Algoritmo de Metropolis–Hastings

Em estatística e física estatística, o algoritmo Metropolis-Hastings é um método de Cadeia de Markov Monte Carlo (MCMC) para obter amostras aleatórias a partir de uma distribuição de probabilidade da qual a amostragem direta é difícil. Essa sequência pode ser usada para aproximar a distribuição (por exemplo, para gerar um histograma) ou para calcular uma integral (por exemplo, um valor esperado). Metropolis – Hastings e outros algoritmos MCMC são geralmente usados para amostragem de distribuições multidimensionais, especialmente quando o número de dimensões é alto. Para distribuições unidimensionais, existem outros métodos (por exemplo, amostragem por rejeição adaptativa) que podem retornar amostras independentes, fugindo do problema de amostras autocorrelacionadas inerente aos métodos MCMC.

História editar

O algoritmo recebeu o nome de Nicholas Metropolis, autor do artigo de 1953 sobre Equação de Cálculos Estatais por Máquinas de Computação Rápida, juntamente com Arianna W. Rosenbluth, Marshall Rosenbluth, Marshall Rosenbluth, Augusta H. Teller e Edward Teller. Este artigo propôs o algoritmo para o caso de distribuições simétricas de propostas, e WK Hastings o estendeu ao caso mais geral em 1970.

Existe alguma controvérsia no que diz respeito ao crédito pelo desenvolvimento do algoritmo. Metropolis cunhou o termo "Monte Carlo" em um artigo anterior com Stanislav Ulam, estava familiarizado com os aspectos computacionais do método e liderou o grupo na Divisão Teórica que projetou e construiu o computador MANIAC I usado nos experimentos em 1952. No entanto, antes de 2003, não havia uma descrição detalhada do desenvolvimento do algoritmo. Então, pouco antes de sua morte, Marshall Rosenbluth participou de uma conferência de 2003 na LANL, comemorando o 50º aniversário da publicação de 1953. Nesta conferência, Rosenbluth descreveu o algoritmo e seu desenvolvimento em uma apresentação intitulada "Gênese do algoritmo de Monte Carlo para mecânica estatística". Esclarecimentos históricos adicionais são feitos por Gubernatis em um artigo de jornal de 2005 recontando a conferência do 50º aniversário. Rosenbluth deixa claro que ele e sua esposa Arianna fizeram o trabalho e que Metropolis não teve nenhum papel no desenvolvimento além de fornecer tempo ao computador.

Intuição editar

O algoritmo Metropolis–Hastings pode coletar amostras de qualquer distribuição de probabilidade $P(x)$ , desde que conheçamos uma função $f(x)$ proporcional à densidade de $P$ e os valores de $f(x)$ possam ser calculados. A exigência de que $f(x)$ deve ser apenas proporcional à densidade, e não exatamente igual a ela, torna o algoritmo Metropolis-Hastings particularmente útil, porque calcular na prática o fator de normalização necessário costuma ser extremamente difícil.

O algoritmo Metropolis-Hastings funciona gerando uma sequência de valores de tal maneira que, à medida que mais e mais valores de amostra são produzidos, a distribuição de valores se aproxima mais da distribuição desejada $P(x)$ . Esses valores de amostra são produzidos iterativamente, com a distribuição da próxima amostra dependendo apenas do valor atual da amostra (transformando a sequência de amostras em uma cadeia de Markov). Especificamente, a cada iteração, o algoritmo escolhe um candidato para o próximo valor de amostra com base no valor atual da amostra. Então, com alguma probabilidade, o candidato é aceito ou rejeitado. Caso seja rejeitado, o valor do candidato é descartado e o valor atual é reutilizado na próxima iteração. A probabilidade de aceitação é determinado comparando os valores da função $f(x)$ dos valores de amostra atuais e candidatos em relação à distribuição desejada $P(x)$ .

Para fins de ilustração, o algoritmo Metropolis, um caso especial do algoritmo Metropolis-Hastings, em que a função da proposta é simétrica, é descrito abaixo.

Algoritmo de Metropolis (distribuição simétrica de propostas)

