1.1 极大似然估计

在进入EM算法之前，先快速回顾一下极大似然估计（Maximum Likelihood Estimation, MLE） 。

假设我们有一个概率模型，其参数为 $\theta$，观测数据为 $X = \{x_1, x_2, ..., x_N\}$。MLE的目标是找到一组参数 $\theta$，使得观测数据出现的概率最大。

似然函数为：

为了计算方便，通常取对数，得到对数似然函数：

MLE的估计值为：

当模型简单时（如单个高斯分布、二项分布等），我们可以直接对 $\ell(\theta)$ 求导并令导数为零，得到解析解。

1.2 隐变量带来的困境

现在考虑一个更复杂的情况：数据中存在隐变量 $Z = \{z_1, z_2, ..., z_N\}$，我们观测不到 $Z$，只能观测到 $X$。

此时，似然函数变成：

或者对于连续隐变量：

对数似然为：

注意看：对数和求和的位置颠倒了。

• 在简单情况下，对 $\log p(X|\theta)$ 求导，其中 $p(X|\theta)$ 是乘积形式，取对数后变成求和，求导很方便。
• 但现在，$\log$ 外面是一个求和（对隐变量求和），而这个求和又无法移到 $\log$ 外面。对数里面是求和，求导变得极其困难。

这就是隐变量带来的核心困境：边际似然函数中的“和的对数”无法直接优化。

1.3 一个直观的例子：二硬币模型

让我们用一个经典例子来感受这个困境。

假设有两枚硬币A和B，它们抛出正面的概率分别是 $\theta_A$ 和 $\theta_B$，但我们不知道。我们做了5轮实验，每轮随机选择一枚硬币，抛10次，记录正面次数。

情况一：我们知道每轮用的是哪枚硬币

这种情况下，没有隐变量。我们可以直接估计：

情况二：我们不知道每轮用的是哪枚硬币

现在，“每轮用的是哪枚硬币”就是隐变量 $Z$。我们只知道观测数据 $X$（每轮的正面次数），却不知道 $Z$。

如果我们想用MLE估计 $\theta_A$ 和 $\theta_B$，需要最大化：$\ell(\theta_A, \theta_B) = \log \sum_{Z} p(X, Z|\theta_A, \theta_B)$

这个式子中的求和（对所有可能的硬币分配组合）使得直接优化变得极其困难。这就是EM算法要解决的问题。

2.1 EM算法的引入

回到二硬币模型。我们陷入了一个循环困境：

• 如果知道 $\theta_A$ 和 $\theta_B$，就能推断每轮用的是哪枚硬币（即估计隐变量 $Z$）
• 如果知道每轮用的是哪枚硬币，就能估计 $\theta_A$ 和 $\theta_B$

EM算法的解决思路为先随便猜一个，然后交替迭代。

1. 随机初始化 $\theta_A$ 和 $\theta_B$（比如都设为0.5）
2. E步：基于当前的 $\theta$，计算隐变量 $Z$ 的概率分布（而不是硬性指定）
3. M步：基于E步得到的隐变量分布，用极大似然估计更新 $\theta$
4. 重复步骤2-3，直到收敛

2.2 E步：期望的计算

E步（Expectation Step，期望步骤）：基于当前的参数估计 $\theta^{(t)}$，计算完整数据对数似然 $\log p(X, Z|\theta)$ 在隐变量后验分布 $p(Z|X, \theta^{(t)})$ 下的期望。

数学上，E步构造如下函数：

也就是说：

这个 $Q$ 函数是EM算法的核心。它不是原始的似然函数，而是完整数据对数似然的期望。

2.3 M步：最大化的步骤

M步（Maximization Step，最大化步骤）：寻找使 $Q$ 函数最大化的新参数：

$Q$ 函数通常比原始的边际似然容易最大化——在 $Q$ 函数中，隐变量 $Z$ 已经被积分掉了（通过期望），剩下的只是关于 $\theta$ 的优化问题，而且往往有解析解。

2.4 EM算法的完整流程

• 输入：观测数据：$X$；隐变量模型：$p(X, Z \mid \theta)$；初始参数：$\theta^{(0)}$
• 输出：“参数估计：$\theta$
• 算法流程：
  1. 初始化：$t = 0$
  2. 重复以下步骤直至收敛： a. E步（期望步） —— 计算联合对数似然在隐变量后验分布下的期望：$Q(\theta, \theta^{(t)}) = \mathbb{E}_{Z \mid X, \theta^{(t)}} \left[ \log p(X, Z \mid \theta) \right]$ b. M步（最大化步） —— 更新参数，使 $Q$ 函数最大化：$\theta^{(t+1)} = \arg\max_{\theta} \, Q(\theta, \theta^{(t)})$ c. 迭代更新：$t = t + 1$
  3. 返回：$\theta^{(t)}$

