聚类算法评估指标之SC、CH轮廓系数介绍

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.metrics import silhouette_score, calinski_harabasz_score

# 1. 生成数据
print("=== 1. 生成数据 ===")
X, y_true = make_blobs(
    n_samples=500,
    centers=4,
    cluster_std=0.5,
    random_state=42
)
print(f"数据形状: {X.shape}")

# 2. SC轮廓系数
print("\n=== 2. SC轮廓系数 ===")

def demo_sc(X):
    """SC轮廓系数评估"""
    print("\n--- SC轮廓系数 ---")
    
    # K值范围（SC至少从2开始）
    k_range = range(2, 11)
    sc_list = []
    
    for k in k_range:
        kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
        labels = kmeans.fit_predict(X)
        sc = silhouette_score(X, labels)
        sc_list.append(sc)
        print(f"K={k}: SC={sc:.4f}")
    
    # 绘制折线图
    plt.figure(figsize=(10, 6))
    plt.plot(k_range, sc_list, 'ro-', linewidth=2, markersize=8)
    plt.xlabel('K值')
    plt.ylabel('SC轮廓系数')
    plt.title('SC轮廓系数选择最佳K值')
    plt.grid(True)
    plt.xticks(k_range)
    plt.savefig('sc_evaluation.png', dpi=100)
    print("SC评估图已保存")
    
    # 找最佳K值
    best_k = k_range[np.argmax(sc_list)]
    print(f"\n最佳K值: {best_k} (SC={max(sc_list):.4f})")
    
    return sc_list

sc_list = demo_sc(X)

# 3. CH轮廓系数
print("\n=== 3. CH轮廓系数 ===")

def demo_ch(X):
    """CH轮廓系数评估"""
    print("\n--- CH轮廓系数 ---")
    
    # K值范围（CH至少从2开始）
    k_range = range(2, 11)
    ch_list = []
    
    for k in k_range:
        kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
        labels = kmeans.fit_predict(X)
        ch = calinski_harabasz_score(X, labels)
        ch_list.append(ch)
        print(f"K={k}: CH={ch:.4f}")
    
    # 绘制折线图
    plt.figure(figsize=(10, 6))
    plt.plot(k_range, ch_list, 'go-', linewidth=2, markersize=8)
    plt.xlabel('K值')
    plt.ylabel('CH轮廓系数')
    plt.title('CH轮廓系数选择最佳K值')
    plt.grid(True)
    plt.xticks(k_range)
    plt.savefig('ch_evaluation.png', dpi=100)
    print("CH评估图已保存")
    
    # 找最佳K值
    best_k = k_range[np.argmax(ch_list)]
    print(f"\n最佳K值: {best_k} (CH={max(ch_list):.4f})")
    
    return ch_list

ch_list = demo_ch(X)

# 4. 手动计算SC
print("\n=== 4. 手动计算SC ===")

def calculate_sc_manual(X, labels):
    """手动计算SC轮廓系数"""
    n_samples = X.shape[0]
    unique_labels = np.unique(labels)
    sc_values = []
    
    for i in range(n_samples):
        # 当前样本的簇
        current_cluster = labels[i]
        
        # a值：到簇内其他样本的平均距离
        same_cluster_mask = (labels == current_cluster)
        same_cluster_points = X[same_cluster_mask]
        if len(same_cluster_points) > 1:
            a = np.mean(np.sqrt(np.sum((X[i] - same_cluster_points) ** 2, axis=1)))
        else:
            a = 0
        
        # b值：到其他簇的最小平均距离
        b_values = []
        for label in unique_labels:
            if label != current_cluster:
                other_cluster_points = X[labels == label]
                b_val = np.mean(np.sqrt(np.sum((X[i] - other_cluster_points) ** 2, axis=1)))
                b_values.append(b_val)
        
        if b_values:
            b = min(b_values)
        else:
            b = 0
        
        # SC值
        if max(a, b) > 0:
            sc = (b - a) / max(a, b)
        else:
            sc = 0
        sc_values.append(sc)
    
    return np.mean(sc_values)

# 使用K=4
kmeans = KMeans(n_clusters=4, random_state=42, n_init=10)
labels = kmeans.fit_predict(X)

sc_manual = calculate_sc_manual(X, labels)
sc_sklearn = silhouette_score(X, labels)
print(f"手动计算SC: {sc_manual:.4f}")
print(f"sklearn SC: {sc_sklearn:.4f}")

# 5. 三种指标对比
print("\n=== 5. 三种指标对比 ===")

