8.吴恩达机器学习-异常检测

一、理论原理

1. 异常检测

  异常检测(Anomaly detection)问题是机器学习算法的一个常见应用。这种算法的一个有趣之处在于：它虽然主要用于非监督学习问题，但从某些角度看，它又类似于一些监督学习问题。

什么是异常检测呢？
  飞机引擎制造商生产的飞机引擎，从生产线上流出时，需要进行 QA(质量控制测试)，而作为这个测试的一部分，你测量了飞机引擎的一些特征变量，比如，你可能测量了 引擎运转时产生的热量 或者 引擎的振动 等等。这样，你就有了一个数据集 $x^{(1)},x^{(2)},\dots ,x^{(m)}$，如果你生产了 m 个引擎的话，也许你会将这些数据绘制成图表，看起来就是这个样子，这里的每个点、每个叉都是你的无标签数据。

  这样异常检测问题可以定义为对于生产线上流出的一个新的飞机引擎，有特征变量 $x_{test}$，那么这个新的飞机引擎是否有某种异常，如果你的新引擎对应的点落在中间的那个绿叉，我们可以直接认为它是正常的。如果你的新飞机引擎对应的点落在外面那个绿叉，那么我们可以认为这是一个异常需要进一步检测。

异常检测问题：我们有一些数据 $x^{(1)},x^{(2)},\dots ,x^{(m)}$，我们通常假定这 m 个样本都是正常的，或者说都不是异常的。然后我们需要一个算法来告诉我们一个新的样本数据 $x_{test}$ 是否是异常的，我们要采取的方法是给定无标签的训练集，对数据建立一个模型 $p(x)$，也就是说我们将对 $x$ 的分布概率建模，其中 $x$ 是这些特征变量，建立了 $x$ 的概率模型之后，对于新的 $x_{test}$，如果概率 $p$ 低于阈值 $ϵ$，那么就将其标记为异常。

  给定图中的这个训练集，如果你建立了一个模型 $p(x)$，将会认为在中心区域的这些点有很大的概率值，而稍微远离中心区域的点概率会小一些。

2. 高斯（正态）分布

  如果 $x$ 的概率分布服从高斯分布，其中均值为 $\mu$，方差为 ${\sigma }^2$，$\sigma$ 也叫标准差。为了表示高斯分布，大写字母 $N$ 表示 $Normal$ (正态)，因为高斯分布就是正态分布。其中 $\mu$ 控制曲线的中心位置，$\sigma$ 控制曲线的宽度，因为这是一个概率分布。因此曲线下的面积，这些阴影区域的积分一定是1，所以当 $\sigma$ 越小的时候，曲线的峰值就越高。

让我们来看参数估计问题：
  这里有一个数据集 $x^{(1)},x^{(2)},\dots ,x^{(m)}$，假设它们都是实数，我们猜测这些样本来自同一个高斯分布的总体。那么我们就要确定参数 $\mu$ 和 ${\sigma }^2$，其中对两参数进行参数估计的标准公式为：$$\mu =\frac{1}{m}\sum _{i=1}^m{x}^{(i)},{\sigma }^2=\frac{1}{m}\sum _{i=1}^m(x^{(i)}-\mu {)}^2$$这里的估计实际就是对 $\mu$ 和 ${\sigma }^2$ 的极大似然估计。

3. 高斯分布开发异常检测算法

  假如说我们有一个无标签的训练集 $x^{(1)},x^{(2)},\dots ,x^{(m)}$，并且训练集里的每一个样本都是 n 维的特征。现在，我们解决异常检测的方法是：我们要从数据中建立一个 $p(x)$ 概率模型，如图所示，这个式子实际上对应于一个从 $x_1$ 到 $x_n$ 上的独立的假设，但实际中无论这些特征是否独立，这个算法的效果也还不错。

估计 $p(x)$ 的分布问题被称为密度估计问题，把所有的结合起来，下面便是我们的异常检测算法：
  ①选择特征，即找出我们认为的具有比较反常的样本特征 $x_i$
  ②给出一个由 m 个无标签数据（n 维）构成的训练集，计算出期望 ${\mu }_1,{\mu }_2,\dots ,{\mu }_n$以及方差值 ${\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_n^2$，其中$${\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)$ 的值，如果 $p(x_{test})<ϵ$（$ϵ$ 为我们设定的阈值），那么你就可以将这一项标注为异常。$$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}{e}^{-\frac{(x_j-{\mu }_j{)}^2}{2{\sigma }_j^2}}$$
