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):
        Bayesian linear regression model with known error variance and
        conjugate Gaussian prior on model parameters.

        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).

        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

        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}\\ $$


$\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 .


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.