GaussianProcesses_Exercise2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "CdtdoyrfJkww"
   },
   "outputs": [],
   "source": [
    "from IPython.display import clear_output\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.spatial.distance import cdist\n",
    "\n",
    "## if you want that your code is reproducible and you get the same results\n",
    "## you should set a seed. \n",
    "SEED = 42\n",
    "np.random.seed(SEED)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5QI-NFoWJkw8"
   },
   "source": [
    "## Gaussian Processes Regression\n",
    "\n",
    "In this exercise session you are going to work and become comfortable with Guassian Processes (GPs), especially in Gaussian Process regression.\n",
    "\n",
    "The notebook is structured as follows:\n",
    "\n",
    "1- You should implement the squared-exponential kernel function and understand the effects of the hyperparameters on the covariance matrix that it creates.\n",
    "\n",
    "2- Then you should implement both noise-free and noisy GP regression for a really simple toy dataset\n",
    "\n",
    "3- As last exercise, you should be ready to apply noisy-GP regression in a more realistic dataset and you should also optimize the parameters of the GP. \n",
    "\n",
    "### Kernel functions\n",
    "\n",
    "From the lecture and the lecture notes, you have seen that the Gaussain Processes assume that the target $y$'s are sampled from a Gaussian distribution. For simplicity, if we condiser our data centered, we assume a zero-mean Gaussian distribution:\n",
    "\n",
    "$$\\begin{bmatrix}y_1 \\\\ y_2\\\\ \\vdots \\\\ y_n \\end{bmatrix} \\sim \\mathcal{N}(0, \\Sigma)$$\n",
    "\n",
    "The trick of GPs is that the covariance matrix between the target $y$'s depends on the input observations $\\mathbf{x}$. This means that each entry of the covariance matrix is given by:\n",
    "\n",
    "$$\\Sigma_{ij} = \\text{Cov}(y_i,y_j)=\\text{Cov}(f(x_i),f(x_j))=k(x_i,x_j)$$\n",
    "\n",
    "where $k(x_i,x_j)$ is the *kernel* or *covariance function*.\n",
    "\n",
    "We have only introduced one kernel, although there are different choices for a kernel and if in the future you are going to work with GPs in the future. We will focus only on the *Squared-Exponential* (SE) kernel, defined by:\n",
    "\n",
    "$$k_{SE}(x_i,x_j)= \\sigma^2 \\exp \\left(-\\frac{\\left\\|x_i - x_k \\right\\|^2}{2l^2}\\right)$$\n",
    "\n",
    "that has two parameters: the lengthscale $l$ and output variance $\\sigma^2$.\n",
    "\n",
    "<font color='blue'> Tasks:\n",
    "\n",
    "<font color='blue'>**1-** You should implement a function `kernel(set1, set2, lengthscale, output_variance)`. Once you call this function it should return a matrix that contains the pairwise evaluation of the kernel. \n",
    "    **HINT**: you should start by computing all the pair-wise distances, a way to get it is to write the expand form of $\\left\\|x_i - x_k \\right\\|^2$, otherwise you should look at the function `scipy.spatial.distance.cdist`. Then you use them to compute the kernel.\n",
    "    \n",
    "<font color='blue'>**2-** Try to apply the function using `lengthscale = 1` and `output_variance = 1` to a toy set and look at the plot. You should expect the value $1$ along the diagonal and that closer points are higher correlated than further points. (If you use `plt.colorbar()`, you will get the mapping from color to numbers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 129
    },
    "colab_type": "code",
    "id": "UjfIPfj0JkxC",
    "outputId": "5677b662-0842-4f4a-8030-4c15b3cf1752"
   },
   "outputs": [],
   "source": [
    "from IPython.display import clear_output\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.spatial.distance import cdist\n",
    "\n",
    "## if you want that your code is reproducible and you get the same results\n",
    "## you should set a seed. \n",
    "SEED = 42\n",
    "np.random.seed(SEED)\n",
    "\n",
    "## IMPLEMENT THE KERNEL FUNCTION\n",
    "def squared_exponential_kernel(x, y, lengthscale, variance):\n",
    "    # COMPUTE ALL THE PAIRWISE DISTANCES, size: NxM\n",
    "    sqdist = cdist(x,y, metric = 'sqeuclidean')\n",
    "    # COMPUTE THE KENEL, this should also be NxM\n",
    "    k = variance * np.exp(-sqdist/(2*lengthscale**2))\n",
    "    return k\n",
    "\n",
    "\n",
    "\n",
    "## TOY DATASET\n",
    "temp_x = np.linspace(0,1,5)\n",
    "\n",
    "# we apply the kernel and get the covariance\n",
    "## DEFINE LENGTHSCALE AND OUTPUT_VARIANCE\n",
    "temp_lengthscale = 1\n",
    "temp_covariance = 0.1\n",
    "\n",
    "temp_cov = squared_exponential_kernel(temp_x.reshape(-1,1), temp_x.reshape(-1,1), temp_lengthscale, temp_covariance)\n",
    "\n",
    "# plot\n",
    "plt.imshow(temp_cov)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "EGf_AvdBJkxM"
   },
   "source": [
    "### Sampling from the prior distribution\n",
    "\n",
    "At this point, as you have seen from the lecture notes, you should know how sample from the prior distribution of our Gaussian Process. We are going to resume the process here. You can start by defining $N_*$ point on the $x$-axis, which can technically be considered as test point since you do not know their target. Since we have no trainign points $(x_i,y_i)_{i=1,\\dotsc,N}$, we compute the covarince matrix $\\Sigma_{**}$ of the distribution using only the $N_*$ points. This will reuslt in a $N_* \\times N_*$ matrix. If you define your prior distribution as $\\mathcal{N}(\\mathbf{0}, \\Sigma_{**})$ and you sample 5 times from it you will get 5 possible values for each point $x_*$ in the test set:\n",
    "\n",
    "$$\\begin{bmatrix}y_{*,1} \\\\ y_{*,2}\\\\ \\vdots \\\\ y_{*,N_*} \\end{bmatrix} \\sim \\mathcal{N}(0, \\Sigma_{**})$$\n",
    "\n",
    "\n",
    "<font color='blue'> Tasks:\n",
    "    \n",
    "<font color='blue'><font color='blue'>**1-** Create the covariance matrix $\\Sigma_{**}$ using the kernel function you implemented before. (Before doing this step be sure that the function you implemented actually works as expected!) As starting point, you can use `lengthscale = 1` and `output_variance = 1` as kernel parameters.\n",
    "\n",
    "<font color='blue'>**2-** Create the distribution $\\mathcal{N}(\\mathbf{0}, \\Sigma_{**})$ and sample $5$ samples from it (**HINT**: you can sample using the following function `np.random.multivariate_normal`). Otherwise if do not wish to rely on this function, sampling from a centered, uncorrelated, homoscedastic Normal with scalar variance of 1 and left multiplying with the matrix squareroot/Cholesky factorization of $\\Sigma$ and adding the mean will do the same (convince yourself that this is true). Create and show the following plots (most code for the plots is given):\n",
    "1. The covariance matrix as you did before\n",
    "2. The prior distribution: for each point $x_*$ plot its mean and the two-standard deviation interval. (**HINT**: to get the standard deviation of each point you should take the square-root of the elements in the diagonal, look at the function `np.diag`. To plot a shaded area for the standard deviation you can use the following function `plt.gca().fill_between(Xtest.flat, mu.reshape(-1) - 2 * var, mu.reshape(-1) + 2 * var,  color='lightblue', alpha=0.5)` where obsviously you should substitute the variables \\texttt{Xtest, mu, var} with the name that you have chosen.\n",
    "3. Plot the different samples. If you use different colors the resuts looks better and clearer.\n",
    "\n",
    "<font color='blue'>**3-** Try to understand the effect of the lengthscale and the output variance. Repeat the steps above selecting different values for the two hyperparameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.105427357601002e-17"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lengthscale = 1\n",
    "output_variance = 1\n",
    "\n",
    "\n",
    "npoints = 50\n",
    "Xtest = np.linspace(-5,5,npoints)\n",
    "\n",
    "# mean and covariance of the distribution\n",
    "# we asusme zero mean, so centered data\n",
    "mu = np.mean(Xtest)\n",
    "mu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 129
    },
    "colab_type": "code",
    "id": "ymaiUuAJJkxP",
    "outputId": "0a3c68ea-3f71-43b9-fcc2-3e4c2e89e9bd",
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5,)\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "Invalid shape () for image data",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_8980/1004972551.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     35\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     36\u001b[0m \u001b[1;31m# we can plot the covariance matrix\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 37\u001b[1;33m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mimshow\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcovariance\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     38\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcolorbar\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     39\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mxlabel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Xtrain elements'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfontsize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m13\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\anaconda3\\lib\\site-packages\\matplotlib\\pyplot.py\u001b[0m in \u001b[0;36mimshow\u001b[1;34m(X, cmap, norm, aspect, interpolation, alpha, vmin, vmax, origin, extent, filternorm, filterrad, resample, url, data, **kwargs)\u001b[0m\n\u001b[0;32m   2901\u001b[0m         \u001b[0mfilternorm\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfilterrad\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m4.0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mresample\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0murl\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2902\u001b[0m         data=None, **kwargs):\n\u001b[1;32m-> 2903\u001b[1;33m     __ret = gca().imshow(\n\u001b[0m\u001b[0;32m   2904\u001b[0m         \u001b[0mX\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcmap\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcmap\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnorm\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnorm\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maspect\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0maspect\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2905\u001b[0m         \u001b[0minterpolation\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0minterpolation\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0malpha\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0malpha\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvmin\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mvmin\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\anaconda3\\lib\\site-packages\\matplotlib\\__init__.py\u001b[0m in \u001b[0;36minner\u001b[1;34m(ax, data, *args, **kwargs)\u001b[0m\n\u001b[0;32m   1359\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0minner\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0max\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1360\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mdata\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1361\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0max\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0mmap\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msanitize_sequence\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1362\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1363\u001b[0m         \u001b[0mbound\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnew_sig\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbind\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0max\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\anaconda3\\lib\\site-packages\\matplotlib\\axes\\_axes.py\u001b[0m in \u001b[0;36mimshow\u001b[1;34m(self, X, cmap, norm, aspect, interpolation, alpha, vmin, vmax, origin, extent, filternorm, filterrad, resample, url, **kwargs)\u001b[0m\n\u001b[0;32m   5607\u001b[0m                               resample=resample, **kwargs)\n\u001b[0;32m   5608\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 5609\u001b[1;33m         \u001b[0mim\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_data\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   5610\u001b[0m         \u001b[0mim\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_alpha\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0malpha\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   5611\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mim\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_clip_path\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\anaconda3\\lib\\site-packages\\matplotlib\\image.py\u001b[0m in \u001b[0;36mset_data\u001b[1;34m(self, A)\u001b[0m\n\u001b[0;32m    707\u001b[0m         if not (self._A.ndim == 2\n\u001b[0;32m    708\u001b[0m                 or self._A.ndim == 3 and self._A.shape[-1] in [3, 4]):\n\u001b[1;32m--> 709\u001b[1;33m             raise TypeError(\"Invalid shape {} for image data\"\n\u001b[0m\u001b[0;32m    710\u001b[0m                             .format(self._A.shape))\n\u001b[0;32m    711\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mTypeError\u001b[0m: Invalid shape () for image data"
     ]
    }
   ],
   "source": [
    "import random\n",
    "\n",
    "\n",
    "# TEST SET (OR RANGE OF INTEREST)\n",
    "lengthscale = 1\n",
    "output_variance = 1\n",
    "\n",
    "\n",
    "npoints = 50\n",
    "Xtest = np.linspace(-5,5,npoints)\n",
    "\n",
    "# mean and covariance of the distribution\n",
    "# we asusme zero mean, so centered data\n",
    "mu = np.mean(Xtest)\n",
    "\n",
    "# since we are using linspace between -5 and 5, our data is already centered, i.e. the mean\n",
    "# of Xtest is zero. You can check with np.mean(Xtest)\n",
    "## COMPUTE THE COVARIANCE MATRIX\n",
    "covariance = np.cov(Xtest)\n",
    "\n",
    "## GET THE STANDARD DEVIATION FROM THE COVARIANCE\n",
    "std = np.std(Xtest)\n",
    "\n",
    "# samples from the prior distribution\n",
    "n_samples = 5\n",
    "\n",
    "## SAMPLE N_SAMPLES POINTS FROM THE PRIOR DISTRIBUTION\n",
    "ys = np.random.choice(Xtest,n_samples)\n",
    "\n",
    "print(ys.shape)\n",
    "\n",
    "\n",
    "#######\n",
    "# the 3 plots we need\n",
    "\n",
    "# we can plot the covariance matrix\n",
    "plt.imshow(covariance)\n",
    "plt.colorbar()\n",
    "plt.xlabel('Xtrain elements', fontsize=13)\n",
    "plt.ylabel('Xtrain elements', fontsize=13)\n",
    "plt.title('Covariance matrix')\n",
    "plt.show()\n",
    "\n",
    "# the prior distribution\n",
    "plt.gca().fill_between(Xtest.flat, mu - 2 * std, mu + 2 * std,  color='lightblue', alpha=0.5, label = r\"$2\\sigma$\")\n",
    "plt.plot(Xtest, mu, 'blue', label = r\"$\\mu$\")\n",
    "plt.xlabel('$x$', fontsize=13)\n",
    "plt.ylabel('$y = f(x)$', fontsize=13)\n",
    "plt.title(('Prior of our Gaussian processes'))\n",
    "plt.legend()\n",
    "plt.axis([-5, 5, -5, 5])\n",
    "plt.show()\n",
    "\n",
    "# samples from the prior\n",
    "plt.figure(figsize=(6, 4))\n",
    "for i in range(n_samples):\n",
    "    plt.plot(Xtest, ys[i], linestyle='-', marker='o', markersize=3)\n",
    "plt.xlabel('$x$', fontsize=13)\n",
    "plt.ylabel('$y = f(x)$', fontsize=13)\n",
    "plt.title((\n",
    "    '5 different function realizations at 50 points\\n'\n",
    "    'sampled from a Gaussian process with squared exponentiated kernel'))\n",
    "plt.axis([-5, 5, -5, 5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "SX1XzYBJJkxW"
   },
   "source": [
    "## Samples from the posterior\n",
    "\n",
    "The goal of Gaussian Process is not to get samples from the prior, but instead, we want to use the knowledge given by the training data about the underlying function to be able to make predictions for new, unseen test points. In the next exercises we are going to understand how to use training points to get information about the posisble regression fucntions that we are trying to model, and as next natural step, how to make predictions for test points.\n",
    "\n",
    "### Noise-free observations\n",
    "We start by considering the training observations we get to be noise-free. This is an unrealistic assumption to make in a real-world scenario, but at the same time, it makes the computation of the posterior easier to understand. In a more analytical way we are considering that the targets we recorded are exactly the underlying function that we want to learn, i.e. $y = f(\\mathbf{x},\\mathbf{w})$, where the usual zero-mean normal error $\\epsilon$ does not appear. \n",
    "\n",
    "Since the prior assumption of having a zero-mean distribution, we should standardize our data. This way we get also that the variance is 1. As you have read in the lecture notes, this would help in preventing numerical issues in computing the posterior distribution. The starting point is to deifne the joint distribution of the training observations $\\mathbf{y}$ and the test outputs $\\mathbf{y}_*$ we want to observe as:\n",
    "\n",
    "$$\\begin{bmatrix} \\mathbf{y} \\\\ \\mathbf{y}_* \\end{bmatrix} \\sim \\left( \\mathbf{0}, \\begin{bmatrix}\n",
    "    K(\\mathbf{X},\\mathbf{X}) & K(\\mathbf{X}, \\mathbf{X_*})\\\\\n",
    "    K(\\mathbf{X_*},\\mathbf{X}) & K(\\mathbf{X_*}, \\mathbf{X_*})\n",
    "    \\end{bmatrix}{}\\right) $$\n",
    "    \n",
    "where $\\mathbf{X}$ is the matrix containing the training points $\\mathbf{x}_i$ and $\\mathbf{X_*}$ is the test set. If there are $N$ training points and $N_*$ test points we have that: \n",
    "\n",
    "- $K(X,X)$ is the $N \\times N$ matrix of the kernel evaluated at all training points pairs, \n",
    "\n",
    "- $K(X,X_*)$ is the $N \\times N_*$ matrix with the kernels evaluated at all pairs created by combining the training points and test points. \n",
    "\n",
    "- $K(X_*,X_*)$ contains the $N_* \\times N_*$ elements obtained by computing the covariance functions to all pairs of test points.\n",
    "\n",
    "We are interested in the conditional probability $\\mathbf{y}_*|\\mathbf{y}$, which is given by $p(\\mathbf{y_*}|\\mathbf{y}, \\mathbf{X}, \\mathbf{X_*})$. We know how to compute the conditional probability of a multivariate Gaussian. We have seen it in the exercise session 1. If we use the same rules here, we get:\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "    p(\\mathbf{y_*}|\\mathbf{y}, \\mathbf{X}, \\mathbf{X_*}) &= \\mathcal{N}(\\mu_{\\mathbf{y_*}|\\mathbf{y}}, \\mathbf{\\Sigma}_{\\mathbf{y_*}|\\mathbf{y}})\\\\\n",
    "    \\mu_{\\mathbf{y_*}|\\mathbf{y}} &= K(\\mathbf{X_*}, \\mathbf{X}) K(\\mathbf{X},\\mathbf{X})^{-1}\\mathbf{y}\\\\\n",
    "    \\mathbf{\\Sigma}_{\\mathbf{y_*}|\\mathbf{y}} &= K(\\mathbf{X_*}, \\mathbf{X_*}) - K(\\mathbf{X_*},\\mathbf{X})K(\\mathbf{X},\\mathbf{X})^{-1}K(\\mathbf{X}, \\mathbf{X_*})\n",
    "\\end{split}\n",
    "\\end{align}\n",
    "\n",
    "To make the posterior more compact, we can define\n",
    "- $\\mathbf{K} = K(\\mathbf{X},\\mathbf{X})$,\n",
    "- $\\mathbf{K_*} = K(\\mathbf{X},\\mathbf{X_*})$ \n",
    "- $\\mathbf{K_{**}}=K(\\mathbf{X_*}, \\mathbf{X_*})$\n",
    "\n",
    "and the posterior becomes\n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "    p(\\mathbf{y_*}|\\mathbf{y}, \\mathbf{X}, \\mathbf{X_*}) &= \\mathcal{N}(\\mu_{\\mathbf{y_*}|\\mathbf{y}}, \\mathbf{\\Sigma}_{\\mathbf{y_*}|\\mathbf{y}})\\\\\n",
    "    \\mu_{\\mathbf{y_*}|\\mathbf{y}} &= \\mathbf{K_*}^T \\mathbf{K}^{-1}\\mathbf{y}\\\\\n",
    "    \\mathbf{\\Sigma}_{\\mathbf{y_*}|\\mathbf{y}} &= \\mathbf{K_{**}} - \\mathbf{K_*}^T\\mathbf{K}^{-1}\\mathbf{K_*}\n",
    "\\end{split}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 231
    },
    "colab_type": "code",
    "id": "4q17qyqSJkxY",
    "outputId": "8190412b-ce47-4802-ebc3-bd01b70f9abe"
   },
   "outputs": [],
   "source": [
    "## TOY DATASET CREATION - DO NOT CHANGE IT\n",
    "## NOISE-FREE GP\n",
    "np.random.seed(SEED)\n",
    "\n",
    "## we start creating a simple dataset\n",
    "# true unknown function we are trying to approximate\n",
    "f = lambda x: np.sin(0.9*x).flatten()\n",
    "\n",
    "N = 10 # number of training points\n",
    "\n",
    "# create the dataset\n",
    "Xtrain = np.random.uniform(-5, 5, size=(N,1))\n",
    "ytrain = f(Xtrain) \n",
    "\n",
    "# we shoul create a test set, which are point in which we want to evaluate the function \n",
    "# with a mean and a variance\n",
    "Xtest = np.linspace(-7,7, npoints).reshape(-1,1)\n",
    "\n",
    "## we can plot the training points\n",
    "plt.plot(Xtrain, ytrain, 'x', label='Training points')\n",
    "plt.plot(Xtest, f(Xtest), '-', label = 'True underline function')\n",
    "plt.xlabel('$x$', fontsize=13)\n",
    "plt.ylabel('$y = f(x)$', fontsize=13)\n",
    "plt.title('Training points')\n",
    "plt.axis([-7, 7, -2, 2])\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zJe8wLz-Jkxf"
   },
   "source": [
    "#### Noise-Free GP regression algorithm\n",
    "To robustly and efficiently caluculate the posterior that we are interested in, you can follow step-by-step this recipe:\n",
    "\n",
    "1. Start computing $\\mathbf{K} = K(\\mathbf{X},\\mathbf{X})$, $\\mathbf{K_*} = K(\\mathbf{X},\\mathbf{X_*})$ and $\\mathbf{K_{**}}=K(\\mathbf{X_*},\\mathbf{X_*})$\n",
    "2. $\\mathbf{L} = \\text{Cholesky}(\\mathbf{K})$\n",
    "3. $\\mathbf{\\alpha}= \\mathbf{L}^T \\backslash (\\mathbf{L} \\backslash \\mathbf{y})$\n",
    "4. $\\mu_{\\mathbf{y_*}|\\mathbf{y}} = \\mathbf{K_*}^T \\mathbf{\\alpha}$\n",
    "5. $\\mathbf{v} = \\mathbf{L}\\backslash \\mathbf{K_*}$\n",
    "6. $\\mathbf{\\Sigma}_{\\mathbf{y_*}|\\mathbf{y}} = K(\\mathbf{X_*},\\mathbf{X_*}) - \\mathbf{v}^T\\mathbf{v}$\n",
    "\n",
    "where the backslash symbol $\\backslash$ means that we have to solve a linear system.\n",
    "\n",
    "\n",
    "<font color='blue'> Tasks:\n",
    "1. Following the algorithm described above, you should compute the posterior distribution. Plot the mean function and the 2-standard deviation interval. Remember that you should choose the hyperparameters for the kernel. (**HINT**: for the Cholesky decomposition look at the NumPy function `np.linalg.cholesky`. To solve a linear system look at `np.linalg.solve`).\n",
    "2. Sample 5 functions from the posterior and plot them.\n",
    "</font>\n",
    "Note that using the NumPy solve operation is not really more efficient with structured matrices (tirangular). For further optimization of efficiency you can look into scipy.linalg.cho_solve and scipy.linalg.cho_factor, however for the dimensions we are working with here this is inconsequential."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 129
    },
    "colab_type": "code",
    "id": "mPzcJByuJkxh",
    "outputId": "8c260ed2-7385-474a-b0d6-1bd3b1204147"
   },
   "outputs": [],
   "source": [
    "# we follow step-by-step the algorithm \n",
    "lengthscale = \n",
    "kernel_variance = \n",
    "\n",
    "\n",
    "# we start computing the kernels (I am using s instead of * in the kernel variable names)\n",
    "K =                        # n_train x n_train\n",
    "Ks =                       # n_train x n_test\n",
    "Kss =                      # n_test x n_test\n",
    "\n",
    "# compute the cholesky decomposition\n",
    "L = \n",
    "\n",
    "# compute alpha\n",
    "alpha = \n",
    "\n",
    "# compute the mean function\n",
    "mu = \n",
    "\n",
    "# compute v\n",
    "v = \n",
    "\n",
    "# compute the covariance\n",
    "covariance = \n",
    "\n",
    "# we get the standard deviation from the covariance matrix\n",
    "std = \n",
    "\n",
    "\n",
    "## now we can plot the mean function with the standard deviation intervals\n",
    "plt.plot(Xtrain, ytrain, 'ro', label='Training points')\n",
    "plt.gca().fill_between(Xtest.flat, mu - 2 * std, mu + 2 * std,  color='lightblue', alpha=0.5, label=r\"$\\mu$\")\n",
    "plt.plot(Xtest, mu, 'blue', label=r\"$2\\sigma$\")\n",
    "plt.axis([-7, 7, -3, 3])\n",
    "plt.legend()\n",
    "plt.show() \n",
    "\n",
    "\n",
    "## SAMPLE FROM THE POSTERIOR\n",
    "samples_from_the_posterior = \n",
    "\n",
    "# PLOT THEM\n",
    "plt.plot(Xtrain, ytrain, 'r+', ms=15, label='Training points')\n",
    "#plt.gca().fill_between(Xtest.flat, mu - 2 * var, mu + 2 * var, color='lightblue', alpha=0.5)\n",
    "#plt.plot(Xtest, mu, 'blue')\n",
    "for sample_id in range(n_samples):\n",
    "    plt.plot(Xtest, samples_from_the_posterior[sample_id], linestyle='-', marker='o', markersize=3)\n",
    "plt.axis([-7, 7, -3, 3])\n",
    "plt.legend()\n",
    "plt.show()  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "E4D4NR9DJkxn"
   },
   "source": [
    "### Noisy observations\n",
    "\n",
    "The next exercise is very similar to the one you just did, but this will be the usual way to use Gaussian Porcess in a real-world setting. In this case we are going to consider that the targets $y$ we observed and registered as training set are noisy. Therefore we have: $y = f(\\mathbf{x},\\mathbf{w})+\\epsilon$, where $\\epsilon \\sim \\mathcal{N}(\\mathbf{0},\\sigma_n^2)$ is the additive independent identically distributed Gaussian noise.\n",
    "\n",
    "Therefore, when we are going to define the joint distribution over the training observations $\\mathbf{y}$ and the test observations $y_*$, we should keep the variance of the additive noise into account. The prior on the noisy training observations becomes:\n",
    "\n",
    "$$\\Sigma_{ij} = \\text{Cov}(y_i,y_j) = k(\\mathbf{x}_i, \\mathbf{x}_j) + \\sigma_n^2 \\delta_{ij}$$\n",
    "\n",
    "where with $\\Sigma_{ij}$ we are denoting the entry $i,j$ of the $K(\\mathbf{X},\\mathbf{X})$ covariance matrix. In addition to that, the $\\delta_{ij}$ is a Kronecker delta which is equal to $1$ if and only if $i=j$, and zero otherwise.\n",
    "\n",
    "Therefore we can write in matrix notation:\n",
    "\n",
    "$$ \\text{Cov}(\\mathbf{y}) = K(\\mathbf{X},\\mathbf{X}) + \\sigma_n^2 \\mathbf{I}$$\n",
    "\n",
    "and using this result, we have that the joint distribution can be defined ad:\n",
    "\n",
    "\\begin{align}\n",
    "    \\begin{bmatrix}\n",
    "    \\mathbf{y}\\\\\n",
    "    \\mathbf{y_*}\n",
    "    \\end{bmatrix} \\sim \\left( \\mathbf{0}, \\begin{bmatrix}\n",
    "    K(\\mathbf{X},\\mathbf{X})+ \\sigma_n^2 \\mathbf{I} & K(\\mathbf{X}, \\mathbf{X_*})\\\\\n",
    "    K(\\mathbf{X_*},\\mathbf{X}) & K(\\mathbf{X_*}, \\mathbf{X_*})\n",
    "    \\end{bmatrix}{}\\right)\n",
    "\\end{align}\n",
    "\n",
    "Since we usually do not know the noise variance of our data, we have that $\\sigma_n^2$ is another hyperparameter. As before, we are interested in the conditional distribution $\\mathbf{y}_*|\\mathbf{y}$, which is given by $p(\\mathbf{y_*}|\\mathbf{y}, \\mathbf{X}, \\mathbf{X_*})$. If we compute it, we get:\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "    p(\\mathbf{y_*}|\\mathbf{y}, \\mathbf{X}, \\mathbf{X_*}) &= \\mathcal{N}(\\mu_{\\mathbf{y_*}|\\mathbf{y}}, \\mathbf{\\Sigma}_{\\mathbf{y_*}|\\mathbf{y}})\\\\\n",
    "    \\mu_{\\mathbf{y_*}|\\mathbf{y}} &= K(\\mathbf{X}, \\mathbf{X_*}) [K(\\mathbf{X},\\mathbf{X})+ \\sigma_n^2 \\mathbf{I}]^{-1}\\mathbf{y}\\\\\n",
    "    \\mathbf{\\Sigma}_{\\mathbf{y_*}|\\mathbf{y}} &= K(\\mathbf{X_*}, \\mathbf{X_*}) - K(\\mathbf{X_*},\\mathbf{X})[K(\\mathbf{X},\\mathbf{X})+ \\sigma_n^2 \\mathbf{I}]^{-1}K(\\mathbf{X}, \\mathbf{X_*})\n",
    "\\end{split}\n",
    "\\end{align}\n",
    "\n",
    "and using the same compact definition we have defined above, we get:\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "    p(\\mathbf{y_*}|\\mathbf{y}, \\mathbf{X}, \\mathbf{X_*}) &= \\mathcal{N}(\\mu_{\\mathbf{y_*}|\\mathbf{y}}, \\mathbf{\\Sigma}_{\\mathbf{y_*}|\\mathbf{y}})\\\\\n",
    "    \\mu_{\\mathbf{y_*}|\\mathbf{y}} &= \\mathbf{K_*}^T [\\mathbf{K} + \\sigma_n^2 \\mathbf{I}]^{-1}\\mathbf{y}\\\\\n",
    "    \\mathbf{\\Sigma}_{\\mathbf{y_*}|\\mathbf{y}} &= \\mathbf{K_{**}} - \\mathbf{K_*}^T[\\mathbf{K}+ \\sigma_n^2 \\mathbf{I}]^{-1}\\mathbf{K_*}\n",
    "\\end{split}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 231
    },
    "colab_type": "code",
    "id": "f03revLWJkxp",
    "outputId": "dc1b26ff-c129-4bbb-e483-5a385495c9d8"
   },
   "outputs": [],
   "source": [
    "## TOY DATASET CREATION - DO NOT CHANGE IT\n",
    "## NOISY-GP regression\n",
    "np.random.seed(SEED)\n",
    "\n",
    "## we start creating a simple dataset\n",
    "# true unknown function we are trying to approximate\n",
    "f = lambda x: np.sin(0.9*x).flatten()\n",
    "\n",
    "N = 10 # number of training points\n",
    "noise_variance =  0.1 # we assume noisy data\n",
    "\n",
    "# create the dataset\n",
    "Xtrain = np.random.uniform(-5, 5, size=(N,1))\n",
    "ytrain = f(Xtrain) + noise_variance*np.random.randn(N)\n",
    "\n",
    "# we shoul create a test set, which are point in which we want to evaluate the function \n",
    "# with a mean and a variance\n",
    "Xtest = np.linspace(-7,7, npoints).reshape(-1,1)\n",
    "\n",
    "## we can plot the training points\n",
    "plt.plot(Xtrain, ytrain, 'x', label='Training points')\n",
    "plt.plot(Xtest, f(Xtest), '-', label = 'True underline function')\n",
    "plt.xlabel('$x$', fontsize=13)\n",
    "plt.ylabel('$y = f(x)$', fontsize=13)\n",
    "plt.title('Training points')\n",
    "plt.axis([-7, 7, -2, 2])\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "tW2itYUYJkxx"
   },
   "source": [
    "#### Noisy GP regression algorithm\n",
    "To compute efficiently the posterior that we are interested in, you can follow step-by-step this recipe:\n",
    "\n",
    "1. Start computing $\\mathbf{K} = K(\\mathbf{X},\\mathbf{X})$, $\\mathbf{K_*} = K(\\mathbf{X},\\mathbf{X_*})$ and $\\mathbf{K_{**}}=K(\\mathbf{X_*},\\mathbf{X_*})$\n",
    "2. $\\mathbf{L} = \\text{Cholesky}([\\mathbf{K} + \\sigma_n^2 \\mathbf{I}])$\n",
    "3. $\\mathbf{\\alpha}= \\mathbf{L}^T \\backslash (\\mathbf{L} \\backslash \\mathbf{y})$\n",
    "4. $\\mu_{\\mathbf{y_*}|\\mathbf{y}} = \\mathbf{K_*}^T \\mathbf{\\alpha}$\n",
    "5. $\\mathbf{v} = \\mathbf{L}\\backslash \\mathbf{K_*}$\n",
    "6. $\\mathbf{\\Sigma}_{\\mathbf{y_*}|\\mathbf{y}} = K(\\mathbf{X_*},\\mathbf{X_*}) - \\mathbf{v}^T\\mathbf{v}$\n",
    "\n",
    "where the backslash symbol $\\backslash$ means that we have to solve a linear system.\n",
    "\n",
    "<font color='blue'> Tasks:\n",
    "1. Following the algorithm described above, you should compute the posterior distribution. Plot the mean function and the 2-standard deviation interval. Remember that you should choose the hyperparameters for the kernel and this time you have also to choose the value for the noise variance. (**HINT**: for the Cholesky decomposition look at the NumPy function `np.linalg.cholesky`. To solve a linear system look at `np.linalg.solve`. To implement $ \\sigma_n^2 \\mathbf{I}$ you should check `np.identity`).\n",
    "2. Sample 5 functions from the posterior and plot them.\n",
    "3. repeat the same procedure for different values of the hyperparameters. What happens?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 129
    },
    "colab_type": "code",
    "id": "TCKmniyPJkxz",
    "outputId": "c5762b7d-193c-4c27-ab47-0eab6e6d59f5"
   },
   "outputs": [],
   "source": [
    "# we follow step-by-step the algorithm \n",
    "lengthscale = \n",
    "kernel_variance = \n",
    "# this time we have to consider a new hyperparameters: the noise variance\n",
    "noise_variance =  \n",
    "\n",
    "# we start computing the kernels (I am using s instead of * in the kernel variable names)\n",
    "K =\n",
    "Ks = \n",
    "Kss = \n",
    "\n",
    "# compute the cholesky decomposition\n",
    "L = \n",
    "\n",
    "# compute alpha\n",
    "alpha = \n",
    "\n",
    "# compute the mean function\n",
    "mu =\n",
    "\n",
    "# compute v\n",
    "v = \n",
    "\n",
    "# compute the covariance\n",
    "covariance =\n",
    "\n",
    "# we get the standard deviation from the covariance matrix\n",
    "std = \n",
    "\n",
    "\n",
    "## now we can plot the mean function with the standard deviation intervals\n",
    "# we get the standard deviation from the covariance matrix\n",
    "plt.plot(Xtrain, ytrain, 'ro', label='Training points')\n",
    "plt.gca().fill_between(Xtest.flat, mu - 2 * std, mu + 2 * std,  color='lightblue', alpha=0.5, label=r\"$\\mu$\")\n",
    "plt.plot(Xtest, mu, 'blue', label=r\"$2\\sigma$\")\n",
    "plt.axis([-7, 7, -3, 3])\n",
    "plt.legend()\n",
    "plt.show() \n",
    "\n",
    "\n",
    "## SAMPLE FROM THE POSTERIOR\n",
    "samples_from_the_posterior = \n",
    "\n",
    "# PLOT THEM\n",
    "plt.plot(Xtrain, ytrain, 'r+', ms=15, label='Training points')\n",
    "#plt.gca().fill_between(Xtest.flat, mu - 2 * var, mu + 2 * var, color='lightblue', alpha=0.5)\n",
    "#plt.plot(Xtest, mu, 'blue')\n",
    "for sample_id in range(n_samples):\n",
    "    plt.plot(Xtest, samples_from_the_posterior[sample_id], linestyle='-', marker='o', markersize=3)\n",
    "plt.axis([-7, 7, -3, 3])\n",
    "plt.legend()\n",
    "plt.show()  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "yZF-e5NbJkx6"
   },
   "source": [
    "## GP-regression and hyperparameters optimization on a real dataset\n",
    "\n",
    "At this point you have the basics of Gaussian Process regression. There is one last thing that you should learn. We have seen that in noisy-GP regression, we have to choose three different hyperparameters: the kernel's lengthscale and output variance, and the noise variance of the data we are trying to fit. A natural way, but not so efficient, is to do a grid search on the values of the three hyperparameters and choose the ones that maximizes the log-likelihood of our training data. This means:\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "\\log p(\\mathbf{y}|\\mathbf{X}) = \\log \\mathcal{N}(\\mathbf{y}|\\mathbf{0}, \\mathbf{K}+\\sigma_n^2 \\mathbf{I}) = - \\frac{1}{2}\\mathbf{y}[\\mathbf{K}+\\sigma_n^2 \\mathbf{I}]^{-1}\\mathbf{y} -\\frac{1}{2}|\\mathbf{K}+\\sigma_n^2 \\mathbf{I}|-\\frac{N}{2}\\log(2\\pi)\n",
    "\\end{split}\n",
    "\\end{align}\n",
    "\n",
    "If we use the Cholesky decomposition, we have to maximize:\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "\\log p(\\mathbf{y}|\\mathbf{X}) = \\log \\mathcal{N}(\\mathbf{y}|\\mathbf{0}, \\mathbf{K}+\\sigma_n^2 \\mathbf{I}) = - \\frac{1}{2}\\mathbf{y}^T\\alpha - \\sum_i \\log L_{ii} -\\frac{N}{2}\\log(2\\pi)\n",
    "\\end{split}\n",
    "\\end{align}\n",
    "\n",
    "In the next exercise we are going to use a dataset that describe the level precipitation over time. We are going to consider only a subsample of it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 425
    },
    "colab_type": "code",
    "id": "UAYXOSsBJkx8",
    "outputId": "9cd080e2-6d00-4ddf-88eb-20263e19c498"
   },
   "outputs": [],
   "source": [
    "## REAL-WORLD dataset example\n",
    "## I use the dataset used by Søren Hauberg\n",
    "from scipy.io import loadmat\n",
    "\n",
    "## We load the data and use only some of them.\n",
    "## We subsample the data, which gives us N pairs of (x, y)\n",
    "M = 1000\n",
    "data = loadmat('weather_dataset.mat')\n",
    "# sampling part \n",
    "Xtrain = np.arange(0, M, 20)\n",
    "ytrain = data['TMPMAX'][Xtrain]\n",
    "N = len(ytrain)\n",
    "Xtrain = Xtrain.reshape(-1,1)\n",
    "\n",
    "# print dataset information\n",
    "print('Xtrain shape', Xtrain.shape)\n",
    "print('ytrain shape', ytrain.shape)\n",
    "\n",
    "# also in this case we standardize the data to have zero mean and unit variance\n",
    "Xtrain = (Xtrain - np.mean(Xtrain)) / np.std(Xtrain)\n",
    "ytrain = (ytrain - np.mean(ytrain)) / np.std(ytrain)\n",
    "\n",
    "# and plot it\n",
    "plt.plot(Xtrain, ytrain, 'o', color='black')\n",
    "plt.ylabel('Normalized level of precipitations')\n",
    "plt.xlabel('Time')\n",
    "plt.show()\n",
    "\n",
    "# we shall also define the test set, that is the range of XTest point we want to \n",
    "# use to compute the mean and the variance\n",
    "Xtest = np.linspace(-2, 2, M).reshape(-1,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RJvAiQ2yJkyC"
   },
   "source": [
    "<font color='blue'> Tasks:\n",
    "1. Choose a value for the three hyperparameters and do a noisy Gaussian Porcess regression on this dataset. Plot the mean and the 2-standard deviation interval. Plot also 10 samples form the posterior. If you think that the fit you get is not the perfect one, try to change the values of the hyperparameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 129
    },
    "colab_type": "code",
    "id": "Sir2hGrlJkyF",
    "outputId": "09600333-be2f-4c85-8380-3f603b5a7a45"
   },
   "outputs": [],
   "source": [
    "## as you did above, you should define the hyperparameters of your kernel \n",
    "## and fit a noisy GP on this dataset\n",
    "\n",
    "lengthscale = \n",
    "kernel_variance = \n",
    "noise_variance = \n",
    "\n",
    "n_samples = 10\n",
    "\n",
    "# we start computing the kernels (I am using s instead of * in the kernel variable names)\n",
    "K = \n",
    "Ks = \n",
    "Kss = \n",
    "\n",
    "# compute the cholesky decomposition\n",
    "L = \n",
    "\n",
    "# compute alpha\n",
    "alpha = \n",
    "\n",
    "# compute the mean function\n",
    "mu = \n",
    "\n",
    "# compute v\n",
    "v = \n",
    "\n",
    "# compute the covariance\n",
    "covariance = \n",
    "\n",
    "\n",
    "# we get the standard deviation from the covariance matrix\n",
    "std = \n",
    "\n",
    "\n",
    "## now we can plot the mean function with the standard deviation intervals\n",
    "plt.plot(Xtrain, ytrain, 'ro', label='Training points')\n",
    "plt.gca().fill_between(Xtest.flat, mu.reshape(-1) - 2 * std, mu.reshape(-1) + 2 * std,  color='lightblue', alpha=0.5, label=r\"$\\mu$\")\n",
    "plt.plot(Xtest, mu, 'blue', label=r\"$2\\sigma$\")\n",
    "plt.legend()\n",
    "plt.show() \n",
    "\n",
    "\n",
    "## SAMPLE FROM THE POSTERIOR\n",
    "samples_from_the_posterior = \n",
    "\n",
    "# PLOT THEM\n",
    "plt.plot(Xtrain, ytrain, 'ro', label='Training points')\n",
    "#plt.gca().fill_between(Xtest.flat, mu - 2 * std, mu + 2 * std, color='lightblue', alpha=0.5)\n",
    "#plt.plot(Xtest, mu, 'blue')\n",
    "for sample_id in range(n_samples):\n",
    "    plt.plot(Xtest, samples_from_the_posterior[sample_id], linestyle='-')\n",
    "plt.legend()\n",
    "plt.show()  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ktWobTdHJkyL"
   },
   "source": [
    "<font color='blue'> Tasks:\n",
    "2. Define an interval for the three hyperparameters. Look at you data and find a reasonable range for each of them that is worth to investigate. Use $\\approx 50$ values for each hyperparameter, because otherwise the process will becomes really time consuming. Substitute the values you obtained in the code you wrote before, and check the fitting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 231
    },
    "colab_type": "code",
    "id": "lZna3YozJkyN",
    "outputId": "8e502ac2-e9db-4493-dfbe-da0e2bd3622e"
   },
   "outputs": [],
   "source": [
    "# maximization of the log-likelihood\n",
    "n_of_tests = 50\n",
    "\n",
    "possible_lengthscales = np.linspace(0.001, 2, n_of_tests)\n",
    "possible_output_variances = np.linspace(0.01,3, n_of_tests)\n",
    "possible_noise_variances = np.linspace(0.001, 1, n_of_tests)\n",
    "\n",
    "best_log_likelihood = -100000\n",
    "best_lengthscale = 0\n",
    "best_output_var = 0\n",
    "best_noise_var = 0\n",
    "\n",
    "## simple grid search: FOR LOOPS!\n",
    "for lengthscale in possible_lengthscales:\n",
    "    for output_var in possible_output_variances:\n",
    "        for noise_var in possible_noise_variances:\n",
    "            ### INSERT YOUR CODE HERE\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "            \n",
    "            if log_likelihood > best_log_likelihood:\n",
    "                # update params\n",
    "                best_log_likelihood = log_likelihood\n",
    "                best_lengthscale = lengthscale\n",
    "                best_output_var = output_var\n",
    "                best_noise_var = noise_var\n",
    "                \n",
    "print(\"Max log-likelihood: \", best_log_likelihood)\n",
    "print(\"Best lengthscale: \", best_lengthscale)\n",
    "print(\"Best output variance: \", best_output_var)\n",
    "print(\"Best noise variance: \", best_noise_var)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "GYuKdwNZJkyS"
   },
   "source": [
    "#### Additional exercise: testing extrapolation of Gaussian Process using only a squared-exponential kernel\n",
    "\n",
    "An interesting things to investigate is if the Gaussian Process with squared-exponential kernel is able to extrapolate out of the range of the training set. We will use the hyperparameters that we found out in the previous step and we define the test set to be from -5 to 5.\n",
    "\n",
    "<font color='blue'> Tasks:\n",
    "1. Define your test set to contain $x$s from -5 to 5 and re-fit the Gaussian Process. Plot the mean and 2-standard deviation interval and also some samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 129
    },
    "colab_type": "code",
    "id": "zpXMU18rJkyT",
    "outputId": "9e465183-0e34-4776-f8a9-cf2e86591473"
   },
   "outputs": [],
   "source": [
    "## what happen if we try to model a bigger test set\n",
    "new_Xtest = np.linspace(-5, 5, M).reshape(-1,1)\n",
    "\n",
    "lengthscale = \n",
    "kernel_variance =\n",
    "noise_variance =\n",
    "\n",
    "# we start computing the kernels (I am using s instead of * in the kernel variable names)\n",
    "K = \n",
    "Ks = \n",
    "Kss = \n",
    "\n",
    "# compute the cholesky decomposition\n",
    "L = \n",
    "\n",
    "\n",
    "# compute alpha\n",
    "alpha = \n",
    "\n",
    "# compute the mean function\n",
    "mu = \n",
    "\n",
    "# compute v\n",
    "v = \n",
    "\n",
    "# compute the covariance\n",
    "covariance = \n",
    "\n",
    "# we get the standard deviation from the covariance matrix\n",
    "std = \n",
    "\n",
    "#plot\n",
    "## now we can plot the mean function with the standard deviation intervals\n",
    "plt.plot(Xtrain, ytrain, 'ro', label='Training points')\n",
    "plt.gca().fill_between(new_Xtest.flat, mu.reshape(-1) - 2 * std, mu.reshape(-1) + 2 * std,  color='lightblue', alpha=0.5, label=r\"$\\mu$\")\n",
    "plt.plot(new_Xtest, mu, 'blue', label=r\"$2\\sigma$\")\n",
    "plt.axis([-5, 5, -3, 3])\n",
    "plt.legend()\n",
    "plt.show() \n",
    "\n",
    "\n",
    "## SAMPLE FROM THE POSTERIOR\n",
    "samples_from_the_posterior = \n",
    "\n",
    "\n",
    "# PLOT THEM\n",
    "plt.plot(Xtrain, ytrain, 'ro', label='Training points')\n",
    "#plt.gca().fill_between(Xtest.flat, mu - 2 * var, mu + 2 * var, color='lightblue', alpha=0.5)\n",
    "#plt.plot(Xtest, mu, 'blue')\n",
    "for sample_id in range(n_samples):\n",
    "    plt.plot(new_Xtest, samples_from_the_posterior[sample_id], linestyle='-')\n",
    "plt.axis([-5, 5, -3, 3])\n",
    "plt.legend()\n",
    "plt.show()  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "kBDOaYPjJkyZ"
   },
   "source": [
    "## Optimization using Gradient Descent\n",
    "\n",
    "Another possibility is to optimize the three hyperparameters using Gradient Descent. Here, there are several options in order to calculate the gradient:\n",
    "1. Do it manually like in the good old days, here checking if you did it correctly by a finite difference approximation would properly be good\n",
    "2. Use automatic diferentiation for example in PyTorch or TensorFlow<br>\n",
    "To use PyTorch, we would have to compute all the operations that involves our hyperparameters using PyTorch, there is a bit of skeleton code for this below.\n",
    "3. Just ignore it and have scipy.optimize automatically estimate the gradient by finite difference a good algorithm for this is L-BFGS.\n",
    "\n",
    "<font color='blue'> Tasks:\n",
    "1. Implement gradient descent for optimizing kernel hyperparameters, by either of the approaches above\n",
    "\n",
    "Experiment a bit with finding hyperparameters do you encounter any problems? - if so what appears to be happening and can you do something to circumvent the problem?\n",
    "</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Oggph6IXJkyb",
    "outputId": "0f740c02-f750-4645-9442-b9ea93c27086"
   },
   "outputs": [],
   "source": [
    "import torch\n",
    "from torch import nn\n",
    "\n",
    "torch.set_default_tensor_type(torch.FloatTensor)\n",
    "\n",
    "def squared_exponential_kernel_torch(x, y, lengthscale, variance):\n",
    "    x = x.squeeze(1).expand(x.size(0), y.size(0))\n",
    "    y = y.squeeze(0).expand(x.size(0), y.size(0))\n",
    "    sqdist = torch.pow(x-y, 2)\n",
    "    k = variance * torch.exp(-0.5 * sqdist * (1/lengthscale**2))  # NxM\n",
    "    return k\n",
    "\n",
    "# tranform our training set in Tensor\n",
    "Xtrain_tensor = torch.from_numpy(Xtrain).float()\n",
    "ytrain_tensor = torch.from_numpy(ytrain).float()\n",
    "\n",
    "\n",
    "### YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "yWOw4j2BJkyx"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "colab": {
   "collapsed_sections": [],
   "name": "Exercise2_solutions.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
}