機器學習之確定最佳聚類數目的10種方法

2021-01-09 雷鋒網

雷鋒網(公眾號：雷鋒網) AI科技評論按，本文作者貝爾塔，原文載於知乎專欄數據分析與可視化，雷鋒網 AI科技評論獲其授權發布。

在聚類分析的時候確定最佳聚類數目是一個很重要的問題，比如kmeans函數就要你提供聚類數目這個參數，總不能兩眼一抹黑亂填一個吧。之前也被這個問題困擾過，看了很多博客，大多泛泛帶過。今天把看到的這麼多方法進行匯總以及代碼實現並儘量弄清每個方法的原理。
數據集選用比較出名的wine數據集進行分析

library(gclus)
data(wine)
head(wine)

Loading required package: cluster

因為我們要找一個數據集進行聚類分析，所以不需要第一列的種類標籤信息，因此去掉第一列。
同時注意到每一列的值差別很大，從1到100多都有，這樣會造成誤差，所以需要歸一化，用scale函數

dataset <- wine[,-1] #去除分類標籤
dataset <- scale(dataset)

去掉標籤之後就可以開始對數據集進行聚類分析了，下面就一一介紹各種確定最佳聚類數目的方法

判定方法

1.mclust包

mclust包是聚類分析非常強大的一個包，也是上課時老師給我們介紹的一個包，每次導入時有一種科技感 :) 幫助文檔非常詳盡，可以進行聚類、分類、密度分析
Mclust包方法有點「暴力」，聚類數目自定義，比如我選取的從1到20，然後一共14種模型，每一種模型都計算聚類數目從1到20的BIC值，最終確定最佳聚類數目，這種方法的思想很直接了當，但是弊端也就顯然易見了——時間複雜度太高，效率低。

library(mclust)
m_clust <- Mclust(as.matrix(dataset), G=1:20) #聚類數目從1一直試到20
summary(m_clust)

Gaussian finite mixture model fitted by EM algorithm


Mclust EVE (ellipsoidal, equal volume and orientation) model with 3 components:


log.likelihood   n  df       BIC       ICL
      -3032.45 178 156 -6873.257 -6873.549


Clustering table:
1  2  3

63 51 64
可見該函數已經把數據集聚類為3種類型了。數目分別為63、51、64。再畫出14個指標隨著聚類數目變化的走勢圖

plot(m_clust, "BIC")

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/9c2e7ed5c2a689fc2caf452112362db6.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;531&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;527&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;531&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic4.zhimg.com/v2-5357b2dfa1d132078fb52f7ec8ea8faf_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/9c2e7ed5c2a689fc2caf452112362db6.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

下表是這些模型的意義

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/92cc4b3be70b912cc2f43a3e8d204f9d.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;715&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;532&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;715&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic2.zhimg.com/v2-5e9833600f6cff8e5233c12e61dc15b9_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/92cc4b3be70b912cc2f43a3e8d204f9d.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

它們應該分別代表著相關性(完全正負相關——對角線、稍強正負相關——橢圓、無關——圓)等參數的改變對應的模型，研究清楚這些又是非常複雜的問題了，先按下表，知道BIC值越大則說明所選取的變量集合擬合效果越好。上圖中除了兩個模型一直遞增，其他的12模型數基本上都是在聚類數目為3的時候達到峰值，所以該算法由此得出最佳聚類數目為3的結論。
mclust包還可以用於分類、密度估計等，這個包值得好好把玩。

注意：此BIC並不是貝葉斯信息準則！！！

最近上課老師講金融模型時提到了BIC值，說BIC值越小模型效果越好，頓時想起這裡是在圖中BIC極大值為最佳聚類數目，然後和老師探討了這個問題，之前這裡誤導大家了，Mclust包裡面的BIC並不是貝葉斯信息準則。

