Python 在科研中的应用 12：聚类分析

Posted on Edited on In Waline:

今天总结了最常见的五种聚类算法。涉及到的算法有：

  • K均值聚类（K-Means Clustering）
  • 层次聚类（Hierarchical Clustering）
  • DBSCAN（Density-Based Spatial Clustering of Applications with Noise）
  • 高斯混合模型（Gaussian Mixture Model，GMM）
  • 谱聚类（Spectral Clustering）

课程作业03 占总成绩25%

构建四组二维空间中的随机分布散点，每组散点符合二维高斯分布，高斯分布的中心点与两个方向的标准差均随机产生，且二维高斯分布的中心点落在[[0,100],[0,100]]的方形区域内。

通过K均值聚类分析方法对上述随机散点执行聚类分析，绘制分类结果，并计算K-means聚类分析的准确率。

K均值聚类（K-Means Clustering）

K均值聚类（K-Means Clustering）是一种常用的无监督学习算法，用于将数据集分成K个簇（clusters）。其目标是将相似的数据点归为同一簇，而将不同的数据点分到不同的簇中，从而使得每个簇内的数据点之间的相似度最大，而不同簇之间的相似度最小。

核心原理

K均值聚类的核心思想是通过迭代优化的方法，最小化簇内点到簇中心的总距离平方和。

基本步骤如下：

  1. 初始化：随机选择K个点作为初始簇中心。

  2. 分配簇：将每个数据点分配到最近的簇中心。

  3. 更新簇中心：重新计算每个簇的中心（即簇内所有点的平均值）。

  4. 重复：重复步骤2和3，直到簇中心不再发生显著变化或达到预定的迭代次数。

核心公式

假设数据集为$x_1,x_2,…,x_n$，每个数据点$x_i$有$d$维。簇的中心为$\mu_1,\mu_2,…,\mu_K$。

  1. 距离计算：通常使用欧氏距离：
    \begin{equation}
    {\bf dist}(x_i,\mu_k) = \Arrowvert x_i-\mu_k\Arrowvert_2 = \sqrt{\sum^d_{j=1} (x_{ij}-\mu_{kj})^2}
    \end{equation}

  2. 簇分配：将每个点分配到最近的簇中心：
    \begin{equation}
    C_k = {x_i:\Arrowvert x_i-\mu_k \Arrowvert_2 \leq \Arrowvert x_i - \mu_j\Arrowvert_2, \forall j,1 \leq j \leq K }
    \end{equation}

  3. 更新簇中心：重新计算每个簇的中心：
    \begin{equation}
    \mu_k = \frac{1}{|C_k|} \sum_{x_i \in C_k} x_i
    \end{equation}

Python示例代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans

# 生成样本数据
X, y = make_blobs(n_samples=500, centers=4, cluster_std=0.60, random_state=0)

# 使用KMeans进行聚类
kmeans = KMeans(n_clusters=4)
kmeans.fit(X)
y_kmeans = kmeans.predict(X)

# 绘制结果图
plt.scatter(X[:, 0], X[:, 1], c=y_kmeans, s=50, cmap='viridis')

centers = kmeans.cluster_centers_
plt.scatter(centers[:, 0], centers[:, 1], c='red', s=200, alpha=0.75)
plt.title("K-Means Clustering")
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")
plt.show()

  1. 生成数据：使用make_blobs生成了500个二维数据点，分为4个簇。

  2. 训练模型：使用KMeans模型进行训练，设定簇数为4。

  3. 预测与绘图：根据训练结果进行预测，并将数据点按簇绘制成不同颜色，同时用红色标注出簇中心。

大家可以清晰地看到K均值聚类的效果，即每个簇内的数据点被分为同一种颜色，而不同簇的数据点被分为不同颜色，同时簇中心用红色圆点表示。

层次聚类（Hierarchical Clustering）

层次聚类（Hierarchical Clustering）是一种常见的无监督学习算法，旨在通过构建层次树状结构将数据点进行分组。它可以分为两种类型：自底向上的聚合（agglomerative）和自顶向下的分裂（divisive）聚类方法。

核心原理

