からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

【R】2.3:ユニグラムモデルの最尤推定の実装【青トピックモデルのノート】

はじめに

 『トピックモデル』(MLPシリーズ)の勉強会資料のまとめです。各種モデルやアルゴリズムを「数式」と「プログラム」を用いて解説します。
 本の補助として読んでください。

 この記事では、カテゴリモデルに対する最尤推定をR言語でスクラッチ実装します。

【前節の内容】

www.anarchive-beta.com

【他の節の内容】

www.anarchive-beta.com

【この節の内容】

2.3 ユニグラムモデルの最尤推定の実装

 ユニグラムモデル(unigram model)に対する最尤推定(maximum likelihood estimation)を実装する。
 ユニグラムモデルの定義や記号については「2.2:ユニグラムモデルの生成モデルの導出【青トピックモデルのノート】 - からっぽのしょこ」、パラメータの計算式については「2.3:ユニグラムモデルの最尤推定の導出【青トピックモデルのノート】 - からっぽのしょこ」を参照のこと。

 利用するパッケージを読み込む。

# 利用パッケージ
library(tidyverse)
library(MCMCpack)

 この記事では、基本的に パッケージ名::関数名() の記法を使うので、パッケージの読み込みは不要である。ただし、作図コードについてはパッケージ名を省略するので、ggplot2 を読み込む必要がある。
 また、ネイティブパイプ演算子 |> を使う。magrittrパッケージのパイプ演算子 %>% に置き換えられるが、その場合は magrittr を読み込む必要がある。

文書データの簡易生成

 まずは、ユニグラムモデルの生成モデルに従って、bag-of-words表現の文書データ(語彙頻度データ)を生成する。
 トイデータの作成については「生成モデルの実装」を参照のこと。

 真の分布などの本来は得られない情報についてはオブジェクト名に true を付けて表す。

処理コード(クリックで展開)

 パラメータ類を設定して、文書データを作成する。

# 文書数を指定
D <- 10

# 語彙数を指定
V <- 6

# 単語分布のハイパーパラメータを指定
true_beta   <- 1
true_beta_v <- rep(true_beta, times = V)

# 単語分布のパラメータを生成
#true_phi_v <- rep(1/V, times = V) # 一様な値の場合
true_phi_v <- MCMCpack::rdirichlet(n = 1, alpha = true_beta_v) |> # 多様な値の場合
  as.vector()

# 各文書の単語数を生成
N_d <- sample(x = 10:20, size = D) # 下限・上限を指定

# 文書データを生成
N_dv <- matrix(NA, nrow = D, ncol = V)
for(d in 1:D) { 
  
  # 各語彙の出現回数を生成
  N_dv[d, ] <- rmultinom(n = 1, size = N_d[d], prob = true_phi_v) |> # (多項乱数)
    as.vector()
  
  # 途中経過を表示
  print(paste0("document: ", d, ", words: ", N_d[d]))
}
[1] "document: 1, words: 17"
[1] "document: 2, words: 11"
[1] "document: 3, words: 13"
[1] "document: 4, words: 14"
[1] "document: 5, words: 18"
[1] "document: 6, words: 20"
[1] "document: 7, words: 15"
[1] "document: 8, words: 19"
[1] "document: 9, words: 10"
[1] "document: 10, words: 12"

 文書数  D、語彙数  V、単語分布(カテゴリ分布)の真のハイパーパラメータ  \beta_{\mathrm{truth}} \gt 0 を指定して、真のパラメータ  \boldsymbol{\phi}_{\mathrm{truth}} を生成する。
 文書ごとに単語数  N_d を決めて、多項分布の乱数により各語彙の出現回数  N_{dv} を生成する。

 文書データを確認する。

# 観測データを確認
head(N_dv); head(N_d)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    5    7    2    1    2
[2,]    2    1    5    2    0    1
[3,]    2    0    4    6    0    1
[4,]    2    6    0    4    0    2
[5,]    1    5    1    8    2    1
[6,]    0    4    1    8    1    6
[1] 17 11 13 14 18 20


 真の単語分布をグラフで確認する。

ユニグラムモデルにおける単語分布

 横軸は語彙番号  v、縦軸は各語彙の出現(生成)確率  \phi_v を表す。

 文書ごとの語彙頻度をグラフで確認する。

ユニグラムモデルにおけるbag-of-words表現の文書データ

 横軸は語彙番号  v、縦軸は文書ごとの各語彙の出現回数  N_{dv} を表すグラフを文書ごとに縦横に並べる。

 全文書での語彙頻度をグラフで確認する。

全文書での頻度データ

 横軸は語彙番号  v、縦軸は全文書での各語彙の出現回数  N_v を表す。

 文書集合  \mathbf{W} から得られる語彙頻度データ  (N_{11}, \cdots, N_{DV}) のみを用いてパラメータ類を推定する。

パラメータの推定

 次は、最尤推定によりパラメータを推定する。

 文書集合に関する値を設定する。

# 文書数を取得
D <- nrow(N_dv)

# 語彙数を取得
V <- ncol(N_dv)

# 全文書の単語数を取得
N <- sum(N_dv)

# 各文書の単語数を取得
N_d <- rowSums(N_dv)

# 各語彙の出現回数を取得
N_v <- colSums(N_dv)
D; V; N; N_d; N_v
[1] 10
[1] 6
[1] 149
 [1] 17 11 13 14 18 20 15 19 10 12
