R: RStanの並列実行

RStan Getting Startedにもあるが、rstan()は、parallelパッケージとくみあわせて並列化できる。

久保本10章の例でやってみる。

library(rstan)
library(parallel)

# read data
d <- read.csv(url("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/hbm/data7a.csv"))

model <- '
  data {
    int<lower=0> N;  // sample size
    int<lower=0> Y[N];  // response variable
  }
  parameters {
    real beta;
    real r[N];
    real<lower=0> sigma;
  }
  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, sigma); // 階層事前分布
	  sigma ~ uniform(0, 1.0e+4); # 無情報事前分布
  }'

data <- list(N = nrow(d), Y = d$y)
pars <- c("beta", "sigma")
seeds <- c(123, 1234, 12345, 123456)

f1 <- stan(model_code = model, data = data, pars = pars,
           chains = 1, iter = 10)
fit.list <- mclapply(1:4,
                     function(i) stan(fit = f1, data = data, pars = pars,
                                      chains = 1, chain_id = i, seed = seeds[i],
                                      warmup = 100, iter = 10100, thin = 10))
fit2 <- sflist2stanfit(fit.list)
print(fit2)

結果

> print(fit2)
Inference for Stan model: model.
4 chains, each with iter=10100; warmup=100; thin=10; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean  sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
beta     0.0     0.0 0.3   -0.6   -0.2    0.0    0.3    0.7  2828    1
sigma    3.0     0.0 0.4    2.4    2.8    3.0    3.3    3.8  3450    1
lp__  -443.6     0.2 9.3 -462.6 -449.7 -443.2 -436.8 -426.7  2836    1

Samples were drawn using NUTS2 at Wed Jul 10 21:40:11 2013.
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).

タグ：R RStan