Rでk-means法とその拡張1 k-means編 - サボタージュ禁止のおさぼり日記

あっ, どーも僕です。

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

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

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)