第 8 章 EM 算法与隐变量模型
前几章讨论的求根、优化和积分，都是统计计算中的基础工具。本章介绍一个在不完全数据和隐变量模型中特别重要的算法：EM 算法。EM 是 Expectation–Maximization 的缩写，通常译为“期望最大化算法”。
在经济统计中，隐变量并不少见。例如，家庭收入调查中有些收入项缺失；消费者可能属于不同但未被直接观测的消费类型；企业信用风险中存在无法直接观测的风险等级；劳动力市场中也可能存在未观测的就业偏好或能力差异。模型中一旦出现这些未观测量，直接最大化观测数据似然往往比较困难。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\) 没有观测到，但可以在当前参数下计算它的条件期望，然后用这个期望替代完整数据对数似然中的未知部分。
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 也有局限：它可能收敛较慢，也可能停在局部极大点。因此实际使用时，仍然需要检查初始值、收敛曲线和最终模型解释是否合理。
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_i\)，若第 \(i\) 个观测属于第 \(k\) 个群体，则 \(z_i=k\)。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\) 个群体解释的概率”。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))
pi_k <- c(0.5, 0.5)
mu_k <- as.numeric(quantile(y, c(0.3, 0.7)))
sigma_k <- c(sd(y), sd(y))
loglik <- numeric(50)
for (iter in seq_along(loglik)) {
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[iter] <- 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])
))))
}
ord <- order(mu_k)
round(data.frame(
component = 1:2,
pi = pi_k[ord],
mu = mu_k[ord],
sigma = sigma_k[ord]
), 3)
#> component pi mu sigma
#> 1 1 0.592 1.673 0.242
#> 2 2 0.408 2.299 0.407
tail(loglik)
#> [1] -322.6 -322.6 -322.6 -322.5 -322.5 -322.5EM 的输出不是把每个家庭硬分到某一类，而是给出每个观测属于不同潜在群体的概率。这一点在经济统计解释中很重要：潜在群体是模型为了刻画异质性而引入的统计结构，不应简单等同于现实中的固定标签。
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 算法的实践建议
使用 EM 算法时，建议关注以下几点：
- 初始值。混合模型对初始值比较敏感，可以尝试多个初始值并比较最终对数似然。
- 收敛检查。记录每次迭代的观测数据对数似然，检查是否单调上升并趋于稳定。
- 局部最优。EM 保证局部改进，不保证找到全局最优。
- 解释约束。隐变量类别是模型结构，不一定对应现实中的天然类别。
- 缺失机制。缺失数据分析必须说明缺失机制假设，不能只报告算法结果。