K均值聚類

哈爾濱工程大學-537

演算法原理：

K均值是發現給定資料集的k$k$個簇的演算法。簇個數k$k$是使用者給定的，每一個簇通過其質心(centroid)，即簇中所有點的中心來描述。

K均值演算法的工作流程是：首先隨機確定k$k$個初始點作為質心。然後將資料集中的每個點分配到一個簇，具體來講，為每個點找距其最近的質心，並將其分配給該之心所對應的簇。這一步完成後，每個簇的質心更新為該簇所有點的平均值。

程式碼實現 ：

一、基礎實現

1）首先匯入所需要的各個庫
%matplotlib inline是為了使接下來的繪圖能夠顯示在Jupyter notebook的瀏覽器上。

import
 
 numpy as np
import math
import matplotlib.pyplot as plt
%matplotlib inline

2）定義載入樣本點資料的函式load_dataset，輸入引數為檔名；
將檔案按行讀取，每一行為一個列表，依次追加到空列表data_mat中，最後以矩陣的形式返回所得。

def load_dataset(file_name):
    data_mat = []
    with open(file_name) as fr:
        lines = fr.readlines()
    for line in lines:
        cur_line = line.strip().split("\t"
 
)
        flt_line = list(map(lambda x:float(x), cur_line))
        data_mat.append(flt_line)
    return np.array(data_mat)

3）呼叫load_dataset函式，引數為將要開啟的檔案的路徑和檔名，得到樣本點的特徵矩陣矩陣data_set。

data_set = load_dataset(r"D:\python_data\MachineLearningInaction\machinelearninginaction\Ch10\testSet.txt")
print(data_set)
 

將其打印出來，得到一個80行，2列的矩陣，其中前9行如下圖：

4）利用python的繪相簿繪製樣本點的散點圖。其中橫座標為”factor1”，縱座標為”factor2”。

point_x = data_set[:,0]
point_y = data_set[:,1]        
fig, ax = plt.subplots(figsize=(10,5))
ax.scatter(point_x, point_y, s=30, c="r", marker="o", label="sample point")
ax.legend()
ax.set_xlabel("factor1")
ax.set_ylabel("factor2")

繪製出的散點圖如下圖所示：

5）接下來定義一個計算兩個向量之間歐氏距離的函式dist_eclud，已知向量X=(x1,x2...xn)$X=(x_1,x_2...x_n)$和Y=(y1,y2...yn)$Y=(y_1,y_2...y_n)$之間的歐氏距離為
(x1−y1)2+(x2−y2)2+...+(xn−yn)2−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−√$\sqrt{(x_1-y_1{)}^2+(x_2-y_2{)}^2+...+(x_n-y_n{)}^2}$

def dist_eclud(vecA, vecB):
    vec_square = []
    for element in vecA - vecB:
        element = element ** 2
        vec_square.append(element)
    return sum(vec_square) ** 0.5

6）構建k$k$個隨機質心。輸入引數為資料集data_set和想要構建的質心個數k。隨機質心必須要在整個資料集的邊界之內。首先可以找到資料集每一維的最小和最大值。然後得到每一維的取值範圍。用0到1之間的隨機數和取值範圍相乘，再用最小值加上該乘積，就可以得到在每一維取值範圍內的隨機數。

# 構建k個隨機質心
def rand_cent(data_set, k):
    n = data_set.shape[1]    
    centroids = np.zeros((k, n))
    for j in range(n):
        min_j = float(min(data_set[:,j]))
        range_j = float(max(data_set[:,j])) - min_j
        centroids[:,j] = (min_j + range_j * np.random.rand(k, 1))[:,0]
    return centroids

7）演算法的主函式，輸入引數為資料集和簇的個數。
首先初始化樣本點的簇分配矩陣，有80行2列，第一列為該樣本點的簇分配索引，第二列為該樣本點到該簇質心的歐氏距離。
當任意一個點的簇分配發生變化時，迭代執行以下操作：遍歷每個樣本點，計算樣本點i到各個質心的距離，找到最小距離，將該質心所在簇編號分配給該樣本點。遍歷完所有樣本點後，重新計算每個簇的質心。直到所有樣本點的簇分配都不再發生變化時迭代停止。
最後返回質心和樣本點的簇分配矩陣

