吴恩达机器学习作业ex8:K 异常检测和推荐系统(Python实现)详细注释

文章目录

1 异常检测
- 1.1 高斯分布
- 1.2 估计高斯参数
- 1.3 选择阈值 ε
- 1.4 高维数据集
2 推荐系统
- 2.1 电影评分数据集
- 2.2 协作过滤学习算法
- - 2.2.1 协同过滤成本函数
  - 2.2.2 梯度协同过滤
  - 2.2.3 Regularized cost function
  - 2.2.4 正则梯度
- 2.3 学习电影推荐
- - 2.3.1 推荐
后记

1 异常检测

在本练习中，您将实施一种异常检测算法，以检测服务器计算机中的异常行为。特征是测量每台服务器的吞吐量（mb/s）和响应延迟（ms）。在服务器运行期间，您收集了 m = 307 个服务器行为的示例，因此有一个未标记的数据集 {x(1), . ，x(m)}。您怀疑这些示例中的绝大多数都是服务器正常运行的 “正常”（非异常）示例，但在这个数据集中也可能存在一些服务器行为异常的示例。首先，您将从一个二维数据集开始，以便直观地了解算法的运行情况。在该数据集上，您将拟合高斯分布，然后找到概率极低的值，这些值可被视为异常值。之后，您将把异常检测算法应用到一个更大、维度更多的数据集。这部分练习将使用 ex8.m。ex8.m 的第一部分将可视化数据集，如图 1 所示。
在这里插入图片描述

1.1 高斯分布

要进行异常检测，首先需要根据数据分布拟合一个模型。
给定训练集 {x(1),…,x(m)}（其中 x(i) ∈ Rn），您需要估计每个特征 xi 的高斯分布。对于每个特征 i = 1 … n，需要找到适合第 i 维数据 {x_i⁽¹⁾，…，x_i^(m)}（每个示例的第 i 维）的参数 µ_i和 σ²_i。
在这里插入图片描述

# 导入必要的库
%matplotlib inline
import pandas as pd  # 导入pandas库，用于数据处理
import numpy as np  # 导入numpy库，用于数值计算
from scipy.io import loadmat  # 导入scipy库中的loadmat函数，用于加载MAT文件
import matplotlib.pyplot as plt  # 导入matplotlib库中的pyplot模块，用于数据可视化# 加载数据文件
mat = loadmat('data/ex8data1.mat')
print(mat.keys())  # 打印加载的MAT文件的所有键
# 输出: dict_keys(['__header__', '__version__', '__globals__', 'X', 'Xval', 'yval'])# 提取数据
X = mat['X']  # 训练集特征数据
Xval, yval = mat['Xval'], mat['yval']  # 验证集特征数据和标签数据
X.shape, Xval.shape, yval.shape  # 打印数据集的形状
# 输出: ((307, 2), (307, 2), (307, 1))# 定义函数用于绘制数据
def plot_data():plt.figure(figsize=(8, 5))  # 创建一个大小为8x5英寸的图形plt.plot(X[:, 0], X[:, 1], 'bx')  # 绘制训练集数据，使用蓝色叉号表示# plt.scatter(Xval[:, 0], Xval[:, 1], c=yval.flatten(), marker='x', cmap='rainbow')  # 可选: 绘制验证集数据，使用彩色表示# 调用函数绘制数据
plot_data()

在这里插入图片描述

