欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

算法模型---聚类分析之DBSCAN密度聚类算法

程序员文章站 2024-01-16 22:51:34
...

来源
DBSCAN(Density-Based Spatial Clustering of Applications with Noise,基于密度的抗噪聚类方法)。和K-Means,BIRCH这些一般只适用于凸样本集的聚类相比,DBSCAN既可以适用于凸样本集,也可以适用于非凸样本集。

1. 密度聚类原理

DBSCAN是一种基于密度的聚类算法,这类密度聚类算法一般假定类别可以通过样本分布的紧密程度决定。同一类别的样本,他们之间紧密相连,也就是说,在该类别任意样本周围不远处一定有同类别的样本存在。
通过将紧密相连的样本划为一类,这样就得到了一个聚类类别。通过将所有各组紧密相连的样本划为各个不同的类别,我们就得到了最终的所有聚类类别结果。

2. DBSCAN密度定义

DBSCAN是基于一组邻域来描述样本集的紧密程度的,参数(ϵ, MinPts)用来描述邻域的样本分布紧密程度。其中,ϵ描述了某一样本邻域的距离阈值(与该样本的距离至少要小于该值才能算该样本领域中的样本),MinPts描述了某一样本的距离为ϵ的邻域中样本个数的阈值(与该样本距离小于ϵ的样本个数超过MinPts才可能构成一个类)。
假设我的样本集是D=(x1,x2,,xm),则DBSCAN具体的密度描述定义如下:
1)ϵ-邻域:对于xjD,其ϵ-邻域包含样本集D中与xj的距离不大于ϵ的子样本集,即Nϵ(xj)={xiD|distance(xi,xj)ϵ}, 这个子样本集的个数记为|Nϵ(xj)|
2) 核心对象:对于任一样本xjD,如果其ϵ-邻域对应的Nϵ(xj)至少包含MinPts个样本,即如果|Nϵ(xj)|MinPts,则xj是核心对象。
3)密度直达:如果xi位于xjϵ-邻域中,且xj是核心对象,则称xixj密度直达。注意反之不一定成立,即此时不能说xjxi密度直达, 除非且xi也是核心对象。
4)密度可达:对于xixj,如果存在样本序列p1,p2,...,pT,满足p1=xi,pT=xj, 且pt+1pt,式中1tT1密度直达,则称xjxi密度可达。也就是说,密度可达满足传递性。此时序列中的传递样本p1,p2,...,pT1均为核心对象,因为只有核心对象才能使其他样本密度直达。密度可达可以理解为本来直达不了,但通过一系列中间点架桥的方式到达了,而且是由核心点可达另一个点。发起方一定是核心对象。
注意密度可达也不满足对称性,这个可以由密度直达的不对称性得出。
5)密度相连:对于xixj,如果存在核心对象样本xk,使xixj均由xk密度可达,则称xixj密度相连。注意密度相连关系是满足对称性的。密度相连可以理解为密度可达的升级版,xi通过一系列的核心对象密度可达,xj也可以通过同一系列对象密度可达,这样就把xixj连起来了。

从下图可以很容易看出上述定义,图中MinPts=5,红色的点都是核心对象,因为其ϵ-邻域至少有5个样本。黑色的样本是非核心对象。所有核心对象密度直达的样本在以红色核心对象为中心的超球体内,如果不在超球体内,则不能密度直达。图中用绿色箭头连起来的核心对象组成了密度可达的样本序列。在这些密度可达的样本序列的ϵ-邻域内所有的样本相互都是密度相连的。
算法模型---聚类分析之DBSCAN密度聚类算法

3. DBSCAN密度聚类思想

DBSCAN的聚类定义很简单:由密度可达关系导出的最大密度相连的样本集合,即为我们最终聚类的一个类别,或者说一个簇。
这个DBSCAN的簇里面可以有一个或者多个核心对象。如果只有一个核心对象,则簇里其他的非核心对象样本都在这个核心对象的ϵ-邻域里;如果有多个核心对象,则簇里的任意一个核心对象的ϵ-邻域中一定有一个其他的核心对象,否则这两个核心对象无法密度可达。这些核心对象的ϵ-邻域里所有的样本的集合组成的一个DBSCAN聚类簇。
那么怎么才能找到这样的簇样本集合呢?DBSCAN使用的方法很简单,它任意选择一个没有类别的核心对象作为种子,然后找到所有这个核心对象能够密度可达的样本集合,即为一个聚类簇。接着继续选择另一个没有类别的核心对象去寻找密度可达的样本集合,这样就得到另一个聚类簇。一直运行到所有核心对象都有类别为止。

