1. 先说清楚为什么K-means解决不了的问题，GMM可以

1.1 K-means的“圆形偏见”

有一次做用户分群，我拿到一张二维散点图：很明显三堆数据，但中间那一堆是斜着拉长的椭圆，另外两堆相对规则。K-means我调了好几次，每次中间那堆椭圆都会被拦腰截断，分出来的两个簇各占半个椭圆，业务上完全没法解释。也就是从那时候开始，我把高斯混合模型的参数估计与聚类应用从头捋了一遍，才发现这类问题从一开始就不该用K-means硬扛。

K-means的思路大家都很熟：初始化K个中心，把每个样本分给最近的中心，重新计算中心，循环。它最小化的是每个样本到所属中心的平方距离和。这个目标函数虽然简单，但背后隐含了极强的几何假设——所有簇在特征空间里看起来都像大小相近的圆球。严格说，K-means等价于一个协方差矩阵为相同各向同性矩阵、混合权重均匀的GMM在这种约束下做“硬聚类”的极限情况。也就是说，它根本不知道某个簇可以在某个方向上更舒展，也不知道某个簇可以比另一个簇更稀疏，更不知道簇和簇之间密度差异有多大。

所以当数据出现方向性、密度差异、或者簇与簇有重叠时，K-means的切分就非常失真。它只会按欧氏距离画垂直平分线，再把空间切成若干块，边界都是直线。我那组椭圆数据就是这么翻车的：椭圆被拦腰截断，两半分别并到了旁边的簇里去，怎么看怎么不对劲。

1.2 GMM给出的三个升级点

高斯混合模型相比K-means，最核心的差异在于三个升级点。

第一，软分配。GMM输出的是每个样本属于每个簇的后验概率，而不是一个非此即彼的标签。交界处的点到底更偏向哪个簇，模型会把不确定性显式表达出来。这个在业务上意义很大，比如用户分群时，边界用户往往是运营重点去研究的对象。

第二，形状表达力。每个簇都有自己的协方差矩阵Σ_k，可以描述“这个簇在哪个方向上更分散、旋转角度是多少”。K-means只关心均值，GMM还关心方差和协方差结构，自然能识别椭圆簇、不同密度的簇。

第三，生成模型能力。GMM拟合完后得到的是一个完整的概率密度函数，后面可以做密度估计、异常检测、随机采样。K-means给不了这些东西，它就是一个分割算法。

这三个升级点不是花架子，实际项目里都会用到：协方差矩阵能解释业务里的“这类用户波动大，那类很稳定”；密度估计是异常检测的基础；概率输出则是风险评估类产品最需要的东西。

1.3 什么时候选GMM，什么时候不选

GMM也不是万能钥匙。如果簇分离得很干净、每个簇都接近圆球，K-means又快又稳，没有必要强行上GMM。但数据出现以下信号时，GMM基本就是首选：簇之间有重叠、簇的形状差异明显、你需要概率输出而不是硬标签、后面还要做密度估计。还有一种典型场景是你需要的不是聚类本身，而是对一组数据的分布建模，中间夹杂几个异常点，GMM也能当密度模型用。

反过来，GMM对环状、月牙形、流形这类非凸几何很无力。每个高斯分量只有一个峰，一个环没法用一个高斯描述清楚。这类数据要么拼很多分量硬凑，要么直接换DBSCAN、谱聚类。很多人在环形数据上硬套GMM然后抱怨不好用，其实是选错了模型。

2. 高斯混合模型到底在算什么：似然、隐变量与EM的直觉

2.1 从“几堆数据混在一起”到混合概率密度

GMM的假设很直观：你看到的样本，来自K个高斯分布按一定比例混合的总体。写成公式就是：

p(x) = Σ_{k=1}^{K} π_k · N(x | μ_k, Σ_k)

其中π_k是第k个分量占所有数据的比例，满足π_k ≥ 0且所有π_k之和为1；N(x | μ_k, Σ_k)是均值为μ_k、协方差矩阵为Σ_k的多元高斯密度。每个高斯分量就是一团密度云，中心是μ_k，形状和朝向由Σ_k决定。

理解生成过程更简单：每来一个样本，先按π_1到π_K掷一次骰子，选中第k个分量，再从这个分量的高斯分布里抽一个点。反复抽，最终得到你手里的数据集。所以“混合”两个字指的不是数据堆在一起，而是产生机制上是多组分叠加。

