6. Algoritmo de Metropolis-Hastings

Este capítulo é o primeiro de uma série de dois sobre métodos de simulação baseados em Cadeias de Markov. Embora o algoritmo de Metropolis-Hastings possa ser visto como um dos algoritmos de Cadeias de Markov Monte Carlo (MCMC) mais gerais, também é um dos mais simples de entender e explicar, tornando-o um algoritmo ideal para começar.

Este capítulo começa com uma rápida atualização das Cadeias de Markov, apenas o básico necessário para entender os algoritmos. Em seguida, definimos o algoritmo Metropolis-Hastings, focando nas versões mais comuns do algoritmo. Acabamos discutindo a calibração do algoritmo via sua taxa de aceitação na Seção 6.5.



6.1 Introdução

Por motivos que ficarão mais claros à medida que prosseguirmos, agora faremos uma mudança fundamental na escolha de nossa estratégia de simulação. Até agora, normalmente geramos variáveis iid diretamente da densidade de interesse \(f\) ou indiretamente no caso de amostragem por importância. Em vez disso, o algoritmo Metropolis-Hastings apresentado abaixo gera variáveis correlacionadas de uma cadeia de Markov.

A razão pela qual optamos por uma mudança tão radical é que as Cadeias de Markov carregam diferentes propriedades de convergência que podem ser exploradas para fornecer propostas mais fáceis nos casos em que a amostragem de importância genérica não se aplica prontamente. Por um lado, os requisitos do alvo \(f\) são mínimos, o que permite configurações onde muito pouco se sabe sobre \(f\). Outra razão, conforme ilustrado no próximo capítulo, é que essa perspectiva de Markov leva a decomposições eficientes de problemas de alta dimensão em uma sequência de problemas menores que são muito mais fáceis de resolver.

Portanto, esteja avisado de que este é um capítulo crucial, pois agora apresentamos uma perspectiva totalmente nova sobre a geração de variáveis aleatórias, que teve um efeito profundo na pesquisa e expandiu a aplicação de métodos estatísticos para resolver problemas mais difíceis e problemas mais relevantes nos últimos vinte anos, embora as origens dessas técnicas estejam ligadas às do método de Monte Carlo no remoto centro de pesquisa de Los Alamos durante a Segunda Guerra Mundial. No entanto, apesar do recurso aos princípios da Cadeia de Markov que são brevemente detalhados na próxima seção, a implementação desses novos métodos não é mais difícil do que os capítulos anteriores, e não há necessidade de se aprofundar mais na teoria das Cadeias de Markov, como você logo descubrirá.

A maior parte do seu tempo e energia serão gastos em projetar e avaliar seus algoritmos MCMC, assim como nos capítulos anteriores, não em estabelecer teoremas de convergência, então vamos com calma!


6.2 Uma olhada na teoria da Cadeias de Markov

Esta seção pretende ser uma atualização minimalista das Cadeias de Markov para definir o vocabulário das Cadeias de Markov, nada mais. Caso você tenha dúvidas ou queira mais detalhes sobre essas noções, é altamente recomendável verificar um tratamento mais completo, como em Cadeias de Markov, Robert and Casella (2004) ou Meyn and Tweedie (1993); pois nenhuma teoria de convergência é fornecida no presente livro.

Uma Cadeia de Markov \(\{X^{(t)}\}\) é uma sequência de variáveis aleatórias dependentes \[ X^{(0)},X^{(1)},X^{(2)},\cdots,X^{(t)},\dots \] tal que a distribuição de probabilidade de \(\{X^{(t)}\}\), dadas as variáveis passadas, depende apenas de \(\{X^{(t-1)}\}\). Essa distribuição de probabilidade condicional é chamada de kernel de transição ou kernel de Markov \(K\); isto é, \[ X^{(t+1)} \, | \, X^{(0)},X^{(1)},X^{(2)},\cdots,X^{(t)} \sim K(X^{(t)},X^{(t+1)})\cdot \]

Por exemplo, uma Cadeia de Markov de caminhada aleatória simples satisfaz \[ X^{(t+1)} = X^{(t)} + \epsilon_t, \] onde \(\epsilon_t\sim N(0,1)\), independentemente de \(X^{(t)}\); portanto, o kernel de Markov \(K(X^{(t)},X^{(t+1)})\) corresponde a uma densidade \(N(X^{(t)},1)\).

Na maioria das vezes, as Cadeias de Markov encontradas nas configurações de Cadeia de Markov Monte Carlo (MCMC) desfrutam de uma propriedade de estabilidade muito forte. De fato, existe uma distribuição de probabilidade estacionária por construção para essas cadeias; isto é, existe uma distribuição de probabilidade \(f\) tal que se \(X^{(t)}\sim f\), então \(X^{(t+1)}\sim f\).

Portanto, formalmente, o kernel e a distribuição estacionária satisfazem a equação \[\begin{equation} \tag{6.1} \int_\mathcal{X} K(x,y)f(x)\mbox{d}x = f(y)\cdot \end{equation}\]

A existência de uma distribuição estacionária ou estacionariedade, impõe uma restrição preliminar em \(K\) chamada irredutibilidade na teoria das Cadeias de Markov, que é que o kernel \(K\) permite movimentos livres em todo o espaço de estados, ou seja, não importa o valor inicial \(X^{(0)}\), a sequência \(\{X^{(t)}\}\) tem uma probabilidade positiva de eventualmente atingir qualquer região do espaço de estados.

Uma condição suficiente é que \(K(x,\cdot)>0\) em todos os lugares. A existência de uma distribuição estacionária tem grandes consequências no comportamento da cadeia \(\{X^{(t)}\}\), sendo uma delas que a maioria das cadeias envolvidas nos algoritmos MCMC são recorrentes, ou seja, retornarão a qualquer conjunto arbitrário não desprezível um número infinito de vezes.

No caso de cadeias recorrentes, a distribuição estacionária também é uma distribuição limite no sentido de que a distribuição limite de \(X^{(t)}\) é \(f\) para quase qualquer valor inicial \(X^{(0)}\). Essa propriedade também é chamada de ergodicidade e obviamente tem grandes consequências do ponto de vista da simulação, pois, se um determinado kernel \(K\) produz uma Cadeia de Markov ergódica com distribuição estacionária \(f\), gerar uma cadeia a partir desse kernel \(K\) acabará por produzir simulações a partir de \(f\).

Em particular, para funções integráveis \(h\), a média padrão \[\begin{equation} \tag{6.2} \dfrac{1}{T}\sum_{t=1}^T h(X^{(t)})-\mbox{E}_f\big(h(X)\big), \end{equation}\] o que significa que a Lei dos Grandes Números, que está na base dos métodos de Monte Carlo (Seção 3.2), também pode ser aplicada em configurações MCMC. Às vezes é chamado de Teorema Ergódico.

Não vamos nos aprofundar mais na teoria da convergência de algoritmos MCMC, confiando na garantia de que as versões padrão desses algoritmos, como o algoritmo Metropolis-Hastings ou o amostrador de Gibbs, são quase sempre teoricamente convergentes. De fato, o problema real com os algoritmos MCMC é que, apesar dessas garantias de convergência, a implementação prática desses princípios pode implicar em um tempo de convergência muito longo ou, pior ainda, pode dar uma impressão de convergência faltando alguns aspectos importantes de \(f\), como discutido em Capítulo 8.

Existe, no entanto, um caso em que a convergência nunca ocorre, nomeadamente quando, num ambiente bayesiano, a distribuição a posteriori não é adequada (Robert, 2001) uma vez que a cadeia não pode ser recorrente. Com o uso de a prioris impróprios \(f(x)\) sendo bastante comum em modelos complexos, existe a possibilidade de que o produto \(\mbox{verossimilhança}\times \mbox{a priori}\), \(\ell(x)\times f(x)\), não seja integrável e que esse problema passe despercebido por causa da complexidade inerente. Em tais casos, as Cadeias de Markov podem ser simuladas em conjunto com o alvo \(\ell(x)\times f(x)\), mas não podem convergir.

Nos melhores casos, as Cadeias de Markov resultantes exibirão rapidamente um comportamento divergente, o que indica que há um problema. Infelizmente, nos piores casos, essas Cadeias de Markov apresentam todos os sinais externos de estabilidade e, portanto, falham em indicar a dificuldade. Mais detalhes sobre este assunto são discutidos na Seção 7.6.4 do próximo capítulo.


6.3 Algoritmo Metropolis-Hasting básico

O princípio de funcionamento dos métodos de Monte Carlo de Cadeias de Markov são bastantes simples de descrever. Dada uma densidade alvo \(f\), construímos um kernel de Markov \(K\) com distribuição estacionária \(f\) e então geramos uma Cadeia de Markov \(\{X^{(t)}\}\) usando este kernel para que a distribuição limite de \(\{X^{(t)}\}\) seja \(f\) e integrais possam ser aproximados de acordo com o Teorema Ergódico (6.2).

A dificuldade deve, portanto, estar em construir um kernel \(K\) associado a uma densidade \(f\) arbitrária. Mas, milagrosamente, existem métodos para derivar tais kerneis (núcleos) que são universais no sentido de que são teoricamente válidos para qualquer densidade \(f\).

O algoritmo Metropolis-Hastings é um exemplo desses métodos. O amostrador de Gibbs, descrita no Capítulo 7, é outro exemplo com potencial igualmente universal. Dada a densidade alvo \(f\), ela está associada a uma densidade condicional funcional \(q(y|x)\) que, na prática, é fácil de simular.

