# Bayesian t-test using pystan

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 `mu`s and `sd`s 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`?

There is no `bridgesampling` package for python yet.

1 Like

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.