7. Amostradores Gibbs

Este capítulo abrange os amostradores Gibbs de dois estágios e de vários estágios. Embora o primeiro seja um caso especial do último, o amostrador de dois estágios tem propriedades de convergência superiores e se aplica naturalmente a uma ampla gama de modelos estatísticos que não exigem a generalidade do amostrador de vários estágios. No entanto, o amostrador Gibbs de vários estágios desfruta de muitas propriedades de otimização e ainda pode ser considerado o burro de carga do mundo MCMC.

Após a introdução na Seção 7.1 com algum histórico, desenvolvemos o amostrador Gibbs de dois estágios na Seção 7.2, passando para o amostrador Gibbs de vários estágios na Seção 7.3. O amostrador de Gibbs é particularmente adequado para lidar com experimentos com dados ausentes e modelos com variáveis latentes, conforme mostrado na Seção 7.4.

Embora façamos uso de modelos hierárquicos ao longo do capítulo, nos concentramos em seu processamento na Seção 7.5. A Seção 7.6 analisa vários tópicos adicionais, como Rao-Blackwellization, reparametrização e o efeito do uso de a prioris impróprias.



7.1 Introdução

O Capítulo 6 descreveu alguns princípios para simulação baseada em Cadeias de Markov, bem como algumas direções de implementação, incluindo o algoritmo genérico de passeio aleatório Metropolis-Hastings. Este capítulo estende o escopo dos algoritmos MCMC estudando outra classe de métodos MCMC agora comuns, chamada amostragem de Gibbs. O apelo desses algoritmos específicos é que primeiro eles coletam a maior parte de sua calibração a partir da densidade do alvo e, segundo, nos permitem quebrar problemas complexos, como distribuições de alvos de alta dimensão, para os quais um algoritmo Metropolis-Hastings de passeio aleatório é quase impossível construir, em uma série de problemas mais fáceis, como uma sequência de alvos de pequena dimensão. Pode haver ressalvas a essa simplificação, pois a sequência de problemas simples pode levar muito tempo para convergir, mas a amostragem de Gibbs é, no entanto, um candidato interessante ao lidar com um novo problema.

O nome amostragem de Gibbs vem do artigo de referência de Geman and Geman (1984), que aplicou pela primeira vez um amostrador de Gibbs em um campo aleatório de Gibbs. Para o bem ou para o mal, ele travou, apesar desse elo fraco. De fato, é um caso especial do algoritmo Metropolis-Hastings conforme detalhado em Robert and Casella (2004). O trabalho de Geman and Geman (1984), construído sobre o de Metropolis et al. (1953), Hastings (1970) e Peskun (1973), influenciaram Gelfand and Smith (1990) a escrever um artigo que despertou um novo interesse em métodos bayesianos, computação estatística, algoritmos e processos estocásticos por meio do uso de algoritmos de computação, como o amostrador Gibbs e o algoritmo Metropolis-Hastings.

É interessante ver, em retrospecto, que artigos anteriores como Tanner and Wong (1987) e Besag and Clifford (1989) propuseram soluções semelhantes, mas não receberam a mesma resposta da comunidade estatística.


7.2 O amostrador Gibbs de dois estágios

O amostrador Gibbs de dois estágios cria uma Cadeia de Markov a partir de uma distribuição conjunta da seguinte maneira. Se duas variáveis aleatórias \(X\) e \(Y\) têm densidade conjunta \(f(x,y)\), com densidades condicionais correspondentes \(f_{Y|X}\) e \(f_{X|Y}\), o amostrador Gibbs de dois estágios gera uma Cadeia de Markov \((X_t,Y_t)\) de acordo com as seguintes etapas:


Algoritmo 7: Amostrador Gibbs de dois estágios
Escolher \(X_0=x_0\)
Para \(t=1,2,\cdots\) gerar
1. \(Y_t\sim f_{Y|X}(\cdot \, | x_{t-1})\);
2. \(X_t\sim f_{X|Y}(\cdot \, | y_{t})\).


O algoritmo 7 é então direto de implementar, desde que a simulação de ambos os condicionais seja viável. Quando \(f(x,y)\) está disponível na forma fechada, até uma constante de normalização, também estão \(f_{Y|X}\) e \(f_{X|Y}\). Portanto, se não for possível simular diretamente a partir desses condicionais, podem ser usadas aproximações Monte Carlo ou MCMC, conforme desenvolvido na Seção 7.6.3.

Também é possível ver por que, se \((X_t,Y_t)\) é distribuído segundo \(f\), então \((X_{t+1},Y_{t+1})\) também é, porque ambas as etapas da iteração \(t\) usam a simulação dos verdadeiros condicionais. A convergência da Cadeia de Markov e, portanto, do algoritmo é garantida, a menos que os suportes das condicionais não sejam conectados.


Exemplo 7.1

Para começar com uma ilustração óbvia, considere o modelo normal bivariado \[\begin{equation} \tag{7.1} (X,Y) \sim N_2\left( 0,\begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix} \right), \end{equation}\] para o qual o amostrador Gibbs é:

Dado \(x_t\), gere \[\begin{equation} \begin{array}{rcl} Y_{t+1} | x_t & \sim & N(\rho \, x_t ,1-\rho^2), \\ X_{t+1} | y_{t+1} & \sim & N(\rho \, y_{t+1},1-\rho^2)\cdot \end{array} \end{equation}\]

A subcadeia \((X_t)_t\) então satisfaz \[\begin{equation} X_{t+1} | X_{t}=x_t \sim N(\rho^2 \, x_{t},1-\rho^4), \end{equation}\] e uma recursão mostra que \[\begin{equation} X_{t} | X_{0}=x_0 \sim N(\rho^{2t} \, x_{0},1-\rho^{4t}), \end{equation}\] que de fato converge para \(N(0,1)\) conforme \(t\) vai para o infinito.


Conforme ilustrado pelo exemplo acima, a sequência \((X_t,Y_t)\), \(t = 1,\cdots,T\), produzido por um amostrador de Gibbs converge para a distribuição conjunta \(f\) e, como consequência, ambas as sequências \((X_t)_t\) e \((Y_t)_t\) convergem para suas respectivas distribuições marginais.

Talvez a principal razão pela qual o amostrador de Gibbs se tornou tão popular na década de 1990 como o algoritmo MCMC de referência é que ele era o complemento computacional perfeito para modelos hierárquicos, que então começavam a ser seriamente investigados. Conforme detalhado e justificado na Seção 7.5, um modelo hierárquico especifica uma distribuição conjunta como camadas sucessivas de distribuições condicionais. O exemplo a seguir dá uma primeira olhada nos modelos hierárquicos.


Exemplo 7.2

Considerando o par de distribuições \[ X \, | \, \theta \sim Binomial(n,\theta), \qquad \theta\sim Beta(\alpha,\beta), \] leva à distribuição conjunta \[ f(x,\theta)=\binom{n}{x}\dfrac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\, x+\alpha-1}(1-\theta)^{\, n-x+\beta-1}\cdot \] 

set.seed(5832)
Nsim=5000 #initial values
n=15
a=3 #alpha
b=7 #beta
X=T=array(0,dim=c(Nsim,1)) #init arrays
T[1]=rbeta(1,a,b) #init chains
X[1]=rbinom(1,n,T[1])
for (i in 2:Nsim){ #sampling loop
 X[i]=rbinom(1,n,T[i-1])
 T[i]=rbeta(1,a+X[i],n-X[i]+b)
}
par(mfrow=c(1,2),mar=c(4,4,2,1))
hist(X, xlab="X", ylab="Marginal density", main = "", freq=FALSE)
grid(); box()
betabi=function(x,a,b,n){
  (gamma(n+1)*gamma(x+a)*gamma(n-x+b)*gamma(a+b))/
    (gamma(x+1)*gamma(n-x+1)*gamma(n+a+b)*gamma(a)*gamma(b))}
curve(betabi(x,a,b,n),col="orange",lwd=3, add=TRUE)
hist(T, xlab=expression(theta), ylim=c(0,3), ylab="Marginal density", main = "", freq=FALSE)
grid(); box()
curve(dbeta(x,a,b), col="orange",lwd=3, add=TRUE)

Figura. 7.1. Histogramas de distribuições marginais do amostrador de Gibbs com base em 5000 iterações do Algoritmo 7 para \(n = 15,\alpha = 3,\beta = 7\). A verdadeira distribuição marginal de \(\theta\) é \(Beta(\alpha,\beta)\) e a distribuição marginal de \(X\) é Beta-Binomial.

A distribuição condicional correspondente de \(X|\theta\) foi mostrada acima, enquanto \[ \theta | X\sim Beta(x + a, n-x + b)\cdot \] O amostrador de Gibbs associado pode ser implementado e sua saída é ilustrada na Figura 7.1 para cada marginal. Como este é um exemplo de brinquedo, as marginais de forma fechada estão disponíveis e, portanto, produzidas no topo dos histogramas mostram um bom ajuste para ambas as amostras Gibbs.



Exemplo 7.3

Considere a distribuição a posteriori de \((\theta,\sigma^2)\), associada ao modelo conjunto \[\begin{equation} \tag{7.2} \begin{array}{rcl} X_i & \sim & N(\theta,\sigma^2), \qquad i=1,\cdots,n, \\ \theta & \sim & N(\theta_0,\tau^2), \quad \sigma^2\sim \mathcal{IG}(\alpha,\beta), \end{array} \end{equation}\] onde \(\mathcal{IG}(\alpha,\beta)\) é a distribuição gama invertida, isto é, a distribuição do inverso de uma variável gama, com densidade \[ \dfrac{1}{x\Gamma(\alpha)}\left(\dfrac{\beta}{x}\right)^\alpha e^{-\beta/x}, \] com \(\theta_0\), \(\tau^2\), \(\alpha\) e \(\beta\) espacificados.

Escrevendo \(x=(x_1,\cdots,x_n)\), a distribuição a posteriori de \((\theta,\sigma^2)\) é dada por \[\begin{equation} \tag{7.3} f(\theta,\sigma^2| x) \approx \dfrac{e^{-\sum_i (x_i-\theta)^2/(2\sigma^2)}}{\big(\sigma^2\big)^{n/2}} \dfrac{e^{-(\theta-\theta_0)^2/(2\tau^2)}}{\tau} \dfrac{e^{-1/(\beta\sigma^2)}}{\big(\sigma^2\big)^{\alpha+1}}, \end{equation}\] a partir do qual podemos obter os condicionais completos de \(\theta\) e \(\sigma^2\).

Observe que esta não é uma configuração conjugada regular em que a integração em \(\theta\) ou \(\sigma^2\) nesta densidade, não produz uma densidade padrão.

Temos então que \[\begin{equation} \tag{7.4} \begin{array}{rcl} \pi(\theta |x,\sigma^2) & \approx & \displaystyle \exp\Big(-\sum_i (x_i-\theta)^2/(2\sigma^2)\Big)\times \exp\Big( -(\theta-\theta_0)^2/(2\tau^2\sigma^2) \Big), \\ \pi(\sigma^2|x,\theta) & \approx & \displaystyle (\sigma^2)^{-\frac{n+2\alpha+3}{2}}\times \exp\left(-\dfrac{1}{2\sigma^2}\Big(\sum_i(x_i-\theta)^2+(\theta-\theta_0)^2/\tau^2+2/\beta\Big)\right)\cdot \end{array} \end{equation}\]

Estas densidades correspondem à \[ \theta|x,\sigma^2 \sim N\left( \dfrac{\sigma^2}{\sigma^2+n\tau^2}\theta_0+\dfrac{n\tau^2}{\sigma^2+n\tau^2}\overline{x},\dfrac{\sigma^2\tau^2}{\sigma^2+n\tau^2}\right) \] e \[ \sigma^2|x,\theta \sim \mathcal{IG}\left( \dfrac{n}{2}+\alpha,\dfrac{1}{2}\sum_i(x_i-\theta)^2+\beta\right), \] onde \(\overline{x}\) é a média empírica das observações, como as distribuições condicionais completas a serem usadas em um amostrador de Gibbs.

Um estudo sobre o metabolismo em mulheres de 15 anos produziu os seguintes dados, denotados por \(x\), correspondentes à ingestão de energia, medida em megajoules, durante um período de 24 horas.

