Clustering Analysis

聚类分析

概述

聚类会将数据划分为由相似观测值组成的组，无需预先定义标签，从而帮助发现数据中的自然模式和结构。

适用场景

• 根据购买行为或人口统计特征细分客户
• 在不了解类别的情况下发现数据中的自然分组
• 识别市场细分，以便开展定向营销活动
• 将大型数据集整理为有意义的类别，供后续分析
• 在基因表达数据或医学影像中发现模式
• 根据相似度对文档、产品或用户分组，用于推荐系统

聚类算法

• K-Means：将数据划分为 k 个簇
• Hierarchical：通过树状图展示嵌套簇
• DBSCAN：基于密度的任意形状聚类
• Gaussian Mixture：概率聚类
• Agglomerative：自底向上的层次聚类方法

关键概念

• Cluster Validation：评估聚类质量的指标
• Optimal Clusters：确定最佳 k 的方法
• Inertia：簇内平方和
• Silhouette Score：衡量簇间分离程度
• Dendrogram：层次聚类的可视化

使用 Python 实现

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans, DBSCAN, AgglomerativeClustering
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import (
    silhouette_score, silhouette_samples, davies_bouldin_score,
    calinski_harabasz_score
)
from scipy.cluster.hierarchy import dendrogram, linkage
import seaborn as sns

# Generate sample data
np.random.seed(42)
n_samples = 300
centers = [[0, 0], [5, 5], [-3, 4]]
X = np.vstack([
    np.random.randn(100, 2) + centers[0],
    np.random.randn(100, 2) + centers[1],
    np.random.randn(100, 2) + centers[2],
])

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

# K-Means with Elbow method
inertias = []
silhouette_scores = []
k_range = range(2, 11)

for k in k_range:
    kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
    kmeans.fit(X_scaled)
    inertias.append(kmeans.inertia_)
    silhouette_scores.append(silhouette_score(X_scaled, kmeans.labels_))

fig, axes = plt.subplots(1, 2, figsize=(14, 4))

axes[0].plot(k_range, inertias, 'bo-')
axes[0].set_xlabel('Number of Clusters (k)')
axes[0].set_ylabel('Inertia')
axes[0].set_title('Elbow Method')
axes[0].grid(True, alpha=0.3)

axes[1].plot(k_range, silhouette_scores, 'go-')
axes[1].set_xlabel('Number of Clusters (k)')
axes[1].set_ylabel('Silhouette Score')
axes[1].set_title('Silhouette Analysis')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Optimal k = 3
optimal_k = 3
kmeans = KMeans(n_clusters=optimal_k, random_state=42, n_init=10)
kmeans_labels = kmeans.fit_predict(X_scaled)

# K-Means visualization
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# K-Means clusters
axes[0].scatter(X[:, 0], X[:, 1], c=kmeans_labels, cmap='viridis', alpha=0.6)
axes[0].scatter(
    kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1],
    c='red', marker='X', s=200, edgecolors='black', linewidths=2
)
axes[0].set_title(f'K-Means (k={optimal_k})')
axes[0].set_xlabel('Feature 1')
axes[0].set_ylabel('Feature 2')

# Silhouette plot
ax = axes[1]
y_lower = 10
silhouette_vals = silhouette_samples(X_scaled, kmeans_labels)

for i in range(optimal_k):
    cluster_silhouette_vals = silhouette_vals[kmeans_labels == i]
    cluster_silhouette_vals.sort()

    size_cluster_i = cluster_silhouette_vals.shape[0]
    y_upper = y_lower + size_cluster_i

    ax.fill_betweenx(np.arange(y_lower, y_upper),
                      0, cluster_silhouette_vals,
                      alpha=0.7, label=f'Cluster {i}')
    y_lower = y_upper + 10

ax.axvline(x=silhouette_score(X_scaled, kmeans_labels), color="red", linestyle="--")
ax.set_xlabel('Silhouette Coefficient')
ax.set_ylabel('Cluster Label')
ax.set_title('Silhouette Plot')

# Hierarchical clustering
linkage_matrix = linkage(X_scaled, method='ward')
dendrogram(linkage_matrix, ax=axes[2], truncate_mode='lastp', p=10)
axes[2].set_title('Dendrogram (Ward)')
axes[2].set_xlabel('Sample Index')

plt.tight_layout()
plt.show()

# Hierarchical clustering
hierarchical = AgglomerativeClustering(n_clusters=optimal_k, linkage='ward')
hier_labels = hierarchical.fit_predict(X_scaled)

# DBSCAN clustering
dbscan = DBSCAN(eps=0.4, min_samples=5)
dbscan_labels = dbscan.fit_predict(X_scaled)
n_clusters_dbscan = len(set(dbscan_labels)) - (1 if -1 in dbscan_labels else 0)
n_noise = list(dbscan_labels).count(-1)

