Introduction

Imagine you want to draw random samples from a distribution over many variables — say, the joint state of all pixels in an image, or the topic assignments for every word in a document. Writing down the full distribution is one thing; drawing from it is quite another.

In high dimensions, the joint density is almost never something you can sample directly. Integration is intractable, and rejection sampling wastes nearly every candidate. You need a cleverer way in.

Gibbs sampling offers a surprising answer: instead of tackling the whole distribution at once, fix all variables but one and sample that one from its conditional. Then move to the next variable, and the next, cycling endlessly. The resulting sequence of states forms a Markov chain, and — under mild conditions — its long-run distribution converges to exactly the joint distribution you wanted.

Named for physicist J. Willard Gibbs and introduced into statistical computation by Stuart Geman and Donald Geman in 1984 in the context of image restoration, Gibbs sampling became the workhorse of Bayesian statistics and probabilistic machine learning. It is one of the clearest examples of how a hard global problem dissolves into a sequence of easy local ones.

Watch the Cloud Fill In

Below is a correlated 2D Gaussian — a distribution where x and y are linked. Direct sampling would need the joint density; Gibbs sampling only ever needs the two conditionals: x given y, and y given x.

<div class="controls">
  <label>Correlation: <span id="rho-val">0.85</span>
    <input type="range" id="rho" min="-0.99" max="0.99" step="0.01" value="0.85">
  </label>
  <div class="btns">
    <button id="step-btn" type="button">Step</button>
    <button id="run-btn" type="button">Run</button>
    <button id="reset-btn" type="button" class="ghost">Reset</button>
  </div>
