企业官网建设流程全解析

如果你翻到课程目录的第九章，看到"EM算法"四个字，多半已经做了点心理准备。我在读研那会儿，第一次看教材里那几页公式，脑子里的画面是：一个没见过数据的老师，非要让一个蒙着眼睛的学生在教室里找黑板擦。后来工作里真正处理带缺失值的数据，才意识到EM其实是个非常朴素的策略：先拿一个不靠谱的初始参数去猜隐藏信息，再用猜出来的信息反过来优化参数，重复到满意为止。它不是什么高深魔法，而是一套把"缺信息"这个问题切成两步的通用框架。这一篇我把整章脉络拆开，讲清楚它为什么这么设计、公式每一步在算什么，再用一个高斯混合模型的小实验把过程跑通，顺便把我在实际建模时踩过的坑也写出来。

1. EM算法：先搞清楚它是来解决什么问题的

1.1 一个看似很简单的估计场景

假设你收到一万条身高数据，想用两个高斯分布去拟合，代表男性和女性的身高分布。理想情况下，如果每条数据都标记了性别，求两个分布的均值和方差就是小菜一碟：对男性样本求均值和方差，对女性样本求均值和方差，收工。这个操作谁都会。

难点在于，性别字段没记录，或者你研究的东西本身就不能直接观测。这时候你面对的数据是X，想推测的类别Z却不可见。如果直接对总似然函数做极大化，会得到一个带 log 内嵌累加的式子。以两个高斯分量为例，目标函数长这样：

log ℓ(θ) = Σᵢ log [ π₁ N(xᵢ|μ₁,Σ₁) + π₂ N(xᵢ|μ₂,Σ₂) ]

问题在于 log 无法穿进加法号，两个子分布的参数缠在一起，求导之后没有漂亮的闭式解。如果Z是可见的，情况完全不同：样本按Z分组后，每组内部的似然可以独立分解，梯度也自然拆开。一旦Z不可见，所有子分布的项都搅合在一个和式里。EM算法本质上就是给这个“纠缠”解套的思路。

1.2 把“一步到位”改成“两步循环”

既然直接求最优θ很难，那就换一种玩法。先随便猜一组初始参数θ₀，这组参数可以很糟糕，但只要有一点点区分能力，就能回答“每条样本更像哪个分布”这个问题。基于θ₀，给每条数据算一个软标签：它属于男生分布的概率是0.7，属于女生分布的概率是0.3。

有了软标签之后，再回头更新参数。这里有个关键点：更新参数时不是按硬分类结果去分组，而是按概率加权。第i个样本对男生组的贡献是0.7，对女生组的贡献是0.3，两边都用到这条数据，只是权重不同。E步和M步交替循环，每轮都让完整数据对数似然的某个下界抬升。用大白话说，E步是在猜身份，M步是在依据身份更新群体画像。

这个过程很像认识新朋友：先凭第一印象判断对方大概属于哪类人，再根据一群人的共性修正印象，然后再重新判断，循环几轮之后，画像会越来越准。

1.3 EM不是无中生有，它是一套通用的“补洞”框架

不少同学学完EM会以为这算法就是为高斯混合模型服务的，这种理解太窄了。任何参数模型，只要似然能够写成“观测变量+隐藏变量”的联合概率积分或求和，并且隐藏变量一旦补齐，模型就变成容易优化的问题，就可以套用EM。

隐马尔可夫模型里的前向-后向训练算法，本质就是EM；LDA主题模型里的变分推断，走的也是类似路径；甚至一些缺失数据填补、推荐系统里的隐因子模型，都能落到完全相同的大框架里。

所以学EM时，我建议大家先吸收建模习惯：判断哪个变量在“偷懒”不肯露面，再把推断分解成E步和M步。只要把模型化成这种完整数据形式，后半段往往只是套用标准求解器的问题。把这一层想明白，第九章的学习曲线会突然变平。

2. 核心公式拆解：E步和M步到底在算什么

