K-means聚类选择最优k值

  最近做了一个数据挖掘的项目,挖掘过程中用到了K-means聚类方法,但是由于根据行业经验确定的聚类数过多并且并不一定是我们获取到数据的真实聚类数,所以,我们希望能从数据自身出发去确定真实的聚类数,也就是对数据而言的最佳聚类数。为此,我查阅了大量资料和博客资源,总结出主流的确定聚类数k的方法有以下两类。

1.手肘法
1.1 理论
手肘法的核心指标是SSE(sum of the squared errors,误差平方和),

                                         

 

其中,Ci是第i个簇,p是Ci中的样本点,mi是Ci的质心(Ci中所有样本的均值),SSE是所有样本的聚类误差,代表了聚类效果的好坏。

       手肘法的核心思想是:随着聚类数k的增大,样本划分会更加精细,每个簇的聚合程度会逐渐提高,那么误差平方和SSE自然会逐渐变小。并且,当k小于真实聚类数时,由于k的增大会大幅增加每个簇的聚合程度,故SSE的下降幅度会很大,而当k到达真实聚类数时,再增加k所得到的聚合程度回报会迅速变小,所以SSE的下降幅度会骤减,然后随着k值的继续增大而趋于平缓,也就是说SSE和k的关系图是一个手肘的形状,而这个肘部对应的k值就是数据的真实聚类数。当然,这也是该方法被称为手肘法的原因。

1.2 实践
  我们对预处理后数据.csv 中的数据利用手肘法选取最佳聚类数k。具体做法是让k从1开始取值直到取到你认为合适的上限(一般来说这个上限不会太大,这里我们选取上限为8),对每一个k值进行聚类并且记下对于的SSE,然后画出k和SSE的关系图(毫无疑问是手肘形),最后选取肘部对应的k作为我们的最佳聚类数。python实现如下:

import pandas as pd
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt
 
df_features = pd.read_csv(r'C:\预处理后数据.csv',encoding='gbk') # 读入数据
'利用SSE选择k'
SSE = []  # 存放每次结果的误差平方和
for k in range(1,9):
    estimator = KMeans(n_clusters=k)  # 构造聚类器
    estimator.fit(df_features[['R','F','M']])
    SSE.append(estimator.inertia_)
X = range(1,9)
plt.xlabel('k')
plt.ylabel('SSE')
plt.plot(X,SSE,'o-')
plt.show()
#有关KMeans及相关用法:
sklearn KMeans用法

画出的k与SSE的关系图如下:

                                         

显然,肘部对于的k值为4,故对于这个数据集的聚类而言,最佳聚类数应该选4

 

2. 轮廓系数法
2.1 理论
该方法的核心指标是轮廓系数(Silhouette Coefficient),某个样本点Xi的轮廓系数定义如下:

                                                           

 

其中,a是Xi与同簇的其他样本的平均距离,称为凝聚度,b是Xi与最近簇中所有样本的平均距离,称为分离度。而最近簇的定义是

                                                    

 

其中p是某个簇Ck中的样本。事实上,简单点讲,就是用Xi到某个簇所有样本平均距离作为衡量该点到该簇的距离后,选择离Xi最近的一个簇作为最近簇。

       求出所有样本的轮廓系数后再求平均值就得到了平均轮廓系数。平均轮廓系数的取值范围为[-1,1],且簇内样本的距离越近,簇间样本距离越远,平均轮廓系数越大,聚类效果越好。那么,很自然地,平均轮廓系数最大的k便是最佳聚类数。

2.2 实践
我们同样使用2.1中的数据集,同样考虑k等于1到8的情况,对于每个k值进行聚类并且求出相应的轮廓系数,然后做出k和轮廓系数的关系图,选取轮廓系数取值最大的k作为我们最佳聚类系数,python实现如下:

import pandas as pd
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
import matplotlib.pyplot as plt
 
df_features = pd.read_csv(r'C:\Users\61087\Desktop\项目\爬虫数据\预处理后数据.csv',encoding='gbk')
Scores = []  # 存放轮廓系数
for k in range(2,9):
    estimator = KMeans(n_clusters=k)  # 构造聚类器
    estimator.fit(df_features[['R','F','M']])
    Scores.append(silhouette_score(df_features[['R','F','M']],estimator.labels_,metric='euclidean'))
X = range(2,9)
plt.xlabel('k')
plt.ylabel('轮廓系数')
plt.plot(X,Scores,'o-')
plt.show()

得到聚类数k与轮廓系数的关系图:

                                

  可以看到,轮廓系数最大的k值是2,这表示我们的最佳聚类数为2。但是,值得注意的是,从k和SSE的手肘图可以看出,当k取2时,SSE还非常大,所以这是一个不太合理的聚类数,我们退而求其次,考虑轮廓系数第二大的k值4,这时候SSE已经处于一个较低的水平,因此最佳聚类系数应该取4而不是2。

       但是,讲道理,k=2时轮廓系数最大,聚类效果应该非常好,那为什么SSE会这么大呢?在我看来,原因在于轮廓系数考虑了分离度b,也就是样本与最近簇中所有样本的平均距离。为什么这么说,因为从定义上看,轮廓系数大,不一定是凝聚度a(样本与同簇的其他样本的平均距离)小,而可能是b和a都很大的情况下b相对a大得多,这么一来,a是有可能取得比较大的。a一大,样本与同簇的其他样本的平均距离就大,簇的紧凑程度就弱,那么簇内样本离质心的距离也大,从而导致SSE较大。所以,虽然轮廓系数引入了分离度b而限制了聚类划分的程度,但是同样会引来最优结果的SSE比较大的问题,这一点也是值得注意的。