这套建模思路对实际问题的映射很直接：比如用户群可能有高活跃、中活跃、低活跃三种人群，每种人群的行为指标各自围绕自己的中心波动，整体分布就是三个高斯的混合。

2.2 为什么最大似然估计在这里“卡壳”

如果每个样本x_n还带着一个标签z_n，告诉你它来自第几个分量，那参数估计就太简单了：把同一个分量的样本归到一起，求均值、求协方差、算占比，都是初等统计。

但聚类场景里z_n看不到，我们能用的只有对数似然：

log p(X) = Σ_n log( Σ_k π_k · N(x_n | μ_k, Σ_k) )

问题出在log里面有一个求和。对μ_k求导之后，这一项依然和别的分量缠在一起，没有办法直接解开得到一个闭式解。换句话说，对数似然的梯度方程是耦合的，不能像单高斯那样一步解出参数。这正是GMM参数估计和拟合单个高斯分布本质不同的地方。

2.3 EM算法的直觉：先猜后改的“软标签”循环

EM算法解决这个问题的思路，说起来很小学生：标签和参数互相纠缠，那就先固定一边猜另一边。

具体说就是：先随便初始化一组参数，然后根据当前参数计算每个样本属于每个分量的概率（这叫E步），再用这些概率作为软标签去重新估计参数（这叫M步）。重复这个过程，每一步都会让对数似然上升，直到收敛。

可以类比成两个瓶子里混着不同甜度的糖，人蒙着眼睛，既看不到每颗糖来自哪个瓶子，也不记得每个瓶子的平均甜度。那就先猜两个平均甜度，然后按甜度距离把每颗糖软性地分给两个瓶子；分完之后，又用这些软归属重新估算每个瓶子的平均甜度；再归一次，再算一次……甜度参数就这么一步步逼近真实值。唯一要记住的是，EM保证对数似然单调不减，但它只能保证收敛到局部最优，不能保证全局最优，所以初始化特别关键。

3. 手写EM做参数估计：公式加代码一起讲透

3.1 初始化策略

EM从一个初始参数开始。直接用随机点当均值，很容易让几个分量跑去做同样的事，或者某个分量完全落空。我自己最常用的初始化方法是先跑一遍K-means：把K-means得到的簇中心作为μ_0，把每个簇里样本相对中心的平均平方偏差作为Σ_0，把簇容量占比作为π_0。这样迭代起点离真实解通常很近，EM基本十几轮就能收敛。sklearn里GaussianMixture的init_params='kmeans'也是这个思路。

3.2 E步与M步公式

E步，对每个样本x_n和每个分量k，计算责任度γ_nk。责任度表示第n个样本属于第k个分量的后验概率：

γ_nk = π_k · N(x_n | μ_k, Σ_k) / Σ_j π_j · N(x_n | μ_j, Σ_j)

分子是这个分量自己贡献的密度，分母是所有分量的密度之和，做一次归一化。这组值加在一起就是样本的软分配，不再是非黑即白的标签。

M步则是把这些责任度当软权重，重新估计三个参数。先把所有样本在分量k上的责任度加起来，得到软计数N_k = Σ_n γ_nk。然后：

π_k = N_k / N

μ_k = (Σ_n γ_nk · x_n) / N_k

Σ_k = (Σ_n γ_nk · (x_n - μ_k)(x_n - μ_k)^T) / N_k

这套公式就是加权版本的占比、均值、协方差。权重不是0或1，而是每个样本属于这个分量的后验概率。

3.3 完整代码：从生成数据到EM迭代

下面这段代码是教学向的完整实现，保留标准EM流程，没有做太多工程优化。生成两簇有明显重叠的二维高斯数据，然后用EM去估计参数。

