《机器学习实战》之二分K-均值聚类算法的python实现

上面博文介绍了K-均值聚类算法及其用python实现，上篇博文中的两张截图，我们可以看到，由于K-均值聚类算法中由于初始质心的选取，会造成聚类的局部最优，并不是全局最优，因此，会造成聚类的效果并不理想，为克服K-均值算法收敛于局部最小值的问题，就有了二分K-均值算法。

二分K-均值聚类算法

二分K均值算法是基本K均值算法的直接扩充，其基本思想是：为了得到K个簇，首先将所有点的集合分裂成两个簇，然后从这些簇中选取一个继续分裂，迭代直到产生K个簇；二分K均值的关键是如何选择待分裂簇，选择策略主要有两种，一是选择包含矢量个数最多的簇进行分裂，二是选择具有最大SSE（误差的平方和）的簇分裂。

二分k均值算法的伪代码如下：

将所有数据点看成一个簇
当簇数目小于k时
            对每一个簇
                    计算总误差
                    在给定的簇上面进行k-均值聚类（k=2）
                    计算将该簇一分为二后的总误差
            选择使得误差最小的那个簇进行划分操作

另一种做法是选择SSE最大的簇进行划分，直到簇数目达到用户指定的数目为止。

python实现如下

环境：python 3.4 使用的库：numpy、matplotlib

biKmeans.py:此文件中的kmeans算法与上篇博文中的代码一致，在上一篇博文中注释写的比较详细，这里就没有在详细注释，若需要看注释，可以看上一篇博文

# 二分k均值算法的伪代码如下：
#***************************************************************
#将所有数据点看成一个簇
#
#当簇数目小于k时
#
# 对每一个簇
#
# 计算总误差
#
# 在给定的簇上面进行k-均值聚类（k=2）
#
# 计算将该簇一分为二后的总误差
#
# 选择使得误差最小的那个簇进行划分操作
#
#*************************************************************** 

from numpy import *  
import time  
import matplotlib.pyplot as plt  


# calculate Euclidean distance 
def euclDistance(vector1, vector2):  
    return sqrt(sum(power(vector2 - vector1, 2)))  

# init centroids with random samples 
def initCentroids(dataSet, k):  
    numSamples, dim = dataSet.shape  
    centroids = zeros((k, dim))  
    for i in range(k):  
        index = int(random.uniform(0, numSamples)) 
        print(index)
        centroids[i, :] = dataSet[index, :]  
    return centroids  

# k-means cluster 
def kmeans(dataSet, k):  
    numSamples = dataSet.shape[0]  
    # first column stores which cluster this sample belongs to, 
    # second column stores the error between this sample and its centroid 
    clusterAssment = mat(zeros((numSamples, 2)))  
    clusterChanged = True  

    ## step 1: init centroids 
    centroids = initCentroids(dataSet, k)  

    while clusterChanged:  
        clusterChanged = False  
        ## for each sample 
        for i in range(numSamples):  
            minDist  = 100000.0  
            minIndex = 0  
            ## for each centroid 
            ## step 2: find the centroid who is closest 
            for j in range(k):  
                distance = euclDistance(centroids[j, :], dataSet[i, :])  
                if distance < minDist:  
                    minDist  = distance  
                    minIndex = j  

            ## step 3: update its cluster 
            if clusterAssment[i, 0] != minIndex:  
                clusterChanged = True  
                clusterAssment[i, :] = minIndex, minDist**2  

        ## step 4: update centroids 
        for j in range(k):  
            pointsInCluster = dataSet[nonzero(clusterAssment[:, 0].A == j)[0]]  
            centroids[j, :] = mean(pointsInCluster, axis = 0)  

    print ('Congratulations, cluster using k-means complete!'  )
    return centroids, clusterAssment  

