异常检测

一、密度估计：给定数据集$x^{\left(1\right)}$,$x^{\left(2\right)}$,…,$x^{\left(m\right)}$，对它进行数据建模p(x)（属于这组数据的可能性），当有新的数据$x_{test}$时：

高斯分布：变量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{{\left(x-\mu \right)}^2}{2{\sigma }^2}\right)$$
应用高斯分布开发异常检测算法：对于给定的数据集$x^{\left(1\right)}$,$x^{\left(2\right)}$,…,$x^{\left(m\right)}$，对每一个特征计算$\mu$和${\sigma }^2$的估计值。$${\mu }_j=\frac{1}{m}\sum _{i=1}^m{x}_j^{\left(i\right)}$$$${\sigma }_j^2=\frac{1}{m}\sum _{j=1}^m{\left(x_j^{\left(i\right)}-{\mu }_j\right)}^2$$
给定新的一个训练实例，根据模型计算p(x)：$$p\left(x\right)=\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)$$
当$p\left(x\right)<\epsilon$ 时，为异常。
开发一个异常检测系统时，从带标记（异常或正常）的数据着手，选择其中一部分正常数据构建训练集，用剩下的正常和异常混合而成的数据构成交叉检验集和测试集。具体评价方法如下：
1、根据测试集数据，估计特征的平均值和方差并构建p(x)函数；
2、对交叉检验集，使用不同的ε值作为阀值，并预测数据是否异常，根据F1值或者查准率与召回率的比例来选择ε；
3、选出ε后，针对测试集进行预测，计算异常检验系统的F1值，或者查准率与召回率之比来评估系统。
异常检测假设特征符合高斯分布，数据不是高斯分布时，最好将其转换成高斯分布，例如进行对数转换（log(x)、log(x+c)），指数转换（$\sqrt{x},x^{\frac{1}{3}}$）。
二、多元高斯分布
假设有两个相关特征，洋红色的线是一般的高斯分布模型获得的判定边界，绿色的X所代表的数据点很可能是异常值，但是其p(x)值仍然在正常范围内。而多元高斯分布将创建像图中蓝色曲线所示的判定边界。

一般的高斯分布是分别计算每个特征对应的概率然后将其累乘起来，而多元高斯分布模型将构建特征的协方差矩阵，用所有特征一起计算p(x)。首先计算所有特征的平均值，然后再计算协方差矩阵：$$\mu =\frac{1}{m}\sum _{i=1}^m{x}^{\left(i\right)}$$$$\sum =\frac{1}{m}\sum _{i=1}^m\left(x^{\left(i\right)}-\mu \right){\left(x^{\left(i\right)}-\mu \right)}^T=\frac{1}{m}{\left(X-\mu \right)}^T\left(X-\mu \right)$$$$p\left(x\right)=\frac{1}{{\left(2\pi \right)}^{\frac{n}{2}}{\left|\sum \right|}^{\frac{1}{2}}}exp\left(-\frac{1}{2}{\left(x-\mu \right)}^T{\Sigma }^{-1}\left(x-\mu \right)\right)$$
一般的高斯分布模型是多元高斯分布模型的一个子集，如果协方差矩阵只在对角线上为非零值时，即为一般的高斯分布模型了。多元高斯分布最大的优点是能描述两个特征变量之间可能存在正相关或者负相关的情况。

一般高斯分布模型和多元高斯分布模型的比较：

一般高斯分布模型	多元高斯分布模型
不能捕捉特征之间的相关性，但可以通过将特征进行组合来解决	自动捕捉特征之间的相关性
计算代价低，能适应大规模的特征	计算代价较高，训练集特征较少时适用
	必须要有m＞n，不然协方差矩阵不可逆，通常需要m＞10n，另外特征冗余也会导致协方差矩阵不可逆

三、以吴恩达机器学习课程练习材料实现，将异常检测算法应用于检测网络中的故障服务器。数据集特征包括每个服务器的吞吐量(mb/s)和响应的延迟(ms)。代码实现来源参考：吴恩达机器学习作业Python实现(八)：异常检测和推荐系统。
基于高斯分布的异常检测相关函数实现代码如下：

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

#计算高斯分布的参数
def getGaussianParams(X, useMultivariate):
    mu = X.mean(axis=0) #求每列的均值
    if useMultivariate:
        sigma2 = (X - mu).T @ (X - mu) / len(X) #多元高斯分布计算协方差矩阵
    else:
        sigma2 = X.var(axis=0, ddof=0)  #一般高斯分布计算方差，自由度为m
    return mu, sigma2

