Bayesian Linear Regression

Simple linear regression¶

We start with a single predictor $x$ and a continuous response $y$.

The deterministic regression line is

with intercept $\beta_0$ and slope $\beta_1$.

To make this into a probabilistic model, we add a noise term:

Equivalently, the likelihood for each observation is

In vector–matrix form, with design matrix

the model is

Likelihood and ordinary least squares (OLS)¶

Under the normal likelihood

the log-likelihood (up to additive constants) is

Maximizing this log-likelihood w.r.t. $\boldsymbol{\beta}$ (for fixed $\sigma$) is equivalent to minimizing the sum of squared residuals:

In matrix form, the OLS solution has a closed form:

provided $X^\top X$ is invertible.

OLS is the maximum likelihood estimator for $\boldsymbol{\beta}$ under the normal linear model with constant variance $\sigma^2$.

Bayesian linear regression¶

In the Bayesian view, the parameters $(\beta_0,\beta_1,\sigma)$ are random variables with a prior distribution.

A common prior choice (conceptually) is:

• Regression coefficients:

  $$\boldsymbol{\beta} \sim \mathcal{N}(\boldsymbol{\mu}_0,\; \Sigma_0),$$

  (10)

• Noise scale:

  $$\sigma \sim \text{some prior on } (0,\infty) \quad\text{(e.g. half-normal, exponential, or half-Cauchy).}$$

  (11)

Given the data $(\mathbf{y}, X)$, the posterior over parameters is

Using Bayes’ theorem explicitly,

For conjugate priors (normal–inverse-gamma), the posterior has an analytic form; in practice we often let PyMC/Bambi do the MCMC sampling from $p(\boldsymbol{\beta},\sigma \mid \mathbf{y},X)$.

Posterior predictive distribution for regression¶

For a new input $x_{\text{new}}$, the model implies

The posterior predictive distribution, which averages over uncertainty in $(\beta_0,\beta_1,\sigma)$, is

In practice, with posterior samples $\{(\beta_0^{(s)},\beta_1^{(s)},\sigma^{(s)})\}_{s=1}^{S}$, we approximate this by simulation:

1. For each posterior draw $s$,

   $$y_{\text{new}}^{(s)} \sim \mathcal{N}\big(\beta_0^{(s)} + \beta_1^{(s)} x_{\text{new}}, (\sigma^{(s)})^2\big).$$

   (16)

2. The sample $\{y_{\text{new}}^{(s)}\}_{s=1}^{S}$ approximates $p(y_{\text{new}} \mid x_{\text{new}}, \mathbf{y}, X)$.

The posterior predictive distribution combines:

• Aleatoric uncertainty: the inherent noise $\varepsilon$ in $y$,

• Epistemic uncertainty: uncertainty about $(\beta_0,\beta_1,\sigma)$ due to finite data.

Epistemic vs aleatoric uncertainty in regression¶

In the linear model

we distinguish:

• Aleatoric uncertainty: variability in $Y$ even if $(\beta_0,\beta_1,\sigma)$ were known exactly. In this model, for a fixed $x$, this is just $\operatorname{Var}(Y \mid X=x,\beta_0,\beta_1,\sigma) = \sigma^2$.

• Epistemic uncertainty: additional variability in predictions due to uncertainty in $(\beta_0,\beta_1,\sigma)$ as captured by the posterior $p(\beta_0,\beta_1,\sigma \mid \mathbf{y},X)$.

For the posterior predictive variance at a given $x_{\text{new}}$, the law of total variance gives

where $\theta = (\beta_0,\beta_1,\sigma)$.

• The first term corresponds to aleatoric variance (expected $\sigma^2$).

• The second term corresponds to epistemic variance (spread of the regression line due to parameter uncertainty).

With many data points covering a range of $x$, the epistemic term often becomes small near the data but remains large when extrapolating far away from observed $x$ values.

Gaussian linear regression¶

Linear regression is perhaps the most frequently used technique in applied statistics for modelling the relationship between set of a predictor variables $\left\lbrace \mathbf{x}_i \right\rbrace_{i=1}^N$ and a response variable $\left\lbrace y_i \right\rbrace_{i=1}^N$. More formally, let $\mathbf{X} \in \mathbb{R}^{N \times D}$ be a design matrix and let $\mathbf{y} \in \mathbb{R}^N$ be the response variables collected into a single vector, then the linear regression model is given by