层次聚类的核心思想是通过递归地合并或分裂簇，构建一个层次树（树状图或树状图），表示数据的嵌套分组。这里我们以自底向上的聚合层次聚类为例：

  1. 初始化：将每个数据点视为一个独立的簇。

  2. 计算距离：计算所有簇之间的距离矩阵。

  3. 合并簇：找到距离最近的两个簇并合并它们。

  4. 更新距离矩阵：更新合并后的簇与其他簇之间的距离。

  5. 重复：重复步骤3和4，直到所有数据点合并成一个簇。

核心公式

假设数据集为$x_1, x_2, …,x_n$，簇之间的距离可以通过多种方式计算，包括最小距离（single linkage）、最大距离（complete linkage）和平均距离（average linkage）等。

  1. 距离计算：假设簇$A$和簇$B$中的点分别为$a_1, a_2,…,a_p$和$b_1,b_2,…b_q$，则不同的距离计算方法如下：

  • 最小距离（单链法）：
    \begin{equation}
    d(A,B) = \min_{i,j} \Arrowvert a_i-b_j \Arrowvert
    \end{equation}

  • 最大距离（全链法）：
    \begin{equation}
    d(A,B) = \max_{i,j} \Arrowvert a_i-b_j \Arrowvert
    \end{equation}

  • 平均距离（平均链法）：
    \begin{equation}
    d(A,B) = \frac{1}{pq} \sum^{p}{i=1} \sum^{q}{j=1} \Arrowvert a_i-b_j \Arrowvert
    \end{equation}

Python案例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import numpy as np
import matplotlib.pyplot as plt
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.datasets import make_blobs

# 生成样本数据
X, y = make_blobs(n_samples=500, centers=3, random_state=0, cluster_std=0.60)

# 进行层次聚类
Z = linkage(X, method='ward')

# 绘制树状图
plt.figure(figsize=(10, 7))
dendrogram(Z, truncate_mode='lastp', p=12, leaf_rotation=90., leaf_font_size=12., show_contracted=True)
plt.title("Hierarchical Clustering Dendrogram")
plt.xlabel("Sample index or (Cluster size)")
plt.ylabel("Distance")
plt.show()

  1. 生成数据：使用make_blobs生成500个二维数据点，分为3个簇。

  2. 层次聚类：使用linkage函数进行层次聚类，采用Ward方法计算簇间距离。

  3. 绘制树状图：使用dendrogram函数绘制树状图，展示聚类的层次结构。

通过树状图，大家可以清晰地看到层次聚类的合并过程，从最底层的单个数据点逐渐合并到最终的一个簇。树状图的高度表示簇之间的距离，合并越早的簇之间距离越近，合并越晚的簇之间距离越远。

DBSCAN（Density-Based Spatial Clustering of Applications with Noise）

DBSCAN（Density-Based Spatial Clustering of Applications with Noise）是一种基于密度的聚类算法，可以发现任意形状的簇，同时能够识别出噪声点。与K均值和层次聚类不同，DBSCAN不需要事先指定簇的数量。

核心原理

DBSCAN通过在数据点周围划定一个$\varepsilon$半径，计算半径内的数据点数量来定义簇。其基本思想是密度足够高的区域形成簇，而密度较低的区域被认为是噪声。主要步骤如下：

  1. 定义核心点：对于数据集中每个点，如果在其ε半径内的点数大于或等于最小点数（minPts），则该点为核心点。

  2. 扩展簇：从核心点出发，将其邻域内的所有点（包括边界点和其他核心点）归入该簇，然后递归地将这些点的邻域内的点也归入该簇。

  3. 标记噪声点：如果一个点既不是核心点也不是边界点，则将其标记为噪声点。

核心公式

  1. 邻域定义：对于点$p$，其邻域定义为：
    \begin{equation}
    N_{\varepsilon}(p) = { q \in D|{\bf dist}(p,q)\leq \varepsilon }
    \end{equation}

