# Regression

### Preliminaries

- Goal 
  - Introduction to Bayesian (Linear) Regression
- Materials        
  - Mandatory
    - These lecture notes
  - Optional
    - Bishop pp. 152-158 
    

### Regression - Illustration


<p style="text-align:center;"></p>

Given a set of (noisy) data measurements, find the 'best' relation between an input variable $x \in \mathbb{R}^D$ and input-dependent outcomes $y \in \mathbb{R}$


### Regression vs Density Estimation


- Observe $N$ IID data **pairs** $D=\{(x_1,y_1),\dotsc,(x_N,y_N)\}$ with $x_n \in \mathbb{R}^M$ and $y_n \in \mathbb{R}$.

- Assume that we are interested in (a model for) the responses $y_n$ for **given inputs** $x_n$?, I.o.w. we are interested in building a model for the conditional distribution $p(y|x)$.

- Note that, since $p(x,y)=p(y|x)\, p(x)$, building a model $p(y|x)$ is similar to density estimation with the assumption that $x$ is drawn from a uniform distribution.

### Bayesian Linear Regression

- Next, we discuss (1) model specification, (2) Inference and (3) a prediction application for a Bayesian linear regression problem. 

#####  1. Model Specification


- In a _regression_ model, we try to 'explain the data' by a purely deterministic term $f(x_n,w)$, plus a purely random term $\epsilon_n$ for 'unexplained noise',

    $$
    y_n  = f(x_n,w) + \epsilon_n
    $$

-  In a _linear regression_ model, i.e., linear w.r.t. the parameters $w$, we assume that 
$$f(x_n,w)= \sum_{j=0}^{M-1} w_j \phi_j(x_n) = w^T \phi(x_n)$$
where $\phi_j(x)$ are called basis functions.
  - For notational simplicity, from now on we will assume $f(x_n,w) = w^T x_n$, with $x_n \in \mathbb{R}^M$.

-  In _ordinary linear regression_ , the noise process $\epsilon_n$ is zero-mean Gaussian with constant variance, i.e.
$$
y_n  = w^T x_n  + \mathcal{N}(0,\beta^{-1}) \,.
$$

- Hence, given a data set  $D=\{(x_1,y_1),\dotsc,(x_N,y_N)\}$, the likelihood for an ordinary linear regression model is 

$$\begin{align*}
p(y\,|\,\mathbf{X},w,\beta) &= \mathcal{N}(y\,|\,\mathbf{X} w,\beta^{-1} \mathbf{I}) \\
  &= \prod_n \mathcal{N}(y_n\,|\,w^T x_n,\beta^{-1}) \tag{B-3.10}
\end{align*}$$

