基於密度的聚類之Dbscan算法

一.算法概述

　　DBSCAN(Density-Based Spatial Clustering of Applications with Noise)是一個比較有代表性的基於密度的聚類算法。與划分和層次聚類方法不同，它將簇定義為密度相連的點的最大集合，能夠把具有足夠高密度的區域划分為簇，並可在噪聲的空間數據庫中發現任意形狀的聚類（筆者認為是因為他不是基於距離的，基於距離的發現的是球狀簇）。

　　該算法利用基於密度的聚類的概念，即要求聚類空間中的一定區域內所包含對象（點或其他空間對象）的數目不小於某一給定閾值。DBSCAN算法的顯著優點是聚類速度快且能夠有效處理噪聲點和發現任意形狀的空間聚類。但是由於它直接對整個數據庫進行操作且進行聚類時使用了一個全局性的表征密度的參數，因此也具有兩個比較明顯的弱點：

　　（1）當數據量增大時，要求較大的內存支持I/O消耗也很大；

　　（2）當空間聚類的密度不均勻、聚類間距差相差很大時，聚類質量較差（有些簇內距離較小，有些簇內距離很大，但是Eps是確定的，所以，大的點可能被誤判斷為離群點或者邊界點，如果Eps太大，那么小距離的醋內，可能會包含一些離群點或者邊界點，KNN的k也存在同樣的問題）。

　　（1）與K-MEANS比較起來，不需要輸入要划分的聚類個數；

　　（2）聚類簇的形狀沒有偏倚（這個不明白啥意思）；

　　（3）可以在需要時輸入過濾噪聲的參數；

二.算法基本定義

三.算法描述

3.1 算法前提

　　DBSCAN算法基於一個事實：一個聚類可以由其中的任何核心對象唯一確定。等價可以表述為：任一滿足核心對象條件的數據對象p，數據庫D中所有從p密度可達的數據對象o所組成的集合構成了一個完整的聚類C，且p屬於C。

3.2 算法流程

 

四.算法實現

%% DBSCAN
clear all;
clc;
%% 導入數據集
% data = load('testData.txt');
data = randn(50,2);
% 定義參數Eps和MinPts
MinPts = 5;
Eps = epsilon(data, MinPts);
[m,n] = size(data);%得到數據的大小
x = [(1:m)' data];
[m,n] = size(x);%重新計算數據集的大小
types = zeros(1,m);%用於區分核心點1，邊界點0和噪音點-1
dealed = zeros(m,1);%用於判斷該點是否處理過,0表示未處理過
dis = calDistance(x(:,2:n));
number = 1;%用於標記類

%% 對每一個點進行處理
for i = 1:m
    %找到未處理的點
    if dealed(i) == 0
        xTemp = x(i,:);
        D = dis(i,:);%取得第i個點到其他所有點的距離
        ind = find(D<=Eps);%找到半徑Eps內的所有點    
        %% 區分點的類型     
        %邊界點
        if length(ind) > 1 && length(ind) < MinPts+1
            types(i) = 0;
            class(i) = 0;
        end
        %噪音點
        if length(ind) == 1
            types(i) = -1;
            class(i) = -1;
            dealed(i) = 1;
        end
        %核心點(此處是關鍵步驟)
        if length(ind) >= MinPts+1
            types(xTemp(1,1)) = 1;
            class(ind) = number;
            
            % 判斷核心點是否密度可達
            while ~isempty(ind)
                yTemp = x(ind(1),:);
                dealed(ind(1)) = 1;
                ind(1) = [];
                D = dis(yTemp(1,1),:);%找到與ind(1)之間的距離
                ind_1 = find(D<=Eps);
                
                if length(ind_1)>1%處理非噪音點
                    class(ind_1) = number;
                    if length(ind_1) >= MinPts+1
                        types(yTemp(1,1)) = 1;
                    else
                        types(yTemp(1,1)) = 0;
                    end
                    
                    for j=1:length(ind_1)
                       if dealed(ind_1(j)) == 0
                          dealed(ind_1(j)) = 1;
                          ind=[ind ind_1(j)];   
                          class(ind_1(j))=number;
                       end                    
                   end
                end
            end
            number = number + 1;
        end
    end
end
% 最后處理所有未分類的點為噪音點
ind_2 = find(class==0);
class(ind_2) = -1;
types(ind_2) = -1;

%% 畫出最終的聚類圖
hold on
for i = 1:m
    if class(i) == -1
        plot(data(i,1),data(i,2),'.r');
    elseif class(i) == 1
        if types(i) == 1
            plot(data(i,1),data(i,2),'+b');
        else
            plot(data(i,1),data(i,2),'.b');
        end
    elseif class(i) == 2
        if types(i) == 1
            plot(data(i,1),data(i,2),'+g');
        else
            plot(data(i,1),data(i,2),'.g');
        end
    elseif class(i) == 3
        if types(i) == 1
            plot(data(i,1),data(i,2),'+c');
        else
            plot(data(i,1),data(i,2),'.c');
        end
    else
        if types(i) == 1
            plot(data(i,1),data(i,2),'+k');
        else
            plot(data(i,1),data(i,2),'.k');
        end
    end
end
hold off

　　么么噠.............

%% 計算矩陣中點與點之間的距離
function [ dis ] = calDistance( x )
    [m,n] = size(x);
    dis = zeros(m,m);
    for i = 1:m
        for j = i:m
            %計算點i和點j之間的歐式距離
            tmp =0;
            for k = 1:n
                tmp = tmp+(x(i,k)-x(j,k)).^2;
            end
            dis(i,j) = sqrt(tmp);
            dis(j,i) = dis(i,j);
        end
    end
end

　　么么噠.............

function [Eps]=epsilon(x,k)
% Function: [Eps]=epsilon(x,k)
%
% Aim: 
% Analytical way of estimating neighborhood radius for DBSCAN
%
% Input: 
% x - data matrix (m,n); m-objects, n-variables
% k - number of objects in a neighborhood of an object
% (minimal number of objects considered as a cluster)

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

　　注意：prod是數組內元素的乘積，A^n是A*A*....*A，A.^n是A中每個元素的n次方。