第 16 章 贝叶斯思维与贝叶斯推断
前面几部分主要从频率学派和计算算法角度讨论统计问题：参数是固定但未知的，样本是随机的，估计量的性质通过重复抽样理解。贝叶斯方法提供了另一种组织不确定性的语言：在看到数据之前，我们用先验分布描述对未知参数的认识；看到数据之后，用后验分布更新这种认识。
对经济统计而言，贝叶斯思想并不神秘。宏观预测会使用历史经验和专家判断；小地区统计会把相邻地区的信息借过来；信用风险模型会把历史违约率与新样本结合；抽样调查中也常常需要在样本信息不足时引入外部信息。贝叶斯方法把这些做法放进统一的概率框架。
本章先从普通条件概率出发，说明贝叶斯更新为什么经常和直觉不同；再进入参数推断，讨论先验、似然、后验、可信区间和几个共轭模型。后续章节会把这些思想推进到更一般的贝叶斯计算问题。
16.1 从条件概率到贝叶斯更新
对两个事件 \(A\) 和 \(B\)，条件概率定义为
\[ P(A\mid B)=\frac{P(A\cap B)}{P(B)}. \]
把 \(A\) 理解为“我们关心的状态”，把 \(B\) 理解为“新观察到的信息”，这个公式就给出了看到信息以后如何修正对状态的判断。贝叶斯公式的事件形式为
\[ P(A\mid B)= \frac{P(B\mid A)P(A)}{P(B\mid A)P(A)+P(B\mid A^c)P(A^c)}. \]
这里 \(P(A)\) 是更新前的概率，\(P(A\mid B)\) 是更新后的概率，\(P(B\mid A)\) 描述新信息在不同状态下出现的可能性。
一个经典例子是诊断检测。设某种检测的灵敏度为 \(P(+\mid H)=0.93\)，特异度为\(P(-\mid H^c)=0.99\)，而人群患病率为 \(P(H)=0.00148\)。如果某人检测为阳性，患病概率是多少？
sensitivity <- 0.93
specificity <- 0.99
prevalence <- 1.48 / 1000
p_positive <- sensitivity * prevalence +
(1 - specificity) * (1 - prevalence)
posterior_after_one_positive <- sensitivity * prevalence / p_positive
c(p_positive = p_positive,
posterior_after_one_positive = posterior_after_one_positive)
#> p_positive posterior_after_one_positive
#> 0.01136 0.12114检测准确率很高，并不意味着阳性以后患病概率也接近 1。原因是患病率很低，假阳性在阳性样本中仍会占据相当比例。贝叶斯更新要求我们同时看“检测有多准”和“事件原本有多稀有”。
如果第一次阳性以后再做一次独立检测，第一次的后验概率就可以作为第二次更新的先验概率：
prior_after_first <- posterior_after_one_positive
posterior_after_two_positives <- sensitivity * prior_after_first /
(sensitivity * prior_after_first +
(1 - specificity) * (1 - prior_after_first))
c(after_first_positive = posterior_after_one_positive,
after_second_positive = posterior_after_two_positives)
#> after_first_positive after_second_positive
#> 0.1211 0.9276这个例子展示了贝叶斯学习的基本结构：先验信息经过数据更新为后验；如果继续收集新数据，当前后验又可以成为下一轮更新的先验。
16.2 概率解释：频率与信念
频率学派通常把概率理解为大量重复试验中的长期频率。若事件 \(E\) 在 \(n\) 次试验中出现 \(n_E\) 次，则
\[ P(E)=\lim_{n\rightarrow\infty}\frac{n_E}{n}. \]
贝叶斯方法允许用概率表达对未知状态或未知参数的不确定性。例如，天气预报中的“明天下雨概率为70%”很难理解为同一个明天重复无穷多次；更自然的解释是，在当前信息和模型下，我们对明天下雨这件事赋予 70% 的概率。类似地，贝叶斯参数推断会把未知参数也放进概率分布中描述。
16.3 贝叶斯公式
设 \(\theta\) 是未知参数，\(y\) 是观测数据。贝叶斯公式写作
\[ p(\theta\mid y) = \frac{p(y\mid\theta)p(\theta)}{p(y)}. \]
其中：
- \(p(\theta)\) 是先验分布，表示看到数据之前对参数的认识；
- \(p(y\mid\theta)\) 是似然函数，表示给定参数时数据出现的可能性；
- \(p(\theta\mid y)\) 是后验分布，表示看到数据之后对参数的更新认识；
- \(p(y)\) 是边际似然，也称证据项，用来保证后验分布积分为 1。
边际似然为
\[ p(y)=\int p(y\mid\theta)p(\theta)\,d\theta. \]
在很多计算中，\(p(y)\) 很难求，但如果只关心后验密度的相对大小，可以写成
\[ p(\theta\mid y)\propto p(y\mid\theta)p(\theta). \]
这就是贝叶斯计算的核心形式：后验正比于似然乘以先验。
16.4 例：睡眠比例的贝叶斯更新
假设我们关心大学生中每天睡眠不少于 8 小时的比例。令 \(\theta\) 表示总体比例，在 \(n\) 名学生中有\(x\) 人睡眠不少于 8 小时。模型为
\[ X\mid\theta\sim \operatorname{Binomial}(n,\theta). \]
16.4.1 离散先验
一种直观做法是先列出 \(\theta\) 的若干可能取值，并给它们分配先验权重。例如，研究者认为\(\theta\) 很可能小于 0.5，于是给较小的比例更高权重。
theta_grid <- seq(0.05, 0.95, by = 0.1)
prior_weight <- c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0)
prior_prob <- prior_weight / sum(prior_weight)
data.frame(theta = theta_grid,
prior = round(prior_prob, 3))
#> theta prior
#> 1 0.05 0.035
#> 2 0.15 0.180
#> 3 0.25 0.277
#> 4 0.35 0.249
#> 5 0.45 0.159
#> 6 0.55 0.073
#> 7 0.65 0.024
#> 8 0.75 0.003
#> 9 0.85 0.000
#> 10 0.95 0.000若样本中有 27 名学生，其中 11 人睡眠不少于 8 小时，则似然函数正比于
\[ L(\theta)\propto \theta^{11}(1-\theta)^{16}. \]
在离散网格上，后验概率只需要把先验概率乘以似然，再归一化：
n_sleep <- 27
x_sleep <- 11
likelihood <- theta_grid^x_sleep * (1 - theta_grid)^(n_sleep - x_sleep)
posterior_prob <- prior_prob * likelihood
posterior_prob <- posterior_prob / sum(posterior_prob)
data.frame(theta = theta_grid,
prior = round(prior_prob, 3),
posterior = round(posterior_prob, 3))
#> theta prior posterior
#> 1 0.05 0.035 0.000
#> 2 0.15 0.180 0.002
#> 3 0.25 0.277 0.129
#> 4 0.35 0.249 0.477
#> 5 0.45 0.159 0.334
#> 6 0.55 0.073 0.056
#> 7 0.65 0.024 0.002
#> 8 0.75 0.003 0.000
#> 9 0.85 0.000 0.000
#> 10 0.95 0.000 0.000这个离散例子很好地展示了贝叶斯更新的机制：先验给出更新前各个取值的权重，似然根据数据重新调整这些权重，后验则是二者结合后的归一化结果。
16.4.2 Beta 先验
离散网格便于理解，但连续比例参数通常使用 Beta 分布建模。选择 Beta 先验
\[ \theta\sim \operatorname{Beta}(a,b). \]
Beta 分布的密度为
\[ p(\theta) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \theta^{a-1}(1-\theta)^{b-1}, \quad 0<\theta<1. \]
其中 \(a\) 和 \(b\) 控制先验形状。结合二项似然后，后验仍为 Beta 分布：
\[ \theta\mid x \sim \operatorname{Beta}(a+x,b+n-x). \]
这称为共轭分析：先验和后验属于同一个分布族。先验中的 \(a\) 和 \(b\) 可以粗略理解为先验成功次数和失败次数的信息量。若研究者认为 \(\theta\) 的中位数大约在 0.3 附近，且 90% 概率小于 0.5，可以用形状参数 \(a=3.26,b=7.19\) 近似表达这一先验判断。
a <- 3.26
b <- 7.19
a_post <- a + x_sleep
b_post <- b + n_sleep - x_sleep
c(prior_mean = a / (a + b),
sample_prop = x_sleep / n_sleep,
posterior_mean = a_post / (a_post + b_post),
posterior_prob_gt_0.5 = 1 - pbeta(0.5, a_post, b_post),
posterior_q05 = qbeta(0.05, a_post, b_post),
posterior_q95 = qbeta(0.95, a_post, b_post))
#> prior_mean sample_prop posterior_mean
#> 0.31196 0.40741 0.38077
#> posterior_prob_gt_0.5 posterior_q05 posterior_q95
#> 0.06902 0.25553 0.51336后验均值位于先验均值和样本比例之间。这里样本量不大，因此先验仍能看出影响；随着样本量增加，似然会逐渐主导后验。
也可以直接从后验分布中抽样，用模拟方式回答后验问题：
16.5 后验区间与解释
贝叶斯推断常用可信区间（credible interval）描述参数不确定性。例如，95% 后验可信区间\([L,U]\) 满足
\[ P(L\leq \theta\leq U\mid y)=0.95. \]
它的解释是：在给定模型、先验和数据后，参数落在该区间内的后验概率为 95%。这与频率学派置信区间的重复抽样解释不同。二者都可以用于不确定性表达，但背后的概率含义不同。
例如，一项民意调查可能报告“95% 置信区间为 60% 到 64%”。在频率学派解释中，未知总体比例是固定的，随机的是抽样过程；如果无限次重复同样的抽样和构造区间程序，约 95% 的区间会覆盖真实比例。对于已经得到的这一个区间，不能严格说“参数有 95% 概率落在其中”。贝叶斯可信区间则是在给定先验、模型和数据以后，直接对参数作概率陈述。
在应用中，贝叶斯区间的解释更接近很多人的直觉，但它依赖先验和模型设定。因此，理解区间时需要同时关注先验选择，并在必要时做敏感性分析。
16.6 Normal–Normal 模型
再看一个连续变量例子。设 \(y_1,\ldots,y_n\) 是家庭月消费支出，近似满足
\[ y_i\mid\mu\sim N(\mu,\sigma^2), \]
其中 \(\sigma^2\) 暂时视为已知，\(\mu\) 是总体平均消费。选择正态先验
\[ \mu\sim N(\mu_0,\tau_0^2). \]
则后验仍为正态分布：
\[ \mu\mid y \sim N(\mu_n,\tau_n^2), \]
其中
\[ \tau_n^2 = \left(\frac{1}{\tau_0^2}+\frac{n}{\sigma^2}\right)^{-1}, \]
\[ \mu_n = \tau_n^2 \left( \frac{\mu_0}{\tau_0^2} + \frac{n\bar y}{\sigma^2} \right). \]
这里 \(\bar y\) 是样本均值。后验均值是先验均值和样本均值的加权平均，权重由各自的不确定性决定。
set.seed(2026)
n_cons <- 80
y <- rnorm(n_cons, mean = 5.8, sd = 1.2)
mu0 <- 5
tau0 <- 1
sigma <- 1.2
tau_n2 <- 1 / (1 / tau0^2 + n_cons / sigma^2)
mu_n <- tau_n2 * (mu0 / tau0^2 + n_cons * mean(y) / sigma^2)
c(sample_mean = mean(y),
posterior_mean = mu_n,
posterior_sd = sqrt(tau_n2),
q025 = qnorm(0.025, mu_n, sqrt(tau_n2)),
q975 = qnorm(0.975, mu_n, sqrt(tau_n2)))
#> sample_mean posterior_mean posterior_sd q025 q975
#> 5.697 5.685 0.133 5.424 5.946这个例子也展示了贝叶斯更新的直觉：样本量越大，后验越集中；先验越不确定，数据权重越大。
16.7 正态模型中的方差推断
课件中另一个重要例子是正态模型的方差推断。设 \(d_1,\ldots,d_n\) 是某个预测误差或得分差异，近似满足
\[ d_i\mid\sigma^2\sim N(0,\sigma^2), \]
其中均值 0 暂时视为已知，未知参数是方差 \(\sigma^2\)。似然函数为
\[ L(\sigma^2) \propto (\sigma^2)^{-n/2} \exp\left\{-\frac{\sum_{i=1}^n d_i^2}{2\sigma^2}\right\}. \]
如果使用常见的非信息先验
\[ p(\sigma^2)\propto \frac{1}{\sigma^2}, \]
则后验为
\[ p(\sigma^2\mid d) \propto (\sigma^2)^{-n/2-1} \exp\left\{-\frac{v}{2\sigma^2}\right\}, \quad v=\sum_{i=1}^n d_i^2. \]
等价地，精度参数 \(1/\sigma^2\) 的后验可以写成
\[ \frac{1}{\sigma^2}\mid d \sim \frac{\chi_n^2}{v}. \]
因此可以直接从 \(\chi^2\) 分布模拟后验样本。下面用一组模拟预测误差演示：
set.seed(17)
d <- rnorm(80, mean = 0, sd = 1.8)
n_d <- length(d)
v <- sum(d^2)
precision_draw <- rchisq(10000, df = n_d) / v
sigma_draw <- sqrt(1 / precision_draw)
c(sample_sd = sd(d),
posterior_median_sigma = median(sigma_draw),
q025 = unname(quantile(sigma_draw, 0.025)),
q975 = unname(quantile(sigma_draw, 0.975)))
#> sample_sd posterior_median_sigma q025
#> 1.994 1.989 1.720
#> q975
#> 2.338这个例子说明，即使后验密度看起来不像常见的均值公式，也可以通过变量变换和随机模拟来总结后验。
16.8 均值和方差同时未知
如果正态总体的均值和方差都未知，设
\[ y_i\mid\mu,\sigma^2\sim N(\mu,\sigma^2), \]
并使用非信息先验
\[ p(\mu,\sigma^2)\propto \frac{1}{\sigma^2}, \]
则联合后验可以分解为两个更容易处理的部分：
\[ \sigma^2\mid y \sim \frac{S}{\chi_{n-1}^2}, \quad S=\sum_{i=1}^n (y_i-\bar y)^2, \]
以及
\[ \mu\mid\sigma^2,y\sim N\left(\bar y,\frac{\sigma^2}{n}\right). \]
因此，模拟联合后验可以分两步进行：先抽 \(\sigma^2\)，再在给定 \(\sigma^2\) 的条件下抽 \(\mu\)。
set.seed(18)
run_time <- rnorm(20, mean = 300, sd = 35)
n_time <- length(run_time)
ybar <- mean(run_time)
S <- sum((run_time - ybar)^2)
sigma2_draw <- S / rchisq(10000, df = n_time - 1)
mu_draw <- rnorm(10000, mean = ybar, sd = sqrt(sigma2_draw / n_time))
rbind(
mu = quantile(mu_draw, c(0.025, 0.5, 0.975)),
sigma = quantile(sqrt(sigma2_draw), c(0.025, 0.5, 0.975))
)
#> 2.5% 50% 97.5%
#> mu 283.61 303.63 324.09
#> sigma 32.41 43.63 62.47在贝叶斯推断中，“总结后验”往往比“给出一个点估计”更重要。我们可以报告后验均值、中位数、可信区间，也可以计算类似 \(P(\mu>300\mid y)\) 这样的概率。
16.9 先验选择与敏感性分析
先验不是随意填写的数字。常见先验来源包括：
- 以往调查或历史数据；
- 专家经验和制度背景；
- 为稳定计算设置的弱信息先验；
- 为表达明确约束设置的结构性先验。
如果结论对先验非常敏感，就应该诚实报告。下面比较 Beta–Binomial 例子中不同先验下的后验均值。
priors <- data.frame(a = c(1, 2, 10, 30),
b = c(1, 8, 10, 70))
priors$prior_mean <- priors$a / (priors$a + priors$b)
priors$post_mean <- (priors$a + x_sleep) /
(priors$a + priors$b + n_sleep)
priors
#> a b prior_mean post_mean
#> 1 1 1 0.5 0.4138
#> 2 2 8 0.2 0.3514
#> 3 10 10 0.5 0.4468
#> 4 30 70 0.3 0.3228当先验信息量很大且与数据差异明显时，后验会受到更强影响。对于本科课程，重要的是理解这种影响，而不是把“客观先验”和“主观先验”的哲学争论展开得过深。
16.10 后验总结的计算路线
本章中的例子大多有解析后验，因而可以直接使用 qbeta()、rbeta()、qnorm()、rnorm() 等函数。实际贝叶斯模型未必如此简单。面对一个后验分布，常见路线包括：
- 如果后验是熟悉分布，直接使用 R 的分布函数；
- 如果参数只有一两个，可以用网格近似后验并可视化；
- 如果后验在众数附近近似正态，可以使用 Laplace 近似等确定性工具；
- 如果模型更复杂，则使用 Monte Carlo、重要性抽样或 MCMC 总结后验。
后续几章会沿着这几条路线展开。核心问题始终是同一个：给定后验分布，我们如何稳定、透明、可复现地计算出需要的后验摘要。