以下是一个异常检测的实例，假如说我们有一系列数据，如左上角绘制的数据集，并估算出两个特征值 $x_1$ 和 $x_2$ 的高斯分布公式：左下角绘制 $p(x)$ 的图像，也就是这两个概率值的乘积。然后我们选取选 $ϵ$ (后面会讲到这个值如何选取)的值为 0.02，那么当给出一个新样本时，我们就能计算出 $p(x_{test}^{(1)})$，来判断它是否异常了。

4. 建立一个异常检测系统

  为了更快地开发出一个异常检测系统，那么最好能找到某种评价异常检测系统的方法。

  为了做到这一点，我们先假定已有了一些带标签的数据。到目前为止我们考虑的异常检测问题是一个非监督问题，使用的是无标签数据。但如果你有一些带标签的数据，能够指明哪些是异常样本，哪些是正常样本，那么这就是我们要找的能够评价异常检测算法的标准方法。

  以飞机引擎为例，假如我们现在有了一些带标签数据，也就是有一些异常样本。同时我们还有一些无异常的样本，我用 y=0 来表示那些完全没有问题的样本，用 y=1 来代表那些异常样本。
那么对于异常检测算法的推导和评价方法，我们先只考虑训练样本，将它们看成无标签的训练集。所以它们都是无异常样本的集合（可能有一些异常的，这也没关系），接下来我们要定义交叉验证集和测试集，通过这两个集合我们将得到异常检测算法。具体来说，我们将假设交叉验证集和测试集中有一些异常的样本。

  假如说：我们有 10000 个无异常的引擎，假设这 10000 个样本中大多数都是好的、没有问题的引擎。从这 10000 个正常的引擎中挑选 6000 个放到无标签的训练集中，我叫它“无标签训练集”。但其实所有这些样本实际上都对应 y=0的情况，我们要用它们来拟合 $p(x)$，参数是 ${\mu }_j,{\sigma }_j^2$，因此我们就是要用这 6000 个样本来估计参数 ${\mu }_1,{\mu }_1,\dots ,{\mu }_n$ 和 ${\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_n^2$。这就是训练集中的好样本，或者说大多数好的样本。然后，我们把剩下的 4000 个样本各挑选 2000 个为交叉验证集和测试集中。同时，我们还有 20 个异常样本，同样也把它们进行一个分割，放 10 个到验证集中，剩下 10 个放入测试集中。

异常检测算法的推导如下：首先我们使用训练样本来拟合模型 $p(x)$，当从交叉验证集和测试集中随机挑选某个测试样本 x 时，假设这个算法对 $p(x)<ϵ$ 的情况做出的预测为 y=1，而 $p(x)\ge ϵ$ 时，算法做出的预测为 y=0。

用什么评价度量好呢？
  因为数据是非常偏斜的，即 y=0 是更常见的，因此分类准确度不是一个好的度量法，我们应该用真阳性、假阳性、假阴性和真阴性的此例来作为评价度量值，也可以算出查准率和召回率或者算出 F1分数，通过一个简单的数字来总结出查准和召回的大小。通过这些方法，你就可以评价你的异常检测算法在交叉验证和测试集样本中的表现。

  当我们需要作出决定时，比如要包括哪些特征或确定参数 $ϵ$ 取多大合适时，我们就可以不断地用交叉验证集来评价这个算法，然后决定我们应该用哪些特征，怎样选择 $ϵ$。

5. 异常检测和监督学习

  既然我们有了这些带标签的数据，那么我们为什么不直接用监督学习的方法呢？为什么不直接用逻辑回归或者神经网络的方法来直接学习这些带标签的数据，从而给出预测 y=1 或 y=0 呢？

下面说明了什么时候用异常检测，什么时候用监督学习。

6. 特征值选择

  在我们的异常检测算法中，我们做的事情之一就是使用高斯分布来对特征向量建模：

  首先画出这些数据，可以用直方图表示数据以确保这些数据在应用异常检测算法前看起来像高斯分布，当然即使你的数据并不是高斯分布，它也可以良好地运行，画柱状图的方法是 plt.hist() 命令。如果我画出来的直方图是左下角这样的话，可以对数据进行一些不同的转换(例如取对数) ，来确保这些数据看起来更像高斯分布，虽然通常来说你不这么做，算法也会运行地很好。但如果你使用一些转换方法，使你的数据更像高斯分布的话，你的算法会工作得更好，除了取对数变换之外，还有别的一些方法也可以用如上图所示。

如何得到异常检测算法的特征变量？
  我们先完整地训练出一个学习算法，然后在一组交叉验证集上运行算法，然后找出那些预测出错的样本，再看看我们能否找到一些其他的特征变量来帮助学习算法，让它在那些交叉验证时判断出错的样本中表现更好。
让我们通过一个例子来解释一下上述过程：
  在异常检测中，我们希望：$p(x)$ 的值对正常样本来说是比较大的，而对异常样本来说，值是很小的。假如数据拟合出的高斯分布是这样的，异常样本中 x 的取值为2.5(绿色)，因此我画出异常样本，它看起来就像被淹没在一堆正常样本中似的。它的概率值很大（即高斯分布中的 $p(x)$ 值很大），由于概率值大，所以我们的算法没能把这个样本判断为异常。

  现在如果说这代表飞机引擎的制造或者别的什么，那么我会看看我的训练样本，看到底是哪一个具体的飞机引擎出错了，通过这个样本能不能启发我想出一个新的特征 $x_2$ 来帮助算法区别出不好的样本和剩下的正确样本(也就是左图中那些红叉)。绿叉的那个样本，如果用特征 $x_1$ 检测，则为正常，而如果用特征 $x_2$ 检测，则为异常，因此我们需要找出这个特征变量 $x_2$。
  通常来说，我想到的选择特征变量的方法是：选那些取值既不会特别特别大、也不会特别特别小的特征变量，比如说数据中心中监控计算机的例子，你有很多台电脑 ，也许上千或者上万台，我们想要知道的是是不是有哪一台机器运作不正常了，这里给出了几种可选的特征变量：包括占用内存、磁盘每秒访问次数、CPU负载、网络流量。现在假如说：我认为在我的数据中，CPU负载和网络流量应该互为线性关系，假如说我怀疑其中一个出错的情形是计算机在执行一个任务时进入了一个死循环，导致其被卡住了，因此CPU负载升高，在这种情况下，要检测出异常，我可以新建一个特征 $x_5$ 等于 CPU负载除以网络流量。
你可以通过不同特征变量的组合捕捉到对应的不寻常现象。

7. 多元高斯分布

  假设我们的无标签数据看起来像这张图一样。使用数据中心的监控机的例子，两个特征变量：$x_1$ 是 CPU 的负载， $x_2$ 是内存使用量。 把 $x_1$ 和 $x_2$ 当做高斯分布来建模，如右图所示，这就是异常检测算法的特征变量建模。测试集中，正常的样本大致如下图左图所示，数据大致呈现一条直线分布， CPU 负载和内存使用量大致呈彼此线性增长的关系。

  现在假如说，有一个样本(下图中左图的绿叉)，它离这里看到的任何数据都很远(下图中左图的红叉)，它应该被当做一个异常数据，它的 CPU 负载和内存使用量不是呈彼此线性增长的关系，它的CPU 负载低，内存使用量高。异常检测算法的做法如图右所示，绿叉的 $x_1$ 和 $x_2$ 都属于比较正常的范畴。所以，异常检测算法不会将这个点标记为异常，它不能察觉到这个蓝色椭圆所表示的好样本概率高的范围，而只能看到图左红色圆的范围。

为了解决这个问题，我们使用多元高斯分布或者多元正态分布。
我们有特征 $x$ ，它是 $n$ 维实数，不要把 $p(x_1),p(x_2),\dots$ 分开，而要建立一个 $p(x)$ 整体的模型，就是一次性建立模型（$|\Sigma |$ 表示 $\Sigma$ 的行列式） $$p(x;\mu ,\Sigma )=\frac{1}{(2\pi {)}^{\frac{n}{2}}|\Sigma {|}^{\frac{1}{2}}}{e}^{-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )}$$ 多元高斯分布的参数包括向量 $\mu$ 和一个 $n\times n$ 的协方差矩阵 $\Sigma$。

