EMアルゴリズム

この記事では最尤推定およびMAP推定を解くための手法「EMアルゴリズム」について解説する。
(一連の記事では、紛らわしくはあるが確率分布を表す\(p\)という記号を異なる確率分布にも使い回して表記を簡素化している。)

関連記事

目次

最尤推定・MAP推定

確率変数\(x\)がパラメーター\(\theta\)の確率分布\(p(x|\theta)\)に従うことがわかっているとする。
\(p\)に従って独立に発生した幾つかのデータ\(X^{D} = \{x^{D}_n\}_{n=1}^{N}\)が既に観測されているとき、次に発生する\(x\)の分布をより正確に推測したい。
そのための基礎的な方針・問題設定が、最尤推定とMAP推定である。

どちらも正規分布や指数分布などの指数型の簡素な分布を扱う場合は解析的に解くこともできるが、より複雑な分布の場合は様々なアルゴリズムを用いなければ解くことが困難である。

最尤推定

最尤推定(Maximum likelihood estimation, MLE)では、尤度を最大化することを目指す。

前提
  • \(x\)および\(x^D_n\)が従う確率分布\(p(x|\theta)\)が関数として既知。
  • 分布\(p\)に従って独立に生成されたデータの集合\(X^D\)が観測されている。
目的

観測データに最も適合する\(\theta\)をただ一つ求める。

パラメーターが\(\theta\)であるときにデータ\(X^D\)が生成される確率\(L(\theta|X^D)\)を尤度と呼ぶ。

\begin{align} L(\theta|X^D) = \prod_{n=1}^N p(x^D_n|\theta) \end{align}

また、尤度の対数\(\ell(\theta|X^D)\)を対数尤度と呼ぶ。

\begin{align} \ell(\theta|X^D) &= \log\left(L(\theta|X^D)\right) \\ &= \sum_n \log p(x^D_n|\theta) \end{align}

最尤推定では尤度を最大化するパラメーターを求めることを目的とする。

\begin{align} \theta_{ML} &= \underset{\theta}{\rm argmax} L(\theta|X^D) \\ &= \underset{\theta}{\rm argmax} \ell(\theta|X^D) \\ \end{align}

MAP推定

MAP推定(Maximum a posteriori estimation)では、\(\theta\)の事前確率が既知であるという前提を追加して、\(\theta\)の事後確率を最大化することを目指す。

前提
  • \(x\)および\(x^D_n\)が従う確率分布\(p(x|\theta)\)が関数として既知。
  • \(\theta\)の(事前)確率分布\(p(\theta)\)が関数として既知。
  • 分布\(p\)に従って独立に生成されたデータの集合\(X^D\)が観測されている。
目的

観測データに最も適合する\(\theta\)をただ一つ求める。

Bayesの定理より、\(\theta\)の事後確率\(p(\theta|X^D)\)は以下のようになる。

\begin{align} p(\theta|X^D) &= \frac{p(\theta, X^D)}{p(X^D)} \\ &= \frac{p(X^D|\theta)p(\theta)}{p(X^D)} \\ \end{align}

MAP推定では\(\theta\)の事後確率を最大化するパラメーターを求めることを目的とする。

\begin{align} \theta_{MAP} &= \underset{\theta}{\rm argmax} p(\theta|X^D) \\ &= \underset{\theta}{\rm argmax} \frac{p(X^D|\theta)p(\theta)}{p(X^D)} \\ &= \underset{\theta}{\rm argmax} p(X^D|\theta)p(\theta) \\ &= \underset{\theta}{\rm argmax} \left( \ell(\theta|X^D) + \log p(\theta) \right) \\ \end{align}

EMアルゴリズム

EMアルゴリズム(expectation–maximization algorithm)は、1977年にDempster, Laird, Rubinらによって考案された最尤推定・MAP推定の反復アルゴリズム。

前提
  • \(p(x, z|\theta), p(z|x, \theta)\)が関数として既知。
  • (MAP推定の場合) \(p(\theta)\)も関数として既知。
  • 分布\(p\)に従って生成された値\((x, z)\)の内、\(x\)だけが観測されていて独立なデータの集合\(X^D\)として値が得られている。
目的

観測データに最も適合する\(\theta\)をただ一つ求める。

\(z\)のように直接観測することはできないが、前提として定めた確率モデルから推測することができる確率変数のことを潜在変数と言う。

アルゴリズムの導出

Jensenの不等式から、対数尤度に対して以下の不等式が成立する。

\begin{align} \ell(\theta|X^D) &= \log p(X^D|\theta) \\ &= \log\left(\int p(X^D, Z|\theta) dZ\right) \\ &= \log\left(\int q(Z) \frac{p(X^D, Z|\theta)}{q(Z)} dZ\right) \\ &\ge \int q(Z) \log\left( \frac{p(X^D, Z|\theta)}{q(Z)}\right) dZ \\ &=: F(\theta, q|X^D) \end{align}

