一、K-means聚类介绍

聚类是一种无监督的学习，它将相似的对象归到同一个簇中。它有点像全自动分类 。聚类方法几乎可以应用于所有对象，簇内的对象越相似，聚类的效果越好。之所以称之为K-均值是因为它可以发现k个不同的簇，且每个簇的中心采用簇中所含值的均值计算而成。

假定有一些数据，现在将相似数据归到一起，簇识别会告诉我们这些簇到底都是些什么。聚类与分类的最大不同在于，分类的目标事先已知，而聚类则不一样。因为其产生的结果与分类相同，而只是类别没有预先定义，聚类有时也被称为无监督分类（unsupervisedclassification）。

K-均值算法的工作流程是这样的。首先，随机确定k个初始点作为质心。然后将数据集中的每个点分配到一个簇中，具体来讲，为每个点找距其最近的质心，并将其分配给该质心所对应的簇。这一步完成之后，每个簇的质心更新为该簇所有点的平均值。算法伪代码如下：

上面提到“最近”质心的说法，意味着需要进行某种距离计算。读者可以使用所喜欢的任意距离度量方法。数据集上K-均值算法的性能会受到所选距离计算方法的影响。

二、算法实现

贴上全文全部代码：

from numpy import *
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans


def load_data_set(file_name):
    """
        Function:
            加载数据
        Parameters:
            file_name - 文件名
        Returns:
            data_mat - 数据矩阵
        Modify:
            2018-11-15
    """
    data_mat = []
    fr = open(file_name)
    for line in fr.readlines():
        cur_line = line.strip().split('\t')
        flt_line = list(map(float, cur_line))
        data_mat.append(flt_line)
    return data_mat


def dist_eclud(vec_a, vec_b):
    """
        Function:
            计算两个向量的欧式距离
        Parameters:
            vec_a - 向量a
            vec_b - 向量b
        Returns:
            欧式距离
        Modify:
            2018-11-15
    """
    return sqrt(sum(power(vec_a - vec_b, 2)))


def rand_cent(data_set, k):
    """
        Function:
            构建一个包含k个随机质心的集合
        Parameters:
            data_set - 数据集
            k - 质心个数
        Returns:
            centroids - 随机质心集合
        Modify:
            2018-11-15
    """
    n = shape(data_set)[1]
    # 初始化为一个(k,n)的矩阵
    centroids = mat(zeros((k, n)))
    # 遍历数据集的每一维度，每维先后选出k个值构成质心
    for j in range(n):
        # 得到该列数据的最小值
        min_j = min(data_set[:, j])
        # 得到该列数据的范围(最大值-最小值)
        range_j = float(max(data_set[:, j] - min_j))
        # rand函数根据给定维度生成[0, 1)之间的数据
        # k个质心向量的第j维数据值随机为位于(最小值，最大值)内的某一值
        centroids[:, j] = min_j + range_j * random.rand(k, 1)
    return centroids


def k_means(data_set, k, dist_meas=dist_eclud, create_cent=rand_cent):
    """
        Function:
            K-均值算法
        Parameters:
            data_set - 数据集
            k - 指定的k个类
            dist_meas - 距离计算方法，默认欧氏距离dist_eclud()
            create_cent - 获得k个质心的方法，默认随机获取rand_cent()
        Returns:
            centroids - k个聚类结果
            cluster_assment - 聚类误差
        Modify:
            2018-11-15
    """
    m = shape(data_set)[0]
    # 初始化一个(m,2)的矩阵
    cluster_assment = mat(zeros((m, 2)))
    # 创建初始的k个质心向量
    centroids = create_cent(data_set, k)
    # 聚类结果是否发生变化的布尔类型
    cluster_changed = True
    # 只要聚类结果一直发生变化，就一直执行聚类算法，直至所有数据点聚类结果不变化
    while cluster_changed:
        cluster_changed = False
        # 遍历数据集每一个样本向量
        for i in range(m):
            # 初始化最小距离最正无穷；最小距离对应索引为-1
            min_dist = inf
            min_index = -1
            # 循环k个类的质心
            for j in range(k):
                # 计算数据点到质心的欧氏距离
                dist_j_i = dist_meas(centroids[j, :], data_set[i, :])
                # 如果距离小于当前最小距离
                if dist_j_i < min_dist:
                    # 当前距离定为当前最小距离
                    min_dist = dist_j_i
                    # 最小距离对应索引对应为j(第j个类)
                    min_index = j
            # 当前聚类结果中第i个样本的聚类结果发生变化：布尔类型置为true，继续聚类算法
            if cluster_assment[i, 0] != min_index:
                cluster_changed = True
            # 更新当前变化样本的聚类结果和平方误差
            cluster_assment[i, :] = min_index, min_dist ** 2
        # print(centroids)
        # 遍历每一个质心
        for cent in range(k):
            # 将数据集中所有属于当前质心类的样本通过条件过滤筛选出来
            pts_in_clust = data_set[nonzero(cluster_assment[:, 0].A == cent)[0]]
            # 计算这些数据的均值（axis=0：求列的均值），作为该类质心向量
            centroids[cent, :] = mean(pts_in_clust, axis=0)
    return centroids, cluster_assment


