K-Means

简介

K-means算法是最为经典的基于划分的聚类方法，是十大经典数据挖掘算法之一。K-means算法的基本思想是：以空间中k个点为中心进行聚类，对最靠近他们的对象归类。通过迭代的方法，逐次更新各聚类中心的值，直至得到最好的聚类结果。

算法大致流程

  1. 随机选取k个点作为种子点(这k个点不一定属于数据集)
  2. 分别计算每个数据点到k个种子点的距离，离哪个种子点最近，就属于哪类
  3. 重新计算k个种子点的坐标(简单常用的方法是求坐标值的平均值作为新的坐标值)
  4. 重复2、3步，直到种子点坐标不变或者循环次数完成

完整计算过程

  1. 设置实验数据
    运行之后，效果如下图所示：

    在图中，ABCDE五个点是待分类点，k1、k2是两个种子点。

  2. 计算ABCDE五个点到k1、k2的距离，离哪个点近，就属于哪个点，进行初步分类。
    结果如图：

    A、B属于k1，C、D、E属于k2

  3. 重新计算k1、k2的坐标。这里使用简单的坐标的平均值，使用其他算法也可以(例如以下三个公式)
    a) Minkowski Distance公式——λ可以随意取值，可以是负数，也可以是正数，或是无穷大。
    b) Euclidean Distance公式——也就是第一个公式λ=2的情况
    c) CityBlock Distance公式——也就是第一个公式λ=1的情况

    采用坐标平均值算法的结果如图：

  4. 重复2、3步，直到最终分类完毕。下面是完整的示例代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
import numpy as np
import matplotlib.pyplot as plt

##样本数据(Xi,Yi)，需要转换成数组(列表)形式
Xn=np.array([2,3,1.9,2.5,4])
Yn=np.array([5,4.8,4,1.8,2.2])

#标识符号
sign_n = ['A','B','C','D','E']
sign_k = ['k1','k2']

def start_class(Xk,Yk):
    ##数据点分类
    cls_dict = {}
    ##离哪个分类点最近，属于哪个分类
    for i in range(len(Xn)):
        temp = []
        for j in range(len(Xk)):
            d1 = np.sqrt((Xn[i]-Xk[j])*(Xn[i]-Xk[j])+(Yn[i]-Yk[j])*(Yn[i]-Yk[j]))
            temp.append(d1)
        min_dis=np.min(temp)
        min_inx = temp.index(min_dis)
        cls_dict[sign_n[i]]=sign_k[min_inx]
    #print(cls_dict)
    return cls_dict
    
##重新计算分类的坐标点
def recal_class_point(Xk,Yk,cls_dict):  
    num_k1 = 0  #属于k1的数据点的个数
    num_k2 = 0  #属于k2的数据点的个数
    x1 =0       #属于k1的x坐标和
    y1 =0       #属于k1的y坐标和
    x2 =0       #属于k2的x坐标和
    y2 =0       #属于k2的y坐标和

    ##循环读取已经分类的数据
    for d in cls_dict:
        ##读取d的类别
        kk = cls_dict[d]
        if kk == 'k1':
            #读取d在数据集中的索引
            idx = sign_n.index(d)
            ##累加x值
            x1 += Xn[idx]
            ##累加y值
            y1 += Yn[idx]
            ##累加分类个数
            num_k1 += 1
        else :
            #读取d在数据集中的索引
            idx = sign_n.index(d)
            ##累加x值
            x2 += Xn[idx]
            ##累加y值
            y2 += Yn[idx]
            ##累加分类个数
            num_k2 += 1
    ##求平均值获取新的分类坐标点
    k1_new_x = x1/num_k1 #新的k1的x坐标
    k1_new_y = y1/num_k1 #新的k1的y坐标

    k2_new_x = x2/num_k2 #新的k2的x坐标
    k2_new_y = y2/num_k2 #新的k2的y坐标

    ##新的分类数组
    Xk=np.array([k1_new_x,k2_new_x])
    Yk=np.array([k1_new_y,k2_new_y])
    return Xk,Yk

def draw_point(Xk,Yk,cls_dict):
    #画样本点
    plt.figure(figsize=(5,4)) 
    plt.scatter(Xn,Yn,color="green",label="数据",linewidth=1)
    plt.scatter(Xk,Yk,color="red",label="分类",linewidth=1)
    plt.xticks(range(1,6))
    plt.xlim([1,5])
    plt.ylim([1,6])
    plt.legend()
    for i in range(len(Xn)):
        plt.text(Xn[i],Yn[i],sign_n[i]+":"+cls_dict[sign_n[i]])
        for i in range(len(Xk)):
            plt.text(Xk[i],Yk[i],sign_k[i])
    plt.show()

if __name__ == "__main__":
    ##种子
    Xk=np.array([3.3,3.0])
    Yk=np.array([5.7,3.2])
    for i in range(3):
        cls_dict =start_class(Xk,Yk)
        Xk_new,Yk_new =recal_class_point(Xk,Yk,cls_dict)
        Xk=Xk_new
        Yk=Yk_new
        draw_point(Xk,Yk,cls_dict)

最终分类结果：

由上图可以看出，C点最终是属于k1类，而不是开始的k2.

K-Means的不足

