《机器学习实战》 Ch10 K-均值聚类

1、K均值聚类算法

　　聚类是一种无监督的学习，它将相似的对象归到同一簇中。聚类的方法几乎可以应用所有对象，簇内的对象越相似，聚类的效果就越好。K-均值是发现给定数据集的k个簇的算法。簇个数k是用户给定的，每一个簇通过其质心(centroid) , 即簇中所有点的中心来描述，即聚类中数据值的均值。
　　K-均值算法的工作流程是这样的。首先，随机确定k个初始点作为质心。然后将数据集中的每个点分配到一个簇中，具体来讲，为每个点找距其最近的质心， 并将其分配给该质心所对应的簇。这一步完成之后，每个簇的质心更新为该簇所有点的平均值。
上述过程的伪代码表示如下：

创建k个点作为起始质心（经常是随机选择）
当任意一个点的簇分配结果发生改变时
  对数据集中的每个数据点
    对每个质心
       计算质心与数据点之间的距离
    将数据点分配到距其最近的簇
  对每一个簇，计算簇中所有点的均值并将均值作为质心

实际代码：

from numpy import *
#K-均值聚类辅助函数
#将文本文件导人到一个列表
def loadDataSet(fileName):      #general function to parse tab -delimited floats
    dataMat = []                #assume last column is target value
    fr = open(fileName)
    for line in fr.readlines():
        curLine = line.strip().split('\t')
        #将每一行的数据映射成float型
        fltLine = map(float,curLine) #map all elements to float()
        dataMat.append(fltLine)
    return dataMat

def distEclud(vecA, vecB):#计算两个向量的欧式距离
    return sqrt(sum(power(vecA - vecB, 2))) #la.norm(vecA-vecB)

def randCent(dataSet, k):
#为给定数据集构建一个包含k个随机质心的集合
    n = shape(dataSet)[1] #得到数据样本的维度
    #初始化为一个(k,n)的矩阵
    centroids = mat(zeros((k,n)))#create centroid mat
    #遍历数据集的每一维度
    for j in range(n):
    #create random cluster centers, within bounds of each dimension
    #构造簇质心
        #得到该列数据的最小值
        minJ = min(dataSet[:,j]) 
        #得到该列数据的范围(最大值-最小值)
        rangeJ = float(max(dataSet[:,j]) - minJ)
        #k个质心向量的第j维数据值随机为位于(最小值，最大值)内的某一值
        centroids[:,j] = mat(minJ + rangeJ * random.rand(k,1))
    #返回初始化得到的k个质心向量
    return centroids

　　K均值算法会创建k个质心，然后将每个点分配到最近的质心，再重新计算质心。这个过程重复数次，直到数据点的簇分配结果不再改变为止。

#k-均值聚类算法
#@dataSet:聚类数据集
#@k:用户指定的k个类
#@distMeas:距离计算方法，默认欧氏距离distEclud()
#@createCent:获得k个质心的方法，默认随机获取randCent()
def kMeans(dataSet, k, distMeas=distEclud, createCent=randCent):
    #获取数据集样本数
    m = shape(dataSet)[0]
    #初始化一个(m,2)的矩阵
    clusterAssment = mat(zeros((m,2)))
    #create mat to assign data points 
    #to a centroid, also holds SE of each point
    #创建初始的k个质心向量               
    centroids = createCent(dataSet, k)
    #聚类结果是否发生变化的布尔类型
    clusterChanged = True
    #只要聚类结果一直发生变化，就一直执行聚类算法，直至所有数据点聚类结果不变化
    while clusterChanged:
        #聚类结果变化布尔类型置为false
        clusterChanged = False
        #遍历数据集每一个样本向量
        for i in range(m):
        #for each data point assign it to the closest centroid
            #初始化最小距离最正无穷；最小距离对应索引为-1
            minDist = inf; minIndex = -1
            #循环k个类的质心
            for j in range(k):  #寻找最近的质心
                #计算数据点到质心的欧氏距离
                distJI = distMeas(centroids[j,:],dataSet[i,:])
                #如果距离小于当前最小距离
                if distJI < minDist:
                    #当前距离定为当前最小距离；最小距离对应索引对应为j(第j个类)
                    minDist = distJI; minIndex = j
            #当前聚类结果中第i个样本的聚类结果发生变化：布尔类型置为true，继续聚类算法
            if clusterAssment[i,0] != minIndex: clusterChanged = True
            #更新当前变化样本的聚类结果和平方误差
            clusterAssment[i,:] = minIndex,minDist**2
        #打印k-均值聚类的质心
        print centroids
        #遍历每一个质心
        for cent in range(k):#recalculate centroids
        #更新质心的位置
            #将数据集中所有属于当前质心类的样本通过条件过滤筛选出来
            ptsInClust = dataSet[nonzero(clusterAssment[:,0].A==cent)[0]]#get all the point in this cluster
            #计算这些数据的均值（axis=0：求列的均值），作为该类质心向量
            centroids[cent,:] = mean(ptsInClust, axis=0) #assign centroid to mean 
    #返回k个聚类，聚类结果及误差
    return centroids, clusterAssment