where $\mathbf{w} \in \mathbb{R}^D$ is the regression weights and $\mathbf{\epsilon} \in \mathbb{R}^N$ is the observation noise vector.

Assuming isotropic Gaussian noise and imposing a multivariate Gaussian prior on $\mathbf{w} \sim \mathcal{N}\left(\mathbf{m}, \mathbf{S}\right)$ gives rise to the following joint distribution

In this notebook, we will use the following simple model as running example:

That is, the parameters are $\mathbf{w} = \left[a, b\right]$, where $a$ and $b$ are the slope and intercept of the line, respectively. Furthermore, we will assume $\mathbf{m} = \mathbf{0}$, $\mathbf{S} = \mathbf{I}$, and $\sigma^2 = 2$.

We will use the following dataset with $N = 8$ data points

The prior, likelihood and posterior¶

This section illustrates the basic building blocks of the linear model.

Reference: Rasmussen and Williams, Section 2.1, for the derivation of these functions (http://www.gaussianprocess.org/gpml/ ).

The multivariate Gaussian density function¶

The multivariate Gaussian density function, implemented using numpy.

• Computes the (log-)density of a multivariate Gaussian $\mathcal{N}(x \mid \mu, \Sigma)$ using a numerically stable Cholesky decomposition of the covariance matrix.

• Adds a small diagonal jitter to $\Sigma$ to ensure positive definiteness and avoid numerical issues during factorization.

• Uses the Cholesky factor $L$ to efficiently solve $v = L^{-1}(x - \mu)$, which is required to compute the quadratic form $(x - \mu)^\top \Sigma^{-1} (x - \mu) = \|v\|^2$.

• Combines the constant term, log-determinant term (computed from $\log \det \Sigma = 2 \sum \log L_{ii}$), and quadratic term to return either the log-density or the density.

The Log density functions¶

The Log density of the (i) prior, (ii) likelihood, and (iii) posterior (likelihood times prior), implemented using the mvn_pdf function.

The Visualization¶

Visualization of the prior density, the likelihood and the posterior as a function of $a$ and $b$.

Analytical expression of the posterior distribution¶

Derivation of the posterior $p(\mathbf{w}\big|\mathbf{y}) = \mathcal{N}\left(\mathbf{w}\big|\mu, \Sigma\right)$ yields:

Reference: Rasmussen and Williams, Section 2.1.1, for the derivation of these functions (http://www.gaussianprocess.org/gpml/ ).

The Implementation¶

Implementation of the function for computing the analytical posterior distribution.

• Computes the inverse of the posterior covariance matrix $\Sigma_{\text{post}}^{-1} = S^{-1} + \frac{1}{\sigma^2} X^\top X$ using the inverse of the prior covariance Sinv = np.linalg.inv(S) and the design matrix X = design_matrix(x).

• Forms the Cholesky factorization of the inverse of the inverse posterior covariance matrix using sc.linalg.cho_factor(Sigmainv), enabling stable and efficient solutions for both the posterior mean and covariance without explicitly inverting matrices.

• Computes the scaled right-hand side $\frac{1}{\sigma^2} X^\top y + S^{-1} m$ using np.dot(X.T, y)/sigma2 + np.linalg.solve(S, m); the call to np.linalg.solve computes $S^{-1}m$ more accurately than explicitly forming $S^{-1}$.

• Solves the linear system $\Sigma_{\text{post}}^{-1}\,\mu_{\text{post}} = \frac{1}{\sigma^2}X^\top y + S^{-1} m$ via sc.linalg.cho_solve(L, scaled_mu), yielding the posterior mean vector mu.

• Computes the posterior covariance matrix $\Sigma_{\text{post}} = (\Sigma_{\text{post}}^{-1})^{-1}$ by solving cho_solve(L, I) where I = np.identity(len(m)); this avoids explicit matrix inversion and ensures numerical stability.

The Application¶

The posterior mean and covariance for the data set $(\mathbf{x}, \mathbf{y})$ are calculated using the function above. The plot of the resulting density matches with the posterior distribution from section 1.3.3.

Combining all of the above¶

All the functions implemented above are used to compute the posterior distribution $p\left(\mathbf{w}|\mathbf{y}\right)$ for various sizes of $N$. The cell below analyses the data set using $N$ data points, where the columns visualise (i) data and sample functions, (ii) the prior, (iii) the likelihood and (iv) posterior.

• generate_mvn_samples samples from a multivariate normal $\mathcal{N}(\mu, \Sigma)$ by computing a Cholesky factor L = np.linalg.cholesky(Sigma + jitter*np.identity(D)), where a small jitter term ensures that $\Sigma$ is numerically positive definite.

• Standard normal samples zs are drawn with np.random.normal(0, 1, size=(D, M)), giving a $D \times M$ matrix of i.i.d. $\mathcal{N}(0, 1)$ variables; multiplying by L linearly transforms them to have covariance $\Sigma$.

• The final samples are constructed as fs = mu + np.dot(L, zs), implementing $f = \mu + L z$; this uses NumPy broadcasting to add the (D)-dimensional mean vector mu to each of the (M) sample columns.

• In compute_f_posterior, the design matrix X = np.column_stack((x, np.ones(len(x)))) encodes a linear model $f(x) = w_1 x + w_0$, with the first column for the input $x$ and the second column for the intercept term.

• Given a posterior over weights $w \sim \mathcal{N}(\mu, \Sigma)$, the induced posterior over function values is $f(x) \sim \mathcal{N}(X\mu,\; X \Sigma X^\top)$; this is implemented by mu_f = np.dot(X, mu) for the mean.

• The pointwise posterior variances $\mathrm{Var}[f(x_i)]$ are the diagonal elements of $X \Sigma X^\top$; these are extracted via np.diag(np.dot(np.dot(X, Sigma), X.T))[:, None], with [:, None] reshaping the vector into an $(N \times 1)$ column.

The Quiz¶

1. How does the blue lines in the first column relate to the blue dots in the fourth column?¶

Each blue line in the first column is a straight line in the form of linear regression: $y = slope * x + intercept$, whose parameters slope and intercept correspond to one blue point in the fourth column, which is a pair of slope and coefficient in the posterior contour.

2. What happens as the number of samples increase?¶

Prior contours stay the same as the number of samples increases.

Likelihood contours become dramatically narrower as the number of samples increase

Posterior contours becomes smaller like likelihood contours, but not so significant as the number of samples increases.

Posterior draws: the range of the posterior draws becomes narrower, which tries to fit the data as much as possible as the number of samples increases.

The reason for this behavior is the following: more data points imply more information, which helps reduce uncertainty in regression and makes likelihood and posterior draws narrower

3. Why does the prior have so little effect in the last row?¶

The amount of data overwhelms the prior: in the last case, the amount of data available is large enough such that the likelihood completely overwhelms the prior, making the prior effect negligible in posterior draws.

Prior is not informative: If the prior is not informative, i.e., it does not provide much useful information about the distribution of the parameters, then it is unlikely to have a significant effect on the posterior. This is true in this case as the prior has zero mean and variance of one, which provideds little prior knowledge about the underlying data structure when there are many datapoints already.

House Price Example¶

Exercises¶

Problem 1.1 : Modeling Credit data with Simple Linear Regression¶

In this exercise, you will perform the basic steps of Bayesian linear regression. To this end, you will work with the Credit dataset that you can load with Pandas from credit_data.csv.

a) Data exploration¶

Load the dataset and display the first five rows. Identify the predictor ($x$) and the response ($y$) variables.

• $x$ = Limit

• $y$ = Balance

b) Data visualization¶

