【异常检测】基于主成分分类器的异常检测方案（文献学习）

A novel anomaly detection scheme based on principal component classifier

Mei-Ling Shyu , Shu-Ching Chen, Kanoksri Sarinnapakorn, LiWu Chang
3rd IEEE International Conference on Data Mining, 353–365.(2003)

一种新的基于主成分分类器的异常检测方法

摘要

提出了一种在训练数据不受监督的情况下，利用鲁棒主成分分类器进行入侵检测的新方案。假设异常可以被视为异常值，从正常情况的主成分和次成分构建入侵预测模型。异常与正常情况之间的度量是主成分空间中的距离。基于占总变化50%的主分量和特征值小于0.20的次分量的距离显示，效果良好。对KDD-Cup-1999数据的实验表明，该方法的查全率达到98.94%，查准率达到97.89%，虚警率为0.92%，优于最近邻方法、基于密度的局部离群点（LOF）方法和基于堪培拉度量的离群点检测算法。
关键词： 异常检测、数据挖掘、入侵检测、异常值、主成分分析、距离度量。

1. 介绍

通信网络使物理距离变得毫无意义。人们可以通过网络进行交流，而不受实际距离的限制。在我们享受连接的便利的同时，人们也认识到，恶意或未经授权的用户从一个地方入侵，可能会对广大地区造成严重损害。计算机网络的安全已成为一个关键问题，因此开发防御入侵的机制显得尤为重要。
Heady等人将入侵定义为“试图破坏信息资源的完整性、机密性或可用性的任何一组行为”。对这组恶意行为的识别称为入侵检测。入侵检测问题已经引起了研究者的极大兴趣。1999年，第三届国际知识发现和数据挖掘工具竞赛与第五届知识发现和数据挖掘国际会议（KDD-99）同时举行。竞赛的任务是建立一个网络入侵检测器从KDD杯1999年的数据，这是一个预测模型，能够区分“坏”连接（称为攻击）和“好”正常连接。
现有的入侵检测方法主要分为两大类：特征识别和异常检测。对于特征识别技术，存储已知攻击的特征码，并根据特征码匹配监控事件。这些技术在有匹配时发出入侵信号。这些技术的一个明显的局限性是它们无法检测到特征码未知的新攻击。另一方面，异常检测则建立正常数据模型，检测观测数据中与正常模型的任何偏差。给定一组要训练的正常数据，并给定一个新的测试数据，目标是确定测试数据是属于“正常”还是属于异常行为。异常检测算法的优点是，它们可以检测出与正常使用不同的新类型入侵。然而，缺点是虚警率高。这是因为以前未出现过但合法的系统行为也可能被认为是异常现象。
有各种入侵检测技术属于异常检测范畴，其中一些是机器学习技术（例如，鲁棒支持向量机），一些是基于统计的。
基于统计的异常检测技术利用正常活动的统计特性来建立一个范数分布，并利用统计检验来确定观测到的活动是否显著偏离范数分布。他们通常假设正态或多元正态分布，这可能是一个缺点。文献中提出了一种基于卡方统计量的技术，该技术建立了信息系统中正常事件的范数分布，并检测到与范数分布的较大偏离是一种可能的入侵。实验证明，该方法具有较低的虚警率和较高的检测率。
Emran和Ye开发了一种基于多元统计的异常检测技术Canberra技术，并将其应用于入侵检测。该方法不受数据正态性假设的影响。然而，他们的实验表明，这种技术只有在所有攻击都放在一起的情况下才能表现得非常好。文献提出了一种检测入侵的多元质量控制技术，该技术建立了信息系统中正常活动的长期轮廓，并使用该轮廓来检测异常。这项技术是基于Hotelling的T检验，检测反关系异常和均值漂移异常。当使用一小组计算机审计数据进行测试时，它的性能导致检测到的所有入侵没有假警报；而对于一个大数据集，92%的入侵被检测到，也没有假警报。
有许多异常检测技术采用离群点检测的概念。Aggarwal和Yu讨论了一种离群点检测技术，该技术通过研究数据集中投影的行为来发现离群点。该方法适用于高维数据集。Breunig等人为每个对象指定了一个异常值的程度，称为对象的局部异常因子（ Local Outlier Factor, LOF）。程度取决于对象相对于周围邻近的孤立程度。Lazarevic等人提出了几种用于检测网络入侵的异常检测方案。对这些方案在DARPA 1998网络连接数据集上的比较研究表明，最有前途的技术是LOF方法。
本文提出了一种基于主成分和离群点检测的异常检测方法。该方法的突出假设是攻击以异常值的形式出现在正常数据中。主成分分类器（Principal Component Classifier, PCC）作为一种孤立点检测方法，在故障检测、传感器检测、统计过程控制、分布式传感器网络等领域有着广泛的应用。首先，它没有任何分布假设。许多基于统计的入侵检测方法采用正态分布，或者使用中心极限定理，要求特征数大于30。其次，这类问题的数据具有典型的高维性。因此，在我们提出的方案中，稳健主成分分析（Principal Component Analysis, PCA）被用来降低维数，得到一个简单的分类器，它是主成分的一个简单函数。主成分分析在不牺牲有价值信息的情况下降低了数据的维数。只有少数的主成分参数需要保留，以备将来检测。该方法的另一个优点是，在检测阶段可以用较少的时间计算统计信息，这使得实时使用该方法成为可能。实验结果表明，该方案具有良好的检测率和较低的虚警率。它优于k近邻法、LOF法和堪培拉度量法。

