Montecarlo hamiltoniano

Algoritmo de muestreo

Muestreo de Monte Carlo hamiltoniano de una distribución de probabilidad bidimensional

El algoritmo Monte Carlo hamiltoniano (originalmente conocido como Monte Carlo híbrido ) es un método Monte Carlo de cadena de Markov para obtener una secuencia de muestras aleatorias cuya distribución converge a una distribución de probabilidad objetivo que es difícil de muestrear directamente. Esta secuencia se puede utilizar para estimar integrales de la distribución objetivo, como valores esperados y momentos .

El Hamiltonian Monte Carlo corresponde a una instancia del algoritmo Metropolis-Hastings , con una evolución de dinámica hamiltoniana simulada utilizando un integrador numérico reversible en el tiempo y que preserva el volumen (normalmente el integrador leapfrog ) para proponer un movimiento a un nuevo punto en el espacio de estados. En comparación con el uso de una distribución de propuesta de paseo aleatorio gaussiano en el algoritmo Metropolis-Hastings, el Hamiltonian Monte Carlo reduce la correlación entre estados muestreados sucesivos al proponer movimientos a estados distantes que mantienen una alta probabilidad de aceptación debido a las propiedades aproximadas de conservación de energía de la dinámica hamiltoniana simulada cuando se utiliza un integrador simpléctico . La correlación reducida significa que se necesitan menos muestras de cadena de Markov para aproximar integrales con respecto a la distribución de probabilidad objetivo para un error de Monte Carlo dado .

El algoritmo fue propuesto originalmente por Simon Duane, Anthony Kennedy, Brian Pendleton y Duncan Roweth en 1987 para cálculos en cromodinámica cuántica de red . ^[1] En 1996, Radford M. Neal mostró cómo el método podría usarse para una clase más amplia de problemas estadísticos, en particular redes neuronales artificiales . ^[2] Sin embargo, la carga de tener que proporcionar gradientes de la red bayesiana retrasó la adopción más amplia del algoritmo en estadística y otras disciplinas cuantitativas, hasta que a mediados de la década de 2010 los desarrolladores de Stan implementaron HMC en combinación con la diferenciación automática . ^[3]

Algoritmo

Supongamos que la distribución objetivo a muestrear es para ( ) y se requiere una cadena de muestras . ${\displaystyle f(\mathbf {x} )}$ ${\displaystyle \mathbf {x} \in \mathbb {R} ^{d}}$ ${\displaystyle d\geq 1}$ ${\displaystyle \mathbf {X} _{0},\mathbf {X} _{1},\mathbf {X} _{2},\ldots }$

Las ecuaciones de Hamilton son

${\displaystyle {\frac {{\text{d}}x_{i}}{{\text{d}}t}}={\frac {\partial H}{\partial p_{i}}}\quad {\text{and}}\quad {\dfrac {{\text{d}}p_{i}}{{\text{d}}t}}=-{\dfrac {\partial H}{\partial x_{i}}}}$

donde y son el componente n del vector de posición y momento respectivamente y es el hamiltoniano. Sea una matriz de masas que es simétrica y definida positiva, entonces el hamiltoniano es ${\displaystyle x_{i}}$ ${\displaystyle p_{i}}$ ${\displaystyle i}$ ${\displaystyle H}$ ${\displaystyle M}$

${\displaystyle H(\mathbf {x} ,\mathbf {p} )=U(\mathbf {x} )+{\dfrac {1}{2}}\mathbf {p} ^{\text{T}}M^{-1}\mathbf {p} }$

¿Dónde está la energía potencial ? La energía potencial de un objetivo se expresa como ${\displaystyle U(\mathbf {x} )}$

${\displaystyle U(\mathbf {x} )=-\ln f(\mathbf {x} )}$

