【统计学习方法死磕系列Ⅱ】马尔科夫链蒙特卡洛法

一、引言：当积分无法解析计算时

贝叶斯推断中的一个核心难题是计算后验分布的各种期望：

$$ \mathbb{E}_{p(\theta|X)}[f(\theta)] = \int f(\theta) p(\theta|X) d\theta $$

这个积分在低维时可以用数值方法（如梯形法则），但在高维时（如 $\theta$ 有数百个维度）会遭遇”维度的诅咒”——积分离散网格数以指数增长。

马尔科夫链蒙特卡洛（Markov Chain Monte Carlo, MCMC）是解决这一难题的利器。它的核心思想是：不必直接计算积分，而是从目标分布中采样。如果用大量样本 $\theta^{(1)}, \theta^{(2)}, \dots, \theta^{(S)}$ 来近似后验分布，那么期望可以直接用样本平均估计：

$$ \mathbb{E}[f(\theta)] \approx \frac{1}{S} \sum_{s=1}^{S} f(\theta^{(s)}) $$

MCMC 方法通过构造一个以目标分布为平稳分布的马尔可夫链来实现采样，已经成为现代贝叶斯统计和机器学习中不可或缺的工具。

二、Monte Carlo 与 Markov Chain

2.1 Monte Carlo 方法

Monte Carlo 方法的核心是：用随机采样来近似确定性计算。

最简单的是独立采样（如 rejection sampling, importance sampling）——直接从目标分布（或某个易采样的建议分布）中独立抽取样本。

但问题是：对于复杂的高维分布，独立采样几乎不可行。你无法从”完全没有解析形式”的后验分布中直接采样。

2.2 马尔可夫链

一个马尔可夫链由状态空间 $\mathcal{X}$ 和转移核 $T(x’|x)$ 定义。其最关键的性质是：如果转移核满足某些条件（不可约、非周期），则存在唯一的平稳分布 $\pi$：

$$ \pi(x') = \int T(x'|x) \pi(x) dx $$

这意味着无论从哪个初始状态出发，经过足够多的转移步后，链的状态分布将收敛到 $\pi$。

2.3 MCMC 的核心思想

MCMC 将上述两个概念结合：

构造一个以目标分布 $p(x)$ 为平稳分布的马尔可夫链
运行这个链，丢弃初始的”burn-in”样本（未收敛）
用收敛后的样本做 Monte Carlo 估计

MCMC 的精妙之处在于：我们只需要知道目标分布的比例（unnormalized density），不需要归一化常数。因为大多数 MCMC 算法只用到 $\frac{p(x’)}{p(x)}$ 这样的比值，归一化常数被约掉了。

在贝叶斯推断中，这尤其方便——后验分布 $p(\theta|X) \propto p(X|\theta)p(\theta)$，我们只需要似然和先验即可采样。

三、Metropolis-Hastings 算法

3.1 动机

Metropolis-Hastings（MH）算法是所有 MCMC 方法中最根本的。它是由 Metropolis 等人（1953）提出，后由 Hastings（1970）推广的。

3.2 算法流程

MH 算法从当前状态 $x$ 出发，通过以下步骤生成下一个状态：