2.1 从“对数似然拆不开”说起

公式还是要看一点，但可以看得轻松一些。设观测数据为X，参数为θ，隐变量为Z。我们希望最大化的是：

ℓ(θ) = log p(X|θ) = log Σ_Z p(X,Z|θ)

麻烦就麻烦在log不能直接钻进求和号里。但如果拥有完整数据(X,Z)，要优化的目标变成 log p(X,Z|θ)。因为联合分布通常会按变量关系分解成若干项的乘积，取对数之后就会变成一堆可独立处理的项相加，极大化难度显著降低。

于是问题被分到两层。log 外与求和号之间有一堵墙，但如果我们引入一个分布q(Z)做加权，再利用Jensen不等式，可以得到：

log p(X|θ) ≥ Σ_Z q(Z) log [ p(X,Z|θ) / q(Z) ]

右侧是一个下界。如果把q选成基于旧参数θ_old的后验分布 p(Z|X,θ_old)，那么在θ_old 这一点，下界与原似然函数相等。接下来固定q不动，去寻找能让下界变得更大的新参数θ。

E步就是在做这件事：拿着旧参数θ_old，把隐藏信息的后验分布算出来。M步则是让q保持不变，去最大化这个下界。因为每一步都让原似然的下界上升，而原似然本身有限，最终会收敛到一个局部极大值。

初次接触的话，别强迫自己从下界公式直接跳到底层代码。先记住“E步给权、M步做加权极大似然”这个宏观流程，公式更重要的作用是解释这个流程为什么不跑偏。

2.2 E步里的“期望”到底是对谁求的

E步中经常出现这样的Q函数：

Q(θ; θ_old) = E_{Z|X, θ_old} [ log p(X,Z|θ) ]

这行式子很容易把心理素质差的人劝退。翻译成人话就是：把完整数据的对数似然看成Z的函数，这里Z服从后验分布 p(Z|X, θ_old)，然后对Z求期望。

如果Z是离散的，所谓的期望就是一个按后验概率加权的和。每个样本的Z不是一个确定值，而是一组概率，也就是软标签。所以E步真正输出的，通常不是把一个样本硬塞给某个簇，而是一张N行K列的权重表。例如第i个样本属于第k个簇的权重是0.82，属于另一个簇的权重是0.18。

这一点是和K-means最大的差异。K-means在分配样本时只保留唯一的归属，EM则保留不确定性。正是这种“软”属性，让GMM在类别边界不清晰的数据上通常会比K-means稳定得多。

2.3 M步为什么往往“白给”

把Q函数写出来以后，M步反而显得很平凡。因为log p(X,Z|θ)展开后，许多项只与某个子参数相关，不相关的部分求导后会被直接消掉。拿高斯混合模型来说，Q函数按类别展开后，均值参数μ_k的更新会变成一个典型的带权平均问题：

μ_k_new = (Σᵢ rᵢₖ xᵢ) / (Σᵢ rᵢₖ)

其中rᵢₖ就是E步算出的第i个样本属于第k个高斯的责任权重。分母是第k个聚类总权重，分子是样本加权求和，这正是加权版本的极大似然闭式解。方差矩阵和混合系数也类似，全部有显式公式。

所以完全可以把各种现成极大似然求解器装进M步口袋里使用。如果M步没有闭式解，还可以用一次梯度上升去近似最大化，这类做法叫广义EM，实际项目中用得非常普遍。不要被M步的符号吓到，它本质上就是做一次加权拟合。

2.4 收敛性只需要一个直觉

严格数学证明容易劝退人，但EM的收敛逻辑可以靠直觉掌握。E步和M步每一步都在优化同一个下界函数。当参数位于θ_old时，E步把下界函数挪到与原似然相切的位置；M步沿着下界往上爬。由于下界在旧点是贴着原函数的，所以爬升高下界后，原似然只能跟着往上走，不会掉下去。

