stringtranslate.com

Algoritmo Metropolis-Hastings

El algoritmo Metropolis-Hastings muestrea 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 de una distribución de probabilidad de la que el muestreo directo es difícil. Se añaden nuevas muestras a la secuencia en dos pasos: primero se propone una nueva muestra basada en la muestra anterior, luego la muestra propuesta se añade 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 el muestreo de distribuciones multidimensionales, especialmente cuando el número de dimensiones es elevado. Para distribuciones unidimensionales, normalmente existen otros métodos (por ejemplo, el muestreo por rechazo adaptativo ) que pueden devolver directamente muestras independientes de la distribución y que están libres del problema de las muestras autocorrelacionadas que es inherente a los métodos MCMC.

Historia

El algoritmo recibe su nombre en parte de Nicholas Metropolis , el primer coautor de un artículo de 1953 titulado Equation of State Calculations by Fast Computing Machines , junto con Arianna W. Rosenbluth , Marshall Rosenbluth , Augusta H. Teller y Edward Teller . Durante muchos años, el algoritmo se conoció simplemente como el algoritmo Metropolis . [1] [2] El artículo propuso el algoritmo para el caso de distribuciones de propuestas simétricas, pero en 1970, W. K. Hastings lo extendió al caso más general. [3] El método generalizado finalmente se identificó con ambos nombres, aunque no está claro el primer uso del término "algoritmo Metropolis-Hastings".

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, había acuñado el término "Monte Carlo" en un artículo anterior con Stanisław Ulam , 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 un relato detallado del desarrollo del algoritmo. Poco antes de su muerte, Marshall Rosenbluth asistió a una conferencia en 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". [4] Gubernatis hace una aclaración histórica adicional en un artículo de revista de 2005 [5] que relata la conferencia del 50 aniversario. Rosenbluth deja claro que él y su esposa Arianna hicieron el trabajo y que Metropolis no jugó ningún papel en el desarrollo más allá de 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] En contraste, el relato detallado de Rosenbluth atribuye a Teller una sugerencia crucial pero temprana de "aprovechar la mecánica estadística y tomar promedios de conjunto 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 John Von Neumann . Arianna Rosenbluth contó (a Gubernatis en 2003) que Augusta Teller comenzó el trabajo de la computadora, pero que Arianna misma lo tomó el control y escribió el código desde cero. En una historia oral registrada poco antes de su muerte, [7] Rosenbluth nuevamente atribuye a Teller el planteamiento del problema original, a él mismo el haberlo resuelto y a Arianna el haber programado 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 que se puedan calcular los valores de . El requisito de que solo 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, lo 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 iterativamente de tal manera que la distribución de la siguiente muestra depende solo del valor de 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 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 siguiente iteración, o es rechazado, 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 muestra actuales y candidatos 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 para es una distribución gaussiana centrada en , de modo que los puntos más cercanos a tengan más probabilidades de ser visitados a continuación, lo que convierte la secuencia de muestras en un paseo aleatorio gaussiano . En el artículo original de Metropolis et al. (1953), se sugirió que fuera una distribución uniforme limitada a una distancia máxima desde . También son posibles funciones de propuesta más complicadas, como las de Monte Carlo hamiltoniano , Monte Carlo de Langevin o Crank-Nicolson preacondicionado .

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

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

Sea una función que es 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 como primera observación en la muestra y elija una función propuesta . En esta sección, se supone que es simétrica; en otras palabras, debe satisfacer .
  2. Para cada iteración t :
    • Proponga un candidato para la siguiente muestra seleccionándolo de la distribución .
    • Calcular la tasa de aceptación , que se utilizará para decidir si se acepta o rechaza al candidato [b] . Como f es proporcional a la densidad de P , tenemos que .
    • Aceptar o rechazar :
      • Generar un número aleatorio uniforme .
      • Si , entonces acepte al candidato estableciendo ,
      • Si , entonces rechaza al candidato y establece en su lugar.

