Entendendo a Inferência Bayesiana

Métodos Computacionalmente Intensivos

Ana Beatriz Tozzo Martins  Edson Antonio Alves da Silva

Orientador: Prof. PhD. Paulo Justiniano Ribeiro Jr.

25 de agosto de 2007

Este trabalho refere-se ao estudo da apostila do curso de Introdução à Inferência Bayesiana ^{1} reproduzindo (muitas vezes textualmente) aspectos teóricos e resolvendo exercícios e exemplos presentes no material. Foi utilizado também a apostila de Métodos Computacionalmente Intensivos em Estatística ^{2}.

1 Monte Carlo Simples

A idéia do método é escrever a integral que se deseja calcular, por exemplo,

I = ∫_{a}^{b} g(θ) dθ

como o valor esperado

I = (b − a) ∫_{a}^{b} g(θ)\frac{1}{b − a} dθ = (b − a)E[g(θ)]

onde θ é uma variável aleatória com distribuição U(a; b). A estimativa de Monte Carlo

Î = (b − a)\frac{1}{n} ∑_{i=1}^{n} g(θ_{i})

pode ser obtida usando o seguinte algoritmo

1. Gerar θ_{1}, …, θ_{n} da distribuição U(a; b);

2. Calcular g(θ_{1}), …, g(θ_{n});

3. Calcular a média amostral ḡ = ∑_{i=1}^{n} \frac{g(θ_{i})}{n};

4. Calcular Î = (b − a) × ḡ.

Exemplo 1: Seja X uma variável aleatória com distribuição exponencial de parâmetro λ = 1, ou seja, f(x) = e^{−x}, x > 0. Calcular P (1 ≤ x ≤ 3) = P (x ≤ 3) − P (x ≤ 1).

Solução:

1Ehlers, R.S. (2005) Introdução à Inferência Bayesiana. Departamento de Estatística, UFPR. Disponível http://leg.est.ufpr.br/∼ ehlers/notas/bayes2006. Acesso em: agosto de 2007.

2Ehlers, R.S. (2003) Métodos Computacionalmente Intensivos em Estatística, UFPR. Disponível http://leg.est.ufpr.br/∼ ehlers/notas/mci.pdf. Acesso em: agosto de 2007.

• Usando a função pexp do R vamos calcular:

P (1 ≤ x ≤ 3) = ∫_{1}^{3} e^{−x}  (1)

> int.exp = pexp(3, 1) - pexp(1, 1)
> int.exp
[1] 0.3180924

• Usando simulação de Monte Carlo vamos calcular a estimativa de Monte Carlo dada por:

(3 − 1)\frac{1}{10} ∑_{i=1}^{10} e^{−x} x : 1, …, 3  (2)

Os passos necessários são:

– Gerar n = 10 valores (λ_{i}) de uma distribuição uniforme no intervalo (1; 3);

– Calcular g(λ_{i});

– Calcular a média ḡ dos g(λ_{i});

– Calcular (3 − 1) × ḡ.

> n = 10
> x = runif(n, 1, 3)
> g = exp(-x)
> (int.exp = (3 - 1) * mean(g))
[1] 0.3156973

Como se pode observar, a precisão na estimativa foi de 2 casas decimais.

• Escrevendo uma função geral,para o R, cujos argumentos serão o número de simulações e limites de integração (intervalo que se pretende calcular a integral).

> int.exp = function(n, a, b) {
+     x = runif(n, a, b)
+     g = exp(-x)
+     int.exp = (b - a) * mean(g)
+     return(int.exp)
+ }

Para o cálculo de P (1 ≤ x ≤ 3) com n = 20 simulações fazemos:

> int.exp(20, 1, 3)
[1] 0.3477491

Observamos que a precisão foi de uma casa decimal em relação aos valores obtidos anteriormente. Uma vantagem em escrever uma função é que podemos repetir facilmente os cálculos. Por exemplo, para obter 20 resultados, cada um com 10 simulações no intervalo (1; 3) fazemos:

> m = NULL
> for (i in 1:20) {
+     m = c(m, int.exp(10, 1, 3))
+ }
> m

 [1] 0.2698334 0.3327875 0.3037126 0.2612778 0.2880306 0.2747677 0.3650898
 [8] 0.2802789 0.3272915 0.3412805 0.3277631 0.3048185 0.3085896 0.3780798
[15] 0.2055650 0.3920031 0.3175781 0.3608137 0.3136632 0.3908758
> summary(m)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
 0.2056  0.2861  0.3156  0.3172  0.3462  0.3920

Logo, aumentando-se o número de simulações observou-se que a precisão da estimativa 0.3172 foi de 2 casas decimais como nos itens (a) e (b). Podemos, também, calcular a esperança de uma função g(x) da v.a. X cuja f.d.p. é p(x). Para isto basta simular valores de p(x) e calcular a E(X).

> n = 1000
> x = rexp(n, 1)
> mean(x)
[1] 1.033371

Iremos novamente calcular P (1 ≤ x ≤ 3)

> sum(x > 1 & x < 3)/n
[1] 0.317

Novamente observamos que a precisão foi de 2 casas decimais. O que foi feito aqui? Simplesmente calculamos a proporção (frequência relativa) dos valores simulados que “cairam” no intervalo (1; 3), ou seja, a probabilidade procurada.

• Variância do Estimador de Monte Carlo. É importante ter uma idéia do erro cometido nesta aproximação, ou seja, o Erro de Monte Carlo. Esse erro é obtido medindo-se a variância empírica do estimador de Monte Carlo dada por:

v = \frac{1}{n^{2}} ∑_{i=1}^{n}(g(x_{i}) − ḡ)^{2} = \frac{1}{n} (∑_{i=1}^{n} \frac{(g(x_{i}) − ḡ)^{2}}{n} ) .

Sabemos que

V ar(X) = E [X − E(X)]^{2} = ∫ (x − E[X])^{2} · f(x) dx

onde (x − E[X])^{2} é uma função da variável aleatória X. Assim,

> (VX = mean((x - mean(x))^2))
[1] 1.089248
> v = mean(VX/n)
> ep = sqrt(v)
> ep
[1] 0.03300376

2 Monte Carlo via Função de Importância

Segundo Ehler(2006), em muitas situações pode ser muito custoso ou mesmo impossível simular valores da distribuição a posteriori. Neste caso, pode-se recorrer à uma função q(θ) que seja de fácil amostragem, usualmente chamada de função de importância. O procedimento é comumente chamado de amostragem por importância. Se q(θ) for uma função de densidade definida no mesmo espaço variação de θ então:

I = ∫ (\frac{g(θ)p(θ)}{q(θ)}) q(θ)dθ = E [\frac{g(θ)p(θ)}{q(θ)}]

onde a esperança agora é com respeito a distribuição q(θ). Assim, se dispomos de uma amostra aleatória θ_{1}, …, θ_{n} tomada da distribuição q(θ) o estimador de Monte Carlo da integral acima fica:

Î = \frac{1}{n} ∑_{i=1}^{n} \frac{g(θ_{i})p(θ_{i})}{q(θ_{i})}

