本次编程作业的实现环境是Python3、Anaconda3(64-bit)、Jupyter Notebook。是在深度之眼“机器学习训练营”作业基础上完成的，个别代码有修改，供交流学习之用。
Anomaly detection（异常检测）?
我们的第一个任务是使用高斯模型来检测数据集中未标记的示例是否应被视为异常。 我们有一个简单的二维数据集开始，以帮助可视化该算法正在做什么。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sb
from scipy.io import loadmat
#读取ex8data1的数据
data = loadmat('data/ex8data1.mat')
#data
for key in data.keys():
print(key, ' ')
X = data['X']
X.shape, data['Xval'].shape, data['yval'].shape
print(X)
#将data['X']的二维数据绘制出来
fig, ax = plt.subplots(figsize = (12, 8))
ax.scatter(X[:,0], X[:,1])
plt.show()
这是一个非常紧密的聚类，几个值远离了聚类。 在这个简单的例子中，这些可以被认为是异常的。 为了弄清楚，我们正在为数据中的每个特征估计高斯分布。 为此，我们将创建一个返回每个要素的均值和方差的函数。
def estimate_gaussian(X): # INPUT：数据X # OUTPUT：数据的均值和方差 # TODO：根据数据计算均值和方差，高斯分布 # STEP：根据数据计算均值和方差 mu = np.mean(X, axis=0) #sigma = np.power(np.std(X, axis=0), 2) sigma = np.var(X, axis=0) return mu, sigma mu, sigma = estimate_gaussian(X) mu, sigma
Out[128]:
(array([14.11222578, 14.99771051]), array([1.83263141, 1.70974533]))
现在我们有了我们的模型参数，我们需要确定概率阈值，这表明一个样本应该被认为是一个异常。 为此，我们需要使用一组标记的验证数据（其中真实异常样本已被标记），并在给出不同阈值的情况下，对模型的性能进行鉴定。
Xval = data['Xval'] yval = data['yval'] ? print(type(Xval), Xval.shape, type(yval), yval.shape) print(Xval, yval)
我们还需要一种计算数据点属于正态分布的概率的方法。 幸运的是SciPy有这个内置的方法。
from scipy import stats dist = stats.norm(mu[0], sigma[0]) #实现正态分布即高斯分布，？换一种方式试试？ dist.pdf(x=15) #计算在x处的概率密度函数，默认均值为0，方差为1
Out[130]:
0.1935875044615038
我们还可以将数组传递给概率密度函数，并获得数据集中每个点的概率密度。
dist.pdf(X[:,0])[0:10]
Out[131]:
array([0.183842 , 0.20221694, 0.21746136, 0.19778763, 0.20858956, 0.21652359, 0.16991291, 0.15123542, 0.1163989 , 0.1594734 ])
我们定义高斯分布函数，计算并保存给定上述的高斯模型参数的数据集中每个值的概率密度。
def Gaussian(X, mu, sigma): px = (1 / (np.sqrt(2 * np.pi * sigma))) * np.exp(-np.power((X - mu), 2) / (2*sigma)) return px[:,0] * px[:,1] ? px = Gaussian(X, mu, sigma) print(pp[0:12])
[0.06470829 0.05030417 0.07245035 0.05031575 0.06368497 0.04245832
0.04790945 0.03651115 0.0186658 0.05068826 0.02651509 0.08471302]
我们还需要为验证集（使用相同的模型参数）执行此操作。 我们将使用与真实标签组合的这些概率来确定将数据点分配为异常的最佳概率阈值。
pxval = Gaussian(Xval, mu, sigma) pxval.min(), pxval.max(), pxval, len(pxval)
接下来，我们需要一个函数，找到给定概率密度值和真实标签的最佳阈值。 为了做到这一点，我们将为不同的epsilon值计算F1分数。 F1是真阳性，假阳性和假阴性的数量的函数。 方程式在练习文本中。
def select_threshold(pval, yval): # INPUT：交叉验证集数据pval，标签yval # OUTPUT：最好的参数epsilon，最高的得分F1 # TODO：寻找最好的参数epsilon，最高的得分F1 # STEP1：初始化变量 best_epsilon = 0 best_f1 = 0 f1 = 0 step = (pval.max() - pval.min()) / 1000 #设定步长 # STEP2：计算F1分数 for epsilon in np.arange(pval.min()+step, pval.max(), step): preds = pval < epsilon tp = yval[preds==1].sum() #真阳性样本数，预测为真实际为真 precision = tp / preds.sum() #查准率，在预测为真的样本里实际为真样本所占比例 recall = tp / yval.sum() #查全率，在实际为真的样本里预测为真样本所占比例 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(pxval, yval) #计算阈值和F1分数 epsilon, f1
输出：(8.990852779269492e-05, 0.8750000000000001)
最后，我们可以将阈值应用于数据集，并可视化结果。
# indexes of the values considered to be outliers
outliers = np.where(pp < epsilon) #获取小于阈值的pp索引号（X数据集）
#outliers = X[np.where(pp < epsilon), :] #获取小于阈值的异常点（X数据集）
outliers
#画出数据集X超过阈值的异常点
fig, ax = plt.subplots(figsize=(9,6))
ax.scatter(X[:,0], X[:,1])
ax.scatter(X[outliers[0],0], X[outliers[0],1], s=100, facecolors='none', edgecolors='r', marker='o')
?
#绘制等高线图
delta = 0.3 # 注意delta不能太小！！！否则会生成太多的数据，导致矩阵相乘会出现内存错误。
x = np.arange(0,30,delta)
y = np.arange(0,30,delta)
# 这部分要转化为X形式的坐标矩阵，也就是一列是横坐标，一列是纵坐标，
# 然后才能传入gaussian中求解得到每个点的概率值
xx, yy = np.meshgrid(x, y) #将x行重复，将y列重复
points = np.c_[xx.ravel(), yy.ravel()] # 按列合并，一列横坐标，一列纵坐标
z = Gaussian(points, mu, sigma)
z = z.reshape(xx.shape) #
cont_levels = [10**h for h in range(-20,0,3)]
plt.contour(xx, yy, z, cont_levels)
plt.title('Gaussian Contours',fontsize=16)
plt.show()
协同过滤?
推荐引擎使用基于项目和用户的相似性度量来检查用户的历史偏好，以便为用户可能感兴趣的新“事物”提供建议。在本练习中，我们将实现一种称为协作过滤的特定推荐系统算法，并将其应用于 电影评分的数据集。
我们首先加载并检查我们将要使用的数据。
data = loadmat('data/ex8_movies.mat')
data
Y是包含从1到5的等级的（数量的电影x数量的用户）数组.R是包含指示用户是否给电影评分的二进制值的“指示符”数组。 两者应该具有相同的维度。
Y = data['Y'] #用户对电影评分 R = data['R'] #用户是否评分 Y.shape, R.shape #((1682, 943), (1682, 943))
我们可以通过平均排序Y来评估电影的平均评级。
Y[1, np.where(R[1,:]==1)[0]].mean() #Out[215]:3.2061068702290076
我们还可以通过将矩阵渲染成图像来尝试“可视化”数据。 我们不能从这里收集太多，但它确实给我们了解用户和电影的相对密度。
fig, ax = plt.subplots(figsize=(12,12))
ax.imshow(Y)
ax.set_xlabel('Users')
ax.set_ylabel('Movies')
fig.tight_layout() #调整坐标轴和标签的位置，使其更清楚
plt.show()
接下来，我们将实施协同过滤的代价函数。 直觉上，“代价”是指一组电影评级预测偏离真实预测的程度。 代价方程在练习文本中给出。 它基于文本中称为X和Theta的两组参数矩阵。 这些“展开”到“参数”输入中，以便稍后可以使用SciPy的优化包。同时在代价和梯度计算中添加正则化。
def cost(params, Y, R, num_features, alambda): # INPUT：参数params，数据与标签Y,R，特征个数num_features，正则化参数alambda # OUTPUT：梯度grad，代价函数J # TODO：计算带正则化的梯度和代价函数 # STEP1：创建变量 Y = np.matrix(Y) # (1682, 943) R = np.matrix(R) # (1682, 943) num_movies = Y.shape[0] num_users = Y.shape[1] # STEP2：重新设置矩阵的维度 # your code here (appro ~ 2 lines) X = np.matrix(np.reshape(params[:num_movies * num_features], (num_movies, num_features))) Theta = np.matrix(np.reshape(params[num_movies * num_features:], (num_users, num_features))) # STEP3：初始化参数 J = 0 X_grad = np.zeros(X.shape) Theta_grad = np.zeros(Theta.shape) # STEP4：计算代价函数 error = np.multiply(X * Theta.T - Y, R) squared_error = np.power(error, 2) J = (1. / 2) * np.sum(squared_error) # STEP5：加入正则项 J = J + alambda / 2 * np.power(X, 2).sum() J = J + alambda / 2 * np.power(Theta, 2).sum() # STEP6：计算包含正则项的梯度 X_grad = error * Theta + alambda * X Theta_grad = error.T * X + alambda * Theta # STEP7：拉直梯度矩阵 grad = np.concatenate((np.ravel(X_grad), np.ravel(Theta_grad))) return J, grad J, grad = cost(params, Y_sub, R_sub, features, 1.5) J, grad
Out[250]:
(31.344056244274217, array([ -0.95596339, 6.97535514, -0.10861109, 0.60308088, 2.77421145, 0.25839822, 0.12985616, 4.0898522 , -0.89247334, 0.29684395, 1.06300933, 0.66738144, 0.60252677, 4.90185327, -0.19747928, -10.13985478, 2.10136256, -6.76563628, -2.29347024, 0.48244098, -2.99791422, -0.64787484, -0.71820673, 1.27006666, 1.09289758, -0.40784086, 0.49026541]))