11.9 聚类K值选取与分析#
在前面几节内容中，我们陆续介绍了3种常见的聚类算法原理及其实现 、4种常见的聚类外部评价指标和3种常见的聚类内部评价指标，对于聚类算法的主体内容算是介绍得差不多了，但是还遗留了最后一个问题——K值的选取。通常来说在实际场景中数据集的标签都是未知的，相反恰好需要通过聚类结果来辅助将各个簇的样本区分开。那到底应该如何选取K值呢？
经过上一节内容的介绍，有读者可能会说，以不同的K值为超参数，以某种内部评价指标为标准，通过交叉验证来选择最佳K值不就行了吗？虽然这种想法看起来有道理，但实际上却行不通。以轮廓系数为例，根据轮廓系数的原理可知我们并不能够严格得出在最佳K值的情况下轮廓系数同样也能够取得最大值。因此，我们常常就需要借助一些其它手段来对K值进行分析。
在聚类分析中，有两种常见的方法用于K值的选择，分别是肘部法和轮廓系数法，下面分别逐一进行介绍。
11.9.1 K值分析肘部法原理#
根据11.1节内容的介绍我们知道，基于Kmeans聚类框架下的聚类算法在最小化目标函数时本质上就是最小化整个簇内距离。根据簇内距离的定义可知，当K值越小时那么簇内距离便会越大；而当K值越大时那么簇内距离便会越小，极端一点在K值等于样本数量时，那么此时的簇内距离可以取到最小值0。
进一步，根据分析可知，在K值由小变大的过程中，随着K值的增大簇内距离便会逐步减小，但是当K值取得最优解后簇内距离便不会出现明显的降幅。此时可以想象这么一个场景，假设某数据集以不同的K值进行聚类处理并同时计算得到对应的簇内距离和，再以不同的K值为横坐标，簇内距离为纵坐标进行可视化，便可以得到类似如图11-19所示的结果。
如图11-19所示，可以发现刚开始随着K值的增大簇内距离开始急剧减小，当K值大于4之后簇内距离的减少幅度便开始明显降低，因此便可以经验性地得出4是该数据集对应K值的最优解。之所以可以凭借这样的经验进行判断是因为，当K值小于最优解时K值每增加1也就代表多了1个簇结构，整个簇内距离和自然会得到大幅下降；而当K值大于最优解之后K值每增加1也就仅仅只是将原本一个正常的簇结构一分为二了，因此即便整个簇内距离和有所降低，但是幅度也不会太大。
基于这样的经验性准则便可以根据类似图11-19中的方法来进行K值的选取。同时，由于图11-19中的变化曲线类似于我们手臂的肘部，因此该方法被称之为肘部法。
11.9.2 K值分析肘部法实现#
在介绍完肘部法的分析原理以后再来看如何实现这一可视化过程，完整示例代码可参见 AllBooKCode/Chapter11/C12_elbow_analysis.py 文件。这里首先需要定义一个函数来根据聚类结果计算簇内距离和并进行可视化，示例代码如下：
1 def elbow_analysis(X, range_n_clusters, all_cluster_labels, all_centers):
2 all_dist = []
3 for n, n_clusters in enumerate(range_n_clusters):
4 cluster_labels = all_cluster_labels[n]
5 centers = all_centers[n]
6 dist = 0
7 for i in range(n_clusters): # 遍历每一个簇，计算当前簇的簇内距离
8 x_data = X[cluster_labels == i]
9 tmp = np.sum((x_data - centers[i]) ** 2, axis=1)
10 dist += np.sum(np.sqrt(tmp)) # 累计当前聚类结果下所有簇的簇内距离和
11 all_dist.append(dist)
12 plt.plot(range_n_clusters, all_dist) # 绘制肘部曲线
13 plt.scatter(range_n_clusters, all_dist) # 绘制各个点
14
15 for i in range(len(range_n_clusters)): # 在图上进行K值标记
16 plt.annotate(f"k = {range_n_clusters[i]}",
17 xy=(range_n_clusters[i], all_dist[i]), fontsize=14,
18 xytext=(range_n_clusters[i] + 0.1, all_dist[i]))
19 plt.xlim(range_n_clusters[0] - 0.5, range_n_clusters[-1] + 0.8) # 调整范围
20 plt.ylim(all_dist[-1] * 0.9, all_dist[0] + all_dist[-1] * 0.1)
21 plt.yticks([]) # 去掉y轴上的刻度显示
22 plt.xlabel("K", fontsize=13)
23 plt.ylabel("distance", fontsize=13)
24 plt.show()在上述代码中，第1行X表示原始的训练样本为一个 [n_samples, n_features]的二维矩阵；range_n_clusters是K值的取值范围，为一个普通的列表；all_cluster_labels是一个普通列表，每个元素为一个一维向量，即某一K值下的聚类标签；all_centers也是一个普通的列表，每个元素为一个[n_clusters,n_features]的二维矩阵，即某一K值下的聚类簇中心。第3~5行开始遍历不同的K取值，并取对应的聚类结果和聚类簇中心。第7~10行是遍历每一个簇并累加得到所有簇的簇内距离和。第12~13行是对计算得到的结果进行可视化。第15~18行是在图中进行标记。第20~24行是调整坐标轴范围等。
最后，可以通过如下代码来根据不同的K值进行聚类，并调用上述函数对结果进行可视化：
1 if __name__ == '__main__':
2 X, y = make_blobs(n_samples=500, n_features=2, centers=4, cluster_std=1,
3 center_box=(-10.0, 10.0), shuffle=True, random_state=1)
4 range_n_clusters = [2, 3, 4, 5, 6]
5 all_cluster_labels, all_centers = [], []
6 for n_clusters in range_n_clusters: # 对不同K值进行聚类处理
7 clusterer = KMeans(n_clusters=n_clusters, random_state=10)
8 cluster_labels = clusterer.fit_predict(X)
9 centers = clusterer.cluster_centers_
10 all_cluster_labels.append(cluster_labels)
11 all_centers.append(centers)
12 elbow_analysis(X, range_n_clusters, all_cluster_labels, all_centers)在上述代码中，第2~3行是得到一个人造数据集。第6~11行是根据不同K值对数据集进行聚类处理并保存对应的聚类结果。第12行则是对结果进行可视化。
进一步，根据类似的做法，还可以对一些常见的数据集进行簇内距离随K值变化的曲线图可视化，如图11-20所示。
如图11-20所示，在簇内距离随K值变化的曲线中我们利用肘部法能够很容易地判断前3个数据集（synthetic、iris和wine）的K值。但是对于最后一个数据集digits（数据集原本一共有10个类别，这里只取了其中若干个类别）来说，K值应该取多少呢？此时还可以借助轮廓系数法来进行分析。
11.9.3 K值分析轮廓系数法原理#
在第11.9.1节内容中，我们介绍了第1种用于K值选取的肘部法。虽然该方法一定程度上能够帮助我们有效地对K值进行选取，但是在某些情况下依旧存在着肘部法失效的情况。如图11-20中最后一幅图所示，随着K值的增大簇内距离和并没有出现明显的骤降情况，而是稳步地在进行减少，此时利用肘部便无法选择正确的K值。下面，开始介绍第2种K值分析方法，轮廓系数法。
轮廓系数法的核心思想是观察每个簇中样本数量的分布情况，以及簇中每个样本点轮廓系数偏离整体轮廓系数的程度来进行K值的分析判断。根据轮廓系数分析方法，图11-19中所列的人造数据集经过不同K值聚类后便可得到如图11-21所示的分析结果。
如图11-21所示，从上到下分别为K取值为2、3、4、5和6时的轮廓系数分析法可视化结果。以最下方K=6为例，左边为不同簇对应的轮廓系数图，右边为簇分布的可视化结果。在左边部分，虚线表示当前K值下聚类结果的平均轮廓系数（也就是聚类结果的整体轮廓系数）；不同颜色的不规则填充图在$y$轴上的宽度表示该簇中对应样本的数量（越宽表示该簇中的样本越多），在$x$轴方向上的长度表示对应簇中每个样本点轮廓系数的分布情况（越长表示该簇中样本的轮廓系数值越大）。
轮廓分析法的基本准则是：①每个簇中样本的分布情况应该尽可能均匀，即每一条填充图的宽度应该尽可能接近；②每个簇中都应该存在部分样本的轮廓系数值大于平均轮廓系数的情况，即每一条填充图的长度都应该超过平均轮廓系数所在的虚线，且每条填充图长度要尽可能保持接近；③当准则①和②都差不多的时候，更多的时候选择平均轮廓系数较大的情况。
基于上述第1条分析准则，我们很容易将K=2和K=6时的情况给排除掉，因为此时填充图之间的宽度差异较大；进一步，根据第2条准则可以将K=3和K=5时的情况排除掉，因为此时存在个别簇中样本的轮廓系数值小于平均轮廓系数值的情况。最终，根据分析当K=4时为最优解，同时从图11-21右侧的可视化结果也可以明显看出不同K值情况下聚类结果。
11.9.4 K值分析轮廓系数法实现#
在介绍完轮廓分析法的原理后再来看如何实现图11-21中的可视化结果，以下完整示例代码可参见AllBooKCode/Chapter11/C14_silhouette_analysis.py 文件。想要实现图11-21中的可视化结果，最重要的便是弄清楚如何绘制两条曲线之间的填充部分，以下为图11-21可视化的完整示例代码：
1 if __name__ == '__main__':
2 X, y = make_blobs(n_samples=500, n_features=2, centers=4, cluster_std=1,
3 center_box=(-10.0, 10.0), shuffle=True, random_state=1)
4 range_n_clusters = [2, 3, 4, 5, 6]
5 for n_clusters in range_n_clusters:
6 fig, (ax1, ax2) = plt.subplots(1, 2) # 创建一个1行2列的子图
7 fig.set_size_inches(12, 5) # 画布大小
8 ax1.set_xlim([-0.1, 1]) # 设置x轴的范围（轮廓系数）
9 ax1.set_ylim([0, len(X) + (n_clusters + 1) * 10]) # 顶端的间隙
10 clusterer = KMeans(n_clusters=n_clusters, random_state=10)
11 cluster_labels = clusterer.fit_predict(X)
12 centers = clusterer.cluster_centers_
13 silhouette_avg = silhouette_score(X, cluster_labels) # 所有样本的轮廓系数均值
14 print(" 当 n_clusters = ", n_clusters, "时，轮廓系数为: ", silhouette_avg)
15 sample_silhouette_values = silhouette_samples(X, cluster_labels)
16 y_lower = 10
17 for i in range(n_clusters): # 遍历每一个簇，取第i个簇中对应所有样本的轮廓系数，并进行排序
18 s_values = sample_silhouette_values[cluster_labels == i]
19 s_values.sort()
20 size_cluster_i = s_values.shape[0] # 得到第i个簇的样本数量
21 y_upper = y_lower + size_cluster_i # 图中每个簇在y轴上的宽度
22 ax1.fill_betweenx(y=np.arange(y_lower, y_upper), x1=0, x2=s_values, alpha=0.7)
23 ax1.text(-0.05, y_lower + 0.5 * size_cluster_i, str(i))# 在y轴右侧标记每个簇的序号
24 y_lower = y_upper + 10
25 fm.fontManager.addfont('../data/SimHei.ttf')
26 plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来 正常显示中文标签
27 ax1.set_title("不同簇对应的轮廓系数图", fontsize=12)
28 ax1.set_xlabel("轮廓系数", fontsize=12)
29 ax1.set_ylabel("聚类簇序号", fontsize=12)
30 ax1.axvline(x=silhouette_avg, color="red", linestyle="--") # 以x=silhouette_avg 画一条平行于y轴的线
31 ax1.set_yticks([]) # 去掉y轴的刻度
32 ax1.set_xticks([-0.1, 0, 0.2, 0.4, 0.6, 0.8, 1]) # 设置x轴的刻度
33 colors = cm.nipy_spectral(cluster_labels.astype(float) / n_clusters) # 定义每个簇的颜色
34 ax2.scatter(X[:, 0], X[:, 1], marker=".", s=30, lw=0, alpha=0.7, c=colors, edgecolor="k")
35 ax2.scatter(centers[:, 0],centers[:, 1], marker="o", c="white", alpha=1, s=200, edgecolor="k", )
36 for i, c in enumerate(centers): # 给每个簇标记序号
37 ax2.scatter(c[0], c[1], marker="$%d$" % i, alpha=1, s=50, edgecolor="k")
38 ax2.set_title("数据集可视化结果", fontsize=12)
39 ax2.set_xlabel("特征维度 1", fontsize=12)
40 ax2.set_ylabel("特征维度 2", fontsize=12)
41 plt.suptitle("轮廓分析法 K = %d"% n_clusters, fontsize=12, fontweight="bold", )
42 plt.tight_layout()
43 plt.show()在上述代码中，第2~4行是制作数据集并枚举可能的K值。第6~9行是创建一个1行2列的画布并设置其大小、范围和间距等。第10~12行是根据不同的K值对数据集进行聚类并预测得到对应标签和簇中心。第13~15行是计算聚类结果的整体轮廓系数和每个样本对应的轮廓系数。第17~19行是遍历每一个簇，取第i个簇中对应所有样本的轮廓系数，并进行排序。第20~21行是得到第i个簇的样本数量并计算每个簇在$y$轴上的宽度。第22行是填充在$y$轴[y_lower, y_upper]区间内，曲线x1=0和曲线x2=s_values之间的区域。第23行是在$y$轴右侧标记每个簇的序号。第25~26行用来指定字体以显示中文。第30行是绘制图3中垂直的红色虚线；第33~35行是绘制图3右侧的样本点和聚类结果。
进一步，根据类似方法，我们还可以得到图11-20中digits数据集对应的轮廓分析结果，如图11-22所示。
如图11-22所示，这里把轮廓分析法的结果和肘部法的分析结果都放到了一张图里以便于观察。从图示的结果可以得出：
① 当K=6时应该排除；因为根据轮廓分析法来看，其存在簇中样本的轮廓系数值小于平均轮廓系数的情况且各个簇中的样本数量差异过大，同时根据肘部法来看从K=5到K=6时簇内距离的降幅并不如前面K=2到K=3或K=3到K=4明显。
② 当K=5时应该排除；因为根据轮廓分析法的结果来看，第4个簇的样本数量明显过少。
③ 当K=2时应该排除；因为根据肘部法来看从K=2到K=3和K=3到K=4时簇内距离产生了最大的降幅，且根据轮廓分析法可知此时两个簇中的样本数量差异较大。
④ 当K=3以及K=4时，这两种情况稍微容易混淆一点；因为根据肘部法来看K=3到K=4的降幅和K=4到K=5的降幅差别并不大，所以实际的K值既可能是3也可能是4；同时根据轮廓分析法前面2个准则来看，虽然当K=4时各个簇中的样本数量相较于K=3时更加均匀，但是K=3时各个簇中的样本数量差距也并不大，因此也不易直接排除；但是根据轮廓分析法的第3个准则来看，当K=4时的平均轮廓系数要大于K=3时的轮廓系数，因此可以排除K=3时的情况。
总的来说不管是肘部法也好还是轮廓系数法也好，都是建立在实际处理经验上的方法论，因此依旧存在两者都不能对K值进行判断的情况。
11.9.5 小结#
在这节内容中，我们首先详细介绍了为什么不能直接用内部指标加交叉验证的组合来选择最佳K值的原因；然后详细介绍了肘部法的原理及其可视化方法；接着详细介绍了轮廓分析法的原理、代码实现和经验性准则；最后通过一个实际的示例来详细介绍了肘部法和轮廓法的混合使用情况。这里需要明白的是，这两种方法都是建立在实际经验上的分析过程，因此并不能够百分之百的肯定最后分析得到的K值就是数据集真实对应的K值，只是说分析得到的K值存在的可能性很大。
引用#
[1]https://towardsdatascience.com/silhouette-method-better-than-elbow-method-to-find-optimal-clusters-378d62ff6891
[2]https://scikit-learn.org/stable/auto_examples/cluster/plot_kmeans_silhouette_analysis.html
[3]https://medium.com/@cmukesh8688/silhouette-analysis-in-k-means-clustering-cefa9a7ad111