1 .Exercícios de simulação

Instruções:

Simule uma amostra de \(\scriptsize Y = (Y_1,...,Y_n)\) do modelo descrito em cada exercício considerando \(\scriptsize n = 1000\). Depois de gerada a amostra, suponha que os parâmetros deste modelo são desconhecidos e utilize a inferência bayesiana para estimá-los.

Considere o seguinte modelo para explicar uma determinada característica (\(\scriptsize Y\)) da população:

\[\scriptsize Y\sim N(\beta_0+\beta_1X,\phi^-1)\]

Considere \(\scriptsize X\) uma covariável do seu modelo e gere esta de uma distribuição normal. Os parâmetros deste modelo são \(\scriptsize \theta=(\beta_0,\beta_1,\phi^-1)\). Assuma a seguinte crença que supõe independência a priori:

\[\scriptsize \beta_i\sim N(a_i, b_i),i= 0,1\] e \[\scriptsize \phi\sim Gama(c,d)\].

Para obtenção da amostra da distribuição a posteriori por meio de técnica de MCMC, utilize dois algoritmos distintos: o amostador de Gibbs e o amostador de Gibbs com um passo de Metropolis-Hastings(para o parâmetro \(\scriptsize \phi\)).

\(\scriptsize Y\sim N(\beta_0 +\beta_1X,\phi{^-1})\) onde \(\scriptsize \theta=(\beta_0 + \beta_1 + \phi^{-1})\)

Supondo independência, temos :

\(\scriptsize p(\theta) = p(\beta_0)* p(\beta_1)* p(\phi^{-1})\)

onde:

\(\scriptsize \beta_0\sim N(a_0, b_0)\)

\(\scriptsize \beta_1\sim N(a_1, b_1)\)

\(\scriptsize \phi^-1 \sim Gama(c,d)\)

library(coda)

1.1 Simulação dos dados

#Tamanho da amostra
n<-1000  

#Verdadeiro valor dos parâmetros
beta0.real<-7 
beta1.real<-1
phi.real<-2

#Simulação dos valores da covariável X
x<- rnorm(n,0,1)

#Calculo da média de Y
mu.real<- beta0.real+(beta1.real*x)

#Simulação os valores de Y
y<-rnorm(n,mu.real,sqrt(1/phi.real))

#Histograma de Y
hist(y,col="darkred",main="",xlab=expression(y),ylab="Frequencia")

1.2 Atribuição dos hiperparâmetros

\(\scriptsize \beta_0\sim N(6,8)\)

\(\scriptsize \beta_1\sim N(2,5)\)

\(\scriptsize \phi^-1 \sim Gama(6,3)\)

# Hiperparametros das distribuições a priori

#Beta_0
a0<-6
b0<-8

#Beta_1
a1<- 2
b1<- 5

#Phi
c<-6
d<-3

par(mfrow=c(2,3))

# Priori para Beta_0
plot(function(x)dnorm(x,a0,sqrt(b0)),-30,30,xlab="",ylab=expression(beta_0),col="darkred")

# Priori para Beta_1
plot(function(x)dnorm(x,a1,sqrt(b1)),-30,30,xlab="",ylab=expression(beta_1),col="darkred")

# Priori para phi
plot(function(x)dgamma(x,c,d),0,10,xlab="",ylab=expression(phi),col="darkred")

1.3 Contas

1.3.1 Distribuição a posteriori

\[\scriptsize p(\theta|Y) \propto p(Y|\theta) p(\theta)\]

\[\tiny p(\theta|Y) \propto (2\pi\phi^{-1})^{-n/2} exp\Bigg\{ -\frac{1}{2\phi^{-1}}\sum^n_{i=1}\Bigg [Y_i - (\beta_0 + \beta_1X_i))^2\Bigg]\Bigg\}(2\pi b_0)^{-1/2}exp\Bigg\{-\frac{1}{2b_0} (\beta_0 -a_0)^2\Bigg\}(2\pi b_1)^{-1/2}exp\Bigg\{-\frac{1}{2b_1} (\beta_1 -a_1)^2\Bigg\}\phi^{c-1}exp\Bigg\{-d\phi\Bigg\}\]

Núcleo da distribuicao a posteriori:

\[\scriptsize p(\theta|Y) \propto (\phi^{-1})^{-n/2} exp\Bigg\{-\frac{1}{2\phi^{-1}} \sum^n_{i=1}\Bigg [Y_i - (\beta_0 + \beta_1X_i))^2\Bigg]\Bigg\} exp\Bigg\{-\frac{1}{2b_0} (\beta_0 -a_0)^2\Bigg\}exp\Bigg\{-\frac{1}{2b_1} (\beta_1 -a_1)^2\Bigg\}\phi^{c-1} exp\Bigg\{-d\phi\Bigg\}\]

