Algoritmo de Metropolis-Hastings

En estadística y física estadística, el algoritmo Metropolis-Hastings es un método de Monte Carlo en cadena de Markov para obtener una secuencia de muestras aleatorias a partir de una distribución de probabilidad a partir de la cual es difícil el muestreo directo. Esta secuencia se puede usar para aproximar la distribución (por ejemplo, para generar un histograma) o para calcular una integral (por ejemplo, un valor esperado). Metropolis-Hastings y otros algoritmos Montecarlo en cadena de Markov se usan generalmente para el muestreo de distribuciones multidimensionales, especialmente cuando el número de dimensiones es alto. Para las distribuciones unidimensionales, generalmente hay otros métodos (por ejemplo, muestreo de rechazo adaptativo) que pueden devolver directamente muestras independientes de la distribución, y estos están libres del problema de las muestras autocorrelacionadas que es inherente a los métodos Montecarlo en cadena de Markov.

Historia

El algoritmo lleva el nombre de Nicholas Metropolis, autor del artículo 1953 Equation of State Calculations by Fast Computing Machines junto con Arianna W. Rosenbluth, Marshall Rosenbluth, Augusta H. Teller y Edward Teller. Este artículo propuso el algoritmo para el caso de distribuciones de propuestas simétricas, y W. K. Hastings lo extendió al caso más general en 1970.

Existe cierta controversia con respecto al crédito para el desarrollo del algoritmo. Metropolis había acuñado el término "Montecarlo" en un artículo anterior con Stanislav Ulam, estaba familiarizado con los aspectos computacionales del método y dirigió el grupo en la División Teórica que diseñó y construyó la computadora MANIAC I utilizada en los experimentos en 1952. Sin embargo, antes de 2003, no había una cuenta detallada del desarrollo del algoritmo. Luego, poco antes de su muerte, Marshall Rosenbluth asistió a una conferencia de 2003 en LANL para conmemorar el 50 aniversario de la publicación de 1953. En esta conferencia, Rosenbluth describió el algoritmo y su desarrollo en una presentación titulada "Génesis del algoritmo de Monte Carlo para la mecánica estadística". Gubernatis hace más aclaraciones históricas en un artículo de revista de 2005 que relata la conferencia del 50 aniversario. Rosenbluth deja en claro que él y su esposa Arianna hicieron el trabajo, y que Metropolis no desempeñó ningún papel en el desarrollo que no sea proporcionar tiempo de computadora.

Esto contradice una versión de Edward Teller, quien afirma en sus memorias que los cinco autores del artículo de 1953 trabajaron juntos durante "días (y noches)". Por el contrario, el relato detallado de Rosenbluth acredita a Teller una sugerencia crucial pero temprana de "aprovechar la mecánica estadística y tomar promedios de conjunto en lugar de seguir una cinemática detallada". Esto, dice Rosenbluth, lo hizo pensar en el enfoque generalizado de Monte Carlo, un tema que, según él, había discutido a menudo con Von Neumann. Arianna contó (a Gubernatis en 2003) que Augusta Teller comenzó el trabajo de la computadora, pero que la propia Arianna se hizo cargo y escribió el código desde cero. En una historia oral registrada poco antes de su muerte, Rosenbluth nuevamente acredita a Teller por plantear el problema original, a él mismo por resolverlo, y a Arianna por programar la computadora. En términos de reputación, hay pocas razones para cuestionar la cuenta de Rosenbluth. En una memoria biográfica de Rosenbluth, Freeman Dyson escribe:

Muchas veces fui con Rosenbluth, le hice una pregunta [...] y recibí una respuesta en dos minutos. Por lo general, me llevaría una semana de arduo trabajo comprender en detalle por qué la respuesta de Rosenbluth era correcta. Tenía una capacidad asombrosa para ver a través de una situación física complicada y llegar a la respuesta correcta mediante argumentos físicos. Enrico Fermi fue el único otro físico que he conocido que era igual a Rosenbluth en su comprensión intuitiva de la física.

Intuición

