Past exam of the mathematics course of the University of Cambridge / 2008 / iii / Paper 47 / 5 / a

A systematic Gibbs sampler updates the coordinates of $\theta$ successively from their full conditional distributions. At iteration $t$, drawFor a single coordinate, if the current state has distribution $\pi$, its untouched coordinates have the correct marginal law and the resampled coordinate has precisely the conditional law under $\pi$. Integrating marginal times conditional therefore recovers $\pi$. Each coordinate kernel preserves this invariant distribution, so their composition also preserves it. This proves the target of the Gibbs sampler; the whole systematic sweep need not be a reversible kernel.

Starting from a suitable state and assuming the required irreducibility and convergence conditions, discard an initial transient and retain subsequent states as a dependent sample from, or approximately from, $\pi$. In finite-state problems irreducibility and aperiodicity suffice; in continuous-state problems the corresponding positive recurrence and ergodicity conditions must hold. Successive draws are generally dependent, so Monte Carlo uncertainty must account for that dependence.

Use the shape-scale convention for the inverse-gamma distribution: $\mathrm{IG}(a,b)$ has density proportional to $s^{-a-1}{e}^{-b/s}$ for $s>0$. Write $X=X_k$, $\beta ={\beta }_k$, $\mu ={\mu }_k$, $\Sigma ={\Sigma }_k$ and $s={\sigma }^2$. The independent normal and inverse-gamma regression priors make the conditional log-likelihood plus log prior in $\beta$ equal, up to a constant, toDefine $V=(X^{T}X/s+{\Sigma }^{-1}{)}^{-1}$ and $m=V(X^{T}y/s+{\Sigma }^{-1}\mu )$. Completing the square turns this expression into $-(\beta -m{)}^T{V}^{-1}(\beta -m)/2$ plus a constant. Thus the conditional posterior distribution is the multivariate normal distributionPositive definiteness of the prior covariance matrix makes this well defined even if the design matrix is rank deficient.

Holding $\beta$ fixed, the factors depending on $s$ areConsequentlyThere is no additional $(k+1)/2$ in the shape: the normal distribution prior for ${\beta }_k$ is independent of ${\sigma }^2$. A prior scaled by ${\sigma }^2$ would be a different model and would change that calculation.

Choose an initial $s^{(0)}>0$. At each iteration, evaluate $V$ and $m$ using $s^{(t)}$ and draw ${\beta }^{(t+1)}$ from the conditional multivariate normal distribution just derived. Next compute $b_t=b+\Vert y-X_k{\beta }^{(t+1)}{\Vert }^2/2$ and drawFor example, generate $G\sim \Gamma (a+n/2,\text{rate }{b}_t)$ and set $s^{(t+1)}=1/G$. A draw of the conditional multivariate normal distribution can be formed as $m+LZ$ with $L{L}^T=V$ a Cholesky decomposition and $Z$ a vector of independent standard normal distributions.

This two-block Gibbs sampler alternates the full conditional distributions of the independent normal and inverse-gamma regression priors posterior. The conditional-integration argument proves that its invariant distribution is the joint posterior distribution. Under the usual convergence conditions, retain the pairs $({\beta }^{(t)},s^{(t)})$ after an initial transient as the requested dependent posterior sample. The new variance update uses the newly drawn coefficients, rather than the previous coefficients.