K-means原理及Python实现

K-means方法是一种非监督学习的算法，它解决的是聚类问题。

1、算法简介：K-means方法是聚类中的经典算法，数据挖掘十大经典算法之一；算法接受参数k,然后将事先输入的n个数据对象划分为k个聚类以便使得所获得的聚类满足聚类中的对象相似度较高，而不同聚类中的对象相似度较小。

2、算法思想：以空间中k个点为中心进行聚类，对最靠近他们的对象归类，通过迭代的方法，逐次更新各聚类中心的值，直到得到最好的聚类结果。

3、算法描述：

（1）适当选择c个类的初始中心； 
（2）在第k次迭代中，对任意一个样本，求其到c各中心的距离，将该样本归到距离最短的那个中心所在的类； 
（3）利用均值等方法更新该类的中心值； 
（4）对于所有的C个聚类中心，如果利用（2）（3）的迭代法更新后，值保持不变，则迭代结束；否则继续迭代。

4、算法举例：

我们假设药物A、B、C、D有两个特征值，分别是药物重量以及PH值。

药物名称	药物重量	药物PH值
A	1	1
B	2	1
C	4	3
D	5	4

现在我们要对这四个药物进行聚类，已知我们要分成两类，那么我们该怎么做呢？

首先我们把上面的数据画到二位坐标系当中A(1,1),B(2,1),C(4,3),D(5,4)：

初始时，我们先假设药物A为聚类1的中心点，B为聚类2的中心点，那么初始时的中心坐标分别为c1=(1,1),c2=(2,1)，矩阵D的第一行代表各个点到中心点c1c1的距离，第二行代表各个点到中心点c2c2的距离；那么初始矩阵D0D0表示成如下： 

矩阵GG代表样本应该归属于哪个聚类，第一行代表各个点是否属于中心c1c1所在的类（0代表不在，1代表在），第二行代表各个点是否属于中心c2c2所在的类（0代表不在，1代表在）；那么此时G0G0表示成如下： 

由矩阵G0G0可知A药物属于一个类，B、C、D属于一类；
然后，利用均值等方法更新该类的中心值。 
c1=（1,1）

c2=（(2+4+5)/3,(1+3+4)/3）=(13/3,8/3)

上图是更新后的坐标图，对应的中心点也发生了变化。

因为中心点跟上次不一样了，所以我们又可以对样本点进行重新划分。划分的方法还是跟以前一模一样，我们先计算出矩阵D1表示成如下： 

此时G1表示成如下： 

由矩阵G1可知A、B药物属于一个类，C、D属于一类；
然后，利用均值等方法再次更新该类的中心值。 
c1=（（1+2）/2,(1+1)/2）=(1.5,1)

c2=（(4+5)/2,(3+4)/2）=(4.5,3.5)

 

上图是更新后的坐标图，对应的中心点也发生了变化。

因为中心点跟上次不一样了，所以我们又可以对样本点进行重新划分。划分的方法还是跟以前一模一样，我们先计算出矩阵D2表示成如下： 

此时G2表示成如下：

 

由矩阵G2可知A、B药物属于一个类，C、D属于一类；
然后，利用均值等方法再次更新该类的中心值。 

c1=（(1+2)/2,(1+1)/2）=(1.5,1)

c2=（(4+5)/2,(3+4)/2）=(4.5,3.5)

因为对应的中心点并没有发生变化，所以迭代停止，计算完毕。

本算法的时间复杂度：O(tkmn)，其中，t为迭代次数，k为簇的数目，m为记录数，n为维数；

空间复杂度：O((m+k)n)，其中，k为簇的数目，m为记录数，n为维数。

适用范围：

K-menas算法试图找到使平凡误差准则函数最小的簇。当潜在的簇形状是凸面的，簇与簇之间区别较明显，且簇大小相近时，其聚类结果较理想。前面提到，该算法时间复杂度为O(tkmn)，与样本数量线性相关，所以，对于处理大数据集合，该算法非常高效，且伸缩性较好。但该算法除了要事先确定簇数K和对初始聚类中心敏感外，经常以局部最优结束，同时对“噪声”和孤立点敏感，并且该方法不适于发现非凸面形状的簇或大小差别很大的簇。

缺点：

1、聚类中心的个数K 需要事先给定，但在实际中这个 K 值的选定是非常难以估计的，很多时候，事先并不知道给定的数据集应该分成多少个类别才最合适； 
2、Kmeans需要人为地确定初始聚类中心，不同的初始聚类中心可能导致完全不同的聚类结果。（可以使用K-means++算法来解决）
 

Python代码实现

# -*- coding:utf-8 -*-
from numpy import *
def loadDataSet(filename):
    dataMat = []
    fr = open(filename)
    for line in fr.readlines():
        curline = line.strip().split('\t')
        fltline = map(float,curline)
        dataMat.append(fltline)
    return dataMat

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

#随机选择中心点
def randCent(dataSet,k):
    n = shape(dataSet)[1]
    centriods = mat(zeros((k,n)))
    for j in range(n):
        minJ = min(dataSet[:,j])
        rangeJ = float(max(array(dataSet)[:,j])-minJ)
        centriods[:,j] = minJ+rangeJ*random.rand(k,j)
    return centriods

def kmeans(dataSet,k,disMea=distE,createCent=randCent):
    m = shape(dataSet)[0]
    clusterA= mat(zeros((m,2)))
    centriods = createCent(dataSet,k)
    clusterC=True
    while clusterC:
        clusterC=False
        for i in range(m):
            minDist = inf
            minIndex = -1
            for j in range(k):
                distJI = disMea(centriods[j,:],dataSet[i,:])
                if distJI<minDist:
                    minDist=distJI;minIndex=j
            if clusterA[i,0] != minIndex:
                clusterC = True

            clusterA[i,:] = minIndex,minDist**2
        print(centriods)

        for cent in range(k):
            ptsInClust = dataSet[nonzero(clusterA[:, 0].A == cent)[0]]  # get all the point in this cluster
            centriods[cent, :] = mean(ptsInClust, axis=0)  # assign centroid to mean
    return centriods,clusterA

def show(dataSet, k, centriods, clusterA):
     import matplotlib.pyplot as plt
     numSamples, dim = dataSet.shape()
     mark = ['or', 'ob', 'og', 'ok', '^r', '+r', 'sr', 'dr', '<r', 'pr']
     for i in range(numSamples):
         markIndex = int(clusterA[i, 0])
         plt.plot(dataSet[i, 0], dataSet[i, 1], mark[markIndex])

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

def main():
    dataMat = mat(loadDataSet('testSet.txt'))
    myCentroids, clustAssing = kmeans(dataMat, 4)
    print(myCentroids)
    show(dataMat, 4, myCentroids, clustAssing)

if __name__ == '__main__':
    main()