高斯混合模型与 EM 算法

一句话定义

高斯混合模型（GMM）假设数据由 K 个高斯分量按混合系数混合生成，因「簇归属」隐变量不可观测而无法直接极大化似然；EM 算法用「E 步估计隐变量后验（责任）→ M 步固定责任更新参数」的交替迭代爬升似然。

为什么重要

EM 是含隐变量模型（GMM、HMM、主题模型、缺失数据）的通用引擎，「为什么不能一步算出、为什么交替会上升」的 ELBO 下界论证是概率机器学习的核心推理范式。GMM 相比 K-means 还给出软分配与簇形状（协方差），是密度估计与异常检测（对数似然低 = 异常）的正规方法。

前置知识

K-means 聚类与 k-means++ 的 K-means（本节给出其概率推广）；学习问题的形式化：假设空间、经验风险与泛化 的极大似然；条件概率与贝叶斯公式。

核心概念

• 混合分布：p(x) = ∑_k=1^K π_k 𝒩(x | μ_k, Σ_k)，π_k 为混合系数（和为 1）。
• 隐变量：z ∈ {1, …, K}，样本来自哪个分量（不可观测）。
• 责任（responsibility）：γ_ik = P(z_i = k | x_i)，E 步输出的软归属。
• E 步/M 步：估计后验 / 固定后验更新 (π, μ, Σ)。
• ELBO：证据下界，EM 每轮令其不降，从而似然单调上升。
• 模型选择：分量数 K 用 BIC/AIC 选，避免纯靠似然（分量越多似然越高）。

直观类比

听多个人同时说话的录音想分开每人的声纹：你不知道每段声音是谁的（隐变量），只好先猜一个归属方案（E 步），按猜测把每人的声纹估计出来（M 步），再用新声纹修正归属——来回迭代直到自洽。K-means 是「每段声音非此即彼」的武断版，GMM 是「七成可能甲、三成可能乙」的温和版。

原理与机制

为什么需要 EM：似然 ∏_i ∑_k π_k 𝒩(x_i | μ_k, Σ_k) 中 log 里套着求和，直接求导无解析解。EM 的策略是引入隐变量使 log 与求和换位：对任意分布 q(z)，log p(x) = ℒ(q, θ) + KL(q  ‖  p(z | x)) ≥ ℒ(q, θ)。E 步取 q = p(z | x, θ^old) 使 KL 为零（下界贴紧）；M 步固定 q 关于 θ 最大化下界——两步都让 log p(x) 不降，迭代保证单调爬升到局部最优。

GMM 的具体两步：E 步按贝叶斯公式算责任

M 步的更新恰似加权版高斯估计：μ_k = ∑_i γ_ik x_iN_k、Σ_k = ∑_i γ_ik(x_i - μ_k)(x_i - μ_k)^⊤N_k、π_k = N_kn，其中 N_k = ∑_i γ_ik。责任全为 0/1 时退化为 K-means——两者的谱系关系在此合拢。

实践要点：协方差需正则（reg_covar）防奇异；对初始化敏感，跑多组随机起点（n_init）取似然最高者；EM 只保证局部最优，且可能收敛到「一个分量塌缩到单点」的退化解。

公式与推导

下界分解（推导骨架）：恒等式 log p(x | θ) = ℒ(q, θ) + KL(q  ‖  p(z | x, θ))，其中

由 KL 非负得 ℒ ≤ log p(x)。E 步选择使 KL 归零的 q；M 步在固定 q 下最大化 ℒ；两步串联：log p(x | θ^old) = ℒ^old ≤ ℒ(θ^new) ≤ log p(x | θ^new)——似然单调不减的完整链条。

图示

log p(x)
 ▲                    ● ← 局部最优（EM 可能停在这）
 │        ●──●       ╱
 │      ╱     ●────●
 │    ●╱  ↑每轮 E贴紧下界、M最大化
 │  ●╱
 │ ● ← 起点（依赖初始化）
 └────────────────────────▶ 参数 θ
 EM = 沿下界爬山的坐标上升

实例或案例

两个拉长的椭圆簇：GMM 学出各向异性协方差，K-means 切错：

from sklearn.datasets import make_blobs
from sklearn.mixture import GaussianMixture
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
import numpy as np

rng = np.random.default_rng(0)
A = rng.multivariate_normal([0, 0], [[6, 3], [3, 2]], 300)      # 相关的椭圆簇
B = rng.multivariate_normal([8, 6], [[2, -1.5], [-1.5, 2]], 300)
X = StandardScaler().fit_transform(np.vstack([A, B]))
gmm = GaussianMixture(2, covariance_type="full", n_init=5, random_state=0).fit(X)
km = KMeans(2, n_init="auto", random_state=0).fit(X)
truth = np.array([0] * 300 + [1] * 300)
print("GMM 一致率:", (gmm.predict(X) == truth).mean().round(3))
print("K-means 一致率:", (km.labels_ == truth).mean().round(3))
print("BIC 曲线:", [GaussianMixture(k, n_init=3, random_state=0).fit(X).bic(X).round(0)
                  for k in (1, 2, 3, 4)])

常见误区

• 用训练似然选 K：分量越多似然越高，必然选出 K = n；必须用 BIC/AIC 或留出似然。
• 以为 EM 有全局保证：只保证单调到局部极值，初始化与多重启（n_init）是必须动作。
• 奇异协方差未处理：某分量塌缩到单点时 Σ → 0、似然爆炸，需 reg_covar 或 covariance_type="tied"。

与其他知识点的关系

K-means 聚类与 k-means++ 是本节的硬分配特例；学习问题的形式化：假设空间、经验风险与泛化 的极大似然是 M 步的目标来源；ELBO 思想是变分推断（深度生成模型 VAE）的起点，深度学习导览：边界与去向 有概念桥接；异常检测可直接用 GMM 的逐点对数似然。

自测题

1. 写出 GMM 的 E 步公式并解释「责任」含义。

- 要点：γ_ik = π_k 𝒩(x_i | μ_k, Σ_k)∑_j π_j 𝒩(x_i | μ_j, Σ_j)，即当前参数下样本 i 来自分量 k 的后验概率（软归属）。

1. 为什么 EM 每轮似然单调不降？仍有什么不保证？

- 要点：E 步把下界贴紧当前似然，M 步抬高下界；KL 项非负保证夹逼。不保证全局最优、不保证分量不退化。

1. K-means 与 GMM 的精确对应关系？

- 要点：GMM 在 Σ_k = ε I、ε → 0、π_k 相等时的硬分配 EM；「最近质心」即责任取 0/1 的极限。

延伸阅读

Dempster 等， Maximum Likelihood from Incomplete Data via the EM Algorithm（JRSS-B, 1977）；Bishop《PRML》第 9 章。