DBSCAB算法介绍（Python3实现）

一、DBSCAB算法简介

1、DBSCAN算法

基于密度的空间聚类的应用(Density-based spatial clustering of applications with noise，DBSCAN)算法是由Martin Ester, Hans-Peter Kriegel, Jörg Sander和Xiaowei Xu于1996年提出的一种聚类分析算法。
其原始论文是在1996年的KDD会议（Knowledge Discovery and Data Mining）上发表的，论文名称为"A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise"。
该算法是一种基于密度的聚类算法，该类算法用样本点密度表示相似度，因此这类算法可以适用于任何类型的数据集。

2、DBSCAN算法基本思想

DBSCAN 算法的基本思想是通过密度可达关系获得最大密度相连的一个簇，该算法具有较强的抗噪性，但是需要手动设置最小样本数和邻域半径作为参数，且聚类结果收到参数影响较大，所以需要选取合适的参数。

DBSCAN算法主要有两个参数：
Eps：指定邻域半径；
MinPts：密度阈值。

DBSCAN算法的基本概念定义如下：
Eps邻域(Eps-neighborhood of a point)：给定对象半径Eps内的邻域，用N_Eps(p）表示点p的Eps半径内的点的集合。

核心对象 (core points)：若某个点的Eps邻域内的样本点数（密度）达到算法设定的阈值MinPts，则其为核心点。(即Eps领域内的样本数量不小于MinPts)
直接密度可达(directly density-reachable)：若某点p在点q的Eps邻域内，且q是核心点，则称对象p从核心对象q是直接密度可达。
密度可达density-reachable)：若有一个点的序列q₀,q₁,…,q_k，对任意q_i到q_i+1是直接密度可达的，则称从q₀到q_k密度可达，这实际上是直接密度可达的“传播”。
密度相连(density-connected)：若从某核心点p出发，点q和点k都是密度可达的，则称点q和点k是密度相连的。
边界点：属于某一个类的非核心点,不能发展下线了。
噪声点(noise)：不属于任何一个类簇的点，从任何一个核心点出发都是密度不可达的。

3、DBSCAN算法步骤

Step1：对一个未访问过的点P，先标记它为已访问。
Step2：如果点P的Eps邻域(即以P为中心，Eps为半径的圆)内的数据点大于等于MinPts阈值，则创建一个新的簇C，并把P加入C。
Step3：对P的Eps邻域内的每个点P’，如果P’未被访问，标记P’为已访问，并且如果P’的Eps邻域内有足够多的点，则将这些点也加入到簇C。
Step4：如果P’不属于任何簇，将P’加入到簇C。
Step5：重复Step2-4，直到所有的点都被访问过。

4、DBSCAN算法伪代码

DBSCAN(D, eps, MinPts) {
   C = 0
   for each unvisited point P in dataset D {
      mark P as visited
      NeighborPts = regionQuery(P, eps)
      if sizeof(NeighborPts) < MinPts
         mark P as NOISE
      else {
         C = next cluster
         expandCluster(P, NeighborPts, C, eps, MinPts)
      }
   }
}

expandCluster(P, NeighborPts, C, eps, MinPts) {
   add P to cluster C
   for each point P' in NeighborPts { 
      if P' is not visited {
         mark P' as visited
         NeighborPts' = regionQuery(P', eps)
         if sizeof(NeighborPts') >= MinPts
            NeighborPts = NeighborPts joined with NeighborPts'
      }
      if P' is not yet member of any cluster
         add P' to cluster C
   }
}

regionQuery(P, eps)
   return all points within P's eps-neighborhood (including P)

5、DBSCAN算法时间复杂度分析

DBSCAN 的基本时间复杂度是 O（n * 找出 Eps 邻域中的点所需要的时间），其中 n 是点的个数。在最坏的情况下，时间复杂度是 O（n^2）。然而，在低维空间，有一些数据结构，如 kd 树，可以有效地检索特定点给定距离内的所有点，时间复杂度可以降低到O（nlogn）。

6、DBSCAN算法优缺点

DBSCAN算法优点：

能够处理任何形状的簇；
能够处理噪声和异常值；
不需要提前指定簇的数量。

DBSCAN算法缺点：

对高维数据效果不好；
对于密度不均匀的数据，聚类效果较差；
对参数敏感，需要选择合适的密度参数，如果Eps、MinPts参数选取不当对结果影响较大。

二、DBSCAB算法实现（Python3）

本文使用的数据集为UCI数据集和人工数据集，分别使用鸢尾花数据集Iris、葡萄酒数据集Wine，和spiral 数据集进行测试，本文从UCI官网上将这三个数据集下载下来，并放入和python文件同一个文件夹内即可。同时由于程序需要，将数据集的列的位置做出了略微改动。数据集具体信息如下表：

