Scikit-learn GaussianMixture 实战：3种协方差矩阵对比与聚类效果可视化

1. 高斯混合模型的核心概念

高斯混合模型（Gaussian Mixture Model, GMM）是一种基于概率的聚类算法，它假设数据是由多个高斯分布组合生成的。与K-means等硬聚类算法不同，GMM属于软聚类方法，能够给出每个样本属于各个簇的概率。

关键数学概念：

• 每个高斯分布由均值向量μ和协方差矩阵Σ定义
• 混合系数α表示各高斯分布的权重，满足∑α=1
• 概率密度函数：p(x) = ∑αᵢN(x|μᵢ,Σᵢ)

在实际应用中，GMM通过EM算法迭代优化这些参数：

1. E步：计算样本属于各分布的概率（责任值）
2. M步：根据责任值更新分布参数

2. 协方差矩阵类型解析

Scikit-learn的GaussianMixture提供了4种协方差矩阵类型，直接影响聚类形状：

from sklearn.mixture import GaussianMixture # 创建不同协方差类型的GMM模型 gmm_full = GaussianMixture(n_components=3, covariance_type='full') gmm_tied = GaussianMixture(n_components=3, covariance_type='tied') gmm_diag = GaussianMixture(n_components=3, covariance_type='diag')

3. 实战案例：鸢尾花数据集聚类

3.1 数据准备与可视化

import matplotlib.pyplot as plt from sklearn.datasets import load_iris import numpy as np # 加载数据 iris = load_iris() X = iris.data[:, :2] # 只使用前两个特征便于可视化 y = iris.target # 可视化原始数据 plt.scatter(X[:, 0], X[:, 1], c=y, s=40, cmap='viridis') plt.title("Iris True Labels") plt.xlabel("Sepal length") plt.ylabel("Sepal width")

3.2 不同协方差矩阵效果对比

我们创建可视化函数来比较三种主要协方差类型：

def plot_gmm(gmm, X, title): # 训练模型 gmm.fit(X) # 预测聚类结果 labels = gmm.predict(X) # 创建网格用于绘制决策边界 x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1 y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1 xx, yy = np.meshgrid(np.linspace(x_min, x_max, 100), np.linspace(y_min, y_max, 100)) Z = gmm.predict(np.c_[xx.ravel(), yy.ravel()]) Z = Z.reshape(xx.shape) # 绘制结果 plt.figure(figsize=(8, 6)) plt.contourf(xx, yy, Z, alpha=0.3) plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis') # 绘制椭圆表示协方差 for i in range(gmm.n_components): # 计算2D协方差椭圆 if gmm.covariance_type == 'full': cov = gmm.covariances_[i][:2, :2] elif gmm.covariance_type == 'tied': cov = gmm.covariances_[:2, :2] elif gmm.covariance_type == 'diag': cov = np.diag(gmm.covariances_[i][:2]) elif gmm.covariance_type == 'spherical': cov = np.eye(2) * gmm.covariances_[i] v, w = np.linalg.eigh(cov) angle = np.arctan2(w[0][1], w[0][0]) angle = 180 * angle / np.pi # 转换为度 v = 2. * np.sqrt(2.) * np.sqrt(v) ell = plt.matplotlib.patches.Ellipse(gmm.means_[i, :2], v[0], v[1], 180 + angle, color='black') ell.set_clip_box(plt.gca().bbox) ell.set_alpha(0.5) plt.gca().add_artist(ell) plt.title(title) plt.xlabel("Sepal length") plt.ylabel("Sepal width") # 比较三种主要类型 plot_gmm(gmm_full, X, "GMM with full covariance") plot_gmm(gmm_tied, X, "GMM with tied covariance") plot_gmm(gmm_diag, X, "GMM with diagonal covariance")

3.3 性能指标对比

我们使用轮廓系数和BIC（贝叶斯信息准则）评估模型：

from sklearn.metrics import silhouette_score models = { 'full': gmm_full, 'tied': gmm_tied, 'diag': gmm_diag, 'spherical': GaussianMixture(n_components=3, covariance_type='spherical') } results = [] for name, model in models.items(): model.fit(X) labels = model.predict(X) sil_score = silhouette_score(X, labels) bic = model.bic(X) results.append({ 'type': name, 'silhouette': sil_score, 'BIC': bic }) # 展示结果对比 import pandas as pd df_results = pd.DataFrame(results) print(df_results.sort_values('silhouette', ascending=False))

典型输出结果：

4. 高级应用技巧

4.1 确定最佳聚类数

使用BIC和轮廓系数确定最优组件数：

n_components = np.arange(1, 8) bics = [] silhouettes = [] for n in n_components: gmm = GaussianMixture(n_components=n, covariance_type='full') gmm.fit(X) bics.append(gmm.bic(X)) labels = gmm.predict(X) silhouettes.append(silhouette_score(X, labels)) # 绘制结果 fig, ax1 = plt.subplots() ax2 = ax1.twinx() ax1.plot(n_components, bics, 'b-', label='BIC') ax2.plot(n_components, silhouettes, 'r-', label='Silhouette') ax1.set_xlabel('Number of components') ax1.set_ylabel('BIC', color='b') ax2.set_ylabel('Silhouette', color='r') plt.title("Optimal Number of Components")

4.2 处理高维数据

当特征维度较高时，建议：

1. 先使用PCA降维
2. 选择'diag'或'spherical'减少参数数量
3. 使用BayesianGaussianMixture自动确定组件数

from sklearn.decomposition import PCA from sklearn.mixture import BayesianGaussianMixture # 高维数据处理示例 X_highdim = iris.data # 使用全部4个特征 # 方法1：PCA降维 pca = PCA(n_components=2) X_pca = pca.fit_transform(X_highdim) gmm_pca = GaussianMixture(n_components=3, covariance_type='full') gmm_pca.fit(X_pca) # 方法2：自动确定组件数 bgmm = BayesianGaussianMixture(n_components=10, covariance_type='diag') bgmm.fit(X_highdim) print(f"实际使用的组件数: {np.sum(bgmm.weights_ > 0.01)}")

5. 常见问题解决方案

问题1：模型收敛警告

解决方案：增加max_iter或tol参数，或尝试不同初始化

gmm = GaussianMixture(n_components=3, covariance_type='full', max_iter=500, tol=1e-4, init_params='kmeans')

问题2：奇异矩阵错误

解决方案：使用更简单的协方差类型，或添加正则化

gmm = GaussianMixture(n_components=3, covariance_type='diag', reg_covar=1e-6) # 添加小的正则化项

问题3：聚类结果不稳定

解决方案：固定random_state，或使用BayesianGaussianMixture

gmm = GaussianMixture(n_components=3, covariance_type='full', random_state=42)