第 19 章 贝叶斯模型比较与后验预测
贝叶斯分析的目标不只是估计参数。很多经济统计问题更关心预测：未来消费是多少？某企业违约概率多大？某项政策下低收入比例会如何变化？本章介绍后验预测、后验预测检验和贝叶斯模型比较的基本思想。
19.1 后验预测分布
设 \(\tilde y\) 表示未来观测值。贝叶斯后验预测分布为
\[ 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)}\) 汇总预测均值和预测区间。
19.2 后验预测示例
继续使用第 18 章的贝叶斯线性回归思想。下面重新生成一组模拟消费数据，并用正态近似得到后验预测。
set.seed(2026)
n <- 160
income <- rlnorm(n, log(8), 0.5)
consumption <- 1.2 + 0.62 * income + rnorm(n, sd = 1.5)
X <- cbind(1, income)
y <- consumption
sigma <- 1.5
b0 <- c(0, 0)
V0 <- diag(100, 2)
Vn <- solve(crossprod(X) / sigma^2 + solve(V0))
bn <- Vn %*% (crossprod(X, y) / sigma^2 + solve(V0, b0))
S <- 4000
L <- chol(Vn)
beta_draw <- matrix(rnorm(S * 2), S, 2) %*% L + matrix(rep(bn, each = S), S)
new_income <- 10
mu_draw <- drop(beta_draw %*% c(1, new_income))
y_tilde <- rnorm(S, mean = mu_draw, sd = sigma)
c(pred_mean = mean(y_tilde),
pred_q025 = quantile(y_tilde, 0.025),
pred_q975 = quantile(y_tilde, 0.975))
#> pred_mean pred_q025.2.5% pred_q975.97.5%
#> 7.492 4.615 10.530后验预测区间可以直接用于表达未来个体消费的不确定性。若研究对象是平均消费，则应使用 mu_draw而不是 y_tilde。
19.3 后验预测检验
后验预测检验的思想是：如果模型合理，那么从模型生成的复制数据应当在关键特征上类似于观测数据。设 \(T(y)\) 是一个检验统计量，例如均值、标准差、最大值或贫困率。对每个后验样本生成复制数据\(y^{\operatorname{rep}}\)，比较
\[ T(y^{\operatorname{rep}}) \]
和
\[ T(y). \]
下面检查模型是否能再现消费数据的样本标准差。
set.seed(1)
T_obs <- sd(y)
T_rep <- numeric(S)
for (s in seq_len(S)) {
y_rep <- rnorm(n, mean = drop(X %*% beta_draw[s, ]), sd = sigma)
T_rep[s] <- sd(y_rep)
}
c(observed_sd = T_obs,
replicated_mean_sd = mean(T_rep),
posterior_predictive_p = mean(T_rep >= T_obs))
#> observed_sd replicated_mean_sd posterior_predictive_p
#> 3.3565 3.3998 0.6165后验预测 p 值不是频率学派检验的 p 值。它更像一个模型检查指标：若观测统计量位于复制分布的极端尾部，说明模型可能遗漏了某些结构，例如异方差、非线性或厚尾误差。
19.4 贝叶斯模型比较
贝叶斯模型比较有多种方式。最经典的是边际似然
\[ 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)}. \]
Bayes 因子衡量数据相对支持哪个模型。但它对先验可能敏感，且边际似然计算并不总是容易。
预测导向的贝叶斯模型比较常使用 WAIC 或留一交叉验证。WAIC 的思想是用逐点对数预测密度衡量预测能力，并对有效复杂度进行惩罚。给定后验样本 \(\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_{\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 越小，通常表示预测表现越好。
19.5 案例：比较两个消费模型
下面比较两个模型：
\[ M_1: C_i=\beta_0+\beta_1Y_i+\varepsilon_i, \]
\[ M_2: C_i=\beta_0+\beta_1\log Y_i+\varepsilon_i. \]
这里用简单训练测试划分展示预测比较。完整贝叶斯项目中可以使用后验预测、WAIC 或 LOO。
set.seed(2)
id_train <- sample(seq_len(n), size = round(0.75 * n))
train <- data.frame(y = y[id_train], income = income[id_train])
test <- data.frame(y = y[-id_train], income = income[-id_train])
fit_m1 <- lm(y ~ income, data = train)
fit_m2 <- lm(y ~ log(income), data = train)
pred_m1 <- predict(fit_m1, newdata = test)
pred_m2 <- predict(fit_m2, newdata = test)
c(RMSE_m1 = sqrt(mean((test$y - pred_m1)^2)),
RMSE_m2 = sqrt(mean((test$y - pred_m2)^2)),
AIC_m1 = AIC(fit_m1),
AIC_m2 = AIC(fit_m2))
#> RMSE_m1 RMSE_m2 AIC_m1 AIC_m2
#> 1.349 1.685 432.671 476.757模型比较不应只看数字。若两个模型预测差异很小，通常应优先选择解释更清楚、与经济机制更一致的模型。若预测是核心目标，则应更重视样本外表现。