机器学习笔记--异常检测

1、基本概念

我们以飞机引擎制造为例，作为质量测试，我们测量了飞机引擎的一些特征变量，比如引擎运转时产生的热量，或者引擎的振动等等。这样就有了一个数据集，绘制成如下图表

给定数据集 $x^{(1)},x^{(2)},\dots ,x^{(m)}$，我们假使数据集是正常的，我们希望知道新的数据 $x_{tset}$ 是不是异常的，即这个测试数据不属于该组数据的几率如何。我们所构建的模型应该能根据该测试数据的位置告诉我们其属于一组数据的可能性 $p(x)$。

上图中，在蓝色圈内的数据属于该组数据的可能性较高，而越是偏远的数据，其属于该组数据的可能性就越低。

这种方法称为密度估计，表达如下：

$$\text{ if }p(x)\{\begin{array}{ll}<\epsilon & \text{ anomaly }\\ >=\epsilon & \text{ normal }\end{array}$$

2、高斯分布

2.1 高斯分布

高斯分布，也称为正态分布。通常如果我们认为变量 x 符合高斯分布 $x\sim N\left(\mu ,{\sigma }^2\right)$ 则其概率密度函数为：

$$p\left(x,\mu ,{\sigma }^2\right)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)$$
其中，

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

注意，机器学习中对于方差我们通常只除以 m 而非统计学中的 (m-1)。在实际使用中，他们的区别甚小，几乎可以忽略不计。

当有多个特征变量时，
$$p(x)=\prod _{j=1}^{n}p\left(x_j;{\mu }_j,{\sigma }_j^2\right)=\prod _{j=1}^n\frac{1}{\sqrt{2\pi }{\sigma }_j}\exp \left(-\frac{{\left(x_j-{\mu }_j\right)}^2}{2{\sigma }_j^2}\right)$$

2.2 多元高斯分布

假使我们有两个相关的特征，而且这两个特征的值域范围比较宽，这种情况下，一般的高斯分布模型可能不能很好地识别异常数据。其原因在于，一般的高斯分布模型尝试的是去同时抓住两个特征的偏差，因此会创造出一个比较大的判定边界。

下图中是两个相关特征，洋红色的线是一般的高斯分布模型获得的判定边界，很明显绿色的 X 所代表的数据点很可能是异常值，但是其 $p(x)$ 值却仍然在正常范围内。多元高斯分布将创建像图中蓝色曲线所示的判定边界。

我们首先计算所有特征的平均值，然后再计算协方差矩阵：

$$\Sigma =\frac{1}{m}\sum _{i=1}^m\left(x^{(i)}-\mu \right){\left(x^{(i)}-\mu \right)}^T=\frac{1}{m}(X-\mu {)}^T(X-\mu )$$
也就是将每一列的方差值求出来后，写成对角矩阵的形式，即为 $\Sigma$ 矩阵。

然后得到多元高斯分布：

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

2.3 二者比较

3、算法流程

1. 根据测试集数据，我们估计特征的平均值和方差并构建 $p(x)$ 函数;

2. 对交叉检验集，我们尝试使用不同的𝜀值作为阀值，并预测数据是否异常，根据 F1 值或者查准率与查全率的比例来选择 𝜀;

3. 选出 𝜀 后，针对测试集进行预测，计算异常检测系统的F1值，或者查准率与查全率之比。

一个技巧：

异常检测假设特征符合高斯分布，如果数据的分布不是高斯分布，异常检测算法也能够工作，但是最好还是将数据转换成高斯分布，例如使用对数函数：$x=\log (x+c)$

4、异常检测与监督学习对比

5、课后习题

# @Time : 2021/11/3 14:50
# @Author : xiao cong
# @Function : 异常检测

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat

data = loadmat("ex8data1.mat")
# print(data)
X = data["X"]
Xval, yval = data["Xval"], data["yval"]  # 注意异常点为 1
print("X.shape = {}, Xval.shape = {}, yval.shape = {}".format(X.shape, Xval.shape, yval.shape))


def plot_data(X):
    plt.figure(figsize=(8, 5))
    plt.scatter(X[:, 0], X[:, 1])
    plt.xlabel("X1")
    plt.ylabel("X2")
    plt.show()


plot_data(X)


# ************************************************************************
# 高斯分布
def estimate_gaussian(X):  # 计算均值和方差
    mu = X.mean(axis=0)
    sigma = X.var(axis=0)
    return mu, sigma


