Markov Chain Monte Carlo/zh

Article
Topic area	Bayesian Methods
Prerequisites	Bayesian Inference, Markov Chain

概述

馬爾可夫鏈蒙特卡羅 (MCMC) 是一類算法,用於從僅在歸一化常數意義下已知的概率分佈中抽取樣本。MCMC 不直接抽樣,而是構造一條馬爾可夫鏈,使其平穩分佈為目標分佈,運行該鏈足夠長的時間以遺忘其初始狀態,並將後續狀態視為近似樣本。然後通過對這些狀態的函數取平均,來估計目標分佈下的期望。 該技術是現代貝葉斯推斷的基礎:形如 $ p(\theta \mid \mathbf{y}) \propto p(\mathbf{y} \mid \theta) p(\theta) $ 的後驗分佈通常難以處理,但其未歸一化的概率密度計算成本很低,這正是 MCMC 所需要的。

MCMC 位於精確抽樣和變分推斷之間:精確抽樣需要可處理的歸一化常數,而變分推斷以固定的近似差距為代價優化一個可處理的替代分佈。只要計算資源充足,MCMC 估計量是漸近無偏的;其主要的實際問題是混合時間、自相關和調參。該框架在 1953 年由Metropolis 算法在統計學中被推廣,並於 1970 年由 Hastings 加以擴展,隨後在 1990 年代,通用採樣器如Gibbs 和 BUGS made使得後驗推斷成為常規做法,從而徹底改變了貝葉斯計算。

直覺

設想要在高維後驗分佈上計算某函數的平均值,其概率密度可以求值但不可逆。直接抽樣行不通,因為不知道歸一化常數;重要性抽樣行不通,因為在高維空間中不存在好的提議分佈。MCMC 通過偏向高密度區域的隨機遊走來迴避這兩個問題:每一步對當前狀態提出一個小的改動,然後接受或拒絕該提議,使得從長期來看,該鏈在每個區域上停留的時間與目標密度成正比。

關鍵事實是:只需要計算目標密度的比值。由於 $ p(\theta) = \tilde{p}(\theta) / Z $,其中 $ Z $ 是未知的歸一化常數,因此比值 $ \tilde{p}(\theta')/\tilde{p}(\theta) $ 等於 $ p(\theta')/p(\theta) $。未知的常數被約去,因此該鏈可以僅使用未歸一化的密度來模擬。

基礎