Exemplo 2: Tomemos uma única observação de uma v.a. X ∼ N(θ; 1) sendo θ desconhecido. A experiência ou conhecimento prévio do parâmetro θ como média da distribuição de X leva a supor que θ ∼ Cauchy(0; 1).

Solução

• Determinação da distribuição posteriori.

– g(x|θ) ∝ e^{−0,5 (x−θ)^{2}}.

– p(θ) = \frac{1}{π (1 + θ^{2})}.

– p(θ|x) = \frac{(\frac{1}{π (1 + θ^{2})}) \frac{1}{\sqrt{2 π}} e^{−0,5 (x−θ)^{2}}}{∫ (\frac{1}{π (1 + θ^{2})}) \frac{1}{\sqrt{2 π}} e^{−0,5 (x−θ)^{2}}d θ}.

– E[θ|x] = ∫ θ p(θ|x)d θ = \frac{∫ θ (\frac{1}{1 + θ^{2}} ) e^{−(x−θ)^{2}}d θ}{∫ (\frac{1}{1 + θ^{2}} ) e^{−(x−θ)^{2}}d θ}.

– V ar(θ|x) = E [θ^{2}|x] − (E [θ|x])^{2}.

• Gráfico da priori e da verossimilhança.

     > options(width = 80)
     > x = rnorm(1, 2, 1)
     > par(mfrow = c(1, 1), mar = c(3.5, 3.5, 0.5, 0.5), mgp = c(2,
     +     0.8, 0))
     > curve(dnorm(x, 2.666545, 1), lty = 1, from = -3, to = 8, ylab = "",
     +     xlab = expression(theta))
     > curve(dcauchy(x, 0, 1), from = -3, to = 8, add = T, lty = 2)
     > legend(4, 0.35, legend = c("priori (Cauchy)", "veross.(Normal)"),
     +     lty = c(2, 1))

• Estimativa pontual de θ (média a posteriori).

Obtenção de E [θ| x]:

a) Gerar θ_{1}, …, θ_{n} (n = 1000), independentes, da distribuição N(x; 1); Figura 1: Distribuição a priori (Cauchy(0;1) e verossimilhança (Normal[2,7 ;1])

> n = 1000
> set.seed(234)
> x = rnorm(1, 2, 1)
> theta = rnorm(n, x, 1)

b) Calcular g_{i} = \frac{θ_{i}}{1 + θ_{i}^{2}} e g_{i}^{∗} = \frac{1}{1 + θ_{i}^{2}} ;

> g.num = theta/(1 + theta^2)
> g.den = 1/(1 + theta^2)

c) Calcular Ê(θ|x) = \frac{∑_{i=1}^{n} g_{i}}{∑_{i=1}^{n} g_{i}^{∗}}.

> media.theta = mean(g.num)/mean(g.den)
> media.theta
[1] 1.877075

• Variância do estimador de θ.

Determina-se aqui V ar(θ|x) = E [θ^{2}|x] − (E [θ|x])^{2}.

a) cálculo de g_{i} = \frac{θ_{i}^{2}}{1 + θ_{i}^{2}} ;

> g.num2 = theta^2/(1 + theta^2)

b) cálculo de E [θ^{2}|x];

> media.theta2 = mean(g.num2)/mean(g.den)

c) cálculo de E [θ^{2}|x] − E^{2} [θ|x].

> var.theta = media.theta2 - (media.theta)^2
> var.theta
[1] 1.065832

• Gráfico da distribuição a posteriori de θ.

a) criar um intervalo de variação para θ;

b) para cada valor θ do intervalo, calcular p(θ|x) na sua forma reduzida (sem a constante normalizadora);

c) construir o gráfico do resultado (θ ∈ (−2; 5)).

> par(mar = c(4, 4, 2, 0.5), mgp = c(3, 0.8, 0))
> curve((1/(pi * (1 + media.theta^2))) * ((1/sqrt(2 * pi)) * exp(-0.5 *
+     (x - media.theta)^2)), from = -2, to = 5, ylab = expression(p(theta/x)),
+     xlab = expression(theta), las = 1)

Figura 2: Distribuição da posteriori na forma analítica sem a constante normalizadora.

• Generalização do problema para k observações.

a) distribuição posteriori;

p(θ|x) ∝ \frac{1}{(1 + θ^{2})} e^{−0,5 k (x̄−θ)^{2}}

b) simulando k = 50 observações de uma população N(θ = 2; σ^{2} = 1);

> x = rnorm(50, 2, 1)

c) simulando 1000 observações de θ ∼ N(x̄; σ^{2}/\sqrt{k} = 1/\sqrt{50});

> n = 1000
> theta = rnorm(n, mean(x), 1/sqrt(50))

d) cálculo de g_{i} = \frac{θ_{i}}{1 + θ_{i}^{2}} e g_{i}^{∗} = \frac{1}{1 + θ_{i}^{2}} ;

> g.num = theta/(1 + theta^2)
> g.den = 1/(1 + theta^2)

e) cálculo de Ê(θ|x) = \frac{∑_{i=1}^{n} g_{i}}{∑_{i=1}^{n} g_{i}^{∗}};

> media.theta = mean(g.num)/mean(g.den)
> media.theta
[1] 2.072455

f) cálculo de E [θ^{2}|x];

> g.num2 = theta^2/(1 + theta^2)
> media.theta2 = mean(g.num2)/mean(g.den)

g) cálculo de E [θ^{2}|x] − E^{2} [θ|x].

> var.theta = media.theta2 - (media.theta)^2
> var.theta
[1] 0.02010174

3 Método de Aceitação-Rejeição

O método consiste consiste em gerar um valor θ^{∗} de uma distribuição auxiliar q(θ) (costuma-se utilizar a priori) e aceitar esse valor como sendo da distribuição a posteriori com probabilidade p(θ|x)/Aq(θ), onde A é uma constante finita tal que p(θ|x) < Aq(θ). Para efeitos práticos, tomase A como sendo o valor máximo da função de verossimilhança, ou seja, A = p(x|θ̂) onde θ̂ é o estimador de máxima verossimilhança de θ. Assim, a probabilidade de aceitação e rejeição se torna p(x|θ^{∗})/p(x|θ̂). O método também funciona se usarmos a versão não normalizada da posteriori dada por: p(θ|x) ∝ p(x|θ) p(θ). Se o processo de maximizar a função de verossimilhança for complexo, sugere-se o método da reamostragem ponderada (a seguir).

O algoritmo para gerar valores a posteriori é:

1. Gerar um valor θ^{∗} da distribuição a priori;

2. Gerar u ∼ U(0; 1);

3. Aceitar θ^{∗} como um valor da posteriori se u < p(x|θ^{∗})/p(x|θ̂), caso contrário, rejeitar θ^{∗} e retornar ao primeiro passo.

Exemplo 3: Tomemos uma única observação de uma v.a. X ∼ N(θ; 1) sendo θ desconhecido. A experiência ou conhecimento prévio do parâmetro θ como média da distribuição de X leva a supor que θ ∼ Cauchy(0; 1). Suponha, agora, que X_{1}, …, X_{n} ∼ N(θ; 1) e assuma uma priori Cauchy(0; 1) para θ.

