高斯混合模型实战：从鸢尾花数据到K值选择
1. 高斯混合模型的核心思想
高斯混合模型（Gaussian Mixture Model, GMM）是一种基于概率分布的聚类方法，它假设数据是由多个高斯分布组合生成的。与K-means等硬聚类方法不同，GMM属于软聚类，能够给出样本属于各个簇的概率。
为什么选择GMM而非其他聚类方法？
- 灵活性：通过调整协方差矩阵，GMM可以识别不同形状（球形、椭圆形）的簇
- 概率输出：不仅给出聚类结果，还提供样本属于各簇的概率
- 密度估计：能够估计整个数据空间的概率密度分布
在Scikit-learn 1.4中，GMM的实现进行了多项优化：
- 内存效率提升
- 计算速度加快
- 新增了收敛检测机制
2. 环境准备与数据加载
首先确保安装了正确版本的Scikit-learn：
pip install scikit-learn==1.4.0加载必要的库和鸢尾花数据集：
import numpy as np import matplotlib.pyplot as plt from sklearn import datasets from sklearn.mixture import GaussianMixture from sklearn.metrics import silhouette_score import pandas as pd # 加载鸢尾花数据集 iris = datasets.load_iris() X = iris.data # 使用所有四个特征 y = iris.target feature_names = iris.feature_names提示：虽然鸢尾花数据集有真实标签，但在聚类任务中我们通常不使用这些标签，仅用于最终评估聚类效果。
3. 模型训练与可视化
3.1 基础模型训练
我们先训练一个简单的GMM模型，假设已知类别数为3：
# 初始化GMM模型 gmm = GaussianMixture(n_components=3, random_state=42) # 训练模型 gmm.fit(X) # 获取聚类结果 labels = gmm.predict(X) probs = gmm.predict_proba(X) # 获取属于各簇的概率3.2 二维可视化
为了直观展示聚类效果，我们选择两个特征进行可视化：
# 选择两个特征进行可视化（花瓣长度和宽度） X_2d = X[:, [2, 3]] # 重新训练模型 gmm_2d = GaussianMixture(n_components=3, random_state=42) gmm_2d.fit(X_2d) labels_2d = gmm_2d.predict(X_2d) # 可视化 plt.figure(figsize=(10, 6)) plt.scatter(X_2d[:, 0], X_2d[:, 1], c=labels_2d, s=50, cmap='viridis') # 绘制高斯分布的均值和椭圆 for i in range(3): # 提取均值和协方差 mean = gmm_2d.means_[i] cov = gmm_2d.covariances_[i] # 绘制椭圆 v, w = np.linalg.eigh(cov) angle = np.degrees(np.arctan2(w[0][1], w[0][0])) v = 2. * np.sqrt(2.) * np.sqrt(v) ell = plt.matplotlib.patches.Ellipse(mean, v[0], v[1], angle, color='red', alpha=0.3) plt.gca().add_patch(ell) # 标记中心点 plt.scatter(mean[0], mean[1], marker='x', s=100, color='red') plt.xlabel(feature_names[2]) plt.ylabel(feature_names[3]) plt.title('GMM聚类结果（二维投影）') plt.colorbar(label='Cluster') plt.show()4. 选择最优K值：AIC与BIC
在实际应用中，我们往往不知道真实的簇数量。GMM提供了两种信息准则帮助选择K值：
- AIC (Akaike Information Criterion)：倾向于选择更复杂的模型
- BIC (Bayesian Information Criterion)：对模型复杂度惩罚更大
# 测试不同K值的AIC和BIC n_components_range = range(1, 8) aics = [] bics = [] for n_components in n_components_range: gmm = GaussianMixture(n_components=n_components, random_state=42) gmm.fit(X) aics.append(gmm.aic(X)) bics.append(gmm.bic(X)) # 可视化结果 plt.figure(figsize=(10, 6)) plt.plot(n_components_range, aics, label='AIC') plt.plot(n_components_range, bics, label='BIC') plt.xlabel('Number of Components (K)') plt.ylabel('Information Criterion') plt.legend() plt.title('AIC和BIC随K值变化曲线') plt.xticks(n_components_range) plt.show()如何解读AIC/BIC曲线？
| K值选择策略 | 适用场景 | 特点 |
|---|---|---|
| AIC最小值 | 数据质量高，噪声少 | 倾向于选择更复杂的模型 |
| BIC最小值 | 样本量较大 | 对模型复杂度惩罚更强 |
| 肘部法则 | 直观判断 | 选择曲线拐点处的K值 |
5. 模型评估与结果分析
5.1 聚类效果评估
虽然聚类是无监督学习，但我们仍有一些评估方法：
# 轮廓系数（越高越好） silhouette_avg = silhouette_score(X, labels) print(f"轮廓系数: {silhouette_avg:.3f}") # 与真实标签比较（仅在有标签时使用） from sklearn.metrics import adjusted_rand_score ari = adjusted_rand_score(y, labels) print(f"调整兰德指数: {ari:.3f}")5.2 参数分析
训练好的GMM模型包含以下关键参数：
# 输出模型参数 print("各簇权重:", gmm.weights_) print("\n均值向量:") for i, mean in enumerate(gmm.means_): print(f"簇{i}: {np.round(mean, 2)}") print("\n协方差矩阵:") for i, cov in enumerate(gmm.covariances_): print(f"簇{i}的协方差矩阵:") print(np.round(cov, 2))5.3 概率输出解析
GMM的一个独特优势是提供样本属于各簇的概率：
# 获取概率最高的前5个样本 prob_df = pd.DataFrame(probs, columns=[f'Cluster_{i}' for i in range(3)]) prob_df['Max_Cluster'] = prob_df.idxmax(axis=1) print(prob_df.head())6. 高级应用与技巧
6.1 协方差矩阵类型选择
Scikit-learn提供了四种协方差矩阵类型：
| 类型 | 参数 | 适用场景 |
|---|---|---|
| 全协方差 | 'full' | 各簇有任意方向的椭圆形状 |
| 对角协方差 | 'diag' | 各特征独立，椭圆轴平行坐标轴 |
| 球形协方差 | 'spherical' | 所有簇是圆形，大小相同 |
| 绑定协方差 | 'tied' | 所有簇共享同一协方差矩阵 |
# 比较不同协方差类型的效果 cov_types = ['full', 'diag', 'spherical', 'tied'] for cov_type in cov_types: gmm = GaussianMixture(n_components=3, covariance_type=cov_type, random_state=42) gmm.fit(X) labels = gmm.predict(X) score = silhouette_score(X, labels) print(f"{cov_type}协方差矩阵的轮廓系数: {score:.3f}")6.2 处理高维数据
当特征维度较高时，可以考虑：
- 特征选择：选择与聚类目标相关的特征
- 降维：使用PCA等降维方法
- 正则化：在协方差矩阵中添加小的常数防止奇异
# PCA降维后聚类 from sklearn.decomposition import PCA pca = PCA(n_components=2) X_pca = pca.fit_transform(X) gmm_pca = GaussianMixture(n_components=3, random_state=42) gmm_pca.fit(X_pca) labels_pca = gmm_pca.predict(X_pca) # 可视化 plt.figure(figsize=(10, 6)) plt.scatter(X_pca[:, 0], X_pca[:, 1], c=labels_pca, cmap='viridis') plt.xlabel('Principal Component 1') plt.ylabel('Principal Component 2') plt.title('PCA降维后的GMM聚类') plt.show()7. 实际应用中的注意事项
数据预处理：
- 标准化：GMM对特征的尺度敏感
- 异常值处理：高斯分布对异常值敏感
初始化策略：
- K-means初始化（默认）
- 随机初始化（多次运行取最优）
收敛判断：
- 设置合理的
max_iter和tol参数 - 监控对数似然的变化
- 设置合理的
# 带预处理和多次初始化的完整流程 from sklearn.preprocessing import StandardScaler from sklearn.pipeline import Pipeline # 创建管道 pipeline = Pipeline([ ('scaler', StandardScaler()), ('gmm', GaussianMixture(n_components=3, n_init=5, # 多次初始化 max_iter=500, tol=1e-4, random_state=42)) ]) pipeline.fit(X) final_labels = pipeline.named_steps['gmm'].predict(X)在实际项目中，我发现GMM对初始化比较敏感，特别是在高维空间中。一个实用的技巧是先用K-means进行粗聚类，然后用其结果初始化GMM：
from sklearn.cluster import KMeans # 使用K-means初始化 kmeans = KMeans(n_clusters=3, random_state=42) initial_means = kmeans.fit(X).cluster_centers_ gmm_smart = GaussianMixture(n_components=3, means_init=initial_means, random_state=42) gmm_smart.fit(X)
)
