論文. Clustering by fast search and find of density peak.?Alex Rodriguez, Alessandro Laio

參考鏈接：Science上發表的超贊聚類算法--

作者(Alex Rodriguez, Alessandro Laio)提出了一種很簡潔優美的聚類算法, 可以識別各種形狀的類簇, 并且其超參數很容易確定.

算法思想

該算法的假設是, 類簇的中心由一些局部密度比較低的點圍繞, 并且這些點距離其他高局部密度的點的距離都比較大.

首先定義兩個值:

????? 局部密度,? 以及到高局部密度點的距離:?

其中：

?? 是一個截斷距離, 是一個超參數.

所以相當于距離點i的距離小于的點的個數. 由于該算法只對的相對值敏感, 所以對的選擇比較魯棒, 一種推薦做法是選擇使得平均每個點的鄰居數為所有點的1%-2%.

? ? ? ? ??

對于密度最大的點, 設置. 注意只有那些密度是局部或者全局最大的點才會有遠大于正常的相鄰點間距.

聚類過程

???? 那些有著比較大的局部密度和很大的的點被認為是類簇的中心. 局部密度較小但是較大的點是異常點.在確定了類簇中心之后, 所有其他點屬于距離其最近的類簇中心所代表的類簇. 圖例如下:

???????左圖是所有點在二維空間的分布, 右圖是以為橫坐標, 以為縱坐標, 這種圖稱作決策圖(decision tree). 可以看到, 1和10兩個點的和都比較大, 作為類簇的中心點. 26, 27, 28三個點的也比較大, 但是較小, 所以是異常點.

聚類分析

??????? 在聚類分析中, 通常需要確定每個點劃分給某個類簇的可靠性. 在該算法中, 可以首先為每個類簇定義一個邊界區域(border region), 亦即劃分給該類簇但是距離其他類簇的點的距離小于的點. 然后為每個類簇找到其邊界區域的局部密度最大的點, 令其局部密度為. 該類簇中所有局部密度大于的點被認為是類簇核心的一部分(亦即將該點劃分給該類簇的可靠性很大), 其余的點被認為是該類簇的光暈(halo), 亦即可以認為是噪音. 圖例如下

???????A圖為生成數據的概率分布, B, C二圖為分別從該分布中生成了4000, 1000個點. D, E分別是B, C兩組數據的決策圖(decision tree), 可以看到兩組數據都只有五個點有比較大的和很大的. 這些點作為類簇的中心, 在確定了類簇的中心之后, 每個點被劃分到各個類簇(彩色點), 或者是劃分到類簇光暈(黑色點). F圖展示的是隨著抽樣點數量的增多, 聚類的錯誤率在逐漸下降, 說明該算法是魯棒的.

???????最后展示一下該算法在各種數據分布上的聚類效果, 非常贊.

數據集DATA來源:http://cs.joensuu.fi/sipu/datasets/ ，測試數據集是幾個文本文件，可以直接看著分析。

Python代碼：原始鏈接：http://www.chinakdd.com/article-Mut2iwV7S3i16cj.html

