K-均值聚类算法对未标注数据分组（1）

1. K-均值聚类算法

优点：容易实现。
缺点：可能收敛到局部最小值，在大规模数据集上收敛较慢。
适用数据类型：数值型数据。

K-均值算法的工作流程是这样的。首先，随机确定&个初始点作为质心。然后将数据集中的每个点分配到一个簇中，具体来讲，为每个点找距其最近的质心，并将其分配给该质心所对应的簇。这一步完成之后，每个簇的质心更新为该簇所有点的平均值。

上述过程的伪代码表示如下：

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

k均值聚类算法：

# -*- coding: utf-8 -*-
from numpy import *# K-均值聚类支持函数
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 distEclud(vecA, vecB):return sqrt(sum(power(vecA - vecB, 2))) # 为给定数据集构建一个包含k个随机质心的集合,是以每列的形式生成的
def randCent(dataSet, k):n = shape(dataSet)[1]centroids = mat(zeros((k,n))) for j in range(n):  minJ = min(dataSet[:,j])  # 找到每一维的最小rangeJ = float(max(dataSet[:,j]) - minJ) # 每一维的最大和最小值之差centroids[:,j] = mat(minJ + rangeJ * random.rand(k,1)) # 生成随机值#print centroids[:,j]return centroids  # 返回随机质心,是和数据点相同的结构# k--均值聚类算法(计算质心--分配--重新计算)
def kMeans(dataSet, k, distMeas=distEclud, createCent=randCent): # k是簇的数目m = shape(dataSet)[0]  # 得到样本的数目clusterAssment = mat(zeros((m,2))) #  创建矩阵来存储每个点的簇分配结果#  第一列：记录簇索引值，第二列：存储误差，欧式距离的平方centroids = createCent(dataSet, k)  # 创建k个随机质心clusterChanged = Truewhile clusterChanged:  # 迭代使用while循环来实现clusterChanged = False  for i in range(m):  # 遍历每个数据点，找到距离每个点最近的质心minDist = inf; minIndex = -1for j in range(k):  # 寻找最近的质心distJI = distMeas(centroids[j,:],dataSet[i,:])if distJI < minDist:minDist = distJI; minIndex = jif clusterAssment[i,0] != minIndex: # 更新停止的条件clusterChanged = TrueclusterAssment[i,:] = minIndex,minDist**2 # minDist**2就去掉了根号         for cent in range(k):  # 更新质心的位置ptsInClust = dataSet[nonzero(clusterAssment[:,0].A==cent)[0]] centroids[cent,:] = mean(ptsInClust, axis=0) # 然后计算均值，axis=0:沿列方向 print 'centroids:',centroidsreturn centroids, clusterAssment # 返回簇和每个簇的误差值，误差值是当前点到该簇的质心的距离# 主函数
datMat=mat(loadDataSet('testSet.txt'))
rand_centroids=randCent(datMat, 2)
print 'Random centroids:',rand_centroids
myCentroids,clustAssing=kMeans(datMat,4)
print 'myCentroids：',myCentroids
#centList,myNewAssments=biKmeans(datMat, 3)

运行结果：

Random centroids: [[ 0.47004012 -0.70046073][ 4.59776536  2.61143197]]
centroids: [[-3.38237045 -2.9473363 ][-2.46154315  2.78737555][ 2.80293085 -2.7315146 ][ 2.6265299   3.10868015]]
myCentroids： [[-3.38237045 -2.9473363 ][-2.46154315  2.78737555][ 2.80293085 -2.7315146 ][ 2.6265299   3.10868015]]

上面的结果可以看出得到了4个质心

2. 二分k-均值算法

在K-均值聚类中簇的数目k是一个用户预先定义的参数，那么用户如何才能知道乂的选择是否正确？如何才能知道生成的簇比较好呢？在包含簇分配结果的矩阵中保存着每个点的误差，即该点到簇质心的距离平方值。

使用后处理来提高聚类性能：

K-均值算法收敛到了局部最小值，而非全局最小值（局部最小值指结果还可以但并非最好结果，全局最小值是可能的最好结果)。
一种用于度量聚类效果的指标是SSE(Sum of Squared Error,误差平方和。SSE值越小表示数据点越接近于它们的质心，聚类效果也越好。因为对误差取了平方，因此更加重视那些远离中心的点。一种肯定可以降低SSE值的方法是增加簇的个数，但这违背了聚类的目标。聚类的目标是在保持族数目不变的情况下提高簇的质量。那么如何进行改进？可以对生成的簇进行后处理，一种方法是将具有最大SSE值的簇划分成两个簇。具体实现时可以将最大簇包含的点过滤出来并在这些点上运行尺-均值算法，其中的K设为2。

为克服k-均值算法收敛于局部最小值的问题：

二分k-均值（bisecting K-means)算法，该算法首先将所有点作为一个簇，然后将该簇一分为二。之后选择其中一个簇继续进行划分，选择哪一个簇进行划分取决于对”其划分是否可以最大程度降低SSE的值。上述基于SSE的划分过程不断重复，直到得到用户指定的簇数目为止。
另一种做法是选择SSE最大的簇进行划分，直到簇数目达到用户指定的数目为止。这个做法听起来并不难实现。