Além disso, \(q\) pode ser quase arbitrário, pois os únicos requisitos teóricos são que a razão \(f(y)/q(y|x)\) seja conhecida até uma constante independente de \(x\) e que \(q(\cdot|x)\) tenha dispersão suficiente para levar a uma exploração de todo o suporte de \(f\). Mais uma vez, enfatizamos a incrível característica do algoritmo Metropolis-Hastings de que, para cada \(q\) dado, podemos então construir um kernel Metropolis-Hastings tal que \(f\) seja sua distribuição estacionária.


6.3.1 Um algoritmo genérico de Cadeias de Markov Monte Carlo (MCMC)

O algoritmo Metropolis-Hastings associado com a densidade objetiva (alvo) \(f\) e a densidade condicional \(q\) produz uma Cadeia de Markov \(\{X^{(t)}\}\) através do seguinte kernel de transição:


Algoritmo 4 Metropolis-Hasting
Dado \(x^{(t)}\),

  1. Gerar, \(Y_t\sim q(y \, | x^{(t)})\);

  2. Aceitar \(X^{(t+1)} = \left\{ \begin{array}{rcl} Y_t & \mbox{ com probabilidade} & \rho(x^{(t)},Y_t), \\ x^{(t)} & \mbox{ com probabilidade} & 1-\rho(x^{(t)},Y_t),\end{array}\right.\),
    onde \(\rho(x,y)=\min\left\{ \dfrac{f(y)}{f(x)}\dfrac{q(x|y)}{q(y|x)},1\right\}\).


Uma implementação R genérica é direta, assumindo um gerador para \(q(y|x)\) está disponível como geneq(x). Se x[t] denota o valor de \(\{X^{(t)}\}\),


y=geneq(x[t])
if (runif(1)<f(y)*q(y,x[t])/(f(x[t])*q(x[t],y)))
  { x[t+1]=y }
  else 
  { x[t+1]=x[t] }
já que o valor \(y\) é sempre aceito quando a razão é maior que um.

A distribuição \(q\) é chamada de distribuição instrumental ou proposta ou candidata e a probabilidade \(\rho(x,y)\) de probabilidade de aceitação Metropolis-Hastings. Deve ser distinguido da taxa de aceitação, que é a média da probabilidade de aceitação ao longo das iterações, \[ \overline{\rho}=\lim_{T\to\infty} \dfrac{1}{T}\sum_{t=0}^T \rho(X^{(t)},Y_t)=\int \rho(x,y)f(x)q(y|x)\mbox{d}y\mbox{d}x\cdot \]

Essa quantidade permite avaliar o desempenho do algoritmo, conforme discutido na Seção 6.5.

Embora, à primeira vista, o Algoritmo 4 não pareça diferir do Algoritmo 2, exceto pela notação, existem duas diferenças fundamentais entre os dois algoritmos. A primeira diferença está na sua utilização, já que o Algoritmo 2 visa maximizar uma função \(h(x)\), enquanto o Algoritmo 4 visa explorar o suporte da densidade \(f\) de acordo com sua probabilidade.

A segunda diferença está em suas propriedades de convergência. Com a escolha adequada de um cronograma de temperatura \(T_t\) no Algoritmo 2, o algoritmo de anelamento simulado converge para os máximos da função \(h\), enquanto o algoritmo Metropolis-Hastings está convergindo para a própria distribuição \(f\). Por fim, modificar a proposta \(q\) ao longo das iterações pode ter consequências drásticas no padrão de convergência desse algoritmo, conforme discutido na Seção 8.5.

O Algoritmo 4 satisfaz a chamada condição de balanço detalhado, \[\begin{equation} \tag{6.3} f(x)K(y|x) = f(y)K(x|y), \end{equation}\] da qual podemos deduzir que \(f\) é a distribuição estacionária da cadeia \(\{X^{(t)}\}\) integrando cada lado da igualdade em \(x\), ver Exercício 6.8.

O fato de o Algoritmo 4 estar naturalmente associado a \(f\) como sua distribuição estacionária vem facilmente como consequência da condição de equilíbrio detalhada para uma escolha arbitrária do par \((f,q)\). Na prática, o desempenho do algoritmo obviamente dependerá fortemente dessa escolha de \(q\), mas considere primeiro um exemplo direto onde o Algoritmo 4 pode ser comparado com a amostragem iid.


Exemplo 6.1

Lembre-se do Exemplo 2.7, onde usamos um algoritmo Accept-Reject para simular uma distribuição beta. Podemos também usar um algoritmo Metropolis-Hastings, onde a densidade alvo \(f\) é a densidade \(Beta(2.7, 6.3)\) e o candidato \(q\) é uniforme em (0,1), o que significa que não depende do valor anterior da cadeia.

Uma amostra Metropolis-Hastings é então gerada com o seguinte código R:

set.seed(6578)
a=2.7; b=6.3 # initial values
Nsim=5000
X=rep(runif(1),Nsim) # initialize the chain
for (i in 2:Nsim){
  Y=runif(1)
  rho=dbeta(Y,a,b)/dbeta(X[i-1],a,b)
  X[i]=X[i-1] + (Y-X[i-1])*(runif(1)<rho)
}
par(mar=c(4,4,1,1))
plot(X,xlab="Iterations",type="l",xlim=c(4500,4800))
grid()

Figura. 6.1. Sequência \(\{X^{(t)}\}\) para \(t = 4500,\cdots,4800\), quando simulada a partir do algoritmo Metropolis-Hastings com proposta uniforme e alvo \(Beta(2.7,6.3)\).

Uma representação da sequência \(\{X^{(t)}\}\) por plot não produz nenhum padrão na simulação, pois a cadeia explora o mesmo intervalo em diferentes períodos. Se ampliarmos o período final, para \(4500\leq t\leq 4800\), a Figura 6.1 exibe algumas características das sequências Metropolis-Hastings, ou seja, para alguns intervalos de tempo, a sequência \(\{X^{(t)}\}\) não muda porque todos \(Y_t\)’s correspondentes são rejeitados.

Observe que essas múltiplas ocorrências do mesmo valor numérico devem ser mantidos na amostra como tal; caso contrário, a validade da aproximação de \(f\) é perdida. De fato, ao considerar toda a cadeia como uma amostra, seu histograma se aproxima adequadamente do alvo \(Beta(2.7,6.3)\).

A Figura 6.2 mostra histogramas e densidades sobrepostas para esta amostra Metropolis-Hastings e para uma amostra iid (exata) desenhada usando o comando rbeta. Os primeiros são bastante semelhantes, e isso pode ser verificado ainda mais usando um teste Kolmogorov-Smirnov de igualdade entre as duas amostras: 

mean(X); var(X)
## [1] 0.3005143
## [1] 0.02083906
ks.test(jitter(X),rbeta(5000,a,b))
## 
##  Asymptotic two-sample Kolmogorov-Smirnov test
## 
## data:  jitter(X) and rbeta(5000, a, b)
## D = 0.0162, p-value = 0.528
## alternative hypothesis: two-sided

que afirma que ambas as amostras são compatíveis com a mesma distribuição.

Uma verificação adicional, se leve, de concordância é fornecida pelos momentos. Por exemplo, como a média e a variância de uma distribuição \(Beta(a,b)\) são \(a/(a+b)\) e \(ab/(a+b)^2(a+b+1)\), respectivamente, podemos comparar \[ \overline{X}=0.3005, \qquad S^2=0.0208, \] com os valores teóricos de

a/(a+b)
## [1] 0.3

para a média e

a*b/((a+b+1)*(a+b)^2)
## [1] 0.021
para a variância. 

par(mfrow=c(1,2),mar=c(4,4,1,1))
hist(X,freq=FALSE,main="",ylab="",ylim=c(0,3),xlab="")
curve(dbeta(x,a,b),col="magenta",add=TRUE,lwd=2)
box();grid()
hist(rbeta(Nsim,a,b),freq=FALSE,main="",ylab="",ylim=c(0,3),xlab="")
curve(dbeta(x,a,b),col="magenta",add=TRUE,lwd=2)
box();grid()

Figura. 6.2. Histogramas de variáveis aleatórias beta \(Beta(2.7,6.3)\) com funções de densidade sobrepostas. No painel esquerdo, as variáveis foram geradas a partir de um algoritmo Metropolis-Hastings com um candidato uniforme, e no painel direito as variáveis aleatórias foram geradas diretamente usando rbeta(n,2.7,6.3).


Embora o MCMC e os resultados da amostragem exata pareçam idênticos na Figura 6.2, é importante lembrar que a amostra Monte Carlo de Cadeia de Markov (MCMC) tem correlação, enquanto a amostra iid não. Isso significa que a qualidade da amostra está necessariamente degradada ou, em outras palavras, que precisamos de mais simulações para obter a mesma precisão. Esta questão é formalizada através da noção de tamanho amostral efetivo para cadeias de Markov (Seção 8.4.3).

No caso simétrico, ou seja, quando \(q(x|y) = q(y|x)\), a probabilidade de aceitação \(\rho(x_t,y_t)\) é determinada pela razão objetiva \(f(y_t)/f(x^{(t)})\) e, portanto, mesmo a probabilidade de aceitação é independente de \(q\). Este caso especial é detalhado na Seção 6.4.1. Mais uma vez, os algoritmos de Metropolis-Hastings compartilham a mesma característica do Algoritmo 2 de otimização estocástica, consulte a Seção 5.5, ou seja, eles sempre aceitam valores de \(y_t\) tais que a razão \(f(y_t)/q(y_t|x^{(t)})\) é aumentado em comparação com o valor anterior \(f(x^{(t)})/q(x^{(t)}|y_t)\).

Alguns valores \(y_t\) tais que a razão é diminuída também podem ser aceitos, dependendo da razão das razões, mas se a diminuição for muito acentuada, o valor \(y_t\) proposto quase sempre será rejeitado. Esta propriedade indica como a escolha de \(q\) pode impactar o desempenho do algoritmo Metropolis-Hastings. Se o domínio explorado por \(q\), seu suporte, for muito pequeno, comparado com o alcance de \(f\), a cadeia de Markov terá dificuldades em explorar este alcance e assim irá convergir muito lentamente, se for o caso para propósitos práticos.

Outra propriedade interessante do algoritmo Metropolis-Hastings que aumenta seu apelo é que ele depende apenas das proporções \[ f(y_t)/f(x^{(t)}) \qquad \mbox{e} \qquad q(x^{(t)}\, | \, y_t)/q(y_t \, | \, x^{(t)})\cdot \] É, portanto, independente das constantes de normalização. Além disso, como tudo o que importa é a capacidade de (a) simular a partir de \(q\) e (b) calcular a razão \(f(y_t)/q(y_t|x^{(t)})\), \(q\) pode ser escolhido de tal forma que as partes intratáveis de \(f\) são eliminados na proporção.

Como \(q(y|x)\) é uma densidade condicional, a integral é um em \(y\) e, como tal, envolve um termo funcional que depende de \(y\) e \(x\), bem como um termo de normalização que depende de \(x\), ou seja, \(q(y|x) = C(x)\widetilde{q}(x,y)\).

Ao observar acima que a probabilidade de aceitação Metropolis-Hastings não depende de constantes de normalização, termos como \(C(x)\) são obviamente excluídos desta observação, pois devem aparecer na probabilidade de aceitação, para não prejudicar a distribuição estacionária da cadeia.


6.3.2 O algoritmo Metropolis-Hastings independente

O algoritmo Metropolis-Hastings da Seção 6.3.1 permite uma distribuição candidata \(q\) que depende apenas do estado atual da cadeia. Se agora exigirmos que o candidato \(q\) seja independente deste estado atual da cadeia, ou seja, \(q(y|x) = g(y)\), obteremos um caso especial do algoritmo original:


Algoritmo 5 Metropolis-Hasting independente
Dado \(x^{(t)}\),

  1. Gerar, \(Y_t\sim q(y \, | x^{(t)})\);

  2. Aceitar \(X^{(t+1)} = \left\{ \begin{array}{rl} Y_t & \mbox{ com probabilidade } \min\left\{ \dfrac{f(Y_t)}{f(x^{(t)}}\dfrac{g(x^{(t)})}{g(Y_t)},1\right\}, \\ x^{(t)} & \mbox{ caso contrário} \end{array}\right.\),


Este método então aparece como uma generalização direta do método Accept-Reject no sentido de que a distribuição instrumental tem a mesma densidade \(g\) do método Accept-Reject. Assim, os valores \(Y_t\) propostos são os mesmos, senão os aceitos.

Os algoritmos Metropolis-Hastings e os métodos Accept-Reject (Seção 2.3), ambos métodos de simulação genéricos, têm semelhanças entre eles que permitem a comparação, embora seja bastante raro considerar o uso de uma solução Metropolis-Hastings quando um algoritmo Accept-Reject está disponível. Em particular, considere o seguinte:

  1. A amostra Accept–Reject é iid, enquanto a amostra Metropolis–Hastings não é. Embora os \(Y_t\) s sejam gerados independentemente, a amostra resultante não é iid, até porque a probabilidade de aceitação de \(Y_t\) depende de \(X^{(t)}\), exceto no caso trivial em que \(f = g\).

  2. A amostra Metropolis–Hastings envolverá ocorrências repetidas do mesmo valor, pois a rejeição de \(Y_t\) leva à repetição de \(X^{(t)}\) no tempo \(t+1\). Isso terá impacto em testes como ks.test que não aceitam empates.

  3. A etapa de aceitação em Accept–Reject requer o cálculo do limite superior \(M\geq \sup_x f(x)/g(x)\), que não é exigido pelo algoritmo Metropolis-Hastings. Esta é uma característica atraente de Metropolis-Hastings, se calcular \(M\) for demorado ou se o \(M\) existente for impreciso e, portanto, induzir um desperdício de simulações.

O Exercício 6.4 fornece uma primeira comparação de Metropolis–Hastings com um algoritmo Accept–Reject já usado no Exercício 2.20 quando ambos os algoritmos são baseados no mesmo candidato.

A Figura 6.3 ilustra o Exercício 6.4 comparando as amostras Accept-Reject e Metropolis-Hastings. Nesse cenário, operacionalmente, o algoritmo Metropolis-Hastings independente tem desempenho muito semelhante ao algoritmo Accept-Reject, que de fato gera variáveis aleatórias perfeitas e independentes. 

a=4.85;nsim=10000;
X1=X2=array(0,dim=c(nsim,1))                  #AR & MH
X1[1]=X2[1]=rgamma(1,a,rate=1)                #initialize the chain
for (i in 2:nsim){
    Y=rgamma(1,floor(a),rate=floor(a)/a)      #candidate
    rhoAR=(exp(1)*Y*exp(-Y/a)/a)^(a-floor(a))
    rhoMH=(dgamma(Y,a,rate=1)/dgamma(X2[i-1],a,rate=1))/(dgamma(Y,floor(a),
    rate=floor(a)/a)/dgamma(X2[i-1],floor(a),rate=floor(a)/a))
    rhoMH=min(rhoMH,1)
    X1[i]=Y*(runif(1)<rhoAR)                  #accepted values
    X2[i]=X2[i-1] + (Y-X2[i-1])*(runif(1)<rhoMH)
}
X1=X1[X1!=0]                                  #The AR sample
mean(X1);var(X1);mean(X2[2500:nsim]);var(X2[2500:nsim])
## [1] 4.865651
## [1] 4.925934
## [1] 4.849637
## [1] 4.91972
par(mfrow=c(2,2),mar=c(4,4,2,2))
hist(X1,col="grey",nclas=125,freq=FALSE,xlab="",main="Accept-Reject",xlim=c(0,15))
curve(dgamma(x, a, rate=1),lwd=2,add=TRUE)
hist(X2[2500:nsim],nclas=125,col="grey",freq=FALSE,xlab="",main="Metropolis-Hastings",xlim=c(0,15))
curve(dgamma(x, a, rate=1),lwd=2,add=TRUE)
acf(X1,lag.max=50,lwd=2,col="red",ylim=c(-0.06,0.06))              #Accept-Reject
acf(X2[2500:nsim],lag.max=50,lwd=2,col="blue",ylim=c(-0.06,0.06))  #Metropolis-Hastings

Figura. 6.3. Histogramas e funções de autocovariância de um algoritmo gama Accept-Reject (painéis à esquerda) e um algoritmo gama Metropolis-Hastings (painéis à direita). O alvo é uma distribuição \(Gama(4.85,1)\) e o candidato é uma distribuição \(Gama(4,4/4.85)\). A função de autocovariância é calculada com a função R acf.

Teoricamente, também é possível usar um par \((f,g)\) tal que não exista um limite \(M\) em \(f/g\) e, portanto, usar Metropolis-Hastings quando Accept-Reject não é possível. No entanto, conforme detalhado em Robert and Casella (2004) e ilustrado no exemplo formal a seguir, o desempenho do algoritmo Metropolis-Hastings é muito ruim, embora seja muito forte desde que \(\sup f/g = M < \infty\).


Exemplo 6.2

Para gerar uma variável aleatória Cauchy, isto é, quando \(f\) corresponde a uma densidade \(Cauchy(0,1)\), formalmente é possível usar um candidato \(N(0,1)\) dentro de um algoritmo Metropolis-Hastings. O seguinte código R fará isso: 

set.seed(105)
Nsim=10^4
X=c(rt(1,1)) # initialize the chain from the stationary
for (t in 2:Nsim){
 Y=rnorm(1)  # candidate normal
 rho=dt(Y,1)*dnorm(X[t-1])/(dt(X[t-1],1)*dnorm(Y))
 X[t]=X[t-1] + (Y-X[t-1])*(runif(1)<rho)
 }

Ao executar este código, às vezes você pode começar com um valor grande para \(X^{(0)}\), digamos 12.788. Nesse caso, dnorm(X[t-1]) é igual a 0 porque, embora 12.788 possa ser formalmente uma realização de uma \(N(0,1)\), ele induz problemas computacionais sob fluxo: 

pnorm(12.78,log=T,low=F)/log(10)
## [1] -36.97455

Fig. 5.4. Comparação de uma maximização numérica e estocástica de uma verossimilhança de Cauchy em termos do tamanho da amostra via (topo) respectivas localizações das avaliações numéricas e estocásticas dos argumentos, plotadas ao longo da diagonal; (abaixo) erro relativo da avaliação estocástica em relação à avaliação numérica em função do tamanho da amostra.

set.seed(547)
nsim=10^4
X=Z=c(rt(1,1))  # initialize the chain from the stationary
for (t in 2:nsim){
    Y=rnorm(1) # candidate normal
    rho=min(1, dt(Y,1)*dnorm(X[t-1])/(dt(X[t-1],1)*dnorm(Y)))
    X[t]=X[t-1] + (Y-X[t-1])*(runif(1)<rho)
    Y=rt(1,.5)  # candidate t
    rho=min(1, dt(Y,1)*dt(Z[t-1],.5)/(dt(Z[t-1],1)*dt(Y,.5)))
    Z[t]=Z[t-1] + (Y-Z[t-1])*(runif(1)<rho)
    }

par(mfrow=c(3,2),mar=c(2,2,1,1))
plot(5000:5800,X[5000:5800],type="l",lwd=2,xlab="",ylab="")
grid()
plot(5000:5800,Z[5000:5800],type="l",lwd=2,xlab="",ylab="")
grid()
hist(X,nclass=100,col="grey",main="",xlab="",ylab="",fre=FALSE,xlim=c(-10,10))
box()
grid()
curve(dt(x,1),col="sienna",lwd=2,add=TRUE)
hist(Z[abs(Z)<10],nclass=100,col="grey",main="",xlab="",ylab="",fre=FALSE,xlim=c(-10,10))
box()
grid()
curve(dt(x,1),col="sienna",lwd=2,add=TRUE)
acf(X,lag.max=50,lwd=2,col="gold3");grid()
acf(Z,lag.max=50,lwd=2,col="gold3",ylim=c(-0.1,0.1));grid()

Figura. 6.4. Comparação de dois esquemas Metropolis-Hastings para um alvo Cauchy ao gerar (esquerda) de uma proposta \(N(0,1)\) e (direita) de uma proposta \(T_{1/2}\) com base em \(10^5\) simulações. (topo) Trecho das cadeias \(\{X^{(t)}\}\); (centro) histogramas das amostras; (abaixo) gráficos de autocorrelação obtidos por acf.

par(mfrow=c(1,1),mar=c(4,4,1,1))
plot(cumsum(X<3)/(1:nsim),lwd=2,ty="l",ylim=c(.85,1),xlab="iterations",ylab="")
lines(cumsum(Z<3)/(1:nsim),lwd=2,col="sienna")
grid()

Figura. 6.5. Exemplo 6.2: gráfico de cobertura acumulada de uma sequência de Cauchy gerada por um algoritmo Metropolis-Hastings baseado em uma proposta \(N(0,1)\) (linhas superiores) e um gerado por um algoritmo Metropolis-Hastings baseado em uma proposta \(T_{1/2}\) (linhas inferiores). Após \(10^5\) iterações, o algoritmo Metropolis-Hastings baseado na proposta normal ainda não convergiu.


Agora veremos um exemplo estatístico um pouco mais realista que corresponde ao cenário geral quando uma proposta independente é derivada de uma estimativa preliminar dos parâmetros do modelo. Por exemplo, ao simular a partir de uma distribuição a posteriori \(\pi(\theta |x)\approx \pi(\theta)f(x|\theta)\), esta proposta independente poderia ser uma distribuição normal ou \(t\)-Student centrada no estimador de máxima verossimilhança \(\widehat{\theta}\) e com matriz de variância-covariância igual à inversa da matriz de informação de Fisher.


Exemplo 6.3

O conjunto de dados cars relaciona a distância de frenagem (\(y\)) à velocidade (\(x\)) em uma amostra de carros. A Figura 6.6 mostra os dados junto com uma curva quadrática ajustada que é dada pela função R lm. O modelo proposto para este conjunto de dados é um modelo quadrático \[ y_{ij}=a+b x_{i}+c x_i^2 +\epsilon_{ij}, \qquad i=1,\cdots,k, \quad j=1,\cdots,n_i, \] onde assumimos que \(\epsilon_{ij} \sim N(0,\sigma^2)\) e são independentes. 

library(datasets)
Nsim = 10^3
x = cars[, 1]
y = cars[, 2]
n = length(x)
x2 = x^2
MSR = 225
SSR = 10809
bhat = coefficients(lm(y ~ x + x2))
loglike = function(a, b, c, s2) {
  -(n/2) * log(s2) - sum((y - a - b * x - c * x2)^2)/(2 * s2)
  }
dcand = function(a, b, c, s2) {
  dnorm(a, bhat[1], sd = 15, log = TRUE) + dnorm(b, bhat[2], sd = 2, log = TRUE)
    +dnorm(c, bhat[3], sd = 0.06, log = TRUE) - (n/2) * log(s2) - SSR/(2 * s2)
  }
b1hat = array(bhat[1], dim = c(nsim, 1))
b2hat = array(bhat[2], dim = c(nsim, 1))
b3hat = array(bhat[3], dim = c(nsim, 1))
s2hat = array(MSR, dim = c(nsim, 1))
for (i in 2:nsim) {
  bcand = c(rnorm(1, mean = bhat[1], sd = 15), 
            rnorm(1, mean = bhat[2], sd = 2), 
            rnorm(1, mean = bhat[3], sd = 0.06), 1/rgamma(1, n/2, rate = SSR/2))
  test = min(exp(loglike(bcand[1], bcand[2], bcand[3], bcand[4]) 
            - loglike(b1hat[i - 1], b2hat[i - 1], b3hat[i - 1], 
            s2hat[i - 1]) + dcand(b1hat[i - 1], b2hat[i - 1], 
            b3hat[i - 1], s2hat[i - 1]) - dcand(bcand[1], 
            bcand[2], bcand[3], bcand[4])), 1)
  rho <- (runif(1) < test)
  b1hat[i] = bcand[1] * rho + b1hat[i - 1] * (1 - rho)
  b2hat[i] = bcand[2] * rho + b2hat[i - 1] * (1 - rho)
  b3hat[i] = bcand[3] * rho + b3hat[i - 1] * (1 - rho)
  s2hat[i] = bcand[4] * rho + s2hat[i - 1] * (1 - rho)
}
par(mar=c(4,4,1,1))
plot(x, b1hat[nsim] + b2hat[nsim] * x + b3hat[nsim] * x2, 
        type = "l", col = "grey", xlab = "", ylab = "", ylim = c(0, 
            120), lwd = 2)
for (i in (nsim - 500):nsim) {
    lines(x, b1hat[i] + b2hat[i] * x + b3hat[i] * x2, col = "grey", lwd = 2)
}
lines(x, bhat[1] + bhat[2] * x + bhat[3] * x2, col = "red", lwd = 2)
points(x, y, pch = 19)
grid()

Figura. 6.6. Dados de frenagem com curva quadrática (escuro) ajustados com a função de mínimos quadrados lm. As curvas cinza representam a amostra de Monte Carlo \((a^{(i)},b^{(i)},c^{(i)})\) e mostram a variabilidade nas linhas ajustadas com base nas últimas 500 iterações de 4000 simulações.

A função de verossimilhança é então proporcional a \[ \left(\dfrac{1}{\sigma^2}\right)^{N/2}\exp\left( -\dfrac{1}{2\sigma^2} \sum_{ij} \big( y_{ij}-a-b x_i-c x_i^2\big)^2\right), \] onde \(N =\sum_i n_i\) é o número total de observações.

Podemos ver essa função de verossimilhança como uma distribuição a a posteriori em \(a\), \(b\), \(c\) e \(\sigma^2\), podemos tentar amostrar dessa distribuição com um algoritmo Metropolis-Hastings, já que essa distribuição padrão pode ser simulada diretamente; consulte o Exercício 6.12. Para começar, podemos obter um candidato gerando coeficientes de acordo com sua distribuição amostral ajustada. 

summary(lm(y ~ x + x2))
## 
## Call:
## lm(formula = y ~ x + x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.720  -9.184  -3.188   4.628  45.152 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  2.47014   14.81716   0.167    0.868
## x            0.91329    2.03422   0.449    0.656
## x2           0.09996    0.06597   1.515    0.136
## 
## Residual standard error: 15.18 on 47 degrees of freedom
## Multiple R-squared:  0.6673, Adjusted R-squared:  0.6532 
## F-statistic: 47.14 on 2 and 47 DF,  p-value: 5.852e-12

Conforme sugerido acima, podemos usar a distribuição normal candidata centrada nos MLEs, \[ a\sim N\big(2.47, (14.82)^2\big), \quad b\sim N\big(0.91, (2.03)^2 \big), \quad c\sim N\big(0.10,(0.06)^2 \big), \] \[ \sigma^{-2} \sim Gama\big(n/2,(n-2)(15.18)^2 \big), \] em um algoritmo Metropolis-Hastings para gerar amostras \((a^{(i)},b^{(i)},c^{(i)})\). A Figura 6.6 ilustra a variabilidade das curvas associadas ao resultado desta simulação.



6.4 Seleção de candidatos

O estudo de algoritmos Metropolis-Hastings independentes é certamente interessante, mas sua implementação prática é mais problemática porque eles são delicados para usar em configurações complexas; porque a construção da proposta é complicada, se estamos usando simulação, muitas vezes é porque derivando estimativas como MLEs é difícil, e porque a escolha da proposta é altamente influente no desempenho do algoritmo.

Em vez de construir uma proposta a partir do zero ou sugerir uma aproximação não paramétrica com base em uma corrida preliminar, porque é improvável que funcione para dimensões moderadas a altas, é mais realista coletar informações sobre o alvo passo a passo, ou seja, explorando a vizinhança do valor atual da cadeia.

Se o mecanismo de exploração tiver energia suficiente para atingir os limites do suporte do alvo \(f\), o método acabará descobrindo a complexidade do alvo. Essa é basicamente a mesma intuição em ação no algoritmo de anelamento simulado da Seção 5.3.3 e no método de gradiente estocástico da Seção 5.3.2.


6.4.1 Passeios aleatórios

Uma abordagem mais natural para a construção prática de uma proposta Metropolis-Hastings é, portanto, levar em conta o valor simulado anteriormente para gerar o seguinte valor; ou seja, considerar uma exploração local da vizinhança do valor atual da Cadeia de Markov.

A implementação desta ideia é simular \(Y_t\) de acordo com \[ Y_t = X^{(t)} + \epsilon_t, \] onde \(\epsilon_t\) é uma perturbação aleatória com distribuição \(g\) independente de \(X^{(t)}\), por exemplo, uma distribuição uniforme ou uma distribuição normal, significando que \(Y_t\sim U(X^{(t)}-\delta, X^{(t)}+\delta)\) ou \(Y_t\sim N(X^{(t)},\tau^2)\) em configurações unidimensionais.

Em termos do algoritmo geral de Metropolis–Hastings, a densidade proposta \(q(y|x)\) está agora na forma \(g(y-x)\). A Cadeia de Markov associada a \(q\) é um passeio aleatório, como descrito na Seção 6.2, quando a densidade \(g\) é simétrica em torno de zero; isto é, satisfazendo \(g(-t) = g(t)\). Mas, devido ao passo de aceitação Metropolis-Hastings adicional, a Cadeia de Markov Metropolis-Hastings \(\{X^{(t)}\}\) não é um passeio aleatório.

Essa abordagem leva ao seguinte algoritmo Metropolis-Hastings, que também é o original proposto por Metrópolis e outros. (1953).


Algoritmo 6 Caminhada aleatória Metrópolis–Hastings
Dado \(x^{(t)}\),

  1. Gerar, \(Y_t\sim g(y- x^{(t)})\);

  2. Aceitar \(X^{(t+1)} = \left\{ \begin{array}{rcl} Y_t & \mbox{ com probabilidade} & \min\{1, f(Y_t)/f(x^{(t)})\} \\ x^{(t)} & \mbox{ caso contrário} & \end{array}\right.\),


Conforme observado acima, a probabilidade de aceitação não depende de \(g\). Isso significa que, para um determinado par \((x^{(t)} , y_t)\), a probabilidade de aceitação é a mesma, quer \(y_t\) seja gerado a partir de uma distribuição normal ou de Cauchy.

Obviamente, a mudança de \(g\) resultará em diferentes intervalos de valores para os \(Y_t\)’s e uma taxa de aceitação diferente, portanto, isso não quer dizer que a escolha de \(g\) não tenha impacto algum no comportamento do algoritmo, mas essa invariância da aceitação probabilidade é digna de nota. Na verdade, está ligada ao fato de que, para qualquer densidade simétrica \(g\), a medida invariante associada ao passeio aleatório é a medida de Lebesgue no espaço correspondente, ver Meyn and Tweedie (1993).


Exemplo 6.4

O exemplo histórico de Hastings (1970) considera o problema formal de geração da distribuição normal \(N(0,1)\) a partir de uma proposta de passeio aleatório igual à distribuição uniforme em \([−\delta, \delta]\).

A probabilidade de aceitação é então \[ \rho(x^{(t)}, y_t) = \exp\left(\big({x^{(t)}}^2 − y_t^2\big)/2\right) \wedge 1\cdot \]

A Figura 6.7 descreve três amostras de 5.000 pontos produzidos por esse método para \(\delta\) = 0.1, 1 e 10 e mostra claramente a diferença nas cadeias produzidas: Candidato muito estreito ou muito largo, isto é, um valor menor ou maior de \(\delta\) resultados em autocovariância mais alta e convergência mais lenta. Observe os padrões distintos para \(\delta\) = 0.1 e \(\delta\) = 10 nos gráficos superiores: No primeiro caso, a Cadeia de Markov se move a cada iteração, mas muito lentamente, enquanto no último ela permanece constante por longos períodos de tempo. 

Nsim = 10^3
a = c(0.1, 1, 10)
na = length(a)
x = array(0, c(na, Nsim))
for (i in 1:na) {
    acc = 0
    for (j in 2:Nsim) {
        y <- x[i, (j - 1)] + runif(1, min = -a[i], max = a[i])
        r = min(exp(-0.5 * ((y^2) - (x[i, (j - 1)]^2))), 1)
        u <- runif(1)
        acc = acc + (u < r)
        x[i, j] <- y * (u < r) + x[i, (j - 1)] * (u > r)
    }
    print(acc/Nsim)
}
## [1] 0.985
## [1] 0.817
## [1] 0.147
par(mfrow = c(3, na), mar = c(4, 4, 2, 1))
for (i in 1:na) plot((Nsim - 500):Nsim, x[i, (Nsim - 500):Nsim], 
    ty = "l", lwd = 2, xlab = "Iterations", ylab = "", main = paste("Rate", 
        (length(unique(x[i, ]))/Nsim), sep = " "))
for (i in 1:na) {
    hist(x[i, ], freq = F, xlim = c(-4, 4), ylim = c(0, 0.4), 
          col = "grey", ylab = "", xlab = "", breaks = 35, main = "")
    curve(dnorm(x), lwd = 2, add = T)
}
for (i in 1:na) acf(x[i, ], main = "")

Figura. 6.7. Resultados dos algoritmos Metropolis–Hastings de passeio aleatório para o Exemplo 6.4. O painel esquerdo tem um candidato \(U(-0.1, 0.1)\), o painel do meio tem \(U(-1, 1)\) e o painel direito tem \(U(-10, 10)\). Os gráficos superiores representam as últimas 500 iterações das cadeias, os gráficos intermediários indicam como os histogramas se ajustam ao alvo e os gráficos inferiores fornecem as respectivas funções de autocovariância.


Conforme observado neste exemplo formal, calibrar a escala \(\delta\) do passeio aleatório é crucial para alcançar uma boa aproximação da distribuição de destino em um número razoável de iterações. Em situações mais realistas, essa calibração se torna uma questão desafiadora, parcialmente abordada na Seção 6.5 e reconsiderada com mais detalhes no Capítulo 8.


Exemplo 6.5

O exemplo de mistura detalhado no Exemplo 5.2 da perspectiva da estimativa de uma estimativa de máxima verossimilhança também pode ser considerada do ponto de vista Bayesiano usando, por exemplo, um prior uniforme U(−2, 5) em µ1 e µ2 . A distribuição posterior na qual estamos interessados é então proporcional à verossimilhança.

A implementação do Algoritmo 6 neste exemplo é surpreendentemente fácil, pois podemos reciclar a maior parte da implementação do Algoritmo 2 de recozimento simulado, já programado no Exemplo 5.2. De fato, o núcleo do código R é muito semelhante, exceto pelo aumento da temperatura, que obviamente não é necessário aqui: 

Niter = 10^4
lange = FALSE
scale = 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]))
    iter = 1
    x = runif(2, -2, 5)
    the = matrix(x, ncol = 2)
    curlike = hval = like(x)
while (iter < Niter) {
          prop = the[iter, ] + rnorm(2) * scale
          if ((max(-prop) > 2) || (max(prop) > 5) || (log(runif(1)) > 
                like(prop) - curlike)) 
                prop = the[iter, ]
            curlike = like(prop)
            hval = c(hval, curlike)
            the = rbind(the, prop)
            iter = iter + 1
        }
        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(unique(the[, 1]), unique(the[, 2]), cex = 0.6, 
            pch = 19)

Ni=10^4
lange=TRUE
scale=.1
mhmix=function(Ni=10^4, lange=TRUE, scale=c(0.1,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]))
    iter = 1
    x = runif(2, -2, 5)
    the = matrix(x, ncol = 2)
    curlike = hval = like(x)
    if (!lange) {
      while (iter < Niter) {
        prop = the[iter, ] + rnorm(2) * scale
        if ((max(-prop) > 2) || (max(prop) > 5) || (log(runif(1)) > 
                like(prop) - curlike)) 
                prop = the[iter, ]
        curlike = like(prop)
        hval = c(hval, curlike)
        the = rbind(the, prop)
        iter = iter + 1
      }
        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 = T)
        points(unique(the[, 1]), unique(the[, 2]), cex = 0.6, 
            pch = 19)
    }
    else {
        gradlike = function(mu) {
            gr = sum(0.2 * (da - mu[1]) * dnorm(da - mu[1])/((0.2 * 
                dnorm(da - mu[1]) + 0.8 * dnorm(da - mu[2]))))
            gr = c(gr, sum(0.8 * (da - mu[2]) * dnorm(da - mu[2])/((0.2 * 
                dnorm(da - mu[1]) + 0.8 * dnorm(da - mu[2])))))
            return(gr)
        }
        curmean = x + 0.5 * scale^2 * gradlike(x)
        while (iter < Niter) {
            prop = curmean + rnorm(2) * scale
            meanprop = prop + 0.5 * scale^2 * gradlike(prop)
            if ((max(-prop) > 2) || (max(prop) > 5) || (log(runif(1)) > 
                like(prop) - curlike - sum(dnorm(prop, curmean, 
                  lo = T)) + sum(dnorm(the[iter, ], meanprop, 
                  lo = T)))) {
                prop = the[iter, ]
                meanprop = curmean
            }
            curlike = like(prop)
            curmean = meanprop
            hval = c(hval, curlike)
            the = rbind(the, prop)
            iter = iter + 1
        }
        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 = T)
        points(unique(the[, 1]), unique(the[, 2]), cex = 0.6, 
            pch = 19)
        lines(unique(the[, 1]), unique(the[, 2]), cex = 0.6, 
            pch = 19)
    }
}


