程序博客网 > 广州数据库开发工程师

k均值聚类，密度聚类，层次聚类

来源：互联网 发布：广州数据库开发工程师 编辑：程序博客网 时间：2024/04/28 06:54

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

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

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

[python] view plain copy
 在CODE上查看代码片派生到我的代码片
  1. #encoding=utf-8  
  2. __author__ = 'freedom'  
  3.   
  4. from numpy import*  
  5. import matplotlib.pyplot as plt  
  6.   
  7. def loadDataSet(fileName):  
  8.     ''''' 
  9.     本函数用于加载数据 
  10.     :param fileName: 数据文件名 
  11.     :return:数据集，具有矩阵形式 
  12.     '''  
  13.     fr = open(fileName)  
  14.     dataSet = []  
  15.     for line in fr.readlines():  
  16.         curLine = line.strip().split('\t')  
  17.         inLine = map(float,curLine) # 利用map广播，是的读入的字符串变为浮点型  
  18.         dataSet.append(inLine)  
  19.     return mat(dataSet)  
  20.   
  21. def getDistance(vecA,vecB):  
  22.     ''''' 
  23.     本函数用于计算欧氏距离 
  24.     :param vecA: 向量A 
  25.     :param vecB: 向量B 
  26.     :return:欧氏距离 
  27.     '''  
  28.     return sqrt(sum(power(vecA-vecB,2)))  
  29.   
  30. def randCent(dataSet,k):  
  31.     ''''' 
  32.     本函数用于生成k个随机质心 
  33.     :param dataSet: 数据集，具有矩阵形式 
  34.     :param k:指定的质心个数 
  35.     :return:随机质心，具有矩阵形式 
  36.     '''  
  37.     n = shape(dataSet)[1] # 获取特征数目  
  38.     centRoids = mat(zeros((k,n)))  
  39.     for j in range(n):  
  40.         minJ = min(dataSet[:,j]) # 获取每个特征的最小值  
  41.         rangeJ = float(max(dataSet[:,j]-minJ)) # 获取每个特征的范围  
  42.         centRoids[:,j] = minJ + rangeJ*random.rand(k,1) # numpy下的rand表示随机生成k*1的随机数矩阵，范围0-1  
  43.     return centRoids  
  44.   
  45. def kMeans(dataSet,k,disMens = getDistance,createCent = randCent):  
  46.     ''''' 
  47.     本函数用于k均值聚类 
  48.     :param dataSet: 数据集，要求有矩阵形式 
  49.     :param k: 指定聚类的个数 
  50.     :param disMens: 求解距离的方式，除欧式距离还可以定义其他距离计算方式 
  51.     :param createCent: 生成随机质心方式 
  52.     :return:随机质心，簇索引和误差距离矩阵 
  53.     '''  
  54.     m = shape(dataSet)[0]  
  55.     clusterAssment = mat(zeros((m,2))) # 要为每个样本建立一个簇索引和相对的误差，所以需要m行的矩阵，m就是样本数  
  56.     centRoids = createCent(dataSet,k) # 生成随机质心  
  57.     clusterChanged = True  
  58.     while clusterChanged:  
  59.         clusterChanged = False  
  60.         for i in range(m): # 遍历所有样本  
  61.             minDist = inf;minIndex = -1 # 初始化最小值  
  62.             for j in range(k): # 遍历所有质心  
  63.                 disJI = disMens(centRoids[j,:],dataSet[i,:])  
  64.                 if disJI < minDist:  
  65.                     minDist = disJI;minIndex = j # 找出距离当前样本最近的那个质心  
  66.             if clusterAssment[i,0] != minIndex: # 更新当前样本点所属于的质心  
  67.                 clusterChanged = True # 如果当前样本点不属于当前与之距离最小的质心，则说明簇分配结果仍需要改变  
  68.                 clusterAssment[i,:] = minIndex,minDist**2  
  69.         for cent in range(k):  
  70.             ptsInClust = dataSet[nonzero(clusterAssment[:,0].A == cent)[0]]  
  71.             # nonzero 返回的是矩阵中所有非零元素的坐标，坐标的行数与列数个存在一个数组或矩阵当中  
  72.             # 矩阵支持检查元素的操作，所有可以写成matrix == int这种形式，返回的一个布尔型矩阵，代表矩阵相应位置有无此元素  
  73.             # 这里指寻找当前质心下所聚类的样本  
  74.             centRoids[cent,:] = mean(ptsInClust,axis = 0) # 更新当前的质心为所有样本的平均值，axis = 0代表对列求平均值  
  75.     return centRoids,clusterAssment  
  76.   
  77. def plotKmens(dataSet,k,clusterMeans):  
  78.     ''''' 
  79.     本函数用于绘制kMeans的二维聚类图 
  80.     :param dataSet: 数据集 
  81.     :param k: 聚类的个数 
  82.     :return:无 
  83.     '''  
  84.     centPoids,assment = clusterMeans(dataSet,k)  
  85.     fig = plt.figure()  
  86.     ax = fig.add_subplot(111)  
  87.     ax.scatter(dataSet[:,0],dataSet[:,1],c = 'blue')  
  88.     ax.scatter(centRoids[:,0],centRoids[:,1],c = 'red',marker = '+',s = 70)  
  89.     plt.show()  
  90.   
  91. def binKMeans(dataSet, k, distMeas = getDistance):  
  92.     ''''' 
  93.     本函数用于二分k均值算法 
  94.     :param dataSet: 数据集，要求有矩阵形式 
  95.     :param k: 指定聚类个数 
  96.     :param distMeas: 求解距离的方式 
  97.     :return:质心，簇索引和误差距离矩阵 
  98.     '''  
  99.     m = shape(dataSet)[0]  
  100.     clusterAssment = mat(zeros((m,2)))  
  101.     centRoids0 = mean(dataSet,axis = 0).tolist()[0] # 初始化一个簇，只有一个质心，分量就是就是所有特征的均值  
  102.     # 注意，tolist函数用于将矩阵转化为一个列表，此列表为嵌套列表  
  103.     #print centRoids0  
  104.     centList = [centRoids0]  
  105.     for j in range(m): # 遍历所有样本，计算所有样本与当前质心的距离作为误差  
  106.         clusterAssment[j,1] = distMeas(mat(centRoids0),dataSet[j,:])**2  
  107.     while (len(centList) < k): # 循环条件为当前质心数目还不够指定数目  
  108.         lowestSSE = inf  
  109.         for i in range(len(centList)): # 遍历所有质心  
  110.             ptsCurrCluster = dataSet[nonzero(clusterAssment[:,0].A == i)[0],:] # 搜索到当前质心所聚类的样本  
  111.             centroidsMat,splitClusterAss = kMeans(ptsCurrCluster,2,distMeas) # 将当前分割成两个簇  
  112.             sseSplit = sum(splitClusterAss[:,1]) # 计算分裂簇后的SSE  
  113.             sseNotSplit = sum(clusterAssment[nonzero(clusterAssment[:,0].A != i)[0],1])  
  114.             # 计算分裂之前的SSE  
  115.             if (sseSplit + sseNotSplit) < lowestSSE: # 如果分裂之后的SSE小，则更新  
  116.                 bestCent2Split = i  
  117.                 bestNewCents = centroidsMat  
  118.                 bestClustAss = splitClusterAss.copy()  
  119.                 lowestSSE = sseSplit+sseNotSplit  
  120.         #重新编制簇的编号，凡是分裂后编号为1的簇，编号为质心列表长度，编号为0的簇，编号为最佳分裂质心的编号，以此更新  
  121.         bestClustAss[nonzero(bestClustAss[:,0].A == 1)[0],0] = len(centList)  
  122.         bestClustAss[nonzero(bestClustAss[:,0].A == 0)[0],0] = bestCent2Split  
  123.         centList[bestCent2Split] = bestNewCents[0,:].tolist()[0] # 添加分裂的质心到质心列表中  
  124.         centList.append(bestNewCents[1,:].tolist()[0])  
  125.         clusterAssment[nonzero(clusterAssment[:,0].A == bestCent2Split)[0],:] = bestClustAss  
  126.     return mat(centList),clusterAssment  
  127.   
  128. def biKmeans(dataSet, k, distMeas=getDistance):  
  129.     m = shape(dataSet)[0]  
  130.     clusterAssment = mat(zeros((m,2)))  
  131.     centroid0 = mean(dataSet, axis=0).tolist()[0]  
  132.     centList =[centroid0] #create a list with one centroid  
  133.     for j in range(m):#calc initial Error  
  134.         clusterAssment[j,1] = distMeas(mat(centroid0), dataSet[j,:])**2  
  135.     while (len(centList) < k):  
  136.         lowestSSE = inf  
  137.         for i in range(len(centList)):  
  138.             ptsInCurrCluster = dataSet[nonzero(clusterAssment[:,0].A==i)[0],:]#get the data points currently in cluster i  
  139.             centroidMat, splitClustAss = kMeans(ptsInCurrCluster, 2, distMeas)  
  140.             sseSplit = sum(splitClustAss[:,1])#compare the SSE to the currrent minimum  
  141.             sseNotSplit = sum(clusterAssment[nonzero(clusterAssment[:,0].A!=i)[0],1])  
  142.             print "sseSplit, and notSplit: ",sseSplit,sseNotSplit  
  143.             if (sseSplit + sseNotSplit) < lowestSSE:  
  144.                 bestCentToSplit = i  
  145.                 bestNewCents = centroidMat  
  146.                 bestClustAss = splitClustAss.copy()  
  147.                 lowestSSE = sseSplit + sseNotSplit  
  148.         bestClustAss[nonzero(bestClustAss[:,0].A == 1)[0],0] = len(centList) #change 1 to 3,4, or whatever  
  149.         bestClustAss[nonzero(bestClustAss[:,0].A == 0)[0],0] = bestCentToSplit  
  150.         print 'the bestCentToSplit is: ',bestCentToSplit  
  151.         print 'the len of bestClustAss is: ', len(bestClustAss)  
  152.         centList[bestCentToSplit] = bestNewCents[0,:].tolist()[0]#replace a centroid with two best centroids  
  153.         centList.append(bestNewCents[1,:].tolist()[0])  
  154.         clusterAssment[nonzero(clusterAssment[:,0].A == bestCentToSplit)[0],:]= bestClustAss#reassign new clusters, and SSE  
  155.     return mat(centList), clusterAssment  
