挑子學習筆記：DBSCAN算法的python實現

轉載請標明出處：https://www.cnblogs.com/tiaozistudy/p/dbscan_algorithm.html

DBSCAN（Density-Based Spatial Clustering of Applications with Noise）聚類算法，是一種基於高密度連通區域的、基於密度的聚類算法，能夠將具有足夠高密度的區域划分為簇(Cluster)，並在具有噪聲的數據中發現任意形狀的簇。DBSCAN算法通過距離定義出一個密度函數，計算出每個樣本附近的密度，從而根據每個樣本附近的密度值來找出那些樣本相對比較集中的區域，這些區域就是我們要找的簇。

1. DBSCAN算法的基本原理

其它聚類方法大都是基於對象之間的距離進行聚類，聚類結果是球狀的簇。DBSCAN 算法利用簇的高密度連通性，尋找被低密度區域分離的高密度區域，可以發現任意形狀的簇，其基本思想是:對於一個簇中的每個對象，在其給定半徑的領域中包含的對象不能少於某一給定的最小數目。

DBSCAN算法中有兩個重要參數：$\varepsilon$表示定義密度時的鄰域半徑，$M$ 表示定義核心點時的閾值。

考慮數據集合$X= \{x^{(1)}, x^{(2)},...,x^{(n)} \}$，引入以下概念與記號。

1. $\varepsilon$鄰域

設$x \in X$，稱
\begin{equation*}
N_{\varepsilon}(x) = \{ y \in X: d(y,x) \le \varepsilon \}
\end{equation*}
為$X$的$\varepsilon$鄰域。顯然，$x \in N_{\varepsilon}(x)$。

為了簡單起見，節點$x^{(i)}$與其下標$i$一一對應，引入記號
\begin{equation*}
N_{\varepsilon}(i) = \{ j: d(y^{(j)},x^{(i)}) \le \varepsilon; \; y^{(j)},x^{(i)} \in X \}
\end{equation*}

2. 密度

設$x \in X$，稱$\rho(x) = |N_{\varepsilon}(x)|$為$x$的密度。密度是一個整數，且依賴於半徑$\varepsilon$。

3. 核心點

設$x \in X$，若$\rho(x) \ge M$，則稱$x$為$X$的核心點。記由$X$中所有核心點構成的集合為$X_c$，並記$X_{nc}=X-X_c$表示$X$中所有非核心點構成的集合。

4. 邊界點

若$x \in X_{nc}$，且$\exists y \in X$，滿足$y \in N_{\varepsilon}(x) \bigcap X_c$，即$X$的非核心點$x$的$\varepsilon$鄰域中存在核心點，則稱$x$ 為$X$的邊界點。記由$X$中所有邊界點構成的集合為$X_{bd}$。

此外，邊界點也可以這么定義，若$x \in X_{nc}$，且$x$落在某個核心點的$\varepsilon$鄰域內，則稱$x$為$X$的邊界點。一個邊界點可能同時落入一個或多個核心點的$\varepsilon$ 鄰域內。

5. 噪聲點

記$X_{noi} = X - (X_c \bigcup X_{bd})$，若$x \in X_{noi}$，則稱$x$為噪音點。

至此，我們嚴格給出了核心點、邊界點和噪音點的數學定義，且滿足$X = X_c \bigcup X_{bd} \bigcup X_{noi}$.

圖1：核心點、邊界點和噪聲點

直觀地說，核心點對應稠密區域內部的點，邊界點對應稠密區域邊緣的點，而噪音點對應稀疏區域中的點。

數據集通過聚類形成的子集是簇。核心點位於簇的內部，它確定無誤地屬於某個特定的簇；噪音點是數據集中的干擾數據，它不屬於任何一個簇；邊界點是一類特殊的點，它位於一個或幾個簇的邊緣地帶，可能屬於一個簇，也可能屬於另外一個簇，其歸屬並不明確。

6. 直接密度可達

設$x,y \in X$. 若滿足$x \in X_c$，則稱$y$是$x$從直接密度可達的。

7. 密度可達

設$p^{(1)}, p^{(2)},..., p^{(m)} \in X$，其中$m \ge 2$。若它們滿足：$p^ {(i+1)}$ 是從$p^{(i)}$直接密度可達的，其中$i = 1,2,...,m-1$，則稱$p^ {(m)}$ 是從$p^{(1)}$ 中密度可達的。

7.1. 當$m = 2$時，密度可達即為直接密度可達。實際上，密度可達是直接密度可達的傳遞閉包。

7.2. 密度可達關系不具有對稱性。若$p^{(m)}$是從$p^{(1)}$密度可達的，那么$p^ {(1)}$ 不一定是從$p^{(m)}$密度可達的。根據上述定義可知，$p^{(1)}, p^{(2)}, ..., p^{(m-1)}$必須為核心點，而$p^{(m)}$可以是核心點，也可以是邊界點。當$p^ {(m)}$是邊界點時，$p^{(1)}$一定不是從$p^{(m)}$密度可達的。

