拟合优度检验

拟合优度检验

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

在数据分析中，我们常常需要判断一组观测数据是否服从某个特定的理论分布。例如，产品的尺寸是否服从正态分布？遗传性状的比例是否符合孟德尔定律？拟合优度检验正是为此而生的强大工具。本节我们将学习由卡尔·皮尔森提出的卡方拟合优度检验，它通过比较观测频数与理论期望频数的差异，来评估数据与理论分布的“拟合”程度。掌握该方法，你将能对分类数据的分布假设进行严格的统计检验，并能将其推广到连续型分布的检验场景。

💡 核心导读

• 核心思想：通过比较观测频数与理论期望频数的差异，构造卡方统计量来评估数据与理论分布的拟合程度。
• 基本形式：针对分类数据，当理论分布完全已知时，检验统计量服从自由度为（类别数-1）的卡方分布。
• 关键推广：当理论分布含有未知参数时，需先用极大似然估计法估计参数，再计算期望频数，此时检验统计量的自由度需减去待估参数的个数。
• 应用场景：从孟德尔豌豆实验的遗传比例检验，到工业产品质量（如滚珠直径）的正态性检验。
• 历史注记：了解皮尔森与费舍尔关于自由度修正的著名争论，体会统计思想的发展脉络。

卡方拟合优度检验的基本思想

拟合优度检验，顾名思义，是检验一个总体的分布与某个指定的理论分布之间“拟合”得有多好。其开山之作是卡尔·皮尔森（Karl Pearson）于1900年提出的卡方拟合优度检验。

它的基本逻辑非常直观：

1. 观测频数：对于分类数据，我们实际观测到每一类别的样本数量，记为 $O_i$（或 $n_i$）。
2. 期望频数：根据我们想要检验的理论分布，可以计算出在假设成立时，每一类别“应该”出现的样本数量，记为 $E_i$。
3. 比较差异：如果总体确实服从该理论分布，那么观测频数 $O_i$ 与期望频数 $E_i$ 之间的差异应该很小。
4. 构造统计量：为了综合衡量所有类别上的差异，并考虑不同类别期望频数大小的影响，皮尔森构造了如下统计量： $$\chi^2 = \sum_{i=1}^{r} \frac{(O_i - E_i)^2}{E_i}$$ 其中，$r$ 为类别总数。这个统计量被称为皮尔森卡方统计量。
5. 做出决策：显然，$\chi^2$ 值越小，说明观测与期望越接近，越支持原假设（即数据服从理论分布）；$\chi^2$ 值越大，则越倾向于拒绝原假设。

那么，这个统计量在原假设成立时服从什么分布呢？皮尔森证明了一个重要定理：在原假设成立且样本量足够大的条件下，上述统计量近似服从自由度为 $r-1$ 的卡方分布，即 $\chi^2 \sim \chi^2(r-1)$。

分类数据的拟合优度检验

我们首先讨论理论分布完全已知的情形，这在历史上是拟合优度检验的起点。

问题的一般形式

设总体按某种属性被分为 $r$ 个互斥的类别 $A_1, A_2, \dots, A_r$。我们关心各类别所占的比例。原假设为： $$H_0: P(A_i) = p_i, \quad i=1,2,\dots,r$$ 其中 $p_i$ 是已知的常数，且满足 $\sum_{i=1}^r p_i = 1$。

现从总体中抽取一个容量为 $n$ 的简单随机样本。记 $O_i$（或 $n_i$）为样本中属于类别 $A_i$ 的观测频数。在原假设 $H_0$ 下，属于类别 $A_i$ 的期望频数为 $E_i = n p_i$。

检验步骤

1. 提出假设：$H_0: P(A_i) = p_i, \quad i=1,2,\dots,r$；$H_1$：至少对某个 $i$，$P(A_i) \neq p_i$。
2. 构造统计量：计算皮尔森卡方统计量 $$\chi^2 = \sum_{i=1}^{r} \frac{(O_i - n p_i)^2}{n p_i}$$
3. 确定分布：在 $H_0$ 成立且 $n$ 较大时，$\chi^2 \sim \chi^2(r-1)$。
4. 拒绝域：给定显著性水平 $\alpha$，拒绝域为 $\chi^2 > \chi^2_{\alpha}(r-1)$，其中 $\chi^2_{\alpha}(r-1)$ 是自由度为 $r-1$ 的卡方分布的 $1-\alpha$ 分位数。