def plot_cluster(data_set, k, centroids, cluster_assment):
    """
        Function:
            绘制聚类结果
        Parameters:
            data_set - 数据集
            k - 指定的k个类
            centroids - 质心集合
            cluster_assment - 聚类误差
        Returns:
            无
        Modify:
            2018-11-16
    """
    figure = plt.figure()
    ax = figure.add_subplot(111)
    ax.scatter(centroids[:, 0].tolist(), centroids[:, 1].tolist(), marker='+', s=60, c='black')
    markers = ['o', 's', 'v', '*'];
    colors = ['blue', 'green', 'yellow', 'red']
    for i in range(k):
        # data_class = data_set[nonzero(cluster_assment[:, 0].A == i)[0]]
        # KMeans().fit().labels_返回的是一维numpy.ndarray，绘制KMeans是用下面这行
        data_class = data_set[cluster_assment == i]
        ax.scatter(data_class[:, 0].tolist(), data_class[:, 1].tolist(), marker=markers[i], c=colors[i], s=20,
                   alpha=0.5)
    plt.show()


def bi_k_means(data_set, k, dist_meas=dist_eclud):
    """
        Function:
            K-均值算法
        Parameters:
            data_set - 数据集
            k - 指定的k个类
            dist_meas - 距离计算方法，默认欧氏距离dist_eclud()
        Returns:
            centroids - k个聚类结果
            cluster_assment - 聚类误差
        Modify:
            2018-11-24
    """
    m = shape(data_set)[0]
    # 将所有的点看成是一个簇
    # cluster_assment存储（所属的中心编号，距中心的距离）的列表
    cluster_assment = mat(zeros((m, 2)))
    # 获取数据集每一列数据的均值，组成一个长为列数的列表
    centroids_0 = mean(data_set, axis=0).tolist()[0]
    # cent_list存储聚类中心
    cent_list = [centroids_0]
    # 计算当前聚为一类时各个数据点距离质心的平方距离
    for j in range(m):
        cluster_assment[j, 1] = dist_meas(mat(centroids_0), data_set[j, :]) ** 2
    # 当簇小于数目k时
    while (len(cent_list) < k):
        lowest_ees = inf
        # 遍历当前每个聚类
        for i in range(len(cent_list)):
            # 通过数组过滤筛选出属于第i类的数据集合
            # nonzero函数是numpy中用于得到数组array中非零元素的位置（数组索引）的函数
            # matrix矩阵名.A矩阵转化为array数组类型
            pts_in_curr_cluster = data_set[nonzero(cluster_assment[:, 0].A == i)[0], :]
            # 对该类利用二分k-均值算法进行划分，返回划分后结果，及误差
            centroid_mat, split_cluster_ass = k_means(pts_in_curr_cluster, 2, dist_meas)
            # 计算该类划分后两个类的误差平方和
            sse_split = sum(split_cluster_ass[:, 1])
            # 计算数据集中不属于该类的数据的误差平方和
            sse_not_split = sum(cluster_assment[nonzero(cluster_assment[:, 0].A != i)[0], 1])
            # 打印这两项误差值
            print('sse_split, and sse_not_split:', sse_split, sse_not_split)
            # 选择使得误差最小的那个簇进行划分
            if (sse_split + sse_not_split) < lowest_ees:
                # 第i类作为本次划分类
                best_cent_to_split = i
                # 第i类划分后得到的两个质心向量
                best_new_cents = centroid_mat
                # 复制第i类中数据点的聚类结果即误差值
                best_clust_ass = split_cluster_ass.copy()
                # 将划分第i类后的总误差作为当前最小误差
                lowest_ees = sse_split + sse_not_split
        # 数组过滤筛选出本次2-均值聚类划分后类编号为1数据点，将这些数据点类编号变为当前类个数+1，作为新的一个聚类
        best_clust_ass[nonzero(best_clust_ass[:, 0].A == 1)[0], 0] = len(cent_list)
        # 将划分数据集中类编号为0的数据点的类编号仍置为被划分的类编号，使类编号连续不出现空缺
        best_clust_ass[nonzero(best_clust_ass[:, 0].A == 0)[0], 0] = best_cent_to_split
        print('the best_cent_to_split is:', best_cent_to_split)
        print('the len of best_clust_ass is:', len(best_clust_ass))
        # 更新质心列表中的变化后的质心向量
        cent_list[best_cent_to_split] = best_new_cents[0, :].tolist()[0]
        # 添加新的类的质心向量
        cent_list.append(best_new_cents[1, :].tolist()[0])
        # 更新cluster_assment列表中参与2-均值聚类数据点变化后的分类编号，及数据该类的误差平方
        cluster_assment[nonzero(cluster_assment[:, 0].A == best_cent_to_split)[0], :] = best_clust_ass
    return mat(cent_list), cluster_assment


