基于密度的聚类算法

普通的聚类算法只能发现一些特定形状的簇群，对于一些特殊形状簇群分类效果就很差。

而基于密度的分类就可以发现任意形状的簇。

先来一张效果图。

一个颜色表示一个簇群

可以看出来算法发现很多不同形状的簇群。

以下就是算法的设计和实现。

算法聚类了各种不同的簇群，所以第一步是要实现一个，可以找出一个簇群的算法。

算法的设计思想就是

找一个点，然后以这个为中心，给定距离为半径画圆，然后看这个圆中有多少个点，如果点的数量大于给定的阈值，那么这个点就是核心点，如果不是不是那么就是边界点或者噪声点。

之后用这个方法去判断每一个刚刚圈中点的，看他是不是边界点，或者噪声点。

边界圈中的其他点就不需要再次判断，噪声点也是，但边界点会直接加入到导致他被圈中的核心点的簇群，而噪声点不会被加入到任何簇群中。

算法设计:

写一个函数，可以统计被圈中的数据的索引标签。

getPointCircle.m

function re_labels = getPointCircle(data,center,radius)

%Author:Q
�te:2013/9/6
�scribe：得到以center为中心,radius为半径的data中的数据点

%参数说明:
%   输入参数:
%       data            输入数据
%       center          中心数据
%       radius          半径
%
%   输出参数:
%       re_labels       标签，圈中数据的标签

distance    = euclideanDistance(data,center);
index       = distance < radius;
re_labels   = index; 

这个函数还使用到了另一个函数，这个函数用于求距离，求一组数据离中心点的欧拉距离（如果以后需要变动，直接变动这个恶函数即可，例如可以变为马氏距离）。

euclideanDistance.m

function re_distance = euclideanDistance(data,center)

%Author:Q
�te:2013/9/1
�scribe:计算欧几里得距离的函数

%参数说明:
%   输入参数:
%       data            需要计算的数据，支持多个数据，行为数据，列为维度
%       center          计算离这个点的欧式距离，默认为原点
%   数据参数:
%       re_distance     计算所得的欧式距离，支持多个数据

[row column] = size(data);

if nargin == 1
    center = zeros(1,column);
end

tmpdis = data;

for i = 1:column
    tmpdis(:,i) = data(:,i) - center(i);
end

tmpdis      = tmpdis .* tmpdis;
re_distance = zeros(row,1);

for i = 1:row
    re_distance(i) = sqrt(sum(tmpdis(i,:)));
end 

最后就是寻找簇群的函数。

getNewCluster.m

function re_lables = getNewCluster(data,center_lab,radius,MinPts)

%Author:Q
�te:2013/9/6
�scribe:找出中心为center范围为radius的簇群
%中心点有资格拉点，边界点没有资格

%参数说明:
%   输入参数:
%       data            待聚类数据，默认为没有进行任何处理的，单纯的n为数据
%       center_lab      中心点的索引，为什么写成这个样子？为什么不直接输入点？这样方便一点
%       radius          Eps领域的半径
%       MinPts          点密度，也就是点的数量
%
%   输出参数:
%       re_labels       分出的新簇群的索引号


[row column] = size(data);

labels = zeros(row,1);          %标记那些数据被标记过
stackmark = zeros(row,1) == 1;

checkpt = [];                   %保存所有要进行查找的点

checkpt = [checkpt ; data(center_lab,:)];
stackmark(center_lab) = true;
checkpos = 1;

judge = 0;

while judge == 0
    [tmprow tmpcolumn] = size(checkpt);
    if tmprow < checkpos
        judge = 1;
    else
        center_pt   = checkpt(checkpos,:);
        Eps_labels  = getPointCircle(data,center_pt,radius);
        
        if sum(Eps_labels) >= MinPts
                                            %这个是核心点，可以拉别的点进来
            mark_label              = cmpVector(data,center_pt);
            labels(mark_label)      = 1;         %核心点赋1，表示处理完毕
            Eps_labels(stackmark)   = false;
            Eps_labels(mark_label)  = false;
            checkpt                 = [checkpt ; data(Eps_labels,:)];
            stackmark(Eps_labels)   = true;
        else                                %这个是边界点，不能拉别的点进来
            mark_label = cmpVector(data,center_pt);
            labels(mark_label) = 1;
        end
        
        checkpos = checkpos + 1;
    end