El algoritmo de Metropolis-Hastings puede extraer muestras de cualquier distribución de probabilidad ${\displaystyle P(x)}$ , siempre que conozcamos una función ${\displaystyle f(x)}$  proporcional a la densidad de ${\displaystyle P}$  y los valores de ${\displaystyle f(x)}$  pueden ser calculados. El único requisito es que ${\displaystyle f(x)}$  debe de ser proporcional a la densidad, en lugar de ser exactamente igual a él, hace que el algoritmo de Metropolis-Hastings sea particularmente útil, porque calcular el factor de normalización necesario a menudo es extremadamente difícil en la práctica.

El algoritmo Metropolis-Hastings funciona generando una secuencia de valores de muestra de tal manera que, a medida que se producen más y más valores de muestra, la distribución de valores se aproxima más a la distribución deseada ${\displaystyle P(x)}$ . Estos valores de muestra se producen de forma iterativa, y la distribución de la siguiente muestra depende solo del valor de muestra actual (convirtiendo así la secuencia de muestras en una cadena de Markov). Específicamente, en cada iteración, el algoritmo elige un candidato para el siguiente valor de muestra basado en el valor de muestra actual. Luego, con cierta probabilidad, el candidato es aceptado (en cuyo caso el valor candidato se usa en la próxima iteración) o rechazado (en cuyo caso el valor candidato se descarta y el valor actual se reutiliza en la próxima iteración): la probabilidad de aceptación se determina comparando los valores de la función ${\displaystyle f(x)}$  de los valores de muestra actuales y candidatos con respecto a la distribución deseada ${\displaystyle P(x)}$ .

Para fines de ilustración, a continuación se describe el algoritmo Metropolis, un caso especial del algoritmo Metropolis-Hastings donde la función de la propuesta es simétrica.

Algoritmo de Metropolis (distribución de propuesta simétrica)

Deja a ${\displaystyle f(x)}$  ser una función que sea proporcional a la distribución de probabilidad deseada ${\displaystyle P(x)}$  (un objetivo de distribución).

