3
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

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

Last updated at Posted at 2019-04-21

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

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

3
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
3
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?