4. Computação Bayesiana

Existem várias formas de resumir a informação descrita na distribuição a posteriori. Esta etapa frequentemente envolve a avaliação de probabilidades ou esperanças.

Neste capítulo serão descritos métodos baseados em simulação, incluindo Monte Carlo simples, Monte Carlo com função de importância, o método do Bootstrap Bayesiano e Monte Carlo via cadeias de Markov (MCMC). O material apresentado é introdutório e mais detalhes sobre os estes métodos podem ser obtidos em Gamerman (1997), Davison e Hinckley (1997) e Robert e Casella (1999). Outros métodos computacionalmente intensivos como técnicas de otimização e integração numérica, bem como aproximações analíticas não serão tratados aqui e uma referência introdutória é Migon e Gamerman (1999).

Todos os algoritmos que serão vistos aqui são não determinísticos, i.e. todos requerem a simulação de números (pseudo) aleatórios de alguma distribuição de probabilidades. Em geral, a única limitação para o número de simulações são o tempo de computação e a capacidade de armazenamento dos valores simulados. Assim, se houver qualquer suspeita de que o número de simulações é insuficiente, a abordagem mais simples consiste em simular mais valores.

4.1 Uma Palavra de Cautela

Apesar da sua grande utilidade, os métodos que serão apresentados aqui devem ser aplicados com cautela. Devido à facilidade com que os recursos computacionais podem ser utilizados hoje em dia, corremos o risco de apresentar uma solução para o problema errado (o erro tipo 3) ou uma solução ruim para o problema certo. Assim, os métodos computacionalmente intensivos não devem ser vistos como substitutos do pensamento crítico sobre o problema por parte do pesquisador.

Além disso, sempre que possível deve-se utilizar soluções exatas, i.e. não aproximadas, se elas existirem. Por exemplo, em muitas situações em que precisamos calcular uma integral múltipla existe solução exata em algumas dimensões, enquanto nas outras dimensões temos que usar métodos de aproximação.

4.2 O Problema Geral da Inferência Bayesiana

A distribuição a posteriori pode ser convenientemente resumida em termos de esperanças de funções particulares do parâmetro $\theta$, i.e.

$$E[g(\theta)\vert\bfx]=\int g(\theta) p(\theta\vert\bfx) d\theta$$

ou distribuições a posteriori marginais quando $\btheta$ for multidimensional, i.e.

$$p(\btheta_1\vert\bfx) =\int p(\btheta\vert\bfx) d\btheta_2$$

onde $\btheta=(\btheta_1,\btheta_2)$.

Assim, o problema geral da inferência Bayesiana consiste em calcular tais valores esperados segundo a distribuição a posteriori de $\theta$. Alguns exemplos são,

1. Constante normalizadora. $g(\theta)=1$ e $p(\theta\vert bfx)=kq(\theta)$, segue que
   $$k=\left[\int q(\theta) d\theta\right]^{-1}.$$

2. Se $g(\theta) = \theta$, então têm-se $\mu = E(\theta\vert\bfx)$, média a posteriori.

3. Quando $g(\theta) = (\theta - \mu)^2$, então $\sigma^2 = var(\theta)=E((\theta - \mu)^2\vert\bfx)$, a variância a posteriori.

4. Se $g(\theta) = I_A(\theta)$, onde $I_A(x) = 1$ se $x \in A$ e zero caso contrário, então $P(A~\vert~\bfx) = \int_A p(\theta\vert\bfx)d\theta$

5. Seja $g(\theta) = p(y\vert\theta)$, onde $y\perp\bfx\vert\theta$. Nestas condições obtemos $E[p(y\vert\bfx)]$, a distribuição preditiva de $y$, uma observação futura.

Portanto, a habilidade de integrar funções, muitas vezes complexas e multidimensionais, é extremamente importante em inferência Bayesiana. Inferência exata somente será possível se estas integrais puderem ser calculadas analiticamente, caso contrário devemos usar aproximações. Nas próximas seções iremos apresentar métodos aproximados baseados em simulação para obtenção dessas integrais.

4.3 Método de Monte Carlo Simples

A idéia do método é justamente escrever a integral que se deseja calcular como um valor esperado. Para introduzir o método considere o problema de calcular a integral de uma função $g(\theta)$ no intervalo $(a,b)$, i.e.

