DBSCAB算法介绍（Python3实现）

helloWorld

3610人浏览 · 2024-04-26 13:28:25

helloWorld · 2024-04-26 13:28:25 发布

本文的代码与数据集已上传至github：https://github.com/bitacron/MachineLearning
数据集地址：https://pan.baidu.com/s/1OYoD06nxvG5kcH35C25cXQ?pwd=p599

一、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	7
Jain	373	2	2
Spiral	312	2	3

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

1、Python3代码实现

"""
DBSCAN.py

DBSCAN聚类算法实现
"""

import time
from collections import namedtuple

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.metrics import (
    accuracy_score,
    adjusted_rand_score,
    f1_score,
    normalized_mutual_info_score,
    rand_score,
)
from sklearn.preprocessing import LabelEncoder, StandardScaler

# ==========================================================
# 数据集配置
# ==========================================================

Dataset = namedtuple(
    "Dataset",
    [
        "path",
        "n_clusters",
        "label_col",
        "id_col",
        "header",
        "sep",
        "skip_rows",
        "comment",
        "eps",
        "min_samples",
    ]
)
DATASETS = {
    # UCI数据集
    "iris": Dataset("dataset/uci/iris.data", 3, -1, None, None, ",", 0, None, 0.14, 8),
    "bezdekIris": Dataset("dataset/uci/bezdekIris.data", 3, -1, None, None, ",", 0, None, 0.14, 8),
    "seeds": Dataset("dataset/uci/seeds_dataset.txt", 3, -1, None, None, r"\s+", 0, None, 0.17, 8),
    "glass": Dataset("dataset/uci/glass.data", 6, -1, 0, None, ",", 0, None, 0.45, 5),
    "wine": Dataset("dataset/uci/wine.data", 3, 0, None, None, ",", 0, None, 0.42, 10),
    "wdbc": Dataset("dataset/uci/wdbc.data", 2, 1, 0, None, ",", 0, None, 0.27, 7),

    # 人工合成数据集
    "flame": Dataset("dataset/synthetic/flame.txt", 2, -1, None, None, ",", 0, None, 0.28, 4),
    "jain": Dataset("dataset/synthetic/jain.txt", 2, -1, None, None, ",", 0, None, 0.315, 4),
    "spiral": Dataset("dataset/synthetic/spiral.txt", 3, -1, None, None, ",", 0, None, 0.45, 4),
    "panelB": Dataset("dataset/synthetic/panelB.txt", 3, None, None, None, ",", 0, None, 0.45, 4),
    "panelC": Dataset("dataset/synthetic/panelC.txt", 3, None, None, None, ",", 0, None, 0.45, 4),
    "aggregation": Dataset("dataset/synthetic/aggregation.txt", 7, -1, None, None, ",", 0, None, 0.18, 4),
    "r15": Dataset("dataset/synthetic/R15.txt", 15, -1, None, None, ",", 0, None, None, None),
    "d31": Dataset("dataset/synthetic/D31.txt", 31, -1, None, None, ",", 0, None, None, None),
    "compound": Dataset("dataset/synthetic/compound.txt", 6, -1, None, None, ",", 0, None, None, None),
    "pathbased": Dataset("dataset/synthetic/pathbased.txt", 3, -1, None, None, ",", 0, None, None, None),
}


# ==========================================================
# 数据加载
# ==========================================================

def load_dataset(dataset_name):
    """
    加载数据集
    """
    config = DATASETS[dataset_name]

    print("加载数据集:   " f"数据集={_dataset_name}")

    df = pd.read_csv(config.path, sep=config.sep, header=config.header, skiprows=config.skip_rows, comment=config.comment)

    if config.label_col is not None:
        labels_true = df.iloc[:, config.label_col].to_numpy()
    else:
        labels_true = None

    drop_cols = []

    if config.id_col is not None:
        drop_cols.append(df.columns[config.id_col])

    if config.label_col is not None:
        drop_cols.append(df.columns[config.label_col])

    if drop_cols:
        features = df.drop(columns=drop_cols).values
    else:
        features = df.values

    print(f"样本数量={features.shape[0]}   " f"属性数量(维度)={features.shape[1]}   " f"真实类别数={config.n_clusters}")
    print(f"eps={config.eps}   min_samples={config.min_samples}")
    return features, labels_true, config.eps, config.min_samples


# ==========================================================
# 常量定义
# ==========================================================

UNCLASSIFIED = 0
NOISE = -1


# ==========================================================
# DBSCAN核心算法
# ==========================================================

def get_distance_matrix(features):
    """
    计算欧氏距离矩阵
    """
    n_samples = features.shape[0]
    dists = np.zeros((n_samples, n_samples), dtype=np.float32)

    for i in range(n_samples):
        for j in range(n_samples):
            vi = features[i]
            vj = features[j]
            dists[i, j] = np.sqrt(np.dot(vi - vj, vi - vj))

    return dists