K-Means算法的不足，都是由初始值引起的：

  1. 初始分类数目k值很难估计，不确定应该分成多少类才最合适(ISODATA算法通过类的自动合并和分裂，得到较为合理的类型数目k。这里不讲这个算法)

  2. 不同的随机种子会得到完全不同的结果(K-Means++算法可以用来解决这个问题，其可以有效地选择初始点)

    K-Means++算法

    算法流程如下：

    1. 在数据集中随机挑选1个点作为种子点
      1
      2
      3
      4
      ##随机挑选一个数据点作为种子点
      def select_seed(Xn):
          idx = np.random.choice(range(len(Xn)))
          return idx

  1. 计算剩余数据点到这个点的距离d(x),并且加入到列表

    1
    2
    3
    4
    5
    6
    7
    ##计算数据点到种子点的距离
    def cal_dis(Xn,Yn,idx):
        dis_list = []
        for i in range(len(Xn)):       
            d = np.sqrt((Xn[i]-Xn[idx])**2+(Yn[i]-Yn[idx])**2)
            dis_list.append(d)
        return dis_list

  2. 再取一个随机值。这次的选择思路是：先取一个能落在上步计算的距离列表求和后(sum(dis_list))的随机值rom，然后用rom -= d(x)，直到rom<=0，此时的点就是下一个“种子点”

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    ##随机挑选另外的种子点
    def select_seed_other(Xn,Yn,dis_list):
        d_sum = sum(dis_list)
        rom = d_sum * np.random.random()
        idx = 0
        for i in range(len(Xn)):
            rom -= dis_list[i]
            if rom > 0 :
                continue
            else :
                idx = i
        return idx
  3. 重复第2步和第3步，直到选出k个种子
  4. 进行标准的K-Means算法。下面完整代码
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    133
    134
    135
    136
    137
    138
    139
    140
    141
    142
    143
    144
    145
    146
    147
    148
    149
    150
    151
    152
    153
    154
    155
    156
    157
    158
    159
    160
    161
    162
    163
    164
    165
    166
    167
    168
    169
    170
    171
    172
    173
    174
    175
    176
    import numpy as np
    import matplotlib.pyplot as plt

    ##样本数据(Xi,Yi)，需要转换成数组(列表)形式
    Xn=np.array([2,3,1.9,2.5,4])
    Yn=np.array([5,4.8,4,1.8,2.2])

    #标识符号
    sign_n = ['A','B','C','D','E']
    sign_k = ['k1','k2']

    ##随机挑选一个数据点作为种子点
    def select_seed(Xn):
        idx = np.random.choice(range(len(Xn)))
        return idx
        
    ##计算数据点到种子点的距离
    def cal_dis(Xn,Yn,idx):
        dis_list = []
        for i in range(len(Xn)):       
            d = np.sqrt((Xn[i]-Xn[idx])**2+(Yn[i]-Yn[idx])**2)
            dis_list.append(d)
        return dis_list

    ##随机挑选另外的种子点
    def select_seed_other(Xn,Yn,dis_list):
        d_sum = sum(dis_list)
        rom = d_sum * np.random.random()
        idx = 0
        for i in range(len(Xn)):
            rom -= dis_list[i]
            if rom > 0 :
                continue
            else :
                idx = i
        return idx

    ##选取所有种子点
    def select_seed_all(seed_count):
         ##种子点
        Xk = []  ##种子点x轴列表
        Yk = []  ##种子点y轴列表
        
        idx = 0  ##选取的种子点的索引
        dis_list = [] ##距离列表
        
        
        ##选取种子点
        #因为实验数据少，有一定的几率选到同一个数据，所以加一个判断
        idx_list = []
        flag = True
        for i in range(seed_count):
            if i == 0:
                 idx = select_seed(Xn)
                 dis_list = cal_dis(Xn,Yn,idx)
                 Xk.append(Xn[idx])
                 Yk.append(Yn[idx])
                 idx_list.append(idx)
            else :
                while flag:
                    idx = select_seed_other(Xn,Yn,dis_list)
                    if idx not in idx_list:
                        flag = False
                    else :
                        continue
                dis_list = cal_dis(Xn,Yn,idx)
                Xk.append(Xn[idx])
                Yk.append(Yn[idx])
                idx_list.append(idx)
                    
        ##列表转成数组       
        Xk=np.array(Xk)
        Yk=np.array(Yk)

        return Xk,Yk
        

    def start_class(Xk,Yk):
        ##数据点分类
        cls_dict = {}
        ##离哪个分类点最近，属于哪个分类
        for i in range(len(Xn)):
            temp = []
            for j in range(len(Xk)):
                d1 = np.sqrt((Xn[i]-Xk[j])*(Xn[i]-Xk[j])+(Yn[i]-Yk[j])*(Yn[i]-Yk[j]))
                temp.append(d1)
            min_dis=np.min(temp)
            min_inx = temp.index(min_dis)
            cls_dict[sign_n[i]]=sign_k[min_inx]
        #print(cls_dict)
        return cls_dict
        
    ##重新计算分类的坐标点
    def recal_class_point(Xk,Yk,cls_dict):  
        num_k1 = 0  #属于k1的数据点的个数
        num_k2 = 0  #属于k2的数据点的个数
        x1 =0       #属于k1的x坐标和
        y1 =0       #属于k1的y坐标和
        x2 =0       #属于k2的x坐标和
        y2 =0       #属于k2的y坐标和

        ##循环读取已经分类的数据
        for d in cls_dict:
            ##读取d的类别
            kk = cls_dict[d]
            if kk == 'k1':
                #读取d在数据集中的索引
                idx = sign_n.index(d)
                ##累加x值
                x1 += Xn[idx]
                ##累加y值
                y1 += Yn[idx]
                ##累加分类个数
                num_k1 += 1
            else :
                #读取d在数据集中的索引
                idx = sign_n.index(d)
                ##累加x值
                x2 += Xn[idx]
                ##累加y值
                y2 += Yn[idx]
                ##累加分类个数
                num_k2 += 1
        ##求平均值获取新的分类坐标点
        k1_new_x = x1/num_k1 #新的k1的x坐标
        k1_new_y = y1/num_k1 #新的k1的y坐标

        k2_new_x = x2/num_k2 #新的k2的x坐标
        k2_new_y = y2/num_k2 #新的k2的y坐标

        ##新的分类数组
        Xk=np.array([k1_new_x,k2_new_x])
        Yk=np.array([k1_new_y,k2_new_y])
        return Xk,Yk

    def draw_point(Xk,Yk,cls_dict):
        #画样本点
        plt.figure(figsize=(5,4)) 
        plt.scatter(Xn,Yn,color="green",label="数据",linewidth=1)
        plt.scatter(Xk,Yk,color="red",label="分类",linewidth=1)
        plt.xticks(range(1,6))
        plt.xlim([1,5])
        plt.ylim([1,6])
        plt.legend()
        for i in range(len(Xn)):
            plt.text(Xn[i],Yn[i],sign_n[i]+":"+cls_dict[sign_n[i]])
            for i in range(len(Xk)):
                plt.text(Xk[i],Yk[i],sign_k[i])
        plt.show()

    def draw_point_all_seed(Xk,Yk):
        #画样本点
        plt.figure(figsize=(5,4)) 
        plt.scatter(Xn,Yn,color="green",label="数据",linewidth=1)
        plt.scatter(Xk,Yk,color="red",label="分类",linewidth=1)
        plt.xticks(range(1,6))
        plt.xlim([1,5])
        plt.ylim([1,6])
        plt.legend()
        for i in range(len(Xn)):
            plt.text(Xn[i],Yn[i],sign_n[i])
        plt.show()

    if __name__ == "__main__":

         ##选取2个种子点
         Xk,Yk = select_seed_all(2)
         ##查看种子点
         draw_point_all_seed(Xk,Yk)
         ##循环三次进行分类
         for i in range(3):
            cls_dict =start_class(Xk,Yk)
            Xk_new,Yk_new =recal_class_point(Xk,Yk,cls_dict)
            Xk=Xk_new
            Yk=Yk_new
            draw_point(Xk,Yk,cls_dict)

