无监督学习与聚类实战

和监督学习不同，无监督学习没有”正确答案”。想象一下，给你一堆没有标签的客户数据，让你自己找出其中的规律——哪些客户比较相似？数据中有哪些隐藏的结构？这就是无监督学习要解决的问题。

无监督学习主要包括两大类任务：聚类（把相似的东西分到一起）和降维（把高维数据压缩到低维，同时保留关键信息）。还有一个重要应用是异常检测——找出那些”不正常”的数据点。让我带你深入了解这些技术。

聚类的核心问题

聚类看似简单——就是把相似的东西归为一类。但实际操作中，有几个关键问题需要解决：

第一，“相似”如何定义？这涉及距离度量。欧氏距离最常用，但有时候余弦相似度更适合文本数据，马氏距离适合特征尺度差异大的数据。

第二，聚成几类？有些算法需要你事先指定K值，有些可以自动推断。选多了会过拟合，选少了会欠拟合。

第三，算法的复杂度。K-Means很快，但只能找到球形簇；层次聚类能发现任意形状的簇，但计算量大；DBSCAN能处理噪声，但参数敏感。

K-Means：简单但强大的经典算法

K-Means大概是名气最大的聚类算法了。它的思想很直观：先把数据分成K个簇，每个簇有个中心点（质心），然后迭代地做两件事——把每个点分配给离它最近的质心，重新计算每个簇的质心。迭代直到质心不再变化（或变化很小）。

算法细节与收敛性

K-Means的数学形式化是：找到K个质心，使得每个点到其所属簇质心的距离平方和最小。这个目标函数是非凸的，所以K-Means可能会收敛到局部最优解。初始化对结果影响很大，运气不好可能收敛到很差的局部最优。

解决初始化问题的主流方法是K-Means++。它的做法是：第一个质心随机选，后续质心按照到已有质心距离的平方成正比的概率来选。这样可以让初始质心分散开来，大幅提高找到好解的概率。scikit-learn默认使用K-Means++初始化。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs, make_moons, make_circles
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score, calinski_harabasz_score, davies_bouldin_score
 
# 生成演示用的数据
# 1. 分离良好的球形簇
X_blobs, y_blobs = make_blobs(
    n_samples=500, 
    centers=4, 
    cluster_std=0.6, 
    random_state=42
)
 
# 2. 半月形数据
X_moons, y_moons = make_moons(n_samples=300, noise=0.1, random_state=42)
 
# 3. 环形数据
X_circles, y_circles = make_circles(n_samples=300, noise=0.05, factor=0.5, random_state=42)
 
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
 
axes[0].scatter(X_blobs[:, 0], X_blobs[:, 1], c=y_blobs, cmap='viridis', alpha=0.7)
axes[0].set_title('Blob Data (K-Means Friendly)')
 
axes[1].scatter(X_moons[:, 0], X_moons[:, 1], c=y_moons, cmap='viridis', alpha=0.7)
axes[1].set_title('Moon Data (K-Means Struggles)')
 
axes[2].scatter(X_circles[:, 0], X_circles[:, 1], c=y_circles, cmap='viridis', alpha=0.7)
axes[2].set_title('Circle Data (K-Means Struggles)')
 
plt.tight_layout()
plt.savefig('clustering_data_types.png', dpi=150)
 
# K-Means实战：对球形数据效果很好
scaler = StandardScaler()
X_blobs_scaled = scaler.fit_transform(X_blobs)
 
# 使用肘部法则确定最佳K值
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_blobs_scaled)
    inertias.append(kmeans.inertia_)
    silhouette_scores.append(silhouette_score(X_blobs_scaled, kmeans.labels_))
 
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
 
axes[0].plot(K_range, inertias, 'bo-')
axes[0].set_xlabel('Number of Clusters (K)')
axes[0].set_ylabel('Inertia (Within-cluster SSE)')
axes[0].set_title('Elbow Method')
axes[0].axvline(x=4, color='r', linestyle='--', label='Optimal K=4')
axes[0].legend()
 
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 Score vs K')
axes[1].axvline(x=4, color='r', linestyle='--', label='Optimal K=4')
axes[1].legend()
 
plt.tight_layout()
plt.savefig('kmeans_k_selection.png', dpi=150)
 
# 可视化K-Means结果
kmeans = KMeans(n_clusters=4, random_state=42, n_init=10)
y_pred = kmeans.fit_predict(X_blobs_scaled)
 
plt.figure(figsize=(8, 6))
plt.scatter(X_blobs_scaled[:, 0], X_blobs_scaled[:, 1], c=y_pred, cmap='viridis', alpha=0.7)
plt.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1], 
            c='red', marker='X', s=200, edgecolors='black', linewidths=2,
            label='Centroids')
plt.title('K-Means Clustering Result')
plt.legend()
plt.savefig('kmeans_result.png', dpi=150)

K-Means的局限性

K-Means有很多局限性需要你了解。首先，它假设簇是球形、大小相似的，这在很多实际数据上不成立。其次，它对噪声和离群点敏感——一个极端的离群点会严重影响质心的位置。第三，它需要预先指定K值，而且对初始化敏感。

对于半月形和环形数据，K-Means表现很差，因为它试图把数据切割成若干个球形区域。下一节介绍的DBSCAN和层次聚类能处理这类非球形簇。

