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

题目概述： 在本练习中，我们将使用高斯模型实现异常检测算法，并将其应用于检测网络上的故障服务器。 我们还将看到如何使用协作过滤构建推荐系统，并将其应用于电影推荐数据集。

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

data = loadmat('E:/shujuji/ex8data1.mat')
X = data['X']
X.shape #(307, 2)
#可以看到是二维的

fig,ax=plt.subplots(figsize=(12,8))
ax.scatter(X[:,0], X[:,1])
plt.show()

计算高斯分布参数：

要做的是，输入一个X矩阵，输出2个n维的向量，mu包含了每一个维度的平均值，sigma包含了每一个维度的方差。

def estimate_gaussian(X):
    mu=X.mean(axis=0) #axis=0,求每列平均值
    sigma=X.var(axis=0)
    return mu, sigma

mu, sigma = estimate_gaussian(X)
mu, sigma

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

选择阈值ε
有了参数后，可以估计每组数据的概率，低概率的数据点更可能是异常的。确定异常点需要先确定一个阈值，我们可以通过验证集来确定这个阈值。

Xval = data['Xval']
yval = data['yval']
Xval.shape, yval.shape #((307, 2), (307, 1))

我们还需要一种计算数据点属于正态分布的概率的方法,幸运的是SciPy有这个内置的方法

from scipy import stats
'''
练习实现正态分布的概率函数
dist=stats.norm(mu[0],sigma[0])
#stats.norm实现正态分布,norm(loc,scale),loc:mean均值,scale:标准差
dist.pdf(X[:,0])[0:50] #pdf概率密度函数
'''
p=np.zeros((X.shape[0],X.shape[1])) #注意两个括号，2维
p[:,0]=stats.norm(mu[0], sigma[0]).pdf(X[:,0])#取出第一列来概率密度。
p[:,1]=stats.norm(mu[1], sigma[1]).pdf(X[:,1])#取出第二列来概率密度。

pval = np.zeros((Xval.shape[0], Xval.shape[1]))
pval[:,0] = stats.norm(mu[0], sigma[0]).pdf(Xval[:,0])
pval[:,1] = stats.norm(mu[1], sigma[1]).pdf(Xval[:,1])
pval.shape #(307, 2)

我们需要一个函数，找到给定概率密度值和真实标签的最佳阈值。 为了做到这一点，我们将为不同的epsilon值计算F1分数。 F1是真阳性，假阳性和假阴性的数量的函数

#在验证集上选择
def select_threshold(pval,yval):
    best_epsilon=0
    best_f1=0
    f1=0
    step=(pval.max()-pval.min())/1000
    
    #epsilon在pval最小最大值遍历,找到best_epsilon
    for epsilon in np.arange(pval.min(), pval.max(), step):
        preds=pval<epsilon  #如果pval<espilon成立，返回布尔值1，反之0
        
        tp=np.sum(np.logical_and(preds==1,yval==1)).astype(float) 
        fp=np.sum(np.logical_and(preds== 1, yval == 0)).astype(float)
        fn=np.sum(np.logical_and(preds== 0, yval == 1)).astype(float)
        #np.logical_and:逻辑与
        #astype：转换数组的数据类型
        
        precision=tp/(tp+fp)
        recall=tp /(tp + fn)
        f1=(2 * precision * recall) / (precision + recall)
        
        if f1>best_f1:
            best_f1 = f1
            best_epsilon = epsilon
            
    return best_epsilon, best_f1           

epsilon, f1 = select_threshold(pval, yval)
epsilon, f1

(0.009566706005956842, 0.7142857142857143)

最后，我们可以将阈值应用于数据集，并可视化结果

outliers = np.where(p < epsilon)#np.where返回索引
outliers

(array([300, 301, 301, 303, 303, 304, 306, 306], dtype=int64),
array([1, 0, 1, 0, 1, 0, 0, 1], dtype=int64))

fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(X[:,0], X[:,1])
ax.scatter(X[outliers[0],0], X[outliers[0],1], s=50, color='r', marker='o')
plt.show()

红色为异常点