如图所示，选择了A、E两点作为种子点。

最终的结果。

补充说明：因为数据量太少，在选取所有种子函数的while阶段有可能陷入死循环，所以需要关闭代码重新运行才可以出结果。

sklearn包中的K-Means算法

  1. 函数：sklearn.cluster.KMeans
  2. 主要参数
    a. n_clusters：要进行的分类的个数，即上文中k值，默认是8
    b. max_iter ：最大迭代次数。默认300
    c. min_iter ：最小迭代次数，默认10
    d. init：有三个可选项
     *  'k-means ++'：使用k-means++算法，默认选项
 *  'random':从初始质心数据中随机选择k个观察值
 *  第三个是数组形式的参数
    e. n_jobs: 设置并行量 （-1表示使用所有CPU）
  3. 主要属性：
    a. cluster_centers_ ：集群中心的坐标
    b. labels_ : 每个点的标签
  4. 官网示例： 
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    >>> from sklearn.cluster import KMeans
    >>> import numpy as np
    >>> X = np.array([[1, 2], [1, 4], [1, 0],
    ...               [4, 2], [4, 4], [4, 0]])
    >>> kmeans = KMeans(n_clusters=2, random_state=0).fit(X)
    >>> kmeans.labels_
    array([0, 0, 0, 1, 1, 1], dtype=int32)
    >>> kmeans.predict([[0, 0], [4, 4]])
    array([0, 1], dtype=int32)
    >>> kmeans.cluster_centers_
    array([[ 1.,  2.],
           [ 4.,  2.]])

AP算法

算法简介

AP(Affinity Propagation)通常被翻译为近邻传播算法或者亲和力传播算法，是在2007年的Science杂志上提出的一种新的聚类算法。AP算法的基本思想是将全部数据点都当作潜在的聚类中心(称之为exemplar)，然后数据点两两之间连线构成一个网络(相似度矩阵)，再通过网络中各条边的消息(responsibility和availability)传递计算出各样本的聚类中心。

相关概念(假如有数据点i和数据点j)

  1. 相似度： 点j作为点i的聚类中心的能力，记为S(i,j)。一般使用负的欧式距离，所以S(i,j)越大，表示两个点距离越近，相似度也就越高。使用负的欧式距离，相似度是对称的，如果采用其他算法，相似度可能就不是对称的。
  2. 相似度矩阵：N个点之间两两计算相似度，这些相似度就组成了相似度矩阵。如图1所示的黄色区域，就是一个5*5的相似度矩阵(N=5)
  3. preference：指点i作为聚类中心的参考度(不能为0)，取值为S对角线的值(图1红色标注部分)，此值越大，最为聚类中心的可能性就越大。但是对角线的值为0，所以需要重新设置对角线的值，既可以根据实际情况设置不同的值，也可以设置成同一值。一般设置为S相似度值的中值。(有的说设置成S的最小值产生的聚类最少，但是在下面的算法中设置成中值产生的聚类是最少的)
  4. Responsibility(吸引度):指点k适合作为数据点i的聚类中心的程度，记为r(i,k)。如图2红色箭头所示，表示点i给点k发送信息，是一个点i选点k的过程。
  5. Availability(归属度):指点i选择点k作为其聚类中心的适合程度，记为a(i,k)。如图3红色箭头所示，表示点k给点i发送信息，是一个点k选点i的过程。
  6. exemplar：指的是聚类中心。
  7. r(i, k)加a(i, k)越大,则k点作为聚类中心的可能性就越大,并且i点隶属于以k点为聚类中心的聚类的可能性也越大

数学公式

  1. 吸引度迭代公式: 说明1：$R_{t+1}(i,k)$表示新的$R(i,k)$，$Rt(i,k)$表示旧的$R(i,k)$，也许这样说更容易理解。其中λ是阻尼系数，取值[0.5,1)，用于算法的收敛

说明2：网上还有另外一种数学公式：

sklearn官网的公式是：

我试了这两种公式之后，发现还是公式一的聚类效果最好。同样的数据都采取S的中值作为参考度，我自己写的算法聚类中心是5个，sklearn提供的算法聚类中心是十三个，但是如果把参考度设置为p=-50，则我自己写的算法聚类中心很多，sklearn提供的聚类算法产生标准的3个聚类中心(因为数据是围绕三个中心点产生的)，目前还不清楚这个p=-50是怎么得到的。

  1. 归属度迭代公式 说明：$A_{t+1}(i,k)$表示新的$A(i,k)$，$At(i,k)$表示旧的$A(i,k)$。其中λ是阻尼系数，取值[0.5,1)，用于算法的收敛

详细的算法流程

  1. 设置实验数据。使用sklearn包中提供的函数，随机生成以[1, 1], [-1, -1], [1, -1]三个点为中心的150个数据。 

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    def init_sample():
        ## 生成的测试数据的中心点
        centers = [[1, 1], [-1, -1], [1, -1]]
        ##生成数据
        Xn, labels_true = make_blobs(n_samples=150, centers=centers, cluster_std=0.5,
                                random_state=0)
        #3数据的长度，即：数据点的个数
        dataLen = len(Xn)

        return Xn,dataLen

  2. 计算相似度矩阵，并且设置参考度，这里使用相似度矩阵的中值

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    def cal_simi(Xn):
        ##这个数据集的相似度矩阵，最终是二维数组
        simi = []
        for m in Xn:
            ##每个数字与所有数字的相似度列表，即矩阵中的一行
            temp = []
            for n in Xn:
                ##采用负的欧式距离计算相似度
                s =-np.sqrt((m[0]-n[0])**2 + (m[1]-n[1])**2)
                temp.append(s)
            simi.append(temp)

        ##设置参考度，即对角线的值，一般为最小值或者中值
        #p = np.min(simi)   ##11个中心
        #p = np.max(simi)  ##14个中心
        p = np.median(simi)  ##5个中心
        for i in range(dataLen):
            simi[i][i] = p
        return simi

  3. 计算吸引度矩阵，即R值。
    如果有细心的同学会发现，在上述求R和求A的公式中，求R需要A，求A需要R，所以R或者A不是一开始就可以求解出的，需要先初始化，然后再更新。(我开始就陷入了这个误区，总觉得公式有问题，囧)

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    ##初始化R矩阵、A矩阵
    def init_R(dataLen):
        R = [[0]*dataLen for j in range(dataLen)] 
        return R

    def init_A(dataLen):
        A = [[0]*dataLen for j in range(dataLen)]
        return A

    ##迭代更新R矩阵
    def iter_update_R(dataLen,R,A,simi):
        old_r = 0 ##更新前的某个r值
        lam = 0.5 ##阻尼系数,用于算法收敛
        ##此循环更新R矩阵
        for i in range(dataLen):
            for k in range(dataLen):
                old_r = R[i][k]
                if i != k:
                    max1 = A[i][0] + R[i][0]  ##注意初始值的设置
                    for j in range(dataLen):
                        if j != k:
                            if A[i][j] + R[i][j] > max1 :
                                max1 = A[i][j] + R[i][j]
                    ##更新后的R[i][k]值
                    R[i][k] = simi[i][k] - max1
                    ##带入阻尼系数重新更新
                    R[i][k] = (1-lam)*R[i][k] +lam*old_r
                else:
                    max2 = simi[i][0] ##注意初始值的设置
                    for j in range(dataLen):
                        if j != k:
                            if simi[i][j] > max2:
                                max2 = simi[i][j]
                    ##更新后的R[i][k]值
                    R[i][k] = simi[i][k] - max2
                    ##带入阻尼系数重新更新
                    R[i][k] = (1-lam)*R[i][k] +lam*old_r
        print("max_r:"+str(np.max(R)))
        #print(np.min(R))
        return R

  4. 计算归属度矩阵，即A值

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    ##迭代更新A矩阵
    def iter_update_A(dataLen,R,A):
        old_a = 0 ##更新前的某个a值
        lam = 0.5 ##阻尼系数,用于算法收敛
        ##此循环更新A矩阵
        for i in range(dataLen):
            for k in range(dataLen):
                old_a = A[i][k]
                if i ==k :
                    max3 = R[0][k] ##注意初始值的设置
                    for j in range(dataLen):
                        if j != k:
                            if R[j][k] > 0:
                                max3 += R[j][k]
                            else :
                                max3 += 0
                    A[i][k] = max3
                    ##带入阻尼系数更新A值
                    A[i][k] = (1-lam)*A[i][k] +lam*old_a
                else :
                    max4 = R[0][k] ##注意初始值的设置
                    for j in range(dataLen):
                        ##上图公式中的i!=k 的求和部分
                        if j != k and j != i:
                            if R[j][k] > 0:
                                max4 += R[j][k]
                            else :
                                max4 += 0

                    ##上图公式中的min部分
                    if R[k][k] + max4 > 0:
                        A[i][k] = 0
                    else :
                        A[i][k] = R[k][k] + max4
                        
                    ##带入阻尼系数更新A值
                    A[i][k] = (1-lam)*A[i][k] +lam*old_a
        print("max_a:"+str(np.max(A)))
        #print(np.min(A))
        return A

  5. 迭代更新R值和A值。终止条件是聚类中心在一定程度上不再更新或者达到最大迭代次数

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    ##计算聚类中心
    def cal_cls_center(dataLen,simi,R,A):
        ##进行聚类，不断迭代直到预设的迭代次数或者判断comp_cnt次后聚类中心不再变化
        max_iter = 100    ##最大迭代次数
        curr_iter = 0     ##当前迭代次数
        max_comp = 30     ##最大比较次数
        curr_comp = 0     ##当前比较次数
        class_cen = []    ##聚类中心列表，存储的是数据点在Xn中的索引
        while True:
            ##计算R矩阵
            R = iter_update_R(dataLen,R,A,simi)
            ##计算A矩阵
            A = iter_update_A(dataLen,R,A)
            ##开始计算聚类中心
            for k in range(dataLen):
                if R[k][k] +A[k][k] > 0:
                    if k not in class_cen:
                        class_cen.append(k)
                    else:
                        curr_comp += 1
            curr_iter += 1
            print(curr_iter)
            if curr_iter >= max_iter or curr_comp > max_comp :
                break
        return class_cen

  6. 根据求出的聚类中心，对数据进行分类
    这个步骤产生的是一个归类列表，列表中的每个数字对应着样本数据中对应位置的数据的分类

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    ##根据聚类中心划分数据
       c_list = []
       for m in Xn:
           temp = []
           for j in class_cen:
               n = Xn[j]
               d = -np.sqrt((m[0]-n[0])**2 + (m[1]-n[1])**2)
               temp.append(d)
           ##按照是第几个数字作为聚类中心进行分类标识
           c = class_cen[temp.index(np.max(temp))]
           c_list.append(c)

  7. 完整代码及效果图

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    133
    134
    135
    136
    137
    138
    139
    140
    141
    142
    143
    144
    145
    146
    147
    148
    149
    150
    151
    152
    153
    154
    155
    156
    157
    158
    159
    160
    161
    162
    163
    164
    165
    166
    167
    168
    169
    170
    171
    172
    173
    174
    175
    176
    177
    178
    179
    180
    181
    182
    183
    184
    185
    186
    187
    188
    189
    190
    191
    192
    193
    194
    195
    196
    197
    198
    199
    200
    201
    202
    203
    204
    205
    from sklearn.datasets.samples_generator import make_blobs
    import numpy as np
    import matplotlib.pyplot as plt
    '''
    第一步：生成测试数据
        1.生成实际中心为centers的测试样本300个，
        2.Xn是包含150个(x,y)点的二维数组
        3.labels_true为其对应的真是类别标签
    '''

    def init_sample():
        ## 生成的测试数据的中心点
        centers = [[1, 1], [-1, -1], [1, -1]]
        ##生成数据
        Xn, labels_true = make_blobs(n_samples=150, centers=centers, cluster_std=0.5,
                                random_state=0)
        #3数据的长度，即：数据点的个数
        dataLen = len(Xn)

        return Xn,dataLen

    '''
    第二步：计算相似度矩阵
    '''
    def cal_simi(Xn):
        ##这个数据集的相似度矩阵，最终是二维数组
        simi = []
        for m in Xn:
            ##每个数字与所有数字的相似度列表，即矩阵中的一行
            temp = []
            for n in Xn:
                ##采用负的欧式距离计算相似度
                s =-np.sqrt((m[0]-n[0])**2 + (m[1]-n[1])**2)
                temp.append(s)
            simi.append(temp)

        ##设置参考度，即对角线的值，一般为最小值或者中值
        #p = np.min(simi)   ##11个中心
        #p = np.max(simi)  ##14个中心
        p = np.median(simi)  ##5个中心
        for i in range(dataLen):
            simi[i][i] = p
        return simi

    '''
    第三步：计算吸引度矩阵，即R
           公式1：r(n+1) =s(n)-(s(n)+a(n))-->简化写法，具体参见上图公式
           公式2：r(n+1)=(1-λ)*r(n+1)+λ*r(n)
    '''

    ##初始化R矩阵、A矩阵
    def init_R(dataLen):
        R = [[0]*dataLen for j in range(dataLen)] 
        return R

    def init_A(dataLen):
        A = [[0]*dataLen for j in range(dataLen)]
        return A

    ##迭代更新R矩阵
    def iter_update_R(dataLen,R,A,simi):
        old_r = 0 ##更新前的某个r值
        lam = 0.5 ##阻尼系数,用于算法收敛
        ##此循环更新R矩阵
        for i in range(dataLen):
            for k in range(dataLen):
                old_r = R[i][k]
                if i != k:
                    max1 = A[i][0] + R[i][0]  ##注意初始值的设置
                    for j in range(dataLen):
                        if j != k:
                            if A[i][j] + R[i][j] > max1 :
                                max1 = A[i][j] + R[i][j]
                    ##更新后的R[i][k]值
                    R[i][k] = simi[i][k] - max1
                    ##带入阻尼系数重新更新
                    R[i][k] = (1-lam)*R[i][k] +lam*old_r
                else:
                    max2 = simi[i][0] ##注意初始值的设置
                    for j in range(dataLen):
                        if j != k:
                            if simi[i][j] > max2:
                                max2 = simi[i][j]
                    ##更新后的R[i][k]值
                    R[i][k] = simi[i][k] - max2
                    ##带入阻尼系数重新更新
                    R[i][k] = (1-lam)*R[i][k] +lam*old_r
        print("max_r:"+str(np.max(R)))
        #print(np.min(R))
        return R
    '''
        第四步：计算归属度矩阵，即A
    '''
    ##迭代更新A矩阵
    def iter_update_A(dataLen,R,A):
        old_a = 0 ##更新前的某个a值
        lam = 0.5 ##阻尼系数,用于算法收敛
        ##此循环更新A矩阵
        for i in range(dataLen):
            for k in range(dataLen):
                old_a = A[i][k]
                if i ==k :
                    max3 = R[0][k] ##注意初始值的设置
                    for j in range(dataLen):
                        if j != k:
                            if R[j][k] > 0:
                                max3 += R[j][k]
                            else :
                                max3 += 0
                    A[i][k] = max3
                    ##带入阻尼系数更新A值
                    A[i][k] = (1-lam)*A[i][k] +lam*old_a
                else :
                    max4 = R[0][k] ##注意初始值的设置
                    for j in range(dataLen):
                        ##上图公式中的i!=k 的求和部分
                        if j != k and j != i:
                            if R[j][k] > 0:
                                max4 += R[j][k]
                            else :
                                max4 += 0

                    ##上图公式中的min部分
                    if R[k][k] + max4 > 0:
                        A[i][k] = 0
                    else :
                        A[i][k] = R[k][k] + max4
                        
                    ##带入阻尼系数更新A值
                    A[i][k] = (1-lam)*A[i][k] +lam*old_a
        print("max_a:"+str(np.max(A)))
        #print(np.min(A))
        return A

    '''
       第5步：计算聚类中心
    '''

    ##计算聚类中心
    def cal_cls_center(dataLen,simi,R,A):
        ##进行聚类，不断迭代直到预设的迭代次数或者判断comp_cnt次后聚类中心不再变化
        max_iter = 100    ##最大迭代次数
        curr_iter = 0     ##当前迭代次数
        max_comp = 30     ##最大比较次数
        curr_comp = 0     ##当前比较次数
        class_cen = []    ##聚类中心列表，存储的是数据点在Xn中的索引
        while True:
            ##计算R矩阵
            R = iter_update_R(dataLen,R,A,simi)
            ##计算A矩阵
            A = iter_update_A(dataLen,R,A)
            ##开始计算聚类中心
            for k in range(dataLen):
                if R[k][k] +A[k][k] > 0:
                    if k not in class_cen:
                        class_cen.append(k)
                    else:
                        curr_comp += 1
            curr_iter += 1
            print(curr_iter)
            if curr_iter >= max_iter or curr_comp > max_comp :
                break
        return class_cen
      
       
    if __name__=='__main__':
        ##初始化数据
        Xn,dataLen = init_sample()
        ##初始化R、A矩阵
        R = init_R(dataLen)
        A = init_A(dataLen)
        ##计算相似度
        simi = cal_simi(Xn)   
        ##输出聚类中心
        class_cen = cal_cls_center(dataLen,simi,R,A)
        #for i in class_cen:
        #    print(str(i)+":"+str(Xn[i]))
        #print(class_cen)

        ##根据聚类中心划分数据
        c_list = []
        for m in Xn:
            temp = []
            for j in class_cen:
                n = Xn[j]
                d = -np.sqrt((m[0]-n[0])**2 + (m[1]-n[1])**2)
                temp.append(d)
            ##按照是第几个数字作为聚类中心进行分类标识
            c = class_cen[temp.index(np.max(temp))]
            c_list.append(c)
        ##画图
        colors = ['red','blue','black','green','yellow']
        plt.figure(figsize=(8,6))
        plt.xlim([-3,3])
        plt.ylim([-3,3])
        for i in range(dataLen):
            d1 = Xn[i]
            d2 = Xn[c_list[i]]
            c = class_cen.index(c_list[i])
            plt.plot([d2[0],d1[0]],[d2[1],d1[1]],color=colors[c],linewidth=1)
            #if i == c_list[i] :
            #    plt.scatter(d1[0],d1[1],color=colors[c],linewidth=3)
            #else :
            #    plt.scatter(d1[0],d1[1],color=colors[c],linewidth=1)
        plt.show()

迭代11次出结果：

补充说明：这个算法重点在讲解实现过程，执行效率不是特别高，有优化的空间。以后我会补充进来

sklearn包中的AP算法

  1. 函数：sklearn.cluster.AffinityPropagation
  2. 主要参数：
    a. damping : 阻尼系数，取值[0.5,1)
    b. convergence_iter ：比较多少次聚类中心不变之后停止迭代，默认15
    c. max_iter ：最大迭代次数
    d. preference :参考度
  3. 主要属性
    a. cluster_centers_indices_ : 存放聚类中心的数组
    b. labels_ :存放每个点的分类的数组
    c. n_iter_ : 迭代次数
  4. 示例
    preference(即p值)取不同值时的聚类中心的数目在代码中注明了。
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    from sklearn.cluster import AffinityPropagation
    from sklearn import metrics
    from sklearn.datasets.samples_generator import make_blobs
    import numpy as np


    ## 生成的测试数据的中心点
    centers = [[1, 1], [-1, -1], [1, -1]]
    ##生成数据
    Xn, labels_true = make_blobs(n_samples=150, centers=centers, cluster_std=0.5,
                                random_state=0)



    simi = []
    for m in Xn:
        ##每个数字与所有数字的相似度列表，即矩阵中的一行
        temp = []
        for n in Xn:
             ##采用负的欧式距离计算相似度
            s =-np.sqrt((m[0]-n[0])**2 + (m[1]-n[1])**2)
            temp.append(s)
        simi.append(temp)

    p=-50   ##3个中心
    #p = np.min(simi)  ##9个中心，
    #p = np.median(simi)  ##13个中心    

    ap = AffinityPropagation(damping=0.5,max_iter=500,convergence_iter=30,
                             preference=p).fit(Xn)
    cluster_centers_indices = ap.cluster_centers_indices_

    for idx in cluster_centers_indices:
        print(Xn[idx])

AP算法的优点

  1. 不需要制定最终聚类族的个数
  2. 已有的数据点作为最终的聚类中心，而不是新生成一个族中心。
  3. 模型对数据的初始值不敏感。
  4. 对初始相似度矩阵数据的对称性没有要求。
  5. 相比与k-centers聚类方法，其结果的平方差误差较小。

AP算法的不足

  1. AP算法需要事先计算每对数据对象之间的相似度，如果数据对象太多的话，内存放不下，若存在数据库，频繁访问数据库也需要时间。
  2. AP算法的时间复杂度较高，一次迭代大概O(N3)
  3. 聚类的好坏受到参考度和阻尼系数的影响。

Mean-shift

概述

Mean-shift（即：均值迁移）的基本思想：在数据集中选定一个点，然后以这个点为圆心，r为半径，画一个圆(二维下是圆)，求出这个点到所有点的向量的平均值，而圆心与向量均值的和为新的圆心，然后迭代此过程，直到满足一点的条件结束。(Fukunage在1975年提出)

后来Yizong Cheng 在此基础上加入了 核函数 和 权重系数 ，使得Mean-shift 算法开始流行起来。目前它在聚类、图像平滑、分割、跟踪等方面有着广泛的应用。

图解过程

为了方便大家理解，借用下几张图来说明Mean-shift的基本过程。

由上图可以很容易看到，Mean-shift 算法的核心思想就是不断的寻找新的圆心坐标，直到密度最大的区域。

Mean-shift 算法函数

  1. 核心函数：sklearn.cluster.MeanShift(核函数：RBF核函数)
    由上图可知，圆心(或种子)的确定和半径(或带宽)的选择，是影响算法效率的两个主要因素。所以在sklearn.cluster.MeanShift中重点说明了这两个参数的设定问题。

  2. 主要参数
    a. bandwidth ：半径(或带宽)，float型。如果没有给出，则使用sklearn.cluster.estimate_bandwidth计算出半径(带宽).（可选）
    b. seeds :圆心（或种子），数组类型，即初始化的圆心。（可选）
    c. bin_seeding ：布尔值。如果为真，初始内核位置不是所有点的位置，而是点的离散版本的位置，其中点被分类到其粗糙度对应于带宽的网格上。将此选项设置为True将加速算法，因为较少的种子将被初始化。默认值：False.如果种子参数(seeds)不为None则忽略。

  3. 主要属性
    a. cluster_centers_ : 数组类型。计算出的聚类中心的坐标。
    b. labels_ :数组类型。每个数据点的分类标签。
  4. 算法示例：代码中有详细讲解内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
from sklearn.datasets.samples_generator import make_blobs
from sklearn.cluster import MeanShift, estimate_bandwidth
import numpy as np
import matplotlib.pyplot as plt
from itertools import cycle  ##python自带的迭代器模块

##产生随机数据的中心
centers = [[1, 1], [-1, -1], [1, -1]]
##产生的数据个数
n_samples=10000
##生产数据
X, _ = make_blobs(n_samples=n_samples, centers= centers, cluster_std=0.6, 
                  random_state =0)

##带宽，也就是以某个点为核心时的搜索半径
bandwidth = estimate_bandwidth(X, quantile=0.2, n_samples=500)
##设置均值偏移函数
ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)
##训练数据
ms.fit(X)
##每个点的标签
labels = ms.labels_
print(labels)
##簇中心的点的集合
cluster_centers = ms.cluster_centers_
##总共的标签分类
labels_unique = np.unique(labels)
##聚簇的个数，即分类的个数
n_clusters_ = len(labels_unique)

print("number of estimated clusters : %d" % n_clusters_)


##绘图
plt.figure(1)
plt.clf()

colors = cycle('bgrcmykbgrcmykbgrcmykbgrcmyk')
for k, col in zip(range(n_clusters_), colors):
    ##根据lables中的值是否等于k，重新组成一个True、False的数组
    my_members = labels == k
    cluster_center = cluster_centers[k]
    ##X[my_members, 0] 取出my_members对应位置为True的值的横坐标
    plt.plot(X[my_members, 0], X[my_members, 1], col + '.')
    plt.plot(cluster_center[0], cluster_center[1], 'o', markerfacecolor=col,
             markeredgecolor='k', markersize=14)
plt.title('Estimated number of clusters: %d' % n_clusters_)
plt.show()
  1. 效果图

Spectral Clustering

概述

Spectral Clustering(SC,即谱聚类)，是一种基于图论的聚类方法,它能够识别任意形状的样本空间且收敛于全局最优解，其基本思想是利用样本数据的相似矩阵进行特征分解后得到的特征向量进行聚类.它与样本特征无关而只与样本个数有关。

基本思路：将样本看作顶点,样本间的相似度看作带权的边,从而将聚类问题转为图分割问题:找到一种图分割的方法使得连接不同组的边的权重尽可能低(这意味着组间相似度要尽可能低),组内的边的权重尽可能高(这意味着组内相似度要尽可能高).

图解过程

如上图所示，断开虚线，六个数据被聚成两类。

Spectral Clustering算法函数

核心函数：sklearn.cluster.SpectralClustering

因为是基于图论的算法，所以输入必须是对称矩阵。

主要参数(参数较多，详细参数)

1. n_clusters：聚类的个数。（官方的解释：投影子空间的维度）
2. affinity：核函数，默认是'rbf'，可选："nearest_neighbors"，"precomputed","rbf"或sklearn.metrics.pairwise_kernels支持的其中一个内核之一。
3. gamma :affinity指定的核函数的内核系数，默认1.0

主要属性

labels_ ：每个数据的分类标签

算法示例：代码中有详细讲解内容

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
from sklearn.datasets.samples_generator import make_blobs
from sklearn.cluster import spectral_clustering
import numpy as np
import matplotlib.pyplot as plt
from sklearn import metrics
from itertools import cycle  ##python自带的迭代器模块

##产生随机数据的中心
centers = [[1, 1], [-1, -1], [1, -1]]
##产生的数据个数
n_samples=3000
##生产数据
X, lables_true = make_blobs(n_samples=n_samples, centers= centers, cluster_std=0.6, 
                  random_state =0)

##变换成矩阵，输入必须是对称矩阵
metrics_metrix = (-1 * metrics.pairwise.pairwise_distances(X)).astype(np.int32)
metrics_metrix += -1 * metrics_metrix.min()
##设置谱聚类函数
n_clusters_= 4
lables = spectral_clustering(metrics_metrix,n_clusters=n_clusters_)

##绘图
plt.figure(1)
plt.clf()

colors = cycle('bgrcmykbgrcmykbgrcmykbgrcmyk')
for k, col in zip(range(n_clusters_), colors):
    ##根据lables中的值是否等于k，重新组成一个True、False的数组
    my_members = lables == k
    ##X[my_members, 0] 取出my_members对应位置为True的值的横坐标
    plt.plot(X[my_members, 0], X[my_members, 1], col + '.')
    
plt.title('Estimated number of clusters: %d' % n_clusters_)
plt.show()

效果图

Hierarchical Clustering

概述

Hierarchical Clustering(层次聚类)：就是按照某种方法进行层次分类，直到满足某种条件为止。
主要分成两类：
a. 凝聚：从下到上。首先将每个对象作为一个簇，然后合并这些原子簇为越来越大的簇，直到所有的对象都在一个簇中，或者某个终结条件被满足。
b. 分裂：从上到下。首先将所有对象置于同一个簇中，然后逐渐细分为越来越小的簇，直到每个对象自成一簇，或者达到了某个终止条件。（较少用）

算法步骤

  1. 将每个对象归为一类, 共得到N类, 每类仅包含一个对象. 类与类之间的距离就是它们所包含的对象之间的距离.
  2. 找到最接近的两个类并合并成一类, 于是总的类数少了一个.
  3. 重新计算新的类与所有旧类之间的距离.
  4. 重复第2步和第3步, 直到最后合并成一个类为止(此类包含了N个对象).

图解过程

Hierarchical Clustering算法函数

  1. sklearn.cluster.AgglomerativeClustering

  2. 主要参数(详细参数)

    a. n_clusters：聚类的个数
    b. linkage：指定层次聚类判断相似度的方法，有以下三种：

    • ward：组间距离等于两类对象之间的最小距离。（即single-linkage聚类）
    • average：组间距离等于两组对象之间的平均距离。（average-linkage聚类）
    • complete：组间距离等于两组对象之间的最大距离。（complete-linkage聚类）

  3. 主要属性
    labels_： 每个数据的分类标签

  4. 算法示例：代码中有详细讲解内容
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    from sklearn.datasets.samples_generator import make_blobs
    from sklearn.cluster import AgglomerativeClustering
    import numpy as np
    import matplotlib.pyplot as plt
    from itertools import cycle  ##python自带的迭代器模块

    ##产生随机数据的中心
    centers = [[1, 1], [-1, -1], [1, -1]]
    ##产生的数据个数
    n_samples=3000
    ##生产数据
    X, lables_true = make_blobs(n_samples=n_samples, centers= centers, cluster_std=0.6, 
                      random_state =0)


    ##设置分层聚类函数
    linkages = ['ward', 'average', 'complete']
    n_clusters_ = 3
    ac = AgglomerativeClustering(linkage=linkages[2],n_clusters = n_clusters_)
    ##训练数据
    ac.fit(X)

    ##每个数据的分类
    lables = ac.labels_

    ##绘图
    plt.figure(1)
    plt.clf()

    colors = cycle('bgrcmykbgrcmykbgrcmykbgrcmyk')
    for k, col in zip(range(n_clusters_), colors):
        ##根据lables中的值是否等于k，重新组成一个True、False的数组
        my_members = lables == k
        ##X[my_members, 0] 取出my_members对应位置为True的值的横坐标
        plt.plot(X[my_members, 0], X[my_members, 1], col + '.')
        
    plt.title('Estimated number of clusters: %d' % n_clusters_)
    plt.show()
  5. 效果图：参数linkage的取值依次为：[‘ward’, ‘average’, ‘complete’]

DBSCAN

概述

DBSCAN（Density-Based Spatial Clustering of Applications with Noise,具有噪声的基于密度的聚类方法）是一种基于密度的空间聚类算法。该算法将具有足够密度的区域划分为簇(即要求聚类空间中的一定区域内所包含对象的数目不小于某一给定阈值)，并在具有噪声的空间数据库中发现任意形状的簇，它将簇定义为密度相连的点的最大集合。

算法步骤（大致非详细）

DBSCAN需要二个参数:扫描半径 (eps)和最小包含点数(min_samples)

  1. 遍历所有点，寻找核心点
  2. 连通核心点，并且在此过程中扩展某个分类集合中点的个数

图解过程

在上图中，第一步就是寻找红色的核心点，第二步就是用绿色箭头联通红色点。图中点以绿色线条为中心被分成了两类。没在黑色圆中的点是噪声点。

DBSCAN算法函数

  1. sklearn.cluster.DBSCAN
  2. 主要参数（详细参数）
    a. eps:两个样本之间的最大距离，即扫描半径
    b. min_samples ：作为核心点的话邻域(即以其为圆心，eps为半径的圆，含圆上的点)中的最小样本数(包括点本身)。
  3. 主要属性
    a. core_sample_indices_:核心样本指数。（此参数在代码中有详细的解释）
    b. labels_:数据集中每个点的集合标签给,噪声点标签为-1。
  4. 算法示例：代码中有详细讲解内容
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    from sklearn.datasets.samples_generator import make_blobs
    from sklearn.cluster import DBSCAN
    import numpy as np
    import matplotlib.pyplot as plt
    from itertools import cycle  ##python自带的迭代器模块
    from sklearn.preprocessing import StandardScaler

    ##产生随机数据的中心
    centers = [[1, 1], [-1, -1], [1, -1]]
    ##产生的数据个数
    n_samples=750
    ##生产数据:此实验结果受cluster_std的影响，或者说受eps 和cluster_std差值影响
    X, lables_true = make_blobs(n_samples=n_samples, centers= centers, cluster_std=0.4, 
                      random_state =0)


    ##设置分层聚类函数
    db = DBSCAN(eps=0.3, min_samples=10)
    ##训练数据
    db.fit(X)
    ##初始化一个全是False的bool类型的数组
    core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
    '''
       这里是关键点(针对这行代码：xy = X[class_member_mask & ~core_samples_mask])：
       db.core_sample_indices_  表示的是某个点在寻找核心点集合的过程中暂时被标为噪声点的点(即周围点
       小于min_samples)，并不是最终的噪声点。在对核心点进行联通的过程中，这部分点会被进行重新归类(即标签
       并不会是表示噪声点的-1)，也可也这样理解，这些点不适合做核心点，但是会被包含在某个核心点的范围之内
    '''
    core_samples_mask[db.core_sample_indices_] = True

    ##每个数据的分类
    lables = db.labels_

    ##分类个数：lables中包含-1，表示噪声点
    n_clusters_ =len(np.unique(lables)) - (1 if -1 in lables else 0)

    ##绘图
    unique_labels = set(lables)
    '''
       1)np.linspace 返回[0,1]之间的len(unique_labels) 个数
       2)plt.cm 一个颜色映射模块
       3)生成的每个colors包含4个值，分别是rgba
       4)其实这行代码的意思就是生成4个可以和光谱对应的颜色值
    '''
    colors = plt.cm.Spectral(np.linspace(0, 1, len(unique_labels)))

    plt.figure(1)
    plt.clf()


    for k, col in zip(unique_labels, colors):
        ##-1表示噪声点,这里的k表示黑色
        if k == -1:
            col = 'k'

        ##生成一个True、False数组，lables == k 的设置成True
        class_member_mask = (lables == k)
        
        ##两个数组做&运算，找出即是核心点又等于分类k的值  markeredgecolor='k',
        xy = X[class_member_mask & core_samples_mask]
        plt.plot(xy[:, 0], xy[:, 1], 'o', c=col,markersize=14)
        '''
           1)~优先级最高，按位对core_samples_mask 求反，求出的是噪音点的位置
           2)& 于运算之后，求出虽然刚开始是噪音点的位置，但是重新归类却属于k的点
           3)对核心分类之后进行的扩展
        '''
        xy = X[class_member_mask & ~core_samples_mask]     
        plt.plot(xy[:, 0], xy[:, 1], 'o', c=col,markersize=6)
        
    plt.title('Estimated number of clusters: %d' % n_clusters_)
    plt.show()

    效果图

如果不进行第二步中的扩展，所有的小圆点都应该是噪声点（不符合第一步核心点的要求）

算法优缺点

  1. 优点
    a. 可以发现任意形状的聚类
  2. 缺点
    a. 随着数据量的增加，对I/O、内存的要求也随之增加。
    b. 如果密度分布不均匀，聚类效果较差

Birch

概述

Birch(利用层次方法的平衡迭代规约和聚类)：就是通过聚类特征(CF)形成一个聚类特征树，root层的CF个数就是聚类个数。

相关概念：

聚类特征(CF)：每一个CF是一个三元组,可以用（N，LS，SS）表示.其中N代表了这个CF中拥有的样本点的数量;LS代表了这个CF中拥有的样本点各特征维度的和向量,SS代表了这个CF中拥有的样本点各特征维度的平方和。

如上图所示：N = 5
LS=(3+2+4+4+3,4+6+5+7+8)=(16,30)
SS =(32+22+42+42+32,42+62+52+72+82)=(54,190)

图解过程

对于上图中的CF Tree,限定了B=7,L=5， 也就是说内部节点最多有7个CF(CF90下的圆),而叶子节点最多有5个CF(CF90到CF94)。叶子节点是通过双向链表连通的。

Birch算法函数

  1. sklearn.cluster.Birch
  2. 主要参数（详细参数）
    a. n_clusters ：聚类的目标个数。（可选）
    b. threshold ：扫描半径（个人理解，官方说法比较绕口），设置小了分类就多。
    c. branches_factor：每个节点中CF子集群的最大数量,默认为50。
  3. 主要属性
    labels_ ：每个数据点的分类

算法示例：代码中有详细讲解内容

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets.samples_generator import make_blobs
from sklearn.cluster import Birch

# X为样本特征，Y为样本簇类别， 共1000个样本，每个样本2个特征，共4个簇，簇中心在[-1,-1], [0,0],[1,1], [2,2]
X, y = make_blobs(n_samples=1000, n_features=2, centers=[[-1,-1], [0,0], [1,1], [2,2]], cluster_std=[0.4, 0.3, 0.4, 0.3], 
                  random_state =9)

##设置birch函数
birch = Birch(n_clusters = None)
##训练数据
y_pred = birch.fit_predict(X)
##绘图
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()

效果图：分别为n_clusters = None 和n_clusters = 4

GaussianMixtureModel(补)

概述

正太分布也叫高斯分布，正太分布的概率密度曲线也叫高斯分布概率曲线。

GaussianMixtureModel(混合高斯模型，GMM)。

聚类算法大多数通过相似度来判断，而相似度又大多采用欧式距离长短作为衡量依据。而GMM采用了新的判断依据：概率，即通过属于某一类的概率大小来判断最终的归属类别。

GMM的基本思想就是：任意形状的概率分布都可以用多个高斯分布函数去近似，也就是说GMM就是有多个单高斯密度分布（Gaussian）组成的，每个Gaussian叫一个”Component”，这些”Component”线性加成在一起就组成了 GMM 的概率密度函数，也就是下面的函数。

数学公式

这里不讲公式的具体推导过程，也不实现具体算法。列出来公式只是方便理解下面的函数中为什么需要那些参数。

  1. K：模型的个数，即Component的个数（聚类的个数）
  2. $\pi_k$为第k个高斯的权重
  3. p(x|k)则为第k个高斯概率密度,其均值为μk,方差为σk
  4. 上述参数，除了K是直接给定之外，其他参数都是通过EM算法估算出来的。(有个参数是指定EM算法参数的)

GaussianMixtureModel 算法函数

  1. from sklearn.mixture.GaussianMixture
  2. 主要参数（详细参数）
    a. n_components ：高斯模型的个数，即聚类的目标个数
    b. covariance_type : 通过EM算法估算参数时使用的协方差类型，默认是”full”
     * full：每个模型使用自己的一般协方差矩阵
 * tied：所用模型共享一个一般协方差矩阵
 * diag：每个模型使用自己的对角线协方差矩阵
 * spherical：每个模型使用自己的单一方差

算法示例：代码中有详细讲解内容

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import matplotlib.pyplot as plt
from sklearn.datasets.samples_generator import make_blobs
from sklearn.mixture import GaussianMixture

# X为样本特征，Y为样本簇类别， 共1000个样本，每个样本2个特征，共4个簇，簇中心在[-1,-1], [0,0],[1,1], [2,2]
X, y = make_blobs(n_samples=1000, n_features=2, centers=[[-1,-1], [0,0], [1,1], [2,2]], cluster_std=[0.4, 0.3, 0.4, 0.3], 
                  random_state = 0)

##设置gmm函数
gmm = GaussianMixture(n_components=4, covariance_type='full').fit(X)
##训练数据
y_pred = gmm.predict(X)

##绘图
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()

效果图

跟图15对比可以看出，虽然使用同样的数据，但是不同的算法的聚类效果是不一样的