密度聚类，基本思路就是将所有密度可达的点都归为一簇。

[python] view plain copy
 在CODE上查看代码片派生到我的代码片
  1. #encoding=utf-8  
  2. import numpy as np  
  3. import kmeans as km  
  4. import matplotlib.pyplot as plt  
  5.   
  6. def createDisMat(dataMat):  
  7.     m = dataMat.shape[0]  
  8.     n = dataMat.shape[1]  
  9.     distMat = np.mat(np.zeros((m,m))) # 初始化距离矩阵，这里默认使用欧式距离  
  10.     for i in range(m):  
  11.         for j in range(m):  
  12.             if i == j:  
  13.                 distMat[i,j] = 0  
  14.             else:  
  15.                 dist = km.getDistance(dataMat[i,:],dataMat[j,:])  
  16.                 distMat[i,j] = dist  
  17.                 distMat[j,i] = dist  
  18.     return distMat  
  19.   
  20. def findCore(dataMat,delta,minPts):  
  21.     core = []  
  22.     m = dataMat.shape[0]  
  23.     n = dataMat.shape[1]  
  24.     distMat = createDisMat(dataMat)  
  25.     for i in range(m):  
  26.         temp = distMat[i,:] < delta # 单独抽取矩阵一行做过滤，凡是小于邻域值的都被标记位True类型  
  27.         ptsNum = np.sum(temp,1) # 按行加和，统计小于邻域值的点个数  
  28.         if ptsNum >= minPts:  
  29.             core.append(i) # 满足条件，增加核心点  
  30.     return core  
  31.   
  32. def DBSCAN(dataMat,delta,minPts):  
  33.     k = 0  
  34.     m = dataMat.shape[0]  
  35.     distMat = createDisMat(dataMat) # 获取距离矩阵  
  36.     core = findCore(dataMat,delta,minPts) # 获取核心点列表  
  37.     unVisit = [1] * m # hash值作为标记，当某一位置的数据位1时，表示还未被访问，为0表示已经被访问  
  38.     Q = []  
  39.     ck = []  
  40.     unVistitOld = []  
  41.     while len(core) != 0:  
  42.         print 'a'  
  43.         unVistitOld = unVisit[:] # 保留原始的未被访问集  
  44.         i = np.random.choice(core) # 在核心点集中随机选择样本  
  45.         Q.append(i) # 加入对列Q  
  46.         unVisit[i] = 0 #剔除当前加入对列的数据，表示已经访问到了  
  47.         while len(Q) != 0:  
  48.             print len(Q)  
  49.             temp = distMat[Q[0],:]<delta # 获取在此核心点邻域范围内的点集  
  50.             del Q[0]  
  51.             ptsNum = np.sum(temp,1)  
  52.             if ptsNum >= minPts:  
  53.                 for j in range(len(unVisit)):  
  54.                     if unVisit[j] == 1 and temp[0,j] == True:  
  55.                         Q.append(j)  
  56.                         unVisit[j] = 0  
  57.         k += 1  
  58.         ck.append([])  
  59.         for index in range(m):  
  60.             if unVistitOld[index] == 1 and unVisit[index] == 0: # 上一轮未被访问到此轮被访问到的点均要加入当前簇  
  61.                 ck[k-1].append(index)  
  62.                 if index in core: # 在核心点集中清除当前簇的点  
  63.                     del core[core.index(index)]  
  64.     return ck  
  65.   
  66. def plotAns(dataSet,ck):  
  67.     fig = plt.figure()  
  68.     ax = fig.add_subplot(111)  
  69.     ax.scatter(dataSet[ck[0],0],dataSet[ck[0],1],c = 'blue')  
  70.     ax.scatter(dataSet[ck[1],0],dataSet[ck[1],1],c = 'red')  
  71.     ax.scatter(dataSet[ck[2],0],dataSet[ck[2],1],c = 'green')  
  72.     ax.scatter(dataSet[ck[3],0],dataSet[ck[3],1],c = 'yellow')  
  73.   
  74.     #ax.scatter(centRoids[:,0],centRoids[:,1],c = 'red',marker = '+',s = 70)  
  75.     plt.show()  
  76.   
  77. if __name__ == '__main__':  
  78.     dataMat = km.loadDataSet("testSet.txt")  
  79.     # distMat = createDisMat(dataMat)  
  80.     # core = findCore(dataMat,1,5)  
  81.     # print distMat  
  82.     # print len(core)  
  83.     ck = DBSCAN(dataMat,2,15)  
  84.     print ck  
  85.     print len(ck)  
  86.     plotAns(dataMat,ck)  

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

