1 Modelo

Seja \[\begin{eqnarray} Y_i &\sim& N\left( X_i \beta, V\right), \; i=1, \ldots, n. \end{eqnarray}\]

Gerando um conjunto de dados considerando o modelo acima:

##########################################################
# Y_i   ~ N(X_i beta, V)
##########################################################

#Gerando uma amostra de Y
n               = 1000
x               = rnorm(n)
beta            = 0.5
V               = 0.1
y               = rnorm(n, x*beta, sqrt(V))
#Analisando a amostra gerada
par(mfrow=c(2,2))
hist(y, freq=F, main="", lwd=2, cex.lab=1.4, cex.axis=1.4, ylab="densidade")
hist(y-x*beta, freq=F, main="", lwd=2, cex.lab=1.4, cex.axis=1.4, ylab="densidade")
plot(x,y)

2 Metropolis-Hastings

Suponha

\(V\) conhecido, \(\beta\) desconhecido e a priori que \[\begin{eqnarray} \beta &\sim& N(0, V_\beta), \end{eqnarray}\]
sendo \(V_\beta\) conhecido;
A distribuição a posteriori \(p(\beta|\boldsymbol{y}, V)\) é a f.d.p. da distribuição \(N\left( \left(\frac{\sum_{i=1}^n{X_i^2}}{V} + \frac{1}{V_{\beta}}\right)^{-1} \left(\frac{\sum_{i=1}^n{Y_i X_i}}{V}\right), \left(\frac{\sum_{i=1}^n{X_i^2}}{V} + \frac{1}{V_{\beta}}\right)^{-1} \right)\). Mas, para exemplificar o método, suponha que ela é desconhecida. Considere então que \[\begin{eqnarray*} p(\beta|\boldsymbol{y}, V) \propto p(\boldsymbol{y}|\beta, V) p(\beta) \end{eqnarray*}\] sendo \(p(\boldsymbol{y}|\beta, V)\) a f.d.p. da \(N(X_i \beta, V)\) e \(p(\beta)\) a f.d.p. da \(N(0, V_\beta)\).
Usaremos então o algoritmo de Metropolis-Hastings para obter uma amostra desta distribuição.

Algoritmo:

Escolha um valor arbitrário para \(\beta^{(1)}\);
Inicializa-se o contador de iteração \(t=2\);
Gera-se \(\beta^{(p)} \sim q(\beta^{(p)}|\beta^{(t-1)})\). Aceita-se o ponto gerado com probabilidade \(\alpha = \min\left\{1,\displaystyle\frac{p(\beta^{(p)})}{q(\beta^{(p)}|\beta^{(t-1)})}\frac{q(\beta^{(t-1)}|\beta^{(p)})}{p(\beta^{(t-1)})} \right\}\), onde \(p(\cdot)\) é a distribuição de interesse; Para decidir sobre a aceitação:
1. calcule \(\alpha\) e gere \(u \sim U(0,1)\);
  1. se \(u<\alpha\), aceite \(\beta^p\) fazendo \(\beta^{(t)}=\beta^p\); caso contrário, rejeite e faça \(\beta^{(t)}=\beta^{(t-1)}\);
Altera-se o contador de \(t\) para \(t+1\);
Repete-se os itens 3 e 4 até que a convergência seja obtida.

2.1 Distribuição proposta: passeio aleatório

#Supondo beta desconhecido, V conhecido 
#supondo que não sabemos gerar da distribuição a posteriori de beta 
#Amostrando beta pelo Metropolis-Hastings

#Definindo alguns parâmetros
naquecimento    = 100 #total de iterações consideradas como período de aquecimento
nespacamento    = 5 #total de iterações consideradas como espaçamento
nite            = (naquecimento + 1)+(1000-1)*nespacamento #para que a amostra tenha 
#nite            = 1500

beta_post   = NULL #inicializando vetor para guardar a amostra da distribuição a posteriori
Vbeta       = 1000 #hiperparâmetro

#Inicializando a cadeia
beta_post[1] = 0

#Inicializando o contador para saber sobre a taxa de aceitação
cont = 0

