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

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.

The order prior and probabilities of choosing upward or downward moves are not specified, so include them explicitly. Let ${\rho }_k$ be the model-order prior probability, $b_k$ the birth-selection probability at order $k$, and $d_{k+1}$ the reverse death-selection probability. Denote the correctly normalized coefficient prior by $p_k(a)$ and the probability density function of $z\sim N(0,{\sigma }_a^2)$ by $q(z)$. Common variance-prior factors cancel since $v={\sigma }^2$ is unchanged.

Put $m=n^{-1}{\sum }_i{x}_i^{k+1}$, $r=y-X_{k}a$, and $w_i=x_i^{k+1}-m$. The centered polynomial birth move in reversible-jump sampling gives $X_{k+1}{a}^{\prime }=X_{k}a+zw$, so the residual-sum-of-squares change isThe map $(a_0,a_1,\dots ,a_k,z)\mapsto (a_0-mz,a_1,\dots ,a_k,z)$ is triangular with determinant one. Thus the reversible-jump Markov chain Monte Carlo ratio isFor an entirely explicit Gaussian expression, define $Q_k=(a-{\mu }_k{)}^T{\Sigma }_k^{-1}(a-{\mu }_k)$ and $Q_{k+1}=(a^{\prime }-{\mu }_{k+1}{)}^T{\Sigma }_{k+1}^{-1}(a^{\prime }-{\mu }_{k+1})$. Dividing the two normalized Gaussian priors and the proposal probability density function givesEqual model-order priors and equal birth/death selection probabilities remove the first ratio, if that is the intended convention. Dropping the order-dependent normal-prior constants instead would give the wrong acceptance rule.

From order $k+1$, set $z=a_{k+1}^{\prime }$, remove that final coefficient, retain $a_i=a_i^{\prime }$ for $1\le i\le k$, and restore the intercept byKeep $v$ fixed. This deterministically recovers the lower-order state and the auxiliary variable that would reproduce the higher-order state under a birth. The inverse map also has absolute Jacobian determinant one. In particular, simply deleting the highest coefficient without restoring the intercept is not the reverse of the supplied move.

For the recovered $a,z$ of the preceding part, compute the forward ratio $R$ from part (i). The reverse probability density function ratio and reverse Jacobian determinant invert every factor, soEquivalently, with $L_k(a;v)$ denoting the Gaussian regression likelihood,The factor $q(z)$ is present even though no random variable is generated during the death proposal: it is the probability density function of the auxiliary variable required by the reverse birth. These reciprocal rules give detailed balance on the paired dimension-matching moves.