Usando o modelo normal acima, com \(\theta\) correspondendo à ingestão média real de energia, o amostrador de Gibbs pode ser implementado, onde theta0, tau2, a e b são valores especificados. 

x=c(91,504,557,609,693,727,764,803,857,929,970,1043,1089,1195,1384,1713)
set.seed(56721)
Nsim=5000
xbar=mean(x)
a=b=3
tau2=10
theta0=5
sh1=(n/2)+a
sigma2=theta=rep(0,Nsim) #init arrays
ra1=(1/2)*(sum((x-theta0)^2))+b
sigma2[1]=1/rgamma(1,shape=sh1,scale=ra1) #init chains
B=sigma2[1]/(sigma2[1]+n*tau2)
theta[1]=rnorm(1,mean=B*theta0+(1-B)*xbar,sd=sqrt(tau2*B))
for (i in 2:Nsim){
 B=sigma2[i-1]/(sigma2[i-1]+n*tau2)
 theta[i]=rnorm(1,mean=B*theta0+(1-B)*xbar,sd=sqrt(tau2*B))
 ra1=(1/2)*(sum((x-theta[i])^2))+b
 sigma2[i]=1/rgamma(1,shape=sh1,scale=ra1)
}
par(mfrow=c(1,2),mar=c(4,4,2,1))
hist(log(theta), xlab=expression(log(theta)), ylab="Marginal density", main = "", freq=TRUE) 
grid(); box()
hist(log(sqrt(sigma2)), xlab=expression(log(sigma)), ylab="Marginal density", main = "", 
     freq=TRUE, col="orangered")
grid(); box()

Figura. 7.2. Histogramas de distribuições a posteriori marginais de \(\log(\theta)\) e \(\log(\sigma)\) do amostrador de Gibbs, com base em 5.000 iterações, com \(a = b = 3\), \(\tau^2 = 10\) e \(\theta_0 = 5\).

mean(theta)
## [1] 870.5
mean(sqrt(sigma2))
## [1] 0.0003088007

As médias a posteriori de \(\theta\) e \(\sigma\) são 870.5 e 0.0002256299. Histogramas das distribuições a posteriori de \(\log(\theta)\) e \(\log(\sigma)\) são dados na Figura 7.2.


Queremos salientar que reconhecer as distribuições condicionais completas de uma distribuição conjunta não é tão difícil. Por exemplo, a distribuição a posteriori proporcional a (7.3) é obtida multiplicando as densidades na especificação (7.2).

Para encontrar uma condicional completa, ou seja, a distribuição condicional de um parâmetro condicional a todos os outros, precisamos apenas selecionar todos os termos da distribuição conjunta que envolvem esse parâmetro.

Por exemplo, de (7.3), vemos que \[ f(\theta|\sigma^2,x) \approx \left( \dfrac{1}{(\sigma^2)^{n/2}} \exp\Big(-\sum_i (x_i-\theta)^2/(2\sigma^2) \Big)\right)\times \left( \dfrac{1}{\tau^2}\exp\Big( -(\theta-\theta_0)^2/(2\tau^2)\Big)\right) \] e \[ f(\sigma^2|\theta,x) \approx \left( \dfrac{1}{(\sigma^2)^{n/2}} \exp\Big(-\sum_i (x_i-\theta)^2/(2\sigma^2) \Big)\right)\times \left( \dfrac{1}{(\sigma^2)^{\alpha+1}}\exp\Big( 1/(\beta\sigma^2)\right)\cdot \]

Deve então ser fácil ver que a distribuição condicional completa de \(\sigma^2\) será uma distribuição gama invertida, conforme definido acima, consulte também o Exercício 7.19. Para \(\theta\), embora haja um pouco mais de álgebra envolvida na derivação, podemos reconhecer que a condicional completa será normal. Veja o Exercício 7.20 para uma ilustração com uma hierarquia maior.


7.3 O amostrador Gibbs de vários estágios

Há uma extensão natural do amostrador Gibbs de dois estágios para o amostrador Gibbs geral de vários estágios. Suponha que, para algum \(p > 1\), a variável aleatória \(X\in\mathcal{X}\) pode ser escrita como \(X = (X_1,\cdots,X_p)\), onde os \(X_i\)’s são componentes unidimensionais ou multidimensionais.

Além disso, suponha que podemos simular a partir das densidades condicionais correspondentes \(f_1,\cdots,f_p\), ou seja, podemos simular \[ X_i \, | \, x_1,x_2,\cdots,x_{i-1},x_{i+1},\cdots,x_p \sim f_i(x_i \, | \, x_1,x_2,\cdots,x_{i-1},x_{i+1},\cdots,x_p), \] para \(i = 1,2,\cdots,p\). O algoritmo de amostragem de Gibbs associado ou amostrador de Gibbs é dado pela seguinte transição de \(X^{(t)}\) para \(X^{(t+1)}\):


Algoritmo 8: Amostrador Gibbs de vários estágios
Na iteração \(t=1,2,\cdots\), dado \(x^{(t)}=(x_1^{(t)},\cdots,x_p^{(t)})\), gerar
1. \(X_1^{(t+1)} \sim f_1(x_1|x_2^{(t)},\cdots,x_p^{(t)})\);
2. \(X_2^{(t+1)} \sim f_2(x_2|x_1^{(t+1)},x_3^{(t)},\cdots,x_p^{(t)})\);
\(\vdots\)
p. \(X_p^{(t+1)} \sim f_p(x_p|x_1^{(t+1)},\cdots,x_{p-1}^{(t+1)})\).


As densidades \(f_1,\cdots,f_p\) são chamados de condicionais completos e uma característica particular do amostrador de Gibbs é que essas são as únicas densidades usadas para simulação. Assim, mesmo em um problema de alta dimensão, todas as simulações podem ser univariadas, o que geralmente é uma vantagem.


Exemplo 7.4

Como uma extensão do Exemplo 7.1, considere a densidade normal multivariada \[\begin{equation} \tag{7.5} (X_1,X_2,\cdots,X_p) \sim N_p\big(0, (1-\rho)\pmb{I}+\rho\pmb{J}\big), \end{equation}\] onde \(\pmb{I}\) é a matriz \(p\times p\) identidade e \(\pmb{J}\) é uma matriz \(p\times p\) de uns.

Este é um modelo para equicorrelação, como \(\mbox{Corr}(X_i,X_j) = \rho\) para cada \(i\) e \(j\). Usando fórmulas padrão para as distribuições condicionais de uma variável aleatória normal multivariada (ver, por exemplo, Johnson and Wichern, 1988) são diretas, mas tedioso, verificar que \[ X_i \, | \, x_{(-i)} \sim N\left( \dfrac{(p-1)\rho}{1+(p-2)\rho}\overline{x}_{(-i)}, \dfrac{1+(p-2)\rho-(p-1)\rho^2}{1+(p-2)\rho}\right), \] onde \(x_{(-i)}=(x_1,x_2,\cdots,x_{i-1},x_{i+1},\cdots,x_p)\) e \(\overline{x}_{(-i)}\) é a média desse vetor.

O amostrador de Gibbs que gera a partir dessas normais univariadas pode então ser facilmente derivado, embora seja inútil para esse problema (Exercício 7.5). É, no entanto, um pequeno passo para considerar a configuração onde as componentes do vetor normal são restritas a um subconjunto de \(\mathbb{R}^p\). Se este subconjunto for um hipercubo, \[ \amalg = \prod_{i=1}^p (a_i,b_i) \] então as condicionais correspondentes simplesmente são as normais acima restritas a \((a_i,b_i)\) para \(i = 1,\cdots,p\). Para restrições mais complexas, um amostrador de Gibbs é quase necessário, pois não existem soluções exatas.

Este amostrador de Gibbs ainda é baseado em condicionais completas normais, que agora estão restritas a subconjuntos da linha real e, portanto, facilmente simuladas como no Exercício 2.22.


Modelos mais complexos que o do Exemplo 7.3 podem ser considerados para o modelo de amostragem normal, como no caso a seguir.


Exemplo 7.5

Uma especificação hierárquica para o modelo normal é o modelo de efeitos aleatórios unidirecionais. Existem diferentes formas de parametrizar este modelo, mas uma possibilidade é a seguinte (veja outras no Exemplo 7.14 e no Exercício 7.24): \[\begin{equation} \tag{7.6} \begin{array}{rcl} X_{ij} & \sim & N(\theta_i,\sigma^2), \qquad i=1,\cdots,k, \quad j=1,\cdots,n_i, \\ \theta_i & \sim & N(\mu,\tau^2), \qquad i=1,\cdots,k, \\ \mu & \sim & N(\mu_0,\sigma^2_\mu), \\ \sigma^2 & \sim & \mathcal{IG}(\alpha_1,\beta_1), \quad \tau \, \sim \, \mathcal{IG}(\alpha_2,\beta_2), \quad \sigma^2_\mu \, \sim \, \mathcal{IG}(\alpha_3,\beta_3)\cdot \end{array} \end{equation}\]

Agora, procedendo como antes e anotarmos a distribuição conjunta dessa hierarquia, podemos derivar o conjunto de condicionais completos: \[\begin{equation} \tag{7.7} \begin{array}{rcl} \theta_i & \sim & N\left(\dfrac{\sigma^2}{\sigma^2+n_i \tau^2}\mu+\dfrac{n_i \tau^2}{\sigma^2+n_i \tau^2}\overline{X}_i, \dfrac{\sigma^2\tau^2}{\sigma^2+n_i \tau^2}\right), \quad i=1,\cdots,k, \\ \mu & \sim & N\left(\dfrac{\tau^2}{\tau^2+k \sigma^2_\mu}\mu_0+\dfrac{k \sigma^2_\mu}{\tau^2+k \sigma^2_\mu}\overline{\theta}, \dfrac{\sigma^2_\mu\tau^2}{\tau^2+k \sigma^2_\mu}\right),\\ \sigma^2 & \sim & \displaystyle \mathcal{IG}\left(\dfrac{n}{2}+\alpha_1, \dfrac{1}{2}\sum_{ij} \big( X_{ij}-\theta_i\big)^2+b_1\right), \\ \tau^2 & \sim & \displaystyle \mathcal{IG}\left(\dfrac{k}{2}+\alpha_2, \dfrac{1}{2}\sum_{i} \big( \theta_{i}-\mu\big)^2+b_2\right), \\ \sigma^2_\mu & \sim & \displaystyle \mathcal{IG}\left(\dfrac{1}{2}+\alpha_3, \dfrac{1}{2}\big( \mu-\mu_0\big)^2+b_3\right), \end{array} \end{equation}\] onde \(n=\sum_i n_i\) e \(\overline{\theta}=\sum_i n_i \theta_i/n\).

Expandindo o estudo do Exemplo 7.3, o conjunto de dados Energy contém dados sobre a ingestão de energia de meninas e meninos, pacote R mcsm. O modelo (7.6) aplica-se, com \(k = 2\) à análise simultânea dos consumos energéticos de moças e rapazes. O resultado do amostrador de Gibbs baseado nas condicionais em (7.7) está resumido na Figura 7.3. 

library(mcsm)
data(Energy)
head(Energy)
##   Girls Boys
## 1    91  457
## 2   504  645
## 3   557  736
## 4   609  790
## 5   693  899
## 6   727  991
nsim = 5000
a = 10
b = 30
x1 = log(Energy$Girls)
x2 = log(Energy$Boys)
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1))
hist(x1, col = "brown", xlab = "Girls", main = "", freq = F)
box(); grid()
hist(x2, col = "orange", xlab = "Boys", main = "", freq = F)
box(); grid()

