EMアルゴリズム

EMアルゴリズム（英: expectation–maximization algorithm）とは、統計学において、確率モデルのパラメータを最尤推定する手法の一つであり、観測不可能な潜在変数に確率モデルが依存する場合に用いられる。EM法、期待値最大化法（きたいちさいだいかほう）^[1][2]とも呼ばれる。その一般性の高さから、機械学習、音声認識、因子分析など、広汎な応用がある^[1]。

EMアルゴリズムは反復法の一種であり、期待値（英: expectation, E）ステップと最大化（英: maximization, M）ステップを交互に繰り返すことで計算が進行する。Eステップでは、現在推定されている潜在変数の分布に基づいて、モデルの尤度の期待値を計算する。Mステップでは、E ステップで求まった尤度の期待値を最大化するようなパラメータを求める。M ステップで求まったパラメータは、次の E ステップで使われる潜在変数の分布を決定するために用いられる。

概要

セッティング・目標

今、2値x、zを取る確率分布があり、その確率分布の確率密度関数${\displaystyle p(x,z|\theta )}$が未知の母数${\displaystyle \theta \in \mathbb {R} ^{m}}$によりパラメトライズされているとする。ここで${\displaystyle \mathbb {R} }$は実数全体の集合を表す。

そして${\displaystyle p(x,z|\theta )}$に従って標本${\displaystyle (x_{1},z_{1}),\ldots ,(x_{n},z_{n})}$を独立に抽出したものの、何らかの事情で${\displaystyle Z=(z_{1},\ldots ,z_{n})}$の値は観測できず、${\displaystyle X=(x_{1},\ldots ,x_{n})}$だけが観測できたとする。実応用上は例えば、${\displaystyle \theta =(\theta _{1},\theta _{2})}$という形をしており、まず観測不能な${\displaystyle z_{i}\sim p_{1}(z|\theta _{1})}$が選ばれた後、${\displaystyle z_{i}}$に依存して観測可能な${\displaystyle x_{i}\sim p_{2}(x|\theta _{2},z_{i})}$が選ばれる、といったケースにEMアルゴリズムが使われる事が多いが、必ずしもこのケースにあてはまらなくてもよい。

簡単の為、記号を混用してX、Zの同時確率分布の確率密度関数も${\displaystyle p(X,Z|\theta )}$と書く。以下ではZが離散変数の場合について説明するが、Zが連続変数の場合も総和を積分に置き換える以外は同様である^[3]。

このような状況において母数θを最尤推定する事が我々の目標である。しかしZを知らない場合の${\displaystyle X=(x_{1},\ldots ,x_{n})}$に関する対数尤度

${\displaystyle \ell (\theta |X):=\log p(X|\theta )=\log \sum _{Z}p(X,Z|\theta )}$

を最大値を直接計算するのは一般には簡単ではない。

EMアルゴリズムは、反復法により、数列${\displaystyle {\hat {\theta }}^{(t)}}$で対数尤度${\displaystyle \ell ({\hat {\theta }}^{(t)}|X)}$が単調非減少であるものを作るアルゴリズムである。最尤推定量を${\displaystyle {\hat {\theta }}_{\mathrm {MLE} }}$とすると、

${\displaystyle \ell ({\hat {\theta }}_{\mathrm {MLE} }|X)\geq \ell ({\hat {\theta }}^{(t)}|X)}$

である事から、${\displaystyle \ell ({\hat {\theta }}_{\mathrm {MLE} }|X)}$が有限であれば${\displaystyle \ell ({\hat {\theta }}^{(t)}|X)}$の単調性より${\displaystyle \ell ({\hat {\theta }}^{(t)}|X)}$は必ず収束する。

アルゴリズム

EMアルゴリズムでは、以下の手順により数列${\displaystyle {\hat {\theta }}^{(0)},{\hat {\theta }}^{(1)},\ldots }$を作る^[3]。

• 初期値${\displaystyle {\hat {\theta }}^{(0)}}$を（何らかの方法で）選ぶ。
• ${\displaystyle t=0,1,\ldots }$に対して以下を実行する
  • E ステップ: ${\displaystyle p(Z|X,{\hat {\theta }}^{(t)})}$を求める。
  • M ステップ: ${\displaystyle {\hat {\theta }}^{(t+1)}={\underset {\theta }{\operatorname {arg\,max} }}\ Q(\theta |{\hat {\theta }}^{(t)})\,}$を求める。

