贝叶斯数据分析（6）——贝叶斯计算导论

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

我们之前提到的，都是一些简单模型，我们可以直接写出后验分布，但当后验分布不能解析写出时，我们如何计算贝叶斯推断需要的量？

在贝叶斯分析里，我们通常关心两类对象：

1. 参数的后验分布：

2. 未来观测值或者未观测量的后验预测分布：

正如我们之前单参数模型和多参数模型提到的，如果我们容易得到后验分布，比如共轭先验下的二项-Beta、Poisson-Gamma、正态模型，我们可以直接写出后验分布。但是一旦模型更复杂，尤其是层级模型、高维模型、非线性模型、缺失数据模型，解析积分几乎就不可行了。因此，贝叶斯计算的本质其实为：积分

我们知道后验均值、后验方差、后验区间、后验概率、预测分布、模型检查量，几乎都可以写成某种积分：

所以贝叶斯计算并不只是从后验抽样，它真正要解决的是：如何在复杂的高维概率分布下，稳定地估计这些积分。

基本概念

归一化密度和未归一化密度

我们把想要模拟或者积分的分布叫做目标分布（target distribution），记为：$p(\theta\mid y)$。但是在实际计算中，我们常常不直接知道归一化后的 $p(\theta\mid y)$，而只知道一个与其成比例的函数：$q(\theta\mid y)$

其中：

这里的 $q(\theta\mid y)$ 就是未归一化密度（unnormalized density）。在贝叶斯定理中，最典型的未归一化密度是：

因为在贝叶斯公式里：

而边际似然：$p(y)=\int p(\theta)p(y\mid \theta)d\theta$ 通常难以计算。但很多算法不需要知道 $p(y)$，因为它只与数据 $y$ 有关，不随 $\theta$ 改变。我们在实际应用时，我们并不追求归一化的密度。

很多抽样算法只需要比较两个参数值 $\theta$ 和 $\theta'$ 哪一个后验密度更高，此时归一化常数会抵消：

$$\frac{p(\theta'\mid y)}{p(\theta\mid y)}=\frac{q(\theta'\mid y)}{q(\theta\mid y)}$$

这就是为什么现代贝叶斯软件通常只要求用户写出 log posterior kernel，即 $\log q(\theta\mid y)$。

使用对数密度

在实际计算里，应该尽可能使用对数密度：$\log q(\theta\mid y)$

原因很简单：概率密度、似然值、后验密度常常是许多小数相乘，直接计算容易出现 underflow；有时也会因为指数项过大出现 overflow。例如，如果：

那么计算对数以后：

乘法变成加法，数值上更稳定，不会出现数值超过计算机的最大长度。如果算法中需要计算两个密度的比值，也先计算 log ratio：

现代计算中 log density 是默认语言。

Stan、PyMC、NumPyro、Turing 等概率编程系统本质上都是让用户定义一个 log density，然后由算法自动完成梯度、采样或者近似。特别是在 Hamiltonian Monte Carlo 中，算法需要 $\nabla_\theta \log q(\theta\mid y)$，所以自动微分和 log density 几乎是现代贝叶斯计算的基础设施。

数值积分

贝叶斯推断中最常见的统计量是后验期望：

这里的 $h(\theta)$ 可以是任意我们关心的函数。例如：

• 如果 $h(\theta)=\theta$，就是后验均值；
• 如果 $h(\theta)=I(\theta>0)$，就是后验概率 $\Pr(\theta>0\mid y)$；
• 如果 $h(\theta)=(\theta-E(\theta\mid y))^2$，就是后验方差；
• 如果 $h(\theta)$ 代表未来预测量，就是预测分布的总结。

模拟方法

如果我们可以从后验分布中抽到 $S$ 个样本：$\theta^1,\theta^2,\ldots,\theta^S\sim p(\theta\mid y)$

那么可以用样本平均估计后验期望：

这里上标 $s$ 表示第 $s$ 次模拟抽样。如果 $\theta^s$ 是从 $p(\theta\mid y)$ 独立抽出的，那么：

也就是说，样本平均是后验期望的无偏估计。它的方差为：

所以 Monte Carlo 误差大约是：

这给了我们一个非常直接的判断：抽样数 $S$ 越大，计算误差越小，但误差下降速度只有 $1/\sqrt S$。

$1/\sqrt S$ 是 Monte Carlo 的温柔和残酷。

它温柔在于：维度再高，只要能从目标分布抽样，样本平均仍然有这个基本收敛速度。它残酷在于：如果想把误差缩小 10 倍，通常需要 100 倍的抽样量。因此真正高效的贝叶斯计算，不是盲目增加 $S$，而是让每一个样本更接近有效样本。

