RStan使ってMCMCの計算をやってみた(1)
久保拓弥先生の「データ分析のための統計モデリング入門」(通称緑本)の9章をトレースしてみました。緑本ではWinBUGSを使用していますが、ここではRStanを使ってStanでやってみます。
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者:久保 拓弥
- 発売日: 2012/05/19
- メディア: 単行本
インストール
Stanのインストールに関しては、以下のブログ記事を参考にさせていただきました。
Stanで統計モデリングを学ぶ(1): まずはStanの使い方のおさらいから - 渋谷駅前で働くデータサイエンティストのブログ
こちらの通りに進めていけば問題なくインストールできるかと思います。
コード
コードの書き方や、結果の可視化に関してはインストールでもお世話になった先ほどの記事と、こちらの記事
http://heartruptcy.blog.fc2.com/blog-entry-89.html
が大変参考になりました。緑本ではBUGSコードが記載されていますが、ここではStanコードを用意します。以下のコードをkubo9.stanという名前で保存します。
data { int<lower=0> N; real X[N]; real MeanX; int<lower=0> y[N]; } parameters { real beta1; real beta2; } model { for (i in 1:N){ y[i] ~ poisson(exp(beta1 + beta2 * (X[i] - MeanX))); } beta1 ~ normal(0, 100); //無情報事前分布 beta2 ~ normal(0, 100); //無情報事前分布 }
$y[i] \sim {\rm poisson}(\lambda)\ $の部分で$\ y\ $が平均$\ \lambda=\exp(\beta_1 + \beta_2x)\ $のポアソン分布に従うことを表しています。次に、以下のコードでStanに渡すデータをlistとして作成して、StanにMCMCサンプリングを実行させます。データは著者のサイト
http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html
からダウンロードできます。
library(rstan) load("d.Rdata") d.dat <- list(N=dim(d)[1], X=d$x, y=d$y, MeanX=mean(d$x)) #Stanに渡すデータ d.fit <- stan(file='kubo9.stan', data=d.dat, iter=1600, warmup=100, thin=3, chains=3)
緑本のn.burnin = 100にあたるのがwarmup = 100です。これを実行すると
> d.fit <- stan(file='kubo9.stan', data=d.dat, iter=1600, warmup=100, thin=3, chains=3) TRANSLATING MODEL 'kubo9' FROM Stan CODE TO C++ CODE NOW. COMPILING THE C++ CODE FOR MODEL 'kubo9' NOW. cygwin warning: MS-DOS style path detected: C:/R/R-30~1.1/etc/x64/Makeconf Preferred POSIX equivalent is: /cygdrive/c/R/R-30~1.1/etc/x64/Makeconf CYGWIN environment variable option "nodosfilewarning" turns off this warning. Consult the user's guide for more details about POSIX paths: http://cygwin.com/cygwin-ug-net/using.html#using-pathnames C:/R/R-3.0.1/library/rstan/include//stansrc/stan/agrad/rev/var_stack.hpp:49:17: warning: 'void stan::agrad::free_memory()' defined but not used [-Wunused-function] C:/R/R-3.0.1/library/rstan/include//stansrc/stan/agrad/rev/chainable.hpp:87:17: warning: 'void stan::agrad::set_zero_all_adjoints()' defined but not used [-Wunused-function] SAMPLING FOR MODEL 'kubo9' NOW (CHAIN 1). Iteration: 1600 / 1600 [100%] (Sampling) Elapsed Time: 0.011 seconds (Warm-up) 0.057 seconds (Sampling) 0.068 seconds (Total) SAMPLING FOR MODEL 'kubo9' NOW (CHAIN 2). Iteration: 1600 / 1600 [100%] (Sampling) Elapsed Time: 0.004 seconds (Warm-up) 0.058 seconds (Sampling) 0.062 seconds (Total) SAMPLING FOR MODEL 'kubo9' NOW (CHAIN 3). Iteration: 1600 / 1600 [100%] (Sampling) Elapsed Time: 0.004 seconds (Warm-up) 0.053 seconds (Sampling) 0.057 seconds (Total)
となり、MCMCサンプリングを3回反復して計算終了となります。
d.fit
とすれば結果が見れます。簡易的に可視化したければ
plot(d.fit) traceplot(d.fit)
で以下の図が得られます。
結果の図示
以下の要領で$\ \beta_1\ $の結果を図示できます。
library(ggplot2) library(reshape2) library(grid) d.ext <- extract(d.fit, permuted=F) b1.1 <- d.ext[1:500,'chain:1','beta1'] b1.2 <- d.ext[1:500,'chain:2','beta1'] b1.3 <- d.ext[1:500,'chain:3','beta1'] b1 <- data.frame(chain1=b1.1, chain2=b1.2, chain3=b1.3) b1.melt <- melt(b1, id=c(), variable.name="chain") #ggplotで扱いやすいようにb1を再形成 b1.melt <- data.frame(b1.melt, iter=1:500) #サンプリング過程 p <- ggplot(b1.melt, aes(x=iter, y=value, group=chain, color=chain)) p <- p + geom_line(size=0.1) p <- p + labs(title="beta1のサンプリング過程", x="Iterations", y="") #事後分布 g <- ggplot(b1.melt, aes(x=value)) g <- g + geom_density() + theme_bw() g <- g + labs(title="beta1の事後分布", x="", y="") #サンプリング過程と事後分布を並べて表示 grid.newpage() pushViewport(viewport(layout=grid.layout(1, 2))) print(p, vp=viewport(layout.pos.row=1, layout.pos.col=1)) print(g, vp=viewport(layout.pos.row=1, layout.pos.col=2))
$\ \beta_2\ $に関しても全く同様にして図示できます。