1.維基上的貝葉斯信息準則定義

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/3b66ee7c7a243bf9c7c3b22a23ad55af.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;292&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;53&amp;amp;amp;quot; class=&amp;amp;amp;quot;content_image&amp;amp;amp;quot; width=&amp;amp;amp;quot;292&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/3b66ee7c7a243bf9c7c3b22a23ad55af.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

與log(likelihood)成反比，極大似然估計是值越大越好，那麼BIC值確實是越小模型效果越好

2.Mclust包中的BIC定義[3]

這是Mclust包裡面作者定義的「BIC值」，此BIC非彼BIC，這裡是作者自己定義的BIC，可以看到，這裡的BIC與極大似然估計是成正比的，所以這裡是BIC值越大越好，與貝葉斯信息準則值越小模型越好的結論並不衝突

2.Nbclust包

Nbclust包是我在《R語言實戰》上看到的一個包，思想和mclust包比較相近，也是定義了幾十個評估指標，然後聚類數目從2遍歷到15（自己設定），然後通過這些指標看分別在聚類數為多少時達到最優，最後選擇指標支持數最多的聚類數目就是最佳聚類數目。

library(NbClust)
set.seed(1234) #因為method選擇的是kmeans，所以如果不設定種子，每次跑得結果可能不同
nb_clust <- NbClust(dataset,  distance = "euclidean",
       min.nc=2, max.nc=15, method = "kmeans",
       index = "alllong", alphaBeale = 0.1)

*** : The Hubert index is a graphical method of determining the number of clusters.
               In the plot of Hubert index, we seek a significant knee that corresponds to a
               significant increase of the value of the measure i.e the significant peak in Hubert
               index second differences plot.

*** : The D index is a graphical method of determining the number of clusters.
               In the plot of D index, we seek a significant knee (the significant peak in Dindex
               second differences plot) that corresponds to a significant increase of the value of
               the measure.

*******************************************************************
* Among all indices:                                                
* 5 proposed 2 as the best number of clusters
* 16 proposed 3 as the best number of clusters
* 1 proposed 10 as the best number of clusters
* 1 proposed 12 as the best number of clusters
* 1 proposed 14 as the best number of clusters
* 3 proposed 15 as the best number of clusters

                  ***** Conclusion *****                            

* According to the majority rule, the best number of clusters is  3

*******************************************************************

barplot(table(nb_clust$Best.nc[1,]),xlab = "聚類數",ylab = "支持指標數")

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/343439c0965ccbde22dbc1d2185b1a13.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;533&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;530&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;533&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic3.zhimg.com/v2-a08e3cf366ada95c9fa43d429a00f16a_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;&amp;lt;a href=&amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/343439c0965ccbde22dbc1d2185b1a13.png&amp;amp;quot;/&amp;gt;&amp;quot; _src=&amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/343439c0965ccbde22dbc1d2185b1a13.png&amp;amp;quot;/&amp;gt;&amp;quot; _href=&amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/343439c0965ccbde22dbc1d2185b1a13.png&amp;amp;quot;/&amp;gt;&amp;quot;&amp;gt;https://static.leiphone.com/uploads/new/article/pic/201710/343439c0965ccbde22dbc1d2185b1a13.png&amp;amp;amp;quot;/&amp;amp;amp;gt;&amp;lt;/a&amp;gt;

可以看到有16個指標支持最佳聚類數目為3，5個指標支持聚類數為2，所以該方法推薦的最佳聚類數目為3.

3. 組內平方誤差和——拐點圖

想必之前動輒幾十個指標，這裡就用一個最簡單的指標——sum of squared error (SSE)組內平方誤差和來確定最佳聚類數目。這個方法也是出於《R語言實戰》，自定義的一個求組內誤差平方和的函數。

wssplot <- function(data, nc=15, seed=1234){
   wss <- (nrow(data)-1)*sum(apply(data,2,var))
   for (i in 2:nc){
       set.seed(seed)
       wss[i] <- sum(kmeans(data, centers=i)$withinss)
       }
   plot(1:nc, wss, type="b", xlab="Number of Clusters",
       ylab="Within groups sum of squares")}

