第 8 章 EM 算法与隐变量模型
前几章讨论的求根、优化和积分，都是统计计算中的基础工具。本章介绍一个在不完全数据和隐变量模型中特别重要的算法：EM 算法。EM 是 Expectation–Maximization 的缩写，通常译为“期望最大化算法”。
在经济统计中，隐变量并不少见。例如，家庭收入调查中有些收入项缺失；消费者可能属于不同但未被直接观测的消费类型；企业信用风险中存在无法直接观测的风险等级；劳动力市场中也可能存在未观测的就业偏好或能力差异。模型中一旦出现这些未观测量，直接最大化观测数据似然往往比较困难。EM算法的作用，就是把一个难以直接优化的问题，转化为反复执行两个较容易的步骤。
本章的目标不是把 EM 算法讲成抽象优化理论，而是让读者掌握三件事：第一，为什么隐变量会让观测数据似然难以直接处理；第二，E 步和 M 步分别在计算什么；第三，如何通过对数似然轨迹、初始值和模型解释来判断 EM 结果是否可信。
8.1 不完全数据与隐变量
设观测数据为 \(y\)，隐变量或缺失数据为 \(z\)，模型参数为 \(\theta\)。如果 \(z\) 也能被观测到，那么完整数据似然为
\[ L_c(\theta)=p(y,z\mid\theta). \]
这里下标 \(c\) 表示 complete data。真实情况下我们只能看到 \(y\)，因此观测数据似然为
\[ L(\theta)=p(y\mid\theta)=\sum_z p(y,z\mid\theta), \]
若 \(z\) 是连续变量，则求和要换成积分。观测数据对数似然为
\[ \ell(\theta)=\log p(y\mid\theta) =\log\left\{\sum_z p(y,z\mid\theta)\right\}. \]
困难在于：对数外面包着求和或积分，直接优化常常不方便。相反，完整数据对数似然
\[ \ell_c(\theta)=\log p(y,z\mid\theta) \]
往往具有更简单的形式。EM 算法正是利用这一点：虽然 \(z\) 没有观测到，但可以在当前参数下计算它的条件期望，然后用这个期望替代完整数据对数似然中的未知部分。
一个典型困难来自“先求和再取对数”。例如两类混合模型的单个观测密度为
\[ p(y_i\mid\theta) = \pi_1 f_1(y_i\mid\theta_1)+\pi_2 f_2(y_i\mid\theta_2). \]
观测数据对数似然包含
\[ \sum_{i=1}^n \log\{\pi_1 f_1(y_i\mid\theta_1)+\pi_2 f_2(y_i\mid\theta_2)\}. \]
这里 \(\pi_k\) 是第 \(k\) 类的混合比例，\(f_k\) 是第 \(k\) 类的条件密度。由于对数函数包住了类别求和，直接对 \(\pi_k\)、\(\theta_k\) 求极大值并不容易。如果知道每个观测属于哪一类，问题会简单得多。EM算法正是把“未知类别”先用条件概率软填充，再在软填充后的完整数据上更新参数。
8.2 EM 算法
设当前参数为 \(\theta^{(t)}\)。EM 算法定义
\[ Q(\theta\mid\theta^{(t)}) = E_{\theta^{(t)}}\{\ell_c(\theta)\mid y\}. \]
这里的期望是在给定观测数据 \(y\) 和当前参数 \(\theta^{(t)}\) 的条件下，对隐变量 \(z\) 的条件分布\(p(z\mid y,\theta^{(t)})\) 取的。\(Q(\theta\mid\theta^{(t)})\) 可以理解为“在当前参数下，对完整数据对数似然的平均版本”。
EM 算法
- 选择初始参数 \(\theta^{(0)}\)，设置迭代次数 \(t=0\)。
- E 步：计算
\[ Q(\theta\mid\theta^{(t)}) = E_{\theta^{(t)}}\{\log p(y,z\mid\theta)\mid y\}. \]
- M 步：更新参数
\[ \theta^{(t+1)} = \arg\max_\theta Q(\theta\mid\theta^{(t)}). \]
- 检查对数似然或参数变化是否足够小。若未收敛，令 \(t\leftarrow t+1\)，返回第 2 步。
EM 算法的一个重要性质是：在一般条件下，每次迭代不会降低观测数据对数似然。也就是说，
\[ \ell(\theta^{(t+1)})\geq \ell(\theta^{(t)}). \]
这个性质使 EM 比很多直接优化方法更稳定。不过，EM 也有局限：它可能收敛较慢，也可能停在局部极大点。因此实际使用时，仍然需要检查初始值、收敛曲线和最终模型解释是否合理。
这个单调性可以用一个直观的不等式理解。对任意给定的隐变量分布 \(q(z)\)，有
\[ \log p(y\mid\theta) = \log\sum_z q(z)\frac{p(y,z\mid\theta)}{q(z)} \geq \sum_z q(z)\log\frac{p(y,z\mid\theta)}{q(z)}. \]
这里的不等式来自 Jensen 不等式，因为 \(\log(\cdot)\) 是凹函数。EM 在第 \(t\) 次迭代中取
\[ q(z)=p(z\mid y,\theta^{(t)}), \]
这会让下界在当前参数 \(\theta^{(t)}\) 处与真实对数似然相切。E 步构造这个下界，M 步最大化这个下界。因此，M 步虽然没有直接最大化原始观测数据对数似然，却能保证对数似然不下降。对本科生而言，记住这个图像就足够了：EM 每一步都在最大化一个当前点处贴住真实目标函数的、较容易处理的下界。
实际计算中常见的停止条件包括：
- 相邻两次对数似然变化很小：\(|\ell^{(t+1)}-\ell^{(t)}|<\epsilon\)；
- 相邻两次参数变化很小：\(\|\theta^{(t+1)}-\theta^{(t)}\|<\epsilon\)；
- 达到最大迭代次数；
- 对数似然出现下降或非有限数值，此时应检查代码、初始值和数值稳定性。
8.3 高斯混合模型
混合模型是理解 EM 算法最清楚的例子之一。设观测到的一维数据 \(y_1,\ldots,y_n\) 来自两个潜在群体。例如，在模拟的家庭收入数据中，一个群体可能对应工资收入较稳定的家庭，另一个群体可能对应经营性收入波动较大的家庭。我们只看到收入结果，却不知道每个家庭属于哪个群体。
两成分高斯混合模型可写为
\[ f(y_i\mid\theta) = \pi_1\phi(y_i;\mu_1,\sigma_1^2) +\pi_2\phi(y_i;\mu_2,\sigma_2^2), \]
其中 \(\phi(y;\mu,\sigma^2)\) 是均值为 \(\mu\)、方差为 \(\sigma^2\) 的正态密度，\(\pi_k\) 是第 \(k\) 个群体的比例，且 \(\pi_1+\pi_2=1\)。参数为
\[ \theta=(\pi_1,\pi_2,\mu_1,\mu_2,\sigma_1^2,\sigma_2^2). \]
引入隐变量 \(z_{ik}\)。若第 \(i\) 个观测属于第 \(k\) 个群体，则 \(z_{ik}=1\)，否则 \(z_{ik}=0\)。如果这些 \(z_{ik}\) 都能被观测到，完整数据对数似然为
\[ \ell_c(\theta) = \sum_{i=1}^n\sum_{k=1}^2 z_{ik} \left[ \log \pi_k+\log \phi(y_i;\mu_k,\sigma_k^2) \right]. \]
这里 \(z_{ik}\) 起到了“开关”的作用：第 \(i\) 个观测属于哪一类，就只贡献那一类的对数密度。由于\(z_{ik}\) 未被观测，E 步要计算它的条件期望，也就是责任度
\[ \tau_{ik} =P(z_i=k\mid y_i,\theta^{(t)}) = \frac{\pi_k^{(t)}\phi(y_i;\mu_k^{(t)},(\sigma_k^2)^{(t)})} {\sum_{j=1}^2\pi_j^{(t)}\phi(y_i;\mu_j^{(t)},(\sigma_j^2)^{(t)})}. \]
这里 \(\tau_{ik}\) 可以理解为“第 \(i\) 个观测由第 \(k\) 个群体解释的概率”。于是\(E(z_{ik}\mid y_i,\theta^{(t)})=\tau_{ik}\)。M 步使用这些概率作为权重：
\[ \pi_k^{(t+1)}=\frac{1}{n}\sum_{i=1}^n\tau_{ik}, \]
\[ \mu_k^{(t+1)} = \frac{\sum_{i=1}^n\tau_{ik}y_i}{\sum_{i=1}^n\tau_{ik}}, \]
\[ (\sigma_k^2)^{(t+1)} = \frac{\sum_{i=1}^n\tau_{ik}(y_i-\mu_k^{(t+1)})^2} {\sum_{i=1}^n\tau_{ik}}. \]
下面用模拟的对数收入数据演示 EM 算法。数据完全由程序生成，不对应真实调查来源。
set.seed(2026)
n <- 600
group <- rbinom(n, size = 1, prob = 0.35) + 1
y <- ifelse(group == 1,
rnorm(n, mean = 1.7, sd = 0.25),
rnorm(n, mean = 2.4, sd = 0.35))
em_normal_mixture2 <- function(y, mu_start, tol = 1e-8, max_iter = 200) {
n <- length(y)
pi_k <- c(0.5, 0.5)
mu_k <- mu_start
sigma_k <- c(sd(y), sd(y))
trace <- data.frame(iter = integer(), loglik = numeric())
for (iter in seq_len(max_iter)) {
dens <- cbind(pi_k[1] * dnorm(y, mu_k[1], sigma_k[1]),
pi_k[2] * dnorm(y, mu_k[2], sigma_k[2]))
tau <- dens / rowSums(dens)
nk <- colSums(tau)
pi_k <- nk / n
mu_k <- colSums(tau * y) / nk
sigma_k <- sqrt(c(
sum(tau[, 1] * (y - mu_k[1])^2) / nk[1],
sum(tau[, 2] * (y - mu_k[2])^2) / nk[2]
))
loglik <- sum(log(rowSums(cbind(
pi_k[1] * dnorm(y, mu_k[1], sigma_k[1]),
pi_k[2] * dnorm(y, mu_k[2], sigma_k[2])
))))
trace <- rbind(trace, data.frame(iter = iter, loglik = loglik))
if (iter > 1 &&
abs(trace$loglik[iter] - trace$loglik[iter - 1]) < tol) {
break
}
}
list(pi = pi_k, mu = mu_k, sigma = sigma_k, tau = tau, trace = trace)
}
fit_mix <- em_normal_mixture2(
y,
mu_start = as.numeric(quantile(y, c(0.3, 0.7)))
)
ord <- order(fit_mix$mu)
round(data.frame(
component = 1:2,
pi = fit_mix$pi[ord],
mu = fit_mix$mu[ord],
sigma = fit_mix$sigma[ord]
), 3)
#> component pi mu sigma
#> 1 1 0.663 1.693 0.253
#> 2 2 0.337 2.392 0.368
tail(fit_mix$trace)
#> iter loglik
#> 195 195 -322.2
#> 196 196 -322.2
#> 197 197 -322.2
#> 198 198 -322.2
#> 199 199 -322.2
#> 200 200 -322.2EM 的输出不是把每个家庭硬分到某一类，而是给出每个观测属于不同潜在群体的概率。这一点在经济统计解释中很重要：潜在群体是模型为了刻画异质性而引入的统计结构，不应简单等同于现实中的固定标签。
根据最终责任度，我们还可以查看分类不确定性。下面计算每个观测最大责任度的几个分位数。
max_tau <- apply(fit_mix$tau, 1, max)
quantile(max_tau, c(0.1, 0.25, 0.5, 0.75, 0.9))
#> 10% 25% 50% 75% 90%
#> 0.7051 0.8517 0.9460 0.9797 0.9952如果很多观测的最大责任度接近 0.5，说明两个潜在群体在数据上难以区分，模型解释应更加谨慎。如果最大责任度普遍接近 1，说明模型能够较清楚地区分两个群体。
8.4 缺失数据下的均值和方差估计
EM 算法也可以用于缺失数据。考虑一个简单例子：总体变量 \(Y\) 服从正态分布
\[ Y\sim N(\mu,\sigma^2), \]
但样本中有一部分 \(Y_i\) 缺失。若缺失机制可以近似看作与未观测值本身无关，则可以用观测数据估计\(\mu\) 和 \(\sigma^2\)。在这个例子中，隐变量就是缺失的 \(Y_i\)。
E 步需要计算缺失值的一阶和二阶条件期望。由于在当前参数下缺失值来自\(N(\mu^{(t)},(\sigma^2)^{(t)})\)，有
\[ E(Y_i\mid \text{missing},\theta^{(t)})=\mu^{(t)}, \]
\[ E(Y_i^2\mid \text{missing},\theta^{(t)}) =(\sigma^2)^{(t)}+\{\mu^{(t)}\}^2. \]
M 步再用这些期望替代缺失值，更新均值和方差。
set.seed(123)
n <- 300
y_full <- rnorm(n, mean = 6, sd = 1.2)
is_miss <- runif(n) < 0.25
y_obs <- y_full
y_obs[is_miss] <- NA
mu <- mean(y_obs, na.rm = TRUE)
sigma2 <- var(y_obs, na.rm = TRUE)
trace <- data.frame(iter = 0, mu = mu, sigma = sqrt(sigma2))
for (iter in 1:30) {
ey <- y_obs
ey2 <- y_obs^2
ey[is_miss] <- mu
ey2[is_miss] <- sigma2 + mu^2
mu_new <- mean(ey)
sigma2_new <- mean(ey2 - 2 * mu_new * ey + mu_new^2)
mu <- mu_new
sigma2 <- sigma2_new
trace <- rbind(trace,
data.frame(iter = iter, mu = mu, sigma = sqrt(sigma2)))
}
tail(trace)
#> iter mu sigma
#> 26 25 6.04 1.141
#> 27 26 6.04 1.141
#> 28 27 6.04 1.141
#> 29 28 6.04 1.141
#> 30 29 6.04 1.141
#> 31 30 6.04 1.141
c(complete_mean = mean(y_full),
observed_mean = mean(y_obs, na.rm = TRUE),
em_mean = mu)
#> complete_mean observed_mean em_mean
#> 6.041 6.040 6.040这个例子有意保持简单。真实调查中的缺失机制可能更复杂，例如高收入家庭更不愿意报告收入，或者企业在经营困难时更容易漏报某些指标。若缺失概率与未观测值本身有关，仅靠 EM 算法不能自动消除偏差，需要把缺失机制纳入模型或进行敏感性分析。
8.5 案例：上限编码收入数据的 EM 估计
收入和资产数据常常会进行上限编码（top-coding）：为了保护隐私，超过某个阈值的收入不报告精确值，只记录为“高于某个水平”。这种数据不是完全缺失，而是被右删失。直接把所有高收入者都当作阈值处理，会低估总体均值和方差。
为了说明 EM 如何处理这类问题，设对数收入
\[ Y=\log(\text{income}) \]
服从正态分布
\[ Y\sim N(\mu,\sigma^2). \]
若某个观测被上限编码，我们只知道 \(Y_i>c\)，其中 \(c\) 是上限编码阈值。在当前参数\(\mu^{(t)},\sigma^{(t)}\) 下，E 步需要计算
\[ E(Y_i\mid Y_i>c,\mu^{(t)},\sigma^{(t)}) = \mu^{(t)}+\sigma^{(t)} \frac{\phi(\alpha^{(t)})}{1-\Phi(\alpha^{(t)})}, \]
其中
\[ \alpha^{(t)}=\frac{c-\mu^{(t)}}{\sigma^{(t)}}. \]
还需要二阶矩
\[ E(Y_i^2\mid Y_i>c,\mu^{(t)},\sigma^{(t)}) = \{\mu^{(t)}\}^2+\{\sigma^{(t)}\}^2 + \sigma^{(t)}\{\mu^{(t)}+c\} \frac{\phi(\alpha^{(t)})}{1-\Phi(\alpha^{(t)})}. \]
这里 \(\phi(\cdot)\) 和 \(\Phi(\cdot)\) 分别为标准正态密度函数和分布函数。M 步把被删失观测的一阶矩和二阶矩代入样本均值与方差公式，更新 \(\mu\) 和 \(\sigma^2\)。
下面使用模拟数据演示。数据完全由程序生成，不对应真实收入调查。
set.seed(2028)
n <- 500
log_income_full <- rnorm(n, mean = 2.1, sd = 0.45)
top_code <- 2.55
is_top <- log_income_full > top_code
log_income_obs <- log_income_full
log_income_obs[is_top] <- top_code
# 初始值：把上限编码值暂时当作真实值
mu <- mean(log_income_obs)
sigma2 <- var(log_income_obs)
trace_top <- data.frame(iter = 0, mu = mu, sigma = sqrt(sigma2))
for (iter in 1:100) {
sigma <- sqrt(sigma2)
alpha <- (top_code - mu) / sigma
lambda <- dnorm(alpha) / pnorm(alpha, lower.tail = FALSE)
ey <- log_income_obs
ey2 <- log_income_obs^2
ey[is_top] <- mu + sigma * lambda
ey2[is_top] <- mu^2 + sigma2 + sigma * (mu + top_code) * lambda
mu_new <- mean(ey)
sigma2_new <- mean(ey2) - mu_new^2
trace_top <- rbind(
trace_top,
data.frame(iter = iter, mu = mu_new, sigma = sqrt(sigma2_new))
)
diff <- max(abs(c(mu_new - mu, sigma2_new - sigma2)))
mu <- mu_new
sigma2 <- sigma2_new
if (diff < 1e-8) break
}
c(topcoded_rate = mean(is_top),
true_mu = mean(log_income_full),
naive_mu = mean(log_income_obs),
em_mu = mu,
true_sigma = sd(log_income_full),
naive_sigma = sd(log_income_obs),
em_sigma = sqrt(sigma2))
#> topcoded_rate true_mu naive_mu em_mu true_sigma
#> 0.1640 2.1202 2.0758 2.1139 0.4516
#> naive_sigma em_sigma
#> 0.3805 0.4409
tail(trace_top)
#> iter mu sigma
#> 9 8 2.114 0.4409
#> 10 9 2.114 0.4409
#> 11 10 2.114 0.4409
#> 12 11 2.114 0.4409
#> 13 12 2.114 0.4409
#> 14 13 2.114 0.4409这个例子展示了 EM 的一个实际用途：它不是简单“补一个值”，而是在当前模型下计算被删失部分的条件期望，再用这些条件期望更新总体参数。需要强调的是，EM 的合理性依赖正态模型和上限编码机制。如果真实收入分布明显偏离对数正态，或者上限编码规则与其他变量有关，就需要更丰富的模型。
8.6 EM 算法的实践建议
使用 EM 算法时，建议关注以下几点：
- 初始值。混合模型对初始值比较敏感，可以尝试多个初始值并比较最终对数似然。
- 收敛检查。记录每次迭代的观测数据对数似然，检查是否单调上升并趋于稳定。
- 局部最优。EM 保证局部改进，不保证找到全局最优。
- 解释约束。隐变量类别是模型结构，不一定对应现实中的天然类别。
- 缺失机制。缺失数据分析必须说明缺失机制假设，不能只报告算法结果。
- 数值稳定性。混合模型中若某个分量方差被估得非常小，似然可能异常增大，应检查是否出现退化解。
- 多初始值比较。对复杂模型可以从多个初始值运行 EM，选择对数似然较高且解释合理的解。