Este algoritmo procede intentando moverse aleatoriamente por el espacio muestral, a veces aceptando los movimientos y a veces permaneciendo en el mismo lugar. Nótese que la razón de aceptación indica qué tan probable es la nueva muestra propuesta con respecto a la muestra actual, de acuerdo con 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 de 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 en la probabilidad, más probable será que rechacemos el nuevo punto. Por lo tanto, tenderemos a permanecer en (y devolver un gran número de muestras de) regiones de alta densidad de , mientras que solo visitaremos 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 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 otra parte, la mayoría de 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 MCMC, no tienen este problema en tal grado, y por lo tanto son a menudo las únicas soluciones disponibles cuando el número de dimensiones de la distribución a muestrear 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 utilizados actualmente en muchas disciplinas.

En distribuciones multivariadas , el algoritmo clásico de Metropolis-Hastings descrito anteriormente implica la elección de un nuevo punto de muestra multidimensional. Cuando el número de dimensiones es alto, puede resultar difícil encontrar la distribución de saltos adecuada para utilizar, ya que las diferentes dimensiones individuales se comportan de formas muy diferentes y el ancho de salto (véase más arriba) debe ser "justo el adecuado" para todas las dimensiones a la vez para evitar una mezcla excesivamente lenta. Un enfoque alternativo que a menudo funciona mejor en tales situaciones, conocido como muestreo de Gibbs , implica la elección de 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 de un espacio potencialmente de alta dimensión se reducirá a una colección de problemas para muestrear de una pequeña dimensionalidad. [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 solo por un pequeño número de otras variables, como es el caso en la mayoría de los modelos jerárquicos típicos . A continuación, se muestrean las variables individuales una a la vez, y cada variable está condicionada por 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 rechazo adaptativo Metropolis, [11] un simple paso Metropolis-Hastings unidimensional o muestreo por cortes .

Derivación formal

El algoritmo Metropolis-Hastings tiene como objetivo generar una colección de estados según una distribución deseada . Para lograrlo, el algoritmo utiliza un proceso de Markov que alcanza asintóticamente una distribución estacionaria única [ ancla rota ] tal que . [12]

Un proceso de Markov se define de forma única por sus probabilidades de transición , es decir, la probabilidad de pasar de un estado dado a 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 exige que cada transición sea reversible: para cada par de estados , la probabilidad de estar en el estado y pasar al estado debe ser igual a la probabilidad de estar en el estado y pasar al 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 el diseño de 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 equilibrio 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 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 ellas:

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 esta tasa de aceptación de Metropolis , ya sea o y, en cualquier caso, se cumple la condición.

El algoritmo de Metropolis-Hastings puede escribirse de la siguiente manera:

  1. Inicializar
    1. Elija un estado inicial .
    2. Colocar .
  2. Iterar
    1. Generar un estado candidato aleatorio de acuerdo a .
    2. Calcular la probabilidad de aceptación .
    3. Aceptar o rechazar :
      1. generar un número aleatorio uniforme ;
      2. si , entonces acepta el nuevo estado y establece ;
      3. si , entonces rechaza el nuevo estado y copia el estado antiguo hacia adelante .
    4. Incremento : establecido .

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 de manera efectiva depende del número de factores, incluida la relación entre y la distribución 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 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 de Metropolis-Hastings es calcular una integral. En concreto, se considera un espacio y una distribución de probabilidad sobre , . Metropolis-Hastings puede estimar una integral de la forma

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 . Suponga 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 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 más probabilidad y así aumentar el número de muestras utilizadas para estimar en las colas. Esto se puede hacer, por ejemplo, utilizando una distribución de muestreo para favorecer a esos estados (por ejemplo, con ).

Instrucciones paso a paso

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

Supongamos que el valor más reciente muestreado 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 razón de probabilidad (por ejemplo, posterior bayesiana) entre la muestra propuesta y la muestra anterior , y

es la relación entre la densidad de propuestas en dos direcciones (de a y viceversa). Es igual a 1 si la densidad de propuestas 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 se "olvida" este estado inicial. Estas muestras, que se descartan, 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 propuestas 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 propuestas gaussiana, el parámetro de varianza debe ajustarse durante el período de rodaje. 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 de alrededor del 50%, disminuyendo a alrededor del 23% para una distribución objetivo gaussiana bidimensional. [15] Estas pautas pueden funcionar bien cuando se muestrea a partir de posteriores bayesianos suficientemente regulares, ya que a menudo siguen una distribución normal multivariante como se puede establecer utilizando el teorema de Bernstein-von Mises . [16]

Si es demasiado pequeño, 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 a ). Por otro lado, si es demasiado grande, la tasa de aceptación será muy baja porque es probable que las propuestas aterricen en regiones de densidad de probabilidad mucho menor, por lo que serán muy pequeñas y, nuevamente, la cadena convergerá muy lentamente. Por lo general, se ajusta la distribución de propuestas para que los algoritmos acepten alrededor del 30% de todas las muestras, de acuerdo con las estimaciones teóricas mencionadas en el párrafo anterior.

