最近、日本の出版社からR本出過ぎ

海外の出版社もですが。


あ、どーも僕です。



上級ハンドブックの目次みれるようになりましたね。
角度データなる初耳の単語があったり、苦手なネットワーク関連があったりしてとても楽しみです。
これ出る前にワンちゃん本読み終わりたいです。

R言語上級ハンドブック

R言語上級ハンドブック




それは置いといて、学会が終わり論文投稿も間に合いそうだし
就活もなんとなかなりそうになったので(ま、どっちもまだ終わっていないけど)
久しぶりにブログ書きます。





いまさらだけど、RでEMアルゴリズム組めたよって話

多変量正規分布を仮定したEMアルゴリズムをRで組みました。

なんですが、EMアルゴリズム [R]とか検索すると山ほど
同じ話題を扱ったブログがでてきますし、
正規分布EMアルゴリズムを扱ったパッケージ*1はたくさんあります。

なんで、ここでは少しだけでも新規性があるように収束の様子をアニメーションしてみました。自分が調べた限りでは収束過程を出力していた方はいらしたのですが、アニメーションにはしている人はいなかったです。
対数尤度が増加する様子のグラフは載せていませんが、下の方にスクリプトはあるのでそれ実行してください。

追加パッケージはmvtnorm*2, RColorBrewer*3, animation*4です。

以下に二変量のときを示します。相関0の二変量正規分布を2つ混ぜています。初期値から徐々にそれらしい正規分布になっているのがわかりますね。


f:id:aaaazzzz036:20130917175123g:plain

スクリプト一部修正しました(日付は忘れた。。)
スクリプトを書き直しました(2013/12/1)

書き直しましたものは, ただクラスタリングをするだけです。動作確認はしたのほうにグラフがあります。計算過程を保存していないので上のグラフが描きたいときは結果を保存するようにスクリプト直してください。

# **********************************************************************************
# 多変量正規EMアルゴリズム
# args 
# - data : matrix (data, variates) , データフレームでもよい. 変数は二変量以上
# - K : 分割するクラスター数
# - epsilon : 収束条件, 尤度の変化量
# - iter.max : 最大収束計算数
# return
# - llh : 対数尤度の記録
# - cluster : データのクラスター
# - iter : 繰り返し計算数
# - mixing.rate : 混合率
# - centers : 各クラスターの平均値ベクトル(variates, cluster)
# - covars : 各クラスターの共分散行列
# **********************************************************************************
myEM.multinormal <- function (data, 
                              K = 2, 
                              epsilon = 1e-5, 
                              iter.max = 1000) {
    # 前処理
    data <- na.omit (as.matrix (data))
    nr   <- nrow (data)
    nc   <- ncol (data)
    K    <- as.integer (K)
    if (nr < nc) return (0)
    if (nc < 2) return (1)
    if (any (!is.numeric (data))) return (2)
    if (K < 2) return (3)
    
    # 追加パッケージ
    require (matrixStats)
    require (mvtnorm)
    
    # 初期クラスターでの尤度計算
    doInitialize <- TRUE
    while (doInitialize) {
        cluster.old  <- sample (1:K, nr, replace = TRUE)
        if (length (unique (cluster.old)) == K & all (table (cluster.old) >= 2)) {
            doInitialize <- FALSE
        }
    }

    cluster.centers <- by (data, cluster.old, colMeans)
    cluster.covars  <- by (data, cluster.old, cov)
    mixing.rate     <- table (cluster.old) / nr
    cluster.density <- sapply (1:K, function (k) {
        center <- cluster.centers[[k]]
        covar  <- cluster.covars[[k]]
        dmvnorm (data, center, covar) * mixing.rate[k]
    })
    mixed.density     <- rowSums (cluster.density)
    cluster.prob      <- sweep (cluster.density, 1, mixed.density, "/")
    logLikelihood.old <- sum (log (mixed.density))
        
    # 尤度の単調増加を確認するために対数尤度だけは記録していく
    logLikelihood.history <- numeric (iter.max)
    logLikelihood.history[1] <- logLikelihood.old

    # 収束計算
    iter   <- 1
    doNEXT <- TRUE
    while (doNEXT) {
        iter <- iter + 1
        if (iter > iter.max) break
        
        # step1 各クラスのデータ数
        cluster.n <- colSums (cluster.prob)
        
        # step2 Mstep
        cluster.centers <- sweep (crossprod (data, cluster.prob),  2, cluster.n, "/")
        cluster.covars  <- lapply (1:K, function (k) {
            data.shifted   <- sweep(data, 2, cluster.centers[, k, drop = TRUE], "-") 
            data.shifted.p <- sweep(data.shifted,  1, cluster.prob[ , k], "*")           
            crossprod(data.shifted.p, data.shifted) / cluster.n[k] 
        })
        
        # step3 事後確率
        mixing.rate <- cluster.n / nr
        cluster.density <- sapply (1:K, function (k) {
            center <- cluster.centers[, k, drop = TRUE]
            covar  <- cluster.covars[[k]]
            dmvnorm (data, center, covar) * mixing.rate[k]
        })

        # step4 Estep
        mixed.density     <- rowSums (cluster.density)
        logLikelihood.new <- sum (log(mixed.density))
        cluster.prob      <- sweep (cluster.density, 1, mixed.density, "/")
        
        # step5 収束判定 
        doNEXT <- abs (logLikelihood.new - logLikelihood.old) > epsilon  
        logLikelihood.old <- logLikelihood.new
        logLikelihood.history[iter] <- logLikelihood.old[1]
        
    }
    logLikelihood.history <- logLikelihood.history [logLikelihood.history != 0]
    
    return (list (llh = logLikelihood.history, 
                  iter = iter, 
                  mixing.rate = mixing.rate, 
                  cluster = max.col (cluster.prob), 
                  centers = cluster.centers, # 列が各クラスの平均値ベクトル
                  covars = cluster.covars))
  
}

