GMM聚类实战：从KMeans的局限到混合高斯模型与EM算法

做机器学习这几年，我经常遇到一个很尴尬的场景：用KMeans把用户分成三群，兴冲冲拿去给业务看，结果人家指着簇边界问，为什么这两个点明明挨得很近，非要被硬切成不同组？这个问题我早期也解释不清楚，直到系统玩了混合高斯模型（GMM），才真正明白KMeans那种“非黑即白”的硬分配，在处理真实数据时有多大的局限性。

GMM作为无监督学习里最经典的概率生成模型之一，做的事情很直白：假设数据来自若干个高斯分布的叠加，然后通过EM算法反推出每个分布的均值、协方差和权重，再算出每个样本属于每个族的概率。样本不再被强制塞进某一个簇，而是每次给出一组“归属概率”，这就是软聚类。这篇文章我会从特征提取、特征融合、特征降维一路讲到GMM建模和调参，完整走一遍实操流程，用Python和scikit-learn实现，既有原理也有可抄的代码，适合正在学聚类算法、准备机器学习相关面试、或者想把用户分群做得更精细的朋友。

1. 项目整体设计：为什么在无监督学习中选GMM而不是KMeans

1.1 三个典型场景让KMeans露馅

先别急着上代码，先把选型逻辑理清楚。KMeans的核心假设是：每个簇是凸的、大小差不多、并且可以用一个质心代表。这个假设在数据简单的时候很好用，但现实里数据往往长这样：

• 场景A：用户消费行为分群。收入分布天生就带重尾，少数高消费用户离群很远，KMeans用均值当中心，对极端值极其敏感，分群结果容易把真正有价值的高端用户甩成孤零零一个簇，甚至直接并入噪声。
• 场景B：图像像素分割。相机拍出来的物体边界常常有颜色渐变，不同物体的颜色范围会有重叠，KMeans的硬分类会把连续的渐变区域切成锯齿状，后期做掩膜特别痛苦。
• 场景C：风控里的异常检测。正常行为本身就不是一个单峰分布，可能包含“低频正常”和“高频正常”两种模式，KMeans只给一个中心，怎么跑都覆盖不了这种多峰结构。

GMM解决了这三类问题，因为它由多个高斯分量组成，每个分量自带一个协方差矩阵。协方差矩阵让簇可以呈椭圆形、可以有大有小，分量之间允许重叠，而EM算法输出的后验概率让每个样本可以“脚踏多只船”。这就是GMM在聚类任务里和KMeans最本质的区别。

1.2 一个完整的GMM聚类流水线设计

我做一个聚类项目，从来不会只调聚类算法，而是先设计一整套流水线。GMM聚类项目我通常这么搭：

原始数据 → 特征提取/融合 → 标准化 → 降维 → GMM建模 → 评估与调参 → 输出分群结果

每个环节都有存在的理由。特征提取是把日志、图像、文本这些原始素材变成数值向量；特征融合是把手头多路信息拼成一个整体；标准化必须做，因为GMM对特征尺度极其敏感，一旦某个特征量纲巨大，距离计算就被它独占了；降维则是减少噪声和计算量，同时也避免高维空间里协方差矩阵奇异的麻烦；最后才是GMM本身。

这条流水线不是GMM专属，但GMM对前几步更挑剔。原因在后面讲协方差类型的时候会提到：GMM要估计的参数数量随维度平方级增长，维度不控制好，后面全是麻烦。这里先记住一句话：聚类效果差，七成问题出在特征，只有三成出在模型。

2. 前端特征工程：提取、融合、降维的三板斧

2.1 特征提取：从原始数据里捞有用的信息

特征提取这件事，不同数据类型套路完全不同，但核心目标就一个：把原始信息变成能算距离的数值向量。我举几个最常见的场景：

• 表格型数据：原始列本身就是特征，但别急着直接用。连续变量可以衍生统计量，比如均值、标准差、偏度、峰度，以及四分位数，这些能帮助GMM区分“高而稳”和“低而波动”这类隐含模式。
• 时间序列：滑动窗口内部的均值、方差、极差，以及FFT之后的主频、幅值，都是经典特征。做设备故障聚类、行为序列聚类时，这些比直接用原始点更抗噪声。
• 文本：TF-IDF向量或者预训练embedding都可以。
• 图像：颜色直方图、纹理特征、深度特征都常见。

