Bayesian Linear Regression

December 16, 2025

On this page

Bayesian Linear Regression: Fundamentals, Ridge, and Lasso Regularization

In linear regression, we treat the coefficients as unknown constants and fit them by maximizing the likelihood. Here, we combine the same likelihood with a prior distribution. Bayes' rule then gives the posterior distribution, which describes uncertainty about the coefficients after observing the data. We derive the Gaussian posterior and show how posterior modes under Gaussian and Laplace priors connect to ridge and lasso estimates.

Recap: The Gaussian likelihood

For the fixed-design model in linear regression, the part of the likelihood that depends on $\boldsymbol\beta$ is

p(\mathbf y\mid\boldsymbol\beta,\mathbf X,\sigma^2) \ \propto\ \exp\!\left[-\frac{\|\mathbf y-\mathbf X\boldsymbol\beta\|_2^2}{2\sigma^2}\right].

Here $\sigma^2>0$ and the observed $(\mathbf X,\mathbf y)$ are fixed. Multiplying by a prior adds its contribution to the log posterior; it does not change this sampling model.

Key Bayesian Regression Framework

Model Setup and Notation

We begin with the familiar linear model:

\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon},

where:

$\mathbf{y}$ is an $n \times 1$ response vector,
$\mathbf{X}$ is an $n \times (p+1)$ design matrix (including a column for the intercept),
$\boldsymbol{\beta}$ is a $(p+1) \times 1$ vector of unknown coefficients, and
the errors $\boldsymbol{\epsilon}$ are assumed to be independent and normally distributed as

\boldsymbol{\epsilon} \sim N(\mathbf{0}, \sigma^2 \mathbf{I}).

Throughout the coefficient derivations, we condition on a fixed noise variance $\sigma^2>0$ and fixed positive prior hyperparameters. We also apply the displayed priors to all $p+1$ coefficients, including the intercept, so that the matrix formulas have a consistent convention. An analysis that estimates the noise variance or gives the intercept a separate prior needs those additional choices stated explicitly.

The design $\mathbf X$ is also fixed. Before observing data, $\boldsymbol\beta$ has the specified prior and $\mathbf Y\mid\boldsymbol\beta,\mathbf X,\sigma^2\sim N(\mathbf X\boldsymbol\beta,\sigma^2\mathbf I)$ .

After observing $\mathbf Y=\mathbf y$ , we keep $\mathbf y$ fixed and average over $\boldsymbol\beta\sim p(\boldsymbol\beta\mid\mathbf y,\mathbf X,\sigma^2)$ to obtain posterior means and covariances. A different calculation, integrating the likelihood over the prior, gives the marginal likelihood. We omit $\mathbf X$ from the densities below, but never average over it.

Standard Bayesian Linear Regression

Likelihood

The likelihood records how the observed data $\mathbf{y}$ depend on the parameters $\boldsymbol{\beta}$ and $\sigma^2$ . Because the response is continuous, we use its probability density, viewed as a function of those parameters. Under the independent Gaussian error model, the likelihood is:

L(\boldsymbol{\beta}, \sigma^2) = p(\mathbf{y} \mid \boldsymbol{\beta}, \sigma^2) = \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left(-\frac{(y_i - \mathbf{x}_i^T \boldsymbol{\beta})^2}{2\sigma^2}\right),

where $\mathbf{x}_i^T$ is the $i$ th row of $\mathbf{X}$ . When the noise variance is held fixed, we may omit factors that do not depend on the coefficients and write:

p(\mathbf{y} \mid \boldsymbol{\beta}, \sigma^2) \propto \exp\!\Bigl(-\frac{1}{2\sigma^2}\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2\Bigr).

Prior

We describe uncertainty about $\boldsymbol{\beta}$ with a prior distribution. A Gaussian prior on $\boldsymbol{\beta}$ is conjugate to the Gaussian likelihood when the noise variance is held fixed: the conditional posterior remains Gaussian, making the update analytically tractable.

\boldsymbol{\beta}\mid\sigma^2 \sim N\Bigl(\mathbf{0},\, \sigma^2 \tau^{-1} \mathbf{I}\Bigr),