import numpy as np from scipy.special import logsumexp from sklearn.cluster import KMeans def log_norm_pdf(x, mu, cov): d = len(mu) sign, logdet = np.linalg.slogdet(cov) diff = x - mu quad = np.einsum('ni,ij,nj->n', diff, np.linalg.inv(cov), diff) return -0.5 * (d * np.log(2 * np.pi) + logdet + quad) def em_gmm(X, K, max_iter=100, tol=1e-4): N, D = X.shape # 用 KMeans 做初始化：均值、协方差、混合权重 kmeans = KMeans(n_clusters=K, n_init=5, random_state=0).fit(X) mu = kmeans.cluster_centers_ labels = kmeans.labels_ pi = np.bincount(labels, minlength=K) / N cov = np.zeros((K, D, D)) for k in range(K): Xk = X[labels == k] if len(Xk) < 2: cov[k] = np.eye(D) else: cov[k] = np.cov(Xk.T) + 1e-6 * np.eye(D) old_loglik = None for it in range(max_iter): # E 步：计算责任度，所有计算放在 log 空间防下溢 log_resp = np.empty((N, K)) for k in range(K): log_resp[:, k] = np.log(pi[k]) + log_norm_pdf(X, mu[k], cov[k]) logp = logsumexp(log_resp, axis=1) resp = np.exp(log_resp - logp[:, None]) # M 步：更新混合权重、均值、协方差 Nk = resp.sum(axis=0) pi = Nk / N mu = (resp.T @ X) / Nk[:, None] for k in range(K): diff = X - mu[k] cov[k] = (resp[:, k][:, None] * diff).T @ diff / Nk[k] cov[k] += 1e-6 * np.eye(D) # ridge 防止奇异 loglik = logp.sum() if old_loglik is not None and abs(loglik - old_loglik) < tol: break old_loglik = loglik return pi, mu, cov, resp np.random.seed(42) n1, n2 = 800, 1200 true_mu1 = np.array([-2.0, 0.0]) true_cov1 = np.array([[1.0, 0.6], [0.6, 1.0]]) true_mu2 = np.array([1.5, 1.5]) true_cov2 = np.array([[0.8, 0.0], [0.0, 1.5]]) X = np.vstack([ np.random.multivariate_normal(true_mu1, true_cov1, n1), np.random.multivariate_normal(true_mu2, true_cov2, n2) ]) pi, mu, cov, resp = em_gmm(X, K=2) print("混合权重估计：", pi) print("均值估计：\n", mu) print("协方差估计：\n", cov)

跑完后可以看到估计结果和真实参数非常接近。因为初始化已经用了K-means贴近真实中心，通常20轮以内就会触发收敛条件。如果换成纯随机初始化，轮次会明显变多，而且更可能收敛到局部最优。

3.4 实现里的数值细节

这段代码里我故意保留了两个工程细节，都是实际手写GMM最容易踩的坑。

第一，所有密度计算都放到log空间，再用logsumexp归一化。如果直接在指数空间算π_k · N(x|μ_k,Σ_k)，维度稍高点、数据点稍多点，乘积很快就下溢成0。一旦在归一化时碰上0/0，得到的责任度全是nan，后续迭代全部报废。

第二，协方差矩阵每次更新后加了一个1e-6的对角ridge。目的是防止迭代过程中某个分量退化到奇异矩阵，尤其是某个分量只捕获到极少数样本的时候。加了ridge之后，矩阵数值稳定性会好很多。教学代码里我直接用np.linalg.inv求逆，批量运算对性能不友好。工程上更推荐先对协方差做Cholesky分解再算密度，只是那样代码会绕，教学场景能看懂逻辑更重要。

4. 用sklearn跑真实聚类：二维合成数据上的效果对比

4.1 生成一组让K-means翻车的数据

手写EM能帮你理解原理，实际干活的时候多半还是用sklearn。我准备造一个三分类数据集：一个斜椭圆、一个正圆、一个细长条，故意让椭圆的长轴方向和坐标轴成45度，让细长条在y方向上发散很大。这种构造就是为了逼出协方差矩阵的表达力。

import numpy as np from sklearn.cluster import KMeans from sklearn.mixture import GaussianMixture np.random.seed(7) X1 = np.random.multivariate_normal([0, 0], [[3.0, 1.5], [1.5, 1.0]], 400) X2 = np.random.multivariate_normal([6, 0], [[1.0, 0.0], [0.0, 1.0]], 300) X3 = np.random.multivariate_normal([3, -3], [[0.4, 0.0], [0.0, 2.0]], 300) X = np.vstack([X1, X2, X3])

4.2 GMM与K-means聚类结果对比

