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