机器学习练习（八）—— 异常检测与推荐

机器学习练习（八）—— 异常检测与推荐

作者：John Wittenauer

翻译：GreatX

源：Machine Learning Exercises In Python, Part 8

这篇文章是一系列 Andrew Ng 在 Coursera 上的机器学习课程的练习的一部分。这篇文章的原始代码，练习文本，数据文件可从这里获得。

Part 1 简单线性回归（Simple Linear Regression）

Part 2 多元线性回归（Multivariate Linear Regression）

Part 3 逻辑回归（Logistic Regression）

Part 4 多元逻辑回归（Multivariate Logistic Regression）

Part 5 神经网络（Neural Networks）

Part 6 支持向量机（Support Vector Machines）

Part 7 K-均值聚类与主成分分析（K-Means Clustering & PCA）

Part 8 异常检测与推荐（Anomaly Detection & Recommendation）

这是本系列的最后一篇文章了，这真是一次有趣的旅程。 Andrew 的课真的很好，而将其用 python 重写也是一次有趣的体验。 在最后一部分中，我们将介绍本课程的最后两个主题 —— 异常检测（anomaly detection）和推荐系统（recommendation system）。 我们将用高斯模型实现一个异常检测算法，并将其应用于检测在网络上的故障服务器。 我们还将了解如何使用协同过滤（collaborative filtering）构建推荐系统，并将其应用到电影推荐数据集。 一如往常，阅读练习文本将有助于我们理解（上传在此）。

Anomaly Detection

我们的第一个任务是使用高斯模型来检测数据集中一个未标记（unlabeled）的样本是否应被视为异常。 我们从一个简单的 2 维数据集开始，因此我们可以很容易地可视化这算法的工作原理。 让我们 导入数据并画出散点图。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sb
from scipy.io import loadmat
%matplotlib inline

data = loadmat('data/ex8data1.mat')
X = data['X']
X.shape

(307L, 2L)

fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(X[:,0], X[:,1])

看起来该聚类的中心十分紧密，只有几个值远离聚类中心。 在这个简单的例子中，这些可以被认为是异常。 为了弄清楚这些，我们的任务是，估计数据中每个特征的高斯分布。 你可能还记得，为了定义概率分布，我们需要两个东西——均值和方差。 为了实现这一点，我们将创建一个简单的函数来计算数据集中每个特征的均值和方差。

def estimate_gaussian(X):
mu = X.mean(axis=0)
sigma = X.var(axis=0)

return mu, sigma

mu, sigma = estimate_gaussian(X)
mu, sigma

(array([ 14.11222578,  14.99771051]), array([ 1.83263141,  1.70974533]))

现在有了我们模型的参数，我们需要确定概率阈值，其指示一个样本是否应该被认为是异常。 为此，我们需要使用一组已被标记的验证数据（其中真正的异常已经被标记出来了），并在不同阈值条件下测试模型识别异常的能力。

Xval = data['Xval']
yval = data['yval']

Xval.shape, yval.shape

((307L, 2L), (307L, 1L))

我们还需要一种方法能够在给定参数下计算数据点属于一个正态分布的概率。幸运的是， SciPy 已经内置了。

from scipy import stats
dist = stats.norm(mu[0], sigma[0])
dist.pdf(X[:,0])[0:50]

array([ 0.183842  ,  0.20221694,  0.21746136,  0.19778763,  0.20858956,
0.21652359,  0.16991291,  0.15123542,  0.1163989 ,  0.1594734 ,
0.21716057,  0.21760472,  0.20141857,  0.20157497,  0.21711385,
0.21758775,  0.21695576,  0.2138258 ,  0.21057069,  0.1173018 ,
0.20765108,  0.21717452,  0.19510663,  0.21702152,  0.17429399,
0.15413455,  0.21000109,  0.20223586,  0.21031898,  0.21313426,
0.16158946,  0.2170794 ,  0.17825767,  0.17414633,  0.1264951 ,
0.19723662,  0.14538809,  0.21766361,  0.21191386,  0.21729442,
0.21238912,  0.18799417,  0.21259798,  0.21752767,  0.20616968,
0.21520366,  0.1280081 ,  0.21768113,  0.21539967,  0.16913173])

在不清楚的情况下，我们只计算了数据集的第一维的前 50 个样本中的每一个属于该分布的概率。 本质上，它计算的是每个实例与均值的距离，以及从均值角度比较这一“典型”距离。

让我们在给定高斯模型参数情况下（上面计算过的）计算并保存数据集中每个值的概率密度.

p = np.zeros((X.shape[0], X.shape[1]))
p[:,0] = stats.norm(mu[0], sigma[0]).pdf(X[:,0])
p[:,1] = stats.norm(mu[1], sigma[1]).pdf(X[:,1])

