Bayesian Linear Regression

In classical linear regression, we learn a single set of weights w that best fits our training data. We find the point estimate that minimizes squared error, and we're done. But this approach hides a profound limitation: it tells us nothing about our uncertainty.

Consider fitting a line through three data points versus fitting the same line through three thousand points. Classical regression gives you the same type of answer in both cases—a single weight vector. But intuitively, you should be far more confident in the second case. Where does that confidence live in the classical framework? It doesn't.

Bayesian linear regression fundamentally reimagines this problem. Instead of asking "What are the best weights?", we ask "What is the probability distribution over all possible weights, given our data?" This shift from point estimates to probability distributions unlocks principled uncertainty quantification, robustness to limited data, and a framework that gracefully incorporates prior knowledge.

The core philosophical departure of Bayesian statistics is treating parameters as random variables rather than fixed but unknown constants. This isn't merely a mathematical trick—it represents a fundamentally different interpretation of what parameters mean.

The Frequentist View (Classical Linear Regression):

In the frequentist paradigm, the true weight vector w* exists as a fixed, deterministic quantity. Our training data is a random sample from some population, and our estimate ŵ is a random variable that approximates w*. The randomness comes from the data, not from the parameters.

This view leads to questions like:

• "What is the sampling distribution of our estimator?"
• "What is the probability that our estimate is within ε of the true value?"

The Bayesian View:

In the Bayesian paradigm, w itself is a random variable with a probability distribution. Before seeing data, we have a prior distribution p(w) representing our beliefs about likely weight values. After seeing data D, we update to a posterior distribution p(w|D) that combines prior beliefs with evidence.

This view leads to different questions:

• "What is the probability that w lies in a certain region, given what we've observed?"
• "How should we update our beliefs as we see more data?"

Before we can define a prior on weights, we need a clear probabilistic model for linear regression. Let's establish the framework carefully.

The Data:

We observe N training examples, each consisting of:

• A feature vector xₙ ∈ ℝᴰ (D features)
• A target value yₙ ∈ ℝ

We collect all features into a design matrix X ∈ ℝᴺˣᴰ and all targets into a vector y ∈ ℝᴺ.

The Generative Model:

We assume that targets are generated by a linear function of features plus Gaussian noise:

$$y_n = \mathbf{w}^\top \mathbf{x}_n + \epsilon_n$$

where:

• w ∈ ℝᴰ is the weight vector we want to learn
• εₙ ~ 𝒩(0, σ²) is independent Gaussian noise with variance σ²

This can be written equivalently as:

$$y_n | \mathbf{x}_n, \mathbf{w}, \sigma^2 \sim \mathcal{N}(\mathbf{w}^\top \mathbf{x}_n, \sigma^2)$$

The Likelihood Function:

Under the assumption that observations are conditionally independent given w, the likelihood of observing all data is:

$$p(\mathbf{y} | \mathbf{X}, \mathbf{w}, \sigma^2) = \prod_{n=1}^{N} \mathcal{N}(y_n | \mathbf{w}^\top \mathbf{x}_n, \sigma^2)$$

Taking the log and expanding:

$$\log p(\mathbf{y} | \mathbf{X}, \mathbf{w}, \sigma^2) = -\frac{N}{2} \log(2\pi\sigma^2) - \frac{1}{2\sigma^2} \sum_{n=1}^{N} (y_n - \mathbf{w}^\top \mathbf{x}_n)^2$$

Maximizing this log-likelihood with respect to w is equivalent to minimizing the sum of squared errors—exactly what ordinary least squares (OLS) does. So OLS is the maximum likelihood estimator for this probabilistic model.

Matrix Notation:

For computational convenience, we write the model in matrix form:

$$\mathbf{y} = \mathbf{X}\mathbf{w} + \boldsymbol{\epsilon}$$

where ε ~ 𝒩(0, σ²Iₙ), giving us:

$$\mathbf{y} | \mathbf{X}, \mathbf{w}, \sigma^2 \sim \mathcal{N}(\mathbf{X}\mathbf{w}, \sigma^2 \mathbf{I}_N)$$

This multivariate Gaussian formulation is crucial for deriving the posterior distribution efficiently.

The prior distribution p(w) encodes our beliefs about the weights before observing any data. This is where we inject domain knowledge, regularization preferences, or express uncertainty.

The Gaussian Prior:

The most common choice is a multivariate Gaussian prior:

$$\mathbf{w} \sim \mathcal{N}(\mathbf{m}_0, \mathbf{S}_0)$$

where:

• m₀ ∈ ℝᴰ is the prior mean (our initial guess for the weights)
• S₀ ∈ ℝᴰˣᴰ is the prior covariance matrix (encoding uncertainty and correlations)

Why Gaussian?

1. Conjugacy: Gaussian priors are conjugate to the Gaussian likelihood, meaning the posterior is also Gaussian. This allows closed-form solutions.

2. Interpretability: The prior mean represents our initial estimate; the prior covariance represents how uncertain we are and how features relate.

3. Flexibility: By choosing m₀ and S₀ appropriately, we can encode a wide range of prior beliefs.

4. Central Limit Theorem: In many settings, the Gaussian emerges as a natural limiting distribution.

The Isotropic Gaussian Prior (Most Common):

Often, we use a zero-mean isotropic prior:

$$\mathbf{w} \sim \mathcal{N}(\mathbf{0}, \alpha^{-1}\mathbf{I}_D)$$

where:

• m₀ = 0: No initial bias toward any particular weights
• S₀ = α⁻¹Iᴰ: Equal, independent uncertainty for all weights
• α is the precision (inverse variance) of the prior

This prior says: "Before seeing data, I believe each weight is equally likely to be positive or negative, with magnitude controlled by α."

The prior precision α (or equivalently, the prior variance α⁻¹) is a critical hyperparameter that controls the strength of our prior beliefs.

High Precision (Large α):

• Prior variance α⁻¹ is small
• Prior concentrated tightly around zero
• Strong belief that weights should be small
• Data must be very convincing to move weights away from zero
• Acts as strong regularization

Low Precision (Small α):

• Prior variance α⁻¹ is large
• Prior spread broadly
• Weak belief about weight magnitudes
• Even weak data signal can dominate
• Acts as weak regularization

The Limit α → 0 (Improper Flat Prior):