2. 多元统计分析

2.1 距离

设$\mathbf{x}=(x_1,x_2,...,x_p{)}^T$且$\mathbf{y}=(y_1,y_2,...,y_p{)}^T$（p为成分个数，即维数）
欧几里得距离（直线距离）
$d(\mathbf{x},\mathbf{y})=\sqrt{{(\mathbf{x}-\mathbf{y})}^T(\mathbf{x}-\mathbf{y})}$
由于每个特征对欧几里德距离的计算贡献相等，因此在许多应用中不需要这种距离。当特征具有非常不同的可变性或在不同的尺度上测量不同的特征时，具有大尺度测量或高可变性的特征的影响将支配具有较小尺度或较少可变性的其他特征。
马氏距离（Mahalanobis distance）
$d^2(\mathbf{x},\mathbf{y})=(\mathbf{x}-\mathbf{y}{)}^T{\mathbf{S}}^{-1}(\mathbf{x}-\mathbf{y})$
其中$\mathbf{S}$是样本协方差矩阵。
作为一种替代方法，可变性的度量可以直接并入距离度量中。
堪培拉度量
$d(\mathbf{x},\mathbf{y})=\sum _{i=1}^p\frac{|x_i-y_i|}{(x_i+y_i)}$
它仅为非负变量定义。

2.2 主成分分析