这里我习惯用一个小函数快速验证特征计算：

python复制import numpy as np
from scipy import stats

def extract_features(sequence):
    return np.array([
        np.mean(sequence),
        np.std(sequence),
        stats.skew(sequence),
        stats.kurtosis(sequence),
        np.percentile(sequence, 25),
        np.percentile(sequence, 75),
    ])

这个方法提取出的6个维度，能把一个任意长度的序列压缩成统一长度的向量，方便后续聚类。特征提取的好坏直接决定聚类天花板，GMM不会帮你凭空变出信息，它只是负责把你给它的向量里的结构找出来。

2.2 特征融合：多渠道特征怎么合而不乱

实际项目里很少有人只用一路特征。用户画像可能是“基础属性 + 行为统计 + 消费记录”三路拼在一起，文本分类可能是“TF-IDF + 词向量”一起用。最常用的融合方式就是直接拼接，一行代码：

python复制X = np.hstack([feature_1, feature_2, feature_3])

直接拼接虽然简单，但有两个坑必须避开。

第一是量纲差异。一个特征是年龄（20~60），另一个特征是月消费频次（0~10000），直接拼在一起，GMM在算距离时基本被频次主导，年龄维度形同虚设。解决办法是拼接之前先对每一路特征单独做标准化。

第二是维度爆炸。每多一路特征，GMM的参数空间就涨一截。我处理这类问题时，会先看各路特征的方差贡献，如果某一路特征全是近常数的，就果断扔掉。特征选择在无监督里没有标签可用，我一般用方差过滤加相关性去重，相关性超过0.9的两个特征只留一个。

2.3 特征降维：PCA为什么是GMM的最佳拍档

特征融合完，通常紧接着就是降维。我几乎每次都会告诉读者：GMM和KMeans在降维这件事上的态度完全不一样。KMeans在原始高维空间也能跑，只是慢；GMM在高维空间是容易直接崩掉，因为每个簇要估计一个d×d的协方差矩阵，参数数量是d(d+1)/2级别。30维的时候，每个簇的协方差矩阵就有465个参数要估计，没有足够样本量根本学不干净，结果就是协方差矩阵奇异，模型直接报错。

PCA的作用是找到方差最大的若干正交方向，把数据投影到低维空间。在GMM流水线里，我一般保留累计方差贡献率90%到95%的主成分：

python复制from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

pca = PCA(n_components=0.95)
X_pca = pca.fit_transform(X_scaled)

这里有一个容易忽略的细节：PCA之前必须先标准化，否则方差最大的方向会被量纲大的特征绑架，降维结果不反映真实结构。

降维之后不只是计算变快，聚类效果往往也变好。高维空间里距离度量会趋向均匀化，GMM对概率密度的估计也容易失真；把有效信息压缩到前几个主成分后，簇与簇之间的区分度反而更明显。

3. 核心原理：混合高斯模型的底层逻辑

3.1 用一枚硬币的故事读懂EM算法

先讲一个经典例子。假设桌上有一枚硬币A和一枚硬币B，它们正面朝上的概率分别是多少我们不知道，只知道每次随机拿一枚硬币抛10次，记录了正面次数，但不知道某次用的具体是哪枚硬币。这个“不知道是哪枚硬币”的信息，就是隐变量。EM算法的思路是：先随便猜一组硬币参数，然后根据当前参数算每次实验更可能来自A还是B，这是E步；再基于这个归属概率重新估算硬币参数，这是M步。反复迭代，参数会收敛到一个稳定值。

GMM的EM训练完全就是这枚硬币故事的连续版本。隐变量是“每个样本来自哪个高斯分量”，我们观察不到，只能从数据里反推。

• E步：根据当前的均值、协方差、权重，计算每个样本属于每个分量的后验概率：

[
\gamma(z_{ik}) = \frac{\pi_k \mathcal{N}(x_i \mid \mu_k, \Sigma_k)}{\sum_{j=1}^{K} \pi_j \mathcal{N}(x_i \mid \mu_j, \Sigma_j)}
]