馬爾可夫鏈是一列隨機變量 $ \theta_0, \theta_1, \ldots $,其中 $ \theta_{t+1} $ 僅依賴於 $ \theta_t $,由一個轉移核 $ T(\theta' \mid \theta) $ 控制。若分佈 $ \pi $ 滿足以下條件,則稱其對 $ T $ 是平穩的:

$ {\displaystyle \pi(\theta') = \int T(\theta' \mid \theta) \pi(\theta)\, d\theta.} $

平穩性的一個充分(但非必要)條件是細緻平衡:

$ {\displaystyle \pi(\theta) T(\theta' \mid \theta) = \pi(\theta') T(\theta \mid \theta'),} $

它表明任意兩個狀態之間的概率流是對稱的。如果該鏈還滿足不可約性(每個狀態都可從任意其他狀態到達)且非周期,則遍歷定理保證:沿該鏈的時間平均收斂到 $ \pi $ 下的期望。 MCMC 算法是用於設計轉移核的方法,使其平穩分佈為所選擇的目標分佈。

Metropolis-Hastings

Metropolis-Hastings (MH) 算法是規範的 MCMC 方法。給定當前狀態 $ \theta $,從提議分佈 $ q(\theta' \mid \theta) $ 中抽取候選 $ \theta' $,並以以下概率接受:

$ {\displaystyle \alpha(\theta, \theta') = \min\!\left(1,\; \frac{\tilde{p}(\theta')\, q(\theta \mid \theta')}{\tilde{p}(\theta)\, q(\theta' \mid \theta)}\right).} $

若提議分佈對稱(即原始 Metropolis 情形),則 $ q $ 比值被約去,接受概率簡化為 $ \min(1, \tilde{p}(\theta')/\tilde{p}(\theta)) $。根據構造,所得的轉移核關於 $ p $ 滿足細緻平衡,因此 $ p $ 是平穩的。提議分佈是唯一的設計槓桿:過窄的提議產生高接受率但探索緩慢;過寬的提議則常被拒絕。對於高維空間中的隨機遊走提議,以約 0.234 的接受率為目標進行步長自適應是漸近最優的。

Gibbs 抽樣

Gibbs 抽樣利用多元目標分佈的結構,每次更新一個坐標(或一組坐標),條件依賴於其他坐標。對於 $ \theta = (\theta_1, \ldots, \theta_d) $,Gibbs 採樣器的一次掃描抽取:

$ {\displaystyle \theta_i^{(t+1)} \sim p\!\left(\theta_i \,\big|\, \theta_{-i}^{(t+1, t)}\right),} $

其中 $ \theta_{-i} $ 表示其他坐標的最新值。每次條件更新本身就是一個接受概率為一的 MH 步,因此當條件分佈可以解析抽樣時無需調參。共軛分層模型、混合模型以及具有可處理馬爾可夫毯的圖模型都是其自然適用對象。當某個條件分佈不可直接抽樣時,可在 Gibbs 內部用 MH 步或切片抽樣動作來替代,從而得到 Metropolis-within-Gibbs 方案。

哈密頓蒙特卡羅

隨機遊走 MH 在高維空間和強相關目標上混合很差。哈密頓蒙特卡羅 (HMC) 通過引入一個輔助動量變量 $ r $ 來擴展狀態,以勢能 $ U(\theta) = -\log \tilde{p}(\theta) $ 和動能 $ K(r) = \tfrac{1}{2} r^\top M^{-1} r $ 模擬哈密頓動力學,從而提議利用梯度信息的長軌跡。 蛙跳積分器保持體積不變且可逆;剩餘的離散化誤差由對聯合能量 $ (\theta, r) $ 進行的單次 MH 接受-拒絕步加以修正。由於軌跡沿梯度方向運動,HMC 抑制了隨機遊走行為,並且在維度上的擴展性更好。No-U-Turn 採樣器 (NUTS) 通過將軌跡延伸到開始折返為止,自動選擇軌跡長度,消除了 HMC 調參的一大難題,從而成為 Stan 和 PyMC 中的默認採樣器。

診斷與收斂

對於有限鏈,MCMC 估計是有偏的;隨着鏈遺忘其初始化,偏差會衰減。要丟棄一段預熱期,然後用自相關、有效樣本量 (ESS) 以及由多條從過度分散點出發的鏈計算的潛在尺度縮減因子 $ \hat{R} $ 來概括剩餘樣本。一種常用的經驗規則是至少運行四條鏈,並將 $ \hat{R} > 1.01 $ 視為未收斂的證據。 軌跡圖、秩圖以及 HMC 的發散診斷能夠捕捉標量診斷難以發現的病態情形,例如漏斗形和多模態。這些工具都無法證明收斂,只能指出失敗。

變體與聯繫

除了規範的採樣器之外,還有多種變體應對特定挑戰。可逆跳躍 MCMC 處理維度可變的模型。並行回火和副本交換在不同溫度下運行多條鏈並交換狀態以逃離模態。序貫蒙特卡羅將重要性抽樣與 MCMC 轉移相結合,非常適合流式和回火目標。諸如 SGLD 和 SGHMC 之類的隨機梯度 MCMC 方法通過用小批量估計代替全數據梯度,以犧牲精確平穩性換取可處理性,從而擴展到大數據集。MCMC 也與其他推斷工具互動:它為評估 Variational Inference 提供漸近無偏的參考後驗,為 Gaussian Processes 提供完全貝葉斯的超參數估計,並在 E 步不可處理時為 Expectation-Maximization Algorithm 提供隨機化的替代方案。

局限性

MCMC 的主要弱點是樣本效率:每次抽樣都需要對未歸一化密度進行一次完整評估(對於 HMC 還需要其梯度),而且相鄰抽樣之間是相關的,因此有效樣本量可能比鏈的長度小若干個數量級。在多模態後驗分佈上,混合性能下降——鏈可能在單一模態中被困任意長的時間;在曲率變化劇烈的幾何上,在某一區域可行的步長在另一區域可能失敗。收斂診斷是必要但非充分的。隨機梯度變體會引入很少被嚴格量化的偏差。對於非常高維或帶潛變量的模型,如果校準的不確定性不是首要目標,通常更傾向於使用變分和攤銷方法,而把 MCMC 留作金標準的參考推斷。