元バイオ系

元バイオウェット系がデータサイエンスやらを勉強していくブログ。 基本自分用のまとめ。

EMアルゴリズムを理解したい【最尤法から】

よく忘れるので一旦ちゃんとEMアルゴリズム*1についてまとめておきます。
最尤法について振り返ってからEMアルゴリズムの説明に入っていきます。
あまり計算に自信はありませんので、間違いを見つけたら教えていただけると幸いです。

最尤法

観測事象の確率密度関数から尤度を構成し、その確率密度関数のパラメータの推測を尤度の最大化によって行う方法です。
さっそく説明をしていきたいところですがその前に、尤度について振り返りましょう。

尤度についてはこちらをどうぞ

hotoke-x.hatenablog.com

それでは最尤法についてですが、実は一度正規分布を例に一度記事を書いています。

hotoke-x.hatenablog.com

ここでは特定の分布を仮定せず、一般的な最尤法について書いていきます。


 \thetaを母数、 x を確率変数 \displaystyle Xの実現値、 p({\boldsymbol x}|{\boldsymbol \theta})を母集団が従う確率密度関数とします。
この時、 n個の観測値 x_1, x_2, \ldots , x_nの同時確率関数は  \displaystyle {\boldsymbol x} = \left(x_1, \ldots, x_n \right)^Tとして、

$$ \begin{align} p({\boldsymbol x} | {\boldsymbol \theta}) = \prod_{i=1}^{n} p(x_i|{\boldsymbol \theta}) \end{align} $$

で与えられます。これを、同時確率分布ではなく、母数 {\boldsymbol \theta}の関数として考えると

$$ \begin{align} L_n ({\boldsymbol \theta} |{\boldsymbol x}) = \prod_{i=1}^{n} p(x_i | {\boldsymbol \theta}) \end{align} $$

と書け、これを尤度関数と呼びます。

最尤法は、式(2)の尤度を最大化する \displaystyle {\boldsymbol \theta}を推定する方法です。この時の \displaystyle {\boldsymbol \theta}最尤推定値と呼び、 \displaystyle \hat {\boldsymbol \theta}と書くのが慣習的となっています。また、 \displaystyle {\boldsymbol x}の関数、 \displaystyle \hat {\boldsymbol \theta} ({\boldsymbol x})最尤推定量と呼びます。忘れがちなのでちゃんと区別しましょう。

現実問題では、最尤推定量を陽に得ることが難しいことが多いため、対数尤度

$$ \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} $$

を数値的に最適化して \displaystyle \hat {\boldsymbol \theta} ({\boldsymbol x})を求めることが多く、その最適化法の一つがEMアルゴリズムがあります。最適化手法はもちろん他にも色々あります。その他についてはこちら。

hotoke-x.hatenablog.com

EMアルゴリズム

実現値が得られない観測に対する確率変数を \displaystyle Yとすると、観測データの確率密度関数

$$ \begin{align} p ({\boldsymbol x}|{\boldsymbol \theta}) = \int p ({\boldsymbol x},{\boldsymbol y}|{\boldsymbol \theta}) \mathrm{d}{\boldsymbol y} \end{align} $$

と書けます。「観測できない点は考えられる点を全部調べる」という気持ちが積分に表れていますね。

なぜ観測できていない点をわざわざ式に埋め込むかというと、不完全データの対数尤度関数より完全データの対数尤度関数の方が扱いやすいからです。

式の上では完全データとはいえ、実際には観測できてない確率変数が混ざっているのでちょっと扱いを工夫します。それがEMアルゴリズムです。

「なぜそれで良いのか」という細かいことは置いといて、まずはEMアルゴリズムの枠組みをおさらいしましょう。

  1. パラメータの初期値 \displaystyle {\boldsymbol \theta}^{(0)}を適当な値に設定(スケールくらいは当たりをつけて置くと良い)。
  2. 観測値 \displaystyle {\boldsymbol x}とパラメータの推定値 \displaystyle {\boldsymbol \theta}^{(0)}から、観測できない \displaystyle Yについて条件付き期待値を算出。これを疑似的な \displaystyle Yの観測値とする。疑似観測 \displaystyle Yから疑似完全データを構成(Eステップ)
  3. 疑似完全データに基づいて対数尤度 \displaystyle l ({\boldsymbol \theta}, {\boldsymbol x})を最大化するパラメータ \displaystyle {\boldsymbol \theta}^{(1)}を推定(Mステップ)
  4.  \displaystyle {\boldsymbol \theta}^{(1)}を新たなパラメータ初期値としてEステップ、Mステップを繰り返す。

さて、いろいろ気になる点があります。

  1. 疑似的な観測値を構成するのに条件付き期待値を使ってよいのか
  2. 疑似的な完全データの尤度を最大化することが、観測データの尤度を最大にしたことになるのか
  3. 疑似的な完全データによる推測と、完全な観測による推測とで乖離はないのか(バイアスがあるのではないかということ)

1つ目は今回の記事では解説しきれそうにないので後日追記ということにして、2つ目の疑問から確認していきます。

疑似的な完全データの尤度を最大化することが、観測データの尤度を最大にしたことになるのか

観測 \displaystyle {\boldsymbol X}={\boldsymbol x}が与えられたときの尤度 \displaystyle l \left({\boldsymbol \theta}, {\boldsymbol x} \right)

$$ \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} $$

となります。しかし困りました。観測できない \displaystyle {\boldsymbol y}が式に含まれています。仕方がないので適当な \displaystyle {\boldsymbol \theta}={\boldsymbol \theta^{(k)}}として完全データの条件付き期待値を計算しましょう(つまり観測していない部分を周辺化するということ)。適当に \displaystyle {\boldsymbol \theta^{(k)}}を使っても問題ない理由もこの後わかります。

$$ \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} $$

これは \displaystyle \boldsymbol \theta^{(k)} \displaystyle \boldsymbol \thetaの関数になっているので、新たに関数 \displaystyle Q, Hを導入して

\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}

と書くことができます。また、 \displaystyle Q, Hは、 \displaystyle \boldsymbol \theta^{(k)}によって定まる関数  p\left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta^{(k)}} \right) を引数にとる汎関数であるという見方もできます。

さて、ここで \displaystyle Qの式を改めて眺めてみます。 \displaystyle Q

$$ \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} $$

です。これは、事前分布に基づいて完全データの対数尤度の期待値を計算しているだけです。完全データの条件付き期待値を用いた対数尤度の最適化ができれば、実際に完全データが得られた場合の対数尤度の値へ収束することがわかります。初期値に適当な \displaystyle {\boldsymbol \theta^{(k)}}を使って問題ないのは、「完全データの条件付き期待値の対数尤度が最適化可能である」という前提からなのでした(実際は局所解へ陥る可能性があることに注意)。

いやいや関数 \displaystyle Hは放っておいていいのかよと思われるかもしれませんが、この後説明します。

疑似的な完全データによる推測と、完全な観測による推測とで乖離はないのか(バイアスがあるのではないかということ)

当然バイアスはあります。実は、疑似データを用いたことによる現実との乖離度合いを表しているのが \displaystyle Hになります。式を眺めてみましょう。

$$ \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)「計算統計学の方法」朝倉書店

qiita.com

*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