异常检测Task03
本次学习参照Datawhale开源学习地址：https://github.com/datawhalechina/team-learning-data-mining/tree/master/AnomalyDetection
本次学习分为五个章节：
一、概述
二、基于统计学的方法
三、线性模型
四、基于邻近度的方法
五、集成方法

本章主要内容包括：

• 线性回归
• 主成分分析

概述

  真实数据集中各特征通常具有一定的相关性，且每个特征本身分布很难转换成高斯分布，在这种情况下上一章中的统计学方法很多时候表现欠佳。因此，我们采用另一种参数化的相关性分析方法——回归建模。
  回归建模主要分为两类：线性回归（通过其他变量预测单独的属性值）， 主成分分析（用一些潜在变量来代表整个数据）。

1、线性回归

  在线性回归中，我们假设不同维度的变量具有一定的相关性，并可以通过一个相关系数矩阵进行衡量，然后通过一系列自变量去预测一个特殊因变量的值。对于特定的观测值，可以通过线性方程组来建模。在实际应用中，观测值的数量往往远大于数据的维度，导致线性方程组是一个超定方程，不能直接求解。因此需要通过优化的方法，最小化模型预测值与真实数据点的误差。

2.1 基于自变量与因变量的线性回归

2.1.1 最小二乘法

自变量与因变量的关系为：

Y = ∑ i = 1 d a i ⋅ X i + a d + 1 Y=\sum_{i=1}^{d} a_{i} \cdot X_{i}+a_{d+1} Y=i=1∑d​ai​⋅Xi​+ad+1​

  变量Y为因变量，也就是我们要预测的值； X 1 . . . X d X_{1}...X_{d} X1​...Xd​为一系列自变量，也就是输入值。系数 a 1 . . . a d + 1 a_{1}...a_{d+1} a1​...ad+1​为要学习的参数。假设维度为 d d d的数据共包含 N N N个样本，第 j j j个样本包含的数据为 x j 1 . . . x j d x_{j1}...x_{jd} xj1​...xjd​和 y j y_{j} yj​，带入式(1)如下式所示：

y j = ∑ i = 1 d a i ⋅ x j i + a d + 1 + ϵ j y_{j}=\sum_{i=1}^{d} a_{i} \cdot x_{j i}+a_{d+1}+\epsilon_{j} yj​=i=1∑d​ai​⋅xji​+ad+1​+ϵj​

  这里 ϵ j \epsilon_{j} ϵj​为第 j j j个样本的误差。以 Y Y Y 代表 N × 1 N \times 1 N×1 的因变量矩阵 ( y 1 . . . y N ) T {(y_{1}...y_{N})}^{T} (y1​...yN​)T，即样本中的真实值；以 U U U代表 N × ( d + 1 ) N \times (d+1) N×(d+1)的自变量矩阵，其中第 j j j行为 ( x j 1 . . . x j d , 1 ) (x_{j1}...x_{jd}, 1) (xj1​...xjd​,1)；以 A A A 代表 ( d + 1 ) × 1 (d+1) \times 1 (d+1)×1 的系数矩阵 ( a 1 . . . a d + 1 ) T (a_{1}...a_{d+1})^{T} (a1​...ad+1​)T。则模型可表示为：
f ( U , A ) = U ⋅ A f(U, A) = U \cdot A f(U,A)=U⋅A

  定义目标函数为：

L ( A ) = 1 2 ∥ Y − U ⋅ A ∥ 2 L(A) = \frac{1}{2}{\left\| {Y - U \cdot A} \right\|^2} L(A)=21​∥Y−U⋅A∥2

  目标函数是关于 A A A的凸函数，其对 A A A求偏导为：

