We've derived the posterior distribution over weights—a complete description of our uncertainty about the model parameters. But in practice, we don't care about weights directly. What we care about is prediction: given a new input x₊, what is the corresponding output y₊?

Classical regression gives a point prediction: ŷ₊ = wᵀx₊. But this ignores all uncertainty. If we're uncertain about the weights, shouldn't we be uncertain about predictions made with those weights?

Bayesian linear regression answers with the predictive distribution: a full probability distribution over the target value for a new input. This distribution captures two sources of uncertainty:

1. Epistemic uncertainty: Uncertainty about the weights (we don't know the true w)
2. Aleatoric uncertainty: Inherent noise in the observation process (even with known w, observations are noisy)

The predictive distribution integrates out weight uncertainty, giving us principled confidence intervals that reflect both what we don't know and what is inherently unknowable.

Given a new input x₊ ∈ ℝᴰ, we want to predict the corresponding target value y₊. In the Bayesian framework, we compute the predictive distribution:

$$p(y_* | \mathbf{x}_*, \mathbf{y}, \mathbf{X})$$

This is the probability distribution over y₊ given:

• The new input x₊
• All observed training data (X, y)

The Key Insight: Marginalization

We don't know the true weights w. But we have a posterior distribution over them. The predictive distribution averages over all possible weights, weighted by their posterior probability:

$$p(y_* | \mathbf{x}*, \mathbf{y}, \mathbf{X}) = \int p(y* | \mathbf{x}_*, \mathbf{w}) \cdot p(\mathbf{w} | \mathbf{y}, \mathbf{X}) , d\mathbf{w}$$

This integral marginalizes out the weights. The result accounts for all sources of uncertainty:

• Each possible w contributes a prediction
• Weights are weighted by how probable they are given the data
• The spread reflects both weight uncertainty and observation noise

Let's compute the integral. We have:

Likelihood for new observation: $$p(y_* | \mathbf{x}*, \mathbf{w}) = \mathcal{N}(y* | \mathbf{w}^\top\mathbf{x}_*, \sigma^2)$$

Posterior over weights: $$p(\mathbf{w} | \mathbf{y}, \mathbf{X}) = \mathcal{N}(\mathbf{w} | \mathbf{m}_N, \mathbf{S}_N)$$

The predictive distribution is: $$p(y_* | \mathbf{x}*, \mathbf{y}, \mathbf{X}) = \int \mathcal{N}(y* | \mathbf{w}^\top\mathbf{x}_*, \sigma^2) \cdot \mathcal{N}(\mathbf{w} | \mathbf{m}_N, \mathbf{S}_N) , d\mathbf{w}$$

This is the convolution of two Gaussians—which is itself Gaussian! The result follows from the general formula for marginalizing Gaussian random variables.

Step 1: Linear Transformation of Weights

Define the latent prediction (without noise): f₊ = wᵀx₊

Since w ~ 𝒩(mₙ, Sₙ), and f₊ is a linear function of w:

$$f_* = \mathbf{w}^\top\mathbf{x}* \sim \mathcal{N}(\mathbf{m}N^\top\mathbf{x}*, \mathbf{x}^\top\mathbf{S}N\mathbf{x})$$

The mean is mₙᵀx₊ and variance is x₊ᵀSₙx₊.

Step 2: Add Observation Noise

The observed y₊ is f₊ plus independent Gaussian noise: $$y_* = f_* + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2)$$

Since f₊ and ε are independent Gaussians, their sum is Gaussian with:

• Mean = mean of f₊ = mₙᵀx₊
• Variance = var(f₊) + var(ε) = x₊ᵀSₙx₊ + σ²

The predictive variance σ²₊ = x₊ᵀSₙx₊ + σ² has a profound interpretation. Let's unpack each term.

Term 1: Epistemic Uncertainty (x₊ᵀSₙx₊)

This term arises because we don't know the true weights. Even with infinite precision observations (σ² = 0), we'd still have prediction uncertainty if we're uncertain about w.