ここでQは対数尤度関数${\displaystyle \log p(X,Z|\theta )}$のZに関する条件付き期待値

${\displaystyle Q(\theta |\theta ^{(t)}):=\operatorname {E} _{Z|X,{\hat {\theta }}^{(t)}}{\big [}\log p(X,Z|\theta ){\big ]}=\sum _{Z}p(Z|X,{\hat {\theta }}^{(t)})\log p(X,Z|\theta )\,}$ 　

である。実応用上は、${\displaystyle {\hat {\theta }}^{(t)}}$の値が十分小さくなったと判定する何らかの条件を事前に定めておき、その条件を満たしたら上述のループを終了する。ループを終了する条件は、パラメータ値や対数尤度関数を使って定められる^[3]。

留意点

EステップとMステップの切れ目は書籍により異なるので注意が必要である。本項では次節の議論と整合性をとる為に文献^[3]の切れ目に従ったが、文献^[4]では${\displaystyle Q(\theta |{\hat {\theta }}^{(t)})}$を計算する所までがEステップであり、${\displaystyle Q(\theta |{\hat {\theta }}^{(t)})}$の${\displaystyle \operatorname {arg\,max} }$を取るところだけがMステップである。

ステップの名称「E」と「M」はそれぞれExpectation（期待値）、Maximization（最大化）の略であり^[4]、文献^[4]のようにEステップで${\displaystyle Q(\theta |{\hat {\theta }}^{(t)})}$を求める為に期待値を計算し、Mステップで${\displaystyle Q(\theta |{\hat {\theta }}^{(t)})}$の${\displaystyle \operatorname {arg\,max} }$を取るところに名称の由来がある。

動作原理

EMアルゴリズムで我々が求めたいのは、${\displaystyle X=(x_{1},\ldots ,x_{n})}$を観測した際における対数尤度

${\displaystyle \ell (\theta |X):=\log p(X|\theta )}$

を最大化する母数${\displaystyle \theta }$であった。EMアルゴリズムの動作原理を説明する為、以下のような汎関数を考える：

${\displaystyle {\mathcal {L}}(q,\theta ):=\sum _{Z}q(Z)\log {p(X,Z|\theta ) \over q(Z)}}$　　...(Eq.1)

ここで${\displaystyle q(Z)}$は任意の確率密度関数である。${\displaystyle p_{X,\theta }(Z):=p(Z|X,\theta )}$とすると、${\displaystyle p(Z|X,\theta )p(X|\theta )=p(X,Z|\theta )}$より、カルバック・ライブラー情報量

${\displaystyle \mathrm {KL} (q||p_{X,\theta })=-\sum _{Z}q(Z)\log {p(Z|X,\theta ) \over q(Z)}}$　

を使って

${\displaystyle {\mathcal {L}}(q,\theta )=\ell (\theta |X)-\mathrm {KL} (q||p_{X,\theta })}$　...(Eq.2)

と書ける事が分かる。カルバック・ライブラー情報量が常に非負である事（ギブスの不等式）から、

${\displaystyle \ell (\theta |X)\geq {\mathcal {L}}(q,\theta )}$　

であるので、${\displaystyle {\mathcal {L}}(q,\theta )}$は${\displaystyle \ell (\theta |X)}$の下限になっている。EMアルゴリズムはこの下限${\displaystyle {\mathcal {L}}(q,\theta )}$を逐次的に改善していくことで、${\displaystyle \ell (\theta |X)}$を可能な限り最大化するアルゴリズムである。すなわち、EステップとMステップは以下のように書き換えられる事を示す事ができる^[3]：

• E ステップ: ${\displaystyle {\hat {q}}^{(t)}={\underset {q}{\operatorname {arg\,max} }}{\mathcal {L}}(q,{\hat {\theta }}^{(t)})}$を求める。
• M ステップ: ${\displaystyle {\hat {\theta }}^{(t+1)}={\underset {\theta }{\operatorname {arg\,max} }}{\mathcal {L}}({\hat {q}}^{(t)},\theta )}$を求める。

