诊断MCMC收敛性和性能

MCMC采样器，例如Metropolis-Hastings、Gibbs，甚至是强大的HMC/NUTS，会从目标后验分布 $p(\theta | D)$ 中生成成千上万、甚至数百万个样本。此时，一个重要问题浮现：如何验证这些样本的质量？又如何确信它们确实代表了所要近似的后验分布？MCMC诊断正是为了解决这些问题。

运行MCMC采样器会产生一系列参数值，称为“链”。核心思想是，在某些理论条件（遍历性）下，如果这条链运行足够长时间，它最终会产生看起来像是从目标平稳分布（我们期望的后验）中抽取的样本。诊断帮助我们评估两个重要方面：

1. 收敛性： 链是否已经“忘记”了它的起始点，并已稳定下来以访问平稳分布的高概率区域？
2. 效率/性能： 一旦链收敛，它访问这个平稳分布的效率如何？样本是否高度相关，还是它们提供了关于后验分布的相对独立信息？

实践中，仅依赖MCMC的理论保证是不够的。有限的运行次数、复杂的后验分布以及采样器调优问题意味着我们必须检查输出。下面我们考察用于此检查的标准工具。

可视化诊断：查看你的链

通常，第一步就是简单地查看输出。可视化检查可以快速显示明显的问题。

轨迹图

轨迹图绘制参数的采样值随迭代次数的变化。对于每个参数（或选定的一些重要参数），你都应该生成一个轨迹图。

• 要看什么（好迹象）： 一个健康的轨迹图应该看起来像平稳的“白噪声”——围绕一个稳定的平均值快速波动，没有长期趋势或漂移。它通常像一只“毛毛虫”在图中蜿蜒爬行。这表明链正在参数空间的一个一致区域中访问。如果你从不同的起始点运行多条链，它们的轨迹图应该重叠并混合在一起。
• 要留意什么（坏迹象）：
  • 初始漂移/趋势： 图中显示有明显的上升或下降趋势，尤其是在早期。这表明链尚未收敛，仍在远离其起始点（“预热”或“燃烧期”）。
  • 混合缓慢： 图看起来像一条平滑、缓慢蜿蜒的河流，而不是模糊的噪声。这意味着连续的样本高度相关，链访问空间效率很低。
  • 停滞： 链在某个特定值或区域停滞多代迭代，然后才跳到其他地方。这可能发生在多峰分布中或采样器调优不佳时。
  • 缺乏重叠（多条链）： 如果来自不同链的轨迹图访问完全不同的区域并且不混合，这强烈表明未收敛，或者链卡在了后验分布的不同众数中。

轨迹图示例。参数A在两条链中都显示出明显的上升趋势，表明它们尚未收敛。参数B显示出平稳行为，两条链都访问了相同的范围，表明该参数已收敛。

自相关图 (ACF)

轨迹图暗示了连续样本之间的相关性。自相关图直接量化了这一点。ACF图显示了链中不同“滞后”（样本之间相隔多少步）下样本之间的相关性。对于参数 $\theta$，滞后 $k$ 自相关是 $\theta_t$ 和 $\theta_{t+k}$ 之间的相关性。

• 要看什么（好迹象）： 随着滞后增加，自相关应迅速降至零。这表明相隔步数较少的样本几乎是独立的，意味着链混合良好且高效地访问参数空间。
• 要留意什么（坏迹象）： 高自相关性持续存在很多滞后。这表明混合缓慢。链正在采取非常小的步长，相邻样本高度冗余。这降低了你拥有的独立样本的有效数量。

ACF图示例。绿线显示自相关性快速下降，表明混合良好。红线显示自相关性缓慢衰减，表明混合不良以及相邻样本之间的高度相关性。

定量诊断：数字会说明问题

虽然可视化检查非常宝贵，但定量指标提供了更客观的评估，特别是对于自动化检查或具有许多参数的复杂模型。

潜在尺度收敛因子 ($\hat{R}$, R-hat)

被称为“R-hat”，这可以说是最受欢迎的收敛诊断方法。由Gelman和Rubin提出，它需要运行多条链 ($M \ge 2$)，这些链从参数空间中分散的起始点初始化。

核心思想是比较每条链内部的方差与链之间的方差。

1. 运行 $M$ 条链，每条链长度为 $N$（丢弃初始预热样本后）。
2. 对于特定参数 $\theta$：
   • 计算每条链 $j$ 内部的均值 $\bar{\theta}_{.j}$ 和方差 $s_j^2$。
   • 计算所有链中所有样本的总体均值 $\bar{\theta}_{..}$。
   • 计算链均值之间的方差： $$B = \frac{N}{M-1} \sum_{j=1}^M (\bar{\theta}_{.j} - \bar{\theta}_{..})^2$$
   • 计算链内部的平均方差： $$W = \frac{1}{M} \sum_{j=1}^M s_j^2$$
3. 将 $\theta$ 的边际后验方差估计为 $W$ 和 $B$ 的加权平均值： $$\hat{\text{Var}}^+(\theta | D) = \frac{N-1}{N} W + \frac{1}{N} B$$
4. 计算潜在尺度收敛因子： $$\hat{R} = \sqrt{\frac{\hat{\text{Var}}^+(\theta | D)}{W}}$$

解释：

• 如果链已收敛到目标分布，它们都应该访问相同的区域。在这种情况下，链内方差 ($W$) 和按比例缩放的链间方差 ($B/N$) 应该相似。估计的后验方差 $\hat{\text{Var}}^+$ 将接近 $W$，使得 $\hat{R}$ 接近1.0。
• 如果链未收敛，起始点仍会产生影响。链均值之间的方差 ($B$) 将明显大于链内平均方差 ($W$)，导致 $\hat{R} > 1.0$。

