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

沒有標簽怎么辦？試試無監督異常檢測

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

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

· 馬氏距離 (Mahalanobis Distance)：

如果數據符合多變量正態分布，可使用樣本到數據集分布中心的馬氏距離檢測異常。利用穩健協方差估計robustcov[2]函數計算馬氏距離，進行離群值檢測。

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

基于觀測點和鄰近樣本之間的相對密度檢測異常點。利用 lof[3] 函數可以創建 LocalOutlierFactor 對象，針對訓練數據，直接返回檢測結果。

· 孤立森林 (Isolation forest)：

通過一組孤立樹模型，將異常值和正常數據點隔離。利用 iforest[4] 函數，創建 IsolationForest 對象，針對訓練數據，直接返回檢測結果。

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

在無監督條件下，訓練支持向量機模型，在變換后的高維空間，將數據點和原點分離。利用 ocsvm[5] 函數創建OneClassSVM 對象，針對訓練數據，直接返回檢測結果。

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

具體實現方式

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

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

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

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

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

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

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

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

利用pdist2函數，計算測試集的馬氏距離后與閾值mah_threshold進行比較：

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

可視化檢測結果：

tiledlayout(3,1)
nexttile
gscatter(X(:,1),X(:,2),tf_robustcov,[],'ox',3)
xlabel('X1')
ylabel('X2')
legend({'正常','異常'})
title("訓練樣本分類結果 (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="測試結果評估-混淆矩陣 (馬氏距離)", Normalization="row-normalized");

? ? ? ?

小結

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

2.局部離群因子

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

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

其中，為k近鄰集合，樣本p關于觀測點o的局部可達距離定義為：

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

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

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

因此，對于正常樣本，一般 LOF 值小于或接近 1，意味著其局部可達密度和近鄰點相近或更高，該樣本和鄰域內的樣本同屬一個簇，當 LOF 值大于 1 時，則可能為異常值，利用 ContaminationFraction 參數可調整 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"]);

可視化檢測結果：

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

隨機選取部分樣本，查看對應 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({'數據點','異常得分'})
hold off
title("訓練樣本異常得分分布")


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

? ? ? ?

小結

·優點：不受數據分布的影響，同時考慮了數據集的局部和全局屬性，比較適用于中等高維的數據集，針對示例數據集的預測準確度比較理想。

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

3.孤立森林

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

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

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

可視化檢測結果：

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


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


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

? ? ? ?

小結

·優點：適合高維表格數據，不需要計算關于距離和密度的指標，具有線性時間復雜度，每個決策樹可獨立采樣，支持并行化處理來實現加速。

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

4.單類支持向量機

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

單類支持向量機，或無監督支持向量機，構建決策邊界，將訓練集中的數據點盡可能劃分為一個類別，位于決策邊界之外的數據則為異常值。

策略是通過核函數將數據映射到新的高維特征空間，在數據與原點間構建超平面（n 維平面）。因為，在低維空間中的非線性特征往往不是線性可分的，在擴展到高維空間后是可分的。在 MATLAB 中，一種實現方法是使用用于構建標準的支持向量機分類模型的 fitcsvm 函數（MATLAB R2022b前），另一種實現方法是使用 ocsvm 函數（MATLAB R2022b 起）。

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

>> openExample('stats/DetectOutliersUsingSVMAndOneClassLearningExample')

ocsvm 函數的求解則是基于 SVM 的原型問題形式，并使用高斯核進行一類學習，以找到決策邊界，針對大規模數據集，求解效率更高。將 KernelScale 設為 "auto"以啟發式地選取合適核函數參數。

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

可視化檢測結果：

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


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


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

? ? ? ?

小結

·優點：可以處理高維數據，適合表格/結構化數據。

·使用限制：分類變量需要先轉換為虛擬變量（啞變量，Dummy Variable），ocsvm函數中可定義相關參數CategoricalPredictors ，以自動進行轉換。

結果對比

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

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

ans = 0.9325

利用馬氏距離/穩健協方差估計的方法，結果不理想，與其他三個算法的結果差異較大：

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

ans = 0.6196

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

關鍵點回顧

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

·如果您有足夠的標注數據（包括異常），可使用有監督學習方法進行異常檢測

·如果您的數據大部分都是正常數據，或者異常數據難以獲取或標記，則可以考慮使用無監督的異常檢測方法

審核編輯：湯梓紅