Sabemos que:

• p(x|θ) ∝ exp \{\frac{n}{2}(x̄ − θ)^{2}\} é a função de verossimilhança;

• θ̂ = x̄ é o EMV de θ;

• a probabilidade de aceitação será:

\frac{p(x|θ^{∗})}{p(x|θ̂)} = \frac{exp \{\frac{n}{2}(x̄ − θ^{∗})^{2}\}}{exp \{\frac{n}{2}(x̄ − x̄)^{2}\}} = exp \{\frac{n}{2}(x̄ − θ^{∗})^{2}\}.

A aplicação do algoritmo fica:

a) simulando uma amostra de 50 observações de X ∼ N(θ = 2; 1) e obtendo o EMV de θ;

     > x = rnorm(50, 2, 1)
     > x.bar = mean(x)

b) simulando 1000 valores de θ^{∗} a partir da distribuição a priori (Cauchy(0; 1));

     > n = 1000
     > theta = rcauchy(n, 0, 1)

c) simulando 1000 valores u ∼ U(0; 1).

     > u = runif(n, 0, 1)

d) calculando o vetor com as probabilidades para aceitação de cada valor de θ∗ simulado

     > prob = exp(-0.5 * 50 * (theta - x.bar)^2)

e) contando os valores de prob que são maiores que o respectivo u simulado;

> soma = sum(u < prob)

f) determinando a taxa de aceitação.

> (taxa = soma/n)
[1] 0.028

g) a figura 3 mostra o gráfico da priori e da verossimilhança.

> options(width = 80)
> par(mar = c(3.5, 3.5, 2, 0.5), mgp = c(2, 0.8, 0))
> curve(dnorm(x, x.bar, 1/sqrt(50)), from = -2, to = 3, ylab = "",
+     xlab = expression(theta))
> curve(dcauchy(x, 0, 1), from = -2, to = 3, add = T, lty = 2)
> legend(-1, 2, legend = c("priori", "veross."), lty = c(2, 1))
> rug(theta)

Figura 3: Distribuição da priori e da verossimilhança

Neste exemplo, a taxa de aceitação foi muito baixa (2.8%), significando que, de 1000 valores de θ^{∗} simulados, apenas 28 foram aceitos.

4 Reamostragem Ponderada

O método consiste em gerar os valores de uma distribuição auxiliar, sem a necessidade da maximização da verossimilhança. A desvantagem do método é que os valores obtidos serão apenas aproximadamente distribuídos segundo a posteriori. O algoritmo consiste basicamente de:

1. Gerar valores θ_{1}, …, θ_{n} da distribuição a priori;

2. Calcular os pesos w_{1}, …, w_{n}. De uma forma geral, esses pesos são calculados como:

w_{i} = \frac{p(θ_{i}|x)/q(θ_{i})}{∑_{j=1}^{n} p(θ_{j}|x)/q(θ_{j})}, i : 1 …, n

Se tomarmos a priori como a densidade auxiliar, ou seja, q(θ) = p(θ), o peso se reduz a:

w_{i} = \frac{p(x|θ_{i})}{∑_{j=1}^{n} p(x|θ_{j})}, i : 1 …, n;

3. Reamostrar valores com probabilidades w_{1}, …, w_{n}.

Exemplo 4: Em um modelo de regressão linear simples temos que y_{i} ∼ N(βx_{i}; 1). Os dados observados são y = (−2, 0, 0, 0, 2) e x = (−2, −1, 0, 1, 2), e usamos uma priori vaga N(0; 4) para β. Faça inferência sobre β obtendo uma amostra da posteriori usando reamostragem ponderada. Compare com a estimativa de máxima verossimilhança β̂ = 0, 8.

Solução:

• A figura 4 mostra o diagrama de dispersão com a reta de regressão linear ajustada por mínimos quadrados:

     > par(mar = c(3.5, 3.5, 0.5, 0.5), mgp = c(2, 0.8, 0))
     > plot(c(-2, -1, 0, 1, 2), c(-2, 0, 0, 0, 2), xlab = "X", ylab = "Y")
     > abline(lm(c(-2, -1, 0, 1, 2) ~ c(-2, 0, 0, 0, 2)))

Figura 4: H-scaterplot de X versus Y e reta de regressão

• Aplicando o algoritmo: a) Entrando com os valores observados de X e Y;

> x = c(-2, -1, 0, 1, 2)
> y = c(-2, 0, 0, 0, 2)

b) Gerando 1000 valores β_{i} ∼ N(0; 2^{2}) da distribuição priori;

> n = 1000
> beta = rnorm(n, 0, 2)

c) Calculando a verossimilhança l(θ_{i}) = p(x|θ_{i});

> l = sapply(beta, function(b) exp(-0.5 * (sum((y - b * x)^2))))

d) Calculando os pesos w_{i};

> w = l/sum(l)

e) Reamostrando 500 β’s com probabilidade w. Serão escolhidos os 500 valores de β correspondentes às 500 maiores probabilidades w;

> m = 500
> beta.resample = sample(beta, size = m, rep = T, prob = w)

f) Visualizando graficamente o resultado da reamostragem (figura 5).

> hist(beta.resample, main = "")

Figura 5: Histograma dos valores reamostrados

Neste exemplo, o estimador β̂ de β será a média dos valores reamostrados dado por β̂ = 0.7665. Podemos ainda visualizar graficamente os valores reamostrados comparados com a distribuição a priori, conforme a figura 6:

> options(width = 80)
> curve(dnorm(x, 0, 2), from = -3, to = 3, ylab = "priori", xlab = expression(beta))
> rug(beta.resample)

Exemplo 5: Para o mesmo modelo do exemplo 4 e os mesmos dados suponha agora que a variância é desconhecida, i.e., y_{i} ∼ N(βx_{i}; σ^{2}). Usamos uma priori hierárquica para (β; σ^{2}), i.e., β|σ^{2} ∼ N(0; σ^{2}) e σ^{−2} ∼ Gama(0, 1; 0, 1).

Figura 6: Valores reamostrados comparados com a distribuição a priori

(a) Obtenha uma amostra da posteriori de (β; σ^{2}) usando reamostragem ponderada. Neste problema, θ = (β; σ^{2}). Assim, a priori hierárquica será dada por p(β; σ^{2}) = p(β|σ^{2}) p(σ^{2}). Será adotado β|σ^{2} ∼ N(0; σ^{2}) e σ^{−2} ∼ Gamma(0, 1; 0, 1).

(b) Estime β e σ^{2} usando uma aproximação para a média a posteriori. Compare com as estimativas de máxima verossimilhança.

(c) Baseado na amostra do item (a), faça um histograma das distribuições marginais de β e σ^{2}.

Solução item (a):

i) Simulando valores dos parâmetros σ^{2} e β|σ^{2} da priori conjunta;

> inv.sigma2 = rgamma(n, 0.1, 0.1)
> sigma2 = 1/inv.sigma2
> beta = rnorm(n, sqrt(sigma2))
> theta = cbind(beta, sigma2)

ii) Calculando a verossimilhança;

