ARTICLE DETAIL

资讯详情

深耕商务建站与企业官网运营的一线实战洞察。

MATLAB中的马氏距离:从原理到实现异常值检测与数据清洗

MATLAB中的马氏距离:从原理到实现异常值检测与数据清洗 简介面向数据预处理与异常检测需求这份MATLAB源码实现了基于马氏距离的异常样本剔除方法。相比欧氏距离马氏距离充分考虑了特征间的相关性在多元统计分析与机器学习建模前清洗异常值时更为可靠。压缩包内含2个文件一个m脚本用于计算均值、协方差矩阵并输出马氏距离一个mat数据文件可直接加载测试整体仅73KB轻量易用。已有2258人学习过该资源适合需要快速上手异常值检测的MATLAB使用者。通过源码演示的完整流程读者可掌握从数据预处理、阈值设定到迭代剔除异常的思路并迁移到自己的数据集中提升模型稳定性。1. 马氏距离为什么是异常值检测的利器做多变量数据清洗时我经常遇到一种尴尬变量两两之间有强相关量纲还差着几个数量级这时候用欧氏距离做异常值筛选结果往往被单位最大的变量牵着走。马氏距离的核心思路是先把数据投影到“标准化”的空间再算距离它同时考虑了变量本身的方差和变量之间的协方差所以对二维平面上一团“斜着的椭圆”数据马氏距离能给出远比欧氏距离合理的异常判定。这个标题提到的“剔除异常样本”和“检测异常值”本质是同一件事的两种说法先用马氏距离给每个样本打分再按一个阈值把尾部样本挑出来。这套方法适合做光谱数据、工业传感器多通道信号、财务指标等场景的预处理也适合刚接触多元统计的 MATLAB 用户快速落地。2. MATLAB 中马氏距离的计算mahal 与手动实现2.1 马氏距离的定义与直觉马氏距离本质上是一个带权重的欧氏距离。给定均值向量mu和协方差矩阵Sigma样本x到总体的马氏距离平方定义为D^2 (x - mu) * inv(Sigma) * (x - mu)2.1.1 公式拆解(x - mu)是把数据中心化inv(Sigma)是对协方差矩阵求逆相当于把椭球形分布压回球形。如果Sigma退化为单位矩阵马氏距离就等于欧氏距离。当变量之间存在相关性时协方差矩阵的非对角元素会改变距离的计算方向——两个变量同步变化不会被视为“异常”只有偏离这个相关性结构时才被突出。这正是它适合异常值检测的根本原因。2.1.2 与欧氏距离的对比很多刚用 MATLAB 的人会用pdist2或直接sqrt(sum((x - mu).^2, 2))算距离。但举个例子一个温度传感器和一个压力传感器温度标准差是 10 度压力标准差是 0.5 MPa欧氏距离会把温度波动当成主要误差源压力通道的微小偏移完全被淹没。马氏距离用协方差做了归一化两个通道的贡献量级一致异常点更容易被识别。2.2 MATLAB 内置函数 mahal 的使用MATLAB 统计与机器学习工具箱里提供了mahal函数这是最省事的路子。2.2.1 最小实现假设你有一个n x p的数据矩阵X想计算每个样本相对整个数据集的马氏距离平方% 生成一个示例数据矩阵200行3列 X randn(200, 3); X(:, 2) X(:, 1) * 0.7 0.3 * randn(200, 1); % 让前两列相关 % 计算每个样本到总体均值的马氏距离平方 D2 mahal(X, X);2.2.2 mahal 返回的是距离平方这里的D2是每个样本的马氏距离平方不是距离本身。为什么返回平方因为平方后服从卡方分布方便直接用chi2inv定阈值。如果你需要距离值自己加一行D sqrt(D2)即可。mahal函数的典型坑有两个一是要求X的列数大于 1二是X的行数必须大于列数否则协方差矩阵不可逆函数会直接报错或者给出NaN。2.3 手动实现马氏距离的代价与收益有些时候你不用内置函数比如需要在 Simulink 里实时计算或者想完全控制协方差的估计方式。手动实现也不复杂mu mean(X, 1); Sigma cov(X); invSigma inv(Sigma); D2_manual zeros(size(X, 1), 1); for i 1:size(X, 1) dx X(i, :) - mu; D2_manual(i) dx * invSigma * dx; end用inv在小规模数据上没什么问题但当p接近样本数时inv(Sigma)极不稳定。实际工程中我通常用pinv求伪逆或者直接改用robustcov这一点在后面章节展开。手动实现的好处是你能在dx * invSigma * dx这一行清楚看到马氏距离的构成也能插入日志调试验证数据形状。代价是循环求值慢数据量大时可以用sum((X - mu) * invSigma .* (X - mu), 2)向量化代替。3. 用马氏距离剔除异常样本的可运行流程3.1 剔除异常样本的完整步骤这里给出一个标准化流程基本适用于大多数表格型数据。第一步是整理数据保证每一行是一个样本每一列是一个变量变量之间必须是连续数值。第二步是估计均值和协方差通常用mean(X)和cov(X)。第三步是计算每个样本的马氏距离平方。第四步是确定阈值推荐使用卡方分布的上分位点。第五步是把距离超过阈值的样本标记为异常然后剔除或替换。3.1.1 数据形状要求如果样本数n小于等于变量数pcov(X)是奇异矩阵马氏距离直接失效。这种情况下需要先降维或者用正则化协方差估计。对很多高维场景比如基因表达谱或高光谱数据直接用mahal是行不通的。一个常见做法是先用 PCA 把维度压到主成分个数小于样本数再对主成分分数计算马氏距离。但需要注意 PCA 本身对异常值敏感异常样本会影响主成分方向。3.1.2 估计协方差矩阵的细节协方差矩阵的估计方法直接影响判别效果。普通cov使用简单算术平均如果样本中存在离群点这些点会“拉大”协方差结果可能让真正的大偏差看起来不极端这就是所谓的掩蔽效应。解决思路是使用稳健协方差估计比如 MCD最小协方差行列式。MATLAB 里robustcov函数就是基于 MCD后面会有例子。3.2 可运行的 MATLAB 函数下面这个函数可以直接复制保存为removeOutliersByMahal.m输入数据X和显著性水平alpha输出剔除后的矩阵和异常索引。function [X_clean, outlierIdx] removeOutliersByMahal(X, alpha) % 输入 % X - n x p 数据矩阵np2 % alpha - 显著性水平默认 0.05 % 输出 % X_clean - 剔除异常后的数据 % outlierIdx - 异常样本的行索引 if nargin 2 || isempty(alpha) alpha 0.05; end % 检查数据形状 [n, p] size(X); if n p error(样本数必须大于变量数当前 cov 矩阵奇异); end % 计算马氏距离平方 D2 mahal(X, X); % 卡方分布阈值自由度等于变量数 p threshold chi2inv(1 - alpha, p); % 标记异常 outlierIdx find(D2 threshold); X_clean X; X_clean(outlierIdx, :) []; % 删除异常行 end这里mahal(X, X)有一个细节第二参数X被当作参考总体函数内部会用mean(X)和cov(X)作为均值向量和协方差矩阵。如果你有一个干净的参考样本集Xref想用它对新的Xnew打分应该写成mahal(Xnew, Xref)这更符合实际生产中的“训练/测试分离”思路。阈值使用chi2inv是因为在多元正态假设下马氏距离平方服从自由度为p的卡方分布。如果数据明显不是正态卡方阈值会偏保守或偏激进这时可以考虑基于经验分布取 99% 分位数作为阈值。3.3 阈值确定卡方分布与经验分位数的取舍3.3.1 为什么要用卡方分布多元正态分布有一个已知结论样本到总体中心的马氏距离平方服从卡方分布自由度是变量数。因此chi2inv(0.95, p)能给出一个理论上的 95% 覆盖范围。这个结论在小样本时并不特别精确当n在 50 以下尤其是p接近n时卡方阈值会低估异常比例导致异常样本漏检。这时候我倾向于用经验分布直接取D2的 97.5% 分位数作为阈值。但对小样本极端值会影响分位数估计所以没有绝对安全的选择。3.3.2 chi2inv 的用法chi2inv是统计工具箱的函数第一个参数是累积概率值第二个参数是自由度。比如chi2inv(0.99, 5)返回 5 个自由度下卡方分布 99% 分位数。注意显著性水平alpha与分位数的关系阈值取1 - alpha的分位数所以alpha0.05等同于 95% 覆盖。工程上常见的alpha是 0.025 或 0.01因为异常值往往是少数拒绝域太大会误删正常点。有一个思路是先用较小的alpha剔除强异常再对剩余数据重新估计协方差这是迭代剔除的雏形。4. 阈值与协方差估计三个影响剔除结果的关键参数4.1 置信度 alpha0.975 还是 0.99alpha是卡方分布的分位数不是实际异常比例。如果你知道数据中大约有 5% 的异常就把阈值设到 95% 分位数附近如果异常比例很低建议用 99% 分位数。实际操作中可以画一下D2的直方图看尾部从哪里开始显著脱离卡方曲线。我自己经常在 0.01 和 0.05 之间做敏感性分析如果剔除结果对alpha剧烈变化说明数据中异常样本还不是明显偏离总体需要回到特征工程层面。4.2 样本量与维度比协方差矩阵的稳定性这是马氏距离最大的一道坎。当n和p的比值小于 2.5 时cov(X)本身噪声太大马氏距离的有效性会快速下降。比如一个 40 行 20 列的数据集协方差矩阵需要估计p(p1)/2个独立参数也就是 210 个值但样本只有 40 个估计结果严重过拟合inv(Sigma)会把微小噪声放大成巨大的距离值。应对方式有三种倾向第一种是做特征选择保留最重要的变量第二种是使用正则化协方差比如 Ledoit-Wolf 收缩估计MATLAB 里cov(X)没有内置参数但可以自己写收缩公式第三种是改用基于马氏距离的稳健版本也就是robustcov它通过子集抽样避免协方差被异常点污染。4.3 稳健估计用 robustcov 解决掩蔽效应当异常值本身数量不多但幅度很大时普通cov估计出的协方差远大于真实总体协方差导致所有点看起来都接近中心马氏距离失效。robustcov基于 MCD 算法它先寻找一个子集使得子集样本的协方差行列式最小再用这个子集的均值和协方差计算距离。这能有效避免掩蔽效应但缺点是计算量大数据量超过几万行时很吃内存。下面是一个对比示例估计方式适用场景缺点推荐用途普通 cov数据干净、异常比例低于 1%对异常敏感可能漏检快速初筛稳健 MCD异常比例 10%-20%且无明显规律计算慢需要统计工具箱正式建模前的清洗收缩估计高维小样本p 接近 n需要选择收缩强度基因、光谱数据代码上用robustcov替换普通协方差通常配合计算稳健马氏距离平方[sigmaRob, muRob, w2, mahDistRob] robustcov(X); % sigmaRob 为稳健协方差muRob 为稳健均值向量 % mahDistRob 为稳健马氏距离平方等价于对 X 中的每行计算注意robustcov的第三个输出w2是每个样本的权重可用作异常程度评分。权重大于 0.5 的样本通常被认为是正常点这个经验值在不少工程场景中有效。使用稳健估计后阈值依然可以用卡方分布但自由度仍然是p因为理论分布没有变。5. 实战MATLAB 多元数据异常检测脚本与 CSV 接入5.1 准备模拟数据与噪声注入为了完整演示剔除流程这里生成一个含相关性的三维数据集并注入少量异常点。实际使用时你可以用readtable或readmatrix把 CSV 数据导入替换这里的模拟部分。模拟数据的关键是让两列之间存在线性关系这样才能体现马氏距离相对于欧氏距离的优势。5.2 完整可运行脚本% detectOutliersDemo.m % 生成带相关性的三维数据注入异常点用马氏距离剔除 rng(1); % 固定随机种子便于复现 n 200; p 3; % 基础数据第一列是标准正态第二列与第一列相关第三列独立 X randn(n, p); X(:, 2) 0.8 * X(:, 1) 0.6 * randn(n, 1); % 注入 10 个异常点把前 10 行的值整体偏移 X(1:10, :) X(1:10, :) [4, 3, 2]; % 导入外部 CSV 的接法如果数据已存在 % data readmatrix(sensor_data.csv); % X data(:, 1:3); % 计算马氏距离平方 D2 mahal(X, X); % 卡方阈值alpha0.02自由度 3 alpha 0.02; threshold chi2inv(1 - alpha, p); % 标记异常 outlierIdx find(D2 threshold); % 剔除异常并输出结果 X_clean X; X_clean(outlierIdx, :) []; fprintf(总样本数%d\n, n); fprintf(检出异常点数%d\n, length(outlierIdx)); fprintf(理论阈值%.2f\n, threshold); % 绘图对比前两个变量散点图正常点与异常点用不同颜色 figure; scatter(X(:, 1), X(:, 2), 20, k, filled); hold on; scatter(X(outlierIdx, 1), X(outlierIdx, 2), 80, r, x); legend({正常样本, 异常样本}, Location, best); xlabel(变量1); ylabel(变量2); title(马氏距离检测异常值结果); grid on;运行后你会看到红色叉号集中在数据云的边缘而且被挤向相关方向上偏出去的区域而不是单纯取变量绝对值的极值。这里的alpha0.02意味着预期误判率约为 2%。如果注入异常偏离强度更大检出率会更高。如果你发现异常点没有被正确分离先检查是否数据中存在缺失值mahal遇到NaN会直接让整个协方差矩阵崩溃。5.3 结果解释与参数调整5.3.1 观察距离排序不要只看阈值把D2从大到小排序取前 20 个索引观察。如果前 10 个恰好是注入的异常后 10 个是正常边界点说明阈值偏严。这种情况下把alpha调到 0.05或者改为取距离排序的后 5% 作为异常都会改变最终清洗后的数据分布。建议在剔除前先保存一份D2变量随后画一个距离分布直方图和卡方概率密度曲线叠加对比目视检查尾部是否一致。5.3.2 导出剔除后的数据writematrix可以避免手工复制writematrix(X_clean, X_clean.csv);注意这里覆盖了原文件内容所以运行时先确认路径。生产环境中我会把异常索引存成outlierIdx.csv保留原始数据而不直接删除方便溯源。马氏距离剔除的局限性在于如果异常是以局部模式出现比如某个传感器只在一段时间内失效那么距离本身难以区分“正常变异性”和“故障偏移”。此时可以考虑对时间序列加滑动窗口在每个窗口内计算局部马氏距离再对距离序列做趋势分析。6. 进阶稳健协方差估计与异常值可视化验证当数据中已经混入一批异常值普通mahal的协方差估计会被污染导致距离分数偏低。此时可以改用robustcov得到稳健距离并配合 Q-Q 图做验证。先看稳健版本的核心调用[sigmaRob, muRob, ~, D2Rob] robustcov(X); thresholdRob chi2inv(0.99, p); outlierRob D2Rob thresholdRob;robustcov的默认方法是用 Fast-MCD 算法它会从样本中随机抽取子集迭代计算因此结果带有随机性。建议设置随机种子以重复实验或者多次运行收集异常索引的并集。robustcov的第四输出已经是稳健马氏距离平方不需要再手动减均值乘逆矩阵。验证手段之一是画卡方 Q-Q 图把距离平方排序后与卡方分布的分位数做散点。正常数据应该大致落在直线附近右上方明显翘起的点就是异常候选。MATLAB 里没有直接的卡方 Q-Q 图函数可以这样生成% 生成理论分位数 p_seq (1:n) / (n 1); theoretical chi2inv(p_seq, p); % 对 D2 排序后画散点 D2_sorted sort(D2); plot(theoretical, D2_sorted, o); hold on; plot(theoretical, theoretical, k--); % 参考对角线 xlabel(卡方理论分位数); ylabel(马氏距离平方排序后);第二点经验是固定异常比例。卡方阈值适合正态数据但工程数据总带偏态我习惯用prctile取 95% 分位数作为阈值这样不用反复调alpha。但注意这种方法输出的异常数量和比例是预设的可能错把边界点圈进来。更保险的做法是对D2做对数变换再对变换后的数据用 3σ 法则因为对数变换后的极端值更接近对称分布。最后一个实操细节mahal和robustcov都会因变量单位不同而得到相同的距离因为协方差矩阵吸收了尺度信息。但这不意味着数据不需要预处理。当某个变量的方差极小比如接近机器精度时协方差矩阵中对应行列接近零逆矩阵放大该维度上的微小偏差本来正常的测量噪声会被误判为异常。处理方法是先剔除方差接近零的变量或者用zscore标准化后再计算马氏距离。两种做法会得到几乎一样的结果但标准化后的协方差矩阵数值上更稳定也能避免mahal因为矩阵病态返回Inf。本文还有配套的精品资源点击获取
返回列表
PREV
查看更多资讯
NEXT
返回资讯列表