簡單易學(xué)的機(jī)器學(xué)習(xí)算法—基于密度的聚類算法DBSCAN

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

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

(流程)
三、實(shí)驗(yàn)仿真
    在實(shí)驗(yàn)中使用了兩個(gè)測試數(shù)據(jù)集，數(shù)據(jù)集的原始圖像如下：

(數(shù)據(jù)集1)

(數(shù)據(jù)集2)
數(shù)據(jù)集1相對比較簡單。顯然我們可以發(fā)現(xiàn)數(shù)據(jù)集1共有兩個(gè)類，數(shù)據(jù)集2有四個(gè)類，下面我們通過DBSCAN算法實(shí)現(xiàn)數(shù)據(jù)點(diǎn)的聚類：

MATLAB代碼
主程序


[plain] view plain copy 在CODE上查看代碼片派生到我的代碼片
%% DBSCAN  
clear all;  
clc;  
 
%% 導(dǎo)入數(shù)據(jù)集  
% data = load('testData.txt');  
data = load('testData_2.txt');  
 
% 定義參數(shù)Eps和MinPts  
MinPts = 5;  
Eps = epsilon(data, MinPts);  
 
[m,n] = size(data);%得到數(shù)據(jù)的大小  
 
x = [(1:m)' data];  
[m,n] = size(x);%重新計(jì)算數(shù)據(jù)集的大小  
types = zeros(1,m);%用于區(qū)分核心點(diǎn)1，邊界點(diǎn)0和噪音點(diǎn)-1  
dealed = zeros(m,1);%用于判斷該點(diǎn)是否處理過,0表示未處理過  
dis = calDistance(x(:,2:n));  
number = 1;%用于標(biāo)記類  
 
%% 對每一個(gè)點(diǎn)進(jìn)行處理  
for i = 1:m  
    %找到未處理的點(diǎn)  
    if dealed(i) == 0  
        xTemp = x(i,:);  
        D = dis(i,:);%取得第i個(gè)點(diǎn)到其他所有點(diǎn)的距離  
        ind = find(D<=Eps);%找到半徑Eps內(nèi)的所有點(diǎn)  
          
        %% 區(qū)分點(diǎn)的類型  
          
        %邊界點(diǎn)  
        if length(ind) > 1 && length(ind) < MinPts+1  
            types(i) = 0;  
            class(i) = 0;  
        end  
        %噪音點(diǎn)  
        if length(ind) == 1  
            types(i) = -1;  
            class(i) = -1;  
            dealed(i) = 1;  
        end  
        %核心點(diǎn)(此處是關(guān)鍵步驟)  
        if length(ind) >= MinPts+1  
            types(xTemp(1,1)) = 1;  
            class(ind) = number;  
              
            % 判斷核心點(diǎn)是否密度可達(dá)  
            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%處理非噪音點(diǎn)  
                    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  
 
% 最后處理所有未分類的點(diǎn)為噪音點(diǎn)  
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  

距離計(jì)算函數(shù)
[plain] view plain copy 在CODE上查看代碼片派生到我的代碼片
%% 計(jì)算矩陣中點(diǎn)與點(diǎn)之間的距離  
function [ dis ] = calDistance( x )  
    [m,n] = size(x);  
    dis = zeros(m,m);  
      
    for i = 1:m  
        for j = i:m  
            %計(jì)算點(diǎn)i和點(diǎn)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函數(shù)
[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);  

最終的結(jié)果

(數(shù)據(jù)集1的聚類結(jié)果)

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