Bayesian regression

Notebook version: 1.0 (Oct 01, 2015)

Author: Jerónimo Arenas García (jarenas@tsc.uc3m.es)

Changes: v.1.0 - First version

Pending changes: * Include regression on the stock data

3. Bayesian regression

In the previous session we tackled the problem of fitting the following model using a LS criterion:

$${\hat s}({\bf x}) = f({\bf x}) = {\bf w}^\top {\bf z}$$

where ${\bf z}$ is a vector with components which can be computed directly from the observed variables. Such model is includes a linear regression problem, where ${\bf z} = [1; {\bf x}]$, as well as any other non-linear model as long as it can be expressed as a "linear in the parameters" model.

The LS solution was defined as the one minimizing the square of the residuals over the training set $\{{\bf x}^{(k)}, s^{(k)}\}_{k=1}^K$. As a result, a single parameter vector ${\bf w}^*$ was obtained, and correspondingly a single regression curve.

In this session, rather than trying to obtain the best single model, we will work with a family of models or functions, and model the problem probabilistically, so that we can assign a probability value to each of the possible functions.

3.1 Maximum Likelihood estimation of the weights

3.1.1 Limitations of the LS approach. The need for assumptions

Consider the same regression task of the previous session. We have a training dataset consisting of 15 points which are given, and depict the regression curves that would be obtained if adding an additional point at a fixed location, depending on the target value of that point:

(You can run this code fragment several times, to check also the changes in the regression curves between executions, and depending also on the location of the training points)

• You can control the degree of the polynomia, and check that when the degree is set to 15 (16 weights) all points will be fitted perfectly
• It seems obvious that we have not solved the problem ...
  • The regression curves overfit the training data
  • The regression curves change a lot when varying the label of just one pattern

The key missing ingredient is assumptions !!

Open questions

• Do we think that all models are equally probable... before we see any data? What does the term model probability mean?

• Do we need to choose a single "best" model or can we consider several simultaneously?

• Perhaps our training targets are contaminated with noise. What to do?

We will start postulating a generative model for the training data that includes the presence of noise contaminating the targets, and work on this model to partly answer the other two questions.

 3.1.2 Generative model

Denoting by $f({\bf x}) = {{\bf w}}^\top {\bf z}$ the true function that we would like to obtain, we could assume that the observations in the training set are obtained as noisy values of the output of such function, i.e.,

$$s^{(k)} = f({\bf x}^{(k)}) + \varepsilon^{(k)}$$

We will further characterize the noise values as i.i.d. and normally distributed, with mean zero, and variance $\sigma_\varepsilon^2$, i.e.,

$$\varepsilon \sim {\cal N}\left(0, \sigma_\varepsilon^2\right)$$

3.1.3 The maximum likelihood solution

• Joint distribution of the noise samples, ${\pmb \varepsilon} = \left[\varepsilon^{(1)}, \dots, \varepsilon^{(K)}\right]^\top$:

  $${\pmb \varepsilon} \sim \left( {\bf 0}, \sigma_{\varepsilon}^2 {\bf I}\right) \;\;\; p({\pmb \varepsilon}) = \left(\frac{1}{\sqrt{2\pi \sigma_{\varepsilon}^2}}\right)^K \exp\left(- \frac{{\pmb \varepsilon}^\top {\pmb \varepsilon}}{2 \sigma_{\varepsilon}^2}\right)$$

• Denoting ${\bf s} = \left[s^{(1)}, \dots, s^{(K)} \right]^\top$ and ${\bf f} = \left[ f({\bf x}^{(1)}), \dots,f({\bf x}^{(K)})\right]^\top$, we have

  $${\bf s} = {\bf f} + {\pmb \varepsilon}$$

• Conditioning on the values of the target function, ${\bf f}$, the pdf of the available targets is obtained as a shifted version of the distribution of the noise. More precisely:

  \begin{align}p({\bf s}|{\bf f}) & = \left(\frac{1}{\sqrt{2\pi \sigma_{\varepsilon}^2}}\right)^K \exp\left(- \frac{\|{\bf s} - {\bf f}\|^2}{2 \sigma_{\varepsilon}^2}\right) \end{align}

