# <span style='font-family:Marker Felt;font-weight:bold'>LINEAR MODELS FOR REGRESSION</span>

---

In this chapter we start to learn some fundamental algorithms behind the most common and simple machine learning tools. To do it, we will use some simple mathematic equations that require basic knowledge of statistics and probability. 

<p align="center">

</p>

## <span style='font-family:Marker Felt;font-weight:bold'>Linear models</span>

---

The goal of a regression problem is to learn a mapping from a given $D$-dimensional vector **x** of input variables to one or more continuous target variables **t**. Some examples of regression problems are the prediction of a stock market price, the prediction of future emissions of $CO_2$, the effect of an actuation in robotics and so on.

The simplest form of regression is given by **linear models** which are models linear in their weights, not necessarily in the inputs. The advantages of these models are many:

* Can be solved analytically. This means, as you will see later, that the associated optimization problem can be solved in _closed form_. We will always obtain the minimum of the approximation loss on the data we have. If the model fails to properly approximate the phenomenon then the problem is for sure related to the hypothesis space, so on the model we are using.

* They are the core concepts of more complex machine learning algorithms.

* Their approximation capability is not limited to linear function: they can approximate any complex function provided the correct input space is considered.

The simplest linear model, is obviously, the one that involves a linear combination of the input variables:

\begin{align}
y(\textbf{x}, \textbf{w}) = w_0 + \sum_{j=1}^{D-1}w_j x_j = \textbf{w}^T\textbf{x}
\end{align}

where $w_0$ is nothing but the offset of the model at zero.

### <span style='font-family:Marker Felt;font-weight:bold'>Loss function</span>

How can we state if this model is sufficiently accurate for our system? We need to define a loss function, also called the _error function_. Generally speaking a loss function can be defined as:

$$
\begin{align}\mathbb{E}[\mathcal{L}] = \int\int L(t, y(\textbf{x}))p(\textbf{x},t)d\textbf{x} dt 
\end{align}
$$

where $p(\textbf{x}, t)$ represents the probability to see that particular combination of input and output; this probability acts as a weight on the loss. Since the distribution of this probability is unknown we have to estimate it with the data points we have.

In regression, a common choice of loss function is the **squared loss function**:

$$
\begin{align}
\mathbb{E}[\mathcal{L}] = \int\int (t - y(\textbf{x}))^2 p(\textbf{x},t)d\textbf{x}dt
\end{align}
$$

Deriving this loss function w.r.t. $y$ and setting the derivative to zero we obtain that the optimal solution is the conditional mean of $t$ conditioned on $\textbf{x}$, also called *regression function*:

$$
\begin{align}
y(\textbf{x}) = \int t p(t|\textbf{x})dt = \mathbb{E}[t|\textbf{x}]
\end{align}
$$

<p align="center">

</p>

There are situations in which the squared loss function doesn't produce any good result, as in the case where $p(t|\textbf{x})$ is multi-modal. The generalization of the squared loss is the *Minkowski loss* defined as:

$$
\begin{align}
\mathbb{E}[\mathcal{L}] = \int\int (t - y(\textbf{x}))^q p(\textbf{x},t)d\textbf{x}dt
\end{align}
$$

Up to now we have only looked at linear functions of weights and inputs. In order to approximate more complex non-linear functions we can use *non-linear basis functions*. In this way the model still remains linear in the weights but it will be non linear in the input features.

$$
\begin{align}
y(\textbf{x}, \textbf{w}) = w_0 + \sum_{j=1}^{M-1}w_j\phi_j(\textbf{x}) = \textbf{w}^T \bf{\phi}(x)
\end{align}
$$

Thanks to $\bf{\phi}$ we transform the original input space into another one. By using nonlinear basis functions, we allow the function $y(\textbf{x}, \textbf{w})$ to be a nonlinear function of the input vector. In this way we can exploit also our domain knowledge by providing a specific transformation of the inputs.

Some examples of basis functions are:

* Polynomial

* Gaussian

* Sigmoidal

