三种聚类方法的简单实现

聚类是机器学习中的无监督学习方法的重要一种，近来看了周志华老师的机器学习，专门研究了有关于聚类的一章，收获很多，对于其中的算法也动手实现了一下。主要实现的包括比较常见的k均值聚类、密度聚类和层次聚类，这三种聚类方法上原理都不难，算法过程也很清晰明白。有关于原理可以参阅周志华老师的机器学习第九章，这里只做一下代码的实现。

运行环境是Python2.7+numpy，说实话，numpy坑还是挺多的，其实用Matlab可能会更简单。

k均值聚类，核心是是不断更新簇样本的质心。

#encoding=utf-8__author__ = 'freedom'from numpy import*import matplotlib.pyplot as pltdef loadDataSet(fileName):    '''    本函数用于加载数据    :param fileName: 数据文件名    :return:数据集，具有矩阵形式    '''    fr = open(fileName)    dataSet = []    for line in fr.readlines():        curLine = line.strip().split('\t')        inLine = map(float,curLine) # 利用map广播，是的读入的字符串变为浮点型        dataSet.append(inLine)    return mat(dataSet)def getDistance(vecA,vecB):    '''    本函数用于计算欧氏距离    :param vecA: 向量A    :param vecB: 向量B    :return:欧氏距离    '''    return sqrt(sum(power(vecA-vecB,2)))def randCent(dataSet,k):    '''    本函数用于生成k个随机质心    :param dataSet: 数据集，具有矩阵形式    :param k:指定的质心个数    :return:随机质心，具有矩阵形式    '''    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] = minJ + rangeJ*random.rand(k,1) # numpy下的rand表示随机生成k*1的随机数矩阵，范围0-1    return centRoidsdef kMeans(dataSet,k,disMens = getDistance,createCent = randCent):    '''    本函数用于k均值聚类    :param dataSet: 数据集，要求有矩阵形式    :param k: 指定聚类的个数    :param disMens: 求解距离的方式，除欧式距离还可以定义其他距离计算方式    :param createCent: 生成随机质心方式    :return:随机质心，簇索引和误差距离矩阵    '''    m = shape(dataSet)[0]    clusterAssment = mat(zeros((m,2))) # 要为每个样本建立一个簇索引和相对的误差，所以需要m行的矩阵，m就是样本数    centRoids = createCent(dataSet,k) # 生成随机质心    clusterChanged = True    while clusterChanged:        clusterChanged = False        for i in range(m): # 遍历所有样本            minDist = inf;minIndex = -1 # 初始化最小值            for j in range(k): # 遍历所有质心                disJI = disMens(centRoids[j,:],dataSet[i,:])                if disJI < minDist:                    minDist = disJI;minIndex = j # 找出距离当前样本最近的那个质心            if clusterAssment[i,0] != minIndex: # 更新当前样本点所属于的质心                clusterChanged = True # 如果当前样本点不属于当前与之距离最小的质心，则说明簇分配结果仍需要改变                clusterAssment[i,:] = minIndex,minDist**2        for cent in range(k):            ptsInClust = dataSet[nonzero(clusterAssment[:,0].A == cent)[0]]            # nonzero 返回的是矩阵中所有非零元素的坐标，坐标的行数与列数个存在一个数组或矩阵当中            # 矩阵支持检查元素的操作，所有可以写成matrix == int这种形式，返回的一个布尔型矩阵，代表矩阵相应位置有无此元素            # 这里指寻找当前质心下所聚类的样本            centRoids[cent,:] = mean(ptsInClust,axis = 0) # 更新当前的质心为所有样本的平均值，axis = 0代表对列求平均值    return centRoids,clusterAssmentdef plotKmens(dataSet,k,clusterMeans):    '''    本函数用于绘制kMeans的二维聚类图    :param dataSet: 数据集    :param k: 聚类的个数    :return:无    '''    centPoids,assment = clusterMeans(dataSet,k)    fig = plt.figure()    ax = fig.add_subplot(111)    ax.scatter(dataSet[:,0],dataSet[:,1],c = 'blue')    ax.scatter(centRoids[:,0],centRoids[:,1],c = 'red',marker = '+',s = 70)    plt.show()def binKMeans(dataSet, k, distMeas = getDistance):    '''    本函数用于二分k均值算法    :param dataSet: 数据集，要求有矩阵形式    :param k: 指定聚类个数    :param distMeas: 求解距离的方式    :return:质心，簇索引和误差距离矩阵    '''    m = shape(dataSet)[0]    clusterAssment = mat(zeros((m,2)))    centRoids0 = mean(dataSet,axis = 0).tolist()[0] # 初始化一个簇，只有一个质心，分量就是就是所有特征的均值    # 注意，tolist函数用于将矩阵转化为一个列表，此列表为嵌套列表    #print centRoids0    centList = [centRoids0]    for j in range(m): # 遍历所有样本，计算所有样本与当前质心的距离作为误差        clusterAssment[j,1] = distMeas(mat(centRoids0),dataSet[j,:])**2    while (len(centList) < k): # 循环条件为当前质心数目还不够指定数目        lowestSSE = inf        for i in range(len(centList)): # 遍历所有质心            ptsCurrCluster = dataSet[nonzero(clusterAssment[:,0].A == i)[0],:] # 搜索到当前质心所聚类的样本            centroidsMat,splitClusterAss = kMeans(ptsCurrCluster,2,distMeas) # 将当前分割成两个簇            sseSplit = sum(splitClusterAss[:,1]) # 计算分裂簇后的SSE            sseNotSplit = sum(clusterAssment[nonzero(clusterAssment[:,0].A != i)[0],1])            # 计算分裂之前的SSE            if (sseSplit + sseNotSplit) < lowestSSE: # 如果分裂之后的SSE小，则更新                bestCent2Split = i                bestNewCents = centroidsMat                bestClustAss = splitClusterAss.copy()                lowestSSE = sseSplit+sseNotSplit        #重新编制簇的编号，凡是分裂后编号为1的簇，编号为质心列表长度，编号为0的簇，编号为最佳分裂质心的编号，以此更新        bestClustAss[nonzero(bestClustAss[:,0].A == 1)[0],0] = len(centList)        bestClustAss[nonzero(bestClustAss[:,0].A == 0)[0],0] = bestCent2Split        centList[bestCent2Split] = bestNewCents[0,:].tolist()[0] # 添加分裂的质心到质心列表中        centList.append(bestNewCents[1,:].tolist()[0])        clusterAssment[nonzero(clusterAssment[:,0].A == bestCent2Split)[0],:] = bestClustAss    return mat(centList),clusterAssmentdef biKmeans(dataSet, k, distMeas=getDistance):    m = shape(dataSet)[0]    clusterAssment = mat(zeros((m,2)))    centroid0 = mean(dataSet, axis=0).tolist()[0]    centList =[centroid0] #create a list with one centroid    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],:]#get the data points currently in cluster i            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) #change 1 to 3,4, or whatever        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]#replace a centroid with two best centroids        centList.append(bestNewCents[1,:].tolist()[0])        clusterAssment[nonzero(clusterAssment[:,0].A == bestCentToSplit)[0],:]= bestClustAss#reassign new clusters, and SSE    return mat(centList), clusterAssment

