经典机器学习算法-第十三章K均值聚类

EDUCATION AND TRAINING

一、算法原理

    所谓聚类算法是指将一堆没有标签的数据自动划分成几类的方法，属于无监督学习方法 。监督学习有特征有标签 。无监督学习只有特征没有标签，两者的区别在这里。

1.1算法思想

    K均值聚类又叫做(k-means算法)是属于无监督学习里的一种最基础最常用聚类算法。所谓聚类即人以类聚、物以群分，将样本按照各自的特点分为不同的类别，所谓无监督即事先不知道任何样本属于哪个类别。如下图所示一些样本被分为了绿色，红色，蓝色的三类。聚类的应用非常广泛包括客户群体的划分，推荐系统，文本聚类中，国家电网用户画像，基于用户位置信息的商业选址等。

    那么K均值聚类究竟是如何把这些样本聚成不同的类别的呢？答案是：样本到聚类中心的距离。(这里使用的就是欧式距离，至于其他的距离度量感兴趣的读者可以问一下度娘，小编在以后也会介绍到)。它的基本思想是，通过迭代方式寻找K个簇的一种划分方案，使得聚类结果对应的代价函数最小。代价函数如下所示可以定义为样本距离所属簇中心点的误差平方和（欧氏距离）。

1.2 算法流程

    下面给出K均值算法的具体步骤：

（1）首先确定一个k值，即我们希望将数据集经过聚类得到k个集合。

（2）从数据集中随机选择k个数据点作为质心。

（3）对数据集中每一个点，计算其与每一个质心的距离（如欧式距离），离哪个质心近，就划分到那个质心所属的集合。

（4）把所有数据归好集合后，一共有k个集合。然后重新计算每个集合的质心（求均值）。

（5）如果新计算出来的质心和原来的质心之间的距离小于某一个设置的阈值（表示重新计算的质心的位置变化不大，趋于稳定，或者说收敛，即损失函数J收敛），我们可以认为聚类已经达到期望的结果，算法终止。

（6）如果新质心和原质心距离变化很大，需要迭代3~5步骤。

    在迭代时，分别先固定簇中心，调整样本所属的类别来让J函数减少，再固定类别，调整簇中心使J函数减少。这两个过程可以交替进行。是不是晕了呢，下面的流程图相信会帮助你理解的：

如果不想看上述公式可以看下图的过程：

（1）,首先空间里的一些样本点(a);

（2）随机初始换两个中心(b)，中心为图中的×;

（3）计算当前样本属于的类别(c);

（4）依据每个类别中的样本重新计算聚类中心(d);

（5）迭代此过程(e),(f)。

1.3 优缺点

    说了这么多，那么真正要掌握K-means算法还要知道其优缺点以及如何调整。

    K均值聚类有很多优点，对于大数据集，K均值聚类是高效的，它的计算复杂度O(NKt)接近于线性，其中N是样本数目，K是聚类簇数，t是迭代次数。

    但是它依然存在一些缺点，该算法受初始值和离群点的影响每次结果不稳定，通常的结果不是全局最优解而是局部最优解，无法很好的解决数据簇分布差别比较大的情况(比如一个类的样本是另一个类的100倍),不太适用于离散类，K值依据人为选取缺乏客观性等。

    针对其缺点调优可以考虑一下步骤：

(1)进行数据归一化和离群点的处理(少量的噪声会对均值产生较大的影响)

(2)选择合理的K值，我们可以采用手肘法，如下图所示纵轴为误差平方和，横轴为不同的K值，这里我们选取K=3，因为K<3曲线下降的熟虑大于k>3。

(3)采用核函数，将空间中的样本点映射到更高的维度进行聚类，该方法又称为核K均值聚类，这样可以增加该样本点线性可分的概率。

    当然还有其他的一些好方法比如Gap statistic方法，感兴趣的读者可以自行百度。

    K-means++算法和ISODATA算法就是基于此改进的模型。K-means++算法的思想是在迭代时选取第n+1个聚类中心时，距离当前n个聚类中心越远的点会有更高的概率会被选为第n+1个聚类中心。

    ISODATA算法的思想也很简单，当某个类别的样本数过少时，把该类去除；当属于某个类别的样本数过多、分散程度较大时，把该类分成两个子类别。

1.4 怎么样设置合理的K值

1.4.1 按需选择

    简单地说就是按照建模的需求和目的来选择聚类的个数。比如说，一个游戏公司想把所有玩家做聚类分析，分成顶级、高级、中级、菜鸟四类，那么K=4；如果房地产公司想把当地的商品房分成高中低三档，那么K=3。按需选择虽然合理，但是未必能保证在做K-Means时能够得到清晰的分界线。