层次聚类：发现数据的层级结构

层次聚类（Hierarchical Clustering）是一种不需要预先指定K值的聚类方法。它通过计算数据点之间的相似度，逐步合并或分裂簇，形成一个树状图（Dendrogram），展示数据的层级结构。

层次聚类有两种策略：凝聚式（Agglomerative）和分裂式（Divisive）。凝聚式从每个点作为一个单独的簇开始，逐步合并最相似的簇；分裂式则相反，从所有点作为一个簇开始，逐步分裂。实际应用中，凝聚式更常用。

链接方法

凝聚式层次聚类的关键是如何定义两个簇之间的”距离”。常见的链接方法有：

单链接（Single Linkage）：两个簇的距离定义为最近两个点之间的距离。优点是能处理任意形状的簇，缺点是容易产生”链式效应”——两个簇之间有一条点组成的链，就会被合并，即使整体形状很奇怪。

完全链接（Complete Linkage）：两个簇的距离定义为最远两个点之间的距离。倾向于产生紧凑、大小相近的簇，但可能分割真正的自然簇。

平均链接（Average Linkage）：两个簇的距离定义为所有点对距离的平均值。是比较中庸的选择，效果通常不错。

Ward链接：最小化合并导致的方差增加。类似于K-Means的目标，能产生紧凑的簇，但对噪声敏感。

from sklearn.cluster import AgglomerativeClustering
from scipy.cluster.hierarchy import dendrogram, linkage
from scipy.spatial.distance import pdist
 
# 层次聚类可视化
fig, axes = plt.subplots(1, 2, figsize=(15, 6))
 
# 计算层次聚类的链接矩阵
Z = linkage(X_blobs_scaled, method='ward')
 
# 绘制树状图
dendrogram(Z, ax=axes[0], truncate_mode='lastp', p=30, 
           leaf_rotation=90, leaf_font_size=8, show_contracted=True)
axes[0].set_title('Hierarchical Clustering Dendrogram (Ward Linkage)')
axes[0].set_xlabel('Sample Index / Cluster Size')
axes[0].set_ylabel('Distance')
 
# 使用不同的链接方法对比
methods = ['single', 'complete', 'average', 'ward']
colors = ['red', 'green', 'blue', 'purple']
 
for method, color in zip(methods, colors):
    Z = linkage(X_blobs_scaled, method=method)
    dendrogram(Z, ax=axes[1], truncate_mode='lastp', p=15, 
               color_threshold=0.7*max(Z[:,2]), no_labels=True)
    axes[1].axhline(y=7, color='gray', linestyle='--', alpha=0.5)
 
axes[1].set_title('Different Linkage Methods Comparison')
axes[1].set_xlabel('Samples')
axes[1].set_ylabel('Distance')
 
plt.tight_layout()
plt.savefig('hierarchical_clustering.png', dpi=150)
 
# 凝聚式层次聚类实战
agg_clustering = AgglomerativeClustering(
    n_clusters=4,
    linkage='ward'  # 也可以是 'single', 'complete', 'average'
)
 
y_pred_agg = agg_clustering.fit_predict(X_blobs_scaled)
 
plt.figure(figsize=(8, 6))
plt.scatter(X_blobs_scaled[:, 0], X_blobs_scaled[:, 1], c=y_pred_agg, cmap='viridis', alpha=0.7)
plt.title('Agglomerative Clustering Result (Ward)')
plt.savefig('agglomerative_result.png', dpi=150)

树状图的使用

树状图是层次聚类的精华。你可以通过观察树状图判断数据应该分成几类：找一条”最长”的垂直线段，这个线段不被任何水平线穿过，代表一次”断裂”，这条线对应的聚类数就是合适的K值。

在scikit-learn中使用层次聚类，你可以不指定n_clusters，而是通过distance_threshold参数来控制聚类：只合并距离小于这个阈值的簇。这样可以自动发现合适的簇数。

# 通过距离阈值自动确定簇数
agg_auto = AgglomerativeClustering(
    distance_threshold=1.5,
    n_clusters=None,
    linkage='ward'
)
 
y_pred_auto = agg_auto.fit_predict(X_blobs_scaled)
n_clusters_auto = len(np.unique(y_pred_auto))
 
print(f"Auto-detected number of clusters: {n_clusters_auto}")
 
plt.figure(figsize=(8, 6))
plt.scatter(X_blobs_scaled[:, 0], X_blobs_scaled[:, 1], c=y_pred_auto, cmap='viridis', alpha=0.7)
plt.title(f'Auto Clustering Result (Distance Threshold=1.5, {n_clusters_auto} clusters)')
plt.savefig('agglomerative_auto.png', dpi=150)

DBSCAN：发现任意形状的簇

DBSCAN（Density-Based Spatial Clustering of Applications with Noise）是密度聚类的代表算法。它的核心思想是把数据分成三类：核心点（在稠密区域内的点）、边界点（在稠密区域边缘的点）、噪声点（远离任何稠密区域的点）。

算法原理