def find_neighbors(point_id, eps, dists):
    """
    寻找以点 point_id 为中心，eps 为半径的圆内的所有点的索引
    """
    indices = np.where(dists[point_id] <= eps)[0]
    return indices.tolist()


def expand_cluster(dists, labels, cluster_id, seeds, eps, min_samples):
    """
    扩展聚类
    """
    i = 0
    while i < len(seeds):
        current_point = seeds[i]

        # 如果该点被标记为NOISE则重新标记为当前聚类
        if labels[current_point] == NOISE:
            labels[current_point] = cluster_id
        # 如果该点未被标记过
        elif labels[current_point] == UNCLASSIFIED:
            # 标记为当前聚类
            labels[current_point] = cluster_id
            # 计算该点的邻域
            new_seeds = find_neighbors(current_point, eps, dists)

            # 如果邻域足够大，则将其加入到seeds队列中
            if len(new_seeds) >= min_samples:
                seeds = seeds + new_seeds

        i += 1

def density_based_clustering(features, eps, min_samples, scale=True):
    """
    DBSCAN聚类算法
    """
    # 数据标准化（默认开启）
    if scale:
        features = StandardScaler().fit_transform(features)
    start = time.time()

    n_samples = features.shape[0]

    # 计算距离矩阵
    dists = get_distance_matrix(features)

    # 初始化标签
    labels = np.full(n_samples, UNCLASSIFIED, dtype=int)

    cluster_id = 0

    # 遍历所有点
    for point_id in range(n_samples):
        # 如果当前点已经处理过，则跳过
        if labels[point_id] != UNCLASSIFIED:
            continue

        # 找到当前点的邻域
        neighbors = find_neighbors(point_id, eps, dists)

        # 如果邻域点数少于min_samples，标记为NOISE
        if len(neighbors) < min_samples:
            labels[point_id] = NOISE
        else:
            # 否则开始一个新的聚类
            cluster_id += 1
            labels[point_id] = cluster_id
            expand_cluster(dists, labels, cluster_id, neighbors, eps, min_samples)

    elapsed_time = time.time() - start
    print(f"耗时={elapsed_time:.4f} s")

    return labels, cluster_id


# ==========================================================
# 标签对齐
# ==========================================================

def align_labels(labels_true, labels_pred, noise_label=None):
    """
    使用匈牙利算法对齐标签。

    参数
    ----------
    labels_true : ndarray
        真实标签。

    labels_pred : ndarray
        预测标签。

    noise_label : int or None
        噪声类别标签。
        如果不为 None，则禁止该类别参与匈牙利匹配。

    返回
    -------
    labels_pred_aligned : ndarray
        对齐后的预测标签。
    """
    from sklearn.metrics import confusion_matrix
    from scipy.optimize import linear_sum_assignment

    labels = np.union1d(labels_true, labels_pred)

    cm = confusion_matrix(
        labels_true,
        labels_pred,
        labels=labels
    )

    row_ind, col_ind = linear_sum_assignment(-cm)

    mapping = {}

    for r, c in zip(row_ind, col_ind):
        pred_label = labels[c]
        true_label = labels[r]

        # 禁止噪声参与匹配
        if noise_label is not None and pred_label == noise_label:
            continue

        mapping[pred_label] = true_label

    labels_pred_aligned = np.array(
        [
            mapping.get(label, label)
            for label in labels_pred
        ]
    )

    return labels_pred_aligned
# ==========================================================
# 聚类评价指标
# ==========================================================

def clustering_indicators(labels_true, labels_pred):
    """
    计算聚类评价指标。

    对于 DBSCAN：
        1. 噪声点参与评价；
        2. 预测为噪声视为错误；
        3. 噪声类别不参与匈牙利匹配。
    """
    if isinstance(labels_true[0], str):
        labels_true = LabelEncoder().fit_transform(labels_true)

    labels_true = labels_true.copy()
    labels_pred = labels_pred.copy()

    noise_label = None

    #
    # 将噪声映射为新的类别
    #
    if np.any(labels_pred == NOISE):
        noise_label = (
            max(
                np.max(labels_true),
                np.max(labels_pred)
            ) + 1
        )

        labels_pred[
            labels_pred == NOISE
        ] = noise_label

    #
    # 标签对齐
    #
    labels_pred_aligned = align_labels(
        labels_true,
        labels_pred,
        noise_label
    )

    #
    # 计算指标
    #
    f_measure = f1_score(
        labels_true,
        labels_pred_aligned,
        average="macro"
    )

    accuracy = accuracy_score(
        labels_true,
        labels_pred_aligned
    )

    normalized_mutual_information = normalized_mutual_info_score(
        labels_true,
        labels_pred
    )

    rand_index = rand_score(
        labels_true,
        labels_pred
    )

    adjusted_rand_index = adjusted_rand_score(
        labels_true,
        labels_pred
    )

    return (
        f_measure,
        accuracy,
        normalized_mutual_information,
        rand_index,
        adjusted_rand_index
    )