经典案例：孟德尔的豌豆实验

孟德尔将豌豆按颜色（黄/绿）和形状（圆/皱）分为四类：黄圆、绿圆、黄皱、绿皱。他观测了 $n=556$ 粒豌豆，得到观测频数为：$O_1=315$， $O_2=108$， $O_3=101$， $O_4=32$。根据遗传学理论，这四类的理论比例应为 $9:3:3:1$，即 $p_1=9/16, p_2=3/16, p_3=3/16, p_4=1/16$。

检验：原假设 $H_0$：比例符合 $9:3:3:1$。 计算期望频数： $E_1 = 556 \times 9/16 = 312.75$， $E_2 = 556 \times 3/16 = 104.25$， $E_3 = 104.25$， $E_4 = 556 \times 1/16 = 34.75$。 计算卡方统计量： $$\chi^2 = \frac{(315-312.75)^2}{312.75} + \frac{(108-104.25)^2}{104.25} + \frac{(101-104.25)^2}{104.25} + \frac{(32-34.75)^2}{34.75} \approx 0.47$$ 在 $\alpha=0.05$ 水平下，查表得 $\chi^2_{0.05}(4-1) = \chi^2_{0.05}(3) = 7.815$。由于 $0.47 < 7.815$，故不拒绝 $H_0$，认为孟德尔的观测数据与 $9:3:3:1$ 的遗传规律是吻合的。

import numpy as np
from scipy import stats

# 孟德尔豌豆实验数据
observed = np.array([315, 108, 101, 32])
n = observed.sum()
# 理论比例
theoretical_props = np.array([9/16, 3/16, 3/16, 1/16])
# 计算期望频数
expected = n * theoretical_props

print(f"观测频数: {observed}")
print(f"期望频数: {expected}")

# 手动计算卡方统计量
chi2_manual = ((observed - expected)**2 / expected).sum()
print(f"手动计算卡方值: {chi2_manual:.4f}")

# 使用scipy函数验证
chi2_stat, p_value = stats.chisquare(f_obs=observed, f_exp=expected)
print(f"scipy计算卡方值: {chi2_stat:.4f}")
print(f"P值: {p_value:.4f}")

# 临界值比较
alpha = 0.05
df = len(observed) - 1
critical_value = stats.chi2.ppf(1 - alpha, df)
print(f"显著性水平 {alpha} 下的临界值 (df={df}): {critical_value:.4f}")
print(f"结论: 由于卡方值({chi2_stat:.4f}) {'<' if chi2_stat < critical_value else '>='} 临界值({critical_value:.4f})，"
      f"{'不拒绝' if chi2_stat < critical_value else '拒绝'}原假设。")

连续型分布的拟合优度检验

对于连续型总体 $X$，我们想检验其分布函数 $F(x)$ 是否等于某个完全已知的分布函数 $F_0(x)$，即 $H_0: F(x) = F_0(x)$。

由于观测值是连续的，无法直接套用分类数据的检验方法。皮尔森的巧妙思路是：将连续样本空间离散化。

1. 划分区间：将实数轴划分为 $r$ 个互不相交的区间 $I_1, I_2, \dots, I_r$。
2. 转化为分类问题：统计样本落在每个区间 $I_i$ 内的观测频数 $O_i$。
3. 计算理论概率：在原假设 $H_0$ 下，样本落在 $I_i$ 内的理论概率为 $p_i = P_{H_0}(X \in I_i) = F_0(b_i) - F_0(a_i)$，其中 $(a_i, b_i]$ 为区间 $I_i$。
4. 应用卡方检验：此时问题转化为检验“观测频数 $(O_1, \dots, O_r)$ 是否来自多项分布 $(n; p_1, \dots, p_r)$”。检验统计量为： $$\chi^2 = \sum_{i=1}^{r} \frac{(O_i - n p_i)^2}{n p_i}$$ 在 $H_0$ 成立时，它近似服从自由度为 $r-1$ 的卡方分布。

重要注意事项：