• M步：用这些后验概率作为软权重，重新加权计算每个分量的均值、协方差和权重：

[
\mu_k = \frac{\sum_i \gamma(z_{ik}) x_i}{\sum_i \gamma(z_{ik})}
]

[
\Sigma_k = \frac{\sum_i \gamma(z_{ik})(x_i - \mu_k)(x_i - \mu_k)^T}{\sum_i \gamma(z_{ik})}
]

[
\pi_k = \frac{\sum_i \gamma(z_{ik})}{N}
]

E步和M步轮流执行，直到对数似然函数的变化量小于某个阈值，比如默认的1e-3。有一个很重要的点必须明确：EM算法保证每次迭代都会提升对数似然的下界，但不保证收敛到全局最优，所以训练时多跑几次不同的初始值很有必要。这是GMM调参的第一大坑，后面章节细说。

3.2 协方差类型的选择：full、tied、diag、spherical怎么选

scikit-learn的GaussianMixture提供了四种协方差类型，这个参数几乎决定了模型复杂度和训练成败，我用一张表来对照：

选择原则其实很实际：数据量少于特征数10倍时，我绝不碰full，出奇异协方差矩阵是大概率事件，每一步都像踩雷。cluster的边界不规整时，优先diag，因为diag只需要估计每个维度的方差，稳定性好很多，训练速度快一到两个量级。只有当样本量充足、并且确认不同簇的形状差别很大时，才考虑full。

3.3 聚类数量怎么定：BIC、AIC和业务常识

KMeans定k可以看肘部图，GMM虽然有手肘法可以用，但更规范的做法是看BIC或AIC。BIC的公式是：

[
BIC = -2 \times \log L + p \times \ln(n)
]

其中logL是模型在数据上的对数似然，p是模型参数个数，n是样本量。BIC越小，代表模型在“拟合优度”和“复杂度”之间找到了更好的平衡，它会给参数过多的模型加惩罚。

实操时遍历候选的n_components：

python复制from sklearn.mixture import GaussianMixture

bics = []
aics = []
K_range = range(1, 10)

for k in K_range:
    gmm = GaussianMixture(n_components=k, covariance_type='diag',
                          random_state=42, n_init=5)
    gmm.fit(X_pca)
    bics.append(gmm.bic(X_pca))
    aics.append(gmm.aic(X_pca))

best_k = np.argmin(bics) + 1

BIC曲线通常会在某个位置出现拐点，之后下降变缓甚至回升，这个拐点就是推荐簇数。但我必须强调一句：BIC算出来的K只是一个技术参考，最后还是要回到业务上验证。我曾经遇到一个数据集，BIC强烈推荐聚类成8簇，但其中两簇的业务画像几乎一样，合并掉效果更好。统计指标帮我们缩小范围，业务常识负责做最终决策。

4. 实操：用Python跑通一个GMM聚类项目

4.1 环境准备和合成数据生成

先把环境列清楚。我用的是Python 3.10，主要依赖库版本：numpy 1.26、scikit-learn 1.5、matplotlib 3.8、scipy 1.12。这些库直接用pip安装即可，不需要额外配置。

为了演示，我自己合成一份三类二维数据，故意让两类在局部重叠，模拟真实场景：

python复制import numpy as np

np.random.seed(42)
n = 500

# 三个高斯分布，第二个和第三个部分重叠
X1 = np.random.multivariate_normal([0, 0], [[1.0, 0.5], [0.5, 1.0]], n)
X2 = np.random.multivariate_normal([3, 3], [[1.2, -0.3], [-0.3, 0.8]], n)
X3 = np.random.multivariate_normal([2.5, -1.0], [[0.6, 0.0], [0.0, 0.6]], n)

X = np.vstack([X1, X2, X3])
y_true = np.array([0] * n + [1] * n + [2] * n)

生成的数据有真实标签，可以方便后面用ARI评估聚类质量。三个类的形状差异明显，有圆有椭圆，正适合展示full协方差的威力。

4.2 完整代码：预处理、降维、GMM拟合