$$I=\int_a^b g(\theta)d\theta.$$

Esta integral pode ser reescrita como

$$I=\int_a^b (b-a)g(\theta)\frac{1}{b-a}d\theta = (b-a)E[g(\theta)]$$

identificando $\theta$ como uma variável aleatória com distribuição $U(a,b)$. Assim, transformamos o problema de avaliar a integral no problema estatístico de estimar uma média, $E[g(\theta)]$. Se dispomos de uma amostra aleatória de tamanho $n$, $\theta_1,\dots,\theta_n$ da distribuição uniforme no intervalo $(a,b)$ teremos também uma amostra de valores $g(\theta_1),\dots,g(\theta_n)$ da função $g(\theta)$ e a integral acima pode ser estimada pela média amostral, i.e.

$$\hat{I}=(b-a)\frac{1}{n}\sum_{i=1}^n g(\theta_i).$$

Não é difícil verificar que esta estimativa é não viesada já que

$$E(\hat{I})=\frac{(b-a)}{n}\sum_{i=1}^n E[g(\theta_i)]= (b-a)E[g(\theta)]=\int_a^b g(\theta)d\theta.$$

Podemos então usar o seguinte algoritmo

1. gere $\theta_1,\dots,\theta_n$ da distribuição $U(a,b)$;
2. calcule $g(\theta_1),\dots,g(\theta_n)$;
3. calcule a média amostral $\overline{g}=\sum_{i=1}^n g(\theta_i)/n$
4. calcule $\hat{I}=(b-a)\overline{g}$

A generalização é bem simples para o caso em que a integral é a esperança matemática de uma função $g(\theta)$ onde $\theta$ tem função de densidade $p(\theta)$, i.e.

$$I=\int_a^b g(\theta)p(\theta)d\theta = E[g(\theta)].$$ (4.1)

Neste caso, podemos usar o mesmo algoritmo descrito acima modificando o passo 1 para gerar $\theta_1,\dots,\theta_n$ da distribuição $p(\theta)$ e calculando $\hat{I}=\overline{g}$.

Uma vez que as gerações são independentes, pela Lei Forte dos Grandes Números segue que $\hat{I}$ converge quase certamente para $I$. Além disso, a variância do estimador pode também ser estimada como

$$v = \frac{1}{n^2}\sum_{i=1}^n (g(\theta_i)-\overline{g})^2,$$

i.e. a aproximação pode ser tão acurada quanto se deseje bastando aumentar o valor de $n$. É importante notar que $n$ está sob nosso controle aqui, e não se trata do tamanho da amostra de dados.

Para $n$ grande segue que

$$\frac{\overline{g}-E[g(\theta)]}{\sqrt{v}}$$

tem distribuição aproximadamente $N(0,1)$. Podemos usar este resultado para testar convergência e construir intervalos de confiança.

No caso multivariado a extensão também é direta. Seja $\bftheta=(\theta_1,\dots,\theta_k)'$ um vetor aleatório de dimensão $k$ com função de densidade $p(\bftheta)$. Neste caso os valores gerados serão também vetores $\bftheta_1,\dots,\bftheta_n$ e o estimador de Monte Carlo fica

$$\hat{I}=\frac{1}{n}\sum_{i=1}^n g(\bftheta_i)$$

4.3.1 Monte Carlo via Função de Importância

Em muitas situações pode ser muito custoso ou mesmo impossível simular valores da distribuição a posteriori. Neste caso, pode-se recorrer à uma função $q(\theta)$ que seja de fácil amostragem, usualmente chamada de função de importância. O procedimento é comumente chamado de amostragem por importância.

Se $q(\theta)$ for uma função de densidade definida no mesmo espaço variação de $\theta$ então a integral (4.1) pode ser reescrita como

$$I=\int\frac{g(\theta)p(\theta)}{q(\theta)}q(\theta)dx = E\left[\frac{g(\theta)p(\theta)}{q(\theta)}\right]$$

onde a esperança agora é com respeito a distribuição $q$. Assim, se dispomos de uma amostra aleatória $\theta_1,\dots,\theta_n$ tomada da distribuição $q$ o estimador de Monte Carlo da integral acima fica