n1 = length(x1)
n2 = length(x2)
n = n1 + n2
x1bar = mean(x1)
x2bar = mean(x2)
xbar = (n1 * x1bar + n2 * x2bar)/n
k = 2
a1 = a2 = a3 = a
b1 = b2 = b3 = b
mu = rep(xbar, nsim)
theta1 = rep(x1bar, nsim)
theta2 = rep(x2bar, nsim)
sigma2mu = sigma2 = tau2 = rep((var(x1) + var(x2))/2, nsim)
mu0 = xbar
for (i in 2:nsim) {
    B1 = sigma2[i - 1]/(sigma2[i - 1] + n1 * tau2[i - 1])
    theta1[i] = rnorm(1, mean = B1 * mu[i - 1] + (1 - B1) * 
        x1bar, sd = sqrt(tau2[i - 1] * B1))
    B2 = sigma2[i - 1]/(sigma2[i - 1] + n2 * tau2[i - 1])
    theta2[i] = rnorm(1, mean = B2 * mu[i - 1] + (1 - B2) * 
        x2bar, sd = sqrt(tau2[i - 1] * B2))
    B = tau2[i - 1]/(k * sigma2mu[i - 1] + tau2[i - 1])
    m = B * mu0 + (1 - B) * (n1 * theta1[i] + n2 * theta2[i])/n
    mu[i] = rnorm(1, mean = m, sd = sqrt(sigma2mu[i - 1] * 
        B))
    sh1 = (n/2) + a1
    ra1 = (1/2) * (sum((x1 - theta1[i])^2) + sum((x2 - theta2[i])^2)) + b1
    sigma2[i] = 1/rgamma(1, shape = sh1, rate = ra1)
    sh2 = (k/2) + a2
    ra2 = (1/2) * (sum((theta1[i] - mu[i])^2) + sum((theta2[i] - mu[i])^2)) + b2
    tau2[i] = 1/rgamma(1, shape = sh2, rate = ra2)
    sh3 = (1/2) + a3
    ra3 = (1/2) * (mu[i] - mu0)^2 + b3
    sigma2mu[i] = 1/rgamma(1, shape = sh3, rate = ra3)
}
par(mfrow = c(2, 3), mar = c(4, 4, 1, 1))
top = max(mu)
bot = min(mu)
hist(mu[(nsim/10):nsim], col = "grey", breaks = 25, xlab = "", 
    main = expression(mu), freq = F)
box(); grid()
hist(theta1[(nsim/10):nsim], col = "grey", breaks = 25, xlab = "", 
    main = expression(theta[1]), freq = F)
box(); grid()
hist(theta2[(nsim/10):nsim], col = "grey", breaks = 25, xlab = "", 
    main = expression(theta[2]), freq = F)
box(); grid()
hist(sqrt(sigma2mu[(nsim/10):nsim]), col = "sienna", breaks = 50, 
    xlim = c(0.5, 4), xlab = "", main = expression(sigma[mu]), freq = F)
box(); grid()
hist(sqrt(tau2[(nsim/10):nsim]), col = "sienna", breaks = 25, 
    xlim = c(0.5, 4), xlab = "", main = expression(tau), freq = F)
box(); grid()
hist(sqrt(sigma2[(nsim/10):nsim]), col = "sienna", breaks = 25, 
    xlim = c(0.5, 2), xlab = "", main = expression(sigma), freq = F)
box(); grid()

Figura. 7.3. Histogramas das distribuições a posteriori marginais do amostrador de Gibbs com base em 5000 iterações. A linha superior apresenta histogramas para a média subjacente \(\mu\) e as médias, \(\theta_1\) e \(\theta_2\), para a energia das meninas e dos meninos. A linha inferior corresponde aos desvios padrão.



7.4 Dados ausentes e variáveis latentes

Começando com o amostrador de Gibbs de dois estágios, trabalhando em uma distribuição conjunta \(f(x,y)\), parece haver uma grande diferença com o algoritmo Metropolis-Hastings que trabalha com uma única distribuição ou, em outras palavras, gera todos os componentes de \((x,y)\) de uma só vez.

Essa diferença no alvo é ilusória porque, uma vez dado \(f(x,y)\), podemos usar o amostrador de Gibbs relevante ou um algoritmo genérico de Metropolis-Hastings, enquanto, se for dada uma densidade marginal \(f_X(x)\), podemos construir \(f_X(x)\) de uma densidade conjunta \(f(x,y)\) para ajudar na simulação, onde a segunda variável \(Y\) é então uma variável auxiliar que pode não ser diretamente relevante do ponto de vista estatístico.

Existem muitas configurações em que existe uma conclusão natural de \(f_X(x)\) em \(f(x,y)\) e, de fato, isso pode levar a um amostrador de Gibbs eficaz. Obviamente, é sempre o caso que qualquer densidade dada \(f_X(x)\) pode ser artificialmente completada em uma densidade conjunta \(f(x,y)\), conforme demonstrado com o amostrador de corte (slice) no final desta seção.

Essas considerações nos trazem de volta ao domínio dos modelos de dados ausentes, conforme descrito na Seção 5.4.2, onde a representação (5.9) \[ g(x | \theta) = \int_\mathcal{Z} f(x,z | \theta)\mbox{d}z, \] foi introduzida.

Conforme discutido no Capítulo 5, \(g(x|\theta)\) é a densidade das observações, ou seja, a verossimilhança e o lado direito representa a densidade conjunta completada. A densidade \(f\) é arbitrária e pode ser escolhida de modo que as funções condicionais completos de \(f\) sejam fáceis de simular e o algoritmo de Gibbs (Algoritmo 8) é implementado em \(f\) em vez de \(g\) e a condicional completo correspondente de dado \((x,z)\).

Dependendo do campo de aplicação, tais representações recebem nomes diferentes. De uma perspectiva matemática, (5.9) é um modelo de mistura. Em estatística, usamos com mais frequência o nome de modelos de dados perdidos, enquanto os econometristas preferem o uso de modelos de variáveis latentes.

Se fatorarmos \(f(x,z|\theta) = f(x|z,\theta)h(z|\theta)\), então (5.9) se torna \[ g(x|\theta) =\int_\mathcal{Z} f(x|z,\theta)h(z|\theta) \mbox{d}z, \] e \(h(z|\theta)\), a distribuição marginal dos dados ausentes \(z\), é claramente uma distribuição mista.

Em uma configuração geral de dados ausentes, \[ g(z) = \int_\mathcal{Z} f(x,z)\mbox{d}z, \] para \(p\geq 2\), escrevemos \(y = (x,z) = (y_1,\cdots,y_p)\) e denotamos as densidades condicionais de \[ f(y) = f(y_1,\cdots,y_p) \] por \[\begin{equation} \begin{array}{rcl} Y_1 | y_2,\cdots,y_p & \sim & f_1(y_1|y_2,\cdots,y_p), \\ Y_2 | y_1,y_3\cdots,y_p & \sim & f_2(y_2|y_1,y_3,\cdots,y_p), \\ \vdots & & \vdots \\ Y_p | y_1,\cdots,y_{p-1} & \sim & f_p(y_p|y_1,\cdots,y_{p-1})\cdot \end{array} \end{equation}\]

Então, aplicar um amostrador de Gibbs de múltiplos estágios como no Algoritmo 8 a esses condicionais completos e assumir que todos podem ser simulados leva a uma Cadeia de Markov \(\{Y^{(t)}\}_t\) que converge para \(f\) e, portanto, a uma subcadeia \(\{X^{(t)}\}_t\) que converge para \(g\).


Exemplo 7.6

Nos Exemplos 5.13 e 5.14, tratamos um modelo de dados censurados como um modelo de dados perdidos. Identificamos \(g(x|\theta)\) com a função de verossimilhança \[ g(x|\theta) = L(\theta|x) \approx \prod_{i=1}^m \exp\big( -(x_i-\theta)^2/2\big) \] e \[ f(x,z|\theta) = L(\theta|x,z)\approx \prod_{i=1}^m \exp\big( -(x_i-\theta)^2/2\big)\prod_{i=m+1}^n \exp\big( -(z-i-\theta)^2/2\big) \] é a verossimilhança dos dados completo.

Dada uma distribuição a priori \(\pi(\theta)\) em \(\theta\), podemos criar um amostrador de Gibbs que itera entre as distribuições condicionais \[ \pi(\theta|x,z) \qquad \mbox{e} \qquad f(z|x,\theta) \] e terá distribuição estacionária \(\pi(\theta,z|x)\), distribuição a posteriori de \((\theta,z)\).

Tomando como distribuição a prior \(\pi(\theta) = 1\), a distribuição condicional de \(\theta|x,z\) é dada por \[ \theta \, | \, x,z\sim N\left( \dfrac{m\overline{x}+(n-m)\overline{z})}{n},\dfrac{1}{n}\right), \] enquanto a distribuição condicional de \(Z|x,\theta\) é o produto das normais truncadas \[ Z_i \, | \, x,\theta \sim \varphi(z-\theta)\Big/ \Big( 1-\Phi(a-\theta)\Big), \] como cada \(Z_i\) deve ser maior que o ponto de truncamento \(a\).

A geração de valores de \(Z\) pode ser feita através da função R rtrun do pacote bayesm, ver Exercícios 7.21 e 7.7. O resultado do amostrador de Gibbs, cujo núcleo R pode ser escrito como a seguir é resumido na Figura 7.4, usando as distribuições a posteriori de \(\theta\) e \(Z\). 

library(bayesm)
set.seed(546)
Nsim=5000
that=zbar=xbar=rep(0,Nsim) #init arrays
a=3.5
for(i in 2:Nsim){
  zbar[i]=mean(rtrun(mu=rep(that[i-1],n-m),sigma=rep(1,n-m),
                     a=rep(a,n-m),b=rep(Inf,n-m)))
  that[i]=rnorm(1,(m/n)*xbar+(1-m/n)*zbar[i],sqrt(1/n))
}
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1))
hist(that, col = "grey", breaks = 50, xlab = "", 
     main = expression(theta), xlim = c(4,7),freq = FALSE)
box(); grid()
hist(zbar, col = "sienna", breaks = 50, xlab = "", 
     main = expression(bar(Z)), xlim = c(4,7), freq = FALSE)
box(); grid()

Figura. 7.4. Histogramas das distribuições a posteriori de \(\theta\) e \(\overline{Z}\). O ponto de truncamento no \(Z_i\) é \(a = 3.5\).



Exemplo 7.7

Lembremos do modelo multinomial do Exemplo 5.16, \[ Multinomial\left(n;\dfrac{1}{2}+\dfrac{\theta}{4};\dfrac{1}{4}(1-\theta),\dfrac{1}{4}(1-\theta),\dfrac{^\theta}{4} \right), \] onde estimamos \(\theta\) usando etapas de EM ou MCEM, introduzindo a variável latente \(Z\) com a desmarginalização \[ \big(z,x_1-z,x_2,x_3,x_4\big) \sim Multinomial\left(n;\dfrac{1}{2}+\dfrac{\theta}{4};\dfrac{1}{4}(1-\theta),\dfrac{1}{4}(1-\theta),\dfrac{^\theta}{4} \right)\cdot \]

Se usarmos um uniforme antes de \(\theta\), os condicionais completos podem ser recuperados como \[ \theta \sim Beta(z+x_4+1, x_2+x_3+1) \qquad \mbox{e} \qquad z\sim Binomial\left(x_1,\dfrac{\theta}{2+\theta} \right), \] levando ao amostrador Gibbs cuja saída é resumida na Figura 7.5. 

Nsim = 10^4
x=c(125,18,20,34)    #data
theta=z=rep(.5,Nsim) #init chain
for (j in 2:Nsim){
  theta[j]=rbeta(1,z[j-1]+x[4]+1,x[2]+x[3]+1)
  z[j]=rbinom(1,x[1],(theta[j]/(2+theta[j])))
}
par(mfrow=c(1,2),mar=c(4,4,1,1))
hist(theta, col = "grey", breaks = 25, xlab = "", main = expression(theta), freq = F)
box(); grid()
hist(z, col = "sienna", breaks = 25, xlab = "", main = expression(z), freq = F)
box(); grid()

Figura. 7.5. Histogramas de distribuições marginais do amostrador de Gibbs do Exemplo 7.7. O interesse principal está na distribuição marginal de \(\theta\).


Este exemplo mostra um caso em que tanto o amostrador EM quanto o Gibbs se aplicam. Como sempre, a abordagem Bayesiana permite uma inferência mais completa que inclui intervalos de confiança.


Exemplo 7.8