KMeans(n_clusters=3, n_init=10)和GaussianMixture(n_components=3, covariance_type='full')分别跑一遍，观察簇边界。

K-means对斜椭圆簇的处理几乎固定是“从腰部连线切开”，因为欧氏距离最近的中心决定了归属，两个新簇各占椭圆一半，椭圆形状信息完全丢失。GMM则通过协方差矩阵学习到了椭圆的主方向，边界会跟着椭圆等密度线走，把斜椭圆完整保留下来。

两类模型簇之间边界的形态也不同：K-means的边界是直线，GMM因为各簇协方差不同，边界是二次曲线。这个差异在高维数据里会更明显，也是为什么GMM对复杂形状簇的还原度比K-means高一个档次。

4.3 协方差类型的选择：一表讲清楚

sklearn的GaussianMixture可以通过covariance_type参数控制协方差矩阵的自由度。如果你不知道这四个选项背后的权衡，建议先看下面这张表。

实在拿不准，我建议在数据量充足时先跑full，然后用BIC比较tied和diag的得分。高维场景下full的参数量会爆炸，容易过拟合，换成diag或tied往往更稳。这里也说明一个问题：GMM不是“一定比K-means复杂就一定更好”，复杂度换来的是表达力，但也要承受过拟合风险。

5. 换个场景实战：图像分割与异常检测

5.1 图像分割：把像素当成样本做软聚类

GMM做图像分割的思路是：把每个像素当成一个D维样本，D是颜色通道数，用GMM把整张图聚成K类，再给每个像素贴上后验概率最高的簇标签。最后用簇中心颜色替换每个像素的原始颜色，图像就变成一种“颜色压缩”版本，肉眼看上去就是分割结果。

实操里有几个细节需要处理。

第一，像素量太大。普通图片几十万甚至上百万个像素点，直接GMM全量拟合会比较慢。稳妥做法是随机抽样10%~20%的像素训练模型，然后把模型apply到全图。如果数据流式读取，也可以直接用MiniBatchGaussianMixture，它支持小批量更新。

第二，颜色空间选择。RGB空间三个通道相关性很强，且欧氏距离不贴合人眼感知。更推荐转换到Lab颜色空间再做聚类，因为Lab空间里的欧氏距离更接近视觉差异。变换后的分割边界通常更干净。

第三，K取值。BIC能给出大致范围，但图片内容影响极大：天空、草地这种简单场景K取3左右就够，复杂纹理可能要到10以上。K取大时容易出现过度分割，同一个物体被劈成好几块。我做过某类组织切片图的前景背景分离，聚类结果比固定阈值好很多，因为组织颜色分布天然就是两团高斯的混合。但稍微把K调大一点，同一个组织就被分成好几块了，最后还是降回K=2才解决问题。

5.2 异常检测：用score_samples揪出低概率样本

GMM本质上就是密度估计，密度越低的区域越可能是异常点。sklearn的GaussianMixture提供了score_samples方法，返回每个样本的对数似然。常规做法是用正常数据拟合GMM，再对测试集算得分，然后取某个分位数作为阈值，低于阈值的标记为异常。

from sklearn.mixture import GaussianMixture gm = GaussianMixture(n_components=2, covariance_type='full', random_state=0) gm.fit(X_normal) scores = gm.score_samples(X_test) threshold = np.percentile(scores, 5) anomaly_mask = scores < threshold

这个方案和IsolationForest、DBSCAN的取舍在于：GMM假设正常数据本身符合高斯混合分布，如果你的正常数据天然有多个聚集区，GMM会给出最自然的密度边界；如果数据模式是任意形状，比如环状、马鞍形，GMM的分类边界会很别扭。另外高维环境下“密度估计”本身会变得非常不可靠，这涉及到维度灾难，通常建议先降维到合适维度再建GMM。我自己做日志指标异常检测时就是这么干的：正常模式的指标空间先PCA降到二维三维，再用GMM定阈值，效果比在高维空间硬算稳定得多。

6. 调参避坑记录：这些坑我不希望你再去踩

6.1 聚类数K怎么定：BIC/AIC不是万能钥匙

GMM的K不能靠看聚类结果拍脑袋。常用方法是BIC和AIC，sklearn里直接调用gm.bic(X)和gm.aic(X)就行。BIC对复杂度惩罚更重，选出的K通常比AIC小。

