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

朝倉書店の「基礎からのベイズ統計学」6章末問題3)をStanを使って解きます。30個の去年のカラオケ大会データから0.90分位数のEAP (Expected A Posteriori)推定値と平均87点、標準偏差5点の人が本選に進める確率を求めます。

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

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

# karaoke score in last year contest
karaoke = [75, 82, 77, 86, 98, 91, 85, 84, 82, 79,
           88, 79, 84, 87, 69, 93, 84, 82, 89, 78,
           90, 74, 75, 84, 89, 81, 74, 79, 82, 84]
stan_data = {
  "n" : len(karaoke),
  "x" : karaoke
}

stan_code = """
  data {
    int<lower=0> n;       # number of data
    real<lower=0> x[n];   # vector of karaoke score
  }
  parameters {
    real<lower=0> mu;     # mean
    real<lower=0> stdev;  # standard deviation
  }
  model {
    x ~ normal(mu, stdev);  # normal distribution model
  }
  generated quantities {
    real<lower=0> borderline;             # 90% percentile
    real<lower=0, upper=1> qualifying;    # probability of qualifying karaoke contest
    borderline <- mu + 1.282 * stdev;
    qualifying <- 1 - normal_cdf(borderline, 87, 5);
  }
"""

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
mu          82.79  6.1e-3   1.22   80.4  81.97  82.79  83.61  85.18  40000    1.0
stdev        6.64  4.6e-3   0.92   5.13   5.99   6.54   7.18   8.72  40000    1.0
borderline   91.3  8.5e-3   1.71  88.36  90.11  91.17  92.33  95.07  40000    1.0
qualifying   0.21  4.4e-4   0.09   0.05   0.14    0.2   0.27   0.39  40000    1.0
lp__       -64.73  5.0e-3   1.01 -67.42 -65.14 -64.43 -64.01 -63.73  40000    1.0

0.90分位数EAP推定値は91.3点となりました。事後標準偏差は1.71点です。95%確信区間は88.36点から95.07点です。
平均87点、標準偏差5点の人が本選に進める確率は21%です。