DBSCAN有两个关键参数：eps（ε，邻域半径）和min_samples（最小样本数）。一个点是核心点，当且仅当它的ε邻域内至少有min_samples个点（包括自己）。两个点是直接密度可达的，当且仅当一个点在另一个点的ε邻域内。一条密度可达的链连接两个点，意味着它们可以通过直接密度可达的关系连接起来。

基于这些定义，DBSCAN的簇是密度可达的最大集合。噪声点不属于任何簇。

from sklearn.cluster import DBSCAN
from sklearn.neighbors import NearestNeighbors
 
# DBSCAN参数选择：使用k-距离图
# 对于每个点，计算它到第k个最近邻的距离，然后画图找拐点
k = 2 * X_blobs_scaled.shape[1] - 1  # 经验法则：2*维度-1
 
nn = NearestNeighbors(n_neighbors=k)
nn.fit(X_blobs_scaled)
distances, indices = nn.kneighbors(X_blobs_scaled)
 
# 取第k个最近邻的距离（即k-距离）
k_distances = distances[:, k-1]
k_distances = np.sort(k_distances)[::-1]
 
plt.figure(figsize=(10, 6))
plt.plot(range(len(k_distances)), k_distances, 'b-')
plt.xlabel('Points (sorted by distance)')
plt.ylabel(f'{k}-Distance')
plt.title('K-Distance Graph for DBSCAN eps Selection')
plt.grid(True, alpha=0.3)
 
# 找一个合理的eps值（拐点位置）
# 通常选择k-距离急剧上升之前的位置
eps_candidates = [0.3, 0.5, 0.7, 1.0]
for eps in eps_candidates:
    n_noise = np.sum(k_distances > eps) / len(k_distances) * 100
    print(f"eps={eps}: ~{n_noise:.1f}% points would be noise")

DBSCAN vs K-Means

DBSCAN最大的优势是可以发现任意形状的簇，包括半月形、环形、螺旋形等K-Means无法处理的形状。另外，DBSCAN不需要指定簇的数量，可以自动检测噪声点。

但DBSCAN也有缺点：它对参数eps和min_samples的选择很敏感，不同的参数可能导致截然不同的结果。在高维数据上，由于”维度诅咒”，密度概念变得模糊，DBSCAN的效果会下降。

# 对比DBSCAN在不同数据上的表现
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
 
datasets = [
    ('Blobs', X_blobs_scaled),
    ('Moons', X_moons),
    ('Circles', X_circles)
]
 
dbscan_params = [
    {'eps': 0.4, 'min_samples': 5},
    {'eps': 0.3, 'min_samples': 5},
    {'eps': 0.3, 'min_samples': 5}
]
 
for row, ((name, X), params) in enumerate(zip(datasets, dbscan_params)):
    # K-Means结果
    kmeans = KMeans(n_clusters=2 if name != 'Blobs' else 4, random_state=42, n_init=10)
    y_kmeans = kmeans.fit_predict(X)
    axes[row, 0].scatter(X[:, 0], X[:, 1], c=y_kmeans, cmap='viridis', alpha=0.7)
    axes[row, 0].set_title(f'{name} - K-Means')
    
    # DBSCAN结果
    dbscan = DBSCAN(eps=params['eps'], min_samples=params['min_samples'])
    y_dbscan = dbscan.fit_predict(X)
    axes[row, 1].scatter(X[:, 0], X[:, 1], c=y_dbscan, cmap='viridis', alpha=0.7)
    n_clusters_dbscan = len(set(y_dbscan)) - (1 if -1 in y_dbscan else 0)
    n_noise = list(y_dbscan).count(-1)
    axes[row, 1].set_title(f'{name} - DBSCAN\nClusters: {n_clusters_dbscan}, Noise: {n_noise}')
    
    # 层次聚类结果
    agg = AgglomerativeClustering(n_clusters=2 if name != 'Blobs' else 4)
    y_agg = agg.fit_predict(X)
    axes[row, 2].scatter(X[:, 0], X[:, 1], c=y_agg, cmap='viridis', alpha=0.7)
    axes[row, 2].set_title(f'{name} - Hierarchical')
 
plt.tight_layout()
plt.savefig('clustering_algorithm_comparison.png', dpi=150)

GMM：软聚类的概率方法

GMM（Gaussian Mixture Model，高斯混合模型）是一种基于概率的聚类方法。与K-Means的”硬分配”不同，GMM给出的是”软分配”——每个点属于每个簇都有一个概率值。

数学原理

GMM假设数据是由K个高斯分布混合生成的。每个高斯分布有自己的均值向量和协方差矩阵。GMM的目标是找到这K个高斯分布的参数，使得生成这些数据的概率最大。

GMM的训练使用EM（Expectation-Maximization）算法。E步计算每个点属于每个簇的后验概率（软分配），M步根据这些概率重新估计高斯分布的参数。迭代直到收敛。

GMM的协方差矩阵可以有不同的约束：有的簇是球形的（只有方差），有的可以拉伸成任意形状的椭球（完整协方差矩阵）。不同的约束影响模型的复杂度和灵活性。

from sklearn.mixture import GaussianMixture
 
# GMM聚类
gmm = GaussianMixture(
    n_components=4,
    covariance_type='full',  # 'full', 'tied', 'diag', 'spherical'
    random_state=42,
    n_init=5  # 多次初始化取最优
)
 
gmm.fit(X_blobs_scaled)
 
