1、先解决一个错误

源代码若以DBSCAN.M为文件名，使用时会出现以下错误：

>> x=[randn(30,2)*.4;randn(40,2)*.5+ones(40,1)*[4 4]];
>> [class,type]=dbscan(x,5,[])
Cannot find an exact (case-sensitive) match for 'dbscan.m'

The closest match is F:\project\clustering\DBSCAN.M

To change the file extension, cd to the file's folder, type:
   movefile DBSCAN.M DBSCAN.m_bad; movefile DBSCAN.m_bad DBSCAN.m
and then cd back.

将DBSCAN.M更名为dbscan.m即可。

2、示例

将dbscan.m放入当前工作区间，

>> x=[randn(3,2)*.4;randn(4,2)*.5+ones(4,1)*[4 4]]
x =
    0.0092   -0.1143
   -0.1048   -0.3325
   -0.7001   -0.3917
    3.4218    4.2600
    3.7332    3.9900
    2.9987    3.9826
    4.4821    3.6009

>> [class,type]=dbscan(x,3,2)
class =
    -1    -1    -1     1     1     1     1
type =
    -1    -1    -1     1     1     1     1

>> [class,type]=dbscan(x,2,2)
class =
     1     1     1     2     2     2     2
type =
     1     1     1     1     1     1     1

>> [class,type]=dbscan(x,2,[])
class =
     1     1     1     2     2     2     2
type =
     1     1     1     1     1     1     1

数据集x含有7个2维的数据点，函数dbscan有三个参数，第一个参数是数据集；第二个参数是核心对象应该满足的其邻域内最少的数据点的个数；第三个参数用来指定邻域半径，如果值为[]，则函数自己判断半径。返回值中class是来标识每个数据点的所属簇的编号，编号从1开始，0表示尚未被分类，-1表示该数据点是噪声，type用来表示数据点是否分类了，0表示未被0表示尚未被分类，-1表示噪声，1表示已经分类。由x本身的特点可以看出，dbscan(x,2,2)和dbscan(x,2,[])都达到了效果；dbscan(x,3,2)因为第二个参数选的不好，效果不好。

3、分析dbscan.m

下面是dbscan.m去掉大部分注释后的内容：

% -------------------------------------------------------------------------
% Function: [class,type]=dbscan(x,k,Eps)
% -------------------------------------------------------------------------
% Input: 
% x - data set (m,n); m-objects, n-variables
% k - number of objects in a neighborhood of an object 
% Eps - neighborhood radius, if not known avoid this parameter or put []
% -------------------------------------------------------------------------
% Output: 
% class - vector specifying assignment of the i-th object to certain 
% cluster (m,1)
% type - vector specifying type of the i-th object 
% -------------------------------------------------------------------------
% Example of use:
% x=[randn(30,2)*.4;randn(40,2)*.5+ones(40,1)*[4 4]];
% [class,type]=dbscan(x,5,[])
% clusteringfigs('Dbscan',x,[1 2],class,type)
% -------------------------------------------------------------------------

function [class,type]=dbscan(x,k,Eps)

[m,n]=size(x);

if nargin<3 | isempty(Eps)
   [Eps]=epsilon(x,k);   %分析可能的邻域半径Eps，关于函数epsilon，见下面的论述
end

