4. MCMC收敛诊断：迹图、自相关图、Gelman-Rubin诊断、有效样本量、收敛前的老化期设置
MCMC跑完了，拿到一堆样本，然后呢？
直接拿去用？别急。我见过太多人把没收敛的链子直接拿去算参数，结果模型一塌糊涂。说白了，MCMC的核心问题就一个：你拿到的样本，到底能不能代表真实的后验分布？
这一章，我们就来聊聊怎么判断链子是否收敛。我会把常用的诊断工具过一遍，顺便讲讲我踩过的坑。
4.1 迹图（Trace Plot）——最直观的第一眼
迹图，就是把采样值按迭代次数画出来。横轴是迭代步数，纵轴是参数值。
我个人习惯，拿到MCMC结果第一件事就是画迹图。为什么？因为一眼就能看出问题。
好的迹图长什么样？
- 像一条“毛虫”，在某个值附近随机波动
- 没有明显的趋势或周期性
- 不同初始值的链子最终混在一起
坏的迹图长什么样？
- 有明显的上升或下降趋势——说明还没收敛
- 长时间停留在某个值——可能卡在局部模式了
- 链子之间分得很开——初始值影响太大
我曾经做过一个波动率模型，跑了5000步，迹图看起来挺平稳。结果我放大一看，前1000步其实还在爬坡。嗯，这里要注意：肉眼判断有时会骗人，尤其是当纵轴范围很大时。
import pymc3 as pm
import matplotlib.pyplot as plt
# 假设 trace 是采样结果
pm.plot_trace(trace, var_names=['mu', 'sigma'])
plt.show()4.2 自相关图（Autocorrelation Plot）——样本到底有多“独立”？
MCMC的样本不是独立的。相邻的样本之间肯定有相关性。自相关图就是看这个相关性有多强。
你想想看，如果自相关降得很慢，说明链子走得很慢，有效信息少。说白了，你采样了10000个点，可能真正有用的只有几百个。
怎么看自相关图？
- 理想情况：自相关在滞后几阶后迅速降到0附近
- 糟糕情况：自相关在滞后几十阶后仍然很高
- 如果自相关长期不衰减，建议增加 thinning（间隔采样）
我记得有一次做随机波动率模型，自相关图显示滞后50阶还有0.3的相关性。我当时就意识到，这链子效率太低了。后来我调整了提议分布，自相关才降下来。
pm.plot_autocorr(trace, var_names=['mu'], max_lag=50)
plt.show()4.3 Gelman-Rubin诊断（R-hat）——多条链子说了算
单条链子可能骗人。但多条链子同时骗人的概率就小多了。
Gelman-Rubin诊断的核心思想：比较链间方差和链内方差。如果链子都收敛到了同一个分布，那么链间方差应该和链内方差差不多。
R-hat值接近1，说明收敛。一般经验是R-hat < 1.1就算可以。但我个人比较严格，我通常要求R-hat < 1.01。
| R-hat值 | 含义 | 建议操作 |
|---|---|---|
| 1.0 - 1.01 | 良好收敛 | 可以放心使用 |
| 1.01 - 1.05 | 可能收敛 | 建议增加迭代次数 |
| 1.05 - 1.1 | 需要警惕 | 检查模型设定或增加链子 |
| > 1.1 | 未收敛 | 重新跑，或者改模型 |
注意：R-hat对多峰分布不敏感。如果后验分布有两个峰，每条链子各占一个峰，R-hat可能还是接近1。这时候要结合迹图一起看。
# PyMC3 自动计算 R-hat
summary = pm.summary(trace, var_names=['mu'])
print(summary[['mean', 'sd', 'r_hat']])4.4 有效样本量（ESS）——你的样本到底值多少？
有效样本量，英文叫Effective Sample Size。它告诉你：你手里的N个样本，相当于多少个独立样本？
如果自相关很高，ESS会远小于实际样本量。比如你采样了10000个点，ESS可能只有500。说白了，你花了10000次计算，只得到了500个有效信息。
ESS的经验法则：
- 对于均值估计：ESS > 400 基本够用
- 对于分位数估计（如5%分位）：ESS > 1000 才靠谱
- 对于尾部概率：ESS > 5000 都不一定够
我在做VaR计算时吃过这个亏。当时ESS只有200多，算出来的95%分位数波动很大。后来我把链子跑长到50000步，ESS到了2000左右，结果才稳定下来。
# 计算有效样本量
ess = pm.effective_n(trace, var_names=['mu'])
print(f"有效样本量: {ess['mu']}")4.5 收敛前的老化期设置（Burn-in）——扔掉“热身”部分
MCMC刚开始时，链子还没走到高概率区域。这部分样本如果保留，会污染整个估计。
老化期，就是扔掉前面一部分样本。问题是：扔多少？
我个人习惯的做法：
- 先跑一个短链（比如2000步），看迹图判断大概什么时候稳定
- 把老化期设为稳定点之前的1.5倍，留点余量
- 正式跑长链，然后重新检查收敛
一个实用的技巧：如果你不确定老化期该设多少，可以多设一些。扔掉1000步和扔掉2000步，如果链子已经收敛，结果应该差不多。如果差别很大，说明链子还没收敛。
我曾经犯过一个错误：老化期设得太短，导致参数估计有偏。后来我养成了一个习惯——每次跑完MCMC，先画迹图，再决定老化期。不要偷懒用默认值。
# 手动设置老化期
burn_in = 1000
trace_after_burn = trace[burn_in:]4.6 知识体系总览
下面这张图，把本章的核心逻辑串起来了。你可以把它当作一个检查清单：
这五个工具，我建议你每次跑MCMC时都用一遍。不要只看一个。它们各有优缺点，组合起来才能全面判断。
我的个人检查清单：
- 画迹图——看趋势
- 画自相关图——看效率
- 算R-hat——多条链子交叉验证
- 算ESS——看有效信息量
- 设老化期——扔掉热身部分
做到这五步，你的MCMC结果基本就靠谱了。剩下的，就是模型本身的问题了——那是另一门课的内容。
公众号:蓝海资料掘金营，微信deep3321