1 概述

1.1 无监督学习与聚类算法

决策树、随机森林、逻辑回归虽然有着不同的功能，但却都属于“有监督学习”的一部分，即是说，模型在训练的时候，既需要特征矩阵$X$，也需要真实标签$y$。在机器学习中，还有很大一部分算法是属于“无监督学习”，无监督的算法在训练的时候只需要特征矩阵$X$，不需要标签。曾经学过的PCA算法就是无监督学习中的一种。聚类算法也是无监督学习的代表算法之一。
聚类算法又叫做“无监督分类”，其目的是将数据划分成有意义或有用的组（或簇）。这种划分可以基于业务需求或建模需求来完成，也可以单纯地帮助探索数据的自然结构和分布。如在商业中，如果有大量的当前和潜在客户的信息，可以使用聚类将客户划分成若干组，以便进一步分析和开展营销活动，最有名的客户价值判断模型RFM，就常常和聚类分析共同使用。再比如，聚类可以用于降维和矢量量化（vector quantization），可以将高维特征压缩到一列当中，常常用于图像、声音、视频等非结构化数据，可以大幅度压缩数据量。

	聚类	分类
核心	将数据分成多个组，探索每个组的数据是否有联系	从已经分组的数据中去学习，把 新数据放到已经分好的组中去
学习类型	无监督，无需标签进行训练	有监督，需要标签进行训练
典型算法	K-Means、DBSCAN、层次聚类、光谱聚类	决策树、贝叶斯、逻辑回归
算法输出	聚类结果是不确定的，不一定总是能够反映数据的真实分类。同样的聚类，根据不同的业务需求，可能是一个好结果，也可能是一个坏结果。	分类结果是确定的，分类的优劣势客观的，不是根据业务或算法需求决定的。

sklearn_10">1.2 sklearn中的聚类算法

聚类算法在sklearn中有两种表现形式，一种是类，需要实例化，训练并使用接口和属性来调用结果。另一种是函数（function），只需要输入特征矩阵和超参数，即可返回聚类的结果和各种指标。

类	含义	输入:[]内代表选择输入,[]外代表必须输入
cluster.AffinityPropagation	执行亲和传播数据聚类	[damping,…]
cluster.AgglomerativeClustering	凝聚聚类	[…]
cluster.Birch	实现Birch聚类算法	[threshold,branching_factor,…]
cluster.DBSCAN	从矢量数组或距离矩阵执行DBSCAN聚类	[eps,min_samples,metric,…]
cluster.FeatureAgglomeration	凝聚特征	[n_clusters,…]
cluster.KMeans	K均值聚类	[n_clusters,init,n_init,…]
cluster.MiniBatchKMeans	小批量K均值聚类	[n_clusters,init,…]
cluster.MeanShift	使用平坦核函数的平均移位聚类	[bandwidth,seeds,…]
cluster.SpectralClustering	光谱聚类，将聚类应用于规范化拉普拉斯的投影	[n_clusters,…]

函数	含义	输入
cluster.affinity_propagation	执行亲和传播数据聚类	S[,…]
cluster.dascan	从矢量数组或距离矩阵执行DBSCAN聚类	X[,eps,min_samples,…]
cluster.estimate_bandwidth	估计要使用均值平移算法的带宽	X[,quantile,…]
cluster.k_means	K均值聚类	X,n_clusters[,…]
cluster.mean_shift	使用平坦核函数的平均移位聚类	X[,bandwidth,seeds,…]
cluster.spectral_clustering	将聚类应用于规范化拉普拉斯的投影	affinity[,…]
cluster.ward_tree	光谱聚类，将聚类应用于规范化拉普拉斯的投影	X[,connectivity,…]

• 输入数据
  该模块中实现的算法可以采用不同类型的矩阵作为输入。所有方法都接受形状[n_samples,n_features]的标准特征矩阵，这些可以从sklearn.feature_extraction模块中的类中获得。对于亲和力传播、光谱聚类和DBSCAN，还可以输入形状[n_samples,n_samples]的相似性矩阵，可以使用skalearn.metrics.pairwise模块中的函数来获取相似性矩阵。

2 K-Means

作为聚类算法的典型代表，K-Means可以说是最简单的聚类算法。

2.1 KMeans是如何工作的

关键概念：簇和质心
KMeans算法将一组N个样本的特征矩阵X划分为K个无交集的簇，直观上来看簇是一组一组聚集在一起的数据，在一个簇中的数据被认为是同一类。簇就是聚类的结果表现。 簇中所有数据的均值${\mu }_j$通常被称为这个簇的“质心”（centroids）。在一个二维平面中，一簇数据点的质心的横坐标就是这一簇数据点的横坐标的均值，质心的纵坐标就是这一簇数据点的纵坐标的均值。同理可推广至高维空间。

在KMeans算法中，簇的个数K是一个超参数，需要人为输入来确定。KMeans的核心任务就是根据设定的K，找出K个最优的质心，并将离这些质心最近的数据分别分配到这些质心代表的簇中去，具体过程总结如下：

1. 随机抽取K个样本作为最初的质心
2. 开始循环（将每个样本点分配到离他们最近的质心，生成K个簇；对于每个簇，计算所有被分到该簇的样本点的平均值作为新的质心）
3. 当质心的位置不再发生变化，迭代停止，聚类完成
   那什么情况质心的位置会不再变化？当找到一个质心，在每次迭代中被分配到这个质心上的样本都是一致的，即每次新生成的簇都是一致的，所有的样本点都不会再从一个簇转移到另一个簇，质心就不会变化了。

2.2 簇内误差平方和的定义和解惑

