Python案例如何用Scikit-learn做GAP统计量

wen python案例 2026-07-16 183

本文目录导读：

GAP统计量简介
完整实现代码
改进版本（使用并行计算）
使用注意事项

我将为您详细介绍如何使用Scikit-learn计算GAP统计量来选择最优聚类数K。

GAP统计量简介

GAP统计量通过比较原始数据的聚类紧密度与随机数据的期望紧密度来评估聚类效果,GAP值最大的K即为最优聚类数。

完整实现代码

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_blobs
def calculate_wcss(data, labels, centers):
    """
    计算Within-Cluster Sum of Squares (WCSS)
    """
    k = len(np.unique(labels))
    wcss = 0
    for i in range(k):
        # 获取属于该簇的数据点
        cluster_points = data[labels == i]
        if len(cluster_points) > 0:
            # 计算该簇中每个点到中心的距离平方和
            cluster_center = centers[i]
            distances = np.sum((cluster_points - cluster_center) ** 2, axis=1)
            wcss += np.sum(distances)
    return wcss
def generate_uniform_data(data, n_references=10):
    """
    在数据范围内生成均匀分布参考数据
    """
    n_samples, n_features = data.shape
    references = []
    for _ in range(n_references):
        # 在原始数据范围内生成随机数据
        ref_data = np.random.uniform(
            low=np.min(data, axis=0),
            high=np.max(data, axis=0),
            size=(n_samples, n_features)
        )
        references.append(ref_data)
    return np.array(references)
def compute_gap_statistic(data, k_max=10, n_references=10, random_state=42):
    """
    计算GAP统计量
    参数:
    data: 输入数据
    k_max: 最大聚类数
    n_references: 参考数据集数量
    random_state: 随机种子
    返回:
    gaps: GAP统计量值
    K: 最优聚类数
    """
    from sklearn.cluster import KMeans
    n_samples = data.shape[0]
    # 标准化数据
    scaler = StandardScaler()
    data_scaled = scaler.fit_transform(data)
    # 生成参考数据集
    reference_data = generate_uniform_data(data_scaled, n_references)
    gaps = []
    wcss_original = []
    wcss_references = []
    for k in range(1, k_max + 1):
        print(f"计算 K={k} 的GAP统计量...")
        # 1. 计算原始数据的WCSS
        kmeans_original = KMeans(n_clusters=k, random_state=random_state, n_init=10)
        labels_original = kmeans_original.fit_predict(data_scaled)
        centers_original = kmeans_original.cluster_centers_
        wcss_orig = calculate_wcss(data_scaled, labels_original, centers_original)
        wcss_original.append(wcss_orig)
        # 2. 计算参考数据的WCSS
        wcss_refs = []
        for ref_idx in range(n_references):
            ref_data = reference_data[ref_idx]
            kmeans_ref = KMeans(n_clusters=k, random_state=random_state + ref_idx, n_init=10)
            labels_ref = kmeans_ref.fit_predict(ref_data)
            centers_ref = kmeans_ref.cluster_centers_
            wcss_ref = calculate_wcss(ref_data, labels_ref, centers_ref)
            wcss_refs.append(wcss_ref)
        wcss_references.append(wcss_refs)
        # 3. 计算GAP统计量
        # gap(k) = E*[log(W_k)] - log(W_k)
        log_wcss_refs = np.log(wcss_refs)
        expected_log_wcss = np.mean(log_wcss_refs)
        gap = expected_log_wcss - np.log(wcss_orig)
        gaps.append(gap)
    # 选择最优K
    gaps = np.array(gaps)
    # 计算gap的标准差
    sd_k = []
    for k_idx in range(len(gaps)):
        if len(wcss_references[k_idx]) > 0:
            log_wcss_refs = np.log(wcss_references[k_idx])
            sd = np.std(log_wcss_refs)
            # 修正因子
            sd_k.append(sd * np.sqrt(1 + 1/n_references))
        else:
            sd_k.append(0)
    sd_k = np.array(sd_k)
    # 找到满足 gap(k) >= gap(k+1) - sd(k+1) 的最小k
    optimal_k = 1
    for k in range(1, len(gaps)):
        if gaps[k-1] >= gaps[k] - sd_k[k]:
            optimal_k = k
            break
    return gaps, sd_k, optimal_k
def plot_gap_statistic(gaps, sd_k, optimal_k):
    """
    可视化GAP统计量
    """
    plt.figure(figsize=(10, 6))
    k_values = range(1, len(gaps) + 1)
    # 绘制GAP值
    plt.plot(k_values, gaps, 'b-o', linewidth=2, markersize=8, label='GAP Statistic')
    # 绘制误差带
    plt.fill_between(k_values, gaps - sd_k, gaps + sd_k, 
                     alpha=0.2, color='blue', label='±1 SD')
    # 标记最优K
    plt.axvline(x=optimal_k, color='r', linestyle='--', 
                label=f'Optimal K = {optimal_k}', linewidth=2)
    plt.xlabel('Number of Clusters (K)', fontsize=12)
    plt.ylabel('GAP Statistic', fontsize=12)
    plt.title('GAP Statistic for Optimal K Selection', fontsize=14)
    plt.legend(fontsize=11)
    plt.grid(True, alpha=0.3)
    # 添加注释
    plt.text(optimal_k + 0.2, gaps[optimal_k-1], 
             f'Optimal K = {optimal_k}', 
             fontsize=12, color='red', fontweight='bold')
    plt.tight_layout()
    plt.show()
