聚类分析法-K-means聚类

K-means聚类法

1.简介

​ 用层次聚类法聚类时，随着聚类样本对象的增多，计算量会迅速增加，而且聚类结果-谱系图会十分复杂，不便于分析，特别是样品的个数很大(如n>=100)的时，层次聚类法的计算量会非常大，将占据大量的计算机内存空间和较多的计算时间.

​ 为了改进上述缺点， 一个自然的想法是先粗略地分一下类，然后按某种最优原则进行修正，直到将类分得比较合理为止. 基于这种思想就产生了动态聚类法,也称逐步聚类法.

​ 动态距离法适用大型数据.动态聚类法有许多种方法，这里介绍一种比较流行的动态聚类法-K-means法.它是一种快速聚类法，该方法得到的结果简单易懂，对计算机的性能要求不高，因而应用广泛，该方法由Macqueen于1967年提出。

​ 算法的思想是假定样本集中的全体样本可分为C类，并选定C个初始聚类中心，然后根据最小距离原则将每个样品分配到某一类中，之后不断迭代计算各类的聚类中心，并依据新的聚类中心调整聚类状况，直到迭代收敛或聚类中心不再改变.

​ K-means algorithms最后将总样本集G划分为C个子集：G1,G2,···Gc,它们满足下面条件:

​ $$1.G_1\cup G_2\cup \cdot \cdot \cdot \cup G_C=G;$$

​ $$2.G_i\cap G_j=\varnothing (1\le i<j\le C);$$

​ $$3.G_i\ne \varnothing ,G_i\ne G(1\le i\le C)$$.
$$设m_i(i=1,\cdot \cdot \cdot ,C)为C个聚类中心，记$$

​ $$J_e=\sum _{i=1}^C\sum _{\omega \in G_i}{||\omega -m_i||}^2$$

$$使J_e最小的聚类是误差平方和准则下的最优结果.$$
K均值聚类算法描述如下:

$$1.初始化.设总样本集G={\omega }_j,j=1,2\cdot \cdot \cdot ,n是n个样品组成的集合，聚类数为C(2\le C\le n),$$

$$将样品集G任意划分为C类,记为G_1,G_2,\cdot \cdot \cdot ,G_c,计算对应的C个初始聚类中心,记为m_1,m_2,\cdot \cdot \cdot ,m_c，并计算J_e.$$

$$2.G_i=\varnothing (i=1,2,\cdot \cdot \cdot ,C),按最小距离原则将样品{\omega }_j(j=1,2,\cdot \cdot \cdot ,n)进行聚类，即若$$

$$d({\omega }_j,G_k)=mi{n}_{1\le i\le C}d({\omega }_j,m_i),则{\omega }_j\in G_k,G_k=G_k\cup {\omega }_j,j=1,2,\cdot \cdot \cdot ,n.$$

重新计算聚类中心

$$m_i=\frac{1}{n_i}\sum _{{\omega }_j\in G_i}^{{\omega }_j},i=1,2,\cdot \cdot \cdot ,C$$

$$上式中，n_i为当前G_i类中的样本数目，并重新计算J_e.$$
$$3.若连续两次迭代的J_e不变，则算法终止，否则算法转(2).$$

值得注意的是，实际计算(编程)时，可以不计算J_e,只要聚类中心不发生变化，算法即可终止

2.栗子：

已知聚类的指标变量为x1,x2，四个样本点的数据分别为

$${\omega }_1=(1,3),{\omega }_2=(1.5,3.2),{\omega }_3=(1.3,2.8),{\omega }_4=(3,1).$$

要求用k-means algorithm把样本点分成两类.

解：
首先先目测随机分为两类，选定两个聚类中心m1,m2

$$G1类:{\omega }_1值，即m1=(1,3);$$

$$G2类:m2=(\frac{1.5+1.3+3}{3},\frac{3.2+2.8+1}{3})=(1.9333,2.3333)$$

计算每个样本点到G1,G2聚类中心的距离

$$d11=||{\omega }_1-m1||=\sqrt{{(1-1)}^2+{(3-3)}^2}=0$$

$$d12=||{\omega }_1-m2||=1.1470$$

$$d21=||{\omega }_2-m1||=0.5385$$

$$d22=||{\omega }_2-m2||=0.9690$$

$$d31=||{\omega }_3-m1||=0.3606$$

$$d32=||{\omega }_3-m2||=0.7867$$

$$d41=||{\omega }_4-m1||=2.8284$$
$$d42=||{\omega }_4-m2||=1.7075$$

$$经过上面的计算，我们得到新的划分为G1={\omega }_1,{\omega }_2,{\omega }_3,G2={\omega }_4$$