Plot the data points of Balance vs. Limit. Describe the relationship you observe. Does a linear model seem appropriate?

The relationship appears to be linear - except for points close to zero. They seem rather arbitrarily distributed.

c) Model definition¶

Consider the linear model $y = ax + b + \varepsilon,\ \ \varepsilon \sim \mathcal{N}(0,\sigma^2).$ Explain the meaning of the parameters $a$, $b$, and $\sigma^2$.

• $a$ is the slope - change in Balance per unit of Limit

• $b$ is the intercept - expected Balance when Limit is zero

• $\sigma^2$ is the noise variance - the unexplained variation in the data

d) Choosing priors¶

Plot the prior distribution using the functions derived in the lectures. Assume a Gaussian prior for $\mathbf{w} = [a,b]^T$: $p(\mathbf{w}) = \mathcal{N}(m,\Sigma_p)$, $\Sigma_p = 10^6 I$. What does this prior represent? How would the model change if $\Sigma_p$ were smaller, such as $\Sigma_p = 10^2 I$? How would that affect the posterior and the predictive region?

Using smaller diagonal values in $\Sigma_p$ means, there is less variance in the parameters. Less variance in turn means that the posterior distribution is more peaked around the true values of $a$ and $b$ - because there is less uncertainty.

e) Computing the posterior¶

