Amostrador de Gibbs

Para um MCMC, a amostra é obtida através de uma estrutura denominada amostrador. Sempre que possível, será utilizado o amostrador de Gibbs, que requer o conhecimento das distribuições condicionais completas. Neste contexto uma distribuição condicional é dita completa para \(\theta_i\) se a densidade \(p(\theta_i|\theta_{-i})\) é conhecida. Seja um vetor de parâmetros \(\theta\in\mathbb{R}^p\) associado a um modelo \(p(x|\theta)\) tal qual as distribuições condicionais completas de \(\theta_i\in\theta\) são conhecidas para todo \(i=1,..,p\). Define-se valores iniciais \(\theta^0 = (\theta_{1}^{0},\theta_{2}^{0},...,\theta_{p}^{0})\), assim sendo, na \(j\)-ésima iteração a cadeia encontra-se no estado \(\theta^{(j)}\) e a transição \(\theta^{(j+1)}\) pode ser descrita como segue:

1. Gere \(\theta_{1}^{(j+1)}\) de \(p(\theta_1|\theta_{2}^{(j)},...,\theta_{p}^{(j)})\)
2. Gere \(\theta_{2}^{(j+1)}\) de \(p(\theta_2|\theta_{1}^{(j+1)},\theta_{3}^{(j)},...,\theta_{p}^{(j)})\)
   
3. Repita os passos anteriores em \(\theta_i\) para \(i=3,...,p\), de forma que a última geração do \(j\)-ésimo passo é gerar \(\theta_{p}^{(j+1)}\) de: \[ p(\theta_p|\theta_{1}^{(j+1)},\theta_{2}^{(j+1)},...,\theta_{p-1}^{(j+1)}) \]

Nota: Distribuição a priori para variância

É proposto pelo exercício uma distribuição a priori Gama Inversa para \(\sigma^2*\) com parâmetros 0,01 e 0,01. No entanto, a distribuição \(%GI(\alpha,\beta)\) tem a média definida apenas para \(\alpha>2\). Outrora, utilizar a proposta \(\alpha=0,01\) e \(\beta=0,01\) conduz a média e variância negativas. Para situações em que a Gama Inversa for utilizada, os hiperparâmetros serão definidos como \(\alpha=3\) e \(\beta=20\), conduzindo a uma média de 10 e uma variância de 100, isto é, uma distribuição pouco informativa.

Distribuições condicionais completas

Distribuição condicional completa para \(\mu|.\)

\[ \mu|. \propto exp\{ -\frac{\phi}{2}\sum(x_i-\mu)^2 \}exp\{ -\frac{t_{0}^{-2}}{2}\mu^2 \}\\ \propto exp\{-\frac{\phi}{2}[\mu^2-2\mu\sum x_i] - -\frac{t_{0}^{-2}}{2}\mu^2 \}\\ \propto exp\{ -\frac{1}{2}[\mu^2(n\phi+t_{0}^{-2})-2\mu\phi n\bar{x}] \}\\ \propto exp\{-\frac{n\phi +t_{0}^{-2}}{2}(\mu^2-2\mu\frac{\phi n\bar{x}}{n\phi + t_{0}^{-2}}) \}\\ \implies N(\frac{\phi n\bar{x}}{\delta^2},\delta^2) \quad ; \quad \delta^2 = (n\phi+t_{0}^{-2})^{-1} \]

Distribuição condicional completa para \(\phi|.\) \[ \phi|. \propto \phi^{\frac{n}{2}}exp\{ -\frac{\phi}{2}[\sum (x_i-\mu)^2] \}\phi^{\frac{a}{2}-1}exp\{\frac{\phi b}{2}\}\\ \propto \phi^{\frac{n+a}{2}+1}exp\{ \phi[\frac{\sum(x_i-\mu)^2+b}{2}] \} \\ \implies \phi|. \sim G(\frac{n+a}{2},\frac{\sum(x_i-\mu)^2+b}{2}) \]

Algoritmo para o modelo normal

1. Defina os hiperparâmetros \(t_0\), \(a\) e \(b\)
2. Defina valores inicias para os parâmetros \(\mu\) e \(\phi\)
3. Gere da distribuição condicional completa de \(\mu|.\)
4. Gere da distribuição condicional completa de \(\phi|.\) utilizando os valores gerados em 3
5. Repita o processo 3 e 4 até obtido a convergência
6. Obtenha \(\sigma^2 = \frac{1}{\phi}\)

Simulação

Uma boa prática para testar um algoritmo é simular um cenário controlado onde todas as quantidades são conhecidas. Suponha a seguinte relação hierárquica.

\[ [\mathbf{X}|\mu,\sigma^2] \sim N(\mu=0;\phi^{-1}=1) \\ \mu \sim N(0;1000)\\ \phi \sim G(0,01;0,01) \]

