期待値最大化アルゴリズム

統計学において、期待値最大化( EM )アルゴリズムは、統計モデルにおけるパラメータの(局所的)最大尤度または最大事後確率 (MAP) 推定値を求める反復法であり、モデルは観測されない潜在変数に依存している。^[1] EM 反復は、パラメータの現在の推定値を使用して評価された対数尤度の期待値に対する関数を作成する期待値 (E) ステップと、E ステップで見つかった期待対数尤度を最大化するパラメータを計算する最大化 (M) ステップを交互に実行する。これらのパラメータ推定値は、次のEステップで潜在変数の分布を決定するために使用される。たとえば、ガウス分布の混合を推定したり、多重線形回帰問題を解いたりするために使用できる。 ^[2]

EMアルゴリズムは、1977年のアーサー・デンプスター、ナン・レアード、ドナルド・ルービンによる古典的な論文で説明され、その名前が付けられました。^[3]彼らは、この方法が以前の著者によって「特別な状況で何度も提案された」ことを指摘しました。最も初期のものの一つは、セドリック・スミスによる対立遺伝子頻度を推定するための遺伝子計数法です。^{[4]もう1つは、1958年に}HOハートリーによって、1977年にハートリーとホッキングによって提案され、デンプスター・レアード・ルービンの論文のアイデアの多くはこの方法から生まれました。^[5]もう1つは、1977年にSK・ン、スリヤムバカム・クリシュナン、GJ・マクラクランによって提案されました。 ^[6]ハートリーのアイデアは、グループ化された任意の離散分布に拡張できます。指数族に対するEM法の非常に詳細な扱いは、Rolf Sundbergによって学位論文といくつかの論文^{[7] [8] [9]}で発表されました。これは、 Per Martin-LöfおよびAnders Martin-Löfとの共同研究に続いて行われました。 ^{[10] [11] [12] [13] [14]} 1977年のDempster–Laird–Rubinの論文では、この方法が一般化され、より広範な問題に対する収束分析の概要が示されました。Dempster–Laird–Rubinの論文は、EM法を統計分析の重要なツールとして確立しました。Mengとvan Dyk（1997）も参照してください。

デンプスター・レアード・ルービン法の収束解析には欠陥があり、正しい収束解析は1983年にCFジェフ・ウーによって発表された。 ^[15] ウーの証明は、デンプスター・レアード・ルービンが主張したように、指数族の外でもEM法の収束を確立した。^[15]

EM アルゴリズムは、方程式を直接解くことができない場合に、統計モデルの(ローカル)最大尤度パラメータを見つけるために使用されます。通常、これらのモデルには、未知のパラメータと既知のデータ観測に加えて、潜在変数が含まれます。つまり、データの中に欠損値が存在するか、または、さらに観測されていないデータ ポイントが存在すると仮定することで、モデルをより簡単に定式化できます。たとえば、混合モデルは、各観測データ ポイントに対応する観測されていないデータ ポイント、つまり各データ ポイントが属する混合コンポーネントを指定する潜在変数があると仮定することで、より簡単に説明できます。

最大尤度解を求めるには、通常、すべての未知の値、パラメータ、潜在変数に関して尤度関数の導関数を取り、同時に結果の方程式を解く必要があります。潜在変数を含む統計モデルでは、これは通常不可能です。代わりに、結果は通常、パラメータの解には潜在変数の値が必要であり、その逆もまた同様である、連動する方程式のセットになりますが、1 つの方程式のセットを他の方程式に代入すると、解けない方程式が生成されます。

EM アルゴリズムは、これら 2 組の方程式を数値的に解く方法があるという観察から出発します。2 組の未知数のうちの 1 組に任意の値を選択し、それを使用して 2 組目の未知数を推定し、これらの新しい値を使用して 1 組目のより適切な推定値を見つけ、結果の値が両方とも固定点に収束するまで 2 組を交互に繰り返します。これが機能するかどうかは明らかではありませんが、このコンテキストでは証明できます。さらに、その時点で尤度の導関数が (任意に) ゼロであることを証明できます。これは、その点が局所的最大値または鞍点のいずれかであることを意味します。^[15]一般に、複数の最大値が発生する可能性があり、グローバル最大値が見つかる保証はありません。一部の尤度には特異点、つまり無意味な最大値もあります。たとえば、混合モデルで EM によって見つかる可能性があるソリューションの 1 つは、コンポーネントの 1 つをゼロ分散に設定し、同じコンポーネントの平均パラメーターをデータ ポイントの 1 つと等しくすることです。

