头歌第一关密度聚类算法

文章目录

• 前言
• 一、DPC是什么？
• 二、算法参数介绍

• 1.局部密度
• 2.距离

• 3.代码部分

• 3.1 导入数据包
• 3.2 计算数据点两两之间的距离
• 3.3找到密度计算的阈值dc
• 3.4 计算每个点的局部密度
• 3.5 找到距离其最近的点的距离
• 3.6 选择聚类中心
• 3.7 对除聚类中心的点进行聚类
• 3.8 代码汇总

• 4.实验

• 4.1人造数据集
• 4.2鸢尾花数据集

• 总结
• 参考文献

---

前言

2014 年 6 月,一篇题为《Clustering by fast search and find of density peaks》的文章刊登在《Science》杂志,该论文介绍了一种新的聚类算法——密度峰值聚类算法,该算法是一种简洁、高效的聚类算法,已成功应用于多个领域。密度峰值聚类算法是综合了样本对象的密度与距离所提出的一种新的聚类算法，该算法的核心步骤主要分为两步完成：第一步依据决策图选择类簇的中心点；第二步采用一步分配策略分配剩余样本点。当确定聚类中心点之后，其余的数据对象点将会分配给大于自身密度值并与之相距最近的类簇，将计算所得的局部密度值按照从大到小的顺序依次更新所有的数据点直至聚类完成，该算法不需要进行复杂的迭代计算，适用性较强，对于任意形状的类簇都可以准确识别。

一、DPC是什么？

DPC 算法具有一种简单且高效的设计思想，即类簇质心点的密度均高于周围邻居的局部密度数值，而且中心点与密度较高的数据点之间的距离都比较远。该聚类算法思想在实现的过程中主要使用了以下两个概念：样本点的局部密度、样本点的距离(这两个概念定义恰好对应算法设定的两个假设)。

二、算法参数介绍

1.局部密度

$\rho_{i}=\sum_{j}\left(-\exp \left(\frac{d_{i j}}{d_{c}}\right)\right)$

2.距离