　　上述清单给出了K-均值算法。kMeans()函数接受4个输人参数。只有数据集及簇的数目是必选参数，而用来计算距离和创建初始质心的函数都是可选的。kMeans()函数一开始确定数据集中数据点的总数，然后创建一个矩阵来存储每个点的簇分配结果。簇分配结果矩阵clusterAssment包含两列：一列记录簇索引值，第二列存储误差。这里的误差是指当前点到簇质心的距离，后边会使用该误差来评价聚类的效果。
　　按照上述方式(即计算质心-分配-重新计算)反复迭代，直到所有数据点的簇分配结果不再改变为止。程序中可以创建一个标志变量clusterChanged，如果该值为True，则继续迭代。上述迭代使用while循环来实现。接下来遍历所有数据找到距离每个点最近的质心，这可以通过对每个点遍历所有质心并计算点到每个质心的距离来完成。计算距离是使用distMeas参数给出的距离函数，默认距离函数是disEclud()。如果任一点的簇分配结果发生改变，则更新clusterChanged标志。
　　最后，遍历所有质心并更新它们的取值。具体实现步骤如下：首先通过数组过滤来获得给定簇的所有点；然后计算所有点的均值，选项axis=0表示沿矩阵的列方向进行均值计算；最后，程序返回所有的类质心与点分配结果。

>>> datMat=mat(kMeans.loadDataSet('testSet.txt'))
>>> myCentroids,clustAssing=kMeans.kMeans(datMat,4)
[[ 2.63373767  1.08576595]
 [-1.95057116  2.19323976]
 [-4.64423117  2.97559065]
 [-0.43325377  3.85644893]]
[[ 2.77226214 -0.46074594]
 [-2.91251273 -0.5788938 ]
 [-4.16112333  1.16232167]
 [ 0.5314535   3.91547925]]
[[ 2.89827912 -1.81544472]
 [-3.34887714 -2.76960033]
 [-3.0002845   2.75265   ]
 [ 1.16289995  3.40272823]]
[[ 2.8692781  -2.54779119]
 [-3.38237045 -2.9473363 ]
 [-2.54951105  2.75812458]
 [ 2.3772111   3.2195035 ]]
[[ 2.80293085 -2.7315146 ]
 [-3.38237045 -2.9473363 ]
 [-2.46154315  2.78737555]
 [ 2.6265299   3.10868015]]
>>> kMeans.showCluster(datMat,4,myCentroids,clustAssing)

　　上面的结果给出了4个质心，可以看出，经过5次迭代之后K-均值算法收敛。其中聚类结果示意代码：

 # 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() 

　　这4个质心以及原始数据的散点图如图：

2、二分K均值算法