where $\tau$ controls the spread or uncertainty of our prior relative to the noise variance $\sigma^2$ .

Conjugacy makes the Gaussian prior convenient, but convenience is not the only consideration. Different priors express different assumptions about plausible coefficient values:

Flat Prior: A constant prior density over $\boldsymbol{\beta}$ , written as $p(\boldsymbol{\beta}) \propto 1$ , assigns the same density to every coefficient value. On an unbounded parameter space, this density cannot integrate to one, so the prior is improper. We must check that the resulting posterior can be normalized to a probability distribution.
Weakly Informative Gaussian Prior: A Gaussian prior with a large variance (small $\tau$ ) permits a wider range of coefficient values. Whether it is weakly informative depends on the scales of the predictors and response.
Laplace Prior (Double-Exponential): A Laplace prior concentrates density near zero and connects posterior mode estimation to lasso variable selection:

$p(\beta_j) = \frac{\lambda}{2} \exp\!\left(-\lambda |\beta_j|\right).$
Heavy-Tailed Priors: A heavy-tailed prior, such as a Student's t-distribution for $\boldsymbol{\beta}$ , allows larger coefficient values than a Gaussian prior with comparable central concentration. This concerns the coefficient distribution; it does not by itself make the observation model robust to outliers.

The prior determines which coefficient values receive more weight before observing the responses. Gaussian and Laplace priors impose different penalties on large coefficients. We show this connection by deriving the maximum a posteriori (MAP) estimate, the coefficient vector at which the posterior density is largest. The prior scales must be chosen in relation to the predictor and response scales; the distribution name alone does not determine a suitable prior.

Posterior: Combine Likelihood and the Prior

For this introductory statement of Bayes' rule, we suppress the fixed design matrix and noise variance from the notation:

p(\boldsymbol{\beta} \mid \mathbf{y}) \;=\; \frac{p(\mathbf{y}, \boldsymbol{\beta})}{p(\mathbf{y})}=\; \frac{p(\mathbf{y} \mid \boldsymbol{\beta}) \;\times\; p(\boldsymbol{\beta})}{p(\mathbf{y})}.

Next, we write the marginal likelihood $p(\mathbf{y})$ by integrating over all possible $\boldsymbol{\beta}$ :

p(\mathbf{y}) \;=\; \int p(\mathbf{y}, \boldsymbol{\beta}) \,\mathrm{d}\boldsymbol{\beta} \;=\; \int p(\mathbf{y} \mid \boldsymbol{\beta}) \;\times\; p(\boldsymbol{\beta}) \,\mathrm{d}\boldsymbol{\beta}.

The marginal likelihood $p(\mathbf{y})$ normalizes the posterior and does not vary with the coefficients. To derive the posterior's form, we can therefore work up to a proportionality constant. For $\boldsymbol{\beta}$ , Bayes' rule becomes:

p(\boldsymbol{\beta} \mid \mathbf{y}) \;\propto\; p(\mathbf{y} \mid \boldsymbol{\beta}) \;\times\; p(\boldsymbol{\beta}).

Therefore, for our setup the posterior distribution is:

p(\boldsymbol{\beta} \mid \mathbf{y}, \sigma^2) \propto p(\mathbf{y} \mid \boldsymbol{\beta}, \sigma^2) \times p(\boldsymbol{\beta}\mid\sigma^2).

With the noise variance held fixed, the Gaussian likelihood and prior give a Gaussian posterior for $\boldsymbol{\beta}$ . This conjugacy lets us describe coefficient uncertainty analytically, rather than approximating the posterior numerically.

Completing the square gives the conditional posterior:

\boldsymbol{\beta}\mid\mathbf{y},\sigma^2 \sim N\!\left(\boldsymbol{m},\,\sigma^2 A^{-1}\right), \qquad A=\mathbf{X}^T\mathbf{X}+\tau\mathbf{I}, \quad \boldsymbol{m}=A^{-1}\mathbf{X}^T\mathbf{y}.