这个循环就像登山时先架一座云梯，再顺势往上走一段，反复架桥、前进，总能到达某个山顶。不过，如果山体有很多峰，最终到的是不是最高峰，取决于初始位置。所以我在实际项目中几乎默认使用多组随机初始化，最后挑似然最大的结果。

3. 实操：手写一个GMM聚类，把流程跑通

3.1 先选一个看得见摸得着的例子

教科书一般用二维高斯混合模型来演示EM，因为能画图、能肉眼判断聚类效果，比一上来就在高维数据上折腾直观得多。我这里构造一个包含三个分量的二维数据，每个分量各300个样本，分量之间有一定重叠，再用EM估计三个高斯的均值、协方差和混合系数。这类问题在无监督聚类、异常检测、图像分割里有大量真实投影。

3.2 代码前的三个准备工作

手写EM之前，先做三件事：确定分量数K、初始化参数、定义好数据结构。

K需要提前指定，在真实场景里可以用BIC或轮廓系数去选，在小实验里我直接设K=3。初始化方式我推荐先用K-means跑几轮，把样本粗略分成K组，再以每组均值作为高斯均值μ，以组内样本的协方差作为初始协方差Σ，以组内样本占比作为混合系数α。

这样初始化比完全随机初始化稳定得多。随机初始化如果不配合多次启动，很容易让两个高斯分量在初始阶段抢同一批数据，最终收敛到很差的结果。

3.3 E步实现：先算权重，别忘归一化

E步的核心是计算责任矩阵 r[i,k]，也就是第i个样本属于第k个高斯的概率。代码上可以用 SciPy 的 multivariate_normal 计算多元高斯密度，也可以用 NumPy 手写公式。

先算score[i,k] = αₖ × pdfₖ(xᵢ)，再把每一行除以该行所有score之和。这里必须提醒一个新手的经典错误：直接用样本到均值向量的欧氏距离作为归属依据，忽略了协方差矩阵的形态和混合系数的先验。那样做等于把GMM降级成了K-means，会丢掉软聚类的核心信息。

3.4 M步实现：三个更新公式各司其职

责任矩阵算好之后，M步只需要做三个更新。

更新混合系数αₖ时，把所有样本对第k个高斯的责任相加，再除以总样本数N，得到的就是新的混合权重。更新第k个高斯的均值向量μₖ时，每个样本用责任r[i,k]作为权重，计算加权平均，再从权重总和归一化。更新协方差矩阵Σₖ时，用带责任的样本构造差异向量外积，并做同样的加权平均。

这三个步骤本质上就是带权重的高斯极大似然估计，只是把传统属于某个类的0/1权重换成了[0,1]之间的连续责任值。循环体写完以后，我要在每一轮迭代里打印对数似然，肉眼确认数值在单调上升。如果发现似然下降，那几乎可以确定是代码里某个归一化或取对数的地方写岔了。

3.5 完整示例代码与运行观察

下面是一个最小可运行版本，只依赖 NumPy 和 SciPy：