def Kmeans(data_set, k):    
    m = data_set.shape[0]   
    cluster_assment = np.zeros((m, 2))  
    centroids = rand_cent(data_set, k)  
    cluster_changed = True        
    while cluster_changed:       
        cluster_changed = False   
        for i in range(m):        
            min_dist = np.inf; min_index = -1   
            for j in range(k):     
                dist_ji = dist_eclud(centroids[j,:], data_set[i,:])
                if dist_ji < min_dist:              
                    min_dist = dist_ji; min_index = j   
            if cluster_assment[i,0] != min_index:    
                cluster_changed = True      
            cluster_assment[i,:] = min_index, min_dist**2   
        for cent in range(k):   
            pts_inclust = data_set[np.nonzero(list(map(lambda x:x==cent, cluster_assment[:,0])))]
            centroids[cent,:] = np.mean(pts_inclust, axis=0)
    return centroids, cluster_assment            

8）呼叫函式Kmean，輸入資料集為data_set，簇數量為4個，得到收斂後的質心和簇分配矩陣(簇分配矩陣略)。

my_centroids, my_cluster_assment = Kmeans(data_set, 4)
print(my_centroids)
# print(my_cluster_assment)

得到的4個質心為：

9）為了更加清晰直觀，用繪相簿繪製散點圖，其中，四個質心用黑色▲標出。

point_x = data_set[:,0]
point_y = data_set[:,1]  
cent_x = my_centroids[:,0]
cent_y = my_centroids[:,1]
fig, ax = plt.subplots(figsize=(10,5))
ax.scatter(point_x, point_y, s=30, c="r", marker="o", label="sample point")
ax.scatter(cent_x, cent_y, s=100, c="black", marker="v", label="centroids")
ax.legend()
ax.set_xlabel("factor1")
ax.set_ylabel("factor2")

得到的散點圖如下圖所示：

二、二分K-均值聚類

上述演算法的聚類簇的數目k$k$是使用者自定義的，那麼使用者如何才能知道k$k$的選擇是否正確？如何才能知道生成的簇比較好呢？在包含簇分配結果的矩陣中儲存著每個點的誤差，即該點到簇質心的距離平方值。

一種用於度量聚類效果的指標使SSE(Sum of Squared Error，誤差平方和)，SSE值越小表示資料點越接近它們的質心，聚類效果也就越好。

於是就可以採用二分K-聚類的方法，具體程式碼實現如下：

1）定義一個bi_Kmeans函式，兩個輸入引數分別為資料集和聚類簇數量。

該函式首先建立一個矩陣cluster_assment來儲存資料集中每個點的簇分配結果及平方誤差，然後計算整個資料集的質心，並使用一個列表cent_list來保留所有的質心。得到上述質心之後，可以遍歷資料集中所有點來計算每個點到質心的誤差值。這些誤差值將會在後面用到。

接下來程式進入while迴圈，該迴圈會不停對簇進行劃分，直到得到想要的簇數目為止。可以通過考察簇列表中的值來獲得當前簇的數目。然後遍歷所有簇來決定最佳的簇進行劃分。

為此需要比較劃分前後的SSE。一開始將最小SSE設定為無窮大，然後遍歷簇列表cent_list中的每一個簇。對每個簇，將該簇中的所有點堪稱一個小的資料集pts_incurr_cluster。將pts_incurr_cluster輸入到函式Kmeans()中進行處理(k=2)。K-均值演算法會生成兩個質心(簇)，同時給出每個簇的誤差值。這些誤差與剩餘資料集的誤差之和作為本次劃分的誤差。如果該劃分的SSE值最小，則儲存本次劃分。

一旦決定了要劃分的簇，接下來就可以實際執行劃分操作。劃分操作很容易，只需要將劃分的簇中所有點的簇分配結果進行修改即可。當使用Kmeans()函式並且制定簇數目為2時，會得到兩個編號分別為0和1的結果簇。需要將這些簇編號修改為被劃分簇以及新加簇的編號，該過程可以通過兩個陣列過濾器來完成。

最後，新的簇分配結果被更新，新的質心會被新增到cent_list中。

