跳到正文
数据统计代码教程

Python 描述性统计与假设检验：SciPy 实用指南

分析数据的第一步不是建模，而是把分布摸清楚。两个变量均值都是 20，可能一个是「所有人都挤在 20 附近」，另一个是「一半人 5、一半人 35」。后者用均值描述就丢掉了全部信息，后续该选哪种检验也完全不同。

这一篇用 sklearn.datasets.load_wine 走一遍完整流程：先把分布看清楚，再决定用哪种检验，最后读懂输出。该数据集有 178 个葡萄酒样本，13 个化学成分变量，target 是三种植物的类别标签。

描述性统计：先看分位数，再看均值

Section titled “描述性统计：先看分位数，再看均值”

describe() 是成本最低的一步，一次给出计数、均值、标准差、最小值、四分位数和最大值。

import numpy as np
import pandas as pd
import scipy.stats as st
from sklearn.datasets import load_wine


wine = load_wine(as_frame=True).frame
wine.columns = [c.replace(" ", "_") for c in wine.columns]
wine["cultivar"] = wine["target"].map({0: "A", 1: "B", 2: "C"})
wine = wine[["cultivar", "alcohol", "malic_acid", "flavanoids", "proline"]]


print(wine[["alcohol", "malic_acid", "proline"]].describe().round(2))
       alcohol  malic_acid  proline
count   178.00      178.00   178.00
mean     13.00        2.34   746.89
std       0.81        1.12   314.91
min      11.03        0.74   278.00
25%      12.36        1.60   500.50
50%      13.05        1.87   673.50
75%      13.68        3.08   985.00
max      14.83        5.80  1680.00

看 alcohol 这一列：均值 13.00，中位数 13.05，两者几乎重合，分布大致对称。malic_acid 不一样，均值 2.34 明显大于中位数 1.87，是右偏（right-skewed）分布——少数高酸样本把均值拉了上去。它的偏度是 1.04，alcohol 只有 −0.05，这个数字比肉眼看表更直接。

需要单独取值时：

print("mean   =", round(wine["alcohol"].mean(), 4))
print("median =", round(wine["alcohol"].median(), 4))
print("mean   =", round(wine["malic_acid"].mean(), 4))
print("median =", round(wine["malic_acid"].median(), 4))
print("skew   =", round(wine["malic_acid"].skew(), 4))
print(wine["alcohol"].quantile([0, 0.25, 0.5, 0.75, 1]).round(2))
mean   = 13.0006
median = 13.05
mean   = 2.3363
median = 1.865
skew   = 1.0397
0.00    11.03
0.25    12.36
0.50    13.05
0.75    13.68
1.00    14.83
Name: alcohol, dtype: float64

quantile() 接受概率序列，返回的索引就是百分位。分位数比均值耐得住极端值：proline 的均值 746.89 与中位数 673.5 差了 73，说明右尾有少量高值，做组间比较时这个差距会体现为标准差偏大。

按分组算，用 groupby 接 describe()：

print(wine.groupby("cultivar")["alcohol"].describe().round(2))
          count   mean   std    min    25%    50%    75%    max
cultivar
A          59.0  13.74  0.46  12.85  13.40  13.75  14.10  14.83
B          71.0  12.28  0.54  11.03  11.92  12.29  12.52  13.86
C          48.0  13.15  0.53  12.20  12.80  13.16  13.50  14.34

三组的均值是 13.74、12.28、13.15，A 组最高。但标准差只有 0.46 到 0.54，组间差距 1.46 是组内标准差的将近三倍——这个比值决定了后面方差分析的 F 值会有多大。

列多的时候用命名聚合（named aggregation），一次算出多个统计量，并且把结果列名定死：

print(wine.groupby("cultivar").agg(
    n=("alcohol", "size"),
    mean_alcohol=("alcohol", "mean"),
    sd_alcohol=("alcohol", "std"),
    mean_proline=("proline", "mean"),
).round(3))
           n  mean_alcohol  sd_alcohol  mean_proline
