Okada(2004) Bayes Factor model

BayesFactorTest.R

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)
fit