$$\hat{I}=\frac{1}{n}\sum_{i=1}^n\frac{g(\theta_i)p(\theta_i)}{q(\theta_i)}.$$

e tem as mesmas propriedades do estimador de Monte Carlo simples.

Em princípio não há restrições quanto à escolha da densidade de importância $q$, porém na prática alguns cuidados devem ser tomados. Pode-se mostrar que a escolha ótima no sentido de minimizar a variância do estimador consiste em tomar $q(\theta)\propto g(\theta)p(\theta)$.

Para uma única observação $X$ com distribuição $N(\theta,1)$, $\theta$ desconhecido, e priori Cauchy(0,1) segue que

$$p(x\vert\theta)\propto\exp[-(x-\theta)^2/2]$$   e$$\quad p(\theta)=\frac{1}{\pi(1+\theta^2)}.$$

Portanto, a densidade a posteriori de $\theta$ é dada por

$$p(\theta\vert x)= \frac{\dfrac{1}{1+\theta^2}\exp[-(x-\theta)^2/2]} {\dint\dfrac{1}{1+\theta^2}\exp[-(x-\theta)^2/2]d\theta}.$$

Suponha agora que queremos estimar $\theta$ usando função de perda quadrática. Como vimos no Capítulo 3 isto implica em tomar a média a posteriori de $\theta$ como estimativa. Mas

$$E[\theta\vert\bfx] = \int\theta p(\theta\vert\bfx)d\theta=\frac {... ...\theta)^2/2]d\theta} {\dint\dfrac{1 }{1+\theta^2}\exp[-(x-\theta)^2/2]d\theta}$$

e as integrais no numerador e denominador não têm solução analítica exata. Uma solução aproximada via simulação de Monte Carlo pode ser obtida usando o seguinte algoritmo,

1. gerar $\theta_1,\dots,\theta_n$ independentes da distribuição $N(x,1)$;
2. calcular $g_i=\dfrac{\theta_i}{1+\theta_i^2}$ e $g_i^*=\dfrac{1}{1+\theta_i^2}$;
3. calcular $\hat{E}(\theta\vert\bfx)=\dfrac{\sum_{i=1}^n g_i}{\sum_{i=1}^n g_i^*}$.

Este exemplo ilustrou um problema que geralmente ocorre em aplicações Bayesianas. Como a posteriori só é conhecida a menos de uma constante de proporcionalidade as esperanças a posteriori são na verdade uma razão de integrais. Neste caso, a aproximação é baseada na razão dos dois estimadores de Monte Carlo para o numerador e denominador.

Exercícios

1. Para uma única observação $X$ com distribuição $N(\theta,1)$, $\theta$ desconhecido, queremos fazer inferência sobre $\theta$ usando uma priori Cauchy(0,1). Gere um valor de $X$ para $\theta=2$, i.e. $x\sim N(2,1)$.
   1. Estime $\theta$ através da sua média a posteriori usando o algoritmo do exemplo 1.
   2. Estime a variância da posteriori.
   3. Generalize o algoritmo para $k$ observações $X_1,\dots,X_k$ da distribuição $N(\theta,1)$.

4.4 Métodos de Reamostragem

Existem distribuições para as quais é muito difícil ou mesmo impossível simular valores. A idéia dos métodos de reamostragem é gerar valores em duas etapas. Na primeira etapa gera-se valores de uma distribuição auxiliar conhecida. Na segunda etapa utiliza-se um mecanismo de correção para que os valores sejam representativos (ao menos aproximadamente) da distribuição a posteriori. Na prática costuma-se tomar a priori como distribuição auxiliar conforme proposto em Smith e Gelfand (1992).

4.4.1 Método de Rejeição

Considere uma densidade auxiliar $q(\theta)$ da qual sabemos gerar valores. A única restrição é que exista uma constante $A$ finita tal que $p(\theta\vert\bfx)<A q(\theta)$. O método de rejeição consiste em gerar um valor $\theta^*$ da distribuição auxiliar $q$ e aceitar este valor como sendo da distribuição a posteriori com probabilidade $p(\theta\vert\bfx)/A q(\theta)$. Caso contrário, $\theta^*$ não é aceito como uma valor gerado da posteriori e o processo é repetido até que um valor seja aceito. O método também funciona se ao invés da posteriori, que em geral é desconhecida, usarmos a sua versão não normalizada, i.e $p(\bfx\vert\theta)p(\theta)$.

