开发和评价一个异常检测系统

解决问题：此次练习是为了检测服务器的吞吐量(throughput)和响应延迟(latency)是否有异常。
问题背景：收集307个训练样本，猜测全都是正常的(但是实际中可能有几个异常点)，所以需要用高斯分布检测异常样本。
可以先用2D散点图查看分布情况(part1图)，用测试机拟合高斯分布然后配合验证集的得到的epision找到异常点，
最后应用到多维度的大数据中。

开发和评价一个异常检测系统
1 part2 ：估计高斯分布参数
用训练集获得均值mu和方差sigma用来构建p(x)函数
2. part3 : 选择阀值 epision
使用带标签的交叉验证集，根据F1值来确定epision
3. part3: 找出异常值outliers
选出 epision 后，针对测试集进行异常值预测(outliers=P(x) < epision )，同时可计算下测试集的F1值，或者召回率与精确率

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns;sns.set()
from sklearn.model_selection import train_test_split
from  scipy.io import loadmat
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline

part1 数据加载及展示

ex8data1 = loadmat('ex8data1.mat')
X = ex8data1['X']
Xval = ex8data1['Xval']
yval = ex8data1['yval']
sns.jointplot(X[:,0], X[:,1], kind="hex", color="purple")

part2 通过高斯分布估计每个特征的均值及方差用于计算概率值

# ================== Part 2: Estimate the dataset statistics ===================
def estimateGaussian(X):
    mu = X.mean(axis=0)
    sigma2 = X.var(axis=0)
    return mu, sigma2
mu, sigma2 = estimateGaussian(X)
mu, sigma2

(array([14.11222578, 14.99771051]), array([1.83263141, 1.70974533]))

def plot_x1_x2():
    # 𝑧轴为根据两个特征的值所估计𝑝(𝑥)值
    delta = 0.25 
    x1 = np.arange(5, 25, delta)
    x2 = np.arange(5, 25, delta)
    X1, X2 = np.meshgrid(x1, x2)
    # x1 代表第一个特征，x2代表第二个特征
    p1 = 1 / (np.sqrt(2*np.pi* sigma2[0] ) ) *np.exp(-(X1- mu[0])**2/(2*sigma2[0]))
    p2 = 1 / (np.sqrt(2*np.pi* sigma2[1] ) ) *np.exp(-(X2 - mu[1])**2/(2*sigma2[1]))
    p = p1 * p2
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_surface(X1, X2, p, cmap='autumn') 
plot_x1_x2()

注：plot_x1_x2中的p就是multivariateGaussian中的p

def multivariateGaussian(X, mu, Sigma2):
    # 多维高斯分布
    Sigma2 = np.diag(Sigma2)   # (2,2)
    n = len(mu)
    p =(2 * np.pi)**(-n/2) * np.linalg.det(Sigma2) ** (-0.5) *\
            np.exp(-0.5 * np.sum((X - mu) @ np.linalg.pinv(Sigma2) * (X-mu), axis=1))
    return p
# 根据X自身的均值mu和方差sigma2来计算多维高斯概率值
p = multivariateGaussian(X, mu, sigma2)

$注：multivariateGaussian中的Sigma就是协方差\Sigma 的主对角线，不过没有考虑特征之间的相关性，为0$

def visualizeFit(X, mu, sigma2, anomaly=False):
    # 通过测试集的mu和sigma2计算2个特征的p值存入Z，画出等高线
    x = np.arange(0,35.1,0.5)
    y = np.arange(0,35.1,0.5)
    X1, X2 = np.meshgrid(x,y)
    Z = multivariateGaussian(np.vstack((X1.ravel(), X2.ravel())).T, mu, sigma2)
    Z = Z.reshape(X1.shape)    
    plt.scatter(X[:,0], X[:,1], s=20)
    plt.contour(X1, X2,Z,levels=[1E-20,1E-17,1E-14,1E-11,1E-8,1E-5,1E-2])
    if anomaly:
        plt.scatter(X[outliers][:,0], X[outliers][:,1], color='', marker='o', edgecolors='g')
visualizeFit(X,  mu, sigma2)

part3 通过带标签的验证集计算F1值确定epision

pval = multivariateGaussian(Xval, mu, sigma2)
def selectThreshold(yval, pval):
    # 根据F1得分获取最合适的epsilon值作为估计是否为异常值（F1值越大越好，epsilon为验证集的高斯概率值）
    # 注意要讲异常值标记为阳性，所以cvPredictions = pval < epsilon
    bestF1 = 0
    stepsize = (max(pval) - min(pval))/1000
    for epsilon in np.arange(min(pval), max(pval)+stepsize, stepsize):
        cvPredictions = pval < epsilon
        tp = sum((cvPredictions == 1) & ((yval == 1).ravel()))
        fp = sum((cvPredictions == 1) & ((yval == 0).ravel()))
        fn = sum((cvPredictions == 0) & ((yval == 1).ravel()))
        rec = tp/(tp+fp)
        prec = tp/(tp+fn)
        F1 = 2 * rec * prec / (rec+prec)
        if F1 > bestF1:
            bestF1 = F1
            bestEpsilon = epsilon
    return bestEpsilon, bestF1
[epsilon, F1] = selectThreshold(yval, pval)
# 找出小于epsilon的异常的点
outliers = np.where( (p<epsilon) == 1)
# 画出异常点
visualizeFit(X, mu, sigma2, anomaly=True)

part4 确定上面part没问题后测试更大更多的数据集

X = loadmat('ex8data2.mat')['X']
Xval = loadmat('ex8data2.mat')['Xval']
yval = loadmat('ex8data2.mat')['yval']
X.shape, Xval.shape, yval.shape

((1000, 11), (100, 11), (100, 1))

mu,sigma2 = estimateGaussian(X)
p = multivariateGaussian(X, mu, sigma2)
pval = multivariateGaussian(Xval, mu, sigma2)
[epsilon,F1]= selectThreshold(yval, pval)
epsilon,F1

(1.377228890761358e-18, 0.6153846153846154)

PS

$其实上面的multivariateGaussian是原高斯分布模型，只要改变Sigma2=(X-mu).T@(X-mu)/len(X)就是多元高斯分布模型$