贝叶斯数据分析（8）——提高马尔科夫链模拟的效率

本笔记基于Bayesian Data Analysis Third Edition， Andrew Gelman et. al. 学习编写，由于是英文教材，可能在学习过程中一些翻译或者内容有误，如有问题或错误，可以发送至邮箱[email protected]反馈。在学习本书中，需要有一定数学基础或者数理统计的基础，并且存在很多需要计算的场景，对于每一条定理，从头开始的证明会让你更明白每一步是如何实现的。

在上一章里，我们学习了 Gibbs 采样器和 Metropolis-Hastings 算法的基础框架。但在实际应用中，基础的 MCMC 方法常常面临效率瓶颈：参数之间的相关性导致 Gibbs 移动缓慢，随机游走行为让 Metropolis 在高维空间中举步维艰。那么，如何让 MCMC 更快、更高效的去贴合后验分布，这时候，我们需要对 Gibbs 和 Metropolis 的工程优化（重参数化、辅助变量、最优调参），和一种从力学角度出发的**哈密顿蒙特卡洛（Hamiltonian Monte Carlo, HMC）**方法来进行。

高效 Gibbs 采样器

重参数化

Gibbs 采样器的效率高度依赖于参数化方式。最理想的情形是参数之间相互独立——这样 Gibbs 的坐标轴方向移动恰好与目标分布的主轴对齐，每一步都走得高效。而当参数高度相关时（理想情况应该为斜线，因为此时参数之间已经不是正交的），Gibbs 只能在水平和垂直方向交替移动，形成短促的锯齿轨迹，链混合极慢。

最直接的重参数化方法是线性变换，也称为白化变换 (Whitening)。假设 $\theta$ 的后验协方差矩阵为 $\Sigma$（对称正定），定义：

其中 $\Sigma^{-1/2}$ 是 $\Sigma^{-1}$ 的矩阵平方根，满足 $\Sigma^{-1/2}\,\Sigma\,\Sigma^{-1/2} = I$。我们验证 $\phi$ 的协方差矩阵：

结果是单位矩阵——$\phi$ 的每个分量方差为 1，任意两个分量之间的协方差为 0，即各分量线性无关。

这对 Gibbs 采样的意义在于几何层面：

• 变换前：若 $\theta_1$ 与 $\theta_2$ 高度正相关，后验密度的等高线是倾斜的椭圆。Gibbs 的每次更新只能沿坐标轴方向移动——水平一步再垂直一步——每一步被椭圆限制得极短，形成锯齿状轨迹。

• 变换后：$\phi_1$ 与 $\phi_2$ 不相关，等高线被"拉圆"。坐标轴方向恰好对齐了分布的自然主轴，每一步的移动距离大幅增加，混合效率显著提升。

特别地，若后验恰好是多元正态分布，不相关等价于独立——此时 $\phi$ 的各分量完全独立，Gibbs 的效率达到理论最优。对于非正态后验，虽然不能保证严格独立，但去除了线性依赖已大幅改善采样效率。这一变换可以借助主成分分析（PCA）来理解。对协方差矩阵做谱分解 $\Sigma = U \Lambda U^\top$，其中 $U$ 的列是主成分方向，$\Lambda = \text{diag}(\lambda_1, \ldots, \lambda_d)$ 是各方差。则：

于是 $\phi = \Sigma^{-1/2}\theta = \Lambda^{-1/2}(U^\top \theta)$ 实质上是两步操作的叠加：

1. $U^\top \theta$（PCA 旋转）：将坐标轴旋转到主成分方向，消去分量间的相关性。
2. $\Lambda^{-1/2}(\cdot)$（方差归一化）：沿每个主轴方向除以其标准差 $\sqrt{\lambda_i}$，将各方差统一缩放到 1。

所以白化变换 = PCA 旋转 + 球形化（sphering），最终将等高线从倾斜椭圆变为单位圆。（如下图所示，由Nano Banana 2生成）

对 Metropolis 跳跃同样适用这一原则：跳跃核的协方差结构应与目标分布匹配。在众数处做正态近似可以估计 $\Sigma$（见第十三章），然后设置跳跃核为 $\mathcal{N}(\theta^{t-1}, c^2\hat{\Sigma})$。

辅助变量方法

辅助变量（auxiliary variables），也称为数据扩充（data augmentation）：通过引入额外的隐变量，将难以直接采样的复杂分布转化为更简单的条件分布。

这一部分主要是引出下面的参数扩充的一节内容。

e.g.：t 分布表示为正态尺度的混合

t 分布是贝叶斯建模中最常用的厚尾分布之一，但它没有共轭的抽样形式。我们可以使用辅助变量方法来解决。

考虑 $n$ 个独立同分布的观测 $y_i \sim t_\nu(\mu, \sigma^2)$，其中自由度 $\nu$ 已知，$(\mu, \sigma)$ 未知。对每个观测引入一个辅助变量 $V_i$，此时，t 似然等价于如下层次模型：

对于$V_i$ ，我们这么想，每个数据点有一个"属于自己的方差"。当 $V_i$ 很小时，该点对似然的贡献集中，表现为正态分布的行为；当 $V_i$ 很大时，该点被"拉宽"，以容许较大的偏离。

等价性的推导如下。边缘化 $V_i$ 应还原为 t 密度：

$$p(y_i \mid \mu, \sigma^2) = \int_0^\infty p(y_i \mid \mu, V_i) \; p(V_i) \; dV_i$$

首先写出两个密度。正态密度和 $\text{Inv-}\chi^2(\nu, \sigma^2)$ 密度分别为：

$$p(y_i \mid \mu, V_i) = \frac{1}{\sqrt{2\pi V_i}} \exp\!\left(-\frac{(y_i - \mu)^2}{2V_i}\right), \qquad p(V_i) = \frac{(\nu\sigma^2/2)^{\nu/2}}{\Gamma(\nu/2)} V_i^{-(\nu/2 + 1)} \exp\!\left(-\frac{\nu\sigma^2}{2V_i}\right)$$

相乘，将指数项合并：

