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

朝倉書店の「基礎からのベイズ統計学」６章末問題1)をStanを使って解きます。１２個のデータから平均のEAP (Expected A Posteriori)推定値と平均が70,000円を超える確率と75,000円を超える確率を求めます。

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

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

uriage = [76230, 73550, 80750, 71500, 75420, 74840,
          71580, 76920, 68450, 76990, 64070, 76200]
stan_data = {
  "n" : len(uriage),
  "x" : uriage
}

stan_code = """
  data {
    int<lower=0> n;         # number of data
    real<lower=0> x[n];   # vector of uriage
  }
  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, upper=1> over0;   # probability of EAP over 70k JPY
    real<lower=0, upper=1> over1;   # probability of EAP over 75k JPY
    over0 <- step(mu-70000);
    over1 <- step(mu-75000);
  }
"""

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     7.4e4    7.53 1506.7  7.1e4  7.3e4  7.4e4  7.5e4  7.7e4  40000    1.0
stdev 5041.7    6.26 1252.8 3274.0 4166.7 4829.6 5655.0 8089.9  40000    1.0
over0   0.99  4.2e-4   0.08    1.0    1.0    1.0    1.0    1.0  40000    1.0
over1   0.21  2.0e-3   0.41    0.0    0.0    0.0    0.0    1.0  40000    1.0
lp__  -87.76  5.4e-3   1.08 -90.63 -88.18 -87.43 -86.99 -86.71  40000    1.0

EAP推定値は74,000円となりました。事後標準偏差は1,506円です。95%確信区間は71,000円から77,000円です。
70,000円を超える確率は99%です。
75,000円を超える確率は21%です。