The matrix $A$ is positive definite because $\tau>0$ . In particular, $\mathbb E[\boldsymbol\beta\mid\mathbf y,\mathbf X,\sigma^2]=\boldsymbol m$ and $\operatorname{Cov}(\boldsymbol\beta\mid\mathbf y,\mathbf X,\sigma^2)=\sigma^2A^{-1}$ ; the posterior mode also equals $\boldsymbol m$ . This is the fixed-variance case of the Gaussian regression posterior derivation.

Posterior Distributions: Analytic and Approximate Forms

The Gaussian calculation is an example of conjugacy: updating the prior with the likelihood keeps the posterior within the same distribution family. Conjugacy is useful because it gives an explicit posterior, but it is not available for every likelihood–prior combination. Logistic regression provides a useful contrast.

Non-Analytic Example: Bayesian Logistic Regression

Model Setup

Likelihood

Suppose $\mathbf{y}\in \{0,1\}^n$ are binary outcomes and $X$ is the design matrix. For logistic regression, each $y_i$ is modeled with

P\bigl(y_i = 1 \mid \mathbf{x}_i,\boldsymbol{\beta}\bigr) \;=\; \sigma\!\bigl(\mathbf{x}_i^\mathsf{T}\,\boldsymbol{\beta}\bigr),

where $\sigma(z) = \tfrac{1}{1+e^{-z}}$ is the logistic function. The likelihood for all $n$ data points is:

p(\mathbf{y} \mid \boldsymbol{\beta}) \;=\; \prod_{i=1}^n \sigma\!\bigl(\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta}\bigr)^{\,y_i} \Bigl(1 - \sigma\!\bigl(\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta}\bigr)\Bigr)^{\,1 - y_i}.

Prior

We can again choose a Gaussian prior for $\boldsymbol{\beta}$ :

\boldsymbol{\beta} \;\sim\; \mathcal{N}\bigl(\boldsymbol{\mu}_0,\;\Sigma_0\bigr).

Posterior Challenge

Form of the Posterior

Using Bayes’ rule, the posterior is:

p(\boldsymbol{\beta} \mid \mathbf{y}) \;\propto\; \underbrace{\prod_{i=1}^n \sigma\!\bigl(\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta}\bigr)^{\,y_i} \bigl[\,1 - \sigma\!\bigl(\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta}\bigr)\bigr]^{1 - y_i}}_{\text{Logistic likelihood}} \;\times\; \underbrace{\exp\!\Bigl(-\tfrac12\bigl(\boldsymbol{\beta}-\boldsymbol{\mu}_0\bigr)^\mathsf{T}\Sigma_0^{-1}\bigl(\boldsymbol{\beta}-\boldsymbol{\mu}_0\bigr)\Bigr)}_{\text{Gaussian prior}}.

Unlike the Gaussian likelihood, this log-likelihood is not quadratic in the coefficients: its logistic terms remain after multiplication by the Gaussian prior. Completing the square therefore does not produce a Gaussian posterior.

No Closed-Form Integral

Unlike the Gaussian–Gaussian case, the integral

p(\mathbf{y}) \;=\; \int p(\mathbf{y},\boldsymbol{\beta})\,d\boldsymbol{\beta} \;=\; \int p(\mathbf{y} \mid \boldsymbol{\beta})\,p(\boldsymbol{\beta})\,d\boldsymbol{\beta}

does not generally yield a closed-form normalizing constant. We therefore need numerical methods for posterior inference.

Markov chain Monte Carlo (MCMC) approximates posterior summaries through samples, whereas variational inference approximates the posterior within a chosen distribution family. Both extend Bayesian inference beyond conjugate models, but their numerical accuracy must be assessed.

Returning to linear regression, the posterior describes uncertainty about each $\beta_j$ under the chosen model and prior. We next examine a particular posterior summary—the mode—to make the connection with regularization explicit.

Bayesian Ridge Regression

In Bayesian Ridge Regression, we use the same Gaussian prior:

\boldsymbol{\beta}\mid\sigma^2 \sim N\Bigl(\mathbf{0},\, \sigma^2 \tau^{-1} \mathbf{I}\Bigr).

Likelihood:
The likelihood function remains:

p(\mathbf{y} \mid \boldsymbol{\beta}, \sigma^2) \propto \exp\!\Bigl(-\frac{1}{2\sigma^2}\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2\Bigr).

