DBSCAN基于密度的聚类算法

**DBSCAN算法和实现

——DBSCAN(Density-Based Spatial Clustering of Applications with Noise)是一个比较有代表性的基于密度的聚类算法，它是一种适应性极强的聚类算法，不同于K-means聚类，它无需确定簇的数量k，可以自己基于密度来实现。首先介绍一下几个基本的概念。

基础概念

1.密度：密度指的是在某距离内含有对象的最小数目。
2.核心对象：如果一个对象的eps邻域内至少包含MinPt个对象，则称该对象是核心对象。
3.直接密度可达：给定一个对象集合D，如果p在q的eps邻域内，并且q是一个核心对象，则p是从q直接密度可达的。（如下图所示），eps可以想象为一个超球体的半径，Minp就是这个超球体内的点，这样比较容易理解。

4.密度可达：对于q和p之间，有一个样本序列p1，p2，p3，……，pn， 其中p1=q， pn=p，若Pi+1是由Pi直接密度可达的（关于eps和MinPts），则对象p是由q关于eps和MinPts密度可达的。（如下图所示）

算法介绍

在DBSCAN算法中，一个簇就是从一个点出发，所有密度可达的点的集合，如果某一个点不属于任何一个簇，也就是说这个点不能由任何点密度可达，那么这个点被称为异常点。
算法流程：
1、输入：样本D={x1, x2, x3,……., xm}，邻域参数（eps，MinPts），距离的度量方式（本文中我们用欧几里得距离来介绍，事实上大部分时候也用的是欧几里得距离方法）
2、输出：簇划分、异常点索引

算法实现：
1、初始化聚类的簇数K=0，已访问样本对象VT为空，分类结果result为空，异常索引noise为空。
2、遍历m的样本对象，判断其是否为核心对象，若不是核心对象，在索引对应的VT中标记为1，表示该对象已访问。如果是核心对象，则聚类簇数K=K+1，并且找到它的领域内的所有对象，标记已访问，遍历每一个邻域中的点，且在result矩阵中，所有密度可达的对象对应的索引位置都赋予K值。
3、对result矩阵取反，就是noise矩阵。也就是说result矩阵中，值为0的位置说明对应的样本不属于任何簇，是异常的。

源代码

总之，说的有点简单，具体见源码如下。

function [result, noncore,noise] = dbscan(SetOfPoints, eps, minPts) 

    cluster = 0;
    sizeOfCell = size(SetOfPoints,1);
    result = zeros(sizeOfCell,1);

%欧氏距离矩阵Dmatrix，对应于各个样本点之间的距离。
    Dmatrix = pdist2(SetOfPoints,SetOfPoints);
%初始化已访问样本和噪声样本。
    visited = false(sizeOfCell,1);
    noncore = false(sizeOfCell,1);


    for index = 1 : sizeOfCell
        if visited(index) == false
            visited(index) = true;%已访问标记。

            Neighbors = regionQuery(index, eps);
            if numel(Neighbors) < minPts
                noncore(index) = true;
            else
                cluster = cluster + 1;
                ExpandCluster(index, Neighbors, cluster, eps, minPts)%找出所有密度可达的点，归于cluster类。
            end  
        end
    end

    function ExpandCluster(index, Neighbors, cluster, eps, minPts)
        result(index) = cluster;
        temp = 1;

        while true
            neighbor = Neighbors(temp); 
            if visited(neighbor) == false
                visited(neighbor) = true;
                NeighborsSecond = regionQuery(neighbor, eps);
                if numel(NeighborsSecond) >= minPts
                   Neighbors = [Neighbors NeighborsSecond]; %密度可达点的索引扩展
                end
            end

            if result(neighbor) == 0
                result(neighbor) = cluster;

            end

            temp = temp + 1;
            if temp > numel(Neighbors)
                break;
            end
        end
    end

    function Neighbors = regionQuery(i, eps)
        Neighbors = find(Dmatrix(i,:)<=eps);%返回Dmatrix矩阵中第i行，小于eps的索引。
    end
 noise=~result;
end