EMアルゴリズム (英 : expectation–maximization algorithm )とは、統計学 において、確率 モデル のパラメータを最尤推定 する手法の一つであり、観測不可能な潜在変数 に確率モデルが依存する場合に用いられる。EM法 、期待値最大化法 (きたいちさいだいかほう)とも呼ばれる。その一般性の高さから、機械学習 、音声認識 、因子分析 など、広汎な応用がある。
EMアルゴリズムは反復法 の一種であり、期待値(英 : expectation, E ) ステップと最大化 (英 : maximization, M )ステップを交互に繰り返すことで計算が進行する。Eステップでは、現在推定されている潜在変数の分布に基づいて、モデルの尤度の期待値を計算する。Mステップでは、E ステップで求まった尤度の期待値を最大化するようなパラメータを求める。M ステップで求まったパラメータは、次の E ステップで使われる潜在変数の分布を決定するために用いられる。
概要
セッティング・目標
今、2値x 、z を取る確率分布があり、その確率分布の確率密度関数
p
(
x
,
z
|
θ θ -->
)
{\displaystyle p(x,z|\theta )}
が未知の母数
θ θ -->
∈ ∈ -->
R
m
{\displaystyle \theta \in \mathbb {R} ^{m}}
によりパラメトライズされているとする。ここで
R
{\displaystyle \mathbb {R} }
は実数全体の集合を表す。
そして
p
(
x
,
z
|
θ θ -->
)
{\displaystyle p(x,z|\theta )}
に従って標本
(
x
1
,
z
1
)
,
… … -->
,
(
x
n
,
z
n
)
{\displaystyle (x_{1},z_{1}),\ldots ,(x_{n},z_{n})}
を独立に抽出したものの、何らかの事情で
Z
=
(
z
1
,
… … -->
,
z
n
)
{\displaystyle Z=(z_{1},\ldots ,z_{n})}
の値は観測できず、
X
=
(
x
1
,
… … -->
,
x
n
)
{\displaystyle X=(x_{1},\ldots ,x_{n})}
だけが観測できたとする。実応用上は例えば、
θ θ -->
=
(
θ θ -->
1
,
θ θ -->
2
)
{\displaystyle \theta =(\theta _{1},\theta _{2})}
という形をしており、まず観測不能な
z
i
∼ ∼ -->
p
1
(
z
|
θ θ -->
1
)
{\displaystyle z_{i}\sim p_{1}(z|\theta _{1})}
が選ばれた後、
z
i
{\displaystyle z_{i}}
に依存して観測可能な
x
i
∼ ∼ -->
p
2
(
x
|
θ θ -->
2
,
z
i
)
{\displaystyle x_{i}\sim p_{2}(x|\theta _{2},z_{i})}
が選ばれる、といったケースにEMアルゴリズムが使われる事が多いが、必ずしもこのケースにあてはまらなくてもよい。
簡単の為、記号を混用してX 、Z の同時確率分布の確率密度関数も
p
(
X
,
Z
|
θ θ -->
)
{\displaystyle p(X,Z|\theta )}
と書く。以下ではZ が離散変数の場合について説明するが、Z が連続変数の場合も総和を積分に置き換える以外は同様である[ 3] 。
このような状況において母数θ を最尤推定する事が我々の目標である。しかしZ を知らない場合の
X
=
(
x
1
,
… … -->
,
x
n
)
{\displaystyle X=(x_{1},\ldots ,x_{n})}
に関する対数尤度
ℓ ℓ -->
(
θ θ -->
|
X
)
:=
log
-->
p
(
X
|
θ θ -->
)
=
log
-->
∑ ∑ -->
Z
p
(
X
,
Z
|
θ θ -->
)
{\displaystyle \ell (\theta |X):=\log p(X|\theta )=\log \sum _{Z}p(X,Z|\theta )}
を最大値を直接計算するのは一般には簡単ではない。
EMアルゴリズムは、反復法により、数列
θ θ -->
^ ^ -->
(
t
)
{\displaystyle {\hat {\theta }}^{(t)}}
で対数尤度
ℓ ℓ -->
(
θ θ -->
^ ^ -->
(
t
)
|
X
)
{\displaystyle \ell ({\hat {\theta }}^{(t)}|X)}
が単調非減少 であるものを作るアルゴリズムである。最尤推定量を
θ θ -->
^ ^ -->
M
L
E
{\displaystyle {\hat {\theta }}_{\mathrm {MLE} }}
とすると、
ℓ ℓ -->
(
θ θ -->
^ ^ -->
M
L
E
|
X
)
≥ ≥ -->
ℓ ℓ -->
(
θ θ -->
^ ^ -->
(
t
)
|
X
)
{\displaystyle \ell ({\hat {\theta }}_{\mathrm {MLE} }|X)\geq \ell ({\hat {\theta }}^{(t)}|X)}
である事から、
ℓ ℓ -->
(
θ θ -->
^ ^ -->
M
L
E
|
X
)
{\displaystyle \ell ({\hat {\theta }}_{\mathrm {MLE} }|X)}
が有限であれば
ℓ ℓ -->
(
θ θ -->
^ ^ -->
(
t
)
|
X
)
{\displaystyle \ell ({\hat {\theta }}^{(t)}|X)}
の単調性より
ℓ ℓ -->
(
θ θ -->
^ ^ -->
(
t
)
|
X
)
{\displaystyle \ell ({\hat {\theta }}^{(t)}|X)}
は必ず収束する 。
アルゴリズム
EMアルゴリズムでは、以下の手順により数列
θ θ -->
^ ^ -->
(
0
)
,
θ θ -->
^ ^ -->
(
1
)
,
… … -->
{\displaystyle {\hat {\theta }}^{(0)},{\hat {\theta }}^{(1)},\ldots }
を作る[ 3] 。
初期値
θ θ -->
^ ^ -->
(
0
)
{\displaystyle {\hat {\theta }}^{(0)}}
を(何らかの方法で)選ぶ。
t
=
0
,
1
,
… … -->
{\displaystyle t=0,1,\ldots }
に対して以下を実行する
E ステップ :
p
(
Z
|
X
,
θ θ -->
^ ^ -->
(
t
)
)
{\displaystyle p(Z|X,{\hat {\theta }}^{(t)})}
を求める。
M ステップ :
θ θ -->
^ ^ -->
(
t
+
1
)
=
a
r
g
m
a
x
θ θ -->
Q
(
θ θ -->
|
θ θ -->
^ ^ -->
(
t
)
)
{\displaystyle {\hat {\theta }}^{(t+1)}={\underset {\theta }{\operatorname {arg\,max} }}\ Q(\theta |{\hat {\theta }}^{(t)})\,}
を求める。
ここでQ は対数尤度関数
log
-->
p
(
X
,
Z
|
θ θ -->
)
{\displaystyle \log p(X,Z|\theta )}
のZ に関する条件付き期待値
Q
(
θ θ -->
|
θ θ -->
(
t
)
)
:=
E
Z
|
X
,
θ θ -->
^ ^ -->
(
t
)
-->
[
log
-->
p
(
X
,
Z
|
θ θ -->
)
]
=
∑ ∑ -->
Z
p
(
Z
|
X
,
θ θ -->
^ ^ -->
(
t
)
)
log
-->
p
(
X
,
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 )\,}
である。実応用上は、
θ θ -->
^ ^ -->
(
t
)
{\displaystyle {\hat {\theta }}^{(t)}}
の値が十分小さくなったと判定する何らかの条件を事前に定めておき、その条件を満たしたら上述のループを終了する。ループを終了する条件は、パラメータ値や対数尤度関数を使って定められる[ 3] 。
留意点
EステップとMステップの切れ目は書籍により異なるので注意が必要である 。本項では次節の議論と整合性をとる為に文献[ 3] の切れ目に従ったが、文献[ 4] では
Q
(
θ θ -->
|
θ θ -->
^ ^ -->
(
t
)
)
{\displaystyle Q(\theta |{\hat {\theta }}^{(t)})}
を計算する所までがEステップであり、
Q
(
θ θ -->
|
θ θ -->
^ ^ -->
(
t
)
)
{\displaystyle Q(\theta |{\hat {\theta }}^{(t)})}
の
a
r
g
m
a
x
{\displaystyle \operatorname {arg\,max} }
を取るところだけがMステップである。
ステップの名称「E」と「M」はそれぞれExpectation(期待値)、Maximization(最大化)の略であり[ 4] 、文献[ 4] のようにEステップで
Q
(
θ θ -->
|
θ θ -->
^ ^ -->
(
t
)
)
{\displaystyle Q(\theta |{\hat {\theta }}^{(t)})}
を求める為に期待値を計算し、Mステップで
Q
(
θ θ -->
|
θ θ -->
^ ^ -->
(
t
)
)
{\displaystyle Q(\theta |{\hat {\theta }}^{(t)})}
の
a
r
g
m
a
x
{\displaystyle \operatorname {arg\,max} }
を取るところに名称の由来がある。
動作原理
EMアルゴリズムで我々が求めたいのは、
X
=
(
x
1
,
… … -->
,
x
n
)
{\displaystyle X=(x_{1},\ldots ,x_{n})}
を観測した際における対数尤度
ℓ ℓ -->
(
θ θ -->
|
X
)
:=
log
-->
p
(
X
|
θ θ -->
)
{\displaystyle \ell (\theta |X):=\log p(X|\theta )}
を最大化する母数
θ θ -->
{\displaystyle \theta }
であった。EMアルゴリズムの動作原理を説明する為、以下のような汎関数を考える:
L
(
q
,
θ θ -->
)
:=
∑ ∑ -->
Z
q
(
Z
)
log
-->
p
(
X
,
Z
|
θ θ -->
)
q
(
Z
)
{\displaystyle {\mathcal {L}}(q,\theta ):=\sum _{Z}q(Z)\log {p(X,Z|\theta ) \over q(Z)}}
...(Eq.1 )
ここで
q
(
Z
)
{\displaystyle q(Z)}
は任意の確率密度関数である。
p
X
,
θ θ -->
(
Z
)
:=
p
(
Z
|
X
,
θ θ -->
)
{\displaystyle p_{X,\theta }(Z):=p(Z|X,\theta )}
とすると、
p
(
Z
|
X
,
θ θ -->
)
p
(
X
|
θ θ -->
)
=
p
(
X
,
Z
|
θ θ -->
)
{\displaystyle p(Z|X,\theta )p(X|\theta )=p(X,Z|\theta )}
より、カルバック・ライブラー情報量
K
L
(
q
|
|
p
X
,
θ θ -->
)
=
− − -->
∑ ∑ -->
Z
q
(
Z
)
log
-->
p
(
Z
|
X
,
θ θ -->
)
q
(
Z
)
{\displaystyle \mathrm {KL} (q||p_{X,\theta })=-\sum _{Z}q(Z)\log {p(Z|X,\theta ) \over q(Z)}}
を使って
L
(
q
,
θ θ -->
)
=
ℓ ℓ -->
(
θ θ -->
|
X
)
− − -->
K
L
(
q
|
|
p
X
,
θ θ -->
)
{\displaystyle {\mathcal {L}}(q,\theta )=\ell (\theta |X)-\mathrm {KL} (q||p_{X,\theta })}
...(Eq.2 )
と書ける事が分かる。カルバック・ライブラー情報量が常に非負である事(ギブスの不等式 )から、
ℓ ℓ -->
(
θ θ -->
|
X
)
≥ ≥ -->
L
(
q
,
θ θ -->
)
{\displaystyle \ell (\theta |X)\geq {\mathcal {L}}(q,\theta )}
であるので、
L
(
q
,
θ θ -->
)
{\displaystyle {\mathcal {L}}(q,\theta )}
は
ℓ ℓ -->
(
θ θ -->
|
X
)
{\displaystyle \ell (\theta |X)}
の下限になっている。EMアルゴリズムはこの下限
L
(
q
,
θ θ -->
)
{\displaystyle {\mathcal {L}}(q,\theta )}
を逐次的に改善していくことで、
ℓ ℓ -->
(
θ θ -->
|
X
)
{\displaystyle \ell (\theta |X)}
を可能な限り最大化するアルゴリズムである。すなわち、EステップとMステップは以下のように書き換えられる事を示す事ができる[ 3] :
E ステップ :
q
^ ^ -->
(
t
)
=
a
r
g
m
a
x
q
L
(
q
,
θ θ -->
^ ^ -->
(
t
)
)
{\displaystyle {\hat {q}}^{(t)}={\underset {q}{\operatorname {arg\,max} }}{\mathcal {L}}(q,{\hat {\theta }}^{(t)})}
を求める。
M ステップ :
θ θ -->
^ ^ -->
(
t
+
1
)
=
a
r
g
m
a
x
θ θ -->
L
(
q
^ ^ -->
(
t
)
,
θ θ -->
)
{\displaystyle {\hat {\theta }}^{(t+1)}={\underset {\theta }{\operatorname {arg\,max} }}{\mathcal {L}}({\hat {q}}^{(t)},\theta )}
を求める。
この事実から対数尤度
ℓ ℓ -->
(
θ θ -->
^ ^ -->
(
t
)
|
X
)
{\displaystyle \ell ({\hat {\theta }}^{(t)}|X)}
の単調非減少性が明らかに従う。
(但し反復法の常として、初期値しだいでは尤度の最大点ではない極大点に到達してそこで停止する可能性がある。)
証明
本節ではEステップ、Mステップが上述のように書き換えられることを示す。本節の証明は文献[ 3] を参考にした。
Eステップの証明
カルバック・ライブラー情報量
K
L
(
q
|
|
p
X
,
θ θ -->
)
{\displaystyle \mathrm {KL} (q||p_{X,\theta })}
が最小値0になるのは
q
=
p
θ θ -->
,
X
{\displaystyle q=p_{\theta ,X}}
の場合だけであった事から、(Eq.2 )より
L
(
q
,
θ θ -->
)
{\displaystyle {\mathcal {L}}(q,\theta )}
は
q
(
Z
)
=
p
(
Z
|
X
,
θ θ -->
)
{\displaystyle q(Z)=p(Z|X,\theta )}
が満たされる場合に最大値を取る。すなわちEMアルゴリズムにおけるEステップは、
θ θ -->
=
θ θ -->
^ ^ -->
(
t
)
{\displaystyle \theta ={\hat {\theta }}^{(t)}}
を固定したままの状態で、
L
(
q
,
θ θ -->
)
{\displaystyle {\mathcal {L}}(q,\theta )}
を最大化する
q
{\displaystyle q}
である
q
^ ^ -->
(
t
)
:=
p
X
,
θ θ -->
^ ^ -->
(
t
)
=
a
r
g
m
a
x
q
L
(
q
,
θ θ -->
^ ^ -->
(
t
)
)
{\displaystyle {\hat {q}}^{(t)}:=p_{X,{\hat {\theta }}^{(t)}}={\underset {q}{\operatorname {arg\,max} }}{\mathcal {L}}(q,{\hat {\theta }}^{(t)})}
を求めるステップである。
Mステップの証明
L
(
q
,
θ θ -->
)
{\displaystyle {\mathcal {L}}(q,\theta )}
の定義式(Eq.1 )に
q
^ ^ -->
(
t
)
=
p
X
,
θ θ -->
^ ^ -->
(
t
)
{\displaystyle {\hat {q}}^{(t)}=p_{X,{\hat {\theta }}^{(t)}}}
を代入すると、
L
(
q
^ ^ -->
(
t
)
,
θ θ -->
)
=
∑ ∑ -->
Z
p
(
Z
|
X
,
θ θ -->
(
t
)
)
log
-->
p
(
X
,
Z
|
θ θ -->
)
p
(
Z
|
X
,
θ θ -->
(
t
)
)
=
Q
(
θ θ -->
|
θ θ -->
(
t
)
)
− − -->
H
X
,
θ θ -->
(
t
)
(
Z
)
{\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)}
が成立し(ここで
H
X
,
θ θ -->
(
t
)
(
Z
)
=
∑ ∑ -->
Z
p
(
Z
|
X
,
θ θ -->
(
t
)
)
log
-->
p
(
Z
|
X
,
θ θ -->
(
t
)
)
{\displaystyle H_{X,\theta ^{(t)}}(Z)=\textstyle \sum _{Z}p(Z|X,\theta ^{(t)})\log p(Z|X,\theta ^{(t)})}
は条件付きエントロピー )、上式右辺第二項はθ に依存しないので、
θ θ -->
^ ^ -->
(
t
+
1
)
=
a
r
g
m
a
x
θ θ -->
Q
(
θ θ -->
|
θ θ -->
^ ^ -->
(
t
)
)
=
a
r
g
m
a
x
θ θ -->
L
(
p
X
,
θ θ -->
^ ^ -->
(
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アルゴリズムでは、期待値の評価において潜在変数のとりうる値すべてを列挙することが必要なため、効率的に扱える分布が限られていた。しかしその後、マルコフ連鎖モンテカルロ法 や変分ベイズ法 (英語版 ) が考案されたことにより、より一般の分布でも現実的な時間での計算が可能になった。
脚注
^ a b c d e f #PRML pp.156, 164-171
^ a b c #ESL pp.316-317.
^ Matsuyama, Yasuo (2003). “The α-EM algorithm: Surrogate likelihood maximization using α-logarithmic information measures”. IEEE Transactions on Information Theory 49 (3): 692-706.
^ Matsuyama, Yasuo (2011). “Hidden Markov model estimation based on alpha-EM algorithm: Discrete and continuous alpha-HMMs”. International Joint Conference on Neural Networks : 808-816.
^ Dempster, A.P., Laird, N.M., Rubin, D.B., (1977). “Maximum Likelihood from Incomplete Data via the EM Algorithm”. Journal of the Royal Statistical Society . Series B (Methodological) 39 (1): 1–38. JSTOR2984875 . MR 0501537 .
参考文献
引用文献
その他の参考文献
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版。
小西貞則 ・越智義道 ・大森裕浩:「計算統計学の方法 ―ブートストラップ,EMアルゴリズム,MCMC―」、朝倉書店(シリーズ:予測と発見の科学、5)、ISBN 978-4-254-12785-0 、2008年3月25日。
関原謙介:「ベイズ信号処理」、共立出版、ISBN 978-4-320-08574-9 、2015年4月。
関原謙介:「ベイズ推論の基礎と信号処理への応用」
Kenneth Lange: "MM Optimization Algorithms", SIAM, ISBN 978-1-611974-39-3 (2016). ※ "MM algorithm " は "EM" アルゴリズムの一般化として提唱されている.
黒田正博:「EMアルゴリズム」、共立出版(シリーズ:統計学One Point、18巻)、ISBN 978-4-320-11269-8 、2020年07月31日。