在使用 MATLAB 進(jìn)行異常檢測（上）中，我們探討了什么是異常值，簡單的一維數(shù)據(jù)異常檢測問題，針對高維數(shù)據(jù)的有監(jiān)督異常檢測方法。 在（下）篇中，我們將和大家一起探討無監(jiān)督異常檢測。

沒有標(biāo)簽怎么辦？試試無監(jiān)督異常檢測

參考文檔頁面：Unsupervised Anomaly Detection[1]

對于沒有標(biāo)簽信息的多變量樣本數(shù)據(jù)，在MATLAB 可以使用以下方法檢測異常值:

· 馬氏距離 (Mahalanobis Distance)：

如果數(shù)據(jù)符合多變量正態(tài)分布，可使用樣本到數(shù)據(jù)集分布中心的馬氏距離檢測異常。利用穩(wěn)健協(xié)方差估計(jì)robustcov[2]函數(shù)計(jì)算馬氏距離，進(jìn)行離群值檢測。

· 局部離群因子 (Local Outlier Factor, LOF)：

基于觀測點(diǎn)和鄰近樣本之間的相對密度檢測異常點(diǎn)。利用 lof[3] 函數(shù)可以創(chuàng)建 LocalOutlierFactor 對象，針對訓(xùn)練數(shù)據(jù)，直接返回檢測結(jié)果。

· 孤立森林 (Isolation forest)：

通過一組孤立樹模型，將異常值和正常數(shù)據(jù)點(diǎn)隔離。利用 iforest[4] 函數(shù)，創(chuàng)建 IsolationForest 對象，針對訓(xùn)練數(shù)據(jù)，直接返回檢測結(jié)果。

· 單類支持向量機(jī) (One Class SVM)：

在無監(jiān)督條件下，訓(xùn)練支持向量機(jī)模型，在變換后的高維空間，將數(shù)據(jù)點(diǎn)和原點(diǎn)分離。利用 ocsvm[5] 函數(shù)創(chuàng)建OneClassSVM 對象，針對訓(xùn)練數(shù)據(jù)，直接返回檢測結(jié)果。

針對測試數(shù)據(jù)進(jìn)行異常檢測時，使用第一種方法的具體檢測步驟將在下文中通過示例說明，如果使用另外三種算法，可以直接調(diào)用檢測模型的對象函數(shù) isanomaly()。

具體實(shí)現(xiàn)方式

1.利用馬氏距離檢測異常值

【定義】馬氏距離：一種衡量樣本和數(shù)據(jù)集分布間相似度的尺度無關(guān)的度量指標(biāo)，例如y到中值點(diǎn)的距離是d2=(y-μ)Σ-1(y-μ)'，其中Σ是多維隨機(jī)變量的協(xié)方差矩陣，μ為樣本均值。在MATLAB中可通過pdist2函數(shù)計(jì)算：

d = pdist2(feat,mean(feat),"mahalanobis");

馬氏距離可以理解為是對歐式距離的一種修正，假設(shè)，下圖中藍(lán)色點(diǎn)和黃色點(diǎn)離樣本均值的歐式距離相近，但是由于樣本整體分布沿 f(x)=x 的方向分布（變量之間具有相關(guān)性），藍(lán)色點(diǎn)更有可能是數(shù)據(jù)集中的點(diǎn)，對應(yīng)的馬氏距離更小，而黃色點(diǎn)更有可能是離群值，對應(yīng)馬氏距離也更大。因此，設(shè)定一個合理的閾值，可以劃分異常樣本和正常樣本。

計(jì)算馬氏距離，首先需要估計(jì) Σ(sig) 和 μ(mu) 。極大似然估計(jì)（Maximum Likelihood Estimation, MLE）對數(shù)據(jù)中的異常值非常敏感，需要采取一種穩(wěn)健的協(xié)方差估計(jì)方法，抵抗數(shù)據(jù)集中存在的異常觀測數(shù)據(jù)。當(dāng)數(shù)據(jù)中存在異常值時，協(xié)方差行列式偏大。使用最小協(xié)方差行列式估計(jì) (Minimum Covariance Determinant, MCD)，從 n 個數(shù)據(jù)樣本中，最多選取h個觀測值，找到協(xié)方差行列式最小的一組觀測子集，計(jì)算其平均值和協(xié)方差，作為估計(jì)量。robustcov 函數(shù)提供了 FAST-MCD、OGK 和 Olive-Hawkins 三種算法供選擇。假設(shè)訓(xùn)練數(shù)據(jù)集中，異常占比為 0.9%：