密度聚类，基本思路就是将所有密度可达的点都归为一簇。

#encoding=utf-8import numpy as npimport kmeans as kmimport matplotlib.pyplot as pltdef createDisMat(dataMat):    m = dataMat.shape[0]    n = dataMat.shape[1]    distMat = np.mat(np.zeros((m,m))) # 初始化距离矩阵，这里默认使用欧式距离    for i in range(m):        for j in range(m):            if i == j:                distMat[i,j] = 0            else:                dist = km.getDistance(dataMat[i,:],dataMat[j,:])                distMat[i,j] = dist                distMat[j,i] = dist    return distMatdef findCore(dataMat,delta,minPts):    core = []    m = dataMat.shape[0]    n = dataMat.shape[1]    distMat = createDisMat(dataMat)    for i in range(m):        temp = distMat[i,:] < delta # 单独抽取矩阵一行做过滤，凡是小于邻域值的都被标记位True类型        ptsNum = np.sum(temp,1) # 按行加和，统计小于邻域值的点个数        if ptsNum >= minPts:            core.append(i) # 满足条件，增加核心点    return coredef DBSCAN(dataMat,delta,minPts):    k = 0    m = dataMat.shape[0]    distMat = createDisMat(dataMat) # 获取距离矩阵    core = findCore(dataMat,delta,minPts) # 获取核心点列表    unVisit = [1] * m # hash值作为标记，当某一位置的数据位1时，表示还未被访问，为0表示已经被访问    Q = []    ck = []    unVistitOld = []    while len(core) != 0:        print 'a'        unVistitOld = unVisit[:] # 保留原始的未被访问集        i = np.random.choice(core) # 在核心点集中随机选择样本        Q.append(i) # 加入对列Q        unVisit[i] = 0 #剔除当前加入对列的数据，表示已经访问到了        while len(Q) != 0:            print len(Q)            temp = distMat[Q[0],:]<delta # 获取在此核心点邻域范围内的点集            del Q[0]            ptsNum = np.sum(temp,1)            if ptsNum >= minPts:                for j in range(len(unVisit)):                    if unVisit[j] == 1 and temp[0,j] == True:                        Q.append(j)                        unVisit[j] = 0        k += 1        ck.append([])        for index in range(m):            if unVistitOld[index] == 1 and unVisit[index] == 0: # 上一轮未被访问到此轮被访问到的点均要加入当前簇                ck[k-1].append(index)                if index in core: # 在核心点集中清除当前簇的点                    del core[core.index(index)]    return ckdef plotAns(dataSet,ck):    fig = plt.figure()    ax = fig.add_subplot(111)    ax.scatter(dataSet[ck[0],0],dataSet[ck[0],1],c = 'blue')    ax.scatter(dataSet[ck[1],0],dataSet[ck[1],1],c = 'red')    ax.scatter(dataSet[ck[2],0],dataSet[ck[2],1],c = 'green')    ax.scatter(dataSet[ck[3],0],dataSet[ck[3],1],c = 'yellow')    #ax.scatter(centRoids[:,0],centRoids[:,1],c = 'red',marker = '+',s = 70)    plt.show()if __name__ == '__main__':    dataMat = km.loadDataSet("testSet.txt")    # distMat = createDisMat(dataMat)    # core = findCore(dataMat,1,5)    # print distMat    # print len(core)    ck = DBSCAN(dataMat,2,15)    print ck    print len(ck)    plotAns(dataMat,ck)