</div>
<canvas id="canvas" width="420" height="300"></canvas>
<div class="info" id="info">Press Step or Run to start sampling.</div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; color: #222; }
.controls { display: flex; flex-wrap: wrap; align-items: center; gap: .6rem; margin-bottom: .5rem; font-size: .9rem; }
label { display: flex; flex-direction: column; gap: .2rem; }
input[type=range] { width: 140px; cursor: pointer; }
.btns { display: flex; gap: .4rem; }
button { font: 600 13px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
canvas { border: 1px solid #cdd9e3; border-radius: 8px; background: #f8fafc; display: block; }
.info { font-size: .85rem; color: #444; margin-top: .4rem; min-height: 1.3em; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Gibbs sampling on a 2D correlated Gaussian.
// Target: (x,y) ~ N(0, Sigma), Sigma = [[1, rho],[rho, 1]]
// Conditionals:
//   x | y ~ N(rho*y, 1 - rho^2)
//   y | x ~ N(rho*x, 1 - rho^2)
// These are exact 1D Gaussians — easy to sample directly.
// No joint density ever needed.

const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');
const W = canvas.width, H = canvas.height;
const CX = W / 2, CY = H / 2;
const SCALE = 50; // pixels per unit std-dev

let rho = 0.85;
let samples = [];
let x = 0, y = 0;
let turn = 0; // 0 = update x, 1 = update y
let running = false;
let runTimer = null;
let totalSteps = 0;

const rhoSlider = document.getElementById('rho');
const rhoVal = document.getElementById('rho-val');
const info = document.getElementById('info');

rhoSlider.addEventListener('input', () => {
  rho = parseFloat(rhoSlider.value);
  rhoVal.textContent = rho.toFixed(2);
});

// Box-Muller transform: sample standard normal
function randn() {
  let u, v;
  do { u = Math.random(); } while (u === 0);
  v = Math.random();
  return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
}

// Convert distribution coords to canvas pixels
function toCanvas(dx, dy) {
  return [CX + dx * SCALE, CY - dy * SCALE];
}

function drawBackground() {
  ctx.clearRect(0, 0, W, H);
  // Draw faint grid lines
  ctx.strokeStyle = '#dde4ea';
  ctx.lineWidth = 0.8;
  for (let i = -4; i <= 4; i++) {
    const [px] = toCanvas(i, 0); const [, py] = toCanvas(0, i);
    ctx.beginPath(); ctx.moveTo(px, 0); ctx.lineTo(px, H); ctx.stroke();
    ctx.beginPath(); ctx.moveTo(0, py); ctx.lineTo(W, py); ctx.stroke();
  }
  // Axes
  ctx.strokeStyle = '#a0adb8'; ctx.lineWidth = 1.2;
  ctx.beginPath(); ctx.moveTo(0, CY); ctx.lineTo(W, CY); ctx.stroke();
  ctx.beginPath(); ctx.moveTo(CX, 0); ctx.lineTo(CX, H); ctx.stroke();
  // Axis labels
  ctx.fillStyle = '#7a8a99'; ctx.font = '11px system-ui';
  ctx.fillText('x', W - 14, CY - 6);
  ctx.fillText('y', CX + 6, 14);

  // Draw true ellipse (1-sigma and 2-sigma contours of the target Gaussian)
  // For N(0,Sigma) the 1-sigma ellipse satisfies x^2 - 2*rho*x*y + y^2 = (1-rho^2)
  // We draw it parametrically: the covariance matrix has eigenvalues 1±rho,
  const lambda1 = 1 + rho, lambda2 = 1 - rho;
  for (const scale of [1, 2]) {
    ctx.beginPath();
    for (let t = 0; t <= 2 * Math.PI + 0.01; t += 0.04) {
      // In eigenbasis: a*cos(t)*sqrt(lambda1), b*sin(t)*sqrt(lambda2)
      const a = scale * Math.sqrt(Math.max(lambda1, 0));
      const b = scale * Math.sqrt(Math.max(lambda2, 0));
      const u1 = a * Math.cos(t), u2 = b * Math.sin(t);
      // Rotate back (eigenvectors at 45 degrees)
      const dx = (u1 + u2) / Math.sqrt(2);
      const dy = (u1 - u2) / Math.sqrt(2);
      const [px, py] = toCanvas(dx, dy);
      t === 0 ? ctx.moveTo(px, py) : ctx.lineTo(px, py);
    }
    ctx.closePath();
    ctx.strokeStyle = scale === 1 ? 'rgba(29,53,87,0.25)' : 'rgba(29,53,87,0.12)';
    ctx.lineWidth = 1;
    ctx.stroke();
  }
}

function drawSamples() {
  const n = samples.length;
  samples.forEach(([sx, sy], i) => {
    const [px, py] = toCanvas(sx, sy);
    const age = (n - i) / Math.max(n, 1);
    const alpha = Math.max(0.12, 1 - age * 0.85);
    ctx.beginPath();
    ctx.arc(px, py, 2.5, 0, 2 * Math.PI);
    ctx.fillStyle = `rgba(230,57,70,${alpha.toFixed(2)})`;
    ctx.fill();
  });
}

function drawCurrent() {
  if (samples.length === 0) return;
  const [px, py] = toCanvas(x, y);
  ctx.beginPath(); ctx.arc(px, py, 5, 0, 2 * Math.PI);
  ctx.fillStyle = '#1d3557'; ctx.fill();
}

function redraw() {
  drawBackground();
  drawSamples();
  drawCurrent();
}

function gibbsStep() {
  const sd = Math.sqrt(Math.max(1 - rho * rho, 1e-9));
  if (turn === 0) {
    // Sample x | y
    x = rho * y + sd * randn();
  } else {
    // Sample y | x
    y = rho * x + sd * randn();
    samples.push([x, y]);
    if (samples.length > 500) samples.shift(); // keep display manageable
  }
  turn = 1 - turn;
  totalSteps++;
}

function step() {
  gibbsStep();
  redraw();
  info.textContent = fmt('Steps: %totalSteps% | Samples: %length% | Current: (%x%, %y%)', { totalSteps: totalSteps, length: samples.length, x: x.toFixed(2), y: y.toFixed(2) });
}

function startRun() {
  if (running) return;
  running = true;
  document.getElementById('run-btn').textContent = 'Pause';
  runTimer = setInterval(() => {
    for (let i = 0; i < 4; i++) gibbsStep();
    redraw();
    info.textContent = fmt('Steps: %totalSteps% | Samples: %length% | Current: (%x%, %y%)', { totalSteps: totalSteps, length: samples.length, x: x.toFixed(2), y: y.toFixed(2) });
  }, 40);
}

function pauseRun() {
  running = false;
  clearInterval(runTimer);
  document.getElementById('run-btn').textContent = 'Run';
}

function reset() {
  pauseRun();
  samples = []; x = 0; y = 0; turn = 0; totalSteps = 0;
  redraw();
  info.textContent = 'Press Step or Run to start sampling.';
}

document.getElementById('step-btn').addEventListener('click', () => { pauseRun(); step(); });
document.getElementById('run-btn').addEventListener('click', () => { running ? pauseRun() : startRun(); });
document.getElementById('reset-btn').addEventListener('click', reset);

// Verify: conditional means must satisfy E[x|y] = rho*y and E[y|x] = rho*x
// Quick numerical check: run 10000 steps, verify sample mean near (0,0)
// and sample correlation near rho. (Logged to console.)
(function verify() {
  let vx = 0, vy = 0, vt = 0, sx2 = 0, sy2 = 0, sxy = 0;
  let tx = 0, ty = 0, tt = 0;
  const vRho = 0.7;
  const vsd = Math.sqrt(1 - vRho * vRho);
  for (let i = 0; i < 20000; i++) {
    if (tt === 0) { tx = vRho * ty + vsd * randn(); }
    else { ty = vRho * tx + vsd * randn(); vx += tx; vy += ty; vt++; sx2 += tx*tx; sy2 += ty*ty; sxy += tx*ty; }
    tt = 1 - tt;
  }
  const mx = vx/vt, my = vy/vt;
  const rx = Math.sqrt(sx2/vt - mx*mx), ry = Math.sqrt(sy2/vt - my*my);
  const corr = (sxy/vt - mx*my) / (rx * ry);
  const ok = Math.abs(mx) < 0.05 && Math.abs(my) < 0.05 && Math.abs(corr - vRho) < 0.05;
  console.log(fmt('[gibbs-verify] mean=(%mx%,%my%) corr=%corr% target_rho=%vRho% PASS=%ok%', { mx: mx.toFixed(3), my: my.toFixed(3), corr: corr.toFixed(3), vRho: vRho, ok: ok }));
})();

redraw();
Result
Console

Press Step to take one Gibbs update (alternating x then y). Press Run to let it go. Each new dot is drawn from a 1D conditional — no joint density ever computed. Watch how the elongated cloud builds up from scratch, with early samples (lighter) leading the chain toward the dense region.

The burn-in phase (the chain wandering to the high-probability region) and mixing (how quickly the chain explores the full distribution) are the two central concerns of any MCMC method. Adjust the correlation slider to make sampling harder and watch mixing slow down.

Convergence and Mixing

Why does endlessly cycling one variable at a time actually converge to the right distribution? The answer comes from Markov chain theory.

  • Detailed balance. Each Gibbs update satisfies a condition called detailed balance: transitions between any two states are balanced so that the target joint distribution is the stationary distribution of the chain. The chain is guaranteed to end up there.
  • Ergodicity. If the chain can reach every part of the distribution from anywhere (irreducibility) and returns to states regularly (aperiodicity), convergence is guaranteed — these conditions usually hold when the conditionals have full support.
  • Mixing time. Convergence is guaranteed, but speed is not. When variables are highly correlated, a single-coordinate update barely moves the chain — each step changes only one dimension while the other holds it back. This is called slow mixing. Near-deterministic relationships between variables can trap the chain for exponentially many steps.
  • Block Gibbs. A classic fix is to group correlated variables and sample them jointly in one block, eliminating the internal bottleneck. The trade-off is that sampling the block may itself be hard.

Gibbs sampling has no formal NP-hard barrier of its own — the per-step cost is polynomial as long as the conditionals are easy to sample. But mixing can be exponentially slow, and distinguishing "converged" from "just stuck" is undecidable in general. In practice, practitioners run multiple chains from different starting points and look for agreement as the main convergence diagnostic.

Gibbs sampling sits alongside other MCMC methods like Metropolis–Hastings in the landscape of Bayesian inference. Where Metropolis–Hastings accepts or rejects proposed moves, Gibbs always accepts — its proposals are already drawn from the exact conditional.

Where It Matters

"Sample from the conditional, one variable at a time" turns out to be useful almost everywhere probabilistic models live:

  • Bayesian inference: the engine behind Bayesian hierarchical models and posterior computation in statistics packages like BUGS, Stan, and PyMC.
  • Latent Dirichlet Allocation (LDA): the topic model that assigns each word to a topic is trained by cycling through every word and resampling its topic from the conditional given all other assignments — pure collapsed Gibbs sampling.
  • Image restoration: Geman & Geman's original 1984 paper used Gibbs sampling to reconstruct noisy images by iteratively updating each pixel's label given its neighbors' current labels.
  • Boltzmann machines: training the early neural energy models relied on Gibbs sampling to estimate expectations under the model distribution.
  • Statistical genetics: genome-wide association studies use Gibbs samplers to impute missing genotype data and estimate haplotype frequencies.
  • Natural language processing: word alignment models in machine translation and coreference resolution both use Gibbs as a core inference engine.

Wherever you have a probabilistic model with many interdependent variables and tractable conditionals, Gibbs sampling turns the problem of global inference into a loop of local coin flips — the same theme as Bayesian inference and k-means clustering.

Conclusion

Gibbs sampling encodes a beautiful insight: a joint distribution that is impossible to sample directly becomes approachable if you are willing to be patient and sequential. Fix everything, update one variable, repeat.

The mathematics guarantees correctness. The practice demands vigilance about mixing — slow-mixing chains produce convincing-looking but wrong samples. And the applications span a remarkable range, from image pixels to document topics to human genomes.

Next time you see a probabilistic model with dozens of latent variables, remember that somewhere underneath it, a Gibbs sampler is probably cycling through them one by one — and that this humble loop is doing the work that no closed-form formula could.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/gibbs-sampling/Content licensed under CC BY-NC 4.0.