はじめに

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

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

【前節の内容】

www.anarchive-beta.com

【他の節の内容】

www.anarchive-beta.com

【この節の内容】

はじめに
2.3 ユニグラムモデルの最尤推定の実装
参考書籍
おわりに

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$ を表す。

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

　横軸は語彙番号 $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ダイバージェンスの推移をグラフで確認する。

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