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

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

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

前回は久保拓弥先生の「データ分析のための統計モデリング入門」(通称緑本)の9章を、Rstanを使って計算しました。
RStan使ってMCMCの計算をやってみた(1) - データサイエンスとプログラミングを学ぶ

今回は緑本の10章をやってみます。10章は

  • 個体差を組み込んだ階層ベイズモデル
  • 個体差と場所差を組み込んだ階層ベイズモデル

の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).

f:id:physyuki:20140310072748p:plain
f:id:physyuki:20140310072757p:plain

結果の図示

切片$\ \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))

f:id:physyuki:20140310162140p:plain
次に、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

f:id:physyuki:20140310162056p:plain