现成调包>底层编写>算法改进——我对k-means由浅入深的理解
2019-12-19 本文已影响0人
交通科研Lab
前言
最近研究中要用到聚类。K-means是第一个去研究的算法，觉得这个简单的算法也很有意思，这期推送把它的来龙去脉记录下，顺便分享！
原理与数据
k-means的原理，相对简单，总结起来就一句话：依据“组内距离最大，组间距离最小”原则把点划分为k类。不赘述，原理可参考：
《机器学习》周志华
《机器学习实战》
网络博客、百度
因为研究数据不公开，本期数据采用随机构造数据，如下代码：
##【代码1】随机产生点##
import numpy as np
import matplotlib.pyplot as plt
def get_data(real_center= [(1, 1), (1, 2), (2, 2), (2, 1)]):
# 先在四个中心点附近产生一堆数据
point_number = 50
points_x = []
points_y = []
for center in real_center:
offset_x, offset_y = np.random.randn(point_number) * 0.3, np.random.randn(point_number) * 0.25
x_val, y_val = center[0] + offset_x, center[1] + offset_y
points_x.append(x_val)
points_y.append(y_val)
points_x = np.concatenate(points_x)
points_y = np.concatenate(points_y)
p_list = np.stack([points_x, points_y], axis=1)
return p_list
real_center = [(1, 1), (1, 2), (2, 2), (2, 1)]
points=get_data(real_center)
plt.scatter(np.array(real_center)[:,0],np.array(real_center)[:,1],c='r',marker='o',s=100)
plt.scatter(points[:,0],points[:,1])
plt.show()
我们以 (1, 1), (1, 2), (2, 2), (2, 1) 四个点为中心产生了随机分布的点，可视化一下，会生成四个可见的簇。
调包
<pre style="background:white">最偷懒的方法就是spss操作、是Python调包、R语言调包了，比如说Python中内置的sklearn.cluster 就提供了k-means算法。程序调用如下：</pre>
##【代码2】sklearn调用##
from sklearn.cluster import KMeans
datas=get_data()
kmeans_model = KMeans(n_clusters=4, n_init=15).fit(datas) # 随机选择初始中心
labels = kmeans_model.labels_
plt.scatter(datas[:, 0], datas[:, 1], c=labels, cmap='rainbow')
plt.show()
这是内置k-means出图的结果，大致上是已经满足聚类出结果的需求了，符合构造数据集的原始特点。但是我想观察它是怎样演化的，并且输出一些自定义的参数，如邓恩指数。这种方法限制了我们对结果的分析，还是自己编写吧！！！
自编写
编写的代码如下（代码还没来得及“美化”，但算法可行的）
##【代码3】k-means自编写##
import numpy as np
import matplotlib.pyplot as plt
def distance(vector1,vector2):
#把距离函数单独拿出来的原因是距离的定义可能会变化的，比如经纬度和普通数字的不同
op1 = np.sqrt(np.sum(np.square(vector1 - vector2)))
op2 = np.linalg.norm(vector1 - vector2)
return op2
def my_Kmeans(points,class_num):
...
:输入待分类的点集、分类数目；输出分类的结果（【1，0，2、、、】），可视化表达
:points的类型为np数组
...
i=1
rand_index = np.random.choice(range(len(points)), size=class_num, replace=False)#replace代表不重复
centre=points[rand_index, :]
centre_diff=(np.sum(centre==0)-class_num)
while centre_diff !=0 and i <=500:
i+=1
class_points=[[]for i in range(class_num)]#用来存放点
class_results=[]#用来存放类别数字，如1，2，
#将全部点归属到类
for point in points:
all_dis=[]
for sta in centre:
dis=distance(point,sta)
all_dis.append(dis)
index=all_dis.index(min(all_dis))
class_results.append(index)
class_points[index].append(point)
#从新计算各类中心
new_centre=[]
for each in class_points:
each_centre=np.array(each).mean(axis=0)
new_centre.append(list(each_centre))
minus =centre-np.array(new_centre)
centre=np.array(new_centre)
centre_diff=(np.sum(minus==0)-class_num)
return class_results
#plt.scatter(X_scatter[:, 0], X_scatter[:, 1], c=lables, cmap='rainbow')
#这种方法参考下，没必要这么麻烦的
def my_plt(points,class_reslts,class_num):
...
:param points: points的类型为np数组
:param class_reslts: 分类结果，是一个列表
...
for clss in range(class_num):
filter=points[np.array(class_reslts)==clss]
print(filter)
plt.scatter(filter[:, 0],filter[:, 1])
plt.show()
if __name__=='__main__':
class_num=3
np.random.seed(0)
a=np.random.randint(low=0,high=20,size=20)#产生10对（x,y）
A=a.reshape(-1,2)
class_reslts=my_Kmeans(A,class_num)
print(class_reslts)
my_plt(A,class_reslts,class_num)
自编写的方便多了，来个演化过程的动图
模型改进Kmedoids
（一） 理论
因为k-means算法是一个最基础的基于点划分的模型，今20年来很多研究者对原模型进行的改进，比如说：
对初始值敏感： K-means++，intelligent K-means， genetic K-means；
对噪声敏感：k-medoids, k-medians,
仅仅适用于数值型：k-modes
不能解决非凸数据：kernel K-means
（二） 实战
改进算法很多，根据我个人实际需求（数据集中有明显离群点），接着解了下k-medoids算法。当存在噪音和孤立点时, K-medoids 比 K-means 更健壮。
k-means与k-medoids之间的差异就是可以理解为对于数据样本的平均值和中位数之间的差异：前者的取值范围可以是连续空间中的任意值，而后者的取值却只能是数据样本范围中的样本。
举个例子说明：在聚类中心更替的时候，左图很明显有两类中心。但是在右图中出现了个异常点，导致中心变到了中心点1，影响后续的聚类。但是k-medois中会选择中心点2作为中心（median oids），显然更加合理。改进算法的代码如下：
##【代码4】k-medoids自编写##
import numpy as np
import matplotlib.pyplot as plt
from sklearn import metrics
def initial(points,k):
#给定初始的点，返回随机的k个初始点（点不重复）
i=0
return_centre=[]
while i<k:
rand_index = np.random.choice(range(len(points)), size=1, replace=False) # replace代表不重复
temp =list(points[rand_index, :][0])
if temp not in return_centre:
return_centre.append(temp)
i+=1
return np.array(return_centre)
def distance(vector1,vector2):
#把距离函数单独拿出来的原因是距离的定义可能会变化的，比如经纬度和普通数字的不同
op1 = np.sqrt(np.sum(np.square(vector1 - vector2)))
op2 = np.linalg.norm(vector1 - vector2)
return op1
def get_better_centre(class_points):
# 输入为二维列表[[],[],[]...],每个子类存放了一类分类点；输出为每个子类的medoids，列表形式[]
new_centre=[]
for each_class in class_points:
evaluation=[]
for each_points in each_class:
diff=np.array(each_class)-np.array(each_points)
diff_square = [[diff[i][j] ** 2 for j in range(len(diff[i]))] for i in range(len(diff))]
diff_sum=sum(sum(i) for i in diff_square)
evaluation.append(diff_sum)
min_index=evaluation.index(min(evaluation))
point=each_class[min_index]
new_centre.append(point)
return new_centre
def my_kmedoids(points,class_num):
'''
输入待分类的点集、分类数目；输出分类的结果（【1，0，2、、、】），可视化表达
points的类型为np数组
'''
i=1
centre=initial(points,class_num)
# print(type(centre))
centre_diff=(np.sum(centre==0)-class_num)
data_dim = points.shape[1] # 每个点数据有多少个维度
# print(centre)
final_label=[[]]
while (centre_diff !=0 and i <=1000):#迭代完500次之后，即使出现不收敛的情况也停止迭代
# global class_
i+=1
class_points=[[]for i in range(class_num)]#用来存放点
class_=[]#用来存放类别数字，如1，2，
#将全部点归属到类
for point in points:
all_dis=[]
for sta in centre:
dis=distance(point,sta)
all_dis.append(dis)
index=all_dis.index(min(all_dis))
class_.append(index)
class_points[index].append(point)
final_label[0]=class_
#从新计算各类中心
new_centre=get_better_centre(class_points)#调用函数，求新的中心点
minus =centre-np.array(new_centre)
centre=np.array(new_centre)
# print(np.sum(minus==0))
centre_diff=(np.sum(minus==0)/data_dim-class_num)
final_class=final_label[0]
# print(i)
return final_class
def k_medoids_(np_points,class_num,repetition=50):#np格式的数据，分类数，重复运行次数
eval=[]
all_labs=[]
for i in range(repetition):
# print('重复了%次',i)
lables=my_kmedoids(np_points,class_num)
# print(np_points)
# print("数据",len(np_points))
# print("分类",len(lables))
silhouettescore = metrics.silhouette_score(np_points, lables, metric='euclidean') # 指标值在[-1,1]之间，越接近1越好
all_labs.append(lables)
eval.append(silhouettescore)
index=eval.index(max(eval))
return all_labs[index]
if __name__ == '__main__':
#构造一个数据集
class_num=3
repetition=20#随机初始运行的次数，最后取最优的结果
np.random.seed(0)
a=np.random.randint(low=0,high=20,size=20)#产生10对（x,y）
A=a.reshape(-1,2)
#调用函数
lables=k_medoids_(A,class_num,repetition)
print(type(lables))
##可视化分类结果
plt.scatter(A[:, 0], A[:, 1], c=lables, cmap='rainbow')
plt.show()
##绘制动图##
k=4
from clustering.k_medoids import k_medoids_
datas=get_data()
labels=k_medoids_(datas,k,repetition=1)
plt.show()
双动图比较，迭代次数比较。
本期内容结束，关注【交通科研Lab】微信公众号，点击“阅读原文”，获取全部代码。欢迎关注下一期《如何科学地选择聚类数K？》
如果您觉得本篇文章对您有用，请不吝点赞！您的鼓励与支持是我创作的最大动力！！
