EM 算法完整深度讲解

1

EM算法是一种用于估计含有隐变量（未观测变量）模型的参数的迭代算法。

1
2
3
4
5

E = Expectation（期望）
M = Maximization（最大化）

整个算法就是在这两步之间反复循环
直到收敛为止

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

数据：观测到100个学生的身高

现象：身高分布呈现两个峰值
可能因为：样本中既有男生也有女生
        男生身高平均175cm
        女生身高平均163cm

问题：
我们知道有两个高斯分布混合
但我们不知道：
1. 每个学生是男生还是女生（隐变量Z）
2. 男生身高的方差是多少？
3. 女生身高的方差是多少？
4. 男女各占多少比例？

如果我们知道每个学生的性别，问题就很简单：
只需要分别对男生和女生计算均值和方差

但我们不知道！所以陷入了困境

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

完整的似然函数（如果知道Z）：
L(θ | X, Z) = ∏_i P(X_i | Z_i, θ)

这很容易最大化（就是简单的回归或MLE）

但观测数据的似然函数（不知道Z）：
L(θ | X) = ∏_i Σ_z P(X_i, Z_i=z | θ)
        = ∏_i Σ_z P(X_i | Z_i=z, θ)·P(Z_i=z | θ)

问题：
这个求和在指数中！
∏与Σ混合在一起，无法直接求导
无法得到闭形解（closed-form solution）

例子：
L = Σ_z [a·z + b·(1-z)]²

这与 [a·E[Z] + b·(1-E[Z])]² 完全不同！

所以不能简单地用条件期望替代
必须用迭代方法

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

与其直接最大化难的函数：
L(θ | X) = ∏_i Σ_z P(X_i | Z_i=z, θ)·P(Z_i=z | θ)

不如：
1. 猜测隐变量的值（或更准确地说，其分布）
2. 假设这个猜测是对的
3. 在这个假设下，最大化似然
4. 根据新的参数，重新更新对隐变量的猜测
5. 重复直到收敛

关键发现：
这个过程会单调地增加观测数据的似然函数
（通过Jensen不等式的巧妙应用）

而且通常会收敛到局部最优解

1
2
3
4
5
6
7
8

梯度下降：
直接对L(θ|X)求导
∇_θ log L(θ|X) = ?
（这很困难，因为有求和符号）

EM算法：
转化为更简单的问题
避免直接求导含有求和的似然函数

1
2
3
4
5
6
7
8

数值优化（如Newton法）：
需要计算目标函数的二阶导数
计算复杂且数值不稳定

EM算法：
只需要在完全数据上的期望
通常有闭形解或简单的更新规则
更稳定，计算更快

1
2
3
4
5
6
7
8
9

EM：频率学派
给出参数的点估计

贝叶斯：
给出参数的后验分布
但计算通常更复杂（需要MCMC）

EM的变种：Bayesian EM
在M步中加入先验信息

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

设有：
- X ∈ ℝ^{n×d}：观测数据（样本大小为n，维度为d）
- Z ∈ ℝ^{n×p}：隐变量（也是n个样本，维度为p）
- θ ∈ Θ：待估计的参数

完全数据：(X, Z)
不完全数据：X

概率模型：
P(X, Z | θ) = 联合分布
P(X | θ) = Σ_Z P(X, Z | θ) = 边际分布（观测似然）

目标：
最大化 log P(X | θ)（边际似然的对数）

困难：
这个最大化很难进行

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

Jensen不等式（对于凹函数f）：
f(E[X]) ≥ E[f(X)]

对于凹函数（如对数函数），
函数值的期望 ≤ 期望的函数值

在我们的问题中：
log是凹函数，所以：

log Σ_z P(X, Z=z | θ)
≥ Σ_z q(Z=z) · log [P(X, Z=z | θ) / q(Z=z)]

其中q(Z=z)是任意概率分布（Σ_z q(z) = 1）

这给了我们一个下界！

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

log P(X_i | θ) = log Σ_{z_i} P(X_i, Z_i=z_i | θ)
                = log Σ_{z_i} q_i(z_i) · [P(X_i, Z_i=z_i | θ) / q_i(z_i)]
                ≥ Σ_{z_i} q_i(z_i) · log [P(X_i, Z_i=z_i | θ) / q_i(z_i)]
                   （Jensen不等式）
                = Σ_{z_i} q_i(z_i) · [log P(X_i, Z_i=z_i | θ) - log q_i(z_i)]
                = Σ_{z_i} q_i(z_i) · log P(X_i, Z_i=z_i | θ)
                  - Σ_{z_i} q_i(z_i) · log q_i(z_i)
                = E_{q_i}[log P(X_i, Z_i | θ)] + H(q_i)

其中H(q_i)是q_i的熵（非负）

所以下界是：
log P(X_i | θ) ≥ E_{q_i}[log P(X_i, Z_i | θ)] + H(q_i)

为了最紧（最接近真实值），我们需要：
等号成立！

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

