# Bayesian equivalent of two sample t-test?

I'm not looking for a plug and play method like BEST in R but rather a mathematical explanation of what are some Bayesian methods I can use to test the difference between the mean of two samples.

• the original BEST paper might be what you are looking for: indiana.edu/~kruschke/BEST/BEST.pdf Dec 26, 2014 at 19:16
• Just to be clear, are we speaking of a two-sample test that is the equivalent of a frequentist test of mean differences in two groups, such as the t-test? or are you interested in tests of the strong null hypothesis for distributional differences like the Kolmogorov-Smirnoff test? Jan 30, 2015 at 22:02
• @Cam.Davidson.Pilon The link is deprecated. Oct 24, 2022 at 3:30
• Updated link: jkkweb.sitehost.iu.edu/articles/Kruschke2013JEPG.pdf Oct 24, 2022 at 15:54

This is a good question, that seems to pop up a lot: link 1, link 2. The paper Bayesian Estimation Superseeds the T-Test that Cam.Davidson.Pilon pointed out is an excellent resource on this subject. It is also very recent, published in 2012, which I think in part is due to the current interest in the area.

I will try to summarize a mathematical explanation of a Bayesian alternative to the two sample t-test. This summary is similar to the BEST paper which assess the difference in two samples by comparing the difference in their posterior distributions (explained below in R).

set.seed(7)

#create samples
sample.1 <- rnorm(8, 100, 3)
sample.2 <- rnorm(10, 103, 7)

#we need a pooled data set for estimating parameters in the prior.
pooled <- c(sample.1, sample.2)
par(mfrow=c(1, 2))

hist(sample.1)
hist(sample.2)


In order to compare the sample means we need to estimate what they are. The Bayesian method to do so uses Bayes' theorem: P(A|B) = P(B|A) * P(A)/P(B) (the syntax of P(A|B) is read as the probability of A given B)

Thanks to modern numerical methods we can ignore the probability of B, P(B), and use the proportional statment: P(A|B) $$\propto$$ P(B|A)*P(A) In Bayesian vernacular the posterior is proportional to the likelihood times the prior

Applying Bayes' theory to our problem where we want to know the means of samples given some data we get $$P(mean.1 | sample.1)$$ $$\propto$$ $$P(sample.1 | mean.1) * P(mean.1)$$. The first term on the right is the likelihood, $$P(sample.1 | mean.1)$$, which is the probability of observing the sample data given mean.1. The second term is the prior, $$P(mean.1)$$, which is simply the probability of mean.1. Figuring out appropriate priors is still a bit of an art and is one of the biggest critisims of Bayesian methods.

Let's put it in code. Code makes everything better.

likelihood <- function(parameters){
mu1=parameters[1]; sig1=parameters[2]; mu2=parameters[3]; sig2=parameters[4]
prod(dnorm(sample.1, mu1, sig1)) * prod(dnorm(sample.2, mu2, sig2))
}

prior <- function(parameters){
mu1=parameters[1]; sig1=parameters[2]; mu2=parameters[3]; sig2=parameters[4]
dnorm(mu1, mean(pooled), 1000*sd(pooled)) * dnorm(mu2, mean(pooled), 1000*sd(pooled)) * dexp(sig1, rate=0.1) * dexp(sig2, 0.1)
}


I made some assumptions in the prior that need to be justified. To keep the priors from prejudicing the estimated mean I wanted to make them broad and uniform-ish over plausible values with the aim of letting the data produce the features of the posterior. I used recommended setting from BEST and distributed the mu's normally with mean = mean(pooled) and a broad standard deviation = 1000*sd(pooled). The standard deviations I set to a broad exponential distribution, because I wanted a broad unbounded distribution.

Now we can make the posterior

posterior <- function(parameters) {likelihood(parameters) * prior(parameters)}


We will sample the posterior distribution using a markov chain monte carlo (MCMC) with Metropolis Hastings modification. Its easiest to understand with code.

#starting values
mu1 = 100; sig1 = 10; mu2 = 100; sig2 = 10
parameters <- c(mu1, sig1, mu2, sig2)