新的聚类中心为：

$$G1类:m1=(\frac{1+1.5+1.3/3}{3},\frac{3+3.2+2.8}{3})=(1.2667,3.0);$$

$$G2类:{\omega }_4值，即m2=(3,1)$$

再次重新计算每个样本点到G1,G2聚类中心的距离

d11,d12,d21,d22,d31,d32,d41,d42

得到的新的划分和我们的上一个最新的假设一样，所以到此聚类结束

当然我们也可以设计程序(利用sklearn库的cluster模块的KMeans函数)进行计算

import numpy as np
from sklearn.cluster import KMeans
a=np.array([[1,3],[1.5,3.2],[1.3,2.8],[3,1]])
md=KMeans(n_clusters=2)	#构建模型 选择参数-构造两个聚类中心的模型
md.fit(a) #求解模型
labels=1+md.labels_ #提取聚类标签 #因为聚类标签一般是从0开始的 
centers=md.cluster_centers_ #提取聚类中心，每一行是一个聚类中心
print(labels,'\n',centers)

3.k-means算法思想揭秘-期望最大化

期望最大化(expectation-maximization,E-M)是一种非常强大的算法，在数据科学中应用广泛.

而K-means就是其一个较为简单且易于理解的应用.

简单的说，期望最大化有以下几个步骤：

1. 猜测一些簇中心点（后面会详细介绍如何猜测）

2. 重复直至收敛

   1. 期望步骤(E-step):将点分配至离其最近的簇中心点

   2. 最大化步骤(M-step):将簇中心点设置为所有点坐标的平均值

E-step不断更新每个点是属于哪一个簇的期望值

M-step计算关于簇中心点的拟合函数值最大化对应坐标(argmax)-一般应用是求每个簇中所有数据点坐标的平均值去得到簇中心点坐标.

简单来说：在典型环境下，每一次重复的E-step和M-step都将会得到更好的聚类效果.

可视化

所以，根据这个思想，我们可以简单地实现一下其源码

from sklearn.metrics import pairwise_distances_argmin
import numpy as np
def find_clusters(X,n_clusters,rseed=2):
#随机选择簇中心点
	rng=np.random.RandomState(rseed)
	i=rng.permutation(X.shape[0])[:n_clusters]
	centers=X[i]

	while True:
#基于最近的中心指定标签
		labels=pairwise_distances_argmin(X,centers)

	#根据点的平均值找到新的中心
		new_centers=np.array([X[labels==i].mean(0) for i in range(n_clusters)])
#确认收敛
		if np.all(centers==new_centers):
			break
		centers=new_centers
	return centers,labels
a=np.array([[1,3],[1.5,3.2],[1.3,2.8],[3,1]])
centers,labels=find_clusters(a,2)

print(labels+1,'\n·······\n',centers)

4.K-means最佳簇数K值的确定

对于K均值聚类来说，如何确定簇数K值是一个至关重要的问题,为了解决这个问题，通过会选用探索法，即给定不同的k值，对比某些评估指标的变动情况，进而选择一个比较合理的k值. 本文将重点介绍两种非常实用的评估方法，即簇内离差平方和拐点法与轮廓系数法.

1.簇内离差平方和拐点法

顾名思义，簇内离差平方和拐点法的思想很简单，就是在不同的k值下计算簇内离差平方和，然后通过可视化的方法找到“拐点”所对应的k值.重点关注的是斜率的变化，当斜率由大突然变小时，并且之后的斜率变化缓慢，则认为突然变换的点就是寻找的目标点，因为继续随着簇数K的增加，聚类效果不再有大的变化.所以，整个核心就是找到”k的剧变点”.有点类似于常见的梯度下降，同样需要计算cost画图看k，不过这个不需要梯度，因为不需要方向.

这里我们随机生成三组二维正态分布数据，首先基于该数据绘制散点图,模拟的数据呈现三个簇.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['axes.unicode_minus']=False
from sklearn.cluster import KMeans
mean=np.array([[-2,-2],[2,2],[6,6]])	#设定三个二维正态分布的均值
cov=np.array([[[0.3,0],[0,0.3]],[[0.4,0],[0,0.4]],[[0.5,0],[0,0.5]]]) #设定三个二维正态分布的协方差矩阵 
X0=[];Y0=[];
for i in range(3):
	x,y=np.random.multivariate_normal(mean[i],cov[i],1000).T
	X0=np.hstack([X0,x]); Y0=np.hstack([Y0,y]) #水平堆叠，方便画图
