ベイズ統計： rstanを使う 〜 正規分布の対応のある2群の平均値差に関する推測

朝倉書店の「基礎からのベイズ統計学」6章末問題5)をStanを使って解きます。あるダイエット方法を20人が試して、ダイエット前後に平均値に差があるか検討します。平均値差のEAP推定値と2群の平均値差が0kgより大きいという研究仮説と、2kgより大きいという研究仮説が正しい確率を求めます。加えて相関係数が0.7より大きい確率も求めます。

今回はpystanではなく、rstanを使います。pystanでは実行中にmemory errorが発生してプログラムが完了しませんでした。stan実行時のchainsパラメータを小さくしましたが変わりませんでした。

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

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

setwd("~/R/repository/bayes_hmc")

stan_file <- "ch0605.stan"
x=read.csv(file="./ch0605.data", header=T)
stan_data <- list(n=nrow(x), x=x)
par<-c("mu","stdev","rho","stdevsq","sigma",
       "delta","over0","over2", "rho_over0", "rho_over07")

fit <- stan(file=stan_file,
            data=stan_data,
            iter=11000,
            warmup=1000,
            thin=1,
            chains=4)

print(fit, digits_summary=3)
plot(fit, pars=par)
traceplot(fit, inc_warmup=F, pars=par)

stanのコード(ch0605.stanの内容)です。

  data {
    int<lower=0>  n;
    vector[2]     x[n];
  }
  parameters {
    vector[2]               mu;     # mean
    vector<lower=0>[2]      stdev;  # standard deviation
    real<lower=-1, upper=1> rho;    # correlation coefficient
  }
  transformed parameters {
    vector<lower=0>[2] stdevsq;   # variance
    matrix[2, 2] sigma;
    stdevsq[1] <- pow(stdev[1], 2);
    stdevsq[2] <- pow(stdev[2], 2);
    sigma[1, 1] <- stdevsq[1];
    sigma[2, 2] <- stdevsq[2];
    sigma[1, 2] <- stdev[1] * stdev[2] * rho;
    sigma[2, 1] <- stdev[1] * stdev[2] * rho;
  }
  model {
    x ~ multi_normal(mu, sigma);    # multi normal distribution model
  }
  generated quantities {
    real delta;
    real<lower=0> over0;
    real<lower=0> over2;
    real<lower=0> rho_over0;
    real<lower=0> rho_over07;
    delta <- mu[1] - mu[2];
    over0 <- step(delta);
    over2 <- step(delta-2);
    rho_over0 <- step(rho);
    rho_over07 <- step(rho-0.7);
  }

データです。

prior, post
53.1, 51.3
51.5, 48.2
45.5, 49.6
55.5, 59.6
49.6, 44.2
50.1, 47.6
59.2, 54.9
54.7, 56.5
53.0, 48.4
48.6, 50.6
55.3, 53.6
54.6, 57.5
51.7, 52.0
48.6, 46.9
56.4, 56.8
58.9, 50.0
53.3, 53.8
42.4, 38.3
51.9, 52.6
39.1, 41.0

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

              mean se_mean     sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
mu[1]       51.625   0.011  1.228  49.197  50.818  51.627  52.423  54.074 13241    1
mu[2]       50.638   0.011  1.312  48.015  49.791  50.640  51.489  53.197 13317    1
stdev[1]     5.434   0.009  0.950   3.956   4.758   5.306   5.962   7.659 12272    1
stdev[2]     5.817   0.009  1.018   4.225   5.095   5.687   6.385   8.192 12752    1
rho          0.744   0.001  0.108   0.482   0.687   0.763   0.823   0.901 14120    1
stdevsq[1]  30.433   0.103 11.255  15.648  22.639  28.153  35.548  58.657 11871    1
stdevsq[2]  34.874   0.117 12.937  17.853  25.963  32.338  40.774  67.103 12215    1
sigma[1,1]  30.433   0.103 11.255  15.648  22.639  28.153  35.548  58.657 11871    1
sigma[1,2]  24.487   0.106 10.469  10.325  17.305  22.467  29.341  50.625  9806    1
sigma[2,1]  24.487   0.106 10.469  10.325  17.305  22.467  29.341  50.625  9806    1
sigma[2,2]  34.874   0.117 12.937  17.853  25.963  32.338  40.774  67.103 12215    1
delta        0.987   0.005  0.899  -0.784   0.400   0.979   1.570   2.765 36824    1
over0        0.872   0.002  0.334   0.000   1.000   1.000   1.000   1.000 29392    1
over2        0.124   0.002  0.329   0.000   0.000   0.000   0.000   1.000 28841    1
rho_over0    1.000   0.000  0.005   1.000   1.000   1.000   1.000   1.000 40000    1
rho_over07   0.714   0.003  0.452   0.000   0.000   1.000   1.000   1.000 18164    1
lp__       -76.711   0.015  1.666 -80.816 -77.580 -76.375 -75.479 -74.500 12078    1

平均値差(delta)は0.987で差があるように見えます。平均値差が0kgより大きい確率(over0)は87.2%で差があることは確かなようです。平均値差が2kgより大きい確率(over2)は12.4%となり、非常に効果があるダイエット方法とは言えません。
また、相関係数が0.7より大きい確率(rho_over07)は71.4%ですので、中程度から強い相関はありそうです。