1.3.2 Distribuição condicional completa para \(\beta_0\):

\[\scriptsize p(\beta_0|.) \propto exp\Bigg\{-\frac{1}{2\phi^{-1}} \sum^n_{i=1}\Bigg [Y_i - (\beta_0 + \beta_1X_i))^2\Bigg]\Bigg\} exp\Bigg\{-\frac{1}{2b_0} (\beta_0 -a_0)^2\Bigg\}\] \[\scriptsize p(\beta_0|.) \propto exp\Bigg\{-\frac{1}{2\phi^{-1}} \sum^n_{i=1}\Bigg [Y_i^2 - 2Y_i (\beta_0 + \beta_1X_i))+ (\beta_0 + \beta_1X_i)^2\Bigg]\Bigg\} exp\Bigg\{-\frac{1}{2b_0} (\beta_0^2 -2\beta_0a_0 + a_0^2)\Bigg\}\] \[\scriptsize p(\beta_0|.) \propto exp\Bigg\{-\frac{1}{2\phi^{-1}}\Bigg[\sum^n_{i=1}Y_i^2 - 2\beta_0\sum^n_{i=1}Y_i - 2\beta_1\sum^n_{i=1}Y_iX_i+ n\beta_0^2 + 2\beta_0\beta_1\sum^n_{i=1}X_i+ \beta_1^2\sum^n_{i=1}X_i^2)\Bigg]\Bigg\} exp\Bigg\{-\frac{1}{2b_0} (\beta_0^2 -2\beta_0a_0)\Bigg\}\] \[\scriptsize p(\beta_0|.)\propto exp\Bigg\{-\frac{1}{2}\Bigg[-2\beta_0\Bigg(\frac{\sum^n_{i=1}Y_i-\beta_1\sum^n_{i=1}X_i}{\phi^{-1}}+\frac{a_0}{b_0}\Bigg)+ \beta_0^2\Bigg(\frac{n}{\phi^{-1}}+\frac{1}{b_0}\Bigg)\Bigg]\Bigg\}\]

\[\scriptsize p(\beta_0|.)\propto exp\Bigg\{-\frac{1}{2}\Bigg[\Bigg(\frac{b_0n +\phi{^-1}}{\phi^{-1}b_0}\Bigg)\Bigg[\beta_o^2-2\beta_0\Bigg(\frac{b_0(\sum^n_{i=1}Y_i-\beta_1\sum^n_{i=1}Xi)+a_0\phi^{-1}}{\phi^{-1}b_0} \Bigg)\Bigg(\frac{b_0n +\phi{^-1}}{\phi^{-1}b_0}\Bigg)^{-1} \Bigg] \Bigg]\Bigg\}\] \[\scriptsize p(\beta_0|.)\propto exp\Bigg\{-\frac{1}{2}\Bigg[\Bigg(\frac{b_0n +\phi{^-1}}{\phi^{-1}b_0}\Bigg)\Bigg[\beta_o^2-2\beta_0\Bigg(\frac{b_0(\sum^n_{i=1}Y_i-\beta_1\sum^n_{i=1}Xi)+a_0\phi^{-1}}{b_0n +\phi{^-1}}\Bigg)\Bigg] \Bigg]\Bigg\}\]

\[\scriptsize \implies\beta_0|.\sim N \Bigg(\Bigg(\frac{b_0(\sum^n_{i=1}Y_i+\beta_1\sum^n_{i=1}Xi)+a_0\phi^{-1}}{\phi^{-1}+ nb_0} \Bigg),\Bigg( \frac{\phi^{-1}b_0}{b_0n+\phi^{-1}}\Bigg)\Bigg)\]

1.3.3 Distribuição condicional completa para \(\beta_1\):

