Q58. Joint Normal conjugacy and exact sampling
[20 marks; approximately 25 minutes]
Let X_1,\ldots,X_n\mid\mu,\sigma^2 \overset{\mathrm{iid}}{\sim}N(\mu,\sigma^2), with both \mu\in\mathbb R and \sigma^2>0 unknown, and n\ge1. Write \theta=(\mu,\sigma^2). Prior densities below are with respect to d\mu\,d\sigma^2. Use the inverse-Gamma convention p(v)\propto v^{-a-1}e^{-b/v} for v>0.
(i) Show the Normal belongs to the exponential family in the form f(x\mid\theta) =\frac{b(x)\exp\{c_1(\theta)T_1(x)+c_2(\theta)T_2(x)\}} {a(\theta)}. Identify b,a,c_1,c_2,T_1,T_2 (signs go in the c_j). Hence write the likelihood of the whole sample in terms of \sum_i x_i and \sum_i x_i^2. [4 marks]
(ii) The conjugate prior for this family is \pi(\theta)\propto \frac{\exp\{c_1(\theta)\alpha_1+c_2(\theta)\alpha_2\}} {a(\theta)^\beta}. Write it out explicitly in \mu and \sigma^2. By completing the square in \mu, show it factorises as \mu\mid\sigma^2\sim N\!\left(m_0,\frac{\sigma^2}{\kappa_0}\right), \qquad\sigma^2\sim\operatorname{InvGamma}(a_0,b_0), and find m_0,\kappa_0,a_0,b_0 in terms of \alpha_1,\alpha_2,\beta. State the conditions on the hyperparameters for a proper distribution. [6 marks]
(iii) Show that the posterior is in the same family, with hyperparameters (\alpha_1+\sum_i x_i,\,\alpha_2+\sum_i x_i^2,\,\beta+n). Then give the posterior m_n,\kappa_n,a_n,b_n in terms of \bar x, s^2=\sum_i(x_i-\bar x)^2, and the prior quantities. Here s^2 denotes the unnormalised sum of squares, not the sample variance. [4 marks]
(iv) Show that \mathbb{E}[\mu\mid{\boldsymbol{x}}] is a weighted average of m_0 and \bar x, and give the weights. What happens as n\to\infty with fixed prior hyperparameters? [3 marks]
(v) Describe how to draw exact (non-MCMC) samples from the joint posterior of (\mu,\sigma^2), and explain why it works. Assess the claim: “With a flat prior on \mu and an independent inverse-Gamma prior on \sigma^2, exact sampling is unavailable and a Gibbs sampler is necessary.” [3 marks]
Definitions. Write v=\sigma^2, \bar x=n^{-1}\sum_i x_i, and s^2=\sum_i(x_i-\bar x)^2. Here s^2 is the unnormalised sum of squares, not the sample variance. N(m,v) uses variance v; \operatorname{InvGamma}(a,b) has kernel v^{-a-1}e^{-b/v} on v>0. A Normal–inverse-Gamma factorisation means p(\mu,v)=p(v)p(\mu\mid v), not independence of \mu and v.
Assumptions. Use the specified conditionally iid Normal model with n\ge1. The joint prior density is with respect to d\mu\,dv and has exactly the kernel stated in the question. Changing to precision, or starting with a density in natural-parameter coordinates, would require a Jacobian; no such change is made here.
(i) Exponential-family form and likelihood. Expanding (x-\mu)^2=x^2-2\mu x+\mu^2 gives f(x\mid\mu,v) =\frac1{\sqrt{2\pi v}} \exp\left(-\frac{x^2}{2v}+\frac{\mu x}{v}-\frac{\mu^2}{2v}\right) =\frac{(2\pi)^{-1/2}\exp\{(\mu/v)x-x^2/(2v)\}} {\sqrt v\exp\{\mu^2/(2v)\}}. One valid identification, with signs in the c_j, is \boxed{\begin{aligned} b(x)&=(2\pi)^{-1/2},& a(\theta)&=\sqrt v\exp\{\mu^2/(2v)\},\\ T_1(x)&=x,&c_1(\theta)&=\frac\mu v,\\ T_2(x)&=x^2,&c_2(\theta)&=-\frac1{2v}. \end{aligned}} Let S_1=\sum_i x_i and S_2=\sum_i x_i^2. Then \boxed{L(\mu,v;{\boldsymbol{x}}) =(2\pi v)^{-n/2} \exp\left[-\frac{S_2-2\mu S_1+n\mu^2}{2v}\right].} Thus the sample enters the likelihood through (S_1,S_2).
(ii) Conjugate prior and parameter matching. Substitution into the given prior formula yields the joint density kernel \pi(\mu,v) \propto v^{-\beta/2} \exp\left[-\frac{\beta\mu^2-2\alpha_1\mu+\alpha_2}{2v}\right]. For \beta>0, complete the square: \beta\mu^2-2\alpha_1\mu+\alpha_2 =\beta\left(\mu-\frac{\alpha_1}{\beta}\right)^2 +\alpha_2-\frac{\alpha_1^2}{\beta}. Set \boxed{m_0=\frac{\alpha_1}{\beta},\qquad \kappa_0=\beta,\qquad a_0=\frac{\beta-3}{2},\qquad b_0=\frac12\left(\alpha_2-\frac{\alpha_1^2}{\beta}\right).} To see the shape parameter, factor the kernel as \pi(\mu,v)\propto \underbrace{v^{-1/2} \exp\left[-\frac{\kappa_0(\mu-m_0)^2}{2v}\right]} _{\text{Normal conditional kernel, including its }v^{-1/2}\text{ factor}} \underbrace{v^{-(\beta-1)/2}e^{-b_0/v}} _{\text{inverse-Gamma marginal kernel}}. The Normal factor has a normalising multiplier \sqrt{\kappa_0/(2\pi)} that does not depend on v. The remaining exponent therefore satisfies a_0+1=(\beta-1)/2, giving a_0=(\beta-3)/2. Hence \mu\mid v\sim N(m_0,v/\kappa_0) and v\sim\operatorname{InvGamma}(a_0,b_0).
The necessary and sufficient conditions for a proper prior are \boxed{\alpha_1\in\mathbb R,\qquad\beta>3,\qquad \alpha_2>\frac{\alpha_1^2}{\beta}.} They ensure \kappa_0,a_0,b_0>0. The specified three-hyperparameter prior is a subfamily of the general Normal–inverse-Gamma family: it imposes \kappa_0=2a_0+3.
(iii) Posterior update. Multiplying the prior kernel by the likelihood gives \pi(\mu,v\mid{\boldsymbol{x}}) \propto v^{-(\beta+n)/2} \exp\left[-\frac{(\beta+n)\mu^2 -2(\alpha_1+S_1)\mu+(\alpha_2+S_2)}{2v}\right]. This has exactly the same form, with \boxed{\alpha_1'=\alpha_1+S_1,\qquad \alpha_2'=\alpha_2+S_2,\qquad\beta'=\beta+n.} Using S_1=n\bar x, S_2=s^2+n\bar x^2, and the prior identities \alpha_1=\kappa_0m_0 and \alpha_2=2b_0+\kappa_0m_0^2, we obtain \boxed{\begin{aligned} \kappa_n&=\kappa_0+n,\\ m_n&=\frac{\kappa_0m_0+n\bar x}{\kappa_0+n},\\ a_n&=a_0+\frac n2,\\ b_n&=b_0+\frac12\left[s^2+ \frac{\kappa_0n}{\kappa_0+n}(\bar x-m_0)^2\right]. \end{aligned}} In particular, the final update follows from \begin{aligned} b_n &=\frac12\left[\alpha_2+S_2- \frac{(\alpha_1+S_1)^2}{\beta+n}\right]\\ &=b_0+\frac12\left[s^2+\kappa_0m_0^2+n\bar x^2 -\frac{(\kappa_0m_0+n\bar x)^2}{\kappa_0+n}\right]\\ &=b_0+\frac12\left[s^2+ \frac{\kappa_0n}{\kappa_0+n}(\bar x-m_0)^2\right]. \end{aligned} Therefore \boxed{\mu\mid v,{\boldsymbol{x}}\sim N(m_n,v/\kappa_n),\qquad v\mid{\boldsymbol{x}}\sim\operatorname{InvGamma}(a_n,b_n).} The v distribution here is the marginal posterior, with \mu integrated out. All updates remain proper when the prior is proper.
(iv) Posterior mean and weights. For n\ge1 and a proper prior, a_n=a_0+n/2>1/2, which ensures that the marginal posterior mean of \mu exists. Iterated expectation gives \mathbb{E}[\mu\mid{\boldsymbol{x}}] =\mathbb{E}\!\left[\mathbb{E}[\mu\mid v,{\boldsymbol{x}}]\mid{\boldsymbol{x}}\right] =m_n. Consequently, \boxed{\mathbb{E}[\mu\mid{\boldsymbol{x}}] =\underbrace{\frac{\kappa_0}{\kappa_0+n}}_{w_0}m_0 +\underbrace{\frac n{\kappa_0+n}}_{w_x}\bar x.} The weights are positive and sum to one. For fixed prior hyperparameters, w_0\to0 and w_x\to1: the prior contribution vanishes and the posterior mean approaches the sample mean. Under iid sampling at fixed true parameters, \bar x\to\mu_{\mathrm{true}} almost surely, so m_n\to\mu_{\mathrm{true}} almost surely as well.
(v) Exact joint-posterior sampling. For each j=1,\ldots,M, independently of all previous pairs:
1. Draw v^{(j)}\sim\operatorname{InvGamma}(a_n,b_n) from the marginal posterior of v=\sigma^2.
2. Draw \mu^{(j)}\sim N(m_n,v^{(j)}/\kappa_n), conditional on that newly drawn variance.
3. Store (\mu^{(j)},v^{(j)}).
For example, step 1 can use G^{(j)}\sim\operatorname{Gamma}(a_n,b_n) in shape–rate form and set v^{(j)}=1/G^{(j)}. With fresh random numbers in each repetition, this produces iid pairs from the joint posterior because p(\mu,v\mid{\boldsymbol{x}})=p(v\mid{\boldsymbol{x}})\,p(\mu\mid v,{\boldsymbol{x}}). The two coordinates within a pair are generally dependent, but distinct pairs are independent. There is no Markov-chain convergence, burn-in, or accept/reject step to manage in this construction.
Assessing the flat-prior claim. The claim is false for a flat mean prior and an independent inverse-Gamma variance prior. To distinguish its hyperparameters from those above, suppose p(\mu)\propto1,\qquad v\sim\operatorname{InvGamma}(a_*,b_*),\qquad a_*,b_*>0. Although this joint prior is improper in \mu, for n\ge1 the posterior is proper. Using \sum_i(x_i-\mu)^2=s^2+n(\mu-\bar x)^2, its kernel is p(\mu,v\mid{\boldsymbol{x}}) \propto v^{-a_*-1-n/2} \exp\left[-\frac{b_*+s^2/2}{v} -\frac{n(\mu-\bar x)^2}{2v}\right]. Integrating over \mu contributes a factor proportional to \sqrt v, since \int_{-\infty}^{\infty} \exp\left[-\frac{n(\mu-\bar x)^2}{2v}\right]d\mu =\sqrt{\frac{2\pi v}{n}}. Hence \boxed{v\mid{\boldsymbol{x}}\sim \operatorname{InvGamma}\!\left(a_*+\frac{n-1}{2},\ b_*+\frac{s^2}{2}\right), \qquad\mu\mid v,{\boldsymbol{x}}\sim N(\bar x,v/n).} One can again draw from the variance marginal and then the Normal conditional, producing independent exact posterior pairs. Gibbs sampling is valid here, but not necessary. The same distinction appears in Question 43 for p(\mu,v)\propto1/v, with its stated propriety conditions.
Key result. Exact direct sampling uses a marginal then a conditional. Gibbs sampling alternates full conditionals. The key issue is whether the needed marginal is available for direct sampling, not merely whether the prior was described as independent.