数据分析 层次 聚类 算法 无监督
开始之前，我们先理解几个变量：
- X 实验样本(n 乘 m的数组)
- n - 样本数量
- m - 样本特征数量
- Z - 集群关系数组（包含层次聚类信息）
- k - 集群数量
模块引入
from matplotlib import pyplot as plt
from scipy.cluster.hierarchy import dendrogram, linkage
import numpy as np
%matplotlib inline
np.set_printoptions(precision=5, suppress=True)
生成实验样本
如果你已经有自己的实验样本了，那么你可以选择跳过此步骤。但是要注意的是，你必须确保你的实验样本是一个 NumPy 矩阵 X，包含 n 个实验样本以及 m 个样本特征，即 X.shape == (n, m)。
# 生成两个集群: 集群 a 有100个数据点, 集群 b 有50个数据点:
np.random.seed(1029)
a = np.random.multivariate_normal([10, 0], [[3, 1], [1, 4]], size=[100,])
b = np.random.multivariate_normal([0, 20], [[3, 1], [1, 4]], size=[50,])
X = np.concatenate((a, b),)
plt.figure(figsize=(25, 10))
plt.scatter(X[:, 0], X[:, 1])
plt.show()
开始聚类
# 建立集群关系数组
Z = linkage(X, 'ward')
# 使用 'ward' 是一个不错的选择。
# 当然还有其他通用的距离算法，例如 'single'、'complete'，'average'。我们可以使用其他的方法进行对比已检验更适合的。
在这一行代码中，根据 SciPy 的 linkage 模块文档所述，'ward' 是一种用于计算集群之间距离的方法。'ward' 表明 linkage() 方法会使用离差平方和算法。
# 你可以使用 `cophenet()` 计算同表象相关系数来判断集群的性能。
# 这个方程非常简单的比对了所有实验样本之间的距离和聚类之后的样本距离。
# 如果这个值越接近于1，集群结果就更加完整的保存了实验样本间的实际距离。在我们的实验中，这个值非常接近于1：
from scipy.cluster.hierarchy import cophenet
from scipy.spatial.distance import pdist
c, coph_dists = cophenet(Z, pdist(X))
print(c)
# 0.9840295755
具体聚类的理解
不管你使用了何种 method 和 metric，linkage() 方程都会利用这些方法计算集群间的距离，并在每一步合并距离最短的两个集群。linkage() 方程会返回一个长度为 n - 1 的数组，包含在每一步合并集群的信息：
print(Z[0])
# [ 124. 141. 0.07372 2. ]
# 这个矩阵的每一行的格式是这样的
# [下标、下标2 、距离、 样本数量]
# [idx, idx2, dist, sample_count]
在第一步中，linkage 算法决定合并集群124和集群141，因为他们之间的距离为0.07372，为当前最短距离。这里的124和141分别代表在原始样本中的数组下标。在这一步中，一个具有两个实验样本的集群诞生了。
print(Z[1])
# [ 131. 139. 0.0753 2. ]
在第二步中，linkage 算法决定合并集群131和集群139，因为他们之间的距离为0.0753，为当前最短距离。在这一步中，另一个具有两个实验样本的集群诞生了。
直到这一步，集群的下标和原本实验样本的下标是互相对应的，但是记住，我们只有150个实验样本，所以我们一共只有下标0至149。让我们来看看前20步：
print(Z[:20])
直到前7步，此算法都在直接合并原有实验样本数据，也不难看出每一步的最短距离都在单调递增。在第8步中，此算法决定合并集群56和集群153。是不是很奇怪，我们只有150个数据点，何来的下标153？此算法中，任何 idx >= len(X)的下标都指向在 Z[idx - len(X)] 中建立的集群。所以，第8步中，样本62与在第四步 Z[153 - 150] 中合并的集群（样本10和样本18）进行了合并。
让我们来看看这些数据点的坐标是不是反应了这一情况：
print(X[[10, 18, 56]])
# 打印结果：
[[ 8.91446 -1.95095]
[ 8.95663 -2.02432]
[ 8.96876 -1.86035]]
让我们再画图看看是不是这样：
idxs = [10, 18, 56]
plt.figure(figsize=(10, 8))
plt.scatter(X[:, 0], X[:, 1])
plt.scatter(X[idxs, 0], X[idxs, 1], c='r')
plt.show()
建立树状图
# 树状图展现了层次聚类中集群合并的顺序以及合并时集群间的距离。
plt.figure(figsize=(50, 10))
plt.title('Hierarchical Clustering Dendrogram')
plt.xlabel('sample index')
plt.ylabel('distance')
dendrogram(Z, leaf_rotation=90., leaf_font_size=8)
plt.show()
在树状图中，x 轴上的标记代表数据点在原有实验样本 X 中的下标，y 轴上表明集群间的距离。图像中，横线所在高度表明集群合并时的距离。
我们可以看到在 y 轴20处我们只有两个集群，但在最后的合并后，集群间的距离直线上升至180左右。让我们看看最后4步的集群间距离：
print(Z[-4:, 2])
# 打印结果
[ 15.2525 16.79548 21.75623 185.07009]
这样的大距离跨度表明了在这一步中虽然有两个集群被合并了，但是它们可能不应该被合并。换句话说，它们可能本来就不属于同一个集群，它们就应该属于两个不同的集群。
在树状图中，我们可以看到绿色集群只包含了大于等于100的下标，而红色集群只包含了小于100的下标。看来我们的算法成功找出了实验样本中的两个集群。
缩减树状图
上一步得到的树状图非常大，然而它却仅仅包含了150个数据点。现实情况中，数据点可能更多。让我们来看看 dendrogram() 方程的其他参数：
plt.title('Hierarchical Clustering Dendrogram (truncated)')
plt.xlabel('sample index')
plt.ylabel('distance')
# 我们可以通过show_contracted=True参数来控制是否显示合并子集个数
# p=12 表示显示最后12次合并的图形
dendrogram(Z, truncate_mode='lastp', p=12, show_leaf_counts=False, leaf_rotation=90., leaf_font_size=12, show_contracted=True)
plt.show()
在树状图上显示距离
def fancy_dendrogram(*args, **kwargs):
max_d = kwargs.pop('max_d', None)
if max_d and 'color_threshold' not in kwargs:
kwargs['color_threshold'] = max_d
annotate_above = kwargs.pop('annotate_above', 0)
ddata = dendrogram(*args, **kwargs)
if not kwargs.get('no_plot', False):
plt.title('Hierarchical Clustering Dendrogram (truncated)')
plt.xlabel('sample index or (cluster size)')
plt.ylabel('distance')
for i, d, c in zip(ddata['icoord'], ddata['dcoord'], ddata['color_list']):
x = 0.5 * sum(i[1:3])
y = d[1]
if y > annotate_above:
plt.plot(x, y, 'o', c=c)
plt.annotate("%.3g" % y, (x, y), xytext=(0, -5), textcoords='offset points', va='top', ha='center')
if max_d:
plt.axhline(y=max_d, c='k')
return ddata
# 方法调用
fancy_dendrogram(Z, truncate_mode='lastp', p=12, leaf_rotation=90., leaf_font_size=12, show_contracted=True, annotate_above=10)
plt.show()
临界距离决定集群数量
就像我们在之前解释的那样，大距离跨度通常是我们所感兴趣的地方。在我们的例子中，我们把距离临界值设为50，因为在这里距离跨度非常明显：
# max_d = 50
# 当临界值为16的时候，我们有4个集群。显然这是不符合我们实际的临界值。
max_d = 16
fancy_dendrogram(Z, truncate_mode='lastp', p=12, leaf_rotation=90., leaf_font_size=12, show_contracted=True, annotate_above=10, max_d=max_d)
plt.show()
获取集群信息
- 如果我们已经通过树状图知道了最大临界值，我们可以通过以下代码获得每个实验样本所对应的集群下标：
from scipy.cluster.hierarchy import fcluster
max_d = 50
clusters = fcluster(Z, max_d, criterion='distance')
print(clusters)
# 打印结果
[2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1]
- 如果我们已经知道最终会有2个集群，我们可以这样获取集群下标：
k = 2
fcluster(Z, k, criterion='maxclust')
# 结果展示
[2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1]
可视化集群
如果你的实验样本特征数量很少，你可以可视化你的集群结果：
plt.figure(figsize=(10, 8))
plt.scatter(X[:, 0], X[:, 1], c=clusters, cmap='prism')
plt.show()
原文地址：原作
