MCMCって大量にサンプルを発生させる技ですが，あー計算しているなー，発生させているなー，というのが目に見えてわかるようになったら楽しいと思うんです。イメージもしやすいだろうし。
ということで，Rでアニメーション画像を作る方法はないかしら，と考えたところ，animationパッケージとimage magickというのを使えば良いらしい，ということがわかりました。この2つのパッケージとアプリが用意されている前提で，次の関数を作りました。
MCMCmovie <- function(dat,n.step=10,dens=FALSE){
L <- length(dat)
B <- nclass.Sturges(dat)
for(i in seq(100,L,n.step)){
hist(dat[1:i],xlim=range(dat),main=paste("MCMC iter at ",i),breaks=B,freq=FALSE)
if(dens){lines(density(dat[1:i]), col = "orange", lwd = 2)}
}
# finish
hist(dat,xlim=range(dat),main=paste("MCMC iter at ",L),breaks=B,freq=FALSE)
if(dens){lines(density(dat[1:i]), col = "orange", lwd = 2)}
}
与えるデータは，rstanが吐き出したstanfitオブジェクトのなかでも，必要な変数のところだけextractしたやつです。
オプションとして，stepが描画するときのサンプルをいくつ飛ばしでやるか(間引き方ですが，thinだとstanの変数とごっちゃになるので名前を変えました)，密度関数の線を引くかどうか(dens=TRUEで書きます。デフォルトはFALSE)，スタートはバーンアウト後100回目のサンプルからにしてます。
試しに，世界一簡単なrstanコードを使って実際に使用してみたいと思います。
require(rstan)
n <- 1000
mu <- 50
sig <- 10
y <- rnorm(n,mu,sig)
stancode <- "
data{
int<lower=0> T;
real N[T];
}
parameters{
real mu;
real<lower=0> s;
}
model{
N ~ normal(mu,sqrt(s));
s~cauchy(0,5);
}
"
datastan <- list(N=y,T=n)
fit <- stan(model_code = stancode,data=datastan,iter=5000,chains=4)
## ## TRANSLATING MODEL 'stancode' FROM Stan CODE TO C++ CODE NOW. ## COMPILING THE C++ CODE FOR MODEL 'stancode' NOW. ## ## SAMPLING FOR MODEL 'stancode' NOW (CHAIN 1). ## ## Iteration: 1 / 5000 [ 0%] (Warmup) ## Iteration: 500 / 5000 [ 10%] (Warmup) ## Iteration: 1000 / 5000 [ 20%] (Warmup) ## Iteration: 1500 / 5000 [ 30%] (Warmup) ## Iteration: 2000 / 5000 [ 40%] (Warmup) ## Iteration: 2500 / 5000 [ 50%] (Warmup) ## Iteration: 2501 / 5000 [ 50%] (Sampling) ## Iteration: 3000 / 5000 [ 60%] (Sampling) ## Iteration: 3500 / 5000 [ 70%] (Sampling) ## Iteration: 4000 / 5000 [ 80%] (Sampling) ## Iteration: 4500 / 5000 [ 90%] (Sampling) ## Iteration: 5000 / 5000 [100%] (Sampling) ## # Elapsed Time: 0.902595 seconds (Warm-up) ## # 0.743175 seconds (Sampling) ## # 1.64577 seconds (Total) ## ## ## SAMPLING FOR MODEL 'stancode' NOW (CHAIN 2). ## ## Iteration: 1 / 5000 [ 0%] (Warmup) ## Iteration: 500 / 5000 [ 10%] (Warmup) ## Iteration: 1000 / 5000 [ 20%] (Warmup) ## Iteration: 1500 / 5000 [ 30%] (Warmup) ## Iteration: 2000 / 5000 [ 40%] (Warmup) ## Iteration: 2500 / 5000 [ 50%] (Warmup) ## Iteration: 2501 / 5000 [ 50%] (Sampling) ## Iteration: 3000 / 5000 [ 60%] (Sampling) ## Iteration: 3500 / 5000 [ 70%] (Sampling) ## Iteration: 4000 / 5000 [ 80%] (Sampling) ## Iteration: 4500 / 5000 [ 90%] (Sampling) ## Iteration: 5000 / 5000 [100%] (Sampling) ## # Elapsed Time: 0.76551 seconds (Warm-up) ## # 0.978063 seconds (Sampling) ## # 1.74357 seconds (Total) ## ## ## SAMPLING FOR MODEL 'stancode' NOW (CHAIN 3). ## ## Iteration: 1 / 5000 [ 0%] (Warmup) ## Iteration: 500 / 5000 [ 10%] (Warmup) ## Iteration: 1000 / 5000 [ 20%] (Warmup) ## Iteration: 1500 / 5000 [ 30%] (Warmup) ## Iteration: 2000 / 5000 [ 40%] (Warmup) ## Iteration: 2500 / 5000 [ 50%] (Warmup) ## Iteration: 2501 / 5000 [ 50%] (Sampling) ## Iteration: 3000 / 5000 [ 60%] (Sampling) ## Iteration: 3500 / 5000 [ 70%] (Sampling) ## Iteration: 4000 / 5000 [ 80%] (Sampling) ## Iteration: 4500 / 5000 [ 90%] (Sampling) ## Iteration: 5000 / 5000 [100%] (Sampling) ## # Elapsed Time: 0.831019 seconds (Warm-up) ## # 0.860767 seconds (Sampling) ## # 1.69179 seconds (Total) ## ## ## SAMPLING FOR MODEL 'stancode' NOW (CHAIN 4). ## ## Iteration: 1 / 5000 [ 0%] (Warmup) ## Iteration: 500 / 5000 [ 10%] (Warmup) ## Iteration: 1000 / 5000 [ 20%] (Warmup) ## Iteration: 1500 / 5000 [ 30%] (Warmup) ## Iteration: 2000 / 5000 [ 40%] (Warmup) ## Iteration: 2500 / 5000 [ 50%] (Warmup) ## Iteration: 2501 / 5000 [ 50%] (Sampling) ## Iteration: 3000 / 5000 [ 60%] (Sampling) ## Iteration: 3500 / 5000 [ 70%] (Sampling) ## Iteration: 4000 / 5000 [ 80%] (Sampling) ## Iteration: 4500 / 5000 [ 90%] (Sampling) ## Iteration: 5000 / 5000 [100%] (Sampling) ## # Elapsed Time: 0.803299 seconds (Warm-up) ## # 0.812294 seconds (Sampling) ## # 1.61559 seconds (Total)
print(fit,digit=3)
## Inference for Stan model: stancode. ## 4 chains, each with iter=5000; warmup=2500; thin=1; ## post-warmup draws per chain=2500, total post-warmup draws=10000. ## ## mean se_mean sd 2.5% 25% 50% 75% ## mu 49.778 0.004 0.312 49.170 49.565 49.780 49.984 ## s 100.270 0.059 4.338 92.089 97.305 100.092 103.174 ## lp__ -2806.161 0.016 0.941 -2808.686 -2806.539 -2805.866 -2805.482 ## 97.5% n_eff Rhat ## mu 50.400 6106 1.000 ## s 109.131 5356 1.001 ## lp__ -2805.234 3644 1.002 ## ## Samples were drawn using NUTS(diag_e) at Wed Jun 17 18:23:28 2015. ## For each parameter, n_eff is a crude measure of effective sample size, ## and Rhat is the potential scale reduction factor on split chains (at ## convergence, Rhat=1).
histdat <- extract(fit,pars="mu")$mu
最後の行でextractsしてます。
さて，このサンプルをアニメーションgifにします。
animationパッケージをrequireして，オプションをいくつか。intervalはパラパラアニメのスピードですが，0.06ぐらいが丁度いいようです。loopオプションは反復するかどうかで，デフォルトではgifが無限ループするので1にしておきます。この辺はanimationパッケージのsaveGIFオプションなので，ヘルプ見てください。
ここでsaveGIF関数の最初に，上で用意した自作関数を入れてやると，どんどんヒストグラムが育っていきます。
require(animation)
ani.options(interval=0.06)
ani.options(loop=1)
saveGIF(MCMCmovie(histdat,dens=T,n.step=20), movie.name="MCMCexample1.gif",width=640, height=480,clean=T)
## [1] FALSE
これを実行して得られる絵がこちら。(動かない人は画像をクリックしてね。)
モコモコしてる！MCMCがMCMCしてる！
・・・しかしここまできたら，ggplot2パッケージをつかってもう少し綺麗に描きたいな。
ということで，書き直したのが次のMCMCmovie2関数。
MCMCmovie2 <- function(dat,nchain=1,n.step=10,dens=F){
require(ggplot2)
dat <- as.data.frame(matrix(dat,ncol=nchain))
L <- nrow(dat)
for(i in seq(100,L,n.step)){
dat.tmp <- as.data.frame(matrix(unlist(dat[1:i,]),ncol=1))
cname <- c()
for(j in 1:nchain){cname <- c(cname,rep(paste("Chain",j),i))}
dat.tmp$cname <- as.factor(cname)
#ggplotting
g <- ggplot(dat.tmp,aes(x=V1,y=..density..,colour=cname,fill=cname))
g <- g + ggtitle(paste("MCMC iter at ",i)) + labs(x="",y="density")
if(dens==TRUE){
g <- g + geom_density(alpha=0,position="identity")
}else{
g <- g + geom_histogram(alpha=0.3,position="identity")
}
print(g)
}
}
使い方として，引数にchain数を入れるようにしてます。そうするとchainごとに色分けしてくれる。chainの区別なしでいいや，という人は何も書かなければ全部を1サンプルとして描いてしまいます。
histdat <- extract(fit,pars="mu")$mu
require(animation)
ani.options(interval=0.06)
ani.options(loop=0)
saveGIF(MCMCmovie2(histdat,nchain=4,dens=TRUE,n.step=10), movie.name="MCMCexample2.gif",width=640, height=480,clean=T)
## Error in .Call("loop_apply", as.integer(n), f, env): "loop_apply" not resolved from current namespace (plyr)
これを実行して得られる絵がこちら。
モコモコしてる！MCMCがMCMCしてる！(；´Д｀)ﾊｱﾊｱ
ちなみに棒グラフにする(dens=FALSE設定)と（なんか警告が出てうるさいですけど）
モコモコしてる！MCMCがMCMCしてる！(；´Д｀)ﾊｱﾊｱ(；´Д｀)ﾊｱﾊｱ
・・・誰かの何かのお役に立てば。