前言
最近研究中要用到聚類。K-means是第一個(gè)去研究的算法池充，覺得這個(gè)簡單的算法也很有意思眠冈，這期推送把它的來龍去脈記錄下略号，順便分享！
原理與數(shù)據(jù)
k-means的原理洋闽，相對簡單玄柠，總結(jié)起來就一句話：依據(jù)“組內(nèi)距離最大，組間距離最小”原則把點(diǎn)劃分為k類诫舅。不贅述羽利，原理可參考：
《機(jī)器學(xué)習(xí)》周志華
《機(jī)器學(xué)習(xí)實(shí)戰(zhàn)》
網(wǎng)絡(luò)博客、百度
因?yàn)檠芯繑?shù)據(jù)不公開刊懈，本期數(shù)據(jù)采用隨機(jī)構(gòu)造數(shù)據(jù)这弧，如下代碼：
##【代碼1】隨機(jī)產(chǎn)生點(diǎn)##
import numpy as np
import matplotlib.pyplot as plt
def get_data(real_center= [(1, 1), (1, 2), (2, 2), (2, 1)]):
# 先在四個(gè)中心點(diǎn)附近產(chǎn)生一堆數(shù)據(jù)
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) 四個(gè)點(diǎn)為中心產(chǎn)生了隨機(jī)分布的點(diǎn)，可視化一下虚汛，會(huì)生成四個(gè)可見的簇匾浪。
調(diào)包
<pre style="background:white">最偷懶的方法就是spss操作、是Python調(diào)包卷哩、R語言調(diào)包了蛋辈，比如說Python中內(nèi)置的sklearn.cluster 就提供了k-means算法。程序調(diào)用如下：</pre>
##【代碼2】sklearn調(diào)用##
from sklearn.cluster import KMeans
datas=get_data()
kmeans_model = KMeans(n_clusters=4, n_init=15).fit(datas) # 隨機(jī)選擇初始中心
labels = kmeans_model.labels_
plt.scatter(datas[:, 0], datas[:, 1], c=labels, cmap='rainbow')
plt.show()
這是內(nèi)置k-means出圖的結(jié)果将谊，大致上是已經(jīng)滿足聚類出結(jié)果的需求了冷溶，符合構(gòu)造數(shù)據(jù)集的原始特點(diǎn)。但是我想觀察它是怎樣演化的尊浓，并且輸出一些自定義的參數(shù)逞频，如鄧恩指數(shù)。這種方法限制了我們對結(jié)果的分析栋齿，還是自己編寫吧Ｃ缯汀！瓦堵！
自編寫
編寫的代碼如下（代碼還沒來得及“美化”基协，但算法可行的）
##【代碼3】k-means自編寫##
import numpy as np
import matplotlib.pyplot as plt
def distance(vector1,vector2):
#把距離函數(shù)單獨(dú)拿出來的原因是距離的定義可能會(huì)變化的，比如經(jīng)緯度和普通數(shù)字的不同
op1 = np.sqrt(np.sum(np.square(vector1 - vector2)))
op2 = np.linalg.norm(vector1 - vector2)
return op2
def my_Kmeans(points,class_num):
...
:輸入待分類的點(diǎn)集谷丸、分類數(shù)目堡掏；輸出分類的結(jié)果（【1，0刨疼，2泉唁、、揩慕、】）亭畜，可視化表達(dá)
:points的類型為np數(shù)組
...
i=1
rand_index = np.random.choice(range(len(points)), size=class_num, replace=False)#replace代表不重復(fù)
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)]#用來存放點(diǎn)
class_results=[]#用來存放類別數(shù)字，如1迎卤，2拴鸵，
#將全部點(diǎn)歸屬到類
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)
#從新計(jì)算各類中心
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數(shù)組
:param class_reslts: 分類結(jié)果，是一個(gè)列表
...
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)#產(chǎn)生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)
自編寫的方便多了劲藐，來個(gè)演化過程的動(dòng)圖
模型改進(jìn)Kmedoids
（一） 理論
因?yàn)閗-means算法是一個(gè)最基礎(chǔ)的基于點(diǎn)劃分的模型八堡，今20年來很多研究者對原模型進(jìn)行的改進(jìn)，比如說：
對初始值敏感： K-means++聘芜，intelligent K-means兄渺， genetic K-means；
對噪聲敏感：k-medoids, k-medians,
僅僅適用于數(shù)值型：k-modes
不能解決非凸數(shù)據(jù)：kernel K-means
（二） 實(shí)戰(zhàn)
改進(jìn)算法很多汰现，根據(jù)我個(gè)人實(shí)際需求（數(shù)據(jù)集中有明顯離群點(diǎn)）挂谍，接著解了下k-medoids算法。當(dāng)存在噪音和孤立點(diǎn)時(shí), K-medoids 比 K-means 更健壯瞎饲。
k-means與k-medoids之間的差異就是可以理解為對于數(shù)據(jù)樣本的平均值和中位數(shù)之間的差異：前者的取值范圍可以是連續(xù)空間中的任意值口叙，而后者的取值卻只能是數(shù)據(jù)樣本范圍中的樣本。
舉個(gè)例子說明：在聚類中心更替的時(shí)候嗅战，左圖很明顯有兩類中心妄田。但是在右圖中出現(xiàn)了個(gè)異常點(diǎn)，導(dǎo)致中心變到了中心點(diǎn)1仗哨，影響后續(xù)的聚類形庭。但是k-medois中會(huì)選擇中心點(diǎn)2作為中心（median oids），顯然更加合理厌漂。改進(jìn)算法的代碼如下：
##【代碼4】k-medoids自編寫##
import numpy as np
import matplotlib.pyplot as plt
from sklearn import metrics
def initial(points,k):
#給定初始的點(diǎn)萨醒，返回隨機(jī)的k個(gè)初始點(diǎn)（點(diǎn)不重復(fù)）
i=0
return_centre=[]
while i<k:
rand_index = np.random.choice(range(len(points)), size=1, replace=False) # replace代表不重復(fù)
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):
#把距離函數(shù)單獨(dú)拿出來的原因是距離的定義可能會(huì)變化的，比如經(jīng)緯度和普通數(shù)字的不同
op1 = np.sqrt(np.sum(np.square(vector1 - vector2)))
op2 = np.linalg.norm(vector1 - vector2)
return op1
def get_better_centre(class_points):
# 輸入為二維列表[[],[],[]...],每個(gè)子類存放了一類分類點(diǎn)苇倡；輸出為每個(gè)子類的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):
'''
輸入待分類的點(diǎn)集、分類數(shù)目旨椒；輸出分類的結(jié)果（【1晓褪，0，2综慎、涣仿、、】）示惊，可視化表達(dá)
points的類型為np數(shù)組
'''
i=1
centre=initial(points,class_num)
# print(type(centre))
centre_diff=(np.sum(centre==0)-class_num)
data_dim = points.shape[1] # 每個(gè)點(diǎn)數(shù)據(jù)有多少個(gè)維度
# print(centre)
final_label=[[]]
while (centre_diff !=0 and i <=1000):#迭代完500次之后好港，即使出現(xiàn)不收斂的情況也停止迭代
# global class_
i+=1
class_points=[[]for i in range(class_num)]#用來存放點(diǎn)
class_=[]#用來存放類別數(shù)字，如1米罚，2钧汹，
#將全部點(diǎn)歸屬到類
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_
#從新計(jì)算各類中心
new_centre=get_better_centre(class_points)#調(diào)用函數(shù)，求新的中心點(diǎn)
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格式的數(shù)據(jù)录择，分類數(shù)拔莱，重復(fù)運(yùn)行次數(shù)
eval=[]
all_labs=[]
for i in range(repetition):
# print('重復(fù)了%次',i)
lables=my_kmedoids(np_points,class_num)
# print(np_points)
# print("數(shù)據(jù)",len(np_points))
# print("分類",len(lables))
silhouettescore = metrics.silhouette_score(np_points, lables, metric='euclidean') # 指標(biāo)值在[-1,1]之間碗降，越接近1越好
all_labs.append(lables)
eval.append(silhouettescore)
index=eval.index(max(eval))
return all_labs[index]
if __name__ == '__main__':
#構(gòu)造一個(gè)數(shù)據(jù)集
class_num=3
repetition=20#隨機(jī)初始運(yùn)行的次數(shù)，最后取最優(yōu)的結(jié)果
np.random.seed(0)
a=np.random.randint(low=0,high=20,size=20)#產(chǎn)生10對（x,y）
A=a.reshape(-1,2)
#調(diào)用函數(shù)
lables=k_medoids_(A,class_num,repetition)
print(type(lables))
##可視化分類結(jié)果
plt.scatter(A[:, 0], A[:, 1], c=lables, cmap='rainbow')
plt.show()
##繪制動(dòng)圖##
k=4
from clustering.k_medoids import k_medoids_
datas=get_data()
labels=k_medoids_(datas,k,repetition=1)
plt.show()
雙動(dòng)圖比較塘秦，迭代次數(shù)比較讼渊。
本期內(nèi)容結(jié)束，關(guān)注【交通科研Lab】微信公眾號(hào)嗤形，點(diǎn)擊“閱讀原文”精偿，獲取全部代碼。歡迎關(guān)注下一期《如何科學(xué)地選擇聚類數(shù)K赋兵？》
如果您覺得本篇文章對您有用，請不吝點(diǎn)贊搔预！您的鼓勵(lì)與支持是我創(chuàng)作的最大動(dòng)力Ｅ凇！
