stringtranslate.com

Algoritmo de Metrópolis-Hastings

El algoritmo de Metropolis-Hastings muestra una distribución de probabilidad posterior unidimensional normal .

En estadística y física estadística , el algoritmo Metropolis-Hastings es un método de Monte Carlo de cadena de Markov (MCMC) para obtener una secuencia de muestras aleatorias a partir de una distribución de probabilidad de la cual el muestreo directo es difícil. Se agregan nuevas muestras a la secuencia en dos pasos: primero se propone una nueva muestra basada en la muestra anterior, luego la muestra propuesta se agrega a la secuencia o se rechaza dependiendo del valor de la distribución de probabilidad en ese punto. La secuencia resultante 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 en parte de Nicholas Metropolis , el primer coautor de un artículo de 1953, titulado Ecuación de cálculos de estado mediante máquinas de computación rápida , con Arianna W. Rosenbluth , Marshall Rosenbluth , Augusta H. Teller y Edward Teller . Durante muchos años el algoritmo fue conocido simplemente como algoritmo de Metropolis . [1] [2] El artículo propuso el algoritmo para el caso de distribuciones de propuestas simétricas, pero en 1970, WK Hastings lo amplió al caso más general. [3] El método generalizado finalmente fue identificado con ambos nombres, aunque el primer uso del término "algoritmo de Metropolis-Hastings" no está claro.

Existe cierta controversia con respecto al crédito por el desarrollo del algoritmo Metropolis. Metropolis, que estaba familiarizado con los aspectos computacionales del método, acuñó 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ó el ordenador MANIAC I utilizado en los experimentos en 1952. Sin embargo, antes de 2003 no existía una descripción detallada del desarrollo del algoritmo. Poco antes de su muerte, Marshall Rosenbluth asistió a una conferencia en LANL en 2003 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 mecánica estadística". [4] Gubernatis hace más aclaraciones históricas en un artículo de revista de 2005 [5] 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)". [6] Por el contrario, el relato detallado de Rosenbluth le da crédito a Teller por una sugerencia crucial pero temprana de "aprovechar la mecánica estadística y tomar promedios conjuntos 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 dice, había discutido a menudo con John 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, [7] Rosenbluth nuevamente le da crédito a Teller por plantear el problema original, a él mismo por resolverlo y a Arianna por programar la computadora.

Descripción

El algoritmo de Metropolis-Hastings puede extraer muestras de cualquier distribución de probabilidad con densidad de probabilidad , siempre que conozcamos una función proporcional a la densidad y se puedan calcular los valores de . El requisito de que sólo debe ser proporcional a la densidad, en lugar de exactamente igual a ella, hace que el algoritmo de Metropolis-Hastings sea particularmente útil, porque elimina la necesidad de calcular el factor de normalización de la densidad, que a menudo es extremadamente difícil en la práctica.

El algoritmo Metropolis-Hastings genera 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 de tal manera que la distribución de la siguiente muestra depende sólo del valor de la muestra actual, lo que hace que la secuencia de muestras sea una cadena de Markov . Específicamente, en cada iteración, el algoritmo propone un candidato para el siguiente valor de muestra en función del 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 de la muestra actual y candidata con respecto a la distribución deseada.

El método utilizado para proponer nuevos candidatos se caracteriza por la distribución de probabilidad (a veces escrita ) de una nueva muestra propuesta dada la muestra anterior . Esto se llama densidad de propuesta , función de propuesta o distribución de salto . Una opción común es una distribución gaussiana centrada en , de modo que es más probable que los puntos más cercanos sean visitados a continuación, lo que convierte la secuencia de muestras en una caminata aleatoria gaussiana . En el artículo original de Metropolis et al. (1953), se sugirió que era una distribución uniforme limitada a una distancia máxima desde . También son posibles funciones de propuesta más complicadas, como las de Hamiltonian Monte Carlo , Langevin Monte Carlo o Crank-Nicolson precondicionadas .

Con fines ilustrativos, 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 de propuestas simétricas)