Prior:
The Gaussian prior on the coefficients is:

p(\boldsymbol{\beta}\mid\sigma^2) \propto \exp\!\Bigl(-\frac{\tau}{2\sigma^2}\boldsymbol{\beta}^T \boldsymbol{\beta}\Bigr).

Posterior:
By Bayes’ rule, the posterior is:

p(\boldsymbol{\beta} \mid \mathbf{y}, \sigma^2) \propto p(\mathbf{y} \mid \boldsymbol{\beta}, \sigma^2) \, p(\boldsymbol{\beta}\mid\sigma^2) \propto \exp\!\left[-\frac{\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 + \tau\|\boldsymbol{\beta}\|^2}{2\sigma^2}\right].

Taking the negative logarithm (and ignoring constants) gives the following objective for Maximum A Posteriori (MAP) estimation:

\min_{\boldsymbol{\beta}} \; \frac{1}{2\sigma^2}\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 + \frac{\tau}{2\sigma^2}\boldsymbol{\beta}^T \boldsymbol{\beta}.

Multiplying through by $2\sigma^2$ (which does not affect the minimizer) leads to:

\min_{\boldsymbol{\beta}} \; \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 + \tau \, \|\boldsymbol{\beta}\|^2.

Derivation of the MAP Estimate:
Taking the derivative with respect to $\boldsymbol{\beta}$ and setting it to zero:

-2\mathbf{X}^T(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) + 2\tau \, \boldsymbol{\beta} = 0.

Rearrange to obtain:

\mathbf{X}^T\mathbf{X}\boldsymbol{\beta} + \tau \, \boldsymbol{\beta} = \mathbf{X}^T\mathbf{y}.

Thus, the MAP estimate for Bayesian Ridge Regression is:

\hat{\boldsymbol{\beta}}_{\text{BayesRidge}} = \bigl(\mathbf{X}^T\mathbf{X} + \tau \, \mathbf{I}\bigr)^{-1}\mathbf{X}^T\mathbf{y}.

The penalty coefficient is $\tau$ because the prior covariance was scaled by $\sigma^2$ . With an unscaled prior covariance $\tau^{-1}\mathbf{I}$ instead, the equivalent RSS penalty would be $\sigma^2\tau$ . These are different prior conventions and should not be mixed within a derivation.

Bayesian Lasso Regression

Bayesian lasso uses a Laplace (double-exponential) prior on the regression coefficients. The connection to sparse estimates concerns the posterior mode; the full posterior still describes uncertainty about those coefficients. We use a fixed-variance, unscaled-prior formulation below. The fully hierarchical model of Park and Casella (2008) instead scales the Laplace prior with the noise standard deviation and also models the noise variance.

Model, Likelihood, and Prior for Bayesian Lasso

Model:
We retain the linear model:

\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}, \quad \text{with } \boldsymbol{\epsilon} \sim N(\mathbf{0}, \sigma^2 \mathbf{I}).

Likelihood:
The likelihood function is:

p(\mathbf{y} \mid \boldsymbol{\beta}, \sigma^2) \propto \exp\!\Bigl(-\frac{1}{2\sigma^2}\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2\Bigr).

Prior:
For this Laplace-prior model, we place independent priors on each coefficient $\beta_j$ , for $j=0,\ldots,p$ . Unlike the Gaussian prior above, the scale used here is not multiplied by the noise standard deviation:

p(\beta_j) = \frac{\lambda}{2} \exp\!\left(-\lambda |\beta_j|\right).

For the entire coefficient vector:

p(\boldsymbol{\beta}) = \prod_{j=0}^{p} \frac{\lambda}{2} \exp\!\left(-\lambda |\beta_j|\right) \propto \exp\!\Bigl(-\lambda \|\boldsymbol{\beta}\|_1\Bigr).

Posterior and MAP Estimation for Bayesian Lasso

Posterior:
Using Bayes’ rule, the posterior distribution is given by:

p(\boldsymbol{\beta} \mid \mathbf{y}, \sigma^2) \propto p(\mathbf{y} \mid \boldsymbol{\beta}, \sigma^2) \, p(\boldsymbol{\beta}) \propto \exp\!\Bigl(-\frac{1}{2\sigma^2}\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 - \lambda \|\boldsymbol{\beta}\|_1\Bigr).

