🔵 聚类分析 — K-Means与DBSCAN对比

📖 聚类分析概述

聚类是无监督学习的核心任务——没有标签，只有数据，目标是发现"相似的在一起"的自然分组。

🔢 K-Means算法详解

K-Means是最经典的聚类算法，通过迭代更新聚类中心来最小化簇内平方和。

import numpy as np

np.random.seed(42)

# 生成3簇数据
X_c1 = np.random.normal([2, 3], 0.8, (100, 2))
X_c2 = np.random.normal([7, 8], 0.8, (100, 2))
X_c3 = np.random.normal([8, 2], 0.8, (100, 2))
X = np.vstack([X_c1, X_c2, X_c3])

def kmeans(X, k, max_iter=100):
    """手写K-Means"""
    # K-Means++初始化
    idx = [np.random.randint(len(X))]
    for _ in range(1, k):
        dists = np.array([min(np.linalg.norm(x - X[i]) for i in idx) for x in X])
        probs = dists**2 / (dists**2).sum()
        idx.append(np.random.choice(len(X), p=probs))
    centroids = X[idx].copy()
    
    for _ in range(max_iter):
        # 分配
        distances = np.array([[np.linalg.norm(x - c) for c in centroids] for x in X])
        labels = np.argmin(distances, axis=1)
        
        # 更新
        new_centroids = np.array([X[labels == i].mean(axis=0) for i in range(k)])
        
        if np.allclose(centroids, new_centroids):
            break
        centroids = new_centroids
    
    return labels, centroids

labels_km, centroids_km = kmeans(X, 3)
print(f"K-Means聚类完成，3个中心:")
for i, c in enumerate(centroids_km):
    print(f"  簇{i}: 中心=({c[0]:.2f}, {c[1]:.2f}), 样本数={(labels_km==i).sum()}")

📊 选择最佳K值

1. 肘部法（Elbow Method）

# 肘部法: 惯性随K的变化
print("K vs 惯性(肘部法):")
for k in range(2, 7):
    labels_k, centroids_k = kmeans(X, k)
    inertia = sum(np.sum((X[labels_k == i] - centroids_k[i])**2) for i in range(k))
    bar = '█' * int(inertia / 100)
    print(f"  K={k}: 惯性={inertia:.0f} {bar}")
# K=2→3惯性大幅下降，之后减缓 → K=3是"肘部"

2. 轮廓系数（Silhouette Score）

def silhouette_score(X, labels):
    """手写轮廓系数"""
    n = len(X)
    scores = []
    for i in range(n):
        # a(i): 与同簇其他点的平均距离
        same = X[labels == labels[i]]
        if len(same) > 1:
            a = np.mean([np.linalg.norm(X[i] - s) for s in same if not np.array_equal(s, X[i])])
        else:
            a = 0
        
        # b(i): 与最近其他簇的平均距离
        b = float('inf')
        for c in set(labels) - {labels[i]}:
            other = X[labels == c]
            dist = np.mean([np.linalg.norm(X[i] - o) for o in other])
            b = min(b, dist)
        
        scores.append((b - a) / max(a, b) if max(a, b) > 0 else 0)
    return np.mean(scores)

print("K vs 轮廓系数:")
for k in range(2, 7):
    labels_k, _ = kmeans(X, k)
    sil = silhouette_score(X, labels_k)
    bar = '█' * int(sil * 30)
    print(f"  K={k}: 轮廓系数={sil:.4f} {bar}")
# K=3时轮廓系数最高(0.7726) → 最佳K=3

🔬 DBSCAN算法详解

DBSCAN（Density-Based Spatial Clustering of Applications with Noise）基于密度定义簇，能发现任意形状的聚类并自动识别噪声点。