where $w = \left(\begin{matrix} w_1 \\ w_2 \\ \vdots \\ w_{M} \end{matrix} \right)$, the $(N\times M)$-dim matrix $\mathbf{X}  = \left(\begin{matrix}x_1^T \\ x_2^T \\ \vdots \\ x_N^T \end{matrix} \right) = \left(\begin{matrix}x_{11},x_{12},\dots,x_{1M}\\ x_{21},x_{22},\dots,x_{2M} \\ \vdots \\ x_{N1},x_{N2},\dots,x_{NM} \end{matrix} \right) $  and $y = \left(\begin{matrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{matrix} \right)$.


- For full Bayesian learning we should also choose a prior $p(w)$:

$$\begin{equation}
p(w\,|\,\alpha) = \mathcal{N}(w\,|\,0,\alpha^{-1}\mathbf{I}) \tag{B-3.52}
\end{equation}$$

- For simplicity, we will assume that $\alpha$ and $\beta$ are known. 


##### 2. Inference

- We'll do Bayesian inference for the parameters $w$. 

$$\begin{align*}
p(w|D) &\propto p(D|w)\cdot p(w) \\
   &= \mathcal{N}(y\,|\,\mathbf{X} w,\beta^{-1} \mathbf{I}) \cdot \mathcal{N}(w\,|\,0,\alpha^{-1}\mathbf{I}) \\
   &\propto \exp \big( -\frac{\beta}{2} \big( {y - \mathbf{X}w } \big)^T \big( {y - \mathbf{X}w } \big)  - \frac{\alpha}{2}w^T w \big)\tag{B-3.55} \\
   &= \exp\big( -\frac{1}{2} w^T\big(\underbrace{\beta \mathbf{X}^T\mathbf{X} + \alpha \mathbf{I}}_{\Lambda_N}\big)w + \big(\underbrace{\beta \mathbf{X}^Ty}_{\eta_N}\big)^T w - \frac{\beta}{2}y^T y \big) \\
     &\propto \mathcal{N}_c\left(w\,|\,\eta_N,\Lambda_N \right)
\end{align*}$$

with natural parameters (see the [natural parameterization of Gaussian](https://nbviewer.jupyter.org/github/bertdv/BMLIP/blob/master/lessons/notebooks/The-Gaussian-Distribution.ipynb#natural-parameterization)):

$$\begin{align*}
\eta_N &= \beta\mathbf{X}^Ty \\
\Lambda_N &= \beta \mathbf{X}^T\mathbf{X} + \alpha \mathbf{I}
\end{align*}$$

- Or equivalently (in the moment parameterization of the Gaussian):
$$\begin{align*}
p(w|D) &= \mathcal{N}\left(w\,|\,m_N,S_N \right) \tag{B-3.49} \\
m_N &= \beta S_N \mathbf{X}^T y  \tag{B-3.53}\\
S_N &= \left(\alpha \mathbf{I} + \beta \mathbf{X}^T \mathbf{X}\right)^{-1} \tag{B-3.54}
\end{align*}$$

- Note that B-3.53 and B-3.54 combine to
$$
m_N = \left(\frac{\alpha}{\beta}\mathbf{I} + \mathbf{X}^T \mathbf{X} \right)^{-1} \mathbf{X}^T y\,.
$$


 - (Bishop Fig.3.7) Illustration of sequential Bayesian learning for a simple linear model of the form $y(x, w) =
w_0 + w_1 x$. (Bishop Fig.3.7, detailed description at Bishop, pg.154.)

<p style="text-align:center;"> </p>

##### 3. Application: predictive distribution


- Assume we are interested in the distribution $p(y_\bullet \,|\, x_\bullet, D)$ for a new input $x_\bullet$. This can be worked out as (exercise B-3.10)

$$\begin{align*}
p(y_\bullet \,|\, x_\bullet, D) &= \int p(y_\bullet \,|\, x_\bullet, w) p(w\,|\,D)\,\mathrm{d}w \\
&= \int \mathcal{N}(y_\bullet \,|\, w^T x_\bullet, \beta^{-1}) \mathcal{N}(w\,|\,m_N,S_N)\,\mathrm{d}w \\
&= \int \mathcal{N}(y_\bullet \,|\, z, \beta^{-1}) \mathcal{N}(z\,|\,x_\bullet^T m_N,x_\bullet^T S_N x_\bullet)\,\mathrm{d}z \quad \text{(sub. }z=x_\bullet^T w \text{)} \\
&= \int \underbrace{\mathcal{N}(z \,|\, y_\bullet, \beta^{-1}) \mathcal{N}(z\,|\,x_\bullet^T m_N,x_\bullet^T S_N x_\bullet)}_{\text{Use Gaussian product formula SRG-6}}\,\mathrm{d}z  \\
&= \int \mathcal{N}\left(y_\bullet\,|\, m_N^T x_\bullet, \sigma_N^2(x_\bullet) \right) \underbrace{\mathcal{N}\left(z\,|\, m_z, S_z\right)}_{\text{integrate this out}} \mathrm{d}z\\
&= \mathcal{N}\left(y_\bullet\,|\, m_N^T x_\bullet, \sigma_N^2(x_\bullet) \right)
\end{align*}$$
with
$$\begin{align*}
\sigma_N^2(x_\bullet) = \beta^{-1} + x^T_\bullet S_N x_\bullet \tag{B-3.59}
\end{align*}$$

- So, the uncertainty $\sigma_N^2(x_\bullet)$ about the output contains comprises both uncertainty about the process ($\beta^{-1}$) and about the model parameters ($x^T_\bullet S_N x_\bullet$).  

### Example Predictive Distribution

- As an example, let's do Bayesian Linear Regression for a synthetic sinusoidal data set and a model with 9 Gaussian basis functions 
$$\begin{align*}
y &=\sum_{m=1}^9 w_m \phi_m(x_n) + \mathcal{N}(0,\beta^{-1}) \\
\phi_m(x_n) &= \exp\left( \frac{x_n-\mu_m}{\sigma^2}\right)
\end{align*}$$

<p style="text-align:center;"></p>

- The predictive distributions for $y$ are shown in the following plots (Bishop, Fig.3.8)
<p style="text-align:center;"> </p>


- And some plots of draws of posteriors for the functions $w^T \phi(x)$ (Bishop, Fig.3.9) 

<p style="text-align:center;"></p> 

###  Maximum Likelihood Estimation for Linear Regression Model

- Recall the posterior mean for the weight vector
$$
m_N = \left(\frac{\alpha}{\beta}\mathbf{I} + \mathbf{X}^T \mathbf{X} \right)^{-1} \mathbf{X}^T y
$$
where $\alpha$ is the prior precision for the weights.

- The Maximum Likelihood solution for $w$ is obtained by letting $\alpha \rightarrow 0$, which leads to 

$$\begin{equation*}
\hat w_{\text{ML}} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T y
\end{equation*}$$

-  The matrix $\mathbf{X}^\dagger \equiv  (\mathbf{X}^T \mathbf{X})^{-1}\mathbf{X}^T$ is also known as the **Moore-Penrose pseudo-inverse** (which is sort-of-an-inverse for non-square matrices).

-  Note that if we have fewer training samples than input dimensions, i.e., if $N<M$, then $\mathbf{X}^T \mathbf{X}$  will not be invertible and maximum likelihood blows up. The Bayesian solution does not suffer from this problem.

###  Least-Squares Regression

-  (You may say that) we don't need to work with probabilistic models. E.g., there's also the deterministic **least-squares** solution: minimize sum of squared errors,

$$\begin{align*} \hat w_{\text{LS}} &= \arg\min_{w} \sum_n {\left( {y_n  - w ^T x_n } \right)} ^2 
  = \arg\min_{w} \left( {y - \mathbf{X}w } \right)^T \left( {y - \mathbf{X} w } \right)
\end{align*}$$

-  Setting the <a id="regression-gradient">gradient</a> 
$ \frac{\partial \left( {y - \mathbf{X}w } \right)^T \left( {y - \mathbf{X}w } \right)}{\partial w} = -2 \mathbf{X}^T \left(y - \mathbf{X} w  \right)
$ to zero yields the so-called  **normal equations** 
$\mathbf{X}^T\mathbf{X} \hat w_{\text{LS}} = \mathbf{X}^T y$  and consequently

$$
\hat w_{\text{LS}} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T y
$$

which is the same answer as we got for the maximum likelihood weights $\hat w_{\text{ML}}$.

- $\Rightarrow$ Least-squares regression ($\hat w_{\text{LS}}$) corresponds to the (probabilistic) maximum likelihood solution ($\hat w_{\text{ML}}$) if the probabilistic model includes the following assumptions:
  1. The observations are independently and identically distributed (**IID**) (this determines how errors are combined), and
  1. The noise signal $\epsilon_n \sim \mathcal{N}(0,\,\beta^{-1})$ is **Gaussian** distributed (determines the error metric)      
  

- If you use the Least-Squares method, you cannot see (nor modify) these assumptions. The probabilistic method forces you to state all assumptions explicitly! 

###  Not Identically Distributed Data

-  Let's do an example regarding changing our assumptions. What if we assume that the variance of the measurement error varies with the sampling index,  $\epsilon_n \sim \mathcal{N}(0,\beta_n^{-1})$?

- The likelihood is now (using $\Lambda \triangleq \mathrm{diag}(\beta_n)$ )
$$
p(y\,|\,\mathbf{X},w,\Lambda) = \mathcal{N}(y\,|\,\mathbf{X} w,\Lambda^{-1} ) \,.
$$

- Combining this likelihood with the prior $p(w) = \mathcal{N}(w\,|\,0,\alpha^{-1}\mathbf{I})$ leads to a posterior

$$\begin{align*}
p(w|D) &\propto p(D|w)\cdot p(w) \\
   &= \mathcal{N}(y\,|\,\mathbf{X} w,\Lambda^{-1} \mathbf{I}) \cdot \mathcal{N}(w\,|\,0,\alpha^{-1}\mathbf{I}) \\
   &\propto \exp \left\{ \frac{1}{2} \left( {y - \mathbf{X}w } \right)^T \Lambda \left( {y - \mathbf{X}w } \right)  + \frac{\alpha}{2}w^T w \right\} \\
   &\propto \mathcal{N}\left(w\,|\,m_N,S_N \right) 
\end{align*}$$
with
$$\begin{align*}
m_N &=  S_N \mathbf{X}^T \Lambda y  \\
S_N &= \left(\alpha \mathbf{I} +  \mathbf{X}^T \Lambda \mathbf{X}\right)^{-1} 
\end{align*}$$

- And maximum likelihood solution 
$$
\hat{w}_{\text{ML}} = \left. m_N\right|_{\alpha \rightarrow 0} = \left(\mathbf{X}^T \Lambda \mathbf{X}\right)^{-1} \mathbf{X}^T \Lambda y
$$

-  This maximum likelihood solution is also called the **Weighted Least Squares** (WLS) solution. (Note that we just stumbled upon it, the crucial aspect is appropriate model specification!)

-  Note also that the dimension of $\Lambda$ grows with the number of data points. In general, models for which the number of parameters grow as the number of observations increase are called **non-parametric models**.

#### CODE EXAMPLE

We'll compare the Least Squares and Weighted Least Squares solutions for a simple linear regression model with input-dependent noise:

$$\begin{align*}
x &\sim \text{Unif}[0,1]\\
y|x &\sim \mathcal{N}(f(x), v(x))\\
f(x) &= 5x - 2\\
v(x) &= 10e^{2x^2}-9.5\\
\mathcal{D} &= \{(x_1,y_1),\ldots,(x_N,y_N)\}
\end{align*}$$


```julia
using Pkg;Pkg.activate("probprog/workspace/");Pkg.instantiate();
IJulia.clear_output();
```


```julia
using PyPlot, LinearAlgebra

# Model specification: y|x ~ 𝒩(f(x), v(x))
f(x) = 5*x .- 2 
v(x) = 10*exp.(2*x.^2) .- 9.5 # input dependent noise variance
x_test = [0.0, 1.0]
plot(x_test, f(x_test), "k--") # plot f(x)

# Generate N samples (x,y), where x ~ Unif[0,1]
N = 50
x = rand(N)
y = f(x) + sqrt.(v(x)) .* randn(N)
plot(x, y, "kx"); xlabel("x"); ylabel("y") # Plot samples

# Add constant to input so we can estimate both the offset and the slope
_x = [x ones(N)]
_x_test = hcat(x_test, ones(2))

# LS regression
w_ls = pinv(_x) * y
plot(x_test, _x_test*w_ls, "b-") # plot LS solution

# Weighted LS regression
W = Diagonal(1 ./ v(x)) # weight matrix
w_wls = inv(_x'*W*_x) * _x' * W * y
plot(x_test, _x_test*w_wls, "r-") # plot WLS solution
ylim([-5,8]); legend(["f(x)", "D", "LS linear regr.", "WLS linear regr."],loc=2);
```


```julia
open("../../styles/aipstyle.html") do f
    display("text/html", read(f, String))
end
```


<!--
This HTML file contains custom styles and some javascript.
Include it a Jupyter notebook for improved rendering.
-->

<!-- Fonts -->
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'>

<!-- Custom style -->
<style>

@font-face {
    font-family: "Computer Modern";
    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}

#notebook_panel { /* main background */
    background: rgb(245,245,245);
}

div.container {
    min-width: 960px;
}

div #notebook { /* centre the content */
    background: #fff; /* white background for content */
    margin: auto;
    padding-left: 0em;
}

#notebook li { /* More space between bullet points */
    margin-top:0.8em;
}

