我最早把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|θ) = Σ_{k=1}^{K} π_k · Π_{d=1}^{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(θ) = Σ_{n=1}^{N} ln( Σ_{k=1}^{K} π_k · Π_{d=1}^{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(θ) = Σ_{n=1}^{N} Σ_{k=1}^{K} [z_{nk} · ( ln π_k + Σ_{d=1}^{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) / Σ_{j=1}^{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) Σ_{n=1}^{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, size=N, p=pi_true) X = np.zeros((N, D), dtype=np.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, seed=0): rng = np.random.default_rng(seed) N, D = X.shape idx = rng.choice(N, K, replace=False) 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]), axis=1 ) log_resp[:, k] = np.log(pi[k]) + logp log_resp_max = log_resp.max(axis=1, keepdims=True) log_resp = log_resp - log_resp_max resp = np.exp(log_resp) resp /= resp.sum(axis=1, keepdims=True) return resp def m_step(X, resp): N, D = X.shape K = resp.shape[1] Nk = resp.sum(axis=0) 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_iter=200, tol=1e-6, seed=0): 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做对比，是验证算法正确性的关键测试。我在一次运行里得到的估计结果大致是：
| 参数 | 真实值 | 估计值（按簇对齐后） |
|---|---|---|
| π_1 | 0.30 | 0.312 |
| π_2 | 0.40 | 0.385 |
| π_3 | 0.30 | 0.303 |
| μ_1 第2维 | 0.80 | 0.763 |
| μ_1 第6维 | 0.90 | 0.884 |
| μ_2 第1维 | 0.80 | 0.772 |
| μ_2 第3维 | 0.70 | 0.714 |
| μ_3 第8维 | 0.90 | 0.891 |
由于EM本身的标签置换性，算法跑出来的簇编号和真实编号不保证一一对应，需要人工按“哪个维度的概率高”进行对齐。但整体数值非常接近，这说明实现是对的。
迭代过程中对数似然的变化趋势也值得观察。收敛曲线通常在前20轮迅速上升，之后进入平台期。如果曲线反复震荡，一般有三个原因：一是K设得太大导致某些分量互相抢样本；二是初始化不当导致后验概率在某个区域反复跳；三是数据本身没有明显的簇结构，模型本身在强行拟合噪声。
4. 实战中躲不开的坑与排查技巧
4.1 空分量问题与重启机制
EM迭代里最常见的问题是空分量：某些簇的后验概率在所有样本上都接近0，导致M步里分母Σ_nγ_{nk}特别小，μ更新后变得极不稳定，甚至下一轮该分量直接消失。
出现空分量的根源一般是初始化不好，或者K设置得比真实簇数大。数据只有3个天然簇，你非要设K=5，那多出来的2个在迭代中大概率会慢慢没人认领。
我常用的处理策略是在每次M步之后做一次检查：统计每个分量的有效样本量Nk，如果Nk小于某个阈值，比如N × 0.01，就把这个分量的μ重新初始化为某条随机样本的拷贝并加扰动，π设为极小值，让它在下一轮重新参与竞争。这个“复活”机制能明显提升高K时的稳定性。
阈值本身要按业务场景调整。如果一条样本私有簇本身就有意义，你甚至可以把阈值设成1；但通常至少设成5~10，否则某个分量只被两三个离群点霸占，会严重过拟合。
4.2 局部最优与多初始化策略
EM只能保证收敛到局部最优，不能保证全局最优。同一个数据集，用不同随机种子初始化，跑完可能得到完全不同的簇结构。
我在文本聚类项目里遇到过非常典型的情况：某次运行分出的簇全都在讲体育，另一轮跑出来却有一条是娱乐新闻，差别完全来自初始化。后来我固定用“多初始化+选对数似然最高”的方案：对同样的K，换5~10组随机种子跑完整EM，最终选对数似然最高的那组结果。
这里有个小技巧：如果两个初始化的对数似然差别很小，而簇结构的业务解释差别很大，优先选业务上更可解释的。因为EM的似然差在百分之零点几以内时，统计意义并不明显，但业务意义很实在。多初始化还能顺带帮你判断K是否合适：如果K=5和K=6的对数似然提升非常微弱，那大概率K=5已经够了，多出来的簇只是在拆解噪声。
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域计算、多初始化的顶层循环、空分量的重初始化，它们比模型本身的精确推导更能决定你这个算法在生产环境里能不能用。
