Amostrador de Gibbs

1 Exemplo 1

Seja \[\begin{eqnarray*} Y_i &\sim& N\left( \mu, V\right), \quad i=1, \ldots, n. \end{eqnarray*}\] A priori, considere que \[\begin{eqnarray*} \mu &\sim& N\left( 0, V_{\mu}\right)\\ V &\sim& GI(a, b). \end{eqnarray*}\]

A distribuição a posteriori de \(\boldsymbol{\theta}=(\mu, V)\) é \[p(\boldsymbol{\theta}|\boldsymbol{y}) \propto V^{-\left(\frac{n}{2}+a+1\right)} \exp\left\{-\frac{1}{2V}\sum_{i=1}^n\left(Y_i-\mu\right)^2 -\frac{b}{V} -\frac{\mu^2}{2V_{\mu}}\right\}.\] Então as distribuições condicionais completas a posteriori são

\[\begin{eqnarray*} \mu|V, \boldsymbol{y} &\sim& N\left( \left(\frac{n}{V} + \frac{1}{V_{\mu}}\right)^{-1} \left(\frac{\sum_{i=1}^n{Y_i}}{V}\right), \left(\frac{n}{V} + \frac{1}{V_{\mu}}\right)^{-1}\right)\\ V|\mu, \boldsymbol{y} &\sim& GI\left( \frac{n}{2} + a, b+\frac{1}{2}\sum_{i=1}^n\left(Y_i-\mu\right)^2\right) \end{eqnarray*}\]

#########################################
#Yi ~ N(mu, V), i=1,...,n
#mu ~ N(0, Vmu)
#V  ~ GI(a,b)
#########################################

#Gerando uma amostra
n   = 1000
mu  = 5
V   = 1
y   = rnorm(n, mu, sqrt(V))
hist(y, freq=F, main="", cex.lab=1.4, cex.axis=1.4, bty="n", ylab="densidade")
curve(dnorm(x, mu, sqrt(V)), col="red", lwd=2, add=T)

#Supondo mu e V desconhecidos
#Amostrando mu e V pelo Amostrador de Gibbs

#Definindo alguns parâmetros
nite    = 1000
mu_post = NULL
V_post  = NULL

#Definindo os hiperparâmetros (=parâmetros da distribuição a priori)
Vmu = 1000
a   = 2
b   = 1

#Inicializando a cadeia
mu_post[1] = 0
V_post[1]  = 10

#Executando o MCMC
for(ite in 2:nite)
{
  #V_post[ite-1] = V #uso isso caso a amostragem de mu esteja com 
          #problema e eu queira verificar se eh algum erro no script
  #Amostrando mu da DCCP 
  variancia     = 1/(n/V_post[ite-1] + 1/Vmu)
  media         = variancia*(sum(y)/V_post[ite-1])
  mu_post[ite]  = rnorm(1, media, sqrt(variancia))

  #mu_post[ite] = mu #uso isso caso a amostragem de V esteja com 
  #problema e eu queira verificar se eh algum erro no script
  #Amostrando V da DCCP
  V_post[ite]   = 1/rgamma(1, n/2+a , b+0.5*sum((y-mu_post[ite])^2))
}

#Analisando a convergencia
par(mfrow=c(2,2))
##Traço da cadeia com toda a amostra obtida
ts.plot(mu_post,main="traço da cadeia",ylab=expression(mu),xlab="iteração")
abline(h=mu, lwd=2, col="red", lty=3)
ts.plot(V_post,main="traço da cadeia",ylab="V",xlab="iteração")
abline(h=V, lwd=2, col="red", lty=3)
##Traço da cadeia removendo o período de aquecimento
ind = seq(50,nite,by=1)
ts.plot(mu_post[ind],main="traço da cadeia",ylab=expression(mu),xlab="iteração")
abline(h=mu, lwd=2, col="red", lty=3)
ts.plot(V_post[ind],main="traço da cadeia",ylab="V",xlab="iteração")
abline(h=V, lwd=2, col="red", lty=3)

par(mfrow=c(2,2))
#com toda a amostra
acf(mu_post,main="Autocorrelação da cadeia", 
    ylab=expression(mu), xlab="defasagem")
