简介：面向数学建模与美赛的数据预处理需求，压缩包内代码基于格拉布斯准则实现异常值判断，用于识别和修正样本中的极端数据，适合参赛选手或数据分析初学者参考。格拉布斯检验以正态分布为前提，计算最大值与均值的偏差并转换为G值，与临界值Gα比较后判定异常值，α通常取0.05或0.01，代码中对此有明确实现。压缩包共3个文件，含一个MATLAB脚本（.m）、一个MATLAB自动保存的备份文件（.asv）及一个txt说明文档，整体仅1KB，结构精简便于查看核心算法。目前已有106人学习下载，适用于快速上手异常值检测任务。代码覆盖了数据读取、缺失值检查、统计量计算、临界值查找、异常值标记与删除替换处理等完整步骤，可帮助读者将格拉布斯准则直接应用于赛题或实际数据集，避免自行推导公式和编写逻辑，显著提升数据预处理效率。

1. 异常值检测为什么先想到格拉布斯准则：从删错一条数据说起

一组测量数据里躺着一条“明显”偏离的值。绝大多数人第一反应是直接删掉，省事。但真实工程里，一条被误删的数据可能让批次放行、让工艺参数调错方向，代价并不比保留异常值小。格拉布斯准则做的事情很明确：在正态假设下，用样本均值与标准差计算当前极值出现的概率，再和预设置信水平比较，给出“删”或“不删”的统计依据，而不是靠肉眼拍脑袋。网上流传的“基于格拉布斯准则判断异常数据代码.rar”这类压缩包，通常就是把这一整套检验流程封装成可运行脚本，改改数据路径就能用。这篇文章不从某个现成脚本出发，而是把准则拆开讲明白，再给出一套能自己复现和改造的 Python 实现，顺带把五个容易翻车的细节捋清楚。适合样品数量不多、理论上近似正态分布的实验测量与传感器数据，也适合想让报告的异常判断多一句统计依据的从业者。

2. 格拉布斯准则的原理与适用边界：检验统计量、临界值表与样本量下限

2.1 检验统计量怎么算：G 值与临界值的关系

格拉布斯准则的出发点非常朴素：先找出数据里离样本均值最远的那个点，计算它跟均值的距离占样本标准差的比例，这个比例越大，说明该点越“不像”这组数据。对应公式是：

G = (max|xᵢ - x̄|) / s

其中 x̄ 是样本均值，s 是样本标准差。注意 s 要用 n − 1 做分母，也就是 ddof=1 的样本标准差；如果用总体标准差 ddof=0，小样本时 G 值会被明显拉大，临界值判断失真。多数教材和国标里讨论的 G 都基于样本标准差，这一点在写代码时要格外留意。

下一个问题是：G 值多大才算异常？标准做法是把 G 和格拉布斯临界值 G_crit 比较。临界值由样本量 n、显著性水平 α，以及 t 分布上侧分位数计算而来。双侧检验的公式为：

G_crit = ((n − 1) / √n) × √( t²_{α/(2n), n−2} / (n − 2 + t²_{α/(2n), n−2}) )

这里 t_{α/(2n), n−2} 表示自由度为 n−2 的 t 分布的上侧 α/(2n) 分位数。为什么要用 α 除以 2n？因为双侧检验要同时考虑上下两个方向，而同一次检验还要对 n 个样本点做“哪个点最可疑”的筛选，相当于多重比较，需要按 n 放大惩罚。单侧检验则把公式里的 α/(2n) 改成 α/n。

公式可以先不背，落实成表更直观。下面是用 Python 计算出的双侧临界值示例，α 取 0.05 和 0.01：

观察这张表能发现两个特点：样本量越小，临界值越接近 1，也就是在小样本下只有离均值特别远的点才会被判为异常；样本量增大后，临界值缓慢上升，单靠“偏差大于三倍标准差”这种拍脑袋阈值，在 n 很大时会频繁误报。日常做质量分析时，我习惯同时把 G 和 G_crit 一起输出，而不是只给一个布尔判断，这样写检验报告的引用依据时可以直接引用原始值。