plt.rc('font',size=16);
plt.rc('font',family='SimHei')
plt.rc('axies');plt.subplot(121)
plt.scatter(X0,Y0,marker='.') #画模拟数据散点图
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['axes.unicode_minus']=False
from sklearn.cluster import KMeans
mean=np.array([[-2,-2],[2,2],[6,6]])
cov=np.array([[[0.3,0],[0,0.3]],[[0.4,0],[0,0.4]],[[0.5,0],[0,0.5]]])
X0=[];Y0=[];
for i in range(3):
	x,y=np.random.multivariate_normal(mean[i],cov[i],1000).T
	X0=np.hstack([X0,x]); Y0=np.hstack([Y0,y])
plt.rc('font',size=16);
plt.rc('font',family='SimHei')
X=np.vstack([X0,Y0]).T
TSSE=[];k=10
for k in range(1,k+1):
    SSE=[]
    md=KMeans(k)
    md.fit(X)
    labels=md.labels_
    centers=md.cluster_centers_
    for label in  set(labels):
        SSE.append(np.sum((X[labels==label,:]-centers[label,])**2)) #labels=label来筛选数据集X,true才会进行计算
    TSSE.append(np.sum(SSE))
plt.subplot(122);plt.style.use('ggplot')
plt.plot(range(1,k+1),TSSE,'b*-')
plt.xlabel('簇的个数');plt.ylabel('簇内离差平方和之和')
plt.show()

经过上面的可视化过程后，我们可以看到，在簇为3个的那里k发生了剧烈的变化，也就是拐点，而后变动都较小，所以合理的k值应该为3.(不浪费算力资源)

2.轮廓系数法

​ 该方法综合考虑了簇的密集性和分散性两个信息，如果数据集被切割为理想的k个簇，那么对应的簇内样本会很密集，而簇间样本会很分散.

​ 如下图所示，假设数据集被拆分为三个簇G1,G2,G3,样本点i对应的 ai 值为所有G1中其他样本点与样本点i的距离平均值;样本点i对应的 bi 值分两步计算，首先计算该点分别到G2和G3中样本点的平均距离，然后将两个平均值中的最小值作为bi的度量.

定义样本点i的轮廓系数

​ $$S_i=\frac{b_i-a_i}{max(a_i,b_i)}$$

k个簇的总轮廓系数定义为所有样本点轮廓系数的平均值.

​ 当总轮廓系数小于0时，说明聚类效果不佳；

​ 当总轮廓系数接近于1时，说明簇内样本的平均距离非常小，而簇间的最近距离非常大，进而表示聚类效果非常理想.

​ 思想虽然较为简单，但实际当样本数据量较大时，运行时间会比较长。

​ 对于轮廓系数的计算，可以直接调用sklearn.metrics中的函数silhouette_score

​ 但是其需要聚类簇数必须大于等于2.

​ 我们同样利用上面的模拟数据来看看其于簇内离差平方和拐点法的效果是否是一样的:

import numpy as np; import matplotlib.pyplot as plt
from sklearn.cluster import KMeans; from sklearn import metrics
X=np.load("Pzdata11_1.npy")
S=[]; K=10
for k in range(2,K+1):
    md = KMeans(k); md.fit(X)
    labels = md.labels_;
    S.append(metrics.silhouette_score(X,labels,metric='euclidean'))  #计算轮廓系数
plt.rc('font',size=16); plt.rc('font',family='SimHei')
plt.plot(range(2,K+1), S, 'b*-')
plt.xlabel('簇的个数'); plt.ylabel('轮廓系数'); plt.show()

可以看到在簇数=3的时候，轮廓系数最接近于1，这与用上一种方法得到的结果是一样的.

所以，在实际情况中，我们可以更为保险地同时使用两种方法来判定结果，也可以只使用一种方法.

5.k均值聚类的应用.

在做k均值聚类时需要注意两点，一个是聚类前必须指定具体的簇数k值，如果k值是已知的，可以直接调用cluster子模块中的KMeans函数，对数据集进行分割; 如果k值是未知的，可以根据行业经验或前面介绍的两种方法确定合理的k值.

另一个是对原始数据做必要的标准化处理.由于k均值的思想是基于点之间的距离实现“物以类聚”的，所以**如果原始数据集存在量纲上的差异，就必须对其进行标准化的预处理.**数据集的标准化预处理可以借助sklearn子模块preprocessing中的scale函数或minmax_scale函数.scale函数的标准化公式为:

​ $$x^{\ast }=\frac{x-u}{\sigma }$$

minmax_scale函数的标准化公式为

​ $$x^{\ast }=\frac{x-x_{min}}{x_{max}-x_{min}}$$

$$其中，u，\sigma ,x_{min},x_{max}分别为x取值的均值、标准值、最小值和最大值.$$