/* draw border around running cells */
div.cell.border-box-sizing.code_cell.running {
    border: 1px solid #111;
}

/* Put a solid color box around each cell and its output, visually linking them*/
div.cell.code_cell {
    background-color: rgb(256,256,256);
    border-radius: 0px;
    padding: 0.5em;
    margin-left:1em;
    margin-top: 1em;
}

div.text_cell_render{
    font-family: 'Alegreya Sans' sans-serif;
    line-height: 140%;
    font-size: 125%;
    font-weight: 400;
    width:800px;
    margin-left:auto;
    margin-right:auto;
}


/* Formatting for header cells */
.text_cell_render h1 {
    font-family: 'Nixie One', serif;
    font-style:regular;
    font-weight: 400;
    font-size: 45pt;
    line-height: 100%;
    color: rgb(0,51,102);
    margin-bottom: 0.5em;
    margin-top: 0.5em;
    display: block;
}

.text_cell_render h2 {
    font-family: 'Nixie One', serif;
    font-weight: 400;
    font-size: 30pt;
    line-height: 100%;
    color: rgb(0,51,102);
    margin-bottom: 0.1em;
    margin-top: 0.3em;
    display: block;
}

.text_cell_render h3 {
    font-family: 'Nixie One', serif;
    margin-top:16px;
    font-size: 22pt;
    font-weight: 600;
    margin-bottom: 3px;
    font-style: regular;
    color: rgb(102,102,0);
}

