第一次被EM算法搞得睡不着是在处理一份用户行为数据的时候。手里几千条样本三分之一标签缺失直接扔掉太可惜用起来又不知道哪些样本可信。网上一搜所有教程都在说“EM算法是期望最大化”然后甩出Q函数、Jensen不等式、收敛证明。说实话第一遍看那些推导我是拒绝的公式每个字都认识连起来不知道在算什么。后来在GMM聚类、缺失值填充、隐马尔可夫模型里反复用它才慢慢总结出一套自己的理解方式。这篇文章不打算复刻教科书推导而是用三样东西把EM讲透一个允许你拿笔手算一遍的硬币例子、一段能在代码里跑通的GMM演示、一堆实际项目里踩出来的初值和退化坑。适合刚接触机器学习、知道极大似然但对隐变量发怵的同学也适合嘴上会说EM、真上手却总是报NaN的老朋友。1. EM算法到底卡在什么问题里1.1 极大似然估计的正常画风先看一个没有隐变量的场景。假设你在测试一枚新发行的纪念币投了100次正面62次让你估计正面概率几乎不需要思考62%。为什么因为结果服从二项分布。如果正面概率是p得到h次正面、t次反面的概率是L(p) p^h × (1-p)^t对p取对数再求导令导数为零d/dp [ h·log(p) t·log(1-p) ] h/p - t/(1-p) 0解得 p h/(ht)。这里没有隐变量因为所有信息都是已知的哪些投掷来自哪枚硬币清清楚楚。所以极大似然估计这条路走得很顺。但一旦出现“样本归属未知”事情就开始麻烦了。1.2 隐变量出现麻烦就来了现在有两枚硬币A和B每轮随机选一枚投10次记录正面次数但不告诉你本轮用的是哪一枚。你观测到的只是5轮实验的正面数5、7、4、8、6。想知道两枚硬币各自的正面概率θ_A和θ_B。这时候数据生成过程变成混合模型。假设每轮选A和B的概率都是0.5那么第i轮观测结果的概率是P(x_i | θ) 0.5 × θ_A^x_i (1-θ_A)^(10-x_i) 0.5 × θ_B^x_i (1-θ_B)^(10-x_i)把所有5轮放一起对数似然是L(θ) Σ_i log( 0.5 × θ_A^x_i (1-θ_A)^(10-x_i) 0.5 × θ_B^x_i (1-θ_B)^(10-x_i) )问题就出在这个结构上如果知道“本轮属于哪枚硬币”这个隐变量Z每个样本项就是单纯的“选A的0.5 A的贡献”取对数后可以拆开分别求导可现在log里面是一个求和求导之后没法把θ_A和θ_B干净地分离开闭式解根本写不出来。这就是EM存在的意义。不是EM想复杂而是隐变量逼着你换一种“先猜后改”的思路。1.3 核心思路先补全再更新循环迭代EM的做法很朴实既然不知道每个样本属于谁那就先用当前参数猜一个“软归属”把数据补全然后用补全后的数据重新估计参数再用新参数把归属猜得更准如此循环。打个比方。你想统计一个餐厅两位厨师的做菜咸淡但每道菜你没看到是谁做的只能尝出来。你先是假设这道菜有70%是师傅A做的于是把这道菜的盐度贡献按70%记到A头上、30%记到B头上。统计完所有菜之后重新估算两位师傅的平均盐度。再用新盐度去修正每道菜的权重再统计一次。旧权重不准不要紧因为只要方向正确每轮估计都会比上一轮更接近真实。这就是EME步是“按当前参数补全数据归属”M步是“按补全后的统计量重新估参数”。下面拿硬币问题把这一步一步算给你看。2. 拿笔算一遍两枚硬币的EM第一次迭代2.1 问题设定与初始参数5轮实验每轮投10次正面数分别是5、7、4、8、6。目标估计θ_A和θ_B。初始值随便给一组比如θ_A0.6θ_B0.5。为什么能随便给因为EM的迭代会自己修正初始误差只要别给得太离谱。但初值也不能完全不讲究后面第五节我会专门讲这个坑。2.2 E步算每条数据来自硬币A的后验概率第一轮E步要回答的问题在当前参数θ_A0.6、θ_B0.5下第一轮的5个正面更可能来自A还是B可能性各有多大如果来自A得到5个正面的概率是L_A 0.6^5 × 0.4^5 ≈ 7.96e-4如果来自BL_B 0.5^5 × 0.5^5 0.5^10 ≈ 9.77e-4那么本轮来自A的后验概率就是w_A L_A / (L_A L_B) ≈ 0.449注意这里不是“非A即B”的硬判定而是软分配。第1轮这10次投掷按0.449的权重算到A头上按0.551的权重算到B头上。后面每一轮都这么算结果列个表轮次正面数L_A(约)L_B(约)w_A(约)w_B(约)157.96e-49.77e-40.4490.551271.79e-39.77e-40.6470.353345.31e-49.77e-40.3520.648482.69e-39.77e-40.7330.267561.19e-39.77e-40.5500.450你会发现一个规律正面数越多它来自那枚“概率更高硬币”的权重就越大。第4轮有8个正面A的权重直接到0.733第3轮只有4个正面B的权重反而占优。这就是E步在做的“软补全”。2.3 M步用软分配结果重新估计参数现在我们把每轮的10次投掷按权重拆成两部分统计A的总正面数和总投掷次数。轮次A分得正面数A分得投掷数B分得正面数B分得投掷数15×0.449≈2.24510×0.449≈4.492.7555.5127×0.647≈4.5296.472.4713.5334×0.352≈1.4083.522.5926.4848×0.733≈5.8647.332.1362.6756×0.550≈3.3005.502.7004.50合计下来A分到大约17.35个正面、27.31次投掷B分到大约12.65个正面、22.69次投掷。于是新参数θ_A_new 17.35 / 27.31 ≈ 0.635θ_B_new 12.65 / 22.69 ≈ 0.558一轮迭代后参数从(0.6, 0.5)变成了(0.635, 0.558)。注意θ_B也上升了因为第一轮E步把正面比例较高的几轮多分了一些给A但B剩下的统计量同样推高了它的估计。这是合理的EM并没有要求两个参数一升一降。2.4 迭代到收敛代码也就十几行后面要做的事就是把E步和M步重复执行用新参数算新权重用新权重再更新参数。手工算太累直接用Python跑。import numpy as np x np.array([5, 7, 4, 8, 6]) # 每轮正面次数 n 10 # 每轮总投掷次数 theta_A, theta_B 0.6, 0.5 eps 1e-6 for it in range(200): # E步计算每条数据来自硬币A的后验概率 l_A theta_A**x * (1 - theta_A)**(n - x) l_B theta_B**x * (1 - theta_B)**(n - x) w_A l_A / (l_A l_B) w_B 1.0 - w_A # M步用软分配的权重累加统计量 A_heads np.sum(w_A * x) A_flips np.sum(w_A * n) B_heads np.sum(w_B * x) B_flips np.sum(w_B * n) new_A A_heads / A_flips new_B B_heads / B_flips delta max(abs(new_A - theta_A), abs(new_B - theta_B)) theta_A, theta_B new_A, new_B if it 10 or it % 20 0: print(fiter {it1:3d}: A{theta_A:.6f}, B{theta_B:.6f}) if delta eps: print(fconverged at iter {it1}) break运行之后你会看到前几轮参数变化很大后面变化越来越小直到满足停止条件。整个过程中每一轮的对数似然只增不减说明EM在稳步逼近某个极大值点。这段代码我建议你亲手跑一下。跑完你会有一种“哦原来EM这么朴素”的感觉。3. EM算法的数学框架Q函数和那套标准流程3.1 完整数据对数似然和观测数据对数似然硬币例子里如果每轮用的是哪枚硬币是已知的完整数据的对数似然写起来非常清爽log P(X, Z | θ) Σ_i log P(x_i, z_i | θ)其中z_i是第i轮选的硬币。可问题在于Z观测不到。我们只能对“所有可能的Z”求和log P(X | θ) log Σ_z P(X, z | θ)这个“log里面套求和”的结构让求导变成了灾难。EM回避硬啃原始目标转而构造一个更容易优化的替代函数。3.2 Q函数到底算什么EM每一轮都在维护一个叫作Q函数的家伙定义是Q(θ | θ^(t) ) Σ_i Σ_z P(z | x_i, θ^(t) ) × log P(x_i, z | θ)读法用上一轮参数θ^(t)算出隐变量z的后验分布然后把完整数据的对数似然在这个分布下求期望。初学者的第一个常见困惑就在这里为什么括号里有两个θ一个θ^(t)用来算后验权重相当于“上一轮的认知”另一个θ是这一轮要优化的变量相当于“待定的新答案”。说白了Q函数就是“把缺失数据用软权重填回去之后再去计算对数似然”。这和硬币例子的做法是逐项对应的。E步算出w_A、w_B那就是P(z | x_i, θ^(t))M步里用权重累加正次数那就是在对完整数据对数似然求期望。数学符号只是把那套话换了个形式。3.3 标准四步EM的标准流程可以浓缩成四步初始化参数 θ^(0)。E步用当前参数θ^(t)计算隐变量的后验分布P(z | x, θ^(t) )构造Q(θ | θ^(t) )。M步求解θ^(t1) argmax θ Q(θ | θ^(t) )。检查停止条件如果对数似然增量或参数变化小于阈值停止否则t加1回到第2步。在实际代码里第4步通常写成两个判断一个是增量阈值比如1e-6另一个是最大迭代次数防止某些场景下收敛太慢跑不出去。3.4 为什么这玩意儿真的收敛一个很多人没想明白的问题为什么m步提升Q函数原始对数似然也会跟着提升直觉是这样的。原始对数似然log P(X | θ)并不好优化但我们可以找到一个不那么紧的下界。EM每轮做的事情是在当前参数附近构造一个和原始对数似然在某点相切的下界函数然后对这个下界做最大化。因为下界永远不超过原始函数所以每次提升下界原始函数即使涨得不多也绝不会比上一轮低。用大白话说原始函数太高难爬你就在脚下垫一块新土把土堆到最高点然后人站上去离目标就更近一点。再垫一块再近一点。反复操作自然收敛到局部极大值附近。注意是局部极大值不是全局最优。这一点后面还会提到。4. GMM高斯混合模型EM最经典的落地场景4.1 从硬币换成高斯分布如果理解了硬币例子GMM几乎不需要额外动脑。硬币是二项分布GMM是高斯分布硬币里“哪一枚硬币”是隐变量GMM里“哪个高斯成分生成这个样本”是隐变量。其余全部相同。GMM假设数据由K个高斯分布按比例混合生成。每个分量有自己的均值μ_k、协方差Σ_k和混合权重π_k。你先以概率π_k选一个分量再从N(μ_k, Σ_k)里采样一个点。我们能观测到的只有采样点不知道它由哪个分量生成。这个设定特别适合做聚类和密度估计。比如客户分群你不想硬指定某个客户属于哪一类而是想知道他“有70%像高价值用户、30%像普通用户”。GMM就能给你这种软归属。4.2 E步在GMM里叫“响应度”GMM的E步核心是计算第i个样本属于第k个分量的后验概率记作r_ikr_ik π_k × N(x_i | μ_k, Σ_k) / Σ_j π_j × N(x_i | μ_j, Σ_j)分子是“在参数为μ_k、Σ_k的高斯分布下看到x_i的密度再乘上这个分量的先验权重”分母是对所有K个分量求和归一化。r_ik就是那个软权重对应硬币例子里的w_A和w_B。所有样本所有分量的r_ik合在一起就是一张N×K的软归属矩阵。每行加起来等于1。4.3 M步在GMM里就是加权平均M步看起来复杂实质只是带权重的平均。先算每个分量的有效样本量N_k Σ_i r_ik然后更新混合权重、均值、协方差π_k N_k / Nμ_k (1 / N_k) Σ_i r_ik × x_iΣ_k (1 / N_k) Σ_i r_ik × (x_i - μ_k)(x_i - μ_k)^T严格来说协方差更新里用的μ_k应该是本轮M步刚算出的新均值不是上一轮的旧均值。很多手写实现踩坑就在这先更新均值再用新均值更新协方差。每个公式的意思都非常直白把r_ik当权重对样本做加权平均。和硬币例子里数正面次数的逻辑一模一样。4.4 不用自己造轮子但要知道sklearn在干什么实际项目里直接用sklearn就行from sklearn.mixture import GaussianMixture from sklearn.preprocessing import StandardScaler scaler StandardScaler() X_scaled scaler.fit_transform(X) model GaussianMixture( n_components3, covariance_typefull, max_iter200, random_state42, init_paramskmeans ) model.fit(X_scaled) labels model.predict(X_scaled) proba model.predict_proba(X_scaled) # 软聚类概率 print(model.means_) print(model.weights_)这里有个细节GMM对数据尺度非常敏感跑之前先把X做标准化否则某个维度数值大协方差就容易出现奇异性。参数covariance_type有几种选择spherical、diag、tied、full。样本量不大时优先用diag而不是full因为full协方差要估计的参数多得多数据不够时很容易退化。关于init_paramskmeans这个非常值得说sklearn默认先用KMeans跑一遍把结果作为GMM的初值。这不是巧合而是EM对随机初值敏感KMeans能给出一个相对靠谱的起点能少跳很多坑。4.5 GMM和KMeans的关系KMeans可以看成GMM的一个极端特例。想象把每个高斯分量的协方差矩阵压缩到无限接近零再把所有混合权重设成相等。这时某个点只可能以接近1的权重属于距离最近的那个中心r_ik退化成0或1M步更新就变成了重算簇均值。所以KMeans是硬聚类GMM是软聚类。GMM比KMeans多做的两件事是给每个样本提供“属于各个簇的概率”以及用协方差刻画簇的形状。数据不是球形分布时GMM的表达能力会强很多。5. 实战中躲不开的初值、退化和监控问题5.1 初值不是随便给的我最早实现GMM时把K个分量的均值全部初始化成同一个向量。结果跑了二十轮所有样本都被分到同一个簇剩下的簇权重趋近于零似然不仅没涨还因为协方差奇异直接崩掉。教训是EM对初值极度敏感。全零、全相等、或者具有完全对称性的初值会让算法停在鞍点或者原地不动。标准做法有两种。一是多组随机初值跑完后选择对数似然最高的那组结果二是先用KMeans聚类把聚类的中心和协方差作为初值这条路在GMM里最常见。5.2 协方差矩阵奇异或退化GMM里最折磨人的报错就是协方差矩阵非正定。出现场景通常很典型某个分量几乎没有样本被分配到或者样本维度比每个簇的有效样本量还高协方差估计出来行列式接近零于是对数似然直接冲到无穷大。应对手段按优先级排数据标准化避免量纲悬殊导致数值问题。协方差选diag或spherical减少待估参数数量。给协方差对角加一个很小的正则项εI相当于强制它保持正定。sklearn里的reg_covar干的就是这件事。每次M步检查N_k如果发现某个分量的有效样本量过小直接把这个分量重置或丢掉防止它拖垮整体。5.3 概率下溢与log-sum-expGMM的E步要计算高斯密度N(x_i | μ_k, Σ_k)。维度一高密度值动辄是1e-50这种级别计算机直接把它当成0处理除法分母变0权重全部NaN。解决方式是在log空间里算然后用减去最大值的技巧归一化# 先计算每个分量的 log 概率 logp [] # 每个元素的长度是 K for k in range(K): logp_k log(pi[k]) - 0.5 * ( d * log(2*pi) log(det(Sigma[k])) (x - mu[k]) inv(Sigma[k]) (x - mu[k]) ) logp.append(logp_k) logp np.array(logp) logp - np.max(logp) # 减最大值防止 exp 溢出 p np.exp(logp) w p / p.sum()这段代码的思路是exp(big)会溢出但exp(big - max)不会。先对log概率做平移再exp归一化结果是一样的。5.4 迭代过程必须监控对数似然我有个雷打不动的习惯每次EM迭代都打印当前对数似然。正常情况下这个值应该单调不减。哪一步突然掉下来几乎可以断定是代码Bug而不是EM算法本身的问题。常见的Bug来源有三个E步的权重忘记归一化M步更新顺序错了比如用旧均值更新协方差或者数据里有NaN导致某些分量退化。先查这三处90%的问题都能定位。停止条件我一般用“参数变化 对数似然增量”双保险同时设置max_iter上限。只设固定迭代次数不行有的数据跑30轮就收敛有的要几百轮只看某一项指标也不够偶尔会出现参数没什么变化但似然还在缓慢爬升的情况。6. EM算法在哪些地方悄悄出现6.1 缺失数据填补EM天生适合缺失数据问题。比如一份用户问卷一部分人没填年龄一部分人没填收入。直接把缺失样本丢了吧浪费填个均值吧又太粗暴。EM的做法是根据现有观测数据估计缺失值的后验分布然后用这个分布去补全再更新模型参数。实际项目里我经常把EM当作数据清洗前的预处理先跑一个简单的混合模型把缺失字段的期望值算出来再进入后面的建模流程。效果比均值填充和删除样本要稳。6.2 HMM训练里的Baum-Welch算法隐马尔可夫模型中你只能观测到输出序列看不见背后的状态序列。这个状态序列就是隐变量。Baum-Welch算法本质就是EME步用前向-后向算法计算每个时刻处于各状态的后验概率M步更新转移概率、发射概率和初始概率。这也解释了为什么语音识别、NLP这类领域里EM是绕不开的基础工具。6.3 和变分推断的关系EM有个硬条件E步要能精确计算后验P(z | x, θ)但很多模型里这个后验根本算不出来。这时就用变分推断找一个简单的分布q(z)去近似它把EM的“精确E步”改成“近似E步”这就是变分EM。你后面如果学变分自编码器VAE会发现它的核心也是这套思路。先理解EM再理解VAE会顺畅很多因为VAE可以理解成“EM在深度模型里的近似版本”。6.4 它的边界在哪EM不是万能的。它只能保证找到局部最优无法保证全局最优每次迭代要扫全量数据样本量大的时候速度会很难看对高维隐变量模型精确后验算不出来必须转向变分推断或随机梯度方法。所以拿到一个实际问题先问自己我面对的到底是不是隐变量问题如果是EM大概率是个好起点如果不是硬套EM只会绕远路。文章写到这如果你真的把那组硬币数据跑了一遍应该已经对EM有了肌肉记忆。最后分享一个小习惯每次跑EM我会把对数似然曲线画出来。看着它单调爬上去再变平心里就知道实现大概没问题哪一步突然掉下来那一定不是EM的问题是你实现的问题。先查E步的除法再查M步的权重有没有归一化八成能找到Bug。这个排查顺序救过我很多次。