哈密 頓蒙特卡洛如何工作？
我製作了下圖來解釋我目前是如何理解 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()
有任何想法嗎？我覺得我已經接近“得到它”了。
在回答有關以直觀方式思考哈密頓蒙特卡洛的問題之前，最好真正牢牢掌握常規 MCMC。讓我們暫時擱置衛星比喻。
當您想要從分佈中獲得無偏樣本時， MCMC很有用到 PDF，但不是 PDF 本身。這齣現在（例如）物理模擬中：PDF 由玻爾茲曼分佈 p ~ exp(-E/kT) 給出，但是您可以為系統的任何配置計算的是 E，而不是 p。比例常數是未知的，因為 exp(-E/kT) 在可能配置的整個空間上的積分通常很難計算。MCMC 通過以特定方式進行隨機遊走來解決該問題，其中採取（“接受”）每一步的概率與 p 值的比率有關（比例常數抵消）。隨著時間的推移，隨機遊走中接受樣本的分佈收斂到我們想要的 PDF，而無需顯式計算 p。
請注意，在上面，任何採取隨機步驟的方法都同樣有效，只要隨機遊走者可以探索整個空間。驗收標准保證所選樣本收斂到真實 PDF。在實踐中，使用當前樣本周圍的高斯分佈（並且可以改變 sigma，以便接受步驟的比例保持相對較高）。從當前樣本周圍的任何其他連續分佈（“跳躍分佈”）採取步驟原則上沒有錯，儘管收斂可能要慢得多。
現在，哈密頓量蒙特卡洛擴展了物理隱喻，特別嘗試朝著比高斯步驟更容易被接受的方向邁出一步。如果它試圖求解勢能為 E 的系統的運動，則這些步驟是跳躍式積分器將採取的步驟。這些運動方程還包括動能項，具有（不是字面意義上的物理）“質量”和“勢頭”。然後，越級積分器在“時間”中採取的步驟作為建議傳遞給 MCMC 算法。
為什麼這行得通？高斯 MC 在每個方向上以相同的概率走相同的距離；唯一使它偏向於 PDF 中人口更密集的區域的是，朝著錯誤方向的步驟更有可能被拒絕。Hamiltonian MC 提出了 E 梯度方向的步驟，以及最近步驟中累積運動的方向（“動量”的方向和大小）。這可以更快地探索空間，也可以更快地到達人口稠密的地區。
現在，衛星比喻：我認為這不是一個非常有用的思考方式。衛星在精確的軌道上運行；你這裡的東西是非常隨機的，更像是容器中的一個氣體粒子和其他粒子。每一次隨機碰撞都會給你一個“台階”；隨著時間的推移，粒子將以相等的概率出現在容器中的任何地方（因為這裡的 PDF 在任何地方都是相等的，除了代表非常高能量/實際上為零的 PDF 的牆壁）。高斯 MCMC 就像一個有效的零質量粒子進行隨機遊走（或相對粘性介質中的非零質量粒子）：它將通過布朗運動到達那裡，但不一定很快。哈密頓量 MC 是一個非零質量的粒子：儘管發生碰撞，它可能會聚集足夠的動量以保持相同的方向運動，因此它有時可能會從容器的一端射到另一端（取決於它的質量與碰撞的頻率/幅度）。當然，它仍然會從牆壁上反彈，但通常會更快地探索。