はじめに
機械学習プロフェッショナルシリーズの『トピックモデル』の勉強時に理解の助けになったことや勉強会用レジュメのまとめです。
この記事では、Rで混合ユニグラムモデルをEMアルゴリズムを用いて最尤推定する方法について書いています。Rの基本的な関数で参考書のアルゴリズム3.1を組むことで、理論をコード・処理の面から理解しようというものです。
パラメータの更新式の導出などはこの記事をご参照ください。
【他の節一覧】
【この節の内容】
・コード全体
・テキスト処理
# 利用パッケージ library(RMeCab) library(tidyverse) ### テキスト処理 # 抽出しない単語を指定 stop_words <- "[a-z]" # 小文字のアルファベットを含む語 # 形態素解析 mecab_df <- docDF("フォルダ名", type = 1) # テキストファイルの保存先を指定する # 文書dの語彙vの出現回数(N_dv)の集合 N_dv <- mecab_df %>% filter(grepl("名詞|形容詞|^動詞", POS1)) %>% # 抽出する品詞(大分類)を指定する filter(grepl("一般|^自立", POS2)) %>% # 抽出する品詞(細分類)を指定する filter(!grepl(stop_words, TERM)) %>% # ストップワードを除く select(-c(TERM, POS1, POS2)) %>% # 数値列のみを残す filter(apply(., 1, sum) >= 5) %>% # 抽出する総出現回数を指定する t() # 転置 # 確認用の行列名 dimnames(N_dv) <- list(paste0("d=", 1:nrow(N_dv)), # 行名 paste0("v=", 1:ncol(N_dv))) # 列名 # 文書dの単語数(N_d)のベクトル N_d <- apply(N_dv, 1, sum) # 行方向に和をとる # 文書数(D) D <- nrow(N_dv) # 総語彙数(V) V <- ncol(N_dv)
・パラメータの初期設定
### パラメータの初期設定 # トピック数(K) K <- 4 # 任意の値を指定する # 負担率(q_dk)の集合 q_dk <- matrix(0, nrow = D, ncol = K, dimnames = list(paste0("d=", 1:D), # 行名 paste0("k=", 1:K))) # 列名 ## トピック分布(theta) # 初期値をランダムに設定 theta_k <- sample(seq(0.1, 10, 0.1), size = K, replace = TRUE) # 正規化：sum_{k=1}^K theta_k = 1 theta_k <- theta_k / sum(theta_k) names(theta_k) <- paste0("k=", 1:K) # 確認用の要素名 ## 単語分布(phi) # 初期値をランダムに設定 phi_kv <- sample(seq(0.1, 10, 0.1), size = K * V, replace = TRUE) %>% matrix(nrow = K, ncol = V, dimnames = list(paste0("k=", 1:K), # 行名 paste0("v=", 1:V))) # 列名 # 正規化：sum_{v=1}^V phi_{kv} = 1 for(k in 1:K) { phi_kv[k, ] <- phi_kv[k, ] / sum(phi_kv[k, ]) }
・最尤推定
### 最尤推定 # 推移の確認用 trace_theta <- as.matrix(theta_k) trace_phi <- as.matrix(phi_kv[1, ]) for(i in 1:50){ # 任意の回数を指定する # 次ステップのパラメータを初期化 new_theta_k <- rep(0, K) names(new_theta_k) <- paste0("k=", 1:K) new_phi_kv <- matrix(0, nrow = K, ncol = V, dimnames = list(paste0("k=", 1:K), paste0("v=", 1:V))) for(d in 1:D){ ## (各文書) for(k in 1:K){ ## (各トピック) # 負担率の計算 tmp_phi_k <- t(phi_kv) ^ N_dv[d, ] %>% apply(2, prod) q_dk[d, k] <- theta_k[k] * tmp_phi_k[k] / sum(theta_k * tmp_phi_k) # トピック分布(theta_k)を更新 new_theta_k[k] <- new_theta_k[k] + q_dk[d, k] for(v in 1:V){ ## (各語彙) # 単語分布(phi_kv)を更新 new_phi_kv[k, v] <- new_phi_kv[k, v] + q_dk[d, k] * N_dv[d, v] } ## (/各語彙) } ## (/各トピック) } ## (/各文書) # パラメータの更新 theta_k <- new_theta_k phi_kv <- new_phi_kv # パラメータの正規化 theta_k <- theta_k / sum(theta_k) for(k2 in 1:K) { phi_kv[k2, ] <- phi_kv[k2, ] / sum(phi_kv[k2, ]) } # 推移の確認用 trace_theta <- cbind(trace_theta, as.matrix(theta_k)) trace_phi <- cbind(trace_phi, as.matrix(phi_kv[1, ])) }
・コードの解説
・利用パッケージ
# 利用パッケージ library(RMeCab) library(tidyverse)
形態素解析のためのRMeCab::docDF()と、グラフ作成時のtidyr::gather()、{ggplot2}関連の関数以外にライブラリの読み込みが必要なのはdplyrのみ。(面倒なのでtidyversパッケージを読み込んでいる)
・テキスト処理
・形態素解析
# 形態素解析 mecab_df <- docDF("フォルダ名", type = 1) # テキストファイルの保存先を指定する
RMeCab::docDF()で、形態素解析を行う。第1引数に解析を行うテキストファイルがあるディレクトリのパスを指定する。引数にtype = 1を指定すると、テキストを単語に分割する。デフォルトの設定は0で、type = 0だと文字に分割する。
ちなみにこの一連の記事で用いているテキストは、ハロー！プロジェクトのグループJuice=Juiceの楽曲の歌詞です！
head(mecab_df[, 1:7])
## TERM POS1 POS2 juice_01.txt juice_02.txt juice_03.txt juice_04.txt ## 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 ## 6 ) 名詞 サ変接続 0 0 0 0
解析結果はこのように、単語とその品詞、文書ごとの頻度の列からなるデータフレームで出力される。(表示は7列目までに省略している)
・特徴語の抽出
# 抽出しない単語を指定 stop_words <- "[a-z]" # 小文字のアルファベットを含む語 # 文書dの語彙vの出現回数(N_dv)の集合 N_dv <- mecab_df %>% filter(grepl("名詞|形容詞|^動詞", POS1)) %>% # 抽出する品詞(大分類)を指定する filter(grepl("一般|^自立", POS2)) %>% # 抽出する品詞(細分類)を指定する filter(!grepl(stop_words, TERM)) %>% # ストップワードを除く select(-c(TERM, POS1, POS2)) %>% # 数値列のみを残す filter(apply(., 1, sum) >= 5) %>% # 抽出する総出現回数を指定する t() # 転置 # 確認用の行列名 dimnames(N_dv) <- list(paste0("d=", 1:nrow(N_dv)), # 行名 paste0("v=", 1:ncol(N_dv))) # 列名
形態素解析の結果から必要な情報を抽出する。
POS1列(品詞大分類)とPOS2列(品詞小分類)の情報を用いて、filter(grepl())に品詞を指定し必要な語彙(行)を取り出す。ここでは、助動詞を含まないようにするために動詞の前に文字列の先頭であることを示す正規表現の記号の^を付けている。同様に非自立を除くため、自立の頭にも付けている。
grepl()の頭に!を付けることで、grepl()に指定した条件以外の行を取り出すことができる。ここでは、小文字のアルファベットのaからzまでを表す正規表現[a-z]を使って、アルファベットを含む語を除く。
apply()関数は、apply(データ, 向き, 関数名)とすることで、行列等のデータに対して指定した向き(1なら行、2なら列)に指定した関数を適応する。第1引数の.はパイプ(%>%)で渡している元のデータを指定する省略記号である。つまり、mecab_dfの各行の要素に対してsum()を行うことで、文書全体での各語彙の総出現回数$N_v$を求めている。その結果が、指定した回数以上である行(語彙)をfilter()で抽出する。
語彙が行、文書が列となっているので、t()で転置して$N_{dv}$とする。
dimnames()で行名と列名を付ける。各次元の名前の要素はlist形式で渡す。ただし、これは目で確認するためのもので、付けなくても処理に影響しない。他の変数についても同様である。
head(N_dv[, 1:15])
## v=1 v=2 v=3 v=4 v=5 v=6 v=7 v=8 v=9 v=10 v=11 v=12 v=13 v=14 v=15 ## d=1 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 ## d=2 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 ## d=3 0 0 0 0 0 0 0 0 0 1 0 0 0 7 0 ## d=4 0 0 0 0 0 0 0 1 0 5 0 0 0 2 4 ## d=5 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 ## d=6 0 1 0 0 0 0 0 1 2 0 0 0 0 1 0
15列目までに省略して表示している。
# 文書dの単語数(N_d)のベクトル N_d <- apply(N_dv, 1, sum) # 行方向に和をとる
apply(N_dv, 1, sum)で、N_dvの各行の要素に対してsum()を行い$N_d = \sum_{v=1}^V N_{dv}$の計算をして、各文書の単語数$N_d$を求める。(この節では利用しない)\
N_d
## d=1 d=2 d=3 d=4 d=5 d=6 d=7 d=8 d=9 d=10 d=11 d=12 d=13 d=14 d=15 ## 45 28 34 34 28 47 8 60 36 59 37 49 24 41 61 ## d=16 d=17 d=18 d=19 d=20 d=21 d=22 d=23 d=24 d=25 d=26 d=27 d=28 d=29 d=30 ## 36 27 25 38 34 28 54 48 47 39 23 32 40 35 20 ## d=31 d=32 d=33 d=34 d=35 d=36 d=37 d=38 d=39 d=40 d=41 d=42 d=43 ## 24 30 33 40 48 45 80 16 47 30 89 45 36
要素が0になった文書があると、推定時にエラーが起きることがあるので、語彙を抽出する条件を緩めるかその文書を除いて次に進む。
# 文書数(D) D <- nrow(N_dv) # 総語彙数(V) V <- ncol(N_dv)
N_dvの行数が文書数$D$に、列数が語彙数$V$に相当する。
・パラメータの初期設定
・トピック数
# トピック数(K) K <- 4 # 任意の値を指定する
任意のトピック数を指定する。
# 負担率(q_dk)の集合 q_dk <- matrix(0, nrow = D, ncol = K, dimnames = list(paste0("d=", 1:D), # 行名 paste0("k=", 1:K))) # 列名
全ての要素が0である$D$行$K$列の行列を作成する。
行列(マトリクス)は、matrix()を使って作る。第1引数に0を指定し、行数を指定する引数nrowに文書数のDを、列数を指定する引数ncolにトピック数のKを指定する。
行名・列名はdimneames =引数にlist形式で指定する。
## トピック分布(theta) # 初期値をランダムに設定 theta_k <- sample(seq(0.1, 10, 0.1), size = K, replace = TRUE) # 正規化：sum_{k=1}^K theta_k = 1 theta_k <- theta_k / sum(theta_k) names(theta_k) <- paste0("k=", 1:K) # 確認用の要素名
初期値をランダムに設定する。
seq(最小値, 最大値, 間隔)でランダムに取り得る範囲を指定する。そこから、sample()を使って、ランダムに値を取る。第2引数のsizeにK(トピック数)を指定し、要素が$K$個のベクトルを作成する。replace = TRUEを指定すると、ランダムに取り得るに重複を許す\
ただし、$\sum_{k=1}^K \theta_k = 1$の条件を満たすために、合計値で割って正規化する。(apply(dat, 1, sum)当たりで上手いことできると思うのだけどとりあえずfor()で…)\
## 単語分布(phi) # 初期値をランダムに設定 phi_kv <- sample(seq(0.1, 10, 0.1), size = K * V, replace = TRUE) %>% matrix(nrow = K, ncol = V, dimnames = list(paste0("k=", 1:K), # 行名 paste0("v=", 1:V))) # 列名 # 正規化：sum_{v=1}^V phi_{kv} = 1 for(k in 1:K) { phi_kv[k, ] <- phi_kv[k, ] / sum(phi_kv[k, ]) }
トピック分布と同様に、sample()で$K * V$個のランダムな値を持つベクトルを作り、matrix()にnrow = K(行数)、ncol = V(列数)を指定して、$K$行$V$列の行列を作る。
$\sum_{v=1}^V \phi_{kv} = 1$の条件を満たすために、各行に対して合計値で割って正規化する。
・最尤推定
以下、for()ループ内の処理の内容。
・パラメータの初期化
# 次ステップのパラメータを初期化 new_theta_k <- rep(0, K) names(new_theta_k) <- paste0("k=", 1:K) new_phi_kv <- matrix(0, nrow = K, ncol = V, dimnames = list(paste0("k=", 1:K), paste0("v=", 1:V)))
この記事の例では、変数の要素は随時置き換えるので初期化する必要はない。しかし初期化しなくとも、ループ処理の前に受け皿としてこの操作で変数を用意しなければならないので、どうせならば参考書のアルゴリズムの説明の通りのこの位置で処理しておく。あるいは、アルゴリズムの理解のために。
・負担率の計算
# 負担率の計算 tmp_phi_k <- t(phi_kv) ^ N_dv[d, ] %>% apply(2, prod) q_dk[d, k] <- theta_k[k] * tmp_phi_k[k] / sum(theta_k * tmp_phi_k)
負担率は
の計算式で求める。
まず$\prod_{v=1}^V \phi_{kv}^{N_dv}$の部分の計算を行い、tmp_phi_kとしておく。マトリクスとベクトルの計算は、ベクトルの要素をマトリクスの1列目の要素の上から1つずつ行われる。ベクトルの要素数よりもマトリクスの要素数の方が多い場合は、ベクトルの要素を繰り返し計算に再利用される。N_dv[d, ]は$V$個のベクトルになるため、phi_kvをt()で転置して列の要素が$V$個のマトリクスとすることで、$\phi_{kv}^{N_dv}$の計算ができる。tmp_phi_kは$K$個のベクトルとなる。
続いてtmp_phi_kを、分子は特定のトピック$k$に関する要素のみ、分母は各トピックごとにtheta_kと掛け合わせる。分母については全体を足し合わせるためsum()している。
・トピック分布の更新
# トピック分布(theta_k)を更新 new_theta_k[k] <- new_theta_k[k] + q_dk[d, k]
$\theta_k$に、今求めた負担率を加えて更新する。繰り返し足して更新していき、全文書分を終えた状態で最後に正規化を行う。
・単語分布の更新
# 単語分布(phi_kv)を更新 new_phi_kv[k, v] <- new_phi_kv[k, v] + q_dk[d, k] * N_dv[d, v]
同様に、単語分布$\phi_{kv}$も各要素について負担率を加えていき更新する。全文書・語彙分を加えた状態、で最後に正規化を行う。
ただし、参考書では語彙$v$(単語を重複せずに扱う集合で出現回数を考慮する)ではなく単語$w_{dn}$(重複を許す単語集合で全ての出現回数が1回となる)を用いて$\phi_{kw_{dn}}$としているが、(そっちの方が面倒なので)語彙ごとに計算している。
単語から語彙に変換するために負担率に語彙の出現回数を掛けている。
ここまでの計算を各文書・各トピックに対して繰り返し行う。また、単語分布の更新は各語彙更新するため$V$回繰り返し処理する。
・パラメータ更新
# パラメータの更新 theta_k <- new_theta_k phi_kv <- new_phi_kv # パラメータの正規化 theta_k <- theta_k / sum(theta_k) for(k2 in 1:K) { phi_kv[k2, ] <- phi_kv[k2, ] / sum(phi_kv[k2, ]) }
次の更新に用いるために、new_theta_kとnew_phi_kvをそれぞれtheta_kとphi_kvに移す。theta_k, phi_kvは次の負担率の計算に用いる。new_theta_k, new_phi_kvは初期化して次の更新式の受け皿となる。
theta_kとphi_kvをそれぞれ$\sum_{k=1}^K \theta_k = \sum_{v=1}^V \phi_{kv} = 1$となるように正規化する。
・パラメータの推定結果の確認
・トピック分布
## トピック分布の推定結果を確認 # データフレームを作成 theta_df <- data.frame(topic = as.factor(1:K), prob = theta_k) # 描画 ggplot(data = theta_df, mapping = aes(x = topic, y = prob, fill = topic)) + # データ geom_bar(stat = "identity", position = "dodge") + # 棒グラフ labs(title = "Mixture of Unigram Models:MLE") # タイトル
どのトピックになる確率が高いのかを確認できる。
・単語分布
## 単語分布の推定結果を確認 # データフレームを作成 phi_df_wide <- cbind(as.data.frame(phi_kv), as.factor(1:K)) # データフレームをlong型に変換 colnames(phi_df_wide) <- c(1:V, "topic") # key用の行名を付与 phi_df_long <- gather(phi_df_wide, key = "word", value = "prob", -topic) # 変換 phi_df_long$term <- phi_df_long$word %>% # 文字列になるため因子に変換 as.numeric() %>% as.factor() # 描画 ggplot(data = phi_df_long, mapping = aes(x = word, y = prob, fill = word)) + # データ geom_bar(stat = "identity", position = "dodge") + # 棒グラフ facet_wrap(~ topic, labeller = label_both) + # グラフの分割 theme(legend.position = "none") + # 凡例 labs(title = "Mixture of Unigram Models:MLE") # タイトル
トピックごとの語彙の出現確率を確認できる。
・推定推移の確認
・推移の記録
# 推移の確認用 trace_theta <- as.matrix(theta_k) trace_phi <- as.matrix(phi_kv[1, ]) # 推移の確認用 trace_theta <- cbind(trace_theta, as.matrix(theta_k)) trace_phi <- cbind(trace_phi, as.matrix(phi_kv[1, ]))
推移を確認するために、パラメータを更新する度に記録しておく。単語分布については、トピック1のみ確認する。
・トピック分布
## トピック分布の推定値の推移を確認 # データフレームを作成 trace_theta_df_wide <- cbind(t(trace_theta), 1:ncol(trace_theta)) %>% # 推定回数 as.data.frame() # データフレームをlong型に変換 colnames(trace_theta_df_wide) <- c(1:K, "n") # key用の行名を付与 trace_theta_df_long <- gather(trace_theta_df_wide, key = "topic", value = "prob", -n) # 変換 # 描画 ggplot(data = trace_theta_df_long, mapping = aes(x = n, y = prob, color = topic)) + geom_line() + # 折れ線グラフ labs(title = "Mixture of Unigram Models:MLE:(theta)") # タイトル
・単語分布
## 単語分布の推定値の推移を確認 # データフレームを作成 trace_phi_df_wide <- cbind(t(trace_phi), 1:ncol(trace_phi)) %>% as.data.frame() # データフレームをlong型に変換 colnames(trace_phi_df_wide) <- c(1:V, "n") # key用の行名を付与 trace_phi_df_long <- gather(trace_phi_df_wide, key = "word", value = "prob", -n) # 変換 trace_phi_df_long$word <- trace_phi_df_long$word %>% # 文字列になるため因子に変換 as.numeric() %>% as.factor() # 描画 ggplot(data = trace_phi_df_long, mapping = aes(x = n, y = prob, color = word)) + geom_line() + # 折れ線グラフ theme(legend.position = "none") + # 凡例 labs(title = "Mixture of Unigram Models:MLE:(phi)") # タイトル
各パラメータが収束するまでの推移を確認できる。
・各トピックの出現確率の上位語の確認
### 各トピックの出現確率の上位語 ## 語彙インデックス(v) # 指定した出現回数以上の単語の行番号 num <- mecab_df %>% select(-c(TERM, POS1, POS2)) %>% apply(1, sum) >= 5 # 抽出する総出現回数を指定する v_index <- mecab_df[num, ] %>% filter(grepl("名詞|形容詞|^動詞", POS1)) %>% # 抽出する品詞を指定する filter(grepl("一般|^自立", POS2)) %>% filter(!grepl(stop_words, TERM)) %>% .[, 1] # 単語の列を抽出する # データフレームを作成 phi_df_wide2 <- cbind(as.data.frame(t(phi_kv)), v_index) %>% as.data.frame() colnames(phi_df_wide2) <- c(paste0("topic", 1:K), "word") # key用の行名を付与 # データフレームの整形 phi_df_long2 <- data.frame() for(i in 1:K) { tmp_df <- phi_df_wide2 %>% select(paste0("topic", i), word) %>% arrange(-.[, 1]) %>% # 降順に並べ替え head(20) %>% # 任意で指定した上位n語を抽出 gather(key = "topic", value = "prob", -word) # long型に変換 phi_df_long2 <- rbind(phi_df_long2, tmp_df) } # 描画 ggplot(data = phi_df_long2, mapping = aes(x = reorder(x = word, X = prob), y = prob, fill = word)) + # データ geom_bar(stat = "identity", position = "dodge") + # 棒グラフ coord_flip() + # グラフの向き facet_wrap( ~ topic, scales = "free") + # グラフの分割 theme(legend.position = "none") + # 凡例 labs(title = "Mixture of Unigram Models:VBE") # タイトル
(えーーと、ストップワードの処理などしていませんので色々お察しですが悪しからず。それとこの処理ももう少し効率よく処理できそう)
参考書籍
- 岩田具治(2015)『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社
おわりに
章が進み数式イジイジにも慣れてきて、内容もふんわり理解したつもりにもなってきたので次は実装です。半月ほどコネコネしてたらループが回るようになりました。やはりR楽しい。
参考書のアルゴリズム3.1を参考にしてるので'for()'多めです。もう少しRっぽいコードでも再現したいです。
数式でもRでも日本語でも説明できればそれなりに理解したと言ってもいいでしょうかね…
2019/08/12：加筆修正しました。
【次節の内容】