數學建模及機器學習算法（一）：聚類-kmeans（Python及MATLAB實現，包括k值選取與聚類效果評估）

本文轉載自查看原文 2019-02-18 01:11 6310

一、聚類的概念

聚類分析是在數據中發現數據對象之間的關系，將數據進行分組，組內的相似性越大，組間的差別越大，則聚類效果越好。我們事先並不知道數據的正確結果(類標)，通過聚類算法來發現和挖掘數據本身的結構信息，對數據進行分簇(分類)。聚類算法的目標是，簇內相似度高，簇間相似度低

二、基本的聚類分析算法

　　1. K均值(K-Means)：
　　　　基於原型的、划分的距離技術，它試圖發現用戶指定個數(K)的簇。

　　2. 凝聚的層次距離：
　　　　思想是開始時，每個點都作為一個單點簇，然后，重復的合並兩個最靠近的簇，直到嘗試單個、包含所有點的簇。

　　3. DBSCAN:
　　　　一種基於密度的划分距離的算法，簇的個數有算法自動的確定，低密度中的點被視為噪聲而忽略，因此其不產生完全聚類。

三、距離量度

不同的距離量度會對距離的結果產生影響，常見的距離量度如下所示：

四、K-Means

在聚類算法中K-Means算法是一種最流行的、使用最廣泛的一種聚類算法，因為它的易於實現且計算效率也高。聚類算法的應用領域也是非常廣泛的，包括不同類型的文檔分類、音樂、電影、基於用戶購買行為的分類、基於用戶興趣愛好來構建推薦系統等。

優點：易於實現
缺點：可能收斂於局部最小值，在大規模數據收斂慢

算法思想：

1 選擇K個點作為初始質心  
2 repeat  
3     將每個點指派到最近的質心，形成K個簇  
4     重新計算每個簇的質心  
5 until 簇不發生變化或達到最大迭代次數

這里的重新計算每個簇的質心，更新過程是：首先找到與每個點距離最近的中心點，構成每個中心點划分的k個點集，然后對於每個點集，計算點的均值代替中心點。

如何計算是根據目標函數得來的，因此在開始時我們要考慮距離度量和目標函數。

考慮歐幾里得距離的數據，使用誤差平方和（Sum of the Squared Error,SSE）作為聚類的目標函數，兩次運行K均值產生的兩個不同的簇集，我們更喜歡SSE最小的那個：

k表示k個聚類中心，ci表示第幾個中心，dist表示的是歐幾里得距離。

前面說的我們更新質心是讓所有的點的平均值，這里就是SSE所決定的：

因此K-Means算法的實現步驟，主要分為四個步驟：

　　1、從樣本集合中隨機抽取k個樣本點作為初始簇的中心。

　　2、將每個樣本點划分到距離它最近的中心點所代表的簇中。

　　3、用各個簇中所有樣本點的中心點代表簇的中心點。

　　4、重復2和3，直到簇的中心點不變或達到設定的迭代次數或達到設定的容錯范圍。

五、k-means代碼實現

本文采用sklearn來實現一個k-means算法的應用，細節的底層實現可見文末第一個鏈接。

　　1.首先使用sklearn的數據集，數據集中包含150個隨機生成的點，樣本點分為三個不同的簇：

 1 from sklearn.datasets import make_blobs
 2 import matplotlib.pyplot as plt
 3  
 4  
 5 if __name__ == "__main__":
 6     '''
 7     n_samples：代表樣本點的個數
 8     n_features：表示每個樣本由兩個特征組成
 9     center：表示樣本點中心的個數(簇)
10     cluster_std：表示每個樣本簇方差的大小
11     '''
12     x,y = make_blobs(n_samples=150,n_features=2,centers=3,
13                      cluster_std=0.5,shuffle=True,random_state=0)
14     #繪點
15     plt.scatter(x[:,0],x[:,1],marker="o",color="blue")
16     #以表格的形式顯示
17     plt.grid()
18     plt.show()

　　效果如下：

　　2.下面使用sklearn內置的KMeans算法來實現對上面樣本點的聚類分析：

 1     from sklearn.cluster import KMeans
 2     '''
 3     n_clusters:設置簇的個數
 4     init:random表示使用Kmeans算法，默認是k-means++
 5     n_init:設置初始樣本中心的個數
 6     max_iter:設置最大迭代次數
 7     tol:設置算法的容錯范圍SSE(簇內誤平方差)
 8     '''
 9     kmeans = KMeans(n_clusters=3,init="random",n_init=10,max_iter=300,
10                     tol=1e-04,random_state=0)
11     y_km = kmeans.fit_predict(x)

