MCMCでモコモコしちゃう

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

これを実行して得られる絵がこちら。(動かない人は画像をクリックしてね。)

MCMCexample1

モコモコしてる!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)

これを実行して得られる絵がこちら。

MCMCexample2

モコモコしてる!MCMCがMCMCしてる!(;´Д`)ハアハア

ちなみに棒グラフにする(dens=FALSE設定)と(なんか警告が出てうるさいですけど)

MCMCexample3

モコモコしてる!MCMCがMCMCしてる!(;´Д`)ハアハア(;´Д`)ハアハア

・・・誰かの何かのお役に立てば。