.text_cell_render h4 {    /*Use this for captions*/
    font-family: 'Nixie One', serif;
    font-size: 14pt;
    text-align: center;
    margin-top: 0em;
    margin-bottom: 2em;
    font-style: regular;
}

.text_cell_render h5 {  /*Use this for small titles*/
    font-family: 'Nixie One', sans-serif;
    font-weight: 400;
    font-size: 16pt;
    color: rgb(163,0,0);
    font-style: italic;
    margin-bottom: .1em;
    margin-top: 0.8em;
    display: block;
}

.text_cell_render h6 { /*use this for copyright note*/
    font-family: 'PT Mono', sans-serif;
    font-weight: 300;
    font-size: 9pt;
    line-height: 100%;
    color: grey;
    margin-bottom: 1px;
    margin-top: 1px;
}

.CodeMirror{
    font-family: "PT Mono";
    font-size: 90%;
}

.boxed { /* draw a border around a piece of text */
  border: 1px solid blue ;
}

h4#CODE-EXAMPLE,
h4#END-OF-CODE-EXAMPLE {
    margin: 10px 0;
    padding: 10px;
    background-color: #d0f9ca !important;
    border-top: #849f81 1px solid;
    border-bottom: #849f81 1px solid;
}

.emphasis {
    color: red;
}

