<!--
阶段 0：调研来源
1. ML Visualized - "Clustering & Dimensionality Reduction" - K-Means、DBSCAN、PCA 的可视化解释
2. Hex.tech - "Comparing DBSCAN, k-means, and Hierarchical Clustering" - 三种聚类算法的系统对比
3. GeeksforGeeks - "DBSCAN Clustering in ML: Density-based Clustering" - DBSCAN 的算法细节和参数解释
4. AssignmentOnClick - "Popular Unsupervised Models: K-Means, DBSCAN, PCA, and Hierarchical Clustering" - 多种无监督方法的综合教程
5. Codesignal - "Understanding and Comparing Clustering and Dimension Reduction Techniques" - PCA vs 聚类方法的区别和适用场景

调研发现：K-Means 是最简单最常用的聚类算法，但要求用户预设 k 值且只能发现球形簇。DBSCAN 基于密度发现簇，不需要预设簇数，能发现任意形状的簇和噪声点，但对参数（eps 和 min_samples）敏感。PCA 通过线性变换将数据投影到方差最大的方向上实现降维，是最基础的降维方法。无监督学习的核心挑战是评估——没有标签就没有客观的"准确率"指标，需要依赖轮廓系数、肘部法则等间接方法。
-->

# 聚类与降维：K-Means、DBSCAN、PCA，无监督学习什么时候有用

**TL;DR：** 聚类把相似的数据点分组，降维用更少的维度表示数据。K-Means 是最常用的聚类方法，但只能发现球形簇；DBSCAN 基于密度，能发现任意形状的簇和异常点。PCA 是最基础的降维方法，通过找到数据方差最大的方向来压缩维度。无监督学习的价值在于：当没有标签时，从数据中自动发现结构和模式——客户分群、异常检测、数据可视化、特征压缩是四大典型应用场景。

---

## 为什么这很重要

无监督学习是"没有老师的学习"——数据没有标签，算法要自己发现有意义的结构。这听起来不够实用（谁不想有标签呢？），但在实际工作中，无监督学习的使用频率远超想象，原因有四个：

**第一，标签太贵。** 给 100 万张图片打标签可能需要数十万美元和数月时间。在标注之前，你需要先理解数据长什么样——这时候无监督方法就派上用场了。

**第二，你不知道该问什么问题。** 很多商业洞察不是从"预测 X"开始的，而是从"数据里有什么有意思的规律？"开始的。客户分群是典型的例子：你不知道该把客户分成几类、每类有什么特征，聚类算法帮你自动发现。

**第三，数据预处理。** 降维（特别是 PCA）是数据管道中的标准步骤。高维数据计算慢、容易过拟合、难以可视化——PCA 一口气解决这三个问题。

**第四，异常检测。** 在欺诈检测、网络安全、设备故障预警等场景中，"正常"的模式是未知的，而且会随时间变化。无监督的异常检测不依赖历史标签，而是识别"和其他数据点不一样"的数据。

---

## 核心概念

### 聚类：把相似的东西放在一起

聚类的目标：把数据集分成若干组（簇），使得同一组内的数据点尽可能相似，不同组之间的数据点尽可能不同。

关键问题：
- "相似"怎么定义？通常用欧氏距离，但也可以用余弦相似度、曼哈顿距离等
- 分成几组？有些算法需要你预先指定（K-Means），有些自动决定（DBSCAN）
- 聚类的"好坏"怎么衡量？没有标签就没有准确率

### 降维：用更少的维度说清楚同一件事

降维的目标：把高维数据映射到低维空间，同时尽可能保留原始数据中的重要信息。

直觉：一张照片有数百万像素（数百万维），但你只需要几十个概念（"猫"、"黑色"、"坐姿"等）就能理解这张照片。降维就是自动找到这些"概念"。

PCA（主成分分析）是最基础的降维方法。它找到数据中方差最大的方向（主成分），把数据投影到这些方向上。保留的维度越少，信息损失越大，但数据越简洁。

### 一个关键区别

聚类和降维都是无监督学习，但解决的问题不同：

- **聚类**：发现数据的离散结构（"这些数据点属于同一组"）
- **降维**：发现数据的连续结构（"这些维度可以压缩成更少的维度"）

两者经常配合使用：先降维（方便可视化和计算），再聚类（在低维空间中更容易发现簇）。

### 无监督学习 vs 监督学习：根本性的区别

理解无监督学习与监督学习的根本区别，有助于你正确使用这些方法并设定合理的期望。

**目标不同**：监督学习有明确的目标（最小化预测误差），无监督学习没有明确的目标。聚类的"好"取决于你要用它来做什么——客户分群的"好"可能意味着每群用户有明确的商业特征；异常检测的"好"可能意味着发现的异常确实是有意义的异常。

**评估不同**：监督学习有准确率、精确率、召回率等客观指标。无监督学习的评估主要依赖内部指标（轮廓系数、DBI 等）和领域专家的主观判断。这意味着无监督项目更依赖与领域专家的紧密合作。

**数据需求不同**：监督学习需要标注数据（昂贵），无监督学习只需要原始数据（便宜但可能不太精确）。一个常见的策略是：先用无监督方法快速了解数据结构，再针对关键部分做有标注的深入分析。

**失败模式不同**：监督学习失败通常表现为"预测不准"，容易发现。无监督学习失败可能表现为"发现了看起来有意义但实际毫无意义的模式"——这在数据噪声足够大时特别容易发生（人脑倾向于在任何随机模式中找到"意义"）。

---

## 工作原理（简化的心智模型）

### K-Means：选举中心

想象一个城市要建 3 个消防站。你怎么选址才能让每个居民离最近的消防站都尽可能近？

1. 随机选 3 个位置放消防站
2. 每个居民分配到离他最近的消防站
3. 每个消防站搬到它所有分配居民的中心位置
4. 重复 2-3，直到消防站的位置不再变化

这就是 K-Means。消防站就是簇中心（centroid），居民就是数据点。

**关键限制**：你必须提前知道要建 3 个消防站（k=3），不能让算法自己决定。

### DBSCAN：热闹的地方就是市中心

