gaussian_process_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook copied from https://github.com/ChakriCherukuri/mlviz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Multi-variate Gaussian Distribution </h2>\n",
    "\n",
    "The multivariate normal distribution of a k-dimensional random vector $ \\large X =(X_{1},\\ldots ,X_{k})^{T}$ is given by:<br>\n",
    "$ \\large  X \\sim  N \\boldsymbol (\\mu, \\boldsymbol \\Sigma) $\n",
    "with k-dimensional mean vector <br>\n",
    "$ \\large \\boldsymbol \\mu = E[X] = [E[X_{1}], E[X_{2}],\\ldots , E[X_{k}]]^{T} $\n",
    "\n",
    "and $ \\large k \\times k $ covariance matrix $\\large \\Sigma $\n",
    "\n",
    "<h2> Gaussian Process </h2>\n",
    "* GP is Gaussian distribution over functions. We are learning a *distribution* from empirical data\n",
    "* Start with prior $\\large f_{prior} \\sim  N  (\\boldsymbol 0, \\boldsymbol K) $\n",
    "* Compute the posterior from training samples by doing a bayesian update<br>\n",
    "$ \\large  f_{posterior} \\sim  N \\boldsymbol (\\mu_{*}, \\Sigma_{*}) $ where <br>\n",
    "$ \\large \\mu_* = K_* K^{-1} y $ <br>\n",
    "$ \\large \\Sigma_* = K_{∗∗} − K_∗ K^{−1} K^T $\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from ipywidgets import *\n",
    "import bqplot.pyplot as plt\n",
    "from bqplot import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def squared_exponential(x1, x2, sigma=1., l=1.):\n",
    "    z = (x1 - x2[:, np.newaxis]) / l\n",
    "    return sigma**2 * np.exp(-.5 * z ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gp_regression(X_train, y_train, X_test,\n",
    "                  kernel=squared_exponential,\n",
    "                  sigma_noise=.1,\n",
    "                  params=dict(sigma=1., l=1.)):\n",
    "    # compute the kernel matrices for train, train_test, test combinations\n",
    "    K = kernel(X_train, X_train, **params)\n",
    "    K_s = kernel(X_train, X_test, **params)\n",
    "    K_ss = kernel(X_test, X_test, **params)\n",
    "    \n",
    "    n, p = len(X_train), len(X_test)\n",
    "    \n",
    "    # compute the posterior mean and cov\n",
    "    mu_s = np.dot(K_s, np.linalg.solve(K + sigma_noise**2 * np.eye(n), y_train))\n",
    "    cov_s = K_ss - np.dot(K_s, np.linalg.solve(K + sigma_noise**2 * np.eye(n), K_s.T))\n",
    "    \n",
    "    # prior and posterior moments\n",
    "    mu_prior, cov_prior = np.zeros(p), K_ss\n",
    "    mu_post, cov_post = mu_s, cov_s + sigma_noise**2\n",
    "    \n",
    "    return dict(prior=(mu_prior, cov_prior), \n",
    "                posterior=(mu_post, cov_post))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel = squared_exponential\n",
    "params = dict(sigma=1., l=1.)\n",
    "\n",
    "X_test = np.arange(-5, 5, .05)\n",
    "p = len(X_test)\n",
    "K_ss = kernel(X_test, X_test, **params)\n",
    "mu_prior, cov_prior = np.zeros(p), K_ss\n",
    "\n",
    "N = 5\n",
    "f_priors = np.random.multivariate_normal(mu_prior, cov_prior, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig_margin=dict(top=60, bottom=40, left=50, right=0)\n",
    "\n",
    "fig = plt.figure(title='Gaussian Process Regression', \n",
    "                 layout=Layout(width='1000px', height='600px'),\n",
    "                 fig_margin=fig_margin)\n",
    "\n",
    "plt.scales(scales={'x': LinearScale(min=-5, max=5),\n",
    "                   'y': LinearScale(min=-5, max=5)})\n",
    "\n",
    "train_scat = plt.scatter([], [], colors=['salmon'], \n",
    "                         enable_move=True,\n",
    "                         interactions={'click': 'add'},\n",
    "                         marker_size=1, marker='square')\n",
    "\n",
    "prior_lines = plt.plot(X_test, f_priors, stroke_width=1, colors=['#ccc'])\n",
    "posterior_lines = plt.plot(X_test, [], stroke_width=1)\n",
    "\n",
    "mean_line = plt.plot(X_test, [], 'm')\n",
    "std_bands = plt.plot(X_test, [],\n",
    "                     fill='between',\n",
    "                     fill_colors=['orange'],\n",
    "                     fill_opacities=[.2], stroke_width=0)\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "\n",
    "# reset btn\n",
    "reset_button = Button(description='Reset Points', button_style='success')\n",
    "reset_button.layout.margin = '20px 0px 0px 70px'\n",
    "\n",
    "# controls for the plot\n",
    "f_priors_cb = Checkbox(description='Display 5 Priors?')\n",
    "f_posteriors_cb = Checkbox(description='Display 5 Posteriors?')\n",
    "std_bands_cb = Checkbox(description='Display Std Bands?')\n",
    "check_boxes = [f_priors_cb, f_posteriors_cb, std_bands_cb]\n",
    "\n",
    "# link widgets\n",
    "_ = jslink((f_priors_cb, 'value'), (prior_lines, 'visible'))\n",
    "_ = jslink((f_posteriors_cb, 'value'), (posterior_lines, 'visible'))\n",
    "_ = jslink((std_bands_cb, 'value'), (std_bands, 'visible'))\n",
    "\n",
    "def update_reg_line(change):\n",
    "    global mu_post, sig_post\n",
    "    \n",
    "    fig.animation_duration = 0\n",
    "    X_train = train_scat.x\n",
    "    y_train = train_scat.y\n",
    "    \n",
    "    fig.animation_duration = 1000\n",
    "\n",
    "    gp_res = gp_regression(X_train, y_train, X_test, sigma_noise=0.04)\n",
    "    mu_post, cov_post = gp_res['posterior']\n",
    "    \n",
    "    # simulate N samples from the posterior distribution\n",
    "    posterior_lines.y = np.random.multivariate_normal(mu_post, cov_post, N)\n",
    "    sig_post = np.sqrt(np.diag(cov_post))\n",
    "\n",
    "    # update the regression line to the mean of the posterior distribution\n",
    "    mean_line.y = mu_post\n",
    "    \n",
    "    # update the std bands to +/- 2 sigmas from the posterior mean\n",
    "    std_bands.y = [mu_post - 2 * sig_post, mu_post + 2 * sig_post]\n",
    "\n",
    "train_scat.observe(update_reg_line, names=['x', 'y'])\n",
    "\n",
    "def reset_points(*args):\n",
    "    with train_scat.hold_trait_notifications():\n",
    "        train_scat.x = []\n",
    "        train_scat.y = []\n",
    "reset_button.on_click(lambda btn: reset_points())\n",
    "\n",
    "fig.on_displayed(update_reg_line)\n",
    "VBox([fig, VBox(check_boxes + [reset_button])])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}