# 获取软分配概率
probs = gmm.predict_proba(X_blobs_scaled)
y_pred_gmm = gmm.predict(X_blobs_scaled)
 
# 可视化GMM结果（用概率作为颜色强度）
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
 
# 左图：硬分配
axes[0].scatter(X_blobs_scaled[:, 0], X_blobs_scaled[:, 1], c=y_pred_gmm, cmap='viridis', alpha=0.7)
axes[0].set_title('GMM Hard Assignment')
 
# 右图：软分配（每个点用4种颜色的混合表示）
# 创建一个混合图：用每个簇的颜色，按概率混合
rgb_colors = plt.cm.viridis(np.linspace(0, 1, 4))
final_colors = np.zeros((len(X_blobs_scaled), 3))
for i in range(4):
    final_colors += probs[:, i:i+1] * rgb_colors[i][:3]
 
axes[1].scatter(X_blobs_scaled[:, 0], X_blobs_scaled[:, 1], c=final_colors, alpha=0.7)
axes[1].set_title('GMM Soft Assignment (Probability Colors)')
 
plt.tight_layout()
plt.savefig('gmm_clustering.png', dpi=150)
 
# 对比不同协方差类型
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
 
cov_types = ['full', 'tied', 'diag', 'spherical']
titles = ['Full (每个簇独立协方差)', 'Tied (所有簇共享协方差)', 
          'Diag (协方差矩阵对角化)', 'Spherical (球形协方差)']
 
for ax, cov_type, title in zip(axes.flatten(), cov_types, titles):
    gmm = GaussianMixture(n_components=4, covariance_type=cov_type, random_state=42, n_init=3)
    y_pred = gmm.fit_predict(X_blobs_scaled)
    
    ax.scatter(X_blobs_scaled[:, 0], X_blobs_scaled[:, 1], c=y_pred, cmap='viridis', alpha=0.7)
    ax.set_title(title)
    
    # 画出高斯分布的等概率轮廓
    from matplotlib.patches import Ellipse
    
    for i in range(4):
        if cov_type == 'spherical':
            cov = np.eye(2) * gmm.covariances_[i]
        elif cov_type == 'diag':
            cov = np.diag(gmm.covariances_[i])
        elif cov_type == 'tied':
            cov = gmm.covariances_
        else:
            cov = gmm.covariances_[i]
        
        mean = gmm.means_[i]
        
        # 绘制2σ椭圆
        eigenvalues, eigenvectors = np.linalg.eigh(cov)
        order = eigenvalues.argsort()[::-1]
        eigenvalues, eigenvectors = eigenvalues[order], eigenvectors[:, order]
        
        theta = np.degrees(np.arctan2(eigenvectors[1, 0], eigenvectors[0, 0]))
        width, height = 2 * 2 * np.sqrt(eigenvalues)  # 2σ
        
        ellipse = Ellipse(xy=mean, width=width, height=height, angle=theta,
                         fill=False, color='red', linewidth=2)
        ax.add_patch(ellipse)
 
plt.tight_layout()
plt.savefig('gmm_covariance_types.png', dpi=150)

聚类评估指标

聚类没有”正确答案”，但我们需要一些指标来比较不同聚类结果的好坏，或者确定最佳的簇数。

内部评估指标

内部指标只使用数据本身，不需要真实标签。

**轮廓系数（Silhouette Score）**是最常用的内部指标之一。对于每个点i，计算：

• a(i)：点i与同簇其他点的平均距离（簇内距离）
• b(i)：点i与最近邻簇的点的平均距离（簇间距离）
• 轮廓系数 s(i) = (b(i) - a(i)) / max(a(i), b(i))

s(i)的取值范围是[-1, 1]。越接近1表示聚类越紧凑、簇间分离越好；越接近-1表示点被分到了错误的簇；接近0表示簇之间有重叠。轮廓系数的平均值反映了整体聚类质量。

CH指数（Calinski-Harabasz Index）：衡量簇间离散度与簇内离散度的比值。越大越好。

DBI（Davies-Bouldin Index）：衡量每个簇与其最相似簇的相似度。越小越好。

from sklearn.metrics import silhouette_score, calinski_harabasz_score, davies_bouldin_score
 
# 计算各种评估指标
def evaluate_clustering(X, labels, name):
    n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
    
    # 排除噪声点（label=-1）计算指标
    mask = labels != -1
    X_filtered = X[mask]
    labels_filtered = labels[mask]
    
    if n_clusters < 2:
        print(f"{name}: 只有1个簇或全部是噪声")
        return
    
    silhouette = silhouette_score(X_filtered, labels_filtered)
    ch = calinski_harabasz_score(X_filtered, labels_filtered)
    dbi = davies_bouldin_score(X_filtered, labels_filtered)
    
    print(f"{name}:")
    print(f"  簇数: {n_clusters}")
    print(f"  轮廓系数: {silhouette:.4f} (越大越好)")
    print(f"  CH指数: {ch:.2f} (越大越好)")
    print(f"  DBI: {dbi:.4f} (越小越好)")
    print()
 
# 对不同算法计算指标
evaluate_clustering(X_blobs_scaled, kmeans.labels_, "K-Means")
evaluate_clustering(X_blobs_scaled, agg_clustering.labels_, "层次聚类")
 
