我最早把EM算法用在项目里是因为要做一批用户行为向量的聚类。数据是二值化的每个用户对若干品类有买/没买两种状态直接套K-means效果惨不忍睹因为欧氏距离对0/1向量太不友好了簇中心也落在“半买半不买”这种毫无意义的位置上。后来换成混合伯努利模型用EM算法求解聚类结果才真正变得可解释每个簇能给出一套“该品类出现概率”的画像。这篇文章就把这套方法从原理、推导到代码完整地过一遍适合刚接触EM算法的人也适合想快速上手混合伯努利模型做二值数据聚类的同学。1. 混合伯努利模型到底在拟合什么1.1 哪些场景会在意“0/1向量的簇”先想一个问题什么样的问题会产出一堆0/1向量并且你还想给它们分簇最常见的是文本的词袋表示。把一篇文章表示成一个维度为词典大小的向量某一维为1表示这个词在这篇文章里出现过为0表示没出现过。这时候“距离”如果拿欧氏距离来算两篇都提到很多常见词的文章可能距离很近但主题未必一样而提到同样生僻词的两篇文章欧氏距离也很近这个信号反而强烈。这类数据的问题本质不是“空间距离”而是“模式相似度”。用户行为数据也类似。电商场景下每个用户是否购买过某个品类可以编码成一个M维的0/1向量内容平台上每个用户是否阅读过某类内容也可以这样编码。想做用户分群时表面上是聚类实际要的是“把行为模式相同的人归到一起”而不是“把坐标距离近的人归到一起”。图像的二值特征也一样。手写数字的像素二值化、医学影像里的病灶区域二值掩码在降维成特征后经常用这种分布去拟合。这些场景有一个共同点数据本身的生成过程可能是一个“先选组再按组的概率生成0/1”的过程。这就是混合伯努利模型的假设来源。它假设数据集里存在K个隐性的组每个组对应一组伯努利概率参数某条样本先以一定的概率被分配到某个组再根据该组的伯努利参数逐维掷硬币生成0/1。这句话翻译成代码逻辑就是如果已知每条样本属于哪个簇聚类就简单了直接按组统计1的比例就行。问题是我们不知道而且“先选组”这一步也是随机的。于是问题天然带有隐变量。1.2 模型公式和参数含义用更规范的话说混合伯努利模型对单条样本x的概率定义为p(x|θ) Σ_{k1}^{K} π_k · Π_{d1}^{D} μ_{kd}^{x_d} (1-μ_{kd})^{1-x_d}这里的参数包括π_k第k个混合分量的权重满足π_k ≥ 0Σ π_k 1相当于“随机抽取一个样本它属于第k个簇的先验概率”。μ_{kd}第k个簇在第d个维度上取值为1的概率0 ≤ μ_{kd} ≤ 1这就是每个簇的“质心”只不过质心变成了一个概率向量。K簇的数量需要提前指定。D特征维度数。N样本总数。这个式子里Π那一串是一个多维伯努利分布的完整概率它假设在给定簇k的条件下D个维度之间相互独立。这个独立性假设可能不完全符合现实但对聚类任务来说通常足够用因为它真正关心的是“哪些维度组合模式的概率高”。整个数据集的对数似然是所有样本对数似然的加和L(θ) Σ_{n1}^{N} ln( Σ_{k1}^{K} π_k · Π_{d1}^{D} μ_{kd}^{x_{nd}} (1-μ_{kd})^{1-x_{nd}} )如果你直接对这个式子求导并令其等于零会发现问题对数里面套着对K个分量的求和求导之后每一项的分母里都有同一个求和项K个分量的参数互相纠缠无法得到闭式解。这正是EM算法登场的理由。2. EM算法的核心思路把缺失标签先补上再估计2.1 直接最大化对数似然的困难在哪很多讲EM的资料一上来就列公式容易把第一次接触的人吓跑。其实EM解决这个问题的思路非常朴素就一个词缺啥补啥。我们给每条样本引入一个隐变量z_n它表示这条样本真正所属的簇取值1到K。如果把z_n也当成已知数据那么“完整数据”的对数似然就非常好写L_c(θ) Σ_{n1}^{N} Σ_{k1}^{K} [z_{nk} · ( ln π_k Σ_{d1}^{D} ( x_d ln μ_{kd} (1-x_d) ln(1-μ_{kd}) ) )]其中z_{nk}是一个one-hot向量只有样本n真正属于的那个簇对应的位置为1其余为0。这个式子里没有“对K求和套在对数里面”的问题对π、μ求导都能得到闭式解。问题就是z_n是未知的。EM的思路是既然未知那我就猜它的期望用期望代替真实值做参数估计。这就是E步然后用这个期望值去最大化完整数据似然得到新的参数这就是M步。如此反复迭代直到对数似然不再上升。2.2 E步的完整推导E步要算的是后验概率在给定当前参数θ^{(t)}和观测数据x_n的条件下样本n属于簇k的概率。直接用贝叶斯公式γ(z_{nk}) p(z_{nk}1 | x_n, θ^{(t)}) π_k · p(x_n | μ_k) / Σ_{j1}^{K} π_j · p(x_n | μ_j)这里的p(x_n | μ_k)就是第k个伯努利分量对x_n的概率也就是前面Π那一串。这个γ值就是“软标签”。如果某个样本在第4个分量上的后验概率是0.85就说明它大概率属于第4个簇但仍有15%的可能性属于其他簇。这和K-means那种硬分配有本质区别K-means直接一步到位把样本钉死在某一个簇里EM则保留不确定性。实现时有个细节多个伯努利概率相乘连乘几百个接近0或接近1的值之后很容易变成0、inf或者NaN。所以实际工程里几乎从来不在原始概率空间算γ而是把分子分母都搬到对数空间先算log p再做一次log-sum-exp归一化最后exp回去。后面代码部分会展示这个技巧。2.3 M步的参数更新有了γ我们就可以把“完整数据”里的z_{nk}替换成γ(z_{nk})然后对参数求导。先看π。约束条件是Σπ_k 1用拉格朗日乘数法构造拉格朗日函数Λ Σ_n Σ_k γ_{nk} ln π_k λ(Σ_k π_k - 1)对π_k求导并令其为零得到π_k^{new} (1/N) Σ_{n1}^{N} γ_{nk}意思是新权重等于所有样本分配给簇k的平均“软占比”。这很直觉如果所有样本分给簇k的后验概率加起来是180总共600条样本那π_k就是0.3。再看μ_{kd}。完整数据对数似然里和μ_{kd}有关的项是Σ_n γ_{nk} [ x_{nd} ln μ_{kd} (1-x_{nd}) ln(1-μ_{kd}) ]求导后令其为零得到μ_{kd}^{new} Σ_n γ_{nk} x_{nd} / Σ_n γ_{nk}这同样非常直觉簇k的第d维概率就是把所有样本在第d维的取值1按“属于簇k的后验概率”加权平均。权重归一化后它就是一个0到1之间的概率值。迭代公式简单到让人怀疑但正是这种“E步期望、M步最大化”的交替保证了每一步对数似然都不下降最终收敛到一个局部最大值。3. 用numpy手写一个可运行的EM求解器3.1 生成一份用于验证的二值混合数据纸上谈兵到这里可以停了。为了确认整个流程真的没问题我先用已知参数生成一份人工数据再用EM把它估计回来。这是验证实现是否正确的最稳妥做法。生成方式就是模拟模型本身先按π采样出每条样本所属的簇再按对应的μ向量逐维掷硬币。import numpy as np np.random.seed(42) N, D, K_true 600, 10, 3 pi_true np.array([0.3, 0.4, 0.3]) mu_true np.array([ [0.2, 0.8, 0.3, 0.1, 0.6, 0.9, 0.2, 0.5, 0.1, 0.7], [0.8, 0.2, 0.7, 0.9, 0.1, 0.3, 0.8, 0.4, 0.9, 0.2], [0.4, 0.5, 0.1, 0.6, 0.8, 0.2, 0.3, 0.9, 0.6, 0.4], ]) z_true np.random.choice(K_true, sizeN, ppi_true) X np.zeros((N, D), dtypenp.int8) for n in range(N): for d in range(D): X[n, d] np.random.binomial(1, mu_true[z_true[n], d])这里我设了3个簇、10个维度。三个簇的核心模式分别偏向了不同的维度组合第一个簇集中在第2、6维第二个簇集中在第1、3、4、7维第三个簇分布比较平均但第8维概率很高。这样生成的X肉眼就能看出结构性。你可以直接改写D、K和mu_true用来模拟更高维、更难分的数据。真实项目里没有这个真值来源这里的真值只是为了验证算法没有写错。3.2 初始化和完整EM循环初始化这一步很容易被忽略但EM对初始值非常敏感。如果初始μ全设成0.5所有分量的初始后验概率就几乎一样迭代起来非常慢如果设成极端值又容易过早掉进坏局部最优。我习惯的做法是从数据里随机抽K条样本把每条样本的取值加上一个小扰动当作μ的初始值π则均匀初始化。完整实现如下核心逻辑都放在e_step和m_step里def init_params(X, K, seed0): rng np.random.default_rng(seed) N, D X.shape idx rng.choice(N, K, replaceFalse) mu X[idx].astype(float) 0.1 * rng.random((K, D)) mu np.clip(mu, 1e-3, 1 - 1e-3) pi np.full(K, 1.0 / K) return pi, mu def e_step(X, pi, mu): N, D X.shape K len(pi) log_resp np.zeros((N, K)) for k in range(K): logp np.sum( X * np.log(mu[k]) (1 - X) * np.log(1 - mu[k]), axis1 ) log_resp[:, k] np.log(pi[k]) logp log_resp_max log_resp.max(axis1, keepdimsTrue) log_resp log_resp - log_resp_max resp np.exp(log_resp) resp / resp.sum(axis1, keepdimsTrue) return resp def m_step(X, resp): N, D X.shape K resp.shape[1] Nk resp.sum(axis0) pi Nk / N mu np.zeros((K, D)) for k in range(K): mu[k] resp[:, k] X / Nk[k] return pi, mu def log_likelihood(X, pi, mu): N, D X.shape K len(pi) ll 0.0 for n in range(N): s 0.0 for k in range(K): p pi[k] * np.prod( mu[k] ** X[n] * (1 - mu[k]) ** (1 - X[n]) ) s p ll np.log(s 1e-300) return ll def em(X, K, max_iter200, tol1e-6, seed0): pi, mu init_params(X, K, seed) ll_prev -np.inf for it in range(max_iter): resp e_step(X, pi, mu) pi, mu m_step(X, resp) ll log_likelihood(X, pi, mu) if abs(ll - ll_prev) tol * abs(ll_prev 1e-12): break ll_prev ll return pi, mu, resp, ll几个细节值得说。E步里我先算每个分量下的对数概率加上logπ再用每行最大值做平移最后exp并归一化。这个log-sum-exp操作能把概率的绝对值放大到可计算范围同时不改变归一化结果。M步里用矩阵乘法resp[:, k] X一句代码就算完了加权求和比显式循环快很多。对10个维度、几千条样本的数据跑几百轮EM耗时都在毫秒级完全够日常调试使用。3.3 从日志和参数对比看收敛跑完上面的代码拿估计出来的mu和真实mu做对比是验证算法正确性的关键测试。我在一次运行里得到的估计结果大致是参数真实值估计值按簇对齐后π_10.300.312π_20.400.385π_30.300.303μ_1 第2维0.800.763μ_1 第6维0.900.884μ_2 第1维0.800.772μ_2 第3维0.700.714μ_3 第8维0.900.891由于EM本身的标签置换性算法跑出来的簇编号和真实编号不保证一一对应需要人工按“哪个维度的概率高”进行对齐。但整体数值非常接近这说明实现是对的。迭代过程中对数似然的变化趋势也值得观察。收敛曲线通常在前20轮迅速上升之后进入平台期。如果曲线反复震荡一般有三个原因一是K设得太大导致某些分量互相抢样本二是初始化不当导致后验概率在某个区域反复跳三是数据本身没有明显的簇结构模型本身在强行拟合噪声。4. 实战中躲不开的坑与排查技巧4.1 空分量问题与重启机制EM迭代里最常见的问题是空分量某些簇的后验概率在所有样本上都接近0导致M步里分母Σ_nγ_{nk}特别小μ更新后变得极不稳定甚至下一轮该分量直接消失。出现空分量的根源一般是初始化不好或者K设置得比真实簇数大。数据只有3个天然簇你非要设K5那多出来的2个在迭代中大概率会慢慢没人认领。我常用的处理策略是在每次M步之后做一次检查统计每个分量的有效样本量Nk如果Nk小于某个阈值比如N × 0.01就把这个分量的μ重新初始化为某条随机样本的拷贝并加扰动π设为极小值让它在下一轮重新参与竞争。这个“复活”机制能明显提升高K时的稳定性。阈值本身要按业务场景调整。如果一条样本私有簇本身就有意义你甚至可以把阈值设成1但通常至少设成5~10否则某个分量只被两三个离群点霸占会严重过拟合。4.2 局部最优与多初始化策略EM只能保证收敛到局部最优不能保证全局最优。同一个数据集用不同随机种子初始化跑完可能得到完全不同的簇结构。我在文本聚类项目里遇到过非常典型的情况某次运行分出的簇全都在讲体育另一轮跑出来却有一条是娱乐新闻差别完全来自初始化。后来我固定用“多初始化选对数似然最高”的方案对同样的K换5~10组随机种子跑完整EM最终选对数似然最高的那组结果。这里有个小技巧如果两个初始化的对数似然差别很小而簇结构的业务解释差别很大优先选业务上更可解释的。因为EM的似然差在百分之零点几以内时统计意义并不明显但业务意义很实在。多初始化还能顺带帮你判断K是否合适如果K5和K6的对数似然提升非常微弱那大概率K5已经够了多出来的簇只是在拆解噪声。4.3 数值下溢把乘法挪到对数域二值向量的概率是许多0到1之间的小数连乘。维度上到几百任何一条样本在所有分量下的概率都可能在10^{-100}量级以下直接算连乘必炸。E步、对数似然这两处都必须在对数域操作。我在e_step里已经展示了log-sum-exp的做法。对数似然函数里我用了np.prod连乘这其实只适合维度低的时候验证用。高维场景下应该改成def log_likelihood_stable(X, pi, mu): N X.shape[0] K len(pi) ll 0.0 for n in range(N): log_p np.array([ np.log(pi[k]) np.sum( X[n] * np.log(mu[k]) (1 - X[n]) * np.log(1 - mu[k]) ) for k in range(K) ]) ll np.log(np.sum(np.exp(log_p - log_p.max()))) log_p.max() return ll这种写法对任意维度都安全。它本质上就是在对数域做了一次“软加”先找到最大值平移再归一化加和最后把最大值加回来。4.4 标签置换与聚类结果的解释EM跑出来的簇编号没有任何顺序含义。同一组数据换一次随机种子很可能簇1和簇3就互换了。这在可视化时特别容易误导人你必须根据μ向量的业务含义做对齐。做法很简单跑完多次EM以后以某一组的簇顺序为基准计算各簇的μ向量与基准簇的相似度按相似度做重新排序。二值场景下相似度可以直接用余弦相似度也可以用负的欧氏距离。对齐后再来比较不同初始化之间的结果差异才有意义。实际输出聚类标签时我也会提醒业务方这个编号只是一个簇的代号不是等级的先后顺序。过去就有同事看到簇编号1就以为是“最重要的用户”其实1只是随机分到的而已。5. EM算法还能用在哪些地方5.1 高斯混合模型连续数据版的EM很多人真正听说EM算法是从高斯混合模型开始的。混合伯努利和高斯混合在EM框架下几乎同构差别只在分量分布GMM的每个分量是一个多维高斯参数是均值向量μ_k和协方差矩阵Σ_k伯努利混合的每个分量是一个多维伯努利参数是概率向量μ_k。E步里GMM计算的是“样本由第k个高斯生成”的后验概率M步里均值和协方差的更新公式也能写成后验概率加权平均的形式。所以如果你把混合伯努利的代码彻底搞懂换到GMM只是替换概率计算函数和参数更新公式。几何上GMM适合连续特征、簇形状可能是椭圆的数据混合伯努利适合离散二值特征、数据天然来自“某个模式的概率激活”的场景。选型时看你的特征类型而不是看哪个模型更“高级”。5.2 隐变量模型大家庭与“EM主要用在哪”搜索引擎里经常有人问“em算法主要用在哪”这个问题背后的困惑是EM看起来不是一个具体的算法倒像是一个思路框架。确实如此。EM的正式适用场景是“概率模型中含有隐变量或缺失数据需要做极大似然估计”的问题。除了混合模型HMM里的Baum-Welch算法本质就是EM的变种E步换成前向-后向算法M步更新转移概率和发射概率因子分析、pLSA、部分深度生成模型的推断也都在用EM的变体。更宽泛地说凡是你能写出完整数据似然、但只能观测到其中一部分的问题都可以尝试用EM。它最大的优点是逻辑清晰、实现简单、内存占用小不像变分推断和马尔可夫链蒙特卡洛那样需要大量采样或者复杂近似。缺点是只能保证局部最优而且在隐变量维度很高时收敛会很慢。5.3 伯努利混合和主题模型的关系最后拎出来聊一聊伯努利混合模型和LDA主题模型的关系因为这两个特别容易被搞混。LDA也是处理词袋数据的也涉及到“每个文档属于哪些主题”看起来和混合伯努利很像。核心区别在于伯努利混合模型假设每条样本只属于一个簇也就是文档只有一个主题LDA假设每条文档可以同时包含多个主题且每个词都有自己的主题分配。前者是硬分配的概率版本后者是更精细的多主题混合。如果你的业务场景是“用户只属于一个细分群体”用混合伯努利就够如果一条文本很可能同时聊到体育、经济和娱乐那你需要的是LDA或更现代的主题模型。判断维度就一句话样本内部是不是同质的。是同质的用混合模型不是用多主题模型。我在实际项目中通常不会只用一次EM就交差而是会把“多初始化-选最优-检查空分量-对齐簇标签-解释μ向量”这套流程固化下来变成一个模板函数。往后遇到任何二值聚类需求直接套用省下大量调试时间。如果你自己实现EM也建议把这三件事放在代码里log域计算、多初始化的顶层循环、空分量的重初始化它们比模型本身的精确推导更能决定你这个算法在生产环境里能不能用。
EM算法与混合伯努利模型：二值数据聚类的原理与NumPy实现
我最早把EM算法用在项目里是因为要做一批用户行为向量的聚类。数据是二值化的每个用户对若干品类有买/没买两种状态直接套K-means效果惨不忍睹因为欧氏距离对0/1向量太不友好了簇中心也落在“半买半不买”这种毫无意义的位置上。后来换成混合伯努利模型用EM算法求解聚类结果才真正变得可解释每个簇能给出一套“该品类出现概率”的画像。这篇文章就把这套方法从原理、推导到代码完整地过一遍适合刚接触EM算法的人也适合想快速上手混合伯努利模型做二值数据聚类的同学。1. 混合伯努利模型到底在拟合什么1.1 哪些场景会在意“0/1向量的簇”先想一个问题什么样的问题会产出一堆0/1向量并且你还想给它们分簇最常见的是文本的词袋表示。把一篇文章表示成一个维度为词典大小的向量某一维为1表示这个词在这篇文章里出现过为0表示没出现过。这时候“距离”如果拿欧氏距离来算两篇都提到很多常见词的文章可能距离很近但主题未必一样而提到同样生僻词的两篇文章欧氏距离也很近这个信号反而强烈。这类数据的问题本质不是“空间距离”而是“模式相似度”。用户行为数据也类似。电商场景下每个用户是否购买过某个品类可以编码成一个M维的0/1向量内容平台上每个用户是否阅读过某类内容也可以这样编码。想做用户分群时表面上是聚类实际要的是“把行为模式相同的人归到一起”而不是“把坐标距离近的人归到一起”。图像的二值特征也一样。手写数字的像素二值化、医学影像里的病灶区域二值掩码在降维成特征后经常用这种分布去拟合。这些场景有一个共同点数据本身的生成过程可能是一个“先选组再按组的概率生成0/1”的过程。这就是混合伯努利模型的假设来源。它假设数据集里存在K个隐性的组每个组对应一组伯努利概率参数某条样本先以一定的概率被分配到某个组再根据该组的伯努利参数逐维掷硬币生成0/1。这句话翻译成代码逻辑就是如果已知每条样本属于哪个簇聚类就简单了直接按组统计1的比例就行。问题是我们不知道而且“先选组”这一步也是随机的。于是问题天然带有隐变量。1.2 模型公式和参数含义用更规范的话说混合伯努利模型对单条样本x的概率定义为p(x|θ) Σ_{k1}^{K} π_k · Π_{d1}^{D} μ_{kd}^{x_d} (1-μ_{kd})^{1-x_d}这里的参数包括π_k第k个混合分量的权重满足π_k ≥ 0Σ π_k 1相当于“随机抽取一个样本它属于第k个簇的先验概率”。μ_{kd}第k个簇在第d个维度上取值为1的概率0 ≤ μ_{kd} ≤ 1这就是每个簇的“质心”只不过质心变成了一个概率向量。K簇的数量需要提前指定。D特征维度数。N样本总数。这个式子里Π那一串是一个多维伯努利分布的完整概率它假设在给定簇k的条件下D个维度之间相互独立。这个独立性假设可能不完全符合现实但对聚类任务来说通常足够用因为它真正关心的是“哪些维度组合模式的概率高”。整个数据集的对数似然是所有样本对数似然的加和L(θ) Σ_{n1}^{N} ln( Σ_{k1}^{K} π_k · Π_{d1}^{D} μ_{kd}^{x_{nd}} (1-μ_{kd})^{1-x_{nd}} )如果你直接对这个式子求导并令其等于零会发现问题对数里面套着对K个分量的求和求导之后每一项的分母里都有同一个求和项K个分量的参数互相纠缠无法得到闭式解。这正是EM算法登场的理由。2. EM算法的核心思路把缺失标签先补上再估计2.1 直接最大化对数似然的困难在哪很多讲EM的资料一上来就列公式容易把第一次接触的人吓跑。其实EM解决这个问题的思路非常朴素就一个词缺啥补啥。我们给每条样本引入一个隐变量z_n它表示这条样本真正所属的簇取值1到K。如果把z_n也当成已知数据那么“完整数据”的对数似然就非常好写L_c(θ) Σ_{n1}^{N} Σ_{k1}^{K} [z_{nk} · ( ln π_k Σ_{d1}^{D} ( x_d ln μ_{kd} (1-x_d) ln(1-μ_{kd}) ) )]其中z_{nk}是一个one-hot向量只有样本n真正属于的那个簇对应的位置为1其余为0。这个式子里没有“对K求和套在对数里面”的问题对π、μ求导都能得到闭式解。问题就是z_n是未知的。EM的思路是既然未知那我就猜它的期望用期望代替真实值做参数估计。这就是E步然后用这个期望值去最大化完整数据似然得到新的参数这就是M步。如此反复迭代直到对数似然不再上升。2.2 E步的完整推导E步要算的是后验概率在给定当前参数θ^{(t)}和观测数据x_n的条件下样本n属于簇k的概率。直接用贝叶斯公式γ(z_{nk}) p(z_{nk}1 | x_n, θ^{(t)}) π_k · p(x_n | μ_k) / Σ_{j1}^{K} π_j · p(x_n | μ_j)这里的p(x_n | μ_k)就是第k个伯努利分量对x_n的概率也就是前面Π那一串。这个γ值就是“软标签”。如果某个样本在第4个分量上的后验概率是0.85就说明它大概率属于第4个簇但仍有15%的可能性属于其他簇。这和K-means那种硬分配有本质区别K-means直接一步到位把样本钉死在某一个簇里EM则保留不确定性。实现时有个细节多个伯努利概率相乘连乘几百个接近0或接近1的值之后很容易变成0、inf或者NaN。所以实际工程里几乎从来不在原始概率空间算γ而是把分子分母都搬到对数空间先算log p再做一次log-sum-exp归一化最后exp回去。后面代码部分会展示这个技巧。2.3 M步的参数更新有了γ我们就可以把“完整数据”里的z_{nk}替换成γ(z_{nk})然后对参数求导。先看π。约束条件是Σπ_k 1用拉格朗日乘数法构造拉格朗日函数Λ Σ_n Σ_k γ_{nk} ln π_k λ(Σ_k π_k - 1)对π_k求导并令其为零得到π_k^{new} (1/N) Σ_{n1}^{N} γ_{nk}意思是新权重等于所有样本分配给簇k的平均“软占比”。这很直觉如果所有样本分给簇k的后验概率加起来是180总共600条样本那π_k就是0.3。再看μ_{kd}。完整数据对数似然里和μ_{kd}有关的项是Σ_n γ_{nk} [ x_{nd} ln μ_{kd} (1-x_{nd}) ln(1-μ_{kd}) ]求导后令其为零得到μ_{kd}^{new} Σ_n γ_{nk} x_{nd} / Σ_n γ_{nk}这同样非常直觉簇k的第d维概率就是把所有样本在第d维的取值1按“属于簇k的后验概率”加权平均。权重归一化后它就是一个0到1之间的概率值。迭代公式简单到让人怀疑但正是这种“E步期望、M步最大化”的交替保证了每一步对数似然都不下降最终收敛到一个局部最大值。3. 用numpy手写一个可运行的EM求解器3.1 生成一份用于验证的二值混合数据纸上谈兵到这里可以停了。为了确认整个流程真的没问题我先用已知参数生成一份人工数据再用EM把它估计回来。这是验证实现是否正确的最稳妥做法。生成方式就是模拟模型本身先按π采样出每条样本所属的簇再按对应的μ向量逐维掷硬币。import numpy as np np.random.seed(42) N, D, K_true 600, 10, 3 pi_true np.array([0.3, 0.4, 0.3]) mu_true np.array([ [0.2, 0.8, 0.3, 0.1, 0.6, 0.9, 0.2, 0.5, 0.1, 0.7], [0.8, 0.2, 0.7, 0.9, 0.1, 0.3, 0.8, 0.4, 0.9, 0.2], [0.4, 0.5, 0.1, 0.6, 0.8, 0.2, 0.3, 0.9, 0.6, 0.4], ]) z_true np.random.choice(K_true, sizeN, ppi_true) X np.zeros((N, D), dtypenp.int8) for n in range(N): for d in range(D): X[n, d] np.random.binomial(1, mu_true[z_true[n], d])这里我设了3个簇、10个维度。三个簇的核心模式分别偏向了不同的维度组合第一个簇集中在第2、6维第二个簇集中在第1、3、4、7维第三个簇分布比较平均但第8维概率很高。这样生成的X肉眼就能看出结构性。你可以直接改写D、K和mu_true用来模拟更高维、更难分的数据。真实项目里没有这个真值来源这里的真值只是为了验证算法没有写错。3.2 初始化和完整EM循环初始化这一步很容易被忽略但EM对初始值非常敏感。如果初始μ全设成0.5所有分量的初始后验概率就几乎一样迭代起来非常慢如果设成极端值又容易过早掉进坏局部最优。我习惯的做法是从数据里随机抽K条样本把每条样本的取值加上一个小扰动当作μ的初始值π则均匀初始化。完整实现如下核心逻辑都放在e_step和m_step里def init_params(X, K, seed0): rng np.random.default_rng(seed) N, D X.shape idx rng.choice(N, K, replaceFalse) mu X[idx].astype(float) 0.1 * rng.random((K, D)) mu np.clip(mu, 1e-3, 1 - 1e-3) pi np.full(K, 1.0 / K) return pi, mu def e_step(X, pi, mu): N, D X.shape K len(pi) log_resp np.zeros((N, K)) for k in range(K): logp np.sum( X * np.log(mu[k]) (1 - X) * np.log(1 - mu[k]), axis1 ) log_resp[:, k] np.log(pi[k]) logp log_resp_max log_resp.max(axis1, keepdimsTrue) log_resp log_resp - log_resp_max resp np.exp(log_resp) resp / resp.sum(axis1, keepdimsTrue) return resp def m_step(X, resp): N, D X.shape K resp.shape[1] Nk resp.sum(axis0) pi Nk / N mu np.zeros((K, D)) for k in range(K): mu[k] resp[:, k] X / Nk[k] return pi, mu def log_likelihood(X, pi, mu): N, D X.shape K len(pi) ll 0.0 for n in range(N): s 0.0 for k in range(K): p pi[k] * np.prod( mu[k] ** X[n] * (1 - mu[k]) ** (1 - X[n]) ) s p ll np.log(s 1e-300) return ll def em(X, K, max_iter200, tol1e-6, seed0): pi, mu init_params(X, K, seed) ll_prev -np.inf for it in range(max_iter): resp e_step(X, pi, mu) pi, mu m_step(X, resp) ll log_likelihood(X, pi, mu) if abs(ll - ll_prev) tol * abs(ll_prev 1e-12): break ll_prev ll return pi, mu, resp, ll几个细节值得说。E步里我先算每个分量下的对数概率加上logπ再用每行最大值做平移最后exp并归一化。这个log-sum-exp操作能把概率的绝对值放大到可计算范围同时不改变归一化结果。M步里用矩阵乘法resp[:, k] X一句代码就算完了加权求和比显式循环快很多。对10个维度、几千条样本的数据跑几百轮EM耗时都在毫秒级完全够日常调试使用。3.3 从日志和参数对比看收敛跑完上面的代码拿估计出来的mu和真实mu做对比是验证算法正确性的关键测试。我在一次运行里得到的估计结果大致是参数真实值估计值按簇对齐后π_10.300.312π_20.400.385π_30.300.303μ_1 第2维0.800.763μ_1 第6维0.900.884μ_2 第1维0.800.772μ_2 第3维0.700.714μ_3 第8维0.900.891由于EM本身的标签置换性算法跑出来的簇编号和真实编号不保证一一对应需要人工按“哪个维度的概率高”进行对齐。但整体数值非常接近这说明实现是对的。迭代过程中对数似然的变化趋势也值得观察。收敛曲线通常在前20轮迅速上升之后进入平台期。如果曲线反复震荡一般有三个原因一是K设得太大导致某些分量互相抢样本二是初始化不当导致后验概率在某个区域反复跳三是数据本身没有明显的簇结构模型本身在强行拟合噪声。4. 实战中躲不开的坑与排查技巧4.1 空分量问题与重启机制EM迭代里最常见的问题是空分量某些簇的后验概率在所有样本上都接近0导致M步里分母Σ_nγ_{nk}特别小μ更新后变得极不稳定甚至下一轮该分量直接消失。出现空分量的根源一般是初始化不好或者K设置得比真实簇数大。数据只有3个天然簇你非要设K5那多出来的2个在迭代中大概率会慢慢没人认领。我常用的处理策略是在每次M步之后做一次检查统计每个分量的有效样本量Nk如果Nk小于某个阈值比如N × 0.01就把这个分量的μ重新初始化为某条随机样本的拷贝并加扰动π设为极小值让它在下一轮重新参与竞争。这个“复活”机制能明显提升高K时的稳定性。阈值本身要按业务场景调整。如果一条样本私有簇本身就有意义你甚至可以把阈值设成1但通常至少设成5~10否则某个分量只被两三个离群点霸占会严重过拟合。4.2 局部最优与多初始化策略EM只能保证收敛到局部最优不能保证全局最优。同一个数据集用不同随机种子初始化跑完可能得到完全不同的簇结构。我在文本聚类项目里遇到过非常典型的情况某次运行分出的簇全都在讲体育另一轮跑出来却有一条是娱乐新闻差别完全来自初始化。后来我固定用“多初始化选对数似然最高”的方案对同样的K换5~10组随机种子跑完整EM最终选对数似然最高的那组结果。这里有个小技巧如果两个初始化的对数似然差别很小而簇结构的业务解释差别很大优先选业务上更可解释的。因为EM的似然差在百分之零点几以内时统计意义并不明显但业务意义很实在。多初始化还能顺带帮你判断K是否合适如果K5和K6的对数似然提升非常微弱那大概率K5已经够了多出来的簇只是在拆解噪声。4.3 数值下溢把乘法挪到对数域二值向量的概率是许多0到1之间的小数连乘。维度上到几百任何一条样本在所有分量下的概率都可能在10^{-100}量级以下直接算连乘必炸。E步、对数似然这两处都必须在对数域操作。我在e_step里已经展示了log-sum-exp的做法。对数似然函数里我用了np.prod连乘这其实只适合维度低的时候验证用。高维场景下应该改成def log_likelihood_stable(X, pi, mu): N X.shape[0] K len(pi) ll 0.0 for n in range(N): log_p np.array([ np.log(pi[k]) np.sum( X[n] * np.log(mu[k]) (1 - X[n]) * np.log(1 - mu[k]) ) for k in range(K) ]) ll np.log(np.sum(np.exp(log_p - log_p.max()))) log_p.max() return ll这种写法对任意维度都安全。它本质上就是在对数域做了一次“软加”先找到最大值平移再归一化加和最后把最大值加回来。4.4 标签置换与聚类结果的解释EM跑出来的簇编号没有任何顺序含义。同一组数据换一次随机种子很可能簇1和簇3就互换了。这在可视化时特别容易误导人你必须根据μ向量的业务含义做对齐。做法很简单跑完多次EM以后以某一组的簇顺序为基准计算各簇的μ向量与基准簇的相似度按相似度做重新排序。二值场景下相似度可以直接用余弦相似度也可以用负的欧氏距离。对齐后再来比较不同初始化之间的结果差异才有意义。实际输出聚类标签时我也会提醒业务方这个编号只是一个簇的代号不是等级的先后顺序。过去就有同事看到簇编号1就以为是“最重要的用户”其实1只是随机分到的而已。5. EM算法还能用在哪些地方5.1 高斯混合模型连续数据版的EM很多人真正听说EM算法是从高斯混合模型开始的。混合伯努利和高斯混合在EM框架下几乎同构差别只在分量分布GMM的每个分量是一个多维高斯参数是均值向量μ_k和协方差矩阵Σ_k伯努利混合的每个分量是一个多维伯努利参数是概率向量μ_k。E步里GMM计算的是“样本由第k个高斯生成”的后验概率M步里均值和协方差的更新公式也能写成后验概率加权平均的形式。所以如果你把混合伯努利的代码彻底搞懂换到GMM只是替换概率计算函数和参数更新公式。几何上GMM适合连续特征、簇形状可能是椭圆的数据混合伯努利适合离散二值特征、数据天然来自“某个模式的概率激活”的场景。选型时看你的特征类型而不是看哪个模型更“高级”。5.2 隐变量模型大家庭与“EM主要用在哪”搜索引擎里经常有人问“em算法主要用在哪”这个问题背后的困惑是EM看起来不是一个具体的算法倒像是一个思路框架。确实如此。EM的正式适用场景是“概率模型中含有隐变量或缺失数据需要做极大似然估计”的问题。除了混合模型HMM里的Baum-Welch算法本质就是EM的变种E步换成前向-后向算法M步更新转移概率和发射概率因子分析、pLSA、部分深度生成模型的推断也都在用EM的变体。更宽泛地说凡是你能写出完整数据似然、但只能观测到其中一部分的问题都可以尝试用EM。它最大的优点是逻辑清晰、实现简单、内存占用小不像变分推断和马尔可夫链蒙特卡洛那样需要大量采样或者复杂近似。缺点是只能保证局部最优而且在隐变量维度很高时收敛会很慢。5.3 伯努利混合和主题模型的关系最后拎出来聊一聊伯努利混合模型和LDA主题模型的关系因为这两个特别容易被搞混。LDA也是处理词袋数据的也涉及到“每个文档属于哪些主题”看起来和混合伯努利很像。核心区别在于伯努利混合模型假设每条样本只属于一个簇也就是文档只有一个主题LDA假设每条文档可以同时包含多个主题且每个词都有自己的主题分配。前者是硬分配的概率版本后者是更精细的多主题混合。如果你的业务场景是“用户只属于一个细分群体”用混合伯努利就够如果一条文本很可能同时聊到体育、经济和娱乐那你需要的是LDA或更现代的主题模型。判断维度就一句话样本内部是不是同质的。是同质的用混合模型不是用多主题模型。我在实际项目中通常不会只用一次EM就交差而是会把“多初始化-选最优-检查空分量-对齐簇标签-解释μ向量”这套流程固化下来变成一个模板函数。往后遇到任何二值聚类需求直接套用省下大量调试时间。如果你自己实现EM也建议把这三件事放在代码里log域计算、多初始化的顶层循环、空分量的重初始化它们比模型本身的精确推导更能决定你这个算法在生产环境里能不能用。