Plot the likelihood and the posterior distribution using the functions derived in the lectures. Compute the posterior mean and covariance of $\mathbf{w}$. Compare and interpret the shapes of the prior, likelihood, and posterior contour plots.

f) Posterior credibility distribution¶

Using the posterior parameters, calculate the posterior predictive mean and posterior variance: $\mu_f(x)$, $\mathrm{Var}[f(x)]$, and plot these together with the scatter plot. Explain what each component (mean and variance) represents. Explain why the posterior credibility region widens toward the edges of the data.

• $\mu_f(x)$ is the most likely linear regression function (the mean of all possible linear functions according to our model)

• $\mathrm{Var}[f(x)]$ is the expected uncertainty at every value of x

g) Predictive region¶

By computing the predictive distribution, plot the predictive region along with the scatter plot. Experiment with different values of $\sigma^2$ (e.g. $10^3, 10^4, 10^5, 10^6$). How does changing $\sigma^2$ affect the posterior credibility region and predictive region? Which value seems most reasonable for this dataset?

Using $\sigma^2 = 1e5$ looks best; it explains most of the variance and its neither too broad nor too narrow.

h) Model performance¶

Compute and interpret the following metrics for your model: RMSE, MAE, and $R^2$. What do these values indicate about model accuracy?

i) Model comparison¶

Suppose we add a polynomial term $x^2$ to the model. How might this affect the results? Would Bayesian regression still prevent overfitting?

A polynomial model is more powerfull, but risks overfitting. Choosing a corresponding Bayesian prior can act as a regularization technique to prevent overfitting.

j) Reflection¶

Summarize how the prior, likelihood, and posterior interact. How do $\Sigma_p$ and $\sigma^2$ influence the model predictions and uncertainty? Why is it important to interpret both the point estimates and the uncertainty bands when communicating regression results?

The prior describes the belief before seeing any data (our parameter initializations), the likelihood describes the probability of observing the given data using our belief (aka. parameters) and finally the posterior is our updated belief (parameters) after we’ve seen the data.

The $\Sigma_p(S)$ describes the uncertainty of our prior (belief): Its diagonal values describe the variance for each parameter (higher variance = higher std. = higher uncertainty).

Lastly, $\sigma^2$ describes the noise in our data (aka. aspects our model cannot explain). This covers the overall predictive uncertainty.

Bayesian linear regression with Bambi¶

Bambi is a high-level interface for Bayesian regression models built on top of PyMC. It uses an R-style formula interface.

For simple linear regression of body fat on BMI, one might specify

import bambi as bmb

model = bmb.Model("BodyFat ~ BMI", data=df)
results = model.fit()

Conceptually, this is equivalent to specifying the normal linear model

with weakly informative priors chosen automatically for $\beta_0$, $\beta_1$, and $\sigma$.

Bambi:

• Builds the corresponding PyMC model and runs MCMC,

• Returns posterior samples in an InferenceData object,

• Provides tools for posterior summaries, HDIs, and predictive intervals.

More complex models (multiple predictors, categorical variables, interactions, non-normal likelihoods) can be specified by extending the formula or choosing a different family (likelihood).

Posterior HDIs for regression coefficients¶

From the posterior samples of regression coefficients, we can compute highest density intervals (HDIs).

For a coefficient $\beta_j$ and a desired probability mass $p$ (e.g. 0.90 or 0.94), the $p$-HDI is an interval $[a,b]$ such that