def gaussian(X, mu, sigma2):'''计算多元高斯分布概率密度函数。参数:X: 样本矩阵，形状为 (m, n)，其中 m 为样本数量，n 为特征数量mu: 均值向量，形状为 (n,)sigma2: 协方差矩阵或方差向量，形状为 (n, n) 或 (n,)输出:一个形状为 (m, ) 的向量，包含每个样本的概率密度值。'''# 获取样本数量 m 和特征数量 nm, n = X.shape# 如果 sigma2 是一个向量，将其转换为对角矩阵if np.ndim(sigma2) == 1:sigma2 = np.diag(sigma2)# 计算归一化系数 norm# norm 是 (2 * pi) 的 n/2 次方乘以 sigma2 行列式的平方根的倒数norm = 1. / (np.power((2 * np.pi), n / 2) * np.sqrt(np.linalg.det(sigma2)))# 初始化存储指数项的数组 expexp = np.zeros((m, 1))# 遍历每个样本for row in range(m):xrow = X[row]  # 取出第 row 个样本# 计算每个样本的指数项exp[row] = np.exp(-0.5 * ((xrow - mu).T).dot(np.linalg.inv(sigma2)).dot(xrow - mu))# 返回概率密度值return norm * exp

1.2 估计高斯参数

您可以使用以下公式估计第 i 个特征的参数 (µi，σ²i)。要估计平均值，可以使用
在这里插入图片描述
以及您将使用的变量：

您的任务是完成 estimateGaussian.m 中的代码。该函数将数据矩阵 X 作为输入，并应输出一个 n 维向量 mu（表示所有 n 个特征的均值）和另一个 n 维向量 sigma2（表示所有特征的方差）。您可以使用对每个特征和每个训练示例的 for 循环来实现这一功能（不过向量化实现可能更有效率；如果您愿意，也可以使用向量化实现）。请注意，在 Octave/MATLAB 中，计算 σ²i 时，var 函数（默认情况下）将使用 1/m-1，而不是 1/m。完成 estimateGaussian.m 中的代码后，ex8.m 的下一部分将可视化拟合高斯分布的轮廓。您应该会得到与图 2 类似的曲线图。从图中可以看出，大部分示例位于概率最高的区域，而异常示例则位于概率较低的区域。
在这里插入图片描述

def getGaussianParams(X, useMultivariate):"""计算数据集 X 的高斯分布参数。参数:X: 数据集，每行是一个 n 维数据点，形状为 (m, n)useMultivariate: 布尔值，是否使用多元高斯分布返回:mu: n 维向量，表示数据集的均值sigma2: 方差或协方差矩阵，根据 useMultivariate 决定"""# 计算数据集 X 的均值向量 mu，形状为 (n,)mu = X.mean(axis=0)# 如果使用多元高斯分布，计算协方差矩阵 sigma2if useMultivariate:# 计算协方差矩阵 sigma2# 公式: sigma2 = (X-mu).T @ (X-mu) / m# 这里的 len(X) 代表样本数量 msigma2 = ((X - mu).T @ (X - mu)) / len(X)else:# 如果不使用多元高斯分布，计算每个特征的方差 sigma2# 使用 numpy 的 var 方法，参数 ddof=0 表示除以样本数 msigma2 = X.var(axis=0, ddof=0)return mu, sigma2

def plotContours(mu, sigma2):"""画出高斯概率分布的等高线图，在三维中是一个上凸的曲面。投影到平面上则是一圈圈的等高线。参数:mu: 高斯分布的均值向量sigma2: 高斯分布的方差向量或协方差矩阵"""# 定义网格步长delta = 0.3  # 注意 delta 不能太小，否则会生成太多的数据，导致矩阵相乘出现内存错误。# 生成x和y坐标范围x = np.arange(0, 30, delta)y = np.arange(0, 30, delta)# 生成网格坐标矩阵xx, yy = np.meshgrid(x, y)# 将网格坐标转换为点集合，每个点包含x和y坐标points = np.c_[xx.ravel(), yy.ravel()]  # 按列合并，一列横坐标，一列纵坐标# 计算每个点的高斯概率值z = gaussian(points, mu, sigma2)# 将概率值转换为网格形式，便于绘制等高线图z = z.reshape(xx.shape)  # 这步骤不能忘# 定义等高线的级别，这里参考作业中的级别cont_levels = [10**h for h in range(-20, 0, 3)]# 绘制等高线图plt.contour(xx, yy, z, cont_levels)  # 这个levels是作业里面给的参考,或者通过求解的概率推出来。# 设置标题plt.title('Gaussian Contours', fontsize=16)

