本文目录导读：
我将为您详细介绍如何使用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()
使用注意事项
- 数据预处理：建议先标准化数据
- 计算复杂度：O(K N d)，K为最大聚类数，N为样本数，d为特征数
- 参考数据集：通常10-20个参考数据集即可
- K值范围：一般设置1到sqrt(N)之间
这个实现提供了完整的GAP统计量计算功能,包括并行计算优化和可视化支持。