• $P(\beta_j \in [a,b] \mid \mathbf{y},X) = p$,

• posterior density inside the interval is higher than outside.

Practical interpretation:

• If a $94\%$ HDI for $\beta_1$ does not include zero, we say that zero is not a plausible value for $\beta_1$ under the model and prior. In informal language: the effect is “significant at the 94% level”.

• If the HDI does include zero, then values on both sides (positive and negative effect) are compatible with data and prior; we cannot confidently say the effect is non-zero.

Bambi and ArviZ compute HDIs directly from the posterior samples.

Bayesian $R^2$ for regression¶

In classical linear regression, the coefficient of determination $R^2$ measures the fraction of variance explained by the model.

In a Bayesian setting, we can define a Bayesian $R^2$ (Gelman et al.), based on posterior samples.

For each posterior draw $s = 1,\dots,S$:

1. Compute the vector of predicted means (fitted values)

   $$\boldsymbol{\mu}^{(s)} = \big(\mu_1^{(s)}, \dots, \mu_n^{(s)}\big)^\top, \quad \mu_i^{(s)} = \mathbb{E}[Y_i \mid x_i,\theta^{(s)}].$$

   (25)

2. Compute residuals

   $$\mathbf{e}^{(s)} = \mathbf{y} - \boldsymbol{\mu}^{(s)}.$$

   (26)

3. Compute the explained variance and residual variance for draw $s$:

   $$V_{\text{fit}}^{(s)} = \operatorname{Var}_i\big(\mu_i^{(s)}\big), \qquad V_{\text{res}}^{(s)} = \operatorname{Var}_i\big(e_i^{(s)}\big),$$

   (27)

   where the variance is over data points $i = 1,\dots,n$.

Define the per-draw $R^2$ as

The Bayesian $R^2$ is then summarized by the posterior distribution of $\{R^{2,(s)}\}_{s=1}^{S}$, e.g. by its mean and an HDI.

This approach acknowledges:

• Uncertainty in the regression function (epistemic),

• Uncertainty in the noise (aleatoric).

Multiple linear regression¶

With $p$ predictors, the linear regression model generalizes to

In matrix notation:

where now $X$ is an $n \times (p+1)$ matrix (including the intercept column of ones).

Bayesian multiple linear regression:

• Place a prior on the coefficient vector $\boldsymbol{\beta}$ and $\sigma$,

• Use MCMC (via PyMC/Bambi) to obtain posterior samples $p(\boldsymbol{\beta},\sigma \mid \mathbf{y},X)$,

• Summarize and interpret parameter and predictive uncertainty.

Adding predictors generally:

• Increases the potential explained variance,

• Increases epistemic uncertainty (more parameters to estimate),

• Risks overfitting if not regularized or if model is too flexible for the data size.

Categorical predictors and interactions¶

Categorical predictors (e.g. Student yes/no, or Ethnicity with several categories) are incorporated via dummy variables.

Example: binary predictor Student (0 = no, 1 = yes):

• $\beta_2$ shifts the mean response up or down for students compared to non-students.

Interactions between predictors allow the effect of one predictor to depend on another.

For two continuous predictors $x_1$ (TV budget) and $x_2$ (radio budget), a model with interaction is

• $\beta_3$ captures synergy (if positive) or redundancy (if negative) between channels.

In Bambi’s formula syntax:

• Main effects: Sales ~ TV + radio

• With interaction: Sales ~ TV * radio

where TV * radio expands to TV + radio + TV:radio.

Robust linear regression with Student’s $t$ likelihood¶

The normal likelihood is sensitive to outliers: large residuals are heavily penalized, which can pull the regression line towards outlying points.

To reduce sensitivity to outliers, we can model residuals with a Student’s $t$ distribution:

where

• $\nu > 0$ is the degrees of freedom parameter (shape),

• smaller $\nu$ yields heavier tails than the normal distribution.

As $\nu \to \infty$, the Student’s $t$ distribution converges to the normal:

A Bayesian robust regression model might use:

• Normal priors for coefficients $\boldsymbol{\beta}$,

• A prior on $\sigma > 0$ (half-normal, exponential, etc.),

