EMアルゴリズム

EMアルゴリズムをJensenの不等式によるELBO導出から解説。EステップとMステップの理論的根拠、古典的な2枚のコインの例をPythonで実行検証し、初期値依存性という落とし穴まで扱います。

9 min read

EMアルゴリズム（Expectation-Maximization Algorithm）は、観測データに隠れた変数（潜在変数）が含まれる統計モデルのパラメータを推定するための反復的なアルゴリズムです。

例えば、フィットネスクラブの利用者の年齢分布を考えた場合、筋トレ目的の20代とメタボ対策目的の50代という2つのグループが存在し、それぞれの年齢分布が正規分布に従うと仮定できます。この場合、各利用者がどちらのグループに属するかは直接観測できない「潜在変数」となります。

EMアルゴリズムの基本的な考え方は、以下の2つのステップを交互に繰り返すことで、モデルのパラメータと潜在変数の推定を同時に行うというものです。

  1. Eステップ (Expectation Step): 現在のモデルパラメータを使って、観測データから潜在変数の期待値（または確率分布）を推定します。上記の例では、各利用者がどちらのグループに属するかの「負担率」を計算します。
  2. Mステップ (Maximization Step): Eステップで推定された潜在変数の期待値を用いて、モデルのパラメータを最大化します。上記の例では、負担率を重みとして、各グループの正規分布の平均や分散といったパラメータを更新します。

このプロセスを繰り返すことで、モデルのパラメータと潜在変数の推定が徐々に改善され、最終的に対数尤度が局所最大値に収束します。本記事では、この「なぜ繰り返すだけで対数尤度が改善するのか」を Jensen の不等式から導出し、GMMとは異なる古典的な例（2枚のコイン問題）で実際にEMアルゴリズムを実行検証します。

Jensenの不等式によるELBOの導出

観測データを \(x\) 、潜在変数を \(z\) 、パラメータを \(\theta\) とします。対数尤度 \(\log p(x|\theta)\) を直接最大化したいのですが、潜在変数を積分消去した周辺尤度 \(p(x|\theta) = \int p(x,z|\theta)dz\) の対数を取ると、対数の中に積分が入ってしまい解析的に扱いにくくなります。

そこで、任意の（\(z\) 上の）確率分布 \(q(z)\) を導入し、恒等的に次のように変形します。

\[ \log p(x|\theta) = \log \int p(x,z|\theta)dz = \log \int q(z)\frac{p(x,z|\theta)}{q(z)}dz = \log \mathbb{E}_{q(z)}\left[\frac{p(x,z|\theta)}{q(z)}\right] \]

対数関数 \(\log(\cdot)\) は凹関数なので、Jensenの不等式 \(\log \mathbb{E}[Y] \ge \mathbb{E}[\log Y]\) が任意の確率変数 \(Y>0\) に対して成り立ちます。これを上式に適用すると、

\[ \log p(x|\theta) = \log \mathbb{E}_{q(z)}\left[\frac{p(x,z|\theta)}{q(z)}\right] \ge \mathbb{E}_{q(z)}\left[\log\frac{p(x,z|\theta)}{q(z)}\right] =: \mathcal{L}(q,\theta) \]

この下界 \(\mathcal{L}(q,\theta)\) を**変分下界（Evidence Lower Bound, ELBO）**と呼びます。対数尤度は常にELBO以上であり、両者の差は次のように厳密に評価できます。

\[ \log p(x|\theta) - \mathcal{L}(q,\theta) = \log p(x|\theta) - \mathbb{E}_{q(z)}\left[\log\frac{p(x,z|\theta)}{q(z)}\right] = \mathbb{E}_{q(z)}\left[\log\frac{q(z)}{p(z|x,\theta)}\right] = KL(q(z) \| p(z|x,\theta)) \]

（\(p(x,z|\theta) = p(z|x,\theta)p(x|\theta)\) を代入し、\(\log p(x|\theta)\) が \(z\) に依存しないため期待値の外に出せることを使っています。）KLダイバージェンスは常に非負なので、この式は「対数尤度 = ELBO + KLダイバージェンス」という厳密な等式であり、ELBOが対数尤度の下界になっていることを裏付けます。さらに、等号が成立するのは \(KL(q(z)\|p(z|x,\theta))=0\) 、すなわち \(q(z) = p(z|x,\theta)\) のときに限られることも分かります。

