I'd like to calculate Bayes Factor for two-sample t-test $H_0: \mu_1=\mu_2$ (model $M_0$) against $H_1: \mu_1\not=\mu_2$ (model $M_1$)

My data are: $x_1,x_2,\ldots, x_{n_1}\sim N(\mu_1, \sigma)$ and $y_1,y_2,\ldots, y_{n_2}\sim N(\mu_2, \sigma)$

BF is a ratio of two marginal likelihoods, so I have to get $$M_H=\int_{\mu\in M_1} f_H(\mu;y)p_H(\mu)d\mu,$$ where $M_i$ is a parameter space under hypothesis $H_0$ or $H_1$, $f_H$ is a probability density function of the data under hypothesis $H$, $p_H$ is a prior distribution on the parameters.

I need to choose prior on $\mu$ (and I do this later) and also let's assume that $\sigma$ is known and equal to both groups (for ease of calculations). The point which I stuck in is how to calculate likelihood for two groups coming from different distributions.

\begin{align} L(\mu_1, \mu_2&| x_1,\ldots, x_{n_1},y_1,\ldots,y_{n_2})\\ &=\prod_{i=1}^{n_1}\frac{1}{\sqrt{2\sigma^2\pi}}\cdot\exp\big\{-\frac{(x_i-\mu_1)^2}{2\sigma^2}\big\}\times\prod_{i=1}^{n_2}\frac{1}{\sqrt{2\sigma^2\pi}}\cdot\exp\big\{-\frac{(y_i-\mu_2)^2}{2\sigma^2}\big\}\\ &= (2\sigma^2\pi)^{\frac{-n_1}{2}}\cdot\exp\big\{{-\frac{1}{2\sigma^2}{\sum_{j=1}^{n_1}(x_j-\mu_1)^2}}\big\}\cdot (2\sigma^2\pi)^{\frac{-n_2}{2}}\cdot\exp\big\{{-\frac{1}{2\sigma^2}{\sum_{j=1}^{n_2}(y_j-\mu_2)^2}}\big\}\\ &= (2\sigma^2\pi)^{-\frac{n_1+n_2}{2}}\cdot\exp\left\{{-\frac{1}{2\sigma^2} ({\sum_{j=1}^{n_1}(x_j-\mu_1)^2}-\frac{1}{2\sigma^2}{\sum_{k=1}^{n_2}(y_k-\mu_2)^2}})\right\} \end{align}

The questions are:

  1. Is it ... Am I horribly wrong or something?
  2. Is it correct likelihood function for two groups form different normal distributions?
  3. How to apply prior to this?

I don't want to reparametrize or place prior on effect size.

  • $\begingroup$ I corrected mistakes in your likelihood, which now stands correct(ed). The derivation of the marginal likelihoods under both hypotheses is processed in Bayesian textbooks when using a conjugate prior on $\mu_1$ and a conjugate prior on $(\mu_1,\mu_2)$. See for instance (!) our book. $\endgroup$
    – Xi'an
    Commented Nov 15, 2016 at 16:56
  • 1
    $\begingroup$ @Xi'an I saw you book (chap.12, example 2.1) but I don't want to use conjugate prior. Because my next step is to use two different distributions (like normal and exponential), where I can't use the mu+eta, mu-eta reparametrization.. $\endgroup$ Commented Nov 15, 2016 at 17:13

1 Answer 1


This is an edited version of Example 2.1 in Bayesian Essentials with R.

When comparing two id normal samples, $$(x_1,\ldots,x_n)\ \text{and }\ (y_1,\ldots, y_n)$$ with respective distributions $\mathscr{N}(\mu_x,\sigma^2)$ and $\mathscr{N}(\mu_y,\sigma^2)$, assume the question is whether or not both means are identical, i.e. $$\mathbf{H_0}:\ \mu_x=\mu_y$$ Further, assume that $\sigma^2$ has a similar meaning under both models and that the same prior $\pi_\sigma(\sigma^2)$ is used under both models. This means that the Bayes factor $$ B^\pi_{21}(\mathscr{D}_n) = \frac{\int\,\ell_2(\mu_x,\mu_y,\sigma|\mathscr{D}_n) \pi(\mu_x,\mu_y)\pi_\sigma(\sigma^2)\,\mathrm{d}\sigma^2\,\mathrm{d}\mu_x\,\mathrm{d}\mu_y} {\int\,\ell_1(\mu,\sigma|\mathscr{D}_n) \pi_\mu(\mu)\pi_\sigma(\sigma^2)\,\mathrm{d}\sigma^2\,\mathrm{d}\mu} $$ does not depend on the normalizing constant of $\pi_\sigma(\sigma^2)$ and thus that we can use an improper prior such as $\pi_\sigma(\sigma^2)=1/\sigma^2$. Furthermore, we can rewrite $\mu_x$ and $\mu_y$ as $$\mu_x=\mu-\xi\ \text{and} \ \mu_y=\mu+\xi$$ respectively, with $\mu$ being a global location parameter, hence a nuisance parameter appearing in both models. Introducing a prior of the form $\pi(\mu,\xi)=\pi_\mu(\mu)\pi_\xi(\xi)$ on the new parameterization, the same prior $\pi_\mu$ can be used underboth models. Once again, an improper Jeffreys prior $\pi_\mu(\mu)=1$ can be used.

However, we do need a proper and well-defined prior on $\xi$, for instance $$\xi\sim\mathscr{N}(0,\tau^2)$$ which leads to \begin{align*} B^\pi_{21}(\mathscr{D}_n) &= \dfrac{\int\, e^{-n\left[(\mu-\xi-\bar x)^2+(\mu+\xi-\bar y)^2+s_{xy}^2\right]/2\sigma^2}\, \sigma^{-2n-2} e^{-\xi^2/2\tau^2}/\tau\sqrt{2\pi}\,\mathrm{d}\sigma^2\,\mathrm{d}\mu\,\mathrm{d}\xi} {\int\,e^{-n \left[(\mu-\bar x)^2+(\mu-\bar y)^2+s_{xy}^2\right]/2\sigma^2}\, \sigma^{-2n-2} \,\mathrm{d}\sigma^2\,\mathrm{d}\mu}\\ %&= \dfrac{2^{2(n+1)/2}}{2^n}\times\\ & = \dfrac{\int\,\left[(\mu-\xi-\bar x)^2+(\mu+\xi-\bar y)^2+s_{xy}^2\right]^{-n} e^{-\xi^2/2\tau^2}/\tau\sqrt{2\pi}\,\mathrm{d}\mu\,\mathrm{d}\xi} {\int\,\left[(\mu-\bar x)^2+(\mu-\bar y)^2+s_{xy}^2\right]^{-n}\,\mathrm{d}\mu}\,, \end{align*} where $s_{xy}^2$ denotes the average $$ s_{xy}^2 = \frac{1}{n}\,\sum_{i=1}^n\,(x_i-\bar x)^2 + \frac{1}{n}\,\sum_{i=1}^n\,(y_i-\bar y)^2 \,. $$

End of the reproduction.

One new item in the OP's question is to avoid reparameterisation, which is easily done by getting back from $(\eta,\xi)$ to $(\mu_x,\mu_y)$ by the change-of-variable lemma. (Hint: the Jacobian is constant.)


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.