ベイズ統計： pystanを使う 〜 正規分布の独立な2群の平均値差に関する推測

朝倉書店の「基礎からのベイズ統計学」6章末問題4)をStanを使って解きます。それぞれ20個のデータをもつ2群において、分散が同一であるという仮定のもとで、平均値に差があるか検討します。

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

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

# karaoke score in last year contest
experiment = [30.86, 29.75, 31.55, 32.29, 29.90, 31.71, 31.35, 29.03, 30.37, 31.55,
              29.26, 32.29, 29.90, 30.18, 30.72, 32.28, 30.72, 29.90, 31.55, 31.55]
control = [31.36, 33.34, 33.16, 31.36, 36.19, 29.80, 31.11, 35.23, 31.36, 31.27,
           31.63, 31.63, 32.00, 31.11, 31.63, 31.36, 31.81, 31.63, 29.21, 33.37]
stan_data = {
  "exN" : len(experiment),
  "ex" : experiment,
  "conN" : len(control),
  "con" : control
}

stan_code = """
  data {
    int<lower=0> exN;
    real<lower=0> ex[exN];
    int<lower=0> conN;
    real<lower=0> con[conN];
  }
  parameters {
    real<lower=0> mu_ex;    # mean for experiment
    real<lower=0> mu_con;   # mean for control
    real<lower=0> stdev;    # standard deviation
  }
  model {
    ex ~ normal(mu_ex, stdev);    # normal distribution model
    con ~ normal(mu_con, stdev);  # normal distribution model
  }
  generated quantities {
    real delta;
    real<lower=0> over0;
    real<lower=0> over1;
    delta <- mu_con - mu_ex;
    over0 <- step(delta);
    over1 <- step(delta-1);
  }
"""

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_ex   30.84  1.6e-3   0.32  30.22  30.63  30.83  31.05  31.46  40000    1.0
mu_con  31.98  1.6e-3   0.31  31.36  31.77  31.97  32.18  32.59  40000    1.0
stdev     1.4  8.3e-4   0.17   1.12   1.28   1.38    1.5   1.77  40000    1.0
delta    1.14  2.2e-3   0.45   0.27   0.84   1.14   1.43   2.02  40000    1.0
over0    0.99  3.8e-4   0.08    1.0    1.0    1.0    1.0    1.0  40000    1.0
over1    0.62  2.4e-3   0.49    0.0    0.0    1.0    1.0    1.0  40000    1.0
lp__   -25.44  6.3e-3   1.26  -28.7 -26.02 -25.12 -24.52 -24.01  40000    1.0

平均の差(delta)のEAP推定値は1.14となりました。事後標準偏差は0.45です。95%確信区間は0.27から2.02です。差がある確率(over0)は99%で、ほぼ差があると見ていいようです。差が1以上ある確率(over1)は62%です。