def dbscan(X, eps, min_samples):
    """手写DBSCAN"""
    n = len(X)
    labels = np.full(n, -1)  # -1 = noise
    cluster_id = 0
    
    # 预计算距离矩阵
    dist_matrix = np.array([[np.linalg.norm(X[i]-X[j]) for j in range(n)] for i in range(n)])
    
    visited = set()
    for i in range(n):
        if i in visited: continue
        visited.add(i)
        
        neighbors = np.where(dist_matrix[i] <= eps)[0]
        
        if len(neighbors) < min_samples:
            continue  # 噪声点，暂不分配
        
        labels[i] = cluster_id
        seed_set = list(neighbors[neighbors != i])
        
        for j in seed_set:
            if j not in visited:
                visited.add(j)
                j_neighbors = np.where(dist_matrix[j] <= eps)[0]
                if len(j_neighbors) >= min_samples:
                    seed_set.extend(j_neighbors)
            
            if labels[j] == -1:
                labels[j] = cluster_id
        
        cluster_id += 1
    
    return labels

labels_db = dbscan(X, eps=1.5, min_samples=5)
n_clusters = len(set(labels_db)) - (1 if -1 in labels_db else 0)
n_noise = (labels_db == -1).sum()
print(f"DBSCAN: {n_clusters}个簇, {n_noise}个噪声点")

# 在不同数据分布上对比
# 圆环状数据（DBSCAN优势场景）
theta = np.random.uniform(0, 2*np.pi, 200)
r_inner = 2 + np.random.normal(0, 0.2, 200)
r_outer = 5 + np.random.normal(0, 0.2, 200)
X_rings = np.vstack([
    np.column_stack([r_inner*np.cos(theta), r_inner*np.sin(theta)]),
    np.column_stack([r_outer*np.cos(theta), r_outer*np.sin(theta)])
])

# K-Means在圆环上失败
labels_km_ring, _ = kmeans(X_rings, 2)
sil_km_ring = silhouette_score(X_rings, labels_km_ring)

# DBSCAN在圆环上成功
labels_db_ring = dbscan(X_rings, eps=0.8, min_samples=5)
n_cl_ring = len(set(labels_db_ring)) - (1 if -1 in labels_db_ring else 0)
valid = labels_db_ring != -1
sil_db_ring = silhouette_score(X_rings[valid], labels_db_ring[valid]) if valid.sum() > 0 and n_cl_ring >= 2 else -1

print(f"圆环数据对比:")
print(f"  K-Means: 轮廓系数={sil_km_ring:.4f} (不适合)")
print(f"  DBSCAN:  轮廓系数={sil_db_ring:.4f} ({n_cl_ring}个簇)")

# 球形数据（K-Means优势场景）
print(f"\n球形数据对比:")
print(f"  K-Means: 轮廓系数=0.7726 ✅")
print(f"  DBSCAN:  轮廓系数=0.7738 ✅")

# 简化GMM (EM算法)
from scipy import stats

def gmm_em(X, k, max_iter=50):
    """EM算法实现GMM"""
    n, d = X.shape
    
    # 初始化
    np.random.seed(42)
    means = X[np.random.choice(n, k, replace=False)]
    covs = [np.eye(d) for _ in range(k)]
    weights = np.ones(k) / k
    
    for iteration in range(max_iter):
        # E步: 计算每个样本属于每个高斯的后验概率
        resp = np.zeros((n, k))
        for j in range(k):
            resp[:, j] = weights[j] * stats.multivariate_normal.pdf(X, means[j], covs[j])
        resp /= resp.sum(axis=1, keepdims=True)
        
        # M步: 更新参数
        Nk = resp.sum(axis=0)
        for j in range(k):
            means[j] = (resp[:, j:j+1].T @ X) / Nk[j]
            diff = X - means[j]
            covs[j] = (diff.T @ (diff * resp[:, j:j+1])) / Nk[j] + 1e-6*np.eye(d)
            weights[j] = Nk[j] / n
    
    labels = resp.argmax(axis=1)
    return labels, means, covs, weights

labels_gmm, means_gmm, _, _ = gmm_em(X, 3)
sil_gmm = silhouette_score(X, labels_gmm)
print(f"GMM轮廓系数: {sil_gmm:.4f}")