Metropolis-Hastings

1 Modelo normal

Considere que a variável de interesse tenha a seguinte distribuição: \[\begin{eqnarray} Y_i|V &\stackrel{iid}\sim& N(0, 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 o parâmetro \(V\). Em seguida, faça um histograma da amostra gerada.

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

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

Suponha que o valor do parâmetro seja desconhecido. E suponha a priori que \[\begin{eqnarray} V &\sim& GI(a, b) \end{eqnarray}\]

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

Suponha que não sabemos gerar da distribuição gama inversa e que usaremos então o Metropolis-Hastings para obter uma amostra desta distribuição a posteriori.

nite          = 5000    #total de iterações
cont          = 0       #contador para calcular a taxa de aceitacao
amostra_V     = NULL    #vetor que guardará a amostra
amostra_V[1]  = 1       #inicializando a amostra

#criando uma função para calcular a densidade da gama inversa
ldgamainversa = function(u, a, b){return(a*log(b) - log(gamma(a))-(a+1)*log(u) - b/u)}
ldgamainversa(3,2,1)
## [1] -3.62917
#comparando com uma função do R
library(invgamma)
dinvgamma(3,2,1, log=T)
## [1] -3.62917
#parâmetros da distribuição a priori
a = 0.1
b = 0.1
#parâmetros da distribuição proposta
ap = 25
bp = 25
for(ite in 2:nite){
  epsilon      = rgamma(1, ap, bp) #gerando da distribuição proposta
  
  #calculando a probabilidade de aceitar o ponto gerado
  p1           = dinvgamma(epsilon, n/2+a, b+0.5*sum(y^2),log=T) - dgamma(epsilon, ap, bp)
  p2           = dinvgamma(amostra_V[ite-1], n/2+a, b+0.5*sum(y^2), log=T) - dgamma(amostra_V[ite-1], ap, bp)
  razao         = exp( p1-p2 )
  alpha         = min(1,razao)
  #decidindo se aceitamos ou não o ponto gerado
  u             = runif(1,0,1)
  if (u<alpha){
    amostra_V[ite]  = epsilon
    cont            = cont+1
  }else{
    amostra_V[ite]  = amostra_V[ite-1]
  }
}
#taxa de aceitação
cont/nite
## [1] 0.288
#Analisando a convergencia
#x11()
par(mfrow=c(2,2))
##Traço da cadeia com toda a amostra obtida
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_V[ind],main="traço da cadeia",ylab="V",xlab="amostra")
abline(h=V, lwd=2, col="red", lty=3)
#autocorrelacao
#com toda a amostra
acf(amostra_V,main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")
#removendo o período de aquecimento
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,1))
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")

#avaliando a convergência com o CODA
library(coda)

resultado = geweke.diag(mcmc(amostra_V))
resultado
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##   var1 
## 0.4176
geweke.plot(mcmc(amostra_V))

Z_V = c(resultado[["z"]][["var1"]])
Z_V
## [1] 0.4176281
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.

1.3 Usando uma outra distribuição como proposta

nite          = 5000    #total de iterações
cont          = 0       #contador para calcular a taxa de aceitacao
amostra_V     = NULL    #vetor que guardará a amostra
amostra_V[1]  = 10       #inicializando a amostra

#parâmetros da distribuição a priori
a = 0.1
b = 0.1
#parâmetros da distribuição proposta
Vp = 0.01

library(stats) #pacote da log-normal

for(ite in 2:nite){
  epsilon      = rlnorm(1, log(amostra_V[ite-1]), sqrt(Vp))
  #calculando a probabilidade de aceitar o ponto gerado
  p1           = dinvgamma(epsilon, n/2+a, b+0.5*sum(y^2),log=T) - dlnorm(epsilon, log(amostra_V[ite-1]) , sqrt(Vp), log = T)
  p2           = dinvgamma(amostra_V[ite-1], n/2+a, b+0.5*sum(y^2),log=T) - dlnorm(amostra_V[ite-1], log(epsilon) , sqrt(Vp), log = T)
  razao         = exp( p1-p2 )
  alpha         = min(1,razao)
  #decidindo se aceitamos ou não o ponto gerado
  u             = runif(1,0,1)
  if (u<alpha){
    amostra_V[ite]  = epsilon
    cont            = cont+1
  }else{
    amostra_V[ite]  = amostra_V[ite-1]
  }
}
#taxa de aceitação
cont/nite
## [1] 0.471
#Analisando a convergencia
#x11()
par(mfrow=c(2,2))
##Traço da cadeia com toda a amostra obtida
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 = 5
ind = seq(aquecimento, nite,by=espacamento)
ts.plot(amostra_V[ind],main="traço da cadeia",ylab="V",xlab="amostra")
abline(h=V, lwd=2, col="red", lty=3)
#autocorrelacao
#com toda a amostra
acf(amostra_V,main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")
#removendo o período de aquecimento
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,1))
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")

#avaliando a convergência com o CODA
library(coda)
resultado = geweke.diag(mcmc(amostra_V))
resultado
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##   var1 
## 0.9313
geweke.plot(mcmc(amostra_V))

Z_V = c(resultado[["z"]][["var1"]])
Z_V
## [1] 0.931283
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.