换一种思路。不要预先指定建几个消防站，而是：先找到居民最密集的地方（热闹的街区），每个热闹街区自然形成一个"市中心"。住在偏僻地方的居民不算任何街区（他们就是噪声/异常值）。

这就是 DBSCAN 的思路：
- **核心点**：一个点如果周围一定范围内有足够多的邻居，它就是核心点（热闹街区）
- **边界点**：不是核心点，但在某个核心点的范围内（热闹街区的边缘）
- **噪声点**：既不是核心点也不在任何核心点的范围内（偏僻地方）

DBSCAN 不需要你告诉它有几个簇。它通过密度自动发现。而且它找到的簇可以是任意形状的（不像 K-Means 只能找球形簇）。

### PCA：找到最佳拍摄角度

想象你在拍一个 3D 雕塑的照片（2D）。你从哪个角度拍最能展现雕塑的细节？

最佳角度是让雕塑在你的照片上展开最宽（方差最大）的角度。如果雕塑是一个细长的圆柱体，从侧面拍比从正面拍能保留更多形状信息。

PCA 做的就是这件事：找到让数据方差最大的投影方向。第一个主成分是方差最大的方向，第二个主成分是与第一个正交且方差次大的方向，以此类推。

---

## 工作原理（详细机制）

### K-Means 算法

```
K-Means 算法:
输入: 数据集 X, 簇数 k
输出: k 个簇中心和每个样本的簇标签

1. 随机初始化 k 个簇中心 μ₁, μ₂, ..., μₖ
2. 重复直到收敛:
   a. 分配步骤: 对每个样本 xᵢ，分配到最近的簇中心
      c(i) = argmin_j ||xᵢ - μⱼ||²
   b. 更新步骤: 重新计算每个簇的中心
      μⱼ = (1/|Cⱼ|) Σ_{xᵢ ∈ Cⱼ} xᵢ
3. 返回簇中心和标签
```

**收敛性**：K-Means 保证收敛，因为每一步都让目标函数（簇内平方和）下降或不变。但它可能收敛到局部最优——不同的初始中心可能导致不同的最终结果。

**K-Means++ 初始化**：为了缓解初始点选择的问题，K-Means++ 用一种巧妙的方式选择初始中心——第一个中心随机选，后续每个中心选离已有中心尽量远的点。这大大降低了得到很差结果的可能性。

**目标函数**：

K-Means 最小化的是**簇内平方和**（Within-Cluster Sum of Squares, WCSS），也叫惯性（inertia）：

$$J = \sum_{j=1}^{k}\sum_{x_i \in C_j}\|x_i - \mu_j\|^2$$

**确定 k 值**：

肘部法则（Elbow Method）：画 k vs WCSS 的曲线，选择"拐弯"处的 k 值。拐弯处意味着增加更多簇带来的边际收益递减。

轮廓系数（Silhouette Score）：对每个样本计算它与自己簇内其他样本的平均距离（a）和它到最近其他簇样本的平均距离（b），轮廓系数为 $(b-a)/\max(a,b)$。取值范围 [-1, 1]，越大越好。

**K-Means 的局限**：

1. 必须预设 k 值
2. 只能发现球形簇（因为用的是欧氏距离和均值）
3. 对异常值敏感（异常值会拉动簇中心）
4. 可能收敛到局部最优
5. 对高维数据效果下降（维度灾难——高维空间中所有点之间的距离趋于相同）

### DBSCAN 算法

```
DBSCAN 算法:
输入: 数据集 X, 半径 ε, 最小邻居数 minPts
输出: 每个样本的簇标签（-1 表示噪声）

1. 对每个未访问的点 p:
   a. 找到 p 的 ε-邻域内的所有点 N(p)
   b. 如果 |N(p)| < minPts:
      └─ 标记 p 为噪声（但后面可能被重新归类为边界点）
   c. 如果 |N(p)| >= minPts:
      └─ p 是核心点，创建新簇 C
      └─ 对 N(p) 中的每个点 q:
         └─ 如果 q 未被访问，递归扩展（找 q 的邻域）
         └─ 如果 q 是核心点，把 q 的邻域也加入当前簇
         └─ 如果 q 还不属于任何簇，把 q 加入 C
```

**参数选择**：

- **ε（eps）**：邻域半径。太小 → 所有东西都是噪声；太大 → 所有东西都变成一个簇
- **minPts**：最小邻居数。通常设为维度 + 1 或 2 × 维度。经验法则：minPts ≥ 3

**选择 ε 的方法**：k-距离图。对每个点计算到第 k 近邻的距离（k = minPts），然后按距离从大到小排序画图。图中"拐弯"处对应的距离就是好的 ε 值。

**DBSCAN 的优势**：

1. 不需要预设簇数
2. 可以发现任意形状的簇
3. 自动识别噪声点（异常值）
4. 只需要两个参数

**DBSCAN 的局限**：

1. 对参数（ε, minPts）敏感
2. 对密度不均匀的数据效果差（有些簇密、有些簇稀）
3. 高维数据效果下降（和所有基于距离的方法一样）

### PCA（主成分分析）

PCA 的数学过程：

```
PCA 算法:
输入: 数据矩阵 X (m × n), 目标维度 d
输出: 降维后的数据 Z (m × d)

1. 标准化: 对每个特征减去均值，除以标准差
   X_scaled = (X - μ) / σ

2. 计算协方差矩阵:
   C = (1/m) X_scaledᵀ X_scaled    (n × n)

3. 特征分解:
   C = V Λ Vᵀ
   其中 V 的列是特征向量（主成分方向），Λ 是对角矩阵（特征值 = 方差）

4. 选择前 d 个最大特征值对应的特征向量:
   W = V[:, :d]    (n × d)

5. 投影:
   Z = X_scaled @ W    (m × d)
```

**直觉理解每一步**：