层次聚类，核心是定义了簇之间的距离衡量，不断寻找距离最近的簇归为一簇。

#encoding=utf-8import numpy as npimport DBSCAN as dbimport kmeans as kmdef calcDistByMin(dataMat,ck1,ck2): # 最小距离点作为簇间的距离    min = np.inf    for vec1 in ck1:        for vec2 in ck2:            dist = km.getDistance(dataMat[vec1,:],dataMat[vec2,:])            if dist <= min:                min = dist    return mindef calcDistByMax(dataMat,ck1,ck2): # 最大距离点作为簇间的距离    max = 0    for vec1 in ck1:        for vec2 in ck2:            dist = km.getDistance(dataMat[vec1,:],dataMat[vec2,:])            if dist >= max:                max = dist    return maxdef createDistMat(dataMat,calcDistType = calcDistByMin): # 生成初始的距离矩阵    m = dataMat.shape[0]    distMat = np.mat(np.zeros((m,m)))    for i in range(m):        for j in range(m):            listI = [i];listJ = [j] # 为配合距离函数的输入参数形式，在这里要列表化一下            distMat[i,j] = calcDistType(dataMat,listI,listJ)            distMat[j,i] = distMat[i,j]    return distMatdef findMaxLoc(distMat,q): # 寻找矩阵中最小的元素并返回其位置，注意，这里不能返回相同的坐标    min = np.inf    I = J = 0    for i in range(q):        for j in range(q):            if distMat[i,j] < min and i != j:                min = distMat[i,j]                I = i                J = j    return I,Jdef ANGES(dataMat,k,calcDistType = calcDistByMax):    m = dataMat.shape[0]    ck = []    for i in range(m):        ck.append([i])    distMat = createDistMat(dataMat,calcDistType)    q = m # 初始化点集个数    while q > k:        i,j = findMaxLoc(distMat,q)        #print i,j        if i > j:            i,j = j,i # 保证i<j，这样做是为了删除的是序号较大的簇        ck[i].extend(ck[j]) # 把序号较大的簇并入序号小的簇        del ck[j] # 删除序号大的簇        distMat = np.delete(distMat,j,0) # 在距离矩阵中删除该簇的数据，注意这里delete函数有返回值，否则不会有删除作用        distMat = np.delete(distMat,j,1)        print distMat.shape        for index in range(0,q-1): # 重新计算新簇和其余簇之间的距离            distMat[i,index] = calcDistType(dataMat,ck[i],ck[index])            distMat[i,index] = distMat[index,i]        q -= 1 # 一个点被分入簇中，自减    return ckif __name__ == '__main__':    dataMat = km.loadDataSet("testSet.txt")    ck = ANGES(dataMat,4)    print ck    db.plotAns(dataMat,ck)