Bayes统计推断

Bayes统计推断

小象实战讲义 · 数据科学的统计基础

在上一节中，我们学习了如何确定参数的先验分布，并利用贝叶斯公式将先验分布与样本信息结合，得到后验概率分布。本节我们将聚焦于贝叶斯统计推断的核心：如何基于后验概率分布对未知参数进行点估计、区间估计和假设检验。学完本节，你将掌握贝叶斯框架下进行统计推断的完整流程，理解其相较于经典（频率学派）方法的优势与直观解释，并能在数据科学实践中灵活运用。

💡 核心导读

• 贝叶斯点估计：基于后验分布，我们可以通过后验均值、后验中位数或后验众数来给出参数的点估计，其中后验均值最为常用。
• 贝叶斯区间估计（可信区间）：给定样本后，参数落在某个区间的后验概率为指定水平（如95%），该区间称为可信区间。其解释比经典置信区间更直观。
• 贝叶斯假设检验：通过计算参数属于原假设和备择假设参数空间的后验概率，并比较其大小，即可做出决策，过程直接明了。
• 核心依赖：所有贝叶斯推断都建立在后验分布之上，而后验分布的计算又依赖于先验分布与样本分布的设定。

贝叶斯点估计

贝叶斯统计推断的基本原则是：一切关于参数的推断都应基于其后验概率分布。参数 $\theta$ 的后验分布 $\pi(\theta \mid \mathbf{x})$ 综合了三种信息：先验信息（$\theta$ 的先验分布）、总体信息（样本分布 $p(\mathbf{x} \mid \theta)$ 的形式）以及样本信息（观测到的数据 $\mathbf{x}$）。

因此，对参数 $\theta$ 进行点估计，自然可以从其后验分布入手。常用的贝叶斯点估计方法有以下三种：

1. 后验众数估计：取后验概率密度函数（或概率质量函数）达到最大值时的 $\theta$ 值作为估计。即： $$\hat{\theta}$${MD} = \arg \max{\theta} \pi(\theta \mid \mathbf{x}) 这被称为众数型贝叶斯估计。

2. 后验中位数估计：取后验分布的中位数作为估计。即满足 $P(\theta \le \hat{\theta}${ME} \mid \mathbf{x}) = 0.5P(θ≤θ^ME​∣x)=0.5 的 $\hat{\theta}${ME}。这被称为中位数型贝叶斯估计。

3. 后验均值估计：取后验分布的期望（均值）作为估计。即： $$\hat{\theta}$${E} = E(\theta \mid \mathbf{x}) = \int{\Theta} \theta , \pi(\theta \mid \mathbf{x}) , d\theta 这被称为期望型贝叶斯估计，也是最常用的一种贝叶斯估计。

对于对称的后验分布（如正态分布），后验众数、中位数和均值三者相等。在后续讨论中，若无特殊说明，我们通常采用后验均值作为 $\theta$ 的贝叶斯估计。

示例1：产品次品率的估计

假设某产品的次品率为 $\theta$。为了估计 $\theta$，我们采用有放回抽样，随机抽取 $n$ 个产品进行检验，发现其中有 $x$ 个次品。

• 样本分布：在 $\theta$ 给定的条件下，次品数 $X$ 服从二项分布 $B(n, \theta)$，即 $p(x \mid \theta) = \binom{n}{x} \theta^x (1-\theta)^{n-x}$。
• 经典估计：经典统计学中的极大似然估计 (MLE) 为 $\hat{\theta}_{MLE} = x/n$。
• 贝叶斯估计：由于对 $\theta$ 的过去信息一无所知，我们采用无信息先验。由于 $\theta \in (0, 1)$，其无信息先验可取为 $(0,1)$ 上的均匀分布 $U(0,1)$，即 $\pi(\theta) = 1$。均匀分布是贝塔分布 $Beta(\alpha, \beta)$ 当 $\alpha=\beta=1$ 时的特例。而二项分布的共轭先验正是贝塔分布，因此后验分布也为贝塔分布： $$\pi(\theta \mid x) \propto p(x \mid \theta) \pi(\theta) = \theta^x (1-\theta)^{n-x} \cdot 1 \propto \theta^{(x+1)-1} (1-\theta)^{(n-x+1)-1}$$ 故 $\theta \mid x \sim Beta(x+1, n-x+1)$。
• 后验均值估计：贝塔分布 $Beta(\alpha, \beta)$ 的均值为 $\alpha/(\alpha+\beta)$。因此，$\theta$ 的贝叶斯估计（后验均值）为： $$\hat{\theta}_{Bayes} = E(\theta \mid x) = \frac{x+1}{(x+1)+(n-x+1)} = \frac{x+1}{n+2}$$