∂ L ( A ) ∂ A = 1 2 ∂ ∥ Y − U ⋅ A ∥ 2 ∂ A = − U T ( Y − U ⋅ A ) \frac{{\partial L(A)}}{{\partial A}} = \frac{1}{2}\frac{{\partial {{\left\| {Y - U \cdot A} \right\|}^2}}}{{\partial A}} = - {U^T}(Y - U \cdot A) ∂A∂L(A)​=21​∂A∂∥Y−U⋅A∥2​=−UT(Y−U⋅A)

  令 ∂ L ( A ) ∂ A = 0 \frac{{\partial L(A)}}{{\partial A}}=0 ∂A∂L(A)​=0，得到最优参数为：

A = ( U T ⋅ U ) − 1 ⋅ ( U T ⋅ Y ) A=\left(U^{T} \cdot U\right)^{-1} \cdot\left(U^{T} \cdot Y\right) A=(UT⋅U)−1⋅(UT⋅Y)

  这种求解线性回归参数的方法也叫最小二乘法。

  最小二乘法要求矩阵 U T ⋅ U U^{T} \cdot U UT⋅U 可逆，即 U T ⋅ U U^{T} \cdot U UT⋅U是满秩的。当 U T ⋅ U U^{T} \cdot U UT⋅U不可逆时可以使用梯度下降法。

2.1.2 梯度下降法

假设用均方误差作为损失函数：

l ( i ) ( w , b ) = 1 2 ( y ^ ( i ) − y ( i ) ) 2 l^{(i)}(\mathbf{w}, b)=\frac{1}{2}\left(\hat{y}^{(i)}-y^{(i)}\right)^{2} l(i)(w,b)=21​(y^​(i)−y(i))2

L ( w , b ) = 1 n ∑ i = 1 n l ( i ) ( w , b ) = 1 n ∑ i = 1 n 1 2 ( w ⊤ x ( i ) + b − y ( i ) ) 2 L(\mathbf{w}, b)=\frac{1}{n} \sum_{i=1}^{n} l^{(i)}(\mathbf{w}, b)=\frac{1}{n} \sum_{i=1}^{n} \frac{1}{2}\left(\mathbf{w}^{\top} \mathbf{x}^{(i)}+b-y^{(i)}\right)^{2} L(w,b)=n1​i=1∑n​l(i)(w,b)=n1​i=1∑n​21​(w⊤x(i)+b−y(i))2
  其中 y ^ \hat{y} y^​ 为预测值， y y y 为真实值。

  当模型和损失函数形式较为简单时，上面的误差最小化问题的解可以直接用梯度下降法公式表达出来。但大多数模型并没有解析解，只能通过优化算法有限次迭代模型参数来尽可能降低损失函数的值，求出数值解。
代码实现

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
"""生成测试数据"""
X = np.array([[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0],
                    [2.0, 4.0, 6.0, 8.0, 18.0, 12.0, 14.0, 5.0, 18.0, 20.0]])


plt.scatter(X[0],X[1])
plt.show()

"""利用训练集拟合线性回归模型"""
lr_train = LinearRegression().fit(X[0].reshape(-1, 1), X[1].reshape([-1,1]).reshape(-1))

xx = np.linspace(0, 10, 100)
yy = lr_train.predict(xx.reshape(-1, 1))
plt.scatter(xx,yy)
plt.show()

"""用残差的1sigma筛选异常值"""
residual = X[1] - lr_train.predict(X[0].reshape(-1, 1))
threshold_l = residual.mean() -  residual.std()
threshold_h = residual.mean() +  residual.std()
error = ((residual>threshold_h) | (residual<threshold_l))*X[0]
print("异常点为：")
for i,j in zip(error,X[1]):
    if i !=0:
        plt.scatter(i,j,c='r')
        plt.show()
        print([i,j])

异常点为：
[5.0, 18.0]
[8.0, 5.0]

2.2 基于异常检测的线性回归

  上一节中我们将一个特定的变量（因变量）认为是特殊的，最优平面是通过最小化该特殊变量的均方误差而确定的。而我们通常所说的异常检测中并不会对任何变量给与特殊对待，异常值的定义是基于基础数据点的整体分布，因此需要采用一种更一般的回归建模：即以相似的方式对待所有变量，通过最小化数据对该平面的投影误差确定最佳回归平面。在这种情况下，假设我们有一组变量 X 1 … X d X_{1}… X_{d} X1​…Xd​， 对应的回归平面如下：