\[\scriptsize p(\beta_1|.) \propto exp\Bigg\{-\frac{1}{2\phi^{-1}} \sum^n_{i=1}\Bigg [Y_i - (\beta_0 + \beta_1X_i))^2\Bigg]\Bigg\} exp\Bigg\{-\frac{1}{2b_1} (\beta_1 -a_1)^2\Bigg\}\] \[\scriptsize p(\beta_1|.) \propto exp\Bigg\{-\frac{1}{2\phi^{-1}} \sum^n_{i=1}\Bigg [Y_i^2 - 2Y_i (\beta_0 + \beta_1X_i))+ (\beta_0 + \beta_1X_i)^2\Bigg]\Bigg\} exp\Bigg\{-\frac{1}{2b_1} (\beta_1^2 -2\beta_1a_1 + a_1^2)\Bigg\}\] \[\scriptsize p(\beta_1|.) \propto exp\Bigg\{-\frac{1}{2\phi^{-1}}\Bigg[\sum^n_{i=1}Y_i^2 - 2\beta_0\sum^n_{i=1}Y_i - 2\beta_1\sum^n_{i=1}Y_iX_i + \beta_0^2 + 2\beta_0\beta_1\sum^n_{i=1}X_i+ \beta_1^2\sum^n_{i=1}X_i^2)\Bigg]\Bigg\} exp\Bigg\{-\frac{1}{2b_1} (\beta_1^2 -2\beta_1a_1)\Bigg\}\] \[\scriptsize p(\beta_1|.)\propto exp\Bigg\{-\frac{1}{2}\Bigg[-2\beta_1\Bigg(\frac{\sum^n_{i=1}Y_iX_i-\beta_0\sum^n_{i=1}X_i}{\phi^{-1}}+\frac{a_1}{b_1}\Bigg)+ \beta_1^2\Bigg(\frac{\sum^n_{i=1}X_i^2}{\phi^{-1}}+\frac{1}{b_1}\Bigg)\Bigg]\Bigg\}\] \[\scriptsize p(\beta_1|.)\propto exp\Bigg\{-\frac{1}{2}\Bigg[-2\beta_1\Bigg(\frac{b_1(\sum^n_{i=1}Y_iX_i-\beta_0\sum^n_{i=1}X_i)+ a_1\phi^{-1}}{\phi^{-1}b_1}\Bigg)+ \beta_1^2\Bigg(\frac{b_1\sum^n_{i=1}X_i^2+ \phi^{-1}}{\phi^{-1}b_1}\Bigg)\Bigg]\Bigg\}\]

\[\scriptsize p(\beta_1|.)\propto exp\Bigg\{-\frac{1}{2}\Bigg[\Bigg(\frac{b_1\sum^n_{i=1}X_i^2 +\phi{^-1}}{\phi^{-1}b_1}\Bigg)\Bigg[\beta_1^2-2\beta_1\Bigg(\frac{b_1(\sum^n_{i=1}Y_iXi-\beta_0\sum^n_{i=1}X_i)+a_1\phi^{-1}}{\phi^{-1}b_1} \Bigg)\Bigg(\frac{b_1\sum^n_{i=1}X_i^2 +\phi{^-1}}{\phi^{-1}b_1}\Bigg)^{-1} \Bigg] \Bigg]\Bigg\}\] \[\scriptsize p(\beta_1|.)\propto exp\Bigg\{-\frac{1}{2}\Bigg[\Bigg(\frac{b_1\sum^n_{i=1}X_i^2 +\phi{^-1}}{\phi^{-1}b_1}\Bigg)\Bigg[\beta_1^2-2\beta_1\Bigg(\frac{b_1(\sum^n_{i=1}Y_iXi-\beta_0\sum^n_{i=1}X_i)+a_1\phi^{-1}}{b_1\sum^n_{i=1}X_i^2 +\phi{^-1}}\Bigg)\Bigg] \Bigg]\Bigg\}\]

\[\scriptsize \implies\beta_1|.\sim N \Bigg(\Bigg(\frac{b_1(\sum^n_{i=1}Y_iXi-\beta_0\sum^n_{i=1}X_i)+a_1\phi^{-1}}{b_1\sum^n_{i=1}X_i^2 +\phi{^-1}} \Bigg),\Bigg( \frac{\phi^{-1}b_1}{b_1\sum^n_{i=1}X_i^2 +\phi{^-1}}\Bigg)\Bigg)\]

1.3.4 Distribuição condicional completa para \(\phi\):

\[\scriptsize p(\phi|.)\propto(\phi^{-1})^{-n/2}exp\Bigg\{-\frac{1}{2\phi^{-1}} \sum^n_{i=1}\Bigg [(Y_i - (\beta_0 + \beta_1X_i))^2\Bigg]\Bigg\}\phi^{c-1} exp\Bigg\{-d\phi\Bigg\}\] \[\scriptsize p(\phi|.)\propto\phi^{(\frac{n}{2}+c-1)}exp\Bigg\{-\phi\Bigg[\frac{\sum^n_{i=1}(Y_i-(\beta_0 +\beta_1X_i))^2}{2}+d\Bigg]\Bigg\}\] \[\scriptsize p(\phi|.)\propto\phi^{(\frac{n}{2}+c-1)}exp\Bigg\{-\phi\Bigg[\frac{\sum^n_{i=1}(Y_i^2- 2Y_i(\beta_0 +\beta_1X_i)+ (\beta_0+\beta_1X_i)^2)}{2}+d\Bigg]\Bigg\}\]

