Solution (source code)

= Solution

For the <Metropolis–Hastings algorithm>, start at a point $x$ with positive target <probability density function>. Draw $Y$ from a proposal <probability density function> $q(y\mid x)$, generate an <independent> <uniform distribution> value $U$, and set the next state to $Y$ if
$$
U\leq\alpha(x,Y),\qquad
\boxed{\alpha(x,y)=\min\!\left\{1,\frac{\pi(y)q(x\mid y)}{\pi(x)q(y\mid x)}\right\}.}
$$
Otherwise retain $x$, including that repeated state in the sample. Where a proposed forward move is possible but its reverse <probability density function> vanishes, its acceptance <probability> is zero. Unknown <normalizing constants> in $\pi$ cancel.

For distinct states, the accepted <probability> flow is
$$
\pi(x)q(y\mid x)\alpha(x,y)
=\min\{\pi(x)q(y\mid x),\pi(y)q(x\mid y)\},
$$
which is symmetric in $x,y$. This <detailed balance> identity, together with the holding <probability>, proves that $\pi$ is invariant. Under appropriate irreducibility and aperiodicity conditions the chain converges to $\pi$; successive states are generally dependent.

The <Gibbs sampler> updates coordinates from their <full conditional distributions>. In a systematic sweep, draw successively
$$
X_j^{(r)}\sim
\pi\!\left(\,\cdot\mid X_1^{(r)},\ldots,X_{j-1}^{(r)},
X_{j+1}^{(r-1)},\ldots,X_k^{(r-1)}\right),
\quad j=1,\ldots,k.
$$
Each update preserves the target joint law: integrating the old coordinate against its conditional <probability density function> and replacing it by an <independent> draw from that same conditional leaves the joint <probability density function> unchanged. The composition of the coordinate updates therefore also preserves $\pi$. A <Gibbs sampler> coordinate move can be regarded as a <Metropolis–Hastings algorithm> move with acceptance <probability> one. Systematic sweeps preserve <stationarity> but need not themselves be reversible. Initialization away from <stationarity> requires convergence before treating the draws as approximately target-distributed, and dependence must be accounted for in <Monte Carlo method> error estimates.