1977 - EM 演算法（Expectation-Maximization）#

案件摘要#

1977 年 Dempster、Laird、Rubin 在《Journal of the Royal Statistical Society B》發表 EM 演算法，正式解決「遺失資料下的極大概似估計」這樁懸案。他們證明了 EM 的單調收斂性，讓隱變數模型（如高斯混合）從理論玩具變成實用工具。

前因 -- 為什麼會有這個案子#

統計學家長期依賴極大概似估計（MLE）：給定資料 \(x\) ，求 \(\hat\theta = \arg\max_\theta L(\theta;x)\) 。
但許多模型含有隱變數（latent variables）或遺失資料（missing data）：例如混合模型中「每筆資料來自哪個成分」是未知的；問卷調查中部分欄位缺失。
直接對觀測資料的概似 \(L(\theta;x)\) 求導常導致難解的非線性方程；若「假裝知道」隱變數 \(z\) ，完整資料概似 \(L(\theta;x,z)\) 卻往往好解得多。
早年已有零散的特例解法（如 1894 年 Pearson 的混合分布手算、遺傳學中的 gene counting），但缺乏統一框架與理論保證——這正是本案要偵破的核心。

線索與推理 -- 數學式、程式、理論#

核心概念：完整資料 vs 觀測資料#

設觀測資料 \(X\) 、隱變數 \(Z\) 、參數 \(\theta\) 。觀測概似可分解為：

\[ L(\theta; X) = p(X|\theta) = \sum_Z p(X, Z|\theta) \]

求和（或積分）使極大化困難。定義完整資料概似 \(p(X, Z|\theta)\) ，它是 \(Z\) 的函數。

E 步（Expectation）#

在目前參數 \(\theta^{(t)}\) 下，計算完整資料對數概似的期望：

\[ Q(\theta \mid \theta^{(t)}) = \mathbb{E}_{Z \sim p(Z \mid X, \theta^{(t})} \left[ \log p(X, Z \mid \theta) \right] \]

直覺：E 步是「補齊線索」——用目前的推論填入缺失的 \(Z\) 。

M 步（Maximization）#

極大化這個期望：

\[ \theta^{(t+1)} = \arg\max_\theta Q(\theta \mid \theta^{(t)}) \]

直覺：M 步是「根據補好的線索更新嫌犯輪廓」。兩步交替，直到收斂。

單調收斂性（本案最關鍵的證據）#

利用 Jensen 不等式可證：每次迭代概似不會變差：

\[ L(\theta^{(t+1)}; X) \ge L(\theta^{(t)}; X) \]

證明骨架：對任意 \(\theta\) ，

\[ \log p(X|\theta) = Q(\theta|\theta^{(t)}) + H(\theta|\theta^{(t)}), \quad H(\theta|\theta^{(t)}) = -\mathbb{E}[\log p(Z|X,\theta)] \]

由 Gibbs 不等式 \(H(\theta|\theta^{(t)}) \ge H(\theta^{(t)}|\theta^{(t)})\) ，故 M 步提升 \(Q\) 就必提升 \(\log p(X|\theta)\) 。這保證 EM 收斂到概似函數的穩定點（可能是局部極大或鞍點）。

經典應用：高斯混合模型（GMM）#

設 \(K\) 個高斯成分，參數 \(\theta = \{(\pi_k, \mu_k, \Sigma_k)\}_{k=1}^K\) ，隱變數 \(z_i \in \{1,\dots,K\}\) 。

E 步：計算責任（responsibility）

\[ \gamma_{ik} = P(z_i = k \mid x_i, \theta^{(t)}) = \frac{\pi_k^{(t)} \, \mathcal{N}(x_i \mid \mu_k^{(t)}, \Sigma_k^{(t)})} {\sum_{j=1}^K \pi_j^{(t)} \, \mathcal{N}(x_i \mid \mu_j^{(t)}, \Sigma_j^{(t)})} \]

M 步：更新參數

\[ \mu_k^{(t+1)} = \frac{\sum_i \gamma_{ik} x_i}{\sum_i \gamma_{ik}}, \quad \Sigma_k^{(t+1)} = \frac{\sum_i \gamma_{ik} (x_i - \mu_k^{(t+1)})(x_i - \mu_k^{(t+1)})^\top}{\sum_i \gamma_{ik}}, \quad \pi_k^{(t+1)} = \frac{\sum_i \gamma_{ik}}{N} \]

Python 實作：二維 GMM 聚類（簡化版 EM）#

import numpy as np

np.random.seed(42)

# 產生二維混合資料：兩群高斯
X = np.vstack([
    np.random.multivariate_normal([0, 0], [[1, .3], [.3, 1]], 200),
    np.random.multivariate_normal([4, 4], [[1, -.2], [.3, 1]], 200),
])
N, d = X.shape
K = 2

# 初始化參數
pi  = np.full(K, 1 / K)
mu  = X[np.random.choice(N, K, replace=False)]
Sigma = np.array([np.eye(d)] * K)

def log_lik(X, pi, mu, Sigma):
    ll = 0.0
    for i in range(N):
        p = sum(pi[k] * _mvn(X[i], mu[k], Sigma[k]) for k in range(K))
        ll += np.log(p + 1e-300)
    return ll

def _mvn(x, m, S):
    diff = x - m
    return np.exp(-0.5 * diff @ np.linalg.solve(S, diff)) / \
           np.sqrt((2 * np.pi) ** d * np.linalg.det(S))

for it in range(50):
    # E 步：計算責任 gamma[i, k]
    gamma = np.zeros((N, K))
    for k in range(K):
        gamma[:, k] = pi[k] * np.array([_mvn(x, mu[k], Sigma[k]) for x in X])
    gamma /= gamma.sum(axis=1, keepdims=True)

    # M 步：更新 pi, mu, Sigma
    Nk = gamma.sum(axis=0)
    pi = Nk / N
    mu = (gamma.T @ X) / Nk[:, None]
    for k in range(K):
        diff = X - mu[k]
        Sigma[k] = (gamma[:, k, None] * diff).T @ diff / Nk[k]

    if it % 10 == 0:
        print(f"iter {it:2d}  log-likelihood = {log_lik(X, pi, mu, Sigma):.2f}")

print("混洗後的均值：\n", mu)
print("混合比例：", pi)

執行後可看到 log-likelihood 單調上升，且 \(\mu\) 收斂到接近 \([0,0]\) 與 \([4,4]\) 。

結案 -- 後果與影響#

EM 成為統計與機器學習的基礎演算法：GMM 聚類、隱馬可夫模型（HMM）的 Baum-Welch 訓練本質上就是 EM、缺失資料填補、因子分析。
理論後續發展：收斂速率分析、變分 EM、與貝葉斯推論結合（MAP-EM）。
局限也被認清：只保證局部極大、對初始值敏感、可能收斂到鞍點（Wu 1983 補充了收斂條件）。
Dempster、Laird、Rubin 的論文成為史上被引用次數最高的統計論文之一，ISI 評選為百大最具影響力的科學論文。

關鍵人物與文獻#

A. P. Dempster（哈佛）：統計推論大師，亦以 Dempster-Shafer 證據理論聞名。
D. B. Laird、D. B. Rubin：Rubin 亦提出「多重插補（multiple imputation）」。
Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. JRSS-B 39(1), 1–38.
Wu, C. F. J. (1983). On the convergence properties of the EM algorithm. Annals of Statistics.
Bishop, C. M. (2006). Pattern Recognition and Machine Learning, Ch. 9.