我们来看一个二维的例子，如果 n=2，两个特征 $x_1$ 和 $x_2$，以下是 $\mu$ 和 $\Sigma$ 发生变化时对高斯分布的影响。$\mu$ 的值可以改变分布的峰值，$\Sigma$ 的值可以改变分布的概率变化范围和高度（积分要为1）。


多元高斯分布的一个优点是你可以用它给数据的相关性建立模型，如果你改变协方差矩阵非对角线上的元素，你会得到一种不同的高斯分布。所以当我将非对角线的元素从 0.5 增加到 0.8 时，我会得到一个沿着 y=x 这条线更窄、更高的分布。当两个特征变量之间可能存在正相关或负相关关系的情况，多元高斯分布可以让你能够描述出来。

8. 利用多元高斯分布进行异常检测

检测步骤：
① 设置参数 $\mu$ 和 $\Sigma$，用数据集 { x ( 1 ) , x ( 2 ) , … , x ( m ) } \{x^{(1)},x^{(2)},\dots,x^{(m)}\} {x(1),x(2),…,x(m)} 拟合模型 $p(x)$：$$\mu =\frac{1}{m}\sum _{i=1}^m{x}^{(i)},\Sigma =\frac{1}{m}\sum _{i=1}^m(x^{(i)}-\mu )(x^{(i)}-\mu {)}^T$$
② 给一个新样本 $x$，计算 $p(x)$，若 $p(x)<ϵ$，则标记为异常：
$$p(x;\mu ,\Sigma )=\frac{1}{(2\pi {)}^{\frac{n}{2}}|\Sigma {|}^{\frac{1}{2}}}{e}^{-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )}$$