聚类算法聚出的类有什么含义？这些类有什么样的性质？可以认为，被分到同一个簇中的数据是有相似性的，而不同簇中的数据是不同的。当聚类完毕之后，要分别研究每个簇中的样本都有什么样的性质，再根据业务需求制定不同的商业或科技策略。和评分卡案例中的“分箱”一样，聚类算法同样追求“簇内差异小，簇间差异大”，而这个“差异”，由样本点到其所在簇的质心的距离来衡量。
对于一个簇来说，所有样本点到质心的距离之和越小，就认为这个簇中的样本越相似，簇内差异就越小。而距离的衡量方法有多种，令$x$表示簇中的一个样本点，$\mu$表示簇中的质心，$n$表示每个样本点中的特征数目，$i$表示组成点$x$的每个特征，则该样本点到质心的距离可以由以下距离来度量：
欧几里得距离：$d(x,\mu )=\sqrt{({\sum }_{i=1}^n(x_i-{\mu }_i{)}^2)}$
曼哈顿距离：$d(x,\mu )={\sum }_{i=1}^n(|x_i-\mu |)$
余弦距离：$cos\theta =\frac{{\sum }_1^n(x_i\ast \mu )}{\sqrt{{\sum }_1^n(x_i{)}^2}\ast \sqrt{{\sum }_1^n(\mu {)}^2}}$
如采用欧几里得距离，则一个簇中所有样本点到质心的距离的平方和为：
$ClusterSumofSquare(CSS)={\sum }_{j=0}^m{\sum }_{i=1}^n(x_i-{\mu }_i{)}^2$
$TotalClusterSumofSquare={\sum }_{l=1}^{k}CS{S}_l$
其中$m$为一个簇中样本的个数，$j$是每个样本的编号，这个公式被称为簇内平方和（cluster Sum of Square），又叫做Inertia。而将一个数据集中的所有簇的簇内平方和相加，就得到了整体平方和（Total Cluster Sum of Square）,又叫做total inertia。Total Inertia越小，代表着每个簇内样本越相似，聚类的效果就越好。因此KMeans追求的是求解能够让Inertia最小化的质心。实际上，在质心不断变化不断迭代的过程中，总体平方和是越来越小的。并且可以使用数学来证明，当整体平方和最小的时候，质心就不会再发生变化了。如此，K-Means的求解过程就变成了一个最优化问题，即需要将某个指标最小化来求解模型中的一部分信息。K-Means中，在一个固定的簇数K下，最小化总体平方和来求解最佳质心，并基于质心的存在去进行聚类，整体距离平方和的最小值可以使用梯度下降来求解。
K-Means有损失函数吗？损失函数本质是用来衡量模型的拟合效果的，只有有着求解参数需求的算法，才会有损失函数。K-Means不求解什么参数，它的模型本质也没有在拟合数据，而是在对数据进行探索。因此K-Means并不存在什么损失函数，Inertia更像是K-Means的模型评估指标。但由于K-Means中的Inertia和逻辑回归中的损失函数确实非常相似，所以，从“求解模型中的某种信息，用于后续模型的使用”这样的功能来看，可以认为Inertia是K-Means中的损失函数，虽然这种说法并不严谨。而对比来看，在决策树中，有衡量分类效果的指标准确度accuracy，不能通过最小化accuracy来求解某个模型中需要的信息。因此决策树、KNN等算法是绝对没有损失函数。
以上是基于欧几里得距离的计算公式得到Inertia，也可以使用其他距离，每个距离都有对应的Inertia。根据经验，总结出不同距离所对应的质心选择方法和Inertia，在K-Means中，只要使用了正确的质心和距离组合，无论使用什么距离，都可以达到不错的聚类效果：

距离度量	质心	Inertia
欧几里得距离	均值	最小化每个样本点到质心的欧式距离之和
曼哈顿距离	中位数	最小化每个样本点到质心的曼哈顿距离之和
余弦距离	均值	最小化每个样本点到质心的余弦距离之和

这些组合都可以由严格的数学证明来推导。在sklearn中，无法选择使用的距离，只能使用欧氏距离。

2.3 KMeans算法的时间复杂度

除了模型本身的效果之外，还使用另一种角度度量算法：算法复杂度。算法的复杂度分为时间复杂度和空间复杂度，时间复杂度是指执行算法所需要的计算工作量，常用大O符号表述；而空间复杂度是指执行这个算法所需要的内存空间。如果一个算法的效果很好，但需要的时间复杂度和空间复杂度都很大，那将会权衡算法的效果和所需的计算成本。
和KNN一样，K-Means算法是一个计算成本很大的算法。接下来，介绍K-Means算法的时间和空间复杂度来加深对K-Means的理解。
K-Means算法的平均复杂度是$O(k\ast n\ast T)$，其中$k$是超参数，所需要输入的簇数，$n$是整个数据集中的样本量，$T$是所需要的迭代次数（相对的，KNN的平均复杂度是$O(n)$）。最坏的情况下，K-Means的复杂度可以写作$O(n^{(k+2)/p})$，其中$n$是整个数据集中的样本量，$p$是特征总数。
在实践中，比起其他聚类算法，k-means算法已经很快了，但它一般找到Inertia的局部最小值。

sklearnclusterKMeans_70">3 sklearn.cluster.KMeans

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 = None, algorithm = ‘auto’)

3.1 重要参数n_clusters

n_clusters是KMeans中的K，表示告诉模型要分几类。这是KMeans中唯一一个必填的参数，默认为8类，但通常聚类结果会是一个小于8的结果。通常，在开始聚类之前，并不知道n_clusters究竟是多少，因此要对其进行探索。

3.1.1 先进行一次聚类看看

当拿到一个数据集，可以通过绘图先观察一下这个数据集的数据分布，一次来为聚类时输入的n_clusters做一个参考。
首先创建一个有标签的数据集。

from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt
#自己创建数据库
x,y = make_blobs(n_samples = 500,n_features = 2,centers = 4,random_state = 1)
fig, ax1 = plt.subplots(1)
ax1.scatter(x[:,0],x[:,1]
           ,marker = 'o'#点的形状
           ,s = 8#点的大小
           )
plt.show()

