简介一套基于格拉布斯准则的异常数据判断代码面向数学建模竞赛和美赛参赛者用于解决数据预处理中的离群点检测问题。该准则通过计算样本最大值与均值的偏离程度并与临界值比较可有效识别正态分布数据中的极端值代码支持显著水平0.01或0.05的选择便于适应不同容错要求。压缩包共3个文件包含txt说明文档、MATLAB主脚本及自动备份文件.asv大小仅1KB逻辑紧凑适合快速调试。目前已有106人学习对新手而言是一份直观的参考实现尤其适合有统计基础的数模爱好者快速理解检验流程。运行后可体验从数据读取、缺失值处理到G值计算、异常值标记的完整流程并可根据需求选择删除或替换策略说明文档还梳理了每一步的计算思路与临界值查找方法方便迁移到其他数据分析项目中长期提升数据质量。1. 格拉布斯准则代码数模预处理里最该先跑的那一步美赛数模里最气人的翻车往往不在模型公式而在数据预处理一组数据里混进一两个离群点均值被拉偏、相关系数虚高、回归系数全变味。格拉布斯准则是解决这类问题的经典工具它用“最远点与均值的距离相当于几个标准差”来判断异常值把“看着不对劲”变成可检验的统计结论。这份基于格拉布斯准则判断异常数据代码.rar就是一套开箱即用的 MATLAB 示例代码解压后直接运行 Untitled2.m改一下数据就能拿到异常值列表和判定依据。适合正在备赛美赛、数模竞赛的队伍也适合做传感器数据清洗的工程人员。代码几十行就能读懂跑通后再按自己的数据改一改比从零手写省一晚上时间。2. 格拉布斯检验的原理统计量怎么构造临界值从哪来用这份代码之前我建议先花十分钟把原理捋清楚。这个方法在网络上流传的公式版本极多有的写法带 (\sqrt{n})有的不带参数稍微用错一点点结果就完全不可信。下面按“统计量怎么构造 → 临界值怎么来 → 什么时候不能用”的顺序过一遍。2.1 统计量的构造逻辑为什么只审问最远的那个点格拉布斯检验的核心逻辑可以压缩成一句话如果一组来自正态总体的样本里最极端的那个点离样本均值的距离远到了“不太可能偶然出现”就把它判为异常值。它构造的统计量是[ G \frac{\max(X) - \bar{X}}{s} ]其中 (\max(X)) 是样本最大值(\bar{X}) 是样本均值(s) 是样本标准差按 (n-1) 计算。这个式子要拆成三层理解。第一层分子“最远点减均值”衡量的是绝对距离第二层除以标准差 (s) 后距离变成了“多少个标准差”也就是标准化第三层标准化后的 (G) 在不同量纲的数据之间可以直接比较所以全世界共用同一套临界值表。和 3σ 准则相比格拉布斯没有把 3 当成固定阈值而是把判定阈值做成了样本量 (n) 和显著性水平 (\alpha) 的函数——小样本时阈值低一点大样本时阈值高一点这比一刀切科学得多。如果怀疑异常出现在最小值端公式对称地写成[ G \frac{\bar{X} - \min(X)}{s} ]实务里我一般两端都算取更大的那个 (G)再和双侧临界值比较。这也是后面代码里的标准写法。2.2 双侧检验与临界值公式alpha 和 n 分别改变了什么很多人拿到格拉布斯公式后最困惑的是临界值 (G_{\alpha}) 怎么算。它不来自常见的 z 表或 t 表而是由 t 分布导出的一个复合形式[ G_{\alpha} \frac{n-1}{\sqrt{n}}\sqrt{\frac{t_{\alpha/(2n),n-2}^{2}}{n-2t_{\alpha/(2n),n-2}^{2}}} ]其中 (t_{\alpha/(2n),n-2}) 是自由度为 (n-2) 的 t 分布上侧分位数。这个公式里有两个需要解释的细节。第一为什么显著性水平要除以 (2n) 而不是直接用 (\alpha)。因为 (G) 检验的是“n 个样本里最极端的那一个”属于多重比较问题(\alpha) 被分摊到 n 个样本和双侧两个端点上这就是 (\alpha/(2n)) 的来源。如果只做单侧检验把 (2n) 改成 (n)。第二为什么自由度是 (n-2)。检验最大值时样本里有一个点被当作候选异常值剩余 (n-1) 个点被用来估计均值再扣掉一个自由度就得到 (n-2)。这个数字在 MATLAB 的 tinv 函数里直接写进去不要顺手改成 (n-1)。先把三种常见方法和格拉布斯放在一起对比选型时更清楚方法判据正态假设一次能剔除的异常值3σ距离 3 倍标准差弱依赖多个IQR 箱线图落在 1.5 倍四分位距之外不依赖多个格拉布斯G 临界值强依赖1 个选型建议数模赛题不会告诉你哪些点是异常IQR 虽然不需要正态假设但写论文时解释力弱3σ 在小样本下阈值固定容易误删格拉布斯胜在“有临界值、有显著性水平、有完整的统计语言”在论文里讲得清楚这是评委能跟上的逻辑。先做正态性检验是使用前提。用 lillietest(data) 或 jbtest(data) 看 p 值p 大于 0.05 再跑格拉布斯如果数据是故障间隔时间、网络延迟这类天然右偏的变量先取对数或做 Box-Cox 变换变换后的分布接近正态再用。2.3 掩蔽效应与适用前提为什么一次只删一个格拉布斯检验在原理上的一个特质是“一轮只审问最远的那个点”。这不是偷懒而是统计上必须如此异常值会把样本均值拉向自己同时把标准差拉大于是其他异常点与均值的标准化距离反而被压缩了看起来像正常点。这就是掩蔽效应。一个很典型的现场是第一轮跑格拉布斯删掉一个点再跑又冒出来两个点继续删又冒出新的。遇到这种情况不用怀疑代码这正是掩蔽效应在起作用。正确做法是每轮只删一个删完重新计算均值和标准差循环直到没有点能超过临界值。这个循环逻辑在第 6 章会封装成函数。使用前提还要补一条样本量。临界值公式里自由度是 (n-2)所以 (n) 至少等于 3 才有意义。而在 (n) 为 38 的小样本区间临界值随 (n) 变化非常陡n3 时是 1.155n4 就直接跳到 1.481此时统计结果对数据极其敏感应当把格拉布斯当“参考”而不是“判决”。网络上很多公式写成 (G (\max(X)-\bar{X})/(s/\sqrt{n}))这是把样本均值当随机变量做标准化的写法不是格拉布斯检验。用这个式子算出的数和本文的临界值表永远对不上这点在第 5 章的避坑部分还会再展开。3. MATLAB 实现把 Untitled2.m 从读数据跑到出结果3.1 解压与文件识别哪个才是主脚本用 WinRAR 或 7-Zip 把压缩包解开后通常会看到这样几个文件文件名类型说明Untitled2.mMATLAB 脚本主代码右键→运行Untitled2.asv自动备份MATLAB 编辑器生成的备份文件可直接忽略实验1文件夹或数据文件可能是原始数据或旧版本脚本www.downma.com.txt文本说明来源说明不参与运行第一次打开时注意双击 .asv 大概率是乱码因为它是编辑器自动保存的中间产物不是给人直接看的。真正的入口只有一个 Untitled2.m。运行前确认环境MATLAB R2016a 及以上即可不需要深度学习工具箱统计工具箱有没有都行——有就用 tinv没有就换成第 4 章的插值表。3.2 核心代码逐段拆解从读数据到标异常我拿到这份代码后一般会先按下面的结构把脚本重写一遍方便替换成自己的数据% 格拉布斯准则判断异常数据 % 用法: 把 data 替换成自己的数据列, 整体运行 data [10.1 10.3 9.8 11.2 10.5 9.9 10.0 29.7 10.2]; % 示例数据 alpha 0.05; % 显著性水平, 可改成0.01 n length(data); % 样本量 mu mean(data); % 样本均值 s std(data); % 样本标准差, MATLAB默认除以(n-1) G_max (max(data) - mu) / s; % 最大值端的G G_min (mu - min(data)) / s; % 最小值端的G G max(G_max, G_min); % 取更极端的一端 % 计算格拉布斯临界值 Gc, 需要统计工具箱 t tinv(1 - alpha/(2*n), n - 2); Gc (n - 1) / sqrt(n) * sqrt(t^2 / (n - 2 t^2)); if G Gc [~, j] max(abs(data - mu)); % 定位离均值最远的点 fprintf(检测到异常值: %.2f (第 %d 个样本), G%.4f Gc%.4f\n, ... data(j), j, G, Gc); else disp(未检测到异常值); end这段代码的逻辑和第 2 章的公式是一一对应的。先说统计量部分G_max 和 G_min 分别检验最大值端和最小值端因为异常点可能出现在任意一端取两者较大值是为了覆盖两端max(abs(data-mu)) 返回的索引 j就是离均值最远的那个点解决了“判出来了却不知道删谁”的问题。再说临界值部分tinv 是 t 分布的分位数函数alpha/(2*n) 是双侧检验的显著性修正自由度 n-2 是格拉布斯统计量的固有属性这两个参数不要动。Gc 的公式就是第 2 章的复合公式拆成两行写是为了可读性。运行这段脚本会在命令窗口输出类似下面的结果检测到异常值: 29.70 (第 8 个样本), G2.65 Gc2.22参数怎么调alpha 取 0.05 是默认选择表示有 5% 的概率把正常点误判为异常如果后续建模对数据完整性要求高可以改成 0.01如果只是想快速排查可疑点放宽到 0.1 也行。data 可以传行向量或列向量脚本里 length 和 mean 都能正确处理如果 data 来自 Excel建议先用 readmatrix(文件名.xlsx) 读进来再赋给 data。n 小于 3 时 tinv 会报错因为自由度变成负数或零。实际使用前先加一句if n 3, error(样本量过少); end。3.3 没有统计工具箱时的替代把临界值表硬编码进脚本很多学校的 MATLAB 许可证不含统计工具箱运行上面的代码会在 tinv 这一行报“未定义函数或变量 tinv”。常见做法是直接把临界值表写进脚本用 interp1 插值代替 tinv% 临界值表(双侧, alpha0.05), 覆盖 n3~50 n_table [3 4 5 6 7 8 9 10 12 15 20 30 50]; g_table [1.155 1.481 1.715 1.887 2.020 2.126 2.215 2.290 ... 2.412 2.549 2.709 2.908 3.128]; if n max(n_table) Gc interp1(n_table, g_table, n, linear, extrap); else warning(n 超过50, 插值结果仅供参考); endinterp1 的 linear 方式做线性插值n 不在表里时也能拿到一个近似临界值extrap 允许外推但 n 超过 50 后临界值变化趋缓线性外推误差会变大只建议把它当应急方案。如果队伍里能装 Python建议还是用第 6 章里的 scipy 版本做精确计算。我先跑完这份代码再去换自己的数据一步步核对输出。这是整个流程里最需要耐心的一段数据列换错、单位没统一、缺失值直接空着传进去都会让结果看起来像“玄学”。把这些细节先排除再谈结果对不对。4. 临界值的三种取法查表、t 分布逆函数与在线工具很多人跑格拉布斯检验时都卡在同一个地方临界值 (G_{\alpha}) 从哪来这一章把三种常用取法都过一遍并给出比赛环境下的选型建议。4.1 常用临界值速查表alpha0.05 与 0.01下面是从标准格拉布斯临界值表里摘出的常用行对应双侧检验数值四舍五入到三位小数nalpha0.05alpha0.0131.1551.15541.4811.49651.7151.76461.8871.97372.0202.13982.1262.27492.2152.387102.2902.482122.4122.636152.5492.806202.7093.001302.9083.236503.1283.483这张表有两个容易忽略的规律。第一同一行里 alpha0.01 的临界值比 0.05 大说明想降低“误杀正常点”的概率判定门槛就要抬高。第二临界值随 n 增大而增大样本量越大极值本身的随机波动越小判定阈值反而更严格。这和小样本直觉刚好相反写论文做解释时注意别写反。查表法的局限在于 n 必须落在表内。n27 这种常用尺寸表里没有线性插值也只能给近似值这在正式比赛里不够严谨。所以我的建议是查表只用来手算验算脚本里用 4.2 的方法。4.2 用 t 分布逆函数直接算一行代码覆盖任意 n精确临界值不需要翻书用 t 分布逆函数几行就能算t tinv(1 - alpha/(2*n), n - 2); Gc (n - 1) / sqrt(n) * sqrt(t^2 / (n - 2 t^2));这是完整公式第一行算出修正后的 t 分位数第二行代入格拉布斯临界值表达式。alpha0.05、n9 时上式代进去 Gc 约为 2.216和速查表的 2.215 基本一致。Python 环境里对应写法是scipy.stats.t.ppf(1 - alpha/(2*n), n-2)自由度同样是 n-2。这条路径不依赖任何统计表n 取任意大于 2 的整数都能算是脚本内联的首选。4.3 三种方式对比与选型建议方式精度依赖适合场景查表受表范围限制无手工验算、课堂演示、快速确认tinv/ppf 公式精确统计工具箱或 SciPy比赛脚本内联推荐在线计算器取决于实现联网临时核对单个 n 值不推荐在比赛里依赖我的习惯是MATLAB 装统计工具箱就用 tinv没装就把 4.1 的表硬编码进脚本interp1 插值应急Python 环境直接用 scipy。绝不为了省事在论文里写“n27 时临界值是……”然后拿表里没有的数字硬凑评委一眼就能看出来。另外补一个容易绕晕的细节如果只检验最大值端临界值要按单侧查同时检验两端才用本文这张双侧表。代码里同时算了 G_max 和 G_min对应双侧临界值两者口径必须保持一致。5. 避坑手册格拉布斯检验在数模里的五个高频翻车点把代码跑通只是第一步真正把它用到自己的数据上坑才会一个个浮出来。以下五条是我在给数据做清洗时踩过或见过的典型问题按“现象→原因→解决”整理。5.1 现象G 明显很大却和临界值表对不上现象自己按网上公式手算G 超过临界值好几倍坚信一定有异常换用脚本一跑结果完全相反。原因网上流传的公式经常把统计量写成 (G (\max(X)-\bar{X})/(s/\sqrt{n}))。这是把“样本均值”当随机变量做标准化的写法本质是 z 检验的形态不是格拉布斯检验。格拉布斯统计量就是 ((\max(X)-\bar{X})/s)没有 (\sqrt{n})。用错公式再去查格拉布斯临界值表相当于拿两把不同的尺子去量同一个物体。解决脚本里删掉 sqrt(n) 相关部分只保留 ((\max(X)-\bar{X})/s)同时确认自己用的是“单侧”还是“双侧”临界值本文代码里是双侧。5.2 现象删掉一个异常值又冒出两个新异常现象第一轮标出 1 个异常剔除后重跑又标出 2 个再删再标好像永远删不完。原因掩蔽效应。异常值会把样本均值拉向自己、把标准差拉大其他离群点的标准化距离因此被压缩看起来像正常点每删掉一个异常均值和标准差就更接近真实分布原来被掩护的点就暴露了。解决接受这个行为采用“每轮只删一个、删完重新计算”的迭代策略。手动一次一次点 Run 容易出错建议直接用第 6 章的封装函数它会自动循环到没有异常为止。5.3 现象解压后打开 .asv 文件内容乱码还报错现象压缩包里有两个 Untitled 相关文件双击 .asv 打开后满屏乱码运行直接报错。原因.asv 是 MATLAB 编辑器自动保存的备份文件保存的是编辑过程中的中间状态格式与 .m 不同内容也可能不完整。它不是给人直接运行的源文件。解决主代码文件是 Untitled2.m.asv 可以直接删除。建议在 MATLAB 的“主页→预设→编辑器/调试器→自动保存”里关掉自动保存或把时间间隔调长避免以后解压出更多 .asv 文件造成混淆。5.4 现象数据明显偏态格拉布斯把正常大值判成异常现象数据是故障间隔时间、网络延迟这类右偏分布明明都是正常样本格拉布斯却连续标出七八个异常值。原因格拉布斯基于正态分布假设。偏态分布的尾巴天然比正态分布厚极值出现的频率远高于正态模型拿正态临界值去判当然处处“显著”。解决先跑 lillietest(data) 或 jbtest(data)p 值小于 0.05 时先取 log 或做 Box-Cox 变换变换后分布接近正态再跑格拉布斯如果变换后仍不满足正态就改用箱线图 IQR 方法或 MAD 中位数绝对偏差论文里说明替换理由即可不是非用某一种方法不可。5.5 现象n 只有 46 个样本结论摇摆不定现象同一组小样本数据alpha 取 0.05 判为异常取 0.01 就不算换掉其中一两个点结论直接反转。原因n 在 38 之间时临界值随 n 变化非常陡峭n3 是 1.155n4 就跳到 1.481统计功效本来就弱加上小样本下均值和标准差对单点极其敏感结果自然不稳定。解决n 小于 8 时把格拉布斯当参考而不是判决。可以结合残差图或业务上下文人工确认如果实在无法确认先保留数据在建模阶段用 Huber 回归这类稳健方法让模型自己对离群点降权。这五条如果全踩一遍一晚上基本就没了。我把它们挂在手边当检查清单每次跑格拉布斯前逐条过一遍能省下大量复盘时间。6. 进阶用法迭代剔除与 Python 交叉验证单次检验在真实清洗里不够用掩蔽效应要求每轮只删一个、删完重算。手动循环容易删错直接封装成函数最省心function [x_clean, r] grubbs_clean(x, alpha) if nargin 2, alpha 0.05; end x_clean x(:); r []; while numel(x_clean) 3 mu mean(x_clean); s std(x_clean); if s 0, break; end [d, j] max(abs(x_clean - mu)); n numel(x_clean); t tinv(1 - alpha/(2*n), n - 2); Gc (n - 1)/sqrt(n) * sqrt(t^2/(n - 2 t^2)); if d/s Gc r [r; x_clean(j)]; x_clean(j) []; else break; end end end每轮循环重新计算 n 和临界值只删离均值最远的一个点删完回到 while 重新判断直到没有点超过临界值。s 等于 0 说明所有样本相同直接退出。美赛常见配置是 MATLAB 建模、Python 画图用下面这段做交叉验证两边结果一致再进论文import numpy as np from scipy import stats def grubbs_clean(x, alpha0.05): x np.asarray(x, dtypefloat).ravel() removed [] while len(x) 3: mu, s x.mean(), x.std(ddof1) if s 0: break j np.argmax(np.abs(x - mu)) t stats.t.ppf(1 - alpha/(2*len(x)), len(x)-2) Gc (len(x)-1)/np.sqrt(len(x)) * np.sqrt(t*t/(len(x)-2t*t)) if np.abs(x[j]-mu)/s Gc: removed.append(x[j]) x np.delete(x, j) else: break return x, np.array(removed)np.std(ddof1) 对应 MATLAB 的 std()两套公式完全一致删掉的点应该逐项对得上。从那以后我拿到赛题数据都会先把原始数据复制一份再在副本上跑格拉布斯迭代同时开一张箱线图对照。这套流程让我再也没做过“删错列、半天找不回原数据”的后悔事。希望帮到你。本文还有配套的精品资源点击获取