数据准备好之后，走一遍完整流水线：

python复制from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.mixture import GaussianMixture

# 标准化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# PCA降维
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)

# 用BIC选聚类数
bics = []
K_range = range(1, 8)
for k in K_range:
    gmm = GaussianMixture(n_components=k, covariance_type='full',
                          random_state=0, n_init=10)
    gmm.fit(X_pca)
    bics.append(gmm.bic(X_pca))

best_k = np.argmin(bics) + 1
print("BIC最优聚类数:", best_k)

# 最终模型
gmm = GaussianMixture(n_components=best_k, covariance_type='full',
                      random_state=0, n_init=10)
gmm.fit(X_pca)

labels = gmm.predict(X_pca)
proba = gmm.predict_proba(X_pca)

这份代码在二维数据上会得到best_k等于3，因为数据本身就是三簇生成，BIC在3处取到最小值。predict返回的是每个样本概率最大的那个簇，predict_proba输出的是n_samples行k列的归属概率矩阵，每行加起来等于1，这就是软聚类的输出。

4.3 可视化技巧：如何画出椭圆置信区间

聚类结果只是看散点图不够，我会把每个高斯簇画成置信椭圆。椭圆的中心和方向由分量的均值和协方差矩阵特征向量决定，长短轴与特征值平方根相关。画95%置信椭圆，可以这样实现：

python复制import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from scipy.stats import chi2

def draw_gmm_ellipse(ax, mean, cov, n_std=2.0):
    eigvals, eigvecs = np.linalg.eigh(cov)
    order = eigvals.argsort()[::-1]
    eigvals, eigvecs = eigvals[order], eigvecs[:, order]
    angle = np.degrees(np.arctan2(eigvecs[1, 0], eigvecs[0, 0]))
    scale = chi2.ppf(0.95, df=2)
    width, height = 2 * np.sqrt(scale * eigvals)
    ellipse = Ellipse(xy=mean, width=width, height=height,
                      angle=angle, alpha=0.3, edgecolor='black', linewidth=1.5)
    ax.add_patch(ellipse)

plt.figure(figsize=(8, 6))
for i in range(best_k):
    mask = labels == i
    plt.scatter(X_pca[mask, 0], X_pca[mask, 1], s=10, alpha=0.6)
    draw_gmm_ellipse(plt.gca(), gmm.means_[i], gmm.covariances_[i])
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.title("GMM Clustering with Confidence Ellipses")
plt.show()

chi2.ppf(0.95, df=2)等于5.991，它把椭圆缩放到了“包含95%概率质量”的形状。实际画出来之后会发现，第三簇的椭圆细长，第一第二簇的椭圆有倾斜，这种形状信息是KMeans完全表达不出来的。

4.4 评估与调参：从轮廓系数到ARI

无监督聚类没有万能评估指标，但有三个维度可以交叉验证。有真实标签时用调整兰德指数（ARI），它的取值在-1到1之间，1表示完全一致，0表示随机分配。轮廓系数不需要真实标签，评估的是簇内紧致和簇间分离的平衡，但它在高维空间里不可靠，所以我只在降维后的空间计算它。第三个是模型本身的对数似然，GMM训练完直接看gmm.score(X_pca)即可，这个值反映的是概率模型对数据的拟合程度，不是聚类质量，但可以用来挑超参数。

我的实操习惯是这样：用ARI确认这个方法在已知答案的数据上没有跑偏，用轮廓系数辅助判断簇内是否紧凑，最后用BIC在多个候选配置之间做决定。

python复制from sklearn.metrics import adjusted_rand_score, silhouette_score

ari = adjusted_rand_score(y_true, labels)
sil = silhouette_score(X_pca, labels)
print("ARI:", ari, "Silhouette:", sil)

在我的经验里，当数据有三个高斯分量且重叠不严重时，ARI一般能超过0.95。如果ARI偏低，问题多半出在预处理而不是GMM本身。

5. 实战中常见的坑与排查技巧

5.1 初始化不稳定：一个结果跑三次不一样

遇到过太多次了：数据没变，代码没变，只是random_state变了，聚类结果完全两样。原因在于EM算法本质上是坐标上升的变体，对初值极其敏感，很容易陷进局部最优。

