はじめに
『トピックモデルによる統計的潜在意味解析』の学習時のメモです。本と併せて読んでください。
この記事では、3.3.8項のLDAの周辺化変分ベイズ法について書いています。図3.5の疑似コードを基にR言語で実装していきます(未完)。
プログラムからアルゴリズムの理解を目指します。できるだけ基本的な(直感的に理解できる)関数を用いて組んでいきます。
【数理編】
【前節の内容】
【他の節一覧】
【この節の内容】
3.3.8 LDAの周辺化変分ベイズ法
図3.5の疑似コードを参考にLDAに対する周辺化変分ベイズ法を実装していきます(未完)。
# 利用パッケージ library(tidyverse)
・簡易文書データ(クリックで展開)
## 簡易文書データ # 文書数 M <- 10 # 語彙数 V <- 20 # 文書ごとの各語彙数 n_dv <- matrix(sample(0:10, M * V, replace = TRUE), M, V)
・パラメータの設定
始めにパラメータ類の初期値を設定していきます。
# イタレーション数を指定 Iter <- 50 # トピック数を指定 K <- 5
試行回数をIter、トピックの数をKとします。
# 事前分布のパラメータを指定 alpha_k <- rep(2, K) beta_v <- rep(2, V)
トピック分布のパラメータ$\boldsymbol{\alpha} = (\alpha_1, \cdots, \alpha_K)$、単語分布のパラメータ$\boldsymbol{\beta} = (\beta_1, \cdots, \beta_V)$の値をそれぞれ指定します。どちらもディリクレ分布のパラメータなので、$\alpha_k > 0,\beta_v > 0$である必要があります。
この例では値を一様に設定します。
# 潜在トピック集合の分布 z_dv_k <- array(0, dim = c(M, V, K)) for(d in 1:M) { for(v in 1:V) { if(n_dv[d, v] > 0) { # ランダムに値を生成 tmp_q_z <- sample(seq(0, 1, by = 0.01), size = K, replace = TRUE) # 正規化 z_dv_k[d, v, ] <- tmp_q_z / sum(tmp_q_z) } } }
各単語に潜在トピックを割り当てる確率$q(\boldsymbol{z})$の初期値をランダムに設定します。
$q(z_{d,i} = k)$の計算式(3.130)より同一単語の確率は同じになることから、ここではz_dv_kとして、文書・(単語$i$ではなく)語彙・トピックの3次元配列として扱うことにします。
## 割り当てられたトピックに関する統計量を計算 tmp_p_n <- array(0, dim = c(M, V, K)) tmp_q_n <- array(0, dim = c(M, V, K)) for(k in 1:K) { tmp_p_n[, , k] <- z_dv_k[, , k] * n_dv tmp_q_n[, , k] <- z_dv_k[, , k] * (1 - z_dv_k[, , k]) * n_dv } # 各文書において各トピックが割り当てられる単語数の期待値と分散 E_n_dk <- apply(tmp_p_n, c(1, 3), sum) V_n_dk <- apply(tmp_q_n, c(1, 3), sum) # 全文書の各語彙において各トピックが割り当てられた語数の期待値と分散 E_n_kv <- apply(tmp_p_n, c(3, 2), sum) V_n_kv <- apply(tmp_q_n, c(3, 2), sum)
$n_{d,k}^{\backslash d,i},\ n_{k,v}^{\backslash d,i}$の期待値と分散を計算式
に従い計算します。
ただしこの段階では、$d, i$要素を含む$n_{d,k},\ n_{k,v}$の期待値と分散としておきます。またこの例では$q(z_{d,i'} = k)$を(重複を考慮しない)語彙として扱っているため、各語彙に対して出現回数n_dvを掛けることで$\sum_{i=1}^{n_d}$の計算を行います。
各語彙の推定の段階で、$d, i$要素に関する値を除きます。
・周辺化変分ベイズ法
続いて推論部分を実装していきます。
・コード全体(クリックで展開)
# 推移の確認用 trace_alpha <- matrix(0, nrow = K, ncol = Iter + 1) trace_beta <- matrix(0, nrow = V, ncol = Iter + 1) # 初期値を代入 trace_alpha[, 1] <- alpha_k trace_beta[, 1] <- beta_v for(I in 1:Iter) { # 動作確認用 start_time <- Sys.time() # 各種統計量を初期化 new_z_dv_k <- array(0, dim = c(M, V, K)) new_E_n_dk <- new_V_n_dk <- matrix(0, M, K) new_E_n_kv <- new_V_n_kv <- matrix(0, K, V) for(d in 1:M) { ## (各文書) for(v in 1:V) { ## (各語彙) if(n_dv[d, v] > 0) { # 各種統計量を移す E_n_dk.di <- E_n_dk V_n_dk.di <- V_n_dk E_n_kv.di <- E_n_kv V_n_kv.di <- V_n_kv # (前回の)d,i要素を除く E_n_dk.di[d, ] <- E_n_dk.di[d, ] - z_dv_k[d, v, ] V_n_dk.di[d, ] <- V_n_dk.di[d, ] - z_dv_k[d, v, ] * (1 - z_dv_k[d, v, ]) E_n_kv.di[, v] <- E_n_kv.di[, v] - z_dv_k[d, v, ] V_n_kv.di[, v] <- V_n_kv.di[, v] - z_dv_k[d, v, ] * (1 - z_dv_k[d, v, ]) # 潜在トピック集合の分布を計算:式(3.130) term1 <- (E_n_kv.di[, v] + beta_v[v]) / apply(t(E_n_kv.di) + beta_v, 2, sum) * (E_n_dk.di[d, ] + alpha_k) term21 <- V_n_kv.di[, v] / (2 * (E_n_kv.di[, v] + beta_v[v])^2) ## (CVB2 ver) term22 <- V_n_dk.di[d, ] / (2 * (E_n_dk.di[d, ] + alpha_k)^2) ## (CVB2 ver) term3 <- apply(V_n_kv.di, 1, sum) / (2 * apply(t(E_n_kv.di) + beta_v, 2, sum)^2) ## (CVB2 ver) tmp_q_z <- term1 * exp(- term21 - term22) * exp(term3) ## (CVB2 ver) # 値を正規化 new_z_dv_k[d, v, ] <- tmp_q_z / sum(tmp_q_z) ## (CVB2 ver) #new_z_dv_k[d, v, ] <- term1 / sum(term1) ## (CVB0 ver) # 各種統計量に(今回の)d,i要素を加える new_E_n_dk[d, ] <- new_E_n_dk[d, ] + new_z_dv_k[d, v, ] * n_dv[d, v] new_V_n_dk[d, ] <- new_V_n_dk[d, ] + new_z_dv_k[d, v, ] * (1 - new_z_dv_k[d, v, ]) * n_dv[d, v] new_E_n_kv[, v] <- new_E_n_kv[, v] + new_z_dv_k[d, v, ] * n_dv[d, v] new_V_n_kv[, v] <- new_V_n_kv[, v] + new_z_dv_k[d, v, ] * (1 - new_z_dv_k[d, v, ]) * n_dv[d, v] } } ## (/各語彙) } ## (/各文書) # 潜在トピック集合の分布を更新 z_dv_k <- new_z_dv_k # 各種統計量を更新 E_n_dk <- new_E_n_dk E_n_kv <- new_E_n_kv # トピック分布のパラメータを更新:式(3.191) alpha_numer <- apply(digamma(t(E_n_dk) + alpha_k) - digamma(alpha_k), 1, sum) alpha_denom <- sum(digamma(apply(t(E_n_dk) + alpha_k, 2, sum)) - digamma(sum(alpha_k))) alpha_k <- alpha_k * alpha_numer / alpha_denom # 単語分布のパラメータを更新:式(3.192) beta_numer <- apply(digamma(t(E_n_kv) + beta_v) - digamma(beta_v), 1, sum) beta_denom <- sum(digamma(apply(t(E_n_kv) + beta_v, 2, sum)) - digamma(sum(beta_v))) beta_v <- beta_v * beta_numer / beta_denom # 推移の確認用 trace_alpha[, I + 1] <- alpha_k trace_beta[, I + 1] <- beta_v # 動作確認 #print(paste0(I, "th try...", round(Sys.time() - start_time, 3))) }
$q(z_{d,i} = k)$の計算前に、各統計量からz_dv_k(前回の値)の$d, i$要素の値を引き、それぞれ変数名の末尾を.diとした別の変数としています。
また各計算結果も、変数名の頭にnew_と付けた別の変数に保存しておきます。それぞれ全ての要素を計算してから、元の変数の値を更新します。
・$q(\boldsymbol{z})$の更新
# 潜在トピック集合の分布を計算:式(3.130) term1 <- (E_n_kv.di[, v] + beta_v[v]) / apply(t(E_n_kv.di) + beta_v, 2, sum) * (E_n_dk.di[d, ] + alpha_k) term21 <- V_n_kv.di[, v] / (2 * (E_n_kv.di[, v] + beta_v[v])^2) ## (CVB2 ver) term22 <- V_n_dk.di[d, ] / (2 * (E_n_dk.di[d, ] + alpha_k)^2) ## (CVB2 ver) term3 <- apply(V_n_kv.di, 1, sum) / (2 * apply(t(E_n_kv.di) + beta_v, 2, sum)^2) ## (CVB2 ver) tmp_q_z <- term1 * exp(- term21 - term22) * exp(term3) ## (CVB2 ver) # 値を正規化 new_z_dv_k[d, v, ] <- tmp_q_z / sum(tmp_q_z) #new_z_dv_k[d, v, ] <- term1 / sum(term1) ## (CVB0 ver)
語彙ごとに、式(3.130)の計算を行い値を更新します。
またこの計算結果を総和が1になるように正規化します。
テイラー展開を0次近似でする場合(CVB0)は、(CVB2 ver)とコメントしている計算をコメントアウトして、(CVB0 ver)のコメントアウトを外してください。
・$q(\boldsymbol{\theta})$の更新
# 各種統計量に(今回の)d,i要素を加える new_E_n_dk[d, ] <- new_z_dv_k[d, v, ] * n_dv[d, v] new_V_n_dk[d, ] <- new_z_dv_k[d, v, ] * (1 - new_z_dv_k[d, v, ]) * n_dv[d, v]
それぞれ計算式(3.124)、(3.125)に従い、値を更新します。
ただし$\sum_{i=1}^{n_d}$の計算について、new_z_dv_k[d, v, ]をループ処理の中でV回繰り返し足していくこと、また各語彙に対して出現回数(n_dv[d, ])を掛けることで行います。
これを文書1から文書$M$まで繰り返し更新し終えてから、最後にハイパーパラメータalpha_kを更新します。
# トピック分布のパラメータを更新:式(3.191) alpha_numer <- apply(digamma(t(E_n_dk) + alpha_k) - digamma(alpha_k), 1, sum) alpha_denom <- sum(digamma(apply(t(E_n_dk) + alpha_k, 2, sum)) - digamma(sum(alpha_k))) alpha_k <- alpha_k * alpha_numer / alpha_denom
$q(\boldsymbol{\theta}_d)$のパラメータの更新式(3.191)
の計算を行い、値を更新します。
・$q(\boldsymbol{\phi})$の更新
# 各種統計量から(今回の)d,i要素を加える new_E_n_kv[, v] <- new_E_n_kv[, v] + new_z_dv_k[d, v, ] * n_dv[d, v] new_V_n_kv[, v] <- new_V_n_kv[, v] + new_z_dv_k[d, v, ] * (1 - new_z_dv_k[d, v, ]) * n_dv[d, v]
それぞれの定義式(3.126)、(3.128)に従い、値を更新します。
これまでと同様に、n_dv[d, v]を掛けることで$\sum_{i=1}^{n_d}$の計算となります。またループ処理の中でnew_z_dv_k[d, v, ]をM回繰り返し足すことで、$\sum_{d=1}^M$の計算となります。
これも文書1から文書$M$まで繰り返し更新し終えてから、最後にハイパーパラメータbeta_vを更新します。
# 単語分布のパラメータを更新:式(3.192) beta_numer <- apply(digamma(t(E_n_kv) + beta_v) - digamma(beta_v), 1, sum) beta_denom <- sum(digamma(apply(t(E_n_kv) + beta_v, 2, sum)) - digamma(sum(beta_v))) beta_v <- beta_v * beta_numer / beta_denom
$q(\boldsymbol{\phi}_k)$のパラメータの更新式(3.192)
の計算を行い、値を更新します。
・推定結果の確認
ggplot2パッケージを利用して、推定結果を可視化しましょう。
ディリクレ分布の期待値(2.10)の計算を行いパラメータの期待値を作図します。
・トピック分布(クリックで展開)
## トピック分布(期待値) # thetaの期待値を計算:式(2.10) theta_k <- alpha_k / sum(alpha_k) # 作図用のデータフレームを作成 theta_df <- tibble( prob = theta_k, topic = as.factor(1:K) ) # 作図 ggplot(theta_df, aes(x = topic, y = prob, fill = topic)) + geom_bar(stat = "identity", position = "dodge") + # 棒グラフ labs(title = "Collapsed Variational Bayes Method for LDA", subtitle = expression(theta)) # ラベル
・単語分布(クリックで展開)
## 単語分布(期待値) # phiの期待値を計算:式(2.10) phi_v <- beta_v / sum(beta_v) # 作図用のデータフレームを作成 phi_df <- tibble( prob = phi_v, word = as.factor(1:V) ) # 作図 ggplot(phi_df, aes(x = word, y = prob, fill = word, color = word)) + geom_bar(stat = "identity", position = "dodge") + # 棒グラフ scale_x_discrete(breaks = seq(0, V, by = 10)) + # x軸目盛 theme(legend.position = "none") + # 凡例 labs(title = "Collapsed Variational Bayes Method for LDA", subtitle = expression(phi)) # ラベル
・更新値の推移の確認
・トピック分布のパラメータ(クリックで展開)
## トピック分布のパラメータ # 作図用のデータフレームを作成 trace_alpha_df <- data.frame() for(I in 1:(Iter + 1)) { # データフレームを作成 tmp_trace_alpha <- data.frame( value = trace_alpha[, I], topic = as.factor(1:K), iteration = I - 1 ) # 結合 trace_alpha_df <- rbind(trace_alpha_df, tmp_trace_alpha) } # 作図 ggplot(trace_alpha_df, aes(x = iteration, y = value, color = topic)) + geom_line() + # 折れ線グラフ labs(title = "Collapsed Variational Bayes Method for LDA", subtitle = expression(alpha)) # ラベル
・単語分布のパラメータ(クリックで展開)
## 単語分布のパラメータ # 作図用のデータフレームを作成 trace_beta_df <- data.frame() for(I in 1:(Iter + 1)) { # データフレームを作成 tmp_trace_beta <- data.frame( value = trace_beta[, I], word = as.factor(1:V), iteration = I - 1 ) # 結合 trace_beta_df <- rbind(trace_beta_df, tmp_trace_beta) } # 作図 ggplot(trace_beta_df, aes(x = iteration, y = value, color = word)) + geom_line(alpha = 0.5) + # 折れ線グラフ theme(legend.position = "none") + # 凡例 labs(title = "Collapsed Variational Bayes Method for LDA", subtitle = expression(beta)) # ラベル
・gif画像を作成
更新ごとの分布の変化をgif画像でみるためのコードです。gganimateパッケージを利用してgif画像を作成します。作図には上で作成したデータフレームを使います。
・コード(クリックで展開)
# 利用パッケージ library(gganimate) ## トピック分布のパラメータ # 作図 graph_alpha <- ggplot(trace_alpha_df, aes(topic, value, fill = topic)) + geom_bar(stat = "identity", position = "dodge") + # 棒グラフ transition_manual(iteration) + # フレーム labs(title = "Collapsed Variational Bayes Method for LDA", subtitle = "I={current_frame}") # ラベル # gif画像を作成 animate(graph_alpha, nframes = (Iter + 1), fps = 10) ## 単語分布のパラメータ # 作図 graph_beta <- ggplot(trace_beta_df, aes(word, value, fill = word, color = word)) + geom_bar(stat = "identity", position = "dodge") + # 棒グラフ scale_x_discrete(breaks = seq(0, V, by = 10)) + # x軸目盛 theme(legend.position = "none") + # 凡例 transition_manual(iteration) + # フレーム labs(title = "Collapsed Variational Bayes Method for LDA", subtitle = "I={current_frame}") # ラベル # gif画像を作成 animate(graph_beta, nframes = (Iter + 1), fps = 10)
transition_manual()に時系列値の列を指定します。
参考文献
- 佐藤一誠『トピックモデルによる統計的潜在意味解析』(自然言語処理シリーズ 8)奥村学監修,コロナ社,2015年.
おわりに
こんな感じで組めませんでした。ちなみに周辺化ギブスサンプリングの方も組めてません。
1か月程粘りましたがこのままではムリそうだったので、区切りとして上げることにしました。誰か教え……実装例が詳しく書いてある本が欲しいのですが、今はゼロつくをやらなきゃだったりで必要不急な外出なのでどーすかなぁと思っております。
まぁそんなこんなで、あと2つ書いたらゼロつくネタが続く予定です(本当はその前に白トピ本(前半)を終わらせたかった…)。Python本も欲しい！
次の分はもう書き終わってます。そっちは実装できたので安心です。また明日！
【次節の内容】
