ここでは，集計済みの行列形式のレイティングデータを読み込んで顧客間の類似度に基づいた協調フィルタリングを実装してみます．まず，Pandas と NumPy ライブラリをインポートします．また，Numpyでの小数点以下の表示桁数を調整するとともに，指数表示にならないように出力形式を設定しておきます．
import pandas as pd
import numpy as np
np.set_printoptions(precision=3)
np.set_printoptions(suppress=True)
次に，GitHub のリポジトリから行列形式のレイティングデータをデータフレームに読み込みます．
url = "https://github.com/rinsaka/sample-data-sets/blob/master/collaborative_filtering_matrix.csv?raw=true"
# url = 'collaborative_filtering_matrix.csv' # ローカルにダウンロードしてから開く場合
df = pd.read_csv(url)
print(df)
id name A B C D E 0 0 Eto 5.0 4.0 3.0 NaN NaN 1 1 Sato NaN 2.0 4.0 2.0 3.0 2 2 Kato 2.0 3.0 4.0 NaN NaN 3 3 Muto 2.0 3.0 2.0 NaN 2.0 4 4 Kito 4.0 3.0 3.0 4.0 NaN 5 5 Goto 2.0 5.0 NaN 2.0 3.0 6 6 Bito 3.0 3.0 2.0 NaN 1.0 7 7 Saito 4.0 4.0 5.0 NaN 5.0 8 8 Naito 5.0 5.0 4.0 2.0 1.0 9 9 Koto 3.0 NaN 2.0 2.0 2.0
目的は顧客 Kato（ と Eto ）に対して商品「D」と商品「E」のどちらを優先的におすすめすべきかを決めることです．顧客の好みが異なれば自ずとおすすめする商品も異なるはずです．
ここでは，（古典的ですが）顧客の類似度に基づいた協調フィルタリングを考えます．この手法ではおおむね次のようなステップで顧客に対するおすすめ商品を決定します．
まずは，「誰」の「どの商品」について推定したいかを指定するパラメータを設定しておきます．以降では，顧客 Kato の商品「D」のスコアを求めてみます．なお，次のコードのいずれかの行を有効にするだけで，別のスコアを得ることができるようになるはずです．
## customer_id は 推定したい顧客（またはユーザ）
## item_id は 推定したい商品（またはアイテム）
# customer_id = 0; item_id = 3;
# customer_id= 0; item_id = 4;
customer_id= 2; item_id = 3;
# customer_id= 2; item_id = 4;
近傍顧客として利用する顧客の人数 \(k\) もここで設定しておきます（実際上は適切な \(k\) がどの程度であるかについて検討する必要があります）．
k = 3
すでに読み込んでいる Pandas のデータフレームからレイティングに関する列だけを取得して NumPy 配列に格納し，レイティング行列 y を生成します．
y = df.iloc[:, 2:].values
print(y)
[[ 5. 4. 3. nan nan] [nan 2. 4. 2. 3.] [ 2. 3. 4. nan nan] [ 2. 3. 2. nan 2.] [ 4. 3. 3. 4. nan] [ 2. 5. nan 2. 3.] [ 3. 3. 2. nan 1.] [ 4. 4. 5. nan 5.] [ 5. 5. 4. 2. 1.] [ 3. nan 2. 2. 2.]]
もちろん，次のような書き方でも動作しますが，列名の指定は面倒かもしれません．
y = df.loc[:, ['A', 'B', 'C', 'D', 'E']].values
顧客 \(i\) と好みが近い顧客を \(k\) 人見つけるために，顧客 \(i\) と \(\ell\) の類似度 \( sim(i, \ell)\) を定式化します．これは顧客 \(i\) と \(\ell\) のレイティングの相関係数です． \begin{eqnarray} sim(i,\ell) &=& \frac{\sum_{j \in \mathcal{J}_{i\ell}}\left(y_{ij} - \overline{y}_{i.} \right)\left( y_{\ell j} - \overline{y}_{\ell.} \right)} { \sqrt{\sum_{j \in \mathcal{J}_{i\ell}} \left( y_{ij} - \overline{y}_{i.} \right)^2 } \sqrt{\sum_{j \in \mathcal{J}_{i\ell}} \left( y_{\ell j} - \overline{y}_{\ell.} \right)^2 } } \end{eqnarray} ここで，\(\mathcal{J}_{i\ell}\) は，顧客 \(i\)と \(\ell\) の両方によって評価された商品の集合を意味します．また，\(y_{ij}\) は 顧客 \(i\) の商品 \(j\) に対するレイティング（評価値）を意味し，\(\overline{y}_{i.}\) は 顧客 \(i\) の平均レイティングですが，これは \(\mathcal{J}_{i\ell}\) での平均とします．つまり，両方の顧客で評価された商品だけでの平均レイティングです．
上の \( sim(i, \ell)\) を求めるコードを記述していきます．まず，上で定義した
customer_id とitem_id を定式化での記号に合わせて内部的には i, と j として利用することにします．
i = customer_id
j = item_id
まず，顧客 \(i\) と顧客 \(\ell\) の両方によって評価された評価の集合を返す関数 get_Y_il() を定義します．
def get_Y_il(y, i, l):
# まずは，yi と yl だけを取り出して２行の行列を生成
yil = np.concatenate((y[i], y[l])).reshape(2,y.shape[1])
# print(yil)
# NaN が含まれる列のインデックスを初期化（あとでNaNが含まれている列は削除したい）
nan_idx = []
for i in range(yil.shape[1]):
if np.isnan(yil[0][i]) or np.isnan(yil[1][i]):
nan_idx.append(i) # どちらかの顧客が未評価 (NaN) なら削除する列のリストに追加する
# NaN が含まれる列を削除する
return np.delete(yil, nan_idx, 1)
例として，顧客 8 と顧客 9 の両方によって評価された評価の集合を表示してみます．行列から計算に必要な部分だけが抽出されていることがわかります．
y_il = get_Y_il(y, 8, 9)
print(y)
print('---')
print(y_il)
[[ 5. 4. 3. nan nan] [nan 2. 4. 2. 3.] [ 2. 3. 4. nan nan] [ 2. 3. 2. nan 2.] [ 4. 3. 3. 4. nan] [ 2. 5. nan 2. 3.] [ 3. 3. 2. nan 1.] [ 4. 4. 5. nan 5.] [ 5. 5. 4. 2. 1.] [ 3. nan 2. 2. 2.]] --- [[5. 4. 2. 1.] [3. 2. 2. 2.]]
顧客 \(i\) と 顧客 \(\ell\) の類似度を返す関数 get_sim_i_l() を定義します．この関数の中で，顧客ごとの平均レイティングを求めていますが，両方の顧客がともに評価した商品だけで平均を求めていることに注意してください．
def get_sim_i_l(y, i, l):
y_il = get_Y_il(y, i, l)
mean0 = np.mean(y_il[0])
mean1 = np.mean(y_il[1])
term1 = np.sum((y_il[0] - mean0) * (y_il[1] - mean1))
term2 = np.sqrt(np.sum((y_il[0] - mean0)**2)) * np.sqrt(np.sum((y_il[1] - mean1)**2))
if term2 == 0:
return 0
return term1 / term2
準備ができたので，例として顧客 8 と顧客 9 の類似度を計算してみます．この結果，顧客 8 と顧客 9 の類似度が 0.730 と高く，2人はよく似た好みを持つ顧客であることがわかりました．
get_sim_i_l(y, 8, 9)
0.7302967433402214
この結果が正しいことを数式を用いて確認してみます．その前に，顧客 8 と 顧客 9 が共通して評価した値を取り出し，その中での平均を求めておきます．
print(y_il)
print(np.mean(y_il, axis=1))
[[5. 4. 2. 1.] [3. 2. 2. 2.]] [3. 2.25]
平均が 3 と 2.25 であることがわかったので，この値を利用して計算してみます． \begin{eqnarray} sim(8,9) &=& \frac{\sum_{j \in \mathcal{J}_{8\ell}}\left(y_{8j} - \overline{y}_{8.} \right)\left( y_{9 j} - \overline{y}_{9.} \right)} { \sqrt{\sum_{j \in \mathcal{J}_{8,9}} \left( y_{8j} - \overline{y}_{8.} \right)^2 } \sqrt{\sum_{j \in \mathcal{J}_{8,9}} \left( y_{9j} - \overline{y}_{9.} \right)^2 } } \\ &=& \frac{ \left(5 - 3\right)\left(3 - 2.25\right) + \left(4 - 3\right)\left(2 - 2.25\right) + \left(2 - 3\right)\left(2 - 2.25\right) + \left(1 - 3\right)\left(2 - 2.25\right) } { \sqrt{ \left( 5-3 \right)^2 + \left( 4-3 \right)^2 + \left( 2-3 \right)^2 + \left( 1-3 \right)^2 } \times \sqrt{ \left( 3-2.25 \right)^2 + \left( 2-2.25 \right)^2 + \left( 2-2.25 \right)^2 + \left( 2-2.25 \right)^2 } } \\ &=& \frac{2} { \sqrt{10} \times \sqrt{0.75} } \\ &\approx& \frac{2}{3.162 \times 0.866} \approx \frac{2}{2.738} \\ &\approx& 0.730 \end{eqnarray}
一方で顧客 6 と顧客 7 の類似度は -1 に近いことから，逆の好みを持つこともわかります．
get_sim_i_l(y, 6, 7)
-0.9045340337332909
すべての顧客の組み合わせについて類似度を計算します．このとき，顧客 \(i\) と 顧客 \(i\) の類似度は当然 1 になることに注意してください．
sim_matrix = np.zeros((y.shape[0], y.shape[0]))
for u in range(sim_matrix.shape[0]):
for v in range(sim_matrix.shape[1]):
sim_matrix[u, v] = get_sim_i_l(y, u, v)
print(sim_matrix)
[[ 1. -1. -1. 0. 0.866 -1. 0.866 -0.866 0.866 1. ] [-1. 1. 1. -0.866 -0.5 -0.189 -0.5 0.866 0. 0. ] [-1. 1. 1. 0. -0.866 1. -0.866 0.866 -0.866 -1. ] [ 0. -0.866 0. 1. -0.5 0.945 0.522 -0.577 0.44 0. ] [ 0.866 -0.5 -0.866 -0.5 1. -1. 0.5 -0.5 -0.408 0.5 ] [-1. -0.189 1. 0.945 -1. 1. 0.189 -0.189 0.343 -0.5 ] [ 0.866 -0.5 -0.866 0.522 0.5 0.189 1. -0.905 0.966 0.866] [-0.866 0.866 0.866 -0.577 -0.5 -0.189 -0.905 1. -0.762 -1. ] [ 0.866 0. -0.866 0.44 -0.408 0.343 0.966 -0.762 1. 0.73 ] [ 1. 0. -1. 0. 0.5 -0.5 0.866 -1. 0.73 1. ]]
上で得られた類似度の高い順に本人を除いて \(n\) 人（本人も含めて \(n+1\)人）を抽出すれば近傍顧客のリストが出来上がります．このために類似度をソートして，要素のインデックスを取得します．このとき，np.argsort() は昇順ソートになるために，逆順にする必要があります．このために [::-1] をつけてソート結果を逆順の降順ソートにします．例えば，i = 2 の Kato さんについて取り出します．この結果の意味は，「0位はインデックス5，1位はインデックス2，2位はインデックス1，3位がインデックス7・・・」となります．つまり，1 位を除いた，0 位から n = 3 位までが近傍顧客となります．
np.argsort(sim_matrix[i])[::-1]
array([5, 2, 1, 7, 3, 8, 6, 4, 9, 0])
上の np.argsort を重ねると順位が得られることになります．つまり，「インデックス5の Goto は0位（Katoに最も近い），インデックス2の Kato（本人）は1位，インデックス1の Sato は2位，インデックス7の Saito は3位・・・」であることを意味しています．
customer_rank = np.argsort(np.argsort(sim_matrix[i])[::-1])
print(customer_rank)
[9 2 1 4 7 0 6 3 5 8]
（類似度が 1 のはずの） Kato自身も含んだ近傍顧客のインデックスを k+1 人分取得します．つまり，インデックス 1, 2, 5, 7 の Sato, Kato, Goto, Saito が近傍顧客です．
idx = np.arange(len(sim_matrix[i]))
near_customer_ids = idx[customer_rank <= k]
print(near_customer_ids)
[1 2 5 7]
スコアを計算するために，類似度の sim ベクトルから近傍顧客の類似度だけを取り出します．
near_sim = sim_matrix[i, near_customer_ids]
print(near_sim)
[1. 1. 1. 0.866]
同様に，レイティング行列から近傍顧客のレイティングだけを取り出します．
near_y = y[near_customer_ids]
print(near_y)
[[nan 2. 4. 2. 3.] [ 2. 3. 4. nan nan] [ 2. 5. nan 2. 3.] [ 4. 4. 5. nan 5.]]
評価したい顧客の変数 i が近傍顧客のみの集団での customer_id のインデックスとなるように i を更新します．これは，近傍顧客の類似度ベクトルやレイティング行列のインデックス 1 がもとの顧客 2 の情報であることを意味しています．
i = np.where(near_customer_ids == i)[0][0]
print(i)
1
顧客 \(i\) の商品 \(j\) に対するスコアは次の式から求められます． \begin{eqnarray} s_{ij} &=& \overline{y}_{i.} + \frac { \sum_{\ell \in \mathcal{I}_{j}(i)}w(i, \ell) \left( y_{\ell j} - \overline{y}_{\ell .}\right) } { \sum_{\ell \in \mathcal{I}_{j}(i)} \left| w(i, \ell) \right| } \end{eqnarray} ただし，\(\mathcal{I}_{j}(i)\) は商品 \(j\) を評価した顧客 \(i\) に類似する顧客の集合です．また \(w(i,\ell)\) は顧客 \(i\) が商品 \(j\) を評価するときの顧客 \(\ell\) の評価に対する重みであり，ここでは \(w(i,\ell) = sim(i, \ell)\) とします．さらに，\(\overline{y}_{i.}\) は顧客 \(i\) の平均レイティングです．
実際に Kato さんの商品「D」についてのスコアは次のようになります．このとき，近傍顧客の中で商品「D」をまだ評価していない Saito さんのデータは利用されないことに注意してください． \begin{eqnarray} s_{Kato,3} &=& 3 + \frac{ 1 \times(2-2.75) + 1 \times (2 -3) }{ 1 + 1 } \\ &=& 3 + \frac{-0.75-1}{2} \\ &=& 3 - 0.875 \\ &=& 2.125 \end{eqnarray}
上の式を用いてスコアを求める関数 get_sij() を次のように定義します．
def get_sij(y, w, i, j):
y_mean = np.nanmean(y, axis=1)
y_mean_i = y_mean[i]
# NaNが含まれている行，および自身の行は削除したい
nan_idx = [i]
for k in range(len(w)):
if np.isnan(y[k][j]):
nan_idx.append(k)
y = np.delete(y, nan_idx, 0)
y_mean = np.delete(y_mean, nan_idx)
w = np.delete(w, nan_idx)
# j 列目だけを取り出す
y = y[:,j]
sij = y_mean_i + np.sum(w * (y - y_mean)) / np.sum(w)
return sij
近傍顧客のレイティング行列，および類似度を与えて，顧客 i=1（customer_id = 2）の商品 j=3（商品「D」） に対するスコアを求めます．数式上で求めたものと同じ結果が得られました．
sij = get_sij(near_y, near_sim, i, j)
print(sij)
2.125
結果を表示してみます．
print(df)
print(f"商品 {item_id} ({df.columns[item_id+2]}) に対する顧客 {customer_id} のスコアは {sij:5.3f} です．")
id name A B C D E 0 0 Eto 5.0 4.0 3.0 NaN NaN 1 1 Sato NaN 2.0 4.0 2.0 3.0 2 2 Kato 2.0 3.0 4.0 NaN NaN 3 3 Muto 2.0 3.0 2.0 NaN 2.0 4 4 Kito 4.0 3.0 3.0 4.0 NaN 5 5 Goto 2.0 5.0 NaN 2.0 3.0 6 6 Bito 3.0 3.0 2.0 NaN 1.0 7 7 Saito 4.0 4.0 5.0 NaN 5.0 8 8 Naito 5.0 5.0 4.0 2.0 1.0 9 9 Koto 3.0 NaN 2.0 2.0 2.0 商品 3 (D) に対する顧客 2 のスコアは 2.125 です．
今の例では，近傍顧客数を n=3 と設定したものの，近傍顧客の中で商品「D」をまだ評価していない顧客がいたため，実際には2名の近傍顧客のデータしか利用しませんでした．大きな顧客数と商品数で構成される現実の評価データでは，巨大な行列の中のごく一部にしか評価値がない（評価値が全体の10%にも満たず，90%以上が欠損値である）ことが少なくありません．このような場合には，商品「D」を評価していない顧客は近傍顧客に含めないようにしたり，そのような顧客はあらかじめ削除してから近傍顧客を求めるなどの工夫が必要になるでしょう．