做聚类分析的人多半遇到过这种场面：数据点画在散点图上，明明是两团斜着叠在一起的椭圆，KMeans 却硬生生从中间切了一刀，把同一个簇劈成两半，又把两个簇的边缘拼成一个。我最早做设备振动信号分群的时候就被这个坑绊了整整两天，后来换成 sklearn 的 GaussianMixture，问题才真正解决。这个模型背后跑的是 EM 算法，它不再回答"这个点属于哪一类"，而是回答"这个点由每个高斯分量生成的概率分别是多少"。这篇内容会从 GMM 到底在拟合什么讲起，把 EM 算法在 sklearn 源码里的执行轨迹拆开，逐个过一遍 GaussianMixture 构造函数的参数，然后用三个任务完整跑通，最后把我踩过的收敛失败、协方差塌陷、标签漂移这些坑原样摆出来。只要你会写 numpy 和 sklearn 的基础代码，就能跟着复现。

1. 为什么簇叠在一起时KMeans会失效：GMM要解决的三个具体问题

1.1 硬分配与软分配的差别到底体现在哪

KMeans 的输出是一个整数标签，样本点到簇心的距离决定了它归谁，这叫硬分配。问题出在"边界"这两个字上——现实数据里两个簇的过渡区域往往是连续渐变的，一个样本可能 60% 像 A 簇、40% 像 B 簇，KMeans 会毫不犹豫地把它判给 A，然后在下一轮迭代里把这个点算进 A 的质心计算，质心被拉偏，边界继续歪，误差就这么累积下去。

GaussianMixture 干的事情完全不同。它假设整份数据是由 K 个多元高斯分布按不同比例混合生成的，每个高斯有自己的均值向量和协方差矩阵，还有一个混合权重。对于每个样本，模型算的是它在每个高斯分量下的概率密度，再乘上分量权重，归一化之后得到一组和为 1 的责任度（responsibility）。这组数字就是软分配，它保留了"这个点有多大程度属于每个簇"的信息。

这个差别在业务上非常实在。我做用户行为分群时，一个刚注册三天、行为特征处于两类用户中间的用户，硬分配会给他贴一个确定的标签，运营策略只能二选一；软分配能告诉我他有 0.55 的倾向属于高活跃群体、0.45 倾向属于观望群体，那就两条策略都触达一次，或者干脆放进一个单独的观察池。这种"模糊地带"的处理能力，是 GMM 相对 KMeans 最直接的价值。

1.2 协方差矩阵带来的"形状自由度"

KMeans 隐含的假设是每个簇都是一个各向同性的球，所有方向上的方差一样大。这个假设在真实数据里经常不成立。传感器采集的温度和湿度两个维度，相关性很强，数据点分布是一条斜着的长条；KMeans 只能用圆形去套，套不住就只能切分。

GMM 的每个分量带一个协方差矩阵，它同时编码了两件事：各个维度上的方差大小（决定椭圆的长短轴），以及维度之间的相关性（决定椭圆倾斜的角度）。一个二维的 full 协方差矩阵有 3 个自由参数，能表达任意方向的椭圆；三维是 6 个；d 维是 d(d+1)/2 个。这种表达能力带来的直接后果是，即使两个簇的均值离得很近，只要它们的朝向和扁长程度不一样，GMM 依然能把它们分开，而 KMeans 在这种情况下基本无解。

代价也很明确：参数多了就容易过拟合。样本量不够的时候，full 协方差会去拟合噪声，把某个簇收得很紧，训练集上的似然很好看，换一批数据就崩。这就是后面 reg_covar 参数存在的理由，也是 covariance_type 需要认真选的原因。

1.3 混合权重：允许簇大小不一样

KMeans 在目标函数里对每个点一视同仁，簇的大小基本由几何距离决定，很难出现"一个大簇 + 一个小簇"的自然结构。GMM 引入了混合权重 π，它满足所有分量权重之和为 1，物理含义是"随机抽一个样本，它来自第 k 个分量的先验概率"。

这个设计对不平衡数据特别友好。比如做线上的异常订单识别，正常订单可能占 97%，异常占 3%，如果强行用 KMeans 分两类，那 3% 的小簇会被大簇吞掉。GMM 的三个分量权重可能是 0.62、0.35、0.03，最后那个分量虽然只占 3%，但它在自己的均值附近依然有独立的方差结构，模型会专门去描述这一小撮数据的形状。