この事実から対数尤度${\displaystyle \ell ({\hat {\theta }}^{(t)}|X)}$の単調非減少性が明らかに従う。 （但し反復法の常として、初期値しだいでは尤度の最大点ではない極大点に到達してそこで停止する可能性がある。）

証明

本節ではEステップ、Mステップが上述のように書き換えられることを示す。本節の証明は文献^[3]を参考にした。

Eステップの証明

カルバック・ライブラー情報量${\displaystyle \mathrm {KL} (q||p_{X,\theta })}$が最小値0になるのは${\displaystyle q=p_{\theta ,X}}$の場合だけであった事から、(Eq.2)より${\displaystyle {\mathcal {L}}(q,\theta )}$は

${\displaystyle q(Z)=p(Z|X,\theta )}$　

が満たされる場合に最大値を取る。すなわちEMアルゴリズムにおけるEステップは、${\displaystyle \theta ={\hat {\theta }}^{(t)}}$を固定したままの状態で、${\displaystyle {\mathcal {L}}(q,\theta )}$を最大化する${\displaystyle q}$である

${\displaystyle {\hat {q}}^{(t)}:=p_{X,{\hat {\theta }}^{(t)}}={\underset {q}{\operatorname {arg\,max} }}{\mathcal {L}}(q,{\hat {\theta }}^{(t)})}$　

を求めるステップである。

Mステップの証明

${\displaystyle {\mathcal {L}}(q,\theta )}$の定義式(Eq.1)に${\displaystyle {\hat {q}}^{(t)}=p_{X,{\hat {\theta }}^{(t)}}}$を代入すると、

${\displaystyle {\mathcal {L}}({\hat {q}}^{(t)},\theta )=\sum _{Z}p(Z|X,\theta ^{(t)})\log {p(X,Z|\theta ) \over p(Z|X,\theta ^{(t)})}=Q(\theta |\theta ^{(t)})-H_{X,\theta ^{(t)}}(Z)}$　

が成立し（ここで${\displaystyle H_{X,\theta ^{(t)}}(Z)=\textstyle \sum _{Z}p(Z|X,\theta ^{(t)})\log p(Z|X,\theta ^{(t)})}$は条件付きエントロピー）、上式右辺第二項はθに依存しないので、

${\displaystyle {\hat {\theta }}^{(t+1)}={\underset {\theta }{\operatorname {arg\,max} }}Q(\theta |{\hat {\theta }}^{(t)})={\underset {\theta }{\operatorname {arg\,max} }}{\mathcal {L}}(p_{X,{\hat {\theta }}^{(t)}},\theta )}$　

が成立する。

一般化

EMアルゴリズムは観測データの対数尤度を、E ステップとM ステップの繰り返しにより最大化するアルゴリズムであるので、正確にはlog-EMアルゴリズムというべきものである。log関数にはα-logとよばれる一般化された対数があるので、それを用いるとlog-EMを特例として含むアルゴリズムを作り上げることができる。ただし、この場合は尤度ではなくてα-log尤度比とαダイバージェンスを用いて基本等式を導くことになる。このようにして得られたものがα-EMアルゴリズム ^[5] であり、log-EMアルゴリズムをサブクラスとして含んでいる。α-EMアルゴリズムは適切なαを選ぶことにより、log-EMアルゴリズムよりも高速になる。また、log-EMが隠れマルコフモデル推定アルゴリズム（Baum-Welchアルゴリズム）を含んでいるように、α-EMアルゴリズムから高速なα-HMMアルゴリズムを得ることができる。 ^[6]

歴史

EMアルゴリズムは、アーサー・デンプスター（英語版）、ナン・レアード（英語版）、ドナルド・ルービンによる1977年の論文^[7]で導入され、その名が付けられた。彼らは、EMアルゴリズムがほかの複数の著者によって「特殊な文脈でなんども提案されてきた」("proposed many times in special circumstances") ことを述べた上で、EMアルゴリズムの一般化を行い、その背後にある理論を追求した。

本来のEMアルゴリズムでは、期待値の評価において潜在変数のとりうる値すべてを列挙することが必要なため、効率的に扱える分布が限られていた。しかしその後、マルコフ連鎖モンテカルロ法や変分ベイズ法（英語版）が考案されたことにより、より一般の分布でも現実的な時間での計算が可能になった^[1][8]。