K-Means算法之K值的选择
K-Means是一个超级简单的聚类方法,说他简单,主要原因是使用它时只需设置一个K值(设置需要将数据聚成几类)。但问题是,有时候我们拿到的数据根本不知道要分为几类,对于二维的数据,我们还能通过肉眼观察法进行确定,超过二维的数据怎么办?今天就一起来学习下。
拍脑袋法
一个非常快速的,拍脑袋的方法是将样本量除以2再平方出来的值作为K值,具体公式为:
肘部法则(Elbow Method)
Elbow Method :Elbow意思是手肘,如下图左所示,此种方法适用于 K 值相对较小的情况,当选择的k值小于真正的时,k每增加1,cost值就会大幅的减小;当选择的k值大于真正的K时, k每增加1,cost值的变化就不会那么明显。这样,正确的k值就会在这个转折点,类似elbow的地方。 如下图:
通过画K与cost function的关系曲线图,如左图所示,肘部的值(cost function开始时下降很快,在肘部开始平缓了)做为K值,K=3。并不是所有的问题都可以通过画肘部图来解决,有的问题如右边的那个图,肘点位置不明显(肘点可以是3,4,5),这时就无法确定K值了。故肘部图是可以尝试的一种方法,但是并不是对所有的问题都能画出如左边那么好的图来确定K值。
Elbow Method公式:
Python实现:
# clustering dataset # determine k using elbow method from sklearn.cluster import KMeans from scipy.spatial.distance import cdist import numpy as np import matplotlib.pyplot as plt x1 = np.array([3, 1, 1, 2, 1, 6, 6, 6, 5, 6, 7, 8, 9, 8, 9, 9, 8]) x2 = np.array([5, 4, 5, 6, 5, 8, 6, 7, 6, 7, 1, 2, 1, 2, 3, 2, 3]) plt.plot() plt.xlim([0, 10]) plt.ylim([0, 10]) plt.title('Dataset') plt.scatter(x1, x2) plt.show() # create new plot and data plt.plot() X = np.array(list(zip(x1, x2))).reshape(len(x1), 2) colors = ['b', 'g', 'r'] markers = ['o', 'v', 's'] # k means determine k distortions = [] K = range(1, 10) for k in K: kmeanModel = KMeans(n_clusters=k).fit(X) kmeanModel.fit(X) distortions.append(sum(np.min(cdist(X, kmeanModel.cluster_centers_, 'euclidean'), axis=1)) / X.shape[0]) # Plot the elbow plt.plot(K, distortions, 'bx-') plt.xlabel('k') plt.ylabel('Distortion') plt.title('The Elbow Method showing the optimal k') plt.show()
间隔统计量 Gap Statistic
根据肘部法则选择最合适的K值有事并不是那么清晰,因此斯坦福大学的Robert等教授提出了Gap Statistic方法。
Gap Statistic的基本思路是:引入参考的测值,这个参考值可以有Monte Carlo采样的方法获得。
B是sampling的次数。为了修正MC带来的误差,我们计算sk也即标准差来矫正Gap Statistic。
选择满足的最小的k作为最优的聚类个数。下图阐释了Gap Statistic的过程。
Python实现:
import scipy from scipy.spatial.distance import euclidean from sklearn.cluster import KMeans as k_means dst = euclidean k_means_args_dict = { 'n_clusters': 0, # drastically saves convergence time 'init': 'k-means++', 'max_iter': 100, 'n_init': 1, 'verbose': False, # 'n_jobs':8 } def gap(data, refs=None, nrefs=20, ks=range(1, 11)): """ I: NumPy array, reference matrix, number of reference boxes, number of clusters to test O: Gaps NumPy array, Ks input list Give the list of k-values for which you want to compute the statistic in ks. By Gap Statistic from Tibshirani, Walther. """ shape = data.shape if not refs: tops = data.max(axis=0) bottoms = data.min(axis=0) dists = scipy.matrix(scipy.diag(tops - bottoms)) rands = scipy.random.random_sample(size=(shape[0], shape[1], nrefs)) for i in range(nrefs): rands[:, :, i] = rands[:, :, i] * dists + bottoms else: rands = refs gaps = scipy.zeros((len(ks),)) for (i, k) in enumerate(ks): k_means_args_dict['n_clusters'] = k kmeans = k_means(**k_means_args_dict) kmeans.fit(data) (cluster_centers, point_labels) = kmeans.cluster_centers_, kmeans.labels_ disp = sum( [dst(data[current_row_index, :], cluster_centers[point_labels[current_row_index], :]) for current_row_index in range(shape[0])]) refdisps = scipy.zeros((rands.shape[2],)) for j in range(rands.shape[2]): kmeans = k_means(**k_means_args_dict) kmeans.fit(rands[:, :, j]) (cluster_centers, point_labels) = kmeans.cluster_centers_, kmeans.labels_ refdisps[j] = sum( [dst(rands[current_row_index, :, j], cluster_centers[point_labels[current_row_index], :]) for current_row_index in range(shape[0])]) # let k be the index of the array 'gaps' gaps[i] = scipy.mean(scipy.log(refdisps)) - scipy.log(disp) return ks, gaps
轮廓系数(Silhouette Coefficient)
Silhouette method 会衡量对象和所属簇之间的相似度——即内聚性(cohesion)。当把它与其他簇做比较,就称为分离性(separation)。该对比通过 silhouette 值来实现,后者在 [-1, 1] 范围内。Silhouette 值接近 1,说明对象与所属簇之间有密切联系;反之则接近 -1。若某模型中的一个数据簇,生成的基本是比较高的 silhouette 值,说明该模型是合适、可接受的。
方法:
1)计算样本i到同簇其他样本的平均距离a(i)。a(i)越小,说明样本i越应该被聚类到该簇。将a(i)称为样本i的簇内不相似度。簇C中所有样本的a(i)均值称为簇C的簇不相似度。
2)计算样本i到其他某簇C(j)的所有样本的平均距离b(ij),称为样本i与簇C(j)的不相似度。定义为样本i的簇间不相似度:b(i) =min{bi1, bi2, …, bik},b(i)越大,说明样本i越不属于其他簇。
3)根据样本i的簇内不相似度a i 和簇间不相似度b i ,定义样本i的轮廓系数:
4)判断:
s(i)接近1,则说明样本i聚类合理
s(i)接近-1,则说明样本i更应该分类到另外的簇
若s(i) 近似为0,则说明样本i在两个簇的边界上
所有样本的s(i )的均值称为聚类结果的轮廓系数,是该聚类是否合理、有效的度量。但是,其缺陷是计算复杂度为O(n^2),需要计算距离矩阵,那么当数据量上到百万,甚至千万级别时,计算开销会非常巨大。
Python实现:
import numpy as np from sklearn.cluster import KMeans from sklearn import metrics import matplotlib.pyplot as plt plt.figure(figsize=(8, 10)) plt.subplot(3, 2, 1) x1 = np.array([1, 2, 3, 1, 5, 6, 5, 5, 6, 7, 8, 9, 7, 9]) x2 = np.array([1, 3, 2, 2, 8, 6, 7, 6, 7, 1, 2, 1, 1, 3]) X = np.array(list(zip(x1, x2))).reshape(len(x1), 2) plt.xlim([0, 10]) plt.ylim([0, 10]) plt.title('Sample') plt.scatter(x1, x2) colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'b'] markers = ['o', 's', 'D', 'v', '^', 'p', '*', '+'] tests = [2, 3, 4, 5, 8] subplot_counter = 1 for t in tests: subplot_counter += 1 plt.subplot(3, 2, subplot_counter) kmeans_model = KMeans(n_clusters=t).fit(X) for i, l in enumerate(kmeans_model.labels_): plt.plot(x1[i], x2[i], color=colors[l], marker=markers[l],ls='None') plt.xlim([0, 10]) plt.ylim([0, 10]) plt.title('K = %s, Silhouette method = %.03f' % (t, metrics.silhouette_score(X, kmeans_model.labels_,metric='euclidean'))) plt.show()
Canopy算法
肘部法则(Elbow Method)和轮廓系数(Silhouette Coefficient)来对k值进行最终的确定,但是这些方法都是属于“事后”判断的,而Canopy算法的作用就在于它是通过事先粗聚类的方式,为k-means算法确定初始聚类中心个数和聚类中心点。
与传统的聚类算法(比如K-Means)不同,Canopy聚类最大的特点是不需要事先指定k值(即clustering的个数),因此具有很大的实际应用价值。与其他聚类算法相比,Canopy聚类虽然精度较低,但其在速度上有很大优势,因此可以使用Canopy聚类先对数据进行“粗”聚类,得到k值,以及大致的k个中心点,再使用K-Means进行进一步“细”聚类。所以Canopy+K-Means这种形式聚类算法聚类效果良好。
Canopy算法解析:
原始数据集合List按照一定的规则进行排序(这个规则是任意的,但是一旦确定就不再更改),初始距离阈值为T1、T2,且T1>T2(T1、T2的设定可以根据用户的需要,或者使用交叉验证获得)。
在List中随机挑选一个数据向量A,使用一个粗糙距离计算方式计算A与List中其他样本数据向量之间的距离d。
根据第2步中的距离d,把d小于T1的样本数据向量划到一个canopy中,同时把d小于T2的样本数据向量从候选中心向量名单(这里可以理解为就是List)中移除。
重复第2、3步,直到候选中心向量名单为空,即List为空,算法结束。
算法原理比较简单,就是对数据进行不断遍历,T2<dis<T1的可以作为中心名单,dis<T2的认为与canopy太近了,以后不会作为中心点,从list中删除,这样的话一个点可能属于多个canopy。
canopy效果图如下:
Canopy算法优势:
Kmeans对噪声抗干扰较弱,通过Canopy对比较小的NumPoint的Cluster直接去掉 有利于抗干扰。
Canopy选择出来的每个Canopy的centerPoint作为Kmeans比较科学。
只是针对每个Canopy的内容做Kmeans聚类,减少相似计算的数量。
Canopy算法缺点:
算法中 T1、T2(T2 < T1) 的确定问题
Python实现:
# -*- coding: utf-8 -*- import math import random import numpy as np import matplotlib.pyplot as plt class Canopy: def __init__(self, dataset): self.dataset = dataset self.t1 = 0 self.t2 = 0 # 设置初始阈值 def setThreshold(self, t1, t2): if t1 > t2: self.t1 = t1 self.t2 = t2 else: print('t1 needs to be larger than t2!') # 使用欧式距离进行距离的计算 def euclideanDistance(self, vec1, vec2): return math.sqrt(((vec1 - vec2)**2).sum()) # 根据当前dataset的长度随机选择一个下标 def getRandIndex(self): return random.randint(0, len(self.dataset) - 1) def clustering(self): if self.t1 == 0: print('Please set the threshold.') else: canopies = [] # 用于存放最终归类结果 while len(self.dataset) != 0: rand_index = self.getRandIndex() current_center = self.dataset[rand_index] # 随机获取一个中心点,定为P点 current_center_list = [] # 初始化P点的canopy类容器 delete_list = [] # 初始化P点的删除容器 self.dataset = np.delete( self.dataset, rand_index, 0) # 删除随机选择的中心点P for datum_j in range(len(self.dataset)): datum = self.dataset[datum_j] distance = self.euclideanDistance( current_center, datum) # 计算选取的中心点P到每个点之间的距离 if distance < self.t1: # 若距离小于t1,则将点归入P点的canopy类 current_center_list.append(datum) if distance < self.t2: delete_list.append(datum_j) # 若小于t2则归入删除容器 # 根据删除容器的下标,将元素从数据集中删除 self.dataset = np.delete(self.dataset, delete_list, 0) canopies.append((current_center, current_center_list)) return canopies def showCanopy(canopies, dataset, t1, t2): fig = plt.figure() sc = fig.add_subplot(111) colors = ['brown', 'green', 'blue', 'y', 'r', 'tan', 'dodgerblue', 'deeppink', 'orangered', 'peru', 'blue', 'y', 'r', 'gold', 'dimgray', 'darkorange', 'peru', 'blue', 'y', 'r', 'cyan', 'tan', 'orchid', 'peru', 'blue', 'y', 'r', 'sienna'] markers = ['*', 'h', 'H', '+', 'o', '1', '2', '3', ',', 'v', 'H', '+', '1', '2', '^', '<', '>', '.', '4', 'H', '+', '1', '2', 's', 'p', 'x', 'D', 'd', '|', '_'] for i in range(len(canopies)): canopy = canopies[i] center = canopy[0] components = canopy[1] sc.plot(center[0], center[1], marker=markers[i], color=colors[i], markersize=10) t1_circle = plt.Circle( xy=(center[0], center[1]), radius=t1, color='dodgerblue', fill=False) t2_circle = plt.Circle( xy=(center[0], center[1]), radius=t2, color='skyblue', alpha=0.2) sc.add_artist(t1_circle) sc.add_artist(t2_circle) for component in components: sc.plot(component[0], component[1], marker=markers[i], color=colors[i], markersize=1.5) maxvalue = np.amax(dataset) minvalue = np.amin(dataset) plt.xlim(minvalue - t1, maxvalue + t1) plt.ylim(minvalue - t1, maxvalue + t1) plt.show() if __name__ == "__main__": dataset = np.random.rand(500, 2) # 随机生成500个二维[0,1)平面点 t1 = 0.6 t2 = 0.4 gc = Canopy(dataset) gc.setThreshold(t1, t2) canopies = gc.clustering() print('Get %s initial centers.' % len(canopies)) showCanopy(canopies, dataset, t1, t2)
参考资料:
https://en.wikipedia.org/wiki/Determining_the_number_of_clusters_in_a_data_set
https://www.quora.com/How-can-we-choose-a-good-K-for-K-means-clustering
https://stackoverflow.com/questions/1793532/how-do-i-determine-k-when-using-k-means-clustering
扫一扫订阅我的微信号:IT技术博客大学习
- 作者:标点符 来源: 标点符
- 标签: K-Means
- 发布时间:2018-06-26 12:36:07
- [55] Oracle MTS模式下 进程地址与会话信
- [55] IOS安全–浅谈关于IOS加固的几种方法
- [54] 如何拿下简短的域名
- [53] android 开发入门
- [52] Go Reflect 性能
- [52] 图书馆的世界纪录
- [49] 【社会化设计】自我(self)部分――欢迎区
- [49] 读书笔记-壹百度:百度十年千倍的29条法则
- [38] 程序员技术练级攻略
- [32] 视觉调整-设计师 vs. 逻辑