01 物以类聚
经过半年的不懈努力，我们已经学习并实践了经典的分类算法和经典的回归算法，下面我们开始学习经典的聚类算法(兴奋~~~)
目前打算对三种聚类算法进行学习和代码实操(俗称“造轮子”)：
KMeans
Apriori
FP-Growth
今天我们学习并实践KMeans聚类算法，分成以下几个部分，跟上节奏燥起来！
KMeans算法理论和代码实现
改进，BiKMeans算法理论和代码实现
实例，上车点规划
抉择，如何挑选最佳的聚类簇数？
其中一二三部分参考了Peter Harrington的《机器学习实战》，第四部分是纯手工打造的轮子，还请多多指教。
02 KMeans理论和算法实现
聚类是一种无监督学习的方法，所谓“无监督”，就是指参与训练的样本没有标签。
KMeans聚类算法过程如下：
1. 对于一组数据集，随机选取k个点作为质心，将数据集中的点归为离其最近的质心一簇，此时数据集被划分为k个簇；
2. 对这k个簇，重新计算各簇的质心(均值)；
3. 根据新的质心，按照step1继续聚类，然后再根据聚类重新计算各簇质心，直到质心不再改变，分类完成。
说白了，就是不断地聚类、划分的过程。
通过KMeans原理，可以看到几个显而易见的缺点：
1. 簇数量k由用户指定，无法预先知道最佳k值 >>解法：分为几簇，最终由轮廓系数S(i)决定，取轮廓系数最大的分类数(05节剧透)
2. 最终质心可能与初始点选择有关 >> 因此KMeans的结果可能收敛到局部最小值，而不是全局最小值 >> 解法：
BiKMeans（03节）
KMeans++（KMeans++ 算法在选择初始质心时并不是随机选择，而是选择尽量相互分离的质心，即，下一个质心点总是离上一个质心点较远的点）
代码实现
def loadDataSet(fileName):
dataList=[]
dataMat=[]
fr=open(fileName) for line in fr.readlines():
curLine=line.strip().split('\t')
fltLine=list(map(float,curLine))
dataList.append(fltLine)
dataMat=mat(dataList) return dataMatdef distEclud(vecA,vecB):
return sqrt(sum(power(vecA-vecB,2))) #欧式距离#为输入数据集构造k个随机中心，中心位置在各特征最大最小值之间def randCent(dataSet,k):
n=shape(dataSet)[1]
center=mat(zeros((k,n))) for j in range(n): #对每个特征
minJ=min(dataSet[:,j])
rangeJ=float(max(dataSet[:,j])-minJ)
center[:,j]=mat(minJ+rangeJ*random.rand(k,1)) #质心第j维坐标在数据集第j维数据之间
return centerdef KMeans(dataSet,k,distMeas=distEclud,createCent=randCent):
m=shape(dataSet)[0]
clusterAssment=mat(zeros((m,2))) #用于记录各样本当前归属于哪个簇以及到该簇质心的欧式距离平方
center=createCent(dataSet,k)
clusterChanged=True
while clusterChanged:
clusterChanged=False
#对每个样本，计算样本到各质心的距离，寻找距离最近的质心，将该样本归为该质心所在簇
for i in range(m):
minDist=inf;minIndex=-1
for j in range(k): #对每个质心,计算到i样本的距离
distJI=distMeas(center[j,:],dataSet[i,:]) if distJI<minDist:
minDist=distJI;minIndex=j #i样本暂属于j簇，到j簇质心距离为minDist
if clusterAssment[i,0]!=minIndex:
clusterChanged=True #若任一样本在本次迭代中改变了簇类，则要进行下一次迭代(即，直到任何样本都不再改变簇类，聚类停止)
clusterAssment[i,:]=minIndex,minDist**2 #记录样本i的簇类情况
#print (center)
#更新质心
for cent in range(k):
ptsInClust=dataSet[nonzero(clusterAssment[:,0].A==cent)[0]] #筛选出属于当前簇类的点
center[cent,:]=mean(ptsInClust,axis=0) #对该簇类各样本的各列求均值，作为新质心
return center,clusterAssment用一组数据来测试一下：
dataMat1=loadDataSet(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\testSet.txt') dataMat2=loadDataSet(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\testSet2.txt') center_testSet1,clusterAssment_testSet1=KMeans(dataMat1,4) center_testSet2,clusterAssment_testSet2=KMeans(dataMat2,3) plt.figure(figsize=(6,6)) plt.scatter(dataMat1[:,0].T.tolist()[0],dataMat1[:,1].T.tolist()[0],c='pink',s=30) plt.scatter(center_testSet1.T[0].tolist()[0],center_testSet1.T[1].tolist()[0],c='blue',s=50)
结果如下，蓝色点为聚类质心
03 BiKMeans理论和算法实现
刚才我们实现了KMeans算法，不过也提到，KMeans有个缺点，就是其聚类的最终质心可能与初始点选择有关，因此KMeans的结果可能收敛到局部最小值，而不是全局最小值。
怎么解决呢？提示：想想决策树是如何分支的？
为了得到全局最优解，BiKMeans算法出现了，它的过程如下（是不是跟决策树分支有点神似？）
1. 将整个数据集看作一个簇，计算初始质心，即所有数据点各特征的均值
2. 遍历各质心，对各质心，将质心所在簇用原始KMeans算法二分，计算二分后整个数据集的SSE(即平方误差和，即簇各点到簇质心距离平方和)，找到二分后整体数据集SSE最小的质心，认为此质心是本次划分的最佳质心，对其进行二分
3. 不断重复step2，直到质心总数=设置的k
说白了，BiKMeans算法过程类似于决策树的分支，通过启发的方法，每次迭代只分裂当前最佳质心，直到簇数量达到k，这样的方法可以保证最终划分得到的质心是全局最优解，而原始KMeans可能会陷入局部最优解。
代码实现
def biKMeans(dataSet,k,distMeas=distEclud):
m=shape(dataSet)[0]
clusterAssment=mat(zeros((m,2))) #记录各样本归属和距离平方
center0=mean(dataSet,axis=0).tolist()[0] #初始质心
centerList=[center0] #用于记录聚类质心坐标
for j in range(m):
clusterAssment[j,1]=distMeas(mat(center0),dataSet[j,:])**2
while len(centerList)<k:
lowestSSE=inf #SSE=Sum of Square Error
#对每个簇，尝试二分，计算二分后整体数据的SSE，若小于lowestSSE，则将簇二分，如此往复，直到分到k个簇为止
for i in range(len(centerList)):
ptsInCluster=dataSet[nonzero(clusterAssment[:,0].A==i)[0],:] #默认质心索引就是数据对应簇类
centerMat,splitClustAss=KMeans(ptsInCluster,2,distMeas)
sseSplit=sum(splitClustAss[:,1]) #被二分后的簇的平方误差和
sseNotSplit=sum(clusterAssment[nonzero(clusterAssment[:,0]!=i)[0],1]) #整体数据中，未被二分的簇的平方误差和
if (sseSplit+sseNotSplit)<lowestSSE:
bestCentToSplit=i
bestNewCents=centerMat
bestClustAss=splitClustAss.copy() #.copy()防止splitClustAss被覆盖时影响到bestClustAss
lowestSSE=sseSplit+sseNotSplit #确定好本次迭代被二分的簇后，将被二分的数据对应的簇类更新
bestClustAss[nonzero(bestClustAss[:,0].A!=0)[0],0]=len(centerList) #更新先后顺序很重要！
bestClustAss[nonzero(bestClustAss[:,0].A==0)[0],0]=bestCentToSplit
print('The bestCentToSplit is:',bestCentToSplit)
print('The number of samples to split is',len(bestClustAss))
centerList[bestCentToSplit]=bestNewCents[0,:].tolist()[0] #将被二分的原簇质心替换为二分后的质心之一
centerList.append(bestNewCents[1,:].tolist()[0]) #将二分后的另一质心添加在质心列表末尾
#将二分后的簇的数据归属更新到总记录中
clusterAssment[nonzero(clusterAssment[:,0].A==bestCentToSplit)[0],:]=bestClustAss return mat(centerList),clusterAssment测试
center_testSet22,clusterAssment_testSet22=biKMeans(dataMat2,3) plt.figure(figsize=(6,6)) plt.scatter(dataMat2[:,0].T.tolist()[0],dataMat2[:,1].T.tolist()[0],c='pink',s=30) plt.scatter(center_testSet22.T[0].tolist()[0],center_testSet22.T[1].tolist()[0],c='blue',s=50)
04 实例：上车点规划
我们在前节实现了KMeans算法和BiKMeans算法，现在让我们来用用这两个轮子吧。
设想，你邀请了70个朋友参加你的party，他们住在城市的各个地方，你需要在party开始前2小时开车去接他们(假设你的车子可以容纳70人，哈哈哈)。
为了时间效率最大化，请问你该如何规划每个人的上车点呢？总不能你一个一个去家门口接吧。
我们用聚类的方法来规划！
现在我们获取了每个朋友住址的经纬度，那么请开始你的表演。
很简单，只需要调用我们早好的轮子BiKMeans
#球面距离计算函数：球面余弦距离(向量夹角*地球半径)#求球面上两向量vecA,vecB的距离def distSLC(vecA,vecB):
a=sin(vecA[0,1]*pi/180)*sin(vecB[0,1]*pi/180) #由于抓取的经纬度为角度，需要通过 *pi/180来转换为弧度
b=cos(vecA[0,1]*pi/180)*cos(vecB[0,1]*pi/180)*cos(pi*(vecA[0,0]-vecB[0,0])/180) return 6371.0*arccos(a+b) #6371为地球半径，单位为英尺def clusterClubs(k=5):
dataList=[] for line in open(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\places.txt').readlines():
lineArr=line.strip().split('\t')
dataList.append([float(lineArr[4]),float(lineArr[3])]) #读取地址的经纬度
dataMat=mat(dataList)
center,clusterAss=biKMeans(dataMat,k,distMeas=distSLC) #将地址按经纬度聚类
#作图
fig=plt.figure(figsize=(10,8))
rect=[0.1,0.1,0.8,0.8] #用于设置坐标轴刻度，[]中前两个值表示左边起始位置，后两个值对应坐标长度
scatterMarkers=['^','o','h','8','p','d','v','s','>','<'] #用于设置散点图点的形状
axprops=dict(xticks=[],yticks=[])
ax0=fig.add_axes(rect,label='ax0',**axprops)
#读图，并将图片显示在设定好的坐标轴中
imgP=plt.imread(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\Portland.png')
ax0.imshow(imgP) #将地址按聚类结果作散点图
ax1=fig.add_axes(rect,label='ax1',frameon=False) for i in range(k):
ptsInCluster=dataMat[nonzero(clusterAss[:,0].A==i)[0],:]
markerStyle=scatterMarkers[i]
ax1.scatter(ptsInCluster[:,0].flatten().A[0],ptsInCluster[:,1].flatten().A[0],\
marker=markerStyle,s=90) #标出质心
ax1.scatter(center[:,0].flatten().A[0],center[:,1].flatten().A[0],marker='+',s=300)
plt.show()好了，我们写的这个“上车点规划器”就可以使用了，只需要输入聚类簇数量即可 ，下面看看几个测试结果：
聚为5类（5个上车点）
聚为9类（9个上车点）
但是到底聚为几类比较合理呢？这个合理是指即不让你的朋友跑太远去坐车，也不用你跑太多地方去接人。
05 如何挑选最佳的聚类簇数
不论是KMeans算法还是BiKMeans算法，都仍有一个缺点没有解决：簇数量k由用户指定，用户指定的k不一定是最佳的簇数量。这也是上一节我们没有解决的问题。
怎么办呢？
使用轮廓系数，最佳k值由轮廓系数S(i)决定，取轮廓系数最大的分类数。
什么是轮廓系数，轮廓系数S(i)用于衡量聚类效果，取值在-1~1之间，越接近1表示聚类效果越好，公式如下，
其中a(i)=i点到同簇各点距离均值；b(i)=min(i点到某个非同簇各点均值)，即i点到其他簇质心距离的最小值，整体数据集的轮廓系数是各点轮廓系数的平均值
基于此，我们可以写一个计算聚类结果轮廓系数的函数，然后看看轮廓系数与聚类数量k的关系，从而找到最佳的聚类数量k。
轮廓系数计算函数
def outlineOfCluster(filename,maxClusterNum,distMeas=distEclud):
dataList=[] for line in open(filename).readlines():
lineArr=line.strip().split('\t')
dataList.append([float(lineArr[4]),float(lineArr[3])]) #读取地址的经纬度
dataMat=mat(dataList) #遍历2~nn个质心，求不同数量的簇的轮廓系数，设置轮廓系数最大的簇数量为k值
sk={} for k in range(2,maxClusterNum+1):
center,clusterAss=biKMeans(dataMat,k,distMeas)
outline=[] #代表聚类为k个簇时各点的轮廓系数列表
for i in range(k): #遍历各簇i
ptsInCluster=clusterAss[nonzero(clusterAss[:,0].A==i)[0],:] for j in range(len(ptsInCluster)): #遍历i簇的各点j
ptsInClusterNotJ=vstack([ptsInCluster[:j,:],ptsInCluster[j+1:,:]])
ajn=[] for n in range(len(ptsInClusterNotJ)): #遍历i簇非j的点
ajn.append(distSLC(ptsInCluster[n],ptsInCluster[j]))
aj=nanmean(ajn)
bjm=[]
centerWithoutI=vstack([center[:i,:],center[i+1:,:]]) for m in range(len(centerWithoutI)): #遍历非i簇质心
bjm.append(distSLC(centerWithoutI[m],ptsInCluster[j]))
bj=min(bjm)
sj=(bj-aj)/max(bj,aj) #i簇中j点的轮廓系数
outline.append(sj) #将i簇中各点的轮廓系数保存在outline中，outline用于存储聚类为k个簇时各点的轮廓系数，在遍历k时需要重置
sk[k]=nanmean(outline) #聚类为k个簇时的轮廓系数是各点轮廓系数的均值
return sk现在我们可以接近上车点规划问题的遗留问题了，规划几个上车点最合理？
sk=outlineOfCluster(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\places.txt',35,distMeas=distSLC)
plt.figure(figsize=(10,5))
plt.plot(list(sk.keys()),list(sk.values()),linewidth=3)
plt.title('聚类簇数-效果',fontsize=18,fontweight='bold')
plt.xlabel('最终聚类簇数量',fontsize=14)
plt.ylabel('轮廓系数',fontsize=14)
可以看到，聚类簇数在达到10个簇之后，轮廓系数的增量就变得非常小了，因此从性价比方面考虑，选择聚类为10个簇的性价比是最高的，即规划10个上车点就好。
06 总结
本文实践了KMeans聚类算法，分成以下几个部分，
KMeans算法理论和代码实现
改进，BiKMeans算法理论和代码实现
实例，上车点规划
抉择，如何挑选最佳的聚类簇数？
下期我们将实践Apriori算法，这是一种挖掘关联规则得到频繁项集的算法，比如可以预测你买了A商品后会买B商品，敬请期待~
07 参考
《机器学习实战》 Peter Harrington Chapter10
作者：邓莎
链接：https://www.jianshu.com/p/605ff743b83b
共同学习，写下你的评论
评论加载中...
作者其他优质文章
