満島さん薄幸似合いすぎ. フォルダ5聴き直すわ.
あっ, どーも僕です.
下痢と肩こりが最近の悩みです.
flexmixとかいうパッケージすげぇ
昨日モデルベースのEMアルゴリズムを使った論文を読みました.
EMアルゴリズムは一年くらい前に勉強して, 結局実装できなくて苦手意識バリバリだったのでうげええええって感じでした. しかも, 勉強したのはモデルベースじゃなかったのでもっと難しくなっているがな(・∀・)って感じです.
再挑戦しようとしてやっぱりなかなか苦戦していたところ, なんとRのflexmix*1パッケージでモデルベースのEMアルゴリズムを提供しているというのでわないですが！！
Rのflexmixパッケージで混合分布モデルによるクラスタ分析を行う。 - Analyze IT.
うほほーいって感じです. ってことでさっそく論文に載っていた例題の正規分布モデルを試してみました. まず, データはこんな感じに三種類の直線モデルに正規分布誤差を与えています. つまり, 潜在的クラスの正しい数は3ですね.
色を抜くとこんな感じです. こんな感じのデータがあったときに, 上のようなモデルの帰属するように分類するにEMアルゴリズムが使えるのです. しかも, 一緒にモデルが求まるというのです.
この図をみたときに, 直感的にはデータが２つの発生源からのものが混在しているように見えますが, ３っつはどうかなって感じですね.
で, 今回はお試しなので, 最初に与えるデータのグループ数は3にしてflexmixパッケージを使ってみました. あっさりできました. 結果が下の図です. 色はEMアルゴリズムがによりわけられたグループの色を使っています. すげええ. まだ, 使い方よくわからないのでパラメータはデフォルトのままなのになかなかよく分類していますね. 線が交差しているところはパラメータをいじればもっとうまくいくのかな？図中に引いてある直線はEMアルゴリズムが推定したモデルです.
ここまでのコードはこんな感じです.
# 架空データの用意 x <- seq (0, 5, length = 50) a1 <- 2; b1 <- 10; sd1 <- 2 a2 <- 5; b2 <- 15; sd2 <- 3 a3 <- 9; b3 <- 5; sd3 <- 2 y1 <- a1 * x + b1 + rnorm (50, m = 0, sd = sd1) y2 <- a2 * x + b2 + rnorm (50, m = 0, sd = sd2) y3 <- a3 * x + b3 + rnorm (50, m = 0, sd = sd3) # 一番目の図 color <- brewer.pal(3, "Set2") def.par <- par (no.readonly = TRUE) ps.options(family = "Japan1") par (mai = c(1.5, 1.5, .5, .5)) plot (0, type = "n", xlim = c (0, 5), ylim = c (0, 60), cex.main = 1.8, xlab = "X", ylab = "Y", tcl = .3, cex.axis = 1.8, cex.lab = 1.8, las = 1) abline (b1, a1, lwd = 1) abline (b2, a2, lwd = 1) abline (b3, a3, lwd = 1) points (x, y1, pch = 19, cex = 2, col = color[3]) points (x, y2, pch = 19, cex = 2, col = color[2]) points (x, y3, pch = 19, cex = 2, col = color[1]) par (def.par) # 二番目の図 def.par <- par (no.readonly = TRUE) ps.options(family = "Japan1") par (mai = c(1.5, 1.5, .5, .5)) plot(X, Y, col = "darkgrey", pch = 19, cex = 2,xlim = c (0, 5), ylim = c (0, 60), cex.main = 1.8, xlab = "X", ylab = "Y", tcl = .3, cex.axis = 1.8, cex.lab = 1.8, las = 1) par(def.par) # ***************** flexmix *************************** library (flexmix) Y <- c(y1, y2, y3) X <- rep (x, times = 3) fit <-FLXMRglm(Y ~ X, family="gaussian") # 線形回帰の正規分布誤差 dat <- data.frame (X = X, Y = Y) # EMアルゴリズムによる推定 rtn <- flexmix(Y ~ X, data=dat, k=3, model=list(fit), control=list(verbose=1, nrep=5)) rtn@cluster ; table (rtn@cluster) rtn@df rtn@logLik rtn@size #1 2 3 #34 66 50 rtn@model rtn@components rtn_coef <- ldply (1:3, function (i) unlist(attr((unlist(rtn@components))[[i]], "parameters"))) # 三番目の図 def.par <- par (no.readonly = TRUE) ps.options(family = "Japan1") par (mai = c(1.5, 1.5, .5, .5)) plot(X, Y, col = color [rtn@cluster], pch = 19, cex = 2,xlim = c (0, 5), ylim = c (0, 60), cex.main = 1.8, xlab = "X", ylab = "Y", tcl = .3, cex.axis = 1.8, cex.lab = 1.8, las = 1) for (i in 1:3) { abline (rtn_coef[i, 1], rtn_coef[i,2]) } par(def.par)
いやーすごいですね. 普段はなるだけ自分でつくるように心がけているのですが, こんなすごいテクニックがこんなにあっさりできてしまうと, パッケージ中毒になりそうです.
おまけ：モデル選択
おまけとして, 今回は最初からグループ数3を決め打ちでしたが実際には自分で試行錯誤しなければなりません. ってことで, 上のような問題設定のもとデータセットを10000用意してグループ数2, 3, 4, 5に対してシュミレーションしAIC, BICを計算してみました.
どちらかと言えばBICの方がいいのかな？
コードはこちら. 結果をまとめるのに少し面倒くさい形で書いているのは, EMアルゴリズムの最適化過程でグールプ数が落ちてしまいデータフレームにまとめられなかったからです. これはパラメータをいじれば解消できるのかな？
# データの用意 X <- rep (seq (0, 5, length = 50), times = 3) a1 <- 2; b1 <- 10; sd1 <- 2 a2 <- 5; b2 <- 15; sd2 <- 3 a3 <- 9; b3 <- 5; sd3 <- 2 Nsim <- 10000 Y <- sapply (seq_len (Nsim) , function (i) return (c(a1 * x + b1 + rnorm (50, m = 0, sd = sd1), # シュミレーション用関数 simMyflex <- function (m, Y) { # m 分割数 # Y 被説明変数 nc <- ncol (Y) tmp <- llply (seq_len (nc), function (i) { y <- Y[ , i, drop = TRUE] dat <- data.frame (X = X, Y = y) fit <- FLXMRglm(Y ~ X, family="gaussian") result <- flexmix (Y ~ X, data = dat, k = m, model = list (fit), control = list (verbose = 1, nrep = 10)) return (result) }) # グループ数が変わってしまうときがあるのでA2はllply A1 <- ldply (tmp, function (X) MyIC (X)) A2 <- llply (tmp, function (X) MyCoef (X)) return (list (A1, A2)) } result.2 <- simMyflex (m = 2, Y) result.3 <- simMyflex (m = 3, Y) result.4 <- simMyflex (m = 4, Y) result.5 <- simMyflex (m = 5, Y) # ******************************************** # サブルーチン的ななにか1 MyIC <- function (flexmix) { n <- sum (flexmix@size) K <- length (flexmix@size) df <- flexmix@df IC <- c (negativeLL = - flexmix@logLik, AIC = - 2 * flexmix@logLik + 2 * df, BIC = - 2 * flexmix@logLik + log (n) * df, K = K ) return (IC) } # サブルーチン的ななにか2 MyCoef <- function (flexmix) { coef <- llply (seq_along (names(flexmix@components)), function (i) unlist(attr((unlist(flexmix@components))[[i]], "parameters"))) return (c(unlist(coef), flexmix@size)) } # **************************************************** # 結果のグラフ化 D <- do.call (rbind, list (result.2[[1]], result.3[[1]], result.4[[1]], result.5[[1]])) # AIC def.par <- par(no.readonly = TRUE) ps.options(family = "Japan1") par(mai = c(1.5, 1.5, 1.5, .5)) boxplot (AIC ~ K, data = D, tcl = -.5, cex.axis = 1.8, cex.lab = 1.8, cex.main = 1.8, main = "AIC", ylab = "AIC", xlab = "グループ数") par(def.par) # BIC def.par <- par(no.readonly = TRUE) ps.options(family = "Japan1") par(mai = c(1.5, 1.5, 1.5, .5)) boxplot (BIC ~ K, data = D, tcl = -.5, cex.axis = 1.8, cex.lab = 1.8, cex.main = 1.8, main = "BIC", ylab = "BIC", xlab = "グループ数") par(def.par)
今回は簡単な正規誤差の線形モデルでしたがブックマーク先やインストールしたときについてくるPDFをみるともっといろいろできそうですね. すごすぎ.
*1:Friedrich Leisch. FlexMix: A general framework for finite mixture models and latent class regression in R. Journal of Statistical Software, 11(8), 1-18, 2004. URL http://www.jstatsoft.org/v11/i08/ Bettina Gruen and Friedrich Leisch. Fitting finite mixtures of generalized linear regressions in R. Computational Statistics & Data Analysis, 51(11), 5247-5252, 2007. doi:10.1016/j.csda.2006.08.014 Bettina Gruen and Friedrich Leisch. FlexMix Version 2: Finite mixtures with concomitant variables and varying and constant parameters Journal of Statistical Software, 28(4), 1-35, 2008. URL http://www.jstatsoft.org/v28/i04/
