异常检测（二）

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

主要内容包括：

• 一元高斯分布
• 多元高斯分布
• 箱线图
• HBOS模型

概述

统计学方法是对数据的正常性做出假定。假定正常的数据对象由一个统计模型产生，违背该模型的数据是异常点。统计学方法的有效性高度依赖于对给定数据所做的统计模型假定是否成立。

根据如何指定和学习模型，异常检测的统计学方法可以划分为两个主要类型：参数方法和非参数方法。

参数方法假定正常的数据对象被一个以$\Theta$为参数的参数分布产生。该参数分布的概率密度函数$f(x,\Theta )$给出对象$x$被该分布产生的概率。该值越小，$x$越可能是异常点。（用样本估计总体分布）

非参数方法并不假定先验统计模型，而是试图从输入数据确定模型。通常假定参数的个数和性质都是不确定的。（把样本当作总体分布）

1、参数方法

1.1 基于正态分布的一元异常点检测

只有一个属性或特征的数据称为一元数据。假定数据由正态分布产生，然后可以由输入数据学习正态分布的参数，并把低概率的点识别为异常点。

假定输入数据集为 { x ( 1 ) , x ( 2 ) , . . . , x ( m ) } \{x^{(1)}, x^{(2)}, ..., x^{(m)}\} {x(1),x(2),...,x(m)}，数据集中的样本服从正态分布，即$x^{(i)}\sim N(\mu ,{\sigma }^2)$，我们可以根据样本求出参数$\mu$和$\sigma$。

$\mu =\frac{1}{m}{\sum }_{i=1}^m{x}^{(i)}$

${\sigma }^2=\frac{1}{m}{\sum }_{i=1}^m(x^{(i)}-\mu {)}^2$

求出参数之后，我们就可以根据概率密度函数计算数据点概率密度。

$p(x)=\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(x-\mu {)}^2}{2{\sigma }^2})$

如果计算出来的概率低于阈值 ε，就可以认为该数据点为异常点。需要注意的是，概率密度函数值和概率值不是一个概念，具体参考吴恩达《机器学习》异常检测章节。

阈值是个经验值，可以选择在验证集上使得评估指标值最大（也就是效果最好）的阈值取值作为最终阈值。

例如常用的3sigma原则中，如果数据点超过范围$(\mu -3\sigma ,\mu +3\sigma )$，那么这些点很有可能是异常点。

代码实现

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import kstest

np.random.seed(0)

"""生成标准正态分布数据"""
data = pd.DataFrame(np.random.randn(1000),columns = ['value'])
# 计算均值
u = data['value'].mean()
# 计算标准差
std = data['value'].std()

"""绘制直方图"""
data.hist(bins=50,alpha = 0.5)
plt.grid()
plt.show()

"""K-S检验是否符合正态分布"""
# 计算p值
p = kstest(data, 'norm', (u, std))[1] 
# 当p<=0.05 服从正态分布
if p<=0.05:
    print('数据服从正态分布')
else:
    print('数据不服从正态分布')

"""3sigma原则检测异常值""" 
if p<=0.05:  
    error = data[np.abs(data['value'] - u) > 3 * std]
    print("异常数据为：{}".format(error))
else:
    print('数据不服从正态分布，不能使用该检测方法')

数据服从正态分布
异常数据为：        value
589 -3.046143

1.2 箱线图
当一元分布不符合高斯分布时，也可以采用箱线图。箱线图对数据分布做了一个简单的统计，利用数据集的上下四分位数（Q1和Q3）、中点等形成。异常点常被定义为小于Q1－1.5IQR或大于Q3+1.5IQR的数据，其中IQR = Q3 - Q1。

代码实现

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

np.random.seed(0)

"""生成标准正态分布数据"""
data = pd.DataFrame(np.random.randn(1000),columns = ['value'])

"""绘制箱线图"""
plt.boxplot(data.value)
plt.show()

"""箱线图检测异常值""" 
# 计算下四分位数和上四分位
Q1 = data.value.quantile(q=0.25)
Q3 = data.value.quantile(q=0.75)

# 计算上下限
threshold_l = Q1 - 1.5 * (Q3 - Q1)
threshold_h = Q3 + 1.5 * (Q3 - Q1)

# 计算异常值  
error = data[(data.value>threshold_h) | (data.value<threshold_l)]
print("异常数据为：{}".format(error))