contaminationFraction = 0.09;
[sig,mu,mah,tf] = robustcov(feat, ...
  OutlierFraction=contaminationFraction);

如果數(shù)據(jù)符合正態(tài)分布的假設(shè)，馬氏距離的平方值，將服從具有 Dim 個自由度的χ2分布，Dim 為原數(shù)據(jù)的維度。默認(rèn)情況下，robustcov 函數(shù)假設(shè)數(shù)據(jù)符合多變量正態(tài)分布，并根據(jù)χ2分布的臨界值，將輸入樣本的 2.5% 作為異常值，如需調(diào)整異常值占比，可以使用 chi2inv 函數(shù)，重新計(jì)算閾值：

mah_threshold = sqrt(chi2inv(1-contaminationFraction,Dim));
tf_robustcov = mah > mah_threshold;

利用pdist2函數(shù)，計(jì)算測試集的馬氏距離后與閾值mah_threshold進(jìn)行比較：

dTest = pdist2(featureTestNoLabels.Variables, ...
mean(featureTestNoLabels.Variables),"mahalanobis");
isanomalyTest = dTest > mah_threshold;
predTest = categorical(isanomalyTest, [1, 0], ["Anomaly", "Normal"]);

可視化檢測結(jié)果：

tiledlayout(3,1)
nexttile
gscatter(X(:,1),X(:,2),tf_robustcov,[],'ox',3)
xlabel('X1')
ylabel('X2')
legend({'正常','異常'})
title("訓(xùn)練樣本分類結(jié)果 (tSNE降維)")


nexttile
plot(d,mah,'o')
line([mah_threshold, mah_threshold], [0, 30], 'color', 'r')
line([0, 6], [mah_threshold, mah_threshold], 'color', 'r')
hold on
plot(d(tf), mah(tf), 'r+')
xlabel('Mahalanobis Distance')
ylabel('Robust Distance')
title('DD Plot')
hold off


nexttile
confusionchart(trueAnomaliesTest, predTest, ...
  Title="測試結(jié)果評估-混淆矩陣 (馬氏距離)", Normalization="row-normalized");

小結(jié)

·以上方法在高維數(shù)據(jù)上應(yīng)用效果不理想，可以看到，測試集中異常數(shù)據(jù)假陰率高。這個方法適用于數(shù)據(jù)符合或接近正態(tài)分布的情況，但是通常情況下，實(shí)際數(shù)據(jù)的分布規(guī)律難以預(yù)估。

2.局部離群因子

參考文檔頁面：Local outlier factor model for anomaly detection[6]

該算法通過計(jì)算樣本p和其周圍 k 個近鄰點(diǎn)的局部可達(dá)密度(local reachability density, lrd)，即觀測樣本 p 到近鄰點(diǎn)的局部可達(dá)距離平均值的倒數(shù)：

其中，為k近鄰集合，樣本p關(guān)于觀測點(diǎn)o的局部可達(dá)距離定義為：

其中 dk(0)為觀測點(diǎn)到其近鄰的第k個最小距離, d(p,0) 為樣本 p 和觀測點(diǎn) o 之間的距離，可參考下圖示意。

再根據(jù)p的局部可達(dá)密度與近鄰點(diǎn)的局部可達(dá)密度比值的平均值，量化每個樣本的離群程度，具體計(jì)算可參考以下公式：

Ird(·)為局部可達(dá)密度函數(shù)，|Nk(p)|為近鄰數(shù)量。

因此，對于正常樣本，一般 LOF 值小于或接近 1，意味著其局部可達(dá)密度和近鄰點(diǎn)相近或更高，該樣本和鄰域內(nèi)的樣本同屬一個簇，當(dāng) LOF 值大于 1 時，則可能為異常值，利用 ContaminationFraction 參數(shù)可調(diào)整 LOF 的閾值。

[mdlLOF,tfLOF,scoresLOF] = lof(feat, ...
  ContaminationFraction=0.09, ...
  NumNeighbors=1000, Distance="mahalanobis");
[isanomalyLOF,~] = isanomaly(mdlLOF, featureTestNoLabels.Variables);
predLOF = categorical(isanomalyLOF, [1, 0], ["Anomaly", "Normal"]);

可視化檢測結(jié)果：

