読者です 読者をやめる 読者になる 読者になる

Rでk-means法とその拡張1 k-means編

R

あっ, どーも僕です。


先週, 「明日から頑張る」と宣言したはずの後輩が, その日から一度も研究室に来てくれません。
まだ, 卒論に書くことないがほとんどはずなのに.....

先輩としてどうしたらいいものか。

k-means法によるクラスリングの拡張

k-means法とその拡張が実装できたのでご紹介。

k-meansとは

k-means法とは, 非階層のハードクラスタリングで, 様々な本で紹介*1されています。
改めて説明する必要もないくらいによく出てきます。

特徴としては次のことがあげられるそうです。

  • 実装が簡単
  • 高速

アルゴリズムはなにか本を参照してください。
ちょっと物足りない点としては次があげられるそうです。

  • 最適なクラスター数が決める基準がない
  • 変数の相関を考慮できない

そんな, k-meansはいくつか拡張されたものがあります

  • Fuzzy c-means
  • x-means*2
  • 変数間の関係を考慮したk-means法(以下, 改良k-means)*3

Fuzzy c-meansはk-meansのソフトクラスタリング版です. Rで学ぶクラスタ解析で解説されています。パッケージでやるならe1071などに入っていますし, わたしも既にスクリプトをRで普通の数値計算に載せてあります。

x-meansではデータを二分割するk-meansを繰り返し, 再度結合することで, クラスタリングを行います。情報量基準(BIC)でクラスター数が決まり, そのクラスター数はアルゴリズムが自動で決めてくれます。石岡さんのHP(Open softwares and its documents)へいくと, Rのスクリプトがあります。パッケージではRWekaにありましたが, 改良される前のx-meansぽかったです。x-meansは自動で, クラスター数が決まるのですが, 同じデータを入れても同じ結果とならないです。

改良k-means法は, 一度k-meansでクラスタリングした後に, 混合正規分布で尤度を考えて, そのクラスタリングを補正します。こちらも, 同じデータを入れても同じ結果とはなりません。ただ, AICとBICが計算できるのでなんどか試して, 最適なものを選択します。この方法は, クラスターを自動で決めることはできませんが, クラスター数を変えながら情報量で最適なクラスター数を決めます。Rのパッケージはなさそうですが, 清水裕士さんが開発しているHADに実装されています。

k-meansの実装

で, 今回はx-meansと改良k-meansをRで実装してみました。

このとき, どちらもまずk-meansが必要となりますので, ここではk-meansの実装だけ載せておきますね。

ただしくできているかはirisで検証しました。その結果, Rのkmeansと同じ結果になりました。
実装したものとRのkmeansを使ったものの比較がこちら。
ちなみに色が同じになっていますが, 色はただクラスターの違いを表すだけですので, 違ってても問題ないです。
f:id:aaaazzzz036:20131127203005p:plain



スクリプトがこちら。
最後の方にある動作確認で結果が完全に一致することがわかります。

