EMアルゴリズム

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

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

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

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

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

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

: \ell(\theta|X):=\log p(X|\theta)=\log \sum_Zp(X,Z|\theta) 

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

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

:\ell(\hat{\theta}_{\mathrm{MLE}}|X)\ge \ell(\hat{\theta}^{(t)}|X) 

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

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

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

ここで*Q*は対数尤度関数 \log p (X,Z|\theta)のに関する条件付き期待値
: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) \,  　
である。実応用上は、\hat{\theta}^{(t)}の値が十分小さくなったと判定する何らかの条件を事前に定めておき、その条件を満たしたら上述のループを終了する。ループを終了する条件は、パラメータ値や対数尤度関数を使って定められる。

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

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

## 動作原理
EMアルゴリズムで我々が求めたいのは、X=(x_1,\ldots,x_n)を観測した際における対数尤度

:\ell(\theta|X):=\log p(X|\theta) 

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

: \mathcal{L}(q,\theta)
:=\sum_Z q(Z) \log {p(X,Z|\theta)\over q(Z)} 　　...()

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

:\mathrm{KL}(q||p_{X,\theta})=-\sum_Z q(Z)\log{p(Z|X,\theta) \over q(Z)} 　

を使って

:\mathcal{L}(q,\theta)=\ell(\theta|X)-\mathrm{KL}(q||p_{X,\theta}) 　...()

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

:\ell(\theta|X)\ge\mathcal{L}(q,\theta) 　

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

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

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

### 証明

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

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

:q(Z)=p(Z|X,\theta) 　

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

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

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

#### Mステップの証明
\mathcal{L}(q,\theta) の定義式()に\hat{q}^{(t)}=p_{X,\hat{\theta}^{(t)}}を代入すると、

:\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) 　

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

:\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アルゴリズム  であり、log-EMアルゴリズムをサブクラスとして含んでいる。α-EMアルゴリズムは適切なαを選ぶことにより、log-EMアルゴリズムよりも高速になる。また、log-EMが隠れマルコフモデル推定アルゴリズム（Baum-Welchアルゴリズム）を含んでいるように、α-EMアルゴリズムから高速なα-HMMアルゴリズムを得ることができる。 

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

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

## 脚注

## 参考文献

### 引用文献
* 
* 
* 

### その他の参考文献
* Robert Hogg, Joseph McKean and Allen Craig. *Introduction to Mathematical Statistics*. pp. 359-364.  Upper Saddle River, NJ: Pearson Prentice Hall, 2005.
* The on-line textbook: ''Information Theory, Inference, and Learning Algorithms'', by David J.C. MacKay includes simple examples of the E-M algorithm such as clustering using the soft K-means algorithm, and emphasizes the variational view of the E-M algorithm.
* A Gentle Tutorial of the EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models, by Jeff Bilmes includes a simplified derivation of the EM equations for Gaussian Mixtures and Gaussian Mixture Hidden Markov Models.
* ''Variational Algorithms for Approximate Bayesian Inference'', by M. J. Beal includes comparisons of EM to Variational Bayesian EM and derivations of several models including Variational Bayesian HMMs.
* ''The Expectation Maximization Algorithm'', by Frank Dellaert, gives an easier explanation of EM algorithm in terms of lowerbound maximization.
* ''The Expectation Maximization Algorithm: A short tutorial'', A self contained derivation of the EM Algorithm by Sean Borman.
* ''The EM Algorithm'', by Xiaojin Zhu.
* Geoffrey J. McLachlan and Thriyambakam Krishnan: *The EM Algorithm and Extensions*, Wiley series in probability and statistics, John Wiley & Sons, Inc., ISBN 0-471-12358-7 (1997).
* Geoffrey J. McLachlan and Thriyambakam Krishnan: *The EM Algorithm and Extensions*, 2nd Edition, Wiley & Sons Inc., ISBN 978-0-471-20170-0 (February 2008).　上記の改訂第2版。
* Kenneth Lange: *MM Optimization Algorithms*, SIAM,  (2016). ※ "EM" アルゴリズムの一般化として　"MM algorithm" を提唱している．
和書
* 小西貞則 ・越智義道 ・大森裕浩:「計算統計学の方法 ―ブートストラップ，EMアルゴリズム，MCMC―」、朝倉書店(シリーズ：予測と発見の科学、5)、ISBN 978-4-254-12785-0　(2008年3月25日)。
* 関原謙介：「ベイズ信号処理」、共立出版、ISBN 978-4-320-08574-9 (2015年4月11日)。
* 関原謙介：「ベイズ推論の基礎と信号処理への応用」
* 黒田正博:「EMアルゴリズム」、共立出版(シリーズ：統計学One Point、18巻)、ISBN 978-4-320-11269-8、(2020年7月31日)。