权重的另一个用处是作为业务指标。跑完拟合之后看gmm.weights_，你立刻就知道这份数据里各个群体的规模比例，这个数字比簇内样本计数更平滑，因为它考虑的是软分配下的期望样本数，不会因为边界上的少数点抖动。

2. EM算法在sklearn源码里是怎么一步步跑起来的

2.1 E步：责任度矩阵的精确计算

EM 算法分两半，E 步负责在给定当前参数的情况下，算出每个样本对每个分量的责任度。数学上写出来是这样的：先算第 i 个样本在第 k 个分量下的加权密度，即 π_k 乘以多元高斯密度函数 N(x_i | μ_k, Σ_k)，然后把这个值除以该样本在所有分量上的加权密度之和。

sklearn 的实现里做了对数化处理，因为直接算密度函数在高维下极容易下溢成 0。源码_estimate_weighted_log_prob计算的是对数加权概率，E 步实际上是在对数域做 softmax，这一步的数值稳定性直接影响模型能不能跑出正确结果。你如果自己手写一遍就会明白为什么不能直接算乘积：

import numpy as np from scipy.stats import multivariate_normal def manual_e_step(gmm, X): n_samples, n_features = X.shape k = gmm.n_components weighted = np.zeros((n_samples, k)) for j in range(k): # 注意：只有 covariance_type='full' 时 covariances_[j] 才是完整矩阵 weighted[:, j] = gmm.weights_[j] * multivariate_normal( mean=gmm.means_[j], cov=gmm.covariances_[j] ).pdf(X) return weighted / weighted.sum(axis=1, keepdims=True) # 与 sklearn 的 predict_proba 对比 resp_manual = manual_e_step(gm, X) resp_sklearn = gm.predict_proba(X) print(np.abs(resp_manual - resp_sklearn).max()) # 正常在 1e-12 量级

这段代码我很推荐每个人都跑一遍，哪怕只为了建立信心。当你看到自己手写的 E 步和predict_proba的输出差在 1e-12 以内，你对这个模型的信任度会和只看文档完全不一样。顺便提醒一个细节：covariances_的形状随covariance_type变化，full是 (K, d, d)，tied是 (d, d)，diag是 (K, d)，spherical是 (K,)，上面这段代码只适用于full。

2.2 M步：三个参数各自的更新公式与推导思路

M 步拿着 E 步算出的责任度矩阵去重新估计参数，三个参数的更新都带一个共同的前置量：第 k 个分量的有效样本数 N_k，它等于所有样本对第 k 个分量的责任度之和。有了这个 N_k，后面就都是加权平均。

均值更新最直观，就是所有样本按责任度加权后的重心：μ_k = (1/N_k) Σ γ_ik x_i。协方差更新是把每个样本到新均值的偏差做外积，再用责任度加权平均：Σ_k = (1/N_k) Σ γ_ik (x_i - μ_k)(x_i - μ_k)^T。权重更新最简单，π_k = N_k / N，也就是这个分量的有效样本数占总样本数的比例。

理解这三个公式的价值在于，你能预判模型在什么情况下会出问题。比如某个分量只"抢"到极少数样本（N_k 很小），它的协方差估计就会非常不稳定，可能变成一个极度扁平的椭圆，甚至接近奇异矩阵。这时候对数似然反而会飙升，因为高斯密度在极窄的峰上能取到很大的值——这就是典型的协方差塌陷，后面第 5 章会专门讲怎么识别和处理。

还有一条性质值得记住：EM 每一轮迭代都不会让对数似然下降。这是由 Jensen 不等式保证的，也是判断你的实现有没有写错的最快方法。如果用 sklearn 看到lower_bound_序列不是单调递增的，那基本可以确定数值上出了问题，通常是精度矩阵求逆那一步的数值误差。

2.3 lower_bound_ 与 converged_：判断收敛的两个信号

跑完fit之后，模型对象上有两个属性必须第一时间看，converged_和n_iter_。converged_是布尔值，表示是否在max_iter之内满足了停止条件；n_iter_是实际迭代轮数。如果converged_是 False 而n_iter_正好等于你设的max_iter，说明模型是被强制截断的，参数还没稳定，这时候拿去做预测风险很大。

