Bayes Factor and likelihood for two sample from different distributions?

Question

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.

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