基本上这就是DBSCAN算法的主要内容了,是不是很简单?但是我们还是有三个问题没有考虑。
第一个是一些异常样本点或者说少量游离于簇外的样本点,这些点不在任何一个核心对象在周围,在DBSCAN中,我们一般将这些样本点标记为噪音点。
第二个是距离的度量问题,即如何计算某样本和核心对象样本的距离。在DBSCAN中,一般采用最近邻思想,采用某一种距离度量来衡量样本距离,比如欧式距离。这和KNN分类算法的最近邻思想完全相同。对应少量的样本,寻找最近邻可以直接去计算所有样本的距离,如果样本量较大,则一般采用KD树或者球树来快速的搜索最近邻。最近邻的思想,距离度量,KD树和球树可参考原文作者的另一篇文章《K近邻法(KNN)原理小结》
第三种问题比较特殊,某些样本可能到两个核心对象的距离都小于ϵ,但是这两个核心对象由于不是密度直达,又不属于同一个聚类簇,那么如果界定这个样本的类别呢?一般来说,此时DBSCAN采用先来后到,先进行聚类的类别簇会标记这个样本为它的类别。也就是说BDSCAN的算法不是完全稳定的算法。

4. DBSCAN聚类算法

下面我们对DBSCAN聚类算法的流程做一个总结。
输入:样本集D=(x1,x2,...,xm),邻域参数(ϵ,MinPts), 样本距离度量方式
输出: 簇划分C. 
1)初始化核心对象集合Ω=, 初始化聚类簇数k=0,初始化未访问样本集合Γ=D, 簇划分C=(初始化的时候是最坏的结果:没有核心对象,没有聚类簇)
2) 对于j=1,2,...m, 按下面的步骤找出所有的核心对象:
a) 通过距离度量方式,找到样本xjϵ-邻域子样本集Nϵ(xj)
b) 如果子样本集样本个数满足|Nϵ(xj)|MinPts, 将样本xj加入核心对象样本集合:Ω=Ω{xj}
3)如果核心对象集合Ω=,则算法结束,否则转入步骤4.(没有核心对象无法聚类)
4)在核心对象集合Ω中,随机选择一个核心对象o,初始化当前簇核心对象队列Ωcur=o, 初始化类别序号k=k+1,初始化当前簇样本集合Ck={o}, 更新未访问样本集合Γ=Γ{o}。(可以从任务一个核心对象开始进行簇集合的搜索)
5)如果当前簇核心对象队列Ωcur=,则当前聚类簇Ck生成完毕, 更新簇划分C=C1,C2,...,Ck, 更新核心对象集合Ω=ΩCk, 转入步骤3。
6)在当前簇核心对象队列Ωcur中取出一个核心对象o,通过邻域距离阈值ϵ找出所有的ϵ-邻域子样本集Nϵ(o),令Δ=Nϵ(o)Γ, 更新当前簇样本集合Ck=CkΔ, 更新未访问样本集合Γ=ΓΔ, 更新Ωcur=Ωcur(Nϵ(o)Ω),转入步骤5.
输出结果为: 簇划分C={C1,C2,...,Ck}

5. DBSCAN小结

DBSCAN的主要优点有:
1) 可以对任意形状的稠密数据集进行聚类,相对的,K-Means之类的聚类算法一般只适用于凸数据集。
2) 可以在聚类的同时发现异常点,对数据集中的异常点不敏感。
3) 聚类结果没有偏倚,相对的,K-Means之类的聚类算法初始值对聚类结果有很大影响。

DBSCAN的主要缺点有:
1)如果样本集的密度不均匀、聚类间距差相差很大时,聚类质量较差,这时用DBSCAN聚类一般不适合。
2) 如果样本集较大时,聚类收敛时间较长,此时可以对搜索最近邻时建立的KD树或者球树进行规模限制来改进。
3) 调参相对于传统的K-Means之类的聚类算法稍复杂,主要需要对距离阈值ϵ,邻域样本数阈值MinPts联合调参,不同的参数组合对最后的聚类效果有较大影响。

6. DBSCAN Python实现

来源
计算过程

Created with Raphaël 2.1.2开始读入数据找一未分类点扩散是否还有未分类数据输出结果结束yesno

样本样例

/* 788points.txt */
15.55,28.65
14.9,27.55
14.45,28.35
14.15,28.8
13.75,28.05
13.35,28.45
13,29.15
13.45,27.5
13.6,26.5
12.8,27.35
12.4,27.85
12.3,28.4
12.2,28.65
13.4,25.1
12.95,25.95

