Simulation Based Calibration

Stan Development Team

2022-11-26

以下は、ベータ二項モデルに対するStanのプログラムである。

data {
  int<lower = 1> N;
  real<lower = 0> a;
  real<lower = 0> b;
}
transformed data { // these adhere to the conventions above
  real pi_ = beta_rng(a, b);
  int y = binomial_rng(N, pi_);
}
parameters {
  real<lower = 0, upper = 1> pi;
}
model {
  target += beta_lpdf(pi | a, b);
  target += binomial_lpmf(y | N, pi);
}
generated quantities { // these adhere to the conventions above
  int y_ = y;
  vector[1] pars_;
  int ranks_[1] = {pi > pi_};
  vector[N] log_lik;
  pars_[1] = pi_;
  for (n in 1:y) log_lik[n] = bernoulli_lpmf(1 | pi);
  for (n in (y + 1):N) log_lik[n] = bernoulli_lpmf(0 | pi);
}

以下の規約を遵守していることに注目してみよう。

• 未知のパラメータの実現値は、transformed data ブロックで描かれ、pi_ のようにアンダースコアが後置される。 これらは、parameters ブロックで宣言された対応するシンボルによって推定される「真の」パラメータとみなされる。これらは、pi のように、末尾のアンダースコアを除いて同じ名前を持っている。

• 未知のパラメータの実測値は、transformed data ブロックの事前予測分布から描画する際に条件付けされる。この場合、int y = binomial_rng(N, pi_); となる。混乱を避けるため、y にはトレーニング用のアンダースコアがない。

• 未知のパラメータの実測値は、generated quantities ブロックの vector にコピーされ、pars_ という名前になる。

• 事前予測分布からの実現値は、generated quantities ブロックの `y_ という名前のオブジェクト（同じ型）にコピーされる。これはオプションである。

• generated quantities ブロックには、ranks_ という名前の整数配列があり、その値は、事後分布からのパラメータの実現値が対応する「真のランク」を超えるかどうかに応じて、0 または 1 だけであるが、後でランクを再構築（間引き）するために使用できる。 " realization, which in this case is ranks_ [1] = {pi > pi_}; . These are not actually" しかし、その後、（間引かれた）ランクを再構築するために使用することができる。

• generated quantities ブロックは log_lik というベクタを含み、その値は各観測による対数尤度への寄与を表する。この場合、「オブザベーション」は、二項尤度に集約される暗黙の成功と失敗である。これはオプションであるが、事後分布が特定のオブザベーションに敏感であるかどうかを示すPareto k形状パラメータ推定値の計算を容易にするものである。

上記が beta_binomial というコード stanmodel にコンパイルされたと仮定すると、sbc 関数を呼び出すことができる。

output <- sbc(beta_binomial, data = list(N = 10, a = 1, b = 1), M = 500, refresh = 0)

この時点で、次に

print(output)

## 0 chains had divergent transitions after warmup

plot(output, bins = 10) # it is best to specify the bins argument yourself