聚类 | KMeans理论与算法实现

01 物以类聚

经过半年的不懈努力，我们已经学习并实践了经典的分类算法和经典的回归算法，下面我们开始学习经典的聚类算法(兴奋~~~)

目前打算对三种聚类算法进行学习和代码实操(俗称“造轮子”)：

1. KMeans
2. Apriori
3. FP-Growth

今天我们学习并实践KMeans聚类算法，分成以下几个部分，跟上节奏燥起来！

1. KMeans算法理论和代码实现
2. 改进，BiKMeans算法理论和代码实现
3. 实例，上车点规划
4. 抉择，如何挑选最佳的聚类簇数？

其中一二三部分参考了Peter Harrington的《机器学习实战》，第四部分是纯手工打造的轮子，还请多多指教。

---

02 KMeans理论和算法实现

聚类是一种无监督学习的方法，所谓“无监督”，就是指参与训练的样本没有标签。

KMeans聚类算法过程如下：
1. 对于一组数据集，随机选取k个点作为质心，将数据集中的点归为离其最近的质心一簇，此时数据集被划分为k个簇；
2. 对这k个簇，重新计算各簇的质心(均值)；
3. 根据新的质心，按照step1继续聚类，然后再根据聚类重新计算各簇质心，直到质心不再改变，分类完成。

说白了，就是不断地聚类、划分的过程。

通过KMeans原理，可以看到几个显而易见的缺点：
1. 簇数量k由用户指定，无法预先知道最佳k值 >>解法：分为几簇，最终由轮廓系数S(i)决定，取轮廓系数最大的分类数(05节剧透)
2. 最终质心可能与初始点选择有关 >> 因此KMeans的结果可能收敛到局部最小值，而不是全局最小值 >> 解法：

• BiKMeans（03节）
• KMeans++（KMeans++ 算法在选择初始质心时并不是随机选择，而是选择尽量相互分离的质心，即，下一个质心点总是离上一个质心点较远的点）

代码实现

def loadDataSet(fileName):
    dataList=[]
    dataMat=[]
    fr=open(fileName)
    for line in fr.readlines():
        curLine=line.strip().split('\t')
        fltLine=list(map(float,curLine))
        dataList.append(fltLine)
        dataMat=mat(dataList)
    return dataMat

def distEclud(vecA,vecB):
    return sqrt(sum(power(vecA-vecB,2))) #欧式距离

#为输入数据集构造k个随机中心，中心位置在各特征最大最小值之间
def randCent(dataSet,k):
    n=shape(dataSet)[1]
    center=mat(zeros((k,n)))
    for j in range(n): #对每个特征
        minJ=min(dataSet[:,j])
        rangeJ=float(max(dataSet[:,j])-minJ)
        center[:,j]=mat(minJ+rangeJ*random.rand(k,1)) #质心第j维坐标在数据集第j维数据之间
    return center

def KMeans(dataSet,k,distMeas=distEclud,createCent=randCent):
    m=shape(dataSet)[0] 
    clusterAssment=mat(zeros((m,2))) #用于记录各样本当前归属于哪个簇以及到该簇质心的欧式距离平方
    center=createCent(dataSet,k)
    clusterChanged=True
    
    while clusterChanged:
        clusterChanged=False     
        #对每个样本，计算样本到各质心的距离，寻找距离最近的质心，将该样本归为该质心所在簇
        for i in range(m): 
            minDist=inf;minIndex=-1
            for j in range(k): #对每个质心,计算到i样本的距离
                distJI=distMeas(center[j,:],dataSet[i,:])
                if distJI<minDist:
                    minDist=distJI;minIndex=j #i样本暂属于j簇，到j簇质心距离为minDist
            if clusterAssment[i,0]!=minIndex:
                clusterChanged=True #若任一样本在本次迭代中改变了簇类，则要进行下一次迭代(即，直到任何样本都不再改变簇类，聚类停止)
            clusterAssment[i,:]=minIndex,minDist**2 #记录样本i的簇类情况
        #print (center)
        #更新质心
        for cent in range(k):
            ptsInClust=dataSet[nonzero(clusterAssment[:,0].A==cent)[0]] #筛选出属于当前簇类的点
            center[cent,:]=mean(ptsInClust,axis=0) #对该簇类各样本的各列求均值，作为新质心
    return center,clusterAssment

