异常检测【Coursera 斯坦福 机器学习】

本文基于Coursera 斯坦福吴恩达机器学习课程

谢绝任何不标记出处的转载

如有问题请联系作者

1. 异常检测 Anomaly Dectection概念

顾名思义，anomaly detection 是用来检测出现异常的情况的。如下图所示，异常（红叉）往往会脱离正常值（蓝圈）范围，变得比较明显。这很好判断，难点就是在如何界定正常值的范围了（类似于确定decision boundary）。这里我们会使用一个threshold ε，当通过模型产生的估计值p(x)<ε时，我们认为该例子为异常的（unusual/anomalous）。当我们侦测到的异常值远超过原本应有的值的数量时，应该适当降低ε进行调节。

异常检测的应用也很广泛，例如：

（1）制造业产品的异常，例如飞机引擎。

（2）诈骗侦测：往往异常用户的一些feature跟普通用户会有差别。

（3）监测指标：例如数据中心的电脑是否正常工作。

2. 异常检测算法 （density estimation）

我们假设每个feature都服从一个均值为μj，标准差为σj的正态分布，那么有：

因为这个思想很简单，所以我用R试了一下：

# 首先set seed

set.seed(7)

# 产生数据，有两个feature x1, x2
# 我们设置数据的后100个全为同类型的异常值（不属于正态分布且全为0。在实际生活中，异常值会是多种多样的，甚至异常的方式都未知），这里注意，我们的x1和x2的均值相差非常大。

x1.1 <- rnorm(900,500,10)
x1.2 <- numeric(100)
x1 <- c(x1.1, x1.2)
x2.1 <- rnorm(900,80,200)
x2.2 <- numeric(100)
x2 <- c(x2.1, x2.2)

data1 <- data.frame(x1,x2)

# 画图观察数据，发现非常imbalanced

plot(data1$x1,data1$x2)

# 我们暂时不对数据做处理，继续进行学习： miu1.1 <- mean(x1) miu1.2 <- mean(x2) sd1.1 <- sd(x1) sd1.2 <- sd(x2) # 设置threshold为0.02 epsilon <- 0.02 #计算p(x) p1.1 <- 1/sqrt(2*pi*sd1.1)*exp(-(x1-miu1.1)^2/2/sd1.1^2) p1.2 <- 1/sqrt(2*pi*sd1.2)*exp(-(x2-miu1.2)^2/2/sd1.2^2) p1 <- p1.1*p1.2 #将结果进行01分类 for(i in 1:length(p1)) { if(p1[i]<epsilon){ p1[i]=1 } else p1[i]=0 } #展示结果，理论上后100全为1，其余全是0 p1

> p1
   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [334] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [371] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [408] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [445] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [482] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [519] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [556] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [593] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [630] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [667] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [704] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [741] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [778] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [815] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [852] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [889] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [926] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [963] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1000] 1

# 结果显示，1000个值全部被判断为异常值，这时我们进行feature scaling再重新学习：

x1 <- (x1-mean(x1))/sd(x1)
x2 <- (x1-mean(x2))/sd(x2)

data1 <- data.frame(x1,x2)

plot(data1$x1,data1$x2)

miu1.1 <- mean(x1) miu1.2 <- mean(x2) sd1.1 <- sd(x1) sd1.2 <- sd(x2) # 设置threshold为0.02 epsilon <- 0.02 #计算p(x) p1.1 <- 1/sqrt(2*pi*sd1.1)*exp(-(x1-miu1.1)^2/2/sd1.1^2) p1.2 <- 1/sqrt(2*pi*sd1.2)*exp(-(x2-miu1.2)^2/2/sd1.2^2) p1 <- p1.1*p1.2 #将结果进行01分类 for(i in 1:length(p1)) { if(p1[i]<epsilon){ p1[i]=1 } else p1[i]=0 } #展示结果，全部分类正确 p1

> p1
   [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [112] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [149] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [186] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [223] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [260] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [297] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [334] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [371] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [408] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [445] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [482] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [519] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [556] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [593] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [630] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [667] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [704] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [741] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [778] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [815] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [852] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [889] 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [926] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [963] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1000] 1

# 尝试增大epsilon

epsilon为0.1, 1时 结果与0.02同，在为2时出现overfit，在为3时全部值被判断为异常值（过拟合）；

# 尝试减小epsilon

epsilon在小至0.0003时结果依然与0.02同，在小于0.0002，所有值均为正常值（欠拟合）。

通过这个例子，大致可以得出epsilon的取值范围在[0.0003,1]之间会比较好。并且，在做异常检测时，一定要对feature进行scaling. 

3. 评估异常检测

依然可以用cross-validation来做，然后做confusion matrix即可。此外，可以用true positive, false positive, false negative ,true negative; precision/recall;F1-score来评估。

但是要注意正态分布是否是skewed的，如果分类skewed，不适合用classfication accuracy来做。

通常异常检测y=1和y=0的情况是非常imbalanced的（不然异常值太多了。。想想就可怕）。

4. Anomaly Detection vs. Supervised Learning

5. feature的选择

（1）要尽量选择属于正态分布的feature，且越符合标准正态分布的形状越好（没有skewed）。如果feature skewed, 可以用log或者sqrt来修正。

例如我们假设feature x3（=x1+x2，x1和x2均为上面例子中的变量，懒得重新做feature了）, 那么直接出来的直方图为：

log一下（加了个3在里面）：

sqrt一下：

（2）除此之外，还要选择当出现异常时，outlier会很突出的feature（比如诈骗中，这个用户一天打出去了500个电话，可能正常人很难做到，这个outlier就很突出；或者监测某一仪器工作时，仪器的指标突然变为0，其他仪器的指标多多少少都在100-500的范围内，这个outlier也很突出）

（3）除了利用现有features，通过比例等形式建造新的feature再使用（1）（2）也是不错的选择。

6. 多元正态分布下的异常检测

上面提到的p(x)=p(x1)*p(x2)*...p(xn)是有局限性的（但是更常用，更简单也更便宜，computationally cheaper）。举个栗子（如下图）。从左图我们很容易就能辨别出来两个红叉为异常值，可是投影到右图，可以看到两个异常值在x1,x2分别的正态分布下还算可以，不算特别奇怪。那么使用p(xj)乘积的形式就不太好探查了（虽然我们有乘积）。这时，可以用多元正态分布一次性将p(x)计算出来。

下图为课程中多元正态分布的示例：

此时的单个均值μj组成的n*1向量μ，ε就变成了一个n*n的矩阵Ε（covariance matrix）。假设各个变量独立（除对角线外都为0），观测值数量远大于feature数量（10倍以上吧），且E不可逆（non-invertible），我们有公式：

|E|被称为“determinant of Ε”。在matlab中用det(Sigma)求得。同样的，p(x)＜ε即为anomaly. 

有一个疑问，这里可以跟AR(1)等correlation structure联系起来么？如果要联系，就要求模型假设为残差独立/features独立，AR(1)是用于自相关的，features独立就是共线性等（也可能不是线性）的问题。

——从课程上来看是可以的，而且可能用了自回归的structure会对模型更好。但是还是不太明白这里的模型假设。

这里记录一下课堂问题，第一次做错的有点离谱，有点绕：