Tomando a priori $p(\theta)$ como densidade auxiliar a constante $A$ deve ser tal que $p(\bfx\vert\theta)<A$. Esta desigualdade é satisfeita se tomarmos $A$ como sendo o valor máximo da função de verossimilhança, i.e. $A=p(\bfx\vert\hat{\theta})$ onde $\hat{\theta}$ é o estimador de máxima verossimilhança de $\theta$. Neste caso, a probabilidade de aceitação se simplifica para $p(\bfx\vert\theta)/p(\bfx\vert\hat{\theta})$.

Podemos então usar o seguinte algoritmo para gerar valores da posteriori

1. gerar um valor $\theta^*$ da distribuição a priori;
2. gerar $u\sim U(0,1)$;
3. aceitar $\theta^*$ como um valor da posteriori se $u < p(\bfx\vert\theta^*)/p(\bfx\vert\hat{\theta})$, caso contrário rejeitar $\theta^*$ e retornar ao item 1.

Um problema técnico associado ao método é a necessidade de se maximizar a função de verossimilhança o que pode não ser uma tarefa simples em modelos mais complexos. Se este for o caso então o método de rejeição perde o seu principal atrativo que é a simplicidade. Neste caso, o método da próxima seção passa a ser recomendado.

Outro problema é que a taxa de aceitação pode ser muito baixa, i.e. teremos que gerar muitos valores da distribuição auxiliar até conseguir um número suficiente de valores da posteriori. Isto ocorrerá se as informações da priori e da verossimilhança forem conflitantes já que neste caso os valores gerados terão baixa probabilidade de serem aceitos.

4.4.2 Reamostragem Ponderada

Estes métodos usam a mesma idéia de gerar valores de uma distribuição auxiliar porém sem a necessidade de maximização da verossimilhança. A desvantagem é que os valores obtidos são apenas aproximadamente distribuidos segundo a posteriori.

Suponha que temos uma amostra $\theta_1,\dots,\theta_n$ gerada da distribuição auxiliar $q$ e a partir dela construimos os pesos

$$w_i=\frac{p(\theta_i\vert\bfx)/q(\theta_i)} {\sum_{j=1}^n p(\theta_j\vert\bfx)/q(\theta_j)},\quad i=1,\dots,n$$

O método consiste em tomar uma segunda amostra (ou reamostra) de tamanho $m$ da distribuição discreta em $\theta_1,\dots,\theta_n$ com probabilidades $w_1,\dots,w_n$. Aqui também não é necessário que se conheça completamente a posteriori mas apenas o produto priori vezes verossimilhança já que neste caso os pesos não se alteram.

Tomando novamente a priori como densidade auxiliar, i.e. $q(\theta)=p(\theta)$ os pesos se simplificam para

$$w_i=\frac{p(\bfx\vert\theta_i)}{\sum_{j=1}^n p(\bfx\vert\theta_j)},\quad i=1,\dots,n$$

e o algoritmo para geração de valores (aproximadamente) da posteriori então fica

1. gerar valores $\theta_1,\dots,\theta_n$ da distribuição a priori;
2. calcular os pesos $w_i$, $i=1,\dots,n$;
3. reamostrar valores com probabilidades $w_1,\dots,w_n$.

Exercícios

1. Em um modelo de regressão linear simples temos que $y_i\sim N(\beta x_i,1)$. Os dados observados são $\bfy=(-2,0,0,0,2)$ e $\bfx=(-2,-1,0,1,2)$, e usamos uma priori vaga $N(0,4)$ para $\beta$. Faça inferência sobre $\beta$ obtendo uma amostra da posteriori usando reamostragem ponderada. Compare com a estimativa de máxima verossimilhança $\hat{\beta}=0,8$.

