第9章：混合模型和EM算法（Mixture Models and EM）
一、导读
9.1 本章学习目标
第9章介绍了混合模型和EM算法——处理包含隐变量问题的强大工具。混合模型假设数据由多个概率分布（组件）混合生成，每个数据点属于哪个组件是未知的（隐变量）。EM算法提供了一种优雅的迭代方法来最大化含有隐变量的模型的似然函数。
本章的核心学习目标包括：
- 理解混合模型的基本概念：什么是混合模型？为什么需要混合模型？
- 掌握高斯混合模型（GMM）：GMM的数学形式、参数估计和应用。
- 深入理解EM算法：EM算法的推导、收敛性保证和几何直觉。
- 掌握EM算法的E步和M步：每个步骤具体在做什么？为什么交替执行能提升似然？
- 理解混合模型的应用：聚类、密度估计、软分类的具体实现。
- 了解EM算法的推广：变分EM、随机EM等变体。
9.2 为什么本章如此重要
混合模型和EM算法是机器学习中最重要的工具之一。在聚类分析中，GMM提供了比K-means更灵活的软聚类方法；在密度估计中，混合模型可以逼近任意复杂的分布；在缺失数据处理中，EM算法是标准的优化方法。
更重要的是，EM算法的思想——通过引入隐变量简化优化问题——贯穿了整个机器学习。从变分推断（第10章）到深度生成模型（VAE、GMM），EM的思想无处不在。理解EM算法，就掌握了处理隐变量问题的核心方法论。
9.3 核心问题引导
在阅读本章时，请带着以下问题思考：
混合模型与单个分布相比有什么优势？它能表示什么样的分布？
EM算法为什么能保证似然函数单调递增？收敛性如何保证？
E步和M步各自在做什么？为什么交替执行比直接优化更容易？
GMM中的隐变量有什么物理含义？后验概率$\gamma(z_{nk})$代表什么？
EM算法的局部最优问题如何解决？有哪些实用的初始化策略？
如何确定混合模型中组件的数量K？有哪些模型选择方法？
二、核心概念详解
9.1 混合模型
9.1.1 基本概念
混合模型假设数据由 $K$ 个概率分布（组件）混合生成：
$$p(\mathbf{x}) = \sum_{k=1}^{K} \pi_k p(\mathbf{x}|C_k)$$
其中 $\pi_k$ 是混合系数（先验概率），满足 $\sum_k \pi_k = 1$，$p(\mathbf{x}|C_k)$ 是第 $k$ 个组件的分布。
混合模型可以表示任意复杂的分布，只要组件数量足够多。
9.1.2 高斯混合模型（GMM）
当组件分布为高斯分布时，得到高斯混合模型：
$$p(\mathbf{x}) = \sum_{k=1}^{K} \pi_k \mathcal{N}(\mathbf{x}|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$$
参数为 $\boldsymbol{\theta} = \{\boldsymbol{\mu}_1, \ldots, \boldsymbol{\mu}_K, \boldsymbol{\Sigma}_1, \ldots, \boldsymbol{\Sigma}_K, \pi_1, \ldots, \pi_K\}$。
GMM可以表示：
- 多峰分布（多个高斯中心）
- 各向异性分布（不同的协方差结构）
- 重叠的簇（软聚类）
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
def generate_gmm_data(n_samples=500, random_state=42):
"""生成GMM数据"""
np.random.seed(random_state)
# 3个高斯组件
means = [np.array([0, 0]), np.array([3, 3]), np.array([6, 0])]
covs = [
np.array([[1, 0.5], [0.5, 1]]),
np.array([[1.5, -0.3], [-0.3, 0.5]]),
np.array([[0.8, 0], [0, 0.8]])
]
weights = [0.3, 0.4, 0.3]
samples = []
labels = []
for k in range(3):
n_k = int(n_samples * weights[k])
x_k = np.random.multivariate_normal(means[k], covs[k], n_k)
samples.append(x_k)
labels.extend([k] * n_k)
X = np.vstack(samples)
y = np.array(labels)
return X, y, means, covs, weights
X_gmm, y_gmm, true_means, true_covs, true_weights = generate_gmm_data()
# 可视化
fig, ax = plt.subplots(figsize=(8, 6))
scatter = ax.scatter(X_gmm[:, 0], X_gmm[:, 1], c=y_gmm, cmap='viridis', alpha=0.5, s=10)
for k in range(3):
ax.scatter(true_means[k][0], true_means[k][1], c='red', s=200, marker='*',
zorder=5, edgecolors='black')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_title('GMM Generated Data')
plt.tight_layout()
plt.savefig('gmm_data.png', dpi=150)
plt.close()9.2 EM算法
9.2.1 隐变量视角
引入隐变量 $\mathbf{z}$，其中 $z_k \in \{0, 1\}$ 且 $\sum_k z_k = 1$，表示数据点属于哪个组件。
联合分布：
$$p(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta}) = \prod_{k=1}^{K} (\pi_k p(\mathbf{x}|C_k))^{z_k}$$
观测数据的对数似然：
$$\ln p(\mathbf{X}|\boldsymbol{\theta}) = \sum_{n=1}^{N} \ln \left\{\sum_{k=1}^{K} \pi_k p(\mathbf{x}_n|C_k)\right\}$$
注意：对数里面有求和，直接优化很困难。
9.2.2 EM算法的推导
EM算法通过引入隐变量的后验分布 $q(\mathbf{Z})$，构造对数似然的下界：
$$\ln p(\mathbf{X}|\boldsymbol{\theta}) = \mathcal{L}(q, \boldsymbol{\theta}) + KL(q || p)$$
其中：
- $\mathcal{L}(q, \boldsymbol{\theta}) = \sum_{\mathbf{Z}} q(\mathbf{Z}) \ln \frac{p(\mathbf{X}, \mathbf{Z}|\boldsymbol{\theta})}{q(\mathbf{Z})}$ 是证据下界（ELBO）
- $KL(q || p) = -\sum_{\mathbf{Z}} q(\mathbf{Z}) \ln \frac{p(\mathbf{Z}|\mathbf{X}, \boldsymbol{\theta})}{q(\mathbf{Z})}$ 是KL散度
由于 $KL(q || p) \geq 0$，所以 $\mathcal{L}(q, \boldsymbol{\theta}) \leq \ln p(\mathbf{X}|\boldsymbol{\theta})$。
EM算法交替执行两步：
E步：固定 $\boldsymbol{\theta}$，优化 $q$。当 $q(\mathbf{Z}) = p(\mathbf{Z}|\mathbf{X}, \boldsymbol{\theta}^{old})$ 时，KL散度为零，下界最紧。
M步：固定 $q$，优化 $\boldsymbol{\theta}$。最大化 $\mathcal{L}(q, \boldsymbol{\theta})$。
9.2.3 GMM的EM算法
对于GMM，EM算法的具体形式为：
E步：计算后验概率（责任度）
$$\gamma(z_{nk}) = \frac{\pi_k \mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}{\sum_{j=1}^{K} \pi_j \mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)}$$
M步：更新参数
$$N_k = \sum_{n=1}^{N} \gamma(z_{nk})$$
$$\boldsymbol{\mu}_k^{new} = \frac{1}{N_k} \sum_{n=1}^{N} \gamma(z_{nk}) \mathbf{x}_n$$
$$\boldsymbol{\Sigma}_k^{new} = \frac{1}{N_k} \sum_{n=1}^{N} \gamma(z_{nk}) (\mathbf{x}_n - \boldsymbol{\mu}_k^{new})(\mathbf{x}_n - \boldsymbol{\mu}_k^{new})^T$$
$$\pi_k^{new} = \frac{N_k}{N}$$
class GaussianMixtureModel:
"""高斯混合模型的EM算法实现"""
def __init__(self, n_components=3, max_iter=100, tol=1e-6):
self.K = n_components
self.max_iter = max_iter
self.tol = tol
def fit(self, X):
N, D = X.shape
# 初始化参数
np.random.seed(42)
indices = np.random.choice(N, self.K, replace=False)
self.means = X[indices].copy()
self.covs = np.array([np.eye(D) for _ in range(self.K)])
self.weights = np.ones(self.K) / self.K
log_likelihoods = []
for iteration in range(self.max_iter):
# E步：计算责任度
responsibilities = np.zeros((N, self.K))
for k in range(self.K):
responsibilities[:, k] = self.weights[k] * \
multivariate_normal.pdf(X, mean=self.means[k], cov=self.covs[k] + 1e-6*np.eye(D))
# 归一化
total = responsibilities.sum(axis=1, keepdims=True)
responsibilities /= total + 1e-300
# 计算对数似然
log_likelihood = np.sum(np.log(total + 1e-300))
log_likelihoods.append(log_likelihood)
# M步：更新参数
Nk = responsibilities.sum(axis=0)
for k in range(self.K):
# 更新均值
self.means[k] = (responsibilities[:, k:k+1].T @ X / Nk[k]).ravel()
# 更新协方差
diff = X - self.means[k]
self.covs[k] = (responsibilities[:, k:k+1] * diff).T @ diff / Nk[k]
self.covs[k] += 1e-6 * np.eye(D) # 正则化
# 更新权重
self.weights[k] = Nk[k] / N
# 检查收敛
if len(log_likelihoods) > 1:
if abs(log_likelihoods[-1] - log_likelihoods[-2]) < self.tol:
print(f"收敛于第 {iteration+1} 次迭代")
break
self.log_likelihoods = log_likelihoods
return self
def predict(self, X):
"""预测每个样本最可能的组件"""
N = len(X)
responsibilities = np.zeros((N, self.K))
for k in range(self.K):
responsibilities[:, k] = self.weights[k] * \
multivariate_normal.pdf(X, mean=self.means[k], cov=self.covs[k] + 1e-6*np.eye(D))
return np.argmax(responsibilities, axis=1)
def predict_proba(self, X):
"""预测每个样本属于各组件的概率"""
N = len(X)
responsibilities = np.zeros((N, self.K))
for k in range(self.K):
responsibilities[:, k] = self.weights[k] * \
multivariate_normal.pdf(X, mean=self.means[k], cov=self.covs[k] + 1e-6*np.eye(D))
responsibilities /= responsibilities.sum(axis=1, keepdims=True)
return responsibilities
def score(self, X):
"""计算对数似然"""
N = len(X)
likelihood = np.zeros(N)
for k in range(self.K):
likelihood += self.weights[k] * \
multivariate_normal.pdf(X, mean=self.means[k], cov=self.covs[k] + 1e-6*np.eye(D))
return np.sum(np.log(likelihood + 1e-300))
# 训练GMM
gmm = GaussianMixtureModel(n_components=3, max_iter=100)
gmm.fit(X_gmm)
print(f"估计的均值:")
for k in range(3):
print(f" 组件{k+1}: {np.round(gmm.means[k], 2)}")
print(f"估计的权重: {np.round(gmm.weights, 3)}")
print(f"最终对数似然: {gmm.log_likelihoods[-1]:.2f}")9.3 EM算法的收敛性
9.3.1 单调性保证
EM算法保证每次迭代后对数似然不减：
$$\ln p(\mathbf{X}|\boldsymbol{\theta}^{(t+1)}) \geq \ln p(\mathbf{X}|\boldsymbol{\theta}^{(t)})$$
证明思路：
E步使下界 $\mathcal{L}$ 最紧（KL散度为零）
M步提升下界 $\mathcal{L}$
由于 $\ln p = \mathcal{L} + KL \geq \mathcal{L}$，对数似然也提升
9.3.2 收敛到局部最优
EM算法保证收敛到对数似然的局部极大值（或鞍点），不保证全局最优。
实际影响：
- 初始化对最终结果影响很大
- 多次随机初始化取最优是常用策略
- K-means++初始化可用于GMM
def em_convergence_demo():
"""演示EM算法的收敛过程"""
np.random.seed(42)
# 不同初始化下的收敛
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左图：对数似然收敛曲线
for seed in [42, 123, 456]:
gmm_temp = GaussianMixtureModel(n_components=3, max_iter=100)
np.random.seed(seed)
gmm_temp.fit(X_gmm)
axes[0].plot(gmm_temp.log_likelihoods, label=f'Seed={seed}', alpha=0.8)
axes[0].set_xlabel('Iteration')
axes[0].set_ylabel('Log Likelihood')
axes[0].set_title('EM Convergence (Different Initializations)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 右图：不同K值的选择
K_values = range(1, 8)
bic_scores = []
aic_scores = []
for K in K_values:
gmm_temp = GaussianMixtureModel(n_components=K, max_iter=100)
gmm_temp.fit(X_gmm)
N, D = X_gmm.shape
n_params = K * (D * (D + 1) / 2 + D) + K - 1
log_lik = gmm_temp.log_likelihoods[-1]
bic = -2 * log_lik + n_params * np.log(N)
aic = -2 * log_lik + 2 * n_params
bic_scores.append(bic)
aic_scores.append(aic)
axes[1].plot(K_values, bic_scores, 'o-', label='BIC', linewidth=2)
axes[1].plot(K_values, aic_scores, 's-', label='AIC', linewidth=2)
axes[1].set_xlabel('Number of Components K')
axes[1].set_ylabel('Score (lower is better)')
axes[1].set_title('Model Selection')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('em_convergence.png', dpi=150)
plt.close()
best_K_bic = K_values[np.argmin(bic_scores)]
print(f"BIC选择的最优K: {best_K_bic}")
em_convergence_demo()9.4 EM与K-means的关系
9.4.1 K-means作为GMM的特殊情况
K-means可以理解为GMM的硬版本（Hard EM）：
- GMM：软分配，每个样本以概率属于各组件
- K-means：硬分配，每个样本确定性地属于最近的中心
数学上，当GMM的所有协方差矩阵为 $\sigma^2 \mathbf{I}$ 且 $\sigma^2 \to 0$ 时，EM退化为K-means。
def compare_gmm_kmeans():
"""比较GMM和K-means"""
from sklearn.cluster import KMeans
# GMM
gmm = GaussianMixtureModel(n_components=3)
gmm.fit(X_gmm)
labels_gmm = gmm.predict(X_gmm)
# K-means
kmeans = KMeans(n_clusters=3, random_state=42, n_init=10)
labels_km = kmeans.fit_predict(X_gmm)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
axes[0].scatter(X_gmm[:, 0], X_gmm[:, 1], c=labels_gmm, cmap='viridis', alpha=0.5, s=10)
axes[0].scatter(gmm.means[:, 0], gmm.means[:, 1], c='red', s=200, marker='*',
edgecolors='black', zorder=5)
axes[0].set_title('GMM (Soft Clustering)')
axes[0].grid(True, alpha=0.3)
axes[1].scatter(X_gmm[:, 0], X_gmm[:, 1], c=labels_km, cmap='viridis', alpha=0.5, s=10)
axes[1].scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1],
c='red', s=200, marker='*', edgecolors='black', zorder=5)
axes[1].set_title('K-means (Hard Clustering)')
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('gmm_vs_kmeans.png', dpi=150)
plt.close()
print("GMM vs K-means:")
print(" GMM: 软聚类，考虑簇的形状和方向")
print(" K-means: 硬聚类，假设簇是球形的")
print(" GMM更灵活，K-means更简单快速")
compare_gmm_kmeans()9.5 混合模型的应用
9.5.1 密度估计
混合模型可以逼近任意复杂的概率密度函数。
9.5.2 软聚类
GMM提供每个样本属于各簇的概率，而不仅仅是硬分配。
9.5.3 异常检测
低概率的样本可能是异常值。
def anomaly_detection_with_gmm():
"""使用GMM进行异常检测"""
np.random.seed(42)
# 正常数据
N_normal = 400
X_normal = np.vstack([
np.random.multivariate_normal([0, 0], [[1, 0], [0, 1]], N_normal//2),
np.random.multivariate_normal([3, 3], [[0.5, 0], [0, 0.5]], N_normal//2)
])
# 训练GMM
gmm = GaussianMixtureModel(n_components=2)
gmm.fit(X_normal)
# 生成测试数据（含异常值）
N_test = 200
X_test = np.vstack([
np.random.multivariate_normal([0, 0], [[1, 0], [0, 1]], N_test//2),
np.random.multivariate_normal([3, 3], [[0.5, 0], [0, 0.5]], N_test//4),
np.random.uniform(-5, 8, (N_test//4, 2)) # 异常值
])
# 计算每个样本的对数概率
log_probs = np.zeros(len(X_test))
for i in range(len(X_test)):
prob = 0
for k in range(2):
prob += gmm.weights[k] * multivariate_normal.pdf(
X_test[i], mean=gmm.means[k], cov=gmm.covs[k] + 1e-6*np.eye(2))
log_probs[i] = np.log(prob + 1e-300)
# 异常检测
threshold = np.percentile(log_probs, 5) # 最低5%为异常
is_anomaly = log_probs < threshold
# 可视化
fig, ax = plt.subplots(figsize=(8, 6))
normal_mask = ~is_anomaly
ax.scatter(X_test[normal_mask, 0], X_test[normal_mask, 1], c='blue', alpha=0.5,
s=10, label='Normal')
ax.scatter(X_test[is_anomaly, 0], X_test[is_anomaly, 1], c='red', s=30,
marker='x', label='Anomaly')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_title('GMM-based Anomaly Detection')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('anomaly_detection.png', dpi=150)
plt.close()
n_detected = np.sum(is_anomaly[N_test//4*3:]) # 后1/4是异常值
n_total_anomaly = N_test // 4
print(f"异常值检测: {n_detected}/{n_total_anomaly} 个异常值被识别")
anomaly_detection_with_gmm()三、实战案例
9.1 完整的GMM密度估计
def gmm_density_estimation():
"""GMM密度估计完整示例"""
np.random.seed(42)
# 生成复杂的2D分布
N = 500
data = np.vstack([
np.random.multivariate_normal([-2, -2], [[0.5, 0], [0, 0.5]], N//4),
np.random.multivariate_normal([2, 2], [[0.3, 0.2], [0.2, 0.3]], N//4),
np.random.multivariate_normal([-2, 2], [[0.8, -0.3], [-0.3, 0.4]], N//4),
np.random.multivariate_normal([2, -2], [[0.4, 0], [0, 0.6]], N//4)
])
# 训练不同K的GMM
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for idx, K in enumerate([2, 4, 8]):
gmm = GaussianMixtureModel(n_components=K, max_iter=100)
gmm.fit(data)
# 绘制密度等值线
x_grid = np.linspace(-5, 5, 100)
y_grid = np.linspace(-5, 5, 100)
xx, yy = np.meshgrid(x_grid, y_grid)
grid_points = np.c_[xx.ravel(), yy.ravel()]
density = np.zeros(len(grid_points))
for k in range(K):
density += gmm.weights[k] * multivariate_normal.pdf(
grid_points, mean=gmm.means[k], cov=gmm.covs[k] + 1e-6*np.eye(2))
density = density.reshape(xx.shape)
axes[idx].contourf(xx, yy, density, levels=20, cmap='viridis', alpha=0.7)
axes[idx].scatter(data[:, 0], data[:, 1], c='white', s=3, alpha=0.3)
axes[idx].scatter(gmm.means[:, 0], gmm.means[:, 1], c='red', s=100,
marker='*', edgecolors='black', zorder=5)
axes[idx].set_title(f'K={K}')
axes[idx].set_xlim(-5, 5)
axes[idx].set_ylim(-5, 5)
plt.tight_layout()
plt.savefig('gmm_density_estimation.png', dpi=150)
plt.close()
gmm_density_estimation()9.2 图像压缩
def gmm_image_compression():
"""使用GMM进行图像颜色量化（压缩）"""
np.random.seed(42)
# 创建简单的渐变图像
L = 50
image = np.zeros((L, L, 3))
for i in range(L):
for j in range(L):
image[i, j] = [i/L*255, j/L*255, (i+j)/(2*L)*255]
# 将图像像素展平
pixels = image.reshape(-1, 3).astype(float)
# 使用GMM进行颜色量化
K = 8
gmm = GaussianMixtureModel(n_components=K, max_iter=50)
gmm.fit(pixels)
# 量化
labels = gmm.predict(pixels)
quantized = gmm.means[labels]
# 计算压缩比
original_size = L * L * 3 * 8 # 每像素3通道，每通道8位
compressed_size = K * 3 * 8 + L * L * np.ceil(np.log2(K)) # K个颜色中心 + 索引
print(f"原始大小: {original_size} bits")
print(f"压缩后大小: {compressed_size:.0f} bits")
print(f"压缩比: {original_size/compressed_size:.1f}x")
# 可视化
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
axes[0].imshow(image.astype(np.uint8))
axes[0].set_title('Original')
axes[0].set_xticks([])
axes[0].set_yticks([])
quantized_image = quantized.reshape(L, L, 3).astype(np.uint8)
axes[1].imshow(quantized_image)
axes[1].set_title(f'GMM Quantized (K={K})')
axes[1].set_xticks([])
axes[1].set_yticks([])
plt.tight_layout()
plt.savefig('gmm_compression.png', dpi=150)
plt.close()
gmm_image_compression()9.3 EM算法的几何直觉
def em_geometric_intuition():
"""EM算法的几何直觉可视化"""
np.random.seed(42)
# 生成简单数据
N = 100
X_em = np.vstack([
np.random.randn(N//2, 2) * 0.5 + [0, 0],
np.random.randn(N//2, 2) * 0.5 + [3, 3]
])
# 运行EM并记录每步的决策边界
gmm_vis = GaussianMixtureModel(n_components=2, max_iter=5)
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# 初始状态
gmm_vis.means = np.array([[0.5, 0.5], [2.5, 2.5]])
gmm_vis.covs = np.array([np.eye(2), np.eye(2)])
gmm_vis.weights = np.array([0.5, 0.5])
for iteration in range(3):
# E步
resp = np.zeros((N, 2))
for k in range(2):
resp[:, k] = gmm_vis.weights[k] * \
multivariate_normal.pdf(X_em, mean=gmm_vis.means[k], cov=gmm_vis.covs[k])
resp /= resp.sum(axis=1, keepdims=True)
# 可视化
ax = axes[iteration]
colors = resp[:, 0]
scatter = ax.scatter(X_em[:, 0], X_em[:, 1], c=colors, cmap='RdBu',
alpha=0.6, s=20)
ax.scatter(gmm_vis.means[:, 0], gmm_vis.means[:, 1], c='black', s=200,
marker='*', zorder=5)
ax.set_title(f'Iteration {iteration+1}')
ax.set_xlim(-2, 5)
ax.set_ylim(-2, 5)
ax.grid(True, alpha=0.3)
# M步
Nk = resp.sum(axis=0)
for k in range(2):
gmm_vis.means[k] = (resp[:, k:k+1].T @ X_em / Nk[k]).ravel()
diff = X_em - gmm_vis.means[k]
gmm_vis.covs[k] = (resp[:, k:k+1] * diff).T @ diff / Nk[k]
gmm_vis.covs[k] += 1e-6 * np.eye(2)
gmm_vis.weights[k] = Nk[k] / N
plt.tight_layout()
plt.savefig('em_geometric.png', dpi=150)
plt.close()
print("EM算法的几何直觉:")
print(" E步: 根据当前参数计算每个样本的责任度（颜色渐变表示）")
print(" M步: 根据责任度更新参数（均值移动到加权中心）")
print(" 交替执行直到收敛")
em_geometric_intuition()四、本章小结
第9章系统介绍了混合模型和EM算法：
混合模型：假设数据由多个概率分布混合生成，可以表示任意复杂的分布。高斯混合模型（GMM）是最常用的混合模型。
EM算法：通过引入隐变量，将困难的优化问题分解为交替的E步和M步。E步计算隐变量的后验分布，M步最大化完全数据的对数似然。
收敛性：EM算法保证对数似然单调递增，收敛到局部极大值。初始化对结果影响很大，多次随机初始化是常用策略。
与K-means的关系：K-means是GMM的硬版本（Hard EM），当协方差趋向零时EM退化为K-means。
应用：密度估计、软聚类、异常检测、图像压缩等。GMM提供了比K-means更灵活的聚类方法。
EM算法是处理隐变量问题的核心方法论。它的思想——通过引入隐变量简化优化——在变分推断、深度生成模型等现代方法中得到了广泛的推广和应用。
五、思考与练习
思考题
为什么EM算法不能保证找到全局最优？有哪些方法可以缓解局部最优问题？
GMM的协方差矩阵退化为零时会发生什么？如何防止这种情况？
如何使用BIC/AIC选择GMM的组件数K？它们的理论依据是什么？
EM算法在什么情况下收敛很慢？有哪些加速方法？
混合模型与核密度估计相比有什么优缺点？
编程练习
带约束的GMM：实现球形协方差（$\Sigma_k = \sigma_k^2 \mathbf{I}$）和对角协方差的GMM，比较不同约束下的聚类效果。
EM初始化策略：比较随机初始化、K-means初始化和K-means++初始化对EM收敛速度和最终结果的影响。
GMM用于说话人识别：使用GMM对不同说话人的语音特征建模，实现简单的说话人识别系统。