半监督密度聚类算法（1）

本算法为小编自己提出，只是实现了初步的想法，但是实际的实验效果并不是很好，在开始半监督话题之前，有几点是读者需要了解的，首先是半监督算法，然后是关于密度聚类算法的一些知识。
什么是半监督算法，半监督算法有什么用？小编比较懒，长话短说，半监督算法就是穷人没钱买数据，最大化的利用极少量的样本，实现和监督算法相近的性能，本篇文章的重点放在密度聚类算法上。本文实现的半监督算法是基于一篇2014年Science文章中的密度聚类算法。该算法的能中science，原因是简单而高效。具体的论文可以从我的博客里面下载论文地址

大体的思想是：通过计算全局（当然也可以优化到局部）所有Sample相互之间的距离，并对某个Sample与其他所有点之间的距离进行排序，通过一个threshold对距离进行分割获取有效后，根据该Sample有效距离内的其他Sample的数量来推算该点的密度，再计算低密度点到高密度点的最小距离（最高密度点取最大距离），通过密度和最小距离这两个参数构建直角坐标系，在这个直角坐标系中最右上角的点即聚类中心点。

该算法有一个缺点就是不能自动地寻找聚类中心，这个方法的想法是很简单而且很新颖的，但是涉及到大量的距离的计算，当数据量比较大的时候或者单个数据的维度很高时，该算法的计算量还是很大的。虽然作者本人在文章中提到了一种自动寻找聚类中心的方法，但是经过小编的实现，该自动寻找聚类中心的方法只能适用于小部分的数据集。

这篇文章的重点是提出了一种密度聚类的思想，算法主要分为两步：首先是通过“决策图”人工选取聚类中心，然后，将剩下的点根据距离和密度的公式分配到相应的聚类中心。那么决策图是怎么得到的呢？在实现决策图之前，需要提到文章中的两个重要的概念，一个是密度，一个是距离。密度的定义如下所示：

其中

dc d c 是截断距离，文中为所有点的相互距离中由小到大排列占总数2%的位置距离数值。什么是截断距离看我下面的解释吧。

dij d i j 是点

xi x i 到

yi y i 的距离，这里的点到点的距离的定义可以有很多种，我们这里使用的距离是欧式距离。假设

xi=(xi1,xi2,xi3,...,xin) x i = ( x i 1 , x i 2 , x i 3 , . . . , x i n ) ,

yi=(yi1,yi2,yi3,...,yin) y i = ( y i 1 , y i 2 , y i 3 , . . . , y i n ) ,那么欧式距离可以有下面的式子得出：

dij=(xi1−yi1)2+(xi2−yi2)2+(xi3−yi3)2+...+(xin−yin)2−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−√ d i j = ( x i 1 − y i 1 ) 2 + ( x i 2 − y i 2 ) 2 + ( x i 3 − y i 3 ) 2 + . . . + ( x i n − y i n ) 2

从上面的式子就可以看出，一旦单个样本的维度很大，距离的计算量也不是一般的大。 χ χ 表示的是二值运算， χ(a) χ ( a ) 当括号中的a值大于0是，则 χ(a)=0 χ ( a ) = 0 ，当a值小于0时， χ(a)=1 χ ( a ) = 1 ，也就是说，在算法的初始阶段，你首先随机选取一个点，然后计算样本中所有的点到该点的距离，如果任意一个点到这个点的距离小于我们设定的截断距离， χ(a) χ ( a ) 的值就会记为1，这样的话，我们遍历整个样本，计算所有点到选取点的距离小于截断距离点的个数就得到了我们随机选取的这个点的密度。密度其实想起来很简单，就是你选定一个点和一个截断距离了，这个截断距离可以理解为这个点的邻域半径，你找到所有样本到选取点距离小于这个领域半径的所有点的个数，就是这个点的密度值。这个点的同样的方法，可以计算出样本中每个点的密度。就像下面这张图：

设定一个邻域半径，这个邻域半径内的点的个数就是该点的密度。从上面的b图可以看出，中间的共色部分的密度比较大，边缘点的密度比较小，也就是说密度越大的点越有可能是聚类中心。接下来看一下下面的这张图，下图中的左图：

