Bridge sampling

R

data

data <- read.csv("~/Williams2017/Williams2017.csv")
X1 <- data[data$C=="1",]$X
N1 <- length(X1)
X2 <- data[data$C=="2",]$X
N2 <- length(X2)

Stan

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

prior

delta_prior <- sqrt(1/2) # 0.707

data for stan

iter <- 5000
datastan <- list(N1=N1, N2=N2, X1=X1, X2=X2, delta_prior=delta_prior)
fit.null <- stan("Hypoth0bs.stan", data=datastan, iter=iter*2, warmup=iter, 
refresh=0, chains=4)
print(fit.null, pars=c("mu", "v", "s"))
fit.test <- stan("Hypoth1bs.stan", data=datastan, iter=iter*2, warmup=iter,
refresh=0, chains=4)
print(fit.test, pars=c("mu", "ma", "delta", "v", "s"))

Bayes factor

library(bridgesampling)
cat("\n** Bridge sampling **", "\n")
bs.null <- bridge_sampler(samples=fit.null, method="warp3", silent=TRUE)
bs.test <- bridge_sampler(samples=fit.test, method="warp3", silent=TRUE)
cat("\nNull Hypothesis\n")
print(bs.null)
cat("\nAlternative Hypothesis\n")
print(bs.test)
cat("\n")
bf10 <- bayes_factor(bs.test, bs.null)
print(bf10,)
cat("\n")

---

Hypoth0bs.stan

data {
  nt<lower=0> N1;
  int<lower=0> N2;
  vector[N1] X1;
  vector[N2] X2;
}
parameters {
  real mu;
  real<lower=0> v;
}
transformed parameters {
  real s = sqrt(v);
}
model {
  target += log(1/v);
  target += normal_lpdf(X1 | mu, s);
  target += normal_lpdf(X2 | mu, s);
}

---

Hypoth1bs.stan

data {
  int<lower=0> N1;
  int<lower=0> N2;
  vector[N1] X1;
  vector[N2] X2;
  real delta_prior;
}
parameters {
  real mu;
  real<lower=0> v;
  real delta;
}
transformed parameters {
  real s = sqrt(v);
  real diff = delta*s;
}
model {
  target += log(1/v);
  target += cauchy_lpdf(delta | 0, delta_prior);
  target += normal_lpdf(X1 | mu - diff/2, s);
  target += normal_lpdf(X2 | mu + diff/2, s);
}
generated quantities {
  array[2] real ma;
  ma[1] = mu - 0.5diff;
  ma[2] = mu + 0.5diff;
}

---

Output

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

mu 4.73 0.00 0.33 4.08 4.51 4.73 4.95 5.39 14259 1
v    4.25 0.01 1.05 2.68 3.51 4.09 4.83 6.75 12340 1
s    2.05 0.00 0.25 1.64 1.87 2.02 2.20 2.60 12844 1

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

mu     4.66 0.00 0.29 4.08 4.47 4.66 4.85 5.23 15012 1
ma[1] 3.78 0.00 0.42 2.94 3.50 3.79 4.07 4.61 15476 1
ma[2] 5.53 0.00 0.40 4.74 5.27 5.53 5.79 6.32 17315 1
delta  0.98 0.00 0.35 0.30 0.75 0.98 1.22 1.67 14200 1
v        3.30 0.01 0.81 2.07 2.73 3.17 3.74 5.20 12342 1
s        1.80 0.00 0.22 1.44 1.65 1.78 1.93 2.28 12575 1

Bridge sampling

Null Hypothesis
Bridge sampling estimate of the log marginal likelihood: -82.67273
Estimate obtained in 4 iteration(s) via method "warp3".

Alternative Hypothesis
Bridge sampling estimate of the log marginal likelihood: -79.3768
Estimate obtained in 4 iteration(s) via method "warp3".

Estimated Bayes factor in favor of x1 over x2: 27.00267

---

BayesFactor: R package
Posterior probability distributions