wssplot(dataset)

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/39ae40d81ccf72c41b48a28e609d661e.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;533&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;531&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;533&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic2.zhimg.com/v2-25b396108e9b5da6094c2097888f2251_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;&amp;lt;a href=&amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/39ae40d81ccf72c41b48a28e609d661e.png&amp;amp;quot;/&amp;gt;&amp;quot; _src=&amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/39ae40d81ccf72c41b48a28e609d661e.png&amp;amp;quot;/&amp;gt;&amp;quot; _href=&amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/39ae40d81ccf72c41b48a28e609d661e.png&amp;amp;quot;/&amp;gt;&amp;quot;&amp;gt;https://static.leiphone.com/uploads/new/article/pic/201710/39ae40d81ccf72c41b48a28e609d661e.png&amp;amp;amp;quot;/&amp;amp;amp;gt;&amp;lt;/a&amp;gt;

隨著聚類數目增多，每一個類別中數量越來越少，距離越來越近，因此WSS值肯定是隨著聚類數目增多而減少的，所以關注的是斜率的變化，但WWS減少得很緩慢時，就認為進一步增大聚類數效果也並不能增強，存在得這個「肘點」就是最佳聚類數目，從一類到三類下降得很快，之後下降得很慢，所以最佳聚類個數選為三

另外也有現成的包(factoextra)可以調用

library(factoextra)
library(ggplot2)
set.seed(1234)
fviz_nbclust(dataset, kmeans, method = "wss") +
   geom_vline(xintercept = 3, linetype = 2)

Loading required package: ggplot2

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/45de0c17eee70394a867f5889f7d3098.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;537&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;534&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;537&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic3.zhimg.com/v2-b80ac3d187e619c5fdaf46a3d3c9361e_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;&amp;lt;a href=&amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/45de0c17eee70394a867f5889f7d3098.png&amp;amp;quot;/&amp;gt;&amp;quot; _src=&amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/45de0c17eee70394a867f5889f7d3098.png&amp;amp;quot;/&amp;gt;&amp;quot; _href=&amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/45de0c17eee70394a867f5889f7d3098.png&amp;amp;quot;/&amp;gt;&amp;quot;&amp;gt;https://static.leiphone.com/uploads/new/article/pic/201710/45de0c17eee70394a867f5889f7d3098.png&amp;amp;amp;quot;/&amp;amp;amp;gt;&amp;lt;/a&amp;gt;

選定為3類為最佳聚類數目
用該包下的fviz_cluster函數可視化一下聚類結果

km.res <- kmeans(dataset,3)
fviz_cluster(km.res, data = dataset)

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/4c6ed79b2deec361f41a79be29eca253.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;532&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;530&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;532&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic2.zhimg.com/v2-9ff52d6a4c1dfeedff83fa03c895e2b9_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/4c6ed79b2deec361f41a79be29eca253.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

4. PAM(Partitioning Around Medoids) 圍繞中心點的分割算法

k-means算法取得是均值，那麼對於異常點其實對其的影響非常大，很可能這種孤立的點就聚為一類，一個改進的方法就是PAM算法，也叫k-medoids clustering
首先通過fpc包中的pamk函數得到最佳聚類數目

library(fpc)
pamk.best <- pamk(dataset)
pamk.best$nc

3

pamk函數不需要提供聚類數目，也會直接自動計算出最佳聚類數，這裡也得到為3
得到聚類數提供給cluster包下的pam函數並進行可視化

library(cluster)
clusplot(pam(dataset, pamk.best$nc))

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/3399bd76503f61d284bdccca0023d899.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;530&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;532&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;530&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic3.zhimg.com/v2-387e44d991195f07e55ffbdeb0c479be_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/3399bd76503f61d284bdccca0023d899.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

5.Calinsky criterion