Algorithm: Metropolis-Hastings
----------------------------------------------
Input: target density p(x) (up to constant), 
       proposal distribution q(x'|x), 
       number of iterations S
Output: samples x^(1), ..., x^(S)

Initialize x^(0)

For t = 0 to S-1:
    1. Propose x' ~ q(x' | x^(t))
    
    2. Compute acceptance probability:
       α = min(1, [p(x') * q(x^(t) | x')] / [p(x^(t)) * q(x' | x^(t))])
    
    3. Accept or reject:
       Draw u ~ Uniform(0, 1)
       If u < α:
           x^(t+1) = x'
       Else:
           x^(t+1) = x^(t)

3.3 接受概率的直观理解

接受概率：

$$ \alpha(x \to x') = \min\left(1, \frac{p(x')}{p(x)} \cdot \frac{q(x|x')}{q(x'|x)}\right) $$

两个部分的乘积：

**$\frac{p(x’)}{p(x)}$**（目标比）：如果 $x’$ 比 $x$ 更”好”（目标密度更高），则倾向于接受。这是算法的”爬山”倾向。
**$\frac{q(x|x’)}{q(x’|x)}$**（Hastings 修正）：如果从 $x$ 很容易到达 $x’$，但从 $x’$ 很难回到 $x$，则需要降低接受概率以保持细致平衡（detailed balance）。

3.4 MH = 接受/拒绝框架

如果把 MH 看作一种”带约束的随机游走”：

建议分布 $q$ 控制探索——向哪个方向走、走多远
接受概率 $\alpha$ 控制利用——只接受”好”的提议

两者配合，确保长期来看，采样点精确地按照 $p(x)$ 分布。

3.5 常用建议分布

| 名称 | $q(x’|x)$ | 特点 |
|——|———–|——|
| 随机游走 MH | $x’ = x + \epsilon, \epsilon \sim \mathcal{N}(0, \sigma^2)$ | 对称，接受率简化为 $p(x’)/p(x)$ |
| 独立 MH | $q(x’|x) = g(x’)$ | 建议独立于当前状态 |
| Adaptive MH | $\sigma$ 随采样自适应调整 | 自动调参 |

3.6 步长的影响

步长太小：几乎所有提议都被接受，但探索缓慢——高度自相关的样本，有效样本量低
步长太大：大多数提议落在低密度区被拒绝——链”粘住”不动
最优接受率：在高斯目标分布下，随机游走 MH 的最优接受率约为 0.234（Roberts et al., 1997）

3.7 Python 实现

import numpy as np

def metropolis_hastings(log_target, x0, sigma, n_iters, n_burnin=1000):
    """
    Random-walk Metropolis-Hastings.
    
    log_target: log of target density (up to const)
    x0: initial state
    sigma: proposal std (scalar or vector)
    """
    x = x0.copy()
    samples = []
    n_accept = 0
    
    for i in range(n_iters):
        # Propose
        x_prop = x + sigma * np.random.randn(*x.shape)
        
        # Accept/reject
        log_alpha = log_target(x_prop) - log_target(x)
        
        if np.log(np.random.rand()) < log_alpha:
            x = x_prop
            n_accept += 1
        
        if i >= n_burnin:
            samples.append(x.copy())
    
    accept_rate = n_accept / n_iters
    return np.array(samples), accept_rate

# Example: sample from a 2D Gaussian
def log_target(x):
    return -0.5 * (x[0]**2 / 2 + x[1]**2 / 0.5)

samples, acc = metropolis_hastings(
    log_target, x0=np.array([0.0, 0.0]), 
    sigma=0.5, n_iters=10000
)
print(f"Acceptance rate: {acc:.3f}")

四、Gibbs 采样

4.1 动机

Gibbs 采样（Geman & Geman, 1984）是 MH 算法的一个特例，专门用于多维分布。当我们可以从每个变量的条件分布 $p(x_i | x_{-i})$ 中直接采样时，Gibbs 采样将 MH 的接受率固定为 1（即总是接受）。

4.2 算法

Algorithm: Gibbs Sampling
----------------------------------------------
Input: conditional distributions p(x_i | x_1, ..., x_{i-1}, x_{i+1}, ..., x_d)
Output: samples

Initialize x^(0) = (x_1^(0), ..., x_d^(0))

For t = 0 to S-1:
    x_1^(t+1) ~ p(x_1 | x_2^(t), x_3^(t), ..., x_d^(t))
    x_2^(t+1) ~ p(x_2 | x_1^(t+1), x_3^(t), ..., x_d^(t))
    ...
    x_d^(t+1) ~ p(x_d | x_1^(t+1), x_2^(t+1), ..., x_{d-1}^(t+1))

Gibbs 采样逐个维度更新，每次从当前所有其他变量条件下的条件分布中采样。

4.3 与 MH 的关系

Gibbs 采样可以看作 MH 的一个特例：提议分布为 $q(x_i’|x) = p(x_i’|x_{-i})$，接受概率恒为 1。

证明：

$$ \alpha = \min\left(1, \frac{p(x_i', x_{-i})}{p(x_i, x_{-i})} \cdot \frac{p(x_i|x_{-i})}{p(x_i'|x_{-i})}\right) $$

由于 $p(x_i’, x_{-i}) = p(x_i’|x_{-i})p(x_{-i})$，代入得 $\alpha = 1$。

4.4 Gibbs 采样的适用条件

Gibbs 采样最适用于以下场景：

可以写出每个变量的条件分布（如共轭模型）
变量之间有较强的依存关系

经典应用包括：LDA 中的主题推断、贝叶斯混合模型、Ising 模型等。

4.5 Gibbs vs MH

方面	Metropolis-Hastings	Gibbs Sampling
接受率	需要调参以达到 ~0.234	恒为 1
条件分布	不需要知道	必须能直接采样
调参	需要调 proposal	无参数
高相关性	可能拒绝很多	每步都移动
模型要求	低（只需密度比值）	高（需要完整的条件分布）

五、Hamiltonian Monte Carlo（HMC）

5.1 动机

MH 和 Gibbs 采样在高维、强相关的分布中表现不佳——它们进行的是随机游走式的探索，有效样本效率低。

Hamiltonian Monte Carlo（HMC, Duane et al. 1987; Neal 2011）通过引入物理学中的哈密尔顿动力学原理，使得采样器能沿着目标分布的等高线进行大幅度跳跃，大幅提高了采样效率。

5.2 物理类比

HMC 将采样问题类比为粒子在势能场中的运动：

位置 $q$：我们想要采样的参数（对应 $\theta$）
动量 $p$：辅助变量，从 $\mathcal{N}(0, M)$ 中采样
势能 $U(q) = -\log p(q)$：目标分布的负对数密度
动能 $K(p) = \frac{1}{2}p^T M^{-1} p$
哈密尔顿量 $H(q, p) = U(q) + K(p)$：总能量

哈密尔顿动力学方程：

$$ \frac{dq}{dt} = \frac{\partial H}{\partial p} = M^{-1}p, \quad \frac{dp}{dt} = -\frac{\partial H}{\partial q} = -\nabla U(q) $$

5.3 Leapfrog 积分

HMC 使用蛙跳法（leapfrog integrator）来数值模拟哈密尔顿动力学：

Leapfrog(q, p, ε, L):
    p = p - (ε/2) * ∇U(q)
    for i = 1..L:
        q = q + ε * M^{-1} * p
        if i < L:
            p = p - ε * ∇U(q)
    p = p - (ε/2) * ∇U(q)
    return q, p

HMC 的接受概率使用 Metropolis 修正：

$$ \alpha = \min(1, \exp(H(q, p) - H(q', p'))) $$

由于数值误差，$H$ 不是完全守恒的，所以需要 Metropolis 步骤修正。

5.4 NUTS：自适应 HMC

No-U-Turn Sampler（NUTS, Hoffman & Gelman 2014）是 HMC 的一种自适应版本，它自动调整蛙跳步数 $L$，当轨迹开始”掉头”时停止。NUTS 是 Stan 概率编程语言中使用的默认采样器。

5.5 HMC vs MH vs Gibbs

方面	MH	Gibbs	HMC
探索方式	随机游走	坐标轴移动	沿梯度”滑行”
高维效率	差	中	好
需要梯度	否	否	是
超参数	proposal 标准差	无	$\varepsilon$ (步长), $L$ (步数)
ESS/iteration	低	中	高

ESS（Effective Sample Size）：有效样本量。由于 MCMC 样本是相关的，$S$ 个 MCMC 样本的信息量相当于 $\text{ESS} \ll S$ 个独立样本。

5.6 Slice Sampling：无需调参的 MCMC

Slice Sampling（Neal, 2003）是另一种经典的 MCMC 方法。它的核心思想是：通过在目标分布 $p(x)$ 的”上方”和”下方”交替采样，自动选择步长，从而消除了对 proposal 标准差的调参需求。

Slice Sampling 引入辅助变量 $y$（”高度”），在增广空间 $(x, y)$ 上采样。联合分布为：

$$ p(x, y) = \begin{cases} 1 / Z & \text{if } 0 \leq y \leq p(x) \\ 0 & \text{otherwise} \end{cases} $$

其中 $Z = \int p(x) dx$ 是归一化常数。$p(x, y)$ 对 $y$ 积分后 $x$ 的边缘分布恰为 $p(x)$。

算法流程：

Algorithm: Slice Sampling
----------------------------------------------
For t = 0 to S-1:
    1. Draw y ~ Uniform(0, p(x^(t)))   // "切片"高度
    2. Find interval (L, R) containing x^(t) 
       where p(x) ≥ y                   // "步出"（stepping out）
    3. Draw x^(t+1) ~ Uniform(L, R)     // "步内采样"
       If p(x^(t+1)) < y: shrink interval and retry

Slice Sampling 的优势：

无 proposal 调参：不需要选择 proposal 分布或步长
对多峰分布适应性好：在步出（stepping out）阶段自动适应局部分布宽度
实现简单

局限：

多维扩展困难：在多维空间中，slice sampling 需要在每个维度独立运行，效率下降
轻尾分布可能效率低：如果目标分布有很长的轻尾，”步出”可能非常耗时

Slice Sampling 在早期贝叶斯计算中广泛使用（如 JAGS、早期的 PyMC），但在 HMC/NUTS 成熟后逐渐被取代。理解它仍然有价值——它的双阶段（步出 + 步内）设计体现了 MCMC 方法设计中的核心权衡：探索（找到高概率区域）vs 利用（在高概率区域内采样）。

5.7 实战：使用 Stan 和 PyMC 进行贝叶斯回归

现代 MCMC 的实际使用通常不手写采样器，而是使用概率编程语言（PPL）。以下是使用 PyMC 进行贝叶斯线性回归的完整示例：

import pymc as pm
import numpy as np

# Simulate data
np.random.seed(42)
X = np.random.randn(100, 2)
true_beta = np.array([2.0, -1.5])
true_sigma = 0.5
y = X @ true_beta + true_sigma * np.random.randn(100)

# Bayesian linear regression with PyMC
with pm.Model() as model:
    # Priors
    beta = pm.Normal('beta', mu=0, sigma=10, shape=2)
    sigma = pm.HalfNormal('sigma', sigma=5)
    
    # Likelihood
    mu = pm.math.dot(X, beta)
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)
    
    # Sample (NUTS is the default sampler)
    trace = pm.sample(draws=2000, tune=1000, chains=4, random_seed=42)