Key properties:

• Depends on x₊: Uncertainty is input-dependent
• Decreases with data: As N → ∞, Sₙ → 0, and this term vanishes
• Highest far from training data: Where we've seen few examples, weight uncertainty is high
• Quadratic form: x₊ᵀSₙx₊ measures how x₊ "projects" onto directions of high weight uncertainty

Term 2: Aleatoric Uncertainty (σ²)

This is irreducible noise—the inherent randomness in the data-generating process. Even with known true weights, observations would be noisy.

Key properties:

• Constant across inputs: Same noise level everywhere
• Cannot be reduced by more data: It's fundamental, not epistemic
• Sets a floor on predictive variance: σ²₊ ≥ σ² always
• Often called "observation noise" or "measurement noise"

The predictive distribution enables principled prediction intervals—ranges that contain the true value with specified probability.

Credible Intervals:

For a Gaussian predictive distribution 𝒩(μ₊, σ²₊), the (1-α) credible interval is:

$$\text{CI}{1-\alpha} = \left[ \mu* - z_{\alpha/2} \sigma_, \quad \mu_ + z_{\alpha/2} \sigma_* \right]$$

where z_{α/2} is the standard normal quantile. Common values:

• 95% CI: z = 1.96 ≈ 2
• 99% CI: z = 2.576
• 68% CI: z = 1 (one standard deviation)

Interpretation (Bayesian):

The 95% credible interval means: "Given our prior and observed data, there's a 95% probability the true y₊ lies in this interval."

This is subtly different from frequentist confidence intervals, which make statements about long-run coverage rather than probability of the true value.

Key Observations from the Visualization:

1. Uncertainty grows away from data: As we move from the training region, the predictive variance increases
2. Narrowest at data-dense regions: Where we have many observations, we're most confident
3. Never zero variance: Even at training points, σ² sets a floor due to observation noise
4. Intervals are honest: The coverage is calibrated—95% intervals really do contain 95% of test points on average

A distinctive feature of Bayesian prediction is that uncertainty varies across the input space. The epistemic component x₊ᵀSₙx₊ depends on where we're predicting.

Geometric Interpretation:

The posterior covariance Sₙ defines an ellipsoid of weight uncertainty. The quadratic form x₊ᵀSₙx₊ measures how x₊ "stretches" this ellipsoid:

• If x₊ aligns with directions of high uncertainty in Sₙ: Large contribution
• If x₊ is orthogonal to uncertain directions: Small contribution

Intuitive Understanding:

Imagine linear regression with two features x₁ and x₂:

• Training data only covers x₁ ∈ [-1, 1], x₂ ∈ [-1, 1]
• The weight for x₁ is well-estimated (small variance)
• The weight for x₂ is less certain (larger variance)

For a test point with large x₂, predictions will be uncertain because the uncertain weight w₂ is multiplied by a large factor.

Extrapolation Warning:

This is why Bayesian methods naturally warn about extrapolation:

• Test points far from training data → large x₊ₜSₙx₊ → wide uncertainty bands
• The model honestly admits "I haven't seen data like this"
• Classical regression gives confident (but unreliable) extrapolations

Often we need predictions at multiple test points simultaneously. Let X₊ ∈ ℝᴺˣᴰ be a matrix of N test inputs. The joint predictive distribution is:

$$p(\mathbf{y}* | \mathbf{X}, \mathbf{y}, \mathbf{X}) = \mathcal{N}(\mathbf{y}_ | \boldsymbol{\mu}*, \boldsymbol{\Sigma}*)$$

where:

Predictive Mean Vector: $$\boldsymbol{\mu}* = \mathbf{X}*\mathbf{m}_N$$

Predictive Covariance Matrix: $$\boldsymbol{\Sigma}* = \mathbf{X}\mathbf{S}N\mathbf{X}^\top + \sigma^2\mathbf{I}{N*}$$

Structure of the Covariance:

The diagonal elements are the individual predictive variances: $$(\boldsymbol{\Sigma}*){ii} = \mathbf{x}_{*i}^\top\mathbf{S}N\mathbf{x}{*i} + \sigma^2$$

The off-diagonal elements capture correlations between predictions: $$(\boldsymbol{\Sigma}*){ij} = \mathbf{x}_{*i}^\top\mathbf{S}N\mathbf{x}{*j}$$

Note: Noise is independent across test points, so it only appears on the diagonal.

Why Correlations Matter:

Predictions at different test points are correlated through shared weight uncertainty. If we're uncertain about w₁, predictions for all test points with large x₁ will be jointly uncertain—they'll tend to be wrong in the same direction.

This correlation structure is crucial for:

• Sampling consistent function realizations
• Propagating uncertainty through downstream computations
• Understanding how errors compound in multi-step predictions

Sampling from the Predictive Distribution:

To generate samples of consistent predictions:

$$\mathbf{y}* \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma}_)$$

Or equivalently:

1. Sample weights: w ~ 𝒩(mₙ, Sₙ)
2. Compute noiseless predictions: f₊ = X₊w
3. Add noise: y₊ = f₊ + ε, where ε ~ 𝒩(0, σ²I)

Frequentist linear regression also provides prediction intervals. How do they compare to Bayesian credible intervals?

Frequentist Prediction Interval:

For OLS with Gaussian errors, the (1-α) prediction interval at x₊ is:

$$\hat{y}* \pm t{n-p, \alpha/2} \cdot \hat{\sigma} \cdot \sqrt{1 + \mathbf{x}*^\top(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{x}*}$$

where:

• t_{n-p, α/2} is the t-distribution quantile (accounts for estimating σ)
• σ̂ is the estimated noise standard deviation
• The term under the square root looks remarkably similar to our Bayesian formula!

Key Differences:

When They Agree:

With a flat/vague prior (α → 0) and known σ²:

• Bayesian predictive variance → Frequentist formula
• The t-correction vanishes with known σ²
• Intervals become identical

When Bayesian Is Better:

1. Small samples: Prior prevents overfitting and unstable variance estimates
2. Near-collinear features: Prior regularization stabilizes (XᵀX)⁻¹
3. Informative prior available: Domain knowledge improves predictions
4. Uncertainty propagation: Full distribution enables proper uncertainty propagation through complex pipelines

When Frequentist Is Better:

1. Unknown prior: If no prior knowledge exists, frequentist avoids arbitrary choices
2. Coverage guarantees: If you need formal frequentist coverage guarantees
3. Computational simplicity: Fewer choices to make, standard implementation

Implementing Bayesian prediction correctly requires attention to several details.

Efficient Variance Computation:

For many test points, computing x₊ᵀSₙx₊ individually is inefficient. Use vectorized computation:

# Slow: loop over test points
var = np.array([x.T @ S_N @ x for x in X_test])

# Fast: vectorized
var = np.sum((X_test @ S_N) * X_test, axis=1)

Numerical Stability:

Ensure Sₙ is computed stably (via Cholesky decomposition). Small negative eigenvalues from numerical error can cause problems when computing variances.

Handling Unknown σ²:

If σ² is estimated from data, the predictive distribution is technically a Student-t, not Gaussian. For practical purposes with reasonable N:

1. Estimate σ² from residuals: σ̂² = ||y - X****mₙ||² / (N - D)
2. Use this estimate in the Gaussian predictive (slight overconfidence for small N)
3. Or: place a prior on σ² and integrate it out (see advanced Bayesian texts)

The predictive distribution is the practical output of Bayesian linear regression—it tells us what we can expect for new inputs, with honest uncertainty quantification.

What's Next:

We've seen that Bayesian linear regression provides uncertainty estimates. But how do we use this uncertainty in practice? When is uncertainty high vs. low? How do we communicate it to stakeholders? The next page explores uncertainty quantification in depth—practical methods for interpreting, visualizing, and acting on predictive uncertainty.