数据集	样本数	属性维度	类别个数
aggregation	240	2	2
jain	373	2	2
spiral	312	2	3

数据集在我主页资源里有，免积分下载，如果无法下载，可以私信我。

1、Python3代码实现

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler, LabelEncoder
from sklearn.metrics import normalized_mutual_info_score  # NMI
from sklearn.metrics import rand_score  # RI
from sklearn.metrics import accuracy_score  # ACC
from sklearn.metrics import f1_score  # F-measure
from sklearn.metrics import adjusted_rand_score  # ARI
# from sklearn.cluster import DBSCAN
from sklearn.decomposition import PCA
# from sklearn.preprocessing import LabelEncoder

# 数据保存在.csv文件中
# iris = pd.read_csv("dataset/iris.csv")  # 鸢尾花数据集 Iris  class=3
# wine = pd.read_csv("dataset/wine.csv")  # 葡萄酒数据集 Wine  class=3
# seeds = pd.read_csv("dataset/seeds.csv")  # 小麦种子数据集 seeds  class=3
# wdbc = pd.read_csv("dataset/wdbc.csv")  # 威斯康星州乳腺癌数据集 Breast Cancer Wisconsin Diagnostic class=2
# glass = pd.read_csv("dataset/glass.csv")  # 玻璃辨识数据集 Glass Identification  class=6
# aggregation = pd.read_csv("dataset/aggregation.csv")  # 人工数据集 Eps=0.18 MinPts=4
# flame = pd.read_csv("dataset/flame.csv")  # 人工数据集 Eps=0.28 MinPts=4
# jain = pd.read_csv("dataset/jain.csv")  # 人工数据集 Eps = 0.315 MinPts = 4
spiral = pd.read_csv("dataset/spiral.csv")  # 人工数据集 Eps=0.45 MinPts=4
df = spiral  # 设置要读取的数据集
# print(df)
columns = list(df.columns)  # 获取数据集的第一行，第一行通常为特征名，所以先取出
features = columns[:len(columns) - 1]  # 数据集的特征名（去除了最后一列，因为最后一列存放的是标签，不是数据）
dataset = df[features]  # 预处理之后的数据，去除掉了第一行的数据（因为其为特征名，如果数据第一行不是特征名，可跳过这一步）
attributes = len(df.columns) - 1  # 属性数量（数据集维度）
original_labels = list(df[columns[-1]])  # 原始标签

UNCLASSIFIED = 0
NOISE = -1


# 计算数据点两两之间的距离
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


#  寻找以点cluster_id 为中心，eps 为半径的圆内的所有点的id
def find_points_in_eps(point_id, eps, dists):
    index = (dists[point_id] <= eps)
    return np.where(index == True)[0].tolist()


# 聚类扩展
# dists ： 所有数据两两之间的距离  N x N
# labs :   所有数据的标签 labs N，
# cluster_id ： 一个簇的标号
# eps ： 密度评估半径
# seeds： 用来进行簇扩展的点
# min_points： 半径内最少的点数
def expand_cluster(dists, labs, cluster_id, seeds, eps, min_points):
    i = 0
    while i < len(seeds):
        # 获取一个临近点
        Pn = seeds[i]
        # 如果该点被标记为NOISE 则重新标记
        if labs[Pn] == NOISE:
            labs[Pn] = cluster_id
        # 如果该点没有被标记过
        elif labs[Pn] == UNCLASSIFIED:
            # 进行标记，并计算它的临近点 new_seeds
            labs[Pn] = cluster_id
            new_seeds = find_points_in_eps(Pn, eps, dists)

            # 如果 new_seeds 足够长则把它加入到seed 队列中
            if len(new_seeds) >= min_points:
                seeds = seeds + new_seeds

        i = i + 1


def dbscan(datas, Eps, MinPts):
    # 计算 所有点之间的距离
    dists = getDistanceMatrix(datas)

    # 将所有点的标签初始化为UNCLASSIFIED
    n_points = datas.shape[0]
    labs = [UNCLASSIFIED] * n_points

    cluster_id = 0
    # 遍历所有点
    for point_id in range(0, n_points):
        # 如果当前点已经处理过了
        if not (labs[point_id] == UNCLASSIFIED):
            continue

        # 没有处理过则计算临近点
        seeds = find_points_in_eps(point_id, Eps, dists)

        # 如果临近点数量过少则标记为 NOISE
        if len(seeds) < MinPts:
            labs[point_id] = NOISE
        else:
            # 否则就开启一轮簇的扩张
            cluster_id = cluster_id + 1
            # 标记当前点
            labs[point_id] = cluster_id
            expand_cluster(dists, labs, cluster_id, seeds, Eps, MinPts)
    return labs, cluster_id


