我制作了下图来解释我目前是如何理解 HMC 算法的。如果这种理解正确或不正确，我希望得到主题专家的验证。下面幻灯片中的文本复制如下，以便于访问：
哈密顿蒙特卡洛：卫星绕行星运行。卫星离地球越近，重力的影响就越大。这意味着，（A）更高的势能和（B）维持轨道所需的更高动能。在离地球更远的地方，同样的动能会将卫星从轨道上弹出。该卫星的任务是收集特定地理区域的照片。卫星绕地球运行的距离越近，它在轨道上的移动速度就越快，它越过该区域的次数越多，它收集的照片就越多。相反，卫星离地球越远，它在轨道上移动的速度越慢，它越过该区域的次数越少，它收集的照片就越少。在采样的上下文中，与行星的距离表示与分布期望的距离。一个低可能性的区域远非预期；当“环绕这种可能性”时，较低的动能意味着在固定时间间隔内收集的样本较少，而当环绕较高的可能性时，意味着在相同的固定时间间隔内收集的样本更多。在给定的轨道上，总能量、动能和势能是恒定的；然而，两者的关系并不简单。哈密顿方程将一个变化与另一个变化联系起来。也就是说，位置相对于时间的梯度等于动量。并且动量相对于时间的梯度等于势能相对于位置的梯度。为了计算卫星将沿着其轨道路径行进多远，必须使用越级积分，迭代更新动量和位置向量。在抽样的背景下，可能性类似于与行星的距离，势能相对于位置的梯度是概率密度函数相对于其输入参数 x 的梯度。该信息允许探索对应于相同可能性 y 的各种输入 X 的轨道路径。
然而，我们不仅仅对探索一种可能性感兴趣，我们还必须探索多条轨道路径。要做到这一点，必须随机增加动量，使卫星离地球更近或更远。这些随机的“动量踢”允许不同的可能性被环绕。幸运的是，汉密尔顿方程确保无论可能性如何，收集的样本数量与可能性成正比，因此收集的样本遵循目标分布的形状。
我的问题是 - 这是思考哈密顿蒙特卡洛如何工作的准确方法吗？
编辑：
根据我对算法的理解，我已经在一些代码中实现了。它适用于 mu=0，sigma=1 的高斯。但是，如果我更改 sigma，它就会中断。任何见解将不胜感激。
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
现在它在这里工作：
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()
但是当我将 sigma 更改为 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()
有任何想法吗？我觉得我已经接近“得到它”了。