Uma generalização do modelo do Exemplo 7.7 é o modelo \[\begin{equation} \tag{7.8} X\sim Multinomial_5 \Big(n; a_1\theta_1 + b_1, a_2\theta_1 + b_2, a_3\theta_2 + b_3, a_4\theta_2 + b_4, c(1-\theta_1-\theta_2)\Big), \end{equation}\] com \(0\leq a_1 + a_2 = a_3 + a_4 = 1-\sum_{i=1}^4 b_i = c\leq 1\), onde os \(a_i,b_i\geq 0\) são conhecidos com base em considerações genéticas, como na Tabela 7.1 que descreve as probabilidades dos quatro tipos de sangue como funções das probabilidades do genótipo, por causa do alelo domínante.

Nosso interesse é estimar as frequências alélicas \(p_{_A}\), \(p_{_B}\) e \(p_{_O}\), que somam 1. Podemos então aumentar os dados com \(Z = (Z_1,Z_2,Z_3,Z_4)\) como \[ X_1 = Z_1 + Z_2, \qquad X_2 = Z_3 + Z_4, \qquad X_3 = Z_5 + Z_6, \qquad X_4 = Z_7 + Z_8, \] que desmarginaliza o modelo para nos permitir amostrar de \[ Y_9 \sim Multinomial_9\Big(n;a_1 \theta_1, b_1, a_2\theta_1, b_2, a_3\theta_3, b_3,a_4\theta_2, c(1-\theta_1-\theta_2) \Big), \] com \(Y = (Z_1, X_1-Z_1,Z_2,X_2-Z_2,Z_3,X_3-Z_3,Z_4,X_4-Z_4,X_5)\). Ver o Exercício 7.23 para uma solução alternativa.

Uma a priori natural em \((\theta_1,\theta_2)\) é a distribuição \(Dirichlet(\alpha_1,\alpha_2,\alpha_3)\), \[ (\theta_1,\theta_2) \approx \theta_1^{\alpha_1-1}\theta_2^{\alpha_2-1} (1-\theta_1-\theta_2)^{\alpha_3-1}, \] o que leva aos condicionais completos \[\begin{equation} \tag{7.9} \begin{array}{rcl} \big(\theta_1,\theta_2,1-\theta_1-\theta_2 \big) \, | \, x,z & \sim & Dirichlet\big(z_1+z_2+\alpha_1,z_3+z_4+\alpha_2,x_5+\alpha_3 \big), \\ Z_i \, | \, x,\theta_1,\theta_2 & \sim & Beta\left( x_1, \dfrac{a_i \theta_1}{a_i \theta_1+b_i}\right), \qquad i=1,3, \\ Z_i \, | \, x,\theta_1,\theta_2 & \sim & Beta\left( x_1, \dfrac{a_i \theta_2}{a_i \theta_2+b_i}\right), \qquad i=5,7, \end{array} \end{equation}\] que podem ser facilmente simulados e, portanto, incluídos em um amostrador Gibbs. A Figura 7.6 mostra as distribuições de cadeias produzidas por tal amostrador.

\[ \begin{array}{c|c|c|c} \mbox{Genotipo} & \mbox{Probabilidade} & \mbox{Observado} & \mbox{Probabilidade} & \mbox{Frequência} \\\hline \mbox{AA} & p_{_A}^2 & A & p_{_A}^2+2p_{_A}p_{_O} & n_{_A}=186 \\ \mbox{AO} & 2p_{_A}p_{_O} & & & \\ \mbox{BB} & p_{_B}^2 & B & p_{_B}^2+2p_{_B}p_{_O} & n_{_B}=38 \\ \mbox{BO} & 2p_{_B}p_{_O} & & & \\ \mbox{AB} & 2p_{_A}p_{_B} & AB & p_{_A}p_{_B} & n_{_{AB}}=13 \\ \mbox{OO} & p_{_O}^2 & O & p_{_O}^2 & n_{_O}=284 \\\hline \end{array} \] Tabela 7.1. Frequências genotípicas observadas em dados de tipo sanguíneo. O efeito de um alelo dominante cria um problema de falta de dados.

set.seed(321)
# Este é um gerador aleatório para distribuições de Dirichlet, baseado na normalização 
# de variáveis aleatórias Gama 
rdirichlet=function (n, shape) 
  {
    d = length(shape)
    x = rgamma(n * d, rep(shape, n))
    x = matrix(x, nrow = d)
    t(sweep(x, 2, apply(x, 2, sum), "/"))
}
nsim=5000
nA=186;nB=38;nAB=13;nO=284          #Blood type data
pA=rep(.25,nsim);pB=array(.05,nsim)
for (i in 2:nsim){
  ZA=rbinom(1,nA,pA[i-1]^2/(pA[i-1]^2+2*pA[i-1]*(1-pA[i-1]-pB[i-1])))
  ZB=rbinom(1,nB,pB[i-1]^2/(pB[i-1]^2+2*pB[i-1]*(1-pA[i-1]-pB[i-1])))
  temp=rdirichlet(1,c(nA+nAB+ZA+1,nB+nAB+ZB+1,nA-ZA+nB-ZB+2*nO+1))
  pA[i]=temp[1];pB[i]=temp[2]
  }
par(mfrow=c(1,3),mar=c(4,4,2,1))
hist(pA,xlab=expression(p[A]),freq=F,col="grey",breaks=25,main="")
box();grid()
hist(pB,xlab=expression(p[B]),freq=F,col="sienna",breaks=25,main="")
box();grid()
hist(1-pA-pB,xlab=expression(p[O]),freq=F,col="gold",breaks=25,main="")
box();grid()

Figura. 7.6. Histogramas das distribuições marginais das probabilidades dos genótipos do amostrador de Gibbs do Exemplo 7.8.


Modelos de mistura finita, que vimos com algum detalhe no Capítulo 5, Exemplo 5.2 e Capítulo 6, Exemplo 6.5, são obviamente candidatos à desmarginalização por meio de variáveis latentes; isto é, como um caso especial de mistura. Conforme já descrito no Exemplo 5.12, dada uma amostra \((x_1,\cdots,x_n)\) de uma distribuição de mistura \[ \sum_{j=1}^k p_j f(x \, | \, \xi_j) \] onde \(\sum_j p_j=1\) e \(f(\cdot|\xi_j)\) é uma densidade parametrizada com parâmetro desconhecido \(\xi_j\), podemos associar a cada observação \(x_i\) uma variável latente \(z_i\in \{1,\cdots,k\}\) que indica qual componente da mistura está associado a \(x_i\).

A conclusão correspondente do modelo de mistura acima é então \[ Z_i \sim Multinomial_k \big(1;p_1,\cdots,p_k \big), \qquad x_i|z_i \sim f(x|\xi_{z_i})\cdot \]

Assim, considerar \(y_i = (x_i, z_i)\), em vez de \(x_i\), elimina totalmente a estrutura da mistura, pois a verossimilhança do modelo completo é \[ \ell(p,\xi \, | \, y_1,\cdots,y_n) \approx \prod_{i=1}^n p_{z_i} f(x_i \, | \, \xi_{z_i}) = \prod_{j=1}^k \prod_{\{i \,: \, z_i=j\}} p_j f(x_i \, | \, \xi_j)\cdot \]

Pode-se perguntar por que a conclusão é útil nesse cenário, uma vez que a verossimilhança observada pode ser calculada de forma fechada, como mostrado, por exemplo, na Figura 5.2, que representa uma verossimilhança de mistura em uma grade de pixels como no Exemplo 6.5, onde produzimos um algoritmo Metropolis-Hastings de caminhada aleatória.

Como no algoritmo EM dos Exemplos 5.12 e 5.13, o uso das variáveis indicadoras latentes produz um algoritmo de simulação geralmente eficiente que se concentra rapidamente na(s) moda(s) da distribuição a posteriori.

As duas etapas do amostrador de Gibbs são então associadas às a posteriori condicionais completas \[ P(Z_i=j \, | \, x,\xi) \approx p_j f(x_i \, | \, \xi_j), \qquad i=1,\cdots,n, \quad j=1,\cdots,k \] e \[ \begin{array}{rcl} \xi_j \, | \, y & \sim & \pi\left(\xi \Big| \, \dfrac{\lambda_j\alpha_j + n_j\overline{x}_j}{\lambda_j+n_j}, \lambda_j+n_j \right),\\ p & \sim & Dirichlet_k\left( \gamma_1+n_1,\cdots,\gamma_k+n_k\right), \end{array} \] onde \[ n_j=\sum_{i=1}^n \pmb{1}_{\{z_i=j\}}, \qquad n_j\overline{x}_j = \sum_{i=1}^n x_i \pmb{1}_{\{z_i=j\}}\cdot \]

Neste amostrador de Gibbs de dois passos, a geração a partir da a posteriori associada à verossimilhança completa não é detalhada, pois ela irá variar dependendo do modelo de amostragem e da a priori utilizada. Na situação padrão que depende de uma família exponencial para \(f(\cdot|\xi)\) e um conjugado a priori de \(\xi\), essa geração é obviamente direta.


Exemplo 7.9

Como ilustração, considere a mesma configuração do Exemplo 5.12, ou seja, uma mistura normal com dois componentes com igual variância conhecida e pesos fixos, \[ p N(\mu_1,\sigma^2) + (1-p)N(\mu_2,\sigma^2)\cdot \] Além disso, assumimos uma distribuição a priori normal \(N(0, \nu^2\sigma^2)\), com \(\nu^2\) conhecido, em ambas as médias \(\mu_1\) e \(\mu_2\).

As variáveis latentes \(z_i\) são as mesmas do Exemplo 5.12, ou seja, \[ P(Z_i=1)=1-P(Z_i=2)=p \qquad \mbox{e} \qquad X_1 \, | \, Z_i=k \sim N(\mu_k,\sigma^2)\cdot \]

A distribuição completa é então \[ \begin{array}{l} \pi(\mu_1,\mu_2,z \, | \, x) \approx \exp\big( -(\mu_1^2+\mu_2)^2/\nu^2\sigma^2\big) \\ \qquad \qquad \qquad \displaystyle \times p\prod_{\{i \, : \, z_i=1\}} \exp\left(-(x_i-\mu_1)^2/2\sigma^2 \right) \times (1-p) \prod_{\{i \, : \, z_i=1\}} \exp\left(-(x_i-\mu_2)^2/2\sigma^2 \right), \end{array} \] para o qual os condicionais completos dos \(\mu_j\)’s são facilmente derivados (Exercício 7.10). 

set.seed(7452)
Niter = 10^4
v = 1
da = sample(c(rnorm(10^2), 2.5 + rnorm(4 * 10^2)))
like = function(mu) sum(log((0.2 * dnorm(da - mu[1]) + 0.8 * dnorm(da - mu[2]))))
mu1 = mu2 = seq(-2, 5, le = 250)
lli = matrix(0, nco = 250, nro = 250)
for (i in 1:250) for (j in 1:250) lli[i, j] = like(c(mu1[i], mu2[j]))
x = prop = runif(2, -2, 5)
the = matrix(x, ncol = 2)
curlike = hval = like(x)
for (i in 2:Niter) {
  pp = 1/(1 + ((0.8 * dnorm(da, mean = the[i - 1, 2]))/(0.2 * 
            dnorm(da, mean = the[i - 1, 1]))))
  z = 2 - (runif(length(da)) < pp)
  prop[1] = (v * sum(da[z == 1]))/(sum(z == 1) * v + 1) + 
    rnorm(1) * sqrt(v/(1 + sum(z == 1) * v))
  prop[2] = (v * sum(da[z == 2]))/(sum(z == 2) * v + 1) + 
    rnorm(1) * sqrt(v/(1 + sum(z == 2) * v))
  curlike = like(prop)
  hval = c(hval, curlike)
  the = rbind(the, prop)
}
par(mar = c(4, 4, 1, 1))
image(mu1, mu2, lli, xlab = expression(mu[1]), ylab = expression(mu[2]))
contour(mu1, mu2, -lli, nle = 100, add = TRUE)
points(the[, 1], the[, 2], cex = 0.6, pch = 19)
lines(the[, 1], the[, 2], cex = 0.6, pch = 19)

Figura. 7.7. Amostrador de Gibbs de 5000 pontos para a mistura a posteriori contra a superfície do logaritmo da a posteriori.

A Figura 7.7 ilustra o comportamento do amostrador de Gibbs correspondente usando um conjunto de dados simulados \(x\) de 500 pontos da distribuição \(0.7N(0,1)+0.3N(2.7,1)\). Esta figura mostra a amostra MCMC após 15000 iterações.