2. Para o mesmo modelo do exercício 1 e os mesmos dados suponha agora que a variância é desconhecida, i.e. $y_i\sim N(\beta x_i,\s)$. Usamos uma priori hierárquica para $(\beta,\s)$, i.e. $\beta\vert\s\sim N(0,\s)$ e $\invs\sim G(0,01,0,01)$.
   1. Obtenha uma amostra da posteriori de $(\beta,\s)$ usando reamostragem ponderada.
   2. Baseado nesta amostra, faça um histograma das distribuições marginais de $\beta$ e $\s$.
   3. Estime $\beta$ e $\s$ usando uma aproximação para a média a posteriori. Compare com as estimativas de máxima verossimilhança.

function()
{
# Obtem uma amostra da distribuicao a posteriori dos parametros em um
# modelo de regressao linear simples via metodo de reamostragem
# ponderada (Exercicio 1, desta seção);
x <- c(-2,-1,0,1,2)
y <- c(-2, 0,0,0,2)
n <- 1000; \# tamanho da amostra da priori
m <- 500 ; \# tamanho da reamostra
par(mfrow = c(2, 2))
beta <- matrix(rnorm(n, 0, 2), nrow = n)
l <- matrix(NA, nrow = n)
for(i in 1:n){
l[i] <- exp(- (1/2) * t(y - beta[i] * x) %*% (y - beta[i] * x))
}
p <- matrix(NA, nrow = n)\\
for(i in 1:n)
{
p[i] <- l[i]/sum(l)
}
resample <- sample(beta, size = m, replace = T, prob = p)
hist(beta, col = 0, prob = T, title="priori")
plot(beta, l, title="")
hist(resample, col = 0, prob = T)
list(beta = summary(beta), resample = summary(resample))
}

4.5 Monte Carlo via cadeias de Markov

Em todos os métodos de simulação vistos até agora obtém-se uma amostra da distribuição a posteriori em um único passo. Os valores são gerados de forma independente e não há preocupação com a convergência do algoritmo, bastando que o tamanho da amostra seja suficientemente grande. Por isso estes métodos são chamados não iterativos (não confundir iteração com interação). No entanto, em muitos problemas pode ser bastante difícil, ou mesmo impossível, encontrar uma densidade de importância que seja simultaneamente uma boa aproximação da posteriori e fácil de ser amostrada.

Os métodos de Monte Carlo via cadeias de Markov (MCMC) são uma alternativa aos métodos não iterativos em problemas complexos. A idéia ainda é obter uma amostra da distribuição a posteriori e calcular estimativas amostrais de características desta distribuição. A diferença é que aqui usaremos técnicas de simulação iterativa, baseadas em cadeias de Markov, e assim os valores gerados não serão mais independentes.

Neste capítulo serão apresentados os métodos MCMC mais utilizados, o amostrador de Gibbs e o algoritmo de Metropolis-Hastings. A idéia básica é simular um passeio aleatório no espaço de $\theta$ que converge para uma distribuição estacionária, que é a distribuição de interesse no problema. Uma discussão mais geral sobre o tema pode ser encontrada por exemplo em gam97.

4.5.1 Cadeias de Markov

Uma cadeia de Markov é um processo estocástico $\{X_0,X_1,\dots\}$ tal que a distribuição de $X_t$ dados todos os valores anteriores $X_0,\dots,X_{t-1}$ depende apenas de $X_{t-1}$. Matematicamente,

$$P(X_t\in A\vert X_0,\dots,X_{t-1})=P(X_t\in A\vert X_{t-1})$$

para qualquer subconjunto $A$. Os métodos MCMC requerem ainda que a cadeia seja,

• homogênea, i.e. as probabilidades de trasição de um estado para outro são invariantes;
• irredutível, i.e. cada estado pode ser atingido a partir de qualquer outro em um número finito de iterações;
• aperiódica, i.e. não haja estados absorventes.

e os algoritmos que serão vistos aqui satisfazem a estas condições.

4.5.2 Algoritmo de Metropolis-Hastings

Os algoritmos de Metropolis-Hastings usam a mesma idéia dos métodos de rejeição vistos no capítulo anterior, i.e. um valor é gerado de uma distribuição auxiliar e aceito com uma dada probabilidade. Este mecanismo de correção garante que a convergência da cadeia para a distribuição de equilibrio, que neste caso é a distribuição a posteriori.

