贝叶斯数据分析（7）——马尔可夫链模拟基础

如之前一章所说，贝叶斯推断的核心任务是计算后验分布 $p(\theta|y)$ 的各种泛函——均值、分位数、区间估计、预测分布等。对于大多数非平凡模型，后验分布没有解析形式，无法直接进行积分或抽样。此时，马尔可夫链蒙特卡洛（Markov chain Monte Carlo, MCMC）会给予我们很多帮助。

MCMC 的核心思想是：我们不去直接对后验分布抽样，而是构造一条马尔可夫链，让它渐进地以后验分布为平稳分布。运行这条链足够长时间后，链的状态就可以看作来自目标分布的（近似）样本，换句话说，MCMC是从一个近似分布出发，逐步修正，直到逼近目标的过程。

MCMC 的三个关键点：

1. 马尔可夫性：第 $t$ 步的抽取值仅依赖于第 $t-1$ 步，即 $\theta^t \sim T_t(\theta^t|\theta^{t-1})$。换句话说，每一个状态只受其前一个状态影响：$p(\theta^t|\theta^{t-1},\cdots,\theta^1)=p(\theta^t|\theta^{t-1})$
2. 收敛性：转移核 $T_t$ 使得每一步的近似分布在某种意义上离目标分布更近，最终收敛到 $p(\theta|y)$。
3. 平稳分布：当链收敛后，从链上采集的点虽然相关，但都来自同一个平稳分布——后验分布。

本章将系统介绍两种最基础的 MCMC 方法：Gibbs采样器和 Metropolis-Hastings 算法，以及如何评估它们是否收敛。

Gibbs采样器

核心思想

Gibbs采样器的思路极其自然：在多元参数空间中，如果联合分布难以直接抽样，但每个分量在给定其他所有分量的条件下容易抽样，那我们就逐分量地、轮番地进行条件抽样。

形式化地说，将参数向量分解为 $d$ 个子向量 $\theta = (\theta_1, \ldots, \theta_d)$。在第 $t$ 轮迭代中，依次从每个分量的满条件分布（full conditional）中抽样：

其中 $\theta_{-j}^{t-1} = (\theta_1^t, \ldots, \theta_{j-1}^t, \theta_{j+1}^{t-1}, \ldots, \theta_d^{t-1})$ 表示除 $\theta_j$ 外所有分量的当前值。注意这里有一个精妙之处：已经被更新的分量用的是第 $t$ 轮的新值，尚未更新的分量仍保留第 $t-1$ 轮的值。应当注意到，Gibbs 是 Metropolis-Hastings 的特例，原因会在后面给出，但我们需要知道的是，Gibbs采样器相当于构造了一个特殊的跳跃分布——只改变当前分量，且跳跃恰好等于该分量的满条件分布。可以证明这时候的接受概率恒为 1，即总是接受。从这个角度看，Gibbs是最聪明的Metropolis：每一步跳得恰到好处，不需要拒绝。

示例：二元正态分布

考虑来自二元正态总体的单个观测 $(y_1, y_2)$，均值 $\theta = (\theta_1, \theta_2)$ 未知，协方差矩阵已知为 $\begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}$。在均匀先验下，后验分布为：

根据多元正态分布的条件分布性质，两个满条件分布为：

$\theta_1 \mid \theta_2, y \sim \mathcal{N}\big(y_1 + \rho(\theta_2 - y_2),\; 1 - \rho^2\big)$

$\theta_2 \mid \theta_1, y \sim \mathcal{N}\big(y_2 + \rho(\theta_1 - y_1),\; 1 - \rho^2\big)$

Gibbs采样器就是在这两个正态分布之间交替抽样。当 $\rho = 0.8$、数据 $(y_1, y_2) = (0, 0)$ 时，Figure 11.2 展示了四条独立链的演化过程。

Gibbs的每次迭代是分量级的跳跃——只能在水平和垂直方向移动，形成锯齿状的轨迹。当参数之间存在较高相关性时（如 $\rho = 0.8$），这种分量级的更新会导致链移动缓慢，因为每一步只能沿坐标轴方向走，这也是 Gibbs 的一个固有局限。

Python代码实现

import numpy as np
import matplotlib.pyplot as plt

def gibbs_bivariate_normal(y1, y2, rho, n_iter, n_chains, init_dispersed=True):
    """
    Gibbs采样器：二元正态分布
    
    参数:
        y1, y2: 观测数据（已知均值）
        rho: 已知相关系数
        n_iter: 每链迭代次数
        n_chains: 独立链数
        init_dispersed: 是否使用过离散初始点
    
    返回:
        chains: shape (n_iter, n_chains, 2) 的数组
    """
    chains = np.zeros((n_iter, n_chains, 2))
    
    # 初始化：过离散的起始点
    if init_dispersed:
        chains[0, :, 0] = np.random.uniform(-5, 5, n_chains)
        chains[0, :, 1] = np.random.uniform(-5, 5, n_chains)
    else:
        chains[0, :, 0] = np.random.randn(n_chains)
        chains[0, :, 1] = np.random.randn(n_chains)
    
    var_cond = 1 - rho**2  # 条件方差
    
    for t in range(1, n_iter):
        for c in range(n_chains):
            theta2_prev = chains[t-1, c, 1]
            
            # 抽样 theta1 | theta2
            mean1 = y1 + rho * (theta2_prev - y2)
            chains[t, c, 0] = np.random.normal(mean1, np.sqrt(var_cond))
            
            # 抽样 theta2 | theta1 (使用刚抽样的 theta1)
            theta1_new = chains[t, c, 0]
            mean2 = y2 + rho * (theta1_new - y1)
            chains[t, c, 1] = np.random.normal(mean2, np.sqrt(var_cond))
    
    return chains

