高斯混合模型实战：从EM参数估计到聚类与异常检测

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) Σ_{k1}^{K} π_k · N(x | μ_k, Σ_k)其中π_k是第k个分量占所有数据的比例满足π_k ≥ 0且所有π_k之和为1N(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_paramskmeans也是这个思路。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_iter100, tol1e-4): N, D X.shape # 用 KMeans 做初始化均值、协方差、混合权重 kmeans KMeans(n_clustersK, n_init5, random_state0).fit(X) mu kmeans.cluster_centers_ labels kmeans.labels_ pi np.bincount(labels, minlengthK) / 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, axis1) resp np.exp(log_resp - logp[:, None]) # M 步更新混合权重、均值、协方差 Nk resp.sum(axis0) 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, K2) 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_clusters3, n_init10)和GaussianMixture(n_components3, covariance_typefull)分别跑一遍观察簇边界。K-means对斜椭圆簇的处理几乎固定是“从腰部连线切开”因为欧氏距离最近的中心决定了归属两个新簇各占椭圆一半椭圆形状信息完全丢失。GMM则通过协方差矩阵学习到了椭圆的主方向边界会跟着椭圆等密度线走把斜椭圆完整保留下来。两类模型簇之间边界的形态也不同K-means的边界是直线GMM因为各簇协方差不同边界是二次曲线。这个差异在高维数据里会更明显也是为什么GMM对复杂形状簇的还原度比K-means高一个档次。4.3 协方差类型的选择一表讲清楚sklearn的GaussianMixture可以通过covariance_type参数控制协方差矩阵的自由度。如果你不知道这四个选项背后的权衡建议先看下面这张表。covariance_type含义参数量级使用经验full每个分量有独立完整协方差矩阵K·D(D1)/2最灵活数据量足够、维度不高时首选tied所有分量共用一个协方差矩阵D(D1)/2 K·D各簇形状一致但位置不同时用省参数diag每个分量协方差限制为对角阵2·K·D特征之间近似独立时够用训练快spherical每个分量只有半径协方差退化为σ²·IK·(D1)最接近K-means低维小样本可考虑实在拿不准我建议在数据量充足时先跑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调大一点同一个组织就被分成好几块了最后还是降回K2才解决问题。5.2 异常检测用score_samples揪出低概率样本GMM本质上就是密度估计密度越低的区域越可能是异常点。sklearn的GaussianMixture提供了score_samples方法返回每个样本的对数似然。常规做法是用正常数据拟合GMM再对测试集算得分然后取某个分位数作为阈值低于阈值的标记为异常。from sklearn.mixture import GaussianMixture gm GaussianMixture(n_components2, covariance_typefull, random_state0) 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和AICsklearn里直接调用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_componentsk, covariance_typefull, random_stateseed) 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_state0和random_state42跑出来的聚类结果可能差别很大而且两个结果看起来都像模像样。这就是局部最优问题。它不像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去相关。第二层给协方差加ridgesklearn里对应reg_covar参数。默认值1e-6够用但如果数据本身比较脏可以放宽到1e-4甚至1e-3。第三层训练后检查每个分量的责任度之和几乎为0的分量直接删掉或者降K重跑。数值问题不解决后面所有分析都建立在沙子上。6.4 分量退化一个簇只吃到一个样本怎么办和奇异矩阵紧密相关的是分量退化。GMM迭代过程中某个分量可能慢慢把所有责任度都收到一个样本上协方差越缩越小最终趋近于零矩阵。这个分量实际上既不代表一个簇也没有任何解释意义只是过拟合了单个噪声点。最直接的观察方式是打印分量责任度之和responsibilities gm.predict_proba(X) print(responsibilities.sum(axis0))正常情况下每个分量的责任度应该和混合权重一致数量级上没有太离谱的偏差。如果有一项明显低于总体样本量的2%说明这个分量在吃噪声。处理办法包括降K、给协方差加扰动、或者对混合权重设置一个最小值约束。实际项目里我更多时候选择降K因为删除一个“废分量”的代价远低于强行保留它然后给业务解释一个空壳簇。6.5 我的常用流程和一个收尾习惯把前面这些经验串成一个固定流程先标准化数据再PCA或者t-SNE看分布形态确定K的候选范围后用多个seed跑BIC选K建模型时init_paramskmeansn_init设5以上固定random_state拟合完打印责任度之和检查是否有退化分量最后画后验概率分布图人工核对聚类边界合不合理。我个人最后一个习惯就是跑完GMM一定把每个样本的后验概率画成直方图看一眼。如果一堆样本的后验概率都堆在0.5以下说明这些簇之间根本没有真正分开K可能选大了。这个检查救过我很多次比单纯看聚类图靠谱得多。GMM这个东西公式推导看起来标准、代码封装也很成熟真正拉开差距的还是这些繁琐的细节和检查动作。