Suponha que a cadeia esteja no estado $\theta$ e um valor $\theta'$ é gerado de uma distribuição proposta $q(\cdot\vert\theta)$. Note que a distribuição proposta pode depender do estado atual da cadeia, por exemplo $q(\cdot\vert\theta)$ poderia ser uma distribuição normal centrada em $\theta$. O novo valor $\theta'$ é aceito com probabilidade

$$\alpha(\btheta,\btheta')= \min\left(1,\frac {\pi(\theta')q(\theta \vert\theta')} {\pi(\theta )q(\theta'\vert\theta )}\right).$$ (4.2)

onde $\pi$ é a distribuição de interesse.

Uma característica importante é que só precisamos conhecer $\pi$ parcialmente, i.e. a menos de uma constante já que neste caso a probabilidade (4.2) não se altera. Isto é fundamental em aplicações Bayesianas aonde não conhecemos completamente a posteriori.

Em termos práticos, o algoritmo de Metropolis-Hastings pode ser especificado pelos seguintes passos,

1. Inicialize o contador de iterações $t=0$ e especifique um valor inicial $\theta^{(0)}$.
2. Gere um novo valor $\theta'$ da distribuição $q(\cdot\vert\theta)$.
3. Calcule a probabilidade de aceitação $\alpha(\theta,\theta')$ e gere $u\sim U(0,1)$.
4. Se $u\le\alpha$ então aceite o novo valor e faça $\theta^{(t+1)}=\theta'$, caso contrário rejeite e faça $\theta^{(t+1)}=\theta$.
5. Incremente o contador de $t$ para $t+1$ e volte ao passo 2.

Uma useful feature of the algorithm is that the target distribution needs only be known up to a constant of proportionality since only the target ratio $\pi(\btheta')/\pi(\btheta)$ is used in the acceptance probability. Note also that the chain may remain in the same state for many iterations and in practice a useful monitoring device is given by the average percentage of iterations for which moves are accepted. Hastings (1970) suggests that this acceptance rate should always be computed in practical applications.

The independence sampler is a Metropolis-Hastings algorithm whose proposal distribution does not depend on the current state of the chain, i.e., $q(\btheta,\btheta')=q(\btheta')$. In general, $q(\cdot)$ should be a good approximation of $\pi(\cdot)$, but it is safest if $q(\cdot)$ is heavier-tailed than $\pi(\cdot)$.

The Metropolis algorithm considers only symmetric proposals, i.e., $q(\btheta,\btheta')=q(\btheta',\btheta)$ for all values of $\btheta$ and $\btheta'$, and the acceptance probability reduces to

$$\alpha(\btheta,\btheta')=\min\left(1,\frac{\pi(\btheta')}{\pi(\btheta)}\right).$$

A special important case is the random-walk Metropolis for which $q(\btheta,\btheta')=q(\vert\btheta-\btheta'\vert)$, so that the probability of generating a move from $\btheta$ to $\btheta'$ depends only on the distance between them. Using a proposal distribution with variance $\s$, very small values of $\s$ will lead to small jumps which are almost all accepted but it will be difficult to traverse the whole parameter space and it will take many iterations to converge. On the other hand, large values of $\s$ will lead to an excessively high rejection rate since the proposed values are likely to fall in the tails of the posterior distribution.

Typically, there will be an optimal value for the proposal scale $\sigma$ determined on the basis of a few pilot runs which lies in between these two extremes (see for example, Roberts, Gelman and Gilks, 1997). We return to this point later and, in particular, discuss an approach for choosing optimal values for the parameters of the proposal distribution for (RJ)MCMC algorithms in Chapter 6.

4.5.3 Amostrador de Gibbs

while in the Gibbs sampler the chain will always move to a new value. Gibbs sampling is an MCMC scheme where the transition kernel is formed by the full conditional distributions, $\pi(\theta_i\vert\btheta_{-i})$, where $\btheta_{-i}=(\theta_1,\dots,\theta_{i-1},\theta_{i+1},\dots,\theta_d)'$. In general, each one of the components $\theta_i$ can be either uni- or multi-dimensional. So, the full conditional distribution is the distribution of the $i$th component of $\btheta$ conditioning on all the remaining components, and it is derived from the joint distribution as follows,