$$p(y_i \mid V_i)\,p(V_i) = \frac{(\nu\sigma^2/2)^{\nu/2}}{\sqrt{2\pi}\,\Gamma(\nu/2)} \cdot V_i^{-\frac{\nu+3}{2}} \cdot \exp\!\left(-\frac{(y_i - \mu)^2 + \nu\sigma^2}{2V_i}\right)$$

注意到逆 Gamma 积分。令 $A = \frac{(y_i - \mu)^2 + \nu\sigma^2}{2}$，$B = \frac{\nu+1}{2}$，则被积函数正比于 $V_i^{-(B+1)} \exp(-A/V_i)$，这是 $\text{Inv-Gamma}(B, A)$ 的核。利用归一化常数 $\int_0^\infty x^{-(\alpha+1)}e^{-\beta/x}dx = \Gamma(\alpha)/\beta^\alpha$：

$$p(y_i \mid \mu, \sigma^2) = \frac{(\nu\sigma^2/2)^{\nu/2}}{\sqrt{2\pi}\,\Gamma(\nu/2)} \cdot \frac{\Gamma\!\left(\frac{\nu+1}{2}\right)}{\left(\frac{(y_i - \mu)^2 + \nu\sigma^2}{2}\right)^{\frac{\nu+1}{2}}}$$

整理。将分母拆为 $\left(\frac{\nu\sigma^2}{2}\right)^{\frac{\nu+1}{2}} \left(1 + \frac{(y_i - \mu)^2}{\nu\sigma^2}\right)^{\frac{\nu+1}{2}}$，与分子中 $(\nu\sigma^2/2)^{\nu/2}$ 相消得到 $\frac{\sqrt{2}}{\sqrt{\nu}\sigma}$：

$$p(y_i \mid \mu, \sigma^2) = \frac{\Gamma\!\left(\frac{\nu+1}{2}\right)}{\Gamma\!\left(\frac{\nu}{2}\right)\sqrt{\nu\pi}\,\sigma} \left(1 + \frac{(y_i - \mu)^2}{\nu\sigma^2}\right)^{-\frac{\nu+1}{2}}$$

这正是位置-尺度 t 分布 $t_\nu(\mu, \sigma^2)$ 的概率密度函数

在扩充模型 $p(\mu, \sigma^2, V \mid y)$ 上实施 Gibbs 采样分三步：

第一步：更新各 $V_i$。 给定 $y$ 和其他参数，每个 $V_i$ 的后验仍为 scaled inverse-$\chi^2$：

$n$ 个 $V_i$ 在条件后验下相互独立，可以并行抽样。

第二步：更新 $\mu$。 给定各 $y_i$ 及其各自的方差 $V_i$，结合均匀先验：

这是加权最小二乘估计的贝叶斯版本——方差越大的点权重越小。

第三步：更新 $\sigma^2$。 $\sigma^2$ 的信息仅来自 $V_i$：

在扩充空间采样后，只需保留 $(\mu, \sigma)$ 的模拟值，就得到了原始 t 模型的后验样本。但应注意的是，我们的 $V_i$ 和 $\sigma$ 是连锁反应，换句话说，每一个抽样都会导致彼此的改变，这样并不是我们想要的，倘若某次迭代抽到了很接近 0 的 $\sigma$，条件分布会让所有 $V_i$ 也接近 0，进而 $\sigma$ 的条件分布进一步集中在 0 附近。于是乎，我们需要想个办法，来解除这样的连锁反应。

参数扩充

参数扩充（parameter expansion）想法是：通过增加一个冗余参数来改善混合。在更大的空间中采样，然后把不关心的维度边缘化掉。

对于 t 模型，引入一个额外的尺度参数 $\alpha > 0$：

这里 $\alpha^2 U_i$ 扮演了之前 $V_i$ 的角色（即引入了一个新的参数 $\alpha$），而 $\alpha\tau$ 扮演了之前 $\sigma$ 的角色。$\alpha$ 本身没有独立含义，赋予对数尺度上的无信息均匀先验。

扩充后的 Gibbs 有四步（增加了更新 $\alpha^2$ 的步骤）：

虽然 $\alpha^2, U, \tau$ 在模型中不被识别（数据不足以区分它们），但 $V_i = \alpha^2 U_i$ 和 $\sigma = \alpha\tau$ 是被识别的。关键在于 $\alpha$ 打破了 $\sigma$ 和 $V_i$ 之间的连锁——即使 $\tau$ 很小，$\alpha$ 可以很大，保持链的自由移动。

Python 实现：t 分布的辅助变量 Gibbs 采样器

import numpy as np
from scipy.stats import invgamma, norm

def gibbs_t_model(y, nu, n_iter, n_chains):
    """
    使用辅助变量 Gibbs 采样器拟合 t 分布（位置-尺度族）

    模型: y_i ~ t_nu(mu, sigma^2)
    辅助变量: V_i ~ Inv-chi2(nu, sigma^2), y_i | V_i ~ N(mu, V_i)

    参数:
        y: 观测数据 (n,)
        nu: 已知自由度
        n_iter: 迭代次数
        n_chains: 独立链数
    """
    n = len(y)
    chains = {
        'mu': np.zeros((n_iter, n_chains)),
        'sigma2': np.zeros((n_iter, n_chains)),
    }
    acceptance = np.zeros(n_chains)

    for c in range(n_chains):
        # 初始化
        mu = np.median(y)
        sigma2 = np.var(y) * (nu - 2) / nu if nu > 2 else np.var(y)
        V = np.ones(n) * sigma2

        for t in range(n_iter):
            # 步骤1: 更新 V_i | mu, sigma2
            shape = (nu + 1) / 2
            for i in range(n):
                scale = (nu * sigma2 + (y[i] - mu)**2) / (nu + 1) / shape
                V[i] = 1.0 / np.random.gamma(shape, scale)

            # 步骤2: 更新 mu | V, sigma2
            prec_sum = np.sum(1.0 / V)
            mu_mean = np.sum(y / V) / prec_sum
            mu_var = 1.0 / prec_sum
            mu = np.random.normal(mu_mean, np.sqrt(mu_var))

            # 步骤3: 更新 sigma^2 | V, mu
            shape_sigma = n * nu / 2
            rate_sigma = 0.5 * nu * np.sum(1.0 / V)
            sigma2 = 1.0 / np.random.gamma(shape_sigma, 1.0 / rate_sigma)

            chains['mu'][t, c] = mu
            chains['sigma2'][t, c] = sigma2

    return chains

