Edited at

緑本の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)

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