第 19 章 贝叶斯模型比较与后验预测
贝叶斯分析的目标不只是估计参数。很多经济统计问题更关心预测和比较：明天某商圈夜间外卖订单量会不会超过运力上限？某个消费模型能否再现数据中的高峰需求？两个模型预测能力相近时，应该选择哪个？本章介绍后验预测、后验预测检验和贝叶斯模型比较的基本思想。
第 18 章重点讨论回归系数的后验分布。本章把视角从“参数是多少”转向“模型能生成怎样的数据”。这是贝叶斯建模中非常重要的一步：一个模型即使有漂亮的系数表，也可能不能合理预测未来，或不能再现观测数据中的关键特征。
本章的学习重点有三个。第一，理解后验预测分布如何同时包含参数不确定性和未来随机性；第二，掌握后验预测检验的基本计算方法；第三，了解 Bayes 因子、WAIC 和样本外预测比较各自回答什么问题。
19.1 后验预测分布
设 \(\tilde y\) 表示未来观测值，\(\theta\) 表示模型参数。贝叶斯后验预测分布为
\[ p(\tilde y\mid y) = \int p(\tilde y\mid\theta)p(\theta\mid y)\,d\theta. \]
其中 \(p(\theta\mid y)\) 描述参数不确定性，\(p(\tilde y\mid\theta)\) 描述给定参数后的未来随机性。这两个不确定性在预测中都会出现。若忽略参数不确定性，预测区间通常会偏窄；若只报告参数可信区间，又不能回答未来观测值的实际波动范围。
如果有后验样本 \(\theta^{(1)},\ldots,\theta^{(S)}\)，后验预测可以按以下步骤模拟：
- 从后验样本中取 \(\theta^{(s)}\)；
- 从 \(p(\tilde y\mid\theta^{(s)})\) 中生成 \(\tilde y^{(s)}\)；
- 用 \(\tilde y^{(1)},\ldots,\tilde y^{(S)}\) 汇总预测均值、预测区间和事件概率。
这个算法与第 17 章中的“用样本归纳后验”完全一致，只是这里归纳的对象从参数变成了未来数据。
19.2 案例：夜间外卖订单预测
下面使用一个模拟案例。某平台想预测一个城市商圈未来夜间外卖订单量，以便安排骑手运力和补贴预算。研究者收集了若干天的夜间订单量、气温、是否下雨、是否周末以及是否发放平台补贴。数据是模拟数据，不对应真实平台记录。
这个案例有两个有趣之处。第一，需求与气温可能不是简单线性关系：天气太冷或太热都会降低夜间订单。第二，平台不仅关心平均预测，还关心“订单量超过 700 单”的概率，因为这会触发额外骑手调度。
set.seed(1901)
n_days <- 220
day_id <- seq_len(n_days)
temp <- runif(n_days, min = 4, max = 34)
temp_c <- temp - 22
rain <- rbinom(n_days, size = 1, prob = plogis(-0.5 + 0.03 * temp))
weekend <- as.integer(day_id %% 7 %in% c(0, 6))
subsidy <- rbinom(n_days, size = 1, prob = plogis(-1.1 + 0.7 * weekend))
mu_true <- 520 + 7 * temp_c - 0.75 * temp_c^2 +
65 * rain + 55 * weekend + 45 * subsidy + 28 * rain * weekend
sigma_true <- 45 + 22 * rain
orders <- round(rnorm(n_days, mean = mu_true, sd = sigma_true))
delivery_dat <- data.frame(
orders = orders,
temp_c = temp_c,
rain = rain,
weekend = weekend,
subsidy = subsidy
)
knitr::kable(head(delivery_dat), row.names = FALSE, digits = 2)| orders | temp_c | rain | weekend | subsidy |
|---|---|---|---|---|
| 514 | 9.35 | 1 | 0 | 0 |
| 220 | -16.45 | 1 | 0 | 0 |
| 488 | -6.90 | 0 | 0 | 0 |
| 495 | -4.71 | 1 | 0 | 0 |
| 591 | 11.42 | 1 | 0 | 0 |
| 626 | 0.58 | 0 | 1 | 1 |
这里 orders 是夜间订单量，temp_c 是相对于 22 摄氏度的温度偏离，rain、weekend 和 subsidy分别表示是否下雨、是否周末和是否发放补贴。使用中心化温度 temp_c 是为了让截距更容易解释：当气温接近 22 摄氏度且其他变量为 0 时，截距近似表示基准订单量。
为保持本章重点集中在预测和比较上，下面继续使用第 18 章的正态线性回归近似。误差标准差 \(\sigma\) 先用训练模型的残差标准差作为插入估计；更完整的处理可以像上一章那样给\(\sigma^2\) 指定先验并一并抽样。
rmvnorm_chol <- function(n, mean, Sigma) {
z <- matrix(rnorm(n * length(mean)), nrow = n)
sweep(z %*% chol(Sigma), 2, mean, FUN = "+")
}
log_mean_exp <- function(x) {
m <- max(x)
m + log(mean(exp(x - m)))
}
fit_bayes_lm <- function(formula, data, S = 4000, sigma = NULL) {
X <- model.matrix(formula, data = data)
y <- model.response(model.frame(formula, data = data))
if (is.null(sigma)) {
sigma <- summary(lm(formula, data = data))$sigma
}
p <- ncol(X)
b0 <- rep(0, p)
b0[1] <- mean(y)
prior_sd <- rep(120, p)
prior_sd[1] <- 500
V0 <- diag(prior_sd^2)
Vn <- solve(crossprod(X) / sigma^2 + solve(V0))
bn <- Vn %*% (crossprod(X, y) / sigma^2 + solve(V0, b0))
beta_draw <- rmvnorm_chol(S, as.vector(bn), Vn)
colnames(beta_draw) <- colnames(X)
list(
formula = formula,
terms = delete.response(terms(formula)),
response = as.character(formula[[2]]),
X = X,
y = y,
sigma = sigma,
beta_draw = beta_draw
)
}
posterior_predict <- function(fit, newdata, include_noise = TRUE) {
X_new <- model.matrix(fit$terms, data = newdata)
mu <- fit$beta_draw %*% t(X_new)
if (!include_noise) {
return(mu)
}
matrix(
rnorm(length(mu), mean = as.vector(mu), sd = fit$sigma),
nrow = nrow(mu)
)
}
pointwise_loglik <- function(fit, data) {
X_eval <- model.matrix(fit$terms, data = data)
y_eval <- data[[fit$response]]
mu <- fit$beta_draw %*% t(X_eval)
matrix(
dnorm(rep(y_eval, each = nrow(mu)),
mean = as.vector(mu),
sd = fit$sigma,
log = TRUE),
nrow = nrow(mu)
)
}19.3 预测均值与未来观测
先拟合一个包含二次温度项的模型：
\[ y_i = \beta_0 +\beta_1\operatorname{temp}_{c,i} +\beta_2\operatorname{temp}_{c,i}^2 +\beta_3\operatorname{rain}_i +\beta_4\operatorname{weekend}_i +\beta_5\operatorname{subsidy}_i +\varepsilon_i. \]
其中 \(y_i\) 是第 \(i\) 天夜间订单量，\(\operatorname{temp}_{c,i}\) 是中心化气温，\(\varepsilon_i\) 是误差项。二次项允许订单量随气温先升后降。
set.seed(1902)
fit_quad <- fit_bayes_lm(
orders ~ temp_c + I(temp_c^2) + rain + weekend + subsidy,
data = delivery_dat,
S = 5000
)
coef_summary <- data.frame(
term = colnames(fit_quad$X),
mean = colMeans(fit_quad$beta_draw),
q025 = apply(fit_quad$beta_draw, 2, quantile, 0.025),
q975 = apply(fit_quad$beta_draw, 2, quantile, 0.975)
)
knitr::kable(coef_summary, row.names = FALSE, digits = 3)| term | mean | q025 | q975 |
|---|---|---|---|
| (Intercept) | 521.074 | 504.622 | 537.774 |
| temp_c | 6.662 | 5.504 | 7.796 |
| I(temp_c^2) | -0.806 | -0.921 | -0.691 |
| rain | 75.879 | 59.831 | 91.522 |
| weekend | 60.712 | 43.059 | 77.962 |
| subsidy | 47.377 | 29.872 | 64.306 |
现在考虑一个未来情景：明天夜间气温比 22 摄氏度高 4 度，下雨，是周末，并且平台发放补贴。我们希望预测订单量，并计算超过 700 单的概率。
future_day <- data.frame(
temp_c = 4,
rain = 1,
weekend = 1,
subsidy = 1
)
mu_future <- posterior_predict(fit_quad, future_day, include_noise = FALSE)[, 1]
y_future <- posterior_predict(fit_quad, future_day, include_noise = TRUE)[, 1]
future_summary <- data.frame(
quantity = c("mean of expected orders",
"q05 of expected orders",
"q95 of expected orders",
"mean of future orders",
"q05 of future orders",
"q95 of future orders",
"Pr(future orders > 700)"),
value = c(mean(mu_future),
quantile(mu_future, 0.05),
quantile(mu_future, 0.95),
mean(y_future),
quantile(y_future, 0.05),
quantile(y_future, 0.95),
mean(y_future > 700))
)
knitr::kable(future_summary, row.names = FALSE, digits = 3)| quantity | value |
|---|---|
| mean of expected orders | 718.79 |
| q05 of expected orders | 702.14 |
| q95 of expected orders | 734.93 |
| mean of future orders | 719.53 |
| q05 of future orders | 622.28 |
| q95 of future orders | 817.60 |
| Pr(future orders > 700) | 0.63 |
这里有两个区间需要区分。expected orders 是条件均值\(E(\tilde y\mid x,\mathbf y)\) 的不确定性，主要来自参数后验不确定性；future orders 是未来实际订单量的预测分布，还包含当天随机波动。因此后者的区间通常更宽。在运营决策中，如果要安排骑手，应更关注未来实际订单量的预测分布。
19.4 后验预测检验
后验预测检验的思想是：如果模型合理，那么从模型生成的复制数据应当在关键特征上类似于观测数据。设 \(T(y)\) 是一个检验统计量，例如均值、标准差、高峰比例或组间差异。对每个后验样本生成复制数据\(y^{\operatorname{rep}}\)，比较
\[ T(y^{\operatorname{rep}}) \]
和
\[ T(y). \]
常用的后验预测 p 值为
\[ p_B = P\{T(y^{\operatorname{rep}})\ge T(y)\mid y\}. \]
这里 \(p_B\) 是模型检查指标，不是频率学派假设检验中的 p 值。若它非常接近 0 或 1，说明观测数据的某个特征位于模型复制分布的极端位置，模型可能遗漏了重要结构。
下面检查二次模型能否再现四个数据特征：平均订单量、订单量标准差、高峰天数比例，以及雨天与非雨天订单波动差异。
set.seed(1903)
y_rep <- posterior_predict(fit_quad, delivery_dat, include_noise = TRUE)
y_obs <- delivery_dat$orders
rain_idx <- delivery_dat$rain == 1
dry_idx <- delivery_dat$rain == 0
T_obs <- c(
mean_orders = mean(y_obs),
sd_orders = sd(y_obs),
high_share = mean(y_obs > 700),
rainy_minus_dry_sd = sd(y_obs[rain_idx]) - sd(y_obs[dry_idx])
)
T_rep <- data.frame(
mean_orders = rowMeans(y_rep),
sd_orders = apply(y_rep, 1, sd),
high_share = rowMeans(y_rep > 700),
rainy_minus_dry_sd = apply(y_rep[, rain_idx], 1, sd) -
apply(y_rep[, dry_idx], 1, sd)
)
ppc_summary <- data.frame(
statistic = names(T_obs),
observed = as.numeric(T_obs),
replicated_mean = colMeans(T_rep),
posterior_predictive_p = sapply(names(T_obs), function(nm) {
mean(T_rep[[nm]] >= T_obs[nm])
})
)
knitr::kable(ppc_summary, row.names = FALSE, digits = 3)| statistic | observed | replicated_mean | posterior_predictive_p |
|---|---|---|---|
| mean_orders | 496.373 | 496.360 | 0.493 |
| sd_orders | 155.462 | 155.874 | 0.531 |
| high_share | 0.077 | 0.065 | 0.259 |
| rainy_minus_dry_sd | 23.375 | 0.913 | 0.002 |
如果均值和高峰比例能够被复制，但雨天与非雨天的波动差异难以复制，就提示当前模型可能仍有缺陷。在本案例中，原因可能是误差方差随天气变化，而我们使用了同方差正态线性回归。后验预测检验的价值正在这里：它不只是给出“模型好坏”的一个数字，而是告诉我们模型在哪些方面可能不够好。
19.5 贝叶斯模型比较的目标
模型比较首先要明确目标。不同准则回答的问题并不相同：
| 方法 | 主要问题 | 适合场景 | 注意事项 |
|---|---|---|---|
| Bayes 因子 | 数据相对支持哪个模型 | 少数候选模型、先验经过认真设定 | 对先验敏感，边际似然可能难算 |
| WAIC | 哪个模型样本内预测能力更好，并惩罚复杂度 | 有后验样本和逐点似然 | 是近似准则，样本量小或后验异常时需谨慎 |
| 留一交叉验证 | 每次留出一个观测时的预测表现 | 预测导向模型比较 | 计算量较大，实际常用 PSIS-LOO |
| 训练测试划分 | 样本外预测误差 | 预测是核心目标且样本量较充足 | 结果受划分影响，可重复划分 |
经典贝叶斯模型比较基于边际似然。给定模型 \(M\)，
\[ p(y\mid M) = \int p(y\mid\theta,M)p(\theta\mid M)\,d\theta. \]
两个模型 \(M_1\) 和 \(M_2\) 的 Bayes 因子为
\[ BF_{12} = \frac{p(y\mid M_1)}{p(y\mid M_2)}. \]
如果模型先验概率相同，\(BF_{12}>1\) 表示数据相对更支持 \(M_1\)。但边际似然会把先验分布也积分进去，因此 Bayes 因子对先验宽窄可能敏感。对于本科阶段的应用课程，更重要的是理解它的含义，而不必把它作为唯一模型比较工具。
19.6 WAIC：基于预测的贝叶斯比较
预测导向的贝叶斯模型比较常使用 WAIC。它利用逐点对数预测密度衡量拟合，同时用后验样本中的变异估计有效复杂度。给定后验样本 \(\theta^{(1)},\ldots,\theta^{(S)}\)，定义
\[ \operatorname{lppd} = \sum_{i=1}^n \log\left[ \frac{1}{S}\sum_{s=1}^S p(y_i\mid\theta^{(s)}) \right]. \]
这里 \(p(y_i\mid\theta^{(s)})\) 是第 \(s\) 个后验样本下第 \(i\) 个观测的似然值。有效参数个数近似为
\[ p_{\operatorname{WAIC}} = \sum_{i=1}^n \operatorname{Var}_{s}\{\log p(y_i\mid\theta^{(s)})\}. \]
于是
\[ \operatorname{WAIC} = -2(\operatorname{lppd}-p_{\operatorname{WAIC}}). \]
WAIC 越小，通常表示预测表现越好。实际比较时常看 \(\Delta\operatorname{WAIC}\)，即某模型 WAIC 与最小 WAIC 的差。
下面比较三个候选模型：
\[ M_1:\ y_i=\beta_0+\beta_1\operatorname{temp}_{c,i} +\beta_2\operatorname{rain}_i +\beta_3\operatorname{weekend}_i +\beta_4\operatorname{subsidy}_i+\varepsilon_i, \]
\[ M_2:\ M_1+\beta_5\operatorname{temp}_{c,i}^2, \]
\[ M_3:\ M_2+\beta_6(\operatorname{rain}_i\times\operatorname{weekend}_i). \]
waic_from_loglik <- function(loglik) {
lppd_i <- apply(loglik, 2, log_mean_exp)
p_waic_i <- apply(loglik, 2, var)
elpd_waic <- sum(lppd_i - p_waic_i)
data.frame(
elpd_waic = elpd_waic,
p_waic = sum(p_waic_i),
waic = -2 * elpd_waic
)
}set.seed(1904)
model_formulas <- list(
linear = orders ~ temp_c + rain + weekend + subsidy,
quadratic = orders ~ temp_c + I(temp_c^2) + rain + weekend + subsidy,
interaction = orders ~ temp_c + I(temp_c^2) + rain * weekend + subsidy
)
model_fits <- lapply(model_formulas, fit_bayes_lm,
data = delivery_dat, S = 4000)
waic_tab <- do.call(
rbind,
lapply(names(model_fits), function(nm) {
out <- waic_from_loglik(pointwise_loglik(model_fits[[nm]], delivery_dat))
data.frame(model = nm, out)
})
)
waic_tab$delta_waic <- waic_tab$waic - min(waic_tab$waic)
raw_weight <- exp(-0.5 * waic_tab$delta_waic)
waic_tab$waic_weight <- raw_weight / sum(raw_weight)
knitr::kable(waic_tab, row.names = FALSE, digits = 3)| model | elpd_waic | p_waic | waic | delta_waic | waic_weight |
|---|---|---|---|---|---|
| linear | -1279 | 5.217 | 2559 | 139.40 | 0.000 |
| quadratic | -1212 | 5.741 | 2424 | 4.83 | 0.082 |
| interaction | -1210 | 6.673 | 2419 | 0.00 | 0.918 |
waic_weight 是由 WAIC 差异构造的相对权重，可用于近似表达模型预测支持度。它不是严格的模型后验概率，但在预测比较中常有参考价值。若两个模型 WAIC 很接近，应避免过度解读微小差异。
19.7 样本外预测比较
WAIC 使用全样本后验和逐点似然，是一种近似的样本外预测准则。更直接的做法是留出测试集，比较模型对未参与拟合数据的预测表现。下面随机留出 25% 的天数作为测试集，并比较三个模型的测试 RMSE、测试对数预测密度和 90% 预测区间覆盖率。
set.seed(1905)
test_id <- sample(seq_len(nrow(delivery_dat)),
size = round(0.25 * nrow(delivery_dat)))
train_dat <- delivery_dat[-test_id, ]
test_dat <- delivery_dat[test_id, ]
holdout_tab <- data.frame()
for (nm in names(model_formulas)) {
fit_train <- fit_bayes_lm(model_formulas[[nm]], train_dat, S = 3500)
mu_test <- posterior_predict(fit_train, test_dat, include_noise = FALSE)
y_test_rep <- posterior_predict(fit_train, test_dat, include_noise = TRUE)
pred_mean <- colMeans(mu_test)
loglik_test <- pointwise_loglik(fit_train, test_dat)
test_lpd <- sum(apply(loglik_test, 2, log_mean_exp))
pred_lower <- apply(y_test_rep, 2, quantile, 0.05)
pred_upper <- apply(y_test_rep, 2, quantile, 0.95)
cover90 <- mean(test_dat$orders >= pred_lower &
test_dat$orders <= pred_upper)
holdout_tab <- rbind(
holdout_tab,
data.frame(
model = nm,
test_rmse = sqrt(mean((test_dat$orders - pred_mean)^2)),
test_lpd = test_lpd,
coverage_90 = cover90
)
)
}
knitr::kable(holdout_tab, row.names = FALSE, digits = 3)| model | test_rmse | test_lpd | coverage_90 |
|---|---|---|---|
| linear | 78.01 | -317.9 | 0.945 |
| quadratic | 61.34 | -304.5 | 0.891 |
| interaction | 60.15 | -303.4 | 0.891 |
RMSE 衡量点预测误差，测试对数预测密度衡量整个预测分布对观测值的支持程度，覆盖率检查预测区间是否大致校准。经济统计应用中，模型选择应与任务目标一致：若目标是调度骑手，预测区间和高峰概率可能比单一 RMSE 更重要；若目标是解释政策效果，模型可解释性也应纳入考虑。
19.8 模型平均与决策
模型比较不一定意味着只能选一个模型。贝叶斯思想中一个自然做法是模型平均：
\[ p(\tilde y\mid y) = \sum_{k=1}^K p(\tilde y\mid y,M_k)P(M_k\mid y). \]
其中 \(M_k\) 是第 \(k\) 个候选模型，\(P(M_k\mid y)\) 是模型后验概率。如果没有计算严格的模型后验概率，也可以用 WAIC 权重构造一个预测导向的近似模型平均。下面用三个模型的 WAIC 权重混合未来情景下的预测分布。
set.seed(1906)
waic_weight <- setNames(waic_tab$waic_weight, waic_tab$model)
future_pred_by_model <- lapply(model_fits, function(fit) {
posterior_predict(fit, future_day, include_noise = TRUE)[, 1]
})
S_avg <- 5000
model_draw <- sample(names(model_fits), size = S_avg, replace = TRUE,
prob = waic_weight[names(model_fits)])
bma_pred <- numeric(S_avg)
for (nm in names(model_fits)) {
id <- which(model_draw == nm)
if (length(id) > 0) {
bma_pred[id] <- sample(future_pred_by_model[[nm]],
size = length(id),
replace = TRUE)
}
}
bma_summary <- data.frame(
quantity = c("model averaged mean",
"model averaged q05",
"model averaged q95",
"Pr(orders > 700)"),
value = c(mean(bma_pred),
quantile(bma_pred, 0.05),
quantile(bma_pred, 0.95),
mean(bma_pred > 700))
)
knitr::kable(bma_summary, row.names = FALSE, digits = 3)| quantity | value |
|---|---|
| model averaged mean | 734.153 |
| model averaged q05 | 634.960 |
| model averaged q95 | 831.230 |
| Pr(orders > 700) | 0.711 |
模型平均的解释要谨慎。它不能修复所有模型共同遗漏的结构，例如所有候选模型都假设同方差正态误差，就仍可能低估雨天需求波动。但当多个合理模型难分高下时，模型平均可以减少“只选一个模型”带来的额外不确定性。