lower_bound_保存的是最后一次迭代的对数似然下界。它在调试中有两个用处：一是横向对比不同n_components的模型，似然越高说明拟合越好（但这个不能直接用来选模型，因为分量越多似然一定越高，得用 BIC）；二是看收敛速度，如果n_iter_只有三四次就 converged 了，通常意味着初始化的 KMeans 已经给出了很好的起点，或者数据本身结构很简单。

停止条件的判定逻辑是相邻两次迭代的lower_bound_差值小于tol，默认tol=1e-3。这个阈值对数量级在 -10^4 的对数似然来说已经相当宽松了，实践中我一般会调到 1e-4 甚至 1e-5，代价是迭代轮数增加，但换来的是参数更稳。数据量在十万级以下的时候，这点额外开销完全可以接受。

3. GaussianMixture构造函数的参数逐个过一遍

3.1 n_components 与 covariance_type 的组合代价

n_components是分量个数，也就是你期望的簇数量，这个参数没有默认值的"正确答案"，必须结合业务和 BIC 一起定。真正容易被低估的是covariance_type，它决定了参数量级，直接影响到需要多少样本、跑得多快、会不会过拟合。

四种取值的参数量对比如下，其中 K 是分量数、d 是特征维度：

举个具体数字，d=10、K=5 的时候，full 需要 275 个协方差参数，tied 只要 55 个，diag 是 50 个，spherical 只有 5 个。参数量的差距在两个数量级以上，这就是为什么高维数据上直接用 full 经常训不动或者训出一堆奇异矩阵。

我的经验法则是：每个分量的有效样本数至少要达到 10d 到 20d 这个量级，才考虑用 full。如果样本达不到，先从 diag 起步，看 BIC 和实际业务效果再决定要不要升到 full。tied 这个选项很多人忽略了，它在"同一批设备的不同工况"这类场景里特别好用——不同工况下信号的波动形态是一样的，只是中心位置不同，共享协方差正好对应这个物理假设。

3.2 reg_covar：那个默认1e-6到底在兜什么底

reg_covar默认是 1e-6，作用是在协方差矩阵的对角线上加一个非负小量，保证矩阵正定、可求逆。听起来像个不起眼的数值技巧，实际上它是 GMM 不崩溃的最后一道防线。

什么时候会用到它？当某个分量分到的样本极少，或者这些样本在某个维度上方差接近 0，协方差矩阵就会退化，求逆的时候直接抛LinAlgError。加了reg_covar之后，对角线被抬起来，矩阵被"撑"成严格正定，计算能继续下去。这个机制在 sklearn 内部是无条件生效的，你不能关掉它，只能调大小。

调多大合适？这里有个很多人踩的坑：reg_covar是一个绝对数值，不是相对比例，所以它对特征量纲极其敏感。如果你的特征标准化到 0 附近、方差为 1，那 1e-6 就是个很温和的正则；但如果你的特征量级在 10^6（比如某些计数型指标），1e-6 相对于量级来说微不足道，等于没加，该塌陷还是塌陷；反过来，如果特征方差在 10^-8 量级（比如某些归一化后的概率值），1e-6 会把所有结构抹平，模型退化成一个球。

我的做法是先做标准化，再按需要把reg_covar调到 1e-3 到 1e-1 之间试试。看到一个分量塌陷成极窄的峰，第一反应就是加大reg_covar，通常加到 1e-3 就能解决大部分问题。

3.3 init_params 与 n_init：初始化决定你会不会掉进局部最优

init_params默认是'kmeans'，也就是先用 KMeans 跑一遍拿到簇心，把簇心当作各分量的初始均值，然后用全体数据的协方差初始化各个分量。另外还有'k-means++'、'random'和'random_from_data'三个选项。

默认用 KMeans 初始化是有道理的：KMeans 虽然假设球形，但它给出的簇心通常比随机初始化的点靠谱得多，能显著降低 EM 掉进糟糕局部最优的概率。'random_from_data'从数据里随机挑 K 个点当初始均值，在数据有明确多峰结构时效果也不错，但方差大。'random'我基本不用，除非在做算法对比实验。

n_init控制的是"整套初始化 + EM 迭代"重复多少遍，最后保留对数似然最高的那次结果。默认值是 1，这一点经常被误解，很多人以为跟某些聚类算法一样默认跑 10 次。实际只有 1 次，意味着你的结果完全取决于那一次初始化的质量。