六、k-means缺陷

k均值算法非常簡單且使用廣泛，但有一些缺點：

1. K值需要預先給定，屬於預先知識，很多情況下K值的估計是非常困難的，因此會有后面的k值確定。
2. K-Means算法對初始選取的聚類中心點是敏感的，不同的隨機種子點得到的聚類結果完全不同。可能只能得到局部的最優解，而無法得到全局的最優解。
3. K-Means算法並不是很所有的數據類型。它不能處理非球形簇、不同尺寸和不同密度的簇。
4. 對離群點的數據進行聚類時，K-Means也有問題。

七、K-Means++

K-Means算法需要隨機選擇初始化的中心點，如果中心點選擇不合適，可能會導致簇的效果不好或產生收斂速度慢等問題。解決這個問題一個比較合適的方法就是，在數據集上多次運行K-Means算法，根據簇內誤差平方和(SSE)來選擇性能最好的模型。除此之外，還可以通過K-Means++算法，讓初始的中心點彼此的距離盡可能的遠，相比K-Means算法，它能夠產生更好的模型。

K-Means++有下面幾個步驟組成：

　　1、初始化一個空的集合M，用於存儲選定的k個中心點

　　2、從輸入的樣本中隨機選擇第一個中心點μ，並將其加入到集合M中

　　3、對於集合M之外的任意樣本點x，通過計算找到與其距離最小的樣本d(x,M)

　　4、使用加權概率分布來隨機來隨機選擇下一個中心點μ

　　5、重復步驟2和3，直到選定k個中心點

　　6、基於選定的中心點執行k-means算法
使用sklearn來實現K-Means++，只需要將init參數設置為"k-means++"，默認設置是"k-means++"。

 1     km = KMeans(n_clusters=3,init="k-means++",n_init=10,max_iter=300,
 2                 tol=1e-04,random_state=0)
 3     #y_km中保存了聚類的結果
 4     y_km = km.fit_predict(x)
 5     #繪制不同簇的點
 6     plt.scatter(x[y_km==0,0],x[y_km==0,1],s=50,c="orange",marker="o",label="cluster 1")
 7     plt.scatter(x[y_km==1,0],x[y_km==1,1],s=50,c="green",marker="s",label="cluster 2")
 8     plt.scatter(x[y_km==2,0],x[y_km==2,1],s=50,c="blue",marker="^",label="cluster 3")
 9     #繪制簇的中心點
10     plt.scatter(km.cluster_centers_[:,0],km.cluster_centers_[:,1],s=250,marker="*",c="red"
11                 ,label="cluster center")
12     plt.legend()
13     plt.grid()
14     plt.show()

通過上面圖可以發現k-means++的聚類效果還不錯，簇的中心點，基本位於球心。

使用技巧：

1.在實際情況中使用k-means++算法可能會遇到，由於樣本的維度太高無法可視化，從而無法設定樣本的簇數。可視化問題可在聚類后顯示前用pca對數據降維顯示。

2.由於k-means算法是基於歐式距離來計算的，所以k-means算法對於數據的范圍比較敏感，所以在使用k-means算法之前，需要先對數據進行標准化，保證k-means算法不受特征量綱的影響。

八、Mini Batch k-Means

在原始的K-means算法中，每一次的划分所有的樣本都要參與運算，如果數據量非常大的話，這個時間是非常高的，因此有了一種分批處理的改進算法。

使用Mini Batch（分批處理）的方法對數據點之間的距離進行計算。
Mini Batch的好處：不必使用所有的數據樣本，而是從不同類別的樣本中抽取一部分樣本來代表各自類型進行計算。n 由於計算樣本量少，所以會相應的減少運行時間n 但另一方面抽樣也必然會帶來准確度的下降。

(其他關於k-means的優化還有很多，可參考鏈接或自行總結)

九、K值的確定

1.根據實際需要

2.肘部法則（Elbow Method）-見后文

3.輪廓系數（Silhouette Coefficient）-見后文

4.層次聚類

層次聚類是通過可視化然后人為去判斷大致聚為幾類，很明顯在共同父節點的一顆子樹可以被聚類為一個類

5.Canopy算法