Jensen不等式中等号成立的条件：
X是常数（几乎处处）

在我们的问题中，意味着：
P(X_i, Z_i=z_i | θ) / q_i(z_i) = 常数

即：q_i(z_i) ∝ P(X_i, Z_i=z_i | θ)

正确的规范化后：
q_i(z_i) = P(Z_i=z_i | X_i, θ) （后验分布！）

这就是E步！

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

初始化：选择初始参数θ^(0)

重复直到收敛（t = 0, 1, 2, ...）：

  E步（Expectation Step）：
    计算隐变量的后验分布
    Q(θ, θ^(t)) = E_{Z|X,θ^(t)}[log P(X, Z | θ)]

    这是一个依赖θ的期望
    （在θ^(t)下计算，但作为θ的函数）

  M步（Maximization Step）：
    最大化Q函数：
    θ^(t+1) = arg max_θ Q(θ, θ^(t))

  检查收敛：
    如果 ||θ^(t+1) - θ^(t)|| < ε，停止
    否则，继续迭代

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

E步：
给定θ^(t)，计算期望的完全数据对数似然：

Q(θ | θ^(t)) = E_{Z|X,θ^(t)}[log P(X, Z | θ)]
              = Σ_z P(Z=z | X, θ^(t)) · log P(X, Z=z | θ)

对于连续隐变量：
Q(θ | θ^(t)) = ∫ P(Z | X, θ^(t)) · log P(X, Z | θ) dZ

解读：
在给定当前参数θ^(t)和观测X的条件下，
计算隐变量Z的后验分布
然后计算"如果参数是θ"时完全数据似然的期望

M步：
θ^(t+1) = arg max_θ Q(θ | θ^(t))

解读：
找到使Q函数最大的θ
（这个最大化通常比原问题简单得多）

关键性质：
log P(X | θ^(t+1)) ≥ log P(X | θ^(t))
观测似然单调增加！

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

为什么似然会单调增加？

记当前似然为：L_t = log P(X | θ^(t))

可以分解为：
L_t = E_{Z|X,θ^(t)}[log P(X, Z | θ^(t))]
      - E_{Z|X,θ^(t)}[log P(Z | X, θ^(t))]
    = Q(θ^(t) | θ^(t)) - H_t

其中第一项是完全数据似然的期望
第二项是后验分布的熵

进行M步得到θ^(t+1)：
L_{t+1} = Q(θ^(t+1) | θ^(t)) - H_{t+1}

现在，由于M步的定义：
Q(θ^(t+1) | θ^(t)) ≥ Q(θ^(t) | θ^(t))

但是H_{t+1}（在新参数下的熵）可能不同

关键证明步骤（使用KL散度）：
H_t - H_{t+1} ≤ Q(θ^(t+1) | θ^(t)) - Q(θ^(t) | θ^(t))

综合起来：
L_{t+1} - L_t = [Q(θ^(t+1) | θ^(t)) - Q(θ^(t) | θ^(t))]
                + [H_{t+1} - H_t]
              ≥ 0

所以L_{t+1} ≥ L_t，似然单调增加！

这保证了算法不会"倒退"
最终会收敛到某个局部最优解

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

假设数据来自K个高斯分布的混合：

X_i ~ Σ_{k=1}^K π_k · N(μ_k, Σ_k)

其中：
- π_k是第k个分量的混合权重（Σ_k π_k = 1）
- μ_k是第k个分量的均值
- Σ_k是第k个分量的协方差矩阵

隐变量：
Z_i ∈ {1, 2, ..., K}表示样本i属于哪个分量
（这是一个分类隐变量）

参数：
θ = {π_1, ..., π_K, μ_1, ..., μ_K, Σ_1, ..., Σ_K}

观测：
独立同分布的n个样本 {X_1, ..., X_n}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

给定当前参数θ^(t) = {π_k^(t), μ_k^(t), Σ_k^(t)}

计算每个样本属于每个分量的后验概率：

γ_{i,k}^(t) = P(Z_i = k | X_i, θ^(t))
             = π_k^(t) · N(X_i; μ_k^(t), Σ_k^(t))
               / Σ_j π_j^(t) · N(X_i; μ_j^(t), Σ_j^(t))

其中N(x; μ, Σ)是在x处的高斯概率密度

解读：
γ_{i,k}是样本i"属于"分量k的"软性"分配
不是硬的0/1（那就知道了），
而是一个概率（0到1之间）

例如：
样本i可能以70%的概率属于第1个分量（男生高斯）
以30%的概率属于第2个分量（女生高斯）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

基于E步计算的γ_{i,k}，更新参数

更新混合权重：
π_k^(t+1) = (1/n) · Σ_i γ_{i,k}^(t)

解读：
第k个分量的新权重 = 所有样本对该分量的后验概率平均

更新均值：
μ_k^(t+1) = Σ_i γ_{i,k}^(t) · X_i
            / Σ_i γ_{i,k}^(t)