# DBSCAN对非球形数据效果更好
dbscan_blobs = DBSCAN(eps=0.4, min_samples=5)
dbscan_moons = DBSCAN(eps=0.3, min_samples=5)
dbscan_circles = DBSCAN(eps=0.3, min_samples=5)
 
y_dbscan_blobs = dbscan_blobs.fit_predict(X_blobs_scaled)
y_dbscan_moons = dbscan_moons.fit_predict(X_moons)
y_dbscan_circles = dbscan_circles.fit_predict(X_circles)
 
evaluate_clustering(X_blobs_scaled, y_dbscan_blobs, "DBSCAN (Blobs)")
evaluate_clustering(X_moons, y_dbscan_moons, "DBSCAN (Moons)")
evaluate_clustering(X_circles, y_dbscan_circles, "DBSCAN (Circles)")

降维技术：从高维到低维

真实世界的数据往往是高维的。图像可能有上万个像素，文本有上万个词汇。高维数据不仅计算量大，还面临”维度诅咒”——在高维空间中，数据变得极其稀疏，距离的概念变得模糊。

降维是把高维数据压缩到低维空间的技术。好的降维方法会保留数据的关键结构和信息，比如相似的点在降维后仍然相近。

PCA：主成分分析

PCA是最经典的降维方法。它的核心思想是找到数据中方差最大的方向（主成分），然后把所有数据投影到这些主成分上。

PCA的数学过程是：对数据做去中心化，计算协方差矩阵，求协方差矩阵的特征值和特征向量，按特征值大小排序，取前K个特征向量作为新的基，把数据投影到这K个基上。特征值的大小表示对应主成分保留了多少信息（方差）。

PCA的缺点是它是线性方法，只能捕捉线性结构。对于非线性结构的数据（如瑞士卷曲面），PCA会”压扁”数据，而不是”展开”它。

from sklearn.decomposition import PCA
from sklearn.datasets import load_iris, load_digits
import seaborn as sns
 
# PCA在鸢尾花数据集上的演示
iris = load_iris()
X_iris = iris.data
y_iris = iris.target
 
# PCA降到2维可视化
pca = PCA(n_components=2)
X_iris_pca = pca.fit_transform(X_iris)
 
print(f"原始维度: {X_iris.shape}")
print(f"降维后维度: {X_iris_pca.shape}")
print(f"主成分解释的方差比例: {pca.explained_variance_ratio_}")
print(f"累计解释方差: {sum(pca.explained_variance_ratio_):.3f}")
 
plt.figure(figsize=(10, 6))
colors = ['red', 'green', 'blue']
for i, (label, color) in enumerate(zip(iris.target_names, colors)):
    mask = y_iris == i
    plt.scatter(X_iris_pca[mask, 0], X_iris_pca[mask, 1], 
                c=color, label=label, alpha=0.7, s=50)
plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.1%} variance)')
plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.1%} variance)')
plt.title('PCA on Iris Dataset')
plt.legend()
plt.savefig('pca_iris.png', dpi=150)
 
# 选择合适的主成分数量（累积解释方差>95%）
pca_full = PCA().fit(X_iris)
cumulative_var = np.cumsum(pca_full.explained_variance_ratio_)
 
plt.figure(figsize=(10, 6))
plt.plot(range(1, len(cumulative_var)+1), cumulative_var, 'bo-')
plt.axhline(y=0.95, color='r', linestyle='--', label='95% threshold')
plt.axhline(y=0.99, color='g', linestyle='--', label='99% threshold')
plt.xlabel('Number of Components')
plt.ylabel('Cumulative Explained Variance')
plt.title('PCA Explained Variance Analysis')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('pca_variance_analysis.png', dpi=150)

t-SNE：非线性降维的可视化利器

t-SNE（t-distributed Stochastic Neighbor Embedding）是专门为可视化设计的非线性降维方法。它特别擅长把高维数据投影到2D或3D空间，同时保留局部结构——相似的点在降维后仍然是邻居。

t-SNE的工作原理是：在高维空间计算每个点对其他点的相似度（用正态分布的概率），在低维空间也做类似计算（用t分布的概率），然后最小化两个分布的KL散度。

t-SNE的主要缺点是计算复杂度高，O(n²)的时间复杂度让它难以处理大数据集；另外它的结果具有随机性，每次运行可能略有不同；最重要的是，t-SNE只适合可视化，不适合作为下游任务的特征。

from sklearn.manifold import TSNE
 
# t-SNE在手写数字数据集上的演示
digits = load_digits()
X_digits = digits.data
y_digits = digits.target
 
print(f"手写数字数据集: {X_digits.shape[0]} 样本, {X_digits.shape[1]} 特征")
 
# 先用PCA降到50维加速t-SNE
pca = PCA(n_components=50, random_state=42)
X_digits_pca = pca.fit_transform(X_digits)
 
# t-SNE降维到2维
tsne = TSNE(
    n_components=2,
    perplexity=30,      # 通常在5-50之间，越大越关注全局结构
    learning_rate='auto',
    init='pca',
    random_state=42,
    n_jobs=-1
)
 
X_digits_tsne = tsne.fit_transform(X_digits_pca)
 