a 1 ⋅ X 1 + … + a d ⋅ X d + a d + 1 = 0 a_{1} \cdot X_{1}+\ldots+a_{d} \cdot X_{d}+a_{d+1}=0 a1​⋅X1​+…+ad​⋅Xd​+ad+1​=0

  为了后续计算的方便，对参数进行如下约束：
∑ i = 1 d a i 2 = 1 \sum\limits_{i = 1}^d {a_i^2 = 1} i=1∑d​ai2​=1
  以 L 2 L_{2} L2​范数作为目标函数：
L = ∥ U ⋅ A ∥ 2 L = {\left\| {U \cdot A} \right\|_2} L=∥U⋅A∥2​

  这样的一个问题可以通过主成分分析方法解决。

2、主成分分析

  PCA方法是通过矩阵变换将高维过程数据投影到正交低维子空间，并保留主要过程信息。也是将各维度相关的原始数据变换成各维度线性无关的表示过程。

2.1 PCA算法过程

• 首先，特征归一化
  
  计算出每个特征的均值( u j uj uj),然后用 x j − u j xj-uj xj−uj来替代，这样归一化后每个特征的均值就为0。

• 其次，求协方差矩阵
  

• 最后，求协方差矩阵的特征值及特征向量
  
  其中，S为对角矩阵，其对角线上的数就是协方差矩阵的特征值，而U就是协方差矩阵的特征向量。而U的前 k k k列就是我们用于代替原来n个特征的新特征。这样就将原始 n n n维数据，经过用变换后变为 k k k维。
  

2.2 PCA用于异常检测

  PCA在异常检测方面的做法，大体有两种思路：一种是将数据映射到低维特征空间，然后在特征空间不同维度上查看每个数据点跟其它数据的偏差；另外一种是将数据映射到低维特征空间，然后由低维特征空间重新映射回原空间，尝试用低维特征重构原始数据，看重构误差的大小。两种思路看似不太一样，其实本质上是差不多的。

• 思路一
  对于某一个特征向量 e j ej ej，数据样本 x i xi xi在该方向上的偏离程度dij可以用如下的公式计算：
  
  这里的 λ j λj λj主要起归一化的作用，这样可以使得不同方向上的偏离程度具有可比性。在计算了数据样本在所有方向上的偏离程度之后，为了给出一个综合的异常得分，最自然的做法是将样本在所有方向上的偏离程度加起来，即：
  
  这个公式只是计算异常得分的一种方式。也有一些算法采取了略微不同的做法，比如，有的只考虑数据在前 k k k个特征向量方向上的偏差，或者只考虑后 r r r个特征向量方向上的偏差，即：
  
  这里的 C 1 C1 C1和 C 2 C2 C2是人为设定的两个阈值，如果得分大于阈值则判断为异常。
    一般而言，前几个特征向量往往直接对应原始数据里的某几个特征，在前几个特征向量方向上偏差比较大的数据样本，往往就是在原始数据中那几个特征上的极值点。而后几个特征向量有些不同，它们通常表示某几个原始特征的线性组合，线性组合之后的方差比较小反应了这几个特征之间的某种关系。在后几个特征方向上偏差比较大的数据样本，表示它在原始数据里对应的那几个特征上出现了与预计不太一致的情况。到底是考虑全部特征方向上的偏差，前几个特征向量上的偏差，还是后几个特征向量上的偏差，在具体使用时可以根据具体数据灵活处理。
