第9章：EM算法及其推广（EM Algorithm and Its Extension）

一、导读

1.1 本章学习目标

第9章介绍EM（Expectation-Maximization）算法——处理含有隐变量模型参数估计的核心方法。当模型中含有无法直接观测的变量时，极大似然估计变得困难，EM算法通过交替执行E步（求期望）和M步（最大化）来迭代求解。

本章的核心学习目标包括：

理解含有隐变量的模型：什么是隐变量？为什么它使参数估计变得困难？
掌握EM算法的原理：E步和M步分别在做什么？
理解EM算法的收敛性：为什么EM算法能保证对数似然单调递增？
掌握混合高斯模型（GMM）：EM在GMM中的具体应用
了解EM算法的推广：GEM算法
实现GMM的EM算法

1.2 为什么本章如此重要

EM算法是统计学习中最重要的算法之一。许多经典模型——混合高斯模型、隐马尔可夫模型、变分推断——的参数学习都依赖EM算法。

EM算法揭示了一个深刻的思想：当直接优化困难时，通过引入隐变量将问题分解为更简单的子问题。这种"分而治之"的思想在机器学习中无处不在。

1.3 核心问题引导

在阅读本章时，请带着以下问题思考：

什么是隐变量？为什么含有隐变量的模型不能直接用极大似然估计？

EM算法的E步和M步分别在做什么？直觉上如何理解？

为什么EM算法能保证收敛？每次迭代对数似然为什么会增加？

GMM中EM算法的具体推导是什么？每个参数如何更新？

EM算法的初始值敏感吗？如何避免陷入局部最优？

EM算法与K-means有什么关系？

二、核心概念详解

9.1 含有隐变量的模型

9.1.1 什么是隐变量

隐变量（Latent Variable）是无法直接观测到的变量。例如：

混合模型中，每个样本属于哪个分量是隐变量
隐马尔可夫模型中，状态序列是隐变量
主题模型中，文档的主题分配是隐变量

9.1.2 为什么隐变量使问题变难

含有隐变量时，观测数据的对数似然为：

L(θ) = log P(X|θ) = log Σ_Z P(X, Z|θ)

由于对数里面有求和，直接对 θ 求导并令其为零变得非常困难。

python

import numpy as np

def hidden_variable_difficulty():
    """演示隐变量使极大似然估计变难"""
    # 假设数据来自两个高斯分布的混合
    # 但我们不知道每个样本来自哪个分布（隐变量）
    
    np.random.seed(42)
    # 生成混合数据
    n1, n2 = 60, 40
    X1 = np.random.randn(n1, 2) + np.array([2, 2])
    X2 = np.random.randn(n2, 2) + np.array([-2, -2])
    X = np.vstack([X1, X2])
    
    # 如果我们知道每个样本属于哪个分量（没有隐变量）
    # 可以直接用MLE估计每个分量的参数
    mu1_known = X1.mean(axis=0)
    mu2_known = X2.mean(axis=0)
    print("如果知道隐变量（每个样本的归属）:")
    print(f"  μ1 = {mu1_known}")
    print(f"  μ2 = {mu2_known}")
    
    # 如果不知道隐变量，需要EM算法
    print("\n如果不知道隐变量，需要EM算法迭代求解")
    print("因为 L(θ) = log[π1*N(x|μ1,Σ1) + π2*N(x|μ2,Σ2)]")
    print("对数里面有求和，无法直接求导求解")

hidden_variable_difficulty()

9.2 EM算法

9.2.1 算法推导

EM算法的核心思想是构造对数似然的下界，然后交替优化下界。

完全数据对数似然：log P(X, Z|θ)（如果Z已知，问题简单）

Q函数：Q(θ, θ^t) = E_Z[log P(X, Z|θ) | X, θ^t]

EM算法的每一步：

E步：计算 Q(θ, θ^t)，即在当前参数 θ^t 下，求完全数据对数似然关于隐变量 Z 的期望
M步：最大化 Q(θ, θ^t)，即 θ^{t+1} = argmax_θ Q(θ, θ^t)