.exercise {
    color: green;
}

.proof {
    color: blue;
}

code {
  padding: 2px 4px !important;
  font-size: 90% !important;
  color: #222 !important;
  background-color: #efefef !important;
  border-radius: 2px !important;
}

/* This removes the actual style cells from the notebooks, but no in print mode
   as they will be removed through some other method */
@media not print {
  .cell:nth-last-child(-n+2) {
    display: none;
  }
}

footer.hidden-print {
    display: none !important;
}

</style>

<!-- MathJax styling -->





```julia

```

fd70e79a9c2f5960dd8dc0d3cf864967a39c9d6b
ipynb
Jupyter Notebook
lessons/notebooks/Regression.ipynb
bertdv/BMLIP
58920c0117efb019f2256c9ef3ca6ed65db69544
CC-BY-3.0
2019-09-14T17:34:14.000Z
2022-01-22T18:29:11.000Z
lessons/notebooks/Regression.ipynb
charlielu05/BMLIP
70c4c3810e0fdea42d611b6c4aab9003506dc243
CC-BY-3.0
2019-08-09T15:49:10.000Z
2021-11-14T10:19:28.000Z
lessons/notebooks/Regression.ipynb
charlielu05/BMLIP
70c4c3810e0fdea42d611b6c4aab9003506dc243
CC-BY-3.0
2020-03-18T14:05:09.000Z
2022-01-04T14:35:32.000Z
Qwen/Qwen-72B
1. YES
2. YES