> options(width = 80)
> l = sapply(1:1000, function(i) (theta[i, 2]^(-length(y)/2)) *
+     exp(-0.5 * sum(((y - theta[i, 1] * x)^2)/theta[i, 2])))

iii) Calculando os pesos w_{i};

> w = l/sum(l)

iv) Reamostrando 500 θ = (β, σ^{2}) com probabilidade w. Serão escolhidos os 500 valores de θ correspondentes às 500 maiores probabilidades w;

     > ind = sample(1:1000, size = m, rep = T, prob = w)
     > theta.resample = theta[ind, ]

Solução item (b):

As estimativas obtidas das médias marginais serão:

> options(width = 80)
> theta.hat = as.numeric(c(round(mean(theta.resample[, 1]), 4),
+     round(mean(theta.resample[, 2]), 4)))

\begin{matrix}β̂ & 0.8092 \\ σ̂^{2} & 0.7166 \end{matrix}

Solução item (c):

> par(mfrow = c(1, 2))
> hist(theta.resample[, 1], main = "", xlab = "")
> hist(theta.resample[, 2], main = "", xlab = "")

Figura 7: Histograma dos valores reamostrados para β (esquerda) e σ^{2} (direita)

5 Monte Carlo via Cadeias de Markov - MCMC

A idéia do método consiste em utilizar técnicas de simulação iterativa (repetidas) baseadas em Cadeias de Markov para gerar valores correlacionados. Os métodos mais utilizados são: algoritmo de Metropolis-Hastings e Amostrador de Gibbs. Esses métodos realizam um “passeio” aleatório no espaço paramétrico até que o parâmetro convirja para uma distribuição de interesse.

5.1 Algoritmo de Metropolis-Hastings

O algoritmo de Metropolis-Hastings pode ser desenvolvido nos seguintes passos:

1. Escolher o número n de etapas (número de simulações ou número de passos da cadeia).

2. Iniciar o contador de iterações t, fazendo t = 1.

3. Especificar um valor inicial para θ^{(t)}.

4. Gerar um valor θ′ de uma distribuição proposta q(θ′|θ). Esta distribuição proposta pode depender do estado atual da cadeia.

5. Gerar u ∼ U(0; 1) (priori).

6. Determinar a probabilidade de aceitação α(θ; θ′) dada por:

α(θ; θ′) = min \{1; \frac{p(x|θ′) p(θ′) q(θ|θ′)}{p(x|θ) p(θ) q(θ′|θ)}\}

onde p(x|θ) é a distribuição de interesse.

7. Se u < α(θ; θ′) então o valor de θ′ é aceito e θ^{(t)} = θ′. Caso contrário θ′ é rejeitado e θ^{(t)} = θ^{(t−1)}.

8. Se t = n (é o último passo) ir para o passo 11.

9. Incrementar o contador t fazendo t = t + 1.

10. Voltar para o passo 4.

11. FIM.

Exemplo 6: Em uma certa população de animais sabe-se que cada animal pode pertencer a uma dentre 4 linhagens genéticas com probabilidades

p_{1} = \frac{1}{2} + \frac{θ}{4}; p_{2} = \frac{1 − θ}{4} ; p_{3} = \frac{1 − θ}{4} ; p_{4} = \frac{θ}{4}

sendo 0 < θ < 1 um parâmetro desconhecido. Para qualquer θ ∈ (0; 1) é fácil verificar que p_{i} > 0, i = 1, 2, 3, 4 e p_{1} + p_{2} + p_{3} + p_{4} = 1. Observando-se n animais dentre os quais y_{i} pertencem à linhagem i, o vetor aleatório Y = (y_{1}, y_{2}, y_{3}, y_{4}) tem distribuição multinomial com parâmetros n, p_{1}, p_{2}, p_{3}, p_{4} e portanto,

\begin{aligned}p(y| θ) & = & \frac{n!}{y_{1}!y_{2}!y_{3}!y_{4}!} p_{1}^{y_{1}}p_{2}^{y_{2}}p_{3}^{y_{3}}p_{4}^{y_{4}} \\ & ∝ & (2 + θ)^{y_{1}} (1 − θ)^{y_{2}+y_{3}} θ^{y_{4}} \end{aligned}

Foram observados 197 animais com a quantidade deles nas diferentes categorias dado por y = (125, 18, 20, 34). Gerar uma cadeia de Markov com 1000 valores do parâmetro θ que definirá as proporções p_{1}, p_{2}, p_{3}, p_{4}. Descartar os 100 primeiros valores (amostra de aquecimento).

Solução:

• Escolhemos a priori θ ∼ U(0; 1). Assim, p(θ) = 1, ∀ θ ∈ (0; 1) e a posteriori será dada por p(θ|y) ∝ (2 + θ)^{y_{1}} (1 − θ)^{y_{2}+y_{3}} θ^{y_{4}}.

Tomando como distribuição porposta q(θ), também, a distribuição U(0; 1), então, q(θ) = 1 ∀ θ. Daí, a probabilidade de aceitação será dada por:

\begin{aligned}α(θ; θ′) & = & min \{1; \frac{p(y|θ′)}{p(y|θ)}\} = \\ & = & min \{1; (\frac{2 + θ′}{2 + θ} )^{y_{1}} (\frac{1 − θ′}{1 + θ} )^{y_{2}+y_{3}} (\frac{θ′}{θ} )^{y_{4}}\} \end{aligned}

• Criando a função para o cálculo da probabilidade de interesse p(y|θ).

> p = function(teta, y) {
+     p = (2 + teta)^y[1] * (1 - teta)^(y[2] + y[3]) * teta^y[4]
+     return(p)
+ }

• Criando a função de implementação do algorítimo de simulação de Metropolis-Hastings (Cadeia de Markov);

> options(width = 80)
> metr = function(n, y, p, start) {
+     theta = matrix(NA, nrow = n)
+     theta[1] = start
+     sucess = 0
+     for (i in 2:n) {
+         theta.prop = runif(1)
+         A = p(theta.prop, y)/p(theta[i - 1], y)
+         prob = min(1, A)
+         theta.u = runif(1)
+         if (theta.u < prob) {
+             theta[i] = theta.prop
+             sucess = sucess + 1
+         }
+         else theta[i] = theta[i - 1]
+     }
+     taxa = sucess/n
+     return(list(theta = theta, taxa = round(taxa, 2)))
+ }

• Simulando para um valor inicial θ = 0.2.

> n = 1000
> y = c(125, 18, 20, 34)
> start = 0.2
> m = metr(n, y, p, start)

• Obtendo a taxa e aceitação e a média, excluindo as 100 primeiras simulações.

> options(width = 80)
> theta.hat = as.numeric(c(round(m$taxa, 2), round(mean(m$theta[101:1000]),
+     4)))

\begin{matrix}Taxa & 0.17 \\ μ̂ & 0.6251 \end{matrix}

Exemplo 7:

Suponha que queiramos simular valores X ∼ N(0; 1) propondo valores Y ∼ N(x; σ^{2}).

Solução:

• Neste exemplo, a probabilidade de aceitação-rejeição fica:

