原本準備自己總結一篇關於KMeans評估的文章，發現網上這篇已經總結的很全面了，所以收藏下。如果有其他方法也歡迎留言交流。原文來自微信公眾號：非凡wang咖
在Kmenas算法中，如何確定簇數K值是一個至關重要的問題，為了解決這個問題，通常會選用探索法，即給定不同的k值下，對比某些評估指標的變動情況，進而選擇一個比較合理的k值。本篇文章總結三種常用的k值選擇方法：簇內離差平方和拐點法，輪廓係數法和間隔統計量法）
拐點法簇內離差平方和拐點法的思想是：在不同的k值下計算簇內離差平方和，然後通過可視化的方法找到「拐點」所對應的k值。隨著簇數量的增加，簇中的樣本量會越來越少，通過可視化方法，重點關注曲線斜率的變化，當斜率由大突然變小時，並且之後的斜率變化緩慢，則認為突然變化的點就是拐點，也就是最佳的k值，因為繼續隨著簇數K的增加，聚類效果不再有大的變化。
接下來我們驗證這個方法，隨機生成三組二元正態分布數據，首先基於該數據繪製散點圖，如下代碼：
np.random.seed(1234)mean1 = [0.5, 0.5]cov1 = [[0.3, 0], [0, 0.3]]x1, y1 = np.random.multivariate_normal(mean1, cov1, 1000).T
mean2 = [0, 8]cov2 = [[0.3, 0], [0, 0.3]]x2, y2 = np.random.multivariate_normal(mean2, cov2, 1000).T
mean3 = [8, 4]cov3 = [[1.5, 0], [0, 1]]x3, y3 = np.random.multivariate_normal(mean3, cov3, 1000).T
plt.scatter(x1, y1)plt.scatter(x2, y2)plt.scatter(x3, y3)plt.show()如上圖，虛擬出來的數據呈現出三個簇，接下來基於這個虛擬數據，使用拐點法繪製簇的個數與總的簇內離差平方和之間的折線圖，確定最終的k值，代碼如下：
def k_SSE(X, clusters): """ 拐點法 繪製不同的k值和對應總的簇內離差平方和的折線圖 """ K = range(1, clusters+1) TSSE = [] for k in K: SSE = [] kmeans = KMeans(n_clusters=k) kmeans.fit(X) labels = kmeans.labels_ centers = kmeans.cluster_centers_ for label in set(labels): SSE.append(np.sum((X.loc[labels == label, ]-centers[label, :])**2)) TSSE.append(np.sum(SSE))
plt.rcParams['font.sans-serif'] = 'SimHei' plt.rcParams['axes.unicode_minus'] =False plt.style.use('ggplot') plt.plot(K, TSSE, 'b*-') plt.xlabel('簇的個數') plt.ylabel('簇內離差平方和之和') plt.show()
X = pd.DataFrame(np.concatenate([np.array([x1, y1]), np.array([x2, y2]), np.array([x3, y3])], axis=1).T)k_SSE(X, 15)如上圖，當簇的個數為3時，形成了一個明顯的「拐點」，因為K值從1到3時，折線的斜率都比較大，但是k值為4時斜率突然就降低了很多，並且之後的簇對應的斜率都變動很小，所以，合理的k值應該為3，與虛擬數據集的三個簇相吻合。
輪廓係數法輪廓係數法綜合考慮了簇的密集性和分散性，如果數據集被分割為理想的K個簇，那麼對應的簇內樣本會很密集，而簇間樣本會很分散。其公式如下：
其中，a(i) 是樣本i與同簇內其他樣本點距離的平均值，體現了簇內的密集性；b(i) 是樣本i與其他非同簇樣本點距離的平均值，然後從平均值中挑選出最小值，反映了簇間的分散性。
通過公式可知當S(i)接近於-1時，說明樣本i分配的不合理，需要將其分配到其他簇中；當S(i)近似為0時，說明樣本i落在了模糊地帶，即簇的邊界處；當S(i)近似為1時，說明樣本i的分配是合理的。
接下來我們就看看如何用輪廓係數解決我們的k取值問題，由於輪廓係數計算較複雜，所以我們直接使用sklearn中的metrics中的silhouette_score方法，需要注意的是該方法需要接受的聚類簇數必須大於等於2。代碼如下：def k_silhouette(X, clusters): """ 輪廓係數法 """ K = range(2, clusters+1) S = [] for k in K: kmeans = KMeans(n_clusters=k) kmeans.fit(X) labels = kmeans.labels_ S.append(metrics.silhouette_score(X, labels, metric='euclidean')) plt.rcParams['font.sans-serif'] = 'SimHei' plt.rcParams['axes.unicode_minus'] =False plt.style.use('ggplot') plt.plot(K, S, 'b*-') plt.xlabel('簇的個數') plt.ylabel('輪廓係數') plt.show()
k_silhouette(X, 15)如上圖，利用之前構造的虛擬數據，繪製了不同K值下對應的輪廓係數圖，當k取值為3時輪廓係數最大，且比較接近於1，說明應該把虛擬數據聚為3類比較合理。
間隔統計量法2000年Hastie等人提出了間隔統計量法（Gap Statistic方法），該方法可以適用與任何聚類算法，公式如下：
詳情參考地址：
https://blog.csdn.net/baidu_17640849/article/details/70769555
接下來我們構造自定義函數，繪製不同K值對應的間隙統計量折線圖：
def short_pair_wise_D(each_cluster): """ 計算簇內任意倆樣本之間的歐式距離Dk """ mu = each_cluster.mean(axis=0) Dk = sum(sum((each_cluster - mu) ** 2 * each_cluster.shape[0])) return Dk
def compute_Wk(data, classfication_result): """ 計算簇內的Wk值 """ Wk = 0 label_set = set(classfication_result) for label in label_set: each_cluster = data[classfication_result == label, :] Wk = Wk + short_pair_wise_D(each_cluster) / (2.0 * each_cluster.shape[0]) return Wk
def gap_statistic(X, B=10, K=range(1, 11), N_init=10): """ 間隔統計法 計算GAP統計量 """ X = np.array(X) shape = X.shape tops = X.max(axis=0) bots = X.min(axis=0) dists = np.matrix(np.diag(tops - bots)) rands = np.random.random_sample(size=(B, shape[0], shape[1])) for i in range(B): rands[i, :, :] = rands[i, :, :] * dists + bots
gaps = np.zeros(len(K)) Wks = np.zeros(len(K)) Wkbs = np.zeros((len(K), B)) for idxk, k in enumerate(K): k_means = KMeans(n_clusters=k) k_means.fit(X) classfication_result = k_means.labels_ Wks[idxk] = compute_Wk(X, classfication_result)
for i in range(B): Xb = rands[i, :, :] k_means.fit(Xb) classfication_result_b = k_means.labels_ Wkbs[idxk, i] = compute_Wk(Xb, classfication_result_b)
gaps = (np.log(Wkbs)).mean(axis=1) - np.log(Wks) sd_ks = np.std(np.log(Wkbs), axis=1) sk = sd_ks * np.sqrt(1 + 1.0 / B) gapDiff = gaps[:-1] - gaps[1:] + sk[1:]
plt.rcParams['font.sans-serif'] = 'SimHei' plt.rcParams['axes.unicode_minus'] = False plt.style.use('ggplot') plt.bar(np.arange(len(gapDiff)) + 1, gapDiff, color='steelblue') plt.xlabel('簇的個數') plt.ylabel('k的選擇標準') plt.show()
gap_statistic(X)如上圖，x軸代表了不同的簇數k，y軸代表k值選擇的判斷指標gapDiff，gapDiff首次出現正值時對應的k為3，所以對於虛擬的數據集來說，將其劃分為三個簇是比較合理的。
完整代碼連結：
https://github.com/jpegbert/MachineLearning/tree/master/kmeans/kmeans_evaluation
歡迎關注，一起學習
參考https://blog.csdn.net/baidu_17640849/article/details/70769555