GaussianProcesses_Exercise1_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "wKx8dXanZSDP"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# we set a seed to make the results reproducible\n",
    "np.random.seed(32)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "GAn8tzTFVm0Z"
   },
   "source": [
    "# Gaussian Processes - Exercise Session 1 - Solutions\n",
    "\n",
    "## Relationship between linear models and conditional Gaussian distribution\n",
    "\n",
    "In this exercise session, you are going to learn the relationship between the predictions given by a linear model and the mean estimate given by a conditional Gaussian distribution. We will derive both formula step-by-step and you should implement them in Python.\n",
    "\n",
    "\n",
    "### Dataset creation\n",
    "\n",
    "We are going to create a toy dataset $D=\\{(\\mathbf{x_1},y_1),...(\\mathbf{x_n},y_n)\\}$, where $\\mathbf{x} \\in \\mathbb{R}^3$ and $y \\in \\mathbb{R}$. The $\\mathbf{x}$ is created in the following way:\n",
    "\n",
    "\n",
    "1.   First we sample a number in $[0,1]$ for each dimension\n",
    "2.   Then we multiply each dimension with a random integer in the range $[1,20]$\n",
    "\n",
    "This way, we now have random $x$'s with very different values. Since we want to have a bias term in our linear model, we use the trick of adding a of 1 in front of our $\\mathbf{x} \\in D$.\n",
    "\n",
    "The next step is to generate the $y$ for each $\\mathbf{x}$. We have to define the real weights vector $\\mathbf{w}\\in \\mathbb{R}^4$ and then compute \n",
    "$$y = \\mathbf{w}^T \\mathbf{x} + \\epsilon$$\n",
    "where $\\epsilon$ is the zero-mean and unit variance Gaussian noise.\n",
    "\n",
    "At this point, we have our datset $D$ ready. We will split it in a training set $\\mathbf{X}$ and a test set $\\mathbf{X_*}$ (we here will adhere to a typical setting were 80% of the data goes into the training set and the remaining 20% will be reserved for testing. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='green'> There was some confusion during this exercise about the role of the bias, which translates into the 1 that we are stacking as first column in our dataset.\n",
    "\n",
    "Suppose we have a dataset $D = \\{(x_0,y_0), \\dotsc, (x_N,y_N)\\}$ where $x_i, y_i \\in \\mathbb{R}^1$. When we are using linear regression, we are assuming that our targets $y$ are given by a deterministic function $f(x,w)$ with additive Gaussian noise:\n",
    "$$y = f(x,w) + \\epsilon $$\n",
    "We also assume that $f(x,w)$ is given by a linear combination of the inputs. In our case, we have that:\n",
    "$$f(x,w) = w x $$\n",
    "\n",
    "As you can seen, the function above is a line, where $w$ is its slope. However, we are forcing this line to go through zero, because when our input $x=0$ we have that we are predicting $f(x,w)=0$, which is not always the case. Indeed, it is difficult to have centered data in a real-word application. To overcome this drawback, we learn also a bias, which represents the intercept of our line. Therefore, we have:\n",
    "\n",
    "\\begin{align}\n",
    "f(x,w) &= b + wx\\\\\n",
    "&= b * 1 + w * x \n",
    "\\end{align}\n",
    "\n",
    "This corresponds to considering our dataset being made of $\\mathbf{x}= [1, x ]$ and we want to learn the following weights $\\mathbf{w}= [b, w ]$. Note that the bias is usually denoted with $w_0$ in most of the books.\n",
    "\n",
    "However, if we center the data, the bias or intercept is not useful anymore, because it forces it to be 0. \n",
    "\n",
    "In the exercise we start by adding the column of 1s in front of the dataset because this is the typical thing to do in linear regression (as centered data is typically not required). But since the GP typically assumes centering we here use centered dataset and a bias term for the linear regression to keep thing compatible.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 102
    },
    "colab_type": "code",
    "id": "nRNyvzCXVk0O",
    "outputId": "96bee2fb-598b-4cbd-b9de-199925bb7719"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Safety check: sizes\n",
      "(80, 4)\n",
      "(80,)\n",
      "(20, 4)\n",
      "(20,)\n"
     ]
    }
   ],
   "source": [
    "## we have to create our dataset\n",
    "real_weights = np.array([3.24, 1.27, -4.52, 1.75])\n",
    "\n",
    "# number of observations\n",
    "N = 100\n",
    "# dimension of X\n",
    "d = 3\n",
    "X = np.random.rand(N,d) * np.random.randint(1,20,(N,d)) # Nxd matrix\n",
    "\n",
    "# we add the bias term adding a column of 1 in front\n",
    "# of our matrix X\n",
    "X = np.hstack((np.ones((N,1)), X)) # Nx(d+1)\n",
    "\n",
    "# compute the y's for each example\n",
    "# noisy --> add random noise\n",
    "y = real_weights @ X.T + np.random.randn(N)# N\n",
    "\n",
    "# now we can split the dataset in training and test set\n",
    "Xtrain = X[:80,:]\n",
    "ytrain = y[:80]\n",
    "\n",
    "Xtest = X[80:,:]\n",
    "ytest = y[80:]\n",
    "\n",
    "print('Safety check: sizes')\n",
    "print(Xtrain.shape)\n",
    "print(ytrain.shape)\n",
    "print(Xtest.shape)\n",
    "print(ytest.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "8zVvrACYr_Fw"
   },
   "source": [
    "### Conditional Gaussian\n",
    "\n",
    "We start by assuming that\n",
    "\n",
    "$$\\begin{bmatrix} y \\\\ \\mathbf{x} \\end{bmatrix} \\sim \\mathcal{N}\\left(\\mathbf{0}, \\begin{bmatrix} \\Sigma_{yy} & \\Sigma_{yx} \\\\ \\Sigma_{xy} & \\Sigma{xx} \\end{bmatrix} \\right) $$\n",
    "\n",
    "where $\\Sigma_{yy} \\in \\mathbb{R}$ and $\\Sigma_{yy} = \\sigma_{y}^2$, i.e. it is the variance of the $y$ because $y \\in \\mathbb{R}$. $\\Sigma_{xx}$ is instead the covariance matrix related to the input vectors $\\mathbf{x}$. \n",
    "\n",
    "<font color='blue'>Tasks:</font>\n",
    "\n",
    "<font color='blue'>**1-** Since we are assuming a zero-mean distribution, we should center or standardize our data. </font>\n",
    "\n",
    "<font color='blue'>**2-** In the distribution above, we are considering the distribution of a vector created by concatenating the training target $y_i$ to each training input $\\mathbf{x}_i$. Have a look at the python function `numpy.hstack` to implement this step.</font>\n",
    "\n",
    "<font color='blue'>**3-** At this point you can compute the sample mean and the sample covariance of this distribution. Be careful and read the documentation of `np.cov`.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "iljs7e25r-Ze",
    "outputId": "9a9f0432-ab1c-4734-8197-f50ff75459d3"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Safety check\n",
      "(80, 5)\n"
     ]
    }
   ],
   "source": [
    "# we assume that the traning data (y,x) ~ N(0,\\Sigma)\n",
    "# to have 0 mean we have to center the data\n",
    "\n",
    "## WRITE THE CODE TO CENTER THE Xs\n",
    "Xtrain_mean = np.mean(Xtrain,0)\n",
    "Xtrain = Xtrain - Xtrain_mean\n",
    "\n",
    "## WRITE THE CODE TO CENTER THE Xs\n",
    "ytrain_mean = np.mean(ytrain)\n",
    "ytrain = ytrain - ytrain_mean\n",
    "\n",
    "# create arrays [y,x1,x2,...,xd] for each array\n",
    "dataset = np.hstack((ytrain.reshape(-1,1), Xtrain))\n",
    "\n",
    "print('Safety check') # you should get a shape of (80, 5)\n",
    "print(dataset.shape)\n",
    "\n",
    "\n",
    "# COMPUTE THE SAMPLE MEAN AND SAMPLE COVARIANCE\n",
    "mu_dataset = np.mean(dataset,0)\n",
    "cov_dataset = np.cov(dataset.T)         # you should get a dxd matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5A7Cl6IcBGIK"
   },
   "source": [
    "At this point we have the distribution of $\\begin{bmatrix} y \\\\ \\mathbf{x} \\end{bmatrix}$ defined above. However, we are interested in the conditional probability $y|\\mathbf{x}$. We know from the lecture that the conditional distribution of a Gaussian distribution is still a Gaussian distribution. \n",
    "\n",
    "Therefore we are looking for:\n",
    "\n",
    "$$y | \\mathbf{x} \\sim \\mathcal{N}(\\mathbf{\\mu}_{y|\\mathbf{x}}, \\Sigma_{y|\\mathbf{x}})$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\\mathbf{\\mu}_{y|\\mathbf{x}} = \\mu_{y} + \\Sigma_{yx}\\Sigma_{xx}^{-1}(\\mathbf{x}-\\mu_{x}) $$\n",
    "\n",
    "$$\\Sigma_{y|\\mathbf{x}} = \\Sigma_{yy} - \\Sigma_{yx}\\Sigma_{xx}^{-1}\\Sigma_{xy}$$\n",
    "\n",
    "Since we have centered our data, we have that both $\\mu_{y} = 0$ and $\\mu_{x}= 0$. Therefore, the equation for computing the $\\mathbf{\\mu}_{y|\\mathbf{x}}$ becomes easier:\n",
    "\n",
    "$$\\mathbf{\\mu}_{y|\\mathbf{x}} = \\Sigma_{yx}\\Sigma_{xx}^{-1}\\mathbf{x} $$\n",
    "\n",
    "We are mostly interested in the mean of this distribution, because it is the one that give us the mean estimate for every input $\\mathbf{x}$. Indeed. if we look at the definition of the mean we can see that it is a linear function of the inputs $\\mathbf{x}$. For simplicity you can write it as $\\mathbf{\\mu}_{y|\\mathbf{x}} = \\mathbf{c} \\mathbf{x}$, where $\\mathbf{c}=\\Sigma_{yx}\\Sigma_{xx}^{-1}$. As we are going to see later, this is really similar to what we get when we are using a linear model for predictions. \n",
    "\n",
    "\n",
    "<font color='blue'>Tasks:</font>\n",
    "\n",
    "<font color='blue'> **1**- Starting from the mean and covariance you have computed at the previous step, you should compute the vector $\\mathbf{c}=\\Sigma_{yx}\\Sigma_{xx}^{-1}$, because if we know its value, we are able to compute the value of the mean $\\mathbf{\\mu}_{y|\\mathbf{x}}^*$ for every test input vectors $\\mathbf{x}_*$. To compute the inverse of a matrix you can use `np.linalg.pinv`.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 139
    },
    "colab_type": "code",
    "id": "WmdtgxqHGbV1",
    "outputId": "54a45870-477d-48d2-f81c-8271eee38687"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Safety check covariance\n",
      "()\n",
      "(4,)\n",
      "(4,)\n",
      "(4, 4)\n",
      "Safety check\n",
      "(4,)\n",
      "The C you find is:  [ 0.          1.26636126 -4.46762863  1.70423914]\n"
     ]
    }
   ],
   "source": [
    "# we are interested in p(y|x1,..xd) which is Gaussian given by\n",
    "# mu_yx = Sigmayx Simgaxx^-1 x\n",
    "\n",
    "# SEPARATE ALL THE PIECES IN THE COVARIANCE MATRIX\n",
    "Sigma_yy = cov_dataset[0,0]\n",
    "Sigma_yx = cov_dataset[0,1:]\n",
    "Sigma_xy = cov_dataset[1:,0]\n",
    "Sigma_xx = cov_dataset[1:,1:]\n",
    "\n",
    "\n",
    "print('Safety check covariance')\n",
    "print(Sigma_yy.shape)\n",
    "print(Sigma_xy.shape)\n",
    "print(Sigma_yx.shape)\n",
    "print(Sigma_xx.shape)\n",
    "\n",
    "## \n",
    "\n",
    "# COMPUTE THE INVERSE AND THEN DO THE MATRIX MULTIPLICATION \\Sigmayx\\Sigmaxx^-1\n",
    "C = Sigma_yx @ np.linalg.pinv(Sigma_xx)\n",
    "\n",
    "print('Safety check')\n",
    "print(C.shape)\n",
    "print('The C you find is: ', C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "343.2479937771417"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Sigma_yy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "MrO-2VK-IhwT"
   },
   "source": [
    "### Linear models\n",
    "\n",
    "<font color='grey'>  Notation alert! We are going to write summations in matrix notation, therefore you should be comfortable with using vectors and matrices. </font>\n",
    "\n",
    "\n",
    "Here we are going to derive again the Normal Equation that you have already seen in your Machine Learning course. In regression, we usually try to learn a model of this form:\n",
    "$$y = \\mathbf{w}^T\\mathbf{x}$$\n",
    "where $\\mathbf{w}=[w_0, w_1,\\dotsc,w_D]$ and $\\mathbf{x} = [1, x_0, x_1, \\dotsc, x_D]$.\n",
    "\n",
    "Given a dataset that contains $N$ training examples $(\\mathbf{x}_i,y_i)_{i=1,\\dotsc,N}$, we learn the weigths or parameters of the model $\\mathbf{w}$ by minimizing the Mean-Squared Error of the dataset given by\n",
    "$$ E(\\mathbf{w}) = \\frac{1}{2} \\sum_{n=1}^{N}(y_n - \\mathbf{w}^T\\mathbf{x}_n)^2$$\n",
    "\n",
    "We can rewrite $E(\\mathbf{w})$ as\n",
    "\\begin{align}\n",
    "E(\\mathbf{w}) &= \\frac{1}{2} \\sum_{n=1}^{N}(y_n - \\mathbf{w}^T\\mathbf{x}_n)^2 \\\\\n",
    "              &= \\frac{1}{2} \\left\\| \\mathbf{X}\\mathbf{w} - \\mathbf{y} \\right\\|^2 \\\\\n",
    "              &= \\frac{1}{2} (\\mathbf{X}\\mathbf{w} - \\mathbf{y})^T(\\mathbf{X}\\mathbf{w} - \\mathbf{y}) \\\\\n",
    "              &= \\frac{1}{2} \\mathbf{w}^T\\mathbf{X}^T\\mathbf{X}\\mathbf{w} - \\frac{1}{2}\\mathbf{w}^T\\mathbf{X}^T\\mathbf{y} - \\frac{1}{2} \\mathbf{y}^T\\mathbf{X}\\mathbf{w} + \\frac{1}{2} \\mathbf{y}^T\\mathbf{y} \\\\\n",
    "              &= \\frac{1}{2} \\mathbf{w}^T\\mathbf{X}^T\\mathbf{X}\\mathbf{w} - \\mathbf{w}^T\\mathbf{X}^T\\mathbf{y} + \\frac{1}{2} \\mathbf{y}^T\\mathbf{y}\n",
    "\\end{align}\n",
    "\n",
    "where in the last line we use the fact that $(\\mathbf{X}\\mathbf{w})^T \\mathbf{y} = \\mathbf{y}^T(\\mathbf{X}\\mathbf{w})$ because the inner product is commutative and the fact that $(\\mathbf{X}\\mathbf{w})^T=\\mathbf{w}^T\\mathbf{X}^T$.\n",
    "\n",
    "Since we want to find the weights $\\mathbf{w}$ that minimize $E(\\mathbf{w})$, from calculus we know that to solve this problem we can take the gradient of $E(\\mathbf{w})$, and then set it to zero and solve the resulting system. \n",
    "$$\\nabla E(\\mathbf{w}) = \\begin{bmatrix} \\frac{\\partial E}{\\partial w_0}\\\\ \\frac{\\partial E}{\\partial w_1}\\\\ \\vdots \\\\ \\frac{\\partial E}{\\partial w_D} \\end{bmatrix} = 0$$ \n",
    "\n",
    "Two compute the gradient we use two results that are important to know:\n",
    "$$ f(\\mathbf{w}) = \\mathbf{w}^T\\mathbf{A}\\mathbf{w} \\implies \\nabla f(\\mathbf{w}) = 2 \\mathbf{A}\\mathbf{w}$$\n",
    "$$ f(\\mathbf{w}) = \\mathbf{a}^T\\mathbf{w}  \\implies \\nabla f(\\mathbf{w}) = \\mathbf{a}$$\n",
    "\n",
    "Therefore we can see that the gradient of $E(\\mathbf{w})$ is given by:\n",
    "$$\\nabla E(\\mathbf{w}) = \\mathbf{X}^T \\mathbf{X} \\mathbf{w} - \\mathbf{X}^T\\mathbf{y} $$\n",
    "and if we put this to zero we get\n",
    "\\begin{align}\n",
    "\\mathbf{X}^T \\mathbf{X} \\mathbf{w} - \\mathbf{X}^T\\mathbf{y} = 0 \\\\\n",
    "\\mathbf{X}^T \\mathbf{X} \\mathbf{w} = \\mathbf{X}^T\\mathbf{y} \\\\\n",
    "\\mathbf{w}  = (\\mathbf{X}^T \\mathbf{X})^{-1}\\mathbf{X}^T\\mathbf{y} \n",
    "\\end{align}\n",
    "\n",
    "<font color='blue'> Tasks:</font>\n",
    "\n",
    "<font color='blue'> **1**- You should code the normal equation in Python. The easiest way to compute it is to write all the matrix multiplication as shown in the formula. Remember, if you want to compute the inverse of a matrix use `numpy.linalg.pinv`. </font>\n",
    "\n",
    "<font color='blue'> (There are other ways you can maybe implement the normal equation, but for this exercise we required you to use the easiest one. You maybe notice that you are trying to solve a linear system, therefore you would try to use `numpy.linalg.solve`. However, in case our system is either over or underconstrained this method will return an error. In those cases you can consider `numpy.linalg.lstsq` which will return the solution that minimizes the squared error.) </font>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "5TUI00zs1FFo",
    "outputId": "e90ade9d-b15d-48fd-bf0c-75a0763e1719"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The weigths obtained by solving the normal equation are: [ 0.          1.26636126 -4.46762863  1.70423914]\n"
     ]
    }
   ],
   "source": [
    "# WRITE THE NORMAL EQUATION\n",
    "\n",
    "w = np.linalg.pinv(Xtrain.T @ Xtrain) @ Xtrain.T @ ytrain\n",
    "\n",
    "print('The weigths obtained by solving the normal equation are:', w)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "h2ZTxI8g6aXu"
   },
   "source": [
    "### Predictions\n",
    "\n",
    "You have already seen in previous courses how to compute the prediction $y_*$ for a test example $\\mathbf{x}_*$:\n",
    "$$y_* = \\mathbf{w}^T \\mathbf{x}_*$$\n",
    "using the weights you have computed in the previous step.\n",
    "\n",
    "Maybe you have not seen before that you can use multivariate Gaussian distribution to make predictions. For predicting the target $y_*$ for a particular test input $x_*$ you should use the conditional Gaussian distribution $y_*|\\mathbf{x}_*$, that you have computed above. Indeed, you want to compute the distribution of the possible targets $y_*$ given the test input $x_*$. A good estimate of $y_*$ would be the mean of this distribution, therefore you should compute:\n",
    "$$y_* = \\mu_{y_*|\\mathbf{x}_*} = \\Sigma_{yx}\\Sigma_{xx}^{-1}\\mathbf{x}_*$$\n",
    "where the $\\Sigma_{yx}$ and $\\Sigma_{xx}$ are those you have calculated from the training set.\n",
    "\n",
    "<font color='blue'>Tasks: </font>\n",
    "\n",
    "<font color='blue'> **1**- For each test point, you should compute the prediction given by the linear model and the mena of the conditional distribution, and plot it. What do you notice? \n",
    "\n",
    "**N.B**: Remember that every modifications you did in the training set should be done also in the test set. This means that  you should center each test set (both the input $\\mathbf{x}_*$ and the true target $y$ using the means you have computed in the training set!) </font>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 129
    },
    "colab_type": "code",
    "id": "XIdKUWL6R40y",
    "outputId": "d2a1c737-e61e-418a-888b-52240c78467c"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Example no.    y (true)    y_hat (lin. reg.)    mu(y|x) (cond. mean)\n",
      "          1       12.88                12.54                   12.54\n",
      "          2       -7.98                -9.02                   -9.02\n",
      "          3      -34.22               -34.70                  -34.70\n",
      "          4       19.27                19.45                   19.45\n",
      "          5       -1.26                -0.49                   -0.49\n",
      "          6      -38.13               -39.66                  -39.66\n",
      "          7       15.70                13.94                   13.94\n",
      "          8       -6.45                -8.91                   -8.91\n",
      "          9       13.65                10.57                   10.57\n",
      "         10        1.72                 1.43                    1.43\n",
      "         11       -1.63                -1.10                   -1.10\n",
      "         12        9.19                 9.32                    9.32\n",
      "         13        2.49                 2.14                    2.14\n",
      "         14       20.84                19.08                   19.08\n",
      "         15        3.24                 2.35                    2.35\n",
      "         16       -7.77                -7.62                   -7.62\n",
      "         17       -2.11                -1.38                   -1.38\n",
      "         18        2.65                 2.21                    2.21\n",
      "         19        2.08                 2.46                    2.46\n",
      "         20       28.25                29.85                   29.85\n"
     ]
    }
   ],
   "source": [
    "# CENTER THE TEST EXAMPLES\n",
    "Xtest -= Xtrain_mean\n",
    "ytest -= ytrain_mean\n",
    "\n",
    "# COMPUTE THE PREDICTION USING THE LINEAR MODEL\n",
    "yhat = Xtest @ w.T\n",
    "\n",
    "# COMPUTE THE MEAN GIVEN BY THE CONDITIONAL GAUSSIAN\n",
    "mu_hat = Xtest @ C\n",
    "\n",
    "# PRINT THE RESULTS\n",
    "columns = ('Example no.', 'y (true)', 'y_hat (lin. reg.)', 'mu(y|x) (cond. mean)')\n",
    "print('    '.join(columns))\n",
    "for i, (y, yh, ym) in enumerate(zip(ytest, yhat, mu_hat)):\n",
    "    print('{:11d}    {:8.2f}    {:17.2f}    {:20.2f}'.format(i+1, y, yh, ym))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "V_ZQzkn_VDMo"
   },
   "source": [
    "### Why do we get this result?\n",
    "\n",
    "<font color='blue'> From the results you get, what can you notice? Try to think at the possible reasons why you observe this findings! </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='green'> ANSWER: To get an answer we should get at how the predictions are computed. For the linear model we have:\n",
    "\\begin{align}\n",
    "y_* &= \\mathbf{w}^T\\mathbf{x}_* \\\\\n",
    "y_* &= ( (\\mathbf{X}^T \\mathbf{X})^{-1}\\mathbf{X}^T\\mathbf{y} )^T\\mathbf{x}_*\n",
    "\\end{align}\n",
    "From this expression we immediately see that $\\mathbf{X}^T \\mathbf{X}$ (with appropriate scaling, $\\tfrac{1}{N-1}$) is the covariance of the training points. Likewise, $\\mathbf{X}^T\\mathbf{y}$ (again with appropriate scaling, $\\tfrac{1}{N-1}$) is the covariance between the training and test data. When we multiply the inverse training covariance and the train/test covariance the scaling factor, $\\tfrac{1}{N-1}$, disappears and we see that the two expressions are identical. While this is sufficient to argue that the results are identical (given familiarity with calculating covariance for matrices/vectors), below a more elaborate derivation is given for the sake of completeness. \n",
    "    \n",
    "To make a little easier, we denote $\\mathbf{A}= (\\mathbf{X}^T \\mathbf{X})$ since it is a $(D+1) \\times (D+1)$ matrix and $\\mathbf{b} = \\mathbf{X}^T\\mathbf{y} $ which is a vector. Therefore, we can write the prediction in a linear model as:\n",
    "\\begin{align}\n",
    "y_* &= (\\mathbf{A}^{-1}\\mathbf{b})^T\\mathbf{x}_* \\\\\n",
    "    &= \\mathbf{b}^T (\\mathbf{A}^{-1})^T\\mathbf{x}_* \\\\\n",
    "    &= \\mathbf{b}^T (\\mathbf{A}^T)^{-1}\\mathbf{x}_*\n",
    "\\end{align}\n",
    "where in the last line we use the fact that $(\\mathbf{A}^{-1})^T = (\\mathbf{A}^T)^{-1}$. </font>\n",
    "\n",
    "<font color='green'> The vector $\\mathbf{b}^T$ is a $1 \\times (D+1)$ vector while the matrix $(\\mathbf{A}^T)^{-1}$ is a $(D+1) \\times (D+1)$ matrix. Their product is therefore a $1 \\times (D+1)$ vector, as expected, given the linear model assumption of having the target being a linear combination of the input variables.</font>\n",
    "\n",
    "<font color='green'>If we look now at the way the mean of the conditional distribution $\\mu_{y_*|\\mathbf{x}_*}$ is computed for a new example $\\mathbf{x}_*$ we can notice that there is a correspondence with the linear models:\n",
    "$$y_* = \\mu_{y_*|\\mathbf{x}_*} = \\Sigma_{yx}\\Sigma_{xx}^{-1}\\mathbf{x}_*$$ </font>\n",
    "\n",
    "<font color='green'>Indeed, also in this case we have a vector $\\Sigma_{yx}$ which is the first row without the first element of the covariance matrix of the distribution of $\\begin{bmatrix} y \\\\ \\mathbf{x} \\end{bmatrix}$ and the inverse of a matrix $\\Sigma_{xx}^{-1}$ which represents the covariance of the distribution of $\\begin{bmatrix} x_1 \\\\ \\vdots \\\\ x_D \\end{bmatrix}$. The vector-matrix multiplication results in a $1 \\times (D+1)$ vector. Therefore also in this case we have that the mean, which we are using as prediction, is a linear combination of the input.  </font>\n",
    "\n",
    "<font color='green'> We can try to see if there is really a link between $\\mathbf{b}^T=\\mathbf{y}^T\\mathbf{X}$ and $\\Sigma_{yx}$ and also between $\\mathbf{A}= (\\mathbf{X}^T \\mathbf{X})$ and $\\Sigma_{xx}$. </font>\n",
    "\n",
    "<font color='green'>To do that, we should go through at how the sample covariance matrix is computed. We consider our vector $\\begin{bmatrix} y \\\\ 1 \\\\ x_1 \\\\ x_2 \\\\ x_3 \\end{bmatrix}$ and for make it more clear we define $\\begin{bmatrix} y \\\\ x_0 \\\\ x_1 \\\\ x_2 \\\\ x_3 \\end{bmatrix}$ we defined the column of 1 as $x_0$. The sample covariance among this $D+2$ variables based on $N$ observation of each, drawn from on unobserved population are given by the $D+2 \\times D+2$ matrix where the entries $q_{jk}$ are computed as:</font>\n",
    "\n",
    "<font color='green'>$$ q_{jk} = \\frac{1}{N-1} \\sum_{i=1}^{N}(x_{j}^{(i)} - \\mu_j)(x_{k}^{(i)} - \\mu_k)$$ </font>\n",
    "\n",
    "<font color='green'> and since in our case we center the data we can write this as</font>\n",
    "\n",
    "<font color='green'>$$ q_{jk} = \\frac{1}{N-1} \\sum_{i=1}^{N}x_{j}^{(i)}x_{k}^{(i)}$$ </font>\n",
    "\n",
    "<font color='green'> We can write the covariance matrix in our case as:</font>\n",
    "\n",
    "<font color='green'> $$ \\Sigma = \\begin{bmatrix} q_{y,y} & q_{y x_0} & q_{y, x_1} & q_{y x_2} & q_{y x_3} \\\\ q_{x_0,y} & q_{x_0,x_0} & q_{x_0,x_1} & q_{x_0, x_2} & q_{x_0, x_3} \\\\ q_{x_1, y} & q_{x_1, x_0} & q_{x_1, x_1} & q_{x_1, x_2} & q_{x_1, x_3} \\\\ q_{x_2, y} & q_{x_2, x_0} & q_{x_2, x_1} & q_{x_2, x_2} & q_{x_2, x_3} \\\\ q_{x_3, y} & q_{x_3, x_0} & q_{x_3, x_1} & q_{x_3, x_2} & q_{x_3, x_3} \\end{bmatrix} $$ </font>\n",
    "\n",
    "<font color='green'> Where we can separates the quantities of our interest, i.e. $\\Sigma_{yx}$ which we know that it is a vector and the sub-matrix $\\Sigma_{xx}$:</font>\n",
    "\n",
    "<font color='green'>$$ \\Sigma_{yx} = [q_{y x_0}, q_{y, x_1}, q_{y x_2},  q_{y x_3}]$$\n",
    "    \n",
    "$$ \\Sigma_{xx} = \\begin{bmatrix} q_{x_0,x_0} & q_{x_0,x_1} & q_{x_0, x_2} & q_{x_0, x_3} \\\\ q_{x_1, x_0} & q_{x_1, x_1} & q_{x_1, x_2} & q_{x_1, x_3} \\\\ q_{x_2, x_0} & q_{x_2, x_1} & q_{x_2, x_2} & q_{x_2, x_3} \\\\  q_{x_3, x_0} & q_{x_3, x_1} & q_{x_3, x_2} & q_{x_3, x_3} \\end{bmatrix} $$</font>\n",
    "\n",
    "<font color='green'>We start considering $\\Sigma_{yx}$. We can see that:</font>\n",
    "\n",
    "\n",
    "<font color='green'>\\begin{align}\n",
    "\\Sigma_{yx} &= [q_{y x_0}, q_{y, x_1}, q_{y x_2},  q_{y x_3}] \\\\\n",
    "            &= [\\frac{1}{N-1} \\sum_{i=1}^{N}y^{(i)}x_{0}^{(i)}, \\frac{1}{N-1} \\sum_{i=1}^{N}y^{(i)}x_{1}^{(i)}, \\frac{1}{N-1} \\sum_{i=1}^{N}y^{(i)}x_{2}^{(i)}, \\frac{1}{N-1} \\sum_{i=1}^{N}y^{(i)}x_{3}^{(i)}] \\\\\n",
    "            &= [\\frac{1}{N-1} (y^{(0)}x_{0}^{(0)}+y^{(1)}x_{0}^{(1)}+\\dotsc+y^{(N)}x_{0}^{(N)}), \\frac{1}{N-1} (y^{(0)}x_{1}^{(0)}+y^{(1)}x_{1}^{(1)}+\\dotsc+y^{(N)}x_{1}^{(N)}), \\frac{1}{N-1} (y^{(0)}x_{2}^{(0)}+y^{(1)}x_{2}^{(1)}+\\dotsc+y^{(N)}x_{2}^{(N)}), \\frac{1}{N-1} (y^{(0)}x_{3}^{(0)}+y^{(1)}x_{3}^{(1)}+\\dotsc+y^{(N)}x_{3}^{(N)})]\\\\\n",
    "            &= \\frac{1}{N-1} [\\mathbf{y}^T \\mathbf{x}_0, \\mathbf{y}^T \\mathbf{x}_1, \\mathbf{y}^T \\mathbf{x}_2, \\mathbf{y}^T \\mathbf{x}_3]\\\\\n",
    "            &= \\frac{1}{N-1} \\mathbf{y}^T\\mathbf{X}\n",
    "\\end{align}</font>\n",
    "\n",
    "<font color='green'>Therefore, we can see that $\\Sigma_{yy}$ is exactly the same as vector $\\mathbf{b}^T=\\mathbf{y}^T\\mathbf{X}$ used in linear models for predicting a test target. </font>\n",
    "\n",
    "\n",
    "<font color='green'>The steps to shows that also $\\Sigma_{xx}$ is equal to $\\mathbf{A}= (\\mathbf{X}^T \\mathbf{X})$ are similar, but a bit longer because it is a matrix. </font>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Exercise1.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}