在（下）篇中，我们将和大家一起探讨无监督异常检测。参考文档页面：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()。
【定义】马氏距离：一种衡量样本和数据集分布间相似度的尺度无关的度量指标，例如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)nexttilegscatter(X(:,1),X(:,2),tf_robustcov,[],'ox',3)xlabel('X1')ylabel('X2')legend({'正常','异常'})
title("训练样本分类结果 (tSNE降维)")
nexttileplot(d,mah,'o')line([mah_threshold, mah_threshold], [0, 30], 'color', 'r')line([0, 6], [mah_threshold, mah_threshold], 'color', 'r')hold onplot(d(tf), mah(tf), 'r+')xlabel('Mahalanobis Distance')ylabel('Robust Distance')title('DD Plot')hold off
nexttileconfusionchart(trueAnomaliesTest, predTest, ... Title="测试结果评估-混淆矩阵 (马氏距离)", Normalization="row-normalized");
· 以上方法在高维数据上应用效果不理想，可以看到，测试集中异常数据假阴率高。这个方法适用于数据符合或接近正态分布的情况，但是通常情况下，实际数据的分布规律难以预估。
局部离群因子
参考文档页面：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)nexttilegscatter(X(:,1),X(:,2),tfLOF,[],'ox',3)xlabel('X1')ylabel('X2')legend({'正常','异常'})title("训练样本分类结果 (tSNE降维)")
nexttileidxes = randi(NumSamples,1,60);scatter(X(idxes,1),X(idxes,2),5,'filled','MarkerFaceColor','k')hold onbubblechart(X(idxes,1),X(idxes,2),scoresLOF(idxes)/100, ...'r','MarkerFaceAlpha',0);legend({'数据点','异常得分'})hold offtitle("训练样本异常得分分布")
nexttileconfusionchart(trueAnomaliesTest, predLOF, ... Title="测试结果评估-混淆矩阵 (LOF)", Normalization="row-normalized");
· 优点：不受数据分布的影响，同时考虑了数据集的局部和全局属性，比较适用于中等高维的数据集，针对示例数据集的预测准确度比较理想。
· 使用限制：对近邻参数较为敏感，由于需要计算数据集中任意两个数据点的距离，算法的时间复杂度较高，在大规模数据集上效率偏低，适合小规模到中等规模的数值型数据。
孤立森林
参考文档页面：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)nexttilegscatter(X(:,1),X(:,2),tfIF,[],'ox',3)xlabel('X1')ylabel('X2')legend({'正常','异常'})title("训练样本分类结果 (tSNE降维)")
nexttilehistogram(scoreTrainIF)xline(mdlIF.ScoreThreshold,"k-",join(["Threshold =" mdlIF.ScoreThreshold]))title("训练样本异常得分分布")
nexttileconfusionchart(trueAnomaliesTest, predIF, ... Title="测试结果评估-混淆矩阵 (孤立森林)", Normalization="row-normalized");
· 优点：适合高维表格数据，不需要计算关于距离和密度的指标，具有线性时间复杂度，每个决策树可独立采样，支持并行化处理来实现加速。
· 使用限制：孤立森林适用于训练集和测试集中，正常样本和异常样本占比接近的情况，且异常样本的特征与正常样本差异很大。
单类支持向量机
参考文档页面：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)nexttilegscatter(X(:,1),X(:,2),tfOCSVM,[],'ox',3)xlabel('X1')ylabel('X2')legend({'正常','异常'})title("训练样本分类结果 (tSNE降维)")
nexttilehistogram(scoreTrainOCSVM, Normalization="probability")xline(mdlSVM.ScoreThreshold,"k-", ... join(["Threshold =" mdlSVM.ScoreThreshold]))title("训练样本异常得分分布")
nexttileconfusionchart(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
综上，各个方法的适用范围不一，或是有特定的使用条件，在使用时需要多加留意，例如马氏距离适合符合正态分布假设的数据集，孤立森林适用于处理正常样本和异常样本差异较大的情况，各个算法计算复杂度有些许区别，可以根据实际情况选择合适的方法。
· 在处理异常检测问题时，首先需要充分了解您的数据
· 如果您有足够的标注数据（包括异常），可使用有监督学习方法进行异常检测
· 如果您的数据大部分都是正常数据，或者异常数据难以获取或标记，则可以考虑使用无监督的异常检测方法
如果您有更多感兴趣的话题或其他观点，欢迎留言告知我们。
[1] Unsupervised Anomaly Detection:https://ww2.mathworks.cn/help/stats/unsupervised-anomaly-detection.html%20/l%20mw_770db8b2-dea5-4ab2-8ab2-287957d26f90
[2]robustcov:https://ww2.mathworks.cn/help/stats/robustcov.html
[3]lof:https://127.0.0.1:31515/static/help/stats/lof.html\
[4]iforest:https://127.0.0.1:31515/static/help/stats/iforest.html
[5]ocsvm:https://127.0.0.1:31515/static/help/stats/ocsvm.html
[6]Local outlier factor model for anomaly detection:https://ww2.mathworks.cn/help/releases/R2022b/stats/localoutlierfactor.html?requestedDomain=cn
[7]Anomaly Detection with Isolation Forest:https://ww2.mathworks.cn/help/releases/R2022b/stats/anomaly-detection-with-isolation-forest.html?requestedDomain=cn
[8]Fit one-class support vector machine (SVM) model for anomaly detection:https://ww2.mathworks.cn/help/releases/R2022b/stats/ocsvm.html