用一组数据来测试一下：

dataMat1=loadDataSet(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\testSet.txt')
dataMat2=loadDataSet(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\testSet2.txt')
center_testSet1,clusterAssment_testSet1=KMeans(dataMat1,4)
center_testSet2,clusterAssment_testSet2=KMeans(dataMat2,3)

plt.figure(figsize=(6,6))
plt.scatter(dataMat1[:,0].T.tolist()[0],dataMat1[:,1].T.tolist()[0],c='pink',s=30)
plt.scatter(center_testSet1.T[0].tolist()[0],center_testSet1.T[1].tolist()[0],c='blue',s=50)

结果如下，蓝色点为聚类质心

---

03 BiKMeans理论和算法实现

刚才我们实现了KMeans算法，不过也提到，KMeans有个缺点，就是其聚类的最终质心可能与初始点选择有关，因此KMeans的结果可能收敛到局部最小值，而不是全局最小值。

怎么解决呢？提示：想想决策树是如何分支的？

为了得到全局最优解，BiKMeans算法出现了，它的过程如下（是不是跟决策树分支有点神似？）
1. 将整个数据集看作一个簇，计算初始质心，即所有数据点各特征的均值
2. 遍历各质心，对各质心，将质心所在簇用原始KMeans算法二分，计算二分后整个数据集的SSE(即平方误差和，即簇各点到簇质心距离平方和)，找到二分后整体数据集SSE最小的质心，认为此质心是本次划分的最佳质心，对其进行二分
3. 不断重复step2，直到质心总数=设置的k

说白了，BiKMeans算法过程类似于决策树的分支，通过启发的方法，每次迭代只分裂当前最佳质心，直到簇数量达到k，这样的方法可以保证最终划分得到的质心是全局最优解，而原始KMeans可能会陷入局部最优解。

代码实现

def biKMeans(dataSet,k,distMeas=distEclud):
    m=shape(dataSet)[0]
    clusterAssment=mat(zeros((m,2))) #记录各样本归属和距离平方
    center0=mean(dataSet,axis=0).tolist()[0] #初始质心
    centerList=[center0] #用于记录聚类质心坐标
    for j in range(m):
        clusterAssment[j,1]=distMeas(mat(center0),dataSet[j,:])**2
    while len(centerList)<k:
        lowestSSE=inf #SSE=Sum of Square Error
        #对每个簇，尝试二分，计算二分后整体数据的SSE，若小于lowestSSE，则将簇二分，如此往复，直到分到k个簇为止
        for i in range(len(centerList)):
            ptsInCluster=dataSet[nonzero(clusterAssment[:,0].A==i)[0],:] #默认质心索引就是数据对应簇类
            centerMat,splitClustAss=KMeans(ptsInCluster,2,distMeas)
            sseSplit=sum(splitClustAss[:,1]) #被二分后的簇的平方误差和
            sseNotSplit=sum(clusterAssment[nonzero(clusterAssment[:,0]!=i)[0],1]) #整体数据中，未被二分的簇的平方误差和
            if (sseSplit+sseNotSplit)<lowestSSE:
                bestCentToSplit=i
                bestNewCents=centerMat
                bestClustAss=splitClustAss.copy() #.copy()防止splitClustAss被覆盖时影响到bestClustAss
                lowestSSE=sseSplit+sseNotSplit
        #确定好本次迭代被二分的簇后，将被二分的数据对应的簇类更新
        bestClustAss[nonzero(bestClustAss[:,0].A!=0)[0],0]=len(centerList) #更新先后顺序很重要！
        bestClustAss[nonzero(bestClustAss[:,0].A==0)[0],0]=bestCentToSplit
        print('The bestCentToSplit is:',bestCentToSplit)
        print('The number of samples to split is',len(bestClustAss))
        centerList[bestCentToSplit]=bestNewCents[0,:].tolist()[0] #将被二分的原簇质心替换为二分后的质心之一
        centerList.append(bestNewCents[1,:].tolist()[0]) #将二分后的另一质心添加在质心列表末尾
        #将二分后的簇的数据归属更新到总记录中
        clusterAssment[nonzero(clusterAssment[:,0].A==bestCentToSplit)[0],:]=bestClustAss
    return mat(centerList),clusterAssment

测试

center_testSet22,clusterAssment_testSet22=biKMeans(dataMat2,3)

plt.figure(figsize=(6,6))
plt.scatter(dataMat2[:,0].T.tolist()[0],dataMat2[:,1].T.tolist()[0],c='pink',s=30)
plt.scatter(center_testSet22.T[0].tolist()[0],center_testSet22.T[1].tolist()[0],c='blue',s=50)

---

04 实例：上车点规划

我们在前节实现了KMeans算法和BiKMeans算法，现在让我们来用用这两个轮子吧。

设想，你邀请了70个朋友参加你的party，他们住在城市的各个地方，你需要在party开始前2小时开车去接他们(假设你的车子可以容纳70人，哈哈哈)。

为了时间效率最大化，请问你该如何规划每个人的上车点呢？总不能你一个一个去家门口接吧。

我们用聚类的方法来规划！

现在我们获取了每个朋友住址的经纬度，那么请开始你的表演。

很简单，只需要调用我们早好的轮子BiKMeans

#球面距离计算函数：球面余弦距离(向量夹角*地球半径)
#求球面上两向量vecA,vecB的距离
def distSLC(vecA,vecB):
    a=sin(vecA[0,1]*pi/180)*sin(vecB[0,1]*pi/180) #由于抓取的经纬度为角度，需要通过 *pi/180来转换为弧度
    b=cos(vecA[0,1]*pi/180)*cos(vecB[0,1]*pi/180)*cos(pi*(vecA[0,0]-vecB[0,0])/180)
    return 6371.0*arccos(a+b) #6371为地球半径，单位为英尺

def clusterClubs(k=5):
    dataList=[]
    for line in open(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\places.txt').readlines():
        lineArr=line.strip().split('\t')
        dataList.append([float(lineArr[4]),float(lineArr[3])]) #读取地址的经纬度
    dataMat=mat(dataList)
    center,clusterAss=biKMeans(dataMat,k,distMeas=distSLC) #将地址按经纬度聚类
    #作图
    fig=plt.figure(figsize=(10,8))
    rect=[0.1,0.1,0.8,0.8] #用于设置坐标轴刻度，[]中前两个值表示左边起始位置，后两个值对应坐标长度
    scatterMarkers=['^','o','h','8','p','d','v','s','>','<'] #用于设置散点图点的形状
    axprops=dict(xticks=[],yticks=[])
    ax0=fig.add_axes(rect,label='ax0',**axprops)
    
    #读图，并将图片显示在设定好的坐标轴中
    imgP=plt.imread(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\Portland.png')
    ax0.imshow(imgP)
    #将地址按聚类结果作散点图
    ax1=fig.add_axes(rect,label='ax1',frameon=False)
    for i in range(k):
        ptsInCluster=dataMat[nonzero(clusterAss[:,0].A==i)[0],:]
        markerStyle=scatterMarkers[i]
        ax1.scatter(ptsInCluster[:,0].flatten().A[0],ptsInCluster[:,1].flatten().A[0],\
                   marker=markerStyle,s=90)
    #标出质心
    ax1.scatter(center[:,0].flatten().A[0],center[:,1].flatten().A[0],marker='+',s=300)
    plt.show()

好了，我们写的这个“上车点规划器”就可以使用了，只需要输入聚类簇数量即可 ，下面看看几个测试结果：

聚为5类（5个上车点） 聚为9类（9个上车点）

但是到底聚为几类比较合理呢？这个合理是指即不让你的朋友跑太远去坐车，也不用你跑太多地方去接人。

---

05 如何挑选最佳的聚类簇数

不论是KMeans算法还是BiKMeans算法，都仍有一个缺点没有解决：簇数量k由用户指定，用户指定的k不一定是最佳的簇数量。这也是上一节我们没有解决的问题。

怎么办呢？

使用轮廓系数，最佳k值由轮廓系数S(i)决定，取轮廓系数最大的分类数。

什么是轮廓系数，轮廓系数S(i)用于衡量聚类效果，取值在-1~1之间，越接近1表示聚类效果越好，公式如下，

其中a(i)=i点到同簇各点距离均值；b(i)=min(i点到某个非同簇各点均值)，即i点到其他簇质心距离的最小值，整体数据集的轮廓系数是各点轮廓系数的平均值

基于此，我们可以写一个计算聚类结果轮廓系数的函数，然后看看轮廓系数与聚类数量k的关系，从而找到最佳的聚类数量k。

轮廓系数计算函数

def outlineOfCluster(filename,maxClusterNum,distMeas=distEclud):
    dataList=[]
    for line in open(filename).readlines():
        lineArr=line.strip().split('\t')
        dataList.append([float(lineArr[4]),float(lineArr[3])]) #读取地址的经纬度
    dataMat=mat(dataList)
    #遍历2~nn个质心，求不同数量的簇的轮廓系数，设置轮廓系数最大的簇数量为k值
    sk={}
    for k in range(2,maxClusterNum+1):
        center,clusterAss=biKMeans(dataMat,k,distMeas)
        outline=[] #代表聚类为k个簇时各点的轮廓系数列表
        for i in range(k): #遍历各簇i
            ptsInCluster=clusterAss[nonzero(clusterAss[:,0].A==i)[0],:]
            for j in range(len(ptsInCluster)): #遍历i簇的各点j
                ptsInClusterNotJ=vstack([ptsInCluster[:j,:],ptsInCluster[j+1:,:]])
                ajn=[]
                for n in range(len(ptsInClusterNotJ)): #遍历i簇非j的点
                    ajn.append(distSLC(ptsInCluster[n],ptsInCluster[j]))
                aj=nanmean(ajn)
                bjm=[]
                centerWithoutI=vstack([center[:i,:],center[i+1:,:]])
                for m in range(len(centerWithoutI)): #遍历非i簇质心
                    bjm.append(distSLC(centerWithoutI[m],ptsInCluster[j]))
                bj=min(bjm)
                sj=(bj-aj)/max(bj,aj) #i簇中j点的轮廓系数
                outline.append(sj) #将i簇中各点的轮廓系数保存在outline中，outline用于存储聚类为k个簇时各点的轮廓系数，在遍历k时需要重置
        sk[k]=nanmean(outline) #聚类为k个簇时的轮廓系数是各点轮廓系数的均值
    return sk

现在我们可以接近上车点规划问题的遗留问题了，规划几个上车点最合理？

sk=outlineOfCluster(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\places.txt',35,distMeas=distSLC)

plt.figure(figsize=(10,5))
plt.plot(list(sk.keys()),list(sk.values()),linewidth=3)
plt.title('聚类簇数-效果',fontsize=18,fontweight='bold')
plt.xlabel('最终聚类簇数量',fontsize=14)
plt.ylabel('轮廓系数',fontsize=14)

可以看到，聚类簇数在达到10个簇之后，轮廓系数的增量就变得非常小了，因此从性价比方面考虑，选择聚类为10个簇的性价比是最高的，即规划10个上车点就好。

---

06 总结

本文实践了KMeans聚类算法，分成以下几个部分，

1. KMeans算法理论和代码实现
2. 改进，BiKMeans算法理论和代码实现
3. 实例，上车点规划
4. 抉择，如何挑选最佳的聚类簇数？

下期我们将实践Apriori算法，这是一种挖掘关联规则得到频繁项集的算法，比如可以预测你买了A商品后会买B商品，敬请期待~

---

07 参考

《机器学习实战》 Peter Harrington Chapter10