signed

QiShunwang

“诚信为本、客户至上”

异常检测(二)

2021/5/14 23:22:35   来源:

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

主要内容包括:

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

概述

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

异常检测的统计学方法的一般思想是:学习一个拟合给定数据集的生成模型,然后识别该模型低概率区域中的对象,把它们作为异常点。

利用统计学方法建立一个模型,然后考虑对象有多大可能符合该模型

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

参数方法假定正常的数据对象被一个以 Θ \Theta Θ为参数的参数分布产生。该参数分布的概率密度函数 f ( x , Θ ) f(x,\Theta) f(x,Θ)给出对象 x x x被该分布产生的概率。该值越小, x 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 ) ∼ N ( μ , σ 2 ) x^{(i)}\sim N(\mu, \sigma^2) x(i)N(μ,σ2),我们可以根据样本求出参数 μ \mu μ σ \sigma σ

μ = 1 m ∑ i = 1 m x ( i ) \mu=\frac 1m\sum_{i=1}^m x^{(i)} μ=m1i=1mx(i)

σ 2 = 1 m ∑ i = 1 m ( x ( i ) − μ ) 2 \sigma^2=\frac 1m\sum_{i=1}^m (x^{(i)}-\mu)^2 σ2=m1i=1m(x(i)μ)2

求出参数之后,我们就可以根据概率密度函数计算数据点服从该分布的概率。正态分布的概率密度函数为

p ( x ) = 1 2 π σ e x p ( − ( x − μ ) 2 2 σ 2 ) p(x)=\frac 1{\sqrt{2\pi}\sigma}exp(-\frac{(x-\mu)^2}{2\sigma^2}) p(x)=2π σ1exp(2σ2(xμ)2)

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

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

例如常用的3sigma原则中,如果数据点超过范围 ( μ − 3 σ , μ + 3 σ ) (\mu-3\sigma, \mu+3\sigma) (μ3σ,μ+3σ),那么这些点很有可能是异常点。
在这里插入图片描述

代码实现

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的那些数据。
在这里插入图片描述
代码实现

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 j j维:

μ j = 1 m ∑ i = 1 m x j ( i ) \mu_j=\frac 1m\sum_{i=1}^m x_j^{(i)} μj=m1i=1mxj(i)

σ j 2 = 1 m ∑ i = 1 m ( x j ( i ) − μ j ) 2 \sigma_j^2=\frac 1m\sum_{i=1}^m (x_j^{(i)}-\mu_j)^2 σj2=m1i=1m(xj(i)μj)2

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

p ( x ) = ∏ j = 1 n p ( x j ; μ j , σ j 2 ) = ∏ j = 1 n 1 2 π σ j e x p ( − ( x j − μ j ) 2 2 σ 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}) p(x)=j=1np(xj;μj,σj2)=j=1n2π σj1exp(2σj2(xjμj)2)
选择一个threshold,即ε,以划定正常与异常的边界。当p(x) >= ε,可认为是正常;当p(x) < ε,可认为是异常。不相关特征的多元高斯分布比较简单,概率密度函数就是各特征的一元概率密度函数乘积,因此省略代码实现。
在这里插入图片描述
对于阈值ε如何选择,有监督学习方法可以采用标签数据进行学习,以得到阈值。
在这里插入图片描述
在这里插入图片描述

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

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

μ = 1 m ∑ i = 1 m x ( i ) \mu=\frac{1}{m}\sum^m_{i=1}x^{(i)} μ=m1i=1mx(i)

∑ = 1 m ∑ i = 1 m ( x ( i ) − μ ) ( x ( i ) − μ ) T \sum=\frac{1}{m}\sum^m_{i=1}(x^{(i)}-\mu)(x^{(i)}-\mu)^T =m1i=1m(x(i)μ)(x(i)μ)T

p ( x ) = 1 ( 2 π ) n 2 ∣ Σ ∣ 1 2 exp ⁡ ( − 1 2 ( x − μ ) T Σ − 1 ( x − μ ) ) 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) p(x)=(2π)2nΣ211exp(21(xμ)TΣ1(xμ))

ps:当多元高斯分布模型的协方差矩阵 ∑ \sum 为对角矩阵,且对角线上的元素为各自一元高斯分布模型的方差时,二者是等价的。
代码实现

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 multi_gussian(x,mu,sigma):
  p=np.zeros((len(x),1))
  n=len(mu)
  if np.ndim(sigma)==1:
    sigma=np.diag(sigma)
    for i in range(len(x)):
        p[i]=(2*np.pi)**(-n/2)*np.linalg.det(sigma)**(-1/2)*np.exp(-0.5*(x[i,:]-mu).T@np.linalg.inv(sigma)@(x[i,:]-mu))   
    return p