异常数据为：        value
271 -2.772593
334 -2.659172
427 -2.739677
494  2.696224
589 -3.046143
685 -2.834555
898  2.594425
943  2.759355

1.3 多个不相关特征，且符合多元高斯分布

涉及两个或多个属性或特征的数据称为多元数据。当各个属性或特征不相关时，可以把多元异常点检测任务转换成一元异常点检测问题。基于正态分布的一元异常点检测扩充到多元情形时，可以求出每一维度的均值和标准差。对于第$j$维：

${\mu }_j=\frac{1}{m}{\sum }_{i=1}^m{x}_j^{(i)}$

${\sigma }_j^2=\frac{1}{m}{\sum }_{i=1}^m(x_j^{(i)}-{\mu }_j{)}^2$

计算概率时的概率密度函数为

$p(x)={\prod }_{j=1}^{n}p(x_j;{\mu }_j,{\sigma }_j^2)={\prod }_{j=1}^n\frac{1}{\sqrt{2\pi }{\sigma }_j}exp(-\frac{(x_j-{\mu }_j{)}^2}{2{\sigma }_j^2})$
选择一个threshold，即ε，以划定正常与异常的边界。当p(x) >= ε，可认为是正常；当p(x) < ε，可认为是异常。不相关特征的多元高斯分布比较简单，概率密度函数就是各特征的一元概率密度函数乘积，代码实现可参见下节相关特征多元高斯分布代码。

1.4 多个相关特征，且符合多元高斯分布

当两个特征相关时（如下图），红色的是上节中高斯分布模型获得的判定边界，很明显绿色X所代表的数据点很可能是异常值，但是其P(x)值却仍然在正常范围内。
对于不相关特征的多元高斯分布模型，我们计算P(x)的方法是：通过分别计算每个特征对应的几率然后将其累乘起来。对于相关特征的多元高斯分布模型，我们需要考虑各特征之间的关系，因此构建特征的协方差矩阵，用所有的特征一起来计算P(x)，获得图中蓝色判定边界。

$\mu =\frac{1}{m}{\sum }_{i=1}^m{x}^{(i)}$

$\sum =\frac{1}{m}{\sum }_{i=1}^m(x^{(i)}-\mu )(x^{(i)}-\mu {)}^T$

$p(x)=\frac{1}{(2\pi {)}^{\frac{n}{2}}|\Sigma {|}^{\frac{1}{2}}}\exp \left(-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )\right)$

代码实现

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)

"""生成二维数据"""
data = np.random.multivariate_normal(mean=[1, 3],
                                      cov=np.array([[1, 0.5], [0.5, 0.8]]),
                                      size=500)

plt.scatter(data[:,0], data[:,1])
plt.xlabel('x1')
plt.ylabel('x2')
plt.show()

"""计算概率密度函数"""
def getGaussianParams(X, useMultivariate):
    mu = X.mean(axis=0)
    if useMultivariate: # 若为多元相关
        sigma2 = ((X-mu).T @ (X-mu)) / len(X)
    else:
        sigma2 = X.var(axis=0,ddof=0) # 若为一元或多元非相关
    return mu,sigma2

def gaussian(X,mu,sigma2):
    m,n = X.shape
    if np.ndim(sigma2) == 1: # 获取数组的维度，若为多元非相关，只有对角线有自相关系数，维度为1
        sigma2 = np.diag(sigma2) #形成一个以一维数组为对角线元素的矩阵
    
    norm = 1 / (np.power((2*np.pi),n/2) * np.sqrt(np.linalg.det(sigma2)))
    exp = np.zeros((m,1))
    for row in range(m):
        xrow = X[row]
        exp[row] = np.exp(-0.5*((xrow-mu).T).dot(np.linalg.inv(sigma2)).dot(xrow-mu))
    return norm*exp

mu, sigma = getGaussianParams(data, useMultivariate=True)

p = gaussian(data, mu, sigma)

"""multi_gussian检测异常值"""
# 假设阈值为0.005（有监督学习可以根据标签确定阈值）
error_local = np.where(p<0.005)
error = data[error_local[0]]
print("异常数据为：{}".format(error))
    
plt.scatter(error[:,0], error[:,1])
plt.xlabel('x1')
plt.ylabel('x2')
plt.show()

