【Datawhale第25期组队学习】Task02：基于统计学的异常检测算法

Task02：基于统计学的异常检测算法

0 写在前面

上篇文章，我们从异常检测的概念、异常检测方法、异常检测应用、异常检测面临的挑战和未来研究方向等方面对异常检测进行了全方位的介绍。

本篇文章，我们将视角回到异常检测算法上来，下面将详细介绍基于统计学的异常检测算法.

1 基于统计学的异常检测算法概述

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

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

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

根据如何指定和学习模型，异常检测的统计学方法可以划分为两个主要类型：参数方法和非参数方法。这两种方法的主要区别在于，前者对给定数据中的潜在分布模型有一个假设，它根据已知数据估计分布模型的参数。后一种方法没有任何假设。

优点：①它们在数学上是可以接受的，一旦模型建立，就有一个快速的评估过程。②模型通常适合定量实值数据集或一些定量有序数据分布。序数数据可以更改为适当的值进行处理，从而缩短复杂数据的处理时间。

缺点：①由于它们的依赖性和参数模型中分布模型的假设，由于缺乏关于潜在分布的先验知识，产生的结果的质量在实际情况和应用中大多是不可靠的。②于大多数模型适用于单变量特征空间，因此它们在多维场景中不适用。在处理多变量数据时，它们会产生很高的计算成本。

2 参数化方法

参数方法假定正常的数据对象被一个以$\Theta$为参数的参数分布产生。该参数分布的概率密度函数$f(x,\Theta )$给出对象$x$被该分布产生的概率。该值越小，$x$越可能是异常点。
这类方法通常会假设给定的数据集服从一个随机分布模型，将与模型不一致的样本视为异常样本。其中最常用的两种分布模型是一元正态分布模型和多元正态分布模型。

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

仅涉及一个属性或变量的数据称为一元数据。我们假定数据由正态分布产生，然后可以由输入数据学习正态分布的参数，并把低概率的点识别为异常点。

算法模型：

在正态分布的假设下，如果有一个新样本，当x的正态分布值小于某个阈值时，就可以任务这个样本是异常的。

在正态分布中，存在3sigma原则，μ-3σ<=x<=μ+3σ的区域包含了绝大部分数据，可以依次为参考，调整ε的值：

现在有一个包含m个一维数据的训练集：

$X=${$x^{(1)},x^{(2)},...,x^{(m)}$}

这里$x^{(i)}$是已知的，μ和σ才是未知的，我们的目的是设法根据训练集求得μ和σ的值，以得到一个确定的函数模型。具体来说，通过最大似然估计量可以得出下面的结果：

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

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

求出参数之后，我们就可以根据概率密度函数计算数据点服从该分布的概率。可以此函数判断一样样本是否是异常的：

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

下面通过数据模拟来看一下基于正态分布的一元异常点检测。

算法实现：

##引用 https://www.cnblogs.com/bigmonkey/p/11431468.html
import numpy as np
import matplotlib.pyplot as plt

def create_data():
    '''
    创建训练数据和测试数据
    :return: X_train:训练集， X_test:测试集
    '''
    np.random.seed(42)  # 设置seed使每次生成的随机数都相等
    m, s = 3, 0.1 # 设置均值和方差
    X_train = np.random.normal(m, s, 100) # 100个一元正态分布数据
    # 构造10测试数据，从一个均匀分布[low,high)中随机采样
    X_test = np.random.uniform(low=m - 1, high=m + 1, size=10)
    return X_train, X_test

def plot_data(X_train, X_test):
    '''
    数据可视化
    :param X_train: 训练集
    :param X_test: 测试集
    :return:
    '''
    fig = plt.figure(figsize=(10, 4))
    plt.subplots_adjust(wspace=0.5)  # 调整子图之间的左右边距
    fig.add_subplot(1, 2, 1)  # 绘制训练数据的分布
    plt.scatter(X_train, [0] * len(X_train), color='blue', marker='x', label='训练数据')
    plt.title('训练数据的分布情况')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend(loc='upper left')

    fig.add_subplot(1, 2, 2)  # 绘制整体数据的分布
    plt.scatter(X_train, [0] * len(X_train), color='blue', marker='x', label='训练数据')
    plt.scatter(X_test, [0] * len(X_test), color='red', marker='^',label='测试数据')
    plt.title('整体数据的分布情况')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend(loc='upper left')

    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 解决中文下的坐标轴负号显示问题
    plt.show()