Um problema que geralmente surge ao usar passeios aleatórios em domínios restritos é se o passeio aleatório deve ou não ser restrito também. A resposta a essa pergunta é não, pois o uso de restrições na proposta modifica a função \(g\) e, portanto, compromete a validade da razão dos alvos encontrada no Algoritmo 6.

Quando são propostos valores \(y_t\) fora do intervalo de \(f\), ou seja, quando \(f(y_t) = 0\), o valor proposto é rejeitado e o valor atual \(X^{(t)}\) é duplicado. Obviamente, escolher uma densidade de passeio aleatório que muitas vezes termina fora do domínio de \(f\) é uma má ideia, pois a cadeia ficará presa na maior parte do tempo! Mas é formalmente correto.


6.4.2 Candidatos alternativos

Enquanto o algoritmo Metropolis-Hastings independente se aplica apenas em situações específicas, o algoritmo Metropolis-Hastings de passeio aleatório geralmente aparece como um algoritmo genérico de Metropolis-Hastings que atende à maioria dos casos.

No entanto, a solução de passeio aleatório não é necessariamente a escolha mais eficiente porque

  1. requer muitas iterações para superar dificuldades como regiões de baixa probabilidade entre regiões modais de \(f\) e

  2. por causa de suas características simétricas, ele gasta aproximadamente metade do tempo de simulação revisitando regiões que já explorou.

Existem alternativas que contornam a simetria perfeita na proposta de passeio aleatório para ganhar em eficiência, embora nem sempre sejam fáceis de implementar ver, por exemplo, Robert and Casella (2004).

Uma dessas alternativas é o algoritmo de Langevin de Roberts and Rosenthal (1998), que tenta favorecer movimentos em direção a valores mais altos do alvo \(f\) incluindo um gradiente na proposta, \[ Y_t=X^{(t)}+\dfrac{\sigma^2}{2}\nabla \log\big( f(X^{(t)})\big)+\sigma \epsilon_t, \qquad \epsilon_t \sim g(\epsilon), \] sendo o parâmetro \(\sigma\) o fator de escala da proposta.

Quando \(Y_t\) é construído dessa maneira, a probabilidade de aceitação Metropolis-Hastings é igual a \[ \rho(x,y)=\min\left\{ \dfrac{f(y)}{f(x)}\dfrac{g\big( (x-y)/\sigma -\sigma \nabla \log(f(y)/2)\big)}{g\big( (y-x)/\sigma -\sigma \nabla \log(f(x)/2)\big)}, 1\right\}\cdot \]

Embora esse esquema possa lembrá-lo das técnicas de gradiente estocástico da Seção 5.3.2, ele difere delas por dois motivos. Uma delas é que a escala é fixa no algoritmo de Langevin, em vez de diminuir no método de gradiente estocástico. Outra é que o movimento proposto para \(Y_t\) não é necessariamente aceito pelo algoritmo de Langevin, garantindo a estacionariedade de \(f\) para a cadeia resultante.


Exemplo 6.6

Com base no mesmo modelo probit do agora conhecido conjunto de dados Pima.tr como no Exemplo 3.10, podemos usar a função de verossimilhança como foi definida: 

set.seed(675)
Niter = 10^4
scale = 0.01
library(MASS)
da = cbind(Pima.tr$type, Pima.tr$bmi)
da[, 1] = da[, 1] - 1
like = function(a, b) {
    sum(pnorm(q = a + outer(X = b, Y = da[, 2], FUN = "*"), 
        log = T) * da[, 1] + pnorm(q = -a - outer(X = b, 
        Y = da[, 2], FUN = "*"), log = T) * (1 - da[, 1]))
}

e calcular o gradiente na forma fechada.

grad = function(a, b) {
    don = pnorm(q = a + outer(X = b, Y = da[, 2], FUN = "*"))
    x1 = sum((dnorm(x = a + outer(X = b, Y = da[, 2], FUN = "*"))/don) * 
        da[, 1] - (dnorm(x = -a - outer(X = b, Y = da[, 2], 
        FUN = "*"))/(1 - don)) * (1 - da[, 1]))
    x2 = sum((dnorm(x = a + outer(X = b, Y = da[, 2], FUN = "*"))/don) * 
        da[, 2] * da[, 1] - (dnorm(x = -a - outer(X = b, 
        Y = da[, 2], FUN = "*"))/(1 - don)) * da[, 2] * (1 - 
        da[, 1]))
    return(c(x1, x2))
}

Ao implementar a iteração básica do algoritmo de Langevin

the = matrix(glm(da[, 1] ~ da[, 2], family = binomial(link = "probit"))$coef, ncol = 2)
curmean = the[1, ] + 0.5 * scale^2 * grad(the[1, 1], the[1, 2])
likecur = like(the[1, 1], the[1, 2])
for (t in 2:Niter) {
    prop = curmean + scale * rnorm(2)
    propmean = prop + 0.5 * scale^2 * grad(prop[1], prop[2])
    if (log(runif(1)) > like(prop[1], prop[2]) - likecur - 
        sum(dnorm(prop, mean = curmean, sd = scale, lo = T)) + 
        sum(dnorm(the[t - 1, ], mean = propmean, sd = scale, lo = T))) {
            prop = the[t - 1, ]
            propmean = curmean
        }
    the = rbind(the, prop)
    curmean = propmean
}

precisamos selecionar uma escala pequena o suficiente porque, caso contrário, grad(prop) retorna NaN dado que


pnorm(q=a+outer(X=b,Y=da[,2],FUN="*"))
é 1 ou 0.

Com escala igual a 0.01, a cadeia explora corretamente a distribuição a posteriori, conforme mostra a Figura 6.9, embora se mova muito lentamente. 

be1 = seq(min(the[, 1]), max(the[, 1]), le = 100)
be2 = seq(min(the[, 2]), max(the[, 2]), le = 130)
li = matrix(0, ncol = 130, nro = 100)
for (i in 1:100) for (j in 1:130) li[i, j] = like(be1[i], be2[j])
par(mar = c(4, 4, 1, 1))
image(be1, be2, li, xlab = expression(beta[1]), ylab = expression(beta[2]))
contour(be1, be2, li, add = T, ncla = 100)
subs = seq(1, Niter, le = 10^3)
points(unique(the[subs, 1]), unique(the[subs, 2]), cex = 0.4, pch = 19)

Figura. 6.9. Repartição da amostra de Langevin correspondente ao probit a posteriori definido no Exemplo 3.10 com base em 20 observações de Pima.tr e \(5\times 10^4\) iterações.


A modificação da proposta de passeio aleatório pode, no entanto, dificultar ainda mais a mobilidade da cadeia de Markov por reforçar a polarização em torno dos modos locais. Por exemplo, quando o alvo é a distribuição a posteriori do modelo de mistura estudado no Exemplo 6.5, a estrutura bimodal do alvo é um empecilho para a implementação do algoritmo de Langevin em que o modo local se torna ainda mais atrativo.


Exemplo 6.7 (Continuação do Exemplo 6.5)

A modificação do algoritmo de passeio aleatório Metropolis-Hastings é direta, pois simplesmente temos que adicionar o desvio do gradiente no código R.

Definindo a função gradiente 

mhmix = function (Niter = 10^4, lange = FALSE, scale = 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]))
    iter = 1
    x = runif(2, -2, 5)
    the = matrix(x, ncol = 2)
    curlike = hval = like(x)
    if (!lange) {
        while (iter < Niter) {
            prop = the[iter, ] + rnorm(2) * scale
            if ((max(-prop) > 2) || (max(prop) > 5) || (log(runif(1)) > 
                like(prop) - curlike)) 
                prop = the[iter, ]
            curlike = like(prop)
            hval = c(hval, curlike)
            the = rbind(the, prop)
            iter = iter + 1
        }
        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 = T)
        points(unique(the[, 1]), unique(the[, 2]), cex = 0.6, 
            pch = 19)
    }
    else {
        gradlike = function(mu) {
            gr = sum(0.2 * (da - mu[1]) * dnorm(da - mu[1])/((0.2 * 
                dnorm(da - mu[1]) + 0.8 * dnorm(da - mu[2]))))
            gr = c(gr, sum(0.8 * (da - mu[2]) * dnorm(da - mu[2])/((0.2 * 
                dnorm(da - mu[1]) + 0.8 * dnorm(da - mu[2])))))
            return(gr)
        }
        curmean = x + 0.5 * scale^2 * gradlike(x)
        while (iter < Niter) {
            prop = curmean + rnorm(2) * scale
            meanprop = prop + 0.5 * scale^2 * gradlike(prop)
            if ((max(-prop) > 2) || (max(prop) > 5) || (log(runif(1)) > 
                like(prop) - curlike - sum(dnorm(prop, curmean, 
                  lo = T)) + sum(dnorm(the[iter, ], meanprop, 
                  lo = T)))) {
                prop = the[iter, ]
                meanprop = curmean
            }
            curlike = like(prop)
            curmean = meanprop
            hval = c(hval, curlike)
            the = rbind(the, prop)
            iter = iter + 1
        }
        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 = T)
        points(unique(the[, 1]), unique(the[, 2]), cex = 0.6, 
            pch = 19)
        lines(unique(the[, 1]), unique(the[, 2]), cex = 0.6, 
            pch = 19)
    }
}
mhmix(scale=0.1)

mhmix()

Figura. 6.10. Exploração das modas no modelo de mistura por um algoritmo de Langevin: representação de duas cadeias de Markov \((\mu_1^{(t)},\mu_2^{(t)})\) no topo da superfície log-posterior com escala igual a 0.1 com base em \(10^4\) simulações e um conjunto de dados simulado de 500 observações.


Ambos os exemplos acima mostram como o ajuste do algoritmo de Langevin pode ser delicado. Isso pode explicar por que não é amplamente implementado, embora seja uma modificação bastante fácil do código básico de passeio aleatório.