# # ***********************************************************
# 二乗ノルムを距離にしたK-平均法
# args
# - data : data.frame (data, variate)
# - - 変数がひとつのときは, iris[, 1, drop = FALSE]のように
# - - drop を入れることで, data.frame落ちを防げる
# - iter.max : 最大繰り返し回数
# - K : クラス数
# return 
# - cluster : 各データのクラス
# - centers : クラスごとの平均値ベクトル
# - iter    : 繰り返し計算数
# - radius  : 最大半径
# *************************************************************
myKmeans <- function (data, K, iter.max = 1000) {
    data <- na.omit (as.data.frame (data))
    nc   <- ncol (data)
    nr   <- nrow (data)
    if (any (sapply (data, class) != "numeric")) {
        return (0)
    }
    K    <- as.integer (K)
    if (K <= 1) {
        return (list (cluster = rep (1, nr)))
    }
    if (nr < K) {
        return (list (cluster = rep (1, nr)))
    }
    if (nc == 1) {
        # 1変数だとcolMeansが動かない
        Mean <- function (x) {mean (unlist (x))}
    } else {
        Mean <- colMeans
    }
    
    # 適当に初期クラスを振る
    # そのうちk-means++を実装する
    doInitialize <- TRUE
    while (doInitialize) {
        cluster.old  <- sample (1:K, nr, replace = TRUE)
        doInitialize <- length (unique (cluster.old)) != K
    }
    cluster.centers <- by (data, cluster.old, Mean)

    
    # 収束計算
    doNEXT   <- TRUE
    iter     <- 1
    while (doNEXT) {
        iter <- iter + 1
        if (iter > iter.max) break
        cluster.dists <- sapply (1:K, function (k) {
            cluster.center <- cluster.centers[[k]]
            data.shifted   <- sweep (data, 2, cluster.center, "-")
            rowSums (data.shifted * data.shifted)
        })
        # 距離が最小のものを選択
        cluster.new     <- max.col (- cluster.dists)
        doNEXT          <- sum (abs (cluster.new - cluster.old)) != 0
        cluster.old     <- cluster.new
        if (length (unique(cluster.old)) != K) {
            # クラスター数が落ちてしまったら振り直す
            doInitialize <- TRUE
            while (doInitialize) {
                cluster.old  <- sample (1:K, nr, replace = TRUE)
                doInitialize <- length (unique (cluster.old)) != K
            }
        }
        cluster.centers <- by (data, cluster.old, Mean)
    }
    names (cluster.new) <- NULL
    cluster.radius <- sapply (1:K, function (k) {
        data.k         <- data[cluster.new == k, ]
        data.k.shifted <- sweep (data.k, 2, cluster.centers[[k]], "-")
        sqrt (max (rowSums (data.k.shifted * data.k.shifted)))
    })
    
    return (list (cluster = cluster.new, 
                  centers = cluster.centers, 
                  iter    = iter,
                  radius  = cluster.radius))
}

# 動作確認
K <- 2
myiris <- iris[, 3:4, drop = FALSE]
myKmeans(myiris, K = K)
print (t1 <- myKmeans (myiris, K = K))
print (t2 <- kmeans (myiris, K))
print (table (t1$cluster, t2$cluster))


# グラフで確認
if (K == 2) {
    col <- c ("skyblue", "yellowgreen")
} else {
    col <- brewer.pal (K, "Set2")
}
layout (matrix (c(1, 2), 1))
plot (myiris, asp = 1, cex = 1.3, col = col[t1$cluster], main = "myKmeans", cex.main = 1.5, pch = 19)
grid (lty = 5, col = "lightgrey")
center <- do.call (rbind, t1$centers)
symbols(center, circles=t1$radius, inches=FALSE, asp=1, col = "lightgrey", add = TRUE)
points (center, pch = c("X"), cex = 2)
plot (myiris, asp = 1, cex = 1.3,col = col[t2$cluster], main = "R Kmeans", cex.main = 1.5, pch = 19)
grid (lty = 5, col = "lightgrey")
points (t2$centers, pch = c("X"), cex = 2)


以上で, k-meansの実装が正しくできているのが確認できたので, 今回はここまで。

k-meansは実装が楽と言われていますが, クラスター数が落ちたりすることに気づかなくて, ちょっと時間がかかってしまいました。


なお, 今回を含めて3回でスクリプトを載せますが, Rで普通の数値解析のフォルダにスクリプトを置いときます。

他にもいままで載せていなかったスクリプトとして
正方行列Aのn乗を計算する関数(myPowerMatrix.r)やグラムシュミットの直交化(myOrthogonalize_Schmidt.r)とかこまいスクリプトも加えているんでよかったらみてください。

*1:例えば, 新納浩幸:Rで学ぶクラスタ解析, オーム社, 2007

*2:石岡恒憲:x-means法改良の一提案ーk-means法の逐次繰り返しとクラスターの再併合ー, 計算機統計学, 第18巻, 第1号, 2006, pp.3-13

*3:豊田秀樹, 池原一哉:変数間の関係性を考慮してクラスター数を決定するk-means法の改良, 心理学研究, 第82号, 第1号, 2011, pp.32-40