¿Cómo funciona Hamiltonian Monte Carlo?
Hice el siguiente gráfico para explicar cómo entiendo actualmente el algoritmo HMC. Me gustaría que un experto en la materia verificara si esta comprensión es correcta o no. El texto de la siguiente diapositiva se copia a continuación para facilitar el acceso:
Montecarlo hamiltoniano: un satélite orbita un planeta. Cuanto más cerca esté el satélite del planeta, mayores serán los efectos de la gravedad. Esto significa, (A) mayor energía potencial y (B) mayor energía cinética necesaria para mantener la órbita. Esa misma energía cinética a una distancia mayor del planeta, expulsaría al satélite de la órbita. El satélite tiene la tarea de recopilar fotos de una región geográfica específica. Cuanto más cerca esté el satélite en órbita alrededor del planeta, más rápido se mueve en órbita, más veces pasa sobre la región, más fotografías recopila. Por el contrario, cuanto más lejos está un satélite del planeta, más lento se mueve en órbita, menos veces pasa sobre la región, menos fotografías recopila. En el contexto del muestreo, la distancia desde el planeta representa la distancia desde la expectativa de distribución. Un área de baja probabilidad está lejos de las expectativas; cuando “orbita esta probabilidad”, menor energía cinética significa menos muestras recolectadas durante un intervalo de tiempo fijo, mientras que cuando orbita una probabilidad más alta significa más muestras recolectadas dado el mismo intervalo de tiempo fijo. En una órbita dada, la energía total, cinética y potencial, es constante; sin embargo, la relación entre los dos no es sencilla. Las ecuaciones hamiltonianas relacionan los cambios de una a otra. Es decir, el gradiente de posición con respecto al tiempo es igual al impulso. Y el gradiente de impulso con respecto al tiempo es igual al gradiente de energía potencial con respecto a la posición. Para calcular qué tan lejos habrá viajado un satélite a lo largo de su trayectoria orbital, se debe utilizar la integración de salto, actualizando iterativamente los vectores de momento y posición. En el contexto del muestreo, la probabilidad es análoga a la distancia al planeta y el gradiente de energía potencial con respecto a la posición es el gradiente de la función de densidad de probabilidad con respecto a su parámetro de entrada, x. Esta información permite explorar la trayectoria orbital alrededor de varias entradas, X, correspondientes a la misma probabilidad, y.
Sin embargo, no estamos simplemente interesados en explorar una probabilidad, debemos explorar múltiples rutas orbitales. Para lograr esto, el impulso debe aumentarse aleatoriamente, acercando o alejando al satélite del planeta. Estas "patadas de impulso" aleatorias permiten orbitar diferentes probabilidades. Afortunadamente, las ecuaciones hamiltonianas garantizan que, independientemente de la probabilidad, el número de muestras recolectadas sea proporcional a la probabilidad, por lo que las muestras recolectadas siguen la forma de la distribución objetivo.
Mi pregunta es: ¿Es esta una forma precisa de pensar en cómo funciona el Hamiltoniano Montecarlo?
Editar:
Lo implementé en un código basado en mi comprensión del algoritmo. Funciona para un gaussiano con mu = 0, sigma = 1. Pero si cambio sigma se rompe. Se agradecería cualquier información.
import numpy as np
import random
import scipy.stats as st
import matplotlib.pyplot as plt
from autograd import grad
def normal(x,mu,sigma):
numerator = np.exp((-(x-mu)**2)/(2*sigma**2))
denominator = sigma * np.sqrt(2*np.pi)
return numerator/denominator
def neg_log_prob(x,mu,sigma):
num = np.exp(-1*((x-mu)**2)/2*sigma**2)
den = sigma*np.sqrt(np.pi*2)
return -1*np.log(num/den)
def HMC(mu=0.0,sigma=1.0,path_len=1,step_size=0.25,initial_position=0.0,epochs=1_000):
# setup
steps = int(path_len/step_size) -1 # path_len and step_size are tricky parameters to tune...
samples = [initial_position]
momentum_dist = st.norm(0, 1)
# generate samples
for e in range(epochs):
q0 = np.copy(samples[-1])
q1 = np.copy(q0)
p0 = momentum_dist.rvs()
p1 = np.copy(p0)
dVdQ = -1*(q0-mu)/(sigma**2) # gradient of PDF wrt position (q0) aka momentum wrt position
# leapfrog integration begin
for s in range(steps):
p1 += step_size*dVdQ/2 # as potential energy increases, kinetic energy decreases
q1 += step_size*p1 # position increases as function of momentum
p1 += step_size*dVdQ/2 # second half "leapfrog" update to momentum
# leapfrog integration end
p1 = -1*p1 #flip momentum for reversibility
#metropolis acceptance
q0_nlp = neg_log_prob(x=q0,mu=mu,sigma=sigma)
q1_nlp = neg_log_prob(x=q1,mu=mu,sigma=sigma)
p0_nlp = neg_log_prob(x=p0,mu=0,sigma=1)
p1_nlp = neg_log_prob(x=p1,mu=0,sigma=1)
# Account for negatives AND log(probabiltiies)...
target = q0_nlp - q1_nlp # P(q1)/P(q0)
adjustment = p1_nlp - p0_nlp # P(p1)/P(p0)
acceptance = target + adjustment
event = np.log(random.uniform(0,1))
if event <= acceptance:
samples.append(q1)
else:
samples.append(q0)
return samples
Ahora funciona aquí:
mu, sigma = 0,1
trial = HMC(mu=mu,sigma=sigma,path_len=2,step_size=0.25)
# What the dist should looks like
lines = np.linspace(-6,6,10_000)
normal_curve = [normal(x=l,mu=mu,sigma=sigma) for l in lines]
# Visualize
plt.plot(lines,normal_curve)
plt.hist(trial,density=True,bins=20)
plt.show()
Pero se rompe cuando cambio sigma a 2.
# Generate samples
mu, sigma = 0,2
trial = HMC(mu=mu,sigma=sigma,path_len=2,step_size=0.25)
# What the dist should looks like
lines = np.linspace(-6,6,10_000)
normal_curve = [normal(x=l,mu=mu,sigma=sigma) for l in lines]
# Visualize
plt.plot(lines,normal_curve)
plt.hist(trial,density=True,bins=20)
plt.show()
¿Algunas ideas? Siento que estoy cerca de "conseguirlo".
Respuestas
Antes de responder a la pregunta sobre una forma intuitiva de pensar sobre el Hamiltoniano Montecarlo, probablemente sea mejor tener una comprensión realmente firme del MCMC regular. Dejemos de lado la metáfora del satélite por ahora.
MCMC es útil cuando desea una muestra sin sesgos de una distribución en la que solo tiene algo disponible que es proporcional al PDF, pero no al PDF en sí. Esto surge en (por ejemplo) simulaciones de física: el PDF viene dado por la distribución de Boltzmann, p ~ exp (-E / kT), pero lo que se puede calcular para cualquier configuración del sistema es E, no p. La constante de proporcionalidad no se conoce, porque la integral de exp (-E / kT) sobre todo el espacio de configuración posible suele ser demasiado difícil de calcular. MCMC resuelve ese problema haciendo una caminata aleatoria de una manera específica, donde la probabilidad de dar ("aceptar") cada paso está relacionada con la razón de los valores p (la constante de proporcionalidad se cancela). Con el tiempo, la distribución de muestras aceptadas del paseo aleatorio converge al PDF que queremos, sin necesidad de calcular explícitamente p.
Tenga en cuenta que en lo anterior, cualquier método para tomar pasos aleatorios es igualmente válido, siempre que el caminante aleatorio pueda explorar todo el espacio. El criterio de aceptación garantiza que las muestras seleccionadas converjan al PDF real. En la práctica, se utiliza una distribución gaussiana alrededor de la muestra actual (y el sigma se puede variar para que la fracción de pasos aceptados se mantenga relativamente alta). En principio, no habría nada de malo en tomar pasos de cualquier otra distribución continua ("distribución de salto") alrededor de la muestra actual, aunque la convergencia puede ser mucho más lenta.
Ahora, el Hamiltoniano Montecarlo extiende la metáfora de la física al tratar específicamente de dar pasos en una dirección que es más probable que sea aceptada que un paso gaussiano. Los pasos son los que tomaría un integrador de salto, si estuviera tratando de resolver el movimiento de un sistema donde la energía potencial fuera E. Estas ecuaciones de movimiento también incluyen un término de energía cinética, con una "masa" (no literalmente física) y "impulso". Los pasos que toma el integrador de salto en "tiempo" se pasan como propuestas al algoritmo MCMC.
¿Por qué funciona esto? El MC gaussiano da pasos de la misma distancia en todas direcciones con igual probabilidad; lo único que lo predispone hacia áreas más densamente pobladas del PDF es que es más probable que se rechacen los pasos en la dirección incorrecta. El Hamiltoniano MC propone pasos tanto en la dirección del gradiente E como en la dirección del movimiento acumulado en los pasos recientes (dirección y magnitud del "momento"). Esto permite una exploración más rápida del espacio y también una mayor probabilidad de llegar más rápido a regiones más densamente pobladas.
Ahora, la metáfora del satélite: creo que esta no es una forma muy útil de pensar en ello. Los satélites se mueven en una órbita exacta; lo que tienes aquí es bastante aleatorio, más como una partícula de gas en un recipiente con otras partículas. Cada colisión aleatoria le da un "paso"; con el tiempo, la partícula estará en todas partes del contenedor con la misma probabilidad (ya que la PDF aquí es igual en todas partes, excepto las paredes que representan una energía muy alta / efectivamente cero PDF). La MCMC gaussiana es como una partícula de masa cero que realiza una caminata aleatoria (o una partícula de masa distinta de cero en un medio relativamente viscoso): llegará allí mediante un movimiento browniano, pero no necesariamente rápido. El Hamiltoniano MC es una partícula con una masa distinta de cero: puede acumular suficiente impulso para seguir en la misma dirección a pesar de las colisiones, por lo que a veces puede dispararse de un extremo del contenedor a otro (dependiendo de su masa frente a la frecuencia / magnitud de las colisiones). Todavía rebotaría en las paredes, por supuesto, pero en general tendería a explorar más rápido.