data1
导入库,读取数据并进行可视化
因为这次的数据是mat文件,需要使用scipy库中的loadmat进行读取数据。
通过对数据类型的分析,发现是字典类型,查看该字典的键,可以发现又X等关键字。
import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt#读取数据
path = "./ex8data1.mat"
data = sio.loadmat(path)
# print(data.keys())
X = data.get("X")
Xval = data.get("Xval")
yval = data.get("yval")
# print(X.shape)#可视化数据
plt.scatter(X[:,0],X[:,1])
plt.show()
开发异常检测系统的步骤:
- 选择合适的特征变量
- 拟合参数。其中,
- 对于给定的计算
- 如果,则认为该样本正常;如果,则认为该样本异常。
我们可以通过可视化数据发现,数据之间并没有线性相关。所以在这,我采用的是一般的高斯分布模型。
首先,计算数据集的平均值以及方差。
#求均值和方差
def fit_parameters(X):m = len(X)means = np.mean(X,axis=0)cov = np.var(X,axis=0)return means,cov
然后,计算密度估计。
#进行密度估计
def probability(X,means,cov):m = len(X)p = np.exp(-np.power(X-means,2)/(2 * cov))/((np.power(2 * np.pi, 0.5) * np.power(cov, 0.5)))p_new = np.ones((m,1))for i in range(p.shape[1]):for j in range(p.shape[0]):p_new[j] = p_new[j] * p[j,i]# p_new = p_new * p[:,i]return p_new
接着,很重要的一步就是阈值的确定。可以尝试不同的作为阈值,使得最后的最大。(通过计算真阳性、真阴性、假阳性和假阴性计算查准率和召回率)
#确定阈值
def get_threshold(yval,p):best_F1= 0best_threshold = 0step = (np.max(p) - np.min(p)) / 1000for th in np.arange(np.min(p), np.max(p), step):y_p = p < thtp = np.sum((yval == 1) & (y_p == 1))fp = np.sum((yval == 0) & (y_p == 1))tn = np.sum((yval == 0) & (y_p == 0))fn = np.sum((yval == 1) & (y_p == 0))if (tp+fp):prec = tp / (tp+fp)else:prec = 0if (tp+fn):rec = tp / (tp+fn)else:rec = 0if (prec+rec):F1 = (2*prec*rec)/(prec + rec)else:F1 = 0if F1 > best_F1:best_F1 = F1best_threshold = threturn best_F1,best_threshold
调用函数,并找到异常点。
means, cov = fit_parameters(X)
# print(cov)
p_new = probability(X,means,cov)
# print(p_new)
best_F1,best_threshold = get_threshold(yval,p_new)
# print(best_F1)
# print(best_threshold)plt.scatter(X[:,0],X[:,1],c="b")
for i in range(len(X)):if p_new[i] < best_threshold:plt.scatter(X[i,0],X[i,1],c="r")
plt.show()
data2
data2对应的是多维的数据
print("-------------data2-------------")#读取数据
path2 = "./ex8data2.mat"
data2 = sio.loadmat(path2)
# print(data2.keys())
X2 = data2.get("X")
Xval2 = data2.get("Xval")
yval2 = data2.get("yval")means2, cov2 = fit_parameters(X2)
# print(cov)
p2_new = probability(Xval2,means2,cov2)
# print(p_new)
best2_F1,best2_threshold = get_threshold(yval2,p2_new)
p2_new = probability(X2,means2,cov2)
print(best_F1)
print(best_threshold)anoms = []
for i in range(len(X2)):if p2_new[i] < best2_threshold:anoms.append(X2[i,:])
print(len(anoms))