其中 $D$ 是数据集，${\bf dist}(p,q)$ 是点 $p$ 和点 $q$ 之间的距离，通常使用欧氏距离。

  1. 核心点：如果点 $p$ 的邻域包含至少 minPts 个点，则 $p$ 是核心点：
    \begin{equation}
    |N_{\varepsilon}(p)| \geq {\bf minPts}
    \end{equation}

Python案例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_moons
from sklearn.cluster import DBSCAN

# 生成样本数据
X, y = make_moons(n_samples=500, noise=0.1, random_state=0)

# 使用DBSCAN进行聚类
dbscan = DBSCAN(eps=0.2, min_samples=5)
y_dbscan = dbscan.fit_predict(X)

# 绘制结果图
plt.scatter(X[:, 0], X[:, 1], c=y_dbscan, s=50, cmap='viridis')

# 标记噪声点
core_samples_mask = np.zeros_like(y_dbscan, dtype=bool)
core_samples_mask[dbscan.core_sample_indices_] = True
noise_mask = (y_dbscan == -1)
plt.scatter(X[noise_mask, 0], X[noise_mask, 1], c='red', s=50, label='Noise')

plt.title("DBSCAN Clustering")
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")
plt.legend()
plt.show()

  1. 生成数据：使用make_moons生成500个二维数据点，这些数据点形成两个半月形簇，并添加了一定的噪声。

  2. 训练模型：使用DBSCAN模型进行训练，设置$\varepsilon$为0.2，minPts为5。

  3. 预测与绘图：根据训练结果进行预测，并将数据点按簇绘制成不同颜色，同时用红色标注出噪声点。

大家可以清晰地看到DBSCAN的聚类效果，两个半月形的簇被正确识别出来，噪声点用红色标注。DBSCAN能够很好地处理复杂形状的簇，并且能够识别出数据中的噪声点。

高斯混合模型（Gaussian Mixture Model，GMM）

高斯混合模型（GMM）是一种概率模型，它假设所有的数据点是由若干个不同的高斯分布（即正态分布）组成的混合生成。GMM常用于聚类分析，特别适合于发现复杂数据分布中的潜在群体结构。

核心原理

GMM通过期望最大化（EM）算法来估计模型参数。EM算法是一个迭代优化算法，包括两个主要步骤：

  1. E步（Expectation）：计算给定当前参数下每个数据点属于每个高斯分布的概率（即责任）。

  2. M步（Maximization）：根据计算出的责任，重新估计每个高斯分布的参数（均值、方差和权重）。

核心公式

  1. 高斯混合模型的概率密度函数：
    \begin{equation}
    p(x) = \sum^{K}_{k = 1}\pi_k \mathcal{N}(x|\mu_k,\Sigma_k)
    \end{equation}

其中，$\pi_k$是第$k$个高斯分布的权重，$\mathcal{N}(x|\mu_k,\sum_k)$是第$k$个高斯分布的概率密度函数，其参数为均值$\mu_k$和协方差矩阵$\Sigma_k$。

  1. E步（计算责任）：
    \begin{equation}
    \gamma_{ik} = \frac{\pi_k\mathcal{N}(x|\mu_k,\Sigma_k)}{\sum_{j = 1}^{K}\pi_j \mathcal{N}(x|\mu_k,\Sigma_k)}
    \end{equation}

其中，$\gamma_{ik}$是数据点$x_i$属于第$k$个高斯分布的概率。

  1. M步（更新参数）：

更新权重：
\begin{equation}
\pi_k = \frac{N_k}{N}
\end{equation}

其中，$N_k = \sum_{i=1}^{N} \gamma_{ik}$是第$k$个高斯分布的有效样本数，$N$是数据点总数。

更新均值：
\begin{equation}
\mu_k = \frac{1}{N_k}\sum_{i=1}^{N} \gamma_{ik}x_i
\end{equation}

更新协方差矩阵：
\begin{equation}
\Sigma_k = \frac{1}{N_k}\sum_{i=1}^{N}\gamma_{ik}(x_i-\mu_k)(x_i-\mu_k)^{T}
\end{equation}

Python案例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.mixture import GaussianMixture

# 生成样本数据
X, y = make_blobs(n_samples=500, centers=3, random_state=42, cluster_std=0.60)

