DBSCAN密度聚类

DBSCAN

经典的密度聚类DBSCAN

基本概念

1. $\epsilon$邻域(eps邻域)
   $\epsilon$是一个距离阈值，对于样本集S中的任意$x_i$样本，与其距离小于等于$\epsilon$阈值的其他样本的集合，叫做样本$x_i$的$\epsilon$邻域
   记为：
   $N_{\epsilon }(x_i)$= {$x_j|dist(x_i,x_j)\le \epsilon$}

2. 核心对象（点）
   若$N_{\epsilon }(x_i)$内样本点数大于等于$MinPts$阈值，即
   |$N_{\epsilon }(x_i)$| $\ge MinPts$
   则$x_i$称为核心点，否则称为边界点

3. 密度直达
   $x_j$在$N_{\epsilon }(x_i)$，且$x_i$为核心点，则$x_j$可由$x_i$密度直达，反之不一定成立($x_j$不一定是核心点)
   只有核心点之间满足对称性

4. 密度可达
   若有一个样本点序列$p_1,p_2,p_3,...p_n$，$p_1=x_i,p_n=x_j$，$p_{i+1}$可由$p_i$密度直达，则$x_j$可有$x_i$密度可达(传递性)

5. 密度相连
   若存在一个核心点$o$，使得$x_i$和$x_j$均由$o$密度可达
   则称$x_i$和$x_j$密度相连
   

   i. 所有密度相连的点分为一簇，同一簇的点必定密度相连
   ii. 一个核心对象所有密度可达的点集合，分为一簇

DBSCAN基本思想

算法流程

python实现

"""
    python 实现 DBSCAN algorithm
    
"""
import numpy as np
from matplotlib import pyplot as plt

import logging
logging.basicConfig(level=logging.INFO,format="%(asctime)s %(message)s",filename="log.txt",filemode="a")


#加载数据并归一化
def load_data(data_path="data.txt"):
    """
        data_path:数据文件路径，str,default-->"data.txt"
        return x_scale:归一化的np 2dim数组
    """
    #加载数据
    x=np.loadtxt(data_path)
    
    #计算数据的mean,std
    x_mean=np.mean(x,axis=0)
    x_std=np.std(x,axis=0)
    
    #归一化
    x_scale=(x-x_mean)/x_std
    
    return x_scale
    #return x


#计算欧式距离
#这里做稍微的改进
def l2(x,y):
    """
    x:样本点向量 1dim 
    y:样本点向量 1dim or 2dim
    return:距离标量/距离1dim 数组
    """
    #y is 1dim
    if y.ndim==1:
        return np.sqrt(np.sum(np.power(x-y,2)))
    
    #y is 2dim
    elif y.ndim==2:
        return np.sqrt(np.sum(np.power(x-y,2),axis=1))
    
    #other 
    else:
        pass
    

#判断是否核心点
def is_kernel(i,S,eps=0.5,min_samples=5):
    """
        判断当前样本x是否核心对象
        i:单个样本索引，int
        S:样本集,2dim数组
        eps:eps阈值,default 0.5
        min_samples:MinPts阈值 default 5
        return (True/False,eps_field)
    """
    #计算当前样本点i与S样本集内所有样本的距离
    d = l2(S[i],S)
    logging.info("Current distance of S[%d] with S:\n%r"%(i,d))
    
    if i == 1:
        print("Current distance of S[%d] with S:\n%r"%(i,d))
    
    #与eps 距离比较
    mask = d <= eps
    
    #获取当前样本i的eps邻域
    eps_field = set(np.nonzero(mask)[0])
    
    logging.info("eps field of current S[%d]:\n%r"%(i,eps_field))
    if i == 1:
        print("eps field of current S[%d]:\n%r"%(i,eps_field))
    
    #判断x是否核心点
    if len(eps_field) >= min_samples:
        return True,eps_field
    
    else:
        return False,eps_field


#检测当前核心点  所有密度可达的核心点
class GetKernels(object):
    
    def __init__(self,kernels_dict):
        """
            kernels_dict: 所有的核心点及其eps 邻域(集合) 的字典
            
        """
        self.kernels_dict = kernels_dict
        
        
    def get_reachable_kernels(self,directs,kernels):
        
        """
            directs:密度直达的核心点 索引集合，set  
            kernels: 所有未分簇的核心点索引集合,set
            return result
                   所有密度可达的核心点集合
        """
        self.all_reachable_kernels = set()
        
        while directs:
            
            #取出一个核心点
            o = directs.pop()
            
            #将o 放入预定的集合
            self.all_reachable_kernels.add(o)
            
            #从所有未分簇的核心点集合中 去除 o核心点
            kernels.remove(o)
            
            #检查o核心点的eps邻域内 是否有其他的核心点
            if not self.kernels_dict.get(o)&kernels :
                continue
            else:
                directs |= self.kernels_dict.get(o)&kernels
        
        return self.all_reachable_kernels
    