p.shape

(307L, 2L)

我们需要对验证集进行相同的操作（使用相同的模型参数）。我们将使用这些概率再结合真标签来确定指定为异常点的概率阈值。

pval = np.zeros((Xval.shape[0], Xval.shape[1]))
pval[:,0] = stats.norm(mu[0], sigma[0]).pdf(Xval[:,0])
pval[:,1] = stats.norm(mu[1], sigma[1]).pdf(Xval[:,1])

接下来，我们需要一个只要给定概率密度值和真标签就能找出最佳阈值的函数。 为此，我们需要计算不同的 ϵ 值的 F1 范数。 F1 范数 是真正（ true positive）、假正（false positive）、假负（false negative）的函数。

def select_threshold(pval, yval):
best_epsilon = 0
best_f1 = 0
f1 = 0

step = (pval.max() - pval.min()) / 1000

for epsilon in np.arange(pval.min(), pval.max(), step):
preds = pval < epsilon

tp = np.sum(np.logical_and(preds == 1, yval == 1)).astype(float)
fp = np.sum(np.logical_and(preds == 1, yval == 0)).astype(float)
fn = np.sum(np.logical_and(preds == 0, yval == 1)).astype(float)

precision = tp / (tp + fp)
recall = tp / (tp + fn)
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(pval, yval)
epsilon, f1

(0.0095667060059568421, 0.7142857142857143)

最后，我们将阈值应用到数据集并观察结果。

# indexes of the values considered to be outliers
outliers = np.where(p < epsilon)

fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(X[:,0], X[:,1])
ax.scatter(X[outliers[0],0], X[outliers[0],1], s=50, color='r', marker='o')

不错！ 红色的点是被标记为异常值的点。 从视觉上来看这很合理。 右上那一点有一些分离（但没有被标记）可能也是一个异常值，这已经很接近正确结果了。 应用其到高维数据集是练习文本的另一个例子，但由于它是二维情形的琐碎扩展，我们将进入最后一节。

Collaborative Filtering

推荐引擎使用项和基于用户的相似性度量来检查用户历史偏好以对用户可能感兴趣的新“事物”提出推荐。在本练习中，我们将实现被称为协同过滤的特定推荐算法，并将其应用于电影评分的数据集。 让我们先加载并检查我们将要使用的数据。

data = loadmat('data/ex8_movies.mat')
data

{'R': array([[1, 1, 0, ..., 1, 0, 0],
[1, 0, 0, ..., 0, 0, 1],
[1, 0, 0, ..., 0, 0, 0],
...,
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0]], dtype=uint8),
'Y': array([[5, 4, 0, ..., 5, 0, 0],
[3, 0, 0, ..., 0, 0, 5],
[4, 0, 0, ..., 0, 0, 0],
...,
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0]], dtype=uint8),
'__globals__': [],
'__header__': 'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Thu Dec  1 17:19:26 2011',
'__version__': '1.0'}

Y 是包含从 1 到 5 评分的（电影数 × 用户数）数组。R 是包含指示用户是否对电影评分的二值数据的数组。 两者应具有相同的形状。

Y = data['Y']
R = data['R']
Y.shape, R.shape

((1682L, 943L), (1682L, 943L))

我们可以通过对 Y 中一行里存在评分的数据进行平均来查看电影的平均评分。

Y[1,R[1,:]].mean()

2.5832449628844114

我们还可以尝试通过渲染该矩阵来“可视化”数据，就像它是一个图像一样。 我们不能收集太多，但它给我们一个用户与电影和评分的相对密度之间关系的想法。

fig, ax = plt.subplots(figsize=(12,12))
ax.imshow(Y)
ax.set_xlabel('Users')
ax.set_ylabel('Movies')
fig.tight_layout()

接下来，我们将实现协同过滤的代价函数。 直观地，“代价”是一组电影评分预测偏离真实预测的程度。 代价方程式在练习文本中已经给出。 它基于两组参数矩阵，在文本中被称为 X 和 Theta。 这些被 “unrolled” 到 “params” 输入中，以便我们以后可以使用 SciPy 的优化包。 注意，我已经在注释中包括数组/矩阵形状，以帮助说明矩阵如何交互工作。

def cost(params, Y, R, num_features):
Y = np.matrix(Y)  # (1682, 943)
R = np.matrix(R)  # (1682, 943)
num_movies = Y.shape[0]
num_users = Y.shape[1]

# reshape the parameter array into parameter matrices
X = np.matrix(np.reshape(params[:num_movies * num_features], (num_movies, num_features)))  #
11aee
(1682, 10)
Theta = np.matrix(np.reshape(params[num_movies * num_features:], (num_users, num_features)))  # (943, 10)