acf(V_post,main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")
#removendo o período de aquecimento
ind = seq(50,nite,by=1)
acf(mu_post[ind],main="Autocorrelação da cadeia", 
    ylab=expression(mu), xlab="defasagem")
acf(V_post[ind],main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")

#Obtendo estimativas pontuais e intervalares com a amostra obtida
#após a remoção do período de aquecimento e o espaçamento
par(mfrow=c(1,2))
hist(mu_post[ind],main="Histograma", freq=F,
    xlab=expression(mu), ylab="densidade", bty="n", lwd=2)
abline(v=mu, lwd=2, col="red", lty=3)
abline(v=quantile(mu_post[ind],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(mu_post[ind],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(mu_post[ind]), lwd=2, lty=4, col="purple")

hist(V_post[ind],main="Histograma",freq=F, 
     xlab="V", ylab="densidade", bty="n", lwd=2)
abline(v=V, lwd=2, col="red", lty=3)
abline(v=quantile(V_post[ind],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(V_post[ind],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(V_post[ind]), lwd=2, lty=4, col="purple")

Rodando mais de uma cadeia para verificar se as cadeias convergirão para a mesma região mesmo inicializando de valores diferentes.

#########################################
#Rodando mais de uma cadeia
#########################################

#Definindo alguns parâmetros
mu_post_2 = mu_post_3 = NULL
V_post_2  = V_post_3 = NULL

#Inicializando a cadeia
mu_post_2[1] = -100
V_post_2[1]  = 100

mu_post_3[1] = 100
V_post_3[1]  = 50

#Executando o MCMC
for(ite in 2:nite)
{
  #2a cadeia
  #Amostrando mu
  variancia     = 1/(n/V_post_2[ite-1] + 1/Vmu)
  media         = variancia*(sum(y)/V_post_2[ite-1])
  mu_post_2[ite]  = rnorm(1, media, sqrt(variancia))
  
  #Amostrando V
  V_post_2[ite]   = 1/rgamma(1, n/2+a , b+0.5*sum((y-mu_post_2[ite-1])^2))
  
  #3a cadeia
  #Amostrando mu
  variancia     = 1/(n/V_post_3[ite-1] + 1/Vmu)
  media         = variancia*(sum(y)/V_post_3[ite-1])
  mu_post_3[ite]  = rnorm(1, media, sqrt(variancia))
  
  #Amostrando V
  V_post_3[ite]   = 1/rgamma(1, n/2+a , b+0.5*sum((y-mu_post_3[ite-1])^2))
}

#Analisando a convergencia com todas as amostras geradas das 3 cadeias
par(mfrow=c(2,3))
ts.plot(mu_post,main="traço da cadeia",ylab=expression(mu),xlab="iteração") #1a cadeia
abline(h=mu, lwd=2, col="red", lty=3) #valor verdadeiro
ts.plot(mu_post_2,main="traço da cadeia",ylab=expression(mu),xlab="iteração")#2a cadeia
abline(h=mu, lwd=2, col="red", lty=3)
ts.plot(mu_post_3,main="traço da cadeia",ylab=expression(mu),xlab="iteração")#3a cadeia
abline(h=mu, lwd=2, col="red", lty=3)

ts.plot(V_post,main="traço da cadeia",ylab="V",xlab="iteração") #1a cadeia
abline(h=V, lwd=2, col="red", lty=3) #valor verdadeiro
ts.plot(V_post_2,main="traço da cadeia",ylab="V",xlab="iteração")#2a cadeia
abline(h=V, lwd=2, col="red", lty=3)
ts.plot(V_post_3,main="traço da cadeia",ylab="V",xlab="iteração")#2a cadeia
abline(h=V, lwd=2, col="red", lty=3)

#Analisando a convergencia com todas as amostras geradas das 3 cadeias
par(mfrow=c(1,2))
ts.plot(mu_post[ind],main="traço da cadeia",ylab=expression(mu),xlab="iteração")
lines(ind,mu_post_2[ind], col="purple")
lines(ind,mu_post_3[ind], col="green")
abline(h=mu, lwd=2, col="red", lty=3)
ts.plot(V_post[ind],main="traço da cadeia",ylab="V",xlab="iteração")
lines(ind,V_post_2[ind], col="purple")
lines(ind,V_post_3[ind], col="green")
abline(h=V, lwd=2, col="red", lty=3)

#Analisando a autocorrelação com todas as amostras geradas das 3 cadeias
par(mfrow=c(2,3))
acf(mu_post,main="Autocorrelação da cadeia", 
    ylab=expression(mu), xlab="defasagem")
acf(mu_post_2,main="Autocorrelação da cadeia", 
    ylab=expression(mu), xlab="defasagem")
acf(mu_post_3,main="Autocorrelação da cadeia", 
    ylab=expression(mu), xlab="defasagem")
acf(V_post,main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")
acf(V_post_2,main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")
acf(V_post_3,main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")

#Analisando a distribuição de interesse e estimativas obtidas
#com base na amostra da 1a cadeia
par(mfrow=c(1,2))
hist(mu_post[ind],main="Histograma", freq=F,
     xlab=expression(mu), ylab="densidade", bty="n", lwd=2)
abline(v=mu, lwd=2, col="red", lty=3)
abline(v=quantile(mu_post[ind],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(mu_post[ind],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(mu_post[ind]), lwd=2, lty=4, col="purple")

hist(V_post[ind],main="Histograma", freq=F,
     xlab="V", ylab="densidade", bty="n", lwd=2)
abline(v=V, lwd=2, col="red", lty=3)
abline(v=quantile(V_post[ind],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(V_post[ind],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(V_post[ind]), lwd=2, lty=4, col="purple")

Exercício 1:

Pesquisar como avaliamos a convergência usando o pacote CODA no R.

#Avaliando a convergência através do pacote CODA
if(!require(coda)){install.packages("coda")}
library(coda)
mu_post           = mcmc(mu_post)
mu_post_2         = mcmc(mu_post_2)
mu_post_3         = mcmc(mu_post_3)
mu_lista          = mcmc.list(mu_post, mu_post_2, mu_post_3)
mu.mcmc           = window(mu_lista, start=50, end=nite)
plot(x = mu.mcmc, trace = TRUE, density = T)

summary(mu.mcmc)
## 
## Iterations = 50:1000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 951 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##      4.9902148      0.0327046      0.0006123      0.0006188 
## 
## 2. Quantiles for each variable:
## 
##  2.5%   25%   50%   75% 97.5% 
## 4.926 4.968 4.990 5.012 5.054
#Este comando fornecerá estatísticas de convergência para cada parâmetro na #análise. Preste atenção especial à estatística Potential scale reduction 
#factor (fator de redução de escala potencial, geralmente denotado como 
#R-hat). Valores de R-hat próximos de 1 indicam convergência adequada, 
#enquanto valores acima de 1 indicam falta de convergência.
gelman.diag(mu.mcmc)
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]          1          1
#A função geweke.diag() retornará os valores Z para cada parâmetro e cada 
#período de tempo. Valores Z dentro do intervalo [-2, 2] são considerados 
#indicativos de convergência adequada.
geweke.diag(mu.mcmc)
## [[1]]
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##    var1 
## -0.3181 
## 
## 
## [[2]]
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##   var1 
## 0.2159 
## 
## 
## [[3]]
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##   var1 
## -0.607
#Além desses diagnósticos, o pacote CODA também oferece outras funções 
#úteis, como autocorr.diag() para verificar a autocorrelação das cadeias e 
#raftery.diag() para estimar o tamanho mínimo da amostra para obter uma 
#estimativa precisa.
autocorr.diag(mu.mcmc)
##                [,1]
## Lag 0   1.000000000
## Lag 1  -0.016962885
## Lag 5  -0.025360828
## Lag 10  0.008901097
## Lag 50 -0.012574016
raftery.diag(mu.mcmc)
## [[1]]
## 
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95 
## 
## You need a sample size of at least 3746 with these values of q, r and s
## 
## [[2]]
## 
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95 
## 
## You need a sample size of at least 3746 with these values of q, r and s
## 
## [[3]]
## 
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95 
## 
## You need a sample size of at least 3746 with these values of q, r and s
V_post   = mcmc(V_post)
V_post_2 = mcmc(V_post_2)
V_post_3 = mcmc(V_post_3)
V_lista  = mcmc.list(V_post, V_post_2, V_post_3)
V.mcmc   = window(V_lista, start=50, end=nite)
plot(x = V.mcmc, trace = TRUE, density = T)

gelman.diag(V.mcmc)
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]          1          1
geweke.diag(V.mcmc)
## [[1]]
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##      var1 
## 0.0005394 
## 
## 
## [[2]]
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##  var1 
## 2.282 
## 
## 
## [[3]]
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##    var1 
## -0.4787
autocorr.diag(V.mcmc)
##               [,1]
## Lag 0  1.000000000
## Lag 1  0.003712328
## Lag 5  0.005782506
## Lag 10 0.027242121
## Lag 50 0.013175670
raftery.diag(V.mcmc)
## [[1]]
## 
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95 
## 
## You need a sample size of at least 3746 with these values of q, r and s
## 
## [[2]]
## 
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95 
## 
## You need a sample size of at least 3746 with these values of q, r and s
## 
## [[3]]
## 
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95 
## 
## You need a sample size of at least 3746 with these values of q, r and s

2 Exemplo 2

Seja \[\begin{eqnarray*} Y_i &\sim& N\left( \beta_0+\beta_1 X_i, V\right), i=1, \ldots, n. \end{eqnarray*}\] A priori, considere que \[\begin{eqnarray*} \beta_0 &\sim& N\left( 0, V_{\beta_0}\right)\\ \beta_1 &\sim& N\left( 0, V_{\beta_1}\right)\\ V &\sim& GI(a, b). \end{eqnarray*}\]

A distribuição a posteriori de \(\boldsymbol{\theta}=(\beta_0, \beta_1, V)\) é \[p(\boldsymbol{\theta}|\boldsymbol{y}) \propto V^{-\left(\frac{n}{2}+a+1\right)} \exp\left\{-\frac{1}{2V}\sum_{i=1}^n\left(Y_i-\beta_0-\beta_1 X_i\right)^2 -\frac{b}{V} -\frac{\beta_0^2}{2V_{\beta_0}}-\frac{\beta_1^2}{2V_{\beta_1}}\right\}.\]

Exercício 2: Calcular as distribuições condicionais completas a posteriori de \(\beta_0\), \(\beta_1\) e \(V\). E implementar o algoritmo Amostrador de Gibbs para obter uma amostra da distribuição a posteriori dada acima.

Então as distribuições condicionais completas a posteriori são

\[\begin{eqnarray*} \beta_0|\beta_1, V, \boldsymbol{y} &\sim& N\left( \left(\frac{n}{V} + \frac{1}{V_{\beta_0}}\right)^{-1} \left(\frac{\sum_{i=1}^n{Y_i-\beta_1 \sum_{i=1}^n{X_i}}}{V}\right), \left(\frac{n}{V} + \frac{1}{V_{\beta_0}}\right)^{-1}\right)\\ \beta_1|\beta_0, V, \boldsymbol{y} &\sim& N\left( \left(\frac{\sum_{i=1}^n{X_i^2}}{V} + \frac{1}{V_{\beta_1}}\right)^{-1} \left(\frac{\sum_{i=1}^n{(Y_i-\beta_0) X_i}}{V}\right), \left(\frac{\sum_{i=1}^n{X_i^2}}{V} + \frac{1}{V_{\beta_1}}\right)^{-1} \right)\\ V|\beta_1, \beta_0, \boldsymbol{y} &\sim& GI\left( \frac{n}{2} + a, b+\frac{1}{2}\sum_{i=1}^n\left(Y_i-\beta_0-\beta_1 X_i\right)^2\right) \end{eqnarray*}\]

#########################################
#Yi ~ N(beta0+beta1*Xi, V), i=1,...,n
#beta0 ~ N(0, Vbeta0)
#beta1 ~ N(0, Vbeta1)
#V  ~ GI(a,b)
#########################################

#Gerando uma amostra
n       = 1000
beta0   = 2
beta1   = 1
V       = 1
x       = rnorm(n, 0, 1) 
mu      = beta0+beta1*x
y       = rnorm(n, mu, sqrt(V))
hist(y, freq=F, main="", cex.lab=1.4, cex.axis=1.4, bty="n", ylab="densidade")

#note que Yi não é identicamente distribuida e por isso esta distribuição não
#tem uma media constante e nem uma variancia constante
plot(x,y)

hist(y-mu, freq=F, main="", cex.lab=1.4, cex.axis=1.4, bty="n", ylab="densidade")
curve(dnorm(x, 0, sqrt(V)), add=T, col="red", lwd=2)

#Supondo beta0, beta1 e V desconhecidos
#Amostrando pelo Amostrador de Gibbs

#Definindo alguns parâmetros
nite        = 1000
beta0_post  = NULL
beta1_post  = NULL
V_post      = NULL

#Definindo os hiperparâmetros (=parâmetros da distribuição a priori)
Vbeta0  = 10000
Vbeta1  = 10000
a       = 2
b       = 1

#Inicializando a cadeia
beta0_post[1] = 0
beta1_post[1] = 0
V_post[1]     = 10

#Executando o MCMC
for(ite in 2:nite)
{
  #Amostrando beta0 da DCCP 
  variancia         = 1/(n/V_post[ite-1] + 1/Vbeta0)
  media             = variancia*(sum(y-beta1_post[ite-1]*x)/V_post[ite-1])
  beta0_post[ite]   = rnorm(1, media, sqrt(variancia))
  
  #Amostrando beta1 da DCCP 
  variancia         = 1/(sum(x^2)/V_post[ite-1] + 1/Vbeta1)
  media             = variancia*(sum((y-beta0_post[ite])*x)/V_post[ite-1])
  beta1_post[ite]   = rnorm(1, media, sqrt(variancia))

  #Amostrando V da DCCP
  V_post[ite]   = 1/rgamma(1, n/2+a , b+0.5*sum((y-beta0_post[ite]-beta1_post[ite]*x)^2))
}

#Analisando a convergencia
par(mfrow=c(2,3))
##Traço da cadeia com toda a amostra obtida
ts.plot(beta0_post,main="traço da cadeia",ylab=expression(beta[0]),xlab="iteração")
abline(h=beta0, lwd=2, col="red", lty=3)
ts.plot(beta1_post,main="traço da cadeia",ylab=expression(beta[1]),xlab="iteração")
abline(h=beta1, lwd=2, col="red", lty=3)
ts.plot(V_post,main="traço da cadeia",ylab="V",xlab="iteração")
abline(h=V, lwd=2, col="red", lty=3)
##Traço da cadeia removendo o período de aquecimento
ind = seq(50,nite,by=1)
ts.plot(beta0_post[ind],main="traço da cadeia",ylab=expression(beta[0]),xlab="iteração")
abline(h=beta0, lwd=2, col="red", lty=3)
ts.plot(beta1_post[ind],main="traço da cadeia",ylab=expression(beta[1]),xlab="iteração")
abline(h=beta1, lwd=2, col="red", lty=3)
ts.plot(V_post[ind],main="traço da cadeia",ylab="V",xlab="iteração")
abline(h=V, lwd=2, col="red", lty=3)

par(mfrow=c(2,3))
#com toda a amostra
acf(beta0_post,main="Autocorrelação da cadeia", 
    ylab=expression(beta[0]), xlab="defasagem")
acf(beta1_post,main="Autocorrelação da cadeia", 
    ylab=expression(beta[1]), xlab="defasagem")
acf(V_post,main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")
#removendo o período de aquecimento
ind = seq(50,nite,by=1)
acf(beta0_post[ind],main="Autocorrelação da cadeia", 
    ylab=expression(beta[0]), xlab="defasagem")
acf(beta1_post[ind],main="Autocorrelação da cadeia", 
    ylab=expression(beta[1]), xlab="defasagem")
acf(V_post[ind],main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")

#Obtendo estimativas pontuais e intervalares com a amostra obtida
#após a remoção do período de aquecimento e o espaçamento
par(mfrow=c(1,3))
hist(beta0_post[ind],main="Histograma", freq=F,
     xlab=expression(beta[0]), ylab="densidade", bty="n", lwd=2)
abline(v=beta0, lwd=2, col="red", lty=3)
abline(v=quantile(beta0_post[ind],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(beta0_post[ind],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(beta0_post[ind]), lwd=2, lty=4, col="purple")

hist(beta1_post[ind],main="Histograma",freq=F, 
     xlab=expression(beta[1]), ylab="densidade", bty="n", lwd=2)
abline(v=beta1, lwd=2, col="red", lty=3)
abline(v=quantile(beta1_post[ind],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(beta1_post[ind],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(beta1_post[ind]), lwd=2, lty=4, col="purple")

hist(V_post[ind],main="Histograma", freq=F,
     xlab="V", ylab="densidade", bty="n", lwd=2)
abline(v=V, lwd=2, col="red", lty=3)
abline(v=quantile(V_post[ind],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(V_post[ind],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(V_post[ind]), lwd=2, lty=4, col="purple")

Exercício 3: - Calcular as distribuições condicionais completas a posteriori de \(\boldsymbol{\beta} = (\beta_0, \beta_1)\) e \(V\). E implementar o algoritmo Amostrador de Gibbs para obter uma amostra da distribuição a posteriori dada acima.

Reescreveremos o modelo usando a forma matricial: \[\begin{eqnarray*} \boldsymbol{Y} &\sim& N_n\left( \boldsymbol{X} \boldsymbol{\beta}, V I_n\right)\\ \boldsymbol{\beta} &\sim& N_2\left( \boldsymbol{0}, \boldsymbol{\Sigma}\right) \end{eqnarray*}\] sendo \(\boldsymbol{Y} = (Y_1 \; Y_2 \; \ldots \; Y_n)^{T}\), \(\boldsymbol{X} = \left(\begin{array}{cc} 1 & X_{12}\\ 1 & X_{22} \\ \vdots& \\ 1 & X_{n2}\end{array} \right)\), \(I_n\) a matriz identidade de ordem \(n\), \(\boldsymbol{0} = (0 \; 0)^T\), \(\boldsymbol{\Sigma}=\left(\begin{array}{cc}V_{\beta_0} & 0 \\ 0 & V_{\beta_1}\end{array}\right)\).

Logo, a distribuição a posteriori de \(\boldsymbol{\theta}=(\boldsymbol{\beta}, V)\) pode ser reescrita da seguinte forma: \[p(\boldsymbol{\theta}|\boldsymbol{y}) \propto V^{-\left(\frac{n}{2}+a+1\right)} \exp\left\{-\frac{1}{2V}\left(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta}\right) -\frac{b}{V} -\frac{1}{2}\left(\boldsymbol{\beta}^T\boldsymbol{\Sigma}^{-1}\boldsymbol{\beta}\right)\right\}.\]

Então as distribuições condicionais completas a posteriori são

\[\begin{eqnarray*} \boldsymbol{\beta}|V, \boldsymbol{y} &\sim& N_2\left( \left(\frac{\boldsymbol{X}^T\boldsymbol{X}}{V} + \boldsymbol{\Sigma}^{-1}\right)^{-1} \left(\frac{\boldsymbol{X}^T\boldsymbol{Y}}{V}\right), \left(\frac{\boldsymbol{X}^T\boldsymbol{X}}{V} + \boldsymbol{\Sigma}^{-1}\right)^{-1}\right)\\ V|\boldsymbol{\beta}, \boldsymbol{y} &\sim& GI\left( \frac{n}{2} + a, b+\frac{1}{2}\sum_{i=1}^n\left(Y_i-\beta_0-\beta_1 X_i\right)^2\right) \end{eqnarray*}\]

#Amostrando pelo Amostrador de Gibbs beta0 e beta1 conjuntamente e V separadamente

#Definindo a matriz X
X     = matrix(1, n, 2)
X[,2] = x

#Definindo alguns parâmetros
nite        = 1000
beta_post   = matrix(NA, nite, 2)
V_post      = NULL

#Definindo os hiperparâmetros (=parâmetros da distribuição a priori)
Vbeta0  = 10000
Vbeta1  = 10000
Sigma   = diag(c(Vbeta0, Vbeta1))
a       = 2
b       = 1

#Inicializando a cadeia
beta_post[1,] = c(0,0)
V_post[1]     = 1

library(MASS)

#Executando o MCMC
for(ite in 2:nite)
{
  #Amostrando beta da DCCP 
  variancia         = solve(t(X)%*%X/V_post[ite-1] + solve(Sigma))
  media             = variancia %*% (t(X)%*%y/V_post[ite-1])
  beta_post[ite,]    = mvrnorm(1, media, variancia)
  
  #Amostrando V da DCCP
  V_post[ite]   = 1/rgamma(1, n/2+a , b+0.5*t(y-X%*%beta_post[ite,])%*%(y-X%*%beta_post[ite,]))
}

#Analisando a convergencia
par(mfrow=c(2,3))
##Traço da cadeia removendo o período de aquecimento
ind = seq(50,nite,by=1)
ts.plot(beta_post[ind,1],main="traço da cadeia",ylab=expression(beta[0]),xlab="iteração")
abline(h=beta0, lwd=2, col="red", lty=3)
ts.plot(beta_post[ind,2],main="traço da cadeia",ylab=expression(beta[1]),xlab="iteração")
abline(h=beta1, lwd=2, col="red", lty=3)
ts.plot(V_post[ind],main="traço da cadeia",ylab="V",xlab="iteração")
abline(h=V, lwd=2, col="red", lty=3)
#removendo o período de aquecimento
acf(beta_post[ind,1],main="Autocorrelação", 
    ylab=expression(beta[0]), xlab="defasagem")
acf(beta_post[ind,2],main="Autocorrelação", 
    ylab=expression(beta[1]), xlab="defasagem")
acf(V_post[ind],main="Autocorrelação", 
    ylab=expression(V), xlab="defasagem")

#Obtendo estimativas pontuais e intervalares com a amostra obtida
#após a remoção do período de aquecimento e o espaçamento
par(mfrow=c(1,3))
hist(beta_post[ind,1],main="Histograma", freq=F,
     xlab=expression(beta[0]), ylab="densidade", bty="n", lwd=2)
abline(v=beta0, lwd=2, col="red", lty=3)
abline(v=quantile(beta_post[ind,1],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(beta_post[ind,1],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(beta_post[ind,1]), lwd=2, lty=4, col="purple")

hist(beta_post[ind,2],main="Histograma", freq=F,
     xlab=expression(beta[1]), ylab="densidade", bty="n", lwd=2)
abline(v=beta1, lwd=2, col="red", lty=3)
abline(v=quantile(beta_post[ind,2],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(beta_post[ind,2],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(beta_post[ind,2]), lwd=2, lty=4, col="purple")

hist(V_post[ind],main="Histograma", freq=F,
     xlab="V", ylab="densidade", bty="n", lwd=2)
abline(v=V, lwd=2, col="red", lty=3)
abline(v=quantile(V_post[ind],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(V_post[ind],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(V_post[ind]), lwd=2, lty=4, col="purple")

3 Exemplo 3

O Amostrador de Gibbs pode ser usada para amostrar de qualquer distribuição de interesse e não somente para distribuições a posteriori. Sendo assim, suponha que queiramos obter uma amostra de \(\theta_1\) e \(\theta_2\) sabendo que \[\left(\begin{array}{c} \theta_1 \\ \theta_2\end{array}\right) \sim N_2\left( \left(\begin{array}{c} \mu_1 \\ \mu_2\end{array}\right) , \left(\begin{array}{c c} V_{11} & V_{12}\\ V_{21} & V_{22}\end{array}\right) \right)\] sendo \(\mu_1, \mu_2, V_{11}, V_{12}, V_{21}, V_{22}\) conhecidos.

3.1 Amostrando da distribuição exata

Pode-se mostrar que, para \(l \neq j\), \[\begin{eqnarray*} \theta_l|\theta_j &\sim& N(\mu_l + V_{ij} V_{jj}^{-1} (\theta_j - \mu_j) \quad , \quad V_{ll} - V_{lj}^2V_{jj}^{-1})\\ \theta_j &\sim& N(\mu_j, V_{jj}) \end{eqnarray*}\] Desta forma, podemos decompor a distribuição conjunta de \((\theta_1, \theta_2)\) no produto de normais da seguinte forma: \[\begin{eqnarray*} p(\theta_1, \theta_2) = p(\theta_1|\theta_2) p(\theta_2) \end{eqnarray*}\] Logo, para obter uma amostra da conjunta, basta amostrar \(\theta_2\) da distribuição normal com média \(\mu_2\) e variância \(V_{22}\) e depois amostrar \(\theta_1\) com média \(\mu_1 + V_{12} V_{22}^{-1} (\theta_2 - \mu_2)\) e variância \(V_{11} - V_{12}^2V_{22}^{-1}\).

#Definindo alguns parâmetros
mu    = c(2, 1)
V     = matrix(c(1, 0.7, 0.7, 1),2,2)

#Definindo os parâmetros do MCMC
n               = 1000
theta_amostra   = matrix(NA, n, 2)

#Obtendo uma amostra
for(ite in 1:n){
  theta_amostra[ite,1]  = rnorm(1, mu[1], sqrt(V[1,1]))
  theta_amostra[ite,2]  = rnorm(1, mu[2]+V[2,1]/V[1,1]*(theta_amostra[ite,1] - mu[1]),
                                sqrt(V[2,2]-V[2,1]^2/V[1,1]))
}

par(mfrow=c(1,2))
hist(theta_amostra[,1],main="Histograma", freq=F,
     xlab=expression(theta[1]), ylab="densidade", bty="n", lwd=2)
hist(theta_amostra[,2],main="Histograma", freq=F,
     xlab=expression(theta[2]), ylab="densidade", bty="n", lwd=2)

3.2 Amostrando usando o pacote do R

library(MASS)
theta_amostra2   = mvrnorm(n, mu, V)

par(mfrow=c(1,2))
hist(theta_amostra2[,1],main="Histograma",freq=F, 
     xlab=expression(theta[1]), ylab="densidade", bty="n", lwd=2)
hist(theta_amostra2[,2],main="Histograma",freq=F, 
     xlab=expression(theta[2]), ylab="densidade", bty="n", lwd=2)

3.3 Amostrando pelo Amostrador de Gibbs

Precisamos calcular as distribuições condicionais completas: \[\begin{eqnarray*} \theta_1| \theta_2 &\sim& N(\mu_1 + V_{12} V_{22}^{-1} (\theta_2 - \mu_2) \quad , \quad V_{11} - V_{12}^2V_{22}^{-1})\\ \theta_2| \theta_1 &\sim& N(\mu_2 + V_{21} V_{11}^{-1} (\theta_1 - \mu_1) \quad , \quad V_{22} - V_{21}^2V_{11}^{-1}) \end{eqnarray*}\]

#Definindo alguns parâmetros

naquecimento    = 50 #total de iterações consideradas como período de aquecimento
nespacamento    = 3 #total de iterações consideradas como espaçamento
nite          = (naquecimento + 1)+(1000-1)*nespacamento #para que a amostra tenha tamanho 1000 apos remover o aquecimento e o espacamento
theta_post    = matrix(NA, nite, 2)

#Inicializando a cadeia
theta_post[1,] = c(0,0)

library(MASS)

#Executando o MCMC
for(ite in 2:nite)
{
  #Amostrando theta da DCCP 
  theta_post[ite,1]         = rnorm(1, mu[1]+V[1,2]/V[2,2]*(theta_post[ite-1,2] - mu[2]), 
                                    sqrt(V[1,1]-V[1,2]^2/V[2,2]))
  theta_post[ite,2]         = rnorm(1, mu[2]+V[2,1]/V[1,1]*(theta_post[ite,1] - mu[1]), 
                                    sqrt(V[2,2]-V[2,1]^2/V[1,1]))
}

#Analisando a convergencia
par(mfrow=c(2,2))
##Traço da cadeia removendo o período de aquecimento
ind = seq(naquecimento+1, nite, nespacamento)
ts.plot(theta_post[ind,1],main="traço da cadeia",ylab=expression(theta[1]),xlab="iteração")
ts.plot(theta_post[ind,2],main="traço da cadeia",ylab=expression(theta[2]),xlab="iteração")
#removendo o período de aquecimento
acf(theta_post[ind,1],main="Autocorrelação", 
    ylab=expression(theta[1]), xlab="defasagem")
acf(theta_post[ind,2],main="Autocorrelação", 
    ylab=expression(theta[2]), xlab="defasagem")

par(mfrow=c(3,2))
hist(theta_amostra[,1],main="Histograma",freq=F, 
     xlab=expression(theta[1]), ylab="densidade", bty="n", lwd=2)
hist(theta_amostra[,2],main="Histograma", freq=F,
     xlab=expression(theta[2]), ylab="densidade", bty="n", lwd=2)
hist(theta_amostra2[,1],main="Histograma", freq=F,
     xlab=expression(theta[1]), ylab="densidade", bty="n", lwd=2)
hist(theta_amostra2[,2],main="Histograma", freq=F,
     xlab=expression(theta[2]), ylab="densidade", bty="n", lwd=2)
hist(theta_post[ind,1],main="Histograma", freq=F,
     xlab=expression(theta[1]), ylab="densidade", bty="n", lwd=2)
hist(theta_post[ind,2],main="Histograma", freq=F,
     xlab=expression(theta[2]), ylab="densidade", bty="n", lwd=2)

##Comparando as amostras obtidas
par(mfrow=c(1,2))
ts.plot(theta_post[ind,1],main="amostras",ylab=expression(theta[1]),xlab="")
lines(theta_amostra[,1], col="red", lwd=2)
lines(theta_amostra[,1], col="purple", lwd=2)
ts.plot(theta_post[ind,2],main="amostras",ylab=expression(theta[2]),xlab="")
lines(theta_amostra[,2], col="red", lwd=2)
lines(theta_amostra[,2], col="purple", lwd=2)

4 Exemplo 4

Seja \[\begin{eqnarray*} Y_i &\stackrel{iid}{\sim}& Pois(\lambda), \; i= 1, 2, \ldots, K \\ Y_i &\stackrel{iid}{\sim}& Pois(\phi), \; i= K+1, 2, \ldots, n. \\ \end{eqnarray*}\]

Para exemplificar o modelo acima, suponha que \(Y_i\) seja o número de pessoas doentes no dia \(i\). Ou seja, estamos considerando que para os \(K\) primeiros dias, em média, \(\lambda\) pessoas ficam doentes e para os \((n-K)\) dias seguintes, esta média muda para \(\phi\) pessoas.

Suponha que \(\lambda, \phi, K\) sejam desconhecidos. Com uma amostra \(\boldsymbol{y} = y_1, \ldots, y_n\), como podemos estimar o valor de \(K, \lambda, \phi\)?

A distribuição das observações condicionada nos parâmetros é \[\begin{eqnarray*} p(y_1, \ldots, y_n | \lambda, \phi, K) &=& p(y_1, \ldots, y_k | \lambda, K) p(y_{k+1}, \ldots, y_n | \phi, K)\\ &=& \frac{1}{\prod_{i=1}^n{y_i}}\lambda^{\sum_{i=1}^{K}{y_i}} \times\\ && \exp(-\lambda K) \phi^{\sum_{i=K+1}^{n}{y_i}} \exp(-\phi (n-K)) \end{eqnarray*}\]

Para usar a inferência bayesiana, precisamos atribuir uma distribuição para \((K, \lambda, \phi)\). Considere que estes parâmetros sejam independentes e que \[\begin{eqnarray*} \lambda &\sim& G(a_\lambda, b_\lambda)\\ \phi &\sim& G(a_\phi, b_\phi)\\ K &\sim& UnifDisc\{1, \ldots, n\}. \end{eqnarray*}\]

Portanto a distribuição conjunta de \((K, \lambda, \phi)\) é proporcional a \[\begin{eqnarray} p(\lambda, \phi, K | y_1, \ldots, y_n) &\propto& p(y_1, \ldots, y_n | \lambda, \phi, K) p(\lambda, \phi, K) \nonumber\\ &\propto& \lambda^{\sum_{i=1}^{K}{y_i}} \exp(-\lambda K) \phi^{\sum_{i=K+1}^{n}{y_i}} \exp(-\phi (n-K)) \times \nonumber\\ && \lambda^{a_\lambda-1} \exp(-b_\lambda\lambda) \phi^{a_\phi-1} \exp(-b_\phi\phi)\label{post1} \end{eqnarray}\]

Como não sabemos amostrar da distribuição acima, recorreremos aos métodos de Monte Carlo via cadeias de Markov e, mais especificamente, ao Amostrador de Gibbs. Para isso, precisamos calcular as distribuições condicionais completas a posteriori:

DCCP de \(\lambda\): \[\begin{eqnarray*} p(\lambda| K, \phi , \boldsymbol{y}) &\propto& \lambda^{\sum_{i=1}^{K}{y_i}+a_\lambda-1} \exp\left\{-\lambda (K+b_\lambda)\right\} \end{eqnarray*}\] e então \(\lambda| K, \phi , \boldsymbol{y} \sim G(\sum_{i=1}^{K}{y_i}+a_\lambda, K+b_\lambda)\).

DCCP de \(\phi\): \[\begin{eqnarray*} p(\phi| K, \lambda , \boldsymbol{y}) &\propto& \phi^{\sum_{i=K+1}^{n}{y_i}+a_\phi-1} \exp\left\{-\phi \left((n-K)+b_\phi\right)\right\} \end{eqnarray*}\] e então \(\phi| K, \lambda , \boldsymbol{y} \sim G(\sum_{i=K+1}^{n}{y_i}+a_\phi, n-K+b_\phi)\), quando \(K = \{1, 2, \ldots, n-1\}\).

DCCP de K:

\[\begin{eqnarray*} P(K| \lambda, \phi , \boldsymbol{y}) &=& \left\{ \begin{array}{cc} c \lambda^{\sum_{i=1}^{K}{y_i}} \phi^{\sum_{i=K+1}^{n}{y_i}} \exp(-\lambda K -\phi (n-K)) & \mbox{ se } K=1, 2, \ldots, n-1 \\ c \lambda^{\sum_{i=1}^{K}{y_i}} \exp(-\lambda K) & \mbox{ se } K=n \end{array}\right . \end{eqnarray*}\] sendo \(c\) uma constante. Para calcular o valor de \(c\), considere que \(q_l\) seja igual a \[\begin{eqnarray*} q_l = \left \{ \begin{array}{cc} \lambda^{\sum_{i=1}^{l}{y_i}} \phi^{\sum_{i=l+1}^{n}{y_i}} \exp(-\lambda l -\phi (n-l)) & \mbox{ se } l=1, \ldots, n-1 \\ \lambda^{\sum_{i=1}^{l}{y_i}} \exp(-\lambda l) & \mbox{ se } l=n \end{array}\right . \end{eqnarray*}\] e então \[c = \frac{1}{q_1 + q_2 + \ldots + q_n}.\]

Logo, para obter uma amostra da distribuição a posteriori dos parâmetros desconhecidos, faremos os seguintes passos:

  • Inicializo \(\boldsymbol{\theta}^{(1)} = (\lambda^{(1)}, \phi^{(1)}, K^{(1)})\).

  • Faço \(j=2\).

  • Gero \(\lambda^{(j)} \sim G(a_\lambda + \sum_{i=1}^{K^{(j-1)}}{y_i}, K^{(j-1)}+b_\lambda)\).

  • Gero \(\phi^{(j)} \sim G(a_\phi + \sum_{i=K^{(j-1)} + 1}^{n}{y_i}, n-K^{(j-1)}+b_\phi)\).

  • Calculo as probabilidades \(p_l=Pr(K=l|\boldsymbol{y}, \lambda, \phi)\) para \(l=1,2, \ldots, n\).

  • Gero \(K^{(j)}\) da distribuição discreta com probabilidade \(p_1, \ldots, p_n\).

  • Faço \(j=j+1\) e repito os 4 últimos passos anteriores até obter convergência.

Problema numérico:

Para calcular as probabilidades \(p_l\) para \(l=1,2, \ldots, n\), quando \(n\) for grande, podemos ter que esta medida assume valores muito altos ou baixos, ocasionando em problemas numéricos. Sendo assim, podemos usar o seguinte artifício: \[\begin{eqnarray*} p_l = \frac{q_l}{q_1 + q_2 + \ldots, q_n} = \exp\left( \log(q_l) - \log(q_1 + q_2 + \ldots, q_n) \right) \end{eqnarray*}\] e para calcular \(\log(q_1 + q_2 + \ldots, q_n)\) podemos fazer o seguinte \[\begin{eqnarray*} q &=& q_1+q_2 \\ \log(q) &=& \log(q_1+q_2) = \log(q_1) + \log\left(1 + q_2/q_1 \right) \\ &=& \log(q_1) + \log\left(1 + \exp(\log(q_2) - \log(q_1) ) \right) \end{eqnarray*}\]

Inicializo o contador \(l = 3\) e considero \(q = \sum_{i=1}^{l-1}{q_l}\). Desta forma, tenho que \[\begin{eqnarray*} \log(q + q_l) &=& \log(q) + \log\left(1 + \exp(\log(q_l) - \log(q) ) \right) \end{eqnarray*}\] Faço \(l = l+1\) e repito a equação anterior até obter \(\log(q_1 + q_2 + \ldots, q_n)\). Depois faço \[p_K = Pr(K=l) = \exp \left\{ \log(q_l) - \log(q_1 + q_2 + \ldots + q_n) \right \}.\]

#########################################
#Yi ~ Poi(lambda), i=1,...,K
#Yi ~ Poi(phi), i=K+1,...,n
#lambda ~ G(alambda, blambda)
#phi    ~ G(aphi, bphi)
#K      ~ UnifDisc({1,2,...,n})
#########################################

#Gerando uma amostra
n           = 100
K           = 20
lambda      = 2
phi         = 4
y           = c( rpois(K,lambda), rpois(n-K, phi) )
par(mfrow=c(1,1))
barplot(table(y))

#Supondo lambda, phi e K desconhecidos
#Amostrando pelo Amostrador de Gibbs

#Definindo alguns parâmetros
naquecimento    = 50 #total de iterações consideradas como período de aquecimento
nespacamento    = 15 #total de iterações consideradas como espaçamento
nite          = (naquecimento + 1)+(1000-1)*nespacamento #para que a amostra tenha tamanho 1000 apos remover o aquecimento e o espacamento
lambda_post   = NULL
phi_post      = NULL
K_post        = NULL

#Definindo os hiperparâmetros (=parâmetros da distribuição a priori)
alambda  = 0.5
blambda  = 0.5
aphi     = 0.5
bphi     = 0.5

#Inicializando a cadeia
lambda_post[1]  = 10
phi_post[1]     = 10
K_post[1]       = n/2

#Executando o MCMC
for(ite in 2:nite)
{
  #Amostrando lambda da DCCP 
  apost             = alambda + sum(y[1:K_post[ite-1]])
  bpost             = blambda + K_post[ite-1]
  lambda_post[ite]  = rgamma(1, apost, bpost)
  
  #Amostrando phi da DCCP 
  if(K_post[ite-1]<n){soma = sum(y[(K_post[ite-1]+1):n])} else{soma=0}
  apost           = aphi + soma
  bpost           = bphi + (n-K_post[ite-1])
  phi_post[ite]   = rgamma(1, apost, bpost)

  #Amostrando K da DCCP
  pl = NULL
  for (l in 1:n){
    t1l = sum(y[1:l])
      t2l = sum(y) - t1l
      pl[l] = (lambda_post[ite]^t1l) * exp(-l*lambda_post[ite]) * (phi_post[ite]^t2l) * exp(-(n-l)*phi_post[ite])
    }   
    K_post[ite] = sample(1:n, size=1, prob=pl)
}

#Analisando a convergencia
par(mfrow=c(2,3))
##Traço da cadeia com toda a amostra obtida
ts.plot(lambda_post,main="traço da cadeia",ylab=expression(lambda),xlab="iteração")
abline(h=lambda, lwd=2, col="red", lty=3)
ts.plot(phi_post,main="traço da cadeia",ylab=expression(phi),xlab="iteração")
abline(h=phi, lwd=2, col="red", lty=3)
ts.plot(K_post,main="traço da cadeia",ylab="K",xlab="iteração")
abline(h=K, lwd=2, col="red", lty=3)
##Traço da cadeia removendo o período de aquecimento e o espacamento
ind = seq(naquecimento,nite,by=nespacamento)
ts.plot(lambda_post[ind],main="traço da cadeia",ylab=expression(lambda),xlab="iteração")
abline(h=lambda, lwd=2, col="red", lty=3)
ts.plot(phi_post[ind],main="traço da cadeia",ylab=expression(phi),xlab="iteração")
abline(h=phi, lwd=2, col="red", lty=3)
ts.plot(K_post[ind],main="traço da cadeia",ylab="K",xlab="iteração")
abline(h=K, lwd=2, col="red", lty=3)

par(mfrow=c(2,3))
#com toda a amostra
acf(lambda_post,main="Autocorrelação da cadeia", 
    ylab=expression(lambda), xlab="defasagem")
acf(phi_post,main="Autocorrelação da cadeia", 
    ylab=expression(phi), xlab="defasagem")
acf(K_post,main="Autocorrelação da cadeia", 
    ylab="K", xlab="defasagem")
#removendo o período de aquecimento e o espaçamento
acf(lambda_post[ind],main="Autocorrelação da cadeia", 
    ylab=expression(lambda), xlab="defasagem")
acf(phi_post[ind],main="Autocorrelação da cadeia", 
    ylab=expression(phi), xlab="defasagem")
acf(K_post[ind],main="Autocorrelação da cadeia", 
    ylab="K", xlab="defasagem")

#Obtendo estimativas pontuais e intervalares com a amostra obtida
#após a remoção do período de aquecimento e o espaçamento
par(mfrow=c(1,4))
hist(lambda_post[ind],main="Histograma", freq=F,
     xlab=expression(lambda), ylab="densidade", bty="n", lwd=2)
abline(v=lambda, lwd=2, col="red", lty=3)
abline(v=quantile(lambda_post[ind],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(lambda_post[ind],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(lambda_post[ind]), lwd=2, lty=4, col="purple")

hist(phi_post[ind],main="Histograma",freq=F, 
     xlab=expression(phi), ylab="densidade", bty="n", lwd=2)
abline(v=phi, lwd=2, col="red", lty=3)
abline(v=quantile(phi_post[ind],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(phi_post[ind],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(phi_post[ind]), lwd=2, lty=4, col="purple")

barplot(table(K_post[ind])/sum(table(K_post[ind])),main="", 
     xlab="K", ylab="frequência relativa", bty="n", lwd=2)
abline(v=K, lwd=2, col="red", lty=3)

hist(K_post[ind],main="Histograma", 
     xlab=expression(K), ylab="densidade", bty="n", lwd=2, freq=F)
abline(v=K, lwd=2, col="red", lty=3)
abline(v=quantile(K_post[ind],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(K_post[ind],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(K_post[ind]), lwd=2, lty=4, col="purple")

mean(K_post[ind])
## [1] 17.526
which.max( table(K_post[ind]) )
## 14 
##  3

Exercício:

  • Alterar o tamanho da amostra \(n\) para 1000 e rodar o script acima.

  • Alterar o script acima implementando a proposta dada para resolver o problema numérico.