\[\scriptsize p(\phi|.)\propto\phi^{(\frac{n}{2}+c-1)}exp\Bigg\{-\phi\Bigg[\frac{\sum^n_{i=1}Y_i^2-2\beta_0\sum^n_{i=1}Y_i-2\beta_1\sum^n_{i=1}Y_iX_i+n\beta_0^2+2\beta_0\beta_1\sum^n_{i=1}X_i+ \beta_1^2\sum^n_{i=1}X_i^2}{2}+d\Bigg]\Bigg\}\]

\[\scriptsize \implies\phi|.\sim GAMA\Bigg(\Bigg(\frac{n}{2}+c\Bigg),\Bigg(\frac{\sum^n_{i=1}(Y_i-(\beta_0 +\beta_1X_i))^2}{2}+d\Bigg)\Bigg)\] ou

\[\scriptsize \implies\phi|.\sim GAMA\Bigg(\Bigg(\frac{n}{2}+c\Bigg),\Bigg(\frac{l}{2}+d\Bigg)\Bigg)\] onde:

\[\scriptsize l =(\sum^n_{i=1}Y_i^2-2\beta_0\sum^n_{i=1}Y_i-2\beta_1\sum^n_{i=1}Y_iX_i+n\beta_0^2+2\beta_0\beta_1\sum^n_{i=1}X_i+ \beta_1^2\sum^n_{i=1}X_i^2) \]

2 Algoritmo para o amostrador de Gibbs

Defina valores inicias para os parâmetros \(\beta_0,\beta_1,\phi\)
Obter um novo valor \(\theta^{(i)} = (\beta_0^{(i)}, \beta_1^{(i)},\phi^{(i)})\) ,através da geração sucessivas das distribuições:
Distribuição condicional completa de \(\beta_0|.\)
Distribuição condicional completa de \(\beta_1|.\), utilizando os valores gerados em 3.
Distribuição condicional completa de \(\phi|.\), utilizando os valores gerados em 3 e 4.
Repita o processo a partir do passo 2 até obter a convergência.

# Valores iniciais
beta0<-NULL
beta1<-NULL
phi<-NULL

beta0[1]<- 5
beta1[1]<- 2
phi[1]<- 4

#Número de interações 
ite<-1000 
n<-length(y)


for(i in 2:ite){
  
  ### Amostrando Beta_0
  
  media<- ((b0*sum(y)) - (b0*beta1[i-1]*sum(x))+ (a0/phi[i-1])) / ((n*b0) +(1/phi[i-1]))
  
  var<- (b0/phi[i-1])/((n*b0) + (1/phi[i-1]))
  
  beta0[i]<-rnorm(1,media,sqrt(var))

  ### Amostrando Beta_1
  
  media1<- ((b1*sum(x*y))-(b1*beta0[i]*sum(x)) + (a1/phi[i-1])) / ((sum(x^2)*b1) + (1/phi[i-1]))
  
  var1<- (b1/phi[i-1]) / (sum(x^2)*b1 + (1/phi[i-1]))
 
  beta1[i]<-rnorm(1,media1,sqrt(var1))

  # Amostrando phi
   
  alpha<-(n/2)+ c
 
  beta_phi<- (sum(y^2)/2) - (beta0[i]*sum(y)) - (beta1[i]*sum(y*x)) + ((n*beta0[i]^2)/2) + (beta0[i]*beta1[i]*sum(x)) + (((beta1[i]^2)*sum(x^2))/2)
  
  beta_phi1<- beta_phi + d
  
  phi[i]<-rgamma(1,alpha,beta_phi1)
}

2.0.1 Análise visual das cadeias antes do tratamento

#Beta0
plot(beta0, type="l", xlab="")
abline(h=beta0.real, lwd=3, col=2)

#Beta1
plot(beta1,type="l",xlab="")
abline(h=beta1.real,lwd=3,col=2)

#Phi
plot(phi,type="l",xlab="")
abline(h=phi.real,lwd=3,col=2)

2.1 Tratamento da cadeia

2.1.1 Amostra para aquecimento

