X-means アルゴリズムを実装してみます。
参考文献
TLに以下の記事の引用ツイートが流れてきたのでやってみます。
スクリプトを全面的に参考にさせていただきました。
aaaazzzz036.hatenablog.com
使用データ
以前に以下の記事で、各都道府県の最低賃金、合計特殊出生率、高齢化率を用意したので流用します。
入門 機械学習による異常検知―Rによる実践ガイド: ノート3 - クッキーの日記
やること
最低賃金、合計特殊出生率、高齢化率の3変数に着目して都道府県をクラスタリングします。
k-means法はクラスタリングのアルゴリズムとしてよく知られていますが、このアルゴリズムではクラスタ数を所与としなければなりません。でも、今回はクラスタ数をどうすべきかよくわかりません。
1つの考え方として、まず適当なクラスタ数でk-means法を一旦適用してみて、次はそのクラスタにさらにk-means法を適用して、もしさらにk-means法を適用した細かいクラスタリングの方がデータをよりよく説明するならそちらを採用した方がよいと考えられます。データが重心の周りに正規分布すると仮定し、ベイズ情報量基準が小さくなる限りk-means法を繰り返し適用するのがX-means法です。
Determining the number of clusters in a data set - Wikipedia
手順
以下のように最低賃金、合計特殊出生率、高齢化率を data.txt にまとめてから読み込みます。
x <- read.table("data.txt", sep=",", header=FALSE) colnames(x) <- c("pref", "region", "wage", "birth", "aging")
3次元プロットしてみると以下です。なんかよくわかりません。
library(scatterplot3d) windowsFonts(MEI=windowsFont("Meiryo")) par(family="MEI", las=2) scatterplot3d(x[,3:5], pch=19)
par(mgp=c(2.7, 0.7, 0)) par(mar=c(4,4,1,1)) layout(matrix(c(1,0,2,3), 2, 2, byrow=TRUE)) plot(x$wage, x$birth, family="MEI", pch=19) plot(x$wage, x$aging, family="MEI", pch=19) plot(x$birth, x$aging, family="MEI", pch=19)
それで、このデータをX-means法でクラスタリングすると以下です。ただし前提条件はコードを参照してください。
x[3:5] <- scale(x[3:5]) # 最低賃金、合計特殊出生率、高齢化率をスケーリング source("xmeans.R") xkm <- xmeans$new() xkm$xmeans(x[3:5]) print(xkm$Cluster) plot(x[3:5], col=ifelse(xkm$Cluster==1, "salmon", "royalblue"), pch=19)
上の結果より、最低賃金、合計特殊出生率、高齢化率をスケーリングしてX-means法でクラスタリングすると都道府県は2クラスタに分割されることがわかります。グラフからわかる各クラスタの特徴は以下です。
最低賃金 | 合計特殊出生率 | 高齢化率 | |
---|---|---|---|
クラスタ1 | 高い | 低い | 低い |
クラスタ2 | 低い | 高い | 高い |
sudillap.hatenablog.com
地図データとレコードの順序を合わせるために、都道府県名をアルファベットにしておく必要があります。
x <- read.table("data.txt", sep=",", header=FALSE) colnames(x) <- c("pref", "region", "wage", "birth", "aging") x <- x[order(x$pref),] # 地図データとの整合性のため県名アルファベット順にソート x[3:5] <- scale(x[3:5]) source("xmeans.R") xkm <- xmeans$new() xkm$xmeans(x[3:5]) print(xkm$Cluster) ##### 日本地図のプロット ##### library(maptools) library(spsurvey) jpn <- readShapePoly("JPN_adm_shp/JPN_adm1.shp") jpn.dbf <- read.dbf("JPN_adm_shp/JPN_adm1.dbf") xlim <- c(126, 146) # 経度 ylim <- c(26, 46) # 緯度 par(mar=c(0,0,0,0)) library(RColorBrewer) plot(jpn, xlim=xlim, ylim=ylim, col=ifelse(xkm$Cluster==1, "salmon", "royalblue"))
上の散布図とクラスタが微妙に違います。元々k-means法が初期値依存なのと、レコード順序も変わったため。
上のスクリプトを何回も実行するとクラスタが1つになったり3つになったりもします。
スクリプト
ほぼ参考文献の記事( Rでk-means法とその拡張2 x-means編 - サボタージュ禁止のおさぼり日記 )のままです。
ただ、k-means法のクラスタ初期化はk-means++法で固定、細分したクラスタをマージするかどうかはするで固定、共分散行列の非対角成分(共分散)は無視で固定しています。元記事のスクリプトは引数で変更できるようになっています。
library(dplyr) kmeans <- setRefClass( Class="kmeans", methods=list( initCluster_kmeanspp=function(dat, NUM_CLUSTER) { ##### クラスタの初期化 (k-means++法による) ##### index <- 1:nrow(dat) center_index <- sample(index, 1) # データから無作為に1点を選ぶ (これが1クラスタ目の中心) center_ <- unlist(dat[center_index,]) min_dist <- rowSums((sweep(dat, 2, center_, "-"))**2) # 各データとcenter_index番目のデータとの2乗距離 while (length(center_index) < NUM_CLUSTER) { # 目的のクラスタ数になるまでデータ中心を追加 # 収集済みのどのデータ中心からも遠いデータほど (距離の2乗に比例した確率で) 選ばれやすいようにする prob <- min_dist[-center_index] / sum(min_dist[-center_index]) new_center <- sample(index[-center_index], 1, prob=prob) # データから1点を選ぶ center_ <- unlist(dat[new_center,]) new_dist <- rowSums((sweep(dat, 2, center_, "-"))**2) min_dist <- pmin(min_dist, new_dist) # 各データと収集済みデータ中心の中で1番近いものとの2乗距離 center_index <- c(center_index, new_center) } cur_center <- list() for (i in seq_len(NUM_CLUSTER)) { cur_center[[i]] <- unlist(dat[center_index[i],]) } dist <- sapply(seq_len(NUM_CLUSTER), function(k) { center_ <- cur_center[[k]] rowSums((sweep(dat, 2, center_, "-"))**2) }) cur_clusters <- max.col(-dist) return(list(cur_center=cur_center, cur_clusters=cur_clusters)) }, updateCenter=function(dat, cur_clusters) { ##### クラスタ分けに従ってデータ中心を更新 ##### cur_center <- cbind(dat, cc=cur_clusters) %>% group_by(cc) %>% do(rtn=colMeans (dplyr::select(., -cc))) %>% as.list() %>% "$"(rtn) return (cur_center) }, kmeans=function(dat_, K_=2) { ##### メイン処理 ##### dat <- na.omit(dat_[, sapply(dat_, is.numeric)]) NUM_CLUSTER <- as.integer(K_) SIZE <- as.integer(nrow(dat)) tmp <- .self$initCluster_kmeanspp(dat, NUM_CLUSTER) # クラスタ初期化 cur_clusters <- tmp$cur_clusters cur_center <- tmp$cur_center old_clusters <- cur_clusters convergence <- FALSE iter <- 1L while (!convergence) { iter <- iter + 1L if (iter > 100L) { cat("iter reached MAX_ITER\n") break } # 各データの所属クラスタを更新 (最短距離のデータ中心へ) distances <- sapply(seq_len(NUM_CLUSTER), function(k) { center_ <- cur_center[[k]] d_ <- sweep(dat, 2, center_, "-") return(rowSums(d_*d_)) }) cur_clusters <- max.col(-distances) if (length(unique(cur_clusters)) != NUM_CLUSTER) { # クラスタが消滅してしまった場合やり直し cur_clusters <- .self$initCluster_kmeanspp(dat, NUM_CLUSTER)$cur_clusters } convergence <- identical(old_clusters, cur_clusters) # 所属クラスタの更新がなければアルゴリズム終了 old_clusters <- cur_clusters cur_center <- .self$updateCenter(dat, cur_clusters) } return(as.integer(cur_clusters)) } ) )
k-means法を実行するときは例えば以下のようにする。
> source("kmeans.R") > km <- kmeans$new() > km$kmeans(x[3:5])
library(matrixStats) library(mvtnorm) source("kmeans.R") xmeans <- setRefClass( Class="xmeans", contains=c("kmeans"), fields=c( "num_split_cluster" ="integer", "stack" ="list", "Cluster" ="integer" ), methods=list( initialize=function() { initFields(num_split_cluster=0L, stack=list()) }, xmeans=function(DATA_) { ##### メイン処理 ##### data_ <- na.omit(DATA_[, sapply(DATA_, is.numeric)]) # 欠損レコード除去 rownames(data_) <- 1:nrow(data_) .self$splitData(data_) # まずデータを分割し切る .self$mergeSplitedData() # もしマージすべきだったらマージする .self$assignCluster() # 各レコードのクラスタ番号を取得する }, xBIC=function(dat_) { ##### ベイズ情報量基準 (分割前) ##### dat <- dat_[, sapply (dat_, is.numeric)] covar <- dat %>% summarise_each(funs(var)) %>% diag() df <- ncol(dat) * 2 center <- colMeans(dat) bic <- -2 * sum(dmvnorm(dat, center, covar, log=TRUE)) + df * log(nrow(dat)) return (bic) }, xBICp=function(dat_, cluster_) { ##### ベイズ情報量基準 (分割後) ##### dat <- dat_[, sapply(dat_, is.numeric)] cluster <- unlist(cluster_) covar_eachCluster <- cbind(dat, cc=cluster) %>% group_by(cc) %>% summarise_each(funs(var)) %>% dplyr::select(., -cc) %>% rowwise() %>% do(rtn=diag(.)) %>% as.list() %>% "$"(rtn) sumDeterminant <- sum(sapply(covar_eachCluster, function(mat) {prod(diag(mat))})) df <- 4 * ncol(dat) center_eachCluster <- cbind(dat, cc=cluster) %>% group_by(cc) %>% do(rtn=colMeans(dplyr::select(., -cc))) %>% as.list() %>% "$"(rtn) beta <- sqrt(c(crossprod(center_eachCluster[[1]] - center_eachCluster[[2]])) / sumDeterminant) logAlpha <- - log(2) - pnorm(beta, log.p=TRUE) logLikelihood_eachCluster <- sapply(seq_along(unique(cluster)), function (k) { center_ <- center_eachCluster[[k]] covar_ <- covar_eachCluster[[k]] dat_ <- dat[cluster==k, ] sum(dmvnorm(dat_, center_, covar_, log=TRUE))}) bicp <- - 2 * (nrow(dat) * logAlpha + sum(logLikelihood_eachCluster)) + df * log(nrow(dat)) return (bicp) }, split2Cluster=function(dat_) { ##### データを2クラスタに分割する ##### km <- kmeans$new() cluster <- km$kmeans(dat_, 2) # k-means で2クラスタにクラスタリング ccounts <- table(cluster) # 各クラスタに所属するデータ数を集計 if (any(ccounts < 3)) { # データ数3未満のクラスタができるようなら分割しない return (list(doSPLIT=FALSE)) } bicp <- .self$xBICp(dat_, cluster=cluster) bic <- .self$xBIC(dat_) if (bic > bicp) { # 分割後のベイズ情報量が小さくなるなら分割する return (list(doSPLIT=TRUE, d1=dat_[cluster==1,], d2=dat_[cluster==2,])) } else { return (list(doSPLIT=FALSE)) } }, merge2Cluster=function(dat_i, dat_j) { ##### 2クラスタをマージする ##### par <- data.frame(len=c(nrow(dat_i), nrow(dat_j)), val=1:2) cluster <- par %>% rowwise() %>% do(rval = rep(.$val, each=.$len)) %>% "$"(rval) %>% unlist() dat_ij <- rbind(dat_i, dat_j) bicp <- .self$xBICp(dat_ij, cluster) bic <- .self$xBIC(dat_ij) if (bic < bicp) { # マージ後のベイズ情報量が小さくなるならマージする return(list(doMERGE=TRUE, D=dat_ij)) } else { return(list(doMERGE=FALSE)) } }, splitData=function(data) { ##### データをクラスタに再帰的に分割する ##### # simple x-means split_data <- split2Cluster(data) if (split_data$doSPLIT) { splitData(split_data$d1) splitData(split_data$d2) } else { # データがもうそれ以上分割されなかった場合は、そのデータを最終的なクラスタとしてメンバに書き込んで終わる num_split_cluster <<- num_split_cluster + 1L stack[[num_split_cluster]] <<- data } }, mergeSplitedData=function() { ##### 分割し切ったクラスタにマージすべきものがあったらマージする ##### ind_order <- .self$getDataSize_order(.self$getDataSize_eachCluster()) num_cluster <- length(ind_order) isMerge_eachCluster <- numeric(num_cluster) for (i in seq_len (num_cluster - 1)) { i_ <- ind_order[sprintf ("%d", i)] if (isMerge_eachCluster[i_] == 1) next for (j in (i+1):num_cluster) { j_ <- ind_order[sprintf ("%d", j)] if (isMerge_eachCluster[j_] == 0) { merge_data <- merge2Cluster(stack[[i_]], stack[[j_]]) if (merge_data$doMERGE) { stack[[i_]] <<- merge_data$D isMerge_eachCluster[j_] <- 1 break } } } } stack <<- stack[isMerge_eachCluster==0] num_split_cluster <<- length(stack) }, getDataSize_eachCluster=function() { out <- sapply(stack, nrow) return (out) }, getDataSize_order=function(DataSize_eachCluster) { out <- seq_len(length(DataSize_eachCluster)) out <- setNames(out, order(DataSize_eachCluster)) return(out) }, assignCluster=function() { out <- numeric(sum(.self$getDataSize_eachCluster())) ind <- lapply(stack, function(d){as.integer(rownames(d))}) for (i in seq_along(ind)) { out[ind[[i]]] <- i } Cluster <<- as.integer(out) } ) )