• 思路二
  PCA提取了数据的主要特征，如果一个数据样本不容易被重构出来，表示这个数据样本的特征跟整体数据样本的特征不一致，那么它显然就是一个异常的样本。对于数据样本 X i Xi Xi, 假设其基于 k k k维特征向量重构的样本为 X ′ i X'i X′i, 则该数据样本的异常得分可以用如下的公式计算：
  
  
  上面的公式考虑了重构使用的特征向量的个数 k k k的影响，将 k k k的所有可能做了一个加权求和，得出了一个综合的异常得分。
    我们在基于低维特征进行数据样本的重构时，舍弃了较小的特征值对应的特征向量方向上的信息。因此，重构误差其实主要来自较小的特征值对应的特征向量方向上的信息。PCA在异常检测上的两种不同思路都会特别关注较小的特征值对应的特征向量，本质上是相似的。
  代码实现

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib.font_manager

from pyod.models.pca import PCA
from pyod.utils.data import generate_data, get_outliers_inliers

"""生成二维数据"""
outlier_fraction = 0.1
X_train, Y_train = generate_data(n_train=1000,train_only=True, contamination=outlier_fraction, n_features=2, offset = 10, random_state = 0)
# 拆分异常数据和正常数据
x_outliers, x_inliers = get_outliers_inliers(X_train,Y_train)

"""训练模型"""
clf = PCA(contamination=outlier_fraction)
clf.fit(x_inliers)

"""预测异常点"""
scores_pred = clf.decision_function(X_train)*-1
y_pred = clf.predict(X_train)

"""训练数据中分错类的点数目"""
n_errors = (y_pred != Y_train).sum()
print('No. of Errors : ','Principal Component Analysis (PCA)', n_errors)

"""可视化"""
xx , yy = np.meshgrid(np.linspace(-10, 10, 200), np.linspace(-10, 10, 200))


plt.figure(figsize=(10, 10))
threshold = stats.scoreatpercentile(scores_pred,100 *outlier_fraction)

Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()]) * -1
Z = Z.reshape(xx.shape)
subplot = plt.subplot(1, 1, 1)
subplot.contourf(xx, yy, Z, levels = np.linspace(Z.min(), threshold, 10),cmap=plt.cm.Blues_r)

a = subplot.contour(xx, yy, Z, levels=[threshold],linewidths=2, colors='red')
subplot.contourf(xx, yy, Z, levels=[threshold, Z.max()],colors='orange')
n_inliers = len(x_inliers)
n_outliers = len(x_outliers)
b = subplot.scatter(X_train[:-n_outliers, 0], X_train[:-n_outliers, 1], c='white',s=20, edgecolor='k') 
c = subplot.scatter(X_train[-n_outliers:, 0], X_train[-n_outliers:, 1], c='black',s=20, edgecolor='k')

subplot.axis('tight')
subplot.legend(
    [a.collections[0], b, c],
    ['decision boundary', 'true inliers', 'true outliers'],
    prop=matplotlib.font_manager.FontProperties(size=10),
    loc='lower right')
subplot.set_title('Principal Component Analysis (PCA)')
subplot.set_xlim((-10, 10))
subplot.set_ylim((-10, 10))

plt.show() 

No. of Errors :  Principal Component Analysis (PCA) 91

3、回归分析的局限性

  为了使回归分析技术有效，数据需要高度相关，并沿着低维子空间对齐。当数据不相关，但在某些区域高度聚集时，这种方法可能不会有效。

  另外，数据中的相关性在本质上可能不是全局性的。子空间相关性可能是特定于数据的特定位置的。在这种情况下，由主成分分析发现的全局子空间对于异常检测是次优的。因此，为了创建更一般的局部子空间模型，有时将线性模型与邻近模型结合起来是有用的。

总结

  当数据不同属性之间具有显著的相关性时，回归建模会是一种非常有效的异常检测方法。

参考资料

[1] 《Outlier Analysis》——Charu C. Aggarwal
[2] Anomaly Detection: A Survey VARUN CHANDOLA, ARINDAM BANERJEE, and VIPIN KUMAR University of Minnesota
[3] Anomaly Detection: A Tutorial
[4] Data Mining Concepts and Techniques Third Edition
[5] http://cs229.stanford.edu/
[6] https://mp.weixin.qq.com/s/DUgvYIL4Nb0UAMaluO-dnA