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 留作金标准的参考推断。