#概率密度函数
def gaussian(X, mu, sigma2):
    m, n = X.shape  #数据集行数、列数
    if np.ndim(sigma2) == 1:
        sigma2 = np.diag(sigma2)    #方差一维时，转化为对角阵
    norm = 1 / (np.power(2 * np.pi, n / 2) * np.sqrt(np.linalg.det(sigma2)))    #概率密度前半部分计算
    exp = np.zeros((m, 1))
    for row in range(m):    
        xrow = X[row]
        exp[row] = np.exp(-0.5 * (xrow - mu).T.dot(np.linalg.inv(sigma2)).dot(xrow - mu))   #概率密度后半部分计算
    return norm * exp

#计算F1值
def computeF1(yval, pval):  
        m = len(yval)   #获得交叉验证集样本数量
        tp = float(len([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    #召回率
        F1 = 2 * prec * rec / (prec + rec) if (prec + rec) else 0   #F1值
        return F1

#寻找最优（F1最大的）的阈值ε
def selectThreshold(yval, pval):
    epsilons = np.linspace(min(pval), max(pval), 1000)
    bestF1, bestEpsilon = 0, 0
    for e in epsilons:
        pval_ = pval < e    #小于阈值的设为1,否则为0
        thisF1 = computeF1(yval, pval_) #计算F1
        if thisF1 > bestF1:
            bestF1 = thisF1
            bestEpsilon = e #获得最大的F1值以及对应最优的阈值
    return bestF1, bestEpsilon

可视化数据点及等概率曲线的代码如下：

#可视化数据点
def plot_data():
    plt.figure(figsize=(8,5))   #初始化图像并设置图的大小
    plt.plot(X[:,0],X[:,1],'bx')    #绘制数据点
    #plt.show()

#绘制等概率曲线
def plotContours(mu,sigma2):
    delta=0.3
    x=np.arange(0,30,delta)
    y=np.arange(0,30,delta)
    xx,yy=np.meshgrid(x,y)  #生成网格点坐标矩阵
    points=np.c_[xx.ravel(),yy.ravel()] #按列合并，一列横坐标，一列纵坐标
    z=gaussian(points,mu,sigma2)    #计算概率密度
    z=z.reshape(xx.shape)
    cont_levels=[10**h for h in range(-20,0,3)] #概率等高线的取值
    plt.contour(xx,yy,z,cont_levels)    #绘制等高线
    plt.title('Gaussian Contours',fontsize=16)  #设置标题
    #plt.show()

数据获取及进行异常检测代码如下：

mat=loadmat('anomaly_de.mat')  #加载数据集
#print(mat.keys())
X=mat['X']  #获得测试集
Xval,yval=mat['Xval'],mat['yval']   #获得交叉验证集
#print(X.shape,Xval.shape,yval.shape)
mu,sigma2=getGaussianParams(X,True) #计算高斯分布的参数
pval=gaussian(Xval,mu,sigma2)   #获得交叉验证集的概率密度
bestF1,bestEpsilon=selectThreshold(yval,pval)   #获得最大的F1值以及对应最优的阈值
print(bestF1,bestEpsilon)
y=gaussian(X,mu,sigma2) #获得测试集的概率密度
xx=np.array([X[i] for i in range(len(y)) if y[i]<bestEpsilon])  #概率小于阈值记为异常点
#print(xx)
plot_data() #绘制数据点
plotContours(mu,sigma2) #绘制等概率曲线
plt.scatter(xx[:,0],xx[:,1],s=80,facecolors='none',edgecolors='r')  #标记异常点
plt.show()

利用一般高斯分布和多元高斯分布检测得到的异常点分别如下：


可以看到两种方法检测得到的异常点是一样的，在该数据集上效果差不多。尽管F1值一致，但是判断异常点的阈值ε并不完全一样，一般高斯分布ε为8.99985263e-05，多元高斯分布ε为9.07484457e-05。

(结语个人日记：很久没有联系之后，前段时间又和高中同学聊了一下天。感慨人跟人之间的关系是很神奇的，在一定的阶段会有一群朋友，当那阶段过去，可能彼此就只变成了通讯录里躺着的名字，渐行渐远，最终留下的好友只有那么一两个。而有的朋友虽然不常联系，但是每次时隔很久聊起天，也总能说个不完。很多东西都被时间推着走，还是庆幸身边存在着一两个可以说个不完的朋友，以及愿意继续陪伴着走下去的人。）