我在任何正式的项目里都会把n_init设成 5 到 10。代价是训练时间线性增长，换来的是结果稳定性大幅提升。判断要不要加大n_init有个简单办法：同一份数据、同一个random_state，连续跑几次不同随机种子的模型，看lower_bound_的波动范围。如果波动超过 1%，说明初始化很敏感，n_init就得往上加。

还有三个手动初始化参数：means_init、weights_init、precisions_init。这里有个大坑——precisions_init要求的是精度矩阵，也就是协方差矩阵的逆，不是协方差本身。我第一次用的时候传了协方差进去，模型不报错，但结果完全不对，查了很久才发现。而且它的形状同样随covariance_type变化，用之前务必对着文档确认。

3.4 tol 与 max_iter 的配合，以及 warm_start 的正确用法

max_iter默认 100，tol默认 1e-3。这两个参数是配套的：tol判得松，max_iter就很少被触顶；tol判得严，max_iter就可能不够用。

我建议的配置思路是：先把tol调到 1e-4，max_iter保持 100 跑一次，看converged_和n_iter_。如果n_iter_在 50 以内就收敛了，说明数据规模不大，这个配置很稳；如果n_iter_逼近 100 还没收敛，先别急着把max_iter调到 500，因为迭代 100 轮还不收敛通常意味着数据本身存在大量重叠结构，或者n_components设多了，加迭代次数只是在给一个不合理的模型硬凑解。

warm_start默认 False。设成 True 之后，再次调用fit会从上次的结果继续迭代替换重新初始化，配合手动传入means_init和weights_init可以做增量式的参数调整。我一般只在做超参网格搜索时用它，比如固定初始化，只改reg_covar，这样可以保证不同配置之间可比。日常训练别开这个，容易因为状态残留导致结果难以复现。

另外verbose和verbose_interval在排查不收敛问题时很有用，设verbose=2, verbose_interval=1，每轮迭代的对数似然都会打出来，你能直观看到似然是不是单调上升、有没有停滞在某处反复横跳。

4. 三个真实任务的完整跑通过程：密度估计、软聚类、异常检测

4.1 密度估计：拟合出来的分布能不能还原真实形状

先从一个可控的实验开始。造三团二维数据，均值分别是 (0,0)、(4,4)、(-4,4)，协方差各不相同，然后让模型去学。这个实验的意义在于，你知道真值，可以直接对比学出来的参数和真实参数。

import numpy as np from sklearn.mixture import GaussianMixture rng = np.random.RandomState(0) true_means = np.array([[0., 0.], [4., 4.], [-4., 4.]]) true_covs = np.array([ [[1.0, 0.6], [0.6, 1.0]], [[0.5, 0.0], [0.0, 2.0]], [[2.0, -0.7], [-0.7, 0.5]], ]) X = np.vstack([ rng.multivariate_normal(true_means[k], true_covs[k], size=400) for k in range(3) ]) gm = GaussianMixture( n_components=3, covariance_type='full', n_init=10, tol=1e-5, max_iter=300, random_state=42 ).fit(X) print("converged:", gm.converged_, "iters:", gm.n_iter_) print("weights:", np.round(gm.weights_, 3)) print("means:\n", np.round(gm.means_, 3)) print("avg log-lik:", round(gm.score(X), 4))

跑完之后大概率会看到权重接近 0.333 的三等分，均值基本对上真值。协方差矩阵因为分量编号是任意的，需要自己按均值距离做一次匹配再比较。这一步很值得做，因为它是验证整条链路正确的唯一办法。

score(X)返回的是平均对数似然，注意是平均，不是总和。这个区别在对比不同大小的数据集时会坑人，如果你想比总和，用score_samples(X).sum()。做可视化的时候，score_samples逐点输出的分数正好可以画等高线，把三团数据的密度分布原样画出来。

还有一点：score_samples返回的是对数密度，np.exp之后才是密度值，但不一定积分为 1，因为它是相对于经验分布的估计。别拿它直接当概率用。

4.2 软聚类：predict_proba 输出的责任度怎么用

predict返回的是predict_proba每行的 argmax，这两者的关系一定要清楚。真正有信息量的是predict_proba的输出矩阵，形状是 (n_samples, n_components)。