這個評估標準定義[5]如下：
其中，k是聚類數，N是樣本數，SSw是我們之前提到過的組內平方和誤差， SSb是組與組之間的平方和誤差，SSw越小，SSb越大聚類效果越好，所以Calinsky criterion值一般來說是越大，聚類效果越好

library(vegan)
ca_clust <- cascadeKM(dataset, 1, 10, iter = 1000)
ca_clust$results

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/a12ebaef7daa5788c09e5eb3ba08c6d0.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;592&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;142&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;592&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic4.zhimg.com/v2-aa28be2976a4f350d80506c45026df97_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/a12ebaef7daa5788c09e5eb3ba08c6d0.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

可以看到該函數把組內平方和誤差和Calinsky都計算出來了，可以看到calinski在聚類數為3時達到最大值。

calinski.best <- as.numeric(which.max(ca_clust$results[2,]))
calinski.best

3

畫圖出來觀察一下

plot(fit, sortg = TRUE, grpmts.plot = TRUE)

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/699416f926e98813797da0ae31c45ce4.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;534&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;533&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;534&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic4.zhimg.com/v2-fc2cda58b5ba4856c681a731bc4906eb_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/699416f926e98813797da0ae31c45ce4.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

注意到那個紅點就是對應的最大值，自帶的繪圖橫軸縱軸取的可能不符合我們的直覺，把數據取出來自己單獨畫一下

calinski<-as.data.frame(ca_clust$results[2,])
calinski$cluster <- c(1:10)
library(ggplot2)
ggplot(calinski,aes(x = calinski[,2], y = calinski[,1]))+geom_line()

Warning message:
"Removed 1 rows containing missing values (geom_path)."

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/d1a5ca803f794122475c83a0f2ca700a.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;533&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;533&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;533&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic2.zhimg.com/v2-7655fc4a181aa14c3614f060e7981575_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/d1a5ca803f794122475c83a0f2ca700a.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

這個看上去直觀多了。這就很清晰的可以看到在聚類數目為3時，calinski指標達到了最大值，所以最佳數目為3

6.Affinity propagation (AP) clustering

這個本質上是類似kmeans或者層次聚類一樣，是一種聚類方法，因為不需要像kmeans一樣提供聚類數，會自動算出最佳聚類數，因此也放到這裡作為一種計算最佳聚類數目的方法。
AP算法的基本思想是將全部樣本看作網絡的節點，然後通過網絡中各條邊的消息傳遞計算出各樣本的聚類中心。聚類過程中，共有兩種消息在各節點間傳遞，分別是吸引度( responsibility)和歸屬度(availability) 。AP算法通過迭代過程不斷更新每一個點的吸引度和歸屬度值，直到產生m個高質量的Exemplar（類似於質心），同時將其餘的數據點分配到相應的聚類中[7]

library(apcluster)
ap_clust <- apcluster(negDistMat(r=2), dataset)
length(ap_clust@clusters)

15

該聚類方法推薦的最佳聚類數目為15，再用熱力圖可視化一下

heatmap(ap_clust)

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/dbc870375b6d45be952ee72e882fbbc0.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;527&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;530&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;527&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic3.zhimg.com/v2-b982fd6c85860755226f4ae347397252_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/dbc870375b6d45be952ee72e882fbbc0.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

選x或者y方向看(對稱)，可以數出來「葉子節點」一共15個

7. 輪廓係數Average silhouette method

輪廓係數是類的密集與分散程度的評價指標。

a(i)是測量組內的相似度,b(i)是測量組間的相似度，s(i)範圍從-1到1，值越大說明組內吻合越高，組間距離越遠——也就是說，輪廓係數值越大，聚類效果越好[9]

require(cluster)
library(factoextra)
fviz_nbclust(dataset, kmeans, method = "silhouette")

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/320ffea3de355eba9f40c34a92397420.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;535&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;531&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;535&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic3.zhimg.com/v2-f24229d24e86ded6219e4636e9849292_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/320ffea3de355eba9f40c34a92397420.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

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