観測データの集合、観測されていない潜在データまたは欠損値、および未知のパラメータのベクトルを生成する統計モデルと尤度関数が与えられた場合、未知のパラメータの最大尤度推定値（MLE）は、観測データの 周辺尤度を最大化することによって決定される。${\displaystyle \mathbf {X} }$ ${\displaystyle \mathbf {Z} }$${\displaystyle {\boldsymbol {\theta }}}$ ${\displaystyle L({\boldsymbol {\theta }};\mathbf {X} ,\mathbf {Z} )=p(\mathbf {X} ,\mathbf {Z} \mid {\boldsymbol {\theta }})}$

${\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} }$

しかし、は観測されておらず、 に達するまでは の分布は不明であるため、この量は扱いにくいことがよくあります。 ${\displaystyle \mathbf {Z} }$${\displaystyle \mathbf {Z} }$${\displaystyle {\boldsymbol {\theta }}}$

EM アルゴリズムは、次の 2 つのステップを繰り返し適用して、周辺尤度の最大尤度推定値を見つけようとします。

期待ステップ（Eステップ）：を、与えられたの現在の条件付き分布とパラメータ の現在の推定値に関して、の対数尤度関数の期待値として定義します。 ${\displaystyle Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})}$${\displaystyle {\boldsymbol {\theta }}}$${\displaystyle \mathbf {Z} }$${\displaystyle \mathbf {X} }$${\displaystyle {\boldsymbol {\theta }}^{(t)}}$

${\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]\,}$

最大化ステップ（Mステップ）：この量を最大化するパラメータを見つけます。

${\displaystyle {\boldsymbol {\theta }}^{(t+1)}={\underset {\boldsymbol {\theta }}{\operatorname {arg\,max} }}\ Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})\,}$

もっと簡潔に言えば、これを 1 つの式で表すことができます。$${\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 が適用される典型的なモデルでは、一連のグループのいずれかのメンバーシップを示す潜在変数として以下を使用します。 ${\displaystyle \mathbf {Z} }$

1. 観測データ ポイントは、離散的(有限または可算無限集合内の値を取る) または連続的(不可算無限集合内の値を取る)です。各データ ポイントには、観測のベクトルが関連付けられる場合があります。${\displaystyle \mathbf {X} }$
2. 欠損値（潜在変数とも呼ばれる）は離散的であり、固定数の値から抽出され、観測単位ごとに 1 つの潜在変数を持ちます。${\displaystyle \mathbf {Z} }$
3. パラメータは連続しており、すべてのデータ ポイントに関連付けられているパラメータと、潜在変数の特定の値に関連付けられているパラメータ (つまり、対応する潜在変数がその値を持つすべてのデータ ポイントに関連付けられているパラメータ) の 2 種類があります。

ただし、EM を他の種類のモデルに適用することは可能です。

動機は次のとおりです。パラメータの値がわかっている場合、潜在変数の値は、のすべての可能な値にわたって対数尤度を最大化することで通常見つけることができます。これは、を反復するか、隠れマルコフモデルのViterbi アルゴリズムなどのアルゴリズムを介して行うだけです。逆に、潜在変数の値がわかっている場合は、通常、観測されたデータ ポイントを関連する潜在変数の値に従ってグループ化し、各グループ内のポイントの値または値の関数を平均するだけで、パラメータの推定値をかなり簡単に見つけることができます。これは、と の両方が不明な場合に反復アルゴリズムを示唆しています。 ${\displaystyle {\boldsymbol {\theta }}}$${\displaystyle \mathbf {Z} }$${\displaystyle \mathbf {Z} }$${\displaystyle \mathbf {Z} }$${\displaystyle \mathbf {Z} }$${\displaystyle {\boldsymbol {\theta }}}$${\displaystyle {\boldsymbol {\theta }}}$${\displaystyle \mathbf {Z} }$

1. まず、パラメータをランダムな値に初期化します。${\displaystyle {\boldsymbol {\theta }}}$
2. が与えられた場合、の各可能な値の確率を計算します。${\displaystyle \mathbf {Z} }$${\displaystyle {\boldsymbol {\theta }}}$
3. 次に、計算したばかりの値を使用して、パラメータのより良い推定値を計算します。${\displaystyle \mathbf {Z} }$${\displaystyle {\boldsymbol {\theta }}}$
4. 収束するまで手順 2 と 3 を繰り返します。

今説明したアルゴリズムは、コスト関数の局所最小値に単調に近づきます。

EM 反復により観測データ (つまり、周辺) 尤度関数は増加しますが、シーケンスが最大尤度推定値に収束するという保証はありません。多峰性分布の場合、これは、開始値に応じて、EM アルゴリズムが観測データ尤度関数の局所的最大値に収束する可能性があることを意味します。局所的最大値を回避するためのさまざまなヒューリスティックまたはメタヒューリスティックなアプローチが存在します。たとえば、ランダム再開ヒルクライミング(いくつかの異なるランダムな初期推定値から開始) やシミュレーテッド アニーリング法の適用などです。 ${\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 とは異なり、このような方法では通常、尤度関数の 1 次導関数および/または 2 次導関数の評価が必要です。

期待最大化は、直接的に改善するのではなく、改善するように機能します。ここでは、前者の改善が後者の改善を意味することが示されています。^[18]${\displaystyle Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})}$${\displaystyle \log p(\mathbf {X} \mid {\boldsymbol {\theta }})}$