Gera-se 100 valores da Normal Padrão proposta. Imaginando que os valores reais são desconhecidos, efetua-se o algoritmo de Gibbs para verificar a efetividade do mesmo, com 10 mil iterações. A figura 3 enuncia que o algoritmo foi proposto corretamente uma vez que a cadeia converge para o ponto verdadeiro de cada parâmetro. A convergência parece ter sido obtida, mas não entraremos em detalhes no momento, assim como o tratamento também não será efetuado no caso simulado.

Figura 3 - Análise das cadeias para \(\mu\) e \(\sigma^2\) para o modelo normal simulado (valor real em vermelho)

Dados reais

A figura 4 apresenta o Amostrador de Gibbs para os dados reais (amostra observada) ainda sem tratamento, com 10 mil iterações. Antes de realizar o tratamento efetivo, é necessário verificar se a cadeia convergiu. Verificou por inspeção gráfica que a convergência foi obtida rapidamente e que existe algum resíduo inicial que pode ser caraterizado como aquecimento.

Figura 4 - Cadeias não tratadas para \(\mu\) e \(\sigma^2\) para os dados reais

A tabela 3 apresenta o diagnóstico de Raftery e Lewis para as cadeias dos dados reais. O número ‘3746’ de iterações é considerado ‘mágico’, e é sempre o mínimo que o diagnóstico do pacote coda retorna. O aquecimento estimado é baixo e não há evidências de autocorrelação.

Tabela 3 - Diagnóstico de Raftery e Lewis

Baseado no critério, garantido que houve a convergência, a cadeia foi tratada retirando mil observações (mais do que o necessário) e não foi tomado espaçamento pois não havia necessidade. A figura 5 apresenta o resultado final para o modelo normal.

Figura 5 - Cadeias tratadas para \(\mu\) e \(\sigma^2\) para os dados reais

Com a cadeia convergida e tratada, é possível conduzir inferência através dos estimados que minimizam o erro quadrático médio e o absoluto, respectivamente sendo a média e a mediana. A tabela 4 apresenta tais estimativas tão bem como o intervalo de credibilidade HDP com probabiliade de 95% para cada parâmetro.

Tabela 4 - Inferência nos parâmetros

Metropolis-Hatings

Como o núcleo da distribuição alvo não é analítico, o algoritmo comumento utilizado é conhecido como Metropolis-Hastings (MH). Em linhas gererais, o algoritmo MH é um processo de aceitação-rejeição baseado em distribuições propostas. A ideia é parecida das práticas estudadadas em simulação de variáveis aleatórias. Temos uma informação parcial (núcleo da distribuição a posteriori) e gostaríamos de atingir a distribuição alvo com o auxílio de uma distribuição proposta.

Para o parâmetro \(\mu\), utilizou-se a proposta do passeio aleatório, já para o parâmetro \(\sigma^2\), não podemos utilizar a mesma abordagem pois o suporte do parâmetro são os positivos, logo uma distribuição lognormal foi utilizada como proposta. As distribuições condicionais completas são as seguintes:

\[ p(\mu|.) \propto \left[\prod(1+(\frac{x-\mu}{\sigma})^2\right]^{\frac{-(v+1)}{2}}\times exp\{-\frac{t_{0}^{-2}}{2}\mu^2\} \\ p(\sigma^2|.) \propto \sigma^{\frac{n}{2}}\left[\prod(1+(\frac{x-\mu}{\sigma})^2\right]^{\frac{-(v+1)}{2}}\times (\sigma^2)^{-(a+1)}exp\{ -b\sigma^2 \} \]

Logo, o algoritmo MH em uma iteração j pode ser escrito como:

1. Gere um valor proposto \(\mu_p\) por uma distribuição normal \(N(\mu_{j-1},\Delta_1)\)
2. Gere um valor proposto \(\sigma^2_p\) por uma distribuição lognormal \(LN(\sigma_{j-1}^{2},\Delta_1)\)
3. Calcule \(\alpha_1 = \min\{1,\frac{p(\mu_{p})}{p(\mu_{j-1})}\}\)
4. Calcule \(\alpha_2 = \min\{1, \frac{p(\sigma_{p}^{2})\times q(\sigma_{j-1}^{2})}{p(\sigma_{j-1}^{2})\times p(\sigma_{p}^{2})} \}\)
5. Gera-se dois valores \(u_1\) e \(u_2\) de uma distribuição uniforme(0,1)
6. Se \(u_1<\alpha_1\) então o valor proposto é aceito como novo valor \(\mu_j\)
7. Se \(u_2<\alpha_2\) então o valor proposto é aceito como novo valor \(\sigma_{j}^{2}\)
8. Para-se quando uma quantidade pré-definida de iterações \(j\) for atingida.

Simulação

Para trienar o algoritmo, como boa prática, suponha a seguinte relação hierárquica com os parâmetros conhecidos:

\[ [\mathbf{X}|\mu,\sigma^2] \sim t_5(0;1) \\ \mu \sim N(0;1000)\\ \sigma^2 \sim GI(3;20) \]

Agora suponha que a locação e a dispersão do modelo t tornaram-se desconhecidos. Tem-se por objetivo inferir a respeito destes parâmetros utilizando o algoritmo MH.

Figura 6 - Análise das cadeias para \(\mu\) e \(\sigma^2\) por Metropolis-Hastings para o modelo t simulado (valor real em vermelho)

Diferentemente do exercício anterior, onde o amostrador de Gibbs foi utilizados, vemos pela figura 6 que existe uma forte autocorrelação. Em linhas gerais, o algoritmo MH foi bem implementado, uma vez que a distribuição a posteriori para cada parâmetro converge para o valor verdadeiro.

Dados reais

A figura 7 apresenta o algoritmo de MH aplicado aos dados reais. Pelo traço da cadeia, há evidências visuais de que se obteve a convergência com a quantidade de iterações propostas. Por outro lado, como é esperado, a autocorrelação é observada para além do lag 120.

Figura 7 - Cadeias não tratadas para \(\mu\) e \(\sigma^2\) para os dados reais

Pela tabela 5, a proposta final para o modelo t irá partir de 601.000 iterações, com burn-in de 1000 e espaçamento de 100, retornando uma amostra efetiva de tamanho 3 mil.

Tabela 5 - Diagnóstico de Raftery e Lewis

A figura 8 apresenta as cadeias tratadas, para a locação, vemos uma concentração entre 5 e 6, já para a dispersão vemos uma curva assimétrica concentrada entre 2 e 4.

Figura 8 - Cadeias tratadas para \(\mu\) e \(\sigma^2\) para os dados reais

Tabela 6 - Inferência nos parâmetros

Modelo normal

Seguindo as mesmas definições apresentadas anteriormente, foi utilizado o Jags para ajustar um modelo normal. Primeiro será apresentado o caso simulado para a normal(0,1). Foram gerados inicialmente 3 cadeias com 10 mil iterações cada para verificação da convergência. A figura 9 apresenta os resultados e indica visualmente que com esta quantidade de iterações, a convergência é obtida.

Figura 9 - Cadeias simuladas para uma normal(0,1), não tratadas para \(\mu\) e \(\sigma^2\) (Valor verdadeiro em vermelho).

Dado que o Jags está bem ajustado para o caso simulado, parte-se para os dados reais. A figura 10 apresenta os resultados não tratados, a cadeia apresenta convergência visual rápida e baixa ou quase nenhuma autocorrelação.

Figura 10 - Cadeias não tratadas para \(\mu\) e \(\sigma^2\) para os dados reais

Pela tabela 7, verificamos que ao conduzir o algoritmo de Gibbs pelo Jags, é necessário apenas 2 de aquecimento e não há a necessidade de espaçamento. Logo, considerou-se as cadeias obtidas na figura 10 como as cadeias já tratadas.

Tabela 7 - Diagnóstico de Raftery e Lewis

Com a cadeia convergida e tratada, efetua-se o processo de inferência. A tabela 8 apresenta as principais medidas resumo, juntando com um intervalo de credebilidade de 95%.

Tabela 8 - Inferência nos parâmetros

Modelo t

O mesmo processo será efetuado também para o modelo supondo função de verossimilhança t.

A figura 12 apresenta evidências de convergência para as cadeias com respeito aos parâmetros investigados da distribuição t.

Figura 12 - Cadeias simuladas para uma \(t_5(0,1)\), não tratadas para \(\mu\) e \(\sigma^2\) (Valor verdadeiro em vermelho).

Dado que o Jags está bem ajustado para o caso simulado, parte-se para os dados reais. A figura 13 apresenta os resultados não tratados, a cadeia apresenta convergência visual rápida e baixa ou quase nenhuma autocorrelação.

Figura 13 - Cadeias não tratadas para \(\mu\) e \(\sigma^2\) para os dados reais

Pela tabela 9, verificamos que ao conduzir o algoritmo de Gibbs pelo Jags, é necessário apenas 2 de aquecimento e não há a necessidade de espaçamento. Logo, considerou-se as cadeias obtidas na figura 13 como as cadeias já tratadas.

Tabela 9 - Diagnóstico de Raftery e Lewis

Com a cadeia convergida e tratada, efetua-se o processo de inferência. A tabela 10 apresenta as principais medidas resumo, juntando com um intervalo de credebilidade de 95%.

Tabela 10 - Inferência nos parâmetros

Gelfand e Gosh

Suponha que seja de interesse calcular o poder preditivo de uma distribuição a posteriori para os parâmetros desconhecidos. Uma forma possível é utilizar o Critério de Gelfand e Gosh como informação. O Critério Preditivo a Posteriori (CPD) pode ser calculado como a soma da variabilidade (incerteza), que pode ser entendida como uma penalização, mais uma medida de bondade de ajuste.

A penalização P é calculada através da variância da previditiva a posteriori. Quanto menor esta quantidade, menos incerteza de previsão teremos.

A bondade do ajuste G é calculado através da distância entre o observação e a média das previsões, pode ser vista como uma medida de erro, quanto menor, melhor.

Logo, CPC = P + G. Formalmente,

\[ P = \sum_{i=1}^{n} V(y_i|z_i,\theta)\\ G = \sum_{i=1}^{n} [z_i - E[y_i|z_i,\theta)]^2 \]

Onde \(z_i\) é observado e \(y_i\) é previsto.

Deviance Information Criteria

Agora suponha que seja de interesse ter uma medida que indique a qualidade do ajuste de uma distribuição a posteriori para os parâmetros desconhecidos. Uma possível forma é através do Deviation Information Criteria (DIC). Para calcular, o DIC pode ser entendido como uma densidade a posteriori penalizada pela quantidade de parâmetros.

Formalmente, é possível calcular da seguinte forma:

\[ DIC = -2*L(\theta) + 2*pD\\ pD = \text{Densidade das médias dos parâmetros - Média das densidades dos parâmetros} \]

Entende-se que \(-2*L(\theta)\) indica a bondade do ajuste, logo, por ser negativo, quanto menor, mais bem ajustado o modelo está. A quantidade pD é uma medida que pode ser interpretada como incerteza, se a densidade das médias dos parâmetros for muito diferentes da média da densidade dos parâmetros, existe mais variabilidade a posteriori, logo penalizando o modelo. A nível de curiosidade, o BIC (Bayesian Information Critério) não é o critério mais utilizado no contexto bayesiano, como o nome sugere. O critério mais utilizado para qualidade de ajuste é o DIC.

A tabela 11 apresenta as quantidades obtidas, para a melhor compreensão, o Critério de Gelfand e Gosh foi dividido entre a Penalidade e o Ajuste. Verificou-se que para todas as abordagens, o Ajuste de Gelfand e Gosh foram numericamente próximos e os DIC’s também encontram-se em patamares simulares. No entanto, quanto a incerteza de previsão, verificou-se que a Normal por Gibbs foi a melhor enquanto o Modelo t por Jags foi o pior.

Tabela 11 - Critérios comparativos

Modelo normal por Gibbs

y <- c(6.08, 4.76, 9.27, 4.31, 4.15, 4.44, 4.87, 4.42, 5.54, 8.34,
       1.66, 7.22, 6.29, 3.81, 3.09, 6.06, 5.58, 2.89, 0.37, 6.42,
       2.34, 6.52, 5.58, 7.26, 3.30, 6.25, 0.30, 2.47, 8.61, 5.86,
       7.76, 5.13, 8.58, 6.71, 3.45, 6.90, 7.25, 2.99, 6.45, 7.51,
       4.48, 4.04, 7.69, 3.61, 7.53, 6.96, 3.84, 4.65, 2.55, 7.13)


n      <- length(y)
ybarra <- mean(y)


#====================================================
# Dados
#====================================================
n      <- length(y)
ybarra <- mean(y)
prec    <- 1/var(y)


#====================================================
# Priori
#====================================================
a0     <- 0.01
b0     <- 0.01
tau02  <- 10

#====================================================
# Amostragem de Gibbs
#====================================================
M      <- 10000      # N?mero de itera??es
mu_1     <- rep(NA,M)
phi    <- rep(NA,M)
mu_1[1]  <-  0         # Valor inicial
phi[1] <-  1         # Valor inicial

for (i in 2:M){
  # Gerando de phi
  somaq  <- sum( (y - mu_1[i-1])^2 )
  phi[i] <- rgamma(1,(n+a0)/2,(somaq+b0)/2)
  
  # Gerando de mu_1
  delta2 <- 1/(n*phi[i]+1/tau02)
  xi     <- n*ybarra*phi[i]*delta2
  mu_1[i]  <- rnorm(1,xi,sqrt(delta2))
}

sigma2_1 <- 1/phi

indc <- 1:M
par(mfrow=c(3,2))
ts.plot(mu_1[indc],xlab="Iteração",ylab=expression(mu))
ts.plot(sigma2_1[indc],xlab="Iteração",ylab=expression(sigma[2]))
hist(mu_1[indc],prob=T,col="darkblue",main="",nclass=50,xlab=expression(mu),ylab="Densidade")
hist(sigma2_1[indc],prob=T,col="darkblue",main="",nclass=50,xlab=expression(sigma[2]),ylab="Densidade")
acf(mu_1[indc],xlab="Defasagem",ylab=expression(paste("ACF:",mu)),main="")
acf(sigma2_1[indc],xlab="Defasagem",ylab=expression(paste("ACF:",sigma[2])),main="")

Modelo t por metropolis

set.seed(02021994)
x  <- c(6.08, 4.76, 9.27, 4.31, 4.15, 4.44, 4.87, 4.42, 5.54, 8.34,
       1.66, 7.22, 6.29, 3.81, 3.09, 6.06, 5.58, 2.89, 0.37, 6.42,
       2.34, 6.52, 5.58, 7.26, 3.30, 6.25, 0.30, 2.47, 8.61, 5.86,
       7.76, 5.13, 8.58, 6.71, 3.45, 6.90, 7.25, 2.99, 6.45, 7.51,
       4.48, 4.04, 7.69, 3.61, 7.53, 6.96, 3.84, 4.65, 2.55, 7.13)


n      <- length(x)
ybarra <- mean(x)


#====================================================
# Dados
#====================================================
n      <- length(x)
xbarra <- mean(x)
prec    <- 1/var(x)

#====================================================
# Priori
#====================================================
a0     <- 3
b0     <- 3
tau02  <- 1000
v      <-  5

#====================================================
# Amostragem de Gibbs
#====================================================
M         <- 200000      # N?mero de itera??es
mu        <- rep(NA,M)
sigma2    <- rep(NA,M)
mu[1]     <-  runif(1)         # Valor inicial
sigma2[1] <-  runif(1)         # Valor inicial

delta     = 3
delta2    = 3
contador  = 0
contador2 = 0
#Algoritmo Gibbs sampling com passo de Metropolis-Hastings
i=2
#Inicializando o vetor de parâmetros e definindo hiperparâmetros




for (i in 2:M) {
   
  #gerando valores aleat?rios
  mu.p  = rnorm(1,mu[i-1],delta) 
  sigma2.p = rlnorm(1,sigma2[i-1],delta2)  
  
  #Gerando amostras de mu via Metropolis Hastings
  numerador   = prod(1/sqrt(sigma2.p)*(1+1/(sigma2.p*v)*(x-mu.p)^2)^(-(v+1)/2))*exp(-1/tau02*mu.p^2)
  denominador = prod(1/sqrt(sigma2[i-1])*(1+1/(sigma2[i-1]*v)*(x-mu[i-1])^2)^(-(v+1)/2))*exp(-1/tau02*mu[i-1]^2)
  alpha       = numerador/denominador
  teste       = runif(1,0,1)
  if (teste<=alpha) {mu[i]=mu.p;contador=contador+1} else {mu[i]=mu[i-1]}
  
  #Gerando amostras de phi via Metropolis Hastings
  
  numerador   = prod(1/sqrt(sigma2.p)*(1+1/(sigma2.p*v)*(x-mu.p)^2)^(-(v+1)/2))*sigma2.p^(-a0-1)*exp(-b0*1/sigma2.p)*dlnorm(sigma2[i-1],sigma2[i-1],delta2)  
  denominador = prod(1/sqrt(sigma2[i-1])*(1+1/(sigma2[i-1]*v)*(x-mu[i-1])^2)^(-(v+1)/2))*sigma2[i-1]^(-a0-1)*exp(-b0*1/sigma2[i-1])*dlnorm(sigma2.p,sigma2[i-1],delta2)  
  alpha       = numerador/denominador
  teste       = runif(1,0,1)
  if (teste<=alpha) {sigma2[i]=sigma2.p;contador2=contador2+1} else {sigma2[i]=sigma2[i-1]}
  
}

tx1 = contador/M
tx2 = contador2/M


indc <- 1:M
par(mfrow=c(3,2))
ts.plot(mu[indc],xlab="Iteração",ylab=expression(mu))
ts.plot(sigma2[indc],xlab="Iteração",ylab=expression(sigma[2]))
hist(mu[indc],prob=T,col="darkblue",main="",nclass=50,xlab=expression(mu),ylab="Densidade")
hist(sigma2[indc],prob=T,col="darkblue",main="",nclass=50,xlab=expression(sigma[2]),ylab="Densidade")
acf(mu[indc],xlab="Defasagem",ylab=expression(paste("ACF:",mu)),main="",lag.max=120)
acf(sigma2[indc],xlab="Defasagem",ylab=expression(paste("ACF:",phi)),main="",lag.max=120)

Modelo normal por Jags

set.seed(02021994)
x = rt(100,df=5)

n = length(x)

model_code = 'model
{
  # Vero
  for (i in 1:n) {
  x[i] ~ dt(mu,phi,5)
  }
  # Prioris
  mu ~ dnorm(0,0.001)
  phi ~ dgamma(0.01,0.01)
  sigma2 <- 1/phi
}'

model_parameters =  c("mu", "sigma2")

model_data = list(n = n, x = x)

model_amostra = jags(data = model_data,
                 parameters.to.save = model_parameters,
                 model.file=textConnection(model_code),
                 n.chains=3, 
                 n.iter=10000,
                 n.burnin=0, 
                 n.thin=1,
                 DIC=FALSE,
                 progress.bar   = "none")

Modelo t por Jags

set.seed(02021994)
x = c(6.08, 4.76, 9.27, 4.31, 4.15, 4.44, 4.87, 4.42, 5.54, 8.34,
       1.66, 7.22, 6.29, 3.81, 3.09, 6.06, 5.58, 2.89, 0.37, 6.42,
       2.34, 6.52, 5.58, 7.26, 3.30, 6.25, 0.30, 2.47, 8.61, 5.86,
       7.76, 5.13, 8.58, 6.71, 3.45, 6.90, 7.25, 2.99, 6.45, 7.51,
       4.48, 4.04, 7.69, 3.61, 7.53, 6.96, 3.84, 4.65, 2.55, 7.13)



n = length(x)

model_code = 'model
{
  # Vero
  for (i in 1:n) {
  x[i] ~ dt(mu,phi,5)
  }
  # Prioris
  mu ~ dnorm(0,0.001)
  phi ~ dgamma(0.01,0.01)
  sigma2 <- 1/phi
}'

model_parameters =  c("mu", "sigma2")

model_data = list(n = n, x = x)

model_amostra = jags(data = model_data,
                 parameters.to.save = model_parameters,
                 model.file=textConnection(model_code),
                 n.chains=1, 
                 n.iter=10000,
                 n.burnin=1000, 
                 n.thin=1,
                 DIC=FALSE,
                 progress.bar   = "none")

Modelo t por Gibbs

set.seed(02021994)
x = rt(100,5)

n      <- length(x)
ybarra <- mean(x)


#====================================================
# Dados
#====================================================
n      <- length(x)
xbarra <- mean(x)
prec    <- 1/var(x)


#====================================================
# Priori
#====================================================
a0     <- 0.01
b0     <- 0.01
tau02  <- 10
v      <-  5


#====================================================
# Amostragem de Gibbs
#====================================================
M      <- 10000      # N?mero de itera??es
mu     <- rep(NA,M)
phi    <- rep(NA,M)
y      <- matrix(NA,nrow=M,ncol=n)
y[1,]  <-  0 
mu[1]  <-  0         # Valor inicial
phi[1] <-  1         # Valor inicial


for (i in 2:M){
  
  # Gerando de y
  for (j in 1:n){
  y[i,j] <- rgamma(1,(v+1)/2,phi[i-1]*(x[j]-mu[i-1])^2/2 + v/2)
  }

  # Gerando de phi
  soma_phi   <- sum(y[i,]*(x-mu[i-1])^2)/2 + b0/2
  phi[i]     <- rgamma(1,(n+a0)/2,soma_phi)
  
  # Gerando de mu
  delta2  <- 1/(sum(y[i,])*phi[i]+1/tau02)
  xi     <- phi[i]*delta2*sum(x*y[i,])
  mu[i]  <- rnorm(1,xi,sqrt(delta2))
  
}

sigma2 <- 1/phi

indc <- 1:M
par(mfrow=c(3,2))
ts.plot(mu[indc],xlab="Iteração",ylab=expression(mu))
abline(h=0,col=2)
ts.plot(phi[indc],xlab="Iteração",ylab=expression(phi))
abline(h=1,col=2)
hist(mu[indc],prob=T,col="darkblue",main="",nclass=50,xlab=expression(mu),ylab="Densidade")
abline(v=0,col=2)
hist(phi[indc],prob=T,col="darkblue",main="",nclass=50,xlab=expression(phi),ylab="Densidade")
abline(v=1,col=2)
acf(mu[indc],xlab="Defasagem",ylab=expression(paste("ACF:",mu)),main="")
acf(phi[indc],xlab="Defasagem",ylab=expression(paste("ACF:",phi)),main="")

Medidas de comparação

set.seed(02021994)

y1.pred= c()
y1=matrix(NA,length(mu_1),length(x))

y2.pred= c()
y2=matrix(NA,length(mu_2),length(x))

y3.pred= c()
y3=matrix(NA,length(mu_3),length(x))

y4.pred= c()
y4=matrix(NA,length(mu_4),length(x))


for (j in 1:length(x)){
for(i in 1:length(mu_1)){y1.pred[i]=rnorm(1,mean = mu_1[i],sd = sqrt(sigma2_1[1]))}
y1[,j] = y1.pred
}

for (j in 1:length(x)){
for(i in 1:length(mu_2)){y2.pred[i]=rnorm(1,mean = mu_2[i],sd = sqrt(sigma2_2[1]))}
y2[,j] = y2.pred
} 

for (j in 1:length(x)){
  for(i in 1:length(mu_3)){y3.pred[i]=rnorm(1,mean = mu_3[i],sd = sqrt(sigma2_3[1]))}
y3[,j] = y3.pred
} 

for (j in 1:length(x)){
for(i in 1:length(mu_4)){y4.pred[i]=rnorm(1,mean = mu_4[i],sd = sqrt(sigma2_4[1]))}
y4[,j] = y4.pred
} 


var1 = sum(apply(y1,2,var))
var2 = sum(apply(y2,2,var))
var3 = sum(apply(y3,2,var))
var4 = sum(apply(y4,2,var))

termo1 = sum((x-apply(y1,2,mean))^2)
termo2 = sum((x-apply(y2,2,mean))^2)
termo3 = sum((x-apply(y3,2,mean))^2)
termo4 = sum((x-apply(y4,2,mean))^2)


gg1 = sum(apply(y1,2,var)) + sum((x-apply(y1,2,mean))^2)
gg2 = sum(apply(y2,2,var)) + sum((x-apply(y2,2,mean))^2)
gg3 = sum(apply(y3,2,var)) + sum((x-apply(y3,2,mean))^2)
gg4 = sum(apply(y4,2,var)) + sum((x-apply(y4,2,mean))^2)

teste = dnorm(x,mu_1,sigma2_1)


dic <- function(data, likelihood, postSamples){
  
  logLikelihood <- function(theta) sum(likelihood(data, theta[1], theta[2], log=T))
  
  thetaBayes <- colMeans(postSamples)
  
  pDIC <- 2*(logLikelihood(thetaBayes) - mean(apply(postSamples, 1, logLikelihood) ))
  dic <- -2*logLikelihood(thetaBayes) + 2*pDIC
  
  return(dic)
}

dic1 = dic(x,dnorm,cbind(mu_1,sigma2_1))
dic2 = dic(x,dnorm,cbind(mu_2,sigma2_2))
dic3 = dic(x,dnorm,cbind(mu_3,sigma2_3))
dic4 = dic(x,dnorm,cbind(mu_4,sigma2_4))

dt = cbind(c(var1,var2,var3,var4),c(termo1,termo2,termo3,termo4),c(dic1,dic2,dic3,dic4))

colnames(dt) = c("Penalidade de Gelfand-Gosh","Ajuste de Gelfand-Gosh","DIC")
rownames(dt) = c("Modelo Normal por Gibbs","Modelo t por Metropolis-Hastings","Modelo Normal por Jags","Modelo t por Jags")

kable(dt) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Monte carlo via Cadeias de Markov com Saltos Reversíveis

y <-c(6.08, 4.76, 9.27, 4.31, 4.15, 4.44, 4.87, 4.42, 5.54, 8.34,
      1.66, 7.22, 6.29, 3.81, 3.09, 6.06, 5.58, 2.89, 0.37, 6.42,
      2.34, 6.52, 5.58, 7.26, 3.30, 6.25, 0.30, 2.47, 8.61, 5.86,
      7.76, 5.13, 8.58, 6.71, 3.45, 6.90, 7.25, 2.99, 6.45, 7.51,
      4.48, 4.04, 7.69, 3.61, 7.53, 6.96, 3.84, 4.65, 2.55, 7.13)
n <- length(y)

#====================================================
# negativo da log-verossimilhancas em theta[1:2] (real)
logflognormal <- function(theta){
    -sum(dnorm(y,theta[1],sqrt(theta[2]),log=T))
}

logfweibull <- function(theta){
    -sum(dst(y,xi = theta[1],omega=theta[2],alpha=0,log=T))
}

#====================================================
# Estimacao de maxima verossimilhanca
# Minimizacao nao linear
#----------------------------------------------------
# Modelo lognormal 
thetalognormal <- nlm(logflognormal,c(3.35,log(1.4)))$estimate

# Modelo Weibull
thetaweibull   <- nlm(logfweibull,log(c(0.85,55)))$estimate

#----------------------------------------------------
# Explorando a verossimilhanca

# Modelo lognormal 
mu     <- seq(2.5,7.5,length=30)
sig2   <- seq(0.5,5.5,length=30)
sig2   <- log(sig2)
likeLN <- matrix(0,30,30)
for(i in 1:30)
    for(j in 1:30)
        likeLN[i,j] <- exp(-logflognormal(c(mu[i],sig2[j])))

# Modelo Weibull
gamma  <- seq(0.5,1.2,length=30)
gamma  <- log(gamma)
delta  <- seq(20,100,length=30)
delta  <- log(delta);
likeWB <- matrix(0,30,30)
for(i in 1:30)
    for(j in 1:30)
        likeWB[i,j] <- exp(-logfweibull(c(gamma[i],delta[j])))

#====================================================
# Amostragem por importancia
# Proposta normal bivariada com media no estimador de
# maxima verossimilhanca e desvio igual ao 'range/4'

set.seed(12345)
M <- 5000 # numero de valores propostos

#----------------------------------------------------
# Para calculo dos pesos (proposta)
dnormbiv <- function(theta,mu,sig2){
    -0.5*sum((theta-mu)^2/sig2)
}
#----------------------------------------------------

# Modelo lognormal 
mediaLNprop <- thetalognormal
sdLNprop    <- c(diff(range(mu)),diff(range(sig2)))/4
thetaLNprop <- matrix(rnorm(M,mediaLNprop,sdLNprop),M,2,byrow=T)

wLN <- rep(NA,M)
for(i in 1:M)
    wLN[i] <- -( logflognormal(thetaLNprop[i,])
                 + dnormbiv(thetaLNprop[i,],mediaLNprop,sdLNprop) )

# Reamostragem 
indLN   <- sample(1:M,replace=T,prob=exp(wLN),size=M)
thetaLN <- thetaLNprop[indLN,]

# Salvando variancia da posteriori
thetaLNmediapost <- apply(thetaLN,2,mean) 
cst              <- 2 # fator para inflar a variancia
thetaLNcovarpost <- cst*cov(thetaLN)

# Modelo Weibull
mediaWBprop <- thetaweibull
sdWBprop    <- c(diff(range(gamma)),diff(range(delta)))/4
thetaWBprop <- matrix(rnorm(M,mediaWBprop,sdWBprop),M,2,byrow=T)

wWB <- rep(NA,M)
for(i in 1:M)
    wWB[i] <- -( logfweibull(thetaWBprop[i,]) + dnormbiv(thetaWBprop[i,],mediaWBprop,sdWBprop) )

# Reamostragem 
indWB   <- sample(1:M,replace=T,prob=exp(wWB),size=M)
thetaWB <- thetaWBprop[indWB,]

# Salvando variancia da posteriori
thetaWBmediapost <- apply(thetaWB,2,mean) 
cst              <- 2 # fator para inflar a variancia
thetaWBcovarpost <- cst*cov(thetaWB)

#====================================================
# Monte Carlo via cadeias de Markov com saltos
# reversiveis
#====================================================
# vetor de medias dos dois modelos
mm    <- rbind(thetaLNmediapost,thetaWBmediapost )

# matrizes da decomposicao de Cholesky da 
# covariancia dos dois modelos
Vchol      <- array(NA,c(2,2,2))
Vchol[1,,] <- t(chol(thetaLNcovarpost))
Vchol[2,,] <- t(chol(thetaWBcovarpost))

# verossimilhanca para cada modelo
loglikeprior <- function(m,theta){
    if(m==1)
        return(-logflognormal(theta))
    if(m==2)
        return(-logfweibull(theta))
}

# Probabilidade para mudar de modelos
J1 <- matrix(0.5,2,2)

# Transicao independente normal multivariada
# com momentos da posteriori
dnorm2 <- function(x,mu,logA,A1){
    -logA-0.5*sum((A1%*%(x-mu))^2)
}
# Entradas da funcao proposta
A1   <- Vchol
logA <- NULL
for(i in 1:2){
    A1[i,,] <- solve(Vchol[i,,])
    logA    <- c(logA,sum(log(diag(Vchol[i,,]))))
}

#----------------------------------------------------
# Inicio do algoritmo Monte Carlo via cadeias
# de Markov com saltos reversiveis

M           <- 5000   # Numero de iteracoes
m           <- 2      # comecando no modelo 2
theta       <- c(0,1) # inicializando os parametros (M_2)
amostra     <- matrix(NA,M,3)
amostra[1,] <- c(m,theta)
for(i in 1:M){
    # escolho modelo entre 1 2 de acordo com matriz J_1
    mprop <- sample(1:2,size=1,prob=J1[m,])

    # proposta: Normal multivariada com media e variancia fixos
    param <- mm[mprop,] + Vchol[mprop,,]%*%rnorm(2) 
    # verossimilhanca no proposto sobre anterior
    vero <- loglikeprior(mprop,param)-loglikeprior(m,theta) 
    prop <- ( dnorm2(theta,mm[m,],logA[m],A1[m,,])
              - dnorm2(param,mm[mprop,],logA[mprop],A1[mprop,,]) ) 
    A    <- vero + prop
    prob <- min(0,A)
    if( log(runif(1)) < prob ){
        m     <- mprop
        theta <- param
    }
    amostra[i,] <- c(m,theta) 
}
proporcao <- table(amostra[,1])/M

modelo1 <- (amostra[,1]==1)
modelo2 <- !modelo1