总结
  从以上两个例子可以看出,轮廓系数法确定出的最优k值不一定是最优的,有时候还需要根据SSE去辅助选取,这样一来相对轮廓系数法就显得有点累赘。因此,如果没有特殊情况的话,我还是建议首先考虑用手肘法。

 

3.Calinski-Harabasz Index(CH指标)

  不像监督学习的分类问题和回归问题,我们的无监督聚类没有样本输出,也就没有比较直接的聚类评估方法。但是我们可以从簇内的稠密程度和簇间的离散程度来评估聚类的效果。常见的方法有轮廓系数Silhouette Coefficient和Calinski-Harabasz Index。个人比较喜欢Calinski-Harabasz Index,这个计算简单直接,得到的Calinski-Harabasz分数值ss越大则聚类效果越好。

    Calinski-Harabasz分数值ss的数学计算公式是:

           

   也就是说,类别内部数据的协方差越小越好,类别之间的协方差越大越好,这样的Calinski-Harabasz分数会高。在scikit-learn中, Calinski-Harabasz Index对应的方法是metrics.calinski_harabaz_score.

   K-Means应用实例

   下面用一个实例来讲解用KMeans类和MiniBatchKMeans类来聚类。我们观察在不同的k值下Calinski-Harabasz分数。完整的代码参见我的github: https://github.com/ljpzzz/machinelearning/blob/master/classic-machine-learning/kmeans_cluster.ipynb首先我们随机创建一些二维数据作为训练集,选择二维特征数据,主要是方便可视化。代码如下:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from sklearn.datasets.samples_generator import make_blobs
# X为样本特征,Y为样本簇类别, 共1000个样本,每个样本2个特征,共4个簇,簇中心在[-1,-1], [0,0],[1,1], [2,2], 簇方差分别为[0.4, 0.2, 0.2]
X, y = make_blobs(n_samples=1000, n_features=2, centers=[[-1,-1], [0,0], [1,1], [2,2]], cluster_std=[0.4, 0.2, 0.2, 0.2], 
                  random_state =9)
plt.scatter(X[:, 0], X[:, 1], marker='o')
plt.show()

从输出图可以我们看看我们创建的数据如下:

                                          

现在我们来用K-Means聚类方法来做聚类,首先选择k=2,代码如下:

from sklearn.cluster import KMeans
y_pred = KMeans(n_clusters=2, random_state=9).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()

k=2聚类的效果图输出如下:

                                        

现在我们来看看我们用Calinski-Harabasz Index评估的聚类分数:

from sklearn import metrics
metrics.calinski_harabaz_score(X, y_pred) 

    输出如下:

3116.1706763322227

    现在k=3来看看聚类效果,代码如下:

from sklearn.cluster import KMeans
y_pred = KMeans(n_clusters=3, random_state=9).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()  

k=3的聚类的效果图输出如下:

                                      

现在我们来看看我们用Calinski-Harabaz Index评估的k=3时候聚类分数:

metrics.calinski_harabaz_score(X, y_pred)  

    输出如下:

2931.625030199556

    可见此时k=3的聚类分数比k=2还差。

    现在我们看看k=4时候的聚类效果:

from sklearn.cluster import KMeans
y_pred = KMeans(n_clusters=4, random_state=9).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()

k=4的聚类的效果图输出如下:

                                      

现在我们来看看我们用Calinski-Harabasz Index评估的k=4时候聚类分数:

metrics.calinski_harabaz_score(X, y_pred)  

    输出如下:

5924.050613480169

    可见k=4的聚类分数比k=2和k=3都要高,这也符合我们的预期,我们的随机数据集也就是4个簇。当特征维度大于2,我们无法直接可视化聚类效果来肉眼观察时,用Calinski-Harabaz Index评估是一个很实用的方法。

    现在我们再看看用MiniBatchKMeans的效果,我们将batch size设置为200. 由于我们的4个簇都是凸的,所以其实batch size的值只要不是非常的小,对聚类的效果影响不大。

for index, k in enumerate((2,3,4,5)):
    plt.subplot(2,2,index+1)
    y_pred = MiniBatchKMeans(n_clusters=k, batch_size = 200, random_state=9).fit_predict(X)
    score= metrics.calinski_harabaz_score(X, y_pred)  
    plt.scatter(X[:, 0], X[:, 1], c=y_pred)
    plt.text(.99, .01, ('k=%d, score: %.2f' % (k,score)),
                 transform=plt.gca().transAxes, size=10,
                 horizontalalignment='right')
plt.show()

对于k=2,3,4,5对应的输出图为:

                                

可见使用MiniBatchKMeans的聚类效果也不错,当然由于使用Mini Batch的原因,同样是k=4最优,KMeans类的Calinski-Harabasz Index分数为5924.05,而MiniBatchKMeans的分数稍微低一些,为5921.45。这个差异损耗并不大。

---------------------
版权声明:本文前两部分为CSDN博主「VultureEye」的原创文章,遵循CC 4.0 by-sa版权协议,转载请附上原文出处链接及本声明。
原文链接:https://blog.csdn.net/qq_15738501/article/details/79036255

第三部分转自:https://www.cnblogs.com/pinard/p/6169370.html

 

posted @ 2019-08-12 15:07  AI_ON  阅读(2524)  评论(0)    收藏  举报