# bisecting k-means cluster 
def biKmeans(dataSet, k):  
    numSamples = dataSet.shape[0]  
    # first column stores which cluster this sample belongs to, 
    # second column stores the error between this sample and its centroid 

    clusterAssment = mat(zeros((numSamples, 2)))  

    # step 1: the init cluster is the whole data set 
    #mean(dataSet,axis=0)返回的是：matrix([[XXX,XXX,XXX]])
    #mean(dataSet, axis = 0).tolist()返回的是：[[XXX,XXX,XXX]],
    #mean(dataSet, axis = 0).tolist()[0],返回的是：[XXX,XXX,XXX]
    centroid0 = mean(dataSet, axis = 0).tolist()[0]  #列表中的第一个列表元素：即全部数据每个属性对应的均值
    centList = [centroid0]  #centList是[[xxx,xxx,xxx]]
    for j in range(numSamples):  
        clusterAssment[j, 1] = (euclDistance(mat(centroid0), dataSet[j, :]))**2 #计算每个样本点与质心之间的距离 

    while (len(centList) < k):  
        # min sum of square error 
        minSSE = inf 
        #numCurrCluster = len(centList) #当前的簇的个数
        # for each cluster 
        #找出numCurrCluster个簇中哪个簇分解得到的误差平方和最小的两个簇
        for i in range(len(centList)):  
            # step 2: get samples in cluster i
            #选取第i个簇的所有数据，然后将其分成两个两个簇
            pointsInCurrCluster = dataSet[nonzero(clusterAssment[:, 0].A == i)[0], :]  

            # step 3: cluster it to 2 sub-clusters using k-means 
            #centroids的元素为每个簇的质心
            #splitClusterAssment第一列为样本所属的类别号，第二列为样本到其所属簇的质心的距离的平方

            #每个2均值聚类时循环N次，去除kmeans的初始值随机选取不当
            N=20
            minSSE_kmeans=100000.0;
            for j in range(N):
                centroids_kmeans, splitClusterAssment_kmeans = kmeans(pointsInCurrCluster, 2)
                if  (sum(splitClusterAssment_kmeans[:, 1]))<minSSE_kmeans:
                    minSSE_kmeans=sum(splitClusterAssment_kmeans[:, 1])
                    splitClusterAssment=splitClusterAssment_kmeans
                    centroids=centroids_kmeans


            # step 4: calculate the sum of square error after split this cluster
            #下面的代码是求误差平方和
            #splitSSE=sum(power(splitClusterAssment[:,1],2))
            splitSSE = sum(splitClusterAssment[:, 1]) 
            #不是标号为第i个簇的误差平方和
            notSplitSSE = sum(clusterAssment[nonzero(clusterAssment[:, 0].A!=i)[0], 1])  
            currSplitSSE = splitSSE + notSplitSSE  #当前所有簇的平方和 
            print('num=%d,%s,%s,%s' %(i,splitSSE,notSplitSSE,currSplitSSE))
            # step 5: find the best split cluster which has the min sum of square error 
            if currSplitSSE < minSSE:  #
                bestCentroidToSplit = i
                bestNewCentroids = centroids
                bestClusterAssment = splitClusterAssment.copy()
                minSSE = currSplitSSE  



        print('minSSE=%s' %minSSE)
        # step 6: modify the cluster index for adding new cluster
        #将新分出来的两个簇的标号一个沿用它父亲的标号，一个用簇的总数来标号。
        bestClusterAssment[nonzero(bestClusterAssment[:, 0].A == 1)[0], 0] = len(centList)  
        bestClusterAssment[nonzero(bestClusterAssment[:, 0].A == 0)[0], 0] = bestCentroidToSplit  

        # step 7: update and append the centroids of the new 2 sub-cluster 

        centList[bestCentroidToSplit] = bestNewCentroids[0, :] #将第一个子簇的质心放在父亲质心的原位置 
        centList.append(bestNewCentroids[1, :]) #将第二个子簇的质心添加在末尾 

        # step 8: update the index and error of the samples whose cluster have been changed
        #由第i个簇分解为j、k两个簇所得到的数据将分解之前的数据替换掉
        clusterAssment[nonzero(clusterAssment[:, 0].A == bestCentroidToSplit)[0], :] = bestClusterAssment  

    print ('Congratulations, cluster using bi-kmeans complete!' ) 
    return mat(centList), clusterAssment  

