ベイズ統計： pystanを使う 〜 正規分布の分散に関する推測

朝倉書店の「基礎からのベイズ統計学」６章末問題2)をStanを使って解きます。20個のデータから分散のEAP (Expected A Posteriori)推定値と分散が1.5以下である確率、1.0以下である確率を求めます。

コードはこちらになります。

# -*- coding: utf-8 -*-
import pystan

juryo = [128.27, 130.00, 130.68, 128.82, 130.08, 129.24, 129.82,
         131.76, 128.19, 131.01, 128.81, 131.53, 129.36, 130.82,
         130.89, 130.35, 130.18, 131.00, 131.03, 128.84]
stan_data = {
  "n" : len(juryo),
  "x" : juryo
}

stan_code = """
  data {
    int<lower=0> n;       # number of data
    real<lower=0> x[n];   # vector of juryo
  }
  parameters {
    real<lower=0> stdev;  # standard deviation
  }
  transformed parameters {
    real<lower=0> stdevsq;  # variance
    stdevsq <- pow(stdev, 2);
  }
  model {
    x ~ normal(130.0, stdev);  # normal distribution model
  }
  generated quantities {
    real<lower=0, upper=1> under0;   # probability of EAP under 1.5
    real<lower=0, upper=1> under1;   # probability of EAP under 1.0
    under0 <- step(1.5-stdevsq);
    under1 <- step(1.0-stdevsq);
  }
"""

def main():
  fit = pystan.stan(model_code=stan_code, 
                    data=stan_data,
                    iter=11000,
                    warmup=1000,
                    thin=1,
                    chains=4)
  print(fit)
  fit.plot()

if __name__ == "__main__":
  main()

実行結果はこうなりました。

          mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
stdev     1.12  9.5e-4   0.19   0.82   0.98   1.09   1.23   1.56  40000    1.0
stdevsq   1.29  2.3e-3   0.46   0.67   0.97    1.2   1.51   2.43  40000    1.0
under0    0.75  2.2e-3   0.44    0.0    0.0    1.0    1.0    1.0  40000    1.0
under1    0.29  2.3e-3   0.45    0.0    0.0    0.0    1.0    1.0  40000    1.0
lp__    -11.41  3.5e-3    0.7 -13.42 -11.57 -11.14 -10.97 -10.92  40000    1.0

分散のEAP推定値は1.29です。この分散の事後標準偏差は0.46です。
1.5以下である確率は75%、1.0以下である確率は29%となりました。