# 使用GMM进行聚类
gmm = GaussianMixture(n_components=3, random_state=42)
gmm.fit(X)
labels = gmm.predict(X)

# 绘制结果图
plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis')
plt.scatter(gmm.means_[:, 0], gmm.means_[:, 1], c='red', s=200, alpha=0.75, marker='X')  # 标记出中心点
plt.title("Gaussian Mixture Model Clustering")
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")
plt.show()

  1. 生成数据：使用make_blobs生成500个二维数据点，分为3个簇，并添加一定的噪声。

  2. 训练模型：使用GaussianMixture模型进行训练，设定簇数为3。

  3. 预测与绘图：根据训练结果进行预测，并将数据点按簇绘制成不同颜色，同时用红色标注出每个高斯分布的均值。

大家可以看到GMM的聚类效果，每个簇用不同的颜色表示，红色的X标记出每个高斯分布的中心。GMM能够很好地处理复杂的簇结构，并且可以估计每个簇的概率分布。

谱聚类（Spectral Clustering）

谱聚类（Spectral Clustering）是一种基于图论的聚类方法，通过利用数据点之间的相似性构建图，并利用图的谱（特征值和特征向量）进行聚类。它特别适合处理非凸形状的簇，并且能有效处理高维数据。

核心原理

谱聚类的核心思想是将数据点看作图中的节点，通过计算节点之间的相似性构建加权图，然后通过图的拉普拉斯矩阵的特征向量将数据映射到低维空间，再在低维空间中进行K均值聚类。

主要步骤如下：

  1. 构建相似度矩阵：计算数据点之间的相似度，构建相似度矩阵。

  2. 构建拉普拉斯矩阵：根据相似度矩阵构建图的拉普拉斯矩阵。

  3. 计算特征向量：计算拉普拉斯矩阵的前k个特征向量，将数据点映射到低维空间。

  4. 聚类：在低维空间中对映射后的数据点进行K均值聚类。

核心公式

  1. 相似度矩阵$W$：
    \begin{equation}
    W_{ij} = \exp\left(-\frac{\Arrowvert x_i-x_j \Arrowvert^2}{2\sigma^2}\right) \quad {\bf if} \quad \Arrowvert x_i - x_j \Arrowvert \leq \varepsilon, {\bf else} \quad 0
    \end{equation}

其中，$\sigma$是高斯核函数的参数，$\varepsilon$是邻域范围。

  1. 度矩阵$D$：
    \begin{equation}
    D_{ij} = \sum_{j} W_{ij}
    \end{equation}

  2. 拉普拉斯矩阵$L$：
    \begin{equation}
    L = D-W
    \end{equation}

或归一化的拉普拉斯矩阵：
\begin{equation}
L_{sym} = D^{-1/2}LD^{-1/2}
\end{equation}

  1. 特征向量：计算$L$或$L_{sym}$的前$k$个特征向量$u_1,u_2,…,u_k$。

  2. 聚类：将每个数据点$x_i$映射到特征向量组成的低维空间，再对这些向量进行$K$均值聚类。

Python案例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_moons
from sklearn.cluster import SpectralClustering

# 生成样本数据
X, y = make_moons(n_samples=500, noise=0.1, random_state=42)

# 使用谱聚类进行聚类
spectral = SpectralClustering(n_clusters=2, affinity='nearest_neighbors', assign_labels='kmeans', random_state=42)
labels = spectral.fit_predict(X)

# 绘制结果图
plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis')
plt.title("Spectral Clustering")
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")
plt.show()

  1. 生成数据：使用make_moons生成500个二维数据点，形成两个半月形簇，并添加一定的噪声。

  2. 训练模型：使用SpectralClustering模型进行训练，设定簇数为2，使用最近邻相似度计算，标签分配使用K均值。

  3. 预测与绘图：根据训练结果进行预测，并将数据点按簇绘制成不同颜色。

可以看到谱聚类的效果，两个半月形的簇被正确识别出来，每个簇用不同的颜色表示。谱聚类通过利用数据点之间的相似性构建图结构，并通过图的谱（特征向量）进行聚类，有效地处理了复杂形状的簇结构。