Bayesian Linear Regression

Bayesian linear regression extends classical OLS by placing probability distributions on the regression coefficients and error variance, yielding full posterior distributions rather than point estimates.

The Bayesian Regression Model

Conjugate Prior: Normal-Inverse-Gamma

Posterior Inference for Coefficients

The marginal posterior of is a multivariate -distribution. For large , this approaches a Gaussian:

Credible Intervals vs. Confidence Intervals

Bayesian Prediction

Under the Normal-Inverse-Gamma prior, the posterior predictive is a -distribution:

The inside the parentheses accounts for parameter uncertainty — this is wider than the OLS prediction interval when the prior is weak and is small.

Bayesian Model Comparison

Python Implementation

import numpy as np
import pymc as pm
import arviz as az
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(42)

# --- Generate synthetic data ---
n = 80
X = np.random.randn(n, 2)
X = np.column_stack([np.ones(n), X])  # Add intercept
true_beta = np.array([1.5, -2.0, 0.8])
sigma_true = 1.5
y = X @ true_beta + np.random.randn(n) * sigma_true

# --- Bayesian Regression with PyMC ---
with pm.Model() as bayes_reg:
    # Priors
    beta = pm.Normal('beta', mu=0, sigma=10, shape=3)
    sigma = pm.HalfNormal('sigma', sigma=5)
    
    # Likelihood
    mu = pm.math.dot(X, beta)
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)
    
    # Posterior sampling
    trace = pm.sample(2000, tune=1000, chains=4, return_inferencedata=True)

# --- Posterior summaries ---
print(az.summary(trace, var_names=['beta', 'sigma']))

# --- Credible intervals ---
for i, name in enumerate(['Intercept', 'X1', 'X2']):
    post_samples = trace.posterior['beta'].values[:, :, i].flatten()
    ci = np.percentile(post_samples, [2.5, 97.5])
    print(f"{name}: {ci[0]:.3f} to {ci[1]:.3f} (true={true_beta[i]:.3f})")

# --- Trace plots ---
az.plot_trace(trace, var_names=['beta', 'sigma'])
plt.tight_layout()
plt.savefig('bayesian_regression_trace.png', dpi=150)
plt.show()

# --- Posterior predictive ---
with bayes_reg:
    ppc = pm.sample_posterior_predictive(trace)
    
y_pred = ppc.posterior_predictive['y_obs'].values.reshape(-1, n)
y_mean = y_pred.mean(axis=0)
y_low = np.percentile(y_pred, 2.5, axis=0)
y_high = np.percentile(y_pred, 97.5, axis=0)

plt.figure(figsize=(10, 6))
order = np.argsort(X[:, 1])
plt.scatter(X[:, 1], y, alpha=0.5, label='Observed')
plt.plot(X[order, 1], y_mean[order], 'r-', label='Posterior mean')
plt.fill_between(X[order, 1], y_low[order], y_high[order], alpha=0.3, label='95% credible band')
plt.xlabel('X₁')
plt.ylabel('y')
plt.title('Bayesian Linear Regression — Posterior Predictive')
plt.legend()
plt.savefig('bayesian_regression_ppc.png', dpi=150)
plt.show()

Related Topics

• See Simple Linear Regression for classical OLS foundations
• See Hypothesis Testing for frequentist inference
• See Bayesian Statistics for the Bayesian–frequentist debate
• See Hierarchical Bayesian Models for multilevel extensions
• See MCMC Diagnostics for ensuring convergence of posterior samples

Key Takeaways