2.3 項目パラメタ推定
2.3.3 EM アルゴリズムによる周辺尤度の最大化
そこでBock & Aitkin (1981) はEMアルゴリズム (Dempster, Laird & Rubin, 1977) を用いる計
算方法を提案することにより計算効率を大幅に向上させている。EM アルゴリズム
(Expected-Maximum algorithm) は潜在変数や欠測値を含む統計モデルとデータにおけるパラメタ推定方法
の一般的な枠組みである。EMアルゴリズムにとよる最適化では潜在変数をデータから欠測した 情報とみなし，欠測データの情報を確率分布で暫定的に補った尤度関数を最大化する。Bock &
Aitkin (1981) が提案したEM アルゴリズムによる周辺最尤推定法は一般的なEM アルゴリズム
の導出とはやや異なる。
一般的な説明のためデータ行列を𝐔，潜在変数ベクトルを𝜽，推定したいパラメタベクトルを
𝛅とおく。まず一般的な EM アルゴリズムでは対数尤度関数ln𝐿(𝛅, 𝜽|𝐔)について暫定的に𝜽が観
測されたとして得られる完全データ対数尤度関数 (complete data likelihood function)，
ln𝐿𝐶(𝛅|𝐔, 𝜽), (2.31)
について考える。さらにこの関数に初期値𝛅𝟎を与えたときの条件付き期待値，
E[ln𝐿𝐶(𝛅|𝐔, 𝜽)|𝜹𝟎], (2.32)
を求めるステップが E ステップである。この期待値は潜在変数ベクトル𝜽の事後分布𝑔(𝜽|𝐮, 𝛅𝟎) により，
E[ln𝐿𝑪(𝛅|𝐔, 𝜽)|𝐔; 𝜹𝟎] = ∫ ln𝐿C(𝛅|𝐔, 𝜽)𝑔(𝜽|𝐔, 𝛅𝟎)d𝜽
∞
−∞
≡ 𝑄(𝜹|𝜹𝟎), (2.33)
と計算される。ただし𝜽の積分区間は−∞から∞である。この関数は𝛅𝟎と𝐔によって与えられる潜 在変数𝜽の事後分布によって，本来は観測され得ない潜在変数に関する情報が暫定的に補われた，
𝛅についての関数であるとみなすことができる。これを一般にQ関数，あるいは期待対数完全デ
ータ尤度関数 (expected log complete likelihood function) と呼ぶ。
期待値を計算し，Q関数を導出できれば，次はQ関数を最大化するような𝛅を何らかの手段で 推定し，更新された推定値を得る。このステップをMステップと呼ぶ。Q関数はもとの潜在変 数を含む対数尤度関数に比べていくらか簡単な表現になっていることが多いため，ときにこの 最大化の問題は解析的に解くことができる場合がある。解析的に解けない場合にはニュートン・
ラフソン法などの反復計算により最適解を計算する必要があるが，その場合に EM アルゴリズ ムが最適解にたどり着くことができるかどうかは，M ステップが適切に最適化されているかど うかに依存してしまう。
EM アルゴリズムによる反復が期待対数完全データ尤度関数を扱っているのにもかかわらず，
もとの対数尤度関数を単調増加させることができる理論的根拠を，樺島・上田 (2003)に則って 示す。まず，t回目の反復における対数尤度は，𝛅𝐭をその時点までに得られている暫定的な推定
値とすれば，潜在変数の事後分布を乗算，除算することで以下のように変形できる。
ln𝐿(𝜹|𝐔) = ln ∫ 𝑔(𝜽|𝐔, 𝛅𝐭)𝐿(𝜹, 𝜽|𝐔)
𝑔(𝜽|𝐔, 𝛅𝐭)d𝜽 = ln𝐸 [𝐿(𝜹, 𝜽|𝐔)
𝑔(𝜽|𝐔, 𝛅𝐭) ∣∣∣ 𝜹𝒕] . (2.34)
これは事後分布𝑔(𝜽|𝐔, 𝛅𝐭)についての期待値と見なせる。ここでJensen不等式
ln𝐸[𝑓(𝑥)] ≥ 𝐸[ln𝑓(𝑥)]
を利用すると，
ln𝐸 [𝐿(𝜹, 𝜽|𝐔)
𝑔(𝜽|𝐔, 𝛅𝐭) ∣∣∣ 𝜹𝒕] ≥ 𝐸 [ ln𝐿(𝜹, 𝜽|𝐔)
𝑔(𝜽|𝐔, 𝛅𝐭) ∣∣∣ 𝜹𝒕] = ∫ 𝑔(𝜽|𝐔, 𝛅𝐭) ln𝐿(𝜹, 𝜽|𝐔)
𝑔(𝜽|𝐔, 𝛅𝐭)d𝜽, (2.35)
というように不等式で表現できる。このとき事後分布𝑔(𝜽|𝐔, 𝛅𝐭)はベイズの定理により自然な形 として，
𝑔(𝜽|𝐔, 𝛅𝐭) = 𝐿(𝜃|𝐔, 𝜹)𝑞(𝜽)
∫ 𝐿(𝜃|𝐔, 𝜹)𝑞(𝜽)d𝜽, (2.36)
と得られる。以上により，Jensen不等式から対数尤度の下限値を代数的に表現することができた。
式 (2.35) の最右辺を変分下限と呼ぶこともある。この下限値の式を𝛅𝐭を所与としたときの𝛅の式 𝐹(𝜹|𝜹𝒕)とおいて，対数尤度ln𝐿(𝜹|𝐔)と𝐹(𝜹|𝜹𝒕)の差をとると，
ln𝐿(𝜹|𝐔) − 𝐹(𝜹|𝜹𝒕) = ln𝐿(𝜹|𝐔) − ∫ 𝑔(𝜽|𝐔, 𝛅𝐭) ln 𝐿(𝜽|𝐔, 𝜹)
𝑔(𝜽|𝐔, 𝛅𝐭)d𝜽, (2.37)
となり，右辺第一項に∫ 𝑔(𝜽|𝐔, 𝛅𝐭)d𝜽 = 1（確率分布なので）を乗じ，第二項の対数の分子を（尤 度だが数式としては確率𝑝(𝐔, 𝜽|𝛅) と同じなので）確率の乗法定理により分解すると式 (2.37) は，
= ln𝐿(𝜹|𝐔) ∫ 𝑔(𝜽|𝐔, 𝛅𝐭)d𝜽 − ∫ 𝑔(𝜽|𝐔, 𝛅𝐭) ln𝐿(𝜽|𝐔, 𝜹)𝐿(𝜹|𝐔)
𝑔(𝜽|𝐔, 𝛅𝐭) d𝜽, (2.38)
となるが，第一項のln𝐿(𝜹|𝐔)は𝜽を含まない関数であるため積分の中に戻してもよい。さらに第 二項の対数記号を分数中のそれぞれの関数に分配すると，さらに，
∫ 𝑔(𝜽|𝐔, 𝛅𝐭)ln𝐿(𝜹|𝐔)d𝜽 − ∫ 𝑔(𝜽|𝐔, 𝛅𝐭) [ln𝐿(𝜽|𝐔, 𝜹) + ln𝐿(𝜹|𝐔) − ln𝑔(𝜽|𝐔, 𝛅𝐭)]d𝜽, (2.39)
と変形される。ここで積分記号の中をまとめると式 (2.39) は，
∫ 𝑔(𝜽|𝐔, 𝛅𝐭) [ln𝐿(𝜹|𝐔) − {ln𝐿(𝜽|𝐔, 𝜹) + ln𝐿(𝜹|𝐔) − ln𝑔(𝜽|𝐔, 𝛅𝐭)}]d𝜽
= ∫ 𝑔(𝜽|𝐔, 𝛅𝐭) [ln𝐿(𝜽|𝐔, 𝜹) − ln𝑔(𝜽|𝐔, 𝛅𝐭)]d𝜽, (2.40)
と変形することができる。ここで確率分布の差異についての指標であるKLダイバージェンス，
KL(𝑃, 𝑄) = ∫ 𝑃(𝜽) ln𝑃(𝜽)
𝑄(𝜽)d𝜽, (2.41)
を利用してKL項で式 (2.40) を置き換えると式 (2.37) は，
ln𝐿(𝜹|𝐔) − 𝐹(𝜹|𝜹𝒕) = KL(𝑔(𝜽|𝐔, 𝛅𝐭), 𝐿(𝜽|𝐔, 𝜹)), (2.42)
となる。さらに変形すると，
ln𝐿(𝜹|𝐔) = 𝐹(𝜹|𝜹𝒕) + KL(𝑔(𝜽|𝐔, 𝛅𝐭), 𝐿(𝜽|𝐔, 𝜹)), (2.43)
が得られる。このとき式 (2.43) で𝜹 = 𝜹𝒕とおいたものとの差をとると，𝜹 = 𝜹𝒕のときのKL項は log1 = 0より消去され，
ln𝐿(𝜹|𝐔) − ln𝐿(𝜹𝒕|𝐔) = 𝐹(𝜹|𝜹𝒕) − 𝐹(𝜹𝒕|𝜹𝒕) + KL(𝑔(𝜽|𝐔, 𝛅𝐭), 𝐿(𝜽|𝐔, 𝜹)), (2.44)
となるが，KL項は非負であるので，
𝐹(𝜹|𝜹𝒕) ≥ 𝐹(𝜹𝒕|𝜹𝒕)
という条件さえ成り立てば，
ln𝐿(𝜹|𝐔) ≥ ln𝐿(𝜹𝒕|𝐔)
が保証される。ここで式 (2.33) および式 (2.34) より変分下限𝐹(𝜹|𝜹𝒕)は，
𝐹(𝜹|𝜹𝒕) = ∫ 𝑔(𝜽|𝐔, 𝛅𝐭) ln𝐿(𝜹, 𝜽|𝐔)d𝜽 − ∫ 𝑔(𝜽|𝐔, 𝛅𝐭) ln𝑔(𝜽|𝐔, 𝛅𝐭)d𝜽, (2.45)
と変形される。第一項は式 (2.33) の Q 関数そのものであり，第二項はパラメタ𝛅とは無関係の 項であるため，この式によりQ関数の最大化が𝐹(𝜹|𝜹𝒕)の最大化と等価であることが示せたため，
Q関数の最大化をすることで対数尤度ln𝐿(𝜹|𝐔)も最大化できることが証明できた。
EMアルゴリズムはEステップとMステップを，適当な収束基準に達するまで繰り返す推定 アルゴリズムである。一般にニュートン・ラフソン法などの計算方法に比べて収束に至るまでの 反復回数は非常に多いため，いくつかの高速化の手法が提案されている (たとえば Ramsay (1975) など)。