1. 信息损失：这种检验实际上检验的是“数据落在各区间内的概率是否等于理论概率”，而非严格检验整个分布函数 $F(x) = F_0(x)$。不同的 $F(x)$ 可能在被划分的区间上产生相同的概率 $p_i$，从而导致检验失效。因此，分点的选择会影响检验的功效。
2. 样本量要求：为保证卡方近似分布的有效性，通常要求每个区间的期望频数 $E_i = n p_i$ 不小于5，且总样本量 $n$ 不小于30。
3. 分点选择：分点不宜过少（损失信息）也不宜过多（导致期望频数过小）。常见的做法是使每个区间的期望频数大致相等，或基于样本分位数进行划分。

含有未知参数的拟合优度检验

实际问题中，我们更常遇到的是检验总体是否服从某一分布族，但该分布族的参数未知。例如，检验数据是否服从正态分布 $N(\mu, \sigma^2)$，但 $\mu$ 和 $\sigma^2$ 未知。

此时，$F_0(x)$ 的形式已知（如正态分布），但依赖于 $k$ 个未知参数 $\theta = (\theta_1, \dots, \theta_k)$，即 $H_0: F(x) = F_0(x; \theta)$。

错误的做法与历史的修正

皮尔森在1900年的论文中意识到了这个问题，并建议：先用样本估计出参数 $\hat{\theta}$，然后用 $F_0(x; \hat{\theta})$ 代替 $F_0(x)$ 来计算理论概率 $p_i(\hat{\theta})$，最后代入卡方统计量： $$\chi^2 = \sum_{i=1}^{r} \frac{(O_i - n p_i(\hat{\theta}))^2}{n p_i(\hat{\theta})}$$ 他最初认为，这个统计量仍然渐近服从 $\chi^2(r-1)$。

然而，费舍尔（R.A. Fisher）在20世纪20年代指出这是一个错误。他证明，当使用参数的极大似然估计（MLE） 时，上述统计量的极限分布是自由度为 $r - 1 - k$ 的卡方分布，其中 $k$ 是待估参数的个数。自由度减少了 $k$，这是因为我们用数据估计了参数，从而“消耗”了一部分自由度。

这个修正具有历史意义，也体现了统计学的严谨性。并非所有的估计方法都适用，但极大似然估计可以满足要求。

检验步骤（以MLE估计参数为例）

1. 提出假设：$H_0$：总体服从分布 $F_0(x; \theta)$，其中 $\theta$ 为 $k$ 维未知参数。
2. 估计参数：基于样本 $X_1, \dots, X_n$，求出 $\theta$ 的极大似然估计 $\hat{\theta}_{MLE}$。
3. 划分区间：划分 $r$ 个区间（通常基于样本或固定分点）。
4. 计算理论概率：计算 $p_i(\hat{\theta}${MLE}) = P{H_0}(X \in I_i | \theta = \hat{\theta}_{MLE})。
5. 构造统计量：计算 $\chi^2 = \sum_{i=1}^{r} \frac{(O_i - n p_i(\hat{\theta}${MLE}))^2}{n p_i(\hat{\theta}{MLE})}。
6. 确定分布与拒绝域：在 $H_0$ 成立下，$\chi^2$ 近似服从 $\chi^2(r-1-k)$。拒绝域为 $\chi^2 > \chi^2_{\alpha}(r-1-k)$。

案例：色盲与性别的遗传模型检验

某地区抽取1000人，按性别（男/女）和色盲（是/否）分为四类，观测频数为：正常男(442)，色盲男(38)，正常女(514)，色盲女(6)。遗传学模型规定各类别的概率为： $p_1 = \frac{1}{2}p$, $p_2 = \frac{1}{2}q$, $p_3 = \frac{1}{2}p^2+pq$, $p_4 = \frac{1}{2}q^2$，其中 $p+q=1$，$p$ 为正常等位基因频率（$q$ 为色盲等位基因频率），是未知参数。

检验：