# Usaremos a seguinte distribuição proposta:
# beta^p ~ N(beta^{t-1}, Vprop)
Vprop         = 0.001

#Executando o MCMC
for(ite in 2:nite)
{
    beta_prop     = rnorm(1, beta_post[ite-1], sqrt(Vprop))
    razao.l.num   = sum(dnorm(y, x*beta_prop, sqrt(V), log=T)) + dnorm(beta_prop, 0, sqrt(Vbeta), log=T) + dnorm(beta_post[ite-1], beta_prop, sqrt(Vprop), log=T)
    razao.l.den   = sum(dnorm(y, x*beta_post[ite-1], sqrt(V), log=T)) + dnorm(beta_post[ite-1], 0, sqrt(Vbeta), log=T) + dnorm(beta_prop, beta_post[ite-1], sqrt(Vprop), log=T)
    #print(dnorm(beta_post[ite-1], beta_prop, sqrt(Vprop), log=T))
    #print(dnorm(beta_prop, beta_post[ite-1], sqrt(Vprop), log=T))
    #razao.l.num   = sum(dnorm(y, x*beta_prop, sqrt(V), log=T)) + dnorm(beta_prop, 0, sqrt(Vbeta), log=T) 
    #razao.l.den   = sum(dnorm(y, x*beta_post[ite-1], sqrt(V), log=T)) + dnorm(beta_post[ite-1], 0, sqrt(Vbeta), log=T) 
    razao       = exp( razao.l.num-razao.l.den )
    alpha       = min(1,razao)
    u           = runif(1,0,1)
    if (u<alpha){
        beta_post[ite]  = beta_prop
        cont            = cont+1
    }else{
        beta_post[ite]  = beta_post[ite-1]
    }
}

#taxa de aceitação
cont/nite

## [1] 0.3647959

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

par(mfrow=c(1,1))
hist(beta_post[ind],main="Histograma", freq=F,
    xlab=expression(beta), ylab="densidade", bty="n", lwd=2)
abline(v=beta,lwd=2, col="red", lty=3)
abline(v=mean(beta_post[ind]),lwd=2, col="purple", lty=3)
abline(v=quantile(beta_post[ind],c(0.025, 0.975)),lwd=2, col="purple", lty=5)
#consideramos que a distribuição a posteriori de beta era desconhecida apenas
#para ilustrar o método mas ela é conhecida. Logo, plotaremos no histograma a curva
#da distribuição verdadeira
variancia = 1/(sum(x^2)/V + 1/Vbeta)
media     =  variancia * sum(y*x)/V
curve(dnorm(x, media, sqrt(variancia)), add=T, lwd=2, col="red")

#summary(glm(y~x-1))

2.2 Distribuição proposta II

#Supondo beta desconhecido, V conhecido e a 
#supondo que não sabemos gerar da distribuição a posteriori de beta 
#Amostrando beta pelo Metropolis-Hastings

#Definindo alguns parâmetros
naquecimento    = 100 #total de iterações consideradas como período de aquecimento
nespacamento    = 5 #total de iterações consideradas como espaçamento
nite            = (naquecimento + 1)+(1000-1)*nespacamento #para que a amostra tenha 
#nite            = 1500

beta_post   = NULL #inicializando vetor para guardar a amostra da distribuição a posteriori
Vbeta       = 1000 #hiperparâmetro

#Inicializando a cadeia
beta_post[1] = 0

#Inicializando o contador para saber sobre a taxa de aceitação
cont = 0

# Usaremos a seguinte distribuição proposta:
# beta^p ~ N(mprop, Vprop)
Vprop         = 0.1
mprop         = 0.5

#Executando o MCMC
for(ite in 2:nite)
{
    beta_prop     = rnorm(1, mprop, sqrt(Vprop))
    razao.l.num   = sum(dnorm(y, x*beta_prop, sqrt(V), log=T)) + dnorm(beta_prop, 0, sqrt(Vbeta), log=T) + dnorm(beta_post[ite-1], mprop, sqrt(Vprop), log=T)
    razao.l.den   = sum(dnorm(y, x*beta_post[ite-1], sqrt(V), log=T)) + dnorm(beta_post[ite-1], 0, sqrt(Vbeta), log=T) + dnorm(beta_prop, mprop, sqrt(Vprop), log=T)
    razao       = exp( razao.l.num-razao.l.den )
    alpha       = min(1,razao)
    u           = runif(1,0,1)
    if (u<alpha){
        beta_post[ite]  = beta_prop
        cont            = cont+1
    }else{
        beta_post[ite]  = beta_post[ite-1]
    }
}

