DESCAN聚类 +Python代码实现

前言

​ DESCAN是一种密度聚类算法，其主要思想是密度和传销。本文不做详细的术语讲解。而是通过案例进行讲解，笔者认为这种方式更能通俗易懂。如果想看严谨的原理解释，那本文可能并不适合你，请移步其他博文。

算法流程

​ ①选定一个半径r和最小密度点数MinPoint

​ ②随机选择数据中的一个点记作核心点。

​ ③从该核心点计算到其他点的距离（比如欧氏距离），选出其中距离小于半径r的点。

​ ④如果这些点的数量（包括核心点）大于等于我们设定的MinPoint，我们就说这些点是密度可达点。并将这些密度可达点记作与核心点同一个簇（类别）。

​ ⑤将得到的密度可达点作为新的核心点。

​ ⑥循环③、④、⑤步，直到不再能够找到密度可达点。

​ ⑦第⑥步结束后开启新一轮的循环，即重复②-⑥步，直到所有的点都被设定为核心点，也就是所有的点都计算完。

实现案例

​ 先来看一下图：

​ 上面每个小圆代表数据。我们的目标就是要将这些数据聚成两类。

​ 根据步骤，我们选定一个半径为r=1，Minpoint=3。

​ 第二步随机选取一个点作为核心点，假设选中x1，那么计算距离

​ 可以看到，计算距离之后，对于所有的点，距离小于r的有{x0、x1、x3、x5}。此时的距离小于r的数目为4，因为4>MinPoint。所以此时可以说{x0、x3、x5}（除去x1核心点）是密度可达点。并且{x0、x3、x5}和核心点x1归为一簇，即{x0、x1、x3、x5}为一簇。

​ 接下来就将{x0、x3、x5}作为新的核心点，分别计算其是否密度可达点。最终算出来

​ x0∈{x0、x1、x3、x5}，点数（4）大于MinPoint，所以可以将{x1、x3、x5}与x0归为一簇，因为x0是从x1的密度可达点，所以x0与x1必须是同一类。又因为{x0、x1、x3、x5}包含了x0里的{x0、x1、x3、x5}，所以不做修改。

​ x5∈{x0、x1、x5}，点数（3）等于Minpoint，所以可以将{x0、x1}与x5归为一簇，因为x5是从x1的密度可达点，所以x5与x1必须是同一类。又因为{x0、x1、x3、x5}包含了x5里的{x0、x1、x5}，所以不做修改。

​ 对于x3

​ x3∈{x0、x1、x3、x4}，点数（4）大于MinPoint，所以可以将{x0、x1、x4}与x3归为一簇，因为x3是从x1的密度可达点，所以x3与x1必须是同一类。又因为{x0、x1、x3、x5}没有包含{x4}，所以修改原始的簇为{x0、x1、x3、x4、x5}。

​ x3中包含有新的密度可达点x4。那么就将x4作为新的核心点。计算后结果也没有新的密度可达点。所以将{x0、x1、x3、x4、x5}归为一簇。

​ 接下来再随机选择一个没有被遍历过的点，假设选中x11。

​ 发现其仅仅包含一个数据点。数据点数目小于3，所以不做分类。

​ 因此直接进行下一步。随机选中一个没有被迭代过的点作为核心点。假设选中x10

​ 可以看到点数等于Minpoint，因此将其余作为密度可达点继续迭代。最终得到x{6、7，8，9，10}作为一簇。

​ 最后完成迭代，结束算法。

代码实现

​ 注意两点。

​ ①因为要遍历所有的点，因此代码中用self.is_check代表是否某个点被遍历过。