# 首先绘制单变量高斯分布的等高线：
plot_data()
useMV = False
plotContours(*getGaussianParams(X, useMV))# 然后绘制多变量高斯分布的等高线：
plot_data()
useMV = True
# *表示解包元组参数
plotContours(*getGaussianParams(X, useMV))

在这里插入图片描述

从上面的图可以看到，一元高斯模型仅在横向和纵向上有变化，而多元高斯模型在斜轴上也有相关变化，对应着特征间的相关关系。

1.3 选择阈值 ε

现在您已经估算出了高斯参数，您可以研究在这种分布下，哪些示例的概率非常高，哪些示例的概率非常低。低概率示例更有可能是我们数据集中的异常示例。确定哪些示例是异常示例的一种方法是根据交叉验证集选择一个阈值。在这部分练习中，您将使用交叉验证集上的 F1 分数来实现选择阈值 ε 的算法。为此，我们将使用交叉验证集 {(x⁽¹⁾cv , y⁽¹⁾cv ), … . ,(x^(mcv)cv , y^(mcv)cv )}，其中 y = 1 对应异常示例，y = 0 对应正常示例。对于每个交叉验证示例，我们将计算 p(x(i)cv )。所有这些概率的向量 p(x⁽¹⁾cv ), … . , p(x^(mcv)cv ) 会以向量 pval 的形式传递给 selectThreshold.m。相应的标签 y⁽¹⁾cv , … . 函数 selectThreshold.m 应返回两个值：第一个是所选阈值 ε。如果示例 x 的低概率 p(x) < ε，则认为该示例异常。函数还应该返回 F1 分数，它可以告诉您在给定阈值的情况下，找到地面实况异常点的效果如何。对于许多不同的 ε 值，您将通过计算当前阈值正确分类和错误分类的示例数量来计算 F1 分数：
在这里插入图片描述
计算精确度和召回率的方法是

其中

tp 是真阳性的数量：实际标签显示为异常，而我们的算法正确地将其归类为异常。
fp 是假阳性的数量：实际标签显示它不是异常点，但我们的算法却错误地将它归类为异常点。
fn 是假阴性的数量：实际标签显示为异常，但我们的算法却错误地将其归类为非异常。

def selectThreshold(yval, pval):"""选择最佳阈值 epsilon 以最大化 F1 分数参数：yval：验证集的真实标签，形状为 (m, 1)pval：验证集上的概率估计值，形状为 (m, 1)返回：bestF1：最佳 F1 分数bestEpsilon：最佳阈值 epsilon"""def computeF1(yval, pval_):"""计算给定阈值下的 F1 分数参数：yval：验证集的真实标签，形状为 (m, 1)pval_：阈值过滤后的概率估计值，形状为 (m, 1)返回：F1：当前阈值下的 F1 分数"""m = len(yval)tp = float(len([i for i in range(m) if pval_[i] and yval[i]]))  # 真阳性fp = float(len([i for i in range(m) if pval_[i] and not yval[i]]))  # 假阳性fn = float(len([i for i in range(m) if not pval_[i] and yval[i]]))  # 假阴性prec = tp / (tp + fp) if (tp + fp) else 0  # 精度rec = tp / (tp + fn) if (tp + fn) else 0  # 召回率F1 = 2 * prec * rec / (prec + rec) if (prec + rec) else 0  # F1 分数return F1epsilons = np.linspace(min(pval), max(pval), 1000)  # 生成 1000 个阈值bestF1, bestEpsilon = 0, 0for e in epsilons:pval_ = pval < e  # 将概率估计值与阈值比较，生成布尔向量thisF1 = computeF1(yval, pval_)  # 计算当前阈值下的 F1 分数if thisF1 > bestF1:  # 更新最佳 F1 分数和阈值bestF1 = thisF1bestEpsilon = ereturn bestF1, bestEpsilon  # 返回最佳 F1 分数和阈值