Taking the negative logarithm (ignoring constants) yields the MAP objective:

\min_{\boldsymbol{\beta}} \; \frac{1}{2\sigma^2}\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 + \lambda \|\boldsymbol{\beta}\|_1.

Multiplying through by $2\sigma^2$ for clarity, we have:

\min_{\boldsymbol{\beta}} \; \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 + 2\sigma^2 \lambda \, \|\boldsymbol{\beta}\|_1.

Note:
The $L_1$ norm is not differentiable at zero, so ordinary gradient equations do not apply there. In general, we compute the MAP estimate of $\boldsymbol{\beta}$ numerically, for example by coordinate descent. Sampling methods such as Gibbs sampling answer a different question: they approximate the full posterior rather than only its mode.

Hyperparameter Interpretation

The Gaussian prior covariance decreases as $\tau$ increases, so a larger $\tau$ penalizes departures from zero more strongly. For the Laplace prior, increasing $\lambda$ strengthens the corresponding $L_1$ penalty and can produce a sparser MAP estimate. Both choices also change the full posterior, not just the location of its mode.

Summary of Bayesian Linear Regression Methods

Standard Bayesian Linear Regression:
- Model:
  $\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}, \quad \boldsymbol{\epsilon} \sim N(\mathbf{0}, \sigma^2 \mathbf{I}).$
- Notation:
  $\mathbf{y}$ is an $n \times 1$ response vector, $\mathbf{X}$ is an $n \times (p+1)$ design matrix, $\boldsymbol{\beta}$ is a $(p+1) \times 1$ coefficient vector.
- Prior:
  $\boldsymbol{\beta}\mid\sigma^2 \sim N\Bigl(\mathbf{0},\, \sigma^2 \tau^{-1} \mathbf{I}\Bigr)$ .
- Posterior:
  $p(\boldsymbol{\beta} \mid \mathbf{y}, \sigma^2) \propto \exp\!\Bigl(-\frac{1}{2\sigma^2}\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 - \frac{\tau}{2\sigma^2}\boldsymbol{\beta}^T \boldsymbol{\beta}\Bigr)$ , a Gaussian distribution with the mean and covariance derived above.
Bayesian Ridge Regression:
- Prior:
  $\boldsymbol{\beta}\mid\sigma^2 \sim N\Bigl(\mathbf{0},\, \sigma^2 \tau^{-1} \mathbf{I}\Bigr)$ .
- MAP Estimator:
  $\hat{\boldsymbol{\beta}}_{\text{BayesRidge}} = \bigl(\mathbf{X}^T\mathbf{X} + \tau \, \mathbf{I}\bigr)^{-1}\mathbf{X}^T\mathbf{y}.$
- Benefits:
  Provides shrinkage of coefficients and quantifies uncertainty via the full posterior distribution.
Bayesian Lasso Regression:
- Prior:
  Laplace prior $p(\boldsymbol{\beta}) \propto \exp\!\Bigl(-\lambda \|\boldsymbol{\beta}\|_1\Bigr)$ .
- MAP Objective:
  $\min_{\boldsymbol{\beta}} \; \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 + 2\sigma^2 \lambda \, \|\boldsymbol{\beta}\|_1.$
- Benefits:
  Can produce a sparse MAP estimate for variable selection. The full posterior and its mode serve different purposes; numerical methods are generally used to compute the required summary.

Practical Considerations

The computation depends on the model and the posterior summary we need. Conjugate Gaussian models allow analytic calculations, while more complex models may require MCMC methods such as Gibbs sampling or Metropolis–Hastings. For posterior simulation, we check whether the chain has stabilized and explored the relevant parts of the posterior. For variational inference, we assess what the chosen family of distributions can represent. Both approaches still require us to assess the model and prior.

The distinction between a posterior and a point estimate is central to this note. Gaussian and Laplace priors connect Bayesian analysis to familiar regularization objectives, but the posterior contains more information than its mode. Keeping those roles separate helps us choose both the prior and the summary that answer the scientific question.