# 動作確認
require (mvtnorm)
col <- c("skyblue", "yellowgreen")
x <- rbind(rmvnorm(100, mean = c(0, 0), sigma = diag(2)), rmvnorm (100, mean = c(3, 3), sigma = diag(2)))
print (t1 <- myEM.multinormal (x, epsilon = 1e-18, K = 2))
layout (matrix (c(1, 2), 1))
plot (x, pch = t1$cluster, cex = 2, asp = 1, col = col[t1$cluster], 
      xlab = "x", ylab = "y")
points (t(t1$centers), pch = "X", cex = 2)
grid(lty = 5, col = "darkgrey")
plot (1:t1$iter, t1$llh, type = "b", col = "orange", 
      xlab = "iterations", ylab = "loglikelihood", cex.main =2)

書き直しましたものでは次のようなグラフで動作確認できます。左がクラスタリングの結果で, 右が対数尤度が単調増加していることの確認です。

f:id:aaaazzzz036:20131201175646p:plain

いいわけ

いまさら感がある話題ですが、実際わたしも取り組んだのは結構前でした。
そのときはなぜかいい値に収束しているのに、対数尤度が増えなくて困っていたのですが
どうやら対数尤度を間違って計算していました。
最近、なんとなく未解決プログラムを見なおしていたらそのことに気づいたいので
嬉しくてつい載せてしまいました。

このプログラムを書くための参考文献は小西先生の本です。
プログラム例が載っているわけでもないのですが、シンプルでわかりやすいです。

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

*1:例えば、Tatiana Benaglia, Didier Chauveau, David R. Hunter, Derek Young (2009). mixtools: An R Package for Analyzing Finite Mixture Models. Journal of Statistical Software, 32(6), 1-29. URL http://www.jstatsoft.org/v32/i06/.

*2: Alan Genz, Frank Bretz, Tetsuhisa Miwa, Xuefei Mi, Friedrich Leisch, Fabian Scheipl, Torsten Hothorn (2013). mvtnorm: Multivariate Normal and t Distributions. R package version 0.9-9995. URL http://CRAN.R-project.org/package=mvtnorm

*3: Erich Neuwirth (2011). RColorBrewer: ColorBrewer palettes. R package version 1.0-5. http://CRAN.R-project.org/package=RColorBrewer

*4: Yihui Xie (2013). animation: An R Package for Creating Animations and Demonstrating Statistical Methods. Journal of Statistical Software, 53(1), 1-27. URL http://www.jstatsoft.org/v53/i01/.