k-means算法

k-means算法是一种非常简单的聚类算法，也就是把信息分类的一种算法。
具体算法思路为：
（这里以二维平面坐标系为例）

1. 在平面上随机投入k个点种子点，表示要把二维平面上的点集分为k个部分。
2. 对二维平面上面的所有点求到这k个点的距离，对于二维平面上的所有信息点，找出它距离这k个种子点中最近的那个，把此信息点归类到该种子点的分类当中。
3. 把这k个种子点移动到它所属于的“点群”的中心。
4. 重复上面的第二、第三步。直到没有点被移动位置，说明最后分好类了。

如图，在点集当中放入5个种子点后，点的移动情况。

图片来自http://www.csdn.net/article/2012-07-03/2807073-k-means

简单粗暴的解决方法虽然很好理解，但是缺点通常也很明显。
k-means算法的缺点如下

• 对异常或者错误的信息值非常敏感
• 需要直到k的值，也就是要直到将信息分成多少部分。
• 经常陷入到局部最优解，要做多次实现

优点如下：

• 速度快
• 方法简单

描述k-means算法划分出来的结果可以用SSE （Sum of squared errors 误差平方和）作为评定标准。

对于k-means算法经常陷入到局部最优解的问题，可以用如下的解决方案。
1.多次使用k-means算法，找到SSE最小的那个值。
2.使用二分k均值、k-means++算法等等

SSE是误差平方和，可以从图像当中得知，如果种子点找的很好，也就是分出的k类聚簇到这k个聚簇的种子点距离的平方和最小，那么就很可能是最优的（可能局部最优）

使用多次计算寻找SSE最小的方案的代码

#!/usr/bin/python
# -*- coding: utf-8 -*-

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

dataSet= [] 
fileIn = open("in.txt")
for line in fileIn.readlines():
    lineArr = line.strip().split(' ')
    dataSet.append([int(lineArr[0]), int(lineArr[1])]) 

dataSet=mat(dataSet)

#计算两点之间距离
def euclDistance(vector1, vector2):
    return sqrt(sum(power(vector2 - vector1, 2)))

#初始化种子点
def initCentroids(dataSet, k):
    numSamples, dim = dataSet.shape #行数 列数 
    centroids = zeros((k, dim))  

    for i in range(k):  
        index = int(random.uniform(0, numSamples))  
        centroids[i, :] = dataSet[index, :]  
    return centroids



def kmeans(dataSet, k):  
    numSamples = dataSet.shape[0]  
    clusterAssment = mat(zeros((numSamples, 2)))  
    clusterChanged = True
    SSE =  zeros(numSamples)#计算 SSE

    centroids = initCentroids(dataSet, k)  

    while clusterChanged:  
        clusterChanged = False  

        for i in xrange(numSamples):  
            minDist  = inf  
            minIndex = 0  
            #对于每个种子点寻找距离它最近的点集
            for j in range(k):  
                distance = euclDistance(centroids[j, :], dataSet[i, :])  
                if distance < minDist:  
                    minDist  = distance  
                    minIndex = j  

            if clusterAssment[i, 0] != minIndex:  
                clusterChanged = True  
                clusterAssment[i, :] = minIndex, minDist**2  
            SSE [i] = euclDistance(dataSet[i], centroids[minIndex])**2

        #更新种子点位置
        for j in range(k):  
            pointsInCluster = dataSet[nonzero(clusterAssment[:, 0].A == j)[0]]  
            centroids[j, :] = mean(pointsInCluster, axis = 0)  
    return centroids, clusterAssment, sum(SSE)


#输出结果图
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! please contact Zouxy"  
        return 1  

    # draw all samples  
    for i in xrange(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() 

k=3


best_clusterAssement = None
best_centriods = None
best_assement = inf

#计算10次，找出SSE最小的那个解
for i in range(10):
    centroids, clusterAssment, SSE=kmeans(dataSet, k)
    if SSE < best_assement :
        best_assement = SSE
        best_centriods = centroids
        best_clusterAssement = clusterAssment

showCluster(dataSet, k, best_centriods, best_clusterAssement)

fileIn.close()

运行效果图

使用二分k-means算法

二分k-means算法的大致思路就是，把上一次分好的所有聚类当中，挑选出SSE最大的。然后然后把挑选出来的聚类再次划分成两个聚类。计算SSE值，如果SSE值变小，说明可以继续”分堆”，如果变大，说明这个聚类已经分好了。

使用k-means++算法

k-means++算法解决初始时种子点位置选取的问题，不同的位置会造成解决的不同。
基本思路如下。

1.在所有点集当中随机找一个当种子点。
2.对于数据点中的所有点计算到距离它最近的种子点的距离，记为D(x)
3.把所有算出的距离长度求和，记为S.
4.按照类似于轮盘选赌的方法，D(x)值越大的点，被选中当种子点的概率越高。
5.重复上面的步骤，直到找到k个种子点，划分成k个聚类为止。

具体代码，稍后补上。

in.txt的数据

769086 392127
769203 362026
767856 387596
772142 360027
777968 404006
752502 343589
773724 342866
703608 385847
769086 392127
731388 381936
753989 341346
780265 39745
733019 394619
782848 391645
781523 396512
780161 398074
795446 427024
805407 397092
779634 405445
781514 391607
769147 387951
794752 403675
769086 392127
747727 400584
796195 385612
780883 401244
805407 397092
779658 391057
753228 339148
786474 414332
786474 414332
773582 39683
78673 400728
750978 339193
754848 338616
684158 396323
754848 338616
733019 394619
767856 387596
772142 360027
761365 332257
779131 404529
808645 432791
783634 399152
781763 39759
765153 358793
765148 358783
786474 414332
761365 332257
772746 366318
754848 338616
791093 396366
786474 414332
760878 37672
781266 390082
780883 401244
756873 40768
753267 339429
73368 396794
780921 394587
769203 362026
773659 394787
77832 406484
794752 403675
747727 400584
731071 399401
737278 410394
775402 333306
786474 414332
684158 396323
777968 404006
78073 401444
786474 414332
786474 414332
786474 414332
815248 426843
745256 407037
795446 427024
754848 338616
743049 413688
744344 404774
782706 401768
684158 396323
73357 384082
787616 387903
73357 384082
768709 391412
785568 419207
782609 398889

代码和内容引用来自
http://blog.csdn.net/zouxy09/article/details/17589329
http://www.jianshu.com/p/5314834f9f8e
http://www.jb51.net/article/49395.htm