确定性积分方法

确定性数值积分不是随机抽样，而是在若干点 $\theta_s$ 上计算积分函数，并用权重 $w_s$ 加权：

如果使用未归一化密度，则更常见的形式是：

低维情况下，网格积分、Simpson 法、高斯求积等方法可能很有效。但在高维空间中，网格点数量会迅速爆炸。假设每一维取 $N$ 个网格点，维度为 $d$，那么总网格点数就是：$N^d$.

如果 $N=100,d=10$，那么需要 $100^{10}=10^{20}$ 个点，这显然不可行，出现维度爆炸的问题。在大数据的时代，这看着很鸡肋，但这放在现在仍然很有用，例如 Laplace approximation、积分嵌套拉普拉斯近似（INLA）、变分推断、稀疏网格、贝叶斯求积等。但它们往往依赖特殊结构，或者作为 MCMC 的初始化、加速、诊断工具。

分布近似

分布近似的想法是：如果真正的后验分布 $p(\theta\mid y)$ 太复杂，我们可以用一个更简单的分布 $g(\theta)$ 去近似它，然后在 $g(\theta)$ 下进行积分、抽样或者初始化。

最经典的是正态近似。设：

也就是后验众数。对 $\log q(\theta\mid y)$ 在 $\hat\theta$ 附近做二阶 Taylor 展开：

其中：

如果 $H$ 是正定矩阵，则：

这就是 Laplace approximation 的基本形式。归一化常数也可以近似为：

分布近似的优点和代价。

正态近似快、清楚、容易解释，但是它会压平偏态、厚尾和多峰结构。对于简单、数据量大、后验接近正态的问题，它可能很好；对于层级模型、弱识别模型、混合模型，它可能非常误导。因此分布近似常常在一开始接触数据的时候使用，可能更多的看成一个预实验的角度出发去检验数据分布。

另外，这个在参数估计方法里也经常用到，比如Delta法

粗略估计：先获得后验的位置感

在使用复杂算法之前，先用简单方法得到目标分布的大致位置。

例如：

• 在层级模型中，先粗略估计超参数 $\phi$，再估计条件后验中的参数 $\gamma\mid \phi,y$；
• 在缺失数据问题中，先用简单插补获得初始值；
• 在复杂回归模型中，先拟合简化模型，观察参数数量级和方向；
• 在高维模型中，先找到后验众数或者近似均值，作为抽样算法的起点。

直接模拟

在简单模型中，我们可以直接从后验分布抽样。例如在共轭模型中：

此时只需要调用对应的随机数生成器。

更复杂一点的情况是例如多参数模型里把联合后验分解为若干部分。例如：

如果我们可以先从 $p(\phi\mid y)$ 抽样，再从 $p(\theta\mid \phi,y)$ 抽样，那么就可以得到联合后验抽样：

这个思想在后续的层级模型中非常重要。它本质上是利用模型的条件结构，把一个难抽的联合分布拆成更容易处理的部分。

网格近似和离散抽样

对于低维参数，可以在一组网格点上计算目标密度：

然后用离散分布近似连续后验：

如果网格点不是等间距，更合理的做法是把每个点代表的体积也放进权重：

一旦有了离散概率，就可以通过**逆 CDF **方法抽样：先抽一个 $U\sim \operatorname{Uniform}(0,1)$，再找到使累计概率超过 $U$ 的网格点。

网格法适合学习，不适合高维通用计算。

网格法适合帮助我们理解后验形状，尤其是 1 维或 2 维参数。但它不适合高维模型，因为维度灾难会让计算量迅速爆炸。现代实践中，网格更常用于可视化、低维的敏感性分析，或者用于检查更复杂算法是否给出合理结果。

后验预测分布的模拟

一旦我们已经从后验分布得到样本：

后验预测分布就很容易模拟。对每一个 $\theta^s$，从条件预测分布中抽一个：

于是：$\tilde y^1,\tilde y^2,\ldots,\tilde y^S$ 就可以看作来自 $p(\tilde y\mid y)$ 的模拟样本，其背后的积分是：

最后一个等号使用的是条件独立：给定参数 $\theta$ 后，未来数据 $\tilde y$ 与已经观测到的数据 $y$ 独立。如果我们只用 $\hat\theta$ 去预测，那么忽略了参数不确定性。后验预测分布会把两层不确定性都放进去：一层是 $\theta$ 的后验不确定性，另一层是即使 $\theta$ 已知，未来观测 $\tilde y$ 也有随机波动。