1. $H_0$：数据服从上述遗传模型。
2. 基于多项分布的似然函数，可求得 $p$ 的MLE为 $\hat{p}=0.91$（$\hat{q}=0.09$）。
3. 代入公式计算理论概率：$\hat{p}_1, \hat{p}_2, \hat{p}_3, \hat{p}_4$。
4. 计算卡方统计量 $\chi^2 \approx 3.056$。
5. 自由度 $df = r - 1 - k = 4 - 1 - 1 = 2$。在 $\alpha=0.05$ 水平下，$\chi^2_{0.05}(2)=5.991$。
6. 由于 $3.056 < 5.991$，不拒绝 $H_0$，认为该地区数据符合遗传学规律。

案例：滚珠直径的正态性检验

某工厂随机抽取50个滚珠测量直径（数据略），问直径是否服从正态分布。

1. $H_0$：直径 $X \sim N(\mu, \sigma^2)$，$\mu, \sigma^2$ 未知。
2. 用样本计算 $\mu, \sigma^2$ 的MLE（即样本均值和样本方差）。
3. 将数据范围划分为7个区间，统计观测频数 $O_i$。
4. 用估计出的 $\hat{\mu}, \hat{\sigma}$ 计算正态分布下落入各区间的理论概率 $\hat{p}_i$ 和期望频数 $E_i = 50\hat{p}_i$。
5. 计算卡方统计量 $\chi^2 \approx 1.7284$。
6. 自由度 $df = r - 1 - k = 7 - 1 - 2 = 4$。$\chi^2_{0.05}(4)=9.488$。
7. 由于 $1.7284 < 9.488$，不拒绝 $H_0$，即没有证据表明滚珠直径不服从正态分布。

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# 模拟滚珠直径数据（假设已知其服从正态分布）
np.random.seed(321)
true_mu, true_sigma = 15.0, 0.5
sample_data = np.random.normal(true_mu, true_sigma, 50)
sample_data.sort()

print(f"样本均值: {sample_data.mean():.4f}, 样本标准差: {sample_data.std(ddof=0):.4f}")

# 步骤1&2：估计参数（MLE）
mu_hat = sample_data.mean()
sigma_hat = sample_data.std(ddof=0)  # ddof=0 给出的是MLE估计量

# 步骤3：划分区间（这里使用样本的近似分位数划分7个区间）
percentiles = np.linspace(0, 100, 8)  # 生成0%, 14.3%, 28.6%, ..., 100%
cut_points = np.percentile(sample_data, percentiles)
# 确保第一个和最后一个点包含所有数据
cut_points[0] = -np.inf
cut_points[-1] = np.inf

print(f"划分区间端点: {cut_points[1:-1]}")

# 计算观测频数
observed, _ = np.histogram(sample_data, bins=cut_points)

# 步骤4：计算理论概率和期望频数
# 计算每个区间在估计的正态分布下的概率
expected_probs = np.diff(stats.norm.cdf(cut_points, loc=mu_hat, scale=sigma_hat))
expected = 50 * expected_probs

print(f"观测频数: {observed}")
print(f"期望频数: {expected}")

# 步骤5：计算卡方统计量（手动）
chi2_stat = ((observed - expected)**2 / expected).sum()
print(f"\n卡方统计量: {chi2_stat:.4f}")

# 步骤6：确定分布与决策
k = 2  # 估计了mu和sigma两个参数
r = len(observed)
df = r - 1 - k
alpha = 0.05
critical_value = stats.chi2.ppf(1 - alpha, df)

print(f"自由度 df = {r} - 1 - {k} = {df}")
print(f"在 α={alpha} 水平下的临界值: {critical_value:.4f}")

if chi2_stat > critical_value:
    print(f"结论: 拒绝原假设 (χ²={chi2_stat:.4f} > {critical_value:.4f})")
else:
    print(f"结论: 不拒绝原假设 (χ²={chi2_stat:.4f} <= {critical_value:.4f})")

# 可视化：观测 vs 期望
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
# 左图：频数对比
x_pos = np.arange(1, r+1)
width = 0.35
ax[0].bar(x_pos - width/2, observed, width, label='观测频数', alpha=0.8)
ax[0].bar(x_pos + width/2, expected, width, label='期望频数', alpha=0.8)
ax[0].set_xlabel('区间编号')
ax[0].set_ylabel('频数')
ax[0].set_title('观测频数与期望频数对比')
ax[0].legend()
ax[0].grid(True, linestyle='--', alpha=0.5)