tiledlayout(3,1)
nexttile
gscatter(X(:,1),X(:,2),tfLOF,[],'ox',3)
xlabel('X1')
ylabel('X2')
legend({'正常','異常'})
title("訓(xùn)練樣本分類結(jié)果 (tSNE降維)")

隨機(jī)選取部分樣本，查看對應(yīng) LOF 值/異常得分

nexttile
idxes = randi(NumSamples,1,60);
scatter(X(idxes,1),X(idxes,2),5,'filled','MarkerFaceColor','k')
hold on
bubblechart(X(idxes,1),X(idxes,2),scoresLOF(idxes)/100, ...
'r','MarkerFaceAlpha',0);
legend({'數(shù)據(jù)點(diǎn)','異常得分'})
hold off
title("訓(xùn)練樣本異常得分分布")


nexttile
confusionchart(trueAnomaliesTest, predLOF, ...
  Title="測試結(jié)果評估-混淆矩陣 (LOF)", Normalization="row-normalized");

小結(jié)

·優(yōu)點(diǎn)：不受數(shù)據(jù)分布的影響，同時考慮了數(shù)據(jù)集的局部和全局屬性，比較適用于中等高維的數(shù)據(jù)集，針對示例數(shù)據(jù)集的預(yù)測準(zhǔn)確度比較理想。

·使用限制：對近鄰參數(shù)較為敏感，由于需要計(jì)算數(shù)據(jù)集中任意兩個數(shù)據(jù)點(diǎn)的距離，算法的時間復(fù)雜度較高，在大規(guī)模數(shù)據(jù)集上效率偏低，適合小規(guī)模到中等規(guī)模的數(shù)值型數(shù)據(jù)。

3.孤立森林

參考文檔頁面：Anomaly Detection with Isolation Forest[7]

孤立森林算法中，集成了多個決策樹模型，訓(xùn)練時，每個決策樹對一個不放回采樣的數(shù)據(jù)子集進(jìn)行分裂，以試圖將每一個觀測樣本劃分到一個對應(yīng)的葉節(jié)點(diǎn)上。假設(shè)異常點(diǎn)與其他正常數(shù)據(jù)差異較大，從根節(jié)點(diǎn)到對應(yīng)葉節(jié)點(diǎn)需要經(jīng)過的路徑長度(path length) 相對較短，對于每個樣本，將孤立森林中的多個決策樹路徑長度的平均值，定義為對應(yīng)樣本的異常得分(anomaly score)。

[mdlIF,tfIF,scoreTrainIF] = iforest(feat, ContaminationFraction=0.09);
[isanomalyIF,~] = isanomaly(mdlIF, featureTestNoLabels.Variables);
predIF = categorical(isanomalyIF, [1, 0], ["Anomaly", "Normal"]);

可視化檢測結(jié)果：

tiledlayout(3,1)
nexttile
gscatter(X(:,1),X(:,2),tfIF,[],'ox',3)
xlabel('X1')
ylabel('X2')
legend({'正常','異常'})
title("訓(xùn)練樣本分類結(jié)果 (tSNE降維)")


nexttile
histogram(scoreTrainIF)
xline(mdlIF.ScoreThreshold,"k-",join(["Threshold =" mdlIF.ScoreThreshold]))
title("訓(xùn)練樣本異常得分分布")


nexttile
confusionchart(trueAnomaliesTest, predIF, ...
  Title="測試結(jié)果評估-混淆矩陣 (孤立森林)", Normalization="row-normalized");

小結(jié)

·優(yōu)點(diǎn)：適合高維表格數(shù)據(jù)，不需要計(jì)算關(guān)于距離和密度的指標(biāo)，具有線性時間復(fù)雜度，每個決策樹可獨(dú)立采樣，支持并行化處理來實(shí)現(xiàn)加速。

·使用限制：孤立森林適用于訓(xùn)練集和測試集中，正常樣本和異常樣本占比接近的情況，且異常樣本的特征與正常樣本差異很大。

4.單類支持向量機(jī)

參考文檔頁面：Fit one-class support vector machine (SVM) model for anomaly detection[8]

單類支持向量機(jī)，或無監(jiān)督支持向量機(jī)，構(gòu)建決策邊界，將訓(xùn)練集中的數(shù)據(jù)點(diǎn)盡可能劃分為一個類別，位于決策邊界之外的數(shù)據(jù)則為異常值。