可以看到左边图显示的就是点的分布图，还没有经过聚类，我们是需要通过算法寻找聚类中心，左图的标号从小到大依次表示的是密度的从大到小，从左图的标号可以看出，大部分的密度大的点都集中在左簇中，右簇中密度最大的点是10，如果按照我们的之前的想法，密度越大的点越有可能是聚类中心，那么左簇中小于10的点是聚类中心的可能性比右簇中的所有点是聚类中心的可能性都大，但是我们可以看到，左簇只可能有一个聚类中心，所以单纯地以密度为度量来选取聚类中心是不够，所以在该论文中引入了距离 δ δ , δ δ 表示的是所有比 i i 点密度高的所有点的最近距离表示 δi δ i ,最近距离，最近，距离，重要的事情说三遍。 对于最大密度的点其为所有点与点距离的最大值 δi=maxj(dij) δ i = m a x j ( d i j ) 。
论文中给出了一个示例，如图6(a)所示，图中一共有28个样本点，样本点按照密度降序排列。从图中大致可以观察到有两个类簇，剩下的26、27和28号样本点可被视为离群点。在图6(b)中，分别以对判定是否为聚类中的最关键信息ρ和δ为横纵坐标绘制决策图(decision graph)，看到1和10号样本点位于决策图的最右上角。9和10号样本点虽然密度值ρ非常接近，但是δ值却相差很大；被孤立的26、27和28号样本点虽然δ值较大，但是ρ值很小。综上可知，只有ρ值很高并且δ相对较大的样本点才会是聚类中心

　在找出聚类中心后，接下来就是将所有剩下的点划分到比其密度更高且最近的样本点所属的类簇中，当然经过这一步之后暂时会为噪声点也分配到类簇中。

talk is cheap，直接上代码吧…

from sklearn.datasets import make_blobs
import numpy as np   
from sklearn.datasets import make_blobs  
test_data = make_blobs(n_samples = 375,n_features = 2,centers = 3,cluster_std = 1)
test_data = test_data[0]


def distanceNorm(Norm,D_value): #距离规范化
    # initialization 


    # Norm for distance 
    if Norm == '1':  
        counter = np.absolute(D_value)
        counter = np.sum(counter)
    elif Norm == '2':  
        counter = np.power(D_value,2)
        counter = np.sum(counter)
        counter = np.sqrt(counter)
    elif Norm == 'Infinity':  
        counter = np.absolute(D_value)
        counter = np.max(counter)
    else:  
        raise Exception('We will program this later......')

    return counter


def chi(x):  
    if x < 0:  
        return 1#截距范围内的点的个数
    else:  
        return 0


def fit(features,t,distanceMethod = '2'):  
    # initialization 
    labels = list(np.arange(test_data.shape[0]))
    distance = np.zeros((len(labels),len(labels))) #初始化计算距离的矩阵
    distance_sort = list()  #产生一个距离的列表用于距离大小的排列
    density = np.zeros(len(labels)) #初始化密度矩阵，一个点对应了一个密度，所以density
    #的个数也就是labels的个数，labels 的个数就是样本的个数
    distance_higherDensity = np.zeros(len(labels)) #密度比较高的点的初始化矩阵
    ''' １．我们首先要计算的是每个点的密度，要计算每个点的密度，首先就要知道截距dc,然后用所有点到该点 的距离值减去这个截距，如果这个减去截距之后结果小于0的点的个数。 2.每个点的密度求出来之后，求出每个点的delta，为每一个点计算Delta距离，该距离的定义是取比该点 局部密度大的所有点的最小距离，如果这个点已经是局部密度最大的点，那么Delta赋值为别的所有点到它 的最大距离。 3.最后，相对应的聚类中心就是相对Delta值比较大的点，如果好几个相近点的Delta值都比较大， 而且比较接近，那么任何一个都可以当作该类的中心 '''


    # compute distance 产生距离矩阵
    for index_i in range(len(labels)):
        for index_j in range(index_i+1,len(labels)):
            D_value = features[index_i] - features[index_j]
            distance[index_i,index_j] = distanceNorm(distanceMethod,D_value)#距离矩阵中的元素用行纵坐标来索引
            distance_sort.append(distance[index_i,index_j])
    distance += distance.T

    # compute optimal cutoff 　选出最优的截距
    distance_sort = np.array(distance_sort)
    cutoff = distance_sort[int(np.round(len(distance_sort) * t))] #t = 0.02
    ''' 选择的距离是按照全部距离排序（从小到大排序），然后按照前2％的那个位置所在的距离作为最优截距， t为什么选择是0.02，这个在其他的文章中是有证明的，此代码中选出的最优截距是经过四舍五入round之后得到的 '''

    # computer density 产生密度矩阵
    for index_i in range(len(labels)):
        distance_cutoff_i = distance[index_i] - cutoff
        '''这里的distance[index_i]是将distance矩阵按照行的方式一行一行遍历下来，得到的distance_cutoff_i 是一个包含一整行的行向量，'''
        for index_j in range(1,len(labels)):
            '''下面的代码的作用就是将上面得到的行向量distance_cutoff_i遍历一遍 如果distance_cutoff_i的值是小于0的，也就是该点在另外一点的规定的截距范围内，截距范围内点的个数就是该 点的密度，这里面密度的个数就是点的个数，就是样本的个数'''
            density[index_i] += chi(distance_cutoff_i[index_j])

    # search for the max density 
    Max = np.max(density) #　得到密度最大的点
    MaxIndexList = list()
    for index_i in range(len(labels)):
        if density[index_i] == Max:  
            MaxIndexList.extend([index_i])
        '''遍历所有点的密度，如果某点的密度==最大密度，将改点的索引加入到一个列表中， 这样这个列表中就包含了所有密度最大点的索引index'''

    # computer distance_higherDensity 
    Min = 0
    for index_i in range(len(labels)):
        ''' 遍历所有的点，一单遇到之前存在列表里面的索引（密度最大点的索引），就将索引对应的点于所有点距离的行向量拿出来 也就是距离矩阵中密度最大那几行的航中的值都拿出来，取最大值，就是delta的取值。 delta的取值方式有两种： １．为每一个点计算Delta距离，该距离的定义是取比该点局部密度大的所有点的最小距离， ２．如果这个点已经是局部密度最大的点，那么Delta赋值为别的所有点到它的最大距离。 下面的if中求delta的方法是第二种方法 '''
        if index_i in MaxIndexList:  
            distance_higherDensity[index_i] = np.max(distance[index_i])
            continue



        else:
            Min = np.max(distance[index_i])
        for index_j in range(1,len(labels)):
            if density[index_i] < density[index_j] and distance[index_i,index_j] < Min:  
                Min = distance[index_i,index_j]
            else:  
                continue
        distance_higherDensity[index_i] = Min
        ''' else中用到的是求delta的第一种方法： １．先将每个距离行向量中的的距离的最大值拿出来， ２．当处于i点的时候，比较该点和所有点的密度，找出密度比该点大的点，并且距离小于最大距离的点，这个距离赋值给min 但是这时候的min还不是所有密度比该点大并且距离是最小值的min,通过不断的比较将最小值赋给min，得到最小的min 以上就是求delta的方法 '''

    return density,distance_higherDensity,cutoff


######### 这部分是自己的一些想法，用于实现自动寻找聚类中心的 #################
dense,dis_higher,cutoff = fit(test_data,t = 0.02)
gamma = np.sort(dense * dis_higher,axis = 0)
gamma_sorted = gamma[::-1]
P = 2
dn = int(test_data.shape[0] * P / 100)
sum = 0
for i in range(dn + 1):
    if i == 0 or i == dn:
        sj = 0
    else:
        sj = abs(abs(gamma_sorted[i - 1] - gamma_sorted[i]) - abs(gamma_sorted[i] - gamma_sorted[i + 1]))
        sum += sj
theta = 1 / (dn - 2) * sum
print('this is theta:',theta)
arr = []
for i in range(2,dn):
    if abs(abs(gamma_sorted[i - 1] - gamma_sorted[i]) - abs(gamma_sorted[i] - gamma_sorted[i + 1])) > theta:
        print(abs(abs(gamma_sorted[i - 1] - gamma_sorted[i]) - abs(gamma_sorted[i] - gamma_sorted[i + 1])))
        arr.append(abs(abs(gamma_sorted[i - 1] - gamma_sorted[i]) - abs(gamma_sorted[i] - gamma_sorted[i + 1])))
    else:
        print('this is abs:',
            abs(abs(gamma_sorted[i - 1] - gamma_sorted[i]) - abs(gamma_sorted[i] - gamma_sorted[i + 1])))
        arr.append(0)
print('this is arr:',arr)
ap = np.argmax(arr) + 2                                                                                                                                                                                                                                                                                                                                                                  
print(ap)
print('this is gamma_sorted[ap]:',gamma_sorted[ap - 1],gamma_sorted[ap],gamma_sorted[ap + 1])


plt.figure()
plt.title('GAMMA')
plt.xlabel('x')
plt.ylabel('gamma')
plt.scatter(np.arange(dn),gamma_sorted[:dn],s = 20,c = 'black',alpha = 1)
plt.show()
plt.scatter(np.arange(1,gamma_sorted.shape[0] + 1),gamma_sorted,s = 20,c = 'r',alpha = 1)
plt.show()

下面是我做的一些结果图，从结果图可以看出，这个算法还是能寻找到一些数据集的聚类中心，但是对有些数据集的聚类中心的自动寻找还是有点困难。