cultivar
A         59        13.745       0.462      1115.712
B         71        12.279       0.538       519.507
C         48        13.154       0.530       629.896

写法是 新列名=("原列名", "函数名")。输出的表可以直接写进报告，不需要再改列名——这一点比字典写法省事得多，字典写法返回的列名就是原列名，看不出是均值还是最大值。

标准差与标准误不是一回事

Section titled “标准差与标准误不是一回事”

这两个量经常被混用，但用途完全不同。标准差描述样本内部的离散程度，标准误描述「样本均值」这个估计量的不确定度。

al = wine["alcohol"]
print("pandas std (ddof=1) =", round(al.std(), 6))
print("numpy  std (ddof=0) =", round(float(np.std(al)), 6))
print("pandas sem          =", round(al.sem(), 6))
print("std / sqrt(n)       =", round(al.std() / np.sqrt(len(al)), 6))
print("n =", len(al))
pandas std (ddof=1) = 0.811827
numpy  std (ddof=0) = 0.809543
pandas sem          = 0.060849
std / sqrt(n)       = 0.060849
n = 178

pandas.Series.std() 默认除以 n−1（样本标准差），numpy.std() 默认除以 n（总体标准差）。同一个数组，两个库给出的小数点后第三位就不一样了。多数统计分析要的是样本标准差，用 pandas 的 .std() 是对的；用 numpy 时要显式写 np.std(x, ddof=1)。

.sem() 等于 std() / sqrt(n)，上面两行输出完全一致。论文里写「均值 ± 数值」时，用的是标准误，因为它随样本量减小而减小，反映的是估计精度。误用标准差会让读者以为数据的离散程度比实际小。

缺失值：pandas 默认跳过，但静默

Section titled “缺失值：pandas 默认跳过，但静默”

pandas 处理缺失值的默认行为和 R 相反。R 的 mean() 碰到 NA 直接返回 NA，逼你先决定怎么处理；pandas 默认跳过。静默跳过意味着代码不会报错，也不会提醒你少算了样本。

w = wine.copy()
w.loc[w.index[:10], "flavanoids"] = np.nan


print("缺失个数：", int(w["flavanoids"].isna().sum()))
print("mean() 默认        =", round(w["flavanoids"].mean(), 4))
print("mean(skipna=False) =", w["flavanoids"].mean(skipna=False))
print("nanmean            =", round(float(np.nanmean(w["flavanoids"])), 4))
print(w.groupby("cultivar")["flavanoids"].agg(["size", "count", "mean"]).round(3))
缺失个数： 10
mean() 默认        = 1.9727
mean(skipna=False) = nan
nanmean            = 1.9727
          size  count   mean
cultivar
A           59     49  2.983
B           71     71  2.081
C           48     48  0.781

size 数行数，count 数非缺失值。A 组两者差 10，说明这 10 个缺失全落在 A 组，均值 2.983 是用 49 个样本算出来的，不是 59 个。写报告时说「A 组 59 个样本的类黄酮均值是 2.983」就是错的。

跳过缺失值还隐含一个假设：缺失与观测值本身无关（missing at random）。如果类黄酮含量低的样本正好更容易测失败，跳过就会把均值抬高。判断这一点需要看缺失的分布，方法见 Pandas数据清洗。

正态性检验：shapiro

Section titled “正态性检验：shapiro”

t 检验和方差分析都假设数据来自正态分布，样本量不大时值得查一下。

for col in ["alcohol", "malic_acid", "proline"]:
    W, p = st.shapiro(wine[col])
    ps = f"{p:.4f}" if p >= 0.001 else f"{p:.2e}"
    print(f"{col:11s} W = {W:.4f}   p = {ps}")
alcohol     W = 0.9818   p = 0.0200
malic_acid  W = 0.8888   p = 2.95e-10
proline     W = 0.9312   p = 1.74e-07