mu, sigma = estimate_gaussian(X)
print("mu = {}, sigma = {}".format(mu, sigma))


def gaussian(X, mu, sigma):  # 定义多元高斯分布函数
    m, n = X.shape
    Sigma = np.diag(sigma)  # 写成对角线形式，以便后续计算
    term1 = 1 / (np.power(2 * np.pi, n / 2) * np.sqrt(np.linalg.det(Sigma)))
    term2 = np.zeros((m, 1))

    for i in range(m):
        term2[i] = np.exp(-0.5 * (X[i] - mu).dot(np.linalg.inv(Sigma)).dot((X[i] - mu).T))

    return term1 * term2


def plot_contours(mu, sigma):  # 画高斯分布轮廓线
    """
    画出高斯概率分布的图，在三维中是一个上凸的曲面。投影到平面上则是一圈圈的等高线。
    """
    delta = 0.3
    x = np.arange(0, 30, delta)
    y = np.arange(0, 30, delta)
    xx, yy = np.meshgrid(x, y)  # 生成网格, xx维度(100, 100)
    points = np.c_[xx.flatten(), yy.flatten()]  # 分别展开，然后合并为两列
    z = gaussian(points, mu, sigma)
    z = z.reshape(xx.shape)

    plt.figure(figsize=(8, 5))
    cont_levels = [10 ** h for h in range(-20, 0, 3)]  # 用来确定轮廓线/区域的数量和位置
    plt.contour(xx, yy, z, cont_levels)
    plt.title('Gaussian Contours')
    plt.scatter(X[:, 0], X[:, 1])
    plt.xlabel("X1")
    plt.ylabel("X2")


plot_contours(mu, sigma)
plt.show()

# **********************************************************************************
# 选择阈值（即找到使F1最大的阈值）
def select_threshold(yval, pval):
    def compute_F1(yval, pval):
        m = len(yval)
        tp = float(len([i for i in range(m) if pval[i] and yval[i]]))
        # print([i for i in range(m) if pval[i] and yval[i]])   # 列表类型，返回符合条件的索引
        fp = float(len([i for i in range(m) if pval[i] and not yval[i]]))
        fn = float(len([i for i in range(m) if not pval[i] and yval[i]]))

        prec = tp / (tp + fp) if (tp + fp) else 0
        rec = tp / (tp + fn) if (tp + fn) else 0          # 防止分母为 0
        F1 = 2 * prec * rec / (prec + rec) if (prec + rec) else 0
        return F1

    epsilons = np.linspace(min(pval), max(pval), 1000)
    bestF1, bestEpsilon = 0, 0
    for eps in epsilons:
        pval_ = (pval < eps)  # 当预测概率小于阈值时，即为异常点(异常点为 1)，说明预测正确
        thisF1 = compute_F1(yval, pval_)
        if thisF1 > bestF1:
            bestF1 = thisF1
            bestEpsilon = eps

    return bestF1, bestEpsilon


pval = gaussian(Xval, mu, sigma)
bestF1, bestEpsilon = select_threshold(yval, pval)
print("bestF1 = {}, bestEpsilon = {}".format(bestF1, bestEpsilon))

# 在训练集中找出离群点
predictions = gaussian(X, mu, sigma)
xx = np.array([X[i] for i in range(len(predictions)) if predictions[i] < bestEpsilon])
# 注意要将得到的列表转换为数组
plot_contours(mu, sigma)
plt.scatter(xx[:, 0], xx[:, 1], c='r')
plt.show()


# ***************************************************
# 高维数据
data2 = loadmat("ex8data2.mat")
# print(data2)
X2 = data2["X"]
Xval2, yval2 = data2["Xval"], data2["yval"]
print("X2.shape = {}, Xval2.shape = {}, yval2.shape = {}".format(X2.shape, Xval2.shape, yval2.shape))

mu2, sigma2 = estimate_gaussian(X2)
pval2 = gaussian(Xval2, mu2, sigma2)
bestF12, bestEpsilon2 = select_threshold(yval2, pval2)

predictions2 = gaussian(X2, mu2, sigma2)
anoms = [X2[i] for i in range(X2.shape[0]) if predictions2[i] < bestEpsilon2]
print("bestEpsilon2 = {}, len(anoms) = {}".format(bestEpsilon2, len(anoms)))