R語言-kmeans聚類理論篇K的選擇（輪廓係數）
kmeans是最簡單的聚類演算法之一，但是運用十分廣泛。最近在工作中也經常遇到這個演算法。kmeans一般在資料分析前期使用，選取適當的k，將資料分類後，然後分類研究不同聚類下資料的特點。
本文記錄學習kmeans演算法相關的內容，包括演算法原理，收斂性，效果評估聚，最後帶上R語言的例子，作為備忘。
演算法原理
kmeans的計算方法如下：
1 隨機選取k箇中心點
2 遍歷所有資料，將每個資料劃分到最近的中心點中
3 計算每個聚類的平均值，並作為新的中心點
4 重複2-3，直到這k箇中線點不再變化（收斂了），或執行了足夠多的迭代
時間複雜度：O(I*n*k*m)
空間複雜度：O(n*m)
其中m為每個元素欄位個數，n為資料量，I為跌打個數。一般I,k,m均可認為是常量，所以時間和空間複雜度可以簡化為O(n)，即線性的。
演算法收斂
從kmeans的演算法可以發現，SSE其實是一個嚴格的座標下降（Coordinate Decendet）過程。設目標函式SSE如下：
SSE(,
,…,
)
=
採用歐式距離作為變數之間的聚類函式。每次朝一個變數的方向找到最優解，也就是求偏倒數，然後等於0，可得
c_i= 其中m是c_i所在的簇的元素的個數
也就是當前聚類的均值就是當前方向的最優解（最小值），這與kmeans的每一次迭代過程一樣。所以，這樣保證SSE每一次迭代時，都會減小，最終使SSE收斂。
由於SSE是一個非凸函式（non-convex function），所以SSE不能保證找到全域性最優解，只能確保區域性最優解。但是可以重複執行幾次kmeans，選取SSE最小的一次作為最終的聚類結果。
0-1規格化
由於資料之間量綱的不相同，不方便比較。舉個例子，比如遊戲使用者的線上時長和活躍天數，前者單位是秒，數值一般都是幾千，而後者單位是天，數值一般在個位或十位，如果用這兩個變數來表徵使用者的活躍情況，顯然活躍天數的作用基本上可以忽略。所以，需要將資料統一放到0~1的範圍，將其轉化為無量綱的純數值，便於不同單位或量級的指標能夠進行比較和加權。具體計算方法如下：
其中屬於A。
輪廓係數
輪廓係數（Silhouette Coefficient）結合了聚類的凝聚度（Cohesion）和分離度（Separation），用於評估聚類的效果。該值處於-1~1之間，值越大，表示聚類效果越好。具體計算方法如下：
-
對於第i個元素x_i，計算x_i與其同一個簇內的所有其他元素距離的平均值，記作a_i，用於量化簇內的凝聚度。
-
選取x_i外的一個簇b，計算x_i與b中所有點的平均距離，遍歷所有其他簇，找到最近的這個平均距離,記作b_i，用於量化簇之間分離度。
-
對於元素x_i，輪廓係數s_i = (b_i – a_i)/max(a_i,b_i)
-
計算所有x的輪廓係數，求出平均值即為當前聚類的整體輪廓係數
從上面的公式，不難發現若s_i小於0，說明x_i與其簇內元素的平均距離小於最近的其他簇，表示聚類效果不好。如果a_i趨於0，或者b_i足夠大，那麼s_i趨近與1，說明聚類效果比較好。
K值選取
在實際應用中，由於Kmean一般作為資料預處理，或者用於輔助分類貼標籤。所以k一般不會設定很大。可以通過列舉，令k從2到一個固定值如10，在每個k值上重複執行數次kmeans(避免區域性最優解)，並計算當前k的平均輪廓係數，最後選取輪廓係數最大的值對應的k作為最終的叢集數目。
實際應用
下面通過例子（R實現，完整程式碼見附件）講解kmeans使用方法，會將上面提到的內容全部串起來
| 1 2 3 |
library(fpc)# install.packages("fpc")
data(iris)
head(iris)
|
載入實驗資料iris，這個資料在機器學習領域使用比較頻繁，主要是通過畫的幾個部分的大小，對花的品種分類，實驗中需要使用fpc庫估計輪廓係數，如果沒有可以通過install.packages安裝。
| 1 2 3 4 5 6 7 8 9 |
# 0-1 正規化資料
min.max.norm <-function(x){
(x-min(x))/(max(x)-min(x))
}
raw.data <- iris[,1:4]
norm.data <- data.frame(sl = min.max.norm(raw.data[,1]),
sw = min.max.norm(raw.data[,2]),
pl = min.max.norm(raw.data[,3]),
pw = min.max.norm(raw.data[,4]))
|
對iris的4個feature做資料正規化，每個feature均是花的某個不為的尺寸。
| 1 2 3 4 5 6 7 8 9 10 11 12 13 |
# k取2到8，評估K
K <- 2:8
round <- 30# 每次迭代30次，避免區域性最優
rst <- sapply(K,function(i){
print(paste("K=",i))
mean(sapply(1:round,function(r){
print(paste("Round",r))
result <- kmeans(norm.data, i)
stats <- cluster.stats(dist(norm.data), result$cluster)
stats$avg.silwidth
}))
})
plot(K,rst,type='l',main='輪廓係數與K的關係',
ylab='輪廓係數')
|
評估k，由於一般K不會太大，太大了也不易於理解，所以遍歷K為2到8。由於kmeans具有一定隨機性，並不是每次都收斂到全域性最小，所以針對每一個k值，重複執行30次，取並計算輪廓係數，最終取平均作為最終評價標準，可以看到如下的示意圖，
當k取2時，有最大的輪廓係數，雖然實際上有3個種類。
| 1 2 3 4 5 6 7 8 |
# 降緯度觀察
old.par <- par(mfrow = c(1,2))
k = 2# 根據上面的評估 k=2最優
clu <- kmeans(norm.data,k)
mds = cmdscale(dist(norm.data,method="euclidean"))
plot(mds, col=clu$cluster, main='kmeans聚類 k=2',
pch = 19)
plot(mds, col=iris$Species, main='原始聚類',
pch = 19)
par(old.par)
|
聚類完成後，有源原始資料是4緯，無法視覺化，所以通過多維定標(Multidimensional scaling)將緯度將至2為，檢視聚類效果，如下
可以發現原始分類中和聚類中左邊那一簇的效果還是擬合的很好的，右測原始資料就連在一起，kmeans無法很好的區分，需要尋求其他方法。
kmeans最佳實踐
1. 隨機選取訓練資料中的k個點作為起始點
2. 當k值選定後，隨機計算n次，取得到最小開銷函式值的k作為最終聚類結果，避免隨機引起的區域性最優解
3. 手肘法選取k值：繪製出k--開銷函式閃點圖，看到有明顯拐點（如下）的地方，設為k值，可以結合輪廓係數。
4. k值有時候需要根據應用場景選取，而不能完全的依據評估引數選取。