解读：
第k个分量的新均值 = 所有样本的加权平均
权重是γ_{i,k}（样本属于该分量的概率）

更新协方差：
Σ_k^(t+1) = Σ_i γ_{i,k}^(t) · (X_i - μ_k^(t+1))(X_i - μ_k^(t+1))^T
            / Σ_i γ_{i,k}^(t)

解读：
第k个分量的新协方差 = 所有样本方差的加权平均
权重同样是γ_{i,k}

直觉：
如果某个样本很可能属于分量k（γ_{i,k}接近1），
它对该分量参数的估计贡献更大；
如果很不可能（γ_{i,k}接近0），
贡献很小

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

初始化：
K = 2（假设有两个高斯）
π_1^(0) = 0.5，π_2^(0) = 0.5
μ_1^(0) = 180，μ_2^(0) = 160
Σ_1^(0) = 100，Σ_2^(0) = 100

数据：
学生1身高175cm
学生2身高168cm
...（共100个学生）

迭代0→1：

E步计算γ：
对学生1（175cm）：
分量1的高斯密度 ∝ exp(-(175-180)²/200) = exp(-0.125) ≈ 0.882
分量2的高斯密度 ∝ exp(-(175-160)²/200) = exp(-1.125) ≈ 0.324

γ_{1,1} = 0.5 × 0.882 / (0.5×0.882 + 0.5×0.324) ≈ 0.731
γ_{1,2} ≈ 0.269

解读：学生1有73.1%的概率是男生，26.9%是女生

M步更新：
假设对所有100个学生计算了γ后
Σ_i γ_{i,1} = 65（平均来说，65个学生被分配给男生）
Σ_i γ_{i,2} = 35

π_1^(1) = 65/100 = 0.65
π_2^(1) = 35/100 = 0.35

（比例更新了，现在认为样本中男生占65%）

μ_1^(1) = Σ_i γ_{i,1}^(0) × X_i / 65
        （用加权平均计算新的男生身高均值）

迭代1→2, 2→3, ...继续，直到参数稳定

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

监测：
对数似然函数：
log L^(t) = Σ_i log[Σ_k π_k · N(X_i; μ_k, Σ_k)]

每次迭代后计算，应该单调增加

收敛条件：
||θ^(t+1) - θ^(t)|| < ε
（参数变化很小）

或者：
|log L^(t+1) - log L^(t)| < ε
（似然增加很小）

典型行为：
第1-5次迭代：快速增加
第5-50次迭代：逐渐减速
第50次迭代后：基本稳定

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

序列模型：时间序列 {X_1, X_2, ..., X_T}

假设存在隐状态序列 {Z_1, Z_2, ..., Z_T}
每个Z_t ∈ {1, 2, ..., K}

模型结构：

1. 初始分布：
   P(Z_1 = k) = π_k

2. 转移分布（Markov性）：
   P(Z_t = j | Z_{t-1} = i) = A_{ij}
   （从状态i转移到j的概率矩阵A）

3. 观测分布：
   P(X_t | Z_t = k) = B_k(X_t)
   （在隐状态k下观测X_t的概率）

完整模型：
P(X_{1:T}, Z_{1:T} | θ)
= P(Z_1) · ∏_{t=1}^{T-1} P(Z_{t+1}|Z_t) · ∏_t P(X_t|Z_t)

参数：
θ = {π, A, B}（初始分布、转移矩阵、观测分布）

观测数据：
只有X_{1:T}，不知道Z_{1:T}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

语音识别：
- X_t = 第t个音频特征（MFCC等）
- Z_t = 第t个音素（音素状态）
- A = 音素转移概率（某个音素后面接什么音素的概率）
- B_k = 音素k的声学模型

天气预测：
- X_t = 观测到的温度/湿度
- Z_t = 真实的天气状态（晴/阴/雨）
- A = 天气转移（今天晴明天还是晴的概率）
- B_k = 天气k下的观测分布（晴天时温度的分布）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

E步（Forward-Backward算法）：

计算前向变量α_t(i)：
P(X_{1:t}, Z_t=i | θ^(t))

计算后向变量β_t(i)：
P(X_{t+1:T} | Z_t=i, θ^(t))

组合得到：
γ_t(i) = P(Z_t=i | X_{1:T}, θ^(t))
       = α_t(i) · β_t(i) / P(X_{1:T})

ξ_t(i,j) = P(Z_t=i, Z_{t+1}=j | X_{1:T}, θ^(t))
         （两个连续时刻的联合后验）

M步：

π_k^(t+1) = γ_1(k)
           （初始分布 = 第1时刻是状态k的概率）

A_{ij}^(t+1) = Σ_t ξ_t(i,j) / Σ_t γ_t(i)
             （转移概率 = 转移次数 / 离开总次数）

B_k^(t+1) = 对于参数分布B_k的最大似然估计
           （通常基于状态k的加权观测）

直观性：
前向算法：从左到右，计算"到达这个点的概率"
后向算法：从右到左，计算"从这个点出发的概率"
两者结合：得到完整的后验分布

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