# 获取高斯分布参数：均值和方差
mu, sigma2 = getGaussianParams(X, useMultivariate=False)# 使用验证集和高斯分布参数计算概率估计值
pval = gaussian(Xval, mu, sigma2)# 使用验证集上的真实标签和概率估计值选择最佳阈值
bestF1, bestEpsilon = selectThreshold(yval, pval)  
# 输出：最佳F1分数和对应的阈值 (0.8750000000000001, 8.999852631901397e-05)# 使用原始数据集和高斯分布参数计算每个样本的概率
y = gaussian(X, mu, sigma2)  # 找到概率小于最佳阈值的样本（离群点）
xx = np.array([X[i] for i in range(len(y)) if y[i] < bestEpsilon])
xx  # 离群点# 可视化原始数据和高斯分布的等高线
plot_data()
plotContours(mu, sigma2)# 在图中标出离群点，用红色圈表示
plt.scatter(xx[:,0], xx[:,1], s=80, facecolors='none', edgecolors='r')

在这里插入图片描述

1.4 高维数据集

脚本 ex8.m 的最后一部分将在更现实、更困难的数据集上运行您实施的异常检测算法。运行异常检测算法。在这个数据集中，每个示例由 11 个特征描述，捕获了计算服务器的更多属性。脚本将使用您的代码来估计高斯参数（μi 和 σ²i），评估训练数据 X 的概率（您根据训练数据 X 估计了高斯参数），并对交叉验证集 Xval 进行评估。最后，它将使用 selectThreshold 找到最佳阈值 ε。你应该会看到 epsilon 值约为 1.38e-18，并发现 117 个异常点。

# 加载数据集
mat = loadmat('/data/ex8data2.mat')
# X2 是训练集，包含 1000 个样本，每个样本有 11 个特征。
# Xval2 和 yval2 分别是交叉验证集的特征和标签。
X2 = mat['X']
Xval2, yval2 = mat['Xval'], mat['yval']
X2.shape  # (1000, 11)
# 计算高斯分布的参数
mu, sigma2 = getGaussianParams(X2, useMultivariate=False)
# 计算训练集每个样本的概率
ypred = gaussian(X2, mu, sigma2)
# 计算交叉验证集每个样本的概率
yval2pred = gaussian(Xval2, mu, sigma2)
# 使用交叉验证集选择最佳的阈值 epsilon
bestF1, bestEps = selectThreshold(yval2, yval2pred)
# 找出训练集中所有概率值小于最佳阈值 epsilon 的样本（即异常点）
anoms = [X2[i] for i in range(X2.shape[0]) if ypred[i] < bestEps]
bestEps, len(anoms)
# (1.378607498200024e-18, 117)

2 推荐系统

在这部分练习中，您将实施协同过滤学习算法，并将其应用于一个电影评分数据集2。该数据集有 n_u = 943 个用户，n_m = 1682 部电影。在这部分练习中，您将使用脚本 ex8 cofi.m。在接下来的练习中，您将实现计算协同拟合目标函数和梯度的函数 cofiCostFunc.m。在实现代价函数和梯度后，您将使用 fmincg.m 学习协同过滤的参数。

2.1 电影评分数据集

脚本 ex8 cofi.m 的第一部分将加载数据集 ex8 movies.mat，并在 Octave/MATLAB 环境中提供变量 Y 和 R。矩阵 Y（num movies × num users matrix）存储评分 y(i,j)（从 1 到 5）。矩阵 R 是二值指示矩阵，如果用户 j 给电影 i 打了分，则 R(i, j) = 1，否则 R(i, j) = 0。协同过滤的目的是预测用户尚未评分的电影的评分，即 R(i, j) = 0 的条目。为了帮助您理解矩阵 Y，脚本 ex8 cofi.m 将计算第一部电影（《玩具总动员》）的平均电影评分，并将平均评分输出到屏幕上：
在这里插入图片描述
X 的第 i 行对应第 i 部电影的特征向量 x(i)，Theta 的第 j 行对应第 j 个用户的参数向量 θ(j)。x(i) 和 θ(j) 都是 n 维向量。本练习将使用 n = 100，因此 x(i) ∈ R¹⁰⁰，θ(j) ∈ R¹⁰⁰。相应地，X 是 n_m× 100 矩阵，Theta 是 n_u × 100 矩阵。