# 计算聚类指标
def clustering_indicators(labels_true, labels_pred):
    if type(labels_true[0]) != int:
        labels_true = LabelEncoder().fit_transform(df[columns[len(columns) - 1]])  # 如果标签为文本类型，把文本标签转换为数字标签
    f_measure = f1_score(labels_true, labels_pred, average='macro')  # F值
    accuracy = accuracy_score(labels_true, labels_pred)  # ACC
    normalized_mutual_information = normalized_mutual_info_score(labels_true, labels_pred)  # NMI
    rand_index = rand_score(labels_true, labels_pred)  # RI
    ARI = adjusted_rand_score(labels_true, labels_pred)
    return f_measure, accuracy, normalized_mutual_information, rand_index, ARI


# 绘图    
def draw_cluster(datas, labs, n_cluster):
    plt.cla()
    colors = np.array(["red", "blue", "green", "orange", "purple", "cyan", "magenta", "beige", "hotpink", "#88c999"])
    if attributes > 2:
        datas = PCA(n_components=2).fit_transform(datas)  # 如果属性数量大于2，降维
    # plt.scatter(datas[:, 0], datas[:, 1], s=7., color="black")
    # plt.show()
    for i, lab in enumerate(labs):
        if lab == NOISE:
            plt.scatter(datas[i, 0], datas[i, 1], s=7., color=(0, 0, 0))
        else:
            plt.scatter(datas[i, 0], datas[i, 1], s=7., color=colors[lab - 1])

    plt.show()


if __name__ == "__main__":
    # 设置DBSCAN参数
    Eps = 0.45
    MinPts = 4

    datas = np.array(dataset).astype(np.float32)
    # 数据正则化
    datas = StandardScaler().fit_transform(datas)
    label, cluster_id = dbscan(datas, Eps=Eps, MinPts=MinPts)  # 执行DBSCAN
    print(original_labels)  # 输出原始标签
    print("label of my dbscan")
    print(label)  # 输出聚类结果
    F_measure, ACC, NMI, RI, ARI = clustering_indicators(original_labels, label)  # 计算聚类指标
    print("F_measure:", F_measure, "ACC:", ACC, "NMI", NMI, "RI", RI, "ARI", ARI)
    draw_cluster(datas, label, cluster_id)  # 画散点图
    # db = DBSCAN(Eps=Eps, MinPts=MinPts).fit(datas)  # 调用sk的dbscan
    # skl_labels = db.labels_
    # print("label of sk-DBSCAN")
    # print(skl_labels)
    # draw_cluster(db, skl_labels, cluster_id)

2、聚类结果分析

本文选择了F值（F-measure，FM）、准确率（Accuracy，ACC）、标准互信息（Normalized Mutual Information，NMI）和兰德指数（Rand Index，RI）作为评估指标，其值域为[0,1]，取值越大说明聚类结果越符合预期。

F值结合了精度（Precision）与召回率（Recall）两种指标，它的值为精度与召回率的调和平均，其计算公式见公式：

$Precision=\frac{TP}{TP+FP}$

$Recall=\frac{TP}{TP+FN}$

$F-measure=\frac{2Recall \times Precision}{Recall+Precision}$

ACC是被正确分类的样本数与数据集总样本数的比值，计算公式如下：

$ACC=\frac{TP+TN}{TP+TN+FP+FN}$

其中，TP（True Positive）表示将正类预测为正类数的样本个数，TN （True Negative）表示将负类预测为负类数的样本个数，FP（False Positive）表示将负类预测为正类数误报的样本个数，FN（False Negative）表示将正类预测为负类数的样本个数。

NMI用于量化聚类结果和已知类别标签的匹配程度，相比于ACC，NMI的值不会受到族类标签排列的影响。计算公式如下：

$NMI=\frac{I\left(U,V\right)}{\sqrt{H\left(U\right)H\left(V\right)}}$

其中H(U)代表正确分类的熵，H(V)分别代表通过算法得到的结果的熵。

其具体实现代吗如下：
由于数据集中给定的正确标签可能为文本类型而不是数字标签，所以在计算前先判断数据集的标签是否为数字类型，如果不是，则转化为数字类型

def clustering_indicators(labels_true, labels_pred):
    if type(labels_true[0]) != int:
        labels_true = LabelEncoder().fit_transform(df[columns[len(columns) - 1]])  # 如果标签为文本类型，把文本标签转换为数字标签
    f_measure = f1_score(labels_true, labels_pred, average='macro')  # F值
    accuracy = accuracy_score(labels_true, labels_pred)  # ACC
    normalized_mutual_information = normalized_mutual_info_score(labels_true, labels_pred)  # NMI
    rand_index = rand_score(labels_true, labels_pred)  # RI
    return f_measure, accuracy, normalized_mutual_information, rand_index


F_measure, ACC, NMI, RI = clustering_indicators(labels_number, labels)
print("F_measure:", F_measure, "ACC:", ACC, "NMI", NMI, "RI", RI)