Gibbs Sampling For Gaussian Regression

Why Use Gibbs Sampling for Gaussian Regression?

In Gaussian regression with a normal-inverse-gamma prior, the posterior distributions are available in closed form. That makes direct posterior calculation possible, but it does not mean that every calculation of interest is convenient, especially in hierarchical models or when a larger algorithm needs repeated draws. Gibbs sampling generates a sequence of parameter values by drawing from simpler conditional posterior distributions in turn. For this regression model, each update is a standard distribution and can be computed from the data and the current value of the other parameter block.

The goal here is to construct those updates correctly and use a simulated chain to estimate posterior summaries. The draws in a Gibbs chain are generally dependent, so they are not a set of independent observations from the posterior. That distinction matters when interpreting summaries and their numerical precision.

Model and Prior

Model and prior: Let \(y\) be an \(n\times 1\) outcome vector and \(X\) an \(n\times k\) design matrix. The Gaussian regression model is \(y\mid\beta,\sigma^2\sim N(X\beta,\sigma^2 I_n)\). The normal-inverse-gamma prior is \(\beta\mid\sigma^2\sim N(b_0,\sigma^2 V_0)\) and \(\sigma^2\sim \operatorname{IG}(a_0,d_0)\). Here \(V_0\) is positive definite, \(a_0>0\), and \(d_0>0\). The inverse-gamma convention used here has density proportional to \((\sigma^2)^{-a_0-1}\exp(-d_0/\sigma^2)\).

The prior parameters \(b_0\) and \(V_0\) describe the prior mean and variance structure for the coefficients, conditional on the error variance. The parameter \(d_0\) is the scale in the stated inverse-gamma convention; it is not the rate of a gamma distribution for \(\sigma^2\). The formulas below use the entire design matrix, including any intercept column.

Full Conditional Distributions

A Gibbs sampler updates one block of unknowns at a time from its posterior distribution conditional on the data and the current values of the other block. For this model, the coefficient update is normal and the variance update is inverse-gamma.

Theorem: Define \(A=V_0^{-1}+X^\top X\) and \(m=A^{-1}(V_0^{-1}b_0+X^\top y)\). Then the full conditional distributions are \(\beta\mid\sigma^2,y\sim N(m,\sigma^2 A^{-1})\) and \(\sigma^2\mid\beta,y\sim \operatorname{IG}(a_0+(n+k)/2,\ d_0+\tfrac12[(y-X\beta)^\top(y-X\beta)+(\beta-b_0)^\top V_0^{-1}(\beta-b_0)])\).

The coefficient update combines the information in the likelihood, \(X^\top X\), with the prior precision, \(V_0^{-1}\). The variance update combines the squared regression residuals with the prior-weighted squared distance of \(\beta\) from \(b_0\). The latter term is required because the prior distribution of \(\beta\) itself depends on \(\sigma^2\).

Why These Updates Have These Forms

For the coefficient conditional, hold \(\sigma^2\) fixed and multiply the likelihood kernel by the prior kernel for \(\beta\). Ignoring factors that do not depend on \(\beta\), the exponent is \(-\tfrac{1}{2\sigma^2}[(y-X\beta)^\top(y-X\beta)+(\beta-b_0)^\top V_0^{-1}(\beta-b_0)]\). Expanding the quadratic terms in \(\beta\) gives the quadratic precision \(A=V_0^{-1}+X^\top X\) and linear term \(V_0^{-1}b_0+X^\top y\). Since \(A\) is positive definite, completing the square gives

\beta^\top A\beta-2\beta^\top(V_0^{-1}b_0+X^\top y) =(\beta-m)^\top A(\beta-m)-m^\top A m.

The final term on the right does not depend on \(\beta\). The remaining kernel is that of a normal distribution with mean \(m\) and covariance \(\sigma^2 A^{-1}\), which establishes the first conditional.

For the variance conditional, hold \(\beta\) fixed. The likelihood contributes a factor \((\sigma^2)^{-n/2}\) and an exponential involving the residual sum of squares. The conditional normal prior for \(\beta\) contributes \((\sigma^2)^{-k/2}\) as a density in \(\beta\), as well as its prior quadratic term. The prior for \(\sigma^2\) contributes \((\sigma^2)^{-a_0-1}\exp(-d_0/\sigma^2)\). Collecting the powers gives the inverse-gamma shape \(a_0+(n+k)/2\). The coefficient-prior quadratic contributes to the scale, while its normalization contributes \(k/2\) to the shape. The terms in the exponent combine as the stated scale divided by \(\sigma^2\). This is why the coefficient-prior penalty appears in the scale update.