#taxa de aceitação
cont/nite

## [1] 0.04415228

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

par(mfrow=c(1,1))
hist(beta_post[ind],main="Histograma", freq=F,
    xlab=expression(beta), ylab="densidade", bty="n", lwd=2)
abline(v=beta,lwd=2, col="red", lty=3)
abline(v=mean(beta_post[ind]),lwd=2, col="purple", lty=3)
abline(v=quantile(beta_post[ind],c(0.025, 0.975)),lwd=2, col="purple", lty=5)
#consideramos que a distribuição a posteriori de beta era desconhecida apenas
#para ilustrar o método mas ela é conhecida. Logo, plotaremos no histograma a curva
#da distribuição verdadeira
variancia = 1/(sum(x^2)/V + 1/Vbeta)
media     =  variancia * sum(y*x)/V
curve(dnorm(x, media, sqrt(variancia)), add=T, lwd=2, col="red")

#summary(glm(y~x-1))

3 Amostrador de Gibbs com passos de Metropolis-Hastings

Suponha

\(\beta\) e \(V\) desconhecidos e a priori que \[\begin{eqnarray} \beta &\sim& N(0, V_\beta), \\ V &\sim& GI(a, b). \end{eqnarray}\]
Considere que \[\begin{eqnarray*} p(\beta, V|\boldsymbol{y}) \propto p(\boldsymbol{y}|\beta, V) p(V) p(\beta) \end{eqnarray*}\] sendo \(p(\boldsymbol{y}|\beta, V)\) a f.d.p. da \(N(X_i \beta, V)\), \(p(\beta)\) a f.d.p. da \(N(0, V_\beta)\) e \(p(V)\) a f.d.p. da Gama Inversa com parâmetros \(a\) e \(b\).
Suponha que a distribuição condicional completa a posteriori (DCCP) de \(\beta\) é desconhecida mas que a DCCP de V é conhecida.
Usaremos então o algoritmo do Amostrador de Gibbs com passos de Metropolis-Hastings para obter uma amostra desta distribuição.

Algoritmo:

Escolha um valor arbitrário para \(\beta^{(1)}\) e \(V^{(1)}\);
Inicializa-se o contador de iteração \(t=2\);
Gera-se \(\beta^{(p)} \sim q(\beta^{(p)}|\beta^{(t-1)})\). Aceita-se o ponto gerado com probabilidade \(\alpha = \min\left\{1,\displaystyle\frac{p(\beta^{(p)})}{q(\beta^{(p)}|\beta^{(t-1)})}\frac{q(\beta^{(t-1)}|\beta^{(p)})}{p(\beta^{(t-1)})} \right\}\), onde \(p(\cdot)\) é a distribuição de interesse; Para decidir sobre a aceitação:
1. calcule \(\alpha\) e gere \(u \sim U(0,1)\);
  1. se \(u<\alpha\), aceite \(\beta^p\) fazendo \(\beta^{(t)}=\beta^p\); caso contrário, rejeite e faça \(\beta^{(t)}=\beta^{(t-1)}\);
Gera-se \(V\) da distribuição condicional completa a posteriori \(p(V|\boldsymbol{y}, \beta)\) que é uma \(GI\left( \frac{n}{2} + a, b+\frac{1}{2}\sum_{i=1}^n\left(Y_i-X_i\beta\right)^2\right)\)
Altera-se o contador de \(t\) para \(t+1\);
Repete-se os itens 3, 4 e 5 até que a convergência seja obtida.

3.1 Distribuição proposta: passeio aleatório

#Supondo beta desconhecido, V desconhecido e 
#supondo que não sabemos gerar da DCCP de beta 
#mas sabemos gerar da DCCP de V
#Amostrando beta pelo Metropolis-Hastings V diretamente da sua DCCP

