企业官网建设流程全解析

简介：面向数据预处理与异常检测需求，这份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小于等于变量数p，cov(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 数据矩阵，n>=p+2 % 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的分位数，所以alpha=0.05等同于 95% 覆盖。工程上常见的alpha是 0.025 或 0.01，因为异常值往往是少数，拒绝域太大会误删正常点。有一个思路是先用较小的alpha剔除强异常，再对剩余数据重新估计协方差，这是迭代剔除的雏形。

4. 阈值与协方差估计：三个影响剔除结果的关键参数

4.1 置信度 alpha：0.975 还是 0.99

alpha是卡方分布的分位数，不是实际异常比例。如果你知道数据中大约有 5% 的异常，就把阈值设到 95% 分位数附近；如果异常比例很低，建议用 99% 分位数。实际操作中可以画一下D2的直方图，看尾部从哪里开始显著脱离卡方曲线。我自己经常在 0.01 和 0.05 之间做敏感性分析：如果剔除结果对alpha剧烈变化，说明数据中异常样本还不是明显偏离总体，需要回到特征工程层面。

4.2 样本量与维度比：协方差矩阵的稳定性

这是马氏距离最大的一道坎。当n和p的比值小于 2.5 时，cov(X)本身噪声太大，马氏距离的有效性会快速下降。比如一个 40 行 20 列的数据集，协方差矩阵需要估计p(p+1)/2个独立参数，也就是 210 个值，但样本只有 40 个，估计结果严重过拟合，inv(Sigma)会把微小噪声放大成巨大的距离值。应对方式有三种倾向：第一种是做特征选择，保留最重要的变量；第二种是使用正则化协方差，比如 Ledoit-Wolf 收缩估计，MATLAB 里cov(X)没有内置参数，但可以自己写收缩公式；第三种是改用基于马氏距离的稳健版本，也就是robustcov，它通过子集抽样避免协方差被异常点污染。

4.3 稳健估计：用 robustcov 解决掩蔽效应

当异常值本身数量不多但幅度很大时，普通cov估计出的协方差远大于真实总体协方差，导致所有点看起来都接近中心，马氏距离失效。robustcov基于 MCD 算法，它先寻找一个子集，使得子集样本的协方差行列式最小，再用这个子集的均值和协方差计算距离。这能有效避免掩蔽效应，但缺点是计算量大，数据量超过几万行时很吃内存。下面是一个对比示例：

代码上，用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); % 卡方阈值，alpha=0.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;

运行后，你会看到红色叉号集中在数据云的边缘，而且被挤向相关方向上偏出去的区域，而不是单纯取变量绝对值的极值。这里的alpha=0.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。

本文还有配套的精品资源，点击获取