Sea una función que sea proporcional a la función de densidad de probabilidad deseada (también conocida como distribución objetivo) [a] .

  1. Inicialización: elija un punto arbitrario para que sea la primera observación en la muestra y elija una función de propuesta . En esta sección, se supone que es simétrico; en otras palabras, debe satisfacer .
  2. Para cada iteración t :
    • Proponga un candidato para la siguiente muestra eligiendo de la distribución .
    • Calcule el índice de aceptación , que se utilizará para decidir si aceptar o rechazar al candidato [b] . Como f es proporcional a la densidad de P , tenemos eso .
    • Aceptar o rechazar :
      • Genera un número aleatorio uniforme .
      • Si , entonces acepte al candidato configurando ,
      • Si es así , rechace al candidato y configúrelo en su lugar.

Este algoritmo procede intentando moverse aleatoriamente por el espacio muestral, a veces aceptando los movimientos y otras veces permaneciendo en el lugar. Nótese que el índice de aceptación indica qué tan probable es la nueva muestra propuesta con respecto a la muestra actual, según la distribución cuya densidad es . 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 correspondiente a un ), 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 probabilidades tendremos de rechazar el nuevo punto. Por lo tanto, tenderemos a permanecer en (y devolver una gran cantidad de muestras de) regiones de alta densidad de , mientras que solo visitaremos ocasionalmente regiones de baja densidad. Intuitivamente, es por eso que este algoritmo funciona y devuelve muestras que siguen la distribución deseada con densidad .

En comparación con un algoritmo como el muestreo de rechazo adaptativo [8] que genera directamente muestras independientes a partir de una distribución, Metropolis-Hastings y otros algoritmos MCMC tienen una serie de desventajas:

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 van a muestrear es alto. Como resultado, los métodos MCMC son a menudo los métodos elegidos 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 describe anteriormente, 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, ya que las diferentes dimensiones individuales se comportan de maneras muy diferentes y el ancho de salto (ver arriba) debe ser "perfecto" 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. De esa manera, el problema de muestrear desde un espacio potencialmente de alta dimensión se reducirá a una colección de problemas para muestrear desde una dimensionalidad pequeña. [10] 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 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 se muestrean las variables individuales una a la vez, condicionando cada variable a los valores más recientes de todas las demás. Se pueden utilizar 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 , [8] el algoritmo de muestreo de Metropolis de rechazo adaptativo, [11] un Metropolis unidimensional simple. "Paso de Hastings, o muestreo de cortes" .

Derivación formal

El propósito del algoritmo Metropolis-Hastings es generar una colección de estados según una distribución deseada . Para lograr esto, el algoritmo utiliza un proceso de Markov , que asintóticamente alcanza una distribución estacionaria única tal que . [12]

Un proceso de Markov se define únicamente por sus probabilidades de transición , la probabilidad de pasar de cualquier estado dado a cualquier otro estado dado . Tiene una distribución estacionaria única cuando se cumplen las dos condiciones siguientes: [12]

  1. 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 sea reversible: para cada par de estados , la probabilidad de estar en un estado y pasar a un estado debe ser igual a la probabilidad de estar en un estado y pasar a un estado .
  2. Unicidad de la distribución estacionaria : la distribución estacionaria debe ser única. 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 recurrente positivo: el número esperado de pasos para regresar al mismo estado es finito.

El algoritmo de 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 se elija que su distribución estacionaria sea . La derivación del algoritmo comienza con la condición de saldo detallado :

que se reescribe como

El enfoque consiste en separar la transición en dos subpasos; la propuesta y la aceptación-rechazo. La distribución de 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 la condición anterior. Una opción común es la opción Metropolis:

Para este índice de aceptación de Metropolis , ya sea o y, de cualquier manera, se cumple la condición.

