4. Metropolis-Hastings算法原理

Metropolis-Hastings（简称MH）算法是MCMC方法族中最具代表性的一种构造性采样方法。它的核心思想是：当我们无法直接从目标分布 p(x) 中抽样时，可以借助一个容易采样的提议分布 q(x'|x)，构造一条马尔可夫链，使该链的平稳分布恰好等于目标分布。理解MH的关键在于三件事：提议分布如何选、接受率如何推导、算法如何落地实现。

一、从"提议"到"接受"：MH的构造逻辑

设当前状态为 x，我们从提议分布 q(x'|x) 中生成候选状态 x'。如果直接接受所有候选，链的平稳分布将是 q 而非 p，因此必须引入一个接受概率 α(x, x') 来"修正"提议分布带来的偏差。MH的接受率定义为：

α(x, x') = min{ 1, [p(x') · q(x | x')] / [p(x) · q(x' | x)] }

这个式子的直观含义是：如果候选状态在目标分布下更"合理"（p(x') ≥ p(x)），则一定接受；如果候选状态更"不合理"，则以一定概率接受，接受概率由目标分布比值与提议分布反向比值的乘积决定。注意这里出现了 q(x|x') 与 q(x'|x) 的比值，这正是MH相比原始Metropolis算法（要求提议分布对称）的推广之处——它允许使用非对称提议分布。

二、接受率的推导思路

接受率的推导依赖细致平衡条件（detailed balance）：p(x)·T(x→x') = p(x')·T(x'→x)，其中 T(x→x') = q(x'|x)·α(x,x') 是转移核。将转移核代入细致平衡式，得到：

p(x) · q(x'|x) · α(x, x') = p(x') · q(x | x') · α(x', x)

若令 α(x,x') = min{1, R(x,x')}，其中 R = [p(x')q(x|x')]/[p(x)q(x'|x)]，可以验证当 R ≥ 1 时 α=1，反向接受率 α(x',x)=1/R，等式成立；当 R < 1 时 α=R，反向接受率为1，等式同样成立。因此该接受率恰好满足细致平衡，从而保证 p 是链的平稳分布。实际实现时，由于 α 中只出现 p(x')/p(x) 的比值，目标分布中的归一化常数会被约掉——这正是MH能处理未归一化密度的根本原因。

三、算法流程与实现要点

将上述逻辑整理为可执行的迭代流程：初始化 x₀；对第 t 步，从 q(·|x_{t-1}) 采样候选 x'；计算接受率 α；生成均匀随机数 u~U(0,1)；若 u ≤ α 则接受 x_t = x'，否则 x_t = x_{t-1}（注意拒绝时也要把当前状态重复记入样本链，这是保证平稳分布正确的关键，不能跳过）。

实现时有几个容易踩坑的细节。第一，提议分布的标准差（步长）直接决定接受率与混合速度：步长过小则接受率高但链移动缓慢、自相关强；步长过大则大量候选被拒、链停滞。随机游走MH的经验目标是接受率落在0.2～0.5之间。第二，必须丢弃前若干步的"预热期"（burn-in）样本，因为链从初始点出发尚未收敛到平稳分布。第三，为降低自相关，可对链进行稀释（thinning），每隔若干步取一个样本。

四、手写实现一维目标分布采样

下面以目标分布 p(x) ∝ exp(-x²/2)·(1 + 0.5·sin(3x)) 为例（一个被轻微调制的标准正态，未归一化），用对称正态随机游走提议 q(x'|x)=N(x, σ²) 实现MH。由于提议对称，q(x|x')/q(x'|x)=1，接受率简化为 min{1, p(x')/p(x)}。

import numpy as np

def target_logpdf(x):
    # 未归一化对数密度，避免数值下溢
    return -0.5 * x**2 + np.log(1 + 0.5 * np.sin(3 * x))

def mh_sample(n_iter=20000, sigma=1.2, x0=0.0, burn_in=2000, thin=5, seed=42):
    rng = np.random.default_rng(seed)
    x = x0
    samples = []
    accepted = 0
    for t in range(n_iter):
        x_prop = x + rng.normal(0, sigma)          # 对称提议
        log_r = target_logpdf(x_prop) - target_logpdf(x)
        if np.log(rng.random()) < log_r:           # 等价于 u <= alpha
            x = x_prop
            accepted += 1
        if t >= burn_in and (t - burn_in) % thin == 0:
            samples.append(x)
    return np.array(samples), accepted / n_iter

samples, acc_rate = mh_sample()
print("接受率:", round(acc_rate, 3))
print("样本均值:", round(samples.mean(), 3), "样本标准差:", round(samples.std(), 3))

这里用对数密度之差与 log(u) 比较，等价于 u ≤ α，但避免了极小概率下的下溢问题，是工程实现的标准做法。以 σ=1.2 运行时，接受率通常落在0.4～0.6附近（示例参数下的合理预期，实际数值随随机种子略有波动）。若接受率明显偏高（如0.9以上），说明步长偏小，应增大 σ；若低于0.15，说明步长过大，应减小 σ。

五、结果判读与收敛检查

得到样本后，不能直接当作目标分布的样本使用，需要做基本诊断。可将样本直方图与目标密度曲线叠加对比，观察形状是否吻合；绘制样本轨迹图（trace plot），看链是否在均值附近稳定波动而无明显漂移；计算样本自相关函数，若滞后若干步后自相关迅速衰减到接近0，说明混合良好。对于本例，样本均值应接近0附近，直方图应呈现带轻微波动的钟形。若轨迹图显示链长期卡在某一区域，往往意味着提议步长不当或目标分布存在多峰，此时可考虑调整步长或改用自适应MH。