• For the particular parametric selection of $f({\bf x})$, ${\bf f} = {\bf Z} {\bf w}$, conditioning on ${\bf f}$ is equivalent to conditioning on ${\bf w}$, so that:

  $$p({\bf s}|{\bf f}) = p({\bf s}|{\bf w}) = \left(\frac{1}{\sqrt{2\pi \sigma_{\varepsilon}^2}}\right)^K \exp\left(- \frac{\|{\bf s} - {\bf Z}{\bf w}\|^2}{2 \sigma_{\varepsilon}^2}\right)$$

• The previous expression represents the probability of the observed targets given the weights, and is also known as the likelihood of the weights for a particular training set.

The maximum likelihood solution is then given by:

$${\bf w}_{ML} = \arg \max_{\bf w} p({\bf s}|{\bf w}) = \arg \min_{\bf w} \|{\bf s} - {\bf Z}{\bf w}\|^2$$

3.1.4 Multiple explanations of the data

With an additive Gaussian independent noise model, the maximum likelihood and the least squares solutions are the same. We have not improved much ...

However, we have already formulated the problem in a probabilistic way. This opens the door to reasoning in terms of a set of possible explanations, not just one. We believe more than one of our models could have generated the data.

• We do not believe all models are equally likely to have generated the data
• We may believe that a simpler model is more likely than a complex one

3.2 Bayesian Inference

3.2.1 Posterior distribution of weights

• If we express our a priori belief of models using a prior distribution $p({\bf f})$, then we can infer the a posteriori distribution using Bayes' rule:

  $$p({\bf f}|{\bf s}) = \frac{p({\bf s}|{\bf f})~p({\bf f})}{p({\bf s})}$$

  In the previous expression:

  • $p({\bf s}|{\bf f})$: is the likelihood function
  • $p({\bf f})$: is the prior distribution of the models (assumptions are needed here)
  • $p({\bf s})$: is the marginal distribution of the observed data, which could be obtained integrating the numerator over all possible models. However, we normally do not need to explicitly compute $p({\bf s})$

• For the parametric model ${\bf f} = {\bf Z} {\bf w}$, the previous expressions become:

  $$p({\bf w}|{\bf s}) = \frac{p({\bf s}|{\bf w})~p({\bf w})}{p({\bf s})}$$

  Where:

  • $p({\bf s}|{\bf w})$: is the likelihood function
  • $p({\bf w})$: is the prior distribution of the weights (assumptions are needed here)
  • $p({\bf s})$: is the marginal distribution of the observed data, which could be obtained integrating

3.2.2 Maximum likelihood vs Bayesian Inference. Making predictions

• Following a ML approach, we retain a single model, ${\bf w}_{ML} = \arg \max_{\bf w} p({\bf s}|{\bf w})$. Then, the predictive distribution of the target value for a new point would be obtained as:

  $$p({s^*}|{\bf w}_{ML},{\bf x}^*) $$

  For the generative model of Section 3.1.2 (additive i.i.d. Gaussian noise), this distribution is:

  $$p({s^*}|{\bf w}_{ML},{\bf x}^*) = \frac{1}{\sqrt{2\pi\sigma_\varepsilon^2}} \exp \left(-\frac{\left(s^* - {\bf w}_{ML}^\top {\bf z}^*\right)^2}{2 \sigma_\varepsilon^2} \right)$$

  • The mean of $s^*$ is just the same as the prediction of the LS model, and the same uncertainty is assumed independently of the observation vector (i.e., the variance of the noise of the model).

  • If a single value is to be kept, we would probably keep the mean of the distribution, which is equivalent to the LS prediction.

• Using Bayesian inference, we retain all models. Then, the inference of the value $s^* = s({\bf x}^*)$ is carried out by mixing all models, according to the weights given by the posterior distribution.

  \begin{align}p({s^}|{\bf x}^,{\bf s}) & = \int p({s^},{\bf w}~|~{\bf x}^,{\bf s}) d{\bf w} \

                  & = \int p({s^*}~|~{\bf w},{\bf x}^*,{\bf s}) p({\bf w}~|~{\bf x}^*,{\bf s}) d{\bf w} \\
                  & = \int p({s^*}~|~{\bf w},{\bf x}^*) p({\bf w}~|~{\bf s}) d{\bf w}\end{align}
  
  

  where:

  • $p({s^*}|{\bf w},{\bf x}^*) = \displaystyle\frac{1}{\sqrt{2\pi\sigma_\varepsilon^2}} \exp \left(-\frac{\left(s^* - {\bf w}^\top {\bf z}^*\right)^2}{2 \sigma_\varepsilon^2} \right)$
  • $p({\bf w}~|~{\bf s})$: Is the posterior distribution of the weights, that can be computed using Bayes' Theorem.

 3.2.3 Example: Selecting a Gaussian prior for the weights

 Prior distribution of the weights

