简介面向数学建模与美赛的数据预处理需求压缩包内代码基于格拉布斯准则实现异常值判断用于识别和修正样本中的极端数据适合参赛选手或数据分析初学者参考。格拉布斯检验以正态分布为前提计算最大值与均值的偏差并转换为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 做分母也就是 ddof1 的样本标准差如果用总体标准差 ddof0小样本时 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样本量 nα 0.05双侧α 0.01双侧31.1551.15551.7151.749102.2902.410152.5492.705202.7092.884302.9083.145503.1283.483观察这张表能发现两个特点样本量越小临界值越接近 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(ddof1) # 样本标准差分母 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(fG统计量: {G:.4f})这段代码的运行结果会很直观16.8 这个读数离均值最远G 值大约在 1.9 左右。注意ddof1是必须写的numpy.std()默认 ddof0 算总体标准差偏差会被低估把 G 值虚高。先运行这一段的目的不是判断异常而是让你看到每个中间量尤其是标准差对结果的敏感程度把样本数据任意一个值改大标准差变大G 值未必同步变大这就是格拉布斯在多个极端点同时存在时容易失灵的根源。有pandas环境时可以把data换成df[col].to_numpy()后续判断直接输出一个标记列。常见做法是不把这段逻辑散落在报告里而是封装成一个小函数便于在不同脚本间复用。3.2 完整实现查临界值、单次判断与循环剔除接下来给出一版可以直接放进项目里用的完整函数。它同时处理临界值计算、单次判断和迭代剔除三个动作import numpy as np from scipy import stats def grubbs_test(data, alpha0.05, sidetwo, max_remove0.2): 基于格拉布斯准则判断并剔除异常数据。 参数 data: 一维数组-like 数据至少 3 个观测点 alpha: 显著性水平常用 0.05 side: two 双侧检验high 只查上端low 只查下端 max_remove: 最多剔除比例防止迭代时误删正常数据 arr np.asarray(data, dtypefloat) 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(ddof1) 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, alpha0.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 更小判定更严格。选择的关键不是数据长什么样而是业务上是否接受“两个方向都算异常”。解决开工前把业务规则写进脚本注释例如“此工序只允许正向偏差使用 sidehigh”。拿不准时维持双侧并在输出里同时保留原始值和被剔除值方便回查。把方向决策留给最懂业务的人不要由写代码的人顺手拍板。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, alpha0.05): df df.copy() df[outlier] False for name, group in df.groupby(group_col, sortFalse): result grubbs_test(group[value_col], alphaalpha) 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_coltemperature, group_colbatch) df[df[outlier]].to_csv(flagged_outliers.csv, indexFalse)跑完后我会习惯性检查三件事剔除比例是否接近 20% 上限、每组 n 是否都大于 8、被标记的数据在原始记录里是否有特殊备注。这三个检查比调整 alpha 本身更能挡掉低级错误。格拉布斯准则真正的价值不是“替代经验判断”而是给经验判断一个可查证的统计依据。我把输出里的 G、G_crit、n 一并写进日志而不是只记录“删除了三条”这样业务方要求把 0.05 改成 0.01 时原始统计量还在不用重跑整条链路。希望这套思路和代码模板能帮你少走几次回头路。本文还有配套的精品资源点击获取