二分-k均值算法

首先我们都知道k均值算法有一个炒鸡大的bug，就是在很多情况下他只会收敛到局部最小值而不是全局最小值，为了解决这个问题，很多学者提出了很多的方法，我们在这里介绍一种叫做2分k均值的方法。

该算法首先将所有点作为一个簇，然后将该簇一分为二。之后选择其中一个簇继续进行划分，选择哪一个簇进行划分取决于哪个簇的sse是最大值。上述基于sse的划分过程不断重复，直到得到用户指定的簇数目为止。

将所有的点看成一个簇，当粗的数目小于k时，对每一个簇计算总误差，在给定的粗上进行k均值聚类（k=2），计算将该粗一分为二之后的总误差。最后选择sse最大的簇进行划分。重复执行若干次直到簇的数目等于k。

 

1.首先贴出k均值函数

# coding=utf-8

from numpy import *
import matplotlib
import matplotlib.pyplot as plt
import operator
from os import listdir
import time


def distEclud(vecA, vecB):
    return sqrt(sum(power(vecA - vecB, 2)))  # la.norm(vecA-vecB)


def randCent(dataSet, k):
    n = shape(dataSet)[1]
    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)
        centroids[:, j] = mat(minJ + rangeJ * random.rand(k, 1))
    return centroids


def kMeans(dataSet, k, distMeas=distEclud, createCent=randCent):
    m = dataSet.shape[0]
    clusterAssment = zeros((m, 2))


    centroids = createCent(dataSet, k)
    # print centroids
    show(centroids)
    clusterChanged = True
    while clusterChanged:
        clusterChanged = False

        for i in range(m):
            point = dataSet[i, :]  # 遍历每个点
            mindist = inf
            minindex = -1

            for n in range(k):
                heart = centroids[n, :]  # 遍历每个质心
                distance = distMeas(point, heart)  # 求点与质心距离
                if distance < mindist:
                    mindist = distance  # 更新最小距离mindist
                    minindex = n  # 更新最小距离的质心序号

            if clusterAssment[i, 0] != minindex: clusterChanged = True
            clusterAssment[i, :] = minindex, mindist ** 2  # 方差
        # print clusterAssment
        for cent in range(k):
            ptsInClust = dataSet[(clusterAssment[:, 0] == cent)]  # get all the point in this cluster
            # print ptsInClust
            if len(ptsInClust):
                centroids[cent, :] = mean(ptsInClust, axis=0)  # assign centroid to mean
            else:
                centroids[cent, :] = array([[0, 0]])
        show(centroids,color='green')

    return centroids, clusterAssment



def show(data,color=None):
    if not color:
        color='green'
    group=createDataSet()
    fig = plt.figure(1)
    axes = fig.add_subplot(111)
    axes.scatter(group[:, 0], group[:, 1], s=40, c='red')
    axes.scatter(data[:, 0], data[:, 1], s=50, c=color)
    plt.show()



def createDataSet():
    group = array([[1.0, 1.1], [1.0, 1.0],
                   [0, 0], [0, 0.1],
                   [2, 1.0], [2.1, 0.9],
                   [0.3, 0.0], [1.1, 0.9],
                   [2.2, 1.0], [2.1, 0.8],
                   [3.3, 3.5], [2.1, 0.9],
                   [2, 1.0], [2.1, 0.9],
                   [3.5, 3.4], [3.6, 3.5]
                   ])


    return group
# centroids, clusterAssment=kMeans(createDataSet(),4)
# show(centroids,color='yellow')

　　

在此基础上我们加上二分的算法

def biKmeans(dataSet, k, distMeas=distEclud):
    m = shape(dataSet)[0]     #点数
    clusterAssment = mat(zeros((m,2)))   #空矩阵
    centroid0 = mean(dataSet, axis=0).tolist()  #数据集的平均值，tolist是转换为列表
    # print centroid0
    centList = [centroid0]  # create a list with one centroid  对每一个质心建立一个列表容器
    for j in range(m):     #遍历数据集
        clusterAssment[j,1]=distMeas(mat(centroid0),dataSet[j,:])**2  #求出数据集中每一个点到先前选定质心的距离平方
                                                      #并将其赋值给clusterAssment这个矩阵的对应列的第二个值。
                                                    #而第一个值全都赋值为0表示当前只有一个簇
        # print mat(centroid0),dataSet[j,:]
    # print clusterAssment

    while (len(centList) < k):  #簇的数量小于4时
        lowestSSE = inf            #初始化lowersse为正无穷
        for i in range(len(centList)):    # 遍历当前存在的每一个质心,整个遍历过程是一个找质心的过程，但即便是2分，得到的结果也是不确定的
                          #这个循环的目的只是得到划分哪个质心可以得到最大效益，而不考虑如何划分，主要的原因是k均值存在很
                         # 大偶然性，不能得到确切结果。
            print i

            ptsInCurrCluster = dataSet[nonzero(clusterAssment[:, 0].A == i)[0], :] # 每一个簇所拥有的所有数据集
            # print ptsInCurrCluster               #当前for 循环中只有一个簇
            centroidMat, splitClustAss = kMeans(ptsInCurrCluster, 2, distMeas)       #k-均值算法，k=2
            # print centroidMat
            # axes.scatter(centroidMat[:, 0], centroidMat[:, 1], s=40, c='blue')         #将得到的两个质心描绘出来
            # plt.show()
            # print splitClustAss                      #点分配到两个质心的分配方式矩阵
            sseSplit = sum(splitClustAss[:, 1])           #分配后的点方差之和   误差sse
            # print sseSplit


            sseNotSplit = sum(clusterAssment[nonzero(clusterAssment[:, 0].A != i)[0], 1])  #不在当前分配簇中的点方差之和

            print "sseSplit, and notSplit: ", sseSplit, sseNotSplit
            if (sseSplit + sseNotSplit) < lowestSSE:   # 两种方差和小于先前的lowersse，则说明这种分配方式减小了误差率
                bestCentToSplit = i                         #暂时将当前划分方式设为最佳
                bestNewCents = centroidMat            #暂时将当前的划分质心设为最好
                bestClustAss = splitClustAss.copy()     #复制当前划分的   点分配到两个质心的分配方式矩阵
                lowestSSE = sseSplit + sseNotSplit      #更新最小lowersse

          # 将得到的两个最佳质心描绘出来
        # show(bestNewCents)
        # print bestClustAss
        bestClustAss[nonzero(bestClustAss[:, 0]== 1)[0], 0] = len(centList)  # 2分k聚类返回系数0，或1  ，需要把1换成当前簇数目，以免造成重复
        bestClustAss[nonzero(bestClustAss[:, 0] == 0)[0], 0] = bestCentToSplit  #把0换成别切分的簇，或者与上面的交换赋值也可以
        print 'the bestCentToSplit is: ', bestCentToSplit
        print 'the len of bestClustAss is: ', len(bestClustAss)

        centList[bestCentToSplit] = bestNewCents[0, :].tolist()[0]  # 将centlist指定位置上的质心换为分割后的质心
        centList.append(bestNewCents[1, :].tolist()[0])                #将另一个质心添加上去

        clusterAssment[nonzero(clusterAssment[:, 0].A == bestCentToSplit)[0],:] = bestClustAss  # 将划分后的新质心及点分布赋值给结果矩阵
        # print clusterAssment
        show(mat(centList),color="blue")
    return mat(centList), clusterAssment

cent,clusterAssment= biKmeans(createDataSet(), 4)
show(cent,color='yellow')

这样我们就结束了2分k均值的编写，效果得到了很大的优化。