# show your cluster only available with 2-D data 
def showCluster(dataSet, k, centroids, clusterAssment):  
    numSamples, dim = dataSet.shape  
    if dim != 2:  
        print ("Sorry! I can not draw because the dimension of your data is not 2!"  )
        return 1  

    mark = ['or', 'ob', 'og', 'ok', '^r', '+r', 'sr', 'dr', '<r', 'pr']  
    if k > len(mark):  
        print ("Sorry! Your k is too large!")  
        return 1  

    # draw all samples 
    for i in range(numSamples):  
        markIndex = int(clusterAssment[i, 0])  
        plt.plot(dataSet[i, 0], dataSet[i, 1], mark[markIndex])  

    mark = ['Dr', 'Db', 'Dg', 'Dk', '^b', '+b', 'sb', 'db', '<b', 'pb']  
    # draw the centroids 
    for i in range(k):  
        plt.plot(centroids[i, 0], centroids[i, 1], mark[i], markersize = 12)  

    plt.show() 

测试文件test.py:与K均值聚类这篇博文中的代码相同，这里也就没有再注释

################################################# 
# kmeans: k-means cluster 
# Author :wojiushimogui 
# Date :2015年7月28日11:12:44
# HomePage :http://write.blog.csdn.net/postlist
# Email : 
################################################# 
from numpy import *  
import time  
import matplotlib.pyplot as plt  
import biKMeans
from biKMeans import *
## step 1: load data 
print ("step 1: load data..." ) 
dataSet = []  
fileIn = open('D:/xuepython/BiKmeans/testSet.txt')  
for line in fileIn.readlines():  
    lineArr = line.strip().split('\t')  
    dataSet.append([float(lineArr[0]), float(lineArr[1])])  

## step 2: clustering... 
print ("step 2: clustering...")  
dataSet = mat(dataSet)  
k = 4  
centroids, clusterAssment = biKmeans(dataSet, k)  

## step 3: show the result 
print ("step 3: show the result..."  )
showCluster(dataSet, k, centroids, clusterAssment) 

运行结果截图如下：

是乃遗憾，运行的结果也不是我们想象的那样，聚类的结果也是随机的，不是那么的理想，但是，按照《机器学习实战》这本书上面提供的二分K-均值聚类算法的思想：找到分解使得总的SSE最小（即找到那个簇使得分解后能够最大的降低SSE））来看，应该是不会出现上面的问题的，但是今天无论是我参考《机器学习实战》这本书上面的源代码，还是这篇博客中的源码，都不能解决这个问题。

由于在参考《数据挖掘导论》这本书上面的二分K均值聚类时看到如下的伪代码。

初始化簇表，使之包含由所有点组成的簇
repeat
    从簇表中取出一个簇
    {对选定的簇进行多次二分“实验”}
    for i= 1 to 实验次数 do
         使用基本K均值，二分选定的簇。
    end for
    从二分试验中选择具有最小总SSE的两个簇
    将这两个簇添加到簇表中。
until 簇表中包含K个簇

因此，按照上面的思路，我在《机器学习实战》上面的源码中添加了如下代码。

#每个2均值聚类时循环N次，去除kmeans的初始值随机选取不当 N=20 minSSE_kmeans=100000.0; for j in range(N): centroids_kmeans, splitClusterAssment_kmeans = kmeans(pointsInCurrCluster, 2) if (sum(splitClusterAssment_kmeans[:, 1]))<minSSE_kmeans: minSSE_kmeans=sum(splitClusterAssment_kmeans[:, 1]) splitClusterAssment=splitClusterAssment_kmeans centroids=centroids_kmeans

发现，依然还是没有办法解决这个问题。

马上就要去跑步了，今天就思考到这里，明天我将会换一种方法来解决这个问题：通过找到最大的SSE的那个簇，然后将其分解；按照这样一个思路我相信应该可以解决这个问题。

源码和测试数据在这里可以获取

身体是自己的，科研是老板的，哈哈
跑步去了