代码实现

# -*- coding: utf-8 -*-
__author__ = 'Wsine'

import numpy as np
import matplotlib.pyplot as plt
import math
import time
import sys

UNCLASSIFIED = False
NOISE = 0

def loadDataSet(fileName, splitChar='\t'):
    """
    输入:文件名
    输出:数据集
    描述:从文件读入数据集
    """
    dataSet = []
    with open(fileName) as fr:
        for line in fr.readlines():
            curline = line.strip().split(splitChar)
            fltline = list(map(float, curline))
            dataSet.append(fltline)
    return dataSet

def dist(a, b):
    """
    输入:向量A, 向量B
    输出:两个向量的欧式距离
    """
    return math.sqrt(np.power(a - b, 2).sum())

def eps_neighbor(a, b, eps):
    """
    输入:向量A, 向量B
    输出:是否在eps范围内
    """
    return dist(a, b) < eps

def region_query(data, pointId, eps):
    """
    输入:数据集, 查询点id, 半径大小
    输出:查询点eps范围内的点的其他点id(包括查询点本身)
    """
    nPoints = data.shape[1]
    seeds = []
    for i in range(nPoints):
        if eps_neighbor(data[:, pointId], data[:, i], eps):#如果距离小于eps即被记录
            if i !=pointId:
                seeds.append(i)
    return seeds

def expand_cluster(data, clusterResult, pointId, clusterId, eps, minPts):
    """
    输入:数据集, 分类结果, 待分类点id, 簇id, 半径大小, 最小点个数
    输出:能否成功分类
    #这个核心函数中,会判断某个点是不同是否是核心对象,如果不是,暂时将其判断噪声点,可能会误判,但会在其他点的判断中得到纠正
    #如果是核心对象,则会以此点为基础生成一个聚类,并将其周围eps距离内的点标识为同一类;在此基础上,寻找该核心对象eps距离内的其他核心对象,将另一核心对象及另一核心对象周围点划分原始
    #核心对象同一类,并不判断扩展,直至找不到核心对象
    """
    seeds = region_query(data, pointId, eps)
    if len(seeds) < minPts: # 不满足minPts条件的为噪声点(应该是非核心对象)
        clusterResult[pointId] = NOISE#某点不是核心对象,暂时判其为噪声点(类别用0来表示),但如果该点虽然自己不是核心对象,但在其他点判断时,如果其他点是核心对象,
        # 而它又在另一点的eps距离内,它仍然会被重新分到另一类中,因而这里不用担心被误判
        return False
    else:
        clusterResult[pointId] = clusterId # 划分到该簇(由核心对象来代表该簇)
        for seedId in seeds:
            clusterResult[seedId] = clusterId#将周围的点一同划分到该簇

        while len(seeds) > 0: # 通过判断周围的点是否为核心对象,持续扩张
            currentPoint = seeds[0]
            if clusterResult[currentPoint]!=0:
                queryResults = region_query(data, currentPoint, eps)
                #这里可以优化,因为如果之前已经判断为非核心对象则对应clusterResult为0,没必要再算一次,从来没判断过的其实是False,两者还是有区别的
                #同是在seed中应该排除原始核心对象,否则在扩展时会重复算;同时利他其他核心对象扩展时,也就排除其他核心对象。修改的方法应该是region_query函数的返回值就不应该包括查询点本身。
                if len(queryResults) >= minPts:#eps距离内的某一点是核心对象
                    for i in range(len(queryResults)):
                        resultPoint = queryResults[i]
                        if clusterResult[resultPoint] == UNCLASSIFIED:#因为两个距离在eps内的核心对象会将彼此周围的点连成一片,即距离可达,因而这些点被判断为原始核心对象同一类
                            #并且将这些点也归到原始核心对象的周围点中去,从而实现不同扩展,这是整个程序中最巧妙和关键的地方。
                            seeds.append(resultPoint)
                            clusterResult[resultPoint] = clusterId
                        elif clusterResult[resultPoint] == NOISE:#将另一核心对象周围的非核心对象标识为原始核心对象同一类,
                            clusterResult[resultPoint] = clusterId
            seeds = seeds[1:]#不断更新,已经判断过点会被移除,直至seed为空;即原始核心对象eps距离内的点都进行了是否为核心对象的判断。
        return True