def fit(X_train):
    '''
    拟合数据，训练模型
    :param X_train: 训练集
    :return:  mu:均值, sigma:方差
    '''
    global mu, sigma
    mu = np.mean(X_train)  # 计算均值μ
    sigma = np.var(X_train) # 计算方差 σ^2

def gaussian(X):
    '''
    计算正态分布
    :param X: 数据集
    :return: 数据集的密度值
    '''
    return np.exp(-((X - mu) ** 2) / (2 * sigma)) / (np.sqrt(2 * np.pi) * np.sqrt(sigma))

def get_epsilon(n=3):
    ''' 调整ε的值，默认ε=3σ '''
    return np.sqrt(sigma) * n

def predict(X):
    '''
    检测训练集中的数据是否是正常数据
    :param X: 待预测的数据
    :return: P1:数据的密度值, P2:数据的异常检测结果，True正常，False异常
    '''
    P1 = gaussian(X) # 数据的密度值
    epsilon = get_epsilon()
    P2 = [p > epsilon for p in P1] # 数据的异常检测结果，True正常，False异常
    return P1, P2

def plot_predict(X):
    '''可视化异常检测结果 '''
    epsilon = get_epsilon()
    xs = np.linspace(mu - epsilon, mu + epsilon, 50)
    ys = gaussian(xs)
    plt.plot(xs, ys, c='g', label='拟合曲线')  # 绘制正态分布曲线

    P1, P2 = predict(X)
    normals_idx = [i for i, t in enumerate(P2) if t == True] # 正常数据的索引
    plt.scatter([X[i] for i in normals_idx], [P1[i] for i in normals_idx],
                color='blue', marker='x', label='正常数据')
    outliers_idx = [i for i, t in enumerate(P2) if t == False] # 异常数据的索引
    plt.scatter([X[i] for i in outliers_idx], [P1[i] for i in outliers_idx],
                color='red', marker='^', label='异常数据')
    plt.title('检测结果，共有{}个异常数据'.format(len(outliers_idx)))
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend(loc='upper left')
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 解决中文下的坐标轴负号显示问题
    plt.show()

if __name__ == '__main__':
    mu, sigma = 0, 0 # 模型的均值μ和方差σ^2
    X_train, X_test = create_data()
    plot_data(X_train, X_test)
    fit(X_train)
    print('μ = {}, σ^2 = {}'.format(mu, sigma))
    plot_predict(np.r_[X_train, X_test])

结果：

可以看出，大部分训练数据都集中在正态分布的均值区域。

接下来使用fit()方法对异常检测模型进行训练，得到了模型参数后就可以使用目标函数对数据进行预测。gussian(X)实现了正态分布的密度函数；predict(X)将对X中的所有样本进行检测，并返回X对应的检测结果列表。

一元异常点检测的不足

在面对多维数据时，基于一元正态分布的异常检测可以单独抽取某一维度进行检测，通常也能工作的很好，但这里有一个假设—所有维度都符合正态分布，并且各维度都是独立的，如果两个维度之间存在相关性，那么基于一元正态分布的异常检测就可能会出现很大程度的误判。

2.2 基于正态分布的多元异常点检测

涉及两个或多个属性或变量的数据称为多元数据。许多一元异常点检测方法都可以扩充，用来处理多元数据。其核心思想是把多元异常点检测任务转换成一元异常点检测问题。例如基于正态分布的一元异常点检测扩充到多元情形时，可以求出每一维度的均值和标准差。对于第$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})$

这是在各个维度的特征之间相互独立的情况下。如果特征之间有相关性，就要用到多元高斯分布了。

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

$\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)$

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

图源自网络

3 非参数化方法

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

使用直方图检测异常点

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

步骤1：构造直方图。使用输入数据（训练数据）构造一个直方图。该直方图可以是一元的，或者多元的（如果输入数据是多维的）。

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

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

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

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

使用马氏距离 Mahalanobis distance检测异常点