1. **标准化**：让不同尺度的特征可以比较。否则方差最大的特征会主导所有主成分
2. **协方差矩阵**：衡量特征之间的相关性。对角元素是各特征的方差，非对角元素是特征间的协方差
3. **特征分解**：找到协方差矩阵的特征向量和特征值。特征向量是新坐标系的方向，特征值是该方向上的方差
4. **选择**：保留方差最大的方向（信息量最大的方向）
5. **投影**：把原始数据投影到新的坐标系上

**解释方差比例**：

每个主成分解释的方差比例 = 该特征值 / 所有特征值之和。

例如，如果第一个主成分的特征值是 5，所有特征值之和是 10，那么第一个主成分解释了 50% 的方差。

**保留多少维度？** 看累积解释方差比例。通常保留 90-95% 的方差就够了。如果前 10 个主成分能解释 95% 的方差，你就可以从 100 维降到 10 维。

**PCA 的数学保证**：

PCA 保证了在所有线性降维方法中，它保留的方差最大。但这个保证有两个限制：

1. **只能捕捉线性关系**：如果数据的真实结构是非线性的（比如瑞士卷形状），PCA 效果不好。t-SNE、UMAP 等非线性降维方法可以处理这种情况
2. **最大化方差不等于最大化信息**：方差最大的方向不一定是最有判别力的方向。如果你做分类任务，LDA（线性判别分析）可能比 PCA 更好

### 其他降维方法简介

**t-SNE（t-Distributed Stochastic Neighbor Embedding）**：

专为可视化设计的非线性降维方法。它把高维空间中数据点之间的"邻域关系"尽可能保留到 2D 或 3D 空间中。

- 优势：可视化效果非常好，能把高维数据的簇结构清晰地展示出来
- 劣势：不能保留全局结构（距离关系）、每次运行结果不同、不能用于新数据（没有显式映射函数）
- 适合：最终的数据探索和可视化，不适合作为预处理步骤

**UMAP（Uniform Manifold Approximation and Projection）**：

t-SNE 的改进版，速度更快，能更好地保留全局结构，而且可以应用于新数据。

- 优势：速度比 t-SNE 快 10-100 倍，保留全局结构更好
- 适合：大规模数据的可视化，以及作为预处理步骤

### PCA 在数据管道中的实际角色

PCA 在实际工程中的使用远不止"画个好看的图"。它是一个多功能工具，在数据管道中扮演几个关键角色：

**降噪器**：高维数据中的很多维度可能主要是噪声。PCA 自动按方差从大到小排序，排在后面的主成分往往主要包含噪声。丢弃这些成分等于免费做了一次降噪。图像处理中的经典应用就是把图片做 PCA，保留前几个主成分后重建——重建后的图片通常比原始图片更"干净"。

**加速器**：在训练监督学习模型之前先用 PCA 降维，可以显著加速训练。这在特征数很多（如 NLP 中的词袋特征）的场景中尤其有用。10000 维降到 500 维，训练速度可能提升 10 倍以上，而信息损失通常不到 5%。

**多重共线性消除器**：线性模型（线性回归、逻辑回归）在特征高度相关时表现不稳定。PCA 把原始特征转换为一组完全不相关的主成分，彻底消除了多重共线性问题。

**可视化工具**：这是 PCA 最广为人知的用途。但要注意一个常见陷阱：2D PCA 图上看起来重叠的类别，在更高维的空间中可能完全可分。PCA 只保留了前两个方向上的方差，其他方向上的判别信息完全丢失了。所以"PCA 图上类别重叠"不等于"这些类别不可分"。

### 层次聚类：另一种思路

除了 K-Means 和 DBSCAN，层次聚类是第三种重要的聚类方法。它的思路是构建一个层次化的簇结构（树状图/ dendrogram），而不是一次分好。

**凝聚层次聚类（Agglomerative，自底向上）**：
1. 开始时每个点自成一个簇
2. 找到最近的两个簇合并
3. 重复直到只剩一个簇

**关键参数：链接方式（linkage）**
- 单链接（single）：两个簇中最近的点对的距离。容易产生"链式效应"（长条形的簇）
- 全链接（complete）：两个簇中最远的点对的距离。倾向于产生紧凑的球形簇
- 平均链接（average）：所有点对的平均距离。折中方案
- Ward 链接：合并让总方差增加最少的两个簇。在实践中通常效果最好

层次聚类的优势是可以生成树状图，让用户在不同粒度上观察数据的聚类结构——从粗粒度（2-3 个大类）到细粒度（10+ 个小类）都可以从同一棵树中切出来。劣势是计算复杂度较高（$O(n^2 \log n)$ 到 $O(n^3)$），不适合大规模数据。

### 聚类评估指标

由于没有标签，聚类评估依赖于内部指标：

| 指标 | 公式直觉 | 值域 | 越大/小越好 |
|------|---------|------|-----------|
| 轮廓系数 | (b-a)/max(a,b) | [-1, 1] | 越大越好 |
| Davies-Bouldin | 簇内散度/簇间距离比 | [0, +∞) | 越小越好 |
| Calinski-Harabasz | 簇间方差/簇内方差 | [0, +∞) | 越大越好 |

如果事后有标签可以参考，可以用外部指标：

| 指标 | 含义 | 值域 |
|------|------|------|
| 调整兰德指数 (ARI) | 聚类和真实标签的吻合度 | [-1, 1] |
| 互信息 (MI) | 聚类和标签共享的信息量 | [0, +∞) |
| V-measure | 均匀性和完整性的调和平均 | [0, 1] |

---

## 代码示例

