Python機器學習的練習八：異常檢測和推薦系統
阿新 • • 發佈：2022-05-02
在這篇文章中，將會涉及兩個話題——異常檢測和推薦系統，我們將使用高斯模型實現異常檢測演算法並且應用它檢測網路上的故障伺服器。我們還將看到如何使用協同過濾建立推薦系統，並將其應用於電影推薦資料集。
異常檢測
我們的第一個任務是利用高斯模型判斷資料集裡未標記的例子是否應該被認為是異常的。我們可以在簡單的二維資料集上開始，這樣就可以很容易的看到演算法如何工作。
載入資料並繪圖。
import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sb from scipy.ioimport 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 scipyimport 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])接下來我們需要一個函式，在給出的概率密度和真實標籤中找到最佳閾值。為了執行這步操作，我們需要計算不同值的epsilon的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 epsilonin 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')
結果還不錯，紅色的點是被標記為離群值的點，視覺上這些看起來很合理。有一些分離(但沒有被標記)的右上角的點也可能是一個離群值，但這是相當接近的。
協同過濾
推薦系統使用專案和基於使用者的相似性計算，以檢查使用者的歷史偏好，從而為使用者推薦可能感興趣的新“東西”。在這個練習中，我們將實現一種特殊的推薦演算法，稱為協同過濾，並將其應用於電影評分的資料集。首先載入並檢查我們要處理的資料。
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是一個包含從一到五的評分的資料組，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引數矩陣為基礎，這些展開到引數輸入中，以便使用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))) # (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答案與練習文字中的結果匹配。接下來需要實現梯度計算，與練習四中神經網路的實現一樣，我們會擴充套件成本函式計算梯度。
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 linein 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) with4 stars.
Rated Twelve Monkeys (1995) with3 stars.
Rated Usual Suspects, The (1995) with5 stars.
Rated Outbreak (1995) with4 stars.
Rated Shawshank Redemption, The (1994) with5 stars.
Rated While You Were Sleeping (1995) with3 stars.
Rated Forrest Gump (1994) with5 stars.
Rated Silence of the Lambs, The (1991) with2 stars.
Rated Alien (1979) with4 stars.
Rated Die Hard2 (1990) with5 stars.
Rated Sphere (1998) with5 stars.我們可以在資料集中新增自定義評分向量。
R= data['R']
Y= data['Y']
Y= np.append(Y, ratings, axis=1)
R= np.append(R, ratings != 0, axis=1)開始訓練協同過濾模型，我們將通過成本函式、引數向量和輸入的資料矩陣使評分正規化，並且執行優化程式。
from scipy.optimizeimport 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 iin 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”，因此要正確地工作，需要將我們的矩陣重新調整回原來的維度。
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 iin range(10):
j= int(idx[i])
print('Predicted rating of {0} for movie {1}.'.format(str(float(my_preds[j])), movie_idx[j]))
Top10 movie predictions:
Predicted rating of5.00000264002 for movie Prefontaine (1997).
Predicted rating of5.00000249142 for movie Santa with Muscles (1996).
Predicted rating of4.99999831018 for movie Marlene Dietrich: Shadowand Light (1996) .
Predicted rating of4.9999967124 for movie Saint of Fort Washington, The (1993).
Predicted rating of4.99999658864 for movie They Made Me a Criminal (1939).
Predicted rating of4.999992533 for movie Someone Else's America (1995).
Predicted rating of4.99999238336 for movie Great Dayin Harlem, A (1994).
Predicted rating of4.99999149604 for movie Star Kid (1997).
Predicted rating of4.99999018592 for movie Aiqing wansui (1994).
Predicted rating of4.99998642746 for movie Entertaining Angels: The Dorothy Day Story (1996).實際上推薦的電影並沒有很好地符合練習文字中的內容。原因不太清楚，我還沒有找到任何可以解釋的理由，在程式碼中可能有錯誤。不過，即使有一些細微的差別，這個例子的大部分也是準確的。
最後的練習結束了!