拿到这个矩阵之后，我通常做三件事。第一件是筛出"高置信样本"，也就是某一行最大值超过 0.9 的样本，这些是各簇的核心成员，适合用来做簇的特征画像，因为它们不会被边界模糊污染。第二件是筛出"模糊样本"，最大值低于 0.6 的，这些样本单独成池，做人工抽检或者走单独的规则。第三件是把责任度当权重做加权统计，比如算某个簇的平均用户价值，用责任度加权比用硬标签更稳。

一个必须提醒的事实是：责任度不是置信度。0.5 对 0.5 不代表"模型很没把握"这种事需要警惕，它就是数学上两个高斯在该点的密度相等，很正常。反过来，0.99 对 0.01 也不代表这个分配一定符合业务语义，它可能只是那个分量收得特别紧，密度峰值很高。解读的时候要结合协方差看，不能光看数字。

还有标签编号的问题。GMM 的predict给出的 0、1、2 是内部编号，跟你的业务类别没有任何对应关系，而且换一个random_state编号就可能整体重排。生产代码里千万不要硬编码"标签 0 是高价值用户"，必须通过均值、权重这些稳定特征去动态映射。

4.3 异常检测：score_samples 分位数阈值法

GMM 是生成模型，它对"正常数据长什么样"建了模，那么密度低的地方就是异常。这个思路比很多专门的异常检测算法更直观，而且解释性强——你可以告诉业务方"这个样本落在所有已知分布的极低密度区"。

scores = gm.score_samples(X) # 每个样本的对数密度，越大越正常 thr = np.percentile(scores, 1.0) # 取最低的 1% 作为异常 flag = scores < thr print("异常样本数:", flag.sum(), "阈值:", round(thr, 4))

分位数的选择直接决定召率和误报率，没有普适值。我一般的做法是先在历史数据上画出scores的直方图，看左尾的形状。如果左尾有一个明显分离的小峰，那阈值就取两峰之间的谷；如果左尾是平滑延伸的，那就只能按业务能承受的误报率倒推，比如运维团队每天只能处理 20 个告警，那就取分位数让日均告警数落在 20 附近。

有个陷阱需要注意：GMM 对训练数据里的异常值是敏感的，因为它会专门分出一个分量去描述那些离群点。所以做异常检测时，要么先做一轮粗筛把明显离谱的数据剔掉，要么让n_components稍微大一点，多出来的分量自然会去承接异常结构。跑完之后检查weights_，如果某个分量权重特别小（比如 0.005）而且均值离其他分量很远，那它大概率就是在描述异常，这时候你可以直接把"属于该分量概率高"作为异常判据，比全局密度阈值更精准。

4.4 用 sample 做数据增强与压力测试

sample(n_samples)从拟合好的模型里采样，这是生成模型独有的能力。较新版本的 sklearn 返回一个二元组 (X, labels)，早期版本只返回 X，用之前确认一下版本，避免解包报错。

这个能力有两个正经用途。第一个是做小样本类别的数据增强：某个类别的真实样本太少，训不动下游分类器，那就用 GMM 拟合这个类的分布，采样出补充样本。这里要诚实——采出来的样本只是对分布的重采样，不带来新信息，对提升模型泛化能力有限，但对稳定某些对样本量敏感的算法（比如小样本下的树模型）确实有帮助。

第二个用途我觉得更有价值：压力测试。把一个离线训练好的 GMM 拿来采样出十万条数据，灌进你的线上推理链路，看看有没有性能瓶颈、有没有数值异常。这比等真实流量慢慢涨要快得多，而且采样数据天然覆盖了模型认为"可能存在"的区域，包括真实数据里还没出现过的组合。

采样的时候注意random_state要设，否则每次结果不一样，没法做回归测试。另外采出来的数据是无标签的联合分布采样，它不会保持原始数据的样本间相关性，如果你的下游逻辑依赖样本顺序或时间序列结构，这个方法就不适用。

5. 踩过的坑：不收敛、协方差塌陷与结果漂移

5.1 "Converged: False" 出现时的排查顺序

看到converged_是 False，别急着调max_iter，按下面这个顺序查，能省很多时间。

第一步，看n_iter_是不是正好等于max_iter。如果有差距，说明是别的原因提前结束了。第二步，把verbose=2打开，观察对数似然曲线。如果它在稳步上升只是没到阈值，那就是tol太严或者max_iter不够，适当放宽或者加迭代；如果它在某个值附近来回震荡，那问题在数值稳定性，通常是reg_covar太小或者特征量纲差异太大。第三步，检查特征。有没有常量列？有没有极端离群值？常量列的方差是 0，直接导致协方差奇异；极端离群值会把某个分量的均值拽走，让 EM 在两套解之间反复。