1. 完全缺失数据
   一些数据点完全没有观测

   例：问卷调查，某些受访者没有回答某些问题

   EM处理：
   E步：推断缺失值的分布
   M步：用推断的值进行参数更新

2. 部分隐变量
   某些变量导致了我们观测到的模式
   但变量本身不可观测

   例：GMM中的聚类标签，HMM中的隐状态

3. 潜在因子模型
   低维隐因子解释高维观测

   例：因子分析、ICA（独立成分分析）

1
2
3
4
5
6
7

当M步的最大值难以获得时，
只需要增加Q函数的值（不一定最大化）

θ^(t+1) = arg max_θ Q(θ | θ^(t)) 近似为
θ^(t+1) 使得 Q(θ^(t+1) | θ^(t)) ≥ Q(θ^(t) | θ^(t))

收敛性仍然保证（只是较慢）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

将M步分解为多个坐标方向

不是同时更新所有参数
而是逐个更新参数θ_j，其他参数固定

每次更新同样保证似然不减少

优点：
1. 计算简单
2. 有时会更快收敛
3. 易于并行化

缺点：
容易陷入更差的局部最优

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

当数据量很大时，完整EM可能很慢
（E步需要处理所有数据）

解决：在E步中随机抽样

每次迭代：
- 随机抽取一个小批量(mini-batch)数据
- 对这个批次进行E步
- 更新参数（M步）

优点：
计算快，可以处理大数据

缺点：
可能不收敛（只是渐近收敛）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

传统EM假设：
在M步后精确知道参数

某些情况下（如贝叶斯模型）：
参数本身有不确定性（有先验和后验）

解决：变分推断
用更简单的分布近似真实的后验分布
在简化的假设下进行EM

例如：平均场假设
假设所有隐变量条件独立给定观测
这样E步可以分解为独立的小问题

应用：
潜在狄利克雷分配（LDA）
变分自编码机（VAE）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

当参数空间很大时：
M步的最大化可能很困难

结合梯度下降：
E步：同常规
M步：用梯度下降而非解析解

θ^(t+1) = θ^(t) + η·∇_θ Q(θ^(t) | θ^(t))

优点：
处理大参数空间
灵活性高

缺点：
收敛可能更慢
需要选择学习率

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

当隐变量空间太大，无法精确计算Q函数时：
（例如，隐变量连续且维度高）

用蒙特卡罗方法近似期望：

E步近似：
从后验分布P(Z|X,θ^(t))中抽取样本 {Z^(s)}_{s=1}^S
近似 Q(θ|θ^(t)) ≈ (1/S)Σ_s log P(X, Z^(s) | θ)

M步：正常最大化

优点：
可以处理复杂的隐变量分布

缺点：
需要大量样本
可能方差很大

1
2
3
4
5
6
7
8

假设消费者在K个选择之间做决定
我们观测到的是最终的选择（比如购买哪个产品）
但不观测到：
- 考虑集（哪些产品被考虑）
- 决策过程（如何权衡）
- 内心偏好强度

可以用EM估计这些隐变量

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

场景：旅行方式选择
观测：
是否选择开车、乘坐公交、步行

模型：
U_i,j = V_i,j + ε_i,j
P(choice=j | X) = exp(V_j) / Σ_k exp(V_k)  （Logit）

但实际上，消费者有一个"考虑集" C_i
只在考虑集内的选项中选择

完整模型：
P(choice=j | C_i) = exp(V_j) / Σ_{k∈C_i} exp(V_k)

隐变量：C_i（每个消费者的考虑集）

EM估计：

E步：
给定观测的选择j，推断考虑集C_i的分布
（如果选择了j，那么j一定在C_i中）
（但其他选项是否在C_i中不确定）

M步：
基于推断的考虑集，更新参数

结果：
估计出每个选项被纳入考虑集的概率
比如：高铁作为出行方式，被50%的消费者考虑

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

经典案例（Heckman 1979）：
评估教育对工资的影响

观测数据问题：
我们只看到在工作的人的工资
不知道不工作的人会赚多少

不可观测的隐变量：
是否工作（工作vs.不工作）

这导致样本选择偏差
直接用OLS会有偏估计

传统解决：Heckman两步法

现代解决：EM算法

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

建立完整模型：

选择方程（决定是否观测）：
S_i = 1 if Z_i > 0，其中Z_i = W_i'γ + ν_i

结果方程（条件于选择）：
Y_i = X_i'β + ε_i （仅当S_i=1时观测）

隐变量：
对于S_i=0的样本，Y_i的值未知

EM估计：

E步：
对于缺失的Y_i（对应S_i=0），
推断其条件分布
E[Y_i | S_i=0, 数据] = ?

M步：
同时估计β和γ
考虑既观测到的Y（S_i=1）
也通过EM推断的Y（S_i=0）

优点（vs. Heckman）：
1. 不需要强的参数假设
2. 可以容易扩展到多个选择
3. 允许更灵活的相关结构

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

