100-非监督学习之DBSCAN密度聚类

作者: wonphen | 来源:发表于2020-11-10 13:49 被阅读0次

> library(pacman)
> p_load(dplyr, dbscan, ggplot2)

k-means （K均值）和分层聚类衡量行间、及行与中心点的距离。k-Means算法比较适用于簇为球型的，对于非球型的，一般需要基于密度的聚类，比如DBSCAN, OPTICS，利用单位空间的样本量，即密度。基于密度的聚类不偏向球形聚类，可以找到形状各异且复杂的聚类。

DBSCAN(Density-Based Spatial Clustering of Applications with Noise，具有噪声的基于密度的聚类方法)是一种很典型的密度聚类算法，和K-Means，BIRCH这些一般只适用于凸样本集的聚类相比，DBSCAN既可以适用于凸样本集，也可以适用于非凸样本集。通过将紧密相连的样本划为一类，这样就得到了一个聚类类别。通过将所有各组紧密相连的样本划为各个不同的类别，则得到最终的所有聚类类别结果。

r-邻域：给定点半径为r的区域。
核心点：如果一个点的r邻域内最少包含M个点，则该点称为核心点。
直接密度可达：对于核心点P而言，如果另一个点O在P的r邻域内，那么称O为P的直接密度可达点。
密度可达：对于P的直接密度可达点O的r邻域内，如果包含另一个点Q，那么称Q为P的密度可达点。
密度相连：如果Q和N都是核心点P的密度可达点，但是并不在一条直线路径上，那么称两者为密度相连。

1 密度聚类算法思想

1.指定r和M。
2.计算所有的样本点，如果点P的r邻域内有超过M个点，那么创建一个以P为核心点的新簇。
3.反复寻找这些核心点的直接密度可达点（之后可能是密度可达），将其加入到相应的簇，对于核心点发生密度相连的情况加以合并。
4.当没有新的点加入到任何簇中时，算法结束。

2 DBSCAN算法实例

dbscan包函数列表：
dbscan(), 实现DBSCAN算法
optics(), 实现OPTICS算法
hdbscan(), 实现带层次DBSCAN算法
sNNclust(), 实现共享聚类算法
jpclust(), Jarvis-Patrick聚类算法
lof(), 局部异常因子得分算法
extractFOSC()，集群优选框架，可以通过参数化来执行聚类。
frNN(), 找到固定半径最近的邻居
kNN(), 最近邻算法，找到最近的k个邻居
sNN(), 找到最近的共享邻居数量
pointdensity(), 计算每个数据点的局部密度
kNNdist()，计算最近的k个邻居的距离
kNNdistplot()，画图，最近距离
hullplot(), 画图，集群的凸壳

> data(banknote, package = "mclust")
> bn <- as_tibble(banknote)
> str(bn)

## tibble [200 × 7] (S3: tbl_df/tbl/data.frame)
##  $ Status  : Factor w/ 2 levels "counterfeit",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Length  : num [1:200] 215 215 215 215 215 ...
##  $ Left    : num [1:200] 131 130 130 130 130 ...
##  $ Right   : num [1:200] 131 130 130 130 130 ...
##  $ Bottom  : num [1:200] 9 8.1 8.7 7.5 10.4 9 7.9 7.2 8.2 9.2 ...
##  $ Top     : num [1:200] 9.7 9.5 9.6 10.4 7.7 10.1 9.6 10.7 11 10 ...
##  $ Diagonal: num [1:200] 141 142 142 142 142 ...

> DataExplorer::profile_missing(bn)

## # A tibble: 7 x 3
##   feature  num_missing pct_missing
##   <fct>          <int>       <dbl>
## 1 Status             0           0
## 2 Length             0           0
## 3 Left               0           0
## 4 Right              0           0
## 5 Bottom             0           0
## 6 Top                0           0
## 7 Diagonal           0           0

去掉类别列，数据标准化：

> # 因为做聚类，所以去掉类别列
> bn <- bn[, -1] %>% 
+   # 标准化
+   mutate(across(everything(), scale))
> str(bn)