入侵检测问题中的数据具有高维性。为了便于探索和进一步分析，需要降低数据的维数。PCA通常用于此目的。主成分分析是通过几个新的变量来解释一组变量的方差结构，这些新变量是原始变量的线性组合。
主成分是p个随机变量X₁，X₂，…，X_p的特殊线性组合，具有三个重要性质：（1）主成分是不相关的，（2）第一主成分的方差最大，第二主成分的方差次之，依此类推，（3）所有变量的总方差主成分组合等于原始变量X₁，X₂，…，X_p的总变化量。通过对X₁，X₂，…，X_p的协方差矩阵或相关矩阵的特征分析，可以很容易地得到具有这种性质的新变量。
设原始数据X是一个n x p的矩阵，n个观察值包含p个变量(X₁，X₂，…，X_p)，令S是X₁，X₂，…，X_p的协方差矩阵，大小为p x p。如果(λ₁, e₁), (λ₂, e₂), …, (λ_p, e_p) 是矩阵S的p个特征向量对（特征值 - 特征向量），则第i个主成分为：
$y_i={\mathbf{e}}_i^T(\mathbf{x}-\overline{\mathbf{x}})=e_{i1}(x_1-\overline{x_1})+e_{i2}(x_2-\overline{x_2})+...+e_{ip}(x_p-\overline{x_p}),i=1,2,...,p$
其中
${\lambda }_1\ge {\lambda }_1\ge ...\ge {\lambda }_p\ge 0$
第i个主成分具有样本方差λ_i，而任意一对主成分的样本协方差为0.此外，如果变量X_i的样本方差为s_ii，则所有变量X₁，X₂，…，X_p的总样本方差为$\sum _{i=1}^p{s}_{ii}={\lambda }_1+{\lambda }_2+...+{\lambda }_p$
也就是所有主成分的总样本方差。意味着原始数据中的所有变化都是由主成分解释的。
可以以与协方差矩阵相同的方式对变量X₁，X₂，…，X_p的p x p相关矩阵R进行PCA。如果(λ₁, e₁), (λ₂, e₂), …, (λ_p, e_p) 是矩阵R的p个特征向量对（特征值 - 特征向量），则第i个主成分为：
$y_i={\mathbf{e}}_i\mathbf{z}=e_{i1}{z}_1+e_{i2}{z}_2+...+e_{ip}{z}_p,i=1,2,...,p$
其中${\mathbf{z}}^T$为标准化（归一化）观测值后的向量，
定义为$z_k=\frac{(xk-\overline{x_k})}{\sqrt{s_{kk}}},k=1,2,...,p$
样本相关矩阵的主成分仍然具有与以前相同的性质。即第i个主分量具有样本方差λ，任意一对主分量的样本协方差为0，所有主分量的总样本方差为${\lambda }_1+{\lambda }_2+...+{\lambda }_p=p$。
它是所有标准化变量z₁，z₂，…，z_p的总样本方差。
样本协方差矩阵S和样本相关矩阵R的主成分通常不相同。此外，它们也不是其他函数的简单函数。当某些变量的大小比其他变量大得多时，它们将在主成分中获得很重的权重。因此，如果变量是在范围大不相同的尺度上测量的，或者如果测量单位不相称，则最好对样本相关矩阵执行PCA。
大多数数据集包含一个或几个不寻常的观测值，这些观测值似乎不属于其他观测值产生的变异模式。当一个观测值与大多数数据不同或在假设的数据概率模型下不太可能时，它被认为是一个异常值。对于单一特征的数据，不寻常的观测是指相对于其他特征非常大或非常小的观测。如果假设正态分布，任何标准化值大于绝对值（例如，大于3或4）的观测值通常被识别为数据集的异常值。然而，由于有许多特点，情况变得复杂起来。在高维情况下，当分别考虑每个维时，可能会出现异常值，这些异常值不会作为离群观测值出现，因此不会从单变量标准中检测出来。因此，需要使用多元方法将所有特征综合考虑。
令${\mathbf{X}}_1,{\mathbf{X}}_2,...,{\mathbf{X}}_n$为多元分布的随机样本，均值向量为μ，协方差矩阵为Σ，其中，
${\mathbf{X}}_j^T=(X_{j1},X_{j2},...,X_{jp}),j=1,2,...,n$
通常用于检测多元异常值的过程是使用马氏距离测量每个观测点到数据中心的距离。统计量T²基于马氏距离。
$T^2=\frac{n}{n+1}(\mathbf{X}-\overline{\mathbf{X}}{)}^T{\mathbf{S}}^{-1}(\mathbf{X}-\overline{\mathbf{X}})$
分布为$\frac{(n-1)p}{n-p}{F}_{p,n-p}$，其中
$\overline{\mathbf{X}}=\frac{1}{n}\sum _{j=1}^n{\mathbf{X}}_j$ ，
$\mathbf{S}=\frac{1}{n-1}\sum _{j=1}^n({\mathbf{X}}_j-\overline{\mathbf{X}}){(X_j-\overline{X})}^T$
F_p,n-p表示具有p和n-p自由度的F分布随机变量。
T²的值越大，表示观测值X与总体中心的偏差越大。F检验统计量$\frac{(n-1)p}{n-p}{F}_{p,n-p}$可用于检测异常值。
代替马氏距离，我们还可以使用其他距离度量，如欧几里德距离和堪培拉度量。任何距离大于阈值的观测值都被视为异常值。阈值通常由距离的经验分布确定。这是因为即使在正态性假设下，也很难得出这些距离的分布。
长期以来，主成分分析一直被用于多变量离群点检测。对于一个观察值x，样本主成分为y₁,y₂,…,y_p。标准化主成分得分的平方和，
$\sum _{i=1}^p\frac{y_i^2}{{\lambda }_i}=\frac{y_1^2}{{\lambda }_1}+\frac{y_2^2}{{\lambda }_2}+...+\frac{y_p^2}{{\lambda }_p}$
等于观测值x与样本平均值的马氏距离。
前几个主成分的方差较大，解释了总样本方差的最大累积比例。这些主要成分往往与方差和协方差相对较大的特征密切相关。因此，与前几个分量相关的异常值通常对应于一个或多个原始变量的异常值。另一方面，最后几个主成分代表原始变量的最小方差的线性函数。这些成分对与数据的相关结构不一致的观测值很敏感，但不是原始变量的异常值。因此，次要成分上的较大观测值反映了使用基于原始变量大值的标准无法检测到的多元异常值。
由于样本主成分是不相关的，在正态假设和假设样本量较大的情况下，可以得出如下结论：
$\sum _{i=1}^q\frac{y_i^2}{{\lambda }_i}=\frac{y_1^2}{{\lambda }_1}+\frac{y_2^2}{{\lambda }_2}+...+\frac{y_q^2}{{\lambda }_q},q\le p$
具有自由度q的卡方分布。为实现这一点，还必须假设所有特征值都是不同的且为正的，即λ₁>λ₂>λ₃>…>λ_p>0。在给定显著性水平α的情况下，给出离群点检测准则：
如果$\sum _{i=1}^q\frac{y_i^2}{{\lambda }_i}>{\chi }_q^2(\alpha )$，观测值x为离群点。
式中${\chi }_q^2(\alpha )$是自由度为q的卡方分布的上α百分点。α的值表示将正常观测值分类为异常值时的错误概率（误报概率）。检查最后r个成分的平方和$\sum _{i=p-r+1}^q\frac{y_i^2}{{\lambda }_i}$对于确定观测值x的多少变化分布在这些后者分量上也很有用。当最后几个组成部分包含一个观测值中的大部分变化时，就相关结构而言，这表明该观测值是一个异常值。通过检查每个观测值x的$\sum _{i=p-r+1}^q\frac{y_i^2}{{\lambda }_i}$或$\underset{p-r+1\le i\le p}{\max }\left|\frac{y_i}{\sqrt{{\lambda }_i}}\right|$可以确定最后r个成分的相对重要性。

