第 17 章 贝叶斯后验计算
第 16 章中的若干例子有解析后验，可以直接写出 Beta、正态或逆卡方分布。但在真实经济统计模型中，后验分布往往没有简单闭式形式。本章介绍几种基础计算策略：网格近似、MAP 估计、Laplace 近似、重要性抽样和后验模拟。它们与前面章节的数值积分、优化和 Monte Carlo 方法紧密相连。
本章的学习重点有三个。第一，理解贝叶斯计算本质上是在计算后验分布的归一化、积分和概率；第二，掌握一维或低维后验的几种可解释近似方法；第三，能够用后验样本计算均值、区间、概率和简单预测。
17.1 后验计算的目标
给定参数 \(\theta\) 和数据 \(y\)，贝叶斯公式为
\[ p(\theta\mid y) = \frac{p(y\mid\theta)p(\theta)} {\int p(y\mid\theta)p(\theta)\,d\theta}. \]
其中 \(p(y\mid\theta)\) 是似然函数，\(p(\theta)\) 是先验分布，\(p(\theta\mid y)\) 是后验分布。分母
\[ p(y)=\int p(y\mid\theta)p(\theta)\,d\theta \]
称为边际似然或归一化常数。它保证后验密度积分为 1，但在许多模型中难以直接计算。因此，实际计算常从未归一化后验开始：
\[ \tilde p(\theta\mid y)=p(y\mid\theta)p(\theta), \]
并利用
\[ p(\theta\mid y)\propto \tilde p(\theta\mid y). \]
后验计算通常关心以下量：
- 后验均值 \(E(\theta\mid y)\)；
- 后验标准差和可信区间；
- 后验众数，即 MAP；
- 后验概率，例如 \(P(\theta>c\mid y)\)；
- 后验预测分布，例如未来一个月订单数或违约数的分布。
这些量本质上都需要对后验分布做归纳。若无法解析计算，就需要数值方法。
17.2 对数后验与数值稳定
后验密度通常是很多概率或密度的乘积。直接相乘容易下溢，因此贝叶斯计算一般使用对数后验：
\[ \log p(\theta\mid y) = \log p(y\mid\theta)+\log p(\theta)-\log p(y). \]
若只需要比较不同 \(\theta\) 下后验密度的相对大小，可以忽略不依赖 \(\theta\) 的常数，写成
\[ \ell_p(\theta) = \log p(y\mid\theta)+\log p(\theta)+C, \]
其中 \(C\) 是与 \(\theta\) 无关的常数。下面的计算都会先写出这种“差一个常数”的对数后验。
为了把一组对数权重变成概率，需要使用第 3 章提到的稳定化技巧。设\(a_1,\ldots,a_K\) 是对数权重，直接计算 \(\exp(a_k)\) 可能上溢或下溢。可以先减去最大值：
\[ w_k=\frac{\exp(a_k-\max_j a_j)} {\sum_{\ell=1}^K \exp(a_\ell-\max_j a_j)}. \]
下面写两个辅助函数：一个把对数权重标准化，另一个计算加权分位数。
17.3 网格近似
当参数只有一维或二维时，网格近似非常直观。选择一组候选点，在每个点计算未归一化后验密度，再把这些值标准化为离散概率。网格法的价值不仅在于得到数值结果，也在于帮助我们看清后验形状。
下面考虑一个模拟的企业月订单数例子。设
\[ y_i\mid\lambda\sim \operatorname{Poisson}(\lambda), \quad i=1,\ldots,n, \]
其中 \(y_i\) 是第 \(i\) 个月的订单数，\(\lambda\) 是每月平均订单数。为了保证 \(\lambda>0\)，令
\[ \eta=\log\lambda. \]
设先验为
\[ \eta\sim N(m_0,\tau_0^2), \]
其中 \(m_0=1\)，\(\tau_0=1\)。Poisson 似然为
\[ p(y\mid\lambda) = \prod_{i=1}^n \frac{\exp(-\lambda)\lambda^{y_i}}{y_i!}. \]
代入 \(\lambda=\exp(\eta)\)，并忽略与 \(\eta\) 无关的常数，后验对数密度为
\[ \ell_p(\eta) = \left(\sum_{i=1}^n y_i\right)\eta -n\exp(\eta) -\frac{1}{2}\left(\frac{\eta-m_0}{\tau_0}\right)^2. \]
这里 \(\sum_i y_i\) 是订单数总和，\(n\) 是月份数，最后一项来自正态先验。
set.seed(2026)
y <- rpois(60, lambda = 3.5)
n <- length(y)
s <- sum(y)
m0 <- 1
tau0 <- 1
log_post_eta <- function(eta) {
s * eta - n * exp(eta) - 0.5 * ((eta - m0) / tau0)^2
}
grid <- seq(-0.2, 2.6, length.out = 1500)
lp <- log_post_eta(grid)
prob <- normalize_log_weights(lp)
lambda_grid <- exp(grid)
eta_mean <- sum(grid * prob)
eta_sd <- sqrt(sum((grid - eta_mean)^2 * prob))
eta_q <- weighted_quantile(grid, prob, c(0.025, 0.5, 0.975))
lambda_mean <- sum(lambda_grid * prob)
lambda_q <- weighted_quantile(lambda_grid, prob, c(0.025, 0.5, 0.975))
prob_lambda_gt_4_grid <- sum(prob[lambda_grid > 4])
grid_summary <- data.frame(
quantity = c("eta_mean", "eta_sd", "lambda_mean",
"lambda_q025", "lambda_q50", "lambda_q975",
"Pr(lambda > 4)"),
value = c(eta_mean, eta_sd, lambda_mean,
lambda_q[1], lambda_q[2], lambda_q[3],
prob_lambda_gt_4_grid)
)
knitr::kable(grid_summary, row.names = FALSE, digits = 4)| quantity | value |
|---|---|
| eta_mean | 1.1227 |
| eta_sd | 0.0734 |
| lambda_mean | 3.0813 |
| lambda_q025 | 2.6559 |
| lambda_q50 | 3.0782 |
| lambda_q975 | 3.5411 |
| Pr(lambda > 4) | 0.0001 |
这个结果把后验分布压缩成若干可解释摘要。\(\lambda\) 的后验均值是每月平均订单数的后验估计，lambda_q025 和 lambda_q975 构成近似 95% 可信区间，Pr(lambda > 4) 表示在给定模型和数据后，平均订单数超过 4 的后验概率。
网格法需要选择范围和密度。范围太窄会漏掉后验尾部，网格太粗会影响分位数和概率。下面比较不同网格下 \(\lambda\) 后验均值的变化。
grid_specs <- data.frame(
lower = c(0.4, -0.2, -0.8),
upper = c(2.0, 2.6, 3.2),
points = c(300, 1500, 3000)
)
grid_result <- data.frame()
for (i in seq_len(nrow(grid_specs))) {
g <- seq(grid_specs$lower[i], grid_specs$upper[i],
length.out = grid_specs$points[i])
p <- normalize_log_weights(log_post_eta(g))
grid_result <- rbind(
grid_result,
data.frame(
lower = grid_specs$lower[i],
upper = grid_specs$upper[i],
points = grid_specs$points[i],
lambda_mean = sum(exp(g) * p)
)
)
}
knitr::kable(grid_result, row.names = FALSE, digits = 4)| lower | upper | points | lambda_mean |
|---|---|---|---|
| 0.4 | 2.0 | 300 | 3.081 |
| -0.2 | 2.6 | 1500 | 3.081 |
| -0.8 | 3.2 | 3000 | 3.081 |
若结果对网格范围非常敏感，通常说明网格没有覆盖主要后验质量，或参数变换不合适。网格法的缺点是维度升高后迅速不可行。如果每个参数取 1000 个网格点，三个参数就需要 \(10^9\) 个点。
17.4 MAP 估计
MAP（maximum a posteriori）估计是后验众数：
\[ \hat\theta_{\text{MAP}} = \arg\max_\theta p(\theta\mid y) = \arg\max_\theta \{\log p(y\mid\theta)+\log p(\theta)\}. \]
其中 \(\hat\theta_{\text{MAP}}\) 是使后验密度最大的参数值。它可以用第 6 章的优化算法计算。对于上面的 Poisson 模型，我们最大化 log_post_eta()，等价于最小化负对数后验。
fit_map <- optim(par = log(mean(y)),
fn = function(e) -log_post_eta(e),
method = "BFGS",
hessian = TRUE)
eta_map <- as.numeric(fit_map$par)
lambda_map <- exp(eta_map)
map_out <- data.frame(
quantity = c("convergence", "eta_map", "lambda_map"),
value = c(fit_map$convergence, eta_map, lambda_map)
)
knitr::kable(map_out, row.names = FALSE, digits = 4)| quantity | value |
|---|---|
| convergence | 0.000 |
| eta_map | 1.125 |
| lambda_map | 3.081 |
convergence = 0 通常表示 optim() 正常收敛。MAP 是一个点估计，计算方便，也常作为更复杂算法的初始值。但只报告 MAP 会丢失后验不确定性。贝叶斯分析通常还需要后验标准差、区间或后验概率。
17.5 Laplace 后验近似
在 MAP 附近，用二次函数近似对数后验，就得到正态近似。若
\[ H = - \left. \frac{\partial^2 \log p(\theta\mid y)} {\partial\theta\partial\theta^\top} \right|_{\theta=\hat\theta_{\text{MAP}}}, \]
则
\[ \theta\mid y \approx N(\hat\theta_{\text{MAP}},H^{-1}). \]
这里 \(H\) 是负 Hessian 矩阵。在一维情形下，它就是负二阶导数。对于 Poisson 例子，
\[ \frac{\partial^2 \ell_p(\eta)}{\partial\eta^2} = -n\exp(\eta)-\frac{1}{\tau_0^2}. \]
因此
\[ H=n\exp(\hat\eta_{\text{MAP}})+\frac{1}{\tau_0^2}. \]
hessian_neg <- n * exp(eta_map) + 1 / tau0^2
eta_sd_laplace <- sqrt(1 / hessian_neg)
set.seed(1)
S <- 10000
eta_draw_laplace <- rnorm(S, mean = eta_map, sd = eta_sd_laplace)
lambda_draw_laplace <- exp(eta_draw_laplace)
laplace_summary <- data.frame(
method = c("grid", "Laplace"),
lambda_mean = c(lambda_mean, mean(lambda_draw_laplace)),
lambda_q025 = c(lambda_q[1], unname(quantile(lambda_draw_laplace, 0.025))),
lambda_q975 = c(lambda_q[3], unname(quantile(lambda_draw_laplace, 0.975))),
prob_lambda_gt_4 = c(prob_lambda_gt_4_grid,
mean(lambda_draw_laplace > 4))
)
knitr::kable(laplace_summary, row.names = FALSE, digits = 4)| method | lambda_mean | lambda_q025 | lambda_q975 | prob_lambda_gt_4 |
|---|---|---|---|---|
| grid | 3.081 | 2.656 | 3.541 | 1e-04 |
| Laplace | 3.088 | 2.656 | 3.565 | 3e-04 |
Laplace 近似把后验计算转化为优化和 Hessian 计算。它在样本量较大、后验近似单峰对称时效果较好。若后验偏斜、多峰或靠近边界，则需要更谨慎。注意本例在 \(\eta=\log\lambda\) 尺度上做正态近似，再变换回 \(\lambda\)；这比直接在正数参数 \(\lambda\) 上近似更自然。
17.6 重要性抽样近似后验
第 10 章介绍过重要性抽样。后验计算中也经常遇到类似问题：我们知道未归一化后验\(\tilde p(\theta\mid y)\)，但不知道归一化常数。若从建议分布 \(q(\theta)\) 中抽样\(\theta^{(1)},\ldots,\theta^{(S)}\)，可以使用权重
\[ w_s\propto \frac{\tilde p(\theta^{(s)}\mid y)} {q(\theta^{(s)})}. \]
归一化权重为
\[ \bar w_s=\frac{w_s}{\sum_{r=1}^S w_r}. \]
后验期望可近似为
\[ E\{h(\theta)\mid y\} \approx \sum_{s=1}^S \bar w_s h(\theta^{(s)}). \]
下面用以 MAP 为中心的正态分布作为建议分布，对 Poisson 后验做重要性抽样。
set.seed(17)
S <- 10000
proposal_sd <- 1.5 * eta_sd_laplace
eta_prop <- rnorm(S, mean = eta_map, sd = proposal_sd)
log_w_raw <- log_post_eta(eta_prop) -
dnorm(eta_prop, mean = eta_map, sd = proposal_sd, log = TRUE)
w_norm <- normalize_log_weights(log_w_raw)
lambda_prop <- exp(eta_prop)
ess <- 1 / sum(w_norm^2)
is_summary <- data.frame(
quantity = c("lambda_mean", "lambda_q025", "lambda_q975",
"Pr(lambda > 4)", "ESS", "max_weight"),
value = c(sum(w_norm * lambda_prop),
weighted_quantile(lambda_prop, w_norm, 0.025),
weighted_quantile(lambda_prop, w_norm, 0.975),
sum(w_norm[lambda_prop > 4]),
ess,
max(w_norm))
)
knitr::kable(is_summary, row.names = FALSE, digits = 4)| quantity | value |
|---|---|
| lambda_mean | 3.0815 |
| lambda_q025 | 2.6574 |
| lambda_q975 | 3.5386 |
| Pr(lambda > 4) | 0.0001 |
| ESS | 8337.6791 |
| max_weight | 0.0001 |
有效样本量
\[ \operatorname{ESS}=\frac{1}{\sum_{s=1}^S \bar w_s^2} \]
用于衡量权重是否集中在少数样本上。若 ESS 远小于 \(S\)，说明建议分布与后验不匹配，重要性抽样估计可能不稳定。后续 MCMC 方法可以看作在更复杂情形下生成后验样本的一类通用方法。
17.7 后验模拟与归纳
一旦获得后验样本
\[ \theta^{(1)},\ldots,\theta^{(S)}, \]
后验均值、标准差和区间都可以用样本统计量近似：
\[ E(\theta\mid y)\approx \frac{1}{S}\sum_{s=1}^S \theta^{(s)}. \]
后验概率也可以直接用比例估计。例如
\[ P(\lambda>4\mid y) \approx \frac{1}{S}\sum_{s=1}^S I(\lambda^{(s)}>4). \]
如果样本带有重要性权重，则把普通平均改成加权平均。为了得到普通后验样本，也可以按归一化权重进行重抽样。
set.seed(18)
id <- sample(seq_len(S), size = 5000, replace = TRUE, prob = w_norm)
eta_draw <- eta_prop[id]
lambda_draw <- exp(eta_draw)
posterior_summary <- c(
post_mean = mean(lambda_draw),
post_sd = sd(lambda_draw),
q025 = unname(quantile(lambda_draw, 0.025)),
q975 = unname(quantile(lambda_draw, 0.975)),
prob_lambda_gt_4 = mean(lambda_draw > 4)
)
posterior_summary
#> post_mean post_sd q025 q975
#> 3.0839 0.2254 2.6652 3.5401
#> prob_lambda_gt_4
#> 0.0004这些计算看起来简单，却是贝叶斯分析中最常用的后验归纳方式。复杂模型中，后验样本可能来自 MCMC 或Stan 等工具，但总结方式仍然是均值、分位数、概率和图形诊断。
后验样本还可以自然传播到预测。未来一个月订单数 \(\tilde y\) 的后验预测分布为
\[ p(\tilde y\mid y) = \int p(\tilde y\mid\lambda)p(\lambda\mid y)\,d\lambda. \]
用后验样本近似时，只需先抽 \(\lambda^{(s)}\)，再抽\(\tilde y^{(s)}\sim\operatorname{Poisson}(\lambda^{(s)})\)。
set.seed(19)
y_new <- rpois(length(lambda_draw), lambda = lambda_draw)
predictive_summary <- c(
pred_mean = mean(y_new),
pred_q025 = unname(quantile(y_new, 0.025)),
pred_q975 = unname(quantile(y_new, 0.975)),
prob_orders_ge_6 = mean(y_new >= 6)
)
predictive_summary
#> pred_mean pred_q025 pred_q975 prob_orders_ge_6
#> 3.0704 0.0000 7.0000 0.0922这个预测分布同时包含参数不确定性和未来订单数本身的随机性。第 19章会进一步讨论后验预测检验和模型比较。
17.8 方法选择
不同后验计算方法适用于不同场景。下面给出本科课程中可以采用的简单判断。
| 方法 | 适用场景 | 主要输出 |
|---|---|---|
| 解析后验 | 共轭模型或简单模型 | 完整后验摘要 |
| 网格近似 | 一维或二维参数 | 完整但低维的离散近似 |
| MAP | 需要快速点估计或算法初值 | 后验众数 |
| Laplace 近似 | 单峰、近似对称的后验 | 近似均值、方差、区间 |
| 重要性抽样 | 可找到合适建议分布的低维问题 | 加权后验摘要 |
| MCMC/现代软件 | 高维、非共轭或复杂层次模型 | 后验样本与诊断 |