第 18 章 贝叶斯回归模型
回归模型是经济统计中最重要的建模工具之一。前面几章讨论了贝叶斯更新、后验计算和 MCMC，本章把这些思想放入回归模型。与经典回归相比，贝叶斯回归的核心变化是：回归系数不再只用点估计和标准误描述，而是用后验分布描述。这样可以更自然地回答下面这类问题：
- 某项消费券政策使人均消费增加的概率有多大？
- 政策带动的消费额超过财政补贴面额的概率有多大？
- 对一个新家庭或一家新企业，预测值的不确定性有多大？
- 当样本量不大或变量较多时，怎样把合理的先验信息转化为稳定估计？
本章的学习重点有三个。第一，理解贝叶斯线性回归后验公式的含义；第二，能够用后验样本解释回归系数、预测均值和未来观测；第三，了解逻辑回归等非共轭模型为什么需要优化、Laplace 近似或 MCMC。
18.1 贝叶斯线性回归
考虑线性模型
\[ \mathbf y\mid\boldsymbol\beta,\sigma^2 \sim N(\mathbf X\boldsymbol\beta,\sigma^2\mathbf I_n). \]
其中 \(\mathbf y=(y_1,\ldots,y_n)^\top\) 是 \(n\times 1\) 因变量向量，\(\mathbf X\) 是 \(n\times p\)设计矩阵，\(\boldsymbol\beta=(\beta_1,\ldots,\beta_p)^\top\) 是回归系数，\(\sigma^2\) 是误差方差，\(\mathbf I_n\) 是 \(n\) 阶单位矩阵。若 \(\mathbf X\) 第一列为 1，则对应系数是截距项。
为了先看清主线，本节暂把 \(\sigma^2\) 视为已知，并给 \(\boldsymbol\beta\) 指定多元正态先验：
\[ \boldsymbol\beta\sim N(\mathbf b_0,\mathbf V_0). \]
其中 \(\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). \]
这里 \(\mathbf V_n\) 是后验协方差矩阵，\(\mathbf b_n\) 是后验均值向量。注意后验精度矩阵\(\mathbf V_n^{-1}\) 可以写成
\[ \mathbf V_n^{-1} = \frac{1}{\sigma^2}\mathbf X^\top\mathbf X+\mathbf V_0^{-1}. \]
这说明后验信息来自两部分：数据提供的精度 \(\mathbf X^\top\mathbf X/\sigma^2\)，以及先验提供的精度\(\mathbf V_0^{-1}\)。如果先验方差很大，\(\mathbf V_0^{-1}\) 很小，后验会更接近普通最小二乘；如果先验方差较小，系数会更明显地向先验均值收缩。
这种收缩与 Ridge 回归有直接联系。若 \(\mathbf b_0=\mathbf 0\) 且 \(\mathbf V_0=\tau^2\mathbf I\)，后验众数等价于最小化
\[ \frac{1}{2\sigma^2} (\mathbf y-\mathbf X\boldsymbol\beta)^\top (\mathbf y-\mathbf X\boldsymbol\beta) + \frac{1}{2\tau^2}\boldsymbol\beta^\top\boldsymbol\beta. \]
第一项衡量拟合误差，第二项惩罚系数大小。\(\tau^2\) 越小，惩罚越强。
18.2 案例：夜间商圈消费券的带动效应
下面使用一个模拟案例。某城市想评估“夜间商圈消费券”是否真的带来增量消费。消费券面额为 15 元，部分居民随机收到消费券。研究者记录一周内夜间商圈消费额，并同时记录收入指数、到商圈距离和当周雨天数。数据是模拟数据，不对应真实调查。
这个案例有一点有意思的地方：政策部门不仅关心“消费券有没有正效应”，还关心“平均带动消费额是否超过 15 元面额”。第二个问题不能只靠一个回归系数点估计回答，更适合用后验概率表达。
set.seed(1801)
n <- 240
coupon <- rbinom(n, size = 1, prob = 0.5)
income_index <- rnorm(n)
distance <- runif(n, min = 0.3, max = 8)
rainy_days <- rbinom(n, size = 4, prob = 0.35)
baseline <- 62 + 9 * income_index - 1.8 * distance - 3.2 * rainy_days
extra_effect <- 17 + 3 * (income_index < -0.5)
night_spend <- baseline + coupon * extra_effect + rnorm(n, sd = 18)
coupon_dat <- data.frame(
night_spend = night_spend,
coupon = coupon,
income_index = income_index,
distance = distance,
rainy_days = rainy_days
)
head(coupon_dat)
#> night_spend coupon income_index distance rainy_days
#> 1 48.62 0 -0.26605 1.961 0
#> 2 30.01 0 0.03238 1.559 0
#> 3 61.52 0 -0.98459 1.801 1
#> 4 70.67 1 -2.25775 5.541 1
#> 5 51.91 0 -1.13338 4.265 1
#> 6 50.78 1 1.37992 6.570 2设回归模型为
\[ y_i = \beta_0 +\beta_1\operatorname{coupon}_i +\beta_2\operatorname{income}_i +\beta_3\operatorname{distance}_i +\beta_4\operatorname{rain}_i +\varepsilon_i, \]
\[ \varepsilon_i\sim N(0,\sigma^2). \]
其中 \(y_i\) 是第 \(i\) 个居民一周夜间商圈消费额，\(\operatorname{coupon}_i\) 表示是否收到消费券，\(\operatorname{income}_i\) 是标准化收入指数，\(\operatorname{distance}_i\) 是到商圈距离，\(\operatorname{rain}_i\) 是当周雨天数。系数 \(\beta_1\) 表示在控制这些变量后，收到消费券带来的平均消费差异。
为了使计算简洁，先用普通最小二乘残差标准差作为 \(\sigma\) 的插入估计。这里的重点不是把\(\sigma^2\) 当作完全已知，而是展示在线性回归中如何从正态后验得到可解释的经济结论。后面会进一步讨论 \(\sigma^2\) 未知时的 Gibbs 抽样。
X <- model.matrix(night_spend ~ coupon + income_index + distance + rainy_days,
data = coupon_dat)
y <- coupon_dat$night_spend
p <- ncol(X)
ols_fit <- lm(night_spend ~ coupon + income_index + distance + rainy_days,
data = coupon_dat)
sigma <- summary(ols_fit)$sigma
b0 <- c(60, 0, 0, 0, 0)
prior_sd <- c(30, 20, 20, 5, 8)
V0 <- diag(prior_sd^2)
Vn <- solve(crossprod(X) / sigma^2 + solve(V0))
bn <- Vn %*% (crossprod(X, y) / sigma^2 + solve(V0, b0))
post_sd <- sqrt(diag(Vn))
coupon_summary <- data.frame(
term = colnames(X),
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
)
knitr::kable(coupon_summary, row.names = FALSE, digits = 3)| term | mean | sd | q025 | q975 |
|---|---|---|---|---|
| (Intercept) | 62.439 | 3.485 | 55.608 | 69.270 |
| coupon | 19.094 | 2.574 | 14.050 | 24.138 |
| income_index | 8.248 | 1.229 | 5.839 | 10.656 |
| distance | -1.306 | 0.560 | -2.404 | -0.208 |
| rainy_days | -4.326 | 1.290 | -6.855 | -1.797 |
表中的 coupon 行是本案例最关心的政策系数。后验均值可以像普通回归系数一样解释，但可信区间的含义是：在给定模型、数据和先验后，系数位于该区间内的后验概率约为 95%。这与经典置信区间的解释不同，后者是关于重复抽样程序长期覆盖率的陈述。
18.3 后验抽样与政策概率
贝叶斯回归的一个优点是：许多政策问题可以直接转化为后验概率。若消费券面额为 15 元，则政策部门可以关注
\[ P(\beta_1>0\mid \mathbf y) \]
和
\[ P(\beta_1>15\mid \mathbf y). \]
第一个概率衡量消费券带来正向消费增量的证据，第二个概率衡量平均消费带动额超过券面额的证据。因为\(\boldsymbol\beta\mid\mathbf y\) 是多元正态分布，可以直接模拟后验样本。
rmvnorm_chol <- function(n, mean, Sigma) {
z <- matrix(rnorm(n * length(mean)), nrow = n)
sweep(z %*% chol(Sigma), 2, mean, FUN = "+")
}
set.seed(1802)
S <- 6000
beta_draw <- rmvnorm_chol(S, as.vector(bn), Vn)
colnames(beta_draw) <- colnames(X)
coupon_effect <- beta_draw[, "coupon"]
face_value <- 15
policy_prob <- data.frame(
quantity = c("Pr(effect > 0)",
"Pr(effect > face value)",
"posterior mean of effect",
"posterior sd of effect"),
value = c(mean(coupon_effect > 0),
mean(coupon_effect > face_value),
mean(coupon_effect),
sd(coupon_effect))
)
knitr::kable(policy_prob, row.names = FALSE, digits = 3)| quantity | value |
|---|---|
| Pr(effect > 0) | 1.000 |
| Pr(effect > face value) | 0.947 |
| posterior mean of effect | 19.097 |
| posterior sd of effect | 2.536 |
这些概率比“显著或不显著”更贴近政策讨论。如果 \(P(\beta_1>0\mid\mathbf y)\) 很高，但\(P(\beta_1>15\mid\mathbf y)\) 不高，说明消费券可能确实刺激消费，但平均带动额未必足以覆盖券面额。
后验样本也可以用来比较不同居民画像下的平均消费。下面比较两个居民：一个住得较近、收入指数较高；另一个住得较远、收入指数较低且雨天较多。我们分别计算收到和未收到消费券时的预测均值差异。
profile_dat <- data.frame(
coupon = c(0, 1, 0, 1),
income_index = c(0.8, 0.8, -0.8, -0.8),
distance = c(1.2, 1.2, 6.5, 6.5),
rainy_days = c(1, 1, 3, 3)
)
X_profile <- model.matrix(~ coupon + income_index + distance + rainy_days,
data = profile_dat)
mu_profile <- beta_draw %*% t(X_profile)
lift_near <- mu_profile[, 2] - mu_profile[, 1]
lift_far <- mu_profile[, 4] - mu_profile[, 3]
profile_summary <- data.frame(
profile = c("near_high_income", "far_low_income"),
mean_lift = c(mean(lift_near), mean(lift_far)),
q025 = c(quantile(lift_near, 0.025), quantile(lift_far, 0.025)),
q975 = c(quantile(lift_near, 0.975), quantile(lift_far, 0.975))
)
knitr::kable(profile_summary, row.names = FALSE, digits = 3)| profile | mean_lift | q025 | q975 |
|---|---|---|---|
| near_high_income | 19.1 | 14.07 | 24.14 |
| far_low_income | 19.1 | 14.07 | 24.14 |
在当前线性模型中，消费券效应被设为相同的 \(\beta_1\)，所以两个画像的平均增量相同。如果研究者怀疑低收入或远距离居民的消费券反应不同，应加入交互项，例如 coupon:income_index 或 coupon:distance。这也提醒我们：贝叶斯方法不会自动替代模型设定，后验解释仍然依赖回归方程本身。
18.4 先验敏感性分析
贝叶斯回归需要报告先验。一个实用做法是比较不同先验下关键结论是否稳定。下面比较三种先验：
neutral：消费券效应先验标准差较大，表示较弱先验；skeptical：消费券效应先验均值为 0，标准差较小，表示对政策效果较保守；optimistic：消费券效应先验均值为 10，表示先前试点或专家判断认为可能存在正效应。
fit_bayes_lm <- function(b0, prior_sd, X, y, sigma) {
V0 <- diag(prior_sd^2)
Vn <- solve(crossprod(X) / sigma^2 + solve(V0))
bn <- Vn %*% (crossprod(X, y) / sigma^2 + solve(V0, b0))
list(mean = as.vector(bn), cov = Vn)
}
prior_list <- list(
neutral = list(b0 = c(60, 0, 0, 0, 0),
sd = c(30, 20, 20, 5, 8)),
skeptical = list(b0 = c(60, 0, 0, 0, 0),
sd = c(30, 6, 20, 5, 8)),
optimistic = list(b0 = c(60, 10, 0, 0, 0),
sd = c(30, 8, 20, 5, 8))
)
sens <- data.frame()
for (nm in names(prior_list)) {
fit_s <- fit_bayes_lm(prior_list[[nm]]$b0,
prior_list[[nm]]$sd,
X, y, sigma)
idx <- which(colnames(X) == "coupon")
m <- fit_s$mean[idx]
se <- sqrt(fit_s$cov[idx, idx])
sens <- rbind(
sens,
data.frame(
prior = nm,
coupon_mean = m,
coupon_sd = se,
prob_gt_0 = 1 - pnorm(0, mean = m, sd = se),
prob_gt_face_value = 1 - pnorm(face_value, mean = m, sd = se)
)
)
}
knitr::kable(sens, row.names = FALSE, digits = 3)| prior | coupon_mean | coupon_sd | prob_gt_0 | prob_gt_face_value |
|---|---|---|---|---|
| neutral | 19.09 | 2.574 | 1 | 0.944 |
| skeptical | 16.36 | 2.382 | 1 | 0.715 |
| optimistic | 18.52 | 2.469 | 1 | 0.923 |
如果不同先验下政策概率差异很小，说明数据对该结论提供了较强信息。如果差异很大，报告时就应说明结论依赖哪些先验判断。对本科阶段而言，先验敏感性分析比追求某个“唯一正确”的先验更重要。
18.5 未知方差与 Gibbs 抽样
前面为了得到简洁公式，把 \(\sigma^2\) 当作已知。在更完整的贝叶斯线性回归中，误差方差也可以作为未知参数。设
\[ \boldsymbol\beta\sim N(\mathbf b_0,\mathbf V_0), \]
\[ \sigma^2\sim \operatorname{IG}(a_0,d_0), \]
其中 \(\operatorname{IG}(a_0,d_0)\) 表示逆伽马分布，\(a_0\) 是形状参数，\(d_0\) 是尺度参数。这里采用约定：若 \(\sigma^2\sim\operatorname{IG}(a,d)\)，则 \(1/\sigma^2\sim\operatorname{Gamma}(a,d)\)，其中 Gamma 的第二个参数为 rate。
在独立正态先验下，可以用 Gibbs 抽样交替更新：
- 给定当前 \(\sigma^2\)，从\(\boldsymbol\beta\mid\sigma^2,\mathbf y\) 中抽样；
- 给定当前 \(\boldsymbol\beta\)，从\(\sigma^2\mid\boldsymbol\beta,\mathbf y\) 中抽样；
- 重复以上步骤，丢弃前若干次 burn-in 后汇总样本。
给定 \(\sigma^2\) 时，\(\boldsymbol\beta\) 的条件后验仍是正态分布：
\[ \boldsymbol\beta\mid\sigma^2,\mathbf y \sim N(\mathbf b_\sigma,\mathbf V_\sigma), \]
其中
\[ \mathbf V_\sigma = \left( \frac{1}{\sigma^2}\mathbf X^\top\mathbf X+\mathbf V_0^{-1} \right)^{-1}, \]
\[ \mathbf b_\sigma = \mathbf V_\sigma \left( \frac{1}{\sigma^2}\mathbf X^\top\mathbf y+\mathbf V_0^{-1}\mathbf b_0 \right). \]
给定 \(\boldsymbol\beta\) 时，令残差平方和
\[ RSS(\boldsymbol\beta) = (\mathbf y-\mathbf X\boldsymbol\beta)^\top (\mathbf y-\mathbf X\boldsymbol\beta), \]
则
\[ \sigma^2\mid\boldsymbol\beta,\mathbf y \sim \operatorname{IG} \left( a_0+\frac{n}{2}, d_0+\frac{RSS(\boldsymbol\beta)}{2} \right). \]
下面对消费券案例运行一个简短的 Gibbs 抽样。
set.seed(1803)
S_gibbs <- 7000
burn <- 1000
a0 <- 2
d0 <- sigma^2
V0_inv <- solve(V0)
XtX <- crossprod(X)
Xty <- crossprod(X, y)
beta_cur <- coef(ols_fit)
sigma2_cur <- sigma^2
beta_gibbs <- matrix(NA_real_, nrow = S_gibbs, ncol = p)
sigma2_gibbs <- numeric(S_gibbs)
colnames(beta_gibbs) <- colnames(X)
for (s in seq_len(S_gibbs)) {
V_beta <- solve(XtX / sigma2_cur + V0_inv)
m_beta <- V_beta %*% (Xty / sigma2_cur + V0_inv %*% b0)
beta_cur <- as.numeric(rmvnorm_chol(1, as.vector(m_beta), V_beta))
resid <- y - drop(X %*% beta_cur)
shape <- a0 + length(y) / 2
rate <- d0 + sum(resid^2) / 2
sigma2_cur <- 1 / rgamma(1, shape = shape, rate = rate)
beta_gibbs[s, ] <- beta_cur
sigma2_gibbs[s] <- sigma2_cur
}
keep <- (burn + 1):S_gibbs
gibbs_summary <- data.frame(
quantity = c("coupon_mean", "coupon_q025", "coupon_q975",
"Pr(coupon > 15)", "sigma_mean"),
value = c(mean(beta_gibbs[keep, "coupon"]),
quantile(beta_gibbs[keep, "coupon"], 0.025),
quantile(beta_gibbs[keep, "coupon"], 0.975),
mean(beta_gibbs[keep, "coupon"] > face_value),
mean(sqrt(sigma2_gibbs[keep])))
)
knitr::kable(gibbs_summary, row.names = FALSE, digits = 3)| quantity | value |
|---|---|
| coupon_mean | 19.069 |
| coupon_q025 | 14.120 |
| coupon_q975 | 23.976 |
| Pr(coupon > 15) | 0.946 |
| sigma_mean | 19.991 |
Gibbs 抽样的结果通常会比插入固定 \(\sigma\) 的后验略宽，因为它额外考虑了误差方差的不确定性。实际报告时，还应检查链的轨迹、有效样本量和自相关；这些诊断思想已经在第 12 章介绍。
18.6 贝叶斯逻辑回归
许多经济统计问题的因变量是二元变量，例如是否违约、是否就业、是否参加某项培训、是否接受电子支付。对于二元响应变量，可以使用逻辑回归：
\[ y_i\mid\boldsymbol\beta\sim \operatorname{Bernoulli}(p_i), \]
\[ \log\frac{p_i}{1-p_i} = \mathbf x_i^\top\boldsymbol\beta. \]
其中 \(y_i\) 取 0 或 1，\(p_i=P(y_i=1\mid\mathbf x_i,\boldsymbol\beta)\)，\(\mathbf x_i\) 是解释变量向量。若给 \(\boldsymbol\beta\) 正态先验
\[ \boldsymbol\beta\sim N(\mathbf 0,\mathbf V_0), \]
后验密度正比于
\[ p(\boldsymbol\beta\mid \mathbf y) \propto \prod_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i} \times \exp\left( -\frac{1}{2}\boldsymbol\beta^\top\mathbf V_0^{-1}\boldsymbol\beta \right). \]
这个后验一般没有像线性回归那样的闭式正态形式。常见做法包括随机游走 Metropolis、HamiltonianMonte Carlo、变分近似，以及本章使用的 MAP 加 Laplace 近似。
18.7 案例：小微企业违约预警
下面构造一个模拟的小微企业贷款数据。变量包括资产负债率、现金流指数、是否年轻企业和近三个月逾期次数。目标是估计企业违约概率。数据仍然是模拟数据，不对应任何真实银行数据。
set.seed(1804)
n_firm <- 450
debt_ratio <- runif(n_firm, min = 0.1, max = 0.9)
cashflow <- rnorm(n_firm)
young <- rbinom(n_firm, size = 1, prob = 0.35)
late_pay <- rpois(n_firm,
lambda = exp(-1.1 + 1.5 * debt_ratio - 0.5 * cashflow))
eta_true <- -3.0 + 2.6 * debt_ratio - 0.9 * cashflow +
0.7 * young + 0.55 * late_pay
default <- rbinom(n_firm, size = 1, prob = plogis(eta_true))
firm_dat <- data.frame(
default = default,
debt_ratio = debt_ratio,
cashflow = cashflow,
young = young,
late_pay = late_pay
)
mean(firm_dat$default)
#> [1] 0.3为了数值稳定，先定义
\[ \log(1+\exp x) \]
的稳定计算函数。逻辑回归的对数似然可写为
\[ \ell(\boldsymbol\beta) = \sum_{i=1}^n \left[ y_i\eta_i-\log(1+\exp(\eta_i)) \right], \]
其中 \(\eta_i=\mathbf x_i^\top\boldsymbol\beta\)。再加上正态先验的对数密度，就得到对数后验。
X_logit <- model.matrix(default ~ debt_ratio + cashflow + young + late_pay,
data = firm_dat)
y_logit <- firm_dat$default
prior_sd_logit <- c(5, 2.5, 2.5, 2.5, 2.5)
log1pexp <- function(x) {
ifelse(x > 0, x + log1p(exp(-x)), log1p(exp(x)))
}
logpost_logit <- function(beta) {
eta <- drop(X_logit %*% beta)
loglik <- sum(y_logit * eta - log1pexp(eta))
logprior <- sum(dnorm(beta, mean = 0, sd = prior_sd_logit, log = TRUE))
loglik + logprior
}
fit_logit <- optim(rep(0, ncol(X_logit)),
function(b) -logpost_logit(b),
method = "BFGS",
hessian = TRUE)
beta_map <- fit_logit$par
cov_laplace <- solve(fit_logit$hessian)
set.seed(1805)
beta_logit_draw <- rmvnorm_chol(6000, beta_map, cov_laplace)
colnames(beta_logit_draw) <- colnames(X_logit)
logit_summary <- data.frame(
term = colnames(X_logit),
mean = colMeans(beta_logit_draw),
q025 = apply(beta_logit_draw, 2, quantile, 0.025),
q975 = apply(beta_logit_draw, 2, quantile, 0.975)
)
knitr::kable(logit_summary, row.names = FALSE, digits = 3)| term | mean | q025 | q975 |
|---|---|---|---|
| (Intercept) | -2.639 | -3.283 | -1.969 |
| debt_ratio | 2.388 | 1.276 | 3.465 |
| cashflow | -0.917 | -1.200 | -0.635 |
| young | 0.293 | -0.207 | 0.783 |
| late_pay | 0.355 | 0.115 | 0.601 |
逻辑回归系数表示对数赔率的变化。为了便于解释，常把系数指数化为赔率比。若某变量的系数为\(\beta_j\)，则该变量增加 1 个单位时，违约赔率乘以 \(\exp(\beta_j)\)，其他变量保持不变。
odds_summary <- data.frame(
term = colnames(X_logit)[-1],
odds_ratio_mean = colMeans(exp(beta_logit_draw[, -1])),
odds_ratio_q025 = apply(exp(beta_logit_draw[, -1]), 2, quantile, 0.025),
odds_ratio_q975 = apply(exp(beta_logit_draw[, -1]), 2, quantile, 0.975)
)
knitr::kable(odds_summary, row.names = FALSE, digits = 3)| term | odds_ratio_mean | odds_ratio_q025 | odds_ratio_q975 |
|---|---|---|---|
| debt_ratio | 12.639 | 3.581 | 31.988 |
| cashflow | 0.404 | 0.301 | 0.530 |
| young | 1.384 | 0.813 | 2.187 |
| late_pay | 1.438 | 1.121 | 1.824 |
最后计算两个企业画像的违约概率后验分布：一家财务状况较稳健的企业，和一家负债率高、现金流偏弱且近期多次逾期的企业。
risk_dat <- data.frame(
debt_ratio = c(0.35, 0.78),
cashflow = c(0.6, -0.8),
young = c(0, 1),
late_pay = c(0, 3)
)
X_risk <- model.matrix(~ debt_ratio + cashflow + young + late_pay,
data = risk_dat)
risk_draw <- plogis(beta_logit_draw %*% t(X_risk))
risk_summary <- data.frame(
profile = c("stable_firm", "strained_firm"),
mean_default_prob = colMeans(risk_draw),
q025 = apply(risk_draw, 2, quantile, 0.025),
q975 = apply(risk_draw, 2, quantile, 0.975),
prob_default_gt_20pct = colMeans(risk_draw > 0.20)
)
knitr::kable(risk_summary, row.names = FALSE, digits = 3)| profile | mean_default_prob | q025 | q975 | prob_default_gt_20pct |
|---|---|---|---|---|
| stable_firm | 0.089 | 0.057 | 0.132 | 0 |
| strained_firm | 0.783 | 0.656 | 0.880 | 1 |
这个结果比单纯给出“违约/不违约”分类更有信息。信贷风控中，概率本身往往比硬分类更有用，因为银行可以结合风险偏好、贷款利率、抵押品和监管要求设定不同阈值。
18.8 本章小结
贝叶斯回归用后验分布描述回归系数不确定性。在线性回归中，正态先验与正态似然结合可以得到解析正态后验；后验抽样可以直接计算政策概率、可信区间和预测分布。当误差方差未知时，可以用 Gibbs抽样同时处理系数和方差不确定性。对于逻辑回归等非共轭模型，后验通常没有闭式形式，需要优化、Laplace 近似或 MCMC。贝叶斯回归与正则化之间也有密切联系：先验可以看作对系数大小和结构的概率约束。
18.9 练习
- 在消费券案例中，把消费券面额从 15 改为 10 和 20，分别计算 \(P(\beta_1>\text{face value}\mid y)\)。
- 在消费券模型中加入交互项
coupon:income_index，解释低收入和高收入居民的消费券效应是否不同。 - 比较
neutral与skeptical先验下的消费券后验均值和政策概率，并解释为什么会有差异。 - 在 Gibbs 抽样中把 burn-in 改为 500 和 2000，比较消费券系数后验摘要是否稳定。
- 在小微企业违约案例中，计算
late_pay的赔率比后验区间，并解释它的经济含义。 - 说明为什么贝叶斯逻辑回归不能像正态线性回归那样直接得到解析正态后验。