Bayes Factor and likelihood for two sample from different distributions? 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:


*

*Is it ... Am I horribly wrong or something? 

*Is it correct likelihood function for two groups form different
normal distributions?

*How to apply prior to this?


I don't want to reparametrize or place prior on effect size. 
 A: 
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.)
