RでK-meansの最適なクラスタ数をAIC / BICに基づいて求める - 渋谷駅前で働くデータサイエンティストのブログ

これはただの備忘録です。既知の話題ばかりが並べられているので、特に新鮮味のない内容である点予めご容赦ください。

クラスタリング手法として広く知られるK-meansは、その簡便さから非常に広汎に使われていますが、一方で「クラスタ数を恣意的に決め打ちせざるを得ない」という難点があり、「最適なクラスタ数をどうやって決めるか」という課題が長年に渡ってあります。この課題の解決策についてちょっと調べてみたので、以下にまとめてみました。

K-meansにおける「最適なクラスタ数の決め方」として、こちらの記事では伝統的な手法ということで

エルボー法
シルエット分析
X-means（K-meansに情報量規準を適用して再帰的に最適クラスタ数を決める）

の3種類が紹介されています。これらは僕も以前から聞いたことがあるもので、実際K-meansの実装の中には最初からエルボー法などを含んでいるものもあったりします。またそれぞれの方法で参照する規準としてDavies-Bouldin indexとかDunn indexといったものがありますが、「そもそもAIC / BICとか使えないんだろうか？」と思ったのでした。

そこで調べてみたところ、AIC / BICを使うやり方が実際にありました。

ちなみに「本当にfit$tot.withinssからAIC / BICを求めて大丈夫なのか？」という質疑がCross Validatedにあって、「データセットを平均0＆分散1に正規化してあればOK」という説明がされています。

fit$tot.withinssが尤度の代わりに使えるかどうかについては、以下のコメントがCross Validatedにあるところを見ると問題なさそうです。

f:id:TJO:20210415221957p:plain

ということで、UCI ML repositoryの"Seeds"データセットで簡単に試してみます。AICだと推定される最適なKが大きめになるということが幾つかの資料に書かれているので、AIC / BICとも算出しますが主にBICの方を参照することとします。

kmeansAIC <- function(fit){
  
  m <- ncol(fit$centers)
  n <- length(fit$cluster)
  k <- nrow(fit$centers)
  D <- fit$tot.withinss
  return(data.frame(AIC = D + 2*m*k,
                    BIC = D + log(n)*m*k))
}

d <- read.csv('https://archive.ics.uci.edu/ml/machine-learning-databases/00236/seeds_dataset.txt', header = F, sep = '\t')
d <- na.omit(d)
d <- scale(d)

res <- c()
for (i in 2:15){
  fit <- kmeans(d, i)
  tmp <- kmeansAIC(fit)
  res <- rbind(res, cbind(tmp, i))
}
print(res)

#R>         AIC      BIC  i
#R> 1  860.3519 913.0448  2
#R> 2  472.5755 551.6148  3
#R> 3  450.6845 556.0703  4
#R> 4  413.8463 545.5785  5
#R> 5  409.4394 567.5180  6
#R> 6  404.5625 588.9875  7
#R> 7  438.1623 648.9338  8
#R> 8  383.7316 620.8495  9
#R> 9  425.3442 688.8086 10
#R> 10 383.4830 673.2938 11
#R> 11 375.5858 691.7431 12
#R> 12 396.1273 738.6310 13
#R> 13 429.2771 798.1272 14
#R> 14 402.3876 797.5842 15

BICだとk = 5がベストらしい、という結果になりました。念のため、Ward法で階層的クラスタリングした結果のdendrogramを描いてみて、本当にこれが妥当かどうかを確かめてみます。

plot(hclust(dist(d), method = "ward.D2"))

f:id:TJO:20210401175307p:plain

これを見る限りではBICが2番目に小さくなるk = 3が妥当なんじゃないかなぁ……ということで、Rのstats::kmeansでもAIC / BICによる最適クラスタ数の推定（らしきこと）ができるということが分かりました。

追記

コメント欄でリクエストがありましたので、エルボー法も試してみました。

#Elbow Method for finding the optimal number of clusters
set.seed(123)
# Compute and plot wss for k = 2 to k = 15.
k.max <- 15
wss <- sapply(1:k.max, 
              function(k){kmeans(d, k, nstart = 50, iter.max = 15)$tot.withinss})
print(wss)
plot(1:k.max, wss,
     type="b", pch = 19, frame = FALSE, 
     xlab="Number of clusters K",
     ylab="Total within-clusters sum of squares")