8. Gap Statistic

之前我們提到了WSSE組內平方和誤差，該種方法是通過找「肘點」來找到最佳聚類數，肘點的選擇並不是那麼清晰，因此史丹福大學的Robert等教授提出了Gap Statistic方法，定義的Gap值為[9]

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/c62dd70be24d437c31a08065d289ac33.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;300&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;72&amp;amp;amp;quot; class=&amp;amp;amp;quot;content_image&amp;amp;amp;quot; width=&amp;amp;amp;quot;300&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/c62dd70be24d437c31a08065d289ac33.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

取對數的原因是因為Wk的值可能很大
通過這個式子來找出Wk跌落最快的點，Gap最大值對應的k值就是最佳聚類數

library(cluster)
set.seed(123)
gap_clust <- clusGap(dataset, kmeans, 10, B = 500, verbose = interactive())
gap_clust

Clustering Gap statistic ["clusGap"] from call:
clusGap(x = dataset, FUNcluster = kmeans, K.max = 10, B = 500,     verbose = interactive())
B=500 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
--> Number of clusters (method 'firstSEmax', SE.factor=1): 3
         logW   E.logW       gap     SE.sim
[1,] 5.377557 5.863690 0.4861333 0.01273873
[2,] 5.203502 5.758276 0.5547745 0.01420766
[3,] 5.066921 5.697322 0.6304006 0.01278909
[4,] 5.023936 5.651618 0.6276814 0.01243239
[5,] 4.993720 5.615174 0.6214536 0.01251765
[6,] 4.962933 5.584564 0.6216311 0.01165595
[7,] 4.943241 5.556310 0.6130690 0.01181831
[8,] 4.915582 5.531834 0.6162518 0.01139207
[9,] 4.881449 5.508514 0.6270646 0.01169532
[10,] 4.855837 5.487005 0.6311683 0.01198264

library(factoextra)
fviz_gap_stat(gap_clust)

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/dd77400b222298012bc677bcc05cb7a8.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;532&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;534&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;532&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic3.zhimg.com/v2-64bd9c3b553db08a095966397e88e26e_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/dd77400b222298012bc677bcc05cb7a8.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

可以看到也是在聚類數為3的時候gap值取到了最大值，所以最佳聚類數為3

9.層次聚類

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

h_dist <- dist(as.matrix(dataset))
h_clust<-hclust(h_dist)
plot(h_clust, hang = -1, labels = FALSE)
rect.hclust(h_clust,3)

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/ff45e8a530a86955edfe5cbe6478c20d.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;529&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;532&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;529&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic1.zhimg.com/v2-3d1b5ba7f76ce0ff1240ca497a613440_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/ff45e8a530a86955edfe5cbe6478c20d.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

10.clustergram

最後一種算法是Tal Galili[10]大牛自己定義的一種聚類可視化的展示，繪製隨著聚類數目的增加，所有成員是如何分配到各個類別的。該代碼沒有被製作成R包，可以去Galili介紹頁面)裡面的github地址找到原始碼跑一遍然後就可以用這個函數了，因為原始碼有點長我就不放博客裡面了，直接放出運行代碼的截圖。

clustergram(dataset, k.range = 2:8, line.width = 0.004)

Loading required package: colorspace
Loading required package: plyr

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/f31e6e16ca2b9f49e4af337493eb703a.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;527&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;530&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;527&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic4.zhimg.com/v2-9797f9f7e0288fdf2d034531db315663_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/f31e6e16ca2b9f49e4af337493eb703a.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