當while迴圈結束時，同Kmeans()函式一樣，函式返回之心列表與簇分配結果。

# 二分K均值聚類的主函式，輸入引數為資料集和要分的簇的個數
def bi_Kmeans(data_set, k):   
    m = data_set.shape[0]   # 得到data_set的行數，即資料集的個數
    cluster_assment = np.zeros((m,2)) # 初始化樣本點的簇分配矩陣，第一列為簇分配索引，第二列為歐氏距離平方
    centroid0 = np.mean(data_set, axis=0)  # 按列計算均值，即找到初始質心
    cent_list = [centroid0]
    for j in range(m):   # 對於每個樣本點
        cluster_assment[j,1] = dist_eclud(centroid0, data_set[j,:])**2 # 計算該樣本點的誤差平方
    while (len(cent_list) < k):   # 當已有的簇個數小於k時，迭代執行以下程式碼
        lowestSSE = np.inf     # 初始化誤差平方和SSE的最小值

        # 找到對哪個簇進行劃分可以最大程度降低SSE值
        for i in range(len(cent_list)):   # 遍歷每個已有的簇
            # 得到屬於該簇的所有樣本點
            pts_incurr_cluster = \
            data_set[np.nonzero(list(map(lambda x:x==i, cluster_assment[:,0])))]
            # 將該簇的所有樣本點通過函式Kmean進行劃分(k=2),得到劃分後的質心和簇分配矩陣
            centroid_mat, split_clust_ass = Kmeans(pts_incurr_cluster, 2)    
            sse_split = np.sum(split_clust_ass[:,1])  # 得到劃分後的誤差平方和
            # 得到其他樣本點的誤差平方和
            sse_not_split = \
            np.sum(cluster_assment[np.nonzero(list(map(lambda x:x!=i, cluster_assment[:,0]))),1])
            if (sse_split + sse_not_split) < lowestSSE:  # 如果總的誤差平方和小於lowestSSE
                best_cent_to_split = i                   # 則儲存本次劃分
                best_new_cents = centroid_mat
                best_clust_ass = split_clust_ass
                lowestSSE = sse_split + sse_not_split
        # 對最大程度降低SSE值的簇進行劃分

        # 將劃分後得到的編號為0的結果簇的編號修改為原最大簇編號+1，即len(cent_list)
        best_clust_ass[np.nonzero(list(map(lambda x:x==1, best_clust_ass[:,0]))), 0] = len(cent_list)  
         # 將劃分後得到的編號為1的結果簇的編號修改為被劃分的簇的編號
        best_clust_ass[np.nonzero(list(map(lambda x:x==0, best_clust_ass[:,0]))), 0] = best_cent_to_split 
        cent_list[best_cent_to_split] = best_new_cents[0,:] # 更新被劃分的簇的質心
        cent_list.append(best_new_cents[1,:])  # 新增一個新的簇的質心
        cluster_assment[np.nonzero(list(map(lambda x:x==best_cent_to_split, cluster_assment[:,0]))),:] = \
            best_clust_ass  # 將簇分配矩陣中屬於被劃分的簇的樣本點的簇分配進行更新
    return np.array(cent_list), cluster_assment

2）接下來呼叫bi_Kmeans()函式，得到三個質心和簇分配結果(簇分配結果略)，在呼叫bi_Kmeans()函式前，先執行load_dataset()函式載入資料。

data_set2 = load_dataset(r"D:\python_data\MachineLearningInaction\machinelearninginaction\Ch10\testSet2.txt")

cent_list, cluster_assment = bi_Kmeans(data_set2, 3)
print(cent_list)

3）最後，用繪相簿繪製散點圖，橫座標為”factor1”，縱座標為”factor2”，其中三個質心用黑色X標出。

point_x = data_set2[:,0]
point_y = data_set2[:,1] 
cent_x = cent_list[:,0]
cent_y = cent_list[:,1]
fig, ax = plt.subplots(figsize=(10,5))
ax.scatter(point_x, point_y, s=30, c="r", marker="o", label="sample point")
ax.scatter(cent_x, cent_y, s=100, c="black", marker="x", label="centroids")
ax.legend()
ax.set_xlabel("factor1")
ax.set_ylabel("factor2")

所繪製的散點圖如下：