# Gaussian Mixture Model
gmm = GaussianMixture(n_components=optimal_k, random_state=42)
gmm_labels = gmm.fit_predict(X_scaled)
gmm_proba = gmm.predict_proba(X_scaled)

# Clustering algorithm comparison
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

algorithms = [
    (kmeans_labels, 'K-Means'),
    (hier_labels, 'Hierarchical'),
    (dbscan_labels, 'DBSCAN'),
    (gmm_labels, 'Gaussian Mixture'),
]

for idx, (labels, title) in enumerate(algorithms):
    ax = axes[idx // 2, idx % 2]

    # Skip noise points for DBSCAN
    mask = labels != -1
    scatter = ax.scatter(
        X[mask, 0], X[mask, 1], c=labels[mask], cmap='viridis', alpha=0.6
    )

    if title == 'DBSCAN' and n_noise > 0:
        noise_mask = labels == -1
        ax.scatter(X[noise_mask, 0], X[noise_mask, 1], c='red', marker='x', s=100, label='Noise')
        ax.legend()

    ax.set_title(f'{title} (n_clusters={len(set(labels[mask]))})')
    ax.set_xlabel('Feature 1')
    ax.set_ylabel('Feature 2')

plt.tight_layout()
plt.show()

# Cluster validation metrics
validation_metrics = {
    'Algorithm': ['K-Means', 'Hierarchical', 'DBSCAN', 'GMM'],
    'Silhouette Score': [
        silhouette_score(X_scaled, kmeans_labels),
        silhouette_score(X_scaled, hier_labels),
        silhouette_score(X_scaled[dbscan_labels != -1], dbscan_labels[dbscan_labels != -1]) if n_noise < len(X_scaled) else np.nan,
        silhouette_score(X_scaled, gmm_labels),
    ],
    'Davies-Bouldin Index': [
        davies_bouldin_score(X_scaled, kmeans_labels),
        davies_bouldin_score(X_scaled, hier_labels),
        davies_bouldin_score(X_scaled[dbscan_labels != -1], dbscan_labels[dbscan_labels != -1]) if n_noise < len(X_scaled) else np.nan,
        davies_bouldin_score(X_scaled, gmm_labels),
    ],
    'Calinski-Harabasz Index': [
        calinski_harabasz_score(X_scaled, kmeans_labels),
        calinski_harabasz_score(X_scaled, hier_labels),
        calinski_harabasz_score(X_scaled[dbscan_labels != -1], dbscan_labels[dbscan_labels != -1]) if n_noise < len(X_scaled) else np.nan,
        calinski_harabasz_score(X_scaled, gmm_labels),
    ],
}

metrics_df = pd.DataFrame(validation_metrics)
print("Clustering Validation Metrics:")
print(metrics_df)

# Cluster size analysis
sizes_df = pd.DataFrame({
    'K-Means': pd.Series(kmeans_labels).value_counts().sort_index(),
    'Hierarchical': pd.Series(hier_labels).value_counts().sort_index(),
    'GMM': pd.Series(gmm_labels).value_counts().sort_index(),
})

print("\nCluster Sizes:")
print(sizes_df)

# Membership probability (GMM)
fig, ax = plt.subplots(figsize=(10, 6))
membership = gmm_proba.max(axis=1)
scatter = ax.scatter(X[:, 0], X[:, 1], c=membership, cmap='RdYlGn', alpha=0.6, s=50)
ax.set_title('Cluster Membership Confidence (GMM)')
ax.set_xlabel('Feature 1')
ax.set_ylabel('Feature 2')
plt.colorbar(scatter, ax=ax, label='Membership Probability')
plt.show()

# Cluster characteristics
kmeans_centers_original = scaler.inverse_transform(kmeans.cluster_centers_)
cluster_df = pd.DataFrame(X, columns=['Feature 1', 'Feature 2'])
cluster_df['Cluster'] = kmeans_labels

for cluster_id in range(optimal_k):
    cluster_data = cluster_df[cluster_df['Cluster'] == cluster_id]
    print(f"\nCluster {cluster_id} Characteristics:")
    print(cluster_data[['Feature 1', 'Feature 2']].describe())

聚类质量指标

• Silhouette Score：范围为 -1 到 1（越高越好）
• Davies-Bouldin Index：越低越好
• Calinski-Harabasz Index：越高越好
• Inertia：越低越好（仅适用于 KMeans）

算法选择

• K-Means：速度快，适用于球形簇，需要指定 k
• Hierarchical：生成树状图，易于解释
• DBSCAN：支持任意形状，可处理噪声
• GMM：基于概率，可进行软分配

交付成果

• 最优簇数量分析
• 聚类可视化
• 验证指标对比
• 簇特征摘要
• 轮廓图
• 层次聚类树状图
• 成员分配结果