x=[[1:m]' x];  %为每个数据点添加编号，其实是合并两个矩阵
[m,n]=size(x);
type=zeros(1,m);
no=1;
touched=zeros(m,1); %标记某数据点是否已被访问

for i=1:m
    if touched(i)==0;
       ob=x(i,:);
       D=dist(ob(2:n),x(:,2:n));  %使用下面的dist函数生成针对当前点的距离矩阵
       ind=find(D<=Eps);   %见下面的论述

       if length(ind)>1 & length(ind)<k+1       
          type(i)=0;
          class(i)=0;
       end
       if length(ind)==1
          type(i)=-1;
          class(i)=-1;  
          touched(i)=1; 
       end

       if length(ind)>=k+1; 
          type(i)=1;
          class(ind)=ones(length(ind),1)*max(no);

          while ~isempty(ind)
                ob=x(ind(1),:);
                touched(ind(1))=1;
                ind(1)=[];
                D=dist(ob(2:n),x(:,2:n));
                i1=find(D<=Eps);

                if length(i1)>1
                   class(i1)=no;
                   if length(i1)>=k+1;
                      type(ob(1))=1;
                   else
                      type(ob(1))=0;
                   end

                   for i=1:length(i1)
                       if touched(i1(i))==0
                          touched(i1(i))=1;
                          ind=[ind i1(i)];   
                          class(i1(i))=no;
                       end                    
                   end
                end
          end
          no=no+1; 
       end
   end
end

i1=find(class==0);
class(i1)=-1;
type(i1)=-1;

%...........................................
function [Eps]=epsilon(x,k)

% Function: [Eps]=epsilon(x,k)
% Input: 
% x - data matrix (m,n); m-objects, n-variables
% k - number of objects in a neighborhood of an object

[m,n]=size(x);

Eps=((prod(max(x)-min(x))*k*gamma(.5*n+1))/(m*sqrt(pi.^n))).^(1/n);

%............................................
function [D]=dist(i,x)

% function: [D]=dist(i,x)						    
% Input: 
% i - an object (1,n)
% x - data matrix (m,n); m-objects, n-variables	    
%                                                                 
% Output: 
% D - Euclidean distance (m,1)

[m,n]=size(x);
D=sqrt(sum((((ones(m,1)*i)-x).^2)'));

if n==1
   D=abs((ones(m,1)*i-x))';
end

3.1、epsilon函数如何分析出可能的邻域半径r（或者是Eps）

epsilon函数中有调用了4个函数，max和min用来获取矩阵每一列的最大、最小值，prod将矩阵的每一列进行列内相乘。示例：

>> x=[1 3 ; 2 4]
x =
     1     3
     2     4

>> max(x)
ans =
     2     4

>> min(x)
ans =
     1     3

>> prod(x)
ans =
     2    12

gamma函数有两个重要的特性，若n为正整数，那么gamma(n) = (n-1)!，另外一个特征就是它与pi的关系：

gamma(1/2) = pi.^(1/2)
gamma(3/2) = 1/2*(pi.^(1/2))
gamma(5/2) = 3/4*(pi.^(1/2))
...

详细请参考资料[2]。

为了更好的理解，对下面的代码（或者说是等式）进行变换：

Eps=((prod(max(x)-min(x))*k*gamma(.5*n+1))/(m*sqrt(pi.^n))).^(1/n);

把Eps换成r方便些：

r=((prod(max(x)-min(x))*k*gamma(.5*n+1))/(m*sqrt(pi.^n))).^(1/n);

最终变换：

r.^n * sqrt(pi.^n)  / gamma(.5*n+1) * (m / k)=  prod(max(x)-min(x))          %（公式1）

公式1中r代表邻域半径，n是数据的维度，m是数据点的个数，k是核心对象的邻域内至少应具有的数据点数。

在二维空间中，n=2，公式左边r.^n * sqrt(pi.^n) / gamma(.5*n+1)代表半径为r的圆的面积；假定所有核心对象的邻域互相没有交集，m/k则是各个核心对象的个数，也是最终的不想交的由邻域构成的圆的个数。公式右边则代表着这样一个矩形的面积：矩形的各个边平行于横坐标轴或者纵坐标轴，切这个矩形是能涵盖所有数据点的最小矩形。由此可见公式1的合理性。同理，n的值更大时，也有相同的结论。所以，epsilon函数求得的邻域半径是合理的。

3.2、距离矩阵与find函数

dbscan.m中查找数据点邻域中数据点的方法是使用距离矩阵，find函数用来查找符合某个条件的索引，例如：

>> a = [4,2,5,3,9,7]
a =
     4     2     5     3     9     7

>> find(a>7)
ans =
     5

>> find(a>3)
ans =
     1     3     5     6

>> a = [4,2,5 ; 3,9,7]
a =
     4     2     5
     3     9     7

>> find(a>3)
ans =
     1
     4
     5
     6

4、可视化聚类结果

要可视化聚类结果很简单，笔者实现了一个可视化2D数据聚类结果的函数：

% -------------
% 输入：
% data：数据集
% class：每个数据点的分类，-1表示噪点，分类从1开始
% -------------

function color_cluster(data, class)
color_theme = {'ro','g+','bx','ks','r+','go','k*','r+','gv','bs','k+', ...
    'rv','gx','bo','kx','rs','gs','b+','ko',...
    'rx','g+','bv','kv'};
non_class_color = 'b*'; %噪点
[m, n] = size(data);
subplot(1,2,1);
plot(data(:,1),data(:,2),'ko');
title('raw data');
subplot(1,2,2);

for x = 1:m
    if class(x) == -1
        plot(data(x,1),data(x,2) ,non_class_color)
    else
        plot(data(x,1),data(x,2) ,cell2mat(color_theme(class(x))))
    end
    hold on;
end
title('result');

我们测试一下：

x=[randn(30,2)*.4;randn(40,2)*.5+ones(40,1)*[4 4]];
[class,type]=dbscan(x,5,[]);
color_cluster(x,class);

结果如下：

5、资料

[1] DBSCAN的matlab实现 https://code.google.com/p/dmfa07/downloads/detail?name=DBSCAN.M

[2] Γ函数 http://zh.wikipedia.org/wiki/Γ函数