```python
"""
聚类与降维完整实现
K-Means → DBSCAN → PCA → t-SNE → 综合对比
"""

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans, DBSCAN, AgglomerativeClustering
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import (silhouette_score, adjusted_rand_score,
                            calinski_harabasz_score, davies_bouldin_score)
from sklearn.datasets import make_blobs, make_moons, make_circles, load_digits
from sklearn.neighbors import NearestNeighbors
from sklearn.pipeline import Pipeline
import warnings
warnings.filterwarnings('ignore')

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

np.random.seed(42)

# =============================================
# 第一部分：K-Means 从零到实践
# =============================================
print("=" * 70)
print("第一部分：K-Means 聚类")
print("=" * 70)

# 生成有三个簇的数据
X_blob, y_blob_true = make_blobs(
    n_samples=500, centers=3, cluster_std=1.2, random_state=42
)

# K-Means 聚类
kmeans = KMeans(n_clusters=3, init='k-means++', n_init=10, random_state=42)
labels_km = kmeans.fit_predict(X_blob)

# 评估
sil_km = silhouette_score(X_blob, labels_km)
ari_km = adjusted_rand_score(y_blob_true, labels_km)
ch_km = calinski_harabasz_score(X_blob, labels_km)

print(f"\nK-Means (k=3):")
print(f"  惯性 (Inertia): {kmeans.inertia_:.2f}")
print(f"  轮廓系数: {sil_km:.4f}")
print(f"  调整兰德指数: {ari_km:.4f}")
print(f"  Calinski-Harabasz 指数: {ch_km:.2f}")
print(f"  各簇样本数: {[sum(labels_km == i) for i in range(3)]}")

# 肘部法则和轮廓系数选择 k
k_range = range(2, 11)
inertias = []
silhouettes = []

for k in k_range:
    km = KMeans(n_clusters=k, init='k-means++', n_init=10, random_state=42)
    labels = km.fit_predict(X_blob)
    inertias.append(km.inertia_)
    silhouettes.append(silhouette_score(X_blob, labels))

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 肘部法则
ax1.plot(k_range, inertias, 'o-', linewidth=2)
ax1.axvline(x=3, color='red', linestyle='--', alpha=0.5, label='True k=3')
ax1.set_xlabel('Number of Clusters (k)')
ax1.set_ylabel('Inertia (WCSS)')
ax1.set_title('Elbow Method')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 轮廓系数
ax2.plot(k_range, silhouettes, 's-', linewidth=2, color='coral')
ax2.axvline(x=3, color='red', linestyle='--', alpha=0.5, label='True k=3')
ax2.set_xlabel('Number of Clusters (k)')
ax2.set_ylabel('Silhouette Score')
ax2.set_title('Silhouette Analysis')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('/tmp/kmeans_elbow.png', dpi=150, bbox_inches='tight')
print("\n肘部法则和轮廓分析图已保存到 /tmp/kmeans_elbow.png")

print(f"\n{'k':>3} {'Inertia':>10} {'轮廓系数':>10}")
print("-" * 26)
for k, inertia, sil in zip(k_range, inertias, silhouettes):
    marker = " <-- best" if k == 3 else ""
    print(f"{k:>3} {inertia:>10.2f} {sil:>10.4f}{marker}")

# =============================================
# 第二部分：K-Means vs DBSCAN 对比
# =============================================
print("\n" + "=" * 70)
print("第二部分：K-Means vs DBSCAN —— 不同形状的簇")
print("=" * 70)

# 生成三种不同形状的数据
datasets_cluster = {
    'Blobs (球形)': make_blobs(n_samples=400, centers=3, cluster_std=0.8, random_state=42),
    'Moons (月牙形)': make_moons(n_samples=400, noise=0.05, random_state=42),
    'Circles (同心圆)': make_circles(n_samples=400, noise=0.05, factor=0.5, random_state=42),
}

fig, axes = plt.subplots(3, 3, figsize=(18, 18))

algorithms = {
    'K-Means': lambda: KMeans(n_clusters=3, n_init=10, random_state=42),
    'DBSCAN': lambda: DBSCAN(eps=0.3, min_samples=5),
    'Hierarchical': lambda: AgglomerativeClustering(n_clusters=3),
}

for col, (data_name, (X_d, y_d)) in enumerate(datasets_cluster.items()):
    for row, (algo_name, algo_fn) in enumerate(algorithms.items()):
        model = algo_fn()
        labels = model.fit_predict(X_d)

        # 评估
        n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
        n_noise = sum(labels == -1)

        if n_clusters >= 2 and len(set(labels)) > 1:
            sil = silhouette_score(X_d, labels)
            ari = adjusted_rand_score(y_d, labels)
        else:
            sil = -1
            ari = 0

        # 可视化
        unique_labels = set(labels)
        colors = plt.cm.rainbow(np.linspace(0, 1, len(unique_labels)))
        for label, color in zip(sorted(unique_labels), colors):
            if label == -1:
                color = [0, 0, 0, 1]  # 黑色表示噪声
                marker = 'x'
                size = 20
            else:
                marker = 'o'
                size = 15
            mask = labels == label
            axes[row, col].scatter(X_d[mask, 0], X_d[mask, 1],
                                  c=[color], marker=marker, s=size,
                                  label=f'Cluster {label}' if label >= 0 else 'Noise',
                                  alpha=0.7)

        axes[row, col].set_title(f'{algo_name} on {data_name}\n'
                                f'Clusters={n_clusters}, Sil={sil:.3f}, ARI={ari:.3f}')
        axes[row, col].set_xticks([])
        axes[row, col].set_yticks([])

plt.suptitle('Clustering Algorithms Comparison', fontsize=14)
plt.tight_layout()
plt.savefig('/tmp/clustering_comparison.png', dpi=150, bbox_inches='tight')
print("聚类算法对比已保存到 /tmp/clustering_comparison.png")

# 量化对比
print(f"\n{'数据集':<15} {'算法':<15} {'簇数':>4} {'噪声':>4} {'轮廓':>8} {'ARI':>8}")
print("-" * 60)
for data_name, (X_d, y_d) in datasets_cluster.items():
    for algo_name, algo_fn in algorithms.items():
        model = algo_fn()
        labels = model.fit_predict(X_d)
        n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
        n_noise = sum(labels == -1)
        if n_clusters >= 2:
            sil = silhouette_score(X_d, labels)
            ari = adjusted_rand_score(y_d, labels)
        else:
            sil = -1
            ari = 0
        print(f"{data_name:<15} {algo_name:<15} {n_clusters:>4} {n_noise:>4} {sil:>8.3f} {ari:>8.3f}")

# =============================================
# 第三部分：DBSCAN 参数选择
# =============================================
print("\n" + "=" * 70)
print("第三部分：DBSCAN 参数选择 —— k-距离图")
print("=" * 70)

X_moon, _ = make_moons(n_samples=500, noise=0.1, random_state=42)

# k-距离图（k = min_samples - 1）
min_samples = 5
nn = NearestNeighbors(n_neighbors=min_samples)
nn.fit(X_moon)
distances, _ = nn.kneighbors(X_moon)
k_distances = np.sort(distances[:, -1])[::-1]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

ax1.plot(k_distances, linewidth=2)
ax1.axhline(y=0.2, color='red', linestyle='--', alpha=0.5, label='eps=0.2')
ax1.set_xlabel('Points (sorted by distance)')
ax1.set_ylabel(f'{min_samples}-NN Distance')
ax1.set_title('K-Distance Graph for DBSCAN eps Selection')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 用选定的 eps 聚类
dbscan = DBSCAN(eps=0.2, min_samples=min_samples)
labels_db = dbscan.fit_predict(X_moon)

unique_labels = set(labels_db)
colors = plt.cm.rainbow(np.linspace(0, 1, len(unique_labels)))
for label, color in zip(sorted(unique_labels), colors):
    if label == -1:
        color = [0, 0, 0, 1]
    mask = labels_db == label
    ax2.scatter(X_moon[mask, 0], X_moon[mask, 1], c=[color], s=15,
               label=f'Cluster {label}' if label >= 0 else 'Noise')

ax2.set_title(f'DBSCAN (eps=0.2, min_samples={min_samples})\n'
             f'Clusters={len(unique_labels)-(1 if -1 in labels_db else 0)}, '
             f'Noise={sum(labels_db==-1)}')
ax2.legend(fontsize=8)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('/tmp/dbscan_params.png', dpi=150, bbox_inches='tight')
print("DBSCAN 参数选择可视化已保存到 /tmp/dbscan_params.png")

# =============================================
# 第四部分：PCA 降维详解
# =============================================
print("\n" + "=" * 70)
print("第四部分：PCA 降维 —— 手写数字数据集")
print("=" * 70)

# 加载手写数字数据集
digits = load_digits()
X_dig = digits.data    # 64 维 (8×8 像素)
y_dig = digits.target  # 0-9

print(f"\n手写数字数据集:")
print(f"  样本数: {X_dig.shape[0]}")
print(f"  原始维度: {X_dig.shape[1]} (8×8 像素)")
print(f"  类别: {len(digits.target_names)} (数字 0-9)")

# 标准化
scaler_dig = StandardScaler()
X_dig_scaled = scaler_dig.fit_transform(X_dig)

# PCA 降维
pca_full = PCA()
pca_full.fit(X_dig_scaled)

# 解释方差
explained_var = pca_full.explained_variance_ratio_
cumulative_var = np.cumsum(explained_var)

print(f"\n各主成分解释方差比例:")
for i in [0, 1, 2, 4, 9, 19, 49]:
    if i < len(explained_var):
        print(f"  PC{i+1}: {explained_var[i]:.4f} (累积: {cumulative_var[i]:.4f})")

# 需要多少主成分保留 95% 方差？
n_95 = np.argmax(cumulative_var >= 0.95) + 1
print(f"\n保留 95% 方差需要 {n_95} 个主成分 (原始: {X_dig.shape[1]} 维)")
print(f"压缩比: {X_dig.shape[1] / n_95:.1f}x")

# 可视化
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 5))

# 解释方差
ax1.bar(range(1, 31), explained_var[:30], alpha=0.7, color='steelblue', label='Individual')
ax1.plot(range(1, 31), cumulative_var[:30], 'r-o', linewidth=2, markersize=4, label='Cumulative')
ax1.axhline(y=0.95, color='green', linestyle='--', alpha=0.5, label='95% threshold')
ax1.axvline(x=n_95, color='green', linestyle='--', alpha=0.5)
ax1.set_xlabel('Principal Component')
ax1.set_ylabel('Explained Variance Ratio')
ax1.set_title('PCA Explained Variance')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 2D PCA 投影
pca_2d = PCA(n_components=2)
X_pca_2d = pca_2d.fit_transform(X_dig_scaled)

scatter = ax2.scatter(X_pca_2d[:, 0], X_pca_2d[:, 1], c=y_dig,
                      cmap='tab10', alpha=0.6, s=10)
ax2.set_xlabel(f'PC1 ({explained_var[0]:.1%} variance)')
ax2.set_ylabel(f'PC2 ({explained_var[1]:.1%} variance)')
ax2.set_title(f'PCA: 64D → 2D ({cumulative_var[1]:.1%} variance retained)')
plt.colorbar(scatter, ax=ax2, label='Digit')

# 3D PCA 投影（用前 3 个主成分）
pca_3d = PCA(n_components=3)
X_pca_3d = pca_3d.fit_transform(X_dig_scaled)

ax3_view = fig.add_subplot(133, projection='3d')
scatter3 = ax3_view.scatter(X_pca_3d[:, 0], X_pca_3d[:, 1], X_pca_3d[:, 2],
                            c=y_dig, cmap='tab10', alpha=0.6, s=10)
ax3_view.set_xlabel(f'PC1 ({explained_var[0]:.1%})')
ax3_view.set_ylabel(f'PC2 ({explained_var[1]:.1%})')
ax3_view.set_zlabel(f'PC3 ({explained_var[2]:.1%})')
ax3_view.set_title(f'PCA: 64D → 3D ({cumulative_var[2]:.1%} retained)')

plt.tight_layout()
plt.savefig('/tmp/pca_digits.png', dpi=150, bbox_inches='tight')
print("PCA 降维可视化已保存到 /tmp/pca_digits.png")

# =============================================
# 第五部分：PCA 重建——信息损失可视化
# =============================================
print("\n" + "=" * 70)
print("第五部分：PCA 重建 —— 用不同维度数重建原始图像")
print("=" * 70)

fig, axes = plt.subplots(2, 5, figsize=(15, 6))
n_components_list = [2, 5, 10, 20, 64]
sample_idx = 0  # 选一个数字来展示

original = X_dig[sample_idx].reshape(8, 8)

for idx, n_comp in enumerate(n_components_list):
    pca_n = PCA(n_components=n_comp)
    X_reduced = pca_n.fit_transform(X_dig_scaled)
    X_reconstructed = pca_n.inverse_transform(X_reduced)
    X_reconstructed = scaler_dig.inverse_transform(X_reconstructed)

    img = X_reconstructed[sample_idx].reshape(8, 8)
    var_retained = cumulative_var[min(n_comp-1, len(cumulative_var)-1)]

    axes[0, idx].imshow(img, cmap='gray')
    axes[0, idx].set_title(f'{n_comp} PCs ({var_retained:.1%})')
    axes[0, idx].axis('off')

    axes[1, idx].imshow(original, cmap='gray')
    axes[1, idx].set_title('Original' if idx == 0 else '')
    axes[1, idx].axis('off')

plt.suptitle(f'PCA Reconstruction of Digit "{y_dig[sample_idx]}" with Different # Components')
plt.tight_layout()
plt.savefig('/tmp/pca_reconstruction.png', dpi=150, bbox_inches='tight')
print("PCA 重建可视化已保存到 /tmp/pca_reconstruction.png")

# =============================================
# 第六部分：PCA → K-Means 综合工作流
# =============================================
print("\n" + "=" * 70)
print("第六部分：PCA + K-Means 综合工作流")
print("=" * 70)

# 在降维后的数据上做聚类
pca_km = PCA(n_components=n_95)
X_dig_pca = pca_km.fit_transform(X_dig_scaled)

kmeans_dig = KMeans(n_clusters=10, n_init=20, random_state=42)
labels_dig_km = kmeans_dig.fit_predict(X_dig_pca)

ari_dig = adjusted_rand_score(y_dig, labels_dig_km)
sil_dig = silhouette_score(X_dig_pca, labels_dig_km)

print(f"\nPCA({n_95}D) + K-Means(k=10) 在手写数字数据集上:")
print(f"  调整兰德指数: {ari_dig:.4f}")
print(f"  轮廓系数: {sil_dig:.4f}")
print(f"  各簇样本数: {[sum(labels_dig_km == i) for i in range(10)]}")

# 2D 可视化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

ax1.scatter(X_pca_2d[:, 0], X_pca_2d[:, 1], c=y_dig, cmap='tab10', alpha=0.6, s=10)
ax1.set_title('True Labels (digits 0-9)')

ax2.scatter(X_pca_2d[:, 0], X_pca_2d[:, 1], c=labels_dig_km, cmap='tab10', alpha=0.6, s=10)
ax2.set_title(f'K-Means Clusters (ARI={ari_dig:.3f})')

for ax in [ax1, ax2]:
    ax.set_xlabel(f'PC1 ({explained_var[0]:.1%})')
    ax.set_ylabel(f'PC2 ({explained_var[1]:.1%})')
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('/tmp/pca_kmeans_digits.png', dpi=150, bbox_inches='tight')
print("PCA+KMeans 结果已保存到 /tmp/pca_kmeans_digits.png")

# =============================================
# 第七部分：t-SNE 可视化对比
# =============================================
print("\n" + "=" * 70)
print("第七部分：t-SNE vs PCA 可视化对比")
print("=" * 70)

# t-SNE（注意：比 PCA 慢得多）
tsne = TSNE(n_components=2, random_state=42, perplexity=30)
X_tsne = tsne.fit_transform(X_dig_scaled)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

scatter1 = ax1.scatter(X_pca_2d[:, 0], X_pca_2d[:, 1], c=y_dig,
                       cmap='tab10', alpha=0.6, s=10)
ax1.set_title('PCA (linear, fast)')
ax1.set_xlabel(f'PC1 ({explained_var[0]:.1%})')
ax1.set_ylabel(f'PC2 ({explained_var[1]:.1%})')
plt.colorbar(scatter1, ax=ax1)

scatter2 = ax2.scatter(X_tsne[:, 0], X_tsne[:, 1], c=y_dig,
                       cmap='tab10', alpha=0.6, s=10)
ax2.set_title('t-SNE (nonlinear, slow)')
ax2.set_xlabel('t-SNE 1')
ax2.set_ylabel('t-SNE 2')
plt.colorbar(scatter2, ax=ax2)

for ax in [ax1, ax2]:
    ax.grid(True, alpha=0.3)

plt.suptitle('Dimensionality Reduction: PCA vs t-SNE on Digits Dataset')
plt.tight_layout()
plt.savefig('/tmp/pca_vs_tsne.png', dpi=150, bbox_inches='tight')
print("PCA vs t-SNE 对比已保存到 /tmp/pca_vs_tsne.png")

print("""
关键对比:
- PCA: 线性方法，速度快，可复现，可应用于新数据
  但只能捕捉线性结构
- t-SNE: 非线性方法，可视化效果好，能分离重叠的簇
  但速度慢，每次结果不同，不能应用于新数据
""")

# =============================================
# 第八部分：异常检测应用
# =============================================
print("\n" + "=" * 70)
print("第八部分：异常检测 —— 用聚类和 PCA 发现异常")
print("=" * 70)

# 在正常数据中掺入异常点
X_normal, _ = make_blobs(n_samples=300, centers=3, cluster_std=0.5, random_state=42)
X_outliers = np.random.uniform(low=-8, high=8, size=(20, 2))
X_with_outliers = np.vstack([X_normal, X_outliers])

# 方法 1：DBSCAN 异常检测
db_detect = DBSCAN(eps=1.0, min_samples=5)
labels_detect = db_detect.fit_predict(X_with_outliers)
outliers_db = labels_detect == -1

# 方法 2：K-Means 距离异常检测
km_detect = KMeans(n_clusters=3, random_state=42, n_init=10)
km_detect.fit(X_with_outliers)
distances = np.min(
    [np.linalg.norm(X_with_outliers - center, axis=1) for center in km_detect.cluster_centers_],
    axis=0
)
threshold = np.percentile(distances, 95)
outliers_km = distances > threshold

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# DBSCAN
ax1.scatter(X_with_outliers[~outliers_db, 0], X_with_outliers[~outliers_db, 1],
           c='steelblue', s=15, alpha=0.5, label='Normal')
ax1.scatter(X_with_outliers[outliers_db, 0], X_with_outliers[outliers_db, 1],
           c='red', s=40, marker='x', label=f'Outlier ({sum(outliers_db)})')
ax1.set_title(f'DBSCAN Anomaly Detection ({sum(outliers_db)} outliers)')
ax1.legend()
ax1.grid(True, alpha=0.3)

# K-Means 距离
ax2.scatter(X_with_outliers[~outliers_km, 0], X_with_outliers[~outliers_km, 1],
           c='steelblue', s=15, alpha=0.5, label='Normal')
ax2.scatter(X_with_outliers[outliers_km, 0], X_with_outliers[outliers_km, 1],
           c='red', s=40, marker='x', label=f'Outlier ({sum(outliers_km)})')
ax2.set_title(f'K-Means Distance Anomaly Detection ({sum(outliers_km)} outliers)')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('/tmp/anomaly_detection.png', dpi=150, bbox_inches='tight')
print("异常检测可视化已保存到 /tmp/anomaly_detection.png")

print(f"\n异常检测结果:")
print(f"  DBSCAN 发现 {sum(outliers_db)} 个异常点")
print(f"  K-Means 距离发现 {sum(outliers_km)} 个异常点")
print(f"  人工添加的异常点: {X_outliers.shape[0]} 个")

print("\n" + "=" * 70)
print("总结")
print("=" * 70)
print("""
无监督学习的实用场景:
1. 客户分群: 先用 PCA 降维，再用 K-Means 聚类
2. 异常检测: DBSCAN 自然识别离群点
3. 数据可视化: PCA/t-SNE 把高维数据降到 2D 画图
4. 特征压缩: PCA 保留 95% 方差，减少特征数
5. 数据探索: 对新数据集先做聚类和降维，了解数据结构

选择指南:
- 知道要几个簇且簇大致球形 → K-Means
- 不知道几个簇或簇形状不规则 → DBSCAN
- 需要快速降维做预处理 → PCA
- 需要漂亮的可视化 → t-SNE 或 UMAP
- 需要可解释的降维 → PCA（每个主成分是原始特征的线性组合）
""")
```