# 运行Gibbs采样器
np.random.seed(42)
chains = gibbs_bivariate_normal(y1=0, y2=0, rho=0.8, 
                                 n_iter=1000, n_chains=4)

# 可视化收敛后的样本点（后半段）
warmup = 500
samples = chains[warmup:].reshape(-1, 2)

plt.figure(figsize=(10, 8))
plt.scatter(samples[:, 0], samples[:, 1], alpha=0.3, s=10)
plt.xlabel(r'$\theta_1$', fontsize=14)
plt.ylabel(r'$\theta_2$', fontsize=14)
plt.title(f'Gibbs采样器：后500次迭代的点 (ρ=0.8, 4链)', fontsize=16)
plt.axhline(0, color='gray', linestyle='--', alpha=0.5)
plt.axvline(0, color='gray', linestyle='--', alpha=0.5)
plt.axis('equal')
plt.show()

关于Gibbs采样器需要特别注意的一点：它要求我们能够从每个满条件分布中直接抽样。在实践中，这通常意味着模型需要使用共轭先验结构。当某个满条件不是标准分布时，就需要引入Metropolis步骤作为替代。

Metropolis与Metropolis-Hastings算法

Metropolis算法

Metropolis算法的本质是一个带有接受-拒绝规则的随机游走。它在参数空间中试探性地跳跃，然后根据目标密度的相对高低来决定接受还是拒绝这个跳跃。

其算法步骤如下：

1. 选择一个起始点 $\theta^0$，满足 $p(\theta^0|y) > 0$。

2. 对 $t = 1, 2, \ldots$：

   1. 从跳跃分布（proposal distribution）$J_t(\cdot|\theta^{t-1})$ 中抽出一个建议点 $\theta^*$。对基本Metropolis算法而言，跳跃分布必须对称：$J_t(a|b) = J_t(b|a)$ 对所有 $a, b, t$ 成立。（重要前提）

      跳跃性的理解有些疑惑，后面我的理解是： 在交换参数之后，其分布得出来的值不发生变化。 比如正态分布：

      $$J_t(b|a) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(b - a)^2}{2\sigma^2}\right)$$

      由于 $(a-b)^2 = (b-a)^2$，两者恒等。正态分布以当前点为中心时，从 $a$ 走到 $b$ 的概率密度和从 $b$ 走回 $a$ 完全一样，这就叫做对称性。 举个反例，对于对数正态分布：

      $$J_t(\theta^*|\theta^{t-1}) = \text{LogNormal}(\log\theta^{t-1}, \sigma^2)$$

      显然此时 $J_t(a|b) \neq J_t(b|a)$，因为对数正态的概率密度在交换参数后不等价。这时候就需要Metropolis-Hastings算法

   2. 计算密度比：

      $$r = \frac{p(\theta^*|y)}{p(\theta^{t-1}|y)} \\$$

   3. 以概率 $\min(r, 1)$ 接受建议，即：

      $$\theta^t = \begin{cases} \theta^* & \text{以概率 } \min(r, 1) \\ \theta^{t-1} & \text{Others} \end{cases} \\$$

这个接受规则指的是：如果建议点的后验密度更高（$r > 1$），总是接受，因为我们找到了更好的参数值；如果建议点的密度更低（$r < 1$），以概率 $r$ 接受——这意味着密度减半时，接受概率也减半。

这里有一个细节：即使拒绝跳跃，原地不动也计为一次迭代。这意味着链可能在同一点停留多步，特别是在目标分布的尾部区域。

Metropolis算法的有效性

证明分为两步。第一步证明序列构成一个具有唯一平稳分布的马尔可夫链，这要求链满足不可约（irreducible）、非周期（aperiodic）和非瞬态（not transient）的条件——对于定义在合适分布上的随机游走，后两个条件自动满足，不可约性要求跳跃分布最终能以正概率到达任何状态。

第二步证明平稳分布就是目标后验分布。取两个点 $a$ 和 $b$，假设 $p(b|y) \geq p(a|y)$。从 $a$ 转移到 $b$ 的无条件概率密度为：

这里接受概率为1，因为 $b$ 的密度更高。反方向的转移：

由于 $J_t$ 对称（$J_t(b|a) = J_t(a|b)$），两个方向的转移概率相等。联合分布对称意味着 $\theta^t$ 和 $\theta^{t-1}$ 有相同的边际分布，因此 $p(\theta|y)$ 是链的平稳分布。这个论证本质上是细致平衡条件（detailed balance）的应用。

与优化的关系