# 右图：QQ图（另一种正态性检验的直观方法）
stats.probplot(sample_data, dist="norm", plot=ax[1])
ax[1].set_title('正态QQ图')
plt.tight_layout()
plt.show()

📝 动手练一练

1. 分类数据检验：某品牌声称其生产的糖果中，四种口味（草莓、柠檬、橙子、葡萄）的比例是 3:2:2:1。消费者协会随机抽取了 120 颗糖果进行检验，得到四种口味的数量分别为：48, 28, 26, 18。在 0.05 的显著性水平下，检验该品牌的比例声称是否可信。 参考答案： 原假设 $H_0$：比例 $p_1=3/8, p_2=2/8, p_3=2/8, p_4=1/8$。 期望频数 $E = 120 \times [3/8, 2/8, 2/8, 1/8] = [45, 30, 30, 15]$。 计算卡方值 $\chi^2 = \frac{(48-45)^2}{45} + \frac{(28-30)^2}{30} + \frac{(26-30)^2}{30} + \frac{(18-15)^2}{15} \approx 0.2 + 0.133 + 0.533 + 0.6 = 1.466$。 自由度 $df=3$，$\chi^2_{0.05}(3)=7.815$。由于 $1.466 < 7.815$，不拒绝 $H_0$，认为品牌声称可信。

2. 含未知参数的检验：为检验某电子元件寿命（单位：小时）是否服从指数分布 $Exp(\lambda)$，随机抽取 80 个元件进行测试。将寿命数据分为 5 个区间后，观测频数为：[38, 22, 11, 5, 4]。已知指数分布的参数 $\lambda$ 需用样本均值倒数（即MLE）估计，若估计出的 $\hat{\lambda} = 0.02$，请计算用于拟合优度检验的卡方统计量，并指出其近似分布的自由度。 参考答案： 首先，用 $\hat{\lambda}=0.02$ 计算指数分布下落入各区间的理论概率（需已知区间端点，假设为 $[0, 30), [30, 60), [60, 90), [90, 120), [120, \infty)$）。 则 $p_1 = 1 - e^{-0.02$30} \approx 0.4512p1​=1−e−0.02∗30≈0.4512，类似计算 $p_2, p_3, p_4$，$p_5 = e^{-0.02$120} \approx 0.0907。 期望频数 $E_i = 80 \times p_i$。 卡方统计量 $\chi^2 = \sum \frac{(O_i - E_i)^2}{E_i}$。 自由度 $df = r - 1 - k = 5 - 1 - 1 = 3$。（$k=1$ 因为估计了一个参数 $\lambda$）

本章小结

本节深入探讨了统计推断中用于评估分布假设的重要工具——卡方拟合优度检验。

要点回顾：

1. 核心思想：通过比较观测频数与理论期望频数的差异（构造皮尔森卡方统计量）来判断数据是否服从指定分布。
2. 分类数据检验：当理论分布完全已知时，检验统计量 $\chi^2 \sim \chi^2(r-1)$，其中 $r$ 为类别数。
3. 连续分布检验：通过离散化（划分区间）将连续分布检验转化为分类数据检验，但会损失信息且结果依赖于分点选择。
4. 含参分布检验：当理论分布含有 $k$ 个未知参数时，需先用极大似然估计出参数，再计算期望频数。此时检验统计量近似服从 $\chi^2(r-1-k)$ 分布。这是费舍尔对皮尔森原始方法的关键修正。
5. 应用前提：要求样本量足够大，通常每个区间的期望频数不小于5，总样本量不小于30。

行动清单：

1. 动手验证：使用 Python 的 scipy.stats.chisquare 函数，对孟德尔豌豆实验或练习题中的数据重新进行计算，验证卡方值和 P 值，并与临界值法结论对比。
2. 设计实验：尝试用 numpy.random 生成一组来自正态分布 $N(0,1)$ 的样本，然后用本节方法（需估计 $\mu, \sigma$）检验其正态性。再生成一组来自均匀分布或指数分布的样本，用同样的方法检验其正态性，观察检验结果（功效）。
3. 文献阅读：搜索了解“皮尔森-费舍尔争论”的更多历史细节，理解自由度修正背后的统计思想，思考其在现代模型诊断（如逻辑回归的 Hosmer-Lemeshow 检验）中的体现。

— 小象教研组