3. 提出的异常检测方案

如前所述，由于对数据的某些假设，图形化方法可能不适用于实时检测应用。另一方面，在我们提出的方案中，我们感兴趣的是开发一种基于主成分的异常检测方案，该方案可以实时应用，并且不会对数据施加太多的限制，其中PCA的应用目标是降低数据空间的维数，而不是作为异常点检测工具。
根据异常检测方法，我们假设异常在性质上与正常情况不同。也就是说，与已建立的正常模式的较大偏差可以标记为攻击。没有试图区分不同类型的攻击。为了建立检测算法，我们对正常组的相关矩阵进行主成分分析。使用相关矩阵是因为每个特征是在不同的尺度上测量的。训练数据在用于确定检测标准之前必须没有异常值，因为异常值会带来方差、协方差和相关性的大幅增加。这些变异和协变度量的相对大小对主成分解有显著影响，特别是对前几个成分。因此，用相关矩阵的稳健估计开始主成分分析是有价值的。一个简单的方法，以获得稳健估计是多元修剪。首先，我们使用马氏度量来确定要修剪的100γ%极端观测值。从常规估计量$\mathbf{x}$和S开始，计算每个观测值x_i（i = 1,2，…，n）的距离$d_i^2=({\mathbf{x}}_i-\overline{\mathbf{x}}{)}^T{\mathbf{S}}^{-1}({\mathbf{x}}_i-\overline{\mathbf{x}})$。对于给定的γ（在我们的实验中为0.005），确定了对应于γ*n最大值{$d_i^2$，i=1,2,…,n}的观测值。然后根据剩余观测值计算均值和协方差矩阵的新修剪估计量$\mathbf{x}$和S。利用S的元素得到了相关矩阵的鲁棒估计量，可以重复修剪过程，以确保估计量$\mathbf{x}$和S抵抗异常值。只要修剪后剩余的观测数超过p（向量$\mathbf{x}$的维数），多元修剪确定的估计量S将是正定的。
顺带一提，这种稳健的过程使得我们的方法非常适合于无监督异常检测。我们不能期望训练数据总是由正常实例组成。一些可疑的数据或入侵可能隐藏在数据集中。然而，为了使异常检测工作，我们假设正常实例的数量必须远远大于异常的数量。因此，通过如上所述的微调过程，异常将被捕获并从训练数据集中移除。我们提出的方案中，主成分分类器（PCC）由两个主成分得分函数组成，一个来自主要成分$\sum _{i=1}^q\frac{y_i^2}{{\lambda }_i}$，一个来自次要成分$\sum _{i=p-r+1}^q\frac{y_i^2}{{\lambda }_i}$。第一个函数是一个在文献中经常使用的方法，是在一些原始特征上检测具有大值的极端观测值。与其他现有方法不同，我们建议在第一个函数之外使用第二个函数来帮助检测不符合正态相关结构的观测值。与其他方法相比，我们提出的方法的一个明显优点是它保留了关于异常值性质的信息，无论它们是极值还是与正常情况不具有相同的相关结构。
主要成分的数量q将根据这些成分所占训练数据的变化量来确定。基于我们的实验，我们建议使用q个主要成分来解释标准化特征中大约50%的总方差。当原始特征不相关时，相关矩阵的每个主成分的特征值都等于1。因此，PCC中使用的r次分量是方差或特征值小于0.20的分量，表示特征之间的某种关系。
使用PCC的分类方案如下。计算要确定类别的观察值x的主成分得分。
如果$\sum _{i=1}^q\frac{y_i^2}{{\lambda }_i}>c_1$或$\sum _{i=p-r+1}^p\frac{y_i^2}{{\lambda }_i}>c_2$，实例x为异常。
如果$\sum _{i=1}^q\frac{y_i^2}{{\lambda }_i}\le c_1$或$\sum _{i=p-r+1}^p\frac{y_i^2}{{\lambda }_i}\le c_2$，实例x为正常。
其中c₁和c₂是离群阈值，使得分类器将产生指定的虚警率。假设数据分布为多元正态分布，则该分类器的虚警率为α=α₁+α₂-α₁α₂，其中α₁和α₂分别为主要成分函数和次要成分函数的虚警率。在其他情况下，Cauchy-Schwartz不等式和Bonferroni不等式提供了虚警率α的一个下界和一个上界。
${\alpha }_1+{\alpha }_2-\sqrt{{\alpha }_1{\alpha }_2}\le \alpha \le {\alpha }_1+{\alpha }_2$
选择α₁和α₂的值是为了反映我们想要检测的异常值类型的相对重要性。在我们的实验中，我们选择α₁=α₂。例如，为了实现约2%的虚警率，根据上式，我们设置α₁=α₂=0.0101。由于正态性假设可能被违背，我们选择根据经验公式设置离群值阈值，选择基于训练数据中的$\sum _{i=1}^q\frac{y_i^2}{{\lambda }_i}$和$\sum _{i=p-r+1}^q\frac{y_i^2}{{\lambda }_i}$而非卡方分布。也就是说，c₁和c₂分别是$\sum _{i=1}^q\frac{y_i^2}{{\lambda }_i}$和$\sum _{i=p-r+1}^q\frac{y_i^2}{{\lambda }_i}$的0.9899分位数。