Metropolis算法可以看作步进式寻模算法（stepwise mode-finding）的随机版本：总是接受上坡（密度增加），但有时也接受下坡（密度减小）。正是这种偶尔接受较差点的能力，使得算法能够探索整个参数空间而不是被困在局部模式上。这与模拟退火有相似之处，但目的不同——Metropolis的目标是完整刻画分布，而非单纯找到最大值。

当我阅读到这里的时候，我的第一直觉是，这跟优化算法里的梯度下降很相似： 当我们把后验密度的负对数定义为损失函数的话：

$$L(\theta) = -\log p(\theta|y)$$

于是乎，Metropolis 的接受比率就变成了：

$$r = \frac{p(\theta^*|y)}{p(\theta^{t-1}|y)} = \exp\big(L(\theta^{t-1}) - L(\theta^*)\big)$$

立刻可以看出接受规则的含义：

• 如果 $\theta^*$ 的损失更小 → $L(\theta^{t-1}) - L(\theta^*) > 0$ → $r > 1$ → 总是接受
• 如果 $\theta^*$ 的损失更大 → $L(\theta^{t-1}) - L(\theta^*) < 0$ → $r < 1$ → 以概率 $r$ 接受

我们在梯度下降法的时候，我们的损失永远不会变大，每次的迭代严格按照梯度下降的方式移动。而Metropolis则允许在其损失函数的形式下损失增大，这里有个好处在于它可以探索整个分布而不被困在局部最小值。但很明显，这也有个问题，他的效率会更低，因为在参数空间的高维区域，随机方向大概率是几乎正交于高密度方向的，因此，大部分的建议也会被拒绝，使得迭代次数很高。

于是乎，我们就有了Metropolis-Adjusted Langevin Algorithm (MALA) 的算法： 其将跳跃分布设为：

$$\theta^* \sim \mathcal{N}\left(\theta^{t-1} + \frac{\epsilon^2}{2}\nabla\log p(\theta^{t-1}|y),\; \epsilon^2 I\right)$$

这里的均值偏移项 $\frac{\epsilon^2}{2}\nabla\log p$ 是梯度的方向，其利用对数后验的梯度信息把建议点往高概率密度区域推。换句话说，他用了梯度的信息指导跳跃的方向，又结合了Metropolis的接受-拒绝机制，因为我们的目标是采样分布，所以我们仍然需要Metropolis来保证正确的平稳分布

Metropolis-Hastings算法

Metropolis-Hastings 对基本 Metropolis 做了两个重要的推广：

1. 跳跃分布可以不对称：不再要求 $J_t(a|b) = J_t(b|a)$。
2. 接受概率修正：用比率之比来纠正不对称性：

这个公式的含义本身非常漂亮。如果跳跃分布倾向于从某个方向提出建议，那么分母中的 $J_t$ 项会自动下调对应的接受概率，从而修正偏差。直觉上，我们把目标密度除以到达该点的难易程度——容易到达的点需要有一个折扣。

当 $J_t$ 恰好等于后验分布本身时（$J(\theta^*|\theta) \equiv p(\theta^*|y)$），$r$ 恒等于1，所有建议都被接受，$\theta^t$ 就是来自 $p(\theta|y)$ 的独立抽样序列。这个极端情况揭示了Metropolis-Hastings设计的优雅：它的目标就是让跳跃分布尽可能接近目标分布，而接受-拒绝规则可以帮助我们做出决策，我们是跳跃还是不跳跃。

良好的跳跃分布应具备的性质

• 易于从中抽样
• 易于计算比率 $r$
• 每次跳跃能走过参数空间中合理的距离（否则随机游走太慢）
• 跳跃不被过于频繁地拒绝（否则太多时间浪费在原地不动）

这些要求之间存在权衡。大步长意味着更容易到达远处但拒绝率更高；小步长意味着高接受率但移动缓慢。

Python代码实现

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

def metropolis_bivariate_normal(n_iter, n_chains, jump_scale=0.2, init_dispersed=True):
    """
    Metropolis算法：二元标准正态分布
    
    目标分布: θ ~ N(0, I_2)
    跳跃分布: θ* | θ^{t-1} ~ N(θ^{t-1}, (jump_scale)^2 * I_2)
    
    参数:
        n_iter: 每链迭代次数
        n_chains: 独立链数
        jump_scale: 跳跃分布标准差
        init_dispersed: 是否使用过离散初始点
    """
    chains = np.zeros((n_iter, n_chains, 2))
    
    # 目标分布的对数密度（相差一个常数）
    def log_target(theta):
        return -0.5 * np.sum(theta**2, axis=-1)
    
    # 初始化
    if init_dispersed:
        chains[0] = np.random.uniform(-5, 5, (n_chains, 2))
    else:
        chains[0] = np.random.randn(n_chains, 2)
    
    log_density_current = log_target(chains[0])
    n_accept = np.zeros(n_chains, dtype=int)
    
    for t in range(1, n_iter):
        for c in range(n_chains):
            # 提出建议
            theta_proposal = chains[t-1, c] + jump_scale * np.random.randn(2)
            log_density_proposal = log_target(theta_proposal)
            
            # 计算接受概率
            log_r = log_density_proposal - log_density_current[c]
            
            # 接受/拒绝
            if np.log(np.random.random()) < log_r:
                chains[t, c] = theta_proposal
                log_density_current[c] = log_density_proposal
                n_accept[c] += 1
            else:
                chains[t, c] = chains[t-1, c]  # 保持原地
    
    acceptance_rate = n_accept / n_iter
    return chains, acceptance_rate