解决方案很直接：把n_init设成10或者更大，让模型从多个初始点出发，保留对数似然最高的那次结果。另一个技巧是把init='kmeans'，先用KMeans的结果作为GMM的初始化，收敛速度和稳定性都有提升。这个参数几乎是我写GMM的标配。

python复制gmm = GaussianMixture(n_components=3, covariance_type='full',
                      init='kmeans', n_init=10, random_state=42)

5.2 协方差矩阵奇异：Warning一大堆

训练时经常看到类似“ConvergenceWarning: The number of distinct clusters is less than the number of components”这样的提示，核心原因是某个分量的协方差矩阵变得奇异，通常是维度太高、样本太少、或者特征之间高度共线。

我的排查顺序是：先检查特征是否做了标准化和PCA，确认特征数量远小于样本数量；其次增大reg_covar，给协方差矩阵对角线加一个小的正则项，比如1e-6到1e-2；最后实在不行就把covariance_type从full降成diag或spherical，牺牲一点形状灵活性，换来数值稳定性。

5.3 聚类数定了对不上业务

BIC是最小化信息准则的技术指标，业务要的是能解释的分群。如果BIC建议8簇，但两个簇的画像几乎重合，这时候不要硬用8簇。我的做法是先用BIC选一个大致的范围，然后跑完看每个簇的weights_，也就是混合权重。如果某些簇的权重低于0.05，说明它只是一个稀疏的边缘分量，可以考虑合并到最近的簇里。

还有另一种情况：业务明确只需要3类，但数据天然是5簇。这时候可以先把GMM当密度估计器用，用多余的分量拟合背景分布，再把业务关注的几个大类拆出来，而不是反过来改数据。

5.4 特征尺度对GMM的影响

GMM里的“距离感”来自高斯分布里的马氏距离，它综合了协方差矩阵的信息。但这个协方差矩阵是被特征本身的尺度强烈影响的，如果一个特征范围是0到10000，另一个是0到1，第一个特征的方差会压过第二个，导致模型以为第二个特征没有区分度。

在聚类之前对每个特征做零均值单位方差标准化，这一步不能省。做PCA之前也必须先标准化，不然PCA找出来的主成分方向会被大方差特征带偏。这两个顺序很关键：先标准化，再降维，再进GMM。

5.5 GMM当分类器用：从无监督走向半监督

GMM不只是聚类工具，它天然也可以当生成式分类器用。如果是完全有监督场景，可以对每个类别单独拟合一个GMM，新样本进来时分别计算它在各类别GMM下的对数似然，取最大那个作为预测类别。这本质上是朴素贝叶斯的进阶版，只不过类条件概率密度用了混合高斯而不是简单的高斯分布。

另一种更常见的使用方式是半监督：先在一部分有标注的数据上初始化GMM分量，然后把无标注数据也丢进去参与EM迭代，让每个高斯分量对应一个类别。这种方式在标签稀缺的场景下特别实用。我在处理带噪设备采集数据时，经常用这个方法先分群，再把业务上确认的少数样本映射到簇上，完成样本扩展。

最后分享一点个人体会。我在实际项目里用过很多次GMM，最大的感受是：它比KMeans多给了一个“概率输出层”，这个输出层在生产环境里特别值钱。比如做用户精细化运营，KMeans只能告诉你用户属于A群还是B群，而GMM告诉你这个用户有0.85的概率属于高价值群，0.1的概率属于流失风险群，运营可以直接拿这个概率做阈值策略。概率大于0.8走高优先级处理，0.5到0.8走普通券包，小于0.5再观察，这样的分群落地方式比硬标签灵活得多。

另外想提醒一句：GMM不是无所不能，它对数据分布有很强的假设——簇必须是高斯形状。如果你的数据分布是U型的、环形的，或者类似月牙形状，GMM会表现得非常别扭，这时候换成DBSCAN或者谱聚类反而更合适。遇到这种数据，先画个二维分布图看一眼，再决定选什么模型。聚类算法没有银弹，选型靠的是对数据结构的理解和试错，这个基本功比会调参重要得多。