Worked Example: Reading the Coefficient Conditional

Suppose a model has two coefficients, and its prior covariance multiplier is \(V_0=I_2\). The data and prior mean produce \(X^\top X=\begin{pmatrix}4&1\\1&3\end{pmatrix}\), \(X^\top y=\begin{pmatrix}5\\2\end{pmatrix}\), and \(b_0=\begin{pmatrix}0\\0\end{pmatrix}\). Find the coefficient conditional mean and covariance multiplier.

First, \(A=I_2+X^\top X=\begin{pmatrix}5&1\\1&4\end{pmatrix}\). Its determinant is \(5(4)-1(1)=19\), so \(A^{-1}=\frac1{19}\begin{pmatrix}4&-1\\-1&5\end{pmatrix}\).

The conditional mean is \(m=A^{-1}X^\top y\), because \(V_0^{-1}b_0=0\). Thus \(m=\frac1{19}\begin{pmatrix}4&-1\\-1&5\end{pmatrix}\begin{pmatrix}5\\2\end{pmatrix} =\frac1{19}\begin{pmatrix}18\\5\end{pmatrix} =\begin{pmatrix}0.9474\\0.2632\end{pmatrix}\).

Therefore, at a current variance value \(\sigma^2\), draw \(\beta\) from \(N(m,\sigma^2 A^{-1})\). The matrix \(A^{-1}\) is the covariance multiplier; it is not itself the conditional covariance unless \(\sigma^2=1\).

A Gibbs Sweep

A single sweep consists of drawing the coefficients conditional on the current variance, then drawing the variance conditional on the newly updated coefficients. Repeating this procedure yields a Markov chain whose target distribution is the joint posterior. The transition is not an instruction to draw both blocks independently from their marginal posteriors; each draw uses the current value of the other block.

Choose starting values.
Set an initial \(\sigma^{2(0)}>0\) and an initial coefficient vector \(\beta^{(0)}\). These values start the computation; they are not posterior estimates.

Update the coefficients.
At iteration \(t\), draw \(\beta^{(t)}\) from \(N(m,\sigma^{2(t-1)}A^{-1})\).

Update the variance.
Draw \(\sigma^{2(t)}\) from the inverse-gamma conditional using \(\beta^{(t)}\).

Retain draws for summaries.
After the chosen initial iterations, store the parameter draws used to estimate posterior quantities. The retained values remain a dependent chain.

The joint posterior is invariant under each full-conditional update: if the current pair has the joint posterior distribution, drawing one parameter block from its conditional distribution given the other preserves that joint distribution. Alternating the two updates therefore has the posterior as its stationary target. This invariance explains the construction, but does not by itself guarantee that a finite run has explored the target adequately.

Worked Example: One Numerical Gibbs Update

Consider hypothetical data with \(x=(-1,0,1)^\top\), outcomes \(y=(1,2,4)^\top\), and design matrix \(X\) whose columns are an intercept and \(x\). Use \(b_0=(0,0)^\top\), \(V_0=I_2\), \(a_0=2\), and \(d_0=1\). Suppose the current variance is \(\sigma^{2(0)}=1.5\).

The design cross-products are \(X^\top X=\begin{pmatrix}3&0\\0&2\end{pmatrix}\) and \(X^\top y=(7,3)^\top\): the first entry is \(1+2+4=7\), and the second is \((-1)(1)+0(2)+1(4)=3\). Hence \(A=\operatorname{diag}(4,3)\) and \(m=(7/4,3/3)^\top=(1.75,1)^\top\). The coefficient conditional covariance is \(1.5A^{-1}=\operatorname{diag}(0.375,0.5)\).

For illustration, suppose the normal innovations used in this draw are \(z=(1/\sqrt{6},0)^\top\), with independent standard normal components. A draw can be written \(m+\sqrt{\sigma^{2(0)}}A^{-1/2}z\). The intercept increment is \(\sqrt{0.375}/\sqrt6=0.25\), and the slope increment is zero. Thus the realized coefficient draw is \(\beta^{(1)}=(2,1)^\top\).

