Comprensión de MCMC y el algoritmo Metropolis-Hastings

En los últimos días he estado tratando de entender cómo funciona Markov Chain Monte Carlo (MCMC). En particular, he estado tratando de entender e implementar el algoritmo Metropolis-Hastings. Hasta ahora creo que tengo una comprensión general del algoritmo, pero hay algunas cosas que aún no me quedan claras. Quiero usar MCMC para adaptar algunos modelos a los datos. Debido a esto, describiré mi comprensión del algoritmo Metropolis-Hastings para ajustar una línea recta a algunos datos observados :$f(x)=ax$$D$

1) Haga una suposición inicial para . Establezca esto como nuestro actual ( ). También agregue al final de la Cadena de Markov ( ).a a a 0 a $a$$a$$a$$a_0$$a$$C$

2) Repita los pasos a continuación varias veces.

3) evaluar la probabilidad actual ( ) dada y . a 0${\cal L_0}$$a_0$$D$

4) Proponga una nueva ( ) muestreando desde una distribución normal con y . Por ahora, es constante.a 1 μ = a 0 σ = s t e p s i z e s t e p s i z $a$$a_1$$\mu=a_0$$\sigma=stepsize$$stepsize$

5) Evaluar nueva verosimilitud ( ) dada y . a 1${\cal L_1}$$a_1$$D$

6) Si es mayor que , acepte como el nuevo , agréguelo al final de y vaya al paso 2.L 0 a 1 a 0${\cal L_1}$${\cal L_0}$$a_1$$a_0$$C$

7) Si es menor que genere un número ( ) en el rango [0,1] a partir de una distribución uniformeL 0${\cal L_1}$${\cal L_0}$$U$

8) Si es menor que la diferencia entre las dos probabilidades ( - ), acepte como el nuevo , agréguelo al final de y vaya al paso 2.L 1 L 0 a 1 a 0$U$${\cal L_1}$${\cal L_0}$$a_1$$a_0$$C$

9) Si es mayor que la diferencia entre las dos probabilidades ( - ), agregue al final de , siga usando el mismo , vaya al paso 2.L 1 L 0 a 0 C a $U$${\cal L_1}$${\cal L_0}$$a_0$$C$$a_0$

10) Fin de la repetición.

11) Elimine algunos elementos del inicio de (fase de quemado).$C$

12) Ahora tomar el promedio de los valores de . Este promedio es el estimado .$C$$a$

Ahora tengo algunas preguntas sobre los pasos anteriores:

• ¿Cómo construyo la función de probabilidad para pero también para cualquier función arbitraria?$f(x)=ax$
• ¿Es esta una implementación correcta del algoritmo Metropolis-Hastings?
• ¿Cómo la selección del método de generación de números aleatorios en el Paso 7 puede cambiar los resultados?
• ¿Cómo va a cambiar este algoritmo si tengo varios parámetros de modelo? Por ejemplo, si tuviera el modelo .$f(x)=ax+b$

Notas / Créditos: La estructura principal del algoritmo descrito anteriormente se basa en el código de un Taller de Python de MPIA.

Parece haber algunas ideas erróneas sobre el algoritmo Metropolis-Hastings (MH) en su descripción del algoritmo.

En primer lugar, uno tiene que entender que MH es un algoritmo de muestreo. Como se indica en wikipedia

En estadística y en física estadística, el algoritmo Metropolis-Hastings es un método de Monte Carlo en cadena de Markov (MCMC) para obtener una secuencia de muestras aleatorias a partir de una distribución de probabilidad para la cual el muestreo directo es difícil.

$Q(\cdot\vert\cdot)$$f(\cdot)$

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$

Un ejemplo en R se puede encontrar en el siguiente enlace:

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

Este método se emplea en gran medida en las estadísticas bayesianas para el muestreo de la distribución posterior de los parámetros del modelo.

$f(x)=ax$$x$

Robert y Casella (2010), Introducing Monte Carlo Methods with R , Ch. 6, "Algoritmos de Metrópolis-Hastings"

También hay muchas preguntas, con punteros a referencias interesantes, en este sitio que discuten sobre el significado de la función de probabilidad.

Otro indicador de posible interés es el paquete R mcmc, que implementa el algoritmo MH con propuestas gaussianas en el comando metrop().