que proviene del factor de Boltzmann . Nótese que el hamiltoniano es adimensional en esta formulación porque el peso de probabilidad exponencial tiene que estar bien definido. Por ejemplo, en simulaciones a temperatura finita el factor (con la constante de Boltzmann ) se absorbe directamente en y . ${\displaystyle H}$ ${\displaystyle \exp \left(-H\right)}$ ${\displaystyle T}$ ${\displaystyle k_{\text{B}}T}$ ${\displaystyle k_{\text{B}}}$ ${\displaystyle U}$ ${\displaystyle M}$

El algoritmo requiere un entero positivo para el número de pasos de salto y un número positivo para el tamaño del paso . Supongamos que la cadena está en . Sea . Primero, se extrae un momento gaussiano aleatorio de . A continuación, la partícula se moverá bajo dinámica hamiltoniana durante un tiempo , esto se hace resolviendo numéricamente las ecuaciones de Hamilton utilizando el algoritmo de salto . Los vectores de posición y momento después del tiempo utilizando el algoritmo de salto son: ^[4] ${\displaystyle L}$ ${\displaystyle \Delta t}$ ${\displaystyle \mathbf {X} _{n}=\mathbf {x} _{n}}$ ${\displaystyle \mathbf {x} _{n}(0)=\mathbf {x} _{n}}$ ${\displaystyle \mathbf {p} _{n}(0)}$ ${\displaystyle {\text{N}}\left(\mathbf {0} ,M\right)}$ ${\displaystyle L\Delta t}$ ${\displaystyle \Delta t}$

${\displaystyle \mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)=\mathbf {p} _{n}(t)-{\dfrac {\Delta t}{2}}\nabla \left.U(\mathbf {x} )\right|_{\mathbf {x} =\mathbf {x} _{n}(t)}}$

${\displaystyle \mathbf {x} _{n}(t+\Delta t)=\mathbf {x} _{n}(t)+\Delta t\mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)}$

${\displaystyle \mathbf {p} _{n}(t+\Delta t)=\mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)-{\dfrac {\Delta t}{2}}\nabla \left.U(\mathbf {x} )\right|_{\mathbf {x} =\mathbf {x} _{n}(t+\Delta t)}}$

Estas ecuaciones se deben aplicar a y veces para obtener y . ${\displaystyle \mathbf {x} _{n}(0)}$ ${\displaystyle \mathbf {p} _{n}(0)}$ ${\displaystyle L}$ ${\displaystyle \mathbf {x} _{n}(L\Delta t)}$ ${\displaystyle \mathbf {p} _{n}(L\Delta t)}$

El algoritmo de salto de rana es una solución aproximada al movimiento de partículas clásicas que no interactúan. Si es exacta, la solución nunca cambiará la distribución de energía inicial generada aleatoriamente, ya que la energía se conserva para cada partícula en presencia de un campo de energía potencial clásico. Para alcanzar una distribución de equilibrio termodinámico, las partículas deben tener algún tipo de interacción con, por ejemplo, un baño de calor circundante, de modo que todo el sistema pueda asumir diferentes energías con probabilidades de acuerdo con la distribución de Boltzmann.

Una forma de mover el sistema hacia una distribución de equilibrio termodinámico es cambiar el estado de las partículas utilizando el algoritmo de Metropolis-Hastings . Primero se aplica el paso de salto y luego un paso de Metropolis-Hastings.

La transición de a es ${\displaystyle \mathbf {X} _{n}=\mathbf {x} _{n}}$ ${\displaystyle \mathbf {X} _{n+1}}$

${\displaystyle \mathbf {X} _{n+1}|\mathbf {X} _{n}=\mathbf {x} _{n}={\begin{cases}\mathbf {x} _{n}(L\Delta t)&{\text{with probability }}\alpha \left(\mathbf {x} _{n}(0),\mathbf {x} _{n}(L\Delta t)\right)\\\mathbf {x} _{n}(0)&{\text{otherwise}}\end{cases}}}$

dónde