原假设是「数据来自正态分布」，所以 p 值大才说明没有证据反对正态性。alcohol 的 p = 0.0200 落在边缘，malic_acid 的 W 只有 0.8888，显著偏离正态。

两点注意事项。第一，shapiro 在样本量大于 5000 时不可靠，大样本下它对微小偏离也敏感，几乎总能拒绝原假设，这时该看 QQ 图而不是 p 值。第二，格式化 p 值时要留意：用 f"{p:.4f}" 打印 2.95e-10 会得到 0.0000，读起来像精确的零。写成 0.0000 不表示 p 等于 0，只表示小于 0.0001。上面用条件表达式在 p 很小时切到科学计数法，避免这个歧义。

t 检验

Section titled “t 检验”

单样本检验问的是「总体均值是否等于某个给定值」。葡萄酒的酒精度均值 13.0006 和 13.0 有区别吗：

r = st.ttest_1samp(wine["alcohol"], popmean=13.0)
lo, hi = r.confidence_interval()
print(f"t = {r.statistic:.4f}, df = {len(wine) - 1}, p = {r.pvalue:.4f}")
print(f"95% 置信区间 = [{lo:.4f}, {hi:.4f}]")
print(f"样本均值 = {wine['alcohol'].mean():.4f}")
t = 0.0102, df = 177, p = 0.9919
95% 置信区间 = [12.8805, 13.1207]
样本均值 = 13.0006

p = 0.9919，没有证据说明总体均值不是 13.0，置信区间 [12.8805, 13.1207] 也包含 13.0。t 统计量只有 0.0102，因为分子（样本均值减假设值）才 0.0006，而标准误是 0.0608。

confidence_interval() 是 SciPy 结果对象的方法，比从 t.ppf() 手工算省事，也避免自由度取错。注意它返回的是元组，可以直接解包。

双样本比较用 ttest_ind，它有两个关键参数：

a = wine.loc[wine["cultivar"] == "A", "alcohol"].to_numpy()
b = wine.loc[wine["cultivar"] == "B", "alcohol"].to_numpy()


s = st.ttest_ind(a, b, equal_var=True)
w_ = st.ttest_ind(a, b, equal_var=False)
print(f"Student : t = {s.statistic:.4f}, df = {s.df:.1f}, p = {s.pvalue:.3g}")
print(f"Welch   : t = {w_.statistic:.4f}, df = {w_.df:.4f}, p = {w_.pvalue:.3g}")
lo, hi = w_.confidence_interval()
print(f"均值 A = {a.mean():.4f} (n={len(a)})")
print(f"均值 B = {b.mean():.4f} (n={len(b)})")
print(f"差值   = {a.mean() - b.mean():.4f}")
print(f"95% 置信区间 = [{lo:.4f}, {hi:.4f}]")
Student : t = 16.4786, df = 128.0, p = 1.96e-33
Welch   : t = 16.7113, df = 127.8471, p = 5.93e-34
均值 A = 13.7447 (n=59)
均值 B = 12.2787 (n=71)
差值   = 1.4660
95% 置信区间 = [1.2924, 1.6396]

equal_var=True 是经典的 Student t 检验，假设两组方差相等，自由度取 n₁+n₂−2 = 128。equal_var=False 是 Welch 校正，不假设方差相等，自由度算出来是 127.8471，不是整数。

SciPy 的默认值是 equal_var=True，这一点和 R 正好相反——R 的 t.test() 默认跑 Welch。直接用 ttest_ind(a, b) 会得到 Student 的结果，两组方差差别大时这个默认值会给出偏窄的置信区间。样本量不等或者方差差异明显时，显式写 equal_var=False。

A 组比 B 组平均高 1.4660 个酒精度单位，p = 5.93e-34。

配对设计（同一个体前后测量、同一批样本的两种测法）要用 ttest_rel：

