【機器學習】肯定最佳聚類數目的10種方法

時間 2019-11-14

標籤機器學習肯定最佳數目方法简体版

原文原文鏈接

【機器學習】肯定最佳聚類數目的10種方法

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

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

Loading required package: cluster

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

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

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

斷定方法

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種類型了。數目分別爲6三、5一、64。再畫出14個指標隨着聚類數目變化的走勢圖網絡

plot(m_clust, "BIC")

下表是這些模型的意義app

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

注意：此BIC並非貝葉斯信息準則！！！ide

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

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

與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 = "支持指標數")

能夠看到有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)

隨着聚類數目增多，每個類別中數量愈來愈少，距離愈來愈近，所以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

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

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

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

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

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

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

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

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

能夠看到該函數把組內平方和偏差和Calinsky都計算出來了，能夠看到calinski在聚類數爲3時達到最大值。

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

畫圖出來觀察一下

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

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

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)."

這個看上去直觀多了。這就很清晰的能夠看到在聚類數目爲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，再用熱力圖可視化一下

heatmap(ap_clust)

選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")

能夠看到也是在聚類數爲3時輪廓係數達到了峯值，因此最佳聚類數爲3

8. Gap Statistic

以前咱們提到了WSSE組內平方和偏差，該種方法是經過找「肘點」來找到最佳聚類數，肘點的選擇並非那麼清晰，所以斯坦福大學的Robert等教授提出了Gap Statistic方法，定義的Gap值爲[9]

取對數的緣由是由於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)

能夠看到也是在聚類數爲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)

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

隨着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)

mclust效果不好，14種模型只有6種有結果
bclust報錯
SSE能夠運行
fpc包中的pamk函數聚成2類，明顯不行
Calinsky criterion聚成2類
Affinity propagation (AP) clustering 聚成28類，相對靠譜
輪廓係數Average silhouette聚類2類
gap-Statistic跑不出結果

可見上述方法中有的由於數據太大不能運行，有的結果很明顯不對，一個多是數據集的自己的緣由（缺失值太多等）可是也告訴了咱們在肯定最佳聚類數目的時候須要多嘗試幾種方法，並無固定的套路，而後選擇一種可信度較高的聚類數目。
最後再把這10種方法總結一下：

參考文獻
[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)