${\displaystyle \alpha \left(\mathbf {x} _{n}(0),\mathbf {x} _{n}(L\Delta t)\right)={\text{min}}\left(1,{\dfrac {\exp \left[-H(\mathbf {x} _{n}(L\Delta t),\mathbf {p} _{n}(L\Delta t))\right]}{\exp \left[-H(\mathbf {x} _{n}(0),\mathbf {p} _{n}(0))\right]}}\right).}$

Una actualización completa consiste en primero muestrear aleatoriamente los momentos (independientemente de las iteraciones anteriores), luego integrar las ecuaciones de movimiento (por ejemplo, con leapfrog) y, finalmente, obtener la nueva configuración a partir del paso de aceptación/rechazo de Metropolis-Hastings. Este mecanismo de actualización se repite para obtener . ${\displaystyle \mathbf {p} }$ ${\displaystyle \mathbf {X} _{n+1},\mathbf {X} _{n+2},\mathbf {X} _{n+3},\ldots }$

Muestra sin vuelta en U

El muestreador sin giro en U (NUTS) ^[5] es una extensión que se controla automáticamente. El ajuste es fundamental. Por ejemplo, en el caso unidimensional , el potencial es que corresponde al potencial de un oscilador armónico simple . Si es demasiado grande, la partícula oscilará y, por lo tanto, se desperdiciará tiempo computacional. Si es demasiado pequeño, la partícula se comportará como un paseo aleatorio. ${\displaystyle L}$ ${\displaystyle L}$ ${\displaystyle {\text{N}}(0,1/{\sqrt {k}})}$ ${\displaystyle U(x)=kx^{2}/2}$ ${\displaystyle L}$ ${\displaystyle L}$

En términos generales, NUTS ejecuta la dinámica hamiltoniana hacia adelante y hacia atrás en el tiempo de manera aleatoria hasta que se satisface una condición de giro en U. Cuando eso sucede, se elige un punto aleatorio de la trayectoria para la muestra MCMC y el proceso se repite desde ese nuevo punto.

En detalle, se construye un árbol binario para rastrear la ruta de los pasos de salto de rana. Para producir una muestra MCMC, se lleva a cabo un procedimiento iterativo. Se muestrea una variable de corte. Sea y la posición y el momento de la partícula hacia adelante respectivamente. De manera similar, y para la partícula hacia atrás. En cada iteración, el árbol binario selecciona al azar de manera uniforme para mover la partícula hacia adelante hacia adelante en el tiempo o la partícula hacia atrás hacia atrás en el tiempo. Además, para cada iteración, el número de pasos de salto de rana aumenta en un factor de 2. Por ejemplo, en la primera iteración, la partícula hacia adelante se mueve hacia adelante en el tiempo usando 1 paso de salto de rana. En la siguiente iteración, la partícula hacia atrás se mueve hacia atrás en el tiempo usando 2 pasos de salto de rana. ${\displaystyle U_{n}\sim {\text{Uniform}}(0,\exp(-H[\mathbf {x} _{n}(0),\mathbf {p} _{n}(0)]))}$ ${\displaystyle \mathbf {x} _{n}^{+}}$ ${\displaystyle \mathbf {p} _{n}^{+}}$ ${\displaystyle \mathbf {x} _{n}^{-}}$ ${\displaystyle \mathbf {p} _{n}^{-}}$

El procedimiento iterativo continúa hasta que se cumple la condición de giro en U, es decir

${\displaystyle (\mathbf {x} _{n}^{+}-\mathbf {x} _{n}^{-})\cdot \mathbf {p} _{n}^{-}<0\quad {\text{or}}\quad .(\mathbf {x} _{n}^{+}-\mathbf {x} _{n}^{-})\cdot \mathbf {p} _{n}^{+}<0}$

o cuando el hamiltoniano se vuelve inexacto

${\displaystyle \exp \left[-H(\mathbf {x} _{n}^{+},\mathbf {p} _{n}^{+})+\delta \right]<U_{n}}$

o

${\displaystyle \exp \left[-H(\mathbf {x} _{n}^{-},\mathbf {p} _{n}^{-})+\delta \right]<U_{n}}$

donde, por ejemplo, . ${\displaystyle \delta =1000}$

