MLiA笔记_kMeans

# -*- coding:utf-8 -*-

# 10.1 K-均值聚类支持函数
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 distEclud(vecA, vecB):
    return sqrt(sum(power(vecA - vecB, 2)))

# randCent()函数，为给定数据集构建一个包含k个随机质心的集合
def randCent(dataSet, k):
    n = shape(dataSet)[1]
    centroids = mat(zeros((k,n)))
    # 随机质心必须要在整个数据集的边界之内，可通过找到数据集每一维德最小和最大值来完成，然后生成(0,1)之间的随机数并通过取值范围和最小值
    for j in range(n):
        minJ = min(dataSet[:,j])
        rangeJ = float(max(dataSet[:,j]) - minJ)
        centroids[:,j ] = minJ + rangeJ*random.rand(k,1)
    return centroids

# 10.2 K-均值聚类算法
# kMeans()，数据集和簇的数目是必选参数，计算距离和创建初始质心的函数是可选的
def kMeans(dataSet, k, disMeas=distEclud, createCent = randCent):
    # 一开始确定数据集中数据点的总数
    m = shape(dataSet)[0]
    # 创建一个矩阵来存储每个点的簇分配结果。
    # 簇分配结果矩阵包含两列：一列记录索引值，一列记录存储误差
    clusterAssment = mat(zeros((m,2)))
    # 误差指当前点到簇质心的距离
    centroids = createCent(dataSet,k)
    clusterChanged = True
    # while循环反复迭代
    while clusterChanged:
        clusterChanged = False
        for i in range(m):
            minDist = inf
            minIndex = -1
            # 遍历所有数据找到距离每个点最近的质心
            for j in range(k):
                # 计算距离是使用disMeas参数给出的距离函数，默认距离函数是disEclud()
                distJI = disMeas(centroids[j,:],dataSet[i,:])
                if distJI < minDist:
                    minDist = distJI
                    minIndex = j
            # 如果任一点的簇分配结果发生改变，则更新clusterChanged标志
            if clusterAssment[i,0] != minIndex:
                clusterChamged = True
            clusterAssment[i,:] = minIndex, minDist**2
        print centroids
        # 遍历所有质心并更新它们的取值
        for cent in range(k):
            # 首先通过数组过滤来获得给定簇的所有点，
            ptsInclust = dataSet[nonzero(clusterAssment[:,0].A == cent)[0]]
            # 然后计算所有点的均值，axis=0表示沿矩阵列方向进行均值计算
            centroids[cent,:] = mean(ptsInclust, axis=0)
    # 最后程序返回所有的类质心与点分配结果
    return centroids, clusterAssment

# 10.3 二分K-均值算法
# 与10-2中kMeans()的参数相同，在给定数据集、所期望的簇数目和距离计算方法的条件下，函数返回聚类结果
def biKmenas(dataSet, k, disMeas=distEclud):
    m = shape(dataSet)[0]
    # 首先创建矩阵来存储数据集中的每个点的簇分配结果及平方误差
    clusterAssment = mat(zeros((m,2)))
    # 然后计算整个数据集的质心，
    centroid0 = mean(dataSet, axis=0).tolist()[0]
    # 并使用一个列表来保留所有的质心
    centList = [centroid0]
    for j in range(m):
        clusterAssment[j,1] = disMeas(mat(centroid0),dataSet[j,:])**2
    # while循环不停地对簇进行划分，直到得到想要的簇数目位置，为此需要比较划分前后的SSE
    while (len(centList) < k):
        # 一开始将SSE设为无穷大
        lowestSSE = inf
        # 然后遍历所有的簇来决定最佳的簇进行划分
        for i in range(len(centList)):
            # 对每个簇，将该簇中的所有点看成一个小的数据集ptsInCurrCluster
            ptsInCurrCluster = dataSet[nonzero(clusterAssment[:,0].A == i)[0],:]
            # 将ptsInCurrCluster输入到kMeans中进行处理，K-均值算法会生成两个质心簇，同时给出每个簇的误差值
            centroidMat, splitClustAss = kMeans(ptsInCurrCluster, 2, disMeas)
            # 这些误差与剩余数据集的误差之和作为本次划分的误差，
            sseSplit = sum(splitClustAss[:,1])
            sseNotSplit = sum(clusterAssment[nonzero(clusterAssment[:,0].A != i)[0],1])
            print "sseSplit, and notSplit: ",sseSplit, sseNotSplit
            # 如果该划分的SSE值最下，则本次划分被保存
            if (sseSplit + sseNotSplit) < lowestSSE:
                bestCentToSplit = i
                bestNewCents = centroidMat
                bestClustAss = splitClustAss.copy()
                lowestSSE = sseSplit + sseNotSplit
        # 一旦决定了要划分的簇，接下来就要实际执行划分操作。
        # 划分操作：将要划分的簇中所有点的簇分配结果进行修改
        # 当使用kMeans函数并指定簇数为2时，会得到两个编号分别为0和1的结果簇。
        # 需要将这些簇编号修改为划分簇和新加簇的编号，该过程通过两个数组过滤器来完成
        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中
        centList[bestCentToSplit] = bestNewCents[0,:]
        centList.append(bestNewCents[1,:])
        clusterAssment[nonzero(clusterAssment[:,0].A == bestCentToSplit)[0],:] = bestClustAss
    # 当while循环结束，同kMeans函数一样，函数返回质心列表和簇分配结果
    return mat(centList), clusterAssment