[python] view plain copy
 在CODE上查看代码片派生到我的代码片
  1. #encoding=utf-8  
  2. import numpy as np  
  3. import DBSCAN as db  
  4. import kmeans as km  
  5.   
  6.   
  7. def calcDistByMin(dataMat,ck1,ck2): # 最小距离点作为簇间的距离  
  8.     min = np.inf  
  9.     for vec1 in ck1:  
  10.         for vec2 in ck2:  
  11.             dist = km.getDistance(dataMat[vec1,:],dataMat[vec2,:])  
  12.             if dist <= min:  
  13.                 min = dist  
  14.     return min  
  15.   
  16. def calcDistByMax(dataMat,ck1,ck2): # 最大距离点作为簇间的距离  
  17.     max = 0  
  18.     for vec1 in ck1:  
  19.         for vec2 in ck2:  
  20.             dist = km.getDistance(dataMat[vec1,:],dataMat[vec2,:])  
  21.             if dist >= max:  
  22.                 max = dist  
  23.     return max  
  24.   
  25. def createDistMat(dataMat,calcDistType = calcDistByMin): # 生成初始的距离矩阵  
  26.     m = dataMat.shape[0]  
  27.     distMat = np.mat(np.zeros((m,m)))  
  28.     for i in range(m):  
  29.         for j in range(m):  
  30.             listI = [i];listJ = [j] # 为配合距离函数的输入参数形式，在这里要列表化一下  
  31.             distMat[i,j] = calcDistType(dataMat,listI,listJ)  
  32.             distMat[j,i] = distMat[i,j]  
  33.     return distMat  
  34.   
  35. def findMaxLoc(distMat,q): # 寻找矩阵中最小的元素并返回其位置，注意，这里不能返回相同的坐标  
  36.     min = np.inf  
  37.     I = J = 0  
  38.     for i in range(q):  
  39.         for j in range(q):  
  40.             if distMat[i,j] < min and i != j:  
  41.                 min = distMat[i,j]  
  42.                 I = i  
  43.                 J = j  
  44.     return I,J  
  45.   
  46.   
  47. def ANGES(dataMat,k,calcDistType = calcDistByMax):  
  48.     m = dataMat.shape[0]  
  49.     ck = []  
  50.     for i in range(m):  
  51.         ck.append([i])  
  52.     distMat = createDistMat(dataMat,calcDistType)  
  53.     q = m # 初始化点集个数  
  54.     while q > k:  
  55.         i,j = findMaxLoc(distMat,q)  
  56.         #print i,j  
  57.         if i > j:  
  58.             i,j = j,i # 保证i<j，这样做是为了删除的是序号较大的簇  
  59.         ck[i].extend(ck[j]) # 把序号较大的簇并入序号小的簇  
  60.         del ck[j] # 删除序号大的簇  
  61.         distMat = np.delete(distMat,j,0) # 在距离矩阵中删除该簇的数据，注意这里delete函数有返回值，否则不会有删除作用  
  62.         distMat = np.delete(distMat,j,1)  
  63.         print distMat.shape  
  64.         for index in range(0,q-1): # 重新计算新簇和其余簇之间的距离  
  65.             distMat[i,index] = calcDistType(dataMat,ck[i],ck[index])  
  66.             distMat[i,index] = distMat[index,i]  
  67.         q -= 1 # 一个点被分入簇中，自减  
  68.     return ck  
  69.   
  70. if __name__ == '__main__':  
  71.     dataMat = km.loadDataSet("testSet.txt")  
  72.     ck = ANGES(dataMat,4)  
  73.     print ck  
  74.     db.plotAns(dataMat,ck)
