からっぽのしょこ

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

2.3:ユニグラムモデルの最尤推定【『トピックモデル』の勉強ノート】

はじめに

 機械学習プロフェッショナルシリーズの『トピックモデル』の勉強時に自分の理解の助けになったことや勉強会資料のまとめです。トピックモデルの各種アルゴリズムを「数式」と「プログラム」から理解することを目指します。

 この記事は、2.3節「最尤推定」の内容です。ユニグラムモデルにおける最尤推定を説明して、最尤推定値を導出します。そして最後にR言語で実装します。

 この記事で用いる記号類やユニグラムモデルの定義については、前節の記事をご確認ください。

【前節の内容】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】


2.3 最尤推定

 最尤推定では、尤度$p(\mathbf{W} | \boldsymbol{\phi})$を最大化するようにパラメータ$\boldsymbol{\phi}$を推定する。

 ユニグラムモデル(生成モデル)のパラメータを$\boldsymbol{\phi}$、データを$\mathbf{W}$としたとき、尤度は条件付き確率$p(\mathbf{W} | \boldsymbol{\phi})$を$\boldsymbol{\phi}$の関数とみたものになる。

 計算を簡単にするため、対数をとった尤度である対数尤度を最大化する。対数は単調増加するので、対数をとっても尤度が最大化されるパラメータの値は変わらない。

$$ \mathop{\mathrm{argmax}}\limits_{\boldsymbol{\phi}} \log p(\mathbf{W} | \boldsymbol{\phi}) $$

 $\mathop{\mathrm{argmax}}\limits_{\boldsymbol{\phi}} f(x)$は、関数$f(x)$を最大化させる$\boldsymbol{\phi}$を求めることである。

 ラグランジュの未定乗数法を用いて、$\sum_{v=1}^V \phi_v = 1$の制約下で$\log p(\mathbf{W} | \boldsymbol{\phi})$が最大となるパラメータ$\boldsymbol{\phi}$を求める。

$$ \begin{align} F &= \log p(\mathbf{W}|\boldsymbol{\phi}) + \lambda \left(\sum_{v=1}^V \phi_v - 1 \right) \tag{2.2}\\ &= \log \left(\prod_{v=1}^V \phi_v^{N_v} \right) + \lambda \left(\sum_{v=1}^V \phi_v - 1 \right) \\ &= \sum_{v=1}^V N_v \log \phi_v + \lambda \left(\sum_{v=1}^V \phi_v - 1 \right) \end{align} $$

【途中式の途中式】

  1. ラグランジュの未定乗数$\lambda$を用いて、関数$F = f(x) + \lambda g(x)$を立てる。
  2. ユニグラムモデルの定義(2.1)より、具体的な値に置き換える。
  3. $\log (A * B) = \log A + \log B$、$\log x^a = a * \log x$であることから、式を変形する。
$$ \begin{aligned} \log \left( \prod_{v=1}^V \phi_v^{N_v} \right) &= \log ( \phi_1^{N_1} * \cdots * \phi_v^{N_v} * \cdots * \phi_V^{N_V} ) \\ &= \log \phi_1^{N_1} + \cdots + \log \phi_v^{N_v} + \cdots + \log \phi_V^{N_V} \\ &= \sum_{v=1}^V \log \phi_v^{N_v} \\ &= \sum_{v=1}^V N_v \log \phi_v \end{aligned} $$


 この式を$\phi_v$に関して微分して、$\frac{\partial F}{\partial \phi_v} = 0$となる$\phi_v$を求める。

$$ \begin{align} \frac{\partial F}{\partial \phi_v} = \frac{N_v}{\phi_v} + \lambda &= 0 \\ \phi_v &= -\frac{N_v}{\lambda} \tag{2.3} \end{align} $$

【途中式の途中式】

  • 自然対数の微分は$(\log x)' = \frac{1}{x}$である。
  • また$\phi_v$に関しての微分なので、それ以外の$\phi_1$から$\phi_V$については定数として扱うため0になる(これを偏微分と言う)。

 従って前の項は

$$ \begin{aligned} \sum_{v=1}^V N_v \log \phi_v &= N_v \log \phi_1 + \cdots + N_v \log \phi_v + \cdots + N_v \log \phi_V \\ \left( \sum_{v=1}^V N_v \log \phi_v \right)' &= 0 + \cdots + \frac{N_v}{\phi_v} + \cdots + 0 \\ &= \frac{N_v}{\phi_v} \end{aligned} $$

になる。

 同様に後ろの項は