* Fourier basis 

The pictures below show the Gaussian basis functions and the Sigmoid ones.

<p align="center">

</p>

<div class="alert alert-block alert-info">
<span style='color:DodgerBlue;font-weight:bold'>REGRESSION APPROACHES</span>

<span style='font-style:italic;font-weight:bold'>Generative</span>: these are the models that explicitly model the distribution of the inputs as well as the outputs. Probabilistically speaking they model also $p(\textbf{x},t)$

<span style='font-style:italic;font-weight:bold'>Discriminative</span>: they directly model the posterior probability $p(t|\textbf{x})$.

<span style='font-style:italic;font-weight:bold'>Direct</span>: they directly find a mapping function from each input to the output
</span>
</div>

## <span style='font-family:Marker Felt;font-weight:bold'>Minimizing Least Squares</span>

---

Let's consider a data set of inputs $\textbf{X}=\{\textbf{x}_1,...,\textbf{x}_N\}$ with the associated target values $t_1,...,t_N$. The error function is defined as the *sum of squared estimate of errors (SSE)*:

$$
\begin{align}
SSE(\textbf{w}) = \mathcal{L}(\textbf{w}) = \frac{1}{2} \sum_{n=1}^N(y(\textbf{x}_n, \textbf{w}) - t_n)^2
\end{align}
$$

The loss can also be written in terms of the $l_2$ norm of the residual errors, known as *residual sum of squares (RSS)*:

$$
\begin{align}
RSS(\textbf{w}) = \left\lVert  \epsilon \right\rVert ^2_2 = \sum^N_{n=1} \epsilon_n^2
\end{align}
$$

Note that if we have an high number $N$ of samples the approximation will be very good otherwise, if $N$ is small, the loss will be very noisy. This is because, without the knowledge of $p(\textbf x, t)$ we have to approximate it with a finite number of points and, higher is the number of points and better is the approximation.

In matrix form the RSS becomes:

$$
\begin{align}
\mathcal{L}(\textbf{w}) = \frac{1}{2}(\textbf{t} - \bf{\Phi}\textbf{w})^T(\textbf{t} - \bf{\Phi}\textbf{w}
\end{align}
$$

where:

$$
\begin{align}
\bf{\Phi} = (\phi(\textbf{x}_1), ...,\phi(\textbf{x}_N))^T  
\end{align}
$$

To obtain the minimum of the loss we calculate its the first and the second derivatives w.r.t. the weights:

$$
\begin{align}
\frac{\partial \mathcal{L}(\textbf{w})}{\partial \textbf{w}} = -\bf{\Phi}^T(\textbf{t} - \bf{\Phi}\textbf{w}) 
\end{align}
$$

$$
\begin{align}
\frac{\partial^2 \mathcal{L}(\textbf{w})}{\partial \textbf{w}\partial \textbf{w}} =  \bf{\Phi}^T\bf{\Phi}
\end{align}
$$

Providing that the Gram matrix (or Gramian matrix, or *Gramian*) $\bf{\Phi}^T\bf{\Phi}$ is nonsingular, the solution of the *ordinary least squares (OLS)* is:

$$
\begin{align}
\hat{\textbf{w}}_{OLS} = (\bf{\Phi}^T\bf{\Phi})^{-1}\bf{\Phi}^T\textbf{t} 
\end{align}
$$

This set of equations are known as *normal equations* for the least squares problem. The inversion of the Gram matrix can be performed with Singular Values Decomposition (SVD) at the complexity cost of $\mathcal{O}(NM^2 + M^3)$, where $M$ is the number of features.

### <span style='font-family:Marker Felt;font-weight:bold'>Sequential learning</span>

Batch techniques can be computationally expensive when dealing with big data. In order to solve this problem we can rely on sequential algorithms, also known as *online algorithms*. This type of updates are also helpful for real-time applications where the data arrives in a continuous stream. Since the total loss can be expressed as the sum of the loss in each data point we can update the parameters with **stochastic gradient descent (SGD)** as follow:

$$
\begin{align}
\textbf{w}^{k+1} = \textbf{w}^{k} - \alpha^{k+1}\nabla\mathcal{L}(x_n) 
\end{align}
$$

where k represents the iteration number and $\alpha$ is called the learning rate. The value of $\alpha$ needs to satisfy certain conditions for convergence:

$$
\begin{align}
\sum^{\infty}_{k=0}\alpha^k = + \infty 
\end{align}
$$

$$
\begin{align}
\sum^{\infty}_{k=0}(\alpha^k)^2 < + \infty
\end{align}
$$

### <span style='font-family:Marker Felt;font-weight:bold'>Geometric interpretation</span>

Let's consider an $N$-dimensional space which axes are given by the $t_n$. Each basis function can also be represented as a vector in the same space. Define $\varphi_j$ the $j^{th}$ column of $\Phi$. If the number $M$ of basis functions is smaller than the number $N$ of data points, then the $\varphi_j$ will span just a linear subspace of dimensionality $M$. Since $\textbf{y}$ is defined as a linear combination of the vectors $\varphi_j$ it lives in the same space of dimension $M$. The least squares solution for $\textbf{w}$ corresponds to the solution $\hat{\textbf t}$ ($\textbf{y}$ in the picture below) that is closest to $\textbf{t}$ in the subspace.

$$
\begin{align}
\hat{\textbf{t}} = \bf{\Phi}\hat{\textbf{w}} = \bf{\Phi}(\bf{\Phi}^T\bf{\Phi})^{-1}\bf{\Phi}^T \textbf{t}
\end{align}
$$

The solution corresponds to the projection of $\textbf{t}$ in the subspace of dimensionality $M$ through a projection matrix called the *hat matrix*:

$$
\begin{align}
H = \bf{\Phi}(\bf{\Phi}^T\bf{\Phi})^{-1}\bf{\Phi}^T
\end{align}
$$

<p align="center">

</p>

### <span style='font-family:Marker Felt;font-weight:bold'>Maximum likelihood</span>

Let's consider the relationship between the least squares approach and the maximum likelihood approach. We assume that the target $t$ is in a deterministic relationship with $\textbf{x}$ and $\textbf{w}$ contaminated with additive noise:

$$
\begin{align}
t = y(\textbf{x},\textbf{w}) + \epsilon 
\end{align}
$$

We assume a Gaussian model for the noise function $\epsilon \sim \mathcal{N}(0, \sigma^2) = \mathcal{N}(0, \beta^{-1}) $

where $\sigma^2$ is the variance and $\beta$ is the precision. Given N samples independent and identically distributed (iid) we can write the *likelihood function* <a name="likelihood"></a>:

$$
\begin{align}
p(\textbf{t}|\textbf{X},\textbf{w},\sigma^2) = \prod^N_{n=1}\mathcal{N}(t_n|\textbf{w}^T\phi(\textbf{x}_n), \sigma^2)
\end{align}
$$
 
<div class="alert alert-block alert-warning">
<b>NOTE:</b> $$ \begin{align} \mathcal{N}(x| \mu, \sigma^2) = \frac{1}{(2 \pi \sigma)^{1/2}} \exp\Big\{-\frac{1}{2\sigma^2}(x- \mu)^2\Big\} \end{align} $$
</div>

Taking the logarithm of the likelihood function:

$$
\begin{align}
\ell(\textbf{w}) = \ln p(\textbf{t}|\textbf{X},\textbf{w},\sigma^2) = \sum^N_{n=1}\ln p(t_n|\textbf{w},\textbf{x}_n, \sigma^2) \\ =-
\frac{N}{2}\log (2 \pi \sigma^2) -\frac{1}{2 \sigma^2} RSS(\textbf{w}) 
\end{align}
$$

To find the maximum likelihood we just have to set to zero the derivative of the likelihood w.r.t. the weights:

$$
\begin{align}
\nabla \ell(\textbf{w}) = \sum^N_{n=1}t_n \phi(\textbf{x}_n)^T - \textbf{w}^T \sum^N_{n=1} \phi(\textbf{x}_n) \phi(\textbf{x}_n)^T = 0
\end{align}
$$

from which we obtain the *maximum likelihood* of the weights:

$$
\begin{align}
\textbf{w}_{ML} = (\bf{\Phi}^T \bf{\Phi})^{-1} \bf{\Phi}^T \textbf{t} 
\end{align}
$$

Notice that the maximum likelihood values of the parameters don't depend on the variance of the additive noise.

Assuming that the observation $t_i$ are uncorrelated and with constant variance $\sigma^2$ and that the $x_i$ are non random, the variance of the least squares estimate is:

$$
\begin{align}
Var(\hat{ \textbf{w} }_{OLS}) = (\bf{\Phi}^T \bf{\Phi})^{-1} \sigma^2
\end{align}
$$

We can notice that the value of the variance decreases as the eigenvalues of $(\bf{\Phi}^T \bf{\Phi})$ increase. The inverse produces low values when the eigenvalues of $\bf{\Phi}^T \bf{\Phi} $ are high so when the features are linearly independent and the number of samples is high.

<div class="alert alert-block alert-warning">
<b>NOTE:</b> $$ \begin{align} A = diag(\lambda_i) \rightarrow A^{-1} = diag(\lambda_i^{-1}) \end{align} $$
</div>

Usually the variance is estimated as:

$$
\begin{align}
\hat{\sigma}^2  = \frac{1}{N-M-1} \sum^N_{n=1}(t_n - \hat{\textbf{w}}^T \phi(\textbf{x}_n))^2
\end{align}
$$

So, we obtain that:

$$
\begin{align}
\hat{\textbf{w}} \sim \mathcal{N}(\textbf{w}, (\bf{\Phi}^T \bf{\Phi})^{-1} \sigma^2)
\end{align}
$$

<div class="alert alert-block alert-info">
<span style='color:DodgerBlue;font-weight:bold'>GAUSS-MARKOV THEOREM</span>

*The least squares estimate of $\textbf{w}$ has the smallest variance among all linear unbiased estimates.*
</div>

This theorem tells us that the least squares estimator has the lowest Mean Squared Error (MSE) among all the linear estimators with no bias. However, there may exists a biased estimator with smaller MSE. *An unbiased estimator is not always the best choice.*

### <span style='font-family:Marker Felt;font-weight:bold'>Multiple outputs</span>

Let's now consider the case of multiple outputs. In this case we have the possibility to use different sets of basis functions for each output, obtaining different and independent regression problems. However, the most common approach is to use the same set of basis functions to model all the components of the target vector.

$$
\begin{align}
\textbf{y(x,w)} = \textbf{W}^T\phi(\textbf{x}) 
\end{align}
$$

where $\textbf{y}$ is a $K$-dimensional column vector, $\textbf{W}$ is an $M \times K$ matrix and $\phi(\textbf{x})$ is an $M$-dimensional column vector.

With the same procedure seen before, the maximum likelihood approach, we obtain the values of the parameters:

$$
\begin{align}
\hat{\textbf{W}}_{ML} = (\bf{\Phi}^T \bf{\Phi})^{-1} \bf{\Phi}^T \textbf{T}
\end{align}
$$

For each output $t_k$ we have:

$$
\begin{align}
\hat{\textbf{w}}_{k} = (\bf{\Phi}^T \bf{\Phi})^{-1} \bf{\Phi}^T \textbf{t}_k = \bf{\Phi}^{\dagger} \textbf{t}_k 
\end{align}
$$

where $\bf{\Phi}^{\dagger}$ is called the *Moore-Penrose pseudo-inverse* and needs to be computed only once.

## <span style='font-family:Marker Felt;font-weight:bold'>Regularization</span>

---

In order to have a clear idea of what is and why we need regularization the best way is considering an example: the polynomial curve fitting. Let's consider the function $sin(2 \pi x)$ and draw 10 samples of this function with random noise included. Now we want to fit these points with a polynomial function of the form:

$$y(x, \textbf{w}) = w_0 + w_1 x + ... + w_M x^M = \sum_{j=0}^M w_j x^j $$

$M$ represents the order of the polynomial. For various values of $M$ we obtain different fitting of the data points:

<p align="center">

</p>
    
The problem is choosing the value of M: this is a very important concept called **model selection**. With a low order polynomial what we obtain is an **underfitting** of the data while, with an high order polynomial, we obtain an **overfitting**. In this last case the model starts to learn the *idiosyncratic* noise of the data set at hand.

* <span style='color:SlateBlue;font-style:italic;font-weight:bold'> Overfitting: </span>the model is too strong and learn the position of the data points inside the training set without discovering the underlying function.

* <span style='color:SlateBlue;font-style:italic;font-weight:bold'> Underfitting: </span>the model is too weak and is not able to approximate well the points inside the training set.

Since the objective of regression, and in general of ML, is the generalization to unseen samples we need a model in the middle of the two extremes. In order to evaluate the overfitting we generate another set composed by $100$ samples, taken from the same distribution of the $10$ samples generate before. To have some quantitative insight into the generalization performance we evaluate the loss both on the initial data set, called *training set*, and also on the new one, called *test set*. In this case we evaluate the Root Mean Square (RMS) error since it allows to account for the difference in size of the two bunch of points.

$$
\begin{align}
E_{RMS} = \sqrt{\frac{2 \cdot RSS}{N}} 
\end{align}
$$

<p align="center">

</p>

The *test set error* is a measure of how well the model will behave on unseen samples. We can see from the above pictures that for $M=9$ the train error goes to zero but the error on the test set increases considerably. 

<p align="center">

</p>

This is because as $M$ increases  the magnitude of the coefficients gets larger, producing wiggling curves. Wiggling curves can be perfectly tuned to pass through all the points of the training set, as in the case of $M=9$, but that's not we want. We know that the points used to optimize the coefficient derive from a function soiled with noise so, we prefer a function that does not interpolate precisely the training points encouraging a smoother interpolation. Below are reported the values of the parameters for various order of the polynomial model.

<p align="center">

</p>

One possibility to overcome overfitting problems is to introduce a regularization term inside the loss function which prevents the coefficients to reach too high values.

### <span style='font-family:Marker Felt;font-weight:bold'>Ridge regression</span>

In ridge regression the loss is modified as follows:

$$
\begin{align}
\mathcal{L}(\textbf{w}) &= \frac{1}{2} \sum_{n=1}^N(y(\textbf{x}_n, \textbf{w}) - t_n)^2 + \frac{\lambda}{2} \left\lVert \textbf{w} \right\rVert_2^2 = \\ &= \mathcal{L}_D(\textbf{w}) + \lambda \mathcal{L}_W(\textbf{w})
\end{align}
$$

where the first term is nothing but the standard RSS while the second term affects the complexity of the model. In ML this choice of regularizer is known as *weights decay*. In statistics this is an example of *parameters shrinkage* since it shrinks the parameter values toward zero. We see that the magnitude of the weights enter directly in the loss function, since we want to minimize the loss we need to converge to a solution with small values of $\textbf w$.

The advantage of this technique is that the loss function remains quadratic on the weights and a closed form solution is possible:

$$
\begin{align}
\textbf w = (\lambda \textbf I + \bf \Phi ^T \bf \Phi)^{-1} \bf \Phi^T \textbf t 
\end{align}
$$

The effect of $\lambda$ is to virtually reduce the hypothesis space. In this case the matrix is always invertible if the regularization term is not null and it gives a lowerbound on the eigenvalues. We can have two cases:

* $\lambda=0$: we have the OLS which is unbiased but has high variance

* $\lambda \ne 0$: we introduce a bias but the variance is reduced

It is important to remark that the bias is necessary when we don't have enough data in order to prevent from overfitting. We can see from the pictures below that as $\lambda$ increases the shape of the approximating function becomes smoother.

<p align="center">

</p>

We can quantify the effect of the regularization by looking at the loss on the train and on the test set.

<p align="center">

</p>

### <span style='font-family:Marker Felt;font-weight:bold'>Lasso regression</span>

Lasso is another popular technique in which the loss is modified as:

$$
\begin{align}
\mathcal{L}(\textbf{w}) = \frac{1}{2} \sum_{n=1}^N(y(\textbf{x}_n, \textbf{w}) - t_n)^2 + \frac{\lambda}{2} \left\lVert \textbf{w} \right\rVert_1 
\end{align}
$$

The important effect of Lasso is that if $\lambda$ is sufficiently large some parameters goes to zero producing a sparse model (easier to interpret) and can be used also for features selection. We can see this by thinking that the Lasso method can be seen as an unregularized minimization problem subjected to the constraint:

$$
\begin{align}
\sum^M_{j=1} \left\lVert w_j \right\rVert < \eta 
\end{align}
$$

The origin of the sparsity of Lasso can be better understand by looking at the following pictures

<p align="center">

</p>

The red circles are the isolines of the unregularized loss function while the light blue parts are the regions of the weights-space that satisfy the constraints.

In the Ridge regression technique the $l_2$ norm squared distance produces an hyper-sphere region (that is a circle in the 2D space) centered in the origin. Reducing $\lambda$ the circle radius reduces as well.
In the Lasso technique the constraint is a polytope (that is a rhomboid in the 2D space). Reducing $\lambda$ the constraint space will reduce as well, but differently from Ridge Regression, $w^*$ can be located exactly on one axis (so that we will have a weight equal to 0).

Regularization allows to train complex models also with few data without obtaining severe overfitting. The problem of determining the maximum number of the polynomial degree is translated into determining the best value of the regularization parameter.

## <span style='font-family:Marker Felt;font-weight:bold'>Bayesian linear regression</span>

---

What we have seen so far was the frequentist approach to linear regression. We have data and we want to find the parameters that better explain the data. In the Bayesian approach the perspective changes. Here the parameters are **uncertain** and they become random variables. We can specify our knowledge on the parameters by specifying a prior distribution based on what we know about the problem. We start from a **prior** distribution, we observe data and this data will change our prior distribution. Combining prior and data we obtain a new probability over parameters called **posterior**. We change our knowledge by observing the real world. The posterior can be used to make predictions and, differently from the frequentist approach where the output is a single value for each weights, here we obtain a probability distribution over all the parameters. In this case we can account for uncertainties and the presence of the uncertainties helps to make decisions about the results.

The posterior can be found by combining the prior with the [likelihood](#likelihood) of the parameters (given data) through the **Bayes' Rule**:

$$
\begin{align}
p(parameters|data) = \frac{p(data|parameters)p(parameters)}{p(data)} 
\end{align}
$$

$$
\begin{align}
p(parameters|data) = \frac{p(data|parameters)p(parameters)}{p(data)} 
\end{align}
$$

$$
\begin{align}
posterior = \frac{likelihood\times prior}{marginal} 
\end{align}
$$

$$
\begin{align}
p(\textbf w | D) = \frac{p( \mathcal D | \textbf w)p(\textbf w)}{p(\mathcal D)} 
\end{align}
$$

where $\textbf w$ are the weights of the model and $ \mathcal D$ represents the observed data.

If we don't have a priori knowledge about parameters distribution we can provide an uninformative prior. This prior is nothing but a uniform probability over all the parameters space. 

Since the denominator is just a normalizing term we can write that the posterior is proportional to product between the likelihood and the prior:

$$
\begin{align}
posterior \propto likelihood \times prior 
\end{align}
$$

If we want to have a single output from our Bayesian model we can take the most probable value of $\textbf w$ given the data, i.e. the mode of the posterior. This approach is called **maximum a posteriori (MAP)** estimate.

The Bayesian approach is another method that allows to avoid overfitting since the prior on the parameters acts as a regularizer. The stronger the prior, the less chance the parameters have of moving in the design space. We can imagine the prior like a spring centered in the mean of the distribution with a stiffness constant $k$ proportional to the variance. The drawback is that if our prior is wrong we add an high bias to the model.

Assuming a multivariate Gaussian prior distribution $ p(\textbf w) = \mathcal{N}( \textbf w | \textbf w_0, \textbf S_0)$ we are specifying that from what we know the mean of the parameters is given by $\textbf w_0$ with a covariance matrix given by $ \textbf S_0$.

If also the likelihood is Gaussian then the posterior will be Gaussian too  (this is true only for Gaussian probability distributions).

$$
\begin{align}
p(\textbf w | \textbf t , \bf \Phi, \sigma^2) \propto \mathcal N (\textbf w | \textbf w_0, \textbf S_0) \mathcal N (\textbf t| \bf \Phi \textbf w, \sigma^2 \textbf I_N) = \mathcal N (\textbf w | \textbf w_N, \textbf S_N) 
\end{align}
$$

$$
\begin{align}
\textbf w_N = \textbf S_N \Big( \textbf S_0^{-1} \textbf w_0 + 
\frac{\bf \Phi^T \textbf t}{\sigma^2} \Big) 
\end{align}
$$

$$
\begin{align}
\textbf S_N^{-1} = \textbf S_0^{-1} + \frac{\bf \Phi^T \bf \Phi}{\sigma^2} 
\end{align}
$$

Let's analyze the results. If we consider the uninformative prior, so $ \textbf S_0 \rightarrow \infty $, and substitute it into the equations before we obtain that the MAP is equal to the maximum likelihood solution of the previous sections. Since the MAP of a Gaussian is equal to the mean:

$$
\begin{align}
\lim_{\textbf S_0 \to \infty} \textbf w_{MAP} = (\bf \Phi^T \bf \Phi)^{-1} \bf \Phi^T \textbf t = \textbf w_{ML} 
\end{align}
$$

If then someone will provide us other data, as in an online sequential case, the just computed posterior will act as a prior for the new model.

What if the mean is $0$? Consider a case in which the prior is defined by a single parameter $\tau^2$:

$$
\begin{align}
p(\textbf w| \tau^2) = \mathcal N (\textbf w | \textbf 0, \tau^2 \textbf I)
\end{align}
$$

we have that the corresponding posterior will have mean and covariance:

$$
\begin{align}
\textbf m_N = \frac{\textbf S_N \bf \Phi^T \textbf t}{\sigma^2}
\end{align}
$$

$$
\begin{align}
\textbf S_N^{-1} = \tau^{-2} \textbf I + \sigma^{-2} \bf \Phi^T \bf \Phi
\end{align}
$$

from which we can compute the log of the posterior:

$$
\begin{align}
\ln p(\textbf w |\textbf t) = -\frac{1}{2 \sigma^2} \sum_{n=1}^N(t_n - \textbf w^T \phi(\textbf x_n))^2 - \frac{1}{2 \tau^2} \left\lVert \textbf w \right\rVert_2^2 
\end{align}
$$

We obtain a result very similar to the case of linear regression with regularization. Also here we are searching for parameters that are around zero and the role of $\lambda$ is played by the covariance matrix of the prior. The maximization of the log of the posterior distribution w.r.t $\textbf w$ is therefore equivalent to the minimization of the sum-of-squares-error with the addition of a quadratic regularization term with $\lambda = \frac{\sigma^2}{\tau^2}$.

Let's consider the example of a single input variable $x$ and a single target variable $t$ given by a linear model of the form $ y(x, \textbf w) = w_0 + w_1 x$. We can generate synthetic data from the function $ y(x, \textbf a) = a_0 + a_1 x$ with for the example $a_0 = -0.3$ and $a_1 = 0.5$ following three steps:

1. Draw $x_n$ from the uniform distribution $U(x|-1,1)$

2. Evaluate $f(x_n, \textbf a)$

3. Add Gaussian noise to obtain $y(x_n, \textbf a)$

At this point our goal is to recover $a_0$ and $a_1$ from the data and we will do that following a Bayesian approach. The first line of the figure below represents the situation before data are obtained. In the central column is plotted the prior over the parameters while in the right column some functions drawn from the prior. In the second line we see the situation after having observed a single point, the blue dot in the right column. Since now we have data we can plot on the left column the likelihood. By multiplying this likelihood with the prior of the previous line, and normalizing, we obtain the posterior distribution in the central column of the second raw. We go on in the same way for all the data. This is also an example of sequential learning. 

<p align="center">

</p>

We can see in the fourth row that as the number of samples increases, the posterior becomes more narrow and in the limit of an infinite number of points it would be a delta function centered on the true value of the parameters.

By taking only the MAP, we are throwing away all the informations about the uncertainty on the weights which are the the very advantage offered by a Bayesian approach. It is better to consider the **posterior predictive distribution**

$$
\begin{align}
p(t| \textbf x, \mathcal D, \sigma^2) = \int \mathcal N(t| \textbf w^T \phi(\textbf x), \sigma^2) \mathcal N (\textbf w| \textbf w_N, \textbf S_N) d\textbf w \\
= \mathcal N(t| \textbf w_N^T \phi(\textbf x), \sigma^2_N(\textbf x)) 
\end{align}
$$

$$
\begin{align}
\sigma^2_N(\textbf x) = \sigma^2 + \phi(\textbf x)^T \textbf S_N \phi(\textbf x) 
\end{align}
$$

The variance term is composed by the irreducible noise inside the data and a term related to the uncertainty about the parameters. The first term is fixed and we cannot act on it but, as $N \rightarrow \infty$ the second term goes to zero. 

By considering again the example with the sinusoidal function we can illustrate the predictive distribution for Bayesian linear regression. For each plot on the left column below, the solid red curve represents the mean of the Gaussian predictive distribution and the shaded region spans one standard deviation either side of the mean. On the right column we have functions drawn from the posterior distribution. As we can see increasing the number of samples, the shaded region decreases and the functions drawn from the posterior become more compact. The shaded region is also narrow near the blue dots since they are observed point and, the only uncertainty we have is the noise on the data.

<p align="center">
 
</p>
<p align="center">

</p>
<p align="center">

</p>
<p align="center">

</p>

When we face a Bayesian linear model we have two challenges:

1. Specify a suitable model, able to describe well the likelihood

2. Provide a suitable prior

In the frequentist approach we just have to face the first problem. The prior is what makes the Bayesian approach strong and at the same time risky. Another aspect related to the Bayesian treatment is associated to the computational issue of the posterior distribution. We have a closed form solution just with Gaussian distributions and so with unimodal problems. When we need to apply other distributions we lose the closed form solution and we need to rely on:

* Analytical integration

* Gaussian (Laplace) approximation

* Monte Carlo integration

* Variational approximation

So, what we have to remember from this last section?

<p align="center">

</p>

<span style='font-family:Marker Felt;font-size:20pt;color:DarkCyan;font-weight:bold'>References:</span>

1. *Restelli M., Machine Learning - course slides*

2. *Bishop C.M., "Pattern Recognition and Machine Learning", chapters: 1.1, 1.2, 1.3, 3.1, 3.3*

e5f5f46c767f22aa8d7490c4d3ba912768c02ca9
ipynb
Jupyter Notebook
lectures/01_Linear_Models_for_Regression.ipynb
stepyt/machine_learning_notes
f7263d96d64a304a7155e0b81b3d2b69df36df24
MIT
lectures/01_Linear_Models_for_Regression.ipynb
stepyt/machine_learning_notes
f7263d96d64a304a7155e0b81b3d2b69df36df24
MIT
lectures/01_Linear_Models_for_Regression.ipynb
stepyt/machine_learning_notes
f7263d96d64a304a7155e0b81b3d2b69df36df24
MIT
Qwen/Qwen-72B
1. YES
2. YES