下表比较了在不同样本情况下，极大似然估计与贝叶斯估计的结果：

对比分析：

• 当 $x=n$（全为次品）时，MLE 始终给出 $\hat{\theta}=1$，无法区分样本量 $n=5$ 和 $n=20$ 所带来的不同确信程度。而贝叶斯估计则反映出，样本量越大（$n=20$），我们对“次品率很高”的确信度也越高（估计值 0.955 > 0.857）。
• 当 $x=0$（全为正品）时，MLE 始终给出 $\hat{\theta}=0$。贝叶斯估计则表明，样本量越大（$n=20$），我们对“次品率很低”的确信度也越高（估计值 0.045 < 0.143）。
• 贝叶斯估计通过引入“伪计数”（本例中为先验均匀分布带来的 +1 计数），避免了在样本量较小时做出“绝对肯定”（0或1）的极端估计，结果更为稳健和合理。即使我们使用的是无信息先验，贝叶斯估计也往往能展现出比极大似然估计更好的小样本性质。

示例2：正态总体均值的估计（方差已知）

设 $X_1, \ldots, X_n \stackrel{\text{i.i.d.}}{\sim} N(\theta, \sigma^2)$，其中 $\sigma^2$ 已知。取 $\theta$ 的共轭先验为正态分布 $N(\mu_0, \tau_0^2)$。

• 后验分布：由共轭性可知，后验分布 $\pi(\theta \mid \mathbf{x})$ 仍为正态分布，其均值和方差为： $$\begin{aligned} E(\theta \mid \mathbf{x}) &= \frac{\frac{1}{\tau_0^2} \mu_0 + \frac{n}{\sigma^2} \bar{x}}{\frac{1}{\tau_0^2} + \frac{n}{\sigma^2}} = w \mu_0 + (1-w) \bar{x}, \ \operatorname{Var}(\theta \mid \mathbf{x}) &= \left( \frac{1}{\tau_0^2} + \frac{n}{\sigma^2} \right)^{-1} \triangleq \tau_n^2. \end{aligned}$$ 其中 $\bar{x} = \frac{1}{n}\sum_{i=1}^n x_i$，权重 $w = \frac{1/\tau_0^2}{1/\tau_0^2 + n/\sigma^2}$。
• 贝叶斯估计：由于后验分布是正态的（对称分布），其后验众数、中位数、均值三者相同。因此 $\theta$ 的贝叶斯估计即为后验均值 $\hat{\theta}_{Bayes} = E(\theta \mid \mathbf{x})$。
• 估计的解释：后验均值是先验均值 $\mu_0$ 与样本均值 $\bar{x}$ 的加权平均。权重 $w$ 和 $(1-w)$ 分别由先验精度 $1/\tau_0^2$ 和样本信息总精度 $n/\sigma^2$ 决定。
  • 先验方差 $\tau_0^2$ 越小（先验信息越确定），先验精度 $1/\tau_0^2$ 越大，权重 $w$ 越大，后验估计更倾向于先验均值。
  • 样本量 $n$ 越大或总体方差 $\sigma^2$ 越小（样本信息越精确），样本信息总精度 $n/\sigma^2$ 越大，权重 $(1-w)$ 越大，后验估计更倾向于样本均值。
  • 当样本量 $n \to \infty$ 时，$(1-w) \to 1$，贝叶斯估计收敛于样本均值（即 MLE），这体现了贝叶斯估计的大样本性质。

import numpy as np
import scipy.stats as stats

# 示例：正态总体均值的贝叶斯估计计算
np.random.seed(321)
# 真实参数
true_theta = 5.0
sigma = 2.0  # 已知总体标准差
# 先验分布参数
mu_0 = 4.0
tau_0 = 1.5  # 先验标准差
# 生成样本
n = 10
sample = stats.norm.rvs(loc=true_theta, scale=sigma, size=n)
x_bar = np.mean(sample)
print(f"样本均值 x̄ = {x_bar:.3f}")

# 计算后验分布的参数
# 精度：方差的倒数
prior_precision = 1 / (tau_0 ** 2)
sample_precision = n / (sigma ** 2)
posterior_precision = prior_precision + sample_precision
tau_n_sq = 1 / posterior_precision  # 后验方差
mu_n = (prior_precision * mu_0 + sample_precision * x_bar) / posterior_precision

print(f"先验分布: N({mu_0}, {tau_0**2:.3f})")
print(f"后验分布: N({mu_n:.3f}, {tau_n_sq:.3f})")
print(f"θ的贝叶斯估计 (后验均值) = {mu_n:.3f}")
print(f"经典极大似然估计 (样本均值) = {x_bar:.3f}")

