MCMC采样：用马尔可夫链采样

可视化演示

展示马尔可夫链状态转移与平稳分布、MH算法接受拒绝过程、采样直方图收敛诊断，交互切换各维度。

采样困境：看不见全貌的山脉

想象你被派去绘制一座常年大雾的山脉。你想要一张精确的等高线图，于是把山区划分成细密网格，逐格测量海拔。可惜这座山没有现成的“海拔解析式”：你只知道局部的陡峭趋势，不知道如何把整座山的高度统一缩放。在贝叶斯推断里，这座山就是后验分布，而“统一缩放”就是那个令人头疼的归一化常数。

贝叶斯定理告诉我们，后验分布由似然与先验相乘得到：

其中分子可以直接计算，分母却要对所有可能的 θ 做积分，在高维空间中几乎无法精确完成。

传统方法比如拒绝采样，相当于把整座山罩进一个巨大的包围盒，盲目标飞镖。飞镖扎进目标区域就收下，否则重掷。但高概率区域之外的山谷面积太大，飞镖大多浪费；维度越高，浪费越是指数级扩大。与其盲掷，不如让探索者真正“走”起来——脚印密度天然反映山的高度。这正是 MCMC 的思想转变。

马尔可夫链：只记脚下的路

给探索者立一条规矩：他下一步怎么走，只取决于此刻脚踩的位置，完全忘记来时的路。这种“健忘”看似缺点，却是 MCMC 的基石。为什么？因为只要走法设计得好，他会在没有全局视野的情况下，自然地向高处聚集。数学上，这种“无记忆”的行走叫马尔可夫链。

记 X_t 为第 t 步踩到的地点，马尔可夫性质写成：

从 i 一步走到 j 的概率写作 P(i,j)，所有状态两两之间的转移概率排成一张方阵，就是转移矩阵。探索者走很久以后，他在各个地点停留的时间比例会趋于稳定，这个稳定的比例分布 π 满足：

如果我们设计的脚步让平稳分布恰好等于目标分布，那么记录足迹就等于从目标分布采样。关键问题是：怎么设计脚步才保证到达目标？下一节的“细致平衡”就是答案。

细致平衡：往返流量守恒

山间有许多小路连接各个地点。探索者长期游走后，期望上从 i 点出发走向 j 点的次数，应该恰好等于从 j 点走回 i 点的次数。这就像一座古镇里，白天从东门进西门出的人流，与从西门进东门出的人流必须持平，街区才不会越挤越失衡。这个“往返流量守恒”的想法，叫细致平衡。

写成数学式：

把细致平衡条件两边同时对 j 求和，左侧变成 π(i)（因为 Σ_j P(i,j)=1），右侧变成 Σ_j π(j)P(j,i)，恰好推出平稳分布方程。所以细致平衡是平稳分布的充分条件，虽然不是必要条件。MH 算法和吉布斯采样正是从“如何构造细致平衡”出发设计的。

Metropolis-Hastings：提议与裁决

探索者站在 i 点，该怎么迈出下一步？他先闭上眼睛，按一个简单规则“乱指一个方向”——这个方向只和当前位置有关，叫作提议分布。指完后他要做一次裁决：新地点到底值不值得去？如果明显比现在好，就去；如果差一些，掷硬币碰碰运气，偶尔也允许下山。这样既能多在高处停留，又不会永远困在一座山峰上。

Metropolis-Hastings（MH）把裁决写成接受率：

算法流程非常直白：先按 q 提议一个新点 j；再按上面公式算 α；接着抽一个 [0,1] 均匀随机数 u，若 u < α 就跳到 j，否则留在原地；最后无论是否移动，都把当前位置记作一个样本。注意“原地不动”也计入足迹，所以在样本序列里会出现连续重复的脚印——这没关系，它们代表了对当前位置的强调。

MH 是“先乱指、后裁决”。只要提议分布满足基本条件，这个简单的博弈过程就会让探索者的长期时间分配精确对齐目标分布。即便目标分布只能算到未归一化的形状，也没问题，因为分子分母中的归一化常数恰好约掉了。

吉布斯采样：一次只走一个方向

与其闭眼乱指，不如沿规定的山路走：探索者一次只改变一个坐标。在二维平面上，这像走“阶梯路径”——先水平移动到合适位置，再垂直移动到合适位置，反复交替。每一步都严格沿着东西或南北方向，不斜跨。

吉布斯采样要求我们能写出“给定其他所有变量时、某个变量所服从的条件分布”，然后直接从中采样：