肘部法則（Elbow Method）和輪廓系數（Silhouette Coefficient）來對k值進行最終的確定，但是這些方法都是屬於“事后”判斷的，而Canopy算法的作用就在於它是通過事先粗聚類的方式，為k-means算法確定初始聚類中心個數和聚類中心點。

與傳統的聚類算法(比如K-Means)不同，Canopy聚類最大的特點是不需要事先指定k值(即clustering的個數)，因此具有很大的實際應用價值。與其他聚類算法相比，Canopy聚類雖然精度較低，但其在速度上有很大優勢，因此可以使用Canopy聚類先對數據進行“粗”聚類，得到k值，以及大致的k個中心點，再使用K-Means進行進一步“細”聚類。所以Canopy+K-Means這種形式聚類算法聚類效果良好。

Canopy算法解析：

原始數據集合List按照一定的規則進行排序（這個規則是任意的，但是一旦確定就不再更改），初始距離閾值為T1、T2，且T1>T2（T1、T2的設定可以根據用戶的需要，或者使用交叉驗證獲得）。
在List中隨機挑選一個數據向量A，使用一個粗糙距離計算方式計算A與List中其他樣本數據向量之間的距離d。
根據第2步中的距離d，把d小於T1的樣本數據向量划到一個canopy中，同時把d小於T2的樣本數據向量從候選中心向量名單（這里可以理解為就是List）中移除。
重復第2、3步，直到候選中心向量名單為空，即List為空，算法結束。

算法原理比較簡單，就是對數據進行不斷遍歷，T2<dis<T1的可以作為中心名單，dis<T2的認為與canopy太近了，以后不會作為中心點，從list中刪除，這樣的話一個點可能屬於多個canopy。

canopy效果圖如下：

Canopy算法優勢：

Kmeans對噪聲抗干擾較弱，通過Canopy對比較小的NumPoint的Cluster直接去掉有利於抗干擾。
Canopy選擇出來的每個Canopy的centerPoint作為Kmeans比較科學。
只是針對每個Canopy的內容做Kmeans聚類，減少相似計算的數量。

Canopy算法缺點：

算法中 T1、T2（T2 < T1）的確定問題（也有專門的算法去描述，但可以自己多次試錯

Python實現：

 1 # -*- coding: utf-8 -*-
 2 import math
 3 import random
 4 import numpy as np
 5 import matplotlib.pyplot as plt
 6  
 7 class Canopy:
 8     def __init__(self, dataset):
 9         self.dataset = dataset
