Bayesian linear regression: not one best-fit line, but a distribution of lines with honest error bars

Ordinary least squares draws one best-fit line and stops. It never tells you how much to trust that line, or where it is basically guessing. Bayesian linear regression keeps a whole distribution over lines — a mean line, a cloud of plausible lines, and a shaded band that balloons out where you have no data and tightens where points pile in. I built an interactive demo where you click to drop points and watch the band breathe, and every number is solved in closed form.

The model, in one sentence

The line is y = w₀ + w₁·x with Gaussian noise of precision β (noise variance is 1/β). Before seeing data, we put a Gaussian prior on the weights, w ~ N(0, α⁻¹I) — the prior precision α gently pulls the line toward flat. Those two knobs, α and β, are the only sliders in the demo; the entire model lives in them.

import numpy as np
# a straight line is just two basis functions: the constant 1, and x
Phi   = np.c_[np.ones_like(x), x]   # design matrix, one row per point
alpha = 2.0   # prior precision -> prior variance 1/alpha
beta  = 6.0   # noise precision -> noise std 1/sqrt(beta)

The posterior is closed form — no sampling

A Gaussian prior times a Gaussian likelihood is again Gaussian (conjugacy), so there is nothing to sample. The posterior precision is Aₙ = αI + βΦᵀΦ, its covariance is the inverse, and its mean is mₙ = β·Σₙ·Φᵀy. This is exactly what runs on every slider move and every click in the demo.

A = alpha * np.eye(Phi.shape[1]) + beta * Phi.T @ Phi   # posterior precision
S = np.linalg.inv(A)                                    # posterior covariance
m = beta * S @ Phi.T @ y                                # posterior mean weights

The band is the predictive variance

For a new input x* the prediction is itself Gaussian. Its mean is mₙᵀφ(x*); its variance is 1/β + φ(x*)ᵀΣₙφ(x*) — a noise floor plus the model's own uncertainty. That second term is a parabola in x*, which is the whole story of the shaded band: it is narrowest near the centre of mass of your data and flares outward as you extrapolate. Add points and Σₙ shrinks, so the band thins — but never below the noise floor 1/√β, because that noise is irreducible.

def predict(xs):
    P    = np.c_[np.ones_like(xs), xs]
    mean = P @ m
    var  = 1.0/beta + np.einsum('ij,jk,ik->i', P, S, P)  # 1/beta + phiᵀ S phi
    return mean, np.sqrt(var)                            # std = the band half-width

This is the payoff over OLS in one expression: a model that knows where it does not know. Right over your data it is confident; far from it, it flares out and tells you extrapolation is a guess.

Sampling whole lines makes it visceral

The band summarises the uncertainty, but drawing actual lines from the posterior makes it click. Sample weight vectors from N(mₙ, Σₙ) and plot each as a line — where they fan apart is where the model is unsure. These are the grey lines in the demo (drawn there with a 2×2 Cholesky factor). There is also a second view in weight space: every line is a single point (w₀, w₁), the prior is a big ring centred at zero, and the posterior is a tight ellipse that you can watch shrink and slide onto the true weights as you add data.

W = np.random.multivariate_normal(m, S, size=20)   # 20 plausible weight vectors
for w0, w1 in W:
    plt.plot(xs, w0 + w1 * xs, color="grey", alpha=0.3)   # one sampled line each
plt.plot(xs, m[0] + m[1] * xs, color="dodgerblue", lw=2)  # posterior mean

Where OLS, ridge, and full Bayes sit

It is worth seeing how these connect. Least squares is the limit α → 0 (a flat prior): one naked line, no error bars. Ridge regression is the MAP point estimate of this same model (α is the ridge penalty). Full Bayes keeps the whole posterior, so you get calibrated uncertainty everywhere instead of a single point.

from sklearn.linear_model import BayesianRidge
model = BayesianRidge()          # even learns alpha & beta from the data (evidence)
model.fit(X, y)
mean, std = model.predict(X_test, return_std=True)   # predictions WITH error bars

BayesianRidge does everything above and estimates α and β by maximising the marginal likelihood, so there is no hand-tuning. Reach for it whenever a prediction interval matters — forecasting, active learning, anything safety-adjacent — and on small or wide data where it self-regularises. And when you want curved fits with the same honest bands, you swap the basis functions: the exact same posterior math over Gaussian basis functions is a Gaussian Process.

Click anywhere on the chart to drop points, drag the α and β sliders, and watch the band widen where data is thin:
https://dev48v.infy.uk/ml/day36-bayesian-linear-regression.html