Algoritmo de Metropolis-Hastings

Distribución de la propuesta Q propone el siguiente punto al que el paseo al azar podría moverse.

En estadística y física estadística, el algoritmo Metropolis-Hastings es un método de cadena de Markov Monte Carlo (MCMC) para obtener una secuencia de muestras aleatorias a partir de una distribución de probabilidad a partir de la cual 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 el muestreo de distribuciones multidimensionales, especialmente cuando el número de dimensiones es alto. Para distribuciones unidimensionales, por lo general existen otros métodos (por ejemplo, muestreo de rechazo adaptativo) que pueden devolver directamente muestras independientes de la distribución, y no tienen el problema de las muestras autocorrelacionadas que es inherente a los métodos MCMC.

Historia

El algoritmo lleva el nombre de Nicholas Metropolis y W.K. Hastings. Metropolis fue el primer autor en aparecer en la lista de autores 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 durante muchos años se conoció como el "algoritmo Metropolis". En 1970, Hastings lo extendió al caso más general. El método generalizado finalmente recibió ambos nombres, aunque el primer uso del término "algoritmo Metropolis-Hastings" no esta claro; puede haber aparecido por primera vez en una revisión de 1995 de Chib y Greenberg.

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 una descripción detallada del desarrollo del algoritmo. Poco antes de su muerte, Marshall Rosenbluth asistió a una conferencia de 2003 en LANL con motivo del 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". 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 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)'. Por el contrario, el relato detallado de Rosenbluth le da crédito a Teller con 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 dice que había discutido a menudo con John Von Neumann. Arianna Rosenbluth contó (a Gubernatis en 2003) que Augusta Teller inició el trabajo informático, pero que Arianna misma se hizo cargo y escribió el código desde cero. En una historia oral registrada poco antes de su muerte, 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.

Intuición