EMアルゴリズムはELBOの座標降下法

この事実から、EMアルゴリズムの2つのステップは「\(q\) と \(\theta\) を交互に最適化してELBOを最大化する座標降下法」として統一的に理解できます。

  • Eステップ: \(\theta\) を現在の推定値 \(\hat{\theta}\) に固定し、ELBOを最大化する \(q(z)\) を求めます。上記の議論から、これは \(q(z) = p(z|x,\hat{\theta})\) とすること（=事後分布そのものを計算すること）に他なりません。このとき \(KL=0\) となるため、\(\mathcal{L}(q,\hat{\theta}) = \log p(x|\hat{\theta})\) 、つまりELBOと対数尤度が一致します。
  • Mステップ: \(q(z)=p(z|x,\hat{\theta})\) を固定し、ELBOを最大化する \(\theta\) を求めます。
\[ \mathcal{L}(q,\theta) = \mathbb{E}_{q(z)}[\log p(x,z|\theta)] - \mathbb{E}_{q(z)}[\log q(z)] = Q(\theta,\hat{\theta}) + H(q) \]

エントロピー項 \(H(q)\) は \(\theta\) に依存しないため、ELBOを \(\theta\) について最大化することは、Q関数 \(Q(\theta,\hat{\theta}) = \mathbb{E}_{p(z|x,\hat{\theta})}[\log p(x,z|\theta)]\) を最大化することと同値です。これがMステップです。

Mステップで \(\theta\) を更新すると（Q関数が改善する限り）ELBOも改善しますが、更新後の \(q(z)=p(z|x,\hat{\theta})\) はもはや新しいパラメータ \(\theta^{(t+1)}\) の下での事後分布とは一致しないため、\(KL>0\) となり、ELBOは対数尤度を下回った状態になります。しかし「対数尤度 \(\ge\) ELBO」は常に成り立つので、

\[ \log p(x|\theta^{(t+1)}) \ge \mathcal{L}(q,\theta^{(t+1)}) \ge \mathcal{L}(q,\theta^{(t)}) = \log p(x|\theta^{(t)}) \]

という不等式の連鎖が得られ、対数尤度がEMの各反復で単調に非減少であることが証明されます（1つ目の不等号は「対数尤度≥ELBO」、2つ目はMステップがELBOを改善すること、最後の等号はEステップ直後の性質です）。この一般的な証明はGMMに限らず任意の潜在変数モデルに対して成立します。GMM特有の共分散の特異性や収束の数値的な詳細は k-means法とGMMの理論・比較・Python実装 で扱っているので、そちらも参照してください。

混合ガウスモデルとEMアルゴリズム

EMアルゴリズムは、複数のガウス分布が重なり合って観測データが生成されると仮定する混合ガウスモデル (Gaussian Mixture Model, GMM) のパラメータ推定によく用いられます。ここでは、上記の一般論をGMMに当てはめた場合の更新式を確認します（完全データの対数尤度、その期待値がQ関数になることは前節の通りです）。

Q関数

EMアルゴリズムでは、観測データ \(x\) と潜在変数 \(z\) の両方がわかっている仮想的な状況を「完全データ」と呼びます。完全データにおける同時分布 \(p(x, z | \theta)\) の対数尤度を考えます。

\[ Q(\theta, \hat{\theta}) = \mathbb{E}_{p(z|x,\hat{\theta})}[\log p(x,z|\theta)] = \int p(z|x,\hat{\theta})\log p(x,z|\theta)dz \]

EMアルゴリズムによる混合ガウスモデルのパラメータ更新式

混合ガウスモデルのパラメータ（各ガウス分布の重み \(\pi_j\) 、平均 \(\mu_j\) 、分散 \(\sigma_j^2\) ）をEMアルゴリズムで推定する際の更新式は以下の通りです。

EMアルゴリズムによる混合ガウスモデルのパラメータ更新式

  1. 初期化: 各ガウス分布のパラメータ \(\hat{\pi}_j^{(0)}, \hat{\mu}_j^{(0)}, \hat{\sigma}_j^{2(0)}\) をランダムな値で初期化します。
  2. Eステップ: 各データ点 \(x_i\) が、どのガウス分布から生成されたかを示す「負担率」 \(r_{ij}\) を計算します。
