{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gaussian Process\n",
    "\n",
    "In this tutorial, we expose what gaussian processes are, and how to use the [GPy library](http://sheffieldml.github.io/GPy/). We first provide a gentle reminder about Gaussian distributions and their properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import all the important libraries\n",
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm\n",
    "from scipy.stats import multivariate_normal\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.mlab as mlab\n",
    "from ipywidgets import widgets as wg\n",
    "from matplotlib import cm\n",
    "\n",
    "import GPy\n",
    "#%matplotlib inline\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1D Gaussian distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "def plot_gaussian(mu=0, sigma=1):\n",
    "    x = np.linspace(-3, 3, 100)\n",
    "    plt.plot(x, norm.pdf(x, mu, sigma))\n",
    "    plt.xlabel('x')\n",
    "    plt.ylabel('p(x)')\n",
    "\n",
    "wg.interact(plot_gaussian, mu=(-2,2,0.1), sigma=(-2,2,0.1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multivariate Gaussian distribution (2D)\n",
    "\n",
    "The multivariable Gaussian distribution is a generalization of the Gaussian distribution to vectors. See [wikipedia](https://en.wikipedia.org/wiki/Multivariate_normal_distribution) for more info."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# moments\n",
    "mu = np.array([0,0])\n",
    "Sigma = np.array([[1,0], \n",
    "                  [0,1]])\n",
    "Sigma1 = np.array([[1,0.5],\n",
    "                   [0.5,1]])\n",
    "Sigma2 = np.array([[1,-0.5],\n",
    "                   [-0.5,1]])\n",
    "Sigmas = [Sigma, Sigma1, Sigma2]\n",
    "\n",
    "pts = []\n",
    "for S in Sigmas:\n",
    "    pts.append(np.random.multivariate_normal(mu, S, 1000).T)\n",
    "\n",
    "# Plotting\n",
    "width = 16\n",
    "height = 4\n",
    "plt.figure(figsize=(width, height))\n",
    "\n",
    "# make plot\n",
    "for i in range(len(Sigmas)):\n",
    "    plt.subplot(1,3,i+1)\n",
    "    plt.title('Plot '+str(i+1))\n",
    "    plt.ylim(-4,4)\n",
    "    plt.xlim(-4,4)\n",
    "    plt.xlabel('x1')\n",
    "    plt.ylabel('x2')\n",
    "    plt.plot(pts[i][0], pts[i][1],'o')\n",
    "    #plt.scatter(pts[i][0], pts[i][1])\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 1st plot above is described by:\n",
    "\\begin{equation}\n",
    "    \\left[ \\begin{array}{c} x_1\\\\x_2 \\end{array}\\right] \\sim \\mathcal{N} \\left(\\left[ \\begin{array}{c} 0\\\\0 \\end{array}\\right], \\left[ \\begin{array}{cc} 1 & 0\\\\ 0 & 1 \\end{array}\\right] \\right)\n",
    "\\end{equation}\n",
    "\n",
    "The 2nd plot is given by:\n",
    "\\begin{equation}\n",
    "     \\left[ \\begin{array}{c} x_1\\\\x_2 \\end{array}\\right] \\sim \\mathcal{N}\\left(\\left[ \\begin{array}{c} 0\\\\0 \\end{array}\\right], \\left[ \\begin{array}{cc} 1 & 0.5\\\\ 0.5 & 1 \\end{array}\\right]\\right)\n",
    "\\end{equation}\n",
    "\n",
    "Finally, the 3rd plot is given by:\n",
    "\\begin{equation}\n",
    "    \\left[ \\begin{array}{c} x_1\\\\x_2 \\end{array}\\right] \\sim \\mathcal{N}\\left(\\left[ \\begin{array}{c} 0\\\\0 \\end{array}\\right], \\left[ \\begin{array}{cc} 1 & -0.5\\\\ -0.5 & 1 \\end{array}\\right]\\right)\n",
    "\\end{equation}\n",
    "\n",
    "The covariance (and the dot product) measures the similarity.\n",
    "\n",
    "For the 2nd and 3rd plots, $x_1$ is **correlated** with $x_2$, i.e. knowing $x_1$ gives us information about $x_2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Joint distribution $p(x_1,x_2)$\n",
    "\n",
    "The joint distribution $p(x_1, x_2)$ is given by:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\left[ \\begin{array}{c} x_1\\\\x_2 \\end{array}\\right] \\sim \\mathcal{N}\\left(\\left[ \\begin{array}{c} \\mu_1 \\\\ \\mu_2 \\end{array}\\right], \\left[ \\begin{array}{cc} \\Sigma_{11} & \\Sigma_{12} \\\\ \\Sigma_{21} & \\Sigma_{22} \\end{array}\\right] \\right) = \\mathcal{N}(\\pmb{\\mu}, \\pmb{\\Sigma})\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reference: http://stackoverflow.com/questions/38698277/plot-normal-distribution-in-3d\n",
    "\n",
    "# moments\n",
    "mu = np.array([0,0])\n",
    "Sigma = np.array([[1,0], [0,1]])\n",
    "\n",
    "# Create grid and multivariate normal\n",
    "step = 500\n",
    "bound = 10\n",
    "x = np.linspace(-bound,bound,step)\n",
    "y = np.linspace(-bound,bound,step)\n",
    "X, Y = np.meshgrid(x,y)\n",
    "pos = np.empty(X.shape + (2,))\n",
    "pos[:, :, 0] = X; pos[:, :, 1] = Y\n",
    "pdf = multivariate_normal(mu, Sigma).pdf(pos)\n",
    "\n",
    "# Plot\n",
    "fig = plt.figure(figsize=plt.figaspect(0.5)) # Twice as wide as it is tall.\n",
    "\n",
    "# 1st subplot (3D)\n",
    "ax = fig.add_subplot(1, 2, 1, projection='3d')\n",
    "ax.plot_surface(X, Y, pdf, cmap='viridis', linewidth=0)\n",
    "ax.set_xlabel('x1')\n",
    "ax.set_ylabel('x2')\n",
    "ax.set_zlabel('p(x1, x2)')\n",
    "\n",
    "# 2nd subplot (2D)\n",
    "ax = fig.add_subplot(1, 2, 2)\n",
    "ax.contourf(x, y, pdf)\n",
    "#ax.colorbar()\n",
    "ax.set_xlabel('x1')\n",
    "ax.set_ylabel('x2')\n",
    "\n",
    "fig.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalization\n",
    "\n",
    "In order to be a valid probability distribution, the volume under the surface should equal to 1.\n",
    "\n",
    "\\begin{equation}\n",
    "    \\int \\int p(x_1,x_2) dx_1 dx_2 = 1\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# p(x1,x2) = pdf\n",
    "# dx = 2.*bound/step\n",
    "# dx1 dx2 = (2.*bound/step)**2\n",
    "print(\"Summation: {}\".format((2.*bound/step)**2 * pdf.sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conditional distribution $p(x_2|x_1)$\n",
    "\n",
    "What is the mean $\\mu_{2|1}$ and the variance $\\Sigma_{2|1}$ of the conditional distribution $p(x_2|x_1) = \\mathcal{N}(\\mu_{2|1}, \\Sigma_{2|1})$?\n",
    "\n",
    "We know the mean $\\pmb{\\mu}$ and the covariance $\\pmb{\\Sigma}$ of the joint distribution $p(x_1,x_2)$. Using the [Schur complement](https://en.wikipedia.org/wiki/Schur_complement), we obtain:\n",
    "\n",
    "\\begin{align}\n",
    "    \\mu_{2|1} &= \\mu_{2} + \\Sigma_{21}\\Sigma_{22}^{-1}(x_2 - \\mu_2) \\\\\n",
    "    \\Sigma_{2|1} &= \\Sigma_{22} - \\Sigma_{21}\\Sigma_{22}^{-1}\\Sigma_{12}\n",
    "\\end{align}\n",
    "\n",
    "For the demo, check Murphy's book \"Machine Learning: A Probabilistic Perspective\", section 4.3.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=plt.figaspect(0.5)) # Twice as wide as it is tall.\n",
    "\n",
    "x1_value = 0\n",
    "z_max = pdf.max()\n",
    "\n",
    "# 1st subplot\n",
    "ax = fig.add_subplot(1, 2, 1, projection='3d')\n",
    "ax.plot_surface(X, Y, pdf, cmap='viridis', linewidth=0)\n",
    "ax.set_xlabel('x1')\n",
    "ax.set_ylabel('x2')\n",
    "ax.set_zlabel('p(x1, x2)')\n",
    "y1 = np.linspace(-bound,bound,2)\n",
    "z = np.linspace(0,z_max,2)\n",
    "Y1, Z = np.meshgrid(y1,z)\n",
    "ax.plot_surface(x1_value, Y1, Z, color='red', alpha=0.2)\n",
    "#cset = ax.contourf(X, Y, pdf, zdir='x', offset=-bound, cmap=cm.coolwarm)\n",
    "\n",
    "# 2nd subplot\n",
    "ax = fig.add_subplot(1, 2, 2)\n",
    "ax.plot(x, pdf[step//2 + x1_value*step//(2*bound)])\n",
    "ax.set_xlabel('x2')\n",
    "ax.set_ylabel('p(x2|x1)')\n",
    "\n",
    "fig.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Marginal distribution $p(x_1)$ and $p(x_2)$\n",
    "\n",
    "\\begin{align}\n",
    "    p(x_1) &= \\int p(x_1, x_2) dx_2 = \\mathcal{N}(\\mu_1, \\Sigma_{11}) \\\\\n",
    "    p(x_2) &= \\int p(x_1, x_2) dx_1 = \\mathcal{N}(\\mu_2, \\Sigma_{22})\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=plt.figaspect(0.5))\n",
    "plt.subplot(1,2,1)\n",
    "plt.title('By summing')\n",
    "dx = 2. * bound / step\n",
    "plt.plot(x, pdf.sum(0) * dx, color='blue')\n",
    "plt.xlabel('x2')\n",
    "plt.ylabel('p(x2)')\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.title('by using the normal distribution')\n",
    "plt.plot(x, norm.pdf(x, mu[1], Sigma[1,1]), color='red')\n",
    "plt.xlabel('x2')\n",
    "plt.ylabel('p(x2)')\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian Processes (GPs)\n",
    "\n",
    "A Gaussian process is a Gaussian distribution over functions. That is, it is a generalization of the multivariable Gaussian distribution to infinite vectors.\n",
    "\n",
    "It will become clearer with an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [0.5,0.8,1.4]\n",
    "f = [1,2,6]\n",
    "\n",
    "plt.plot(x,f,'o')\n",
    "for i in range(len(x)):\n",
    "    plt.annotate('f'+str(i+1), (x[i],f[i]))\n",
    "plt.xlim(0,2)\n",
    "plt.ylim(0,6.5)\n",
    "plt.ylabel('f(x)')\n",
    "plt.xlabel('x')\n",
    "plt.xticks(x, ['x'+str(i+1) for i in range(len(x))])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "    \\left[ \\begin{array}{c} f_1 \\\\ f_2 \\\\ f_3 \\end{array}\\right] \n",
    "    &\\sim \\mathcal{N}\\left( \\left[ \\begin{array}{c} 0 \\\\ 0 \\\\ 0 \\end{array}\\right], \\left[ \\begin{array}{ccc} K_{11} & K_{12} & K_{13} \\\\ K_{21} & K_{22} & K_{23} \\\\ K_{31} & K_{32} & K_{33} \\end{array}\\right] \\right) \\\\\n",
    "    &\\sim \\mathcal{N}\\left( \\left[ \\begin{array}{c} 0 \\\\ 0 \\\\ 0 \\end{array}\\right], \\left[ \\begin{array}{ccc} 1 & 0.7 & 0.2 \\\\ 0.7 & 1 & 0.6 \\\\ 0.2 & 0.6 & 1 \\end{array}\\right] \\right)\n",
    "\\end{align}\n",
    "\n",
    "Similarity measure: $K_{ij} = \\exp(- ||x_i - x_j||^2) = \\left\\{ \\begin{array}{ll} 0 & ||x_i - x_j|| \\rightarrow \\infty \\\\ 1 & x_i = x_j \\end{array} \\right.$\n",
    "\n",
    "Prediction (noiseless GP regression): given data $\\mathcal{D} = \\{(x_1,f_1), (x_2,f_2), (x_3,f_3)\\}$, and new point $x_*$ (e.g. $x_*$=1.4), what is the value of $f_*$?\n",
    "\n",
    "\\begin{equation}\n",
    "    \\pmb{f} \\sim \\mathcal{N}(\\pmb{0}, \\pmb{K}) \\qquad \\mbox{and} \\qquad f_* \\sim \\mathcal{N}(0, K(x_*,x_*)) = \\mathcal{N}(0, K_{**})\n",
    "\\end{equation}\n",
    "\n",
    "In this case, $K_{**} = K(x_*,x_*) = \\exp(- ||x_* - x_*||^2) = 1$.\n",
    "\n",
    "Now, we can write:\n",
    "\n",
    "\\begin{equation}\n",
    "        \\left[ \\begin{array}{c} \\pmb{f} \\\\ f_* \\end{array} \\right] \\sim \\mathcal{N}\\left(\\pmb{0}, \\left[\\begin{array}{cc} \\left[ \\begin{array}{ccc} K_{11} & K_{12} & K_{13} \\\\ K_{21} & K_{22} & K_{23} \\\\ K_{31} & K_{32} & K_{33} \\end{array} \\right] & \\left[ \\begin{array}{c} K_{1*} \\\\ K_{2*} \\\\ K_{3*} \\end{array} \\right] \\\\ \\left[ \\begin{array}{ccc} K_{*1} & K_{*2} & K_{*3} \\end{array} \\right] & \\left[\\begin{array}{c} K_{**} \\end{array} \\right] \\end{array} \\right] \\right) = \\mathcal{N}\\left(\\pmb{0}, \\left[\\begin{array}{cc} \\pmb{K} & \\pmb{K}_* \\\\ \\pmb{K}_*^T & \\pmb{K}_{**} \\end{array}\\right]\\right)\n",
    "\\end{equation}\n",
    "\n",
    "Using the formula for the conditional probability $p(f_*|f)$, we have:\n",
    "\n",
    "\\begin{align}\n",
    "    \\mu_* &= \\mathbb{E}[f_*] = \\pmb{K}_*^T \\pmb{K}^{-1}\\pmb{f} \\\\\n",
    "    c_* &= K_{**} - \\pmb{K}_*^T \\pmb{K}^{-1}\\pmb{K}_*\n",
    "\\end{align}\n",
    "\n",
    "We can thus predict the mean $\\mu_*$ and the variance $c_*$ for the test point $x_*$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [0.5,0.8,1.4]\n",
    "f = [1,2,6]\n",
    "\n",
    "x_new = 1.3\n",
    "f_new = 5.2\n",
    "\n",
    "plt.plot(x+[x_new],f+[f_new],'o')\n",
    "for i in range(len(x)):\n",
    "    plt.annotate('f'+str(i+1), (x[i],f[i]))\n",
    "plt.errorbar(x_new, f_new, yerr=1)\n",
    "plt.annotate('f*', (x_new+0.02, f_new))\n",
    "plt.xlim(0,2)\n",
    "plt.ylim(0,6.5)\n",
    "plt.ylabel('f(x)')\n",
    "plt.xlabel('x')\n",
    "plt.xticks(x+[x_new], ['x'+str(i+1) for i in range(len(x))]+['x*'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generalization\n",
    "\n",
    "A GP defines a distribution over functions $p(f)$ (i.e. it is the joint distribution over all the infinite function values).\n",
    "\n",
    "Definition: $p(f)$ is a GP if for any finite subset $\\{x_1,...,x_n\\} ⊂ X$, the marginal distribution over that finite subset $p(f)$ has a multivariate Gaussian distribution.\n",
    "\n",
    "Prior on $f$:\n",
    "\\begin{equation}\n",
    "    \\pmb{f}|\\pmb{x} \\sim \\mathcal{GP}(\\pmb{\\mu}(\\pmb{x}), \\pmb{K}(\\pmb{x}, \\pmb{x}))\n",
    "\\end{equation}\n",
    "with\n",
    "\\begin{align*}\n",
    "    \\pmb{\\mu}(\\pmb{x}) &= \\mathbb{E}_f \\lbrack \\pmb{x} \\rbrack \\\\\n",
    "    k(\\pmb{x}, \\pmb{x'}) &= \\mathbb{E}_f \\lbrack (\\pmb{x} - \\pmb{\\mu}(\\pmb{x})) (\\pmb{x'} - \\pmb{\\mu}(\\pmb{x'})) \\rbrack\n",
    "\\end{align*}\n",
    "\n",
    "Often written as:\n",
    "\\begin{equation}\n",
    "    \\pmb{f} \\sim \\mathcal{GP}(\\pmb{0}, \\pmb{K})\n",
    "\\end{equation}\n",
    "\n",
    "Concretely, assume $\\pmb{x} \\in \\mathbb{R}^{50}$, then $\\pmb{K}(\\pmb{x}, \\pmb{x}) \\in \\mathbb{R}^{50 \\times 50}$, then $\\pmb{f} \\sim \\mathcal{GP}(\\pmb{0}, \\pmb{K})$ means:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\left[ \\begin{array}{c} f_1 \\\\ \\vdots \\\\ f_{50} \\end{array}\\right] := \\left[ \\begin{array}{c} f(x_1) \\\\ \\vdots \\\\ f(x_{50}) \\end{array}\\right]  \\sim \\mathcal{N}\\left( \\left[ \\begin{array}{c} 0 \\\\ \\vdots \\\\ 0 \\end{array}\\right], \\left[ \\begin{array}{ccc} k(x_1,x_1) & \\cdots & k(x_1, x_{50}) \\\\ \\vdots & \\ddots & \\vdots \\\\ k(x_{50},x_1) & \\cdots & k(x_{50}, x_{50}) \\end{array} \\right] \\right)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### RBF kernel\n",
    "\n",
    "Let's choose a RBF (a.k.a Squared Exponential, Gaussian) kernel:\n",
    "\n",
    "\\begin{equation}\n",
    "\\pmb{K} = \\left[ \\begin{array}{ccc} k(x_1,x_1) & \\cdots & k(x_1, x_d) \\\\ \\vdots & \\ddots & \\vdots \\\\ k(x_d,x_1) & \\cdots & k(x_d, x_d) \\end{array} \\right]\n",
    "\\end{equation}\n",
    "with\n",
    "\\begin{equation}\n",
    "    k(x_i, x_j) = \\alpha^2 \\exp \\left( - \\frac{(x_i - x_j)^2}{2l} \\right) \\qquad \\mbox{ and hyperparameters } \\pmb{\\Phi} =  \\left\\{ \\begin{array}{l} \\alpha \\mbox{: amplitude} \\\\ l \\mbox{: the lengthscale} \\end{array} \\right.\n",
    "\\end{equation}\n",
    "\n",
    "This function $k$ is infinitely differentiable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reference: https://www.youtube.com/watch?v=4vGiHC35j9s&t=51s\n",
    "\n",
    "# Hyperparameters\n",
    "alpha = 1\n",
    "l = 2\n",
    "\n",
    "# Parameters\n",
    "n = 50 # nb of points\n",
    "n_func = 10 # nb of fct to draw\n",
    "x_bound = 5 # bound on the x axis\n",
    "\n",
    "def RBF_kernel(a,b):\n",
    "    sqdist = np.sum(a**2,1).reshape(-1,1) + np.sum(b**2,1) - 2*np.dot(a,b.T)\n",
    "    return alpha**2 * np.exp(-1/l * sqdist)\n",
    "\n",
    "n = 50\n",
    "X = np.linspace(-x_bound, x_bound, n).reshape(-1,1)\n",
    "K = RBF_kernel(X, X) # dim(K) = n x n\n",
    "\n",
    "L = np.linalg.cholesky(K + 1e-6 * np.eye(n))\n",
    "f_prior = np.dot(L, np.random.normal(size=(n, n_func)))\n",
    "\n",
    "# Plotting\n",
    "width = 16\n",
    "height = 4\n",
    "plt.figure(figsize=(width, height))\n",
    "\n",
    "# plot f_prior\n",
    "plt.subplot(1,3,1)\n",
    "plt.title('GP: prior on f')\n",
    "plt.plot(X, f_prior)\n",
    "plt.plot(X, f_prior.mean(1), linewidth=3, color='black')\n",
    "plt.ylabel('f(x)')\n",
    "plt.xlabel('x')\n",
    "\n",
    "# plot Kernel\n",
    "plt.subplot(1,3,2)\n",
    "plt.title('Kernel matrix')\n",
    "plt.pcolor(K[::-1])\n",
    "plt.colorbar()\n",
    "\n",
    "plt.subplot(1,3,3)\n",
    "plt.title('Kernel function')\n",
    "plt.plot(X, RBF_kernel(X, np.array([[1.0]])))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kernel (prior knowledge)\n",
    "\n",
    "By choosing a specific kernel, we can incorporate prior knowledge that we have about the function $f$, such as, if the function is:\n",
    "* periodic\n",
    "* smooth\n",
    "* symmetric\n",
    "* etc.\n",
    "\n",
    "The hyperparameters for each kernel are also very intuitive/interpretable.\n",
    "\n",
    "Note: kernels can be combined!\n",
    "\n",
    "Indeed, if $k(x,y)$, $k_1(x,y)$ and $k_2(x,y)$ are valid kernels then:\n",
    "* $\\alpha k(x,y) $ with $\\alpha \\geq 0$\n",
    "* $k_1(x,y) + k_2(x,y)$\n",
    "* $k_1(x,y) k_2(x,y)$\n",
    "* $p(k(x,y))$ with $p$ being a polynomial function with non-negative coefficients\n",
    "* $exp(k(x,y))$\n",
    "* $f(x) k(x,y) \\overline{f(y)}$ with $\\overline{f} = $ complex conjugate\n",
    "* $k(\\phi(x),\\phi(y))$\n",
    "\n",
    "are all valid kernels!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Periodic Exponential kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variance = 1.\n",
    "lengthscale = 1.\n",
    "period = 2.*np.pi\n",
    "\n",
    "#K = periodic_kernel(X, X) # dim(K) = n x n\n",
    "kern = GPy.kern.PeriodicExponential(variance=variance, lengthscale=lengthscale, period=period)\n",
    "K1 = kern.K(X)\n",
    "\n",
    "L = np.linalg.cholesky(K1 + 1e-6 * np.eye(n))\n",
    "f_prior = np.dot(L, np.random.normal(size=(n, 1)))\n",
    "\n",
    "# Plotting\n",
    "width = 16\n",
    "height = 4\n",
    "plt.figure(figsize=(width, height))\n",
    "\n",
    "# plot f_prior\n",
    "plt.subplot(1,3,1)\n",
    "plt.title('GP: prior on f')\n",
    "plt.plot(X, f_prior)\n",
    "plt.plot(X, f_prior.mean(1), linewidth=3, color='black')\n",
    "plt.ylabel('f(x)')\n",
    "plt.xlabel('x')\n",
    "\n",
    "# plot Kernel\n",
    "plt.subplot(1,3,2)\n",
    "plt.title('Kernel matrix')\n",
    "plt.pcolor(K1[::-1])\n",
    "plt.colorbar()\n",
    "\n",
    "plt.subplot(1,3,3)\n",
    "plt.title('Kernel function')\n",
    "plt.plot(X, kern.K(X, np.array([[1.0]])))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### addition and multiplication of 2 kernels (SE and PE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "K_add = K + K1\n",
    "\n",
    "L = np.linalg.cholesky(K_add + 1e-6 * np.eye(n))\n",
    "f_prior = np.dot(L, np.random.normal(size=(n, n_func)))\n",
    "\n",
    "# Plotting\n",
    "width = 16\n",
    "height = 8\n",
    "plt.figure(figsize=(width, height))\n",
    "\n",
    "# plot f_prior\n",
    "plt.subplot(2,2,1)\n",
    "plt.title('GP: prior on f with K_add')\n",
    "plt.plot(X, f_prior)\n",
    "plt.plot(X, f_prior.mean(1), linewidth=3, color='black')\n",
    "plt.ylabel('f(x)')\n",
    "plt.xlabel('x')\n",
    "\n",
    "# plot Kernel\n",
    "plt.subplot(2,2,2)\n",
    "plt.title('Kernel matrix: K_add')\n",
    "plt.pcolor(K_add[::-1])\n",
    "plt.colorbar()\n",
    "\n",
    "K_prod = K * K1\n",
    "\n",
    "L = np.linalg.cholesky(K_prod + 1e-6 * np.eye(n))\n",
    "f_prior = np.dot(L, np.random.normal(size=(n, n_func)))\n",
    "\n",
    "# plot f_prior\n",
    "plt.subplot(2,2,3)\n",
    "plt.title('GP: prior on f with K_prod')\n",
    "plt.plot(X, f_prior)\n",
    "plt.plot(X, f_prior.mean(1), linewidth=3, color='black')\n",
    "plt.ylabel('f(x)')\n",
    "plt.xlabel('x')\n",
    "\n",
    "# plot Kernel\n",
    "plt.subplot(2,2,4)\n",
    "plt.title('Kernel matrix: K_prod')\n",
    "plt.pcolor(K_prod[::-1])\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GP Posterior\n",
    "\n",
    "Given $\\mathcal{D}=\\{(x_i, y_i)\\}_{i=1}^{i=N} = (\\pmb{X}, \\pmb{y})$, we have:\n",
    "\n",
    "\\begin{equation}\n",
    "    p(f|\\mathcal{D}) = \\frac{p(\\mathcal{D}|f)p(f)}{p(\\mathcal{D})}\n",
    "\\end{equation}\n",
    "\n",
    "### GP Regression\n",
    "\n",
    "\\begin{equation}\n",
    "    y_i = f(\\pmb{x}_i) + \\epsilon_i \\qquad \n",
    "\t\\left\\{ \\begin{array}{l}\n",
    "\t\tf \\sim \\mathcal{GP}(\\pmb{0}, \\pmb{K}) \\\\\n",
    "\t\t\\epsilon_i \\sim \\mathcal{N}(0, \\sigma^2)\n",
    "\t\\end{array} \\right.\n",
    "\\end{equation}\n",
    "\n",
    "* Prior $f$ is a GP $\\Leftrightarrow p(\\pmb{f}|\\pmb{X}) = \\mathcal{N}(\\pmb{0}, \\pmb{K})$\n",
    "* Likelihood is Gaussian $\\Leftrightarrow p(\\pmb{y}|\\pmb{X},\\pmb{f}) = \\mathcal{N}(\\pmb{f}, \\sigma^2\\pmb{I})$\n",
    "* $\\rightarrow p(f|\\mathcal{D})$ is also a GP.\n",
    "\n",
    "#### Predictive distribution: \n",
    "$$p(\\pmb{y}_*|\\pmb{x}_*, \\pmb{X}, \\pmb{y}) = \\int p(\\pmb{y}_{*}| \\pmb{x}_{*}, \\pmb{f}, \\pmb{X}, \\pmb{y}) p(f|\\pmb{X}, \\pmb{y}) d\\pmb{f} = \\mathcal{N}(\\pmb{\\mu}_*, \\pmb{\\Sigma}_*)$$\n",
    "\\begin{align}\n",
    "    \\pmb{\\mu}_* &= \\pmb{K}_{*N} (\\pmb{K}_N + \\sigma^2 \\pmb{I})^{-1} \\pmb{y} \\\\\n",
    "    \\pmb{\\Sigma}_* &= \\pmb{K}_{**} - \\pmb{K}_{*N} (\\pmb{K}_N + \\sigma^2 \\pmb{I})^{-1} \\pmb{K}_{N*}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning a GP\n",
    "#### Marginal likelihood:\n",
    "\\begin{equation}\n",
    "    p(\\pmb{y}|\\pmb{X}) = \\int p(\\pmb{y}|\\pmb{f},\\pmb{X}) p(\\pmb{f}|\\pmb{X}) d\\pmb{f} = \\mathcal{N}(\\pmb{0}, \\pmb{K} + \\sigma^2\\pmb{I})\n",
    "\\end{equation}\n",
    "\n",
    "By taking the logarithm, and setting $\\pmb{K}_y = (\\pmb{K} + \\sigma^2\\pmb{I})$, we have:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\mathcal{L} = \\log p(\\pmb{y}|\\pmb{X}; \\pmb{\\Phi}) = \\underbrace{-\\frac{1}{2} \\pmb{y}^T \\pmb{K}_y^{-1} \\pmb{y}}_{\\mbox{data fit}} \\underbrace{-\\frac{1}{2} \\log |\\pmb{K}_y^{-1}|}_{\\mbox{complexity penalty}} - \\frac{n}{2} \\log 2\\pi\n",
    "\\end{equation}\n",
    "\n",
    "The marginal likelihood (i.e. ML-II) is used to optimize the hyperparameters $\\pmb{\\Phi}$ that defines the covariance function and thus the GP.\n",
    "\n",
    "\\begin{equation}\n",
    "    \\pmb{\\Phi}^* = argmax_{\\pmb{\\Phi}} \\log p(\\pmb{y}|\\pmb{X}; \\pmb{\\Phi})\n",
    "\\end{equation}\n",
    "\n",
    "Optimizing the marginal likelihood is more robust than the likelihood as it tries to optimize the complexity of the model, and the fitting of this last one to the observed data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GP Regression\n",
    "# Based on the tutorial: https://github.com/SheffieldML/notebook/blob/master/GPy/GPyCrashCourse.ipynb\n",
    "\n",
    "# Create dataset\n",
    "X = np.random.uniform(-3.0, 3.0, (20,1))\n",
    "Y = np.sin(X) + np.random.randn(20,1) * 0.05 \n",
    "\n",
    "# Create the kernel\n",
    "# Reminder 1: The sum of valid kernels gives a valid kernel.\n",
    "# Reminder 2: The product of valid kernels gives a valid kernel.\n",
    "# Available kernels: RBF, Exponential, Matern32, Matern52, Brownian, Bias, Linear, PeriodicExponential, White.\n",
    "kernel = GPy.kern.RBF(input_dim=1, variance=1.0, lengthscale=1.0)\n",
    "\n",
    "# Create the model\n",
    "gp_model = GPy.models.GPRegression(X, Y, kernel)\n",
    "\n",
    "# Display and plot\n",
    "print(\"Before optimization: \", gp_model)\n",
    "gp_model.plot()\n",
    "plt.show()\n",
    "\n",
    "# Optimize the model (that is find the 'best' hyperparameters of the kernel matrix)\n",
    "# By default, the optimizer is a 2nd order algo: lbfgsb. Others are available such as the scg, ...\n",
    "gp_model.optimize(messages=False)\n",
    "\n",
    "# Display and plot\n",
    "print(\"After optimization: \", gp_model)\n",
    "gp_model.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian Process Latent Variable Model (GP-LVM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GPLVM\n",
    "# Based on the tutorials: \n",
    "# http://nbviewer.jupyter.org/github/SheffieldML/notebook/blob/master/GPy/MagnificationFactor.ipynb\n",
    "# https://github.com/SheffieldML/notebook/blob/master/lab_classes/gprs/lab4-Copy0.ipynb\n",
    "\n",
    "# Create dataset\n",
    "N = 100\n",
    "k1 = GPy.kern.RBF(5, variance=1, lengthscale=1./np.random.dirichlet(np.r_[10,10,10,0.1,0.1]), ARD=True)\n",
    "k2 = GPy.kern.RBF(5, variance=1, lengthscale=1./np.random.dirichlet(np.r_[0.1,10,10,10,0.1]), ARD=True)\n",
    "X = np.random.normal(0, 1, (N,5))\n",
    "A = np.random.multivariate_normal(np.zeros(N), k1.K(X), 10).T\n",
    "B = np.random.multivariate_normal(np.zeros(N), k2.K(X), 10).T\n",
    "\n",
    "Y = np.vstack((A,B))\n",
    "\n",
    "# latent space dimension\n",
    "latent_dim = 2\n",
    "\n",
    "# Create the kernel\n",
    "kernel = GPy.kern.RBF(input_dim=latent_dim, variance=1.0, lengthscale=1.0)\n",
    "\n",
    "# Create the GPLVM model\n",
    "gplvm_model = GPy.models.GPLVM(Y, latent_dim, init='PCA', kernel=kernel)\n",
    "\n",
    "# Display and plot\n",
    "print(\"Before optimization: \", gplvm_model)\n",
    "gplvm_model.plot_latent()\n",
    "plt.show()\n",
    "\n",
    "# Optimize the model (that is find the 'best' hyperparameters of the kernel matrix)\n",
    "# By default, the optimizer is a 2nd order algo: lbfgsb. Others are available such as the scg, ...\n",
    "gplvm_model.optimize(messages=False)\n",
    "\n",
    "# Display and plot\n",
    "print(\"After optimization: \", gplvm_model)\n",
    "gplvm_model.plot_latent()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}