El algoritmo de Metropolis-Hastings puede extraer muestras de cualquier distribución de probabilidad con densidad de probabilidad ${displaystyle P(x)}$, siempre que conozcamos una función ${displaystyle f(x)}$ proporcional a la densidad ${displaystyle P}$ y los valores de ${displaystyle f(x)}$ se puede calcular. El requisito de que ${displaystyle f(x)}$ debe ser sólo proporcional a la densidad, en lugar de exactamente igual a ella, hace que el algoritmo de 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 de 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 de cerca a la distribución deseada. Estos valores de muestra se producen iterativamente, ya que la distribución de la siguiente muestra depende únicamente del valor de muestra actual (por lo que la secuencia de muestras se convierte en una cadena 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. Entonces, con cierta probabilidad, el candidato es aceptado (en cuyo caso el valor candidato se utiliza en la próxima iteración) o rechazado (en cuyo caso se descarta el valor candidato, 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.

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 Metropolis (distribución de propuesta simétrica)

Vamos ${displaystyle f(x)}$ ser una función proporcional a la función de densidad de probabilidad deseada ${displaystyle P(x)}$ (a.k.a. a target distribution).

1. Iniciación: Elija un punto arbitrario ${displaystyle x_{t}$ ser la primera observación en la muestra y elegir una densidad de probabilidad arbitraria ${displaystyle g(xmid y)}$ (a veces escritos ${displaystyle Q(xmid 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, ${displaystyle g}$ se supone que es simétrico; en otras palabras, debe satisfacer ${displaystyle g(xmid y)=g(ymid x)}$. Una opción habitual es dejar ${displaystyle g(xmid y)}$ ser una distribución Gausiana centrada en ${displaystyle y}$, por lo que apunta más cerca ${displaystyle y}$ es más probable que se visite a continuación, haciendo la secuencia de muestras en un paseo aleatorio. La función ${displaystyle g}$ se conoce como el propuesta densidad o distribución de saltos.
2. Para cada iteración t:
   • Generar un candidato ${displaystyle x'}$ para la siguiente muestra seleccionando de la distribución ${displaystyle g(x'mid x_{t}}$.
   • Cálculo el ratio de aceptación ${displaystyle alpha =f(x')/f(x_{t}}$, que se utilizará para decidir si aceptar o rechazar al candidato. Porque... f es proporcional a la densidad de P, tenemos eso ${displaystyle alpha =f(x')/f(x_{t})=P(x')/P(x_{t}}$.
   • Aceptar o rechazar:
     • Generar un número aleatorio uniforme ${displaystyle uin [0,1]}$.
     • Si ${displaystyle uleq alpha }$, entonces aceptar el candidato estableciendo ${displaystyle x_{t+1}=x}$,
     • Si ${displaystyle u confianzaalpha }$, entonces Rechazo el candidato y establecer ${displaystyle x_{t+1}=x_{t}$ en su lugar.

Este algoritmo procede al intentar moverse aleatoriamente sobre el espacio de la muestra, a veces aceptando los movimientos y a veces permaneciendo en su lugar. Note that the acceptance ratio ${displaystyle alpha }$ indica cuán probable es la nueva muestra propuesta con respecto a la muestra actual, según la distribución cuya densidad es ${displaystyle P(x)}$. Si intentamos pasar a un punto más probable que el punto existente (es decir, un punto en una región de densidad superior) ${displaystyle P(x)}$ correspondiente a un ${displaystyle alpha >1geq u}$Siempre aceptaremos el movimiento. Sin embargo, si intentamos pasar a un punto menos probable, a veces rechazaremos el movimiento, y cuanto más grande sea la caída relativa en probabilidad, más probable es que rechacemos el nuevo punto. Así, tendremos a permanecer en (y devolver un gran número de muestras de) regiones de alta densidad de ${displaystyle P(x)}$, mientras que sólo ocasionalmente visitan regiones de baja densidad. Intuitivamente, por eso este algoritmo funciona y devuelve muestras que siguen la distribución deseada con densidad ${displaystyle P(x)}$.

En comparación con un algoritmo como el muestreo de rechazo adaptativo 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 sí siguen correctamente ${displaystyle P(x)}$, un conjunto de muestras cercanas se correlacionarán entre sí y no reflejarán correctamente la distribución. Esto significa que los tamaños de muestras eficaces pueden ser significativamente menores que el número de muestras realmente tomadas, lo que conduce a errores grandes.
• Aunque la cadena 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. As a result, a quemador el período suele ser necesario, donde se arroja un número inicial de muestras.

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 hasta tal punto y, por lo tanto, a menudo son las únicas soluciones disponibles cuando el número de dimensiones de la distribución que se va 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 que se utilizan hoy en día en muchas disciplinas.

En las 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 elevado, puede resultar difícil encontrar la distribución de saltos adecuada, 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 evitar 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 de muestrear de un espacio potencialmente de alta dimensión se reducirá a una colección de problemas para muestrear de pequeña dimensionalidad. 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 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. Las variables individuales luego se muestrean una a la vez, 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, el algoritmo de muestreo de Metropolis de rechazo adaptativo, un paso Metropolis-Hastings unidimensional simple o muestreo por sectores..

Derivación formal

El propósito del algoritmo de Metropolis-Hastings es generar una colección de estados según una distribución deseada ${displaystyle P(x)}$. Para lograr esto, el algoritmo utiliza un proceso de Markov, que asintomáticamente alcanza una distribución estacionaria única ${displaystyle pi (x)}$ tales que ${displaystyle pi (x)=P(x)}$.

Un proceso de Markov está definido por sus probabilidades de transición ${displaystyle P(x'mid x)}$, la probabilidad de transición de cualquier estado dado ${displaystyle x}$ a cualquier otro estado dado ${displaystyle x'}$. Tiene una distribución estacionaria única ${displaystyle pi (x)}$ cuando se cumplen las dos condiciones siguientes:

1. Existencia de distribución estacionaria: debe existir una distribución estacionaria ${displaystyle pi (x)}$. Una condición suficiente pero no necesaria es un equilibrio detallado, que requiere que cada transición ${displaystyle xto x}$ es reversible: para cada par de estados ${displaystyle x,x}$, la probabilidad de estar en estado ${displaystyle x}$ y transición al Estado ${displaystyle x'}$ debe ser igual a la probabilidad de estar en estado ${displaystyle x'}$ y transición al Estado ${displaystyle x}$, ${displaystyle pi (x)P(x'mid x)=pi (x')P(xmid x')}$.
2. Unicidad de la distribución estacionaria: la distribución estacionaria ${displaystyle pi (x)}$ debe 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 recidivante 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 (construyendo probabilidades de transición) que cumpla las dos condiciones anteriores, de tal manera que su distribución estacionaria ${displaystyle pi (x)}$ es elegido para ser ${displaystyle P(x)}$. La derivación del algoritmo comienza con la condición de balance detallado:

${displaystyle P(x'mid x)P(x)=P(xmid x')P(x'),}$

que se reescribe como

${displaystyle {frac {P(x'mid x)}{P(xmid x')}={frac {P(x')}{P(x)}}}}$

El enfoque es separar la transición en dos sub-pasos; la propuesta y el rechazo de la aceptación. Distribución de la propuesta ${displaystyle g(x'mid x)}$ es la probabilidad condicional de proponer un estado ${displaystyle x'}$ dado ${displaystyle x}$, y la distribución de aceptación ${displaystyle A(x',x)}$ es la probabilidad de aceptar el estado propuesto ${displaystyle x'}$. La probabilidad de transición se puede escribir como producto de ellos:

${displaystyle P(x'mid x)=g(x'mid x)A(x',x). }$

Insertando esta relación en la ecuación anterior, tenemos

${displaystyle {frac {A(x',x)}{A(x,x')}={frac {P(x'}{P(x)}}}{frac {g(xmid x')}{g(x'mid x)}}}}}}}$

El siguiente paso en la derivación es elegir un índice de aceptación que cumpla la condición anterior. Una elección común es la elección de Metropolis:

${displaystyle A(x',x)=min left(1,{frac {P(x'}{P(x)}}{frac {g(xmid x')}{g(x'mid x)}right). }$

Para esta relación de aceptación de Metropolis ${displaystyle A}$, o ${displaystyle A(x',x)=1}$ o ${displaystyle A(x,x')=1}$ y, de cualquier manera, la condición está satisfecha.

El algoritmo Metropolis-Hastings se puede escribir de la siguiente manera:

1. Inicial
   1. Escoge un estado inicial ${displaystyle x_{0}$.
   2. Set ${displaystyle t=0}$.
2. Iterate
   1. Generar un estado candidato aleatorio ${displaystyle x'}$ según ${displaystyle g(x'mid x_{t}}$.
   2. Cálculo la probabilidad de aceptación ${displaystyle A(x',x_{t})=min left(1,{frac {P(x')}{P(x_{t}}}}{frac {g(x_{t}mid x')}{g(x'mid x_{t}}}}right)}}}}}right)}$.
   3. Aceptar o rechazar:
      1. generar un número aleatorio uniforme ${displaystyle uin [0,1]}$;
      2. si ${displaystyle uleq A(x',x_{t}}$, entonces aceptar el nuevo estado y establecer ${displaystyle x_{t+1}=x}$;
      3. si ${displaystyle u confianzaA(x',x_{t}}$, entonces Rechazo el nuevo estado, y copia el viejo estado adelante ${displaystyle x_{t+1}=x_{t}$.
   4. Incremento: set ${displaystyle t=t+1}$.

Siempre que se cumplan las condiciones específicas, la distribución empírica de estados salvados ${displaystyle x_{0},ldotsx_{T}$ enfoque ${displaystyle P(x)}$. El número de iteraciones (${displaystyle T}$) requerido para calcular eficazmente ${displaystyle P(x)}$ depende del número de factores, incluyendo la relación entre ${displaystyle P(x)}$ y la distribución de propuestas y la exactitud deseada de la estimación. Para distribución en espacios estatales discretos, tiene que ser del orden del tiempo de autocorrelación del proceso de Markov.

Es importante notar que no está claro, en un problema general, qué distribución ${displaystyle g(x'mid x)}$ uno debe utilizar o el número de iteraciones necesarias para una estimación adecuada; ambos son parámetros libres del método, que debe ajustarse al problema particular en la mano.

Uso en integración numérica

Un uso común del algoritmo de Metropolis-Hastings es calcular una integral. Específicamente, considere un espacio ${displaystyle Omega subset mathbb {R}$ y una distribución de probabilidad ${displaystyle P(x)}$ sobre ${displaystyle Omega }$, ${displaystyle xin Omega }$. Metropolis–Hastings puede estimar una parte integral de la forma de

${displaystyle P(E)=int _{Omega }A(x)P(x),dx,}$

Donde ${displaystyle A(x)}$ es una función de interés (mesurable).

Por ejemplo, considere una estadística ${displaystyle E(x)}$ y su distribución de probabilidad ${displaystyle P(E)}$, que es una distribución marginal. Supongamos que el objetivo es estimar ${displaystyle P(E)}$ para ${displaystyle E}$ en la cola ${displaystyle P(E)}$. Formalmente, ${displaystyle P(E)}$ puede ser escrito como

${displaystyle P(E)=int _{Omega }P(Emid x)P(x),dx=int _{\ Omega.$

y, por consiguiente, estimación ${displaystyle P(E)}$ se puede lograr estimando el valor esperado de la función indicadora ${displaystyle A_{E}(x)equiv mathbf {1} _{E}(x)}$, que es 1 cuando ${displaystyle E(x)in [E,E+Delta E]}$ y cero de lo contrario. Porque... ${displaystyle E}$ está en la cola ${displaystyle P(E)}$, la probabilidad de dibujar un estado ${displaystyle x}$ con ${displaystyle E(x)}$ en la cola ${displaystyle P(E)}$ es proporcional a ${displaystyle P(E)}$, que es pequeña por definición. El algoritmo de Metrópolis-Hastings se puede utilizar aquí para mostrar (rare) estados más probables y así aumentar el número de muestras utilizadas para estimar ${displaystyle P(E)}$ en la cola. Esto se puede hacer por ejemplo utilizando una distribución de muestreo ${displaystyle pi (x)}$ para favorecer a esos estados (por ejemplo, ${displaystyle pi (x)propto e^{aE}$ con ${displaystyle a confía0}$).

Instrucciones paso a paso

Supongamos que el valor más reciente muestra ${displaystyle x_{t}$. Para seguir el algoritmo de Metropolis-Hastings, a continuación dibujamos un nuevo estado de propuesta ${displaystyle x'}$ con densidad de probabilidad ${displaystyle g(x'mid x_{t}}$ y calcula un valor

${displaystyle a=a_{1}a_{2},}$

dónde

${displaystyle - ¿Qué?$

es la relación de probabilidad (por ejemplo, Bayesian posterior) entre la muestra propuesta ${displaystyle x'}$ y la muestra anterior ${displaystyle x_{t}$, y

${displaystyle a_{2}={frac {g(x_{t}}{g(x'mid x_{t}}}}}$

es la proporción de la densidad de la propuesta en dos direcciones (desde ${displaystyle x_{t}$ a ${displaystyle x'}$ y a la inversa). Esto es igual a 1 si la densidad de la propuesta es simétrica. Entonces el nuevo estado ${displaystyle x_{t+1}$ es elegido según las reglas siguientes.

Si ${displaystyle ageq 1{:}$

${displaystyle x_{t+1}=x',}$

más:

${displaystyle x_{t+1}={begin{cases}x' afectados{text{with with probability }a,x_{t} {text{with probability {fnMicrosoft Sans Serif}$

La cadena Markov se inicia desde un valor inicial arbitrario ${displaystyle x_{0}$, y el algoritmo se ejecuta para muchas iteraciones hasta que este estado inicial es "olvidado". Estas muestras, que se descartan, se conocen como quemador. El conjunto restante de valores aceptados ${displaystyle x}$ representa una muestra de la distribución ${displaystyle P(x)}$.

El algoritmo funciona mejor si la densidad de la propuesta coincide con la forma de la distribución del objetivo ${displaystyle P(x)}$, desde el cual el muestreo directo es difícil, es decir ${displaystyle g(x'mid x_{t})approx P(x')}$. Si una densidad de propuesta Gausiana ${displaystyle g}$ se utiliza, el parámetro de variación ${displaystyle sigma ^{2}$ tiene que ser sintonizado durante el período de encendido. Esto se hace generalmente calculando Tasa de aceptación, que es la fracción de muestras propuestas que se acepta en una ventana de la última ${displaystyle N}$ muestras. La tasa de aceptación deseada depende de la distribución de destino, sin embargo se ha demostrado teóricamente que la tasa de aceptación ideal para una distribución gaussiiana unidimensional es de alrededor del 50%, disminuyendo a alrededor del 23% para un ${displaystyle N}$- distribución de objetivos gaussianos. Estas pautas pueden funcionar bien cuando se muestren de los posteriores Bayesianos suficientemente regulares, ya que a menudo siguen una distribución normal multivariable como se puede establecer utilizando el teorema Bernstein-von Mises.

Si ${displaystyle sigma ^{2}$ es demasiado pequeña, la cadena será mezcla lentamente (es decir, la tasa de aceptación será alta, pero muestras sucesivas se moverán alrededor del espacio lentamente, y la cadena convergerá sólo lentamente a ${displaystyle P(x)}$). Por otro lado, si ${displaystyle sigma ^{2}$ es demasiado grande, la tasa de aceptación será muy baja porque las propuestas probablemente aterrizarán en regiones de densidad de probabilidad mucho menor, por lo que ${displaystyle A_{1}$ será muy pequeño, y de nuevo la cadena convergerá muy lentamente. Una típicamente sintoniza la distribución de la propuesta para que los algoritmos acepten en el orden del 30% de todas las muestras – de acuerdo con las estimaciones teóricas mencionadas en el párrafo anterior.

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