import numpy as np from scipy.stats import multivariate_normal # 生成三簇模拟数据 np.random.seed(42) true_means = [np.array([-2, 0]), np.array([2, 0]), np.array([0, 3])] true_covs = [np.array([[0.6, 0.2], [0.2, 0.5]]), np.array([[0.6, -0.2], [-0.2, 0.5]]), np.array([[0.5, 0.1], [0.1, 0.4]])] X = [] for m, c in zip(true_means, true_covs): X.append(np.random.multivariate_normal(m, c, 300)) X = np.vstack(X) N, D = X.shape K = 3 # 用K-means初始化参数 from sklearn.cluster import KMeans km = KMeans(n_clusters=K, random_state=0).fit(X) labels = km.labels_ alpha = np.array([np.mean(labels == k) for k in range(K)]) mu = np.array([X[labels == k].mean(axis=0) for k in range(K)]) cov = np.array([np.cov(X[labels == k], rowvar=False) + 1e-6 * np.eye(D) for k in range(K)]) def log_likelihood(X, alpha, mu, cov): n = X.shape[0] log_ll = 0.0 for i in range(n): row = np.array([ alpha[k] * multivariate_normal.pdf(X[i], mean=mu[k], cov=cov[k]) for k in range(K) ]) log_ll += np.log(row.sum()) return log_ll for it in range(100): # E-step: 计算责任矩阵 r[i, k] r = np.zeros((N, K)) for k in range(K): r[:, k] = alpha[k] * multivariate_normal.pdf(X, mean=mu[k], cov=cov[k]) r = r / r.sum(axis=1, keepdims=True) # M-step: 更新参数 Nk = r.sum(axis=0) alpha_new = Nk / N mu_new = (r.T @ X) / Nk[:, None] cov_new = np.zeros((K, D, D)) for k in range(K): diff = X - mu_new[k] cov_new[k] = (r[:, k][:, None] * diff).T @ diff / Nk[k] cov_new[k] += 1e-6 * np.eye(D) alpha, mu, cov = alpha_new, mu_new, cov_new if it % 10 == 0: print("iter", it, "log-likelihood:", log_likelihood(X, alpha, mu, cov))

运行的时候，你会看到似然函数在前面几次迭代快速爬升，随后增量变小，最终趋于平缓。由于这个实验的数据是从已知参数生成的，收敛得到的中心会和真实中心非常接近。我实测时还发现，如果完全换成随机初始化，偶尔会让两个高斯分量在数据中互相纠缠，最后收敛到局部最优，导致某个簇被切分。这个问题也是下一节要聊的重头戏。

代码里写逐样本循环，是为了让公式对应关系直白，数据量大了以后需要改成矩阵化计算，不然会很慢。但作为学习版本，清晰度优先。

4. 实战中必踩的坑与排查手段

4.1 局部最优：每次结果都不一样

EM对初始参数非常敏感，尤其当K偏大、数据分布重叠度高时，同一份数据跑十次，很可能得到几种截然不同的聚类结果。

我做项目的标准操作是：先用K-means初始化，这是降低随机性的第一步。第二步再做多次启动，比如10次甚至20次，每次固定不同随机种子，最后取对数似然最大的那一组参数。第三步，如果业务对稳定性有硬性要求，还可以先做标准化或PCA降维，让数据分布更接近各项同性，EM跑起来会稳一些。

如果发现不同启动的结果差异极大，优先怀疑初始化太差，而不是代码写错了。这是EM的固有特性，不需要慌。

4.2 协方差矩阵奇异：大多是某个分量样本太少

运行GMM时最常见的报错是"LinAlgError: Singular matrix"，源头往往是某一个高斯分量的责任权重总和接近0或为0，导致协方差矩阵无法求逆。当K设得比数据真实分布多，或者初始化把某一簇里样本扔得太少时，就会出现这种问题。

最简单的应对是给每次M步更新后的协方差矩阵加一个很小的正则项，例如1e-6乘以单位矩阵。这个操作在代码里已经列出来了，虽然增加了极小的偏差，但能稳定矩阵求逆。如果加了正则仍然报错，就要检查Nk是不是掉到了0附近。若是这样，减少K，或者用更保守的初始化会更有效。

4.3 收敛慢或震荡：先检查是不是用了硬标签

收敛过程如果出现似然震荡，而不是稳定上升，通常是在E步做了多余操作。比如有人觉得软标签不如硬标签清晰，手动把责任矩阵里最大值置1、其余置0，拿硬EM去跑。那样做不是完全不行，但它已经不是标准EM，更容易出现震荡。

还有一种常见问题是没有对责任矩阵做行归一化。如果责任矩阵每一行的和不等于1，M步的若干加权平均就会失去概率意义，参数更新方向会变得奇怪。排查时打印几行责任矩阵，确认每个元素都在0到1之间，并且行和严格等于1。

4.4 分量标签翻转和可视化陷阱