確率がゼロでない任意の に対して、次のように書くことができる。 ${\displaystyle \mathbf {Z} }$${\displaystyle p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }})}$

${\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 }}).}$

現在のパラメータ推定値の下での未知のデータの可能な値に対する期待値を求めるには、両辺に を掛けて を合計（または積分）します。左辺は定数の期待値なので、次の式が得られます。 ${\displaystyle \mathbf {Z} }$${\displaystyle \theta ^{(t)}}$${\displaystyle p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }}^{(t)})}$${\displaystyle \mathbf {Z} }$

${\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}}}$

ここで、はそれが置き換える和の負の数で定義されます。この最後の式はを含むのあらゆる値に対して成り立ちます。 ${\displaystyle H({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})}$${\displaystyle {\boldsymbol {\theta }}}$${\displaystyle {\boldsymbol {\theta }}={\boldsymbol {\theta }}^{(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)}),}$

この最後の式を前の式から引くと、

${\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)}).}$

しかし、ギブスの不等式によれば、 ${\displaystyle H({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})\geq H({\boldsymbol {\theta }}^{(t)}\mid {\boldsymbol {\theta }}^{(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 }}}$${\displaystyle Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})}$${\displaystyle \log p(\mathbf {X} \mid {\boldsymbol {\theta }})}$

EMアルゴリズムは、2つの交互の最大化ステップ、つまり座標降下法の例として見ることができます。^{[19] [20]}次の関数を考えてみましょう。

${\displaystyle F(q,\theta ):=\operatorname {E} _{q}[\log L(\theta ;x,Z)]+H(q),}$

ここで、qは観測されていないデータz上の任意の確率分布であり、H(q)は分布qのエントロピーです。この関数は次のように書くことができます。