　　为克服K-均值算法收敛于局部最小值的问题，有人提出了另一个称为二分K-均值（bisecting K-means)的算法。该算法首先将所有点作为一个簇，然后将该簇一分为二。之后选择其中一个簇继续进行划分，选择哪一个簇进行划分取决于对其划分是否可以最大程度降低SSE的值。上述基于SSE的划分过程不断重复，直到得到用户指定的簇数目为止。
　　一种用于度量聚类效果的指标是SSE(Sum of Squared Error，误差平方和），对应clusterAssment矩阵的第二列之和。SSE值越小表示数据点越接近于它们的质心，聚类效果也越好。因为对误差取了平方，因此更加重视那些远离中心的点。一种肯定可以降低SSE值的方法是增加簇的个数，但这违背了聚类的目标。聚类的目标是在保持族数目不变的情况下提高簇的质量。
　　二分&-均值算法的伪代码形式为：

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

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

#二分K-均值聚类算法
#@dataSet:待聚类数据集
#@k：用户指定的聚类个数
#@distMeas:用户指定的距离计算方法，默认为欧式距离计算
def biKmeans(dataSet, k, distMeas=distEclud):
    #获得数据集的样本数
    m = shape(dataSet)[0]
    #初始化一个元素均值0的(m,2)矩阵
    clusterAssment = mat(zeros((m,2)))#创建一个初始簇
    #获取数据集每一列数据的均值，组成一个长为列数的列表
    centroid0 = mean(dataSet, axis=0).tolist()[0]
    #当前聚类列表为将数据集聚为一类
    centList =[centroid0] #create a list with one centroid
    #遍历每个数据集样本
    for j in range(m):#calc initial Error
        #计算当前聚为一类时各个数据点距离质心的平方距离
        clusterAssment[j,1] = distMeas(mat(centroid0), dataSet[j,:])**2
    #循环，直至二分k-均值达到k类为止
    while (len(centList) < k):
        #将当前最小平方误差置为正无穷
        lowestSSE = inf
        #遍历当前每个聚类
        for i in range(len(centList)):#尝试划分每一簇
            #通过数组过滤筛选出属于第i类的数据集合
            ptsInCurrCluster = dataSet[nonzero(clusterAssment[:,0].A==i)[0],:]
            #get the data points currently in cluster i
            #对该类利用二分k-均值算法进行划分，返回划分后结果，及误差
            centroidMat, splitClustAss = kMeans(ptsInCurrCluster, 2, distMeas)
            #计算该类划分后两个类的误差平方和
            sseSplit = sum(splitClustAss[:,1])#compare the SSE to the currrent minimum
            #计算数据集中不属于该类的数据的误差平方和
            sseNotSplit = sum(clusterAssment[nonzero(clusterAssment[:,0].A!=i)[0],1])
            #打印这两项误差值
            print "sseSplit, and notSplit: ",sseSplit,sseNotSplit
            #划分第i类后总误差小于当前最小总误差
            if (sseSplit + sseNotSplit) < lowestSSE:
                #第i类作为本次划分类
                bestCentToSplit = i
                #第i类划分后得到的两个质心向量
                bestNewCents = centroidMat
                #复制第i类中数据点的聚类结果即误差值
                bestClustAss = splitClustAss.copy()
                #将划分第i类后的总误差作为当前最小误差
                lowestSSE = sseSplit + sseNotSplit
                #更新簇的分配结果
        #数组过滤筛选出本次2-均值聚类划分后类编号为1数据点，将这些数据点类编号变为
        #当前类个数+1，作为新的一个聚类
        bestClustAss[nonzero(bestClustAss[:,0].A == 1)[0],0] = len(centList) #change 1 to 3,4, or whatever
        #同理，将划分数据集中类编号为0的数据点的类编号仍置为被划分的类编号，使类编号
        #连续不出现空缺
        bestClustAss[nonzero(bestClustAss[:,0].A == 0)[0],0] = bestCentToSplit
        #打印本次执行2-均值聚类算法的类
        print 'the bestCentToSplit is: ',bestCentToSplit
        #打印被划分的类的数据个数
        print 'the len of bestClustAss is: ', len(bestClustAss)
        #更新质心列表中的变化后的质心向量
        centList[bestCentToSplit] = bestNewCents[0,:].tolist()[0]
        #replace a centroid with two best centroids 
        #添加新的类的质心向量
        centList.append(bestNewCents[1,:].tolist()[0])
        #更新clusterAssment列表中参与2-均值聚类数据点变化后的分类编号，及数据该类的误差平方
        clusterAssment[nonzero(clusterAssment[:,0].A == bestCentToSplit)[0],:]= bestClustAss
        #reassign new clusters, and SSE
    #返回聚类结果
    return mat(centList), clusterAssment

　　该函数首先创建一个矩阵来存储数据集中每个点的簇分配结果及平方误差，然后计算整个数据集的质心，并使用一个列表来保留所有的质心。得到上述质心之后，可以遍历数据集中所有点来计算每个点到质心的误差值。这些误差值将会在后面用到。
　　一开始将最小SSE置设为无穷大，然后遍历簇列表centList中的每一个簇 。对每个簇，将该簇中的所有点看成一个小 的数据集ptsInCurrCluster。将ptsInCurrCluster输人到函数kMeans()中进行处理 (k=2)。K-均值算法会生成两个质心(簇)，同时给出每个簇的误差值。这些误差与剩余数据集的误差之和作为本次划分的误差。如果该划分的SSE值最小，则本次划分被保存。一旦决定了要划分的簇，接下来就要实际执行划分操作。划分操作很容易，只需要将要划分的簇中所有点的簇分配结果进行修改即可。
　　当使用kMeans()函数并且指定簇数为2时，会得到两个编号分别为0和1的结果簇。需要将这些簇编号修改为划分簇及新加簇的编号，该过程可以通过两个数组过滤器来完成。最后，新的簇分配结果被更新，新的质心会被添加到centList中。

>>> datMat3=mat(kMeans.loadDataSet('testSet2.txt'))
>>> centlist,myNewAssments=kMeans.biKmeans(datMat3,3)
[[-0.15316589  1.20439814]
 [-4.73761553  0.52399209]]
[[ 0.78986789  0.65484736]
 [-3.58362738  3.28784469]]
[[ 1.19084544  0.24642754]
 [-3.06779095  3.33769884]]
[[ 1.23710375  0.17480612]
 [-2.94737575  3.3263781 ]]
sseSplit, and notSplit:  570.722757425 0.0
the bestCentToSplit is:  0
the len of bestClustAss is:  60
[[ 2.28522325  0.1678996 ]
 [ 2.26581725 -3.85684349]]
[[ 2.41135835  2.56370726]
 [-0.35159365 -3.05723659]]
[[ 2.93386365  3.12782785]
 [-0.45965615 -2.7782156 ]]
sseSplit, and notSplit:  68.6865481262 38.0629506357
[[-1.11915506  3.9373216 ]
 [-1.51175944  2.23825046]]
[[-2.90993709  3.99469745]
 [-2.99313411  2.50954333]]
[[-3.1349821  4.0830388]
 [-2.7597694  2.5697174]]
[[-3.25187989  4.160664  ]
 [-2.698236    2.64378055]]
[[-3.39413187  4.2372255 ]
 [-2.64953833  2.7191465 ]]
[[-3.61124957  4.26650957]
 [-2.58990523  2.82015346]]
[[-4.095738    4.4204886 ]
 [-2.56458833  2.9616746 ]]
sseSplit, and notSplit:  21.2908596794 532.659806789
the bestCentToSplit is:  0
the len of bestClustAss is:  40

质心结果：

>>> centlist
matrix([[ 2.93386365,  3.12782785],
        [-2.94737575,  3.3263781 ],
        [-0.45965615, -2.7782156 ]])

　　上述函数可以运行多次，聚类会收敛到全局最小值，而原始的别的kMeans()函数偶尔会陷人局部最小值。