高斯混合模型本身不区分哪个分量排第一、哪个排第二。同一组聚类结果，把分量1和分量2互换，似然不会变化。这会给代码调试带来困惑：两次运行聚类效果一样，但中心点看起来对不上，甚至标签名字都对不上。

遇到这种情形，不要用分量序号直接与真实簇类别做比较，而是使用调整兰德指数、归一化互信息等不受标签顺序影响的指标。如果必须对比参数，就给分量做一个按位置排列的后处理，比如按均值向量在空间中的位置排序，否则对比结果很容易自我怀疑。

4.5 决定K值：不只看图，还要算账

GMM的K值是需要预先指定的。如果数据本身有四类，你偏用K=3，算法会把两类糊成一团。如果你用K=5，多余的分量又会去抢数据，甚至某个高斯形状会被拉成细长条来拟合离群点。

一套比较可复现的决策方式是画BIC或AIC曲线，看指标下降的拐点。BIC的计算并不麻烦，在最优参数时给对数似然加一个与自由参数数量的惩罚项。比起纯粹人眼看散点图，BIC更客观，也更容易写进流程。

不过我要多提醒一句：无监督场景里统计上的最优K，不一定等于业务想要的分类数。之前做个客服文本聚类项目，统计结果偏好K=8，但业务侧只需要5个标签，那就拿K=8的结果做归类合并，把细节留到后续人工打标，而不是机械地追求统计最优。

5. 走出课本：EM在真实项目里的变体与体会

5.1 不止是聚类，缺失值填补也能用

日常数据处理时，我们经常遇到部分字段缺失的情况。设备传感器可能在某段时间断传，用户画像里某个字段可能成片缺失。直接删掉缺失行浪费信息，用均值填充又太粗糙。如果能将缺失值视作隐变量，将已观测字段视作已知部分，完全可以用EM框架估计出不带主观先验的填补结果。

对于待填补的连续变量，E步要计算缺失变量的条件期望，M步基于期望后的补全结果重新拟合回归或分类参数。这套逻辑比独立插值更能捕捉变量间的关联结构。不过要记牢，普通EM默认缺失机制可以被忽略，如果字段是否缺失与字段值本身有关，例如高收入用户更不愿意填收入，那么直接用EM会有偏差，需要引入更复杂的模型才能处理。

5.2 从GMM到HMM、LDA，方法异曲同工

EM骨架可以顺滑地搬到许多经典模型上。隐马尔可夫模型的Baum-Welch算法本质上就是EM：E步用前向-后向算法估计每个时刻落在某个隐状态的后验概率，M步再重估状态转移矩阵和发射概率。LDA主题模型里的训练过程也与EM同源，唯一的区别在于它需要用变分分布去近似真正的后验分布，不然E步可能不可解。

在第九章掌握的那套思路，能帮你快速理解其他高级模型的训练步骤。如果真的只是为了用库，多数时候我们不需要手写EM，库内部已经完成了相关计算。但理解原理仍然能帮我们判断参数设置是否合理，以及当模型结果异常时，排查问题出在E步还是M步。

5.3 手推小例子和调包之间的分寸

我自己的习惯是，学习任何迭代式算法，先在小数据集上手动迭代三到五轮，把每轮参数的变化写下来。EM尤其适合这种练法：用两个高斯分量、几十个样本，把责任矩阵逐轮打印出来，你会非常直观地看到软标签怎样从先验均衡走向分化，均值向量怎样一步一步移动到聚类中心。

等这个流程做完，再看网上的高速矩阵版本，或者看框架源码里的加速技巧，都会有底气。工作后我用现成高斯混合库的机会很多，但每当效果异常，例如某个类始终不饱满，我第一反应仍然是调出中间责任矩阵看几眼，看看能量是不是被某个分量吸走了。这种“打开黑盒”的直觉，来自当初那几页手推笔记。所以我的建议很朴素：第九章别急着跳到sklearn里的GaussianMixture，先花半小时手写一个不加加速的最小版本，收益比刷十遍教学视频都高。