burnin<-100
beta0_aq<- beta0[(burnin+1):1000]
beta1_aq<- beta1[(burnin+1):1000]
phi_aq<- phi[(burnin+1):1000]

2.1.2 Analisando a autocorrelação

Podemos observar nos gráficos da autocorrelação, que não é necessario retirar espaçamento nas cadeias de \(\scriptsize\beta_0, \beta_1,\phi\). Também foi gerado uma amostra de tamanho 4000 para realizar o teste de diagnóstico de Raftery e Lewis (1992), que também confirmou o resultado.

##Autocorrelação
acf(beta0_aq , main = "Autocorrelação Beta 0")

acf(beta1_aq , main = "Autocorrelação Beta 1")

acf(phi_aq , main = "Autocorrelação Phi")

2.1.3 Método de Gelman e Rubin

Portanto, para realizar o teste foi gerada uma segunda amostra do algoritmo de Gibbs, com diferentes valores iniciais.

# Valores iniciais
beta0_2<-NULL
beta1_2<-NULL
phi_2<-NULL

beta0_2[1]<- 10
beta1_2[1]<- 3
phi_2[1]<- 1

#Número de interações 
ite<-1000 
n<-length(y)

# MCMC 
for(i in 2:ite){
  
  ### Amostrando Beta_0
  
  media2<- ((b0*sum(y)) - (b0*beta1_2[i-1]*sum(x))+ (a0/phi_2[i-1])) / ((n*b0) +(1/phi_2[i-1]))
  
  var2<- (b0/phi_2[i-1])/(n*b0 + (1/phi_2[i-1]))
  
  beta0_2[i]<-rnorm(1,media2,sqrt(var2))
  
  
  ### Amostrando Beta_1
  
  media12<- ((b1*sum(x*y))-(b1*beta0_2[i]*sum(x)) + (a1/phi_2[i-1])) / ((sum(x^2)*b1) + (1/phi_2[i-1]))
  
  var12<- (b1/phi_2[i-1]) / (sum(x^2)*b1 + (1/phi_2[i-1]))
 
  beta1_2[i]<-rnorm(1,media12,sqrt(var12))
  
  
  # Amostrando phi_2
   
  alpha2<-(n/2)+ c
 
  beta_phi_2<- (sum(y^2)/2) - (beta0_2[i]*sum(y)) - (beta1_2[i]*sum(y*x)) + ((n*beta0_2[i]^2)/2) + (beta0_2[i]*beta1_2[i]*sum(x)) + (((beta1_2[i]^2)*sum(x^2))/2)
  
  beta_phi_12<- beta_phi_2 + d
  
  #beta_phi12<- ((sum(( y - (beta0_2[i] + beta1_2[i] * x))^2)) /2 )+ d
  
  phi_2[i]<-rgamma(1,alpha2,beta_phi_12)
  
}

Análise visual das cadeias antes do tratamento

#Beta0_2
plot(beta0_2, type="l", xlab="")
abline(h=beta0.real, lwd=3, col=2)

#Beta1_2
plot(beta1_2, type="l", xlab="")
abline(h=beta1.real, lwd=3, col=2)

#phi_2
plot(phi_2, type="l", xlab="")
abline(h=phi.real, lwd=3, col=2)

Amostra para aquecimento

burnin<-100
beta0_aq2<- beta0_2[(burnin+1):1000]
beta1_aq2<- beta1_2[(burnin+1):1000]
phi_aq2<- phi_2[(burnin+1):1000]

Analisando a autocorrelação

Como pode ser visto nos gráficos da autocorrelação, não é necessario retirar espaçamento nas cadeias de \(\beta_0, \beta_1,\phi\).

## Autocorrelação
acf(beta0_aq2 , main = "Autocorrelação Beta 0")

acf(beta1_aq2 , main = "Autocorrelação Beta 1")

acf(phi_aq2 , main = "Autocorrelação Phi")

O resultado do teste de Gelman e Rubin (1992) mostra que as cadeias convergem. Como podemos observar todos os valores do “Point est.” estão próximos de 1.

cadeia.beta0<-mcmc.list(as.mcmc(beta0_aq),as.mcmc(beta0_aq2))
cadeia.beta1<-mcmc.list(as.mcmc(beta1_aq),as.mcmc(beta1_aq2))
cadeia.phi<-mcmc.list(as.mcmc(phi_aq),as.mcmc(phi_aq2))

gelman.diag(cadeia.beta0)

## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]      0.999          1

gelman.diag(cadeia.beta1)

## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]       1.01       1.03

gelman.diag(cadeia.phi)

## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]          1          1

