k-means算法算是经典的机器学习算法，记得应该是大三的课上学过，后面就一直没在接触过该算法对应的问题，现在就来回顾记录一下kmeans算法。

K-means 有一个著名的解释：牧师—村民模型：

有四个牧师去郊区布道，一开始牧师们随意选了几个布道点，并且把这几个布道点的情况公告给了郊区所有的村民，于是每个村民到离自己家最近的布道点去听课。
听课之后，大家觉得距离太远了，于是每个牧师统计了一下自己的课上所有的村民的地址，搬到了所有地址的中心地带，并且在海报上更新了自己的布道点的位置。
牧师每一次移动不可能离所有人都更近，有的人发现A牧师移动以后自己还不如去B牧师处听课更近，于是每个村民又去了离自己最近的布道点……
就这样，牧师每个礼拜更新自己的位置，村民根据自己的情况选择布道点，最终稳定了下来。

我们可以看到该牧师的目的是为了让每个村民到其最近中心点的距离和最小。

所以 K-means 的算法步骤为:

1. 选择初始化的\( \mathrm{k} \)个样本作为初始聚类中心 \(a=a_{1}, a_{2}, \ldots a_{k}\) ；
2. 针对数据集中每个样本 \(x_{i}\) 计算它到 \(\mathrm{k}\) 个聚类中心的距离并将其分到距离最小的聚类中心所对 应的类中；
3. 针对每个类别 \(a_{j}\) ，重新计算它的聚类中心 \(a_{j}=\frac{1}{\left|c_{i}\right|} \sum_{x \in c_{i}} x\) (即属于该类的所有样本的质心）；
4. 重复上面 23 两步操作，直到达到某个中止条件（迭代次数、最小误差变化等）。

K-means算法Python实现代码如下：

# -*- coding:utf-8 -*-
import numpy as np
from matplotlib import pyplot


class K_Means(object):
    # k是分组数；tolerance‘中心点误差’；max_iter是迭代次数
    def __init__(self, k=2, tolerance=0.0001, max_iter=300):
        self.k_ = k
        self.tolerance_ = tolerance
        self.max_iter_ = max_iter

    def fit(self, data):
        self.centers_ = {}
        for i in range(self.k_):
            self.centers_[i] = data[i]

        for i in range(self.max_iter_):
            self.clf_ = {}
            for i in range(self.k_):
                self.clf_[i] = []
            # print("质点:",self.centers_)
            for feature in data:
                # distances = [np.linalg.norm(feature-self.centers[center]) for center in self.centers]
                distances = []
                for center in self.centers_:
                    # 欧拉距离
                    # np.sqrt(np.sum((features-self.centers_[center])**2))
                    distances.append(np.linalg.norm(feature - self.centers_[center]))
                classification = distances.index(min(distances))
                self.clf_[classification].append(feature)

            # print("分组情况:",self.clf_)
            prev_centers = dict(self.centers_)
            for c in self.clf_:
                self.centers_[c] = np.average(self.clf_[c], axis=0)

            # '中心点'是否在误差范围
            optimized = True
            for center in self.centers_:
                org_centers = prev_centers[center]
                cur_centers = self.centers_[center]
                if np.sum((cur_centers - org_centers) / org_centers * 100.0) > self.tolerance_:
                    optimized = False
            if optimized:
                break

    def predict(self, p_data):
        distances = [np.linalg.norm(p_data - self.centers_[center]) for center in self.centers_]
        index = distances.index(min(distances))
        return index


if __name__ == '__main__':
    x = np.array([[1, 2], [1.5, 1.8], [5, 8], [8, 8], [1, 0.6], [9, 11]])
    k_means = K_Means(k=2)
    k_means.fit(x)
    print(k_means.centers_)
    for center in k_means.centers_:
        pyplot.scatter(k_means.centers_[center][0], k_means.centers_[center][1], marker='*', s=150)

    for cat in k_means.clf_:
        for point in k_means.clf_[cat]:
            pyplot.scatter(point[0], point[1], c=('r' if cat == 0 else 'b'))

    predict = [[2, 1], [6, 9]]
    for feature in predict:
        cat = k_means.predict(predict)
        pyplot.scatter(feature[0], feature[1], c=('r' if cat == 0 else 'b'), marker='x')

    pyplot.show()

2.1 优点

• 容易理解，聚类效果不错，虽然是局部最优， 但往往局部最优就够了；
• 处理大数据集的时候，该算法可以保证较好的伸缩性；
• 当簇近似高斯分布的时候，效果非常不错；
• 算法复杂度低。

2.2 缺点

• K 值需要人为设定，不同 K 值得到的结果不一样；
• 对初始的簇中心敏感，不同选取方式会得到不同结果；
• 对异常值敏感；
• 样本只能归为一类，不适合多分类任务；
• 不适合太离散的分类、样本类别不平衡的分类、非凸形状的分类。

针对 K-means 算法的缺点，我们可以有很多种调优方式：如数据预处理（去除异常点），合理选择 K 值，高维映射等。

数据预处理

K-means 的本质是基于欧式距离的数据划分算法，均值和方差大的维度将对数据的聚类产生决定性影响。所以未做归一化处理和统一单位的数据是无法直接参与运算和比较的。常见的数据预处理方式有：数据归一化，数据标准化。

此外，离群点或者噪声数据会对均值产生较大的影响，导致中心偏移，因此我们还需要对数据进行异常点检测。

合理选择 K 值

K 值的选取对 K-means 影响很大，这也是 K-means 最大的缺点，常见的选取 K 值的方法有：手肘法、Gap statistic 方法。

手肘法 (Elbow Method) 本质上也是一种间接的观察法。当我们对数据 \(\left\{x_{1}, \ldots, x_{N}\right\}\) 进行K均值聚 类后，我们将得到 \(K\) 个聚类的中心点 \(\mu_{k}, k=1,2, \ldots, K\) ，以及数据点 \(x_{i}\) 所对应的簇 \(C_{k}, k=1,2, \ldots, K\) ，每个簇中有 \(n_{k}\) 个数据点。我们定义每个簇中的所有点相互之间的距离的和 为

$$D_{k}=\sum_{x_{i} \in C_{k}} \sum_{x_{j} \in C_{k}}\left\|x_{i}-x_{j}\right\|^{2}$$

其中 \(\|\cdot\|\) 为2范数。则对于聚类个数为 \(K\) 时，我们可以定义一个测度量为

$$W_{K}=\sum_{k=1}^{K} \frac{1}{2 n_{k}} D_{k}$$

对于不同的 \(K\) ，经过 \(\mathrm{K}-\mathrm{means}\) 算法后我们会得到不同的中心点和数据点所属的簇，从而得到不同的 度量 \(W_{K}\) 。将聚类个数 \(K\) 作为横坐标， \(W_{K}\) 作为纵坐标，我们可以得到类似下面的图像:

图像中图形很像人的手肘，该方法的命名就是从这而来。从图像中我们可以观察到，K=1 到 4 下降很快，K=4 之后趋于平稳，因此 K=4 是一个拐点，手肘法认为这个拐点就是最佳的聚类个数 K。

 Gap statistic 方法，这个方法出自斯坦福大学的几个学者的论文：Estimating the number of clusters in a data set via the gap statistic

$$
\operatorname{Gap}(K)=\mathrm{E}\left(\log D_{k}\right)-\log D_{k}
$$
其中 \(D_{k}\) 为损失函数，这里\(E\left(\log D_{k}\right)\) 指的是 \(\log D_{k}\) 的期望。这个数值通常通过蒙特卡洛 模拟产生，我们在样本里所在的区域中按照均匀分布随机产生和原始样本数一样多的随机样本，并 对这个随机样本做 K-Means，从而得到一个 \(D_{k}\) 。如此往复多次，通常 20 次，我们可以得到 20 个 \(\log D_{k}\)。对这 20 个数值求平均值，就得到了 \(E\left(\log D_{k}\right)\) 的近似值。最终可以计算 Gap Statisitc。而 Gap statistic 取得最大值所对应的 K 就是最佳的 K。

python实现gap：

import numpy as np


def calculate_Wk(data, centroids, cluster):
    K = centroids.shape[0]
    wk = 0.0
    for k in range(K):
        data_in_cluster = data[cluster == k, :]
        center = centroids[k, :]
        num_points = data_in_cluster.shape[0]
        for i in range(num_points):
            wk = wk + np.linalg.norm(data_in_cluster[i, :]-center, ord=2) ** 2

    return wk


def bounding_box(data):
    dim = data.shape[1]
    boxes = []
    for i in range(dim):
        data_min = np.amin(data[:, i])
        data_max = np.amax(data[:, i])
        boxes.append((data_min, data_max))

    return boxes


def gap_statistic(data, max_K, B, cluster_algorithm):
    num_points, dim = data.shape
    K_range = np.arange(1, max_K, dtype=int)
    num_K = len(K_range)
    boxes = bounding_box(data)
    data_generate = np.zeros((num_points, dim))

    ''' 写法1
    log_Wks = np.zeros(num_K)
    gaps = np.zeros(num_K)
    sks = np.zeros(num_K)
    for ind_K, K in enumerate(K_range):
        cluster_centers, labels, _ = cluster_algorithm(data, K)
        log_Wks[ind_K] = np.log(calculate_Wk(data, cluster_centers, labels))

        # generate B reference data sets
        log_Wkbs = np.zeros(B)
        for b in range(B):
            for i in range(num_points):
                for j in range(dim):
                    data_generate[i][j] = \
                        np.random.uniform(boxes[j][0], boxes[j][1])
            cluster_centers, labels, _ = cluster_algorithm(data_generate, K)
            log_Wkbs[b] = \
                np.log(calculate_Wk(data_generate, cluster_centers, labels))
        gaps[ind_K] = np.mean(log_Wkbs) - log_Wks[ind_K]
        sks[ind_K] = np.std(log_Wkbs) * np.sqrt(1 + 1.0 / B)
    '''

    ''' 写法2
    '''
    log_Wks = np.zeros(num_K)
    for indK, K in enumerate(K_range):
        cluster_centers, labels, _ = cluster_algorithm(data, K)
        log_Wks[indK] = np.log(calculate_Wk(data, cluster_centers, labels))

    gaps = np.zeros(num_K)
    sks = np.zeros(num_K)
    log_Wkbs = np.zeros((B, num_K))

    # generate B reference data sets
    for b in range(B):
        for i in range(num_points):
            for j in range(dim):
                data_generate[i, j] = \
                    np.random.uniform(boxes[j][0], boxes[j][1])
        for indK, K in enumerate(K_range):
            cluster_centers, labels, _ = cluster_algorithm(data_generate, K)
            log_Wkbs[b, indK] = \
                np.log(calculate_Wk(data_generate, cluster_centers, labels))

    for k in range(num_K):
        gaps[k] = np.mean(log_Wkbs[:, k]) - log_Wks[k]
        sks[k] = np.std(log_Wkbs[:, k]) * np.sqrt(1 + 1.0 / B)

    return gaps, sks, log_Wks


采用核函数
基于欧式距离的 K-means 假设了了各个数据簇的数据具有一样的的先验概率并呈现球形分布，但这种分布在实际生活中并不常见。面对非凸的数据分布形状时我们可以引入核函数来优化，这时算法又称为核 K-means 算法，是核聚类方法的一种。核聚类方法的主要思想是通过一个非线性映射，将输入空间中的数据点映射到高位的特征空间中，并在新的特征空间中进行聚类。非线性映射增加了数据点线性可分的概率，从而在经典的聚类算法失效的情况下，通过引入核函数可以达到更为准确的聚类结果。

 K-means++算法 ：改进的算法

原始K-means算法最开始随机选取数据集中K个点作为聚类中心，而K-means++按照如下的思想选取K个聚类中心：假设已经选取了n个初始聚类中心(0<n<K)，则在选取第n+1个聚类中心时：距离当前n个聚类中心越远的点会有更高的概率被选为第n+1个聚类中心。在选取第一个聚类中心(n=1)时同样通过随机的方法。可以说这也符合我们的直觉：聚类中心当然是互相离得越远越好。这个改进虽然直观简单，但是却非常得有效。

1. 随机选取一个中心点 \(a_{1}\) ；
2. 计算数据到之前 \(\mathrm{n}\) 个聚类中心最远的距离 \(D(x)\) ，并以一定概率 \(\frac{D(x)^{2}}{\sum D(x)^{2}}选择新中心点 a_{i} ;\)
3. 重复第二步。

但是这个算法的缺点在于，难以并行化。所以 k-means II 改变取样策略，并非按照 k-means++ 那样每次遍历只取样一个样本，而是每次遍历取样 k 个，重复该取样过程 log(n)次，则得到 klog(n)个样本点组成的集合，然后从这些点中选取 k 个。当然一般也不需要 log(n)次取样，5 次即可。

K-means与ISODATA：ISODATA的全称是迭代自组织数据分析法。在K-means中，K的值需要预先人为地确定，并且在整个算法过程中无法更改。而当遇到高维度、海量的数据集时，人们往往很难准确地估计出K的大小。ISODATA就是针对这个问题进行了改进，它的思想也很直观：当属于某个类别的样本数过少时把这个类别去除，当属于某个类别的样本数过多、分散程度较大时把这个类别分为两个子类别。

K-means与Kernel K-means：传统K-means采用欧式距离进行样本间的相似度度量，显然并不是所有的数据集都适用于这种度量方式。参照支持向量机中核函数的思想，将所有样本映射到另外一个特征空间中再进行聚类，就有可能改善聚类效果。