吴恩达机器学习作业八：异常检测和推荐系统
在本练习中，您将实现异常检测算法，并将其应用于检测网络上出现故障的服务器。在第二部分中，您将使用协作过滤来构建电影推荐系统。
1 Anomaly detection
在本练习中，您将实现异常检测算法来检测服务器计算机中的异常行为。这些特性测量每个服务器的吞吐量（mb/s）和响应延迟（ms）。当您的服务器运行时，您收集了m=307个它们的行为示例，因此有一个未标记的数据集{x（1），x（m）}。您怀疑这些示例中的绝大多数都是服务器正常运行的“正常”（非异常）示例，但也可能有一些服务器在此数据集中异常运行的示例。您将使用高斯模型来检测数据集中的异常示例。您将首先从一个2D数据集开始，该数据集将允许您可视化算法正在执行的操作。在该数据集上，您将拟合高斯分布，然后找到概率非常低的值，因此可以认为是异常。然后，将异常检测算法应用于具有多个维度的较大数据集。
path = "D:\编程\ex8data1.mat"
data = loadmat(path)
print(data.keys())
#dict_keys(['__header__', '__version__', '__globals__', 'X', 'Xval', 'yval'])
X = data['X']
Xval = data['Xval']
yval = data['yval']
print(X.shape,Xval.shape,yval.shape)
#(307, 2) (307, 2) (307, 1)
def plotData(X):
plt.figure(figsize = (5,5),dpi = 80)
plt.scatter(X[:,0],X[:,1],marker='x',c = 'b')
plt.xlabel('Latency (ms)')
plt.ylabel('Throughput (mb/s)')
plt.title('The first dataset')
plotData(X)
plt.show()
1.1 Gaussian distribution
要执行异常检测，首先需要根据数据的分布调整模型。
给定一个训练集{x（1），…，x（m）}（其中x（i）∈ Rn），你要估计每个特征的高斯分布。对于每个特征i=1…n，你需要找到参数µ和σ2将数据拟合到第i维{x（1）i，…，x（m）i}（每个例子的第i维）。
高斯分布由下式给出：
其中，µ是指平均值，σ2控制方差。
def Gaussian_distribution(X,mu,sigma2):
if np.ndim(sigma2) == 1:
sigma2 = np.diag(sigma2)
m,n = X.shape
X = X - mu
first = np.power(2*np.pi,-n/2)*np.power(np.linalg.det(sigma2),-1/2)
second = np.exp(-1/2*X@np.linalg.inv(sigma2)@X.T)
#(307,307)
second = np.diag(second)
#(307,)
p = first*second
p = p.reshape(-1,1)
# 以下为另一种方法：
# m = len(X)
# second = np.zeros((m,1))
# first = np.power(2 * np.pi, -n / 2) * np.power(np.linalg.det(sigma2), -1 / 2)
# for row in range(m):
# second[row] =np.exp(X[row].T@np.linalg.inv(sigma2)@X[row])
# # （1,2)*(2,2)*(2,1)
# p = first*second
return p
1.2 Estimating parameters for a Gaussian
您的任务是完成代码estimateGaussian。此函数将数据矩阵X作为输入，并应输出一个n维向量mu（保存所有n个特征的平均值）和另一个n维向量sigma2（保存所有特征的方差）。您可以在每个特性和每个训练示例上使用for循环来实现这一点（尽管向量化实现可能更有效；如果您愿意，可以随意使用矢量化实现）。
一旦您完成了estimateGaussian代码，下一部分将可视化拟合高斯分布的轮廓。您应该得到下图。从你的图中，你可以看到大多数的例子在概率最高的区域，而异常的例子在概率较低的区域。
# 计算均值和方差
def estimateGaussian(X,isCovariance):
mu = np.mean(X,axis = 0)
# 如果具有相关性，则计算协方差；反之计算方差
if isCovariance :
sigma2 = ((X-mu).T@(X-mu))/len(X)
else:
sigma2 = np.var(X,axis = 0,ddof=0)
return mu,sigma2
def plotContours(mu,sigma2):
"""
画出高斯概率分布的图，在三维中是一个上凸的曲面。投影到平面上则是一圈圈的等高线。
"""
delta = 0.3
x = np.arange(0,30,delta)
y = np.arange(0,30,delta)
xx,yy = np.meshgrid(x,y)
link = np.c_[xx.ravel(),yy.ravel()]
z = Gaussian_distribution(link,mu,sigma2)
z = z.reshape(xx.shape)
level = [10**h for h in range(-20,0,3)]
plt.contour(xx,yy,z,level)
plt.title('Gaussion Contours')
运行入口：
if __name__ =="__main__":
# 有相关性
mu,sigma2 = estimateGaussian(X,isCovariance=True)
plotData(X)
plotContours(mu,sigma2)
plt.show()
if __name__ =="__main__":
# 无相关性
mu,sigma2 = estimateGaussian(X,isCovariance=False)
plotData(X)
plotContours(mu,sigma2)
plt.show()
从上面的图可以看到，一元高斯模型仅在横向和纵向上有变化，而多元高斯模型在斜轴上也有相关变化，对应着特征间的相关关系。而一元高斯模型就是多元高斯模型中协方差矩阵为对角矩阵的结果，即协方差都为0，不考虑协方差，只考虑方差，故一元高斯模型不会有斜轴上的变化。
从上面的图我们可以清晰的看到，哪些样本的概率高，哪些样本的概率低，概率低的样本很大程度上就是异常值。
1.3 Selecting the threshold, ε
现在您已经估计了高斯参数，您可以研究哪些示例在给定此分布的情况下具有非常高的概率，哪些示例具有非常低的概率。低概率的例子更可能是我们数据集中的异常。确定哪些示例是异常的一种方法是基于交叉验证集选择阈值。在这部分练习中，您将实现一个选择阈值的算法ε 使用交叉验证集上的分数。
现在应该完成代码selectThreshold。为此，我们将使用交叉验证集
，其中标签y=1对应于异常示例，y=0对应于正常示例。对于每个交叉验证示例，我们将计算p（x（i）cv）。所有这些概率的向量p（x（1）cv），p（x（mcv）cv）传递给向量pval中的selectThreshold。相应的标签y（1）cv，y（mcv）cv传递给向量yval中的相同函数。
函数selectThreshold应该返回两个值；第一个是选定的阈值ε. 如果示例x具有低概率p（x）<ε, 那么它就被认为是一个异常。函数还应该返回F1score，它告诉您在给定某个阈值的情况下，您在查找地面真值异常方面做得有多好。对于许多不同的值ε, 您将通过计算当前阈值正确分类和错误分类的示例数来计算结果分数。
在代码selectThreshold中，已经有一个循环将尝试ε 来选择最好的ε 根据分数。现在应该完成代码selectThreshold。您可以在所有交叉验证示例上使用for循环来实现F1分数的计算（以计算值tp、fp、fn）。你应该看到一个大约8.99e-05的ε值。
def selectThreshold(yval,pval):
bestEpsilon = 0
bestF1 = 0
epsilons = np.linspace(min(pval),max(pval),1000)
for e in epsilons:
pval_ = pval < e
thisF1 = computeF1(yval,pval_)
if thisF1 > bestF1:
bestF1 = thisF1
bestEpsilon = e
print(pval_)
return bestF1,bestEpsilon
def computeF1(yval,pval):
# pval为预测值，yval是实际值
length = len(yval)
#也可以这么写tp = np.sum((yval == 1) & (pval == 1))
tp = float(len([i for i in range(length) if yval[i] and pval[i]]))
fp =float(len([j for j in range(length) if not yval[j] and pval[j]]))
fn =float(len([k for k in range(length) if yval[k] and not pval[k]]))
# 以下由于做除法，分母不能等于0
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
return F1
if __name__ =="__main__":
mu,sigma2 = estimateGaussian(X,isCovariance=False)
pval = Gaussian_distribution(Xval,mu,sigma2)
bestF1,bestEpsilon = selectThreshold(yval,pval)
print("bestF1 = {}".format(bestF1))
print("bestEpsilon = {}".format(bestEpsilon))
#bestF1 = 0.8750000000000001
#bestEpsilon = [8.99985263e-05]
def plotAnomalies(X,bestEpsilon):
# 画出异常点
p = Gaussian_distribution(X, mu, sigma2) # X的概率
xx = np.array([X[i] for i in range(len(p)) if p[i] < bestEpsilon])# 离群点
plt.scatter(xx[:,0],xx[:,1],marker='o',edgecolors = 'r',facecolors='none')
if __name__ =="__main__":
mu,sigma2 = estimateGaussian(X,isCovariance=False)
pval = Gaussian_distribution(Xval,mu,sigma2)
bestF1,bestEpsilon = selectThreshold(yval,pval)
plotData(X)
plotContours(mu,sigma2)
plotAnomalies(X, bestEpsilon)
plt.show()
1.4 High dimensional dataset
最后一部分将运行您在更真实、更困难的数据集上实现的异常检测算法。在这个数据集中，每个示例都由11个特性描述，它们捕获了计算服务器的更多属性。
脚本将使用您的代码来估计高斯参数(µi,σ2 i），评估用于估计高斯参数的训练数据X的概率，并对交叉验证集Xval进行评估。最后，它将使用selectThreshold来找到最佳阈值ε. 你应该看到ε值约为1.38e-18，发现117个异常。
path2 = "D:\编程\ex8data2.mat"
data2 = loadmat(path2)
print(data2.keys())
X2 = data2['X']
Xval2 = data2['Xval']
yval2 = data2['yval']
print(X2.shape,Xval2.shape,yval2.shape)
if __name__ =="__main__":
mu,sigma2 = estimateGaussian(X2,isCovariance=False)
pred = Gaussian_distribution(X2,mu,sigma2)
pval2 = Gaussian_distribution(Xval2,mu,sigma2)
bestF1,bestEpsilon = selectThreshold(yval2,pval2)
# 用测试集训练出bestF1,bestEpsilon，然后将其应用于训练集，计算出异常点个数
Anomalies = [X2[i] for i in range(len(pred)) if pred[i] < bestEpsilon]
# 离群点
print(bestEpsilon,len(Anomalies))
# (1.378607498200024e-18, 117)
完整代码：
#encoding=utf-8
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.io import loadmat
path = "D:\编程\ex8data1.mat"
data = loadmat(path)
#print(data.keys())
X = data['X']
Xval = data['Xval']
yval = data['yval']
print(X.shape,Xval.shape,yval.shape)
def plotData(X):
plt.figure(figsize = (5,5),dpi = 80)
plt.scatter(X[:,0],X[:,1],marker='x',c = 'b')
plt.xlabel('Latency (ms)')
plt.ylabel('Throughput (mb/s)')
plt.title('The first dataset')
#plotData(X)
#plt.show()
def Gaussian_distribution(X,mu,sigma2):
if np.ndim(sigma2) == 1:
sigma2 = np.diag(sigma2)
m,n = X.shape
X = X - mu
first = np.power(2*np.pi,-n/2)*np.power(np.linalg.det(sigma2),-1/2)
second = np.exp(-1/2*X@np.linalg.inv(sigma2)@X.T)
#(307,307)
second = np.diag(second)
#(307,)
p = first*second
p = p.reshape(-1,1)
# 以下为另一种方法：
# m = len(X)
# second = np.zeros((m,1))
# first = np.power(2 * np.pi, -n / 2) * np.power(np.linalg.det(sigma2), -1 / 2)
# for row in range(m):
# second[row] =np.exp(X[row].T@np.linalg.inv(sigma2)@X[row])
# # （1,2)*(2,2)*(2,1)
# p = first*second
return p
def estimateGaussian(X,isCovariance):
mu = np.mean(X,axis = 0)
if isCovariance :
sigma2 = ((X-mu).T@(X-mu))/len(X)
else:
sigma2 = np.var(X,axis = 0,ddof=0)
return mu,sigma2
# if __name__ == "__main__":
# mu,sigma2 = estimateGaussian(X,isCovariance = False)
# print(mu,sigma2)
# p = Gaussian_distribution(X,mu,sigma2)
# print(p)
def plotContours(mu,sigma2):
"""
画出高斯概率分布的图，在三维中是一个上凸的曲面。投影到平面上则是一圈圈的等高线。
"""
delta = 0.3
x = np.arange(0,30,delta)
y = np.arange(0,30,delta)
xx,yy = np.meshgrid(x,y)
link = np.c_[xx.ravel(),yy.ravel()]
z = Gaussian_distribution(link,mu,sigma2)
z = z.reshape(xx.shape)
level = [10**h for h in range(-20,0,3)]
plt.contour(xx,yy,z,level)
plt.title('Gaussion Contours')
# if __name__ =="__main__":
# mu,sigma2 = estimateGaussian(X,isCovariance=False)
# plotData(X)
# plotContours(mu,sigma2)
# plt.show()
def selectThreshold(yval,pval):
bestEpsilon = 0
bestF1 = 0
epsilons = np.linspace(min(pval),max(pval),1000)
for e in epsilons:
pval_ = pval < e
thisF1 = computeF1(yval,pval_)
if thisF1 > bestF1:
bestF1 = thisF1
bestEpsilon = e
return bestF1,bestEpsilon
def computeF1(yval,pval):
# pval为预测值，yval是实际值
length = len(yval)
#也可以这么写tp = np.sum((yval == 1) & (pval == 1))
tp = float(len([i for i in range(length) if yval[i] and pval[i]]))
fp =float(len([j for j in range(length) if not yval[j] and pval[j]]))
fn =float(len([k for k in range(length) if yval[k] and not pval[k]]))
# 以下由于做除法，分母不能等于0
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
return F1
# if __name__ =="__main__":
# mu,sigma2 = estimateGaussian(X,isCovariance=False)
# pval = Gaussian_distribution(Xval,mu,sigma2)
# bestF1,bestEpsilon = selectThreshold(yval,pval)
# print("bestF1 = {}".format(bestF1))
# print("bestEpsilon = {}".format(bestEpsilon))
# def plotAnomalies(X,bestEpsilon):
# p = Gaussian_distribution(X, mu, sigma2) # X的概率
# xx = np.array([X[i] for i in range(len(p)) if p[i] < bestEpsilon])# 离群点
# plt.scatter(xx[:,0],xx[:,1],marker='o',edgecolors = 'r',facecolors='none')
# if __name__ =="__main__":
# mu,sigma2 = estimateGaussian(X,isCovariance=False)
# pval = Gaussian_distribution(Xval,mu,sigma2)
# bestF1,bestEpsilon = selectThreshold(yval,pval)
# plotData(X)
# plotContours(mu,sigma2)
# plotAnomalies(X, bestEpsilon)
# plt.show()
path2 = "D:\编程\ex8data2.mat"
data2 = loadmat(path2)
print(data2.keys())
X2 = data2['X']
Xval2 = data2['Xval']
yval2 = data2['yval']
print(X2.shape,Xval2.shape,yval2.shape)
if __name__ =="__main__":
mu,sigma2 = estimateGaussian(X2,isCovariance=False)
pred = Gaussian_distribution(X2,mu,sigma2)
pval2 = Gaussian_distribution(Xval2,mu,sigma2)
bestF1,bestEpsilon = selectThreshold(yval2,pval2)
Anomalies = [X2[i] for i in range(len(pred)) if pred[i] < bestEpsilon] # 离群点
print(bestEpsilon,len(Anomalies))
2 Recommender Systems
在本练习的这一部分中，您将实现协作过滤学习算法，并将其应用于电影分级数据集。此数据集由分级为1到5的分级组成。数据集有nu=943个用户，nm=1682个电影。在本部分练习中，您将使用脚本ex8cofi。在本练习的下一部分中，您将实现计算协作拟合目标函数和梯度的函数cofiCostFunc。在实现了代价函数和梯度之后，您将使用fmincg来学习协作过滤的参数。
2.1 Movie ratings dataset
# 推荐系统
#encoding=utf-8
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.io import loadmat
path = "D:\编程\ex8_movies.mat"
data = loadmat(path)
print(data.keys())
Y = data['Y']
R = data['R']
nm,nu = Y.shape
print(Y.shape,R.shape)
print(nm,nu)
plt.figure(figsize=(6,6))
plt.imshow(Y,cmap='rainbow')
#cmap控制用于显示值的颜色映射
plt.colorbar()
plt.xlabel("Users(%d)" % nu)
plt.ylabel("Movies(%d)" % nm)
plt.show()
2.2 Collaborative filtering learning algorithm
现在，您将开始实现协作过滤学习算法。您将从实现成本函数开始（没有正则化）。
您将在cofiCostFunc中完成代码，以计算协作过滤的代价函数和梯度。请注意，函数的参数（即您尝试学习的值）是X和θ。为了使用fmincg等现成的最小化方法，建立了代价函数，将参数展开为单个向量参数。您以前在神经网络编程练习中使用过相同的向量展开方法。
path2 = "D:\编程\ex8_movieParams.mat"
mat = loadmat(path2)
print(mat.keys())
X = mat['X']
theta = mat['Theta']
nu = mat['num_users']
nm = mat['num_movies']
nf = mat['num_features']
print(X.shape,theta.shape,nu,nm,nf)
#dict_keys(['__header__', '__version__', '__globals__', 'X', 'Theta', 'num_users', 'num_movies', 'num_features'])
#(1682, 10) (943, 10) [[943]] [[1682]] [[10]]
# 由于nu,nm,nf都是矩阵形式，为了方便后续使用，将其转成int型
nu = int(nu)
nm = int(nm)
nf = int(nf)
2.2.1 Collaborative filtering cost function
请注意，只有当R(i,j)=1时，才应该累计用户j和电影i的成本。完成函数后应该期望看到22.22的输出。
我们强烈建议您使用矢量化实现来计算J，因为优化包fmincg稍后会多次调用它。像往常一样，首先编写一个非矢量化的实现（以确保您有正确的答案）并将其修改为矢量化的实现（检查矢量化步骤是否不会更改算法的输出）可能是最简单的。要实现一个矢量化的实现，下面的技巧可能会有所帮助：可以使用R矩阵将所选条目设置为0。例如，R.*M将在M和R之间进行元素相乘；由于R只有值为0或1的元素，因此只有当R中的相应值为0时，才会将M的元素设置为0。因此，sum（sum（R.*M））是M的所有元素的和，其中R中的相应元素等于1。
现在，您应该将正则化添加到成本函数J的原始计算中。完成之后运行正则化的代价函数，您应该看到大约31.34的代价。
# 序列化
def serialize(X,theta):
serial = np.append(X.flatten(), theta.flatten())
return serial
# 解序列化
def deserialize(serial,nm,nu,nf):
X = serial[:nm*nf].reshape(nm,nf)
theta = serial[nm*nf:].reshape(nu,nf)
return X,theta
# 代价函数
# def CostFunction(serial,Y,nm,nu,nf,lam,):
# X,theta = deserialize(serial,nm,nu,nf)
# first = 0
# for i in range(nm):
# for j in range(nu):
# if R[i,j] == 1:
# fir = np.power((X[i] @ theta[j].T - Y[i,j]),2)
# first = first + 0.5 * fir
# second = 0
# for j in range(nu):
# for k in range(nf):
# sec = (lam / 2) * np.power(theta[j,k],2)
# second = second + sec
# third = 0
# for i in range(nm):
# for k in range(nf):
# thi = (lam / 2) * np.power(X[i,k],2)
# third = third + thi
# costF = first + second + third
# return costF
# 上面的计算代价函数实在太复杂啦，因此，我们能要使用矢量化的方式去做！
def CostFunction(serial,nm,nu,nf,Y,R,lam):
X,theta = deserialize(serial,nm,nu,nf)
error = 0.5 * np.square((X @ theta.T - Y) * R).sum()
reg1 = (lam/2) * np.square(theta).sum()
reg2 = (lam/2) * np.square(X).sum()
costF = error + reg1 + reg2
return costF
users = 4
movies = 5
features = 3
X_sub = X[:movies,:features]
theta_sub = theta[:users,:features]
Y_sub = Y[:movies,:users]
R_sub = R[:movies,:users]
# 第一种代价函数输出
#costf = CostFunction(serialize(X_sub,theta_sub),Y_sub,movies,users,features,lam=0)
# 第二种代价函数输出
costf = CostFunction(serialize(X_sub,theta_sub),movies,users,features,Y_sub,R_sub,lam = 0)
print(costf)
#22.224603725685675
2.2.2 Collaborative filtering gradient
现在，您应该实现梯度（没有正则化）。具体来说，应该在cofiCostFunc中完成代码，以返回变量X grad和θgrad。请注意，X grad应该是与X大小相同的矩阵，类似地，θgrad是与θ大小相同的矩阵。成本函数的梯度如下所示：
请注意，该函数通过将两组变量展开为单个向量来返回它们的梯度。完成计算渐变的代码后，运行渐变检查（checkCostFunction）以数字方式检查渐变的实现。如果实现正确，您应该会发现分析渐变和数字渐变非常匹配。
这意味着您只需要添加 λx(i) 添加到前面描述的x grad(i,:)变量中，并将λθ(j）到前面描述的θ_grad(j,:)变量。在完成计算渐变的代码之后，将运行另一个渐变检查（checkCostFunction），以数字方式检查渐变的实现。
# 梯度下降
#计算X和theta的梯度，并将其进行序列化输出
def Gradient(serial,nm,nu,nf,Y,R,lam):
X,theta = deserialize(serial,nm,nu,nf)
# *R的目的是将r(i,j) = 0的过滤掉
X_grad = ((X @ theta.T - Y) * R) @ theta + lam * X
theta_grad = ((X @ theta.T - Y) * R).T @ X + lam * theta
ser = serialize(X_grad, theta_grad)
return ser
# ser = Gradient(serialize(X_sub,theta_sub),movies,users,features,Y_sub,R_sub,lam = 1)
# print(ser)
# 提醒：不要将内部函数和变量名同名，容易报错！
2.3 Learning movie recommendations
完成协作过滤成本函数和梯度的实现后，现在可以开始训练算法，为自己推荐电影。在下一部分中，您可以输入自己的电影首选项，以便稍后当算法运行时，您可以获得自己的电影推荐！我们已经根据自己的喜好填写了一些数值，但您应该根据自己的喜好进行更改。数据集中所有电影及其编号的列表可以在文件movie_idx.txt中找到。
2.3.1 Recommendations
将额外的评分添加到数据集后，脚本将继续训练协作过滤模型。这将学习参数X和θ。为了预测用户j对电影i的评价，需要计算(θ(j))T X(i)。下一部分计算所有电影和用户的收视率，并根据脚本前面输入的收视率显示它推荐的电影。注意，由于不同的随机初始化，您可能会获得不同的预测集。
# 添加一个新用户
my_ratings = np.zeros((nm,1))
my_ratings[9] = 5
my_ratings[12] = 4
my_ratings[56] = 3
my_ratings[78] = 5
my_ratings[278] = 4
my_ratings[568] = 5
my_ratings[980] = 3
my_ratings[1090] = 2
Y = np.c_[Y,my_ratings]
# R为二值的0或者1
R = np.c_[R,my_ratings != 0]
print(Y.shape,R.shape)
nm,nu = Y.shape
# 均值归一化
# 以下方法是将每一行看作整体计算的
# def normalizeRatings(Y,R):
# Y_norm = np.zeros(Y.shape)
# Y_means = np.zeros((Y.shape[0],1))
# for row in range(Y.shape[0]):
# mean = np.sum(Y[row,:]) / np.sum(R[row,:])
# Y_norm[row,:] = Y[row,:] - mean
# Y_means[row,:] = mean
# return Y_norm,Y_means
# 第二种方法，将矩阵看作一个整体进行计算
def normalizeRatings(Y,R):
#axis=1按行的方向相加，返回每个行的值；axis=0按列相加，返回每个列的值。
Y_mean = (np.sum(Y,axis=1) / np.sum(R,axis=1)).reshape(-1,1)
Y_norm = (Y - Y_mean) * R
# 这里要注意不要归一化未评分的数据
return Y_norm,Y_mean
Y_norm,Y_mean = normalizeRatings(Y,R)
# 参数初始化
X = np.random.random((nm,nf))
theta = np.random.random((nu,nf))
serial = serialize(X,theta)
lam = 5
# 模型训练
from scipy.optimize import minimize
res = minimize(fun=CostFunction,
x0 = serial,
args = (nm,nu,nf,Y_norm,R,lam),
method='TNC',
jac = Gradient,
options={'maxiter':100})
serialize_fit = res.x
fit_X,fit_theta = deserialize(serialize_fit,nm,nu,nf)
# 预测
Y_pred = fit_X @ fit_theta.T
y_pred = Y_pred[:,-1] + Y_mean.flatten()
# y预测值 = y预测归一化值 + 平均值
# 排序，使之从大到小排列
#np.argsort(-y_pred)或者np.argsort(y_pred)[::-1]
index = np.argsort(y_pred)[::-1]
# 输出推荐给我的排名前十的电影
print(index[:10])
movies = []
with open('D:\编程\movie_ids.txt','r',encoding = 'gbk') as f:
for line in f:
tokens = line.strip().split(' ')
movies.append(''.join(tokens[1:]))
print(len(movies))
for i in range(10):
# 预测的排名前十的电影索引，电影名字，以及评分
print(index[i],movies[index[i]],y_pred[index[i]])
完整代码如下：
# 推荐系统
#encoding=utf-8
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.io import loadmat
path = "D:\编程\ex8_movies.mat"
data = loadmat(path)
print(data.keys())
Y = data['Y']
R = data['R']
nm,nu = Y.shape
print(Y.shape,R.shape)
print(nm,nu)
# plt.figure(figsize=(6,6))
# plt.imshow(Y,cmap='rainbow')
# #cmap控制用于显示值的颜色映射
# plt.colorbar()
# plt.xlabel("Users(%d)" % nu)
# plt.ylabel("Movies(%d)" % nm)
# plt.show()
print("*********************************")
path2 = "D:\编程\ex8_movieParams.mat"
mat = loadmat(path2)
print(mat.keys())
X = mat['X']
theta = mat['Theta']
nu = mat['num_users']
nm = mat['num_movies']
nf = mat['num_features']
print(X.shape,theta.shape,nm,nu,nf)
# 由于nu,nm,nf都是矩阵形式，为了方便后续使用，将其转成int型
nu = int(nu)
nm = int(nm)
nf = int(nf)
# 序列化
def serialize(X,theta):
serial = np.append(X.flatten(), theta.flatten())
return serial
# 解序列化
def deserialize(serial,nm,nu,nf):
X = serial[:nm*nf].reshape(nm,nf)
theta = serial[nm*nf:].reshape(nu,nf)
return X,theta
# 代价函数
# def CostFunction(serial,Y,nm,nu,nf,lam,):
# X,theta = deserialize(serial,nm,nu,nf)
# first = 0
# for i in range(nm):
# for j in range(nu):
# if R[i,j] == 1:
# fir = np.power((X[i] @ theta[j].T - Y[i,j]),2)
# first = first + 0.5 * fir
# second = 0
# for j in range(nu):
# for k in range(nf):
# sec = (lam / 2) * np.power(theta[j,k],2)
# second = second + sec
# third = 0
# for i in range(nm):
# for k in range(nf):
# thi = (lam / 2) * np.power(X[i,k],2)
# third = third + thi
# costF = first + second + third
# return costF
# 上面的计算代价函数实在太复杂啦，因此，我们能要使用矢量化的方式去做！
def CostFunction(serial,nm,nu,nf,Y,R,lam):
X,theta = deserialize(serial,nm,nu,nf)
error = 0.5 * np.square((X @ theta.T - Y) * R).sum()
reg1 = (lam/2) * np.square(theta).sum()
reg2 = (lam/2) * np.square(X).sum()
costF = error + reg1 + reg2
return costF
users = 4
movies = 5
features = 3
X_sub = X[:movies,:features]
theta_sub = theta[:users,:features]
Y_sub = Y[:movies,:users]
R_sub = R[:movies,:users]
# 第一种代价函数输出
#costf = CostFunction(serialize(X_sub,theta_sub),Y_sub,movies,users,features,lam=0)
# 第二种代价函数输出
costf = CostFunction(serialize(X_sub,theta_sub),movies,users,features,Y_sub,R_sub,lam = 0)
print(costf)
# 梯度下降
#计算X和theta的梯度，并将其进行序列化输出
def Gradient(serial,nm,nu,nf,Y,R,lam):
X,theta = deserialize(serial,nm,nu,nf)
# *R的目的是将r(i,j) = 0的过滤掉
X_grad = ((X @ theta.T - Y) * R) @ theta + lam * X
theta_grad = ((X @ theta.T - Y) * R).T @ X + lam * theta
ser = serialize(X_grad, theta_grad)
return ser
# ser = Gradient(serialize(X_sub,theta_sub),movies,users,features,Y_sub,R_sub,lam = 1)
# print(ser)
# 添加一个新用户
my_ratings = np.zeros((nm,1))
my_ratings[9] = 5
my_ratings[12] = 4
my_ratings[56] = 3
my_ratings[78] = 5
my_ratings[278] = 4
my_ratings[568] = 5
my_ratings[980] = 3
my_ratings[1090] = 2
Y = np.c_[Y,my_ratings]
# R为二值的0或者1
R = np.c_[R,my_ratings != 0]
print(Y.shape,R.shape)
nm,nu = Y.shape
# 均值归一化
# 以下方法是将每一行看作整体计算的
# def normalizeRatings(Y,R):
# Y_norm = np.zeros(Y.shape)
# Y_means = np.zeros((Y.shape[0],1))
# for row in range(Y.shape[0]):
# mean = np.sum(Y[row,:]) / np.sum(R[row,:])
# Y_norm[row,:] = Y[row,:] - mean
# Y_means[row,:] = mean
# return Y_norm,Y_means
# 第二种方法，将矩阵看作一个整体进行计算
def normalizeRatings(Y,R):
Y_mean = (np.sum(Y,axis=1) / np.sum(R,axis=1)).reshape(-1,1)
Y_norm = (Y - Y_mean) * R
# 这里要注意不要归一化未评分的数据
return Y_norm,Y_mean
Y_norm,Y_mean = normalizeRatings(Y,R)
# 参数初始化
X = np.random.random((nm,nf))
theta = np.random.random((nu,nf))
serial = serialize(X,theta)
lam = 5
# 模型训练
from scipy.optimize import minimize
res = minimize(fun=CostFunction,
x0 = serial,
args = (nm,nu,nf,Y_norm,R,lam),
method='TNC',
jac = Gradient,
options={'maxiter':100})
serialize_fit = res.x
fit_X,fit_theta = deserialize(serialize_fit,nm,nu,nf)
# 预测
Y_pred = fit_X @ fit_theta.T
y_pred = Y_pred[:,-1] + Y_mean.flatten()
# y预测值 = y预测归一化值 + 平均值
# 排序，使之从大到小排列
#np.argsort(-y_pred)或者np.argsort(y_pred)[::-1]
index = np.argsort(y_pred)[::-1]
# 输出推荐给我的排名前十的电影
print(index[:10])
movies = []
with open('D:\编程\movie_ids.txt','r',encoding = 'gbk') as f:
for line in f:
tokens = line.strip().split(' ')
movies.append(''.join(tokens[1:]))
print(len(movies))
for i in range(10):
# 预测的排名前十的电影索引，电影名字，以及评分
print(index[i],movies[index[i]],y_pred[index[i]])
988
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