# initializations
J = 0

# compute the cost
error = np.multiply((X * Theta.T) - Y, R)  # (1682, 943)
squared_error = np.power(error, 2)  # (1682, 943)
J = (1. / 2) * np.sum(squared_error)

return J

为了测试运行情况，我们提供了一组预训练参数来评估结果。 为了减少评估时间，我们仅看数据的一个小子集。

users = 4
movies = 5
features = 3

params_data = loadmat('data/ex8_movieParams.mat')
X = params_data['X']
Theta = params_data['Theta']

X_sub = X[:movies, :features]
Theta_sub = Theta[:users, :features]
Y_sub = Y[:movies, :users]
R_sub = R[:movies, :users]

params = np.concatenate((np.ravel(X_sub), np.ravel(Theta_sub)))

cost(params, Y_sub, R_sub, features)

22.224603725685675

这个答案符合练习文本中给出的结果。 接下来，我们需要实现梯度计算。 就像我们在练习 4 神经网络中实现的一样，我们将扩展代价函数使其也能计算梯度。

def cost(params, Y, R, num_features):
Y = np.matrix(Y)  # (1682, 943)
R = np.matrix(R)  # (1682, 943)
num_movies = Y.shape[0]
num_users = Y.shape[1]

# reshape the parameter array into parameter matrices
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)

# initializations
J = 0
X_grad = np.zeros(X.shape)  # (1682, 10)
Theta_grad = np.zeros(Theta.shape)  # (943, 10)

# compute the cost
error = np.multiply((X * Theta.T) - Y, R)  # (1682, 943)
squared_error = np.power(error, 2)  # (1682, 943)
J = (1. / 2) * np.sum(squared_error)

# calculate the gradients
X_grad = error * Theta
Theta_grad = error.T * X

# unravel the gradient matrices into a single array
grad = np.concatenate((np.ravel(X_grad), np.ravel(Theta_grad)))

return J, grad

J, grad = cost(params, Y_sub, R_sub, features)
J, grad

(22.224603725685675,
array([ -2.52899165,   7.57570308,  -1.89979026,  -0.56819597,
3.35265031,  -0.52339845,  -0.83240713,   4.91163297,
-0.76677878,  -0.38358278,   2.26333698,  -0.35334048,
-0.80378006,   4.74271842,  -0.74040871, -10.5680202 ,
4.62776019,  -7.16004443,  -3.05099006,   1.16441367,
-3.47410789,   0.        ,   0.        ,   0.        ,
0.        ,   0.        ,   0.        ]))

我们的下一步是将正则化添加到代价和梯度的计算中。 我们将创建函数的最终的正则化版（注意，此版本包括一个额外的学习率参数，称为 “lambda”）。

def cost(params, Y, R, num_features, learning_rate):
Y = np.matrix(Y)  # (1682, 943)
R = np.matrix(R)  # (1682, 943)
num_movies = Y.shape[0]
num_users = Y.shape[1]

# reshape the parameter array into parameter matrices
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)

# initializations
J = 0
X_grad = np.zeros(X.shape)  # (1682, 10)
Theta_grad = np.zeros(Theta.shape)  # (943, 10)

# compute the cost
error = np.multiply((X * Theta.T) - Y, R)  # (1682, 943)
squared_error = np.power(error, 2)  # (1682, 943)
J = (1. / 2) * np.sum(squared_error)

# add the cost regularization
J = J + ((learning_rate / 2) * np.sum(np.power(Theta, 2)))
J = J + ((learning_rate / 2) * np.sum(np.power(X, 2)))

# calculate the gradients with regularization
X_grad = (error * Theta) + (learning_rate * X)
Theta_grad = (error.T * X) + (learning_rate * Theta)

# unravel the gradient matrices into a single array
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

(31.344056244274221,
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]))

这个结果再次与练习文本代码的预期输出相匹配，所以看起来正则化是正常工作的。 在我们训练模型之前，我们还有一个最后一步。 我们的任务是创建我们自己的电影评分，以便我们可以使用该模型生成个性化的推荐。 我们有一个将电影索引和其标题相连接的文件。 让我们将文件加载到字典中，并使用练习中提供的一些样本评分。

movie_idx = {}
f = open('data/movie_ids.txt')
for line in f:
tokens = line.split(' ')
tokens[-1] = tokens[-1][:-1]
movie_idx[int(tokens[0]) - 1] = ' '.join(tokens[1:])

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]))))