Algoritmos de passeio aleatório Metropolis-Hastings também se aplicam a alvos de suporte discretos. Embora isso soe mais como uma configuração combinatória ou de processamento de imagem, uma vez que a maioria dos problemas estatísticos envolve espaços de parâmetros contínuos, uma exceção é o caso da escolha do modelo ver, por exemplo, Robert (2001), onde o índice do modelo a ser selecionado é o parâmetro de interesse.


Exemplo 6.8

Dada uma regressão linear comum com \(n\) observações, \[ y \, | \, \beta,\sigma^2,X \sim N_n(X\beta,\sigma^2\pmb{I}_n), \] onde \(X\) é uma matriz \(n\times p\), a verossimilhança é \[ \ell(\beta,\sigma^2 \, | \, y,X) = \left(\dfrac{1}{2\pi\sigma^2}\right)^{n/2}\exp\left( -\dfrac{1}{2\sigma^2} \big( y-X\beta\big)^\top \big( y-X\beta\big)\right) \] e, sob o chamado \(g\)-a priori de Zellner (1986), \[ \beta \, | \, \sigma^2,X \sim N_{k+1}\big(\widetilde{\beta},n\sigma^2(X^\top X)^{-1} \big) \qquad \mbox{e} \qquad \pi(\sigma^2 \, | \, X) \approx \sigma^2, \] onde a constante \(g\) é escolhida igual a \(n\), a distribuição marginal de \(y\) é uma distribuição \(t\)-Student multivariada, \[ \begin{array}{l} m(y\, | \, X)=\dfrac{1}{(n+1)^{(k+1)/2}}\dfrac{1}{\pi^{n/2}}\Gamma(n/2)\times \\ \qquad \qquad \qquad \qquad \qquad \qquad \left(y^\top y -\dfrac{n}{n+1}y^\top X(X^\top X)^{-1}X^\top y-\dfrac{1}{n+1}\widetilde{\beta}^\top X^\top X^\top X\widetilde{\beta} \right)^{-n/2}\cdot \end{array} \]

Como ilustração, consideramos o conjunto de dados swiss, onde o logaritmo da fertilidade em 47 distritos da Suíça por volta de 1888 é a variável \(y\) a ser explicada por alguns indicadores socioeconômicos,

library(datasets)
data("swiss")
y=log(as.vector(swiss[,1]))
X=as.matrix(swiss[,2:6])

A matriz de covariáveis \(X\) envolve cinco variáveis explicativas

names(swiss)
## [1] "Fertility"        "Agriculture"      "Examination"      "Education"       
## [5] "Catholic"         "Infant.Mortality"
que são explicados por ?swiss e queremos comparar os \(2^5\) modelos correspondentes a todos os possíveis subconjuntos de covariáveis.

Neste exemplo, o número de modelos é pequeno o suficiente para permitir o cálculo de todos os modelos marginais e, portanto, as verdadeiras probabilidades de todos os modelos em comparação. Seguindo Marin and Robert (2007), indexamos todos os modelos por vetores \(\gamma\) de indicadores binários onde \(\gamma_i = 0\) indica que a coluna correspondente de \(X\) é usada na regressão.

Observe que, adotando a convenção de Marin and Robert (2007), sempre incluímos o intercepto num modelo. Usando a função de matriz inversa rápida

inv=function(X){
  EV=eigen(X)
  EV$vector%*%diag(1/EV$values)%*%t(EV$vector)
}

então calculamos o logaritmo da densidade marginal correspondente ao modelo \(\gamma\), agora denotado como $m(y , | ,X,), como

t1 = function(gam) {
  p = length(gam)
  if (sum(gam) != 0) 
    order(gam)[(p - sum(gam) + 1):p]
  else 0
}
lpostw = function(gam, y, X, betatilde) {
  n = length(y)
  qgam = sum(gam)
  t1gam = t1(gam)
  Xt1 = cbind(rep(1, n), X[, t1gam])
  if (qgam != 0) 
    P1 = Xt1 %*% inv(t(Xt1) %*% Xt1) %*% t(Xt1)
  else P1 = matrix(0, n, n)
  -(qgam + 1)/2 * log(n + 1) - n/2 * log(t(y) %*% y - n/(n + 
  1) * t(y) %*% P1 %*% y - 1/(n + 1) * t(betatilde) %*% 
  t(cbind(rep(1, n), X)) %*% P1 %*% cbind(rep(1, n), X) %*% betatilde)
}

A exploração do espaço dos modelos pode resultar de um algoritmo Metropolis-Hastings que se move em torno dos modelos alterando um indicador de modelo por vez; ou seja, dado o vetor indicador atual \(\gamma^{(t)}\), a proposta Metropolis-Hastings escolhe uma das \(p\) coordenadas, digamos \(i\), e escolhe entre manter \(\gamma_i^{(t)}\) e mudar para \(1-\gamma_i^{(t)}\) com probabilidades proporcionais às marginais associadas.

A probabilidade de aceitação Metropolis-Hastings do modelo proposto \(\gamma^*\) é então igual a \[ \min\left\{ \dfrac{m(y\, | \, X,\gamma^*)}{m(y\, | \, X,\gamma^{(t)})}\dfrac{m(y\, | \, X,\gamma^{(t)})}{m(y\, | \, X,\gamma^*)},1\right\}=1 \] uma vez que as constantes de normalização se cancelam.

Isso significa que não temos que considerar a rejeição do modelo proposto \(\gamma^*\) porque é sempre aceito na etapa Metropolis-Hastings.

Executando a função R

set.seed(2976)
mochoice=function (Niter = 10^4) 
{
    library(datasets)
    y = log(as.vector(swiss[, 1]))
    X = as.matrix(swiss[, 2:6])
    inv = function(X) {
        EV = eigen(X)
        EV$vector %*% diag(1/EV$values) %*% t(EV$vector)
    }
    t1 = function(gam) {
        p = length(gam)
        if (sum(gam) != 0) 
            order(gam)[(p - sum(gam) + 1):p]
        else 0
    }
    lpostw = function(gam, y, X, betatilde) {
        n = length(y)
        qgam = sum(gam)
        t1gam = t1(gam)
        Xt1 = cbind(rep(1, n), X[, t1gam])
        if (qgam != 0) 
            P1 = Xt1 %*% inv(t(Xt1) %*% Xt1) %*% t(Xt1)
        else P1 = matrix(0, n, n)
        -(qgam + 1)/2 * log(n + 1) - n/2 * log(t(y) %*% y - n/(n + 
            1) * t(y) %*% P1 %*% y - 1/(n + 1) * t(betatilde) %*% 
            t(cbind(rep(1, n), X)) %*% P1 %*% cbind(rep(1, n), X) %*% betatilde)
    }
    integerate = function(gamma) {
        lga = length(gamma)
        vale = gamma[1]
        for (i in 2:lga) vale = vale + gamma[i] * 2^(i - 1)
        vale
    }
    back2bin = function(k) {
        gam = k%/%2^4
        res = k%%2^4
        for (i in 3:0) {
            logam = res%/%2^i
            res = res%%2^i
            gam = c(logam, gam)
        }
        gam
    }
    gocho = function(niter, y, X) {
        lga = dim(X)[2]
        betatilde = lm(y ~ X)$coeff
        gamma = matrix(0, niter, lga)
        gamma[1, ] = sample(c(0, 1), lga, rep = T)
        for (i in 1:(niter - 1)) {
            j = sample(1:lga, 1)
            gam0 = gam1 = gamma[i, ]
            gam1[j] = 1 - gam0[j]
            pr = lpostw(gam0, y, X, betatilde)
            pr = c(pr, lpostw(gam1, y, X, betatilde))
            pr = exp(pr - max(pr))
            gamma[i + 1, ] = gam0
            if (sample(c(0, 1), 1, prob = pr)) 
                gamma[i + 1, ] = gam1
        }
        gamma
    }
    out = gocho(Niter, y, X)
    inde = apply(out, 1, integerate)
    cole = rep(0, length(unique(inde)))
    for (i in 1:length(unique(inde))) cole[i] = (sum(inde == unique(inde)[i]))
    besm = matrix(0, ncol = 5, nrow = 5)
    for (i in 1:5) besm[i, ] = back2bin(unique(inde)[order(-jitter(cole))[i]])
    list(model = out, top = besm, means = apply(out, 2, mean))
}
então produz uma amostra (aproximadamente) distribuída a partir da distribuição a posteriori no conjunto de indicadores; ou seja, na coleção de possíveis submodelos.

Com base no resultado

resultado=mochoice(10^5)
resultado$model[1:10,]
##       [,1] [,2] [,3] [,4] [,5]
##  [1,]    0    1    1    0    1
##  [2,]    0    1    1    0    1
##  [3,]    0    1    1    0    1
##  [4,]    0    1    1    0    1
##  [5,]    0    1    1    0    1
##  [6,]    0    1    1    1    1
##  [7,]    0    1    1    1    1
##  [8,]    0    1    1    1    1
##  [9,]    0    1    1    1    1
## [10,]    0    1    1    1    1
resultado$top
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0    1    1    1
## [2,]    0    0    1    1    1
## [3,]    1    1    1    1    1
## [4,]    1    0    1    1    0
## [5,]    0    1    1    1    1
o modelo mais provável corresponde à exclusão da variável Agriculture, ou seja, \(\gamma= (1,0,1,1,1)\), com probabilidade estimada 0.4995, enquanto a probabilidade verdadeira é 0.4997. Este modelo também é o indicado por lm(y ~ X).

Da mesma forma, o segundo modelo mais provável é \(\gamma= (0,0,1,1,1)\), com uma probabilidade estimada de 0.237 versus uma probabilidade verdadeira de 0.234. A probabilidade de cada variável ser incluída no modelo também é fornecida por