隨著K的增加，從最開始的兩類到最後的八類，圖肯定是越到後面越密集。通過這個圖判斷最佳聚類數目的方法應該是看隨著K每增加1，分出來的線越少說明在該k值下越穩定。比如k=7到k=8，假設k=7是很好的聚類數，那分成8類時應該可能只是某一類分成了兩類，其他6類都每怎麼變。反應到圖中應該是有6簇平行線，有一簇分成了兩股，而現在可以看到從7到8，線完全亂了，說明k=7時效果並不好。按照這個分析，k=3到k=4時，第一股和第三股幾本沒變，就第二股拆成了2類，所以k=3是最佳聚類數目

方法匯總與比較

wine數據集我們知道其實是分為3類的，以上10種判定方法中：

層次聚類和clustergram方法、肘點圖法，需要人工判定，雖然可以得出大致的最佳聚類數，但算法本身不會給出最佳聚類數

除了Affinity propagation (AP) clustering 給出最佳聚類數為15，剩下6種全都是給出最佳聚類數為3

選用上次文本挖掘的矩陣進行分析(667*1623)

可見上述方法中有的因為數據太大不能運行，有的結果很明顯不對，一個可能是數據集的本身的原因（缺失值太多等），但是也告訴了我們在確定最佳聚類數目的時候需要多嘗試幾種方法，並沒有固定的套路，然後選擇一種可信度較高的聚類數目。
最後再把這10種方法總結一下：

&amp;amp;amp;lt;img src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/e62c8e97de526716f340c418282308ab.png&amp;amp;amp;quot; data-rawwidth=&amp;amp;amp;quot;534&amp;amp;amp;quot; data-rawheight=&amp;amp;amp;quot;532&amp;amp;amp;quot; class=&amp;amp;amp;quot;origin_image zh-lightbox-thumb&amp;amp;amp;quot; width=&amp;amp;amp;quot;534&amp;amp;amp;quot; data-original=&amp;amp;amp;quot;https://pic3.zhimg.com/v2-d43a3c6a70d48edb2add13722f224ac2_r.png&amp;amp;amp;quot; _src=&amp;amp;amp;quot;https://static.leiphone.com/uploads/new/article/pic/201710/e62c8e97de526716f340c418282308ab.png&amp;amp;amp;quot;/&amp;amp;amp;gt;

參考文獻
[1]R語言實戰第二版
[2]Partitioning cluster analysis: Quick start guide - Unsupervised Machine Learning
[3]BIC：http://www.stat.washington.edu/raftery/Research/PDF/fraley1998.pdf
[4]Cluster analysis in R: determine the optimal number of clusters
[5]Calinski-Harabasz Criterion：Calinski-Harabasz criterion clustering evaluation object
[6]Determining the optimal number of clusters: 3 must known methods - Unsupervised Machine Learning
[7] affinity-propagation：聚類算法Affinity Propagation(AP)
[8]輪廓係數https://en.wikipedia.org/wiki/Silhouette(clustering))
[9]gap statistic-Tibshirani R, Walther G, Hastie T. Estimating the number of clusters in a data set via the gap statistic[J]. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2001, 63(2): 411-423.

[10]ClustergramsClustergram: visualization and diagnostics for cluster analysis (R code)

雷鋒網版權文章，未經授權禁止轉載。詳情見轉載須知。