python

import numpy as np

def em_intuition():
    """EM算法的直觉理解"""
    print("EM算法的直觉理解:")
    print()
    print("E步（Expectation）:")
    print("  固定参数θ，推断隐变量Z的后验分布 P(Z|X,θ)")
    print("  即：给定当前模型，每个样本属于哪个分量的概率是多少？")
    print()
    print("M步（Maximization）:")
    print("  固定隐变量的后验分布，更新参数θ使完全数据似然最大")
    print("  即：根据样本的归属概率，重新估计每个分量的参数")
    print()
    print("交替执行E步和M步，直到收敛")
    print("每次迭代都保证对数似然不减少")

em_intuition()

9.2.2 收敛性证明

EM算法保证每次迭代后对数似然不减少：

L(θ^{t+1}) ≥ L(θ^t)

证明基于Jensen不等式。关键步骤：

利用Jensen不等式构造对数似然的下界

M步最大化这个下界

因此对数似然也增大

python

import numpy as np

def em_convergence_demo():
    """演示EM算法的收敛性"""
    np.random.seed(42)
    
    # 生成混合高斯数据
    n = 200
    X1 = np.random.randn(n // 2, 2) + np.array([3, 3])
    X2 = np.random.randn(n // 2, 2) + np.array([-2, -2])
    X = np.vstack([X1, X2])
    
    # 简化的EM（只估计均值）
    mu1 = np.array([2.0, 2.0])
    mu2 = np.array([-1.0, -1.0])
    sigma = np.eye(2)
    
    log_likelihoods = []
    
    for iteration in range(20):
        # E步：计算后验概率
        p1 = np.exp(-0.5 * np.sum((X - mu1) ** 2, axis=1))
        p2 = np.exp(-0.5 * np.sum((X - mu2) ** 2, axis=1))
        total = p1 + p2
        gamma1 = p1 / total
        gamma2 = p2 / total
        
        # 计算对数似然
        ll = np.sum(np.log(total + 1e-300))
        log_likelihoods.append(ll)
        
        # M步：更新参数
        mu1 = np.sum(gamma1[:, np.newaxis] * X, axis=0) / np.sum(gamma1)
        mu2 = np.sum(gamma2[:, np.newaxis] * X, axis=0) / np.sum(gamma2)
    
    print("EM收敛过程:")
    print(f"{'迭代':>5} {'对数似然':>12} {'增量':>10}")
    print("-" * 30)
    for i, ll in enumerate(log_likelihoods):
        delta = ll - log_likelihoods[i-1] if i > 0 else 0
        print(f"{i:>5} {ll:>12.4f} {delta:>10.4f}")
    
    print(f"\n最终均值估计: μ1={mu1}, μ2={mu2}")
    print("对数似然单调递增，验证了EM的收敛性")

em_convergence_demo()

9.3 混合高斯模型（GMM）

9.3.1 模型定义

GMM假设数据由K个高斯分布混合生成：

P(x) = Σ_k π_k * N(x | μ_k, Σ_k)

其中 π_k 是混合权重（Σπ_k = 1），μ_k 和 Σ_k 是第k个高斯分量的均值和协方差。

隐变量 Z：每个样本属于哪个分量，P(Z=k) = π_k。

9.3.2 GMM的EM算法

E步：计算后验概率（responsibility）

γ_{nk} = π_k * N(x_n | μ_k, Σ_k) / Σ_j π_j * N(x_n | μ_j, Σ_j)

M步：更新参数

N_k = Σ_n γ_{nk}
μ_k = (1/N_k) Σ_n γ_{nk} * x_n
Σ_k = (1/N_k) Σ_n γ_{nk} * (x_n - μ_k)(x_n - μ_k)^T
π_k = N_k / N

python

import numpy as np
from scipy.stats import multivariate_normal

class GMM:
    def __init__(self, n_components=3, n_iter=100, tol=1e-6):
        self.K = n_components
        self.n_iter = n_iter
        self.tol = tol
        self.log_likelihoods = []
    
    def fit(self, X):
        N, D = X.shape
        
        # 初始化
        np.random.seed(42)
        indices = np.random.choice(N, self.K, replace=False)
        self.mu = X[indices].copy()
        self.sigma = np.array([np.eye(D) for _ in range(self.K)])
        self.pi = np.ones(self.K) / self.K
        
        for iteration in range(self.n_iter):
            # E步
            gamma = np.zeros((N, self.K))
            for k in range(self.K):
                gamma[:, k] = self.pi[k] * multivariate_normal.pdf(
                    X, self.mu[k], self.sigma[k] + 1e-6 * np.eye(D))
            
            total = gamma.sum(axis=1, keepdims=True)
            gamma /= total + 1e-300
            
            # 对数似然
            ll = np.sum(np.log(total.flatten() + 1e-300))
            self.log_likelihoods.append(ll)
            
            # 检查收敛
            if len(self.log_likelihoods) > 1:
                if abs(self.log_likelihoods[-1] - self.log_likelihoods[-2]) < self.tol:
                    print(f"收敛于第 {iteration + 1} 次迭代")
                    break
            
            # M步
            Nk = gamma.sum(axis=0)
            for k in range(self.K):
                self.mu[k] = (gamma[:, k:k+1].T @ X) / Nk[k]
                diff = X - self.mu[k]
                self.sigma[k] = (diff.T * gamma[:, k]) @ diff / Nk[k]
                self.sigma[k] += 1e-6 * np.eye(D)  # 防止奇异
            self.pi = Nk / N
        
        return self
    
    def predict(self, X):
        gamma = np.zeros((len(X), self.K))
        for k in range(self.K):
            gamma[:, k] = self.pi[k] * multivariate_normal.pdf(
                X, self.mu[k], self.sigma[k])
        return np.argmax(gamma, axis=1)
    
    def predict_proba(self, X):
        gamma = np.zeros((len(X), self.K))
        for k in range(self.K):
            gamma[:, k] = self.pi[k] * multivariate_normal.pdf(
                X, self.mu[k], self.sigma[k])
        gamma /= gamma.sum(axis=1, keepdims=True)
        return gamma

# 测试
np.random.seed(42)
X1 = np.random.randn(100, 2) + np.array([3, 3])
X2 = np.random.randn(100, 2) + np.array([-2, -2])
X3 = np.random.randn(100, 2) + np.array([3, -2])
X = np.vstack([X1, X2, X3])

gmm = GMM(n_components=3, n_iter=100)
gmm.fit(X)

print(f"\n估计的均值:")
for k in range(3):
    print(f"  分量{k}: {gmm.mu[k]}")
print(f"估计的混合权重: {gmm.pi}")

9.4 EM与K-means的关系

EM算法在GMM中有一个有趣的特殊情况：当所有高斯分量的协方差趋近于0时，GMM的EM退化为K-means算法。

python

import numpy as np

def em_vs_kmeans():
    """对比GMM（EM）和K-means"""
    np.random.seed(42)
    
    # 生成椭圆形簇的数据
    X1 = np.random.randn(100, 2) @ np.array([[2, 0], [0, 0.5]]) + np.array([3, 3])
    X2 = np.random.randn(100, 2) @ np.array([[0.5, 0], [0, 2]]) + np.array([-2, -2])
    X = np.vstack([X1, X2])
    
    # GMM可以捕捉椭球形簇
    gmm = GMM(n_components=2, n_iter=50)
    gmm.fit(X)
    
    print("GMM（EM）:")
    print(f"  估计均值: {gmm.mu}")
    print(f"  协方差矩阵的对角线:")
    for k in range(2):
        print(f"    分量{k}: diag = {np.diag(gmm.sigma[k])}")
    print(f"  → GMM可以捕捉不同形状的簇")
    
    print("\nK-means:")
    print("  → 只能捕捉球状簇（等方差假设）")
    print("  → GMM是K-means的概率化推广")

em_vs_kmeans()

9.5 EM算法的推广：GEM

GEM（Generalized EM）算法放宽了M步的要求：不要求完全最大化Q函数，只需要使Q函数增大即可。

Q(θ^{t+1}, θ^t) ≥ Q(θ^t, θ^t)

这使得M步可以使用梯度上升等迭代方法，适用于Q函数难以解析最大化的情况。

三、实战案例

3.1 GMM聚类可视化

python

import numpy as np

def gmm_clustering_demo():
    """GMM聚类演示"""
    np.random.seed(42)
    
    # 生成三个不同形状的高斯簇
    X1 = np.random.randn(80, 2) @ np.array([[1, 0.5], [0.5, 1]]) + np.array([3, 3])
    X2 = np.random.randn(80, 2) @ np.array([[2, 0], [0, 0.3]]) + np.array([-3, -2])
    X3 = np.random.randn(80, 2) + np.array([0, -3])
    X = np.vstack([X1, X2, X3])
    
    # 拟合GMM
    gmm = GMM(n_components=3, n_iter=100)
    gmm.fit(X)
    
    # 聚类结果
    labels = gmm.predict(X)
    
    print("GMM聚类结果:")
    print(f"  混合权重: {gmm.pi}")
    print(f"  均值:")
    for k in range(3):
        print(f"    分量{k}: {gmm.mu[k]}")
    
    # 计算聚类准确率（需要匹配标签）
    true_labels = np.array([0]*80 + [1]*80 + [2]*80)
    from itertools import permutations
    best_acc = 0
    for perm in permutations(range(3)):
        mapped = np.zeros_like(labels)
        for i, p in enumerate(perm):
            mapped[labels == i] = p
        acc = np.mean(mapped == true_labels)
        best_acc = max(best_acc, acc)
    print(f"  聚类准确率: {best_acc:.4f}")

gmm_clustering_demo()

3.2 初始化对EM的影响

python

import numpy as np

def initialization_sensitivity():
    """演示EM对初始值的敏感性"""
    np.random.seed(42)
    
    X1 = np.random.randn(100, 2) + np.array([3, 3])
    X2 = np.random.randn(100, 2) + np.array([-3, -3])
    X = np.vstack([X1, X2])
    
    results = []
    for seed in range(5):
        np.random.seed(seed)
        gmm = GMM(n_components=2, n_iter=50)
        gmm.fit(X)
        final_ll = gmm.log_likelihoods[-1]
        results.append((seed, final_ll))
    
    print("不同初始化下的最终对数似然:")
    for seed, ll in results:
        print(f"  seed={seed}: {ll:.4f}")
    
    print("\n不同的初始化可能导致不同的局部最优")
    print("实际中常用K-means++初始化或多次运行取最优")

initialization_sensitivity()

四、本章小结

第9章介绍了EM算法及其推广：

隐变量：无法直接观测的变量使极大似然估计变得困难。

EM算法：交替执行E步（推断隐变量后验）和M步（更新参数），保证对数似然单调递增。

GMM：混合高斯模型是EM最重要的应用之一，可以捕捉复杂的数据分布。

收敛性：EM保证收敛到局部最优，但不保证全局最优。

EM与K-means：GMM是K-means的概率化推广。

GEM：放宽M步要求，适用于Q函数难以解析最大化的情况。

EM算法的核心启示是：当直接优化困难时，引入隐变量将问题分解为交替优化的子问题。

五、思考与练习

思考题

为什么EM算法不能保证找到全局最优？如何缓解这个问题？

EM算法的E步和M步各自在做什么？能否用一句话概括？

GMM中的混合权重 π_k 有什么约束？为什么需要这个约束？

EM算法与K-means有什么联系和区别？

在什么情况下应该使用GMM而不是K-means？

编程练习

从零实现GMM的EM算法，在二维数据上进行聚类。

对比不同初始化下EM的收敛结果，验证其对初始值的敏感性。

在鸢尾花数据集上对比GMM和K-means的聚类效果。