Entendendo o MCMC e o algoritmo Metropolis-Hastings

Nos últimos dias, tenho tentado entender como o Markov Chain Monte Carlo (MCMC) funciona. Em particular, tenho tentado entender e implementar o algoritmo Metropolis-Hastings. Até agora, acho que tenho uma compreensão geral do algoritmo, mas há algumas coisas que ainda não estão claras para mim. Eu quero usar o MCMC para ajustar alguns modelos aos dados. Por causa disso, descreverei minha compreensão do algoritmo Metropolis-Hastings para ajustar uma linha reta $f(x)=ax$ a alguns dados observados $D$ :

1) Faça um palpite inicial para . Defina esta um como o nosso atual um ( a 0 ). Adicione também um no final da cadeia de Markov ( C ).$a$$a$$a$$a_0$$a$$C$

2) Repita as etapas abaixo várias vezes.

3) avaliar a probabilidade actual ( ) dado um 0 e D .${\cal L_0}$$a_0$$D$

4) Propor um novo ( a 1 ) amostrando a partir de uma distribuição normal com μ = a 0 e σ = s t e p s i z e . Por agora, s t e p s i z e é constante.$a$$a_1$$\mu=a_0$$\sigma=stepsize$$stepsize$

5) Avaliar nova probabilidade ( ) deu um 1 e D .${\cal L_1}$$a_1$$D$

6) Se for maior que L ${\cal L_1}$${\cal L_0}$ , aceitar como o novo um 0 , acrescentá-la na extremidade de C e ir para o passo 2.$a_1$$a_0$$C$

7) se for menor que L 0, gere um número (U) no intervalo [0,1] a partir de uma distribuição uniforme${\cal L_1}$${\cal L_0}$$U$

8) Se for menor que a diferença entre as duas probabilidades ( L 1 - L $U$${\cal L_1}$${\cal L_0}$ ), aceitar como o novo um 0 , acrescentá-la na extremidade de C e ir para o passo 2.$a_1$$a_0$$C$

9) Se for maior que a diferença entre as duas probabilidades ( L 1 - L $U$${\cal L_1}$${\cal L_0}$ ), acrescente a no final de C , continue usando a mesma a 0 , vá para a etapa 2.$a_0$$C$$a_0$

10) Fim da repetição.

11) Remova alguns elementos do início de (fase de queima).$C$

12) Agora pegue a média dos valores em . Essa média é a estimada a .$C$$a$

Agora eu tenho algumas perguntas sobre as etapas acima:

• Como eu construo a função de probabilidade para mas também para qualquer função arbitrária?$f(x)=ax$
• Esta é uma implementação correta do algoritmo Metropolis-Hastings?
• Como a seleção do método de geração de número aleatório na Etapa 7 pode alterar os resultados?
• Como esse algoritmo vai mudar se eu tiver vários parâmetros de modelo? Por exemplo, se eu tivesse o modelo .$f(x)=ax+b$

Notas / Créditos: A estrutura principal do algoritmo descrito acima é baseada no código de um Workshop MPIA Python.

Parece haver alguns conceitos errados sobre o que é o algoritmo Metropolis-Hastings (MH) na sua descrição do algoritmo.

Primeiro de tudo, é preciso entender que o MH é um algoritmo de amostragem. Como afirmado em wikipedia

Em estatística e em física estatística, o algoritmo Metropolis-Hastings é um método de Monte Carlo da cadeia de Markov (MCMC) para obter uma sequência de amostras aleatórias a partir de uma distribuição de probabilidade para a qual a amostragem direta é difícil.

$Q(\cdot\vert\cdot)$$f(\cdot)$ , o algoritmo MH pode ser implementado da seguinte maneira:

1. $x_0$ .
2. $x^{\star}$$Q(\cdot\vert x_0)$ .
3. $\alpha=f(x^{\star})/f(x_0)$ .
4. $x^{\star}$$f$$\alpha$
5. $x^{\star}$

$N$$k$

Um exemplo em R pode ser encontrado no seguinte link:

http://www.mas.ncl.ac.uk/~ndjw1/teaching/sim/metrop/metrop.html

Este método é amplamente empregado nas estatísticas bayesianas para amostragem a partir da distribuição posterior dos parâmetros do modelo.

$f(x)=ax$$x$

Robert & Casella (2010), Introduzindo Métodos de Monte Carlo com R , cap. 6, "Algoritmos de Metropolis-Hastings"

Também há muitas perguntas, com dicas para referências interessantes, neste site, discutindo sobre o significado da função de probabilidade.

Outro ponteiro de possível interesse é o pacote R mcmc, que implementa o algoritmo MH com propostas gaussianas no comando metrop().