拒绝采样

拒绝采样（rejection sampling）用于从目标分布 $p(\theta\mid y)$ 中产生独立样本。它需要一个容易抽样的近似函数 $g(\theta)$，并要求存在一个常数 $M$，使得：

直觉上，$Mg(\theta)$ 是覆盖目标密度 $p(\theta\mid y)$ 的上包络函数。且应注意，$\int g(\theta) d\theta = C$

图中上方曲线是 $Mg(\theta)$，下方曲线是目标密度 $p(\theta\mid y)$。在某个候选点处，接受概率等于下方高度与上方高度的比值。

拒绝采样算法可以写成：

1. 从与 $g(\theta)$ 成比例的分布中抽取候选值 $\theta$；

2. 以概率 $\frac{p(\theta\mid y)}{Mg(\theta)} \\$ 接受这个 $\theta$；

3. 如果拒绝，就回到第一步。

接受后的样本服从目标分布，假设 $g(\theta)$ 本身已经是归一化密度。候选值为 $\theta$ 且被接受的联合密度正比于：

因此：

即证明了拒绝采样的正确性。

如果 $g(\theta)$ 只是未归一化函数，设：

那么从 $g$ 成比例的分布抽样就是从 $r$ 抽样，推导只是多了一个常数 $G$，最终仍然会抵消。

拒绝采样的效率

拒绝采样的接受率为：

这里假设 $g$ 是归一化密度。如果 $M$ 很大，接受率就很低。理想情况下，$g(\theta)$ 与 $p(\theta\mid y)$ 成比例，此时可以用很小的 $M$，几乎所有候选点都会被接受。

拒绝采样的优点是样本独立，且算法“自我监控”：如果接受率很低，我们马上知道近似函数 $g$ 不好。缺点是高维问题中很难找到合适的包络函数，$M$ 往往巨大，导致大量候选样本被浪费。

如果用一个近似分布 $g$ 来帮助计算，$g$ 必须覆盖目标分布的重要区域，尤其不能尾部太薄。

重要性采样

重要性采样（importance sampling）不是直接从目标分布抽样，而是从一个容易抽样的近似分布 $g(\theta)$ 中抽样，然后用权重修正。

我们希望计算：

如果从 $g(\theta)$ 抽样，那么：

其中重要性权重为：

从 $g(\theta)$ 中抽取 $\theta^1,\ldots,\theta^S\sim g(\theta)$

则有自归一化重要性采样估计量：

归一化重要性权重为：

则估计量可以写成更直观的加权平均：

重要性采样的有效样本量

如果所有权重都相等，那么每个样本贡献一样，重要性采样接近普通 Monte Carlo。如果少数样本权重极大，那么估计几乎由少数点决定，实际有效样本量会很低。

有效样本量近似：

若所有权重相等，$\tilde w_s=1/S$，则：（这里等权重会导致信息丢失和纳入重复抽样，但即便如此，这也很方便进行计算）

如果只有一个样本权重几乎为 1，其余几乎为 0，则：

也就是说，虽然你抽了 $S$ 个样本，但实际上只有一个样本在发挥作用。

权重尾部问题

重要性采样最怕的是：大多数权重很小，但偶尔出现极大权重。此时估计量可能看似稳定，实际上漏掉了罕见但极重要的区域。

尤其当目标分布 $q(\theta\mid y)$ 的尾部比近似分布 $g(\theta)$ 更厚时，重要性权重：

可能出现极端大值。此时估计方差可能很大，甚至不存在。

重要性分布的尾部必须足够厚。

如果 $g$ 的尾部比目标分布更薄，那么它很少抽到目标分布尾部的重要区域；但一旦抽到，权重会巨大。一个实用原则是：作为 proposal 的 $g$ 应该覆盖目标分布的所有重要区域，尾部宁可厚一点，不要薄。(即给予尾部一定的密度)

重要性重采样

重要性采样得到的是带权样本。如果我们想得到等权样本，可以做重要性重采样（importance resampling，也叫 sampling-importance resampling, SIR）。

步骤如下：

1. 从 $g(\theta)$ 中抽取候选样本：$\theta^1,\ldots,\theta^S$

2. 计算权重：

3. 按照概率：

从候选样本中重采样。

如果不放回抽取 $k<S$ 个样本，可以减少少数高权重点被重复抽中的问题；如果放回抽样，则更接近传统的重采样思想。

PSIS 是重要性采样的现代版本之一。

