在正课第10讲中,我们介绍了常用的聚类算法:Kmeans聚类、系统聚类和Dbscan聚类。
其中,
Dbscan聚类算法最为特殊,它是
一
种基于密度的聚类方法,聚类前不需要预先指定聚类的
个数,下图来自我们的课件:
在课程的最后,我们给出了MATLAB实现Dbscan聚类的代码,该代码下载于MATLAB官网:
https://ww2.mathworks.cn/matlabcentral/fileexchange/52905-dbscan-clustering-algorithm
该代码中借助了pdist2函数,该函数可用来计算两组观测值之间的成对距离。使用该函数计算时,可能会有一个隐藏的问题:当输入矩阵X的行数特别大时,MATLAB可能没有足够的内存计算并报错。
例如,当X的行数为20万时,D矩阵的大小为20万乘以20万,一般的计算机是无法在内存中保存这么大的数据的。
因此,上面链接中给的代码只能用于数据量不大的情况,如果数据量较多则无法借助一般的家用电脑完成计算。
那我们应该如何解决这个问题呢?
从MATLAB2019版本开始,MATLAB强大的统计与机器学习工具箱中内置了dbscan函数,我们不需要下载第三方的代码即可方便的使用dbscan算法进行聚类,其调用的方法可在MATLAB官网找到:
https://ww2.mathworks.cn/help/stats/density-based-spatial-clustering.html
该内置函数的使用方法和课堂上讲解的DBSCAN函数非常类似,大家可以参考正课第10讲。
该函数的核心代码是闭源的,MATLAB官网也没有介绍内部是如何实现的,但我个人推测其内部应该用到了特殊的数据结构来巧妙的节省内存消耗,例如可能用到了并查集来合并多个类别。
经过测试,该函数在面对几十万条数据时也能计算出聚类的结果,因此大家在面对大型数据的dbscan聚类问题时可以选用内置的这个函数,前提是你的MATLAB版本要高于2019且安装好了
统
计与
机器学习工具箱( Statistics and Machine LearningToolbox)
。
(1)在1至500的正方形网格点上随机生成3000个点,使用x和y分别表示这3000个点的横纵坐标。
rng(100) % 随机数种子
idx = randperm(500^2,3000);
[x,y]=ind2sub([500,500],idx);
x = x'; y = y';
这就是随机生成的3000个点的横纵坐标,例如第一个点的坐标为(352,272),第二个点坐标为(93,140),依次类推。
figure('Name','散点分布图','NumberTitle','off');
scatter(x,y,0.5,'k') % 0.5表示散点大小
axis([0,500,0,500]) % 设置x轴y轴范围
(3)对这三千个点进行聚类,规则如下:将距离小于10的点聚成一类,例如A和B距离为8,则将A和B聚类一类,且该聚类具有传递性,假设C和B的距离为5.21,则C也和A、B归于同一类。
epsilon = 10; % 半径
minpts = 1; % 只要包含一个点就可以认为是一类,即不含噪音点
idx = dbscan([x,y],epsilon,minpts);
length(unique(idx)) % 聚类的类别数
从结果可以看出,最终一共分为了211类。(因为数据是随机生成的,你的结果可能和我的有所差异,这是正常现象)
idx中保存了聚类的结果,它是一个长度为3000的向量,数值相同的位置表示该样本处于相同的一类中。
例如我们可以找到第一类中包含的样本,它们位于原样本的第1行、第392行、第425行等共计26个样本点。
另外,我们还可以求出每一个类别中包含的样本点的个数,这里需要借助MATLAB中分组求和的函数:
结果显示,第一个类别中包含26个样本点,第二个类别中包含334个样本点,依次类推。
当然,你也可以按照样本点的个数来对这些类别进行排序:
结果表明,第6类中样本点最多,有4
91个样本点;第二多的是第2类,有334个样本点;依次类推。
figure('Name','DBSCAN聚类结果','NumberTitle','off');
gscatter(x,y,idx,[],[],1,'doleg','off')
xlabel('x坐标'); ylabel('y坐标')
这里我们用到了gcatter函数,它能根据分类结果来绘制散点图,但要注意的是,由于我们的类别太多,MATLAB内置的颜色不足,因此可能有某些不同的类别使用的颜色也相同,因此下面我们可以考虑:
每次随机的取出k个类别画图,将其他的类别全部绘制为黑色:
k = 5;
p = randperm(length((unique(idx))),k);
[id1,pp] = ismember(idx,p); id0 = ~id1;
x1 = x(id1); y1 = y(id1); % 作者:数学建模清风
x0 = x(id0); y0 = y(id0);
figure('Name',"随机显示"+k+"个类别的效果图",'NumberTitle','off');
gscatter(x1,y1,pp(pp~=0),[],'.',8,'doleg','off')
hold on
scatter(x0,y0,1,'.k')
hold off
xlabel('x坐标'); ylabel('y坐标')
注意,这是随机取的k个类别(上图k等于5),因此每次运行的结果可能也不同,下图是k取8时随机生成的一个图形,大家也可以自己测试完成。