可以证明，吉布斯采样是 MH 的特例：它的提议分布恰好是目标分布的全条件分布，因此接受率恒等于 1。它高效、无需调参，但前提是条件分布容易采样；如果某个维度的条件分布太复杂，就只好回到 MH。工程中常把两种策略混用，能用吉布斯的维度就用吉布斯，不能用的交给 MH，这辆车既能走阶梯路，也能跳跃前进。

收敛诊断：探索者真的逛遍了吗

探索者走了一天又一天，脚印铺满山坡。可你心里打鼓：他真的逛遍了整片山，还是困在某座偏峰上反复转圈？光看脚印多不多不够，得用三个工具自查。

第一是烧入期（burn-in）。探索者是被随机空投的，开头那段路受起点位置影响极大，不能代表山脉的真实地形。通常把最前面若干步丢掉，就像锅炉先烧热再读数，预热阶段的数字不算数。第二是自相关。相邻脚印不可能独立：上一步在哪，下一步八九不离十。相隔 k 步的脚印之间的相似程度用 ρ_k 衡量，自相关越强，样本携带的新信息越少。有效样本量（ESS）把这个损失量化出来：

第三是多链对比。从几个相隔很远的起点同时放出多支探索队伍，比较“队伍之间的差异”与“单支队伍内部的差异”。若所有队伍最终都汇合到同一片区域，说明链已收束；工程上常用 Gelman-Rubin 的 R-hat 指标，它约接近 1 越好。烧入期清掉预热垃圾，ESS 告诉你别高估脚印的独立性，R-hat 判断多支队伍有没有会师——三项都过关，足迹地图才敢放心用。

贝叶斯应用：从足迹到地图

探索者终于逛完了山，现在要把脚印整理成一张正式地图。在贝叶斯工作流里，这张地图就是后验分布，而“整理”的核心动作是用后验样本做预测。拿到一串后验样本 θ^(1),…,θ^(N) 后，对任意新观测 x̃，预测分布要把模型在每一个可能参数下的预报按后验权重平均：

概率编程库（PyMC、Stan、NumPyro）正是这么工作的：我们只写模型代码描述先验与似然，库内部自动调度 MCMC 得到后验样本。返回值不是单点估计，而是一堆样本——可以画直方图、算分位数、汇报最高密度区间。在医疗、金融、自动驾驶等领域，这种“不仅给预测，还给置信度”的能力无可替代。

山还是那座山，但我们已用足迹读懂它的轮廓。MCMC 把数学上算不出的后验积分，变成一次次的行走与记录；足迹不是地图本身，但已经足够让我们看清山的形状。

AI 应用场景

Python 演示

import numpy as npimport matplotlib.pyplot as plt def target_logpdf(x, y):    g1 = -((x - 0.35) ** 2 + (y - 0.65) ** 2) / 0.012    g2 = -((x - 0.70) ** 2 + (y - 0.30) ** 2) / 0.012    return np.logaddexp(g1, g2) rng = np.random.default_rng(42)x, y = 0.5, 0.5samples = []n_steps = 20000burn_in = 2000accepted = 0 for t in range(n_steps):    x_prop = x + rng.normal(0, 0.05)    y_prop = y + rng.normal(0, 0.05)    log_alpha = target_logpdf(x_prop, y_prop) - target_logpdf(x, y)    if np.log(rng.random()) < log_alpha:        x, y = x_prop, y_prop        accepted += 1    if t >= burn_in:        samples.append((x, y)) samples = np.array(samples)print("accept rate: {:.2%}".format(accepted / n_steps))print("sample mean: ({:.3f}, {:.3f})".format(samples[:, 0].mean(), samples[:, 1].mean())) plt.figure(figsize=(6, 5))plt.scatter(samples[:, 0], samples[:, 1], s=2, alpha=0.4, color="#818CF8")plt.title("MH sampler on bimodal distribution")plt.xlim(0, 1)plt.ylim(0, 1)plt.show()

1. return np.logaddexp(g1, g2)用 log-sum-exp 合并两座山峰的对数高度，避免概率乘到极小后数值下溢。
2. x_prop = x + rng.normal(0, 0.05)对称高斯提议：在当前位置附近随机迈一小步，q(i,j)=q(j,i)。
3. log_alpha = target_logpdf(x_prop, y_prop) - target_logpdf(x, y)在对数域计算接受率，等价于 min(1, π(j)/π(i))。
4. if np.log(rng.random()) < log_alpha:取随机数再取对数后比较，等价于以概率 min(1, e^{log_alpha}) 接受移动。
5. if t >= burn_in: samples.append((x, y))丢弃 burn-in 前的预热足迹，只保留稳定段样本用于统计。

常见误区

面试题