## tibble [200 × 6] (S3: tbl_df/tbl/data.frame)
##  $ Length  : num [1:200, 1] -0.255 -0.786 -0.255 -0.255 0.276 ...
##   ..- attr(*, "scaled:center")= num 215
##   ..- attr(*, "scaled:scale")= num 0.377
##  $ Left    : num [1:200, 1] 2.43 -1.17 -1.17 -1.17 -1.44 ...
##   ..- attr(*, "scaled:center")= num 130
##   ..- attr(*, "scaled:scale")= num 0.361
##  $ Right   : num [1:200, 1] 2.83 -0.635 -0.635 -0.882 -0.635 ...
##   ..- attr(*, "scaled:center")= num 130
##   ..- attr(*, "scaled:scale")= num 0.404
##  $ Bottom  : num [1:200, 1] -0.289 -0.912 -0.497 -1.327 0.68 ...
##   ..- attr(*, "scaled:center")= num 9.42
##   ..- attr(*, "scaled:scale")= num 1.44
##  $ Top     : num [1:200, 1] -1.184 -1.433 -1.308 -0.312 -3.675 ...
##   ..- attr(*, "scaled:center")= num 10.7
##   ..- attr(*, "scaled:scale")= num 0.803
##  $ Diagonal: num [1:200, 1] 0.448 1.056 1.49 1.316 1.143 ...
##   ..- attr(*, "scaled:center")= num 140
##   ..- attr(*, "scaled:scale")= num 1.15

使用kNN()函数，计算数据集中每个值最近的5个点。

> nn <- kNN(bn, k = 5)
> head(nn$dist)

##              1         2         3         4         5
## [1,] 2.3979290 2.6547499 2.7136683 2.7506359 2.7758450
## [2,] 0.8114257 0.9311021 0.9597960 1.0923877 1.1280840
## [3,] 0.6325241 0.8114257 0.8644072 0.9798160 1.1217324
## [4,] 0.7787735 0.8736941 0.8803796 0.8863102 0.9570431
## [5,] 1.9323265 1.9571331 2.0636176 2.1836362 2.2670993
## [6,] 1.6226807 1.6650844 1.6863721 1.8445092 1.8529403

行为每个点的索引，列为最近邻的５个点。
查找与第79号点（79行）距离最近的5个点：

> nn$dist[79, ]

##         1         2         3         4         5 
## 0.8936911 0.9257298 1.0673777 1.1936619 1.2228467

画图：

> # 79号点使用红色，其他黑色
> cols <- ifelse(1:nrow(bn) %in% nn$id[79, ], "red", "black")
> # 近邻点使用蓝色
> cols[79] <- "blue"
> 
> # 避免拥挤，只画3个特征
> plot(bn[, 1:3], pch = 19, col = cols)

近邻点示意图

三个特征两两组合的二维平面图中，可以看出红色点确实围绕在蓝色点周围，但同时，因为是二维平面，所以有些点被遮盖了。
选取两个列，画出最近邻前5连接路径：

> plot(nn, bn)

最近邻连接路径

通过连接路径，可以看到最近邻的分组过程，能够连接在一起的就组成了一个聚类，没有连接在一起的就聚为了不同的类。

DBSACN算法函数语法：
eps：搜索半径，设置得非常小，则意味着没有点是核心样本，可能会导致所有点被标记为噪声；设置得非常大，可能会导致所有点形成单个簇。
minPts：成为聚类的最少的行数，要成为核心对象所需要的 ϵ-邻域的样本数阈值，默认为5。
weights，数据点的权重，仅用于加权聚类。
borderPoints，边界点是否为噪声，默认为TRUE。

> dbscan(x, eps = 0.42, minPts = 5)

寻找最优的参数

方法一

eps，可以使用绘制k-距离曲线(k-distance graph)方法得到，在k-距离曲线图明显拐点位置为较好的参数。若参数设置过小，大部分数据不能聚类；若参数设置过大，多个簇和大部分对象会归并到同一个簇中。
minPts，通常让minPts≥dim+1，其中dim表示数据集聚类数据的维度。若该值选取过小，则稀疏簇中结果由于密度小于minPts，从而被认为是边界点；若该值过大，则密度较大的两个邻近簇可能被合并为同一簇。
本例中数据集为6维，所以选择k=6+1=7。

> kNNdistplot(bn, k = 7)
> abline(h = 2, col = "red", lty = 2)

K-距离曲线

kNNdistplot()会计算点矩阵中的k=7的最近邻的距离，然后按距离从小到大排序后，以图形进行展示。x轴为距离的序号，y轴为距离的值。图中黑色的线，从左到右y值越来越大。
通过绘制k-距离曲线，寻找eps，即明显拐点位置为对应较好的参数。本例中eps为2。
最后使用参数聚类：