#this is the MCMC /w Metropolis method
n.iter <- 10000
results <- matrix(0, nrow=n.iter, ncol=4)
results[1, ] <- parameters
for (iteration in 2:n.iter){
candidate <- parameters + rnorm(4, sd=0.5)
ratio <- posterior(candidate)/posterior(parameters)
if (runif(1) < ratio) parameters <- candidate #Metropolis modification
results[iteration, ] <- parameters
}


The results matrix is a list of samples from the posterior distribution for each parameter which we can use to answer our original question: Is sample.1 different than sample.2? But first to avoid affects from the starting values we will "burn-in" the first 500 values of the chain.

#burn-in
results <- results[500:n.iter,]


Now, is sample.1 different than sample.2?

mu1 <- results[,1]
mu2 <- results[,3]

hist(mu1 - mu2)


mean(mu1 - mu2 < 0)
[1] 0.9953689


From this analysis I would conclude there is a 99.5% chance that the mean for sample.1 is less than the mean for sample.2.

An advantage of the Bayesian approach, as pointed out in the BEST paper, is that it can make strong theories. E.G. what is the probability that sample.2 is 5 units bigger than sample.1.

mean(mu2 - mu1 > 5)
[1] 0.9321124


We would conclude that there is a 93% chance that the mean of sample.2 is 5 unit greater than sample.1. An observant reader would find that interesting because we know the true populations have means of 100 and 103 respectively. This is most likely due to the small sample size, and choice of using a normal distribution for the likelihood.

I will end this answer with a warning: This code is for teaching purposes. For a real analysis use RJAGS and depending on your sample size fit a t-distribution for the likelihood. If there is interest I will post a t-test using RJAGS.

EDIT: As requested here is a JAGS model.

model.str <- 'model {
for (i in 1:Ntotal) {
y[i] ~ dt(mu[x[i]], tau[x[i]], nu)
}
for (j in 1:2) {
mu[j] ~ dnorm(mu_pooled, tau_pooled)
tau[j] <- 1 / pow(sigma[j], 2)
sigma[j] ~ dunif(sigma_low, sigma_high)
}
nu <- nu_minus_one + 1
nu_minus_one ~ dexp(1 / 29)
}'

# Indicator variable
x <- c(rep(1, length(sample.1)), rep(2, length(sample.2)))

cpd.model <- jags.model(textConnection(model.str),
data=list(y=pooled,
x=x,
mu_pooled=mean(pooled),
tau_pooled=1/(1000 * sd(pooled))^2,
sigma_low=sd(pooled) / 1000,
sigma_high=sd(pooled) * 1000,
Ntotal=length(pooled)))
update(cpd.model, 1000)
chain <- coda.samples(model = cpd.model, n.iter = 100000,
variable.names = c('mu', 'sigma'))
rchain <- as.matrix(chain)
hist(rchain[, 'mu[1]'] - rchain[, 'mu[2]'])
mean(rchain[, 'mu[1]'] - rchain[, 'mu[2]'] < 0)
mean(rchain[, 'mu[2]'] - rchain[, 'mu[1]'] > 5)

• Just wondering if there is a reasonable solution to use the Bayesian two sample comparison with this type of data sets. stackoverflow.com/q/57503523/7288088 Aug 15, 2019 at 9:05
• I don't think I followed the "burn-in" part, anywhere i can read an explanation about that please?
– baxx
Aug 15, 2020 at 14:25
• @baxx burn-in is removing samples from the beginning of an MCMC chain that are influenced by the starting points and thus don't represent the parameters you want to estimate see e.g. users.stat.umn.edu/~geyer/mcmc/burn.html#what. Also called "warm up" Jan 20, 2021 at 21:05
• An implementation of the same test using PyMC3 with the NUTS sampler: docs.pymc.io/en/v3/pymc-examples/examples/case_studies/… Jul 6, 2022 at 11:12

The excellent answer by user1068430 implemented in Python

import numpy as np
from pylab import plt

def dnorm(x, mu, sig):
return 1/(sig * np.sqrt(2 * np.pi)) * np.exp(-(x - mu)**2 / (2 * sig**2))

