机器学习练习八:用Python实现异常检测和协同过滤算法

981 阅读10分钟

在本练习中,我们将使用高斯模型实现异常检测算法,并将其应用于检测网络上的故障服务器。我们还将看到如何使用协作过滤构建推荐系统,并将其应用于电影推荐数据集。

1.异常检测模型

我们的第一个任务是使用高斯模型来检测数据集中未标记的示例是否应被视为异常。我们有一个简单的二维数据集开始,以帮助可视化该算法正在做什么。

1.1 展示原始数据

在本练习中,您将实现一个异常检测算法来检测服务器计算机中的异常行为。这些特性度量每个服务器的吞吐量(mb/s)和响应延迟(ms)。

# 展示原始数据
plt.figure(figsize=(15, 8))
plt.scatter(X[:, 0], X[:, 1])
plt.title('Source Data')
plt.xlabel('响应延迟(ms)')
plt.ylabel('吞吐量(mb/s)')
plt.show()

1.2 高斯分布

即正态分布,建立关于样本特征值的高斯分布模型,找到概率密度低于阙值的异常点。

概率密度公式:

第一步是找到适合模型的期望与方差。

代码如下:

def estimate_uv(X):
    """
    计算数据的均值与方差
    """
    # 以列为轴计算
    mu = X.mean(axis=0)
    sigma = X.var(axis=0)
    
    return mu, sigma

1.3 绘制3D图和等高线图

利用1.2计算出的均值与方差绘制3D图和等高线图。

1.3.1 绘制3D图

关于绘制3D图的技巧详见之前的文章机器学习练习一:用python实现线性回归第1.6节。

另外,在本次实验中使用了多元高斯分布,与普通的高斯分布相比,它能够体现自变量之间的相关性。详见:Scipy

# 绘制3D图
from mpl_toolkits.mplot3d import Axes3D
# 创建对象
fig = plt.figure(figsize=(10, 5))
ax = Axes3D(fig)
# 作为X轴和Y轴
X0, X1 = np.meshgrid(
        # 随机两组数
        np.linspace(5, 25, num=1000).reshape(-1, 1),
        np.linspace(5, 25, num=1000).reshape(-1, 1),
    )
# ravel()方法将高维数组降为一维数组,c_[]将两个数组以列的形式拼接起来,形成矩阵
X_grid_matrix = np.c_[X0.ravel(), X1.ravel()]

# 这一部分是计算高斯分布概率
multi_normal = stats.multivariate_normal(mu, sigma)
# 需要形状和X0与X1一致
Z = multi_normal.pdf(X_grid_matrix).reshape(X0.shape)

# 设置标签
ax.set_xlabel("响应延迟")
ax.set_ylabel("吞吐量")
ax.set_zlabel("概率密度")

# 绘制3D曲面图
ax.plot_surface(X0, X1, Z, rstride=20, cstride=20, cmap=plt.get_cmap('rainbow'))
plt.show()

结果如图:

1.3.2 绘制等高线图

等高线图与3D图的绘制原理相似,不再赘述。

# 作为X轴和Y轴
X0, X1 = np.meshgrid(
        # 随机两组数
        np.linspace(10, 18, num=300).reshape(-1, 1),
        np.linspace(10, 20, num=300).reshape(-1, 1),
    )

# ravel()方法将高维数组降为一维数组,c_[]将两个数组以列的形式拼接起来,形成矩阵
X_grid_matrix = np.c_[X0.ravel(), X1.ravel()]

# 这一部分是计算高斯分布概率
multi_normal = stats.multivariate_normal(mu, sigma)
# 需要形状和X0与X1一致
Z = multi_normal.pdf(X_grid_matrix).reshape(-1, 300)
# 创建对象
fig, ax = plt.subplots(figsize=(15, 8))

# 绘制等高线图
ax.contourf(X0, X1, Z, cmap='Reds')
# 展示原始数据
ax.scatter(X[:, 0], X[:, 1])
plt.title('Source Data')
plt.xlabel('响应延迟(ms)')
plt.ylabel('吞吐量(mb/s)')
plt.show()

结果如图:

可以看到异常点在概率密度为零的区域。

1.4 确定阙值

接下来,我们需要一个函数,找到给定概率密度值和真实标签的最佳阈值。 为了做到这一点,我们将为不同的ε值计算F1-score。 F1是真阳性(tp),假阳性(fp)和假阴性(fn)的数量的函数。

首先是获取训练集与验证集的概率密度,我们计算并保存给定上述的高斯模型参数的数据集中每个值的概率密度,还需要为验证集(使用相同的模型参数)执行此操作。我们将使用交叉验证的方法来确定将数据点分配为异常的最佳概率阈值。:

# 多元高斯分布
multi_normal = stats.multivariate_normal(mu, sigma)
# 训练集概率密度
p = multi_normal.pdf(X)
p.shape
# 获取验证集
Xval, yval = data['Xval'], data['yval']
Xval.shape, yval.shape
# 验证集概率密度
pval = multi_normal.pdf(Xval)
pval.shape

选取阙值:

注意yval是矩阵形式,需要展开为一维数组,否则逻辑与的结果将会出现错误。

def select_threshold(pval, yval):
    """
    选取阙值
    """
    # 用来存储最佳值
    best_epsilon = 0
    best_f1 = 0
    f1 = 0
    
    # 从验证集数值范围内抽取1000个预设ε值
    epsilon_list = np.linspace(pval.min(), pval.max(), 1000)
    
    # 循环对比最优的F1-score,其对应的ε值也是最优值
    for epsilon in epsilon_list:
        # 将验证集概率密度与ε值对比,返回预测值,异常点为正样本
        preds = pval < epsilon
        # 计算tp,fp,fn的值,注意yval是矩阵形式,需要展开为一维数组
        tp = np.sum(np.logical_and(preds == 1, yval.ravel() == 1)).astype(float)
        fp = np.sum(np.logical_and(preds == 1, yval.ravel() == 0)).astype(float)
        fn = np.sum(np.logical_and(preds == 0, yval.ravel() == 1)).astype(float)
        # 计算预测率与召回率
        precision = tp / (tp + fp)
        recall = tp / (tp + fn)
        # 计算f1-score
        f1 = (2 * precision * recall) / (precision + recall)
        # 保存最优的值
        if f1 > best_f1:
            best_f1 = f1
            best_epsilon = epsilon
    
    return best_epsilon, best_f1

最佳ε值:

1.5 展示被标记为异常的点

注意outliers是一个元组

# 选取被标记为异常的点
outliers = np.where(p < epsilon)
# 显示异常值
plt.figure(figsize=(15, 8))
plt.scatter(X[:, 0], X[:, 1])
plt.scatter(X[outliers, 0][0], X[outliers, 1][0], color='r')
plt.title('Source Data')
plt.xlabel('响应延迟(ms)')
plt.ylabel('吞吐量(mb/s)')
plt.show()

结果如图:

可以看到所有的异常点都被标记为红色了。

2.基于协同过滤算法的推荐系统

推荐引擎使用基于项目和用户的相似性度量来检查用户的历史偏好,以便为用户可能感兴趣的新“事物”提供建议。在本部分练习中,您将实现协同过滤学习算法并将其应用于电影评级数据集。评级由1到5的等级组成。数据集有n_u = 943个用户,n_m = 1682个电影。

2.1 读取数据