Esta simulação está de fato em claro acordo com a superfície a posteriori. Embora possa parecer muito concentrado em torno de uma moda, você deve levar em conta o fato de que a segunda moda representada neste gráfico é muito menor, pois há uma diferença de pelo menos 50 nos valores a posteriori.


Como um último exemplo de um amostrador de Gibbs variável latente, olhamos para o amostrador de fatias (slice sampler), que aparece mais como um tipo genérico de desmarginalização. Ver Neal (2003), para um tratamento abrangente. De fato, podemos alternativamente considerar o amostrador de Gibbs como sendo derivado do amostrador de fatia; ver Robert and Casella (2004).

Dada uma densidade de interesse \(f_{_X}(x)\), podemos sempre representá-la como a densidade marginal da densidade da conjunta \[ f(x,u)=\pmb{1}_{\{0< u<f_{_X}(x)\}}, \] já que integrar a função acima em \(u\) retorna \(f_{_X}\). As densidades condicionais associadas são \[ f_{_{X|U}}(x|u)=\dfrac{\pmb{1}_{\{0< u<f_{_X}(x)\}}}{\displaystyle \int \pmb{1}_{\{0< u<f_{_X}(x)\}} \mbox{d}x}, \qquad f_{_{U|X}}(u|x)=\dfrac{\pmb{1}_{\{0< u<f_{_X}(x)\}}}{\displaystyle \int \pmb{1}_{\{0< u<f_{_X}(x)\}} \mbox{d}u}, \] o que significa que ambos são uniformes. Esses dois condicionais então definem o amostrador de fatia como o amostrador de Gibbs associado.


Algoritmo 9: Amostrador 2D amostrador de fatias (slice sampler)
Na iteraçõa \(t\), simular
1. \(U^{(t+1)}\sim U_{[0,f(x^{(t)})]}\);
2. \(X^{(t+1)}\sim U_{A^{(t+1)}}\), com
\[ A^{(t+1)} = \{x \, : \, f(x)\geq u^{(t+1)}\}\cdot \]


O apelo desse algoritmo é que ele se aplica formalmente a qualquer densidade conhecida até uma constante multiplicativa sem restrição em sua forma ou dimensão. Obviamente, sua implementação pode ser dificultada pela simulação uniforme sobre o conjunto \(A^{(t)}\).


Exemplo 7.10

Considere a densidade \(f(x)=\frac{1}{2}e^{-\sqrt{x}}\) definida para \(x>0\). Embora possa ser simulada diretamente, ela também cede facilmente ao amostrador de fatias. De fato, aplicando as fórmulas acima, temos \[ U \, | \, x \sim U\left(0,\frac{1}{2}e^{-\sqrt{x}}\right), \qquad X \, | \, u \sim U\big(0,\log^2(2u)\big)\cdot \]

Implementamos o amostrador para gerar 5.000 variáveis e mostrá-las junto com a densidade na Figura 7.8, o que mostra que a concordância é muito boa. O painel direito mostra algumas autocorrelações fortes, o que é típico do amostrador de fatias. 

fff=function(y)exp(-sqrt(y))/2;
nsim=10^3;x=rep(pi,nsim);x[1]=-log(runif(1))
for (i in 2:nsim){
  w=runif(1,min=0,max=fff(x[i-1]));x[i]=runif(1,min=0,max=(-log(2*w))^2)
}
par(mfrow=c(1,2),mar=c(4,4,2,1))
hist(x,xlab="Slice Sampler",main="",xlim=c(0,40),ylim=c(0,.25),
     freq=FALSE,col="grey",breaks=250)
box();grid()
curve(fff,add=TRUE,lwd=3,col="sienna",xlim=c(0,40))
acf(x,xlab="Autocorrelation",ylab="",lwd=3)
grid()

Figura. 7.8. Um histograma de amostrador de fatias e densidade alvo para o Exemplo 7.10 usando 5.000 iterações. O painel esquerdo é o histograma com a densidade real sobreposta e o painel direito mostra a função de autocorrelação.


Há uma extensão óbvia para o amostrador de fatias 2D acima, semelhante à extensão de vários estágios para o amostrador Gibbs de dois estágios.

Se a densidade alvo for escrita como um produto de funções, \[ f(x)=\prod_{i=1}^n g_i(x), \] como, por exemplo, no caso de uma distribuição a posteriori associada a uma amostra de \(n\) observações, onde os \(g_i\) são então as densidades componentes, uma conclusão associada é \[ f(x,u_1,\cdots,u_n)=\prod_{i=1}^n \pmb{1}_{\{0<u_i<g_i(x)\}}, \] o que leva a um amostrador de fatias com \((n + 1)\) etapas, \(X^{(t)}\) sendo gerado uniformemente sobre o conjunto \[ A^{(t)}=\bigcap_{i=1}^n \left\{ x \, : \, g_i(x)>u_i^{(t)} \right\}\cdot \]


Exemplo 7.11

Lembremos da regressão logística, que vimos pela primeira vez no Exemplo 4.11 e ajustada com um algoritmo de Metropolis-Hastings no Exercício 6.13.

O modelo é \[ Y_i \sim Bernoulli\big( p(x_i)\big), \qquad p(x) = \dfrac{\exp(\alpha+\beta x)}{1+\exp(\alpha+\beta x)}, \] onde \(p(x)\) é a probabilidade de sucesso e \(x\) é uma covariável unidimensional. A verossimilhança associada com a amostra \((y,x)=(y_1,x_1),\cdots,(y_n,x_n)\) é \[ L(\alpha,\beta \, | \, y) \approx \prod_{i=1}^n \left( \dfrac{\exp(\alpha+\beta x_i)}{1+\exp(\alpha+\beta x_i)}\right)^{y_i} \left( \dfrac{1}{1+\exp(\alpha+\beta x_i)}\right)^{1-y_i}\cdot \]

Usando uma a prior não informativa em \((a,b)\), a distribuição a posteriori pode ser associada a um amostrador de fatias (slice sampler) baseado nas variáveis uniformes \[ U-i \sim U\left(0,\dfrac{\exp\big(y_i(\alpha+\beta x_i)\big)}{1+\exp(\alpha+\beta x_i)} \right)\cdot \]

Gerando uma distribuição uniforme sobre o conjunto \[ \left\{ (a,b) \, : \, y_i(a+be x_i)>\log\big(u_i/(1-u_i)\big)\right\} \] sendo bastante pesado, podemos decompor ainda mais a simulação uniforme simulando consecutivamente \[ a^{(t)} \sim U\left( \max_{\{i\, : \, y_i=1\}} \Big(\log\big(u_i^{(t)}/(1-u_i^{(t)})\big)-b^{(t-1)}x_i\Big), \min_{\{i\, : \, y_i=0\}}\Big( \log\big(u_i^{(t)}/(1-u_i^{(t)})\big)-b^{(t-1)}x_i\Big)\right) \] e \[ b^{(t)} \sim U\left( \max_{\{i\, : \, y_i=1\}} \left(\log\left(\dfrac{u_i^{(t)}}{1-u_i^{(t)}}-a^{(t)}\right) \Big/ x_i\right), \min_{\{i\, : \, y_i=0\}} \left(\log\left(\dfrac{u_i^{(t)}}{1-u_i^{(t)}}-b^{(t-1)}\right)\Big/x_i \right) \right), \] se assumirmos, sem perda de generalidade, que todos os \(x_i\) são positivos.

No entanto, executar o amostrador de fatia correspondente no conjunto de dados challenger descrito no Exercício 6.13 exibe um comportamento de passeio aleatório na cadeia \(\{(a^{(t)},b^{(t)})\}\), conforme mostrado na Figura 7.9. Portanto, introduzimos \(N(0,\sigma^2)\) a prioris normais em \(a\) e \(b\). A modificação no amostrador de fatias é mínima, pois ambas as distribuições uniformes acima são substituídas por normais \(N(0,\sigma^2)\) truncadas, sendo os intervalos de truncamento os usados acima. 

challenger=read.csv2("http://leg.ufpr.br/~lucambio/SSim/challenger.csv", sep=",")
set.seed(15191)
Nsim = 10^4
rtrun = function(mu = 0, sigma = 1, minn = -Inf, maxx = Inf) {
  qnorm(pnorm(minn, mean = mu, sd = sigma) + (pnorm(maxx, 
      mean = mu, sd = sigma) - pnorm(minn, mean = mu, sd = sigma)) * runif(1))
}
x = challenger[, 2]
y = challenger[, 1]
n = length(x)
a = b = rep(0, Nsim)
a[1] = glm(challenger[, 1] ~ challenger[, 2])$coef[1]
b[1] = glm(challenger[, 1] ~ challenger[, 2])$coef[2]
meana = 0
sigmaa = 5
meanb = 0
sigmab = 5/sd(x)
for (t in 2:Nsim) {
  uni = runif(n) * exp(y * (a[t - 1] + b[t - 1] * x))/(1 + exp(a[t - 1] + b[t - 1] * x))
  mina = max(log(uni[y == 1]/(1 - uni[y == 1])) - b[t - 1] * x[y == 1])
  maxa = min(-log(uni[y == 0]/(1 - uni[y == 0])) - b[t - 1] * x[y == 0])
  a[t] = rtrun(0, sigmaa, mina, maxa)
  uni = runif(n) * exp(y * (a[t] + b[t - 1] * x))/(1 + exp(a[t] + b[t - 1] * x))
  minb = max((log(uni[y == 1]/(1 - uni[y == 1])) - a[t])/x[y == 1])
  maxb = min((-log(uni[y == 0]/(1 - uni[y == 0])) - a[t])/x[y == 0])
  b[t] = rtrun(0, sigmab, minb, maxb)
}
plot(a[(Nsim-100):Nsim],b[(Nsim-100):Nsim],type="l",xlab="a",ylab="b")
grid()

Figura. 7.9. Evolução da cadeia \(\{(a^{(t)},b^{(t)})\}\) ao longo de \(10^3\) iterações finais para o conjunto de dados challenger sob um flat prior.



7.5 Estruturas hierárquicas

Vimos o amostrador de Gibbs de vários estágios aplicado a vários exemplos, a maioria decorrente de estruturas de dados ausentes. No entanto, é igualmente adequado para amostrar de maneira direta a partir de qualquer modelo hierárquico.

Um modelo hierárquico é definido por uma sequência de distribuições condicionais como, por exemplo, na hierarquia genérica de dois níveis \[ \begin{array}{rcl} X_i & \sim & f_i(x|\theta), \quad i=1,\cdots,n, \quad \theta=(\theta_1,\cdots,\theta_p), \\ \theta_j & \sim & \pi_j(\theta|\gamma), \quad j=1,\cdots,p, \quad \gamma=(\gamma_1,\cdots,\gamma_s), \\ \gamma_k & \sim & g(\gamma), \quad k=1,\cdots,s\cdot \end{array} \]

A distribuição conjunta desta hierarquia é \[ \prod_{i=1}^n f_i(x_i|\theta) \prod_{j=1}^p \pi_j(\theta_j|\gamma)\prod_{k=1}^s g(\gamma_k)\cdot \]

Assumindo que os \(x_i\)’s são observações, a distribuição a posteriori correspondente de \((\theta,\gamma)\) está associada às condicionais a posteriori completas \[ \begin{array}{rcl} \theta_j & \approx & \displaystyle \pi_j(\theta_j|\gamma)\prod_{i=1}^n f_i(x_i|\theta), \quad j=1,\cdots,p, \\ \gamma_k & \approx & \displaystyle g(\gamma_k)\prod_{j=1}^p \pi_j(\theta_j|\gamma), \quad k=1,\cdots,s\cdot \end{array} \]

Em hierarquias padrão, essas densidades são fáceis de simular e, portanto, são naturalmente associadas a um amostrador de Gibbs. Em hierarquias mais complexas, podemos precisar usar métodos mais sofisticados, como uma etapa de Metropolis-Hastings ou outro amostrador de fatia (slice), para amostrar a partir dos condicionais, conforme explicado na Seção 7.6.3.

No entanto, nossa principal mensagem aqui é que os condicionais completos são bastante fáceis de escrever dada a especificação hierárquica, enquanto reduzem consideravelmente a dimensão das variáveis aleatórias a serem simuladas em cada etapa.

Quando um condicional completo em um amostrador de Gibbs não pode ser simulado diretamente, é suficiente executar, em vez disso, uma única etapa de qualquer algoritmo MCMC associado a esse condicional completo. A validação teórica é a mesma de qualquer amostrador MCMC. Caso seja utilizado um amostrador slice para este fim, a variável auxiliar é simplesmente adicionada ao vetor de parâmetros.


Exemplo 7.12

Um exemplo hierárquico de referência na literatura de amostragem de Gibbs descreve falhas múltiplas de dez bombas em uma usina nuclear, com os dados fornecidos na Tabela 7.2.

\[ \mbox{Tabela 7.2. Número de falhas e tempos de observação de dez bombas em uma usina nuclear.}\\ \mbox{Fonte: Gaver and O'Muircheartaigh (1987).}\\ \begin{array}{ccccccc} \hline \mbox{Pump} & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\\hline \mbox{Failures} & 5 & 1 & 5 & 14 & 3 & 19 & 1 & 1 & 1 & 22 \\ \mbox{Time} & 94.32 & 15.72 & 62.88 & 125.76 & 5.24 & 31.44 & 1.05 & 1.05 & 2.10 & 10.48 \\\hline \end{array} \]

A modelagem é baseada na suposição de que o número de falhas da \(i\)-ésima bomba segue um processo de Poisson com parâmetro \(\lambda_i\), \(1\leq \lambda_i\leq 10\).

Para um tempo de observação \(t_i\), o número de falhas \(X_i\) é, portanto, uma variável aleatória \(Poisson(\lambda_i t_i)\). As distribuições a priori padrão são distribuições gama, que levam ao modelo hierárquico \[ \begin{array}{rcl} X_i & \sim & Poisson(\lambda_i t_i) \quad i=1,\cdots,10, \\ \lambda_i & \sim & Gama(\alpha,\beta), \quad i=1,\cdots,10, \\ \beta & \sim & Gama(\gamma,\delta)\cdot \end{array} \]

A distribuição conjunta é então \[ \pi(\lambda_1,\cdots,\lambda_{10},\beta \,| \, t_1,\cdots,t_{10},p_1,\cdots,p_{10}) \approx \prod_{i=1}^{10} \Big((\lambda_i t_i)^{x_i} e^{-\lambda_i t_i}\lambda_i^{\alpha-1}e^{-\beta \lambda_i} \Big)\beta^{10\alpha}\beta^{\gamma-1}e^{-\delta \beta} \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \approx \prod_{i=1}^{10} \Big(\lambda_i^{x_i+\alpha-1} e^{-(t_1+\beta)\lambda_i} \Big)\beta^{10\alpha+\gamma-1}e^{-\delta \beta}, \] levando às distribuições condicionais completas \[ \lambda_i \, | \, \beta,t_i,x_i \sim Gama(x_i+\alpha,t_i+\beta), \quad i=1,\cdots,10 \\ \beta \, | \, \lambda_1,\cdots,\lambda_{10} \sim Gama\left(\gamma+10\alpha,\delta+\sum_{i=1}^{10}\lambda_i\right)\cdot \]

O amostrador de Gibbs associado é bastante direto, com código R:

Nsim=5000
xdata=c(5,1,5,14,3,19,1,1,4,22)
Time=c(94.32,15.72,62.88,125.76,5.24,31.44,1.05,1.05,2.10,10.48)
alpha=1.8
gamma=0.01
delta=1
nx=10
lambda=matrix(0,nrow=Nsim,ncol=nx)
beta=matrix(0,nrow=Nsim,ncol=1)
for(i in 2:Nsim){
  for(j in 1:nx) lambda[i,j]=rgamma(1,sh=xdata[j]+alpha,ra=Time[j]+beta[i-1])
  beta[i]=rgamma(1,sh=gamma+nx*alpha,ra=delta+sum(lambda[i,]))}
par(mfrow=c(2,3),mar=c(4,4,1,1))
hist(lambda[,1], col = "grey", breaks = 25, xlab = "", main = expression(lambda[1]), freq = F)
box(); grid()
hist(lambda[,2], col = "grey", breaks = 25, xlab = "", main = expression(lambda[2]), freq = F)
box(); grid()
hist(beta, col = "sienna", breaks = 50, xlab = "", main = expression(beta), freq = F)
box(); grid()
acf(lambda[,1][-1])
grid()
acf(lambda[,2])
grid()
acf(beta)
grid()
Figura. 7.10. Histogramas das distribuições marginais de \(\lambda_1\), \(\lambda_2\) e \(\beta\) dos dados de falha da bomba. Os painéis inferiores correspondentes são gráficos de autocorrelação. Os valores dos hiperparâmetros são \(\alpha= 1.8\), \(\gamma= 0.01\) e \(\delta= 1\).

O resultado de uma execução de 5.000 iterações é mostrado na Figura 7.10.



7.6 Outras considerações

Nesta última seção, examinamos alguns problemas que podem surgir na implementação de um amostrador de Gibbs.


7.6.1 Reparametrização

Muitos fatores contribuem para as propriedades de convergência de um amostrador de Gibbs. Por exemplo, o desempenho da convergência pode ser bastante afetado pela escolha das coordenadas ou seja, pela parametrização. Se a matriz de covariância do alvo tiver uma ampla faixa de autovalores, o amostrador de Gibbs pode ser muito lento para explorar toda a faixa de suporte do alvo.


Exemplo 7.13

Lembre-se do Exemplo 7.1, onde vimos um primeiro amostrador de Gibbs para a normal bivariada em (7.1). Para essa distribuição normal bivariada, a Figura 7.11 mostra a autocorrelação para \(\rho\) = 0.3, 0.6, 0.9.

A correlação mais alta resulta em um amostrador que terá mais problemas para explorar todo o espaço e, portanto, exigirá mais iterações. Também é interessante notar que não importa qual seja o valor de \(\rho\), \(X+Y\) e \(X-Y\) são independentes e, portanto, alterar as coordenadas de \((x,y)\) para \((x+y,x-y)\) levaria a um algoritmo de Gibbs de convergência imediata. 

set.seed(74748)
Nsim=5000 #initial values
#
rho=0.3
X=Y=array(0,dim=c(Nsim,1)) #init arrays
X[1]=rnorm(1,0,1)
Y[1]=rnorm(1,rho*X[1],1-rho^2)
for (i in 2:Nsim){ #sampling loop
  Y[i]=rnorm(1,rho*X[i-1],1-rho^2)
  X[i]=rnorm(1,rho*Y[i],1-rho^2)
}
par(mfrow=c(1,3),mar=c(4,4,1,1))
acf(Y)
text(10,0.7,labels=expression(paste(rho," = 0.3")), cex=1.1)
grid()
rho=0.6
X=Y=array(0,dim=c(Nsim,1)) #init arrays
X[1]=rnorm(1,0,1)
Y[1]=rnorm(1,rho*X[1],1-rho^2)
for (i in 2:Nsim){ #sampling loop
  Y[i]=rnorm(1,rho*X[i-1],1-rho^2)
  X[i]=rnorm(1,rho*Y[i],1-rho^2)
}
acf(Y)
grid()
text(10,0.7,labels=expression(paste(rho," = 0.6")), cex=1.1)
rho=0.9
X=Y=array(0,dim=c(Nsim,1)) #init arrays
X[1]=rnorm(1,0,1)
Y[1]=rnorm(1,rho*X[1],1-rho^2)
for (i in 2:Nsim){ #sampling loop
  Y[i]=rnorm(1,rho*X[i-1],1-rho^2)
  X[i]=rnorm(1,rho*Y[i],1-rho^2)
}
acf(Y)
grid()  
text(10,0.7,labels=expression(paste(rho," = 0.9")), cex=1.1)

Figura. 7.11. Autocorrelações em uma marginal de uma normal bivariada gerada a partir de um amostrador de Gibbs para \(\rho\) = 0.3 (esquerda), 0.6 (meio) e 0.9 (direita).


A convergência dos algoritmos de amostragem de Gibbs e Metropolis-Hastings pode, portanto, sofrer de uma má escolha de parametrização. Como resultado disso, a literatura MCMC tem considerado mudanças na parametrização de um modelo como forma de acelerar a convergência em um amostrador de Gibbs. Parece, no entanto, que a maior parte dos esforços se concentrou no aperfeiçoamento de modelos específicos, resultando na falta de uma metodologia geral para a escolha de uma parametrização “adequada”.

No entanto, o conselho geral é tentar fazer com que os componentes “o mais independente possível” e usar várias parametrizações simultaneamente para misturar os condicionais.


Exemplo 7.13

Lembre-se do Exemplo 7.1, onde vimos um primeiro amostrador de Gibbs para a normal bivariada em (7.1). Para essa distribuição normal bivariada, a Figura 7.11 mostra a autocorrelação para \(\rho= 0.3\); 0.6 e 0.9. A correlação mais alta resulta em um amostrador que terá mais problemas para explorar todo o espaço e, portanto, exigirá mais iterações.

Também é interessante notar que não importa qual seja o valor de \(\rho\), \(X +Y\) e \(X-Y\) são independentes e, portanto, alterar as coordenadas de \((x, y)\) para \((x + y, x-y)\) levaria a um algoritmo de Gibbs imediatamente convergente. 

nsim=10^4
X=Y=rep(0,nsim)
rho=c(.3,.6,.9)
par(mfrow=c(1,3),mar=c(4,4,3,1))
for (t in 1:3){
  Std=sqrt(1-rho[t]^2)
  X[1]=rnorm(1)         #initialize the chain
  Y[1]=rnorm(1)         #initialize the chain
  for(i in 2:nsim){
     X[i]=rnorm(1,mean=rho[t]*Y[i-1],sd=Std)
     Y[i]=rnorm(1,mean=rho[t]*X[i],sd=Std)
     }
  acf(X,lwd=3);grid()
  title(expression(rho==rho[t]),cex=2,col="sienna")
  }

Figura. 7.11. Autocorrelações marginai de uma normal bivariada gerada a partir de um amostrador de Gibbs para \(\rho= 0.3\) (esquerda), \(\rho= 0.6\) (meio) e \(\rho= 0.9\) (direita).


A convergência dos algoritmos de amostragem de Gibbs e Metropolis-Hastings pode, portanto, sofrer de uma má escolha de parametrização. Como resultado disso, a literatura MCMC tem considerado mudanças na parametrização de um modelo como forma de acelerar a convergência em um amostrador de Gibbs. Parece, no entanto, que a maior parte dos esforços se concentrou na melhoria de modelos específicos, resultando na falta de uma metodologia geral para a escolha de uma parametrização adequada.

No entanto, o conselho geral é tentar fazer com que as componentes “tão independente quanto possível” e usar várias parametrizações simultaneamente para misturar as condicionais.


Exemplo 7.14

Uma reparametrização do efeito aleatório unidirecional do Exemplo 7.5 é introduzir a média geral no nível de observação, como em \[\begin{equation} \tag{7.10} \begin{array}{rcl} X_{ij} & \sim & N(\mu+\theta_i, \sigma^2), \qquad i=1,\cdots,k, \quad j=1,\cdots,n_i \\ \theta_i & \sim & N(0,\tau^2), \qquad i=1,\cdots,k, \\ \mu & \sim & N(\mu_0,\sigma^2_\mu)\cdot \end{array} \end{equation}\]

Embora a hierarquia pareça a mesma, os condicionais são diferentes (Exercício 7.14) e as propriedades do amostrador de Gibbs correspondente também. Quando aplicado ao conjunto de dados strong>Energy, o novo amostrador Gibbs não é tão bom. 

set.seed(7493)
nsim = 10^4;a = 10;b = 10
x1 = c(91, 504, 557, 609, 693, 727, 764, 803, 857, 929, 970, 1043, 1089, 1195, 1384, 1713)
x2 = c(457, 645, 736, 790, 899, 991, 1104, 1154, 1203, 1320, 1417, 1569, 1685, 1843, 2296, 2710)
x1 = log(x1)
x2 = log(x2)
n1 = length(x1)
n2 = length(x2)
n = n1 + n2
x1bar = mean(x1)
x2bar = mean(x2)
xbar = (n1 * x1bar + n2 * x2bar)/n
k = 2
a1 = a2 = a3 = a
b1 = b2 = b3 = b
mu = array(xbar, c(nsim, 1))
theta1 = array(x1bar, c(nsim, 1))
theta2 = array(x2bar, c(nsim, 1))
sigma2mu = array((var(x1) + var(x2))/2, c(nsim, 1))
sigma2 = array((var(x1) + var(x2))/2, c(nsim, 1))
tau2 = array((var(x1) + var(x2))/2, c(nsim, 1))
muR = array(xbar, c(nsim, 1))
theta1R = array(x1bar, c(nsim, 1))
theta2R = array(x2bar, c(nsim, 1))
sigma2muR = array((var(x1) + var(x2))/2, c(nsim, 1))
sigma2R = array((var(x1) + var(x2))/2, c(nsim, 1))
tau2R = array((var(x1) + var(x2))/2, c(nsim, 1))
mu0 = xbar
for (i in 2:nsim) {
  B1 = sigma2[i - 1]/(sigma2[i - 1] + n1 * tau2[i - 1]) 
    theta1[i] = rnorm(1, mean = B1 * mu[i - 1] + (1 - B1) * x1bar, sd = sqrt(tau2[i - 1] * B1))
    B2 = sigma2[i - 1]/(sigma2[i - 1] + n2 * tau2[i - 1])
        theta2[i] = rnorm(1, mean = B2 * mu[i - 1] + (1 - B2) * x2bar, sd = sqrt(tau2[i - 1] * B2))
    B1R = sigma2R[i - 1]/(sigma2R[i - 1] + n1 * tau2R[i - 1])
        theta1R[i] = rnorm(1, mean = (1 - B1R) * x1bar, sd = sqrt(tau2R[i - 1] * B1R))
    B2R = sigma2R[i - 1]/(sigma2R[i - 1] + n2 * tau2R[i - 1])
        theta2R[i] = rnorm(1, mean = (1 - B2R) * x2bar, sd = sqrt(tau2R[i - 1] * B2R))
    B = tau2[i - 1]/(k * sigma2mu[i - 1] + tau2[i - 1])
    m = B * mu0 + (1 - B) * (n1 * theta1[i] + n2 * theta2[i])/n
    mu[i] = rnorm(1, mean = m, sd = sqrt(sigma2mu[i - 1] * B))
    BR = sigma2muR[i - 1]/(sigma2muR[i - 1] + sigma2R[i - 1]/n)
    m = (1 - BR) * mu0 + BR * (xbar - (n1 * theta1R[i] + n2 * theta2R[i])/n)
    muR[i] = rnorm(1, mean = m, sd = sqrt(sigma2R[i - 1] * BR/n))
    sh1 = (n/2) + a1
    ra1 = (1/2) * (sum((x1 - theta1[i])^2) + sum((x2 - theta2[i])^2)) + b1
    sigma2[i] = 1/rgamma(1, shape = sh1, rate = ra1)
    ra1 = (1/2) * (sum((x1 - muR[i] - theta1R[i])^2) + sum((x2 - muR[i] - theta2R[i])^2)) + b1
    sigma2[i] = 1/rgamma(1, shape = sh1, rate = ra1)
    sh2 = (k/2) + a2
    ra2 = (1/2) * (sum((x1bar - theta1[i] - mu[i])^2) + sum((theta2[i] - mu[i])^2)) + b2
    tau2[i] = 1/rgamma(1, shape = sh2, rate = ra2)
    ra2 = (1/2) * (theta1R[i]^2 + theta2R[i]^2) + b2
    tau2[i] = 1/rgamma(1, shape = sh2, rate = ra2)
    sh3 = (1/2) + a3
    ra3 = (1/2) * (mu[i] - mu0)^2 + b3
    sigma2mu[i] = 1/rgamma(1, shape = sh3, rate = ra3)
}
ra3 = (1/2) * (muR[i] - mu0)^2 + b3
sigma2mu[i] = 1/rgamma(1, shape = sh3, rate = ra3)
lim = c(-0.05, 0.1)
par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))
acf(mu, ylim = lim, lwd = 2, main = expression(mu), cex = 3);grid()
acf(theta1, ylim = lim, lwd = 2, main = expression(theta[1]));grid()
acf(theta2, ylim = lim, lwd = 2, main = expression(theta[2]));grid()
acf(muR, ylim = lim, lwd = 2, main = expression(mu));grid()
acf(theta1R, ylim = lim, lwd = 2, main = expression(theta[1]));grid()
acf(theta2R, ylim = lim, lwd = 2, main = expression(theta[2]));grid()

