EMアルゴリズム
この記事では最尤推定およびMAP推定を解くための手法「EMアルゴリズム」について解説する。
(一連の記事では、紛らわしくはあるが確率分布を表すpという記号を異なる確率分布にも使い回して表記を簡素化している。)
関連記事
- EMアルゴリズム (本記事)
- 変分Bayes
- VAE: 変分オートエンコーダー
目次
最尤推定・MAP推定
確率変数xがパラメーターθの確率分布p(x|θ)に従うことがわかっているとする。
pに従って独立に発生した幾つかのデータXD={xDn}Nn=1が既に観測されているとき、次に発生するxの分布をより正確に推測したい。
そのための基礎的な方針・問題設定が、最尤推定とMAP推定である。
どちらも正規分布や指数分布などの指数型の簡素な分布を扱う場合は解析的に解くこともできるが、より複雑な分布の場合は様々なアルゴリズムを用いなければ解くことが困難である。
最尤推定
最尤推定(Maximum likelihood estimation, MLE)では、尤度を最大化することを目指す。
- xおよびxDnが従う確率分布p(x|θ)が関数として既知。
- 分布pに従って独立に生成されたデータの集合XDが観測されている。
観測データに最も適合するθをただ一つ求める。
パラメーターがθであるときにデータXDが生成される確率L(θ|XD)を尤度と呼ぶ。
L(θ|XD)=N∏n=1p(xDn|θ)
また、尤度の対数ℓ(θ|XD)を対数尤度と呼ぶ。
ℓ(θ|XD)=log(L(θ|XD))=∑nlogp(xDn|θ)
最尤推定では尤度を最大化するパラメーターを求めることを目的とする。
θML=argmaxθL(θ|XD)=argmaxθℓ(θ|XD)
MAP推定
MAP推定(Maximum a posteriori estimation)では、θの事前確率が既知であるという前提を追加して、θの事後確率を最大化することを目指す。
- xおよびxDnが従う確率分布p(x|θ)が関数として既知。
- θの(事前)確率分布p(θ)が関数として既知。
- 分布pに従って独立に生成されたデータの集合XDが観測されている。
観測データに最も適合するθをただ一つ求める。
Bayesの定理より、θの事後確率p(θ|XD)は以下のようになる。
p(θ|XD)=p(θ,XD)p(XD)=p(XD|θ)p(θ)p(XD)
MAP推定ではθの事後確率を最大化するパラメーターを求めることを目的とする。
θMAP=argmaxθp(θ|XD)=argmaxθp(XD|θ)p(θ)p(XD)=argmaxθp(XD|θ)p(θ)=argmaxθ(ℓ(θ|XD)+logp(θ))
EMアルゴリズム
EMアルゴリズム(expectation–maximization algorithm)は、1977年にDempster, Laird, Rubinらによって考案された最尤推定・MAP推定の反復アルゴリズム。
- p(x,z|θ),p(z|x,θ)が関数として既知。
- (MAP推定の場合) p(θ)も関数として既知。
- 分布pに従って生成された値(x,z)の内、xだけが観測されていて独立なデータの集合XDとして値が得られている。
観測データに最も適合するθをただ一つ求める。
zのように直接観測することはできないが、前提として定めた確率モデルから推測することができる確率変数のことを潜在変数と言う。
アルゴリズムの導出
Jensenの不等式から、対数尤度に対して以下の不等式が成立する。
ℓ(θ|XD)=logp(XD|θ)=log(∫p(XD,Z|θ)dZ)=log(∫q(Z)p(XD,Z|θ)q(Z)dZ)≥∫q(Z)log(p(XD,Z|θ)q(Z))dZ=:F(θ,q|XD)
ここで、q(Z)はZの任意の確率分布。F(θ,q|XD)を変分下限(Evidence lower bound, ELBO)と呼ぶ。
変分下限を再び調べると、以下のようになる。
F(θ,q|XD)=∫q(Z)log(p(XD,Z|θ)q(Z))dZ=∫q(Z)log(p(Z|XD,θ)p(XD|θ)q(Z))dZ=∫q(Z)logp(XD|θ)dZ+∫q(Z)log(p(Z|XD,θ)q(Z))dZ=p(XD|θ)−∫q(Z)log(q(Z)p(Z|XD,θ))dZ=ℓ(θ|XD)−DKL(q‖
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}
更新アルゴリズム
\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)})を容易に計算できるモデルでなければ利用できない。