def compare_metrics(X):
    """对比三种评估指标"""
    print("\n--- 三种指标对比 ---")
    
    k_range = range(1, 11)
    sse_list = []
    sc_list = []
    ch_list = []
    
    for k in k_range:
        kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
        labels = kmeans.fit_predict(X)
        
        sse = kmeans.inertia_
        sse_list.append(sse)
        
        if k >= 2:  # SC和CH至少需要2个簇
            sc = silhouette_score(X, labels)
            ch = calinski_harabasz_score(X, labels)
            sc_list.append((k, sc))
            ch_list.append((k, ch))
    
    # 绘制对比图
    fig, axes = plt.subplots(1, 3, figsize=(18, 5))
    
    # SSE
    axes[0].plot(k_range, sse_list, 'bo-')
    axes[0].set_xlabel('K值')
    axes[0].set_ylabel('SSE')
    axes[0].set_title('SSE（越小越好）')
    axes[0].grid(True)
    
    # SC
    sc_k = [x[0] for x in sc_list]
    sc_v = [x[1] for x in sc_list]
    axes[1].plot(sc_k, sc_v, 'ro-')
    axes[1].set_xlabel('K值')
    axes[1].set_ylabel('SC')
    axes[1].set_title('SC轮廓系数（越大越好）')
    axes[1].grid(True)
    
    # CH
    ch_k = [x[0] for x in ch_list]
    ch_v = [x[1] for x in ch_list]
    axes[2].plot(ch_k, ch_v, 'go-')
    axes[2].set_xlabel('K值')
    axes[2].set_ylabel('CH')
    axes[2].set_title('CH轮廓系数（越大越好）')
    axes[2].grid(True)
    
    plt.tight_layout()
    plt.savefig('metrics_comparison.png', dpi=100)
    print("三种指标对比图已保存")
    
    # 最佳K值
    best_k_sse = 4  # 肘部法
    best_k_sc = sc_k[np.argmax(sc_v)]
    best_k_ch = ch_k[np.argmax(ch_v)]
    
    print(f"\n最佳K值:")
    print(f"  SSE肘部法: K={best_k_sse}")
    print(f"  SC轮廓系数: K={best_k_sc} (SC={max(sc_v):.4f})")
    print(f"  CH轮廓系数: K={best_k_ch} (CH={max(ch_v):.4f})")

compare_metrics(X)

# 6. SC公式说明
print("\n=== 6. SC公式说明 ===")
print("""
SC轮廓系数公式:

SC = (b - a) / max(a, b)

参数说明:
- a: 样本i到簇内其他样本的平均距离（内聚）
  - a值越小越好（内聚越高）
  - 极端情况: a=0

- b: 样本i到最近其他簇所有样本的平均距离（分离）
  - b值越大越好（分离越高）
  - 极端情况: b=1

SC值范围: [-1, 1]
- SC=1: 最佳（a=0, b=1）
- SC=0: 一般（a=b）
- SC<0: 最差（a>b）
""")

# 7. 完整流程
def sc_ch_evaluation_pipeline():
    """SC和CH评估完整流程"""
    print("=" * 50)
    print("SC和CH评估 - 完整流程")
    print("=" * 50)
    
    # 1. 生成数据
    print("\n1. 生成数据")
    X, _ = make_blobs(n_samples=500, centers=4, cluster_std=0.5, random_state=42)
    
    # 2. 计算SC
    print("\n2. 计算SC轮廓系数")
    k_range = range(2, 11)
    sc_list = []
    for k in k_range:
        kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
        labels = kmeans.fit_predict(X)
        sc = silhouette_score(X, labels)
        sc_list.append(sc)
        print(f"   K={k}: SC={sc:.4f}")
    
    # 3. 计算CH
    print("\n3. 计算CH轮廓系数")
    ch_list = []
    for k in k_range:
        kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
        labels = kmeans.fit_predict(X)
        ch = calinski_harabasz_score(X, labels)
        ch_list.append(ch)
        print(f"   K={k}: CH={ch:.4f}")
    
    # 4. 选择最佳K值
    print("\n4. 选择最佳K值")
    best_k_sc = list(k_range)[np.argmax(sc_list)]
    best_k_ch = list(k_range)[np.argmax(ch_list)]
    print(f"   SC最佳K值: {best_k_sc}")
    print(f"   CH最佳K值: {best_k_ch}")
    
    print("\n" + "=" * 50)
    print("评估完成!")
    print("=" * 50)

sc_ch_evaluation_pipeline()