Figura. 7.12. Gráficos de autocovariância para o amostrador de Gibbs associado ao modelo (7.7) e o amostrador de Gibbs associado à sua reparametrização (7.10). A linha superior fornece as autocovariâncias para \(\mu\), \(\theta_1\) e \(\theta_2\) (da esquerda para a direita) para o modelo (7.7), e a linha inferior os fornece para o modelo (7.10).

Por exemplo, a Figura 7.12 mostra as autocorrelações, que, a olho nu, parecem apenas ligeiramente melhores para o primeiro modelo. No entanto, se olharmos para a matriz de covariância da subcadeia \(\big(\mu^{(t)},\theta_1^{(t)},\theta_2^{(t)}\big)\), sua estimativa é

print(cov(cbind(mu, theta1, theta2)))
##             [,1]         [,2]         [,3]
## [1,] 0.743021406 0.0066109797 0.0089859869
## [2,] 0.006610980 0.0384451022 0.0001737657
## [3,] 0.008985987 0.0001737657 0.0395439286
print(cov(cbind(muR, theta1R, theta2R)))
##              [,1]          [,2]          [,3]
## [1,]  0.022742682 -9.014282e-03 -9.034942e-03
## [2,] -0.009014282  1.948148e-02 -7.692061e-05
## [3,] -0.009034942 -7.692061e-05  1.913127e-02
para o modelo (7.7) e modelo (7.10), respectivamente, portanto as variâncias e covariâncias são maiores para o modelo reparametrizado. Assim, claramente devemos usar a parametrização do modelo (7.10).



7.6.2 Rao-Blackwellization

Já vimos a Rao-Blackwellização na Seção 4.6, onde o condicionamento em um subconjunto das variáveis simuladas pode produzir melhorias consideráveis no estimador empírico padrão em termos de variância por uma simples “reciclagem” das variáveis rejeitadas.

No entanto, como o amostrador Gibbs aceita todos os valores simulados, esse tipo de reciclagem não pode ser aplicado. No entanto, Gelfand and Smith (1990) propõem um tipo de condicionamento que chamaremos de Rao-Blackwellization paramétrico para diferenciá-la da forma estudada na Seção 4.6.