def dexp(x, l):
return l * np.exp(- l*x)

def like(parameters):
[mu1, sig1, mu2, sig2] = parameters
return dnorm(sample1, mu1, sig1).prod()*dnorm(sample2, mu2, sig2).prod()

def prior(parameters):
[mu1, sig1, mu2, sig2] = parameters
return dnorm(mu1, pooled.mean(), 1000*pooled.std()) * dnorm(mu2, pooled.mean(), 1000*pooled.std()) * dexp(sig1, 0.1) * dexp(sig2, 0.1)

def posterior(parameters):
[mu1, sig1, mu2, sig2] = parameters
return like([mu1, sig1, mu2, sig2])*prior([mu1, sig1, mu2, sig2])

#create samples
sample1 = np.random.normal(100, 3, 8)
sample2 = np.random.normal(100, 7, 10)

pooled= np.append(sample1, sample2)

plt.figure(0)
plt.hist(sample1)
plt.hold(True)
plt.hist(sample2)
plt.show(block=False)

mu1 = 100
sig1 = 10
mu2 = 100
sig2 = 10
parameters = np.array([mu1, sig1, mu2, sig2])

niter = 10000

results = np.zeros([niter, 4])
results[1,:] = parameters

for iteration in np.arange(2,niter):
candidate = parameters + np.random.normal(0,0.5,4)
ratio = posterior(candidate)/posterior(parameters)
if np.random.uniform() < ratio:
parameters = candidate
results[iteration,:] = parameters

#burn-in
results = results[499:niter-1,:]

mu1 = results[:,1]
mu2 = results[:,3]

d = (mu1 - mu2)
p_value = np.mean(d > 0)

plt.figure(1)
plt.hist(d,normed = 1)
plt.show()
$$$$
`

With a Bayesian analysis you have more things to specify (that is actually a good thing, since it gives much more flexibility and ability to model what you believe the truth to be). Are you assuming normals for the likelihoods? Will the 2 groups have the same variance?

One straight forward approach is to model the 2 means (and 1 or 2 variances/dispersions) then look at the posterior on the difference of the 2 means and/or the Credible Interval on the difference of the 2 means.

• Could you provide some more details on this? I'm not sure how to model 2 means and look at posteriors.
– John
Jan 6, 2015 at 21:03
• I think what Greg meant is that you look at the the posterior on the difference and see how far is it from 0 probabilistically. Apr 7, 2021 at 7:37

a mathematical explanation of what are some Bayesian methods I can use to test the difference between the mean of two samples.

There are several approaches to "testing" this. I'll mention a couple:

• If you want an explicit decision you could look at decision theory.

• A pretty simple thing that's sometimes done is to find an interval for the difference in the means and consider whether it includes 0 or not. That would involve starting with a model for the observations, priors on the parameters and computation of the posterior distribution of the difference in means conditional on the data.

You'd need to say what your model is (e.g. normal, constant variance), and then (at least) some prior for the difference in means and a prior for the variance. You might have priors on the parameters of those priors in turn. Or you might not assume constant variance. Or you might assume something other than normality.

First you should consider whether you would like to do this at all.

A good start would be to read this even more recent (2022) paper by Costello and Watts to learn about the various theoretical pitfalls associated with hypothesis testing in a Bayesian framework.

One of the most serious problems they see with the usual t-test null hypothesis of $$\mathcal{H}_0: \Delta\mu = 0$$ is that this is a "point-form null hypothesis". As $$\Delta\mu$$ is a real number, and its distribution is continuous, the probability $$\mathrm{Pr}(\Delta\mu = 0)$$ will unfortunately be exactly 0. Now, "if the null hypothesis is always false, what’s the big deal about rejecting it?" -- they ask in the paper, quoting Cohen (2016).

After showing that the currently used "Bayesian t-test" is equivalent to its "classical"/frequentist counterpart (thus offering no obvious advantages), Costello and Watts also propose a different testing procedure which is supposed to be more informative. See the paper for details.

Disclaimer: I haven't done such an analysis yet (just reading myself into the subject). I would be rather curious to know if there's someone who has some practical experience with what Costello and Watts propose.