MCMC收敛诊断

1 引言

MCMC（马尔可夫链蒙特卡洛）方法的核心假设是链在足够长时间后收敛到目标分布。收敛诊断是验证这一假设是否成立的关键技术。¹

本章系统介绍MCMC收敛诊断的理论与实践，包括：

Gelman-Rubin统计量（R̂）
有效样本量（ESS）
Geweke检验
谱密度分析
实用诊断流程

2 为什么需要收敛诊断

2.1 MCMC的双重极限

MCMC样本来自两个渐近过程：

链内： $N \to \infty$ （样本量增加）
链间： $M \to \infty$ （独立链数增加）

问题：实践中只能运行有限长度链。

2.2 收敛失败的类型

类型	描述	诊断方法
未充分预热	初始状态远离后验	轨迹图、预热期分析
局部模式陷阱	链困在局部后验模式	多链比较
慢混合	链移动缓慢	自相关分析、ESS
周期性	链有规律振荡	自相关函数

2.3 诊断原则

关键洞察：没有任何单一诊断能完全保证收敛。

推荐策略：

使用多种诊断方法
结合可视化分析
保守的参数选择
充足的预热和采样

3 Gelman-Rubin统计量（R̂）

3.1 基本思想

R̂通过比较多条独立链的链内方差和链间方差来检测收敛。²

核心假设：

如果链已收敛，所有链来自同一分布
链间方差应与链内方差相当

3.2 正式定义

设运行 $m$ 条独立链，每条链有 $n$ 个样本（排除预热期）。

链内方差（Within-chain variance）：

W = \frac{1}{m} j = 1 \sum m s_{j}^{2}

其中 $s_{j}^{2} = \frac{1}{n - 1} \sum_{i = 1}^{n} (θ_{ij} - \overset{ˉ}{θ}_{j \cdot})^{2}$ 是第 $j$ 链的样本方差。

链间方差（Between-chain variance）：

B = \frac{n}{m - 1} j = 1 \sum m (\overset{ˉ}{θ}_{j \cdot} - \overset{ˉ}{θ}_{\cdot\cdot})^{2}

其中 $\overset{ˉ}{θ}_{j \cdot}$ 是第 $j$ 链的均值， $\overset{ˉ}{θ}_{\cdot\cdot}$ 是总均值。

总体方差估计（后验方差的无偏估计）：

\hat{Var} (θ ∣ y) = \frac{n - 1}{n} W + \frac{1}{n} B

3.3 R̂计算

潜在尺度缩减因子（PSRF）：

\hat{R} = \frac{Var ^ ( θ ∣ y )}{W} = \frac{\frac{n - 1}{n} W + \frac{1}{n} B}{W}

直观理解：

$\hat{R} \approx 1$ ：链间方差与链内方差相当 → 可能收敛
$\hat{R} > 1.1$ ：链间方差显著大于链内方差 → 未充分混合
$\hat{R} > 1.2$ ：严重未收敛问题

3.4 多变量R̂

对于多参数问题，分别计算每个参数的R̂可能不够。

修正R̂：

\hat{R} = \frac{V ^}{W}

其中 $\hat{V}$ 是参数向量 $θ = (θ_{1}, \dots, θ_{p})$ 的协方差矩阵估计。

3.5 代码实现

import numpy as np
 
def gelman_rubin(chains):
    """
    计算Gelman-Rubin统计量
    
    Parameters:
    chains: shape (m, n) 的数组，m条链，每链n个样本
    
    Returns:
    R_hat: 每个参数的R̂值
    """
    m, n = chains.shape[0], chains.shape[1]
    
    # 链内方差
    chain_means = np.mean(chains, axis=1)
    chain_vars = np.var(chains, axis=1, ddof=1)  # 使用n-1
    W = np.mean(chain_vars, axis=0)
    
    # 链间方差
    grand_mean = np.mean(chain_means, axis=0)
    B = n / (m - 1) * np.sum((chain_means - grand_mean)**2, axis=0)
    
    # 总体方差估计
    var_theta = (n - 1) / n * W + 1 / n * B
    
    # R̂
    R_hat = np.sqrt(var_theta / W)
    
    return R_hat
 
def multi_R_hat(chains):
    """多变量R̂"""
    m, n = chains.shape[0], chains.shape[1]
    p = chains.shape[2]  # 参数维度
    
    # 展平为 (m*n, p)
    flat_chains = chains.reshape(-1, p)
    
    # 计算协方差矩阵
    # ...
    
    return R_hat_multivariate

3.6 阈值选择

R̂范围	解释	建议
< 1.01	理想收敛	继续采样
1.01 - 1.05	可接受	继续或增加迭代
1.05 - 1.1	边缘	增加迭代、检查混合
> 1.1	问题	增加预热、检查链行为
> 1.2	严重问题	重新设计模型/采样器

