吴恩达老师机器学习-ex8

news2024/9/25 1:20:53

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()

开发异常检测系统的步骤：

选择合适的特征变量 $x_{j}$
拟合参数 $(\mu _{1},\mu _{2},\mu _{3},...,\mu _{n},\sigma _{1}^{2},\sigma _{2}^{2},\sigma _{3}^{2},...,\sigma _{n}^{2})$ 。其中 $\mu _{j}=\frac{1}{m}\sum_{i=1}^{m}x_{j}^{(i)}$ ， $\sigma ^{2}=\frac{1}{m}\sum_{i=1}^{m}(x_{j}^{(i)}-\mu _{j})^{2}$
对于给定的 $x_{test}$ 计算 $P(x_{test})=\prod_{j=1}^{n}P(x_{test};\mu ,\sigma ^{2})=\prod_{j=1}^{n}\frac{1}{\sqrt{2\pi }\sigma _{j}}exp(-\frac{(x_{i}-\mu _{i})^{2}}{2\sigma ^{2}})$
如果 $P(x_{test})\geq \varepsilon$ ，则认为该样本正常；如果 $P(x_{test})< \varepsilon$ ，则认为该样本异常。

我们可以通过可视化数据发现，数据之间并没有线性相关。所以在这，我采用的是一般的高斯分布模型。

首先，计算数据集的平均值以及方差。

#求均值和方差
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

接着，很重要的一步就是阈值 $\varepsilon$ 的确定。可以尝试不同的 $\varepsilon$ 作为阈值，使得最后的 $F1-score$ 最大。（通过计算真阳性、真阴性、假阳性和假阴性计算查准率和召回率）

#确定阈值
def get_threshold(yval,p):
    best_F1= 0
    best_threshold = 0
    step = (np.max(p) - np.min(p)) / 1000
    for th in np.arange(np.min(p), np.max(p), step):
        y_p = p < th
        tp = 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 = 0
        if (tp+fn):
            rec = tp / (tp+fn)
        else:
            rec = 0
        if (prec+rec):
            F1 = (2*prec*rec)/(prec + rec)
        else:
            F1 = 0
        if F1 > best_F1:
            best_F1 = F1
            best_threshold = th
    return 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))

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/1981635.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！