第 18 章 贝叶斯回归模型
回归模型是经济统计中最重要的建模工具之一。本章把贝叶斯思想应用到回归模型中。与经典回归相比，贝叶斯回归的核心变化是：回归系数不再只用点估计和标准误描述，而是用后验分布描述。这样可以更自然地进行区间估计、预测不确定性传播和先验信息整合。
18.1 贝叶斯线性回归
考虑线性模型
\[ \mathbf y\mid\boldsymbol\beta,\sigma^2 \sim N(\mathbf X\boldsymbol\beta,\sigma^2\mathbf I), \]
其中 \(\mathbf y\) 是 \(n\times 1\) 因变量向量，\(\mathbf X\) 是 \(n\times p\) 设计矩阵，\(\boldsymbol\beta\) 是回归系数，\(\sigma^2\) 是误差方差。为了保持本科阶段的计算简洁，本节先把\(\sigma^2\) 视为已知，给 \(\boldsymbol\beta\) 指定正态先验：
\[ \boldsymbol\beta\sim N(\mathbf b_0,\mathbf V_0). \]
根据贝叶斯公式，后验仍为多元正态分布：
\[ \boldsymbol\beta\mid \mathbf y \sim N(\mathbf b_n,\mathbf V_n), \]
其中
\[ \mathbf V_n = \left( \frac{1}{\sigma^2}\mathbf X^\top\mathbf X + \mathbf V_0^{-1} \right)^{-1}, \]
\[ \mathbf b_n = \mathbf V_n \left( \frac{1}{\sigma^2}\mathbf X^\top\mathbf y + \mathbf V_0^{-1}\mathbf b_0 \right). \]
这个公式与 Ridge 回归有相似之处：正态先验会把系数向先验均值收缩。
18.2 案例：消费支出回归
下面使用模拟家庭数据。目标是估计收入和就业状态对消费支出的关系。数据为模拟数据，不对应真实调查。
set.seed(2026)
n <- 180
income <- rlnorm(n, log(8), 0.5)
employed <- rbinom(n, 1, plogis(-0.4 + 0.18 * income))
consumption <- 1.2 + 0.60 * income + 0.7 * employed + rnorm(n, sd = 1.5)
X <- cbind(1, income, employed)
y <- consumption
p <- ncol(X)设先验均值为 0，先验方差较大：
\[ \boldsymbol\beta\sim N(\mathbf 0,10^2\mathbf I). \]
误差标准差暂取 \(\sigma=1.5\)。
sigma <- 1.5
b0 <- rep(0, p)
V0 <- diag(100, p)
Vn <- solve(crossprod(X) / sigma^2 + solve(V0))
bn <- Vn %*% (crossprod(X, y) / sigma^2 + solve(V0, b0))
post_sd <- sqrt(diag(Vn))
out <- data.frame(
term = c("intercept", "income", "employed"),
mean = as.vector(bn),
sd = post_sd,
q025 = as.vector(bn) + qnorm(0.025) * post_sd,
q975 = as.vector(bn) + qnorm(0.975) * post_sd
)
out[, -1] <- round(out[, -1], 3)
out
#> term mean sd q025 q975
#> intercept 1.506 0.271 0.976 2.037
#> income income 0.578 0.025 0.530 0.627
#> employed employed 0.525 0.268 0.000 1.051后验均值可以像普通回归系数一样解释，但区间是后验可信区间。在给定模型和先验下，收入系数的后验分布描述了边际消费倾向的不确定性。
18.3 后验抽样与预测
因为后验是多元正态分布，可以直接从中抽取系数样本。这里用 Cholesky 分解生成多元正态随机数。
set.seed(1)
S <- 5000
L <- chol(Vn)
beta_draw <- matrix(rnorm(S * p), S, p) %*% L + matrix(rep(bn, each = S), S)
new_x <- c(1, income = 10, employed = 1)
mu_pred <- drop(beta_draw %*% new_x)
y_pred <- rnorm(S, mean = mu_pred, sd = sigma)
c(mean_prediction = mean(y_pred),
q025 = quantile(y_pred, 0.025),
q975 = quantile(y_pred, 0.975))
#> mean_prediction q025.2.5% q975.97.5%
#> 7.833 4.968 10.950mu_pred 反映平均消费的不确定性，y_pred 进一步加入个体扰动，因此预测区间通常比均值区间更宽。这是贝叶斯预测的一个优点：参数不确定性和未来观测噪声可以自然合并。
18.4 未知方差与 Gibbs 抽样
如果 \(\sigma^2\) 也未知，可以给它指定先验，例如逆伽马先验。此时后验通常可以通过 Gibbs 抽样：
- 给定 \(\sigma^2\)，从 \(\boldsymbol\beta\mid\sigma^2,\mathbf y\) 中抽样；
- 给定 \(\boldsymbol\beta\)，从 \(\sigma^2\mid\boldsymbol\beta,\mathbf y\) 中抽样；
- 重复以上步骤，得到联合后验样本。
本书不在这里展开完整推导，因为核心思想已经在第 12 章介绍。重要的是认识到：贝叶斯回归的计算通常不是一个公式结束，而是通过解析更新、优化近似或 MCMC 抽样完成。
18.5 贝叶斯逻辑回归
对于二元响应变量，可以使用贝叶斯逻辑回归：
\[ y_i\mid\boldsymbol\beta\sim \operatorname{Bernoulli}(p_i), \]
\[ \log\frac{p_i}{1-p_i} = \mathbf x_i^\top\boldsymbol\beta. \]
若给 \(\boldsymbol\beta\) 正态先验，后验一般没有闭式形式。可以使用第 12 章的随机游走Metropolis，也可以使用 MAP 加 Laplace 近似，或使用第 20 章介绍的 Stan、brms 等工具。
下面只给出 MAP 的计算框架。数据是模拟违约数据。
set.seed(2)
n <- 300
x_income <- as.numeric(scale(rlnorm(n, log(8), 0.55)))
eta <- -1 - 1.0 * x_income
default <- rbinom(n, 1, plogis(eta))
logpost_logit <- function(beta) {
eta <- beta[1] + beta[2] * x_income
loglik <- sum(default * eta - log1p(exp(eta)))
logprior <- sum(dnorm(beta, 0, 5, log = TRUE))
loglik + logprior
}
fit <- optim(c(0, 0), function(b) -logpost_logit(b),
method = "BFGS", hessian = TRUE)
fit$par
#> [1] -0.8746 -0.8634这个例子与第 12 章的 MCMC 例子对应。MAP 给出后验众数，而 MCMC 可以进一步给出完整后验不确定性。