# 示例：从 t_4(0, 5^2) 生成数据并拟合
np.random.seed(42)
y_true = 0 + 5 * np.random.standard_t(4, size=200)

chains = gibbs_t_model(y_true, nu=4, n_iter=2000, n_chains=4)
warmup = 1000

print(f"mu 后验均值: {chains['mu'][warmup:].mean():.3f}")
print(f"sigma2 后验均值: {chains['sigma2'][warmup:].mean():.3f}")
print(f"mu 95% 可信区间: [{np.percentile(chains['mu'][warmup:], 2.5):.3f}, "
      f"{np.percentile(chains['mu'][warmup:], 97.5):.3f}]")

高效 Metropolis 跳跃规则

基础的 Metropolis 算法几乎有无穷多种实现方式——跳跃分布 $J_t$ 的选择决定了算法的成败。好的跳跃规则让链轻盈地穿梭于参数空间；差的跳跃规则让链原地踏步或频繁被拒。

随机游走 Metropolis 的最优缩放

这是 Metropolis 调参中最重要的理论结果。假设有 $d$ 个参数，经过适当变换后后验分布近似为多元正态 $\mathcal{N}(\theta_0, \Sigma)$。如果使用与目标分布同形状的随机游走跳跃核：

这一结论的理论基础是 扩散极限分析：令 $c = \ell/\sqrt{d}$，当 $d \to \infty$ 时，链的单个坐标分量经时间加速后收敛到朗之万扩散过程，其扩散速度由 $h(\ell) = \ell^2 \cdot a(\ell)$ 决定——即步长平方与接受概率的乘积。对标准正态目标，log 接受比的渐近分布为 $\mathcal{N}(-\ell^4/4,\; \ell^4/2)$，由此导出接受概率 $a(\ell) = 2\Phi(-\ell^2/2)$。最大化 $h(\ell)$ 得 $\ell_{\text{opt}} \approx 2.38$，对应 $c_{\text{opt}} \approx 2.38/\sqrt{d}$ 及渐近接受率 $a \approx 0.234$。一维情形（$d = 1$）需单独分析，结果为最优接受率约 44%。

在此最优设置下，Metropolis 的效率约为 $0.3/d$。作为对比，若 $d$ 个参数后验独立，Gibbs 的效率为 $1/d$。

直观理解 23%：在高维空间中，大多数随机方向几乎正交于高密度区域的"脊线"。步长太小时接受率虽高但链原地踏步；步长太大时多数跳跃被拒绝。23% 精确均衡了两者——最大化 $\ell^2 \times$ 接受概率的乘积。

自适应 Metropolis 算法

基于上述理论，一个实用的自适应算法框架如下：

1. 初始阶段：使用固定的初始算法，比如 Gibbs 采样器，或者使用在众数处估计的协方差矩阵 $\hat{\Sigma}$ 并缩放 $c = 2.4/\sqrt{d}$ 的随机游走 Metropolis。

2. 自适应阶段：运行一定迭代后，根据已有的模拟值调整跳跃规则：

   • (a) 将跳跃协方差更新为与已估计的后验协方差矩阵成比例：$J \propto \hat{\Sigma}_{\text{posterior}}$
   • (b) 根据接受率调整尺度：过高则增大 $c$，过低则减小 $c$，目标接受率设为 0.44（一维）或 0.23（多维联合更新）

关键安全规则：自适应只能在预热（warm-up）阶段进行。一旦开始为推断而保留模拟值，跳跃规则必须固定。原因是依赖于历史状态的适应会破坏细致平衡，使链不再收敛到目标分布。

假设你根据当前链的表现增大步长，使得算法在平坦区域走得更快——这听起来合理，但后果是链会不成比例地避免平坦区域、集中在变化剧烈的区域，最终得到的样本分布不再是目标后验。

Python 实现：自适应 Metropolis

import numpy as np
from scipy.stats import multivariate_normal

def adaptive_metropolis(log_target, n_iter, n_chains, d,
                        warmup_ratio=0.5, target_accept=0.23):
    """
    自适应 Metropolis 算法（仅使用 warmup 后的固定核进行推断）

    参数:
        log_target: 对数目标密度函数，签名 f(theta) -> float
        n_iter: 总迭代次数
        n_chains: 链数
        d: 参数维度
        warmup_ratio: 预热比例
        target_accept: 目标接受率（高维 ≈0.23）
    """
    n_warmup = int(n_iter * warmup_ratio)
    chains = np.zeros((n_iter, n_chains, d))
    n_accept = np.zeros(n_chains, dtype=int)

    # 初始化：过离散起始点
    chains[0] = np.random.randn(n_chains, d) * 5

    # 初始跳跃尺度
    c = 2.4 / np.sqrt(d)
    Sigma = np.eye(d)

    log_density = np.array([log_target(chains[0, c]) for c in range(n_chains)])

    for t in range(1, n_iter):
        # Warmup 期间每 50 步自适应一次
        if t < n_warmup and t % 50 == 0 and t > 100:
            recent = chains[max(0, t-200):t].reshape(-1, d)
            Sigma = np.cov(recent.T) + 1e-6 * np.eye(d)
            # 调整尺度以接近目标接受率
            recent_ar = n_accept / t
            if np.mean(recent_ar) > target_accept + 0.05:
                c *= 1.2
            elif np.mean(recent_ar) < target_accept - 0.05:
                c *= 0.8

        for ch in range(n_chains):
            # 提出建议
            proposal = chains[t-1, ch] + c * multivariate_normal.rvs(
                mean=np.zeros(d), cov=Sigma
            )
            log_prop = log_target(proposal)

            # MH 比率（数值稳定版）
            log_r = log_prop - log_density[ch]

            if np.log(np.random.random()) < log_r:
                chains[t, ch] = proposal
                log_density[ch] = log_prop
                n_accept[ch] += 1
            else:
                chains[t, ch] = chains[t-1, ch]

    acceptance_rate = n_accept / n_iter
    return chains, acceptance_rate, c, Sigma