# Posterior summary
summary = pm.summary(trace)
print(summary)

# Check convergence
print(f"R-hat: {pm.rhat(trace)}")
print(f"ESS: {pm.ess(trace)}")

# Posterior predictive
with model:
    post_pred = pm.sample_posterior_predictive(trace)

使用 Stan 的等价代码（在 .stan 文件中）：

data {
    int<lower=0> N;
    int<lower=0> K;
    matrix[N, K] X;
    vector[N] y;
}
parameters {
    vector[K] beta;
    real<lower=0> sigma;
}
model {
    beta ~ normal(0, 10);
    sigma ~ normal(0, 5);
    y ~ normal(X * beta, sigma);
}

关键实践要点：

tune（warm-up）：Stan/PyMC 默认使用 warm-up 阶段自动调整 HMC/NUTS 的步长和 mass matrix。这是 MCMC 自动化的重要一步。
Mass Matrix 自适应：NUTS 在 warm-up 阶段会学习参数之间的相关性（通过估计 posterior 协方差矩阵），从而设置合适的 HMC mass matrix。这解释了为什么 NUTS 在高相关性分布中仍然高效。
Target Acceptance Rate：NUTS 默认目标接受率为 0.8（对于 HMC，最优接受率高于随机游走 MH 的 0.234）。
Divergences：如果出现 divergent transitions，通常意味着步长太大或模型参数化不当。可以使用非中心化参数化（non-centered parameterization）解决。