#如果想要观察这些点的分布
color = ["red","pink","orange","gray"]
fig,ax1 = plt.subplots(1)
for i in range(4):
    ax1.scatter(x[y == i,0],x[y == i,1]
               ,marker = 'o'
               ,s = 8
               ,c = color[i]
               )
plt.show()

基于这个分布，使用KMeans进行聚类。首先，猜一下聚成几类。

from sklearn.cluster import KMeans
n_clusters = 3
cluster = KMeans(n_clusters = n_clusters, random_state = 0).fit(x)
#重要属性labels_，查看聚好的类别，每个样本所对应的类
y_pred = cluster.labels_
y_pred
#KMeans因为并不需要建立模型或者预测结果，因此只需要fit就能够得到聚类结果
#KMeans也有接口predict和fit_predict，表示学习数据X并对X的类进行预测
#但所得到的的结果和不调用predict直接fit之后调用属性labels一模一样
pre = cluster.fit_predict(x)
pre == y_pred
#什么时候需要predict呢？当数据量太大的时候
#其实不必使用所有的数据来寻找质心，少量的数据就可以帮助我们确定质心
#当数据量非常大的时候，可以使用部分数据来确定质心
#剩下的数据的聚类结果，使用predict来调用
cluster_smallsub = KMeans(n_clusters = n_clusters,random_state = 0).fit(x[:200])
y_pred_ = cluster_smallsub.predict(x)
y_pred_ == y_pred
#但这样的结果肯定与直接fit全部数据会不一致。当不要求那么精确，或者数据量实在太大，那么可以使用这样的方法。
#重要属性cluster_centers_，查看质心
centroid = cluster.cluster_centers_
centroid
#结果：array([[-7.09306648, -8.10994454],
#       	 [-1.54234022,  4.43517599],
#       	 [-8.0862351 , -3.5179868 ]])
centroid.shape
#结果：(3, 2)
#重要属性inertia_，查看总距离平方和
inertia = cluster.inertia_
inertia
#结果：1903.4503741659241
color = ["red","pink","orange","gray"]
fig,ax1 = plt.subplots(1)
for i in range(n_clusters):
    ax1.scatter(x[y_pred == i,0],x[y_pred == i,1]
               ,marker = "o"
               ,s = 8
               ,c = color[i]
               )
ax1.scatter(centroid[:,0],centroid[:,1]
           ,marker = "x"
           ,s = 15
           ,c = "black")
plt.show()

#如果把猜测的簇数换成4，Inertia会怎样？
n_clusters = 4
cluster_ = KMeans(n_clusters = n_clusters,random_state = 0).fit(x)
inertia_ = cluster_.inertia_
inertia_
#结果：908.3855684760603
n_clusters = 5
cluster_ = KMeans(n_clusters = n_clusters,random_state = 0).fit(x)
inertia_ = cluster_.inertia_
inertia_
#结果：811.0841324482416
n_clusters = 500
cluster_ = KMeans(n_clusters = n_clusters,random_state = 0).fit(x)
inertia_ = cluster_.inertia_
inertia_
#结果：0.0

3.1.2 聚类算法的模型评估指标

不同于分类模型和回归，聚类算法的模型评估不是一件简单的事。在分类中，有直接结果（标签）的输出，并且分类的结果有正误之分，所以使用预测的准确度、混淆矩阵、ROC曲线等等指标来进行评估。但无论如何评估，都是在评价“模型找到正确答案”的能力。而回归中，由于要拟合数据，有SSE均方误差，有损失函数来衡量模型的拟合程度。但这些衡量指标都不能使用与聚类。

如何衡量聚类算法的效果？
聚类模型的结果不是某种标签输出，并且聚类的结果是不确定的，其优劣由业务需求或者算法需求来决定，并且没有固定的正确答案。

那如何衡量聚类的效果？KMeans的目标是确保“簇内差异小，簇外差异大”，可以通过衡量簇内差异来衡量聚类的效果。Inertia是用距离来衡量簇内差异的指标，因此，是否可以使用Inertia来作为聚类的衡量指标呢？Inertia越小模型越好。
可以，但是这个指标的缺点和极限太大。首先它不是有界的，只知道Inertia是越小越好，是0最好，但是不知道一个较小的Inertia究竟有没有达到模型的极限，能否继续降低。第二，其计算太容易受到特征数目的影响，数据维度很大的时候，Inertia的计算量会陷入维度诅咒之中，计算量会爆炸，不适合用来一次次评估模型。第三，其会受到超参数K的影响，可以发现，随着K越大，Inertia注定会越来越小，但这并不代表模型的效果越来越好了。第四，Inertia对数据的分布有假设，它假设数据满足凸分布（即数据在二维平面图像上看起来是一个凸函数的样子），并且它假设数据是各向同性的（isotropic），即数据的属性在不同方向上代表着相同的含义。但是现实数据往往不是这样。所以使用Inertia作为评估指标，会让聚类算法在一些细长簇、环形簇或者不规则状的流形时表现不佳：

3.1.2.1 当真实标签已知的时候

虽然在聚类中不输入真实标签，但这不代表拥有的数据一定不具有真实标签，或者一定没有任何参考信息。当然，在现实中，拥有真实标签的情况非常少见（几乎是不可能的）。如果拥有真实标签，更倾向于使用分类算法，但也不排除依然使用聚类算法的可能性。如果有样本真实聚类情况的数据，可以根据聚类算法的结果和真实结果来衡量聚类的效果，常用的有以下三种方法：

3.1.2.2 当真实标签未知的时候：轮廓系数

