KMeans算法实践

　　这几天学习了无监督学习聚类算法Kmeans，这是聚类中非常简单的一个算法，它的算法思想与监督学习算法KNN（K近邻算法）的理论基础一样都是利用了节点之间的距离度量，不同之处在于KNN是利用了有标签的数据进行分类，而Kmeans则是将无标签的数据聚簇成为一类。接下来主要是我对《机器学习实战》算法示例的代码实现和理解。

     首先叙述下算法项目《对地图上的俱乐部进行聚类》的要求：朋友Drew希望让我们带她去城里庆祝生日，由于其他一些朋友也会过来，所以需要提供一个大家都可行的计划，Drew给出了她希望去的地址列表，有69个位置。这些位置都是俄勒冈州的波特兰地区。

     在给出的地址列表中有69行数据，每行数据分别是名称、具体地址和地区。这么多地址如果要在一晚上全部到达的话，需要合理的安排，在示例中的思想值得学习，将这些位置点在地图上标记出来，通过聚类算法找到数据点的若干个质心，安排交通工具到达这些簇的质心，然后步行到达每个簇内的地址。

1     数据预处理

　　这里只给出了想去位置的具体信息，但是这些数据无法直接用于计算两点之间的距离，因此在这里首先要获取到这些位置的经纬度信息，然后再通过经纬度计算出距离。原文中获取经纬度的方法是通过Yahoo! PlaceFinder API收集数据，但是在学习的过程中发现直接利用基于Python的地理编码库geopy更加的方便，而且该API提供了直接利用经纬度计算距离的方法，非常方便的获得两个位置之间的球面距离。进入到Python的根目录，在命令提示符下输入 pip install geopy安装geopy包。在geopy包下利用地理编码可以获得具体的地址和地区经纬度，同时也可以利用反地理编码将经纬度信息转换成具体位置信息。在这里我定义了一个函数用于获取location

1 def getPlaceThroughGeo(address,city):
2     geolocator = Nominatim()
3     location = geolocator.geocode(address,city)
4     return location

同时定义一个批量处理地址列表portlandClubs.txt的函数

 1 def getPlaceLocation(path):
 2     fw = open('E:/place.txt','w')    #将获取到的经纬度信息写入到place.txt
 3     for line in open(path).readlines():    #遍历解析每一个具体位置
 4         line = line.strip()
 5         lineArr = line.split('\t')
 6         loc = getPlaceThroughGeo(lineArr[1],lineArr[2])   #调用函数得到location
 7         if loc == None:                                   #如果没有解析出位置，loc将为None
 8             print 'Not Fetch the address!'
 9         else:                                             #解析出后，写入到文件中
10             fw.write("%s\t%f\t%f\n" % (lineArr[0],loc[0].latitude,loc[0].longitude))
11     fw.close()

在处理这个函数的时候，由于网络的问题，频繁的调用API可能会出现响应时间延迟的错误信息，我使用了sleep(1、10)仍然会有出错误信息，我就10个为一组，依次运行，最后得到了所有位置的location，显示其中的一部分

2  Kmeans聚类算法

　　这里直接使用二分聚类，这是对于直接聚类算法的改进，因为直接聚类算法在最初的时候随机选择k个质心点，然后做聚类，这种方式首先会受到随机选择的质心点的影响，其次比较容易到达局部最小值而不是全局最小值。因此在此基础上就使用了二分聚类：首先将所有节点作为一个簇，然后将其一分为二，然后选择能够最大程度降低SSE(误差平方和)的簇继续划分，直到簇的个数达到设定的K值。首先根据算法原理写出普通kmeans算法，在此基础上由biKmeans直接调用就好。

def kMeans(dataSet, k, distMeas=getPlaceThroughGeo, createCent=randCent):
    m = shape(dataSet)[0]                     #获取数据集的样本数
    clusterAssment = mat(zeros((m,2)))        #用于存储样本的聚类结果和距离误差
    centroids = createCent(dataSet, k)        #在dataSet中创建随机创建k个质心
    clusterChanged = True                     # 观察标记
    while clusterChanged:
        clusterChanged = False                #如果簇质心不在更改就为false,否则接下来会设为True
        for i in range(m):                    #遍历每个样本
            minDist = inf; minIndex = -1
            for j in range(k):                #将每个样本点划分至最近的质心节点所在簇
                distJI = distMeas(centroids[j,:],dataSet[i,:])
                if distJI < minDist:
                    minDist = distJI; minIndex = j
            if clusterAssment[i,0] != minIndex: clusterChanged = True
            clusterAssment[i,:] = minIndex,minDist**2
        print centroids
        for cent in range(k):                    #更新质心节点，不断迭代
            ptsInClust = dataSet[nonzero(clusterAssment[:,0].A==cent)[0]]#get all the point in this cluster
            centroids[cent,:] = mean(ptsInClust, axis=0) #assign centroid to mean 
    return centroids, clusterAssment