多元高斯分布与原始模型的联系：
原始模型其实就是多元高斯分布的一种特殊情况，即多元高斯分布的 $$\Sigma =\left[\begin{array}{cccc}{\sigma }_1^2& 0& \dots & 0\\ 0& {\sigma }_2^2& \dots & 0\\ ⋮& ⋮& \ddots & 0\\ 0& 0& \dots & {\sigma }_n^2\end{array}\right]$$

如何选择原始模型与多元高斯分布：

二、代码实现

import matplotlib.pyplot as plt
import seaborn as sns
sns.set(context="notebook", style="white", palette=sns.color_palette("RdBu"))

import numpy as np
import pandas as pd
import scipy.io as sio
from scipy import stats
from sklearn.model_selection import train_test_split

# 服务器响应的延迟(Latency)和吞吐量(Throughput)数据集
mat = sio.loadmat('../data/ex8data1.mat')

划分训练集、交叉验证集和测试集

X = mat.get('X')

# 将原始验证集数据分为验证集和测试集
Xval, Xtest, yval, ytest = train_test_split(mat.get('Xval'),
                                            mat.get('yval').ravel(),  # 将yval变成一维数组
                                            test_size=0.5)            # 测试集占50%

可视化训练数据

sns.regplot('Latency', 'Throughput', 
            data=pd.DataFrame(X, columns=['Latency', 'Throughput']),
            fit_reg=False,    # 不拟合直线
            scatter_kws={"s":20, "alpha":0.5})      # 对散点样式的设置

plt.show()

估计多元高斯分布的参数 $\mu$ 和 ${\sigma }^2$

# 设定初始 mu 值为均值
mu = X.mean(axis=0)     # [14.11222578 14.99771051]
# 计算协方差矩阵
cov = np.cov(X.T)       # [[ 1.83862041 -0.22786456],[-0.22786456  1.71533273]]

# 创建多元高斯模型
multi_normal = stats.multivariate_normal(mu, cov)

# 创建一个 2d 网格
x, y = np.mgrid[0:30:0.01, 0:30:0.01]
pos = np.dstack((x, y))

fig, ax = plt.subplots()

# 画概率密度，multi_normal.pdf(pos)：概率密度函数pdf对数据pos处理后得到满足设置的mean和cov的值
ax.contourf(x, y, multi_normal.pdf(pos), cmap='Blues')

# 画原始数据点
sns.regplot('Latency', 'Throughput',
            data=pd.DataFrame(X, columns=['Latency', 'Throughput']),
            fit_reg=False,
            ax=ax,
            scatter_kws={"s":10, "alpha":0.4})
plt.show()