# 加载数据集
mat = loadmat('/data/ex8_movies.mat')
print(mat.keys())
# dict_keys(['__header__', '__version__', '__globals__', 'Y', 'R'])
'''
Y 是评分矩阵，R 是指示矩阵，表明用户是否对电影评分。
nm 是电影数量，nu 是用户数量，nf 是特征数量（假设为100）。
Y.shape, R.shape 打印评分矩阵和指示矩阵的形状，均为 (1682, 943)。
'''
Y, R = mat['Y'], mat['R']
nm, nu = Y.shape  # 获取电影数量和用户数量
nf = 100  # 特征数量
Y.shape, R.shape
# ((1682, 943), (1682, 943))
Y[0].sum() / R[0].sum()  # 分子代表第一个电影的总分数，分母代表这部电影有多少评分数据
# 可视化评分矩阵
fig = plt.figure(figsize=(8, 8*(1682./943.)))  # 设置图形大小
plt.imshow(Y, cmap='rainbow')  # 显示评分矩阵，使用 'rainbow' 颜色映射
plt.colorbar()  # 添加颜色条
plt.ylabel('Movies (%d)' % nm, fontsize=20)  # 添加 y 轴标签
plt.xlabel('Users (%d)' % nu, fontsize=20)  # 添加 x 轴标签
plt.show()  # 显示图形

2.2 协作过滤学习算法

现在，您将开始执行协同过滤学习算法。在电影推荐设置中，协同过滤算法考虑了一组 n 维参数向量 x(1), …, x(nm) 和 θ(1), …, θ(nu)，其中模型预测用户 j 对电影 i 的评分为 y(i,j) = (θ(j))T x(i)。给定的数据集由一些用户对一些电影的评分组成，您希望学习能产生最佳拟合效果（使平方误差最小）的参数向量 x(1)、…、x(nm)、θ(1)、…、θ(nu)。您将完成 cofiCostFunc.m 中的代码，以计算协同过滤的成本函数和梯度。请注意，函数的参数（即您试图学习的值）是 X 和 Theta。为了使用 fmincg 等现成的最小化器，我们设置了代价函数，将参数展开为单个向量 params。您之前在神经网络编程练习中使用过相同的向量展开方法。

2.2.1 协同过滤成本函数

协同过滤成本函数（无正则化）为
在这里插入图片描述
请注意，只有当 R(i, j) = 1 时，您才会累计用户 j 和影片 i 的成本。完成函数后，脚本 ex8 cofi.m 将运行您的成本函数。您应该会看到 22.22 的输出结果。

# 加载电影参数数据集
mat = loadmat('data/ex8_movieParams.mat')
X = mat['X']
Theta = mat['Theta']
'''
X 是电影特征矩阵，每行代表一个电影，每列代表一个特征。
Theta 是用户参数矩阵，每行代表一个用户，每列代表一个特征。
nu 是用户数量，nm 是电影数量，nf 是特征数量。
'''
nu = int(mat['num_users'])
nm = int(mat['num_movies'])
nf = int(mat['num_features'])
# 为了更快运行，缩小数据集大小
nu = 4  # 将用户数量减少到4
nm = 5  # 将电影数量减少到5
nf = 3  # 将特征数量减少到3# 调整数据集大小
X = X[:nm, :nf]  # 取前5部电影的前3个特征
Theta = Theta[:nu, :nf]  # 取前4个用户的前3个特征
Y = Y[:nm, :nu]  # 取前5部电影和前4个用户的评分
R = R[:nm, :nu]  # 取前5部电影和前4个用户的评分指示矩阵
X.shape, Theta.shape
#((5, 3), (4, 3))