#coding=utf-8 import math import pylab as pldef getDistance(pt1, pt2):tmp = pow(pt1[0]-pt2[0],2) + pow(pt1[1]-pt2[1],2)#dis = pow(tmp,0.5)dis = math.sqrt(math.fabs(tmp) )return disdef ChooseDc(dc_percent,points,dis,distance):avgNeighbourNum = dc_percent*len(points)maxd =0for i in range(0,len(points)):for j in range(i ,len(points)):pt1 = points[i]pt2 = points[j]d = getDistance(pt1,pt2)dis.append(d)distance[i,j]= ddis.append(d)distance[j,i]= dif d>maxd:maxd = d#print disdis.sort()#print disprint avgNeighbourNumprint len(points)*2dc = dis[ int ( avgNeighbourNum* len(points)*2)]#print dcreturn dcdef drawOriginGraph(pl,points,cl,colorNum):x =[xx for(xx,yy)in points]y =[yy for(xx,yy)in points]cm = pl.get_cmap("RdYlGn")#print cl#for i in range(len(points)):# if cl[i]==0:# pl.plot(x[i],y[i],'o', color =cm(cl[i]*1.0/colorNum))# else:# pl.plot(x[i],y[i],'*', color =cm(cl[i]*1.0/colorNum))for i in range(len(points)):pl.plot(x[i],y[i],'o', color =cm( math.fabs(cl[i]) *1.0/colorNum))#pl.plot(x[i],y[i],'o', color =cm( cl[i] *1.0/colorNum))def drawDecisionGraph(pl,rho, delta,cl,colorNum):cm = pl.get_cmap("RdYlGn")for i in range(len(rho)):pl.plot(delta[i],rho[i], 'o',color=cm( math.fabs(cl[i]) *1.0/colorNum ))#pl.xlabel(r'$rho$')#pl.ylabel(r'$delta$')pl.xlabel('wishchin')pl.ylabel('delta')def Cluster(dc_percent): #=========Load Data=========InputFileName = "flame"InputFileName = "Compound"InputFileName = "spiral"FolderName ="E:\Develope\EclipseWorks\MeachinLearning/Ch19_ScineceCluster/"InputFileName = FolderName + InputFileNameOutputFileName= InputFileName +"_out"suffix =".txt"Fin= open(InputFileName+ suffix,"r")Fout= open(OutputFileName+ suffix,"w")points =[]for line in Fin.readlines():data = line.split()if len(data)==3:a =float(data[0])b =float(data[1])points.append((a,b))#=========Calculating=========#-----choose dc-----#dc_percent = 0.015 #0.5 #0.015dis =[]distance ={}dc =ChooseDc( dc_percent, points, dis, distance)print("dc:" ,str(dc))#-----cal rho:"Cut off" kernel#''' # rho = [0 for i in range(len(points))] # for i in range(0,len(points)): # for j in range(i 1,len(points)): # dij = getDistance(points[i],points[j]) # if dij<dc: # rho[i] = 1 # rho[j] = 1 # '''#-----cal rho:"Gaussian Kernel"rho =[i for i in range(len(points))]for i in range(0,len(points)):for j in range(i ,len(points)):dij = getDistance(points[i],points[j])#print (dij, dc)#高斯核！！rho[i] = math.exp(-(dij/dc)*(dij/dc))rho[j] = math.exp(-(dij/dc)*(dij/dc))rho_list =[(rho[i],i)for i in range(len(rho))]rho_sorted = sorted(rho_list, reverse=1)print("Highest rho:",rho_sorted[0])maxd = dis[-1]delta =[maxd for i in range(len(points))]nneigh =[-1for i in range(len(points))]for ii in range(1,len(rho_sorted)):for jj in range(0,ii):id_p1 = rho_sorted[ii][1]#get point1's idid_p2 = rho_sorted[jj][1]#get point2's idif(distance[id_p1,id_p2]<delta[id_p1]):delta[id_p1]= distance[id_p1,id_p2]nneigh[id_p1]= id_p2#assignmentcl =[-1 for i in range(len(points))]colorNum =0for ii in range(len(rho_sorted)):id_p = rho_sorted[ii][1]if(cl[id_p]==-1 and delta[id_p]>7.0):cl[id_p]= colorNumcolorNum =1else:if(cl[id_p]==-1 and cl[nneigh[id_p]!=-1]):cl[id_p]= cl[nneigh[id_p]]#print(colorNum)#import pylab as plfig1 = pl.figure(1)pl.subplot(121)#pl.subplot(211)drawOriginGraph(pl,points,cl,colorNum)pl.subplot(122)drawDecisionGraph(pl,rho,delta,cl,colorNum)pl.show()for i in range(len(points)):Fout.write(str(i) +"," +str(rho[i])+ "," +str(delta[i])+ "\n")#Assign ClusterFin.close()Fout.close()#if __name__=="__main__":#Cluster()

測試：

import Cluster Cluster.Cluster()

后記：發表在 Science 上的一種新聚類算法-

?

總結

以上是生活随笔為你收集整理的Science上发表的超赞聚类算法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。