2.2 Traços das cadeia

#Beta0
plot(beta0_aq, type="l", xlab=" ", main="Traço da cadeia de Beta0")
abline(h=beta0.real, lwd=3, col=2)

#Beta1
plot(beta1_aq,type="l",xlab=" ", main="Traço da cadeia de Beta1")
abline(h=beta1.real,lwd=3,col=2)

#Phi
plot(phi_aq,type="l",xlab=" ", main="Traço da cadeia de phi")
abline(h=phi.real,lwd=3,col=2)

2.3 Estimativas a posteriori

resumo<- cbind(round(mean(beta0_aq),3),beta0.real, round(mean(beta1_aq),3),beta1.real,round(mean(phi_aq),3),phi.real)
colnames(resumo)<- c("Beta0 Estimado","Beta0 Real","Beta1 Estimado","Beta1 Real","Phi Estimado","Phi Real")

DT::datatable(resumo, rownames = FALSE)

#Beta0
intervalo_beta0 <- cbind(
round(HPDinterval(as.mcmc(beta0_aq)),3),beta0.real
)
colnames(intervalo_beta0)<- c("Intervalo HPD Beta0 95% (min)","Intervalo HPD Beta0 95% (máx)","Beta0 Real")

DT::datatable(intervalo_beta0, rownames = FALSE)

#Beta1
intervalo_beta1 <- cbind(
round(HPDinterval(as.mcmc(beta1_aq)),3),beta1.real
)
colnames(intervalo_beta1)<- c("Intervalo HPD Beta1 95% (min)","Intervalo HPD Beta1 95% (máx)","Beta1 Real")

DT::datatable(intervalo_beta1, rownames = FALSE)

#Phi
intervalo_phi <- cbind(
round(HPDinterval(as.mcmc(phi_aq)),3),phi.real
)
colnames(intervalo_phi)<- c("Intervalo HPD Phi 95% (min)","Intervalo HPD Phi 95% (máx)","Phi Real")

DT::datatable(intervalo_phi, rownames = FALSE)

3 Algoritmo para o amostrador de Gibbs com um passo de Metropolis-Hatings

3.0.1 Distribuição Proposta

