Past exam of the mathematics course of the University of Cambridge / 2007 / iii / Paper 48 / 5

Both the Metropolis–Hastings algorithm and the Gibbs sampler construct a Markov chain with a specified target distribution, usually a posterior distribution known only up to normalization. For a proposal density $q(y\mid x)$ and target density $\pi$, the Metropolis–Hastings acceptance probability isAfter drawing $y$, move there with this probability; otherwise stay at $x$. The unknown normalizing constant of $\pi$ cancels. The original Metropolis algorithm uses a symmetric proposal, $q(y\mid x)=q(x\mid y)$, so the proposal-density ratio disappears. The general Metropolis–Hastings algorithm permits nonsymmetric proposals and needs that correction.

A Gibbs sampler instead chooses a coordinate or block and draws directly from its full conditional distribution, leaving the other coordinates fixed. It is exactly a conditional Metropolis–Hastings update whose proposal is that full conditional. For a change from $x_i$ to $y_i$ with $x_{-i}$ fixed, the acceptance ratio isThus a Gibbs update is an always-accepted Metropolis–Hastings update with the full conditional as proposal. This statement concerns each update: a systematic sweep of several Gibbs kernels preserves the target but need not itself be reversible.

When the full conditional distributions are easy to sample, a Gibbs sampler avoids proposal tuning and rejection, and is particularly convenient for conjugate priors and data augmentation. When a conditional cannot be sampled directly, a Metropolis–Hastings algorithm can use only its evaluable density; such an update can also be inserted into a Gibbs sweep. Block or joint Metropolis–Hastings proposals can move efficiently along strongly correlated directions for which coordinate-wise Gibbs updates mix slowly. Conversely, a badly chosen Metropolis–Hastings proposal can reject most moves or make only tiny moves. Acceptance probability one by itself does not guarantee rapid exploration. Both methods require appropriate irreducibility and convergence checks before their long-run averages are useful posterior estimates.

Write $S_m={\sum }_{i=1}^m{x}_i$ and $S=S_N$. Use shape-rate notation for the Gamma distribution, as in the authoritative PDF: shape $a$ and rate $b$ give a density proportional to $z^{a-1}{e}^{-bz}$. Conditional independence of the Poisson distributions gives a likelihood function proportional toThe product of factorial terms depends only on the data and can be absorbed into the normalizing constant. Multiplying by the two independent Gamma distribution prior kernels and the uniform prior on the possible indices gives the Poisson count change-point model with gamma priors:Because the hyperparameters are positive, each possible index gives finite positive gamma integrals, and summing over the finitely many indices yields a proper posterior distribution. At $m=N$ the second likelihood factor is one, so the posterior law of its unused rate $\varphi$ remains its prior $\mathrm{Gamma}(\gamma ,\delta )$.

The posterior kernel gives two gamma full conditional distributions directly:with rates as the second parameters. For the index, terms independent of the proposed value $j$ cancel. Its conditional weights areDefine ${\ell }_j=j(\varphi -\lambda )+S_j\log (\lambda /\varphi )$ and ${\ell }_{\max }={\max }_j{\ell }_j$. The normalized categorical probabilities areSubtracting the largest log weight prevents needless numerical overflow, an application of the log-sum-exp function.

Precompute the cumulative counts $S_j$. Initialize any $m\in \{1,\dots ,N\}$ and positive rates. At each iteration, draw $\lambda$ from its gamma conditional given the current index, draw $\varphi$ from its gamma conditional given that index, and draw a new index from the categorical probabilities using the newly drawn rates. Repeat and retain the parameter triples after initialization effects have diminished. These are the Gibbs updates for a Poisson count change point.

To verify invariance rather than just identify a sampler, a coordinate update integrates the old joint density as $\pi (\lambda ,\varphi ,m\mid x)=\pi (\lambda \mid \varphi ,m,x)\pi (\varphi ,m\mid x)$ and replaces the first factor by an independent draw from the same conditional. The joint posterior distribution is therefore unchanged. The same argument applies to the other two updates, so their composition preserves the posterior. All gamma densities and all categorical weights are positive on their supports. The full-sweep transition density is consequently positive across the joint support, giving the usual irreducibility and aperiodicity conditions for this proper-posterior Gibbs sampler. This supports posterior long-run averages under the ergodic theorem for a positive Harris recurrent Markov chain.

The event (m=N) represents no change within the observed sequence: every observation belongs to the first-rate segment, and the second-rate segment is empty. It is not the event $\lambda =\varphi$; under continuous independent rate priors that exact equality has posterior probability zero. The endpoint is instead a distinct index value with positive prior mass.

For retained Markov chain Monte Carlo draws $({\lambda }^{(b)},{\varphi }^{(b)},m^{(b)})$, the direct estimate isIts consistency follows from the ergodic theorem for a positive Harris recurrent Markov chain. The draws are correlated, so an independent-binomial standard error would generally be incorrect; uncertainty should account for the Markov chain Monte Carlo asymptotic variance.

A conditional Monte Carlo alternative averages the conditional probability from the preceding solution:The tower property of conditional expectation shows that its posterior expectation is the same target. This replaces an indicator by its conditional expectation; for correlated chains one should assess efficiency rather than automatically assume its entire asymptotic variance is smaller.

An exact finite-sum check is available by integrating out both rates. The marginal posterior weights of a Poisson count change point areThese follow from the Gamma function integral ${\int }_0^{\infty }{z}^{a-1}{e}^{-bz}dz=\Gamma (a)/b^a$. All omitted prior-normalization and factorial factors are independent of $j$, so they cancel from this ratio. This exact calculation provides a useful independent check on either simulation estimate.