EMアルゴリズムを理解したい【最尤法から】
よく忘れるので一旦ちゃんとEMアルゴリズム*1についてまとめておきます。
最尤法について振り返ってからEMアルゴリズムの説明に入っていきます。
あまり計算に自信はありませんので、間違いを見つけたら教えていただけると幸いです。
最尤法
観測事象の確率密度関数から尤度を構成し、その確率密度関数のパラメータの推測を尤度の最大化によって行う方法です。
さっそく説明をしていきたいところですがその前に、尤度について振り返りましょう。
尤度についてはこちらをどうぞ
それでは最尤法についてですが、実は一度正規分布を例に一度記事を書いています。
ここでは特定の分布を仮定せず、一般的な最尤法について書いていきます。
を母数、 を確率変数の実現値、を母集団が従う確率密度関数とします。
この時、個の観測値の同時確率関数は として、
$$ \begin{align} p({\boldsymbol x} | {\boldsymbol \theta}) = \prod_{i=1}^{n} p(x_i|{\boldsymbol \theta}) \end{align} $$
で与えられます。これを、同時確率分布ではなく、母数の関数として考えると
$$ \begin{align} L_n ({\boldsymbol \theta} |{\boldsymbol x}) = \prod_{i=1}^{n} p(x_i | {\boldsymbol \theta}) \end{align} $$
と書け、これを尤度関数と呼びます。
最尤法は、式(2)の尤度を最大化するを推定する方法です。この時のを最尤推定値と呼び、と書くのが慣習的となっています。また、の関数、を最尤推定量と呼びます。忘れがちなのでちゃんと区別しましょう。
現実問題では、最尤推定量を陽に得ることが難しいことが多いため、対数尤度
$$ \begin{align} l \left({\boldsymbol \theta}, {\boldsymbol x} \right) = \log L \left({\boldsymbol \theta}, {\boldsymbol x} \right) = \sum_{i=1}^{n} \log p \left( x_i | {\boldsymbol \theta} \right) \end{align} $$
を数値的に最適化してを求めることが多く、その最適化法の一つがEMアルゴリズムがあります。最適化手法はもちろん他にも色々あります。その他についてはこちら。
EMアルゴリズム
実現値が得られない観測に対する確率変数をとすると、観測データの確率密度関数は
$$ \begin{align} p ({\boldsymbol x}|{\boldsymbol \theta}) = \int p ({\boldsymbol x},{\boldsymbol y}|{\boldsymbol \theta}) \mathrm{d}{\boldsymbol y} \end{align} $$
と書けます。「観測できない点は考えられる点を全部調べる」という気持ちが積分に表れていますね。
なぜ観測できていない点をわざわざ式に埋め込むかというと、不完全データの対数尤度関数より完全データの対数尤度関数の方が扱いやすいからです。
式の上では完全データとはいえ、実際には観測できてない確率変数が混ざっているのでちょっと扱いを工夫します。それがEMアルゴリズムです。
「なぜそれで良いのか」という細かいことは置いといて、まずはEMアルゴリズムの枠組みをおさらいしましょう。
- パラメータの初期値を適当な値に設定(スケールくらいは当たりをつけて置くと良い)。
- 観測値とパラメータの推定値から、観測できないについて条件付き期待値を算出。これを疑似的なの観測値とする。疑似観測から疑似完全データを構成(Eステップ)
- 疑似完全データに基づいて対数尤度を最大化するパラメータを推定(Mステップ)
- を新たなパラメータ初期値としてEステップ、Mステップを繰り返す。
さて、いろいろ気になる点があります。
- 疑似的な観測値を構成するのに条件付き期待値を使ってよいのか
- 疑似的な完全データの尤度を最大化することが、観測データの尤度を最大にしたことになるのか
- 疑似的な完全データによる推測と、完全な観測による推測とで乖離はないのか(バイアスがあるのではないかということ)
1つ目は今回の記事では解説しきれそうにないので後日追記ということにして、2つ目の疑問から確認していきます。
疑似的な完全データの尤度を最大化することが、観測データの尤度を最大にしたことになるのか
観測が与えられたときの尤度は
$$ \begin{align} l \left({\boldsymbol \theta}, {\boldsymbol x} \right) &= \log p \left({\boldsymbol x}|{\boldsymbol \theta} \right) \\ &= \log p \left({\boldsymbol x, \boldsymbol y}|{\boldsymbol \theta} \right) - \log p \left({\boldsymbol y}|{\boldsymbol x, \boldsymbol \theta} \right) \end{align} $$
となります。しかし困りました。観測できないが式に含まれています。仕方がないので適当なとして完全データの条件付き期待値を計算しましょう(つまり観測していない部分を周辺化するということ)。適当にを使っても問題ない理由もこの後わかります。
$$ \begin{align} l \left({\boldsymbol \theta}, {\boldsymbol x} \right) &= \int p \left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta^{(k)}} \right) \left\{ \log p \left({\boldsymbol x, \boldsymbol y}|{\boldsymbol \theta} \right) - \log p \left({\boldsymbol y}|{\boldsymbol x, \boldsymbol \theta} \right) \right\} \mathrm{d}{\boldsymbol y} \\ &= \mathrm{E}_{\boldsymbol \theta^{(k)} } \left[ \log p \left({\boldsymbol x, \boldsymbol y}|{\boldsymbol \theta} \right) \right] - \mathrm{E}_{\boldsymbol \theta^{(k)} } \left[ \log p \left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta} \right) \right] \end{align} $$
これはとの関数になっているので、新たに関数を導入して
\begin{align} l \left({\boldsymbol \theta}, {\boldsymbol x} \right) &= Q \left({\boldsymbol \theta}, {\boldsymbol \theta^{(k)}} \right) - H \left( {\boldsymbol \theta}, {\boldsymbol \theta^{(k)}} \right) \end{align}
と書くことができます。また、は、によって定まる関数を引数にとる汎関数であるという見方もできます。
さて、ここでの式を改めて眺めてみます。は
$$ \begin{align} Q \left({\boldsymbol \theta}, {\boldsymbol \theta^{(k)}} \right) &= \int p \left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta^{(k)}} \right) \log p \left({\boldsymbol x, \boldsymbol y}|{\boldsymbol \theta} \right) \mathrm{d}{\boldsymbol y} \end{align} $$
です。これは、事前分布に基づいて完全データの対数尤度の期待値を計算しているだけです。完全データの条件付き期待値を用いた対数尤度の最適化ができれば、実際に完全データが得られた場合の対数尤度の値へ収束することがわかります。初期値に適当なを使って問題ないのは、「完全データの条件付き期待値の対数尤度が最適化可能である」という前提からなのでした(実際は局所解へ陥る可能性があることに注意)。
いやいや関数は放っておいていいのかよと思われるかもしれませんが、この後説明します。
疑似的な完全データによる推測と、完全な観測による推測とで乖離はないのか(バイアスがあるのではないかということ)
当然バイアスはあります。実は、疑似データを用いたことによる現実との乖離度合いを表しているのがになります。式を眺めてみましょう。
$$ \begin{align} H \left( {\boldsymbol \theta}, {\boldsymbol \theta^{(k)}} \right) &= - \int p \left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta^{(k)}} \right) \log p \left({\boldsymbol y}|{\boldsymbol x, \boldsymbol \theta} \right) \mathrm{d}{\boldsymbol y} \\ \end{align} $$
ある程度勉強してきた方ならピンとくるのではないでしょうか。クロスエントロピーの式になっていますね。つまり、理想的な最適解が求まったとすればこの値は自己エントロピーに一致し最小となります。
対数尤度を最大にするのにこの項は最小になってよいのかと思われる方がいるかもしれませんが、その分、完全データの条件付き期待値の対数尤度が大きくなります。そもそも「疑似データを用いたことによる現実との乖離度合いを表している」だけなのです。
疑似的な観測値を構成するのに条件付き期待値を使ってよいのか(pending)
まず、これは個人の見解であり正確な解釈ではない可能性があることをご理解下さい。
実はEMアルゴリズムは暗に指数分布族と呼ばれる確率密度関数を使うことを仮定しているのだと思います。
指数分布で考えると、対数尤度をとるときにexponentialをきれいに打ち消せるのです。
すると全体が見えてきて、十分統計量をその条件付き期待値で置き換えて考える問題に帰着します。
今回は実用重視の記事にしたいので、数理統計の細かいところは次の機会もしくは追記という形にしたいと思います。
参考
小西 貞則、越智 義道、大森 裕浩(2008)「計算統計学の方法」朝倉書店
*1:Dempster, A. P., Laird, N. M. and Rubin, D. B. (1977), Maximum Likelihood from Incomplete Data Via the EM Algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39: 1-22. doi:10.1111/j.2517-6161.1977.tb01600.x