六、MCMC 诊断与最佳实践

6.1 收敛诊断

$\hat{R}$ 统计量（Gelman-Rubin 统计量）：比较多个独立链的链间方差和链内方差：

$$ \hat{R} = \sqrt{\frac{\widehat{\text{var}}^+(\psi)}{W}} $$

其中 $\widehat{\text{var}}^+$ 是 $\text{pooled}$ 后验方差的估计，$W$ 是链内方差。

当 $\hat{R} < 1.01$ 时，链通常被认为已经收敛。

迹图（Trace plot）：画出参数值随时间的变化，目测是否有趋势或周期。

自相关函数（ACF）：检查样本的相关性衰减速度。

6.2 Burn-in 和 Thinning

Burn-in：丢弃前 $B$ 个样本，以消除初始状态的影响。一般在迹图稳定后即可判定 burn-in 结束
Thinning：每隔 $m$ 个样本保留 1 个，以减少存储和自相关。注意：thinning 不会增加有效信息量，纯属存储优化

6.3 实践清单

MCMC 最佳实践清单：
□ 运行至少 4 条独立链，从不同初始点出发
□ 确保 R̂ < 1.01（最好 < 1.005）
□ 检查 ESS（应 ≥ 100 每链用于基本推断）
□ 目检 trace plots 确保没有趋势
□ 检查是否有发散转移（divergent transitions）
□ 对高维模型考虑使用 HMC/NUTS

七、面试高频问题

Q1：MCMC 与普通的 Monte Carlo 有什么区别？

