无监督学习 · 05

马尔可夫链蒙特卡罗法

当积分算不动、后验写不出，MCMC 用一条精心设计的随机游走链，从任意复杂分布中采出样本——蒙特卡罗遇上马尔可夫链的魔法。

18 min read

为什么需要采样

上一篇讨论了 EM 算法——它需要对隐变量的后验做期望。当后验分布足够简单（比如 GMM 的离散责任值），期望可以解析地算出来。但更多时候，这个期望是一个高维积分，根本写不出闭合形式。

考虑贝叶斯推断的核心问题：给定数据 DD，计算参数 θ\theta 的后验：

p(θ∣D)=p(D∣θ) p(θ)p(D)=p(D∣θ) p(θ)∫p(D∣θ) p(θ) dθp(\theta \mid D) = \frac{p(D \mid \theta) \, p(\theta)}{p(D)} = \frac{p(D \mid \theta) \, p(\theta)}{\int p(D \mid \theta) \, p(\theta) \, d\theta}

分母 p(D)p(D) 是一个对所有可能 θ\theta 的积分。参数空间一旦超过几维（神经网络有百万维），这个积分就彻底算不了。

蒙特卡罗法的思路：既然积分算不动，不如从 p(θ∣D)p(\theta \mid D) 里抽一批样本 {θ(1),…,θ(N)}\{\theta^{(1)}, \dots, \theta^{(N)}\}，然后用样本均值代替期望：

Ep(θ∣D)[f(θ)]≈1N∑i=1Nf(θ(i))\mathbb{E}_{p(\theta|D)}[f(\theta)] \approx \frac{1}{N} \sum_{i=1}^{N} f(\theta^{(i)})

大数定律保证：只要样本足够多，这个近似可以任意精确。问题变成了——如何从一个只知道未归一化密度 p(D∣θ)p(θ)p(D|\theta)p(\theta) 的复杂分布中采样？

简单的拒绝采样（rejection sampling）在高维空间效率指数级下降。均匀网格（quadrature）在 10 维以上就不可行。这正是 MCMC 登场的地方。

马尔可夫链基础

MCMC 的"MC"之一是马尔可夫链（Markov Chain）。直觉很简单：一条链在状态空间上游走，下一步只取决于当前位置，不关心历史。

定义：随机过程 {X0,X1,X2,… }\{X_0, X_1, X_2, \dots\} 满足马尔可夫性：

P(Xt+1=j∣Xt=i,Xt−1,…,X0)=P(Xt+1=j∣Xt=i)=TijP(X_{t+1} = j \mid X_t = i, X_{t-1}, \dots, X_0) = P(X_{t+1} = j \mid X_t = i) = T_{ij}

矩阵 T=[Tij]T = [T_{ij}] 叫做转移矩阵——它完整描述了链的行为。

平稳分布：如果存在分布 π\pi 满足

πT=π\pi T = \pi

即 π\pi 是 TT 的左特征向量（特征值为 1），那么一旦链达到 π\pi，之后每一步的边际分布都不再变化。π\pi 就叫做链的平稳分布。

核心定理：如果链是不可约（从任何状态都能到达任何状态）且非周期（不会在几个状态间固定轮转）的，那么：

  1. 平稳分布 π\pi 唯一存在。
  2. 无论从哪个初始状态出发，链都会收敛到 π\pi。

这意味着：只要设计一条转移矩阵，让它的平稳分布恰好等于我们的目标分布 p(θ∣D)p(\theta|D)，然后沿着链跑足够多步，后面的样本就近似来自目标分布。

细致平衡条件（detailed balance）是构造这种链的常用充分条件：

π(i) Tij=π(j) Tji∀ i,j\pi(i) \, T_{ij} = \pi(j) \, T_{ji} \quad \forall \, i, j

直觉：从 ii 流向 jj 的"概率流"等于从 jj 流回 ii 的。满足细致平衡的链，π\pi 自动是平稳分布（对 jj 求和即可验证 πT=π\pi T = \pi）。Metropolis-Hastings 和吉布斯抽样都利用细致平衡来保证正确性。

Metropolis-Hastings 算法

Metropolis-Hastings（MH）是最通用的 MCMC 算法。它的精妙之处在于：只需要知道目标分布的未归一化密度，就能构造出满足细致平衡的转移核。

算法步骤：

  1. 初始化 x(0)x^{(0)}，选择提议分布 q(x′∣x)q(x' \mid x)（通常是以当前位置为中心的高斯）。
  2. 在第 tt 步：
    • 从 q(⋅∣x(t))q(\cdot \mid x^{(t)}) 中采一个候选点 x′x'。
    • 计算接受概率：
α(x(t),x′)=min⁡(1,p(x′) q(x(t)∣x′)p(x(t)) q(x′∣x(t)))\alpha(x^{(t)}, x') = \min\left(1, \frac{p(x') \, q(x^{(t)} \mid x')}{p(x^{(t)}) \, q(x' \mid x^{(t)})}\right)
  • 掷一枚均匀硬币 u∼U(0,1)u \sim U(0,1)：
    • 若 u<αu < \alpha，接受：x(t+1)=x′x^{(t+1)} = x'。
    • 否则拒绝：x(t+1)=x(t)x^{(t+1)} = x^{(t)}（停在原地）。

