簡單易學的機器學習算法—基于密度的聚類算法DBSCAN

一、基于密度的聚類算法的概述

    我想了解下基于密度的聚類算法，熟悉下基于密度的聚類算法與基于距離的聚類算法，如K-Means算法之間的區別。
    基于密度的聚類算法主要的目標是尋找被低密度區域分離的高密度區域。與基于距離的聚類算法不同的是，基于距離的聚類算法的聚類結果是球狀的簇，而基于密度的聚類算法可以發現任意形狀的聚類，這對于帶有噪音點的數據起著重要的作用。
二、DBSCAN算法的原理
1、基本概念
    DBSCAN(Density-Based Spatial Clustering of Application with Noise)是一種典型的基于密度的聚類算法，在DBSCAN算法中將數據點分為一下三類：
核心點。在半徑Eps內含有超過MinPts數目的點
邊界點。在半徑Eps內點的數量小于MinPts，但是落在核心點的鄰域內
噪音點。既不是核心點也不是邊界點的點
在這里有兩個量，一個是半徑Eps，另一個是指定的數目MinPts。
    一些其他的概念
Eps鄰域。簡單來講就是與點p的距離小于等于Eps的所有的點的集合，可以表示為。
直接密度可達。如果p在核心對象q的Eps鄰域內，則稱對象p從對象q出發是直接密度可達的。
密度可達。對于對象鏈：，是從關于Eps和MinPts直接密度可達的，則對象是從對象關于Eps和MinPts密度可達的。
2、算法流程

(流程)
三、實驗仿真
    在實驗中使用了兩個測試數據集，數據集的原始圖像如下：

(數據集1)

(數據集2)
數據集1相對比較簡單。顯然我們可以發現數據集1共有兩個類，數據集2有四個類，下面我們通過DBSCAN算法實現數據點的聚類：

MATLAB代碼
主程序


[plain] view plain copy 在CODE上查看代碼片派生到我的代碼片
%% DBSCAN  
clear all;  
clc;  
 
%% 導入數據集  
% data = load('testData.txt');  
data = load('testData_2.txt');  
 
% 定義參數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  

距離計算函數
[plain] view plain copy 在CODE上查看代碼片派生到我的代碼片
%% 計算矩陣中點與點之間的距離  
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  

epsilon函數
[plain] view plain copy 在CODE上查看代碼片派生到我的代碼片
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);  

最終的結果

(數據集1的聚類結果)

(數據集2的聚類結果)
在上面的結果中，紅色的點代表的是噪音點，點代表的是邊界點，十字代表的是核心點。不同的顏色代表著不同的類。