データサイエンスとプログラミングを学ぶ

新米データアナリストです。

RStan使ってMCMCの計算をやってみた(1)

久保拓弥先生の「データ分析のための統計モデリング入門」(通称緑本)の9章をトレースしてみました。緑本ではWinBUGSを使用していますが、ここではRStanを使ってStanでやってみます。

インストール

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)

で以下の図が得られます。
f:id:physyuki:20140215222838p:plain
f:id:physyuki:20140215222744p:plain

結果の図示

以下の要領で$\ \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\ $に関しても全く同様にして図示できます。
f:id:physyuki:20140215222255p:plain
f:id:physyuki:20140215222303p:plain