最近、日本の出版社からR本出過ぎ
海外の出版社もですが。
あ、どーも僕です。
上級ハンドブックの目次みれるようになりましたね。
角度データなる初耳の単語があったり、苦手なネットワーク関連があったりしてとても楽しみです。
これ出る前にワンちゃん本読み終わりたいです。
- 作者: 荒引健,石田基広,高橋康介,林真広,二階堂愛
- 出版社/メーカー: シーアンドアール研究所
- 発売日: 2013/09/25
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (10件) を見る
それは置いといて、学会が終わり論文投稿も間に合いそうだし
就活もなんとなかなりそうになったので(ま、どっちもまだ終わっていないけど)
久しぶりにブログ書きます。
いまさらだけど、RでEMアルゴリズム組めたよって話
なんですが、EMアルゴリズム [R]とか検索すると山ほど
同じ話題を扱ったブログがでてきますし、
正規分布のEMアルゴリズムを扱ったパッケージ*1はたくさんあります。
なんで、ここでは少しだけでも新規性があるように収束の様子をアニメーションしてみました。自分が調べた限りでは収束過程を出力していた方はいらしたのですが、アニメーションにはしている人はいなかったです。
対数尤度が増加する様子のグラフは載せていませんが、下の方にスクリプトはあるのでそれ実行してください。
追加パッケージはmvtnorm*2, RColorBrewer*3, animation*4です。
以下に二変量のときを示します。相関0の二変量正規分布を2つ混ぜています。初期値から徐々にそれらしい正規分布になっているのがわかりますね。
スクリプト一部修正しました(日付は忘れた。。)
スクリプトを書き直しました(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)
書き直しましたものでは次のようなグラフで動作確認できます。左がクラスタリングの結果で, 右が対数尤度が単調増加していることの確認です。
いいわけ
いまさら感がある話題ですが、実際わたしも取り組んだのは結構前でした。
そのときはなぜかいい値に収束しているのに、対数尤度が増えなくて困っていたのですが
どうやら対数尤度を間違って計算していました。
最近、なんとなく未解決プログラムを見なおしていたらそのことに気づいたいので
嬉しくてつい載せてしまいました。
このプログラムを書くための参考文献は小西先生の本です。
プログラム例が載っているわけでもないのですが、シンプルでわかりやすいです。
- 作者: 小西貞則
- 出版社/メーカー: 岩波書店
- 発売日: 2010/01/27
- メディア: 単行本(ソフトカバー)
- 購入: 14人 クリック: 347回
- この商品を含むブログ (8件) を見る
*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/.