def serialize(X, Theta):'''展开参数'''# 将矩阵X和Theta展开成一维数组，并将它们连接在一起return np.r_[X.flatten(), Theta.flatten()]def deserialize(seq, nm, nu, nf):'''提取参数'''# 从一维数组中提取参数，重塑为原来的矩阵形状X = seq[:nm*nf].reshape(nm, nf)Theta = seq[nm*nf:].reshape(nu, nf)return X, Thetadef cofiCostFunc(params, Y, R, nm, nu, nf, l=0):"""params : 拉成一维之后的参数向量(X, Theta)Y : 评分矩阵 (nm, nu)R ：0-1矩阵，表示用户对某一电影有无评分nu : 用户数量nm : 电影数量nf : 自定义的特征的维度l : lambda for regularization"""# 反序列化params，恢复为矩阵X和ThetaX, Theta = deserialize(params, nm, nu, nf)# (X @ Theta.T) 计算用户对电影的预测评分# (X @ Theta.T - Y) 计算预测评分和实际评分之间的误差# (X @ Theta.T - Y) * R 将误差矩阵中没有评分的数据置为0error = 0.5 * np.square((X @ Theta.T - Y) * R).sum()# 计算正则化项reg1 = 0.5 * l * np.square(Theta).sum()reg2 = 0.5 * l * np.square(X).sum()# 返回代价函数的值return error + reg1 + reg2

cofiCostFunc(serialize(X,Theta),Y,R,nm,nu,nf),cofiCostFunc(serialize(X,Theta),Y,R,nm,nu,nf,1.5)
# (22.224603725685675, 31.34405624427422)

2.2.2 梯度协同过滤

现在，您应该实现梯度计算（无正则化）。具体来说，您应该完成 cofiCostFunc.m 中的代码，以返回变量 X grad 和 Theta grad。请注意，X grad 应该是与 X 大小相同的矩阵，同样，Theta grad 也是与 Theta 大小相同的矩阵。代价函数的梯度取值为
在这里插入图片描述
请注意，该函数通过将两组变量展开为一个单一向量来返回它们的梯度。完成梯度计算代码后，脚本 ex8 cofi.m 将运行梯度检查（checkCostFunction），对梯度的实现进行数值检查。


def cofiGradient(params, Y, R, nm, nu, nf, l=0):"""计算X和Theta的梯度，并序列化输出。"""# 反序列化params，恢复为矩阵X和ThetaX, Theta = deserialize(params, nm, nu, nf)# 计算X的梯度# (X @ Theta.T - Y) * R 计算误差并掩盖掉没有评分的数据# (误差矩阵 @ Theta) 加上正则化项 l * XX_grad = ((X @ Theta.T - Y) * R) @ Theta + l * X# 计算Theta的梯度# (X @ Theta.T - Y) * R 计算误差并掩盖掉没有评分的数据# 转置误差矩阵 (误差矩阵.T @ X) 加上正则化项 l * ThetaTheta_grad = ((X @ Theta.T - Y) * R).T @ X + l * Theta# 序列化梯度并返回return serialize(X_grad, Theta_grad)

2.2.3 Regularized cost function

正则化协同过滤的成本函数为
在这里插入图片描述
完成后，脚本 ex8 cofi.m 将运行您的正则化成本函数，您将看到约 31.34 的成本。