#Definindo alguns parâmetros
naquecimento    = 100 #total de iterações consideradas como período de aquecimento
nespacamento    = 5 #total de iterações consideradas como espaçamento
nite            = (naquecimento + 1)+(1000-1)*nespacamento #para que a amostra tenha 
#nite            = 1500

#inicializando vetor para guardar a amostra da distribuição de interesse
beta_post   = NULL 
V_post      = NULL 
#hiperparâmetros
Vbeta       = 1000 
a           = 0.5
b           = 0.5

#Inicializando a cadeia
beta_post[1]  = 0
V_post[1]     = 5

#Inicializando o contador para saber sobre a taxa de aceitação
cont = 0

# Usaremos a seguinte distribuição proposta:
# beta^p ~ N(beta^{t-1}, Vprop)
Vprop         = 0.001

#Executando o MCMC
for(ite in 2:nite)
{
  #Amostrando beta
    beta_prop     = rnorm(1, beta_post[ite-1], sqrt(Vprop))
    razao.l.num   = sum(dnorm(y, x*beta_prop, sqrt(V_post[ite-1]), log=T)) + dnorm(beta_prop, 0, sqrt(Vbeta), log=T) + dnorm(beta_post[ite-1], beta_prop, sqrt(Vprop), log=T)
    razao.l.den   = sum(dnorm(y, x*beta_post[ite-1], sqrt(V_post[ite-1]), log=T)) + dnorm(beta_post[ite-1], 0, sqrt(Vbeta), log=T) + dnorm(beta_prop, beta_post[ite-1], sqrt(Vprop), log=T)
    #print(dnorm(beta_post[ite-1], beta_prop, sqrt(Vprop), log=T))
    #print(dnorm(beta_prop, beta_post[ite-1], sqrt(Vprop), log=T))
    #razao.l.num   = sum(dnorm(y, x*beta_prop, sqrt(V), log=T)) + dnorm(beta_prop, 0, sqrt(Vbeta), log=T) 
    #razao.l.den   = sum(dnorm(y, x*beta_post[ite-1], sqrt(V), log=T)) + dnorm(beta_post[ite-1], 0, sqrt(Vbeta), log=T) 
    razao       = exp( razao.l.num-razao.l.den )
    alpha       = min(1,razao)
    u           = runif(1,0,1)
    if (u<alpha){
        beta_post[ite]  = beta_prop
        cont            = cont+1
    }else{
        beta_post[ite]  = beta_post[ite-1]
    }

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

#taxa de aceitação
cont/nite

## [1] 0.3720565

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

#com toda a amostra
acf(beta_post,main="Autocorrelação da cadeia", 
    ylab=expression(beta), xlab="defasagem")
#removendo o período de aquecimento
acf(beta_post[ind],main="Autocorrelação da cadeia", 
    ylab=expression(beta), xlab="defasagem")
#com toda a amostra
acf(V_post,main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")
#removendo o período de aquecimento
acf(V_post[ind],main="Autocorrelação da cadeia", 
    ylab=expression(V), xlab="defasagem")

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

4 Exercício 1:

Seja \[\begin{eqnarray} \theta &\sim& GI\left( a, b\right). \end{eqnarray}\]

Suponha que não sabemos amostrar desta distribuição. Obtenha uma amostra desta distribuição usando o Metropolis-Hastings.

5 Exercício 2:

Seja \[\begin{eqnarray} Y_i &\sim& N\left( \beta_0 + X_i \beta_1, V\right), \; i=1, \ldots, n.\\ \end{eqnarray}\]
A priori, considere que \[\begin{eqnarray} \beta_0 &\sim& N(0, V_{\beta_0}), \\ \beta_1 &\sim& N(0, V_{\beta_1}),\\ V &\sim& GI(a, b). \end{eqnarray}\]
Implemente o algoritmo amostrador de Gibbs usando o Metropolis-Hastings para obter uma amostra de V e amostrando diretamente da DCCP de \(\beta_0\) e \(\beta_1\).

Inferência Bayesiana - GET/IME/UFF

Patrícia Lusié Velozo da Costa

2023-1