第10.8节 聚类内部评估指标
各位朋友大家好，欢迎来到月来客栈，我是掌柜空字符。
本期推送内容目录如下，如果本期内容对你有所帮助，欢迎点赞、转发支持掌柜！
10.8 聚类内部评估指标 10.8.1 Silhouette Coefficient Index 10.8.2 Calinski-Harabasz Index 10.8.3 Davies-Bouldin Index 10.8.4 小结 引用
10.8 聚类内部评估指标
在第10.6节内容中，笔者详细介绍了聚类算法中几种常见的评估指标，包括纯度（准确率）、精确率、召回率、兰德系数和F值等。虽然这些评价指标都能很好的评估聚类结果的优劣，但是它们都有着一个共同的缺点，需要知道训练样本的真实标签（Ground Truth），而这在现实情况中却很难做到。
对于聚类结果的评价方法一般可以分为内部评估法（Internal Evaluation）与外部评估方法（External Evaluation）。所谓外部评估方法是指在知道真实标签的情况下来评估聚类结果的好坏，即笔者上面提到几种评估指标。一般来说在做论文，或者是有少量的标注数据时，都可以用外部评估法选择一个相对最优的聚类模型，然后再将其应用到其它未被标记的数据中。所谓内部评估法则是不借助于外部信息，仅仅只是依靠聚类结果和样本本身的属性（空间分布特性等）来进行评估的办法。常见的内部评估法有轮廓系数（Silhouette Coefficient）、Calinski-Harabasz Score和Davies-Bouldin Score等。在接下来的这节内容中，笔者将会分别就这3种常见的内部评价指标进行介绍。
10.8.1 Silhouette Coefficient Index
所谓轮廓系数（Silhouette Coefficient）[1]，它本质上衡量的是每个样本点到其簇内样本的距离与其最近簇结构之间距离的比值。如果该比值越小，则说明该样本点所在的簇结构与其最近簇结构之间的距离越远，从而说明此时的聚类结果越好。
1. 原理
在介绍轮廓系数原理之前，首先需要明白的是对于聚类结果的轮廓系数，它实际上是每个样本点轮廓系数的平均值，也就是说聚类样本中的每个样本点对应都有一个轮廓系数值，而总的轮廓系数则是所有样本点轮廓系数的平均值。
如图10-12所示，数据样本中一共包含有3个簇结构，对于最左边这个簇中（圆形）的样本点来说，距离其最近的簇为中间（方形）这个簇。现在定义样本到簇中每个样本距离的均值为，到最近簇中每个样本距离的均值为，即此时根据图10-12所示的结果有：
这里需要注意的一点是，对于同一个簇中的每个样本点来说，距离自己最近的簇可能并不是同一个；同时，在寻找距离当前样本点最近的簇结构时，计算的是当前样本点到各个簇中心的最短距离，而不是计算当前样本点所在簇的簇中心到其它每个簇中心的最短距离。
此时，样本的轮廓系数定义为
从式(10.35)可以看出，当越小而越大时，此时对应的情况是样本所在的簇中所有样本点之间距离都比较近（簇内距离较小），样本所在的簇距离其最近簇的距离较远（因为此时样本所在的簇中簇内距离较小，所以样本到其它簇的距离可以近似地看做样本所在簇到其它簇之间的距离），所以此时也就越接近于1，即聚类效果越好。
类似地，当越小而越大时，此时对应的情况是样本所在的簇中所有样本点之间距离都比较远（簇内距离较大），样本所在的簇距离其最近簇的距离较近，所以此时也就越接近于-1，即聚类效果越差。
最后，当样本所在的簇中只有样本一个样本时，那么对于样本到底属于哪个簇就变得不清楚了，此时可以粗略地认为属于离它最近的簇中，所以此时有，即表示一种中立的情况。
由此可知的取值范围应该是，且进一步式(10.35)可以改为为
最终，整个聚类结果的轮廓系数结果为
其中表示样本点的个数。
2. 实现
在介绍完轮廓系数的原理和计算方式后，下面笔者就来介绍如何编程实现这一过程。具体实现如下所示：
1 def get_silhouette_coefficient(X, labels):
2 n_clusters = np.unique(labels).shape[0]
3 s = []
4 for k in range(n_clusters): # 遍历每一个簇
5 index = (labels == k) # 取对应簇所有样本的索引
6 x_in_cluster = X[index] # 去对应簇中的所有样本
7 for sample in x_in_cluster: # 计算每个样本的轮廓系数
8 a = ((sample - x_in_cluster) ** 2).sum(axis=1)
9 a = np.sqrt(a).sum() / (len(a) - 1) # 去掉当前样本点与当前样本点的组合计数
10 nearest_cluster_id = None
11 min_dist2 = np.inf
12 for c in range(n_clusters): # 寻找距离当前样本点最近的簇
13 if k == c:
14 continue
15 centroid = X[labels == c].mean(axis=0)
16 dist2 = ((sample - centroid) ** 2).sum()
17 if dist2 < min_dist2:
18 nearest_cluster_id = c
19 min_dist2 = dist2
20 x_nearest_cluster = X[labels == nearest_cluster_id]
21 b = ((sample - x_nearest_cluster) ** 2).sum(axis=1)
22 b = np.sqrt(b).mean()
23 s.append((b - a) / np.max([a, b]))
24 return np.mean(s)
在上述代码中，第1行X和labels分别是原始训练集和聚类得到的标签；第2行是取聚类结果的簇数量；第4~6行是遍历每个簇并取对应中所有的样本；第7行开始是遍历每一个样本点来计算轮廓系数；第8~9行是计算当前样本距离其簇中每个样本点之间的平均距离，其中第9行分母减去1是因为要去掉当前样本与当前样本之间距离为0的情况；第10~19行是计算得到里当前样本点最近的簇，其中15行是计算簇中心；第20行是取距离当前样本点最近的簇对应的所有样本点；第21~22行是计算当前样本点到其最近簇中所有样本的平均距离；第23行是计算当前样本点对应的轮廓系数并保存；第24行是返回所有样本对应的轮廓系数的均值。
完成上述代码实现之后，便可以通过如下方式进行使用：
1 import numpy as np
2 from sklearn.cluster import KMeans
3 from sklearn.datasets import load_iris
4 from sklearn.metrics import silhouette_score
5
6 def test_silhouette_score():
7 x, y = load_iris(return_X_y=True)
8 model = KMeans(n_clusters=3)
9 model.fit(x)
10 y_pred = model.predict(x)
11 print(f"轮廓系数 by sklearn: {silhouette_score(x, y_pred)}")
12 print(f"轮廓系数 by ours: {get_silhouette_coefficient(x, y_pred)}")
13
14
15 if __name__ == '__main__':
16 test_silhouette_score()
上述代码运行结束后的结果如下：
1 轮廓系数 by sklearn: 0.5528190123564091
2 轮廓系数 by ours: 0.5528190123564101
10.8.2 Calinski-Harabasz Index
在介绍完轮廓系数后我们再来看第2种内部评价指标Calinski-Harabasz指数[2]。由于Calinski-Harabasz指数的本质是簇间距离与簇内距离的比值，且整体计算过程与方差计算方式类似，所以又将其称之为方差比准则。
为你认可的知识付费，欢迎订阅本专栏阅读更多优质内容！
1. 原理
方差比准则的计算过程相较于轮廓系数更加简单，我们只需要计算出所有簇的簇内距离和然后再比上所有簇的簇间距离和即可，原理如图10-13所示。
如图10-13所示一共包含有3个聚类簇，每个簇中均有6个样本点，且每个簇内部的白色圆点表示对应簇的簇中心。表示第个簇对应的簇内距离和，即每个样本点减去簇中心的平方和；表示第个簇的簇中心到全局簇中心（即所有样本点的中心）的加权距离，而方差比则是所有簇内距离和与簇间距离和的比值，即
其中表示簇数量，表示第个簇中的所有样本，表示第个簇对应的簇中心，表示第个簇中样本的数量，表示所有样本的数量，表示全局簇中心。
值得一提的是，由于每个簇的簇中心到全局簇中心的距离只会计算一次，在计算方差比时为了消除与簇内距离上的差异所以还乘上了每个簇对应的样本总数；同时，在计算方差时采用是方差的无偏估计。
根据式(10.38)可知，如果越大，则表示相对于越大，此时的簇间距离便远大于簇内距离，即簇与簇之间相距较远聚类结果好；相反，如果越小，则表示簇与簇之间相距较近聚类效果较差。同时，可以看出方差比的取值范围应该是。
2. 实现
在清楚了方差比的计算原理后，实现起来就比较容易了，整体实现代码如下：
1 def get_calinski_harabasz(X, labels):
2 n_samples = X.shape[0]
3 n_clusters = np.unique(labels).shape[0]
4 betw_disp = 0. # 所有的簇间距离和
5 within_disp = 0. # 所有的簇内距离和
6 global_centroid = np.mean(X, axis=0) # 全局簇中心
7 for k in range(n_clusters): # 遍历每一个簇
8 x_in_cluster = X[labels == k] # 取当前簇中的所有样本
9 centroid = np.mean(x_in_cluster, axis=0) # 计算当前簇的簇中心
10 within_disp += np.sum((x_in_cluster - centroid) ** 2)
11 betw_disp += len(x_in_cluster) * np.sum((centroid - global_centroid) ** 2)
12
13 return (1. if within_disp == 0. else
14 betw_disp * (n_samples - n_clusters) /
15 (within_disp * (n_clusters - 1.)))
在上述代码中在上述代码中，第1行X和labels分别是原始训练集和聚类得到的标签；第6行是计算全局簇中心；第7~9行是开始遍历每一个簇并计算得到每个簇对应的簇中心；第10行是累计每个簇对应簇内距离的总和；第11行是累计每个簇中心到全局中心的距离总和；第13~15行则是返回最终计算得到方差比结果。
完成上述代码实现之后，便可以通过如下方式进行使用：
1 def test_calinski_harabasz_score():
2 x, y = load_iris(return_X_y=True)
3 model = KMeans(n_clusters=3)
4 model.fit(x)
5 y_pred = model.predict(x)
6 print(f"方差比 by sklearn: {calinski_harabasz_score(x, y_pred)}")
7 print(f"方差比 by ours: {get_calinski_harabasz(x, y_pred)}")
8
9 if __name__ == '__main__':
10 test_calinski_harabasz_score()
上述代码运行结束后的结果如下：
1 方差比 by sklearn: 561.62775662962
2 方差比 by ours: 561.62775662962
10.8.3 Davies-Bouldin Index
在介绍完前面两种聚类内容部评价指标后我们再来看第3种评价方法Davies-Bouldin Index（DB指数）[3]。DB指数的核心思想是计算每个簇与之最相似簇之间相似度，然后再通过求得所有相似度的平均值来衡量整个聚类结果的优劣。如果簇与簇之间的相似度越高（DB指数偏高），也就说明簇与簇之间的距离越小（直观上理解只有相距越近的事物才会约相似），那么此时聚类结果就越差，反之亦然。
1. 原理
在DB指数中，簇与簇之间的相似度被定义为两个簇的簇内直径和与簇间距离的比值，即图10-14中与的比值。
从公式(10.40)可以看出，对于第个簇的相似度来说，其等于中的最大值；而DB指数则是所有簇相似度的平均值。虽然计算过程很清晰明了也容易实现，但是为什么要取最大值呢，有什么含义？如图10-14所示，对于第个簇来说其相似度值是还是呢？
根据式(10.39)可知，在分别计算和时是一直保持不变的。同时，如果要取得最大值，那么要尽可能大，而则要尽可能小。想一想大，小又意味着什么呢？意味着簇中的样本点分布离散，且簇与簇之间相距较近（相似度较高），对应的便是聚类结果较差。因此，越大，反映的便是这两个簇相距越近，越相似。所以，DB指数越大反映的是聚类结果越差，越小则是聚类结果越好。可以看出，DB指数有一点先求最大，再找最小的味道。
最后，根据式(10.39)和式(10.40)的分析可知，DB指数的取值范围在之间，越小表示聚类结果越好。
2. 实现
在清楚了方差比的计算原理后，实现起来就比较容易了，整体实现代码如下：
1 def get_davies_bouldin(X, labels):
2 n_clusters = np.unique(labels).shape[0]
3 centroids = np.zeros((n_clusters, len(X[0])), dtype=float)
4 s_i = np.zeros(n_clusters)
5 for k in range(n_clusters): # 遍历每一个簇
6 x_in_cluster = X[labels == k] # 取当前簇中的所有样本
7 centroids[k] = np.mean(x_in_cluster, axis=0) # 计算当前簇的簇中心
8 s_i[k] = pairwise_distances(x_in_cluster, [centroids[k]]).mean() #
9 centroid_distances = pairwise_distances(centroids) # [K,K]
10 combined_s_i_j = s_i[:, None] + s_i # [K,k]
11 centroid_distances[centroid_distances == 0] = np.inf
12 scores = np.max(combined_s_i_j / centroid_distances, axis=1)
13 return np.mean(scores)
在上述代码中，第1行X和labels分别是原始训练集和聚类得到的标签；第3~4行是分别初始化一个簇中心矩阵和簇内直径向量；第5~7行是遍历每一个簇，并计算得到对应的簇中心；第8行是同时计算当前簇中所有样本x_in_cluster到对应簇中心centroids[k]的距离，然后取平均得到；第9行是计算centroids中任意两个簇中心之间的距离，得到的将是一个形状为[K,K]的对称矩阵，对角线均为0（自己与自己的距离）；第10行是构造得到一个形状同为[K,K]的矩阵便于后续计算。
例如：
1 centroid_distances = [[0. 3.35693455 5.01756852]
2 [3.35693455 0. 1.7971818 ]
3 [5.01756852 1.7971818 0. ]]
4 combined_s_i_j = [[0.96341047 1.21985761 1.20154379]
5 [1.21985761 1.47630474 1.45799092]
6 [1.20154379 1.45799092 1.4396771 ]]
上述结果中，centroid_distances[i][j]表示第i个簇与第j个簇之间的距离；combined_s_i_j[i][j]表示s_i[i]+s_i[j]的结果。
同时，第11行是将centroid_distances中对角线的0填充为无穷大；第12行是计算得到所有的的情况并取最大值；第13行则是返回最后DB的计算结果。
完成上述代码实现之后，便可以通过如下方式进行使用：
1 def test_davies_bouldin_score():
2 x, y = load_iris(return_X_y=True)
3 model = KMeans(n_clusters=3)
4 model.fit(x)
5 y_pred = model.predict(x)
6 print(f"db_score by sklearn: {davies_bouldin_score(x, y_pred)}")
7 print(f"db_score by ours: {get_davies_bouldin(x, y_pred)}")
8
9 if __name__ == '__main__':
10 test_davies_bouldin_score()
上述代码运行结束后的结果如下：
1 db_score by sklearn: 0.6619715465007511
2 db_score by ours: 0.6619715465007511
10.8.4 小结
在这节内容中，笔者首先详细介绍了轮廓系数的思想原理与实现过程，明确了聚类结果的轮廓系数实质上就是所有样本点轮廓系数的均值；然后介绍了Calinski-Harabasz指数的原理与实现，其大致上可以看做是簇间距离与簇内距离的比值，因其计算过程与方差类似所以又被称之为方差比；最后介绍了Davies-Bouldin评价指标的原理与实现，DB指数可以看做是通过簇与簇之间相似度的大小来衡量聚类结果的优劣。在下一节内容中，笔者将会详细介绍如果通过轮廓系数等方法来选择聚类算法中的参数K值。
引用
[1] Rousseeuw P J. Silhouettes: a graphical aid to the interpretation and validation of cluster analysis[J]. Journal of computational and applied mathematics, 1987, 20: 53-65.
[2] Caliński T, Harabasz J. A dendrite method for cluster analysis[J]. Communications in Statistics-theory and Methods, 1974, 3(1): 1-27.
[3] Davies, David L.; Bouldin, Donald W. (1979). “A Cluster Separation Measure” IEEE Transactions on Pattern Analysis and Machine Intelligence. PAMI-1 (2): 224-227.
[4] https://en.wikipedia.org/wiki/Davies%E2%80%93Bouldin_index
[5] https://scikit-learn.org/stable/modules/clustering.html#clustering-evaluation