$\left\{\begin{array}{l} \delta_{i}=\min _{j: \rho_{j}>\rho i}\left(d_{i j}\right) \\ \delta_{i}=\max _{j \in c}\left(d_{i j}\right) \end{array}\right.$
如果数据点${x}_{i}$的局部密度取值最大时，那么该样本点的距离长度为在所有数据中与数据点${x}_{i}$之间距离最远的数据点到数据点${x}_{i}$的距离；否则，该样本点的距离表示的是该样本点与其距离最近，且密度高于自身的样本距离。也就是说，局部密度越大的样本其距离越大。

3.代码部分

3.1 导入数据包

import numpy as np
import matplotlib.pyplot as plt
import csv

3.2 计算数据点两两之间的距离

def getDistanceMatrix(datas):
    N,D = np.shape(datas)
    dists = np.zeros([N,N])
    
    for i in range(N):
        for j in range(N):
            vi = datas[i,:]
            vj = datas[j,:]
            dists[i,j]= np.sqrt(np.dot((vi-vj),(vi-vj)))
    return dists

这里的距离指的是欧式距离，得到的结果为一个N*N的矩阵。

3.3找到密度计算的阈值dc

这里dc有着多种选择方法，最简单的一种是手动进行选择，还有一种要求平均每个点周围距离小于dc的点的数目占总点数的1%-2%。

def select_dc(dists):    
    N = np.shape(dists)[0]
    tt = np.reshape(dists,N*N)
    percent = 2.0
    position = int(N * (N - 1) * percent / 100)
    dc = np.sort(tt)[position  + N]

3.4 计算每个点的局部密度

def get_density(dists,dc,method=None):
    N = np.shape(dists)[0]
    rho = np.zeros(N)
    
    for i in range(N):
        if method == None:
            rho[i]  = np.where(dists[i,:]<dc)[0].shape[0]-1
        else:
            rho[i] = np.sum(np.exp(-(dists[i,:]/dc)**2))-1
    return rho

3.5 找到距离其最近的点的距离

def get_deltas(dists,rho):
    N = np.shape(dists)[0]
    deltas = np.zeros(N)
    nearest_neiber = np.zeros(N)
    # 将密度从大到小排序
    index_rho = np.argsort(-rho)
    for i,index in enumerate(index_rho):
        # 对于密度最大的点
        if i==0:
            continue
            
        # 对于其他的点
        # 找到密度比其大的点的序号    
        index_higher_rho = index_rho[:i]
        # 获取这些点距离当前点的距离,并找最小值
        deltas[index] = np.min(dists[index,index_higher_rho])
        
        #保存最近邻点的编号
        index_nn = np.argmin(dists[index,index_higher_rho])
        nearest_neiber[index] = index_higher_rho[index_nn].astype(int)
    
    deltas[index_rho[0]] = np.max(deltas)   
    return deltas,nearest_neiber

3.6 选择聚类中心

选取 rho与delta乘积较大的点作为聚类中心

def find_centers_K(rho,deltas,K):
    rho_delta = rho*deltas
    centers = np.argsort(-rho_delta)
    return centers[:K]

3.7 对除聚类中心的点进行聚类

def cluster_PD(rho,centers,nearest_neiber):
    K = np.shape(centers)[0]
    if K == 0:
        print("can not find centers")
        return
    
    N = np.shape(rho)[0]
    labs = -1*np.ones(N).astype(int)
    
    # 首先对几个聚类中进行标号
    for i, center in enumerate(centers):
        labs[center] = i
   
    # 将密度从大到小排序
    index_rho = np.argsort(-rho)
    for i, index in enumerate(index_rho):
        # 从密度大的点进行标号
        if labs[index] == -1:
            # 如果没有被标记过
            # 那么聚类标号与距离其最近且密度比其大
            # 的点的标号相同
            labs[index] = labs[int(nearest_neiber[index])]
    return labs

3.8 代码汇总

import numpy as np
import matplotlib.pyplot as plt
import csv

# 计算数据点两两之间的距离
def getDistanceMatrix(datas):
    N,D = np.shape(datas)
    dists = np.zeros([N,N])
    
    for i in range(N):
        for j in range(N):
            vi = datas[i,:]
            vj = datas[j,:]
            dists[i,j]= np.sqrt(np.dot((vi-vj),(vi-vj)))
    return dists
    
# 找到密度计算的阈值dc
# 要求平均每个点周围距离小于dc的点的数目占总点数的1%-2%
def select_dc(dists):    
    '''算法1'''
    N = np.shape(dists)[0]
    tt = np.reshape(dists,N*N)
    percent = 2.0
    position = int(N * (N - 1) * percent / 100)
    dc = np.sort(tt)[position  + N]
    
    ''' 算法 2 '''
    # N = np.shape(dists)[0]
    # max_dis = np.max(dists)
    # min_dis = np.min(dists)
    # dc = (max_dis + min_dis) / 2
   
    # while True:
        # n_neighs = np.where(dists<dc)[0].shape[0]-N
        # rate = n_neighs/(N*(N-1))
        
        # if rate>=0.01 and rate<=0.02:
            # break
        # if rate<0.01:
            # min_dis = dc
        # else:
            # max_dis = dc
            
        # dc = (max_dis + min_dis) / 2
        # if max_dis - min_dis < 0.0001:
            # break
    return dc
    
# 计算每个点的局部密度    
def get_density(dists,dc,method=None):
    N = np.shape(dists)[0]
    rho = np.zeros(N)
    
    for i in range(N):
        if method == None:
            rho[i]  = np.where(dists[i,:]<dc)[0].shape[0]-1
        else:
            rho[i] = np.sum(np.exp(-(dists[i,:]/dc)**2))-1
    return rho
    
# 计算每个数据点的密度距离
# 即对每个点，找到密度比它大的所有点
# 再在这些点中找到距离其最近的点的距离
def get_deltas(dists,rho):
    N = np.shape(dists)[0]
    deltas = np.zeros(N)
    nearest_neiber = np.zeros(N)
    # 将密度从大到小排序
    index_rho = np.argsort(-rho)
    for i,index in enumerate(index_rho):
        # 对于密度最大的点
        if i==0:
            continue
            
        # 对于其他的点
        # 找到密度比其大的点的序号    
        index_higher_rho = index_rho[:i]
        # 获取这些点距离当前点的距离,并找最小值
        deltas[index] = np.min(dists[index,index_higher_rho])
        
        #保存最近邻点的编号
        index_nn = np.argmin(dists[index,index_higher_rho])
        nearest_neiber[index] = index_higher_rho[index_nn].astype(int)
    
    deltas[index_rho[0]] = np.max(deltas)   
    return deltas,nearest_neiber
        
# 通过阈值选取 rho与delta都大的点
# 作为聚类中心    
def find_centers_auto(rho,deltas):
    rho_threshold = (np.min(rho) + np.max(rho))/ 2
    delta_threshold  = (np.min(deltas) + np.max(deltas))/ 2
    N = np.shape(rho)[0]
    
    centers = []
    for i in range(N):
        if rho[i]>=rho_threshold and deltas[i]>delta_threshold:
            centers.append(i)
    return np.array(centers)

# 选取 rho与delta乘积较大的点作为
# 聚类中心   
def find_centers_K(rho,deltas,K):
    rho_delta = rho*deltas
    centers = np.argsort(-rho_delta)
    return centers[:K]


def cluster_PD(rho,centers,nearest_neiber):
    K = np.shape(centers)[0]
    if K == 0:
        print("can not find centers")
        return
    
    N = np.shape(rho)[0]
    labs = -1*np.ones(N).astype(int)
    
    # 首先对几个聚类中进行标号
    for i, center in enumerate(centers):
        labs[center] = i
   
    # 将密度从大到小排序
    index_rho = np.argsort(-rho)
    for i, index in enumerate(index_rho):
        # 从密度大的点进行标号
        if labs[index] == -1:
            # 如果没有被标记过
            # 那么聚类标号与距离其最近且密度比其大
            # 的点的标号相同
            labs[index] = labs[int(nearest_neiber[index])]
    return labs
        
def draw_decision(rho,deltas,name="0_decision.jpg"):       
    plt.cla()
    for i in range(np.shape(datas)[0]):
        plt.scatter(rho[i],deltas[i],s=16.,color=(0,0,0))
        plt.annotate(str(i), xy = (rho[i], deltas[i]),xytext = (rho[i], deltas[i]))
        plt.xlabel("rho")
        plt.ylabel("deltas")
    plt.savefig(name)
    plt.show()

def draw_cluster(datas,labs,centers, dic_colors, name="0_cluster.jpg"):     
    plt.cla()
    K = np.shape(centers)[0]
    
    for k in range(K):
        sub_index = np.where(labs == k)
        sub_datas = datas[sub_index]
        # 画数据点
        plt.scatter(sub_datas[:,0],sub_datas[:,1],s=16.,color=dic_colors[k])
        # 画聚类中心
        plt.scatter(datas[centers[k],0],datas[centers[k],1],color="k",marker="+",s = 200.)
    plt.savefig(name)
    plt.show()

    
        
    
    


if __name__== "__main__":
    
    dic_colors = {0:(.8,0,0),1:(0,.8,0),
                  2:(0,0,.8),3:(.8,.8,0),
                  4:(.8,0,.8),5:(0,.8,.8),
                  6:(0,0,0)}
    file_name = "jain"    #spiral jain  flame D31
    # with open(file_name+".txt","r",encoding="utf-8") as f:
    #     lines = f.read().splitlines()
    # lines = [line.split("\t")[:-1] for line in lines]


    with open(r'G:\desktop\PR\iris2.csv ', 'r') as fc:   #dataset   BreastTissue.csv iris2.csv
        reader = csv.reader(fc)
        lines1 = []
        for line in reader:
            lines1.append(line)
    lines = lines1[1:]

    datas = np.array(lines).astype(np.float32)
    # 计算距离矩阵
    dists = getDistanceMatrix(datas)
    # 计算dc
    dc = select_dc(dists)
    print("dc",dc)
    # 计算局部密度 
    rho = get_density(dists,dc,method="Gaussion")
    # 计算密度距离
    deltas, nearest_neiber= get_deltas(dists,rho)
  
    # 绘制密度/距离分布图
    draw_decision(rho,deltas,name=file_name+"_decision.jpg")
    
    # 获取聚类中心点
    centers = find_centers_K(rho,deltas,3)
    # centers = find_centers_auto(rho,deltas)
    print("centers",centers)
    
    labs = cluster_PD(rho,centers,nearest_neiber)
    draw_cluster(datas,labs,centers, dic_colors, name=file_name+"_cluster.jpg")

4.实验

4.1人造数据集

4.2鸢尾花数据集

上方为决策图，下方为聚类结果。

总结

DPC 算法第一个最为突出的贡献就是提出决策图的思想，根据决策图来人工手动选择类簇质心。利用数据点的局部密度与数据点的距离这两个数值建立决策图，构建所有样本对象的 坐标点从而形成一个二维空间，通过对各个数据点局部密度和距离的计算统计，很容易就可以构建完成数据集的决策图。DPC 算法第二个最为突出的贡献就是高效的一步分配策略机制，将剩余的样本点分配到距离最近且局部密度高于自身的样本所在的类簇中。 一般来说，类簇质心选择的是在数据点中局部密度值很大且距离也较大的，从决策图的角度来看，类簇质心位于右上角。