α(x; y) = min \{1; \frac{π(y)}{π(x)} \frac{q(x|y)}{q(y|x)}\}

α(x; y) = min \{1; \frac{exp\{−0.5 y^{2}\}}{exp\{−0.5 x^{2}\}} \frac{exp \{(x − y)^{2}\}}{exp \{(y − x)^{2}\}}\}

α(x; y) = min \{1; exp \{−\frac{1}{2} (y^{2} − x^{2})\}\}  (3)

Como podemos notar na equação 3, houve uma simplificação no quociente das “propostas”, devido à simetria da função q(x|y) e q(y|x). Esta é a situação de um processo denominado random walk no qual o algoritmo dá saltos aleatórios cobrindo mais eficientemente o espaço paramétrico em torno desse novo centro. Neste caso, taxas de aceitação entre 20 e 30

• Definindo uma função para executar o algoritmo.

> options(width = 80)
> metrop = function(n, sigma) {
+     x = matrix(NA, nrow = n)
+     x[1] = 0
+     sucess = 0
+     for (i in 2:n) {
+         y = rnorm(1, x[i - 1], sigma)
+         prob = min(1, exp(-0.5 * (y^2 - x[i - 1]^2)))
+         u = runif(1)
+         if (u < prob) {
+             x[i] = y
+             sucess = sucess + 1
+         }
+         else x[i] = x[i - 1]
+     }
+     return(list(x = x, taxa = round(sucess/n, 4)))
+ }

Nesta função, n representa o número de simulações, sigma é o desvio padrão da “proposta”, y é simulado da distribuição da “proposta” com média igual ao valor de x obtido no “passo” anterior e os x são os valores da distribuição que desejamos simular.

• Simulando 1000 valores com diferentes valores de sigma da “proposta q()”.

> n = 1000
> sigma = c(0.05, 2, 4)
> m1 = metrop(n, sigma[1])
> m2 = metrop(n, sigma[2])
> m3 = metrop(n, sigma[3])
> taxa.r = as.numeric(c(m1$taxa, m2$taxa, m3$taxa))
> cat("\n Taxa de aceitação", taxa.r, "\n")
 Taxa de aceitação 0.988 0.504 0.306

• Construindo gráficos dos valores simulados.

> options(width = 80)
> par(mfrow = c(1, 3))
> plot(m1$x, type = "l", xlab = "iterações", ylab = expression(theta))
> acf(m1$x, main = "", xlab = "defasagens", ylab = "autocorrelações")
> qqnorm(m1$x[seq(100, n, 5)], pch = "*", main = "", xlab = "Quantil teorico",
+     ylab = "Quantil amostral")
> qqline(m1$x[seq(100, n, 5)], lwd = 2)

Figura 8: Iterações (esquerda), correlações (centro) de 1000 valores simulados usando o algoritmo de Metropolis-Hastings com σ = 0.05 e taxa de aceitação de 98.8% e qqplot (direita) com burning dos 100 primeiros valores e thinning de 5 em 5.

A figura 8 à esquerda mostra a ausência de estacionariedade sugerindo a não-convergência da cadeia. A mesma figura, ao centro, indica uma alta autocorrelação serial.

O thinning é um procedimento que faz uma espécie de amostragem sistemática no conjunto de valores simulados da distribuição a posteriori do(s) parâmetro(s) de interesse. A idéia central consiste em tomar sequencialmente um valor em um intervalo de k simulados. Neste exemplo, os valores simulados da distribuição a posteriori deveriam corresponder valores de uma distribuição N(0; 1). Podemos verificar se a suposição é correta contruindo um gráfico qqnorm. A figura 8 (direita) ilustra o respectivo qqnorm de θ|y com um thinning de 5 simulações e um burning das 100 primeiras observações onde notamos um afastamento da normalidade em ambas as caudas da distribuição.

> par(mfrow = c(1, 3))
> plot(m2$x, type = "l", xlab = "iterações", ylab = expression(theta))
> acf(m2$x, main = "", xlab = "defasagens", ylab = "autocorrelações")
> qqnorm(m2$x[seq(100, n, 5)], pch = "*", main = "", xlab = "Quantil teorico",
+     ylab = "Quantil amostral")
> qqline(m2$x[seq(100, n, 5)], lwd = 1)

A figura 9 à esquerda mostra estacionariedade sugerindo a convergência da cadeia. A mesma figura, ao centro, ainda indica uma autocorrelação serial, mas somente para os valores relativamente próximos (menores que o 10º valor anterior simulado).

> options(width = 80)
> par(mfrow = c(1, 3))
> plot(m3$x, type = "l", xlab = "iterações", ylab = expression(theta))
> acf(m3$x, main = "", xlab = "defasagens", ylab = "autocorrelações")
> qqnorm(m3$x[seq(100, n, 5)], pch = "*", main = "", xlab = "Quantil teorico",
+     ylab = "Quantil amostral")
> qqline(m3$x[seq(100, n, 5)], lwd = 1)

Figura 9: Iterações (esquerda), correlações (centro) de 1000 valores simulados usando o algoritmo de Metropolis-Hastings com σ = 2 e taxa de aceitação de 50.4% e qqplot (direita) com burning dos 100 primeiros valores e thinning de 5 em 5.

Figura 10: Iterações (esquerda), correlações (centro) de 1000 valores simulados usando o algoritmo de Metropolis-Hastings com σ = 4 e taxa de aceitação de 30.6% e qqplot (direita) com burning dos 100 primeiros valores e thinning de 5 em 5.

A figura 10 à esquerda mostra estacionariedade indicando a convergência da cadeia. A mesma figura, à direita, ainda indica uma alta autocorrelação serial.

> par(mfrow = c(1, 3))
> hist(m1$x[seq(100, n, 5)], main = "", xlab = "")
> hist(m2$x[seq(100, n, 5)], main = "", xlab = "")
> hist(m3$x[seq(100, n, 5)], main = "", xlab = "")

• Determinando as estimativas da posteriori.

As estimativas de θ em relação a cada variância atribuída à “proposta” com “burning” das 100 primeiras observações resulta em:

> theta.hat = as.numeric(c(round(mean(m1$x[100:n]), 4), round(mean(m2$x[100:n]),
+     4), round(mean(m3$x[100:n]), 4)))

\begin{matrix}σ_{Y} & 0.05 & 2 & 4 \\ μ̂_{X} & -0.4872 & 0.0989 & 0.0057 \end{matrix}

Figura 11: Distribuição da posteriori para σ_{1} = 0.05 (esquerda), σ_{2} = 2 (centro) e σ_{3} = 4 (direita).

Como esperado, as estimativas obtidas estão próximas de zero já que X ∼ N(0; 1).

Exemplo 8: Suponha que queiramos simular valores X ∼ N(0; 1) propondo valores Y ∼ Cauchy(0; σ).

(a) Para a distribuição proposta (Cauchy) selecione experimentalmente o valor de σ que maximiza a taxa de aceitação.

(b) Para este valor de σ faça os gráficos dos valores simulados da cadeia ao longo das iterações e verifique se há indicação de convergência.

(c) Repita os itens anteriores com a distribuição proposta Cauchy(θ; σ).

