This notebook was developed as a part of the Graduate Bayesian Econometrics course at Insper in 2020. The question to be solved can be found here http://hedibert.org/wp-content/uploads/2020/02/takehome-midterm-exam-bayesianeconometrics-2020.pdf
Parte A¶
Nessa questão queremos que marginalizando a distribuição de $(\epsilon_i,\lambda_i)| \sigma^2$ em relação a $\lambda_i$, chegamos que $\epsilon_i | \sigma^2 \sim \text{Laplace}(0,\sigma)$. Ou seja, queremos mostrar que:
\begin{align} \\ p(\epsilon_i|\sigma^2)&=\int p(\epsilon_i|\lambda_i,\sigma^2)p(\lambda_i| \sigma^2)\mathrm{d}\lambda_i \\ \\ &=\int \text{N}(\epsilon_i|0,\lambda_i \sigma^2)\cdot\text{Exp}(\lambda_i|1/2)\mathrm{d}\lambda_i \\ \\ &= \text{Laplace}(\epsilon_i|0,\sigma) \\ \\ \end{align}Sendo que, por um pequeno abuso de notação, $\text{N}(\epsilon_i|0,\lambda_i \sigma^2)$ denota a função densidade condicional de $\epsilon_i|\lambda_i, \sigma^2$, $\text{Exp}(\lambda_i|1/2)$ denota a função de densidade de $\lambda_i$ e $\text{Laplace}(\epsilon_i|0,\sigma)$ a distribuição condicional de $\epsilon_i|\sigma$, que é onde queremos chegar. Vamos para as contas:
\begin{align} \\ p(\epsilon_i|\sigma^2)&=\int \text{N}(\epsilon_i|0,\lambda_i \sigma^2)\cdot\text{Exp}(\lambda_i|1/2)\mathrm{d}\lambda_i \\ \\ &= \int_{0}^{\infty} \frac{1}{\sqrt{2\pi\lambda_i\sigma^2}}\text{exp}\left ( -\frac{1}{2}\frac{\epsilon_i^2}{\lambda_i \sigma^2} \right ) \cdot \frac{1}{2} \text{exp}\left ( -\frac{\lambda_i}{2} \right ) \mathrm{d}\lambda_i \\ \\ &\propto \int_{0}^{\infty} \frac{1}{\sqrt{\lambda_i}}\text{exp}\left [ -\frac{1}{2}\left ( \frac{\epsilon_i^2}{\lambda_i \sigma^2} +\lambda_i \right ) \right ]\mathrm{d}\lambda_i \propto \star \\ \\ \end{align}Aqui faremos uma substituição de variáveis, tomando $\theta=\sqrt{\lambda_i}$. Logo temos o seguinte diferencial $\mathrm{d}\lambda_i=2\theta \mathrm{d}\theta$ e, consequentemente:
\begin{align} \\ \star \propto \int_{0}^{\infty}\text{exp}\left \{ -\frac{1}{2}\left [ \left (\frac{\epsilon}{\sigma} \right )^2 \theta^{-2} + \theta^2\right ] \right \}\mathrm{d}\theta \propto \text{exp}\left ( -\frac{|\epsilon|}{\sigma} \right ) \\ \\ \end{align}Sendo que a última proporcionalidade é dado pela dica do enunciado. É possível perceber que chegamos ao núcleo da $\text{Laplace}(0,\sigma)$, concluindo nossa prova.
Parte B¶
Nesta parte, o nosso objetivo é obter as distribuições condicionais completas para todos os parâmetros, pois assim conseguiremos aplicar a técnica Amostrador de Gibbs a fim de amostrador a distribuição a posteriori dos parâmetros do nosso modelo. Adotaremos as seguintes distirbuições a priori para os nossos parâmetros:
$\beta \sim \text{N}(\beta_0, V_0)$
$\sigma^2 \sim \text{IG} \left ( \nu_0/2, \sigma^2_0 \nu_0/2 \right )$
$\lambda_1, ...,\lambda_n \overset{iid}{\sim} \text{Exp}(1/2)$
Primeiramente, veja que pelo Teorema de Bayes, a distribuição a posteriori dos parâmetros a proporcional à distribuição conjunta dos parâmetros e da nossa variável resposta, dado as variáveis explicativas:
\begin{align} \\ &p(\beta, \sigma^2, \lambda_1, ...,\lambda_n | y, X)\propto \\ \\ & \propto p(\beta, \sigma^2, \lambda_1, ...,\lambda_n, y| X)\\ \\ &=p(y|X,\beta, \sigma^2, \lambda_1, ...,\lambda_n)\cdot p(\beta) \cdot p(\sigma^2) \cdot \prod_{i=1}^{n}p(\lambda_i) \\ \\ &=\text{N}(y|X\beta, \sigma^2 \Lambda)\cdot\text{N}(\beta|\beta_0, V_0)\cdot \text{IG}\left (\sigma^2 | \nu_0/2, \sigma^2_0 \nu_0/2 \right ) \cdot \prod_{i=1}^{n} \text{Exp}(\lambda_i|1/2) \\ \\ &\propto \left(\sigma^{2n} \prod_{i=1}^{n} \lambda_i \right)^{-\frac{1}{2}} \text{exp}\left[-\frac{1}{2\sigma^2} \left(y-X\beta\right)^\top \Lambda^{-1} \left(y-X\beta\right) \right] \cdot\\ \\ &~~ \cdot \text{exp} \left[-\frac{1}{2} \left(\beta-\beta_0\right)^\top V_0^{-1} \left(\beta-\beta_0\right) \right] \cdot (\sigma^2)^{-(\nu_0/2+1)} \text{exp} \left(- \frac{ \sigma^2_0 \nu_0/2}{\sigma^2} \right) \cdot \\ \\ &~~ \cdot \text{exp} \left(- \frac{1}{2} \sum_{i=1}^{n} \lambda_i \right) \end{align}Sendo que $\Lambda=\text{diag}(\lambda_1,...,\lambda_n)$. Primeiramente vamos chegar a distirbuição condicional completa de $\beta$:
\begin{align} \\ &p(\beta |\sigma^2, \lambda_1, ...,\lambda_n, y, X) \propto \\ \\ &\propto \text{exp}\left[-\frac{1}{2\sigma^2} \left(y-X\beta\right)^\top \Lambda^{-1} \left(y-X\beta\right) \right] \cdot \text{exp} \left[-\frac{1}{2} \left(\beta-\beta_0\right)^\top V_0^{-1} \left(\beta-\beta_0\right) \right] \\ \\ &\propto \text{exp}\left\{-\frac{1}{2}\left[-2\beta^\top \left(\frac{X^\top \Lambda^{-1} y}{\sigma^2} + V_0^{-1} \beta_0 \right)+\beta^\top \left(\frac{X^\top \Lambda^{-1} X}{\sigma^2} + V_0^{-1}\right)\beta \right]\right\} \\ \\ &\propto \text{exp}\left\{-\frac{1}{2}\left[-2\beta^\top V_1^{-1}\beta_1+\beta^\top V_1^{-1}\beta \right]\right\} \Rightarrow \beta |\sigma^2, \lambda_1, ...,\lambda_n, y, X \sim \text{N}(\beta_1,V_1) \end{align}Sendo que $V_1=\left(\frac{X^\top \Lambda^{-1} X}{\sigma^2} + V_0^{-1}\right)^{-1}$ e $\beta_1=V_1 \left(\frac{X^\top \Lambda^{-1} y}{\sigma^2} + V_0^{-1} \beta_0 \right)$. Agora que já temos a condicional completa para $\beta$, vamos obter a condicional completa para $\sigma^2$:
\begin{align} \\ &p(\sigma^2 |\beta, \lambda_1, ...,\lambda_n, y, X) \propto \\ \\ &\propto \left(\sigma^{2n} \right)^{-\frac{1}{2}} \text{exp}\left[-\frac{1}{2\sigma^2} \left(y-X\beta\right)^\top \Lambda^{-1} \left(y-X\beta\right) \right] \cdot (\sigma^2)^{-(\nu_0/2+1)} \text{exp} \left(- \frac{ \sigma^2_0 \nu_0/2}{\sigma^2} \right) \\ \\ &=(\sigma^2)^{-\left(\frac{\nu_0+n}{2}+1\right)} \text{exp} \left[\frac{1}{\sigma^2} \left(-\frac{1}{2} \left(y-X\beta\right)^\top \Lambda^{-1} \left(y-X\beta\right) -\sigma_0^2 \frac{\nu_0}{2} \right) \right] \\ \\ &=(\sigma^2)^{-\left(\frac{\nu_1}{2}+1\right)} \text{exp} \left(- \frac{ \sigma^2_1 \nu_1/2}{\sigma^2} \right) \Rightarrow \sigma^2 |\beta, \lambda_1, ...,\lambda_n, y, X \sim \text{IG} \left(\frac{\nu_1}{2}, \sigma^2_1 \nu_1/2\right) \end{align}Sendo que $\nu_1=\nu_0+n$ e $\sigma^2_1=\left(y-X\beta\right)^\top \Lambda^{-1} \left(y-X\beta\right) +\sigma_0^2\nu_0$. Agora só falta sabermos a condicional para $\lambda_i,~i=1,...,n$:
\begin{align} \\ &p( \lambda_i |\beta,\sigma^2, \lambda_1,..., \lambda_{i-1}, \lambda_{i+1}, ...,\lambda_n, y, X) \propto \\ \\ & \propto \text{exp}\left[-\frac{1}{2\sigma^2} \left(y-X\beta\right)^\top \Lambda^{-1} \left(y-X\beta\right) \right] \cdot \text{exp} \left(- \frac{1}{2} \sum_{i=1}^{n} \lambda_i \right) \\ \\ & \propto \frac{1}{\sqrt{\lambda_i}}\text{exp}\left(-\frac{1}{2\sigma^2}\sum_{i=1}^{n}\frac{(y_i-x_i^\top\beta)^2}{\lambda_i}\right)\cdot \text{exp}\left(-\frac{\lambda_i}{2} \right)\\ \\ & \propto \lambda_i^{-\frac{1}{2}} \text{exp}\left\{ -\frac{1}{2} \left[\lambda_i^{-1} \left(\frac{y_i-x_i^\top\beta}{\sigma}\right)^2 + \lambda_i \right]\right\}=g(\lambda_i,x_i,y_i,\sigma^2) \end{align}É possível deduzir que $\lambda_i^{-1}|\beta,\sigma^2, \lambda_1,..., \lambda_{i-1}, \lambda_{i+1}, ...,\lambda_n, y, X \sim \text{InvGaussian}(1,|\sigma_i/(y_i-x_i^\top\beta)|)$.
Parte C¶
Fixando uma semente:
set.seed(42)
Carregando pacotes necessários:
library(MASS)
library(LaplacesDemon)
library(asbio)
library(statmod)
Para fazer ambas as simulações, primeiramente precisamos gerar nossos dados:
n=200
β=c(0,1,2,3)
q=3
σ2=1
Λ=diag(rexp(n, rate=.5))
X=cbind(1, matrix(rnorm(3*n), c(n,3)))
y=mvrnorm(n = 1, X%*%β, σ2*Λ)
Definindo os parâmetros das distribuições a priori:
#beta
β_0=c(rep(0, q+1))
V_0=diag(rep(1, q+1))
#sigma2
ν_0=100
σ2_0=1
Simulação por Gibbs-Sampler/Metropolis-Hastings¶
Nesta primeira simulação, amostraremos $\beta$ e $\sigma^2$ de suas condicionais completas e amostraremos $\lambda_i,~i=1,...,n$ com a ajuda do algoritmo Metropolis-Hastings clássico da seguinte maneira. O algoritmo completo que utilizei, para amostrar uma cadeia, é o seguinte:
- Inicializamos $\beta^{(1)}$, $\sigma^{2{(1)}}$ e $\left\{\lambda^{(1)}\right\}$;
Para $l$ indo de de $2$ a $L$:
Amostramos $\beta^{(l)}$ de $p(\beta |\sigma^{2{(l-1)}}, \lambda_1^{(l-1)}, ...,\lambda_n^{(l-1)}, y, X)$;
Amostramos $\sigma^{2{(l)}}$ de $p(\sigma^2 |\beta^{(l)}, \lambda_1^{(l-1)}, ...,\lambda_n^{(l-1)}, y, X)$;
Para $i$ indo de $1$ a $n$:
Amostrar $\lambda_i^* \sim \text{Lognormal}(\lambda_i^{(l-1)}, .1)$, sendo que $f(.|\mu, \omega)$ representa a função densidade dessa variável aleatória;
Calcular $\alpha=\text{min}\left[1, \frac{g(\lambda_i^*,x_i,y_i,\sigma^2)}{g(\lambda_i^{(l-1)},x_i,y_i,\sigma^2)}\frac{f(\lambda_i^{(l-1)}|\lambda_i^*, .1)}{f(\lambda_i^*|,\lambda_i^{(l-1)}, .1)} \right];$
Amostrar $V \sim \text{U}([0,1])$;
Se $V\leq\alpha$:
- $\lambda_i^{(l)}=\lambda_i^*$;
Se $V>\alpha$:
- $\lambda_i^{(l)}=\lambda_i^{(l-1)}$;
Ao final, descartamos as primeiras $B$ amostras de cada um dos parâmetros como o período de 'Burn-in'.
OBS: Vamos amostrar de duas cadeias com inícios diferentes afim de comparar os resultados e estarmos mais confiantes a respeito da convergência para uma distribuição de equilíbrio.
Criando arrays para armazenar amostras:
L=50000 #Total de simulações
chains=1
βs_mh=array(rep(0, L*chains*(q+1)), c(L, chains, q+1))
σ2s_mh=array(rep(0, L*chains), c(L, chains))
λs_mh=array(rep(0, L*chains*n), c(L, chains, n))
Inicializando primeiros valores:
βs_mh[1,1,]=mvrnorm(1, c(0,1,2,3), .01*diag(q+1))
#βs_mh[1,2,]=mvrnorm(1, c(0,1,2,3), .01*diag(q+1))
σ2s_mh[1,]=array(rlnorm(chains*1, 0, .01))
λs_mh[1,,]=array(rexp(chains*n, .5), c(chains, n))
Definindo função que utilizaremos no passo do Metropolis-Hastings:
g=function(λ,β,σ2,x,y) λ^-.5 * exp(-.5*(λ + (y-x%*%β)^2/σ2 * λ^-1))
Simulando:
start.time <- Sys.time()
for(c in 1:chains){
for(l in 2:L){
#sigma2
ν_1=ν_0+n
σ2_1=((ν_0*σ2_0) + t(y-X%*%βs_mh[l-1,c,])%*%solve(diag(λs_mh[l-1,c,]))%*%(y-X%*%βs_mh[l-1,c,]))/(ν_1)
σ2s_mh[l,c]=1/rgamma(1, ν_1/2, σ2_1*ν_1/2)
#beta
V_1=solve((t(X)%*%solve(diag(λs_mh[l-1,c,]))%*%X)/σ2s_mh[l,c] + solve(V_0))
β_1=V_1%*%((t(X)%*%solve(diag(λs_mh[l-1,c,]))%*%y)/σ2s_mh[l,c] + solve(V_0)%*%β_0)
βs_mh[l,c,]=mvrnorm(n = 1, β_1, V_1)
#lambdas (M-H)
for(i in 1:n){
λ_star=rlnorm(1, λs_mh[l-1, c, i], .5)
frac1=g(λ_star,βs_mh[l, c,],σ2s_mh[l, c],X[i,],y[i])/g(λs_mh[l-1, c, i],βs_mh[l, c,],σ2s_mh[l, c],X[i,],y[i])
frac2=dlnorm(λs_mh[l-1, c, i], λ_star, .5)/dlnorm(λ_star, λs_mh[l-1, c, i], .5)
α=min(frac1*frac2, 1)
if(runif(1) < α){
λs_mh[l, c, i]=λ_star
}
else{
λs_mh[l, c, i]=λs_mh[l-1, c, i]
}
}
}
}
end.time <- Sys.time()
time_mh <- as.numeric(end.time - start.time)
Descartaremos as B primeiras amostras de cada cadeia como Burn-in:
B=45000
σ2s_mh=σ2s_mh[((B+1):L),]
βs_mh=βs_mh[((B+1):L),,]
λs_mh=λs_mh[((B+1):L),,]
Plotando traceplots do $\sigma^2$:
plot(σ2s_mh, col=2, type='l', xlab="Sigma^2", main="Traceplots do 𝜎^2")
#lines(σ2s_mh[,2], col=4)
Plotando traceplots dos Betas:
for(j in 1:(q+1)){
plot(βs_mh[,j], col=2, type='l', xlab="Beta", main=paste("Traceplots do Beta",j,sep=" "))
#lines(βs_mh[,2,j], col=4)
}
Para trazermos mais uma evidência de que nossas cadeias convergiram, olharemos para as distribuições das estatísticas $\hat{R}$, propostas por Gelman and Rubin (1992). Se as estatísticas estão próximas de 1, temos evidências de que as cadeias convergiram. Primeiramente olharemos para os betas e sigma2:
rhat_mh=c(R.hat(σ2s_mh, burn.in = 0))
for(j in 1:(q+1)){
rhat_mh=c(rhat_mh, R.hat(βs_mh[,,j], burn.in = 0))
}
rhat_mh
Tivemos valores muito próximos a 1. Agora vamos olhar os lambdas:
rhat_mh=c()
for(j in 1:n){
rhat_mh=c(rhat_mh, R.hat(λs_mh[,,j], burn.in = 0))
}
summary(rhat_mh)
boxplot(rhat_mh, main="R_hat dos lambdas")
É possível ver que temos um outlier e apenas dois pontos um pouco acima de 1. O cenário é otimista para a grande maioria dos indivíduos.
Simulação por Gibbs-Sampler (Full)¶
Nesta segunda maneira de simular, amostraremos $\beta$ e $\sigma^2$ e $\lambda_i^{-1},~i=1,...,n$ de suas condicionais completas. O algoritmo completo que utilizei é o seguinte:
- Inicializamos $\beta^{(1)}$, $\sigma^{2{(1)}}$ e $\left\{\lambda^{(1)}\right\}$;
Para $l$ indo de de $2$ a $L$:
Amostramos $\beta^{(l)}$ de $p(\beta |\sigma^{2{(l-1)}}, \lambda_1^{(l-1)}, ...,\lambda_n^{(l-1)}, y, X)$;
Amostramos $\sigma^{2{(l)}}$ de $p(\sigma^2 |\beta^{(l)}, \lambda_1^{(l-1)}, ...,\lambda_n^{(l-1)}, y, X)$;
- Amostramos $(\lambda_i^{(l)})^{-1}$ de $p(\lambda_i^{-1}|\beta,\sigma^2, \lambda_1^{(l)},..., \lambda_{i-1}^{(l)}, \lambda_{i+1}^{(l-1)}, ...,\lambda_n^{(l-1)}, y, X)$ e invertemos o valor amostrado para armazená-lo.
Ao final, descartamos as primeiras $B$ amostras de cada um dos parâmetros como o período de 'Burn-in'.
OBS: Vamos amostrar de duas cadeias com inícios diferentes afim de comparar os resultados e estarmos mais confiantes a respeito da convergência para uma distribuição de equilíbrio.
Vamos codificar uma função que amostra da Inverse Gaussian:
inv_gaussian=function(n, mu, psi){
out=c()
for(i in 1:n){
v=rnorm(1)^2
x1=mu+(mu^2*v)/(2*psi)-mu/(2*psi)*sqrt(4*mu*psi*v+mu^2*v^2)
x2=mu^2/x1
p1=mu/(mu+x1)
if(runif(1)<p1) out=c(out,x1)
else out=c(out,x2)
}
return(out)
}
Comparando densidades gerada pela minha função e pela função de um pacote do R:
sample1=inv_gaussian(10000, 1, 2)
sample2=rinvgauss(10000, 1, 2)
par(mfrow=c(1,2))
hist(sample1, main='Minha função')
hist(sample2, main='Função do R')
É possível ver que estamos conseguindo amostrar da Inverse Gaussian.
Agora vou criar arrays para armazenar amostras do MCMC:
L=50000 #Total de simulações por cadeia
chains=1
βs_gib=array(rep(0, L*chains*(q+1)), c(L, chains, q+1))
σ2s_gib=array(rep(0, L*chains), c(L, chains))
λs_gib=array(rep(0, L*chains*n), c(L, chains, n))
Inicializando primeiros valores:
βs_gib[1,1,]=mvrnorm(1, c(0,1,2,3), .01*diag(q+1))
#βs_gib[1,2,]=mvrnorm(1, c(0,1,2,3), .01*diag(q+1))
σ2s_gib[1,]=array(rlnorm(chains*1, 0, .01))
λs_gib[1,,]=array(rexp(chains*n, .5), c(chains, n))
Simulando:
start.time <- Sys.time()
for(c in 1:chains){
for(l in 2:L){
#sigma2
ν_1=ν_0+n
σ2_1=((ν_0*σ2_0) + t(y-X%*%βs_gib[l-1,c,])%*%solve(diag(λs_gib[l-1,c,]))%*%(y-X%*%βs_gib[l-1,c,]))/(ν_1)
σ2s_gib[l,c]=1/rgamma(1, ν_1/2, σ2_1*ν_1/2)
#beta
V_1=solve((t(X)%*%solve(diag(λs_gib[l-1,c,]))%*%X)/σ2s_gib[l,c] + solve(V_0))
β_1=V_1%*%((t(X)%*%solve(diag(λs_gib[l-1,c,]))%*%y)/σ2s_gib[l,c] + solve(V_0)%*%β_0)
βs_gib[l,c,]=mvrnorm(n = 1, β_1, V_1)
#lambdas
for(i in 1:n){ #rinvgauss
λs_gib[l,c,i]=1/inv_gaussian(1, 1, abs(σ2s_gib[l,c]/(y[i]-X[i,]%*%βs_gib[l,c,])))
}
}
}
end.time <- Sys.time()
time_gib <- as.numeric(end.time - start.time)
time_gib/60
time_mh/60
Descartaremos as B primeiras amostras de cada cadeia como Burn-in:
B=45000
σ2s_gib=σ2s_gib[((B+1):L),]
βs_gib=βs_gib[((B+1):L),,]
λs_gib=λs_gib[((B+1):L),,]
Plotando traceplots do $\sigma^2$:
plot(σ2s_gib[,1], col=2, type='l', xlab="Sigma^2", main="Traceplots do 𝜎^2")
lines(σ2s_gib[,2], col=4)
Plotando traceplots dos Betas:
for(j in 1:(q+1)){
plot(βs_gib[,1,j], col=2, type='l', xlab="Beta", main=paste("Traceplots do Beta",j,sep=" "))
lines(βs_gib[,2,j], col=4)
}
Para trazermos mais uma evidência de que nossas cadeias convergiram, olharemos para as distribuições das estatísticas $\hat{R}$, propostas por Gelman and Rubin (1992). Se as estatísticas estão próximas de 1, temos evidências de que as cadeias convergiram. Primeiramente olharemos para os betas e sigma2:
rhat_gib=c(R.hat(σ2s_gib, burn.in = 0))
for(j in 1:(q+1)){
rhat_gib=c(rhat_gib, R.hat(βs_gib[,,j], burn.in = 0))
}
rhat_gib
Tivemos valores muito próximos a 1. Agora vamos olhar os lambdas:
rhat_gib=c()
for(j in 1:n){
rhat_gib=c(rhat_gib, R.hat(λs_gib[,,j], burn.in = 0))
}
summary(rhat_gib)
boxplot(rhat_gib, main="R_hat dos lambdas")
Apesar de o cenário não ser muito alarmante quando utilizamos M-H para os lambdas, é claro que temos um resultado muito melhor agora.
Medidas resumo sobre das distribuições a posteriori dos parâmetros¶
Comparando as duas maneiras de amostragem¶
Afim de comparar os dois algoritmos de amostragem, utilizaremos a razão ESS$/s$ como uma medida de eficiência ao amostrar dos dois algoritmos, sendo ESS o Effective Sample Size e s os segundos que o código demorou para rodar. O ESS, que é calculado para uma cadeia específica, é dado pela fórmula abaixo (Gamerman and Lopes, 2006). Para estimar o ESS, utilizei a implementação do pacote "LaplaceDemons", que utiliza a densidade espectral na frequência zero.
\begin{align} ESS(\theta)=\frac{L}{1+2\sum_{k=1}^\infty\rho_k(\theta)} \end{align}Na realidade não nos ateremos aos ESS/s propriamente ditos, vamos olhar para a razão dos ESS/s para amostras dos dois algoritmos. Calcularemos o ESS para cada uma das cadeias, de cada parâmetro de cada um dos algoritmos e tiraremos a razão, sendo que se a razão for maior que 1, temos que o algoritmo Full Gibbs é mais eficiente ao inferir sobre um parâmetro específico em uma das duas cadeias. Começaremos por $\beta$ e $\sigma^2$:
ess_mh=c()
ess_gib=c()
for(c in 1:chains){
ess_mh=c(ess_mh, ESS(σ2s_mh[,c]))
ess_gib=c(ess_gib, ESS(σ2s_gib[,c]))
}
for(c in 1:chains){
for(j in 1:(q+1)){
ess_mh=c(ess_mh, ESS(βs_mh[,c,j]))
ess_gib=c(ess_gib, ESS(βs_gib[,c,j]))
}
}
efic_mh=ess_mh/time_mh
efic_gib=ess_gib/time_gib
summary(efic_gib/efic_mh)
É possível ver que o nosso segundo algoritmo é em média 11% menos eficiente que o primeiro algoritmo ao inferir sobre $\beta$ e $\sigma^2$. Agora nos ateremos aos lambdas:
ess_mh=c()
ess_gib=c()
for(c in 1:chains){
for(j in 1:n){
ess_mh=c(ess_mh, ESS(λs_mh[,c,j]))
ess_gib=c(ess_gib, ESS(λs_gib[,c,j]))
}
}
efic_mh=ess_mh/time_mh
efic_gib=ess_gib/time_gib
summary(efic_gib/efic_mh)
Em relação aos lambdas, temos que o agoritmo Full Gibss é muito mais eficiente.
Referências:¶
Gamerman, Dani, and Hedibert F. Lopes. Markov chain Monte Carlo: stochastic simulation for Bayesian inference. Chapman and Hall/CRC, 2006.
Gelman, A. and D. B. Rubin (1992) Inference from iterative simulation using multiple sequences (with discussion). Statistical Science, 7:457-511