# yahoo！PlaceFinder API
import urllib
import json

# 函数geoGrab()从yahoo返回一个字典
def geoGrab(stAddress, city):
    # 首先为yahooApi设置apiStem
    apiStem = 'http://where.yahooapis.com/geocode?'
    # 然后创建一个字典
    params = {}
    # 为字典设置不同值，以便返回JSON格式的结果
    params['flags'] = 'J'
    params['appid'] = 'ppp68N8t'
    params['location'] = '%s %s' % (stAddress. city)
    # 接下来使用urllib的urlencode()将创建的字典转换为可以通过URL进行传递的字符串格式
    url_params = urllib.urlencode(params)
    # 最后打开URL读取返回值
    yahooApi = apiStem + url_params
    print yahooApi
    c = urllib.urlopen(yahooApi)
    return json.loads(c.read())

from time import sleep

# massPlaceFind()将所有这些封装起来并且将相关信息保存到文件中
def massPlaceFind(fileName):
    # 该函数打开一个tab分隔de文本文件，
    fw =  open('places.txt','w')
    for line in open(fileName).readlines():
        line = line.strip()
        lineArr = line.split('\t')
        # 获取第2列和第3咧结果，这些值被输入到函数geoGrab()中
        retDict = geoGrab(lineArr[1],lineArr[2])
        # 然后检查geoGrab()的输出字典判断有没有错误，如果有错误就不需要去抽取纬度和经度，如果没有错误就可以读取
        if retDict['ResultSet']['Error'] == 0:
            lat = float(retDict['ResultSet']['Results'][0]['latitude'])
            lng = float(resDict['ResultSet']['Results'][0]['longtitude'])
            print "%S\t$s\t%f" % (lineArr[0].lat,lng)
            # 这些值被添加到原来对应的行上，同时写道一个新的文件中
            fw.write('%s\t$f\t$f\n' % (line,lat,lng))
        else: print("error fetching")
        # 最后调用sleep()函数将massPalceFind()函数延迟1秒，确保不要再短时间内过于频繁地调用API
        sleep(1)
    fw.close()


# 10.5 球面距离计算及簇绘图函数
# distSLC()返回地球表面两点之间的距离
# 给定两个点的经纬度，可以使用球面余弦定理来计算两点的距离
def distSLC(vecA, vecB):
    a = sin(vecA[0,1]*pi/180) * sin(vecB[0,1]*pi/180)
    b = cos(vecA[0,1]*pi/180) * cos(vecB[0,1]*pi/180) * cos(pi*  (vecB[0.0]-vecA[0,0])/180)
    return arccos(a+b)*6371.0

import matplotlib
import matplotlib.pyplot as plt
# clusterClubs()将文本文件中的俱乐部进行聚类并画出结果，只有一个参数，即所希望得到的簇数目
# 该函数将文本文件的解析、聚类以及画图都封装在一起
def clusterClubs(numClust=5):
    # 首先创建一个空列表
    datList = []
    # 打开文件获取第4列he第5列，这两列分别对应纬度和经度
    for line in open('places.txt').readlines():
        lineArr = line.split('\t')
        datList.append([float(lineArr[4]),float(lineArr[3])])
    # 基于这些经纬度对的列表创建一个矩阵
    datMat = mat(datList)
    # 在这些数据点上运行biKmeans()，并使用distSLC()作为聚类中使用的距离计算方法
    myCentroids, clustAssing = biKmenas(datMat, numClust, disMeas=distSLC)
    # 最后将簇以及簇质心画在图上
    fig = plt.figure()
    # 为了画出这些簇，首先创建一幅图和一个矩形
    rect = [0.1,0.1,0.8,0.8]
    scatterMarkers = ['s','c','^','8','p','d','v','h','>','<']
    # 然后使用该矩形来决定绘制图的那一部分
    # 构建一个标记形状的列表用于绘制散点图
    axprops = dict(xticks=[], yticks=[])
    ax0 = fig.add_axes(rect,label='ax0',**axprops)
    # 下一步使用imread()函数基于一幅图像来创建矩阵
    imgP = plt.imread('Portlan.ong')
    # 然后使用imshow()绘制该矩形
    ax0.imshow(imgP)
    ax1 = fig.add_axes(rect, label='ax1',frameon = False)
    for i in range(numClust):
        ptsIncurrCluster = datMat[nonzero(clustAssing[:,0].A==i)[0],:]
        markerStyle = scatterMarkers[i % len(scatterMarkers)]
        ax1.scatter(ptsIncurrCluster[:,0].flatten().A[0],ptsIncurrCluster[:,1].flatten().A[0],marker= markerStyle, s=90)
    ax1.scatter(myCentroids[:,0].flatten().A[0], myCentroids[:,1].flatten().A[0], marker='+', s=300)
    plt.show()