$$\pi(\theta_i\vert\btheta_{-i})=\frac{\pi(\btheta)}{\dint\pi(\btheta)d\theta_i}.$$

If generation schemes to draw a sample directly from $\pi(\btheta)$ are costly, complicated or simply unavailable but the full conditional distributions are completely known and can be sampled from, then Gibbs sampling proceeds as follows,

1. Initialize the iteration counter of the chain $t=1$ and set initial values
   $\btheta^{(0)}=(\theta_1^{(0)},\dots,\theta_d^{(0)})'$.
2. Obtain a new value of $\btheta^{(t)}$ from $\btheta^{(t-1)}$ through successive generation of values
   

   $$\theta_1^{(t)} \sim \pi(\theta_1\vert\theta_2^{(t-1)},\theta_3^{(t-1)},\dots,\theta_d^{(t-1)})$$

   $$\theta_2^{(t)} \sim \pi(\theta_2\vert\theta_1^{(t)},\theta_3^{(t-1)},\dots,\theta_d^{(t-1)})$$

   $$\vdots$$

   $$\theta_d^{(t)} \sim \pi(\theta_d\vert\theta_1^{(t)},\theta_2^{(t)},\dots,\theta_{d-1}^{(t)})$$

   
   

3. Increment the counter $t$ to $t+1$ and return to step 2 until convergence is reached.

So, each iteration is completed after $d$ moves along the coordinates axes of the components of $\btheta$. When convergence is reached, the resulting value $\btheta$ is a draw from $\pi(\btheta)$. It is worth noting that, even in a high-dimensional problem, all of the simulations may be univariate, which is usually a computational advantage.

However, the Gibbs sampler does not apply to problems where the number of parameters varies because of the lack of irreducibility of the resulting chain. When the length of $\btheta$ is not fixed and its elements need not have a fixed interpretation across all models, to resample some components conditional on the remainder would rarely be meaninful.

Note also that the Gibbs sampler is a special case of the Metropolis-Hastings algorithm, in which individual elements of $\btheta$ are updated one at a time (or in blocks), with the full conditional distribution as the candidate generating function and acceptance probabilities uniformly equal to 1.

4.5.4 Updating strategies

In the above scheme, all the components of $\btheta$ are updated in the same deterministic order at every iteration. However, other scanning or updating strategies are possible for visiting the components of $\btheta$. Geman and Geman (1984) showed in a discrete setting that any updating scheme that guarantees that all components are visited infinitely often when the chain is run indefinitely, converges to the joint distribution of interest, i.e., $\pi(\btheta)$. For example, Zeger and Karim (1991) describe a Gibbs sampling scheme where some components are visited only every $k$th iteration, which still guarantees that every component is updated infinitely often for finite, fixed $k$.

Roberts and Sahu (1997) consider a random permutation scan where at each iteration a permutation of $\{1,\dots,d\}$ is chosen and components are visited in that order. In particular, they showed that when $\pi$ is multivariate normal, convergence for the deterministic scan is faster than for the random scan if the precision matrix is tridiagonal ($\theta_i$ depends only on $\theta_{i-1}$ and $\theta_{i+1}$) or if it has non-negative partial correlations.

4.5.5 Blocking

In principle, the way the components of the parameter vector $\btheta$ are arranged in blocks of parameters is completely arbitrary and includes blocks formed by scalar components as special cases. However, the structure of the Gibbs sampler imposes moves according to the coordinate axes of the blocks, so that larger blocks allow moves in more general directions. This can be very beneficial, although more computationally demanding, in a context where there is high correlation between individual components since these higher dimensional moves incorporate information about this dependence. Parameter values are then generated from the joint full conditional distribution for the block of parameters considered.

Roberts and Sahu (1997) showed that for a multivariate normal $\pi$ and random scans, convergence improves as the number of blocks decreases. They also proved that blocking can hasten convergence for non-negative partial correlation distributions and even more as the partial correlation of the components in the block gets larger. However, they also provided an example where blocking worsens convergence.

4.5.6 Completion

Even when every full conditional distribution associated with the target distribution $\pi$ is not explicit there can be a density $\pi^*$ for which $\pi$ is a marginal density, i.e.,

