简介面向MATLAB用户与数据分析学习者这份资源聚焦蒙特卡洛模拟在异常值剔除中的具体应用解决传统统计方法难以处理复杂数据分布时识别偏离样本的难题。压缩包共19个文件包含15个m脚本、2个中英文讲义ppt、1个mat格式数据文件及1个说明文档整体仅388KB脚本涵盖蒙特卡洛模拟演示、方差缩减、湖面面积、投资组合等经典案例并配有随机数生成与统计检验代码便于直接运行与二次修改。目前已有1080人学习使用。通过对照英文与法文PPT框架配合readme说明与可执行脚本读者可系统理解基于随机抽样构建分布模型、识别并剔除异常值的完整流程同时掌握boxplot可视化、normfit参数估计及Z-score、IQR等常用检验方法的MATLAB实现适合需要快速上手蒙特卡洛方法或改进异常值处理流程的初、中级研究者。1. 蒙特卡洛异常值剔除为什么随机抽样能定位离群点一份带异常值的数据集摆到面前时直接算均值、拟合回归结果往往被少数离群点牵着走。传统手段用 3σ 规则或箱线图也能筛但它们默认数据服从正态分布且阈值完全由样本自身决定——当样本里已经混入异常值时均值和标准差本身就被污染了筛完可能漏掉真异常也可能误杀正常点。蒙特卡洛思路的不同之处在于它先用正常部分的分布假设生成大量随机样本构建“如果数据真是干净的话应该长什么样”的经验包络再拿原始点去对包络偏离超过分位数的才判为异常。这相当于给数据清洗加了一个可复现的参照系。下面用 MATLAB 把这套流程完整走一遍覆盖参数估计、模拟样本生成、迭代剔除与验证四个环节落地到可以直接跑的脚本。2. 蒙特卡洛剔除异常值的统计原理与分布假设2.1 异常值的统计定义与蒙特卡洛的定位先说清楚什么是异常值。异常值是明显偏离主体观测的数据点来源可能是测量误差、录入错误也可能是真实存在的极端事件。统计意义上异常值在给定分布模型下的出现概率极低。这里的难点是“概率低”是相对哪个分布而言——分布错了异常判定就错了。蒙特卡洛在异常值剔除中扮演的角色是用随机抽样近似一个参照分布再把观测值和该分布的分位数比较。换言之蒙特卡洛不负责“找异常”它负责生成“如果没有异常数据应该长什么样”的对照。这是它和确定性阈值方法最大的区别确定性方法用样本自身估计参数再切阈值而样本被污染时均值、标准差都偏离真实值这叫掩蔽效应。蒙特卡洛抽样并不能免疫掩蔽效应但它给了我们一个反复迭代的框架每剔除一轮就重新抽样、重新估计逐步逼近干净分布。我拿到 MonteCarlo.rar 这类资源时一般会先看目录结构MyMC 目录通常放自定义的蒙特卡洛核心函数VarReduction 是方差缩减实现PortSim 偏投资组合模拟LakeArea 常用于面积型随机模拟readme.txt 记录运行环境和依赖。异常值剔除不是这些模块的主线但蒙特卡洛模拟框架完全可以复用来构建异常判定的参照分布。2.2 三条可落地的异常判定路径实际用蒙特卡洛做异常值剔除常见三条路径判定逻辑都建立在模拟样本的分位数上而不是直接查统计表。2.2.1 基于模拟分位数的 Z-score 路径传统 Z-score 判定计算 z (x - μ) / σ看 |z| 是否超过 3。当 μ 和 σ 用被污染的样本估计时掩蔽效应会让真实异常点的 z 值变小。蒙特卡洛思路是先假设数据来自某个分布用模拟样本估计 z 的经验分布取 0.025 和 0.975 分位数作为正常范围。这样即使原分布不正态也能得到对应的分位阈值。这条路径的优点是直观、容易可视化缺点是对初始分布假设敏感。2.2.2 IQR 路径与模拟结合箱线图里的 IQR 法则不需要正态假设但 1.5 倍 IQR 这个经验阈值在重尾分布下会把大量正常点标成异常。可以在蒙特卡洛框架内对模拟样本计算 IQR观察 IQR 的经验分布再用高分位数放大或收缩判定带宽。也就是说把 IQR 的倍数从固定值变成由模拟分布决定的随机量。2.2.3 Grubbs 检验路径Grubbs 检验专门针对单变量正态数据逐点检测最大偏离点统计量 G max|x_i - μ| / σ。先用模拟样本生成 G 统计量的经验分布再把原始数据里最偏离点拿进去比 p 值。相比查表模拟版能适应非标准样本量得到更贴近实际的临界值。较新的 Statistics Toolbox 自带 grubbs 函数旧版本需要手写 G 统计量并用 tinv 查临界值。三条路径的参数对应关系如下判定路径分布假设需要设置的参数MATLAB 入手函数Z-score 分位数任意可模拟分位水平 αnormfit、quantileIQR 经验带宽任意可模拟IQR 倍数或分位数iqr、quantileGrubbs 经验临界值近似正态显著性水平max、tinv、grubbs我一般先用 Z-score 分位数做初筛再用 Grubbs 交叉验证。两个方法同时标记的点剔除置信度明显更高。2.3 适用边界小样本与分布误设蒙特卡洛路径不是万能的。样本量小于 30 时MLE 估计出的分布参数本身波动大模拟分位数的方差也大容易把正常点误判成异常。分布假设错误是另一个更隐蔽的坑真实数据是重尾分布却硬套正态模型蒙特卡洛包络会偏窄尾部正常点全被标红。我通常先用分布拟合优度检验对比正态、对数正态、t location-scale 三种候选分布再看模拟包络的宽度变化决定到底用哪个分布做抽样。这一步比单纯增大模拟样本量更关键。3. MATLAB 脚本实战数据预处理与模拟样本生成3.1 数据读取与探索性分析拿到数据第一步不是急着写蒙特卡洛循环而是把分布形态看清楚。我一般先读入数据计算基本统计量画直方图和箱线图重点看偏度、峰度以及箱线图外的点。数据文件格式上CSV 或 Excel 都行readmatrix 对两类文件都支持如果源数据是 Excel 里多个 sheet用 readmatrix 的 Sheet 参数指定避免读错表。% 读取单列数值型数据假定第一列为观测值 data readmatrix(obs_data.csv); data data(:, 1); % 只保留第一列 n length(data); % 探索性统计量 mu0 mean(data); sig0 std(data); fprintf(原始样本: n%d, mean%.3f, std%.3f\n, n, mu0, sig0); % 直方图与箱线图 figure(Color, w); subplot(1, 2, 1); histogram(data, 30, FaceColor, [0.3 0.6 0.9]); title(原始数据直方图); xlabel(观测值); ylabel(频数); subplot(1, 2, 2); boxplot(data, Symbol, r); title(原始数据箱线图); ylabel(观测值);逻辑是先读数据并固定到单列算均值和标准差作为后续对比基准直方图看分布形态箱线图快速暴露离群点。readmatrix 在较新版本 MATLAB 中可用旧版本换成 csvread 或 importdata。boxplot 的 Symbol 参数控制离群点标记默认红色加号能直观看到哪些点落在箱线图外。注意这里的初始统计量是被污染过的不能直接当作剔除阈值的依据只用于对照“如果不处理统计量会偏到哪里”。真正用于判定的参数要在干净样本上迭代估计。3.2 正态假设下的参数估计与蒙特卡洛样本生成蒙特卡洛样本生成的核心是先用极大似然估计得到分布参数再用随机数生成函数构造与原数据量级一致的模拟样本。MATLAB 的标准做法如下% 极大似然估计正态分布参数 [mu_hat, sigma_hat] normfit(data); fprintf(MLE: mu%.3f, sigma%.3f\n, mu_hat, sigma_hat); % 固定随机种子保证结果可复现 rng(2024); % 构造分布对象并生成蒙特卡洛样本 pd makedist(Normal, mu, mu_hat, sigma, sigma_hat); n_mc 10000; mc random(pd, n_mc, 1); % 对模拟样本标准化用于后续分位数对照 z_mc (mc - mu_hat) / sigma_hat;normfit 返回均值和标准差的最大似然估计它比直接 meanstd 受到样本污染的收缩稍稳。makedist 将参数封装成分布对象random 从该分布中抽出 10000 个样本。rng(2024) 固定随机种子保证同一份数据每次运行结果一致这个习惯在调试和出报告时非常有用。对模拟样本做标准化是为了与后续原始数据的 Z-score 对齐尺度。如果数据明显非正态把 makedist 的 Normal 换成 tLocationScale 或 Lognormal 即可其余代码不用动。用分布对象而不是直接调 random(Normal, mu, sigma, n_mc, 1) 的好处也在这里切换分布时只需改 makedist 一行参数估计和标准化代码完全复用。3.3 模拟样本量怎么定模拟样本量 n_mc 决定了分位数估计的稳定性。样本太少尾部分位数每次运行都在变样本太多计算成本线性上升收益却很小。我常用的参考量级如下n_mc0.5% 分位数相对波动适用场景1000±15% 左右快速探索不追求稳定10000±4% ~ ±6%单变量数据通用推荐100000±1.5% ~ ±2%多维或高显著水平这里的相对波动是不同随机种子下同一样本量分位数变化的经验范围。对单变量异常值剔除10000 个模拟样本足够如果原始数据本身只有几百个点模拟样本再多也无法提升参数估计精度瓶颈在 MLE 而不在抽样。如果数据包含明确的上下物理边界比如传感器量程我会把模拟样本直接截断到量程内再计算分位数这样包络不会延伸到物理上不可能的区域。4. 异常点识别、迭代剔除与收敛判据4.1 用蒙特卡洛包络识别潜在异常点有了模拟样本的经验分布下一步把原始观测值与模拟包络比较。这里采用双尾 1% 分位水平即 0.5% 和 99.5% 分位数比常规 3σ 更严格适合剔除阶段使用。alpha 0.01; lo quantile(z_mc, alpha / 2); hi quantile(z_mc, 1 - alpha / 2); % 原始数据标准化 z_obs (data - mu_hat) / sigma_hat; % 超出包络即为潜在异常 flag (z_obs lo) | (z_obs hi); idx_out find(flag); fprintf(首次识别潜在异常点 %d 个\n, length(idx_out)); % 包络可视化 figure(Color, w); hold on; scatter(1:n, data, 20, [0.5 0.5 0.5], filled); plot(1:n, mu_hat lo * sigma_hat, r--, LineWidth, 1.5); plot(1:n, mu_hat hi * sigma_hat, r--, LineWidth, 1.5); xlabel(样本序号); ylabel(观测值); legend({原始数据, MC 下界, MC 上界}, Location, best); hold off;核心逻辑是计算模拟样本标准化值的分位数作为包络边界再代入原始观测值超出边界的标记为潜在异常。红色虚线是包络不是回归线表示“在给定分布假设下正常观测值应该落在这个范围内”。alpha 直接影响识别结果。0.01 偏保守适合样本量大、异常占比低的数据0.05 更激进适合异常占比高、需要强清洗的场景。quantile 接受 0 到 1 的比例prctile 接受 0 到 100这两个函数经常被混用我习惯统一用 quantile改分位水平时不用换算。4.2 迭代剔除每剔一点就重新估计参数一次性剔除所有超包络点是危险动作异常值参与均值和标准差估计包络本身被拉宽部分真异常可能藏在包络内部。常见做法是逐轮剔除每轮只用当前干净数据重新做 MLE 和蒙特卡洛模拟再剔除超出新包络的点直到没有新异常或达到迭代上限。data_clean data; removed_idx []; max_iter 10; n_mc 20000; for iter 1:max_iter n_now length(data_clean); if n_now 10 warning(剩余样本过少停止迭代); break; end % 当前干净集上的分布估计 [mu_i, sigma_i] normfit(data_clean); if sigma_i 1e-8 break; end % 生成蒙特卡洛样本并计算双尾阈值 mc_i random(... makedist(Normal, mu, mu_i, sigma, sigma_i), n_mc, 1); z_i (mc_i - mu_i) / sigma_i; lo_i quantile(z_i, alpha / 2); hi_i quantile(z_i, 1 - alpha / 2); % 检查当前数据 z_cur (data_clean - mu_i) / sigma_i; bad find(z_cur lo_i | z_cur hi_i); if isempty(bad) fprintf(第 %d 轮收敛无新异常\n, iter); break; end removed_idx [removed_idx; bad]; data_clean(bad) []; fprintf(第 %d 轮: 剔除 %d 个点剩余 %d 个\n, ... iter, length(bad), length(data_clean)); end循环结构是估计参数 → 生成模拟样本 → 计算包络 → 标记越界点 → 有点则剔除并进入下一轮。注意每轮都要重新拟合不能沿用最初的 mu_hat 和 sigma_hat否则掩蔽效应会一直存在。注意迭代循环里不要复用第一轮估计的 mu_hat 和 sigma_hat 作为每轮判定参数掩蔽效应会随着异常值残留累积导致永远收不了敛。停止条件有三个没有新异常被标记、剩余样本少于 10 个、达到最大迭代次数。前两个是数据层面的收敛第三个是保护性上限防止异常占比极高时循环把样本清空。迭代轮数超过 5 轮时我会打一条日志看每轮剔除量是否递减如果第二轮剔除数量比第一轮还多说明分布假设可能有问题先停手排查分布。4.3 剔除前后统计量与模型稳定性对比剔除完成后需要验证。验证对象不是“剔得多干净”而是“统计量是否稳定下来”。我对比剔除前后的均值、标准差、偏度和峰度并重新画箱线图。mu_after mean(data_clean); sig_after std(data_clean); fprintf(剔除前: n%d, mean%.3f, std%.3f, skew%.2f, kurt%.2f\n, ... n, mu0, sig0, skewness(data), kurtosis(data)); fprintf(剔除后: n%d, mean%.3f, std%.3f, skew%.2f, kurt%.2f\n, ... length(data_clean), mu_after, sig_after, ... skewness(data_clean), kurtosis(data_clean)); figure(Color, w); subplot(1, 2, 1); boxplot(data, Symbol, r); title(剔除前); subplot(1, 2, 2); boxplot(data_clean, Symbol, r); title(剔除后); % 交叉验证Grubbs 检验新版 Statistics Toolbox 提供 if exist(grubbs, file) [h_grubbs, idx_grubbs] grubbs(data, 0.05); fprintf(Grubbs 标记异常点 %d 个\n, sum(h_grubbs)); end判定要落到量级上均值变化不大但标准差显著缩小说明异常值主要拉宽了分布没有造成系统性偏移均值也明显变化时说明异常值的分布不对称剔除后的分析结论需要重写。比如一组模拟数据的对比指标剔除前剔除后变化方向样本量500486减少 2.8%均值12.7412.31下降 0.43标准差3.912.68下降 31.5%偏度1.820.34趋近对称峰度8.753.21趋近正态标准差下降 30% 以上是常见信号说明原本数据的变异主要来自离群点而不是主体分布。grubbs 函数在较新的 Statistics Toolbox 中提供旧版本可以手写 G 统计量并用 tinv 查临界值。蒙特卡洛和 Grubbs 同时标记的点剔除置信度更高只被一方标记的点保留下来进 review 名单。5. 蒙特卡洛剔除异常值的阈值调优与业务收口蒙特卡洛剔除异常值的参数里alpha 是最敏感的一个。我通常先跑一轮 alpha0.05 的初筛统计潜在异常占比若占比超过 10%说明数据质量比预想差不能单靠蒙特卡洛要回去查采集链路若占比在 2% 到 5% 之间再改用 alpha0.01 跑正式迭代。这个两段式策略比直接跑一个固定阈值稳定得多。分布假设不成立时优先换 t location-scale 分布。makedist 支持 tLocationScale拟合函数是 fitdist(data, tLocationScale)。相比正态分布它的尾部更厚模拟包络更宽误杀率会明显下降。另一个思路是做分布无关的 bootstrap对干净样本做有放回重采样每次计算均值或其他统计量用重采样分布替代解析分布。这类方法在 LakeArea 这种面积模拟、PortSim 这类投资组合模拟里也常用——先剔除异常观测再做蒙特卡洛预测预测区间的稳定性会好很多。迭代剔除在小样本下存在一个隐蔽风险每轮剔除后样本量变小MLE 方差变大包络可能收缩下一轮误杀更多点。我一般会加一个保护条件单轮剔除比例超过当前样本量的 5% 时停止迭代并输出警告让业务方确认这批数据是否真的异常占比这么高。数学上这是合理的保守条件工程上是防呆设计。最后给一个落地技巧不要把剔除结果直接覆盖原数据而是在原表上加一列 mc_label取值 keep、remove 或 review。蒙特卡洛包络外且与 Grubbs 交叉验证通过的点标为 remove包络边界附近z 值落在 0.9% 到 1% 分位区间标为 review其余标为 keep。这个三分类方案比一味删除更适合真实数据分析流程——保留原始记录、标注去向、留给下游业务判断才是蒙特卡洛异常值剔除在生产环境里的正确形态。本文还有配套的精品资源点击获取
蒙特卡洛异常值剔除：基于MATLAB的统计模拟与数据清洗实践
简介面向MATLAB用户与数据分析学习者这份资源聚焦蒙特卡洛模拟在异常值剔除中的具体应用解决传统统计方法难以处理复杂数据分布时识别偏离样本的难题。压缩包共19个文件包含15个m脚本、2个中英文讲义ppt、1个mat格式数据文件及1个说明文档整体仅388KB脚本涵盖蒙特卡洛模拟演示、方差缩减、湖面面积、投资组合等经典案例并配有随机数生成与统计检验代码便于直接运行与二次修改。目前已有1080人学习使用。通过对照英文与法文PPT框架配合readme说明与可执行脚本读者可系统理解基于随机抽样构建分布模型、识别并剔除异常值的完整流程同时掌握boxplot可视化、normfit参数估计及Z-score、IQR等常用检验方法的MATLAB实现适合需要快速上手蒙特卡洛方法或改进异常值处理流程的初、中级研究者。1. 蒙特卡洛异常值剔除为什么随机抽样能定位离群点一份带异常值的数据集摆到面前时直接算均值、拟合回归结果往往被少数离群点牵着走。传统手段用 3σ 规则或箱线图也能筛但它们默认数据服从正态分布且阈值完全由样本自身决定——当样本里已经混入异常值时均值和标准差本身就被污染了筛完可能漏掉真异常也可能误杀正常点。蒙特卡洛思路的不同之处在于它先用正常部分的分布假设生成大量随机样本构建“如果数据真是干净的话应该长什么样”的经验包络再拿原始点去对包络偏离超过分位数的才判为异常。这相当于给数据清洗加了一个可复现的参照系。下面用 MATLAB 把这套流程完整走一遍覆盖参数估计、模拟样本生成、迭代剔除与验证四个环节落地到可以直接跑的脚本。2. 蒙特卡洛剔除异常值的统计原理与分布假设2.1 异常值的统计定义与蒙特卡洛的定位先说清楚什么是异常值。异常值是明显偏离主体观测的数据点来源可能是测量误差、录入错误也可能是真实存在的极端事件。统计意义上异常值在给定分布模型下的出现概率极低。这里的难点是“概率低”是相对哪个分布而言——分布错了异常判定就错了。蒙特卡洛在异常值剔除中扮演的角色是用随机抽样近似一个参照分布再把观测值和该分布的分位数比较。换言之蒙特卡洛不负责“找异常”它负责生成“如果没有异常数据应该长什么样”的对照。这是它和确定性阈值方法最大的区别确定性方法用样本自身估计参数再切阈值而样本被污染时均值、标准差都偏离真实值这叫掩蔽效应。蒙特卡洛抽样并不能免疫掩蔽效应但它给了我们一个反复迭代的框架每剔除一轮就重新抽样、重新估计逐步逼近干净分布。我拿到 MonteCarlo.rar 这类资源时一般会先看目录结构MyMC 目录通常放自定义的蒙特卡洛核心函数VarReduction 是方差缩减实现PortSim 偏投资组合模拟LakeArea 常用于面积型随机模拟readme.txt 记录运行环境和依赖。异常值剔除不是这些模块的主线但蒙特卡洛模拟框架完全可以复用来构建异常判定的参照分布。2.2 三条可落地的异常判定路径实际用蒙特卡洛做异常值剔除常见三条路径判定逻辑都建立在模拟样本的分位数上而不是直接查统计表。2.2.1 基于模拟分位数的 Z-score 路径传统 Z-score 判定计算 z (x - μ) / σ看 |z| 是否超过 3。当 μ 和 σ 用被污染的样本估计时掩蔽效应会让真实异常点的 z 值变小。蒙特卡洛思路是先假设数据来自某个分布用模拟样本估计 z 的经验分布取 0.025 和 0.975 分位数作为正常范围。这样即使原分布不正态也能得到对应的分位阈值。这条路径的优点是直观、容易可视化缺点是对初始分布假设敏感。2.2.2 IQR 路径与模拟结合箱线图里的 IQR 法则不需要正态假设但 1.5 倍 IQR 这个经验阈值在重尾分布下会把大量正常点标成异常。可以在蒙特卡洛框架内对模拟样本计算 IQR观察 IQR 的经验分布再用高分位数放大或收缩判定带宽。也就是说把 IQR 的倍数从固定值变成由模拟分布决定的随机量。2.2.3 Grubbs 检验路径Grubbs 检验专门针对单变量正态数据逐点检测最大偏离点统计量 G max|x_i - μ| / σ。先用模拟样本生成 G 统计量的经验分布再把原始数据里最偏离点拿进去比 p 值。相比查表模拟版能适应非标准样本量得到更贴近实际的临界值。较新的 Statistics Toolbox 自带 grubbs 函数旧版本需要手写 G 统计量并用 tinv 查临界值。三条路径的参数对应关系如下判定路径分布假设需要设置的参数MATLAB 入手函数Z-score 分位数任意可模拟分位水平 αnormfit、quantileIQR 经验带宽任意可模拟IQR 倍数或分位数iqr、quantileGrubbs 经验临界值近似正态显著性水平max、tinv、grubbs我一般先用 Z-score 分位数做初筛再用 Grubbs 交叉验证。两个方法同时标记的点剔除置信度明显更高。2.3 适用边界小样本与分布误设蒙特卡洛路径不是万能的。样本量小于 30 时MLE 估计出的分布参数本身波动大模拟分位数的方差也大容易把正常点误判成异常。分布假设错误是另一个更隐蔽的坑真实数据是重尾分布却硬套正态模型蒙特卡洛包络会偏窄尾部正常点全被标红。我通常先用分布拟合优度检验对比正态、对数正态、t location-scale 三种候选分布再看模拟包络的宽度变化决定到底用哪个分布做抽样。这一步比单纯增大模拟样本量更关键。3. MATLAB 脚本实战数据预处理与模拟样本生成3.1 数据读取与探索性分析拿到数据第一步不是急着写蒙特卡洛循环而是把分布形态看清楚。我一般先读入数据计算基本统计量画直方图和箱线图重点看偏度、峰度以及箱线图外的点。数据文件格式上CSV 或 Excel 都行readmatrix 对两类文件都支持如果源数据是 Excel 里多个 sheet用 readmatrix 的 Sheet 参数指定避免读错表。% 读取单列数值型数据假定第一列为观测值 data readmatrix(obs_data.csv); data data(:, 1); % 只保留第一列 n length(data); % 探索性统计量 mu0 mean(data); sig0 std(data); fprintf(原始样本: n%d, mean%.3f, std%.3f\n, n, mu0, sig0); % 直方图与箱线图 figure(Color, w); subplot(1, 2, 1); histogram(data, 30, FaceColor, [0.3 0.6 0.9]); title(原始数据直方图); xlabel(观测值); ylabel(频数); subplot(1, 2, 2); boxplot(data, Symbol, r); title(原始数据箱线图); ylabel(观测值);逻辑是先读数据并固定到单列算均值和标准差作为后续对比基准直方图看分布形态箱线图快速暴露离群点。readmatrix 在较新版本 MATLAB 中可用旧版本换成 csvread 或 importdata。boxplot 的 Symbol 参数控制离群点标记默认红色加号能直观看到哪些点落在箱线图外。注意这里的初始统计量是被污染过的不能直接当作剔除阈值的依据只用于对照“如果不处理统计量会偏到哪里”。真正用于判定的参数要在干净样本上迭代估计。3.2 正态假设下的参数估计与蒙特卡洛样本生成蒙特卡洛样本生成的核心是先用极大似然估计得到分布参数再用随机数生成函数构造与原数据量级一致的模拟样本。MATLAB 的标准做法如下% 极大似然估计正态分布参数 [mu_hat, sigma_hat] normfit(data); fprintf(MLE: mu%.3f, sigma%.3f\n, mu_hat, sigma_hat); % 固定随机种子保证结果可复现 rng(2024); % 构造分布对象并生成蒙特卡洛样本 pd makedist(Normal, mu, mu_hat, sigma, sigma_hat); n_mc 10000; mc random(pd, n_mc, 1); % 对模拟样本标准化用于后续分位数对照 z_mc (mc - mu_hat) / sigma_hat;normfit 返回均值和标准差的最大似然估计它比直接 meanstd 受到样本污染的收缩稍稳。makedist 将参数封装成分布对象random 从该分布中抽出 10000 个样本。rng(2024) 固定随机种子保证同一份数据每次运行结果一致这个习惯在调试和出报告时非常有用。对模拟样本做标准化是为了与后续原始数据的 Z-score 对齐尺度。如果数据明显非正态把 makedist 的 Normal 换成 tLocationScale 或 Lognormal 即可其余代码不用动。用分布对象而不是直接调 random(Normal, mu, sigma, n_mc, 1) 的好处也在这里切换分布时只需改 makedist 一行参数估计和标准化代码完全复用。3.3 模拟样本量怎么定模拟样本量 n_mc 决定了分位数估计的稳定性。样本太少尾部分位数每次运行都在变样本太多计算成本线性上升收益却很小。我常用的参考量级如下n_mc0.5% 分位数相对波动适用场景1000±15% 左右快速探索不追求稳定10000±4% ~ ±6%单变量数据通用推荐100000±1.5% ~ ±2%多维或高显著水平这里的相对波动是不同随机种子下同一样本量分位数变化的经验范围。对单变量异常值剔除10000 个模拟样本足够如果原始数据本身只有几百个点模拟样本再多也无法提升参数估计精度瓶颈在 MLE 而不在抽样。如果数据包含明确的上下物理边界比如传感器量程我会把模拟样本直接截断到量程内再计算分位数这样包络不会延伸到物理上不可能的区域。4. 异常点识别、迭代剔除与收敛判据4.1 用蒙特卡洛包络识别潜在异常点有了模拟样本的经验分布下一步把原始观测值与模拟包络比较。这里采用双尾 1% 分位水平即 0.5% 和 99.5% 分位数比常规 3σ 更严格适合剔除阶段使用。alpha 0.01; lo quantile(z_mc, alpha / 2); hi quantile(z_mc, 1 - alpha / 2); % 原始数据标准化 z_obs (data - mu_hat) / sigma_hat; % 超出包络即为潜在异常 flag (z_obs lo) | (z_obs hi); idx_out find(flag); fprintf(首次识别潜在异常点 %d 个\n, length(idx_out)); % 包络可视化 figure(Color, w); hold on; scatter(1:n, data, 20, [0.5 0.5 0.5], filled); plot(1:n, mu_hat lo * sigma_hat, r--, LineWidth, 1.5); plot(1:n, mu_hat hi * sigma_hat, r--, LineWidth, 1.5); xlabel(样本序号); ylabel(观测值); legend({原始数据, MC 下界, MC 上界}, Location, best); hold off;核心逻辑是计算模拟样本标准化值的分位数作为包络边界再代入原始观测值超出边界的标记为潜在异常。红色虚线是包络不是回归线表示“在给定分布假设下正常观测值应该落在这个范围内”。alpha 直接影响识别结果。0.01 偏保守适合样本量大、异常占比低的数据0.05 更激进适合异常占比高、需要强清洗的场景。quantile 接受 0 到 1 的比例prctile 接受 0 到 100这两个函数经常被混用我习惯统一用 quantile改分位水平时不用换算。4.2 迭代剔除每剔一点就重新估计参数一次性剔除所有超包络点是危险动作异常值参与均值和标准差估计包络本身被拉宽部分真异常可能藏在包络内部。常见做法是逐轮剔除每轮只用当前干净数据重新做 MLE 和蒙特卡洛模拟再剔除超出新包络的点直到没有新异常或达到迭代上限。data_clean data; removed_idx []; max_iter 10; n_mc 20000; for iter 1:max_iter n_now length(data_clean); if n_now 10 warning(剩余样本过少停止迭代); break; end % 当前干净集上的分布估计 [mu_i, sigma_i] normfit(data_clean); if sigma_i 1e-8 break; end % 生成蒙特卡洛样本并计算双尾阈值 mc_i random(... makedist(Normal, mu, mu_i, sigma, sigma_i), n_mc, 1); z_i (mc_i - mu_i) / sigma_i; lo_i quantile(z_i, alpha / 2); hi_i quantile(z_i, 1 - alpha / 2); % 检查当前数据 z_cur (data_clean - mu_i) / sigma_i; bad find(z_cur lo_i | z_cur hi_i); if isempty(bad) fprintf(第 %d 轮收敛无新异常\n, iter); break; end removed_idx [removed_idx; bad]; data_clean(bad) []; fprintf(第 %d 轮: 剔除 %d 个点剩余 %d 个\n, ... iter, length(bad), length(data_clean)); end循环结构是估计参数 → 生成模拟样本 → 计算包络 → 标记越界点 → 有点则剔除并进入下一轮。注意每轮都要重新拟合不能沿用最初的 mu_hat 和 sigma_hat否则掩蔽效应会一直存在。注意迭代循环里不要复用第一轮估计的 mu_hat 和 sigma_hat 作为每轮判定参数掩蔽效应会随着异常值残留累积导致永远收不了敛。停止条件有三个没有新异常被标记、剩余样本少于 10 个、达到最大迭代次数。前两个是数据层面的收敛第三个是保护性上限防止异常占比极高时循环把样本清空。迭代轮数超过 5 轮时我会打一条日志看每轮剔除量是否递减如果第二轮剔除数量比第一轮还多说明分布假设可能有问题先停手排查分布。4.3 剔除前后统计量与模型稳定性对比剔除完成后需要验证。验证对象不是“剔得多干净”而是“统计量是否稳定下来”。我对比剔除前后的均值、标准差、偏度和峰度并重新画箱线图。mu_after mean(data_clean); sig_after std(data_clean); fprintf(剔除前: n%d, mean%.3f, std%.3f, skew%.2f, kurt%.2f\n, ... n, mu0, sig0, skewness(data), kurtosis(data)); fprintf(剔除后: n%d, mean%.3f, std%.3f, skew%.2f, kurt%.2f\n, ... length(data_clean), mu_after, sig_after, ... skewness(data_clean), kurtosis(data_clean)); figure(Color, w); subplot(1, 2, 1); boxplot(data, Symbol, r); title(剔除前); subplot(1, 2, 2); boxplot(data_clean, Symbol, r); title(剔除后); % 交叉验证Grubbs 检验新版 Statistics Toolbox 提供 if exist(grubbs, file) [h_grubbs, idx_grubbs] grubbs(data, 0.05); fprintf(Grubbs 标记异常点 %d 个\n, sum(h_grubbs)); end判定要落到量级上均值变化不大但标准差显著缩小说明异常值主要拉宽了分布没有造成系统性偏移均值也明显变化时说明异常值的分布不对称剔除后的分析结论需要重写。比如一组模拟数据的对比指标剔除前剔除后变化方向样本量500486减少 2.8%均值12.7412.31下降 0.43标准差3.912.68下降 31.5%偏度1.820.34趋近对称峰度8.753.21趋近正态标准差下降 30% 以上是常见信号说明原本数据的变异主要来自离群点而不是主体分布。grubbs 函数在较新的 Statistics Toolbox 中提供旧版本可以手写 G 统计量并用 tinv 查临界值。蒙特卡洛和 Grubbs 同时标记的点剔除置信度更高只被一方标记的点保留下来进 review 名单。5. 蒙特卡洛剔除异常值的阈值调优与业务收口蒙特卡洛剔除异常值的参数里alpha 是最敏感的一个。我通常先跑一轮 alpha0.05 的初筛统计潜在异常占比若占比超过 10%说明数据质量比预想差不能单靠蒙特卡洛要回去查采集链路若占比在 2% 到 5% 之间再改用 alpha0.01 跑正式迭代。这个两段式策略比直接跑一个固定阈值稳定得多。分布假设不成立时优先换 t location-scale 分布。makedist 支持 tLocationScale拟合函数是 fitdist(data, tLocationScale)。相比正态分布它的尾部更厚模拟包络更宽误杀率会明显下降。另一个思路是做分布无关的 bootstrap对干净样本做有放回重采样每次计算均值或其他统计量用重采样分布替代解析分布。这类方法在 LakeArea 这种面积模拟、PortSim 这类投资组合模拟里也常用——先剔除异常观测再做蒙特卡洛预测预测区间的稳定性会好很多。迭代剔除在小样本下存在一个隐蔽风险每轮剔除后样本量变小MLE 方差变大包络可能收缩下一轮误杀更多点。我一般会加一个保护条件单轮剔除比例超过当前样本量的 5% 时停止迭代并输出警告让业务方确认这批数据是否真的异常占比这么高。数学上这是合理的保守条件工程上是防呆设计。最后给一个落地技巧不要把剔除结果直接覆盖原数据而是在原表上加一列 mc_label取值 keep、remove 或 review。蒙特卡洛包络外且与 Grubbs 交叉验证通过的点标为 remove包络边界附近z 值落在 0.9% 到 1% 分位区间标为 review其余标为 keep。这个三分类方案比一味删除更适合真实数据分析流程——保留原始记录、标注去向、留给下游业务判断才是蒙特卡洛异常值剔除在生产环境里的正确形态。本文还有配套的精品资源点击获取
