Amostrador de Gibbs

1 Modelo normal

Considere que a variável de interesse tenha a seguinte distribuição: \[\begin{eqnarray} Y_i|\mu, V &\stackrel{iid}\sim& N(\mu, V), \quad i=1, \ldots, n \end{eqnarray}\]

1.1 Simulando o dado

Simule um conjunto de dados desse modelo. Para isso, considere que a amostra terá tamanho 1000. Atribua algum valor para os parâmetros \(\mu\) e \(V\). Em seguida, faça um histograma da amostra gerada.

#########################################
#Yi ~ N(mu, V), i=1,...,n
#########################################

#Gerando uma amostra
n     = 1000
mu    = 10
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)

1.2 Estimando os parâmetros (1 amostra)

Suponha que os valores dos parâmetros sejam desconhecidos. E suponha a priori que \[\begin{eqnarray} \mu &\sim& N(0,1) \\ V &\sim& GI(a, b) \end{eqnarray}\]

Note que a distribuição a posteriori de \((\mu, V)\) é proporcional a \[\begin{eqnarray} p(\mu, V|y_1, \ldots, y_n) &\propto& V^{-n/2-a-1} \exp\left\{-\frac{1}{2V}(y_i - \mu)^2\right\} \exp\left\{-\frac{b}{V}\right\} \exp\left\{-\frac{\mu^2}{2}\right\}. \end{eqnarray}\]

Obtenha uma amostra da distribuição a posteriori de \(\mu\) e \(V\) usando o Amostrador de Gibbs. Para isso, calcule a distribuição condicional completa a posteriori de \(\mu\) e a DCCP de \(V\).

#########################################
# Amostrador de Gibbs
#########################################

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

#Supondo mu desconhecido

#Definindo alguns parâmetros
nite        = 2000
amostra_mu  = NULL
amostra_V   = NULL

#Definindo os hiperparâmetros (=parâmetros da distribuição a priori)
a   = 0.1
b   = 0.1

#Inicializando a cadeia
amostra_mu[1] = 0
amostra_V[1]  = 10

#Executando o MCMC
for(ite in 2:nite)
{
  #amostra_V[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/amostra_V[ite-1] + 1)
  media             = variancia*(sum(y)/amostra_V[ite-1])
  amostra_mu[ite]   = rnorm(1, media, sqrt(variancia))
  
  #amostra_mu[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
  amostra_V[ite]   = 1/rgamma(1, n/2+a , 
              b+0.5*sum((y-amostra_mu[ite])^2))
}

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

par(mfrow=c(2,2))
#com toda a amostra
acf(amostra_mu,main="Autocorrelação da cadeia", 
    ylab=expression(mu), xlab="defasagem")
