岩波データサイエンスvol1のいくつかの例題をStanでやってみた

岩波データサイエンスは従来の書籍とは異なり、以下のサポートページの異様な充実がウリの一つです。 https://sites.google.com/site/iwanamidatascience/

岩波データサイエンス Vol.1

また、StanとRStanのはじめの一歩と言える使い方が載っています。また、作図を含めたソースコードがGithubのリポジトリにアップされていますので参考になるでしょう。

今回は久保さんや伊東さんの記事にある例題をいくつかStanで書きなおしてみました。

「階層ベイズ　最初の一歩」の例題から

モデル1のベイズ版

data {
   int N;
   real Mean_Y[N];
   real X[N];
   real Age[N];
}

parameters {
   real beta[3];
   real<lower=0> sigma;
}

transformed parameters {
   real mu[N];
   for (i in 1:N)
      mu[i] <- beta[1] + (beta[2] + beta[3] * X[i]) * Age[i];
}

model {
   for (i in 1:N)
      Mean_Y[i] ~ normal(mu[i], sigma);
}

気をつけるべきは本にもちょろっと書いてありますが、21行目のnormal関数です。JAGSなどのBUGS言語ではdnorm(平均, 精度)で渡すのに対し、Stanではnormal(平均, 標準偏差)で渡します。また、本にも書きましたがStanでは事前分布の記述をしなければ自動で十分に広い一様分布が使われるので、積極的に使っています。

これを実行するRコードは以下のとおり。

library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

load('data.RData')

data <- list(
   N=nrow(d),
   Mean_Y=d$mean.Y,
   X=d$X,
   Age=d$Age
)

stanmodel <- stan_model(file='model/model1.stan')
fit <- sampling(stanmodel, data=data, iter=1000, warmup=200, seed=1)
save.image('post1.RData')

• 2～3行目は並列計算させる時にこれらを書きます。うまくいくまではコメントアウトして実行しましょう。

モデル2の場合

data {
   int N_r;
   int N_pref;
   real Mean_Y[N_r, N_pref];
   real Sigma_Y[N_r, N_pref];
   real X[N_r, N_pref];
   real Age[N_r, N_pref];
}

parameters {
   real beta[3];
   real s_r[N_r];
   real r[N_r, N_pref];
}

transformed parameters {
   real mu[N_r, N_pref];
   for (i in 1:N_r)
      for (j in 1:N_pref)
         mu[i, j] <- beta[1] + r[1, j] + (beta[2] + r[2, j] + beta[3] * X[i, j]) * Age[i, j];
}

model {
   for (i in 1:N_r)
      for (j in 1:N_pref)
         Mean_Y[i, j] ~ normal(mu[i, j], Sigma_Y[i, j]);

   for (i in 1:N_r)
      for (j in 1:N_pref)
         r[i, j] ~ normal(0, s_r[i]);
}

• 4～7行目:データを(測定時点, 県)の2次元配列にして扱っています。
• 30行目:階層ベイズになっている部分です。ゆるくしばっています。

キックするRコードは以下です。

library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

load('data.RData')
to2d <- function(v) t(matrix(v, N.pref, 2))

data <- list(
   N_r=2,
   N_pref=length(unique(d$pref)),
   Mean_Y=to2d(d$mean.Y),
   Sigma_Y=to2d(d$sd.Y/sqrt(d$N)),
   X=to2d(d$X),
   Age=to2d(d$Age)
)

stanmodel <- stan_model(file='model/model2.stan')
fit <- sampling(stanmodel, data=data, iter=1000, warmup=200, seed=1)
save.image('post2.RData')

• 6行目:この関数は久保さんのソースコードにあったものを使っています。

「時系列・空間データのモデリング」の年輪データの例題から

data {
   int<lower=0> N;
   real<lower=0> Y[N];
}

parameters {
   real alpha[N];
   real<lower=0> s_Y;
   real<lower=0> s_a;
}

model {
   for(i in 1:N)
      Y[i] ~ normal(alpha[i], s_Y);
      # Y[i] ~ lognormal(alpha[i], s_Y);

   for(i in 3:N)
      alpha[i] ~ normal(2*alpha[i-1] - alpha[i-2], s_a);
}

• 14行目:観測モデル。
• 18行目:システムモデル。「Stan入門」の空間構造のあるベイズモデルであえてincrement_log_prob()関数を使っているのはこちらのシステムモデルの方。
• 15行目:伊東さんの記事にある対数正規分布を使う場合は14行目の代わりにこちらを使います。こういうのを探す時はマニュアルで「lognormal」とかで検索します。

キックするRコードは以下です。

library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

X <- 1961:1990
Y <- c(4.71, 7.70, 7.97, 8.35, 5.70,
       7.33, 3.10, 4.98, 3.75, 3.35,
       1.84, 3.28, 2.77, 2.72, 2.54,
       3.23, 2.45, 1.90, 2.56, 2.12,
       1.78, 3.18, 2.64, 1.86, 1.69,
       0.81, 1.02, 1.40, 1.31, 1.57)

data <- list(N=length(X), Y=Y)
stanmodel <- stan_model(file='model/model3.stan')
fit <- sampling(stanmodel, data=data, iter=1700, warmup=200, thin=3, seed=1)
save.image('post3.RData')

thin=1だとちょっと自己相関が高くてn_effが稼げずRhatが1.1ぐらいだったのでthin=3にしました。

よくひっかかる注意点は書籍とサポートページにちゃんと書いてあるので1行1句見落としのないように読むことをオススメします。