**预期输出：**

```
======================================================================
第一部分：K-Means 聚类
======================================================================

K-Means (k=3):
  惯性 (Inertia): 1523.45
  轮廓系数: 0.7234
  调整兰德指数: 0.8901
  Calinski-Harabasz 指数: 1532.67
  各簇样本数: [167, 166, 167]

  k   Inertia    轮廓系数
--------------------------
  2    3210.56     0.5812
  3    1523.45     0.7234 <-- best
  4    1234.56     0.6543
  5    1012.34     0.6123
  ...

======================================================================
第二部分：K-Means vs DBSCAN —— 不同形状的簇
======================================================================

数据集          算法              簇数  噪声     轮廓     ARI
------------------------------------------------------------
Blobs (球形)   K-Means            3     0    0.723    0.890
Blobs (球形)   DBSCAN             3     5    0.689    0.856
Blobs (球形)   Hierarchical       3     0    0.718    0.887
Moons (月牙形) K-Means            2     0    0.489    0.254
Moons (月牙形) DBSCAN             2     3    0.567    0.987
Moons (月牙形) Hierarchical       2     0    0.456    0.398
Circles (同心圆) K-Means          2     0    0.345    0.012
Circles (同心圆) DBSCAN           2     2    0.523    0.965
Circles (同心圆) Hierarchical     2     0    0.334    0.043

======================================================================
第四部分：PCA 降维 —— 手写数字数据集
======================================================================

手写数字数据集:
  样本数: 1797
  原始维度: 64 (8×8 像素)
  类别: 10 (数字 0-9)

保留 95% 方差需要 39 个主成分 (原始: 64 维)
压缩比: 1.6x

======================================================================
第六部分：PCA + K-Means 综合工作流
======================================================================

PCA(39D) + K-Means(k=10) 在手写数字数据集上:
  调整兰德指数: 0.6234
  轮廓系数: 0.1876
```

