聚类算法之DBSCAN

K-Means算法 和Mean Shift算法都是基于距离的聚类算法，基于距离的聚类算法的聚类结果是球状的簇，当数据集中的聚类结果是非球状结构时，基于距离的聚类算法的聚类效果并不好。

与基于距离的聚类算法不同的是，基于密度的聚类算法可以发现任意形状的聚类。在基于密度的聚类算法中，通过在数据集中寻找被低密度区域分离的高密度区域，将分离出的高密度区域作为一个独立的类别。DBSCAN（Density-Based Spatial Clustering of Application with Noise）是一种典型的基于密度的聚类算法。

DBSCAN算法原理

DBSCAN（Density-Based Spatial Clustering of Applications with Noise，具有噪声的基于密度的聚类方法）是一种基于密度的空间聚类算法。该算法将具有足够密度的区域划分为簇，并在具有噪声的空间数据库中发现任意形状的簇，它将簇定义为密度相连的点的最大集合。

在DBSCAN算法中将数据点分为三类：

• 核心点（Core point）。若样本的邻域内至少包含了MinPts个样本，即 ，则称样本点为核心点。
• 边界点（Border point）。若样本的邻域内包含的样本数目小于MinPts，但是它在其他核心点的邻域内，则称样本点为边界点。
• 噪音点（Noise）。既不是核心点也不是边界点的点

在这里有两个量，一个是半径Eps（），另一个是指定的数目MinPts。

在DBSCAN算法中，还定义了如下一些概念：

• 密度直达(directly density-reachable)：我们称样本点 p 是由样本点 q 对于参数 {Eps,MinPts} 密度直达的，如果它们满足 p∈NEps(q) 且 |NEps(q)|≥MinPts （即样本点 q 是核心点）
• 密度可达(density-reachable)：我们称样本点 p 是由样本点 q 对于参数{Eps,MinPts}密度可达的，如果存在一系列的样本点 p1,…,pn（其中 p1=q,pn=p）使得对于i=1,…,n−1，样本点 pi+1 可由样本点 pi 密度可达
• 密度相连(density-connected)：我们称样本点 p 与样本点 q 对于参数 {Eps,MinPts} 是密度相连的，如果存在一个样本点 o，使得 p 和 q 均由样本点 o 密度可达。

基于密度的聚类算法通过寻找被低密度区域分离的高密度区域，并将高密度区域作为一个聚类的“簇”。在DBSCAN算法中，聚类“簇”定义为：由密度可达关系导出的最大的密度连接样本的集合。

DBSCAN算法流程

在DBSCAN算法中，有核心对象出发，找到与该核心对象密度可达的所有样本形成“簇”。DBSCAN算法的流程为：

• 根据给定的邻域参数Eps和MinPts确定所有的核心对象
• 对每一个核心对象
  • 选择一个未处理过的核心对象，找到由其密度可达的的样本生成聚类“簇”
• 重复以上过程

伪代码：

(1) 首先将数据集D中的所有对象标记为未处理状态  
(2) for（数据集D中每个对象p） do  
(3)    if （p已经归入某个簇或标记为噪声） then  
(4)         continue;  
(5)    else  
(6)         检查对象p的Eps邻域 NEps(p) ；  
(7)         if (NEps(p)包含的对象数小于MinPts) then  
(8)                  标记对象p为边界点或噪声点；  
(9)         else  
(10)                 标记对象p为核心点，并建立新簇C, 并将p邻域内所有点加入C  
(11)                 for (NEps(p)中所有尚未被处理的对象q)  do  
(12)                       检查其Eps邻域NEps(q)，若NEps(q)包含至少MinPts个对象，则将NEps(q)中未归入任何一个簇的对象加入C；  
(13)                 end for  
(14)        end if  
(15)    end if  
(16) end for

Python实现：

# -*- coding: utf-8 -*-
import numpy as np
 
 
def distance(data):
    '''计算样本点之间的距离
    :param data(mat):样本
    :return:dis(mat):样本点之间的距离
    '''
    m, n = np.shape(data)
    dis = np.mat(np.zeros((m, m)))
    for i in range(m):
        for j in range(i, m):
            # 计算i和j之间的欧式距离
            tmp = 0
            for k in range(n):
                tmp += (data[i, k] - data[j, k]) * (data[i, k] - data[j, k])
            dis[i, j] = np.sqrt(tmp)
            dis[j, i] = dis[i, j]
    return dis
 
 