Inferencia bayesiana

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

Véase también

Referencias

  1. ^ Kalos, Malvin H.; Whitlock, Paula A. (1986). Métodos de Monte Carlo, volumen I: Fundamentos . Nueva York: Wiley. págs. 78–88.
  2. ^ Tierney, Luke (1994). "Cadenas de Markov para explorar distribuciones posteriores". Anales de 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. Bibcode :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 mecánica estadística". Actas de la conferencia AIP . 690 : 22–30. Código Bibliográfico :2003AIPC..690...22R. doi :10.1063/1.1632112.
  5. ^ JE Gubernatis (2005). "Marshall Rosenbluth y el algoritmo Metropolis". Física de plasmas . 12 (5): 057303. Bibcode :2005PhPl...12e7303G. doi :10.1063/1.1887186.
  6. ^ Teller, Edward. Memorias: un viaje del siglo XX por la ciencia y la política . Perseus Publishing , 2001, pág. 328
  7. ^ Rosenbluth, Marshall. "Transcripción de la historia oral". Instituto Americano de Física
  8. ^ ab Gilks, WR; Wild, P. (1992-01-01). "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.
  9. ^ Análisis de datos bayesianos . Gelman, Andrew (2.ª ed.). Boca Raton, Fla.: Chapman & Hall / CRC. 2004. ISBN 978-1584883883.OCLC 51991499  .{{cite book}}: CS1 maint: others (link)
  10. ^ Lee, Se Yoon (2021). "Inferencia variacional mediante el muestreador de Gibbs y el 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; Best, NG ; Tan, KKC (1 de enero de 1995). "Muestreo de metrópolis con 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.
  12. ^ de Robert, Christian; Casella, George (2004). Métodos estadísticos de Monte Carlo. Springer. ISBN 978-0387212395.
  13. ^ Raftery, Adrian E. y Steven Lewis. "¿Cuántas iteraciones hay en el muestreador de Gibbs?", en Bayesian Statistics 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, GO; Gelman, A.; Gilks, WR (1997). "Convergencia débil y escalamiento óptimo de algoritmos de recorrido aleatorio de Metropolis". Ann. Appl. Probab. 7 (1): 110–120. CiteSeerX 10.1.1.717.2582 . doi :10.1214/aoap/1034625254.  
  16. ^ Schmon, Sebastian M.; Gagnon, Philippe (15 de abril de 2022). "Escalamiento óptimo de algoritmos Metropolis de recorrido aleatorio utilizando asintóticos bayesianos 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 tomó como 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), en realidad era la distribución de Boltzmann , tal como se aplicaba a los 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 dada en equilibrio térmico). En consecuencia, la razón de aceptación era en sí misma un exponencial de la diferencia en los parámetros del numerador y el denominador de esta razón.

Lectura adicional