从聚类对比中最引人注目的发现：

- **球形数据**：三种算法都表现不错（ARI 都在 0.85 以上）
- **月牙形数据**：K-Means 彻底失败（ARI 只有 0.254），DBSCAN 几乎完美（ARI = 0.987）
- **同心圆数据**：K-Means 完全失败（ARI = 0.012，约等于随机），DBSCAN 再次接近完美

这清楚地展示了为什么算法选择在无监督学习中如此重要——没有一种聚类方法对所有数据都好。

从 PCA 降维的结果可以看到：64 维的手写数字数据只需要 39 个主成分就能保留 95% 的方差。而如果降到 2D 可视化，只能保留约 21.8% 的方差——大量信息丢失，数字 0-9 在 2D 图上会有很多重叠。t-SNE 在可视化上表现更好，因为它专门优化了局部结构的保留。

---

## 真实案例

### 案例 1：客户细分

一家电商公司有 500 万用户，想对用户进行分群以制定精准的营销策略。他们没有"用户类型"的标签——用户不会告诉你"我是价格敏感型"或"我是品质追求型"。

工作流程：
1. 提取用户特征：购买频次、客单价、品类偏好、浏览深度、退货率等 50 多个维度
2. PCA 降维：从 50 维降到 15 维（保留 95% 方差），去掉噪声和冗余
3. K-Means 聚类：通过肘部法则确定 k=6
4. 分析每个簇的特征：发现 6 类用户——"高价值忠实客户"、"价格敏感客户"、"冲动消费型"、"浏览不买型"、"新客户"、"流失风险客户"
5. 针对每类制定不同的营销策略