# 示例：30 维相关正态分布
np.random.seed(42)
d = 30
true_cov = np.eye(d)
for i in range(d):
    for j in range(d):
        true_cov[i, j] = 0.7 ** abs(i - j)

def log_target_30d(theta):
    return -0.5 * theta @ np.linalg.solve(true_cov, theta)

chains, ar, final_c, final_sigma = adaptive_metropolis(
    log_target_30d, n_iter=5000, n_chains=4, d=d
)
print(f"最终跳跃尺度 c: {final_c:.3f}")
print(f"理论最优 c ≈ 2.4/√d: {2.4/np.sqrt(d):.3f}")
print(f"接受率: {ar}")
print(f"目标接受率: 0.23")

Gibbs 与 Metropolis 的进一步扩展

切片采样

切片采样（slice sampling）提供了一个优雅的视角。从 $d$ 维目标分布 $p(\theta \mid y)$ 中抽样，等价于在 $(d+1)$ 维空间中，从密度函数曲线下的均匀分布中抽样：

在实际操作中，我们直接在高维连续区域内做均匀采样很难，这时候我们需要用到Gibbs采样的思路，将过程拆解为纵向采样和横向采样。以最常见的一维参数 $\theta$ 的更新过程为例：

**1.纵向采样：确定切片高度：**采样辅助变量 u。

• 在当前位置 $\theta^{(t)}$ 处计算函数值 $f(\theta^{(t)})$。
• 从 Uniform$(0, f(\theta^{(t)}))$ 中随机抽取一个高度 $u$。
• 这个高度 $u$ 定义了一个水平“切片区域” $S_u = \{\theta : f(\theta) \ge u\}$。

**2.扩展区间：外推定位 (Stepping-out)：**找到包含切片区域的初始区间。

由于通常无法直接解析求出 $S_u$ 的精确边界，算法会围绕 $\theta^{(t)}$ 构建一个估计区间 $[L, R]$：

• 设定一个初始步长 $w$。
• 随机将区间 $[L, R]$ 放在 $\theta^{(t)}$ 左右两侧（如 $L = \theta^{(t)} - w \cdot r$，$R = L + w$，其中 $r \sim \text{Uniform}(0, 1)$）。
• 向左扩展：若 $f(L) \ge u$，说明 $L$ 还在切片内，将 $L$ 向左移动 $w$（$L \leftarrow L - w$），直到 $f(L) < u$。
• 向右扩展：若 $f(R) \ge u$，说明 $R$ 还在切片内，将 $R$ 向右移动 $w$（$R \leftarrow R + w$），直到 $f(R) < u$。

**3.横向采样与收缩：收缩拒绝 (Shrinkage)：**采样并自适应缩小区间。

在当前的区间 $[L, R]$ 内均匀抽取候选点 $\theta^* \sim \text{Uniform}(L, R)$。

• 检查候选点：接受：若 $f(\theta^*) \ge u$，说明 $\theta^*$ 在切片 $S_u$ 内部，接受 $\theta^{(t+1)} = \theta^*$，本次迭代结束。

• 拒绝与收缩：若 $f(\theta^*) < u$，说明 $\theta^*$ 落在目标曲线外。此时不丢弃整个流程，而是利用 $\theta^*$ 缩小区间：如果 $\theta^* < \theta^{(t)}$，更新左边界 $L = \theta^*$。如果 $\theta^* > \theta^{(t)}$，更新右边界 $R = \theta^*$。缩小区间后，重新在新的 $[L, R]$ 中抽取 $\theta^*$，重复此步骤直至被接受。

可逆跳转采样：跨越不同维度的空间

有时我们希望马尔可夫链能在不同维度的参数空间之间跳转。典型场景包括：

• 模型平均：链在不同预测变量集合的回归模型之间移动，每个模型的参数维度不同 （e.g. 比较线性模型（2个参数）与二次模型（3个参数））
• 有限混合模型：混合成分数量本身是未知的，需要在不同 $K$ 之间跳转

