はじめに
機械学習プロフェッショナルシリーズの『トピックモデル』の勉強時に自分の理解の助けになったことや勉強会資料のまとめです。トピックモデルの各種アルゴリズムを「数式」と「プログラム」から理解することを目指します。
この記事は、2.4節「最大事後確率推定」の内容です。ユニグラムモデルにおけるMAP推定を説明して、MAP推定値を導出します。そして最後にR言語で実装します。
この記事で用いる記号類やユニグラムモデルの定義については、2.1-2：ユニグラムモデル【『トピックモデル』の勉強ノート】 - からっぽのしょこをご確認ください。
【前節の内容】
【他の節一覧】
【この節の内容】
2.4 最大事後確率推定
最大事後確率推定(MAP推定)では、事後確率$p(\boldsymbol{\phi} | \mathbf{W}, \beta)$が最大となるようにパラメータ$\boldsymbol{\phi}$を推定する。
文書データ$\mathbf{W}$を観測する前の単語分布$\boldsymbol{\phi}$の分布を、事前分布$p(\boldsymbol{\phi} | \beta)$とする(事前分布に従って$\boldsymbol{\phi}$が生成されるとする)。$\beta$は、事前分布のパラメータである。事前分布(パラメータ)のパラメータをハイパーパラメータと呼ぶ。事前分布と尤度(生成モデル)$p(\mathbf{W} | \boldsymbol{\phi})$にベイズの定理(1.4)を用いて、事後分布$p(\boldsymbol{\phi} | \mathbf{W}, \beta)$を求める。
分母の項は、分子を$\boldsymbol{\phi}$について周辺化したもので、事後分布の総和を1とするための正規化項である。
ただしこの項は$\boldsymbol{\phi}$に依存せず、最大化したい$\boldsymbol{\phi}$と無関係なのでここでは無視できる。
従って、分子を最大化するパラメータがMAP推定値になる。最尤推定のときと同様に、計算を簡単にするため対数をとって求めていく。これを式で表すと次のように書く。
MAP推定では、事前分布が必要である。事後分布が事前分布と同じ形の分布になる共役事前分布を事前分布として用いると計算が簡単になる。カテゴリ分布の共役事前分布はディリクレ分布なので、事前確率$p(\boldsymbol{\phi} | \beta)$をパラメータ$\beta, \cdots, \beta$を持つディリクレ分布とする。
尤度(2.1)と事前分布(2.7)を用いて、最適化問題(2.6)を解く。
【途中式の途中式】
- 尤度(2.1)、事前分布(2.7)より、具体的な値に置き換える。
- 対数をとり式を変形する。
$\log \left(\frac{A}{B} C\right) = \log A - \log B + \log C$より、それぞれ項を変形する。前の項は
となる。後の項についても同様に変形できる。
- $\log \Gamma(\beta V)$と$V \log \Gamma(\beta)$は、$\boldsymbol{\phi}$に影響しないので省ける。また$\log x^{a} = a \log x$である。
- $\log \phi_v$の項をまとめる。
最尤推定と同様に、ラグランジュの未定乗数法を用いて、$\sum_{v=1}^V \phi_v = 1$の制約の下で$\sum_{v=1}^V (N_v + \beta - 1) \log \phi_v$を最大化する$\boldsymbol{\phi}$を求める。
$\phi_v$に関して微分して、$\frac{\partial F}{\partial \phi_v} = 0$になる$\phi_v$を求める。
この式の両辺で$v$に関して和をとる。
【途中式の途中式】
- $\sum_{v=1}^V \phi_v = 1$なので置き換える。
- $\sum_{v=1}^V N_v$は、各単語の出現回数$N_v$の総和なので、総単語数$N$である。
これを$\phi_v = \frac{N_v + \beta - 1}{- \lambda}$に代入すると
が得られる。これがMAP推定値である。
観測がない語彙に対して確率が負にならないためには、ディリクレ分布のパラメータを$\beta \geq 1$で設定する必要がある。また$\beta = 2$のときラプラススムージングと呼ばれる。
最尤推定値$\frac{N_v}{N}$とは分子に$\beta - 1$が加えられているところが異なる。分母の$(\beta - 1) V$は、分子に$\beta - 1$が加えたとき総和が1になるのに必要な項である。これは観測が少ないときにも、事前分布(として設定した)の情報が反映されることを意味する。
・Rでやってみる
簡単な例として、文書が1つの場合について考えます。文書は1つですが、次章以降への繋がりを考えてこの文書を文書$d$と呼ぶことにします。またそれぞれの過程をイメージしやすいように、サイコロに例えながら進めます。
まずはパラメータを指定し、そのパラメータに従ってランダムに単語を生成します。生成された単語データ(文書)を用いて、MAP推定によりパラメータの推定値を求めます。
# 利用パッケージ library(tidyverse)
利用するパッケージを読み込みます。
文書$d$(文書全体)の単語数と語彙数とハイパーパラメータを指定します。
# 単語数(サイコロを振る回数)を指定 N_d <- 10 # 語彙数(サイコロの目の数)を指定 V = 6 # 単語分布のパラメータを指定 beta <- 2
単語数はサイコロを振る回数と言えます。
語彙数はサイコロの目の数と言えます。この例では、一般的なサイコロとして6を指定します。
事前分布のパラメータ$\beta \geq 1$を指定します。
真の単語分布を指定します。
# 真の単語分布(カテゴリ分布のパラメータ)を指定 phi_turth <- rep(1, V) / V phi_turth
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
sum(phi_turth)
## [1] 1
この例では(正六面体の)サイコロのように一様に設定します。一様でなくとも問題なく推定を行えますが、総和が1となるように設定する必要があります。
この値は本来知り得ない値です。ここで指定した値をMAP推定により推定することが目的になります。
指定したパラメータ類に従って、単語をランダムに生成し文書を作成します。
# 文書を生成(サイコロを振る) w_dn <- sample(x = 1:V, size = N_d, replace = TRUE, prob = phi_turth) w_dn
## [1] 1 6 6 5 5 6 2 4 3 4
サイコロの目に応じた語彙が設定してあるイメージです。
語彙ごとの出現回数をカウントします。
# 各語彙の出現回数を集計 doc_df1 <- tibble(v = w_dn) %>% group_by(v) %>% # 出目でグループ化 summarise(N_v = n()) # 出目ごとにカウント # 出ない目があったとき用の対策 doc_df2 <- tibble(v = 1:V) %>% # 1からVまでの受け皿(行)を作成 left_join(doc_df1, by = "v") %>% # 結合 mutate(N_v = replace_na(N_v, 0)) # NAを0に置換
1度も出ない語(目)があると、値が0ではなくデータフレームに含まれません。なのでその対策をします。
MAP推定を行います。
# MAP推定 map_df <- doc_df2 %>% mutate(prob = (N_v + beta - 1) / (sum(N_v) + (beta - 1) * V)) head(map_df)
## # A tibble: 6 x 3 ## v N_v prob ## <int> <dbl> <dbl> ## 1 1 1 0.125 ## 2 2 1 0.125 ## 3 3 1 0.125 ## 4 4 2 0.188 ## 5 5 2 0.188 ## 6 6 3 0.25
$\phi_v = \frac{N_v + \beta - 1}{N + (\beta - 1) V}$の計算を行い、各語彙の出現確率$\boldsymbol{\phi}$の推定値を求めます。
推定結果を可視化します。
ggplot(map_df, aes(x = v, y = prob)) + # データ geom_bar(stat = "identity", fill = "#00A968") + # 棒グラフ scale_x_continuous(breaks = 1:V, labels = 1:V) + # 軸目盛 labs(title = "Unigram Model:Maximum A Posteriori Estimation", subtitle = paste0("N=", N_d)) # ラベル
この例ではどの値も0.167としたので、上手く推定できているとは言えませんね。
ただし最尤推定値と比較すると、語彙の出現回数$N_v$に$\beta - 1$を加えることで、観測値に極端に依存しない推定値になっている。
試行回数が増やしたときの推定結果の変化を、最尤推定の場合と比較しながら確認します。
・最尤推定とMAP推定との比較
試行回数が1回・10回・100回のときの最尤推定とMAP推定との推定結果を比較する。
・コード(クリックで展開)
# 単語数を指定 N_d_vec <- c(1, 10, 100) # パラメータ生成 res_df <- tibble() for(i in seq_along(N_d_vec)){ # 文書を生成 doc_df <- tibble( v = sample(x = 1:V, size = N_d_vec[i], replace = TRUE, prob = phi_turth) # 単語を生成 ) %>% group_by(v) %>% # 語彙ごとにグループ化 summarise(N_v = n()) %>% # 各語彙の出現回数を集計 right_join( tibble( v = 1:V, N = as.factor(N_d_vec[i]) # 作図用のラベル ), by = "v" ) %>% # 結合 mutate(N_v = replace_na(N_v, 0)) # NAを0に置換 # 推定 tmp_df <- doc_df %>% mutate(ml = N_v / sum(N_v)) %>% # 最尤推定:式(2.3) mutate(map = (N_v + beta - 1) / (sum(N_v) + (beta - 1) * V)) # MAP推定 # 推定結果を結合 res_df <- rbind(res_df, tmp_df) } # wide型データフレームをlong型に変換 res_df_long <- pivot_longer( res_df, cols = c(ml, map), # 変換する列を指定 names_to = "method", # 現列名を値として格納する列名 names_ptypes = list(method = factor()), # 現列名を値とする際の型を指定 values_to = "prob" # 現セルの値を格納する列前 ) # 作図 ggplot(res_df_long, aes(x = v, y = prob, fill = method)) + # データの設定 geom_bar(stat = "identity") + # 棒グラフ facet_wrap(method ~ N, nrow = 2, labeller = label_both) + # グラフの分割 scale_x_continuous(breaks = 1:V, labels = 1:V) + # 軸目盛 scale_fill_manual(values = c("#00A968", "Orange")) + # 塗りつぶし色 labs(title = "Unigram Model") # ラベル
観測データ数が多くなるに従って、$\beta$に比べて$N_v$が大きくなるため、最尤推定値とMAP推定値が近づいていくことが確認できる。
参考書籍
- 岩田具治(2015)『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社
おわりに
前節の記事を加筆修正する際に分割したものになります。
2020.07.01：加筆修正しました。
【次節の内容】