はじめに

pythonでトピックモデルの最尤推定を実装しました．
教科書としてトピックモデルを使いました．

本記事の構成

はじめに
トピックモデル
確率的潜在意味解析
実装
結果
おわりに

トピックモデル

トピックモデルは，1つの文書が複数のトピックを持つと仮定します．

例
「医療に関する法案の国会審議」に関する記事 => 「医療」と「政治」の2つのトピック
「オリンピックの経済効果」に関する記事 => 「スポーツ」と「経済」の2つのトピック

1つの文書が複数のトピックを持つという仮定をうまく表現するため，
文書ごとにトピック分布 $\boldsymbol \theta_d = (\theta_{d1}, \cdots, \theta_{dK})$ を設定します．
$\theta_{dk} = p(k \ | \ \boldsymbol \theta_d)$ は，文書 $d$ の単語にトピック $k$ が割り当てられる確率です．
トピック分布 $\boldsymbol \theta_d$ に従って文書 $d$ の各単語にトピック $z_{dn}$ が割り当てられます．
割り当てられたトピックの単語分布 $\boldsymbol \phi_{z_{dn}}$ に従って単語が生成されます．
$\boldsymbol \phi_k = (\phi_{k1}, \cdots, \phi_{kV})$ は，トピック $k$ で語彙 $v$ が生成される確率を表しています．

具体的な文書集合の生成過程を下図に示します．

単語ごとにトピックが割り当てられ，そのトピックに応じた単語分布から単語が生成されます．
図の中段では，「スポーツ」と「経済」のトピックが割り当てられる確率が高いため，
この2つのトピックに関連した単語が多く観測されます．
経済トピックの単語分布を見ると，生成される確率の高い順に
「経済」「株価」「景気」「円高」と並んでいます．

確率的潜在意味解析

トピックモデルを最尤推定する手法は，確率的潜在意味解析(PLSA)と呼ばれます．
ちなみにトピックモデルのベイズ推定は，潜在ディリクレ配分モデル(LDA)という手法です．
※ LDAは全く理解していないので解説できません．

PLSAでは，EMアルゴリズムを使ってトピックモデルの最尤推定をします．
以下の対数尤度を最大にするパラメータ $\boldsymbol \Theta = (\boldsymbol \theta_1, \cdots, \boldsymbol \theta_D)$ および $ \boldsymbol \Phi = (\boldsymbol \phi_1, \cdots, \boldsymbol \phi_K)$ を求めます．

\begin{align}
L &= \sum^D_{d = 1} \sum^{N_d}_{n = 1} \log \sum^K_{k = 1} \theta_{dk}\phi_{kw_{dn}} \\
&= \sum^D_{d = 1} \sum^{N_d}_{n = 1} \log \sum^K_{k = 1} q_{dnk} \frac{\theta_{dk}\phi_{kw_{dn}}}{q_{dnk}} \\
&\geq \sum^D_{d = 1} \sum^{N_d}_{n = 1} \sum^K_{k = 1} q_{dnk} \log \frac{\theta_{dk}\phi_{kw_{dn}}}{q_{dnk}} \\
&= \sum^D_{d = 1} \sum^{N_d}_{n = 1} \sum^K_{k = 1} q_{dnk} \log \theta_{dk}\phi_{kw_{dn}} - \sum^D_{d = 1} \sum^{N_d}_{n = 1} \sum^K_{k = 1} q_{dnk} \log q_{dnk} \equiv F \tag{1}
\end{align}

下限 $F$ を最大化するためにEMアルゴリズムを使います．
以下 E step と M step に分けて式展開していきます．

※ ここで，使用される文字について一度整理しておきます．

$D$：文書数
$N_d$：文書 $d$ に含まれる単語数（文書長）
$V$：全文書の中で現れる単語の種類数（語彙数）
$\boldsymbol W$：文書集合
$\boldsymbol w_d$：文書 $d$
$w_{dn}$：文書 $d$ の $n$ 番目の単語
$K$：トピック数
$\theta_{dk}$：文書 $d$ でトピック $k$ が割り当てられる確率
$\phi_{kv}$：トピック $k$ で語彙 $v$ が生成される確率
$q_{dnk}$：文書 $d$ において $n$ 番目の単語にトピック $k$ が割り当てられる確率（負担率）

E step

ラグランジュの未定乗数法を用いるために $F_q$ を以下のように定義します．

F_q \equiv \sum^D_{d = 1} \sum^{N_d}_{n = 1} \sum^K_{k = 1} q_{dnk} \log \theta_{dk}\phi_{kw_{dn}} - \sum^D_{d = 1} \sum^{N_d}_{n = 1} \sum^K_{k = 1} q_{dnk} \log q_{dnk} + \lambda \bigl(\sum^K_{k = 1} q_{dnk} - 1 \bigr) \tag{2}