mu=data.mean(axis=0)
sigma=data.std(axis=0)
p=multi_gussian(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()

在这里插入图片描述

异常数据为:[[ 3.08549385  5.23758135]
 [-0.50833363  0.58886536]
 [-1.5623398   1.66217866]
 [-1.51564778  2.01708909]
 [-1.05790246  1.27879517]
 [ 2.04720929  5.52586327]
 [-0.6958563   0.75362435]
 [ 3.20172631  5.29466879]
 [ 3.47608451  4.20409127]
 [ 3.00196877  4.84432947]
 [-1.44656498  0.93459232]
 [ 2.83375548  5.20115891]
 [ 2.89755716  4.87988951]
 [ 3.16287674  4.78955739]
 [-0.30589789  0.62841779]
 [-0.40398886  0.03560377]
 [ 3.8286177   4.38694095]
 [-1.2223719   0.85178016]
 [-0.64831769  0.61560726]]

1.5 使用混合参数分布

在许多情况下假定数据是由正态分布产生的。当实际数据很复杂时,这种假定过于简单,可以假定数据是被混合参数分布产生的。

2、非参数方法

在异常检测的非参数方法中,“正常数据”的模型从输入数据学习,而不是假定一个先验。通常,非参数方法对数据做较少假定,因而在更多情况下都可以使用。

例子:使用直方图检测异常点。

直方图是一种频繁使用的非参数统计模型,可以用来检测异常点。该过程包括如下两步:

步骤1:构造直方图。使用输入数据(训练数据)构造一个直方图。该直方图可以是一元的,或者多元的(如果输入数据是多维的)。

尽管非参数方法并不假定任何先验统计模型,但是通常确实要求用户提供参数,以便由数据学习。例如,用户必须指定直方图的类型(等宽的或等深的)和其他参数(直方图中的箱数或每个箱的大小等)。与参数方法不同,这些参数并不指定数据分布的类型。

步骤2:检测异常点。为了确定一个对象是否是异常点,可以对照直方图检查它。在最简单的方法中,如果该对象落入直方图的一个箱中,则该对象被看作正常的,否则被认为是异常点。

对于更复杂的方法,可以使用直方图赋予每个对象一个异常点得分。例如令对象的异常点得分为该对象落入的箱的容积的倒数。

使用直方图作为异常点检测的非参数模型的一个缺点是,很难选择一个合适的箱尺寸。一方面,如果箱尺寸太小,则许多正常对象都会落入空的或稀疏的箱中,因而被误识别为异常点。另一方面,如果箱尺寸太大,则异常点对象可能渗入某些频繁的箱中,因而“假扮”成正常的。

2.1 HBOS

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

HBOS算法流程:

1.为每个数据维度做出数据直方图。对分类数据统计每个值的频数并计算相对频率。对数值数据根据分布的不同采用以下两种方法:

  • 静态宽度直方图:标准的直方图构建方法,在值范围内使用k个等宽箱。样本落入每个桶的频率(相对数量)作为密度(箱子高度)的估计。时间复杂度: O ( n ) O(n) O(n)

  • 2.动态宽度直方图:首先对所有值进行排序,然后固定数量的 N k \frac{N}{k} kN个连续值装进一个箱里,其中N是总实例数,k是箱个数;直方图中的箱面积表示实例数。因为箱的宽度是由箱中第一个值和最后一个值决定的,所有箱的面积都一样,因此每一个箱的高度都是可计算的。这意味着跨度大的箱的高度低,即密度小,只有一种情况例外,超过k个数相等,此时允许在同一个箱里超过 N k \frac{N}{k} kN值。

    时间复杂度: O ( n × l o g ( n ) ) O(n\times log(n)) O(n×log(n))

2.对每个维度都计算了一个独立的直方图,其中每个箱子的高度表示密度的估计。然后为了使得最大高度为1(确保了每个特征与异常值得分的权重相等),对直方图进行归一化处理。最后,每一个实例的HBOS值由以下公式计算:

H B O S ( p ) = ∑ i = 0 d log ⁡ ( 1 hist i ( p ) ) H B O S(p)=\sum_{i=0}^{d} \log \left(\frac{1}{\text {hist}_{i}(p)}\right) HBOS(p)=i=0dlog(histi(p)1)

推导过程

假设样本pi 个特征的概率密度为 p i ( p ) p_i(p) pi(p) ,则p的概率密度可以计算为:
P ( p ) = P 1 ( p ) P 2 ( p ) ⋯ P d ( p ) P(p)=P_{1}(p) P_{2}(p) \cdots P_{d}(p) P(p)=P1(p)P2(p)Pd(p)
两边取对数:
log ⁡ ( P ( p ) ) = log ⁡ ( P 1 ( p ) P 2 ( p ) ⋯ P d ( p ) ) = ∑ i = 1 d log ⁡ ( P i ( p ) ) \begin{aligned} \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{aligned} log(P(p))=log(P1(p)P2(p)Pd(p))=i=1dlog(Pi(p))
概率密度越大,异常评分越小,为了方便评分,两边乘以“-1”:
− log ⁡ ( P ( p ) ) = − 1 ∑ i = 1 d log ⁡ ( P t ( p ) ) = ∑ i = 1 d 1 log ⁡ ( P i ( p ) ) -\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)} log(P(p))=1i=1dlog(Pt(p))=i=1dlog(Pi(p))1
最后可得:
H B O S ( p ) = − log ⁡ ( P ( p ) ) = ∑ i = 1 d 1 log ⁡ ( P i ( p ) ) H B O S(p)=-\log (P(p))=\sum_{i=1}^{d} \frac{1}{\log \left(P_{i}(p)\right)} HBOS(p)=log(P(p))=i=1dlog(Pi(p))1
代码实现

# 导入相关依赖模块
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] http://cs229.stanford.edu/
[2] https://blog.csdn.net/Sirow/article/details/112692357