机器学习 | 吴恩达机器学习第九周编程作业(Python版）

最新推荐文章于 2022-08-13 14:22:06 发布

CoreJT

最新推荐文章于 2022-08-13 14:22:06 发布

阅读量1.1k

点赞数 4

分类专栏：林轩田机器学习文章标签：吴恩达机器学习实验异常检测高斯分布电影推荐系统协同过滤

本文链接：https://blog.csdn.net/sdu_hao/article/details/84338471

版权

林轩田机器学习专栏收录该内容

44 篇文章 47 订阅

订阅专栏

实验指导书下载密码:bgi7

本篇博客主要讲解，吴恩达机器学习第九周的编程作业，主要包含异常检测实验和电影推荐系统实验两部分。原始实验使用Matlab实现，本篇博客提供Python版本。

1.实验包含的文件

文件名	含义
ex8.py	异常检测实验主程序
ex8_cofi.py	推荐系统实验主程序
ex8data1.mat	异常检测第一个示例数据集
ex8data2.mat	异常检测第二个示例数据集
ex8_movies.mat	电影评分数据集
ex8_movieParams.mat	调试参数
multivariateGaussian.py	多元高斯分布程序
visualizeFit.py	可视化训练集及其概率分布
checkCostFunction.py	协同过滤梯度检测
computeNumericalGradient.py	计算每个参数的近似梯度
loadMovieList.py	加载电影列表
movie_ids.txt	电影id列表
normalizeRatings.py	协同过滤均值规范化
estimateGaussian.py	估计高斯分布参数
selectThreshold.py	选择异常检测的阈值
cofiCostFunc.py	实现协同过滤的代价函数

完成红色部分程序的关键代码。

2.异常检测实验

打开异常检测实验主程序ex8.py

'''第1部分 加载示例数据集'''

#先通过一个小数据集进行异常检测  便于可视化

# 数据集包含两个特征 
# 一些机器的等待时间和吞吐量  实验目的找出其中可能有异常的机器


print('Visualizing example dataset for outlier detection.')


data = scio.loadmat('ex8data1.mat')
X = data['X']#训练集样本特征矩阵
Xval = data['Xval']  #验证集样本特征矩阵
yval = data['yval'].flatten() #验证集样本标签 异常/正常 

# 可视化样例训练集
plt.figure()
plt.scatter(X[:, 0], X[:, 1], c='b', marker='x', s=15, linewidth=1)
plt.axis([0, 30, 0, 30])
plt.xlabel('Latency (ms)')  #x1等待时间
plt.ylabel('Throughput (mb/s') #x2吞吐量

样例训练集可视化效果

估计训练集的概率分布

'''第2部分 估计训练集的分布'''
# 假设数据集的各个特征服从高斯分布

print('Visualizing Gaussian fit.')

# 参数估计  
mu, sigma2 = eg.estimate_gaussian(X)

# 计算训练集的概率分布
p = mvg.multivariate_gaussian(X, mu, sigma2)
#可视化训练集的概率分布  画出等高线图
vf.visualize_fit(X, mu, sigma2)
plt.xlabel('Latency (ms)')
plt.ylabel('Throughput (mb/s')

编写参数估计程序estimateGaussian.py

def estimate_gaussian(X):
   
    m, n = X.shape #样本数m 特征数n
    

    mu = np.zeros(n)
    sigma2 = np.zeros(n)

  
    mu=np.mean(X,axis=0) #计算每个特征的均值
    sigma2=np.var(X,axis=0) #对于基于单元高斯分布的异常检测模型 直接计算每个特征的方差 返回一维数组
    #sigma2=(1/m)*(X.T.dot(X)) #对于基于多元高斯分布的异常检测模型 计算协方差矩阵
    return mu, sigma2

查看计算概率分布的程序multivariateGaussian.py

注意当使用基于单元高斯分布的异常检测模型时，因为当协方差矩阵是对角矩阵时，对角线上的元素是每个特征的方差，此时基于单元高斯分布的模型和基于多元的高斯分布模型是等价的，所以此时计算出每个特征的方差，把它转换为对角矩阵作为协方差矩阵，对原始特征矩阵进行均值标准化，再代入多元高斯分布公式计算训练集概率分布即可，而不是使用多个单元高斯分布函数累乘的公式。

如果想使用基于多元高斯分布的异常检测模型，此时协方差矩阵需要基于训练集X进行计算： $\Sigma = 1/m X^{T}X$ ,再对原始特征矩阵进行均值规范化，再代入多元高斯分布公式计算训练集概率分布即可。


def multivariate_gaussian(X, mu, sigma2):
    k = mu.size  #特征数
    
    #如果是基于单元高斯分布的模型  将其sigma2转换为对角矩阵 作为协方差矩阵 代入多元高斯分布公式
    #此时单元模型和多元模型是等价的
    #如果是基于多元高斯分布的模型 直接将计算的协方差矩阵sigma2代入多元高斯分布公式
    if sigma2.ndim == 1 or (sigma2.ndim == 2 and (sigma2.shape[1] == 1 or sigma2.shape[0] == 1)):
        sigma2 = np.diag(sigma2)
    
    x = X - mu #对原始特征矩阵进行均值规范化
    p = (2 * np.pi) ** (-k / 2) * np.linalg.det(sigma2) ** (-0.5) * np.exp(-0.5*np.sum(np.dot(x, np.linalg.pinv(sigma2)) * x, axis=1))

    return p

查看可视化训练集及其概率分布程序 visualizeFit.py

#可视化训练集的概率分布  画出等高线图
def visualize_fit(X, mu, sigma2):
    grid = np.arange(0, 35.5, 0.5) #生成网格点
    x1, x2 = np.meshgrid(grid, grid)

    Z = mvg.multivariate_gaussian(np.c_[x1.flatten('F'), x2.flatten('F')], mu, sigma2) #得到每个网格点的概率
    Z = Z.reshape(x1.shape, order='F') 

    plt.figure() #画出训练集样本
    plt.scatter(X[:, 0], X[:, 1], marker='x', c='b', s=15, linewidth=1)

    if np.sum(np.isinf(X)) == 0:  #画出训练集概率分布的等高线图
        lvls = 10 ** np.arange(-20, 0, 3).astype(np.float)
        plt.contour(x1, x2, Z, levels=lvls, colors='r', linewidths=0.7)

可视化效果

基于验证集选择一个最好的阈值参数

'''第3部分 基于验证集 得到一个最好的概率分布阈值'''
pval = mvg.multivariate_gaussian(Xval, mu, sigma2) #根据训练集的概率分布 得到验证集样本的概率

epsilon, f1 = st.select_threshold(yval, pval)   #选择合适的概率阈值
print('Best epsilon found using cross-validation: {:0.4e}'.format(epsilon))
print('Best F1 on Cross Validation Set: {:0.6f}'.format(f1))
print('(you should see a value epsilon of about 8.99e-05 and F1 of about 0.875)')

# 标出训练集中的异常值
outliers = np.where(p < epsilon)
plt.scatter(X[outliers, 0], X[outliers, 1], marker='o', facecolors='none', edgecolors='r')

编写阈值选择程序selectThreshold.py

def select_threshold(yval, pval):
    f1 = 0 #f1-score

    best_eps = 0  #最好的阈值参数
    best_f1 = 0#最好的阈值参数对应的f1-score
    y=np.zeros(yval.size) #存放预测值
    for epsilon in np.linspace(np.min(pval), np.max(pval), num=1001): #尝试不同的阈值
        
        #y=1 异常 ；y=0正常
        y[pval<epsilon]=1  #预测为异常的样本
        tp=np.sum([yval[x] for x in range(len(y)) if y[x]])  #true positive
        
        precision=tp/np.sum(y) #查准率
        recall=tp/np.sum(yval)   #召回率
        
        f1=2*precision*recall/(precision+recall)  #f1-score

        if f1 > best_f1: #得到对应最高f1-score的阈值
            best_f1 = f1
            best_eps = epsilon
        y=np.zeros(yval.size)

    return best_eps, best_f1

基于最好的阈值，标出训练集中的异常样本

验证程序正确性：

在大数据集上进行异常检测

异常检测的训练过程是无监督的过程，得到训练集的概率分布，训练集一般全为正常样本，不过有异常样本存在，也没关系。将验证集样本代入训练得到的概率分布和验证集真实label进行比较，来进行模型选择，比如选择阈值参数或是否应该添加新特征等。

'''第4部分 基于大数据集 进行异常检测（特征数很多）'''
data = scio.loadmat('ex8data2.mat')
X = data['X'] #训练集样本特征矩阵
Xval = data['Xval'] #验证集样本特征矩阵
yval = data['yval'].flatten() #验证集样本标签 1异常 0正常

#参数估计
mu, sigma2 = eg.estimate_gaussian(X)

# 计算训练集的概率分布
p = mvg.multivariate_gaussian(X, mu, sigma2)

# 得到验证集每个样本的概率
pval = mvg.multivariate_gaussian(Xval, mu, sigma2)

# 选择一个最好的阈值
epsilon, f1 = st.select_threshold(yval, pval)

#验证程序正确性
print('Best epsilon found using cross-validation: {:0.4e}'.format(epsilon))
print('Best F1 on Cross Validation Set: {:0.6f}'.format(f1))
print('# Outliers found: {}'.format(np.sum(np.less(p, epsilon)))) #训练集上的异常样本数量
print('(you should see a value epsilon of about 1.38e-18, F1 of about 0.615, and 117 outliers)')

验证程序正确性

3.异常检测实验完整代码

下载链接下载密码:3u81

4.推荐系统实验

打开推荐系统实验主程序ex8_cofi.py

'''第1部分 加载电影评分数据集'''
print('Loading movie ratings dataset.')


data = scio.loadmat('ex8_movies.mat')
Y = data['Y']  #1682*943  1682部电影 943个用户 评分（1-5）   评分矩阵
R = data['R']  #1682*943  1682部电影 943个用户 当且仅当R(i,j)=1时 Y(i,j)有评分



# 计算第一部电影的平均评分
print('Average ratings for movie 0(Toy Story): {:0.6f}/5'.format(np.mean(Y[0, np.where(R[0] == 1)])))

# 可视化评分矩阵
plt.figure()
plt.imshow(Y)
plt.colorbar()
plt.xlabel('Users')
plt.ylabel('Movies')

可视化评分矩阵

计算协同过滤的代价

'''第2部分 计算协同过滤的代价'''

data = scio.loadmat('ex8_movieParams.mat') #加载预训练好的协同过滤的参数
X = data['X']    #电影特征向量构成的矩阵
theta = data['Theta']  #用户喜好向量构成的矩阵
num_users = data['num_users']  #用户数量
num_movies = data['num_movies']  #电影数量
num_features = data['num_features']  #特征数量

#缩减数据集 使运行速度更快
num_users = 4
num_movies = 5
num_features = 3
X = X[0:num_movies, 0:num_features] #选择前5部电影的特征向量 其中每个特征向量只取前3个特征
theta = theta[0:num_users, 0:num_features]#选择前4位用户的喜好向量 其中每个向量只取前3个特征
Y = Y[0:num_movies, 0:num_users] #缩减后的评分矩阵
R = R[0:num_movies, 0:num_users] #缩减后的R矩阵

# 根据训练好的参数向量  计算此时的代价
cost, grad = ccf.cofi_cost_function(np.concatenate((X.flatten(), theta.flatten())), Y, R, num_users, num_movies, num_features, 0)

#验证程序正确性
print('Cost at loaded parameters: {:0.2f}\n(this value should be about 22.22)'.format(cost))

编写计算代价的程序cofiCostfunc.py

不带正则化(lmd=0)：

def cofi_cost_function(params, Y, R, num_users, num_movies, num_features, lmd):
    #将一维向量 转为矩阵形式
    X = params[0:num_movies * num_features].reshape((num_movies, num_features)) 
    theta = params[num_movies * num_features:].reshape((num_users, num_features))

    cost = 0
    X_grad = np.zeros(X.shape)  #存放电影特征向量梯度的矩阵
    theta_grad = np.zeros(theta.shape) #存放用户喜好向量梯度的矩阵


    term1=np.sum((X.dot(theta.T)-Y)*(X.dot(theta.T)-Y))/2
    term1=term1*R    #只考虑有评分的数据
    term2=(lmd/2)*((X.dot(X.T)).diagonal().sum()+(theta.dot(theta.T)).diagonal().sum())#正则化惩罚项
    cost=term1+term2
    
    X_grad=(X.dot(theta.T)-Y)*R.dot(theta)+lmd*X
    theta_grad=((X.dot(theta.T)-Y)*R).T.dot(X)+lmd*theta
    grad = np.concatenate((X_grad.flatten(), theta_grad.flatten())) #把所有参数的梯度放在一个一维向量中

    return cost, grad

验证程序正确性（lmd=0）：

进行梯度检查

'''第3部分 进行梯度检查'''

print('Checking gradients (without regularization) ...')

#构建小数据集  对每一个参数进行梯度检查 （可以理解为 弦的斜率是否近似于切线的斜率）
cf.check_cost_function(0)

查看梯度检查程序checkCostFunction.py

def check_cost_function(lmd):

    #构建一个小数据集 进行梯度检查
    x_t = np.random.rand(4, 3)  #构建电影的特征向量矩阵 4部电影 每部电影3个特征
    theta_t = np.random.rand(5, 3)#构建用户的喜好向量矩阵 5位用户 对应的3个特征

    
    Y = np.dot(x_t, theta_t.T)  # 4x5 得到评分矩阵
    Y[np.random.rand(Y.shape[0], Y.shape[1]) > 0.5] = 0 #去掉一部分评分值（将一部分随机生成的评分值置0）
    R = np.zeros(Y.shape)  #R矩阵与Y同大小
    R[Y != 0] = 1  #Y矩阵有评分的位置 R[i][j]=1


    x = np.random.randn(x_t.shape[0], x_t.shape[1]) #存放电影特征向量梯度的矩阵
    theta = np.random.randn(theta_t.shape[0], theta_t.shape[1]) #存放用户喜好向量梯度的矩阵
    num_users = Y.shape[1]  #5
    num_movies = Y.shape[0]  #4
    num_features = theta_t.shape[1] #3

    def cost_func(p):
        return ccf.cofi_cost_function(p, Y, R, num_users, num_movies, num_features, lmd)

    #返回每个参数的近似梯度  可以理解为2维空间中参数加减一个很小的数得到的弦的斜率（梯度）
    numgrad = cng.compute_numerial_gradient(cost_func, np.concatenate((x.flatten(), theta.flatten()))) #把两个梯度矩阵转型为一个一维向量 计算每一个分量的近似梯度
     
    #返回每个参数的真实梯度  可以理解为2维空间中参数的切线斜率（梯度）
    cost, grad = ccf.cofi_cost_function(np.concatenate((x.flatten(), theta.flatten())), Y, R, num_users, num_movies, num_features, lmd)
    #对2者进行比较
    print(np.c_[numgrad, grad])
    print('The above two columns you get should be very similar.\n'
          '(Left-Your Numerical Gradient, Right-Analytical Gradient')
    #如果该数值非常小 即2者非常接近 说明没有问题
    diff = np.linalg.norm(numgrad - grad) / np.linalg.norm(numgrad + grad)
    print('If you backpropagation implementation is correct, then\n'
          'the relative difference will be small (less than 1e-9).\n'
          'Relative Difference: {:0.3e}'.format(diff))

查看计算参数近似梯度的程序computeNumericalGradient.py

def compute_numerial_gradient(cost_func, theta):
    #theta存放所有参数 包括所有特征向量和喜好向量的每一个分量
    numgrad = np.zeros(theta.size)
    perturb = np.zeros(theta.size)

    e = 1e-4
    #用弦的斜率（梯度）与切线的斜率（梯度）进行比较  来进行梯度检查 如果差不多 就说明没问题
    #对每个参数都进行梯度检查 
    for p in range(theta.size): 
        perturb[p] = e
        #返回弦两个端点的代价
        loss1, grad1 = cost_func(theta - perturb) 
        loss2, grad2 = cost_func(theta + perturb)

        numgrad[p] = (loss2 - loss1) / (2 * e) #弦的斜率（梯度）
        perturb[p] = 0

    return numgrad

编写cofiCostFunc.py 计算协同过滤每个参数的梯度

def cofi_cost_function(params, Y, R, num_users, num_movies, num_features, lmd):
    #将一维向量 转为矩阵形式
    X = params[0:num_movies * num_features].reshape((num_movies, num_features)) 
    theta = params[num_movies * num_features:].reshape((num_users, num_features))

    cost = 0
    X_grad = np.zeros(X.shape)  #存放电影特征向量梯度的矩阵
    theta_grad = np.zeros(theta.shape) #存放用户喜好向量梯度的矩阵


    term1=np.sum((X.dot(theta.T)-Y)*(X.dot(theta.T)-Y))/2
    term1=term1*R    #只考虑有评分的数据
    term2=(lmd/2)*((X.dot(X.T)).diagonal().sum()+(theta.dot(theta.T)).diagonal().sum())#正则化惩罚项
    cost=term1+term2
    
    X_grad=(X.dot(theta.T)-Y)*R.dot(theta)+lmd*X
    theta_grad=((X.dot(theta.T)-Y)*R).T.dot(X)+lmd*theta
    grad = np.concatenate((X_grad.flatten(), theta_grad.flatten())) #把所有参数的梯度放在一个一维向量中

    return cost, grad

在构建的小数据集上验证计算梯度的程序和梯度检查程序的正确性：

计算协同过滤的代价（带正则化）


'''第4部分 计算协同过滤的代价 带正则化'''
#lmd=1.5
cost, _ = ccf.cofi_cost_function(np.concatenate((X.flatten(), theta.flatten())), Y, R, num_users, num_movies, num_features, 1.5)

#验证程序正确性
print('Cost at loaded parameters (lambda = 1.5): {:0.2f}\n'
      '(this value should be about 31.34)'.format(cost))

验证程序正确性：

进行梯度检查（带正则化）

'''第5部分 进行梯度检查 （带正则化）'''
#计算带正则化的近似梯度和梯度 并进行比较
print('Checking Gradients (with regularization) ...')

#lmd=1.5
cf.check_cost_function(1.5)

验证程序正确性：

增加一个新用户，并对所有电影进行评分

'''第6部分 增加一个新用户 并对所有的1682部电影进行评分'''
movie_list = lm.load_movie_list() #读取电影名放在列表里

# 初始新用户对每部电影的评分为0 表示新用户都没看过
my_ratings = np.zeros(len(movie_list))

#新用户对其中看过的部分电影进行评分（1-5）
#比如新用户看过第一部电影 并给他的评分为4 其他类似
my_ratings[0] = 4

my_ratings[97] = 2

my_ratings[6] = 3
my_ratings[11] = 5
my_ratings[53] = 4
my_ratings[63] = 5
my_ratings[65] = 3
my_ratings[68] = 5
my_ratings[182] = 4
my_ratings[225] = 5
my_ratings[354] = 5

print('New user ratings:\n')
for i in range(my_ratings.size):
    if my_ratings[i] > 0:
        print('Rated {} for {}'.format(my_ratings[i], movie_list[i]))

在原始评分数据集上增加新用户的评分重新训练协同过滤算法

'''第7部分 在原始电影评分数据集上，增加新用户的评分 即增加一列数据 重新训练协同过滤算法'''
print('Training collaborative filtering ...\n'
      '(this may take 1 ~ 2 minutes)')


# 加载原始数据
data = scio.loadmat('ex8_movies.mat')
Y = data['Y']
R = data['R']

# Y is a 1682x943 matrix, containing ratings (1-5) of 1682 movies by
# 943 users
#
# R is a 1682x943 matrix, where R[i,j] = 1 if and only if user j gave a
# rating to movie i

# 加上新用户的数据 构建新的电影评分数据集
Y = np.c_[my_ratings, Y]
R = np.c_[(my_ratings != 0), R]

# 规范化评分矩阵 计算每部电影的平均评分 原始评分矩阵每一行减去各自的平均评分 返回新的评分矩阵
Ynorm, Ymean = nr.normalize_ratings(Y, R)

num_users = Y.shape[1] #用户数
num_movies = Y.shape[0] #电影数
num_features = 10  #特征向量取10个特征

# 初始化所有参数为很小的随机值
X = np.random.randn(num_movies, num_features)
theta = np.random.randn(num_users, num_features)

initial_params = np.concatenate([X.flatten(), theta.flatten()]) #把所有参数放在一个一维向量中

lmd = 10 #正则化系数


def cost_func(p):  #返回协同过滤的代价
    return ccf.cofi_cost_function(p, Ynorm, R, num_users, num_movies, num_features, lmd)[0]


def grad_func(p):  #返回协同过滤的梯度
    return ccf.cofi_cost_function(p, Ynorm, R, num_users, num_movies, num_features, lmd)[1]

#调用高级优化算法 进行训练 得到最优的参数
theta, *unused = opt.fmin_cg(cost_func, fprime=grad_func, x0=initial_params, maxiter=100, disp=False, full_output=True)

# 把训练好的参数 再重新转型为矩阵形式
X = theta[0:num_movies * num_features].reshape((num_movies, num_features))
theta = theta[num_movies * num_features:].reshape((num_users, num_features))

#打印用户的喜好向量矩阵
print('Recommender system learning completed')
print(theta)

查看均值规范化评分矩阵的程序 normalizeRatings.py

def normalize_ratings(Y, R):
    m, n = Y.shape  
    Ymean = np.zeros(m)
    Ynorm = np.zeros(Y.shape)
    for i in range(m):
        idx = np.where(R[i] == 1)
        Ymean[i] = np.mean(Y[i, idx])  #计算每一部电影的平均评分 仅对有评分的数据计算 
        Ynorm[i, idx] = Y[i, idx] - Ymean[i]  #得到新的评分矩阵   对有评分的数据减去对应电影的平均评分

    return Ynorm, Ymean

打印训练好的944位用户喜好向量构成的矩阵（944*10）：

根据训练好的参数，预测评分矩阵

'''第8部分 根据训练好的参数 预测评分矩阵的缺失值'''

p = np.dot(X, theta.T) #X训练好的每部电影的特征向量构成的矩阵 theta训练好的每位用户的喜好向量构成的矩阵  得到预测的评分矩阵
my_predictions = p[:, 0] + Ymean #得到预测的第一位新添加用户对每部电影的评价 别忘了加上减去的每部电影的平均评分

indices = np.argsort(my_predictions)[::-1] #得到预测的新用户对每部电影的评价的排序 从大到小  返回电影的索引
print('\nTop recommendations for you:')  #为新用户推荐的前10位的电影
for i in range(10):
    j = indices[i]
    print('Predicting rating {:0.1f} for movie {}'.format(my_predictions[j], movie_list[j])) #打印预测评分+电影名称

print('\nOriginal ratings provided:') #新用户对部分电影的原始评分+对应的电影名
for i in range(my_ratings.size):
    if my_ratings[i] > 0:
        print('Rated {} for {}'.format(my_ratings[i], movie_list[i]))

5.推荐系统实验完整代码

下载链接下载密码:uhae

CoreJT

关注

4
点赞
踩
13

收藏

觉得还不错? 一键收藏
0
评论
机器学习 | 吴恩达机器学习第九周编程作业(Python版）

实验指导书下载密码:bgi7本篇博客主要讲解，吴恩达机器学习第九周的编程作业，主要包含异常检测实验和电影推荐系统实验两部分。原始实验使用Matlab实现，本篇博客提供Python版本。目录1.实验包含的文件2.异常检测实验3.异常检测实验完整代码4.推荐系统实验5.推荐系统实验完整代码1.实验包含的文件文件名含义 ex8.py ...
复制链接

扫一扫