机器学习数学基础 120 章

MCMC、Metropolis–Hastings 与 Gibbs 采样

层级：C｜建议先修：05-22、08-05、08-07

当目标分布只能计算到比例常数、状态空间又太大而无法精确求和时，可以构造以目标分布为平稳分布的 Markov 链，用链上的样本近似期望。这就是 Markov chain Monte Carlo（MCMC）。

1. 目标

设目标密度

π(x)=π~(x)Z,\pi(x)=\frac{\tilde\pi(x)}{Z},

其中 ZZ 难以计算，但未归一化密度 π~(x)\tilde\pi(x) 可求。我们希望估计

Eπ[f(X)]=∫f(x)π(x)dx.\mathbb E_\pi[f(X)] =\int f(x)\pi(x)dx.

若能生成近似来自 π\pi 的样本 X1,…,XTX_1,\ldots,X_T，就用

I^=1T∑t=1Tf(Xt)\hat I=\frac1T\sum_{t=1}^Tf(X_t)

估计该期望。

2. MCMC 的核心设计

构造转移核 P(x′∣x)P(x'\mid x)，使 π\pi 是其平稳分布，并尽量满足不可约、非周期等遍历条件。运行足够久后，链的边缘分布接近 π\pi，时间平均收敛到目标期望。

样本之间相关并不妨碍一致性，但会增大 Monte Carlo 方差。

3. Metropolis–Hastings 算法

当前状态为 xx：

  1. 从提议分布 q(x′∣x)q(x'\mid x) 生成候选 x′x'；
  2. 计算接受概率
α(x,x′)=min{1,π~(x′)q(x∣x′)π~(x)q(x′∣x)};\alpha(x,x')=min\left\{1, \frac{\tilde\pi(x')q(x\mid x')} {\tilde\pi(x)q(x'\mid x)} \right\};
  1. 以概率 α\alpha 接受候选，否则保持在 xx。

归一化常数 ZZ 在比值中相消，这是算法最有价值的性质。

4. 为什么它正确

对 x≠x′x\ne x'，转移概率满足细致平衡：

π(x)q(x′∣x)α(x,x′)=π(x′)q(x∣x′)α(x′,x).\pi(x)q(x'\mid x)\alpha(x,x') =\pi(x')q(x\mid x')\alpha(x',x).

因此 π\pi 是平稳分布。在加上遍历条件后，链可从一般初始状态趋近目标分布。

5. 随机游走 Metropolis

若提议对称，例如

x′=x+ε,qquadε∼N(0,σ2I),x'=x+\varepsilon,qquad \varepsilon\sim\mathcal N(0,\sigma^2I),

则 q(x′∣x)=q(x∣x′)q(x'\mid x)=q(x\mid x')，接受率简化为

α=min⁡{1,π~(x′)π~(x)}.\alpha=\min\left\{1,\frac{\tilde\pi(x')}{\tilde\pi(x)}\right\}.

步长太小会导致高接受率但移动缓慢；步长太大会导致低接受率、长时间停留。接受率本身不是唯一目标，应关注有效样本量和混合情况。

6. 独立提议

若 q(x′∣x)=q(x′)q(x'\mid x)=q(x') 与当前状态无关，则

α=min⁡{1,π~(x′)q(x)π~(x)q(x′)}.\alpha=\min\left\{1, \frac{\tilde\pi(x')q(x)}{\tilde\pi(x)q(x')} \right\}.

提议分布必须充分覆盖目标分布尾部，否则链可能很难访问重要区域。

7. Gibbs 采样

对 x=(x1,ldots,xd)\mathbf x=(x_1,ldots,x_d)，若各个完整条件分布易采样，Gibbs 采样依次更新：

x1∼p(x1∣x2,ldots,xd),x_1\sim p(x_1\mid x_2,ldots,x_d), x2∼p(x2∣x1,x3,ldots,xd),x_2\sim p(x_2\mid x_1,x_3,ldots,x_d),

直到全部坐标。每次更新使用本轮已经更新的最新值。

Gibbs 可视为接受率恒为 1 的特殊 MH，因为候选直接来自条件目标分布。

8. 分块 Gibbs

若变量强相关，逐坐标更新可能移动缓慢。可以将相关变量组成块，从联合条件分布中同时采样。块越大，单步计算越贵，但可能显著改善混合。

9. burn-in 与初始化

链开始阶段仍受初始状态影响，常丢弃前一段样本，称为 burn-in 或 warm-up。但“丢弃固定若干步”不是收敛证明。应结合多链、轨迹图、诊断统计量和领域知识判断。

10. 自相关与有效样本量

MCMC 样本相关，样本均值的方差可用积分自相关时间描述。粗略地，

Neff≈T1+2∑k≥1ρk,N_{\mathrm{eff}} \approx\frac{T}{1+2\sum_{k\ge1}\rho_k},

其中 ρk\rho_k 是滞后 kk 的自相关。有效样本量可能远小于迭代次数。

11. 收敛诊断

常见检查包括：

  • 多个分散初值的链是否混合到相同区域；
  • 轨迹图是否稳定且频繁往返；
  • R^\hat R 是否接近 1；
  • 有效样本量是否足够；
  • 后验统计量的 Monte Carlo 标准误是否可接受。

诊断只能发现部分问题，不能有限步严格证明已收敛。

12. 多峰与高维困难

局部随机游走可能困在一个模式中。改进方法包括重参数化、温度方法、Hamiltonian Monte Carlo、slice sampling、并行回火以及更合适的提议分布。

13. 易错点

  1. 接受率高不代表链好，若步长极小，有效移动仍很少。
  2. 拒绝候选时必须把当前状态再记一次，它仍是链的一个样本。
  3. thinning 会减少存储，却通常不会凭空增加总信息量；优先延长有效运行或改善采样器。
  4. 单链看起来平稳不代表访问了所有模式。

常见问答

Q1：MCMC 样本不是独立的，为什么还能估计期望？
遍历定理保证在适当条件下时间平均仍收敛，只是相关性会降低有效样本量。

Q2：为什么无需知道配分函数？
MH 接受率只用密度比，同一个未知归一化常数在分子分母中相消。

Q3：Gibbs 采样一定比 MH 好吗？
不一定。条件分布可能难采样，强相关时逐坐标 Gibbs 也会很慢；算法选择依赖目标几何结构。

Q4：burn-in 应该设多少？
没有通用固定值。应结合模型、初始化、多链诊断和目标精度评估，而不是机械取迭代数的某个百分比。

练习

  1. 对称提议下推导 MH 接受率的简化形式。
  2. 说明拒绝候选仍需保留当前状态的原因。
  3. 写出二元联合分布 p(x,y)p(x,y) 的一轮 Gibbs 更新。
  4. 若样本自相关很强，有哪些改进方向？

答案与提示

  1. 提议密度比为 1，只剩 min⁡{1,π~(x′)/π~(x)}\min\{1,\tilde\pi(x')/\tilde\pi(x)\}。
  2. 自环概率是保证转移核和目标平稳性的组成部分，删除会改变链的分布。
  3. 先采 xnew∼p(x∣yold)x^{new}\sim p(x\mid y^{old})，再采 ynew∼p(y∣xnew)y^{new}\sim p(y\mid x^{new})。
  4. 调整提议尺度与形状、重参数化、分块更新，或使用 HMC、温度方法等更适合的采样器。