大家在学习算法的时候会学习到关于Kmeans的算法，但是网络和很多机器学习算法书中关于Kmeans的算法理论核心一样，但是代码实现过于复杂，效率不高，不方便阅读。这篇文章首先列举出Kmeans核心的算法过程，并且会给出如何最大限度的在不用for循环的前提下，利用numpy, pandas的高效的功能来完成Kmeans算法。这里会用到列表解析，它是相当于速度更快的for循环，标题里指出的无for loop指的是除了列表解析解析以外不用for循环，来完成Kmeans算法。
一般在pythonshujuqingxi/' style='color:#000;font-size:inherit;'>python数据清洗中，数据量大的情况下，for循环的方法会使的数据处理的过程特别慢，效率特别低。一个很好的解决方法就是使用numpy，pandas自带的高级功能，不仅可以使得代码效率大大提高，还可以使得代码方便理解阅读。这里在介绍用numpy，pandas来进行Kmeans算法的同时，也是带大家复习一遍numpy，pandas用法。
创建k个点作为初始质⼼心（通常是随机选择）
当任意一个点的簇分配结果发生改变时:
对数据集中的每个点:
对每个质⼼:
计算质⼼与数据点之间的距离
将数据点分配到据其最近的簇
对每个簇，计算簇中所有点的均值并将均值作为新的质⼼点
直到簇不再发⽣变化或者达到最大迭代次数
SSE = \sum_{i=1}^k\sum_{x\in C_{i}}(c_{i} - x)^2SSE=i=1∑kx∈Ci∑(ci−x)2
C_{i}指的是第i个簇, x是i个簇中的点，c_{i}是第i个簇的质心Ci指的是第i个簇,x是i个簇中的点，ci是第i个簇的质心
import numpy as np import pandas as pd import matplotlib as mpl import matplotlib.pyplot as plt from sklearn.datasets import make_blobs import seaborn as sns
#r = np.random.randint(1,100)
r = 4
#print(r)
k = 3
x , y = make_blobs(n_samples = 51,
cluster_std = [0.3, 0.3, 0.3],
centers = [[0,0],[1,1],[-1,1]]
,random_state = r
)
sim_data = pd.DataFrame(x, columns = ['x', 'y'])
sim_data['label'] = y
sim_data.head(5)
data = sim_data.copy()
plt.scatter(sim_data['x'], sim_data['y'], c = y)
上图是一个随机生成的2维的数据，可以用来尝试完成Kmeans的代码。
实际过程中，Kmeans需要能运行在多维的数据上，所以下面的代码部分，会考虑多维的数据集，而不是仅仅2维的数据。
这里的严格意义上不是随机的生成k个质心点，而是取出每个特征的最大值最小值，在最大值和最小值中取出一个随机数作为质心点的一个维度
def initial_centers(datasets, k = 3):
#首先将datasets的特征名取出来，这里需要除去label那一列
cols = datasets.columns
data_content = datasets.loc[:, cols != 'label']
#直接用describe的方法将每一列的最小值最大值取出来
range_info = data_content.describe().loc[['min','max']]
#用列表解析的方法和np.random.uniform的方法生成k个随机的质心点
#np.random.uniform(a, b, c) 随机生成在[a,b)区间里的3个数
#对每个特征都做此操作
k_randoms = [np.random.uniform(range_info[i]['min'],
range_info[i]['max'], k)
for i in range_info.columns]
centers = pd.DataFrame(k_randoms, index = range_info.columns)
return centers.T
centers = initial_centers(data, k = 3) centers
xy00.1225750.0217621-0.9225961.3675042-0.677202-0.411821
将每一个中心点取出来，然后使用pandas的广播的功能，可以直接将所有的实例和其中一个质心点相减。如下图，下图中是给出相加的例子，而我们的例子是减法。
所以对于一个DataFrame来说，比如说这里的只包含x和y的data，假设我们的质心是c = [1,1]，可以用以下的方式来给出所有的实例点的x和y和点(1,1)之间的差值。注意，这里的c可以是list，也可以是numpy array，甚至可以是元组。
$$
$$
算出每个实例的每个特征和质心点的差距之后，则需要将所有的数平方一下，然后按每一行加起来则给出了每一个实例点到质心的距离了
$$
$$
用的方法就是使用np.power(data - c, 2).sum(axis = 1)
def cal_distant(dataset, centers):
#选出不是label的那些特征列
data = dataset.loc[:, dataset.columns != 'label']
#使用列表解析式的格式，对centers表里的每一行也就是每一个随机的质心点，都算一遍所有的点到该质心点的距离，并且存入一个list中
d_to_centers = [np.power(data - centers.loc[i], 2).sum(axis = 1)
for i in centers.index]
#所有的实例点到质心点的距离都已经存在了list中，则可以直接带入pd.concat里面将数据拼起来
return pd.concat(d_to_centers, axis = 1)
d_to_centers = cal_distant(data, centers) d_to_centers.head(5)
01200.1533653.9355460.52828611.9878790.0880062.46244420.0279772.3617530.79500430.5434105.1832830.56569641.5055142.2482644.031165
当每个实例点都和中心点计算好距离后，对于每个实例点找出最近的那个中心点,可以用np.where的方法，但是pandas已经提供更加方便的方法，用idxmin和idxmax,这2个函数可以直接给出DataFrame每行或者每列的最小值和最大值的索引，设置axis = 1则是想找出对每个实例点来说，哪个质心点离得最近。
curr_group = d_to_centers.idxmin(axis=1)
这个时候，每个点都有了新的group，这里我们则需要开始更新我们的3个中心点了。对每一个临时的簇来说，算出X的平均， 和Y的平均，就是这个临时的簇的中心点。
centers = data.loc[:, data.columns != 'label'].groupby(curr_group).mean() centers
xy00.5484680.5234741-1.0036801.0449552-0.125490-0.475373
这样我们新的质心点就得到了，只是这个时候的算法还是没有收敛的，需要将上面的步骤重复多次。
Kmeans代码迭代部分就完成了，将上面的步骤做成一个函数，做成函数后，方便展示Kmeans的中间过程。
def iterate(dataset, centers):
#计算所有的实例点到所有的质心点之间的距离
d_to_centers = cal_distant(dataset, centers)
#得出每个实例点新的类别
curr_group = d_to_centers.idxmin(axis=1)
#算出当前新的类别下每个簇的组内误差
SSE = d_to_centers.min(axis = 1).sum()
#给出在新的实例点类别下，新的质心点的位置
centers = dataset.loc[:, dataset.columns != 'label'].groupby(curr_group).mean()
return curr_group, SSE, centers
curr_group, SSE, centers = iterate(data,centers)
centers, SSE
( x y 0 0.892579 0.931085 1 -1.003680 1.044955 2 0.008740 -0.130172, 19.041432436034352)
最后需要判断什么时候迭代停止，可以判断SSE差值不变的时候，算法停止
#创建一个空的SSE_list,用来存SSE的，第一个位置的数为0，无意义，只是方便收敛时最后一个SSE和上一个SSE的对比
SSE_list = [0]
#初始化质心点
centers = initial_centers(data, k = 3)
#开始迭代
while True:
#每次迭代中得出新的组，组内误差，和新的质心点，当前的新的质心点会被用于下一次迭代
curr_group, SSE, centers = iterate(data,centers)
#检查这一次算出的SSE和上一次迭代的SSE是否相同，如果相同，则收敛结束
if SSE_list[-1] == SSE:
break
#如果不相同，则记录SSE，进入下一次迭代
SSE_list.append(SSE)
SSE_list
[0, 37.86874675507244, 11.231524142566894, 8.419267088238051]
算法完成了，将所有的代码整合在一起
def initial_centers(datasets, k = 3):
cols = datasets.columns
data_content = datasets.loc[:, cols != 'label']
range_info = data_content.describe().loc[['min','max']]
k_randoms = [np.random.uniform(range_info[i]['min'],
range_info[i]['max'], k)
for i in range_info.columns]
centers = pd.DataFrame(k_randoms, index = range_info.columns)
return centers.T
def cal_distant(dataset, centers):
data = dataset.loc[:, dataset.columns != 'label']
d_to_centers = [np.power(data - centers.loc[i], 2).sum(axis = 1)
for i in centers.index]
return pd.concat(d_to_centers, axis = 1)
def iterate(dataset, centers):
d_to_centers = cal_distant(dataset, centers)
curr_group = d_to_centers.idxmin(axis=1)
SSE = d_to_centers.min(axis = 1).sum()
centers = dataset.loc[:, dataset.columns != 'label'].groupby(curr_group).mean()
return curr_group, SSE, centers
def Kmeans_regular(data, k = 3):
SSE_list = [0]
centers = initial_centers(data, k = k)
while True:
curr_group, SSE, centers = iterate(data,centers)
if SSE_list[-1] == SSE:
break
SSE_list.append(SSE)
return curr_group, SSE_list, centers
上面的函数已经完成，当然这里推荐大家尽量写成class的形式更好，这里为了方便观看，则用简单的函数完成。
最后的函数是Kmeans_regular函数，这个函数里面包含了上面所有的函数。现在需要测试Kmeans_regular代码对于多特征的数据集鸢尾花数据集，是否也能进行Kmeans聚类算法
from sklearn.datasets import load_iris data_dict = load_iris() iris = pd.DataFrame(data_dict.data, columns = data_dict.feature_names) iris['label'] = data_dict.target
curr_group, SSE_list, centers = Kmeans_regular(iris.copy(), k = 3)
np.array(SSE_list)
array([ 0. , 589.73485975, 115.8301874 , 83.29216169,
79.45325846, 78.91005674, 78.85144143])
pd.crosstab(iris['label'], curr_group)
col_0012label0500010482201436
np.diag(pd.crosstab(iris['label'], curr_group)).sum() / iris.shape[0]
0.8933333333333333
最后可以看出我们的代码是可以适用于多特征变量的数据集，并且对于鸢尾花数据集来说，对角线上的数是预测正确的个数，准确率大约为90%。
在完成代码后，还是需要讨论一下，为什么我们的代码的算法是那样的，这个算法虽然看起来很有逻辑，但是它到底是从哪里来的。
这个时候，我们就需要从Kmeans的损失函数出发来解释刚才提出的问题。对于无监督学习算法来说，也是有一个损失函数。而我们的Kmeans的中间过程的逻辑，就是从最小化Kmeans的损失函数的过程。
假设我们有一个数据集{x_1, x_2, ..., x_N}x1,x2,...,xN， 每个样本实例点x有多个特征。我们的目标是将这个数据集通过某种方式切分成K份，或者说我们最后想将每个样本点标上一个类别(簇)，且总共有K个类别，使得每个样本点到各自的簇中心点的距离最小，并且u_kuk来表示各个簇的中心点。
我们还需要一些其他的符号，比如说r_{nk}rnk, 它的值是0或者1。下标k代表的是第k个簇，下标n表示的是第n个样本点。
举例说明，加入当前K=3，k的可取1，2，3。对于第一个实例点n = 1来说它属于第3个簇，所以
r_{n=1, k = 1} = 0rn=1,k=1=0
r_{n=1, k = 2} = 0rn=1,k=2=0
r_{n=1, k = 3} = 1rn=1,k=3=1
这个也可以把想象成独热编码。
将上面的符号解释完了后，以下就是损失函数。这里是使用了求和嵌套了求和的公式，并且也引入了刚才所提到了r_{nk}rnk。这个损失函数其实很好理解，在给定的k个中心点u_kuk以及分配好了各个实例点属于哪一个簇之后，计算各个实例点到各自的簇中心点的距离，距离平方以下并且相加起来，就是损失函数。这个公式其实也就是在算簇内误差和。
C = \sum_{n=1}^N\sum_{k=1}^K r_{nk} (x_n - u_k)^2C=n=1∑Nk=1∑Krnk(xn−uk)2
那怎么来最小化这个损失函数呢，用的就是EM算法，这个算法总体来说分2个步骤，Expectation和Maximization，对Kmeans来说M应该说是Minimization
Expection:
保持u_{k}uk不变，也就是各个簇的中心点的位置不变，计算各个实例点到哪个u_{k}uk最近，将各个实例点划分到离各自最近的那个簇里面去，从而使得整体SSE降低。
Minimization:
保持当前实例点的簇的类别不变，为了整体降低损失函数，可以对每个簇内的损失函数公式做微分。由于当前我们的各个点的类别是不变的，变的是u_{k}uk，所以做的微分是基于u_{k}uk的
\frac{d}{du_{k}}\sum_{k=1}^K r_{nk} (x_n - u_k)^2 = 0dukdk=1∑Krnk(xn−uk)2=0
-2\sum_{k=1}^K r_{nk} (x_n - u_k) = 0−2k=1∑Krnk(xn−uk)=0
u_{k} = \frac{\sum_{n} r_{nk} x_{n}}{\sum_{n} r_{nk}}uk=∑nrnk∑nrnkxn
得出来的u_{k}uk其实就是在算各个簇内的新的中心点，也就是对各个簇内所有的实例点的各个特征做平均数。
这时候得到新的中心点u_{k}uk, 紧接着再到E阶段，保持u_{k}uk更新簇类别，再到M阶段，保持簇类别不变更新u_{k}uk，不断的迭代知道SSE不变位置。这个就是Kmeans的算法过程。下面将用plotly可视化，动态展示Kmeans的过程。
使用之前写好的函数，然后将数据的中间过程通过plotly展示出来。因为代码比较长，所以这里就不展示代码了。由于当前是一个markdown，这里放入一个gif图片用来展示最后的Kmeans中间过程。
对于这个数据集来看的话，我们的Kmeans算法可以使得每一个点最终可以找到各自的簇，但是这个算法也是有缺陷的，比如以下例子。
假如说现在有4个簇的话，Kmeans算法不一定能使最后的SSE最小。对于2列的数据集来说，我们取2组随机的质心点来做对比。
第一组为设置seed为5的时候，以下为演示的结果。
从上面的动图可以看出一共用了8次迭代，才收敛。那加入我们的seed为1的话，随机的质心点的分布会变的很离谱，会导致下面的结果。这里我们加快动画的速度。
这里用34次，数据才迭代收敛，并且可以看出，在迭代的过程中，差点陷入了一个局部最小的一个情况。所以对于复杂的数据来说的话，我们最后看到迭代的次数会明显的增加。
假如说我们的数据集再变的集中一点，其中的2个簇，稍微近一点，我们会看到以下的结果。
所以在这次迭代的过程中，我们明显看到其中有个质心点消失了，原因就是因为由于点的分布的原因和初始质心点的原因，最开始随机生成的一个离所有的点都最远的质心点，由于它离所有的点都最远，所以导致了在迭代的过程中，没有任何一个点属于这个质心点，最后导致这个点消失了。所以这个就是Kmeans算法的缺陷，那怎么来优化这个算法了，我们可以利用BiKmeans算法。
数据分析咨询请扫描二维码
