异常监测的要点:1. 适用于数据集符合某种分布,能够转换为某种分布也算,比如车的航行轨迹,就不能用这招。 2. 或者使用阈值设定,结合逻辑回归设定异常,也可以。3. 在数据集中,异常数据点非常少,1%都算多。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.io as spio
def display_2d_data(X,marker):
#plt.figure(figsize=(10,8))
plt.plot(X[:,0],X[:,1],marker)
return plt
def estimateGaussian(X): #求均值与方差
m,n = X.shape
mu = np.zeros((n,1))
sigma2 = np.zeros((n,1))
mu = np.mean(X,axis = 0)
sigma2 = np.var(X,axis = 0) * (m-1) / m
# 在概率论中计算sigma2时,除以的是(1-m),机器学习中,除以(1-m)和除以m的差别不大
return mu,sigma2
def multivariateGaussian(X,mu,Sigma2):
k = len(mu)
if (Sigma2.shape[0]>1):
Sigma2 = np.diag(Sigma2)
'''多元高斯分布函数'''
X = X-mu
argu = (2*np.pi)**(-k/2)*np.linalg.det(Sigma2)**(-0.5)
p = argu*np.exp(-0.5*np.sum(np.dot(X,np.linalg.inv(Sigma2))*X,axis=1))
# axis表示每行
return p
def visualizeFit(X,mu,sigma2):
x = np.arange(0, 36, 0.5) # 0-36,步长0.5
y = np.arange(0, 36, 0.5)
X1,X2 = np.meshgrid(x,y) # 要画等高线,所以meshgird
Z = multivariateGaussian(np.hstack((X1.reshape(-1,1),X2.reshape(-1,1))), mu, sigma2) # 计算对应的高斯分布函数
Z = Z.reshape(X1.shape) # 调整形状
plt.figure(figsize=(10,8))
plt.plot(X[:,0],X[:,1],'bx')
if np.sum(np.isinf(Z).astype(float)) == 0: # 如果计算的为无穷,就不用画了
#plt.contourf(X1,X2,Z,10.**np.arange(-20, 0, 3),linewidth=.5)
CS = plt.contour(X1,X2,Z,10.**np.arange(-20, 0, 3),color='black',linewidth=.5)
# 画等高线,Z的值在10.**np.arange(-20, 0, 3)
#plt.clabel(CS)
plt.show()
# 选择最优的epsilon,即:使F1Score最大
def selectThreshold(yval,pval):
'''初始化所需变量'''
bestEpsilon = 0.
bestF1 = 0.
F1 = 0.
step = (np.max(pval)-np.min(pval))/1000
'''计算'''
for epsilon in np.arange(np.min(pval),np.max(pval),step):
cvPrecision = pval<epsilon
tp = np.sum((cvPrecision == 1) & (yval == 1).ravel()).astype(float) # sum求和是int型的,需要转为float
fp = np.sum((cvPrecision == 1) & (yval == 0).ravel()).astype(float)
fn = np.sum((cvPrecision == 0) & (yval == 1).ravel()).astype(float)
precision = tp/(tp+fp) # 精准度
recision = tp/(tp+fn) # 召回率
F1 = (2*precision*recision)/(precision+recision) # F1Score计算公式
if F1 > bestF1: # 修改最优的F1 Score
bestF1 = F1
bestEpsilon = epsilon
return bestEpsilon,bestF1
def AnomalyDetection2():
data = spio.loadmat("data2.mat")
X = data['X']
Xval = data['Xval']
yval = data['yval']
#print(pd.DataFrame(X))
mu,sigma2 = estimateGaussian(X)
#print(mu,sigma2)
p = multivariateGaussian(X,mu,sigma2)
#print(pd.DataFrame(p))
pval = multivariateGaussian(Xval,mu,sigma2)
epsilon,F1 = selectThreshold(yval,pval)
print("the best epsilon is ",epsilon)
print("the best F1 is ",F1)
print("Outliers found",np.sum(p < epsilon))
AnomalyDetection2()
def AnomalyDetection():
data = spio.loadmat("data.mat")
#print(data)
X = data['X']
Xval = data['Xval']
yval = data['yval'] # y = 1 为异常
#plt.plot(X[:,0],X[:,1],'x')
plt = display_2d_data(X,'x')
plt.title("Origin data")
plt.show()
mu,sigma2 = estimateGaussian(X)
#print(mu,sigma2)
p = multivariateGaussian(X,mu,sigma2)
#print(p)
visualizeFit(X,mu,sigma2) # 显示图像
# 选择异常点
pval = multivariateGaussian(Xval,mu,sigma2)
epsilon,F1 = selectThreshold(yval,pval)
print(u'在CV上得到的最好的epsilon是:%e'%epsilon)
print(u'对应的F1Score值为:%f'%F1)
outliers = np.where(p<epsilon) # 找到小于临界值的异常点,并作图
#plt.figure(figsize=(10,8))
plt.plot(X[outliers,0],X[outliers,1],'o',markeredgecolor='r',markerfacecolor='w',markersize=10.)
plt = display_2d_data(X, 'bx')
plt.show()
if __name__ == "__main__":
AnomalyDetection()