緑本をStanで書く

はじめに

短くて恐縮ですが, 自分用のメモです．

緑本("データ解析のための統計モデリング入門"(久保拓也著))の9章以降では既に開発が終了しているWinBUGSが使われていますが、入門を兼ねてStanで書き直してみました．

教科書ではギブスサンプリングでMCMCを実行していますが, Stanではハミルトニアンモンテカルロ法を採用しているため,パラメータの事後分布が異なる点に注意してください．

9.4 ベイズ統計モデルの事後分布の推定

• Stanファイルの作成

data{
  int N;
  int Y[N];
  real X[N];
}

parameters{
  real beta1;
  real beta2;
}

transformed parameters{
  real lambda[N];
  for(i in 1:N){
    lambda[i] = exp(beta1 + beta2 * X[i]);
  }  
}

model{
  for(i in 1:N){
    Y[i] ~ poisson(lambda[i]);
  }
  beta1 ~ normal(0, 100);
  beta2 ~ normal(0, 100);
}

• 実行方法 (サンプリングパラメータは教科書に準拠しています)

set.seed(123)
d = load.csv('d.RData') 
fit = sampling(model, data=list(N=20, X=d$x, Y=d$yiter=1600, warmup=100, thin=3, chains=3))
model = stan_model("poisson.stan")
fit

• 結果 (lambdaパラメータの記載は省略)

             mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
beta1        1.55    0.01 0.37   0.82   1.32   1.56   1.81   2.25   654    1
beta2        0.09    0.00 0.07  -0.05   0.04   0.08   0.13   0.22   661    1

教科書に記載の値 (beta1 mean: 1.973, beta2 mean: 0.082)とは異なりますが,より真の値(beta2: 0.1)に近いですね． 95%信頼区間の中にゼロが入っちゃってますが、、

10.3 階層ベイズモデルの推定・予測

• Stanファイルの作成

data{
  int N;
  int Y[N];
}

parameters{
  real beta;
  real r[N];
  real 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]);
    r[i] ~ normal(0, s);
  }
  beta ~ normal(0, 100);
  s ~ uniform(0, 10000);
}

• 結果 (一部のみ抜粋)

          mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
beta      0.03    0.01 0.34   -0.64   -0.19    0.03    0.25    0.71   621 1.01
r[1]     -3.82    0.03 1.74   -7.80   -4.77   -3.59   -2.59   -1.11  2920 1.00
r[2]     -1.19    0.02 0.88   -3.05   -1.74   -1.15   -0.60    0.43  3017 1.00
r[3]      1.99    0.02 1.09    0.13    1.22    1.90    2.68    4.42  3212 1.00
s         3.02    0.01 0.36    2.39    2.77    2.99    3.24    3.84  1200    1

教科書に数値の記載がないため正しく推定できているか怪しいですが, 図10.3に一致していそうです．

10.5

こちらの記事 を参照ください．

11章

こちらの記事 を参照ください．