This is my first time doing a Bayesian analysis, so I'm not sure whether what I did makes perfect sense. I'm trying to tell if two samples come from the same distribution, more specifically, if they have the same mean.

I have very few samples, and no subjective beliefs about them. I am modelling the observed distribution as a Gaussian. Using pyMC:

import pymc as pm

sample1 = [8.254828927, 8.485694524,10.58058423,5.221356325]
sample2 = [7.921107724, 9.744301571, 5.750874082,8.421883012,11.09813068]

sample_mean = sum(sample1+sample2)

#Empirical prior for the mean
normal_mean_prior = pm.Normal("mean_prior", mu = sample_mean, tau = 1)

#Uninformative prior for the precision
normal_precision_prior = pm.Uniform("precision_prior", lower = 0, upper = 100)

#Normal observation
normal_observation_sample1 = pm.Normal("normal1", mu = normal_mean_prior, tau = normal_precision_prior, observed = True, value = sample1)
normal_observation_sample2 = pm.Normal("normal2", mu = normal_mean_prior, tau = normal_precision_prior, observed = True, value = sample2)

With this model, I would sample it using MCMC and check if the posterior distribution of the means are different or the same.

Is my analysis correct? Have I made any significant blunder? Thank you.

  • 1
    $\begingroup$ This line: sample_mean = sum(sample1+sample2). Did you want to divide by sample size? $\endgroup$ – Cam.Davidson.Pilon Sep 19 '14 at 21:32
  • $\begingroup$ Yes, that is what I wanted. Interestingly, the posterior mean was reasonable, close to a frequentist answer, so I couldn't realize there was a mistake until you saw it. (p.s. I love your book, it's what motivated me to start exploring pyMC. Thank you!) $\endgroup$ – Pedro Tabacof Sep 22 '14 at 13:05

I would say that you have an good start, but your analysis is incomplete. How did you pick your hyper-parameters? How will you check if the posterior distribution of the means are the same? The devil is in the details.

@TomMinka has answered with a reference to a paper by Kruschke that proceeds in the general direction you have been going and fills in the details. Here is the PyMC implementation of something along those lines (but using the normal instead of t for the likelihood to keep it simple):

def model():
    sample_0 = [8.254828927, 8.485694524, 10.58058423, 5.221356325]
    sample_1 = [7.921107724, 9.744301571, 5.750874082, 8.421883012, 11.09813068]
    pooled_data = np.hstack((sample_0, sample_1))

    # hyper-parameters
    S = 1000 * np.std(pooled_data)
    M = np.mean(pooled_data)

    L = .001 * np.std(pooled_data)
    H = 1000 * np.std(pooled_data)

    # prior
    mu = pm.Normal('mu_1', M, S**-2, value=[M,M])
    sigma = pm.Uniform('sigma_1', L, H, value=[np.std(pooled_data), np.std(pooled_data)])

    # likelihood
    obs_1 = pm.Normal('obs_1', mu=mu[0], tau=sigma[0]**-2, value=sample_0, observed=True)
    obs_2 = pm.Normal('obs_2', mu=mu[1], tau=sigma[1]**-2, value=sample_1, observed=True)

    # derived quantities of interest
    d_mu = mu[0] - mu[1]
    d_sigma = sigma[0] - sigma[1]
    def effect_size(d_mu=d_mu, sigma=sigma):
        return d_mu / np.sqrt(.5*(sigma[0]**2 + sigma[1]**2))
    def outside_rope(effect_size=effect_size):
        return -.1 < effect_size < .1

    return locals()

And here is an IPython Notebook to sample from the approximation posterior distribution with MCMC to see that the probability that the two distributions have different means is 0.1.

  • $\begingroup$ I noticed a weird thing about this Kruschke approach, in the absence of data, the posterior probability of outside_rope is 2%. Maybe I made a mistake in my implementation, but if not... is that low? $\endgroup$ – Abraham D Flaxman Sep 19 '14 at 23:40
  • $\begingroup$ In the absence of data, Kruschke's rule for setting the hyperparameters doesn't apply. So you must have done something different. $\endgroup$ – Tom Minka Sep 20 '14 at 14:29
  • $\begingroup$ Good point. I used 0 instead of the mean and 1e-6 instead of the standard deviation. Here is the notebook for it. $\endgroup$ – Abraham D Flaxman Sep 20 '14 at 21:43
  • $\begingroup$ I see a small mistake on the code: "pooled_data = np.hstack((sample_0, sample_0))" is doubling sample_0 and ignoring sample_1. In any case, I don't understand something about your code: Why do you set the values of the prior distributions? Does it influence the posterior sampling? Thank you. $\endgroup$ – Pedro Tabacof Sep 22 '14 at 13:29
  • 1
    $\begingroup$ @PedroTabacof: I see. For non-observed stochastics, the value parameter is used only as the initial point for MCMC. It is not necessary to provide it, and if you do not, PyMC will choose an initial value randomly from the stochastic's prior distribution. Well-chosen initial values can speed convergence. $\endgroup$ – Abraham D Flaxman Sep 22 '14 at 19:22

Your prior on the mean is too narrow. Also, you don't just want samples for the two means, you want samples of the difference between the means. The posterior distribution of the difference is what you want to look at. For more details, see the paper Bayesian estimation supersedes the t test.

  • $\begingroup$ How can I know whether the prior is too narrow or too wide? This was the part I had most difficulty with, so I ended up choosing the precision rather arbitrarily. I'm gonna read Kruschke's paper to understand what I'm doing better, thanks for the link. $\endgroup$ – Pedro Tabacof Sep 22 '14 at 13:18
  • 3
    $\begingroup$ You said that you had no subjective beliefs, which I interpret to mean that you don't want the prior to influence the results. Choosing tau=1 will definitely influence the results, since the likelihood function for the mean has precision around 1 for both datasets. $\endgroup$ – Tom Minka Sep 22 '14 at 18:11
  • $\begingroup$ The method BEST discussed in detail in the paper is also one of the exampes for PyMC. Almost as easy as plug in your data and interpret the results. pymc-devs.github.io/pymc3/BEST $\endgroup$ – Vladislavs Dovgalecs Apr 12 '16 at 21:55

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.