4. 实验

通过与基于密度的局部离群点（LOF）方法和其他两种基于距离的入侵检测方法（Canberra度量和欧氏距离）的比较，研究了PCC方法的性能。基于欧氏距离的方法实际上是k-近邻法。我们选择k=1和5进行比较研究。
实验在以下框架下进行：
1） 所有的异常值阈值都是由训练数据确定的。我们把误报率从1%调整到10%。对于PCC方法，阈值的选择应确保α₁=α₂。
2） 训练和测试数据均来自KDD’99训练数据集。
3） 每个训练数据集由5000个正常连接组成，通过系统抽样从KDD’99数据中的所有正常连接中随机选择。
4） 为了评估分类器的准确性，我们用5个不同的训练样本进行了5个独立的实验。在每一个实验中，分类器用从KDD’99数据中随机选择的92279个正常连接和39674个攻击连接的测试集进行测试。

4.1 KDD’99数据

KDD CUP 1999数据(http://kdd.ics.uci.edu/databases/kddcup99/kddcup99.html) 是第三届国际知识发现和数据挖掘工具竞赛使用的数据集。训练数据集包含494021条连接记录，测试数据集包含311029条与训练数据不来自同一概率分布的记录。由于概率分布是不同的，在我们的实验中，我们只从训练数据集中采样数据，并在训练和测试阶段使用。
连接是包含41个特征值的TCP数据包序列，标记为正常或攻击，只有一种特定的攻击类型。训练数据中有22种攻击类型。然而，在本研究中，我们将他们视为同一攻击组。这41个特征可分为三组：第一组是单个TCP连接的基本特征，第二组是由领域知识建议的连接内的内容特征，第三组是使用两秒时间窗口计算的流量特征。在41个特征中，34个是数字特征，7个是符号特征。实验中只使用了34个数字特征。KDD CUP 1999数据中列出了完整的功能和详细信息。

4.2 性能指标

分类结果通常以一个称为混淆矩阵的矩阵表示。分类器的准确度是由其错误分类率或正确分类率来衡量的（TP、FN、FP、TN）。
另外两个性能指标，精确率和召回率也值得关注。
精确率Precision=TP/（TP+FP）
召回率Recall=TP/（TP+FN）
评估异常检测方案的另一个有价值的工具是接受者操作特征曲线（Receiver OperatingCharacteristic, ROC）曲线，它是检测率与虚警率之间的关系图。方案的ROC曲线越靠近左上角，方案的性能越好。如果不同方案的roc相互叠加，则这些方案具有相同的性能。

5. 实验结果与讨论

从以上实验结果来看，本文提出的基于主成分的异常检测方案能够有效地识别攻击。在我们的研究中，唯一可比较的竞争对手是LOF方法，但只有当虚警率为4%或更高时。该方案不仅具有良好的查准率和查全率，而且能够将虚警保持在期望的水平。在KDD’99数据的实验中，PCC具有足够的灵敏度来检测攻击。另外，与马氏距离不同，PCC通过使用两种不同的主成分函数提供了更多关于攻击性质的信息。PCC的另一个好处是，在检测阶段，统计数据可以在较少的时间内计算出来，这使得实时使用该方法成为可能。这是因为PCC中只使用了三分之一的主成分，5个主成分解释了34个特征总变化的50%，6-7个次成分的特征值小于0.20。

6. 结论

随着越来越多的组织变得容易受到各种各样的网络威胁，因此有一种能够区分攻击和正常连接的有效算法是非常重要的。在异常检测方法的基础上，研究了鲁棒PCA在异常检测中的应用，并将其应用于异常检测问题。预测模型由正态连接主成分的两个函数组成，其中主要成分约占总变异的50%，次要成分的特征值小于0.20。这种方法的另一个好处是，它能够区分异常的性质，无论它们在极值或不同的相关结构方面是否不同于正常情况。对KDD’99数据的实验表明，本文提出的异常检测方案比其他方法具有更好的性能。其检测率接近99%，虚警率低至1%。由于其鲁棒性，我们提出的方案也适用于无监督的训练数据。

附：方法复现

方法复现数据和源码已上传至Github仓库(https://github.com/xzydn/anomaly-detection-example/tree/master/PrincipalComponentClassifier(PCC))
首先载入数据

import pandas as pd
import numpy as np

# 准备训练数据
origin_training_data = pd.read_csv(r'kddcup.data.corrected', header=None)
print(origin_training_data)

# 删除符号列   
origin_training_data_1 = origin_training_data.drop(columns=[1, 2, 3, 6, 11, 20, 21])
print(origin_training_data_1)

# 筛选正常行
origin_training_data_1_normal = origin_training_data_1.loc[origin_training_data_1.loc[:, 41] == 'normal.']
print(origin_training_data_1_normal)

# 筛选异常行
origin_training_data_1_abnormal = origin_training_data_1.loc[origin_training_data_1.loc[:, 41] != 'normal.']
print(origin_training_data_1_abnormal)

# 转为numpy矩阵
origin_training_data_2_abnormal = origin_training_data_1_abnormal.drop(columns=[41])
X_train_abnormal = np.array(origin_training_data_2_abnormal)
print(X_train_abnormal)

为了避免变量范围尺度相差较大，对样本进行标准化（归一化）

mu_normal = np.mean(X_train_normal, axis=0) # 均值
print(mu_normal)
sigma_normal = np.std(X_train_normal, axis=0) # 标准差
print(sigma_normal)
Z_train_normal = (X_train_normal - mu_normal)/(sigma_normal + np.spacing(1)) # 归一化
print(Z_train_normal)
print(np.mean(Z_train_normal, axis=0))
print(np.std(Z_train_normal, axis=0))

在这里基于numpy实现PCA

#基于numpy实现PCA
import numpy as np

covMat_Z_train_normal = np.cov(Z_train_normal,rowvar=False) #计算协方差矩阵，rowvar=0表示数据的每一列代表一个feature
print(covMat_Z_train_normal)
(featValue_Z_train_normal, featVec_Z_train_normal) = np.linalg.eig(covMat_Z_train_normal) #计算协方差矩阵的特征值和特征向量
print("协方差矩阵的特征值：")
print(featValue_Z_train_normal)
print("协方差矩阵的特征向量：")
print(featVec_Z_train_normal)

index_Z_train_normal=np.argsort(featValue_Z_train_normal) #将特征值按从小到大排序，index保留的是对应原featValue中的下标
print(index_Z_train_normal)

featValue_Z_train_normal_sorted = featValue_Z_train_normal[index_Z_train_normal[-1::-1]] # 从大到小排序后的特征值
print("特征值：")
print(featValue_Z_train_normal_sorted)
featVec_Z_train_normal_sorted = featVec_Z_train_normal[:, index_Z_train_normal[-1::-1]] # 从大到小排序后的特征向量
print("特征向量：")
print(featVec_Z_train_normal_sorted)

# 特征值之和等于方差之和
sum_featValue_Z_train_normal = np.sum(featValue_Z_train_normal_sorted) 
print(sum_featValue_Z_train_normal)
print(np.sum(np.std(Z_train_normal, axis=0))) # 方差已归一化

print(np.sum(featValue_Z_train_normal_sorted[:5])/sum_featValue_Z_train_normal)

计算主要成分和次要成分

# 主要成分
n_featValue_Z_train_normal_sorted = featValue_Z_train_normal_sorted[:5] # 取前5个成分
print("特征值：")
print(n_featValue_Z_train_normal_sorted)
n_featVec_Z_train_normal_sorted = featVec_Z_train_normal_sorted[:, :5]
print("特征向量：")
print(n_featVec_Z_train_normal_sorted)

Z_train_normal_pc = np.dot(Z_train_normal,n_featVec_Z_train_normal_sorted)
print(Z_train_normal_pc)

Z_train_normal_scores = np.power(Z_train_normal_pc, 2)
Z_train_normal_scores = Z_train_normal_scores/n_featValue_Z_train_normal_sorted
Z_train_normal_scores = np.sum(Z_train_normal_scores, axis=1)

print("主要成分得分：")
print(Z_train_normal_scores)

index_Z_train_normal_scores = np.argsort(Z_train_normal_scores)
Z_train_normal_scores_sorted = Z_train_normal_scores[index_Z_train_normal_scores[-1::-1]]
print(Z_train_normal_scores_sorted)

Z_train_normal_scores_choice = np.int(np.size(Z_train_normal_scores_sorted)*0.01)
print(Z_train_normal_scores_choice)
Z_train_normal_scores_threshold = Z_train_normal_scores_sorted[Z_train_normal_scores_choice]
print("主要组成成分阈值：", Z_train_normal_scores_threshold)

# 次要成分
n1_featValue_Z_train_normal_sorted = featValue_Z_train_normal_sorted[-1:-4:-1] # 取后3个成分
print("特征值：")
print(n1_featValue_Z_train_normal_sorted)
n1_featVec_Z_train_normal_sorted = featVec_Z_train_normal_sorted[:, -1:-4:-1]
print("特征向量：")
print(n1_featVec_Z_train_normal_sorted)

Z_train_normal_pc1 = np.dot(Z_train_normal,n1_featVec_Z_train_normal_sorted)
print(Z_train_normal_pc1)

Z_train_normal_scores1 = np.power(Z_train_normal_pc1, 2)
Z_train_normal_scores1 = Z_train_normal_scores1/(n1_featValue_Z_train_normal_sorted + np.spacing(1))
Z_train_normal_scores1 = np.sum(Z_train_normal_scores1, axis=1)

print("次要成分得分：")
print(Z_train_normal_scores1)

index_Z_train_normal_scores1 = np.argsort(Z_train_normal_scores1)
Z_train_normal_scores1_sorted = Z_train_normal_scores1[index_Z_train_normal_scores1[-1::-1]]
print(Z_train_normal_scores1_sorted)

Z_train_normal_scores1_choice = np.int(np.size(Z_train_normal_scores1_sorted)*0.01)
print(Z_train_normal_scores1_choice)
Z_train_normal_scores1_threshold = Z_train_normal_scores1_sorted[Z_train_normal_scores1_choice]
print("次要组成成分阈值：", Z_train_normal_scores1_threshold)

计算异常检出率

Z_train_abnormal = (X_train_abnormal - mu_normal)/(sigma_normal + np.spacing(1))
Z_train_abnormal_pc1 = np.dot(Z_train_abnormal,n1_featVec_Z_train_normal_sorted)
Z_train_abnormal_scores1 = np.power(Z_train_abnormal_pc1, 2)
Z_train_abnormal_scores1 = Z_train_abnormal_scores1/(n1_featValue_Z_train_normal_sorted + np.spacing(1))
Z_train_abnormal_scores1 = np.sum(Z_train_abnormal_scores1, axis=1)

print(Z_train_abnormal_scores1 > Z_train_normal_scores1_threshold)

Z_train_abnormal_res1 = Z_train_abnormal_scores[np.logical_or((Z_train_abnormal_scores > Z_train_normal_scores_threshold) , (Z_train_abnormal_scores1 > Z_train_normal_scores1_threshold))]
print(np.size(Z_train_abnormal_res1))

print("检出率：", np.size(Z_train_abnormal_res1)/np.size(Z_train_abnormal_scores1))

计算正常误报率

Z_train_normal = (X_train_normal - mu_normal)/(sigma_normal + np.spacing(1))
Z_train_normal_pc = np.dot(Z_train_normal,n_featVec_Z_train_normal_sorted)
Z_train_normal_scores = np.power(Z_train_normal_pc, 2)
Z_train_normal_scores = Z_train_normal_scores/n_featValue_Z_train_normal_sorted
Z_train_normal_scores = np.sum(Z_train_normal_scores, axis=1)
Z_train_normal_res = Z_train_normal_scores[Z_train_normal_scores > Z_train_normal_scores_threshold]
print(np.size(Z_train_normal_res))

Z_train_normal_pc1 = np.dot(Z_train_normal,n1_featVec_Z_train_normal_sorted)
Z_train_normal_scores1 = np.power(Z_train_normal_pc1, 2)
Z_train_normal_scores1 = Z_train_normal_scores1/(n1_featValue_Z_train_normal_sorted + np.spacing(1))
Z_train_normal_scores1 = np.sum(Z_train_normal_scores1, axis=1)
Z_train_normal_res = Z_train_normal_scores[Z_train_normal_scores1 > Z_train_normal_scores1_threshold]
print(np.size(Z_train_normal_res))

Z_train_normal_res1 = Z_train_normal_scores[np.logical_or((Z_train_normal_scores > Z_train_normal_scores_threshold) , (Z_train_normal_scores1 > Z_train_normal_scores1_threshold))]
print("误报率：", np.size(Z_train_normal_res1)/np.size(Z_train_normal_scores1))

复现程序仅供参考。检出率99.36%，虚警率1.93%，与文献结果差异较小。