10         self.t1 = 0
11         self.t2 = 0
12  
13     # 設置初始閾值
14     def setThreshold(self, t1, t2):
15         if t1 > t2:
16             self.t1 = t1
17             self.t2 = t2
18         else:
19             print('t1 needs to be larger than t2!')
20  
21     # 使用歐式距離進行距離的計算
22     def euclideanDistance(self, vec1, vec2):
23         return math.sqrt(((vec1 - vec2)**2).sum())
24  
25     # 根據當前dataset的長度隨機選擇一個下標
26     def getRandIndex(self):
27         return random.randint(0, len(self.dataset) - 1)
28  
29     def clustering(self):
30         if self.t1 == 0:
31             print('Please set the threshold.')
32         else:
33             canopies = []  # 用於存放最終歸類結果
34             while len(self.dataset) != 0:
35                 rand_index = self.getRandIndex()
36                 current_center = self.dataset[rand_index]  # 隨機獲取一個中心點，定為P點
37                 current_center_list = []  # 初始化P點的canopy類容器
38                 delete_list = []  # 初始化P點的刪除容器
39                 self.dataset = np.delete(
40                     self.dataset, rand_index, 0)  # 刪除隨機選擇的中心點P
41                 for datum_j in range(len(self.dataset)):
42                     datum = self.dataset[datum_j]
43                     distance = self.euclideanDistance(
44                         current_center, datum)  # 計算選取的中心點P到每個點之間的距離
45                     if distance < self.t1:
46                         # 若距離小於t1，則將點歸入P點的canopy類
47                         current_center_list.append(datum)
48                     if distance < self.t2:
49                         delete_list.append(datum_j)  # 若小於t2則歸入刪除容器
50                 # 根據刪除容器的下標，將元素從數據集中刪除
51                 self.dataset = np.delete(self.dataset, delete_list, 0)
52                 canopies.append((current_center, current_center_list))
53         return canopies
54  
55  
56 def showCanopy(canopies, dataset, t1, t2):
57     fig = plt.figure()
58     sc = fig.add_subplot(111)
59     colors = ['brown', 'green', 'blue', 'y', 'r', 'tan', 'dodgerblue', 'deeppink', 'orangered', 'peru', 'blue', 'y', 'r',
60               'gold', 'dimgray', 'darkorange', 'peru', 'blue', 'y', 'r', 'cyan', 'tan', 'orchid', 'peru', 'blue', 'y', 'r', 'sienna']
61     markers = ['*', 'h', 'H', '+', 'o', '1', '2', '3', ',', 'v', 'H', '+', '1', '2', '^',
62                '<', '>', '.', '4', 'H', '+', '1', '2', 's', 'p', 'x', 'D', 'd', '|', '_']
63     for i in range(len(canopies)):
64         canopy = canopies[i]
65         center = canopy[0]
66         components = canopy[1]
67         sc.plot(center[0], center[1], marker=markers[i],
68                 color=colors[i], markersize=10)
69         t1_circle = plt.Circle(
70             xy=(center[0], center[1]), radius=t1, color='dodgerblue', fill=False)
71         t2_circle = plt.Circle(
72             xy=(center[0], center[1]), radius=t2, color='skyblue', alpha=0.2)
73         sc.add_artist(t1_circle)
74         sc.add_artist(t2_circle)
75         for component in components:
76             sc.plot(component[0], component[1],
77                     marker=markers[i], color=colors[i], markersize=1.5)
78     maxvalue = np.amax(dataset)
79     minvalue = np.amin(dataset)
80     plt.xlim(minvalue - t1, maxvalue + t1)
81     plt.ylim(minvalue - t1, maxvalue + t1)
82     plt.show()
83  
84  
85 if __name__ == "__main__":
86     dataset = np.random.rand(500, 2)  # 隨機生成500個二維[0,1)平面點
87     t1 = 0.6
88     t2 = 0.4
89     gc = Canopy(dataset)
90     gc.setThreshold(t1, t2)
91     canopies = gc.clustering()
92     print('Get %s initial centers.' % len(canopies))
93 showCanopy(canopies, dataset, t1, t2)

View Code

6.間隔統計量 Gap Statistic

根據肘部法則選擇最合適的K值有事並不是那么清晰，因此斯坦福大學的Robert等教授提出了Gap Statistic方法。

這里我們要繼續使用上面的 $D_k$ 。Gap Statistic的定義為:

$Gap_n(k)=E^*_n(log(D_k))-logD_k$

這里 ${E}(\log D_k)$ 指的是 $log D_k$ 的期望。這個數值通常通過蒙特卡洛模擬產生，我們在樣本里所在的矩形區域中（高維的話就是立方體區域）按照均勻分布隨機地產生和原始樣本數一樣多的隨機樣本，並對這個隨機樣本做K-Means，從而得到一個 $D_k$ 。如此往復多次，通常20次，我們可以得到20個 $log D_k$ 。對這20個數值求平均值，就得到了 ${E}(\log D_k)$ 的近似值。最終可以計算Gap Statisitc。而Gap statistic取得最大值所對應的K就是最佳的K。

Gap Statistic的基本思路是：引入參考的測值，這個參考值可以有Monte Carlo采樣的方法獲得。

$E^*_n(log(D_k)) =(1/B)\sum_{b=1}^{B}log(D_{kb}^*)$

B是sampling的次數。為了修正MC帶來的誤差，我們計算sk也即標准差來矯正Gap Statistic。

$w' = (1/B)\sum_{b=1}^{B}log(D^*_{kb})$