3.1 Jensen不等式：建立下界

EM算法的数学基础是 Jensen不等式。对于凹函数 $f$（如 $\log$ 函数），Jensen不等式告诉我们：

现在，我们想最大化 $\log p(X|\theta)$。引入任意一个关于隐变量 $Z$ 的分布 $q(Z)$：

关于上述变换的说明

• $q(Z)$：任意一个关于隐变量 $Z$ 的概率分布（满足 $\sum_Z q(Z)=1$，且 $q(Z)>0$）。

• 数学技巧：原式是 $\sum_Z p(X,Z)$，乘以一个 $1 = \frac{q(Z)}{q(Z)}$，拆成了 $\sum_Z q(Z) \cdot \frac{p}{q}$。

• 目的：这一步是为了把求和伪装成数学期望。因为 $q(Z)$ 正好是一个概率权重，所以 $\sum_Z q(Z) \cdot [\cdot]$ 就等于“在 $q(Z)$ 这个分布下，对括号里的值求期望”。

由于 $\log$ 是凹函数，应用Jensen不等式：

右边这个量被称为 ELBO（Evidence Lower Bound） ——证据下界：

3.2 ELBO与KL散度的关系

ELBO与对数似然之间有一个重要的关系：

其中 $\text{KL}(q \| p)$ 是KL散度，衡量两个分布之间的差异，恒为非负。

ELBO与对数似然之间的重要关系数学推导

这个等式不是凭空冒出来的，它是贝叶斯公式的简单代数变形。

• 联合分布可以拆解为：$p(X, Z|\theta) = p(Z|X, \theta) \cdot p(X|\theta)$

• 两边同时取对数：$\log p(X, Z|\theta) = \log p(Z|X, \theta) + \log p(X|\theta)$

• 两边同时在分布 $q(Z)$ 下取期望 $\mathbb{E}_q[\cdot]$：$\mathbb{E}_q[\log p(X, Z|\theta)] = \mathbb{E}_q[\log p(Z|X, \theta)] + \log p(X|\theta)$

• 把 $\log p(X|\theta)$ 单独拎到左边，整理一下：$\log p(X|\theta) = \mathbb{E}_q[\log p(X, Z|\theta)] - \mathbb{E}_q[\log p(Z|X, \theta)]$

在右边巧妙地加一项、减一项 $\mathbb{E}_q[\log q(Z)]$（恒等变形）：

$\log p(X|\theta) = \mathbb{E}_q[\log \frac{p(X, Z|\theta)}{q(Z)}] + \mathbb{E}_q[\log \frac{q(Z)}{p(Z|X, \theta)}]$

第一项正好就是 ELBO，第二项正好就是 $KL(q \parallel p)$。

于是，得到ELBO与对数似然之间的重要关系式： $\boxed{\log p(X|\theta) = \text{ELBO}(q, \theta) + D_{KL}(q(Z) \parallel p(Z|X, \theta))}$

因此：

$\log p(X|\theta) \geq \text{ELBO}(q, \theta)$

当且仅当 $q(Z) = p(Z|X, \theta)$ 时，KL散度为0，ELBO等于对数似然。

3.3 EM算法的ELBO视角

从ELBO的视角来看，EM算法的两步实际上是在交替优化两个变量：

• E步：固定 $\theta$，选择 $q(Z) = p(Z|X, \theta)$，使ELBO等于对数似然（KL散度为0）
• M步：固定 $q$，最大化ELBO（等价于最大化 $Q$ 函数），提升对数似然

这个视角揭示了EM算法的本质：它通过交替优化ELBO，逐步提升对数似然的下界，从而间接提升对数似然本身。

3.4 收敛性证明

EM算法最优雅的性质是它保证每次迭代后观测数据的对数似然不会下降。即 $\boxed{\log p(X|\theta^{(t)}) \leq \log p(X|\theta^{(t+1)})}$

证明概要：

• 在E步中，选择 $q(Z) = p(Z|X, \theta^{(t)})$，此时ELBO等于对数似然：$\log p(X|\theta^{(t)}) = \text{ELBO}(q, \theta^{(t)})$

• 在M步中，最大化 $Q$ 函数（即ELBO）：$\text{ELBO}(q, \theta^{(t+1)}) \geq \text{ELBO}(q, \theta^{(t)})$