$$ \begin{aligned} \lambda \left( \sum_{v=1}^V \phi_v - 1 \right) &= (\lambda \phi_1 - \lambda) + \cdots + (\lambda \phi_v - \lambda) + \cdots + (\lambda \phi_V - \lambda) \\ \left\{ \lambda \left( \sum_{v=1}^V \phi_v - 1 \right) \right\}' &= 0 + \cdots + \lambda * 1 + \cdots + 0 \\ &= \lambda \end{aligned} $$

になる。



 この式の両辺で$v$に関して和をとると

$$ \begin{aligned} \sum_{v=1}^V \phi_v &= - \sum_{v=1}^V \frac{N_v}{\lambda} \\ 1 &= - \frac{N}{\lambda} \\ \lambda &= -N \end{aligned} $$

【途中式の途中式】

  • 単語分布$\boldsymbol{\phi}$の定義より、$\sum_{v=1}^V \phi_v = 1$である。
  • $\sum_{v=1}^V N_v$は各語彙の出現回数$N_v$の総和なので、総単語数$N$になる。


となる。これを式(2.3)に代入すると

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

が得られる。

 これは語彙$v$の出現回数$N_v$を総単語数$N$で割ったものである。つまり文書全体で語彙$v$が現れた割合が、$\phi_v$の最尤推定値となる。

・Rでやってみる

 簡単な例として、文書が1つの場合について考えます。文書は1つですが、次章以降への繋がりを考えてこの文書を文書$d$と呼ぶことにします。またそれぞれの過程をイメージしやすいように、サイコロに例えながら進めます。

 まずはパラメータを指定し、そのパラメータに従ってランダムに単語を生成します。生成された単語データ(文書)を用いて、最尤推定によりパラメータの推定値を求めます。

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

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

 文書$d$(文書全体)の単語数と語彙数を指定します

# 単語数(サイコロを振る回数)を指定
N_d <- 10

# 語彙数(サイコロの目の数)を指定
V <- 6

 単語数はサイコロを振る回数と言えます。

 語彙数はサイコロの目の数と言えます。この例では、一般的なサイコロとして6を指定します。

 真の単語分布を指定します。

# 真の単語分布(カテゴリ分布のパラメータ)を指定
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となるように値を設定する必要があります。

 この値は本来知り得ない値です。ここで指定した値を最尤法により推定することが目的になります。

 指定したパラメータ類に従って、単語をランダムに生成し文書を作成します。

# 文書を生成(サイコロを振る)
w_dn <- sample(x = 1:V, size = N_d, replace = TRUE, prob = phi_turth)
w_dn
##  [1] 3 4 5 2 1 3 2 1 2 2

 サイコロの目に応じた語彙が設定してあるイメージです。

 語彙ごとの出現回数をカウントします。

# 各語彙の出現回数を集計
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ではなくデータフレームに含まれません。なのでその対策をします。

 最尤推定を行います。

# 最尤推定
ml_df <- doc_df2 %>% 
  mutate(prob = N_v / sum(N_v)) # 式(4.2)
head(ml_df)
## # A tibble: 6 x 3
##       v   N_v  prob
##   <int> <dbl> <dbl>
## 1     1     2   0.2
## 2     2     4   0.4
## 3     3     2   0.2
## 4     4     1   0.1
## 5     5     1   0.1
## 6     6     0   0

 式(2.4)の計算を行い、各語彙の出現確率$\boldsymbol{\phi}$の推定値を求めます。

 推定結果を可視化します。

# 作図
ggplot(ml_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 Likelihood Estimation", 
       subtitle = paste0("N=", N_d)) # ラベル

最尤推定値:$N=10$

 この例ではどの値も0.167としたので、上手く推定できているとは言えませんね。

 単語数(試行回数)を増やしてみましょう。

・単語数が100の場合

最尤推定値:$N=100$

 先ほどよりは、真の値に近い値が推定されています。

・単語数が1000の場合

最尤推定値:$N=1000$

 かなり良い推定値となりました。

 最尤推定では、試行回数が少ないと過学習のおそれがあることが分かる。例えば、観測されなかった要素については確率を0と推定する。あるいは、サイコロを1回振って4が出たという観測データから推定すると、必ず4の目が出るサイコロであると推定することになる。このように、試行回数によっては事実と異なる推定結果となってしまう。
 そこで次節の最大事後確率推定を用いることで、この問題を緩和できる。

参考書籍

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

おわりに

2019.07.14:加筆修正の上、記事を一部分割しました。

2020.07.01:加筆修正の上、記事を一部分割しました。

【次節の内容】

www.anarchive-beta.com