Iterative method for finding maximum likelihood estimates in statistical models
統計学 において 、 期待最大化 ( EM ) アルゴリズムは、 統計モデル における パラメータ の(局所的) 最大尤度 または 最大事後尤度 (MAP)推定値を求める 反復法 であり 、モデルは観測されない 潜在変数 に依存する。 [1] EM反復は、パラメータの現在の推定値を使用して評価された対数尤度の期待 値の関数を作成する期待値(E)ステップと、 E ステップで見つかった期待対数尤度を最大化するパラメータを計算する最大化(M)ステップを交互に実行する。これらのパラメータ推定値は、次のEステップで潜在変数の分布を決定するために使用される。例えば、 ガウス 分布の混合を推定したり、多重線形回帰問題を解いたりするのに 使用できる。 [2]
オールドフェイスフル 噴火データのEMクラスタリング 。ランダムな初期モデル(軸のスケールが異なるため、非常に平坦で幅の広い2つの楕円のように見える)を観測データに当てはめる。最初の反復ではモデルは大きく変化するが、その後、 間欠泉 の2つのモードに収束する。ELKIを用いて可視 化 。
歴史
EMアルゴリズムは、 アーサー・デンプスター 、 ナン・レアード 、 ドナルド・ルービン による1977年の古典的な論文で説明され、その名前が付けられました。 [3]彼らは、この方法が以前の著者によって「特別な状況下で何度も提案された」ことを指摘しました。最も初期のものの一つは、 セドリック・スミス による対立遺伝子頻度を推定するための遺伝子計数法です 。 [4]もう1つは、1958年に HOハートレー によって、1977年にハートレーとホッキングによって提案され 、デンプスター・レアード・ルービンの論文のアイデアの多くはこの方法に由来しています。 [5] もう1つは、1977年にSK・ン、スリヤムバカム・クリシュナン、GJ・マクラクランによって提案されました。 [6] ハートレーのアイデアは、グループ化された任意の離散分布に拡張できます。指数族に対するEM法の非常に詳細な扱いは、ロルフ・サンドバーグが ペル・マーティン=レーフ および アンダース・マーティン=レーフ との共同研究を経て、学位論文およびいくつかの論文 [7] [8] [9] で発表しました。 [10] [11] [12] [13] [14] 1977年のデンプスター・レアード・ルビンの論文では、この手法が一般化され、より広範な問題に対する収束解析の概要が示されました。デンプスター・レアード・ルビンの論文は、EM法を統計解析の重要なツールとして確立しました。メンとヴァン・ダイク(1997)も参照してください。
デンプスター・レアード・ルビン法の収束解析には欠陥があり、正しい収束解析は 1983年に CFジェフ・ウーによって発表された。 [15]ウーの証明は、デンプスター・レアード・ルビンが主張したように、 指数族
の外側でもEM法が収束することを証明した 。 [15]
導入
EMアルゴリズムは、 方程式を直接解くことができない場合に、 統計モデル の(局所的) 最大尤度パラメータを求めるために使用されます。通常、これらのモデルには、未知の パラメータ と既知のデータ観測に加えて、 潜在変数が 含まれます。つまり、データの中に 欠損値が 存在するか、あるいは、さらに観測されていないデータポイントが存在すると仮定することで、モデルをより単純に定式化できます。例えば、 混合モデルは 、各観測データポイントに対応する観測されていないデータポイント、つまり潜在変数があり、各データポイントが属する混合成分を特定すると仮定することで、より単純に記述できます。
最大尤度解を求めるには、通常、尤度関数をすべての未知数、パラメータ、潜在変数について微分し、得られた方程式を同時に解く必要があります 。 潜在 変数 を含む統計モデルでは、これは通常不可能です。その代わりに、結果として得られる方程式は通常、パラメータの解を求めるには潜在変数の値が必要であり、その逆もまた同様である、互いに連動した方程式の集合となります。しかし、ある方程式の集合を別の方程式の集合に代入すると、解けない方程式が生成されます。
EMアルゴリズムは、これら2組の方程式を数値的に解く方法があるという観察から出発する。2組の未知数のうちの1組に任意の値を選び、それを用いて2組目の未知数を推定し、次にこれらの新しい値を用いて1組目の未知数をより正確に推定し、そして得られた値が両方とも固定点に収束するまで、2組を交互に繰り返すだけでよい。これがうまくいくかどうかは明らかではないが、この文脈では証明できる。さらに、その点における尤度の導関数が(任意に)ゼロに近いことが証明できる。これは、その点が局所的最大値または鞍点のいずれかであることを意味する 。 [ 15] 一般に、複数の最大値が発生する可能性があり、必ずしも全体的最大値が見つかるとは限らない。尤度の中には 特異点 、すなわち無意味な最大値を持つものもある。例えば、混合モデルにおいてEMによって見つかる可能性のある 解 の1つは、成分の1つをゼロ分散に設定し、同じ成分の平均パラメータをデータ点の1つと等しく設定することである。
説明
シンボル
観測データの集合 、観測されていない潜在データまたは 欠損値 の集合、未知のパラメータのベクトル、および 尤度関数を生成する 統計モデル が与えられた場合 、未知のパラメータの 最大 尤度推定 値(MLE)は、観測データの
周辺尤度 を最大化することによって決定される。
X
{\displaystyle \mathbf {X} }
Z
{\displaystyle \mathbf {Z} }
θ
{\displaystyle {\boldsymbol {\theta }}}
L
(
θ
;
X
,
Z
)
=
p
(
X
,
Z
∣
θ
)
{\displaystyle L({\boldsymbol {\theta }};\mathbf {X} ,\mathbf {Z} )=p(\mathbf {X} ,\mathbf {Z} \mid {\boldsymbol {\theta }})}
L
(
θ
;
X
)
=
p
(
X
∣
θ
)
=
∫
p
(
X
,
Z
∣
θ
)
d
Z
=
∫
p
(
X
∣
Z
,
θ
)
p
(
Z
∣
θ
)
d
Z
{\displaystyle L({\boldsymbol {\theta }};\mathbf {X} )=p(\mathbf {X} \mid {\boldsymbol {\theta }})=\int p(\mathbf {X} ,\mathbf {Z} \mid {\boldsymbol {\theta }})\,d\mathbf {Z} =\int p(\mathbf {X} \mid \mathbf {Z} ,{\boldsymbol {\theta }})p(\mathbf {Z} \mid {\boldsymbol {\theta }})\,d\mathbf {Z} }
しかし、 は観測されておらず、に達するまで の の分布は不明である ため、この量は扱いにくいことがよくあります 。
Z
{\displaystyle \mathbf {Z} }
Z
{\displaystyle \mathbf {Z} }
θ
{\displaystyle {\boldsymbol {\theta }}}
EMアルゴリズム
EM アルゴリズムは、次の 2 つのステップを繰り返し適用して、周辺尤度の最大尤度推定値を見つけようとします。
期待値ステップ(Eステップ) :を、 与えられた の現在 の 条件付き分布 とパラメータの現在の推定値 に関して、 の対数 尤度関数 の 期待値 として定義します 。
Q
(
θ
∣
θ
(
t
)
)
{\displaystyle Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})}
θ
{\displaystyle {\boldsymbol {\theta }}}
Z
{\displaystyle \mathbf {Z} }
X
{\displaystyle \mathbf {X} }
θ
(
t
)
{\displaystyle {\boldsymbol {\theta }}^{(t)}}
Q
(
θ
∣
θ
(
t
)
)
=
E
Z
∼
p
(
⋅
|
X
,
θ
(
t
)
)
[
log
p
(
X
,
Z
|
θ
)
]
:=
∫
log
p
(
X
,
Z
|
θ
)
p
(
Z
|
X
,
θ
(
t
)
)
d
Z
{\displaystyle Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})=\operatorname {E} _{\mathbf {Z} \sim p(\cdot |\mathbf {X} ,{\boldsymbol {\theta }}^{(t)})}\left[\log p(\mathbf {X} ,\mathbf {Z} |{\boldsymbol {\theta }})\right]:=\int \log p(\mathbf {X} ,\mathbf {Z} |{\boldsymbol {\theta }})\,p(\mathbf {Z} |\mathbf {X} ,{\boldsymbol {\theta }}^{(t)})\,d\mathbf {Z} \,}
最大化ステップ(Mステップ) :この量を最大化するパラメータを見つけます。
θ
(
t
+
1
)
=
a
r
g
m
a
x
θ
Q
(
θ
∣
θ
(
t
)
)
{\displaystyle {\boldsymbol {\theta }}^{(t+1)}={\underset {\boldsymbol {\theta }}{\operatorname {arg\,max} }}\ Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})\,}
もっと簡潔に言えば、これを 1 つの式で表すことができます。
θ
(
t
+
1
)
=
a
r
g
m
a
x
θ
E
Z
∼
p
(
⋅
|
X
,
θ
(
t
)
)
[
log
p
(
X
,
Z
|
θ
)
]
{\displaystyle {\boldsymbol {\theta }}^{(t+1)}={\underset {\boldsymbol {\theta }}{\operatorname {arg\,max} }}\operatorname {E} _{\mathbf {Z} \sim p(\cdot |\mathbf {X} ,{\boldsymbol {\theta }}^{(t)})}\left[\log p(\mathbf {X} ,\mathbf {Z} |{\boldsymbol {\theta }})\right]\,}
変数の解釈
EM が適用される典型的なモデルでは、 一連のグループのいずれかのメンバーシップを示す潜在変数として以下を使用します。
Z
{\displaystyle \mathbf {Z} }
観測データポイントは、 離散的 (有限または可算無限集合内の値を取る)または 連続的 (非可算無限集合内の値を取る)です 。各データポイントには、観測値のベクトルが関連付けられる場合があります。
X
{\displaystyle \mathbf {X} }
欠損 値 ( 潜在変数 とも呼ばれる)は 離散的 であり 、固定数の値から抽出され、観測単位ごとに 1 つの潜在変数を持ちます。
Z
{\displaystyle \mathbf {Z} }
パラメータは連続的であり、すべてのデータ ポイントに関連付けられているパラメータと、潜在変数の特定の値に関連付けられているパラメータ (つまり、対応する潜在変数がその値を持つすべてのデータ ポイントに関連付けられているパラメータ) の 2 種類があります。
ただし、EM を他の種類のモデルに適用することは可能です。
動機は次のとおりです。パラメータの値が わかっている場合、通常、潜在変数の値は、 のすべての可能な値にわたって対数尤度を最大化することで見つけることができます。 これは、 を単純に反復するか、 隠れマルコフモデル に対する ビタビアルゴリズム などのアルゴリズムを介して行います 。逆に、潜在変数の値がわかっている場合 、パラメータの推定値をかなり簡単に見つけることができます。通常は、観測されたデータポイントを関連する潜在変数の値に従ってグループ化し、各グループ内のポイントの値、または値の何らかの関数を平均化するだけで済みます。これは、とが両方 とも 未知の
場合に反復アルゴリズムを使用することを示唆しています。
θ
{\displaystyle {\boldsymbol {\theta }}}
Z
{\displaystyle \mathbf {Z} }
Z
{\displaystyle \mathbf {Z} }
Z
{\displaystyle \mathbf {Z} }
Z
{\displaystyle \mathbf {Z} }
θ
{\displaystyle {\boldsymbol {\theta }}}
θ
{\displaystyle {\boldsymbol {\theta }}}
Z
{\displaystyle \mathbf {Z} }
まず、パラメータを ランダムな値に初期化します。
θ
{\displaystyle {\boldsymbol {\theta }}}
が与えられた場合、 の各可能な値の確率を計算します 。
Z
{\displaystyle \mathbf {Z} }
θ
{\displaystyle {\boldsymbol {\theta }}}
次に、 の計算された値を使用して 、パラメータのより正確な推定値を計算します 。
Z
{\displaystyle \mathbf {Z} }
θ
{\displaystyle {\boldsymbol {\theta }}}
収束するまで手順 2 と 3 を繰り返します。
今説明したアルゴリズムは、コスト関数の局所的最小値に単調に近づきます。
プロパティ
EM反復は観測データ(すなわち周辺)尤度関数を増加させますが、その系列が 最大尤度推定値 に収束するという保証はありません。 多峰性分布 の場合、これはEMアルゴリズムが初期値に応じて観測データ尤度関数の 局所的最大値 に収束する可能性があることを意味します。局所的最大値を回避するための様々なヒューリスティックまたは メタヒューリスティックな アプローチが存在します。例えば、ランダムリスタート ヒルクライミング (複数の異なるランダムな初期推定値から開始する )や シミュレーテッドアニーリング 法の適用などが挙げられます。
θ
(
t
)
{\displaystyle {\boldsymbol {\theta }}^{(t)}}
EMは、尤度が 指数族 である場合に特に有用です。包括的な説明については、Sundberg(2019、第8章)を参照してください。 [16] Eステップは 十分な統計量 の期待値の合計となり 、Mステップは線形関数の最大化を伴います。このような場合、通常、 Sundbergの公式 [17] ( Per Martin-Löf と Anders Martin-Löf の未発表の結果に基づいてRolf Sundbergによって証明および発表された)を使用して、各ステップの 閉じた形式の式の 更新を導くことが可能です。 [8] [9] [11] [12] [13] [14]
EM メソッドは、 Dempster、Laird、および Rubin による元の論文で、
ベイズ推論の 最大事後確率 (MAP) 推定値を計算するように変更されました。
最尤推定値を求める他の方法としては、 勾配降下法 、 共役勾配法 、 ガウス・ニュートン法 の変種などがあります。EMとは異なり、これらの方法では通常、尤度関数の一次導関数および/または二次導関数の評価が必要です。
正しさの証明
期待最大化は、 直接的に改善するのではなく、改善するように機能します 。ここでは、前者の改善が後者の改善を意味することが示されています。 [18]
Q
(
θ
∣
θ
(
t
)
)
{\displaystyle Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})}
log
p
(
X
∣
θ
)
{\displaystyle \log p(\mathbf {X} \mid {\boldsymbol {\theta }})}
確率がゼロでない 任意の に対して 、次のように書くことができる。
Z
{\displaystyle \mathbf {Z} }
p
(
Z
∣
X
,
θ
)
{\displaystyle p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }})}
log
p
(
X
∣
θ
)
=
log
p
(
X
,
Z
∣
θ
)
−
log
p
(
Z
∣
X
,
θ
)
.
{\displaystyle \log p(\mathbf {X} \mid {\boldsymbol {\theta }})=\log p(\mathbf {X} ,\mathbf {Z} \mid {\boldsymbol {\theta }})-\log p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }}).}
現在のパラメータ推定値における 未知データの可能な値に対する期待値を求めるには、 両辺に を掛け 、 について合計(または積分)します 。左辺は定数の期待値なので、以下の式を得ます。
Z
{\displaystyle \mathbf {Z} }
θ
(
t
)
{\displaystyle \theta ^{(t)}}
p
(
Z
∣
X
,
θ
(
t
)
)
{\displaystyle p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }}^{(t)})}
Z
{\displaystyle \mathbf {Z} }
log
p
(
X
∣
θ
)
=
∑
Z
p
(
Z
∣
X
,
θ
(
t
)
)
log
p
(
X
,
Z
∣
θ
)
−
∑
Z
p
(
Z
∣
X
,
θ
(
t
)
)
log
p
(
Z
∣
X
,
θ
)
=
Q
(
θ
∣
θ
(
t
)
)
+
H
(
θ
∣
θ
(
t
)
)
,
{\displaystyle {\begin{aligned}\log p(\mathbf {X} \mid {\boldsymbol {\theta }})&=\sum _{\mathbf {Z} }p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }}^{(t)})\log p(\mathbf {X} ,\mathbf {Z} \mid {\boldsymbol {\theta }})-\sum _{\mathbf {Z} }p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }}^{(t)})\log p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }})\\&=Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})+H({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)}),\end{aligned}}}
ここで 、 はそれが置き換える和の負の数によって定義されます。この最後の式は を含むのあらゆる値に対して成り立ちます 。
H
(
θ
∣
θ
(
t
)
)
{\displaystyle H({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})}
θ
{\displaystyle {\boldsymbol {\theta }}}
θ
=
θ
(
t
)
{\displaystyle {\boldsymbol {\theta }}={\boldsymbol {\theta }}^{(t)}}
log
p
(
X
∣
θ
(
t
)
)
=
Q
(
θ
(
t
)
∣
θ
(
t
)
)
+
H
(
θ
(
t
)
∣
θ
(
t
)
)
,
{\displaystyle \log p(\mathbf {X} \mid {\boldsymbol {\theta }}^{(t)})=Q({\boldsymbol {\theta }}^{(t)}\mid {\boldsymbol {\theta }}^{(t)})+H({\boldsymbol {\theta }}^{(t)}\mid {\boldsymbol {\theta }}^{(t)}),}
この最後の式を前の式から引くと
log
p
(
X
∣
θ
)
−
log
p
(
X
∣
θ
(
t
)
)
=
Q
(
θ
∣
θ
(
t
)
)
−
Q
(
θ
(
t
)
∣
θ
(
t
)
)
+
H
(
θ
∣
θ
(
t
)
)
−
H
(
θ
(
t
)
∣
θ
(
t
)
)
.
{\displaystyle \log p(\mathbf {X} \mid {\boldsymbol {\theta }})-\log p(\mathbf {X} \mid {\boldsymbol {\theta }}^{(t)})=Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})-Q({\boldsymbol {\theta }}^{(t)}\mid {\boldsymbol {\theta }}^{(t)})+H({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})-H({\boldsymbol {\theta }}^{(t)}\mid {\boldsymbol {\theta }}^{(t)}).}
しかし、 ギブスの不等式 によれば 、
H
(
θ
∣
θ
(
t
)
)
≥
H
(
θ
(
t
)
∣
θ
(
t
)
)
{\displaystyle H({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})\geq H({\boldsymbol {\theta }}^{(t)}\mid {\boldsymbol {\theta }}^{(t)})}
log
p
(
X
∣
θ
)
−
log
p
(
X
∣
θ
(
t
)
)
≥
Q
(
θ
∣
θ
(
t
)
)
−
Q
(
θ
(
t
)
∣
θ
(
t
)
)
.
{\displaystyle \log p(\mathbf {X} \mid {\boldsymbol {\theta }})-\log p(\mathbf {X} \mid {\boldsymbol {\theta }}^{(t)})\geq Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})-Q({\boldsymbol {\theta }}^{(t)}\mid {\boldsymbol {\theta }}^{(t)}).}
言い換えれば、 改善する ことを選択すると 、少なくとも同じだけ改善がもたらされます。
θ
{\displaystyle {\boldsymbol {\theta }}}
Q
(
θ
∣
θ
(
t
)
)
{\displaystyle Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})}
log
p
(
X
∣
θ
)
{\displaystyle \log p(\mathbf {X} \mid {\boldsymbol {\theta }})}
最大化-最大化手順として
EMアルゴリズムは、2つの交互の最大化ステップ、つまり 座標降下法 の例として見ることができます。 [19] [20] 次の関数を考えてみましょう。
F
(
q
,
θ
)
:=
E
q
[
log
L
(
θ
;
x
,
Z
)
]
+
H
(
q
)
,
{\displaystyle F(q,\theta ):=\operatorname {E} _{q}[\log L(\theta ;x,Z)]+H(q),}
ここで、 q は観測されていないデータ z 上の任意の確率分布であり、 H(q) は分布 qの エントロピー である 。この関数は次のように書ける。
F
(
q
,
θ
)
=
−
D
K
L
(
q
∥
p
Z
∣
X
(
⋅
∣
x
;
θ
)
)
+
log
L
(
θ
;
x
)
,
{\displaystyle F(q,\theta )=-D_{\mathrm {KL} }{\big (}q\parallel p_{Z\mid X}(\cdot \mid x;\theta ){\big )}+\log L(\theta ;x),}
ここで 、 は観測データが与えられた場合の観測されていないデータの条件付き分布であり 、は Kullback–Leibler 情報 です 。
p
Z
∣
X
(
⋅
∣
x
;
θ
)
{\displaystyle p_{Z\mid X}(\cdot \mid x;\theta )}
x
{\displaystyle x}
D
K
L
{\displaystyle D_{KL}}
EM アルゴリズムのステップは次のように考えられます。
期待ステップ : 最大化するように選択 :
q
{\displaystyle q}
F
{\displaystyle F}
q
(
t
)
=
a
r
g
m
a
x
q
F
(
q
,
θ
(
t
)
)
{\displaystyle q^{(t)}=\operatorname {arg\,max} _{q}\ F(q,\theta ^{(t)})}
最大化ステップ : 最大化を選択 :
θ
{\displaystyle \theta }
F
{\displaystyle F}
θ
(
t
+
1
)
=
a
r
g
m
a
x
θ
F
(
q
(
t
)
,
θ
)
{\displaystyle \theta ^{(t+1)}=\operatorname {arg\,max} _{\theta }\ F(q^{(t)},\theta )}
アプリケーション
フィルタリングとスムージングEMアルゴリズム
カルマンフィルタ は 一般的にオンライン状態推定に使用され、最小分散スムーザはオフラインまたはバッチ状態推定に使用されます。ただし、これらの最小分散解法では、状態空間モデルパラメータの推定値が必要です。EMアルゴリズムは、状態推定とパラメータ推定を同時に行う問題に使用できます。
フィルタリングとスムージングの EM アルゴリズムは、次の 2 段階の手順を繰り返すことによって生成されます。
Eステップ
現在のパラメータ推定値を使用して設計されたカルマン フィルタまたは最小分散スムージングを操作して、更新された状態推定値を取得します。
Mステップ
最大尤度計算内でフィルタリングまたは平滑化された状態推定値を使用して、更新されたパラメータ推定値を取得します。
加法性白色ノイズを含む単入力単出力システムの測定データに カルマンフィルタ または最小分散平滑化フィルタを適用すると仮定する。測定ノイズ分散 推定
値は、 最尤法によって更新することができる。
σ
^
v
2
=
1
N
∑
k
=
1
N
(
z
k
−
x
^
k
)
2
,
{\displaystyle {\widehat {\sigma }}_{v}^{2}={\frac {1}{N}}\sum _{k=1}^{N}{(z_{k}-{\widehat {x}}_{k})}^{2},}
ここで 、はN個のスカラー測定値からフィルタまたは平滑化器によって計算されたスカラー出力推定値である 。上記の更新は、ポアソン測定のノイズ強度の更新にも適用できる。同様に、1次自己回帰プロセスの場合、更新されたプロセスノイズ分散推定値は次のように計算できる。
x
^
k
{\displaystyle {\widehat {x}}_{k}}
z
k
{\displaystyle z_{k}}
σ
^
w
2
=
1
N
∑
k
=
1
N
(
x
^
k
+
1
−
F
^
x
^
k
)
2
,
{\displaystyle {\widehat {\sigma }}_{w}^{2}={\frac {1}{N}}\sum _{k=1}^{N}{({\widehat {x}}_{k+1}-{\widehat {F}}{\widehat {x}}_{k})}^{2},}
ここで 、および はフィルタまたは平滑化器によって計算されるスカラー状態推定値である。更新されたモデル係数推定値は次のように得られる。
x
^
k
{\displaystyle {\widehat {x}}_{k}}
x
^
k
+
1
{\displaystyle {\widehat {x}}_{k+1}}
F
^
=
∑
k
=
1
N
(
x
^
k
+
1
−
F
^
x
^
k
)
2
∑
k
=
1
N
x
^
k
2
.
{\displaystyle {\widehat {F}}={\frac {\sum _{k=1}^{N}{({\widehat {x}}_{k+1}-{\widehat {F}}{\widehat {x}}_{k})}^{2}}{\sum _{k=1}^{N}{\widehat {x}}_{k}^{2}}}.}
上記のようなパラメータ推定値の収束は十分に研究されている。 [26] [27] [28] [29]
変種
EMアルゴリズムの収束が遅い場合があり、これを加速するために共役勾配法 や修正 ニュートン法 (ニュートン・ラプソン法) など、いくつかの手法が提案されている。 [30] また、EMは制約付き推定法と組み合わせて用いることもできる。
パラメータ拡張期待値最大化(PX-EM) アルゴリズムは、多くの場合、「補完された完全データで捕捉された追加情報を活用して、Mステップの分析を修正するための『共分散調整』を使用する」ことで速度向上を実現します。 [31]
期待値条件付き最大化(ECM) は、各Mステップを、他のパラメータが固定されたままの条件付きで各パラメータ θ i が個別に最大化される条件付き最大化(CM)ステップのシーケンスに置き換えます。 [32]これは 期待値条件付き最大化(ECME) アルゴリズムに拡張できます 。 [33]
この考え方は、一般化期待最大化(GEM) アルゴリズムにさらに拡張され、 最大化-最大化手順のセクションで説明されているように、EステップとMステップの両方で 目的関数 Fの増加のみが求められます。 [19] GEMは分散環境でさらに開発され、有望な結果を示しています。 [34]
EMアルゴリズムをMM (コンテキストに応じてMajorize/MinimizeまたはMinorize/Maximize)アルゴリズム のサブクラスと見なすことも可能であり、 [35] したがって、より一般的なケースで開発された任意のメカニズムを使用することもできる。
α-EMアルゴリズム
EMアルゴリズムで使用されるQ関数は対数尤度に基づいているため、log-EMアルゴリズムとみなされる。対数尤度の使用は、α対数尤度比の使用に一般化できる。そして、観測データのα対数尤度比は、α対数尤度比のQ関数とαダイバージェンスを用いることで、等式として正確に表現できる。このQ関数を求めることは一般化されたEステップであり、その最大化は一般化されたMステップである。このペアはα-EMアルゴリズム [36]と呼ばれ、そのサブクラスとしてlog-EMアルゴリズムを含む。したがって、 松山康雄
によるα-EMアルゴリズムは、 log-EMアルゴリズムの正確な一般化である。勾配やヘッセ行列の計算は不要である。α-EMは、適切なαを選択することで、log-EMアルゴリズムよりも収束が速い。α-EMアルゴリズムは、隠れマルコフモデル推定アルゴリズムα-HMMの高速版となる
[37] 。
変分ベイズ法との関係
EM は部分的に非ベイジアンな最大尤度法です。その最終結果は、潜在変数の 確率分布 (ベイジアン スタイル) と θ の点推定値( 最尤推定値 または事後モード) を示します。これを完全にベイジアン化したもので、 θ と潜在変数の確率分布が得られるものが必要な場合があります。ベイジアンによる推論のアプローチは、単に θ を 別の潜在変数として扱うことです。このパラダイムでは、 E ステップと M ステップの区別はなくなります。上記のように因数分解された Q 近似 ( 変分ベイズ ) を使用する場合、解法は各潜在変数 (今度は θ を含む) を反復処理し、一度に 1 つずつ最適化できます。ここで、反復ごとに k ステップが必要です。k は 潜在変数の数です。 グラフィカル モデルの場合、各変数の新しい Q はその マルコフ ブランケット のみに依存する ため、これは簡単に実行でき 、ローカル メッセージ パッシングを 使用して効率的な推論を行うことができます。
幾何学的解釈
情報幾何学 では 、E ステップと M ステップは 、e 接続と m 接続と呼ばれる双対 アフィン接続による射影として解釈されます。 カルバック-ライブラー ダイバージェンス もこれらの用語で理解できます。
例
ガウス混合
ELKI で可視化された人工データにおける k-means法 とEM法の比較 。EM法は分散を用いて正規分布を正確に記述できるのに対し、k-means法はデータを ボロノイセル に分割する。クラスター中心は、より明るい大きなシンボルで示されている。
EMアルゴリズムを用いて2成分ガウス混合モデルを Old Faithful データセットに フィッティングする様子を示すアニメーション 。このアルゴリズムは、ランダムな初期化から収束までを段階的に実行します。
を次元の2つの 多変量正規分布 の 混合 からの独立した観測値 のサンプルとし 、 を 観測の起源となる成分を決定する潜在変数とする。 [20]
x
=
(
x
1
,
x
2
,
…
,
x
n
)
{\displaystyle \mathbf {x} =(\mathbf {x} _{1},\mathbf {x} _{2},\ldots ,\mathbf {x} _{n})}
n
{\displaystyle n}
d
{\displaystyle d}
z
=
(
z
1
,
z
2
,
…
,
z
n
)
{\displaystyle \mathbf {z} =(z_{1},z_{2},\ldots ,z_{n})}
X
i
∣
(
Z
i
=
1
)
∼
N
d
(
μ
1
,
Σ
1
)
{\displaystyle X_{i}\mid (Z_{i}=1)\sim {\mathcal {N}}_{d}({\boldsymbol {\mu }}_{1},\Sigma _{1})}
そして
X
i
∣
(
Z
i
=
2
)
∼
N
d
(
μ
2
,
Σ
2
)
,
{\displaystyle X_{i}\mid (Z_{i}=2)\sim {\mathcal {N}}_{d}({\boldsymbol {\mu }}_{2},\Sigma _{2}),}
どこ
P
(
Z
i
=
1
)
=
τ
1
{\displaystyle \operatorname {P} (Z_{i}=1)=\tau _{1}\,}
そして
P
(
Z
i
=
2
)
=
τ
2
=
1
−
τ
1
.
{\displaystyle \operatorname {P} (Z_{i}=2)=\tau _{2}=1-\tau _{1}.}
目的は、ガウス分布間の 混合 値とそれぞれの平均および共分散を表す未知のパラメータを推定することです。
θ
=
(
τ
,
μ
1
,
μ
2
,
Σ
1
,
Σ
2
)
,
{\displaystyle \theta ={\big (}{\boldsymbol {\tau }},{\boldsymbol {\mu }}_{1},{\boldsymbol {\mu }}_{2},\Sigma _{1},\Sigma _{2}{\big )},}
ここで、不完全データ尤度関数は
L
(
θ
;
x
)
=
∏
i
=
1
n
∑
j
=
1
2
τ
j
f
(
x
i
;
μ
j
,
Σ
j
)
,
{\displaystyle L(\theta ;\mathbf {x} )=\prod _{i=1}^{n}\sum _{j=1}^{2}\tau _{j}\ f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{j},\Sigma _{j}),}
そして、完全データ尤度関数は
L
(
θ
;
x
,
z
)
=
p
(
x
,
z
∣
θ
)
=
∏
i
=
1
n
∏
j
=
1
2
[
f
(
x
i
;
μ
j
,
Σ
j
)
τ
j
]
I
(
z
i
=
j
)
,
{\displaystyle L(\theta ;\mathbf {x} ,\mathbf {z} )=p(\mathbf {x} ,\mathbf {z} \mid \theta )=\prod _{i=1}^{n}\prod _{j=1}^{2}\ [f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{j},\Sigma _{j})\tau _{j}]^{\mathbb {I} (z_{i}=j)},}
または
L
(
θ
;
x
,
z
)
=
exp
{
∑
i
=
1
n
∑
j
=
1
2
I
(
z
i
=
j
)
[
log
τ
j
−
1
2
log
|
Σ
j
|
−
1
2
(
x
i
−
μ
j
)
⊤
Σ
j
−
1
(
x
i
−
μ
j
)
−
d
2
log
(
2
π
)
]
}
,
{\displaystyle L(\theta ;\mathbf {x} ,\mathbf {z} )=\exp \left\{\sum _{i=1}^{n}\sum _{j=1}^{2}\mathbb {I} (z_{i}=j){\big [}\log \tau _{j}-{\tfrac {1}{2}}\log |\Sigma _{j}|-{\tfrac {1}{2}}(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{j})^{\top }\Sigma _{j}^{-1}(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{j})-{\tfrac {d}{2}}\log(2\pi ){\big ]}\right\},}
ここで 、 は 指標関数 であり、は 多変量正規分布の
確率密度関数 です。
I
{\displaystyle \mathbb {I} }
f
{\displaystyle f}
最後の等式では、各 i について、1つの指標 は0に等しく、もう1つの指標は1に等しい。したがって、内部和は1つの項に簡約される。
I
(
z
i
=
j
)
{\displaystyle \mathbb {I} (z_{i}=j)}
Eステップ
パラメータ θ ( t )の現在の推定値を考えると、 ベイズの定理 によって Z i の条件付き分布は τ で重み付けされた 正規 密度 の比例高さとして決定されます。
T
j
,
i
(
t
)
:=
P
(
Z
i
=
j
∣
X
i
=
x
i
;
θ
(
t
)
)
=
τ
j
(
t
)
f
(
x
i
;
μ
j
(
t
)
,
Σ
j
(
t
)
)
τ
1
(
t
)
f
(
x
i
;
μ
1
(
t
)
,
Σ
1
(
t
)
)
+
τ
2
(
t
)
f
(
x
i
;
μ
2
(
t
)
,
Σ
2
(
t
)
)
.
{\displaystyle T_{j,i}^{(t)}:=\operatorname {P} (Z_{i}=j\mid X_{i}=\mathbf {x} _{i};\theta ^{(t)})={\frac {\tau _{j}^{(t)}\ f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{j}^{(t)},\Sigma _{j}^{(t)})}{\tau _{1}^{(t)}\ f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{1}^{(t)},\Sigma _{1}^{(t)})+\tau _{2}^{(t)}\ f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{2}^{(t)},\Sigma _{2}^{(t)})}}.}
これらは「メンバーシップ確率」と呼ばれ、通常は E ステップの出力と見なされます (ただし、これは以下の Q 関数ではありません)。
この E ステップは、Q に対してこの関数を設定することに対応します。
Q
(
θ
∣
θ
(
t
)
)
=
E
Z
∣
X
=
x
;
θ
(
t
)
[
log
L
(
θ
;
x
,
Z
)
]
=
E
Z
∣
X
=
x
;
θ
(
t
)
[
log
∏
i
=
1
n
L
(
θ
;
x
i
,
Z
i
)
]
=
E
Z
∣
X
=
x
;
θ
(
t
)
[
∑
i
=
1
n
log
L
(
θ
;
x
i
,
Z
i
)
]
=
∑
i
=
1
n
E
Z
i
∣
X
i
=
x
i
;
θ
(
t
)
[
log
L
(
θ
;
x
i
,
Z
i
)
]
=
∑
i
=
1
n
∑
j
=
1
2
P
(
Z
i
=
j
∣
X
i
=
x
i
;
θ
(
t
)
)
log
L
(
θ
j
;
x
i
,
j
)
=
∑
i
=
1
n
∑
j
=
1
2
T
j
,
i
(
t
)
[
log
τ
j
−
1
2
log
|
Σ
j
|
−
1
2
(
x
i
−
μ
j
)
⊤
Σ
j
−
1
(
x
i
−
μ
j
)
−
d
2
log
(
2
π
)
]
.
{\displaystyle {\begin{aligned}Q(\theta \mid \theta ^{(t)})&=\operatorname {E} _{\mathbf {Z} \mid \mathbf {X} =\mathbf {x} ;\mathbf {\theta } ^{(t)}}[\log L(\theta ;\mathbf {x} ,\mathbf {Z} )]\\&=\operatorname {E} _{\mathbf {Z} \mid \mathbf {X} =\mathbf {x} ;\mathbf {\theta } ^{(t)}}[\log \prod _{i=1}^{n}L(\theta ;\mathbf {x} _{i},Z_{i})]\\&=\operatorname {E} _{\mathbf {Z} \mid \mathbf {X} =\mathbf {x} ;\mathbf {\theta } ^{(t)}}[\sum _{i=1}^{n}\log L(\theta ;\mathbf {x} _{i},Z_{i})]\\&=\sum _{i=1}^{n}\operatorname {E} _{Z_{i}\mid X_{i}=x_{i};\mathbf {\theta } ^{(t)}}[\log L(\theta ;\mathbf {x} _{i},Z_{i})]\\&=\sum _{i=1}^{n}\sum _{j=1}^{2}P(Z_{i}=j\mid X_{i}=\mathbf {x} _{i};\theta ^{(t)})\log L(\theta _{j};\mathbf {x} _{i},j)\\&=\sum _{i=1}^{n}\sum _{j=1}^{2}T_{j,i}^{(t)}{\big [}\log \tau _{j}-{\tfrac {1}{2}}\log |\Sigma _{j}|-{\tfrac {1}{2}}(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{j})^{\top }\Sigma _{j}^{-1}(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{j})-{\tfrac {d}{2}}\log(2\pi ){\big ]}.\end{aligned}}}
和の中の の期待値は 、確率密度関数 に関して求められますが、この 関数は訓練セット ごとに異なる場合があります 。Eステップでは、 を除き 、ステップ実行前に全てが既知です。 は、Eステップセクションの冒頭にある式に従って計算されます。
log
L
(
θ
;
x
i
,
Z
i
)
{\displaystyle \log L(\theta ;\mathbf {x} _{i},Z_{i})}
P
(
Z
i
∣
X
i
=
x
i
;
θ
(
t
)
)
{\displaystyle P(Z_{i}\mid X_{i}=\mathbf {x} _{i};\theta ^{(t)})}
x
i
{\displaystyle \mathbf {x} _{i}}
T
j
,
i
{\displaystyle T_{j,i}}
この完全な条件付き期待値は、 τ と μ / Σ が 別々の線形項で表示され、独立して最大化できる
ため、1 つのステップで計算する必要はありません。
Mステップ
Q
(
θ
∣
θ
(
t
)
)
{\displaystyle Q(\theta \mid \theta ^{(t)})}
が二次形式であるため、 を最大化する値を決定するのは 比較的簡単です。また、 、 、 は すべて別々の線形項として現れるため、それぞれ独立に最大化できます。
θ
{\displaystyle \theta }
τ
{\displaystyle \tau }
(
μ
1
,
Σ
1
)
{\displaystyle ({\boldsymbol {\mu }}_{1},\Sigma _{1})}
(
μ
2
,
Σ
2
)
{\displaystyle ({\boldsymbol {\mu }}_{2},\Sigma _{2})}
まず、 制約 を持つ を考えます 。
τ
{\displaystyle \tau }
τ
1
+
τ
2
=
1
{\displaystyle \tau _{1}+\tau _{2}=1}
τ
(
t
+
1
)
=
a
r
g
m
a
x
τ
Q
(
θ
∣
θ
(
t
)
)
=
a
r
g
m
a
x
τ
{
[
∑
i
=
1
n
T
1
,
i
(
t
)
]
log
τ
1
+
[
∑
i
=
1
n
T
2
,
i
(
t
)
]
log
τ
2
}
.
{\displaystyle {\begin{aligned}{\boldsymbol {\tau }}^{(t+1)}&={\underset {\boldsymbol {\tau }}{\operatorname {arg\,max} }}\ Q(\theta \mid \theta ^{(t)})\\&={\underset {\boldsymbol {\tau }}{\operatorname {arg\,max} }}\ \left\{\left[\sum _{i=1}^{n}T_{1,i}^{(t)}\right]\log \tau _{1}+\left[\sum _{i=1}^{n}T_{2,i}^{(t)}\right]\log \tau _{2}\right\}.\end{aligned}}}
これは二項分布 の最大尤度推定値と同じ形をしており 、
τ
j
(
t
+
1
)
=
∑
i
=
1
n
T
j
,
i
(
t
)
∑
i
=
1
n
(
T
1
,
i
(
t
)
+
T
2
,
i
(
t
)
)
=
1
n
∑
i
=
1
n
T
j
,
i
(
t
)
.
{\displaystyle \tau _{j}^{(t+1)}={\frac {\sum _{i=1}^{n}T_{j,i}^{(t)}}{\sum _{i=1}^{n}(T_{1,i}^{(t)}+T_{2,i}^{(t)})}}={\frac {1}{n}}\sum _{i=1}^{n}T_{j,i}^{(t)}.}
次の推定値については :
(
μ
1
,
Σ
1
)
{\displaystyle ({\boldsymbol {\mu }}_{1},\Sigma _{1})}
(
μ
1
(
t
+
1
)
,
Σ
1
(
t
+
1
)
)
=
a
r
g
m
a
x
μ
1
,
Σ
1
Q
(
θ
∣
θ
(
t
)
)
=
a
r
g
m
a
x
μ
1
,
Σ
1
∑
i
=
1
n
T
1
,
i
(
t
)
{
−
1
2
log
|
Σ
1
|
−
1
2
(
x
i
−
μ
1
)
⊤
Σ
1
−
1
(
x
i
−
μ
1
)
}
.
{\displaystyle {\begin{aligned}({\boldsymbol {\mu }}_{1}^{(t+1)},\Sigma _{1}^{(t+1)})&={\underset {{\boldsymbol {\mu }}_{1},\Sigma _{1}}{\operatorname {arg\,max} }}\ Q(\theta \mid \theta ^{(t)})\\&={\underset {{\boldsymbol {\mu }}_{1},\Sigma _{1}}{\operatorname {arg\,max} }}\ \sum _{i=1}^{n}T_{1,i}^{(t)}\left\{-{\tfrac {1}{2}}\log |\Sigma _{1}|-{\tfrac {1}{2}}(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{1})^{\top }\Sigma _{1}^{-1}(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{1})\right\}\end{aligned}}.}
これは正規分布の重み付き最尤推定と同じ形をしており、
μ
1
(
t
+
1
)
=
∑
i
=
1
n
T
1
,
i
(
t
)
x
i
∑
i
=
1
n
T
1
,
i
(
t
)
{\displaystyle {\boldsymbol {\mu }}_{1}^{(t+1)}={\frac {\sum _{i=1}^{n}T_{1,i}^{(t)}\mathbf {x} _{i}}{\sum _{i=1}^{n}T_{1,i}^{(t)}}}}
そして
Σ
1
(
t
+
1
)
=
∑
i
=
1
n
T
1
,
i
(
t
)
(
x
i
−
μ
1
(
t
+
1
)
)
(
x
i
−
μ
1
(
t
+
1
)
)
⊤
∑
i
=
1
n
T
1
,
i
(
t
)
{\displaystyle \Sigma _{1}^{(t+1)}={\frac {\sum _{i=1}^{n}T_{1,i}^{(t)}(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{1}^{(t+1)})(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{1}^{(t+1)})^{\top }}{\sum _{i=1}^{n}T_{1,i}^{(t)}}}}
そして対称性により、
μ
2
(
t
+
1
)
=
∑
i
=
1
n
T
2
,
i
(
t
)
x
i
∑
i
=
1
n
T
2
,
i
(
t
)
{\displaystyle {\boldsymbol {\mu }}_{2}^{(t+1)}={\frac {\sum _{i=1}^{n}T_{2,i}^{(t)}\mathbf {x} _{i}}{\sum _{i=1}^{n}T_{2,i}^{(t)}}}}
そして
Σ
2
(
t
+
1
)
=
∑
i
=
1
n
T
2
,
i
(
t
)
(
x
i
−
μ
2
(
t
+
1
)
)
(
x
i
−
μ
2
(
t
+
1
)
)
⊤
∑
i
=
1
n
T
2
,
i
(
t
)
.
{\displaystyle \Sigma _{2}^{(t+1)}={\frac {\sum _{i=1}^{n}T_{2,i}^{(t)}(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{2}^{(t+1)})(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{2}^{(t+1)})^{\top }}{\sum _{i=1}^{n}T_{2,i}^{(t)}}}.}
終了
事前に設定されたしきい値を下回る 場合、反復プロセスを終了します 。
E
Z
∣
θ
(
t
)
,
x
[
log
L
(
θ
(
t
)
;
x
,
Z
)
]
≤
E
Z
∣
θ
(
t
−
1
)
,
x
[
log
L
(
θ
(
t
−
1
)
;
x
,
Z
)
]
+
ε
{\displaystyle E_{Z\mid \theta ^{(t)},\mathbf {x} }[\log L(\theta ^{(t)};\mathbf {x} ,\mathbf {Z} )]\leq E_{Z\mid \theta ^{(t-1)},\mathbf {x} }[\log L(\theta ^{(t-1)};\mathbf {x} ,\mathbf {Z} )]+\varepsilon }
ε
{\displaystyle \varepsilon }
一般化
上記のアルゴリズムは、2 つ以上の 多変量正規分布 の混合に対して一般化できます。
打ち切り回帰分析
EMアルゴリズムは 、ある量の変動を説明する 線形回帰モデルが存在するものの、実際に観測された値はモデルで表現された値の打ち切りまたは切り捨てられたバージョンである場合に実装されています。 [38]このモデルの特殊なケースとしては、1つの 正規分布 からの打ち切りまたは切り捨てられた観測値が含まれます 。 [38]
代替案
EMは典型的には局所最適解に収束しますが、必ずしも大域最適解には収束しません。収束速度には一般に上限がありません。高次元では収束性が任意に低下する可能性があり、局所最適解の数が指数関数的に増加する可能性があります。そのため、特に高次元環境において、学習を保証する代替手法が求められています。EMの代替手法として、一貫性の保証がより高く、 モーメントベースアプローチ [39] またはいわゆる スペクトル手法 [40] [41] と呼ばれる手法があります。 確率モデルのパラメータを学習するモーメントベースアプローチは、局所最適解に陥ってしまう問題に悩まされることが多いEMとは異なり、特定の条件下での大域収束などの保証があります。学習を保証するアルゴリズムは、混合モデルやHMMなどの多くの重要なモデルに対して導出可能です。これらのスペクトル手法では、偽の局所最適解は発生せず、いくつかの正則性条件下で真のパラメータを一貫して推定できます。 [ 要出典 ]
参照
参考文献
^ Meng, X.-L.; van Dyk, D. (1997). 「EMアルゴリズム ― 古き良きフォークソングを新しいテンポで歌い上げる」 J. Royal Statist. Soc. B. 59 ( 3): 511– 567. doi : 10.1111/1467-9868.00082 . S2CID 17461647.
^ ジョンヨル・クォン、コンスタンティン・カラマニス『
第23回国際人工知能・統計会議議事録』 、PMLR 108:1727-1736、2020年。
^ Dempster, AP ; Laird, NM ; Rubin, DB (1977). 「EMアルゴリズムによる不完全データからの最大尤度」. Journal of the Royal Statistical Society, Series B. 39 ( 1): 1– 38. doi :10.1111/j.2517-6161.1977.tb01600.x. JSTOR 2984875. MR 0501537.
^ Ceppelini, RM (1955). 「ランダム交配集団における遺伝子頻度の推定」. Ann. Hum. Genet . 20 (2): 97– 115. doi :10.1111/j.1469-1809.1955.tb01360.x. PMID 13268982. S2CID 38625779.
^ ハートレー、ハーマン・オットー (1958). 「不完全データからの最大尤度推定」. バイオメトリクス . 14 (2): 174– 194. doi :10.2307/2527783. JSTOR 2527783.
^ Ng, Shu Kay; Krishnan, Thriyambakam; McLachlan, Geoffrey J. (2011-12-21)、「EMアルゴリズム」、 計算統計ハンドブック 、ベルリン、ハイデルベルク:Springer Berlin Heidelberg、pp. 139– 172、 doi :10.1007/978-3-642-21551-3_6、 ISBN 978-3-642-21550-6 、 S2CID 59942212
^ Sundberg, Rolf (1974). 「指数関数族からの不完全データに対する最大尤度理論」. Scandinavian Journal of Statistics . 1 (2): 49– 58. JSTOR 4615553. MR 0381110.
^ ab
Rolf Sundberg. 1971. 「指数族変数の関数を観測する際に生成される分布に対する最大尤度理論と応用」 ストックホルム大学数理統計研究所博士論文。
^ ab Sundberg, Rolf (1976). 「指数族からの不完全データに対する尤度方程式の反復解法」. Communications in Statistics – Simulation and Computation . 5 (1): 55– 64. doi :10.1080/03610917608812007. MR 0443190.
^ 3、5、11ページのDempster、Laird、Rubinによる謝辞を参照。
^ ab Per Martin-Löf . 1966. 統計力学の観点から見た統計学 . 講義ノート、オーフス大学数学研究所。(「Sundbergの公式」はAnders Martin-Löfによる)。
^ ab Per Martin-Löf 。 1970. Statistiska Modeller (統計モデル): Anteckningar från seminarier läsåret 1969–1970 (講義ノート 1969-1970)、Rolf Sundberg の協力を得て。 ストックホルム大学。
^ ab Martin-Löf, P. 冗長性の概念と、統計仮説と観測データ集合間の偏差を定量的に測る尺度としてのその利用。F. Abildgård、 A.P. Dempster 、 D. Basu 、 D.R. Cox 、 AWF Edwards 、D.A. Sprott、 G.A. Barnard 、O. Barndorff-Nielsen、J.D. Kalbfleisch、 G. Rasch による議論と著者による返答付き。 統計推論における基礎的問題に関する会議議事録 (オーフス、1973年)、1~42頁。回顧録、第1号、オーフス大学数学研究所理論統計学科、オーフス、1974年。
^ ab Martin-Löf, Per (1974). 「冗長性の概念と、統計仮説と観測データ集合間の乖離を定量的に測る尺度としてのその利用」 Scand. J. Statist . 1 (1): 3– 18.
^ abc Wu, CF Jeff (1983年3月). 「EMアルゴリズムの収束特性について」 Annals of Statistics . 11 (1): 95–103 . doi : 10.1214/aos/1176346060 . JSTOR 2240463. MR 0684867.
^ サンドバーグ、ロルフ(2019年) 『指数族による統計モデリング』 ケンブリッジ大学出版局、 ISBN 978-1-108-70111-2 。
^ レアード、ナン (2006). 「サンドバーグの公式」. 統計科学百科事典 . Wiley. doi :10.1002/0471667196.ess2643.pub2. ISBN 0-471-66719-6 。
^ リトル、ロデリック JA; ルービン、ドナルド B. (1987). 欠損データを含む統計分析 . ワイリー確率・数理統計シリーズ. ニューヨーク: ジョン・ワイリー・アンド・サンズ. pp. 134–136. ISBN 978-0-471-80254-9 。
^ ab Neal, Radford; Hinton, Geoffrey (1999). 「増分法、スパース法、その他の変種を正当化するEMアルゴリズムの視点」 Michael I. Jordan (編). 『グラフィカルモデルにおける学習』 (PDF) . マサチューセッツ州ケンブリッジ: MIT Press. pp. 355– 368. ISBN 978-0-262-60032-3 . 2009年3月22日 閲覧 。
^ ab ハスティー, トレバー ; ティブシラニ, ロバート ; フリードマン, ジェローム (2001). 「8.5 EMアルゴリズム」. 統計学習の要素 . ニューヨーク: シュプリンガー. pp. 236–243. ISBN 978-0-387-95284-0 。
^ リンドストロム, メアリー J; ベイツ, ダグラス M (1988). 「反復測定データのための線形混合効果モデルのためのニュートン-ラプソン法とEMアルゴリズム」 アメリカ統計学会誌 . 83 (404): 1014. doi :10.1080/01621459.1988.10478693.
^ Van Dyk, David A (2000). 「効率的なEM型アルゴリズムを用いた混合効果モデルのフィッティング」. 計算およびグラフィカル統計ジャーナル . 9 (1): 78– 98. doi :10.2307/1390614. JSTOR 1390614.
^ Diffey, S. M; Smith, A. B; Welsh, A. H; Cullis, B. R (2017). 「線形混合モデルのための新しいREML(パラメータ拡張)EMアルゴリズム」. Australian & New Zealand Journal of Statistics . 59 (4): 433. doi : 10.1111/anzs.12208 . hdl : 1885/211365 .
^ Matarazzo, TJ, Pakzad, SN (2016). 「期待値最大化を用いた構造同定のためのSTRIDE:モーダル同定のための反復出力のみの手法」Journal of Engineering Mechanics.http://ascelibrary.org/doi/abs/10.1061/(ASCE)EM.1943-7889.0000951
^ Kreer, Markus; Kizilersu, Ayse; Thomas, Anthony W. (2022). 「混合体に対する打ち切り期待値最大化アルゴリズム:トレード間待ち時間への応用」 . Physica A: 統計力学とその応用 . 587 (1) 126456. Bibcode :2022PhyA..58726456K. doi :10.1016/j.physa.2021.126456. ISSN 0378-4371. S2CID 244198364.
^ Einicke, GA; Malos, JT; Reid, DC; Hainsworth, DW (2009年1月). 「慣性航法アライメントのためのリカッチ方程式とEMアルゴリズムの収束」. IEEE Trans. Signal Process . 57 (1): 370– 375. Bibcode :2009ITSP...57..370E. doi :10.1109/TSP.2008.2007090. S2CID 1930004.
^ Einicke, GA; Falco, G.; Malos, JT (2010年5月). 「ナビゲーションのためのEMアルゴリズム状態行列推定」. IEEE Signal Processing Letters . 17 (5): 437– 440. Bibcode :2010ISPL...17..437E. doi :10.1109/LSP.2010.2043151. S2CID 14114266.
^ Einicke, GA; Falco, G.; Dunn, MT; Reid, DC (2012年5月). 「反復スムーザーベースの分散推定」. IEEE Signal Processing Letters . 19 (5): 275– 278. Bibcode :2012ISPL...19..275E. doi :10.1109/LSP.2012.2190278. S2CID 17476971.
^ Einicke, GA (2015年9月). 「ポアソンノイズを含む測定値の反復フィルタリングと平滑化」. IEEE Transactions on Aerospace and Electronic Systems . 51 (3): 2205– 2011. Bibcode :2015ITAES..51.2205E. doi :10.1109/TAES.2015.140843. S2CID 32667132.
^ ジャムシディアン、モルタザ;ロバート I ジェンリッヒ (1997)。 「準ニュートン法を使用した EM アルゴリズムの高速化」。 王立統計協会ジャーナル、シリーズ B 。 59 (2): 569–587 。 土井 :10.1111/1467-9868.00083。 MR 1452026。S2CID 121966443 。
^ Liu, C (1998). 「EMを加速するためのパラメータ拡張:PX-EMアルゴリズム」. Biometrika . 85 (4): 755– 770. CiteSeerX 10.1.1.134.9617 . doi :10.1093/biomet/85.4.755.
^ Meng, Xiao-Li; Rubin, Donald B. (1993). 「ECMアルゴリズムによる最大尤度推定:一般的な枠組み」. Biometrika . 80 (2): 267– 278. doi :10.1093/biomet/80.2.267. MR 1243503. S2CID 40571416.
^ Liu, Chuanhai; Rubin, Donald B (1994). 「ECMEアルゴリズム:EMとECMのシンプルな拡張による高速単調収束」 Biometrika . 81 (4): 633. doi :10.1093/biomet/81.4.633. JSTOR 2337067.
^ Jiangtao Yin、Yanfeng Zhang、Lixin Gao (2012). 「頻繁な更新による期待値最大化アルゴリズムの高速化」 (PDF) . IEEE International Conference on Cluster Computing の議事録 .
^ Hunter DRとLange K(2004)、MMアルゴリズムのチュートリアル、アメリカ統計学者、58:30–37
^ 松山康夫 (2003). 「α-EMアルゴリズム:α対数情報量を用いた代理尤度最大化」. IEEE Transactions on Information Theory . 49 (3): 692– 706. doi :10.1109/TIT.2002.808105.
^ 松山康夫 (2011). 「α-EMアルゴリズムに基づく隠れマルコフモデル推定:離散型および連続型α-HMM」. 国際ニューラルネットワーク合同会議 : 808–816 .
^ ab Wolynetz, MS (1979). 「限定的かつ打ち切り正規分布データを用いた線形モデルによる最大尤度推定」. Journal of the Royal Statistical Society, Series C. 28 ( 2): 195– 206. doi :10.2307/2346749. JSTOR 2346749.
^ ピアソン, カール (1894). 「進化の数学的理論への貢献」. ロンドン王立協会哲学論文集 A. 185 : 71–110 . Bibcode : 1894RSPTA.185...71P. doi : 10.1098/rsta.1894.0003 . ISSN 0264-3820. JSTOR 90667.
^ Shaban, Amirreza; Mehrdad, Farajtabar; Bo, Xie; Le, Song; Byron, Boots (2015). 「外点法を用いたスペクトル解の改善による潜在変数モデルの学習」 (PDF) . UAI : 792–801 . 2016年12月24日時点の オリジナル (PDF)からアーカイブ。 2019年6月12日 閲覧 。
^ バレ、ボルハ・クアットーニ、アリアドナ・カレーラス、ザビエル (2012-06-27)。 演算子モデルにおける局所損失の最適化: スペクトル学習への新たな洞察 。 OCLC 815865081。 {{cite book }}: CS1 maint: multiple names: authors list (link )
^ ランゲ、ケネス. 「MMアルゴリズム」 (PDF) .
さらに読む
ホッグ, ロバート; マッキーン, ジョセフ; クレイグ, アレン (2005). 『数理統計学入門 』 アッパーサドルリバー, ニュージャージー州: ピアソン・プレンティス・ホール. pp. 359– 364.
デラート, フランク (2002年2月). 期待値最大化アルゴリズム (PDF) (技術レポート番号 GIT-GVU-02-20). ジョージア工科大学コンピューティング学部. 下限最大化に関する EM アルゴリズムをより簡単に説明します。
ビショップ、クリストファー・M. (2006). パターン認識と機械学習 . シュプリンガー. ISBN 978-0-387-31073-2 。
Gupta, MR; Chen, Y. (2010). 「EMアルゴリズムの理論と利用」. 信号処理の基礎と動向 . 4 (3): 223– 296. CiteSeerX 10.1.1.219.6830 . doi :10.1561/2000000034. GMM、HMM、ディリクレの EM の詳細な導出を含む、EM に関するよく書かれた短い本です。
ビルメス、ジェフ (1997). EMアルゴリズムとそのガウス混合モデルおよび隠れマルコフモデルのパラメータ推定への応用に関する分かりやすいチュートリアル (技術レポート TR-97-021). 国際コンピュータサイエンス研究所. ガウス混合分布とガウス混合隠れマルコフモデルの EM 方程式の簡略化された導出が含まれています。
マクラクラン, ジェフリー・J.; クリシュナン, スリヤンバカム (2008). 『EMアルゴリズムとその拡張 (第2版)』 ホーボーケン: ワイリー. ISBN 978-0-471-20170-0 。
外部リンク
SOCR アクティビティとアプレットの組み合わせには、EMと混合モデリングを組み合わせた様々な1D、2D、3Dデモンストレーションが用意されています 。これらのアプレットとアクティビティは、多様な設定におけるパラメータ推定におけるEMアルゴリズムの特性を経験的に示します。
ガウス混合分布を含むC++ のクラス階層 (GPL)
オンライン教科書「情報理論、推論、学習アルゴリズム」( David JC MacKay 著) には、ソフト k 平均法アルゴリズムを使用したクラスタリングなどの EM アルゴリズムの簡単な例が含まれており、バージョン 7.2 (第 4 版) の第 33.7 章で説明されているように、EM アルゴリズムの変分ビューが強調されています。
MJ Beal 著の「近似ベイズ推論のための変分アルゴリズム」には、EM と変分ベイズ EM の比較と、変分ベイズ HMM (章) を含むいくつかのモデルの導出が含まれています。
期待値最大化アルゴリズム: 短いチュートリアル、Sean Borman による EM アルゴリズムの自己完結型導出。
EM アルゴリズム、Xiaojin Zhu 著。
EMアルゴリズムとその変種:Alexis Rocheによる非公式チュートリアル。EMとその多くの興味深い変種について、簡潔かつ非常に分かりやすく解説しています。