sleep = pd.DataFrame({
    "subject": list(range(1, 11)) * 2,
    "group": ["1"] * 10 + ["2"] * 10,
    "extra": [0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0.0, 2.0,
              1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4],
})
g1 = sleep.loc[sleep["group"] == "1", "extra"].to_numpy()
g2 = sleep.loc[sleep["group"] == "2", "extra"].to_numpy()


rp = st.ttest_rel(g1, g2)
lo, hi = rp.confidence_interval()
print(f"paired      : t = {rp.statistic:.4f}, df = {len(g1) - 1}, p = {rp.pvalue:.6f}")
print(f"95% 置信区间 = [{lo:.6f}, {hi:.6f}]")
print(f"均值差 = {(g1 - g2).mean():.4f}")
d = g1 - g2
print(f"差值 sd = {d.std(ddof=1):.4f}, sem = {d.std(ddof=1) / np.sqrt(len(d)):.4f}")
print(f"t = 均值差 / sem = {d.mean() / (d.std(ddof=1) / np.sqrt(len(d))):.4f}")
ri = st.ttest_ind(g1, g2)
print(f"independent : t = {ri.statistic:.4f}, df = {ri.df:.1f}, p = {ri.pvalue:.4f}")
paired      : t = -4.0621, df = 9, p = 0.002833
95% 置信区间 = [-2.459886, -0.700114]
均值差 = -1.5800
差值 sd = 1.2300, sem = 0.3890
t = 均值差 / sem = -4.0621
independent : t = -1.8608, df = 18.0, p = 0.0792

数据是 10 位受试者在两种条件下的睡眠时长增量，每位受试者贡献一对观测。

配对检验的 t 统计量等于「差值均值除以差值标准误」：−1.58 / 0.389 = −4.0621，和 ttest_rel 的输出一致。它的自由度是 9（对数减 1），不是 18。同一批数据按独立样本处理，t = −1.8608、p = 0.0792，落在 0.05 边缘。配对检验先减掉了个体间差异带来的噪声，把真实效应从背景波动里提了出来。方法选错，显著的结果会被判成不显著。

ttest_rel 要求两个数组长度相同且顺序一一对应。如果两份数据来自不同的排序，先按样本编号 sort_values 对齐，否则算出来的是两组随机配对的结果，没有意义。

方差齐性：levene

Section titled “方差齐性：levene”

经典 t 检验和方差分析都要求各组方差相等。检验方差齐性用 Levene 检验：

lv = st.levene(a, b)
print(f"Levene W = {lv.statistic:.4f}, p = {lv.pvalue:.4f}")
print(f"方差 A = {a.var(ddof=1):.4f}  方差 B = {b.var(ddof=1):.4f}  比值 = {a.var(ddof=1) / b.var(ddof=1):.4f}")
Levene W = 0.3125, p = 0.5771
方差 A = 0.2136  方差 B = 0.2894  比值 = 0.7379

p = 0.5771，没有证据说明两组方差不齐。但这不等于「两组方差相等」——不拒绝原假设只是证据不足，别把 p > 0.05 读成「假设成立」。

另外，Levene 检验本身也受样本量影响。样本量小时它检不出真实的方差差异，样本量大时又会因为微小差异而拒绝。把「先检验方差齐性、再决定用哪种 t 检验」当成固定流程，会让整个分析的实际错误率失控。多数情况下直接用 Welch（equal_var=False）更稳妥，代价只是自由度不再是整数。

非参数替代：mannwhitneyu 与 wilcoxon

Section titled “非参数替代：mannwhitneyu 与 wilcoxon”

数据明显非正态，或者只有序数尺度时，用基于秩的检验：