我遇到过一次特别隐蔽的情况：数据里有几列是整数计数，量级在 10^5，另外几列是比例，量级在 10^-2。协方差矩阵的条件数直接爆掉，EM 每轮都在两个解之间跳，永远不收敛。解决办法就一句话——标准化。做完StandardScaler之后同一个模型 12 轮就收敛了。

还有一种情况是n_components设得比真实结构多太多。比如数据其实只有 2 个簇，你设了 10，多出来的 8 个分量会去争抢边界上的残余变异，互相推挤，收敛极慢甚至不收敛。这时候 BIC 会明确告诉你答案是 2。

5.2 单样本簇导致的协方差奇异

这个坑的表现很典型：模型能跑完，converged_是 True，但weights_里有一个分量权重接近 0，对应的协方差矩阵行列式接近 0，predict_proba在这块区域输出极端值（0 或者 1）。

原因就是 M 步里那个 N_k。如果某个分量在 E 步几乎没抢到样本，N_k 极小，它估计出的协方差就完全由那一两个点决定，方向任意、尺度极小。下一轮 E 步，这个极窄的高斯在它均值附近给出巨大的密度，又会把附近的点全部吸过来，然后又散掉——如此反复。

处理办法有三个层次。最直接的是加大reg_covar，从 1e-6 提到 1e-3，把那个窄峰撑开，让它没法独占样本。第二个是减少n_components，从根上避免多余分量的存在。第三个是后处理：跑完之后检查weights_，把权重低于某个阈值（比如 1/(10K)）的分量标记出来，如果它的均值也远离其他分量，那它就是个伪分量，要么合并要么丢弃。

我更倾向于在训练前就把这件事防住：先用 BIC 定一个合理的n_components，再用reg_covar=1e-3起步，跑完检查权重分布是否合理。三步下来，塌陷基本不会出现。

5.3 每次运行结果不一样：随机种子与标签置换

同一份数据、同样的参数，两次运行结果不同，原因通常是random_state没设。GMM 的随机性来自两处：KMeans 初始化（KMeans 本身也有随机性）和n_init>1时的多次重启选择。设了random_state之后结果就固定了。

但设了种子之后还有一个更微妙的问题：标签编号不稳定。你把n_components从 3 改成 4，或者换个随机种子，簇的编号可能整体重排，原本的"簇 0"变成了"簇 2"。如果你的下游代码里有基于编号的硬编码逻辑，就会出大问题。

解决办法是建立一套稳定的簇标识规则，用可解释的特征去命名，而不是用编号。常用的规则有：按权重从大到小排序、按某个关键维度的均值排序、按簇心到原点的距离排序。我一般会在模型外面包一层，拟合完之后算一个从内部编号到业务名称的映射字典，所有下游逻辑都走这个字典。这样即使重新训练导致编号变化，业务语义也不会乱。

顺便说一个相关的坑：n_init大于 1 的时候，最终选中的是似然最高的那次，但不同次的运行可能选中不同的局部最优，编号自然也不同。如果你需要严格可比的结果（比如做 A/B 实验对比两个模型），除了固定random_state，还要固定n_init=1并用means_init手动传同一组初始中心。

6. 上线前的收尾工作：模型选择、持久化与算法分工

6.1 BIC/AIC 选 n_components 的实操细节

n_components怎么定？最常用的方法是扫一遍候选值，看 BIC 曲线。BIC 的定义里带了参数个数的惩罚项，所以它不会像似然那样单调变好，而是在某个点之后开始回升，那个拐点就是答案。

bics, aics = {}, {} for k in range(1, 9): m = GaussianMixture( n_components=k, covariance_type='full', n_init=5, random_state=42, reg_covar=1e-3 ).fit(X) bics[k] = m.bic(X) aics[k] = m.aic(X) best_k = min(bics, key=bics.get) print("best k by BIC:", best_k)

几个实操细节值得说。第一，扫描的时候n_init和random_state必须固定，否则 BIC 的差异里混入了初始化噪声，曲线会很难看。第二，BIC 倾向于选较小的模型，AIC 倾向选较大的，如果你的业务需要捕捉更细的结构，可以两个都看，BIC 给的 k 作为下界，AIC 给的作为上界。第三，bic和aic都是方法不是属性，必须传 X 进去，而且是负数还是正数别记混——越小越好。

