X-means アルゴリズムを実装してみます。
手順
以下のように最低賃金、合計特殊出生率、高齢化率を data.txt にまとめてから読み込みます。
北海道,北海道,764,1.29,0.281
青森,東北,695,1.43,0.290
岩手,東北,695,1.50,0.
296
宮城,東北,726,1.31,0.246
（以下略）
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)
2次元に投影したプロットもしておきます。
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
クラスタに分割されることがわかります。グラフからわかる各
クラスタの特徴は以下です。
追記
上の結果を日本地図にプロットします。以下の記事を全面的に参考にさせていただきました。
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++法で固定、細分したクラスタをマージするかどうかはするで固定、共分散行列の非対角成分（共分散）は無視で固定しています。元記事のスクリプトは引数で変更できるようになっています。
kmeans.R
library(dplyr)
kmeans <- setRefClass(
Class="kmeans",
methods=list(
initCluster_kmeanspp=function(dat, NUM_CLUSTER) {
index <- 1:nrow(dat)
center_index <- sample(index, 1)
center_ <- unlist(dat[center_index,])
min_dist <- rowSums((sweep(dat, 2, center_, "-"))**2)
while (length(center_index) < NUM_CLUSTER) {
prob <- min_dist[-center_index] / sum(min_dist[-center_index])
new_center <- sample(index[-center_index], 1, prob=prob)
center_ <- unlist(dat[new_center,])
new_dist <- rowSums((sweep(dat, 2, center_, "-"))**2)
min_dist <- pmin(min_dist, new_dist)
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])
xmeans.R
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_) {
km <- kmeans$new()
cluster <- km$kmeans(dat_, 2)
ccounts <- table(cluster)
if (any(ccounts < 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) {
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) {
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)
}
)
)