相關焦點

  • Python機器學習10:機器學習中的六種分類算法及實現(上)
    在機器學習中，可以使用多種算法模型解決同一個問題，那麼如何從中選擇出最佳的算法模型呢？當然，這個問題沒有一種固定的答案，需要根據不同的問題，嘗試使用多種機器學習方法，比較各種算法模型在該問題上的效果，最終才能決定究竟選擇哪一種模型。
  • 分享最適合新手入門的10種機器學習算法
    最常見的機器學習類型是學習映射Y=f（X），用它來預測Y的值。這被稱為預測建模或預測分析，我們的目標是做出最準確的預測。 對於想了解機器學習基礎知識的新手，以下是數據科學家最常用的10種機器學習算法。 1.線性回歸 線性回歸也許是數據科學和機器學習中最知名、最好理解的算法了吧。
  • 每個數據科學家都應該知道的10種機器學習方法
    作者 | CDA數據分析師 10 machine learning methods that every data scientist should know機器學習是研究和工業中的熱門話題，新方法一直在發展。該領域的速度和複雜性使得即使對於專家而言也難以跟上新技術 - 並且對於初學者而言可能是壓倒性的難度。
  • 入門| 機器學習新手必看10大算法
    最常見的機器學習算法是學習映射 Y = f(X) 來預測新 X 的 Y。這叫做預測建模或預測分析，我們的目標是儘可能作出最準確的預測。 對於想了解機器學習基礎知識的新手，本文將概述數據科學家使用的 top 10 機器學習算法。 1.
  • 機器學習和深度學習最佳作品
    來源：CPDA數據分析師網 機器學習算法的超參數優化：理論與實踐 機器學習算法已廣泛應用於各種應用程式和領域，為了使機器學習模型適應不同的問題，必須調整其超參數，為機器學習模型選擇最佳的超參數配置會直接影響模型的性能，它通常需要對機器學習算法和適當的超參數優化技術有深入的了解
  • 入門機器學習之線性回歸
    什麼是回歸分析在客觀世界中普通存在著變量之間的關係，變量之間的關係一般來說可以分成確定性關係和不確定關係，確定性關係是說變量之間的關係是可以用函數關係來表示的，另一種不確定性關係即所謂相關關係。這也是處理高維數據集的方法之一。
  • 數據科學和機器學習的最佳Python庫
    數據科學和機器學習是該時代最需求的技術，這一需求促使每個人都學習不同的庫和軟體包以實現它們。這篇博客文章將重點介紹用於數據科學和機器學習的Python庫。這些是您掌握市場上最被炒作的兩項技能的庫。更少的代碼：實施數據科學和機器學習涉及無數的算法。得益於Python對預定義包的支持，我們不必編寫算法。為了使事情變得更容易，Python提供了一種「在編碼時檢查」的方法，從而減輕了測試代碼的負擔。
  • 機器學習頂會ICML 2018:復旦大學副教授獲最佳論文亞軍 騰訊清華...
    智東西(公眾號:zhidxcom)文 | 心緣智東西7月12日消息，第35屆機器學習國際大會ICML 2018在7月10日至15日期間登陸瑞典斯特哥爾摩。ICML官網上提前公布了最佳論文名單，來自MIT和UC Berkeley的研究人員摘得最佳論文的桂冠。
  • 機器學習備忘錄 | AUC值的含義與計算方法
    內容導讀在機器學習領域， AUC 值經常用來評價一個二分類模型的訓練效果，對於許多機器學習或者數據挖掘的從業者或在校學生來說， AUC 值的概念也經常被提起，但由於大多數時候我們都是藉助一些軟體包來訓練機器學習模型，模型評價指標的計算往往被軟體包所封裝，因此我們常常會忽略了它們的具體意義，這在有些時候會讓我們對於手頭正在進行的任務感到困惑。
  • 五分鐘了解機器學習十大算法
    今天，我們將簡要介紹 10 種最流行的機器學習算法，這樣你就可以適應這個激動人心的機器學習世界了！讓我們言歸正傳！1. 線性回歸線性回歸（Linear Regression）可能是最流行的機器學習算法。線性回歸就是要找一條直線，並且讓這條直線儘可能地擬合散點圖中的數據點。它試圖通過將直線方程與該數據擬合來表示自變量（x 值）和數值結果（y 值）。
  • 謝國彤:疾病預測的機器學習、深度學習和經典回歸方法
    數據不均衡：很多疾病的發病率都不高，比如房顫患者發生腦卒中的平均概率是 10%，腦卒中患者出院後導致殘疾的平均概率是 4%。造成數據中正例相對較少，很不均衡，對機器學習算法的要求更高。疾病預測的主要方法疾病預測的主要方法可以簡單的分為經典回歸方法、機器學習方法和深度學習方法三大類。下面分別用三篇論文舉例介紹一下。
  • 機器學習的統計方法 貝葉斯決策理論入門
    引言無論你是在建立機器學習模型還是在日常生活中做決定，我們總是選擇風險最小的方案。作為人類，我們天生就採取任何有助於我們生存的行動；然而，機器學習模型最初並不是基於這種理解而建立的。這些算法需要經過訓練和優化，以選擇風險最小的最優方案。此外，很重要的一點在於，我們必須明白，如果某些高風險的決定做的不正確，將會導致嚴重的後果。我們以癌症診斷為例。
  • 路基土的最大幹密度和最佳含水量確定方法
    【學員問題】路基土的最大幹密度和最佳含水量確定方法？　　【解答】路基受到的荷載應力，隨深度而迅速減少，所以路基上部的壓實度應高一些；另外，道路等級高，其路面等級也高，對路基強度的要求則相應提高，所以對路基壓實度的要求也應高一些。
  • 物理所等利用機器學習方法預測材料性能獲進展
    近二十年來，機器學習方法的發展為我們的生活帶來許多便利。智能網絡搜索、語音識別，乃至無人超市、無人駕駛汽車等，依託於機器學習方法的新事物正迅速地在生活中普及。Alpha Go的橫空出世更讓世界驚嘆於人工智慧的潛在價值。在科研領域，大數據的理念正在改變著科研人員對未知世界的探索方式。
  • 機器學習的創新之舉：解答未知
    通過機器學習，研究人員可以多次測試同一假設以找到最佳答案。越來越多的生物學家開始使用計算模型來解釋複雜系統，比如：在神經科學領域，研究人員正在嘗試調試用於預測天氣或過濾垃圾郵件的算法，以了解大腦神經網絡的信息處理機制。
  • 解答未知——機器學習的創新之舉
    通過機器學習，研究人員可以多次測試同一假設以找到最佳答案。越來越多的生物學家開始使用計算模型來解釋複雜系統，比如：在神經科學領域，研究人員正在嘗試調試用於預測天氣或過濾垃圾郵件的算法，以了解大腦神經網絡的信息處理機制。
  • 分享一個Python中機器學習的特徵選擇工具
    我對臨時的特徵選擇方法感到很失望，但是在解決機器學習問題時又反覆用到了這些方法，所以就創建了一個關於特徵選擇的Python類，該類可以在GitHub上找到。　　特徵選擇器有5種查找要刪除特徵的方法，找到被標識的特徵並手動地將它們從數據中刪除，或者使用特徵選擇器中的remove函數。
  • 僅需10分鐘:開啟你的機器學習之路
    有幾種方法可以實現這種技術。在很多時候這都是常識。我們以登船地數據為例——這是用 Q、S 或 C 填充的數據。Python 庫不能處理這個，因為它只能處理數字。所以你需要用所謂的獨熱向量化（One Hot Vectorization）來處理，它可以把一列變成三列。
  • 量子機器學習入門科普:解讀量子力學和機器學習的共生關係
    雖然最終能找到一個結果，但這種方法相當耗時。對比之下，量子計算機「解鎖了神奇的並行性。它們同時探尋玉米迷宮中的每一條路。」因此，量子計算機可能指數級減少解決問題的步驟。這種並行性正是起源於量子物理中「量子位（qubit）」、「疊加（superposition）」和「糾纏（entanglement）」等理論。這其中的神奇之處遠不止如此，還得繼續往下看。
  • 機器學習最佳統計書籍推薦
    幾乎每一個機器學習項目中都離不開統計方法。所以需要基礎的統計相關基礎知識和方法對機器學習和AI必不可少。我們說AI不是神話，AI是數學算法，說明紮實數學基礎是做AI的必須，而統計學知識尤其是如此。這意味著重要的是要牢固掌握統計學主要發現的基礎和相關統計方法的工作知識。