8. 密度相連

設$x,y,z \in X$，若$y$和$z$均是從$x$密度可達的，則稱$y$和$z$是密度相連的。顯然，密度相連具有對稱性。

9. 簇（cluster）

非空集合$C \subset X$，如果$C$滿足：對於$x,y \in X$

• 若$x \in C$，且$y$是從$x$密度可達的，則$y \in C$，
• 若$x \in C$，$y \in C$，則$x,y$是密度相連的。

則稱$C$是$X$的一個簇。

DBSCAN 算法基於以下一個基本事實：對於任一核心點$x$，數據集$X$中所有從$x$ 密度可達的數據點可以構成一個完整的簇$C$，且$x \in C$。其核心思想描述如下:從某個選定的核心點出發，不斷向密度可達的區域擴張，從而得到一個包含核心點和邊界點的最大化區域，區域中任意兩點密度相連。

2. DBSCAN算法的實現

《數據挖掘概念與技術》給出的算法偽代碼如下:

考慮數據集合$X= \{x^{(1)}, x^{(2)},..., x^{(n)} \}$。DBSCAN算法的目標是將數據集合$X$分成$K$個簇及噪聲點集合，其中$K$也是由算法得到，為此，引入簇的標記數組
\begin{equation*}
m_i =
\begin{cases}
j, & \text{若}x^{(i)}\text{屬於第}j\text{個簇}; \\
-1, & \text{若}x^{(i)}\text{為噪聲點}
\end{cases}
\end{equation*}
DBSCAN算法的目標就是生成標記數組$m_i, \; i=1,...,n$.

為了保證可以更有效地實現算法1中第3句隨機選擇一個unvisited對象$p$，設計了一個數據結構visitlist，其中包含兩個列表visitedlist和unvisitedlist，分別用於存儲已訪問的點和未訪問的點，每次從unvisitedlist 中取點可以保證每次取到的點都是未訪問過的點，實現代碼如下：

代碼1：visitlist數據結構

# visitlist類用於記錄訪問列表
# unvisitedlist記錄未訪問過的點
# visitedlist記錄已訪問過的點
# unvisitednum記錄訪問過的點數量
class visitlist:
    def _init_(self, count=0):
        self.unvisitedlist=[i for i in range(count)]
        self.visitedlist=list()
        self.unvisitednum=count

    def visit(self, pointId):
        self.visitedlist.append(pointId)
        self.unvisitedlist.remove(pointId)
        self.unvisitednum -= 1

DBSCAN算法實現代碼如下:

代碼2：DBSCAN算法實現

import numpy as np
import matplotlib.pyplot as plt
import math
import random

def  dist(a, b):
    # 計算a,b兩個元組的歐幾里得距離
    return math.sqrt(np.power(a-b, 2).sum())

def my_dbscanl(dataSet, eps, minPts):
    # numpy.ndarray的 shape屬性表示矩陣的行數與列數
    nPoints = dataSet.shape[0]
    # (1)標記所有對象為unvisited
    # 在這里用一個類vPoints進行買現
    vPoints = visitlist(count=nPoints)
    # 初始化簇標記列表C,簇標記為 k
    k = -1
    C = [-1 for i in range(nPoints)]
    while(vPoints.unvisitednum > 0):
        # (3)隨機上選擇一個unvisited對象p
        P = random.choice(vPoints.unvisitedlist)
        # (4)標記p為visited
        vPoints.visit(p)
        # (5)if p的$\varepsilon$-鄰域至少有MinPts個對象
        # N是p的$\varepsilon$-鄰域點列表
        N = [i for i in range(nPoints) if dist(dataSet[i], dataSet[p])<= eps]
        if  len(N) >= minPts:
            # (6)創建個新簇C，並把p添加到C
            # 這里的C是一個標記列表，直接對第p個結點進行賦植
            k += 1
            C[p]=k
            # (7)令N為p的ε-鄰域中的對象的集合
            # N是p的$\varepsilon$-鄰域點集合
            # (8) for N中的每個點p'
            for p1 in N:
                # (9) if p'是unvisited
                if p1 in vPoints.unvisitedlist:
                    # (10)標記p’為visited
                    vPoints.visit(p1)
                    # (11) if p'的$\varepsilon$-鄰域至少有MinPts個點，把這些點添加到N
                    # 找出p'的$\varepsilon$-鄰域點，並將這些點去重添加到N
                    M=[i for i in range(nPoints) if dist(dataSet[i], \
                        dataSet[p1]) <= eps]
                    if len(M) >= minPts:
                        for i in M:
                            if i not in N:
                                N.append(i)
                    # (12) if p'還不是任何簇的成員，把P'添加到C
                    # C是標記列表，直接把p'分到對應的簇里即可
                    if  C[p1] == -1:
                        C[p1]= k
        # (15)else標記p為噪聲
        else:
            C[p]=-1

    # (16)until沒有標t己為unvisitedl內對象
    return C

 利用sklearn生成數據集，共2500條數據，並利用matplotlib畫出散點圖，代碼如下:

代碼3：生成數據集

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

X1, Y1 = datasets.make_circles(n_samples=2000, factor=0.6, noise=0.05,
                               random_state=1)
X2, Y2 = datasets.make_blobs(n_samples=500, n_features=2, centers=[[1.5,1.5]],
                             cluster_std=[[0.1]], random_state=5)

X = np.concatenate((X1, X2))
plt.figure(figsize=(12, 9), dpi=80)
plt.scatter(X[:,0], X[:,1], marker='.')
plt.show()

 

圖2：數據集散點圖

設置參數Eps=0.1, MinPts=10，聚類結果如下圖：

 

圖3：聚類結果

3. 利用KD樹進行優化

 

KD樹（K-Dimensional Tree），是一種分割k維數據空間的數據結構，是二叉搜索樹在多維條件下的推廣。主要應用於多維空間關鍵數據的搜索。KD樹的介紹見：https://www.jianshu.com/p/ffe52db3e12b，不贅述。

利用scipy實現KD樹的構造和查詢，對代碼2的算法進行改進，代碼如下:

代碼4：DBSCAN算法的優化實現

import numpy as np
import matplotlib.pyplot as plt
import math
import random
from scipy.spatial import KDTree

def my-dbscan2(dataSet, eps, minPts):
    # numpy.ndarray的 shape屬性表示矩陣的行數與列數
    # 行數即表小所有點的個數
    nPoints = dataSet.shape[0]
    # (1) 標記所有對象為unvisited
    # 在這里用一個類vPoints進行實現
    vPoints = visitlist(count=nPoints)
    # 初始化簇標記列表C，簇標記為 k
    k = -1
    C = [-1 for i in range(nPoints)]
    # 構建KD-Tree，並生成所有距離<=eps的點集合
    kd = KDTree(X)
    while(vPoints.unvisitednum>0):
        # (3) 隨機選擇一個unvisited對象p
        p = random.choice(vPoints.unvisitedlist)
        # (4) 標t己p為visited
        vPoints.visit(p)
        # (5) if p 的$\varepsilon$-鄰域至少有MinPts個對象
        # N是p的$\varepsilon$-鄰域點列表
        N = kd.query_ball_point(dataSet[p], eps)
        if len(N) >= minPts:
            # (6) 創建個一個新簇C，並把p添加到C
            # 這里的C是一個標記列表，直接對第p個結點進行賦值
            k += 1
            C[p] = k
            # (7) 令N為p的$\varepsilon$-鄰域中的對象的集合
            # N是p的$\varepsilon$-鄰域點集合
            # (8) for N中的每個點p'
            for p1 in N:
                # (9) if p'是unvisited
                if p1 in vPoints.unvisitedlist:
                    # (10) 標記p'為visited
                    vPoints.visit(p1)
                    # (11) if p'的$\varepsilon$-鄰域至少有MinPts個點，把這些點添加到N
                    # 找出p'的$\varepsilon$-鄰域點，並將這些點去重新添加到N
                    M = kd.query_ball_point(dataSet[p1], eps)
                    if len(M) >= minPts:
                        for i in M:
                            if i not in N:
                                N.append(i)
                    # (12) if p'還不是任何簇的成員，把p'添加到c
                    # C是標記列表，直接把p'分到對應的簇里即可
                    if C[p1] == -1
                        C[p1] = k
        # (15) else標記p為噪聲
        else:
            C[p1] = -1

    # (16) until沒有標記為unvisited的對象
    return C

 

 以代碼3中生成的2500條數據作為測試，比較優化前后的算法性能

import time
start = time.time()
C1 = my_dbscanl(X, 0.1, 10)
end = time.time()
print "`運行時間`:", end - start
plt.scatter(X[:, 0], X[:, 1], c=C1, marker='.')
plt.show()
>>> `運行時間：`29.1249849796

圖4：優化前算法結果

import time
start = time.time()
C2 = my_dbscan2(X, 0.1, 10)
end = time.time()
print "運行時間:", end - start
plt.scatter(X[:, 0], X[:, 1], c=C2, marker='.')
plt.show()
>>> 運行時間：4.72340583801

 

 

圖5：優化后算法結果

可以看到優化后的算法運行時間從29.12s降到了4.72s，優化的效果非常明顯。

4. 后記

上文僅僅是對DBSCAN算法的思想與實現進行了簡略摘要，是學習算法的一個過程。算法的學習還比較粗劣和淺層，在實踐應用中上述代碼並不實用。如果需要使用DBSCAN的算法求解聚類問題，建議使用sklearn自帶的DBSCAN函數。以代碼3中生成數據為例：

# DBSCAN eps = 0.1, MinPts = 10
import time
from sklearn.cluster import DBSCAN
start = time.time()
C = DBSCAN(eps=0.1, min_pts=10).
end = time.time()
print "運行時間：", end - start
plt.scatter(X[:, 0], X[:, 1], c=C, marker='.')
plt.show()
>>> 運行時間：0.0240921974182