$F_q$ を $q_{dnk}$ で偏微分します．

\frac{\partial F_q}{\partial q_{dnk}} = \log \theta_{dk}\phi_{kw_{dn}} - \log q_{dnk} - 1 + \lambda = 0 \tag{3}

$\frac{\partial F_q}{\partial q_{dnk}} = 0$ より，

q_{dnk} = \exp(\lambda - 1) \theta_{dk}\phi_{kw_{dn}} \tag{4}

$\sum^K_{k = 1} q_{dnk} = 1$ より，

\begin{align}
& \sum^K_{k = 1} \exp(\lambda - 1) \theta_{dk}\phi_{kw_{dn}} = 1 \\
& \exp(\lambda - 1) = \frac{1}{\displaystyle \sum^K_{k = 1} \theta_{dk}\phi_{kw_{dn}}} \tag{5}
\end{align}

式(5)を式(4)に代入して $q_{dnk}$ を得ます．

q_{dnk} = \frac{\theta_{dk}\phi_{kw_{dn}}}{\displaystyle \sum^K_{k' = 1} \theta_{dk'}\phi_{k'w_{dn}}} \tag{6}

$q_{dnk}$ は負担率と呼ばれ，文書 $d$ において $n$ 番目の単語にトピック $k$ が割り当てられる確率を表します．

M step

M stepでは，パラメータ $\boldsymbol \Theta = (\boldsymbol \theta_1, \cdots, \theta_D)$ と $ \boldsymbol \Phi = (\phi_1, \cdots, \phi_K)$ を分けて考えます．
どちらも E stepと同様にラグランジュの未定乗数法を用います．

トピック分布

ラグランジュの未定乗数法を用いるために $F_\theta$ を以下のように定義します．

F_\theta \equiv \sum^D_{d = 1} \sum^{N_d}_{n = 1} \sum^K_{k = 1} q_{dnk} \log \theta_{dk}\phi_{kw_{dn}} - \sum^D_{d = 1} \sum^{N_d}_{n = 1} \sum^K_{k = 1} q_{dnk} \log q_{dnk} + \lambda \bigl(\sum^K_{k = 1} \theta_{dk} - 1 \bigr) \tag{7}

$F_\theta$ を $\theta_{dk}$ で偏微分します．

\frac{\partial F_\theta}{\partial \theta_{dk}} = \sum^{N_d}_{n = 1} \frac{q_{dnk}}{\theta_{dk}} + \lambda \tag{8}

$\frac{\partial F_\theta}{\partial \theta_{dk}} = 0$ より，

\theta_{dk} = - \sum^{N_d}_{n = 1} \frac{q_{dnk}}{\lambda} \tag{9}

$\sum^K_{k = 1} \theta_{dk} = 1$ より，

\begin{align}
& - \sum^K_{k = 1} \sum^{N_d}_{n = 1} \frac{q_{dnk}}{\lambda} = 1 \\
& \lambda = - \sum^K_{k = 1} \sum^{N_d}_{n = 1} q_{dnk} \tag{10}
\end{align}

式(10)を式(9)に代入して $\theta_{dk}$ を得ます．

\theta_{dk} = \frac{\displaystyle \sum^{N_d}_{n = 1} q_{dnk}}{\displaystyle \sum^K_{k' = 1} \sum^{N_d}_{n = 1} q_{dnk'}} \tag{11}

単語分布

ラグランジュの未定乗数法を用いるために $F_\phi$ を以下のように定義します．

F_\phi \equiv \sum^D_{d = 1} \sum^{N_d}_{n = 1} \sum^K_{k = 1} q_{dnk} \log \theta_{dk}\phi_{kw_{dn}} - \sum^D_{d = 1} \sum^{N_d}_{n = 1} \sum^K_{k = 1} q_{dnk} \log q_{dnk} + \lambda \bigl(\sum^V_{v = 1} \phi_{kv} - 1 \bigr) \tag{12}

$F_\phi$ を $\phi_{kv}$ で偏微分します．

\frac{\partial F_\phi}{\partial \phi_{kv}} = \sum^D_{d = 1} \sum^{}_{n:w_{dn} = v} \frac{q_{dnk}}{\phi_{kv}} + \lambda \tag{13}

$\frac{\partial F_\phi}{\partial \phi_{kv}} = 0$ より，

\phi_{kv} = - \sum^D_{d = 1} \sum^{}_{n:w_{dn} = v} \frac{q_{dnk}}{\lambda} \tag{14}

$\sum^V_{v = 1} \phi_{kv} = 1$ より，

\begin{align}
& - \sum^V_{v = 1} \sum^D_{d = 1} \sum^{}_{n:w_{dn} = v} \frac{q_{dnk}}{\lambda} = 1 \\
& \lambda = - \sum^V_{v = 1} \sum^D_{d = 1} \sum^{}_{n:w_{dn} = v} q_{dnk} \tag{15}
\end{align}

式(15)を式(14)に代入して $\phi_{kv}$ を得ます．

\phi_{kv} = \frac{\displaystyle \sum^D_{d = 1} \sum^{}_{n:w_{dn} = v} q_{dnk}}{\displaystyle \sum^V_{v' = 1} \sum^D_{d = 1} \sum^{}_{n:w_{dn} = v'} q_{dnk}}

実装

pythonでトピックモデルの最尤推定のトイプログラムを実装しました．

import numpy as np

def normalize(ndarray, axis):
	return ndarray / ndarray.sum(axis = axis, keepdims = True)
	
def normalized_random_array(d0, d1):
	ndarray = np.random.rand(d0, d1)
	return normalize(ndarray, axis = 1)

if __name__ == "__main__":
	
	# initialize parameters
	D, K, V = 1000, 3, 6
	theta = normalized_random_array(D, K)
	phi = normalized_random_array(K, V)
	theta_est = normalized_random_array(D, K)
	phi_est = normalized_random_array(K, V)

	# for generate documents
	_theta = np.array([theta[:, :k+1].sum(axis = 1) for k in range(K)]).T
	_phi = np.array([phi[:, :v+1].sum(axis = 1) for v in range(V)]).T

	# generate documents
	W, Z = [], []
	N = np.random.randint(100, 300, D)
	for (d, N_d) in enumerate(N):
		Z.append((np.random.rand(N_d, 1) < _theta[d, :]).argmax(axis = 1))
		W.append((np.random.rand(N_d, 1) < _phi[Z[-1], :]).argmax(axis = 1))
	
	# estimate parameters
	q = np.zeros((D, np.max(N), K))
	T = 300
	likelihood = np.zeros(T)
	for t in range(T):
		print(t)

		# E step
		for (d, N_d) in enumerate(N):
			q[d, :N_d, :] = normalize(theta_est[d, :] * phi_est[:, W[d]].T, axis = 1)

		# M step
		theta_est[:, :] = normalize(q.sum(axis = 1), axis = 1)
		q_sum = np.zeros((K, V))
		for (d, W_d) in enumerate(W):
			v, index, count = np.unique(W_d, return_index= True, return_counts = True)
			q_sum[:, v[:]] += count[:] * q[d, index[:], :].T
		phi_est[:, :] = normalize(q_sum, axis = 1)

		# likelihood
		for (d, W_d) in enumerate(W):
			likelihood[t] += np.log(theta_est[d, :].dot(phi_est[:, W_d])).sum()

	# perplexity
	perplexity = np.exp(-likelihood[:] / N.sum())

結果

実装したコードを動かした結果を示します．
繰り返しになりますが，各パラメータは以下になります．

文書数 $D = 1000$
トピック数 $K = 3$
語彙数 $V = 6$
学習回数 $T = 300$

真の単語分布と推定した単語分布を以下に示します．

phi
[[ 0.06837655  0.2294022   0.29048815  0.01610669  0.31437584  0.08125057]
 [ 0.19902324  0.19283392  0.0251427   0.23391767  0.10622823  0.24285424]
 [ 0.04165762  0.30678289  0.2579833   0.00882557  0.25580282  0.12894781]]

phi_est
[[ 0.04572686  0.20003413  0.33589344  0.00282668  0.36368718  0.0518317 ]
 [ 0.20342652  0.16884355  0.04403442  0.23473554  0.12226066  0.22669932]
 [ 0.05439782  0.37614278  0.18788228  0.01364525  0.18248741  0.18544446]]

対数尤度のグラフを示します．
対数尤度は以下の式で求めています．

\begin{align}
likelihood &= \sum^D_{d = 1} \log p(\boldsymbol w_d \ | \ M) \\
&= \sum^D_{d = 1} \log \prod^{N_d}_{n = 1} \sum^K_{k = 1} \theta_{dk} \phi_{kw_{dn}}
\end{align}

パープレキシティのグラフを示します．
※ 評価データに学習データを使っているのであまり意味はないですが．．．
パープレキシティは以下の式で求めています．

perplexity = \exp\Biggl(- \frac{\displaystyle \sum^D_{d = 1} \log p(\boldsymbol w_d \ | \ M)}{\displaystyle \sum^D_{d = 1} N_d}\Biggr)

おわりに

トピックモデルの最尤推定を実装できました．
次は LDA に挑戦しようと思います．