矩阵Y(1682×943)存储评级Y(i;j) = [1,5],矩阵R为二值指标矩阵,其中R(i;j)(如果用户j对电影i进行评级,则j = 1, 否则R(i;j) = 0。协同过滤的目的是预测用户尚未打分的电影的评分,即含有R(i;j) = 0。这将允许我们向用户推荐预测收视率最高的电影。

接下来,我们将实施协同过滤的代价函数。 直觉上,“代价”是指一组电影评级预测偏离真实预测的程度。它基于文本中称为X和Theta的两组参数矩阵,其中参数θ是用户的偏好权重,参数X是电影的类型特征。参数需要展开为一维数组,以便稍后可以使用SciPy的优化包。

展开元素的代码:

def serialize(a, b):
    """
    展开元素
    """
    return np.concatenate((np.ravel(a), np.ravel(b)))

2.2 协同过滤的损失函数和梯度下降

协同过滤算法的本质是将关于用户喜好的参数θ和关于电影类型的参数X当做特征值,将评分值当做目标值,预测的评分值与真实值之差当做损失计算梯度。梯度包括两部分:即关于参数Θ的梯度和关于参数X的梯度。

函数公式如图:

代码实现:

显然在一个函数中先计算了损失,然后利用损失计算了梯度。

# 在Cost函数中计算
def cost(params, Y, R, num_features, lambd):
    """
    正则化的代价函数
    同时计算正则化梯度
    """
    # 将Y和R转换为矩阵类型
    Y = np.matrix(Y)  # (1682, 943)
    R = np.matrix(R)  # (1682, 943)
    # 获取电影样本数和用户数
    num_movies = Y.shape[0]
    num_users = Y.shape[1]
    
    # 将展开后的X和Theta恢复成原来的形状,并转换为矩阵类型
    X = np.matrix(np.reshape(params[:num_movies * num_features], (num_movies, num_features)))  # (1682, 10)
    Theta = np.matrix(np.reshape(params[num_movies * num_features:], (num_users, num_features)))  # (943, 10)
    
    # 计算损失值
    error = np.multiply((X * Theta.T) - Y, R)  # (1682, 943)
    J = (1. / 2) * np.sum(np.power(error, 2))
    J_regu = (lambd / 2) * (np.sum(np.power(Theta, 2)) + np.sum(np.power(X, 2)))
    J += J_regu
    
    # 计算梯度
    X_grad = error * Theta + lambd * X
    Theta_grad = error.T * X + lambd * Theta
    
    # 展开
    grad = serialize(X_grad, Theta_grad)
    return J, grad

2.3 自定义一个新用户

在我们训练模型之前,我们有一个最后步骤,我们的任务是创建自己的电影评分(即在数据集Y和R中加入一个新用户),以便我们可以使用该模型来生成个性化的推荐。提供一个连接电影索引到其标题的文件,接着我们将电影名加载到字典中。

2.3.1 创建一个存储电影名的字典

目的是方便可视化数据。创建的字典如图:

代码如下:

# 用来存储电影名的字典
movie_idx = {}
# 打开文件
f = open('data/movie_ids.txt',encoding= 'gbk')
# 每次取一行(以换行符结尾)
for line in f:
    # 将字符串以空格分隔为列表
    line_list = line.split(' ')
    # 将最后一个字符去掉换行符
    line_list[-1] = line_list[-1][:-1]
    # 去掉序号,存储到字典中
    movie_idx[int(line_list[0]) - 1] = ' '.join(line_list[1:])

2.3.2 构造一个新用户

  • 对电影评分
# 存储新用户评分的数组
ratings = np.zeros((1682, 1))

ratings[0] = 4
ratings[6] = 3
ratings[11] = 5
ratings[53] = 4
ratings[63] = 5
ratings[65] = 3
ratings[68] = 5
ratings[97] = 2
ratings[182] = 4
ratings[225] = 5
ratings[354] = 5

print('Rated {0} with {1} stars.'.format(movie_idx[0], str(int(ratings[0]))))
print('Rated {0} with {1} stars.'.format(movie_idx[6], str(int(ratings[6]))))
print('Rated {0} with {1} stars.'.format(movie_idx[11], str(int(ratings[11]))))
print('Rated {0} with {1} stars.'.format(movie_idx[53], str(int(ratings[53]))))
print('Rated {0} with {1} stars.'.format(movie_idx[63], str(int(ratings[63]))))
print('Rated {0} with {1} stars.'.format(movie_idx[65], str(int(ratings[65]))))
print('Rated {0} with {1} stars.'.format(movie_idx[68], str(int(ratings[68]))))
print('Rated {0} with {1} stars.'.format(movie_idx[97], str(int(ratings[97]))))
print('Rated {0} with {1} stars.'.format(movie_idx[182], str(int(ratings[182]))))
print('Rated {0} with {1} stars.'.format(movie_idx[225], str(int(ratings[225]))))
print('Rated {0} with {1} stars.'.format(movie_idx[354], str(int(ratings[354]))))

输出如图:

  • 将新用户加入到数据集中
R = data['R']
Y = data['Y']
# 将新用户的评级分数加入数据中
Y = np.append(Y, ratings, axis=1)
R = np.append(R, ratings != 0, axis=1)

2.4 初始化变量和均值归一化

为了进行模型训练,需要初始化一些变量并随机初始化参数X和Theta。

对于无法判断喜好的用户(即未对任何电影评分的用户),为了避免预测的评分均为0,需要对评分数据Y进行均值归一化。

2.4.1 初始化变量

# 初始化变量
movie_num = Y.shape[0]  # 1682
user_num = Y.shape[1]  # 944
features = 10
lambd = 10.
# 随机初始化X,Theta
X = np.random.random(size=(movie_num, features))
Theta = np.random.random(size=(user_num, features))
params = serialize(X, Theta)

2.4.2 均值归一化

# 初始化每部电影评分的均值
Ymean = np.zeros((movie_num, 1))
# 初始化均值归一化后的评分数据集
Ynorm = np.zeros((movie_num, user_num))
# 每次对一行样本操作
for i in range(movie_num):
    # 提取有评分的列
    idx = np.where(R[i,:] == 1)[0]
    # 求这一行的均值
    Ymean[i] = Y[i, idx].mean()
    # 这一行的评分减去均值
    Ynorm[i, idx] = Y[i, idx] - Ymean[i]

Ynorm为归一化之后的评分值。

2.5 模型训练与预测

2.5.1 模型训练

使用Scipy中的包。

# 模型训练
from scipy.optimize import minimize
fmin = minimize(fun=cost, x0=params, args=(Ynorm, R, features, lambd), 
                method='CG', jac=True, options={'maxiter': 100})
fmin

结果如下:

2.5.2 预测

我们要使用的是模型训练出的参数X和Theta,而预测值y_predict = X * θ.T。

  • 获取预测结果
# 将X和Theta恢复成原来的形状
X = np.matrix(np.reshape(fmin.x[:movie_num * features], (movie_num, features)))  # (1682, 10)
Theta = np.matrix(np.reshape(fmin.x[movie_num * features:], (user_num, features)))  # (943, 10)
# 根据训练好的用户喜好Theta和电影类型分析X,可以得到预测结果
predictions = X * Theta.T  # (1682, 10) * (10, 944)
  • 查看自定义用户的预测结果
# 取自定义的列
my_preds = predictions[:, -1] + Ymean
# 按评分从高到低获取索引
idx = np.argsort(my_preds, axis=0)[::-1]
# 展示预测评分最高的前十个电影
for i in idx[:10, :]:
    print('Predicted {0} with {1} stars.'.format(movie_idx[i[0, 0]], str(int(my_preds[i]))))

结果如图: