Inferência Bayesiana - GET/IME/UFF
Patrícia Lusié Velozo da Costa
2023-1
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*}\]
#########################################
# Amostrador de Gibbs
#########################################
#########################################
#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] = 1
#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="amostra")
abline(h=mu, lwd=2, col="red", lty=3)
ts.plot(V_post,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
ind = seq(50,nite,by=1)
ts.plot(mu_post[ind],main="traço da cadeia",ylab=expression(mu),xlab="amostra")
abline(h=mu, lwd=2, col="red", lty=3)
ts.plot(V_post[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(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",
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",
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 duas cadeias
#########################################
#########################################
#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="amostra") #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="amostra")#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="amostra")#3a cadeia
abline(h=mu, lwd=2, col="red", lty=3)
ts.plot(V_post,main="traço da cadeia",ylab="V",xlab="amostra") #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="amostra")#2a cadeia
abline(h=V, lwd=2, col="red", lty=3)
ts.plot(V_post_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(mu_post[ind],main="traço da cadeia",ylab=expression(mu),xlab="amostra")
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="amostra")
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",
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",
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.
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.
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.