Bayesian Linear Regression

Bayesian Linear Regression

The Bayesian linear regression model is
$\mathbf{y}\mid \mathbf{\beta },{\sigma }^2\sim N(\mathbf{X}\mathbf{\beta },{\sigma }^2{\mathbf{I}}_n),\mathbf{\beta }\mid {\sigma }^2\sim N({\mathbf{\beta }}_0,{\sigma }^2{\mathbf{\Lambda }}_0^{-1}),{\sigma }^2\sim \text{InvGamma}(a_0,b_0),$
where $\mathbf{y}\in {\mathbb{R}}^n$, $\mathbf{X}\in {\mathbb{R}}^{n\times d}$, and the prior is the Normal-Inverse-Gamma (NIG) family scaled by ${\sigma }^2$. The conjugate posterior is again NIG with
${\mathbf{\Lambda }}_n={\mathbf{X}}^{\top }\mathbf{X}+{\mathbf{\Lambda }}_0,{\mathbf{\beta }}_n={\mathbf{\Lambda }}_n^{-1}\left({\mathbf{X}}^{\top }\mathbf{y}+{\mathbf{\Lambda }}_0{\mathbf{\beta }}_0\right),$
$a_n=a_0+\frac{n}{2},b_n=b_0+\frac{1}{2}\left({\mathbf{y}}^{\top }\mathbf{y}+{\mathbf{\beta }}_0^{\top }{\mathbf{\Lambda }}_0{\mathbf{\beta }}_0-{\mathbf{\beta }}_n^{\top }{\mathbf{\Lambda }}_n{\mathbf{\beta }}_n\right).$
Marginal posterior for $\mathbf{\beta }$ (integrating out ${\sigma }^2$):
$\mathbf{\beta }\mid \mathbf{y}\sim t_{2a_n}\left({\mathbf{\beta }}_n,\frac{b_n}{a_n}{\mathbf{\Lambda }}_n^{-1}\right),$
a multivariate Student-$t$ centered at ${\mathbf{\beta }}_n$ with $2a_n$ degrees of freedom. The scale matrix $b_n{\mathbf{\Lambda }}_n^{-1}/a_n$ inherits both the design-matrix geometry (through ${\mathbf{\Lambda }}_n$) and the residual scale ${\hat{\sigma }}^2=b_n/a_n$.

Posterior predictive for new inputs ${\mathbf{X}}_{\star }$: again multivariate Student-$t$,
${\mathbf{y}}_{\star }\mid \mathbf{y}\sim t_{2a_n}\left({\mathbf{X}}_{\star }{\mathbf{\beta }}_n,\frac{b_n}{a_n}(\mathbf{I}+{\mathbf{X}}_{\star }{\mathbf{\Lambda }}_n^{-1}{\mathbf{X}}_{\star }^{\top })\right).$
The predictive covariance splits into noise variance ($\mathbf{I}$) and parameter-uncertainty variance (${\mathbf{X}}_{\star }{\mathbf{\Lambda }}_n^{-1}{\mathbf{X}}_{\star }^{\top }$) — the latter grows with distance from the training design, producing natural uncertainty inflation for extrapolation.

Ridge-regression correspondence. With prior $\mathbf{\beta }\sim N(0,{\sigma }^2{\lambda }^{-1}\mathbf{I})$, the posterior mean is
${\mathbf{\beta }}_n=({\mathbf{X}}^{\top }\mathbf{X}+\lambda \mathbf{I}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y},$
exactly the ridge regression estimator with penalty $\lambda$. Ridge is the Bayesian MAP under a Gaussian prior; the Bayesian interpretation adds a principled way to choose $\lambda$ (via marginal likelihood) and a credible interval on each coefficient.

