★ FEATURED ARTICLE EM算法大概是机器学习里最容易被名字吓住的模型之一。我第一次接触“期望最大化”这几个字时满脑子都是抽象的期望和迭代符号直到后来把经典的“三硬币模型”从头到尾手算了一遍才真正理解它到底在忙什么。三硬币模型之所以被当成EM算法入门的第一课是因为它把所有核心要素都浓缩在一个足够简单的情景里有一组只能看到正反面的观测数据背后藏着一个看不见的选择过程我们需要估计三枚硬币各自的正面概率。这个过程看起来小但把隐变量、后验概率、E步、M步全都串起来了。这篇文章会从零开始拆解三硬币模型给出完整的数学推导、手算数值全过程、一份可以直接跑的NumPy实现以及我在实际使用中踩过的坑和边界条件。适合正在自学机器学习、准备面试或者跟EM算法纠缠许久但总觉得差点意思的读者。就算你现在对“期望最大化”只有模糊概念只要耐心跟着算一遍应该也能把这块硬骨头啃下来。1. 三硬币模型先明确我们到底在估什么1.1 游戏规则假设桌子上有三枚硬币分别叫A、B、C。每一次试验这样进行先抛硬币A。如果A正面朝上这一轮就选用硬币B如果A反面朝上就选用硬币C。再用选中的那枚硬币B或C抛一次记录结果是正面还是反面。重复上述过程 n 次最终我们只记录到一串“正面/反面”序列。问题来了我们看不到每一轮到底用了B还是C只能看到最终结果。请问三枚硬币各自的正面概率分别是多少这个设定就是三硬币模型的标准形态。A硬币的正面概率记为 pB硬币的正面概率记为 qC硬币的正面概率记为 r。我们要估计的参数就是 θ (p, q, r)。观测数据是 y₁, y₂, ..., yₙ每个 yᵢ 取值为 1正面或 0反面。真正决定每一轮用哪枚硬币的变量 zᵢ 看不到zᵢ 1 表示选了Bzᵢ 0 表示选了C这个 zᵢ 就是EM算法里常说的“隐变量”。为了后面手算方便固定一组观测数据n 101, 1, 0, 1, 0, 0, 1, 1, 1, 0也就是6次正面、4次反面。整篇文章会反复用这组数据做演示。1.2 把数学形式写出来给定参数 θ 后某一轮观测到正面 yᵢ 1 的概率是P(yᵢ 1 | θ) p · q (1 - p) · r意思很直白先由A决定走B还是C再由对应硬币给出正面概率。如果把 yᵢ 同时允许取0和1写成一个统一形式P(yᵢ | θ) p · q^yᵢ · (1-q)^(1-yᵢ) (1-p) · r^yᵢ · (1-r)^(1-yᵢ)这个式子并不复杂当 yᵢ 1 时后面变成 q 和 r当 yᵢ 0 时后面变成 (1-q) 和 (1-r)。两边的结构都一样只是被 p 加权后相加。隐变量 zᵢ 的分布也顺便写出来。在不知道观测结果的情况下P(zᵢ 1) pP(zᵢ 0) 1 - p。但一旦看到了 yᵢ我们就能根据贝叶斯公式算出P(zᵢ 1 | yᵢ, θ) p · q^yᵢ · (1-q)^(1-yᵢ) / P(yᵢ | θ)这个后验概率就是EM算法E步的核心对象。1.3 为什么一个“选硬币”也算一枚硬币初学者常问这里明明只需要两枚硬币B和C为什么非要加一枚A关键在于如果没有A每次直接抛B或C而且知道抛的是哪一枚那就没有隐变量参数的估计无非是数正面次数除以总次数小学生都能算。真正让问题变难的是“每一轮来自哪个分布不确定”这个不确定性就是由A引入的。所以A并不是可有可无的道具。它的正面概率 p 实际上控制着混合权重观测数据以概率 p 来自B分布以概率 1-p 来自C分布。这个结构在更复杂的混合模型里一模一样比如高斯混合模型 GMM 中每个样本以一定概率属于某个高斯分量那个概率就对应这里的 p。2. 直接做极大似然会卡在哪一步2.1 完整数据场景下可以闭式求解先做一个思想实验假设我们作弊每一轮用B还是C全都看得见。也就是说隐变量 zᵢ 是已知的。这时可以写出“完整数据”的对数似然。用 Z_B 表示选了B的总次数Z_C 表示选了C的总次数用 Z_BY 表示选了B且正面朝上的次数用 Z_CY 表示选了C且正面朝上的次数。那么log L(θ) Z_B · log p Z_C · log(1 - p) Σ zᵢ [ yᵢ log q (1-yᵢ) log(1-q) ] Σ (1-zᵢ) [ yᵢ log r (1-yᵢ) log(1-r) ]注意这个式子可以拆成三个互不纠缠的部分。分别对 p、q、r 求导并令其为零马上能得到闭式解p Z_B / n q Z_BY / Z_B r Z_CY / Z_C逻辑完全符合直觉p 是选B的比例q 是“选B时正面次数占选B总次数的比例”r 同理。如果隐变量可见这个估计过程没有任何悬念。2.2 缺失信息让对数似然变成“log(和)”现实是我们看不到 zᵢ。对单个样本只能把“走B”和“走C”两种情况都算出来再相加P(yᵢ | θ) p · B(yᵢ | q) (1-p) · B(yᵢ | r)其中 B(y|q) q^y · (1-q)^(1-y)。整个观测数据的对数似然是log L(θ) Σᵢ log [ p · B(yᵢ | q) (1-p) · B(yᵢ | r) ]麻烦就麻烦在 log 里面有一个加法结构。求导时分母会同时混入 p、q、r 三个参数导致没法像完整数据那样把 p、q、r 一个一个干净地解出来。2.3 拆不出闭式解的本质原因很多人第一次看到这里会问那就硬解呗最多是个非线性方程组用数值方法照样能做。确实能做但我们要理解EM算法存在的理由完整数据下的极大似然估计太简单了而观测数据的似然函数因为 log(和) 的结构变得复杂。一个自然的想法是能不能先用某种方式“猜”出缺失的隐变量或者更准确地说猜出隐变量的分布然后回到完整数据那个简单问题上EM算法的思路就是这么一个循环先用当前参数猜隐变量然后假装猜对了去做完整数据的极大似然更新参数。猜完再更新更新完再猜反复迭代。这里有一个数学上容易混淆的地方log(a b) 不等于 log(a) log(b)。很多初学推导的人会不自觉地把对数分配到和里面结果整条推导都错了。三硬币模型就是用来强化这个印象的最好例子。3. E步和M步的分工逻辑3.1 Q函数把“猜测”放进期望里EM算法的每次迭代分成两步分别是E步和M步。但先要定义Q函数。Q函数的形式是Q(θ | θᵗ) E_{Z | Y, θᵗ} [ log P(Y, Z | θ) ]翻译成人话给定当前参数 θᵗ 和观测数据 Y先算出隐变量 Z 的后验分布然后在这个后验分布下求“完整数据对数似然”的期望。为什么要这么做因为真实的后验分布我们不知道但可以拿当前参数近似。既然不知道每个样本到底来自B还是C那就给每个样本一个“软标签”比如有 0.36 的概率来自B、0.64 的概率来自C而不是直接硬性指定它来自哪个硬币。入门阶段不必急着证明EM算法单调收敛。你只需要理解直接最大化观测数据的对数似然很难但最大化这个 Q 函数在结构上跟完整数据极大似然很像每次迭代都能让真实对数似然不下降所以迭代下去会收敛。3.2 E步计算隐变量的后验概率E步要做的事非常具体对每个样本算出它来自B的后验概率。用上标 t 表示当前迭代参数μᵢ 表示第 i 个样本来自B的后验概率μᵢ P(zᵢ 1 | yᵢ, θᵗ) pᵗ · qᵗ^yᵢ · (1-qᵗ)^(1-yᵢ) / [ pᵗ · qᵗ^yᵢ · (1-qᵗ)^(1-yᵢ) (1-pᵗ) · rᵗ^yᵢ · (1-rᵗ)^(1-yᵢ) ]因为观测 yᵢ 只有0和1两种取值所以在三硬币模型里 μᵢ 只有两种值一种是“正面样本来自B的概率”一种是“反面样本来自B的概率”。E步就是把这 n 个 μᵢ 算出来。3.3 M步用加权统计更新参数拿到 μᵢ 之后M步要做的最大似然更新是p^(t1) (1/n) · Σᵢ μᵢ q^(t1) Σᵢ μᵢ · yᵢ / Σᵢ μᵢ r^(t1) Σᵢ (1-μᵢ) · yᵢ / Σᵢ (1-μᵢ)这三个公式从哪来其实把 Q 函数分别对 p、q、r 求导并令其为零就好。对 p 求导时只挑出含有 p 的项Σᵢ [ μᵢ log p (1-μᵢ) log(1-p) ]求导后得到 Σᵢ μᵢ / p - Σᵢ (1-μᵢ) / (1-p) 0整理就是 p 的更新公式。对 q 求导时只看含有 q 的项Σᵢ μᵢ [ yᵢ log q (1-yᵢ) log(1-q) ]求导后得到 q Σ μᵢ yᵢ / Σ μᵢ。分母是所有“软性来自B”的样本权重之和分子是其中“软性来自B且正面朝上”的加权次数。r 同理只是把 μᵢ 换成 1-μᵢ。这里必须强调一个容易写错的点q 的分母除的是 Σ μᵢ不是 n。因为 q 只跟来自B的样本有关来自C的样本哪怕权重很低也不能用来估计 q。同理r 只能用 Σ(1-μᵢ) 做分母。3.4 一个生活化的类比可以把 EM 想象成班主任在统计两个平行班的成绩。本来每个学生属于哪个班是清楚的但突然发现档案丢了只知道每个学生的分数不知道班级归属。E步拿着上一轮粗略了解到的班级平均分反推每个学生属于A班的概率。 M步用这个概率作为权重重新统计两个班的平均分。 第二轮再用更新的平均分反推归属反推完再更新平均分。几轮下来班级归属估计和平均分估计一起趋于稳定。三硬币模型的 E步 就是“反推每个观测来自B还是C的概率”M步就是“用软概率重新统计三枚硬币的正面率”。整套逻辑就是这个循环。4. 手算一轮EM数据、初值、迭代过程全走一遍4.1 观测数据与初始参数手算能极大提升对EM算法的理解。我们固定前面那组数据Y 1, 1, 0, 1, 0, 0, 1, 1, 1, 0统计一下正面 6 次反面 4 次。取一组经典初始值p⁰ 0.4, q⁰ 0.6, r⁰ 0.7这个初值和很多教材里的三硬币例子一致算起来比较顺。4.2 第一轮E步计算先算正面样本的后验概率。所有 yᵢ 1 的样本来自B的概率相同μ(正面) 0.4 × 0.6 / (0.4 × 0.6 0.6 × 0.7) 0.24 / 0.66 ≈ 0.3636解释一下分母0.24 是“走B且正面”的概率0.42 是“走C且正面”的概率加起来就是观测到正面的总概率 0.66。再算反面样本的后验概率μ(反面) 0.4 × 0.4 / (0.4 × 0.4 0.6 × 0.3) 0.16 / 0.34 ≈ 0.4706这里 0.16 是“走B且反面”的概率0.18 是“走C且反面”的概率分母是观测到反面的总概率 0.34。所以E步得到的 μ 序列是0.3636, 0.3636, 0.4706, 0.3636, 0.4706, 0.4706, 0.3636, 0.3636, 0.3636, 0.4706六个 0.3636四个 0.4706。4.3 第一轮M步计算先算所有 μᵢ 的总和Σμ 6 × 0.3636 4 × 0.4706 2.1816 1.8824 4.0640于是p¹ Σμ / n 4.0640 / 10 0.4064接下来算 q¹。分子是“加权后的正面计数”Σ μᵢ yᵢ 6 × 0.3636 2.1816因为反面样本 y0不会贡献到分子里。所以q¹ 2.1816 / 4.0640 ≈ 0.5368r¹ 的分子是“来自C的软权重中正面样本加权次数”Σ (1-μᵢ) yᵢ 6 × (1 - 0.3636) 6 × 0.6364 3.8184 Σ (1-μᵢ) 10 - 4.0640 5.9360 r¹ 3.8184 / 5.9360 ≈ 0.6432第一轮迭代结果θ¹ (0.4064, 0.5368, 0.6432)很多教材里最后展示的三硬币结果就是这个附近的结果比如 (0.4064, 0.5368, 0.6438)略差异来自四舍五入。4.4 第二轮验证为什么这个例子一次迭代就接近收敛正常EM算法是迭代很多次才收敛的但三硬币模型这里有个特殊现象第二轮E步算出来的 μ 会基本保持不变。验证一下。用 θ¹ (0.4064, 0.5368, 0.6432) 重新计算对正面样本 分子 0.4064 × 0.5368 ≈ 0.2182 分母 0.2182 0.5936 × 0.6432 ≈ 0.2182 0.3818 0.6000 μ(正面) ≈ 0.2182 / 0.6000 ≈ 0.3636对反面样本 分子 0.4064 × (1 - 0.5368) ≈ 0.4064 × 0.4632 ≈ 0.1883 分母 0.1883 0.5936 × (1 - 0.6432) ≈ 0.1883 0.2118 0.4001 μ(反面) ≈ 0.1883 / 0.4001 ≈ 0.4706μ 没变那么下一轮M步算出的参数也不会变迭代就在 θ¹ 附近停住。这不是巧合。三硬币模型的数据只有0/1两类所有正面样本共享同一个后验概率所有反面样本也共享同一个后验概率。当样本类别只有两种时EM在这种混合伯努利模型里有可能第一轮就跳到自洽解。所以如果你看到某些资料里三硬币模型只迭代了一两轮就出结果别觉得是偷工减料这是该模型的固有特征。但要注意这绝不代表EM算法通常一轮收敛。真正需要多轮迭代的是GMM这类每个样本观测值都不同的连续模型。三硬币模型只是适合用来算清公式它隐藏的这个特点反而容易让初学者产生错误预期。5. 初值敏感、不可辨识与常见误区5.1 换一组初值会发生什么EM算法对初值敏感这在三硬币模型里可以直观看到。把初值换成p⁰ 0.3, q⁰ 0.5, r⁰ 0.8第一轮E步μ(正面) 0.3 × 0.5 / (0.3 × 0.5 0.7 × 0.8) 0.15 / 0.71 ≈ 0.2113 μ(反面) 0.3 × 0.5 / (0.3 × 0.5 0.7 × 0.2) 0.15 / 0.29 ≈ 0.5172第一轮M步Σμ 6 × 0.2113 4 × 0.5172 1.2678 2.0688 3.3366 p¹ 3.3366 / 10 0.3337 Σμy 1.2678 q¹ 1.2678 / 3.3366 ≈ 0.3800 Σ(1-μ)y 6 × 0.7887 ≈ 4.7322 Σ(1-μ) 10 - 3.3366 6.6634 r¹ 4.7322 / 6.6634 ≈ 0.7102第二组的收敛结果是(0.3337, 0.3800, 0.7102)和第一组 (0.4064, 0.5368, 0.6432) 完全是不同的参数。这就是初值敏感性的直观体现。5.2 三硬币模型中的“自洽解”不止一个更值得思考的是这两个自洽解给出的观测分布几乎一样。第一组参数算出的正面总概率0.4064 × 0.5368 0.5936 × 0.6432 ≈ 0.2182 0.3818 0.6000第二组参数算出的正面总概率0.3337 × 0.3800 0.6663 × 0.7102 ≈ 0.1268 0.4732 0.6000两组参数不同但都对应“正面概率0.6反面概率0.4”的观测分布所以它们的对数似然基本一样。这说明三硬币模型存在不可辨识性从观测数据里你只能识别出最终的正面概率却无法唯一确定 (p, q, r) 三者分别是什么。原因也很简单观测数据只有一类概率分布正面概率这一个自由度却要估计三个参数。方程不够自然有无穷多组解。另外还存在分量标签交换问题把B和C互换把 p 换成 1-p观测分布完全不变。初学者看到算法跑出不同结果就以为实现出错其实很多时候是模型本身就不能唯一辨识。5.3 初值敏感与随机重启策略既然EM对初值敏感实际使用中几乎不会只跑一组初值。标准做法叫“随机重启”随机生成多组初始参数分别跑完整EM迭代记录每组收敛后的对数似然值最后选对数似然最高的一组作为最终结果。要注意的是EM的不动点未必是全局最优。它只保证每轮迭代后对数似然不下降但可能停在局部极大值或鞍点。多组初值同时启动就是为了尽量避开差的局部极值。另外初值不要取恰好是 0 或 1。比如 p⁰ 1 时算法会认为所有样本都来自B那么 r 将完全无法被估计分量C直接坍缩。这是混合模型里最常见的“分量坍缩”现象初值远离边界能有效降低发生概率。5.4 几个容易踩的认知坑第一个坑把E步理解成先给每个样本硬性分类再统计。如果是硬分类最后得到的就不是EM而是类似K-Means的硬分配聚类流程。EM用的是软概率 μᵢ每个样本对两类贡献都不为0这是关键区别。第二个坑以为EM保证得到全局最优解。它只保证“不下降”和“收敛到稳定点”这个稳定点可能是局部最优。实际项目里一定要配合多初值或先验知识。第三个坑在三硬币例子里看到一轮收敛就以为EM一般只需要迭代一两轮。三硬币模型由于观测数据只有两类后验概率只有两个取值才会出现快速固定的特殊情况。换成连续观测模型比如GMM通常要几十到上百轮。第四个坑忽略模型不可辨识性。三硬币模型本身不能唯一确定三个参数如果你拿着EM跑出的某组参数去解释业务含义要先想想是否有其他等价的参数组合。6. 用NumPy把三硬币EM跑起来6.1 完整代码实现理论说再多不如把代码跑一遍。用NumPy实现三硬币模型的EM只需要十几行import numpy as np Y np.array([1, 1, 0, 1, 0, 0, 1, 1, 1, 0], dtypefloat) p, q, r 0.4, 0.6, 0.7 max_iter 200 tol 1e-9 for t in range(max_iter): # E step分别计算来自B和来自C的概率 prob_b np.where(Y 1, p * q, p * (1 - q)) prob_c np.where(Y 1, (1 - p) * r, (1 - p) * (1 - r)) mu prob_b / (prob_b prob_c) # M step软权重统计 p_new mu.mean() q_new (mu * Y).sum() / mu.sum() r_new ((1 - mu) * Y).sum() / (1 - mu).sum() if max(abs(p - p_new), abs(q - q_new), abs(r - r_new)) tol: break p, q, r p_new, q_new, r_new print(fIterations: {t 1}) print(fp {p:.6f}, q {q:.6f}, r {r:.6f})运行后会出现我们手算的结果先是很少几次迭代就停下p 约 0.4064q 约 0.5368r 约 0.6432。你也可以把tol改小一点会发现迭代次数几乎不增加这正对应之前说的“一轮接近自洽”的三硬币模型特性。6.2 多初值自动筛选工程实践里推荐先跑多组随机初值再挑似然最好的结果。可以这样写rng np.random.default_rng(42) best_loglik -np.inf best_theta None for _ in range(50): p0, q0, r0 rng.uniform(0.01, 0.99, 3) p, q, r p0, q0, r0 for t in range(200): prob_b np.where(Y 1, p * q, p * (1 - q)) prob_c np.where(Y 1, (1 - p) * r, (1 - p) * (1 - r)) mu prob_b / (prob_b prob_c) p_new mu.mean() q_new (mu * Y).sum() / mu.sum() r_new ((1 - mu) * Y).sum() / (1 - mu).sum() if max(abs(p - p_new), abs(q - q_new), abs(r - r_new)) 1e-12: break p, q, r p_new, q_new, r_new # 计算观测数据的对数似然 prob_b np.where(Y 1, p * q, p * (1 - q)) prob_c np.where(Y 1, (1 - p) * r, (1 - p) * (1 - r)) loglik np.sum(np.log(prob_b prob_c)) if loglik best_loglik: best_loglik loglik best_theta (p, q, r) print(Best theta:, best_theta) print(Best log-likelihood:, best_loglik)你会发现不同初值收敛到的自洽解可能各不相同但被挑出来的最高对数似然值往往很接近。这就是模型本身不可辨识带来的结果并不一定是谁写错了。6.3 从三硬币到GMM和K-Means的迁移把三硬币模型跑通之后再去学GMM会轻松很多。GMM里的隐变量是“每个样本属于哪个高斯分量”E步算样本属于每个分量的后验概率M步按后验概率加权更新每个分量的均值、协方差和混合系数。你只要把这里 B(y|q) 的二项分布换成高斯密度函数几乎整个框架都能平移过去。K-Means也可以看成EM的一个极端特例E步把每个样本硬分配到一个簇相当于后验概率变成0/1M步更新簇中心。所以理解了软分配和硬分配的区别再去看K-Means和各种软聚类方法之间的演变脉络会清晰很多。在我实际学习过程中的体会是三硬币模型不要只看不练一定要自己完整算一遍E步和M步。算完之后再回头看EM算法的收敛性证明会发现那些看似绕的期望和Jensen不等式其实都是在说一件很朴素的事——用当前的猜测补全缺失信息然后重新统计再补全再统计。做项目时碰见缺失数据或者混合分布结构我也习惯先在类似的小模型上把迭代逻辑跑通再上完整数据能省掉很多调试时间。
阅读完成 · 觉得有帮助？