• 由ELBO的定义，对于任意 $\theta$：$\log p(X|\theta) \geq \text{ELBO}(q, \theta)$

• 特别地：$\log p(X|\theta^{(t+1)}) \geq \text{ELBO}(q, \theta^{(t+1)}) \geq \text{ELBO}(q, \theta^{(t)}) = \log p(X|\theta^{(t)})$

• 因此：$\log p(X|\theta^{(t+1)}) \geq \log p(X|\theta^{(t)})$

Important

EM算法保证收敛到局部最优（似然函数的局部最大值），但不保证收敛到全局最优。不同的初始值可能导致不同的局部最优解。

4.1 问题设定

回到二硬币模型。我们有5轮实验，每轮抛10次，但不知道用的是哪枚硬币：

目标是估计 $\theta_A$ 和 $\theta_B$。

4.2 EM算法求解

• 初始化：$\theta_A^{(0)} = 0.6$，$\theta_B^{(0)} = 0.5$

• 第1轮：

1. E步：对每一轮，计算它来自硬币A和硬币B的概率。

Tip

以第一轮（5正5反）为例：

• 如果是硬币A：$p = 0.6^5 \times 0.4^5 = 0.000796$
• 如果是硬币B：$p = 0.5^5 \times 0.5^5 = 0.000977$

归一化后：$P(A|\text{第一轮}) = \frac{0.000796}{0.000796 + 0.000977} \approx 0.45 \quad P(B|\text{第一轮}) = \frac{0.000977}{0.000796 + 0.000977} \approx 0.55$

同理计算其他轮次。

2. M步：基于E步得到的概率，计算每枚硬币的“期望正面次数”和“期望总次数”。

Tip

以硬币A为例，它“贡献”的正面次数 = 各轮正面次数 × 该轮来自A的概率之和。

计算得到新参数：$\theta_A^{(1)} = \frac{21.3}{21.3 + 8.6} \approx 0.71 \quad \theta_B^{(1)} = \frac{11.7}{11.7 + 8.4} \approx 0.58$

• 迭代：重复E步和M步，直到 $\theta$ 收敛。

4.3 二硬币模型的Python实现

1
import numpy as np
2

3
def em_coin(observations, theta_A, theta_B, max_iter=100, tol=1e-6):
4
    """
5
    二硬币模型的EM算法
6
    observations: 每轮实验的正面次数（假设每轮抛10次）
7
    theta_A, theta_B: 初始参数
8
    """
9
    n_rounds = len(observations)
10
    n_flips = 10  # 每轮抛10次
11

12
    for t in range(max_iter):
13
        # E步：计算每轮来自A和B的概率
14
        prob_A = np.zeros(n_rounds)
15
        prob_B = np.zeros(n_rounds)
16

17
        for i, heads in enumerate(observations):
18
            tails = n_flips - heads
19
            # 在给定参数下，出现该观测结果的概率
20
            p_A = (theta_A ** heads) * ((1 - theta_A) ** tails)
21
            p_B = (theta_B ** heads) * ((1 - theta_B) ** tails)
22
            prob_A[i] = p_A / (p_A + p_B)
23
            prob_B[i] = p_B / (p_A + p_B)
24

25
        # M步：更新参数
26
        # 硬币A的期望正面次数和期望总次数
27
        expected_heads_A = np.sum(prob_A * observations)
28
        expected_total_A = np.sum(prob_A * n_flips)
29
        theta_A_new = expected_heads_A / expected_total_A
30

31
        expected_heads_B = np.sum(prob_B * observations)
32
        expected_total_B = np.sum(prob_B * n_flips)
33
        theta_B_new = expected_heads_B / expected_total_B
34

35
        # 检查收敛
36
        if abs(theta_A_new - theta_A) < tol and abs(theta_B_new - theta_B) < tol:
37
            print(f"收敛于第 {t+1} 轮迭代")
38
            break
39

40
        theta_A, theta_B = theta_A_new, theta_B_new
41
        print(f"第 {t+1} 轮: θA={theta_A:.4f}, θB={theta_B:.4f}")
42

43
    return theta_A, theta_B
44

45
# 运行EM算法
46
observations = [5, 9, 4, 4, 5]  # 每轮正面次数
47
theta_A, theta_B = em_coin(observations, 0.6, 0.5)
48
print(f"\n最终结果: θA={theta_A:.4f}, θB={theta_B:.4f}")

5.1 高斯混合模型的定义

