当前位置：首页 > news >正文

【ML】EM算法：从三硬币到高斯混合模型的参数估计之旅

news 2026/6/3 18:34:02

1. 从三硬币问题理解EM算法的本质

我第一次接触EM算法是在研究生阶段，当时被这个看似简单实则精妙的思想震撼到了。让我们从一个经典的三硬币问题开始，逐步揭开EM算法的神秘面纱。

想象你面前有三枚特殊的硬币A、B、C，它们出现正面的概率分别是π、p、q。现在进行这样的实验：先掷硬币A，如果正面朝上就掷硬币B，反面朝上就掷硬币C；然后记录下第二次掷币的结果（1表示正面，0表示反面）。重复这个实验10次，得到观测序列：1,1,0,1,0,0,1,0,1,1。

这里的关键在于：我们只能看到最终的掷币结果，无法知道每次实验时到底是掷了B还是C。这种情况下，如何估计三个硬币的参数(π,p,q)？这就是著名的三硬币问题，也是理解EM算法最直观的案例。

在实际项目中，我经常遇到类似的情况。比如分析用户行为数据时，我们能看到用户最终点击了哪个广告，但不知道用户在选择过程中经历了哪些心理活动。这与三硬币问题何其相似！

2. EM算法的核心：E步与M步的舞蹈

2.1 E步：基于当前参数的"猜测"

E步全称Expectation Step，也就是期望步骤。在这个阶段，我们会基于当前的参数估计，对缺失的数据（隐变量）进行"最佳猜测"。

回到三硬币问题，假设我们已经有了参数(π,p,q)的初始值（可以随机设定），那么对于每个观测结果x_i，我们可以计算它来自硬币B的概率：

P(z_i=1|x_i) = π * p^x_i * (1-p)^(1-x_i) / [π * p^x_i * (1-p)^(1-x_i) + (1-π) * q^x_i * (1-q)^(1-x_i)]

这个计算本质上是在说："基于当前参数，我认为这个观测结果有百分之多少的可能性来自硬币B"。

2.2 M步：基于猜测的参数优化

M步全称Maximization Step，即最大化步骤。有了E步的"猜测"后，我们就可以更新参数了。具体来说：

更新π：所有实验中硬币A出现正面的平均概率
更新p：所有"可能"来自硬币B的结果中，正面的比例
更新q：所有"可能"来自硬币C的结果中，正面的比例

在实际编程实现时，我发现一个技巧：E步计算出的概率可以看作"软分配"，即一个观测结果可能同时属于硬币B和C，只是权重不同。这与K-means等"硬分配"算法形成鲜明对比。

3. 从三硬币到高斯混合模型

3.1 高斯混合模型简介

高斯混合模型(GMM)是EM算法最典型的应用场景之一。想象你面前有一组数据点，它们实际上来自多个不同的高斯分布，但你不知道每个点属于哪个分布。GMM的目标就是识别出这些隐藏的高斯分布。

我在图像分割项目中就使用过GMM。比如一张风景照中的像素点，天空、山脉、树木的像素值可能分别服从不同的高斯分布。通过GMM，我们可以自动发现这些潜在的分布。

3.2 GMM中的EM算法实现

在GMM中，EM算法的实现步骤如下：

初始化：随机设定各个高斯分布的参数(均值μ，协方差Σ)和混合系数π

E步：计算每个数据点属于各个高斯分布的概率（称为责任系数）

# 伪代码示例 def e_step(X, mus, sigmas, pis): responsibilities = [] for x in X: # 计算x对每个高斯分布的"贡献" probs = [pi * gaussian(x, mu, sigma) for pi, mu, sigma in zip(pis, mus, sigmas)] total = sum(probs) responsibilities.append([p/total for p in probs]) return responsibilities

M步：基于当前的责任系数重新估计参数

# 伪代码示例 def m_step(X, responsibilities): new_mus = [] new_sigmas = [] new_pis = [] for k in range(n_components): # 计算加权均值和协方差 resp_k = [r[k] for r in responsibilities] total_resp = sum(resp_k) new_mu = sum(r*x for r,x in zip(resp_k, X)) / total_resp new_sigma = sum(r*(x-new_mu)**2 for r,x in zip(resp_k, X)) / total_resp new_pi = total_resp / len(X) new_mus.append(new_mu) new_sigmas.append(new_sigma) new_pis.append(new_pi) return new_mus, new_sigmas, new_pis

