1. 为什么需要MCMC：从贝叶斯推断的计算瓶颈说起
做贝叶斯推断的人迟早会撞上一堵墙：后验分布往往没有解析解。你写好先验、写好似然，贝叶斯公式一摆，分母那个归一化常数——也就是对全参数空间积分得到的边缘似然——直接把路堵死了。参数少的时候还能数值积分凑合，参数一多，维度一高，网格积分的计算量呈指数爆炸，别说跑模型，连内存都放不下。
我在早期做A/B测试的贝叶斯建模时就吃过这个亏。当时参数只有四五个，我天真地以为数值积分能扛住，结果把网格从100格加密到200格，计算时间翻了不止一倍，精度却没提升多少。后来才彻底转向MCMC（Markov Chain Monte Carlo，马尔可夫链蒙特卡洛）这条路。
MCMC的核心思想其实很朴素：既然我没法直接算出后验分布的解析形式，那我就构造一条马尔可夫链，让它的平稳分布恰好等于我想要的後验分布。然后我从某个初始点出发，让链跑足够长的时间，等它“忘记”初始状态、进入平稳状态后，把链上访问过的样本收集起来，就相当于从后验分布里抽样了。有了这些样本，均值、方差、分位数、可信区间，全都能估计出来。
这个专题是贝叶斯网系列的第8篇，前面几篇我们聊了贝叶斯网的结构、条件独立、推断算法。这一篇专门啃MCMC的抽样原理。我会把Metropolis-Hastings、Gibbs采样、收敛诊断这些核心内容拆开讲透，配上可复现的代码和参数选择的经验。适合已经了解贝叶斯基本概念、想搞明白MCMC到底怎么运转的读者。如果你之前只是调包跑过PyMC或Stan，但说不清里面发生了什么，这篇就是写给你的。
2. MCMC的核心机制拆解：平稳分布与细致平衡
2.1 马尔可夫链为什么能“记住”目标分布
马尔可夫链是一个状态序列，下一个状态只依赖当前状态，跟更早的历史无关。这个性质叫马尔可夫性。听起来像是“没有记忆”，但恰恰是这种结构让链的长期行为变得可分析。
一条马尔可夫链如果满足一定条件（不可约、非周期、正常返），它会有一个平稳分布。所谓平稳分布，就是链跑到后面，处于各个状态的概率不再随时间变化。MCMC的巧妙之处在于：我不需要知道后验分布的具体形式，我只需要设计一个转移规则，让这条链的平稳分布正好是目标后验分布。
打个比方。你想测量一个湖里鱼的长度分布，但湖水太浑，你没法一次性捞上来统计。你能做的是：站在湖边，每次随机撒一网，记录捞到的鱼的长度，然后把鱼放回去。如果你撒网的方式设计得当，长期下来，你记录到的长度分布就会逼近湖里真实的长度分布。MCMC里的“撒网规则”就是转移核，“长期记录”就是链的迭代。
2.2 细致平衡条件：构造转移核的黄金法则
要让目标分布π(x)成为链的平稳分布，一个充分条件是满足细致平衡条件：
π(x) T(x→y) = π(y) T(y→x)
这个式子的意思是：从x转移到y的概率通量，等于从y转移到x的概率通量。如果每条边都满足这个条件，那么整体上分布就不会变化，π就是平稳分布。
细致平衡的好处在于，它把“验证平稳分布”这个全局问题，简化成了“验证每对状态之间的通量平衡”这个局部问题。Metropolis-Hastings算法就是围绕这个条件设计的。
具体来说，我们有一个建议分布q(y|x)，它负责提出候选状态。但候选状态不能照单全收，得按一定概率接受。接受率α(x,y)设计成：
α(x,y) = min(1, [π(y)q(x|y)] / [π(x)q(y|x)])
这个式子的精妙之处在于，π(x)和π(y)的比值消掉了归一化常数。我们不需要知道后验分布的具体数值，只需要知道它正比于什么。先验乘以似然，这个乘积是能算出来的，归一化常数就让它待在分母里不用管。
注意：接受率公式里的π只需要知道到常数倍的程度。这是MCMC能在贝叶斯推断里大杀四方的根本原因。你永远不需要计算那个该死的边缘似然。
2.3 建议分布的选择：随机游走与独立采样
建议分布q的选择直接影响链的效率和混合速度。常见的有两种：
随机游走建议分布：q(y|x) = q(|y-x|)，比如以当前点为中心的正态分布。这是最常用的方案，实现简单，但步长需要调。步长太小，链走得慢，自相关高；步长太大，接受率低，链老是原地踏步。
独立采样建议分布：q(y|x) = q(y)，建议分布不依赖当前状态。这种方案理论上可以做到接受率很高，但如果q跟目标分布差太远，链会卡在某个区域出不来。
我个人的经验是，高维问题优先用随机游走，配合自适应步长；低维问题可以试试独立采样，但要对建议分布做足够的调优。下面这张表总结了两种方案的适用场景：
| 建议分布类型 | 优点 | 缺点 | 适用场景 |
|---|---|---|---|
| 随机游走 | 实现简单，对目标分布形状要求低 | 步长敏感，高维混合慢 | 中高维参数空间，目标分布复杂 |
| 独立采样 | 接受率可调高，样本独立性好 | 需要好的建议分布，否则链卡死 | 低维问题，已有近似后验可用 |
| 混合方案 | 兼顾探索与利用 | 实现复杂，参数多 | 多峰分布，需要跨峰探索 |
3. Metropolis-Hastings算法实操：从伪代码到可运行实现
3.1 算法流程的逐步拆解
Metropolis-Hastings（简称MH）是MCMC的基石算法。它的流程可以拆成五步：
- 初始化：选一个起点x₀，要求目标分布在x₀处的密度大于零。
- 提议：从建议分布q(y|xₜ)中抽一个候选点y。
- 计算接受率：α = min(1, [π(y)q(xₜ|y)] / [π(xₜ)q(y|xₜ)])。
- 接受/拒绝：从均匀分布U(0,1)抽一个数u，如果u < α，则xₜ₊₁ = y；否则xₜ₊₁ = xₜ。
- 重复：回到第2步，直到收集到足够多的样本。
这个流程里，第3步是核心。如果建议分布是对称的，即q(y|x) = q(x|y)，那么接受率简化为α = min(1, π(y)/π(x))。这就是原始的Metropolis算法。
3.2 用Python实现一个完整的MH采样器
下面这段代码实现了一个针对一维正态后验的MH采样器。目标分布是N(3, 2²)，建议分布是以当前点为中心、步长为1的正态分布。
import numpy as np import matplotlib.pyplot as plt def target_logpdf(x): """目标分布的对数密度：N(3, 2^2)""" return -0.5 * ((x - 3) / 2) ** 2 def metropolis_hastings(n_samples, init_x, step_size): samples = np.zeros(n_samples) x_current = init_x n_accepted = 0 for i in range(n_samples): # 提议：对称随机游走 y = x_current + np.random.normal(0, step_size) # 计算对数接受率 log_alpha = target_logpdf(y) - target_logpdf(x_current) # 接受/拒绝 if np.log(np.random.uniform()) < log_alpha: x_current = y n_accepted += 1 samples[i] = x_current accept_rate = n_accepted / n_samples return samples, accept_rate # 运行采样 np.random.seed(42) samples, acc_rate = metropolis_hastings(n_samples=20000, init_x=0.0, step_size=1.0) print(f"接受率: {acc_rate:.3f}") print(f"样本均值: {samples[5000:].mean():.3f} (理论值: 3.000)") print(f"样本标准差: {samples[5000:].std():.3f} (理论值: 2.000)")跑下来你会发现，接受率大概在0.7左右，样本均值和标准差都接近理论值。但这里有个细节：我丢弃了前5000个样本，这就是所谓的“预热期”（burn-in）。预热期的样本还没收敛到平稳分布，不能用来做推断。
3.3 步长调优：接受率不是越高越好
很多人第一次调MH，看到接受率低就拼命缩小步长，结果接受率上去了，但样本自相关极高，有效样本量反而下降。这里有个反直觉的权衡。
步长太小：链每次只挪一点点，相邻样本高度相关，虽然接受率高，但探索效率低。步长太大：候选点经常跳到低密度区，被拒绝，链原地踏步，接受率低，探索效率也低。
理论上的最优接受率，对于一维对称随机游走是0.44左右，高维情况下降到0.23左右。这是Roberts和Rosenthal在1997年给出的经典结果。我实测下来，一维问题把接受率控制在0.4到0.5之间，高维问题控制在0.2到0.3之间，混合效果通常不错。
实操心得：不要死盯接受率。更靠谱的做法是看迹图（trace plot）和自相关函数（ACF）。迹图应该像一条毛线，快速上下波动，没有明显的趋势或周期性。ACF应该快速衰减到零附近。如果ACF衰减很慢，说明步长偏小；如果迹图经常长时间停在同一个值，说明步长偏大。
3.4 预热期的处理与样本收集
预热期要丢多少样本？这个问题没有标准答案，取决于初始点和目标分布的差距。我的经验是：
- 如果初始点是通过某种近似方法（比如拉普拉斯近似）选的，预热期可以短一些，几千个样本就够。
- 如果初始点是随便选的，预热期可能要几万甚至几十万。
- 更稳妥的做法是跑多条链，从不同初始点出发，看它们是否收敛到同一个区域。如果多条链的迹图混在一起分不清，说明预热期够了。
收集样本时，还有一个技巧叫“稀释”（thinning）：每隔k个样本取一个，降低自相关。但现代观点认为，如果你有足够的计算资源，不如直接多跑样本，而不是稀释。稀释会丢弃信息，除非存储空间真的不够。
4. Gibbs采样：当条件分布可以解析采样时
4.1 Gibbs采样的核心思想
Gibbs采样是MH的一个特例，适用于多维参数问题。它的巧妙之处在于：不需要设计建议分布，也不需要计算接受率。它每次只更新一个参数，从该参数在给定其他参数下的条件分布中直接采样。
假设参数向量θ = (θ₁, θ₂, ..., θₚ)，Gibbs采样的一轮迭代是这样的：
- 从p(θ₁ | θ₂, ..., θₚ, 数据)中抽θ₁
- 从p(θ₂ | θ₁, θ₃, ..., θₚ, 数据)中抽θ₂
- ...
- 从p(θₚ | θ₁, ..., θₚ₋₁, 数据)中抽θₚ
每一轮都更新所有参数。因为每次都是从精确的条件分布中采样，所以接受率恒为1，所有样本都被接受。
4.2 条件分布推导的实操案例
用一个简单的贝叶斯线性回归来演示。模型是：
yᵢ = β₀ + β₁xᵢ + εᵢ，εᵢ ~ N(0, σ²)
先验设成：
- β₀ ~ N(0, 100)
- β₁ ~ N(0, 100)
- σ² ~ InverseGamma(1, 1)
在这个设定下，β₀和β₁的联合条件分布是正态分布，σ²的条件分布是逆伽马分布。推导过程涉及配平方，这里不展开，直接给结果：
β的条件后验：β | σ², 数据 ~ N(μ_β, Σ_β) 其中Σ_β = (XᵀX/σ² + I/100)⁻¹，μ_β = Σ_β Xᵀy/σ²
σ²的条件后验：σ² | β, 数据 ~ InverseGamma(1 + n/2, 1 + ||y - Xβ||²/2)
4.3 Gibbs采样的Python实现
import numpy as np def gibbs_sampling(X, y, n_samples=10000, burn_in=2000): n, p = X.shape samples = np.zeros((n_samples, p + 1)) # 最后一列存sigma^2 # 初始化 beta = np.zeros(p) sigma2 = 1.0 for i in range(n_samples + burn_in): # 更新beta Sigma_beta = np.linalg.inv(X.T @ X / sigma2 + np.eye(p) / 100) mu_beta = Sigma_beta @ X.T @ y / sigma2 beta = np.random.multivariate_normal(mu_beta, Sigma_beta) # 更新sigma^2 residuals = y - X @ beta alpha_post = 1 + n / 2 beta_post = 1 + residuals @ residuals / 2 sigma2 = 1.0 / np.random.gamma(alpha_post, 1.0 / beta_post) if i >= burn_in: samples[i - burn_in, :p] = beta samples[i - burn_in, p] = sigma2 return samples # 生成模拟数据 np.random.seed(123) n = 100 X = np.column_stack([np.ones(n), np.random.normal(0, 1, n)]) y = 2.0 + 3.0 * X[:, 1] + np.random.normal(0, 1, n) samples = gibbs_sampling(X, y) print(f"beta0后验均值: {samples[:, 0].mean():.3f} (真实值: 2.0)") print(f"beta1后验均值: {samples[:, 1].mean():.3f} (真实值: 3.0)") print(f"sigma2后验均值: {samples[:, 2].mean():.3f} (真实值: 1.0)")跑出来结果应该跟真实值很接近。Gibbs采样的优势在这里体现得很明显：不需要调步长，不需要算接受率，代码写起来也直观。
4.4 Gibbs采样的局限与应对
Gibbs采样不是万能的。它最大的问题是：当参数之间存在强相关时，条件分布的移动会非常慢。想象一个狭长的椭圆后验分布，Gibbs采样每次只能沿着坐标轴方向移动，要走到椭圆另一端需要很多很多步。
应对方案有几个：
- 重参数化：把强相关的参数变换成近似独立的参数。比如中心化预测变量，可以大幅降低β₀和β₁的相关性。
- 分组更新：把强相关的参数放在一组，联合采样。如果联合条件分布可以解析采样，就用块Gibbs；如果不行，就在块内用MH。
- 退火与温度：在预热期使用高温链，让链更容易跨过低密度区域，然后逐渐降温到目标分布。
注意：Gibbs采样要求每个参数的条件分布都能解析采样。如果某个条件分布没法直接采样，就得用MH步骤替代，这就变成了Metropolis-within-Gibbs。实际建模中，纯Gibbs和混合方案都很常见。
5. 收敛诊断：怎么判断链已经“熟了”
5.1 迹图与自相关图：肉眼诊断的第一道关
迹图是最直观的诊断工具。把每次迭代的样本值画成折线图，观察几个特征：
- 链是否围绕一个稳定的均值波动，没有明显的上升或下降趋势。
- 链是否充分混合，没有长时间卡在某个区域。
- 多条链的迹图是否重叠在一起，分不清哪条是哪条。
自相关图（ACF）展示的是样本与其滞后版本的相关性。理想情况下，ACF应该快速衰减到零。如果衰减很慢，说明样本之间高度相关，有效样本量低。
def plot_diagnostics(samples, param_idx=0): fig, axes = plt.subplots(1, 2, figsize=(12, 4)) # 迹图 axes[0].plot(samples[:, param_idx], alpha=0.7, linewidth=0.5) axes[0].set_xlabel('迭代次数') axes[0].set_ylabel('参数值') axes[0].set_title('迹图') # 自相关图 from statsmodels.tsa.stattools import acf acf_vals = acf(samples[:, param_idx], nlags=100) axes[1].bar(range(len(acf_vals)), acf_vals) axes[1].set_xlabel('滞后') axes[1].set_ylabel('自相关') axes[1].set_title('自相关函数') plt.tight_layout() plt.show()5.2 Gelman-Rubin诊断：多链比较的量化指标
Gelman-Rubin诊断（也叫R-hat）是比较多条链的经典方法。思路是：如果多条链都收敛了，那么链间方差应该跟链内方差差不多。R-hat就是这两个方差的比值开根号。
计算步骤：
- 跑m条链，每条链长度n，丢弃前一半作为预热。
- 计算每条链的均值，以及所有链的总均值。
- 计算链间方差B = n * Var(链均值)。
- 计算链内方差W = 各链方差的均值。
- 估计边际后验方差：V̂ = (n-1)/n * W + B/n。
- R-hat = sqrt(V̂ / W)。
R-hat接近1说明收敛。经验阈值是R-hat < 1.1可以接受，< 1.01更好。如果R-hat大于1.1，说明链还没收敛，需要跑更长时间或者检查模型设定。
| 诊断指标 | 判断标准 | 含义 | 应对措施 |
|---|---|---|---|
| R-hat | < 1.01 | 收敛良好 | 可以开始推断 |
| R-hat | 1.01 ~ 1.1 | 勉强可接受 | 增加迭代次数 |
| R-hat | > 1.1 | 未收敛 | 检查模型、增加预热、重参数化 |
| 有效样本量 | > 1000 | 足够 | 正常推断 |
| 有效样本量 | < 100 | 不足 | 增加迭代或降低自相关 |
5.3 有效样本量：你真正有多少独立信息
MCMC样本是相关的，所以n个样本的信息量不等于n个独立样本。有效样本量（ESS）衡量的是：这些相关样本相当于多少个独立样本。
ESS的计算公式是：ESS = n / (1 + 2 * Σρₖ)，其中ρₖ是滞后k的自相关。实际计算中，通常用Geyer的初始正序列方法或者Stan的改进方法。
ESS低意味着你需要跑更长的链才能得到稳定的后验估计。我的经验是，对于后验均值估计，ESS至少要到1000；对于95%可信区间的边界估计，ESS最好到10000。
5.4 收敛诊断的常见陷阱
陷阱一：只看一条链。单链的迹图看起来平稳，不代表收敛了。链可能卡在某个局部区域，看起来稳定，但跟真实后验差很远。至少跑4条链，从不同初始点出发。
陷阱二：预热期丢太少。有些人只丢几百个样本就开始推断，结果后验估计有偏。预热期的长度应该通过诊断确定，而不是拍脑袋。
陷阱三：忽略自相关。接受率高不代表混合好。如果ACF衰减很慢，即使接受率99%，有效样本量也可能很低。
陷阱四：多峰分布。如果后验是多峰的，单条链可能只探索了一个峰。这时候需要特殊处理，比如并行回火或者更激进的建议分布。
6. 从理论到落地：MCMC实战中的经验与避坑
6.1 高维问题的维度灾难与应对
MCMC在高维问题下会变慢，这是众所周知的。维度升高时，随机游走MH的接受率会急剧下降，因为高维空间里大部分体积都在“壳”上，随机方向很容易跳到低密度区。
应对策略：
- 自适应步长：在预热期动态调整步长，让接受率维持在目标水平。Stan和PyMC都内置了这种机制。
- 利用梯度信息：Hamiltonian Monte Carlo（HMC）和No-U-Turn Sampler（NUTS）利用后验的梯度来指导提议方向，在高维问题下效率远高于随机游走MH。这是Stan的默认采样器。
- 重参数化：把参数变换到无约束空间，让后验更接近正态分布。比如方差参数取对数，相关系数用Fisher变换。
6.2 多峰分布的跨峰探索
多峰后验是MCMC的噩梦。链一旦进入一个峰，就很难跳到另一个峰，因为中间的谷底密度太低，提议点几乎总被拒绝。
解决方案：
- 并行回火：同时跑多条不同温度的链，高温链更容易跨峰，低温链精细探索。定期交换链的状态。
- 跳跃提议：设计能直接从一个峰跳到另一个峰的建议分布。这需要对问题结构有深入了解。
- 增加链数：从多个初始点出发，每个峰至少有一条链。但这样得到的样本需要加权处理。
6.3 计算资源与时间的权衡
MCMC的计算成本不低。每次迭代都要计算后验密度（先验乘似然），对于大数据集，这个计算量可能很大。
优化思路：
- 似然函数的向量化：用NumPy或JAX把循环改成矩阵运算，速度能提升几十倍。
- 子采样：大数据集下，每次迭代只用一部分数据计算似然。这引入了额外的噪声，但可以大幅加速。随机梯度MCMC就是基于这个思路。
- 并行化：多条链可以并行跑，充分利用多核CPU。Stan和PyMC都支持并行链。
6.4 常见问题速查表
| 问题现象 | 可能原因 | 排查方法 | 解决方案 |
|---|---|---|---|
| 接受率过低 | 步长太大 | 查看接受率 | 减小步长 |
| 接受率过高但混合差 | 步长太小 | 查看ACF | 增大步长 |
| R-hat > 1.1 | 未收敛 | 多链诊断 | 增加迭代、重参数化 |
| 迹图有趋势 | 预热不足 | 观察迹图前段 | 增加预热期 |
| 链卡住不动 | 多峰或强相关 | 多初始点测试 | 并行回火、分组更新 |
| ESS极低 | 自相关高 | 计算ACF | 稀释、增加迭代、改进提议 |
| 后验估计有偏 | 预热期丢弃不当 | 对比不同预热长度 | 增加预热、多链验证 |
6.5 我踩过的几个坑
第一个坑是盲目相信默认参数。早期用某个概率编程库，直接调默认采样器，结果后验估计跟理论值差很远。后来才发现默认迭代次数太少，R-hat都超过1.5了。从那以后，我养成了每次跑完先看R-hat和ESS的习惯。
第二个坑是忽略参数的相关性。有一次做层次模型，中心化参数和非中心化参数混用，导致后验出现漏斗形状，Gibbs采样直接卡死。后来改用非中心化参数化，问题迎刃而解。这个经验告诉我，参数化方式对MCMC效率的影响，有时候比算法选择还大。
第三个坑是预热期丢太多。有段时间我迷信“预热期越长越好”，结果把一半样本都丢了，有效样本量反而不够。后来才明白，预热期的长度应该通过诊断确定，而不是拍脑袋定一个固定比例。
第四个坑是在多峰问题上用单链。有一次做混合模型，后验有两个峰，单链只探索了一个，得出的结论完全错误。后来跑多条链才发现问题。多峰问题一定要多链，或者用能跨峰的采样器。
6.6 工具选型建议
如果你不想从零实现MCMC，这几个工具值得考虑：
- PyMC：Python生态里最成熟的贝叶斯建模库，支持NUTS、Metropolis、Slice等多种采样器，诊断工具完善。
- Stan：通过PyStan或CmdStanPy调用，NUTS实现非常高效，适合复杂模型。
- NumPyro/JAX：基于JAX，支持自动微分和GPU加速，适合大规模问题。
- TensorFlow Probability：适合已经用TensorFlow的团队，支持HMC和变分推断。
我的建议是：学习和理解原理时，自己手写MH和Gibbs；实际项目里，优先用Stan或PyMC，把精力放在模型设计上，而不是采样器调优上。
7. 一个完整的实战案例：贝叶斯逻辑回归的MCMC实现
7.1 问题设定与数据生成
用一个二分类问题来串起整个流程。假设有n个观测，每个观测有d个特征，响应变量y是0或1。模型是：
P(yᵢ=1) = sigmoid(Xᵢβ)
先验设β ~ N(0, 10²I)。后验分布没有解析形式，需要用MCMC采样。
import numpy as np def sigmoid(z): return 1 / (1 + np.exp(-z)) def log_posterior(beta, X, y): """对数后验（未归一化）""" log_prior = -0.5 * np.sum(beta**2) / 100 log_lik = np.sum(y * np.log(sigmoid(X @ beta) + 1e-10) + (1 - y) * np.log(1 - sigmoid(X @ beta) + 1e-10)) return log_prior + log_lik def mh_logistic(X, y, n_samples=50000, burn_in=10000, step_size=0.1): n, d = X.shape samples = np.zeros((n_samples, d)) beta = np.zeros(d) n_accepted = 0 for i in range(n_samples + burn_in): beta_proposal = beta + np.random.normal(0, step_size, d) log_alpha = log_posterior(beta_proposal, X, y) - log_posterior(beta, X, y) if np.log(np.random.uniform()) < log_alpha: beta = beta_proposal if i >= burn_in: n_accepted += 1 if i >= burn_in: samples[i - burn_in] = beta accept_rate = n_accepted / n_samples return samples, accept_rate # 生成模拟数据 np.random.seed(456) n, d = 500, 3 X = np.column_stack([np.ones(n), np.random.normal(0, 1, (n, d-1))]) true_beta = np.array([0.5, -1.0, 2.0]) y = np.random.binomial(1, sigmoid(X @ true_beta)) samples, acc_rate = mh_logistic(X, y) print(f"接受率: {acc_rate:.3f}") print(f"后验均值: {samples.mean(axis=0)}") print(f"真实参数: {true_beta}")7.2 结果解读与诊断
跑完这个例子，你应该能看到后验均值跟真实参数比较接近。接受率大概在0.2到0.3之间，对于三维问题稍微偏低，可以适当减小步长。
接下来做诊断：
# 计算R-hat（需要跑多条链） def compute_r_hat(chains): """chains: list of arrays, each shape (n_samples, d)""" m = len(chains) n = chains[0].shape[0] chain_means = np.array([c.mean(axis=0) for c in chains]) overall_mean = chain_means.mean(axis=0) B = n * np.var(chain_means, axis=0, ddof=1) W = np.mean([np.var(c, axis=0, ddof=1) for c in chains], axis=0) V_hat = (n - 1) / n * W + B / n r_hat = np.sqrt(V_hat / W) return r_hat # 跑4条链 chains = [] for seed in [1, 2, 3, 4]: np.random.seed(seed) chain, _ = mh_logistic(X, y, n_samples=20000, burn_in=5000) chains.append(chain) r_hat = compute_r_hat(chains) print(f"R-hat: {r_hat}")如果R-hat都在1.01以下，说明收敛良好。如果某个参数的R-hat偏高，可能需要增加迭代次数或者调整步长。
7.3 后验预测与决策
有了后验样本，就可以做后验预测了。对于新的输入x_new，预测概率的分布是：
P(y_new=1 | x_new, 数据) = ∫ sigmoid(x_new β) p(β | 数据) dβ
用MCMC样本近似这个积分：
def posterior_predictive(samples, X_new): """对每个后验样本计算预测概率，然后平均""" probs = sigmoid(X_new @ samples.T) # shape: (n_new, n_samples) return probs.mean(axis=1), probs.std(axis=1) X_new = np.array([[1.0, 0.5, -0.5], [1.0, -1.0, 1.0]]) mean_prob, std_prob = posterior_predictive(samples, X_new) print(f"预测概率均值: {mean_prob}") print(f"预测概率标准差: {std_prob}")这个预测分布不仅给了点估计，还给了不确定性量化。这是贝叶斯方法相比频率学派的一大优势。
8. 从MH到现代采样器：MCMC的演进脉络
8.1 Hamiltonian Monte Carlo的核心直觉
HMC把参数空间想象成一个物理系统，参数是位置，后验密度的负对数是势能。系统还有一个动量变量，动量从标准正态分布中抽取。然后模拟哈密顿动力学：势能转化为动能，动能转化为势能，粒子在势能面上滚动。
HMC的提议方向由梯度决定，所以它总是朝着高密度区域走，接受率远高于随机游走MH。在几百维的问题上，HMC的效率可以比MH高几个数量级。
代价是需要计算后验的梯度。对于复杂模型，梯度推导很麻烦，但自动微分工具（如JAX、PyTorch、Stan）解决了这个问题。
8.2 NUTS：不用手动调轨迹长度
HMC需要指定轨迹长度（模拟多少步动力学）。太短，探索不充分；太长，浪费计算。NUTS（No-U-Turn Sampler）自动决定轨迹长度：当轨迹开始“掉头”往回走时，就停止。
NUTS是Stan的默认采样器，也是目前最先进的MCMC算法之一。它几乎不需要手动调参，对复杂后验的适应能力很强。
8.3 变分推断与MCMC的取舍
变分推断（VI）用优化代替采样，速度快但精度低。MCMC精度高但速度慢。实际项目中怎么选？
- 如果只需要点估计和粗略的不确定性，VI够用。
- 如果需要精确的后验推断，特别是可信区间和模型比较，MCMC更可靠。
- 大数据场景下，随机梯度VI（如ADVI）可以扩展到海量数据，MCMC则受限于每次迭代的全数据似然计算。
我个人的做法是：先用VI快速探索，确定模型大致合理后，再用MCMC做最终推断。这样既节省时间，又保证精度。
9. 写在最后：一些个人体会
MCMC是贝叶斯统计里最迷人的部分之一。它的核心思想——用马尔可夫链的平稳分布来逼近目标分布——既优雅又实用。但真正用好MCMC，光懂原理不够，还得有调参的经验和诊断的直觉。
我刚开始学MCMC的时候，总觉得接受率越高越好，后来才发现混合效率才是关键。再后来，又发现参数化方式比采样器选择更重要。这些经验都是踩坑踩出来的，书本上不会写。
如果你正在学MCMC，我的建议是：先手写一个MH采样器，跑通一个简单的一维例子，把迹图和ACF画出来，感受一下步长对混合的影响。然后手写一个Gibbs采样器，跑一个贝叶斯线性回归，体会条件分布采样的便利。最后，用PyMC或Stan跑一个稍微复杂的模型，对比手写实现和成熟工具的结果。
这个过程中，你会遇到各种奇怪的现象：链卡住不动、R-hat爆表、ESS低得可怜。别慌，这些都是正常的。每一个问题背后都有原因，找到原因的过程就是进步的过程。
最后分享一个小技巧：跑MCMC的时候，把随机种子固定下来，这样结果可复现。但诊断的时候，一定要换不同的种子跑多条链，看看结果是否稳定。单链的结果再漂亮，也不如多链的一致性好。
