ベイズ勉強会 Part 9 非負値行列因子分解
ベイズ推論による非負値行列因子分解
ベイズ勉強会資料は『ベイズ推論による機械学習入門』1を元に、途中式計算をできるだけ省略せずに行ったものです。
今回は5.2の内容を扱います。
import numpy as np
np.random.seed(42)
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
データのimportと前処理
Kaggleにあるポケモン(第6世代まで)のデータセットを使う。
英語名だとわかりにくいので日本語訳する。 対応表はポケモン徹底攻略から拝借した。
日本語訳したものをgistに置いてある。
作成過程
#collapse-hide
pokemon = pd.read_csv('data/for0714/pokemon.csv')
# メガシンカは前に、フォルムチェンジは後ろについているので、メガシンカのみ2列目を使う
split = pokemon['Name'].str.split(' ', expand=True)
# 効率的な書き方ではないが行数少ないので素朴に
English = []
for i in range(800):
if split.iloc[i,0] == 'Mega':
English.append(split.iloc[i,1])
else:
English.append(split.iloc[i,0])
pokemon['English'] = English
# 日本語訳
translate = pd.read_csv('data/for0714/japanese.csv')
pokemon = pd.merge(pokemon, translate, on = 'English', how = 'left')
# 以下特殊事例の処理
# オコリザル
pokemon.loc[62,'Name']='Primeape'
pokemon.loc[62,'Pokedex_No']=57
pokemon.loc[62,'Japanese']='オコリザル'
# バリヤード
pokemon.loc[131,'Pokedex_No']=122
pokemon.loc[131,'Japanese']='バリヤード'
# ホウオウ
pokemon.loc[270,'Pokedex_No']=250
pokemon.loc[270,'Japanese']='ホウオウ'
# ゲンシカイキ
pokemon.loc[422,'Pokedex_No']=382
pokemon.loc[422,'Japanese'] = 'カイオーガ'
pokemon.loc[424,'Pokedex_No']=383
pokemon.loc[424,'Japanese'] = 'グラードン'
# デオキシス
pokemon.loc[429,'Pokedex_No']=386
pokemon.loc[429,'Japanese']='デオキシス'
# マネネ
pokemon.loc[487,'Pokedex_No']=439
pokemon.loc[487,'Japanese']='マネネ'
# ロトム
pokemon.loc[532:537,'Pokedex_No']=479
pokemon.loc[532:537,'Japanese']='ロトム'
# フラベベ
pokemon.loc[737,'Pokedex_No']=669
pokemon.loc[737,'Japanese']='フラベベ'
# 合計種族値の列も作る
pokemon['sum'] = pokemon.iloc[:,4:10].sum(axis=1)
# 出力
# pokemon.drop(['#','English'],axis=1).to_csv('pokemon_data_g6.csv',index=False)
pokemon = pd.read_csv('https://gist.githubusercontent.com/Vintea01/'
'9f83dd0b54dd0360067a08fad167fe98/raw/'
'b7f7f1253515f0151cf5428330fcab0ed26900c8/pokemon_data_g6.csv')
# 種族値だけ取り出す
X = pokemon.iloc[:,3:9]
X
vmin = 0
vmax = np.max(np.max(X))
sns.heatmap(X, vmin=vmin, vmax=vmax)
from sklearn.decomposition import NMF
def fitNMF(X,M):
model = NMF(n_components=M, init='random', random_state=0)
W = model.fit_transform(X)
H = model.components_
return (W,H)
# M=1でNMF
WH1 = fitNMF(X,1)
fig, axes = plt.subplots(2,1, figsize=(6, 4))
# データセットでの通し番号を横軸にしてプロット
axes[0].plot(range(800),WH1[0])
# 図鑑番号を横軸にしてプロット
axes[1].plot(pokemon['Pokedex_No'],WH1[0])
# 各世代の最初と終盤(御三家や600族、伝説がいるあたり)で高くなっているのがわかる
# 特徴量1つなので、合計種族値にほぼ比例する
plt.scatter(pokemon.iloc[:,3:9].sum(axis=1),WH1[0])
pokemon['W1_1'] = WH1[0]
pokemon.loc[np.argsort(- WH1[0], axis=0).flatten()[:10],:]
# ヒートマップで確認
fig, axes = plt.subplots(1,2, figsize=(10,4))
sns.heatmap(X, ax=axes[0],vmin=vmin,vmax=vmax)
sns.heatmap(np.dot(WH1[0],WH1[1]), ax=axes[1], vmin=vmin, vmax=vmax)
# M=2
WH2 = fitNMF(X,2)
pokemon['W2_1'] = WH2[0][:,0]
pokemon['W2_2'] = WH2[0][:,1]
pokemon['W2_norm'] = np.sqrt(WH2[0][:,0]**2 + WH2[0][:,1]**2) # 原点からの距離
fig, axes = plt.subplots(3,1, figsize=(6, 6))
axes[0].plot(pokemon['Pokedex_No'],WH2[0][:,0])
axes[1].plot(pokemon['Pokedex_No'],WH2[0][:,1])
axes[2].plot(pokemon['Pokedex_No'],pokemon['W2_norm'])
# W2_1の1位はメガハガネール
pokemon.loc[np.argsort(- WH2[0][:,0])[:10],:]
# W2_2の1位はデオキシスアタックフォルム
pokemon.loc[np.argsort(- WH2[0][:,1])[:10],:]
上は防御系、下は攻撃系が考慮されていることがわかる。
# 原点からの距離上位10
pokemon.loc[np.argsort(- np.sqrt(WH2[0][:,0]**2 + WH2[0][:,1]**2))[:10],:]
# つよいポケモン よわいポケモン
# そんなの ひとの かって
# ほんとうに つよい トレーナーなら
# すきなポケモンで かてるように がんばるべき
pokemon.loc[np.argsort(np.sqrt(WH2[0][:,0]**2 + WH2[0][:,1]**2))[:10],:]
# ヒートマップで確認
fig, axes = plt.subplots(1,3, figsize=(15,4))
sns.heatmap(X, ax=axes[0],vmin=vmin,vmax=vmax)
sns.heatmap(np.dot(WH1[0],WH1[1]), ax=axes[1],vmin=vmin,vmax=vmax)
sns.heatmap(np.dot(WH2[0],WH2[1]), ax=axes[2],vmin=vmin,vmax=vmax)
レコメンドシステム用の関数も書いてみる。
# ポケモン同士の距離を計算
def calc_distance(data, pokemonA, pokemonB):
A = np.array(data.loc[pokemonA, ['W2_1','W2_2']]).flatten()
B = np.array(data.loc[pokemonB, ['W2_1','W2_2']]).flatten()
return np.linalg.norm(A-B)
# 距離が近い上位n種を推薦
def recommendation(data, base_pokemon, n):
distances = np.array([calc_distance(data,base_pokemon, p) for p in range(800)])
index = np.argsort(distances)[:n+1]
return data.loc[index]
# フシギダネに近いポケモン
recommendation(pokemon,0,10)
種族値しか考慮されていないのでレコメンドシステムとしては使い物にならないが
ベイズ推論で解く
散々遊んだのでNMFのベイズ推論による解法の1例を見てみる。ポケモンの種族値は全て自然数なので、ポアソン分布を使うことができる。
モデル定義
目標
行列${\bf X} \in \mathbb{N}^{D \times N}$を2つの行列${\bf W} \in \mathbb{R}^{+ D \times M}$と${\bf H} \in \mathbb{R}^{+ M \times N}$の積に近似する。
$${\bf X} \approx {\bf W H}$$
Notes:ここでの$\mathbb{N}$は0を含んだ非負の整数である。
登場する変数
- 元の行列: ${\bf X} \in \mathbb{N}^{D \times N}$
- 分解先の行列: ${\bf W} \in \mathbb{R}^{+ D \times M}$,${\bf H} \in \mathbb{R}^{+ M \times N}$
- 潜在変数: ${\bf S} \in \mathbb{N}^{D \times M \times N}$
ここで、${\bf X}$を成分で見たとき
$$ \begin{eqnarray} X_{d,n} &=& \Sigma_{m=1}^M S_{d,m,n} \\ &\approx& \Sigma_{m=1}^M W_{d,m} H_{m,n} \end{eqnarray} $$となるような潜在変数${\bf S} \in \mathbb{N}^{D \times M \times N}$をおいた。このような中間変数を導入することで効率的な近似アルゴリズムを導出できることが知られている。
それぞれの確率分布
${\bf S}$の分布
潜在変数${\bf S}$は${\bf W},{\bf H}$の成分の積をパラメータに持つポアソン分布に従って発生するとする。
$$ \begin{eqnarray} p({\bf S}|{\bf W},{\bf H}) &=& \Pi_{d=1}^D \Pi_{m=1}^M \Pi_{n=1}^N p(S_{d,m,n}|W_{d,m},H_{m,n}) \\ &=& \Pi_{d=1}^D \Pi_{m=1}^M \Pi_{n=1}^N \mathrm{Poi}(S_{d,m,n}|W_{d,m} H_{m,n}) \end{eqnarray} $$${\bf X}$の分布
${\bf X}$を${\bf S}$に条件づけられたデルタ分布に従うとする。
$$ \begin{eqnarray} p({\bf X}|{\bf S}) &=& \Pi_{d=1}^D \Pi_{n=1}^N p(X_{d,n}|S_{d,:,n}) \\ &=& \Pi_{d=1}^D \Pi_{n=1}^N \mathrm{Del}(X_{d,n}|\Sigma_{m=1}^M S_{d,m,n}) \end{eqnarray} $$デルタ分布は次のような離散値に対する確率分布である。
$$ \mathrm{Del}(x|z) = \begin{cases} 1 \mathrm{if} x=z \\ 0 \mathrm{otherwise} \end{cases} $$$\Sigma_{m=1}^M S_{d,m,n}$を$X_{d,n}$にそのまま写す確率分布である。
${\bf W},{\bf H}$の分布
各成分がポアソン分布の共役事前分布であるガンマ分布に独立に従うとする。
$$ \begin{eqnarray} p({\bf W}) &=& \Pi_{d=1}^D \Pi_{m=1}^M p(W_{d,m}) \\ &=& \Pi_{d=1}^D \Pi_{m=1}^M \mathrm{Gam}(W_{d,m}|a_W,b_W) \end{eqnarray} $$$$ \begin{eqnarray} p({\bf H}) &=& \Pi_{m=1}^M \Pi_{n=1}^N p(H_{m,n}) \\ &=& \Pi_{m=1}^M \Pi_{n=1}^N \mathrm{Gam}(H_{m,n}|a_H,b_H) \end{eqnarray} $$$a_W,b_W,a_H,b_H$は正の実数であり、ガンマ分布のハイパーパラメータである。
同時分布
確率変数の関係をまとめるとNMFに用いるモデルは次のようになる。
$$ p({\bf X},{\bf S},{\bf W},{\bf H}) = p({\bf X}|{\bf S})p({\bf S}|{\bf W},{\bf H})p({\bf W})p({\bf H}) $$これをDAGで表すと次のようになる。
Important: 事後分布に対する近似公式 $$\begin{eqnarray} \ln q(z_i) &=& \langle \ln p(z_1,\dots,z_M|\mathcal{D}) \rangle_{q({\bf Z}_{\verb|\| i} ) }+const. \\ &=& \langle \ln p(\mathcal{D},z_1,\dots,z_M) \rangle_{q({\bf Z}_{\verb|\| i} ) } + const. \end{eqnarray}$$
近似公式を用いて更新式を書く。
$$ \begin{eqnarray} \ln q({\bf W}) &=& \langle \ln p({\bf S}|{\bf W},{\bf H}) \rangle_{q({\bf S})q({\bf H})} + \ln p({\bf W}) + const. \\ &=& \Sigma_{d=1}^D \Sigma_{m=1}^M \{ \Sigma_{n=1}^N \langle \ln p(S_{d,m,n}|W_{d,m},H_{m,n}) \rangle_{q({\bf S})q({\bf H})} + \ln p(W_{d,m}) \} + const. \\ \ln q({\bf H}) &=& \langle \ln p({\bf S}|{\bf W},{\bf H}) \rangle_{q({\bf S})q({\bf W})} + \ln p({\bf H}) + const. \\ &=& \Sigma_{m=1}^M \Sigma_{n=1}^N \{ \Sigma_{d=1}^D \langle \ln p(S_{d,m,n}|W_{d,m},H_{m,n}) \rangle_{q({\bf S})q({\bf W})} + \ln p(H_{m,n}) \} + const. \\ \ln q({\bf S}) &=& \ln p({\bf X}|{\bf S}) + \langle \ln p({\bf S}|{\bf W},{\bf H}) \rangle_{q({\bf W})q({\bf H})} + const. \\ &=& \Sigma_{d=1}^D \Sigma_{n=1}^N \{ \ln p(X_{d,n}|\Sigma_{m=1}^M S_{d,m,n}) \Sigma_{m=1}^M \langle \ln p(S_{d,m,n}|W_{d,m},H_{m,n}) \rangle_{q({\bf W})q({\bf H})} \} + const. \end{eqnarray} $$${\bf W}$の近似事後分布
$\ln q({\bf W})$は$D \times M$個の分布の対数の和で表されているので、$q({\bf W})$は$D \times M$個の独立な分布の積であり、要素$W_{d,m}$に関する近似分布を求めれば導出として十分である。$W_{d,m}$に無関係な項は無視して期待値の部分を計算すると、
$$ \begin{eqnarray} \langle \ln p({\bf S}|{\bf W},{\bf H}) \rangle_{q({\bf S})q({\bf H})} &=& \langle \ln \mathrm{Poi} (S_{d,m,n}|W_{d,m} H_{m,n}) \rangle_{q({\bf S})q({\bf H})} \\ &=& \langle S_{d,m,n} \ln W_{d,m} H_{m,n} - \ln S_{d,m,n} ! - W_{d,m} H_{m,n} \rangle_{q({\bf S})q({\bf H})} \\ &=& \langle S_{d,m,n} \rangle \ln W_{d,m} - \langle H_{m,n} \rangle W_{d,m} + const. \end{eqnarray} $$また、ガンマ事前分布$p(W_{d,m})$の対数をとると
$$ \ln p(W_{d,m}|a_W,b_W) = (a_W - 1) \ln W_{d,m} - b_W W_{d,m} + const. $$となるので
$$ \begin{eqnarray} \ln q(W_{d,m}) &=& (\Sigma_{n=1}^N \langle S_{d,m,n} \rangle + a_W - 1) \ln W_{d,m} - (\Sigma_{n=1}^N \langle H_{m,n} \rangle + b_W) W_{d,m} + const. \end{eqnarray} $$となり、これはガンマ分布に対数をとったものである。よって次のように書ける。
$$ \begin{eqnarray} q(W_{d,m}) &=& \mathrm{Gam} (W_{d,m}| \hat{a}_W^{(d,m)}, \hat{b}_W^{(m)}) \\ ただし \hat{a}_W^{(d,m)} &=& \Sigma_{n=1}^N \langle S_{d,m,n} \rangle + a_W \\ \hat{b}_W^{(m)} &=& \Sigma_{n=1}^N \langle H_{m,n} \rangle + b_W \end{eqnarray} $$${\bf H}$の近似事後分布
${\bf W}$との対称性から次のようになる。
$$ \begin{eqnarray} q(H_{m,n}) &=& \mathrm{Gam} (H_{m,n}| \hat{a}_H^{(m,n)}, \hat{b}_H^{(m)}) \\ ただし \hat{a}_H^{(m,n)} &=& \Sigma_{d=1}^D \langle S_{d,m,n} \rangle + a_H \\ \hat{b}_H^{(m)} &=& \Sigma_{d=1}^D \langle W_{d,m} \rangle + b_H \end{eqnarray} $$${\bf S}$の近似事後分布
${\bf S}$の分布に対数をとったものも同様にdとnに関する和で分解できる。期待値の部分を整理すると、
$$ \begin{eqnarray} \Sigma_{m=1}^M \langle \ln p(S_{d,m,n}|W_{d,m},H_{m,n}) \rangle_{q({\bf W})q({\bf H})} &=& \Sigma_{m=1}^M \langle \ln \mathrm{Poi}(S_{d,m,n}|W_{d,m}H_{m,n}) \rangle_{q({\bf W})q({\bf H})} \\ &=& - \Sigma_{m=1}^M \ln S_{d,m,n}! + \Sigma_{m=1}^M S_{d,m,n} (\langle \ln W_{d,m} \rangle + \langle \ln H_{m,n} \rangle) + const. \end{eqnarray} $$となる。これは多項分布と同じ関数形式であり、またデルタ分布$p({\bf X}|{\bf S})$が$\Sigma_{m=1}^M S_{d,m,n} = X_{d,n}$を満たすことを要求しているので近似分布$q({\bf S})$は次のような試行回数を$X_{d,n}$としたM項分布として表現できる。
$$ \begin{eqnarray} q(S_{d,:,n}) &=& \mathrm{Mult}(S_{d,:,n}|\hat{\pi}_{d,n},X_{d,n}) \\ ただし \hat{\pi}_{d,n}^{(m)} &\propto& \exp (\langle \ln W_{d,m} \rangle + \langle \ln H_{m,n} \rangle) \\ (\mathrm{s.t.} \Sigma_{m=1}^M \hat{\pi}_{d,n}^{(m)} &=& 1) \end{eqnarray} $$期待値
各近似分布の更新で必要になる期待値は次のように解析的に計算できる。
$$ \begin{eqnarray} \langle S_{d,m,n} \rangle &=& X_{d,n} \hat{\pi}_{d,n}^{(m)} \\ \langle W_{d,m} \rangle &=& \frac{\hat{a}_{W}^{(d,m)} }{\hat{b}_{W}^{(m)} } \\ \langle H_{m,n} \rangle &=& \frac{\hat{a}_{H}^{(m,n)} }{\hat{b}_{H}^{(m)} } \\ \langle \ln W_{d,m} \rangle &=& \psi (\hat{a}_{W}^{(d,m)}) - \ln \hat{b}_{W}^{(m)} \\ \langle \ln H_{m,n} \rangle &=& \psi (\hat{a}_{H}^{(m,n)}) - \ln \hat{b}_{H}^{(m)} \end{eqnarray} $$※$\psi()$はディガンマ関数
from scipy.special import digamma
class NMF_model():
def __init__(self, X, M):
self.X = X
self.D, self.N = X.shape
self.M = M
self.E_S = np.zeros((self.D,self.M,self.N))
self.a_W = np.random.random((self.D, self.M)) + 1e-5
self.b_W = np.random.random((self.M)) + 1e-5
self.a_H = np.random.random((self.M, self.N)) + 1e-5
self.b_H = np.random.random((self.M)) + 1e-5
self.W = np.zeros((self.D,self.M))
self.H = np.zeros((self.M,self.N))
def calc_E_Sdmn(self):
for d in range(self.D):
for n in range(self.N):
E_lnWd = digamma(self.a_W[d,:]) - np.log(self.b_W)
E_lnHn = digamma(self.a_H[:,n]) - np.log(self.b_H)
pidn = np.exp(E_lnWd + E_lnHn)
pidn = pidn / np.sum(pidn)
E_Sdn = self.X[d,n] * pidn
self.E_S[d,:,n] = E_Sdn
def update_W(self):
self.a_W += np.sum(self.E_S, axis=2)
self.b_W += np.sum(self.a_H/self.b_H[:,np.newaxis], axis=1)
def update_H(self):
self.a_H += np.sum(self.E_S, axis=0)
self.b_H += np.sum(self.a_W/self.b_W[np.newaxis,:], axis=0)
def fit(self, n_iter):
for i in range(n_iter):
self.calc_E_Sdmn()
self.update_W()
self.update_H()
def get_W(self):
for d in range(self.D):
for m in range(self.M):
self.W[d,m] = np.random.gamma(self.a_W[d,m], 1/self.b_W[m])
return self.W
def get_H(self):
for m in range(self.M):
for n in range(self.N):
self.H[m,n] = np.random.gamma(self.a_H[m,n], 1/self.b_H[m])
return self.H
# numpy配列に
X_for_bayes = np.array(X)
# M=1
bayes_NMF1 = NMF_model(X_for_bayes, 1)
bayes_NMF1.fit(1000)
W1 = bayes_NMF1.get_W()
H1 = bayes_NMF1.get_H()
# 特徴量1つなので、合計種族値にほぼ比例する
plt.scatter(pokemon.iloc[:,3:9].sum(axis=1),W1)
pokemon['W1_bayes'] = W1
pokemon.loc[np.argsort(- W1,axis=0).flatten()[:10],:]
# ヒートマップで確認
fig, axes = plt.subplots(1,2, figsize=(10,4))
sns.heatmap(X, ax=axes[0],vmin=vmin,vmax=vmax)
sns.heatmap(np.dot(W1,H1), ax=axes[1], vmin=vmin, vmax=vmax)
# M=2
bayes_NMF2 = NMF_model(X_for_bayes, 2)
bayes_NMF2.fit(3000)
W2 = bayes_NMF2.get_W()
H2 = bayes_NMF2.get_H()
pokemon['W2_1_bayes'] = W2[:,0]
pokemon.loc[np.argsort(- W2[:,0])[:10],:]
pokemon['W2_2_bayes'] = W2[:,1]
pokemon.loc[np.argsort(- W2[:,1])[:10],:]
# つよいポケモン よわいポケモン
# そんなの ひとの かって
# ほんとうに つよい トレーナーなら
# すきなポケモンで かてるように がんばるべき
pokemon.loc[np.argsort(np.sqrt(W2[:,0]**2 + W2[:,1]**2))[:10],:]
# ヒートマップで確認
fig, axes = plt.subplots(1,3, figsize=(15,4))
sns.heatmap(X, ax=axes[0],vmin=vmin,vmax=vmax)
sns.heatmap(np.dot(W1,H1), ax=axes[1],vmin=vmin,vmax=vmax)
sns.heatmap(np.dot(W2,H2), ax=axes[2],vmin=vmin,vmax=vmax)
sklearnと比べると性能がイマイチですが、実装はできましたね。