
简介面向数据分析与统计建模学习者这份资料以MATLAB为工具系统演示如何用蒙特卡洛模拟识别并剔除异常值。内容涵盖数据预处理、随机样本生成、统计检验Z-score、IQR等到迭代剔除与结果验证的完整流程配有可直接运行的m脚本可用于科研数据处理、实验数据清洗等场景。压缩包共19个文件其中15个MATLAB脚本构成主要代码主体2个PPT分别以英文和法文讲解蒙特卡洛模拟原理另含1个mat数据文件和1个txt说明文档整体仅388KB轻量便于下载。已有1080人学习下载。资源附加了MyMC、方差缩减、湖泊面积估计、港口模拟等演示示例帮助读者在具体任务中体会蒙特卡洛方法的应用技巧PPT与脚本配合既适合教师备课引用也适合自学时对照演练。1. 蒙特卡洛异常值剔除为什么 3σ 和 IQR 不够用一组传感器数据里混进几个离群点第一反应多半是 3σ 或 IQR 看一眼。这两个方法在单峰对称分布上有效可生产数据很少长成那样流量曲线有峰谷、交易量有尖峰、模型的训练样本里往往混杂半正常半噪声的多模态分布。条件一偏离3σ 会把正常业务波动剔除掉IQR 对偏态分布则基本无能为力。蒙特卡洛异常值剔除的思路完全不同——不假定数据服从某个参数族而是对样本做大量随机重采样把均值、协方差这类统计量的经验分布估出来再基于这个分布计算每个样本的马氏距离和判定分位点。对没有先验分布、却需要可解释剔除结果的数据清洗任务这比硬套正态公式稳得多。数据工程、算法、量化分析相关岗位的人都能在这套流程里找到自己需要的步骤下面的 Python 实现可以直接改来用。2. 蒙特卡洛剔除异常样本的原理重采样、稳健统计量与距离度量2.1 重采样为什么能代替分布假设蒙特卡洛异常值剔除的内核是经验分布思想。对数据集做 B 次有放回抽样每个子样本都独立计算中心统计量B 次之后就得到了这些统计量的经验分布。你不需要提前假设数据是正态的或对称的只要重采样次数足够bootstrap 分布就会稳定逼近真实抽样分布的形状判定阈值完全从数据自身来。传统 3σ 用原始数据的全局均值和全局标准差卡一刀异常点一旦存在均值和标准差自身就被带偏出现 Masking Effect掩蔽效应。换到蒙特卡洛框架里每次采样出来的子样本相当于一个带随机噪声的估计B 次之后再做统计量聚合异常点对中心的贡献被稀释到极小。蒙特卡洛算法家族里重采样路径在处理异常样本剔除问题时能站住脚就是这个原因。要顺带区分的是动力学蒙特卡洛那种处理状态转移轨迹的蒙特卡洛算法跟这里完全无关异常值剔除只用“重复抽样 统计量聚合”这件事。最小实现长这样import numpy as np def bootstrap_center(X, B1000, ratio0.8, seed42): 用蒙特卡洛重采样估计稳健中心和协方差矩阵。 rng np.random.default_rng(seed) n X.shape[0] m max(1, int(n * ratio)) mus np.zeros((B, X.shape[1])) sigmas np.zeros((B, X.shape[1], X.shape[1])) for i in range(B): idx rng.integers(0, n, m) sub X[idx] mus[i] sub.mean(axis0) sigmas[i] np.cov(sub.T) # 逐元素取中位数后加正则化项防止协方差矩阵退化 sigma np.median(sigmas, axis0) 1e-8 * np.eye(X.shape[1]) return np.median(mus, axis0), sigmaratio 控制在 0.6 到 0.9 之间0.8 是常用起点子样本太小统计量方差偏大子样本太大异常点进入子样本的概率升高中心估计的稳健性就会被削弱。B 在 1000 到 5000 之间1000 次经验分布的尾部分位点已经足够平滑再往上主要增加运行时间。seed 固定写死结果才可复现线上任务里建议把 seed 值写进日志。2.2 中位数聚合为什么稳健中心是干净数据的前提知道要对 mu 取中位数的人不少忽略对协方差做同样抗扰的人更多。如果直接把各子样本协方差取平均一个恰好包含极端异常点的子样本会把协方差整体拉大马氏距离的分母被撑大极端样本反而显得“正常”整个剔除动作失效。常见的做法有两个方向一是像上面代码那样逐元素取中位数再给对角线加一个小值二是用几何中位数定义协方差。逐元素中位数简单、计算快但可能破坏正定性所以加正则化项兜底。中心和协方差是一体的这一环不稳健后续距离判定就会来回抖去掉的点换个 seed 可能又回来了。2.3 马氏距离与卡方分位点判定边界从哪来距离度量上欧氏距离假设各特征独立且同方差特征尺度一不同就完全没法用。马氏距离把协方差结构考虑进去相当于把数据白化再算距离相关性冗余被消掉公式是md_i sqrt((x_i - mu)^T Sigma^{-1} (x_i - mu))如果数据来自多维正态总体距离平方近似服从 p 维卡方分布所以判定可以写成from scipy.stats import chi2 cutoff np.sqrt(chi2.ppf(1 - alpha, dfp))这里的 alpha 是名义显著性水平不是异常样本占比。把三种方法的差异摆在一起看方法中心估计尺度估计判定依据对异常单点的敏感性3σ全局均值全局标准差|z|3极敏感掩蔽效应强IQR中位数四分位距1.5×IQR 越界中等但只看单特征蒙特卡洛马氏bootstrap 稳健中心bootstrap 稳健协方差卡方分位点低联合分布一致3σ 的均值和标准差会被异常点拉偏这是掩蔽效应的根源。IQR 用中位数和四分位距缓解了中心估计问题但判定只沿单特征方向做相关特征之间的联合异常基本漏检。蒙特卡洛异常样本剔除把中心估计、协方差估计、分位点判定揉成一个整体相关结构场景下的表现比固定阈值稳定。3. 用 Python 写蒙特卡洛异常值剔除最小可复现代码3.1 完整实现一个能直接接进流水线的剔除器把前面的原理组装成一个类方便在各处就地调用import numpy as np from scipy.stats import chi2 def mahalanobis_distance(X, mu, sigma): inv np.linalg.pinv(sigma) dx X - mu return np.sqrt(np.einsum(ij,jk,ik-i, dx, inv, dx)) class MonteCarloOutlierRemover: 基于蒙特卡洛重采样的异常样本剔除器。 def __init__(self, B1000, alpha0.05, ratio0.8, max_iter5, seed42): self.B B self.alpha alpha self.ratio ratio self.max_iter max_iter self.seed seed def _center(self, X): rng np.random.default_rng(self.seed) n, p X.shape m max(1, int(n * self.ratio)) mus np.zeros((self.B, p)) sigmas np.zeros((self.B, p, p)) for i in range(self.B): idx rng.integers(0, n, m) sub X[idx] mus[i] sub.mean(axis0) sigmas[i] np.cov(sub.T) mu np.median(mus, axis0) sigma np.median(sigmas, axis0) np.eye(p) * 1e-8 return mu, sigma def fit(self, X): X np.asarray(X, dtypefloat) mask np.ones(X.shape[0], dtypebool) for _ in range(self.max_iter): cur X[mask] mu, sigma self._center(cur) md mahalanobis_distance(cur, mu, sigma) cutoff np.sqrt(chi2.ppf(1 - self.alpha, dfcur.shape[1])) new_mask mask.copy() new_mask[mask] md cutoff if new_mask.sum() mask.sum(): break mask new_mask return maskfit 的流程是第一轮全量估计之后迭代剔除。每轮重新算一次稳健中心和协方差再算马氏距离和卡方阈值直到某轮没有样本被标记为异常就收敛退出。设计成迭代是为了处理中心偏移异常点离群明显时第一轮剔除留下的数据中心还会漂移一次性定阈值容易误伤。3.2 代码逻辑与参数说明fit 里的参数按这张表选参数推荐范围作用调大影响调小影响B1000-5000重采样次数分位点更平滑但更慢偶然抖动增大alpha0.01-0.10名义误判概率剔除范围扩大剔除范围收缩ratio0.6-0.9每次抽样占比更接近全量分布统计量方差更大max_iter3-8最大迭代轮数收敛更充分但耗时聚团异常可能剔不净数据量在 500 条以下时B 降到 500、alpha 用 0.01 更稳否则 5% 的误判率在样本少时会误删很多有效数据。数据量过万时 B1000 就够ratio 提到 0.9运行时间主要耗在 np.covB 加大的边际收益很小。3.3 对照实验3σ、IQR 与蒙特卡洛在生成数据上的表现写个小实验验证算法逻辑rng np.random.default_rng(0) normal rng.multivariate_normal([0, 0], [[1, 0.5], [0.5, 1]], 500) outliers rng.multivariate_normal([8, 8], [[0.5, 0], [0, 0.5]], 20) data np.vstack([normal, outliers]) # 3σ mu, sig data.mean(0), data.std(0) flag_3sigma np.all(np.abs(data - mu) 3 * sig, axis1) # IQR q1, q3 np.percentile(data, [25, 75], axis0) iqr q3 - q1 flag_iqr np.all((data q1 - 1.5 * iqr) (data q3 1.5 * iqr), axis1) # 蒙特卡洛 remover MonteCarloOutlierRemover(B1000, alpha0.05, seed1) flag_mc remover.fit(data) for name, flag in [(3σ, flag_3sigma), (IQR, flag_iqr), (蒙特卡洛, flag_mc)]: hits int((~flag[500:]).sum()) print(f{name}: 剔除 {int((~flag).sum())} 个召回异常 {hits}/20)这个数据里异常点聚集在 (8,8) 附近与正常数据呈相关结构。3σ 的均值被拉到中心附近、标准差被放大异常点几乎全被漏检IQR 各维独立判断能捡回部分两个方向都突出的点但对沿相关方向的联合异常仍然漏蒙特卡洛加马氏距离用卡方分位点做边界召回最稳误伤数量也最低。这个实验的价值不在对比本身而在验证剔除流程的收敛性——mask 停止变化后结果不再依赖 max_iter。4. 实战踩坑蒙特卡洛异常样本剔除的参数边界4.1 三个必调参数B、alpha、max_iter 的真实含义第一坑是 B。B 小于 300 时经验分布尾部抖动明显同一个点在不同 seed 下标记结果可能不一致。出现这种情况不是算法错而是统计量还没稳定。线上任务固定 seed部署前用两三个不同 seed 跑一遍剔除集合如果相差超过 1%就得加大 B。B 超过 5000 边际收益接近零计算时间却线性上升。大数据集用 ratio0.9、B500效果接近 ratio0.8、B1000还省一半开销。第二坑是 alpha。alpha 不是数据里异常样本占比而是“把正常误判为异常”的概率。数据完全干净时 alpha0.05 也会剔掉约 5% 的边缘样本。生产上我一般把 alpha 调到 0.01然后把剔除索引存下来人工复核宁可少剔也不误伤。alpha 大于 0.1 时边缘分布会被成片剪掉损失的信息比异常还严重。第三坑是 max_iter。异常点聚团时单轮剔除只能摘掉离中心最远的点下一轮中心移动后又露出新远点。max_iter 小于 3聚团异常剔不干净迭代超过 8 还没收敛往往不是 max_iter 的事而是数据本身是多模态单个马氏距离拟合不了。4.2 偏态、多模态和高维数据形态决定模型选择偏态分布是最常见的坑。右偏数据里的正常长尾样本会被卡方分位点误杀因为卡方假设本身带正态遗传。解决办法是先把数据做对数或 Box-Cox 变换压掉尾部变换做完再进剔除器剔除完要把索引映射回原始空间不要拿变换后的值当成最终数据使用。多模态则要提前分裂。两个相距很远的正常簇被单一中心和协方差定义成异常几乎是必然的。常见做法是先聚类K 簇后每簇单独跑一遍异常样本剔除。这和全局异常剔除的目标不一样前者抓的是“该簇内部既然出现的离群点”后者抓的是“不属于任何一个簇的样本”得分开处理。高维情况需要额外小心。马氏距离要求样本量远大于维数样本量与维度比低于 10:1 时协方差估计方差剧增伪逆引入的数值噪声会盖过真实异常信号。把 np.cov 换成 LedoitWolf 收缩估计或先用 PCA 压到 20 维以内比硬调 B 和 alpha 有效得多。这个 10:1 的经验值在工程上是很好的预警线低于它就不要指望卡方分布能给出准确阈值。4.3 蒙特卡洛异常值剔除与交叉验证的配合交叉验证里最容易犯的是泄露错误。剔除器必须在训练折内拟合只用训练折剔剩下的数据训模型验证折原封不动拿来评估。from sklearn.model_selection import KFold kf KFold(n_splits5, shuffleTrue, random_state0) for fold, (train_idx, val_idx) in enumerate(kf.split(X)): remover MonteCarloOutlierRemover(B1000, alpha0.05, seed0) clean train_idx[remover.fit(X[train_idx])] model.fit(X[clean], y[clean]) print(ffold {fold}: 剔除 {len(train_idx) - clean.sum()} 个)如果先在整个 X 上做剔除再切折验证集信息被中心估计间接用过了验证分数会乐观得没有参考意义。常见做法是时间窗口滚动每周重写一次 B 和 alpha而不是一次性定死再无更新。5. 验证蒙特卡洛异常值剔除效果AUC 评估与滑动窗口进阶5.1 用 ROC/AUC 检验剔除器的排序质量有真实异常标签时盯召回率还不够还得看误伤。马氏距离本身是个排序分数用 roc_auc_score 衡量正合适from sklearn.metrics import roc_auc_score mu, sigma remover._center(data[mask]) score mahalanobis_distance(data, mu, sigma) auc roc_auc_score(y_true, score) print(fAUC: {auc:.3f})AUC 高于 0.9 说明距离排序对异常和正常的区分力很强0.8 到 0.9 是可用区间低于 0.7 基本要回头查数据形态看是多模态还是偏态。比对过程中最重要的一点是标签和分数必须来自同一个时间窗口跨窗口验证会被概念漂移干扰造成评估失真。AUC 只反映排序质量真正剔除多少还要结合固定 alpha 看一版具体结果。5.2 三种进阶用法第一把剔除从布尔判断改成带分位数的输出。每轮蒙特卡洛子样本算出的距离可以累积成经验分布留给下游业务定自己的阈值而不是依赖卡方假设。这样剔除结果的解释从“它是异常”变成“它在 90 分位以上离 99 分位还有多少”可解释性高一个档次。第二对模型预测残差做蒙特卡洛异常值剔除。无论回归模型还是时序预测训练结束后都能拿残差向量做重采样。把“输入不像总体”转化为“预测误差不像总体”常用于抓取突变样本和概念漂移前的异常记录。第三滑动窗口的在线部署。窗口长度 1000、每步都跑 B1000 的重采样热路径撑不住。实际线路上把这套剔除器做成小时级批任务把剔除索引和距离分数推给实时推理端用离线重采样换实时稳定性是性价比最高的应用形态。本文还有配套的精品资源点击获取