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

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

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

算法思想

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

首先定義兩個值:

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

其中：

?? 是一個截斷距離, 是一個超參數(shù).

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

? ? ? ? ??

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

聚類過程

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

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

聚類分析

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

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

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

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

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()

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

?

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。