end

re_lables = labels == 1;

最后一个测试代码查看函数是否正确工作。

clear all
close all
clc

datanum = 80;
data = abs(rand(datanum,2))*100;

%function re_lables = getNewCluster(data,center_lab,radius,MinPts)

center_lab = zeros(datanum,1);

center_lab(25) = 1;

labels = center_lab == 1;

re_labels = getNewCluster(data,labels,10,2);

tmpdata = data(re_labels,:);

figure
plot(data(:,1),data(:,2),'*b');
hold on
plot(tmpdata(:,1),tmpdata(:,2),'*','Color',[0.5 0.5 0.5]); 

测试结果为：

可见函数分类正确找出了同密度的簇群，而且形状不规则。

之后就只需要重复调用该函数对没有分类的数据分类，就可以得到最后的结果

函数代码:

dbScan.m

function [re_labels num] = dbScan(data,MinPts,Eps)

%Author:Q
�te:2013/9/6
�scribe:密度聚类算法

%参数说明:
%   输入参数:
%       data        待聚类数据
%       MinPts      领域的中心点
%       Eps         邻域半径
%
%   输出参数:
%       re_labels   发分类之后的标签
%       num         有多少个簇群

mark = 1;               %分类标记

[row column] = size(data);

labels = zeros(row,1);          %数据的标签

index_deal = labels == 0;       %未分类数据的索引号

for i = 1:row
    if (labels(i) == 0)         %表示数据未被分类
        center_lab = cmpVector(data,data(i,:));
        labels_ind = getNewCluster(data,center_lab,Eps,MinPts);
        
        
        if sum(labels_ind) == 1
            labels(labels_ind) = -1;            %噪声点为-1
            index_deal(labels_ind) = false;
        else                                    %一个普通的簇群点
            labels(labels_ind) = mark;
            mark = mark + 1;
            index_deal(labels_ind) = false;
        end
    end
end

re_labels = labels;
num = mark; 

这个函数使用了另一个cmpVector

因为整个过程都是同一份数据，所有对数据的操作都是索引号。

所以用到一个cmpVector去寻找数据的索引

cmpVector.m

function re_labels = cmpVector(data,vector)

%Author:Q
�te:2013/9/6
�scribe:找出Data中行为vector的矩阵的索引
%目的是为了解决不用直接data == vector这个问题，要找出符合条件的索引

%参数说明:
%   输入参数:
%       data        待比较的矩阵
%       vector      待比较的向量
%
%   输出参数:
%       re_labels   返回找出之后的索引

[row column] = size(data);

labels = [];

for i = 1:row
    tmplabe = 0;
    
    if sum(data(i,:) == vector) ~= 0
        tmplabe = 1;
    end
    
    labels = [labels ; tmplabe];
end

re_labels = labels == 1;

最后来一个测试demo

clear all
close all
clc

datanum = 200;
data = abs(rand(datanum,2))*100;

[re_labels labels_num]= dbScan(data,2,8)


figure
index = re_labels == -1;
plot(data(index,1),data(index,2),'.k');

for i = 1:labels_num
    index = re_labels == i;
    hold on
    plot(data(index,1),data(index,2),'*','Color',[abs(rand(1,1)) abs(rand(1,1)) abs(rand(1,1))]);
end

运行查看效果

分类效果还是很不错的。

刚刚设置的是只要圈中两个点就可以算为核心点

改变参数，设置为3个点为核心点(目地是测试函数是否正常工作)。

可以发现两个点的簇群消失，这些点全部被标记为噪声点。