mw = st.mannwhitneyu(a, b)
print(f"mannwhitneyu U = {mw.statistic}, p = {mw.pvalue:.3g}")
n1, n2 = len(a), len(b)
print(f"秩双列相关 r = {1 - 2 * mw.statistic / (n1 * n2):.4f}")
wl = st.wilcoxon(g1, g2)
print(f"wilcoxon W = {wl.statistic}, p = {wl.pvalue:.6f}")
mannwhitneyu U = 4079.5, p = 1.67e-20
秩双列相关 r = -0.9477
wilcoxon W = 0.0, p = 0.003906

mannwhitneyu 对应两组独立样本，wilcoxon 对应配对样本，对应关系和 ttest_ind / ttest_rel 一样。

这两个函数检验的是位置移动（location shift），不是均值差。A、B 两组的酒精度结论和 t 检验一致（都显著），但 p 值不同，因为它比较的是秩而不是均值和方差。

SciPy 的 mannwhitneyu 默认 alternative="two-sided"，早期版本默认是单侧，升级后行为变了，看旧代码时留意。默认还会做连续性校正（method="auto"），想在并列值很多时关掉，传 method="asymptotic" 或 method="exact" 显式指定。

U 统计量本身不好解读，转成秩双列相关（rank-biserial correlation）就有界了：r = −0.9477，取值范围 [−1, 1]，绝对值接近 1 表示两组几乎完全分离。报告非参数检验时给出这个效应量，比只报 p 值有信息量。

配对的 wilcoxon 输出 W = 0.0，因为有 9 个差值非零、且全部同号，秩和的下界就是 0。注意有一个差值为 0，wilcoxon 默认丢弃零差值（zero_method="wilcox"），有效样本量变成 9，p = 0.003906 = 2/512。

卡方检验：两个分类变量是否独立

Section titled “卡方检验：两个分类变量是否独立”

两个分类变量之间有没有关联，用 pd.crosstab 建列联表，再交给 chi2_contingency。把酒精度按中位数切成高、低两类，看它和品种的关系：

wine["high_alcohol"] = np.where(wine["alcohol"] > wine["alcohol"].median(), "high", "low")
tab = pd.crosstab(wine["cultivar"], wine["high_alcohol"])
print(tab)
chi2, p, dof, exp = st.chi2_contingency(tab)
print(f"X-squared = {chi2:.4f}, df = {dof}, p-value = {p:.3g}")
n = tab.to_numpy().sum()
print(f"Cramer's V = {np.sqrt(chi2 / (n * (min(tab.shape) - 1))):.4f}")
high_alcohol  high  low
cultivar
A               53    6
B                5   66
C               27   21
X-squared = 90.4222, df = 2, p-value = 2.32e-20
Cramer's V = 0.7127

p 值远小于 0.05，酒精度高低和品种不是独立的。Cramér’s V = 0.7127 是效应量，取值在 0 到 1 之间，0.71 属于强关联。只看 p 值会漏掉这件事：p 值随样本量变小，效应量不随样本量变。样本量足够大时，Cramér’s V = 0.05 也能得到 p < 0.001。

Cramér’s V 的算法是卡方统计量除以「样本量 ×（行列数较小者减 1）」，再开平方。min(tab.shape) 取行列数里较小的那个，3×2 的表减 1 得 1。

chi2_contingency 返回四元组：统计量、p 值、自由度、期望频数。前三个通常用得上，第四个在检查前提假设时才用。

卡方检验有个容易被忽略的前提：每个格子的期望频数最好都不小于 5。把酸度分成四档再看：

wine["acid_level"] = pd.cut(wine["malic_acid"], bins=[0, 1.5, 3.0, 4.0, 10],
                            labels=["very low", "low", "high", "very high"])