def dbscan(X,eps=0.5,min_samples=5):
    """
        X:样本集，2dim数组
        eps: eps邻域距离阈值，default 0.5
        min_samples: MinPts阈值，default 5
        
        return (y_pred,kernel_points)
        y_pred: 聚类的样本标记，一维数组
        kernels: 核心点集合
    """
    #获取数据行，列
    m,n = X.shape
    
    #1 初始化
    #核心点及其eps邻域
    kernels_dict = {}
    
    #标记值
    k = 0
    
    #初始化的标记
    y_pred = np.ones(m)*(-1)
    
    #2 计算出所有的核心点
    for i in range(m):
        
        #判断是否核心点
        bool_,eps_field=is_kernel(i,X,eps,min_samples)
        
        if bool_:
            logging.info("S[{i}] is kernel point".format(i=i))
            kernels_dict[i]=eps_field
        else:
            logging.info("S[{i}] not kernel point".format(i=i))
            continue
        
    #所有未分簇的核心点
    kernels = set(kernels_dict.keys())
    
    
    #3 取未分簇的核心点，并找到它所有的密度可达的核心点
    #实例化
    get_reachable = GetKernels(kernels_dict)
    while kernels:
        
        #取一个未分簇的核心点
        o_i = kernels.pop()
        
        #获取当前核心点的eps邻域
        c_i = kernels_dict.get(o_i)
        
        #当前核心点的密度直达的核心点
        directs = c_i & kernels
        
        #当前核心点 eps邻域内有其他核心点
        if directs:
            #获取当前核心点所有密度可达的核心点
            all_reachable = get_reachable.get_reachable_kernels(directs,kernels)
            
            for i in all_reachable:
                
                #所有核心点的eps邻域 取并集
                c_i |= kernels_dict.get(i)
            
        
        #分簇
        y_pred[np.array(list(c_i))] = k
        
        k += 1
        
    return y_pred,list(kernels_dict.keys())
        
        


if __name__=="__main__":
    
    X = load_data("data.txt")
#    from sklearn.datasets import make_circles
#    X,y = make_circles(n_samples=300,noise=0.1,factor=0.3)
    y_pred,kernel_points = dbscan(X,eps=0.5,min_samples=2)
    print("cluster result:",np.unique(y_pred))
    print("kernels:",kernel_points)
    
    
    
    #可视化
    #
    plt.rcParams["font.sans-serif"] = ["SimHei"]
    plt.rcParams["axes.unicode_minus"] = False
    
    fig = plt.figure()
    plt.title("DBSCAN_python聚类",fontsize=20)
    plt.scatter(X[:,0],X[:,1],s=300,c="lightblue",label="Samples",cmap="cool")
    
    plt.scatter(X[kernel_points,0],X[kernel_points,1],s=100,linewidths=3,marker="*",label="kernels")
    plt.scatter(X[y_pred==-1,0],X[y_pred==-1,1],s=100,marker="s",label="noise")
    
    #聚类中心点
    c0 = np.mean(X[y_pred==0],axis=0)
    c1 = np.mean(X[y_pred==1],axis=0)
    c_1 = np.mean(X[y_pred==-1],axis=0)
    
    plt.scatter(c0[0],c0[1],s=300,c="white",edgecolors=None,linewidths=2,marker="o")
    plt.scatter(c0[0],c0[1],s=300,c="r",edgecolors=None,linewidths=2,marker="$0$")
    
    plt.scatter(c1[0],c1[1],s=300,c="white",edgecolors=None,linewidths=2,marker="o")
    plt.scatter(c1[0],c1[1],s=300,c="r",edgecolors=None,linewidths=2,marker="$1$")
    
    plt.scatter(c_1[0],c_1[1],s=500,c="white",edgecolors=None,linewidths=2,marker="o")
    plt.scatter(c_1[0],c_1[1],s=500,c="r",edgecolors=None,linewidths=2,marker="$-1$")
    plt.grid()
    plt.legend(loc="best",frameon=True,framealpha=0.7)
    plt.show()