plt.figure(figsize=(12, 10))
scatter = plt.scatter(X_digits_tsne[:, 0], X_digits_tsne[:, 1], 
                      c=y_digits, cmap='tab10', alpha=0.6, s=10)
plt.colorbar(scatter, ticks=range(10))
plt.title('t-SNE on Handwritten Digits (64D -> 2D)')
plt.savefig('tsne_digits.png', dpi=150)
 
# 对比PCA和t-SNE在同一数据上的效果
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
 
pca_2d = PCA(n_components=2)
X_digits_pca_2d = pca_2d.fit_transform(X_digits)
 
axes[0].scatter(X_digits_pca_2d[:, 0], X_digits_pca_2d[:, 1], 
                c=y_digits, cmap='tab10', alpha=0.6, s=10)
axes[0].set_title('PCA (线性降维)')
 
axes[1].scatter(X_digits_tsne[:, 0], X_digits_tsne[:, 1], 
                c=y_digits, cmap='tab10', alpha=0.6, s=10)
axes[1].set_title('t-SNE (非线性降维)')
 
plt.savefig('pca_vs_tsne.png', dpi=150)

UMAP：更快更好的选择

UMAP（Uniform Manifold Approximation and Projection）是近年来很火的降维方法。相比t-SNE，它有几个优势：速度更快（用近似最近邻算法）、保留更多全局结构、结果更稳定、支持降维到更高维度。

UMAP的理论基础是黎曼几何和代数拓扑，它假设数据分布在一个流形上，目标是找到这个流形的低维表示。从实际效果看，UMAP在大多数任务上都能媲美甚至超越t-SNE。

# UMAP的使用（需要安装umap-learn）
try:
    import umap
    
    reducer = umap.UMAP(
        n_components=2,
        n_neighbors=15,      # 局部邻居数，越大越关注全局
        min_dist=0.1,        # 最小距离，越小嵌入越紧密
        metric='euclidean',
        random_state=42
    )
    
    X_digits_umap = reducer.fit_transform(X_digits_pca)
    
    plt.figure(figsize=(12, 10))
    scatter = plt.scatter(X_digits_umap[:, 0], X_digits_umap[:, 1], 
                          c=y_digits, cmap='tab10', alpha=0.6, s=10)
    plt.colorbar(scatter, ticks=range(10))
    plt.title('UMAP on Handwritten Digits')
    plt.savefig('umap_digits.png', dpi=150)
    
    print("UMAP可视化完成")
    
except ImportError:
    print("UMAP未安装，跳过UMAP演示。请运行: pip install umap-learn")

异常检测：找出”不正常”的数据

异常检测（Anomaly Detection）是另一个重要的无监督学习任务。它的目标是找出那些偏离正常模式的”异常”数据点。应用场景包括：金融欺诈检测、网络入侵检测、设备故障预警、质量控制等。

Isolation Forest

Isolation Forest是一种基于随机森林思想的异常检测算法。它的核心假设是：异常点更容易被”隔离”——在随机分割过程中，异常点通常会更快地被分离出来，因为它们通常位于稀疏区域。

具体做法是：构建多棵随机决策树，每棵树随机选择一个特征和一个分割值。如果一个点在少数几次分割后就被单独隔离出来（落在叶子节点），它就更可能是异常点。最终的异常分数基于点在所有树中被隔离所需的平均分割次数。

from sklearn.ensemble import IsolationForest
from sklearn.neighbors import LocalOutlierFactor
from sklearn.svm import OneClassSVM
 
# 创建包含异常点的数据
np.random.seed(42)
n_normal = 500
n_outlier = 50
 
# 正常数据：从多元正态分布采样
mean = [0, 0]
cov = [[1, 0], [0, 1]]
X_normal = np.random.multivariate_normal(mean, cov, n_normal)
 