Por tanto, el algoritmo de Metropolis-Hastings se puede escribir de la siguiente manera:

  1. Inicializar
    1. Elija un estado inicial .
    2. Colocar .
  2. Iterar
    1. Genere un estado candidato aleatorio según .
    2. Calcule la probabilidad de aceptación .
    3. Aceptar o rechazar :
      1. generar un número aleatorio uniforme ;
      2. if , entonces acepte el nuevo estado y configúrelo ;
      3. Si es así , rechace el nuevo estado y copie el estado anterior hacia adelante .
    4. Incremento : establecer .

Siempre que se cumplan las condiciones especificadas, la distribución empírica de los estados guardados se aproximará a . El número de iteraciones ( ) necesarias para estimar eficazmente depende del número de factores, incluida la relación entre la distribución de la propuesta y la precisión deseada de la estimación. [13] Para la distribución en espacios de estados discretos, tiene que ser del orden del tiempo de autocorrelación del proceso de Markov. [14]

Es importante notar que no está claro, en un problema general, qué distribución se debe usar 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 de Metropolis-Hastings es calcular una integral. Específicamente, considere un espacio y una distribución de probabilidad sobre , . Metropolis-Hastings puede estimar una integral de la forma de

donde es una función de interés (medible).

Por ejemplo, considere una estadística y su distribución de probabilidad , que es una distribución marginal . Supongamos que el objetivo es estimar para en la cola de . Formalmente, se puede escribir como

y, por lo tanto, la estimación se puede lograr estimando el valor esperado de la función indicadora , que es 1 cuando y cero en caso contrario. Debido a que está en la cola de , la probabilidad de dibujar un estado en la cola de es proporcional a , que es pequeña por definición. El algoritmo Metropolis-Hastings se puede utilizar aquí para muestrear estados (raros) con mayor probabilidad y así aumentar el número de muestras utilizadas para estimar las colas. Esto se puede hacer, por ejemplo, utilizando una distribución muestral para favorecer a esos estados (por ejemplo, con ).

Instrucciones paso a paso

Tres cadenas de Markov que se ejecutan en la función 3D Rosenbrock utilizando el algoritmo Metropolis-Hastings. Las cadenas convergen y se mezclan en la región donde la función es alta. Se ha iluminado la posición aproximada del máximo. Los puntos rojos son los que quedan después del proceso de quemado. Los anteriores han sido descartados.

Supongamos que el valor muestreado más reciente es . Para seguir el algoritmo de Metropolis-Hastings, a continuación dibujamos un nuevo estado de propuesta con densidad de probabilidad y calculamos un valor

dónde

es la relación de probabilidad (por ejemplo, posterior bayesiana) 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 viceversa). Esto es igual a 1 si la densidad de la propuesta es simétrica. Luego se elige el nuevo estado 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 se "olvida" este estado inicial. Estas muestras, que son desechadas, se conocen como burn-in . El conjunto restante de valores aceptados de representa una muestra de la distribución .

El algoritmo funciona mejor si la densidad de la propuesta coincide con la forma de la distribución objetivo , de la cual el muestreo directo es difícil, es decir . Si se utiliza una densidad de propuesta gaussiana , el parámetro de varianza debe ajustarse durante el período de preparación. 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 las últimas muestras. 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 el 50 %, disminuyendo a aproximadamente el 23 % para una distribución objetivo gaussiana bidimensional. [15] Estas pautas pueden funcionar bien cuando se toman muestras de posteriores bayesianos suficientemente regulares, ya que a menudo siguen una distribución normal multivariada, como se puede establecer utilizando el teorema de Bernstein-von Mises . [dieciséis]

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 ). 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ña 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.

Inferencia bayesiana

MCMC se puede utilizar para extraer muestras de la distribución posterior de un modelo estadístico. La probabilidad de aceptación está dada por: donde está la probabilidad , la densidad de probabilidad previa y la probabilidad de propuesta (condicional).

Ver también