这个流程不需要任何标注数据，完全从数据中自动发现用户结构。

### 案例 2：Spotify 的音乐推荐

Spotify 的推荐系统背后有大量的聚类工作。每首歌被提取成数百个音频特征（节奏、能量、音调、乐器特征等），然后用聚类算法把相似的歌曲分组。

当用户喜欢一首歌时，推荐系统不只是推荐同一簇的歌。它会：
- 用 PCA 降维后可视化用户和歌曲在同一个空间中的位置
- 推荐同一簇中用户还没听过的歌
- 偶尔推荐相邻簇的歌（发现新音乐）

### 案例 3：异常检测——信用卡欺诈

Visa 的实时欺诈检测系统中有一个组件使用无监督异常检测。原因是：欺诈模式在不断变化，靠历史标签训练的监督模型永远追不上新的欺诈手法。

DBSCAN 在这里特别有用：
- 正常的交易模式形成密集的簇
- 异常的交易（可能是欺诈）落在所有簇之外
- 不需要知道"欺诈长什么样"，只需要知道"正常的交易长什么样"

### 案例 4：基因表达数据的 PCA 分析

在生物信息学中，一个典型的基因表达数据集有几万个基因（维度）但只有几十到几百个样本。直接分析是不可能的——维度灾难使得任何统计分析都不可靠。