\(\scriptsize c= log(\phi)\sim N(log(\phi^{(i)},var)\)

A razão de metropolis é dada por:

\(\scriptsize RM= \frac{\pi{\phi^{p}}}{\pi{\phi^{(i)}}} \frac{\phi^p}{\phi^{(i)}}\)

\(\scriptsize log\big(\frac{\pi{\phi^{p}}}{\pi{\phi^{(i)}}} \frac{\phi^p}{\phi^{(i)}}\big)= log(\pi{\phi^{p}})- log(\pi{\phi^{(i)}})+ log(\phi^p) - log(\phi^{(i)})\)

onde:

\[\scriptsize \pi{\phi}= \phi^{(\frac{n}{2}+c-1)}exp\Bigg\{-\phi\Bigg[\frac{\sum^n_{i=1}Y_i^2-2\beta_0\sum^n_{i=1}Y_i-2\beta_1\sum^n_{i=1}Y_iX_i+n\beta_0^2+2\beta_0\beta_1\sum^n_{i=1}X_i+ \beta_1^2\sum^n_{i=1}X_i^2}{2}+d\Bigg]\Bigg\}\] sendo

\[\scriptsize l =(\sum^n_{i=1}Y_i^2-2\beta_0\sum^n_{i=1}Y_i-2\beta_1\sum^n_{i=1}Y_iX_i+n\beta_0^2+2\beta_0\beta_1\sum^n_{i=1}X_i+ \beta_1^2\sum^n_{i=1}X_i^2) \] \[\scriptsize \pi{\phi}= \phi^{(\frac{n}{2}+c-1)}exp\Bigg\{-\phi\Bigg[\frac{l}{2}+d\Bigg]\Bigg\}\]

3.1 Algoritmo para o amostrador de Gibbs com passo de Metropolis_Hastings

Defina valores inicias para os parâmetros \(\scriptsize \beta_0,\beta_1,\phi\)
Obter um novo valor \(\scriptsize \theta^{(i)} = (\beta_0^{(i)}, \beta_1^{(i)},\phi^{(i)})\) ,através da geração sucessivas das distribuições:
Distribuição condicional completa de \(\scriptsize \beta_0|.\)
Distribuição condicional completa de \(\scriptsize \beta_1|.\), utilizando os valores gerados em 3.
O auxiliar c com uma distribuição \(\scriptsize c\sim N(log(\phi^{(i)},var)\) onde o phi proposto é \(\scriptsize \phi^p= exp(c)\), utilizando os valores gerados em 3 e 4.
Calcular \(\scriptsize \alpha=min(1,\frac{\pi{\phi^{p}}}{\pi{\phi^{(i)}}} \frac{\phi^p}{\phi^{(i)}})\)
Gerar um valor u da distribuição uniforme(0,1)
Se \(\scriptsize u<=\alpha\) então o valor proposto é aceito como novo valor de phi proposto, \(\scriptsize \phi^{(i+1)}= \phi^p\). Caso contrário, \(\scriptsize \phi^{(i+1)}= \phi^{(i)}\)
Repita o processo a partir do passo 2 até obter a convergência.

O algoritmo de Gibbs com passo de Metropolis Hastings possui taxa de aceitação de aproximadamente 40%, considerada boa(para passeio aleatório).

# Valores iniciais
beta0<-NULL
beta1<-NULL
phi<-NULL

beta0[1]<- 5
beta1[1]<- 2
phi[1]<- 4

#Número de interações 
ite<-1000 
n<-length(y)

#Contador 
contador <- 0

#Variância 
variancia <- 0.01

#Alpha
alpha <- 0


for(i in 2:ite){
  
  ### Amostrando Beta_0
  
  media<- ((b0*sum(y)) - (b0*beta1[i-1]*sum(x))+ (a0/phi[i-1])) / ((n*b0) +(1/phi[i-1]))
  
  var<- (b0/phi[i-1])/((n*b0) + (1/phi[i-1]))
  
  beta0[i]<-rnorm(1,media,sqrt(var))
 
  
  ### Amostrando Beta_1
  
  media1<- ((b1*sum(x*y))-(b1*beta0[i]*sum(x)) + (a1/phi[i-1])) / ((sum(x^2)*b1) + (1/phi[i-1]))
  
  var1<- (b1/phi[i-1]) / (sum(x^2)*b1 + (1/phi[i-1]))
 
  beta1[i]<-rnorm(1,media1,sqrt(var1))
  
 
   ## Phi com passo Metropolis
  
  phi_proposto1<- rnorm(1,log(phi[i-1]),sqrt(variancia))
  phi_proposto<- exp(phi_proposto1)
  
  
  l<- (((sum(y^2)/2) - (beta0[i]*sum(y)) - (beta1[i]*sum(y*x)) + ((n*beta0[i]^2)/2) + (beta0[i]*beta1[i]*sum(x)) + (((beta1[i]^2)*sum(x^2))/2))) + d
  
  
  log_numerador<-   ((n/2)+c-1)*log(phi_proposto) - phi_proposto*l + log(phi_proposto)
  log_denominador<- ((n/2)+c-1)*log(phi[i-1]) - phi[i-1]*l + log(phi[i-1])
  log_alpha<- log_numerador- log_denominador
  
  uniforme<- runif(1,0,1)
  
  if(log(uniforme) <= min(1,log_alpha)){
    phi[i] = phi_proposto
    contador = contador + 1
  }else{
    phi[i] = phi[i-1]
  }

}
aceitacao = contador/ite
aceitacao

## [1] 0.488

3.1.1 Análise visual das cadeias antes do tratamento

#Beta0
plot(beta0, type="l", xlab="")
abline(h=beta0.real, lwd=3, col=2)

#Beta1
plot(beta1,type="l",xlab="")
abline(h=beta1.real,lwd=3,col=2)

#Phi
plot(phi,type="l",xlab="")
abline(h=phi.real,lwd=3,col=2)

3.2 Tratamento da cadeia

3.2.1 Amostra para aquecimento

burnin<-100
beta0_aq<- beta0[(burnin+1):1000]
beta1_aq<- beta1[(burnin+1):1000]
phi_aq<- phi[(burnin+1):1000]

3.2.2 Analisando a autocorrelação

Podemos observar nos gráficos da autocorrelação, que na cadeia relacionada ao \(\phi\) temos lag 6. Logo, foi realizado um espaçamento de 7 em 7.

##Autocorrelação
acf(beta0_aq , main = "Autocorrelação Beta 0")

acf(beta1_aq , main = "Autocorrelação Beta 1")

acf(phi_aq , main = "Autocorrelação Phi")

## Tirando o espaçamento
thin<-7
beta0_aq_esp<-beta0_aq[seq(1,length(beta0_aq),by=thin)]
beta1_aq_esp<-beta1_aq[seq(1,length(beta1_aq),by=thin)]
phi_aq_esp<-phi_aq[seq(1,length(phi_aq),by=thin)]

Gráfico da autocorrelação depois do espaçamento.

## Analisando autocorrelação
acf(beta0_aq_esp , main = "Autocorrelação Beta 0")

acf(beta1_aq_esp , main = "Autocorrelação Beta 1")

acf(phi_aq_esp , main = "Autocorrelação Phi")

3.2.3 Método de Gelman e Rubin

Para verificar a convergência da distribuição a posteriori, foi realizado o método de Gelman e Rubin(1992), onde podemos observar que em muitos casos a convergência pode ser facilmente determinada provinda de sequências múltiplas independentes em paralelo. Portanto, para realizar o teste foi gerada uma segunda amostra do algoritmo de Gibbs com passo de Metropolis, com diferentes valores iniciais.

# Valores iniciais
beta0_2<-NULL
beta1_2<-NULL
phi_2<-NULL

beta0_2[1]<- 9
beta1_2[1]<- 3
phi_2[1]<- 1

#Número de interações 
ite<-1000 
n<-length(y)

#Contador 
contador <- 0

#Variância 
variancia <- 0.01

#Alpha
alpha2 <- 0


for(i in 2:ite){
  
  ### Amostrando Beta_0
  
  media2<- ((b0*sum(y)) - (b0*beta1_2[i-1]*sum(x))+ (a0/phi_2[i-1])) / ((n*b0) +(1/phi_2[i-1]))
  
  var2<- (b0/phi_2[i-1])/((n*b0) + (1/phi_2[i-1]))
  
  beta0_2[i]<-rnorm(1,media2,sqrt(var2))
 
  
  ### Amostrando Beta_1
  
  media12<- ((b1*sum(x*y))-(b1*beta0_2[i]*sum(x)) + (a1/phi_2[i-1])) / ((sum(x^2)*b1) + (1/phi_2[i-1]))
  
  var12<- (b1/phi_2[i-1]) / (sum(x^2)*b1 + (1/phi_2[i-1]))
 
  beta1_2[i]<-rnorm(1,media12,sqrt(var12))
 
 
   ## phi com passo Metropolis
  
  phi_proposto1<- rnorm(1,log(phi_2[i-1]),sqrt(variancia))
  phi_proposto<- exp(phi_proposto1)
  
  
  l<- (((sum(y^2)/2) - (beta0_2[i]*sum(y)) - (beta1_2[i]*sum(y*x)) + ((n*beta0_2[i]^2)/2) + (beta0_2[i]*beta1_2[i]*sum(x)) + (((beta1_2[i]^2)*sum(x^2))/2))) +d
  
  log_numerador<-   ((n/2)+c-1)*log(phi_proposto) - phi_proposto*l + log(phi_proposto)
  log_denominador<- ((n/2)+c-1)*log(phi_2[i-1]) - phi_2[i-1]*l + log(phi_2[i-1])
  log_alpha<- log_numerador- log_denominador
  
  uniforme<- runif(1,0,1)
  
  if(log(uniforme) <= min(1,log_alpha)){
    phi_2[i] = phi_proposto
    contador = contador + 1
  }else{
    phi_2[i] = phi_2[i-1]
  }
}

aceitacao = contador/ite
aceitacao

## [1] 0.456

Análise visual das cadeias antes do tratamento

#Beta0_2
plot(beta0_2, type="l", xlab="")
abline(h=beta0.real, lwd=3, col=2)

#Beta1_2
plot(beta1_2, type="l", xlab="")
abline(h=beta1.real, lwd=3, col=2)

#phi_2
plot(phi_2, type="l", xlab="")
abline(h=phi.real, lwd=3, col=2)

Amostra para aquecimento e espaçamento

burnin<-100
beta0_aq2<- beta0_2[(burnin+1):1000]
beta1_aq2<- beta1_2[(burnin+1):1000]
phi_aq2<- phi_2[(burnin+1):1000]

## Tirando o espaçamento
thin<-7
beta0_aq_esp2<-beta0_aq2[seq(1,length(beta0_aq2),by=thin)]
beta1_aq_esp2<-beta1_aq2[seq(1,length(beta1_aq2),by=thin)]
phi_aq_esp2<-phi_aq2[seq(1,length(phi_aq2),by=thin)]

Analisando a autocorrelação

## Autocorrelação
acf(beta0_aq2 , main = "Autocorrelação Beta 0")

acf(beta1_aq2 , main = "Autocorrelação Beta 1")

acf(phi_aq2 , main = "Autocorrelação Phi")