也可以用距离可以直接利用数学球面余弦定理计算得到：
　　R•arc cos[cosβ1cosβ2cos（α1-α2）+sinβ1sinβ2]    其中R为球体半径，此处为地球的半径6371km

二分Kmeans:

def biKmeans(dataSet, k, distMeas=getPlaceThroughGeo):
    m = shape(dataSet)[0]
    clusterAssment = mat(zeros((m,2)))
    centroid0 = mean(dataSet, axis=0).tolist()[0]                        #原始质心，将所有节点作为同一个簇
    centList =[centroid0]                                                #创建质心列表，其长度可以表示簇的个数
    for j in range(m):#calc initial Error
        clusterAssment[j,1] = distMeas(mat(centroid0), dataSet[j,:])**2  #每个节点距离初始质心的距离误差平方
    while (len(centList) < k):
        lowestSSE = inf
        for i in range(len(centList)):
            ptsInCurrCluster = dataSet[nonzero(clusterAssment[:,0].A==i)[0],:]#得到所在簇的所有节点
            centroidMat, splitClustAss = kMeans(ptsInCurrCluster, 2, distMeas) #对数据集进行二分类

            #如何选择要二分的类别，通过遍历每一个类别，计算该类别下进行二分的结果误差之和是否是当前
            #最小的误差，如果是就将该类别进行二分

            sseSplit = sum(splitClustAss[:,1])#compare the SSE to the currrent minimum
            sseNotSplit = sum(clusterAssment[nonzero(clusterAssment[:,0].A!=i)[0],1])
            print "sseSplit, and notSplit: ",sseSplit,sseNotSplit
            if (sseSplit + sseNotSplit) < lowestSSE:
                bestCentToSplit = i
                bestNewCents = centroidMat
                bestClustAss = splitClustAss.copy()
                lowestSSE = sseSplit + sseNotSplit
        bestClustAss[nonzero(bestClustAss[:,0].A == 1)[0],0] = len(centList)      #将二分的类别更新簇号，一个为新分配的类别标号
        bestClustAss[nonzero(bestClustAss[:,0].A == 0)[0],0] = bestCentToSplit   # 一个原标号
        print 'the bestCentToSplit is: ',bestCentToSplit
        print 'the len of bestClustAss is: ', len(bestClustAss)
        centList[bestCentToSplit] = bestNewCents[0,:].tolist()[0]          #更新簇质心的坐标，即更新经纬度
        centList.append(bestNewCents[1,:].tolist()[0])
        clusterAssment[nonzero(clusterAssment[:,0].A == bestCentToSplit)[0],:]= bestClustAss#reassign new clusters, and SSE
    return mat(centList), clusterAssment

3.绘图

根据matplotlib.pyplot可以在图片上绘制出每个点的位置，并将质心点画出。最后运行结果为：

这是给定的标准簇个数为k = 5时的运行结果图，从图中可以看出聚类效果还不错，我又分别在k=4,6,8下测试了算法，结果如下

k=4时波特兰地区夜店分布图

k=6

 

 

 k=8

在k=5时分割之后的误差和与不分割的误差平方和为    sseSplit, and notSplit:  58125.7363242     102958044.466

k=6时   sseSplit, and notSplit:  13667.7740819     41266437.6316

因此最后我觉得在k = 6时结果也许会比较好。

　　这是我的第一篇博客，因为刚刚开始工作主要做关于机器学习的算法，所以现在算法基础还比较薄弱，因此开通这个博客将我学习的算法进行总结，一方面方便以后的复习，另一方面也可以与各位大牛进行沟通探讨提升自己。目前刚开始学习机器学习，所以博客内容会以机器学习为主，同时会写一点其他算法的理解，以后会逐渐写关于深度学习和神经网络相关的学习笔记。还希望能得到各位朋友的建议和批评，以此不断进步！