def dbscan(data, eps, minPts):
    """
    输入:数据集, 半径大小, 最小点个数
    输出:每个样本的类别,和总的类别数
    """
    clusterId = 1#类别号从1开始,用什么来标识类型本来是无所谓的,但如果从1开始,得到聚类结果的同时也就得到簇的个数。
    nPoints = data.shape[1]
    clusterResult = [UNCLASSIFIED] * nPoints#初始每个样本的类别为False,即还未参考聚类过程,最终聚类后的结果为类别标识,即上面的clusterId
    for pointId in range(nPoints):#逐个点进行类别判定
        point = data[:, pointId]
        if clusterResult[pointId] == UNCLASSIFIED:#如果还未被聚到某一类(虽然是逐个点判断,但只要之前已经生成了该点类别就会跳过去)
            if expand_cluster(data, clusterResult, pointId, clusterId, eps, minPts):#这里没有定义全局变量,expand_cluster也只返回True和False但clusterResult却得到了更新
                clusterId = clusterId + 1#以某一点为核心对象进行扩展,连成一片,之后就是一个类别自然就加1了。
    return clusterResult, clusterId - 1

def plotFeature(data, clusters, clusterNum):
    nPoints = data.shape[1]
    matClusters = np.mat(clusters).transpose()
    fig = plt.figure()
    scatterColors = ['black', 'blue', 'green', 'yellow', 'red', 'purple', 'orange', 'brown']
    ax = fig.add_subplot(111)
    for i in range(clusterNum + 1):
        colorSytle = scatterColors[i % len(scatterColors)]#如果类别很多,就可能存在颜色重复
        subCluster = data[:, np.nonzero(matClusters[:, 0].A == i)]
        ax.scatter(subCluster[0, :].flatten().A[0], subCluster[1, :].flatten().A[0], c=colorSytle, s=50)

def main():
    filePath="%s/788points.txt"%(sys.path[0])
    dataSet = loadDataSet(filePath, splitChar=',')
    dataSet = np.mat(dataSet).transpose()#将数据变成一列一个样本
    # print(dataSet)
    clusters, clusterNum = dbscan(dataSet, 2, 14)
    print("cluster Numbers = ", clusterNum)
    # print(clusters)
    plotFeature(dataSet, clusters, clusterNum)


if __name__ == '__main__':
    start = time.clock()
    main()
    end = time.clock()
    print('finish all in %s' % str(end - start))
    plt.show()

结果如下:
算法模型---聚类分析之DBSCAN密度聚类算法

7、scikit-learn中的DBSCAN类

来源
官网
在scikit-learn中,DBSCAN算法类为sklearn.cluster.DBSCAN。要熟练的掌握用DBSCAN类来聚类,除了对DBSCAN本身的原理有较深的理解以外,还要对最近邻的思想有一定的理解。集合这两者,就可以玩转DBSCAN了。

7.1、DBSCAN类重要参数