重复E步和M步直到收敛

4. EM算法的数学原理与实现细节

4.1 为什么EM算法有效：Jensen不等式

EM算法的有效性建立在Jensen不等式的基础上。对于凹函数（如对数函数），有：

E[log(X)] ≤ log(E[X])

在EM算法的E步中，我们实际上是在构造对数似然函数的一个下界；在M步中，我们通过最大化这个下界来提升原始的对数似然函数。

我曾经在推导这个数学关系时遇到困难，直到画出了函数图像才恍然大悟。关键在于理解：通过引入隐变量的分布，我们能够"绕过"直接优化困难的原始问题。

4.2 实际实现中的技巧

在真实项目中应用EM算法时，有几个实用技巧：

初始化策略：随机初始化可能导致收敛到局部最优。我通常会尝试多次随机初始化，选择似然函数值最大的结果。
收敛判断：不要简单地比较参数变化，而是监控对数似然函数的变化。我常用的阈值是1e-6。
数值稳定性：在计算高斯密度时，加入小的正则项防止协方差矩阵奇异。
加速技巧：可以使用K-means的结果作为初始值，或者采用变分EM等改进算法。

5. EM算法的应用场景与局限性

5.1 典型应用场景

除了GMM，EM算法还广泛应用于：

隐马尔可夫模型(HMM)的参数学习
主题模型(如LDA)的参数估计
缺失数据填补
计算机视觉中的图像分割

我曾经在一个自然语言处理项目中使用EM算法来估计词性标注模型的参数，效果出奇地好。

5.2 局限性及解决方案

EM算法并非万能，存在以下局限：

收敛到局部最优：可以通过多次随机初始化缓解
收敛速度慢：特别是接近最优解时
需要知道隐变量的结构：即需要预先设定分布形式

在实践中，我通常会结合其他优化方法。比如先用EM算法快速找到一个不错的解，再用梯度下降法进行精细调整。

6. 从理论到实践：一个完整的Python实现

让我们用Python实现一个完整的GMM模型。这个实现包含了我在实际项目中学到的各种技巧：

import numpy as np from scipy.stats import multivariate_normal class GMM: def __init__(self, n_components, max_iter=100, tol=1e-6): self.n_components = n_components self.max_iter = max_iter self.tol = tol def fit(self, X): # 初始化参数 n_samples, n_features = X.shape self.weights_ = np.ones(self.n_components) / self.n_components self.means_ = X[np.random.choice(n_samples, self.n_components, replace=False)] self.covariances_ = [np.eye(n_features) for _ in range(self.n_components)] prev_log_likelihood = None for iteration in range(self.max_iter): # E步：计算责任系数 responsibilities = np.zeros((n_samples, self.n_components)) for k in range(self.n_components): responsibilities[:, k] = self.weights_[k] * \ multivariate_normal.pdf(X, mean=self.means_[k], cov=self.covariances_[k]) responsibilities /= responsibilities.sum(axis=1, keepdims=True) # M步：更新参数 Nk = responsibilities.sum(axis=0) self.weights_ = Nk / n_samples self.means_ = np.dot(responsibilities.T, X) / Nk[:, np.newaxis] for k in range(self.n_components): diff = X - self.means_[k] self.covariances_[k] = (responsibilities[:, k, np.newaxis] * diff).T @ diff / Nk[k] # 加入正则项确保数值稳定 self.covariances_[k] += 1e-6 * np.eye(n_features) # 检查收敛 log_likelihood = self._compute_log_likelihood(X) if prev_log_likelihood is not None and \ abs(log_likelihood - prev_log_likelihood) < self.tol: break prev_log_likelihood = log_likelihood def _compute_log_likelihood(self, X): likelihood = np.zeros((X.shape[0], self.n_components)) for k in range(self.n_components): likelihood[:, k] = self.weights_[k] * \ multivariate_normal.pdf(X, mean=self.means_[k], cov=self.covariances_[k]) return np.log(likelihood.sum(axis=1)).sum()

这个实现包含了几个关键点：