RStan使ってMCMCの計算をやってみた(2)
前回は久保拓弥先生の「データ分析のための統計モデリング入門」(通称緑本)の9章を、Rstanを使って計算しました。
RStan使ってMCMCの計算をやってみた(1) - データサイエンスとプログラミングを学ぶ
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者:久保 拓弥
- 発売日: 2012/05/19
- メディア: 単行本
- 個体差を組み込んだ階層ベイズモデル
- 個体差と場所差を組み込んだ階層ベイズモデル
の2つの例題があります。今回は前者をやってみます。
MCサンプリング
まずはStanのコードです。以下のコードをkubo10_1.stanという名前で保存します。
data { int<lower=0> N; //データ数 int<lower=0> y[N]; //生存種子数 } parameters { real beta; //切片 real r[N]; //個体差 real<lower=0> s; //超パラメータ } transformed parameters { real q[N]; for (i in 1:N) { q[i] <- inv_logit(beta + r[i]); //生存確率 } } model { for (i in 1:N) { y[i] ~ binomial(8, q[i]); //二項分布 } beta ~ normal(0, 100); //無情報事前分布 r ~ normal(0, s); //階層事前分布 s ~ uniform(0, 1.0e+4); //超事前分布 }
生存種子数$\ y \ $は二項分布に従うとします。リンク関数はロジットリンク関数、線形予測子は$\ \beta + r_i\ $です。$\ r_i\ $が個体差を表しています。次に、以下のコードでStanに渡すデータをlistとして作成して、StanにMCサンプリングを実行させます。データは前回と同様に著者のサイト
http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html
からダウンロードできます。
library(rstan) d <- read.csv("data7a.csv") d.dat <- list(N=dim(d)[1], y=d$y) d.fit <- stan(file='kubo10_1.stan', data=d.dat, iter=10100, warmup=100, thin=10, chains=3)
前回よりも計算に時間がかかります。切片$\ \beta\ $と超パラメータ$\ s\ $のサンプリング結果を簡易的に見たければ
print(d.fit, pars=c("beta", "s"), digits_summary=3) plot(d.fit, pars=c("beta", "s")) traceplot(d.fit, pars=c("beta", "s"))
あたりでいいでしょう。
> print(d.fit,pars=c("beta", "s"), digits_summary=3) Inference for Stan model: kubo10_1. 3 chains, each with iter=10100; warmup=100; thin=10; post-warmup draws per chain=1000, total post-warmup draws=3000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat beta 0.055 0.007 0.343 -0.609 -0.172 0.051 0.283 0.733 2781 1 s 3.031 0.007 0.366 2.391 2.772 3.006 3.260 3.824 2980 1 Samples were drawn using NUTS(diag_e) at Mon Mar 10 07:14:43 2014. 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).
結果の図示
切片$\ \beta\ $と超パラメータ$\ s\ $の事後分布を図示します。extract関数でサンプリング結果を抽出して、ggplotで図示します。
library(ggplot2) library(grid) #サンプリング結果の抽出 d.ext <- extract(d.fit, permuted=T) beta <- data.frame(val=d.ext$beta) s <- data.frame(val=d.ext$s) #事後分布を並べて表示 p1 <- ggplot(beta, aes(x=val)) + geom_density() + theme_bw() + labs(title="betaの事後分布", x="beta") p2 <- ggplot(s, aes(x=val)) + geom_density() + theme_bw() + labs(title="sの事後分布", x="s") grid.newpage() pushViewport(viewport(layout=grid.layout(1, 2))) print(p1, vp=viewport(layout.pos.row=1, layout.pos.col=1)) print(p2, vp=viewport(layout.pos.row=1, layout.pos.col=2))
次に、231ページの図10.4を描画してみます。生存種子数$\ y\ $の確率分布は、二項分布$\ p(y|\beta,r)\ $と正規分布$\ p(r|s)\ $を使って$$p(y|\beta, s) = \int_{-\infty}^\infty p(y|\beta,r)p(r|s) dr$$と表せます。これを計算するにあたって、先ほどの著者のサイトに作図用Rコードがあるので、参考にさせてもらいました。
logistic <- function(z) 1/(1+exp(-z)) #ロジスティック関数 n <- nrow(d) #データ数 size <- 8 #各個体から採集する種子数 #被積分関数 gaussian.binom <- function(r, y, size, beta, sd) dbinom(y, size, logistic(beta + r))*dnorm(r, 0, sd) #積分を実行する関数 integrate.gaussian.binom <- function(y, size, beta, sd) sapply( y, function(y) integrate( f=gaussian.binom, lower=-sd*10, upper=sd*10, # for gaussian.binom y=y, size=size, beta=beta, sd=sd )$value ) #MCサンプリングで得た各beta,sで積分を実行 suv.prob <- sapply( 1:nrow(d.ext$beta), function(i) { if (i %% 100 == 0) cat(".") integrate.gaussian.binom( 0:size, size, beta=beta$val[i], sd=s$val[i] ) } ) #積分して求まった確率分布を使って生存種子数をサンプリング sample <- apply(suv.prob, 2, function(prob) summary( factor( sample(0:size, n, replace=TRUE, prob=prob), levels = 0:size ) ) ) #結果の図示 #予測分布の中央値 median <- data.frame(suv=0:size, number=apply(suv.prob*n, 1, median)) p <- ggplot(median, aes(x=suv,y=number)) + geom_line() + geom_point(shape=1, size=4) + theme_bw() #観測値 obs <- data.frame(suv=0:size, number=summary(as.factor(d$y))) p <- p + geom_point(data=obs, aes(x=suv, y=number), size=4) #予測分布の95%区間 pre.val <- apply(sample, 1, quantile, probs=c(0.025, 0.975)) pre.interval <- data.frame(suv=c(0:size, size:0), number=c(pre.val[1,], rev(pre.val[2,]))) p <- p + geom_polygon(data=pre.interval, aes(x=suv, y=number), alpha=0.2) p <- p + labs(x="生存種子数", y="個体数") p
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\ $に関しても全く同様にして図示できます。