はじめに
機械学習プロフェッショナルシリーズの『トピックモデル』の勉強時に自分の理解の助けになったことや勉強会資料のまとめです。トピックモデルの各種アルゴリズムを「数式」と「プログラム」から理解することを目指します。
この記事は、3.3節「EMアルゴリズム」の内容です。
EMアルゴリズを用いた最尤推定による混合ユニグラムモデルにおけるパラメータ推定を実装します。
【数式編】
【他の節一覧】
【この節の内容】
3.3 EMアルゴリズム
アルゴリズム3.1を参考にして、混合ユニグラムモデル(混合カテゴリ分布)に対するEMアルゴリズムを用いた最尤推定を実装する。
この記事では、3重のループを使って実装する。「【R】3.3：混合ユニグラムモデルの最尤推定(EMアルゴリズム)の実装：ループなし版【『トピックモデル』の勉強ノート】 - からっぽのしょこ」では、行列計算を使ってループを使わずに実装する。簡潔に計算(処理)したい場合は「ループなし版」を参照のこと。
利用するパッケージを読み込む。
# 利用パッケージ library(tidyverse) library(gganimate) library(gtools)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はない。ただし、作図コードに関してはパッケージ名を省略するので、ggplot2を読み込む必要がある。
また、magrittrパッケージのパイプ演算子%>%ではなく、ネイティブパイプ演算子|>を使う。%>%に置き換えても処理できるが、その場合はmagrittrを読み込む必要がある。
文書データの簡易生成
混合ユニグラムモデルでは、トピック分布(カテゴリ分布)に従い文書ごとにトピックが割り当てられ、さらにトピックに応じた単語分布(カテゴリ分布)に従い単語が出現していると仮定する。詳しくは「3.1-2：混合ユニグラムモデル【『トピックモデル』の勉強ノート】 - からっぽのしょこ」を参照のこと。
まずは、生成モデル(真のトピック分布と単語分布)を設定して、簡易的に文書データ(bag-of-words)を生成する。
真の分布の設定
文書データに関する値を設定する。
# 文書数を指定 D <- 20 # 語彙数を指定 V <- 26 # トピック数を指定 K <- 4
生成する文書数$D$、語彙(単語の種類)数$V$、トピック数$K$を指定する。
真のトピック分布を作成する。
・コード(クリックで展開)
# 真のトピック分布を生成 theta_true_k <- MCMCpack::rdirichlet(n = 1, alpha = rep(1, times = K)) |> as.vector() theta_true_k
## [1] 0.2120760 0.2882704 0.1771578 0.3224958
ディリクレ分布からカテゴリ分布(のパラメータ)を生成して、真のトピック分布$\boldsymbol{\theta}^{(\mathrm{true})} = (\theta_1^{(\mathrm{true})}, \cdots, \theta_K^{(\mathrm{true})})$、$0 \leq \theta_k^{(\mathrm{true})} \leq 1$、$\sum_{k=1}^K \theta_k^{(\mathrm{true})} = 1$とする。または、値を直接指定する。
ディリクレ分布の乱数は、MCMCpackパッケージのrdirichlet()で生成できる。各トピックの割り当て確率を生成するため、サンプルサイズの引数nに1、パラメータの引数alphaにK個の正の値(ハイパーパラメータ)を指定する。この例では、一様な確率で生成するため、全ての値を1とする。
作図用に、生成した値をデータフレームに格納する。
# 作図用のデータフレームを作成 true_topic_df <- tibble::tibble( topic = factor(1:K), # トピック番号 probability = theta_true_k # 割り当て確率 ) true_topic_df
## # A tibble: 4 × 2 ## topic probability ## <fct> <dbl> ## 1 1 0.212 ## 2 2 0.288 ## 3 3 0.177 ## 4 4 0.322
因子型の1からKの整数をトピック番号列として、対応する確率とあわせて格納する。
真のトピック分布をグラフで確認する。
# 真のトピック分布を作図 ggplot() + geom_bar(data = true_topic_df, mapping = aes(x = topic, y = probability, fill = topic), stat = "identity", show.legend = FALSE) + # トピック分布 labs(title = "Topic Distribution", subtitle = "truth", x = "k", y = expression(theta[k]))
x軸をトピック番号、y軸を各トピックとなる確率として棒グラフを描画する。
同様に、真の単語分布を作成する。
・コード(クリックで展開)
# 真の単語分布を生成 phi_true_kv <- MCMCpack::rdirichlet(n = K, alpha = rep(1, times = V)) phi_true_kv[, 1:5]
## [,1] [,2] [,3] [,4] [,5] ## [1,] 0.0568399975 0.07531482 0.01279492 0.07577412 0.016774667 ## [2,] 0.0006434284 0.04475967 0.02762827 0.01839845 0.001307481 ## [3,] 0.0012722861 0.07577773 0.03325202 0.04024892 0.007597065 ## [4,] 0.0492145704 0.04913936 0.02513589 0.00503566 0.012337878
ディリクレ分布からカテゴリ分布(のパラメータ)を生成して、K個の真の単語分布$\boldsymbol{\Phi}^{(\mathrm{true})} = \{\boldsymbol{\phi}_1^{(\mathrm{true})}, \cdots, \boldsymbol{\phi}_K^{(\mathrm{true})}\}$、$\boldsymbol{\phi}_k^{(\mathrm{true})} = (\phi_{k1}^{(\mathrm{true})}, \cdots, \phi_{kV}^{(\mathrm{true})})$、$0 \leq \phi_{kv}^{(\mathrm{true})} \leq 1$、$\sum_{v=1}^V \phi_{kv}^{(\mathrm{true})} = 1$とする。または、値を直接指定する。
こちらは、トピックごとに各語彙の出現確率を生成するため、n引数にトピック数K、alpha引数にV個の正の値(ハイパーパラメータ)を指定する。
生成した値をデータフレームに格納する。
# 作図用のデータフレームを作成 true_word_df <- phi_true_kv |> tibble::as_tibble() |> tibble::add_column( topic = factor(1:K) # トピック番号 ) |> tidyr::pivot_longer( cols = !topic, names_to = "vocabulary", # 列名を語彙番号に変換 names_prefix = "V", names_ptypes = list(vocabulary = factor()), values_to = "probability" # 出現確率列をまとめる ) true_word_df
## # A tibble: 104 × 3 ## topic vocabulary probability ## <fct> <fct> <dbl> ## 1 1 1 0.0568 ## 2 1 2 0.0753 ## 3 1 3 0.0128 ## 4 1 4 0.0758 ## 5 1 5 0.0168 ## 6 1 6 0.104 ## 7 1 7 0.00690 ## 8 1 8 0.0149 ## 9 1 9 0.00799 ## 10 1 10 0.0514 ## # … with 94 more rows
phi_true_kvはK行V列のマトリクスなので、as_tibble()でデータフレームに変換して、トピック番号列を追加する。
さらに、pivot_longer()で語彙ごと(V列)の出現確率列を1列にまとめる。データフレームの変換時に列名が"V"に通し番号を付けた文字列になるので、names_prefix引数で先頭の"V"を取り除き、names_ptypes引数で因子型に変換して、語彙番号列とする。(列名のVはベクトルの頭文字？)
トピックごとの真の単語分布をグラフで確認する。
# 真の単語分布を作図 ggplot() + geom_bar(data = true_word_df, mapping = aes(x = vocabulary, y = probability, fill = vocabulary), stat = "identity", show.legend = FALSE) + # 単語分布 facet_wrap(topic ~ ., labeller = label_bquote(k==.(topic)), scales = "free_x") + # 分割 labs(title = "Word Distribution", subtitle = "truth", x = "v", y = expression(phi[kv]))
x軸を語彙番号、y軸を各語彙となる確率として棒グラフを描画する。
facet_wrap()にトピック番号列を指定して、トピックごとにグラフを分割して描画する。
文書の生成
設定した生成モデルに従って、文書を生成する。
・コード(クリックで展開)
# 文書を生成 z_d <- rep(NA, times = D) N_d <- rep(NA, times = D) N_dv <- matrix(NA, nrow = D, ncol = V) for(d in 1:D) { # 文書dのトピックを生成 one_hot_k <- rmultinom(n = 1, size = 1, prob = theta_true_k) |> as.vector() k <- which(one_hot_k == 1) # トピック番号を抽出 z_d[d] <- k # 単語数を決定 N_d[d] <- sample(100:150, size = 1) # 範囲を指定 # トピックkに従い単語を生成 N_dv[d, ] <- rmultinom(n = 1, size = N_d[d], prob = phi_true_kv[k, ]) |> as.vector() } z_d[1:5]; N_dv[1:5, 1:5]
## [1] 2 4 3 4 2 ## [,1] [,2] [,3] [,4] [,5] ## [1,] 0 17 10 5 1 ## [2,] 7 17 6 1 7 ## [3,] 1 17 5 6 2 ## [4,] 10 10 5 0 1 ## [5,] 0 9 6 6 0
for()を使って、文書ごとに処理する。
トピック分布(カテゴリ分布)に従いトピックを割り当てる。
カテゴリ分布の乱数は、多項分布の乱数生成関数rmultinom()の試行回数の引数sizeを1にすることで生成できる。サンプルサイズの引数nに1、パラメータの引数probにトピックの割り当て確率theta_true_kを指定する。
one-hotベクトルのトピック(割り当てられたトピック番号に対応する要素が1でそれ以外の要素が0のベクトル)が出力されるので、which()でトピック番号を抽出してkとする。
また確認用に、真のトピックとしてN_dに保存する。
単語数$N_d$をランダムに決めてN_dに保存する。全ての文書で同じ単語数を指定してもよい。
単語ごとに、単語分布(カテゴリ分布)に従い語彙を生成する。この例では、多項分布に従い$N_d$個の単語を1度に生成する。
rmultinom()のsize引数に文書dの単語数N_d[d]、prob引数にトピックkの単語分布phi_true_kv[k, ]を指定する。語彙ごとの出現度数$N_{dv}\ (v = 1, \dots, V)$(総和が$N_d$となるV個の値)が出力されるので、N_dvに保存する。
作図用に、文書ごとの出現度数をデータフレームに格納する。
# 作図用のデータフレームを作成 freq_df <- N_dv |> tibble::as_tibble() |> tibble::add_column( document = factor(1:D), # 文書番号 word_count = N_d, # 単語数 topic = factor(z_d) # 割り当てトピック ) |> tidyr::pivot_longer( cols = !c(document, word_count, topic), names_to = "vocabulary", # 列名を語彙番号に変換 names_prefix = "V", names_ptypes = list(vocabulary = factor()), values_to = "frequency" # 出現度数列をまとめる ) |> dplyr::mutate( relative_frequency = frequency / N_d[d] # 相対度数 ) freq_df
## # A tibble: 520 × 6 ## document word_count topic vocabulary frequency relative_frequency ## <fct> <int> <fct> <fct> <int> <dbl> ## 1 1 289 2 1 0 0 ## 2 1 289 2 2 17 0.0821 ## 3 1 289 2 3 10 0.0483 ## 4 1 289 2 4 5 0.0242 ## 5 1 289 2 5 1 0.00483 ## 6 1 289 2 6 0 0 ## 7 1 289 2 7 5 0.0242 ## 8 1 289 2 8 18 0.0870 ## 9 1 289 2 9 3 0.0145 ## 10 1 289 2 10 23 0.111 ## # … with 510 more rows
単語分布のときと同様にして、D行V列のN_dvをデータフレームに変換する。
因子型の1からDの整数を文書番号列、文書ごとの単語数列と真のトピック列を追加する。
pivot_longer()で語彙ごと(V列)の出現度数列を1列にまとめて、列名を語彙番号列に変換する。
文書ごとに、各語彙の度数を単語数で割って相対度数を求める。
各文書の語彙ごとの出現度数をグラフで確認する。
# 度数を作図 ggplot() + geom_bar(data = freq_df, mapping = aes(x = vocabulary, y = frequency, fill = vocabulary), stat = "identity", show.legend = FALSE) + # 出現度数 facet_wrap(document ~ ., labeller = label_bquote(list(d==.(document), N[d]==.(N_d[document]))), scales = "free_x") + # 分割 labs(title = "Word Frequency", x = "v", y = expression(N[dv]))
x軸を語彙番号、y軸を各語彙の度数として、文書ごとにグラフを分割して描画する。
・コード(クリックで展開)
確認用に、割り当てられた真のトピックに応じて、単語分布を複製する。
# 真のトピック分布を複製 tmp_word_df <- tidyr::expand_grid( document = factor(1:D), # 文書番号 true_word_df ) |> # 文書数分にトピック分布を複製 dplyr::filter(topic == z_d[document]) # 割り当てられたトピックを抽出 tmp_word_df
## # A tibble: 520 × 4 ## document topic vocabulary probability ## <fct> <fct> <fct> <dbl> ## 1 1 2 1 0.000643 ## 2 1 2 2 0.0448 ## 3 1 2 3 0.0276 ## 4 1 2 4 0.0184 ## 5 1 2 5 0.00131 ## 6 1 2 6 0.00562 ## 7 1 2 7 0.0271 ## 8 1 2 8 0.0633 ## 9 1 2 9 0.0150 ## 10 1 2 10 0.0843 ## # … with 510 more rows
文書番号と単語分布true_word_dfの行の全ての組み合わせをexpand_grid()で作成することで、true_word_dfをD個分に複製する。複製した単語分布(の行)から、真のトピックの単語分布(トピック番号topicが各文書の真のトピックz_d[document]と一致する行)のみを取り出す。
語彙ごとの相対度数を真の単語分布と重ねて描画する。
# 度数と単語分布を比較 ggplot() + geom_bar(data = freq_df, mapping = aes(x = vocabulary, y = relative_frequency, fill = vocabulary, linetype = "sample"), stat = "identity", alpha = 0.6) + # 出現度数 geom_bar(data = tmp_word_df, mapping = aes(x = vocabulary, y = probability, color = vocabulary, linetype = "truth"), stat = "identity", alpha = 0) + # 単語分布 facet_wrap(document ~ ., labeller = label_bquote(list(d==.(document), z[d]==.(z_d[document]))), scales = "free_x") + # 分割 scale_x_discrete(breaks = 1:V, labels = LETTERS[1:V]) + # x軸目盛:(雰囲気用の演出) scale_linetype_manual(breaks = c("sample", "truth"), values = c("solid", "dashed"), name = "distribution") + # 線の種類:(凡例表示用) guides(color = "none", fill = "none", linetype = guide_legend(override.aes = list(color = c("black", "black"), fill = c("white", "white")))) + # 凡例の体裁:(凡例表示用) labs(title = "Word Frequency", x = "v", y = expression(frac(N[dv], N[d])))
相対度数を塗りつぶし、真の出現確率を破線で示す。
x軸が語彙に対応しているのを表現するために、語彙番号ではなくアルファベットにした。
文書ごとにトピックが割り当てられており、それぞれに対応する単語分布に従って単語(語彙)が生成されているのを確認できる。
本来は、真のトピック分布と単語分布は分からないあるいは存在せず、割り当てられたトピックとトピック数も分からない。以降は、文書ごとの語彙の度数N_dvのみが手元にあるとして、トピック分布と単語分布を推定する。
推論処理
生成した文書データを用いて、トピック分布と単語分布を推定する。
パラメータの初期化
トピック数を指定して、文書数・語彙数を設定する。
# トピック数を指定 K <- 4 # 文書数と語彙数を取得 D <- nrow(N_dv) V <- ncol(N_dv) D; V
## [1] 20 ## [1] 26
トピック数$K$を指定する。
文書数$D$はN_dvの行数、語彙数$V$は列数に対応する。
負担率の受け皿を作成する。
# 負担率qの初期化 q_dk <- matrix(0, nrow = D, ncol = K) q_dk[1:5, ]
## [,1] [,2] [,3] [,4] ## [1,] 0 0 0 0 ## [2,] 0 0 0 0 ## [3,] 0 0 0 0 ## [4,] 0 0 0 0 ## [5,] 0 0 0 0
全ての値を0とするD行K列のマトリクスを$\mathbf{q} = \{\mathbf{q}_1, \cdots, \mathbf{q}_D\}$、$\mathbf{q}_d = (q_{d1}, \cdots, q_{dK})$の初期値とする。
トピック分布の初期値を作成する。
# トピック分布θの初期化 tmp_theta_k <- runif(n = K, min = 0, max = 1) theta_k <- tmp_theta_k / sum(tmp_theta_k) # 正規化 tmp_theta_k; theta_k
## [1] 0.880183187 0.960191366 0.003786403 0.156169213 ## [1] 0.440018953 0.480016440 0.001892889 0.078071718
K個の一様乱数をrunif()で生成し、総和で割って正規化して、トピック分布$\boldsymbol{\theta} = (\theta_1, \cdots, \theta_K)$、$0 \leq \theta_k \leq 1$、$\sum_{k=1}^K \theta_k = 1$の初期値とする。
「真の分布の設定」のときと同様にして、トピック分布の初期値をグラフで確認する。
同様に、単語分布の初期値を作成する。
# 単語分布Φの初期化 tmp_phi_kv <- runif(n = K*V, min = 0, max = 1) |> matrix(nrow = K, ncol = V) phi_kv <- tmp_phi_kv / rowSums(tmp_phi_kv) # 正規化 tmp_phi_kv[, 1:5]; phi_kv[, 1:5]
## [,1] [,2] [,3] [,4] [,5] ## [1,] 0.5504015 0.3243743 0.2130970 0.06740453 0.69602915 ## [2,] 0.1492291 0.6474596 0.7413559 0.13016941 0.03093222 ## [3,] 0.4249492 0.4766022 0.2187806 0.60908219 0.48546104 ## [4,] 0.3961126 0.3157212 0.2495844 0.67992133 0.47149214 ## [,1] [,2] [,3] [,4] [,5] ## [1,] 0.04337584 0.02556318 0.01679367 0.005311993 0.054852418 ## [2,] 0.01473401 0.06392638 0.07319714 0.012852167 0.003054066 ## [3,] 0.03477138 0.03899788 0.01790168 0.049838029 0.039722753 ## [4,] 0.03548589 0.02828399 0.02235910 0.060910989 0.042238788
K×V個の一様乱数を生成してK行V列のマトリクスに整形し、行ごとの和で割って正規化して、K個の単語分布$\boldsymbol{\Phi} = \{\boldsymbol{\phi}_1, \cdots, \boldsymbol{\phi}_K\}$、$\boldsymbol{\phi}_k = (\phi_{k1}, \cdots, \phi_{kV})$、$0 \leq \phi_{kv} \leq 1$、$\sum_{v=1}^V \phi_{kv} = 1$の初期値とする。
単語分布の初期値をグラフで確認する。
ここまでで、推論処理の準備ができた。続いて、推論処理を行う。
コード全体
試行回数を指定して、トピック分布と単語分布の更新を繰り返す。
# 試行回数を指定 MaxIter <- 20 # 推移の確認用の受け皿を作成 trace_theta_ki <- matrix(NA, nrow = K, ncol = MaxIter+1) trace_phi_kvi <- array(NA, dim = c(K, V, MaxIter+1)) # 初期値を保存 trace_theta_ki[, 1] <- theta_k trace_phi_kvi[, , 1] <- phi_kv # EMアルゴリズムによる最尤推定 for(i in 1:MaxIter){ # 次ステップのパラメータの初期化 new_theta_k <- rep(0, K) new_phi_kv <- matrix(0, nrow = K, ncol = V) for(d in 1:D){ ## (各文書) for(k in 1:K){ ## (各トピック) # 負担率qの計算:式(3.3) log_term_k <- log(theta_k + 1e-7) + colSums(N_dv[d, ] * log(t(phi_kv) + 1e-7)) # 分子 log_term_k <- log_term_k - min(log_term_k) # アンダーフロー対策 log_term_k <- log_term_k - max(log_term_k) # オーバーフロー対策 q_dk[d, k] <- exp(log_term_k[k]) / sum(exp(log_term_k)) # 正規化 # トピック分布θの計算:式(3.7)の分子 new_theta_k[k] <- new_theta_k[k] + q_dk[d, k] for(v in 1:V){ ## (各語彙) # 単語分布Φの計算:式(3.8)の分子 new_phi_kv[k, v] <- new_phi_kv[k, v] + q_dk[d, k] * N_dv[d, v] } ## (各語彙) } ## (各トピック) } ## (各文書) # パラメータの正規化と更新 theta_k <- new_theta_k / sum(new_theta_k) phi_kv <- new_phi_kv / rowSums(new_phi_kv) # i回目の更新値を保存 trace_theta_ki[, i+1] <- theta_k trace_phi_kvi[, , i+1] <- phi_kv # 動作確認 message("\r", i, "/", MaxIter, appendLF = FALSE) }
推移の確認用に、トピック分布と単語分布の初期値と更新値をそれぞれtrace_***に格納していく。
コードの解説
繰り返し行う更新処理をそれぞれ確認する。
Eステップ
Eステップでは、前ステップで更新した(または初期値の)トピック分布$\boldsymbol{\theta}$と単語分布$\boldsymbol{\Phi}$を用いて、負担率$\mathbf{q}$を更新する。
次の処理によって、負担率を更新(計算)する。
# 負担率qの初期化 q_dk <- matrix(0, nrow = D, ncol = K) for(d in 1:D) { for(k in 1:K) { # 負担率qの計算:式(3.3) log_term_k <- log(theta_k + 1e-7) + colSums(N_dv[d, ] * log(t(phi_kv + 1e-7))) # 分子 log_term_k <- log_term_k - min(log_term_k) # アンダーフロー対策 log_term_k <- log_term_k - max(log_term_k) # オーバーフロー対策 q_dk[d, k] <- exp(log_term_k[k]) / sum(exp(log_term_k)) # 正規化 } } q_dk[1:5, ]
## [,1] [,2] [,3] [,4] ## [1,] 2.646944e-74 1.237776e-193 1.122263e-88 1.000000e+00 ## [2,] 1.000000e+00 1.717865e-184 4.303884e-95 9.330041e-59 ## [3,] 8.336747e-59 1.000000e+00 8.139674e-132 7.833705e-82 ## [4,] 1.000000e+00 2.605016e-120 1.918943e-77 1.587688e-48 ## [5,] 4.908026e-52 6.313179e-186 1.488565e-81 1.000000e+00
負担率の各要素の更新式は、次の式である。
この式のまま計算すると、指数の計算においてアンダーフローの可能性がある。そこで、この式の分子を$\tilde{q}_{dk}$とおき、対数をとって計算する。
さらに、最小値を引くことでアンダーフロー、最大値を引くことでオーバーフローが起きにくくする。また、$\log 0$は計算できないので、結果に影響しないほどの小さい値1e-7を加えて計算する。
$x = \exp(\log x)$より、$\log$を外して総和で割って正規化する。
$\log \tilde{q}_{dk}\ (k = 1, \dots, K)$をlog_term_kとして、文書番号dとトピック番号kに関するfor()によって、q_dkのD×K個の要素を1回ずつ更新する。
Mステップ
Mステップでは、Eステップで更新した負担率$\mathbf{q}$を用いて、トピック分布$\boldsymbol{\theta}$と単語分布$\boldsymbol{\Phi}$を更新する。
次の処理によって、トピック分布を更新する。
# 次ステップのパラメータの初期化 new_theta_k <- rep(0, K) for(d in 1:D) { for(k in 1:K) { # トピック分布θの計算:式(3.7)の分子 new_theta_k[k] <- new_theta_k[k] + q_dk[d, k] } } new_theta_k
## [1] 6 3 3 8
トピック分布の各要素の更新式は、次の式である。
この式の分子を$\tilde{\theta}_k$とおいて計算する。
$\tilde{\theta}_k\ (k = 1, \dots, K)$をnew_theta_kとして、ステップごとに全ての要素を0に初期化しておく。トピック番号kに関するループにより、new_theta_kのK個の要素を計算する。
また、文書番号dに関するループにより、new_theta_k[k]にq_dk[, k]のＤ個の値を順番に足すことで、$\sum_d$の計算を行う。
全ての要素を計算できたら、正規化して、トピック分布の値を上書きする。
# パラメータの正規化と更新 res_theta_k <- new_theta_k / sum(new_theta_k) res_theta_k
## [1] 0.30 0.15 0.15 0.40
$\tilde{\theta}_k$を総和で割って正規化する。
(次にやる作図処理に影響しないように変数名をres_***にしている。)
続いて、次の処理によって、単語分布を更新する。
# 次ステップのパラメータの初期化 new_phi_kv <- matrix(0, nrow = K, ncol = V) for(d in 1:D) { for(k in 1:K) { for(v in 1:V) { # 単語分布Φの計算:式(3.8)の分子 new_phi_kv[k, v] <- new_phi_kv[k, v] + q_dk[d, k] * N_dv[d, v] } } } new_phi_kv[, 1:5]
## [,1] [,2] [,3] [,4] [,5] ## [1,] 77 87 36 6 23 ## [2,] 1 46 22 30 6 ## [3,] 61 57 5 49 7 ## [4,] 2 98 51 34 2
単語分布の各要素の更新式は、次の式である。
この式の分子を$\tilde{\phi}_{kv}$とおいて計算する。
$\tilde{\phi}_{kv}\ (k = 1, \dots, K,\ v = 1, \dots, V)$をnew_phi_kvとして、ステップごとに全ての要素を0に初期化しておく。トピック番号kと語彙番号vについてのループによって、new_phi_kvのK×V個の要素を計算する。
また、文書番号dのループにより、new_phi_kv[k, v]にq_dk[, k] * N_dv[, v]のＤ個の値を順番に足すことで、$\sum_d$の計算を行う。
全ての要素を計算できたら、正規化して、単語分布の値を上書きする。
# パラメータの正規化と更新 res_phi_kv <- new_phi_kv / rowSums(new_phi_kv) res_phi_kv[, 1:K]
## [,1] [,2] [,3] [,4] ## [1,] 0.049358974 0.05576923 0.023076923 0.003846154 ## [2,] 0.001386963 0.06380028 0.030513176 0.041608877 ## [3,] 0.090504451 0.08456973 0.007418398 0.072700297 ## [4,] 0.001050972 0.05149764 0.026799790 0.017866527
$\tilde{\phi}_{kv}$を$v$(行)についての和で割って正規化する。
以上が、1ステップでの更新処理である。この処理を指定した回数繰り返し行う。
推論結果の可視化
推定したトピック分布と単語分布を可視化する。
推定した分布の確認
「真の分布の設定」のときと同様にして、推定したトピック分布のグラフを作成する。
推定した単語分布のグラフを作成する。
ここまでで、トピック分布と単語分布を推定できた。続いて、更新推移の様子や真の分布との対応関係を確認する。
更新推移の確認
・コード(クリックで展開)
各試行におけるトピック分布の更新値をデータフレームに格納する。
# 作図用のデータフレームを作成 trace_topic_df <- trace_theta_ki |> tibble::as_tibble() |> tibble::add_column( topic = factor(1:K) # トピック番号 ) |> tidyr::pivot_longer( cols = !topic, names_to = "iteration", # 列名を試行番号に変換 names_prefix = "V", names_transform = list(iteration = as.numeric), values_to = "probability" # 割り当て確率列をまとめる ) |> dplyr::mutate( iteration = iteration - 1 ) trace_topic_df
## # A tibble: 84 × 3 ## topic iteration probability ## <fct> <dbl> <dbl> ## 1 1 0 0.440 ## 2 1 1 0.182 ## 3 1 2 0.308 ## 4 1 3 0.3 ## 5 1 4 0.3 ## 6 1 5 0.3 ## 7 1 6 0.3 ## 8 1 7 0.3 ## 9 1 8 0.3 ## 10 1 9 0.3 ## # … with 74 more rows
K行MaxIter + 1列のマトリクスtrace_theta_kiをデータフレームに変換して、トピック番号列を追加する。
試行ごとの割り当て確率列を1列にまとめて、列名を試行番号に変換する。初期値の試行番号を0で表すために、試行番号列から1を引き、0からMaxIterの整数になるようにする。
トピック分布の更新値の推移の折れ線グラフを作成する。
# トピック分布の推移を作図 ggplot() + geom_line(data = trace_topic_df, mapping = aes(x = iteration, y = probability, color = topic), alpha = 0.5, size = 1) + # 更新推移 labs(title = "Topic Distribution", subtitle = "EM Algorithm", color = "k", x = "iteration", y = expression(theta[k]))
x軸を試行回数、y軸を推定した確率(更新値)として折れ線グラフを描画する。
・コード(クリックで展開)
各試行における単語分布の更新値をデータフレームに格納する。
# 作図用のデータフレームを作成 trace_word_df <- trace_phi_kvi |> tibble::as_tibble() |> tibble::add_column( topic = factor(1:K) # トピック番号 ) |> tidyr::pivot_longer( cols = !topic, names_to = "vi", # 列名を語彙番号×試行番号に変換 names_prefix = "V", names_transform = list("vi" = as.numeric), values_to = "probability" # 出現確率列をまとめる ) |> dplyr::arrange(vi, topic) |> dplyr::bind_cols( tidyr::expand_grid( iteration = 0:MaxIter, # 試行番号 vocabulary = factor(1:V), # 語彙番号 tmp_topic = 1:K ) |> # 試行ごとに語彙番号を複製 dplyr::select(!tmp_topic) ) trace_word_df
## # A tibble: 2,184 × 5 ## topic vi probability iteration vocabulary ## <fct> <dbl> <dbl> <int> <fct> ## 1 1 1 0.0434 0 1 ## 2 2 1 0.0147 0 1 ## 3 3 1 0.0348 0 1 ## 4 4 1 0.0355 0 1 ## 5 1 2 0.0256 0 2 ## 6 2 2 0.0639 0 2 ## 7 3 2 0.0390 0 2 ## 8 4 2 0.0283 0 2 ## 9 1 3 0.0168 0 3 ## 10 2 3 0.0732 0 3 ## # … with 2,174 more rows
K行V列のマトリクスをMaxIter + 1個分並べた3次元配列trace_phi_kviをデータフレームに変換すると、trace_phi_kvi[, , 1]の隣の列にtrace_phi_kvi[, , 2]のように並べたデータフレームになる。
トピック番号を追加して、試行と語彙ごとの出現確率列を1列にまとめて、試行番号と語彙番号を追加する。
単語分布の更新値の推移の折れ線グラフを作成する。
# 単語分布の推移を作図 ggplot() + geom_line(data = trace_word_df, mapping = aes(x = iteration, y = probability, color = vocabulary), alpha = 0.5, size = 1) + # 更新推移 facet_wrap(topic ~ ., labeller = label_bquote(k==.(topic))) + # 分割 labs(title = "Topic Distribution", subtitle = "EM Algorithm", color = "v", x = "iteration", y = expression(phi[kv]))
トピックごとにグラフを分割して、x軸を試行回数、y軸を推定した確率(更新値)として折れ線グラフを描画する。
値が収束する様子を確認できる。
・コード(クリックで展開)
トピック分布のアニメーションを作成する。
# 推定トピック分布のアニメーションを作図 anime_topic_graph <- ggplot() + geom_bar(data = trace_topic_df, mapping = aes(x = topic, y = probability, fill = topic), stat = "identity", show.legend = FALSE) + # トピック分布 gganimate::transition_manual(frames = iteration) + # フレーム labs(title = "Topic Distribution", subtitle = "iteration = {current_frame}", x = "k", y = expression(theta[k])) # gif画像を作成 gganimate::animate(plot = anime_topic_graph, nframes = MaxIter+1, fps = 5, width = 600, height = 450)
gganimateパッケージを利用して、アニメーション(gif画像)を作成する。
transition_manual()のフレーム制御の引数framesに試行回数列iterationを指定して、グラフを作成する。
animate()のplot引数にグラフオブジェクト、nframes引数に初期値分を加えた試行回数MaxIter + 1を指定して、gif画像を作成する。また、fps引数に1秒当たりのフレーム数を指定できる。
同様に、単語分布のアニメーションを作成する。
# 推定単語分布のアニメーションを作図 anime_word_graph <- ggplot() + geom_bar(data = trace_word_df, mapping = aes(x = vocabulary, y = probability, fill = vocabulary), stat = "identity", show.legend = FALSE) + # 単語分布 gganimate::transition_manual(frames = iteration) + # フレーム facet_wrap(topic ~ ., labeller = label_bquote(k==.(topic)), scales = "free_x") + # 分割 labs(title = "Word Distribution", subtitle = "iteration = {current_frame}", x = "v", y = expression(phi[kv])) # gif画像を作成 gganimate::animate(plot = anime_word_graph, nframes = MaxIter+1, fps = 5, width = 800, height = 600)
ランダムな初期値が、徐々に偏りを持つように変化するのを確認できる。
真の分布との比較
最後に、文書データの生成モデル(真の分布)と推定した分布との対応関係を考える。詳しくは「KL情報量が最小となるように確率分布の次元を入れ替えたい - からっぽのしょこ」を参照のこと。
・コード(クリックで展開)
真の分布と推定した分布のKL情報量を計算する。
# インデックスの組み合わせを作成してKL情報量を計算 kl_df <- tidyr::expand_grid( group = 1:gamma(K+1), # 組み合わせ番号 k = factor(1:K) # 真の分布・並べ替え後のトピック番号 ) |> # トピック番号を複製 tibble::add_column( j = gtools::permutations(n = K, r = K, v = 1:K) |> # トピック番号の順列を作成 t() |> as.vector() |> factor() # 元のトピック番号・並べ替え用のインデックス ) |> dplyr::group_by(group) |> # 組み合わせごとの計算用 dplyr::mutate( # KL情報量を計算 kl_topic = sum(theta_true_k * (log(theta_true_k) - log(theta_k[j]))), # トピック分布のKL情報量 kl_word = sum(phi_true_kv * (log(phi_true_kv) - log(phi_kv[j, ]))), # 単語分布のKL情報量 kl_sum = kl_topic + kl_word/V ) |> dplyr::ungroup() |> dplyr::arrange(kl_sum, k) # 当て嵌まりが良い順に並べ替え kl_df
## # A tibble: 96 × 6 ## group k j kl_topic kl_word kl_sum ## <int> <fct> <fct> <dbl> <dbl> <dbl> ## 1 18 1 3 0.0318 0.820 0.0633 ## 2 18 2 4 0.0318 0.820 0.0633 ## 3 18 3 2 0.0318 0.820 0.0633 ## 4 18 4 1 0.0318 0.820 0.0633 ## 5 13 1 3 0.0220 2.49 0.118 ## 6 13 2 1 0.0220 2.49 0.118 ## 7 13 3 2 0.0220 2.49 0.118 ## 8 13 4 4 0.0220 2.49 0.118 ## 9 24 1 4 0.107 2.43 0.200 ## 10 24 2 3 0.107 2.43 0.200 ## # … with 86 more rows
単語分布のKL情報量$\mathrm{KL}(\boldsymbol{\theta}^{(\mathrm{truth})}, \boldsymbol{\theta}) = \sum_{k=1}^K \theta_k^{(\mathrm{truth})} \log \frac{\theta_k^{(\mathrm{truth})}}{\theta_k}$とトピック分布のKL情報量$\mathrm{KL}(\boldsymbol{\Phi}^{(\mathrm{truth})}, \boldsymbol{\Phi}) = \sum_{k=1}^K \sum_{v=1}^V \phi_{kv}^{(\mathrm{truth})} \log \frac{\phi_{kv}^{(\mathrm{truth})}}{\phi_{kv}}$を足してkl_sum列とする。ただし単語分布に関して、$v$についての和の影響を取り除くために、$V$で割った値を足すことにする。
真の分布を推定(近似)できている組み合わせが上にくるように、KL情報量が小さい順に行を並べ替える。
KL情報量が最小になるインデックスを取り出す。
# KL情報量が最小となるトピック番号を抽出 adapt_idx <- kl_df |> dplyr::slice_head(n = K) |> # 最小の組み合わせを抽出 dplyr::pull(j) |> # 列をベクトルに変換 as.numeric() # 数値に変換 adapt_idx
## [1] 3 4 2 1
kl_dfの上からK行分のj列をベクトルとして取り出してadapt_idxとする。
推定したトピック分布のトピック番号を入れ替えてデータフレームに格納する。
# 作図用のデータフレームを作成 topic_df <- tibble::tibble( k = factor(adapt_idx), # 推定時のトピック番号 topic = factor(1:K), # 真の分布に対応したトピック番号 probability = theta_k[adapt_idx] # 割り当て確率 ) topic_df
## # A tibble: 4 × 3 ## k topic probability ## <fct> <fct> <dbl> ## 1 3 1 0.15 ## 2 4 2 0.4 ## 3 2 3 0.15 ## 4 1 4 0.3
adapt_idxを使ってtheta_kの要素を入れ替えてから、これまでと同様にして処理する。
トピック分布に真の分布を重ねてグラフを作成する。
# トピック分布を作図 ggplot() + geom_bar(data = topic_df, mapping = aes(x = topic, y = probability, fill = k, linetype = "result"), stat = "identity", alpha = 0.5) + # 推定トピック分布 geom_bar(data = true_topic_df, mapping = aes(x = topic, y = probability, color = topic, linetype = "truth"), stat = "identity", alpha = 0) + # 真のトピック分布 scale_linetype_manual(breaks = c("result", "truth"), values = c("solid", "dashed"), name = "distribution") + # 線の種類:(凡例表示用) guides(color = "none", fill = "none", linetype = guide_legend(override.aes = list(color = c("black", "black"), fill = c("white", "white")))) + # 凡例の体裁:(凡例表示用) labs(title = "Topic Distribution", subtitle = "EM Algorithm", x = "k", y = expression(theta[k]))
推定した分布を塗りつぶし、真の分布を破線で示す。
推定した単語分布のトピック番号を入れ替えてデータフレームに格納する。
# 作図用のデータフレームを作成 word_df <- phi_kv[adapt_idx, ] |> tibble::as_tibble() |> tibble::add_column( k = factor(adapt_idx), # 推定時のトピック番号 topic = factor(1:K) # 真の分布に対応したトピック番号 ) |> tidyr::pivot_longer( cols = !c(k, topic), names_prefix = "V", names_to = "vocabulary", # 列名を語彙番号に変換 names_ptypes = list(vocabulary = factor()), values_to = "probability" # 出現確率列をまとめる ) word_df
## # A tibble: 104 × 4 ## k topic vocabulary probability ## <fct> <fct> <fct> <dbl> ## 1 3 1 1 0.0905 ## 2 3 1 2 0.0846 ## 3 3 1 3 0.00742 ## 4 3 1 4 0.0727 ## 5 3 1 5 0.0104 ## 6 3 1 6 0.128 ## 7 3 1 7 0.00742 ## 8 3 1 8 0.0134 ## 9 3 1 9 0.00445 ## 10 3 1 10 0.0534 ## # … with 94 more rows
adapt_idxを使ってphi_kvの行を入れ替えてから、これまでと同様にして処理する。
単語分布を真の分布に重ねてのグラフを作成する。
# 単語分布を作図 ggplot() + geom_bar(data = word_df, mapping = aes(x = vocabulary, y = probability, fill = vocabulary, linetype = "result"), stat = "identity", alpha = 0.5) + # 推定単語分布 geom_bar(data = true_word_df, mapping = aes(x = vocabulary, y = probability, color = vocabulary, linetype = "truth"), stat = "identity", alpha = 0) + # 真の単語分布 facet_wrap(topic ~ ., labeller = label_bquote(k==.(topic)), scales = "free_x") + # 分割 scale_linetype_manual(breaks = c("result", "truth"), values = c("solid", "dashed"), name = "distribution") + # 線の種類:(凡例表示用) guides(color = "none", fill = "none", linetype = guide_legend(override.aes = list(color = c("black", "black"), fill = c("white", "white")))) + # 凡例の体裁:(凡例表示用) labs(title = "Word Distribution", subtitle = "EM Algorithm", x = "v", y = expression(phi[kv]))
真の分布を推定(近似)できているかを確認できる。
・コード(クリックで展開)
各試行におけるトピック分布の更新値をデータフレームに格納する。
# 作図用のデータフレームを作成 trace_topic_df <- trace_theta_ki[adapt_idx, ] |> tibble::as_tibble() |> tibble::add_column( k = factor(adapt_idx), # 推定時のトピック番号 topic = factor(1:K) # 真の分布に対応付けたトピック番号 ) |> tidyr::pivot_longer( cols = !c(k, topic), names_to = "iteration", # 列名を試行番号に変換 names_prefix = "V", names_transform = list(iteration = as.numeric), values_to = "probability" # 割り当て確率列をまとめる ) |> dplyr::mutate( iteration = iteration - 1 ) trace_topic_df
## # A tibble: 84 × 4 ## k topic iteration probability ## <fct> <fct> <dbl> <dbl> ## 1 3 1 0 0.00189 ## 2 3 1 1 0.363 ## 3 3 1 2 0.0957 ## 4 3 1 3 0.15 ## 5 3 1 4 0.15 ## 6 3 1 5 0.15 ## 7 3 1 6 0.15 ## 8 3 1 7 0.15 ## 9 3 1 8 0.15 ## 10 3 1 9 0.15 ## # … with 74 more rows
adapt_idxを使ってtrace_theta_kiの行を入れ替えてから、「更新推移の可視化」のときと同様にして処理する。
真のトピック分布をフレーム数分に複製する。
# 真のトピック分布を複製 anime_true_topic_df <- tidyr::expand_grid( iteration = 0:MaxIter, # 試行番号 true_topic_df ) anime_true_topic_df
## # A tibble: 84 × 3 ## iteration topic probability ## <int> <fct> <dbl> ## 1 0 1 0.212 ## 2 0 2 0.288 ## 3 0 3 0.177 ## 4 0 4 0.322 ## 5 1 1 0.212 ## 6 1 2 0.288 ## 7 1 3 0.177 ## 8 1 4 0.322 ## 9 2 1 0.212 ## 10 2 2 0.288 ## # … with 74 more rows
試行番号と真のトピック分布true_topic_dfの行の全ての組み合わせを作成することで、MaxIter + 1個に複製する。
トピック分布のアニメーションを作成する。
# トピック分布のアニメーションを作図 anime_topic_graph <- ggplot() + geom_bar(data = trace_topic_df, mapping = aes(x = topic, y = probability, fill = k, linetype = "result"), stat = "identity", alpha = 0.5) + # 推定トピック分布 geom_bar(data = anime_true_topic_df, mapping = aes(x = topic, y = probability, color = topic, linetype = "truth"), stat = "identity", alpha = 0) + # 真のトピック分布 gganimate::transition_manual(frames = iteration) + # フレーム scale_linetype_manual(breaks = c("result", "truth"), values = c("solid", "dashed"), name = "distribution") + # 線の種類:(凡例表示用) guides(color = "none", fill = "none", linetype = guide_legend(override.aes = list(color = c("black", "black"), fill = c("white", "white")))) + # 凡例の体裁:(凡例表示用) labs(title = "Topic Distribution", subtitle = "iteration = {current_frame}", x = "k", y = expression(theta[k])) # gif画像を作成 gganimate::animate(plot = anime_topic_graph, nframes = MaxIter+1, fps = 5, width = 700, height = 450)
更新を繰り返すことで真の分布に近付く様子を確認できる。
・コード(クリックで展開)
各試行における単語分布の更新値をデータフレームに格納する。
# 作図用のデータフレームを作成 trace_word_df <- trace_phi_kvi[adapt_idx, , ] |> tibble::as_tibble() |> tibble::add_column( k = factor(adapt_idx), # 推定時のトピック番号 topic = factor(1:K) # 真の分布に対応付けたトピック番号 ) |> tidyr::pivot_longer( cols = !c(k, topic), names_to = "vi", # 列名を語彙番号×試行番号に変換 names_prefix = "V", names_transform = list("vi" = as.numeric), values_to = "probability" # 出現確率列をまとめる ) |> dplyr::arrange(vi, topic) |> dplyr::bind_cols( tidyr::expand_grid( iteration = 0:MaxIter, # 試行番号 vocabulary = factor(1:V), # 語彙番号 tmp_topic = 1:K ) |> # 試行ごとに語彙番号を複製 dplyr::select(!tmp_topic) ) trace_word_df
## # A tibble: 2,184 × 6 ## k topic vi probability iteration vocabulary ## <fct> <fct> <dbl> <dbl> <int> <fct> ## 1 3 1 1 0.0348 0 1 ## 2 4 2 1 0.0355 0 1 ## 3 2 3 1 0.0147 0 1 ## 4 1 4 1 0.0434 0 1 ## 5 3 1 2 0.0390 0 2 ## 6 4 2 2 0.0283 0 2 ## 7 2 3 2 0.0639 0 2 ## 8 1 4 2 0.0256 0 2 ## 9 3 1 3 0.0179 0 3 ## 10 4 2 3 0.0224 0 3 ## # … with 2,174 more rows
adapt_idxを使ってtrace_phi_kviの行を入れ替えてから、「更新推移の可視化」のときと同様にして処理する。
真の単語分布をフレーム数分に複製する。
# 真の単語分布を複製 anime_true_word_df <- tidyr::expand_grid( iteration = 0:MaxIter, # 試行番号 true_word_df ) anime_true_word_df
## # A tibble: 2,184 × 4 ## iteration topic vocabulary probability ## <int> <fct> <fct> <dbl> ## 1 0 1 1 0.0568 ## 2 0 1 2 0.0753 ## 3 0 1 3 0.0128 ## 4 0 1 4 0.0758 ## 5 0 1 5 0.0168 ## 6 0 1 6 0.104 ## 7 0 1 7 0.00690 ## 8 0 1 8 0.0149 ## 9 0 1 9 0.00799 ## 10 0 1 10 0.0514 ## # … with 2,174 more rows
単語分布のアニメーションを作成する。
# 単語分布のアニメーションを作図 anime_word_graph <- ggplot() + geom_bar(data = trace_word_df, mapping = aes(x = vocabulary, y = probability, fill = vocabulary, linetype = "result"), stat = "identity", alpha = 0.5) + # 推定単語分布 geom_bar(data = anime_true_word_df, mapping = aes(x = vocabulary, y = probability, color = vocabulary, linetype = "truth"), stat = "identity", alpha = 0) + # 真の単語分布 gganimate::transition_manual(frames = iteration) + # フレーム facet_wrap(topic ~ ., labeller = label_bquote(k==.(topic)), scales = "free_x") + # 分割 scale_linetype_manual(breaks = c("result", "truth"), values = c("solid", "dashed"), name = "distribution") + # 線の種類:(凡例表示用) guides(color = "none", fill = "none", linetype = guide_legend(override.aes = list(color = c("black", "black"), fill = c("white", "white")))) + # 凡例の体裁:(凡例表示用) labs(title = "Word Distribution", subtitle = "iteration = {current_frame}", x = "v", y = expression(phi[kv])) # gif画像を作成 gganimate::animate(plot = anime_word_graph, nframes = MaxIter+1, fps = 5, width = 900, height = 600)
この記事では、EMアルゴリズムによる最尤推定を実装した。次の記事では、変分ベイズ推論を実装する。
参考書籍
- 岩田具治(2015)『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社
おわりに
章が進み数式イジイジにも慣れてきて、内容もふんわり理解したつもりにもなってきたので次は実装です。半月ほどコネコネしてたらループが回るようになりました。やはりR楽しい。
参考書のアルゴリズム3.1を参考にしてるので'for()'多めです。もう少しRっぽいコードでも再現したいです。
数式でもRでも日本語でも説明できればそれなりに理解したと言ってもいいでしょうかね…
2019/08/12：加筆修正しました。
- 2023/02/07：加筆修正しました。
数式編の修正から随分空いてしまいましたが、プログラム編もようやく加筆修正できました(し始めました)。
主な変更は次の4点です。
- 推定できてるのか確認できるように、人工データを使うようにしました。
- 真の分布を設定したことにより、可視化を充実させました。
- 負担率の計算において、アンダーフロー(による0除算)対策をしました。
- 作図処理において、forループを関数による処理に置き換えました。
作図に関して、ループ処理をtidyverse関数によって置き替えたのでモダンなRっぽいコードになったのではないでしょうか。ついでに、推論処理でもループを使わずに実装したのを別記事で書きました。どっちの方が分かりやすいかは人によると思うので、そっちも読んでみてください。
【次節の内容】