def checkGradient(params, Y, myR, nm, nu, nf, l=0.):"""Let's check my gradient computation real quick"""print('Numerical Gradient \t cofiGrad \t\t Difference')# 分析出来的梯度grad = cofiGradient(params, Y, myR, nm, nu, nf, l)# 用微小的e来计算数值梯度。e = 0.0001nparams = len(params)e_vec = np.zeros(nparams)# Choose 10 random elements of param vector and compute the numerical gradient# 每次只能改变e_vec中的一个值，并在计算完数值梯度后要还原。for i in range(10):idx = np.random.randint(0, nparams)e_vec[idx] = e# 计算数值梯度loss1 = cofiCostFunc(params - e_vec, Y, myR, nm, nu, nf, l)loss2 = cofiCostFunc(params + e_vec, Y, myR, nm, nu, nf, l)numgrad = (loss2 - loss1) / (2 * e)# 恢复e_vec的值e_vec[idx] = 0# 计算差异diff = np.linalg.norm(numgrad - grad[idx]) / np.linalg.norm(numgrad + grad[idx])# 输出数值梯度、分析梯度和差异print('%0.15f \t %0.15f \t %0.15f' % (numgrad, grad[idx], diff))

在这里插入图片描述

2.2.4 正则梯度

现在，您已经实现了正则化代价函数，接下来应该对梯度进行正则化。您应该在 cofiCostFunc.m 中添加正则化项，以返回正则化梯度。请注意，正则化代价函数的梯度由以下公式给出：
在这里插入图片描述
这意味着您只需将 λx(i) 添加到前面描述的 X grad(i,:) 变量中，并将λθ(j) 添加到前面描述的 Theta grad(j,:) 变量中。在您完成计算梯度的代码后，脚本 ex8 cofi.m 将运行另一个梯度检查（checkCostFunction），对梯度的实现进行数值检查。

'''
输出结果包括解析梯度、数值梯度和两者之间的差异。
如果解析梯度和数值梯度的差异很小（接近于0），则说明梯度计算是正确的。
'''
print("Checking gradient with lambda = 0...")
checkGradient(serialize(X,Theta), Y, R, nm, nu, nf)
print("\nChecking gradient with lambda = 1.5...")
checkGradient(serialize(X,Theta), Y, R, nm, nu, nf, l=1.5)

在这里插入图片描述

2.3 学习电影推荐

完成协同过滤成本函数和梯度的实现后，现在就可以开始训练算法，为自己推荐电影了。在 ex8 cofi.m 脚本的下一部分，你可以输入自己的电影偏好，这样以后算法运行时，你就可以得到自己的电影推荐了！我们根据自己的喜好填写了一些值，但您也可以根据自己的口味进行修改。所有电影的列表及其在数据集中的编号可以在文件 movie idx.txt 中找到。

# 初始化电影列表
movies = []  # 存储所有电影的名称
with open('data/movie_ids.txt', 'r', encoding='latin-1') as f:for line in f:# 读取每行并去除行首行尾的空白符，通过空格分割后，仅提取电影名称部分movies.append(' '.join(line.strip().split(' ')[1:]))# 初始化用户评分矩阵
my_ratings = np.zeros((1682, 1))  # 初始化一个 1682x1 的零矩阵，表示1682部电影的评分# 添加用户对电影的评分
my_ratings[0]   = 4   # Toy Story
my_ratings[97]  = 2   # Silence of the Lambs, The (1991)
my_ratings[6]   = 3   # Twelve Monkeys
my_ratings[11]  = 5   # Usual Suspects
my_ratings[53]  = 4   # Outbreak
my_ratings[63]  = 5   # Shawshank Redemption
my_ratings[65]  = 3   # While You Were Sleeping
my_ratings[68]  = 5   # Forrest Gump
my_ratings[182] = 4   # Alien
my_ratings[225] = 5   # Die Hard 2
my_ratings[354] = 5   # Sphere# 输出用户评分过的电影及其评分
for i in range(len(my_ratings)):if my_ratings[i] > 0:  # 过滤出用户评分过的电影print(my_ratings[i], movies[i])# 加载电影数据
mat = loadmat('data/ex8_movies.mat')  # 读取包含电影评分数据的MAT文件
Y, R = mat['Y'], mat['R']  # 提取电影评分矩阵Y和评分标记矩阵R
Y.shape, R.shape  # 输出矩阵的形状，分别为(1682, 943)# 添加用户的评分数据到矩阵中
Y = np.c_[Y, my_ratings]  # 将用户的评分数据添加到电影评分矩阵Y的最后一列，形状变为(1682, 944)
R = np.c_[R, my_ratings != 0]  # 更新评分标记矩阵R，形状变为(1682, 944)# 重新定义参数
nm, nu = Y.shape  # 更新后的电影评分矩阵Y的形状，表示1682部电影和944个用户
nf = 10  # 使用10维的特征向量表示电影和用户的特征

