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$ – Cam.Davidson.Pilon Dec 26 '14 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 '15 at 22:02

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)
| cite | improve this answer | |
  • $\begingroup$ 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 $\endgroup$ – pyring Aug 15 '19 at 9:05
  • $\begingroup$ I don't think I followed the "burn-in" part, anywhere i can read an explanation about that please? $\endgroup$ – baxx Aug 15 at 14:25

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)
| cite | improve this answer | |

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.

| cite | improve this answer | |
  • 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 '15 at 21:03

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.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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