#    bool_,eps_field=is_kernel(x[0],x)
#    if bool_:
#        logging.info("x[{i}] is kernel point".format(i=0))
#    else:
#        logging.info("x[{i}] not kernel point".format(i=0))
#    

sklearn库的实现

from sklearn.cluster import DBSCAN
import numpy as np
samples = np.loadtxt("data.txt")

clustering = DBSCAN(eps=6, min_samples=3).fit(samples)

import matplotlib.pyplot as plt
plt.scatter(samples[:,0],samples[:,1],c=clustering.labels_+1.5,linewidths=np.power(clustering.labels_+1.5, 2))
plt.show()

DBSCAN优点及缺点

优点

1. 适用凸、非凸的稠密数据集
   不像KMeans/Birch之类的，只适用凸集
2. 不需事先指定聚类簇数$k$
3. 聚类完成可以识别噪声点 ，对异常点不敏感

缺点

1. 非稠密数据，距离差异大的样本集，都不适用DBSCAN
2. 样本集大时，聚类收敛时间较长，使用kd_tree,ball_tree增加空间复杂度
3. $\epsilon 和MinPts$需联合调参，比较复杂，DBSCAN对这俩个参数比较敏感

$\epsilon 和MinPts$选参

1. 交叉验证
   取$\epsilon 和MinPts$合理的数值组合，使DBSCAN聚类的评价指标最优，比如轮廓系数，DB指数等
2. 网格搜索
3. 手肘法，根据MinPts，确定$\epsilon$值
   i. MinPts取初始值$k=n+1$，即特征数加$1$
   ii.迭代每个样本$x_i$，计算与样本集S内其他样本的距离，取第$k$近的距离$k_{dist}$，如$[0,1,2,2,3,4,5]$中第三近的距离为3，0表示样本与自己的距离。将所有样本的$k_{dist}$放入列表$L$
   iii. 对$L$从大到小排序，作线性可视化，取拐点处的$k_{dist}$作为$\epsilon$值
   python 代码手肘法实现：

"""
    decide 'eps' according to 'MinPts'
    'MinPts'=n+1
"""
import numpy as np
from matplotlib import pyplot as plt
from sklearn.preprocessing import StandardScaler

class SelectEPS(object):
    
    def __init__(self,datapath):
        self.x=np.loadtxt(datapath)
        self.x_scale=StandardScaler().fit_transform(self.x)
        self.m,self.n=x.shape
        
    def select_eps(self):
        
        #初始化MinPts
        min_pts=self.n+1
        k_dist_list=[]
        #遍历样本点
        for i in range(self.m):
            p=self.x[i]
            p_to_other=self.l2(p)
        
            print("Current sample:",p)
            print("Distance of p with others:",p_to_other)
            k_dist_list.append(self.get_k_th_dist(min_pts,p_to_other))
            
        print(k_dist_list)
        
        k_dist_list_sort=np.sort(k_dist_list)[::-1]
        plt.scatter(range(self.m),k_dist_list_sort,marker="o")
        
        plt.xlabel("index")
        plt.ylabel("k_dist        ",rotation="horizontal")
        plt.title("k_dist of all samples")
        plt.show()
    
    def l2(self,x):
        """
            计算单个样本与样本集S内其他样本的距离
            x:single sample,1dim
            return:distance,1dim
        """
        return np.sqrt(np.sum((x-self.x)**2,axis=1))
        
    def get_k_th_dist(self,k,dist_list):
        """
            从包含0的几个距离中取最近的k个
        """
        #排序
        sort_dist=np.sort(dist_list)
        
        #取出第k近的值
        i=0
        while i<k:
            if sort_dist[i]==0:
                k+=1
                i+=1
                if i>len(sort_dist)-1:
                    return "no result"
                continue
            elif sort_dist[i]==sort_dist[i-1]:
                k+=1
                i+=1
                if i>len(sort_dist)-1:
                    return "no result"
                continue
            else:
                i+=1
                if i>len(sort_dist)-1:
                    return "no result"
        
        #取到除0外的第k个最近距离
        return sort_dist[i-1]
            
          
if __name__=="__main__":
    eps=SelectEPS("data.txt")
    eps.select_eps()
    
#    d=eps.get_k_th_dist(1,[0,0,0,0,0,1,1,1,2,2,2,3,4,4])
#    print(d)