在99%的情况下，是对没有真实标签的数据进行探索，也就是对不知道真正答案的数据进行聚类。这样的聚类是完全依赖于评价簇内的稠密程度（簇内差异小）和簇间的离散程度（簇间差异大）来评估聚类的效果。其中，轮廓系数是最常用的聚类算法的评价指标。它是对每个样本来定义的，它能够同时衡量：
1）样本与其自身所在的簇中的其他样本的相似度$\text{a}$，等于样本与同一个簇中所有其他点之间的平均距离
2）样本与其他簇中的样本相似度$\text{b}$，等于样本与下一个最近的簇中的所有点之间的平均距离
根据聚类的要求“簇内差异小，簇间差异大”，希望b永远大于a，并且大得越多越好。
单个样本的轮廓系数计算公式为：
$s=\frac{b-a}{max(a,b)}$
这个公式可以被解析为：
$s=\{\begin{array}{l}1-a/b，ifa<b\\ 0，ifa=b\\ b/a-1，ifa>b\end{array}$
很容易理解轮廓系数范围是(-1,1)，其中值越接近1表示样本与自己所在簇中的样本很相似，并且与其他簇中的样本不相似，当样本点与簇外的样本更相似的时候，轮廓系数就为负。当轮廓系数为0时，则代表两个簇中的样本相似度一致，两个簇本应该是一个簇。
如果一个簇中的大多数样本具有比较高的轮廓系数，则簇会有较高的总轮廓系数，则整个数据集的平均轮廓系数越高，说明聚类是合适的。如果许多样本点具有低轮廓系数甚至负值，则聚类是不合适的，聚类的超参数K可能设定得太大或者太小。
在sklearn中，使用模块metrics中的类silhouette_score来计算轮廓系数，它返回的是一个数据集中所有样本的轮廓系数的均值。但还有同在metrics模块中的silhouette_sample，它的参数与轮廓系数一致，但返回的是数据集中每个样本的轮廓系数。

from sklearn.metrics import silhouette_score
from sklearn.metrics import silhouette_samples
x.shape
y_pred
silhouette_score(x,y_pred)
#结果：0.5882004012129721
silhouette_score(x,cluster_.labels_)#分4簇的确比分3簇效果要好
#结果：0.6505186632729437
silhouette_score(x,cluster_.labels_)#分5簇
#结果：0.5746932321727458
silhouette_score(x,cluster_.labels_)#分6簇
#结果：0.5150064498560357

轮廓系数有很多优点，它在有限空间中取值，使得对模型的聚类效果有一个“参考”。并且，轮廓系数对数据的分布没有假设，因此在很多数据集上都表现良好。但它在每个簇的分割比较清晰时表现最好。但轮廓系数也有缺陷，它在凸型的类上表现会虚高，如基于密度进行的聚类、或通过DBSCAN获得的聚类结果。如果使用轮廓系数来衡量，则会表现出比真实聚类效果更高的分数。

3.1.2.3 当真实标签未知时的其他评价指标

除了轮廓系数是最常用的，还有卡林斯基-哈拉巴斯指数（Calinski-Harabaz Index，简称CHI，也被称为方差比标准），戴维斯-布尔丁指数（Davies-Bouldin）以及权变矩阵（Contingency Matrix）可以使用。

标签未知时的评估指标
卡林斯基-哈拉巴斯指数：sklearn.metrics.calinski_harabasz_score(x,y_pred)
戴维斯-布尔丁指数：sklearn.metrics.davies_bouldin_score(x,y_pred)
权变矩阵：sklearn.metrics.cluster.contingency_matrix(x,y_pred)

重点了解卡林斯基-哈拉巴斯指数。Calinski-Harabaz指数越高越好。对于有k个簇的聚类而言，Calinski-Harabaz指数s(k)写作如下公式：
$s(k)=\frac{Tr(B_k)}{Tr(W_k)}\ast \frac{N-k}{k-1}$
其中$N$为数据集中的样本量，$k$为簇的个数（即类别的个数），$B_k$是组间离散矩阵，即不同簇之间的协方差矩阵，$W_k$是簇内离散矩阵，即一个簇内数据的协方差矩阵，而tr表示矩阵的迹。在线性代数中，一个n*n矩阵A的主对角线（从左上方至右下方的对角线）上各个元素的总和被称为矩阵A的迹（或迹数），一般记作$tr(A)$。数据之间的离散程度越高，协方差的迹就会越大。组内离散程度低，协方差矩阵的迹就会越大。组内离散程度低，协方差的迹就会越小，$Tr(W_k)$也就越小，同时，组间离散程度大，协方差的迹也会越大，$Tr(B_k)$就越大，因此Calinski-Harabaz指数越高越好。

from sklearn.metrics import calinski_harabasz_score
#注意新版本这里是calinski_harabasz_score，旧版是calinski_harabaz_score
x
y_pred
calinski_harabasz_score(x,y_pred)
#结果：1809.991966958033

虽然Calinski-Harabaz指数没有界，在凸型的数据上的聚类也会表现虚高。但是比起轮廓系数，它有一个巨大的优点，就是计算非常快速。之前使用过魔法命令%%timeit来计算一个命令的运算时间，本文使用时间戳计算运行时间。

from time import time
#time()：记下每一次time()这一行命令时的时间戳
#时间戳是一行数字，用来记录此时此刻的时间
t0 = time()
calinski_harabasz_score(x,y_pred)
time()-t0
#结果：0.0010018348693847656
t0 = time()
silhouette_score(x,y_pred)
time()-t0
#结果：0.00950932502746582
#时间戳可以通过datetime中的函数fromtimestamp转换成真正的时间格式
import datetime
datetime.datetime.fromtimestamp(t0).strftime("%Y-%m-%d %H:%M:%S")
#结果：'2021-12-14 16:17:19'

可以看出，Calinski-Harabaz指数比轮廓系数的计算快了不止一倍。如果是以万计的数据量，轮廓系数会大大拖慢模型的运行速度。

3.1.3 案例：基于轮廓系数来选择n_cluster

通常会绘制轮廓系数分布图和聚类后的数据分布图来选择最佳的n_clusters.

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_samples, silhouette_score
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import pandas as pd
#基于轮廓系数来选择最佳的n_clusters
#需要知道每个聚出来的类的轮廓系数，以及各个类之间的轮廓系数的对比
#需要知道聚类完毕之后图像的分布
#先设定要分成的簇数
n_clusters = 4
#创建一个画布，画布上共有一行两列两个图
fig,(ax1,ax2) = plt.subplots(1,2)
#画布尺寸
#plt.figure(figsize=[18,7])
fig.set_size_inches(18,7)
#第一个图是轮廓系数图像，是由各个簇的轮廓系数组成的横向条形图
#横向条形图的横坐标是轮廓系数取值，走坐标是每个样本，因为轮廓系数是对于每一个样本进行计算的
#首先设定横坐标
#轮廓系数的取值范围在[-1,1]之间，但至少是希望轮廓系数要大于0
#太长的横坐标不利于可视化，所以只设定X轴的取值在[-0.1,1]之间
ax1.set_xlim([-0.1,1])
#接下来设定纵坐标，通常来说，纵坐标是从0开始，最大值取到x.shape[0]的取值
#但是希望每个簇能够排在一起，不同的簇之间能够有一定的空隙，以便看到不同的条形图聚合成的块，理解其实对应了哪一个簇
#因此在设定纵坐标的取值范围的时候，在x.shape[0]的基础上，加上一个距离(n_clusters+1)*10，留作间隔用
ax1.set_ylim([0,x.shape[0]+(n_clusters+1)*10])

#开始建模，调用聚类好的标签
clusterer = KMeans(n_clusters = n_clusters,random_state = 10).fit(x)
cluster_labels = clusterer.labels_
#调用轮廓系数分数，注意，silhouette_score生成的是所有样本点的轮廓系数均值
#两个需要输入的参数是：特征矩阵x和聚类完毕后的标签
silhouette_avg = silhouette_score(x,cluster_labels)
#用print输出一下结果，现在的簇数量下，整体的轮廓系数是多少
print("For n_cluster = ",n_clusters
     ,"The average silhouette_score is :",silhouette_avg)
#调用silhouette_samples，返回每个样本点的轮廓系数，这就是横坐标
sample_silhouette_values = silhouette_samples(x,cluster_labels)
#设定y轴上的初始取值
y_lower = 10

#接下来，对每一个簇进行循环
for i in range(n_clusters):
    #从每个样本的轮廓系数结果中抽取出第i个簇的轮廓系数，并对它进行排序
    ith_cluster_silhouette_values = sample_silhouette_values[cluster_labels == i]
    #注意，.sort()这个命令会直接改掉原数据的顺序
    ith_cluster_silhouette_values.sort()
    #查看这一个簇中究竟有多少个样本
    size_cluster_i = ith_cluster_silhouette_values.shape[0]
    #这一个簇在y轴上的取值，应该是由初始值（y_lower）开始，到初始值加上这个簇中的样本数量结束（y_upper）
    y_upper = y_lower + size_cluster_i
    #在colormap库中，使用小数来调用颜色的函数
    #在nipy_spectral([输入任意小数来代表一个颜色])
    #在这里，希望每个簇的颜色是不同的，需要的颜色种类刚好是循环的个数的种类
    #在这里，只要能够确保，每次循环生成的小数是不同的，可以使用任意方式来获取小数
    #在这里，是用i的浮点数除以n_clusters，在不同的i下，自然生成不同的小数，以确保所有的簇会有不同的颜色
    color = cm.nipy_spectral(float(i)/n_clusters)
    #开始填充子图1中的内容
    #fill_between是让一个范围中的柱状图都统一颜色的函数
    #fill_betweenx的范围是在纵坐标上
    #fill_betweeny的范围是在横坐标上
    #fill_betweenx的参数应该输入（纵坐标的下限，纵坐标的上限，柱状图的颜色）
    ax1.fill_betweenx(np.arange(y_lower,y_upper)
                     ,ith_cluster_silhouette_values
                     ,facecolor = color
                     ,alpha = 0.7
                     )
    #为每个簇的轮廓系数写上簇的编号，并且让簇的编号显示在坐标轴上每个条形图的中间位置
    #text的参数为(要显示编号的位置的横坐标，要显示编号的位置的纵坐标，要显示的编号内容)
    ax1.text(-0.05
            ,y_lower + 0.5 * size_cluster_i
            ,str(i))
    #为下一个簇计算新的y轴上的初始值，是每一次迭代之后，y的上限再加上10
    #以此来保证，不同的簇的图像之间显示有空隙
    y_lower = y_upper + 10
    
#给图1加上标题，横坐标轴，纵坐标轴的标签
ax1.set_title("The silhouette plot for the various clusters.")
ax1.set_xlabel("The silhouette coefficient values")
ax1.set_ylabel("Cluster label")

#把整个数据集上的轮廓系数的均值以虚线的形式放入图中
ax1.axvline(x = silhouette_avg,color = "red",linestyle = "--")
#让y轴不显示任何刻度
ax1.set_yticks([])
#让x轴上的刻度显示为规定的列表
ax1.set_xticks([-0.1,0,0.2,0.4,0.6,0.8,1])

#开始对第二个图进行处理，首先获取新颜色，由于这里没有循环，因此需要一次性生成多个小数来获取多个颜色
colors = cm.nipy_spectral(cluster_labels.astype(float)/n_clusters)
ax2.scatter(x[:,0],x[:,1]
           ,marker = 'o'#点的形状
           ,s = 8#点的大小
           ,c = colors
           )
#把生成的质心放到图像中去
centers = clusterer.cluster_centers_
#Draw white circles at cluster centers
ax2.scatter(centers[:,0],centers[:,1],marker = 'x',c = 'red',alpha = 1,s = 200)
#为图二设置标题，横坐标标题，纵坐标标题
ax2.set_title("The visualization of the clustered data.")
ax2.set_xlabel("Feature space for the 1st feature")
ax2.set_ylabel("Feature space for the 2nd feature")
#为整个图设置标题
plt.suptitle(("Silhouette analysis for KMeans clustering on sample data with n_clusters = %d" % n_clusters)
            ,fontsize = 14, fontweight = 'bold')
plt.show()

将上述这一过程包装成一个循环，可以得到：

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_samples,silhouette_score
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np

for n_clusters in [2,3,4,5,6,7]:
    n_clusters = n_clusters
    fig,(ax1,ax2) = plt.subplots(1,2)
    fig.set_size_inches(18,7)
    ax1.set_xlim([-0.1,1])
    ax1.set_ylim([0,x.shape[0] + (n_clusters + 1) * 10])
    clusterer = KMeans(n_clusters = n_clusters,random_state = 10).fit(x)
    cluster_labels = clusterer.labels_
    silhouette_avg = silhouette_score(x,cluster_labels)
    print("For n_cluster = ",n_clusters
         ,"The average silhouette_score is :",silhouette_avg)
    sample_silhouette_values = silhouette_samples(x,cluster_labels)
    y_lower = 10
    for i in range(n_clusters):
        ith_cluster_silhouette_values = sample_silhouette_values[cluster_labels == i]
        ith_cluster_silhouette_values.sort()
        size_cluster_i = ith_cluster_silhouette_values.shape[0]
        y_upper = y_lower + size_cluster_i
        color = cm.nipy_spectral(float(i) / n_clusters)
        ax1.fill_betweenx(np.arange(y_lower,y_upper)
                         ,ith_cluster_silhouette_values
                         ,facecolor = color
                         ,alpha = 0.7
                         )
        ax1.text(-0.05
                ,y_lower + 0.5 * size_cluster_i
                ,str(i))
        y_lower = y_upper + 10
    
    ax1.set_title("The silhouette plot for the various clusters.")
    ax1.set_xlabel("The silhouette coefficient values")
    ax1.set_ylabel("Cluster label")
    ax1.axvline(x = silhouette_avg,color = "red",linestyle = "--")
    ax1.set_yticks([])
    ax1.set_xticks([-0.1,0,0.2,0.4,0.6,0.8,1])

    colors = cm.nipy_spectral(cluster_labels.astype(float)/n_clusters)
    ax2.scatter(x[:,0],x[:,1]
               ,marker = 'o'
               ,s = 8
               ,c = colors
               )
    centers = clusterer.cluster_centers_
    ax2.scatter(centers[:,0],centers[:,1],marker = 'x',c = 'red',alpha = 1,s = 200)
    ax2.set_title("The visualization of the clustered data.")
    ax2.set_xlabel("Feature space for the 1st feature")
    ax2.set_ylabel("Feature space for the 2nd feature")
    plt.suptitle(("Silhouette analysis for KMeans clustering on sample data with n_clusters = %d" % n_clusters)
                ,fontsize = 14, fontweight = 'bold')
    plt.show()

3.2 重要参数init & random_state & n_init：初始质心怎么放好？

在K-Means中有一个重要的环节，就是放置初始质心。如果有足够的时间，K-Means一定会收敛，但Inertia可能收敛到局部最小值。Inertia是否能够收敛到真正的最小值很大程度上取决于质心的初始化，init就是用来帮助决定初始化方式的参数。
初始质心放置的位置不同，聚类的结果很可能也会不一样，一个好的质心选择可以让K-Means避免更多的计算，让算法收敛稳定且更快。之前讲解初始质心的放置时，使用“随机”的方法在样本点中抽取k个样本作为初始质心，这种方法显然不符合“稳定且更快”的需求。为此，可以使用random_state参数来控制每次生成的初始质心都在相同位置，甚至可以画出学习曲线来确定最优的random_state是哪个整数值。
一个random_state对应一个质心随机初始化的随机数种子。如果不指定随机种子，则sklearn中的K-Means并不会只选择一个随机模式扔出结果，而会在每个随机种子下运行多次，并使用结果最好的一个随机数种子来作为初始质心。可以使用参数n_init来选择每个随机数种子下运行的次数。这个参数不常用到，默认为10次，如果希望运行的结果更加精确，可以增加参数n_init的值来增加每个随机数种子下运行的次数。然而这种方法依然是具有随机性的。
为了优化选择初始质心的方法，2007年Arthur、David and Sergei Vassilvitskii三人发表了论文“k-means++: The advantages of careful seeding”（点击可跳转下载页面），他们开发了“k-means++”初始化方案，使得初始质心（通常）彼此远离，以此来引导出比随机初始化更可靠的结果。在sklearn中，使用参数init = ‘k-means++’来选择使用k-means++作为质心初始化的方案。通常来说，建议保留默认的“k-means++”的方法。

参数	说明
init	可输入“k-means++”、“random”或者一个n维数组。初始化质心的方法，默认为“k-means++”。输入“k-means++”：一种为K均值聚类选择初始聚类中心的聪明方法，以加速收敛。如果输入了n维数组，数组的形状应该是(n_clusters, n_features)并给出初始质心
random_state	控制每次质心随机初始化的随机数种子
n_init	整数，默认为10。使用不同的质心随机初始化的种子来运行k-means算法的次数。最终结果会是基于Inertia来计算的n_init次连续运行后的最佳输出。

plus = KMeans(n_clusters = 10).fit(x)
plus.n_iter_
#结果：15
random = KMeans(n_clusters = 10,init = "random",random_state = 420).fit(x)
random.n_iter_
#结果：19
t0 = time()
plus = KMeans(n_clusters = 10).fit(x)
plus.n_iter_
time() - t0
#结果：0.07305192947387695
t0 = time()
random = KMeans(n_clusters = 10,init = "random",random_state = 420).fit(x)
random.n_iter_
time() - t0
#结果：0.03700757026672363

k-means++的迭代次数更少，但运行所需的时间却更长。

3.3 重要参数max_iter & tol：让迭代停下来

当质心不再移动，K-Means算法就会停下来。但在完全收敛之前，也可以使用最大迭代次数max_iter，或者两次迭代之间Inertia下降的量tol，这两个参数来让迭代提前停下来。有时当n_clusters选择不符合数据的自然分布，或者为了业务需求必须要填入与数据的自然分布不合的n_clusters，提前让迭代停下来反而能够提升模型的表现。

参数	说明
max_iter	整数，默认为300。单次运行的k-means算法的最大迭代次数
tol	浮点数，默认为1e-4。两次迭代间Inertia下降的量，如果两次迭代之间Inertia下降的值小于tol所设定的值，迭代就会停下

random = KMeans(n_clusters = 10,init = "random",max_iter = 10,random_state = 420).fit(x)
y_pred_max10 = random.labels_
silhouette_score(x,y_pred_max10)
#结果：0.39525864440341574
random = KMeans(n_clusters = 10,init = "random",max_iter = 20,random_state = 420).fit(x)
y_pred_max20 = random.labels_
silhouette_score(x,y_pred_max20)
#结果：0.3401504537571701

3.4 重要属性与重要接口

属性	含义
cluster_centers_	收敛到的质心。如果算法在完全收敛之前就已经停下来（受到参数max_iter和tol的控制），所返回的内容将与labels_属性中反映出来的聚类结果不一致
labels_	每个样本点对应的标签
inertia_	每个样本点到距离他们最近的簇心的均方距离，又叫做“簇内平方和”
n_iter_	实际的迭代次数

接口	输入	功能&返回
fit	训练特征矩阵X,[训练用标签， sample_weight]	拟合模型，计算K均值的聚类结果
fit_predict	训练特征矩阵X,[训练用标签， sample_weight]	返回每个样本所对应的簇的索引，计算质心并且为每个样本预测所在的簇的索引，功能相当于先fit再predict
fit_transform	训练特征矩阵X,[训练用标签， sample_weight]	返回新空间中的特征矩阵，进行聚类并且将特征矩阵X转换到簇距离空间当中，功能相当于先fit再transform
get_params	不需要任何输入	获取该类的参数
predict	测试特征矩阵X,[sample_weight]	预测每个测试集X中的样本的所在簇，并返回每个样本所对应的簇的索引。在矢量量化的相关文献中，cluster_centers_被称为代码簿，而predict返回的每个值是代码簿中最接近的代码的索引
score	测试特征矩阵X,[训练用标签,sample_weight]	返回聚类后的Inertia，即簇内平方和的负数。簇内平方和是Kmeans常用的模型评价指标，簇内平方和越小越好，最佳值为0
set_params	需要新设定的参数	为建立好的类重设参数
transform	任意特征矩阵X	将X转换到簇距离空间中，在新空间中，每个维度（即每个坐标轴）是样本点到集群中心的距离。请注意，即使X是稀疏的，变换返回的数组通常也是密集的

3.5 函数cluster.k_means

sklearn.cluster.k_means(X, n_clusters, sample_weight = None, init = ‘k-means++’, precompute_distances = ‘auto’, n_init = 10, max_iter = 300, verbose = False, tol = 0.0001, random_state = None, copy_x = True, n_jobs = None, algorithm = ‘auto’, return_n_iter = False)
函数k_means的用法和类非常相似，不过函数是输入一系列值，而直接一次性返回所有结果。函数k_means会依次返回质心，每个样本对应的簇的标签，inertia以及最佳迭代次数。

from sklearn.cluster import k_means
k_means(x,4,return_n_iter = True)
#参数return_n_iter默认为False，调整为True后就可以看到返回的最佳迭代次数

4 案例：聚类算法用于降维，KMeans的矢量量化应用

K-Means聚类最重要的应用之一是非结构化数据（图像、声音）上的矢量量化（VQ）。非结构化数据往往占有更多的储存空间，文件本身也会比较大，运算非常缓慢。希望能够在保证数据质量的前提下，尽量地缩小非结构化数据的大小，或者简化非结构化数据的结构，而矢量量化就可以帮助实现这个目的。K-Means聚类的本质是一种降维运用，但与之前学过的任何一种降维算法的思路都不相同。特征选择的降维是直接选择对模型贡献大的特征，PCA降维是聚合信息，而矢量量化的降维是在同等样本量上压缩信息的大小，即不改变特征的数目也不改变样本的数目，只改变在这些特征下的样本上的信息量。
用K-Means聚类中获得的质心来代替原有的数据，可以把数据上的信息量压缩到非常小，但又不损失太多信息。以下通过一张图片的矢量量化来看看K-Means如何实现压缩数据大小，却不损失大多信息量。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.metrics import pairwise_distances_argmin
#对两个序列中的点进行距离匹配的函数
from sklearn.datasets import load_sample_image
#导入图片数据所用的库
from sklearn.utils import shuffle
#实例化，导入颐和园的图片
china = load_sample_image("china.jpg")
china
#查看数据类型
china.dtype
#结果：dtype('uint8')
china.shape#长度*宽度*像素（三个数决定的颜色，三原色）
#结果：(427, 640, 3)
china[0][0]
#这里的图像使用RGB编码，每个像素由三个数字来定义
#结果：array([174, 201, 231], dtype=uint8)
#图像中包含的不同颜色的种类
newimage = china.reshape((427*640,3))
import pandas as pd
pd.DataFrame(newimage).drop_duplicates().shape
#现在有9w+种颜色
#结果：(96615, 3)
#图像可视化
plt.figure(figsize = [15,15])
plt.imshow(china)#导入三维数组形成的图片

#查看模块中的另一张图片
flower = load_sample_image("flower.jpg")
plt.figure(figsize = [15,15])
plt.imshow(flower)

经过图像探索，可以了解到图像中有9w多种颜色。可以试试能否使用K-Means将颜色压缩到64种，还不严重损耗图像的质量。为此，使用K-Means来将9w种颜色聚类成64类，然后使用64个簇的质心来替代全部的9w种颜色。质心有着这样的性质：簇中的点都是离质心最近的样本点。为了比较，还要画出随机压缩到64种颜色的矢量量化图像，需要随机选取64个样本点作为随机质心，计算原数据中每个样本到它们的距离来找出离每个样本最近的随机质心，然后用每个样本所对应的随机质心来替换原本的样本。两种状况下，观察图像可视化之后的状况，以查看图片信息的损失。
首先需要把数据处理成sklearn中的K-Means类能够接受的数据类型。

n_clusters = 64
#plt.imshow()在浮点数上表现非常优异，在这里把china中的数据转换为浮点数，压缩到[0,1]之间
china = np.array(china,dtype = np.float64) / china.max()
#把china从图像格式，转换成矩阵格式
w,h,d = original_shape = tuple(china.shape)
w
#结果：427
h
#结果：640
d
#结果：3
assert d == 3
#assert相当于raise error if not，表示为“不为True就报错”
#要求d必须等于3，如果不等于就报错

#展示assert的功能
d_ = 5
assert d_ == 3,"一个格子中的特征数不等于3"

image_array = np.reshape(china,(w * h,d))#reshape是改变结构
image_array
#结果：
#array([[0.68235294, 0.78823529, 0.90588235],
#       [0.68235294, 0.78823529, 0.90588235],
#       [0.68235294, 0.78823529, 0.90588235],
#       ...,
#       [0.16862745, 0.19215686, 0.15294118],
#       [0.05098039, 0.08235294, 0.02352941],
#       [0.05882353, 0.09411765, 0.02745098]])
image_array.shape
#结果：(273280, 3)

#np.reshape(a,newshape,order = 'C'),reshape函数的第一个参数a是要改变结构的对象，第二个参数是要改变成的新结构
#展示np.reshape的效果
a = np.random.random((2,4))
a
#结果：array([[0.21274756, 0.41660165, 0.99374756, 0.76289839],
#            [0.40746363, 0.62590664, 0.11708159, 0.98503634]])
a.reshape((4,2))
#结果：array([[0.21274756, 0.41660165],
#	         [0.99374756, 0.76289839],
#	         [0.40746363, 0.62590664],
#	         [0.11708159, 0.98503634]])
a.reshape((4,2)).shape
#结果：(4, 2)
np.reshape(a,(4,2))
#结果：array([[0.21274756, 0.41660165],
#	         [0.99374756, 0.76289839],
#	         [0.40746363, 0.62590664],
#	         [0.11708159, 0.98503634]])
p.reshape(a,(2,2,2))
#结果：array([[[0.21274756, 0.41660165],
#	          [0.99374756, 0.76289839]],
#	  
#	         [[0.40746363, 0.62590664],
#	          [0.11708159, 0.98503634]]])
#无论有几维，只要维度之间相乘后的总数据量不变，维度可以随意变换
a.shape
#结果：(2,4)

#对数据进行K-Means的矢量量化
#首先，先使用1000个数据来找出质心
image_array_sample = shuffle(image_array,random_state = 0)[:1000]
kmeans = KMeans(n_clusters = n_clusters,random_state = 0).fit(image_array_sample)
kmeans.cluster_centers_.shape
#结果：(64, 3)
#找出质心之后，按照已存在的质心对所有数据进行聚类
labels = kmeans.predict(image_array)
labels.shape
#结果：(273280,)
set(labels)
#使用质心来替换所有样本
image_kmeans = image_array.copy()
image_kmeans#27w个样本点，9w多种不同的颜色（像素点）
#结果：
#array([[0.68235294, 0.78823529, 0.90588235],
#       [0.68235294, 0.78823529, 0.90588235],
#       [0.68235294, 0.78823529, 0.90588235],
#       ...,
#       [0.16862745, 0.19215686, 0.15294118],
#       [0.05098039, 0.08235294, 0.02352941],
#       [0.05882353, 0.09411765, 0.02745098]])
labels#这27w多个样本点所对应的簇的质心的索引
#结果：array([62, 62, 62, ...,  1,  6,  6])
for i in range(w*h):
    image_kmeans[i] = kmeans.cluster_centers_[labels[i]]
#查看生成的新图片信息
image_kmeans
#结果：
#array([[0.73524384, 0.82021116, 0.91925591],
#       [0.73524384, 0.82021116, 0.91925591],
#       [0.73524384, 0.82021116, 0.91925591],
#       ...,
#       [0.15546218, 0.1557423 , 0.12829132],
#       [0.07058824, 0.0754637 , 0.0508744 ],
#       [0.07058824, 0.0754637 , 0.0508744 ]])
pd.DataFrame(image_kmeans).drop_duplicates().shape
#结果：(64, 3)
#恢复图片的结构
image_kmeans = image_kmeans.reshape(w,h,d)
image_kmeans.shape
#结果：(427, 640, 3)

#对数据进行随机的矢量化
centroid_random = shuffle(image_array, random_state = 0)[:n_clusters]#随机选择n_clusters个样本点作为质心
labels_random = pairwise_distances_argmin(centroid_random,image_array,axis = 0)
#函数pairwise_distances_argmin(x1,x2,axis)，其中x1和x2分别是序列
#用来计算x2中的每个样本到x1中的每个样本点的距离，并返回和x2相同形状的x1中对应的最近的样本点的索引
labels_random.shape
#结果：(273280,)
len(set(labels_random))
#结果：64
#使用随机质心来替换所有样本
image_random = image_array.copy()
for i in range(w*h):
    image_random[i] = centroid_random[labels_random[i]]
#恢复图片的结构
image_random = image_random.reshape(w,h,d)
image_random.shape
#结果：(427, 640, 3)

#将原题，按KMeans矢量量化和随机矢量量化的图像绘制出来
plt.figure(figsize = [10,10])
plt.axis("off")
plt.title("original image (96,615 colors)")
plt.imshow(china)

plt.figure(figsize = [10,10])
plt.axis("off")
plt.title("Quantized image (64 colors, K-Means)")
plt.imshow(image_kmeans)

plt.figure(figsize = [10,10])
plt.axis("off")
plt.title("Quantized image (64 colors, Random)")
plt.imshow(image_random)
plt.show()

参考：https://www.bilibili.com/video/BV1WL4y1H7rD?spm_id_from=333.999.0.0