一个经常被忽略的点：BIC 的绝对数值没有意义，只有相对比较有意义。你拿两个不同特征工程版本跑出来的 BIC 不能直接比，因为样本维度变了。要在同一套预处理下比较。

还有，如果 BIC 曲线出现了两个局部极小值，那说明数据里可能真的存在两种尺度的结构。这时候不要盲选，把两个 k 都训练出来，看它们在业务指标上的表现，往往小的那个对应粗粒度分群，大的那个对应细分，两个都能用，取哪个取决于你要回答什么问题。

6.2 模型持久化与推理阶段的对齐

GMM 模型的持久化用 joblib 就够，整个对象序列化下来一般几百 KB 到几 MB，取决于n_components和维度。真正要小心的不是序列化本身，而是训练和推理两端的预处理必须完全一致。

我见过的最常见事故是：训练时对特征做了标准化，用的是训练集拟合出来的StandardScaler，但线上推理时直接用原始数据喂给模型，或者用线上实时数据的均值方差重新标准化。结果就是score_samples的输出整体偏移，异常阈值完全失效。正确做法是把 scaler 和 GMM 打包成一个 pipeline 一起存，推理时整体加载。

import joblib from sklearn.pipeline import Pipeline from sklearn.preprocessing import StandardScaler pipe = Pipeline([ ("scaler", StandardScaler()), ("gmm", GaussianMixture(n_components=3, covariance_type="full", n_init=10, reg_covar=1e-3, random_state=42)), ]) pipe.fit(X) joblib.dump(pipe, "gmm_pipeline.joblib") # 推理 loaded = joblib.load("gmm_pipeline.joblib") prob = loaded.predict_proba(X_new) # 注意 Pipeline 的 predict_proba 会走两步

用 Pipeline 的一个小副作用是，访问内部属性要写成pipe.named_steps["gmm"].weights_，而不是pipe.weights_，这一点在写监控代码时容易漏。另外注意 Pipeline 只有最后一个 estimator 的predict和predict_proba方法会被暴露，score_samples在 Pipeline 上不可直接用，需要取出来单独调用。

特征维度变化也是个大坑。如果你的特征集在迭代中加了列、删了列或者改了顺序，旧模型直接报形状错误还好，最怕的是顺序变了但维度没变，模型照常跑，结果全错。我给的做法是在 pipeline 里加一层列顺序校验，用ColumnTransformer或者干脆在入口处做 assert，把列名列表写死。

6.3 和 KMeans、DBSCAN 的分工

GMM 不是万能的，知道它什么时候该让位，比知道它怎么用更重要。

KMeans 的优势是快、简单、结果可解释。数据量在千万级、簇形状接近球形、只需要一个粗粒度标签的时候，它依然是最好的选择。GMM 的训练成本明显更高，特别是covariance_type='full'加n_init>1的组合，十万样本、十维特征、五个分量，跑起来也要几秒到十几秒。

DBSCAN 解决的是另一类问题：任意形状、带噪声、不需要预先指定簇数。如果数据是环形、月牙形或者有大量孤立的噪声点，GMM 那套椭圆假设就不成立，怎么调都不对。这种情况下 DBSCAN 或者它的变体更合适。

我自己的选择顺序是这样：先画图看数据分布形状，如果是明显的椭圆团块结构，直接上 GMM；如果形状诡异但有密度分层，用 DBSCAN；如果只是要一个快的粗分，用 KMeans。还有一个组合用法我挺常用——先用 GMM 算出每个样本的对数密度，把密度极低的样本剔掉，再对剩余样本跑 KMeans。这样既避免了离群点污染质心，又保留了 KMeans 的速度。

最后一个实操心得是：GMM 的predict_proba输出在现代业务流程里其实比硬标签更有价值，因为它可以直接对接需要概率的决策逻辑（比如期望收益计算、风险敞口加权）。如果你的下游系统只接受标签，那用 KMeans 就行，硬套 GMM 反而增加了不必要的复杂度。我在实际项目里踩过这个坑——模型用 GMM 训得很好，结果下游只吃整数标签，软分配的信息全部浪费，最后又把责任度矩阵单独存了一份做加权统计，绕了一大圈。