# 异常数据：从远离正常区域的地方采样
X_outlier = np.random.uniform(low=[-5, -5], high=[-3, -3], size=(n_outlier // 4, 2))
X_outlier = np.vstack([X_outlier, np.random.uniform(low=[3, 3], high=[5, 5], size=(n_outlier // 4, 2))])
X_outlier = np.vstack([X_outlier, np.random.uniform(low=[-5, 3], high=[-3, 5], size=(n_outlier // 4, 2))])
X_outlier = np.vstack([X_outlier, np.random.uniform(low=[3, -5], high=[5, -3], size=(n_outlier // 4, 2))])
 
X_combined = np.vstack([X_normal, X_outlier])
y_true = np.concatenate([np.zeros(n_normal), np.ones(n_outlier)])
 
# 异常检测
fig, axes = plt.subplots(2, 2, figsize=(14, 14))
 
# 1. Isolation Forest
iso_forest = IsolationForest(
    n_estimators=100,
    contamination=0.1,  # 假设10%的异常点
    random_state=42
)
y_pred_iso = iso_forest.fit_predict(X_combined)
y_score_iso = -iso_forest.decision_function(X_combined)  # 异常分数（越大越异常）
 
axes[0, 0].scatter(X_combined[:, 0], X_combined[:, 1], c=y_pred_iso, cmap='coolwarm', alpha=0.6)
axes[0, 0].set_title('Isolation Forest')
 
# 2. Local Outlier Factor
lof = LocalOutlierFactor(
    n_neighbors=20,
    contamination=0.1
)
y_pred_lof = lof.fit_predict(X_combined)
 
axes[0, 1].scatter(X_combined[:, 0], X_combined[:, 1], c=y_pred_lof, cmap='coolwarm', alpha=0.6)
axes[0, 1].set_title('Local Outlier Factor (LOF)')
 
# 3. One-Class SVM
ocsvm = OneClassSVM(
    kernel='rbf',
    gamma='auto',
    nu=0.1  # 异常比例的上界
)
y_pred_ocsvm = ocsvm.fit_predict(X_combined)
 
axes[1, 0].scatter(X_combined[:, 0], X_combined[:, 1], c=y_pred_ocsvm, cmap='coolwarm', alpha=0.6)
axes[1, 0].set_title('One-Class SVM')
 
# 4. 用异常分数热力图可视化
axes[1, 1].scatter(X_combined[:, 0], X_combined[:, 1], c=y_score_iso, cmap='YlOrRd', alpha=0.6)
axes[1, 1].set_title('Isolation Forest Anomaly Scores')
plt.colorbar(axes[1, 1].collections[0], ax=axes[1, 1], label='Anomaly Score')
 
plt.tight_layout()
plt.savefig('anomaly_detection_comparison.png', dpi=150)
 
# 评估异常检测效果
from sklearn.metrics import precision_score, recall_score, f1_score
 
for name, y_pred in [('Isolation Forest', y_pred_iso), 
                     ('LOF', y_pred_lof), 
                     ('One-Class SVM', y_pred_ocsvm)]:
    # 转换为1（异常）和-1（正常）的标签格式
    precision = precision_score(y_true, (y_pred == -1).astype(int))
    recall = recall_score(y_true, (y_pred == -1).astype(int))
    f1 = f1_score(y_true, (y_pred == -1).astype(int))
    print(f"{name}: Precision={precision:.3f}, Recall={recall:.3f}, F1={f1:.3f}")

Local Outlier Factor (LOF)

LOF是一种基于密度的异常检测方法。它的核心思想是比较一个点与其邻居的密度：如果一个点的密度显著低于其邻居的密度，它就是异常点。

具体计算：对于每个点，计算它到第K个最近邻的距离（K-distance），然后计算”局部可达密度”——点是它K个邻居的K距离的平均值的倒数。最后，LOF分数是该点的局部可达密度与邻居的局部可达密度之比的平均值。LOF接近1表示正常点，远大于1表示异常点。

LOF的优势是能检测局部异常——在稠密簇中的稀疏点也会被识别出来。但它对参数K的选择敏感，在高维数据上效果下降。

实战案例：客户分群

最后用一个完整的实战案例把知识串起来。假设你是一个电商公司的数据分析师，需要对客户进行分群，以便制定针对性的营销策略。

import pandas as pd
from sklearn.cluster import KMeans, AgglomerativeClustering, DBSCAN
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn.decomposition import PCA
from sklearn.metrics import silhouette_score
import warnings
warnings.filterwarnings('ignore')
 
# 模拟电商客户数据
np.random.seed(42)
n_customers = 1000
 
customers = pd.DataFrame({
    'customer_id': range(n_customers),
    'age': np.random.randint(18, 70, n_customers),
    'annual_income': np.random.lognormal(10.5, 0.5, n_customers),
    'spending_score': np.random.lognormal(4, 0.5, n_customers),
    'purchase_frequency': np.random.poisson(10, n_customers),
    'avg_basket_size': np.random.lognormal(3.5, 0.8, n_customers),
    'days_since_last_purchase': np.random.exponential(30, n_customers) + 1,
    'num_categories': np.random.poisson(3, n_customers) + 1,
    'discount_sensitivity': np.random.beta(2, 5, n_customers)  # 0-1，越高越敏感
})
 
# 特征工程
features = ['age', 'annual_income', 'spending_score', 'purchase_frequency', 
            'avg_basket_size', 'days_since_last_purchase', 'num_categories', 
            'discount_sensitivity']
 
X_customers = customers[features].values
 
# 标准化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_customers)
 
# 降维到2维用于可视化
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)
 
print(f"PCA累计解释方差: {sum(pca.explained_variance_ratio_):.2%}")
 
# 使用肘部法则和轮廓系数确定最佳K
K_range = range(2, 11)
inertias = []
silhouettes = []
 
for k in K_range:
    kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
    kmeans.fit(X_scaled)
    inertias.append(kmeans.inertia_)
    silhouettes.append(silhouette_score(X_scaled, kmeans.labels_))
 
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
 
axes[0].plot(K_range, inertias, 'bo-')
axes[0].set_xlabel('K')
axes[0].set_ylabel('Inertia')
axes[0].set_title('Elbow Method')
 
axes[1].plot(K_range, silhouettes, 'go-')
axes[1].set_xlabel('K')
axes[1].set_ylabel('Silhouette Score')
axes[1].set_title('Silhouette Score vs K')
optimal_k = K_range[np.argmax(silhouettes)]
axes[1].axvline(x=optimal_k, color='r', linestyle='--', label=f'Optimal K={optimal_k}')
axes[1].legend()
 
plt.tight_layout()
plt.savefig('customer_clustering_k_selection.png', dpi=150)
 
# 选定K=5进行聚类
optimal_k = 5
kmeans_final = KMeans(n_clusters=optimal_k, random_state=42, n_init=10)
customers['cluster'] = kmeans_final.fit_predict(X_scaled)
 
# 分析各簇的特征
cluster_summary = customers.groupby('cluster')[features].agg(['mean', 'std'])
print("\n===== 客户分群结果 =====\n")
 
cluster_profiles = []
for cluster_id in range(optimal_k):
    cluster_data = customers[customers['cluster'] == cluster_id]
    profile = {
        '簇ID': cluster_id,
        '人数': len(cluster_data),
        '占比': f"{len(cluster_data)/len(customers):.1%}",
        '平均年龄': f"{cluster_data['age'].mean():.1f}",
        '平均年收入': f"${cluster_data['annual_income'].mean():,.0f}",
        '消费评分': f"{cluster_data['spending_score'].mean():.1f}",
        '购买频率': f"{cluster_data['purchase_frequency'].mean():.1f}/月",
        '折扣敏感度': f"{cluster_data['discount_sensitivity'].mean():.2f}"
    }
    cluster_profiles.append(profile)
    
    print(f"簇 {cluster_id} ({len(cluster_data)}人, {len(cluster_data)/len(customers):.1%}):")
    print(f"  年龄: {cluster_data['age'].mean():.1f}, 年收入: ${cluster_data['annual_income'].mean():,.0f}")
    print(f"  消费评分: {cluster_data['spending_score'].mean():.1f}, 购买频率: {cluster_data['purchase_frequency'].mean():.1f}/月")
    print(f"  上次购买距今: {cluster_data['days_since_last_purchase'].mean():.1f}天, 折扣敏感度: {cluster_data['discount_sensitivity'].mean():.2f}")
    print()
 
# 可视化聚类结果
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
 
# PCA降维可视化
scatter = axes[0].scatter(X_pca[:, 0], X_pca[:, 1], c=customers['cluster'], 
                          cmap='viridis', alpha=0.6)
axes[0].set_xlabel('PC1')
axes[0].set_ylabel('PC2')
axes[0].set_title('Customer Segments (PCA)')
plt.colorbar(scatter, ax=axes[0], label='Cluster')
 
# 收入 vs 消费评分
for cluster_id in range(optimal_k):
    cluster_data = customers[customers['cluster'] == cluster_id]
    axes[1].scatter(cluster_data['annual_income'], cluster_data['spending_score'], 
                   label=f'Cluster {cluster_id}', alpha=0.6)
axes[1].set_xlabel('Annual Income ($)')
axes[1].set_ylabel('Spending Score')
axes[1].set_title('Income vs Spending Score by Cluster')
axes[1].legend()
 
# 各簇的关键指标雷达图
from math import pi
 
categories = ['Age', 'Income', 'Spending', 'Frequency', 'Recency']
N = len(categories)
 
# 归一化各簇的指标
cluster_means = customers.groupby('cluster')[features[:5]].mean()
cluster_means_norm = (cluster_means - cluster_means.min()) / (cluster_means.max() - cluster_means.min())
 
angles = [n / float(N) * 2 * pi for n in range(N)]
angles += angles[:1]
 
ax = axes[2]
ax = plt.subplot(133, polar=True)
 
colors = plt.cm.viridis(np.linspace(0, 1, optimal_k))
 
for cluster_id in range(optimal_k):
    values = cluster_means_norm.loc[cluster_id].values.tolist()
    values += values[:1]
    ax.plot(angles, values, 'o-', linewidth=2, label=f'Cluster {cluster_id}', color=colors[cluster_id])
    ax.fill(angles, values, alpha=0.1, color=colors[cluster_id])
 
ax.set_xticks(angles[:-1])
ax.set_xticklabels(categories)
ax.set_title('Cluster Profiles (Normalized)')
ax.legend(loc='upper right', bbox_to_anchor=(1.3, 1.0))
 
plt.tight_layout()
plt.savefig('customer_segmentation_results.png', dpi=150)
 
# 保存结果
customers.to_csv('customer_segmentation_results.csv', index=False)
print("\n聚类结果已保存到 customer_segmentation_results.csv")

总结

无监督学习是机器学习中非常有趣的一个领域。没有”正确答案”的约束，反而让算法能够发现数据中隐藏的结构和规律。

聚类算法各有优劣：K-Means简单高效，适合球形簇；层次聚类能发现任意形状的簇，还能展示层级结构；DBSCAN不需要指定簇数，能自动检测噪声；GMM提供概率分配，适合软聚类场景。

降维技术让我们能够可视化和理解高维数据。PCA是线性降维的经典方法，t-SNE和UMAP擅长非线性降维的可视化。

异常检测是另一个重要的无监督任务。Isolation Forest、LOF、One-Class SVM各有适用场景，选择时要考虑数据的特性和检测需求。

实际应用中，很多情况下会结合使用多种方法。比如先用PCA降维加速聚类，再用聚类结果做客户分群。关键是理解每种方法的假设和局限性，选择最适合你数据特点的方案。

---

本文为机器学习实战指南系列文章，主要涵盖无监督学习的核心算法和实践技巧。