[1] 10 35 32 47  7 18

 観測データ N_dv から文書数  D、語彙数  V、単語数  N = \sum_{d=1}^D \sum_{v=1}^V N_{dv}、各文書の単語数  (N_1, \cdots, N_D) N_d = \sum_{v=1}^V N_{dv}、各語彙の単語数  (N_1, \cdots, N_V) N_v = \sum_{d=1}^D N_{dv} を取得する。

 単語分布のパラメータの最尤解を計算する。

# 単語分布のパラメータを計算:式(2.4)
phi_v <- N_v / N
phi_v; sum(phi_v)
[1] 0.06711409 0.23489933 0.21476510 0.31543624 0.04697987 0.12080537
[1] 1

 単語分布のパラメータ  \boldsymbol{\phi} の各項を、次の式で計算する。

 \displaystyle
\phi_v
    = \frac{N_v}{N}
\tag{2.4}


 以上で、最尤推定によりパラメータを推定した。

推定結果の可視化

 続いて、推定したパラメータや分布のグラフを作成する。

分布の図

 単語分布の描画用のデータフレームを作成する。

# 推定した単語分布を格納
phi_df <- tibble::tibble(
  v    = 1:V, 
  prob = phi_v
)
phi_df
# A tibble: 6 × 2
      v   prob
  <int>  <dbl>
1     1 0.0671
2     2 0.235 
3     3 0.215 
4     4 0.315 
5     5 0.0470
6     6 0.121 

 語彙番号  v とパラメータ  \phi_v を格納する。

 単語分布のグラフを作成する。

# 真の単語分布を格納
true_phi_df <- tibble::tibble(
  v    = 1:V, 
  prob = true_phi_v
)

# ラベル用の文字列を作成
param_label <- paste0(
  "N == ", N
)

# 凡例用の設定を格納
linetype_lt <- list(
  fill  = c("gray", NA), 
  color = c(NA, "red")
)

# 単語分布を作図
ggplot() + 
  geom_bar(data = phi_df, 
           mapping = aes(x = v, y = prob, fill = factor(v), linetype = "estimated"), 
           stat = "identity") + # 推定した分布
  geom_bar(data = true_phi_df, 
           mapping = aes(x = v, y = prob, linetype = "true"), 
           stat = "identity", show.legend = FALSE, 
           fill = NA, color = "red", linewidth = 1) + # 真の分布
  scale_linetype_manual(breaks = c("estimated", "true"), 
                        values = c("solid", "dashed"), 
                        labels = c("estimated", "true"), 
                        name = "distribution") + # 凡例の表示用
  guides(fill = "none", 
         linetype = guide_legend(override.aes = linetype_lt)) + # 凡例の体裁
  labs(title = "word distribution (maximum likelihood estimation)", 
       subtitle = parse(text = param_label), 
       x = expression(vocabulary ~ (v)), 
       y = expression(probability ~ (phi[v])))

ユニグラムモデルにおける単語分布

 横軸は語彙番号  v、縦軸は各語彙の出現確率  \phi_v を表す。
 真の分布を赤色の破線、推定した分布を塗りつぶしで示す。

推定への影響の可視化

 最後は、パラメータなどの推定結果への影響をグラフで確認する。
 アニメーションなどの作図コードは「map_estimation.R · anemptyarchive/Topic-Models-AO · GitHub」を参照のこと。

データ数の影響

 データ数の影響をアニメーションで確認する。

 新たに観測された語彙をオレンジ色の点で示す。

  N が大きくなるほど  \beta の影響が小さくなり、真の分布の形状に近付いていくのを確認できる。

 KLダイバージェンスの推移をグラフで確認する。

KLダイバージェンスの推移

 横軸はデータ数(全文書の単語数)  N、縦軸は真のパラメータ  \boldsymbol{\phi}_{\mathrm{truth}} とデータ数  N のときの最尤推定値  \boldsymbol{\phi}_{\mathrm{ML}} のKLダイバージェンスを表す。
 ただし、観測されないカテゴリ(語彙  v の出現回数が  N_v = 0 のとき)の推定確率が  \phi_v = 0 になり、 - \log 0 = \infty なので、データ数が少ないとKLダイバージェンスを計算できないことがある。また、ggplot2の作図では、Inf は描画領域外の端を意味する。

 KLダイバージェンスについては「1.1.8-10:カルバック・ライブラー・ダイバージェンスとイェンゼンの不等式【『トピックモデル』の勉強ノート】 - からっぽのしょこ」を参照のこと。

 この記事では、ユニグラムモデルに対する最尤推定を実装してパラメータを推定した。次の記事では、MAP推定によるパラメータの計算式を導出する。

参考書籍

  • 岩田具治(2015)『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社

おわりに

  • 2024.05.20:加筆修正の際に「ユニグラムモデルの最尤推定の導出」から記事を分割しました。

 推定処理自体は1行で済むのですが、図を充実させて1記事分の内容になったと思います。Python版も書いたときの構成上の都合もあります。
 修正前は本に合わせて、1つのサイコロを連続して振った結果のイメージで解説しました。修正後は次章以降と対応するように、複数の文書データのイメージで解説しました。

【次節の内容】

  • 数式読解編

 ユニグラムモデルに対するMAP推定を数式で確認します。

www.anarchive-beta.com


  • スクラッチ実装編

 ユニグラムモデルに対するMAP推定をプログラムで確認します。

www.anarchive-beta.com