${\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 情報です。 ${\displaystyle p_{Z\mid X}(\cdot \mid x;\theta )}$${\displaystyle x}$${\displaystyle D_{KL}}$

EM アルゴリズムの手順は次のように考えることができます。

期待ステップ:最大化を選択: ${\displaystyle q}$${\displaystyle F}$

${\displaystyle q^{(t)}=\operatorname {arg\,max} _{q}\ F(q,\theta ^{(t)})}$

最大化ステップ:最大化を選択: ${\displaystyle \theta }$${\displaystyle F}$

${\displaystyle \theta ^{(t+1)}=\operatorname {arg\,max} _{\theta }\ F(q^{(t)},\theta )}$

• EMは混合モデルのパラメータ推定によく使用され、^{[21] [22]}特に量的遺伝学でよく使用されます。^[23]
• 心理測定学において、EM は項目反応理論モデルの項目パラメータと潜在能力を推定するための重要なツールです。
• 欠損データに対処し、未確認の変数を観察する能力により、EM はポートフォリオの価格設定とリスク管理に役立つツールになりつつあります。^[要出典]
• EM アルゴリズム (およびその高速版である順序付きサブセット期待値最大化) は、医療用画像の再構成、特に陽電子放出断層撮影、単一光子放出コンピューター断層撮影、X 線コンピューター断層撮影で広く使用されています。EM の他の高速版については、以下を参照してください。
• 構造工学において、期待値最大化を用いた構造同定（STRIDE）^[24]アルゴリズムは、センサーデータを用いて構造システムの固有振動特性を識別する出力のみの方法です（実用モード解析を参照）。
• EM はデータ クラスタリングにも使用されます。自然言語処理では、このアルゴリズムの代表的な例として、隠れマルコフ モデルのBaum-Welch アルゴリズムと、確率的文脈自由文法の教師なし誘導のinside-outside アルゴリズムが挙げられます。
• 取引間の待ち時間、つまり証券取引所における株式の連続取引間の時間を分析する場合、EMアルゴリズムは非常に有用であることが証明されている。^[25]

カルマン フィルタは通常、オンライン状態推定に使用され、最小分散スムーザはオフラインまたはバッチ状態推定に使用されます。ただし、これらの最小分散ソリューションでは、状態空間モデル パラメータの推定が必要です。EM アルゴリズムは、状態とパラメータの結合推定問題を解決するために使用できます。

フィルタリングとスムージングの EM アルゴリズムは、次の 2 段階の手順を繰り返すことによって生成されます。

Eステップ

現在のパラメータ推定値を使用して設計されたカルマン フィルタまたは最小分散スムージングを操作して、更新された状態推定値を取得します。

Mステップ

最大尤度計算内でフィルタリングまたは平滑化された状態推定値を使用して、更新されたパラメータ推定値を取得します。

カルマンフィルタまたは最小分散スムーザが、加法性ホワイトノイズを持つ単一入力単一出力システムの測定に作用すると仮定する。更新された測定ノイズ分散推定値は、最大尤度計算 から得られる。

${\displaystyle {\widehat {\sigma }}_{v}^{2}={\frac {1}{N}}\sum _{k=1}^{N}{(z_{k}-{\widehat {x}}_{k})}^{2},}$

ここで、はN個のスカラー測定値からフィルタまたはスムーザーによって計算されたスカラー出力推定値です。上記の更新は、ポアソン測定ノイズ強度の更新にも適用できます。同様に、1次自己回帰プロセスの場合、更新されたプロセスノイズ分散推定値は次のように計算できます。 ${\displaystyle {\widehat {x}}_{k}}$${\displaystyle z_{k}}$

${\displaystyle {\widehat {\sigma }}_{w}^{2}={\frac {1}{N}}\sum _{k=1}^{N}{({\widehat {x}}_{k+1}-{\widehat {F}}{\widehat {x}}_{k})}^{2},}$

ここで、およびはフィルタまたはスムーザによって計算されたスカラー状態推定値である。更新されたモデル係数推定値は次のように得られる。 ${\displaystyle {\widehat {x}}_{k}}$${\displaystyle {\widehat {x}}_{k+1}}$

${\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 アルゴリズムで使用される 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 接続と呼ばれる双対アフィン接続による射影として解釈されます。カルバック・ライブラー情報もこれらの用語で理解できます。

を次元 の2つの多変量正規分布の混合からの独立した観測値のサンプルとし、を観測値の起源となる成分を決定する潜在変数とする。^[20]${\displaystyle \mathbf {x} =(\mathbf {x} _{1},\mathbf {x} _{2},\ldots ,\mathbf {x} _{n})}$${\displaystyle n}$${\displaystyle d}$${\displaystyle \mathbf {z} =(z_{1},z_{2},\ldots ,z_{n})}$

${\displaystyle X_{i}\mid (Z_{i}=1)\sim {\mathcal {N}}_{d}({\boldsymbol {\mu }}_{1},\Sigma _{1})}$そして${\displaystyle X_{i}\mid (Z_{i}=2)\sim {\mathcal {N}}_{d}({\boldsymbol {\mu }}_{2},\Sigma _{2}),}$

どこ

${\displaystyle \operatorname {P} (Z_{i}=1)=\tau _{1}\,}$そして${\displaystyle \operatorname {P} (Z_{i}=2)=\tau _{2}=1-\tau _{1}.}$

目的は、ガウス分布間の混合値とそれぞれの平均および共分散を 表す未知のパラメータを推定することです。

${\displaystyle \theta ={\big (}{\boldsymbol {\tau }},{\boldsymbol {\mu }}_{1},{\boldsymbol {\mu }}_{2},\Sigma _{1},\Sigma _{2}{\big )},}$

ここで、不完全データ尤度関数は

${\displaystyle L(\theta ;\mathbf {x} )=\prod _{i=1}^{n}\sum _{j=1}^{2}\tau _{j}\ f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{j},\Sigma _{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)},}$

または

${\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\},}$

ここで、は指標関数であり、は多変量正規分布の確率密度関数です。 ${\displaystyle \mathbb {I} }$${\displaystyle f}$

最後の等式では、各iについて、1 つの指標は0 に等しく、1 つの指標は 1 に等しくなります。したがって、内部の合計は 1 つの項に減ります。 ${\displaystyle \mathbb {I} (z_{i}=j)}$

