本文的代码与数据地址已上传至github:https://github.com/helloWorldchn/MachineLearning
一、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内的邻域,用NEps(p)表示点p的Eps半径内的点的集合。
- 核心对象 (core points):若某个点的Eps邻域内的样本点数(密度)达到算法设定的阈值MinPts,则其为核心点。(即Eps领域内的样本数量不小于MinPts)
- 直接密度可达(directly density-reachable):若某点p在点q的Eps邻域内,且q是核心点,则称对象p从核心对象q是直接密度可达。
- 密度可达density-reachable):若有一个点的序列q0,q1,…,qk,对任意qi到qi+1是直接密度可达的,则称从q0到qk密度可达,这实际上是直接密度可达的“传播”。
- 密度相连(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代码实现
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler, LabelEncoder
from sklearn.metrics import normalized_mutual_info_score # NMI
from sklearn.metrics import rand_score # RI
from sklearn.metrics import accuracy_score # ACC
from sklearn.metrics import f1_score # F-measure
from sklearn.metrics import adjusted_rand_score # ARI
# from sklearn.cluster import DBSCAN
from sklearn.decomposition import PCA
# from sklearn.preprocessing import LabelEncoder
# 数据保存在.csv文件中
# iris = pd.read_csv("dataset/iris.csv") # 鸢尾花数据集 Iris class=3
# wine = pd.read_csv("dataset/wine.csv") # 葡萄酒数据集 Wine class=3
# seeds = pd.read_csv("dataset/seeds.csv") # 小麦种子数据集 seeds class=3
# wdbc = pd.read_csv("dataset/wdbc.csv") # 威斯康星州乳腺癌数据集 Breast Cancer Wisconsin Diagnostic class=2
# glass = pd.read_csv("dataset/glass.csv") # 玻璃辨识数据集 Glass Identification class=6
# aggregation = pd.read_csv("dataset/aggregation.csv") # 人工数据集 Eps=0.18 MinPts=4
# flame = pd.read_csv("dataset/flame.csv") # 人工数据集 Eps=0.28 MinPts=4
# jain = pd.read_csv("dataset/jain.csv") # 人工数据集 Eps = 0.315 MinPts = 4
spiral = pd.read_csv("dataset/spiral.csv") # 人工数据集 Eps=0.45 MinPts=4
df = spiral # 设置要读取的数据集
# print(df)
columns = list(df.columns) # 获取数据集的第一行,第一行通常为特征名,所以先取出
features = columns[:len(columns) - 1] # 数据集的特征名(去除了最后一列,因为最后一列存放的是标签,不是数据)
dataset = df[features] # 预处理之后的数据,去除掉了第一行的数据(因为其为特征名,如果数据第一行不是特征名,可跳过这一步)
attributes = len(df.columns) - 1 # 属性数量(数据集维度)
original_labels = list(df[columns[-1]]) # 原始标签
UNCLASSIFIED = 0
NOISE = -1
# 计算数据点两两之间的距离
def getDistanceMatrix(datas):
N, D = np.shape(datas)
dists = np.zeros([N, N])
for i in range(N):
for j in range(N):
vi = datas[i, :]
vj = datas[j, :]
dists[i, j] = np.sqrt(np.dot((vi - vj), (vi - vj)))
return dists
# 寻找以点cluster_id 为中心,eps 为半径的圆内的所有点的id
def find_points_in_eps(point_id, eps, dists):
index = (dists[point_id] <= eps)
return np.where(index == True)[0].tolist()
# 聚类扩展
# dists : 所有数据两两之间的距离 N x N
# labs : 所有数据的标签 labs N,
# cluster_id : 一个簇的标号
# eps : 密度评估半径
# seeds: 用来进行簇扩展的点
# min_points: 半径内最少的点数
def expand_cluster(dists, labs, cluster_id, seeds, eps, min_points):
i = 0
while i < len(seeds):
# 获取一个临近点
Pn = seeds[i]
# 如果该点被标记为NOISE 则重新标记
if labs[Pn] == NOISE:
labs[Pn] = cluster_id
# 如果该点没有被标记过
elif labs[Pn] == UNCLASSIFIED:
# 进行标记,并计算它的临近点 new_seeds
labs[Pn] = cluster_id
new_seeds = find_points_in_eps(Pn, eps, dists)
# 如果 new_seeds 足够长则把它加入到seed 队列中
if len(new_seeds) >= min_points:
seeds = seeds + new_seeds
i = i + 1
def dbscan(datas, Eps, MinPts):
# 计算 所有点之间的距离
dists = getDistanceMatrix(datas)
# 将所有点的标签初始化为UNCLASSIFIED
n_points = datas.shape[0]
labs = [UNCLASSIFIED] * n_points
cluster_id = 0
# 遍历所有点
for point_id in range(0, n_points):
# 如果当前点已经处理过了
if not (labs[point_id] == UNCLASSIFIED):
continue
# 没有处理过则计算临近点
seeds = find_points_in_eps(point_id, Eps, dists)
# 如果临近点数量过少则标记为 NOISE
if len(seeds) < MinPts:
labs[point_id] = NOISE
else:
# 否则就开启一轮簇的扩张
cluster_id = cluster_id + 1
# 标记当前点
labs[point_id] = cluster_id
expand_cluster(dists, labs, cluster_id, seeds, Eps, MinPts)
return labs, cluster_id
# 计算聚类指标
def clustering_indicators(labels_true, labels_pred):
if type(labels_true[0]) != int:
labels_true = LabelEncoder().fit_transform(df[columns[len(columns) - 1]]) # 如果标签为文本类型,把文本标签转换为数字标签
f_measure = f1_score(labels_true, labels_pred, average='macro') # F值
accuracy = accuracy_score(labels_true, labels_pred) # ACC
normalized_mutual_information = normalized_mutual_info_score(labels_true, labels_pred) # NMI
rand_index = rand_score(labels_true, labels_pred) # RI
ARI = adjusted_rand_score(labels_true, labels_pred)
return f_measure, accuracy, normalized_mutual_information, rand_index, ARI
# 绘图
def draw_cluster(datas, labs, n_cluster):
plt.cla()
colors = np.array(["red", "blue", "green", "orange", "purple", "cyan", "magenta", "beige", "hotpink", "#88c999"])
if attributes > 2:
datas = PCA(n_components=2).fit_transform(datas) # 如果属性数量大于2,降维
# plt.scatter(datas[:, 0], datas[:, 1], s=7., color="black")
# plt.show()
for i, lab in enumerate(labs):
if lab == NOISE:
plt.scatter(datas[i, 0], datas[i, 1], s=7., color=(0, 0, 0))
else:
plt.scatter(datas[i, 0], datas[i, 1], s=7., color=colors[lab - 1])
plt.show()
if __name__ == "__main__":
# 设置DBSCAN参数
Eps = 0.45
MinPts = 4
datas = np.array(dataset).astype(np.float32)
# 数据正则化
datas = StandardScaler().fit_transform(datas)
label, cluster_id = dbscan(datas, Eps=Eps, MinPts=MinPts) # 执行DBSCAN
print(original_labels) # 输出原始标签
print("label of my dbscan")
print(label) # 输出聚类结果
F_measure, ACC, NMI, RI, ARI = clustering_indicators(original_labels, label) # 计算聚类指标
print("F_measure:", F_measure, "ACC:", ACC, "NMI", NMI, "RI", RI, "ARI", ARI)
draw_cluster(datas, label, cluster_id) # 画散点图
# db = DBSCAN(Eps=Eps, MinPts=MinPts).fit(datas) # 调用sk的dbscan
# skl_labels = db.labels_
# print("label of sk-DBSCAN")
# print(skl_labels)
# draw_cluster(db, skl_labels, cluster_id)
2、聚类指标
本文选择了F值(F-measure,FM)、准确率(Accuracy,ACC)、标准互信息(Normalized Mutual Information,NMI)和兰德指数(Rand Index,RI)作为评估指标,其值域为[0,1],取值越大说明聚类结果越符合预期。
F值结合了精度(Precision)与召回率(Recall)两种指标,它的值为精度与召回率的调和平均,其计算公式见公式:
P r e c i s i o n = T P T P + F P Precision=\frac{TP}{TP+FP} Precision=TP+FPTP
R e c a l l = T P T P + F N Recall=\frac{TP}{TP+FN} Recall=TP+FNTP
F − m e a s u r e = 2 R e c a l l × P r e c i s i o n R e c a l l + P r e c i s i o n F-measure=\frac{2Recall \times Precision}{Recall+Precision} F−measure=Recall+Precision2Recall×Precision
ACC是被正确分类的样本数与数据集总样本数的比值,计算公式如下:
A C C = T P + T N T P + T N + F P + F N ACC=\frac{TP+TN}{TP+TN+FP+FN} ACC=TP+TN+FP+FNTP+TN
其中,TP(True Positive)表示将正类预测为正类数的样本个数,TN (True Negative)表示将负类预测为负类数的样本个数,FP(False Positive)表示将负类预测为正类数误报的样本个数,FN(False Negative)表示将正类预测为负类数的样本个数。
NMI用于量化聚类结果和已知类别标签的匹配程度,相比于ACC,NMI的值不会受到族类标签排列的影响。计算公式如下:
N M I = I ( U , V ) H ( U ) H ( V ) NMI=\frac{I\left(U,V\right)}{\sqrt{H\left(U\right)H\left(V\right)}} NMI=H(U)H(V)I(U,V)
其中H(U)代表正确分类的熵,H(V)分别代表通过算法得到的结果的熵。
其具体实现代吗如下:
由于数据集中给定的正确标签可能为文本类型而不是数字标签,所以在计算前先判断数据集的标签是否为数字类型,如果不是,则转化为数字类型
def clustering_indicators(labels_true, labels_pred):
if type(labels_true[0]) != int:
labels_true = LabelEncoder().fit_transform(df[columns[len(columns) - 1]]) # 如果标签为文本类型,把文本标签转换为数字标签
f_measure = f1_score(labels_true, labels_pred, average='macro') # F值
accuracy = accuracy_score(labels_true, labels_pred) # ACC
normalized_mutual_information = normalized_mutual_info_score(labels_true, labels_pred) # NMI
rand_index = rand_score(labels_true, labels_pred) # RI
return f_measure, accuracy, normalized_mutual_information, rand_index
F_measure, ACC, NMI, RI = clustering_indicators(labels_number, labels)
print("F_measure:", F_measure, "ACC:", ACC, "NMI", NMI, "RI", RI)
如果需要计算出聚类分析指标,只要将以上代码插入实现代码中即可。
3、聚类结果散点图
- aggregation数据集
原图:
聚类效果图(Eps=0.18,MinPts=4):
- jain数据集
原图:
聚类效果图(Eps=0.315,MinPts=4):
- Spiral数据集
原图:
聚类效果图(Eps=0.45,MinPts=4):