异常数据为：[[-0.50833363  0.58886536]
 [ 1.7314507   1.36040277]
 [ 2.04720929  5.52586327]
 [ 1.21888318  0.96678171]
 [ 1.25947239  1.06794516]
 [-1.44656498  0.93459232]
 [ 1.02330845  0.95222086]
 [ 2.56651417  1.82406223]
 [ 2.72521972  2.125056  ]
 [-0.40398886  0.03560377]
 [ 3.8286177   4.38694095]
 [-1.2223719   0.85178016]
 [-0.53442685  3.97062306]
 [-0.64831769  0.61560726]]

1.5 使用混合参数分布

在许多情况下假定数据是由正态分布产生的。实际数据不符合高斯分布时，可以假定数据是被混合参数分布产生的。

2、非参数方法

在非参数方法中，模型从输入数据学习，而不是假定一个先验。

例如：使用直方图方法检测异常点。

2.1 HBOS

HBOS全名为：Histogram-based Outlier Score。它是一种单变量方法的组合，不能对特征之间的依赖关系进行建模，但是计算速度较快，对大数据集友好。其基本假设是数据集的每个维度相互独立。然后对每个维度进行区间(bin)划分，区间的密度越高，异常评分越低。

HBOS算法流程：

1. 为每个数据维度做出数据直方图。对分类数据统计每个值的频数并计算相对频率。

2. 对直方图进行归一化处理。

3. 计算每一个实例的HBOS值：

$$HBOS(p)=\sum _{i=0}^d\log \left(\frac{1}{{\text{hist}}_i(p)}\right)$$

推导过程：

假设样本p第 i 个特征的概率密度为$p_i(p)$ ，则p的概率密度可以计算为：
$$P(p)=P_1(p)P_2(p)\cdots P_d(p)$$
两边取对数：
$$\begin{array}{rl}\log (P(p))& =\log \left(P_1(p)P_2(p)\cdots P_d(p)\right)=\sum _{i=1}^d\log \left(P_i(p)\right)\end{array}$$
概率密度越大，异常评分越小，为了方便评分，两边乘以“-1”：
$$-\log (P(p))=-1\sum _{i=1}^d\log \left(P_t(p)\right)=\sum _{i=1}^d\frac{1}{\log \left(P_i(p)\right)}$$
最后可得：
$$HBOS(p)=-\log (P(p))=\sum _{i=1}^d\frac{1}{\log \left(P_i(p)\right)}$$
代码实现

# 导入相关依赖模块
from pyod.utils.data import evaluate_print,generate_data
from pyod.models.hbos import HBOS
from pyod.utils.example import visualize

contamination = 0.05  # percentage of outliers
n_train = 1000  # number of training points
n_test = 500  # number of testing points
X_train, y_train, X_test, y_test = generate_data(n_train=n_train, n_test=n_test, contamination=contamination)

# 初始化HBOS模型
clf_name = 'HBOS'
clf = HBOS() 
clf.fit(X_train)

# get the prediction labels and outlier scores of the training data
y_train_pred = clf.labels_  # binary labels (0: inliers, 1: outliers)
y_train_scores = clf.decision_scores_  # raw outlier scores

# get the prediction on the test data
y_test_pred = clf.predict(X_test)  # outlier labels (0 or 1)
y_test_scores = clf.decision_function(X_test)  # outlier scores

# evaluate and print the results
print("\nOn Training Data:")
evaluate_print(clf_name, y_train, y_train_scores)
print("\nOn Test Data:")
evaluate_print(clf_name, y_test, y_test_scores)

# 可视化展示训练集、测试集异常检测结果
visualize(clf_name, X_train, y_train, X_test, y_test, y_train_pred,
          y_test_pred, show_figure=False, save_figure=True)

3、总结

• 使用统计方法做异常检测时，数据需满足所使用方法的统计假定。

• HBOS是一种单变量方法的组合，有点类似参数方法中多个不相关特征的多元高斯分布方法。它不能对特征之间的依赖关系进行建模，但是计算速度较快，对大数据集友好。

参考资料

[1] https://github.com/datawhalechina/team-learning-data-mining/tree/master/AnomalyDetection
[2] http://cs229.stanford.edu/
[3] https://zhuanlan.zhihu.com/p/354872685
[4] https://blog.csdn.net/Sirow/article/details/112692357
[5] https://blog.csdn.net/weixin_42032907/article/details/112726848