def dist_s_l_c(vec_a, vec_b):
    """
        Function:
            计算地球表面两点之间的距离
        Parameters:
            vec_a - 数据集
            vec_b - 指定的k个类
        Returns:
            地球表面两点之间的距离，单位英里
        Modify:
            2018-11-24
    """
    # sin()和cos()以弧度未输入，将float角度数值转为弧度，即*pi/180
    a = sin(vec_a[0, 1] * pi / 180) * sin(vec_b[0, 1] * pi / 180)
    b = cos(vec_a[0, 1] * pi / 180) * cos(vec_b[0, 1] * pi / 180) * cos(pi * (vec_b[0, 0] - vec_a[0, 0]) / 180)
    return arccos(a + b) * 6371.0


def cluster_clubs(num_clust=5):
    """
        Function:
            簇绘图函数
        Parameters:
            num_clust - 指定簇的个数，默认5个
        Returns:
            无
        Modify:
            2018-11-24
    """
    data_list = []
    for line in open('./machinelearninginaction/Ch10/places.txt').readlines():
        line_arr = line.split('\t')
        data_list.append([float(line_arr[4]), float(line_arr[3])])
    data_mat = mat(data_list)
    # 利用2-均值聚类算法进行聚类
    my_centroids, clust_assing = bi_k_means(data_mat, num_clust, dist_meas=dist_s_l_c)
    # 对聚类结果进行绘图
    fig = plt.figure()
    rect = [0.1, 0.1, 0.8, 0.8]
    scatter_markers = ['s', 'o', '^', '8', 'p', 'd', 'v', 'h', '>', '<']
    axprops = dict(xticks=[], yticks=[])
    # 创建一幅图
    ax0 = fig.add_axes(rect, label='ax0', **axprops)
    img_p = plt.imread('./machinelearninginaction/Ch10/Portland.png')
    ax0.imshow(img_p)
    # 创建矩形
    ax1 = fig.add_axes(rect, label='ax1', frameon=False)
    for i in range(num_clust):
        pts_in_cruu_cluster = data_mat[nonzero(clust_assing[:, 0].A == i)[0], :]
        marker_style = scatter_markers[i % len(scatter_markers)]
        ax1.scatter(pts_in_cruu_cluster[:, 0].flatten().A[0], pts_in_cruu_cluster[:, 1].flatten().A[0],
                    marker=marker_style, s=90)
    ax1.scatter(my_centroids[:, 0].flatten().A[0], my_centroids[:, 1].flatten().A[0], marker='+', s=300)
    plt.show()


if __name__ == '__main__':
    # data_mat = mat(load_data_set('./machinelearninginaction/Ch10/testSet.txt'))
    # # 测试k个随机质心集合
    # t = rand_cent(data_mat, 2)
    # # 测试距离计算方法
    # m = dist_eclud(data_mat[0], data_mat[1])
    # print('k个随机质心集合：', t)
    # print('距离计算：', m)

    # data_mat = mat(load_data_set('./machinelearninginaction/Ch10/testSet.txt'))
    # my_centroids, clust_assing = k_means(data_mat, 4)
    # print(my_centroids)
    #
    # plot_cluster(data_mat, 4, my_centroids, clust_assing)

    # data_mat = mat(load_data_set('./machinelearninginaction/Ch10/testSet2.txt'))
    # my_centroids, clust_assing = bi_k_means(data_mat, 3)
    # print(my_centroids)
    # print(clust_assing)
    #
    # plot_cluster(data_mat, 3, my_centroids, clust_assing)

    # cluster_clubs(5)

    data_mat = mat(load_data_set('./machinelearninginaction/Ch10/testSet2.txt'))
    k = 3
    y_pred = KMeans(n_clusters=k).fit(data_mat)
    centroids = y_pred.cluster_centers_
    clust_assing = y_pred.labels_
    print(clust_assing)
    print(centroids)
    plot_cluster(data_mat, k, centroids, clust_assing)

下面用这个数据来测试实现

选取出来的质心和聚类结果：

2.1 使用后处理来提高聚类性能