ここで、\(q(Z)\)は\(Z\)の任意の確率分布。\(F(\theta, q|X^D)\)を変分下限(Evidence lower bound, ELBO)と呼ぶ。
変分下限を再び調べると、以下のようになる。

\begin{align} F(\theta, q|X^D) &= \int q(Z) \log\left( \frac{p(X^D, Z|\theta)}{q(Z)}\right) dZ \\ &= \int q(Z) \log\left( \frac{p(Z|X^D, \theta)p(X^D|\theta)}{q(Z)}\right) dZ \\ &= \int q(Z) \log p(X^D|\theta) dZ + \int q(Z) \log\left( \frac{p(Z|X^D, \theta)}{q(Z)}\right) dZ \\ &= p(X^D|\theta) - \int q(Z) \log\left( \frac{q(Z)}{p(Z|X^D, \theta)}\right) dZ \\ &= \ell(\theta|X^D) - D_{KL}\left(q\|p(\cdot|X^D, \theta)\right) \\ \end{align}

KL情報量の最小値は0なので、尤度の最大値と変分下限の最大値は一致する。

\begin{align} \max_{\theta, q} F(\theta, q|X^D) = \max_\theta \ell(\theta|X^D) \\ \end{align}

EMアルゴリズムでは\(F(\theta, q|X^D)\)を\(\theta\)と\(q\)に対して交互に最大化することで、\(\theta\)を近似させて最尤推定を行う。

\(q\)について最大化

\(F(\theta^{(t)}, q|X^D)\)を\(q\)について最大化する。

\(q\)と\(p(\cdot|X^D, \theta^{(t)})\)が一致するときKL情報量は最小値の0となるので、

\begin{align} q^{(t)} &= \underset{q}{\rm argmax} F(\theta^{(t)}, q|X^D) \\ &= \underset{q}{\rm argmax} \left( \ell(\theta^{(t)}|X^D) - D_{KL}\left(q\|p(\cdot|X^D, \theta^{(t)})\right) \right) \\ &= \underset{q}{\rm argmin} D_{KL}\left(q\|p(\cdot|X^D, \theta^{(t)})\right) \\ &= p(\cdot|X^D, \theta^{(t)}) \end{align}

\(\theta\)について最大化

\(F(\theta, q^{(t)}|X^D)\)を\(\theta\)について最大化する。

\begin{align} \theta^{(t+1)} &= \underset{\theta}{\rm argmax} F(\theta, q^{(t)}|X^D) \\ &= \underset{\theta}{\rm argmax} \int q^{(t)}(Z) \log\left( \frac{p(X^D, Z|\theta)}{q^{(t)}(Z)}\right) dZ \\ &= \underset{\theta}{\rm argmax} \int p(Z|X^D, \theta^{(t)}) \log\left( \frac{p(X^D, Z|\theta)}{p(Z|X^D, \theta^{(t)})}\right) dZ \\ &= \underset{\theta}{\rm argmax} \int p(Z|X^D, \theta^{(t)}) \log p(X^D, Z|\theta) dZ \\ &=: \underset{\theta}{\rm argmax} Q(\theta, \theta^{(t)}) \\ \end{align}

Qは以下で定義される関数。

\begin{align} Q(\theta, \theta^{(t)}) = \int p(Z|X^D, \theta^{(t)}) \log p(X^D, Z|\theta) dZ \end{align}

MAP推定の場合は最大化の対象が\(\ell(\theta|X^D)\)ではなく\(\ell(\theta|X^D)+\log p(\theta)\)なので、Q関数を以下のように置き換える。

\begin{align} Q_{MAP}(\theta, \theta^{(t)}) = \int p(Z|X^D, \theta^{(t)}) \log p(X^D, Z|\theta) dZ + \log p(\theta) \end{align}

更新アルゴリズム

EMアルゴリズム

\(\theta^{(0)}\)を適当に初期化し、以下の更新式で\(\theta^{(t)}\)を更新していく。

  • 最尤推定の場合

\begin{align} \theta^{(t+1)} &= \underset{\theta}{\rm argmax} Q(\theta, \theta^{(t)}) \\ &= \underset{\theta}{\rm argmax} \left( \int p(Z|X^D, \theta^{(t)}) \log p(X^D, Z|\theta) dZ \right) \\ \end{align}

  • MAP推定の場合

\begin{align} \theta^{(t+1)} &= \underset{\theta}{\rm argmax} Q_{MAP}(\theta, \theta^{(t)}) \\ &= \underset{\theta}{\rm argmax} \left( \int p(Z|X^D, \theta^{(t)}) \log p(X^D, Z|\theta) dZ + \log p(\theta) \right) \\ \end{align}

特徴

  • 潜在変数を含む確率モデルの最尤推定・MAP推定に適用することができる。
  • \(\theta^{(t)}\)は収束するが、\(\ell(\theta^{(t)}|X^D))\)の極限が最大値になるとは限らない。
  • \(\underset{\theta}{\rm argmax} Q(\theta, \theta^{(t)})\)を容易に計算できるモデルでなければ利用できない。

参考