Referencias

  1. ^ Kalos, Malvin H.; Whitlock, Paula A. (1986). Métodos Monte Carlo Volumen I: Conceptos básicos . Nueva York: Wiley. págs. 78–88.
  2. ^ Tierney, Lucas (1994). "Cadenas de Markov para explorar las distribuciones posteriores". Los anales de la estadística . 22 (4): 1701-1762. doi : 10.1214/aos/1176325750.
  3. ^ 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.
  4. ^ MN Rosenbluth (2003). "Génesis del algoritmo de Monte Carlo para la mecánica estadística". Actas de la conferencia AIP . 690 : 22–30. Código Bib : 2003AIPC..690...22R. doi : 10.1063/1.1632112.
  5. ^ 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.
  6. ^ Cajero, Eduardo. Memorias: un viaje por la ciencia y la política del siglo XX . Editorial Perseo , 2001, pág. 328
  7. ^ Rosenbluth, Marshall. "Transcripción de historia oral". Instituto Americano de Física
  8. ^ ab Gilks, WR; Salvaje, P. (1 de enero de 1992). "Muestreo de rechazo adaptativo para muestreo de Gibbs". Revista de la Real Sociedad de Estadística. Serie C (Estadística Aplicada) . 41 (2): 337–348. doi :10.2307/2347565. JSTOR  2347565.
  9. ^ Análisis de datos bayesianos . Gelman, Andrew (2ª ed.). Boca Ratón, Florida: Chapman & Hall / CRC. 2004.ISBN 978-1584883883. OCLC  51991499.{{cite book}}: CS1 maint: others (link)
  10. ^ Lee, Se Yoon (2021). "Muestreador de Gibbs e inferencia variacional de ascenso de coordenadas: una revisión de la teoría de conjuntos". Comunicaciones en Estadística - Teoría y Métodos . 51 (6): 1549-1568. arXiv : 2008.01006 . doi :10.1080/03610926.2021.1921214. S2CID  220935477.
  11. ^ Gilks, WR; Mejor, NG ; Bronceado, KKC (1 de enero de 1995). "Muestreo de metrópolis de rechazo adaptativo dentro del muestreo de Gibbs". Revista de la Real Sociedad de Estadística. Serie C (Estadística Aplicada) . 44 (4): 455–472. doi :10.2307/2986138. JSTOR  2986138.
  12. ^ ab Robert, cristiano; Casella, George (2004). Métodos estadísticos de Montecarlo. Saltador. ISBN 978-0387212395.
  13. ^ Raftery, Adrian E. y Steven Lewis. "¿Cuántas iteraciones hay en el muestreador de Gibbs?" En Estadística bayesiana 4 . 1992.
  14. ^ Newman, MEJ; Barkema, GT (1999). Métodos de Monte Carlo en Física Estadística . Estados Unidos: Oxford University Press. ISBN 978-0198517979.
  15. ^ Roberts, VAMOS; Gelman, A.; Gilks, WR (1997). "Convergencia débil y escalado óptimo de algoritmos de paseo aleatorio Metropolis". Ana. Aplica. Probablemente. 7 (1): 110–120. CiteSeerX 10.1.1.717.2582 . doi : 10.1214/aoap/1034625254.  
  16. ^ Schmon, Sebastián M.; Gagnon, Philippe (15 de abril de 2022). "Escalado óptimo de algoritmos de Metropolis de paseo aleatorio utilizando asintóticas bayesianas de muestras grandes". Estadística y Computación . 32 (2): 28. doi :10.1007/s11222-022-10080-8. ISSN  0960-3174. PMC 8924149 . PMID  35310543. 

Notas

  1. ^ En el artículo original de Metropolis et al. (1953), se consideró la distribución de Boltzmann ya que la aplicación específica considerada fue la integración de Monte Carlo de ecuaciones de estado en química física ; la extensión de Hastings se generalizó a una distribución arbitraria .
  2. ^ En el artículo original de Metropolis et al. (1953), era en realidad la distribución de Boltzmann , tal como se aplicaba a sistemas físicos en el contexto de la mecánica estadística (por ejemplo, una distribución de entropía máxima de microestados para una temperatura determinada en equilibrio térmico). En consecuencia, el índice de aceptación era en sí mismo un exponencial de la diferencia en los parámetros del numerador y denominador de este índice.

Otras lecturas