1.4.2 观察法

    就是用肉眼看，看这些点大概聚成几堆。这个方法虽然简单，但是同时也模棱两可。

    左上角是原始点。右上角分成了两类。左下角是三类，左下角是四类。至于K到底是选3还是选4，可能每个人都有不同的选择。

    观察法的另一个缺陷就是：原始数据维数要低，一般是两维（平面散点）或者三维（立体散点），否则人类肉眼则无法观察。对于高维数据，我们通常利用PCA降维，然后再进行肉眼观察。

1.4.3 手肘法

    将所有的的点 按照K 个分类 计算距离总和 。总和越大 分类越不合理。手肘法本质上也是一种间接的观察法。我们将得到K个聚类的中心点Mi, i=1,2,⋯,K。以及每个原始点所对应的聚类Ci,i=1,2,⋯,K。我们通常采用所有样本点到它所在的聚类的中心点的距离的和作为模型的度量，记为DK。

    这里距离可以采用欧式距离。

    对于不同的K，最后我们会得到不同的中心点和聚类，所有会有不同的度量。

    我们把上面的例子用不同的K去计算，会得到不同的结果。把K作为横坐标，DK作为纵坐标，我们可以得到下面的折线。

K 为 分类的个数 DK 为 每一分类 到质心的距离之和 再相加。

    很显然K越大，距离和越小。但是我们注意到K=3是一个拐点，就像是我们的肘部一样，K=1到3下降很快，K=3之后趋于平稳。手肘法认为这个拐点就是最佳的K。

手肘法是一个经验方法，而且肉眼观察也因人而异，特别是遇到模棱两可的时候。相比于直接观察法，手肘法的一个优点是，适用于高维的样本数据。有时候人们也会把手肘法用于不同的度量上，如组内方差组间方差比。

1.4.4 Gap Statistic方法

    这里我们要继续使用上面的DK。Gap Statistic的定义为：

    这里E(logDk)指的是logDk的期望。对于一个聚类个数k，首先利用k-means聚类将样本聚成k类，然后计算k类中各类内各点与类中心的距离加和W(ki)，进而计算k类的距离加和W(k)=sum(W(k1),…,W(ki),…,W(kk))；根据原始数据的特点产生B个均匀分布的参考数据集，对于每个数据集都计算W(sk)，计算B个数据集的平均E.W(k)=mean(W(1k),…,W(sk),…,W(bk))；那么对于每个k就有：gap(k)=log(E.W(k))-log(W(k))；然后选取最小的k，使得gap(k)为局部最大值，并且超出了其邻居1个标准差，即gap(k)-gap(k+1)>0.25*sd(W(s(k+1)))

    用上图的例子，我们计算了K=1,2,…9对应的Gap Statisitc.

二、Numpy手写代码

import numpy as np# 定义欧式距离def euclidean_distance(x1, x2): distance = 0 # 距离的平方项再开根号 for i in range(len(x1)): distance += pow((x1[i] - x2[i]), 2) return np.sqrt(distance)
print(euclidean_distance(X[0], X[4]))

# 定义中心初始化函数def centroids_init(k, X):    m, n = X.shape    centroids = np.zeros((k, n))    for i in range(k):        # 每一次循环随机选择一个类别中心        centroid = X[np.random.choice(range(m))]        centroids[i] = centroid    return centroids

# 定义样本的最近质心点所属的类别索引def closest_centroid(sample, centroids): closest_i = 0 closest_dist = float('inf') for i, centroid in enumerate(centroids): # 根据欧式距离判断，选择最小距离的中心点所属类别 distance = euclidean_distance(sample, centroid) if distance < closest_dist: closest_i = i closest_dist = distance return closest_i

# 定义构建类别过程def build_clusters(centroids, k, X):    clusters = [[] for _ in range(k)]    for x_i, x in enumerate(X):        # 将样本划分到最近的类别区域        centroid_i = closest_centroid(x, centroids)        clusters[centroid_i].append(x_i)    return clusters

# 根据上一步聚类结果计算新的中心点def calculate_centroids(clusters, k, X): n = X.shape[1] centroids = np.zeros((k, n)) # 以当前每个类样本的均值为新的中心点 for i, cluster in enumerate(clusters): centroid = np.mean(X[cluster], axis=0) centroids[i] = centroid return centroids

# 获取每个样本所属的聚类类别def get_cluster_labels(clusters, X):    y_pred = np.zeros(X.shape[0])    for cluster_i, cluster in enumerate(clusters):        for X_i in cluster:            y_pred[X_i] = cluster_i    return y_pred