普通 Monte Carlo 方法尝试从目标分布独立采样（如 rejection sampling, importance sampling）。这在低维可行，但在高维遭遇维度灾难。

MCMC 不要求独立采样，而是构造一个马尔可夫链，其平稳分布等于目标分布。样本虽然是相关的，但只要链足够长，样本平均仍然收敛到真实期望。MCMC 在高维问题中通常更可行。

Q2：MH 算法为什么需要接受/拒绝步骤？不能直接接受所有提议吗？

如果直接接受所有提议，链的平稳分布将由提议分布 $q$ 决定，而非目标分布 $p$。例如，如果 $q$ 倾向于往右走，直接接受所有提议会导致链越来越右，采样不到左侧的区域。

接受/拒绝步骤通过精确的 Metropolis 修正（$\alpha = \min(1, \frac{p(x’)q(x|x’)}{p(x)q(x’|x)})$）保证了细致平衡条件（detailed balance）：

$$ p(x) T(x'|x) = p(x') T(x|x') $$

细致平衡蕴含 $p$ 是链的平稳分布。

Q3：Gibbs 采样为什么总是接受？它的秘密是什么？

Gibbs 采样的提议分布是 $q(x_i’|x) = p(x_i’|x_{-i})$。代入 MH 的接受概率：

$$ \alpha = \min\left(1, \frac{p(x')}{p(x)} \cdot \frac{p(x_i|x_{-i})}{p(x_i'|x_{-i})}\right) $$

而 $p(x’) = p(x_i’|x_{-i}) p(x_{-i})$，$p(x) = p(x_i|x_{-i}) p(x_{-i})$，代入得：

$$ \alpha = \min\left(1, \frac{p(x_i'|x_{-i})p(x_{-i})}{p(x_i|x_{-i})p(x_{-i})} \cdot \frac{p(x_i|x_{-i})}{p(x_i'|x_{-i})}\right) = \min(1, 1) = 1 $$

Gibbs 不”神奇”，它只是选择了使接受率恒为 1 的特殊提议分布。

Q4：HMC 为什么比随机游走 MH 更高效？原理是什么？

随机游走 MH 的探索是无方向的——$\sqrt{n}$ 步后只走了 $\approx \sigma \sqrt{n}$ 的距离（随机游走的均方根位移）。在 $d$ 维中，这极其缓慢。

HMC 利用目标分布的梯度信息，沿着分布的”山谷”（高概率密度方向）进行长距离移动。在物理类比中，粒子沿着等势能面的切线方向滑行（维持了总能量），能快速穿越状态空间——这等价于利用了分布的几何结构。

代价是需要计算梯度 $\nabla \log p(x)$，但对于可微分的模型，现代自动微分工具（如 Stan, PyMC, TensorFlow Probability）使这变得简单。

Q5：如何判断 MCMC 链是否收敛？有哪些诊断手段？

$\hat{R}$（Gelman-Rubin）：比较链间方差和链内方差，$\hat{R} < 1.01$ 表示收敛
迹图：多条链的迹图混合在一起，没有明显趋势
自相关图：lag 增大后自相关快速衰减
有效样本量（ESS）：$\text{ESS} > 100$ per chain
发散转移检测：HMC/NUTS 中，如果蛙跳积分出现数值不稳定（步长太大），会标记为 divergent transition

综合使用这些诊断手段，而非依赖单一指标。

八、总结

MCMC 是现代贝叶斯统计的引擎。它将一个看似不可能的高维积分问题，转化为构造和模拟马尔可夫链的问题。从 Metropolis-Hastings 到 Gibbs 采样再到 Hamiltonian Monte Carlo，MCMC 方法在理论与计算效率上不断进化，使得越来越复杂的贝叶斯模型能够在可接受的时间内进行推断。

理解 MCMC 不仅是学习统计方法的必要步骤，也是理解概率编程语言（如 Stan, PyMC, JAGS）内部工作原理的关键。在后续的学习中，当你看到这些工具快速地对复杂贝叶斯模型进行采样时，你应当意识到：它们背后运行的正是本文介绍的这些优雅算法。

参考文献：

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. (1953). Equation of state calculations by fast computing machines. Journal of Chemical Physics, 21(6), 1087-1092.
Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1), 97-109.
Geman, S., & Geman, D. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6(6), 721-741.
Neal, R. M. (2011). MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo.
Hoffman, M. D., & Gelman, A. (2014). The No-U-Turn Sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. JMLR, 15(1), 1593-1623.