Kmeans聚类时K值选择的方法

1.K-means算法

（1）简单介绍

聚类属于非监督学习，K均值聚类是最基础常用的聚类算法。它的基本思想是，通过迭代寻找K个簇（Cluster）的一种划分方案，使得聚类结果对应的损失函数最小。其中，损失函数可以定义为各个样本距离所属簇中心点的误差平方和：

其中代表第个样本，是所属的簇，代表簇对应的中心点，是样本总数。

（2）具体步骤

KMeans的核心目标是将给定的数据集划分成K个簇（K是超参），并给出每个样本数据对应的中心点。具体步骤非常简单，可以分为4步：

（1）数据预处理。主要是标准化、异常点过滤。

（2）随机选取K个中心，记为

（3）定义损失函数：

（4）令t=0,1,2,… 为迭代步数，重复如下过程直到收敛：

（4.1）对于每一个样本，将其分配到距离最近的中心

（4.2）对于每一个类中心k，重新计算该类的中心

KMeans最核心的部分就是先固定中心点，调整每个样本所属的类别来减少 J；再固定每个样本的类别，调整中心点继续减小J 。两个过程交替循环， J单调递减直到最（极）小值，中心点和样本划分的类别同时收敛。

2.K值的选取方法

（1）手肘法

手肘法将簇间误差平方和看成是类簇数量k的函数。随着k的增加，每个类簇内的离散程度越小，总距离平方和也就在不断减小，并且减小的程度越来越不明显。极限情况是当k=N时，每个类簇只有一个点，这时总的误差平方和为0。手肘法认为我们应该选择这样的k：当k继续增大时，总误差平方和减少的趋势不再明显，也就是“拐点”处。具体过程如下：

选择一个聚类算法（例如K-means），计算不同k时的聚类结果，例如k可以取为0～10。
对每个k，计算总的簇间距离平方和。
画出总簇间距离平方和随k值增加的变化趋势。
图中弯曲的“拐点”处对应的k就是最合适的类簇数量

# 手肘法调研了一下基本是画出图片以后，采取目测的方式选择合适的K,
# 这里我自己写了一个获取K的方法，好像有点不准
import matplotlib.pyplot as plt
from sklearn.cluster import KMeansdef get_k_value(distortions, start_class_num = 1):"""通过手肘法计算最优的k值Args:border_entity_info: Returns:k: 最优的k值"""k = 0for i in range(len(distortions) - 1):if distortions[i] - distortions[i+1] < 1:k = i + start_class_numbreakreturn kdef elbow_method_K(data, range_K, pro_num):K = range(1, range_K + 1)meandistortions = []for k in K:kmeans = KMeans(n_clusters=k)kmeans.fit(data)meandistortions.append(kmeans.inertia_)best_k = get_k_value(meandistortions)plt.plot(K, meandistortions, 'bx-')plt.xlabel('k')plt.ylabel('Average Dispersion')plt.title('Selecting k with the Elbow Method')plt.savefig(f'/Users/cecilia/Desktop/K_图片/{pro_num}_elbow_K.jpg')plt.cla()return best_k
# 这个函数是我自己使用的时候封装的
# data是需要进行聚类的数据，可以是多维矩阵
# range_K是类别的可选择范围
# pro_num是名称，没有实际意义是为了将图片保存下来，不想保存图片可以直接使用plt.show()

（2）Gap Statistic

是斯坦福大学的三位教授在2001年的一篇论文中(R. Tibshirani, G. Walther, and T. Hastie, 2001)提出来的，可用于任何的聚类方法。Gap Statistic的主要思想是比较不同k时原始数据的簇内偏差总和与数据在均匀分布推断下的簇内偏差总和。使Gap Statistic这个统计量达到最大值意味着此时的聚类结果结构与随机均匀分布产生的数据的聚类结果差别最大，此时的k就是最优的k。算法如下：

将原始的观测数据进行聚类，k=0,…, k_max，计算不同k值对应的簇内偏离和W_k。
通过随机的均匀分布产生B个推断数据，对这些推断数据进行聚类，k=0,…, k_max。计算不同k值对应的在B个推断数据上的平均簇内偏离和W_kb。
计算gap statistic：W_k与W_kb的log偏差Gap(k)。同时计算这个偏差的标准差sd_k，然后令s_k = sprt(1+1/B*sd_k)。
选择一个最小的k，这样的k满足Gap(k) > Gap(k+1) - s_k+1。
流行的做法是直接选择最大的Gap(k)所对应的k作为最优k，也就是忽略上述的第四步。需要注意的是当B=500时，W_kb是非常精确的，在下一次迭代中基本保持不变。
注意⚠️：使用这个需要安装一个库，具体信息可以看Gap Statistic