2.2 单侧还是双侧：方向没有想清楚，结论会反过来

很多人第一次用格拉布斯准则时忽略了一个前提：检验是双侧的还是单侧的。双侧检验不预设异常方向，既管上端极值也管下端极值，适用于“偏高和偏低都异常”的场景，比如化学实验的平行样浓度、传感器校准值。单侧检验只盯一个方向，适用于业务上只担心超上限或只担心超下限的场景。

选错方向最直接的后果是临界值不一样。以 n = 10、α = 0.05 为例：双侧检验在计算临界值时用 α/(2n)，单侧检验用 α/n，后者的 t 分位数更小，得到的 G_crit 更小，判定阈值更严。也就是说，如果业务上明确“只关心上限”，却选了双侧检验，可能会把一些轻微偏离但不足以判异常的数值放过去；反过来，如果上下限都重要却选了单侧，则会把一侧的正常极值误杀。

落地时我的建议是：不确定方向就默认双侧，这是大多数统计软件和网上示例代码的默认行为；只有在工艺规范、安全标准里明确写了“只允许单向偏差”时，才改成单侧。判断口诀是“业务先说方向，代码再定参数”，不要为了临界值更小而去选单侧。

2.3 适用边界：正态假设与样本量的硬约束

格拉布斯准则最容易被滥用的一点是正态性假设。它的 G 值本质上是正态总体下极差与标准差的比值统计量，数据必须近似来自正态分布。如果原始数据是偏态分布，比如反应时间、浓度这种往往右偏的数据，直接套格拉布斯会频繁把正常的高值误判为异常，因为右偏分布本身就天然拥有更多高值。处理方式有两种：一是先做对数变换或 Box-Cox 变换，把数据拉回近似正态再检验；二是干脆换成不依赖正态假设的方法，比如基于中位数和 MAD（绝对中位差）的离群值识别。常见做法是在业务上把两种方法并行跑一遍，结论一致再剔除。

样本量同样有硬约束。格拉布斯准则要求 n ≥ 3，但 n = 3 或 4 的时候检验几乎没有区分度：三个点里最远的那个只要偏离同伴一点点，就可能被判为“显著异常”，这个结果毫无实际意义。经验上，少于 8 个样本时我不太信任格拉布斯的最终结论，只把它当佐证；真正做剔除决策至少要有 8 到 10 个以上的重复测量值。另外，随着 n 变大，异常点对样本标准差的影响越来越大，这就是后面要详细说的遮蔽效应。

还有一类场景需要提前打个预防针：量化交易策略回测时也会遇到单日极端收益，有人拿格拉布斯去清洗收益率序列。这个方向不是不行，但必须先确认序列不是重尾分布，否则会把正常风险事件当脏数据洗掉，回测结果看着漂亮，实盘立刻打脸。

3. 用 Python 实现格拉布斯检验：从统计量到循环剔除的完整示例代码

3.1 最小复现：5 行代码算出 G 值与可疑点索引

先把最核心的统计量算出来，不急着谈临界值。假设有一组样本数据，来自某次环境监测仪器的 10 次读数：

import numpy as np data = np.array([ 15.2, 14.9, 15.6, 15.1, 16.8, 14.8, 15.3, 15.0, 15.4, 15.2 ]) mean = data.mean() std = data.std(ddof=1) # 样本标准差，分母 n-1 diff = np.abs(data - mean) idx = np.argmax(diff) G = diff[idx] / std print(f"均值: {mean:.4f}") print(f"样本标准差: {std:.4f}") print(f"最可疑点: index={idx}, 值={data[idx]:.4f}") print(f"G统计量: {G:.4f}")

这段代码的运行结果会很直观：16.8 这个读数离均值最远，G 值大约在 1.9 左右。注意ddof=1是必须写的，numpy.std()默认 ddof=0 算总体标准差，偏差会被低估，把 G 值虚高。先运行这一段的目的不是判断异常，而是让你看到每个中间量，尤其是标准差对结果的敏感程度：把样本数据任意一个值改大，标准差变大，G 值未必同步变大，这就是格拉布斯在多个极端点同时存在时容易失灵的根源。