高斯混合模型（Gaussian Mixture Model, GMM） 是EM算法最经典的应用场景。

GMM假设数据由 $K$ 个高斯分布混合而成。每个数据点 $x_i$ 的生成过程是：

• 以概率 $\pi_k$ 选择一个高斯分量 $k$
• 从该高斯分布 $\mathcal{N}(\mu_k, \Sigma_k)$ 中采样一个点

其中 $\pi_k$ 是混合系数（$\sum_k \pi_k = 1$），$\mu_k$ 和 $\Sigma_k$ 是第 $k$ 个高斯分量的均值和协方差矩阵。

隐变量是每个数据点 $x_i$ 来自哪个高斯分量——这个信息无法观测。

5.2 GMM的EM算法推导

完整数据的对数似然为：

其中 $z_{ik} \in \{0, 1\}$ 表示第 $i$ 个点是否属于第 $k$ 个分量。

E步：计算后验概率

M步：更新参数

5.3 GMM的Python实现

1
import numpy as np
2
import matplotlib.pyplot as plt
3
from scipy.stats import multivariate_normal
4

5
class GMM:
6
    """高斯混合模型（使用EM算法）"""
7

8
    def __init__(self, n_components=3, max_iter=100, tol=1e-6):
9
        self.n_components = n_components
10
        self.max_iter = max_iter
11
        self.tol = tol
12
        self.pi = None      # 混合系数
13
        self.mu = None      # 均值
14
        self.sigma = None   # 协方差矩阵
15
        self.gamma = None   # 后验概率（责任）
16

17
    def _initialize(self, X):
18
        """初始化参数（使用K-Means++思想）"""
19
        n_samples, n_features = X.shape
20

21
        # 初始化均值：随机选择K个样本
22
        indices = np.random.choice(n_samples, self.n_components, replace=False)
23
        self.mu = X[indices].copy()
24

25
        # 初始化协方差：单位矩阵
26
        self.sigma = np.array([np.eye(n_features) for _ in range(self.n_components)])
27

28
        # 初始化混合系数：均匀分布
29
        self.pi = np.ones(self.n_components) / self.n_components
30

31
    def _e_step(self, X):
32
        """E步：计算后验概率"""
33
        n_samples = X.shape[0]
34
        self.gamma = np.zeros((n_samples, self.n_components))
35

36
        for k in range(self.n_components):
37
            # 计算每个点属于第k个分量的概率密度
38
            rv = multivariate_normal(mean=self.mu[k], cov=self.sigma[k])
39
            self.gamma[:, k] = self.pi[k] * rv.pdf(X)
40

41
        # 归一化
42
        self.gamma = self.gamma / np.sum(self.gamma, axis=1, keepdims=True)
43

44
    def _m_step(self, X):
45
        """M步：更新参数"""
46
        n_samples, n_features = X.shape
47

48
        # 更新混合系数
49
        N_k = np.sum(self.gamma, axis=0)
50
        self.pi = N_k / n_samples
51

52
        # 更新均值
53
        for k in range(self.n_components):
54
            self.mu[k] = np.sum(self.gamma[:, k:k+1] * X, axis=0) / N_k[k]
55

56
        # 更新协方差
57
        for k in range(self.n_components):
58
            diff = X - self.mu[k]
59
            weighted_diff = self.gamma[:, k:k+1] * diff
60
            self.sigma[k] = (weighted_diff.T @ diff) / N_k[k]
61
            # 加小值保证正定
62
            self.sigma[k] += 1e-6 * np.eye(n_features)
63

64
    def fit(self, X):
65
        """训练GMM"""
66
        self._initialize(X)
67

68
        prev_log_likelihood = -np.inf
69

70
        for t in range(self.max_iter):
71
            # E步
72
            self._e_step(X)
73

74
            # M步
75
            self._m_step(X)
76

77
            # 计算对数似然
78
            log_likelihood = 0
79
            for i in range(X.shape[0]):
80
                likelihood = 0
81
                for k in range(self.n_components):
82
                    rv = multivariate_normal(mean=self.mu[k], cov=self.sigma[k])
83
                    likelihood += self.pi[k] * rv.pdf(X[i])
84
                log_likelihood += np.log(likelihood + 1e-10)
85

86
            if abs(log_likelihood - prev_log_likelihood) < self.tol:
87
                print(f"GMM收敛于第 {t+1} 轮迭代")
88
                break
89

90
            prev_log_likelihood = log_likelihood
91

92
        return self
93

94
    def predict(self, X):
95
        """预测每个点最可能属于的簇"""
96
        self._e_step(X)