np.random.seed(42)
chains_small, ar_small = metropolis_bivariate_normal(
    n_iter=1000, n_chains=5, jump_scale=0.2
)
print(f"小跳跃尺度 (0.2) 接受率: {ar_small}")

# 使用更合理的跳跃尺度
chains_good, ar_good = metropolis_bivariate_normal(
    n_iter=1000, n_chains=5, jump_scale=1.0
)
print(f"合理跳跃尺度 (1.0) 接受率: {ar_good}")

# 可视化：对比两种跳跃尺度
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

for ax, chains, title in zip(
    axes,
    [chains_small, chains_good],
    [f'跳跃尺度 = 0.2 (接受率: {ar_small.mean():.2f})',
     f'跳跃尺度 = 1.0 (接受率: {ar_good.mean():.2f})']
):
    warmup = 500
    for c in range(chains.shape[1]):
        samples = chains[warmup:, c, :]
        ax.plot(samples[:, 0], samples[:, 1], alpha=0.5, linewidth=0.5)
        ax.scatter(samples[0, 0], samples[0, 1], s=30, marker='o')
    ax.set_xlim(-4, 4); ax.set_ylim(-4, 4)
    ax.set_xlabel(r'$\theta_1$', fontsize=12)
    ax.set_ylabel(r'$\theta_2$', fontsize=12)
    ax.set_title(title, fontsize=13)
    ax.set_aspect('equal')

plt.suptitle('Metropolis算法：跳跃尺度对混合效率的影响', fontsize=15)
plt.tight_layout()
plt.show()

注意到小跳跃尺度（0.2）虽然接受率高（约90%），但链移动极其缓慢，在空间中形成短促的随机游走轨迹。而尺度为1.0时，接受率虽然下降到约50%，但每次跳跃走得更远，混合效率反而更高。

Gibbs与Metropolis的组合使用

实际应用中，很少有模型能纯粹用Gibbs或纯粹用Metropolis来解决。更常见的情况是混合使用：Gibbs和Metropolis是积木，可以按需拼装。

三种混合策略

1. Gibbs优先，Metropolis补充：对条件共轭的参数用Gibbs直接抽样，对非共轭的参数用一维Metropolis跳跃。。

2. 分块更新（block updating）：将参数分成若干块，对每块内部做联合的Metropolis跳跃，块之间仍然保持Gibbs式的条件更新。当块内参数高度相关时，联合更新比逐个分量更新效率高得多。

3. 近似Gibbs加Metropolis修正：当满条件分布 $p(\theta_j|\theta_{-j}, y)$ 难以直接抽样时，用一个近似分布 $g(\theta_j|\theta_{-j})$ 代替。此时跳跃分布变为：

   $$J_{j,t}(\theta^*|\theta^{t-1}) = \begin{cases} g(\theta_j^*|\theta_{-j}^{t-1}) & \text{若 } \theta_{-j}^* = \theta_{-j}^{t-1} \\ 0 & \text{否则} \end{cases} \\$$

   然后使用Metropolis-Hastings比率来决定接受还是拒绝，以纠正近似带来的偏差。

Gibbs与Metropolis-Hastings的关系

还记得我们前面说的，Gibbs是Metropolis-Hastings的一种特例，现在我们来证明他

这个观点的证明非常漂亮。将 Gibbs 采样的第 $j$ 步看作一个 Metropolis-Hastings 步骤，其跳跃分布为：

即只能沿第 $j$ 个分量方向跳跃，且跳跃恰好来自该分量的满条件分布。

代入 Metropolis-Hastings 比率：

由于 $\theta^*$ 和 $\theta^{t-1}$ 仅在 $j$ 分量上不同（$\theta_{-j}^* = \theta_{-j}^{t-1}$），运用条件概率公式 $p(\theta) = p(\theta_j|\theta_{-j}) \, p(\theta_{-j})$：

因此 Gibbs 的每一步接受概率恒为 1。在这个意义上，Gibbs可以被理解为一个永不拒绝的 Metropolis-Hastings 算法——它总能提出一个恰好满足细致平衡的建议。

混合策略的一个实际考量

关于 Gibbs 与 Metropolis 混合使用时的一个经验法则：如果某个参数的满条件分布有共轭形式，优先用 Gibbs；如果有现成的 Metropolis 实现但 Gibbs 更慢，衡量维度——高维空间中 Metropolis 的接受率会自然下降，可能需要更精细的调参。

推断与收敛性评估

迭代模拟的两大挑战

从迭代模拟（而非独立抽样）中做推断面临两个独特问题：

1. 收敛不足：如果迭代次数不够，模拟值可能与目标分布相差甚远。即使大致收敛，早期迭代仍反映起始分布的影响。