贝叶斯区间估计（可信区间）

在经典统计学中，区间估计的结果是置信区间，其解释为：在重复抽样下，随机区间以 $(1-\alpha)$ 的概率覆盖固定的真实参数 $\theta$。这种基于频率的解释对于单次观测得到的区间而言，不够直观。

贝叶斯区间估计则基于后验分布，概念上更为直接。给定样本 $\mathbf{x}$ 和显著性水平 $\alpha$（通常取 0.05 或 0.01），若存在区间 $[a, b]$ 使得： $$P(a \le \theta \le b \mid \mathbf{x}) = 1 - \alpha,$$ 则称 $[a, b]$ 为参数 $\theta$ 的可信水平为 $1-\alpha$ 的贝叶斯可信区间（Credible Interval）。$a$ 和 $b$ 分别称为可信下限和可信上限。

可信区间与置信区间的本质区别：

1. 解释不同：对于可信区间 $[a, b]$，我们可以直接说“在给定当前样本下，参数 $\theta$ 落在该区间内的概率为 $1-\alpha$”。因为贝叶斯观点中 $\theta$ 是随机变量。而对于置信区间，只能说“这个随机区间以 $1-\alpha$ 的概率覆盖 $\theta$”，因为频率观点中 $\theta$ 是固定常数。
2. 求解难度不同：经典置信区间的构造常依赖于寻找枢轴量（其分布与 $\theta$ 无关），这有时很困难。而贝叶斯可信区间的求解，在得到后验分布后，通常只需找到该分布的 $\alpha/2$ 和 $1-\alpha/2$ 分位数即可，过程更简便。

示例：正态总体均值的可信区间（方差已知）

沿用上一个示例的设定：$X_i \stackrel{\text{i.i.d.}}{\sim} N(\theta, \sigma^2)$，$\sigma^2$ 已知，先验 $\theta \sim N(\mu_0, \tau_0^2)$。后验分布为 $\theta \mid \mathbf{x} \sim N(\mu_n, \tau_n^2)$。

则 $\theta$ 的可信水平为 $1-\alpha$ 的可信区间为： $$\left[ \mu_n - z_{1-\alpha/2} \cdot \tau_n, \quad \mu_n + z_{1-\alpha/2} \cdot \tau_n \right],$$ 其中 $z_{1-\alpha/2}$ 是标准正态分布的 $1-\alpha/2$ 分位数（即上侧 $\alpha/2$ 分位数），$\tau_n = \sqrt{\tau_n^2}$ 是后验标准差。

数值案例：儿童智商测验

假设儿童智商测验得分 $X$ 服从 $N(\theta, 100)$，其中 $\theta$ 为真实智商。根据历史经验，$\theta$ 的先验分布为 $N(100, 225)$。现有一名儿童测得 $x=115$。

• 样本量：$n=1$。
• 计算后验分布：代入公式，其中 $\sigma^2=100$, $\mu_0=100$, $\tau_0^2=225$。 $$\begin{aligned} \tau_n^2 &= \left( \frac{1}{225} + \frac{1}{100} \right)^{-1} = \left( \frac{100+225}{22500} \right)^{-1} = \frac{22500}{325} \approx 69.23, \ \mu_n &= \frac{\frac{1}{225} \cdot 100 + \frac{1}{100} \cdot 115}{\frac{1}{225} + \frac{1}{100}} = \frac{\frac{100}{225} + \frac{115}{100}}{\frac{325}{22500}} \approx 110.38. \end{aligned}$$ 故后验分布为 $\theta \mid x=115 \sim N(110.38, 69.23)$。
• 计算95%可信区间：$z_{0.975} \approx 1.96$，区间为 $110.38 \pm 1.96 \times \sqrt{69.23} \approx 110.38 \pm 16.31$，即 $[94.07, 126.69]$。
• 对比经典95%置信区间：经典方法下，$\theta$ 的置信区间为 $x \pm z_{0.975} \cdot \sigma$，即 $115 \pm 1.96 \times 10 = [95.4, 134.6]$。

对比分析：

1. 精度：贝叶斯可信区间 $[94.07, 126.69]$ 的长度约为 32.62，而经典置信区间 $[95.4, 134.6]$ 的长度约为 39.2。贝叶斯区间更短，精度更高，因为它利用了先验信息（认为智商普遍在100左右）。
2. 解释：对于可信区间，我们可以说“基于本次测验得分115和先验信息，该儿童真实智商 $\theta$ 落在 $[94.07, 126.69]$ 内的概率是95%”。对于置信区间，则只能说“如果用同样的方法对许多儿童进行测验，那么计算出的区间中大约有95%会包含其真实智商”。