from gap_statistic import OptimalKef gap_statistic_K(data, range_K, pro_num):K = np.arange(1, range_K)optimalK = OptimalK(n_jobs=1, parallel_backend='joblib')n_clusters = optimalK(data, cluster_array=K)# Gap values plotplt.plot(optimalK.gap_df.n_clusters, optimalK.gap_df.gap_value, linewidth=3)plt.scatter(optimalK.gap_df[optimalK.gap_df.n_clusters == n_clusters].n_clusters,optimalK.gap_df[optimalK.gap_df.n_clusters == n_clusters].gap_value, s=250, c='r')plt.grid(True)plt.xlabel('Cluster Count')plt.ylabel('Gap Value')plt.title('Gap Values by Cluster Count')plt.savefig(f'/Users/cecilia/Desktop/K_图片/{pro_num}_gap_values_K.jpg')plt.cla()# plt.show()# # diff plot# plt.plot(optimalK.gap_df.n_clusters, optimalK.gap_df["diff"], linewidth=3)# plt.grid(True)# plt.xlabel("Cluster Count")# plt.ylabel("Diff Value")# plt.title("Diff Values by Cluster Count")# # plt.show()# plt.savefig(f'/Users/cecilia/Desktop/K_图片/{pro_num}_diff_2.jpg')# plt.cla()# Gap* plot# max_ix = optimalK.gap_df[optimalK.gap_df["gap*"] == optimalK.gap_df["gap*"].max()].index[0]# plt.plot(optimalK.gap_df.n_clusters, optimalK.gap_df["gap*"], linewidth=3)# plt.scatter(#     optimalK.gap_df.loc[max_ix]["n_clusters"],#     optimalK.gap_df.loc[max_ix]["gap*"],#     s=250,#     c="r",# )# plt.grid(True)# plt.xlabel("Cluster Count")# plt.ylabel("Gap* Value")# plt.title("Gap* Values by Cluster Count")# plt.savefig(f'/Users/cecilia/Desktop/K_图片/{pro_num}_Gap*_3.jpg')# plt.cla()# plt.show()# # diff* plot# plt.plot(optimalK.gap_df.n_clusters, optimalK.gap_df["diff*"], linewidth=3)# plt.grid(True)# plt.xlabel("Cluster Count")# plt.ylabel("Diff* Value")# plt.title("Diff* Values by Cluster Count")# plt.savefig(f'/Users/cecilia/Desktop/K_图片/{pro_num}_diff*_4.jpg')# plt.cla()# plt.show()return n_clusters

（3）平均轮廓系数法

平均轮廓系数方法衡量了聚类结果的质量，即衡量每个点被放到当前类簇有多合适，平均轮廓系数很高意味着聚类的结果很好。这种方法计算不同k值下，所有点的平均轮廓系数，能够使平均轮廓系数最大的k就是最优的类簇数量（Kaufman and Rousseeuw 1990）。
具体的过程与手肘法是相似的：

选择一个聚类算法（例如K-means），计算不同k时的聚类结果，例如k可以取为0～10。
对于每个k，计算所有观测点的平均轮廓系数。
画出这个指标随着k变化的曲线。
曲线中最高点对应的k就是最优聚类数量。

from sklearn.metrics import silhouette_score
from sklearn.cluster import KMeansdef get_silhouette_K(data, range_K, pro_num):K = range(2, range_K)Scores = [] for k in K:kmeans = KMeans(n_clusters=k)kmeans.fit(data)Scores.append(silhouette_score(data, kmeans.labels_, metric='euclidean'))max_idx = Scores.index(max(Scores))best_k = K[max_idx]plt.plot(K, Scores, 'bx-')plt.xlabel('k')plt.ylabel('silhouette')plt.title('Selecting k with the silhouette Method')plt.savefig(f'/Users/cecilia/Desktop/K_图片/{pro_num}_silhouette_K.jpg')plt.cla()return best_k
# 这个函数是我自己使用的时候封装的
# data是需要进行聚类的数据，可以是多维矩阵
# range_K是类别的可选择范围
# pro_num是名称，没有实际意义是为了将图片保存下来，不想保存图片可以直接使用plt.show()

注意⚠️：在使用轮廓系数法时，遇到一个问题就是K值的选择必须从2开始，最多只能选择n_samples -1(最大K候选就是样本数量-1)，不然会报错的，具体没有细细研究。