$$\int \pi^*(\btheta,\bfz) d\bfz = \pi(\btheta)$$

and such that all the full conditionals associated with $\pi^*$ are easy to simulate from. Then the Gibbs sampler can be implemented in $\pi^*$ instead of $\pi$ and this is called the completion Gibbs sampler because $\pi^*$ is a completion of $\pi$. The required sample from the target distribution is obtained by marginalizing again, i.e., integrating $\bfz$ out.

This approach was actually one of the first appearances of the Gibbs sampler in Statistics with the introduction of data augmentation by Tanner and Wong (1987). It is also worth noting that, in principle, this Gibbs sampler does not require that the completion of $\pi$ into $\pi^*$ and of $\btheta$ into $(\btheta,\bfz)$ should be related to the problem of interest and the vector $\bfz$ might have no meaning from a statistical point of view.

4.5.7 The Slice Sampler

This is a very general version of the Gibbs sampler which applies to most distributions and is based on the simulation of specific uniform random variables. In its simplest version when only one variable is being updated, if $\pi$ can be written as a product of functions, i.e.,

$$\pi(\theta)=\prod_{i=1}^k f_i(\theta),$$

where $f_i$ are positive functions but not necessarily densities then $f$ can be completed (or demarginalised) into

$$\prod_{i=1}^k I_{0<z_i<f_i(\theta)}.$$

The slice sampler consists of generating $(z_1,\dots,z_k,\theta)$ from their full conditional distributions, i.e.,

• generate $z_i^{(t+1)}$ from $U[0,f_i(\theta^{(t)})]$, $i=1,\dots,k$ and
• generate $\theta^{(t+1)}$ from a uniform distribution in $A^{(t+1)}=\{y:f_i(y)>z_i^{(t+1)}, i=1,\dots,k\}$.

Roberts and Rosenthal (1998) study the slice sampler and show that it usually enjoys good theoretical properties. In practice there may be problems as $d$ increases since the determination of the set $A^{(t+1)}$ may get increasingly complex.

Further details about the Gibbs sampler and related algorithms are given, for example, in Gamerman (1997, Chapter 5) and Robert and Casella (1999, Chapter 7).

4.6 Posterior Model Probabilities

The posterior model probability is obtained as

$$p(k\vert\bfy)=\frac{p(\bfy\vert k)p(k)}{p(\bfy)}$$

where the term $p(\bfy\vert k)$ is sometimes referred to as the marginal likelihood for model $k$ and is calculated as

$$p(\bfy\vert k)=\int p(\bfy\vert\btheta^{(k)},k)p(\btheta^{(k)}\vert k)d\btheta^{(k)}.$$

Also, $1/p(\bfy\vert k)$ is the normalisation constant for $p(\btheta^{(k)}\vert k,\bfy)$, the posterior density of $\btheta$ within model $k$.

Hence, the posterior probability of a certain model is proportional to the product of the prior probability and the marginal likelihood for that model. It is also worth noting that, in practice, $p(\bfy)$ is unknown so that typically the model probabilities are known only up to a normalisation constant.

The above integral is commonly analytically intractable but may be approximated in a number of ways by observing that it can be regarded as the expected value of the likelihood with respect to the prior distribution $p(\btheta^{(k)}\vert k)$. In terms of simulation techniques, the simplest estimate consists of simulating $n$ values $\btheta_1,\dots,\btheta_n$ from the prior, evaluating the likelihood at those values and computing the Monte Carlo estimate

$$\hat{p}(\bfy\vert k)=\frac{1}{n}\sum_{i=1}^n p(\bfy\vert\btheta_i^{(k)},k).$$

This estimator has high variance with possibly few terms contributing substantially to the sum in cases of disagreement between prior and likelihood. Various alternative estimators are reviewed in Gamerman (1997, Chapter 7) and analytical approximations supported by asymptotic normal theory might also be used. Other alternatives will be explored in the next section.

Having obtained the posterior model probabilities these may be used for either selecting the model with the highest probability or highest Bayes factor from the list of candidate models (model selection), or estimating some quantity under each model and then averaging the estimates according to how likely each model is, that is, using these probabilities as weights (model averaging). In the next section, we present MCMC methods that take into account different models simultaneously.