簡介面向數(shù)據(jù)預(yù)處理與異常檢測需求這份MATLAB源碼實(shí)現(xiàn)了基于馬氏距離的異常樣本剔除方法。相比歐氏距離馬氏距離充分考慮了特征間的相關(guān)性在多元統(tǒng)計(jì)分析與機(jī)器學(xué)習(xí)建模前清洗異常值時(shí)更為可靠。壓縮包內(nèi)含2個(gè)文件一個(gè)m腳本用于計(jì)算均值、協(xié)方差矩陣并輸出馬氏距離一個(gè)mat數(shù)據(jù)文件可直接加載測試整體僅73KB輕量易用。已有2258人學(xué)習(xí)過該資源適合需要快速上手異常值檢測的MATLAB使用者。通過源碼演示的完整流程讀者可掌握從數(shù)據(jù)預(yù)處理、閾值設(shè)定到迭代剔除異常的思路并遷移到自己的數(shù)據(jù)集中提升模型穩(wěn)定性。1. 馬氏距離為什么是異常值檢測的利器做多變量數(shù)據(jù)清洗時(shí)我經(jīng)常遇到一種尷尬變量兩兩之間有強(qiáng)相關(guān)量綱還差著幾個(gè)數(shù)量級(jí)這時(shí)候用歐氏距離做異常值篩選結(jié)果往往被單位最大的變量牽著走。馬氏距離的核心思路是先把數(shù)據(jù)投影到“標(biāo)準(zhǔn)化”的空間再算距離它同時(shí)考慮了變量本身的方差和變量之間的協(xié)方差所以對(duì)二維平面上一團(tuán)“斜著的橢圓”數(shù)據(jù)馬氏距離能給出遠(yuǎn)比歐氏距離合理的異常判定。這個(gè)標(biāo)題提到的“剔除異常樣本”和“檢測異常值”本質(zhì)是同一件事的兩種說法先用馬氏距離給每個(gè)樣本打分再按一個(gè)閾值把尾部樣本挑出來。這套方法適合做光譜數(shù)據(jù)、工業(yè)傳感器多通道信號(hào)、財(cái)務(wù)指標(biāo)等場景的預(yù)處理也適合剛接觸多元統(tǒng)計(jì)的 MATLAB 用戶快速落地。2. MATLAB 中馬氏距離的計(jì)算mahal 與手動(dòng)實(shí)現(xiàn)2.1 馬氏距離的定義與直覺馬氏距離本質(zhì)上是一個(gè)帶權(quán)重的歐氏距離。給定均值向量mu和協(xié)方差矩陣Sigma樣本x到總體的馬氏距離平方定義為D^2 (x - mu) * inv(Sigma) * (x - mu)2.1.1 公式拆解(x - mu)是把數(shù)據(jù)中心化inv(Sigma)是對(duì)協(xié)方差矩陣求逆相當(dāng)于把橢球形分布?jí)夯厍蛐巍Ｈ绻鸖igma退化為單位矩陣馬氏距離就等于歐氏距離。當(dāng)變量之間存在相關(guān)性時(shí)協(xié)方差矩陣的非對(duì)角元素會(huì)改變距離的計(jì)算方向——兩個(gè)變量同步變化不會(huì)被視為“異?！敝挥衅x這個(gè)相關(guān)性結(jié)構(gòu)時(shí)才被突出。這正是它適合異常值檢測的根本原因。2.1.2 與歐氏距離的對(duì)比很多剛用 MATLAB 的人會(huì)用pdist2或直接sqrt(sum((x - mu).^2, 2))算距離。但舉個(gè)例子一個(gè)溫度傳感器和一個(gè)壓力傳感器溫度標(biāo)準(zhǔn)差是 10 度壓力標(biāo)準(zhǔn)差是 0.5 MPa歐氏距離會(huì)把溫度波動(dòng)當(dāng)成主要誤差源壓力通道的微小偏移完全被淹沒。馬氏距離用協(xié)方差做了歸一化兩個(gè)通道的貢獻(xiàn)量級(jí)一致異常點(diǎn)更容易被識(shí)別。2.2 MATLAB 內(nèi)置函數(shù) mahal 的使用MATLAB 統(tǒng)計(jì)與機(jī)器學(xué)習(xí)工具箱里提供了mahal函數(shù)這是最省事的路子。2.2.1 最小實(shí)現(xiàn)假設(shè)你有一個(gè)n x p的數(shù)據(jù)矩陣X想計(jì)算每個(gè)樣本相對(duì)整個(gè)數(shù)據(jù)集的馬氏距離平方% 生成一個(gè)示例數(shù)據(jù)矩陣200行3列 X randn(200, 3); X(:, 2) X(:, 1) * 0.7 0.3 * randn(200, 1); % 讓前兩列相關(guān) % 計(jì)算每個(gè)樣本到總體均值的馬氏距離平方 D2 mahal(X, X);2.2.2 mahal 返回的是距離平方這里的D2是每個(gè)樣本的馬氏距離平方不是距離本身。為什么返回平方因?yàn)槠椒胶蠓目ǚ椒植挤奖阒苯佑胏hi2inv定閾值。如果你需要距離值自己加一行D sqrt(D2)即可。mahal函數(shù)的典型坑有兩個(gè)一是要求X的列數(shù)大于 1二是X的行數(shù)必須大于列數(shù)否則協(xié)方差矩陣不可逆函數(shù)會(huì)直接報(bào)錯(cuò)或者給出NaN。2.3 手動(dòng)實(shí)現(xiàn)馬氏距離的代價(jià)與收益有些時(shí)候你不用內(nèi)置函數(shù)比如需要在 Simulink 里實(shí)時(shí)計(jì)算或者想完全控制協(xié)方差的估計(jì)方式。手動(dòng)實(shí)現(xiàn)也不復(fù)雜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在小規(guī)模數(shù)據(jù)上沒什么問題但當(dāng)p接近樣本數(shù)時(shí)inv(Sigma)極不穩(wěn)定。實(shí)際工程中我通常用pinv求偽逆或者直接改用robustcov這一點(diǎn)在后面章節(jié)展開。手動(dòng)實(shí)現(xiàn)的好處是你能在dx * invSigma * dx這一行清楚看到馬氏距離的構(gòu)成也能插入日志調(diào)試驗(yàn)證數(shù)據(jù)形狀。代價(jià)是循環(huán)求值慢數(shù)據(jù)量大時(shí)可以用sum((X - mu) * invSigma .* (X - mu), 2)向量化代替。3. 用馬氏距離剔除異常樣本的可運(yùn)行流程3.1 剔除異常樣本的完整步驟這里給出一個(gè)標(biāo)準(zhǔn)化流程基本適用于大多數(shù)表格型數(shù)據(jù)。第一步是整理數(shù)據(jù)保證每一行是一個(gè)樣本每一列是一個(gè)變量變量之間必須是連續(xù)數(shù)值。第二步是估計(jì)均值和協(xié)方差通常用mean(X)和cov(X)。第三步是計(jì)算每個(gè)樣本的馬氏距離平方。第四步是確定閾值推薦使用卡方分布的上分位點(diǎn)。第五步是把距離超過閾值的樣本標(biāo)記為異常然后剔除或替換。3.1.1 數(shù)據(jù)形狀要求如果樣本數(shù)n小于等于變量數(shù)pcov(X)是奇異矩陣馬氏距離直接失效。這種情況下需要先降維或者用正則化協(xié)方差估計(jì)。對(duì)很多高維場景比如基因表達(dá)譜或高光譜數(shù)據(jù)直接用mahal是行不通的。一個(gè)常見做法是先用 PCA 把維度壓到主成分個(gè)數(shù)小于樣本數(shù)再對(duì)主成分分?jǐn)?shù)計(jì)算馬氏距離。但需要注意 PCA 本身對(duì)異常值敏感異常樣本會(huì)影響主成分方向。3.1.2 估計(jì)協(xié)方差矩陣的細(xì)節(jié)協(xié)方差矩陣的估計(jì)方法直接影響判別效果。普通cov使用簡單算術(shù)平均如果樣本中存在離群點(diǎn)這些點(diǎn)會(huì)“拉大”協(xié)方差結(jié)果可能讓真正的大偏差看起來不極端這就是所謂的掩蔽效應(yīng)。解決思路是使用穩(wěn)健協(xié)方差估計(jì)比如 MCD最小協(xié)方差行列式。MATLAB 里robustcov函數(shù)就是基于 MCD后面會(huì)有例子。3.2 可運(yùn)行的 MATLAB 函數(shù)下面這個(gè)函數(shù)可以直接復(fù)制保存為removeOutliersByMahal.m輸入數(shù)據(jù)X和顯著性水平alpha輸出剔除后的矩陣和異常索引。function [X_clean, outlierIdx] removeOutliersByMahal(X, alpha) % 輸入 % X - n x p 數(shù)據(jù)矩陣np2 % alpha - 顯著性水平默認(rèn) 0.05 % 輸出 % X_clean - 剔除異常后的數(shù)據(jù) % outlierIdx - 異常樣本的行索引 if nargin 2 || isempty(alpha) alpha 0.05; end % 檢查數(shù)據(jù)形狀 [n, p] size(X); if n p error(樣本數(shù)必須大于變量數(shù)當(dāng)前 cov 矩陣奇異); end % 計(jì)算馬氏距離平方 D2 mahal(X, X); % 卡方分布閾值自由度等于變量數(shù) p threshold chi2inv(1 - alpha, p); % 標(biāo)記異常 outlierIdx find(D2 threshold); X_clean X; X_clean(outlierIdx, :) []; % 刪除異常行 end這里mahal(X, X)有一個(gè)細(xì)節(jié)第二參數(shù)X被當(dāng)作參考總體函數(shù)內(nèi)部會(huì)用mean(X)和cov(X)作為均值向量和協(xié)方差矩陣。如果你有一個(gè)干凈的參考樣本集Xref想用它對(duì)新的Xnew打分應(yīng)該寫成mahal(Xnew, Xref)這更符合實(shí)際生產(chǎn)中的“訓(xùn)練/測試分離”思路。閾值使用chi2inv是因?yàn)樵诙嘣龖B(tài)假設(shè)下馬氏距離平方服從自由度為p的卡方分布。如果數(shù)據(jù)明顯不是正態(tài)卡方閾值會(huì)偏保守或偏激進(jìn)這時(shí)可以考慮基于經(jīng)驗(yàn)分布取 99% 分位數(shù)作為閾值。3.3 閾值確定卡方分布與經(jīng)驗(yàn)分位數(shù)的取舍3.3.1 為什么要用卡方分布多元正態(tài)分布有一個(gè)已知結(jié)論樣本到總體中心的馬氏距離平方服從卡方分布自由度是變量數(shù)。因此chi2inv(0.95, p)能給出一個(gè)理論上的 95% 覆蓋范圍。這個(gè)結(jié)論在小樣本時(shí)并不特別精確當(dāng)n在 50 以下尤其是p接近n時(shí)卡方閾值會(huì)低估異常比例導(dǎo)致異常樣本漏檢。這時(shí)候我傾向于用經(jīng)驗(yàn)分布直接取D2的 97.5% 分位數(shù)作為閾值。但對(duì)小樣本極端值會(huì)影響分位數(shù)估計(jì)所以沒有絕對(duì)安全的選擇。3.3.2 chi2inv 的用法chi2inv是統(tǒng)計(jì)工具箱的函數(shù)第一個(gè)參數(shù)是累積概率值第二個(gè)參數(shù)是自由度。比如chi2inv(0.99, 5)返回 5 個(gè)自由度下卡方分布 99% 分位數(shù)。注意顯著性水平alpha與分位數(shù)的關(guān)系閾值取1 - alpha的分位數(shù)所以alpha0.05等同于 95% 覆蓋。工程上常見的alpha是 0.025 或 0.01因?yàn)楫惓Ｖ低巧贁?shù)拒絕域太大會(huì)誤刪正常點(diǎn)。有一個(gè)思路是先用較小的alpha剔除強(qiáng)異常再對(duì)剩余數(shù)據(jù)重新估計(jì)協(xié)方差這是迭代剔除的雛形。4. 閾值與協(xié)方差估計(jì)三個(gè)影響剔除結(jié)果的關(guān)鍵參數(shù)4.1 置信度 alpha0.975 還是 0.99alpha是卡方分布的分位數(shù)不是實(shí)際異常比例。如果你知道數(shù)據(jù)中大約有 5% 的異常就把閾值設(shè)到 95% 分位數(shù)附近如果異常比例很低建議用 99% 分位數(shù)。實(shí)際操作中可以畫一下D2的直方圖看尾部從哪里開始顯著脫離卡方曲線。我自己經(jīng)常在 0.01 和 0.05 之間做敏感性分析如果剔除結(jié)果對(duì)alpha劇烈變化說明數(shù)據(jù)中異常樣本還不是明顯偏離總體需要回到特征工程層面。4.2 樣本量與維度比協(xié)方差矩陣的穩(wěn)定性這是馬氏距離最大的一道坎。當(dāng)n和p的比值小于 2.5 時(shí)cov(X)本身噪聲太大馬氏距離的有效性會(huì)快速下降。比如一個(gè) 40 行 20 列的數(shù)據(jù)集協(xié)方差矩陣需要估計(jì)p(p1)/2個(gè)獨(dú)立參數(shù)也就是 210 個(gè)值但樣本只有 40 個(gè)估計(jì)結(jié)果嚴(yán)重過擬合inv(Sigma)會(huì)把微小噪聲放大成巨大的距離值。應(yīng)對(duì)方式有三種傾向第一種是做特征選擇保留最重要的變量第二種是使用正則化協(xié)方差比如 Ledoit-Wolf 收縮估計(jì)MATLAB 里cov(X)沒有內(nèi)置參數(shù)但可以自己寫收縮公式第三種是改用基于馬氏距離的穩(wěn)健版本也就是robustcov它通過子集抽樣避免協(xié)方差被異常點(diǎn)污染。4.3 穩(wěn)健估計(jì)用 robustcov 解決掩蔽效應(yīng)當(dāng)異常值本身數(shù)量不多但幅度很大時(shí)普通cov估計(jì)出的協(xié)方差遠(yuǎn)大于真實(shí)總體協(xié)方差導(dǎo)致所有點(diǎn)看起來都接近中心馬氏距離失效。robustcov基于 MCD 算法它先尋找一個(gè)子集使得子集樣本的協(xié)方差行列式最小再用這個(gè)子集的均值和協(xié)方差計(jì)算距離。這能有效避免掩蔽效應(yīng)但缺點(diǎn)是計(jì)算量大數(shù)據(jù)量超過幾萬行時(shí)很吃內(nèi)存。下面是一個(gè)對(duì)比示例估計(jì)方式適用場景缺點(diǎn)推薦用途普通 cov數(shù)據(jù)干凈、異常比例低于 1%對(duì)異常敏感可能漏檢快速初篩穩(wěn)健 MCD異常比例 10%-20%且無明顯規(guī)律計(jì)算慢需要統(tǒng)計(jì)工具箱正式建模前的清洗收縮估計(jì)高維小樣本p 接近 n需要選擇收縮強(qiáng)度基因、光譜數(shù)據(jù)代碼上用robustcov替換普通協(xié)方差通常配合計(jì)算穩(wěn)健馬氏距離平方[sigmaRob, muRob, w2, mahDistRob] robustcov(X); % sigmaRob 為穩(wěn)健協(xié)方差muRob 為穩(wěn)健均值向量 % mahDistRob 為穩(wěn)健馬氏距離平方等價(jià)于對(duì) X 中的每行計(jì)算注意robustcov的第三個(gè)輸出w2是每個(gè)樣本的權(quán)重可用作異常程度評(píng)分。權(quán)重大于 0.5 的樣本通常被認(rèn)為是正常點(diǎn)這個(gè)經(jīng)驗(yàn)值在不少工程場景中有效。使用穩(wěn)健估計(jì)后閾值依然可以用卡方分布但自由度仍然是p因?yàn)槔碚摲植紱]有變。5. 實(shí)戰(zhàn)MATLAB 多元數(shù)據(jù)異常檢測腳本與 CSV 接入5.1 準(zhǔn)備模擬數(shù)據(jù)與噪聲注入為了完整演示剔除流程這里生成一個(gè)含相關(guān)性的三維數(shù)據(jù)集并注入少量異常點(diǎn)。實(shí)際使用時(shí)你可以用readtable或readmatrix把 CSV 數(shù)據(jù)導(dǎo)入替換這里的模擬部分。模擬數(shù)據(jù)的關(guān)鍵是讓兩列之間存在線性關(guān)系這樣才能體現(xiàn)馬氏距離相對(duì)于歐氏距離的優(yōu)勢。5.2 完整可運(yùn)行腳本% detectOutliersDemo.m % 生成帶相關(guān)性的三維數(shù)據(jù)注入異常點(diǎn)用馬氏距離剔除 rng(1); % 固定隨機(jī)種子便于復(fù)現(xiàn) n 200; p 3; % 基礎(chǔ)數(shù)據(jù)第一列是標(biāo)準(zhǔn)正態(tài)第二列與第一列相關(guān)第三列獨(dú)立 X randn(n, p); X(:, 2) 0.8 * X(:, 1) 0.6 * randn(n, 1); % 注入 10 個(gè)異常點(diǎn)把前 10 行的值整體偏移 X(1:10, :) X(1:10, :) [4, 3, 2]; % 導(dǎo)入外部 CSV 的接法如果數(shù)據(jù)已存在 % data readmatrix(sensor_data.csv); % X data(:, 1:3); % 計(jì)算馬氏距離平方 D2 mahal(X, X); % 卡方閾值alpha0.02自由度 3 alpha 0.02; threshold chi2inv(1 - alpha, p); % 標(biāo)記異常 outlierIdx find(D2 threshold); % 剔除異常并輸出結(jié)果 X_clean X; X_clean(outlierIdx, :) []; fprintf(總樣本數(shù)%d\n, n); fprintf(檢出異常點(diǎn)數(shù)%d\n, length(outlierIdx)); fprintf(理論閾值%.2f\n, threshold); % 繪圖對(duì)比前兩個(gè)變量散點(diǎn)圖正常點(diǎn)與異常點(diǎn)用不同顏色 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(馬氏距離檢測異常值結(jié)果); grid on;運(yùn)行后你會(huì)看到紅色叉號(hào)集中在數(shù)據(jù)云的邊緣而且被擠向相關(guān)方向上偏出去的區(qū)域而不是單純?nèi)∽兞拷^對(duì)值的極值。這里的alpha0.02意味著預(yù)期誤判率約為 2%。如果注入異常偏離強(qiáng)度更大檢出率會(huì)更高。如果你發(fā)現(xiàn)異常點(diǎn)沒有被正確分離先檢查是否數(shù)據(jù)中存在缺失值mahal遇到NaN會(huì)直接讓整個(gè)協(xié)方差矩陣崩潰。5.3 結(jié)果解釋與參數(shù)調(diào)整5.3.1 觀察距離排序不要只看閾值把D2從大到小排序取前 20 個(gè)索引觀察。如果前 10 個(gè)恰好是注入的異常后 10 個(gè)是正常邊界點(diǎn)說明閾值偏嚴(yán)。這種情況下把a(bǔ)lpha調(diào)到 0.05或者改為取距離排序的后 5% 作為異常都會(huì)改變最終清洗后的數(shù)據(jù)分布。建議在剔除前先保存一份D2變量隨后畫一個(gè)距離分布直方圖和卡方概率密度曲線疊加對(duì)比目視檢查尾部是否一致。5.3.2 導(dǎo)出剔除后的數(shù)據(jù)writematrix可以避免手工復(fù)制writematrix(X_clean, X_clean.csv);注意這里覆蓋了原文件內(nèi)容所以運(yùn)行時(shí)先確認(rèn)路徑。生產(chǎn)環(huán)境中我會(huì)把異常索引存成outlierIdx.csv保留原始數(shù)據(jù)而不直接刪除方便溯源。馬氏距離剔除的局限性在于如果異常是以局部模式出現(xiàn)比如某個(gè)傳感器只在一段時(shí)間內(nèi)失效那么距離本身難以區(qū)分“正常變異性”和“故障偏移”。此時(shí)可以考慮對(duì)時(shí)間序列加滑動(dòng)窗口在每個(gè)窗口內(nèi)計(jì)算局部馬氏距離再對(duì)距離序列做趨勢分析。6. 進(jìn)階穩(wěn)健協(xié)方差估計(jì)與異常值可視化驗(yàn)證當(dāng)數(shù)據(jù)中已經(jīng)混入一批異常值普通mahal的協(xié)方差估計(jì)會(huì)被污染導(dǎo)致距離分?jǐn)?shù)偏低。此時(shí)可以改用robustcov得到穩(wěn)健距離并配合 Q-Q 圖做驗(yàn)證。先看穩(wěn)健版本的核心調(diào)用[sigmaRob, muRob, ~, D2Rob] robustcov(X); thresholdRob chi2inv(0.99, p); outlierRob D2Rob thresholdRob;robustcov的默認(rèn)方法是用 Fast-MCD 算法它會(huì)從樣本中隨機(jī)抽取子集迭代計(jì)算因此結(jié)果帶有隨機(jī)性。建議設(shè)置隨機(jī)種子以重復(fù)實(shí)驗(yàn)或者多次運(yùn)行收集異常索引的并集。robustcov的第四輸出已經(jīng)是穩(wěn)健馬氏距離平方不需要再手動(dòng)減均值乘逆矩陣。驗(yàn)證手段之一是畫卡方 Q-Q 圖把距離平方排序后與卡方分布的分位數(shù)做散點(diǎn)。正常數(shù)據(jù)應(yīng)該大致落在直線附近右上方明顯翹起的點(diǎn)就是異常候選。MATLAB 里沒有直接的卡方 Q-Q 圖函數(shù)可以這樣生成% 生成理論分位數(shù) p_seq (1:n) / (n 1); theoretical chi2inv(p_seq, p); % 對(duì) D2 排序后畫散點(diǎn) D2_sorted sort(D2); plot(theoretical, D2_sorted, o); hold on; plot(theoretical, theoretical, k--); % 參考對(duì)角線 xlabel(卡方理論分位數(shù)); ylabel(馬氏距離平方排序后);第二點(diǎn)經(jīng)驗(yàn)是固定異常比例?？ǚ介撝颠m合正態(tài)數(shù)據(jù)但工程數(shù)據(jù)總帶偏態(tài)我習(xí)慣用prctile取 95% 分位數(shù)作為閾值這樣不用反復(fù)調(diào)alpha。但注意這種方法輸出的異常數(shù)量和比例是預(yù)設(shè)的可能錯(cuò)把邊界點(diǎn)圈進(jìn)來。更保險(xiǎn)的做法是對(duì)D2做對(duì)數(shù)變換再對(duì)變換后的數(shù)據(jù)用 3σ 法則因?yàn)閷?duì)數(shù)變換后的極端值更接近對(duì)稱分布。最后一個(gè)實(shí)操細(xì)節(jié)mahal和robustcov都會(huì)因變量單位不同而得到相同的距離因?yàn)閰f(xié)方差矩陣吸收了尺度信息。但這不意味著數(shù)據(jù)不需要預(yù)處理。當(dāng)某個(gè)變量的方差極小比如接近機(jī)器精度時(shí)協(xié)方差矩陣中對(duì)應(yīng)行列接近零逆矩陣放大該維度上的微小偏差本來正常的測量噪聲會(huì)被誤判為異常。處理方法是先剔除方差接近零的變量或者用zscore標(biāo)準(zhǔn)化后再計(jì)算馬氏距離。兩種做法會(huì)得到幾乎一樣的結(jié)果但標(biāo)準(zhǔn)化后的協(xié)方差矩陣數(shù)值上更穩(wěn)定也能避免mahal因?yàn)榫仃嚥B(tài)返回Inf。本文還有配套的精品資源點(diǎn)擊獲取