Using this updated vector, the fitted values are \(1,2,3\), so the residuals are \(0,0,1\) and the residual sum of squares is \(1\). The prior penalty is \((\beta^{(1)}-b_0)^\top V_0^{-1}(\beta^{(1)}-b_0)=2^2+1^2=5\). Therefore the variance conditional has shape \(2+(3+2)/2=4.5\) and scale \(1+\tfrac12(1+5)=4\). The next step is a draw from \(\operatorname{IG}(4.5,4)\); its particular realized value depends on the random-number draw.

Estimating Posterior Summaries From a Chain

For a scalar parameter \(\theta\), a posterior expectation can be approximated by the average of its retained draws. If \(\theta^{(1)},\ldots,\theta^{(M)}\) are retained values, the Monte Carlo estimate is \(\widehat{E}[\theta\mid y]=M^{-1}\sum_{t=1}^{M}\theta^{(t)}\). Posterior probabilities can be estimated by the fraction of draws satisfying the event. Empirical quantiles of the draws approximate posterior quantiles and can define credible intervals.

These are simulation-based approximations, not exact values. Because successive Gibbs draws can be correlated, the number of retained draws alone does not tell how much independent information the chain contains. A short chain may give unstable summaries even when its arithmetic is correct.

Worked Example: Summarizing Simulated Slope Draws

Suppose a hypothetical Gibbs run retains the following ten draws of a regression slope, in outcome units per one-unit increase in the predictor: \(0.72, 0.91, 1.08, 0.83, 1.16, 0.95, 0.67, 1.02, 0.88, 1.28\).

Their sum is \(9.50\), so the estimated posterior mean slope is \(9.50/10=0.950\) outcome units per predictor unit. All ten draws are positive, giving a chain-based estimate of \(10/10=1.00\) for the posterior probability that the slope is positive. This is an estimate from this finite chain, not proof that the exact posterior probability equals one.

Sorted, the draws are \(0.67,0.72,0.83,0.88,0.91,0.95,1.02,1.08,1.16,1.28\). Using the nearest-rank empirical quantile convention, the 10th percentile is the first ordered value and the 90th percentile is the ninth. The resulting approximate 80% equal-tail credible interval is \([0.67,1.16]\) outcome units per predictor unit. Under the specified model and prior, this interval is intended to contain 80% of the posterior probability, subject to Monte Carlo approximation; it is not a statement that there is an 80% frequentist probability that the fixed slope lies in this already computed interval.

Common Mistakes and Practical Checks

Dropping the prior penalty from the variance update. The scale includes both the residual sum of squares and the prior quadratic \((\beta-b_0)^\top V_0^{-1}(\beta-b_0)\). Omitting the second term changes the full conditional.
Using the wrong current value. In a sweep, the variance update uses the newly drawn \(\beta^{(t)}\), not \(\beta^{(t-1)}\). The coefficient update uses \(\sigma^{2(t-1)}\).
Confusing covariance with covariance multiplier. The coefficient conditional covariance is \(\sigma^2 A^{-1}\), not just \(A^{-1}\).
Treating chain draws as independent. A retained Gibbs sequence is generally autocorrelated. Its sample average estimates a posterior mean, but the usual independent-sample precision calculation does not automatically apply.
Interpreting a finite-run fraction as exact. If every retained draw has a positive coefficient, the estimated posterior probability is 1.00 for that run; the exact posterior probability need not be exactly one.

Gibbs sampling is especially useful when the conditional distributions are easy to sample even if a direct draw from the joint posterior is inconvenient. The normal-inverse-gamma structure makes each regression update explicit and provides checks on implementations: verify matrix dimensions, ensure \(A\) is positive definite, use the stated inverse-gamma parameterization, and recompute the residual and prior quadratic terms from the current coefficient draw.

Check Your Understanding

Use the model and notation above to answer the following questions.

Which quantities determine the conditional mean \(m\) of the coefficient update, and which additional quantity determines its covariance?
Why does the conditional variance scale include a quadratic penalty involving \(\beta-b_0\)?
In a Gibbs sweep, which coefficient vector is used to calculate the variance update at iteration \(t\)?
If the retained slope draws are all positive, what does the chain-based fraction estimate, and what does it not establish?
Explain why \(A^{-1}\) alone is not the conditional covariance of \(\beta\) unless \(\sigma^2=1\).