$sd(k) = \sqrt{(1/B)\sum_b(logD^*_{kb}-w')^2}$

$s_k = \sqrt{\frac{1+B}{B}}sd(k)$

選擇滿足 $Gap_k >= Gap_{k+1}-s_{k+1}$ 的最小的k作為最優的聚類個數。

Python實現：

 1 import scipy
 2 from  scipy.spatial.distance import euclidean
 3 from sklearn.cluster import KMeans as k_means
 4  
 5 dst = euclidean
 6  
 7 k_means_args_dict = {
 8     'n_clusters': 0,
 9     # drastically saves convergence time
10     'init': 'k-means++',
11     'max_iter': 100,
12     'n_init': 1,
13     'verbose': False,
14     # 'n_jobs':8
15 }
16  
17  
18 def gap(data, refs=None, nrefs=20, ks=range(1, 11)):
19     """
20     I: NumPy array, reference matrix, number of reference boxes, number of clusters to test
21     O: Gaps NumPy array, Ks input list
22  
23     Give the list of k-values for which you want to compute the statistic in ks. By Gap Statistic
24     from Tibshirani, Walther.
25     """
26     shape = data.shape
27  
28     if not refs:
29         tops = data.max(axis=0)
30         bottoms = data.min(axis=0)
31         dists = scipy.matrix(scipy.diag(tops - bottoms))
32         rands = scipy.random.random_sample(size=(shape[0], shape[1], nrefs))
33         for i in range(nrefs):
34             rands[:, :, i] = rands[:, :, i] * dists + bottoms
35     else:
36         rands = refs
37  
38     gaps = scipy.zeros((len(ks),))
39  
40     for (i, k) in enumerate(ks):
41         k_means_args_dict['n_clusters'] = k
42         kmeans = k_means(**k_means_args_dict)
43         kmeans.fit(data)
44         (cluster_centers, point_labels) = kmeans.cluster_centers_, kmeans.labels_
45  
46         disp = sum(
47             [dst(data[current_row_index, :], cluster_centers[point_labels[current_row_index], :]) for current_row_index
48              in range(shape[0])])
49  
50         refdisps = scipy.zeros((rands.shape[2],))
51  
52         for j in range(rands.shape[2]):
53             kmeans = k_means(**k_means_args_dict)
54             kmeans.fit(rands[:, :, j])
55             (cluster_centers, point_labels) = kmeans.cluster_centers_, kmeans.labels_
56             refdisps[j] = sum(
57                 [dst(rands[current_row_index, :, j], cluster_centers[point_labels[current_row_index], :]) for
58                  current_row_index in range(shape[0])])
59  
60         # let k be the index of the array 'gaps'
61         gaps[i] = scipy.mean(scipy.log(refdisps)) - scipy.log(disp)
62  
63     return ks, gaps

View Code

十、聚類效果評估

1.簇內誤方差(SSE)

在對簇的划分中，我們就使用了SSE作為目標函數來划分簇。當KMeans算法訓練完成后，我們可以通過使用inertia屬性來獲取簇內的誤方差，不需要再次進行計算。

 1     #用來存放設置不同簇數時的SSE值
 2     distortions = []
 3     for i in range(1,11):
 4         km = KMeans(n_clusters=i,init="k-means++",n_init=10,max_iter=300,tol=1e-4,random_state=0)
 5         km.fit(x)
 6         #獲取K-means算法的SSE
 7         distortions.append(km.inertia_)
 8     #繪制曲線
 9     plt.plot(range(1,11),distortions,marker="o")
10     plt.xlabel("簇數量")
11     plt.ylabel("簇內誤方差(SSE)")
12     plt.show()

可以使用圖形工具肘方法，根據簇的數量來可視化簇內誤方差。通過圖形可以直觀的觀察到k對於簇內誤方差的影響。也可以用來確定K值。

通過上圖可以發現，當簇數量為3的時候出現了肘型，這說明k取3是一個不錯的選擇。但不一定所有的問題都能用肘部法則來解決，如下圖右圖中，肘部不明顯。因此肘部法則只是一種可嘗試的方法。

2、輪廓圖定量分析聚類質量

輪廓分析(silhouette analysis)，使用圖形工具來度量簇中樣本的聚集程度，除k-means之外也適用於其他的聚類算法。通過三個步驟可以計算出當個樣本的輪廓系數(silhouette coefficient)：

　　1、將樣本x與簇內的其他點之間的平均距離作為簇內的內聚度a
　　2、將樣本x與最近簇中所有點之間的平均距離看作是與最近簇的分離度b
　　3、將簇的分離度與簇內聚度之差除以二者中比較大的數得到輪廓系數，計算公式如下

輪廓系數的取值在-1到1之間。當簇內聚度與分度離相等時，輪廓系數為0。當b>>a時，輪廓系數近似取到1，此時模型的性能最佳。

 1     km = KMeans(n_clusters=3,init="k-means++",n_init=10,max_iter=300,tol=1e-4,random_state=0)
 2     y_km = km.fit_predict(x)
 3     import numpy as np
 4     from matplotlib import cm
 5     from sklearn.metrics import silhouette_samples
 6     #獲取簇的標號
 7     cluster_labels = np.unique(y_km)
 8     #獲取簇的個數
 9     n_clusters = cluster_labels.shape[0]
10     #基於歐式距離計算輪廓系數
11     silhoutte_vals = silhouette_samples(x,y_km,metric="euclidean")
12     #設置y坐標的起始位置
13     y_ax_lower,y_ax_upper=0,0
14     yticks=[]
15     for i,c in enumerate(cluster_labels):
16         #獲取不同簇的輪廓系數
17         c_silhouette_vals = silhoutte_vals[y_km == c]
18         #對簇中樣本的輪廓系數由小到大進行排序
19         c_silhouette_vals.sort()
20         #獲取到簇中輪廓系數的個數
21         y_ax_upper += len(c_silhouette_vals)
22         #獲取不同顏色
23         color = cm.jet(i / n_clusters)
24         #繪制水平直方圖
25         plt.barh(range(y_ax_lower,y_ax_upper),c_silhouette_vals,
26                  height=1.0,edgecolor="none",color=color)
27         #獲取顯示y軸刻度的位置
28         yticks.append((y_ax_lower+y_ax_upper) / 2)
29         #下一個y軸的起點位置
30         y_ax_lower += len(c_silhouette_vals)
31     #獲取輪廓系數的平均值
32     silhouette_avg = np.mean(silhoutte_vals)
33     #繪制一條平行y軸的輪廓系數平均值的虛線
34     plt.axvline(silhouette_avg,color="red",linestyle="--")
35     #設置y軸顯示的刻度
36     plt.yticks(yticks,cluster_labels+1)
37     plt.ylabel("簇")
38     plt.xlabel("輪廓系數")
39     plt.show()

通過輪廓圖，我們能夠看出樣本的簇數以及判斷樣本中是否包含異常值。為了評價聚類模型的性能，可以通過評價輪廓系數，也就是圖中的紅色虛線進行評價。

類似SSE，也可以做出不同k值下的效果圖：

可以看到也是在聚類數為3時輪廓系數達到了峰值，所以最佳聚類數為3

十一、MATLAB實現

最后針對使用MATLAB的給出代碼，細節與上文類似：

代碼一：

生成隨機二維分布圖形，三個中心

 1 % 使用高斯分布（正態分布）
 2 % 隨機生成3個中心以及標准差
 3 s = rng(5,'v5normal');
 4 mu = round((rand(3,2)-0.5)*19)+1;
 5 sigma = round(rand(3,2)*40)/10+1;
 6 X = [mvnrnd(mu(1,:),sigma(1,:),200); ...
 7      mvnrnd(mu(2,:),sigma(2,:),300); ...
 8      mvnrnd(mu(3,:),sigma(3,:),400)];
 9 % 作圖
10 P1 = figure;clf;
11 scatter(X(:,1),X(:,2),10,'ro');
12 title('研究樣本散點分布圖')

分層聚類：

1 eucD = pdist(X,'euclidean');
2 clustTreeEuc = linkage(eucD,'average');
3 cophenet(clustTreeEuc,eucD);
4 
5 P3 = figure;clf;
6 [h,nodes] =  dendrogram(clustTreeEuc,20);
7 set(gca,'TickDir','out','TickLength',[.002 0],'XTickLabel',[]);

調用k-means,分成三類

1 [cidx3,cmeans3,sumd3,D3] = kmeans(X,3,'dist','sqEuclidean');
2 P4 = figure;clf;
3 [silh3,h3] = silhouette(X,cidx3,'sqeuclidean');

結果顯示：

 1 P5 = figure;clf
 2 ptsymb = {'bo','ro','go',',mo','c+'};
 3 MarkFace = {[0 0 1],[.8 0 0],[0 .5 0]};
 4 hold on
 5 for i =1:3
 6     clust = find(cidx3 == i);
 7     plot(X(clust,1),X(clust,2),ptsymb{i},'MarkerSize',3,'MarkerFace',MarkFace{i},'MarkerEdgeColor','black');
 8     plot(cmeans3(i,1),cmeans3(i,2),ptsymb{i},'MarkerSize',10,'MarkerFace',MarkFace{i});
 9 end
10 hold off

分別用等高線、分布圖、熱能圖和概率圖展示結果

 1 % 等高線
 2 options = statset('Display','off');
 3 gm = gmdistribution.fit(X,3,'Options',options);
 4 
 5 P6 = figure;clf
 6 scatter(X(:,1),X(:,2),10,'ro');
 7 hold on
 8 ezcontour(@(x,y) pdf(gm,[x,y]),[-15 15],[-15 10]);
 9 hold off
10 
11 P7 = figure;clf
12 scatter(X(:,1),X(:,2),10,'ro');
13 hold on
14 ezsurf(@(x,y) pdf(gm,[x,y]),[-15 15],[-15 10]);
15 hold off
16 view(33,24)
17 
18  
19 
20 熱能圖
21 cluster1 = (cidx3 == 1);
22 cluster3 = (cidx3 == 2);
23 % 通過觀察，K均值方法的第二類是gm的第三類
24 cluster2 = (cidx3 == 3);
25 % 計算分類概率
26 P = posterior(gm,X);
27 P8 = figure;clf
28 plot3(X(cluster1,1),X(cluster1,2),P(cluster1,1),'r.')
29 grid on;hold on
30 plot3(X(cluster2,1),X(cluster2,2),P(cluster2,2),'bo')
31 plot3(X(cluster3,1),X(cluster3,2),P(cluster3,3),'g*')
32 legend('第 1 類','第 2 類','第 3 類','Location','NW')
33 clrmap = jet(80); colormap(clrmap(9:72,:))
34 ylabel(colorbar,'Component 1 Posterior Probability')
35 view(-45,20);
36 % 第三類點部分概率值較低，可能需要其他數據來進行分析。
37 
38 % 概率圖
39 P9 = figure;clf
40 [~,order] = sort(P(:,1));
41 plot(1:size(X,1),P(order,1),'r-',1:size(X,1),P(order,2),'b-',1:size(X,1),P(order,3),'y-');
42 legend({'Cluster 1 Score' 'Cluster 2 Score' 'Cluster 3 Score'},'location','NW');
43 ylabel('Cluster Membership Score');
44 xlabel('Point Ranking');

View Code

AIC准則尋找最優分類

1 AIC = zeros(1,4);
2 NlogL = AIC;
3 GM = cell(1,4);
4 for k = 1:4
5     GM{k} = gmdistribution.fit(X,k);
6     AIC(k)= GM{k}.AIC;
7     NlogL(k) = GM{k}.NlogL;
8 end
9 [minAIC,numComponents] = min(AIC);

代碼二（簡單版）：

%隨機初始化中心點，可以隨機取數據集中的任意k個點
function centroids = kMeansInitCentroids(X, K)
    centroids = zeros(K, size(X, 2));
    randidx = randperm(size(X, 1));
    centroids = X(randidx(1 : K), :);
end

%對每個點，找到其所屬的中心點，即距離其最近的中心點。 返回一個向量，為每個點對應的中心點id。
function idx = findClosestCentroids(X, centroids)
    K = size(centroids, 1);
    for i = 1 : size(X, 1)
        min_dis = +inf;
        for j = 1 : K
            now_dis = sum((X(i, :) - centroids(j, :)).^2);
            if now_dis < min_dis
                min_dis = now_dis;
                idx(i) = j;
            end
        end
    end
end

%用每個中心點統領的那類點的均值，更新中心點。
function centroids = computeCentroids(X, idx, K)
    [m n] = size(X);
    cnt = zeros(K, n);
    for i = 1 : m
        cnt(idx(i), :) = cnt(idx(i), :) + 1;
        centroids(idx(i), :) = centroids(idx(i), :)  + X(i, :);
    end
    cnt = 1 ./ cnt;
    centroids = cnt .* centroids;
end

centroids = kMeansInitCentroids(X, k);
for iter = 1 : iterations 
    idx = findClosestCentroids(X, centroids);
    centroids = computeMeans(X, idx, K);
end

十二、參考鏈接：

深入理解K-Means聚類算法

K-Means算法詳細介紹(SSE、輪廓分析)

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 python機器學習——kmeans聚類算法機器學習-KMeans聚類 K值以及初始類簇中心點的選取機器學習-KMeans聚類 K值以及初始類簇中心點的選取【機器學習】KMeans 聚類算法原理與實現機器學習：Python實現聚類算法(一)之K-Means 機器學習——聚類算法的評估指標【機器學習】：Kmeans均值聚類算法原理(附帶Python代碼實現) 機器學習（十）—聚類算法（KNN、Kmeans、密度聚類、層次聚類）機器學習之--kmeans聚類簡單算法實例機器學習——詳解經典聚類算法Kmeans