Bayesian t-test using pystan

Modeling
1

Hi all,

I’ve been following Ch. 7.4 of the Bayesian Methods for Hackers book, which is about implementing a t-test (Bayesian Estimation Supersedes the t-test or BESt), and there are some examples out there (including here) on this using pymc but I’ve playing around with how one would do this using pystan. There’s a very cool vignette on this using rstan, but that takes a slightly different tack or approach.

Would anybody know if there’s a way to do this using pystan?

Take this example from pymc:

dataset_a = np.random.normal(10, 2, size=100)
dataset_b = np.random.normal(8, 3, size=100)

# the pymc part
pooled_sd = np.r_[dataset_a, dataset_b].std()
tau = 1./np.sqrt(1000.*pooled_sd) # pymc precision param, 1/sigma**2

# mu from normal
mu_a = pm.Normal('mu_a', np.r_[dataset_a, dataset_b].mean(), tau)
mu_b = pm.Normal('mu_b', np.r_[dataset_a, dataset_b].mean(), tau)

# sd from uniform
sd_a = pm.Uniform('sd_a', pooled_sd/1000., 1000.*pooled_sd)
sd_b = pm.Uniform('sd_b', pooled_sd/1000., 1000.*pooled_sd)

# nu from exponential
nu_minus_1 = pm.Exponential('nu-1', 1./29)

# putting it together
obs_a = pm.NoncentralT('obs_a', mu_a, 1.0/sd_a**2, nu_minus_1 + 1, observed=True, value=dataset_a)
obs_b = pm.NoncentralT('obs_b', mu_a, 1.0/sd_a**2, nu_minus_1 + 1, observed=True, value=dataset_b)

mcmc = pm.MCMC([obs_a, obs_b, mu_a, mu_b, sd_a, sd_b, nu_minus_1])
mcmc.sample(25000, 10000)

Then one can look at the traces of mus and sds and compare them.

Now, it seems to me, from the vignette, that the Stan equivalent would be to use…

data {
  int<lower=1> n;
  vector[n] y;
  real<lower=0> r;
}
parameters {
  real delta;
  real<lower=0> sigma2;
}
model {
  target += cauchy_lpdf(delta | 0, r);
  target += log(1/sigma2);
  target += normal_lpdf(y | delta*sqrt(sigma2), sqrt(sigma2));
}

…and get mu and sd from delta and sigma2. That said, is there a pystan version of bridgesampling::bridge_sampler and bridgesampling::error_measures?

2

There is no bridgesampling package for python yet.

1 Like
3

Thanks for that, @ahartikainen Ari. That said, would you be aware of some instructive guides as to how one might go about calculating Bayes Factors with pystan? I found this thread, and I agree that it may not be the best tool, but it would still be cool to be able to do that.