在BDA3 里提示了 Pareto-smoothed importance resampling。现在更常用的是 Pareto smoothed importance sampling（PSIS）：它会对最大的一批重要性权重拟合 generalized Pareto distribution，并平滑极端权重。Vehtari et al. 的 PSIS 论文已在 2024 年发表于 JMLR，并提出 Pareto $\hat k$ 作为权重尾部稳定性的诊断量。直观地说，$\hat k$ 越大，说明权重尾部越危险，重要性采样估计越不可靠。

Pareto-$\hat k$ diagnostic

重要性采样的问题是右尾太危险：大部分权重很小，少数权重极大，最后期望估计几乎被少数样本支配。PSIS 的想法是：不要只看权重的均值或者方差，而是直接检查最大一批权重的尾部形状。

在极值理论中，超过高阈值的尾部超额分布常常可以用 generalized Pareto distribution（GPD）近似。PSIS 正是对最大的若干个重要性权重拟合 GPD，并用估计出来的形状参数 $\hat k$ 判断尾部有多厚。

标准化的 GPD 可以写作：

对应密度为：

图中右侧使用 log-scale 展示尾部概率。可以看到，$k$ 越大，尾部下降越慢，也就越容易出现极端大的重要性权重。如果尾部近似满足 GPD 且 $k>0$，那么阶数为 $m$ 的矩存在大致要求 $m<1/k$，等价于 $k<1/m$。因为尾部密度满足类似幂律衰减：

因此可以得到几个非常重要的判断：

• $k<1$：权重分布的均值有限；
• $k<1/2$：权重分布的方差有限；
• $k$ 越接近或者超过 1，重要性采样对期望的估计越不稳定，甚至均值本身都可能无法良好定义。

这也是为什么 Pareto-$\hat k$ diagnostic 对所有依赖期望的估计都很重要。重要性采样、PSIS-LOO、prior/likelihood sensitivity 这类方法，最后通常都要计算某个加权期望。如果权重尾部太厚，那么“期望”这个目标本身就会变得非常脆弱。

早期实践中经常使用固定经验阈值：$\hat k<0.5$ 很好，$0.5<\hat k<0.7$ 还可以谨慎使用，$\hat k>0.7$ 通常要警惕。现在 loo 文档采用更细的样本量相关阈值：

$$\hat k < \min\left(1-\frac{1}{\log_{10}(S)},0.7\right)$$

时，PSIS 估计和相应 MCSE 通常被认为可靠。若 $\hat k$ 低于 0.7 但高于样本量相关阈值，增加有效样本量可能有帮助；若 $\hat k\ge 0.7$，估计常有较大偏差，单纯增加样本未必解决问题；若 $\hat k\ge 1$，目标权重分布估计为均值不有限，PSIS 估计和 MCSE 都不应被信任。 当 $\hat k$ 接近 0.5 时，更多有效抽样通常有助于降低 $\hat k$ 估计本身的不确定性，从而判断真实尾部是否确实处在 $k<0.5$ 的方差有限区间；但如果真实尾部很厚，增加样本只会更清楚地暴露问题，而不会从根本上修复 proposal 与目标分布不匹配。

在 PSIS-LOO 中的用途

PSIS-LOO 想用完整数据后验 $p(\theta\mid y)$ 的样本，近似留一后验 $p(\theta\mid y_{-i})$。对第 $i$ 个观测，原始重要性权重大致与下式成比例：

如果某个观测 $y_i$ 对后验影响很大，那么完整数据后验和留一后验差异会很大，$r_i^s$ 的尾部就会很厚，Pareto $\hat k_i$ 会升高。于是 $\hat k_i$ 不只是一个数值诊断，它也可以解释为观测影响力诊断：哪些点一旦被留出，模型后验会明显改变。

当 PSIS-LOO 中出现高 $\hat k$ 时，常见处理方式包括：

• 使用 moment matching 改善重要性采样；
• 对问题观测做 exact LOO 或者 K-fold cross-validation；
• 检查是否有异常值、高影响点、模型错设；
• 使用更稳健的 likelihood 或更合理的层级结构。

模拟样本数量计算

后验模拟通常用于报告：

• 后验中位数；
• 后验均值；
• 50% 或 95% 后验区间；
• 某个事件的后验概率；
• 后验预测分布；
• 后验预测检查量。

不同目标需要的模拟样本数不同。

后验均值的 Monte Carlo 误差

设标量参数 $\theta$ 的后验标准差为 $\sigma_\theta$，用 $S$ 个独立后验样本估计后验均值：

则：

实际中用样本标准差 $s_\theta$ 代替 $\sigma_\theta$：

如果我们关心的是参数本身的不确定性，后验标准差是 $\sigma_\theta$；如果再加上有限模拟导致的计算误差，总标准差近似为：