In this section, we consider a particular example in which we assume the following prior for the weights:

$${\bf w} \sim {\cal N}\left({\bf 0},{\pmb \Sigma}_{p} \right)$$

The following figure shows functions which are generated by drawing points from this distribution

Likelihood of the weights

According to the generative model, ${\bf s} = {\bf Z}{\bf w} + {\pmb \varepsilon}$

$${\bf s}~|~{\bf w} \sim {\cal N}\left({\bf Z}{\bf w},\sigma_\varepsilon^2 {\bf I} \right)$$

Posterior distribution of the weights

$$p({\bf w}|{\bf s}) = \frac{p({\bf s}|{\bf w})~p({\bf w})}{p({\bf s})}$$

Since both $p({\bf s}|{\bf w})$ and $p({\bf w})$ follow a Gaussian distribution, we know also that the joint distribution and the posterior distribution of ${\bf w}$ given ${\bf s}$ are also Gaussian. Therefore,

$${\bf w}~|~{\bf s} \sim {\cal N}\left({\bar{\bf w}} , {\pmb\Sigma}_{\bf w}\right)$$

where the mean and the covariance matrix of the distribution are to be determined.

Exercise:

Show that the posterior mean and posterior covariance matrix of ${\bf w}$ given ${\bf s}$ are:

$${\pmb\Sigma}_{\bf w} = \left[\frac{1}{\sigma_\varepsilon^2} {\bf Z}^{\top}{\bf Z} + {\pmb \Sigma}_p^{-1}\right]^{-1}$$$${\bar{\bf w}} = {\sigma_\varepsilon^{-2}} {\pmb\Sigma}_{\bf w} {\bf Z}^\top {\bf s}$$

The following fragment of code draws random vectors from $p({\bf w}|{\bf s})$, and plots the corresponding regression curve along with the training points. Compare these curves with those extracted from the prior distribution of ${\bf w}$.

Posterior distribution of the target

• Since $f^* = f({\bf x}^*) = [{\bf x}^*]^\top {bf w}$, $f^*$ is also a Gaussian variable whose posterior mean and variance can be calculated as follows:

  $$\mathbb{E}\{{{\bf x}^*}^\top {\bf w}~|~{\bf s}, {\bf x}^*\} = {{\bf x}^*}^\top \mathbb{E}\{{\bf w}|{\bf s}\} = {\sigma_\varepsilon^{-2}} {{\bf x}^*}^\top {\pmb\Sigma}_{\bf w} {\bf Z}^\top {\bf s}$$

  $$\text{Cov}\left[{{\bf x}^*}^\top {\bf w}~|~{\bf s}, {\bf x}^*\right] = {{\bf x}^*}^\top \text{Cov}\left[{\bf w}~|~{\bf s}\right] {{\bf x}^*} = {{\bf x}^*}^\top {\pmb \Sigma}_{\bf w} {{\bf x}^*}$$

• Therefore, $f^*~|~{\bf s}, {\bf x}^* \sim {\cal N}\left({\sigma_\varepsilon^{-2}} {{\bf x}^*}^\top {\pmb\Sigma}_{\bf w} {\bf Z}^\top {\bf s}, {{\bf x}^*}^\top {\pmb \Sigma}_{\bf w} {{\bf x}^*} \right)$

• Finally, for $s^* = f^* + \varepsilon^*$, the posterior distribution is $s^*~|~{\bf s}, {\bf x}^* \sim {\cal N}\left({\sigma_\varepsilon^{-2}} {{\bf x}^*}^\top {\pmb\Sigma}_{\bf w} {\bf Z}^\top {\bf s}, {{\bf x}^*}^\top {\pmb \Sigma}_{\bf w} {{\bf x}^*} + \sigma_\varepsilon^2\right)$

Not only do we obtain a better predictive model, but we also have confidence intervals (error bars) for the predictions.