2. 序列自相关：收敛后的抽取虽然同分布，但不独立。自相关降低了有效信息量——2500个相关抽取可能只相当于远少于2500个独立抽取的有效信息。

应注意，一个容易产生的误解是认为相关性本身会让推断有偏。事实是：在收敛的前提下，序列相关的抽取仍然是 $p(\theta|y)$ 的同分布样本，做推断时不考虑顺序即可——均值、分位数等汇总统计仍然一致。相关性的影响仅在于效率，即需要更多迭代才能达到独立抽样同等的精度。

Warm-up（预热）期

书中建议丢弃每条链的前半段作为预热期（warm-up），这是一个保守但稳妥的做法。文本还特别解释了为什么不再使用 burn-in 这个术语——后者暗示某种工业化的应力测试逻辑，而 MCMC 的预热更接近一种逐步逼近的过程。

例如：模拟 200 次迭代后舍弃前 100 次；如果仍未收敛，再追加 200 次，此时之前的全部 200 次都视为预热并舍弃。

对于Warm-up的比例，我们通常会选择1/2。然后计算后半段的$\widehat{R}$ 和 $n_{\text{eff}}$，如果 $\widehat{R} > 1.1$，说明收敛还不充分——此时应该追加迭代，而不是加大丢弃比例。新跑的迭代需要把之前所有的迭代（包括之前保留的那些）都视为 warm-up 并一并丢弃，重复直到诊断通过

多条序列与过离散起点

单一序列的收敛诊断极不可靠。Figure 11.3给出了两个警示案例：左侧图中，两条链各自看起来稳定，但放在一起看显然尚未混合到共同分布；右侧图中，两条链覆盖了相同的分布范围，但各自不平稳。这说明收敛既要求各链之间的充分混合，也要求每条链内部的平稳性。

因此推荐的方法是：至少运行两条独立链，从过离散分布中抽取起始点。过离散（overdispersed）意味着起始点的散布范围应该大于目标分布的估计范围——宁可始于太分散，也不能始于太集中而错过整个区域。

Split-$\widehat{R}$ 统计量

这是本章最核心的实用工具。这个统计量将将每条链对半劈开，比较组间方差与组内方差。

记号：令 $\psi_{ij}$ 表示第 $j$ 条链（$j = 1, \ldots, m$）中第 $i$ 个（$i = 1, \ldots, n$）保存的迭代值。$\psi$ 是我们关心的任意标量——可以是单个参数、参数的比值、对数后验密度，或预测值。

计算步骤：

1. 计算链内均值和总体均值：

   链内均值：$\displaystyle \psi_{\cdot j} = \frac{1}{n}\sum_{i=1}^{n} \psi_{ij}$

   总体均值：$\displaystyle \psi_{\cdot\cdot} = \frac{1}{m}\sum_{j=1}^{m} \psi_{\cdot j}$

2. 组间方差：

   $$B = \frac{n}{m-1} \sum_{j=1}^{m} (\psi_{\cdot j} - \psi_{\cdot\cdot})^2$$

   注意 $B$ 中的因子 $n$——因为 $\psi_{\cdot j}$ 是 $n$ 个值的均值，它的方差本应是总体方差的 $1/n$，乘以 $n$ 恢复为总体尺度。

3. 组内方差：

   $$W = \frac{1}{m}\sum_{j=1}^{m} s_j^2, \quad s_j^2 = \frac{1}{n-1}\sum_{i=1}^{n} (\psi_{ij} - \psi_{\cdot j})^2$$

4. 边际后验方差的混合估计：

   $$\widehat{\text{var}}^+(\psi|y) = \frac{n-1}{n}W + \frac{1}{n}B$$

   这个估计量的巧妙之处在于：起始分布过离散时，它高估边际方差（保守）；但在平稳条件下（起始分布等于目标分布时）它是无偏的；在 $n \to \infty$ 的极限下，它收敛到真实的总体方差 $\text{var}(\psi|y)$。

5. 潜在尺度缩减因子：

   $$\widehat{R} = \sqrt{\frac{\widehat{\text{var}}^+(\psi|y)}{W}}$$

当 $\widehat{R}$ 接近 1 时，说明组间方差和组内方差已经相当，链已经充分混合。实际应用中，$\widehat{R}$ < 1.1 通常被视为收敛的合理阈值。需要注意，截断值 1.1 不是显著性检验——没有 p 值，只是实用标准。

Split 操作的意义：将每条链劈成前半和后半，假设有 5 条原始链，劈开后得到 $m=10$ 条半链，每半长度为 $n=250$（假设原始长度 1000，丢弃前半预热）。这样操作同时检验了混合（不同链之间）和平稳性（同链的前后半之间）。

Python代码实现：$\widehat{R}$ 统计量

import numpy as np

