緑本のBUGSをStanで書いてみた

緑本のBugsをRstanで書いてみた

みんな大好き緑本（データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)こちらに出てくる、BUGSのコードをRstanで実装したので、今日はそれのご紹介。

「10.5 個体差＋場所差の階層ベイズモデル」をRstanで書きました。

data {
    int N;
    int P; #Pot数（1～10に変換）
    int F; #施肥処理の有無（0-1に変換）
    int y[N];　＃個体差
    int pot[N];　＃Pot差
    int f[N];　＃施肥処理差
}

parameters{
    real beta1;
    real beta2;
    real id_est[N];　＃個体差
    real pot_est[P];　＃Pot差
    real<lower = 0> s1;
    real<lower = 0> s2;
}    

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

model {
    for(i in 1:N){
      y[i] ~ poisson(lambda[i]);
    }
    beta1 ~ normal(0, 100);
    beta2 ~ normal(0, 100);
    for(i in 1:N){
      id_est[i] ~ normal(0, s1);
    }
    for(p in 1:P){
      pot_est[p] ~ normal(0, s2);
    }
    s1 ~ normal(0, 1000);
    s2 ~ normal(0, 1000);
}

書籍のBUGSのコードでは、lambdaはlog(lambda[i] <- beta1 + beta2*F[i] + r[i] + rp[Pot[i]]としていますが、このコードではtransformed parametersのところで、exponentialに式変形しています。

書籍では、「beta2の（周辺）事後分布の95%区間は-2.45%から0.75ぐらい」とあり、このStanで再現した場合のbeta2の95%区間は-2.49 ~ 0.74となりましたので、再現できているかと思います。(iteration 2000, burn in 400, chain 3)

もし間違ってるよ、等あればご連絡ください。。