在这里插入图片描述

def normalizeRatings(Y, R):"""归一化电影评分矩阵Y，计算每部电影的平均评分，并用这个平均评分来调整评分矩阵。参数：Y - 电影评分矩阵 (nm, nu)R - 评分标记矩阵 (nm, nu)，0-1矩阵，表示用户对某一电影有无评分输出：Ynorm - 归一化后的评分矩阵 (nm, nu)Ymean - 每部电影的平均评分 (nm, 1)"""# 计算每部电影的平均评分，只有那些评分过的电影才参与计算Ymean = (Y.sum(axis=1) / R.sum(axis=1)).reshape(-1,1)# 对评分进行归一化处理，未评分的电影不参与归一化Ynorm = (Y - Ymean) * Rreturn Ynorm, Ymean

# 归一化后的评分矩阵和平均评分的形状
Ynorm, Ymean = normalizeRatings(Y, R)
print(Ynorm.shape, Ymean.shape)
# 输出：((1682, 944), (1682, 1))# 随机初始化电影特征矩阵X和用户特征矩阵Theta
X = np.random.random((nm, nf))
Theta = np.random.random((nu, nf))# 序列化参数
params = serialize(X, Theta)# 正则化参数
l = 10

import scipy.optimize as opt# 使用优化函数来最小化协同过滤的成本函数
res = opt.minimize(fun=cofiCostFunc,       # 目标函数，计算协同过滤的成本x0=params,              # 初始参数，展开后的X和Thetaargs=(Y, R, nm, nu, nf, l),  # 传递给目标函数的额外参数method='TNC',           # 使用的方法，这里选择的是信赖区域牛顿共轭梯度法jac=cofiGradient,       # 梯度函数，计算协同过滤的梯度options={'maxiter': 100}  # 优化选项，这里设置最大迭代次数为100
)# 最优参数
ret = res.x# 反序列化参数，恢复为电影特征矩阵和用户特征矩阵
fit_X, fit_Theta = deserialize(ret, nm, nu, nf)

2.3.1 推荐

将额外的评分添加到数据集后，脚本将继续训练协同过滤模型。这将学习参数 X 和 Theta。要预测用户 j 对电影 i 的评分，需要计算 (θ(j))T x(i)。脚本的下一部分将计算所有电影和用户的评分，并根据脚本前面输入的评分显示其推荐的电影（图 4）。请注意，由于随机初始化的不同，您可能会得到一组不同的预测结果。在这里插入图片描述

# 所有用户的剧场分数矩阵
pred_mat = fit_X @ fit_Theta.T
# 最后一个用户的预测分数， 也就是我们刚才添加的用户
pred = pred_mat[:,-1] + Ymean.flatten()
# 排序并翻转，使之从大到小排列
pred_sorted_idx = np.argsort(pred)[::-1]
print(pred_mat)
print("")
print(pred)
print("")
print(pred_sorted_idx)

在这里插入图片描述

# 打印前10个推荐的电影
for i in range(10):print('Predicting rating %0.1f for movie %s.' \%(pred[pred_sorted_idx[i]], movies[pred_sorted_idx[i]]))# 打印原始评分
print("\nOriginal ratings provided:")
for i in range(len(my_ratings)):if my_ratings[i] > 0:print('Rated %d for movie %s.' % (my_ratings[i], movies[i]))