雑記: X-means法

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 にまとめてから読み込みます。

北海道,北海道,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)

f:id:cookie-box:20170301131830p:plain:w440
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)

f:id:cookie-box:20170301164808p:plain:w690

それで、このデータを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)

f:id:cookie-box:20170302000449p:plain:w640

上の結果より、最低賃金合計特殊出生率、高齢化率をスケーリングして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"))

f:id:cookie-box:20170302124428p:plain:w500

上の散布図とクラスタが微妙に違います。元々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) { ##### クラスタの初期化 (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])

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_) { ##### データを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)
    }
  )
)