其中S为二元高斯分布的协方差矩阵，[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-AvTHvjOz-1620980071315)(https://www.zhihu.com/equation?tex=%5Cbar%7Bx%7D+)]为所有数据对象的均值。Mahalanobis距离越大，数据对象偏离均值越大，也就越异常。

4 基于角度的方法

基于角度的方法的主要思想是：数据边界上的数据很可能将整个数据包围在一个较小的角度内，而内部的数据点则可能以不同的角度围绕着他们。如下图所示，其中点A是一个异常点，点B位于数据内部。

5 HBOS

HBOS算法基于多维数据各个维度的独立性假设，对于单个数据维度，先做出数据直方图。对于categroy 值，统计每个值出现的次数，并计算相对频率。对于数值特征，可以用两种方法：

• 静态跨度的柱状图：将值域分成K个等宽的桶，落入每个桶的值的频数作为密度的估计（桶的高度）
• 动态宽度柱状图：先将所有值排序，然后将连续的N/k个值装进一个桶里，其中N是所有的样例数，k是桶的个数，是一个超参数；柱状图的面积对应桶中的样例数。因为桶的宽度是有桶中第一个值和最后一个值决定的，所有桶的面积都一样，所以，每一个桶的高度可以被计算出来。这意味着跨度大的桶的高度低，即密度小，只有一种情况例外，超过k个数相等，此时一个桶里允许超过N/k个值。

之所以HBOS提供上述两种方法是因为真实的数据分布种类繁多，尤其是值域跨度大的时候。我们更倾向于使用第二种方法，尤其是值得分布位置或者长尾的情形。
按照上面的方法，每一个维度的数据都生成一个柱状图，每个桶的高度代表数据的密度，使用归一化操作保证桶最高的高度是1.这样可以保证每个特征的权重相同。最终每一个样本的HBOS值按照下面的公示计算：

HBOS分值越高，样本越异常。

算法演示：

# -*- coding: utf-8 -*-
"""Example of using Histogram- based outlier detection (HBOS) for
outlier detection
"""
# # Author: Yue Zhao <zhaoy@cmu.edu>
# # License: BSD 2 clause


from pyod.models.hbos import HBOS
from pyod.utils.data import generate_data
from pyod.utils.data import evaluate_print
from pyod.utils.example import visualize

if __name__ == "__main__":
    contamination = 0.1  # 离群点比重
    n_train = 200  # 训练集样本数
    n_test = 100  # 测试集样本数

    # 生成采样数据
    X_train, y_train, X_test, y_test = \
        generate_data(n_train=n_train,
                      n_test=n_test,
                      n_features=2,
                      contamination=contamination,
                      random_state=42)

    # 训练HBOS离群点检测器
    clf_name = 'HBOS'
    clf = HBOS()
    clf.fit(X_train)

    # 得到训练集样本点的预测标签 和 outlier得分 得分越高点的很可能就是离群点
    y_train_pred = clf.labels_  # 2分类标签 (0: 正常, 1: 异常)
    y_train_scores = clf.decision_scores_  # 原始数据的outlier scores

    # 测试集结果
    y_test_pred = clf.predict(X_test)  # 预测标签 (0 or 1)
    y_test_scores = clf.decision_function(X_test)  # 测试集数据的outlier scores

    # 评估
    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=True, save_figure=True)

模型优化

HBOS(n_bins=10, alpha=0.1, tol=0.5, contamination=0.1)调整四个参数即可
1、n_bins：分箱的数量
2、alpha：用于防止边缘溢出的正则项
3、tol：用于设置当数据点落在箱子外时的宽容度
4、contamination：用于设置异常点的比例
根据可视化结果，来调整4个输入参数，直到得出更优参数。"""

本模型的参数为

print(HBOS.get_params(clf))

{'alpha': 0.1, 'contamination': 0.1, 'n_bins': 10, 'tol': 0.5}

6 总结

• 异常检测的统计学方法由数据学习模型，以区别正常的数据对象和异常点。使用统计学方法的一个优点是，异常检测可以是统计上无可非议的。当然，仅当对数据所做的统计假定满足实际约束时才为真。

• HBOS在全局异常检测问题上表现良好，但不能检测局部异常值。但是HBOS比标准算法快得多，尤其是在大数据集上。

参考文献

[1]Datawhale开源课程.异常检测.2021年.https://github.com/datawhalechina/team-learning-data-mining/tree/master/AnomalyDetection

[2]呼广跃.异常检测之HBOS算法.2021年.https://zhuanlan.zhihu.com/p/84587517

[3]CSDN.马氏距离与异常点检测.2021年.https://blog.csdn.net/weixin_30405421/article/details/98026796

[4]博客园.异常检测(2)——基于概率统计的异常检测（1）.2021年.https://www.cnblogs.com/bigmonkey/p/11431468.html