A distribuição Cauchy deriva de uma distribuição t-student com 1 grau de liberdade. Considerando a distribuição t-student:

f(t|μ; σ^{2}) = \frac{Γ (\frac{n + 1}{2}) n^{n/2}}{Γ (\frac{n}{2})} [n + \frac{(t − μ)^{2}}{σ^{2}} ]^{− \frac{(n + 1)}{2}} t ∈ ℝ

então, para n = 1 grau de liberdade, temos:

Cauchy(y|μ; σ) = (π σ [1 + (\frac{y − μ}{σ})^{2}])^{−1}

que representa a distribuição de Cauchy com parâmetro de locação μ e parâmetro de escala σ.

Solução item (a):

• Obtendo a expressão da probabilidade de aceitação/rejeição.

\begin{aligned}α(x; y) & = & min \{1; \frac{π(y)}{π(x)} \frac{q(x|y)}{q(y|x)}\} \\ α(x; y) & = & min \{1; \frac{exp\{−0.5 y^{2}\}}{exp\{−0.5 x^{2}\}} \begin{matrix}(π σ [1 + (\frac{x − 0}{σ})^{2}]) & −1 \\ (π σ [1 + (\frac{y − 0}{σ})^{2}]) & −1 \end{matrix} \} \\ α(x; y) & = & min \{1; exp \{−\frac{1}{2} (y^{2} − x^{2})\} (\frac{σ^{2} + y^{2}}{σ^{2} + x^{2}} )\} \end{aligned}

• Definindo a função para executar o algoritmo.

> metrop = function(n, sigma) {
+     x = matrix(NA, nrow = n)
+     x[1] = 0
+     sucess = 0
+     for (i in 2:n) {
+         y = rcauchy(1, 0, sigma)
+         prob = min(1, exp(-0.5 * (y^2 - x[i - 1]^2)) * (y^2 +
+             sigma^2)/(x[i - 1]^2 + sigma^2))
+         u = runif(1)
+         if (u < prob) {
+             x[i] = y
+             sucess = sucess + 1
+         }
+         else x[i] = x[i - 1]
+     }
+     return(list(x = x, taxa = round(sucess/n, 4)))
+ }

Nesta função, n representa o número de simulações, sigma é o parâmetro de escala da “proposta”, y é simulado da distribuição da “proposta” com média igual zero e σ como parâmetro de escala.

• Simulando 1000 valores com diferentes valores de sigma da “proposta q()”.

> n = 1000
> sigma = c(0.545, 2, 3)
> m1 = metrop(n, sigma[1])
> m2 = metrop(n, sigma[2])
> m3 = metrop(n, sigma[3])
> taxa.r = as.numeric(c(m1$taxa, m2$taxa, m3$taxa))
> cat("\n Taxa de aceitação", taxa.r, "\n")
 Taxa de aceitação 0.693 0.439 0.32

Atribuindo valores ao vetor σ observamos, após várias e sucessivas tentativas, que o valor máximo ou maior taxa de aceitação, 69.3%, ocorreu para σ = 0.545.

Solução item (b):

• Construindo gráficos dos valores simulados para σ_{1} = 0.545.

> par(mfrow = c(1, 3))
> plot(m1$x, type = "l", xlab = "iterações", ylab = expression(theta))

> acf(m1$x, main = "", xlab = "defasagens", ylab = "autocorrelações")
> qqnorm(m1$x[seq(100, n, 5)], pch = "*", main = "", xlab = "Quantil teorico",
+     ylab = "Quantil amostral")
> qqline(m1$x[seq(100, n, 5)], lwd = 1)

Figura 12: Iterações (esquerda), correlações (centro) de 1000 valores simulados usando o algoritmo de Metropolis-Hastings com σ = 0.545 e taxa de aceitação de 69.3% e qqplot (direita) com burning dos 100 primeiros valores e thinning de 5 em 5.

• Construindo gráficos dos valores simulados para σ_{2} = 2

> par(mfrow = c(1, 3))
> plot(m2$x, type = "l", xlab = "iterações", ylab = expression(theta))
> acf(m2$x, main = "", xlab = "defasagens", ylab = "autocorrelações")
> qqnorm(m2$x[seq(100, n, 5)], pch = "*", main = "", xlab = "Quantil teorico",
+     ylab = "Quantil amostral")
> qqline(m2$x[seq(100, n, 5)], lwd = 1)

Figura 13: Iterações (esquerda), correlações (centro) de 1000 valores simulados usando o algoritmo de Metropolis-Hastings com σ = 2 e taxa de aceitação de 43.9% e qqplot (direita) com burning dos 100 primeiros valores e thinning de 5 em 5.

• Construindo gráficos dos valores simulados para σ_{3} = 3

> options(width = 80)
> par(mfrow = c(1, 3))
> plot(m3$x, type = "l", xlab = "iterações", ylab = expression(theta))

> acf(m3$x, main = "", xlab = "defasagens", ylab = "autocorrelações")
> qqnorm(m3$x[seq(100, n, 5)], pch = "*", main = "", xlab = "Quantil teorico",
+     ylab = "Quantil amostral")
> qqline(m3$x[seq(100, n, 5)], lwd = 1)

Figura 14: Iterações (esquerda), correlações (centro) de 1000 valores simulados usando o algoritmo de Metropolis-Hastings com σ = 3 e taxa de aceitação de 32% e qqplot (direita) com burning dos 100 primeiros valores e thinning de 5 em 5.

As figuras 12, 13 e 14, à esquerda, mostram estacionariedade sugerindo a convergência das cadeias. Ao centro, os gráficos revelam alta autocorrelação e a normalidade dos dados pode ser verificada pelos gráficos à direita.

• Distribuição da posteriori para σ_{1} = 0.545, 2 e 3.

> par(mfrow = c(1, 3))
> hist(m1$x[seq(100, n, 5)], main = "", xlab = "")
> hist(m2$x[seq(100, n, 5)], main = "", xlab = "")
> hist(m3$x[seq(100, n, 5)], main = "", xlab = "")

Figura 15: Distribuição da posteriori para σ = 0.545 (esquerda), σ_{2} = 2 (centro) e σ_{3} = 3 (direita).

• Determinando as estimativas da posteriori. As estimativas de θ em relação a cada variância atribuída a “proposta” com “burning” das 100 primeiras observações resulta em:

> theta.hat = as.numeric(c(round(mean(m1$x[100:n]), 4), round(mean(m2$x[100:n]),
+     4), round(mean(m3$x[100:n]), 4)))

\begin{matrix}σ_{Y} & 0.545 & 2 & 3 \\ μ̂_{X} & 0.048 & -0.0763 & -0.0396 \end{matrix}

Novamente, sendo X ∼ N(0; 1), as estimativas se aproximaram de zero.

• Densidade de uma distribuição Cauchy com parâmetros de dispersão iguais a σ =0.545, σ = 2 e σ = 3.