​ ②代码计算距离的公式用欧氏距离

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import imageio
plt.ion()
a=0
fig,ax=plt.subplots()
image=[]
class DESCAN():
    def __init__(self,x):
        self.x=x
        #假定数据对应的标签y为-1，训练过程会更改
        self.y=np.ones(shape=(self.x.shape[0],1))*-1
        #计算某个点到其他点的距离
        self.distance=np.ones(shape=(self.x.shape[0],self.x.shape[0]))
        for index,i in enumerate(self.x):
            #比如第i行代表第i个数据到其他点的距离。
            self.distance[index,:]=np.linalg.norm(self.x-i,ord=2,axis=1).reshape(-1)
    def fit(self,r,numpoit):
        self.r=r#半径
        self.numpoint=numpoit#最小点数
        # 生成一个矩阵用来记录该点是否被遍历过，事先是定全为1，被遍历过后设置为0.
        self.is_check = np.ones(shape=self.x.shape[0])
        class_num=0#标签值
        #设置循环直到所有的数据都被遍历过
        while True:
            # 计算概率。因为被遍历过的数据点对应的is_check会被置为0，也就是概率为0，不会让下面随机选的时候选到它
            p=self.is_check/sum(self.is_check)
            # 随机选择一个点，返回对应索引，p是用来排除已经被计算过的点
            index=np.random.choice(self.x.shape[0],1,p=p)
            self.is_check[index]=0#令被选中的点更改为已被遍历过
            self.cycle(index,class_num,True)#代入循环式
            if (self.is_check == 0).all():  # 判断所有的数据是否都被遍历过，是则跳出循环
                plt.ioff() #绘图函数，与模型无关
                plt.show() #绘图函数，与模型无关
                break
            # 判断是否标签值y中是否有新的类别。意为如果self.y中只要包含有一个我们传过去的class_num，说明有满足条件的簇.
            #因此我们就要寻找下一个簇。因此class_num+=1
            if (self.y==class_num).any():
                class_num+=1
        # imageio.mimsave("result.gif",image) #保存图像函数，与模型无关
    def cycle(self,index,class_num,First=False):
        '''
        :param index: 核心点的索引
        :param class_num: 类别标签
        :param First: 用于判断是从fit函数过来的，还是从cycle循环本身循环过来的
        :return:
        '''
        #因此index在第二次循环的时候会有可能有多个index，因此作循环。
        for i in index:#迭代选取一个核心点的索引，寻找密度可达点
            # 取出一个核心点
            kernel_point = self.x[i, :]
            # 取出该核心点和其他点的距离。也就是第i行
            distances=self.distance[i,:].squeeze()
            # 找出距离小于r且没有被遍历过的数据点。如果没有没被遍历过的点那就没有意义了。
            new_index = np.argwhere(((distances <= self.r) & (self.is_check == 1)))
            plot_figure(self.x,self.y,kernel_point,self.r)#绘图函数，与模型无关
            # 判断下一步距离小于r的数据点数能否达到我们预设的最少点数
            if new_index.size != 0:  # 判断是否有没有被遍历过的数据点，有的话则将该这些点作为新的核心点。
                # 找出距离小于r的点的对应的索引，然后判断是否设定的最小点数。
                point_index = np.argwhere((distances <= self.r))
                if point_index.shape[0]>=self.numpoint:
                    ############################
                    # 将这些点赋给对应的标签值，因此此处只更新新的没有遍历过的点。但是一开始从fit传过来的点没有赋予标签值。
                    # 因此用First判断，如果是从fit函数过来的，就要赋给标签值进行额外更新。
                    if First:
                        self.y[index, :] = class_num
                    self.y[new_index] = class_num
                    ############################
                    self.is_check[new_index]=0#设定这些核心点为被计算过
                    self.cycle(new_index,class_num)#将这些核心点放到下一轮循环找出其他的直接密度可达点
def plot_figure(x,y,kernel_point,r):
    global a
    '''绘图函数，与模型无关'''
    map_color={-1:"k",0:"r",1:"g",2:"b",3:"y",4:"c"}
    color_initial=[map_color[i] for i in y.squeeze()]
    a+=1
    kernel_point = kernel_point.squeeze()
    plt.cla()
    plt.scatter(x[:,0],x[:,1],c=color_initial)
    circle=plt.Circle((kernel_point[0],kernel_point[1]),radius=r,fill=False)
    plt.axis("equal")
    ax.add_artist(circle)
    # plt.savefig("1.png") #用于保存图像
    # image.append(imageio.v3.imread("1.png")) #用于保存图像
    plt.pause(0.01)
def main():
    #生成第一类数据g
    x1=stats.norm.rvs(0,1,(100,2))
    #生成第二类数据r
    x2=stats.norm.rvs(10,1,(100,2))
    #生成第三类数据b
    x3=stats.norm.rvs(20,1,(100,2))
    #合并数据
    x=np.concatenate([x1,x2,x3],axis=0)
    model=DESCAN(x)
    model.fit(1,3)
if __name__ == '__main__':
    main()

结束语

​ 文章并不严谨。有任何错误的地方。还望指出。阿里嘎多。