# 示例用法
if __name__ == "__main__":
    # 生成示例数据
    n_samples = 500
    # 生成3个聚类
    X, y_true = make_blobs(n_samples=n_samples, centers=3, 
                           cluster_std=0.8, random_state=42)
    print("数据形状:", X.shape)
    print("实际聚类数:", len(np.unique(y_true)))
    # 计算GAP统计量
    gaps, sd_k, optimal_k = compute_gap_statistic(
        X, 
        k_max=10, 
        n_references=10,
        random_state=42
    )
    print(f"\n最优聚类数 K = {optimal_k}")
    print(f"实际聚类数 = {len(np.unique(y_true))}")
    # 可视化结果
    plot_gap_statistic(gaps, sd_k, optimal_k)
    # 打印详细结果
    print("\n详细GAP统计量:")
    for k in range(1, len(gaps) + 1):
        print(f"K={k}: GAP={gaps[k-1]:.4f} ± {sd_k[k-1]:.4f}")

改进版本（使用并行计算）

import numpy as np
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from joblib import Parallel, delayed
import warnings
warnings.filterwarnings('ignore')
class GapStatistics:
    """
    GAP统计量计算类，支持并行计算
    """
    def __init__(self, n_references=10, random_state=42, n_jobs=-1):
        self.n_references = n_references
        self.random_state = random_state
        self.n_jobs = n_jobs
        self.gaps_ = None
        self.sd_k_ = None
        self.optimal_k_ = None
    def _compute_wcss_for_k(self, k, data_scaled, reference_data):
        """计算单个K值的WCSS"""
        n_references = len(reference_data)
        # 原始数据的WCSS
        kmeans = KMeans(n_clusters=k, random_state=self.random_state, n_init=10)
        labels = kmeans.fit_predict(data_scaled)
        centers = kmeans.cluster_centers_
        # 计算WCSS
        wcss_orig = 0
        for i in range(k):
            cluster_points = data_scaled[labels == i]
            if len(cluster_points) > 0:
                distances = np.sum((cluster_points - centers[i]) ** 2)
                wcss_orig += distances
        # 参考数据的WCSS
        wcss_refs = []
        for ref_idx in range(n_references):
            ref_data = reference_data[ref_idx]
            kmeans_ref = KMeans(n_clusters=k, 
                               random_state=self.random_state + ref_idx, 
                               n_init=10)
            labels_ref = kmeans_ref.fit_predict(ref_data)
            centers_ref = kmeans_ref.cluster_centers_
            wcss_ref = 0
            for i in range(k):
                cluster_points = ref_data[labels_ref == i]
                if len(cluster_points) > 0:
                    distances = np.sum((cluster_points - centers_ref[i]) ** 2)
                    wcss_ref += distances
            wcss_refs.append(wcss_ref)
        return wcss_orig, wcss_refs
    def fit(self, data, k_max=10):
        """
        计算GAP统计量
        """
        # 标准化数据
        scaler = StandardScaler()
        data_scaled = scaler.fit_transform(data)
        n_samples, n_features = data_scaled.shape
        # 生成参考数据集
        reference_data = []
        for _ in range(self.n_references):
            ref_data = np.random.uniform(
                low=np.min(data_scaled, axis=0),
                high=np.max(data_scaled, axis=0),
                size=(n_samples, n_features)
            )
            reference_data.append(ref_data)
        # 并行计算各K值的WCSS
        results = Parallel(n_jobs=self.n_jobs)(
            delayed(self._compute_wcss_for_k)(k, data_scaled, reference_data)
            for k in range(1, k_max + 1)
        )
        # 处理结果
        wcss_original = [r[0] for r in results]
        wcss_references = [r[1] for r in results]
        # 计算GAP统计量
        gaps = []
        sd_k = []
        for k_idx in range(k_max):
            log_wcss_refs = np.log(wcss_references[k_idx])
            expected_log_wcss = np.mean(log_wcss_refs)
            gap = expected_log_wcss - np.log(wcss_original[k_idx])
            gaps.append(gap)
            sd = np.std(log_wcss_refs)
            sd_k.append(sd * np.sqrt(1 + 1/self.n_references))
        self.gaps_ = np.array(gaps)
        self.sd_k_ = np.array(sd_k)
        # 选择最优K
        for k in range(1, len(gaps)):
            if gaps[k-1] >= gaps[k] - sd_k[k]:
                self.optimal_k_ = k
                break
        else:
            self.optimal_k_ = len(gaps)
        return self
    def plot(self, figsize=(10, 6)):
        """可视化GAP统计量"""
        import matplotlib.pyplot as plt
        plt.figure(figsize=figsize)
        k_values = range(1, len(self.gaps_) + 1)
        plt.plot(k_values, self.gaps_, 'b-o', linewidth=2, markersize=8, 
                label='GAP Statistic')
        plt.fill_between(k_values, 
                        self.gaps_ - self.sd_k_, 
                        self.gaps_ + self.sd_k_, 
                        alpha=0.2, color='blue', label='±1 SD')
        plt.axvline(x=self.optimal_k_, color='r', linestyle='--', 
                   label=f'Optimal K = {self.optimal_k_}', linewidth=2)
        plt.xlabel('Number of Clusters (K)', fontsize=12)
        plt.ylabel('GAP Statistic', fontsize=12)
        plt.title('GAP Statistic for Optimal K Selection', fontsize=14)
        plt.legend(fontsize=11)
        plt.grid(True, alpha=0.3)
        plt.tight_layout()
        plt.show()
# 使用示例
if __name__ == "__main__":
    from sklearn.datasets import make_blobs
    # 生成示例数据
    X, y = make_blobs(n_samples=500, centers=3, cluster_std=0.8, random_state=42)
    # 计算GAP统计量
    gap_stats = GapStatistics(n_references=10, random_state=42, n_jobs=-1)
    gap_stats.fit(X, k_max=10)
    print(f"最优聚类数: K = {gap_stats.optimal_k_}")
    gap_stats.plot()