一个非常容易被忽略的坑是：同一个K用不同random_state跑，BIC分数可能差别很大。如果只跑一次就下结论，很容易因为初始化撞到局部最优而选错K。我一般每个K跑5到10个随机种子，取BIC中位数比较。

best_k, best_bic = None, np.inf for k in range(1, 11): vals = [] for seed in range(5): gm = GaussianMixture(n_components=k, covariance_type='full', random_state=seed) gm.fit(X) vals.append(gm.bic(X)) median_bic = np.median(vals) if median_bic < best_bic: best_bic, best_k = median_bic, k

选完之后还不是终点。BIC只是统计指标，业务解释更关键。如果某个分量只占了总体样本量的2%，哪怕BIC分数很好看，这个簇在业务上通常也很难自圆其说，这时候宁愿降K。

6.2 初始化与局部最优：random_state带来的“幸存者偏差”

EM的收敛点完全依赖起点。同样一份数据，random_state=0和random_state=42跑出来的聚类结果可能差别很大，而且两个结果看起来都像模像样。这就是局部最优问题。它不像K-means那样可以通过多次重启直观解决，因为每次重启EM的代价更高。

我的做法有两条。第一，init_params用'kmeans'，先用K-means找到较好起点，实测比'random'稳定得多。第二，n_init至少设成5到10，模型内部会跑多组初始化，最终选对数似然最高的结果。运行时固定random_state，方便复现和调参。

还有一个容易被忽略的关联问题：喂给GMM的数据是否标准化，会直接影响K-means初始化的合理性。特征量纲差异大时，K-means会被高方差维度带偏，进而带偏GMM。所以我的流程都是先标准化，再跑初始化，再建GMM。

6.3 数值稳定性：奇异协方差矩阵的预防

GMM训练时最经典的问题就是协方差矩阵变成奇异矩阵。症状有两种：要么某次迭代后对数似然突然变成nan，要么某个分量的协方差矩阵行列式变成负数或者0。常见诱因有三个：某个分量只捕获到一个样本；数据特征高度共线；两个分量均值几乎重合且样本量极少。

对策也分三层。第一层，数据预处理时检查有没有常数特征，有就删掉。标准化之后再看特征相关性，严重时做PCA去相关。第二层，给协方差加ridge，sklearn里对应reg_covar参数。默认值1e-6够用，但如果数据本身比较脏，可以放宽到1e-4甚至1e-3。第三层，训练后检查每个分量的责任度之和，几乎为0的分量直接删掉或者降K重跑。数值问题不解决，后面所有分析都建立在沙子上。

6.4 分量退化：一个簇只吃到一个样本怎么办

和奇异矩阵紧密相关的是分量退化。GMM迭代过程中，某个分量可能慢慢把所有责任度都收到一个样本上，协方差越缩越小，最终趋近于零矩阵。这个分量实际上既不代表一个簇，也没有任何解释意义，只是过拟合了单个噪声点。

最直接的观察方式是打印分量责任度之和：

responsibilities = gm.predict_proba(X) print(responsibilities.sum(axis=0))

正常情况下每个分量的责任度应该和混合权重一致，数量级上没有太离谱的偏差。如果有一项明显低于总体样本量的2%，说明这个分量在吃噪声。处理办法包括：降K、给协方差加扰动、或者对混合权重设置一个最小值约束。实际项目里我更多时候选择降K，因为删除一个“废分量”的代价远低于强行保留它然后给业务解释一个空壳簇。

6.5 我的常用流程和一个收尾习惯

把前面这些经验串成一个固定流程：先标准化数据，再PCA或者t-SNE看分布形态；确定K的候选范围后，用多个seed跑BIC选K；建模型时init_params='kmeans'，n_init设5以上，固定random_state；拟合完打印责任度之和，检查是否有退化分量；最后画后验概率分布图，人工核对聚类边界合不合理。

我个人最后一个习惯，就是跑完GMM一定把每个样本的后验概率画成直方图看一眼。如果一堆样本的后验概率都堆在0.5以下，说明这些簇之间根本没有真正分开，K可能选大了。这个检查救过我很多次，比单纯看聚类图靠谱得多。GMM这个东西，公式推导看起来标准、代码封装也很成熟，真正拉开差距的还是这些繁琐的细节和检查动作。