策略是通過核函數(shù)將數(shù)據(jù)映射到新的高維特征空間，在數(shù)據(jù)與原點(diǎn)間構(gòu)建超平面（n 維平面）。因?yàn)?，在低維空間中的非線性特征往往不是線性可分的，在擴(kuò)展到高維空間后是可分的。在 MATLAB 中，一種實(shí)現(xiàn)方法是使用用于構(gòu)建標(biāo)準(zhǔn)的支持向量機(jī)分類模型的 fitcsvm 函數(shù)（MATLAB R2022b前），另一種實(shí)現(xiàn)方法是使用 ocsvm 函數(shù)（MATLAB R2022b 起）。

fitcsvm 函數(shù)的求解是基于 SVM 的對偶問題形式，需要求解每對樣本的格拉姆矩陣（Gram Matrix），相關(guān)示例可在 MATLAB 命令行輸入以下指令打開：

>> openExample('stats/DetectOutliersUsingSVMAndOneClassLearningExample')

ocsvm 函數(shù)的求解則是基于 SVM 的原型問題形式，并使用高斯核進(jìn)行一類學(xué)習(xí)，以找到?jīng)Q策邊界，針對大規(guī)模數(shù)據(jù)集，求解效率更高。將 KernelScale 設(shè)為 "auto"以啟發(fā)式地選取合適核函數(shù)參數(shù)。

[mdlSVM,tfOCSVM,scoreTrainOCSVM] = ocsvm(feat, ...
    ContaminationFraction=0.09, ...
    StandardizeData=true,KernelScale="auto");
[isanomalyOCSVM,~] = isanomaly(mdlSVM, featureTestNoLabels.Variables);
predOCSVM = categorical(isanomalyOCSVM, [1, 0], ["Anomaly", "Normal"]);

可視化檢測結(jié)果：

tiledlayout(3,1)
nexttile
gscatter(X(:,1),X(:,2),tfOCSVM,[],'ox',3)
xlabel('X1')
ylabel('X2')
legend({'正常','異常'})
title("訓(xùn)練樣本分類結(jié)果 (tSNE降維)")


nexttile
histogram(scoreTrainOCSVM, Normalization="probability")
xline(mdlSVM.ScoreThreshold,"k-", ...
  join(["Threshold =" mdlSVM.ScoreThreshold]))
title("訓(xùn)練樣本異常得分分布")


nexttile
confusionchart(trueAnomaliesTest, predOCSVM,...
  Title="測試結(jié)果評估-混淆矩陣 (OCSVM)", Normalization="row-normalized");

小結(jié)

·優(yōu)點(diǎn)：可以處理高維數(shù)據(jù)，適合表格/結(jié)構(gòu)化數(shù)據(jù)。

·使用限制：分類變量需要先轉(zhuǎn)換為虛擬變量（啞變量，Dummy Variable），ocsvm函數(shù)中可定義相關(guān)參數(shù)CategoricalPredictors ，以自動進(jìn)行轉(zhuǎn)換。

結(jié)果對比

在該數(shù)據(jù)集的離群值檢測問題中，孤立森林、局部離群因子與單類支持向量機(jī)的結(jié)果比較接近，各自的準(zhǔn)確度都比較理想，預(yù)測結(jié)果的重合度也超過了90%：

mean((predIF==predLOF) & (predLOF==predOCSVM))

ans = 0.9325

利用馬氏距離/穩(wěn)健協(xié)方差估計(jì)的方法，結(jié)果不理想，與其他三個算法的結(jié)果差異較大：

mean((predIF==predLOF) & (predLOF==predOCSVM) & (predOCSVM==predTest))

ans = 0.6196

綜上，各個方法的適用范圍不一，或是有特定的使用條件，在使用時需要多加留意，例如馬氏距離適合符合正態(tài)分布假設(shè)的數(shù)據(jù)集，孤立森林適用于處理正常樣本和異常樣本差異較大的情況，各個算法計(jì)算復(fù)雜度有些許區(qū)別，可以根據(jù)實(shí)際情況選擇合適的方法。

關(guān)鍵點(diǎn)回顧

·在處理異常檢測問題時，首先需要充分了解您的數(shù)據(jù)

·如果您有足夠的標(biāo)注數(shù)據(jù)（包括異常），可使用有監(jiān)督學(xué)習(xí)方法進(jìn)行異常檢測

·如果您的數(shù)據(jù)大部分都是正常數(shù)據(jù)，或者異常數(shù)據(jù)難以獲取或標(biāo)記，則可以考慮使用無監(jiān)督的異常檢測方法

審核編輯：湯梓紅