吴恩达机器学习第八次编程作业-异常检测

机器学习

Programming Exercise 8: Anomaly Detection

目标： 检测服务器的异常行为

工具： Pycharm，Python3.6

参考资料

  1. 吴恩达机器学习作业Python实现(八)：异常检测和推荐系统
  2. Andrew-NG-Meachine-Learning
  3. 斯坦福大学2014机器学习教程中文笔记目录

完整代码：

  1. Anomaly Detection

绘制数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.io as sio

def Get_data(path):
    data=sio.loadmat(path)
    # for key in data:
    #     print(key)
    return data

def Plot_Data(X):
    plt.plot(X.T[0],X.T[1],'bx')
    plt.xlabel('Latency (ms)',fontsize=10)
    plt.ylabel('Throughput (mb/s)',fontsize=10)
    plt.xlim(0,30)
    plt.ylim(0,30)
    plt.show()

data=Get_data('./ex8/ex8data1.mat')
X,Xval,yval=data['X'],data['Xval'],data['yval']
#print(X.shape,Xval.shape,yval.shape) #  (307, 2) (307, 2) (307, 1)
Plot_Data(X)

Figure_1.png

Gaussian distribution

_TP`3~L_C~WEWMK0JUP_F@2.png

Estimating parameters for a Gaussian

_PVWV_T`~_TMWO_A_AG81IS.png

这个地方花了很多时间，没有搞明白∑和σ之间的关系，文档所给的矩阵运算貌似有问题，参考别人的作业才能实现。
在课程中原始高斯概率模型和多元高斯概率模型区别听的不是很懂，而所给的文档中使用的Octave实现的代码用的是多元高斯模型

求μ和σ
1
2
3
4
5
6
7
8
9
def Estimate_Gaussian(X):
    mu=X.mean(axis=0)
    sigma2=X.var(axis=0,ddof=0)
    '''
    # 当使用多元高斯概率模型时
    sigma2 = ((X-mu).T @ (X-mu)) /X.shape[0]
    # 这个地方也不明白
    '''
    return mu,sigma2
多元高斯概率分布
1
2
3
4
5
6
7
def Multivariate_Gaussian(X,mu,sigma2):
    # print(X.shape,mu.shape,sigma2.shape)
    if np.ndim(sigma2)==1:
        sigma2=np.diag(sigma2) # Extract a diagonal or construct a diagonal array
    norm=1/(np.power(2*np.pi,X.shape[1]/2)*np.power(np.linalg.det(sigma2),0.5))
    p=np.exp(-1/2*(X-mu).dot(np.linalg.inv(sigma2)).dot((X-mu).T)) # 这里和视频中所给的公式不太一样，可能默认的矩阵形状不同，所给资料矩阵是(307，2）
    return norm*np.diag(p) # 为什么取对角线
绘制高斯概率分布图
1
2
3
4
5
6
7
8
9
10
11
def Plot_Contour(X,mu,sigma2):
    plt=Plot_Data(X)
    x=np.arange(0,30,0.5)
    y=np.arange(0,30,0.5)
    xx,yy=np.meshgrid(x,y)
    points=np.c_[xx.ravel(),yy.ravel()] # 按列拼接
    z=Multivariate_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.show()

Figure_2.png

Selecting the threshold ε

选择不同的 ε，根据F1的值来判定阈值是否合适
9`UJ7F`K1~0SCP_MNE_X_LU.png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
def Select_Threshold(yval,pval):
    best_epsilon=0
    best_F1=0
    F1=0
    epsilons=np.linspace(min(pval),max(pval),1000)
    for e in epsilons:
        pval_tmp=pval<e # True or False vector
        Tp=np.sum(pval_tmp&yval) # correctly classified
        Fp=np.sum(pval_tmp&(yval^1)) # incorrectly classified
        Fn=np.sum((pval_tmp^1)&yval) # incorrectly classified
        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
        if F1>best_F1:
            best_F1=F1
            best_epsilon=e
    return best_epsilon,best_F1
pval=Multivariate_Gaussian(Xval,mu,sigma2)
epsilon,F1 = Select_Threshold(yval.ravel(), pval) # 注意维度
print(epsilon,F1)  # 8.999852631901397e-05 0.8750000000000001

High dimensional dataset

这里主要就是使用之前的代码跑高维度特征
核心是Multivariate_Gaussian函数，根据估计的参数 μ 和 σ ，传入数据 X ，计算每个样本的高斯概率
使用验证集计算高斯概率 pval，然后和yval对比选择最佳 ε
最后计算训练样本的异常数量

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def High_Dimensional_Dataset(path):
    data=Get_data(path)
    # print(data.keys()) # dict_keys(['__header__', '__version__', '__globals__', 'X', 'Xval', 'yval'])
    X,Xval,yval=data['X'],data['Xval'],data['yval']
    # print(X.shape,Xval.shape,yval.shape) # (1000, 11) (100, 11) (100, 1)
    mu,sigma2=Estimate_Gaussian(X)
    pval=Multivariate_Gaussian(Xval,mu,sigma2)
    epsilon,F1=Select_Threshold(yval.ravel(),pval)
    print('Best epsilon found using cross-validation: %e' %epsilon)
    print('Best F1 on Cross Validation Set:  %f' %F1)
    pval=Multivariate_Gaussian(X,mu,sigma2)
    print('Outliers found:  %d' %(sum(pval<epsilon)))
    '''
    Best epsilon found using cross-validation: 1.378607e-18
    Best F1 on Cross Validation Set:  0.615385
    Outliers found:  117
    '''