有pandas环境时可以把data换成df['col'].to_numpy()，后续判断直接输出一个标记列。常见做法是不把这段逻辑散落在报告里，而是封装成一个小函数，便于在不同脚本间复用。

3.2 完整实现：查临界值、单次判断与循环剔除

接下来给出一版可以直接放进项目里用的完整函数。它同时处理临界值计算、单次判断和迭代剔除三个动作：

import numpy as np from scipy import stats def grubbs_test(data, alpha=0.05, side='two', max_remove=0.2): """ 基于格拉布斯准则判断并剔除异常数据。 参数： data: 一维数组-like 数据，至少 3 个观测点 alpha: 显著性水平，常用 0.05 side: 'two' 双侧检验，'high' 只查上端，'low' 只查下端 max_remove: 最多剔除比例，防止迭代时误删正常数据 """ arr = np.asarray(data, dtype=float) n0 = len(arr) if n0 < 3: raise ValueError("样本量至少为 3 才有检验意义") if not (0 < max_remove < 0.5): raise ValueError("max_remove 应设置在 0 到 0.5 之间") removed = [] arr_work = arr.copy() max_iter = max(1, int(n0 * max_remove)) for _ in range(max_iter): n = len(arr_work) if n < 3: break mean = arr_work.mean() std = arr_work.std(ddof=1) if std == 0: break # 数据全部相同，不存在异常值 diff = np.abs(arr_work - mean) idx = np.argmax(diff) G = diff[idx] / std if side == 'two': t_crit = stats.t.ppf(1 - alpha / (2 * n), n - 2) elif side in ('high', 'low'): t_crit = stats.t.ppf(1 - alpha / n, n - 2) else: raise ValueError("side 参数只能是 two / high / low") G_crit = (n - 1) / np.sqrt(n) * np.sqrt( t_crit**2 / (n - 2 + t_crit**2) ) if G > G_crit: removed.append({ 'value': arr_work[idx], 'index_in_original': int( np.where(arr == arr_work[idx])[0][0] ), 'G': G, 'G_crit': G_crit, 'n': n }) arr_work = np.delete(arr_work, idx) else: break return { 'cleaned': arr_work, 'removed': removed, 'n_total': n0, 'n_removed': len(removed) }

逻辑说明：每次迭代先计算当前数据集的 n、均值、样本标准差和最大偏差点，再依据 side 用scipy.stats.t.ppf拿 t 分布分位数，换算成 G_crit。G 大于 G_crit 才删除该点，然后进入下一轮；否则直接退出。这里有两个容易被忽略的细节。第一是自由度 n−2：检验统计量依赖均值和标准差两个估计量，删除可疑点后剩余点的独立信息量是 n−2。第二是剔除后的 index 不会自动映射回原始数组，函数里通过np.where(arr == arr_work[idx])找回原始位置，实际使用中这个位置信息比数值本身更重要，因为要回查实验记录。

注意：上面原始索引的查找代码在数据存在重复值时只会返回第一个匹配位置。真实项目里建议先给每条记录加一列自增 row_id，用 row_id 代替数值做定位，避免两个人拿着相同读数时不知道该删的是哪一行。

参数说明：alpha 默认 0.05 是多数质量检验的常规水平；如果处理的是安全相关参数，我会改到 0.01，让剔除更保守。max_remove 的默认值 0.2 是经验值，避免循环剔除把正常数据一批批洗掉——格拉布斯准则对单个异常点设计，多轮迭代后检验水平已经和名义 alpha 不等价，必须用比例硬控。

调用方式和输出如下：

data = np.array([12.3, 12.5, 12.7, 13.1, 15.2]) result = grubbs_test(data, alpha=0.05) print(result['removed']) # 输出示例：[{'value': 15.2, 'index_in_original': 4, 'G': 1.968, 'G_crit': 1.715, 'n': 5}]

看到 G > G_crit 时再确认一次业务上是否说得通，比如 15.2 是否来自仪器未校准时段。统计检验能给你拒绝的理由，但不能替你决定该不该信这个值。

3.3 没有 SciPy 怎么办：查表回退与 Excel 里的手工做法

如果你的环境不允许安装 scipy，临界值可以直接从上一章那种表里线性插值，或者在代码里维护一份关键 n 值的临界值表。工程上我更推荐保留 scipy，因为插值在小样本区间会带来不可忽略的误差，而 t 分布分位数计算是一行代码的事。

Excel 用户在报表里也有替代路径：格拉布斯临界值本质上来自 t 分布，可以在单元格里用T.INV函数配合公式手工算。步骤是：先用AVERAGE、STDEV.S算均值和样本标准差，找最大偏差并算 G 值；再在另一列用 T 分布函数得到分位数，套入临界值公式；最后用IF(G > G_crit, "异常", "正常")输出标记。如果确实想自动化，写一段 VBA 代码也能套同样的 t 分位数公式，但维护成本比 Python 高不少，只适合偶尔复核几列数据的情况。数据量到几十列、上百列，我还是建议用 Python 函数统一跑，结果落成一个 CSV，比在 Excel 里拉公式可复查性好得多。

4. 格拉布斯准则实战避坑：五个让检验结果翻车的场景

4.1 遮蔽效应：两个异常值互相掩护，单次检验直接停手

现象：数据里其实有两个离群点，一个偏大一个偏小，格拉布斯单次检验却报告“未发现异常”，或者只删了较大的那个，剩下那个再也报不出来了。只要数据整体标准差被两个极值拉得很宽，无论删哪个，下一步的统计量都会变得相对“正常”。

原因：G 值用的是全体样本的均值和标准差，异常点本身在抬高标准差。这就是经典的遮蔽效应：极端点越多，它们越“互坑”，单次检验的检出能力越差。

解决：先把 G 值接近但未超过临界值的点做一个排序散点图看整体分布，再用 3.2 的循环剔除函数跑一遍，限制最多剔除 20%。如果图上明显存在两个方向相反的极端点，或者循环第一次剔除后 G 值又快速上升，就基本可以确认遮蔽效应，此时应改用广义 ESD 检验，它显式支持多个异常点，虽然 scipy 没有内置实现，按算法步骤写一遍并不算难。

4.2 单侧双侧选错：业务方向与临界值不匹配

现象：同一组数据，用双侧检验判为正常，换单侧就判为异常，同事拿着两份结论来问你哪份对。这在传感器超限分析里很常见：只关心上限超限，却用了默认的双侧检验，极限值被卡在临界值内侧。

原因：双侧检验把显著性水平按两个方向平分，使用的 t 分位数更大，G_crit 更大，判定更宽松；单侧检验用全量 α，G_crit 更小，判定更严格。选择的关键不是数据长什么样，而是业务上是否接受“两个方向都算异常”。

解决：开工前把业务规则写进脚本注释，例如“此工序只允许正向偏差，使用 side='high'”。拿不准时维持双侧，并在输出里同时保留原始值和被剔除值，方便回查。把方向决策留给最懂业务的人，不要由写代码的人顺手拍板。

4.3 偏态数据硬套正态假设：误报率和想象中完全不同

现象：一组反应时间数据，最小值 2.3 秒，最大值 5.1 秒，格拉布斯检验把 5.1 秒判为异常。但从业务角度这组数据本来就右偏，5.1 秒只是分布的高值尾巴，不是测量错误。

原因：格拉布斯准则的 G 统计量在数据近似正态时才服从推导所用的分布。偏态分布中，极值出现的天然频率更高，检验会把分布的自然高值误认为离群点。

解决：先判断原始数据的形态。最简单的做法是看偏度，scipy.stats.skew绝对值超过 1 的，优先做对数变换再检验；或者换用不依赖正态假设的 MAD 方法。需要注意的是，变换后的结论只能说明“变换空间里的异常”，回写成业务报告时要注明用的是对数空间，否则复核的人会对不上号。

4.4 迭代剔除没有硬上限：好数据被一把把洗掉

现象：循环剔除脚本不设比例上限，跑完后原本 30 条数据只剩 20 条，而且每次删除后剩下的数据越来越“齐”，G 值始终显著。最后得到的一小撮数据确实干净，但已经不是原来的样本。

原因：格拉布斯准则默认假设数据里最多只有一个离群点。多轮检验时每轮显著性水平都会偏移，随着样本量变小，标准差被不断压缩，原本正常但偏边缘的点会逐步被踢出。

解决：给循环加两个硬条件——最大剔除比例不超过 20%，且每轮删除后用标记确认“删的是不是同一实验条件下的重复测量”。我一般会把 20% 当成硬写死的参数，而不是靠人工停。真出现接近 20% 还收不住的情况，检查数据是否混入了不同批次，这是质量问题，不是统计问题。

4.5 下载的 rar 代码包跑不通：先分清是算法问题还是环境问题

现象：从一个分享链接里下载了“基于格拉布斯准则判断异常数据代码.rar”，解压后双击运行，报 ModuleNotFoundError，或者中文注释变成乱码，甚至解压时提示需要密码。此时第一反应不应该是怀疑算法，九成是环境或共享包自身的问题。

原因：共享压缩包常见三个坑。一是解压软件对文件名的编码处理不一，中文文件名在部分工具里解出来直接乱码，脚本里如果硬编码了路径就会立刻报错；二是部分分享者会给 rar 包做伪加密标记，头部标了加密字段但数据区没有真正加密，解压工具会提示输入密码，实际上文件可以正常释放；三是脚本依赖的 pandas、scipy 等库版本和你的环境不一致，换台机器就缺模块。

解决：先把解压目录改成纯英文路径，再换一款解压工具或直接在命令行模式下重试，很多伪加密包在命令行模式下能正常释放；解压后优先检查import段，缺库就在虚拟环境里补装；中文乱码则用文本编辑器把脚本另存为 UTF-8 编码再运行。环境问题排除干净后，如果 G 值计算结果和预期不一致，才轮到去怀疑临界值公式和自由度实现。

5. 把格拉布斯检验沉淀成通用检测函数：一个可直接带进数据流程的模板

日常处理结构化数据时，数据往往不是一维数组，而是带时间戳、批次号的宽表。一个实用的做法是把格拉布斯检验按分组跑，而不是全表一起跑，因为不同批次方差结构不同，混在一起会被最大方差组主导。下面这个模板会按分组给 DataFrame 打上异常标记：

import pandas as pd def flag_outliers_by_group(df, value_col, group_col, alpha=0.05): df = df.copy() df['outlier'] = False for name, group in df.groupby(group_col, sort=False): result = grubbs_test(group[value_col], alpha=alpha) vals = [r['value'] for r in result['removed']] mask = df.index.isin(group.index[group[value_col].isin(vals)]) df.loc[mask, 'outlier'] = True return df

逻辑说明：先按 group_col 分组，每个组单独调用前面写好的grubbs_test，再把结果里被剔除的数值映射回原表并打上outlier标记。分组的意义在于不同批次有各自的均值和离散度，混在一起检验会把高均值批次的正常低值、低均值批次的正常高值都误判成异常。这个模板省去了重复值的精确行号映射，因为格拉布斯判断关心的是数值本身是否离群；如果项目要求逐行可追溯，给源表加一列自增 row_id，再按 row_id 回写即可。

调用方式很直接：

df = pd.read_csv('sensor_log.csv') df = flag_outliers_by_group(df, value_col='temperature', group_col='batch') df[df['outlier']].to_csv('flagged_outliers.csv', index=False)

跑完后我会习惯性检查三件事：剔除比例是否接近 20% 上限、每组 n 是否都大于 8、被标记的数据在原始记录里是否有特殊备注。这三个检查比调整 alpha 本身更能挡掉低级错误。

格拉布斯准则真正的价值不是“替代经验判断”，而是给经验判断一个可查证的统计依据。我把输出里的 G、G_crit、n 一并写进日志，而不是只记录“删除了三条”，这样业务方要求把 0.05 改成 0.01 时，原始统计量还在，不用重跑整条链路。希望这套思路和代码模板能帮你少走几次回头路。

本文还有配套的精品资源，点击获取