经验法则： $\hat{R}$ 值低于1.01（或有时放宽到1.05或1.1）通常被认为表明该参数收敛可接受。然而，检查所有你关心的参数的 $\hat{R}$ 很重要。任何参数的 $\hat{R}$ 值高都表明存在问题。请记住，$\hat{R} \approx 1$ 是收敛的必要但不充分条件。

有效样本量 (ESS)

MCMC样本通常是自相关的。长度为 $N$ 的链不包含 $N$ 个关于后验的独立信息。有效样本量 (ESS) 估计链所代表的独立样本的等效数量。

高自相关（混合缓慢）导致低ESS，而低自相关（混合良好）导致ESS接近实际样本数量 $N$。

对于基于 $M$ 条预热后长度为 $N$ 的链的参数 $\theta$ 的ESS，可以利用自相关结构估计，或从计算 $\hat{R}$ 时使用的量中推导出来：

$$\text{ESS} = \frac{MN}{1 + 2 \sum_{k=1}^{\infty} \rho_k} \approx \frac{MN}{\hat{R}} \quad (\text{使用涉及 } \hat{R} \text{ 的简化关系})$$

（注意：软件程序包中通常使用涉及自相关性的更精确计算）。

解释：

• 低ESS（例如，小于100或小于几百，取决于应用所需的精度）表明，尽管进行了多次MCMC迭代，但你关于该参数后验分布的独立信息相对较少。你的后验汇总（如均值或可信区间）可能不可靠。
• 高ESS表明链混合良好，你的后验估计很可能基于足够的信息。

经验法则： 对于每个你关心的参数，争取ESS至少达到几百，以确保后验均值和分位数估计的稳定。如果ESS较低，你需要运行更长时间的链或提高采样器的效率（例如，通过模型重参数化或使用像NUTS这样更好的采样器，如果你尚未这样做）。

接受率

对于Metropolis-Hastings、HMC和NUTS等采样器，接受率（提议步数中被接受的比例）是一个有用的诊断方法，特别是用于调优。

• 接受率非常低： 采样器提议的移动经常被拒绝，通常是因为提议的步长过大，落在了低概率区域。链会卡住并且混合不良。
• 接受率非常高： 采样器提议的移动几乎总是被接受。这通常意味着提议的步长太小。链虽然移动，但在空间中移动非常缓慢且效率低下（高自相关）。

最佳范围： 理论结果表明，最佳接受率取决于采样器和问题的维度。

• 随机游走Metropolis： 1维目标约为0.44，多维目标约为0.234。
• HMC/NUTS： 目标接受率通常更高，通常在0.65到0.8之间，甚至更高，具体取决于具体的实现和自适应策略。请查阅你的概率编程语言（如PyMC、Stan、Turing.jl等）的文档。

偏离这些范围表明采样器的调优参数（例如HMC/NUTS中的步长或MH中的提议方差）可能需要调整。像NUTS这样的现代采样器通常会自动调整这些参数，但检查接受率仍然是一个很好的检查方法。

常见问题及处理方法

诊断通常会显示潜在问题：

• 未收敛 ($\hat{R} \gg 1$，轨迹图不重叠)：
  • 让链运行更长时间（增加 $N$）。
  • 丢弃更多初始样本（增加预热/燃烧期）。
  • 检查模型规范：模型是否可识别？先验是否合理？
  • 尝试不同、更分散的起始值。
• 混合缓慢 (高ACF，低ESS)：
  • 让链运行更长时间。
  • 模型重参数化：有时改变参数的定义方式可以显著改善采样器的几何特性（例如，非中心化参数化）。
  • 使用更高效的采样器：对于复杂、连续的后验分布，HMC/NUTS通常比随机游走Metropolis好得多。
  • 调整采样器调优参数（如果适用且未自动化）。
• 多峰性（链卡在不同位置，如果链没有访问其他众数，$\hat{R}$ 可能较低）：
  • 这具有挑战性。需要非常长的运行时间、专用采样器（例如，并行退火）或仔细的模型检查/重新思考。如果所有链都卡在同一个局部众数中，诊断可能会错误地表明收敛。

汇总：最佳实践

诊断MCMC输出对于可靠的贝叶斯推断来说很重要。

1. 始终运行多条链（通常 $M \ge 4$）。 这是计算 $\hat{R}$ 所必需的，并有助于可视化评估收敛性。
2. 从分散的起始点初始化链。 使用先验或其他方法生成不同的起始位置。
3. 丢弃初始预热/燃烧期样本。 这些不代表平稳分布。概率编程语言通常默认丢弃总运行量的一半。
4. 使用多种诊断方法： 综合使用轨迹图、$\hat{R}$ 和ESS来检查所有重要参数。不要只信任一个指标。检查采样器特有的诊断方法，如接受率或树深度（对于NUTS）。
5. 对于你关心的参数，目标是 $\hat{R} < 1.01$ 且ESS > 几百。
6. 运用概率编程语言功能： PyMC、Stan、CmdStanPy、ArviZ、Turing.jl等库提供了自动计算和汇总这些诊断的函数。
7. 请记住： 诊断检查的是未收敛。它们不能证明收敛。通过所有诊断可以增加信心，但保持警惕（并思考你的模型和后验分布的几何特性）仍然很重要。

通过仔细应用这些诊断工具，你可以获得信心，相信你的MCMC模拟提供了后验分布的忠实近似，从而实现可靠的贝叶斯推断。