选择阈值 $ϵ$
  1.使用训练集 $X$ 建立多元高斯模型
  2.使用交叉验证集(𝑋𝑣𝑎𝑙,𝑦𝑣𝑎𝑙)通过寻找最好的F1值来找到最好的 𝜖
$$F_1=2\cdot \frac{precision\cdot recall}{precision+recall}$$

def select_threshold(X, Xval, yval):
    """use CV data to find the best epsilon
    Returns:
        e: best epsilon with the highest f-score
        f-score: such best f-score
    """
    # 使用训练集数据来创建多元模型
    mu = X.mean(axis=0)
    cov = np.cov(X.T)
    multi_normal = stats.multivariate_normal(mu, cov)   # 设置模型的 mu 和 cov
    
    # 计算验证集数据的概率密度
    pval = multi_normal.pdf(Xval)
    
    # 设置候选的 epsilon
    epsilon = np.linspace(np.min(pval), np.max(pval), num=10000)   # 在min(pval) 和 max(pval) 中产生 10000 个候选值
    
    # 计算 f1-score
    fs = []
    for e in epsilon:
        y_pred = (pval <= e).astype('int')
        fs.append(f1_score(yval, y_pred))
        
    # 找到最好的 f1-score
    argmax_fs = np.argmax(fs)     # 返回最大值的索引
    
    return epsilon[argmax_fs], fs[argmax_fs]

from sklearn.metrics import f1_score, classification_report
e, fs = select_threshold(X, Xval, yval)
print(f'Best epsilon: {e}\nBest F-score on validation data: {fs}')

使用学到的 𝜖 可视化对验证集的预测
  1.使用验证数据集寻找最好的 𝜖
  2.使用训练集和验证集建立模型
  3.对测试集进行预测

# 1.使用验证数据集寻找最好的 epsilon
def select_threshold(X, Xval, yval):
    """use CV data to find the best epsilon
    Returns:
        e: best epsilon with the highest f-score
        f-score: such best f-score
    """
    # 使用训练集数据来创建多元模型
    mu = X.mean(axis=0)
    cov = np.cov(X.T)
    multi_normal = stats.multivariate_normal(mu, cov)
    
    # 关键：使用验证集数据微调超参数
    pval = multi_normal.pdf(Xval)
    
    # 设置候选的 epsilon
    epsilon = np.linspace(np.min(pval), np.max(pval), num=10000)
    
    # 计算 f1-score
    fs = []
    for e  in epsilon:
        y_pred = (pval <= e).astype('int')
        fs.append(f1_score(yval, y_pred))
        
    # 找到最好的 f1-score
    argmax_fs = np.argmax(fs)
    
    return epsilon[argmax_fs], fs[argmax_fs]

# 3.对测试集进行预测
def predict(X, Xval, e, Xtest, ytest):
    """with optimal epsilon, combine X, Xval and predict Xtest
    Returns:
        multi_normal: multivariate normal model
        y_pred: prediction of test data
    """
    Xdata= np.concatenate((X, Xval), axis=0)    # 数据拼接
    
    mu = Xdata.mean(axis=0)
    cov = np.cov(Xdata.T)
    multi_normal = stats.multivariate_normal(mu, cov)
    
    # 计算测试集的概率
    pval = multi_normal.pdf(Xtest)
    y_pred = (pval <= e).astype('int')
    
    print(classification_report(ytest, y_pred))
    
    return multi_normal, y_pred

multi_normal, y_pred = predict(X, Xval, e, Xtest, ytest)

# 搭建测试数据框架
data = pd.DataFrame(Xtest, columns=['Latency', 'Throughput'])
data['y_pred'] = y_pred

# 创建网格方便画图
x, y = np.mgrid[0:30:0.01, 0:30:0.01]
pos = np.dstack((x, y))

fig, ax = plt.subplots()

# 画概率密度
ax.contourf(x, y, multi_normal.pdf(pos), cmap='Blues')

# 画原始验证集的点
sns.regplot('Latency', 'Throughput',
            data=data,
            fit_reg=False,
            ax=ax,
            scatter_kws={"s":10,"alpha":0.4})

# 标记验证集上预测异常的数据，我们应该有一个测试集来做这件事
anamoly_data = data[data['y_pred']==1]
ax.scatter(anamoly_data['Latency'], anamoly_data['Throughput'], marker='x', s=50)
plt.show()