注意 α\alpha 的分子分母里 pp 出现的是比值 p(x′)/p(x(t))p(x')/p(x^{(t)})——归一化常数完美抵消了。这就是 MH 能在只知道未归一化密度时工作的原因。

如果提议分布是对称的（q(x′∣x)=q(x∣x′)q(x'|x) = q(x|x')，比如零均值高斯），接受概率简化为：

α=min⁡(1,p(x′)p(x(t)))\alpha = \min\left(1, \frac{p(x')}{p(x^{(t)})}\right)

直觉：往高概率区域走一定接受；往低概率区域走有一定概率接受——概率越低拒绝越频繁。这让链最终在高概率区域聚集，同时偶尔探索低概率区域。

提议分布的方差是关键调参：

  • 方差太大：候选点经常跑到低概率区域，频繁拒绝，链原地踏步。
  • 方差太小：每步移动太小，需要极多步才能遍历整个分布。
  • 经验法则：接受率在 20%--50% 之间通常效率最高。

下面的交互演示了 MH 采样的过程——调节提议方差，观察接受率和轨迹的变化：

0.8
总采样: 0接受: 0接受率: 0%接受拒绝
起点 (0, 0)
Metropolis-Hastings 采样：目标分布是两个高斯的混合（青色热力图）。点击「采样」逐步观察链的游走轨迹。调节提议方差 σ 观察接受率的变化——太大会频繁拒绝，太小会原地踏步。

吉布斯抽样

吉布斯抽样（Gibbs Sampling）是 MH 的一个特例，适用于条件分布容易采样的场景——这在图模型和贝叶斯模型中非常常见。

思路：假设 x=(x1,x2,…,xd)x = (x_1, x_2, \dots, x_d)，不直接从联合分布 p(x)p(x) 采样，而是逐个分量从条件分布中采：

  1. 初始化 x(0)=(x1(0),x2(0),…,xd(0))x^{(0)} = (x_1^{(0)}, x_2^{(0)}, \dots, x_d^{(0)})。
  2. 在第 tt 步，依次更新每个分量：
x1(t+1)∼p(x1∣x2(t),x3(t),…,xd(t))x_1^{(t+1)} \sim p(x_1 \mid x_2^{(t)}, x_3^{(t)}, \dots, x_d^{(t)}) x2(t+1)∼p(x2∣x1(t+1),x3(t),…,xd(t))x_2^{(t+1)} \sim p(x_2 \mid x_1^{(t+1)}, x_3^{(t)}, \dots, x_d^{(t)}) ⋮\vdots xd(t+1)∼p(xd∣x1(t+1),x2(t+1),…,xd−1(t+1))x_d^{(t+1)} \sim p(x_d \mid x_1^{(t+1)}, x_2^{(t+1)}, \dots, x_{d-1}^{(t+1)})

每一步只采一个维度，固定其余维度。注意后面的分量用的是已经更新过的值，这让收敛更快。

吉布斯是 MH 的特例：把"从条件分布采样"看成一个提议，可以证明这个提议的接受概率恒为 1——每一步都接受。这意味着吉布斯抽样没有浪费的步骤，效率通常比一般的 MH 高。

证明很简洁：对分量 xjx_j，提议分布就是条件分布 q(xj′∣x−j)=p(xj′∣x−j)q(x_j' \mid x_{-j}) = p(x_j' \mid x_{-j})。代入 MH 的接受概率：

α=min⁡(1,p(xj′,x−j)⋅p(xj∣x−j)p(xj,x−j)⋅p(xj′∣x−j))=min⁡(1,p(xj′∣x−j)p(x−j)⋅p(xj∣x−j)p(xj∣x−j)p(x−j)⋅p(xj′∣x−j))=1\alpha = \min\left(1, \frac{p(x_j', x_{-j}) \cdot p(x_j \mid x_{-j})}{p(x_j, x_{-j}) \cdot p(x_j' \mid x_{-j})}\right) = \min\left(1, \frac{p(x_j' \mid x_{-j}) p(x_{-j}) \cdot p(x_j \mid x_{-j})}{p(x_j \mid x_{-j}) p(x_{-j}) \cdot p(x_j' \mid x_{-j})}\right) = 1

局限：当分量之间高度相关时，逐维采样的效率很低——链在一个窄长的"山谷"里缓慢蠕动。这时需要分块（block Gibbs）或切换到其他方法。

MCMC 与机器学习

MCMC 在机器学习中的典型应用场景：

1. LDA 话题模型

潜在狄利克雷分配（LDA）为每篇文档分配话题分布，为每个话题分配词分布。精确后验不可解，标准训练方法就是吉布斯抽样——交替采样「每个词属于哪个话题」这个离散隐变量。条件分布是一个简单的多项式分布，所以吉布斯特别合适。

2. 贝叶斯神经网络

标准神经网络只学一组点估计的权重。贝叶斯神经网络（BNN）对权重维护一个完整的后验分布。推断时对后验采样，多组权重的预测取平均——天然给出不确定性。MCMC（特别是 Hamiltonian Monte Carlo，一种利用梯度信息加速的 MH 变体）是训练 BNN 的经典方法。

3. 贝叶斯超参数调优

把超参数看作随机变量，用 MCMC 从它们的后验中采样，比网格搜索高效得多。贝叶斯优化（Bayesian Optimization）的代理模型——高斯过程——本身的推断也依赖 MCMC。

4. 统计物理与能量模型

玻尔兹曼机（Boltzmann Machine）和受限玻尔兹曼机（RBM）的训练核心就是从能量函数定义的分布中采样。对比散度（Contrastive Divergence）算法本质上是"只跑几步"的吉布斯抽样。

Burn-in 与实践注意事项：

链的前若干步还没有收敛到平稳分布，这段叫做 burn-in，采到的样本应该丢弃。判断收敛的常用方法：

  • 跑多条独立链，看它们的统计量是否一致（Gelman-Rubin 诊断，R^≈1\hat{R} \approx 1）。
  • 观察 trace plot：参数的采样值应该像噪声一样上下波动，而不是有明显趋势。

此外，相邻样本有相关性（因为链的马尔可夫性），有效样本量小于总采样量。常用做法是间隔采样（thinning）：每隔 kk 步取一个样本。

这个想法在前沿里