# 根据上述各流程定义kmeans算法流程def kmeans(X, k, max_iterations): # 1.初始化中心点 centroids = centroids_init(k, X) # 遍历迭代求解 for _ in range(max_iterations): # 2.根据当前中心点进行聚类 clusters = build_clusters(centroids, k, X) # 保存当前中心点 prev_centroids = centroids # 3.根据聚类结果计算新的中心点 centroids = calculate_centroids(clusters, k, X) # 4.设定收敛条件为中心点是否发生变化 diff = centroids - prev_centroids if not diff.any(): break # 返回最终的聚类标签 return get_cluster_labels(clusters, X)

# 测试数据X = np.array([[0,2],[0,0],[1,0],[5,0],[5,2]])# 设定聚类类别为2个，最大迭代次数为10次labels = kmeans(X, 2, 10)# 打印每个样本所属的类别标签print(labels)

三、scikit-learn集成方法

class sklearn.cluster.KMeans (n_clusters=8, init='k-means++', n_init=10, max_iter=300, tol=0.0001, precompute_distances='auto', verbose=0, random_state=None, copy_x=True, n_jobs=1, algorithm='auto')

参数：

n_clusters：整形，缺省值=8 【生成的聚类数，即产生的质心（centroids）数。】

max_iter：整形，缺省值=300 执行一次k-means算法所进行的最大迭代数。

n_init：整形，缺省值=10 用不同的质心初始化值运行算法的次数，最终解是在inertia意义下选出的最优结果。

init：有三个可选值：’k-means++’， 'random’，或者传递一个ndarray向量。此参数指定初始化方法，默认值为 'k-means++’。

（１）'k-means++’ 用一种特殊的方法选定初始质心从而能加速迭代过程的收敛（即上文中的k-means++介绍）

（２）'random’ 随机从训练数据中选取初始质心。

（３）如果传递的是一个ndarray，则应该形如 (n_clusters, n_features) 并给出初始质心。

precompute_distances：三个可选值，'auto’，True 或者 False。预计算距离，计算速度更快但占用更多内存。

（１）'auto’：如果 样本数乘以聚类数大于 12million 的话则不预计算距离。This corresponds to about 100MB overhead per job using double precision.

（２）True：总是预先计算距离。

（３）False：永远不预先计算距离。

tol：float形，默认值= 1e-4　与inertia结合来确定收敛条件。

n_jobs：整形数。　指定计算所用的进程数。内部原理是同时进行n_init指定次数的计算。

（１）若值为 -1，则用所有的CPU进行运算。若值为1，则不进行并行运算，这样的话方便调试。

（２）若值小于-1，则用到的CPU数为(n_cpus + 1 + n_jobs)。因此如果 n_jobs值为-2，则用到的CPU数为总CPU数减1。

random_state：整形或 numpy.RandomState 类型，可选用于初始化质心的生成器（generator）。如果值为一个整数，则确定一个seed。此参数默认值为numpy的随机数生成器。

copy_x：布尔型，默认值=True，当我们precomputing distances时，将数据中心化会得到更准确的结果。如果把此参数值设为True，则原始数据不会被改变。如果是False，则会直接在原始数据上做修改并在函数返回值时将其还原。但是在计算过程中由于有对数据均值的加减运算，所以数据返回后，原始数据和计算前可能会有细小差别。

属性：

cluster_centers_：向量，[n_clusters, n_features] (聚类中心的坐标)

Labels_: 每个点的分类

inertia_：float形

每个点到其簇的质心的距离之和。

Notes：

　　这个k-means运用了 Lioyd’s 算法,平均计算复杂度是 O(k*n*T)，其中n是样本量，T是迭代次数。

　　计算复杂读在最坏的情况下为 O(n^(k+2/p))，其中n是样本量，p是特征个数。(D. Arthur and S. Vassilvitskii, 'How slow is the k-means method?’ SoCG2006）

　　在实践中，k-means算法时非常快的，属于可实践的算法中最快的那一类。但是它的解只是由特定初始值所产生的局部解。所以为了让结果更准确真实，在实践中要用不同的初始值重复几次才可以。

Methods：

fit(X[,y]):

　计算k-means聚类。

fit_predictt(X[,y]):

　计算簇质心并给每个样本预测类别。

fit_transform(X[,y])：

计算簇并 transform X to cluster-distance space。

get_params([deep])：

　取得估计器的参数。

predict(X):predict(X)

　给每个样本估计最接近的簇。

score(X[,y]):

　计算聚类误差

set_params(**params):

　为这个估计器手动设定参数。

transform(X[,y]): 将X转换为群集距离空间。

　在新空间中，每个维度都是到集群中心的距离。请注意，即使X是稀疏的，转换返回的数组通常也是密集的。