理工学—算法模型—聚类分析之DBSCAN密度聚类算法

大家好，欢迎来到IT知识分享网。

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

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

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

样本样例

/* 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() 

• eps： DBSCAN算法参数，即我们的$ϵ$-邻域的距离阈值，和样本距离超过ϵ的样本点不在$ϵ$-邻域内。默认值是0.5。一般需要通过在多组值里面选择一个合适的阈值。eps过大，则更多的点会落在核心对象的$ϵ$-邻域，此时我们的类别数可能会减少， 本来不应该是一类的样本也会被划为一类。反之则类别数可能会增大，本来是一类的样本却被划分开。
• min_samples： DBSCAN算法参数，即样本点要成为核心对象所需要的$ϵ$-邻域的样本数阈值。默认值是5. 一般需要通过在多组值里面选择一个合适的阈值。通常和eps一起调参。在eps一定的情况下，min_samples过大，则核心对象会过少，此时簇内部分本来是一类的样本可能会被标为噪音点，类别数也会变多。反之min_samples过小的话，则会产生大量的核心对象，可能会导致类别数过少。
• metric：最近邻距离度量参数。可以使用的距离度量较多，一般来说DBSCAN使用默认的欧式距离（即p=2的闵可夫斯基距离）就可以满足我们的需求。可以使用的距离度量参数有：
  • a) 欧式距离 “euclidean”: $\sqrt{{\sum }_{i=1}^n(x_i-y_i{)}^2}$
  • b) 曼哈顿距离 “manhattan”：$ \sum_{i=1}^n |x_i−y_i|$
  • c) 切比雪夫距离“chebyshev”😒 max|x_i−y_i|(i=1,2,…n)$
  • d) 闵可夫斯基距离 “minkowski”: $\sqrt[p]{{\sum }_{i=1}^n|x_i-y_i{|}^p}$，p=1为曼哈顿距离， p=2为欧式距离。
  • e) 带权重闵可夫斯基距离 “wminkowski”: $\sqrt[p]{{\sum }_{i=1}^{n}w\ast |x_i-y_i{|}^p}$，其中$w$为特征权重
  • f) 标准化欧式距离 “seuclidean”: 即对于各特征维度做了归一化以后的欧式距离。此时各样本特征维度的均值为0，方差为1.
  • g) 马氏距离“mahalanobis”：$(x-y{)}^T{S}^{-1}(x-y)$，其中，$S^{-1}$为样本协方差矩阵的逆矩阵。当样本分布独立时， 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为欧式距离。如果使用默认的欧式距离不需要管这个参数。

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() 

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() 

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

怎么办？看来我们需要对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() 

可以看到聚类效果有了改进，至少边上的那个簇已经被发现出来了。此时我们需要继续调参增加类别，有两个方向都是可以的，一个是继续减少eps，另一个是增加min_samples。我们现在将min_samples从默认的5增加到10，代码如下： “`python y_pred = DBSCAN(eps = 0.1, min_samples = 10).fit_predict(X) plt.scatter(X[:, 0], X[:, 1], c=y_pred) plt.show() “` 输出的效果图如下：