轮廓系数、方差比、DB指数（三种常见的聚类内部评价指标）-526互联

1 引言

在之前的一篇文章（https://www.cnblogs.com/emanlee/p/17742869.html）中掌柜详细介绍了聚类算法中几种常见的评估指标，包括纯度（准确率）、精确率、召回率、兰德系数和F值等。虽然这些评价指标都能很好的评估聚类结果的优劣，但是它们都有着一个共同的缺点，需要知道训练样本的真实标签（ground truth），而这在现实情况中却很难做到。

对于聚类结果的评价方法一般可以分为内部评估法（internal evaluation）与外部评估方法（external evaluation）。所谓外部评估方法是指在知道真实标签（ground truth ）的情况下来评估聚类结果的好坏，也就是上面提到几种评估指标。一般来说在做论文，或者是有少量的标注数据时，都可以用外部评估法选择一个相对最优的聚类模型，然后再应用到其它未被标记的数据中。所谓内部评估法则是不借助于外部信息，仅仅只是依靠聚类结果和样本本身的属性来进行评估的办法。常见的内部评估法有轮廓系数（Silhouette Coefficient）、Calinski-Harabasz Score和Davies-Bouldin score等。在接下来的这篇文章中，将会分别就这3中内部评价指标进行介绍。

2 Silhouette Coefficient Index

所谓轮廓系数（Silhouette Coefficient）[1]，它本质上衡量的是每个样本点到其簇内样本的距离与其最近簇结构之间距离的比值。如果该比值越小，则说明该样本点所在的簇结构与其最近簇结构之间的距离越远，从而说明此时的聚类结果越好。

2.1 原理

在介绍轮廓系数原理之前，首先需要明白的是对于聚类结果的轮廓系数，它实际上是每个样本点轮廓系数的平均值。这也就是说聚类样本中的每个样本点对应都有一个轮廓系数值，而总的轮廓系数则是所有样本点轮廓系数的平均值。

图 1. 轮廓系数原理图

如图1所示，数据样本中一共包含有3个簇结构，对于最左边这个簇中（圆形）的样本点 i 来说，距离其最近的簇为中间（方形）这个簇。现在定义样本 i

这里需要注意的一点是，对于同一个簇中的每个样本点来说，距离自己最近的簇可能并不是同一个；同时，在寻找距离当前样本点最近的簇结构时，计算的是当前样本点到各个簇中心的最短距离，而不是计算当前样本点所在簇的簇中心到其它每个簇中心的最短距离。

此时，样本 i 定义为

从式(2)

s(i) 此时也就越接近于1，即聚类效果越好。

类似地，当 b(i)

此时 s(i) 也就越接近于-1，即聚类效果越差。

最后，当样本 i

s(i) =0 表示一种中立的情况。

由此可知 s(i)

可以改为为

最终，整个聚类结果的轮廓系数结果为

其中 n 表示样本点的个数。

2.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行是返回所有样本对应的轮廓系数的均值。

所有完整代码均可从本仓库获取：https://github.com/moon-hotel/MachineLearningWithMe

完成上述代码实现之后，便可以通过如下方式进行使用：

 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

3 Calinski-Harabasz Index

在介绍完轮廓系数后我们再来看第2种内部评价指标Calinski-Harabasz指数[2]。由于Calinski-Harabasz指数的本质是簇间距离与簇内距离的比值，且整体计算过程与方差计算方式类似，所以又将其称之为方差比准则。

如果还有不知道怎么有效入门机器学习的朋友，掌柜强烈推荐阅读下面这篇文章，里面详细介绍了掌柜所推荐的通过五个步骤和三个阶段的方法来高效学习机器学习的路线，同时也按顺序整理了客栈了所有与机器学习相关的文章。

月来客栈：初学者如何入门机器学习？49 赞同 · 14 评论文章

3.1 原理

方差比准则的计算过程相较于轮廓系数更加简单，我们只需要计算出所有簇的簇内距离和然后再比上所有簇的簇间距离和即可，原理如图2所示。

如图2所示一共包含有3个聚类簇，每个簇中均有6个样本点，且每个簇内部的白色圆点表示对应簇的簇中心。

其中K

表示全局簇中心。

值得一提的是，由于每个簇的簇中心到全局簇中心的距离只会计算一次，在计算方差比时为了消除与簇内距离上的差异所以还乘上了每个簇对应的样本总数；同时，在计算方差时采用是方差的无偏估计。

根据式

。

3.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

4 Davies-Bouldin Index

在介绍完前面两种聚类内容部评价指标后我们再来看第3种评价方法Davies-Bouldin Index（DB指数）[3]。DB指数的核心思想是计算每个簇与之最相似簇之间相似度，然后再通过求得所有相似度的平均值来衡量整个聚类结果的优劣。如果簇与簇之间的相似度越高（DB指数偏高），也就说明簇与簇之间的距离越小（直观上理解只有相距越近的事物才会约相似），那么此时聚类结果就越差，反之亦然。

4.1 原理

在DB指数中，簇与簇之间的相似度被定义为两个簇的簇内直径和与簇间距离的比值，即图3中

的比值。

如图3所示一共包含有3个聚类簇，每个簇中均有6个样本点，且每个簇内部的白色圆点表示对应簇的簇中心。

之间的相似度为

进一步，DB指数的计算公式为

其中

表示总的簇数量。

从公式

呢？

根据式

越大，反映的便是这两个簇相距越近，越相似。所以，DB指数越大反映的是聚类结果越差，越小则是聚类结果越好。可以看出，DB指数有一点先求最大，再找最小的味道。

最后，根据式

之间，越小表示聚类结果越好。

4.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]的矩阵便于后续计算。

所有完整代码均可从本仓库获取：https://github.com/moon-hotel/MachineLearningWithMe

例如：

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

5 总结

在这篇文章中，掌柜首先详细介绍了轮廓系数的思想原理与实现过程，明确了聚类结果的轮廓系数实质上就是所有样本点轮廓系数的均值；然后介绍了Calinski-Harabasz指数的原理与实现，其大致上可以看做是簇间距离与簇内距离的比值，因其计算过程与方差类似所以又被称之为方差比；最后介绍了Davies-Bouldin评价指标的原理与实现，DB指数可以看做是通过簇与簇之间相似度的大小来衡量聚类结果的优劣。

REF

https://zhuanlan.zhihu.com/p/530944459

系数轮廓coefficient silhouette