4 有效样本量（ESS）

4.1 核心概念

有效样本量（Effective Sample Size, ESS）衡量MCMC样本的信息量。³

问题：自相关样本的信息量低于独立样本。

定义：与当前链等效的独立样本数量。

ESS = \frac{N}{1 + 2 \sum _{k = 1}^{\infty} ρ _{k}}

其中 $ρ_{k}$ 是滞后 $k$ 的自相关函数。

4.2 自相关函数

样本自相关函数（ACF）：

ρ_{k} = \frac{\sum _{i = 1}^{N - k} ( θ _{i} - θ ˉ ) ( θ _{i + k} - θ ˉ )}{\sum _{i = 1}^{N} ( θ _{i} - θ ˉ ) ^{2}}

4.3 ESS估计

批量均值法（Batch Means）：

将链划分为 $b$ 个批量，每批 $a = N / b$ 个样本：

ESS = \frac{N}{τ ^}

其中 $\overset{τ}{^} = a \cdot \frac{Var ( θ ˉ _{b} )}{Var ( θ )} / a$ 是估计的时间窗口。

def effective_sample_size(samples, max_lag=None):
    """
    计算有效样本量
    """
    n = len(samples)
    mean = np.mean(samples)
    var = np.var(samples, ddof=1)
    
    if max_lag is None:
        max_lag = n // 2
    
    # 计算自相关
    acf = []
    for k in range(max_lag + 1):
        if k == 0:
            acf.append(1.0)
        else:
            cov = np.mean((samples[:-k] - mean) * (samples[k:] - mean))
            acf.append(cov / var)
    
    # 使用直到ACF变为负的滞后
    negative_idx = next((i for i, x in enumerate(acf) if x < 0), len(acf))
    acf_truncated = acf[:negative_idx]
    
    # ESS = N / (1 + 2 * sum(acf))
    ess = n / (1 + 2 * sum(acf_truncated[1:]))
    
    return int(ess)

4.4 尾部ESS

标准ESS可能低估复杂后验的信息量。尾部ESS关注分布尾部的采样质量：

def tail_ess(samples, alpha=0.05):
    """
    计算尾部有效样本量
    alpha: 尾部概率（如0.05表示5%和95%分位数）
    """
    n = len(samples)
    
    # 批量均值
    batch_size = int(np.sqrt(n))
    n_batches = n // batch_size
    
    # 计算每个分位数的ESS
    lower_q = np.percentile(samples, 100 * alpha)
    upper_q = np.percentile(samples, 100 * (1 - alpha))
    
    # ...

4.5 ESS阈值

ESS值	解释
< 100	不充足
100 - 500	可接受
500 - 1000	良好
> 1000	非常良好

4.6 效率指标

ESS/n（每迭代的ESS）：

Efficiency = \frac{ESS}{N}

高效率意味着链快速混合。

5 Geweke检验

5.1 基本思想

Geweke检验比较链的两个不同部分来判断是否来自同一分布。⁴

假设：

$H_{0}$ ：两部分样本来自同一分布
$H_{1}$ ：两部分样本来自不同分布

5.2 检验统计量

选择链的开始部分 $A$ 和结尾部分 $B$ ：

z = \frac{θ ˉ _{A} - θ ˉ _{B}}{\frac{S _{A}}{n _{A}} + \frac{S _{B}}{n _{B}}}

其中 $S_{A}, S_{B}$ 是使用谱密度估计的方差。

渐近分布：在 $H_{0}$ 下， $z \to N (0, 1)$ 。

5.3 代码实现

from scipy import stats
 
def geweke_test(samples, first_frac=0.1, last_frac=0.5, alpha=0.05):
    """
    Geweke收敛诊断检验
    
    Parameters:
    samples: MCMC样本
    first_frac: 开头部分的比例
    last_frac: 结尾部分的比例
    alpha: 显著性水平
    
    Returns:
    z_score: 检验统计量
    p_value: p值
    converged: 是否收敛
    """
    n = len(samples)
    
    # 选择两部分
    nA = int(n * first_frac)
    nB = int(n * last_frac)
    
    # 避免重叠
    nB = min(nB, n - nA)
    
    samples_A = samples[:nA]
    samples_B = samples[-nB:]
    
    # 计算均值和谱密度方差
    mean_A = np.mean(samples_A)
    mean_B = np.mean(samples_B)
    
    # 使用Bartlett窗估计谱密度
    var_A = spectral_density_variance(samples_A)
    var_B = spectral_density_variance(samples_B)
    
    # 检验统计量
    z = (mean_A - mean_B) / np.sqrt(var_A / nA + var_B / nB)
    
    # p值（双侧检验）
    p_value = 2 * (1 - stats.norm.cdf(abs(z)))
    
    # 收敛判定
    converged = p_value > alpha
    
    return z, p_value, converged
 
