前言

最近研究中要用到聚類。K-means是第一個去研究的算法，覺得這個簡單的算法也很有意思，這期推送把它的來龍去脈記錄下，順便分享！

原理與數(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)]):
    # 先在四個中心點(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) 四個點(diǎn)為中心產(chǎn)生了隨機(jī)分布的點(diǎn)，可視化一下，會生成四個可見的簇。

調(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é)果的分析，還是自己編寫吧?。。?/p>

自編寫

編寫的代碼如下（代碼還沒來得及“美化”，但算法可行的）

##【代碼3】k-means自編寫##
import numpy as np
import matplotlib.pyplot as plt
def distance(vector1,vector2):
    #把距離函數(shù)單獨(dú)拿出來的原因是距離的定義可能會變化的，比如經(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é)果，是一個列表
    ...
    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)

自編寫的方便多了，來個演化過程的動圖

模型改進(jìn)Kmedoids

（一） 理論

因?yàn)閗-means算法是一個最基礎(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ù)我個人實(shí)際需求（數(shù)據(jù)集中有明顯離群點(diǎn)），接著解了下k-medoids算法。當(dāng)存在噪音和孤立點(diǎn)時, K-medoids 比 K-means 更健壯。

k-means與k-medoids之間的差異就是可以理解為對于數(shù)據(jù)樣本的平均值和中位數(shù)之間的差異：前者的取值范圍可以是連續(xù)空間中的任意值，而后者的取值卻只能是數(shù)據(jù)樣本范圍中的樣本。

舉個例子說明：在聚類中心更替的時候，左圖很明顯有兩類中心。但是在右圖中出現(xiàn)了個異常點(diǎn)，導(dǎo)致中心變到了中心點(diǎn)1，影響后續(xù)的聚類。但是k-medois中會選擇中心點(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個初始點(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ú)拿出來的原因是距離的定義可能會變化的，比如經(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):
    # 輸入為二維列表[[],[],[]...],每個子類存放了一類分類點(diǎn)；輸出為每個子類的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]  # 每個點(diǎn)數(shù)據(jù)有多少個維度
    # 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)造一個數(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()

##繪制動圖##
k=4
from clustering.k_medoids import k_medoids_
datas=get_data()
labels=k_medoids_(datas,k,repetition=1)
plt.show()

雙動圖比較，迭代次數(shù)比較。

本期內(nèi)容結(jié)束，關(guān)注【交通科研Lab】微信公眾號，點(diǎn)擊“閱讀原文”，獲取全部代碼。歡迎關(guān)注下一期《如何科學(xué)地選擇聚類數(shù)K？》

如果您覺得本篇文章對您有用，請不吝點(diǎn)贊！您的鼓勵與支持是我創(chuàng)作的最大動力！！