下面就来看一下该算法的实际效果。

# -*- coding: utf-8 -*-
from numpy import *# K-均值聚类支持函数
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 distEclud(vecA, vecB):return sqrt(sum(power(vecA - vecB, 2))) # 为给定数据集构建一个包含k个随机质心的集合,是以每列的形式生成的
def randCent(dataSet, k):n = shape(dataSet)[1]centroids = mat(zeros((k,n))) for j in range(n):  minJ = min(dataSet[:,j])  # 找到每一维的最小rangeJ = float(max(dataSet[:,j]) - minJ) # 每一维的最大和最小值之差centroids[:,j] = mat(minJ + rangeJ * random.rand(k,1)) # 生成随机值#print centroids[:,j]return centroids  # 返回随机质心,是和数据点相同的结构# k--均值聚类算法(计算质心--分配--重新计算)
def kMeans(dataSet, k, distMeas=distEclud, createCent=randCent): # k是簇的数目m = shape(dataSet)[0]  # 得到样本的数目clusterAssment = mat(zeros((m,2))) #  创建矩阵来存储每个点的簇分配结果#  第一列：记录簇索引值，第二列：存储误差，欧式距离的平方centroids = createCent(dataSet, k)  # 创建k个随机质心clusterChanged = Truewhile clusterChanged:  # 迭代使用while循环来实现clusterChanged = False  for i in range(m):  # 遍历每个数据点，找到距离每个点最近的质心minDist = inf; minIndex = -1for j in range(k):  # 寻找最近的质心distJI = distMeas(centroids[j,:],dataSet[i,:])if distJI < minDist:minDist = distJI; minIndex = jif clusterAssment[i,0] != minIndex: # 更新停止的条件clusterChanged = TrueclusterAssment[i,:] = minIndex,minDist**2 # minDist**2就去掉了根号         for cent in range(k):  # 更新质心的位置ptsInClust = dataSet[nonzero(clusterAssment[:,0].A==cent)[0]] centroids[cent,:] = mean(ptsInClust, axis=0) # 然后计算均值，axis=0:沿列方向 print 'centroids:',centroidsreturn centroids, clusterAssment # 返回簇和每个簇的误差值，误差值是当前点到该簇的质心的距离# 二分k--均值聚类算法
def biKmeans(dataSet, k, distMeas=distEclud):m = shape(dataSet)[0] clusterAssment = mat(zeros((m,2))) # 存储数据集中每个点的簇分配结果及平方误差centroid0 = mean(dataSet, axis=0).tolist()[0] # 计算整个数据集的质心：1*2的向量centList =[centroid0] # []的意思是使用一个列表保存所有的质心,簇列表,[]的作用很大for j in range(m):  # 遍历所有的数据点，计算到初始质心的误差值，存储在第1列clusterAssment[j,1] = distMeas(mat(centroid0), dataSet[j,:])**2while (len(centList) < k):  # 不断对簇进行划分，直到klowestSSE = inf  # 初始化SSE为无穷大for i in range(len(centList)): # 遍历每一个簇print 'i:',i               # 数组过滤得到所有的类别簇等于i的数据集ptsInCurrCluster = dataSet[nonzero(clusterAssment[:,0].A==i)[0],:]# 得到2个簇和每个簇的误差，centroidMat：簇矩阵  splitClustAss：[索引值,误差]centroidMat, splitClustAss = kMeans(ptsInCurrCluster, 2, distMeas) # centroidMat是矩阵sseSplit = sum(splitClustAss[:,1])  # 求二分k划分后所有数据点的误差和     # 数组过滤得到整个数据点集的簇中不等于i的点集print nonzero(clusterAssment[:,0].A!=i)[0]sseNotSplit = sum(clusterAssment[nonzero(clusterAssment[:,0].A!=i)[0],1])# 所有剩余数据集的误差之和print "sseSplit and notSplit: ",sseSplit,',',sseNotSplitif (sseSplit + sseNotSplit) < lowestSSE: # 划分后的误差和小于当前的误差，本次划分被保存print 'here..........'bestCentToSplit = i  # i代表簇数bestNewCents = centroidMat  # 保存簇矩阵print 'bestNewCents',bestNewCentsbestClustAss = splitClustAss.copy() # 拷贝所有数据点的簇索引和误差lowestSSE = sseSplit + sseNotSplit  # 保存当前误差和# centList是原划分的簇向量，bestCentToSplit是i值print 'len(centList) and  bestCentToSplit ',len(centList),',',bestCentToSplit# 数组过滤得到的是新划分的簇类别是1的数据集的类别簇重新划为新的类别值为最大的类别数bestClustAss[nonzero(bestClustAss[:,0].A == 1)[0],0] = len(centList) # 数组过滤得到的是新划分的簇类别是0的数据集的类别簇重新划为新的类别值为ibestClustAss[nonzero(bestClustAss[:,0].A == 0)[0],0] = bestCentToSplitprint 'the bestCentToSplit is: ',bestCentToSplit   # 代表的是划分的簇个数-1print 'the len of bestClustAss is: ', len(bestClustAss) # 数据簇的数据点个数# 新划分簇矩阵的第0簇向量新增到当前的簇列表中centList[bestCentToSplit] = bestNewCents[0,:].tolist()[0] print 'centList[bestCentToSplit]:',centList[bestCentToSplit]# 新划分簇矩阵的第1簇向量添加到当前的簇列表中centList.append(bestNewCents[1,:].tolist()[0]) # centList是列表的格式print 'centList',centList# 数组过滤得到所有数据集中簇类别是新簇的数据点clusterAssment[nonzero(clusterAssment[:,0].A == bestCentToSplit)[0],:]= bestClustAssreturn mat(centList), clusterAssment # 返回质心列表和簇分配结果# 主函数
datMat=mat(loadDataSet('testSet.txt'))
centList,myNewAssments=biKmeans(datMat, 3)