def find_eps(distance_D, eps):
    '''找到距离≤eps的样本的索引
    :param distance_D(mat):样本i与其他样本之间的距离
    :param eps(float):半径的大小
    :return: ind(list):与样本i之间的距离≤eps的样本的索引
    '''
    ind = []
    n = np.shape(distance_D)[1]
    for j in range(n):
        if distance_D[0, j] <= eps:
            ind.append(j)
    return ind
 
 
def dbscan(data, eps, MinPts):
    '''DBSCAN算法
    :param data(mat):需要聚类的数据集
    :param eps(float):半径
    :param MinPts(int):半径内最少的数据点数
    :return:
        types(mat):每个样本的类型：核心点、边界点、噪音点
        sub_class(mat):每个样本所属的类别
    '''
    m = np.shape(data)[0]
    # 在types中，1为核心点，0为边界点，-1为噪音点
    types = np.mat(np.zeros((1, m)))
    sub_class = np.mat(np.zeros((1, m)))
    # 用于判断该点是否处理过，0表示未处理过
    dealt = np.mat(np.zeros((m, 1)))
    # 计算每个数据点之间的距离
    dis = distance(data)
    # 用于标记类别
    number = 1
 
    # 对每一个点进行处理
    for i in range(m):
        # 找到未处理的点
        if dealt[i, 0] == 0:
            # 找到第i个点到其他所有点的距离
            D = dis[i,]
            # 找到半径eps内的所有点
            ind = find_eps(D, eps)
            # 区分点的类型
            # 边界点
            if len(ind) > 1 and len(ind) < MinPts + 1:
                types[0, i] = 0
                sub_class[0, i] = 0
            # 噪音点
            if len(ind) == 1:
                types[0, i] = -1
                sub_class[0, i] = -1
                dealt[i, 0] = 1
            # 核心点
            if len(ind) >= MinPts + 1:
                types[0, i] = 1
                for x in ind:
                    sub_class[0, x] = number
                # 判断核心点是否密度可达
                while len(ind) > 0:
                    dealt[ind[0], 0] = 1
                    D = dis[ind[0],]
                    tmp = ind[0]
                    del ind[0]
                    ind_1 = find_eps(D, eps)
 
                    if len(ind_1) > 1:  # 处理非噪音点
                        for x1 in ind_1:
                            sub_class[0, x1] = number
                        if len(ind_1) >= MinPts + 1:
                            types[0, tmp] = 1
                        else:
                            types[0, tmp] = 0
 
                        for j in range(len(ind_1)):
                            if dealt[ind_1[j], 0] == 0:
                                dealt[ind_1[j], 0] = 1
                                ind.append(ind_1[j])
                                sub_class[0, ind_1[j]] = number
                number += 1
 
    # 最后处理所有未分类的点为噪音点
    ind_2 = ((sub_class == 0).nonzero())[1]
    for x in ind_2:
        sub_class[0, x] = -1
        types[0, x] = -1
 
    return types, sub_class

DBSCAN的参数选择

MinPts

这个参数建议根据数据量及具体的业务进行自行设定

Eps

《Python机器学习算法》这本书上给出了一个计算公式，但是没有解释中间的原因，并不清楚理论依据是什么，算法如下：

def epsilon(data, MinPts):
    '''计算最佳半径
    input:  data(mat):训练数据
            MinPts(int):半径内的数据点的个数
    output: eps(float):半径
    '''
    m, n = np.shape(data)
    xMax = np.max(data, 0)
    xMin = np.min(data, 0)
    eps = ((np.prod(xMax - xMin) * MinPts * math.gamma(0.5 * n + 1)) / (m * math.sqrt(math.pi ** n))) ** (1.0 / n)
    return eps

其他参考资料：

• https://stackoverflow.com/questions/12893492/choosing-eps-and-minpts-for-dbscan-r
• https://blog.csdn.net/dwf_android/article/details/78059260
• https://en.wikipedia.org/wiki/OPTICS_algorithm

Scikit-learn 中的 DBSCAN 的使用

主要函数介绍：

DBSCAN(eps=0.5, min_samples=5, metric='euclidean', algorithm='auto', leaf_size=30, p=None, n_jobs=1)

核心参数：

• eps: float，ϵ-邻域的距离阈值
• min_samples ：int，样本点要成为核心对象所需要的 ϵ-邻域的样本数阈值

其他参数：

