RStan 2.0.0 をためす

RStanは表面上はあまりかわっていないような。

data(Nile)を状態空間モデルで平滑化してみる。

Stanコード(ar1.stan)

data {
  int<lower=0> n;
  real y[n];
}
parameters {
  real x[n];
  real<lower=0> sigma_x;
  real<lower=0> sigma_y;
  real y1;
}
model {
  x[1] ~ normal(y1, sigma_x);
  y[1] ~ normal(x[1], sigma_y);

  for (i in 2:n) {
    x[i] ~ normal(x[i - 1], sigma_x);
    y[i] ~ normal(x[i], sigma_y);
  }

  y1 ~ normal(0.0, 1.0e+4);
  sigma_x ~ uniform(0.0, 1.0e+4);
  sigma_y ~ uniform(0.0, 1.0e+4);
}

Rコード

library(rstan)

data(Nile)

nile.stan <- stan("ar1.stan",
                  data = list(y = as.vector(Nile), n = length(Nile), y1 = Nile[1]),
                  pars = c("x", "sigma_x", "sigma_y"),
                  chains = 4, iter = 4000, warmup = 1000, thin = 3)
st <- start(Nile)[1]
ed <- end(Nile)[1]
len <- ed - st + 1
plot(Nile, xlab = "Year", ylab = "Flow", las = 1)
x <- apply(as.matrix(nile.stan), 2, mean)
loc.x1 <- match("x[1]", labels(x))
x <- x[loc.x1:len]
x.95 <- apply(as.matrix(nile.stan), 2, 
              quantile, probs = c(0.025, 0.975))[, loc.x1:len]
lines(st:ed, x, col = 2)
lines(st:ed, x.95[1, ], col = 2, lt = 2)
lines(st:ed, x.95[2, ], col = 2, lt = 2)

グラフ

mclapply()による並列化でRStan実行。

library(parallel)

nile.model <- stan_model("ar1.stan")
seed <- c(1234, 12345, 123456, 1234567)
nile.fit <- mclapply(1:4, function(i) 
                            sampling(nile.model,
                                     data = list(y = as.vector(Nile),
                                                 n = length(Nile),
                                                 y1 = Nile[1]),
                                     pars = c("x", "sigma_x", "sigma_y"),
                                     chain_id = i, seed = seed[i],
                                     chains = 1, iter = 4000, warmup = 1000,
                                     thin = 3),
                     mc.cores = 4)
nile.stan2 <- sflist2stanfit(nile.fit)
plot(nile.stan2)

※ 2013-10-20: "chain_id = i"をパラメータに追加。

結果

タグ：STAn RStan