ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

用BIC确定GMM聚类簇数:从原理到工程实践

用BIC确定GMM聚类簇数:从原理到工程实践 简介针对高斯混合模型GMM聚类时“到底该分几类”这一关键难题资源以贝叶斯信息准则BIC为核心提供一套自动确定最优聚类簇数的Python实现。使用时只需准备数据集脚本会自动遍历多个候选簇数逐一训练GMM并计算BIC最终返回使BIC取最小值的簇数——该值在拟合精度与模型复杂度之间取得平衡可有效替代人工目测或反复试参。适合正在学习聚类与模型选择的机器学习初学者也可供数据分析人员快速借用或嵌入自身流程。资源为一个轻量py脚本完整压缩包仅2KB大小代码基于sklearn的GaussianMixture类编写包含循环候选簇数、采用score_samples计算对数似然、代入BIC公式比较等核心步骤逻辑清晰、易读易改还可扩展AIC等其他准则进行交叉验证。脚本已吸引996人学习下载是理解BIC准则与GMM实战结合的良好样例。1. 为什么聚类数不能拍脑袋BIC 与 GMM 之间的双人舞给一堆没有标签的数据做高斯混合模型GMM聚类最难回答的问题往往不是“怎么分”而是“分成几簇”。K-Means 时代还能画个肘部图硬看换成 GMM簇数 K 直接决定模型复杂度——选小欠拟合选大过拟合定错一个后面所有下游分析都得重跑。BICBayesian Information Criterion贝叶斯信息准则就是回答“K 取多少”的量化工具它把模型拟合优度和参数惩罚放进同一个指标BIC 最小的 K 就是统计意义上的候选最优簇数。常用做法是在 K1 到 K15 之间逐个拟合 GMM画出 BIC 曲线找谷底。下面把原理、代码、参数设置和真实项目里容易踩的坑一并讲清适合已经跑通 GMM 基础代码、正卡在“K 值怎么定”上的从业者。2. GMM 如何计算 BIC从似然函数到信息准则的推导2.1 高斯混合模型参数与 EM 算法GMM 假设数据由 K 个高斯分布加权叠加产生每个分量有自己的均值 μ_k 和协方差 Σ_k另有权重 π_k所有权重之和为 1。聚类输出是软标签每个样本得到一组后验概率表示它属于每个簇的可能性。这种软分配特性让 GMM 在面对类别边界重叠的数据时比 K-Means 更从容——K-Means 只有一个距离归属GMM 天生能表达“这个样本七成属于 A 簇、三成属于 B 簇”。如果你在头歌这类机器学习实训平台上刷过 GMM 实验题大概率见过那种三团点云靠在一起的二维示意图那就是一个三高斯的混合分布。实际调用 sklearn.mixture.GaussianMixture 时只要传 n_components 设定簇数再 fit 就能完成训练。但 EM 算法的大致流程还是值得理解E 步根据当前参数计算每个样本的责任度即后验概率M 步用责任度加权更新均值、协方差和权重反复迭代直到对数似然不再明显提升。EM 收敛性是最早要面对的玄学点。当数据维度较高、分量重叠严重或初始化位置偏离时EM 容易落入局部最优同一份数据、同一个 K换一个 random_stateBIC 数值可能差出几百甚至更多。很多人忽略这一点把 sklearn 默认的 n_init1 直接当成最终答案这往往是后来选簇翻车的根源。固定 random_state、把 n_init 提高到 10 左右是应该在项目一开始就养成的习惯而不是等结果不对时才补救。GMM 的迭代停止由 tol 和 max_iter 控制。tol 默认 1e-3表示对数似然提升小于该阈值就停max_iter 默认 100数据量大或分量数多时可能不够。每次 fit 完建议检查 gmm.converged_ 和 gmm.n_iter_若 converged_ 中存在 False说明 EM 超过最大迭代仍未收敛这时的 BIC 只能当作参考不能作为决策依据。2.2 BIC 的数学形式与 GMM 的参数计数BIC 的标准形式是BIC -2 * ln(L_hat) k * ln(n)L_hat 是模型在最大似然估计下的似然值k 是自由参数个数n 是样本量。第一项衡量拟合优度似然越高 BIC 越低第二项对复杂度做惩罚模型参数越多扣分越重。ln(n) 随样本量增大而变大因此样本越多BIC 对复杂模型的惩罚越严厉。这一条在 GMM 场景里尤其重要因为高斯混合模型的参数数量可以因为 covariance_type 的选择而剧烈变化。GMM 的 k 由三部分组成每个分量的均值向量、协方差矩阵和混合权重。均值参数 K × D权重参数 K - 1协方差参数取决于 covariance_type。full 模式每个分量有 D(D1)/2 个独立参数diag 模式每个分量 D 个tied 模式全部分量共享一个完整协方差矩阵因此共 D(D1)/2 个spherical 模式每个分量只有 1 个方差标量。我用一张表来演示 D5、K4 时的参数计数差异covariance_type每个分量协方差参数协方差总参数均值参数权重参数总参数 kfullD(D1)/2 154×15 604×5 204-1 383diagD 54×5 2020343tied15共享1520338spherical14×1 420327同样是 1000 个样本full 与 spherical 的惩罚项差别大约是 ln(1000) × (83-27) ≈ 386。这个量级的差距足以吞掉似然项带来的差异所以把 covariance_type 当作和 K 一样的搜索维度不是可做可不做的优化而是必要的决策项。即便换成 D50、K10 的 full 模型协方差参数会高达 10×50×51/212750 个这种参数爆炸会让 BIC 惩罚项压过一切此时直接选 diag 或先降维才是合理路径。sklearn 的 GaussianMixture.bic(X) 已经内置了这套计算直接用就好。想核对内部数值时可以打印 gmm._n_parameters() 看自由参数总数再对照上表逻辑确认自己理解正确。手算值和接口返回值对不上时多半是某个 cov_type 的参数计数公式写错了。2.3 为什么选 BIC 而不是 AIC 或肘部法则AIC 和 BIC 看起来非常像只是把惩罚项从 k×ln(n) 换成常数 2k。样本量一大BIC 对复杂模型的惩罚明显更重因此更保守选出的 K 通常更小。聚类项目的下游往往要对每个簇做业务解释BIC 的保守通常不容易出问题AIC 在高样本量下容易把噪声也拆成独立簇解释成本没准会翻倍。这就像买设备AIC 偏向把每个小毛病都当作独立故障处理BIC 则倾向于用一套更简洁的方案解释现象。肘部法则在 K-Means 里还能用但 GMM 场景不适用。GMM 的 BIC 曲线很少出现类似 SSE 那种尖锐拐点更多是一条逐渐走低的平滑曲线人工找拐点全凭感觉。层次聚类也类似sklearn 里 AgglomerativeClustering 画完树状图后切口位置取决于观察者换个人看可能就换个结果。BIC 至少把这个环节自动化了不依赖人眼拿出来的数字也更容易在评审时说清楚。需要提醒的是BIC 不是“真实簇数探测器”。它只比较同一批候选模型的相对优劣如果数据生成过程根本不是高斯混合BIC 选出的数字也就没有物理意义。所以正式分析前跑一次 t-SNE 或 UMAP 做目检确认数据确实呈现近似高斯的团状结构再谈用 BIC 定簇数。我过手的数据集里标准化前后跑 GMMBIC 绝对值会有变化但最优 K 通常不变如果变化明显说明某个特征主导了似然计算先用 StandardScaler 处理再重新评估。3. 用 BIC 确定 GMM 聚类簇数完整可复现的 Python 流程3.1 搭建最小实验枚举 K 并绘制 BIC 曲线先给一个能在本地直接跑通的最小脚本。我习惯把选簇阶段和建模阶段分离选簇阶段只关心 BIC建模阶段再精调其他参数。下面代码生成三簇模拟数据对 K1 到 10 逐个拟合 GMM记录 BIC 并画曲线。import numpy as np import matplotlib.pyplot as plt from sklearn.mixture import GaussianMixture from sklearn.datasets import make_blobs # 生成三簇模拟数据样本量 1000二维簇间标准差 0.8 X, y_true make_blobs( n_samples1000, centers3, cluster_std0.8, random_state42 ) K_range range(1, 11) bic_scores [] for k in K_range: gmm GaussianMixture( n_componentsk, covariance_typefull, # 先用 full后面会对比其他类型 n_init5, # 多次初始化降低 EM 局部最优影响 random_state42 # 固定随机种子保证可复现 ) gmm.fit(X) bic_scores.append(gmm.bic(X)) # sklearn 内置 BIC 接口 best_k K_range[int(np.argmin(bic_scores))] print(fBIC 最优簇数: {best_k}) plt.plot(K_range, bic_scores, markero) plt.xlabel(Number of components (K)) plt.ylabel(BIC) plt.title(BIC vs K for GMM) plt.grid(True) plt.show()代码逻辑很简单先构造三簇高斯点云然后对每个 K 训练一个 full 协方差 GMM。gmm.bic(X) 返回该模型的 BIC 值数值越小代表“拟合与复杂度”综合表现越好argmin 直接给出统计最优 K。n_init5 表示 EM 从 5 个不同起点开始训练取对数似然最高的一次random_state42 固定随机种子确保运行结果可复现。这两个参数组合起来能明显减少 EM 局部最优带来的 BIC 抖动。真实项目里一般不会只打印一个 best_k 就结束。我会把 BIC 曲线保存成图片连同一个按 K 排列的数值表一起归档方便后续排障和汇报。如果 K 可能在 15 以上把 K_range 拉长到 20如果数据量很大n_init 从 5 降到 2 或 3曲线大体趋势仍然保留。样本量小时不要盲目拉高 K比如 n200 时 K10 的 GMM 每个分量的有效样本量已经太少参数估计方差会很大BIC 曲线也会出现很多伪谷底。3.2 把 covariance_type 纳入搜索避免单一参数结构误导只跑 full 一种协方差结构不够因为数据各维度的相关结构未必支持 full 的复杂假设。我通常把四种 covariance_type 全部跑一遍生成对比曲线再做横向比较。下面这段代码直接输出每种结构下的最优 K 和最小 BIC。cov_types [full, tied, diag, spherical] K_range np.arange(1, 11) result {} for cov_type in cov_types: bics [] for k in K_range: gmm GaussianMixture( n_componentsk, covariance_typecov_type, n_init5, random_state42 ) gmm.fit(X) bics.append(gmm.bic(X)) result[cov_type] bics best_k K_range[int(np.argmin(bics))] print(f{cov_type}: 最优 K{best_k}, 最小 BIC{np.min(bics):.2f})这段代码遍历四种协方差结构每个结构内部再遍历 K1 到 10。输出的四对结果可以让你直观看到“协方差类型一变最优 K 就跟着变”有多普遍。如果 full 和 diag 都指向同一个 K这个结果可信度较高如果四个结构指向四个不同的 K说明数据本身对“簇”的假设不敏感需要回到业务层面确认到底要分多细或者先做特征工程。这里还要提一个容易忽视的参数 reg_covar。GMM 拟合时协方差矩阵可能变得奇异sklearn 会在对角线加一个小值保证可逆默认是 1e-6。当你发现某些 K 下 BIC 突然异常或者拟合直接报错可以尝试把 reg_covar 调到 1e-5 或 1e-4但不要超过 1e-3否则协方差结构被严重扭曲BIC 就失去了比较意义。数据经过标准化或白化后reg_covar 触发概率会明显降低。3.3 记录收敛状态并核对参数数量别让未收敛模型混进来最后这一步容易被忽略却是选簇阶段最重要的“质量门禁”。GaussianMixture 训练完成后converged_ 记录每次初始化的收敛结果n_iter_ 记录实际迭代次数_n_parameters() 返回参数总数。在选簇循环里把这些信息打印出来能挡住大量无效结果。from sklearn.mixture import GaussianMixture gmm GaussianMixture( n_components5, covariance_typefull, n_init10, max_iter200, tol1e-3, random_state42 ) gmm.fit(X) print(收敛状态:, gmm.converged_) # 每次初始化的收敛布尔值 print(实际迭代次数:, gmm.n_iter_) # 最优模型的迭代次数 print(参数总个数:, gmm._n_parameters()) # 核对 BIC 里的 k这段代码把 n_init 提高到 10并打印收敛状态、迭代次数和参数总数。如果某个 K 的 converged_ 全是 False说明 EM 到 max_iter 上限都没收敛对应的 BIC 不能参与比较。n_iter_ 如果总是贴着 max_iter200 的上限说明模型难收敛优先检查数据是否做了标准化特征之间量纲差异过大是 EM 难以收敛的常见原因。_n_parameters() 是私有属性不建议在生产逻辑里依赖但可以作为排查工具当你怀疑某个 K 的 BIC 看起来离谱时先对一下这个数和手算的 k 是否一致。我一般会把 StandardScaler 和数据降维放进选簇流程的前置步骤。GMM 对特征尺度敏感一个量纲特别大的特征会主导似然计算BIC 曲线因此失真。标准化之后重新跑一遍最优 K 如果没变结论更稳如果变了说明原结果很大程度只是量纲效应的产物。还有一点血泪经验K 的上限不要无脑拉到 20K 超过某个阈值后每个簇的样本量迅速减少BIC 曲线会出现一堆伪谷底。先用业务常识设定一个最大簇数比如“最多分 8 类”再去枚举比纯数据驱动靠谱得多。4. BIC 选簇避坑指南5 个高发问题与排查方法4.1 坑一协方差类型一变最优 K 就跟着变现象同一份数据用 covariance_typefull 时 BIC 给出 K4换成 diag 后最优 K 变成 K6。原因不同协方差结构对应不同的模型复杂度惩罚项差异在 BIC 里占比很大数据维度越高full 与 diag 的参数数量差距越大两个模型下的 BIC 曲线形态可能完全不同。解决不要在选簇阶段只试一种协方差结构把四种 cov_type 全部跑一遍看它们指向的 K 是否一致。如果一致选那个公共 K如果不一致按 3.2 节的方式画出四条曲线综合判断。另一种思路是干脆把 cov_type 也作为搜索参数选 BIC 最低的那个组合。4.2 坑二EM 陷入局部最优BIC 每次运行结果不同现象固定 K、固定数据BIC 连续跑三次得到三个不同值BIC 曲线也完全变形。原因GMM 的似然面在高维空间有大量局部最优EM 对初始化很敏感random_state 一变就掉进不同的坑。解决把 n_init 从默认的 1 提到 10 或 20让 sklearn 从多个起始点出发再用一个固定的 random_state 保证实验可复现。如果这样还是不稳定考虑做标准化或 PCA 降维把问题约束在一个条件更好的空间里。调完参数之后回头检查 gmm.converged_确保模型真的收敛不要把未收敛模型的 BIC 拿来比较。4.3 坑三样本量太小BIC 惩罚项失灵现象数据只有一两百条BIC 曲线剧烈抖动无法找到清晰谷底甚至选出的 K 等于 1。原因BIC 的惩罚项是 k×ln(n)n 很小时惩罚力度弱而似然项本身受噪声影响大曲线不稳定。解决先看样本量是否支持你想要的聚类粒度。一个经验法则是每个簇至少要有几十个有效样本且样本量要远大于特征维度。如果样本量不足可以先用 PCA 降到低维再聚类或者改用 BayesianGaussianMixture它给分量权重加了狄利克雷过程先验会在拟合过程中自动把多余的簇“杀死”本质上达到了类似 BIC 的模型选择效果在小样本下更稳健。4.4 坑四BIC 曲线单调递减找不到极小值现象K 从 1 增大到 15BIC 一直下降没有谷底。原因数据本身不满足“几个高斯簇”的假设例如有流形结构、长尾分布或者特征空间大到让更多分量总能带来似然提升。解决别硬用 GMM。可以先跑一遍 t-SNE 或 UMAP 看数据低维投影判断数据形态。若是流形结构用谱聚类或 DBSCAN若是长尾考虑对特征做变换。如果业务上仍然要 GMM可以给 K 加业务上界比如“最多分 6 类”然后只看 K1 到 6 范围内的相对最优不和全局最优死磕。4.5 坑五高维数据下 full 协方差参数爆炸现象特征维度几十上百时full 协方差的 BIC 在 K2 就到底选不出有意义的簇。原因full 协方差的参数个数是 O(D²)维度高时每个分量的参数数量爆炸BIC 惩罚项迅速压过似然项导致模型被推向极小的 K。解决维度低于 10 且各维度相关性明显时才考虑 full维度较高时优先用 diag假设各特征独立或 tied所有分量共享一个协方差矩阵。更稳妥的做法是先 PCA 降到 5 到 10 维再跑 GMM既保留主要结构又控制 BIC 里 k 的数量让曲线重新出现可判断的谷底。5. 验证 BIC 选出的簇数是否可信三种佐证方法5.1 用轮廓系数对照 BIC 决策BIC 是从统计模型角度做选择轮廓系数是从“簇内紧密、簇间分离”的几何角度做评价。虽然 GMM 是软聚类但可以用最大后验概率把标签硬化再计算轮廓系数。如果 BIC 选出的 K 在轮廓系数下也明显优于相邻 K这个结果的把握就大了很多。from sklearn.metrics import silhouette_score def best_k_by_bic(X, K_range, cov_typefull, seed42): bics [] for k in K_range: gmm GaussianMixture(n_componentsk, covariance_typecov_type, n_init10, random_stateseed) gmm.fit(X) bics.append(gmm.bic(X)) return K_range[int(np.argmin(bics))], bics k_bic, _ best_k_by_bic(X, range(2, 11)) sil_scores [] for k in range(2, 11): gmm GaussianMixture(n_componentsk, covariance_typefull, n_init10, random_state42) labels gmm.fit_predict(X) sil_scores.append(silhouette_score(X, labels)) best_k_sil range(2, 11)[int(np.argmax(sil_scores))] print(fBIC 最优 K: {k_bic}, 轮廓系数最优 K: {best_k_sil})说明轮廓系数取值范围 -1 到 1越高说明簇间分离越明显。它天然偏向紧凑的凸形簇对狭长簇不友好所以只能做佐证不能单独替代 BIC。当两个指标指向同一个 K 时这个 K 通常会在大样本项目中表现稳定当两者不一致时优先看数据形态——如果数据不是标准团状轮廓系数的参考价值要打折。5.2 用交叉验证的对数似然检验泛化能力BIC 对全量数据做一次拟合后给出评分但一个信息准则毕竟不是“未来数据”的直接指标。K 折交叉验证的做法是把样本拆成训练和验证两部分在训练集上拟合 GMM在验证集上计算对数似然最后取平均。验证集对数似然越高代表模型泛化能力越好。如果交叉验证的最优 K 与 BIC 的最优 K 基本一致说明这个簇数不是过拟合出来的。from sklearn.model_selection import KFold kf KFold(n_splits5, shuffleTrue, random_state42) cv_scores [] for k in range(1, 11): val_lls [] for train_idx, val_idx in kf.split(X): gmm GaussianMixture(n_componentsk, covariance_typefull, n_init5, random_state42) gmm.fit(X[train_idx]) val_lls.append(gmm.score(X[val_idx])) # score 返回平均对数似然 cv_scores.append(np.mean(val_lls)) best_k_cv range(1, 11)[int(np.argmax(cv_scores))] print(f交叉验证最优 K: {best_k_cv}, 平均对数似然: {max(cv_scores):.2f})说明gmm.score(X) 返回的是平均对数似然越大越好。这里每折都重新拟合一次 GMM耗时约为单次拟合的 K 倍但换来一个更可信的选簇依据。注意验证集要独立于训练集绝不能用全量数据的 BIC 和交叉验证结果做同一批样本的比较否则逻辑上是循环论证。交叉验证的结果也能暴露出一个隐含问题如果每折的最优 K 互相打架说明数据本身不稳定这时候更要依赖业务判断。5.3 用 Bootstrap 冲击检验看 BIC 的稳定性BIC 是一个点估计它没有自带“置信区间”。Bootstrap 的思想是对原始数据做有放回的重采样在每次重采样上重新执行 BIC 选簇统计最优 K 的分布。如果最优 K 每次都稳定出现说明选择是可信的如果分布很散说明数据本身没有提供足够的簇结构信号BIC 的结果更多是噪声。rng np.random.default_rng(42) K_range np.arange(1, 11) best_ks [] for i in range(200): idx rng.integers(0, X.shape[0], X.shape[0]) X_boot X[idx] bics [] for k in K_range: gmm GaussianMixture(n_componentsk, covariance_typefull, n_init3, random_state42) gmm.fit(X_boot) bics.append(gmm.bic(X_boot)) best_ks.append(int(K_range[int(np.argmin(bics))])) counts np.bincount(best_ks, minlengthK_range[-1] 1) for k in K_range: print(fK{k}: {counts[k]} 次 ({counts[k] / len(best_ks) * 100:.1f}%))说明这段代码跑 200 次 bootstrap每次从原始样本中重采样同等数量的样本重新选 K。当数据包含明显的三个簇时K3 的出现比例会压倒性领先如果结果在两个 K 上各占一半说明 BIC 的“最优”很不稳定你需要认真对待业务侧的解释。Bootstrap 的代价是计算量成倍上涨实践中可以先跑 50 次预检曲线稳定了再放大到 200 次。6. 两阶段决策法把 BIC 选簇做成可落地的工程流程在真实项目里我很少直接拿 BIC 最小值当最终答案。信息准则回答的是“统计上最合理”不是“业务上最合理”。我的习惯是走一个两阶段流程。第一阶段把 K1 到 15、四种 covariance_type 的网格全部跑一遍画出全部 BIC 曲线圈出每个协方差结构下的局部极小点得到候选 K 集合。第二阶段用轮廓系数和交叉验证对数似然做两轮筛选把候选 K 收窄到 1 到 2 个。如果两个候选在业务上都说得通选较少的那个——更简单的模型更容易维护对噪声的容忍度也更高。几个工程细节值得记下来。数据量超过十万条时EM 迭代明显变慢可以把 n_init 降到 3 或 5BIC 曲线的主要趋势依然存在维度超过 20 时默认只跑 diag 和 tiedfull 的参数数量太夸张。每次实验固定 random_state把 BIC、AIC、轮廓系数和交叉验证得分写成一个字典并保存为 CSV后续汇报时直接把图和表丢出来别人需要复现时也用同一套脚本跑。我早期做聚类分析吃过一次亏当时 BIC 明确指向 K5我直接给业务方交付了五类分群结果业务方回来说“第五类是我们统计口径下的重复样本”。后来跑 Bootstrap 重采样检查发现最优 K 在 4 和 5 之间几乎各占一半说明数据根本没有足够信号把第五簇稳定分离。从那以后我给自己设了一条规矩BIC 只是提案者不是裁决者。任何信息准则给出的簇数都要经过至少一个独立验证和一轮业务解释才能进入生产流程。最后分享一个实用技巧把第 5 章的三种验证方法封装成一个 auto_select_k 函数对每个新数据集输出一张多指标对照表列出每个 K 下的 BIC、轮廓系数、交叉验证对数似然。这张表既能帮你积累不同数据集的选簇经验也能在评审的时候用一张表说清楚“为什么是这个 K”。希望这套流程能帮你在自己的聚类项目里少走几个弯路。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进