EM算法大概是机器学习里最容易被名字吓住的模型之一。我第一次接触“期望最大化”这几个字时，满脑子都是抽象的期望和迭代符号，直到后来把经典的“三硬币模型”从头到尾手算了一遍，才真正理解它到底在忙什么。三硬币模型之所以被当成EM算法入门的第一课，是因为它把所有核心要素都浓缩在一个足够简单的情景里：有一组只能看到正反面的观测数据，背后藏着一个看不见的选择过程，我们需要估计三枚硬币各自的正面概率。这个过程看起来小，但把隐变量、后验概率、E步、M步全都串起来了。
这篇文章会从零开始拆解三硬币模型，给出完整的数学推导、手算数值全过程、一份可以直接跑的NumPy实现，以及我在实际使用中踩过的坑和边界条件。适合正在自学机器学习、准备面试，或者跟EM算法纠缠许久但总觉得差点意思的读者。就算你现在对“期望最大化”只有模糊概念，只要耐心跟着算一遍，应该也能把这块硬骨头啃下来。
1. 三硬币模型：先明确我们到底在估什么
1.1 游戏规则
假设桌子上有三枚硬币，分别叫A、B、C。每一次试验这样进行：
- 先抛硬币A。如果A正面朝上，这一轮就选用硬币B；如果A反面朝上，就选用硬币C。
- 再用选中的那枚硬币（B或C）抛一次，记录结果是正面还是反面。
- 重复上述过程 n 次，最终我们只记录到一串“正面/反面”序列。
问题来了：我们看不到每一轮到底用了B还是C，只能看到最终结果。请问三枚硬币各自的正面概率分别是多少？
这个设定就是三硬币模型的标准形态。A硬币的正面概率记为 p，B硬币的正面概率记为 q，C硬币的正面概率记为 r。我们要估计的参数就是 θ = (p, q, r)。观测数据是 y₁, y₂, ..., yₙ，每个 yᵢ 取值为 1（正面）或 0（反面）。真正决定每一轮用哪枚硬币的变量 zᵢ 看不到，zᵢ = 1 表示选了B，zᵢ = 0 表示选了C，这个 zᵢ 就是EM算法里常说的“隐变量”。
为了后面手算方便，固定一组观测数据，n = 10：
1, 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) = p，P(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^(t+1) = (1/n) · Σᵢ μᵢ q^(t+1) = Σᵢ μᵢ · yᵢ / Σᵢ μᵢ r^(t+1) = Σᵢ (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因为反面样本 y=0，不会贡献到分子里。所以：
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], dtype=float) 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(f"Iterations: {t + 1}") print(f"p = {p:.6f}, q = {q:.6f}, r = {r:.6f}")运行后会出现我们手算的结果：先是很少几次迭代就停下，p 约 0.4064，q 约 0.5368，r 约 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/1；M步更新簇中心。所以理解了软分配和硬分配的区别，再去看K-Means和各种软聚类方法之间的演变，脉络会清晰很多。
在我实际学习过程中的体会是：三硬币模型不要只看不练，一定要自己完整算一遍E步和M步。算完之后再回头看EM算法的收敛性证明，会发现那些看似绕的期望和Jensen不等式，其实都是在说一件很朴素的事——用当前的猜测补全缺失信息，然后重新统计，再补全，再统计。做项目时碰见缺失数据或者混合分布结构，我也习惯先在类似的小模型上把迭代逻辑跑通，再上完整数据，能省掉很多调试时间。