> bn.dbscan <- dbscan(bn, eps = 2, minPts = 7)
> bn.dbscan

## DBSCAN clustering for 200 objects.
## Parameters: eps = 2, minPts = 7
## The clustering contains 1 cluster(s) and 4 noise points.
## 
##   0   1 
##   4 196 
## 
## Available fields: cluster, eps, minPts

从结果可知，整个数据集聚为了一个类，其中4个点为噪声点。但实际数据中只有两个类，所以说明eps或者minPts参数设置得过大。

方法二

使用参数网格寻找最优的参数。

> # 参数网格
> # eps从1.0到3.0，每次0.2增加
> # minPts从5到10
> # 一共11*6=66组
> grid.dbscan <- expand.grid(eps = seq(1.0, 3.0, 0.2),
+                            MinPts = 5:10)
> 
> dim(grid.dbscan)

## [1] 66  2

> head(grid.dbscan)

##   eps MinPts
## 1 1.0      5
## 2 1.2      5
## 3 1.4      5
## 4 1.6      5
## 5 1.8      5
## 6 2.0      5

> # 训练模型
> # 将grid.dbscan中的每一行传给dbscan()函数，bn为dbscan的数据框
> bn.dbs <- purrr::pmap(.l = grid.dbscan, .f = dbscan, bn)
> 
> # 返回所有模型结果保存在list中
> class(bn.dbs)

## [1] "list"

> # 查看其中一个模型的输出
> bn.dbs[[1]]

## DBSCAN clustering for 200 objects.
## Parameters: eps = 1, minPts = 5
## The clustering contains 3 cluster(s) and 81 noise points.
## 
##  0  1  2  3 
## 81 55  4 60 
## 
## Available fields: cluster, eps, minPts

> # 查看其中一个聚类的结果
> bn.dbs[[11]]$cluster

##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [54] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [107] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [160] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

> bn.clust <- bn.dbs %>% 
+   # 提取所有的聚类结果
+   purrr::map_dfc(~.$cluster) %>% 
+   # 设置列名
+   setNames(paste0("cluster", 1:66))
> bn.clust[1:6, 1:8]

## # A tibble: 6 x 8
##   cluster1 cluster2 cluster3 cluster4 cluster5 cluster6 cluster7 cluster8
##      <int>    <int>    <int>    <int>    <int>    <int>    <int>    <int>
## 1        0        0        0        0        0        0        0        1
## 2        1        1        1        1        1        1        1        1
## 3        1        1        1        1        1        1        1        1
## 4        1        1        1        1        1        1        1        1
## 5        0        0        0        0        0        1        1        1
## 6        0        0        0        0        1        1        1        1

超参数结果的可视化比较：（因为模型太多，一共66个，所以图形看不清，这里将不展示）

> bn.clust %>% 
+   # 合并原数据框中的两列
+   bind_cols(right = bn$Right,
+             diagonal = bn$Diagonal) %>%
+   # 宽表变长表，只变换聚类结果列
+   tidyr::pivot_longer(names_to = "permutation", 
+                       values_to = "cluster", cols = c(1:66)) %>% 
+   # 画图
+   ggplot(aes(right, diagonal, col = as.factor(cluster))) +
+   geom_point() +
+   facet_wrap(~ permutation) +
+   theme_bw() +
+   theme(legend.position = "none")

超参数结果的指标比较：

> p_load(clusterSim, purrr)
> 
> # 计算分类效果函数
> msr_cluster <- function(data, clusters, dist_mat) {
+   # 计算戴维森堡丁指数，DB越小意味着类内距离越小，同时类间距离越大
+   list(db = index.DB(data, clusters)$DB,
+        # Calinski-Harabasz伪F统计量，CH越大代表着类自身越紧密，类与类之间越分散
+        G1 = index.G1(data, clusters),
+        # 聚类数量
+        clusters = length(unique(clusters)))
+ }
> 
> # 十折交叉验证
> bn.boot <- map(1:10, ~ {
+   # 有放回抽样10次
+   sample_n(bn, size = nrow(bn), replace = T)
+ })
> 
> class(bn.boot)

## [1] "list"

> length(bn.boot)

## [1] 10