acf(amostra_V,main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")
#removendo o período de aquecimento
acf(amostra_mu[ind],main="Autocorrelação da cadeia", 
    ylab=expression(mu), xlab="defasagem")
acf(amostra_V[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(amostra_mu[ind],main="Histograma", 
     xlab=expression(mu), ylab="densidade", bty="n", lwd=2)
abline(v=mu, lwd=2, col="red", lty=3)
abline(v=quantile(amostra_mu[ind],c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(amostra_mu[ind],c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(amostra_mu[ind]), lwd=2, lty=4, col="purple")

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

1.3 Estimando os parâmetros (3 amostras)

Obtenha mais 2 amostras da distribuição a posteriori de \(\mu\) e \(V\) usando o Amostrador de Gibbs.

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

#Definindo alguns parâmetros
amostra_mu_2 = amostra_mu_3 = NULL
amostra_V_2  = amostra_V_3 = NULL

#Inicializando a cadeia
amostra_mu_2[1] = -100
amostra_V_2[1]  = 100

amostra_mu_3[1] = 100
amostra_V_3[1]  = 50

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

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

ts.plot(amostra_V,main="traço da cadeia",ylab="V",xlab="amostra") #1a cadeia
abline(h=V, lwd=2, col="red", lty=3) #valor verdadeiro
ts.plot(amostra_V_2,main="traço da cadeia",ylab="V",xlab="amostra")#2a cadeia
abline(h=V, lwd=2, col="red", lty=3)
ts.plot(amostra_V_3,main="traço da cadeia",ylab="V",xlab="amostra")#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(amostra_mu[ind],main="traço da cadeia",ylab=expression(mu),xlab="amostra")
lines(ind,amostra_mu_2[ind], col="purple")
lines(ind,amostra_mu_3[ind], col="green")
abline(h=mu, lwd=2, col="red", lty=3)
ts.plot(amostra_V[ind],main="traço da cadeia",ylab="V",xlab="amostra")
lines(ind,amostra_V_2[ind], col="purple")
lines(ind,amostra_V_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(amostra_mu,main="Autocorrelação da cadeia", 
    ylab=expression(mu), xlab="defasagem")
acf(amostra_mu_2,main="Autocorrelação da cadeia", 
    ylab=expression(mu), xlab="defasagem")
acf(amostra_mu_3,main="Autocorrelação da cadeia", 
    ylab=expression(mu), xlab="defasagem")
acf(amostra_V,main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")
acf(amostra_V_2,main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")
acf(amostra_V_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(c(amostra_mu[ind],amostra_mu_2[ind],amostra_mu_3[ind]),main="Histograma", 
     xlab=expression(mu), ylab="densidade", bty="n", lwd=2)
abline(v=mu, lwd=2, col="red", lty=3)
abline(v=quantile(c(amostra_mu[ind],amostra_mu_2[ind],amostra_mu_3[ind]),c(0.025, 0.975)), lwd=2, col="blue", lty=3)
abline(v=quantile(c(amostra_mu[ind],amostra_mu_2[ind],amostra_mu_3[ind]),c(0.5)), lwd=2, col="blue", lty=4)
abline(v=mean(c(amostra_mu[ind],amostra_mu_2[ind],amostra_mu_3[ind])), lwd=2, lty=4, col="purple")

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

##Avaliando a convergência usando o pacote CODA

#inicializando o pacote
library(coda)

#criando um objeto do tipo MCMC
cadeia1 = mcmc(data=data.frame(mu=amostra_mu[ind], V=amostra_V[ind]))
cadeia2 = mcmc(data=data.frame(mu=amostra_mu_2[ind], V=amostra_V_2[ind]))
cadeia3 = mcmc(data=data.frame(mu=amostra_mu_3[ind], V=amostra_V_3[ind]))
summary(cadeia1)
## 
## Iterations = 1:1901
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 1901 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##      Mean      SD  Naive SE Time-series SE
## mu 9.9736 0.03124 0.0007164      0.0007164
## V  0.9586 0.04306 0.0009877      0.0009877
## 
## 2. Quantiles for each variable:
## 
##      2.5%    25%    50%    75%  97.5%
## mu 9.9119 9.9531 9.9740 9.9947 10.035
## V  0.8754 0.9297 0.9564 0.9874  1.045
summary(cadeia2)
## 
## Iterations = 1:1901
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 1901 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##      Mean      SD  Naive SE Time-series SE
## mu 9.9729 0.02997 0.0006874      0.0006874
## V  0.9577 0.04223 0.0009685      0.0009685
## 
## 2. Quantiles for each variable:
## 
##      2.5%    25%    50%   75%  97.5%
## mu 9.9131 9.9530 9.9730 9.993 10.031
## V  0.8781 0.9282 0.9571 0.984  1.045
summary(cadeia3)
## 
## Iterations = 1:1901
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 1901 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##      Mean      SD  Naive SE Time-series SE
## mu 9.9738 0.03120 0.0007156      0.0007156
## V  0.9591 0.04457 0.0010223      0.0010223
## 
## 2. Quantiles for each variable:
## 
##      2.5%    25%    50%    75%  97.5%
## mu 9.9148 9.9532 9.9748 9.9941 10.035
## V  0.8776 0.9285 0.9572 0.9875  1.053
plot(cadeia1)

plot(cadeia2)

plot(cadeia3)

amostra_mcmc   = mcmc.list(cadeia1, cadeia2, cadeia3)
plot(amostra_mcmc)

plot(amostra_mcmc, type = "p")

summary(amostra_mcmc)
## 
## Iterations = 1:1901
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 1901 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##      Mean      SD  Naive SE Time-series SE
## mu 9.9735 0.03081 0.0004079      0.0004080
## V  0.9585 0.04330 0.0005733      0.0005734
## 
## 2. Quantiles for each variable:
## 
##      2.5%   25%    50%    75%  97.5%
## mu 9.9133 9.953 9.9739 9.9940 10.033
## V  0.8772 0.929 0.9568 0.9863  1.048

1.3.1 Gelman-Rubin

O método de Gelman-Rubin, também conhecido como fator de redução da escala potencial (PSRF), é uma medida usada para avaliar a convergência de múltiplas cadeias de Markov Chain Monte Carlo (MCMC). Ele foi proposto por Donald B. Rubin e Andrew Gelman como uma extensão do trabalho anterior de Brooks e Gelman.

A ideia fundamental do método é comparar a variabilidade intra-cadeias com a variabilidade inter-cadeias. Se as cadeias convergiram para a mesma distribuição alvo, esperamos que a variabilidade intra-cadeias seja similar à variabilidade inter-cadeias.

O PSRF é calculado da seguinte forma:

Variabilidade intra-cadeias (dentro): É a média das variâncias estimadas para cada cadeia individual. A variância é estimada considerando todas as iterações da cadeia.

Variabilidade inter-cadeias (entre): É a variância média entre as médias das cadeias. Calcula-se a média das médias de todas as cadeias e, em seguida, calcula-se a variância dessas médias.

Fator de redução da escala potencial (PSRF): É a razão da variabilidade total (intra + inter) para a variabilidade intra-cadeias. Um PSRF próximo a 1 sugere convergência.

O PSRF é calculado para cada parâmetro de interesse. Um PSRF grande (maior que 1,1, por exemplo) indica que as cadeias ainda estão convergindo. Um PSRF próximo de 1 sugere convergência.

gelman.diag(amostra_mcmc)
## Potential scale reduction factors:
## 
##    Point est. Upper C.I.
## mu          1          1
## V           1          1
## 
## Multivariate psrf
## 
## 1

1.3.2 Geweke

A medida de Geweke é outro método para avaliar a convergência de cadeias de Markov Chain Monte Carlo (MCMC). Ela compara as médias de duas subamostras de iterações da cadeia para identificar se há evidências de que as amostras finais e iniciais provêm da mesma distribuição estacionária. A abordagem de Geweke é frequentemente usada em conjunto com outras técnicas de diagnóstico de convergência.

O teste de Geweke é executado comparando estatísticas z padronizadas calculadas a partir de duas subamostras da cadeia. A estatística de teste é definida como:

\[Z = \frac{\bar{\theta}_1-\bar{\theta}_2}{\sqrt{S_1^2+S_2^2}}\] sendo \(\bar{\theta}_1\) e \(\bar{\theta}_2\) as médias das 2 subamostras e \(S_1^2\) e \(S_2^2\) as variâncias das 2 subamostras.

Se as iterações estão convergindo, a estatística de teste deve se aproximar de uma distribuição normal padrão, ou seja, \(Z \sim N(0,1)\).

Valores Z próximos de zero (\(-1,96 < Z < 1,96\)): sugere convergência para a distribuição estacionária. Valores próximos de zero indicam que as médias das subamostras inicial e final não são significativamente diferentes.

Valores Z distantes de zero (\(Z < -1.96 ou Z > 1.96\)): Indicam falta de convergência. Valores significativamente diferentes de zero sugerem que as médias das subamostras inicial e final são estatisticamente diferentes.

resultado = geweke.diag(amostra_mcmc)
resultado
## [[1]]
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##      mu       V 
## -0.5868  0.7321 
## 
## 
## [[2]]
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##     mu      V 
##  0.573 -1.816 
## 
## 
## [[3]]
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##     mu      V 
## 1.5325 0.6426
geweke.plot(amostra_mcmc)

Z_mu = c(resultado[[1]][["z"]][["mu"]], resultado[[2]][["z"]][["mu"]], resultado[[3]][["z"]][["mu"]])
#Z_mu2 = unlist(sapply(resultado, function(res) res[["z"]][["mu"]]))
Z_mu
## [1] -0.5867952  0.5729842  1.5324756
#Z_mu2
Z_V = c(resultado[[1]][["z"]][["V"]], resultado[[2]][["z"]][["V"]], resultado[[3]][["z"]][["V"]])
Z_V
## [1]  0.7321190 -1.8162394  0.6425647
if (any(Z_mu > 1.96)) {
  cat("Possível falta de convergência para mu. Considere mais iterações ou ajustes no modelo.\n")
} else {
  cat("Convergência aparente para mu.\n")
}
## Convergência aparente para mu.
if (any(Z_V > 1.96)) {
  cat("Possível falta de convergência para V. Considere mais iterações ou ajustes no modelo.\n")
} else {
  cat("Convergência aparente para V.\n")
}
## Convergência aparente para V.