97
        return np.argmax(self.gamma, axis=1)
98

99
    def predict_proba(self, X):
100
        """预测每个点属于每个簇的概率"""
101
        self._e_step(X)
102
        return self.gamma
103

104
# 生成数据
105
np.random.seed(42)
106
n_samples = 300
107

108
# 三个高斯分量
109
X1 = np.random.multivariate_normal([0, 0], [[1, 0.5], [0.5, 1]], n_samples // 3)
110
X2 = np.random.multivariate_normal([5, 5], [[1, -0.3], [-0.3, 1]], n_samples // 3)
111
X3 = np.random.multivariate_normal([0, 5], [[0.5, 0], [0, 0.5]], n_samples // 3)
112
X = np.vstack([X1, X2, X3])
113

114
# 训练GMM
115
gmm = GMM(n_components=3, max_iter=100)
116
gmm.fit(X)
117
labels = gmm.predict(X)
118

119
# 可视化
120
plt.figure(figsize=(10, 6))
121
plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis', alpha=0.6)
122
plt.scatter(gmm.mu[:, 0], gmm.mu[:, 1], c='red', marker='X', s=200, label='簇中心')
123
plt.title('GMM聚类结果（EM算法）')
124
plt.legend()
125
plt.axis('equal')
126
plt.show()

6.1 K-Means是EM的硬版本

在后续的文章中，我们会详细介绍K-Means算法。我们可以从EM的视角审视它 —— K-Means可以看作是EM算法的一个特例。

关键区别：K-Means做的是硬分配（Hard Assignment） —— 每个数据点只能属于一个簇；而GMM-EM做的是软分配（Soft Assignment） —— 每个数据点可以以一定的概率属于多个簇。

从数学上看，K-Means等价于GMM的一个特殊情形：当所有高斯分量的协方差矩阵 $\Sigma_k = \sigma^2 I$（各向同性且相同）且 $\sigma^2 \to 0$ 时，后验概率 $\gamma_{ik}$ 退化为 one-hot 向量，GMM-EM 退化为 K-Means。

6.2 硬分类 vs 软分类

6.3 从硬到软：软分类的优势

在真实数据中，簇与簇之间往往存在重叠。一个数据点可能部分属于簇A、部分属于簇B——比如一个身高175cm的人，既有可能是男生也有可能是女生。

K-Means的硬分配无法表达这种不确定性。而GMM的软分配通过概率自然地表达了这种不确定性，这种软分配不仅更符合真实情况，还避免了硬分配中“边界点”的尴尬 —— 边界点不会被迫完全属于某一个簇，而是以概率分布在多个簇之间。

总结

概念	核心内容
隐变量困境	边际似然 $\log \sum_Z p(X,Z\|\theta)$ 中的“和的对数”难以直接优化
E步	计算 $Q(\theta,\theta^{(t)}) = \mathbb{E}_{Z\|X,\theta^{(t)}}[\log p(X,Z\|\theta)]$
M步	$\theta^{(t+1)} = \arg\max_\theta Q(\theta,\theta^{(t)})$
ELBO	$\log p(X
收敛性	EM保证每一步似然函数不下降，但只能收敛到局部最优
K-Means	EM的硬分配特例（$\Sigma_k = \sigma^2 I, \sigma^2 \to 0$）
GMM	EM的软分配经典应用，可处理任意椭圆形状的簇

核心要点回顾

1. 隐变量的困境：当模型包含隐变量时，边际似然 $\log \sum_Z p(X,Z|\theta)$ 中的“和的对数”无法直接优化。EM算法通过迭代的方式绕开了这个困难。
2. E步与M步：E步计算完整数据对数似然在隐变量后验分布下的期望（即 $Q$ 函数），M步最大化这个期望来更新参数。两步交替进行，逐步逼近最优解。
3. ELBO与收敛性：EM算法的数学基础是Jensen不等式，它保证ELBO始终是边际对数似然的下界。每一步M步都在提升ELBO，从而间接提升对数似然。EM算法保证收敛，但只能收敛到局部最优。
4. K-Means是EM的特例：K-Means相当于GMM-EM在硬分配极限下的退化版本。K-Means做硬分类（每个点100%属于一个簇），GMM做软分类（每个点以概率属于各簇）。软分类能更好地处理簇重叠和不确定性。
5. EM算法的本质：EM与其说是一种具体的算法，不如说是一种解决问题的框架。它不限定具体的模型——可以是GMM、HMM、LDA等——只要模型包含隐变量，都可以用EM框架来求解。