Para \((X,Y)\sim f(x,y)\), a Rao–Blackwellização paramétrica é baseada na identidade de marginalização, esperança iterada \[ \mbox{E}(X)=\mbox{E}\big(\mbox{E}(X \, | \, Y) \big)\cdot \] Definindo \(\delta (Y) = \mbox{E}([X \, | \, Y\big)\), temos \(\mbox{E}\big(\delta(Y)\big) = \mbox{E}(X)\) e \(\mbox{Var}\big(\delta(Y)\big)\leq \mbox{Var}(X)\), mostrando que \(\delta(Y)\) é melhor estimador, desde que possa ser calculado.


Exemplo 7.15 (Continuação do Exemplo 7.1)

No caso em que o alvo é a distribuição normal bivariada, as condicionais completas são \[ X_{t+1} \, | \, y_t \sim N(\rho y_t, 1-\rho^2), \\ Y_{t+1} \, | \, x_{t+1} \sim N(\rho x_{t+1} , 1-\rho^2), \] e assim segue que \(\mbox{E}\big(X \, | \, Y\big) = \rho Y\). Como \(X\) e \(Y\) têm a mesma distribuição marginal, a variância da versão Rao–Blackwellizada é obviamente reduzida por um fator \(\rho^2\).


Infelizmente, a redução da variância do uso de \(\delta_Y\) não é válida em geral, devido à correlação na amostra MCMC. No entanto, Liu et al. (1994) mostraram que, em particular, a melhoria vale para qualquer amostrador Gibbs de dois estágios.

Agora, examinamos outro exemplo de Rao–Blackwellization em um amostrador Gibbs de dados ausentes para uma ocorrência comum em que ocorrem possíveis ganhos.


Exemplo 7.16

Para 360 unidades de tempo consecutivas, considere registrar o número de passagens de indivíduos por unidade de tempo por algum sensor. Isso pode ser, por exemplo, o número de carros observados em uma encruzilhada.

Os resultados hipotéticos são

\[ \begin{array}{lccccc}\hline \mbox{Número de passagens} & 0 & 1 & 2 & 3 & \mbox{4 ou mais} \\\hline \mbox{Número de observações} & 139 & 128 & 55 & 25 & 13 \\\hline \end{array} \]

Os dados envolvem um agrupamento das observações com quatro passagens ou mais. Isso pode ser tratado como um modelo de dados perdidos, onde assumimos que as observações não agrupadas são \(X_i\sim Poisson(\lambda)\). A verossimilhança do modelo é \[ \ell(\lambda \, | \, x_1,\cdots,x_5) \approx e^{-347 \lambda} \lambda^{128+55+25+13} \left( 1-e^{-\lambda} \sum_{i=0}^3 \lambda^i/i!\right)^{13}, \] para \(x_1=139\),\(\cdots\),\(x_5=13\).

Para \(\pi(\lambda)=1/\lambda\) e \(z=(z_1,\cdots,z_{13})\) o vetor das 13 unidades maiores que 4, podemos derivar um amostrador de Gibbs de conclusão dos condicionais completos \[ Z_i^{(t)}\sim Poisson \big(\lambda^{(t-1)} \big)\pmb{I}_{\{y\geq 4\}} \] para \(i=1,\cdots,13\) e \[ \lambda^{(t)}\sim Gama\left(313+\sum_{i=1}^{13} Z_i^{(t)}, 360 \right)\cdot \]

A estimativa Rao–Blackwellizada de \(\lambda\) é então dada por \[ \mbox{E}\left( \lambda \, | \, z_1^{(t)},\cdots,z_{13}^{(t)}\right) = \dfrac{1}{360 T}\sum_{t=1}^T \left( 313+\sum_{i=1}^{13} y_i^{(t)}\right), \] e a evolução deste estimador, juntamente com a média empírica, é apresentada na Figura 7.13. Ele exibe uma enorme redução de variância. 

set.seed(824)
nsim=10^3
lam=RB=rep(313/360,nsim)
z=rep(0,13)
for (j in 2:nsim){
  top=round(lam[j -1]+6*sqrt(lam[j -1]))
  prob=dpois(c(4:top),lam[j -1])
  cprob=cumsum(prob/sum(prob))
  for(i in 1:13){z[i] = 4+sum(cprob<runif(1))}
  RB[j]=(313+sum(z))/360
  lam[j]=rgamma(1,360*RB[j],scale=1/360);
}
par(mfrow=c(1,3),mar=c(4,4,2,1))
hist(lam,col="grey",breaks=25,xlab="",main="Empirical Average");box()
plot(cumsum(lam)/c(1:nsim),ylim=c(1,1.05),type="l",lwd=1.5,ylab="")
lines(cumsum(RB)/c(1:nsim),col="sienna",lwd=1.5);grid()
hist(RB,col="sienna",breaks=62,xlab="",main="Rao-Blackwell",xlim=c(1,1.05));box()

Figura. 7.13. Para os dados de contagem do Exemplo 7.16, o histograma de \(\lambda\) (esquerda) e sua esperança condicional \(\mbox{E}(\lambda \, | \, Z) = 313 + \sum_{i=1}^{13} z_i\) (direita). Observe a diferença na escala dos histogramas. O painel central mostra a evolução das médias acumuladas da média empírica (preto) e Rao–Blackwellization (cinza).


Também existem estimadores Rao–Blackwellizados não paramétricos em configurações de variáveis ausentes. Ao considerar uma aproximação para a distribuição marginal \(f_{_X}\) associada a \(f(x,y_1,\cdots,y_p)\), um estimador Rao–Blackwellizado associado à cadeia de Gibbs \(\{(x^{(t)},y^{(t)})\}\) é dado por \[\begin{equation} \tag{7.11} \widehat{f}_{_X}(x)=\dfrac{1}{T}\sum_{t=1}^T f(x \, | \, y^{(t)}), \end{equation}\] que converge a uma velocidade paramétrica para \(f_{_X}(x)\). Este estimador dá uma aproximação suave ao marginal, que pode ser mostrado nos valores finais.

Conforme estudado nos Exercícios 3.15, 4.1 e 4.2, a aproximação dos fatores de Bayes exige soluções específicas. Chib (1995) propõe uma abordagem alternativa baseada em Rao-Blackwellization que é muito mais eficiente quando pode ser implementado.


7.6.3 Metropolis dentro de Gibbs e estratégias híbridas

Um ponto que vale a pena enfatizar sobre a implementação de um amostrador de Gibbs é que ele pode ser facilmente estendido para configurações onde alguns dos condicionais completos não podem ser simulados por geradores aleatórios padrão

Se, dentro de um conjunto de condicionais completos \(f_1,\cdots,f_p\), alguma densidade \(f_i\) não é convencional, por exemplo (5.15) no Exemplo 5.17, isso não prejudica o amostrador de Gibbs resultante no sentido de que a seguinte estratégia Metropolis-dentro-Gibbs pode ser adotada.

Em vez de simular \[ X_i^{(t+1)}\sim f_i\big( x_i \, | \, x_1^{t+1},\cdots,x_{i-1}^{t+1},x_{i+1}^{(t)},\cdots,x_p^{(t)}\big), \] você pode executar uma única etapa de qualquer esquema MCMC associado à distribuição estacionária \[ f_i\big( x_i \, | \, x_1^{t+1},\cdots,x_{i-1}^{t+1},x_{i+1}^{(t)},\cdots,x_p^{(t)}\big)\cdot \]

Uma solução simples é, por exemplo, usar um algoritmo Metropolis de passeio aleatório centrado em \(x_i^{(t)}\). Embora a princípio isso pareça uma aproximação grosseira, já que a condicional completa não é exatamente simulada, a validade do algoritmo resultante é exatamente a mesma do amostrador de Gibbs original, pois a distribuição conjunta \(f\) permanece a distribuição estacionária da Cadeia de Markov correspondente.

Você pode então se perguntar qual é o sentido de usar um amostrador de Gibbs se as simulações de componente precisam ser substituídas por etapas Metropolis-Hastings, pois usar um algoritmo Metropolis-Hastings direcionado à distribuição conjunta \(f\) é mais natural. Embora não haja nada restritivo impedir que você use um algoritmo conjunto de Metropolis-Hastings, na maioria das vezes, projetar um algoritmo de Metropolis-Hastings em um alvo de grande dimensão é desafiador ou até mesmo impossível.

O ganho fundamental em usar uma estrutura do tipo Gibbs é que ele divide um modelo complexo em um grande número de alvos menores e mais simples, onde os algoritmos locais de Metropolis-Hastings podem ser projetados com pouco custo.


Exemplo 7.17

Consideremos a distribuição alvo (5.15), mencionada no Exemplo 5.17, sendo que esta não é uma distribuição padrão.

Enquanto Booth and Hobert (1999) projetaram um algoritmo específico de Accept-Reject para simular de (5.15), uma proposta de caminhada aleatória em cada interface do usuário, como em


for (i in 1:n){
  mu=u[i]
  u[i]=factor*sigma[iter-1]*rnorm(1)+mu
  if (log(runif(1))>gu(u[i],i,beta[iter-1],sigma[iter-1])-
      + gu(mu,i,beta[iter-1],sigma[iter-1])){
    u[i]=mu
    }
  }
produz uma amostra de \(u_i\)’s na iteração iter condicional nos valores atuais dos parâmetros e a amostra de \(u_i\)’s na iteração iter-1.

No amostrador Gibbs geral, os parâmetros são então simulados por


sigma=c(sigma,1/sqrt(2*rgamma(1,0.5*n)/sum(u^2)))
tau=sigma[iter]/sqrt(sum(as.vector(x^2)*pro(beta[iter-1],u)))
betaprop=beta[iter-1]+rnorm(1)*factor*tau
if (log(runif(1))>likecomp(betaprop,sigma[iter],u)-
    likecomp(beta[iter-1],sigma[iter],u))
  betarop=beta[iter-1]
beta=c(beta,betaprop)
de maneira direta. Veja o Exemplo 8.1 para a implementação completa.

O termo de calibração factor pode ainda ser ajustado em relação à taxa de aceitação da Seção 6.5, conforme descrito na Seção 8.5.


Embora permaneçamos próximos a essa ideia de incorporar os passos de Metropolis-Hastings quando a simulação direta não é possível, podemos também sinalizar a possível extensão para estratégias híbridas. As estratégias híbridas não devem ser confundidas com o Monte Carlo híbrido (Neal, 1999), também chamado de MCMC Hamiltoniano, que é uma forma de implementação de Langevin que visa reduzir o desperdício de simulação em propostas de passeio aleatório.

O conceito é mais uma vez baseado na estacionaridade da distribuição de destino correta, mesmo que a intuição possa discordar. Quando dado um alvo, univariado ou multivariado, onde vários esquemas MCMC naturais estão disponíveis, um algoritmo híbrido mescla esses diferentes esquemas completamente.

Esquematicamente, se as funções de atualização MCMC local ou global, significando componentes ou conjuntos, mcmc.1(x,y), \(\cdots\), mcmc.q(x,y) estiverem disponíveis, o kernel de transição definido por


mcmc(x,y)=function(x,y){
  switch(sample(1:p,1),
  mcmc.1(x,y)
    ...
  mcmc.p(x,y))
}
permanece uma função de atualização MCMC válida na mesma distribuição de destino.

Embora isso pareça uma ideia ridícula porque os esquemas ruins são misturados com os bons, a mistura cega de todas as estratégias disponíveis é, no entanto,

  1. válida da perspectiva de produzir a distribuição estacionária correta e

  2. livre de riscos no sentido de que se a lista de funções contiver um único algoritmo de bom desempenho, a versão híbrida terá pelo menos o mesmo desempenho, exigindo simplesmente uma extensão do tempo de computação.

Por exemplo, se várias estratégias de bloqueio ou reparametrização estiverem disponíveis simultaneamente, elas podem ser todas incorporadas no mesmo algoritmo. Essa solução pode parecer uma perda de tempo de computação, mas nosso conselho sobre esse assunto é que, a menos que algumas das funções mcmc.i funcionem claramente, o tempo gasto (desperdiçado) executando a solução híbrida é o tempo economizado em projetar e selecionar as funções mcmc.i mais eficientes. Em outras palavras, é mais eficiente deixar o computador classificar entre as soluções disponíveis do que executar testes preliminares para classificar essas soluções manualmente.


7.6.4 A prioris impróprias

Esta seção discute um perigo específico resultante do uso descuidado do amostrador Gibbs. Sabemos que o amostrador de Gibbs é baseado em distribuições condicionais derivadas da distribuição conjunta.

No entanto, o que é particularmente insidioso é que essas distribuições condicionais podem ser bem definidas e podem ser simuladas, mas podem não corresponder a nenhuma distribuição conjunta!

Este problema não é um defeito do amostrador de Gibbs ou mesmo um problema de simulação, mas sim um problema de uso inadvertido do amostrador de Gibbs em uma situação em que as suposições subjacentes são violadas. No entanto, é importante alertar o usuário de algoritmos MCMC contra esse perigo, pois corresponde a uma situação frequentemente encontrada em modelos bayesianos não informativos ou “padrão”.

A construção do amostrador de Gibbs diretamente das distribuições condicionais é um forte incentivo para ignorar a verificação da propriedade da a posteriori, especialmente em configurações complexas. Mas essa verificação é essencial, como mostra o exemplo simples a seguir.


Exemplo 7.18

O seguinte modelo foi utilizado por Casella and George (1992) para apontar a dificuldade de avaliar a impropriedade de uma distribuição a posteriori por meio das distribuições condicionais.

O par de densidades condicionais \[\begin{equation} \tag{7.12} X \, | \, y \sim Exp(y) \qquad Y \, | \, x \sim Exp(x), \end{equation}\] são distribuições condicionais bem definidas, mas essas distribuições condicionais não correspondem a nenhuma distribuição de probabilidade conjunta.

A Figura 7.14 mostra um histograma e uma média cumulativa para uma amostra gerada usando o amostrador de Gibbs correspondente a essas condicionais. Os gráficos são extremamente curiosas e na verdade são um lixo absoluto! Esta não é uma Cadeia de Markov recorrente. De fato, a única função que poderia ser a distribuição conjunta é \[ f(x,y) \approx \exp(-xy), \] que não possui integral finita. 

nsim=10^3
X=Y=rep(0,nsim)
X[1]=rexp(1)            #initialize the chain
Y[1]=rexp(1)            #initialize the chain
for (i in 2:nsim){
  X[i]=rexp(1,rate=Y[i-1])
  Y[i]=rexp(1,rate=X[i])
}
st=0.1*nsim
par(mfrow=c(1,2),mar=c(4,4,2,1))
hist(X,col="grey",breaks=25,xlab="",main="")
box()/grid()
## numeric(0)
plot(cumsum(X)[(st+1):nsim]/c(1:(nsim-st)),type="l",lwd=1.5,ylab="",xlab="Iterations")
grid()

Figura. 7.14. Histograma e média cumulativa da variável \(X\) do amostrador de Gibbs em (7.12). Observe os intervalos nos gráficos que sinalizam problemas de convergência.


Dados os resultados do Exemplo 7.18, pode parecer que um simples monitoramento gráfico é suficiente para exibir um comportamento desviante do amostrador de Gibbs. No entanto, este não é o caso em geral e há muitos exemplos, alguns dos quais são publicados (ver Casella, 1996), onde a saída do amostrador de Gibbs aparentemente não difere de uma cadeia de Markov convergente.

Freqüentemente, esse fenômeno ocorre quando a divergência da densidade a posteriori ocorre “em 0”; isto é, em um ponto específico cuja vizinhança imediata raramente é visitada pela cadeia, como no exemplo de efeitos aleatórios a seguir. A única maneira de garantir que o amostrador de Gibbs que você está usando é válido é verificar se a distribuição conjunta tem uma integral finita.


Exemplo 7.19

Considere um modelo de efeitos aleatórios, \[ Y_{ij}=\beta+U_i+\epsilon_{ij}, \quad i=1,\cdots,I, \quad j=1,\cdots,J, \] onde \(U_i\sim N(0,\sigma^2)\) e \(\epsilon_{ij}\sim N(0,\tau^2)\).

À a priori de Jeffreys, imprópria para os parâmetros \(\beta\), \(\sigma\) e \(\tau\) é \[ \pi(\beta,\sigma^2,\tau^2) = \dfrac{1}{\sigma^2\tau^2}\cdot \]

As distribuições condicionais \[ \begin{array}{rcl} U_i \, | \, y,\beta,\sigma^2,\tau^2 & \sim & N\left( \dfrac{J(\overline{y}_i-\beta)}{J+\tau^2/\sigma^{2}}, \big( J/\tau^2 +1/\sigma^2 \big)^{-1}\right), \\ \beta \, | \, u,y,\sigma^2,\tau^2 & \sim & N\left( \overline{y}-\overline{u}, \tau^2/JI\right), \\ \sigma^2 \, | \, u,y,\beta,\tau^2 & \sim & \displaystyle \mathcal{IG}\left( I/2, \big(1/2\big) \sum_{i} u_i^2\right), \\ \tau^2 \, | \, u,y,\beta,\sigma^2 & \sim & \displaystyle \mathcal{IG}\left( IJ/2, \big(1/2\big) \sum_{ij} \big( y_{ij}-u_i-\beta\big)^2\right), \end{array} \] são bem definidas e um amostrador Gibbs pode ser facilmente implementado nessa configuração.

No entanto, não há distribuição conjunta adequada que corresponda a essas condicionais! E, em muitos casos, como você pode verificar por si mesmo, isso é impossível de detectar monitorando a saída.


Se a priori impróprios forem usados em um amostrador Gibbs, a posteriori deve sempre ser verificada quanto à propriedade. No entanto, muitas vezes, as priori impróprias causam mais problemas sobre as variâncias do que sobre as médias.


7.7 Exercícios