tab2 = pd.crosstab(wine["acid_level"], wine["cultivar"])
print(tab2)
chi2b, pb, dofb, expb = st.chi2_contingency(tab2)
print(f"X-squared = {chi2b:.4f}, df = {dofb}, p-value = {pb:.3g}")
print("期望频数：")
print(pd.DataFrame(expb, index=tab2.index, columns=tab2.columns).round(2))
print("最小期望频数 =", round(expb.min(), 3), "   期望频数 < 5 的格子：", int((expb < 5).sum()), "/", expb.size)
cultivar     A   B   C
acid_level
very low     5  27   3
low         47  34  15
high         6   6  18
very high    1   4  12
X-squared = 66.2559, df = 6, p-value = 2.39e-12
期望频数：
cultivar        A      B      C
acid_level
very low    11.60  13.96   9.44
low         31.82  38.29  25.89
high         9.94  11.97   8.09
very high    5.63   6.78   4.58
最小期望频数 = 4.584    期望频数 < 5 的格子： 1 / 12

「very high」品种 C 那一格的期望频数是 4.58，低于 5 的门槛。这种情况下卡方近似不再可靠，应该改用 Fisher 精确检验，或者把相邻类别合并。

SciPy 不会像 R 的 chisq.test() 那样自动给出警告。R 会打印 “Chi-squared approximation may be incorrect”，SciPy 只把期望频数放在返回值里，需要自己检查。写一个 (exp < 5).any() 的判断放进分析脚本，比事后返工便宜。

2×2 表还有连续性校正的差异：

wine["acid2"] = pd.cut(wine["malic_acid"], bins=2, labels=["low", "high"])
tab3 = pd.crosstab(wine["acid2"], wine["high_alcohol"])
print(tab3)
c1, p1 = st.chi2_contingency(tab3)[:2]
c0, p0 = st.chi2_contingency(tab3, correction=False)[:2]
print(f"correction=True  : X-squared = {c1:.4f}, p = {p1:.4f}")
print(f"correction=False : X-squared = {c0:.4f}, p = {p0:.4f}")
print("fisher_exact     :", f"p = {st.fisher_exact(tab3).pvalue:.4f}")
high_alcohol  high  low
acid2
low             65   74
high            20   19
correction=True  : X-squared = 0.1011, p = 0.7505
correction=False : X-squared = 0.2493, p = 0.6175
fisher_exact     : p = 0.7173

SciPy 的 correction 参数默认是 True，对 2×2 表自动加 Yates 连续性校正，这一点和 R 一致。校正后统计量从 0.2493 降到 0.1011，p 值从 0.6175 升到 0.7505。自由度大于 1 的表不受影响，校正只对 2×2 生效。

三个结论都不显著，方向也一致。样本量小时（每格期望频数接近或低于 5）优先用 fisher_exact，它给出的是精确概率，不依赖近似。

p 值到底在说什么

Section titled “p 值到底在说什么”

几种常见的误读：

  • p 值不是「原假设为真的概率」，也不是「结果由随机造成的概率」。
  • p < 0.05 不表示效应有 95% 的概率存在。它的含义是：如果原假设成立，出现当前或更极端结果的概率小于 5%。
  • p 值大小和效应大小无关。上面的卡方检验里，p = 2.32e-20 和 Cramér’s V = 0.7127 说的是两件事，前者关于证据强度，后者关于关联强度。
  • p = 0.06 和 p = 0.04 之间没有质的差别，别把一个说成「无效应」、另一个说成「有效应」。
  • f"{p:.4f}" 打印出的 0.0000 不表示 p 等于 0，只是小于 0.0001。

配合置信区间一起读，比盯着 p 值本身可靠。SciPy 的结果对象都有 confidence_interval() 方法，取用成本和打印 p 值差不多。

下一步

Section titled “下一步”

基础检验覆盖的是「两组比较」，但研究里更常见的是三组以上，或者要同时控制多个变量。前者用方差分析，后者用线性回归，两者背后是同一个线性模型框架。

三组以上的均值比较见 /python/modeling/anova，连续变量之间的线性关系见 /python/modeling/linear-regression。画分布图和箱线图先看形状，再决定用哪种检验，方法见 /python/visualization/seaborn。R 语言的 描述性统计与假设检验 是同一套检验的另一种实现，函数名不同，判断标准完全一致。