import numpy as np
import scipy.stats as stats

# 儿童智商测验案例：可信区间 vs 置信区间
np.random.seed(321)
# 已知条件
x = 115  # 观测值
sigma = 10  # 已知测量标准差
mu_0 = 100  # 先验均值
tau_0 = 15  # 先验标准差 (sqrt(225))
alpha = 0.05

# 1. 贝叶斯方法：计算后验分布及可信区间
n = 1
prior_precision = 1 / (tau_0 ** 2)
sample_precision = n / (sigma ** 2)
posterior_precision = prior_precision + sample_precision
tau_n_sq = 1 / posterior_precision
mu_n = (prior_precision * mu_0 + sample_precision * x) / posterior_precision
tau_n = np.sqrt(tau_n_sq)
z_half_alpha = stats.norm.ppf(1 - alpha/2)
bayes_ci_lower = mu_n - z_half_alpha * tau_n
bayes_ci_upper = mu_n + z_half_alpha * tau_n

print("=== 贝叶斯可信区间 ===")
print(f"后验分布: θ | x ~ N({mu_n:.2f}, {tau_n_sq:.2f})")
print(f"{int((1-alpha)*100)}% 可信区间: [{bayes_ci_lower:.2f}, {bayes_ci_upper:.2f}]")
print(f"区间长度: {bayes_ci_upper - bayes_ci_lower:.2f}\n")

# 2. 经典方法：计算置信区间
classical_ci_lower = x - z_half_alpha * sigma
classical_ci_upper = x + z_half_alpha * sigma
print("=== 经典置信区间 ===")
print(f"{int((1-alpha)*100)}% 置信区间: [{classical_ci_lower:.2f}, {classical_ci_upper:.2f}]")
print(f"区间长度: {classical_ci_upper - classical_ci_lower:.2f}")

贝叶斯假设检验

经典假设检验通过控制第一类错误概率，并寻求在给定第一类错误下第二类错误最小的检验（如 UMP 检验）。过程相对复杂。

贝叶斯假设检验则非常直接。考虑检验问题： $$H_0: \theta \in \Theta_0 \quad vs. \quad H_1: \theta \in \Theta_1,$$ 其中 $\Theta_0 \cup \Theta_1 = \Theta$，且 $\Theta_0 \cap \Theta_1 = \varnothing$。

在获得后验分布 $\pi(\theta \mid \mathbf{x})$ 后，我们只需计算两个后验概率： $$\alpha_0 = P(\theta \in \Theta_0 \mid \mathbf{x}), \quad \alpha_1 = P(\theta \in \Theta_1 \mid \mathbf{x}) = 1 - \alpha_0.$$

决策规则很简单：

• 若 $\alpha_0 > \alpha_1$（或等价地，$\alpha_0 / \alpha_1 > 1$），则接受 $H_0$。
• 若 $\alpha_0 < \alpha_1$（或 $\alpha_0 / \alpha_1 < 1$），则拒绝 $H_0$，接受 $H_1$。

示例：儿童智商检验（续）

沿用上一个儿童智商的案例。现在要检验： $$H_0: \theta \le 100 \quad vs. \quad H_1: \theta > 100.$$ 观测值仍为 $x=115$，后验分布为 $\theta \mid x \sim N(110.38, 69.23)$。

计算后验概率： $$\begin{aligned} \alpha_0 &= P(\theta \le 100 \mid x=115) = \Phi\left( \frac{100 - 110.38}{\sqrt{69.23}} \right) \approx \Phi(-1.25) \approx 0.106, \ \alpha_1 &= 1 - \alpha_0 \approx 0.894. \end{aligned}$$ 其中 $\Phi(\cdot)$ 是标准正态分布函数。

由于 $\alpha_0 \approx 0.106 < 0.894 \approx \alpha_1$，即 $\alpha_0 / \alpha_1 \ll 1$，因此我们拒绝原假设 $H_0$，认为该儿童的真实智商 $\theta$ 大于100。

贝叶斯检验的优势：即使在小样本（本例 $n=1$）下，只要先验信息合理，我们仍然可以进行假设检验，并得到直观的“支持原假设或备择假设的概率”。而经典方法在小样本下可能难以进行或功效很低。

import numpy as np
import scipy.stats as stats

# 贝叶斯假设检验示例：儿童智商是否大于100
np.random.seed(321)
# 使用前面计算的后验分布参数
mu_n = 110.38
tau_n = np.sqrt(69.23)  # 后验标准差