运行结果：

i: 0
centroids: [[-1.12630774  1.13042276][ 2.59258145 -2.78274655]]
[]
sseSplit and notSplit:  1000.74868136 , 0.0
here..........
bestNewCents [[-1.12630774  1.13042276][ 2.59258145 -2.78274655]]
len(centList) and  bestCentToSplit  1 , 0
the bestCentToSplit is:  0
the len of bestClustAss is:  80
centList[bestCentToSplit]: [-1.1263077413793101, 1.13042275862069]
centList [[-1.1263077413793101, 1.13042275862069], [2.592581454545454, -2.782746545454545]]
i: 0
centroids: [[ 1.88660209  3.23434291][-3.10621991 -0.25215334]]
[ 2  6 10 ..., 72 74 78]
sseSplit and notSplit:  390.512447704 , 95.5368926046
here..........
bestNewCents [[ 1.88660209  3.23434291][-3.10621991 -0.25215334]]
i: 1
centroids: [[ 4.5110462  -1.0349174 ][ 2.02832712 -3.29681394]]
[ 0  1  3 ..., 76 77 79]
sseSplit and notSplit:  51.9548039494 , 905.211788759
len(centList) and  bestCentToSplit  2 , 0
the bestCentToSplit is:  0
the len of bestClustAss is:  58
centList[bestCentToSplit]: [1.8866020869565217, 3.2343429130434784]
centList [[1.8866020869565217, 3.2343429130434784], [2.592581454545454, -2.782746545454545], [-3.1062199142857145, -0.2521533428571429]]

可以看出算法运行的中间过程，其中最容易迷惑的地方就是数组过滤器。可参考：
tolist()的用法
nonzero()

要注意的几点：

（1）centList =[centroid0] # []的意思是使用一个列表保存所有的质心,簇列表,[]的作用很大

In [12]: mean(datMat, axis=0).tolist()[0]
Out[12]: [-0.10361321250000004, 0.05430119999999998]In [13]: centroid0=mean(datMat, axis=0).tolist()[0]In [14]: centList =[centroid0]In [15]: centList
Out[15]: [[-0.10361321250000004, 0.05430119999999998]]

（2）dataSet[nonzero(clusterAssment[:,0].A==i)[0],:]

>>> from numpy import *
>>> clusterAssment = mat(zeros((5,2)))
>>> clusterAssment[:,0].A==0
array([[ True],[ True],[ True],[ True],[ True]], dtype=bool)
>>> clusterAssment[:,0]==0
matrix([[ True],[ True],[ True],[ True],[ True]], dtype=bool)
>>> clusterAssment[:,0].A==1
array([[False],[False],[False],[False],[False]], dtype=bool)
>>> clusterAssment
matrix([[ 0.,  0.],[ 0.,  0.],[ 0.,  0.],[ 0.,  0.],[ 0.,  0.]])
>>> mm=transpose(nonzero(clusterAssment[:,0]==0))
>>> mm
array([[0, 0],[1, 0],[2, 0],[3, 0],[4, 0]], dtype=int64)
>>> nonzero(clusterAssment[:,0].A==0)[0]
array([0, 1, 2, 3, 4], dtype=int64)
>>> transpose(nonzero(clusterAssment[:,0].A==0))
array([[0, 0],[1, 0],[2, 0],[3, 0],[4, 0]], dtype=int64)
>>> nonzero(clusterAssment[:,0].A==0)
(array([0, 1, 2, 3, 4], dtype=int64), array([0, 0, 0, 0, 0], dtype=int64))
>>> dataSet=mat([[1,2],[2,3],[3,4],[4,5],[5,6]])
>>> dataSet[nonzero(clusterAssment[:,0].A==0)[0],:]
matrix([[1, 2],[2, 3],[3, 4],[4, 5],[5, 6]])
>>>

（3）bestNewCents[0,:].tolist()[0]

In [17]:a=mat([[1,2],[2,3]])In [18]:a
Out[18]:
matrix([[1, 2],[2, 3]])In [19]:a[0,:]
Out[19]: matrix([[1, 2]])In [20]:a[0,:].tolist()
Out[20]: [[1, 2]]In [21]:a[0,:].tolist()[0]
Out[21]: [1, 2]