En estadística y 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 de la que el muestreo directo es difícil. Esta secuencia se puede utilizar 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 MCMC se utilizan generalmente para muestrear distribuciones multidimensionales, especialmente cuando el número de dimensiones es alto. Para distribuciones unidimensionales, generalmente existen 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 MCMC.
Historia
El algoritmo lleva el nombre de Nicholas Metropolis , autor del artículo de 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 WK Hastings lo extendió al caso más general en 1970. [1]
Existe cierta controversia con respecto al crédito por el desarrollo del algoritmo. Metropolis, que estaba familiarizado con los aspectos computacionales del método, había acuñado el término "Monte Carlo" en un artículo anterior con Stanisław Ulam , y dirigió el grupo de 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 descripción detallada del desarrollo del algoritmo. Luego, poco antes de su muerte, Marshall Rosenbluth asistió a una conferencia de 2003 en LANL que marcaba 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". [2] Gubernatis hace una aclaración histórica adicional en un artículo de revista de 2005 [3] que relata la conferencia del 50º aniversario. Rosenbluth deja en claro que él y su esposa Arianna hicieron el trabajo, y que Metropolis no jugó ningún papel en el desarrollo más que proporcionar tiempo de computadora.
Esto contradice un relato de Edward Teller, quien afirma en sus memorias que los cinco autores del artículo de 1953 trabajaron juntos durante "días (y noches)". [4] En contraste, el relato detallado de Rosenbluth acredita a Teller con una sugerencia crucial pero temprana de "aprovechar la mecánica estadística y tomar promedios de conjuntos en lugar de seguir la cinemática detallada". Esto, dice Rosenbluth, lo hizo pensar en el enfoque generalizado de Monte Carlo, un tema que dice haber discutido a menudo con Von Neumann . Arianna Rosenbluth contó (a Gubernatis en 2003) que Augusta Teller inició el trabajo con 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, [5] Rosenbluth nuevamente le da crédito a Teller por plantear el problema original, a él mismo por resolverlo ya Arianna por programar la computadora. En términos de reputación, hay pocas razones para cuestionar el relato de Rosenbluth. En una memoria biográfica de Rosenbluth, Freeman Dyson escribe: [6]
Muchas veces vine a Rosenbluth, le hice una pregunta [...] y recibí una respuesta en dos minutos. Entonces, 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 habilidad 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 era 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 Metropolis-Hastings puede extraer muestras de cualquier distribución de probabilidad , siempre que conozcamos una función proporcional a la densidad de y los valores de se puede calcular. El requisito de que debe ser solo proporcional a la densidad, en lugar de exactamente igual a ella, hace que el algoritmo Metropolis-Hastings sea particularmente útil, porque calcular el factor de normalización necesario es a menudo 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. . Estos valores de muestra se producen de forma iterativa, y la distribución de la siguiente muestra depende únicamente del valor de muestra actual (lo que convierte 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 basándose en el valor de muestra actual. Luego, con cierta probabilidad, el candidato se acepta (en cuyo caso el valor candidato se usa en la siguiente iteración) o se rechaza (en cuyo caso el valor candidato se descarta y el valor actual se reutiliza en la siguiente iteración): la probabilidad de aceptación se determina comparando los valores de la función de los valores muestrales actuales y candidatos con respecto a la distribución deseada .
A modo de ilustración, a continuación se describe el algoritmo Metropolis, un caso especial del algoritmo Metropolis-Hastings donde la función propuesta es simétrica.
Algoritmo de Metropolis (distribución simétrica de propuestas)
Dejar ser una función que es proporcional a la distribución de probabilidad deseada (también conocido como una distribución de destino).
- Inicialización: elija un punto arbitrario ser la primera muestra y elegir una densidad de probabilidad arbitraria (a veces escrito ) que sugiere un candidato para el siguiente valor de muestra , dado el valor de la muestra anterior . En esta sección,se supone que es simétrico; en otras palabras, debe satisfacer. Una elección habitual es dejarser una distribución gaussiana centrada en, de modo que apunte más cerca de es más probable que se visiten a continuación, lo que convierte la secuencia de muestras en un recorrido aleatorio [a] . La funciónse denomina densidad de propuesta o distribución de salto .
- Para cada iteración t :
- Genera un candidato para la siguiente muestra seleccionando de la distribución .
- Calcule la tasa de aceptación , que se utilizará para decidir si acepta o rechaza al candidato. Como f es proporcional a la densidad de P , tenemos que.
- Aceptar o rechazar :
- Genera un número aleatorio uniforme .
- Si , luego acepte al candidato configurando,
- Si , luego rechace al candidato y establezca en lugar de.
Este algoritmo procede al intentar moverse aleatoriamente por el espacio muestral, a veces aceptando los movimientos y, a veces, permaneciendo en su lugar. Tenga en cuenta que la tasa de aceptación indica qué tan probable es la nueva muestra propuesta con respecto a la muestra actual, de acuerdo con la distribución . Si intentamos movernos a un punto que es más probable que el punto existente (es decir, un punto en una región de mayor densidad de correspondiente a un ), siempre aceptaremos la mudanza. Sin embargo, si intentamos movernos a un punto menos probable, a veces rechazaremos el movimiento, y cuanto mayor sea la caída relativa en la probabilidad, más probable será que rechacemos el nuevo punto. Por lo tanto, tendremos a permanecer en (y devolver un gran número de muestras de) regiones de alta densidad de, mientras que solo visita ocasionalmente regiones de baja densidad. Intuitivamente, esta es la razón por la que este algoritmo funciona y devuelve muestras que siguen la distribución deseada..
En comparación con un algoritmo como el muestreo de rechazo adaptativo [7] que genera directamente muestras independientes a partir de una distribución, Metropolis-Hastings y otros algoritmos de MCMC tienen una serie de desventajas:
- Las muestras están correlacionadas. Aunque a largo plazo sigan correctamente, un conjunto de muestras cercanas se correlacionarán entre sí y no reflejarán correctamente la distribución. Esto significa que si queremos un conjunto de muestras independientes, que tienen que tirar la mayoría de las muestras y sólo tomar cada n ésima muestra, para un cierto valor de n (normalmente determinado por el examen de la auto-correlación entre muestras adyacentes). La autocorrelación se puede reducir 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 dará lugar a una cadena de Markov de mezcla lenta , es decir, un conjunto de muestras altamente correlacionadas, de modo que se necesitará una gran cantidad 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, normalmente es necesario un período de quemado , [8] en el que se desecha un número inicial de muestras (por ejemplo, las primeras 1000 aproximadamente).
Por otro lado, la mayoría de los métodos de muestreo de rechazo 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 MCMC, 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 que se muestrea es alto. Como resultado, los métodos MCMC 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 que se utilizan hoy en día en muchas disciplinas.
En distribuciones multivariadas , el algoritmo clásico de Metropolis-Hastings como se describe arriba implica elegir un nuevo punto de muestra multidimensional. Cuando el número de dimensiones es alto, puede resultar difícil encontrar la distribución de salto adecuada para usar, ya que las diferentes dimensiones individuales se comportan de formas muy diferentes, y el ancho de salto (ver arriba) debe ser "justo" para todas las dimensiones a la vez para Evite una mezcla excesivamente lenta. 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. De esa manera, el problema del muestreo de un espacio potencialmente de alta dimensión se reducirá a una colección de problemas para muestrear a partir de una dimensión pequeña. [9] Esto es especialmente aplicable cuando la distribución multivariante se compone de un conjunto de variables aleatorias individuales en las que cada variable está condicionada a sólo un pequeño número de otras variables, como es el caso en la mayoría de los modelos jerárquicos típicos . Luego, las variables individuales se muestrean una por 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 multivariante: algunas posibilidades son los métodos de muestreo de rechazo adaptativo , [7] el algoritmo de muestreo de rechazo adaptativo de Metropolis, [10] un simple Metropolis unidimensional– Paso de Hastings o muestreo por rebanadas .
Derivación formal
El propósito del algoritmo Metropolis-Hastings es generar una colección de estados de acuerdo con una distribución deseada . Para lograr esto, el algoritmo utiliza un proceso de Markov , que alcanza asintóticamente una distribución estacionaria única tal que . [11]
Un proceso de Markov se define de forma única por sus probabilidades de transición , la probabilidad de pasar de cualquier estado dado a cualquier otro estado dado . Tiene una distribución estacionaria únicacuando se cumplen las dos condiciones siguientes: [11]
- Existencia de distribución estacionaria : debe existir una distribución estacionaria. Una condición suficiente pero no necesaria es el equilibrio detallado , que requiere que cada transición es reversible: para cada par de estados , la probabilidad de estar en estado y la transición al estado debe ser igual a la probabilidad de estar en estado y la transición al estado , .
- Unicidad de la distribución estacionaria : la distribución estacionariadebe ser único. Esto está garantizado por la ergodicidad del proceso de Markov, que requiere que cada estado debe (1) ser aperiódico: el sistema no regresa al mismo estado a intervalos fijos; y (2) ser positivo recurrente: el número esperado de pasos para volver al mismo estado es finito.
El algoritmo Metropolis-Hastings implica diseñar un proceso de Markov (mediante la construcción de probabilidades de transición) que cumpla las dos condiciones anteriores, de modo que su distribución estacionaria es elegido para ser . La derivación del algoritmo comienza con la condición de equilibrio detallado :
que está reescrito como
El enfoque consiste en separar la transición en dos subpasos; la propuesta y la aceptación-rechazo. La distribución de la propuesta es la probabilidad condicional de proponer un estado dado , y la distribución de aceptación es la probabilidad de aceptar el estado propuesto . La probabilidad de transición se puede escribir como el producto de ellos:
Insertando esta relación en la ecuación anterior, tenemos
El siguiente paso en la derivación es elegir un índice de aceptación que cumpla con la condición anterior. Una opción común es la opción Metropolis:
Para esta tasa de aceptación de Metropolis , ya sea o y, de cualquier manera, se cumple la condición.
El algoritmo Metropolis-Hastings, por tanto, consiste en lo siguiente:
- Inicializar
- Elige un estado inicial .
- Colocar .
- Iterar
- Genera un estado candidato aleatorio de acuerdo a .
- Calcule la probabilidad de aceptación.
- Aceptar o rechazar :
- generar un número aleatorio uniforme ;
- Si , luego acepte el nuevo estado y configure;
- Si , luego rechace el nuevo estado y copie el estado anterior hacia adelante.
- Incremento : establecer.
Siempre que se cumplan las condiciones especificadas, la distribución empírica de los estados guardados se acercará . El número de iteraciones () necesarios para estimar de forma eficaz depende del número de factores, incluida la relación entre y la distribución de la propuesta y la precisión deseada de la estimación. [12] Para la distribución en espacios de estados discretos, tiene que ser del orden del tiempo de autocorrelación del proceso de Markov. [13]
Es importante notar que no está claro, en un problema general, qué distribución se debe utilizar o el número de iteraciones necesarias para una estimación adecuada; ambos son parámetros libres del método, que deben ajustarse al problema particular en cuestión.
Uso en integración numérica
Un uso común del algoritmo Metropolis-Hastings es calcular una integral. Específicamente, considere un espacio y una distribución de probabilidad encima , . Metropolis-Hastings puede estimar una integral de la forma de
dónde es una función (medible) de interés.
Por ejemplo, considere una estadística y su distribución de probabilidad , que es una distribución marginal . Suponga que el objetivo es estimar por en la cola de . Formalmente, Se puede escribir como
y, por tanto, estimar se puede lograr estimando el valor esperado de la función del indicador , que es 1 cuando y cero en caso contrario. Porque está en la cola de , la probabilidad de dibujar un estado con en la cola de es proporcional a , que es pequeño por definición. El algoritmo Metropolis-Hastings se puede utilizar aquí para muestrear estados (raros) con mayor probabilidad y, por lo tanto, aumentar la cantidad de muestras utilizadas para estimaren las colas. Esto se puede hacer, por ejemplo, utilizando una distribución de muestreo. para favorecer a esos estados (p. ej. con ).
Instrucciones paso a paso
Suponga que el valor más reciente muestreado es . Para seguir el algoritmo Metropolis-Hastings, dibujamos a continuación un nuevo estado de propuesta con densidad de probabilidad y calcula un valor
dónde
es la razón de probabilidad (por ejemplo, posterior bayesiano) entre la muestra propuesta y la muestra anterior , y
es la relación de la densidad de la propuesta en dos direcciones (de a y por el contrario). Esto es igual a 1 si la densidad de la propuesta es simétrica. Entonces el nuevo estado se elige de acuerdo con las siguientes reglas.
- Si
- demás:
La cadena de Markov se inicia a partir de un valor inicial arbitrario , y el algoritmo se ejecuta durante muchas iteraciones hasta que este estado inicial se "olvida". Estas muestras, que se descartan, se conocen como quemaduras . El conjunto restante de valores aceptados derepresentar una muestra de la distribución.
El algoritmo funciona mejor si la densidad de la propuesta coincide con la forma de la distribución de destino. , de la cual el muestreo directo es difícil, es decir . Si una propuesta de densidad gaussiana se utiliza, el parámetro de varianza tiene que ser sintonizado durante el período de quemado. Esto generalmente se hace calculando la tasa de aceptación , que es la fracción de muestras propuestas que se acepta en una ventana de la últimamuestras. La tasa de aceptación deseada depende de la distribución objetivo; sin embargo, se ha demostrado teóricamente que la tasa de aceptación ideal para una distribución gaussiana unidimensional es aproximadamente del 50%, disminuyendo a aproximadamente el 23% para una distribución gaussiana unidimensional.-distribución objetivo gaussiana dimensional. [14]
Si es demasiado pequeña, la cadena se mezclará lentamente (es decir, la tasa de aceptación será alta, pero las muestras sucesivas se moverán lentamente por el espacio y la cadena convergerá solo lentamente para). Por otro lado, si es demasiado grande, la tasa de aceptación será muy baja porque es probable que las propuestas lleguen a regiones con una densidad de probabilidad mucho menor, por lo que será muy pequeño, y nuevamente la cadena convergerá muy lentamente. Por lo general, se ajusta la distribución de la propuesta para que los algoritmos acepten del orden del 30% de todas las muestras, de acuerdo con las estimaciones teóricas mencionadas en el párrafo anterior.
Ver también
- Saldo detallado
- Algoritmos genéticos
- Muestreo de Gibbs
- Métodos de partículas de campo medio
- Algoritmo de Langevin ajustado por Metropolis
- Transporte ligero Metropolis
- Metrópolis de múltiples intentos
- Templado paralelo
- Algoritmo Crank-Nicolson preacondicionado
- Montecarlo secuencial
- Recocido simulado
Referencias
- ^ Hastings, WK (1970). "Métodos de muestreo de Monte Carlo utilizando cadenas de Markov y sus aplicaciones". Biometrika . 57 (1): 97–109. Código bibliográfico : 1970Bimka..57 ... 97H . doi : 10.1093 / biomet / 57.1.97 . JSTOR 2334940 . Zbl 0219.65008 .
- ^ MN Rosenbluth (2003). "Génesis del algoritmo de Monte Carlo para la mecánica estadística". Actas de la conferencia AIP . 690 : 22-30. doi : 10.1063 / 1.1632112 .
- ^ JE Gubernatis (2005). "Marshall Rosenbluth y el algoritmo de Metropolis" . Física de Plasmas . 12 (5): 057303. Código Bibliográfico : 2005PhPl ... 12e7303G . doi : 10.1063 / 1.1887186 .
- ^ Teller, Edward. Memorias: Un viaje del siglo XX en ciencia y política . Perseus Publishing , 2001, pág. 328
- ^ Rosenbluth, Marshall. "Transcripción de historia oral" . Instituto Americano de Física
- ^ F. Dyson (2006). "Marshall N. Rosenbluth". Actas de la American Philosophical Society . 250 : 404.
- ^ a b Gilks, WR; Wild, P. (1 de enero de 1992). "Muestreo de rechazo adaptativo para muestreo de Gibbs". Revista de la Royal Statistical Society. Serie C (Estadística aplicada) . 41 (2): 337–348. doi : 10.2307 / 2347565 . JSTOR 2347565 .
- ^ Análisis de datos bayesianos . Gelman, Andrew (2ª ed.). Boca Raton, Fla .: Chapman & Hall / CRC. 2004. ISBN 978-1584883883. OCLC 51991499 .CS1 maint: otros ( enlace )
- ^ Lee, Se Yoon (2021). "Inferencia variacional de ascenso y muestreo de Gibbs: una revisión teórica de conjuntos". Comunicaciones en estadística: teoría y métodos . arXiv : 2008.01006 . doi : 10.1080 / 03610926.2021.1921214 .
- ^ Gilks, WR; Mejor, NG ; Tan, KKC (1 de enero de 1995). "Muestreo de metrópolis de rechazo adaptativo dentro del muestreo de Gibbs". Revista de la Royal Statistical Society. Serie C (Estadística aplicada) . 44 (4): 455–472. doi : 10.2307 / 2986138 . JSTOR 2986138 .
- ^ a b Robert, Christian; Casella, George (2004). Métodos estadísticos de Monte Carlo . Saltador. ISBN 978-0387212395.
- ^ Raftery, Adrian E. y Steven Lewis. "¿Cuántas iteraciones en el muestreador de Gibbs?" En Estadísticas Bayesianas 4 . 1992.
- ^ Newman, MEJ; Barkema, GT (1999). Métodos de Monte Carlo en Física Estadística . Estados Unidos: Oxford University Press. ISBN 978-0198517979.
- ^ Roberts, GO; Gelman, A .; Gilks, WR (1997). "Convergencia débil y escalamiento óptimo de algoritmos de Metropolis de paseo aleatorio" . Ana. Apl. Probab. 7 (1): 110–120. CiteSeerX 10.1.1.717.2582 . doi : 10.1214 / aoap / 1034625254 .
Notas
- ^ Sin embargo, en el documento original,era en realidad la distribución de Boltzmann , ya que se aplicó a los sistemas físicos en el contexto de la mecánica estadística
Otras lecturas
- Bernd A. Berg . Simulaciones de Monte Carlo de la cadena de Markov y su análisis estadístico . Singapur, World Scientific , 2004.
- Siddhartha Chib y Edward Greenberg: "Comprensión del algoritmo Metropolis-Hastings". Estadístico estadounidense , 49 (4), 327–335, 1995
- David DL Minh y Do Le Minh. "Comprensión del algoritmo de Hastings". Comunicaciones en estadística: simulación y computación, 44: 2332-349, 2015
- Bolstad, William M. (2010) Comprensión de las estadísticas bayesianas computacionales , John Wiley & SonsISBN 0-470-04609-0