EM 算法

EM（Expectation-Maximization）算法用于含有隐变量的概率模型参数估计。它常用于高斯混合模型、隐马尔可夫模型等场景。

概念详解

假设观测变量为 $X$，隐变量为 $Z$，模型参数为 $\theta$。我们希望最大化观测数据的对数似然：

$$\log p(X|\theta )=\log \sum _{Z}p(X,Z|\theta )$$

由于求和在对数内部，直接优化通常困难。EM 算法通过两步迭代解决：

• E 步：用当前参数估计隐变量后验分布。
• M 步：在该后验分布下最大化完全数据似然的期望。

数学公式及其推导

引入任意分布 $q(Z)$：

$$\log p(X|\theta )=\log \sum _{Z}q(Z)\frac{p(X,Z|\theta )}{q(Z)}$$

由 Jensen 不等式：

$$\log p(X|\theta )\ge \sum _{Z}q(Z)\log \frac{p(X,Z|\theta )}{q(Z)}$$

右侧称为证据下界：

$$\mathcal{L}(q,\theta )=\sum _{Z}q(Z)\log p(X,Z|\theta )-\sum _{Z}q(Z)\log q(Z)$$

E 步令 $q(Z)=p(Z|X,{\theta }^{old})$，M 步最大化：

$$Q(\theta ,{\theta }^{old})=E_{Z|X,{\theta }^{old}}[\log p(X,Z|\theta )]$$$${\theta }^{new}=\arg \underset{\theta }{max}Q(\theta ,{\theta }^{old})$$

应用代码：一维高斯混合模型

import numpy as np

np.random.seed(0)
X = np.r_[np.random.normal(-2, 0.6, 100), np.random.normal(2, 0.8, 100)]
K = 2
pi = np.ones(K) / K
mu = np.array([-1.0, 1.0])
sigma = np.ones(K)

def normal_pdf(x, mu, sigma):
    return np.exp(-0.5 * ((x - mu) / sigma) ** 2) / (np.sqrt(2 * np.pi) * sigma)

for _ in range(50):
    gamma = np.zeros((len(X), K))
    for k in range(K):
        gamma[:, k] = pi[k] * normal_pdf(X, mu[k], sigma[k])
    gamma /= gamma.sum(axis=1, keepdims=True)

    Nk = gamma.sum(axis=0)
    pi = Nk / len(X)
    mu = (gamma * X[:, None]).sum(axis=0) / Nk
    sigma = np.sqrt((gamma * (X[:, None] - mu) ** 2).sum(axis=0) / Nk)

print("pi:", pi)
print("mu:", mu)
print("sigma:", sigma)

小结

EM 的本质是“先估计隐变量分布，再更新模型参数”。它保证观测似然不下降，但可能收敛到局部最优。