def rhat(chains, warmup_ratio=0.5):
    """
    计算 split-R-hat 统计量
    
    参数:
        chains: shape (n_iter, n_chains) 或 (n_iter, n_chains, n_params) 的数组
        warmup_ratio: 预热期比例
    
    返回:
        R_hat: 每条参数链的潜在尺度缩减因子
    """
    # 丢弃预热期
    n_iter = chains.shape[0]
    warmup = int(n_iter * warmup_ratio)
    post_warmup = chains[warmup:]
    
    n_half = post_warmup.shape[0] // 2
    n_chains = post_warmup.shape[1]
    
    # 将每条链劈成两半，重新组织为 m = 2 * n_chains 条链
    first_half = post_warmup[:n_half]
    second_half = post_warmup[n_half:2*n_half]
    
    # 堆叠: (n_half, 2 * n_chains) 或更高维
    split_chains = np.concatenate([first_half, second_half], axis=1)
    
    n = split_chains.shape[0]   # 每条半链的长度
    m = split_chains.shape[1]   # 半链数量
    
    # 处理多维情况
    if split_chains.ndim == 2:
        split_chains = split_chains[..., np.newaxis]
    
    n_params = split_chains.shape[2]
    R_hats = np.zeros(n_params)
    
    for p in range(n_params):
        psi = split_chains[:, :, p]  # (n, m)
        
        # 链内均值
        psi_dot_j = np.mean(psi, axis=0)      # (m,)
        psi_dot_dot = np.mean(psi_dot_j)       # 标量
        
        # 组间方差 B
        B = n / (m - 1) * np.sum((psi_dot_j - psi_dot_dot)**2)
        
        # 组内方差 W
        s2_j = np.var(psi, axis=0, ddof=1)     # (m,)
        W = np.mean(s2_j)
        
        # 混合方差估计
        var_plus = (n - 1) / n * W + B / n
        
        # R-hat
        R_hats[p] = np.sqrt(var_plus / W)
    
    return R_hats.squeeze()

# 在Metropolis示例上测试R-hat
np.random.seed(42)
chains_test, _ = metropolis_bivariate_normal(
    n_iter=1000, n_chains=5, jump_scale=0.2
)

# 分别计算theta_1和theta_2的R-hat
rhat_values = rhat(chains_test)
print(f"跳跃尺度=0.2时 R-hat: theta1 = {rhat_values[0]:.3f}, theta2 = {rhat_values[1]:.3f}")

# 合理跳跃尺度
chains_good, _ = metropolis_bivariate_normal(
    n_iter=5000, n_chains=5, jump_scale=1.0
)
rhat_good = rhat(chains_good)
print(f"跳跃尺度=1.0时 R-hat: theta1 = {rhat_good[0]:.3f}, theta2 = {rhat_good[1]:.3f}")

关于非正态分布的稳健处理

$\widehat{R}$ 统计量基于均值和方差，对近乎正态的边际后验最有效。当参数分布远离正态时——例如正值参数的对数正态、有界参数（如概率）、或重尾分布——书中建议先做变换再计算 $\widehat{R}$。具体策略：

• 对正值参数：取对数
• 对有界参数（0到1）：取 logit
• 对重尾分布：用秩变换（rank transformation）

变换的目的是让分布更接近正态，从而使基于一阶矩和二阶矩的诊断更可靠，而不是对原尺度上的偏态分布产生误判。

有效样本量

自相关影响估计精度

当 MCMC 链收敛后，每个抽取 $\psi_{ij}$ 都是来自后验的无偏样本，但它们不是独立的。问题出在我们关心的汇总统计上——用链均值 $\bar{\psi}_{\cdot\cdot}$ 估计后验期望 $\mathbb{E}(\psi|y)$ 时，正自相关会放大其方差。

对于一个平稳的、长度为 $n$ 的相关序列，其均值的渐近方差为：（其中 $m$ 为链的数量）

其中 $\rho_t$ 是序列在滞后 $t$ 处的自相关系数。

推导过程： 考虑一条长度为 $n$ 的平稳序列 $\psi_1, \psi_2, \ldots, \psi_n$，均值为 $\bar{\psi} = \frac{1}{n}\sum_{i=1}^n \psi_i$。对这条单链而言：

$$\text{var}(\bar{\psi}) = \text{var}\!\left(\frac{1}{n}\sum_{i=1}^n \psi_i\right) = \frac{1}{n^2} \; \text{var}\!\left(\sum_{i=1}^n \psi_i\right)$$

那么，使用方差分解

$$\text{var}\left(\sum_{i=1}^n \psi_i\right) = \sum_{i=1}^n \text{var}(\psi_i) + 2\sum_{i=1}^{n-1}\sum_{j=i+1}^n \text{cov}(\psi_i, \psi_j)$$

其中 $\text{cov}(\psi_i,\psi_j)$ 为协方差

此时，还记得马尔可夫链的平稳性，对于所有的 $\psi_i$ 有相同方差 $\sigma^2 = \text{var}(\psi|y)$，且$\psi_i$ 和 $\psi_j$ 的协方差只取决于它们的间距 $|i-j|$，即 $\text{cov}(\psi_i, \psi_j) = \rho_{|i-j|} \cdot \sigma^2$，代入原式：

$$\text{var}\left(\sum_{i=1}^n \psi_i\right) = n\sigma^2 + 2\sigma^2\sum_{i=1}^{n-1}\sum_{j=i+1}^n \rho_{j-i}$$

把内层和按滞后 $t = j-i$ 重新组织。对于固定的 $t$，序列中有 $n-t$ 组点对间隔恰好为 $t$：原式则：