Rated Toy Story (1995) with 4 stars.
Rated Twelve Monkeys (1995) with 3 stars.
Rated Usual Suspects, The (1995) with 5 stars.
Rated Outbreak (1995) with 4 stars.
Rated Shawshank Redemption, The (1994) with 5 stars.
Rated While You Were Sleeping (1995) with 3 stars.
Rated Forrest Gump (1994) with 5 stars.
Rated Silence of the Lambs, The (1991) with 2 stars.
Rated Alien (1979) with 4 stars.
Rated Die Hard 2 (1990) with 5 stars.
Rated Sphere (1998) with 5 stars.

我们可以将这个自定义评分向量添加到数据集，以便它能够包含在模型中。

R = data['R']
Y = data['Y']

Y = np.append(Y, ratings, axis=1)
R = np.append(R, ratings != 0, axis=1)

我们现在准备训练协同过滤模型。 我们将归一化电影评分，然后使用我们的成本函数，参数向量和输入端的数据矩阵等进入优化例程。

from scipy.optimize import minimize

movies = Y.shape[0]
users = Y.shape[1]
features = 10
learning_rate = 10.

X = np.random.random(size=(movies, features))
Theta = np.random.random(size=(users, features))
params = np.concatenate((np.ravel(X), np.ravel(Theta)))

Ymean = np.zeros((movies, 1))
Ynorm = np.zeros((movies, users))

for i in range(movies):
idx = np.where(R[i,:] == 1)[0]
Ymean[i] = Y[i,idx].mean()
Ynorm[i,idx] = Y[i,idx] - Ymean[i]

fmin = minimize(fun=cost, x0=params, args=(Ynorm, R, features, learning_rate),
method='CG', jac=True, options={'maxiter': 100})
fmin

status: 1
success: False
njev: 149
nfev: 149
fun: 38953.88249706676
x: array([-0.07177334, -0.08315075,  0.1081135 , ...,  0.1817828 ,
0.16873062,  0.03383596])
message: 'Maximum number of iterations has been exceeded.'
jac: array([ 0.01833555,  0.07377974,  0.03999323, ..., -0.00970181,
0.00758961, -0.01181811])

由于为了让优化程序正常工作一切都被 “unrolled”了 ，我们需要 reshape 我们的矩阵使它们回到原始尺寸。

X = np.matrix(np.reshape(fmin.x[:movies * features], (movies, features)))
Theta = np.matrix(np.reshape(fmin.x[movies * features:], (users, features)))

X.shape, Theta.shape

((1682L, 10L), (944L, 10L))

我们训练后的参数现在在 X 和 Theta 中了。 我们可以使用这些来为我们之前添加的用户给出一些推荐。

predictions = X * Theta.T
my_preds = predictions[:, -1] + Ymean
sorted_preds = np.sort(my_preds, axis=0)[::-1]
sorted_preds[:10]

matrix([[ 5.00000264],
[ 5.00000249],
[ 4.99999831],
[ 4.99999671],
[ 4.99999659],
[ 4.99999253],
[ 4.99999238],
[ 4.9999915 ],
[ 4.99999019],
[ 4.99998643]]

这给出一个有序的评分列表，但我们失去了这些评分的索引。 实际上我们需要使用 argsort，以便让我们知道预测评分对应的电影。

idx = np.argsort(my_preds, axis=0)[::-1]
print("Top 10 movie predictions:")
for i in range(10):
j = int(idx[i])
print('Predicted rating of {0} for movie {1}.'.format(str(float(my_preds[j])), movie_idx[j]))

Top 10 movie predictions:
Predicted rating of 5.00000264002 for movie Prefontaine (1997).
Predicted rating of 5.00000249142 for movie Santa with Muscles (1996).
Predicted rating of 4.99999831018 for movie Marlene Dietrich: Shadow and Light (1996) .
Predicted rating of 4.9999967124 for movie Saint of Fort Washington, The (1993).
Predicted rating of 4.99999658864 for movie They Made Me a Criminal (1939).
Predicted rating of 4.999992533 for movie Someone Else's America (1995).
Predicted rating of 4.99999238336 for movie Great Day in Harlem, A (1994).
Predicted rating of 4.99999149604 for movie Star Kid (1997).
Predicted rating of 4.99999018592 for movie Aiqing wansui (1994).
Predicted rating of 4.99998642746 for movie Entertaining Angels: The Dorothy Day Story (1996).

推荐的电影与练习文本中给出的不一致。 原因是什么还不太清楚，我没有找到任何理由来解释它，可能在某处的代码中有错误。 如果有人发现错误并指出会有奖励。 不过，这仅是一些细微差别，这个例子的大部分是准确的。

最后一次练习结束了！ 当我开始这个系列时，我的目标是在 Andrew 的课程中更熟练的掌握 python，以及凝练机器学习知识。 我确信我实现了这一目标。 我希望这些文章是有价值的阅读材料，因为它是为我创造的。