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

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.

A reversible-jump Markov chain Monte Carlo calculation compares normalized model-specific densities. The fixed-order posterior written only up to proportionality in part (a) omits normalizing factors that depend on $k$ and therefore cannot simply be canceled between dimensions. Moreover, the printed question does not specify the model-order prior distribution or the probabilities of selecting birth and death moves. Denote those by $p_k$, $b_k$ and $d_k$, respectively, retaining them in the answer.

Put $s={\sigma }^2$, $c=n^{-1}{\sum }_i{x}_i^{k+1}$, $v_i=x_i^{k+1}-c$, and $r=y-X_k\beta$. The centered polynomial birth move in reversible-jump sampling changes the fitted vector by $zv$, so its new residual is $r^{\prime }=r-zv$. ThereforeThe map from $({\beta }_0,\dots ,{\beta }_k,z)$ to $({\beta }_0^{\prime },\dots ,{\beta }_{k+1}^{\prime })$ has ones down the diagonal and a single additional entry $-c$ in the intercept row. Its Jacobian determinant is one.

Write $Q_k=(\beta -{\mu }_k{)}^T{\Sigma }_k^{-1}(\beta -{\mu }_k)$ and $Q_{k+1}=({\beta }^{\prime }-{\mu }_{k+1}{)}^T{\Sigma }_{k+1}^{-1}({\beta }^{\prime }-{\mu }_{k+1})$. The normalized multivariate normal distribution prior ratio isDivide also by the forward proposal density $q(z)=(2\pi {\sigma }_{\beta }^2{)}^{-1/2}{e}^{-z^2/(2{\sigma }_{\beta }^2)}$. The common inverse-gamma distribution prior for $s$ and the common normal-likelihood normalization cancel, since $s$ stays fixed. The forward acceptance ratio is consequentlyHere ${\sigma }_{\beta }>0$ is the proposal standard deviation. If equal model priors and symmetric birth/death selection are explicitly adopted, their ratios become one; they are not determined by the PDF. The expression assumes the same hyperparameters for the variance prior in both models, as in part (a). A dimension-dependent variance prior would contribute its own ratio.

For the reverse centered polynomial birth move in reversible-jump sampling, recover $z={\beta }_{k+1}^{\prime }$ from the current coefficient vector. Remove that highest-order coefficient and undo the intercept shift:No auxiliary random number is drawn for this death move: the removed coefficient supplies the auxiliary value needed to reconstruct the birth move. The reverse Jacobian determinant is also one. The inverse prediction change is $-zv$, hence the residual in the lower-order model is $r=r^{\prime }+zv$.

Reversible-jump Markov chain Monte Carlo requires the reciprocal density/proposal ratio for the matched reverse transformation, so the death acceptance probability isExpressed directly in the current residual $r^{\prime }=y-X_{k+1}{\beta }^{\prime }$ and the recovered lower-order parameters, this isIndeed $r=r^{\prime }+zv$ makes this exponent the negative of the birth exponent, and all prefactors invert. Thus $R_{-}{R}_{+}=1$ for matched states. This reciprocity, including the normalized prior distributions and auxiliary proposal density, gives the required detailed balance between models.