第 17 章 贝叶斯后验计算
第 16 章中的两个例子都有解析后验。但在真实模型中，后验分布往往没有简单闭式形式。本章介绍几种基础计算策略：网格近似、MAP 估计、Laplace 近似和后验模拟。它们与前面章节的数值积分、优化和 Monte Carlo 方法紧密相连。
17.1 后验计算的目标
给定后验分布
\[ p(\theta\mid y)\propto p(y\mid\theta)p(\theta), \]
我们通常关心以下量：
- 后验均值 \(E(\theta\mid y)\)；
- 后验标准差和可信区间；
- 后验众数，也称 MAP；
- 后验预测分布；
- 参数之间的相关性和不确定性传播。
这些量本质上都需要对后验分布做归纳。若无法解析计算，就需要数值方法。
17.2 网格近似
当参数只有一维或二维时，网格近似非常直观。选择一组候选点，在每个点计算未归一化后验密度，再把这些值标准化为离散概率。
下面考虑模拟企业月订单数。设
\[ y_i\mid\lambda\sim \operatorname{Poisson}(\lambda), \]
为了保证 \(\lambda>0\)，令 \(\eta=\log\lambda\)，并设先验
\[ \eta\sim N(1,1^2). \]
后验对数密度忽略常数后为
\[ \log p(\eta\mid y) = \left(\sum_{i=1}^n y_i\right)\eta -n\exp(\eta) -\frac{1}{2}(\eta-1)^2. \]
set.seed(2026)
y <- rpois(60, lambda = 3.5)
n <- length(y)
s <- sum(y)
log_post_eta <- function(eta) {
s * eta - n * exp(eta) - 0.5 * (eta - 1)^2
}
grid <- seq(0, 2.5, length.out = 1000)
lp <- log_post_eta(grid)
w <- exp(lp - max(lp))
prob <- w / sum(w)
eta_mean <- sum(grid * prob)
lambda_mean <- sum(exp(grid) * prob)
c(eta_mean = eta_mean,
lambda_mean = lambda_mean,
lambda_q025 = exp(grid[which(cumsum(prob) >= 0.025)[1]]),
lambda_q975 = exp(grid[which(cumsum(prob) >= 0.975)[1]]))
#> eta_mean lambda_mean lambda_q025 lambda_q975
#> 1.123 3.081 2.654 3.539网格法的优点是透明，缺点是维度升高后迅速不可行。如果每个参数取 1000 个网格点，三个参数就需要\(10^9\) 个点。
17.3 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)\}. \]
它可以用第 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 <- fit_map$par
lambda_map <- exp(eta_map)
c(eta_map = eta_map, lambda_map = lambda_map)
#> eta_map lambda_map
#> 1.125 3.081MAP 是一个点估计。它计算方便，也常作为更复杂算法的初始值。但只报告 MAP 会丢失后验不确定性。
17.4 Laplace 后验近似
在 MAP 附近，用二次函数近似对数后验，就得到正态近似。若
\[ H = - \left. \frac{\partial^2 \log p(\theta\mid y)} {\partial\theta^2} \right|_{\theta=\hat\theta_{\text{MAP}}}, \]
则
\[ \theta\mid y \approx N(\hat\theta_{\text{MAP}},H^{-1}). \]
在 Poisson 例子中，\(\eta\) 的二阶导数为
\[ \frac{\partial^2 \log p(\eta\mid y)}{\partial\eta^2} = -n\exp(\eta)-1. \]
hessian_neg <- n * exp(eta_map) + 1
eta_sd_laplace <- sqrt(1 / hessian_neg)
set.seed(1)
eta_draw <- rnorm(5000, mean = eta_map, sd = eta_sd_laplace)
lambda_draw <- exp(eta_draw)
c(lambda_mean_laplace = mean(lambda_draw),
lambda_q025 = quantile(lambda_draw, 0.025),
lambda_q975 = quantile(lambda_draw, 0.975))
#> lambda_mean_laplace lambda_q025.2.5% lambda_q975.97.5%
#> 3.089 2.650 3.573Laplace 近似把后验计算转化为优化和 Hessian 计算。它在样本量较大、后验近似单峰对称时效果较好。若后验偏斜、多峰或靠近边界，则需要更谨慎。
17.5 后验模拟与归纳
一旦获得后验样本
\[ \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). \]
c(post_mean = mean(lambda_draw),
post_sd = sd(lambda_draw),
prob_lambda_gt_4 = mean(lambda_draw > 4))
#> post_mean post_sd prob_lambda_gt_4
#> 3.0893 0.2326 0.0004这种“用样本归纳后验”的思想会贯穿后续贝叶斯回归、后验预测和现代贝叶斯软件。