# ==========================================================
# 可视化
# ==========================================================
def draw_cluster(features, labels_pred, eps=None, min_samples=None):
    """
    绘制聚类结果
    """
    features = np.asarray(features)

    if features.shape[1] > 2:
        pca = PCA(n_components=2)
        features_2d = pca.fit_transform(features)
    else:
        features_2d = features

    plt.figure(figsize=(8, 6))

    # 分离噪声点和聚类点
    noise_mask = labels_pred == NOISE
    cluster_mask = ~noise_mask

    # 绘制聚类点
    if np.any(cluster_mask):
        # 获取聚类标签（排除噪声）
        cluster_labels = labels_pred[cluster_mask]
        unique_labels = np.unique(cluster_labels)

        for label in unique_labels:
            mask = (labels_pred == label) & cluster_mask
            plt.scatter(features_2d[mask, 0], features_2d[mask, 1], s=7, cmap="nipy_spectral", label=f"Cluster {label}")

    # 绘制噪声点（黑色）
    if np.any(noise_mask):
        plt.scatter(features_2d[noise_mask, 0], features_2d[noise_mask, 1], color="black", s=7, alpha=0.6, label="Noise")
    title = "DBSCAN Clustering Result"
    if eps is not None and min_samples is not None:
        title += f" (eps={eps}, min_samples={min_samples})"
    plt.title(title)
    if np.any(noise_mask):
        plt.legend()
    plt.show()

# ==========================================================
# 主程序
# ==========================================================

if __name__ == "__main__":

    # 选择数据集
    _dataset_name = "spiral"

    _features, _labels_true, _eps, _min_samples = load_dataset(_dataset_name)

    # 执行DBSCAN聚类
    _labels_pred, _num_clusters = density_based_clustering(_features, eps=_eps, min_samples=_min_samples, scale=True)

    print(f"聚类数量={_num_clusters}")

    if _labels_true is not None:
        # 有标签：计算聚类指标
        F1, ACC, NMI, RI, ARI = clustering_indicators(_labels_true, _labels_pred)
        print("聚类指标: "  f"F1={F1:.6f}  "  f"ACC={ACC:.6f}  "  f"NMI={NMI:.6f}  "  f"RI={RI:.6f}  "  f"ARI={ARI:.6f}")

    draw_cluster(_features, _labels_pred, eps=_eps, min_samples=_min_samples)

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 align_labels(labels_true, labels_pred, noise_label=None):
    """
    使用匈牙利算法对齐标签。

    参数
    ----------
    labels_true : ndarray
        真实标签。

    labels_pred : ndarray
        预测标签。

    noise_label : int or None
        噪声类别标签。
        如果不为 None，则禁止该类别参与匈牙利匹配。

    返回
    -------
    labels_pred_aligned : ndarray
        对齐后的预测标签。
    """
    from sklearn.metrics import confusion_matrix
    from scipy.optimize import linear_sum_assignment

    labels = np.union1d(labels_true, labels_pred)

    cm = confusion_matrix(
        labels_true,
        labels_pred,
        labels=labels
    )

    row_ind, col_ind = linear_sum_assignment(-cm)

    mapping = {}

    for r, c in zip(row_ind, col_ind):
        pred_label = labels[c]
        true_label = labels[r]

        # 禁止噪声参与匹配
        if noise_label is not None and pred_label == noise_label:
            continue

        mapping[pred_label] = true_label

    labels_pred_aligned = np.array(
        [
            mapping.get(label, label)
            for label in labels_pred
        ]
    )

    return labels_pred_aligned
# ==========================================================
# 聚类评价指标
# ==========================================================

def clustering_indicators(labels_true, labels_pred):
    """
    计算聚类评价指标。

    对于 DBSCAN：
        1. 噪声点参与评价；
        2. 预测为噪声视为错误；
        3. 噪声类别不参与匈牙利匹配。
    """
    if isinstance(labels_true[0], str):
        labels_true = LabelEncoder().fit_transform(labels_true)

    labels_true = labels_true.copy()
    labels_pred = labels_pred.copy()

    noise_label = None

    #
    # 将噪声映射为新的类别
    #
    if np.any(labels_pred == NOISE):
        noise_label = (
            max(
                np.max(labels_true),
                np.max(labels_pred)
            ) + 1
        )

        labels_pred[
            labels_pred == NOISE
        ] = noise_label

    #
    # 标签对齐
    #
    labels_pred_aligned = align_labels(
        labels_true,
        labels_pred,
        noise_label
    )

    #
    # 计算指标
    #
    f_measure = f1_score(
        labels_true,
        labels_pred_aligned,
        average="macro"
    )

    accuracy = accuracy_score(
        labels_true,
        labels_pred_aligned
    )

    normalized_mutual_information = normalized_mutual_info_score(
        labels_true,
        labels_pred
    )

    rand_index = rand_score(
        labels_true,
        labels_pred
    )

    adjusted_rand_index = adjusted_rand_score(
        labels_true,
        labels_pred
    )

    return (
        f_measure,
        accuracy,
        normalized_mutual_information,
        rand_index,
        adjusted_rand_index
    )