$$\begin{aligned} \text{var}\left(\sum_{i=1}^n \psi_i\right) &= n\sigma^2 + 2\sigma^2\sum_{t=1}^{n-1} (n-t)\rho_t \\ &= \sigma^2\left[n + 2\sum_{t=1}^{n-1} (n-t)\rho_t\right] \end{aligned}$$

除以 $n^2$ 得到样本均值的方差：

$$\text{var}(\bar{\psi}) = \frac{\sigma^2}{n}\left[1 + 2\sum_{t=1}^{n-1} \left(1 - \frac{t}{n}\right)\rho_t\right] \Rightarrow\\ n \cdot \text{var}(\bar{\psi}) = \sigma^2\left[1 + 2\sum_{t=1}^{n-1} \left(1 - \frac{t}{n}\right)\rho_t\right]$$

取 $n \to \infty$，因子 $\left(1 - \frac{t}{n}\right) \to 1$（对任意固定的 $t$），于是：

$$\lim_{n\to\infty} n \cdot \text{var}(\bar{\psi}) = \sigma^2\left(1 + 2\sum_{t=1}^{\infty} \rho_t\right)$$

最后，有 $m$ 条独立链时，总均值 $\bar{\psi}_{\cdot\cdot}$ 的方差是单链的 $1/m$。两边乘以 $m$：

$$\lim_{n\to\infty} mn \cdot \text{var}(\bar{\psi}_{\cdot\cdot}) = \left(1 + 2\sum_{t=1}^{\infty} \rho_t\right) \text{var}(\psi|y)$$

Q.E.D

如果抽取是独立的（$\rho_t = 0$ 对所有 $t>0$），那么

样本量就是 $mn$。

由此定义有效样本量：

这个有效样本量的具体含义是：这个相关样本的信息量相当于多少个独立抽取？

这时候我们就可以开始分析，这个有效样本量的变化因素有哪些：

$\rho_t > 0$ ： 此时分母 >0, 且大于1，均值的方差被放大了。换句话说，如果当前点高于均值，下一点大概率也高于均值。这导致样本均值的波动更大，因为整条链可能在某个区域逗留很久，漂移缓慢。实际获得的信息少于 $mn$，有效样本量自然就打折扣。

$\rho_t < 0$：此时分母<1, 均值的方差反而比独立情况更小。这是因为负自相关意味着值在交替振荡，使得局部平均更快地收敛到总体均值。不过在 MCMC 的实践中，负自相关极少出现，几乎总是面对正自相关。

除此之外，我们还很清晰的知道，$\sum \rho_t$ 必须收敛，如果发散的话，我们将不能求得一个稳定的有效样本量。换句话说，这个条件等价于要求序列是短记忆的（short memory）——自相关衰减得足够快，使得无穷的和收敛于一个确定值。对于实际中的 MCMC 链，只要链满足几何遍历性（geometric ergodicity），$\rho_t$ 呈指数衰减，这个条件是成立的。

几何遍历性，geometric ergodicity： 用 $P^t(\theta, \cdot)$ 表示从初始点 $\theta$ 出发、经过 $t$ 步后链所在状态的分布。用 $\pi(\cdot)$ 表示目标平稳分布 $p(\cdot|y)$。两者之间的差异用全变差距离（total variation distance）来度量：

$$\|P^t(\theta, \cdot) - \pi(\cdot)\|_{\text{TV}}$$

几何遍历性：存在一个常数 $r \in (0, 1)$ 和一个函数 $M(\theta)$，使得对任意起始点 $\theta$ 和任意 $t$：

$$\|P^t(\theta, \cdot) - \pi(\cdot)\|_{\text{TV}} \leq M(\theta) \cdot r^t$$

• $r^t$：指数衰减。每一步，当前分布和平稳分布之间的差距被一个小于 1 的因子 $r$ 压缩。
• $M(\theta)$：起始点的惩罚因子。从不同的 $\theta$ 出发，初始的差距不同，$M(\theta)$ 把这个初始差距考虑进去。从参数空间的极端角落出发，$M(\theta)$ 很大；从高概率密度区域出发，$M(\theta)$ 很小。
• 全变差距离：衡量两个概率分布之间最坏情况差异的标准度量:

$$\|P - Q\|_{\text{TV}} = \sup_{A}\; \big|P(A) - Q(A)\big|$$

右侧的 $\sup$ 取遍所有可测集合 $A$。遍历参数空间中所有可能的子集，找出使得让两个分布的概率差距最大的子集，那个最大的差距就是全变差距离。等价形式：

$$\|P - Q\|_{\text{TV}} = \frac{1}{2} \int |p(\theta) - q(\theta)| \, d\theta$$

变异函数（Variogram）估计自相关

直接估计 $\rho_t$ 对大的 $t$ 很不稳定（样本相关性在大滞后下噪声极大）。在这里介绍一个巧妙的方法：先计算变异函数 $V_t$，再反推 $\rho_t$。

变异函数衡量间隔 $t$ 步的两个值之间的期望平方差：