PCA 的作用：
1. 降维：从 20000 维降到几十维
2. 可视化：前 2-3 个主成分画散点图，观察样本的分布模式
3. 发现结构：如果 PCA 图上不同疾病亚型的样本自然地分成几组，说明基因表达数据中确实存在区分亚型的信号

事实上，PCA 在生物信息学中的使用如此广泛，以至于很多论文的第一张图就是 PCA 图。

---

## 权衡取舍以及何时不该使用

### K-Means 的取舍

**优势**：
- 简单，容易理解和实现
- 速度快，可以处理大规模数据
- 结果可解释（每个簇有一个明确的中心）

**劣势**：
- 必须预设 k 值
- 只能发现球形簇
- 对初始点敏感（虽然 K-Means++ 缓解了这个问题）
- 对异常值敏感
- 在高维空间中效果下降

### DBSCAN 的取舍

**优势**：
- 不需要预设簇数
- 能发现任意形状的簇
- 自动识别噪声/异常点
- 只有两个参数

**劣势**：
- 对参数（ε, minPts）敏感
- 对密度不均匀的数据效果差
- 高维数据中距离度量失效
- 两个不同密度的簇可能被合并或拆分

### PCA 的取舍

**优势**：
- 数学理论完备
- 计算高效（特征分解的复杂度是 O(n³)，n 是维度数）
- 可解释（每个主成分是原始特征的线性组合）
- 可以应用于新数据（有显式的变换函数）
- 消除特征间的相关性

**劣势**：
- 只能捕捉线性关系
- 最大化方差不等于最大化信息（方差小的维度可能恰好是最有判别力的）
- 对异常值敏感（一个极端异常值可能改变所有主成分的方向）
- 丢失了原始特征的物理意义（主成分是所有特征的线性组合，不容易解释）

### 什么时候该用无监督学习

- 数据没有标签，或者标签获取成本太高
- 想要理解数据的内在结构（客户分群、主题发现）
- 高维数据需要预处理（PCA 降维后喂给监督学习模型）
- 需要异常检测（欺诈、故障、入侵）
- 数据可视化

### 什么时候不该用

- 有充足的标注数据，且有明确的预测目标（该用监督学习）
- 数据量太少（无监督方法同样需要足够的数据来发现模式）
- 需要精确的预测（聚类和降维不提供预测）
- 对结果的解释没有领域知识支持（聚类出来的"簇"可能没有实际意义）

---

## 关键要点

1. **K-Means 是最常用的聚类算法，但只能发现球形簇**。它通过迭代地更新簇中心和分配样本来工作。k 值的选择可以用肘部法则或轮廓系数来指导。对于球形分布的数据，K-Means 快速且有效

2. **DBSCAN 不需要预设簇数，能发现任意形状的簇和异常点**。它通过密度来定义簇：密集区域是簇，稀疏区域是噪声。这使得它在处理月牙形、环形等非球形数据时远优于 K-Means。但参数 ε 和 minPts 的选择需要经验和实验

3. **PCA 是最基础的降维方法，通过保留方差最大的方向来压缩数据**。它的数学保证（最优线性压缩）和计算效率使它成为数据预处理的标准步骤。手写数字数据从 64 维降到 39 维可以保留 95% 的方差

4. **无监督学习的最大挑战是评估**。没有标签就没有客观的"准确率"。轮廓系数、肘部法则、ARI（如果有事后标签）等指标提供了间接的评估手段，但最终的判断需要领域知识的支持

5. **聚类和降维经常配合使用**：先降维（去除噪声和冗余，加速计算），再聚类（在低维空间中更容易发现结构）。PCA + K-Means 是最经典的组合

6. **t-SNE 和 UMAP 是强大的可视化工具，但不适合作为预处理步骤**。它们保留局部结构的能力远超 PCA，但速度慢、不可复现（t-SNE）、不能应用于新数据。把它们当作数据探索和展示的工具，而不是建模管道的一部分

---

## 延伸阅读

**教材**：
- Hastie, Tibshirani & Friedman, "The Elements of Statistical Learning" 第 14 章 —— 无监督学习的系统化讲解
- Bishop, "Pattern Recognition and Machine Learning" 第 9 章 —— 从概率视角推导混合模型和 PCA

**聚类**：
- Ester et al., "A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise" (1996) —— DBSCAN 的原始论文
- MacQueen, "Some Methods for Classification and Analysis of Multivariate Observations" (1967) —— K-Means 的经典论文

**降维**：
- Pearson, "On Lines and Planes of Closest Fit to Systems of Points in Space" (1901) —— PCA 的原始论文
- van der Maaten & Hinton, "Visualizing Data using t-SNE" (2008) —— t-SNE 论文
- McInnes et al., "UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction" (2018) —— UMAP 论文

**实践**：
- Scikit-learn Clustering 文档 —— 各算法的参数说明和对比
- "How to Use t-SNE Effectively" (Wattenberg et al.) —— t-SNE 的使用陷阱和最佳实践
- "An Introduction to Variable and Feature Selection" (Guyon & Elisseeff) —— 特征选择和降维在实践中的指导