• metric ：度量方式，默认为欧式距离，可以使用的距离度量参数有：
  • 欧式距离 “euclidean”
  • 曼哈顿距离 “manhattan”
  • 切比雪夫距离“chebyshev”
  • 闵可夫斯基距离 “minkowski”
  • 带权重闵可夫斯基距离 “wminkowski”
  • 标准化欧式距离 “seuclidean”
  • 马氏距离“mahalanobis”
  • 自己定义距离函数
• algorithm：近邻算法求解方式，有四种：
  • “brute”蛮力实现
  • “kd_tree” KD树实现
  • “ball_tree”球树实现
  • “auto”上面三种算法中做权衡，选择一个拟合最好的最优算法。
• leaf_size：使用“ball_tree”或“kd_tree”时,停止建子树的叶子节点数量的阈值
• p：只用于闵可夫斯基距离和带权重闵可夫斯基距离中p值的选择，p=1为曼哈顿距离， p=2为欧式距离。如果使用默认的欧式距离不需要管这个参数。
• n_jobs ：CPU并行数，若值为 -1，则用所有的CPU进行运算

属性：

• core_sample_indices_ : 核心点的索引，因为labels_不能区分核心点还是边界点，所以需要用这个索引确定核心点
• components_：训练样本的核心点
• labels_：每个点所属集群的标签，-1代表噪声点

使用示例：

import numpy as np
from sklearn.cluster import DBSCAN
import matplotlib.pyplot as plt
from math import radians, sin, cos, asin, sqrt
 
 
def haversine(latlon1, latlon2):
    """
    计算两经纬度之间的距离
    """
    if (latlon1 - latlon2).all():
        lat1, lon1 = latlon1
        lat2, lon2 = latlon2
        lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
        dlon = lon2 - lon1
        dlat = lat2 - lat1
        a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
        c = 2 * asin(sqrt(a))
        r = 6370996.81  # 地球半径
        distance = c * r
    else:
        distance = 0
    return distance
 
 
if __name__ == "__main__":
    data = []
    f = open("k_means_sample_data.txt", 'r')
    for line in f:
        data.append([float(line.split(',')[0]), float(line.split(',')[1])])
    data = np.array(data)
    MinPts = int(data.shape[0] / 100)
    eps = 2000
 
    db = DBSCAN(eps=eps, min_samples=MinPts, metric=haversine).fit(data)
    core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
    core_samples_mask[db.core_sample_indices_] = True
    labels = db.labels_
    n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)
 
    unique_labels = set(labels)
    colors = ['r', 'b', 'g', 'y', 'c', 'm', 'orange']
    for k, col in zip(unique_labels, colors):
        if k == -1:
            col = 'k'
        class_member_mask = (labels == k)
        xy = data[class_member_mask & core_samples_mask]
        plt.plot(xy[:, 0], xy[:, 1], 'o', markerfacecolor=col, markeredgecolor='w', markersize=10)
 
        xy = data[class_member_mask & ~core_samples_mask]
        plt.plot(xy[:, 0], xy[:, 1], 'o', markerfacecolor=col, markeredgecolor='w', markersize=3)
 
    plt.title('Estimated number of clusters: %d' % n_clusters_)
    plt.show()

执行结果：

DBSCAN优缺点总结

优点：

• 相比K-Means，DBSCAN 不需要预先声明聚类数量。
• 可以对任意形状的稠密数据集进行聚类，相对的，K-Means之类的聚类算法一般只适用于凸数据集。
• 可以在聚类的同时发现异常点，对数据集中的异常点不敏感。
• 聚类结果没有偏倚，相对的，K-Means之类的聚类算法初始值对聚类结果有很大影响。

缺点：

• 当空间聚类的密度不均匀、聚类间距差相差很大时，聚类质量较差，因为这种情况下参数MinPts和Eps选取困难。
• 如果样本集较大时，聚类收敛时间较长，此时可以对搜索最近邻时建立的KD树或者球树进行规模限制来改进。
• 在两个聚类交界边缘的点会视乎它在数据库的次序决定加入哪个聚类，幸运地，这种情况并不常见，而且对整体的聚类结果影响不大（DBSCAN*变种算法，把交界点视为噪音，达到完全决定性的结果。）
• 调参相对于传统的K-Means之类的聚类算法稍复杂，主要需要对距离阈值eps，邻域样本数阈值MinPts联合调参，不同的参数组合对最后的聚类效果有较大影响。