Question 9. Gibbs and random-walk Metropolis–Hastings

Topic: Markov chains, Gibbs sampling, and Metropolis–Hastings · more questions on this topic
‹ Question 6 Question 13 ›

Consider the joint density on x>0, y>0,

f(x,y)\propto x^{3}ye^{-xy-2x-3y}.

(a) Describe a Gibbs sampler for this density. [10 marks]

(b) Describe a random-walk Metropolis–Hastings algorithm with increment distribution \operatorname{N}(0,\sigma^{2}I_2). [10 marks]

(c) Describe how to obtain a 95\% confidence interval for \mathbb{E}(Y). [5 marks]

Worked answer (click to show)

Definitions. The target is the normalised version of q(x,y)=x^3y\exp(-xy-2x-3y) on the positive quadrant. A \operatorname{Gamma}(a,b) full conditional has shape a and rate b. Gibbs sampling redraws coordinates from exact full conditionals; Metropolis–Hastings proposes a move and corrects it with an acceptance probability. In part (c), the target is the fixed expectation \mu_Y=\mathbb{E}_f(Y), and the interval measures numerical error in estimating it.

(a) Given x, the factors in y are y\,e^{-(x+3)y}, so

Y\mid x\sim\operatorname{Gamma}(2,\,x+3).

Given y, the factors in x are x^{3}e^{-(y+2)x}, so

X\mid y\sim\operatorname{Gamma}(4,\,y+2).

Each draw uses the newest available value and comes from an exact full conditional, so it needs no acceptance step. The pairs across sweeps are generally dependent. Both gamma distributions use rates, not scales.

(b) From the current point (x,y), draw an increment Z\sim \operatorname{N}(0,\sigma^{2}I_2) and set (x',y')=(x,y)+Z. This proposal is symmetric. If x'\leq 0 or y'\leq 0, reject. Otherwise accept with probability

\min\Bigl(1,\frac{f(x',y')}{f(x,y)}\Bigr) = \min\Bigl(1,\frac{(x')^{3}y'}{x^{3}y} \exp\bigl(-(x'y'-xy)-2(x'-x)-3(y'-y)\bigr)\Bigr).

Draw an independent U\sim\operatorname{Uniform}(0,1) and accept the valid proposal if U is no greater than this probability; otherwise retain and record the current point. A proposal outside the support is rejected, not redrawn until positive, because repeated redrawing would change the proposal kernel. The target’s unknown normalising constant cancels from the ratio.

Symmetry means equality of the forward and reverse density values: a centred normal increment has the same density at Z and -Z. Hence the proposal ratio in the general Metropolis–Hastings formula is one. The scale \sigma^2 affects exploration, not the target distribution: very small moves can mix slowly, while very large moves can cause many rejections. Repeated states after rejection are part of the chain and must be retained in averages.

For a general proposal density g(u'\mid u), the acceptance ratio includes g(u\mid u')/g(u'\mid u). The following algorithm uses the symmetric proposal in this question, so that factor cancels. Define the log kernel h(x,y)=3\log x+\log y-xy-2x-3y on the positive quadrant. Evaluating acceptance on the log scale avoids exponentiating a potentially extreme density ratio.

An invalid proposal or failed acceptance test leaves the current state unchanged. In particular, do not keep proposing until a move is accepted, and do not save only accepted moves: either action changes the sampling procedure. Check traces, exploration and autocorrelation before using the retained states in part (c).

(c) After checking convergence and allowing suitable warm-up, retain N=ab draws of Y in a non-overlapping batches of length b. Let \bar Y_N be their overall mean and \bar Y_k^{\mathrm{batch}} the mean in batch k. With

s_B^2=\frac{1}{a-1}\sum_{k=1}^a \bigl(\bar Y_k^{\mathrm{batch}}-\bar Y_N\bigr)^2, \qquad \widehat{\operatorname{MCSE}}(\bar Y_N)=\frac{s_B}{\sqrt a},

an approximate 95\% confidence interval for the fixed target \mathbb{E}(Y) is \bar Y_N\pm t_{0.975,a-1}s_B/\sqrt a (or 1.96 in place of the t quantile for many batches). Batches must be long enough to weaken within-chain dependence. The iid standard error s_Y/\sqrt N ignores autocorrelation, and quantiles of the Y-draws describe Y, not uncertainty in its mean.

The reason for batching is that Monte Carlo uncertainty depends on the whole autocovariance sequence. Under a stationary Markov-chain central limit theorem,

\sqrt N(\bar Y_N-\mu_Y)\ \overset{d}{\longrightarrow}\ \operatorname{N}(0,\tau_Y^2), \qquad \tau_Y^2=\gamma_0+2\sum_{k=1}^{\infty}\gamma_k,

where \gamma_k=\operatorname{Cov}(Y_t,Y_{t+k}). Batch means estimates \tau_Y^2 by bs_B^2, giving \sqrt{bs_B^2/N}=s_B/\sqrt a. There must be enough batches to estimate their variance and enough observations within each batch to account for dependence. Check that the reported error is reasonably stable as the simulation and batch length are increased; discarding warm-up alone does not remove autocorrelation.