Past exam of the mathematics course of the University of Cambridge / 2002 / iii / Paper 40 / 3 / a

For a target probability density function $\pi ({\theta }_1,\dots ,{\theta }_p)$, one systematic sweep of the Gibbs sampler updates coordinates successively:Each update leaves the target invariant: keeping the other coordinates with their target marginal and replacing one coordinate by its conditional distribution reconstructs the target joint probability density function. Therefore composition of the coordinate updates also preserves $\pi$. It defines a Markov chain whose consecutive draws are generally dependent. Under irreducibility, aperiodicity and the usual recurrence conditions it converges to the target law from an arbitrary valid start; a stationary start gives exact target marginals immediately. Invariance alone does not prove these ergodicity conditions, and a systematic full sweep need not itself be reversible.

Write $v={\sigma }^2$, $d=k+1$, $X=X_k$, $a=a_k$, $\mu ={\mu }_k$ and $\Sigma ={\Sigma }_k$, with $\Sigma$ positive-definite. For fixed $v$, the terms involving $a$ in the posterior exponent areComplete the square. With $V_k=(X_k^T{X}_k/v+{\Sigma }_k^{-1}{)}^{-1}$ and $m_k=V_k(X_k^{T}y/v+{\Sigma }_k^{-1}{\mu }_k)$, this givesFor fixed $a_k$, collect the powers and exponential terms involving $v$:Hence, in the paper's shape-scale inverse-gamma distribution convention,There is no extra $(k+1)/2$ in this shape, because the normal coefficient prior is independent of $v$. These are the independent normal and inverse-gamma regression priors conditionals.

The Gaussian normalization printed in the PDF is incorrect. The proper $d$-dimensional prior probability density function is $(2\pi {)}^{-d/2}|{\Sigma }_k{|}^{-1/2}\exp [-(a_k-{\mu }_k{)}^T{\Sigma }_k^{-1}(a_k-{\mu }_k)/2]$. For fixed $k$ this correction is constant in $a_k,v$ and leaves the two conditionals unchanged, but it must be retained when comparing different model orders in part (b).

Start with $v^{(0)}>0$ and any coefficient vector. At sweep $m$, draw $a_k^{(m+1)}$ from the multivariate normal distribution above using $v^{(m)}$, then draw $v^{(m+1)}$ from the inverse-gamma distribution using the newly drawn coefficients. For the second step one may generate $1/v$ from a gamma distribution with the same shape and rate equal to the displayed inverse-gamma scale.

Each block update uses an exact conditional distribution, so it preserves the joint posterior distribution. The resulting pairs form a dependent Gibbs sampler sample. With proper priors, positive scale and positive-definite coefficient covariance, these conditionals have full probability support on ${\mathbb{R}}^{k+1}\times (0,\infty )$ and the usual posterior ergodicity applies. A burn-in reduces dependence on initialization; retained sweeps are still serially dependent and Monte Carlo errors should account for that dependence.