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

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

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

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