{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gaussian Processes "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`conda install plotly`\n",
    "\n",
    "Consider a simple function\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sklearn.metrics\n",
    "\n",
    "np.random.seed(2)\n",
    "def f(x):\n",
    "    \"\"\"The function to predict.\"\"\"\n",
    "    return x * np.sin(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<ErrorbarContainer object of 3 artists>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XeclNX1x/HPoSlgiUZABSk2FBHbqrjrIioaQYqx62pI1BeJRmOJJtYYCzHFKIk1CCpRUFnssUPELkpRCRIr0oS46g8LIEXu74+zKGUXdneemTvzzPf9eu1r2GF8nrOye/bOufeeayEEREQkPRrFDkBERJKlxC4ikjJK7CIiKaPELiKSMkrsIiIpo8QuIpIySuwiIimjxC4ikjJK7CIiKdMkxk232GKL0LFjxxi3FhEpWJMmTfo0hNBqfa+Lktg7duzIxIkTY9xaRKRgmdnMurxOpRgRkZRRYhcRSRkldhGRlEkksZvZuWY2zcz+Y2b3mNmGSVxXRETqL+PEbmZtgV8BJSGErkBj4PhMrysiIg2TVCmmCdDczJoALYCPE7quiIjUU8aJPYQwF7gWmAXMA74IITyd6XVFRKRhMl7HbmabAQOATsACoNLMTgoh3L3G6wYBgwDat2+f6W1FpBDNmgXPP++PANtsA927ww47xI0rZZLYoNQLmBFCqAIwsweAUmC1xB5CGAoMBSgpKdFBqyLFYsUKuP9++Mtf4PXXa37NnnvCRRfBkUdCIy3Wy1QS/wdnAd3NrIWZGXAwMD2B64pIoXvnHTjgADj2WFiwwJP7m2/CokWweDFMmwZDhsDChXDMMdCzJ3z0UeyoC14SNfYJwBhgMjC1+ppDM72uiBS4Bx+Evfby5D18OEyfDuefD926QfPmsOGG0KULnH22v2bYMHjjDdhtNxg3Lnb0BS2R9zwhhMtDCDuFELqGEE4OISxJ4roiUqBuvNHLKrvsAlOnwimnQOPGtb++cWM49VQfzXfoAL17Q2Vl7uJNGRWzRCRZN90EZ50FAwbAc89B27Z1/287dfLJ1X33hRNPhCeeyF6cKabELiLJueceOPNMT+qjR3u5pb5+8AN47DHYdVc46iiYPDn5OFNOiV1EkvHaa/Czn0F5Odx3HzRr1vBrbbKJj9Y339yT++efJxdnEVBiF5HMVVXBEUfAVlv50sYNNsj8mm3awJgxMHeu/8IIWiVdV0rsIpKZEHzi87PP4KGHoNV6D/ipu+7d4Zpr4JFH4K67krtuyimxi0hmhg6FRx+FP/3Jlyom7ZxzoKzMl0XOnZv89VNIiV1EGu6DD+Dcc+GQQ+BXv8rOPRo3hjvugCVLPMnLeimxi0jDhABnnAFNmnjizWYrgB128JYDY8bAs89m7z4pocQuIg1TWQlPPw1XX12/teoNdf750LGjl2SWL8/+/QqYEruI1N+XX3pZZM89fdSeC82bw7XX+k7WO+7IzT0LlBK7iNTfNdfAvHlwyy1eismVI4/0lTJXXgnffJO7+xYYJXYRqZ/Zs70jY0UF7LNPbu9t5qWfOXPgtttye+8CosQuIvVz+eXeY/3qq+Pc/6CDvL3v4MHe/lfWosQuInU3dSqMGOFNvjp2jBODGVxxBfzvf3D77XFiyHNK7CJSd5dc4n1cLr44bhzl5bDffnDddVohUwMldhGpm8mTfYfp+ed7c66YzOCCC2DGDO9NI6tRYheRurn6am+pe+aZsSNx/fvDjjv6cXtqELYaJXYRWb+33vKj7s4+GzbdNHY0rnFj+PWvYdIkP9BDvqPELiLrN3gwbLyxJ/Z8cvLJXha66abYkeSVRBK7mf3AzMaY2X/NbLqZ7ZfEdUUkD7zzjrcPOOss2Gyz2NGsrnlz79X+0EPw8cexo8kbSY3Y/wY8GULYCdgNmJ7QdUUktiFD/DSkfButr/SLX/jKmGHDYkeSNzJO7Ga2CdADGA4QQlgaQliQ6XVFJA98+qmvWz/5ZGjdOnY0Ndt+ezj0UO8Lr6WPQDIj9m2BKuAOM5tiZsPMrGUC1xWR2P7xD1i8OP/7oJ9xhh/C8eijsSPJC0kk9ibAnsAtIYQ9gIXAhWu+yMwGmdlEM5tYVVWVwG1FJKuWLIEbb4TDDoNddokdzbodfjhsvbW6PlZLIrHPAeaEECZUfz4GT/SrCSEMDSGUhBBKWiV5JqKIZMe998L8+XDeebEjWb8mTeCkk+Dxx73VQJHLOLGHEOYDs82sc/VTBwNvZ3pdEYkoBLj+eujaFXr1ih1N3QwcCN9+C6NGxY4kuqRWxZwFjDSzt4DdgT8kdF0RieHll+HNN/0cU7PY0dRNly6w995w552xI4kukcQeQnijuszSLYRwRAjh/5K4rohEcuut3uzrxBNjR1I/Awf6Ltk33ogdSVTaeSoiq/v0Uxg9Gn7yE2hZYAvcjj/e19yPGBE7kqiU2EVkdXfeCUuX+safQvPDH0K/fjByZFGvaVdiF5HvrVjhZZjy8vxf4libE06AqioYPz52JNEosYvI98aOhQ8+gNNPjx1Jw/XpAxttBPfdFzuSaJTYReR7t94KrVrBkUfGjqThmjf3Xu0PPADLlsWOJgoldhFx8+fDI494t8QNNogdTWaOOw4+/xzGjYsdSRRK7CLi7r7bN/icckrsSDL3ox/5gSBFWo5RYhcR32l6xx1+QHTnzut/fb7bYAM44gg/9WnJktjR5JwSu4jA66/D2297GSYtjj0WvvgCnnkmdiQ5p8QuIj5ab97ca9Np0auX75596KHYkeScErtIsVu8GO65B446yhNhWjRr5u18H3nE5w6KiBK7SLF76CEvWaSpDLPSgAG+WemVV2JHklNK7CLF7o47oEMH6NkzdiTJ690bmjYtunKMErtIMZs1y3ebDhwIjVKYDjbZBA4+2BN7CLGjyZkU/kuKSJ3ddZcnvJ/+NHYk2XPEEd4mYdq02JHkjBK7SLEKwRN7jx7QqVPsaLKnf39/LKJyjBK7SLGaPBneeQcqKmJHkl1bbQXdu8PDD8eOJGeU2EWK1ahRPrF49NGxI8m+AQNg4kSYMyd2JDmhxC5SjL791teu9+4Nm28eO5rs69vXH598Mm4cOZJYYjezxmY2xcz+ldQ1RSRLxo+HefPSX4ZZaZddoF07ePzx2JHkRJIj9rOB6QleT0SyZdQoP4yiX7/YkeSGmR/AMXasH/uXcokkdjNrBxwODEvieiKSRd98A2PG+GEazZvHjiZ3+vSBr76Cl16KHUnWJTViHwL8BliR0PVEJFseewy+/LJ4yjArHXywTxYXQTkm48RuZn2BT0IIk9bzukFmNtHMJlZVVWV6WxFpqFGjoE0bOOig2JHk1kYb+Zr9J56IHUnWJTFiLwP6m9lHwL3AQWZ295ovCiEMDSGUhBBKWrVqlcBtRaTeFiyAf/0Ljj8emjSJHU3u9enjO1BnzowdSVZlnNhDCBeFENqFEDoCxwP/DiGclHFkIpK8++/3ycMTT4wdSRx9+vhjykftWscuUkxGjYLtt4e9944dSRydO0PHjkrs9RFCGB9C6JvkNUUkIfPn+/r1E07w5X/FaNVlj998EzuarNGIXaRYPPAArFiRruPvGuKww2DRInj55diRZI0Su0ixGD0aunTxXZjFrGdPnzgeOzZ2JFmjxC5SDObNg+efh2OPjR1JfBtv7N0en3kmdiRZo8QuUgzuv9/7rx9zTOxI8kOvXjBpEnz+eexIsiJVib1nz3Qe2yiSsdGjoWtXL8UIHHKI/6J79tnYkWRFqhK7iNRg7lx48UWVYVa1995ekklpOUaJXSTtVIZZW9OmcOCBqZ1AVWIXSbvRo2HXXWGnnWJHkl969fJDrmfMiB1J4pTYRdJszhxvU6syzNoOOcQfUzhqV2IXSbMxY/xRZZi1de4Mbdumss6uxC6SZpWVsNtunsRkdWY+ah83znfkpogSu0hazZ7t2+ZVhqldr16+ln3KlNiRJEqJvQBofb40iMow69erlz+mrM6uxC6SVqNHwx57wA47xI4kf7Vp45u2xo+PHUmilNhF0mjmTHj1VZVh6uLAA30D17JlsSNJTGoS+8iR/n383HPeR3/kyNgRiUSkMkzd9ewJX3/tvWNSIhWJfeRIGDQIlizxz2fO9M+V3KVojR4Ne+0F220XO5L816OHP6aoHJOKxH7JJd43f1WLFvnzIkXno4/gtddUhqmr1q29R70Se36ZNat+z4ukWmWlP6oMU3cpq7NnnNjNbBsze9bMppvZNDM7O4nA6qN9+/o9L5Jqo0d798JOnWJHUjh69oSFC2HixNiRJCKJEfty4NchhJ2B7sAvzSynTZ8HD4YWLVZ/rkULf15WpzXxKTdjhicnjdbr54AD/DEl5ZiME3sIYV4IYXL1n78CpgNtM71ufVRUwNChsMEG/nmHDv55RUUuoxDJAytXwxx9dNw4Cs0WW3gHzJQcvNEkyYuZWUdgD2BCkteti4oKuO02/3NKfumK1F9lJZSUqAzTED17wvDhsHQpNGsWO5qMJDZ5amYbAfcD54QQvqzh7weZ2UQzm1hVVZXUbVNP6/Olzj76CF5/XWWYhurZ05fTpaDOnkhiN7OmeFIfGUJ4oKbXhBCGhhBKQgglrVq1SuK2qefr84PW50vdaFNSZlauZ09BOSaJVTEGDAemhxCuyzykIvXFF/D443DRRdCvH3TuzCUnzWTRIlvtZYsWwSU/+xhOOAFuuAHefNOPPROprPRNSSrDNMwWW0C3bqmo5SZRYy8DTgammtkb1c9dHEJ4PIFrp9usWX4eZWWl11tC8LMYd9oJunVj1rs1r9ectWxLX3N7773+RKdOPll22mmw4445/AIkb8yc6ZuS/vjH2JEUtp49fbKuwOvsSayKeTGEYCGEbiGE3as/lNRrs3w5PPCAtwvt0AHOOw8WL4bLLoN//xsWLIC33oLKStp3sBov0b5DI++1PXMmDBvmhyhcf70/9unjR6FJcVEZJhkHHOA/jwVeZ0/FztOC8PXX8Oc/ezI/6ih491248kp/nDIFrrjCd7+tsiB/vevz27eHU0+FJ57wRP/733sjo/33hwED4O23c/blSWSVlbDnnrDttrEjKWz77++PL7wQN44MKbFn2xdfeCbu2BF++1vYeWd4+GHfSHLZZevslV2v9flbbgmXXw4ffuj3Gz/ej0S7/PLvu6NJwanThrJZs2DCBI3Wk9C6tb/zVWKXGi1b5pOb220Hl14K3bvDK6/4SS39+0PjxnW6TEWF/6cHHOCr2da76aplS7j4YvjgAzj+eH9XsNdePskq6aQyTLLKy72cWcDnoCqxJy0EeOQR6NoVfvUrHzW//jr861+eoXNliy3grrvgscf8TMfu3Rk5aLzWxKdRZaWflKQWvckoL/e5rv/8J3YkDZaqxD5+fOSVSu+9B4ce6vXtRo3g0Ud9hF5SEi+mPn1gyhRGbnspg27bW2vi02b2bF9RpdF6csrL/bGAyzGpSuzRLFkCV13lvSZefx1uvNFXtvTtC1bzypacatOGS76+mEW0XO1p9axPAZVhktexI7RtW9CJPdFeMUXp+efh5z+H//7Xa9rXX+8TmXlm1uyaf8GoZ32Bq6yE3XeH7bePHUl6mPmo/fnnvbSaD4OzetKIvaE++8yXGh5wAHzzje8aveeevEzqsI6e9W2X5zYQSc7s2T4hr9F68srL4eOPffVaAVJir68QfFJyp51gxAhfwjhtGvTuHTuydapxTTwLGfzthTBvXpygJDP33++PSuzJK/A6uxJ7fbz3HhxyCPzkJ/7Wd/Jk38K9ZsbMQzWuib90FhVf3goHHeTvQKSwVFb6qqt17IWQBtplF9hsMyX2VFt1cnTiRLjlFl/n2q1b7MjqZa018Vft7CWkGTN8onfhwtghSl3NmQMvv6zRerY0agRlZUrsqfX88z459bvf+TLG6dPhF7/wf/g06NHD5wZeew2OOy41h/mmnsow2Vde7i0//ve/2JHUW0qyUxbUNDl6332w1VY5DyXr6/N//GO4+WbfzHTGGWoDnCfWechKZaW/Y1Q3z+xZWWd/8cW4cTSAEvuaVqzwSdGddoJ//rNgJkcz9vOfeyuCYcN8Hb5E5YesUPOGsrlzvRSo0Xp27bUXNG9ekOUYJfZVvfoqlJbCT3/qE1IFNDmaiKuu8j42557LebuNW3/zKcmaSy7xDWSr+m5DmcowudGsGey7rxJ7wZo7F04+Gfbbz3fs3Hmnv/3addfYkeVWo0a+lLNzZy5/+1i2Wvxh7IiKVm0bx2bNwg9Y2XVX70Io2VVeDm+8AV99FTuSeinuxP7ZZ3DhhV6nrKz0UsS778LAgemZHK2vTTaBRx7BCFw17cd+6IDkXK0byrZe7puSTjghtwEVq/JyL8++8krsSOqlOLPXggW+yqVTJz/8YuVql8GDYaONYkcX33bbMXinu9l+4Vt+wpPkXK2HrOz/hH9y/PG5D6oYde/ug7wCm0AtrsQ+Y4YnqvbtvZ78ox/B1KkwapQOAF7DhB/24Z52F8Ctt8Lo0bHDKTq1HrLy7uVe99X3a25svLFvAnv55diR1Esiid3MDjOzd8zsfTO7MIlrJmb5cl+qeNRRvlv0hhugXz8/jq6y0neYSY2GdRrsI5bTToP3348dTtFZa0NZyTv+favRem6VlvrCiuWF01cp48RuZo2Bm4DeQBfgBDPrkul1M7JsmS/8PvdcaNcODj/cFwP/9rf+EzJypG86KkL1WRP/baOmPlHXuLFvXlq6NJuhyfrcd593GtRqmNwqK/Nd2VOnxo6kzpJo27sP8H4I4UMAM7sXGADk7iTlL7/0kczEib6+d+xYn8Vu1sxH5yef7OvQmzXLWUip0aED3H47HHmkH7j93UnaklMh+A7hHj28V7jkTmmpP770kp9UVQCSSOxtgdmrfD4H2DeB665t1ChP2osW+W/Q+fN9BP7pp9+/pmNHXzHQuzccfLDXyCQzP/4x/Oxnvqa/X7/cHvEn7q23vOf/2WfHjqT4tG/vv0xfegnOPDN2NHWSRGKvqQv9WnvSzWwQMAigfW1rudZn6lR4+mk/sLllSz9RvKTEk3m3br5TrHXrhl1b1m3IEPj3v72z5ZQp/v9fcmdlSezoo2NHUnzMfNReQBOoSST2OcA2q3zeDvh4zReFEIYCQwFKSkoa1ozkmmv8Q3Jvk01849aBB/pchdoO5E4IntgPOcQPKZfcKyvzxRZz5vi8XZ5LYlXM68AOZtbJzJoBxwOPJHBdiaTW5lM9e/qE9E03+Tun6qfUeiC7dv7qNS85ajVMPGVl/lggo/aME3sIYTlwJvAUMB0YHUKYlul1JY51Np8CnzzdeWdfAllg26wL1UGf3OsT/0ccETuU4rXbbr5D7KWXYkdSJ4msYw8hPB5C2DGEsF0IQcsmCtg6m0+Bd7sbPtzfkn73pGRLo/AtB1bdB336wKabxg6neDVtCvvsUzwjdkmXdTafWmm//Xx1wI030uWLwuqhUWh2XzCeLZbOUxkmH5SW+sKBAjhpTIldVlNr86k1nx88GLbZhgvePY2mK5ZkPa5idd0ed/nEdf/+sUORsjL49ls/bSzPKbHLamptPrVmgW3jjeHWW+m06G1OnPXHnMVXVBYt8t7rRx/tJTCJa+X+jQIoxyixy2pqbT5VUcOLe/fmmdYVnDRrsJ8yJcl66CH4+mvfOS3xbb45dOlSEBOoSuyylrWaT9WU1KvduN31LGyyqS+dWbEiZzEWhbvu8hpYjx6xI5GVSku9N3uef68rsUtGvmjWilu3/Yu/PR0xInY46TF/vu8VqKgo3kNf8lFZmZ/nMH167EjWSd8xkrGn2vwE9t8fLrjAT6WSzN17r48KVYbJLwWyUUmJXTIWrBHcfLOPZC6+OHY46XDXXd77aOedY0ciq9p+e2jVKu/r7Ers0mCrtR7otysjfzTCZ1pffTV2aIXt7bdh8mSN1vPRyoZgSuySRjW2Hnj2REb+4Jdw+ukFddpM3rnrLu/kqE1J+am01E8U++ST2JHUSoldGqTG1gOLjUua/AneeMNLM1J/K1b4b81DD4U2bWJHIzUpgDq7Ers0SK2tBz5r4YeEX3opzJuX26DSYNw4mD3b+95LftprL2/KpsQuaVN76wHzXu1Ll8J55+U2qDQYPhw220ydHPPZhht6cs/jOrsSuzTIOlsPbL+9H8Zx773w7LNR4itIn30GDz4IJ53kyUPyV2mpn7G8JD/7JCmxS43Gj/eP2qy39cCFF0KnTvDLX8KyZVmONiVGjvR3OqeeGjsSWZ+yMv+3mjQpdiQ1UmKXBltn64HmzeHvf/cdekOGxAqxcITgZZi99vJDHSS/lZb6Y56WY5TYJXv69oV+/eCKK/xgDqndpEnw1lsarReKNm1gu+3ydgJViV2y629/8x7Wv/517Ejy2/DhXlc/4YTYkUhdlZV5Yg8hdiRrUWKX7OrUCS66CEaPhrFjY0eTnxYtglGjvO/6D34QOxqpq9JS36T0wQexI1mLErtk329+A9tu68fpLV0aO5r8c//98OWXKsMUmpV19jwsx2SU2M3sL2b2XzN7y8weNDMNN2RtG27oE6nvvAPXXx87mvxz882w444+Cy2FY5dd/NjCPJxAzXTE/gzQNYTQDXgXuCjzkCSVDj8cBgyAK6/0nZXiJk/2pmlnnOENpqRwNGrkB7unbcQeQng6hLCy29OrQLvMQ5LUGjLEJ5q0I/V7t9ziO7sGDowdiTREWZkfC7lgQexIVpNkjf0U4Ina/tLMBpnZRDObWFVVleBtpWB07Ojdw8aM8dOBit2CBb4pqaJCk6aFqrTUByt51qp6vYndzMaa2X9q+BiwymsuAZYDI2u7TghhaAihJIRQ0qpVq2Sil8Jz/vnecuDMM/N2O3bOjBgBixd7GUYK0777ekkmz8oxTdb3ghBCr3X9vZkNBPoCB4eQhws6JavW1XagRhtsADfcAL17w1//Wq8Tl3r2bOA981EIPmm6336w++6xo5GG2mgj3ymcZxOoma6KOQz4LdA/hLBofa8XAeCww+DII+Hqq/2EjmI0bhy8+65G62lQWgoTJuTV4TKZ1thvBDYGnjGzN8zs1gRikmKwctnjuefGjSOWIUOgdWs45pjYkUimyspg4UJvCZEnMl0Vs30IYZsQwu7VH79IKjBJufbt4bLLvE3tE7XOuafT9Onw2GPe+XJle0wpXHm4UUk7TyWe887zjTlnnQXffBM7mty5/nrftHX66bEjkSS0bw9t2yqxiwA+Wr3xRu+1ce21saPJjU8+gX/+09eta3VYOpj5qD2PJlCV2CWuQw7xOvPgwd7UvZ569vx+tUxBuPlmX+Z5zjmxI5EklZX5QcB50p5aiV3iu+46aNw4/clu8WJP7H37wk47xY5GkpRndXYldomvXTv43e/g4Yd9UjGthg+HqirfpCXpsvvufmqYErvIKs45x0exZ53lI9u0WbIE/vQn2H9/6NEjdjSStKZNYZ998qbOrsQu+aFZM59InTHDE2CBWW+tf8QIr79edpm6OKZVaSlMmeJr2iNTYpf8cfDBcNxxcM01vtZ7FSNHep+l557zXmIja+1KlIeWLfOvaZ99fLJY0qmszI+BnDgxdiRK7JJn/vY32HhjOOUU/yHBk/igQd/3DJs50z8vmOR+992+4kej9XTbbz9/zINyjBK75Jc2bfy0pVdf9W33eKffRWt0Ilq0yJ/Pe8uWwR/+AHvs4YeNSHptvjnsvHNeTKAqsUv+OeEE6N8fLr0U3nuPWbNqflltz+eV4cPh/ffh97/XaL0YlJZ6Yl+xImoYSuySf8z8ZKENN4RTT6X9NjV3g27fPsdx1dfChXDFFb4Spl+/2NFILpSVwf/9n5/vG5ESu+Snrbf2niovvMDgA56iRYvV/7pFC9+smteGDIH5832Vj0brxWHlRqXIdXYldslfAwdC795UjDmSob+f+10jxA4dYOhQ/3PerpT59FNP6Ecc8f0Pu6TfjjvCD38Yvc6uxC75y8xr1C1bUjGqL+X7LOGAA75vKZPXK2Uuu8xnePP+bYUkamVDMCV2kXXYaiu44w544w1Om/H9Mpi8XikzaRL84x9+rmuXLrGjkVwrLfUa+6efRgtBiV3yX9++cMYZHDfnr5R8/jRQ+4qY6CtlVqzwAzRatfKJUyk+ZWX++Mor0UJQYpfCcO21zGjRhYveGQj/+1+tK2JirJRZbVds60WMnLAt/PnPsOmmuQ9G4isp8d4xESdQldilMDRvzpU730vL5V/Acccx+MrlebFSZq1dsZ9txKBGwxnZ+OTcBiL5o3lz2HPPqHX2RBK7mZ1vZsHMtkjieiI1mbHRrly7423w3HNUTLmAoUNZa6VMRUVuY6qx1r+iOZdcqjFTUSsthddfh6VLo9w+4+8+M9sGOASIXd2UIjC2TQWcfTYMGUKFjaJ7d75bKZPrpA55XOuXuMrK/BzfKVOi3D6JYcX1wG+AmrcHiiTtL3/xnuannUbnr+J20sunWr/kkZUNwSKVYzJK7GbWH5gbQnizDq8dZGYTzWxiVVVVJreVYte0KYweDW3acM3Uw9lq8YfRQhl8xXJaNFr9YJCC2BUr2bX11r5rLtIE6noTu5mNNbP/1PAxALgE+F1dbhRCGBpCKAkhlLTS6eySqTZt4MknaRyW8+eph/mRcxFUTP41Q1ecyqZN/HCFWLV+yUNlZZ7YQ+6LGetN7CGEXiGErmt+AB8CnYA3zewjoB0w2cy2zG7IUqzGj/eP73TuzMVdH6X1ktm+1v2LL3Ib0LBh8Pe/U/GrVuxe1jJqrV/yUGmp9wpauVU6hxpcigkhTA0htA4hdAwhdATmAHuGEOYnFp3IeoycUcpPN7wPJk+GQw+FBQtyc+OnnoJf/AJ+9CO49trc3FMKy8qNShHq7FqTJQXviab94f77fQVCr17w2WfZveGECXDMMdC1K1RWes1fZE1du/ppYBHq7Ikl9uqRe7zmCFLc+veHBx+EqVP9LfB772XnPi+/7OeWtm4Njz3mP7giNWncGLp314hdJCOHHw7jxsHnn/sP1NixyV7/4Ye93LPVVt4/oG3bZK8v6VNa6oONL7/M6W2V2CVd9t/fG7dsuaUn4d/+NvO4/dEzAAAGfklEQVTdf8uXw1VXeW/1XXZRUpe6KyvzxnATJuT0tkrsUtB2390/VrPddr6de9Agb8a1554+km+It9/2Xxa/+x2cdJIvy9lSC7+kjvbdFxo1ynmdXYld0qlFC7j1Vnj0UW/m0quX18afeqpuBw1Pnw6nnAK77ur1+nvugX/+0xs8idTVJptAt27w4os5vW2TnN5NJNf69vWkfsMNfobqYYd5jbx/f9h7bz/KbJNNYNkymDPHl00++aSP+DfYwPvSXHSR91cXaYjycj8JbNmynK2gUmKX9NtwQ7jgAk/SDzwAY8Z4v91//GPt1zZqBHvtBX/9q+80atMm9/FKupSX+8Bi8mQvzeSAErsUj2bN4Pjj/WPFCvjwQ//4+mto0sT7e6wcwdfTajtiRVZVXu6PL7ygxC6SVY0awfbb+4dINm25pX+fvfACnH9+Tm6pyVMpWKsdSdfRPxfJSz16+ARqXSbuE6DELgVprSPpZvrnSu6Sl8rLfePc22/n5HZK7FKQajySbpE/L5J3Vq2z54ASuxQkHUknBWXbbX2ZrRK7SO10JJ0UFDMftT//fE4O3lBil4I0eLBvLl2VjqSTvFZeDnPn5uTgDSV2KUgVFX4E3QYb+Oc6kk7yXo8evl9i2rSs30rr2KVgVVTAbbf5n7VBSPJe165+wlfLllm/lUbsIiK50KhRTpI6KLGLiKSOEruISMpknNjN7Cwze8fMppnZn5MISkREGi6jyVMzOxAYAHQLISwxs9bJhCUiIg2V6Yj9dOCPIYQlACGETzIPSUREMpFpYt8RKDezCWb2nJntXdsLzWyQmU00s4lVVVUZ3lZERGqz3lKMmY0Fajq995Lq/34zoDuwNzDazLYNYe09syGEocBQgJKSkuzvqRURKVLrTewhhF61/Z2ZnQ48UJ3IXzOzFcAWgIbkIiKRZLrz9CHgIGC8me0INAM+zTgqkTrSjlORtWWa2G8Hbjez/wBLgYE1lWFERCR3MkrsIYSlwEkJxSIiIgnQzlMRkZRRYhcRSRkldhGRlFFiFxFJGSV2EZGUUWIXEUkZJXYRkZSxGPuJzKwKmNnA/3wLim93q77m4qCvuThk8jV3CCG0Wt+LoiT2TJjZxBBCSew4cklfc3HQ11wccvE1qxQjIpIySuwiIilTiIl9aOwAItDXXBz0NReHrH/NBVdjFxGRdSvEEbuIiKxDQSV2MzvMzN4xs/fN7MLY8WSbmW1jZs+a2XQzm2ZmZ8eOKRfMrLGZTTGzf8WOJRfM7AdmNsbM/lv9b71f7JiyzczOrf6e/o+Z3WNmG8aOKWlmdruZfVJ9XsXK5zY3s2fM7L3qx82yce+CSexm1hi4CegNdAFOMLMucaPKuuXAr0MIO+Pnyv6yCL5mgLOB6bGDyKG/AU+GEHYCdiPlX7uZtQV+BZSEELoCjYHj40aVFXcCh63x3IXAuBDCDsC46s8TVzCJHdgHeD+E8GH1AR/3AgMix5RVIYR5IYTJ1X/+Cv+Bbxs3quwys3bA4cCw2LHkgpltAvQAhoMfXhNCWBA3qpxoAjQ3syZAC+DjyPEkLoTwPPD5Gk8PAEZU/3kEcEQ27l1Iib0tMHuVz+eQ8iS3KjPrCOwBTIgbSdYNAX4DrIgdSI5six/+fkd1+WmYmbWMHVQ2hRDmAtcCs4B5wBchhKfjRpUzbUII88AHbkDrbNykkBK71fBcUSzpMbONgPuBc0IIX8aOJ1vMrC/wSQhhUuxYcqgJsCdwSwhhD2AhWXp7ni+q68oDgE7A1kBLM9MRmwkqpMQ+B9hmlc/bkcK3b2sys6Z4Uh8ZQnggdjxZVgb0N7OP8FLbQWZ2d9yQsm4OMCeEsPKd2Bg80adZL2BGCKEqhLAMeAAojRxTrvzPzLYCqH78JBs3KaTE/jqwg5l1MrNm+GTLI5FjyiozM7z2Oj2EcF3seLIthHBRCKFdCKEj/u/77xBCqkdyIYT5wGwz61z91MHA2xFDyoVZQHcza1H9PX4wKZ8wXsUjwMDqPw8EHs7GTZpk46LZEEJYbmZnAk/hs+i3hxCmRQ4r28qAk4GpZvZG9XMXhxAejxiTJO8sYGT1gOVD4GeR48mqEMIEMxsDTMZXfk0hhTtQzeweoCewhZnNAS4H/giMNrNT8V9wx2Tl3tp5KiKSLoVUihERkTpQYhcRSRkldhGRlFFiFxFJGSV2EZGUUWIXEUkZJXYRkZRRYhcRSZn/B6zQ6ylNUeYHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xmax = 10\n",
    "X_space = np.linspace(0.001, xmax - 0.001, 1000)\n",
    "\n",
    "# Get the \"real\" value\n",
    "y_real = f(X_space)\n",
    "\n",
    "n = 10\n",
    "# Uniform sampling and add Gaussian noise\n",
    "x_sample = np.random.random(n)*10\n",
    "\n",
    "sigma = 1.0\n",
    "noise = np.random.normal(0, sigma, n)\n",
    "y_meas = f(x_sample) + noise\n",
    "\n",
    "plt.plot(X_space,y_real,\"r\")\n",
    "plt.errorbar(x_sample,y_meas,yerr=sigma,fmt=\"bo\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probabilistic interpretation\n",
    "Let our samples be $x_i, y_i$, $i \\in [0,n]$, then the probability of value $\\hat{y}_i$ at position $x_i$ given the noisy $y_i$ is:\n",
    "$$ Pr(Y_i = \\hat{y}_i | y_i) = \\mathcal{N}(\\hat{y}_i; y_i,\\sigma) = \\frac{1}{\\sqrt{2\\pi}\\sigma}exp\\left(\\frac{-(y_i - \\hat{y}_i)^2}{2\\sigma^2}\\right) $$\n",
    "\n",
    "If we assume that every sample is **independent** of the other (which are not!), then the probability becomes trivialy a **multivariate Gaussian distribution** with diagonal covariance\n",
    "$$ \\mathbf{Y} \\sim \\prod_i Pr(Y_i = \\hat{y}_i | y_i, \\sigma) = \\mathcal{N}(\\mathbf{\\hat{y}}; \\mathbf{y},\\sigma I)$$\n",
    "\n",
    "Following this idea, we model the $\\mathbf{Y}$ random vector as a Multivariate Gaussian Distribution\n",
    "\n",
    "$$\\mathcal{N}\\left( \\mathbf{\\hat{y}}; \\boldsymbol{\\mu},\\Sigma \\right) = \\frac{1}{\\sqrt{(2\\pi )^{k}|\\Sigma}|} exp\\left(-\\frac{1}{2}(\\boldsymbol{\\mu} - \\mathbf{\\hat{y}})^{\\mathrm{T}}\\Sigma^{-1}(\\boldsymbol{\\mu} - \\mathbf{\\hat{y}})\\right)$$\n",
    "\n",
    "where $\\Sigma_{ij}$ is a value that indicates how the sample $i$ varies w.r.t. $j$. \n",
    "\n",
    "If we interpret this model as an stochastic **process** (i.e., $x$ represent time) one would expect that samples that are closer to each other are more correlated than those that are far away, and that only depends on $x$, not $y$. One can define then a function $\\Sigma_{ij} = K(x_i,x_j)$, called **Kernel**, that defines this covariance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian Processes\n",
    "\n",
    "In the real world we just have the data and we want to predict a new value.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD8CAYAAACSCdTiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEmpJREFUeJzt3X9sXWd9x/HPx21C5tC1aRoWiOs466oG2nlQXVBopf0IMHUE0TUDqdUdTbZKXkW7wtSoBYw2pCkFjQxWafyQoT+QeilCkEAFrEBYJTTYEE7p8oOkS8ti11BWk2QpaqhI6Xd/HGd2ghPbOU/uOffx+yVZ1+ex73O+vUo+ffKc5zzHESEAQJ66qi4AAHD2EPIAkDFCHgAyRsgDQMYIeQDIGCEPABkj5AEgY4Q8AGSMkAeAjJ1bxUkvuuii6Ovrq+LUANCxduzY8bOIWDaX91QS8n19fRoeHq7i1ADQsWyPzPU9TNcAQMYIeQDIGCEPABmrZE4eANrp2LFjGhsb0/PPP191KbOyaNEi9fT0aMGCBaX7IuQBZG9sbEznnXee+vr6ZLvqck4rInTw4EGNjY1p1apVpfur3XRNqyX19UldXcVrq1V1RQA63fPPP6+lS5fWPuAlybaWLl2a7F8dtRrJt1rSwIB09GhxPDJSHEtSs1ldXQA6XycE/HEpa63VSH5wcDLgjzt6tGgHAMxdrUJ+dHRu7QCA06tVyPf2zq0dAM6GnK4N1irkN2+WurtPbOvuLtoBoB2OXxscGZEiJq8Nlg36t73tbVq4cKHGx8clSbfddptsa9++fQmqPrVahXyzKQ0NSStXSnbxOjTERVcA7XO2rg2+853v1LFjx/TAAw8oIvSlL31Jr33ta7V69epyHc+gVqtrpCLQCXUAVTlb1wbXrl2r1atX67777tNVV12lp556SnfccUe5TmehViN5AKja2bw2ePPNN2vXrl16//vfrwULFuiGG24o3+kMSoe87YttP2J7r+09tt+VojAAqMLZvDa4YcMGdXd3a/v27Vq3bp2WLl1avtMZpBjJvyDp9oh4paQ1km6x/aoE/QJA253Na4MXXHCBrr/+eknSjTfeWL7DWSg9Jx8RT0t6euL7n9veK2mFpB+W7RsAqnC2rg0+8sgj2r9/v5YvX65169alP8E0kl54td0n6TWSvpeyXwDIwdq1a7Vs2TJ96lOf0sKFC9tyzmQhb/ulkr4o6d0R8ew0Px+QNCBJvdzdBGAeioi2nzPJ6hrbC1QEfCsitk73OxExFBGNiGgsWzan59ACAM5QitU1lnSPpL0R8ZHyJc0fOd06DaCeUozkr5b0DklrbT828fXmBP1m7WzdOg0AU6VYXfNvkjpno+aaON2t09zxCyAV7nitCNsqA2gHQr4ibKsMoB0I+YqwrTKAdiDkK8K2ysD8tGvXLi1fvly7d+9uy/kI+Qo1m9KBA9KLLxavBDyQv7vuukvf/e53ddddd7XlfLXbTx4AauH++4vXjRuTdvvggw9Kkj772c8m7fdUGMkDQMYIeQA42c6d0rZt0j33SB/4QHFc0q5du3T11Vf///Gjjz6qtWvXlu53JoQ8AEy1c6e0ZUtxd+KFF0qHDxfHJYP+8ssv15NPPqlf/epXkqTbb79dW7ZsSVHxaTEnDwBTbd0qLVkiPTuxme6SJZPt/f1n3G1XV5cuv/xy7dmzR/v371dvb6+uvPLKBAWfHiEPAFONjko9PSe2nX9+ktvR16xZo+985zv6+Mc/rocffrh0f7NByAPAVL29xRTNVEeOJLkdfc2aNdq4caNuueUWrVixonR/s8GcfGbYvhgoaf36IuSfe67YIvbw4eJr/frSXa9evVoveclLdOeddyYodHYI+YywfTGQQH+/tGlTsc/IoUPFnPymTaXm44+7++679cEPflCLFy9OUOjsEPIZOd32xQDmoL9fuu466aabiiWUJQP+ySef1OrVq/WLX/xCGzZsSFPjLDEnnxG2LwYSSnin6yWXXKJ9+/Yl628uGMlnhO2LAZws1YO877X9jO32bKuGabF9MYCTpRrJ3y/pmkR94QyxfTGAkyWZk4+Ib9vuS9EXymk2CXUAk5iTBzAvRETVJcxaylrbFvK2B2wP2x4eHx9v12mB2bv//sk9xJGVRYsW6eDBgx0R9BGhgwcPatGiRUn6a9sSyogYkjQkSY1Go/6fNIBs9PT0aGxsTJ0ywFy0aJF6Tt4/5wyxTh5A9hYsWKBVq1ZVXUYlUi2hfFDSv0u6zPaY7ZtS9AsAKCfV6pobUvQDAEiL1TUAkDFCHgAyRsgDHYDnBOBMsboGqLnjzwk4vo308ecESNzdjJkxkgdqjucEoAxCHqg5nhOAMgh5oOZ4TgDKIOSBmuM5ASiDkAdqjucEoAxW1wAdgOcE4EwxkgeAjBHyAJAxQh4AMkbIA0DGCHkAyBghDwAZI+QBIGOpHv93je3HbT9h+z0p+gQAlFf6Zijb50j6mKQ3SRqT9H3bD0XED8v2DbTNzp3Stm3SoUPSgQPS+vVSf3/VVQGlpRjJv07SExHxo4j4paTPSbo2Qb9Ae+zcKW3ZUuzfe+GF0uHDxfHOnVVXBpSWIuRXSHpqyvHYRBvQGbZulZYskRYvLjaHWbKk+Nq6terKgNJShLynaYtf+yV7wPaw7eHx8fEEpwUSGR2Vzj//xLbzz2fDdmQhRciPSbp4ynGPpJ+c/EsRMRQRjYhoLFu2LMFpgUR6e6UjR05sO3KEDduRhRQh/31Jl9peZXuhpOslPZSgX6A91q8v5uGfe06KKL4/fLhoBzpc6ZCPiBck3Srp65L2Svp8ROwp2y/QNv390qZNxZM4Dh0q5uM3bWJ1DbKQZD/5iPiapK+l6AuoRH+/dN11xfcbN1ZaCpASd7zOY62W1NcndXUVr61W1RUBSI0nQ81TrZY0MFAsDZekkZHiWOIJREBOGMnPU4ODkwF/3NGjRTuAfBDy89SploCzNBzICyE/T51qCThLw4G8EPLz1ObNxYrBqbq7i3YA+SDk56lmUxoaklauLLZrWbmyOOaiK5AXVtfMY80moQ7kjpE8AGSMkAeAjBHyAJAxQh4AMkbIA0DGCHkAyBghDwAZI+QBIGOEPABkrFTI23677T22X7TdSFUUACCNsiP53ZLWS/p2gloAAImV2rsmIvZKku001QAAkmJOHgAyNuNI3vZ2Scun+dFgRHx5tieyPSBpQJJ6eTIF6mjjxqorAJKbcSQfEW+MiCum+Zp1wE/0MxQRjYhoLFu27MwrRiVaLamvT+rqKl5braorAjAb7CePGbVa0sDA5IO/R0aKY4n96IG6K7uE8jrbY5JeL+mrtr+epizUyeDgZMAfd/Ro0Q6g3squrtkmaVuiWlBTo6NzawdQH6yuwYxOdZ2c6+dA/RHymNHmzVJ394lt3d1FO4B6I+Qxo2ZTGhqSVq6U7OJ1aIiLrkAnYHUNZqXZJNSBTsRIHgAyRsgDQMYIeQDIGCEPABkj5AEgY4Q8AGSMkAeAjBHyAJAxQh4AMkbIA0DGCHkAyBghDwAZI+QBIGOEPABkrOwzXj9se5/tnba32b4gVWEAgPLKjuS/KemKiOiX9F+S3lu+JABAKqVCPiK+EREvTBz+h6Se8iUBAFJJOSf/l5L+5VQ/tD1ge9j28Pj4eMLTAgBOZcbH/9neLmn5ND8ajIgvT/zOoKQXJLVO1U9EDEkakqRGoxFnVC0AYE5mDPmIeOPpfm57g6S3SHpDRBDeAFAjpR7kbfsaSXdK+oOIOJqmJABAKmXn5P9Z0nmSvmn7MdufTFATACCRUiP5iPidVIUAANLjjlcAyBghDwAZI+QBIGOEPABkjJAHgIwR8gCQMUIeADJGyANAxgh5AMgYIQ8AGSPkASBjhDwAZIyQB4CMEfIAkDFCHgAyRsgDQMYIeQDIWKmQt/33tndOPPrvG7ZfkaowAEB5ZUfyH46I/oh4taSvSPrbBDUBABIpFfIR8eyUw8WSolw5AICUSj3IW5Jsb5Z0o6Qjkv6odEUAgGRmHMnb3m579zRf10pSRAxGxMWSWpJuPU0/A7aHbQ+Pj4+n+y9A1lotqa9P6uoqXlutqisCOosj0syw2F4p6asRccVMv9toNGJ4eDjJeZGvVksaGJCOHp1s6+6WhoakZrO6uoCq2N4REY25vKfs6ppLpxy+VdK+Mv0BUw0OnhjwUnE8OFhNPUAnKjsn/yHbl0l6UdKIpJvLlwQURkfn1g7g15UK+Yj4s1SFACfr7ZVGRqZvBzA73PGK2tq8uZiDn6q7u2gHMDuEPGqr2Swusq5cKdnFKxddgbkpvU4eOJuaTUIdKIORPABkjJAHgIwR8gCQMUIeHYetDoDZ48IrOsrJWx2MjBTHEhdogekwkkdHYasDYG4IeXQUtjoA5oaQR0c51ZYGbHUATI+QR0dhqwNgbgh5dBS2OgDmhtU16DhsdQDMHiN5IBHW76OOGMkDCbB+H3XFSB5IgPX7qKskIW97k+2wfVGK/oBOw/p91FXpkLd9saQ3SeKPM+Yt1u+jrlKM5D8q6Q5JkaAvoCOxfh91VSrkbb9V0o8j4j8T1QN0JNbvo65mXF1je7uk5dP8aFDS+yT98WxOZHtA0oAk9fJvWGSI9fuoI0ec2SyL7d+V9C1Jx9cU9Ej6iaTXRcRPT/feRqMRw8PDZ3ReAJivbO+IiMZc3nPG6+QjYpekl005+QFJjYj42Zn2CQBIi3XyAJCxZHe8RkRfqr4AAGkwkgeAjBHyAJAxQh4AMkbIA0DGCHkAyBghDwAZI+QBIGOEPABkjJAHgIwR8gCQMUIeADJGyANAxgh5AMgYIQ8AGSPkASBjhDwAZIyQB4CMlQp52x+w/WPbj018vTlVYQCA8lI8/u+jEbElQT8AgMSYrgGAjKUI+Vtt77R9r+0lCfoDACQyY8jb3m579zRf10r6hKRLJL1a0tOS/vE0/QzYHrY9PD4+nuw/AABwao6INB3ZfZK+EhFXzPS7jUYjhoeHk5wXAOYL2zsiojGX95RdXfPyKYfXSdpdpj9gNlotqa9P6uoqXlutqisC6qvs6pp/sP1qSSHpgKS/Kl0RcBqtljQwIB09WhyPjBTHktRsVlcXUFfJpmvmgukanKm+viLYT7ZypXTgQLurAdqr7dM1QLuNjs6tHZjvCHl0lN7eubUD8x0hj46yebPU3X1iW3d30Q7g1xHy6CjNpjQ0VMzB28Xr0BAXXYFTSbF3DdBWzSahDswWI3kAyBghDwAZI+QBIGOEPABkjJAHgIxVsq2B7Z9LerztJ66niyT9rOoiaoLPYhKfxSQ+i0mXRcR5c3lDVUsoH5/r/gu5sj3MZ1Hgs5jEZzGJz2KS7Tlv+sV0DQBkjJAHgIxVFfJDFZ23jvgsJvFZTOKzmMRnMWnOn0UlF14BAO3BdA0AZKytIW/7GtuP237C9nvaee46sX2x7Uds77W9x/a7qq6parbPsf0D21+pupYq2b7A9hds75v48/H6qmuqiu2/mfj7sdv2g7YXVV1TO9m+1/YztndPabvQ9jdt7594XTJTP20LedvnSPqYpD+R9CpJN9h+VbvOXzMvSLo9Il4paY2kW+bxZ3HcuyTtrbqIGrhb0sMRsVrS72mefia2V0i6TVIjIq6QdI6k66utqu3ul3TNSW3vkfStiLhU0rcmjk+rnSP510l6IiJ+FBG/lPQ5Sde28fy1ERFPR8SjE9//XMVf5BXVVlUd2z2S1kn6dNW1VMn2b0r6fUn3SFJE/DIi/rfaqip1rqTfsH2upG5JP6m4nraKiG9LOnRS87WSPjPx/Wck/elM/bQz5FdIemrK8ZjmcbAdZ7tP0mskfa/aSir1T5LukPRi1YVU7LcljUu6b2Lq6tO2F1ddVBUi4seStkgalfS0pCMR8Y1qq6qF34qIp6VisCjpZTO9oZ0h72na5vXSHtsvlfRFSe+OiGerrqcKtt8i6ZmI2FF1LTVwrqQrJX0iIl4j6TnN4p/jOZqYa75W0ipJr5C02PafV1tVZ2pnyI9JunjKcY/m2T+/prK9QEXAtyJia9X1VOhqSW+1fUDFFN5a2w9UW1JlxiSNRcTxf9V9QUXoz0dvlPTfETEeEcckbZV0VcU11cH/2H65JE28PjPTG9oZ8t+XdKntVbYXqriI8lAbz18btq1i3nVvRHyk6nqqFBHvjYieiOhT8WfiXyNiXo7YIuKnkp6yfdlE0xsk/bDCkqo0KmmN7e6Jvy9v0Dy9CH2ShyRtmPh+g6Qvz/SGtm1QFhEv2L5V0tdVXCm/NyL2tOv8NXO1pHdI2mX7sYm290XE1yqsCfXw15JaEwOhH0n6i4rrqUREfM/2FyQ9qmI12g80z+58tf2gpD+UdJHtMUl/J+lDkj5v+yYV/yN8+4z9cMcrAOSLO14BIGOEPABkjJAHgIwR8gCQMUIeADJGyANAxgh5AMgYIQ8AGfs/PPPtMYuL9vUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sigma = 1.0\n",
    "y_p = 0.0\n",
    "x_p = 5.0\n",
    "plt.plot(x_sample,y_meas,\"bo\",label=\"$\\mathbf{y}$\")\n",
    "plt.errorbar(x_p,y_p,yerr=sigma,fmt=\"ro\",label=\"$\\hat{y}$\",alpha=0.5)\n",
    "plt.legend()\n",
    "plt.xlim(0,10)\n",
    "ylim=plt.ylim()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We model the random vector (a priori) as $$\\begin{bmatrix}\\hat{Y} \\\\ \\mathbf{Y} \\end{bmatrix} \\sim Pr\\left(\\begin{bmatrix}\\hat{Y} \\\\ \\mathbf{Y} \\end{bmatrix} = \\begin{bmatrix}\\hat{y} \\\\ \\mathbf{y} \\end{bmatrix}\\right) = \\mathcal{N}\\left(\\begin{bmatrix}0 \\\\ \\mathbf{0} \\end{bmatrix}, \\begin{bmatrix} K(\\hat{x},\\hat{x}) & K(\\hat{x},\\mathbf{x}) \\\\ K(\\mathbf{x},\\hat{x}) & K(\\mathbf{x},\\mathbf{x}) \\end{bmatrix}\\right) = \\mathcal{N}\\left(\\begin{bmatrix}0 \\\\ \\mathbf{0} \\end{bmatrix}, \\begin{bmatrix} a & \\mathbf{b}^T \\\\  \\mathbf{b} & C \\end{bmatrix}\\right) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-5.570975251039949, 2.7362768371568604)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD8CAYAAACSCdTiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEJtJREFUeJzt3X2MHPV9x/HPx4/Xc6ic2k6d+DifQxEXTFMeFmRqqbgObSmOoLQggVyw2lRXQgikIgoQW22lChM1NClSSipDIJFsXFU8i7Q8mEJQoUDWhGA7NgGnGA6ccjgioBonBr79Y9e9w9zd3t78dufud++XtPLN7Mz3953R7cdzs7M7jggBAPI0rewGAACtQ8gDQMYIeQDIGCEPABkj5AEgY4Q8AGSMkAeAjBHyAJAxQh4AMjajjEHnz58fPT09ZQwNAJPW1q1bX4+IBc2sU0rI9/T0qFqtljE0AExatvc0uw6nawAgY4Q8AGSMkAeAjJVyTh4A2ungwYPq7+/XgQMHym5lTDo6OtTV1aWZM2cWrkXIA8hef3+/jjjiCPX09Mh22e2MKiK0b98+9ff3a8mSJYXrcboGQPYOHDigefPmTfiAlyTbmjdvXrK/Ogh5AFPCZAj4Q1L2SsgDwDBWrKg9JjtCHgAyRsgDwGE2bZKeeEL63veknp7a9GRFyAPAEJs2SX190i9+UZves6c2XTTozz33XM2aNUsDAwOSpMsuu0y2tWvXroIdj46QB4Ah1q6V9u9//7z9+2vzi7jkkkt08OBBbdy4URGhu+66SyeffLJ6e3uLFW6AkAeAIV56qbn5Y7Vy5Ur19vbqlltu0VNPPaWXX35ZF110UbGiY0DIA8AQ3d3NzW/GxRdfrG3btmndunWaOXOmLrjgguJFGygc8raPtP2w7Z22d9i+PEVjAFCGa66ROjvfP6+zsza/qDVr1qizs1NbtmzRqlWrNG/evOJFG0hxJP+OpCsi4hOSlkn6nO1jE9QFgLZbvVrasEGaPbs2vXhxbXr16uK1586dq/PPP1+S2nKqRkrw3TURsVfS3vrPb9neKWmRpB8VrQ0AZVi9WrrxxtrPjzySru7DDz+s559/XgsXLtSqVavSFR5F0i8os90j6QRJT6asCwDtljLcD1m5cqUWLFigG2+8UbNmzUo/wDCShbztD0m6XdIXIuLNYZ7vk9QnSd0p3sEAgEkmIto+ZpKra2zPVC3gN0XEHcMtExEbIqISEZUFC5q6Dy0AYJxSXF1jSd+StDMivla8JQBAKimO5JdLulDSStvP1B9nJqgLACgoxdU1/ylp8nxRMwBMIXziFQAyRsgDQMYIeQDIGCEPAG20bds2LVy4UNu3b2/LeIQ8ALTR+vXr9fjjj2v9+vVtGS/p1xoAQDYO3cU78fcbbN68WZJ06623Jq07Eo7kASBjhDwAHK4Fd/Letm2bli9f/v/TTz/9tFauXFm4biOEPAAM1aI7eS9dulS7d+/Wu+++K0m64oordN111xXttiFCHgCGatGdvKdNm6alS5dqx44duv3229Xd3a0TTzyxUM2x4I1XABiqVXfylrRs2TI99thjuuGGG3TfffcVrjcWHMkDwFAtvJP3smXLtG7dOp1zzjlatGhR4XpjQcgDwFAtvJN3b2+vZs+erSuvvLJwrbEi5AFgqBbeyfv666/Xtddeqzlz5hSuNVaEPAAcbvVqadky6bTTpBdfLBzwu3fvVm9vr95++22tWbMmTY9jxBuvADCchJ90Peqoo7Rr165k9ZrBkTwAZCzVjbxvtv2a7fZ8rRoAYExSHcl/W9IZiWoBABJJEvIR8aikn6WoBQBIh3PyAKaEiCi7hTFL2WvbQt52n+2q7erAwMCoy65YMfhVzrlrxbY2U7NVy6Y03Lhz59YeLR+o/SVKqZ27jo4O7du3b1IEfURo37596ujoSFKvbZdQRsQGSRskqVKpTPw9DSAbXV1d6u/vV6MDzImio6NDXV1dSWpxnTyA7M2cOVNLliwpu41SpLqEcrOk/5J0jO1+259JURcAUEySI/mIuCBFHQBAWlxdAwAZI+QBIGOEPABkjJAHgIwR8gCQMUIeADJGyANAxgh5AMgYIQ8AGSPkASBjhDwAZIyQB4CMEfIAkDFCHgAyRsgDQMYIeQDIGCEPABlLdfu/M2w/Z/sF21elqAkASCAiCj0kTZe0W9LHJc2S9ENJx462zkknnRQj2bgxYvbsCCli8eLadK5asa3N1GzVsikNN+5nP1ubliKmT69Nt2Sg9pcopTYmD0nVaDajm13hAwWkUyXdP2T6aklXj7bOSCG/cWNEZ+fgC1iqTef4C92KbW2mZquWTWm4cWfMeP/0oUehoE+wga3cR1PpdYHRlRXy50q6acj0hZK+Mdo6I4X84sXDv4AXL061iyaOVmxrMzVbtWxKI4073GP69BYM1MQGtnIfTaXXBUY3npB3bb3xs32epD+IiL+oT18o6ZSI+Pxhy/VJ6pOk7u7uk/bs2fOBWtOm1X59PziG9N57hdqccFqxrc3UbNWyKY007kjG/aucYANbuY+m0usCo7O9NSIqzayT4o3XfklHDpnukvTq4QtFxIaIqEREZcGCBcMW6u4efoCR5k9mrdjWZmq2atmUmqk/fXoLBmqigVbuo6n0ukB6KUL++5KOtr3E9ixJ50u6ZzyFrrlG6ux8/7zOztr83LRiW5up2aplUxpu3Bkzhl+2ry/xQE1uYCv30VR6XaAFmj2/M9xD0pmSfqzaVTZrGy3P1TU1XF3TGFfXtL42Jg+VcU5+PCqVSlSr1RGfX7Gi9u8jj7SlnVK1YlubqdmqZVMabty5c2v/vvFGiwdqf4lSamNyKOucPABggiLkASBjhDwAZIyQB4CMEfIAkDFCHgAyRsgDQMYIeQDIGCEPABkj5AEgY4Q8AGSMkAeAjBHyAJAxQh4AMkbIA0DGCHkAyBghDwAZKxTyts+zvcP2e7abulsJAKD1ih7Jb5f0x5IeTdALACCxGUVWjoidkmQ7TTcAgKQ4Jw8AGXNEjL6AvUXSwmGeWhsRd9eXeUTSFyOiOkqdPkl9ktTd3X3Snj17xtszAExJtrdGRFPvfzY8XRMRp4+/pffV2SBpgyRVKpXR/2cBACTB6RoAyFjRSyjPsd0v6VRJ37V9f5q2AAApFL265k5JdybqBQCQGKdrACBjhDwAZIyQB4CMEfIAkDFCHgAyRsgDQMYIeQDIGCEPABkj5AEgY4Q8AGSMkAeAjBHyAJAxQh4AMkbIA0DGCHkAyBghDwAZI+QBIGOEPABkrOg9Xr9qe5ftZ23faXtuqsYAAMUVPZJ/UNJxEfFJST+WdHXxlgAAqRQK+Yh4ICLeqU8+IamreEsAgFRSnpP/c0n/PtKTtvtsV21XBwYGEg4LABjJjEYL2N4iaeEwT62NiLvry6yV9I6kTSPViYgNkjZIUqVSiXF1CwBoSsOQj4jTR3ve9hpJn5b0qYggvAFgAmkY8qOxfYakKyWdFhH707QEAEil6Dn5b0g6QtKDtp+x/c8JegIAJFLoSD4ifiNVIwCA9PjEKwBkjJAHgIwR8gCQMUIeADJGyANAxgh5AMgYIQ8AGSPkASBjhDwAZIyQB4CMEfIAkDFCHgAyRsgDQMYIeQDIGCEPABkj5AEgY4Q8AGSsUMjb/jvbz9Zv/feA7Y+lagwAUFzRI/mvRsQnI+J4SfdK+usEPQEAEikU8hHx5pDJOZKiWDsAgJQK3chbkmxfI+kiST+X9LuFOwIAJNPwSN72Ftvbh3mcLUkRsTYijpS0SdKlo9Tps121XR0YGEi3BQCAETkizRkW24slfTcijmu0bKVSiWq1mmRcAJgqbG+NiEoz6xS9uuboIZNnSdpVpB4AIK2i5+S/YvsYSe9J2iPp4uItAQBSKRTyEfEnqRoBAKTHJ14BIGOEPABkjJAHgIwR8gCQMUIeADJGyANAxgh5AMgYIQ8AGSPkASBjhDwAZIyQB4CMEfIAkDFCHgAyRsgDQMYIeQDIGCEPABkj5AEgY0lC3vYXbYft+SnqAQDSKBzyto+U9HuSXireDgAgpRRH8l+X9CVJkaAWACChQiFv+yxJr0TEDxP1AwBIaEajBWxvkbRwmKfWSvqypN8fy0C2+yT1SVJ3d3cTLQIAxssR4zvLYvs3JT0kaX99VpekVyWdEhE/HW3dSqUS1Wp1XOMCwFRle2tEVJpZp+GR/EgiYpukjwwZ/EVJlYh4fbw1AQBpcZ08AGRs3Efyh4uInlS1AABpcCQPABkj5AEgY4Q8AGSMkAeAjBHyAJAxQh4AMkbIA0DGCHkAyBghDwAZI+QBIGOEPABkjJAHgIwR8gCQMUIeADJGyANAxgh5AMgYIQ8AGSsU8rb/1vYrtp+pP85M1RgAoLgUt//7ekRcl6AOACAxTtcAQMZShPyltp+1fbPtDyeoBwBIpGHI295ie/swj7MlfVPSUZKOl7RX0j+MUqfPdtV2dWBgINkGAABG5ohIU8jukXRvRBzXaNlKpRLVajXJuAAwVdjeGhGVZtYpenXNR4dMniNpe5F6AIC0il5d8/e2j5cUkl6U9JeFOwIAJFMo5CPiwlSNAADS4xJKAMgYIQ8AGSPkASBjhDwAZIyQB4CMEfIAkDFCHgAyRsgDQMYIeQDIGCEPABlL9i2UTQ1qvyXpubYPPDHNl/R62U1MEOyLQeyLQeyLQcdExBHNrJDi9n/j8VyzX5eZK9tV9kUN+2IQ+2IQ+2KQ7aa/o53TNQCQMUIeADJWVshvKGnciYh9MYh9MYh9MYh9MajpfVHKG68AgPbgdA0AZKytIW/7DNvP2X7B9lXtHHsisX2k7Ydt77S9w/blZfdUNtvTbf/A9r1l91Im23Nt32Z7V/3349SyeyqL7b+qvz62295su6PsntrJ9s22X7O9fci8X7P9oO3n6/9+uFGdtoW87emS/knSH0o6VtIFto9t1/gTzDuSroiIT0haJulzU3hfHHK5pJ1lNzEBXC/pvojolfRbmqL7xPYiSZdJqkTEcZKmSzq/3K7a7tuSzjhs3lWSHoqIoyU9VJ8eVTuP5E+R9EJE/CQifinpXySd3cbxJ4yI2BsRT9d/fku1F/Kicrsqj+0uSask3VR2L2Wy/auSfkfStyQpIn4ZEW+U21WpZkj6FdszJHVKerXkftoqIh6V9LPDZp8t6Tv1n78j6Y8a1WlnyC+S9PKQ6X5N4WA7xHaPpBMkPVluJ6X6R0lfkvRe2Y2U7OOSBiTdUj91dZPtOWU3VYaIeEXSdZJekrRX0s8j4oFyu5oQfj0i9kq1g0VJH2m0QjtD3sPMm9KX9tj+kKTbJX0hIt4su58y2P60pNciYmvZvUwAMySdKOmbEXGCpP/VGP4cz1H9XPPZkpZI+pikObb/tNyuJqd2hny/pCOHTHdpiv35NZTtmaoF/KaIuKPsfkq0XNJZtl9U7RTeStsby22pNP2S+iPi0F91t6kW+lPR6ZL+OyIGIuKgpDsk/XbJPU0E/2P7o5JU//e1Riu0M+S/L+lo20tsz1LtTZR72jj+hGHbqp133RkRXyu7nzJFxNUR0RURPar9TvxHREzJI7aI+Kmkl20fU5/1KUk/KrGlMr0kaZntzvrr5VOaom9CH+YeSWvqP6+RdHejFdr2BWUR8Y7tSyXdr9o75TdHxI52jT/BLJd0oaRttp+pz/tyRPxbiT1hYvi8pE31A6GfSPqzkvspRUQ8afs2SU+rdjXaDzTFPvlqe7OkFZLm2+6X9DeSviLpX21/RrX/CM9rWIdPvAJAvvjEKwBkjJAHgIwR8gCQMUIeADJGyANAxgh5AMgYIQ8AGSPkASBj/we2NeEnP0z7YgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.errorbar(x_sample,np.zeros((n,1)),fmt=\"bo\",label=\"$\\mathbf{y}$\",yerr=sigma)\n",
    "plt.errorbar(x_p,0,yerr=sigma,fmt=\"ro\",label=\"$\\hat{y}$\")\n",
    "plt.legend()\n",
    "plt.xlim(0,10)\n",
    "plt.ylim(ylim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exponential Quadratic Kernel\n",
    "We define a Kernel function:\n",
    "\n",
    "$$ \\Sigma_{ij} = K(x_j,x_i) = \\sigma^2 exp\\left(-\\frac{(x_j - x_i)^2}{2\\mathcal{l}^2} \\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5.        , 4.35994902, 0.25926232, 5.49662478, 4.35322393,\n",
       "       4.20367802, 3.30334821, 2.04648634, 6.19270966, 2.99654674,\n",
       "       2.66827275])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_comb = np.append(x_p,x_sample)\n",
    "x_comb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x1a18cc28d0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS0AAAD8CAYAAAAi9vLQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFNZJREFUeJzt3X+s3XV9x/Hnq7eUliLTUnXaVilakaZxYm7wB4niwFHRgC5uAaNTZ+xcBBR1is6ocdkynXO6hDA7RWUyiSIZjemoRgG3RbGFErCtjU3d4EKxBaUwSmnvva/9cU7h3F/nnt7z6f1+v3evR/KN95zz5XU+7YW3n8/n+/l+vrJNRERTzKu6ARERRyNFKyIaJUUrIholRSsiGiVFKyIaJUUrIholRSsijhlJV0vaK+nnU3wuSf8oaZekuyS9bLrMFK2IOJa+Dqzt8vnrgVXtYx1w1XSBKVoRcczY/jHwmy6nXAhc45afAk+X9JxumfNLNnA6S5cM+JQVxxXL2zO8sFgWwKPby9XwZ695vFgWwIHR48vmjSwomjf8QNk8PXKgaN4TyxYXy1r40EixLAAODxeLenzkUQ6NPq5+Ms577WI/9Jve/oy33/XENuBgx1vrba8/iq9bBtzb8Xqo/d6eqf6BWS1ap6w4jp9tWlEs728ePK1YFsCtL1lULOuyG39RLAtg64FTiubduX950bx9nz21aN6im+4omrfr8sFiWauuebhYFoB+3a0jcnR+8uB3+s546Dcj/GzT83o6d+A5vzxou5+/3MkKbNd7C2e1aEVE/RkYZXS2vm4I6OzJLAfu7/YPZE4rIsYw5rBHejoK2AD8Sfsq4iuA/banHBpCeloRMYlSPS1J3wLOBpZKGgI+BRwHYPufgI3A+cAu4ADwrukyU7QiYgxjRgptWWX74mk+N/C+o8lM0YqICUa7z4VXqq85LUlrJe1sr2a9olSjIqI6BkZwT0cVZly0JA0AV9Ja0boauFjS6lINi4jqjOKejir0Mzw8E9hlezeApOtorW7dXqJhEVENA4drvA17P8PDqVayjiFpnaQtkrbsK72SOCKKc49Dw8YND+lxJavt9bYHbQ8+8+SBPr4uImaFYaTHowr9DA+PeiVrRNRfa0V8ffVTtDYDqyStBO4DLgLeWqRVEVEhMTLpQKoeZly0bA9LugTYBAwAV9veVqxlEVGJ1kT8HCxaALY30lqGHxFzRGud1hwtWhExN43O1Z5WRMw96WlFRKMYMVLjXatStCJiggwP2/YMLyy6RfLHl+4slgVwKy8tlnXppncUywLY/eYvF8174uS7i+ad/pb3ls3bcnLRvJU3PlEsa2jtkmJZAEt2nFQsa/TW/vfqN+KQ67sQPD2tiBijtbg0w8OIaJBMxEdEY9hixOlpRUSDjKanFRFN0ZqIr29pqG/LIqISmYiPiMYZyTqtiGiKrIiPiMYZzdXDiGiK1g3TKVoR0RBGHM5tPBHRFDZZXBoRTaIsLo2I5jDpaUVEw2QiPiIawyibAEZEc7QeIVbf0lDflkVERebow1ojYm4yWRH/pEe3z+PWlywqlldyT3eATfffWSzrvOcWi2rlva/sn7W0VfPK/d0BHDy77J93/o9uL5a1fPPCYlkAOnFxsax5j5XZC7/OPa36ltOIqIQtRj2vp6MXktZK2ilpl6QrJvn8eZJulrRV0l2Szu+Wl+FhRIzRmogvcxuPpAHgSuB1wBCwWdIG29s7TvsE8G3bV0laDWwETpkqM0UrIsYpukf8mcAu27sBJF0HXAh0Fi0DR56j9jvA/d0CU7QiYozWRHzPc1pLJW3peL3e9vqO18uAezteDwEvH5fxaeD7ki4FFgPndvvCFK2ImOAoVsQ/aHuwy+eTVT+Pe30x8HXbfy/plcC/SFpje3SywBStiBij8Ir4IWBFx+vlTBz+vRtYC2D7J5IWAkuBvZMFznjgKmlFe8Z/h6Rtkt4/06yIqJdR5vV09GAzsErSSkkLgIuADePOuQc4B0DS6cBCYN9Ugf30tIaBD9m+Q9LTgNsl/WDcVYGIaBgbDo+WmYi3PSzpEmATMABcbXubpM8AW2xvAD4E/LOky2kNHd9pe/wQ8kkzLlq29wB72j8/KmkHrUm3FK2IBmsND8st4bS9kdYyhs73Ptnx83bgrF7zisxpSToFOAO4bZLP1gHrABZyQomvi4hjrM4r4vsuWpJOBL4LfMD2I+M/b1/+XA9wkpZM2eWLiHo4yiUPs66voiXpOFoF61rbN5RpUkRUq+zwsLQZFy1JAr4K7LD9hXJNioiqzdU94s8C3g7cLenILf4fb0+6RURDta4ezsFHiNn+TyZf7RoRDZbtliOicebq8DAi5qA5ffUwIuamOXn1cCaeveZxLrvxF8XyLt30jmJZUHaL5JJbNwOcesOfFc07fm/ZidZTrxkqmkfB7ZEB1txe7j/CHeccXywLYHT/o8WyPDLpxghHl2ExnKIVEU2S4WFENEbmtCKicVK0IqIxsk4rIhon67QiojFsGC60CeCxkKIVERNkeBgRjZE5rYhoHKdoRUSTZCI+IhrDzpxWRDSKGMnVw4hoksxpRURj5N7DiGgWt+a16ipFKyImyNXDiGgMZyI+Ipomw8OIaJRcPWw7MHo8Ww+cUixv95u/XCwL4Lz3vbRYVuk93Xf/Ydk/62GPFM077QXvKZp3+kefKJp31wfKPQDggXWLimUBLNkxXCzLt97af4ZTtCKiYbLkISIaJXNaEdEYRozm6mFENEmNO1rUt5xGRDXaE/G9HL2QtFbSTkm7JF0xxTl/LGm7pG2S/rVbXnpaETFRoa6WpAHgSuB1wBCwWdIG29s7zlkFfAw4y/ZvJT2rW2bfPS1JA5K2Svpev1kRUQ8Fe1pnArts77Z9CLgOuHDcOe8BrrT929Z3e2+3wBLDw/cDOwrkREQNGBgdVU8HsFTSlo5j3bi4ZcC9Ha+H2u91ehHwIkn/JemnktZ2a19fw0NJy4E3AH8NfLCfrIioCQO9r9N60PZgl88nCxo/+JwPrALOBpYD/yFpje2HJwvst6f1ReAjwOhUJ0had6QKP/bbQ31+XUTMBru3owdDwIqO18uB+yc550bbh23/CthJq4hNasZFS9Ibgb22b+92nu31tgdtDy5+xoKZfl1EzCb3eExvM7BK0kpJC4CLgA3jzvk34LUAkpbSGi7uniqwn+HhWcAFks4HFgInSfqm7bf1kRkRlet9OcN0bA9LugTYBAwAV9veJukzwBbbG9qf/YGk7cAI8Be2H5oqc8ZFy/bHaF2mRNLZwIdTsCLmiIKrS21vBDaOe++THT+b1px4T/PiWacVEWMZPDrHb5i2fQtwS4msiKiDOV60ImKOqfHNhylaETFRilZENMbRLS6ddSlaETFBNgFsOzCygDv3Ly+W98TJdxfLKu34vQNF80rv6X6cyrbvWc98pGgeixYWjZv/yMFiWQeXlm3b40vK/S5G5xfqIc31q4cRMbcoPa2IaIzeb9GpRIpWRIyjTMRHRMOkpxURjTLlZlPVS9GKiLGyTisimiZXDyOiWWpctPLcw4holPS0ImKCDA8jojlMbuOJiIZJTysimiTDw4holhStiGiUFK2IaAo5w8OIaJpcPYyIJklPKyKaJUWrZfiBBez77KnF8k5/y3uLZQGsmndnsaxTrxkqlgVw2gveUzSv9J7uP33p9UXzTnv3nxfNW/npzcWybr7olmJZANfuP6NY1pU/K/B7zZxWRDROilZENIlqvAlgdnmIiEZJTysiJsrwMCIaIxPxEdE4NS5afc1pSXq6pOsl/ULSDkmvLNWwiKiQezwq0G9P60vATbbfImkBcEKBNkVEhcQcvXoo6STg1cBXAWwfsv1wqYZFREX81E3T0x29kLRW0k5JuyRd0eW8t0iypMFuef0MD08F9gFfk7RV0lckLZ6kIeskbZG05fChx/r4uoiYNYWGh5IGgCuB1wOrgYslrZ7kvKcBlwG3TZfZT9GaD7wMuMr2GcBjwIQqanu97UHbg8ctmFDTIqKOys1pnQnssr3b9iHgOuDCSc77K+BzwMHpAvspWkPAkO0jlfF6WkUsIhruKIaHS4+MpNrHunFRy4B7O14Ptd976rukM4AVtr/XS9tmPBFv+wFJ90o6zfZO4Bxg+0zzIqJGer8y+KDtbnNQk23M9WS6pHnAPwDv7PUL+716eClwbfvK4W7gXX3mRUTVXPTq4RCwouP1cuD+jtdPA9YAt0gC+F1gg6QLbG+ZLLCvomX7TqDrTH9ENFC5NVibgVWSVgL3ARcBb33ya+z9wNIjryXdAnx4qoIFuWE6IiZRasmD7WHgEmATsAP4tu1tkj4j6YKZtC238UTERAVXu9veCGwc994npzj37OnyUrQiYqwKb9HpRYpWRIwhssvDk/TIARbddEexvNO3nFwsC+Dg2S8tF/aj28tlAad/9ImieSxaWDSu9J7uO//0qqJ5532i3O/2/DvK7tf/zhdOuwi8Zyr05K8UrYholhStiGiUFK2IaIzsXBoRjZOiFRFNUudNAFO0ImKCDA8jojmyuDQiGidFKyKaIiviI6JxNFrfqpWiFRFjZU4rIpomw8OIaJYUrYhokvS0IqJZUrQiojHKPo2nuBStiBgj67Qionlc36qVohURE6Sn1fbEssXsurzcs11X3lh23/T5Bfd1X3N72UdK3vWB5xbNm//IwaJ5Kz+9uWheyT3dATbdf2exrPOWFdqIvW2TTyqWtd8D/YdkcWlENE0m4iOiUVK0IqI5TCbiI6JZMhEfEc1S46LV1yUuSZdL2ibp55K+JansY4sjYtYdWVzay1GFGRctScuAy4BB22uAAeCiUg2LiIrYaLS3owr9Dg/nA4skHQZOAO7vv0kRUbm5ODy0fR/weeAeYA+w3/b3x58naZ2kLZK2jDz22MxbGhGzZq4OD58BXAisBJ4LLJb0tvHn2V5ve9D24MDixTNvaUTMDgOj7u2oQD8T8ecCv7K9z/Zh4AbgVWWaFRGVco9HBfopWvcAr5B0giQB5wA7yjQrIqpUcngoaa2knZJ2Sbpiks8/KGm7pLsk/VDS87vl9TOndRtwPXAHcHc7a/1M8yKiPkpdPZQ0AFwJvB5YDVwsafW407bSWoXwElo15XPdMvtap2X7U7ZfbHuN7bfbLrvtQkTMvl6Hhr31tM4EdtnebfsQcB2tufCnvs6+2faB9sufAsu7BWZFfESM0Vpc2vOE1VJJWzper7fdOeJaBtzb8XoIeHmXvHcD/97tC1O0ImKi3nd5eNB2t03yJtt8bNKK2F59MAi8ptsXpmhFxARH0dOazhCwouP1ciZZhC7pXOAvgddMN81UdnvNiGi+snNam4FVklZKWkDrVr8NnSdIOgP4MnCB7b3TBc5qT2vhQyOsuubhYnlDa5cUywJYvrnc/d47zjm+WBbAA+sWFc07uLTsve03X3RL0bzz73hP0bySWyRvum9rsSyAXx3+32JZb3rDowVSyt1XaHtY0iXAJlr3J19te5ukzwBbbG8A/g44EfhOa/UU99i+YKrMDA8jYqKCmwDa3ghsHPfeJzt+Pvdo8lK0ImKsPKw1Ihon2y1HRKPUt2alaEXERBqt7/gwRSsixjJHs7h01qVoRcQYwiUXlxaXohURE6VoRUSjpGhFRGNkTisimiZXDyOiQZzhYUQ0iEnRioiGqe/oMEUrIibKOq2IaJYUrYhoDBtG6js+TNGKiInS04qIRknRajs8jH79m2JxS3acVCwLQCcuLpY1ur/EXt1PWbJjuGje40sGiuZdu/+MonnvfOFtRfM2udy/KyX3dAdYedyJxbKO12/7DzFQaI/4YyE9rYgYx+DMaUVEU5hMxEdEw2ROKyIaJUUrIpojN0xHRJMYqPHWNPOmO0HS1ZL2Svp5x3tLJP1A0i/b//uMY9vMiJhVdm9HBaYtWsDXgbXj3rsC+KHtVcAP268jYk5o38bTy1GBaYuW7R8D41eEXgh8o/3zN4A3FW5XRFTFYI/2dFRhpnNaz7a9B8D2HknPmupESeuAdQAL55Vb+RsRx1CNV8T3Mjzsi+31tgdtDy6Yt+hYf11ElFDjOa2Z9rR+Lek57V7Wc4C9JRsVERWym331cAobgHe0f34HcGOZ5kRELTS5pyXpW8DZwFJJQ8CngL8Fvi3p3cA9wB8dy0ZGxGwyHhmpuhFTmrZo2b54io/OKdyWiKiDbE0TEY2TrWkioikMOD2tiGgMZxPAiGiYOk/Ey7N42VLSPuB/ejh1KfDgMW7OTNW5bVDv9tW5bVDv9vXatufbfmY/XyTppvb39eJB2+PvTT6mZrVo9UrSFtuDVbdjMnVuG9S7fXVuG9S7fXVu22w75rfxRESUlKIVEY1S16K1vuoGdFHntkG921fntkG921fnts2qWs5pRURMpa49rYiISaVoRUSj1KpoSVoraaekXZJqte+8pBWSbpa0Q9I2Se+vuk3jSRqQtFXS96puy3iSni7pekm/aP8dvrLqNh0h6fL27/Tnkr4laWHF7cnDZLqoTdGSNABcCbweWA1cLGl1ta0aYxj4kO3TgVcA76tZ+wDeD+youhFT+BJwk+0XA79HTdopaRlwGTBoew0wAFxUbavyMJlualO0gDOBXbZ32z4EXEfrARq1YHuP7TvaPz9K6z+6ZdW26imSlgNvAL5SdVvGk3QS8GrgqwC2D9l+uNpWjTEfWCRpPnACcH+VjcnDZLqrU9FaBtzb8XqIGhWFTpJOAc4Abqu2JWN8EfgIUMc7XU8F9gFfaw9fvyJpcdWNArB9H/B5WptZ7gH22/5+ta2a1JiHyQBTPkxmrqtT0dIk79VuPYakE4HvAh+w/UjV7QGQ9EZgr+3bq27LFOYDLwOusn0G8Bg1Gd6054YuBFYCzwUWS3pbta2KbupUtIaAFR2vl1NxN308ScfRKljX2r6h6vZ0OAu4QNJ/0xpW/76kb1bbpDGGgCHbR3qm19MqYnVwLvAr2/tsHwZuAF5VcZsm8+v2Q2T4//4wmToVrc3AKkkrJS2gNRm6oeI2PUmSaM3J7LD9harb08n2x2wvt30Krb+3H9muTW/B9gPAvZJOa791DrC9wiZ1ugd4haQT2r/jc6jJRYJx8jCZttrsp2V7WNIlwCZaV3Cutr2t4mZ1Ogt4O3C3pDvb733c9sYK29QklwLXtv8PaTfwrorbA4Dt2yRdD9xB6wrxViq+ZSYPk+kut/FERKPUaXgYETGtFK2IaJQUrYholBStiGiUFK2IaJQUrYholBStiGiU/wNr+4fqQpIKqAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "l = 1.0\n",
    "X = np.atleast_2d(x_comb).T\n",
    "M=sklearn.metrics.pairwise_distances(X,X)\n",
    "S=sigma*sigma*np.exp(M*M/(-2*l*l))\n",
    "plt.imshow(S)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we want to do is to do is to **update** this model with our data:\n",
    "\n",
    "$$\\hat{y} \\sim Pr(\\hat{Y}=\\hat{y}| \\mathbf{Y} = \\mathbf{y}) = \\mathcal{N}\\left(\\hat{y};\\mathbf{b}^T C^{-1} \\mathbf{y},  a - \\mathbf{b}^T C^{-1} \\mathbf{b} \\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a18e9cf98>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAG+RJREFUeJzt3X90VeWd7/H3NwkJBn8AMYAmQFCpKNOqGL1UXK4O2FZlWnBG7sJhSsZyb5ZTb2vveJc/qp3eWdNaO7drrK5bXQu1ipdox+UvGEu1Fn+NP7AGRREB+aFAACEiYivyI8n3/vE8pyeEQE7ICSf4fF5r7bX3fp5n7/OcnXM+e5/nnJNj7o6IiKShqNAdEBGRw0ehLyKSEIW+iEhCFPoiIglR6IuIJEShLyKSEIW+iEhCFPoiIglR6IuIJKSk0B0AOP74472mpqbQ3RAROaIsXrz4Q3ev7M42fSL0a2pqaGxsLHQ3RESOKGa2rrvbdDm8Y2anmtmSdtMnZvZ9MxtsZk+b2ao4HxTbm5ndbmarzewtMxt3KHdGRETyr8vQd/eV7n6mu58JnA3sBB4DrgcWuvtoYGFcB7gYGB2neuDO3ui4iIh0X3ffyJ0ErHH3dcAUYE4snwNMjctTgPs9WAQMNLMT8tJbERHpke6O6U8HHozLQ919M4C7bzazIbG8CtjQbpumWLa5Jx0VEcmnvXv30tTUxK5duwrdlS7179+f6upq+vXr1+N95Rz6ZlYKfBO4oaumnZTt90/7zayeMPzDiBEjcu2GiEheNDU1ccwxx1BTU4NZZ7HVN7g727Zto6mpiVGjRvV4f90Z3rkYeN3dt8T1LZlhmzjfGsubgOHttqsGNnXcmbvPdvdad6+trOzWJ45Eel9DA9TUQFFRmDc0FLpHkme7du2ioqKiTwc+gJlRUVGRt1ck3Qn9y8kO7QDMB+rich0wr135zPgpnvHAjswwkMgRoaEB6uth3TpwD/P6egX/51BfD/yMfPYzp9A3s3Lgq8Cj7YpvAb5qZqti3S2xfAGwFlgN3AV8J2+9FTkcbrwRdu7ct2znzlAucoTLaUzf3XcCFR3KthE+zdOxrQNX5aV3IoWwfn33ykWOIPrfOyIdHeiDBfrAQdI+L2/zKPRFOvrJT6C8fN+y8vJQLknqrbd5LrvsMkpLS2lubgbge9/7HmbGihUr8tDrzin0RTqaMQNmz4aRI8EszGfPDuWSpN56m+c73/kOe/fuZe7cubg7jz/+OOeccw5jxozp2Y4Pok/8wzWRPmfGDIW8/Flvvc0zceJExowZw7333st5553Hhg0buPbaa3u20y7oSl9EpAu9+TbPlVdeydKlS7npppvo168fl19+ec93ehAKfRGRLvTm2zx1dXWUl5fz+9//nsmTJ1NRUdH1Rj2g0BcR6UJvvs0zcOBApk+fDsDMmTN7vsMuaExfRCQHvfU2z7PPPsuqVasYNmwYkydPzv8NdKDQFxEpoIkTJ1JZWcldd91FaWlpr9+eQl9EpIDCPzE4fDSmLyKSEIW+iEhCFPoiIglR6IuIJEShLyKSEIW+iEhCFPoiIglR6IuIFNjSpUsZNmwYb7/9dq/flkJfRKTAbr75Zl5++WVuvvnmXr8tfSNXRCRXX/lKmD/3XF53++CDDwLwwAMP5HW/ncnpSt/MBprZw2a2wsyWm9mXzWywmT1tZqvifFBsa2Z2u5mtNrO3zGxc794FERHJVa7DO7cBT7r7GOAMYDlwPbDQ3UcDC+M6wMXA6DjVA3fmtcciIoXQ0ACLFsHzz+ftl9GXLl3KhAkT/rz++uuvM3HixB7v92C6DH0zOxa4ALgHwN33uPvHwBRgTmw2B5gal6cA93uwCBhoZifkveciIodL5pfRd+8O63n6ZfSxY8eyZs0aWltbAbjmmmv4+c9/3tPeHlQuV/onAc3AvWb2hpndbWYDgKHuvhkgzofE9lXAhnbbN8UyEZEjUy/9MnpRURFjx45l2bJlPPLII4wYMYJx43p3RDyXN3JLgHHAd939VTO7jexQTmesk7L9/neomdUThn8YkY8fmhQR6S299cvowPjx43nppZe44447ePLJJ3u8v67kcqXfBDS5+6tx/WHCSWBLZtgmzre2az+83fbVwKaOO3X32e5e6+61lZWVh9p/EZHe14u/jD5+/HhuuukmLr30Uqqqen9QpMvQd/cPgA1mdmosmgS8A8wH6mJZHTAvLs8HZsZP8YwHdmSGgUREjki9+MvoY8aMoaysjOuuu67H+8pFrp/T/y7QYGalwFrgCsIJ4yEzmwWsB6bFtguAS4DVwM7YVkTkyJX5cdxZs8KbuSNHhsDPw4/m3nbbbfz0pz9lwIABPd5XLnIKfXdfAtR2UjWpk7YOXNXDfomI9C0zZsBdd4XlPHw5a82aNUyePJkJEyZQV1fX9QZ5om/kiojkKo/fxD355JNZsWJF3vaXK/3vHRGRhCj0RUQSotAXEUmIQl9EJCEKfRFJVviwYd+Xz34q9EUkSf3792fbtm19PvjdnW3bttG/f/+87E8f2RSRJFVXV9PU1ERzc3Ohu9Kl/v37U11dnZd9KfRFJEn9+vVj1KhRhe7GYafhHRGRhCj0RUQSotAXEUmIQl9EJCEKfRGRhCj0RUQSotAXEUmIQl9EJCEKfRGRhCj0RUQSotAXEUlITqFvZu+b2VIzW2JmjbFssJk9bWar4nxQLDczu93MVpvZW2Y2rjfvgIiI5K47V/p/6e5nunttXL8eWOjuo4GFcR3gYmB0nOqBO/PVWRER6ZmeDO9MAebE5TnA1Hbl93uwCBhoZif04HZERCRPcg19B35nZovNrD6WDXX3zQBxPiSWVwEb2m3bFMv2YWb1ZtZoZo1Hwv+zFhH5PMj1/+lPcPdNZjYEeNrMVhykrXVStt9P07j7bGA2QG1tbd/+6RoRkc+JnK703X1TnG8FHgPOBbZkhm3ifGts3gQMb7d5NbApXx0WEZFD12Xom9kAMzsmswx8DXgbmA/UxWZ1wLy4PB+YGT/FMx7YkRkGEhGRwspleGco8JiZZdo/4O5PmtlrwENmNgtYD0yL7RcAlwCrgZ3AFXnvtYiIHJIuQ9/d1wJndFK+DZjUSbkDV+WldyIiklf6Rq6ISEIU+iIiCVHoi4gkRKEvIpIQhb6ISEIU+iIiCVHoi4gkRKEvIpIQhb6ISEIU+iIiCVHoi4gkRKEvIpIQhb6ISEIU+iIiCVHoi4gkRKEvIpIQhb6ISEIU+iIiCVHoi4gkJOfQN7NiM3vDzJ6I66PM7FUzW2Vm/25mpbG8LK6vjvU1vdN1ERHpru5c6V8NLG+3/jPgVncfDWwHZsXyWcB2dz8FuDW2ExGRPiCn0DezamAycHdcN2Ai8HBsMgeYGpenxHVi/aTYXkRECizXK/1fANcCbXG9AvjY3VviehNQFZergA0AsX5HbC8iIgXWZeib2V8BW919cfviTpp6DnXt91tvZo1m1tjc3JxTZ0VEpGdyudKfAHzTzN4Hfk0Y1vkFMNDMSmKbamBTXG4ChgPE+uOAjzru1N1nu3utu9dWVlb26E6IiEhuugx9d7/B3avdvQaYDjzj7jOAZ4HLYrM6YF5cnh/XifXPuPt+V/oiInL49eRz+tcB/2hmqwlj9vfE8nuAilj+j8D1PeuiiIjkS0nXTbLc/Tngubi8Fji3kza7gGl56JuIiOSZvpErIpIQhb6ISEIU+iIiCVHoi4gkRKEvIpIQhb6ISEIU+iIiCVHoi4gkRKEvIpIQhb6ISEIU+iIiCVHoi4gkRKEvIpIQhb6ISEIU+iIiCVHoi4gkRKEvIpIQhb6ISEIU+iIiCVHoi4gkpMvQN7P+ZvYHM3vTzJaZ2T/H8lFm9qqZrTKzfzez0lheFtdXx/qa3r0LIiKSq1yu9HcDE939DOBM4CIzGw/8DLjV3UcD24FZsf0sYLu7nwLcGtuJiEgf0GXoe/CnuNovTg5MBB6O5XOAqXF5Slwn1k8yM8tbj0VE5JDlNKZvZsVmtgTYCjwNrAE+dveW2KQJqIrLVcAGgFi/A6joZJ/1ZtZoZo3Nzc09uxciIpKTnELf3Vvd/UygGjgXOK2zZnHe2VW971fgPtvda929trKyMtf+iohID3Tr0zvu/jHwHDAeGGhmJbGqGtgUl5uA4QCx/jjgo3x0VkREeiaXT+9UmtnAuHwUcCGwHHgWuCw2qwPmxeX5cZ1Y/4y773elLyIih19J1004AZhjZsWEk8RD7v6Emb0D/NrMfgy8AdwT298D/D8zW024wp/eC/0WEZFD0GXou/tbwFmdlK8ljO93LN8FTMtL70REJK/0jVwRkYQo9EVEEqLQFxFJiEJfRCQhCn0RkYQo9EVEEqLQFxFJiEJfRCQhCn0RkYQo9EVEEqLQFxFJiEJfRCQhCn0RkYQo9EVEEqLQFxFJiEJfRCQhCn0RkYQo9EVEEqLQFxFJSJehb2bDzexZM1tuZsvM7OpYPtjMnjazVXE+KJabmd1uZqvN7C0zG9fbd0JERHKTy5V+C3CNu58GjAeuMrPTgeuBhe4+GlgY1wEuBkbHqR64M++9FhGRQ9Jl6Lv7Znd/PS7/EVgOVAFTgDmx2RxgalyeAtzvwSJgoJmdkPeei4hIt3VrTN/MaoCzgFeBoe6+GcKJARgSm1UBG9pt1hTLRESkwHIOfTM7GngE+L67f3Kwpp2UeSf7qzezRjNrbG5uzrUbIiLSAzmFvpn1IwR+g7s/Gou3ZIZt4nxrLG8ChrfbvBrY1HGf7j7b3WvdvbaysvJQ+y8iIt2Qy6d3DLgHWO7u/9auaj5QF5frgHntymfGT/GMB3ZkhoFERKSwSnJoMwH4FrDUzJbEsh8AtwAPmdksYD0wLdYtAC4BVgM7gSvy2mMRETlkXYa+u79I5+P0AJM6ae/AVT3sl4iI9AJ9I1dEJCEKfRGRhCj0RUQSotAXEUmIQl9EJCEKfRGRhCj0RUQSotAXEUmIQl9EJCEKfRGRhCj0RUQSotAXEUmIQl9EJCEKfRGRhCj0RUQSotAXEUmIQl9EJCEKfRGRhCj0RUQSotAXEUlIl6FvZr8ys61m9na7ssFm9rSZrYrzQbHczOx2M1ttZm+Z2bje7LyIiHRPLlf69wEXdSi7Hljo7qOBhXEd4GJgdJzqgTvz000REcmHLkPf3V8APupQPAWYE5fnAFPbld/vwSJgoJmdkK/OiohIzxzqmP5Qd98MEOdDYnkVsKFdu6ZYJiIifUC+38i1Tsq804Zm9WbWaGaNzc3Nee6GiIh05lBDf0tm2CbOt8byJmB4u3bVwKbOduDus9291t1rKysrD7EbIiLSHYca+vOBurhcB8xrVz4zfopnPLAjMwwkIiKFV9JVAzN7EPgKcLyZNQE/Am4BHjKzWcB6YFpsvgC4BFgN7ASu6IU+i4jIIeoy9N398gNUTeqkrQNX9bRTIiLSO/SNXBGRhCj0RUQSotAXEUmIQl9EJCEKfRGRhCj0RUQSotAXEUmIQl9EJCEKfRGRhCj0RUQSotAXEUmIQl9EJCEKfRGRhCj0RUQSotAXEUmIQl9EJCEKfRGRhCj0RUQSotAXEUmIQl9EJCG9EvpmdpGZrTSz1WZ2fU/21dAANTVQVBTmDQ356WNfle/72539pXasj2R98m/V1ga7doF7WP/4Y1i7FlasgKVLYckSWLw4237VKnj+eXjmGfjd7+C3v4WnnsrWNzbC44+H6bHHwrRgQbb+5ZdD2bx58B//EepeeCFb/+abYf3FF2HRonDb776brf/gA9i0CT78EHbsgJ07oaWld45NH1KS7x2aWTHwS+CrQBPwmpnNd/d3uruvhgaorw9/C4B168I6wIwZ+epx35Hv+9ud/XX3thsa4MYbYf16GDECfvKTw/M36ex2oTB96UpvHaO8PU7cYc+eMO3enZ1XVUFpaQjFtWv3rduzByZPhqOOCqH88svZur17w35vuAHKyuAPfwj1HZ11VjhbrVwZ9lFUBMXFYSorg69/PbRrbob33gOz7LYDBmSXt20Ld76tLUytrTBwIFxwQah/4w14//19b/vEE+ELXwjL8+bB5s371o8cCVdcEZbnzg0HuawM+vcP8+pqqK0N9cuWhT6Xl2en/v3D/enDzDNn5Xzt0OzLwP9296/H9RsA3P2nB9qmduxYb5wzJ6xk+jNoEDUXnsK6dVDLa5QQzsCGM2woPPriEDjllND2P/8zbJfZ1j38cU45JTwQnn02W56ZTjoJRo8OVyYLF+5b5w6nnx7q//SncPXRsf7ss8P+P/ooXKG0v213OP/8cBubN2fr209f+1q4RFu3Dn7zG3Dnn37obN/uGM6j/DUbqeYLrOTvBv2GH94Y919UFJ4E06fDsGHwzjvw3HOhzCxbP20aNWcOZMC6ZdTSiGO0UUQbRQw+vpj/u+4b4UH67ruwahV13y5m09ZiWgnTS0yglRLGV23glfnN2SdlcTHzf1PM5f80mp2fGcfxMUfxGWVHFfN/bi1h2vRiKCnJPjnd933S9kDHsIOQTS0t4Tmf0a8f3HtvYYO/s76Wl8Ps2T3oV1sb7NnDyaeWsHZ9CWXsYhgfUMoeStlDzbDdPPHoHjj1VBg8OFzFvvLK/qH+N38Twm/JknAV3dGVV4bH1muvhccmhANdWhqCb+ZMOO44WL48PP4y5Zl5bW34I2zbFp4/paVhPfMYOuaY8Jhobc0+ZnuDe/aE0NISTkruoe8AGzaEP9DeveHY7N0LRx8NY8eG+gULYPv2kBG7d4f5SSfB1Kmh/pZbQll7Z5wBl14alu+7L9z39ieF4cPDicU9vMIoLw8n0EM8Bma22N1ru7NN3q/0gSpgQ7v1JuC/dGxkZvVAPcCYQYPgiSf2bTB2LOvXh1C/kN9Txu7stluAd8ZlQ3/hwv17MX58qG9rCw/8TChmpqOOCqHe2hoeuB3rq6rCflpaYOPG/eszf+y2Nvj00/3rMylUVBRCsGN9STz0/fpBRQWYsWK70YbhGHsoBeBPHM3i7SfD6BicmQdyWVlYLysL22dOJm1tYV5UxPr1cDpOMa0YHiO/Df+wNXuSam6GZcuo2trK8Bj5hvMKX6YVOG7jOzD/lX0O7dpfwM7PfgTABbzAOF6Hz2DjdcBmwgP9Bz8IjR99NHtFVFIS5scem700feopaGrat37gQLj44lC/aFEYJigu5olrSjh7ZzE7OI63OAOAkXvepYzdtFJMCyW0Usynewdw9dXDQrh+9FHYT2bfmavK0nB82Z19XP355JT5m0F4fLQ/tm1t2e3dw7BApjzTprycG288hl07WzmZ9yihJUw7W/j1/2phxl+eGEL3s8/gpZfCYywTSnv2wLhx4bHZ3AwPPJC9Go9X0kev/2vgSwxlC3/PfdnufwA8RTh+gweH+7ZxYzaMBwwI5Zn7VlUVLj7aB3ppKQwaFOrPPBO+9KVQ1tmJ+7TTwnQgFRVhOpDi4gPX5YNZ9kTTr194zrc3fPjBt7/kkoPXX3llOGl89lmY79yZvb+treFx9Mkn4RVTZujo/PND6O/aBb/8ZbafRx0V/j7nnRdeCe3aFV4llZeH8sxJY9Cg8GqiB3rjSn8a8HV3/29x/VvAue7+3QNtUztunDc+/3xmB38OxZpTy1i3DsoIAeuEB96IEcaqNR2emO23zSwfYTIX/h2NHLn/q9R87699W6MtHmtj9IjdvLt0dzjGcaoa1somTgRgCFsYxHaKaaWEVt5cHMdEx40L82XLwoO+tTU86FtbQ4hkXsI/80wI/fb1xx4Lf/u3ob6hIYyPtLbyzz9swYGNVHE3/x2AK7mToWzZ5768xyjupy6c1267LVyttTdmTHilBPCv/7rv5Tjse7X24x/vP857zjlhiKO1Ff7lX/Y/wOefT9HXLqTMP+M6frZPlQE/en5iGILYsQNuvz0EUklJmMrKQjB88YshMBYuzIZynMb911N4o6mSMnZxApvjdX4pw4aXsvTdGNx9fIghSZlXGqWlYXnFiuzJ4tNPw/yLXwwn0uZmuOOO7MVZxje+EUYZPvgAfvtb7Nvf7vaVft8Y3qmt9cbGxv3Ke+Ulch+W7/vbnf11p22+T065qhnpbFgfXq200A+AY/iEUvbEk04LxbSyh1K2MCw8X1auDFdNmRNKW1u4Wjr11LDT117LjkVnnguVldlx3xdfzA5RFRWFaehQGDUqlC9Zki3PTBUV1Jw7hPXr2qhiY+Y6nxZKOHF4CcvW9A9Bf4hSe14ky33fVxGffhqG3QYNCkNDb72FTZrU7dDH3fM6EYaM1gKjgFLgTWDswbY5++yz/UDmznUfOdLdLMznzj1g08+FfN/f7uwv17Zz57qXl+/7JkV5ee//bTq73QNNFRW925dD6Ws+j1FqzwvpHNDo3c3o7m6Q007hEuBdYA1wY1ftDxb60jcVKnQ63u4//IN7aem+4Vpa2jdCUMEsve1QQj/vwzuH4kDDOyK5KNTHR0UKra98ekfksJoxQyEvkiu9xS8ikhCFvohIQhT6IiIJUeiLiCREoS8ikpA+8ZFNM/sjsLLQ/egjjgc+LHQn+ggdiywdiywdi6xT3f2Y7mzQVz6yubK7nzX9vDKzRh2LQMciS8ciS8ciy8y6/QUnDe+IiCREoS8ikpC+EvqzC92BPkTHIkvHIkvHIkvHIqvbx6JPvJErIiKHR1+50hcRkcOg4KFvZheZ2UozW21m1xe6P4ViZsPN7FkzW25my8zs6kL3qZDMrNjM3jCzJ7pu/flmZgPN7GEzWxEfH18udJ8Kwcz+Z3xuvG1mD5pZz3438AhjZr8ys61m9na7ssFm9rSZrYrzQV3tp6Chb2bFwC+Bi4HTgcvN7PRC9qmAWoBr3P00YDxwVcLHAuBqYHmhO9FH3AY86e5jgDNI8LiYWRXwPaDW3f8CKAamF7ZXh919wEUdyq4HFrr7aGBhXD+oQl/pnwusdve17r4H+DUwpcB9Kgh33+zur8flPxKe2FWF7VVhmFk1MBm4u9B9KTQzOxa4ALgHwN33uPvHhe1VwZQAR5lZCVAObCpwfw4rd38B+KhD8RRgTlyeA0ztaj+FDv0qYEO79SYSDbr2zKwGOAt4tbA9KZhfANcCbYXuSB9wEtAM3BuHu+42swGF7tTh5u4bgZ8D64HNwA53/11he9UnDHX3zRAuHIEhXW1Q6NC3TsqS/jiRmR0NPAJ8390/KXR/Djcz+ytgq7svLnRf+ogSYBxwp7ufBXxKDi/hP2/iWPUUwm9vnwgMMLO/K2yvjkyFDv0mYHi79WoSe8nWnpn1IwR+g7s/Wuj+FMgE4Jtm9j5huG+imc0tbJcKqglocvfMq76HCSeB1FwIvOfuze6+F3gUOK/AfeoLtpjZCQBxvrWrDQod+q8Bo81slJmVEt6YmV/gPhWEmRlh3Ha5u/9boftTKO5+g7tXu3sN4fHwjLsne0Xn7h8AG8zs1Fg0CXingF0qlPXAeDMrj8+VSST4hnYn5gN1cbkOmNfVBgX9h2vu3mJm/wN4ivBu/K/cfVkh+1RAE4BvAUvNbEks+4G7Lyhgn6Rv+C7QEC+M1gJXFLg/h527v2pmDwOvEz7p9gaJfTPXzB4EvgIcb2ZNwI+AW4CHzGwW4cQ4rcv96Bu5IiLpKPTwjoiIHEYKfRGRhCj0RUQSotAXEUmIQl9EJCEKfRGRhCj0RUQSotAXEUnI/wcLj8W8sbupcQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = S[0,0]\n",
    "b = S[0,1:]\n",
    "C = S[1:,1:]\n",
    "Cinv = np.linalg.inv(C)\n",
    "mu_p = b.dot(Cinv).dot(y_meas)\n",
    "sig_p = a - b.dot(Cinv).dot(b)\n",
    "plt.plot(x_sample,y_meas,\"bo\",label=\"$\\mathbf{y}$\")\n",
    "plt.errorbar(x_p,mu_p,yerr=sig_p,fmt=\"ro\",label=\"$\\hat{y}$\")\n",
    "plt.legend()\n",
    "plt.xlim(0,10)\n",
    "ylim=plt.ylim()\n",
    "plt.plot(X_space,y_real,\"r--\",alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x1a18e27208>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASYAAAEDCAYAAACYpalfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEQFJREFUeJzt3X+MpVV9x/H3ZwcFGaRg1wbdWQtaoxJai5mgSGItYItK4J82wRbT2iYbm6porFvRJjb9qwlWIaklmSKaViK2SFNiNoAGSdtEKcOPqrBSEC0MYGHFX9lWdnfm0z/u3XRcZ+c+wzw/zp3zeSVP2HvvM2e+THY+e855znke2SYioiTbhi4gIuJICaaIKE6CKSKKk2CKiOIkmCKiOAmmiChOgiliC5N0raQnJX2jwbkvkfRlSfdI+pqkt/RR41oSTBFb26eBCxqe+2fAP9g+E7gE+JuuipokwRSxhdn+F+Dp1e9JepmkmyXdJelfJb3y8OnAieM//xzweI+l/pRjhvrGETGYBeCdth+U9FpGPaNzgT8HbpX0bmAWOH+oAhNMERWRdALweuAfJR1++9jxf98GfNr2X0k6G/h7SWfYXum7zgRTRF22AT+w/atrfPaHjOejbH9F0nHAduDJHusDMscUURXbPwK+Lem3ATTy6vHHjwDnjd9/FXAc8NQQdSp3F4jYuiR9Fngjo57PfwMfAW4DrgZeBDwHuN72X0g6Hfhb4ARGE+G7bd86SN0JpogoTYZyEVGcTia/X/CCbZ6bm2m93Qe+f0rrbXbl2KX9nbT7zM7ZTtolHeepcujpp1nev1+Tzzy63/z1WX/v6eVG5971tWdusd10oeamdRJMc3Mz3LRne+vtnvu597XeZlde9oGvdNLug7tf20m7Mz/ppvPsTf3q9EwdpXMHP4THrvz4ptv43tPL/PstL2l07syLHmz/F3odWS4QUSkDK/S+RKmRBFNEpYw56GZDub4lmCIqlh5TRBTFmOVClwslmCIqtlLo5dgEU0SlDCwXGkyNrhFLukDSA5IekvTBrouKiH6s4EZH3yb2mCTNAJ8A3gQsAXdKusn2/V0XFxHdMXCw0DmmJj2ms4CHbD9s+wBwPXBxt2VFRNeMWW549K1JMO0AHl31emn83k+RtEvSoqTF7z1d5iXIiFjFsNzw6FuTYFprPf3PlGp7wfa87fmff0H2BkeUbrTyu9nRtyZX5ZaAnatezzHgTcojoi1iec1+x/CaBNOdwMslnQY8xuixLr/TaVUR0bnR5PeUBpPtQ5LeBdwCzADX2r6v88oiolOjdUxTGkwAtvcAezquJSJ6tjKtPaaI2JqmvscUEVuPEcuF3l07wRRRsQzlIqIoRhxw+/fmb0OCKaJSowWW7Q3lxvtqF4HHbF+4mbY6CaYHvn/KVD04oAvfuuLsTtrd9kwnzU7XQwO6UuEPoeXJ78uAvcCJm22ozJmviOicLZa9rdExiaQ54K3ANW3UlqFcRMVW2usxXQnsBp7fRmMJpohKjSa/G0fAdkmLq14v2F4AkHQh8KTtuyS9sY3aEkwRldrg5Pc+2/NH+ewc4CJJbwGOA06U9Bnblz7b2jLHFFGxZavRsR7bl9ues30qo03+t20mlCA9pohqZeV3RBRppcEVt42wfTtw+2bbSTBFVGq0iTc9pogoiBEHsyUlIkpi02jx5BASTBHVUpsLLFuVYIqolEmPKSIKlMnviCiKUW4UFxFlGT2+qcwIKLOqiOjBdD/wMiK2INP+yu+2JJgiKpYeU0QUxVZ6TBFRltHkd7akRERRlAWWEVGW0eR35pgiojBZ+R0RRcnK74goUptP4m1TgimiUjYcXEkwRURBRkO5BFNEFCYrvyOiKCUvF5jYj5O0U9KXJe2VdJ+ky/ooLCK6NhrKNTn61qTHdAh4v+27JT0fuEvSF23f33FtEdGxqb3nt+0ngCfGf/6xpL3ADiDBFDHFRlfltsBeOUmnAmcCd6zx2S5gF8DMySe3UFpEdKnkBZaNB4+STgA+D7zX9o+O/Nz2gu152/Mzs7Nt1hgRHVkZP8Jp0tG3Rj0mSc9hFErX2b6x25Iiog8lX5WbGEySBHwS2Gv7Y92XFBF9meYFlucAbwe+Lune8Xsfsr2nu7Iiomu2ODStwWT736DQa4oRsSlTO5SLiK1pqueYImLrSjBFRFFKXseUYIqo2NRuSXk2jl3az8s+8JXW2/3WFWe33mZXHvrdqztp9xWf+qNO2l0+1p20S6H/Iq9FXf0ItnXU8CbZcKiFG8VJ2gn8HXAKsAIs2L5qM22mxxRRsZaGcq1v9E8wRVSqrTmmLjb6J5giKubmwbRd0uKq1wu2F448ab2N/huRYIqo2AYmv/fZnl/vhEkb/TciwRRRKbu9dUxtb/RPMEVUSyy3c1Wu9Y3+Ze7gi4he2Gp0THB4o/+5ku4dH2/ZTF3pMUVUqq29cl1s9E8wRdTKo3mmEiWYIipW1ZaUiCifW5r87kKCKaJiGcpFRHE2sPK7VwmmiErZCaaIKFBuFBcRxckcU0QUxYiVXJWLiNIU2mFKMEVUK5PfEVGkQrtMCaaIilXVY3pm5ywP7n5t6+1ue6b1JjuTp5lMn85+BIX+bA2srJRZW3pMEbUyxYZmgimiYlnHFBHlSTBFRFka3TZ3EAmmiJqlxxQRRTE4V+UiojxlBlPjHXySZiTdI+kLXRYUET1yw6NnG9lafBmwt6tCImIA0xxMkuaAtwLXdFtORPTm8ALLJkfPms4xXQnsBp5/tBMk7QJ2AcycfNLmK4uIzpW6wHJij0nShcCTtu9a7zzbC7bnbc/PnHBCawVGRIdW1OzoWZMe0znAReNnkR8HnCjpM7Yv7ba0iOiaprXHZPty23O2TwUuAW5LKEVsAU0nvgcIr6xjiqjWMBPbTWwomGzfDtzeSSUR0b9Ch3LpMUXUbGXoAtaWYIqoVW4UFxElKvWqXIIpomaFBlOZj+GMiKp102MyzPyk/cwrdDi8pjzNJKZBW0M5SRcAVwEzwDW2/3Iz7aXHFFEr08qWFEkzwCeANwOnA2+TdPpmSkswRdSsnZXfZwEP2X7Y9gHgeuDizZSVYIqomNzsALZLWlx17FrVzA7g0VWvl8bvPWu5KhdRs+ZzTPtszx/ls7XGepuavUowRdSsncnvJWDnqtdzwOObaTBDuYhKNR3GNbhydyfwckmnSXouo7uQ3LSZ2tJjiqhZCzeBs31I0ruAWxgtF7jW9n2baTPBFFGxttYx2d4D7GmntQRTRN0K3ZKSYIqoVbP5o0EkmCJqlmCKiNKo0BvFZblARBQnPaaImmUoFxFFyeR3RBQpwRQRxUkwRURJRLlX5RJMEbXKHFNEFCnBFBHFqS2Yqn+YR/U/gJgGGcpFRHkSTBFRFOeqXESUKD2miChN5pgiojwJpogoSrOn7A4iwRRRKVHuUK7RjeIknSTpBknflLRX0tldFxYR3WvpuXKta9pjugq42fZvjR9od3yHNUVEXwrtMU0MJkknAm8Afh/A9gHgQLdlRUQvCg2mJkO5lwJPAZ+SdI+kayTNHnmSpF2SFiUtLu/f33qhEdGy9h4R3romwXQM8BrgattnAvuBDx55ku0F2/O252dmfya3IqJEbnj0rEkwLQFLtu8Yv76BUVBFxJTTSrOjbxODyfZ3gUclvWL81nnA/Z1WFRG9KHUo1/Sq3LuB68ZX5B4G3tFdSRHRi2lfYGn7XmC+41oiom/THEwRsfWUvPI7wRRRMa2UmUwJpohaTfscU0RsTaUO5Rpt4o2ILaqHBZaSrhjfAOBrkv5J0kmTvibBFFGxntYxfRE4w/avAP8JXD7pCxJMETXrocdk+1bbh8YvvwrMTfqazDFF1GpjT0nZLmlx1esF2wvP4rv+AfC5SSclmCIqtcF1TPtsH3WRtaQvAaes8dGHbf/z+JwPA4eA6yZ9swRTRM3czmU52+ev97mk3wMuBM6zJ3/TBFNExfpYLiDpAuBPgV+z/T9NvibBFFGr/hZY/jVwLPBFSQBftf3O9b4gwRRRsT7utWT7lzb6NQmmiIoNcRO4JhJMEbUyrU1+ty3BFFGxUvfKJZgiapZgioiS5EZxEVEeOzeKi4gClZlLCaaImmUoFxFlMZChXEQUp8xcSjBF1CxDuYgoTq7KRURZqnx8Uxd9RKv9NjvSVRd5in4EUbjRAssykyk9poia5e4CEVGa9JgioixVzjFFROGyVy4iSpShXEQUZWMPvOxVgimiZoX2mLY1OUnS+yTdJ+kbkj4r6biuC4uIHrjh0bOJwSRpB/AeYN72GcAMcEnXhUVE97Sy0ujoW9Oh3DHA8yQdBI4HHu+upIjohSl2geXEHpPtx4CPAo8ATwA/tH3rkedJ2iVpUdLi8v797VcaEa0SRm529K3JUO5k4GLgNODFwKykS488z/aC7Xnb8zOzs+1XGhHts5sdPWsy+X0+8G3bT9k+CNwIvL7bsiKiF4UGU5M5pkeA10k6Hvhf4DxgsdOqIqJ7Bc8xTQwm23dIugG4GzgE3AMsdF1YRHRviCtuTTS6Kmf7I8BHOq4lIno1zDCtiaz8jqiVSTBFRIHKHMklmCJqlhvFRUR5EkwRURQblsscy3UXTJU/zsPb8piUmAI99pgk/QlwBfBC2/vWOzc9poia9RRMknYCb2K0YHuiRvdjiogtyMCKmx2b93FgNw3v7pQeU0S1DG48x7Rd0uqtaAu2G+0AkXQR8Jjt/5CaTUUkmCJqZTYy+b3P9vzRPpT0JeCUNT76MPAh4Dc2UlqCKaJmLc0x2T5/rfcl/TKjWyYd7i3NAXdLOsv2d4/WXoIpomYdT37b/jrwC4dfS/oOo9t056pcRKwlm3gjojQGer7tie1Tm5yXYIqoWXpMEVGWGrekRETZDG6+jqlXCaaImrWzqrt1CaaImmWOKSKKYvd+Va6pBFNEzdJjioiyGC8vD13EmhJMEbU6fNuTAiWYImqW5QIRURIDTo8pIoriDd0orlcJpoiKlTr5LXdwuVDSU8B/NTh1O7DufVkKM031TlOtMF31llDrL9p+4WYakHQzo/+XJvbZvmAz328jOgmmxt9cWlzvdp2lmaZ6p6lWmK56p6nWaZWnpEREcRJMEVGcoYOp0eNfCjJN9U5TrTBd9U5TrVNp0DmmiIi1DN1jioj4GQmmiCjOYMEk6QJJD0h6SNIHh6pjEkk7JX1Z0l5J90m6bOiampA0I+keSV8Yupb1SDpJ0g2Svjn+GZ89dE3rkfS+8d+Db0j6rKTjhq5pKxokmCTNAJ8A3gycDrxN0ulD1NLAIeD9tl8FvA7444JrXe0yYO/QRTRwFXCz7VcCr6bgmiXtAN7D6IGNZwAzwCXDVrU1DdVjOgt4yPbDtg8A1wMXD1TLumw/Yfvu8Z9/zOgXZ8ewVa1P0hzwVuCaoWtZj6QTgTcAnwSwfcD2D4ataqJjgOdJOgY4Hnh84Hq2pKGCaQfw6KrXSxT+yw4g6VTgTOCOYSuZ6EpgN1DmDs3/91LgKeBT42HnNZJmhy7qaGw/BnwUeAR4Avih7VuHrWprGiqYtMZ7Ra9bkHQC8HngvbZ/NHQ9RyPpQuBJ23cNXUsDxwCvAa62fSawHyh5vvFkRj3704AXA7OSLh22qq1pqGBaAnauej1HwV1iSc9hFErX2b5x6HomOAe4SNJ3GA2Rz5X0mWFLOqolYMn24R7oDYyCqlTnA9+2/ZTtg8CNwOsHrmlLGiqY7gReLuk0Sc9lNIF400C1rEuSGM2B7LX9saHrmcT25bbnxs+IvwS4zXaR/6rb/i7wqKRXjN86D7h/wJImeQR4naTjx38vzqPgyfppNsj9mGwfkvQu4BZGVzautX3fELU0cA7wduDrku4dv/ch23sGrGkreTdw3fgfqIeBdwxcz1HZvkPSDcDdjK7W3kO2p3QiW1IiojhZ+R0RxUkwRURxEkwRUZwEU0QUJ8EUEcVJMEVEcRJMEVGc/wMp3OIjbF7lAwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(Cinv)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x106babf98>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD8CAYAAACSCdTiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl01VWW6PHvzkSYQWYJGRgjKAgVEUUcIo5YattWKcYquqxXPMsqpXpVr64B13q93lvY9frZ/aq6u1orZWtXa0qfs7blhBZOoGiYAwEBJRBECVEQmUKS8/7YiSGQ+f7uPb/fvfuzFusmN8nvt7m52ffc89tnH3HOYYwxJjml+Q7AGGNM/FiSN8aYJGZJ3hhjkpgleWOMSWKW5I0xJolZkjfGmCRmSd4YY5KYJXljjEliluSNMSaJZfg46dChQ11+fr6PUxtjTGStWrVqn3NuWHd+xkuSz8/Pp7y83MepjTEmskSkqrs/Y9M1xhiTxCzJG2NMErMkb4wxSSyQOXkR+WvgvwEO2AB8zzl3tDvHOH78ONXV1Rw92q0f8yI7O5ucnBwyMzN9h2KMMR2KOcmLyGjgLmCyc+6IiDwO3Az8R3eOU11dTf/+/cnPz0dEYg0rbpxz1NbWUl1dTUFBge9wjDGmQ0FN12QAvUUkA+gDfNLdAxw9epQhQ4aEOsEDiAhDhgyJxDsOY4yJOck753YD9wI7gT3AAefcqz05VtgTfLOoxGmMMUFM1wwGrgMKgP3AEyJyq3PukZO+byGwECA3NzfW0xpjktmePbBtGzgH/frBqFEwYgSkWa1IdwVx4XUu8LFzrgZARJ4GzgdaJXnnXClQClBUVGQbyxpjWtTVQXk5TJoEQ4bA/v3w+uutv2fAALjlFhg50k+MERVEkt8JzBKRPsAR4FIg7stZy8pg8WLYuRNyc2HJEigpifdZjTGBcg4qKuDVV+HgQcjI0CQ/caL+gYvo/bt2wUcfwbCmFf0HDmjSt6nTTgUxJ78SeBJYjZZPptE0Yo+XsjJYuBCqqvQ5UlWln5eVxXbcG2+8kaysLGpqagC46667EBE2b94cQNTGmFaOHYOnn4annoKBA+G222DmTP1aejpkZmrSHzwYpk6F66/X+48cgd//Hh5/XN8BmA4FMsHlnPsfzrlC59yZzrnvOOeOBXHc9ixeDIcPt77v8GG9PxZ33HEHx48f55FHHsE5x7PPPss555xDYWFhbAc2xpxq2TLYuBGKizXBd/VaXXY2nH8+bN4MDz2kI33Trkhexdi5s3v3d1VxcTGFhYU89NBDvP/+++zatYvvfve7sR3UGNO24mJYsAAuvLB7F1RFNMnPnw+1tfCf/wmHDsUvzoiLZJJv7wU/iKKd22+/nQ0bNnD33XeTmZnJ/PnzYz+oMUZ99RU88wwcPQpZWZCX1/NjTZyoF2L374c33ggsxGQTySS/ZAn06dP6vj599P5YLViwgD59+vDaa68xb948hgwZEvtBjTE6B//ww1BZCZ9/Hswx8/P13cDllwdzvCQUySRfUgKlpToIENHb0tJgqmsGDRrEzTffDGBTNcYExTm9yFpTAzfdBKefHtyxc3L0Iu2RI7B+fXDHTRJeNg0JQklJfEomly1bxtatWxk5ciTz5s0L/gTGpKI334QtW+Cqq2DcuPicY8UKeOcdLa20nee+FsmRfDwVFxezefNm7r//frKysnyHY0z0HTsGa9bAtGktJZLxMGeOlls2z/kbIMIj+XhxzhbjGhOoXr3g9tu15j2ei5eysuCGG+DBB+Hll7Wu3thI3hgTRxs3QkMD9O6t8+bxlpMDF1wAa9fCxx/H/3wRYEneGBMflZXwxBOwalVizztnDpxzDpx2WmLPG1I2XWOMCd6xY/DSS9pM7BvfSOy5MzPBiia+ZiN5Y0zw3nxT2w1cc432m/Ghtlb726T4alhL8saYYH32Gbz3HsyYoXPkvjin/W3eestfDCFgSd4YE7yxY+HSS/3GMHQoTJ+ufeq/+MJvLB5ZkjfGBGvECLj11lN7j/hw0UVatrlsme9IvLEk34YNGzYwcuRIKioqfIdiTHQ0NmqjsDDNgQ8YAOeeCxs26Bx9CrIk34Z77rmHFStWcM899/gOxZjoWLtWk3ysPb+Ddt552s44DO8sPIh2CeXFF+ttwG1GH330UQD++Mc/BnpcY5LW8eM6JTJmDIRtk51+/eCSS3xH4Y2N5I0xsVu5UksmL7ssvPuubtmicaaY6Cb5sjIt03rzTe04F+sGr+hc/OzZs7/+fPXq1RQXF8d8XGOSWl2ddoAcPz6YnXvipbISXnstXNcMEiCaSb55J+9jTVvJBrST95QpU9i+fTsNDQ0A/PSnP+Xee++NNVpjkltdnSb4iy7yHUnHZs/WaaVEt1nwLJpJPk47eaelpTFlyhQ2btzIU089RW5uLjNmzIjpmMYkvX79tPvjmDG+I+nYsGHay/6DD7RpWooIJMmLyCAReVJENotIpYicF8Rx2xWvnbyBWbNmsXz5cv7u7/7OqmuM6czWrbrCNSrOPVevHVRW+o4kYYIayf8GeNk5VwhMA+L7CMZxJ+9Zs2Zx99138xd/8ReMHj065uMZk7Tq6+H557V3e1RMmKBTS7766XgQc5IXkQHAhcC/Azjn6pxz+2M9bofiuJN3YWEhvXr14mc/+1nMxzImqa1bp6PiOXN8R9J1Iroa94wzfEeSMEGM5McCNcBDIrJGRB4Qkb4BHLd9zTt59+qlnwe4k/dvfvMb/v7v/56+feP7XzAm0pyDd9+FUaOgoMB3NN1XVxe+RVtxEkSSzwBmAPc556YDh4Cfn/xNIrJQRMpFpLympib2s5aUwKxZekV/x46YE/z27dspLCzkyJEjLFiwIPb4jElmW7fCvn26mjSsdfEdWboUHn64pUIviQWx4rUaqHbONa8yeJI2krxzrhQoBSgqKgpmI9UAV7qOGzeOzZs3B3Y8Y5LawYNarTJliu9IembaNK2yqahI/KYmCRbzSN459ymwS0QmNd11KbAp1uMaY0LsG9+AO+6I7gXM0aNh+HBYvdp3JHEXVHXNnUCZiKwHzgas9tCYZLV3r87JR3GappmIbmqye3e0SkB7IJAk75xb65wrcs5Ndc5d75xL3Q79xiSzL7+E+++H5ct9RxK7qVP1ncjWrb4jiatod6E0xiTWqlU6io/qXPyJ+vSBO++EQYN8RxJXoWpr4Fww12PjLSpxGhOohgZN8hMmwODBvqMJRpIneAhRks/Ozqa2tjb0CdQ5R21tLdnZ2b5DMSaxKivhq6/gnHN8RxKsV16B557zHUXchGa6Jicnh+rqagKpoY+z7OxscnzuQm+MDxUVOoIfP953JMFqaNDtAa+8smWBZRIJTZLPzMykIIor54xJFd/6FnzxRbSratoydSq8/76+Uzn7bN/RBC400zXGmBBzTitRhg71HUnwRo+G006D9et9RxIXluSNMR07dgzuuw8+/NB3JPEhAmedBR9/rCt5k0xopmuMMSG1bp0ugErmpn1Tp2rTsiRkSd4Y0z7ndOn/qFE6rZGshgyBK67wHUVc2HSNMaZ9e/bAp59qC4Bk19ioUzZffeU7kkBZkjfGtG/NGsjI0DnrZLd/P/zhD1pOmUQsyRtj2ldYCHPnQios/jvtNBg5EjYlVxNdS/LGmPaNG6eb86SKyZNh1y5txJYkLMkbY9q2apUufkolkyfrbWWl3zgCZEneGHOq2lr4r//SVgapZOhQGDECtm/3HUlgrITSGHOqNWt0kVASLvPv1Pz5MGCA7ygCY0neGNNaY6MugJo4Efr39x1N4iVZ+2GbrjHGtLZtmy7vnz7ddyT+LF8OL7zgO4pAWJI3xrRWW6vTFRMm+I7En4MHYe3apGh1YEneGNPaeefBokXadTJVTZoE9fVJcQHWkrwxpkV9vd6mcoIHyM2F3r1hyxbfkcQssCQvIukiskZEkmMiy5hU9PDD8MwzvqPwLz1dp6s+/FAvREdYkNU1i4BKIHlqj4xJJV98AVVVybe9X09NmaK3x47pqD6iAhnJi0gOMA94IIjjGWM8aN4ZaepUv3GExaRJcMMNkU7wENx0za+BvwXafV8jIgtFpFxEyqOwWbcxKcU5rY3Pz4eBA31HEx7OabVRhMWc5EXkGmCvc25VR9/nnCt1zhU554qGDRsW62mNMUHavRs+/xymTfMdSbisWwf/8i+wb5/vSHosiJH8bOBaEdkBPAYUi8gjARzXGJMoQ4fCvHktDbqMKijQ261b/cYRg5iTvHPuF865HOdcPnAz8Gfn3K0xR2aMSZzsbDjnHOjVy3ck4TJwIAwbltpJ3hgTcVVVUF7eUiNvWhs/Xh+jiK5+DTTJO+fecM5dE+QxjTFx9t578MYbkGZjvjZNmAANDbr/awTZb9WYVHbkiC74OfNMS/Ltyc2Fb30L8vJ8R9Ij1mrYmFRWWamj1FTYqLunMjJaFkZFkL10G5PKKip0A+vTT/cdSbgdOqTthyO4HaIleWNSVUMDHD+uo3gR39GE2/HjsHRpJBuW2XSNMakqPR2+/31d1Wk6NmiQriXYuhVmzfIdTbfYSN6YVNVcEmij+K6ZMEFLKY8f9x1Jt1iSNyYV7d8P//APsGmT70iiY/x4XUsQsVJKS/LGpKKKCk1YdsG16/LydGVwxBqW2Zy8MalowwYYM0bnmk3XZGTA3/yN3kaIjeSNSTU1NfDZZ7oAynRPxBI8WJI3JvVs2KAXWyO8wMebw4fhwQe1BXFERO9lyRgTm2nTdAFUv36+I4me3r217/62bZHpvW9J3phUM2SI/jPdJwJjx8L27bq+IALlpzZdY0wqWb9eG5KZnhs7VtscfPaZ70i6xJK8MamisVGX5q9e7TuSaBs7Vm8/+shvHF1kSd6YVLFzJxw8aFU1sRowAGbMgMGDfUfSJaFL8mVlumF8WprelpX5jsiYJLFhA2RlwcSJviOJvmuvhTPO8B1Fl4TqwmtZGSxcqFVKoG0iFi7Uj0tK/MVlTOQ1NGgLg0mTNNGb2B0+rBdf+/b1HUmHQjWSX7y4JcE3O3xY7zfGxGD/ft2k26ZqgnH8OPzjP+rWiSEXqiS/c2f37jfGdNGQIbBokU3VBCUzE0aPjsTF15iTvIiMEZFlIlIpIhtFZFFPj5Wb2737jTFd0NCg/0QiUdcdGWPHwief6D65IRbESL4e+Klz7gxgFvAjEZnckwMtWQJ9+rS+r08fvd8Y00ObN+vUQsS6J4be2LE6J79jh+9IOhRzknfO7XHOrW76+CBQCYzuybFKSqC0VDt6iuhtaalddDUmJhUVugtUREr+ImP0aJ22CXl/+UCra0QkH5gOrOzpMUpKLKkbE5ijR3XLuqIirUs2wUlPhxtvhGHDfEfSocB+6yLSD3gK+Ilz7ss2vr5QRMpFpLympiao00aarQkwcbd5s24OYlU18TFpkjZ7C7FAkryIZKIJvsw593Rb3+OcK3XOFTnnioaF/JUvER578DC/+ME+6qo+YYTbw8GqWn7yg0OUPWKbKpsAbdig0zSjezSDajpTXw9r18Lu3b4jaVfM0zUiIsC/A5XOuX+KPaQkc+QI7NoF1dWwZw/ccguIsPTnr3PbkVWtvrXhSDqLF99Nya3oPOrRo1BQYB0DTc9deKE+B62qJj5E4MUXte1wSF9Ig5iTnw18B9ggImub7vulc+7FAI4dXdu3w5tvaoJ3Tudkhg+HY8cgO5tXa6YziHzqyEJwZFFHGo3s3NX0x7h2rfasBp3zO+MM7Zdh27WZ7sjL8x1BcktP18c4xBdfY07yzrl3ABsmHDsGGzfqL3zIEO34V1enI6mxY3XD5MzMr789PS+HiqqcUw6T17wmoKQEvvhCE31lJbz9tr4TsKvSpqvee08v9owc6TuS5Jafrxe3v/oqlBuxhKp3TSTt3w8rV2r71mPH4NJLYc4cmDBB/7VjyZLWfXrgpDUBInpBZ+ZM/XfggL5oNJ/zhRfg4osh59QXCmM4cABefhmKiy3Jx1t+vt7u2BHKC9yW5HvKOXjuOd2EAfSXO3Nml+flmgfkixdr24bcXE3w7Q7UBw5s+fjzz3VU/8ADMH06zJ0b+iZJJsEqKvQ2hEkn6YwapX2BQro5ujiX+GqOoqIiV15envDzBuLLL7WfNMArr+iI+9xzWyfhRKirg7feghUr9Ak2b14on2DGk9/9Tq8D/eAHviNJDUePQnZ23E8jIqucc0Xd+RlbHdFVBw7oyP3Xv9aLqQBXXAGXX574BA/aLnbuXPjhD/UaQFOjJKu9N+zbp+/07EU/cRKQ4HvKpms6U18P776ro+bGRp2SCdPih2HD4LbboKGBsjL45Q9qOHLE4Rhu/fhT1b590Ls3TJniO5LUceQIPP88TJ0aus1ELMl3xDl46CFd6FBYqCP3MPb/SEuDtDQWL4biIy8wij08zQ1sofDrfvyW5FNIYaFe9E9P9x1J6ujVS8soe/e2JB8JBw9qKZSIjtz79OmwUiYsdu6Ep/hLbuL/cTOPsZTLWMFs68efSurrNblbgk+s5vnRENbL25z8ierrtR79n/+5pWpm2rRIJHjQCp2DDOA/+CsqOJPLWMpclpI7xlolpIxly+D++7V/vEms/Hxd23LggO9IWrEk32zbNrjvPnj9dRg/PpIrBZv78deTydPcwAecw/isXdzzv+wPPiU4p6WTAwbYSN6HE+vlQ8SmawBeekkXNA0ZArfeqkk+glrX3qexcczVlPzPem75bobuSZmRYT1Mkll1tY4ii4t9R5KaRozQXlMZ4Uqr4Yomkerr9TYjA8aN0zn4884L3S+ou1r34xcgU9+6P/KIvoh985uW6JPVhg36/C0s9B1JahKBBQt8R3GK1Jyu2boV/u3fYPly/XziRG1FEPEE367mJkqrV2vHPA8L4EycNTZq76SJE7XSw/jT0NAyiAyB1EryX3wBjz6qK4TS0mDMGN8RJc4ll8D558MHH2h3TBMpnS5ycw6uukrfjRp/amvhV7/SpoIhkaRD1zasWQN/+pP+lcydq38MqXRxSgQuu0wXbbzxhrYsPvts31GZLigra93Mrs1FbunptsI1DAYP1hyzcyecdZbvaIBkT/LO6dumzEztxDd5sia6/v19R+aHCFxzjd6GdIMDc6rFi1t3KwVaL3Krr9fCgalTU/e5HRbNMwRVVb4j+VryTtfs3g1/+IOO3kE7xd1wg/0RpKfDtddqOwTnTs0eJnTaW8z29f3btsHSpdoF0fiXlwd794bmbyv5knxtLTz+OPz+91BTY/3WO/L66/o4heTJaNqWm9vJ/RUVukCioCBhMZkONK+xCclS8+RK8uvXw29/qyObiy+Gu+6Com515UwthYXaOvnJJ7U6w4RS8yK3E329wUxdHWzZolORqXSNKcxOP13XKgwb5jsSIBmSfG2tvjUCLTuYORMWLdIkb6VkHcvJ0Tn6jz6C117zHY1pR0kJlJbqAFFEb0tLm+bjKyt1odvUqb7DNM0yMnTbzyFDfEcCRPXCq3P6VmjlSn2Sjx+vz/gBA+DKK31HFy3Tp2vv8RUr9LpFSCoCTGutF7md4PPPtfV1KpUDR0FdnV58zcvTvR88il6S37BBm4jt3auN+mfPhlmzfEcVbVdcocvhPT8ZTQ9ccoku5LNVzOFSXa21ryFokxJIkheRK4HfAOnAA865XwVxXEC31frwQ+3RnJmpySgjA667TuuCMzMDO1XKSk+H+fN9R2G6q75e/xaSdaV2lOXkaDllVVX0k7yIpAO/BS4DqoEPROR559ymHh3w+HF9Fayqgu3btRSysRFuukkT/fnnwwUXxBq2ac+77+oLqU17hd+DD2oyufpq35GYk2Vl6fRnCOrlgxgCzAS2Oec+AhCRx4DrgPaT/NGj2mfj+HE4dEiTytixWu1x4IDWt4voVeoLLtB+7s2lkGnRv1Ycal9+Ce+9p/V5kyf7jsa0p6YGPvnELriGWV6eXjdsfsflSRBnHg3sOuHzauDck79JRBYCCwHOGDgQnnii5YvZ2XrRFFra/Y4erVtpmcSaO1cvaj//vL7IDhrkOyLTlvXrdSBkrQzCKy9PCxp27/a6P0UQSb6tKz6ntDl0zpUCpQBF06Y57rhD59N7926907mI9zmslJaeDjfeqLsLPfustk61i3rh4pwWIDS3yDbhVFAAt98Ow4d7DSOIuY9q4MT6rRzgkw5/IjNT/+ODB7dO8Cah2u1sOHiwzslXVYVm1Z45wc6dsH+/TdWEXVaW9szyPMUcxEj+A2CCiBQAu4GbgVsCOK6Jo047G559tk6ZeR6FmDYMHaovwrY5SPjt3g3r1unvy1Oyj/mszrl64MfAK0Al8LhzbmOsxzXx1VFnQ0CnaJoT/K5dtjF0mPTtq2tDbF1D+O3fD++/rwsOPQnkpcU596JzbqJzbpxzbkkQxzTx1Wlnw2Y1NVqqZxuNhMPOnbB2rb3oRkVzFzmPpZRWj5iiOu1s2GzYMJ26efttXb9g/FqxQvsM2cXwaOjfXysGLcmbROuws+HJrrxSS1yfey5Ue1emnEOHdPX31KneL+aZbsjL03dgnvZWtmdKiuqws+HJevXSbpU1NTqiN36sX6+rv23bxmjJzdXrJwcPejm9Nb1IYe12NmzLhAnam79v37jGZNrhnM7FW8VT9Eyb5vWF2ZK86bprrvEdQeo6elTfctkoPno8Xz+x6RrTPc2rLT/4wHckqaV3b109aTudRdPKlbrVpod5eUvypvs2bYJXXtFduUz8NTTAsWP6sVXVRFNami6M2r8/8adO+BlNtIloa9uMDG1i5qliIKVs2QL33tuyzaWJnubaZA9tQizJmy5p1efmrP682HCF1v6uWeM7tOS3Zo1O1wwd6jsS01PDh+vv0EO9vCV506nmPjdVVTpwr6qCby05m3c/yYNXX4UjR3yHmLwOHoRt27RCw2rjo0tE9+H1kOStusZ0qs0+N0eEO176JmuW7rO+//G0bp2+slpVTfRNmaIbvTQ2JvQF25K86VR704jrdg+FwqYphAQ/cVOCczpVk5urS+NNtE2bpv8SzP4qTac67XOzahX87nfW8iBoIrq38eWX+47EBKWxMeErXy3Jm0512udm0CD47DNreRAPw4e37G9sou/RR+GPf0zoKS3Jm0512udm3DhtmvXOO7Bvn9dYk8ahQ/D00/Z4JpvTT4dPP21Z95AAluRNl5SUwI4d+m5zx442et5cfrlu6/jCC1Y7H4Q1a7QhmT2WySUvT3+nu3Yl7JSW5E0w+vWDyy7TErHdu31HE23O6XWOvDzt52+SR06OFigksJTSqmtMcGbM0Fpg65IYm48/hi++gOJi35GYoGVlwahRCV35akneBOfEfWFra63sr6fKy/XK9hln+I7ExMMll0B6esJOZ9M1JniVlfCv/6qT96b7hg6F887T/kAm+YwfDwUFCTtdTEleRP6PiGwWkfUi8oyIDAoqMBNh48drWeULL1jtfE8UF8OcOb6jMPH08ccJm5ePdSS/FDjTOTcV+BD4RewhmcjLzIR587T8b/ly39FER0MDbN1qFTWp4MUXE7auJKYk75x71TnXPFR7D7BVG0aNHw9nnqlPZOs73zWbNmk3uI8+8h2Jibe8PC2jbGyM+6mCnJO/DXgpwOOZqLviCr2AaH3Qu2blSr1YPXas70hMvOXm6oKozz6L+6k6vbIjIq8BI9v40mLn3HNN37MYqAfKOjjOQmAhQG57zVBMcunfHxYtSmglQWTt3g3V1XDVVbb7UyrIy9PbnTu1pDKOOk3yzrm5HX1dRBYA1wCXOtf+ZKJzrhQoBSgqKrJJx1SRnq5zzBUVOoVjbYnbtnIl9OplLYVTxcCB+q+6Gs49N66niqlGS0SuBH4GXOScO9zZ95sU9fnn8MwzmsCuvdZ3NOHT0KB/7NOna6I3qeGv/koTfZzFWoj7r0AvYKnoW8z3nHO3xxyVSS5Dhmjd9/Ll2k+7+a2qUenp8OMfw/HjviMxiTR4cEJOE2t1zXjn3Bjn3NlN/yzBm7ZddFFL7XxDg+9owqO+Xv+lpdkoPtUcPw4vvaQbtceRrXg1iZGVBVdfDTU1sGKF72jCY/Vq+PWvE76RhAmBjAzYsEFXiMeRJXmTOBMn6kUma2CmGhr0BW/wYK1EMqlFREsp47zy1ZK8SayrroJJk3xHEQ4bN8L+/XDBBb4jMb7k5WnH0Ti+k7MkbxKvoQHeektXeKYq5/RC9LBh+g7HpKbmNUNxHM1bmzuTeCKweTN8+aWu7szO9h1R4lVV6WrH66+3xU+pbNQorT6rq4vbKWwkbxIvLQ2++U3dx/SVV3xH40deHixYAGed5TsS41NaGtx5p264E69TxO3IxnRk1CiYPVv3Mt22zXc0ieWcjt4LCqzlg2kRp+6jluSNPxdfrHPSf/pT6tTOOwcPPwzvvOM7EhMW+/ZpGe3WrXE5vM3JG38yMuCGGzTxpcqIdvt2bSU8ebLvSExYDByo1TVVVXG5CG8jeePXqFFw+un68bFjp3y5rAzy83XqMj9fP48s52DZMv2jnj7ddzQmLDIz9W8gTpt7W5I34fDWW3Dffa0SfVkZLFyoAxzn9Hbhwggn+q1btaXwhRemzjsX0zW5ufDJJ3HpX2RJ3oRDQQEcOABLl3591+LFcPik3qaHD+v9keMc/PnPcNpp1k7YnCovT69L7d4d+KFtTt6Ew5gx2qlyxQooLITx49t99xqnd7XxJQJ/+Zf6KmWjeHOy3Fz4xjfist+CjeRNeBQXa1+bZ5+FQ4dobwOxyG0s1lwaN2yYtVk2bcvO1rUjI0YEfmhL8iY8MjJ0tNvYCJ9+ypIlukXsifr0gSVL/ITXY2+9BU88kTploqZnnNNV0AFv7m1J3oTLiBHwk5/AuHGUlEBpqQ5+RfS2tFS/LTIVN19+qT1qUqlM1PTMpk1afLBnT6CHtSRvwicrS5PimjWUXF7Djh06uNmxQ78cqYqbV1/V4Od2uFWyMS3zkAFfdLIkb8Lp6FGttHnyyVZlZZGquNm+XTcwnzNHq2qM6Uj//rq3QMAdKS3Jm3Dq3VtXw+7dq1sGNl28jFTFzbJlmtxnz/YdiYmKvDx9Mgdd/1lLAAAJNklEQVTYx8aSvAmv8eN1b9h162DVKqD9ypowVNycvDr3MZkP3/62XlA2pityc/Wt6b59gR3SkrwJt4su0mT/8stw8GBoK25OXJ3b1x1kZ1Uj37+rL2Wvj/QbmImWiRPhllu09UVAAknyIvI3IuJEZGgQxzPmayI6bfPtb0P//u1W3JSU+A2z+VpBOvV8h4e5kSfDe63AhFe/fpros7ICO2TM7yNFZAxwGRDGWVGTDPr0aenOV11NybdHUFKS6TemkzRfE7iYNxjOXl5jbqv7jemyvXu1U+msWYEcLoiR/P8F/haIT8d7Y5odOAAPPQTPPRe3DRZ6KjcXxrKd2SxnDdPZysSv7zemWz7+WKcnDxwI5HAxJXkRuRbY7ZxbF0g0xnRk4EBtfVBRoc2+QuR/L/6S+ZlPUcMwXuIqIBzXCkwEBby5d6fTNSLyGtDW1aPFwC+By7tyIhFZCCwEyLXhjemp88+Hzz+Ht9/W+ctzz/UdEQA3XXeUgncH8N+X3kj97izycjXB+75WYCJoxAjo1Uvn+qZOjflw4nr4tldEzgJeB5qXpuQAnwAznXOfdvSzRUVFrry8vEfnNYbGRu0FU1kJt93md06k+e9HpGXvVmNi9cgjOl3zox+1ultEVjnnirpzqB5P1zjnNjjnhjvn8p1z+UA1MKOzBG9MzNLStJHZtddqi2Kfli/XxVqNjZbgTXDy8mD//jZ3S+suq5M30ZSRATNmaGLdtw9Wr058DGvWwGuvQV2dJXgTrJkz4Wc/02mbGAW2FK9pNG9M4r37rq6I/eor7ROTiIS7di08/7wu1LruOkvyJlgBJPdmtt7aRN/VV0N9vVbc1NbCNdfo5sjxsmqVTtEUFMBNN1nbAhMf5eW6HeB118V0GJuuMdGXng7XXw+XXKJ9bh588NRWlUEaNAjOOAPmz4/vi4lJbQcO6PO5ri6mw1iSN8lBRPvc3HKLdrDMzg72+AcO6BQNwLhx2mbBEryJp7w8vaBfXR3TYex9pkkuEyfChAma9A8f1mmVCy6A00/v2fEaG3V65vXXW45/coc0Y+JhzBh9HldVwdixPT6MJXmTfJovgn72mW4ntWmTTq+ce25LZ7PONDToz739tvYSKSjQjZYtwZtE6dULRo6MuQGSJXmTvAoKYNEirWX/4ANdPDV0KPzwhzqPf/iw/iGlp+vuUwcP6vznyJFw5Ag8+6zu1HPTTVBYaBU0JvEmTdJighhYkjfJrVcv7XczZ472vKmpadlQ+7HHdJTUvFoVdKT/ve9py4SFC2H4cEvuxp+LL475EJbkTWrIzITp01vfN3OmXkRtaND+3X37wqhRLV8fMSKxMRrTnoaGlsFJN1mSN6nrzDN9R2BM58rKtJXH/Pk9+nEroTSRc/JeqmVlviMyJo76949pc29L8iZSTtxL1Tm9XbgwgERvrxwmrPLytBBg794e/bgleRMpzXupnijmvVTj9sphTACaW2n3sJTSkryJlPae5zGVEsfllcOYgAwaBAMG9HinKLvwaiIlN7ft53pM+4bE5ZXDmICIwIUX9nghno3kTaQsWXLqcz3mvVTbe4WwbSpNWBQVweTJPfpRS/ImUkpKoLS0pTtBXp5+HtNeqnF55TAmQM7piuwesOkaEzklJQFvkN18sMWLdYom13bhNiHUw+kaS/LGQBxeOYwJkEiPV7zadI0xxiQxS/LGGJPELMkbY0wSiznJi8idIrJFRDaKyD8EEZQxxphgxHThVUQuAa4DpjrnjonI8GDCMsYYE4RYR/I/BH7lnDsG4JzrWQcdY4wxcRFrkp8IzBGRlSLypoic0943ishCESkXkfKampoYT2uMMaYrOp2uEZHXgJFtfGlx088PBmYB5wCPi8hY505tfOycKwVKAYqKinrWGNkYY0y3SBv5uOs/LPIyOl3zRtPn24FZzrkOh+oichDY0uMTJ5ehwD7fQYSEPRYt7LFoYY9Fi0nOuf7d+YFYV7w+CxQDb4jIRCCLrv0ytjjnimI8d1IQkXJ7LJQ9Fi3ssWhhj0ULESnv7s/EmuQfBB4UkQqgDljQ1lSNMcYYP2JK8s65OuDWgGIxxhgTMF8rXks9nTeM7LFoYY9FC3ssWthj0aLbj0VMF16NMcaEm/WuMcaYJJbQJC8iVzb1udkmIj9P5LnDRETGiMgyEals6vmzyHdMvolIuoisEZEXfMfik4gMEpEnRWRz0/PjPN8x+SIif93091EhIo+KSLbvmBJJRB4Ukb1NhS3N950mIktFZGvT7eDOjpOwJC8i6cBvgauAycB8EenZpoXRVw/81Dl3BrqQ7Ecp/Fg0WwRU+g4iBH4DvOycKwSmkaKPiYiMBu4CipxzZwLpwM1+o0q4/wCuPOm+nwOvO+cmAK83fd6hRI7kZwLbnHMfNVXlPIY2N0s5zrk9zrnVTR8fRP+QR/uNyh8RyQHmAQ/4jsUnERkAXAj8O2j1mnNuv9+ovMoAeotIBtAH+MRzPAnlnHsL+Pyku68D/tD08R+A6zs7TiKT/Ghg1wmfV5PCia2ZiOQD04GVfiPx6tfA3wKNvgPxbCxQAzzUNHX1gIj09R2UD8653cC9wE5gD3DAOfeq36hCYYRzbg/oYBHotPNvIpO8tHFfSpf2iEg/4CngJ865L33H44OIXAPsdc6t8h1LCGQAM4D7nHPTgUN04e14Mmqaa74OKABOB/qKiK3J6YFEJvlqYMwJn+eQYm+/TiQimWiCL3POPe07Ho9mA9eKyA50Cq9YRB7xG5I31UC1c675Xd2TaNJPRXOBj51zNc6548DTwPmeYwqDz0RkFEDTbaft3ROZ5D8AJohIgYhkoRdRnk/g+UNDRASdd610zv2T73h8cs79wjmX45zLR58Tf3bOpeSIzTn3KbBLRCY13XUpsMljSD7tBGaJSJ+mv5dLSdGL0Cd5HljQ9PEC4LnOfiDW3jVd5pyrF5EfA6+gV8ofdM5tTNT5Q2Y28B1gg4isbbrvl865Fz3GZMLhTqCsaSD0EfA9z/F44ZxbKSJPAqvRarQ1pNjKVxF5FLgYGCoi1cD/AH6FtnT/PvpC+K1Oj2MrXo0xJnnZildjjEliluSNMSaJWZI3xpgkZkneGGOSmCV5Y4xJYpbkjTEmiVmSN8aYJGZJ3hhjktj/B0ZmVmBfRpuWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "lamb = 0.01\n",
    "Cinv = np.linalg.inv(C + lamb*np.eye(C.shape[0],C.shape[1]))\n",
    "mu_p = b.dot(Cinv).dot(y_meas)\n",
    "sig_p = a - b.dot(Cinv).dot(b)\n",
    "plt.plot(x_sample,y_meas,\"bo\",label=\"$\\mathbf{y}$\")\n",
    "plt.errorbar(x_p,mu_p,yerr=sig_p,fmt=\"ro\",label=\"$\\hat{y}$\")\n",
    "plt.legend()\n",
    "plt.xlim(0,10)\n",
    "ylim=plt.ylim()\n",
    "plt.plot(X_space,y_real,\"r--\",alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x1a19468748>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASwAAAD8CAYAAADNNJnuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEqBJREFUeJzt3X2MXXWdx/H3pzOtlELDQ1eBttKabdCuWcXMIkJidq3JViHgH5rg07IuSf9YC1XZKOgf/rNmn4wPybqsDUhIbGRJwUAMy8Mi/GEilQKNWgaWpi7tSLEUn6A8tDPz2T/u7TqrnbmnvefMub/ezys56T33Xr7nO5c73/k9nXNkm4iIEixoO4GIiKpSsCKiGClYEVGMFKyIKEYKVkQUIwUrIoqRghURxUjBiohipGBFRDFGmwi67IwRn7uy/tDjz72+9pgAmqo/5sgLB+sPCmjhwkbiHjpzUSNxNd1IWKwmgjYQk2Y+g0Mv/pLJVw/29Sn85V8s8Qu/rPblf/THr91re30/x6tDIwXr3JWjPHzPitrjjv3DxtpjAowerP+besbNP6w9JsDoWcsbifvMR89tJO7oy42EZbqBut3EHy6A0Vfq/349dftX+47xwi+n+NG9b6z03pGzn17W9wFr0EjBiojBZ2CahprADUnBihhSxhx2Q83KhqRgRQyxtLAiogjGTBV2eaksa4gYYtO40taLpNMkbZX0pKRxSe+SdIak+yU93f339H7zTcGKGFIGpnClrYKvA/fYfjPwNmAcuA54wPYa4IHufl8qFSxJ6yU9JWmXpL4PGhGDoY4WlqSlwLuBmwBsH7L9a+By4Jbu224BPtBvvj0LlqQR4BvA+4C1wIclre33wBHRLgOH7UpbD28CngdulvS4pBslLQHeYHsfQPffvld+V2lhXQDssr3b9iHgVjqVMyIK5ordwW6XcJmk7TO2DTNCjQLvAG6wfT5wkBq6f0dTZZZwObB3xv4E8M7ff1P3B9gA8MblmXyMGHiGqeqThAdsj83y2gQwYXtbd38rnYL1C0ln294n6Wxgf1/5Uq2FdbTzlf7gx7S92faY7bFlZ2YsP2LQdVa6V9vmjGM/B+yVdF73qXXAE8BdwJXd564E7uw35ypNoQlg5Yz9FcCz/R44Itompo7aHjkuVwNbJC0CdgOfoNMguk3SVcAe4EP9HqRKwXoEWCNpNfBz4ArgI/0eOCLa1Rl0r6dg2d4BHK3LuK6WA3T1LFi2JyVtBO4FRoBv2d5ZZxIRMf8667CauE5PcyqNjtu+G7i74VwiYp5NN3JhseZkOi9iSJ2wLayIOPEYMVXY2XkpWBFDLF3CiCiCEYc80nYaxyQFK2JIdRaOpkvI+HOvb+SGEdMN/TGYXFJ/s3j/xotqjwkN3S0GWPhSM3GbuhPNSEHXnZtcXP//NNdUZzLoHhFFsMVUXZVvnqRgRQyx6bSwIqIEnUH3skpAWdlGRG0y6B4RRZnKOqyIKEFWukdEUaYzSxgRJeic/JyCFREFMOJwTs2JiBLYZOFoRJRCWTgaEWUwaWFFREEy6B4RRTDKBfwiogyd23yVVQLKyjYialTrjVTnRQpWxJAyWekeEQVJCysiimArLayIKENn0D2n5kREEXJNdwA0BaMH67+tSRN3twE4tLT+mCc/18xtXV4+q5nPYOS1RsKy4HAzcZv4fk0vqj0kAAum6/9/php+/M6ge8awIqIQWekeEUXISveIKEpuQhERRbDh8HRZBausbCOiNp0u4YJKWxWSRiQ9Lul73f3VkrZJelrSf0jqe1ojBStiiE11zyfstVW0CRifsf9PwFdtrwF+BVzVb749C5aklZIelDQuaaekTf0eNCLad2RZQ5WtF0krgEuAG7v7At4DbO2+5RbgA/3mXGUMaxK41vZjkk4FHpV0v+0n+j14RLSp1lNzvgZ8Fji1u38m8Gvbk939CWB5vwfpma3tfbYf6z5+kU6Tr+8DR0T7prvXde+1AcskbZ+xbTgSQ9KlwH7bj84IfbRmWd/LXY9pllDSKuB8YFu/B46IdnVmCSufS3jA9tgsr10MXCbp/cBJwFI6La7TJI12W1krgGf7zblye1DSKcDtwKds//Yor284Un0nXz3Yb14R0bAjC0f7HcOyfb3tFbZXAVcA37f9UeBB4IPdt10J3NlvzpUKlqSFdIrVFtt3zJL0ZttjtsdGT1rSb14RMQ+OoUt4PD4HfEbSLjpjWjf1m2/PLmF3tP8mYNz2V/o9YEQMhiZOfrb9EPBQ9/Fu4II641dpYV0MfBx4j6Qd3e39dSYREe2oc+HofOjZwrL9A44+4h8RBbPF5AAVoypyLmHEEMvVGiKiCLmAX0QUJQUrIoqQC/hFRFH6WGPVikYK1sgLBznj5h/WHnf/xotqjwnN3DDikb+/ofaYAJf8WTMrSvZ8ZFUjcfs/e+zoDi9t4Bdtuv6QAKOvNvQh9MmGycIu4JcWVsQQS5cwIoqQMayIKIpTsCKiFBl0j4gi2BnDiohiiKnMEkZEKTKGFRFFyLmEEVEOd8axSpKCFTHEMksYEUVwBt0joiTpEkZEMTJLGBFFsFOwIqIgWdYQEcXIGFZEFMGI6cwSRkQpCmtgpWBFDK0MukdEUQprYqVgRQyxtLAALVzI6FnLa4/b1Gf78ln1B27q7jbffnhrI3E/9FcbG4n7qzWvayQur9UfcnJJM1+wJr5f0zX85hqYnk7BiogSmOZaAQ1JwYoYYlmHFRHlSMGKiDIog+4RUZDCWlhlrcuPiPoYPK1K21wkrZT0oKRxSTslbeo+f4ak+yU93f339H5TTsGKGGqquM1pErjW9luAC4FPSloLXAc8YHsN8EB3vy+VC5akEUmPS/pevweNiAHhittcIex9th/rPn4RGAeWA5cDt3TfdgvwgX7TPZYW1qZuIhFxoqhesJZJ2j5j23C0cJJWAecD24A32N4HnaIGvL7fdCsNuktaAVwCfAn4TL8HjYgBcGwLRw/YHpvrDZJOAW4HPmX7t1L9M5BVW1hfAz4LTM/2BkkbjlTfQ9Ov1JJcRDTLrrb1ImkhnWK1xfYd3ad/Iens7utnA/v7zbdnwZJ0KbDf9qNzvc/2ZttjtscWLVjcb14RMR+mVW2bgzpNqZuAcdtfmfHSXcCV3cdXAnf2m26VLuHFwGWS3g+cBCyV9G3bH+v34BHRLtWzDuti4OPATyTt6D73eeAfgdskXQXsAT7U74F6Fizb1wPXA0j6c+DvUqwiTgAVZgArhbF/wOxrH9b1f4TfyUr3iKGlE/tqDbYfAh5qJJOImH+FnZqTFlbEMJt13n8wpWBFDKtcwC8iSlLTLOG8ScGKGGaFFaxcrSEiitFIC+vQmYt45qPn1h534Uu1hwRgpIE7sOz5yKr6g9Lc3W2mPv9CI3EP3Vf/3ZMAFr5Uf9PgnG/u6P2m47D36rfXHrOurly6hBFRBtPztJtBk4IVMczSwoqIUqRLGBHlSMGKiGKkYEVECeR0CSOiJJkljIhSpIUVEeVIwYqIImQMKyKKkoIVEaVQYRfwy9UaIqIYaWFFDLN0CSOiCBl0j4iipGBFRDFSsCKiBKK8WcIUrIhhlTGsiChKClZEFCMFq9MvHn25gcANfbgLDjcQtKFcf7XmdY3EberuNlPNpMvoy/VfFmXvNfXf3QbgzJ2Ttcfc+0o9X7B0CSOiHClYEVEEZ5YwIkqSFlZElKK0MaxcrSFimLni1oOk9ZKekrRL0nVNpZuCFTGsqharHgVL0gjwDeB9wFrgw5LWNpFypYIl6TRJWyU9KWlc0ruaSCYi5o/43a2+em09XADssr3b9iHgVuDyJnKuOob1deAe2x+UtAg4uYlkImJ+1TSGtRzYO2N/AnhnLZF/T8+CJWkp8G7grwG6FfRQE8lExDyrXrCWSdo+Y3+z7c3dx0dbxdvIcH6VFtabgOeBmyW9DXgU2GT7YBMJRcQ8ql5WDtgem+W1CWDljP0VwLN9ZDWrKmNYo8A7gBtsnw8cBP5gFkDSBknbJW2ffCW1LGLgVRy/qtBtfARYI2l1d8joCuCuJlKuUrAmgAnb27r7W+kUsP/H9mbbY7bHRhcvqTPHiGhKDbOEtieBjcC9wDhwm+2dTaTbs0to+zlJeyWdZ/spYB3wRBPJRMT8quvUHNt3A3fXE212VWcJrwa2dJt7u4FPNJdSRMyX0la6VypYtncAsw24RUSJKq5iHyQ5lzBimKVgRUQJjqx0L0kKVsQQ03RZFSsFK2JYZQwrIkqSLmFElCMFCyyYXlh/3JGGPtzRg/UHPry0/ru6APBaM2EXvtTMh9vE3W0AJhu4Xkgjd3oCfrO6/l+zqYfr+VzTwoqIcqRgRUQRcteciChF1mFFRFlcVsVKwYoYYmlhRUQZsnA0IkqSQfeIKEYKVkSUwWTQPSLKkUH3iChHClZElCALRyOiHHYu4BcRBSmrXqVgRQyzdAkjogwG0iWMiGKUVa9SsCKGWbqEEVGMzBJGRBlytYYug6YaidyI6UVNBG0gJjC5pJmbOpzzzR2NxN17zdsbidvEDSO+e+0/1x8U+Nvz3lt7zGde7f8D6CwcLatipYUVMcxytYaIKEVaWBFRhoxhRUQ5ci5hRJSksC7hgrYTiIiWdG+kWmXrh6R/kfSkpB9L+q6k02a8dr2kXZKekvSXvWKlYEUMM7va1p/7gbfa/lPgv4HrASStBa4A/gRYD/ybpJG5AlUqWJI+LWmnpJ9K+o6kk/pKPyIGgytu/RzCvs/2ZHf3YWBF9/HlwK22X7P9M2AXcMFcsXoWLEnLgWuAMdtvBUboVMWIKJympyttwDJJ22dsG47zkH8D/Gf38XJg74zXJrrPzarqoPsosFjSYeBk4NljTDIiBo05loWjB2yPzfaipP8CzjrKS1+wfWf3PV8AJoEtR/6zWbKaVc+CZfvnkr4M7AFeAe6zfd9REt4AbABYeOrpvcJGRMuEa1s4anvO848kXQlcCqyz/++gE8DKGW9bQY/GUJUu4el0+pqrgXOAJZI+dpSEN9sesz02snhJr7ARMQjmYdBd0nrgc8BltmeeBHkXcIWk10laDawBfjRXrCqD7u8Ffmb7eduHgTuAi44v9YgYKPMzS/ivwKnA/ZJ2SPr3zqG9E7gNeAK4B/ik7Tkvm1BlDGsPcKGkk+l0CdcB2/tIPiIGwbGNYR3/Yew/nuO1LwFfqhqryhjWNklbgcfoDJg9DmyueoCIGFzdGcBiVJoltP1F4IsN5xIR86qW7t68yrmEEcPKpGBFREHK6hGmYEUMs1zALyLKkYIVEUWwYaqsPmEjBUvTMPpK/ZV7cnEzd4xZMF1/3NFXm/nL9fJZzXwGe69u5u42Z+6c7P2m4/Cb1fV/dZu4uw3AkvtPqT3mgqvmvApLdWlhRUQxUrAioggGck33iCiDwRnDiogSmAy6R0RBMoYVEcVIwYqIMuTk54gohYET8fIyEXGCSgsrIsqQU3MiohQGZx1WRBQjK90johgZw4qIItiZJYyIgqSFFRFlMJ6a876lAycFK2JY5fIyEVGULGuIiBIYcFpYEVEE5wJ+EVGQ0gbd5QamNSU9DzxT4a3LgAO1J9CckvItKVcoK99ByPVc23/UTwBJ99D5Wao4YHt9P8erQyMFq/LBpe22x1pL4BiVlG9JuUJZ+ZaU64lmQdsJRERUlYIVEcVou2Btbvn4x6qkfEvKFcrKt6RcTyitjmFFRByLtltYERGVtVawJK2X9JSkXZKuayuPXiStlPSgpHFJOyVtajunKiSNSHpc0vfazmUukk6TtFXSk93P+F1t5zQXSZ/ufg9+Kuk7kk5qO6dh0krBkjQCfAN4H7AW+LCktW3kUsEkcK3ttwAXAp8c4Fxn2gSMt51EBV8H7rH9ZuBtDHDOkpYD1wBjtt8KjABXtJvVcGmrhXUBsMv2btuHgFuBy1vKZU6299l+rPv4RTq/UMvbzWpuklYAlwA3tp3LXCQtBd4N3ARg+5DtX7ebVU+jwGJJo8DJwLMt5zNU2ipYy4G9M/YnGPAiACBpFXA+sK3dTHr6GvBZYNBPFHsT8Dxwc7f7eqOkJW0nNRvbPwe+DOwB9gG/sX1fu1kNl7YKlo7y3EBPV0o6Bbgd+JTt37adz2wkXQrst/1o27lUMAq8A7jB9vnAQWCQxzNPp9MTWA2cAyyR9LF2sxoubRWsCWDljP0VDHDTWtJCOsVqi+072s6nh4uByyT9D52u9nskfbvdlGY1AUzYPtJi3UqngA2q9wI/s/287cPAHcBFLec0VNoqWI8AayStlrSIzsDlXS3lMidJojPGMm77K23n04vt622vsL2Kzuf6fdsD2Qqw/RywV9J53afWAU+0mFIve4ALJZ3c/V6sY4AnCU5ErVxexvakpI3AvXRmWr5le2cbuVRwMfBx4CeSdnSf+7ztu1vM6URyNbCl+4drN/CJlvOZle1tkrYCj9GZPX6crHqfV1npHhHFyEr3iChGClZEFCMFKyKKkYIVEcVIwYqIYqRgRUQxUrAiohgpWBFRjP8FL5tdU/zduj4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(Cinv)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now consider that we want to do the same but for a vector of $\\mathbf{\\hat{y}}$ values for the $\\mathbf{\\hat{x}}$ points.\n",
    "\n",
    "$$\\begin{bmatrix}\\mathbf{\\hat{Y}} \\\\ \\mathbf{Y} \\end{bmatrix} \\sim Pr\\left(\\begin{bmatrix}\\mathbf{\\hat{Y}} \\\\ \\mathbf{Y} \\end{bmatrix} = \\begin{bmatrix}\\mathbf{\\hat{y}} \\\\ \\mathbf{y} \\end{bmatrix}\\right) = \\mathcal{N}\\left(\\begin{bmatrix}\\mathbf{0} \\\\ \\mathbf{0} \\end{bmatrix}, \\begin{bmatrix} K(\\mathbf{\\hat{x}},\\mathbf{\\hat{x}}) & K(\\mathbf{\\hat{x}},\\mathbf{x}) \\\\ K(\\mathbf{x},\\mathbf{\\hat{x}}) & K(\\mathbf{x},\\mathbf{x}) + \\lambda I \\end{bmatrix}\\right) = \\mathcal{N}\\left(\\begin{bmatrix}\\mathbf{0} \\\\ \\mathbf{0} \\end{bmatrix}, \\begin{bmatrix} A & B^T \\\\  B & C \\end{bmatrix}\\right) $$\n",
    "\n",
    "Now\n",
    "$$\\mathbf{\\hat{y}} \\sim Pr(\\mathbf{\\hat{Y}}=\\mathbf{\\hat{y}}| \\mathbf{Y} = \\mathbf{y}) = \\mathcal{N}\\left(\\mathbf{\\hat{y}} ; B^T C^{-1} \\mathbf{y},  A - B^T C^{-1} B \\right)$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.atleast_2d(x_sample).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact\n",
    "X_space = np.linspace(0.1, 9.9, 1000)\n",
    "X_space = np.atleast_2d(X_space).T\n",
    "\n",
    "def gen_data(n = 10,sigma =1.0):\n",
    "    # Get the \"real\" value\n",
    "    x_sample = np.random.random(n)*10\n",
    "    noise = np.random.normal(0, sigma, n)\n",
    "    y_meas = f(x_sample) + noise\n",
    "    X = np.atleast_2d(x_sample).T\n",
    "    return(X,y_meas)\n",
    "\n",
    "def predict_GMM_plot(lamb=0.1,sigma=1.0,l=1.0,bars=False,instances=False):\n",
    "    dA=sklearn.metrics.pairwise_distances(X_space,X_space)\n",
    "    dB=sklearn.metrics.pairwise_distances(X,X_space)\n",
    "    dC=sklearn.metrics.pairwise_distances(X,X)\n",
    "    A=sigma*sigma*np.exp(dA*dA/(-2*l*l))\n",
    "    B=sigma*sigma*np.exp(dB*dB/(-2*l*l))\n",
    "    C=sigma*sigma*np.exp(dC*dC/(-2*l*l))\n",
    "    sI = lamb*np.eye(C.shape[0],C.shape[1])\n",
    "    Cinv = np.linalg.inv(C + sI)\n",
    "    mu_p = B.T.dot(Cinv).dot(y_meas)\n",
    "    cov_p = A - B.T.dot(Cinv).dot(B)\n",
    "    plt.plot(X,y_meas,\"bo\",label=\"$\\mathbf{y}$\")\n",
    "    if bars:\n",
    "        plt.errorbar(X_space,mu_p,yerr=1.96*np.diag(cov_p),fmt=\"r\",label=\"$\\mathbf{\\hat{y}}$\",alpha=0.2)\n",
    "    else:\n",
    "        plt.plot(X_space,mu_p,\"r\",label=\"$\\mathbf{\\hat{y}}$\",alpha=0.8)\n",
    "    plt.plot(X_space,y_real,\"k--\",alpha=0.5,label=\"f(x)\")\n",
    "    if instances:\n",
    "        #L=np.linalg.cholesky(cov_p)\n",
    "        for i in range(10):\n",
    "            inst=np.random.multivariate_normal(mu_p, cov_p)\n",
    "            #z = np.random.normal(0, 1.0, mu_p.length)\n",
    "            #inst = mu_p + L.dot(z)\n",
    "            plt.plot(X_space,inst,alpha=0.5)\n",
    "    \n",
    "    plt.legend()\n",
    "    plt.xlim(0,max(X_space))\n",
    "    ylim=plt.ylim() \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1a7f68ab93ea46f78b52449b1b0f94dd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.1, description='lamb', max=1.0, step=0.01), FloatSlider(value=1.0, d…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.predict_GMM_plot(lamb=0.1, sigma=1.0, l=1.0, bars=False, instances=False)>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "interact(predict_GMM_plot,lamb=(0.0,1.0,0.01),sigma=(0.1,2.0,0.1),l=(0.1,5.0,0.1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "453b87e64d60493ba4be1abecf7c3a69",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.1, description='lamb', max=1.0, step=0.01), FloatSlider(value=1.0, d…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.predict_GMM_plot(lamb=0.1, sigma=1.0, l=1.0, bars=False, instances=False)>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X,y_meas = gen_data(n=25,sigma=1.0)\n",
    "interact(predict_GMM_plot,lamb=(0.0,1.0,0.01),sigma=(0.1,4.0,0.1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a877ab634c894b3d9244e05588d8ae48",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.1, description='lamb', max=1.0, step=0.01), FloatSlider(value=1.0, d…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.predict_GMM_plot(lamb=0.1, sigma=1.0, l=1.0, bars=False, instances=False)>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_space = np.linspace(0.1, 30.0, 1000)\n",
    "X_space = np.atleast_2d(X_space).T\n",
    "y_real = f(X_space)\n",
    "X,y_meas = gen_data(n=25,sigma=1.0)\n",
    "interact(predict_GMM_plot,lamb=(0.0,1.0,0.01),sigma=(1,40.0,1),l=(0.1,4.0,0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Periodic Kernel\n",
    "\n",
    "$K(x,x') = \\sigma^2 exp\\left(−\\frac{2 sin^2\\left(\\frac{\\pi(x− x′)}{p}\\right)}{l^2}\\right)$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict_GMMper_plot(lamb=0.1,sigma=1.0,l=1.0,p=10.0,bars=False,instances=False):\n",
    "    dA=sklearn.metrics.pairwise_distances(X_space,X_space)\n",
    "    dB=sklearn.metrics.pairwise_distances(X,X_space)\n",
    "    dC=sklearn.metrics.pairwise_distances(X,X)\n",
    "    sA=np.sin(np.pi*dA/p)\n",
    "    sB=np.sin(np.pi*dB/p)\n",
    "    sC=np.sin(np.pi*dC/p)\n",
    "    A=sigma*sigma*np.exp(2*sA*sA/(-l*l))\n",
    "    B=sigma*sigma*np.exp(2*sB*sB/(-l*l))\n",
    "    C=sigma*sigma*np.exp(2*sC*sC/(-l*l))\n",
    "    sI = lamb*np.eye(C.shape[0],C.shape[1])\n",
    "    Cinv = np.linalg.inv(C + sI)\n",
    "    mu_p = B.T.dot(Cinv).dot(y_meas)\n",
    "    cov_p = A - B.T.dot(Cinv).dot(B)\n",
    "    plt.plot(X,y_meas,\"bo\",label=\"$\\mathbf{y}$\")\n",
    "    if bars:\n",
    "        plt.errorbar(X_space,mu_p,yerr=1.96*np.diag(cov_p),fmt=\"r\",label=\"$\\mathbf{\\hat{y}}$\",alpha=0.2)\n",
    "    else:\n",
    "        plt.plot(X_space,mu_p,\"r\",label=\"$\\mathbf{\\hat{y}}$\",alpha=0.8)\n",
    "    plt.plot(X_space,y_real,\"k--\",alpha=0.5,label=\"f(x)\")\n",
    "    if instances:\n",
    "        #L=np.linalg.cholesky(cov_p)\n",
    "        for i in range(10):\n",
    "            inst=np.random.multivariate_normal(mu_p, cov_p)\n",
    "            #z = np.random.normal(0, 1.0, mu_p.length)\n",
    "            #inst = mu_p + L.dot(z)\n",
    "            plt.plot(X_space,inst,alpha=0.5)\n",
    "    \n",
    "    plt.legend()\n",
    "    plt.xlim(0,max(X_space))\n",
    "    ylim=plt.ylim() \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9fb60daeb9fc4572b007aab12f1b8237",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.1, description='lamb', max=1.0, step=0.01), FloatSlider(value=1.0, d…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.predict_GMMper_plot(lamb=0.1, sigma=1.0, l=1.0, p=10.0, bars=False, instances=False)>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "interact(predict_GMMper_plot,lamb=(0.0,1.0,0.01),sigma=(1,40.0,1),l=(0.1,4.0,0.1),p=(0.1,30.0,0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Periodic Kernel + Exponential Quadratic\n",
    "\n",
    "$K(x,x') = \\sigma^2 exp\\left(−\\frac{2 sin^2\\left(\\frac{\\pi(x− x′)}{p}\\right)}{l^2}\\right) \\sigma^2 exp\\left(-\\frac{(x_j - x_i)^2}{2\\mathcal{l}^2} \\right) $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict_GMMcomb_plot(lamb=0.1,sigma=1.0,l=1.0,p=10.0,bars=False,instances=False):\n",
    "    dA=sklearn.metrics.pairwise_distances(X_space,X_space)\n",
    "    dB=sklearn.metrics.pairwise_distances(X,X_space)\n",
    "    dC=sklearn.metrics.pairwise_distances(X,X)\n",
    "    sA=np.sin(np.pi*dA/p)\n",
    "    sB=np.sin(np.pi*dB/p)\n",
    "    sC=np.sin(np.pi*dC/p)\n",
    "    A=sigma*sigma*np.exp(2*sA*sA/(-l*l))*np.exp(dA*dA/(-2*l*l))\n",
    "    B=sigma*sigma*np.exp(2*sB*sB/(-l*l))*np.exp(dB*dB/(-2*l*l))\n",
    "    C=sigma*sigma*np.exp(2*sC*sC/(-l*l))*np.exp(dC*dC/(-2*l*l))\n",
    "    sI = lamb*np.eye(C.shape[0],C.shape[1])\n",
    "    Cinv = np.linalg.inv(C + sI)\n",
    "    mu_p = B.T.dot(Cinv).dot(y_meas)\n",
    "    cov_p = A - B.T.dot(Cinv).dot(B)\n",
    "    plt.plot(X,y_meas,\"bo\",label=\"$\\mathbf{y}$\")\n",
    "    if bars:\n",
    "        plt.errorbar(X_space,mu_p,yerr=1.96*np.diag(cov_p),fmt=\"r\",label=\"$\\mathbf{\\hat{y}}$\",alpha=0.2)\n",
    "    else:\n",
    "        plt.plot(X_space,mu_p,\"r\",label=\"$\\mathbf{\\hat{y}}$\",alpha=0.8)\n",
    "    plt.plot(X_space,y_real,\"k--\",alpha=0.5,label=\"f(x)\")\n",
    "    if instances:\n",
    "        #L=np.linalg.cholesky(cov_p)\n",
    "        for i in range(10):\n",
    "            inst=np.random.multivariate_normal(mu_p, cov_p)\n",
    "            #z = np.random.normal(0, 1.0, mu_p.length)\n",
    "            #inst = mu_p + L.dot(z)\n",
    "            plt.plot(X_space,inst,alpha=0.5)\n",
    "    \n",
    "    plt.legend()\n",
    "    plt.xlim(0,max(X_space))\n",
    "    ylim=plt.ylim() \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact(predict_GMMcomb_plot,lamb=(0.0,1.0,0.01),sigma=(1,40.0,1),l=(0.1,4.0,0.1),p=(0.1,30.0,0.1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}