\[ r_{ij} = p(z_{ij}=1 \mid x_i, \hat{\theta}^{(t)}) = \frac{\hat{\pi}_j^{(t)} \mathcal{N}(x_i \mid \hat{\mu}_j^{(t)}, \hat{\sigma}_j^{2(t)})}{\sum_{k=1}^K \hat{\pi}_k^{(t)} \mathcal{N}(x_i \mid \hat{\mu}_k^{(t)}, \hat{\sigma}_k^{2(t)})} \]

ここで \(z_{ij}=1\) はデータ点 \(x_i\) が \(j\) 番目のガウス分布に属することを示します。

  1. Mステップ: Eステップで計算された負担率 \(r_{ij}\) を用いて、新しいパラメータ \(\hat{\theta}^{(t+1)}\) を計算します。
\[ N_j = \sum_{i=1}^N r_{ij}, \qquad \hat{\pi}_j^{(t+1)} = \frac{N_j}{N} \] \[ \hat{\mu}_j^{(t+1)} = \frac{1}{N_j} \sum_{i=1}^N r_{ij} x_i, \qquad \hat{\sigma}_j^{2(t+1)} = \frac{1}{N_j} \sum_{i=1}^N r_{ij} (x_i - \hat{\mu}_j^{(t+1)})^2 \]
  1. 収束判定: パラメータの変化が十分に小さくなるか、最大反復回数に達するまでステップ2と3を繰り返します。

古典例：2枚のコイン問題（Do & Batzoglou, 2008）

GMMは連続値データの例ですが、EMアルゴリズムはもっと単純な離散的設定でも同じ理屈で動きます。ここでは教育目的でよく使われる「2枚のコイン問題」を実装し、実際に収束の様子を数値で確認します。

設定: コインAとコインBがあり、それぞれ未知の確率 \(\theta_A, \theta_B\) で表が出ます。5回の試行があり、各試行では10回コインを投げますが、どちらのコインを使ったか（潜在変数）は記録されていません。観測データ（各試行の表の枚数）だけから \(\theta_A, \theta_B\) を推定します。

import numpy as np

# Do & Batzoglou (2008) の5試行×10投げのデータ
data = np.array([
    [1,0,0,0,1,1,0,1,0,1],
    [1,1,1,1,0,1,1,1,1,1],
    [1,0,1,1,1,1,1,0,1,1],
    [1,0,1,0,0,0,1,1,0,0],
    [0,1,1,1,0,1,1,1,0,1],
])
n_heads = data.sum(axis=1)
n_tails = data.shape[1] - n_heads

def bernoulli_loglik(k, n, p):
    return k*np.log(p) + (n-k)*np.log(1-p)

def em_two_coins(n_heads, n_tails, theta_a0, theta_b0, n_iter=10):
    theta_a, theta_b = theta_a0, theta_b0
    n = n_heads + n_tails
    loglik_history = []
    for t in range(n_iter):
        # Eステップ: 各試行がコインAである負担率
        loglik_a = bernoulli_loglik(n_heads, n, theta_a)
        loglik_b = bernoulli_loglik(n_heads, n, theta_b)
        m = np.maximum(loglik_a, loglik_b)
        wa, wb = np.exp(loglik_a - m), np.exp(loglik_b - m)
        r_a = wa / (wa + wb)
        r_b = 1 - r_a

        # 観測データの対数尤度（コインの選択は五分五分の混合と仮定）
        mix = 0.5*np.exp(loglik_a) + 0.5*np.exp(loglik_b)
        loglik_history.append(np.log(mix).sum())

        # Mステップ: 負担率で重み付けした最尤推定
        theta_a = (r_a*n_heads).sum() / (r_a*n).sum()
        theta_b = (r_b*n_heads).sum() / (r_b*n).sum()
    return theta_a, theta_b, loglik_history

theta_a, theta_b, loglik_history = em_two_coins(n_heads, n_tails, 0.6, 0.5, n_iter=10)
for i, ll in enumerate(loglik_history):
    print(f"iter {i+1:2d}: log-likelihood = {ll:.6f}")
print(f"converged theta_A = {theta_a:.4f}, theta_B = {theta_b:.4f}")

実行結果は以下の通りです。