> par(mfrow = c(1, 1))
> x <- seq(-4, 4, l = 501)
> C1 <- dcauchy(x, 0, sigma[1])
> C2 <- dcauchy(x, 0, sigma[2])
> C3 <- dcauchy(x, 0, sigma[3])
> plot(x, C1, ty = "l")
> lines(x, C2, lty = 2)
> lines(x, C3, lty = 3)

Figura 16: Densidade de uma distribuição Cauchy com parâmetros de dispersão iguais a σ =0.545 (linha contínua), σ = 2 (linha pontilhada maior) e σ = 3 (linha pontilhada menor)

5.2 Algoritmo Amostrador de Gibbs

No método amostrador de Gibbs, as transações de estado são feitas de acordo com as distribuições condicionais completas, dadas por:

π(θ_{i}|θ_{−i}) = \frac{π(θ)}{∫ π(θ) d θ_{i}}  (4)

onde θ_{−i} = 𝕀_{−i}θ, sendo θ = (θ_{1}, θ_{2}, …, θ_{p})′, i = 1, 2, …, p, ou seja:

θ_{−i} = [ θ_{1} ; θ_{2} ; ⋮ ; θ_{i−1} ; θ_{i+1} ; ⋮ ; θ_{p}]_{(p−1)×1} = [\begin{matrix}1 & 0 & … & 0 & 0 & 0 & … & 0 \\ 0 & 1 & … & 0 & 0 & 0 & … & 0 \\ ⋮ & ⋮ & ... & ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & … & 1 & 0 & 0 & … & 0 \\ 0 & 0 & … & 0 & 0 & 1 & … & 0 \\ ⋮ & ⋮ & ... & ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & … & 0 & 0 & 0 & … & 1 \end{matrix}]_{(p−1)×p} [θ_{1} ; θ_{2} ; ⋮ ; ⋮ ; θ_{p}]_{p×1}

A matriz 𝕀_{−i} corresponde a uma matriz identidade excluída a i-ésima linha. Em outras palavras, podemos dizer que a distribuição condicional completa é a distribuição da i-ésima componente de θ condicionada em todas as outras componentes. A equação 4 também pode ser escrita como:

π(θ) = π(θ_{i}|θ_{−i})π_{θ_{−i}}(θ_{−i}).

Queremos gerar uma amostra de π(θ) onde a transição se dá da seguinte forma:

• especificamos valores para um vetor inicial θ^{(0)};

• calcula-se a probabilidade de θ_{1} condicionada aos θ′s restantes com os valores especificados inicialmente;

• nas transições de θ_{r−1} para θ_{r}, r = 2, 3, …, p, calcula-se, para θ_{r}, a probabilidade π(θ_{r}|θ_{−r}) substituindo em θ_{−r} o valor de θ_{r−1} calculado no passo anterior.

Segundo Paulino, Turkman e Murteira (2003)^{3} este procedimento é o adotado pelo algoritmo de Metropolis-Hastings com q(θ; θ′) e q(θ′; θ) positivas se, e somente se, θ e θ′ diferirem no máximo em uma coordenada.

Se as distribuições condicionais completas forem conhecidas, então o algoritmo será desenvolvido como:

1. Inicie o contador de iterações em t = 0 e defina o número n de iterações;

2. Especifique valores iniciais para θ^{(0)};

3. Avance t fazendo t = t + 1 e obtenha θ^{(t)} a partir de θ^{(t−1)} por geração sucessiva de valores como:

θ_{1}^{(t)} ∼ π(θ_{1}| θ_{2}^{(t−1)}, θ_{3}^{(t−1)}, θ_{4}^{(t−1)}, …, θ_{p}^{(t−1)})
θ_{2}^{(t)} ∼ π(θ_{2}| θ_{1}^{(t)},  θ_{3}^{(t−1)}, θ_{4}^{(t−1)}, …, θ_{p}^{(t−1)})
θ_{3}^{(t)} ∼ π(θ_{3}| θ_{1}^{(t)},  θ_{2}^{(t)},  θ_{4}^{(t−1)}, …, θ_{p}^{(t−1)})
⋮
θ_{p}^{(t)} ∼ π(θ_{p}| θ_{1}^{(t)},  θ_{2}^{(t)},  θ_{3}^{(t)},  …, θ_{p−1}^{(t)})

4. Se t = n, encerrar;

5. Voltar a etapa 3.

Exemplo 9: Em um processo de contagem no qual foram observados Y_{1}, Y_{2}, …, Y_{m}, Y_{m+1}, …, Y_{n}, suspeita-se que houve um ponto de mudança m tal que:

\begin{aligned}Y_{i} & ∼ & P oisson(λ), & i = 1, …, m \\ Y_{i} & ∼ & P oisson(ϕ), & i = (m + 1), …, n. \end{aligned}

3Paulino, Carlos Daniel; Turkman, M. Antónia Amaral; Murteira, Bento. Estatística Bayesiana. Lisboa: Fundação Calouste Gulbenkian, 2003. ISBN 972-31-1043-1

Assumindo-se as distribuições a priori independentes λ ∼ Gamma(a; b), ϕ ∼ Gamma(c; d) e p(m) = 1/n, estimar o ponto de mudança m e os parâmetros λ e ϕ (CARLIN, et all)^{4}.

Solução

• Calculando a posteriori.

\begin{aligned}p(λ; ϕ; m| y) & ∝ & (∏_{i=1}^{m} e^{−λ}λ^{y_{i}}) (∏_{i=1}^{m} e^{−ϕ}ϕ^{y_{i}}) λ^{a−1}e^{−bλ}ϕ^{c−1}e^{−d ϕ} (\frac{1}{n}) \\ p(λ; ϕ; m| y) & ∝ & λ^{a+t_{1}−1}e^{−(b+m)λ}ϕ^{c+t_{2}−1}e^{−(d+n−m)ϕ} (\frac{1}{n}) \end{aligned}

onde t_{1} = ∑_{i−1}^{m} y_{i} e t_{2} = ∑_{i=m+1}^{n} y_{i}.

• Calculando as distribuições condicionais completas.

(a) Para o parâmetro λ:

p(λ|ϕ; m; y) = \frac{p(λ; ϕ; m; y)}{p(ϕ; m; y)} = \frac{p(y|λ; ϕ; m)p(λ; ϕ; m)}{p(ϕ; m; y)}

p(λ|ϕ; m; y) = \frac{p(y|λ; ϕ; m)p(λ) p(ϕ) p(m)}{p(ϕ; m; y)}

Como p(m), p(ϕ) e p(ϕ; m; y) não envolvem o parâmetro λ de interesse, então:

p(λ|ϕ; m; y) ∝ p(y|λ; ϕ; m)p(λ)

p(λ|ϕ; m; y) ∝ (∏_{i=1}^{m} e^{−λ}λ^{y_{i}}) (∏_{i=1}^{m} e^{−ϕ}ϕ^{y_{i}}) λ^{a−1}e^{−bλ}

Novamente ∏_{i=1}^{m} e^{−ϕ}ϕ^{y_{i}} não envolve λ, então:

p(λ|ϕ; m; y) ∝ (∏_{i=1}^{m} e^{−λ}λ^{y_{i}}) λ^{a−1}e^{−bλ}