Reference prior limit. With ${\mathbf{\Lambda }}_0\to 0$ and $a_0,b_0\to 0$ (improper prior $p(\mathbf{\beta },{\sigma }^2)\propto 1/{\sigma }^2$), the posterior mean equals the OLS estimate ${\hat{\mathbf{\beta }}}_{\text{OLS}}=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}$, and the marginal posterior is $\mathbf{\beta }\mid \mathbf{y}\sim t_{n-d}({\hat{\mathbf{\beta }}}_{\text{OLS}},{\hat{\sigma }}^2({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1})$ with ${\hat{\sigma }}^2=\mathrm{RSS}/(n-d)$ — numerically identical to the frequentist sampling distribution.

Prerequisites

• 3.1 Bayesian Framework — Bayesian mechanics.
• 3.5.2 Normal-Normal Model — the univariate NIG conjugate analysis this generalizes.
• 3.8 Bayesian Regression — parent index.
• 1.2. Linear Regression — frequentist linear regression for comparison.

Intuition

The conjugate analysis works because everything in the model is Gaussian. Writing the log-joint as a function of $\mathbf{\beta }$, both the prior and the likelihood contribute quadratic forms; their sum is another quadratic, completed to identify a Gaussian posterior for $\mathbf{\beta }$ conditional on ${\sigma }^2$. The variance, ${\sigma }^2$, has its own Inverse-Gamma conjugate prior, and the two pieces combine into the NIG hierarchy. Integrating ${\sigma }^2$ out of the joint posterior yields a Student-$t$ marginal for $\mathbf{\beta }$, mirroring the univariate Normal-Inverse-Gamma story in 3.5.2 Normal-Normal Model.

The matrix ${\mathbf{\Lambda }}_n={\mathbf{X}}^{\top }\mathbf{X}+{\mathbf{\Lambda }}_0$ is the posterior precision: the data contribute their Gram matrix ${\mathbf{X}}^{\top }\mathbf{X}$ and the prior contributes ${\mathbf{\Lambda }}_0$. For ridge, ${\mathbf{\Lambda }}_0=\lambda \mathbf{I}$ regularizes directions where the data are information-poor (small eigenvalues of ${\mathbf{X}}^{\top }\mathbf{X}$). The prior effectively adds phantom observations that pull $\mathbf{\beta }$ toward ${\mathbf{\beta }}_0$, and the ratio of prior-precision eigenvalues to data-precision eigenvalues controls how much each coordinate of $\mathbf{\beta }$ is shrunk.

The predictive variance is the crucial Bayesian payoff. Plug-in prediction ${\mathbf{X}}_{\star }\hat{\mathbf{\beta }}$ forgets that $\hat{\mathbf{\beta }}$ itself was estimated; the Bayesian predictive distribution correctly inflates variance by ${\mathbf{X}}_{\star }{\mathbf{\Lambda }}_n^{-1}{\mathbf{X}}_{\star }^{\top }$, which grows for test inputs far from the training design. A regression model is therefore automatically more uncertain about extrapolation — no post-hoc confidence-inflation tricks needed.

Properties

(Ridge as MAP Under Gaussian Prior) The ridge estimator is the posterior mode under $\mathbf{\beta }\sim N(0,{\sigma }^2/\lambda )$.

Setting ${\mathbf{\beta }}_0=0$ and ${\mathbf{\Lambda }}_0=\lambda \mathbf{I}$ yields posterior mean ${\mathbf{\beta }}_n=({\mathbf{X}}^{\top }\mathbf{X}+\lambda \mathbf{I}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}$, the ridge estimator. The Bayesian view explains why ridge shrinks: the Gaussian prior puts mass near zero, so coefficients without strong data support are pulled toward zero. The choice of $\lambda$ corresponds to choosing the prior’s concentration; cross-validation and Type-II MLE (evidence maximization, see 3.7 Marginal Likelihood) give two principled ways to set it.

(Lasso as MAP Under Laplace Prior) A Laplace (double-exponential) prior on $\mathbf{\beta }$ yields the LASSO posterior mode.

Taking ${\mathbf{\beta }}_j\stackrel{\text{iid}}{\sim }\text{Laplace}(0,1/\lambda )$ with density $\propto \exp (-\lambda |{\beta }_j|)$, the log-posterior is $\mathbf{y}$-fit plus $-\lambda \Vert \mathbf{\beta }{\Vert }_1$; maximizing it recovers the LASSO. The Laplace prior is not conjugate to the Gaussian likelihood — the posterior has no closed form, but it admits a data-augmented Gibbs sampler via a scale mixture representation. Practical Bayesian LASSO fits typically use the Park-Casella Gibbs sampler.

(Multivariate Student-t Marginal for $\mathbf{\beta }$) Integrating ${\sigma }^2$ out of the NIG posterior yields a Student-$t$ marginal for the coefficient vector.

$\mathbf{\beta }\mid \mathbf{y}\sim t_{2a_n}\left({\mathbf{\beta }}_n,\frac{b_n}{a_n}{\mathbf{\Lambda }}_n^{-1}\right).$
Degrees of freedom $2a_n$ grow linearly with $n$; the tails thicken for small $n$ and approach Gaussian as $n\to \infty$. The $j$-th marginal is univariate Student-$t$:
${\beta }_j\mid \mathbf{y}\sim t_{2a_n}\left({\beta }_{n,j},\frac{b_n}{a_n}[{\mathbf{\Lambda }}_n^{-1}{]}_{jj}\right),$
whose $(1-\alpha )$ credible interval is ${\beta }_{n,j}\pm t_{2a_n,1-\alpha /2}\sqrt{b_n[{\mathbf{\Lambda }}_n^{-1}{]}_{jj}/a_n}$. Under the reference prior limit this collapses to the classical $t$-interval.

(Posterior Predictive is Multivariate Student-t) Predictions at new inputs inherit heavy tails and design-dependent uncertainty.

${\mathbf{y}}_{\star }\mid \mathbf{y}\sim t_{2a_n}\left({\mathbf{X}}_{\star }{\mathbf{\beta }}_n,\frac{b_n}{a_n}(\mathbf{I}+{\mathbf{X}}_{\star }{\mathbf{\Lambda }}_n^{-1}{\mathbf{X}}_{\star }^{\top })\right).$
The predictive variance decomposes as noise variance $\mathbf{I}$ plus parameter-uncertainty variance ${\mathbf{X}}_{\star }{\mathbf{\Lambda }}_n^{-1}{\mathbf{X}}_{\star }^{\top }$. For a test input far from training, ${\mathbf{X}}_{\star }{\mathbf{\Lambda }}_n^{-1}{\mathbf{X}}_{\star }^{\top }$ is large; predictive intervals widen automatically. This is the discrete-grid analog of the increasing posterior variance in GP regression away from data.

(Reference-Prior Posterior Equals OLS Sampling Distribution) The Jeffreys prior limit recovers frequentist inference exactly.

With ${\mathbf{\Lambda }}_0=0$ and $a_0,b_0\to 0$,
${\mathbf{\beta }}_n=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}={\hat{\mathbf{\beta }}}_{\text{OLS}},\mathbf{\beta }\mid \mathbf{y}\sim t_{n-d}({\hat{\mathbf{\beta }}}_{\text{OLS}},{\hat{\sigma }}^2({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}),$
with ${\hat{\sigma }}^2=\mathrm{RSS}/(n-d)$. The Bayesian $(1-\alpha )$ credible interval coincides numerically with the frequentist $(1-\alpha )$ confidence interval, though the interpretations differ (see 3.2 Credible Intervals).

(Conjugate Gibbs Sampler) In non-conjugate extensions, the conditional for $\mathbf{\beta }$ given ${\sigma }^2$ stays Gaussian, enabling Gibbs.

If ${\sigma }^2$ has any prior (not necessarily Inverse-Gamma), the conditional $\mathbf{\beta }\mid {\sigma }^2,\mathbf{y}\sim N({\mathbf{\beta }}_n,{\sigma }^2{\mathbf{\Lambda }}_n^{-1})$ remains Gaussian — a full conditional ready for Gibbs. This conditional conjugacy is why Bayesian regression is a canonical Gibbs target even when the full joint is not conjugate. See 3.11.2 Gibbs Sampling.

(Evidence and Model Comparison) The marginal likelihood has closed form, enabling Bayesian model selection.

$p(\mathbf{y})=\frac{|{\mathbf{\Lambda }}_0{|}^{1/2}}{|{\mathbf{\Lambda }}_n{|}^{1/2}}\cdot \frac{b_0^{a_0}}{b_n^{a_n}}\cdot \frac{\Gamma (a_n)}{\Gamma (a_0)}\cdot (2\pi {)}^{-n/2}.$
This lets Bayes factors compare nested linear models (whether to include a covariate, polynomial order, etc.) without resorting to BIC. The closed form also enables Type-II MLE tuning of ${\mathbf{\Lambda }}_0$ (evidence-based choice of ridge penalty). See 3.7 Marginal Likelihood.

(Sensitivity to Prior Scale) Regression coefficient priors must be calibrated to covariate scale.

A prior ${\beta }_j\sim N(0,1)$ has very different regularizing force when the covariate $x_{ij}$ ranges over $[-1,1]$ than when it ranges over $[-1000,1000]$. Standard practice is to standardize covariates before regression, making a unit-scale prior the default. Alternately, use a scale-invariant hierarchical prior that learns an appropriate per-coefficient scale.

Examples

Example 1: Closed-Form Posterior for a Single-Predictor Regression

Data $\{(x_i,y_i){\}}_{i=1}^n$ with $x_i$ and $y_i$ centered (so intercept drops out) and $y_i=\beta x_i+{\epsilon }_i$, ${\epsilon }_i\sim N(0,{\sigma }^2)$. Prior $\beta \sim N(0,{\sigma }^2/\lambda )$ and improper $p({\sigma }^2)\propto 1/{\sigma }^2$.

Design quantities: $S_{xx}=\sum x_i^2$, $S_{xy}=\sum x_i{y}_i$, $S_{yy}=\sum y_i^2$.
Posterior for $\beta \mid {\sigma }^2$: $N({\beta }_n,{\sigma }^2/(S_{xx}+\lambda ))$ with ${\beta }_n=S_{xy}/(S_{xx}+\lambda )$ — the ridge estimator.
Posterior for ${\sigma }^2$: $\text{InvGamma}(n/2,(S_{yy}-{\beta }_n{S}_{xy})/2+\lambda {\beta }_n^2/2)$ under the improper prior.
Marginal posterior for $\beta$: $t_n({\beta }_n,{\hat{\sigma }}^2/(S_{xx}+\lambda ))$ with ${\hat{\sigma }}^2=b_n/a_n$.

Setting $\lambda =0$ recovers the frequentist OLS. Increasing $\lambda$ shrinks ${\beta }_n$ toward zero and narrows the credible interval around the shrunk value — a variance-bias tradeoff made explicit.

Example 2: Predictive Uncertainty for Extrapolation

Suppose a simple linear regression is fit with $n=50$ training points $x_i\in [-1,1]$. Consider predicting at $x_{\star }=2$ (outside the training range).

Plug-in variance: ${\hat{\sigma }}^2=$ residual variance, giving predictive standard deviation $\sqrt{{\hat{\sigma }}^2}$ regardless of $x_{\star }$ — wrong.
Bayesian predictive variance: ${\hat{\sigma }}^2(1+{\mathbf{x}}_{\star }^{\top }{\mathbf{\Lambda }}_n^{-1}{\mathbf{x}}_{\star })$. For centered $x$ and $\lambda =0$:
$\text{Var}[{\tilde{y}}_{\star }\mid \mathbf{y}]={\hat{\sigma }}^2\left(1+\frac{1}{n}+\frac{x_{\star }^2}{\sum x_i^2}\right).$
With $n=50$, $\sum x_i^2=16$ (say), and $x_{\star }=2$, the variance inflates by a factor $1+1/50+4/16=1.27$ over the plug-in. At $x_{\star }=10$, the factor is $1+0.02+6.25=7.27$ — a 2.7-fold increase in predictive standard deviation. Extrapolation uncertainty is automatic.

Example 3: Ridge via Marginal-Likelihood Tuning

With prior $\mathbf{\beta }\sim N(0,{\sigma }^2{\lambda }^{-1}\mathbf{I})$ and ${\sigma }^2\sim \text{InvGamma}(a_0,b_0)$, the marginal likelihood $p(\mathbf{y}\mid \lambda )$ is a smooth function of $\lambda$. Maximizing it gives an empirical-Bayes ridge penalty.

Recipe: for each candidate $\lambda$ on a log-grid, compute the closed-form evidence using the formula above and pick the maximizer. This works out to a GCV-like criterion but is principled (the evidence already penalizes complexity — see 3.7 Marginal Likelihood’s BIC property). For a design with $n=100,d=10$ and well-conditioned ${\mathbf{X}}^{\top }\mathbf{X}$, the evidence-optimal $\hat{\lambda }$ typically lands between cross-validation’s choice and the LOO-CV choice, with similar predictive performance.

Derivation

Goal: derive the NIG posterior for $(\mathbf{\beta },{\sigma }^2)$.

The likelihood times prior is
$p(\mathbf{y}\mid \mathbf{\beta },{\sigma }^2)p(\mathbf{\beta }\mid {\sigma }^2)p({\sigma }^2)\propto ({\sigma }^2{)}^{-n/2}\exp \left(-\frac{1}{2{\sigma }^2}\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }^2\right)$
$\times ({\sigma }^2{)}^{-d/2}|{\mathbf{\Lambda }}_0{|}^{1/2}\exp \left(-\frac{1}{2{\sigma }^2}(\mathbf{\beta }-{\mathbf{\beta }}_0{)}^{\top }{\mathbf{\Lambda }}_0(\mathbf{\beta }-{\mathbf{\beta }}_0)\right)\times ({\sigma }^2{)}^{-a_0-1}\exp (-b_0/{\sigma }^2).$
Collecting terms quadratic in $\mathbf{\beta }$ inside the exponent:
$\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }^2+(\mathbf{\beta }-{\mathbf{\beta }}_0{)}^{\top }{\mathbf{\Lambda }}_0(\mathbf{\beta }-{\mathbf{\beta }}_0)=(\mathbf{\beta }-{\mathbf{\beta }}_n{)}^{\top }{\mathbf{\Lambda }}_n(\mathbf{\beta }-{\mathbf{\beta }}_n)+\text{const}(\mathbf{\beta }),$
by completing the square, with ${\mathbf{\Lambda }}_n={\mathbf{X}}^{\top }\mathbf{X}+{\mathbf{\Lambda }}_0$ and ${\mathbf{\beta }}_n={\mathbf{\Lambda }}_n^{-1}({\mathbf{X}}^{\top }\mathbf{y}+{\mathbf{\Lambda }}_0{\mathbf{\beta }}_0)$. The constant-in-$\mathbf{\beta }$ piece is
$\text{const}(\mathbf{\beta })={\mathbf{y}}^{\top }\mathbf{y}+{\mathbf{\beta }}_0^{\top }{\mathbf{\Lambda }}_0{\mathbf{\beta }}_0-{\mathbf{\beta }}_n^{\top }{\mathbf{\Lambda }}_n{\mathbf{\beta }}_n.$
Substituting:
$p(\mathbf{\beta },{\sigma }^2\mid \mathbf{y})\propto ({\sigma }^2{)}^{-(n+d)/2-a_0-1}\exp \left(-\frac{1}{2{\sigma }^2}(\mathbf{\beta }-{\mathbf{\beta }}_n{)}^{\top }{\mathbf{\Lambda }}_n(\mathbf{\beta }-{\mathbf{\beta }}_n)\right)\exp \left(-\frac{b_n}{{\sigma }^2}\right),$
which is NIG$({\mathbf{\beta }}_n,{\mathbf{\Lambda }}_n,a_n,b_n)$ with $a_n,b_n$ as stated.

