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.

  • 16
    $\begingroup$ the original BEST paper might be what you are looking for: indiana.edu/~kruschke/BEST/BEST.pdf $\endgroup$ Dec 26, 2014 at 19:16
  • 4
    $\begingroup$ 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? $\endgroup$
    – AdamO
    Jan 30, 2015 at 22:02
  • 1
    $\begingroup$ @Cam.Davidson.Pilon The link is deprecated. $\endgroup$
    – Galen
    Oct 24, 2022 at 3:30
  • 1
    $\begingroup$ Updated link: jkkweb.sitehost.iu.edu/articles/Kruschke2013JEPG.pdf $\endgroup$ Oct 24, 2022 at 15:54

5 Answers 5


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


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


enter image description here

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.

results <- results[500:n.iter,]

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

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

hist(mu1 - mu2)

enter image description here

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),
                                  tau_pooled=1/(1000 * sd(pooled))^2,
                                  sigma_low=sd(pooled) / 1000,
                                  sigma_high=sd(pooled) * 1000,
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)

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)


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

results = results[499:niter-1,:]

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

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

plt.hist(d,normed = 1)

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.

  • 1
    $\begingroup$ Could you provide some more details on this? I'm not sure how to model 2 means and look at posteriors. $\endgroup$
    – John
    Jan 6, 2015 at 21:03
  • $\begingroup$ 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. $\endgroup$
    – Steve
    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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.