resultado$means
## [1] 0.66234 0.19184 0.99999 0.91099 0.94865
o que, novamente, indica que as três últimas variáveis de swiss são as mais significativas nesta análise.


O fato da probabilidade de aceitação ser sempre igual a 1 no Exemplo 6.8 se deve ao uso da verdadeira probabilidade alvo em um subconjunto dos valores possíveis do indicador do modelo.


6.5 Taxas de aceitação

Existem infinitas escolhas para a distribuição candidata \(q\) em um algoritmo Metropolis-Hastings, e aqui discutimos a possibilidade de alcançar uma escolha ótima. Obviamente, esse não é um conceito bem definido, pois a escolha ótima de \(q\) é tomar \(q = f\), a distribuição de destino, ao raciocinar em termos de velocidade de convergência.

Este é obviamente um resultado formal que não tem relevância na prática. Em vez disso, precisamos adotar um critério prático que permita a comparação dos núcleos propostos em situações em que quase nada se sabe sobre \(f\). Um desses critérios é a taxa de aceitação do algoritmo Metropolis-Hastings correspondente, pois pode ser facilmente calculado ao executar esse algoritmo por meio da frequência empírica de aceitação.

Em contraste com o Capítulo 2, onde a calibração de um algoritmo Accept-Reject foi baseada na taxa de aceitação máxima, simplesmente otimizar a taxa de aceitação não resultará necessariamente no melhor algoritmo em termos de mistura e convergência.


Exemplo 6.9

Em um algoritmo Accept-Reject gerando uma amostra \(N(0,1)\) de uma distribuição exponencial dupla \(\mathcal{L}(\alpha)\) com densidade \(g(x|\alpha) = (\alpha/2) \exp(−\alpha |x|)\), a escolha \(\alpha = 1\) otimiza a taxa de aceitação (Exercício 2.19). Podemos usar essa distribuição como um candidato independente \(q\) em um algoritmo Metropolis-Hastings.

A Figura 6.11 compara o comportamento desse candidato \(\mathcal{L}(1)\) junto com uma distribuição \(\mathcal{L}(3)\), que, para esta simulação, produz um resultado inferior no sentido de ter autocovariâncias maiores e, consequentemente, convergência mais lenta.

Obviamente, uma análise mais profunda seria necessária para validar essa afirmação, mas nosso ponto aqui é que a taxa de aceitação estimada para \(\alpha = 1\) é quase duas vezes maior, 0.84, do que a taxa de aceitação estimada para \(\alpha = 3\), 0.499. 

set.seed(9462)
Nsim = 10^3
a = 3
dl = function(x, a) (a/2) * exp(-a * abs(x))
X1 = X2 = array(0, dim = c(Nsim, 1))
X1[1] = X2[1] = rnorm(1)
acc1 = acc2 = 0
for (i in 2:Nsim) {
  Y1 = ifelse(runif(1) > 0.5, 1, -1) * rexp(1)
  Y2 = Y1/a
  rho1 = min(1, dnorm(Y1) * dl(X1[i - 1], 1)/(dnorm(X1[i - 1]) * dl(Y1, 1)))
  rho2 = min(1, dnorm(Y2) * dl(X2[i - 1], a)/(dnorm(X2[i - 1]) * dl(Y2, a)))
  test1 = (runif(1) < rho1)
  acc1 = acc1 + test1
  test2 = (runif(1) < rho2)
  acc2 = acc2 + test2
  X1[i] = X1[i - 1] + (Y1 - X1[i - 1]) * test1
  X2[i] = X2[i - 1] + (Y2 - X2[i - 1]) * test2
}
print(c(acc1/Nsim, acc2/Nsim))
## [1] 0.840 0.524
print(c(mean(X1[(Nsim/2):Nsim]), var(X1[(Nsim/2):Nsim])))
## [1] -0.0226247  1.0047279
print(c(mean(X2[(Nsim/2):Nsim]), var(X2[(Nsim/2):Nsim])))
## [1] -0.1009395  0.9831432
par(mfrow = c(1, 3), mar = c(4, 4, 2, 1))
plot(cumsum(X1)/(1:Nsim), type = "l", ylim = c(-0.25, 0.25), 
        lwd = 2, xlab = "Iterations", ylab = "", col = "blue")
lines(cumsum(X2)/(1:Nsim), lwd = 2, xlab = "", ylab = "")
grid()
acf(X1)
grid()
acf(X2)
grid()
Figura. 6.11. Gráfico de média acumulada (esquerda) de um algoritmo de Metropolis-Hastings usado para gerar uma variável aleatória \(N(0,1)\) de uma distribuição proposta exponencial dupla \(\mathcal{L}(1)\) (mais clara) e \(\mathcal{L}(3)\) (preto). Os painéis central e esquerdo mostram a autocovariância para as propostas \(\mathcal{L}(1)\) e \(\mathcal{L}(3)\), respectivamente.


Embora os algoritmos Metropolis-Hastings independentes possam de fato ser otimizados ou pelo menos comparados por meio de sua taxa de aceitação, porque isso reduz o número de réplicas na cadeia \(\{X^{(t)}\}\) e, portanto, o nível de correlação na cadeia, isso não é verdade para outros tipos de algoritmos de Metropolis-Hastings, em primeiro lugar a versão de passeio aleatório.

A versão de caminhada aleatória do algoritmo Metropolis-Hastings, apresentada na Seção 6.4.1, de fato requer uma abordagem diferente para as taxas de aceitação, dada a dependência da distribuição do candidato no estado atual da cadeia. De fato, como já visto no Exemplo 6.4, uma alta taxa de aceitação não indica necessariamente que o algoritmo está se comportando satisfatoriamente, pois pode corresponder ao fato de que a cadeia está se movendo muito lentamente na superfície de \(f\).

Quando \(x^{(t)}\) e \(y_t\) são próximos, no sentido de que \(f(x^{(t)})\) e \(f(y_t)\) são aproximadamente iguais, o algoritmo Metropolis–Hastings de passeio aleatório leva à aceitação de \(y_t\) com probabilidade \[ \min\left\{ f(y_t)/f(x^{(t)}),1\right\} \approx 1\cdot \]

Uma alta taxa de aceitação pode, portanto, sinalizar um padrão de convergência ruim, pois os movimentos no suporte de \(f\) são mais limitados. Obviamente, nem sempre é esse o caso. Por exemplo, quando \(f\) é quase plano, altas taxas de aceitação não são indicativas de nenhum comportamento errado. Mas, a menos que \(f\) seja completamente plano, ou seja, corresponda a um alvo uniforme, há partes do domínio a serem exploradas onde \(f\) assume valores menores e, portanto, onde as probabilidades de aceitação devem ser pequenas. Uma alta taxa de aceitação indica que essas partes do domínio não são frequentemente ou nem um pouco, exploradas pelo algoritmo Metropolis-Hastings.

Em contraste, se a taxa média de aceitação for baixa, os valores sucessivos de \(f(y_t)\) muitas vezes são pequenos quando comparados com \(f(x^{(t)})\), o que corresponde ao cenário em que o passeio aleatório se move rapidamente na superfície de \(f\) desde frequentemente atinge as bordas do suporte de \(f\) ou pelo menos quando o passeio aleatório explora regiões com baixa probabilidade sob \(f\).

Novamente, uma baixa taxa de aceitação não significa que a cadeia explore todo o suporte de \(f\). Mesmo com uma pequena taxa de aceitação, pode faltar um modo importante, mas isolado, de \(f\). No entanto, uma baixa taxa de aceitação é um problema menor, exceto do ponto de vista do tempo de computação, porque indica explicitamente que um número maior de simulações é necessário. Usar o tamanho efetivo da amostra como um indicador de convergência, ver Seção 8.4.3, sinalizaria claramente esse requisito.


Exemplo 6.10. (Continuação do Exemplo 6.4)

Os três algoritmos de caminhada aleatória Metropolis–Hastings da Figura 6.7 têm taxas de aceitação iguais a 0.977, 0.8 e 0.171, respectivamente.

Observando o ajuste do histograma, vemos que a taxa de aceitação média se sai melhor, mas que a taxa de aceitação mais baixa ainda se sai melhor do que a mais alta.


A questão é, então, decidir sobre uma taxa de aceitação de ouro contra a qual calibrar os algoritmos de passeio aleatório Metropolis-Hastings, a fim de evitar taxas de aceitação “muito altas” e “muito baixas”. Roberts et al. (1997) recomendam o uso de distribuições instrumentais com taxas de aceitação próximas a 1/4 para modelos de alta dimensão e iguais a 1/2 para os modelos de dimensão 1 ou 2.

Esta é a regra adotada no pacote amcmc descrito na Seção 8.5.2. Embora esta regra não seja universal, no sentido de que foi projetada principalmente para um ambiente Gaussiano, nós a defendemos como uma meta de calibração padrão sempre que pode ser alcançado, o que nem sempre é o caso. Por exemplo, se considerarmos o algoritmo Metropolis–Hastings do Exemplo 6.8, não há taxa de aceitação, pois a probabilidade de aceitação é sempre igual a 1. No entanto, como a proposta inclui o valor atual em seu suporte, a cadeia \(\{\gamma^{(t)}\}\) tem valores idênticos em uma linha e, portanto, uma taxa de aceitação ou renovação implícita.

Esta taxa é igual a 0.1805, muito abaixo da meta de 0.25, e o algoritmo não pode ser facilmente modificado, por exemplo, observando mais movimentos alternativos em torno do modelo atual, para atingir essa taxa de aceitação de token.


6.6 Exercícios