# Understand the implement detail of BayesianLinearRegression in python

I am learning the implement of BayesianLinearRegression through numpy-ml project

I copy the code here

class BayesianLinearRegressionKnownVariance:
def __init__(self, b_mean=0, b_sigma=1, b_V=None, fit_intercept=True):
r"""
Bayesian linear regression model with known error variance and
conjugate Gaussian prior on model parameters.

Notes
-----
Uses a conjugate Gaussian prior on the model coefficients. The
posterior over model parameters is

.. math::

b \mid b_{mean}, \sigma^2, b_V \sim \mathcal{N}(b_{mean}, \sigma^2 b_V)

Ridge regression is a special case of this model where :math:b_{mean}
= 0, :math:\sigma = 1 and b_V = I (ie., the prior on b is a
zero-mean, unit covariance Gaussian).

Parameters
----------
b_mean : :py:class:ndarray <numpy.ndarray> of shape (M,) or float
The mean of the Gaussian prior on b. If a float, assume b_mean is
np.ones(M) * b_mean. Default is 0.
b_sigma : float
A scaling term for covariance of the Gaussian prior on b. Default
is 1.
b_V : :py:class:ndarray <numpy.ndarray> of shape (N,N) or (N,) or None
A symmetric positive definite matrix that when multiplied
element-wise by b_sigma^2 gives the covariance matrix for the
Gaussian prior on b. If a list, assume b_V = diag(b_V). If None,
assume b_V is the identity matrix. Default is None.
fit_intercept : bool
Whether to fit an intercept term in addition to the coefficients in
b. If True, the estimates for b will have M + 1 dimensions, where
the first dimension corresponds to the intercept. Default is True.
"""
# this is a placeholder until we know the dimensions of X
b_V = 1.0 if b_V is None else b_V

if isinstance(b_V, list):
b_V = np.array(b_V)

if isinstance(b_V, np.ndarray):
if b_V.ndim == 1:
b_V = np.diag(b_V)
elif b_V.ndim == 2:
fstr = "b_V must be symmetric positive definite"
assert is_symmetric_positive_definite(b_V), fstr

self.posterior = {}
self.posterior_predictive = {}

self.b_V = b_V
self.b_mean = b_mean
self.b_sigma = b_sigma
self.fit_intercept = fit_intercept

def fit(self, X, y):
"""
Compute the posterior over model parameters using the data in X and
y.

Parameters
----------
X : :py:class:ndarray <numpy.ndarray> of shape (N, M)
A dataset consisting of N examples, each of dimension M.
y : :py:class:ndarray <numpy.ndarray> of shape (N, K)
The targets for each of the N examples in X, where each target
has dimension K.
"""
# convert X to a design matrix if we're fitting an intercept
if self.fit_intercept:
X = np.c_[np.ones(X.shape[0]), X]

N, M = X.shape
self.X, self.y = X, y

if is_number(self.b_V):
self.b_V *= np.eye(M)

if is_number(self.b_mean):
self.b_mean *= np.ones(M)

b_V = self.b_V
b_mean = self.b_mean
b_sigma = self.b_sigma

b_V_inv = np.linalg.inv(b_V)
L = np.linalg.inv(b_V_inv + X.T @ X)
R = b_V_inv @ b_mean + X.T @ y

# (b_v^{-1} + X^{\top}X)^{-1} @ (b_v^{-1}@b_mean + X^{\top}y)

mu = L @ R
cov = L * b_sigma ** 2

# posterior distribution over b conditioned on b_sigma
self.posterior["b"] = {"dist": "Gaussian", "mu": mu, "cov": cov}


Pull out the fit algorithm :

$$\mu_w^{'} = (b_v^{-1} + X^{\top}X)^{-1} (b_v^{-1}b_{mean} + X^{\top}y) = (b_v^{-1} + X^{\top}X)^{-1}b_v^{-1}b_{mean} + (X^{\top}X+ b_v^{-1} )^{-1}X^{\top}y\\ \\ \Sigma_w^{'} = b_{sigma}^2(X^{\top}X+ b_v^{-1} )^{-1}$$

This is not what I learned .

# Below are my math notes:

We need calculate the postiror $$P(w \mid \operatorname{Data} ) \propto P(Y \mid x, w) \cdot P(w)$$

Let $$P(w \mid Data) = N(\mu_w, \Sigma_w)$$, then we need $$\mu_w , \Sigma_w$$

Since we have $$\begin{array}{l} P(Y|X, w)=\prod_{i=1}^{N} N\left(w^{\top} x_{i}, \sigma^{2}\right) = \frac{1}{\sqrt{2\pi}\sigma^2}e^{-\frac{1}{2\sigma^2} \sum_{i=1}^{N}(y_i - w^{\top}x_i ) } = \frac{1}{\sqrt{2\pi}\sigma^2}e^{-\frac{1}{2}(Y-w^{\top}X)^{\top}\sigma^{-2}I(Y-w^{\top}X) } \\ P(w)=N\left(0, \Sigma_{p}\right) = = \frac{1}{\sqrt{2\pi}\Sigma_p^2}e^{-\frac{1}{2}(w^{\top}\Sigma_p^{-1}w)^2} \end{array}$$

Therefore $$P(w \mid Data) \propto N(Xw,\sigma^{-2}I) \cdot N(0, \Sigma_p)$$

We can derive

$$A = \sigma^{-2}X^{\top}X+\Sigma_p^{-1} \\ \mu_w = \sigma^{-2}A^{-1}X^{\top}Y = (X^{\top}X+\Sigma_p^{-1})^{-1}X^{\top}Y\\ \Sigma_w= A^{-1} = \sigma^{2}(X^{\top}X+\Sigma_p^{-1} )^{-1}\\$$

Question:

$$\Sigma_w^{'}$$ is much like $$\Sigma_w$$, $$b_{sigma}^2$$ is variance of error, $$\Sigma_p$$ is assigned by us . I can understand here .

But the $$\mu_w$$ is quite different from the implement of numpy-ml , how do I understand it ?

$$\mu_w^{'} = (b_v^{-1} + X^{\top}X)^{-1}b_v^{-1}b_{mean} + (X^{\top}X+ b_v^{-1} )^{-1}X^{\top}y\\ \mu_w = (X^{\top}X+\Sigma_p^{-1})^{-1}X^{\top}Y$$

Here $$(b_v^{-1} + X^{\top}X)^{-1}b_v^{-1}b_{mean}$$ looks like a regulizasition term, but refer to the code docurment, it says

Ridge regression is a special case of this model where :math:b_{mean} = 0, :math:\sigma = 1 and b_V = I

As we know Ridge has L2 term, but $$b_{mean}=0$$ would exterminate $$(b_v^{-1} + X^{\top}X)^{-1}b_v^{-1}b_{mean}$$, no term ? I can't figure it out .