可逆跳转 MCMC（reversible jump MCMC, RJMCMC）通过引入辅助随机变量 $u \sim q(u)$**[低维转高维]**来实现维度匹配，使得$\dim(\theta_k) + \dim(u) = \dim(\theta_{k'})$，然后通过双射函数合成高维参数：

其核心是保持 Metropolis-Hastings 所需的平衡条件在跨维度时依然成立。

从高维到低维的时候，我们并没有凭空采样新的辅助变量，而是直接对高维参数应用逆映射 $g^{-1}$：

因此，我们通过一个表达式表式如下：

其中，$\mathbb{R}^{d_k}$ 为低维空间，$\mathbb{R}^{d_k‘}$ wei

设 $\mathcal{M}_k$ 为第 $k$ 个候选模型，$\theta_k$ 为其 $d_k$ 维参数向量。不失一般性，假设当前处于低维模型，试图跳转到高维模型（即 $d_k < d_{k^*}$）：

低维到高维跳转（$d_k \to d_{k^*}$，升维）

• 提出新模型与采样辅助变量**：**

  • 从当前模型 $\mathcal{M}_k$ 以概率 $J_{k,k^*}$ 提出新模型 $\mathcal{M}_{k^*}$。
  • 从生成密度 $J(u \mid \theta_k, k, k^*)$ 中采样一个维度为 $d_u = d_{k^*} - d_k$ 的辅助变量 $u$。

• 维度匹配**：**

  • 通过可逆的一一对应映射 $\theta_{k^*} = g_{k,k^*}(\theta_k, u)$，构建高维参数 $\theta_{k^*}$。

• 接受概率：

高维到低维跳转（$d_{k^*} \to d_k$，降维）

• 提出新模型与逆向映射：

  • 从当前高维模型 $\mathcal{M}_{k^*}$ 以概率 $J_{k^*,k}$ 提出低维模型 $\mathcal{M}_k$。
  • 无需采样新的辅助变量，直接应用逆映射 $(\theta_k, u) = g_{k,k^*}^{-1}(\theta_{k^*})$，直接解构得到低维参数 $\theta_k$ 及对应的残差变量 $u$。

• 接受概率：

  $$r = \frac{p(y \mid \theta_k, \mathcal{M}_k)\,p(\theta_k \mid \mathcal{M}_k)\,p(\mathcal{M}_k)}{p(y \mid \theta_{k^*}, \mathcal{M}_{k^*})\,p(\theta_{k^*} \mid \mathcal{M}_{k^*})\,p(\mathcal{M}_{k^*})} \cdot \frac{J_{k,k^*}\,J(u \mid \theta_k, k, k^*)}{J_{k^*,k}} \cdot \left\vert{}\frac{\partial g_{k,k^*}(\theta_k, u)}{\partial (\theta_k, u)}\right\vert{}^{-1}$$

• 然后以概率 $\min(r,1)$ 接受新模型。

三项因子的具体含义：第一项为后验比（Posterior Ratio）： 目标分布在新的状态与当前状态的密度之比。第二项为跳转提议修正比（Proposal Ratio）： 包含两部分：模型转移概率比值（如 $J_{k^*,k} / J_{k,k^*}$）和 辅助变量的生成密度（降维时不需要生成 $u$，因此分子或分母不再包含 $u^*$ 项），有点类似于Metropolis-Hastings算法里的校正跳跃概率。第三项为雅可比行列式（Jacobian Determinant）： 修正由空间映射变换导致的的微元体积变化。降维跳转时的雅可比行列式即为升维映射雅可比行列式的倒数。

模拟回火与并行回火

当后验分布有两个（或多个）被极低密度区域隔开的模式时，标准 MCMC 面临严重困难：链在某个模式内随机游走，任何试图跨越"低谷"的跳跃都因密度比极小而被拒绝。结果：链可能在一个模式内困住成千上万步，从此不再访问其他模式——即使两个模式的后验概率相当。

回火方法的思想来自统计力学：升温让分布变平，低谷消失，跨模式移动变得容易。

温度阶梯

定义 $K+1$ 个分布，对应温度序列 $1 = T_0 < T_1 < \cdots < T_K$：

各项含义：

• $T_0 = 1$：$q_0(\theta) = p(\theta \mid y)$，即原始目标后验。这是我们最终关心的分布。
• $T_k > 1$：指数 $1/T_k < 1$ 压缩了密度的动态范围——高峰被压低，低谷被抬高。极限 $T \to \infty$ 下，$q_\infty(\theta) \to p_0(\theta)$，一个预先选定的、方差很大的"基分布"（如扩散正态）。
• $p_0(\theta)$：基分布，通常取高方差、易于采样的分布。它确保即使最热的分布也有定义良好的形式。

一个具体的数值直觉：若 $p(\theta \mid y)$ 在两个模式处的密度比为 $100:1$，取 $T=10$ 后变为 $100^{1/10} : 1^{1/10} \approx 1.58:1$——几乎被抹平。

模拟回火

模拟回火维护一条链，链的状态是 $(\theta, s)$，其中 $s \in \{0,1,\ldots,K\}$ 指示当前使用的温度层。每步迭代包含两个子步骤：

子步骤 1：在当前温度内更新。 固定 $s_t$，使用以 $q_{s_t}$ 为平稳分布的采样器（如 Metropolis 或 HMC）对 $\theta$ 做一步更新，得到 $\theta^{t+1}$。这一步就是普通的 MCMC——在当前的"热度"下探索。

子步骤 2：提议切换温度层。 以概率 $J_{s_t, j}$ 提议跳到相邻层 $j$（通常 $j = s_t \pm 1$，仅允许相邻层之间跳转）。接受概率为 Metropolis-Hastings 比率：

其中 $c_k$ 是 $q_k$ 的归一化常数的倒数的估计值：$c_k \approx 1 / \int q_k(\theta)\,d\theta$。这些常数之所以需要，是因为各 $q_k$ 的归一化常数不同，而 MH 比率需要未归一化密度的比值。$c_k$ 在预热阶段通过各层被访问的频率自适应估计——目标是让链在各层花费大致相等的时间，既不只在冷层也不只在热层逗留。

最终，仅保留来自 $s=0$（$T=1$，原始目标）的 $\theta$ 样本用于推断。热层的样本全部丢弃——它们只在帮助链跨越模式时发挥作用。

并行回火

并行回火维护 $K+1$ 条链并行运行，第 $k$ 条链在温度 $T_k$ 下演化。它与模拟回火的核心区别在于：各温度层始终在运行，通过"交换状态"而非"跳转"来传递信息。

算法步骤如下：

1. 独立更新：每条链各自用其温度的采样器独立跑一步（这 $K+1$ 步可以并行计算）。

2. 交换提议：周期性地（如每步或每隔几步），对相邻的温度对 $(k, k+1)$ 提议交换状态。即，链 $k$ 当前的 $\theta^{(k)}$ 与链 $k+1$ 当前的 $\theta^{(k+1)}$ 尝试互换。接受概率为：

注意这里不再需要 $c_k$——因为两条链都在各自的 $q_k$ 下，归一化常数在两边的乘积中消去，是并行回火相对于模拟回火的一个重要简化。

交换的直觉：当冷链（$T=1$）困在模式 A 时，相邻的稍热链可能已通过更平坦的分布跳到了模式 B。如果交换被接受，模式 B 的信息就传递到了冷链——就像接力赛中的交接棒。

3. 推断：收敛后，仅使用链 0（$T=1$）的样本。

粒子滤波与遗传算法

这些方法维护一群并行链：

• 粒子滤波：周期性地杀死低概率区域的链，分裂高概率区域的链，使得样本集合更快地向高密度区域集中
• 重要性加权：从不正确的分布 $g(\theta)$ 采样，然后用权重 $p(\theta \mid y)/g(\theta)$ 修正
• 遗传算法：在粒子滤波的基础上加入链之间的"突变"和"交叉"操作

这些方法源自数值优化领域，但在以分布为目标的贝叶斯框架下同样有效。

哈密顿(Hamiltonian)蒙特卡洛（HMC）

Gibbs 和随机游走 Metropolis 的主要效率瓶颈是随机游走行为：每次迭代只在随机方向上试探性地移动一小步，已经走过的方向信息不会被继续利用。随着参数维数升高，链往往需要更多步才能穿过后验分布的典型区域。HMC 的核心改进是引入一个临时的动量变量，利用后验密度的梯度构造一条连续轨迹，从而一次提出距离较远、但仍大概率位于典型区域内的新状态。

动量(Momentum)

设模型有 $d$ 个连续参数，记为 $\theta=(\theta_1,\ldots,\theta_d)^\top$。HMC 为每个参数 $\theta_j$ 配置一个辅助动量 $\phi_j$，组成动量向量 $\phi=(\phi_1,\ldots,\phi_d)^\top$。可以把 $\theta$ 想象成位置，$\phi$ 是动量；速度为 $M^{-1}\phi$。

HMC 把原来的后验分布扩展成位置与动量的联合分布：

动量与参数相互独立，其分布通常取为

其中 $M$ 称为质量矩阵（mass matrix）。定义势能、动能和总能量（哈密顿量）为

这里忽略了与 $\theta$ 和 $\phi$ 无关的归一化常数。于是联合密度也可以写成

因此，保持哈密顿量 $H$ 不变，就等价于沿着联合密度大致不变的轨迹运动。这正是 HMC 能够提出远距离候选点、同时维持较高接受率的原因。

质量矩阵

若 $\theta$ 有 $d$ 个分量，$M$ 就是一个 $d\times d$ 的对称正定矩阵。定义 $\phi\sim\mathcal N(0,M)$，矩阵元素满足

例如，二维参数对应的质量矩阵为

其中，对角元素 $m_{11}$ 和 $m_{22}$ 分别是两个动量分量的方差；非对角元素 $m_{12}$ 是二者的协方差。正定性保证动能始终非负，也保证多元正态动量分布有定义。

实践中常见两种选择：

• 对角质量矩阵 $M=\operatorname{diag}(m_1,\ldots,m_d)$：所有非对角元素均为 0，只校正不同参数的尺度，存储和计算成本较低。
• 稠密质量矩阵：允许非对角元素非零，除尺度外还能校正参数之间的线性相关性，但估计和矩阵运算的成本更高。

质量矩阵的作用可以从位置更新 $d\theta/dt=M^{-1}\phi$ ，即速度，其作用是 HMC 用来调整参数空间几何形状的算法参数。若后验协方差近似为 $\Sigma_\theta$，一个理想化的尺度匹配是

，即参数后验协方差矩阵 (\Sigma_\theta) 的逆矩阵，如果只使用对角矩阵，并且 $\theta_j$ 的后验标准差约为 $s_j$，则可取 $m_j\approx 1/s_j^2$。直观上，后验较宽的方向允许更快地移动，后验较窄的方向移动得更谨慎；经过这种标准化后，各方向完成一次振荡所需的时间尺度更接近。

HMC 的每次迭代由三部分组成：

步骤 1：刷新动量

从辅助动量分布中独立抽取 $\phi$：

这步本质上是 Gibbs 步骤——从联合分布 $p(\theta, \phi \mid y)$ 中给定 $\theta$ 抽取 $\phi$。

步骤 2：蛙跳积分

这是 HMC 的核心。给定当前的 $(\theta,\phi)$，哈密顿动力学由两条微分方程描述：

第一条方程说明动量如何改变位置，第二条说明后验 log 密度的梯度如何改变动量。直接进行解析求解通常不可行，因此让 $(\theta,\phi)$ 经历 $L$ 步“蛙跳”（leapfrog）数值积分，每步步长为 $\epsilon$。它把一次动量更新拆成前后两个半步，从而得到可逆、保持体积且数值上较稳定的离散轨迹。

具体而言，重复 $L$ 次：

(a) 动量的半步更新（离散化第二条哈密顿方程）：

(b) 位置的整步更新（离散化第一条哈密顿方程）：

(c) 使用新位置处的梯度，再完成动量的半步更新：

除了第一步和最后一步外，(c) 和下一轮的 (a) 可以合并为整步。所以实际实现中：先做半个 $\phi$ 步，然后交替 $L-1$ 轮整步的 $\theta$ 和 $\phi$，最后再做半个 $\phi$ 步。

蛙跳积分器的物理解释：

• 当链走向低密度区域时，$\frac{d\log p}{d\theta}$ 在该方向为负，则下次更新的时候 $\phi$ 变小，动量被减速，链逐渐停下然后转向
• 当链走向高密度区域时，$\frac{d\log p}{d\theta}$ 在该方向为正，则下次更新的时候 $\phi$ 变大，动量被加速

连续的哈密顿动力学精确保持 $H(\theta,\phi)$，因而也保持联合密度不变。蛙跳积分在 $\epsilon\to0$ 时趋近精确解；当 $\epsilon$ 有限时会产生小的能量误差，但其对称、辛积分结构通常使误差保持有界振荡，而不是像普通数值积分那样持续单向漂移。

步骤 3：接受-拒绝

设 $(\theta^{t-1}, \phi^{t-1})$ 为蛙跳前的状态，$(\theta^*, \phi^*)$ 为 L 步蛙跳后的状态。计算：

等价地，接受概率可以直接写成哈密顿量之差：

以概率 $\alpha$ 接受 $(\theta^*,\phi^*)$，否则保持原来的 $\theta$。如果能够精确求解哈密顿方程，前后能量完全相同，候选点总会被接受；实际的接受—拒绝步骤是在修正有限步长造成的数值能量误差，使算法仍以正确的后验分布为平稳分布。

处理受限参数

标准 HMC 最适合在连续、可微且 log 密度及其梯度可计算的无约束空间中运行。若参数带有边界，例如标准差 $\sigma>0$，直接在原始空间积分会遇到边界处梯度未定义或数值积分越界的问题。概念上有三种处理方式：

1. 直接拒绝：如果蛙跳过程中密度变为零（如越界），放弃整条轨迹，在当前位置多待一轮。这保持了细致平衡，但浪费迭代。

2. 弹性反弹（bouncing）：当触及边界时，反转动量方向使轨迹反弹回可行域。这通常比直接拒绝更高效。

3. 变换法（实践中的常规做法）：将受限参数变换到无约束空间。例如：

   • 正参数 $\sigma \to \log\sigma$
   • $(0,1)$ 内的概率 $p \to \text{logit}(p) = \log\frac{p}{1-p}$

   变换后需在 log 后验密度中加上 Jacobian 项（$\log|\partial\theta_{\text{orig}}/\partial\theta_{\text{new}}|$），梯度也需相应调整。

HMC 的调参

HMC 有三个调节参数：

如何理解接受率：在高维、近似独立同分布目标和蛙跳积分等理想化渐近条件下，经典理论给出的效率最优接受率约为 65%。它是理解步长量级的参考值，而不是所有 HMC 或 NUTS 实现都必须精确达到的固定目标。实际软件可能采用更保守的目标接受率，并结合发散诊断共同调节。

• 接受率太低或经常发散 → 通常首先说明 $\epsilon$ 太大，应减小 $\epsilon$。
• 接受率非常高但每个有效样本耗时很长 → $\epsilon$ 可能过小，可以在没有发散的前提下适当增大。
• 调整 $\epsilon$ 后，若想维持相近的轨迹长度，可以反向调整 $L$；例如将 $\epsilon$ 减半时把 $L$ 加倍，使总积分时间 $T=\epsilon L$ 大致不变。

在简单的实现中，若后验尺度已由 $M$ 大致标准化，可以把 $\epsilon L\approx1$ 当作轨迹长度的粗略起点；它依赖参数化和质量矩阵，并不是通用常数。质量矩阵和步长的自适应应限制在预热阶段，正式采样时固定自适应结果。也可以按与当前状态无关的预定分布随机抖动 $\epsilon$ 或 $L$，以避免轨迹反复落在特殊周期上。

无 U 型转弯采样器

HMC 最棘手的调参是 $L$：步数太少则轨迹太短，移动不充分；步数太多则轨迹可能绕回来（U 型转弯），浪费计算且可能降低接受率。无 U 型转弯采样器（No-U-Turn Sampler, NUTS）通过在每次迭代中自动决定轨迹长度来解决这个问题。其核心直觉是比较从轨迹起点到当前位置的位移与当前运动方向；当二者的内积表明轨迹开始折返时停止扩展。在 $M=I$ 的简单情形，运动方向可直接用动量表示；一般情形则应结合由质量矩阵决定的速度 $M^{-1}\phi$ 理解。

单纯的"走到转弯处就停"会破坏细致平衡。完整的 NUTS 算法采用更复杂的过程：沿轨迹前后两个方向同时扩展，使用一种满足细致平衡的方式从整条轨迹上的点中抽样。详细的实现不在本章范围内，但核心直觉很简单：每条轨迹应该走到不能再往前走为止。

黎曼自适应 HMC

另一种扩展是让描述局部几何的度量矩阵随位置 $\theta$ 变化，使积分器能够适应后验曲率在不同区域的变化。这时动能、哈密顿方程和数值积分器都会出现额外项，不能简单地把固定的 $M$ 替换成 $M(\theta)$。黎曼流形 HMC 能更精细地处理强烈变化的局部几何，但每一步通常需要昂贵的矩阵分解和隐式积分，因此在高维模型中计算成本很高。

HMC 与 Gibbs 的组合

两种组合方式：

1. 分块 HMC：在分层模型中，按组分别更新 $\theta^{(j)}$ 和超参数 $\tau$——每次 HMC 只涉及一组参数，降低了每次迭代的计算量。参数扩充可以进一步改善跨组的混合。

2. 离散参数的 Gibbs 步骤：哈密顿动力学仅对连续分布有定义。离散参数（如混合成分指示变量、零膨胀模型中的零/非零状态）通过 Gibbs 或一维 Metropolis 更新，连续参数通过 HMC 更新。

分层模型的 HMC 实例

本节以后面讲述的分层模型中的第五章的教育测试八校模型为例，完整演示 HMC 的实现和调参过程。

参数向量有 $d = 10$ 维：$\theta = (\theta_1, \ldots, \theta_8, \mu, \tau)$。

梯度计算

HMC 需要解析梯度。八校模型的 log 后验梯度为：

调试建议：总是用有限差分（扰动 $\pm 0.0001$）检验解析梯度的正确性。两者应对到小数点后多位一致。

质量矩阵与初始化

粗略估计各参数的后验标准差约为 15。按照本文 $\phi\sim\mathcal N(0,M)$ 的约定，使用对角近似时应设置

这里的 $1/225$ 是每个动量分量的方差。四条链的参数起始值仍从独立的 $\mathcal N(0,15^2)$ 中抽取，以形成过离散初始点。参数初始值的分布与辅助动量的分布作用不同：前者决定链从哪里出发，后者决定每次 HMC 轨迹的运动尺度，不能把两者的方差混为一谈。

实际调参过程

这是一个有价值的案例研究：

1. 初始设置：$\epsilon_0 = 0.1, L_0 = 10$，以 $\epsilon L = 1$ 作为教学示例中的初始轨迹长度。每次迭代从均匀分布中随机抖动 $\epsilon \sim U(0, 2\epsilon_0), L \sim U[1, 2L_0]$。4 条链跑 20 步确认程序不崩溃。

2. 第一轮 100 步：接受率分别为 0.23, 0.59, 0.02, 0.57，其中两条链的接受率极低，且 $\widehat{R}>2$。这首先提示步长 $\epsilon$ 太大，蛙跳积分的能量误差过高。

3. 调参：$\epsilon_0$ 减半至 0.05，同时将 $L_0$ 加倍至 20，使总积分时间 $T=\epsilon L$ 大致保持为 1。再跑 100 步，接受率为 0.72, 0.87, 0.33, 0.55。仍未收敛，但数值稳定性有所改善。

4. 延长运行：1000 步后 $\widehat{R} < 1.2$，接受率趋于稳定（0.52, 0.68, 0.75, 0.51）。

5. 最终运行：10,000 步后所有参数的 $\widehat{R} < 1.1$，近似收敛。

变换到log尺度

$\tau > 0$ 的约束可以通过变换 $\xi = \log\tau$ 来处理。此时需要在 log 后验密度中加 Jacobian 项 $\log\tau = \xi$，梯度也需要调整：

Python 实现：简化的 HMC

import numpy as np
from scipy.stats import multivariate_normal

def hmc_sampler(log_prob, grad_log_prob, n_iter, n_chains, d,
                M_diag=None, epsilon=0.1, L=10, randomize_params=True):
    """
    基本的哈密顿蒙特卡洛采样器

    参数:
        log_prob: 对数目标密度 (未归一化), f(theta) -> float
        grad_log_prob: 对数目标密度的梯度, f(theta) -> (d,) array
        n_iter: 迭代次数
        n_chains: 链数
        d: 参数维度
        M_diag: 动量协方差矩阵 M 的对角元 (d,), 默认为全1
        epsilon: 蛙跳步长
        L: 蛙跳步数
        randomize_params: 是否随机抖动 epsilon 和 L
    """
    if M_diag is None:
        M_diag = np.ones(d)
    M_inv = 1.0 / M_diag

    chains = np.zeros((n_iter, n_chains, d))
    n_accept = np.zeros(n_chains, dtype=int)

    # 初始化：过离散起始点
    chains[0] = np.random.randn(n_chains, d) * 5
    log_prob_current = np.array([log_prob(chains[0, c]) for c in range(n_chains)])

    for t in range(1, n_iter):
        for c in range(n_chains):
            theta = chains[t-1, c].copy()
            logp_theta = log_prob_current[c]

            # 步骤1: 按 phi ~ N(0, M) 刷新动量
            phi = np.random.randn(d) * np.sqrt(M_diag)

            # 抖动调参
            eps = epsilon * np.random.uniform(0.5, 1.5) if randomize_params else epsilon
            L_actual = np.random.randint(max(1, L//2), L*2) if randomize_params else L

            # 步骤2: 蛙跳积分
            theta_prop = theta.copy()
            phi_prop = phi.copy()

            # 初始半步步
            grad = grad_log_prob(theta_prop)
            phi_prop = phi_prop + 0.5 * eps * grad

            # L-1 轮整步
            for _ in range(L_actual - 1):
                # d theta / dt = M^{-1} phi
                theta_prop = theta_prop + eps * M_inv * phi_prop
                grad = grad_log_prob(theta_prop)
                phi_prop = phi_prop + eps * grad

            # 最后一步：theta 整步 + phi 半步步
            # d theta / dt = M^{-1} phi
            theta_prop = theta_prop + eps * M_inv * phi_prop
            grad = grad_log_prob(theta_prop)
            phi_prop = phi_prop + 0.5 * eps * grad

            # 步骤3: 接受-拒绝
            logp_prop = log_prob(theta_prop)
            # log p(phi) = -K(phi) + const，其中 K = phi^T M^{-1} phi / 2
            logp_phi_current = -0.5 * np.sum(phi**2 / M_diag)
            logp_phi_prop = -0.5 * np.sum(phi_prop**2 / M_diag)

            log_r = (logp_prop + logp_phi_prop) - (logp_theta + logp_phi_current)

            if np.isfinite(log_r) and np.log(np.random.random()) < log_r:
                chains[t, c] = theta_prop
                log_prob_current[c] = logp_prop
                n_accept[c] += 1
            else:
                chains[t, c] = theta

    return chains, n_accept / n_iter


# 示例：八校模型的简化版本
def eight_schools_log_prob_and_grad():
    """
    构建八校模型的 log 后验和梯度

    参数向量: (theta1..theta8, mu, tau)
    theta_j ~ N(mu, tau^2), y_j ~ N(theta_j, sigma_j^2)
    p(mu) propto 1, p(tau) propto 1
    """
    y = np.array([28, 8, -3, 7, -1, 1, 18, 12])
    sigma = np.array([15, 10, 16, 11, 9, 11, 10, 18])
    J = 8

    def log_prob(theta):
        theta_j = theta[:J]
        mu = theta[J]
        tau = theta[J+1]

        if tau <= 0:
            return -np.inf

        lp = -0.5 * np.sum((theta_j - y)**2 / sigma**2)
        lp += -0.5 * np.sum((theta_j - mu)**2 / tau**2)
        lp += -J * np.log(tau)  # 层级正态密度中的尺度项；平坦先验不另加项
        return lp

    def grad_log_prob(theta):
        theta_j = theta[:J]
        mu = theta[J]
        tau = theta[J+1]

        if tau <= 0:
            return np.full_like(theta, np.nan)

        grad = np.zeros_like(theta)
        grad[:J] = -(theta_j - y) / sigma**2 - (theta_j - mu) / tau**2
        grad[J] = -np.sum((mu - theta_j) / tau**2)
        grad[J+1] = -J / tau + np.sum((mu - theta_j)**2) / tau**3

        return grad

    return log_prob, grad_log_prob, J + 2

np.random.seed(42)
log_prob, grad_log_prob, d = eight_schools_log_prob_and_grad()

chains, ar = hmc_sampler(
    log_prob, grad_log_prob,
    n_iter=5000, n_chains=4, d=d,
    M_diag=np.ones(d) * (1/15**2),  # 后验标准差约为15时，M_diag ≈ 1/15^2
    epsilon=0.05, L=20
)

print(f"HMC 接受率: {ar}")
print(f"mu 后验均值: {chains[2500:, :, 8].mean():.2f}")
print(f"tau 后验均值: {chains[2500:, :, 9].mean():.2f}")

方法总结

在方法上我们大概率可以分为5类：一类是基础的MCMC，然后考虑到效率问题，我们进一步提出了一些提高效率的方法；后面根据扩散极限分析的理论，我们提出了自适应的Metropolis，同时还提及了一些其他的扩展，如切片采样和模拟/并行回火，最后提出了现在MCMC最常用的哈密顿算法。

接受数值

在这一章节里，我们有很多的讨论到接受率的数值，且他们使用方法和范围也不尽一样。通过最优接受率，我们知道如何调节参数，使得效率更高。

HMC 调参

其实在实际应用中，我们会更偏向于直接使用 Stan 或其他自动化的 HMC 实现，理解底层原理的价值在于——当链不收敛时，你知道应该调整什么、如何重参数化、以及诊断输出了什么问题。