iter  1: log-likelihood = -33.093863
iter  2: log-likelihood = -31.859258
iter  3: log-likelihood = -31.723116
iter  4: log-likelihood = -31.627827
iter  5: log-likelihood = -31.584947
iter  6: log-likelihood = -31.573058
iter  7: log-likelihood = -31.570678
iter  8: log-likelihood = -31.570275
iter  9: log-likelihood = -31.570212
iter 10: log-likelihood = -31.570202
converged theta_A = 0.7967, theta_B = 0.5197

対数尤度が \(-33.09 \to -31.57\) と単調に増加し、10回の反復で \(\theta_A \approx 0.797\) , \(\theta_B \approx 0.520\) に収束しました。これは元論文（Do & Batzoglou, 2008）が報告する \(\theta_A\approx0.80\) , \(\theta_B\approx0.52\) とよく一致します。前節で証明した「対数尤度の単調非減少性」が、GMMとは異なる離散モデルでも数値的に成立することが確認できます。

対数尤度と \(\theta_A, \theta_B\) の収束を可視化すると、以下のようになります。

2枚のコイン問題におけるEMアルゴリズムの収束。左：対数尤度の単調増加。右：θA・θBの推定値の収束

エッジケース・注意点

  • 初期値への依存性（局所最適解）: EMアルゴリズムが保証するのは対数尤度の単調非減少性のみで、大域最適解への収束は保証されません。同じデータに対して初期値を変えて実行すると、以下のように挙動が変わります。
for (ia, ib) in [(0.6, 0.5), (0.5, 0.6), (0.9, 0.1), (0.5, 0.5)]:
    ta, tb, llh = em_two_coins(n_heads, n_tails, ia, ib, n_iter=15)
    print(f"init=({ia},{ib}) -> converged=({ta:.4f},{tb:.4f}), final loglik={llh[-1]:.4f}")

init=(0.6,0.5) -> converged=(0.7968,0.5196), final loglik=-31.5702
init=(0.5,0.6) -> converged=(0.5196,0.7968), final loglik=-31.5702
init=(0.9,0.1) -> converged=(0.7968,0.5196), final loglik=-31.5702
init=(0.5,0.5) -> converged=(0.6600,0.6600), final loglik=-32.0518

初期値 \((0.6,0.5)\) と \((0.5,0.6)\) は「コインAとBのラベルを入れ替えただけ」の対称な解に収束します（ラベルスイッチング：潜在変数モデルは一般にラベルの並べ替えに関して非識別的です）。一方、完全に対称な初期値 \((0.5,0.5)\) から始めると、\(\theta_A=\theta_B=0.66\) という劣った鞍点的な解（対数尤度 \(-32.05\) ）に留まってしまい、真の最適解（対数尤度 \(-31.57\) ）に到達できません。実務上は複数の初期値からEMを実行し、最も対数尤度が高い解を採用する（マルチスタート）のが定石です。

  • 収束判定の基準: パラメータの変化量、Q関数の変化量、対数尤度の変化量のいずれで収束を判定するかによって、停止タイミングがわずかに変わります。上記の例では10反復で対数尤度の変化が \(10^{-4}\) 以下になっています。
  • GMM特有の落とし穴: 共分散が特異行列に近づく発散（対数尤度が原理上\(+\infty\) に発散しうる問題）は、この記事の一般的な議論だけでは説明できないGMM特有の病理です。詳しくは k-means法とGMMの理論・比較・Python実装 を参照してください。

発展：EM収束理論の最近の進展

EMアルゴリズムの収束速度の解析は、対数尤度が単調非減少であることの証明（本記事で導出した内容）よりも難しい問題として、現在も研究が続いています。例えば Caprio & Johansen (2025, Biometrika) は、対数ソボレフ不等式（logarithmic Sobolev inequality）という条件の下でEMアルゴリズムの収束速度を非漸近的に評価する枠組みを提案しています。これはEMをKLダイバージェンスに関する自由エネルギー汎函数の座標ごとの最小化とみなす視点に基づいており、本記事のJensenの不等式による導出とも整合する現代的な理論の発展です。

参考

この記事に関連する無料ツール

DevToolBox Moving Average Calculator SMA/EMA/WMA をブラウザで即座に計算 DevToolBox JSON Formatter JSONを整形・圧縮・バリデーション DevToolBox Regex Tester 正規表現をリアルタイムテスト

関連記事