1. Inicialización: elija un punto arbitrario ${\displaystyle x_{t}}$  ser la primera muestra y elegir una densidad de probabilidad arbitraria ${\displaystyle g(x\mid y)}$  (a veces escrita ${\displaystyle Q(x\mid y)}$ ) que sugiere un candidato para el siguiente valor de muestra ${\displaystyle x}$ , dado el valor de muestra anterior ${\displaystyle y}$ . En esta sección, se asume que el valor de ${\displaystyle g}$  es simétrico; en otras palabras, este debe de satisfacer ${\displaystyle g(x\mid y)=g(y\mid x)}$ . Una opción habitual es dejar que ${\displaystyle g(x\mid y)}$  sea una distribución gaussiana centrada en ${\displaystyle y}$ , así los puntos cercanos a ${\displaystyle y}$  son más propensos a ser visitados a continuación, convirtiendo la secuencia de muestras en una caminata aleatoria. La función ${\displaystyle g}$  se denomina densidad de propuesta o distribución de saltos.
2. Para cada iteración t:
   • Genera un candidato ${\displaystyle x'}$  para la siguiente muestra seleccionandolo de la distribución ${\displaystyle g(x'\mid x_{t})}$ .
   • Calculate the acceptance ratio ${\displaystyle \alpha =f(x')/f(x_{t})}$ , que se utilizará para decidir si acepta o rechaza al candidato. Ya que f es proporcional a la densidad de P, podemos deducir que ${\displaystyle \alpha =f(x')/f(x_{t})=P(x')/P(x_{t})}$ .
   • Aceptar o rechazar:
     • Genera un número aleatorio uniforme ${\displaystyle u\in [0,1]}$ .
     • Si ${\displaystyle u\leq \alpha }$ , entonces acepta el candidato de forma que ${\displaystyle x_{t+1}=x'}$ ,
     • Si ${\displaystyle u>\alpha }$ , entonces rechaza el candidato y coloca ${\displaystyle x_{t+1}=x_{t}}$  en su lugar.

Este algoritmo continúa intentando moverse aleatoriamente por el espacio muestral, a veces aceptando los movimientos y otras permaneciendo en su lugar. Nótese que la relación de aceptación ${\displaystyle \alpha }$  indica qué tan probable es la nueva muestra propuesta con respecto a la muestra actual, de acuerdo con la distribución ${\displaystyle P(x)}$ . Si intentamos movernos a un punto que sea más probable que el punto existente (es decir, un punto en una región de mayor densidad de ${\displaystyle P(x)}$ ), siempre aceptaremos el movimiento. Sin embargo, si intentamos movernos a un punto menos probable, a veces rechazaremos el movimiento, y cuanto mayor sea la caída relativa de la probabilidad, más probable es que rechacemos el nuevo punto. Por lo tanto, tendremos que permanecer en (y devolver grandes cantidades de muestras de) regiones de alta densidad de ${\displaystyle P(x)}$ , mientras visita ocasionalmente regiones de baja densidad. Intuitivamente, esta es la razón por la cual este algoritmo funciona y devuelve muestras que siguen la distribución deseada ${\displaystyle P(x)}$ .

En comparación con un algoritmo como el muestreo de rechazo adaptativo^[1]​ que genera directamente muestras independientes de una distribución, Metropolis-Hastings y otros algoritmos Montecarlo en cadena de Markov tienen una serie de desventajas:

• Las muestras están correlacionadas. Aunque a largo plazo siguen correctamente ${\displaystyle P(x)}$ , un conjunto de muestras cercanas se correlacionaron entre sí y no reflejarán correctamente la distribución. Esto significa que si queremos un conjunto de muestras independientes, tenemos que tirar la mayoría de las muestras y solo se seleccionan cada cierto número de muestra, para algún valor de n (típicamente determinado al examinar la autocorrelación entre muestras adyacentes). La autocorrelación puede reducirse aumentando el ancho del salto (el tamaño promedio de un salto, que está relacionado con la variación de la distribución del salto), pero esto también aumentará la probabilidad de rechazo del salto propuesto. Un tamaño de salto demasiado grande o demasiado pequeño conducirá a una cadena de Markov de mezcla lenta, es decir, un conjunto altamente correlacionado de muestras, por lo que se necesitará una cantidad muy grande de muestras para obtener una estimación razonable de cualquier propiedad deseada de la distribución.
• Aunque la cadena de Markov finalmente converge a la distribución deseada, las muestras iniciales pueden seguir una distribución muy diferente, especialmente si el punto de partida está en una región de baja densidad. Como resultado, generalmente es necesario un período de quemado, en el que se desecha un número inicial de muestras (por ejemplo, los primeros 1000).

Por otro lado, los métodos de muestreo de rechazo más simples sufren la "maldición de la dimensionalidad", donde la probabilidad de rechazo aumenta exponencialmente en función del número de dimensiones. Metropolis-Hastings, junto con otros métodos de Montecarlo en cadena de Markov, no tienen este problema en tal grado y, por lo tanto, a menudo son las únicas soluciones disponibles cuando el número de dimensiones de la distribución a muestrear es alto. Como resultado, los métodos de este tipo son a menudo los métodos de elección para producir muestras a partir de modelos bayesianos jerárquicos y otros modelos estadísticos de alta dimensión utilizados hoy en día en muchas disciplinas.

En distribuciones multivariadas, el algoritmo clásico de Metropolis-Hastings como se describió anteriormente implica elegir un nuevo punto de muestra multidimensional. Cuando el número de dimensiones es alto, encontrar la distribución de salto adecuada para usar puede ser difícil, ya que las diferentes dimensiones individuales se comportan de maneras muy diferentes, y el ancho de salto (ver arriba) debe ser "justo" para todas las dimensiones a la vez. evite mezclar excesivamente lento. Un enfoque alternativo que a menudo funciona mejor en tales situaciones, conocido como muestreo de Gibbs, implica elegir una nueva muestra para cada dimensión por separado de las demás, en lugar de elegir una muestra para todas las dimensiones a la vez. Esto es especialmente aplicable cuando la distribución multivariada se compone de un conjunto de variables aleatorias individuales en las que cada variable está condicionada a un pequeño número de otras variables, como es el caso en la mayoría de los modelos jerárquicos típicos. Las variables individuales se muestran de una en una, con cada variable condicionada a los valores más recientes de todas las demás. Se pueden usar varios algoritmos para elegir estas muestras individuales, dependiendo de la forma exacta de la distribución multivariada: algunas posibilidades son los métodos de muestreo de rechazo adaptativo,^[1]​^[2]​^[3]​^[4]​ el algoritmo de muestreo Metrópolis de rechazo adaptativo^[5]​ o sus mejoras,^[6]​^[7]​ Un simple paso unidimensional de Metrópolis-Hastings, o muestreo de corte.

 

El resultado de tres cadenas de Markov que se ejecutan en la función 3D Rosenbrock utilizando el algoritmo Metropolis-Hastings. El algoritmo toma muestras de regiones donde la probabilidad posterior es alta, y las cadenas comienzan a mezclarse en estas regiones. La posición aproximada del máximo ha sido iluminada. Tenga en cuenta que los puntos rojos son los que quedan después del proceso de quemado. Los primeros han sido descartados.

Véase también

• Balance detallado
• Algoritmo Genético
• Muestreo de Gibbs
• Métodos de partículas de campo medio
• Algoritmo Langevin ajustado por Metropolis
• Transporte de la luz de Metropolis
• Metrópolis de intentos múltiples
• Temple en paralelo
• Algoritmo de Crank-Nicolson precondicionado
• Filtro de partículas
• Algoritmo de recocido simulado

Referencias

1. Gilks, W. R.; Wild, P. (1 de enero de 1992). «Adaptive Rejection Sampling for Gibbs Sampling». Journal of the Royal Statistical Society. Series C (Applied Statistics) 41 (2): 337-348. JSTOR 2347565. doi:10.2307/2347565. 
2. Görür, Dilan; Teh, Yee Whye (1 de enero de 2011). «Concave-Convex Adaptive Rejection Sampling». Journal of Computational and Graphical Statistics 20 (3): 670-691. ISSN 1061-8600. doi:10.1198/jcgs.2011.09058. 
3. Hörmann, Wolfgang (1 de junio de 1995). «A Rejection Technique for Sampling from T-concave Distributions». ACM Trans. Math. Softw. 21 (2): 182-193. ISSN 0098-3500. doi:10.1145/203082.203089.  Parámetro desconocido |citeseerx= ignorado (ayuda)
4. Martino, Luca; Míguez, Joaquín (25 de agosto de 2010). «A generalization of the adaptive rejection sampling algorithm». Statistics and Computing 21 (4): 633-647. ISSN 0960-3174. doi:10.1007/s11222-010-9197-9.  Parámetro desconocido |hdl= ignorado (ayuda)
5. Gilks, W. R.; Best, N. G.; Tan, K. K. C. (1 de enero de 1995). «Adaptive Rejection Metropolis Sampling within Gibbs Sampling». Journal of the Royal Statistical Society. Series C (Applied Statistics) 44 (4): 455-472. JSTOR 2986138. doi:10.2307/2986138. 
6. Martino, L.; Read, J.; Luengo, D. (1 de junio de 2015). «Independent Doubly Adaptive Rejection Metropolis Sampling Within Gibbs Sampling». IEEE Transactions on Signal Processing 63 (12): 3123-3138. Bibcode:2015ITSP...63.3123M. ISSN 1053-587X. arXiv:1205.5494. doi:10.1109/TSP.2015.2420537. 
7. Meyer, Renate; Cai, Bo; Perron, François (15 de marzo de 2008). «Adaptive rejection Metropolis sampling using Lagrange interpolation polynomials of degree 2». Computational Statistics & Data Analysis 52 (7): 3408-3423. doi:10.1016/j.csda.2008.01.005. 

•   Datos: Q910810