# 计算后验概率
theta_0 = 100
alpha_0 = stats.norm.cdf(theta_0, loc=mu_n, scale=tau_n)
alpha_1 = 1 - alpha_0

print("=== 贝叶斯假设检验 ===")
print(f"原假设 H0: θ ≤ {theta_0}")
print(f"备择假设 H1: θ > {theta_0}")
print(f"后验概率 α0 = P(θ ∈ Θ0 | x) = {alpha_0:.3f}")
print(f"后验概率 α1 = P(θ ∈ Θ1 | x) = {alpha_1:.3f}")
print(f"比值 α0 / α1 = {alpha_0/alpha_1:.3f}")
if alpha_0 > alpha_1:
    print("结论：接受 H0 (α0 > α1)")
else:
    print("结论：拒绝 H0，接受 H1 (α0 < α1)")

📝 动手练一练

1. 产品抽检的贝叶斯估计：某生产线产品的合格率 $\theta$ 未知。采用无信息先验（$Beta(1,1)$）。现抽检15个产品，发现12个合格。 a) 计算 $\theta$ 的后验分布。 b) 给出 $\theta$ 的贝叶斯点估计（后验均值）。 c) 计算 $\theta$ 的95%贝叶斯可信区间。

2. 测量系统的贝叶斯推断：假设某物理量 $\theta$ 的测量值 $X$ 服从 $N(\theta, 0.1^2)$。根据仪器校准报告，$\theta$ 的先验分布为 $N(5.0, 0.2^2)$。现进行一次测量，得到 $x=5.1$。 a) 求 $\theta$ 的后验分布。 b) 对假设 $H_0: 4.95 \le \theta \le 5.05$ 与 $H_1: \theta < 4.95 \text{ 或 } \theta > 5.05$ 进行贝叶斯检验。

参考答案：

1. a) 后验分布为 $Beta(12+1, 15-12+1) = Beta(13, 4)$。 b) 后验均值 = $13/(13+4) \approx 0.7647$。 c) $Beta(13,4)$ 分布的0.025和0.975分位数可由统计软件计算。使用 scipy.stats.beta.ppf([0.025, 0.975], 13, 4) 得到近似区间 $[0.544, 0.927]$。

2. a) 已知 $\sigma^2=0.01$, $\mu_0=5.0$, $\tau_0^2=0.04$, $n=1$, $\bar{x}=x=5.1$。 后验方差：$\tau_n^2 = \left( \frac{1}{0.04} + \frac{1}{0.01} \right)^{-1} = (25 + 100)^{-1} = \frac{1}{125} = 0.008$。 后验均值：$\mu_n = \frac{25 \times 5.0 + 100 \times 5.1}{125} = \frac{125 + 510}{125} = \frac{635}{125} = 5.08$。 故后验分布为 $\theta \mid x \sim N(5.08, 0.008)$。 b) 计算后验概率： $\alpha_0 = P(4.95 \le \theta \le 5.05 \mid x) = \Phi\left(\frac{5.05-5.08}{\sqrt{0.008}}\right) - \Phi\left(\frac{4.95-5.08}{\sqrt{0.008}}\right) \approx \Phi(-0.34) - \Phi(-1.46) \approx 0.366 - 0.072 = 0.294$。 $\alpha_1 = 1 - 0.294 = 0.706$。 由于 $\alpha_0 < \alpha_1$，我们拒绝 $H_0$，认为 $\theta$ 不太可能落在 $[4.95, 5.05]$ 区间内。

本章小结

本节系统介绍了贝叶斯统计推断的三大任务：点估计、区间估计和假设检验。其核心思想高度统一：一切推断均基于后验概率分布。

• 点估计：常用后验均值，它是先验信息与样本信息的加权平均，在小样本下通常比经典估计更稳健。
• 区间估计：得到的是可信区间，其解释——“参数落在区间内的概率”——比经典置信区间更直观、更贴近人们对不确定性的认知。
• 假设检验：通过比较参数属于原假设和备择假设空间的后验概率大小来决策，过程直接，并能给出支持假设的“概率”，而非二元结论。

行动清单

1. 理解核心：反复体会“基于后验分布进行推断”这一贝叶斯统计的核心原则。
2. 掌握计算：熟练推导正态总体（方差已知）在共轭先验下的后验分布，并能据此进行点估计、区间估计和假设检验。
3. 代码实践：运行并理解本节提供的Python代码，尝试修改参数（如先验均值、方差、样本量），观察后验分布及推断结果如何变化，加深对贝叶斯推断过程的感性认识。

— 小象教研组