> # 对每一次抽样进行聚类，并计算结果
> metric.bn <- map_df(bn.boot, function(boot) {
+   clust.out = pmap(grid.dbscan, dbscan, boot)
+   
+   map_df(clust.out, function(permutation) {
+     clust = permutation$cluster %>%
+       as_tibble() %>% 
+       bind_cols(boot) %>% 
+       # 去掉噪声点
+       filter(value != 0)
+     
+   d = dist(dplyr::select(clust, -value))
+     
+   msr_cluster(data = clust, 
+               clusters = clust$value,
+               dist_mat = d)
+   })
+ })
> 
> str(metric.bn)

## tibble [660 × 3] (S3: tbl_df/tbl/data.frame)
##  $ db      : num [1:660] 0.628 0.793 0.926 NaN NaN ...
##  $ G1      : num [1:660] 89.3 87.4 95.9 NaN NaN ...
##  $ clusters: int [1:660] 5 4 3 1 1 1 1 1 1 1 ...

> summary(metric.bn)

##        db               G1           clusters    
##  Min.   :0.4627   Min.   : 10.7   Min.   :1.000  
##  1st Qu.:0.7133   1st Qu.:112.4   1st Qu.:1.000  
##  Median :0.8091   Median :138.5   Median :1.000  
##  Mean   :0.8180   Mean   :142.0   Mean   :1.645  
##  3rd Qu.:0.9527   3rd Qu.:165.3   3rd Qu.:2.000  
##  Max.   :1.0715   Max.   :270.0   Max.   :8.000  
##  NA's   :463      NA's   :463

结果中有463个缺失值，最少的聚为了一个类，最多的聚为了8个类。

> metric.sum <- metric.bn %>% 
+   # 增加三列
+   mutate(bootstrap = factor(rep(1:10, each = 66)),
+          eps = factor(rep(grid.dbscan$eps, times = 10)),
+          MinPts = factor(rep(grid.dbscan$MinPts, times = 10))) %>% 
+   # 将Na转换为Inf（无穷大）
+   mutate(across(where(is.numeric), ~ na_if(., Inf))) %>% 
+   # 去掉有缺失值的行
+   tidyr::drop_na()
> str(metric.sum)

## tibble [197 × 6] (S3: tbl_df/tbl/data.frame)
##  $ db       : num [1:197] 0.628 0.793 0.926 0.596 0.78 ...
##  $ G1       : num [1:197] 89.3 87.4 95.9 114 114.2 ...
##  $ clusters : int [1:197] 5 4 3 4 3 2 5 2 2 5 ...
##  $ bootstrap: Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ eps      : Factor w/ 11 levels "1","1.2","1.4",..: 1 2 3 1 2 3 1 2 3 1 ...
##  $ MinPts   : Factor w/ 6 levels "5","6","7","8",..: 1 1 1 2 2 2 3 3 3 4 ...

此时只需要找到db最小同时G1最大（即db/G1最小）时的参数即可。
画图查看结果：

> metric.sum %>% 
+   mutate(value = round(db / G1, 3)) %>% 
+   dplyr::select(eps, MinPts, clusters, value) %>% 
+   # 将clusters转换为因子型
+   transform(clusters = as.factor(.$clusters)) %>% 
+   # 将1、2、3列宽变长
+   tidyr::pivot_longer(names_to = "metrics", values_to = "dbs", cols = 1:3) %>% 
+   # 找出value值最小的行
+   filter(value == min(value)) %>% 
+   ggplot(aes(dbs, value)) +
+   geom_point(size = 2, col = "blue") +
+   facet_wrap(~ metrics) +
+   theme_bw() +
+   labs(x = "")

最优参数可视化

最终选择的聚类数量为2，此时eps=1，MinPts=9或者10。

3 密度聚类的优缺点

优点：
（1）聚类速度快且能够有效处理噪声点和发现任意形状的空间聚类；
（2）与K-MEANS比较起来，不需要输入要划分的聚类个数；
（3）聚类簇的形状没有偏倚；
（4）可以在需要时输入过滤噪声的参数。

缺点：
（1）当数据量增大时，要求较大的内存，I/O消耗也很大；
（2）当空间聚类的密度不均匀、聚类间距差相差很大时，聚类质量较差，因为这种情况下参数MinPts和eps选取困难
（3）聚类效果依赖于距离公式选取，实际应用中常用欧式距离，对于高维数据，存在“维数灾难”。

100-非监督学习之DBSCAN密度聚类

1 密度聚类算法思想

2 DBSCAN算法实例

3 密度聚类的优缺点

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

R炒面