Bayesian Linear Regression

Compute the posterior over the weights of a Bayesian linear regression.

Ordinary least squares gives one weight vector. The Bayesian version gives a distribution over weight vectors, which also says how sure the model is. Use these assumptions:

• prior: w ~ N(0, (1 / alpha) * I)
• likelihood: y = X @ w + noise, with noise ~ N(0, 1 / beta) for each row

Here alpha is the precision of the prior and beta is the precision of the noise. The posterior is Gaussian, w ~ N(m_N, S_N), with

S_N = inverse(alpha * I + beta * X^T X)
m_N = beta * S_N @ X^T @ y

Return (m_N, S_N).

Function signature

def bayes_linreg_posterior(X, y, alpha, beta):

Example

Input

X = np.array([[1.0], [2.0]])
y = np.array([1.0, 2.0])
bayes_linreg_posterior(X, y, alpha=1.0, beta=1.0)

Output

(array([0.83333333]), array([[0.16666667]]))

Reasoning

X^T X = 1 + 4 = 5, so S_N = 1 / (1 + 5) = 1/6. X^T y = 1 + 4 = 5, so m_N = 5/6. Least squares alone would give w = 1. The prior pulls the mean towards zero, like ridge regression with penalty alpha / beta.

Constraints and notes

• X has shape (n, d) and y has shape (n,).
• m_N has shape (d,) and S_N has shape (d, d).
• The posterior covariance is symmetric.