Una vez que se cumple la condición de U-Turn, la siguiente muestra MCMC, , se obtiene muestreando uniformemente la trayectoria de salto de rana trazada por el árbol binario que satisface ${\displaystyle \mathbf {x} _{n+1}}$ ${\displaystyle \{\mathbf {x} _{n}^{-},\ldots ,\mathbf {x} _{n}(-\Delta t),\mathbf {x} _{n}(0),\mathbf {x} _{n}(\Delta t),\ldots ,\mathbf {x} _{n}^{+}\}}$

${\displaystyle U_{n}<\exp \left[-H(\mathbf {x_{n+1}} ,\mathbf {p_{n+1})} \right]}$

Generalmente esto se cumple si los parámetros HMC restantes son sensatos.

Véase también

• Método dinámico de Monte Carlo
• Software para modelado molecular de Monte Carlo
• Stan , un lenguaje de programación probabilístico que implementa HMC.
• PyMC , un lenguaje de programación probabilística que implementa HMC.
• Algoritmo de Langevin ajustado a la metrópolis

Referencias

1. Duane, Simon; Kennedy, Anthony D.; Pendleton, Brian J.; Roweth, Duncan (1987). "Monte Carlo híbrido". Physics Letters B . 195 (2): 216–222. Código Bibliográfico :1987PhLB..195..216D. doi :10.1016/0370-2693(87)91197-X.
2. Neal, Radford M. (1996). "Implementación de Monte Carlo". Aprendizaje bayesiano para redes neuronales . Apuntes de clase sobre estadística. Vol. 118. Springer. págs. 55–98. doi :10.1007/978-1-4612-0745-0_3. ISBN . 0-387-94724-8.
3. Gelman, Andrew; Lee, Daniel; Guo, Jiqiang (2015). "Stan: un lenguaje de programación probabilístico para la inferencia y optimización bayesiana". Revista de estadística educativa y conductual . 40 (5): 530–543. doi :10.3102/1076998615606113. S2CID  18351694.
4. Betancourt, Michael (15 de julio de 2018). "Una introducción conceptual al método Monte Carlo hamiltoniano". arXiv : 1701.02434 [stat.ME].
5. Hoffman, Matthew D; Gelman, Andrew (2014). "El muestreador sin giro en U: configuración adaptativa de longitudes de ruta en el modelo Monte Carlo hamiltoniano". Journal of Machine Learning Research . 15 (1): 1593–1623 . Consultado el 28 de marzo de 2024 .

Lectura adicional

• Betancourt, Michael; Girolami, Mark (2015). "Monte Carlo hamiltoniano para modelos jerárquicos". En Upadhyay, Satyanshu Kumar; et al. (eds.). Tendencias actuales en la metodología bayesiana con aplicaciones . CRC Press. págs. 79–101. ISBN 978-1-4822-3511-1.
• Betancourt, Michael (2018). "Una introducción conceptual al método Monte Carlo hamiltoniano". arXiv : 1701.02434 [stat.ME].
• Barbu, Adrian; Zhu, Song-Chun (2020). "Monte Carlo hamiltoniano y de Langevin". Métodos de Monte Carlo . Singapur: Springer. págs. 281–326. ISBN 978-981-13-2970-8.
• Neal, Radford M (2011). "MCMC usando dinámica hamiltoniana" (PDF) . En Steve Brooks; Andrew Gelman; Galin L. Jones; Xiao-Li Meng (eds.). Manual de Monte Carlo de cadenas de Markov . Chapman y Hall/CRC. ISBN. 9781420079418.

Enlaces externos

• Betancourt, Michael. "Inferencia bayesiana eficiente con el método Monte Carlo hamiltoniano". MLSS Islandia 2014 – vía YouTube .
• McElreath, Richard. "Monte Carlo con cadenas de Markov". Replanteamiento estadístico 2022 – vía YouTube .
• Hamiltoniano Monte Carlo desde cero
• Optimización y métodos de Monte Carlo