0 0
广州数据库开发工程师 广州数据库开发工程师
原创粉丝点击
热门IT博客
linux acl权限控制linux acl权限控制
联通电信网络转换联通电信网络转换
深圳岂凡网络 怎么样
发票网络认证系统
js 放入localstorage[]
软件开发培训收费
spring boot sql脚本
nba2k17pc版优化
淘宝宝贝排名在线查询
简明python教程电子书
mac有没有抢网速的软件
济宁市网络问政平台app
深圳网络诈骗报警
情感寄语一叶落而知秋
windows samba 客户端
数码宝贝网络侦探存档
农村淘宝村小二工资
java的switch里的
网站域名收费
c语言求质因数
热门问题 老师的惩罚 人脸识别 我在镇武司摸鱼那些年 重生之率土为王 我在大康的咸鱼生活 盘龙之生命进化 天生仙种 凡人之先天五行 春回大明朝 姑娘不必设防，我是瞎子 学化妆师 韩国免税店化妆品价格 打扮 那个彩妆好 彩妆学校哪里好 怎么样护肤 后护肤 好用的护肤产品 妆后 彩妆都有哪些 韩妆 化彩妆的正确步骤 彩妆有那些 如何护肤 专业化妆培训学校 日化产品 彩妆有哪些 化妆师好吗 什么护肤好 日妆 什么彩妆好用 妆容 每日护肤 好用的彩妆 护肤过程 初学者如何学化妆 怎么学会给自己化妆 画淡妆的步骤及化妆用品 慈溪哪里可以学化妆 彩妆教程教你如何化妆 化学符号 化学周期表 化学 化学方程式 氨气化学式 化学实验 明矾化学式 酒精化学式 甲烷化学式 中国化学 纯碱化学式

程序博客网，程序员的互联网技术博客家园。csdn论坛精品 msdn技术资料都在这里