As α → 0, the prior variance → ∞, and we approach an "uninformative" or flat prior. In this limit, the Bayesian posterior mode matches the OLS solution. However, truly flat priors are technically improper (don't integrate to 1) and can cause issues.

The Limit α → ∞ (Infinitely Informative Prior):

As α → ∞, the prior collapses to a point mass at zero. No matter what data we see, the weights stay at zero. The prior completely dominates.

The Regularization Interpretation:

We'll explore this deeply in a later page, but the key insight is:

$$\text{Posterior Mode} = \arg\max_\mathbf{w} \left[ \log p(\mathbf{y}|\mathbf{X}, \mathbf{w}) + \log p(\mathbf{w}) \right]$$

For a Gaussian prior with precision α and Gaussian likelihood with noise variance σ²:

$$\text{Posterior Mode} = \arg\min_\mathbf{w} \left[ |\mathbf{y} - \mathbf{X}\mathbf{w}|^2 + \frac{\alpha\sigma^2}{1} |\mathbf{w}|^2 \right]$$

This is exactly Ridge Regression with regularization parameter λ = ασ². The Bayesian prior directly corresponds to the regularization term!

This connection reveals that every regularized method implicitly assumes some prior. Bayesian inference makes this assumption explicit and principled.

While the isotropic zero-mean prior is common, real problems often benefit from more sophisticated prior specifications.

Non-Zero Mean:

If we have prior knowledge that certain features should have positive or negative effects, we can encode this:

$$\mathbf{w} \sim \mathcal{N}(\mathbf{m}_0, \mathbf{S}_0)$$

with m₀ ≠ 0. For example, in a housing price model, we might believe:

• Prior mean for "square footage" weight: positive (larger homes cost more)
• Prior mean for "distance to city center" weight: negative (farther is cheaper)

Diagonal (Anisotropic) Prior:

Different features may warrant different prior variances:

$$\mathbf{S}_0 = \text{diag}(\sigma_1^2, \sigma_2^2, \ldots, \sigma_D^2)$$

This is useful when:

• Some features are known to be more reliable predictors (smaller variance)
• Some features are exploratory and should have high uncertainty
• Feature scales differ substantially

Full Covariance Prior:

For correlated features, a full covariance matrix captures prior beliefs about relationships:

$$\mathbf{S}_0 = \begin{pmatrix} \sigma_1^2 & \rho\sigma_1\sigma_2 \ \rho\sigma_1\sigma_2 & \sigma_2^2 \end{pmatrix}$$

Positive correlation in S₀ says: "If weight 1 is large, weight 2 is likely large too." This can encode structural knowledge about feature relationships.

While Gaussian priors are convenient and interpretable, other prior families encode different assumptions about weight structure.

Laplace Prior (Sparse Weights):

$$p(w_j) \propto \exp(-\lambda |w_j|)$$

The Laplace (double exponential) prior has a sharp peak at zero and heavy tails. This encourages sparsity—many weights exactly or nearly zero, with a few large weights. The Laplace prior corresponds to L1 regularization (Lasso).

Spike-and-Slab Prior (Explicit Sparsity):

$$p(w_j) = \pi \cdot \delta(w_j) + (1-\pi) \cdot \mathcal{N}(w_j | 0, \sigma^2)$$

A mixture of a point mass at zero (spike) and a diffuse Gaussian (slab). This explicitly models the belief that each weight is either exactly zero or drawn from a continuous distribution. Powerful for feature selection but computationally challenging.

Horseshoe Prior (Adaptive Shrinkage):

$$w_j | \lambda_j \sim \mathcal{N}(0, \lambda_j^2), \quad \lambda_j \sim \text{Half-Cauchy}(0, \tau)$$

A hierarchical prior where each weight has its own scale parameter λⱼ with a heavy-tailed distribution. This provides adaptive shrinkage: truly zero weights are shrunk aggressively, while large true weights experience less shrinkage.

Automatic Relevance Determination (ARD):

$$w_j \sim \mathcal{N}(0, \alpha_j^{-1})$$

Each weight has its own precision αⱼ, learned from data. If αⱼ → ∞ during learning, the corresponding feature is deemed irrelevant. ARD provides automatic feature selection within a Gaussian framework.

One of the most profound insights from Bayesian linear regression is that every regularization scheme corresponds to some prior distribution. This isn't just a mathematical curiosity—it provides deep insight into what regularization actually means.

The MAP Estimator:

The Maximum A Posteriori (MAP) estimate is the mode of the posterior:

$$\mathbf{w}{\text{MAP}} = \arg\max\mathbf{w} p(\mathbf{w} | \mathbf{y}, \mathbf{X}) = \arg\max_\mathbf{w} \left[ p(\mathbf{y} | \mathbf{X}, \mathbf{w}) \cdot p(\mathbf{w}) \right]$$

Taking logs:

$$\mathbf{w}{\text{MAP}} = \arg\min\mathbf{w} \left[ -\log p(\mathbf{y} | \mathbf{X}, \mathbf{w}) - \log p(\mathbf{w}) \right]$$

The first term is the negative log-likelihood (data fit). The second term is the negative log-prior (regularization).

Specific Correspondences:

Why This Matters:

1. Principled Hyperparameter Interpretation: Regularization strength λ isn't arbitrary—it's the ratio of prior precision to noise precision. Setting λ = 1 means you trust the prior equally to one data point.

2. Informed Prior Design: Want L1-like sparsity? Use a Laplace prior. Want smooth solutions? Use a prior that penalizes weight differences.

3. Unified Framework: All regularization methods become instances of Bayesian inference with different priors. This unified view clarifies when each method is appropriate.

4. Beyond Point Estimates: Once you recognize the prior, you can do full Bayesian inference with it—getting posterior distributions, not just MAP estimates.

Setting priors in practice requires balancing domain knowledge, computational tractability, and robustness. Here are key principles:

1. Feature Scaling and Priors:

If features are on different scales, an isotropic prior may be inappropriate. A weight of 0.01 for a feature measured in millions differs from 0.01 for a binary feature. Options:

• Standardize features to zero mean, unit variance
• Use feature-specific prior variances (ARD)
• Design priors in the standardized space

2. Prior Predictive Checks:

Before seeing data, sample from the prior and simulate predictions: $$\mathbf{w} \sim p(\mathbf{w}), \quad \mathbf{y}_{\text{simulated}} = \mathbf{X}\mathbf{w} + \boldsymbol{\epsilon}$$

Are the simulated predictions plausible? If your prior produces predictions like "house prices of -$10 million," the prior is poorly calibrated.

3. Weakly Informative Priors:

If you lack strong prior knowledge, use weakly informative priors that:

• Exclude implausible regions (e.g., infinite weights)
• Don't commit strongly to any particular value
• Provide just enough regularization for stable inference

A common choice: 𝒩(0, s²) where s is set to cover plausible weight magnitudes.

4. Sensitivity Analysis:

Run inference with different prior settings. If conclusions change dramatically with minor prior adjustments, the data may be insufficient or the model misspecified.

We've established the foundational concepts for Bayesian linear regression by understanding how to place prior distributions on weights. Let's consolidate the key insights:

What's Next:

With the prior established, we're ready to see what we actually learn when data arrives. The next page derives the posterior distribution over weights—the result of combining our prior beliefs with the evidence from observed data. This posterior is the central object in Bayesian inference, from which all predictions and uncertainty estimates flow.