パラメータθ ^{( t )の現在の推定値を考えると、}ベイズの定理により、 Z _iの条件付き分布は、τで重み付けされた正規密度の比例高さになるように決定されます。

${\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 に対してこの関数を設定することに対応します。

${\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 ステップ セクションの冒頭の式に従って計算されます。 ${\displaystyle \log L(\theta ;\mathbf {x} _{i},Z_{i})}$${\displaystyle P(Z_{i}\mid X_{i}=\mathbf {x} _{i};\theta ^{(t)})}$${\displaystyle \mathbf {x} _{i}}$${\displaystyle T_{j,i}}$

この完全な条件付き期待値は、 τとμ / Σ が別々の線形項で表示され、独立して最大化できる ため、1 つのステップで計算する必要はありません。

${\displaystyle Q(\theta \mid \theta ^{(t)})}$が二次形式であるということは、 を最大化する値を決定することが比較的簡単であることを意味します。また、、およびはすべて別々の線形項として現れるため、すべて独立して最大化できます。 ${\displaystyle \theta }$${\displaystyle \tau }$${\displaystyle ({\boldsymbol {\mu }}_{1},\Sigma _{1})}$${\displaystyle ({\boldsymbol {\mu }}_{2},\Sigma _{2})}$

まず、制約 を持つを考えます。 ${\displaystyle \tau }$${\displaystyle \tau _{1}+\tau _{2}=1}$

${\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}}}$

これは二項分布の最大尤度推定と同じ形をしており、

${\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)}.}$

次の推定値については: ${\displaystyle ({\boldsymbol {\mu }}_{1},\Sigma _{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}}.}$

これは正規分布の加重最大尤度推定と同じ形式なので、

${\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)}}}}$そして${\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)}}}}$

そして対称性により、

${\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)}}}}$そして${\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)}}}.}$

事前に設定されたしきい値を下回る 場合は、反復プロセスを終了します。${\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 などの重要なモデルの多くで導出できます。これらのスペクトル手法では、偽の局所最適値は発生せず、いくつかの規則性条件下で真のパラメータを一貫して推定できます。^[要出典]

• 混合分布
• 複合分布
• 密度推定
• 主成分分析
• 全吸収分光法
• EMアルゴリズムは、Majorize-Minimization（MM）アルゴリズムの特殊なケースとして見ることができます。^[42]

• ホッグ、ロバート; マッキーン、ジョセフ;クレイグ、アレン(2005)。『数理統計学入門』アッパーサドルリバー、ニュージャージー: ピアソン・プレンティス・ホール。pp. 359–364。
• Dellaert, Frank (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 に関するよく書かれた短い本です。
• Bilmes, Jeff (1997)。EM アルゴリズムの簡単なチュートリアルと、ガウス混合モデルおよび隠れマルコフ モデルのパラメータ推定への応用 (技術レポート TR-97-021)。国際コンピュータ サイエンス研究所。ガウス混合分布とガウス混合隠れマルコフモデルの EM 方程式の簡略化された導出が含まれています。
• McLachlan, Geoffrey J.; Krishnan, Thriyambakam (2008)。EMアルゴリズムと拡張(第 2 版)。ホーボーケン: Wiley。ISBN 978-0-471-20170-0。

• EM と混合モデリングのさまざまな 1D、2D、3D デモンストレーションが、SOCRアクティビティとアプレットのペアの一部として提供されます。これらのアプレットとアクティビティは、さまざまな設定でのパラメータ推定のための EM アルゴリズムの特性を経験的に示します。
• ガウス混合分布を含むC++のクラス階層(GPL)
• オンライン教科書「情報理論、推論、学習アルゴリズム」( David JC MacKay著) には、ソフトk平均法アルゴリズムを使用したクラスタリングなどの EM アルゴリズムの簡単な例が含まれており、バージョン 7.2 (第 4 版) の第 33.7 章で説明されているように、EM アルゴリズムの変分観点が強調されています。
• MJ Beal 著の「近似ベイズ推論のための変分アルゴリズム」には、EM と変分ベイズ EM の比較、および変分ベイズ HMM (章) を含むいくつかのモデルの導出が含まれています。
• 期待値最大化アルゴリズム: 短いチュートリアル、Sean Borman による EM アルゴリズムの自己完結型導出。
• Xiaojin Zhu による EM アルゴリズム。
• EM アルゴリズムとその変種: Alexis Roche による非公式チュートリアル。EM と多くの興味深い変種について簡潔かつ非常にわかりやすく説明しています。