習作；ベイズファクター

岡田(2014)のベイズファクターの計算をStanでやってみた件。

これについては，次のスライドに詳しい。

ただ，途中のStanコードがちょっとわかりにくいのと，計算をstan内部でgenerated quantitiesを使ってやったらいいんじゃないかと思ったので，そのように書き直して見た習作。

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
model <- '
data{
  int<lower=0> N1;
  int<lower=0> N2;
  int<lower=0> N3;
  real y1[N1];
  real y2[N2];
  real y3[N3];
}
parameters{
  real mu1;
  real mu2;
  real mu3;
  real<lower=0> sigma;
}
model{
  y1 ~ normal(mu1,sigma);
  y2 ~ normal(mu2,sigma);
  y3 ~ normal(mu3,sigma);
  mu1 ~ normal(0,100);
  mu2 ~ normal(0,100);
  mu3 ~ normal(0,100);
  sigma ~ inv_gamma(0.01,0.01);
}
generated quantities{
  real f1a;
  real f1b;
  real f2a;
  real f2b;
  real f1;
  real f2;
  real BF1;
  real BF2;
  real PMP1;
  real PMP2;
  f1a = (mu2>mu1? 1:0);
  f1b = (mu3>mu2? 1:0);
  f1 = f1a * f1b;
  f2a = (mu2>mu1? 1:0);
  f2b = (mu3>mu1? 1:0);
  f2 = f2a * f2b;
  BF1 = f1 / (1.0/6.0);
  BF2 = f2 / (1.0/3.0);
  PMP1 = BF1 / (BF1 + BF2 + 1);
  PMP2 = BF2 / (BF1 + BF2 + 1);
}
'

library(Lock5Data)
dat  <- subset(LightatNight,select=c("Light","BMGain"))
y1=subset(dat,dat$Light=="LD")
y2=subset(dat,dat$Light=="DM")
y3=subset(dat,dat$Light=="LL")
datastan <- list(N1=nrow(y1),N2=nrow(y2),N3=nrow(y3),
                 y1=y1$BMGain,y2=y2$BMGain,y3=y3$BMGain)
fit <- stan(model_code=model,data=datastan,
                iter=21000,warmup=1000,chains=5)

推定結果は論文とほぼ一致した。BF1,BF2,PMP1,PMP2を参照。

Inference for Stan model: 635fb9451ae02a2b78fba0033bf61ec5.
5 chains, each with iter=21000; warmup=1000; thin=1; 
post-warmup draws per chain=20000, total post-warmup draws=1e+05.

        mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
mu1     5.93    0.00 0.91   4.12   5.34   5.93   6.52   7.74 81771    1
mu2     7.86    0.00 0.82   6.24   7.33   7.85   8.39   9.47 87166    1
mu3    11.01    0.00 0.86   9.30  10.45  11.00  11.57  12.71 88478    1
sigma   2.68    0.00 0.39   2.05   2.40   2.63   2.90   3.57 72728    1
f1a     0.94    0.00 0.23   0.00   1.00   1.00   1.00   1.00 76653    1
f1b     0.99    0.00 0.08   1.00   1.00   1.00   1.00   1.00 80118    1
f2a     0.94    0.00 0.23   0.00   1.00   1.00   1.00   1.00 76653    1
f2b     1.00    0.00 0.01   1.00   1.00   1.00   1.00   1.00 89506    1
f1      0.94    0.00 0.24   0.00   1.00   1.00   1.00   1.00 72399    1
f2      0.94    0.00 0.23   0.00   1.00   1.00   1.00   1.00 76683    1
BF1     5.63    0.01 1.45   0.00   6.00   6.00   6.00   6.00 72399    1
BF2     2.83    0.00 0.69   0.00   3.00   3.00   3.00   3.00 76683    1
PMP1    0.56    0.00 0.14   0.00   0.60   0.60   0.60   0.60 72399    1
PMP2    0.29    0.00 0.08   0.00   0.30   0.30   0.30   0.30 81240    1
lp__  -39.67    0.01 1.46 -43.41 -40.36 -39.32 -38.62 -37.93 44375    1

Samples were drawn using NUTS(diag_e) at Tue Dec 13 07:14:28 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).