DBSCAN类的重要参数也分为两类,一类是DBSCAN算法本身的参数,一类是最近邻度量的参数,下面我们对这些参数做一个总结。

  • eps: DBSCAN算法参数,即我们的ϵ-邻域的距离阈值,和样本距离超过ϵ的样本点不在ϵ-邻域内。默认值是0.5。一般需要通过在多组值里面选择一个合适的阈值。eps过大,则更多的点会落在核心对象的ϵ-邻域,此时我们的类别数可能会减少, 本来不应该是一类的样本也会被划为一类。反之则类别数可能会增大,本来是一类的样本却被划分开。

  • min_samples: DBSCAN算法参数,即样本点要成为核心对象所需要的ϵ-邻域的样本数阈值。默认值是5. 一般需要通过在多组值里面选择一个合适的阈值。通常和eps一起调参。在eps一定的情况下,min_samples过大,则核心对象会过少,此时簇内部分本来是一类的样本可能会被标为噪音点,类别数也会变多。反之min_samples过小的话,则会产生大量的核心对象,可能会导致类别数过少。

  • metric:最近邻距离度量参数。可以使用的距离度量较多,一般来说DBSCAN使用默认的欧式距离(即p=2的闵可夫斯基距离)就可以满足我们的需求。可以使用的距离度量参数有:

    • a) 欧式距离 “euclidean”: i=1n(xiyi)2
    • b) 曼哈顿距离 “manhattan”:i=1n|xiyi|
    • c) 切比雪夫距离“chebyshev”:max|xiyi|(i=1,2,...n)
    • d) 闵可夫斯基距离 “minkowski”: i=1n|xiyi|pp,p=1为曼哈顿距离, p=2为欧式距离。
    • e) 带权重闵可夫斯基距离 “wminkowski”: i=1nw|xiyi|pp,其中w为特征权重
    • f) 标准化欧式距离 “seuclidean”: 即对于各特征维度做了归一化以后的欧式距离。此时各样本特征维度的均值为0,方差为1.
    • g) 马氏距离“mahalanobis”:(xy)TS1(xy),其中,S1为样本协方差矩阵的逆矩阵。当样本分布独立时, S为单位矩阵,此时马氏距离等同于欧式距离。
    • 还有一些其他不是实数的距离度量,一般在DBSCAN算法用不上,这里也就不列了。
  • algorithm:最近邻搜索算法参数,算法一共有三种,第一种是蛮力实现,第二种是KD树实现,第三种是球树实现。这三种方法在K近邻法(KNN)原理小结中都有讲述,如果不熟悉可以去复习下。对于这个参数,一共有4种可选输入,‘brute’对应第一种蛮力实现,‘kd_tree’对应第二种KD树实现,‘ball_tree’对应第三种的球树实现, ‘auto’则会在上面三种算法中做权衡,选择一个拟合最好的最优算法。需要注意的是,如果输入样本特征是稀疏的时候,无论我们选择哪种算法,最后scikit-learn都会去用蛮力实现‘brute’。个人的经验,一般情况使用默认的 ‘auto’就够了。 如果数据量很大或者特征也很多,用”auto”建树时间可能会很长,效率不高,建议选择KD树实现‘kd_tree’,此时如果发现‘kd_tree’速度比较慢或者已经知道样本分布不是很均匀时,可以尝试用‘ball_tree’。而如果输入样本是稀疏的,无论你选择哪个算法最后实际运行的都是‘brute’。

  • leaf_size:最近邻搜索算法参数,为使用KD树或者球树时, 停止建子树的叶子节点数量的阈值。这个值越小,则生成的KD树或者球树就越大,层数越深,建树时间越长,反之,则生成的KD树或者球树会小,层数较浅,建树时间较短。默认是30. 因为这个值一般只影响算法的运行速度和使用内存大小,因此一般情况下可以不管它。
  • p: 最近邻距离度量参数。只用于闵可夫斯基距离和带权重闵可夫斯基距离中p值的选择,p=1为曼哈顿距离, p=2为欧式距离。如果使用默认的欧式距离不需要管这个参数。

以上就是DBSCAN类的主要参数介绍,其实需要调参的就是两个参数eps和min_samples,这两个值的组合对最终的聚类效果有很大的影响。

7.2、scikit-learn DBSCAN聚类实例

生成一组随机数据,为了体现DBSCAN在非凸数据的聚类优点,我们生成了三簇数据,两组是非凸的。代码如下:

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
X1, y1=datasets.make_circles(n_samples=5000, factor=.6,noise=.05)
X2, y2 = datasets.make_blobs(n_samples=1000, n_features=2, centers=[[1.2,1.2]], cluster_std=[[.1]],random_state=9)

X = np.concatenate((X1, X2))
plt.scatter(X[:, 0], X[:, 1], marker='o')
plt.show()

可以直观看看我们的样本数据分布输出:
算法模型---聚类分析之DBSCAN密度聚类算法
首先我们看看K-Means的聚类效果,代码如下:

from sklearn.cluster import KMeans
y_pred = KMeans(n_clusters=3, random_state=9).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()

K-Means对于非凸数据集的聚类表现不好,从上面代码输出的聚类效果图可以明显看出,输出图如下:
算法模型---聚类分析之DBSCAN密度聚类算法
那么如果使用DBSCAN效果如何呢?我们先不调参,直接用默认参数,看看聚类效果,代码如下:

from sklearn.cluster import DBSCAN
y_pred = DBSCAN().fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()

发现输出让我们很不满意,DBSCAN居然认为所有的数据都是一类!输出效果图如下:
算法模型---聚类分析之DBSCAN密度聚类算法

怎么办?看来我们需要对DBSCAN的两个关键的参数eps和min_samples进行调参!从上图我们可以发现,类别数太少,我们需要增加类别数,那么我们可以减少ϵ-邻域的大小,默认是0.5,我们减到0.1看看效果。代码如下:

y_pred = DBSCAN(eps = 0.1).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()

算法模型---聚类分析之DBSCAN密度聚类算法
可以看到聚类效果有了改进,至少边上的那个簇已经被发现出来了。此时我们需要继续调参增加类别,有两个方向都是可以的,一个是继续减少eps,另一个是增加min_samples。我们现在将min_samples从默认的5增加到10,代码如下:

y_pred = DBSCAN(eps = 0.1, min_samples = 10).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()

输出的效果图如下:
算法模型---聚类分析之DBSCAN密度聚类算法

相关标签: 算法 clustering