模型：
Y_i | Z_i=k ~ N(X_i'β_k, σ_k²)，有概率π_k

隐变量Z_i：样本i属于哪个子群体

观测：(Y_i, X_i)，但不知道Z_i

应用场景：

1. 消费行为的异质性
   有些消费者对价格很敏感，有些不敏感
   我们观测到购买决定，
   但不知道消费者属于哪一类

2. 企业生产率的异质性
   有些企业技术先进，有些落后
   我们观测到产出，但不知道企业类型

3. 交易数据分析
   有些交易是知情者，有些不是
   我们观测到价格冲击，但不知道交易类型

EM估计步骤：

E步：
计算每个样本属于每个子群体的概率
γ_{i,k} = P(Z_i=k | Y_i, X_i)
        = π_k · φ(Y_i; X_i'β_k, σ_k²)
          / Σ_j π_j · φ(Y_i; X_i'β_j, σ_j²)

M步：
对每个子群体k：
β_k^{new} = (Σ_i γ_{i,k} X_i X_i')^{-1} Σ_i γ_{i,k} X_i Y_i
（加权最小二乘，权重是γ_{i,k}）

σ_k^{2,new} = Σ_i γ_{i,k} (Y_i - X_i'β_k)² / Σ_i γ_{i,k}

π_k^{new} = (1/n) Σ_i γ_{i,k}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

传统问题：
Y_{it} = α_i + X_{it}'β + ε_{it}

如果直接加虚拟变量（dummies），参数太多（维度灾难）

如果用within转换消除α_i，会有偏估计（Nickell偏差，动态面板）

EM处理：

思想：
把α_i视为"缺失数据"的一种
用EM推断α_i

E步：
给定当前的β估计和ε的分布
计算α_i的后验分布E[α_i | 数据]

M步：
用α̂_i代替α_i，进行标准回归

优点：
1. 同时估计β和α_i
2. 获得α_i的完整分布，不仅是点估计
3. 可以进行推断（比如，预测新的α_j）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

策略1：随机初始化
随机选择参数起点
运行多次EM（多启动）
选择最好的结果（最高似然）

风险：可能陷入不同的局部最优

策略2：K-means初始化（对于GMM）
首先用K-means聚类
用聚类结果初始化：
- π_k = 聚类k的大小
- μ_k = 聚类k的中心
- Σ_k = 聚类k的样本方差

优点：通常更好

策略3：启发式初始化
利用领域知识
比如，对于语音识别，用已知的音素特征初始化

策略4：分层初始化
从简单模型开始（如1个高斯）
逐步增加复杂性（K=2,3,...）
用前一个K的结果初始化K+1

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

监测量1：对数似然
log L = log P(X | θ)

应该单调增加
如果下降，有bug

监测量2：参数变化
||θ^(t+1) - θ^(t)||

通常应该逐渐减小

监测量3：后验变化
max_i ||P(Z_i|X_i,θ^(t+1)) - P(Z_i|X_i,θ^(t))||

如果后验稳定，参数可能也稳定了

停止准则：
- 相对参数变化 < ε（如0.0001）
- 相对似然变化 < ε
- 最大迭代次数（如1000）

取决于应用的精度要求

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

直接计算会导致数值下溢（underflow）：
L = ∏_i P(X_i) 非常小，接近0

改进：用对数似然
log L = Σ_i log P(X_i)

但即使这样，仍有问题

对于混合模型：
log P(X_i) = log Σ_k π_k·P(X_i|θ_k)

问题：如果一个π_k很大，其他很小，
那么 Σ_k 被大项主导，导致数值不稳定

解决：log-sum-exp 技巧
log(a + b) = log a + log(1 + exp(log b - log a))
           = max(log a, log b) + log(1 + exp(-|log a - log b|))

代码：
def logsumexp(x):
    x_max = max(x)
    return x_max + log(sum(exp(x_i - x_max) for x_i in x))

这样避免了指数的下溢

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

计算协方差矩阵时：
Σ_k = (1/n_k) Σ_i γ_{i,k} (X_i - μ_k)(X_i - μ_k)^T

问题1：如果γ_{i,k}中有些非常小，数值精度丧失
解决：添加小的正则项 Σ_k + λ I

问题2：协方差矩阵可能不是正定的
解决：在计算前进行Cholesky分解
或者进行特征值分解，删除小的负特征值

问题3：当样本量小相对维数大时，协方差矩阵奇异
解决：
- 降维（PCA）
- 使用对角协方差矩阵（假设变量独立）
- 使用球形协方差（Σ_k = σ_k² I）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

一次EM迭代的计算：

对于GMM（K个分量，n个样本，d维数据）：

E步：
- 对每个样本和每个分量计算高斯密度：O(nKd)
- 规范化计算γ：O(nK)
合计：O(nKd)

M步：
- 计算均值：O(nKd)
- 计算协方差：O(nKd²)  （最昂贵）
合计：O(nKd²)

每次迭代：O(nKd²)
总共T次迭代：O(T·nKd²)

对于大规模问题（大n, 大d）：
可能很慢

优化策略：
1. 使用对角或球形协方差（避免d²因子）
2. 并行化（对样本并行）
3. 小批量EM（处理子集）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

需要存储：
- 数据X：O(nd)
- 参数θ：O(Kd²)（因为有K个协方差矩阵）
- 后验γ：O(nK)

总共：O(nd + Kd² + nK)

对于n很大：O(nd)主导
对于d很大：O(Kd²)主导

优化：
1. 不存储完整的X，流式处理
2. 对角协方差：O(Kd)
3. 不存储γ，边计算边使用

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

相同点：
都是最大化似然函数
都给出参数的点估计

不同点：
MLE：直接对观测似然P(X|θ)求导
    对于简单模型有效
    但对于有隐变量的模型很困难

EM：通过引入隐变量的期望
   将困难问题转化为简单问题
   迭代求解

何时用哪个：
- 简单的显式模型（直接MLE更快）
- 复杂的隐变量模型（EM更实用）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

MCMC（贝叶斯方法）：
从参数的后验分布中抽样
P(θ|X) ∝ P(X|θ)·P(θ)

EM（频率学派）：
点估计参数
θ̂ = arg max P(X|θ)

对比：

精度：
MCMC给出完整的后验分布
EM只给出点估计

不确定性：
MCMC：有后验方差，可以直接做推断
EM：需要另外估计方差（比如bootstrap）

计算速度：
EM：通常更快，迭代数从10-100
MCMC：需要更多迭代（1000-100000）

可扩展性：
EM：对大数据更友好
MCMC：对高维参数困难（维度诅咒）

应用：
需要点估计，数据很大 → EM
需要不确定性量化，参数少 → MCMC
两者的杂交（变分推断）逐渐流行

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

梯度下降：
直接计算∇ log P(X|θ)
θ^(t+1) = θ^(t) + α·∇_θ log P(X|θ)

问题（对于隐变量模型）：
∇ log P(X|θ) 的计算涉及 Σ_Z
求和和梯度不交换，很复杂

EM：
避免直接求导
通过E步简化问题
通常有闭形解（M步）

对比：

梯度下降优点：
- 通用，对任何可导函数有效
- 有现成的优化库
- 变种很多（SGD, Adam等）

EM优点：
- 对隐变量模型自然
- 通常有闭形解
- 保证似然不减少（稳定性好）
- 通常收敛更快（对隐变量模型）

何时用哪个：
- 简单的损失函数 → 梯度下降
- 隐变量模型 → EM
- 超大规模数据 → 随机梯度下降
- 结合使用 → EM做初始化，然后梯度下降精调

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

EM保证的是：单调增加似然
但不保证全局最优

可能陷入局部最优而不是全局最优

例子（GMM）：
初始化不当可能导致某个分量"坍缩"
（所有数据都被分配给一个分量）

这也是一个（局部）最优解
但不是我们想要的

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

1. 多启动（Multiple Restarts）
   从不同的初始点运行EM多次
   比较最终的似然
   选择似然最高的

   数量：K*K*... （对于K个分量，运行K²次）

2. 随机初始化 + 好的初始化
   结合两种方法

3. 初始化约束
   确保各个分量有足够的样本
   避免分量坍缩

4. 正则化
   添加先验（贝叶斯EM）
   增加参数的过度选择成本

5. 动态模型
   先从K=1开始，逐渐增加K
   用K-1的结果初始化K
   （这样找到的是全局最优的好近似）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

当模型复杂度（分量数K）增加时：
训练样本的似然会持续增加
但这可能不是真实的模型结构
而是过度拟合

例子：
真实数据来自2个高斯
但我们指定K=10
EM会努力找到10个分量的最优混合
即使多数分量是虚假的
似然会很高，但模型无法泛化

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

1. 模型选择标准

AIC（Akaike Information Criterion）：
AIC = -2·log L + 2·p
其中p是参数数量

BIC（Bayes Information Criterion）：
BIC = -2·log L + p·log(n)
其中n是样本数

选择使AIC或BIC最小的K

2. 交叉验证
将数据分成训练和验证集
对各种K值：
- 用训练集的EM估计参数
- 在验证集上评估对数似然
- 选择验证似然最高的K

3. 正则化
通过添加参数的先验惩罚复杂模型

4. 可视化诊断
检查学到的分量是否有意义
比如，在聚类结果上叠加真实的分量
看是否匹配

1
2
3
4
5

某些EM应用需要很多次迭代才能收敛
特别是当：
1. 隐变量高度不确定（后验分布很平缓）
2. 模型复杂，参数间相互依赖
3. 初始化不好

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

1. 更好的初始化
直接影响初始似然距离最优有多远
好初始化能显著减少迭代次数

2. 加速技巧

Aitken加速（Aitken Acceleration）：
猜测序列会线性收敛到极限
用Aitken公式外推
提前一步

数学：
如果θ^(t) → θ* 且收敛速率是r
Aitken加速给出：
θ̃^(t) ≈ θ* （更接近真实极限）
然后继续EM迭代

3. 共轭梯度法（Conjugate Gradient）
在M步中使用共轭梯度加速
（当M步本身也是优化问题时）

4. 一阶梯度EM
不进行完整M步的最大化
而只做一步梯度上升
（计算更快，但收敛不同）

5. 随机EM
处理小批量而非全部数据
通常更快达到"足够好"的解
（可能不会精确收敛）

1
2
3
4
5
6
7
8
9

在GMM中常见：
某个分量中的样本太少
该分量的方差趋向于0
似然趋向于无穷

例如：
数据中有一个离群点
EM可能会创建一个非常小方差的分量围绕该点
这使似然非常高，但模型无用

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

1. 方差下界
设置σ_k² ≥ σ_min
防止方差变得太小

2. 最小分量大小
要求每个分量至少有n_min个样本
否则合并或删除该分量

3. 正则化
添加先验使参数远离退化的值
比如，高斯先验：
P(σ_k²) ∝ exp(-(σ_k² - σ_prior)²)

4. 重新初始化
如果检测到退化（方差很小），
重新随机初始化该分量

5. 约束优化
在M步中加入约束
θ^(t+1) = arg max_θ Q(θ|θ^(t)) s.t. 约束条件
（如协方差必须满足某种条件）

1
2
3
4

例如，在社交网络分析中：
- 节点之间的连接是隐变量
- 网络太大，无法枚举所有可能的配置
- 精确的E步不可行

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

使用蒙特卡罗期望修正（MC-EM）：

E步修改为：
1. 从后验分布中抽样：Z^(s) ~ P(Z|X, θ^(t))
2. 近似期望：Q̂(θ|θ^(t)) ≈ (1/S) Σ_s log P(X, Z^(s) | θ)

这样避免了精确计算期望，只需要采样

M步：照常

优点：
可以处理复杂的隐变量空间

缺点：
需要大量样本S
如果S太小，估计有偏

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

VAE是神经网络的生成模型：

编码器（推断）：q_φ(z|x)
解码器（生成）：p_θ(x|z)

训练目标是最大化ELBO（下界）：
ELBO = E_{q_φ}[log p_θ(x|z)] - KL(q_φ(z|x) || p(z))

这与EM的思想完全对应：
- 第一项：类似M步（最大化似然）
- 第二项：类似E步的正则化（保持q接近先验）

联系：
EM是VAE的"概率"版本
VAE是EM的"深度学习"版本

1
2
3
4
5
6
7
8
9

GAN与EM的关系：
GAN：对抗性训练，没有显式的似然
EM：最大化似然

最近的研究试图结合两者：
- 学习隐变量的分布（类似E步）
- 学习数据分布（类似M步）

（这是一个活跃的研究方向）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

对于超大规模数据，可以使用分布式EM：

Mapper（映射）：
- 对分配到的数据子集进行E步
- 计算局部统计量（如Σ γ_{i,k}, Σ γ_{i,k} X_i等）

Shuffleand Sort（洗牌排序）：
- 按参数分组

Reducer（化简）：
- 合并局部统计量：全局 Σ γ_{i,k} = Σ_mapper Σ_local γ_{i,k}
- 进行M步计算全局参数更新

循环：
对每次迭代重复上述过程

优点：
可以处理PB级数据

缺点：
通信开销大
需要迭代同步（MapReduce 2-3个版本之间延迟）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

问题：有隐变量或缺失数据？

      ↓是

建模：参数化模型易否指定？

   ↙是          ↘否
  ↓             ↓
EM可用         考虑其他方法
           （可能问题本身不适合EM）

EM内：

能否闭形求解参数？

   ↙是                ↘否
  ↓                   ↓
标准EM              梯度EM/MCEM
（快）             （柔灵活）

数据规模？

   ↙小（<1M)        ↘大(>1M)
  ↓                 ↓
批EM               小批EM/随机EM
（精确）          （快速）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

适合EM：
☑ 有清晰的隐变量或缺失值
☑ 能够指定参数模型P(X,Z|θ)
☑ 参数空间中等大小（<10000）
☑ 隐变量空间不是太大
☑ 需要点估计（不关心完整的后验分布）

不适合EM：
☒ 隐变量空间极其复杂（>10维连续）
☒ 模型无法参数化（非参数方法更好）
☒ 需要参数的完整后验（用贝叶斯/MCMC）
☒ 高维离散选择（组合爆炸，用其他方法）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

□ 模型设定
  ☐ 隐变量清晰定义？
  ☐ 参数化分布选择合理？
  ☐ 边际化是否导致难以处理的形式？
  ☐ 是否有更简单的替代方法？

□ E步实现
  ☐ 后验分布有闭形吗？
  ☐ 否则，数值计算是否稳定？
  ☐ log-sum-exp是否用于避免下溢？
  ☐ 后验是否出现异常值？

□ M步实现
  ☐ 参数更新有闭形吗？
  ☐ 否则，数值优化收敛吗？
  ☐ 是否有参数约束（如协方差正定）？
  ☐ 是否需要正则化防止退化？

□ 算法设置
  ☐ 初始化策略清晰？
  ☐ 收敛标准合理？
  ☐ 最大迭代数足够？
  ☐ 是否进行多启动？

□ 诊断和验证
  ☐ 对数似然单调增加吗？
  ☐ 似然值合理吗？
  ☐ 参数估计合理吗？
  ☐ 后验分布形状合理吗？
  ☐ 敏感性分析：初始值的影响？

□ 结果报告
  ☐ 参数估计值和标准误？
  ☐ 对数似然最终值？
  ☐ 收敛用了多少迭代？
  ☐ 初始值敏感性？
  ☐ 模型选择（如K的选择）依据？

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

背景：
观测5000个基因在100个细胞样本中的表达水平
目标：识别不同的细胞类型

假设：
基因表达来自K个细胞类型的混合
每个细胞类型有不同的基因表达特征

模型：
X_{ij} ~ Σ_k π_k · N(μ_{jk}, σ_{jk}²)

其中i是细胞，j是基因

隐变量：每个细胞属于哪个类型

EM估计：

E步：
计算p_{ik} = P(细胞i属于类型k | X_i, θ)

M步：
μ_{jk} = Σ_i p_{ik} X_{ij} / Σ_i p_{ik}
σ_{jk}² = Σ_i p_{ik} (X_{ij} - μ_{jk})² / Σ_i p_{ik}
π_k = Σ_i p_{ik} / n

结果：
识别出4个细胞类型
计算每个细胞的"类型概率"
识别最能区分类型的基因

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

背景：
1000名工人的教育（年）和工资（万元）调查
缺失数据：
- 50名工人没有报告教育水平
- 200名工人没有报告工资

目标：估计教育对工资的影响，尽管数据缺失

模型：
Wage_i = α + β·Educ_i + ε_i

其中ε_i ~ N(0, σ²)

隐变量：缺失的Wage和Educ值

EM估计：

E步：
给定当前的(α, β, σ)
对于缺失的Educ：
E[Educ_i | Wage_i, 模型] = ...（贝叶斯推断）

对于缺失的Wage：
E[Wage_i | Educ_i, 模型] = α + β·Educ_i

M步：
β̂ = Cov(Wage, Educ) / Var(Educ)
（现在使用了推断的缺失值）

结果：
完整案例的OLS可能有偏
EM方法通过利用缺失机制给出更好的估计

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

优点：
✓ 对隐变量模型自然高效
✓ 保证似然不减少（稳定）
✓ M步通常有闭形解（计算快）
✓ 易于理解和实现
✓ 对部分缺失数据的处理很优雅

缺点：
✗ 可能陷入局部最优
✗ 收敛可能较慢（早期快，后期慢）
✗ 不直接给出参数的不确定性
✗ 只适合特定问题结构
✗ 需要能指定参数模型

何时是最好的选择：
- 中等规模数据（几千到百万）
- 明确的隐变量或缺失数据
- 参数化模型
- 需要高效的点估计

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

入门（1周）：
1. 理解混合高斯的直觉
2. 学习一次EM迭代的数学
3. 实现简单的GMM

进阶（2周）：
1. 理解Jensen不等式的作用
2. 学习隐马尔可夫模型
3. 实现HMM的Baum-Welch算法

高阶（3周）：
1. 学习EM的变种（Online EM, MCEM等）
2. 理解EM与变分推断的关系
3. 在实际数据上应用EM
4. 优化大规模实现

研究前沿（持续）：
1. EM在深度学习中的应用
2. 分布式EM
3. 自适应EM
4. Bayesian EM与不确定性量化

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

经典教材：
1. "Mixture Models and EM Algorithm"
   (Bishop, Pattern Recognition and Machine Learning)
2. "The EM Algorithm and Extensions"
   (McLachlan & Krishnan)

论文：
1. Dempster, Laird & Rubin (1977) - 原始EM论文
2. 各领域的应用论文

代码资源：
- scikit-learn中的GaussianMixture
- 各种R包（mixtools, mclust）
- PyMC/Stan中的变分方法

实践建议：
从简单的GMM开始
逐步理解更复杂的模型
重视可视化和诊断

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

1. EM算法的本质：
   通过迭代期望和最大化
   来处理隐变量或缺失数据

2. 关键性质：
   似然单调增加（由Jensen不等式保证）
   收敛到局部最优解

3. 两个步骤的含义：
   E步：给定参数，推断隐变量
   M步：给定推断的隐变量，更新参数

4. 实际应用的关键：
   正确的初始化
   仔细的数值实现
   充分的诊断检查

5. 何时使用：
   有隐变量和参数模型时
   是解决这类问题的优雅方法