Shapiro-Wilk与Shapiro-Francia正态性检验:R和Python实战指南
简介这份资源提供复合正态性的 Shapiro-Wilk 参数假设检验实现适用于样本量 3≤n≤5000 的统计分析场景面向需要判断数据是否服从正态分布的研究者、数据分析人员与统计学习者。其核心基于 Royston R94 算法并对 platykurtic低峰态样本额外执行 Shapiro-Francia 正态性检验兼顾两种经典方法的互补性便于在偏态或峰态异常时交叉验证结论。压缩包为 zip 格式仅含 1 个 m 文件体积约 3KB属于轻量级脚本可直接在 MATLAB 环境中调用无需额外依赖。目前已有 946 人学习下载说明该实现具备一定的实用参考价值。读者可获得一套可直接运行的检验脚本理解 Shapiro-Wilk 与 Shapiro-Francia 的算法流程、样本量适用范围及低峰态处理逻辑并据此快速完成正态性诊断与结果解读。1. 正态性检验的两种武器为什么 Shapiro-Wilk 和 Shapiro-Francia 总被一起提跑 t 检验、方差分析、线性回归之前很多人会先做一步正态性检验。R 里shapiro.test()默认走的就是 Shapiro-WilkPython 的scipy.stats.shapiro也是同一套。但如果你翻过旧版统计教材或者某些遗传统计软件的输出会看到另一个名字Shapiro-Francia。它俩不是替代关系而是针对不同样本量区间的互补工具。Shapiro-Wilk 在 n 介于 3 到 5000 时表现稳健Shapiro-Francia 则在小样本到中等样本大约 5 到 100时对尾部偏离更敏感。实际项目中我经常两个都跑一遍对比 p 值和 W 统计量判断数据是否真的偏离正态。这篇文章不扯理论推导只讲怎么在 R 和 Python 里落地这两个检验参数怎么设结果怎么读以及那些让我翻过车的坑。如果你手头有一批数据要决定用参数检验还是非参数检验这篇笔记能直接抄作业。2. Shapiro-Wilk 实操从统计量到 p 值的完整链路2.1 检验原理与选型理由Shapiro-Wilk 的核心思想是把样本顺序统计量与正态分布下顺序统计量的期望值做回归回归的斜率平方就是 W 统计量。W 越接近 1说明样本越像正态分布。它的零假设是“样本来自正态分布”所以 p 值小于显著性水平通常 0.05时拒绝正态性。为什么选它而不是 Kolmogorov-Smirnov因为 K-S 检验对参数估计后的复合假设处理不够精细而 Shapiro-Wilk 专门针对正态性做了优化功效更高。但要注意Shapiro-Wilk 在样本量超过 5000 时scipy会直接报错R 的shapiro.test也会限制在 5000 以内。这时候要么抽样要么换用 Shapiro-Francia 或 Anderson-Darling。另一个选型理由是Shapiro-Wilk 对偏度和峰度都敏感但如果你只关心尾部行为比如极值是否异常Shapiro-Francia 的 W 统计量更直接。我一般会先跑 Shapiro-Wilk 看整体如果 p 值在 0.05 附近徘徊再跑 Shapiro-Francia 交叉验证。2.2 R 语言实现shapiro.test 与参数细节R 基础包里的shapiro.test用起来最简单但有几个参数和限制必须知道。# 生成一组模拟数据正态分布和偏态分布 set.seed(123) normal_data - rnorm(100, mean 50, sd 10) skewed_data - rexp(100, rate 0.5) # Shapiro-Wilk 检验 sw_normal - shapiro.test(normal_data) sw_skewed - shapiro.test(skewed_data) print(sw_normal) print(sw_skewed) # 提取统计量和 p 值 cat(Normal W:, sw_normal$statistic, p-value:, sw_normal$p.value, \n) cat(Skewed W:, sw_skewed$statistic, p-value:, sw_skewed$p.value, \n)逻辑说明shapiro.test只接受一个数值向量返回 W 统计量和 p 值。normal_data的 p 值应该大于 0.05skewed_data的 p 值会远小于 0.05。参数方面R 版本没有额外可调参数但样本量 n 必须满足 3 ≤ n ≤ 5000。如果数据有重复值过多W 统计量的计算会受影响R 会给出警告“ties should not be present for the Shapiro-Wilk test”。这时候要么加微小的随机扰动要么改用其他检验。注意R 的shapiro.test不接受数据框直接输入必须用$提取向量或者用with()。如果数据里有 NA需要加na.rm TRUE但shapiro.test本身没有这个参数得先手动剔除。2.3 Python 实现scipy.stats.shapiro 的边界与替代Python 里scipy.stats.shapiro的接口和 R 类似但返回的是(statistic, p_value)元组。样本量限制同样是 3 到 5000。import numpy as np from scipy import stats np.random.seed(42) normal_data np.random.normal(loc50, scale10, size100) skewed_data np.random.exponential(scale2, size100) # Shapiro-Wilk 检验 stat_normal, p_normal stats.shapiro(normal_data) stat_skewed, p_skewed stats.shapiro(skewed_data) print(fNormal: W{stat_normal:.4f}, p{p_normal:.4f}) print(fSkewed: W{stat_skewed:.4f}, p{p_skewed:.4f}) # 当样本量超过 5000 时scipy 会报错 large_data np.random.normal(size6000) try: stats.shapiro(large_data) except Exception as e: print(Error:, e)逻辑说明stats.shapiro返回的统计量就是 Wp 值越小越拒绝正态性。当 n 5000 时scipy抛出ValueError提示“样本量太大”。这时候常见做法是随机抽样 5000 个点或者改用stats.normaltest基于 DAgostino-Pearson 的偏度峰度检验。但normaltest对样本量更敏感大样本下容易拒绝所以更推荐用 Shapiro-Francia 或者直接看 Q-Q 图。参数方面scipy.stats.shapiro没有nan_policy参数数据里有 NaN 会直接返回 NaN。需要先用np.nan_to_num或者dropna处理。另外shapiro对重复值的容忍度比 R 稍好但重复值超过 20% 时结果也不可靠。3. Shapiro-Francia 实操小样本下的替代方案与手动实现3.1 为什么需要 Shapiro-FranciaShapiro-Francia 是 Shapiro-Wilk 的简化版本用顺序统计量的期望值与样本值的相关系数平方作为统计量 W。它的优势在于计算更简单对小样本n 50的尾部偏离更敏感而且没有 5000 的上限。但它的临界值表不如 Shapiro-Wilk 那么普及很多软件不直接提供。R 的nortest包里有sf.testPython 没有现成函数需要手动实现或者用scipy.stats的shapiro替代。我一般在样本量小于 30 且怀疑有离群值时会优先跑 Shapiro-Francia。3.2 R 中 nortest 包的 sf.test 用法nortest包提供了sf.test用法和shapiro.test几乎一样。# 安装并加载 nortest 包 if (!require(nortest)) { install.packages(nortest) } library(nortest) # 使用之前的数据 sf_normal - sf.test(normal_data) sf_skewed - sf.test(skewed_data) print(sf_normal) print(sf_skewed) # 对比 Shapiro-Wilk 和 Shapiro-Francia 的 p 值 cat(SW p:, sw_normal$p.value, SF p:, sf_normal$p.value, \n) cat(SW p:, sw_skewed$p.value, SF p:, sf_skewed$p.value, \n)逻辑说明sf.test返回的统计量是 Wp 值含义相同。对于正态数据两个检验的 p 值应该都大于 0.05对于偏态数据Shapiro-Francia 的 p 值通常更小说明它更敏感。参数方面sf.test没有额外参数但样本量建议在 5 到 100 之间。如果 n 太小比如小于 5临界值表不准确结果不可靠。注意nortest包里的sf.test对重复值的处理比shapiro.test更宽松但重复值超过 30% 时依然会警告。如果数据有大量重复建议先做频数统计或者改用lillie.testKolmogorov-Smirnov 的 Lilliefors 修正。3.3 Python 手动实现 Shapiro-FranciaPython 没有现成的 Shapiro-Francia 函数但可以手动计算 W 统计量然后用近似公式求 p 值。下面是一个可复现的实现。import numpy as np from scipy import stats def shapiro_francia(data): 手动实现 Shapiro-Francia 正态性检验 参数 data: 一维数组 返回: W 统计量, p 值 n len(data) if n 5: raise ValueError(样本量至少为 5) # 排序 x np.sort(data) # 计算正态顺序统计量的期望值 m_i # 使用 Blom 近似: m_i Phi^{-1}((i - 0.375) / (n 0.25)) i np.arange(1, n 1) m stats.norm.ppf((i - 0.375) / (n 0.25)) # 计算相关系数平方 numerator np.sum(m * x) denominator np.sqrt(np.sum(m**2) * np.sum(x**2)) W_prime (numerator / denominator) ** 2 # 计算 p 值使用 Royston 近似 # 变换到正态分布 mu 0.0038915 * np.log(n)**3 - 0.083751 * np.log(n)**2 - 0.31082 * np.log(n) - 1.5861 sigma np.exp(0.0030302 * np.log(n)**2 - 0.082676 * np.log(n) - 0.4803) z (np.log(1 - W_prime) - mu) / sigma p_value 1 - stats.norm.cdf(z) return W_prime, p_value # 测试 W_prime_normal, p_normal_sf shapiro_francia(normal_data) W_prime_skewed, p_skewed_sf shapiro_francia(skewed_data) print(fNormal SF: W{W_prime_normal:.4f}, p{p_normal_sf:.4f}) print(fSkewed SF: W{W_prime_skewed:.4f}, p{p_skewed_sf:.4f})逻辑说明这段代码先排序然后用 Blom 近似计算正态顺序统计量的期望值再求相关系数平方得到 W。p 值部分用了 Royston 提出的近似公式把 W 变换到标准正态分布。参数方面n必须大于等于 5否则期望值计算不稳定。m的计算公式里 0.375 和 0.25 是 Blom 的经验常数也可以用 0.5 和 0.5 替代但 0.375 在小样本下更准。如果数据有重复值排序后相关系数会偏高导致 W 偏大p 值偏大容易漏掉非正态性。这时候建议先做 jitter 或者用scipy.stats.shapiro交叉验证。注意手动实现的 p 值只是近似和 R 的sf.test结果可能有微小差异。如果对精度要求高建议用 R 的nortest包或者用scipy.stats.shapiro的结果作为参考。4. 避坑与排查正态性检验里那些让我翻车的细节4.1 样本量超过 5000 时 shapiro 报错现象Python 的scipy.stats.shapiro在 n 5000 时抛出ValueErrorR 的shapiro.test也限制在 5000 以内。 原因Shapiro-Wilk 的临界值表只算到 5000超过后统计量的分布不再准确。 解决随机抽样 5000 个点或者改用 Shapiro-Francia无上限或者用stats.normaltest配合 Q-Q 图。我一般会先画 Q-Q 图如果点基本在直线上就不纠结 p 值了。4.2 重复值过多导致警告现象R 的shapiro.test报“ties should not be present”Python 的shapiro返回的 p 值异常大。 原因重复值会让顺序统计量的期望值计算出现偏差W 统计量被高估。 解决如果重复值比例低于 20%可以忽略警告如果超过 20%先做频数合并或者改用lillie.test。我遇到过一批问卷数据1 到 5 的 Likert 量表重复值极多最后直接用了 K-S 检验的 Lilliefors 修正。4.3 p 值在 0.05 附近时的决策困境现象p 0.048 或 p 0.052到底算不算正态 原因p 值只是证据强度不是非黑即白。样本量越大越容易拒绝正态性即使偏离很小。 解决结合 W 统计量和 Q-Q 图。如果 W 0.98 且 Q-Q 图接近直线即使 p 0.05也可以认为近似正态。我通常会把显著性水平设为 0.01减少假阳性。4.4 Shapiro-Francia 的 p 值计算偏差现象手动实现的 Shapiro-Francia 和 R 的sf.test结果不一致。 原因p 值近似公式不同Royston 的公式在小样本下误差较大。 解决以 R 的nortest为准或者用scipy.stats.shapiro的结果做交叉验证。如果两者结论矛盾优先相信 Shapiro-Wilk因为它的临界值表更精确。4.5 数据里有 NaN 或 Inf现象Python 的shapiro返回 NaNR 的shapiro.test报错“missing values”。 原因两个检验都不处理缺失值。 解决先剔除 NaN 和 Inf或者用np.nan_to_num替换。但替换会引入偏差最好直接删除。我一般会在检验前加一行data data[np.isfinite(data)]。5. 进阶技巧用 Q-Q 图和 Monte Carlo 模拟验证检验结果5.1 Q-Q 图比 p 值更直观的判据p 值只能告诉你“是否拒绝”Q-Q 图能告诉你“哪里偏离”。我习惯把 Shapiro-Wilk 和 Shapiro-Francia 的 p 值放在一起再画一张 Q-Q 图。如果 Q-Q 图的尾部点明显偏离直线即使 p 值大于 0.05也要警惕。import matplotlib.pyplot as plt fig, axes plt.subplots(1, 2, figsize(12, 5)) # 正态数据的 Q-Q 图 stats.probplot(normal_data, distnorm, plotaxes[0]) axes[0].set_title(fNormal Data\nSW p{p_normal:.4f}, SF p{p_normal_sf:.4f}) # 偏态数据的 Q-Q 图 stats.probplot(skewed_data, distnorm, plotaxes[1]) axes[1].set_title(fSkewed Data\nSW p{p_skewed:.4f}, SF p{p_skewed_sf:.4f}) plt.tight_layout() plt.show()逻辑说明probplot会画出样本分位数和理论分位数的散点图红线是参考线。正态数据的点应该紧贴红线偏态数据的尾部会明显偏离。参数方面distnorm指定正态分布也可以换成其他分布。如果点呈现 S 形说明尾部厚如果呈现倒 S 形说明尾部薄。5.2 Monte Carlo 模拟验证检验的功效如果你不确定某个样本量下 Shapiro-Wilk 和 Shapiro-Francia 哪个更可靠可以跑一个 Monte Carlo 模拟比较它们的拒绝率。def monte_carlo_power(n, n_sim1000, alpha0.05): 比较 SW 和 SF 在偏态分布下的拒绝率 sw_reject 0 sf_reject 0 for _ in range(n_sim): data np.random.exponential(scale1, sizen) _, p_sw stats.shapiro(data) _, p_sf shapiro_francia(data) if p_sw alpha: sw_reject 1 if p_sf alpha: sf_reject 1 return sw_reject / n_sim, sf_reject / n_sim # 测试不同样本量 for n in [10, 20, 50, 100]: sw_power, sf_power monte_carlo_power(n) print(fn{n}: SW power{sw_power:.3f}, SF power{sf_power:.3f})逻辑说明这段代码生成指数分布的样本分别跑 Shapiro-Wilk 和 Shapiro-Francia统计拒绝零假设的比例。拒绝率越高说明检验越能识别出非正态性。参数方面n_sim是模拟次数1000 次足够稳定alpha是显著性水平。从结果看小样本下 Shapiro-Francia 的拒绝率通常更高但样本量超过 50 后两者差距缩小。注意Monte Carlo 模拟耗时较长建议在 Jupyter Notebook 里跑或者用并行加速。如果只是做一次决策直接看 Q-Q 图更快。5.3 我的固定检查流程从那以后我每次做正态性检验都强制走一遍这个流程先看样本量n 5000 就跑 Shapiro-Wilkn 5000 就抽样或者换 Shapiro-Francia然后检查重复值和缺失值该剔除的剔除接着画 Q-Q 图肉眼确认尾部行为最后如果 p 值在 0.05 附近再跑一次 Monte Carlo 模拟看功效。这套流程帮我省掉了不少后悔药。希望帮到你。本文还有配套的精品资源点击获取