p(λ|ϕ; m; y) ∝ e^{−mλ}λ^{t_{1}}λ^{a−1}e^{−bλ}

p(λ|ϕ; m; y) ∝ λ^{a+t_{1}−1}e^{−(m+b)λ} ∼ Gamma(a + t_{1}; b + m)

(b) Para o parâmetro ϕ: Com um desenvolvimento análogo ao item (a) temos:

p(ϕ|λ; m; y) ∝ ϕ^{c+t_{2}−1}e^{−(d+n−m)ϕ} ∼ Gamma(c + t_{2}; d + n − m)

(c) Para o parâmetro m: Com um desenvolvimento análogo ao item (a) temos:

p(m|λ; ϕ; y) ∝ λ^{t_{1}}e^{−mλ}ϕ^{t_{2}}e^{−(n−m)ϕ}, m = 1, …, n

• Escrevendo a função de implementação do algoritmo.

> rm(list = ls(all = TRUE))
> Gibbs = function(a, b, c, d, y, n, nburn) {
+     N = length(y)
+     lambda = matrix(0, nrow = n)
+     phi = matrix(0, nrow = n)
+     m = matrix(0, nrow = n)
+     lambda[1] = 1
+     phi[1] = 1
+     m[1] = 10
+     for (i in 2:n) {

4.... localizar ...

+         t1 = sum(y[1:m[i - 1]])
+         t2 = 0
+         if (m[i - 1] < N) {
+             t2 = sum(y[(m[i - 1] + 1):N])
+         }
+         lambda[i] = rgamma(1, (a + t1), (b + m[i - 1]))
+         phi[i] = rgamma(1, (c + t2), (d + N - m[i - 1]))
+         prob = NULL
+         for (j in 1:N) {
+             t1 = sum(y[1:j])
+             t2 = 0
+             if (j < N) {
+                 t2 = sum(y[(j + 1):N])
+             }
+             aux = (lambda[i]^t1) * exp(-j * lambda[i]) * (phi[i]^t2) *
+                 exp(-(N - j) * phi[i])
+             prob = c(prob, aux)
+         }
+         soma = sum(prob)
+         probm = prob/soma
+         m[i] = sample(x = (1:N), size = 1, prob = probm)
+     }
+     return(theta = list(lambda = lambda, phi = phi, m = m))
+ }

• Especificando dados de entrada da função Gibbs.

> a = 0.1
> b = 0.1
> c = 0.1
> d = 0.1
> y = c(rpois(23, 2), rpois(17, 5))
> n = 1000
> nburn = 100

• Passando valores para a função.

> g = Gibbs(a, b, c, d, y, n, nburn)

• Determinando a distribuição de probabilidade dos 1000 valores simulados de m (ponte de mudança), descontados os 100 primeiros valores de Burn in simulados.

> foo1 = prop.table(table(g$m[nburn + 1:n]))
> foo2 = data.frame(m = as.integer(names(foo1)), Freq = as.numeric(foo1))
> library(xtable)
> xtable(foo2, caption = "Tabela de frequências relativas dos valores de corte simulados",
+     "tab:one")

A tabela 1 indica que o valor de m=22 aparece com uma frequência de 34.22% dos resultados, indicando ser o valor do ponto de mudança com maior probabilidade.

• Gerando gráficos das cadeias simuladas, após um período de “aquecimento” (Burn in) para o parâmetro m.

> par(mar = c(3.5, 3.5, 0.5, 0.5), mgp = c(2, 0.8, 0), mfrow = c(1,
+     3))
> range = ((nburn + 1):n)
> plot(g$m[range], type = "l", ylab = "m", lwd = 0.5, xlab = "Iteracao")

m Freq
1 12 0.00
2 18 0.01
3 19 0.00
4 20 0.14
5 21 0.20
6 22 0.34
7 23 0.23
8 24 0.08
9 25 0.00
10 26 0.00
11 28 0.00

Tabela 1: Tabela de frequências relativas dos valores de corte simulados

> acf(g$m[range], main = "", xlab = "defasagem", ylab = "autocorrelações")
> plot(foo1, ylab = "Frequencia relativa", xlab = "Ponto de mudança")

Figura 17: Iterações, autocorrelações e distribuição de frequências relativa do ponto de mudança do processo de Poisson.

• Gerando gráficos das cadeias simuladas, após um período de “aquecimento” (Burn in) para o parâmetro λ.

> par(mar = c(3.5, 3.5, 0.5, 0.5), mgp = c(2, 0.8, 0), mfrow = c(1,
+     3))
> range = ((nburn + 1):n)
> plot(g$lambda[range], type = "l", ylab = expression(lambda),
+     lwd = 0.5, xlab = "Iteracao")
> acf(g$lambda[range], main = "", xlab = "defasagem", ylab = "autocorrelações")
> plot(density(g$lambda[range]), cex = 0.7, lwd = 0.5, main = "",
+     xlab = expression(lambda))

• Gerando gráficos das cadeias simuladas, após um período de “aquecimento” (Burn in) para o parâmetro ϕ.

> par(mar = c(3.5, 3.5, 0.5, 0.5), mgp = c(2, 0.8, 0), mfrow = c(1,
+     3))
> range = ((nburn + 1):n)
> plot(g$phi[range], type = "l", lwd = 0.5, ylab = expression(phi),
+     xlab = "Iteracao")
> acf(g$phi[range], main = "", xlab = "defasagem", ylab = "autocorrelações")
> plot(density(g$phi[range]), cex = 0.7, lwd = 0.5, main = "",
+     xlab = expression(phi))

Figura 18: Iterações, autocorrelações e densidade (distribuição de probabilidades) do parâmetro λ.

Figura 19: Iterações, autocorrelações e densidade (distribuição de probabilidades) do parâmetro ϕ.

• Comparando ϕ com λ pela diferença entre eles.

> par(mar = c(3.5, 3.5, 0.5, 0.5), mgp = c(2, 0.8, 0), mfrow = c(1,
+     3))
> dif = g$phi[range] - g$lambda[range]
> par(mfrow = c(1, 3))
> plot(dif, type = "l", ylab = expression(phi - lambda), xlab = "Iteracao")
> acf(dif, main = "")
> plot(density(dif), main = "", lwd = 0.5, xlab = expression(phi -
+     lambda))
> rug(dif)

Figura 20: Iterações, autocorrelações e densidade da diferença dos parâmetros ϕ e λ.

• Estimação das contagens médias dos processos de Poisson condicionado nos valores de m =22 com maior probabilidade a posteriori.

> xp = matrix(0, nrow = n - nburn, ncol = 3)
> xp[, 1] = g$lambda[range]

> xp[, 2] = g$phi[range]
> xp[, 3] = g$m[range]
> lambda.22 = xp[, 1][xp[, 3] == max.m]
> phi.22 = xp[, 2][xp[, 3] == max.m]
> par(mfrow = c(1, 2))
> plot(lambda.22, type = "l", ylab = "", xlab = "Iteração")
> plot(phi.22, type = "l", ylab = "", xlab = "Iteração")

Figura 21: À esquerda, densidade de (λ|m = 22) e a direita, densidade de (ϕ|m = 22).