def spectral_density_variance(samples):
    """估计谱密度方差（用于Geweke检验）"""
    n = len(samples)
    
    # 计算自相关
    acf = []
    for k in range(n // 2):
        acf.append(np.mean((samples[:-k] - np.mean(samples)) * 
                          (samples[k:] - np.mean(samples))))
    
    acf = np.array(acf)
    
    # 截断自相关（使用保守估计）
    M = int(np.sqrt(n))
    acf_truncated = acf[:M]
    
    # 谱密度方差
    return np.var(samples) + 2 * np.sum(acf_truncated)

5.4 解读

结果	解释
p > 0.05	不能拒绝 $H_{0}$ → 可能收敛
p < 0.05	拒绝 $H_{0}$ → 显著差异，未收敛

6 谱密度分析

6.1 蒙特卡洛标准误

蒙特卡洛采样的标准误（Standard Error）：

SE (\hat{I}) = \frac{σ}{N} \cdot \frac{1 + 2 \sum _{k = 1}^{\infty} ρ _{k}}{1} = \frac{σ}{ESS}

关键洞察：高自相关增大标准误，减小ESS。

6.2 功率谱密度

功率谱密度（PSD）描述不同频率成分的功率分布：

S (ω) = k = - \infty \sum \infty ρ_{k} e^{- iωk}

特征频率：

低频： $ω \approx 0$ ，对应慢混合
高频： $ω \approx π$ ，对应噪声

6.3 诊断应用

诊断慢混合：

功率集中在低频 → 链移动缓慢
均匀分布 → 快速混合

import matplotlib.pyplot as plt
from scipy import signal
 
def diagnose_mixing(samples):
    """诊断链的混合性质"""
    n = len(samples)
    
    # 计算功率谱密度
    frequencies, psd = signal.periodogram(samples - np.mean(samples), 
                                          fs=2*np.pi, 
                                          window='hann')
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 1. 轨迹图
    axes[0, 0].plot(samples)
    axes[0, 0].set_title('Trace Plot')
    axes[0, 0].set_xlabel('Iteration')
    axes[0, 0].set_ylabel('Value')
    
    # 2. 自相关函数
    acf = np.correlate(samples - np.mean(samples), 
                        samples - np.mean(samples), mode='full')
    acf = acf[len(acf)//2:]
    acf /= acf[0]
    axes[0, 1].plot(acf[:50])
    axes[0, 1].set_title('Autocorrelation Function')
    axes[0, 1].set_xlabel('Lag')
    axes[0, 1].set_ylabel('ACF')
    
    # 3. 功率谱密度
    axes[1, 0].semilogy(frequencies[1:], psd[1:])
    axes[1, 0].set_title('Power Spectral Density')
    axes[1, 0].set_xlabel('Frequency')
    axes[1, 0].set_ylabel('Power')
    
    # 4. 直方图
    axes[1, 1].hist(samples, bins=50, density=True, alpha=0.7)
    axes[1, 1].set_title('Histogram')
    axes[1, 1].set_xlabel('Value')
    axes[1, 1].set_ylabel('Density')
    
    plt.tight_layout()
    plt.savefig('mcmc_diagnostics.png')
    
    return {'frequencies': frequencies, 'psd': psd}

7 实用诊断流程

7.1 推荐流程

开始采样
   │
   ├─ 阶段1：预热期诊断
   │    │
   │    ├─ 可视化轨迹图（检查趋势）
   │    ├─ 调整参数（步长、迭代数）
   │    └─ 重复直到稳定
   │
   ├─ 阶段2：收敛诊断
   │    │
   │    ├─ 运行多条独立链（≥4）
   │    ├─ 计算R̂（< 1.1）
   │    ├─ 计算ESS（足够大）
   │    └─ Geweke检验（p > 0.05）
   │
   ├─ 阶段3：质量评估
   │    │
   │    ├─ 谱密度分析（无异常模式）
   │    ├─ 蒙特卡洛标准误（足够小）
   │    └─ 后验汇总（合理范围）
   │
   └─ 结束

7.2 诊断检查清单

检查项	方法	阈值
R̂	gelman_rubin	< 1.1
ESS	effective_sample_size	> 100/参数
ESS效率	ESS/N	> 1%
Geweke	geweke_test	p > 0.05
自相关	acf图	在合理滞后内衰减
轨迹图	可视化	无趋势/周期性

7.3 常见问题与解决方案

问题	诊断信号	解决方案
未充分预热	轨迹有趋势	增加预热期
慢混合	高自相关、低ESS	减小步长、增大学习率
多峰分布	多链在不同模式	增加迭代、使用重参数化
局部模式	单链在单一模式	尝试不同初始化
参数相关	高自相关	重参数化（正交化）

7.4 完整诊断工具

class MCMCDiagnostics:
    """完整的MCMC诊断工具"""
    
    def __init__(self, chains, warmup=None):
        """
        Parameters:
        chains: shape (m, n) 或 (m, n, p) 的数组
        warmup: 预热期长度
        """
        self.chains = chains
        self.warmup = warmup
        self.n_chains = chains.shape[0]
        self.n_samples = chains.shape[1]
        
    def run_diagnostics(self):
        """运行所有诊断"""
        results = {}
        
        # 1. R̂
        results['R_hat'] = gelman_rubin(self.chains)
        
        # 2. ESS
        results['ESS'] = self.compute_ess_all_params()
        
        # 3. ESS效率
        total_samples = self.n_chains * self.n_samples
        results['ESS_efficiency'] = results['ESS'] / total_samples
        
        # 4. Geweke检验
        results['geweke'] = self.geweke_all_params()
        
        # 5. 诊断总结
        results['summary'] = self.summarize()
        
        return results
    
    def compute_ess_all_params(self):
        """计算所有参数的ESS"""
        # 对每个参数计算ESS
        ess = {}
        for i in range(self.chains.shape[2]):
            param_samples = self.chains[:, :, i]
            ess[i] = effective_sample_size(param_samples.flatten())
        return ess
    
    def summarize(self):
        """生成诊断总结"""
        r_hat = gelman_rubin(self.chains)
        max_r_hat = np.max(r_hat)
        
        ess = self.compute_ess_all_params()
        min_ess = np.min(list(ess.values()))
        
        if max_r_hat < 1.1 and min_ess > 100:
            return "✓ 收敛良好"
        elif max_r_hat < 1.2:
            return "⚠ 边缘收敛，建议增加迭代"
        else:
            return "✗ 收敛问题，需要检查"

8 与相关内容的联系

8.1 采样方法

收敛诊断适用于所有MCMC方法：

mcmc-methods — MCMC方法基础
mcmc-advanced — NUTS、切片采样
hamiltonian-monte-carlo-geometric-theory — HMC几何理论

8.2 贝叶斯推断

收敛诊断是贝叶斯推断流程的一部分：

bayesian-inference — 贝叶斯推断基础
em-algorithm — EM算法与收敛性

8.3 变分推断

变分推断的收敛诊断不同：

variational-inference-deep-dive — 变分推断
variational-vs-sampling-comparison — 变分vs采样对比

9 总结

本章系统介绍了MCMC收敛诊断的方法：

Gelman-Rubin统计量：比较链间/链内方差
有效样本量：衡量独立样本的信息量
Geweke检验：比较链的不同部分
谱密度分析：诊断混合速度
实用流程：综合多种诊断方法

关键原则：没有任何单一诊断能保证收敛。使用多种方法、保守的参数选择、充足的迭代是最佳实践。

参考文献

Cowles, M.K. & Carlin, B.P. (1996). Markov chain Monte Carlo convergence diagnostics. Journal of the American Statistical Association, 91(434), 883-904. ↩
Gelman, A. & Rubin, D.B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457-472. ↩
Geyer, C.J. (1992). Practical Markov Chain Monte Carlo. Statistical Science, 7(4), 473-483. ↩
Geweke, J. (1992). Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. Bayesian Statistics, 4, 169-193. ↩

Metaphor

探索

MCMC收敛诊断

MCMC收敛诊断

1 引言

2 为什么需要收敛诊断

2.1 MCMC的双重极限

2.2 收敛失败的类型

2.3 诊断原则

3 Gelman-Rubin统计量（R̂）

3.1 基本思想

3.2 正式定义

3.3 R̂计算

3.4 多变量R̂

3.5 代码实现

3.6 阈值选择

4 有效样本量（ESS）

4.1 核心概念

4.2 自相关函数

4.3 ESS估计

4.4 尾部ESS

4.5 ESS阈值

4.6 效率指标

5 Geweke检验

5.1 基本思想

5.2 检验统计量

5.3 代码实现

5.4 解读

6 谱密度分析

6.1 蒙特卡洛标准误

6.2 功率谱密度

6.3 诊断应用

7 实用诊断流程

7.1 推荐流程

7.2 诊断检查清单

7.3 常见问题与解决方案

7.4 完整诊断工具

8 与相关内容的联系

8.1 采样方法

8.2 贝叶斯推断

8.3 变分推断

9 总结

参考文献

Footnotes

关系图谱

目录