对于平稳序列，有 $\mathbb{E}[(\psi_i - \psi_{i-t})^2] = 2(1 - \rho_t)\text{var}(\psi)$，因此：

其中总体方差 $\widehat{\text{var}}^+$ ，：

$B$ 为组间方差， $W$ 为组内方差

截断策略：不是把所有 $\hat{\rho}_t$ 都加起来，因为远距离的估计噪声会更大，所以我们从滞后0开始累加，直到连续两个滞后的自相关之和为负时停止。即找到第一个奇正整数 $T$，使得 $\hat{\rho}_{T+1} + \hat{\rho}_{T+2} < 0$。这个规则来自 Geyer (1992)，直觉是：当自相关估计的噪声超过了信号，就不再累加。

最终的估计量为：

有效样本量取值

一个问题是：$n_{\text{eff}}$ 多大可以满足我们的计算，得到一个可信的结果。一个非常实用的经验法则是：每条链至少需要 $n_{\text{eff}} \geq 10$（即每半链至少 5 个有效样本）。对于大多数后验均值和方差的估计，100 个独立抽取已经足够——如果 $n_{\text{eff}} = 10$，标准误的膨胀因子仅为 $\sqrt{1 + 1/10} = 1.05$，几乎可以忽略。

但也有例外：如果关心的是后验尾部的极值（如 97.5% 分位数），就需要更大的有效样本量。

Python代码实现：$n_{\text{eff}}$ 的计算

import numpy as np

def effective_sample_size(chains, warmup_ratio=0.5):
    """
    计算有效样本量 n_eff
    
    参数:
        chains: shape (n_iter, n_chains) 或 (n_iter, n_chains, n_params)
        warmup_ratio: 预热比例
    
    返回:
        n_eff: 有效样本量
    """
    n_iter = chains.shape[0]
    warmup = int(n_iter * warmup_ratio)
    post_warmup = chains[warmup:]
    
    # Split chains
    n_half = post_warmup.shape[0] // 2
    first_half = post_warmup[:n_half]
    second_half = post_warmup[n_half:2*n_half]
    split_chains = np.concatenate([first_half, second_half], axis=1)
    
    n = split_chains.shape[0]  # 半链长度
    m = split_chains.shape[1]  # 半链数量
    
    if split_chains.ndim == 2:
        split_chains = split_chains[..., np.newaxis]
    
    n_params = split_chains.shape[2]
    n_effs = np.zeros(n_params)
    
    for p in range(n_params):
        psi = split_chains[:, :, p]  # (n, m)
        
        # 计算 var^+（与R-hat中的相同）
        psi_dot_j = np.mean(psi, axis=0)
        psi_dot_dot = np.mean(psi_dot_j)
        B = n / (m - 1) * np.sum((psi_dot_j - psi_dot_dot)**2)
        s2_j = np.var(psi, axis=0, ddof=1)
        W = np.mean(s2_j)
        var_plus = (n - 1) / n * W + B / n
        
        # 计算变异函数和自相关
        rho_hat = np.zeros(n)
        rho_hat[0] = 1.0  # 滞后0
        
        for t in range(1, n):
            # 变异函数 V_t
            V_t = 0.0
            for j in range(m):
                for i in range(t, n):
                    V_t += (psi[i, j] - psi[i-t, j])**2
            V_t /= (m * (n - t))
            
            # 自相关估计
            rho_hat[t] = 1.0 - V_t / (2.0 * var_plus)
        
        # 截断规则：找到第一个奇数T使得ρ_{T+1}+ρ_{T+2}<0
        T = 1
        while T + 2 < n:
            if rho_hat[T+1] + rho_hat[T+2] < 0:
                break
            T += 2
        
        # 有效样本量
        rho_sum = np.sum(rho_hat[1:T+1])
        if 1 + 2 * rho_sum <= 0:
            n_effs[p] = m * n  # 退化情况
        else:
            n_effs[p] = m * n / (1 + 2 * rho_sum)
    
    return n_effs.squeeze()

# 测试
np.random.seed(42)
chains_good, ar = metropolis_bivariate_normal(
    n_iter=5000, n_chains=5, jump_scale=1.0
)
n_eff = effective_sample_size(chains_good)
print(f"跳跃尺度1.0的有效样本量: {n_eff}")
print(f"原始迭代数: {chains_good.shape[0] * chains_good.shape[1]}")

停止准则

停止规则可以总结为：

1. 对所有感兴趣的标量汇总，计算 split-$\widehat{R}$，确保都小于 1.1。
2. 对所有这些量，确保 $\hat{n}_{\text{eff}} \geq 5m$（半链数量 $m$ 的5倍，即每条原始链至少有10个有效独立抽取）。
3. 两个条件都满足后，将所有预热后的迭代合并，当作来自目标分布的样本使用。
4. 如果不满足，可能需要调整算法本身（如改变跳跃尺度、重参数化），而非仅仅增加迭代次数。

最后需要注意的是：即使 $\widehat{R}$ 和 $n_{\text{eff}}$ 都通过了，也不能百分之百保证收敛——如果目标分布存在起始分布未覆盖且算法难以到达的重要区域，当前的诊断无法发现。这种情形下需要更深入的模型理解和探索性分析。

总结