Deixe $f(x)$ ser uma função proporcional à distribuição de probabilidade desejada $P(x)$ (também conhecida como distribuição de destino).

Inicialização: escolha um ponto arbitrário $x_{t}$ para ser a primeira amostra e escolher uma densidade de probabilidade arbitrária $g(x\mid y)$ (às vezes escrito $Q(x\mid y)$ ) que sugere um candidato para o próximo valor da amostra $x$ , dado o valor da amostra anterior $y$ . Nesta secção, $g$ é assumido como simétrico; em outras palavras, ele deve satisfazer $g(x\mid y)=g(y\mid x)$ . Uma escolha usual é deixar $g(x\mid y)$ ser uma distribuição gaussiana centrada em $y$ , para que aponte para mais perto $y$ são mais propensos a serem visitados a seguir, transformando a sequência de amostras em uma caminhada aleatória . A função $g$ é chamado de densidade da proposta ou distribuição de salto .
Para cada iteração t :
- Gere um candidato $x'$ para a próxima amostra, escolhendo na distribuição $g(x'\mid x_{t})$ .
- Calcule a taxa de aceitação $\alpha =f(x')/f(x_{t})$ , que será usado para decidir se aceita ou rejeita o candidato. Como f é proporcional à densidade de P, temos que $\alpha =f(x')/f(x_{t})=P(x')/P(x_{t})$ .
- Aceite ou rejeite :
  - Gere um número aleatório uniforme $u\in [0,1]$ .
  - E se $u\leq \alpha$ , aceite o candidato definindo $x_{t+1}=x'$ ,
  - E se $u>\alpha$ , rejeite o candidato e defina $x_{t+1}=x_{t}$ em vez de.

Esse algoritmo prossegue ao tentar mover-se aleatoriamente no espaço da amostra, às vezes aceitando as movimentações e às vezes permanecendo no local. Observe que a taxa de aceitação $\alpha$ indica a probabilidade da nova amostra proposta em relação à amostra atual, de acordo com a distribuição $P(x)$ . Se tentarmos avançar para um ponto mais provável que o ponto existente (isto é, um ponto em uma região de maior densidade de $P(x)$ ), sempre aceitaremos a mudança. No entanto, se tentarmos avançar para um ponto menos provável, às vezes rejeitaremos a mudança e, quanto mais a queda relativa na probabilidade, maior a probabilidade de rejeitarmos o novo ponto. Assim, tenderemos a permanecer (e retornar um grande número de amostras de) regiões de alta densidade de $P(x)$ , enquanto visitam apenas ocasionalmente regiões de baixa densidade. Intuitivamente, é por isso que esse algoritmo funciona e retorna amostras que seguem a distribuição desejada $P(x)$ .

Comparado com um algoritmo como a amostragem por rejeição adaptativa^[1] que gera diretamente amostras independentes de uma distribuição, o Metropolis – Hastings e outros algoritmos do MCMC têm algumas desvantagens:

As amostras estão correlacionadas. Embora, a longo prazo, sigam corretamente $P(x)$ , um conjunto de amostras próximas será correlacionado entre si e não refletirá corretamente a distribuição. Isso significa que, se queremos um conjunto de amostras independentes, temos que jogar fora a maioria das amostras e coletar apenas a n- amostra, por algum valor de n (geralmente determinado examinando a autocorrelação entre amostras adjacentes). A autocorrelação pode ser reduzida aumentando a largura do salto (o tamanho médio de um salto, que está relacionado à variância da distribuição do salto), mas isso também aumentará a probabilidade de rejeição do salto proposto. Um tamanho de salto muito grande ou muito pequeno levará a uma cadeia de Markov de mistura lenta, ou seja, um conjunto de amostras altamente correlacionado, de modo que um número muito grande de amostras será necessário para obter uma estimativa razoável de qualquer propriedade desejada da distribuição.
Embora a cadeia de Markov eventualmente converja para a distribuição desejada, as amostras iniciais podem seguir uma distribuição muito diferente, especialmente se o ponto de partida estiver em uma região de baixa densidade. Como resultado, normalmente é necessário um período de queima,^[2] onde um número inicial de amostras (por exemplo, as primeiras 1000) são jogadas fora.

Por outro lado, os métodos mais simples de amostragem de rejeição sofrem com a " maldição da dimensionalidade ", onde a probabilidade de rejeição aumenta exponencialmente em função do número de dimensões. Metropolis – Hastings, juntamente com outros métodos MCMC, não têm esse problema a esse grau e, portanto, são frequentemente as únicas soluções disponíveis quando o número de dimensões da distribuição a ser amostrada é alto. Como resultado, os métodos MCMC são frequentemente escolhidos para produzir amostras a partir de modelos hierárquicos bayesianos e demais modelos estatísticos de alta dimensionalidade.

Nas distribuições multivariadas, o algoritmo clássico Metropolis-Hastings, conforme descrito acima, envolve a escolha de um novo ponto de amostra multidimensional. Quando o número de dimensões é alto, é difícil encontrar a distribuição de salto adequada a ser usada, pois as diferentes dimensões individuais se comportam de maneiras muito diferentes, e a largura do salto (veja acima) deve ser "perfeita" para todas as dimensões de uma só vez para evitar misturas excessivamente lentas. Uma abordagem alternativa que geralmente funciona melhor nessas situações, conhecida como amostragem de Gibbs, envolve a escolha de uma nova amostra para cada dimensão separadamente das outras, em vez de escolher uma amostra para todas as dimensões de uma só vez. Isso é especialmente aplicável quando a distribuição multivariada é composta por um conjunto de variáveis aleatórias individuais , nas quais cada variável é condicionada apenas a um pequeno número de outras variáveis, como é o caso dos modelos hierárquicos mais comuns. As variáveis individuais são então amostradas uma de cada vez, com cada variável condicionada aos valores mais recentes de todas as outras. Vários algoritmos podem ser usados para escolher essas amostras individuais, dependendo da forma exata da distribuição multivariada: algumas possibilidades são os métodos de amostragem por rejeição adaptativa,^[1]^[3]^[4]^[5] o algoritmo de amostragem por metrópole de rejeição adaptativa^[6] ou suas melhorias^[7]^[8] (consulte o código do matlab), uma etapa unidimensional simples de Metropolis – Hastings ou amostragem de fatias.

Derivação formal editar

O objetivo do algoritmo Metropolis – Hastings é gerar uma coleção de estados conforme uma distribuição desejada $P(x)$ . Para conseguir isso, o algoritmo usa um processo de Markov, que atinge assintoticamente uma distribuição estacionária exclusiva $\pi (x)$ de tal modo que $\pi (x)=P(x)$ .

Um processo de Markov é definido exclusivamente por suas probabilidades de transição $P(x'\mid x)$ , a probabilidade de fazer a transição de qualquer estado $x$ para qualquer outro estado $x'$ . O processo possui uma distribuição estacionária única $\pi (x)$ quando as duas condições a seguir são atendidas:

Existência de distribuição estacionária : deve existir uma distribuição estacionária $\pi (x)$ . Uma condição suficiente (mas não necessária) é o equilíbrio detalhado, que exige que cada transição $x\to x'$ seja reversível: para cada par de estados $x,x'$ , a probabilidade de estar no estado $x$ e transitar para o estado $x'$ deve ser igual à probabilidade de estar no estado $x'$ e transitar para o estado $x$ , $\pi (x)P(x'\mid x)=\pi (x')P(x\mid x')$ .
Exclusividade da distribuição estacionária : a distribuição estacionária $\pi (x)$ deve ser única. Isso é garantido pela ergodicidade do processo de Markov, que exige que todo estado (1) seja aperiódico (o sistema não retorna ao mesmo estado em intervalos fixos); e (2) seja positivo recorrente (o número esperado de etapas para retornar ao mesmo estado é finito).

O algoritmo Metropolis-Hastings envolve o desenho de uma cadeia de Markov (construindo probabilidades de transição) que atenda às duas condições acima, de modo que sua distribuição estacionária $\pi (x)$ seja $P(x)$ . A derivação do algoritmo começa com a condição de equilíbrio detalhado :

P(x'\mid x)P(x)=P(x\mid x')P(x'),

que é reescrita como

{\frac {P(x'\mid x)}{P(x\mid x')}}={\frac {P(x')}{P(x)}}.

A abordagem é separar a transição em duas sub-etapas; a proposta e a aceitação-rejeição. A distribuição da proposta $g(x'\mid x)$ é a probabilidade condicional de caminhar para um estado $x'$ dado $x$ e a distribuição de aceitação $A(x',x)$ é a probabilidade de aceitar o estado proposto $x'$ e, de fato, permanecer neste. A probabilidade de transição pode ser escrita como o produto deles:

P(x'\mid x)=g(x'\mid x)A(x',x).

Inserindo essa relação na equação anterior, temos

{\frac {A(x',x)}{A(x,x')}}={\frac {P(x')}{P(x)}}{\frac {g(x\mid x')}{g(x'\mid x)}}.

O próximo passo na derivação é escolher uma taxa de aceitação que atenda à condição acima. Uma escolha comum é a escolha de Metropolis:

A(x',x)=\min \left(1,{\frac {P(x')}{P(x)}}{\frac {g(x\mid x')}{g(x'\mid x)}}\right).

Para essa taxa de aceitação do Metropolis $A$ , ou $A(x',x)=1$ ou $A(x,x')=1$ e, de qualquer forma, a condição é satisfeita.

O algoritmo Metropolis – Hastings consiste, portanto, no seguinte:

Inicialize
1. Escolha um estado inicial $x_{0}$ .
2. Coloque $t=0$ .
Itere
1. Gere um estado candidato aleatório $x'$ de acordo com $g(x'\mid x_{t})$ .
2. Calcule a probabilidade de aceitação $A(x',x_{t})=\min \left(1,{\frac {P(x')}{P(x_{t})}}{\frac {g(x_{t}\mid x')}{g(x'\mid x_{t})}}\right)$ .
3. Aceite ou rejeite :
  1. gere um número aleatório uniforme $u\in [0,1]$ ;
  2. E se $u\leq A(x',x_{t})$ , aceite o novo estado e defina $x_{t+1}=x'$ ;
  3. E se $u>A(x',x_{t})$ , rejeite o novo estado e copie o estado antigo para a frente $x_{t+1}=x_{t}$ .
4. Incremente : coloque $t=t+1$ .

Desde que as condições especificadas sejam atendidas, a distribuição empírica dos estados salvos $x_{0},\ldots ,x_{T}$ vai se aproximar de $P(x)$ . O número de iterações ( $T$ ) necessário para estimar efetivamente $P(x)$ depende de vários fatores, como a relação entre $P(x)$ e a distribuição da proposta e a precisão desejada.^[9] Para distribuição em espaços de estado discretos, deve ser da ordem do tempo de autocorrelação do processo de Markov.

É importante notar que, em um problema geral, não é óbvio qual distribuição $g(x'\mid x)$ deve-se usar ou o número de iterações necessárias para uma estimativa adequada; ambos são parâmetros livres do método, que devem ser ajustados ao problema específico em questão.

Uso na integração numérica editar

Um uso comum do algoritmo Metropolis-Hastings é calcular uma integral. Especificamente, considere um espaço $\Omega \subset \mathbb {R}$ e uma distribuição de probabilidade $P(x)$ sobre $\Omega$ , $x\in \Omega$ . Metropolis-Hastings pode estimar uma integral da forma de

P(E)=\int _{\Omega }A(x)P(x)\,dx,

Onde $A(x)$ é uma função (mensurável) de interesse.

Por exemplo, considere uma estatística $E(x)$ e sua distribuição de probabilidade $P(E)$ , que é uma distribuição marginal . Suponha que o objetivo seja estimar $P(E)$ para $E$ na cauda de $P(E)$ . Formalmente, $P(E)$ pode ser escrito como

P(E)=\int _{\Omega }P(E\mid x)P(x)\,dx=\int _{\Omega }\delta {\big (}E-E(x){\big )}P(x)\,dx=E{\big (}P(E\mid X){\big )}

e, assim, pode-se estimar $P(E)$ estimando primeiro o valor esperado da função indicadora $A_{E}(x)\equiv \mathbf {1} _{E}(x)$ , que é 1 quando $E(x)\in [E,E+\Delta E]$ e 0 caso contrário. Como $E$ está na cauda de $P(E)$ , a probabilidade de tomar um estado $x$ com $E(x)$ na cauda de $P(E)$ é proporcional a $P(E)$ , que é pequeno por definição. O algoritmo Metropolis-Hastings pode ser usado aqui para amostrar estados (raros) com maior probabilidade e, assim, aumentar o número de amostras usadas para estimar $P(E)$ nas caudas. Isso pode ser feito, por exemplo, usando uma distribuição de amostragem $\pi (x)$ favorecer esses estados (por exemplo, $\pi (x)\propto e^{aE}$ com $a>0$ )

Instruções passo a passo editar

Suponha que o valor mais recente amostrado seja $x_{t}$ . Para seguir o algoritmo Metropolis-Hastings, pegamos um novo estado proposto $x'$ com densidade de probabilidade $g(x'\mid x_{t})$ e calculamos um valor

a=a_{1}a_{2},

onde

a_{1}={\frac {P(x')}{P(x_{t})}}

é a razão de probabilidade (por exemplo, a distribuição posterior bayesiana) entre a amostra proposta $x'$ e a amostra anterior $x_{t}$ e

a_{2}={\frac {g(x_{t}\mid x')}{g(x'\mid x_{t})}}

é a razão da densidade da proposta em duas direções (de $x_{t}$ para $x'$ e vice-versa). Isso é igual a 1 se a densidade da proposta for simétrica. Então o novo estado $x_{t+1}$ é escolhido de acordo com as seguintes regras.

Se

a\geq 1{:}

x_{t+1}=x',

De outro modo:

x_{t+1}={\begin{cases}x'&{\text{com probabilidade }}a,\\x_{t}&{\text{com probabilidade }}1-a.\end{cases}}

A cadeia de Markov é iniciada a partir de um valor inicial arbitrário $x_{0}$ , e o algoritmo é executado para muitas iterações até que esse estado inicial seja "esquecido". Essas amostras, que são descartadas, são conhecidas como burn-in . O conjunto restante de valores aceitos de $x$ representam uma amostra de valores da distribuição $P(x)$ .

O algoritmo funciona melhor se a densidade da proposta corresponder à forma da distribuição de destino $P(x)$ , da qual a amostragem direta é difícil (ou seja, $g(x'\mid x_{t})\approx P(x')$ ) . Se uma densidade de proposta gaussiana $g$ é usado, o parâmetro de variância $\sigma ^{2}$ precisa ser ajustado durante o período de queima. Isso geralmente é feito calculando a taxa de aceitação, que é a fração de amostras propostas que é aceita em uma janela da última $N$ amostras. A taxa de aceitação desejada depende da distribuição de destino. No entanto, foi demonstrado teoricamente que a taxa de aceitação ideal para uma distribuição Gaussiana unidimensional é de cerca de 50%, diminuindo para cerca de 23% para uma distribuição Gaussiana $N$ -dimensional.

Se $\sigma ^{2}$ for muito pequeno, a cadeia se misturará lentamente (ou seja, a taxa de aceitação será alta, mas amostras sucessivas se moverão lentamente pelo espaço, e a cadeia convergirá apenas lentamente para $P(x)$ ) Por outro lado, se $\sigma ^{2}$ for muito grande, a taxa de aceitação será muito baixa porque é provável que as propostas cheguem a regiões com uma densidade de probabilidade muito menor. $a_{1}$ será muito pequena e, novamente, a cadeia convergirá muito lentamente. Geralmente, a distribuição da proposta é ajustada para que os algoritmos aceitem na ordem de 30% de todas as amostras - de acordo com as estimativas teóricas mencionadas no parágrafo anterior.

O resultado de três cadeias de Markov em execução na função Rosenbrock 3D usando o algoritmo Metropolis–Hastings. O algoritmo coleta amostras de regiões onde a probabilidade posterior é alta e as cadeias começam a se misturar nessas regiões. A posição aproximada do máximo foi iluminada. Observe que os pontos vermelhos são os que permanecem após o processo de gravação. Os anteriores foram descartados.

Ver também editar

Referências

↑ ^a ^b Gilks. «Adaptive Rejection Sampling for Gibbs Sampling». Journal of the Royal Statistical Society. Series C (Applied Statistics). 41: 337–348. JSTOR 2347565. doi:10.2307/2347565
↑ Bayesian data analysis. Chapman & Hall / CRC 2nd ed. Boca Raton, Fla.: [s.n.] 2004. ISBN 978-1584883883. OCLC 51991499
↑ Görür. «Concave-Convex Adaptive Rejection Sampling». Journal of Computational and Graphical Statistics. 20: 670–691. ISSN 1061-8600. doi:10.1198/jcgs.2011.09058
↑ Hörmann. «A Rejection Technique for Sampling from T-concave Distributions». ACM Trans. Math. Softw. 21: 182–193. CiteSeerX 10.1.1.56.6055 . ISSN 0098-3500. doi:10.1145/203082.203089
↑ Martino. «A generalization of the adaptive rejection sampling algorithm». Statistics and Computing. 21: 633–647. ISSN 0960-3174. doi:10.1007/s11222-010-9197-9
↑ Gilks. «Adaptive Rejection Metropolis Sampling within Gibbs Sampling». Journal of the Royal Statistical Society. Series C (Applied Statistics). 44: 455–472. JSTOR 2986138. doi:10.2307/2986138
↑ Martino. «Independent Doubly Adaptive Rejection Metropolis Sampling Within Gibbs Sampling». IEEE Transactions on Signal Processing. 63: 3123–3138. Bibcode:2015ITSP...63.3123M. ISSN 1053-587X. arXiv:1205.5494 . doi:10.1109/TSP.2015.2420537
↑ Meyer. «Adaptive rejection Metropolis sampling using Lagrange interpolation polynomials of degree 2». Computational Statistics & Data Analysis. 52: 3408–3423. doi:10.1016/j.csda.2008.01.005
↑ Raftery, Adrian E., and Steven Lewis. "How Many Iterations in the Gibbs Sampler?" In Bayesian Statistics 4. 1992.

[:0-1] Gilks. «Adaptive Rejection Sampling for Gibbs Sampling». Journal of the Royal Statistical Society. Series C (Applied Statistics). 41: 337–348. JSTOR 2347565. doi:10.2307/2347565

[2] Bayesian data analysis. Chapman & Hall / CRC 2nd ed. Boca Raton, Fla.: [s.n.] 2004. ISBN 978-1584883883. OCLC 51991499

[3] Görür. «Concave-Convex Adaptive Rejection Sampling». Journal of Computational and Graphical Statistics. 20: 670–691. ISSN 1061-8600. doi:10.1198/jcgs.2011.09058

[4] Hörmann. «A Rejection Technique for Sampling from T-concave Distributions». ACM Trans. Math. Softw. 21: 182–193. CiteSeerX 10.1.1.56.6055 . ISSN 0098-3500. doi:10.1145/203082.203089

[5] Martino. «A generalization of the adaptive rejection sampling algorithm». Statistics and Computing. 21: 633–647. ISSN 0960-3174. doi:10.1007/s11222-010-9197-9

[6] Gilks. «Adaptive Rejection Metropolis Sampling within Gibbs Sampling». Journal of the Royal Statistical Society. Series C (Applied Statistics). 44: 455–472. JSTOR 2986138. doi:10.2307/2986138

[7] Martino. «Independent Doubly Adaptive Rejection Metropolis Sampling Within Gibbs Sampling». IEEE Transactions on Signal Processing. 63: 3123–3138. Bibcode:2015ITSP...63.3123M. ISSN 1053-587X. arXiv:1205.5494 . doi:10.1109/TSP.2015.2420537

[8] Meyer. «Adaptive rejection Metropolis sampling using Lagrange interpolation polynomials of degree 2». Computational Statistics & Data Analysis. 52: 3408–3423. doi:10.1016/j.csda.2008.01.005

[9] Raftery, Adrian E., and Steven Lewis. "How Many Iterations in the Gibbs Sampler?" In Bayesian Statistics 4. 1992.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]