有时候当我们观察聚类的结果图时，发现聚类的效果没有那么好，如上图所示，K-means算法在k值选取为3时的聚类结果，我们发现，算法能够收敛但效果较差。显然，这种情况的原因是，算法收敛到了局部最小值，而并不是全局最小值，局部最小值显然没有全局最小值的结果好。既然知道了算法已经陷入了局部最小值，如何才能够进一步提升K-means算法的效果呢？

一种用于度量聚类效果的指标是SSE，即误差平方和， 为所有簇中的全部数据点到簇中心的误差距离的平方累加和。SSE的值如果越小，表示数据点越接近于它们的簇中心，即质心，聚类效果也越好。因为，对误差取平方后，就会更加重视那些远离中心的数据点。

显然，一种改善聚类效果的做法就是降低SSE，那么如何在保持簇数目不变的情况下提高簇的质量呢？

一种方法是：我们可以将具有最大SSE值得簇划分为两个簇（因为，SSE最大的簇一般情况下，意味着簇内的数据点距离簇中心较远），具体地，可以将最大簇包含的点过滤出来并在这些点上运行K-means算法，其中k设为2。同时，当把最大的簇(上图中的下半部分)分为两个簇之后，为了保证簇的数目是不变的，我们可以再合并两个簇。

合并两个簇有两种可以量化的办法：合并最近的质心，或者合并两个使得SSE增幅最小的质心。第一种思路通过计算所有质心之间的距离，然后合并距离最近的两个点来实现。第二种方法需要合并两个簇然后计算总SSE值。必须在所有可能的两个簇上重复上述处理过程，直到找到合并最佳的两个簇为止。这样，就可以满足簇的数目不变。

上面，是对已经聚类完成的结果进行改善的方法，在不改变k值的情况下，上述方法能够起到一定的作用，会使得聚类效果得到一定的改善。我们可以用二分k-均值克服算法收敛于局部最小值问题的K-means算法。

2.2 二分 K-均值算法

二分K-均值（bisecting K-means）算法首先将所有点作为一个簇，然后将该簇一分为二。之后选择其中一个簇继续进行划分，选择哪一个簇进行划分取决于对其划分是否可以最大程度降低SSE的值。上述基于SSE的划分过程不断重复，直到得到用户指定的簇数目为止。

当然，也可以选择SSE最大的簇进行划分，知道簇数目达到用户指定的数目为止。下面就来看一下该算法的实际效果。

通过上述算法，之前陷入局部最小值的的这些数据，经过二分K-means算法多次划分后，逐渐收敛到全局最小值，从而达到了令人满意的聚类效果。

三、示例：对地图上的点进行聚类

现在有一个存有70个地址、城市名、对应经纬度的文本数据，而没有这些地点的距离信息。因为在北极每走几米的经度变化可能达到数十度，而沿着赤道附近走相同的距离，带来的经度变化可能是零，所以我们使用球面余弦定理来计算两个经纬度之间的实际距离。想要对这些地点进行聚类，找到每个簇的质心地点，从而可以安排合理的行程，即质心之间选择交通工具抵达，而位于每个质心附近的地点就可以采取步行的方法抵达。显然，K-means算法可以为我们找到一种更加经济而且高效的出行方式。

四、应用scikit-learn构建KMeans聚类器

数据使用上面用过的testSet2.txt数据集。

从图中看出，scikit-learn的KMeans聚类效果和二分 K-均值聚类效果基本一样。

五，小结

聚类是一种无监督的学习方法。聚类区别于分类，即事先不知道要寻找的内容，没有预先设定好的目标变量。

聚类将数据点归到多个簇中，其中相似的数据点归为同一簇，而不相似的点归为不同的簇。相似度的计算方法有很多，具体的应用选择合适的相似度计算方法。

K-means聚类算法，是一种广泛使用的聚类算法，其中k是需要指定的参数，即需要创建的簇的数目，K-means算法中的k个簇的质心可以通过随机的方式获得，但是这些点需要位于数据范围内。在算法中，计算每个点到质心得距离，选择距离最小的质心对应的簇作为该数据点的划分，然后再基于该分配过程后更新簇的质心。重复上述过程，直至各个簇的质心不再变化为止。

K-means算法虽然有效，但是容易受到初始簇质心的情况而影响，有可能陷入局部最优解。为了解决这个问题，可以使用另外一种称为二分K-means的聚类算法。二分K-means算法首先将所有数据点分为一个簇；然后使用K-means（k=2）对其进行划分；下一次迭代时，选择使得SSE下降程度最大的簇进行划分；重复该过程，直至簇的个数达到指定的数目为止。实验表明，二分K-means算法的聚类效果要好于普通的K-means聚类算法。