Marginal posterior for $\mathbf{\beta }$: integrating out ${\sigma }^2$ against the Inverse-Gamma marginal of the NIG. The resulting density is $\propto {\left[1+(\mathbf{\beta }-{\mathbf{\beta }}_n{)}^{\top }(b_n{\mathbf{\Lambda }}_n^{-1}/a_n{)}^{-1}(\mathbf{\beta }-{\mathbf{\beta }}_n)/(2a_n)\right]}^{-(2a_n+d)/2}$, the multivariate-$t$ kernel.

Posterior predictive: ${\mathbf{y}}_{\star }={\mathbf{X}}_{\star }\mathbf{\beta }+{\mathbf{\epsilon }}_{\star }$ with $\mathbf{\beta }\mid \mathbf{y}\sim t_{2a_n}({\mathbf{\beta }}_n,\cdot )$ and ${\mathbf{\epsilon }}_{\star }\mid {\sigma }^2\sim N(0,{\sigma }^2\mathbf{I})$. Summing two Student-$t$ dependence structures preserves degrees of freedom and adds scales.

See Also

• 3.8 Bayesian Regression — parent index.
• 3.5.2 Normal-Normal Model — univariate conjugate analysis this generalizes.
• 3.9 Laplace Approximation — non-conjugate regressions require this.
• 3.11.2 Gibbs Sampling — conditional conjugacy makes Gibbs the natural sampler for non-NIG priors (Laplace, horseshoe, etc.).
• 3.8.2 Bayesian AutoRegression — linear regression with lagged responses.
• 3.8.3 Bayesian Logistic Regression — non-conjugate binary-response counterpart.
• 3.8.4 Bayesian Poisson Regression — non-conjugate count-response counterpart.
• 3.10 Gaussian Processes — generalization to infinite-dimensional linear models via kernels.
• 1.2. Linear Regression — frequentist analog.