• A prior on $\nu$ (e.g. $\nu \sim \operatorname{Gamma}(s,r)$).

Compared to the normal-likelihood model, the $t$-likelihood:

• Gives less weight to extreme residuals,

• Produces regression lines less distorted by outliers,

• Often yields broader predictive intervals, reflecting uncertainty about outliers.

Polynomial regression¶

To capture simple nonlinear relationships between a single predictor $x$ and response $y$, we can use polynomial regression.

For a degree-3 polynomial:

This is still linear in the parameters $\boldsymbol{\beta}$, so it fits into the linear regression framework by defining transformed predictors:

However:

• Polynomial terms can be strongly correlated (e.g. $x$ and $x^2$), causing numerical issues.

• High-degree polynomials can overfit heavily, especially with few data points.

Bambi’s poly(x, d) construct can create orthogonal polynomials, which reduce the correlations between polynomial terms and improve sampling stability.

Overfitting and epistemic uncertainty in polynomial models¶

With few data points and a highly flexible model (e.g. high-degree polynomial), the model can fit the training data very closely, but:

• The regression function becomes wiggly and unstable,

• Posterior predictive intervals can become very wide away from data points,

• Epistemic uncertainty increases, reflecting that many different polynomials are compatible with the data.

Bayesian inference:

• Makes overfitting visible in the posterior distribution of functions and predictive intervals,

• Encourages the use of priors (shrinkage) to regularize complex models,

• Allows the use of model comparison (e.g. Bayes factors, ELPD) to decide whether increased flexibility is justified.

Shrinkage via normal priors and L2 regularisation¶

Consider a linear regression model with predictors $X$ (possibly including polynomial terms) and coefficients $\boldsymbol{\beta}$. Suppose we place an independent normal prior on each coefficient:

and a normal likelihood

Ignoring constants, the log posterior is

Up to additive constants, this becomes

The maximum a posteriori (MAP) estimate $\hat{\boldsymbol{\beta}}_{\text{MAP}}$ maximizes this log posterior, i.e. minimizes

Up to a constant factor, this is equivalent to ridge regression (L2 regularisation):

with penalty parameter $\lambda = \sigma^2 / \tau^2$.

Thus:

• L2 regularisation corresponds to a normal prior on coefficients centered at zero.

• Stronger prior shrinkage (smaller $\tau$) corresponds to larger $\lambda$ and stronger penalization of large coefficients.

L1 regularisation and Laplace prior¶

Analogously, L1 regularisation (lasso) corresponds to a Laplace (double-exponential) prior on the coefficients.

The Laplace prior for a scalar coefficient $\beta$ with scale $b > 0$ has density

Ignoring constants, the log prior is

so for a vector of coefficients $\boldsymbol{\beta}$ with independent Laplace priors, the negative log prior adds a term proportional to $\|\boldsymbol{\beta}\|_1$.

The MAP estimator then minimizes

which is the lasso objective.

In summary:

• Normal prior on $\boldsymbol{\beta}$ $\Rightarrow$ ridge (L2) penalty at MAP.

• Laplace prior on $\boldsymbol{\beta}$ $\Rightarrow$ lasso (L1) penalty at MAP.

Regularisation can therefore be understood as a natural consequence of using informative priors in Bayesian models.

Model selection versus regularisation¶

When dealing with flexible models (e.g. polynomials of various degrees), there are two complementary strategies:

1. Model selection
   Compare several models (e.g. linear, quadratic, cubic) using:

   • Bayes factors (marginal likelihoods),

   • ELPD / LOO,

   • Posterior predictive checks,

   • Predictive RMSE/MAE.

   Then choose the model that best balances fit and complexity.

2. Regularisation (shrinkage)
   Use a richer model (e.g. higher-degree polynomial) but place strong priors that shrink higher-order coefficients towards zero. This keeps the model flexible but penalizes unnecessary complexity.

Often, a combination of both is used:

• Start with a simple model,

• Add complexity if strongly supported by posterior predictive checks and model comparison,

• Use priors to regularize complex models and avoid overfitting.

In Bayesian regression, regularisation is not an ad-hoc fix but an integral part of the modelling philosophy through the choice of priors.