当 $S=100$ 时：

因此对于粗略后验总结，100 个独立样本已经可能够用。但是样本必须接近独立，且我们关心的是普通后验总结，而不是极端尾部或罕见事件。

后验概率的 Monte Carlo 误差

如果我们关心事件 $A$ 的后验概率：

用模拟估计：

那么：

所以：

当 $p$ 接近 0.5 时，$S=100$ 的标准误约为：

如果希望标准误约为 0.01，则需要：

分位数和罕见事件需要更多样本

后验分位数的 Monte Carlo 误差通常比均值更难稳定。设 $\xi_\alpha$ 是后验的 $\alpha$ 分位数，后验密度在该点为 $f(\xi_\alpha)$，则样本分位数的近似标准误为：

如果 $\alpha=0.025$ 或 $0.975$，并且尾部密度 $f(\xi_\alpha)$ 很小，那么分位数估计会比中位数更不稳定。

罕见事件也类似。如果某事件真实概率是 $0.0003$，那么 $S=200$ 个样本很可能一个都看不到，于是简单估计会给出 0。这不是事件不可能发生，而是模拟样本太少。

抽样数和有效样本量

BDA3 书本里为了说明概念，说很多例子中 100 到 2000 个模拟样本已经够表达主要推断。但在现代 MCMC 实践中，样本之间通常相关，所以应该看 effective sample size（ESS）、Monte Carlo standard error（MCSE）、$\hat R$、能量诊断、divergence 等指标，而不是只看总 draws 数。ArviZ 当前诊断接口就包括 ESS、rank-normalized split $\hat R$、MCSE 和 BFMI 等指标。

编程计算

Stan

Stan 的设计目标是用更高效的通用算法拟合连续参数贝叶斯模型。它使用 Hamiltonian Monte Carlo（HMC）及其自适应版本 No-U-Turn Sampler（NUTS），并通过自动微分计算 log density 的梯度。

现代概率编程生态

现在的贝叶斯计算生态已经比 BDA3 出版时丰富很多。常见工具包括：

• Stan：强调稳定的 HMC/NUTS、静态模型语言、C++ 后端；
• PyMC：Python 生态，常与 ArviZ 一起做诊断和可视化；
• NumPyro 和 BlackJAX：基于 JAX，强调自动微分、加速和可组合采样内核；
• Turing.jl：Julia 生态中的概率编程框架；
• INLA：适合某些 latent Gaussian models 的快速确定性近似；
• 变分推断工具：用于大数据或需要快速近似的场景。

贝叶斯计算的调试

贝叶斯计算需要调试。不是只有代码会出 bug，模型也会出 bug，数据也会出 bug，算法也可能被几何结构卡住。

假数据调试

假数据调试的基本步骤是：

1. 选择一个合理的真实参数值 $\theta$；
2. 如果是层级模型，先选择超参数，再从条件先验中生成下层参数；
3. 从数据生成模型中模拟假数据：

4. 用同一个模型对 $y^{\operatorname{fake}}$ 做后验推断；
5. 比较后验推断是否能覆盖或者恢复真实 $\theta$。

如果模型和计算都正确，那么后验区间应当有合适的覆盖性质。例如一个 50% 后验区间，在重复模拟下应该大约 50% 的时候覆盖真实参数。

这个性质可以这样理解。设 $C_\alpha(y)$ 是一个后验概率为 $\alpha$ 的区间，即：

如果 $\theta\sim p(\theta)$ 且 $y\sim p(y\mid \theta)$，则：

残差图和覆盖率检查

如果模型参数很多，比如层级模型中有许多 $\theta_j$，可以对每个参数比较：

然后画 $\text{error}_j$ 对后验均值的图。如果计算正确，误差会围绕 0，且没有明显的聚集。

更直接的方法是检查：真实值落在 50%、80%、95% 后验区间中的比例是否接近对应概率。

SBC 是假数据调试的现代形式。

Simulation-Based Calibration（SBC）为假数据检查系统化。基本思想是：从先验抽真实参数 $\theta^{(0)}$，再模拟数据 $y$，然后从 $p(\theta\mid y)$ 抽后验样本。如果计算正确，真实参数在后验样本中的 rank 应该近似均匀分布。Talts et al. 将 SBC 作为验证贝叶斯推断算法和模型实现的一般程序。2025 年还有 posterior SBC 的新进展，尝试在已观测数据条件下检查推断是否可靠；这类方法仍在发展中，但方向很重要。

这一节的计算方法可以整理成一张方法表格：