Question 18. Bayesian model selection
Suppose Y_1,\ldots,Y_n are independent observations, with observed data {\boldsymbol{y}}=(y_1,\ldots,y_n). A researcher wishes to compare two competing models.
Under model M_1,
Y_i\mid\theta,M_1\stackrel{\mathrm{iid}}{\sim}\operatorname{N}(\theta,1),
where the unknown mean has prior distribution
\theta\mid M_1\sim \operatorname{N}(0,\tau^2).
Under model M_2,
Y_i\mid\mu,\sigma^2,M_2\stackrel{\mathrm{iid}}{\sim}\operatorname{N}(\mu,\sigma^2),
where the unknown mean and variance have prior distributions
\mu\mid\sigma^2,M_2\sim \operatorname{N}\Bigl(0,\frac{\sigma^2}{\kappa}\Bigr), \qquad \sigma^2\mid M_2\sim\operatorname{InvGamma}(\alpha,\beta).
Assume that n\geq 2, and \tau^2,\kappa,\alpha,\beta are known positive constants.
The inverse-gamma distribution is parameterised by
f(v;a,b)=\frac{b^{a}}{\Gamma(a)}v^{-a-1}\exp\Bigl(-\frac{b}{v}\Bigr), \qquad v>0.
Define
\bar y=\frac1n\sum_{i=1}^{n}y_i, \qquad S=\sum_{i=1}^{n}(y_i-\bar y)^{2}.
You may use the identity
\sum_{i=1}^{n}(y_i-\eta)^{2}=S+n(\bar y-\eta)^{2}
and the result
\int_{0}^{\infty}v^{-a-1}e^{-b/v}\,dv=\frac{\Gamma(a)}{b^{a}}, \qquad a,b>0.
(a) Write down the likelihood functions L_1(\theta;{\boldsymbol{y}}) and L_2(\mu,\sigma^2;{\boldsymbol{y}}) under models M_1 and M_2, respectively. [2 marks]
(b) Using the prior distribution of \theta, derive an expression for the marginal likelihood p({\boldsymbol{y}}\mid M_1). The derivation should explicitly show the integration over \theta, including any completion of the square required. [5 marks]
(c) The marginal likelihood under model M_2 is given by
p({\boldsymbol{y}}\mid M_2) = \int_0^{\infty}\int_{-\infty}^{\infty} p({\boldsymbol{y}}\mid\mu,\sigma^2,M_2)\, p(\mu\mid\sigma^2,M_2)\, p(\sigma^2\mid M_2)\,d\mu\,d\sigma^2.
(i) Integrate out \mu to obtain an expression for p({\boldsymbol{y}}\mid\sigma^2,M_2). [4 marks]
(ii) Hence integrate out \sigma^2 to derive the marginal likelihood p({\boldsymbol{y}}\mid M_2). [4 marks]
(d)
(i) Obtain an expression for the Bayes factor
B_{12}=\frac{p({\boldsymbol{y}}\mid M_1)}{p({\boldsymbol{y}}\mid M_2)}.
(ii) Suppose the prior model probabilities are P(M_1)=0.25 and P(M_2)=0.75. Express the posterior odds in favour of M_1 in terms of B_{12}.
[3 marks]
(e)
(i) Explain why normalising constants that do not depend on the unknown parameters must nevertheless be retained when calculating marginal likelihoods.
(ii) State the conditions under which the Bayes factor is well-defined. Explain briefly why replacing the proper parameter priors with unspecified improper priors can cause difficulties in Bayesian model comparison.
[2 marks]
Definitions. \operatorname{N}(\eta,v) has mean \eta and variance v. Write \bar y=n^{-1}\sum_i y_i and S=\sum_i(y_i-\bar y)^2. Model M_1 has known variance one and unknown mean \theta; M_2 has unknown mean and variance with a dependent normal–inverse-gamma prior. The quantity B_{12}=p({\boldsymbol{y}}\mid M_1)/p({\boldsymbol{y}}\mid M_2) compares prior predictive densities of the same data, and each parameter prior is distinct from the prior probability of its model.
(a) By the given sum-of-squares identity, for \theta\in\mathbb{R},
\boxed{ L_1(\theta;{\boldsymbol{y}}) = (2\pi)^{-n/2} \exp\Bigl(-\tfrac12\bigl[S+n(\bar y-\theta)^{2}\bigr]\Bigr). }
For \mu\in\mathbb{R} and \sigma^{2}>0,
\boxed{ L_2(\mu,\sigma^{2};{\boldsymbol{y}}) = (2\pi\sigma^{2})^{-n/2} \exp\Bigl(-\frac{S+n(\bar y-\mu)^{2}}{2\sigma^{2}}\Bigr). }
(b) The \operatorname{N}(0,\tau^{2}) prior has density (2\pi\tau^{2})^{-1/2}\exp(-\theta^{2}/(2\tau^{2})). Thus
\begin{aligned} p({\boldsymbol{y}}\mid M_1) &= \int_{-\infty}^{\infty}L_1(\theta;{\boldsymbol{y}})\,(2\pi\tau^{2})^{-1/2}e^{-\theta^{2}/(2\tau^{2})}\,d\theta\\ &= (2\pi)^{-n/2}(2\pi\tau^{2})^{-1/2}e^{-S/2} \int_{-\infty}^{\infty} \exp\Bigl(-\tfrac12\Bigl[n(\bar y-\theta)^{2}+\frac{\theta^{2}}{\tau^{2}}\Bigr]\Bigr)\,d\theta. \end{aligned}
The expression in brackets is
n(\bar y-\theta)^{2}+\frac{\theta^{2}}{\tau^{2}} = \Bigl(n+\frac1{\tau^{2}}\Bigr)\theta^{2}-2n\bar y\,\theta+n\bar y^{2}.
Set A=n+1/\tau^{2}. Then A>0, because n\geq 2 and \tau^{2}>0, and \tau^{2}A=1+n\tau^{2}. Completing the square,
\begin{aligned} A\theta^{2}-2n\bar y\,\theta+n\bar y^{2} &= A\Bigl(\theta-\frac{n\bar y}{A}\Bigr)^{2} +n\bar y^{2}\frac{A-n}{A}\\ &= A\Bigl(\theta-\frac{n\bar y}{A}\Bigr)^{2} +\frac{n\bar y^{2}}{1+n\tau^{2}}, \end{aligned}
where A-n=1/\tau^{2}. The integral equals
\exp\Bigl(-\frac{n\bar y^{2}}{2(1+n\tau^{2})}\Bigr) \int_{-\infty}^{\infty} \exp\Bigl(-\frac A2\Bigl(\theta-\frac{n\bar y}{A}\Bigr)^{2}\Bigr)\,d\theta = \exp\Bigl(-\frac{n\bar y^{2}}{2(1+n\tau^{2})}\Bigr) \sqrt{\frac{2\pi}{A}},
since A>0. The remaining normalising factors simplify by \tau^{2}A=1+n\tau^{2}:
(2\pi\tau^{2})^{-1/2}\sqrt{2\pi/A}=(1+n\tau^{2})^{-1/2}.
\boxed{ p({\boldsymbol{y}}\mid M_1) = (2\pi)^{-n/2}(1+n\tau^{2})^{-1/2} \exp\Bigl(-\frac12\Bigl[S+\frac{n\bar y^{2}}{1+n\tau^{2}}\Bigr]\Bigr). }
(c)
(i) The conditional prior of \mu has density (\kappa/(2\pi\sigma^{2}))^{1/2}\exp(-\kappa\mu^{2}/(2\sigma^{2})). Therefore
\begin{aligned} p({\boldsymbol{y}}\mid\sigma^{2},M_2) &= (2\pi\sigma^{2})^{-n/2}\Bigl(\frac{\kappa}{2\pi\sigma^{2}}\Bigr)^{1/2}e^{-S/(2\sigma^{2})}\\ &\qquad\times \int_{-\infty}^{\infty} \exp\Bigl(-\frac{n(\bar y-\mu)^{2}+\kappa\mu^{2}}{2\sigma^{2}}\Bigr)\,d\mu. \end{aligned}
Since n+\kappa>0,
\begin{aligned} n(\bar y-\mu)^{2}+\kappa\mu^{2} &= (n+\kappa)\Bigl(\mu-\frac{n\bar y}{n+\kappa}\Bigr)^{2} +n\bar y^{2}\frac{\kappa}{n+\kappa}. \end{aligned}
The integral is
\exp\Bigl(-\frac{n\kappa\bar y^{2}}{2\sigma^{2}(n+\kappa)}\Bigr) \sqrt{\frac{2\pi\sigma^{2}}{n+\kappa}}.
Multiplying by the constant prefactor, the factors (2\pi\sigma^{2})^{\pm 1/2} cancel, and
\boxed{ p({\boldsymbol{y}}\mid\sigma^{2},M_2) = (2\pi\sigma^{2})^{-n/2}\Bigl(\frac{\kappa}{n+\kappa}\Bigr)^{1/2} \exp\Bigl(-\frac{Q}{2\sigma^{2}}\Bigr), }
where Q=S+n\kappa\bar y^{2}/(n+\kappa).
(ii) By (i) and the inverse-gamma prior,
\begin{aligned} p({\boldsymbol{y}}\mid M_2) &= (2\pi)^{-n/2}\Bigl(\frac{\kappa}{n+\kappa}\Bigr)^{1/2} \frac{\beta^{\alpha}}{\Gamma(\alpha)}\\ &\qquad\times \int_0^{\infty} (\sigma^{2})^{-(\alpha+n/2)-1} \exp\Bigl(-\frac{\beta+Q/2}{\sigma^{2}}\Bigr)\,d\sigma^{2}. \end{aligned}
Here \alpha+n/2>0 because \alpha>0 and n\geq 2, and \beta+Q/2>0 because \beta>0 and Q\geq 0. The given integral with these values yields
\boxed{ p({\boldsymbol{y}}\mid M_2) = (2\pi)^{-n/2}\Bigl(\frac{\kappa}{n+\kappa}\Bigr)^{1/2} \frac{\beta^{\alpha}\,\Gamma(\alpha+n/2)}{\Gamma(\alpha)\,(\beta+Q/2)^{\alpha+n/2}}, }
with Q as in (i).
The two integrations under M_2 have different roles. Integrating over \mu combines the sampling information about the mean with its conditional normal prior, producing the adjustment n\kappa\bar y^2/(n+\kappa). Integrating over \sigma^2 then averages over uncertainty in the common variance. Neither step replaces an unknown by its maximum likelihood estimate. The factor \sqrt{\kappa/(n+\kappa)} is one of the normalisation factors that would be lost by keeping only an exponential kernel.
(d)
(i) By (b) and (c)(ii), the factors (2\pi)^{-n/2} cancel, so
\boxed{ B_{12} = (1+n\tau^{2})^{-1/2} \exp\Bigl(-\frac12\Bigl[S+\frac{n\bar y^{2}}{1+n\tau^{2}}\Bigr]\Bigr) \Bigl(\frac{n+\kappa}{\kappa}\Bigr)^{1/2} \frac{\Gamma(\alpha)\,(\beta+Q/2)^{\alpha+n/2}}{\beta^{\alpha}\,\Gamma(\alpha+n/2)}, }
where Q=S+n\kappa\bar y^{2}/(n+\kappa).
(ii)
\frac{\mathbb{P}(M_1\mid{\boldsymbol{y}})}{\mathbb{P}(M_2\mid{\boldsymbol{y}})} = B_{12}\,\frac{\mathbb{P}(M_1)}{\mathbb{P}(M_2)} = B_{12}\,\frac{0.25}{0.75} = \frac{B_{12}}{3}.
\therefore \text{the posterior odds in favour of }M_1\text{ are }B_{12}/3.
Therefore the posterior favours M_1 only if B_{12}>3, rather than merely B_{12}>1. If the two models exhaust the alternatives, \mathbb{P}(M_1\mid{\boldsymbol{y}})=B_{12}/(B_{12}+3). A Bayes factor above one says that the data shift the odds towards M_1; it need not make M_1 the more probable model after accounting for the prior odds.
(e)
(i) A constant free of the parameters cancels in the posterior distribution within one model. It does not cancel between M_1 and M_2, because the two marginal likelihoods need not contain the same constant. Those factors remain in B_{12}.
(ii) B_{12} is well-defined when 0<p({\boldsymbol{y}}\mid M_k)<\infty for k=1,2. That holds here: both parameter priors are proper, n\geq 2, and \tau^{2},\kappa,\alpha,\beta>0, so the integrals in (b) and (c) converge. An unspecified improper prior is determined only up to an arbitrary positive multiple. That multiple remains in the marginal likelihood, so the Bayes factor is not uniquely defined.
A marginal likelihood is a prior-weighted average of likelihood values, so it depends on both fit and how much prior mass lies in parameter regions that fit the data. Broadening a proper prior can change this average even when the posterior within a model changes little. Sensitivity to the stated prior scales is therefore part of interpreting the model comparison; a proper posterior by itself does not justify a Bayes factor based on an unspecified improper prior.