sampling_from_gaussian_by_hand.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Sampling from Gaussian Process by Hand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "This example is taken from [here](https://blog.dominodatalab.com/fitting-gaussian-process-models-python/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "In the function-space view of gaussian processes (see section 2.2 [here](http://www.gaussianprocess.org/gpml/chapters/RW2.pdf)) you look at gaussian processes from the perspective of sampling functions in a similar way as you normally sample random numbers from a distribution, e.g.:\n",
    "$$f(x)\\sim\\mathcal{GP}(m(x), k(x, x'))$$\n",
    "where\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "m(x)&=&\\mathbb{E}[f(x)],\\\\\n",
    "k(x,x')&=&\\mathbb{E}[(f(x)-m(x))(f(x')-m(x'))]\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "Usually the mean function $m(x)$ is taken to be zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "If you look at it from that perspective then it is interesting to draw a few samples and plot them. The purpose of this notebook is to show how that can be done."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Sampling from a Gaussian Process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np, matplotlib.pyplot as plt, seaborn as sns\n",
    "\n",
    "# https://stackoverflow.com/questions/37149933/how-to-set-max-output-width-in-numpy\n",
    "np.set_printoptions(edgeitems=10)\n",
    "np.core.arrayprint._line_width = 180\n",
    "\n",
    "sns.set()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "SEED = 42\n",
    "np.random.seed(SEED)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "For this example we take the SE (squared exponential) kernel (a.k.a. covariance function; have a look [here](http://www.gaussianprocess.org/gpml/chapters/RW4.pdf) p. 83). Have a look at the excellent Phd thesis [Automatic Model Construction with Gaussian Processes](https://github.com/duvenaud/phd-thesis/) for an overview of different choice of kernels and their properties.\n",
    "\n",
    "$$k_{SE}(r)=\\sigma^2\\cdot e^{-\\frac{r^2}{2\\ell^2}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "class GPCovarianceFunctionSquaredExponential:\n",
    "    def __init__(self, l, sigma):\n",
    "        self.l = float(l)\n",
    "        self.sigma = float(sigma)\n",
    "\n",
    "    def covariance(self, x1, x2):\n",
    "        return (self.sigma**2) * np.exp( -0.5 * (np.subtract.outer(x1, x2)**2) / (self.l ** 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "xpts = np.arange(-3, 3, step=0.01)\n",
    "k = GPCovarianceFunctionSquaredExponential(l=np.sqrt(0.1), sigma=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We will assume a zero function as the mean, so we can plot a band that represents one standard deviation from the mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigma_0 = k.covariance(0, 0)\n",
    "sigma_0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-3.0, 3.0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8YAAAIICAYAAACsBPexAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAMTQAADE0B0s6tTgAAGVpJREFUeJzt3WFo3PX9wPHP5dKy0oZNk04r0kelOnGlgxano2wWHyQd\nCn9kig8KWoQho6yRabBbReqaYmsn2NWilDm0tBooI7SWSDEocw/GqBuBlYKZUCqIJ2lTS6stSe7/\nwH/zd3Ntc7lL7q6f1wvE5JL73qfpx+I797u0UC6XywEAAABJtdR7AAAAAKgnYQwAAEBqwhgAAIDU\nhDEAAACpCWMAAABSE8YAAACkJowBAABIrbWaO09MTMT27dvj6NGj0draGu3t7bF169ZYsGBBreYD\nAACAGVXVM8YffPBBlEql6Ovri3379sW8efPizTffrNVsAAAAMOOqesZ4xYoVsWLFioiIuHjxYpRK\npfjJT35Si7kAAABgVtTkNcbbtm2L1atXx5IlS2LNmjVX/fzx8YlaPCwAAABUrVAul8u1OOjixYvR\n09MT3//+92PdunVX/NzPPjsbhUItHnXmFAoR7e1tMTJyNmrzFeJaZ2eolJ2hUnaGStkZpsPeUKlm\n2ZmOjrbLfqyqS6k//PDDGB8fj1tvvTXmzp0bnZ2dceDAgauGcUQ09Bfs68rl5pmVxmBnqJSdoVJ2\nhkrZGabD3lCpZt6Zqi6lHh4ejs2bN8fY2FhEfPXDuJYsWVKTwQAAAGA2VPWMcWdnZ/zzn/+Mhx56\nKIrFYnR0dMSWLVtqNRsAAADMuKrCuFAoxK9+9atazQIAAACzriY/lRoAAACalTAGAAAgNWEMAABA\nasIYAACA1IQxAAAAqQljAAAAUhPGAAAApCaMAQAASE0YAwAAkJowBgAAIDVhDAAAQGrCGAAAgNSE\nMQAAAKkJYwAAAFITxgAAAKQmjAEAAEhNGAMAAJCaMAYAACA1YQwAAEBqwhgAAIDUhDEAAACpCWMA\nAABSE8YAAACkJowBAABITRgDAACQmjAGAAAgNWEMAABAasIYAACA1IQxAAAAqQljAAAAUhPGAAAA\npCaMAQAASE0YAwAAkJowBgAAIDVhDAAAQGrCGAAAgNSEMQAAAKkJYwAAAFITxgAAAKQmjAEAAEhN\nGAMAAJCaMAYAACA1YQwAAEBqwhgAAIDUhDEAAACpCWMAAABSE8YAAACkJowBAABITRgDAACQmjAG\nAAAgNWEMAABAasIYAACA1IQxAAAAqQljAAAAUhPGAAAApCaMAQAASE0YAwAAkJowBgAAIDVhDAAA\nQGrCGAAAgNSEMQAAAKkJYwAAAFITxgAAAKQmjAEAAEhNGAMAAJCaMAYAACA1YQwAAEBqwhgAAIDU\nhDEAAACpCWMAAABSE8YAAACkJowBAABIrbXaA1555ZV4++23o1gsxuLFi6O3tzfmzp1bi9kAAABg\nxlX1jPHRo0fj4MGD8cYbb0RfX19cuHAh+vv7azUbAAAAzLiqnjFevnx57N+/P+bMmRMREdddd12c\nPn16SvctFKp55Jl3ab5Gn5PGYWeolJ2hUnaGStkZpsPeUKlrYWcK5XK5XIuDTpw4EWvXro19+/bF\nzTfffMXPHR+fiGLRy5sBAACov6pfYxwRcfz48Vi/fn309vZeNYojIk6dOtfw300oFCLa29tiZORs\n1OZbB1zr7AyVsjNUys5QKTvDdNgbKtUsO9PR0XbZj1UdxseOHYsNGzbE9u3bY/ny5VO+XyN/wb6u\nXG6eWWkMdoZK2RkqZWeolJ1hOuwNlWrmnakqjM+fPx/d3d2xc+fOuOWWW2o1EwAAAMyaqsL40KFD\nMTo6Gr/97W8nb7vrrrviscceq3owAAAAmA1VhfEDDzwQDzzwQK1mAQAAgFnnR0MDAACQmjAGAAAg\nNWEMAABAasIYAACA1IQxAAAAqQljAAAAUhPGAAAApCaMAQAASE0YAwAAkJowBgAAIDVhDAAAQGrC\nGAAAgNSEMQAAAKkJYwAAAFITxgAAAKQmjAEAAEhNGAMAAJCaMAYAACA1YQwAAEBqwhgAAIDUhDEA\nAACpCWMAAABSE8YAAACkJowBAABITRgDAACQmjAGAAAgNWEMAABAasIYAACA1IQxAAAAqQljAAAA\nUhPGAAAApCaMAQAASE0YAwAAkJowBgAAIDVhDAAAQGrCGAAAgNSEMQAAAKkJYwAAAFITxgAAAKQm\njAEAAEhNGAMAAJCaMAYAACA1YQwAAEBqwhgAAIDUhDEAAACpCWMAAABSE8YAAACkJowBAABITRgD\nAACQmjAGAAAgNWEMAABAasIYAACA1IQxAAAAqQljAAAAUhPGAAAApCaMAQAASE0YAwAAkJowBgAA\nIDVhDAAAQGrCGAAAgNSEMQAAAKkJYwAAAFITxgAAAKQmjAEAAEhNGAMAAJCaMAYAACA1YQwAAEBq\nwhgAAIDUhDEAAACpCWMAAABSE8YAAACkVnUYf/7557Fhw4b40Y9+VIt5AAAAYFZVHcaPP/543HHH\nHbWYBQAAAGZda7UHvPDCC3HmzJn4/e9/X9H9CoVqH3lmXZqv0eekcdgZKmVnqJSdoVJ2humwN1Tq\nWtiZqsO4ra0tzpw5U9F9rr9+fhSLjf3y5v958mCMjU9E6//NOTY+ERHxb+9f7u0rfd5MfMz5zne+\n8xvxsZ3vfOf7b9v5znd+Y50/U4/9p233RkREe3tbNKuqw3g6Tp061wTfTSj/x7//8/YrvT3bH3O+\n853fvOfX87Gd35znT/XMRp3f+bP/2F+/rRl+bc53vvOb7bFHRs5Ge3tbjIycjfJ/+yOoQXR0XD7c\n6xLGEdHQXzAAAACm5lLblcvN23kt9R4AAAAA6qmqZ4xHR0dj/fr1ceHChThz5kysXbs2li5dGps2\nbarVfAAAADCjqgrj73znO/H666/XahYAAACYdS6lBgAAIDVhDAAAQGrCGAAAgNSEMQAAAKkJYwAA\nAFITxgAAAKQmjAEAAEhNGAMAAJCaMAYAACA1YQwAAEBqwhgAAIDUhDEAAACpCWMAAABSE8YAAACk\nJowBAABITRgDAACQmjAGAAAgNWEMAABAasIYAACA1IQxAAAAqQljAAAAUhPGAAAApCaMAQAASE0Y\nAwAAkJowBgAAIDVhDAAAQGrCGAAAgNSEMQAAAKkJYwAAAFITxgAAAKQmjAEAAEhNGAMAAJCaMAYA\nACA1YQwAAEBqwhgAAIDUhDEAAACpCWMAAABSE8YAAACkJowBAABITRgDAACQmjAGAAAgNWEMAABA\nasIYAACA1IQxAAAAqQljAAAAUhPGAAAApCaMAQAASE0YAwAAkJowBgAAIDVhDAAAQGrCGAAAgNSE\nMQAAAKkJYwAAAFITxgAAAKQmjAEAAEhNGAMAAJCaMAYAACA1YQwAAEBqwhgAAIDUhDEAAACpCWMA\nAABSE8YAAACkJowBAABITRgDAACQmjAGAAAgNWEMAABAasIYAACA1IQxAAAAqQljAAAAUhPGAAAA\npCaMAQAASE0YAwAAkFprtQe8/PLLceTIkSgWi7Fs2bLYuHFjFAqFWswGAAAAM66qZ4yHhobi0KFD\nsXfv3ti/f38MDw/HkSNHajUbAAAAzLhCuVwuT/fOO3fujLGxseju7o6IiL6+vvjHP/4Rvb29V7zf\nZ5+djUZ/Unnd1sGY9hcGAAAggUJEvLpxdbS3t8XIyNmYfl3OvI6Otst+rKpLqUulUtxyyy2T7y9c\nuDA+/fTTq97v+uvnR7Ho5c0AAADNrr297d/+3Yyqfo3x1031yedTp841/DPGxWIhxsbL0Vr8atCx\n8a9+bV9//3JvX+nzZuJjzne+85v5/JaIKPuzw/kVfGxqO9O48zt/9h/7/3emGX5tzne+8xvlz46p\nfaxYLMTIyNnczxjfeOONUSqVJt//5JNP4qabbprSfRv5CwYAAMDUXGq7crl5O6+q65nvvvvueOed\nd+KLL76IsbGxOHz4cNxzzz21mg0AAABmXFXPGN92223xs5/9LNauXRstLS1x5513xo9//ONazQYA\nAAAzrurXGD/88MPx8MMP12AUAAAAmH1+NDQAAACpCWMAAABSE8YAAACkJowBAABITRgDAACQmjAG\nAAAgNWEMAABAasIYAACA1IQxAAAAqQljAAAAUhPGAAAApCaMAQAASE0YAwAAkJowBgAAIDVhDAAA\nQGrCGAAAgNSEMQAAAKkJYwAAAFITxgAAAKQmjAEAAEhNGAMAAJCaMAYAACA1YQwAAEBqwhgAAIDU\nhDEAAACpCWMAAABSE8YAAACkJowBAABITRgDAACQmjAGAAAgNWEMAABAasIYAACA1IQxAAAAqQlj\nAAAAUhPGAAAApCaMAQAASE0YAwAAkJowBgAAIDVhDAAAQGrCGAAAgNSEMQAAAKkJYwAAAFITxgAA\nAKQmjAEAAEhNGAMAAJCaMAYAACA1YQwAAEBqwhgAAIDUhDEAAACpCWMAAABSE8YAAACkJowBAABI\nTRgDAACQmjAGAAAgNWEMAABAasIYAACA1IQxAAAAqQljAAAAUhPGAAAApCaMAQAASE0YAwAAkJow\nBgAAIDVhDAAAQGrCGAAAgNSEMQAAAKkJYwAAAFITxgAAAKQmjAEAAEhNGAMAAJCaMAYAACA1YQwA\nAEBqVYfx4OBg/PCHP4z9+/fXYh4AAACYVVWF8eDgYLz11luxcuXKWs0DAAAAs6qqMF65cmXs2LEj\n5s+fX/F9C4XG/gcAAICru9RP9W64ahqvtZovQFtb27Tud/3186NYbOyXN/9p2331HgEAAKBptLdP\nrw8bwZTCeGBgIHbv3v2N2/v7+6f1oKdOnWv4Z2ULha9+Y0dGzka5XO9paAZ2hkrZGSplZ6iUnWE6\n7A2Vapad6ei4fLhPKYw7Ozujs7OzZgNFREN/wb6uXG6eWWkMdoZK2RkqZWeolJ1hOuwNlWrmnWns\n65kBAABghlX1GuNXX301BgcH46OPPoq///3vcfjw4ejp6Ynbb7+9VvMBAADAjKoqjB955JF45JFH\najULAAAAzDqXUgMAAJCaMAYAACA1YQwAAEBqwhgAAIDUhDEAAACpCWMAAABSE8YAAACkJowBAABI\nTRgDAACQmjAGAAAgNWEMAABAasIYAACA1IQxAAAAqQljAAAAUhPGAAAApCaMAQAASE0YAwAAkJow\nBgAAIDVhDAAAQGrCGAAAgNSEMQAAAKkJYwAAAFITxgAAAKQmjAEAAEhNGAMAAJCaMAYAACA1YQwA\nAEBqwhgAAIDUhDEAAACpCWMAAABSE8YAAACkJowBAABITRgDAACQmjAGAAAgNWEMAABAasIYAACA\n1IQxAAAAqQljAAAAUhPGAAAApCaMAQAASE0YAwAAkJowBgAAIDVhDAAAQGrCGAAAgNSEMQAAAKkJ\nYwAAAFITxgAAAKQmjAEAAEhNGAMAAJCaMAYAACA1YQwAAEBqwhgAAIDUhDEAAACpCWMAAABSE8YA\nAACkJowBAABITRgDAACQmjAGAAAgNWEMAABAasIYAACA1IQxAAAAqQljAAAAUhPGAAAApCaMAQAA\nSE0YAwAAkJowBgAAIDVhDAAAQGrCGAAAgNSEMQAAAKkJYwAAAFITxgAAAKQmjAEAAEittZo7nzhx\nIjZt2hQTExNx/vz56O7ujlWrVtVqNgAAAJhxVYXxM888Ew8++GD89Kc/jePHj8fPf/7zeO+992o1\nGwAAAMy4qsJ4586dMW/evIiIaG9vj9HR0SiXy1EoFK563yl8Sl1dmq/R56Rx2BkqZWeolJ2hUnaG\n6bA3VOpa2JlCuVwu1+KgzZs3R0TE008/fdXPHR+fiGLRy5sBAACovyk9YzwwMBC7d+/+xu39/f1R\nLpdjy5YtcfLkydi1a9eUHvTUqXMN/92EQiGivb0tRkbORm2+dcC1zs5QKTtDpewMlbIzTIe9oVLN\nsjMdHW2X/diUwrizszM6Ozu/cXu5XI6nnnoqCoVCvPTSSzFnzpwpD9XIX7CvK5ebZ1Yag52hUnaG\nStkZKmVnmA57Q6WaeWeqeo3xa6+9FhMTE7Ft27ZazQMAAACzqqow3rNnTyxcuDDWrl07edvzzz8f\nN9xwQ9WDAQAAwGyoKoz//Oc/12oOAAAAqAs/GhoAAIDUhDEAAACpCWMAAABSE8YAAACkJowBAABI\nTRgDAACQmjAGAAAgNWEMAABAasIYAACA1IQxAAAAqQljAAAAUhPGAAAApCaMAQAASE0YAwAAkJow\nBgAAIDVhDAAAQGrCGAAAgNSEMQAAAKkJYwAAAFITxgAAAKQmjAEAAEhNGAMAAJCaMAYAACA1YQwA\nAEBqwhgAAIDUhDEAAACpCWMAAABSE8YAAACkJowBAABITRgDAACQmjAGAAAgNWEMAABAasIYAACA\n1IQxAAAAqQljAAAAUhPGAAAApCaMAQAASE0YAwAAkJowBgAAIDVhDAAAQGrCGAAAgNSEMQAAAKkJ\nYwAAAFITxgAAAKQmjAEAAEhNGAMAAJCaMAYAACA1YQwAAEBqwhgAAIDUhDEAAACpCWMAAABSE8YA\nAACkJowBAABITRgDAACQmjAGAAAgNWEMAABAasIYAACA1IQxAAAAqQljAAAAUhPGAAAApCaMAQAA\nSE0YAwAAkJowBgAAIDVhDAAAQGrCGAAAgNSEMQAAAKkJYwAAAFITxgAAAKQmjAEAAEhNGAMAAJCa\nMAYAACC11mruPDQ0FM8991y0tLTEl19+GevWrYuurq5azQYAAAAzrqow3rNnTzzxxBOxfPnyOHny\nZNx3333R2dkZhUKhVvMBAADAjKoqjF988cXJtz/++ONYtGjRlKO40dv50nyNPieNw85QKTtDpewM\nlbIzTIe9oVLXws4UyuVyuZoDhoeHo6enJ06fPh27du2K733ve7WaDQAAAGbclMJ4YGAgdu/e/Y3b\n+/v7J98eGhqKDRs2xMGDB2P+/Pm1nRIAAABmyLSfMZ6YmIiBgYHo6uqavHz63nvvjS1btsSyZctq\nOiQAAADMlGn/dU0tLS2xa9eueP/99yMiolQqRalUisWLF9dsOAAAAJhpVb3G+Pjx4/Hss89GS0tL\nnDt3Lh599NFYs2ZNLecDAACAGVX1D98CAACAZjbtS6kBAADgWiCMAQAASE0YAwAAkFprvQdoZEND\nQ/Hcc89FS0tLfPnll7Fu3bro6uqq91g0sBMnTsSmTZtiYmIizp8/H93d3bFq1ap6j0WDGxwcjI0b\nN8Yvf/nLeOihh+o9Dg3q5ZdfjiNHjkSxWIxly5bFxo0bJ/+6RLiczz//PJ5++un429/+Fn/5y1/q\nPQ5N4JVXXom33347isViLF68OHp7e2Pu3Ln1HosGNTExEdu3b4+jR49Ga2trtLe3x9atW2PBggX1\nHq1injG+gj179sQTTzwRr7/+evzud7+LjRs3hp9VxpU888wz8eCDD8bevXujt7c3fvOb39R7JBrc\n4OBgvPXWW7Fy5cp6j0IDGxoaikOHDsXevXtj//79MTw8HEeOHKn3WDSBxx9/PO644456j0GTOHr0\naBw8eDDeeOON6OvriwsXLkR/f3+9x6KBffDBB1EqlaKvry/27dsX8+bNizfffLPeY02LML6CF198\nMZYvXx4RER9//HEsWrTId+e5op07d0ZnZ2dERLS3t8fo6KhvpnBFK1eujB07dsT8+fPrPQoN7L33\n3ovVq1fHt771rWhpaYmurq5499136z0WTeCFF15w5RJTtnz58ti/f3/MmTMnIiKuu+66OH36dJ2n\nopGtWLEiduzYERERFy9ejFKpFIsWLarzVNMjjK9ieHg47r///vj1r389+ZsOl7NgwYIoFosREbF7\n9+64//77fTOFK2pra6v3CDSBUqkUCxcunHx/4cKF8emnn9ZxIpqFP2OoRLFYnLwE9sSJE/Huu+/G\nmjVr6jwVzWDbtm2xevXqWLJkSdPujNcYR8TAwEDs3r37G7f39/fHkiVL4sCBAzE0NBS/+MUv4uDB\ng57Z4Yo7Uy6XY8uWLXHy5MnYtWtXHaajEV1pZ6BSrkQBZtLx48dj/fr10dvbGzfffHO9x6EJPPnk\nk7Fhw4bo6emJP/zhD7Fu3bp6j1QxYRwRnZ2dk5e/XjIxMRGHDx+Orq6uKBQKsWzZspg/f37861//\nimXLltVpUhrFf9uZiK/+Z/Wpp56KQqEQL7300uSlSHC5nYGpuPHGG6NUKk2+/8knn8RNN91Ux4mA\na9WxY8diw4YNsX379smXFMLlfPjhhzE+Ph633nprzJ07Nzo7O+PAgQNNGcYupb6MlpaW2LVrV7z/\n/vsR8dVlbKVSKRYvXlznyWhkr732WkxMTMTWrVtFMVAzd999d7zzzjvxxRdfxNjYWBw+fDjuueee\neo8FXGMu/Y0aO3fuFMVMyfDwcGzevDnGxsYi4qsfxrVkyZI6TzU9hbLrsS7r+PHj8eyzz0ZLS0uc\nO3cuHn300aa9Zp7ZsWrVqli4cOG/XW7//PPPxw033FDHqWhkr776agwODsZHH30UCxYsiO9+97vR\n09MTt99+e71Ho8H88Y9/jEOHDkVLS0vceeed0d3dXe+RaHCjo6Oxfv36uHDhQhw7dix+8IMfxNKl\nS2PTpk31Ho0G1dfXFzt27IilS5dO3nbXXXfFY489VsepaGTlcjl27NgRf/3rX6NYLEZHR0ds2bIl\nvv3tb9d7tIoJYwAAAFJzKTUAAACpCWMAAABSE8YAAACkJowBAABITRgDAACQmjAGAAAgNWEMAABA\nasIYAACA1P4XDj9O2pi8lJEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4cc5594e80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(15, 8), dpi= 80, facecolor='w', edgecolor='k')\n",
    "ax = plt.subplot(1, 1, 1)\n",
    "ax.errorbar(xpts, np.zeros(len(xpts)), yerr=sigma_0, capsize=0)\n",
    "ax.set_ylim(-3.0,3.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Let's select an arbitrary starting point to sample, say $x=1$. Since there are no prevous points, we can sample from an unconditional Gaussian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.4967141530112327]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = [1.]\n",
    "y = [np.random.normal(scale=sigma_0)]\n",
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigma_1 = k.covariance(x, x)\n",
    "sigma_1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "In order to create new function points `y_new` at locations `x_new` after we created a set of initial ones at locations `x` with their corresponding `y` values we have to use the conditional probability as described [here](http://www.gaussianprocess.org/gpml/chapters/RWA.pdf) in equation A.6.\n",
    "\n",
    "Even while I've taken the example from [here](https://blog.dominodatalab.com/fitting-gaussian-process-models-python/) I will change the notation to fit to the [Gaussian Processes for Machine Learning](http://www.gaussianprocess.org/gpml/chapters/) notation. I will use $y_*$ instead of `y_new`, though.\n",
    "\n",
    "Given a joint distribution:\n",
    "$$\\left[\\begin{matrix}y_*\\\\ y\\end{matrix}\\right]\\sim\\mathcal{N}\\left(\n",
    "\\left[\\begin{matrix}\\bar y_*\\\\\\bar y\\end{matrix}\\right],\n",
    "\\left[\\begin{matrix}A&C\\\\ C^T&B\\end{matrix}\\right]\n",
    "\\right)$$\n",
    "\n",
    "We get the conditional distribution analytically as follows:\n",
    "\n",
    "$$p(y_*\\,|\\, X, y, X_*)\\sim\\mathcal{N}(\\bar y_*+CB^{-1}(y-\\bar y), A-CB^{-1}C^T)$$\n",
    "\n",
    "And as we set the mean function $m(x)=0$ this reduces to:\n",
    "\n",
    "$$p(y_*\\,|\\, X, y, X_*)\\sim\\mathcal{N}(CB^{-1}y, A-CB^{-1}C^T)$$\n",
    "\n",
    "Please also keep in mind that both, $y_*$ and $y$ can be vectors and do not need to be scalars."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "In case that $A$, $B$ and $C$ are defined by the covariance function and we only look at predicting one new point $(x_*,y_*)$ we have the following:\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "A&=&k(x_*, x_*)\\space\\hbox{is a scalar}\\\\\n",
    "B&=&k(x_*, x)=k_*\\space\\hbox{is a vector}\\\\\n",
    "C&=&k(X,X)\\space\\hbox{is a matrix}\\\\\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "See also [here](http://www.gaussianprocess.org/gpml/chapters/RW2.pdf) equations 2.25 and 2.26 for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def conditional(x_new, x, y, cov, sigma_n=0):\n",
    "    total_covariance_function = cov\n",
    "\n",
    "    A = total_covariance_function.covariance(x_new, x_new)\n",
    "    C = total_covariance_function.covariance(x_new, x)\n",
    "    B = total_covariance_function.covariance(x, x) + np.power(sigma_n,2)*np.diag(np.ones(len(x)))\n",
    "\n",
    "    mu = [np.linalg.inv(B).dot(C.T).T.dot(y).squeeze()]\n",
    "    sigma = [(A - C.dot(np.linalg.inv(B).dot(C.T))).squeeze()]\n",
    "\n",
    "    return (mu, sigma)\n",
    "\n",
    "def predict(x_new, x, y, cov, sigma_n=0):\n",
    "    l_y_pred, l_sigmas = conditional(x_new, x, y, cov=cov, sigma_n=sigma_n)\n",
    "    if len(l_sigmas[0].shape) > 1:\n",
    "        return l_y_pred, [np.diagonal(ls) for ls in l_sigmas]\n",
    "    else:\n",
    "        return l_y_pred, l_sigmas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Now we can draw the predicted function value together with the sigmas via the predict method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "y_pred, sigmas         = predict(xpts, x, y, cov=k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-3.0, 3.0)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8YAAAIICAYAAACsBPexAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAMTQAADE0B0s6tTgAAIABJREFUeJzt3W+MXFeZJ+D3VrWtRI5nCG0PhM1E+8GTAAKvGcVEQYqY\nNmjUzixIIwQoWmUVDAGhUUScpePBjCMUxs5iE5AIjkWTJQxkEjCKRl6byMiLkyzhA0SBlaWJvIsH\nKUp2Ec04ceKJ86/r3v1QrnZ1ubvd1VXVdavO80iJq+vPqePqU9f1q/ecc7OiKIoAAACARFX63QEA\nAADoJ8EYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSRjp5cJ7nsWfPnnjqqadi\nZGQkRkdH46677opLLrmkW/0DAACAnuqoYvyrX/0qpqamYv/+/fHggw/GxRdfHD/84Q+71TcAAADo\nuY4qxldffXVcffXVERHx+uuvx9TUVPzFX/xFN/oFAAAAy6Ira4x3794dmzZtinXr1sX1119/wfvX\nank3nhYAAAA6lhVFUXSjoddffz22bdsW7373u2PLli0L3vcPfzgdWdaNZ+2dLIsYHV0dJ0+eju68\nQgw7Y4Z2GTO0y5ihXcYMS2Hc0K5BGTNr1qye97aOplL/5je/iVqtFm9/+9tj5cqVMT4+Hg8//PAF\ng3FElPoFa1YUg9NXysGYoV3GDO0yZmiXMcNSGDe0a5DHTEdTqU+cOBF33nlnTE9PR0R9M65169Z1\npWMAAACwHDqqGI+Pj8c///M/xw033BDVajXWrFkTO3fu7FbfAAAAoOc6CsZZlsXnP//5bvUFAAAA\nll1XdqUGAACAQSUYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRj\nAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYY\nAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTB\nGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJ\nxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRN\nMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBp\ngjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJ\nE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABI\nmmAMAABA0kY6bWBycjJ+8pOfRLVajSuuuCJ27doVK1eu7EbfAAAAoOc6qhg/9dRTcfDgwfjBD34Q\n+/fvj9deey0OHDjQrb4BAABAz3VUMd6wYUM89NBDsWLFioiIuPTSS+OFF15Y1GOzrJNn7r1G/8re\nT8rDmKFdxgztMmZolzHDUhg3tGsYxkxWFEXRjYaeeeaZuPHGG+PBBx+Myy+/fMH71mp5VKuWNwMA\nANB/Ha8xjog4fvx43HLLLbFr164LhuKIiOeff7n03yZkWcTo6Oo4efJ0dOerA4adMUO7jBnaZczQ\nLmOGpTBuaNegjJk1a1bPe1vHwfjpp5+OW2+9Nfbs2RMbNmxY9OPK/II1K4rB6SvlYMzQLmOGdhkz\ntMuYYSmMG9o1yGOmo2B85syZ2Lp1a9xzzz1x1VVXdatPAAAAsGw6CsaHDh2KU6dOxd///d/PXPe+\n970vPvvZz3bcMQAAAFgOHQXjj33sY/Gxj32sW30BAACAZWdraAAAAJImGAMAAJA0wRgAAICkCcYA\nAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAG\nAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIx\nAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROM\nAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpg\nDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIE\nYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJIm\nGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0\nwRgAAICkCcYAAAAkTTAGAAAgaSP97gAAwFw+vefRmK4VERExUs0iImK6Vsy6PN9tERGTE2PL2V0G\nwKf3PDpzeSljy5iC4SUYAwDLajGBtzngLsV0rYgt//WooJy4xljr1tha7Ng11mDwdByMX3rppbjj\njjviySefjJ///Ofd6BMAMKAWW5FbbtO1YqZvwsvwmiu4lqUfZjRAuXUcjG+77bb4wAc+EE8++WQ3\n+gMAlMxCVbcyhd92Nf+9hJTB1To+B0HrjIaF3lMRQjQsh46D8de//vV48cUX45vf/GZbj8tKfuxq\n9K/s/aQ8jBnaZczQrqWMmZt3Ly3UlqHqtlwa1eRv3z584WOYjzM37370wncaAnPNdmhcjlj8+7ed\n8T3M44beGIYx03EwXr16dbz44ottPebNb14V1Wq5N8T+69sPxnQtj5Gz/Zyu5RERs36e7/JC9+vF\nbdrXvva1X8bn1n5Z2q9/QI5ofFopWi4v9bZutFGu9m/e/djMa/dPuz8Uw2R0dHW/u9A1jc9oEefG\nf9nHVr/bn67lcfPux2YuR5Th2KT9frTfq+duHDMH+VjTl823nn/+5QH4NqFo+bP1+oUuL/dt2te+\n9ge3/X4+t/YHs/3FtlnW/g9G+399+3+PiPaqbGWUZfUPqidPno6i9a88YOauEA/e2NK+9vvbfm+e\n++TJ0wNxrFmzZv7g3rddqcv8ggEATNeKuHn3o0OxvrMoBvuz1yCuI4aUNI4vg3ysqVz4LgAAaWpe\n38ny+/SeR73+wLLoqGJ86tSpuOWWW+K1116LF198MW688ca48sorY8eOHd3qHwBA3zXC2TBUjweF\nKjGwnDoKxm9605vi+9//frf6AgBQWo3qsXDce83nAQZYDn1bYwwAMIhUj3vHtGmgX6wxBgBok7XH\n3adKDPSTYAwAAEDSBGMAgCVQNe4OO08DZSAYAwB0QLBbOtOngbIQjAEAOqR63B5fJgBlIxgDAACQ\nNMEYAKBLVEIvzPRpoIwEYwCALjKten5CMVBWI/3uAADAMGqE48mJsT73pP98UQCUnYoxAAAASROM\nAQB6xLRq06eBwSAYAwD0UKrh2EZkwCARjAEAAEiaYAwAsAxSqqCaPg0MGsEYAGCZpDCtWigGBpFg\nDAAAQNKcxxgAYJkN4zmOh70SDgw3FWMAgD4YpmnVpk8Dg04wBgAAIGmCMQBAnwx61TilnbaB4SYY\nAwD02SAGTNOngWEiGAMAAJA0wRgAoAQGaVq1ajEwbJyuCQCgRMp8KqdBCe4A7VIxBgAomTJWj1WJ\ngWEmGAMAAJA0wRgAoITKUjUexB2zAdolGAMAlFg/g6np00AqBGMAgJJb7uqxKjGQGsEYAGBALEdg\nVSUGUiQYAwAMkHaqxyO//EXEP/xD/c9FEIqBVDmPMQDAgGkNx63nPK4892z80U3/KVYc+18REfGm\niHhj/YZ46bv/GPnlfzrrvo0wPFLNet5vgLISjAEABlgjJN972/vj9Jk34vSZ16O4bUfk/7Yi3viz\na2O6MhIj+XSsPP16jGz9YmT33hurV62M//LNJyLLhGGACMEYAGAg5XkRjUnP9XD82Lkb//w/R/z5\nPA/8b79s+qGIRjTOi3OXAVIjGAMAlFhR1ONvnhdRNC3/zedZCvymM6fiqv/3v+OiN16NFbU3oprn\n8UZ1JKarI3Fm5cXxfy67Kl66+I/Otd9oL5/dTp4XoaAMpEIwBgAoqVpTGG4NwlkWkZ29vlqprzPO\nsixGfvmLuPQ/3jRvmy8cOhLT770mbt59NCLqwTjP6201P0VenLuilhdREZKBIWZXagCAksiLImpN\nCbi5QpxlEZWmT27VShaVShYj1Sy+ffummfXC0++9Jt5Yv2HO9t/4D++J6fdeExER3759U3z79k1R\nOfu4ajWbtQFXc7W4KCJqZyvK9bBu52pguAjGAAB9VBRF5HljuvTsMFzJ6tXgiLNBuGVu8+TE2Hk7\nUkdEvPTdfzwvHL+xfkO8dP8D5913cmJszh2pq5Vs5rkr84TkRv8BBp2p1AAAy6w5TNZa1vZWsnPT\npisziXR2+BypZnMG4ob88j+NU//jf8aKJ38Rb/rD/41Ta/9dvLHxmnnv32ir9fzI9Sp0EZVKFllR\nRC2vV5Kbs3D9OuEYGGyCMQDAMsrzYs71wsXZtcJZlkVe607QnH7vNRFrVsf0v55uzdZta4TkaiWL\n4mxIbmgE5VqtmDXdG2BQOHQBAPTYrOnSLQG1WqlPW46IRZ1X+ELV4k7MN626VXM/62G+frkI06yB\nwSQYAwD0WC2fHYib1+wuJgw3zLemuJsWG44bsiybCfatO1fX/97CMVB+gjEAQJcVRTFvIKxWmtcO\nD5dK04ZdDY3zI+eF3ayB8rLGGACgSxrB77wNtSr18wTXN6pqPxT3cvr0XObbjGsxGmuRI2ZvJJbn\n8z8GoN9UjAEAuqB1Q6pmlSxbUiAedHNVkBtUj4EyUTEGAOjAXFXiSlbfkGq+oNyO5awUz/XcS6ka\nNzu3m3V99+pGBbmWR2SdbpUN0CUqxgAAS1TLi/MCcUS9UppihXghWZadt7a6EYtrufXHQH8JxgAA\nbWgOcM1Zrtubai33uuL59Gon7FmneSrOVdcFZKAfBGMAgEVqrRBnWcysoVUhbk/zaZ6aX7paHjPn\nfAZYLoIxAMACiqKYCWqtxcxqj6ZMl6Va3Kzd8xu3o9qySVdzLlZBBpaDYAwAsIBafi6oZRHz7rJM\nZxpfMFTPntqqof76C8dAbzm0AwA0aa4QNzSCWrVqU61eq2/SNfu6xjmQBWSgV5yuCQAg5j7tUkRE\npVI/D/F0bXlCWdmmUDdr9K3TUzhdSOMUTxH1nb4b31PkeUThFE9AD6gYAwDJy4ti3nMOV1SI+6rS\nsv64ORZbfwx0i2AMACSpOVTlLeci7sc64jJuuDWfXm7ENZfm9ceVlh2sa3awBrpAMAYAklPMUSFu\nBK5Kj3aapnP19cezfzeN7zdqeaGCDCyZYAwAJKEoipnNm1pDcbUS5wUuyq9aOXcO5KI4//cKsFg2\n3wIAhtq8m2o1berU7wrxoEyhbjY5MdbzTbguJMuyqGYR07Uismz2eaZreRG+6wAWS8UYABhaeT73\nploqxMOn2rpJV1MFOS9MswYWJhgDAENl1qZaTVmoeVOtfleIGwZpw625TE6Mlar/zb/X5l9xnp8L\nyQIyMBdTqQGAoVHfgGn2dZVKPRidqxALRimoVrI5N1mr5RFZZgwAswnGAMBAy/NipjLcGoobpxTK\nheEk1SvIRVQr9a9D8pmq8bn75EUR5Zg/APSTqdQAwEApWtaLzjdduuwGfRp1s+U+r3G7siyLStO5\nkE2zBlqpGAMAA2G+3aWziMhMl2aRmneyblUfW8YOpEgwBgBKqbl6N1eIaaiaLk0HGjMM5tq9vHb2\nNFDA8BOMAYDSOFcVPn8TrYaFgsygGJYp1M0af6d+n9u4Xed2sp5jLXKcW488XXNeZBhmgjEA0DfN\nVeFarZip+baG4uYw3BxkoJuyLIsszs0+qFTqY7ExHpvXs9dym3bBMBmQ7SkAgGHQunFWrXnjo6b7\ntW6ilWVZac493Ilh2nBrPmXfiKsdlSyLamXuTbuK4lxQnq4VUTv7Q+sYBwaDijEA0FN5fq4SvND0\n50qlvpFWzSZalFDrpl2VbPZU68afrWPc6aBgMAjGAEDH6lWy+uXW9cH5Atm2XoXLzq7fFIYZHI0v\nb6ZrZ9cmF3OP9bwpKE+3bOalsgzlIRgDAPOafb7gYiaztobf5ipZ62f9LKtXgvNidhCu35ZOLW3Y\np1A3m5wYG7hNuDpRn+ofkdeKGKlmURTFzHsii9lf9cz3vqnVimguLTdXmouiSOq9Av3QcTD+1re+\nFUeOHIlqtRrr16+P7du3e+MCQMk0B9yiKGZ9UJ811bk2+7bmD+75AuE3oh6Ai2L2lOiImFmjmdd8\nuCcN9XFef5M0Tid2ocpyMfO/uub3W/P5lc8L0PnsL6+a32Eq0rB4HQXjY8eOxaFDh+JHP/pRrFy5\nMj75yU/GkSNH4i//8i+71b++aRxI8pYj16yDzzyXl/s27Wtf+4Pe/jD/3bS/2PZrTT/XFrxt9gLG\n5iDbfK7f1vP+1s77kN3Ux6a7zn7UbI1MO1f4HWn68G9KNMytubIcEbOqy62hufFFU6ti5n91ze/f\nlsPDrPd66zGh+eda6+WmdH2hY1PjrmU9tmq/9889LF/AdBSMH3/88di0aVNcdNFFERGxefPmeOyx\nxxYVjMv8hfG/nXlj5kDSMgZmH3zmubzct2lf+9of5PaLRd2vN8+t/TK13zwUWj9jLHhbtC87+79G\nW419rvJi4cBbrQi/SzFSzeLbt/dnGnXj19SPz12Nv/PNux9d/icfII3qcmtorrasYY6IcwE6YiYE\nN28CttR3ZdF6eZ5jzlzHpsZVZT22ar/3z13LI06+9EqsWbO61BnvQjoKxlNTU3HVVVfN/Lx27dr4\n/e9/f8HHvfnNq6JaLe+ZokaL+j/4eVHM+nY8Yva35fNdXuh+vbhN+9rXvvbL+NzL3X4jqOVFMevy\nUm/rRhtdbT+LyPNiZsOfxrf1jQ/PtbyY93JExMjZf3ena/nMv8HTZ9NvpVL/Oa/lUcka/z43Sk1n\nX+goFri80P2WetswtZ/FmjWro59GR/v5/FmU93czGO1nTe/LLKucPdfy7PfvXO/tWe/7lmPCfMeL\n1mPMUo9bZfm3Qfu9f+5KlsWf/fu1EdHvY01nurr51mLL6M8///LMi1lWlUpEXpv9bV1E67d3c19e\n6H69uE372tf+ILdfiYjCsaML7Z/9bBh5bfblpd7WjTa63n40PogW52Jr07+n812uK+a5vNTbutFG\nKu0X8a//ejr6IcvqH1RPnjwdi/yY1gNFy5+t1y90uRe3pdl+6zFhvuPFzAySxs9LPG6V5d8G7ff+\nuSuViBde+LcSHGsubKEvKTsq2771rW+NqampmZ9/97vfxdve9rZFPbYoyv0fAECnRqpZTE6M9f0z\nTT+ff3JibGZKPjCcynCs6TTjdRSMx8bG4qc//Wm88sorMT09HY888kh88IMf7KRJAAAAWFYdTaV+\n5zvfGR/96EfjxhtvjEqlEtdee228//3v71bfAAAAoOc6XmN80003xU033dSFrgAADI/JibF+d6E0\nGq/Fp/c82ueeAMytvFtDAwAAwDIQjAEAAEiaYAwA0EWNnag5nx2qgbISjAEAAEiaYAwAAEDSOt6V\nGgCAOlOoL2xyYszu1EDpqBgDAACQNMEYAACApAnGAAAdshN1eyYnxrxeQKkIxgAAACTN5lsAAB1Q\n+Vw6G3EBZaFiDAAAQNIEYwAAAJImGAMALIENt7rDRlxAGQjGAAAAJE0wBgAAIGmCMQBAm0yj7r7J\nibEYqWb97gaQKMEYAACApAnGAAAAJG2k3x0AABgkplD3TuO1/fSeR/vcEyA1KsYAAAAkTTAGAFgE\nG24tHxtxActNMAYAACBpgjEAAABJs/kWAMAFmEK9/CYnxmzCBSwbFWMAAACSJhgDAMzDhlv9NTkx\n5vUHloVgDAAAQNIEYwAAAJImGAMAzME06vJwXmOg1wRjAAAAkuZ0TQAALVSKy6fxO3EKJ6AXVIwB\nAABImmAMAHCWdcXlZ70x0AuCMQAAAEkTjAEAAEiazbcAAMKGW4NkcmLMJlxAV6kYAwAAkDTBGABI\nmg23BtPkxJjfG9A1gjEAAABJE4wBgGSpFg8+p28CukEwBgAAIGl2pQYAkqRSPDwav0s7VQNLpWIM\nAABA0gRjACAp1hUPL+uNgaUSjAEAAEiaNcYAQDJUioff5MSYtcZA21SMAYChZ/p0WiYnxvy+gbYI\nxgAAACTNVGoAYKipHKbLtGpgsVSMAQAASJpgDAAMJeuKibDeGFgcwRgAGDpCMa2c4xhYiDXGAMBQ\nEYiZT2NsWHcMtFIxBgAAIGmCMQAwFEyfZrFMqwZamUoNAAw8gZh2OZUT0EzFGAAYWKrEdMKO1UCD\nYAwADCShmG4xtRowlRoAGDgCMd1mx2pIm4oxADAwVInpNdVjSJOKMQAwEARilouNuSA9KsYAQKmp\nEtMPNuaCtAjGAEBpCcX0m6nVkAZTqQGAUhKIKQsbc8HwE4wBgNIQhikzARmGl6nUAEDfmTLNILH+\nGIaPijEA0BfCMIOuefyqIsNgE4wBgGUlDDOMhGQYbIIxANAzqsKkqHXMC8pQfoIxANA1gjCcb66g\nPF0r+tQbYC4dB+OjR4/G9u3b43Of+1zccMMN3egTADAABGBYmvneOyrL0D8dBeOjR4/Gj3/849i4\ncWO3+gMAlITqLyyvhd5vQjP0VkfBeOPGjbFp06b427/927Yfm2WdPDMAsBgj1Sy+fbtwm6rG5y2f\nuwZfJ+/jm3cL1fTWMBxrOgrGq1evXtLj3vzmVVGtlvsUyv+0+8P97gIAQFeMji7tMxvDwedalssg\nH2sWFYwPHz4c+/btO+/6AwcOLOlJn3/+5dJ/m5Bl9V/syZOno7A3AotgzNAuY4Z2GTO0y5hhKYwb\n2jUoY2bNmvmD+6KC8fj4eIyPj3etQxFR6hesWVEMTl8pB2OGdhkztMuYoV3GDEth3NCuQR4z5Z7P\nDAAAAD3W0Rrj+++/P44ePRq//e1v49e//nU88sgjsW3btnjXu97Vrf4BAABAT3UUjD/xiU/EJz7x\niW71BQAAAJadqdQAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgA\nAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYA\nAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAG\nAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIx\nAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROM\nAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpg\nDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIE\nYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJIm\nGAMAAJA0wRgAAICkjXTy4GeeeSZ27NgReZ7HmTNnYuvWrXHdddd1q28AAADQcx0F4y996Uvx8Y9/\nPP7qr/4qjh8/Hp/5zGfi8ccf71bfAAAAoOc6Csb33HNPXHzxxRERMTo6GqdOnYqiKCLLsgs+dhF3\n6atG/8reT8rDmKFdxgztMmZolzHDUhg3tGsYxkxWFEXRjYbuvPPOiIi44447LnjfWi2PatXyZgAA\nAPpvURXjw4cPx759+867/sCBA1EURezcuTOeffbZ2Lt376Ke9PnnXy79twlZFjE6ujpOnjwd3fnq\ngGFnzNAuY4Z2GTO0y5hhKYwb2jUoY2bNmtXz3raoYDw+Ph7j4+PnXV8URXzhC1+ILMvi3nvvjRUr\nViy6U2V+wZoVxeD0lXIwZmiXMUO7jBnaZcywFMYN7RrkMdPRGuPvfe97ked57N69u1v9AQAAgGXV\nUTC+7777Yu3atXHjjTfOXPfVr3413vKWt3TcMQAAAFgOHQXjn/3sZ93qBwAAAPSFraEBAABImmAM\nAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRj\nAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYY\nAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTB\nGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJ\nxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRN\nMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBp\ngjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJ\nE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpI508+NixY/GVr3wl\nKpVKvPrqq7Fly5bYvHlzt/oGAAAAPddRML7vvvtiYmIiNmzYEM8++2x8+MMfjvHx8ciyrFv9AwAA\ngJ7qKBj06JXGAAAF40lEQVR/4xvfmLn83HPPxWWXXbboUFz27NzoX9n7SXkYM7TLmKFdxgztMmZY\nCuOGdg3DmMmKoig6aeDEiROxbdu2eOGFF2Lv3r3xjne8o1t9AwAAgJ5bVDA+fPhw7Nu377zrDxw4\nMHP52LFjceutt8bBgwdj1apV3e0lAAAA9MiSK8Z5nsfhw4dj8+bNM9OnP/ShD8XOnTtj/fr1Xe0k\nAAAA9MqST9dUqVRi79698cQTT0RExNTUVExNTcUVV1zRtc4BAABAr3W0xvj48ePx5S9/OSqVSrz8\n8svxqU99Kq6//vpu9g8AAAB6quPNtwAAAGCQLXkqNQAAAAwDwRgAAICkCcYAAAAkbaTfHSizY8eO\nxVe+8pWoVCrx6quvxpYtW2Lz5s397hYl9swzz8SOHTsiz/M4c+ZMbN26Na677rp+d4uSO3r0aGzf\nvj0+97nPxQ033NDv7lBS3/rWt+LIkSNRrVZj/fr1sX379pnTJcJ8XnrppbjjjjviySefjJ///Of9\n7g4DYHJyMn7yk59EtVqNK664Inbt2hUrV67sd7coqTzPY8+ePfHUU0/FyMhIjI6Oxl133RWXXHJJ\nv7vWNhXjBdx3330xMTER3//+9+NrX/tabN++PexVxkK+9KUvxcc//vF44IEHYteuXfF3f/d3/e4S\nJXf06NH48Y9/HBs3bux3VyixY8eOxaFDh+KBBx6Ihx56KE6cOBFHjhzpd7cYALfddltcc801/e4G\nA+Kpp56KgwcPxg9+8IPYv39/vPbaa3HgwIF+d4sS+9WvfhVTU1Oxf//+ePDBB+Piiy+OH/7wh/3u\n1pIIxgv4xje+ERs2bIiIiOeeey4uu+wy386zoHvuuSfGx8cjImJ0dDROnTrlyxQWtHHjxrj77rtj\n1apV/e4KJfb444/Hpk2b4qKLLopKpRKbN2+Oxx57rN/dYgB8/etfN3OJRduwYUM89NBDsWLFioiI\nuPTSS+OFF17oc68os6uvvjruvvvuiIh4/fXXY2pqKi677LI+92ppBOMLOHHiRHzkIx+JL37xizO/\ndJjPJZdcEtVqNSIi9u3bFx/5yEd8mcKCVq9e3e8uMACmpqZi7dq1Mz+vXbs2fv/73/exRwwKxxja\nUa1WZ6bAPvPMM/HYY4/F9ddf3+deMQh2794dmzZtinXr1g3smLHGOCIOHz4c+/btO+/6AwcOxLp1\n6+Lhhx+OY8eOxd/8zd/EwYMHVXZYcMwURRE7d+6MZ599Nvbu3duH3lFGC40ZaJeZKEAvHT9+PG65\n5ZbYtWtXXH755f3uDgPg9ttvj1tvvTW2bdsW3/nOd2LLli397lLbBOOIGB8fn5n+2pDneTzyyCOx\nefPmyLIs1q9fH6tWrYp/+Zd/ifXr1/epp5TFXGMmov5h9Qtf+EJkWRb33nvvzFQkmG/MwGK89a1v\njampqZmff/e738Xb3va2PvYIGFZPP/103HrrrbFnz56ZJYUwn9/85jdRq9Xi7W9/e6xcuTLGx8fj\n4YcfHshgbCr1PCqVSuzduzeeeOKJiKhPY5uamoorrriizz2jzL73ve9Fnudx1113CcVA14yNjcVP\nf/rTeOWVV2J6ejoeeeSR+OAHP9jvbgFDpnFGjXvuuUcoZlFOnDgRd955Z0xPT0dEfTOudevW9blX\nS5MV5mPN6/jx4/HlL385KpVKvPzyy/GpT31qYOfMszyuu+66WLt27azp9l/96lfjLW95Sx97RZnd\nf//9cfTo0fjtb38bl1xySfzJn/xJbNu2Ld71rnf1u2uUzHe/+904dOhQVCqVuPbaa2Pr1q397hIl\nd+rUqbjlllvitddei6effjre8573xJVXXhk7duzod9coqf3798fdd98dV1555cx173vf++Kzn/1s\nH3tFmRVFEXfffXf84he/iGq1GmvWrImdO3fGH//xH/e7a20TjAEAAEiaqdQAAAAkTTAGAAAgaYIx\nAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICk/X9/Z6jLiTypGwAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4cbd4974e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(15, 8), dpi= 80, facecolor='w', edgecolor='k')\n",
    "ax = plt.subplot(1, 1, 1)\n",
    "ax.errorbar(xpts, y_pred[0], yerr=sigmas[0], capsize=0)\n",
    "ax.plot(x, y, \"ro\")\n",
    "ax.set_ylim(-3.0,3.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Now let's take a second arbitrary location on the $x$-axis, e.g. $-0.7$ and let's repeat the process from before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[array(2.633608839077876e-07)]\n",
      "[array(0.9999999999997189)]\n"
     ]
    }
   ],
   "source": [
    "new_dp = [-0.7]\n",
    "m, s = predict(new_dp, x, y, cov=k)\n",
    "print(m)\n",
    "print(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Drawing a random number with the mean function and standard deviation (sigma):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[array(2.633608839077876e-07)]\n",
      "[array(0.9999999999997189)]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-0.1382640378102619"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y2 = np.random.normal(m[0], s[0])\n",
    "print(m)\n",
    "print(s)\n",
    "y2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Now we add the new x-point to the array of x-values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.0, -0.7]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# x = [1.0]\n",
    "x += new_dp\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "and the new y-point to the array of y-values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.4967141530112327, -0.1382640378102619]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y.append(y2)\n",
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "And repeat the process from above to predict how the function should look like along the x-axis given the two points we just fixed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "y_pred, sigmas         = predict(xpts, x, y, cov=k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-3.0, 3.0)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8YAAAIICAYAAACsBPexAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAMTQAADE0B0s6tTgAAIABJREFUeJzt3X+spHV9L/DPM3PYSGCtuLtVvEhubrZgG11XLz+iCdXd\n0maXVhNjlZCGBqloTMMVuN1dxUIMulAXqYkUiAu3WKWgGGK2IMFQV2j1DyRYu0kJ97o1IXCv8VRg\ngfL7zDz3jzlzzjNzzpydOTNz5nnm+3olysyZmed895zv85zn/Xy+3++T5XmeBwAAACSqNukGAAAA\nwCQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkjYzzIebzWZce+218cgjj8TM\nzExs2LAhrrnmmjj++ONH1T4AAAAYq6Eqxj/96U9jdnY27rzzzrj99tvj2GOPjW9/+9ujahsAAACM\n3VAV49NOOy1OO+20iIh49dVXY3Z2Nt7//vePol0AAACwJkYyx3jfvn2xffv22Lx5c5xzzjlHfX+j\n0RzFtwUAAIChZXme56PY0Kuvvhp79uyJd7zjHXHhhReu+N7/+I/nI8tG8V3HJ8siNmxYH0899XyM\n5ifEtNNnGJQ+w6D0GQalz7Aa+g2Dqkqf2bhxfc/XhhpK/fOf/zwajUa87W1vi3Xr1sWOHTvirrvu\nOmowjohS/8CK8rw6baUc9BkGpc8wKH2GQekzrIZ+w6Cq3GeGGkp9+PDhuOqqq2Jubi4iWotxbd68\neSQNAwAAgLUwVMV4x44d8W//9m9x3nnnRb1ej40bN8bevXtH1TYAAAAYu6GCcZZl8Rd/8RejagsA\nAACsuZGsSg0AAABVJRgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDS\nBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACS\nJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQ\nNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACA\npAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAA\nJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAA\nIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAA\nAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEA\nAEiaYAwAAEDSZobdwP79++P73/9+1Ov1OPnkk+Pqq6+OdevWjaJtAAAAMHZDVYwfeeSRuPvuu+Nb\n3/pW3HnnnfHKK6/EgQMHRtU2AAAAGLuhKsZbt26NO+64I4455piIiDjhhBPimWee6euzWTbMdx6/\ndvvK3k7KQ59hUPoMg9JnGJQ+w2roNwxqGvpMlud5PooNPf7443H++efH7bffHieddNKK7200mlGv\nm94MAADA5A09xzgi4rHHHouLL744rr766qOG4oiIp59+ofRXE7IsYsOG9fHUU8/HaC4dMO30GQal\nzzAofYZB6TOshn7DoKrSZzZuXN/ztaGD8aOPPhqXXHJJXHvttbF169a+P1fmH1hRnlenrZSDPsOg\n9BkGpc8wKH2G1dBvGFSV+8xQwfjFF1+MSy+9NK6//vo49dRTR9UmAAAAWDNDBeN77rknjhw5El/8\n4hcXvvbe9743PvWpTw3dMAAAAFgLQwXjj370o/HRj350VG0BAACANWdpaAAAAJImGAMAAJA0wRgA\nAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYA\nAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAG\nAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIx\nAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROM\nAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpg\nDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIE\nYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJIm\nGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaTOTbgAAwHI+ce0PY66RR0TETD2LiIi5Rt7xuNdrERH7\nd21by+ZSAZ+49ocLj1fTt/QpmF6CMQCwpvoJvMWAuxpzjTwu/KuDgnLi2n1tVH2r376rr0H1DB2M\nn3vuubjyyivj4Ycfjh//+MejaBMAUFH9VuTW2lwjX2ib8DK9lguuZWmHEQ1QbkMH48suuyx+7/d+\nLx5++OFRtAcAKLHlKnAR5Qi/gyr+W4SU6uruk1XQPaJhpX0qQoiGtTB0MP7KV74Szz77bPzN3/zN\nQJ/LSn7sarev7O2kPPQZBqXPMKhx9JmL9nVWeKch8A6qXU2+eff0hY9pPs4U++40W260Q/txRH9z\npCNioP49zf2G8ZiGPjN0MF6/fn08++yzA33mjW88Lur1ci+I/aHdd8dcoxkz8+2cazQjIjqe93q8\n0vvG8Zrt236Zth8R8d19Hwj6t2HD+kk3gYpZqc8M+ver+FpEHhFZ4XF0Pe/1eKX3rfa1td/+Rfse\nWPj5TNtxbJqOM+0+HqHv9ruNuUYzLtr3wMLjiOk577D9cnzv9jGzyseaiSy+9fTTL1TgakLe9d/u\nr6/0eK1fs33bL8/25xp5fGj3Pyw8Xosr21WVZa0/IE899Xzk3T9OWEarQpZFRN5HRbdax46ybb99\nHKv6sWiajjPLV4ir17ds3/Ynu/3xfO+nnnq+EseajRt7B/eJrUpd5h8YMFlzjXzhBCiFW2bkuWMi\ni8q6eFWK2seiaTjeVP04U8V5xJCS9vGlyscat2sCKqnXLTim4QSWtIz6djKMVnt+p2PLZBQvFAGM\n01DB+MiRI3HxxRfHK6+8Es8++2ycf/75ccopp8QVV1wxqvYB9K24QEmEVTwpp16rOlNu7WOL48ra\nUSUG1tJQwfgNb3hDfPOb3xxVWwBGxj1LKRMn+NNB9XjtFO8DDLAWDKUGkqCazFoThqeX6vH4GDoN\nTEq575kEMAbdIRlGTbVr+jmOjJ79BpgkFWMgWcWKnsoPw1IhBoDqUjEGkqfywzA+ce0P9Z9EOXaM\nhn0IKAPBGGCekzMGZegnEY4dw7APAWUhGAMUqADRLyf0FDl2DMbFBKBszDEGWIZVZ+nFyTwATB8V\nY4AeVIDopkrM0aiEHp39CCgjwRjgKJzoog8wCBfVehOKgbISjAH64EQ3XU7kWS0XVBb5WQBlJxgD\n9Ek4To9QDABpEIwBBqTyMf38jhkVF9RcYAKqQTAGWAUnu9PLSTyjlurxwgUmoEoEYwAAAJLmPsYA\nQ3C/4+mhssW4pXS8aI+8mKlnk24KQF9UjAGGlOowyWli+DRrJYXjhf0JqCLBGAAAgKQZSg0wAu0q\nUApDJKfJtFfuKK9pHFZtfwKqTMUYgCQZ7smkTdOwavsTUHWCMcAIuT0JAED1CMYAIzZNVaBppbpF\nWVT9eOFiIDAtBGMAgAmrYsB0gQmYJhbfAhiTaVxcp+qqFjwAgLWhYgwAUAJVGlatWgxMG8EYYIyq\ndKI77ZzIUxVlHlZd5rYBDEMwBgAomTJeVHNxCZhmgjHAmJXxBDclKlwAwNEIxgAAJVSWi2ouLgEp\nEIwB1oiTy7Vn6CfTYJLHDvsQkArBGACg5Na6euxCHpAawRhgDZVlaGQKVLqYRmsRWO07QIoEYwCA\nChnkAtvMTx6K+Lu/a/23D0IxkKqZSTcAIEXtk9rv7vvghFsyfVTkSUF3ON6/a1vH67Unn4jXX/An\nccyhn0VExBsi4rUtW+O5r/99NE96a8d722F4pp6Nvd0AZSUYAwBUWDsk33jZ++L5F1+L5198NfLL\nrojmfx4Tr/3We2KuNhMzzblY9/yrMXPp5yK78cZYf9y6+J9/86PIMmEYIEIwBpiYuUYeH9p9d9y8\n+/2TbsrUUPkiJc1mHu1Bz61w/MDii+/+04h39/jg//pJ4Uke7b2lmS8+BkiNYAwAUGJ53oq/zWYe\neWH6b7PHVOA3vHgkTv1//zte99rLcUzjtag3m/FafSbm6jPx4rpj4/+ceGo8d+zrF7ff3l6zczvN\nZh4KykAqBGOACbtoX2ueYPccQfpnXjHTqlEIw91BOMsisvmv12utY0iWZTHzk4fihD+6oOc2n7nn\n/pg748y4aN/BiGgF42azta3it2jmi19oNPOoCcnAFLMqNQBASTTzPBqFBFysEGdZRK1w5lavZVGr\nZTFTz+Lm3dsX5gvPnXFmvLZl67Lbf+2d74q5M86MiIibd2+Pm3dvj9r85+r1rGMaQrFanOcRjfmK\nciusW7kamC6CMUAJuL/x6rm9DFWX53k0m+3h0p1huJa1qsER80G4a2zz/l3blh1t8tzX/35JOH5t\ny9Z47tbblrx3/65ty87Lr9eyhe9d6xGS2+0HqDpDqQEA1lgxTDa65vbWssVh07WFRNoZPmfq2YrT\nL5onvTWO/OM/xTEPPxRv+I//G0c2/Zd47fQze76/va3uC3StKnQetVoWWZ5Ho9mqJBezcOtrwjFQ\nbYIxAMAaajbzZecL5/NzhbMsi+aIRkHMnXFmxMb1Mffr57uz9cDaIbleyyKfD8lt7aDcaOQdw70B\nqsKhC6AkDKcenGHUVEXHcOmuLluvtYYtR0Rf9xU+WrV4GL2GVXcrtrMV5luP8zDMGqgmwRgAYMwa\nzc5AXJyz208Ybus1p3iU+g3HbVmWLQT77pWrW/9u4RgoP0OpAUqmXTV2+6beVNYpuzzPe45cHvVw\n6TIpzkVua98fuZnn4Y5PQFkJxgAAI9IeOrxkQa1a6z7BrYWqBo+H4xw+vZxei3H1oz0XOaJzIbFm\ns/dnACbNUGqAEjLfuDfziimr7gWpimpZtqpAXHW1wi2fupl/DJSJijEAwBCWqxLXstaCVL2C8iAm\nOa1i/65tQ1+kW1zNurV6dbuC3GhGZMMulQ0wIirGAACr1GjmSwJxxPxc2wQrxCvJsqxwX+aWdixu\nNHMVZGCiVIwBSsxCXIsMLacsigGumOVGvajWWs8r7mWY+cYrqdda1eN8/n/tH1ue5y4qAGtOxRgA\noE/dFeIsi4U5tMLcYIq3eSr+6BrNWLjnM8BaEYwBSs5CXBbcYrLyPF8Iat2jfetjGjJdlmpx0aD3\nNx5EvWuRrmIuNsQaWAuCMQDAChrNxaCWRfRcZZnhtC8w1OdvbdXW+vkLx8B4ObQDVECqVeNPXPvD\nJP/dTFaxQtzWDmr1ukW1xq21SFfn19r3QBaQgXGx+BYAQCx/26WIiFqtdR/itRrOX7Yh1EXjWoir\nW/sWTxGtlb7b1ymazYjcLZ6AMVAxBgCS18zznvccrqkQT1Sta/5xMRabfwyMimAMUCEpDS224Bbj\nVgxVza57EU9iHnEZF9zqZZwLcS2nOP+41rWCdcMK1sAICMYAQHLyZSrE7cBVG9NK0wyvNf+483fT\nvr7RaOYqyMCqCcYAFZPCQlyqxYxDnucLizd1h+J6LZYELsqvXlu8B3KeL/29AvTL4lsAwFTruahW\nYVGnSVeIqzKEumj/rm0Tv0iXZVnUs9YFwyzrvM90o5mHax1Av1SMAYCp1Wwuv6iWCvH0qXcv0lWo\nIDdzw6yBlakYA1RUu1JTxUpTL5OuPjEdOhbVKmShWtYadttoTr5C3FalBbeWs1a3b+pX8TZPxQpy\ncXG1PM9L8/sHykMwBgCmRmsBps6v1WqtYLRYIVY5TEG9li27yFrrwog+AHQylBqgwqZpIS4LbrFa\nzWa+0He6Q/FMPXMf4oR13Oapa5h1m2HWQIRgDABUTN4VZLqHS0/iHsSrUfVh1EVrfV/jQWXZ4gWS\n4krWEa3RBO2qsoAM6TKUGqDi2lXjqp5gT0vFm/Hrtbp0FhGZ4dL0qbiSdbdW39J3IEWCMQBQeisN\ns6/PVyqbAg2r0B5hsNzq5Y3520AB008wBpgSVVyluj2vuMxDMFlbi1XhpYtota0UZKqiSvtpv8q2\nQnW/FleozqNea9WL26tY57E4H3mu4b7IMM0EYwBgohbCcCNfqPl2h+JiGC4GGRilLMsii8XRB7Va\nqy8u3Pap0OUazTzkZJgeFVmeAoB+VGmVaqtQp6l74ay5xuLtdIq9IetaRCvLsqm49+w0LbjVS9kX\n4hpELcuiXlt+0a48XwzKrX7cXhndKtdQRSrGAMBYtW6H03q80vDnWtYKHo1mLIQRVWHKonvRrlrW\nOdS6Vx9v5irLUAWCMcAUKvN846pUtBlMni8Og240844821whDBeHSFtRmipp99e5xvzc5LxzqHVb\nsf/PdS3mpbIM5SEYAwB9WRJ+C4pVsu5z/Sxr3VKpmbeHo2YLVbeU5guX8ULVuOzftS2pi2Ctof4R\nzfnFBPN8cYpAFp29O++Yp1x43MijWFouVprzPJ+KqQRQZkPPMf7a174Wf/zHfxznnntu7N2715Uv\ngJIo43xj84rLoXsOZDPPozkfdJvNvCP0zjXyhd9Zo1lYrTdfJgDP/7fWNT+4XssWqmtO7klBsZ/X\n69nCnOt6LXqubF0clh3R2tfawbnRXBzC3Wh07qPNZmH/7dq3zXeG/g1VMT506FDcc8898Z3vfCfW\nrVsXf/Znfxb3339//MEf/MGo2jcx7YNIs+uKeLPrQLTc45G+1nXl8GiPx/Ga7Vdz+91/CLv/UA76\nvkFfc/ILa2PJSfD842Yz76hSFS9ILJkDWXi+3FDQonb1q5a1nrQ/2z7xn2vkhkRDD8XKckR0VJe7\nh2Nn2dKLTxHze1UxQK8wbaHRNYy7qOOY0P248Ce8GMKXjhRZrGpP/JzZ9if2vafl4stQwfjBBx+M\n7du3x+te97qIiNi5c2c88MADfQXjMp8z/+eLry0cSLpPEDoOPj0ej/S1Hge7Xo/H8ZrtV3P73Se+\njR6v9fu+wV9rdd5+/xC3Ptf7ALxSKKe3dtX45t2THcJ50b4fTvT7V0W7bxcDbkRnyG00mh2v9dov\njxZwi7L5/8vzxWpW+/PFOcDLhd+m8HtUM/VsYvtg+3xrEudd7X+z/X9lrQvJ+ZLQXO+awxwRiwE6\nFv/mFxcBW+0lqbz7ceELeY/H7ecLF+Mmfc5s+xP73o1mxFPPvRQbN64vdcY7mqGC8ezsbJx66qkL\nzzdt2hS/+tWvjvq5N77xuKjXy3unqA15HrUsa83tmP/ltg8Exee9Hq/0vnG8Zvsj3H7MP+96vNJr\n/b6vitsft+7v0TGErOvFXif/3cG7M1xHx1Xv1u+8Xe8qy2+g/dp4t79x4/qYrOJfymnYO/p/rXto\nZFv3AlUrVnQL71tp3+wMuPPBdb4B7ZPsRjOPei2bn+fb+kbtv8lzjWbUaq3HzfnXsqz997q4Q1Xn\n51+O7U9+H9ywYZLfv2zH3eptv7gfZllt/l7LrX20vc/ONZod+3JExEzhefEYELH0mNB+vHDRa/59\nxWNJ93Flpdcqcd6X4PbH8b1rWRa/9V83RcSkjzXDGeniW/1WcZ5++oXSX02o1SKajc6rdRHdV++W\nf7zS+8bxmu2PcPuFakjx8Uqv9fu+Km5/puu1mVW8VrzK3fr5Lz5f7qp3RCsELHfRop8jTGe47gwe\njcIwjEaj2ZHVmnmzsMhJc4UFgfIej8fx2ui2/6Hd/xARMZGq1UX7ftjRL9pt6jTZn89q3tcxR7eQ\neJdWdAvv6xj5ECsq9vnivtE6bi0u7DPTY3+eP1eOZmN+e4Uff+vxsD/X6vT/yW8/j1//+vmYhCxr\nnag+9dTzR+1z45N3/bf76ys9HsdraW6/+9x76TGhZeF4035eOJZ0H1dWeq0S530Jbn8c37tWi3jm\nmf8swbHm6Fa6SDlUMH7zm98cs7OzC89/+ctfxlve8pa+PlvmHxgwOt0Bs/g8y7L5q96t1xauUh/l\nosVKwXtJuM6XniYstKYYonsMBS/On4ro/wJgWVW8+Wuq+3fdvSBVUd8V3SwWVmdeOIksXDRqr9bc\nfeGpuG+ktIrzNJipZ7F/17aJ73t5Prn9P7UVqiFF7ePLJI81wxpqPPO2bdviBz/4Qbz00ksxNzcX\n9957b5x99tmjahvAUbXmZC1G11ph9dt6LVsI1K3nixXrem0xmLS2s/z2iwuhRLRC88IqoV0rg5Z9\n9c9JrFJdhVWo8zxfqOQ2mvnC3Pfi77r1vt5/7LNs8aJMrdDPIlrBaHFF2sX+WcuyhaGHrW306IQA\nwNgNVTH+nd/5nfjIRz4S559/ftRqtXjPe94T73vf+0bVNoCRWqla3WuRk+6FiIq6q84r3Y+yTIG5\nHY7HeU/VslWHltx/t8fvbaVfU3EV5qVTBBYrujUVXQConKHnGF9wwQVxwQUXjKApAOVQDNDFIawR\nnbfWqNUiIh8uNDfzPLIJBKh29Xgc4bhdJe6cUzx+xYsPg4bf9nze7vDbHtocER2rMBvOTD/GefGp\nato/i7JdNANoG+niWwApaK9sXZufMDpMaO6+LVUxNJepylwWnQtedd7SaDXhN6L3LYgMbQaAdAjG\nACM0aGjuXm27GJqX3iu6cz7zKPQzrHrmJw9F/ReHo/HfNsfcGWcedVvD6r6Hb6OZLwTdo92nd9Dw\nCwAQMeTiWwD0r12BrGWFBcLqtY4hx7Xa4rzm7nplcfGn7oWhiguBNZt5122BVg6AvRblqj35RLzh\n7N+NE/7o9+P1/+NTccIf/X684ezfjdqTTyx5bz+LbBUXJ2vmeUe1vNHIFz7faHauEt6r+bWscwG1\nmfriYmu1WufCVrCW2itRs9T+XdvWfJoFQD9UjAFKpFhp7r7HdL+V5sLtmiOi8/ZTc43F2091337o\n41862HqSRVz6ka1x/J4r4phnmhFvfUdkkUeWR8TTr8XcpZ+Lp//qK/HVu/41moUFyrq32Wj0Hurc\n7KqGd2ff4r+vVms9bzRVfgGA8RCMASpipdBcr7WCZLO5eOupXpXW9pe7X1/ItHnEdd/+WcS7/zTi\n3T0a851/Xfq5rm0e7X6+xfcXV3lecg9fC10BAGMmGANMgeLtp4q3noronN9cDNDdt6Iqjjz+78e+\nHDOHfhbt+nKeZZFHFnkWkeURtXe8PX7y0rGR550ht7jN7lsaFQNvdxut8sy0MIT66Pbv2mZ1aqB0\nBGOABLQXBSsG6O5bUdVri/MiZ37yUJzwxb/qub1nLro/PnnGmR0ntx3DmxtuaQQAVIdgDJC45RYK\nmjvjzHhty9Y45tDPlrz/tXe+a2F16uLnVIAAgKoSjAESttKwz+e+/vfx+gv+pCMcv7Zlazx3620r\nbktAJkVWoh6M4wVQNoIxQIL6OYlvnvTWOPKP/9T3fYzbzB8EAKpGMAZgRXNnnNlXIIZUqRSvngtp\nQFkIxgCJWYuTeMMkAYAqqU26AQCsjUnMgdy/a9vCbZoAAMpKxRggAZMc6mmoJNPKglujYYQJUAYq\nxgAAACRNMAaYYmWpaO3fta0U7QAAWI5gDDClyhKKi8w5ZlqUcf+qOscHYJIEYwAAAJJm8S2AKVTm\nSpaFdgCAshGMAQAGUOYLT1XnwhkwKYZSA0yRKs17NJ8QACgLwRgAoA9VuvBUdS6cAWvNUGqAKVHF\nE/b9u7YZMgkATJyKMUDFVb2K5R7HAMCkqRgDAByFizdrz4gSYC2pGANUWNWrxUXmFAIAkyIYAwD0\nME0Xn6rIVAtgrRhKDVBR03iy6B6mAMAkqBgDAACQNMEYoGJSGNppvjFlkMK+VhWOCcC4CcYAAAAk\nzRxjgApJqXrlVi1MUkr7WlVYgwAYJxVjAAAAkiYYA1RAqnMd3aqFtZbqvlYl5hsD4yAYAwAAkDTB\nGKDkVLBUiACA8bL4FgBAWHCrSizOB4yaYAxQYk7UF1mRFgAYF0OpAYCkma5QTRbnA0ZJMAYoISfq\nvZlvDACMmmAMACTLRajqc7EMGAVzjAFKxkn60ZlvDACMkmAMACTJRajp4WIZMCxDqQFKwpDOwRlC\nCQCMgmAMACTFRajp5WIZsFqCMUAJOFFfPSfCAMCwzDEGAJLhAtT0279rm7nGwMAEY4AJu3n3tsjz\nSbei2iy8w9EYlZEWxwRgUIZSAwAAkDTBGGBCZupZfHffBybdjKlivjHL2b9rm2pxohwTgH4ZSg0w\nAft3bYvMudpYGEIJAAxKxRgAmErmFRNhxADQH8EYYA05UV87hlCmzb5GN8cEYCWGUgMAU0UgphdT\nLYBeBGOANeJkfe25nykA0A9DqQHGzJDOyTK/MB32NfplWDXQTcUYAKg8gZhBGVECFKkYA4yRClZ5\nqBBNJ/sYwzCiBGhTMQYYEydb5WPhnekiFDMqqseAijHAiDlZLz/V4+pT6WPU9ClIm2AMAFSGC0+M\nmwtnkCZDqQFGyAl7dRg6WT32L9aK4wOkR8UYYARUsarJ0MlqsH8xCY4PkBbBGGBITtqrz9DJ8rJ/\nMWmOD5AGQ6kBhuCEfXpYsbp87F+UheMDTD/BGGAVVLGml7mFk2W/oswEZJhegjHAgJy4Tz8nv2vP\nxSaqxDECpo9gDNAnJ+7paVeP5xr5pJsylexTVF2x/wrJUG2CMcBROHFPmxPf0bNPMY0cK6DaBGOA\nHlSz6GaFEbQBAAANC0lEQVT+8eDsR6Sou887bkD5CcYABU7iOZruqpBh1p3sQ7DUckHZsQPKZehg\nfPDgwbj88svj05/+dJx33nmjaBPAmnIiz2oZOmn/gdXotc+kehyBMhgqGB88eDC+973vxemnnz6q\n9gCMnRN5xmGaK0L2F1gbK+1rQjOM11DB+PTTT4/t27fHZz7zmYE/m2XDfGeA/ty8u7wn9O3joOPh\ndOrV9y7aV46T2zLvG4yO48z0GGafLctxh+k1DceaLM/zoS9nf+Yzn4l3vvOdfQ+lbjSaUa/Xhv22\nQMV8aPfdC4/nGs2YmT8OzDWaEREdz2cKx4jv7vvAGrYSAIDU9FUxvu++++Kmm25a8vUDBw6s6ps+\n/fQLpb+akGURGzasj6eeej6Gv3RACvSZo7t59/tX9blf//r50TakJPQZBqXPMCh9htXQbxhUVfrM\nxo3re77WVzDesWNH7NixY2QNiohS/8CK8rw6baUc9BkGpc8wKH2GQekzrIZ+w6Cq3GeMZwYAACBp\nQy2+deutt8bBgwfjF7/4RfzLv/xL3HvvvbFnz554+9vfPqr2AQAAwFgNFYw/9rGPxcc+9rFRtQUA\nAADWnKHUAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnG\nAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0w\nBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmC\nMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkT\njAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEia\nYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDS\nBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACS\nJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQ\nNMEYAACApM0M8+HHH388rrjiimg2m/Hiiy/GpZdeGmedddao2gYAAABjN1Qw/vznPx/nnntu/OEf\n/mE89thj8clPfjIefPDBUbUNAAAAxm6oYHz99dfHscceGxERGzZsiCNHjkSe55Fl2VE/28dbJqrd\nvrK3k/LQZxiUPsOg9BkGpc+wGvoNg5qGPpPleZ6PYkNXXXVVRERceeWVR31vo9GMet30ZgAAACav\nr4rxfffdFzfddNOSrx84cCDyPI+9e/fGE088ETfccENf3/Tpp18o/dWELIvYsGF9PPXU8zGaSwdM\nO32GQekzDEqfYVD6DKuh3zCoqvSZjRvX93ytr2C8Y8eO2LFjx5Kv53ken/3sZyPLsrjxxhvjmGOO\n6btRZf6BFeV5ddpKOegzDEqfYVD6DIPSZ1gN/YZBVbnPDDXH+Bvf+EY0m83Yt2/fqNoDAAAAa2qo\nYHzLLbfEpk2b4vzzz1/42pe//OV405veNHTDAAAAYC0MFYz/+Z//eVTtAAAAgImwNDQAAABJE4wB\nAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAM\nAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRj\nAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYY\nAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTB\nGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJ\nxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRN\nMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBp\ngjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRtZpgPHzp0KL70pS9F\nrVaLl19+OS688MLYuXPnqNoGAAAAYzdUML7lllti165dsXXr1njiiSfigx/8YOzYsSOyLBtV+wAA\nAGCshgrGX/3qVxceP/nkk3HiiSf2HYrLnp3b7St7OykPfYZB6TMMSp9hUPoMq6HfMKhp6DNZnuf5\nMBs4fPhw7NmzJ5555pm44YYb4rd/+7dH1TYAAAAYu76C8X333Rc33XTTkq8fOHBg4fGhQ4fikksu\nibvvvjuOO+640bYSAAAAxmTVFeNmsxn33Xdf7Ny5c2H49Ac+8IHYu3dvbNmyZaSNBAAAgHFZ9e2a\narVa3HDDDfGjH/0oIiJmZ2djdnY2Tj755JE1DgAAAMZtqDnGjz32WHzhC1+IWq0WL7zwQnz84x+P\nc845Z5TtAwAAgLEaevEtAAAAqLJVD6UGAACAaSAYAwAAkDTBGAAAgKTNTLoBZXbo0KH40pe+FLVa\nLV5++eW48MILY+fOnZNuFiX2+OOPxxVXXBHNZjNefPHFuPTSS+Oss86adLMouYMHD8bll18en/70\np+O8886bdHMoqa997Wtx//33R71ejy1btsTll1++cLtE6OW5556LK6+8Mh5++OH48Y9/POnmUAH7\n9++P73//+1Gv1+Pkk0+Oq6++OtatWzfpZlFSzWYzrr322njkkUdiZmYmNmzYENdcc00cf/zxk27a\nwFSMV3DLLbfErl274pvf/Gb89V//dVx++eVhrTJW8vnPfz7OPffcuO222+Lqq6+Ov/zLv5x0kyi5\ngwcPxve+9704/fTTJ90USuzQoUNxzz33xG233RZ33HFHHD58OO6///5JN4sKuOyyy+LMM8+cdDOo\niEceeSTuvvvu+Na3vhV33nlnvPLKK3HgwIFJN4sS++lPfxqzs7Nx5513xu233x7HHntsfPvb3550\ns1ZFMF7BV7/61di6dWtERDz55JNx4oknujrPiq6//vrYsWNHRERs2LAhjhw54mIKKzr99NPjuuuu\ni+OOO27STaHEHnzwwdi+fXu87nWvi1qtFjt37owHHnhg0s2iAr7yla8YuUTftm7dGnfccUccc8wx\nERFxwgknxDPPPDPhVlFmp512Wlx33XUREfHqq6/G7OxsnHjiiRNu1eoIxkdx+PDh+PCHPxyf+9zn\nFn7p0Mvxxx8f9Xo9IiJuuumm+PCHP+xiCitav379pJtABczOzsamTZsWnm/atCl+9atfTbBFVIVj\nDIOo1+sLQ2Aff/zxeOCBB+Kcc86ZcKuogn379sX27dtj8+bNle0z5hhHxH333Rc33XTTkq8fOHAg\nNm/eHHfddVccOnQo/vzP/zzuvvtulR1W7DN5nsfevXvjiSeeiBtuuGECraOMVuozMCgjUYBxeuyx\nx+Liiy+Oq6++Ok466aRJN4cK2L17d1xyySWxZ8+e+Nu//du48MILJ92kgQnGEbFjx46F4a9tzWYz\n7r333ti5c2dkWRZbtmyJ4447Lv793/89tmzZMqGWUhbL9ZmI1snqZz/72ciyLG688caFoUjQq89A\nP9785jfH7OzswvNf/vKX8Za3vGWCLQKm1aOPPhqXXHJJXHvttQtTCqGXn//859FoNOJtb3tbrFu3\nLnbs2BF33XVXJYOxodQ91Gq1uOGGG+JHP/pRRLSGsc3OzsbJJ5884ZZRZt/4xjei2WzGNddcIxQD\nI7Nt27b4wQ9+EC+99FLMzc3FvffeG2efffakmwVMmfYdNa6//nqhmL4cPnw4rrrqqpibm4uI1mJc\nmzdvnnCrVifLjcfq6bHHHosvfOELUavV4oUXXoiPf/zjlR0zz9o466yzYtOmTR3D7b/85S/Hm970\npgm2ijK79dZb4+DBg/GLX/wijj/++PjN3/zN2LNnT7z97W+fdNMoma9//etxzz33RK1Wi/e85z1x\n6aWXTrpJlNyRI0fi4osvjldeeSUeffTReNe73hWnnHJKXHHFFZNuGiV15513xnXXXRennHLKwtfe\n+973xqc+9akJtooyy/M8rrvuunjooYeiXq/Hxo0bY+/evfEbv/Ebk27awARjAAAAkmYoNQAAAEkT\njAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGn/H7uznZLzJ/zTAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4cc4f6e3c8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(15, 8), dpi= 80, facecolor='w', edgecolor='k')\n",
    "ax = plt.subplot(1, 1, 1)\n",
    "ax.errorbar(xpts, y_pred[0], yerr=sigmas[0], capsize=0)\n",
    "ax.plot(x, y, \"ro\")\n",
    "ax.set_ylim(-3.0,3.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "At every location where there is a \"known\" point the standard deviation goes to zero, because the function at this point is not \"unknown\" or \"random\" any longer, but fixed and known. This is the effect of using the conditional probability $p(y_*\\,|\\, X, y, X_*)$ as described above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Let's draw several values at once rather than going single point by point as we did above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ -7.66697664e-06  -5.63595765e-03   4.19316345e-02   2.02471666e-02   6.46090979e-06]\n",
      "[[  9.99999997e-01   1.65296628e-01  -3.73627090e-07   1.19843900e-12   3.82424663e-16]\n",
      " [  1.65296628e-01   9.98338443e-01  -2.74559569e-04   8.80965650e-10   2.81118181e-13]\n",
      " [ -3.73627090e-07  -2.74559569e-04   9.92508018e-01  -3.50450933e-03  -1.12241541e-06]\n",
      " [  1.19843900e-12   8.80965650e-10  -3.50450933e-03   9.98338443e-01   8.62930563e-02]\n",
      " [  3.82424663e-16   2.81118181e-13  -1.12241541e-06   8.62930563e-02   1.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "x_more = [-2.1, -1.5, 0.3, 1.8, 2.5]\n",
    "m, s = conditional(x_more, x, y, cov=k)\n",
    "print(m[0])\n",
    "print(s[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.5128756 ,  0.52371713, -0.13952425, -0.93665367, -1.29343995])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_more = np.random.multivariate_normal(m[0], s[0])\n",
    "y_more"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "x += x_more\n",
    "y += y_more.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "y_pred, sigmas         = predict(xpts, x, y, cov=k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-3.0, 3.0)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8YAAAIICAYAAACsBPexAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAMTQAADE0B0s6tTgAAIABJREFUeJzt3X+8ZWV9H/rvWnsGJTAKzqCiSCNO1KiZjPclEPWaONQm\ngxZ7ja3W66UlhMGXTWkxlcGY6DUaMA4SbyRIBROtv1CstQT0hVUHRHtfhQTDnSaURLShYJQRZIAg\nDpy91v1jzz5n7X1+7rN/rbWf9/v1Us6Zc84+6+z9rLWfz3qe5/tkZVmWAQAAAInKp30AAAAAME2C\nMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApG0Y5oeLooiLLroobrnlltiwYUNs\n3rw53vve98aRRx45quMDAACAsRpqxPhb3/pW7N+/P6666qr49Kc/HYcffnh89rOfHdWxAQAAwNgN\nNWL8ohe9KF70ohdFRMSjjz4a+/fvj5e//OWjOC4AAACYiJGsMd6zZ0+ccsopsXXr1njlK1+56ve3\n28Uofi0AAAAMLSvLshzFAz366KNx/vnnx8/93M/FmWeeueL3/vCHD0WWjeK3jk+WRWzevCnuu++h\nGM0zxKzTZhiUNsOgtBkGpc2wHtoNg2pKm9myZdOyXxtqKvW3v/3taLfb8dznPjcOO+yw2LlzZ3z+\n859fNRhHRK2fsKqybM6xUg/aDIPSZhiUNsOgtBnWQ7thUE1uM0NNpb7jjjvi3e9+d8zNzUVEpxjX\n1q1bR3JgAAAAMAlDjRjv3Lkz/uqv/ire8IY3RKvVii1btsQFF1wwqmMDAACAsRsqGGdZFm9961tH\ndSwAAAAwcSOpSg0AAABNJRgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwA\nAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMA\nAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgD\nAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEY\nAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnG\nAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0w\nBgAAIGmCMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmC\nMQAAAEkTjAEAAEiaYAwAAEDSBGMAAACSJhgDAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkT\njAEAAEiaYAwAAEDSNgz7AJdffnl8+ctfjlarFccff3xceOGFcdhhh43i2AAAAGDshhoxvuWWW+Ka\na66Jz3zmM3HVVVfFwYMH4+qrrx7VsQEAAMDYDTVivH379rjyyitj48aNERFx9NFHx/3337+mn82y\nYX7z+HWPr+7HSX1oMwxKm2FQ2gyD0mZYD+2GQc1Cm8nKsixH8UB33nlnnH766fHpT386jjvuuBW/\nt90uotWyvBkAAIDpG3qNcUTE7bffHuecc05ceOGFq4biiIgf/ejh2t9NyLKIzZs3xX33PRSjuXXA\nrNNmGJQ2w6C0GQalzbAe2g2Dakqb2bJl07JfGzoY33bbbXHuuefGRRddFNu3b1/zz9X5Casqy+Yc\nK/WgzTAobYZBaTMMSpthPbQbBtXkNjNUMP7xj38cb3nLW+KSSy6J5zznOaM6JgAAAJiYoYLxtdde\nGwcOHIjf+73fm/+3l7zkJfHmN7956AMDAACASRgqGL/uda+L173udaM6FgAAAJg4paEBAABImmAM\nAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRj\nAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYY\nAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTB\nGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJ\nxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRN\nMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBp\ngjEAAABJE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJ\nE4wBAABImmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpQwfjBx98MM4999x4\n6UtfOorjAQAAgIkaOhj/5m/+Zpx88smjOBYAAACYuA3DPsAHPvCBeOCBB+KP/uiPBvq5LBv2N49X\n9/jqfpzUhzbDoLQZBqXNMChthvXQbhjULLSZoYPxpk2b4oEHHhjoZ570pCOi1WrG8ubNmzdN+xBo\nGG2GQWkzDEqbYVDaDOuh3TCoJreZoYPxevzoRw/X/m5ClnVe2PvueyjKctpHQxNoMwxKm2FQ2gyD\n0mZYD+2GQTWlzWzZsnxwn0owjohaP2EREbv2XB8RWVyx++W1P1bqpSzr376pF22GQWkzDEqbYT20\nGwbV5DYztWDcBHPt4lBAjphrd17hDa0sLj9vxzQPCwAAgBEaKhgfOHAgzjnnnDh48GA88MADcfrp\np8ezn/3seMc73jGq46ulsy+6PubaZWxoLcwHF5ahmfrP56U+dn4DAMy2oYLxUUcdFZ/4xCdGdSyN\nNdcu4+yLrp//XCca6q16vnZng6z2/UIyAMDsMpV6xObaZZz5+3t1oKGmlprxsVZuggEAzCbBeIyM\nMkF9VAPtKHRDsnMbAKD5BOMJEZJheoYZJV7LY0cYPQYAaLJ82geQmv6pmMB4dUMxAAAsx4jxlBhB\nhvGa5A0o06oBAJrNiHENnH3R9UaRYYSmMUpsNggAQHMJxjWhUw3Dc5MJAID1EIxrRDiG9avLWmLh\nHACgeQTjGtKxhmZzkwsAoFkE45rSsYa1q8toMQAAzaQqdc3ZIxWWV/ebR85fAIBmMGLcAEaPYTGj\nxAAAjIpgDDBGbmwBANSfYNwQOtfQoTgdAACjJhg3jFBAypo6fdqNLQCAehOMG0gnG2D6zr7o+jjz\n9/e6HgPADFCVGqi9WQkeqlQ3V7UNzrXL2NDKej7vfr3/a15rAGgGI8YNZlo1KWjq9GlmQ3dUeD1t\n0OweAGgOI8YN1+14GZWAZnDONsOoAm1/OPa6z47uTbvuDIHlPo7wugM0gWAM1JbRYqahP/CMipsi\nzdc/pX4t5tplnPn7e2NDK/PaA9SYYDwjrF1klph+yjRMqt25XjfTKG6YeO0B6kswBpiCbgf5C3te\nPeUjIWJ8o8TLMXrcHKO+YWJ6PUA9Kb41QxR6YRaYPs0kKWLISsZ9PfK+DVAfgvGM8SZLU6UYUOba\nZbxm9zXTPgymKMV23xSTukmX2vt2t9K7PcCBujGVGpi6paaxlmVvh7T/cxjWpKdPL8e06nqZVlib\nxfXHK+3/Xf0eVbyBOjBiPKOMQtBERVFG+9AITbvo/K+r+vlcu4x20fm+WQjMu/Y4XyetblP2Uxs1\nZGmz0g7s/w00kRFjYGoemyvizRffEIcy7vx/16Kbh9tFRJbVJ+BQbzrcrKQuswiaOHo86ueu+nhN\neh6A5jJiPMPcdaWuyrKMs/bsjTe9/4aeMJxlEfmhq9KGVtbTwap+3soj8qz6eAsfF0U5E6PIpMls\nH7qa8h4+7jbblOcBaD4jxsDYbbj5pmh9945on7A17nrWC+JPvvQ/oigWf1/rUNotYuVgm2VZZFlE\n0S6jlXeCcc+o86GPy7KMLJvuyM9aWWc6fnUZDaR+6hq86r61k3MKmCWCcQKaOCWL2ZDffVc84Yw3\nxsZ9t0YZEV99/j+MD7/iTXGwdVjn63nEFeftiDe9/4Z1r/eshuR+7SKilRs9plncJJmspoS7urSL\n/oJak/690/77gdklGANj0w3Fj7Y2xh/+yjlx43N/MSIiTr3rpvg//nB3bNzQiojejs6wIzf9I8jt\nJUamSUtdRwNXIwiwlGm2i2nfRKjLzQFgNlljnAhrdJi0DTffFBv33Rp//7gj4p2vfVfc+NxfjE2P\nPBi/858viH/1uffG4d/68yV/7vLzdgzV6cqyLPJDU7L7Z1HXfe2x83T06lZ9mnppavsY5Fqx4eab\nIv7Df+j8dx26FabrdG2yFh/qY5bORyPGwFi0vntH3Hvk5njXr74z7tzyD+Jp938vfvfzvxtPfXD/\n/NfnTjp5yZ/tjgYMP3qcRVmW86PG7SIiV8GaBjFCxnJWW39cXcoSEXFURDy2bXs8+LFPRXHcM1Z8\n7GmPDANMg2AMjMW9Tz8hLnzdBXHPUU+Nn/nB38T//YXfiyc+8uD819snbF31MS4/b8fQIzqd4lsL\nP180oDCXKbSj0dTRQMZvVkY3ImL+5t+uPXvnK/T/05dvjdZHLo/DN/6DePzPPTmOPPhwPOWBe+Jp\nf/3teMIZb4wDX71x0eNMa+3werlpBNM3azfRBOOEeBNhUv7+kcdiz9/kcc9RT41t/2tf/M7VF8bh\nj/1k/uuP/fwLlx0t7tcNx6PQyqNn9FhhLprAtZuq6pKQaoCt1lO46vo7Ip51SsSzFv/8UQ/fH8/6\n46/HX9zbjqLsXBezLJupzi3AegjGwEj95NG5+MPP/X/xvR8+HBvbj8X537m2NxRv2x4PfvSTUzm2\nLMuilfdOrWb2zNJoIKPXxBGObhhut8tlN7PL5v8v4lef+FAcec1/ikc3PC4e2Xh4PHj4prjniU+J\n7x91bPzoyCfFLT9sz/9cu4jIGrrExE0jYJQE4wSZpsm4FGUZv/EHN0YZEU/bckS87Y3/Wxz87V+J\n+yv7GK91pLhqVGuOIxamVufZwrTqiPpNrdbho1+3/X9hz6unfCRMUjUM98fX1qESqu0iotXK5q8X\nG26+KY5+5xeWfLz7f+qJ8eeXfiYu+cuD81Ovu/+da3eujQArmdUb0KpSAyOza8/1UUbE0ZseF//u\n9dvjyMM3RkTE3Eknx8F//sZ1heKqYStWV+V5FnnlCtgu6l+1mtVZV8xymlI5tSzLKCp37aqtOc8X\nwnBEdx/3LDZUQnFE55r72LbtSz7+kT9zQmx/1Uvjj88/JVp5dyr1wtf7bxg2QVNeW6DeBGNgJL71\nNz+cH3X4V695QRy96XHTPaA1yPtGiE2tpu7m2mW8Zvc10z4MxmB+unTRG07zbCEM54eCcNXl5+1Y\ncmbJgx/71KJw3L+U5Yrdp8QVu0+J1qFh4v7R4s6xNCMcA5MxyzegTaVOlGmajNKuPXvnQ2WeRTzr\naU8c2+8a5bTq5dRpWrWlD2tjtIiV1H1dcVGUPWG4Wks/n0+rvR3R/lHiRY953DPiwFdvjI1/dlMc\n9cPvxYFjnh6Pnbj8FnndcyirbHHXObZDv71G18Wl6NcAwxKMgaFUQ/HLX/j0+Be/8pzpHtA6tfLO\nKE1Z2u+Y+tu1xw2TWVINxdUq0ctZLRRXzZ10csSWTTF370OLFykvobrFXbUWg+siMOtMpU6cdTkM\n4+yLru8ZWXjDP/yZif3uUa43juh0BvNl1tlRb7M8rYvZVJZltPvabLfmwbRGZZeakp3nWc+a5uo+\n8HWlXwPjM+vvt0aMgYF1Ox3VztG7zzwpNm5o9r22pSpWt4vpV2k1RZClaBerq1tAqq4jruqOEhdr\nGNKd9OvdvS727wOf2wcemDHN7sUCU1NW1qHlWcRxTz5yugc0QtWK1WVp9BgYXtm3drc6MLyWUeJB\npk8PY7nZOJ194Bc+LxQrhGSkMhNDMGZ+1AHWojuNpqdQzJRGVJerxjoK1YrVNZ41mLxZn9a1Gtfv\n5dWtbVRDcWebpPoWslopHEcsrl5dt6nVzgtgPUylBtakv5PR7Qe18s6WH7OqOn0wYnqVWU2bhebp\nD4xZ1rl2DnoNqdt5n+dZT/XqdhHRMrUaaDgjxsBAqh29PJteoZiqURfiquqfPtgu6jc6kqJUpnXR\nXEVRNmqUeFD91/667QNv1BhGo26zb8ZJMCYivIGwsupFsazBFOpJq1MHUCBkKdpFr2l25Lo3zvpr\nE6znJuKk1hUvZZClKtW/rHDjEGgoU6mBFfWG4oW1xbM+hXotpjWtOnXdNjmuWQKwXv0FtvI8jSJV\neWUf+KKICNOqgQYyYgysWW/BrXqFknFOp65q5QujI6ZVUydm/kxv5HypbZhaeW8Rv0FNc7S4ai3X\n1kX7wNfoZoDZFLA+KZ47Rozp0T0B6vBmzHQtdTGUAQ91APPeojOTlmohrtTeoGmGsix7gmB3H/S6\n3Twct6X2gQdoEsEYWGSlqap51qlIWsdQ1j2mcQeobgewW2EW6iTVG5zTmmLfP3U6z7IohlzfXLfX\nbpBra34oGXfDcVGWQ42cA0yKqdTAQPRvFuRZ9E0flJIhZQJgR165MBbF9AtyWWYAg0mpEnWVYMwi\n3kDStdR6kp7tmfL0pgeupDOtutIBLCcbjlNa/5PqmzT10ymwtdAWsyx6tnRjsTqtOQZYjks5sKKe\n7ZlisC08pmVShbiW0u0vK8rFNKV0g3PSN4g6RfcWPm/l2chuGNal4NZy1nNtrT41rotAnVljDETE\n0uvzqtszRRgtXk11a5aijMhDJ3BYqYQ76m/a04GbKs8iiujcTGgXEa0pbuWU6vp7GETKM7SMGLOs\nlKZpsrQm9wOnMbKdZ9n8lMr5/TzHLKWRQQY369fxSXXgimKh8nQW45k63YTZOOvRv5XTNCr5A6yF\nYAws27lUS2pw1VF1T99wUr5rTb1Ur4Wp11pYT4BP+fkCmkMwBlbU7c7Ufe1b3fSPKI1zbZ1RY1ai\nfaxP/znbPaeFvPWb5HVxOc4HWNqszzBaC8GYFXkDmW3LXQT7K1E32bQKcWWVadUR3YI9Rj9hVMY5\no6Aoyp4pv60xjhI39abj+gpx9V4XXRKBOml4lxcYh25nJQujI8Pof+6srVs706hHy03OtenevOpf\nRuI6ODrV57IoFTUD6kNVakjUUlWoI3orUXdHi5s6olEnWbZww6Eoyp79j0dFxVVSMK6A39mfeOHz\napX5cXFtnd4ex66XsGC5PmFqjBgDPXr2LZ6RUZJp7mvclWcL+3kWZSccs5g1TkxL/9TpfEauf+M0\n7LXVHsdAnRgxZk3cWZ0dq4UOeW08siyLVhbz04O7z3NZljHK/nd3yqxzlaW4lvfqD2PdmR2TuCno\nNejcMGzPz6SJyKe4xzGAEWNgEeMk49VblKs0UnKIdcWsZpRtZNHU6SyiNYYlDiyvWoyrjMkX47L2\nntSZpdVLMIaErLVTmVWuDOvZs7KO6vR3LFWxGiZFGOjoP+/Gse4/BcNeW/uLcQFMi6nUrJkpmunQ\nPRy/TmewjM7/L+hMrR7+FWjSlFkhjdWMqo1UiwtGxKLzbxIU3FpanvUvMfFOBEyWEWOgR6dIlA7J\npOR5FtWBqqJUhIbJSHEKXWc/8c7Hedb8fdpnSZYt3JS17zswDd4SIAGrdYCrHZBZz8R1qFBdlWVZ\nzxTOshzd1OomTJm1rpjVDNtGyrJcsgp8nmcTvwlYpyUdddO5Fi58Pqlp1U24TsI4eP9dTDAGeqcW\nHuoomu43fYXCXIxZCqGgKBaucVnWW/yO0RnFTcfqjQqXPmDSvD0wkBQ6Uakpy1IHpEZaef9+x9M9\nnnFyt5rVjKKNVH+6NYVRYganDhowDYIxzLjVOpYpVgGt83TGTsXqyqjJkI/nZharmbU20tmGqbI8\nJKY/SmwGzmCyrHdZzyRmzqS45p40aevLG7oq9Yc//OH4yle+Eq1WK7Zt2xZvf/vb3Y2FBjFaXF+t\nvDNi3H2J2kU5EyMp3pAZh254WrwN00IVeMarG/6HPcezLIs8yuje0y3L2a9/AUzfUMF43759ce21\n18bnPve5OOyww+LXf/3X4ytf+Ur88i//8qiOD5iA/i5jKiMbl5+3o9bTeTvFaMr5jn5ZRqWjONh2\nJk3avonpqFMbGTRYdQpsLXxevaZN+2a90eL1qd7MKMqIzF1cYMyGCsZf//rX45RTTonHP/7xERFx\n6qmnxg033LCmYOzOX7N1Oy1X7PZm39Vt03Vq27v2LB/6qlPT8rx3lKVOf0Pqqp3Dame/XUS08sE7\nitN+bbttsk6VwenVxDZSvX7leedcGVV191GY5nNax/em9ZjU69n052lUZqXdMDmz0GaGCsb79++P\n5zznOfOfH3PMMXHPPfes+nNPetIR0Zr2gp9VLGSG7qtbLvH5ch+v9H3j+Nq0Hj9iy5ZNQa/Nm+v0\nnPSMmxz6b+c1Lec/Xvx9ab2u3b992uf26t+X553Xrbv1TG9HcfXfPdcuYteeG+ILe06L6Vnq+Y4Y\n7bVp1h+/+/noH78pbaQsFy8DybLOv+VZt39RLPq56by+9XivnOR70xf2vDpes/uamGsXMYpzo9ra\nO9Oq13bNHPS1r0f7r5d69WkYVue8LGPDfA4b3XWx21aa3GaGXmNctdbiCD/60cO1vpvwyMG5+cId\nc+2irwBEUZmWVf17+//2SX5teo9/770PBR1Z1rkY3HffQzVat7v8a1oeauN51n93r0zqdb1i98tj\n157rY3rn10pf6/2+LOu8DRWxEAK6iqJ6rVr58V+z+08jYrIzPjrP8XLH1P95/a990338tT7m+o9x\nGteAQdpIu69qYOvQOuLODJn6vL7dc2ya19TpvTeVff/t/3ilr/V+X3VWU2cZyVoeY+2PX4f2Xzf1\n7NMwvNXOy/Vf++6776FGtJmVblQOFYyf+tSnxv79++c///73vx9Pe9rT1vSzdX7CWnnvaHb1WDsX\n5s4/FEX/xTktc+0ydu253tqpPkuNZEzaauvzyrJcuN9XacPdtXDTPn5W18qzQ9V3O58XlTyw1vXH\nXmeW0w2pdbq+968jrpr2OuLl1Okcq8N703p16y1096Qe93rjpj5P49DkdsNkleXCf5vaZoaaz7xj\nx4742te+Fo888kjMzc3Fl770pXjFK14xqmObmo0b8vm1Ta18+f30inLhDma7KOenN0LdVZtqXTuU\nrK762lVfxXYRq16PJrlFT50LnLG8urWRdqVCe55NfwumlSi41bmpMspaAnnlejfu9ca2s4E0DTVi\n/LznPS/+2T/7Z3H66adHnufx4he/OH7pl35pVMdWC1mWRZZFFIfesLtvxO2idypjWfZOKEh9NJl6\na+qdvHEZ1RYj09RqZVFURtSquXjQCtajJBSzXkVZRtm/9VLWadv5/B1rbQtgLbwfr27oNcZnnHFG\nnHHGGSM4lGaori9u5d1iDZ29RctYCBzVaY1CMpPUvfCtdqe+WrqE2ZBnWRRRRp5HlD37H0dkWeez\n/pBcpy16qJ/uqPG42kf/zaiyLBfeR/tCcXcdcVHzjp1zabxaeb0qjgOzo8YTkZolz7P5oBzRu26z\nOuU6ItZcpKwpTDlqjmrby5z9i4x66t+05FkWrb6/o/vSt4uFwkXd9jCOKbOzcF3ohLRyyc8H+dos\nGMfr2b2JV32+2kXvbIfqUibLPog4tN64pyjqbJ1rwPToGo9JK8/mp133v5evZf0fjFt/F9OauNnV\nXythqZC8UjjecPNN8bjPfCo23HzTmn5fHadrVcNXZ8p5pcJxu5w/3rnKx53nZ+Exqp8P8rWex2yX\n88950Vebou6BetAbKCu1m1179i77PHebaiuvTpmmiS4/b8dY3lf6Bx/GYZJr7GGcZuFG9aSMdLsm\nenX32euvHhuxcCE3zZpRWctFr9vuOls0mUydimqthFbeaQfV+ghd7XYZZ71vb2RZxL9/68tj4999\nL55wxhtj475b57/nsW3b48GPfSqK456x6PdM4423GiSLSrGHdtFb+KHn+ts3DXM9Z8F6V7iWlR/q\n79D3BO122XP3qijLRTezpmW16ff53Xct2W7uueIT8Ruf+86KFUur2y41bYTYzcXJqb5/leWhcx9g\nSILxhFQv4nkevQVyXM+ZgOravYb1NydqFgpxrSTLsmhlndGQVt7tVHa+Nl8noYzYteeGePaBu+JZ\nx5wcz9z2lPjpH/5tPPWBH8RR+26NJ5zxxjjw1Rt7Hneta9vXai2Bt39UuuhZsrL8Y2eH/q/7PXne\n+bd20Vtgsfu3dH9P9fPWCl/r/7nqY+Z5RBx6zrvn4VLHWg3Q/X/bXN8Cy6JYCNGTGnFeae3xE854\nY2zYd2v84AlPjtuf9ty4/WnPif9x7HPjf1757SjzxRPVmhyGqYfltvECGIRgPAV5lkWWl0sWj5hm\n9dhhjLtAC8Or9peb2MYYvf6q+3le2X+wLONvjnpG/M1RvSPDj3vsJ/HkB/bHk674euz7UTvKsnea\ndnXkpj+kVb/Wv5xk1IE3zzofd3+2Gr76Q21eKapYLbA4KtXHzLOsc1zt3gKOEZ1w3Z1d1H/TYqX5\nHdUbrP1TuHumnI44QJdlGbv27J0P/K/9pWfFD/7HnfGDF/yfcfcvnh+PHHZ4z/c//tFHolW045HD\nj4wsWzjWWbgeee9b2eXn7RjbzcbqudHUPhRQD4LxlCxMs+5c0IvK2rRu9VhYq7WM1q20DkunbrFx\nduTqqH8aaPbpT8X9731/fPeYZ8bfHvPMuHPL8XHPE54S927aHHdtOT7uuq89/73VtlX0rbutKpZY\nTtK1nsBbHYmNiJ7A212bWhzqMjels9x9b1i0VWB1tDrvfM9CsOz8d8mR5+prs0yA7j5uV7vvxkT1\n857v65uARNuKAAAgAElEQVSe/rnrv9P55KnPjrxox7H3/108+wffjp/9u9vjuX93e/z0vXfGw394\naRz856+e/7mUzjHGI69UqS7L0c6IctOfphv1bK5ZJxhPWZZlkcVC5y2isuWT9cfAmC3X4Wtt3Ro/\nc8934mfu+U7Pv8/lrbj3yM3xnUs/Hntu/UlnT9nKdnXdVSNlLF6HWw1whwZOF9a9901n7o7uRiwf\neMcxutsE/X93deR50ahzLNyQqL42K+n/+lqe3SyLeNWLfzqeceD78bNve3M87cDfxePmHl30fe0T\ntvZ8Xm1/TQzJ1hVPX5ZlkeedPdyLMiKz3hhYJ8G4ZlqH1h+X0b3AT/uImCV5Nr4KnjTLah36uZNO\njse2be8poBQRsaFox+ZnPj1av/KS+Miv9IaZub7pwf1Tlnu+1g24K05nZj16Rp1j4cZrf4Duftz5\nWudnlwrU/V/rflwdgVhoSyfEUR88Kjbe+7eLjuuxn39hzJ108rLHPevr+xmf7h7uEfY4BtZPMK6Z\nLMui1VrorDQpw5hyVH8KUQ9mVjvqaz1HH/zYp5auSv3RTy75WLP2PKUk67sxUQ3U/V/rfhyxdFta\nS7tZSVPalNHiwTT5erpaJXZgNgjGNVYdPY7oFKexpSP91rpnbBa2aErdoB354rhnxIGv3hgbbr4p\nWt+9I9onbF1xxI90rNSWRtluumv967YvtoBUT63ciDFENPMGVB0IxjXWXTdTLSphGiyDqFaezRbv\nkmLEY43q2jlfi1G8xnMnnbymYNPkESHWbq3taa3tZi2/T5tiLTpF6xa2JlSlGhjEEl1l6qT/gl7d\nkgAGoWvApFx+3g4VMGfQNG+kXX7ejlrcxHMzcXjjvj70bh83tl8DzCDBuEFa+UK4aRe9e4LCUqqd\nAnfNh9PEsCfIMCp1eT2neR7W5TlgZdX3us6+7KPpK3XrqACzy1TqBumfWl3UdB2NIhWTsdobdFmW\nK+4Ny2yry/nX5Gno1JMp+wzCmmNSY+/i9ROMG6ZbPKm67Y79jlnKaqG4LsGpSZrQITfVk1Gqc3ua\n1PlY5+egqawbB+rIVOqGyvMs8kOvXqEoF0vQJtJT16meTZyGTnMC4TjbfVOeA5bWqvRyLT8DViMY\nN1jet46GdJiaOl11C3pN6LzXNbQzO0Z9XmqzzZdl2XwxrqIYzXrjsy+63mg3zChTqWdA/759ddie\noFukQqdiurLMTZNZ14RQTPM0tU1Vj3u9NxCdU7Pl0Aq0iDCTitln4GQ4gvEMyLIsWtWiXGVEHk6K\nVFXviOdZhOvjeNRhvXETO+/WFjIp/SF5JcLw5E3qGppVkrEbxcBKBOMZ0X/hV4QxXdX3/Wq76NIB\nHK1pVF32GjIus9q2ZvFvYjDVoqXDzqwzKw5mk2A8o9wVnU1rubPefe3z+iyBZYRmoSNWh9F2FpuF\ntgXLyfMsyqKzjWFRROS5jhLQS/GtGZRXqzBOeUGNIhWTVd272BZekzOJYlyzOpIHpG2SxQy7N4zL\nMIDAbNHfHg3BeAbl1SqMpS0KUlJ9qaddgC01towZXN2qe6dqVtsX9Ku+LyrEBaPTLsqYazd/Iadg\nPKPyyjzaovntlDVa7Y3e9iPjZ8sYgHobxX3j7jpjSFl38K0sIx58+NEpH83wBOMZZ9BwdijB3xyj\nCLOpjOIZNZ6uVNoZ9TbpG4BqcMDwirKcH3zLs4gnPeHx0z2gEVB8a8Z1qzB249Q09jhWvXGy7F1c\nH4MWmUo1pCjGBUxS/zaX0+gbwah0B04mfZO5OiM1n5G7TYLxjMuyLHJ7HCfF3sX1s1zYzbKILVs2\nxb33PuRmBhOX4k0Y6KpuZ9guIlqqVMPAZm0wxlTqBFTvgpblbDVgFnPXm6YyrRrSNs1rwHqLcakG\nTErKSojIstlbliAYJ0glxtm13AUq1Sm6wGKuB7CYQQNYWVkuzECN6PQ5Z20wxlTqxFSnPFhT0xwr\nFd7qv3sHTWa9MTBp3XosEfpGNMek3yf7d7mZxfPEiHFi8iyi24yLojdUjZNtDcbH3sXAWtn+CxbL\nhuwb6eMwy7rnQwqTKgTjxHSKcXU+LsO06lngNWQWWW8MaZrGzZv+vpFp1bCgej60Zjw5zvifx1L6\ni3HRXGu5q209IeA6ACur9o3ccIaOoiznz4c8n/2ZiYJxwmatktwsWq3apTdvZplRY2DS9I1ogpVq\nz4xSz17FMx6KIwTjpFXX1ERMZr2xbQ1Gy4g/sBqjxTTRtG6M9ff9B+kb6eMwixLIw/ME44RV19RE\nCFlA/SgWBUxSlmU96yj1jUhN/82glGZR2K4pcZ21Ap0ToCgjMu8AjVPdZgKgyk0FGJy+ESkr+opt\nVc+HWWfEmB7tYvXvYTJGsXex0TZmhfXGkJ66nPeD9I1s3cS4jXN9cVGUSW8DKhhDA3UvWlmW3kUL\nWJ11xTB6k6jFAtOU+gxEwZh5vWtqxndmuJs6nLJaOl8mJiFmQEB6pn3eW29MCtY6E3HWWWPMvCzL\nIss6UyjaRUQr9w4A0DRuHsDodAqVllEUh0bT1piOuwMAzkeaoDpSnGcRE9gJqpYEY3pUTwZ3Rqdj\ntdH0Yg1rP0yjZJZdft4Os06AicmzLIpuMS61WJiScb3vVdcVp1Zsq5+p1PSoBq3U1xnUlRsWwFLc\nEGOW1aEQV8pTTJldaxlwSYVgzJKq58W41huffdH1Rn3WKe3LFkx/3SGQnv66Hqv1j9RUoc6sK17M\nVGqW1D+l2glTL1lm5BhY4CYBjF+WZdHKy/ntm/SPaDLrihczYsySOsUmOh8XpS0K6sbexdBRh+mV\nwGTU4f2tf8lZoX/EBIx67+KiXGpdMUaMWVZPsQnX/YnoXvhW6ugv7F3sRYHUWVcM07VaMa7udOpU\nztNqgKv2ZVL5+5ui2m6F4gWCMWvihuh0VUfs7V0MvVSphrR0z/lRjqCtR3VZU1mWyQaM1W7q96+1\nFpLrwbK8xUylZlXjDGIKU6yN6xbQz2gxTFeeLSxtahfpLTtbTxFV/b76MNCymBFjVmXW7vRV32vt\nXQyLddt+Kh2uK3bvcKcfpizLssijnC9a1F5hWnX32jQL79PLTZce9DEiZuP5GKdRvqdV18Onvl/x\ncowYs6pOFcaFz1O7IzoJK911LSsFEgCA+hTeW+pm9az2k0a9zabR48kpy9K64jUQjFmTniqMCU4X\nmiZPNaxdXTrL47KhlcUX9pw27cMA+lQHEGZxN49xrekWjidjpdkMLBCMGVgZow9ro74LOUtUBIfB\n1GFLF2D86nSuV2fXleXS791C4PL0A6kDa4xZF2GtXqwvhtl3+Xk7VtzDHJiu6prNWRkw3rVncmE1\nta2tVrOWLTxX07+u2MjxyowYMzD9ssnTGYbBzcq0aje+YGV1O9dbfb3r/mnVTRkdfc3ua6a+JRbr\nV11X3KmgXp9zpK4EYwaWK8Q1Usut27F3MQA0T3/R0iauOd615/qYm9LwYlNuHNRdu6fY1vSOo0lM\npWZgWZZFlnUqJRdFRJ4362LfFNVnVUl9WJ/Lz9vR6A6W0WJYm7pt2dY/rbopS9Dq8vwxWkaL18aI\nMevSHcEsY3RxTVGKXmVl+stK6lR8BOqoqedIU48bWKzbV+qOHOvzrCz152dUVcD7p/WzMk8X69K/\nfROjVZbl/Juom3wwGnVbhwiMXh3P82o4add4y8txbcm0XqmH4/XoX4ZnpHgwgjFM0fLrixc+dlGD\ntJg+DbOlf81xd+1nXYKfNb2zoSxL64qHZI0xQ1H6fTyashYJmqZu6xD7CcQwvDrWFuiuOa5bxZBR\nbAk0bt3XMoXr4zDttr/vaGBlcEaMGUqnENfoHq8ud0+bwKgSrF/dpls6n2G06rpGP897R/KKooxd\ne/bq+zC0ms7QbxTBmKFVi0MVzsqRsUUTjFfdwjEw+7Is6+s3dWq1THrdcdOmTxs4WZ7tPUdHMGZo\n/YW46lpUok6We0OqPndmwMD41WFUyWgxjE8db4D1T3Eto7Ms7aw9eyfy++tWZIv1q64rzjJ9x2EJ\nxoycXLx+3ecuU0kQJmpanec6BHOYdXU+z1p5RPfKUxQRn/7q38Sjj7XH9vuaHoqbNtI9iPW8NtV1\nxapQD0/xLUZuFIWjUiq00FWW5fxzt5a9i4HRmlRhLucvTEddi3Ll+cKo31f//O746p/fHf/u9dvj\n+c98UkREbLj5pmh9945on7A15k46eV2/p25/N6NhF5PREowZqSwzYrxenjaoh3F1nk2ZhumrQ2X6\npa4DZ190fc+02Is/e2v8wk9virM+/q44+s+/Mf99j23bHg9+7FNRHPeMNf++JlSeZv3yzG4moyIY\nM1J5FtGdBVKWpbtXAyh79p7zvME0VTuuw0w9FIahnro3wMY5rXiQ87/7fUVZxo23/l18/Mt/Hf/t\nbx+Kv3jxb8Svbnx6vObP/3Mc/thPYuO+W+MJZ7wxDnz1xlUfc1ZHibuFuGbl2jro69RfjybPsyga\nPD2+TgRjRqpTbbEzJbhdRLRyJ2q/5e7ceqagnvpD8qA/A9TTes7tpYzyBlieZfHyFz49XvT3/yuu\n+fA18dXnnxJXvvifx3U/98vx2j/7QvzKf/9yPH7frbHh5ptWnFZtlHh29a8rZnQEY0Yuq+xeb1r1\nYFabim70CabL+Qezqf/cXm40eVLXgM3f+278m//yR3Hat66Nj/7iv4y/+OkXxkd2/Hr8x5NeG//k\nW1fHi7/9nThsmWDc9AJbazFro8ZrVRSldcVjJBgzclklGQ+z5iHFi151KjoAMB3T7nu0T9gaERHP\nvPdv493/6Xfjr57+vPjMya+LW396e/yHl/3L+Pj/KqL8/c72Tt1RYSPEs89a4vEaerumvXv3xi/8\nwi/ElVdeOYrjYYa4NA8mC3f+AICIuZNOjse2bZ///Pnfuy3e85/eFRdduTv+9x/89yizhS78XLuM\nItHE1PTtm9Y6ut+/rpjxGCoY7927N774xS/GiSeeOKrjYYbkdsnusdrFO1vl+TKNGgDS8eDHPtUT\njiMinnXMT8VZb/3V+IN//dKe9aXdXNxOOCTPqrIso6gUaLWueHyGmkp94oknximnnBJve9vbBv5Z\ndztmX5ZlkWXl0OuMm9JWusc5yPH23AEc4HcwG9bTZkibNsOgtJnmKp/xjHjgazcuuY/x0RHxx287\nJc56394oo1OfpCxj/uOudlHOfJDqLr27YvdsDh4U5UKB1lbeu2SxTmbhWjNUMN60adO6fu5JTzoi\nWq26Dyd2G1331S0r/x5LfG2t3zeOr9X38TtrZhfWG+fZUs/ryo+xa88NERHxhT2nRRNs3rzceVG9\nUnT+tqJnakze87XFH0ds2bK+c456W77NwNK0GQalzTTYK18REa9Y8kv/+aJ/EhERr9l9Tcy1i8jz\nLMpyYVCiLBdql7Tb/dtoNq9fudLXmtRHes3uaw59tHq/uLfYVt7ztTq9Nt1rTJOvNWsKxtddd11c\ndtlli/796quvXtcv/dGPHm7A3YSy77/9/77Sx5P+Wn0fv/o6F0UZWb7U87q2x7/33oeizrKsczG4\n776Hlhkl7/3Hsiz6vm/55797F7TuzwGDWb3NQC9thkFpM2m4YvfLY9eeQ8u1smx+3Wp1t4vOaPJC\nI5hrF/P9tLIsKqG5vv3Klb72mt1/GhHRkJHjQf7OzvTpzgz5+r429933UCOuNSvdQFlTMN65c2fs\n3LlzZAcUEbV+wiI6J1VnpLLmB9oww7zudW8zXd3pTFVL7Sc4yN/TlL+d9VmqzcBKtBkGpc2kqZUv\nVKzO88XtoPtxu4jo9nmLomzAANbyZqGd9xfbyvMsippvW1KdpdDU18B2TSv4wp7T4t57F9/1aHL1\nu2lTD2KB5wIAGJVugc7l+ql5lkVkMT+a3DoUlPv7I0UZPeNCRVHaamSElhosqSrL8tCNio5ZXyNe\nJ0Mt9P3oRz8ap59+enzjG9+Ij33sY3H66afHX/7lX47q2Grr8vN2xJ+87ZS4/Lwdcfl5O+wZt0bD\n3H3sFlZIjUrUAMAg1to3zbIs8kOpq5V3/tf5997vK8qYr4o81y6jXeM7+7PQX6xWoI6wleckDTVi\n/Gu/9mvxa7/2a6M6lsbqDy5NPyHHpVOIq/NxWfYXgEhHdXpMni++AAIATFJ1fXF1+nXnaxFRVqrv\nVHJxp6BX99/T7duNQrd/WN/bDrPPVOoxqAZlIXlBlmWRZ2UUZWctSytP49Tv37y9+ld7+wAARu3y\n83aMrA9aDcr90687Bb06H7eLiCyrR9+u+7fXbdZdf5+wqjoQ38qjZzo1kyEYj1l/SF7uZEhFdeu1\npi7MH1Z56EKXZyvvRVe3izkAkLYsyyLLYr4QVH9Br9T2UB6VoljYYqvOexXPurpvJjxTqmuTU1Wd\nYjPoEpVZWDdSluXC7m/eLACAMenWwhmnPMvmR5RbeW+hqLJcGPXsBL/JBr069RvPvuj6FY+l2ic2\nHX16BOMpSb1oV6p/uQsfADBJk+pzVot5dT5f+FpRCcllOfmQXFdFpWNodH36TKWeopSnWed57wVy\nFkPiUncG1/I+oBI1ANB0rTyb33ooy/r3TJ6M7qhxHftVRVHOD5jkDdmreNYZMa6J1EaQu2tUIjpV\nmQe5c7jadJS6m8F7AABAjU1iWvVSugMfrTyb3w6qX3sK06wnabnBr+oswtxwcS0YMa6R1KpZd7dv\n6ttHfuZVt60CAEhBt6BUK+/0+4r5mYPj386zrlWq82zwmjuMj2BcU90Td5YDcrXiXip7+WYqUQMA\nU9Ldxmmay/eyLIssIopY2Ce5jls+Dau/D1+WZe8osenTtWMqdc2lNsV6VlTfdKrTg8yUAQBY0D/N\nujqruhhhoa5pV6nuLB1c+Nz06foRjBtg1sPxcmtOVjLti9sg1lKJWsEtAGAS6tivXFiL3DuIUBS9\nxVqbYqkBkuYcfboE44aYVtGESagW4opo1oVvNZ0tCaZ9FAAA9de/5VNVu+jd3mi9JlnEtThUlbtr\nPYNBTI6Xp2HqeJdvFPo3hG+q1+y+JnbtWbjYruX6bbQYAJikpgy4LBpBrvSrhhlIGefMw+pocdEX\nimdxe9JZIhg3UFMuZoOoXiiKcnZGjWfkzwAAmLjqCHLel1raRWerp7rp78N2g71QXH+CcYPN6uhx\nxNqrVDdhT2PXQQCgjprUl8yzbNGxLlSzHrxI16hHjXft2Ru//r69i6ZOK7LVHLZrarhZ3dapfvf/\nBtNfiXq5SvyzNvIPADAJrbw7y7B3L+SiKAcalBhkj+MNN98Ure/eEe0TtsbcSSdHRMQjB+finP/n\nxiWXzxklbhYjxjOiSXf8VpP3lOxvTkTetef6mGsXUVYKLeSZiyIAUF9NXaKXZVm0utOs+9YhV0dt\n19KXXG30OL/7rjjqFb8YR//jfxRP+DdvjqP/8T+K9qteFf/xmr+It37o/50PxVmmwFaTeelmyKyE\n4ywWph939nxb+YJWt62bqoe7XCZWcAsAYDTyyl7I/V2vdhHR7ts6aTnLLdF7whlvjI37bo0HHr8p\n/ssLXhG//U/fHW966bnxpb+6Px45ODf/O1t5ltyAyCz1aU2lnjGzMLU6y7LIo4x22ZlS3aBB44hY\n277FAAB1cvl5Oxrff4woo9XKembvRSws0ev820JHrTrtuhua20VnvfAFu34hfnjgJ3HvX9wW+485\nOf7yX5wZd205fv5nD3/0kXjZX38jdpz1T+IpO168pmNs8vO7lFkJxF2C8YyalYtbxNq2PJq2amn+\niM7dyuUOe9YuIgAAdVLtR7byziDLUv3Jopz/tp4g3S4i3vbh/7bwDz9/akREbHrkwdh213+Pl/31\nN+NF//OWeNzco/HgaT8XB2Ntwbi/D9jUvvosjRJXCcYzbDZGj9c2YtydTl2XkzTPey+wEbN7EQEA\nZsMs9B37ZVkWWRZRtMvY0DeavFo/c/vWLfGUgwfiWR+9JJ73vdviGffdHXnf0Ef7hK3rPrZqv7AJ\nz/ms92MF4wR0R4/nliuNXGMrVXSuk/4q1NU7lQAATXLF7h2xa88NMdd/l38G9I4md+ZRzx0KzdWP\nqyHwqI/+Tmy8765Fj/XYz79wvjr1sOoakmc9DFcJxtRap+Lgwp29oihrtR9c98K1WsEto8UAAPW3\nVJ/twY99ar4AV9dj27bHgx/95FiOoQ4j9yn2XQXjRDR5zXH/euNshTkvg+xFNyplWa5YcCu1iwoA\n0Hxf2HNavGb3n077MCZquT5bcdwz4sBXb1xyH+NJH8+4+/Mp91sF44TU4e7TKNRpVk9/1UMAAJpj\nkCA4d9LJEwnEKxllAa8UR4VXIhgnqMmjx1VlWS65HdKkCnGdfdH1K4ZiFxsAoMmu2L0jyrL5gypL\nmZV+2iz8DXUhGCeqqaPHrUq156KMaE1xuXFRmT/d6qtC7SIFAADNkU/7AJiuy8/bMV+Frwk6xbg6\nH5dlbzidlLMvuj527dk7v654oQp1xIZWHlfsFooBgNnRtP7iamZltJjRMmJM4/QX41rOuApx9a8r\n7lbJ3tDK4gt7Tot7731opL8PAIDREIhZjmBMY9cc53lEMeHCV48cnFvYFD4Wdiq+/LwdS27TBAAw\nC5q6DK/LKDGrMZWaiOhc7Jp2scizLKpbGi81rbpbiGsUHptrx7/+wI0Lvz93kQUA0tLEadVN7Ocy\neUaM6dG00eM8z6JoL0yrLse05njXnr1RFAsjxK084ordp4zldwEAMDwDGAzCiDGLNPmuWll2/9sb\nkM++6Po1B/4NN98Uj/vMp2LDzTdFRMT9Dx2MtlAMABARzRg1FooZlBFjltW00ePqlkntIiLPe8Px\navsb53ffFU84442xcd+t8//2rV96TfzuC/+viLwVERHv/1cviSc94fHj+QMAABiaQMx6CMasqEmF\nFjpbOS1UjO4W5irLcn47pZVUQ/HDh/1UfPKlb4wvbj81yiyP7Vu3xK7TnheHP84pAwBQxz6iUWKG\nYSo1M6UbgKtFudpFRPvQ2uPlinFtuPmm2Ljv1jhw+BPjqpP+aez69X8f177wVbGx/Vicdf0fx1ue\n8WOhGACgT12mVTd5KSD1oKfPmnSnVc+1x1PcatTyPIusst9wdclxuyjjrPftjcgi3vr6F8bfP/JY\n7P/WD+Lbr/3d+KvjnhdzrY0REXHSd26Os274kzj2gR/Eg//zldE++Rem8JcAANTbNEePjRIzKoIx\na9a0Nced0eMyWvmhitXzhbkOFdIqI/Zc+ReHvntTxD/4+Ths7mC8/LYb4rRbvxjP/sG35x+rfcLW\nCR89AAAwKYIxA6njepLVZFkWrSzmR7vzrBOMyzIii4iTn/+UeOrRPxXPv+T3Yvs3/jQeP3ew5+cf\n+/kXxtxJJ0/+wAEAGmRSMwyNEjMO1hizLnVZT7IeeZ5FK89iQyuLP37bKXH2ac+PV//vz4yfufhd\n0Xrez/Z872PbtseDH/3klI4UAKBZxt1HFIoZFyPGrFv1otSU9ccrXUyL454RB756Y2y4+aZoffeO\naJ+w1UgxAMCAxjXDUCBmnARjRqLu648HuZDOnXSyQAwAMKT+QZT1MkrMJAjGjEzd1h+7gAIA1MOg\ngyjCMJMmGDNy0xw9dgEFAKgn/TTqTDBmLCa5/tgdRQAAYBiCMWM3qvUlSz0eAADAsARjJmq5ULtU\nYDYSDAAATIJgTC0IwAAAwLTk0z4AAAAAmCbBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABI\nmmAMAABA0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA\n0gRjAAAAkiYYAwAAkDTBGAAAgKQJxgAAACRNMAYAACBpgjEAAABJE4wBAABImmAMAABA0gRjAAAA\nkiYYAwAAkLQNw/zwnXfeGe94xzuiKIr48Y9/HG95y1viZS972aiODQAAAMZuqGD8rne9K17/+tfH\nq171qrj99tvjTW96U3z9618f1bEBAADA2A0VjC+55JI4/PDDIyJi8+bNceDAgSjLMrIsW/Vn1/At\nU9U9vrofJ/WhzTAobYZBaTMMSpthPbQbBjULbSYry7IcxQO9+93vjoiId77znat+b7tdRKtleTMA\nAADTt6YR4+uuuy4uu+yyRf9+9dVXR1mWccEFF8Rdd90Vl1566Zp+6Y9+9HDt7yZkWcTmzZvivvse\nitHcOmDWaTMMSpthUNoMg9JmWA/thkE1pc1s2bJp2a+tKRjv3Lkzdu7cuejfy7KM3/qt34osy+JD\nH/pQbNy4cc0HVecnrKosm3Os1IM2w6C0GQalzTAobYb10G4YVJPbzFBrjD/+8Y9HURSxZ8+eUR0P\nAAAATNRQwfgjH/lIHHPMMXH66afP/9v73//+eMpTnjL0gQEAAMAkDBWMv/GNb4zqOAAAAGAqlIYG\nAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIx\nAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROM\nAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpg\nDAAAQNIEYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIE\nYwAAAJImGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJIm\nGAMAAJA0wRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0\nwRgAAICkCcYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0wRgAAICk\nCcpjkzoAAAaSSURBVMYAAAAkTTAGAAAgaYIxAAAASROMAQAASJpgDAAAQNIEYwAAAJImGAMAAJA0\nwRgAAICkbRjmh/ft2xfve9/7Is/z+MlPfhJnnnlmnHrqqaM6NgAAABi7oYLxRz7ykTjvvPNi+/bt\ncdddd8WrX/3q2LlzZ2RZNqrjAwAAgLEaKhh/8IMfnP/47rvvjmOPPXbNobju2bl7fHU/TupDm2FQ\n2gyD0mYYlDbDemg3DGoW2kxWlmU5zAPccccdcf7558f9998fl156afzsz/7sqI4NAAAAxm5Nwfi6\n666Lyy67bNG/X3311fMf79u3L84999y45ppr4ogjjhjtUQIAAMCYrHvEuCiKuO666+LUU0+dnz59\n2mmnxQUXXBDbtm0b6UECAADAuKx7u6Y8z+PSSy+Nb37zmxERsX///ti/f38cf/zxIzs4AAAAGLeh\n1hjffvvt8Z73vCfyPI+HH344zjrrrHjlK185yuMDAACAsRq6+BYAAAA02bqnUgMAAMAsEIwBAABI\nmmAMAABA0jZM+wDqbN++ffG+970v8jyPn/zkJ3HmmWfGqaeeOu3DosbuvPPOeMc73hFFUcSPf/zj\neMtb3hIve9nLpn1Y1NzevXvj7W9/e/zbf/tv4w1veMO0D4ea+vCHPxxf+cpXotVqxbZt2+Ltb3/7\n/HaJsJwHH3ww3vnOd8af/dmfxX/9r/912odDA1x++eXx5S9/OVqtVhx//PFx4YUXxmGHHTbtw6Km\niqKIiy66KG655ZbYsGFDbN68Od773vfGkUceOe1DG5gR4xV85CMfifPOOy8+8YlPxB/8wR/E29/+\n9lCrjJW8613vite//vXxyU9+Mi688ML4nd/5nWkfEjW3d+/e+OIXvxgnnnjitA+FGtu3b19ce+21\n8clPfjKuvPLKuOOOO+IrX/nKtA+LBvjN3/zNOPnkk6d9GDTELbfcEtdcc0185jOfiauuuioOHjwY\nV1999bQPixr71re+Ffv374+rrroqPv3pT8fhhx8en/3sZ6d9WOsiGK/ggx/8YGzfvj0iIu6+++44\n9thj3Z1nRZdcckns3LkzIiI2b94cBw4ccDOFFZ144olx8cUXxxFHHDHtQ6HGvv71r8cpp5wSj3/8\n4yPP8zj11FPjhhtumPZh0QAf+MAHzFxizbZv3x5XXnllbNy4MSIijj766Lj//vunfFTU2Yte9KK4\n+OKLIyLi0Ucfjf3798exxx475aNaH8F4FXfccUe89rWvjd/+7d+ef9FhOUceeWS0Wq2IiLjsssvi\nta99rZsprGjTpk3TPgQaYP/+/XHMMcfMf37MMcfEPffcM8UjoilcYxhEq9WanwJ75513xg033BCv\nfOUrp3xUNMGePXvilFNOia1btza2zVhjHBHXXXddXHbZZYv+/eqrr46tW7fG5z//+di3b1/8xm/8\nRlxzzTVGdlixzZRlGRdccEHcddddcemll07h6KijldoMDMpMFGCcbr/99jjnnHPiwgsvjOOOO27a\nh0MD7N69O84999w4//zz40/+5E/izDPPnPYhDUwwjoidO3fOT3/tKooivvSlL8Wpp54aWZbFtm3b\n4ogjjojvfOc7sW3btikdKXWxVJuJ6HRWf+u3fiuyLIsPfehD81ORYLk2A2vx1Kc+Nfbv3z//+fe/\n//142tOeNsUjAmbVbbfdFueee25cdNFF80sKYTnf/va3o91ux3Of+9w47LDDYufOnfH5z3++kcHY\nVOpl5Hkel156aXzzm9+MiM40tv3798fxxx8/5SOjzj7+8Y9HURTx3ve+VygGRmbHjh3xta99LR55\n5JGYm5uLL33pS/GKV7xi2ocFzJjujhqXXHKJUMya3HHHHfHud7875ubmIqJTjGvr1q1TPqr1yUrz\nsZZ1++23x3ve857I8zwefvjhOOussxo7Z57JeNnLXhbHHHNMz3T797///fGUpzxlikdFnX30ox+N\nvXv3xne/+9048sgj48lPfnKcf/758YIXvGDah0bNfOxjH4trr7028jyPF7/4xfGWt7xl2odEzR04\ncCDOOeecOHjwYNx2223xwhe+MJ797GfHO97xjmkfGjV11VVXxcUXXxzPfvaz5//tJS95Sbz5zW+e\n4lFRZ2VZxsUXXxw33XRTtFqt2LJlS1xwwQXxxCc+cdqHNjDBGAAAgKSZSg0AAEDSBGMAAACSJhgD\nAACQNMEYAACApAnGAAAAJE0wBgAAIGmCMQAAAEkTjAEAAEja/w+B0ykaZKR6kgAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4cc4e75d68>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(15, 8), dpi= 80, facecolor='w', edgecolor='k')\n",
    "ax = plt.subplot(1, 1, 1)\n",
    "ax.errorbar(xpts, y_pred[0], yerr=sigmas[0], capsize=0)\n",
    "ax.plot(x, y, \"ro\")\n",
    "ax.set_ylim(-3.0,3.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "You can see how a sample of the function drawn out of the relevant function space (configured via the covariance function) emerges."
   ]
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}