GaussianProcess_tensorflow_sine.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import time\n",
    "import errno\n",
    "import time\n",
    "import timeit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import tqdm\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/data/TimeSeriesForcasting'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cwd = os.getcwd()\n",
    "cwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from keras import backend as K\n",
    "from keras import activations, initializers\n",
    "from keras.layers import Layer\n",
    "\n",
    "from keras.layers import Input\n",
    "from keras.models import Model\n",
    "\n",
    "from keras import callbacks, optimizers\n",
    "\n",
    "import tensorflow as tf\n",
    "import tensorflow_probability as tfp\n",
    "\n",
    "\n",
    "class DenseVariational(Layer):\n",
    "    def __init__(self,\n",
    "                 units,\n",
    "                 kl_weight,\n",
    "                 activation=None,\n",
    "                 prior_sigma_1=1.5,\n",
    "                 prior_sigma_2=0.1,\n",
    "                 prior_pi=0.5, **kwargs):\n",
    "        self.units = units\n",
    "        self.kl_weight = kl_weight\n",
    "        self.activation = activations.get(activation)\n",
    "        self.prior_sigma_1 = prior_sigma_1\n",
    "        self.prior_sigma_2 = prior_sigma_2\n",
    "        self.prior_pi_1 = prior_pi\n",
    "        self.prior_pi_2 = 1.0 - prior_pi\n",
    "        self.init_sigma = np.sqrt(self.prior_pi_1 * self.prior_sigma_1 ** 2 +\n",
    "                                  self.prior_pi_2 * self.prior_sigma_2 ** 2)\n",
    "\n",
    "        super().__init__(**kwargs)\n",
    "\n",
    "    def compute_output_shape(self, input_shape):\n",
    "        return input_shape[0], self.units\n",
    "\n",
    "    def build(self, input_shape):\n",
    "        self.kernel_mu = self.add_weight(name='kernel_mu',\n",
    "                                         shape=(input_shape[1], self.units),\n",
    "                                         initializer=initializers.normal(stddev=self.init_sigma),\n",
    "                                         trainable=True)\n",
    "        self.bias_mu = self.add_weight(name='bias_mu',\n",
    "                                       shape=(self.units,),\n",
    "                                       initializer=initializers.normal(stddev=self.init_sigma),\n",
    "                                       trainable=True)\n",
    "        self.kernel_rho = self.add_weight(name='kernel_rho',\n",
    "                                          shape=(input_shape[1], self.units),\n",
    "                                          initializer=initializers.constant(0.0),\n",
    "                                          trainable=True)\n",
    "        self.bias_rho = self.add_weight(name='bias_rho',\n",
    "                                        shape=(self.units,),\n",
    "                                        initializer=initializers.constant(0.0),\n",
    "                                        trainable=True)\n",
    "        super().build(input_shape)\n",
    "\n",
    "    def call(self, inputs, **kwargs):\n",
    "        kernel_sigma = tf.math.softplus(self.kernel_rho)\n",
    "        kernel = self.kernel_mu + kernel_sigma * tf.random.normal(self.kernel_mu.shape)\n",
    "\n",
    "        bias_sigma = tf.math.softplus(self.bias_rho)\n",
    "        bias = self.bias_mu + bias_sigma * tf.random.normal(self.bias_mu.shape)\n",
    "\n",
    "        self.add_loss(self.kl_loss(kernel, self.kernel_mu, kernel_sigma) +\n",
    "                      self.kl_loss(bias, self.bias_mu, bias_sigma))\n",
    "\n",
    "        return self.activation(K.dot(inputs, kernel) + bias)\n",
    "\n",
    "    def kl_loss(self, w, mu, sigma):\n",
    "        variational_dist = tfp.distributions.Normal(mu, sigma)\n",
    "        return self.kl_weight * K.sum(variational_dist.log_prob(w) - self.log_prior_prob(w))\n",
    "\n",
    "    def log_prior_prob(self, w):\n",
    "        comp_1_dist = tfp.distributions.Normal(0.0, self.prior_sigma_1)\n",
    "        comp_2_dist = tfp.distributions.Normal(0.0, self.prior_sigma_2)\n",
    "        return K.log(self.prior_pi_1 * comp_1_dist.prob(w) +\n",
    "                     self.prior_pi_2 * comp_2_dist.prob(w))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def f(x, sigma):\n",
    "    epsilon = np.random.randn(*x.shape) * sigma\n",
    "    return 10 * np.sin(2 * np.pi * (x)) + epsilon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data_size = 30\n",
    "noise = 1.\n",
    "\n",
    "X = np.linspace(-0.5, 0.5, data_size).reshape(-1, 1)\n",
    "Y = f(X, sigma=noise)\n",
    "Y_true = f(X, sigma=0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEICAYAAAC6fYRZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA2SklEQVR4nO3deXxU1fnH8c+TnZAFAiQhBAhLAgGEAAFlEUEWBdzFIq0L1QpiWxDXorViW391L9JWBat1V6i7KKIgiAqyLyZsIYQlQFgSyEb2nN8fM6EhZiOZyZ2ZPO/XK69kZu7c+72TyZMz5957jhhjUEop5Zm8rA6glFLKebTIK6WUB9Mir5RSHkyLvFJKeTAt8kop5cG0yCullAfTIu+GRGSpiNxqdY6qzidXU+6DiBgR6d4U23IkEXlNRP5qdY7zJSIjRSTdgu3uF5ExTb1dV+djdYDmSET2Ay2ArsaYfPt9vwFuMsaMrOv5xpjxTshkgFhjzN6GruN8cjljHxpLRGKANMDXGFNqcRyP5Yj3moi8BqQbY/7osGAeSlvy1vEBZlkdor5ERBsEbkps3OZvXd9rjuU2v3gP9DRwn4i0qu5BERkqIhtEJNv+fWilx1bZW/6ISHcR+da+3EkRWWS//18i8myVdX4mIndXs63V9h+3iUieiEyu+MgtIg+KSAbwHxFpLSJLROSEiJyy/xxdQ66pIvK9iDxjXzZNRMY3cNkuIrJaRHJFZLl9396q6YUVkftF5KiIHBGR26o8NlFEtohIjogcEpG5lR6ueB1O21+HISLSTUS+EZFM++v7dk2/M/v6n7evN0dENonIxZUemysii0XkDfu+JItIYqXH+4vIZvtji4CAWrbjLSLP2jOlicjv7N1SPpVe38dF5AfgDNC1jvfUOV0d9qxv2X+Osa/7VhE5aN/mw5WWbSG2rqVTIrIDGFRL7vq+16aKyPdVnmvs7/dpwK+AB+zr+KzSYgkist2+j4tEpMbXsNkwxuhXE38B+4ExwIfAX+33/QZYZf85DDgF3IytxT/FfruN/fFVwG/sP78LPIztH3YAMNx+/2DgCOBlv90W2x97RA2ZDNC90u2RQCnwJOCPrXupDXA9EAgEA/8FPq70nMq5pgIlwB2ANzDDnkcasOxa4BnADxgO5ABv1bAflwPHgD5AS+Cdyvtm368L7K9XX/uy19gfi7Ev61Npfd2BsfbXoB22fwTzavnd3mR/nXyAe4EMIMD+2FygEJhg38+/AT/aH/MDDgCzAV9gkv01+WsN27kT2AFEA62B5ZWz21/fg0Bve5YIan9P7QfGVFr/3IrXuNLr8rL9fdAPKALi7Y8/AXyH7X3bEUjC1pVS02tUn/faVOD7mp4HvFb1tbHvw3ogyp5lJ3Cn1X/vVn9pS95afwJ+LyLtqtw/EUgxxrxpjCk1xrwL7AKurGYdJUBnIMoYU2iM+R7AGLMeyAZG25e7Eds/kWPnka8ceNQYU2SMKTDGZBpjPjDGnDHG5AKPA5fU8vwDxpiXjTFlwOtAe2zFpt7LikgnbC3DPxljiu3792kt2/wF8B9jTJKxHe+YW/lBY8wqY8xPxphyY8x2bP8ka9wHY8xeY8zX9tfgBPBcHcu/ZX+dSo0xz2IrWj0qLfK9MeYL+36+ia1gAlyErbjPM8aUGGPeBzbUsZ/PG2PSjTGnsBXaql4zxiQb2/GFcdT/PVWTx+zvg23AtkrZfwE8bozJMsYcAuafxzornPNea8DzK8w3xhwxxmQBnwEJjViXR9AibyFjTBKwBPhDlYeisLXqKjsAdKhmNQ8AAqy3f/yv3D3xOraWJfbvb55nxBPGmMKKGyISKCILROSAiORga9W2EhHvGp6fUfGDMeaM/ceg81w2CsiqdB/AoVoyR1V5/JzXUUQuFJGV9i6nbGwt4rY1rUxEwkXkPRE5bN/nt+pY/l4R2WnvLjgNhFZZPqPSz2eAAHsXSxRw2NibpNVlr2M/q3tNKt93Pu+pmlTNXvG7rPU1r6dz3muNUFPGZkuLvPUexdZNUfmP7Qi21nllnYDDVZ9sjMkwxtxhjIkCpgMvyP9OF3wLuFpE+gHxwMfnma3qEKX3YmuVXmiMCQFG2O+X81zv+TgKhIlIYKX7OtaxfOXHO1V5/B1snwQ6GmNCgZf4X/7qhmT9m/3+vvZ9voka9tfe//4gtpZta2NMK2yfpurz+hwFOohI5WWrZq+6fHSl29W9JpX3p673VD62brgKkbWm/XmW2l7z+qj62p+TR0Sq5tHhc+tJi7zFjO00skXAzEp3fwHEicgvRcRHRCYDvbC1+s8hIjfI/w5+nsL25i+zrzsd20f+N4EP6vgYfAzoWkfcYKAA24HJMGz/oJzKGHMA2AjMFRE/ERlC7V0Mi4GpItLL/o+hasZgbJ8MCkVkMPDLSo+dwNZt0LXK8nnY9rkDcH8t2w7G1rd8AvARkT8BIXXupM1a+3Nn2n/n12E7rlKTxcAsEelgPxD8YB3rr+s9tRW4UUR87QeDJ9Uzd0WWOWI7MB8N/L6O5evzXtsG9BaRBPvB07kNWIdCi7yr+DO2g4QAGGMygSuwtZwzsXXJXGGMOVnNcwcB60QkD1sLdZYxJq3S469jO9BYV1fNXOB1ETktIr+oYZl52A6KnQR+BL6sY52O8itgCLbX4q/Y/ikWVbegMWYptpzfAHvt3yu7C/iziORiOyayuNJzz2A7zvCD/XW4CHgMGICtRf45toPlNVkGLAX2YOuyKKT2rqXKuYuB67AdcDwFTK5jWy8DXwHbgS3Yingp9n/w1ay/rvfUI0A3+7Yfw/aJp74ew7a/afZMjX6vGWP2YPu7WA6kAN9XWeQVoJd9HR+fR9ZmR87tAlSeRkRGYOu2iTHGlFudxxHspxfuMsY4/ZOEuxDbKacvGWOqdsmoZk5b8h5MRHyxXXD1b3cu8CIySGznq3uJyOXA1Zz/8QWPYj83fYK966UDtm6pj6zOpVyPFnkPJSLxwGlspyLOszRM40ViO+87D9vpeTOMMVssTWQ9wdZNcgpbd81ObN1PSp1Du2uUUsqDaUteKaU8mEsNBNS2bVsTExNjdQyllHIrmzZtOmmMqXrlPOBiRT4mJoaNGzdaHUMppdyKiNR4lbF21yillAfTIq+UUh5Mi7xSSnkwl+qTr05JSQnp6ekUFjpigDrlCAEBAURHR+Pr62t1FKVUHVy+yKenpxMcHExMTAznDtCnrGCMITMzk/T0dLp06WJ1HKVUHVy+u6awsJA2bdpogXcRIkKbNm30k1UzN3nBWiYvWGt1DFUPLl/kAS3wLkZ/H0q5D5fvrlFKuY6K1vu6tKxzbi+aPsSyTKp2btGSt1JmZiYJCQkkJCQQGRlJhw4dzt4uLi6u9bkbN25k5syZtS4DMHToUEfFPcfIkSPrvLhs3rx5nDlzptZllFLuS1vydWjTpg1bt24FYO7cuQQFBXHfffedfby0tBQfn+pfxsTERBITE+vcxpo1axyStSHmzZvHTTfdRGBgYN0Lq2avosWuLXj34ZEteWcfFJo6dSr33HMPo0aN4sEHH2T9+vUMHTqU/v37M3ToUHbv3g3AqlWruOKKKwDbP4jbbruNkSNH0rVrV+bP/9+E9kFBQWeXHzlyJJMmTaJnz5786le/omKU0C+++IKePXsyfPhwZs6ceXa9lRUUFHDjjTfSt29fJk+eTEHB/2b7mzFjBomJifTu3ZtHH7XNtTF//nyOHDnCqFGjGDVqVI3LKaXcl7bkG2jPnj0sX74cb29vcnJyWL16NT4+PixfvpyHHnqIDz744GfP2bVrFytXriQ3N5cePXowY8aMn51rvmXLFpKTk4mKimLYsGH88MMPJCYmMn36dFavXk2XLl2YMmVKtZlefPFFAgMD2b59O9u3b2fAgAFnH3v88ccJCwujrKyM0aNHs337dmbOnMlzzz3HypUradu2bY3L9e3b14GvnPIE2oJ3Hx5V5JvyoNANN9yAt7c3ANnZ2dx6662kpKQgIpSUlFT7nIkTJ+Lv74+/vz/h4eEcO3aM6Ojoc5YZPHjw2fsSEhLYv38/QUFBdO3a9ex56VOmTGHhwoU/W//q1avPHgPo27fvOcV58eLFLFy4kNLSUo4ePcqOHTuqLd71XU4p5R48srumKbRseXbebR555BFGjRpFUlISn332WY3nkPv7+5/92dvbm9LS0notcz4Tu1R3emNaWhrPPPMMK1asYPv27UycOLHajPVdTinlPhxS5EXkVRE5LiJJle4LE5GvRSTF/r21I7ZVm0XTh7Bo+hAu7BLGhV3Czt52tuzsbDp06ADAa6+95vD19+zZk3379rF//34AFi1aVO1yI0aM4O233wYgKSmJ7du3A5CTk0PLli0JDQ3l2LFjLF269OxzgoODyc3NrXM5pZR7clRL/jXg8ir3/QFYYYyJBVbYb3ukBx54gDlz5jBs2DDKysocvv4WLVrwwgsvcPnllzN8+HAiIiIIDQ392XIzZswgLy+Pvn378tRTTzF48GAA+vXrR//+/enduze33XYbw4YNO/ucadOmMX78eEaNGlXrckop9+SwOV5FJAZYYozpY7+9GxhpjDkqIu2BVcaYHrWtIzEx0VQ9r3vnzp3Ex8c7JKM7y8vLIygoCGMMv/3tb4mNjWX27NmW5dHfi6qLMYacglLyi0tpHxqgV0o7kYhsMsZUe762Mw+8RhhjjgLYC314DeGmAdMAOnXq5MQ47u3ll1/m9ddfp7i4mP79+zN9+nSrI6lmrLSsnIycQo7lFJKRXURGTiHHcwrJyCkkI9t+f04hhSXlAHRuE8i4XhFc1juS/p1a4+2lBb+pOLMlf9oY06rS46eMMbX2y2tL3n3o76V5KikrZ/HGQ/xjxV4ycs49KO/n7UVEqD+RIQFEhAQQGRJAZGgAPl7Cqj0nWLM3k+KyctoG+TEm3lbwh3RrQ4Cvt0V74zmsaskfE5H2lbprjjtxW0opJyorN3yy9TDzlqdwMOsMAzq14vejuxMV2sJW0EMDaB3oW2OXzNRhXcgtLGHV7hMsS85gyfajvLfhEC39vBnZM5xxvSJ4Y+1+fLy89Bx8B3Nmkf8UuBV4wv79EyduSynlBMYYvkzK4Lmv95ByPI9e7UN4dWoio3qEn3cfe3CAL1f2i+LKflEUlZaxJjWTr5KP8fWOY3y+/SgCtA3yJ6+olCB/j7qEx1IOeSVF5F1gJNBWRNKBR7EV98UicjtwELjBEdtSSjmfMYZVe07w7Fe7STqcQ7d2LfnXLwcwvk8kXg7oT/f38WZUj3BG9Qgn9Xgu4cF+7Diay4m8Igb+5Wu6hwfx+cyLHbAnyiFF3hhT/XX2MNoR61dKNZ11+zJ55qvdbNh/iujWLXh6Ul+u7d8BH2/nXDspIgQH/G94D2Ngx5EcXli1l+kjuulB2kbSK17r0Jihhk+fPs0LL7xw9nblAcuUcjVHThdwy6vrmbzwRw5knuEv1/Thm3tHckNiR6cVePj5RYwbHh7DhAva89SXu7np3+vIyNarrhtDO77q0JihhiuK/F133dUUUZVqsJRjudzy6npyC0t5eEI8Nw/p7JCzXhoyflRooC///GV/LtnUjrmfJnP586t58vq+XNY7stF5miMt8g0wdepUwsLC2LJlCwMGDCA4OPic4t+nTx+WLFnCH/7wB1JTU0lISGDs2LFMnDiRvLw8Jk2aRFJSEgMHDuStt97Si0SUpTYdOMVtr23Az8eLxdOH0CsqxJIclf8RiAi/SOxIYufWzHpvK9Pf3MQvL+zEIxN70cJPT7k8H25V5B/7LJkdR3Icus5eUSE8emXv835e5aGG586dW+0yTzzxBElJSWc/CaxataraoYSHDx/eiD1QquG+2XWMu97eTGRIAG/efiEdwxwzeYyjRoTt2i6ID2YM5dmvdrNg9T7Wp2Ux/8b+lv0jckfaJ99AlYcaPh8VQwl7eXmdHUpYKSt8sCmdO97YRPfwIN6fMdRhBd7R/Hy8mDMhnrduv5CcghImzv+OS55aeV6jszZnbtWSb0iL21kqDzXs4+NDeXn52du1Dc9bn+GGlXKk6lrRC1en8n9f7GJotzYsuHngOWe3OIIzpgkcHtuWpbMuZvRz33Ig6wx/X57CPWPjGr1eT+dWRd5VxcTEsGTJEgA2b95MWloacO4wvkq5gvJywxNf7mLh6n1MvKA9z03uh7+P7ROpq8/bWpHv9BnbpDzzV6TwyZbDfPvAKCtjuTwt8g5w/fXX88Ybb5CQkMCgQYOIi7O1Ltq0acOwYcPo06cP48ePZ+LEiRYnVc1J1X7xG15aQ9rJfE7mFXPLkM48emVvp5+D7sx/GK0DfTmQdYZPth7m6oQOTtuOu3PYAGWOoAOUuQ/9vbi+qkW+VQtfTheUMHtMHDNHdz97VlfV5S7sEga4fov+9dsGc8ur69l84BSvTB3EJXHtLE5mndoGKNMDr0p5qIqLjAZ2akWQvw85hSU8fm0fZo2J9YjTdgN8vfn3rYnERQRz55ub2HzwlNWRXJJ21yjlwQqKy9iZkUtBcRkv3jSAy/u0/9kyzjhI6kyV84UE+PL6bYOZ9NIabnttA/+dPoTYiGAL07ket2jJu1KXktLfh7swxjDnw+0UlJTx71sTqy3wnqBdsD9v3nYhvt5e3PzKeg6fLrA6kktx+SIfEBBAZmamFhYXYYwhMzOTgIAAq6OoOvznh/18vPUI94yJY3R8RJ3LN9XE987QqU0gb9w2mPziUm5+ZR2ZeUVWR3IZLn/gtaSkhPT09FrPPVdNKyAggOjoaHx9HXtutXKctamZ3PTKOi7tGc6CmwY6ZHhgd7A+LYubX1lHj8hg3rnjomYzLn1tB15dvsgrpc7PkdMFXPmP7wkN9OWT3w5z+IVOrm7FzmNMe3MTF3UN49Wpg85eB+DJ9OwapZqJwpIy7nxrE0Wl5Sy8ObHZFXiA0fERPHV9X37Ym8nsRVspK3edhqwVtMgr5SGMMTzycRLb07N57hf96B4eZHUky1w/MJqHJ8TzxU8ZDHtihdVxLKVFXikP8da6g/x3UzozL+3OOB17nTtGdCUixJ+MnCJW7jpudRzLNI+jEkp5uI37s3js02RG9WjH3WN00K6Kc/6P5djOsrnjjY30jQ7lw7uGWRnLEtqSV8rNHcspZMbbm4lu3YJ5N/ZvNmfSnI8yY9h3Ir9ZnoqtRV4pN1ZcWs6MtzaRX1TKgpsTCW3R/A60VqfqvLGPTOzF6YIS3lp30OpoTU6LvFJu7LHPktl88DRPT+pHj0i9nL8mU4fGMCKuHY9/voO9x5vX8N9a5JVyU4s2HOTtdQeZfklXJvb1zCELGquiRe/lJTwzqS+Bfj7Mem8rxaXldT/ZQ2iRV8oNJR3O5pGPk7k4ti0PXNbT6jhuITwkgCeuu4DkIzk8+/Vuq+M0GS3ySrmZkrJyHnh/O6GBvsy/sb/TJ/7wJON6RzJlcCcWrt7HmtSTVsdpElrklXIzL3+3jx1Hc/jL1b1p3dLP6jhu55Er4unSpiX3LNrG6TPFVsdxOi3ySrmR1BN5zFuewuW9Iz126GBnC/Tz4fkb+3Myr4iHP0ry+NMqtcgr5SbKyw1zPviJAB8v/nx1b6vjuLULokO5Z1wcn/90lA82H7Y6jlNpkVfKTby9/iDr92fxx4m9CA/R8fwba/qIbgzuEsajnyRxIDPf6jhOo0VeKTdw5HQBTy7dxbDubbghMdrqOB7B20v4++QEvLyEuxdtpbTMM0+r1CKvlIszxvDHj5MoKzf87dq+HjEJt6vo0KoFj197AVsOnuaSp1daHccptMgr5eI+3XaEb3Yd595xcXRqE2h1HI9zVb8o2gb5cfh0IUmHs62O43A6CqVSLiwrv5jHPttBv46t+PWwLlbH8TgVo1WezLOdSjl54Vp6tw9h8Z1DrYzlUNqSV8qF/fmzZHILS3jq+r560VMTyC8q40SeZ507ry15pVxIRcty0fQhfLPrGB9vPcKs0bE6+JiTLJo+BLC97sYYDJB6Ip/TZ4ppFegZF5ppS14pF5RbWMIfP0oiNjyIu0Z1szpOsyAi/PnqPmQXlPDMV54zto225JVyARUt+HVpWQCMfvZbjucW8eFdQ/H38a5x+YqWqGqcyq/jzRd15vW1+5mc2IkLokMtTOUY2pJXygUdzy0iMsSfAZ1aWx2l2Zk9No42Lf155JMkysvdf8gDp7fkRWQ/kAuUAaXGmERnb1Mpd1PRkrzhxTX8dCSbtkH+LLt7xM+Wq9ri1xa944W28OWhCT25Z/E2/rvpEJMHdbI6UqM0VUt+lDEmQQu8UrU7fLqAwpJy/u/aC2jpr72pVrm2fwcGxbTmiaW73H6kSu2uUcpFHMjM50ReEdf178CIuHbVLlN17tKK28qxKg7C5hSW8vQy9z4I2xRF3gBficgmEZnWBNtTyi397Ytd+Hp78eB4nenJFcS3D+GWIZ15Z/1Bfkp33ythm6LIDzPGDADGA78VkXM6GkVkmohsFJGNJ06caII4SrmeH/dl8mVyBjMu6UZEPUaY1BZ806g4CPtHNz4I6/Qib4w5Yv9+HPgIGFzl8YXGmERjTGK7dtV/RFXKk5WVG/6yZAdRoQHcMaKr1XFUJSEBvjw8sSfbDp1m8cZDVsdpEKcWeRFpKSLBFT8D44AkZ25TKXfzweZ0ko/k8OD4ngT4/vyceGWtaxI6MDgmjCe/dM+DsM5uyUcA34vINmA98Lkx5ksnb1Mpt5FfZDuwl9CxFVf1i7I6jqqGiPDna3q77UFYp56jZYzZB/Rz5jaUcmcvfZvKidwiXrppoI4T78J6RoZw65AY/rMmjcmDOtI3upXVkepNT6FUyiKHTxewcPU+ruoXxcDOemWrq7t7bCxtg/x55GP3OgirRV4pizz15S4APWXSTYQE+PLwhHi2pWezyI0OwmqRV8oCmw+e4pOtR7jj4q50aNXC6jiqnq5OiGJwTBjPfrWH/KJSq+PUixZ5pZqYMbZTJtsF+zNjpA4j7E5EhDkTenIyr4iXv9tndZx60SKvVBP7dNsRthw8zf3jeuj4NG6of6fWTLggkoWr93Eit8jqOHXSIq9UEyosKePJpbvoHRXC9QOjrY6jGuj+y3pSXFrO/BUpVkepkzYjlGpC//5uH0eyC3lucoLO2eqmKoZ3njK4E++sP8ivh8XQtV2Qxalqpi15pZrI8ZxCXliVymW9I7ioaxur46hGmjk6lgAfL5e/QEpb8ko1kaeX7aakrJw54+OtjqIaoOqELb97ZzOtAv1YmpTB5oOnXHYWL23JK9UEkg5n8/7mdH49rAsxbVtaHUc5SPvQANoG+fPEF7swxjUvkNKWvFJOVnHKZOtAP353aXer46gGqhjaueqUi2/+eIBHPk5ixc7jjOkVYVm+mmhLXiknW77zOOvSspg9No6QAF+r4ygHu3FQR7q2bcmTX+6itKzc6jg/o0VeKScqKzfMfHcLAb5eTBnU0eo4ygGqTtji6+3F/Zf1IOV4Hh9sTrcwWfW0yCvlRJ9sPUxBSRnRrVrg461/bp7q8j6R9O/Uiue+3kNBcZnVcc6h7zqlnOSGl9Yw58OfANh7Ip/JC9ae7c9VnkVEmDM+nmM5Rbz6Q5rVcc6hRV4pJzmeW0RRqev10SrnGNwljDHx4by0KpWsfNeZQUqLvFJOcKa4lPyiMgZ3CWNwTGsu7BKmk283Aw9e3pP84lL++c1eq6OcpUVeKSd4bc1+TuYV8cBlPXTGp2YkNiKYGwZ25M0f93Mo64zVcQAt8ko5XPaZEl5alcqlPcNJjNEWfHMze2wc3l7CM1+5xnAHWuSVaqCaDqQuWJ1KTmEp943rYUEqZbXI0ABuG9aFT7YeIelwttVxtMgr5UjHcwv5zw/7uapfFL2iQqyOoyxy58hutA705Ymlu6yOosMaKHW+qg5UVfky939+s5eSsnLuGRtnWT5lvZAAX353aSx/WbKD71NOMjy2rWVZtCWvlIMcyjrDu+sP8otBHXUQMsWXSUfx8/bima92Wzp4mbbklTpPNQ1Udc+irXiJMPPSWMuyKdfhJUKH1gFsPXTa0sHLtMgr5QC7M3L5aOthpl3clcjQAKvjKAtV7c7z9/Fi5ntbSJp7GV4WzAamRV6pBqp8WuQzX+0myM+HOy/pZmEi5YqiW7cg9UQ+XyQd5Yq+UU2+fe2TV6qRthw8xdc7jjFtRFdat/SzOo6yWMV1ERd2CePCLmF8NfsSYsODeO7rPZYMRaxFXqlGenrZbtoG+XHb8C5WR1EuyNtLuHdcHPtO5PPx1iPVLuPMweu0u0apRvg+5SRrUjN59MpetPTXPyf1P5W78y7rHUmfDiE8v2IPV/WLws+n6drX2pJXqoGMMTy9bBcdWrXglxd2sjqOcmEiwr3jenAoq4DFGw+dvb+iBb8uLYt1aVlOadFrkVeqgZYlZ7AtPZtZY2Lx9/G2Oo5ycSPj2pHYuTX/+CaFwpKmm1hEP18q1QBl5YZnv9pDt3Ytua5/B6vjKDdQ0Zqf8vKPvPXjAX5zcdcar7lwJG3JK9UAn207QsrxPO4d10On9VP1NqRbG4Z3b8uLq1LJLyptkm3qu1Op81RSVs685XuIbx/C5b0jrY6j3My94+LIzC/mtTX7z97nzOGotcgrdZ4+3JzO/swz3Ds2zpIrGJV769+pNWPiw1nwbSrZBSVO354WeaXOQ1FpGfNX7KVfx1aMjg+3Oo5yU7PHxpFTWMq/v9vn9G1pkVfqPCzecIjDpwu4b1ycTuunGqx3VCgT+7bn1e/TyMwrcuq2tMgrVU+FJWX845u9DO4SxvDu1o0PrjzD7DFxFJSU8dK3qU7djhZ5perprR8PcDy3iHvHaiteNV738CCu7R/NG2sPcCyn0GnbcXqRF5HLRWS3iOwVkT84e3tKOUN+USkvrErl4ti2XNi1jdVxlIeYNTqWsnLDP75Jcdo2nFrkRcQb+BcwHugFTBGRXs7YlrP7tVTz9tqa/WTlF+u0fsqhOrUJZPKgjizacIhDWWecsg1nt+QHA3uNMfuMMcXAe8DVjt7Iun2ZDHniG9bsPenoVStFdkEJC75NZUx8OP07tbY6jvIwv7u0OyLC8yuc05p3dpHvAByqdDvdft9ZIjJNRDaKyMYTJ040aCP9OrYiLNDP8rkUlWd65fs0cgpLma2teOUE7UNb8ND4nkzs294p63d2ka/u6NQ5VdgYs9AYk2iMSWzXrl2DNhLg683vR3dn88HTrNrdsH8USlXnVH4xr36fxoQLIukdFWp1HOWhpg7rwqgezrnuwtlFPh3oWOl2NFD9qPmN9IvEjnQKC+SZr3ZTXq6teeUYC1bvI7+4lNljtBWv3JOzi/wGIFZEuoiIH3Aj8KkzNuTr7cWs0bEkH8lhWXKGMzahmpnjuYW8tiaNaxI6EBsRbHUcpRrEqUXeGFMK/A5YBuwEFhtjkp21vWv6d6Bbu5Y8+/UeyrQ1rxrpxVWplJQZZo2OtTqKUg3m9PPkjTFfGGPijDHdjDGPO3Nb3l7CPWN7sPd4Hp9uO+zMTSkPdzS7gLd/PMikAdHEtG1pdRylGszjJg0Z3yeSXu1D+PvXKVzRNwpfHetbnYeKyRu6hwdhMPx+dHeLEynVOB5XAb3sM6MfzDrD+5vSrY6j3FBhSRmLNhxiyuBORLcOtDqOUo3icS15gEt7htO/Uyvmr0jh2v4dCPDV+TdV7Spa8OvSsgAQgZ/Ss62MpJRDeFxLHmxzKd43rgdHswt5d/1Bq+MoNxQRHICfj0f+eahmxiNb8gDDurdlSNc2/GvlXiYP6kign8fuqnKAiqnXBvzla06fKWbJzOG0DfK3OJVSjecxTZXJC9ae/chd4b7L4jiZV8zraw5YlEq5k6TD2WTlF9M+NEALvPIYHlPkqzOwcxijerTjpW9TySl0/lyKyr09tWw3rQJ9WXr3CKujKOUwbl/kK1rw69KyWJeW9bMW/b3jepBdUMIr36VZmFK5uh/3ZbJ6zwnuGtmNkABfq+Mo5TBuX+Tr0qdDKOP7RPLK92mcyi+2Oo5yQcYYnl62m8iQAG4ZEmN1HKUcyu2PRlYcMKtovVfcruyesXF8mZzBS6tTmTM+vknzKdf3za7jbDpwiv+79gI93VZ5HI9vyQPERgRzTUIHXl+zn+NOnEtRuZ/yclsrPqZNIDckRlsdRymHc/uWfIXqWvCV3T0mlk+3HeHKf35PTJuWdS6vmofPth9hV0Yu86f01yEwlEdqNu/qzm1a8ovEaI7nFFFUWmZ1HOUCSsrKee7rPcS3D+GKC5wzK49SVvOYlnxdJi9YS1FpGQbYeii71j581Tws2nCIA5lneHVqIl5e1U1ippT7azYteQB/n/8dVDtTrK355qyguIz5K1JI7NzaadOuKeUKmk2RXzR9CIumD2FAp1Z4ixAREqCt+Gbs9bX7OZ5bxAOX90REW/HKczWbIl/B19uL9q0CWL7zGBv2Z1kdR1kgu6CEF1elMrJHOwZ3CbM6jlJO1eyK/KLpQ/h69iVEhPjzf1/sxBidJrC5eXn1PrILSrhvXA+royjldM2uyAO08PNm9pg4thw8rZN+NzMncot49Yc0rujbnj4dQq2Oo5TTNcsiDzBpYDTdw4N46svdlJSVWx1HNZF/rdxLUWk594yNszqKUk2i2RZ5H28vHry8J/tO5rNowyGr46gmcCjrDG+vO8ANA6Pp2i7I6jhKNYlmW+QBxsSHMyimNfOWp5BfVPqzx6sbo165r3nLUxARZo2JtTqKUk2mWRd5EeEP4+M5mVfEv3UoYo+WciyXj7akc8tFnWkf2sLqOEo1mWZzxWtNBnZuzeW9I1m4OpVfXdSJtkH+P5vUWa+OdX9PL9tNoJ8Pd43qbnUUpZpUs27JV3jg8h4UlpYzf0WK1VHUeahvd9ra1Ey+2nGMVi18CWvp1wTJlHIdzb4lD9C1XRBTBnfknXUH+fWwLvUao165h7Jyw18/34GftxftQwOsjqNUk9OWvN2s0XH4+XjxzLLdVkdRdahrysfKxj73LclHciguK2fDgVN6MF01O9qSt2sX7M8dF3fl+RUp3HHoNAkdW2kL3s3lF5Vy6NQZWvp7k1+kA9Kp5kmLfCV3jOjK2+sO8LcvdvLetIt04CoXVd/utJe+TaWkzLBo+oU8uXRXrcsq5am0u6aSIH8fZo2OZV1aFt/sOm51HNUIh08XsHD1Pq7qF8WATq2tjqOUZcSVBuhKTEw0GzdutDRDSVk54/6+Gl9vYemsEXjrZBJuaea7W1iWnME3942kQys9L155NhHZZIxJrO4xbclX4evtxf2X9WDPsTw+2JRudRzVAJsPnuLTbUe44+KuWuBVs6dFvhrj+0SS0LEVz329hwKdQcqtGGP4y5IdtAv2Z8bIblbHUcpyWuSrISLMGd+TjJxCFqxOtTqOOg+fbjvCloOnuX9cD1r663kFSmmRr8GFXdtwZb8oXliVysHMM1bHUfVQWFLGk0t30TsqhOsHRlsdRymXoEW+Fg9PiMfXS3jss2Sro6h6eHn1Po5kF/LIFb30gLlSdlrkaxEZGsCsMbGs2HWc5TuOWR1H1eJYTiEvfpvKZb0juKhrG6vjKOUytMjX4dfDuhAbHsRjS5IpLNGDsK7qmWW2Gb7mjI+3OopSLsVpRV5E5orIYRHZav+a4KxtOZOvtxePXd2bQ1kFvLhKD8K6oqTD2by/OZ2pQ2OIadvS6jhKuRRnt+T/boxJsH994eRtOc3Qbm25ql8UL36byoHMfKvjqEoqTplsHejH7y7VGZ+Uqkq7a+rp4YkVB2F3WB1FVbIs+Rjr0rKYPTaO0Ba+VsdRyuU4u8j/TkS2i8irIlLtACIiMk1ENorIxhMnTjg5TsNFhARw95g4vtGDsC6jqLSMvy3dSVxEEFMGdbQ6jlIuqVFFXkSWi0hSNV9XAy8C3YAE4CjwbHXrMMYsNMYkGmMS27Vr15g4Tjd1WAyx4UHM/UwPwrqCl1bt40DmGR6e2Asfb/1QqlR1GvWXYYwZY4zpU83XJ8aYY8aYMmNMOfAyMNgxka3j6+3Fn6/uQ/qpAl7Qg7CW2nMsl3+uTOHKflFcEufajQOlrOTMs2vaV7p5LZDkrG01pSHd2nB1QhQvfZvK/pN6ENYKZeWGBz/YTpC/D49e2cvqOEq5NGd+xn1KRH4Ske3AKGC2E7fVpB6aEI+ftxePfZaMKw3V3Fy8vmY/Ww6e5tEre9M2yN/qOEq5NKcVeWPMzcaYC4wxfY0xVxljjjprW03NdhA2lpW7T7B8p04u0pQOZZ3h6WW7GdmjHVcnRFkdRymXp0erGujWoTHERQQx99NkHY64iRhjeOijn/ASePzaC3R6RqXqQYt8A1UchD18uoAXV+21Ok6z8MHmw3yXcpIHx/fUyUCUqict8o1wUdc2XJMQxUvf7tODsE52IreIvyzZQWLn1tx0YWer4yjlNrTIN9JDE+Lx8/HikU+SKC/Xg7DOUtEt9sT1ffHSYYSVqjct8o0UHhLAg+N78l3KSV5fu9/qOB5pWXIGn/90lFljYukeHmR1HKXcihZ5B7jpwk5c2jOcvy3dxe6MXKvjeJTsghIe+TiJ+PYhTBvR1eo4SrkdLfIOICI8NakvIQE+zHpviw554EBPLN3Jybwinrq+L746dIFS503/ahykbZA/T0/qx66MXJ5ettvqOG5t8oK1TF6wljWpJ3l3/SHuuLgrF0SHWh1LKbekRd6BRvUM59YhnXnl+zS+S3HdETXdQVm5Yc6HPxHTJpC7x8RZHUcpt+VjdQBPM2dCPGtSM7l38Ta+vHsEYS39rI7kNiYvWAvAurSss/fFRwbTws/bqkhKuT1tyTtYgK83825M4PSZEuZ8uF3HtmmE8GB/QnQiEKUaRYu8A1T0IVfoHRXK/Zf1YFnyMRZtOGRhMveyaPoQXpk6iABfL3y9heX3XsKi6UOsjqWUW9Mi7yS3D+/CsO5teOyzHew7kWd1HLdgjOG+xdsoLCmnW7sgQgK0Fa9UY2mRb4SKFvy6tCzWpWWd06L38hKevSEBf18vZi/aSklZucVpXd8Lq1L5MjmDP06M58u7R1gdRymPoEXeiSJDA/jbtRewLT2b55enWB3Hpa3afZxnvtrN1QlR3D68i9VxlPIYenZNI1T0F1e03qvrPx5/QXt+kRjNv1btZURcOwZ3CWvSjO7gQGY+M9/dQs/IEJ64rq8OIayUA2lLvgk8emVvOocFMnvRVrILSqyO41LOFJcy/c1NeHkJC28eqKdLKuVgWuQdYNH0IbWeBdLS34e/T04gI6eQP33iEVPdOoQxhgfe386eY7nMv7E/HcMCrY6klMfRIt9E+ndqzd2jY/lk6xH+80Oa1XFcwr+/S2PJ9qPcf1lPRsS1szqOUh5Ji3wTumtUd1oH+vLYZzv4ZOvhWpeteu69p/k+5SR/W7qTCRdEcuclOrqkUs6iRb4JeXsJ3dsFERzgw72Lt7Fqd/OZBLzyP61DWWf4/bub6R4exNOT+umBVqWcSM+uaSIVBW7DgVMABPp5c9trG/jvnUMZ2Ln1z5arGL+ltjN33FFhSRl3vrWJ0nLDgpsTaemvb0GlnEn/wizSMzKY5CM59kI/hLiIYKsjOUXVf1rDn/yGk3nFvDo1kS5tW1oZTalmQYt8E6nunPpDWWe4/sU13PLKet6fMYTo1oH1OvfenZ3MK6ZDqxZc2jPC6ihKNQvaJ2+hjmGBvHH7YM4Ul3LLK+vJzCuyOpLDVZxeGhsehABj4iP47oFRVsdSqtkQVxoKNzEx0WzcuNHqGE1u4/4sbnplHXERwbxzx0UEeVg/9TvrDvLQRz8R5O/NmjmjdeAxpRxMRDYZYxKre0xb8i4gMSaMF341gOQjOUx7YyNFpZ4xR6wxhvkrUnjoo58Y1aMdGx4eqwVeqSamRd5FXNozgqcn9WVNaiZ3v7eVsnLX+YTVEOXlhkc/Tea5r/dw3YAOLLwlUYcsUMoCWuRdyHUDovnjxHiWJmXwyCdJbjurVFFpGb9/bwtvrD3AtBFdeWZSP3y99a2mlBU8q/PXA/zm4q5k5RfzwqpUWvh6M2d8T3zcqEDmFZUy/c2N/LA3k4cm9GTaiG5WR1KqWdMi74Luv6wHZ4rLeOX7NLYcPMXzbjJ418m8In79nw3sOJrDszf04/qB0VZHUqrZc58mYjMiIsy9qjfzp/Qn5XgeE57/rs6xburDmePhHMo6ww0vrSXleC4v3zJQC7xSLkKLvAu7ql8UX8y8mB6Rwcx6byv3LNpKbqHrjUe/82gO1724hqz8Yt7+zUV6oZNSLkTPk3cDpWXl/HPlXuavSCG6dSDzp/QnoWOrs4/XdXVs1aEFLrTPTuWIq2m/SznBXW9vpqWfD2/cPthjh2dQypXpefJuzsfbi7vHxLF4+hDKyg2TXlzDv1butfQ0y6TD2Uz9z3pufmU94cH+fHDXUC3wSrkgbcm7meyCEh7+6CeWbD9KcIAP3dq1ZOuhbKDuFvr5jIdT07Ipx3J57us9LE3KoFWgL3de0o1bhnQm0E+P4Stlldpa8vqX6WZCW/jyjyn9GdkjnAfe38ZPh3OaZLsHM88wb/kePt56mBa+3swaHcvtF3fRK1iVcnHakndjaSfzmfXeFranZ+Pv48UtQzozrnckAzq1xturYRNxVO2/T+jYiiOnC8jKL8bbS7h1aAx3XtKNsJZ+DtsPpVTjOK0lLyI3AHOBeGCwMWZjpcfmALcDZcBMY8yyxmxL/VyXti15/86hXDbvW7LyS3htzX5e/i6NtkF+jImPYFzvCIZ2a0uAb8OHE9iWfhoM3DykM78d1Z2IkIBzHvfUIZGV8hSN7a5JAq4DFlS+U0R6ATcCvYEoYLmIxBljPGPkLRfi5+PFyvtsQ/fmFpawcvcJvkrOYMn2o7y34RAt/bwZ2SOccb0jGNUznJAAX0rLyjmZV0xGTiEZ2YUcz7V9z8gpxMdbyMguRASMgUkDopk5OtYtLsZSSv1co4q8MWYnUN0cnVcD7xljioA0EdkLDAY8d2ZqFxAc4MtV/aK4ql8URaVlrEnN5KvkDL7ecYzPfzqKr7cQ1tKPE7lFVD0xx8dLCA/2JyI0gLiIYPKLyggP9ufpG/pVuy1Pn6ZQKU/hrAOvHYAfK91Ot9/3MyIyDZgG0KlTJyfFaX78fbwZ1SOcUT3C+es1hq2HTvFV8jEy84uJDAkgIjSAyJAA+8/+tG3pj1cD+/GVUq6rziIvIsuByGoeetgY80lNT6vmvmqP8BpjFgILwXbgta486vx5ewkDO4cxsHOYw9bp6dMUKuUp6izyxpgxDVhvOtCx0u1o4EgD1qOUUqoRnNVd8ynwjog8h+3Aayyw3knbUhbSFrxSrq1RwxqIyLUikg4MAT4XkWUAxphkYDGwA/gS+K2eWaOUUk2vsWfXfAR8VMNjjwOPN2b9SimlGkcHKFNKKQ+mRV4ppTyYFnmllPJgWuSVUsqDudQolCJyAjhgdY4GaAuctDpEE9N9bh6a2z676/52Nsa0q+4Blyry7kpENtY0zKen0n1uHprbPnvi/mp3jVJKeTAt8kop5cG0yDvGQqsDWED3uXlobvvscfurffJKKeXBtCWvlFIeTIu8Ukp5MC3yDSAiYSLytYik2L+3rmVZbxHZIiJLmjKjo9Vnn0Wko4isFJGdIpIsIrOsyNoYInK5iOwWkb0i8odqHhcRmW9/fLuIDLAipyPVY59/Zd/X7SKyRkSqnxPSjdS1z5WWGyQiZSIyqSnzOZIW+Yb5A7DCGBMLrLDfrsksYGeTpHKu+uxzKXCvMSYeuAj4rX1Sd7cgIt7Av4DxQC9gSjX5x2ObHyEW27SVLzZpSAer5z6nAZcYY/oCf8HND07Wc58rlnsSWNa0CR1Li3zDXA28bv/5deCa6hYSkWhgIvDvponlVHXuszHmqDFms/3nXGz/3Kqd29dFDQb2GmP2GWOKgfew7XdlVwNvGJsfgVYi0r6pgzpQnftsjFljjDllv/kjtpne3Fl9fs8Avwc+AI43ZThH0yLfMBHGmKNgK2xAeA3LzQMeAMqbKJcz1XefARCRGKA/sM750RymA3Co0u3qJqCvzzLu5Hz353ZgqVMTOV+d+ywiHYBrgZeaMJdTOGv6P7dX2wTm9Xz+FcBxY8wmERnpwGhO09h9rrSeIGwtoLuNMTmOyNZE6jMBfb0nqXcT9d4fERmFrcgPd2oi56vPPs8DHjTGlIlUt7j70CJfg9omMBeRYyLS3hhz1P5RvbqPc8OAq0RkAhAAhIjIW8aYm5wUudEcsM+IiC+2Av+2MeZDJ0V1lvpMQO9pk9TXa39EpC+2bsfxxpjMJsrmLPXZ50TgPXuBbwtMEJFSY8zHTZLQgbS7pmE+BW61/3wr8EnVBYwxc4wx0caYGOBG4BtXLvD1UOc+i+0v4hVgpzHmuSbM5igbgFgR6SIifth+b59WWeZT4Bb7WTYXAdkV3Vhuqs59FpFOwIfAzcaYPRZkdLQ699kY08UYE2P/+30fuMsdCzxokW+oJ4CxIpICjLXfRkSiROQLS5M5T332eRhwM3CpiGy1f02wJu75M8aUAr/DdjbFTmCxMSZZRO4UkTvti30B7AP2Ai8Dd1kS1kHquc9/AtoAL9h/pxstiusQ9dxnj6HDGiillAfTlrxSSnkwLfJKKeXBtMgrpZQH0yKvlFIeTIu8Ukp5MC3ySinlwbTIK6WUB/t/bkYY6ni04EkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X, Y, marker='+', label='Training data')\n",
    "plt.plot(X, Y_true, label='Truth')\n",
    "plt.title('Noisy training data and ground truth')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "30"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train_size=len(Y)\n",
    "train_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "batch_size = train_size\n",
    "num_batches = train_size / batch_size\n",
    "\n",
    "kl_weight = 1.0 / num_batches\n",
    "prior_params = {\n",
    "    'prior_sigma_1': 1.5, \n",
    "    'prior_sigma_2': 0.1, \n",
    "    'prior_pi': 0.5 \n",
    "}\n",
    "\n",
    "x_in = Input(shape=(1,))\n",
    "x = DenseVariational(20, kl_weight, **prior_params, activation='relu')(x_in)\n",
    "x = DenseVariational(20, kl_weight, **prior_params, activation='relu')(x)\n",
    "x = DenseVariational(1, kl_weight, **prior_params)(x)\n",
    "\n",
    "model = Model(x_in, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def neg_log_likelihood(y_obs, y_pred, sigma=noise):\n",
    "    dist = tfp.distributions.Normal(loc=y_pred, scale=sigma)\n",
    "    return K.sum(-dist.log_prob(y_obs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- training took: 8.827615022659302 seconds ---\n"
     ]
    }
   ],
   "source": [
    "start_time = time.time()\n",
    "\n",
    "model.compile(loss=neg_log_likelihood, optimizer=tf.optimizers.Adam(lr=0.01), metrics=['mse'])\n",
    "model.fit(X, Y, batch_size=batch_size, epochs=1500, verbose=0);\n",
    "\n",
    "print(\"--- training took: %s seconds ---\" % (time.time() - start_time))    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1000/1000 [03:25<00:00,  4.87it/s]\n"
     ]
    }
   ],
   "source": [
    "X_test = np.linspace(-1.5, 1.5, 1000).reshape(-1, 1)\n",
    "Y_pred_list = []\n",
    "\n",
    "for i in tqdm.tqdm(range(1000)):\n",
    "    Y_pred = model.predict(X_test)\n",
    "    Y_pred_list.append(Y_pred)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Y_preds = np.concatenate(Y_pred_list, axis=1)\n",
    "\n",
    "Y_mean = np.mean(Y_preds, axis=1)\n",
    "Y_sigma = np.std(Y_preds, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEICAYAAAC6fYRZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA/0UlEQVR4nO3de3zU1Z34/9fJ5H4nJAFygQQEuQXCLYJUFKmIVVu1UnR1K2sVcHXd2u+2dbffde3aftvt+rOuthXQtvaytRS1Sr1VUUALCAQICReRW4BACCEhCQkhl5nz++PMNZmQ20zmkvfz8Rjmc5vP53xmwnvOnM/5vI/SWiOEECI8RQS6AEIIIfxHgrwQQoQxCfJCCBHGJMgLIUQYkyAvhBBhTIK8EEKEMQnyQrhRSr2slPqBffoapdTBPu5npVLq331bOiF6T4K8CElKqXKlVLNSqlEpVaWU+rVSKtGXx9Baf6K1vrIHZVmqlPpbh9eu0Fo/5cvyCNEXEuRFKLtVa50ITAdmAf/XfaVSKjIgpRIiiEiQFyFPa30KeBeYrJTSSqmHlVKHgEMASqlblFIlSqk6pdQWpdQUx2uVUtOUUruUUheUUmuAWLd11ymlKtzmc5VSryulqpVSNUqpnymlJgArgTn2XxV19m2dzT72+QeVUoeVUrVKqXVKqSy3dVoptUIpdUgpdV4p9XOllPLbGyYGFQnyIuQppXKBLwG77YtuA64CJiqlpgO/ApYDQ4FVwDqlVIxSKhp4A/gdkAasBb7axTEswFvAcSAPyAb+qLU+AKwAtmqtE7XWqV5eez3wI+BrwAj7Pv7YYbNbML9Gptq3u7F374IQ3kmQF6HsDXvN+W/AJuD/2Zf/SGtdq7VuBh4EVmmtt2mtrVrr3wAtwGz7Iwp4VmvdprV+FdjRxbGKgCzg21rrJq31Ja3137rYtqN7gF9prXdprVuAf8XU/PPctvmx1rpOa30C2AAU9nDfQlyWtFmKUHab1nq9+wJ7K8dJt0WjgPuUUv/ktiwaE7A1cEp7Zuk73sWxcoHjWuv2PpQzC9jlmNFaNyqlajC/Bsrti8+4bX8R8OlFZDF4SU1ehCP3oH0S+KHWOtXtEa+1fgWoBLI7tH+P7GKfJ4GRXVzM7S6V62nMlw0ASqkETNPRqe5ORIj+kiAvwt2LwAql1FXKSFBK3ayUSgK2Au3Ao0qpSKXUHZhmGW+2Y74UfmzfR6xSaq59XRWQY2/j9+YPwD8opQqVUjGYZqVtWutyH52jEF2SIC/Cmta6GNMu/zPgPHAYWGpf1wrcYZ8/DywBXu9iP1bgVuAK4ARQYd8e4CNgH3BGKXXOy2s/BP4deA3zRTEGuMsHpydEt5QMGiKEEOFLavJCCBHGJMgLIUQYkyAvhBBhTIK8EEKEsaC6GSo9PV3n5eUFuhhCCBFSdu7ceU5rneFtXVAF+by8PIqLiwNdDCGECClKqa7u1JbmGiGECGcS5IUQIoxJkBdCiDAWVG3y3rS1tVFRUcGlS5cCXRQRomJjY8nJySEqKirQRRFiwAV9kK+oqCApKYm8vDyUDJYjeklrTU1NDRUVFeTn5we6OEIMuKBvrrl06RJDhw6VAC/6RCnF0KFD5ZegGLSCPsgDEuBFv8jfjxjMQiLICyGE6BsJ8j1gsVgoLCxk8uTJLF68mIsXL/Z5X0uXLuXVV18F4IEHHmD//v1dbrtx40a2bNninF+5ciW//e1v+3xsEYSOfQKN1YEuhQhjEuR7IC4ujpKSEvbu3Ut0dDQrV670WG+1Wvu035deeomJEyd2ub5jkF+xYgVf//rX+3QsEYRam+D4Zji8vvtthegjCfK9dM0113D48GE2btzI/Pnz+bu/+zsKCgqwWq18+9vfZtasWUyZMoVVq1YBpnfHI488wsSJE7n55ps5e/asc1/XXXedM43De++9x/Tp05k6dSoLFiygvLyclStX8tOf/pTCwkI++eQTnnzySZ5++mkOHDhAUZFrlLry8nKmTJkCwM6dO7n22muZMWMGN954I5WVlZ3OYenSpTz00EPMnz+f0aNHs2nTJu6//34mTJjA0qVLndu9//77zJkzh+nTp7N48WIaGxsB+M///E9mzZrF5MmTWbZsGY6BZ6677jq++93vUlRUxLhx4/jkk098++aHm9pjoDWcL4cLVYEujQhTQd+F0sM3vwklJb7dZ2EhPPtsjzZtb2/n3XffZdGiRQBs376dvXv3kp+fz+rVq0lJSWHHjh20tLQwd+5cFi5cyO7duzl48CBlZWVUVVUxceJE7r//fo/9VldX8+CDD/Lxxx+Tn59PbW0taWlprFixgsTERP7lX/4FgA8//BCACRMm0NraytGjRxk9ejRr1qzha1/7Gm1tbfzTP/0Tb775JhkZGaxZs4bvfe97/OpXv+p0LufPn+ejjz5i3bp13HrrrWzevJmXXnqJWbNmUVJSQk5ODj/4wQ9Yv349CQkJ/Nd//RfPPPMMTzzxBI888ghPPPEEAH//93/PW2+9xa233up8j7Zv384777zD97//fdavl1pql+pOuKZP74YrFwWuLCJshVaQD5Dm5mYKCwsBU5P/xje+wZYtWygqKnL2vX7//fcpLS11trfX19dz6NAhPv74Y+6++24sFgtZWVlcf/31nfb/6aefMm/ePOe+0tLSui3T1772Nf70pz/x+OOPs2bNGtasWcPBgwfZu3cvN9xwA2CakUaMGOH19bfeeitKKQoKChg2bBgFBQUATJo0ifLycioqKti/fz9z55qxqltbW5kzZw4AGzZs4Cc/+QkXL16ktraWSZMmOYP8HXfcAcCMGTMoLy/v9jwGtYZTrumqvTDmeojsaixwIfomtIJ8D2vcvuZok+8oISHBOa215vnnn+fGG2/02Oadd97ptguf1rrX3fyWLFnC4sWLueOOO1BKMXbsWMrKypg0aRJbt27t9vUxMTEAREREOKcd8+3t7VgsFm644QZeeeUVj9ddunSJf/zHf6S4uJjc3FyefPJJjz7ojn1ZLBba29t7dU6DSlszNLmN+W1tg7P7IaswYEUS4Una5H3kxhtv5IUXXqCtrQ2Azz//nKamJubNm8cf//hHrFYrlZWVbNiwodNr58yZw6ZNmzh27BgAtbW1ACQlJXHhwgWvxxszZgwWi4WnnnqKJUuWAHDllVdSXV3tDPJtbW3s27evT+cze/ZsNm/ezOHDhwG4ePEin3/+uTOgp6en09jY6PzlInqp4TQAa3eeZO3Ok2ZZZUngyiPCVmjV5IPYAw88QHl5OdOnT0drTUZGBm+88Qa33347H330EQUFBYwbN45rr72202szMjJYvXo1d9xxBzabjczMTD744ANuvfVW7rzzTt58802ef/75Tq9bsmQJ3/72t51fDtHR0bz66qs8+uij1NfX097ezje/+U0mTZrU6/PJyMjg5Zdf5u6776alpQWAH/zgB4wbN44HH3yQgoIC8vLymDVrVq/3LYD6is7LGirNBdikYQNfHhG2lKNnRDCYOXOm7jhoyIEDB5gwYUKASiTCRbD9Hf38p0+Raa3kVF0zANmpcQAs/vJXYNyNl3upEJ0opXZqrWd6WyfNNUIMNJuNNNs57+uq9kJ768CWR4Q1aa4RYqA1VXP39OEAzvb4xTNyzbr2VrkAK3xKavJCDDT3rpPeyAVY4UNSkxdioLkFeWcN3mN9JVw4A0nDB7BQIlxJTV6IgVbfTU0e4HSJ34shBgefBHml1K+UUmeVUnvdlqUppT5QSh2yPw/xxbGECGmtTdB83jXfcBGOnum8XdVeaG8ZuHKJsOWrmvzLQMfEG48DH2qtxwIf2udDTk1NDYWFhRQWFjJ8+HCys7Od862tl+8FUVxczKOPPtrtMa6++mpfFdeDewK0rjz77LP9Sp0sesm9Ft9uhYd+AQ88BxtLPbeztkFV325kE8KdT4K81vpjoLbD4q8Av7FP/wa4zRfHGmhDhw6lpKSEkpISVqxYwWOPPeacj46Ovuyt+zNnzuS5557r9hju6YQHmgT5AeZ+0fVv+6HS/t/mNS9/A6d3myyVQvSDP9vkh2mtKwHsz5neNlJKLVNKFSuliqurfTN4wpJVW1myqvv8LX21dOlSvvWtbzF//ny++93vsn37dq6++mqmTZvG1VdfzcGDBwGTD/6WW24B4Mknn+T+++/nuuuuY/To0R7BPzEx0bn9ddddx5133sn48eO55557nGl833nnHcaPH88XvvAFHn30Ued+3TU3N3PXXXcxZcoUlixZQnNzs3PdQw89xMyZM5k0aRL/8R//AcBzzz3H6dOnmT9/PvPnz+9yO+FD7kH+k32QlgQPLIR9J+BEh7//xrPO9AdC9FXAe9dorVcDq8Hc8Rrg4vTY559/zvr167FYLDQ0NPDxxx8TGRnJ+vXr+bd/+zdee+21Tq/57LPP2LBhAxcuXODKK6/koYceIioqymOb3bt3s2/fPrKyspg7dy6bN29m5syZLF++3JmK+O677/ZaphdeeIH4+HhKS0spLS1l+vTpznU//OEPSUtLw2q1smDBAkpLS3n00Ud55pln2LBhA+np6V1u58hVL/rJZjU9Z8DU0EvLoTAfFs2AX62H93bCsg6tnqd3Q0r2gBdVhA9/BvkqpdQIrXWlUmoEcLbbV/STo/a+7Vitx/ya5XN8fqzFixdjsVgAk1b4vvvu49ChQyilnEnKOrr55puJiYkhJiaGzMxMqqqqyMnJ8dimqKjIuaywsJDy8nISExMZPXq0MxXx3XffzerVqzvt/+OPP3ZeA5gyZYpHcP7Tn/7E6tWraW9vp7Kykv3793sN3j3dTvRBYxXY7M17p2uhpgGm5Jva/Jzx8NedcP8NEGn+rtbuPImN0yy5YgFExQWw4CKU+bO5Zh1wn336PuBNPx5rwLmnGf73f/935s+fz969e/nLX/7ikXrXnXtK365S8Xrbpjf5hbylLD527BhPP/00H374IaWlpdx8881ey9jT7UQfuV903WOSyjHVfHHzpZlwvgm2fubxkgiscKZsgAoowpGvulC+AmwFrlRKVSilvgH8GLhBKXUIuME+71drls9hzfI5XJWfxlX5ac55f6uvryc72/ykfvnll32+//Hjx3P06FHnIBxr1qzxut28efP43//9XwD27t1LaanpsdHQ0EBCQgIpKSlUVVXx7rvvOl/jns74ctsJH6g/6ZrecwyGJMDIDDNfNBaGJsM7xc70w6fqmjlV18yv177OkpWBuzgvQptPmmu01t4biWGBL/Yf7L7zne9w33338cwzz3gd+am/4uLi+MUvfsGiRYtIT0/3GN/V3UMPPcQ//MM/MGXKFAoLC53bTZ06lWnTpjFp0iRGjx7tHO0JYNmyZdx0002MGDGCDRs2dLmd6CetPS+6lh6Dgnxw/PKyWGDRdHhlE3FfLKI5NdG5aaKtgUyrl770QvSApBoOEY2NjSQmJqK15uGHH2bs2LE89thjgS5WyAj431FzHXz6gpk+cx7+7r/hkVvgDrd7JE7Xwr1PU7aoiILv3OaZvCxjHEz+6sCXW4QESTUcBl588UUKCwuZNGkS9fX1LF++PNBFEr3RsRYPrvZ4h6w0zo7JIn/7Z2Czea47dxguNfi3jCIsBbwLpeiZxx57TGruoazjRdfkOMh3jQDlqLUnTRrDonWfsPHV7TAm25XATNt4avXv2RszbUCuM4nwITV5IQaC+0XX0nIoyIOIzv/9Do0fRUtMlKnNd3BF20Es2nv3XCG6IkFeCD+7Z+Um1m7aZWbONcCpGkqGeObrWzwjl8UzchmekUTFtLGM2nuUxVeanjeO3jZV5xuoP17q9zu6xQBrqoHDH/pt9xLkhfCzNGsNYO/gYG+Prx6T1eX2x4omQGs7fLSn07qpEUcln004qT4IO170azI6aZMXwk8ctW116ginIpr5xabDXPPeHq6MjmR/fDx1HYf+c0xPz4G3NsO7O+HLVznXr915kmwu8f/uHiWpDkJdSyNUbIdTu/z+pR16QX7Dj3y7v/n/2u0mFouFgoIC5/xdd93F4493nTl55cqVxMfH8/Wvf93r+o0bNxIdHe2zFMMPPPAA3/rWt5g4caJP9udrzz77LMuWLSM+Pv6y2/XkPN544w3GjRsXtOfqzQhqnNNZJ6qozMlEO9vjFVzxRTi83vUCpeCmGfDzt02u+dEdRog6vUuCfKhqroOT26FyjyvFhZ+FXpAPgLi4OEpKSnq8/YoVKy67fuPGjSQmJvosyL/00ks+2Y8/WK1Wnn32We69995ug3xPzuONN97glltuCYkgv2b5HLDZeOWnfyRSx7F4XDpU11E5bSzZqXGmhp6WD9nT4fjfoM0thcQXC2HVe6Y2//DNgFuN/+xnMGYBRF/+/RRB5GItHN9immW0rfvtfUja5PshLy+P7373uxQVFVFUVMThw4cBk1b46aefBkw634kTJzJlyhTuuusuysvLWblyJT/96U8pLCzkk08+obq6mq9+9avMmjWLWbNmsXnzZud+7rvvPhYuXEheXh6vv/463/nOdygoKGDRokXORGjug4O89957TJ8+nalTp7JgQecbjl9++WUeeeQR5/wtt9zCxo0bAZPy+Hvf+x5Tp05l9uzZVFVVAVBVVcXtt9/O1KlTmTp1qjP//e9//3uKioooLCxk+fLlWK1W536eeOIJrrrqKn74wx/2OJ2x+3l4K8uWLVtYt24d3/72tyksLOTIkSMemTYPHTrEjBkz+vOR+l5TNZGOHjF7j4PWVOePcK0fkgcRFsjuUO6UBJg7AT7YDW0dany2djjTYZAREZwaq2HfG7B9tclBNMABHiTI90hzc7NzNKjCwkKP3DHJycls376dRx55hG9+85udXvvjH/+Y3bt3U1paysqVK8nLy/MYfOSaa67hn//5n3nsscfYsWMHr732Gg888IDz9UeOHOHtt9/mzTff5N5772X+/PmUlZURFxfH22+/7XGs6upqHnzwQV577TX27NnD2rVre3WeTU1NzJ49mz179jBv3jxefPFFAB599FGuvfZa9uzZw65du5g0aRIHDhxgzZo1bN68mZKSEiwWizNvTlNTE5MnT2bbtm088cQTZGVlsWHDBjZs2ACYdMbFxcWUlpayadMmZ46d7spy9dVX8+Uvf5n//u//pqSkhDFjxpCSkuL8lfXrX/+apUuX9uqc/a6+wtlzhrJyiLQw/yvTXbXyFHsW0hFTXSkOHG6aaYYH3HKg835P7+58w5QIDq0XoXwzlL0KO16CswcCerFcmmt64HLNNY7c7nfffbfXm5WmTJnCPffcw2233cZtt93mdR/r169n//79zvmGhgZn0rCbbrqJqKgoCgoKsFqtLFpk8o0XFBQ4E5Y5fPrpp8ybN8+ZkjgtLa03p0l0dLRzMJIZM2bwwQcfAPDRRx/x29/+FjDXJ1JSUvjd737Hzp07mTVrFmC+CDMzM53bfPWrXd+C35N0xl2VpaMHHniAX//61zzzzDOsWbOG7du39+qc/UprqDvumi8thyuzIcY+hkCEBRLt7e2xKaZWX3vMtf2MKyAjxTTZXOu6JgSYtt3ao5B+hR9PQPRKa5Npbz+10wzfGCQkyPeTe2pfb2l+3377bT7++GPWrVvHU089xb59nbtK2Ww2tm7dSlxc55zhjtTDERERREVFOY8RERHRKVWx1tprGdxFRkZic6sBuqcSdt9/V6mQ3Y9133338aMfdb4QHhsb68y135EjnfGOHTsYMmQIS5cu9ZrOuKdl+epXv8r3v/99rr/+embMmMHQoUO7LPOAaz5vusgBtLTB56fgTrekb8lZYHH7Lzi8wDPIWyLgxunwh41QXW8CvrvTuyTIB4PWJji5zR7cB+Ziam9Ic00/OZpu1qxZw5w5nreb22w2Tp48yfz58/nJT35CXV0djY2NHul9ARYuXMjPfvYz53xvLvK6mzNnDps2beLYMRMoams7DrtrriOUlJQ4y9aTmu+CBQt44QWTXMtqtdLQ0MCCBQt49dVXOXv2rPNYx48f9/p6X6Yz7vjexcbGcuONNzozcAYV97tcD5w0A3cX5LmWpXgOGEP6OIiM9ly2aAbYNPx1l8fitTtPsvb9TeaCngiM1iY48pFJPHdiW1AGeAjFmnwPujz6mqNN3mHRokX8+McmPX5LSwtXXXUVNpuNV155xeN1VquVe++9l/r6erTWPPbYY6SmpnLrrbdy55138uabb/L888/z3HPP8fDDDzNlyhTa29uZN28eK1eu7HU5MzIyWL16NXfccQc2m43MzMxOzRxz584lPz+fgoICJk+e7HHhsiv/8z//w7Jly/jlL3+JxWLhhRdeYM6cOfzgBz9g4cKF2Gw2oqKi+PnPf86oUaM6vd6X6YzvuusuHnzwQZ577jleffVVxowZwz333MPrr7/OwoULe7Uvv3PPV1NWbtrcJ7u9P6kjPbe3REHGBNO9ziErDQpHmyabv7u2QyoEbdrmrxgUGb2DQ9slsLVBxY6grbl3JKmG+yEvL4/i4mLn+KgiMJ5++mnq6+t56qmnutwmIH9H21bDRXsf+e/8GmovwEtmeEYiLDD3m51r7nUnYffvPZd9sBt+tBaeeYC1VtOef6rODNKeMSSVtxPu4A8rrvHjiQhaLpjAfnq3aW/3ddyMToC5j/b55ZdLNRx6NXkh3Nx+++0cOXKEjz76KNBF8dTa5ArwVivsOw4Lp7nWJ43oHODBNOHEpZoLqw7zJsNzfzG1+YWzPTaP1i2MbD8GSJD3i4u1pr39TJkZiD0ESZDvh469W8TA+/Of/xzoInjn3lRz5Aw0t3q2xw/p3KwFmCadYZOh/G+uZTFRcP0UeH83ix+5BZLiPAYUWZxor1l2c9Fd9ELjWTixNeDdH30hJC68BlOTkgg9Afn76ZhaGDyDfMf2eHfDJ3de9uWrzE1Rf9jYeV3jWaiv6EMhRSf1FVC6Fnb8Eqr2h3yAhxAI8rGxsdTU1EigF32itaampobY2NiBPbB70C0rhxFDXF0gIyIhOcfrywCIG9K5582YETB3oull02513WDlcGqnz4o+aB3dCLt+BzWHA10Snwr65pqcnBwqKiqorq4OdFFEiIqNjSUn5zJB1desbXDBPvC21ibIF41zrU/J9uwf783wgs6185tmwCf74NOD8IUOuXuqD5qLgzFJ/S7+oHOxFiqKPXs1hZGgD/JRUVHOOziFCAkNp105Sk6eg7ommJznWn+5phqHjPFw6APPTIWzxsLQJHinuHOQ1zY4XQL5cgG2R7Q2TWont5uaexi3FAR9c40QIce9Pb6s3DxPuUz/eG+iYiF9rOcyi8XcAbv9oBlhqqPTu0O2B8iAsdlMFs9dv4Hd/wvnDoV1gAcJ8kL4Xsf2+NQEyDVD+WGJhKSuR4XyMMzLBVjHHbAf7O68rrUJqjuPDTuoaQ2XGqC9FSp2wvZVsO/P0FAZ6JINmKBvrhEipNhsHYL8cXOXq6N7Y3JO9+3xDmn5Jmd860XXspx0mJJn+szfNa9zt8mKYhg2qV+nEBZsVtP98cRWaGsGbfXM1z+ISE1eCF9qOuvKQHiuASpre951sqMIC2R6Cdg3zYSKcyY/fUcNp+F8eW9KHF6sbeaLbtsqOPAXaDpnfuEM0gAPEuSF8K06b+3xea5lXd0E1RVvtfJ5kyE+xtTm7dbuPOm8QYrDH4Z9OzMAl+pd061NcOwT2Ppzc8Hafd0gJ801QvhSx5ug4qLhCvtIUJZIk86gN5KGQ0K6qZE6xEXD/Cnw4R545BYT8N01noW6E73/QgkFWptzO7EVmqqh8B7TQ+ZM2YCNmRpqJMgL4Stad77oOnGk6RUDkJJrmmB6QylTmz+6yXP5TTPg7R3s+M0myosmOBOWOdMdZBSHV5DXGmqOwIktnikjtq0KXJlChAR5IXxgyaqtJNoa+OX4JrOgsRmOVcG1bj1ketMe7y5zYucgPyEXRmWQv/0zyou8ZNc8d8gMWhI3pG/HDBY2G1TbL6A2yg2RfSFBXggfybBWuWbsg3b3+aKru7hUSM31bO9XCm6aSfrKd1mcHoNjNF9nqgOtTaqDK77Yt2MGmtU+WPnJbZ4ZOUWvSZAXoh+WrNoKwLZjtaREHGXtJTNS1uLScoi0wHh7OoW+tMe7GzbZM8gDfLEQXvwrvLcTZnrpU1+5B/KugciYzuuC3d7XzBi2ot+kd40QPpKtXBdHz239nJqsdIi154zvS3u8u4zxnV+flgSzx8P7u1g8NcszYRmYG4DOlPX9mIHU1hToEoQNCfJC9MOa5XNYs3wO146KZeIQm8kOOXk4aRVnqc53q7n3tanGISoWhnoZtPtLM+F8E2w76P11FcWh2Z2yvSXQJQgbfg/ySqlypVSZUqpEKVXc/SuECD0ZVtNMs3bnSTas20WE1cbBzKGu/uv9DfIAY67vfIdr0VhTo3+3i1TDzedDM3Vu++C9ecnXBqomP19rXdjVGIRChLqnF6Y5m0vSj5m8KKdzMwGw0s/2eIe4VNPs485iMcMKfnoQarwkLQMzNmko0Vpq8j4kzTVC+EK9azi+gprz1A9PI31EKotn5HLXF+f0rz3enbc7YG+aaboaflDi/TXnj5sbpIKdzQbnDpuLrqHYxBSkBiLIa+B9pdROpdSyjiuVUsuUUsVKqWIZGESEpLZmVx9uqw32nfBte7y7jPFmZCl3uelQMArevUz7ezDX5i/WmlGZPv05lK01ffyFzwxEkJ+rtZ4O3AQ8rJSa575Sa71aaz1Taz0zIyNjAIojhI+53+V6pBIutnDF9ZNdvV18GeSjYmHomM7Lb5ppBijZd8L766r2m/wugWCzme6Qn73tSp7W3mK6eO7+vblr9fhWaGkMTPnCnN/7yWutT9ufzyql/gwUAR/7+7hCDJg6t8DaMSlZhI/a490Nm2yG+3N37WR4/i+mNj/ZSzoDW7sZOSpvrm/L0hWbDepPmAE6zh10pUu2tUNlqVlmlVwzA8GvQV4plQBEaK0v2KcXAv/pz2MKMeA6jgQ1LNU1aHdPxnPtraFjTI3ePX1uXIxJWvZRqUlaFuflBqjTu2DkbN9dH+jIEdirD5rBS9zz4DtU7ffPsUWX/F2THwb8WZluX5HAH7TW7/n5mEIMnPZWuGBPZ+AYtHuG27B9vmyqcYiwQMYEM9yfu5tmmvFfN5aZ6Y5aGs1AGsO93B3bVy2NcHa/SR7WWGWuT4ig4tcgr7U+Ckz15zGECKiGU65Bu0/VmBuT3PPHJ2f757jDJnYO8hNzYWQGvF1shgns2KceTC6YYZO8r+up9laoOQRV+0xbu/SECWqSu0aI/vA2aLd7UrKk4f45bkouxKZ4Do6hFHz5KvjZW1B8GGaN7fy6xrOm1p3u5e7Zy7FZofYYnN1ner84Rr8SQU/6yQvRH+4XXUvLISXe1KbBpPmNivPPcZWCTC8phm8tMmV4a3vXrz22ybSf94TW8Nk7sOU5072xar8E+BAjQV6IvrK2Q0Ola76sHCbnuZpC/NEe726Yl7b1qEhYOB22HIDaC95f13jWXITtiebzpqvjIB4jNdRJkBeiry6cdg05V9MAp2s92+NTc72+zGcSMyAxs/Pym2eam7Le3915ncORj0zbelfaLpk290Mf9L+cIqCkTV6IvvI2aLd7e3zHPDP+MGxSp5QFa6tbuC5/OBnvFMOSa7xfZLVZ4UKla4hAre3JzI6YhGZ1J1wXlEVIkyAvRF91bI+PjYax9hufYlNMQjF/y5xoUgJ06OFyrGgCGWs2mHJNzff+2rK1kDbafAk0VHpexBVhQ4K8EH1hs0KD+6Ddx2GS26Dd/m6qcYhNNr8Y6k44B/E+VddM1cgspsVGc+p3m8jrKshb2zrfOSvCjrTJC9EXFypdt+U3NsPRM57pBAaiqcZh2MROi9qjozg+bSy5pUdM+cSgJUFeiL6oP+Wa3nvCd4N294V9aMDFM3JZPCOX7NQ4slPjuOK+67C0W+HDPQNXFhF0JMgL0ReX6lyjPu0tB0uEueMUICbR9JEfKFFxpm29o3HZcEUWvL1D7kodxCTIC9EXl9xGYSorNwHVOWh3Tv/SBvSF22Aijho9YLpTHq6EQ6cHtjwiaEiQF6KXlqzaym/W7+BUXTNnzl3Auv8kBzPTXBsMZFONw9ArwBLVefmCqRATZWrzYlCSIC9EL8XZmoi3mQEuhp2uwWK1eY4ENZAXXR0sUZBxZefliXEm1/xHe6D5Mjc/ibAlQV6IXnr59uHOC5wTqmsA+MJX7Kl9o2IhIUAjnGV27mUDwJdmQlMLbCob2PKIoCBBXojecss8mX6sEvIyTVIwMLX4gW6PdxiS5z0hWkGeGQf2neKBLpEIAhLkheihJau2smTVVmc6g8WF2YyoODvwqQy6EmHxnplSKfjSLNh7HI6f7bxehDUJ8kL0QoztElw0TTQcPWOaQTz6xwcwyEPXTTYLp5lunlKbH3QkrYEQ3ViyaisA247VcoWqYG2TvSZ/3N5s4wjylihIHBaAErpJyek8mAjAkESYOwHe3wXfWAjR8l9/sJCavBC9kKPOuWbKyiEz1QzcDWaoP38Nkt1TSkH2dO/rbp4F9RfNGLBi0JCvcyG6sWb5HMDU6L/QdJHF03LNHaQ/Og7Tx7g2DHRTjUOHzJSOxGWLZ1wBo4fD/24w/ectUscbDORTFqKHonQLKbY6M3O61oy8VBCgpGSX48hM2VFEBNw7H06eu/zwgCKsSE1eiB76/Z3ZUJZjZkrLzbOjPT7CAslZgSiWd5kTWPvhZsCkHgZ7jT4xhcXTRsNLf4VpY1zj0YqwJTV5IXqq7rhruqwckuNhlH34vaTh3tMKBErGeLSX/97VTa28vWiOGQv2X35pMlTWX5QEZmFMavJC9FTH4f4KRrlufErJCUiRuhQdz9cWzoOaI25t8rms3XmSi8TAf98PP1xjHgBRFpPjJirS9LxxPF9qheoGsz7SYpY7pzs+R5pnxzL3bR3LRmXCohkBfGMGHwnyQvREews0Vpnp2gtwqgZuLXKtD5b2eHeZE8yYrUB1Ywtrd550Nd0Qh1pxO3dGtZsslecbobUNWtuhzWp/bod2G8y+ElQEtNvXtVs7P7e3my8E5/L2zts1XTLl+tUHkBhrvlSio8xzTKRrut0Kx6rMF2i025dHVKSZd0zXXoCKGvNa931FR7r2Fx0JJUehrgkiI8zIXVEW8xwZAZGR9mcLpCRAaoK5IB2hzDUMS0Tn+Zx00yXVfZ0lwrU+QkHNBThZbfbrPJbFrHc8WyJMZaHmAkRGw63x8MADPv8zkCAvRE/UV7iaNDq2x0Pw1eQB0sdBRKSzBt+RtkSYdvlpY7y82A+aW+H3G6C6HlrazKO1zQT/Wvt0S7tZnpUGqYmuL4vmFtOs5PiiaWsHm4YxI0xAd7zu4iU43+7ad0ub+UKYM971BdJmBavVfIE5p63m5raWNrDZwGqzP2vXtE37p1krIQZiY2D4DgnyQgRMx/b42CgYa7/QmpDuPWdMoEXGwNDRUP25M7+8e9PNgIuLhgdvHPjj+lJrGxyqNF827l8AVvuXgGNaAWOzzS8Jq83+a8fqmnY8x8fYv6gSYe6jfimyBHkhesKjPf44TMg1P7shOJtqHDInQvXngS5F+IiOMgO2hxAJ8kJ0p70VLpwx002X4Gil6W/uEIxNNQ6OwUSsbUCAavAioKQLpRDdqT8J2mam950wP8un5LnWB8udrt5YokzbvBi0JMgL0Z26E67psnLTK2KC/Sd7bLJJCBbMuspMKQYFCfJCdMdtkBBKy80F1zj7oN3J2QEpUq+k5ZsRq8SgJEFeiMtpb4WGSjPd2g6fVQTPICE9FWGBjPGBLoUIEL8HeaXUIqXUQaXUYaXU4/4+nhA+1VDhao8/WGH6Z4dKe7w7abIZtPwa5JVSFuDnwE3AROBupZT8tYnQ0bE9HmBynnmOjIb49IEuUd+kjoSYpECXQgSAv2vyRcBhrfVRrXUr8EfgK34+phC+4x7kS8tN7hXHoN3JOeZW9lCgFAyfHOhSiADw919oNuB+P3WFfZmTUmqZUqpYKVVcXV3t5+II0Qvu7fFWG+w73iF/fBD3j/dm+JRAl0AEgL+DvPKyzCP5g9Z6tdZ6ptZ6ZkaG5LYWQcS9Pb68qvOg3aEW5OPTYEheoEshBpi/g3wF4H5lKgc47edjCuEbHZtqwBXkVURwDRLSUzkzA10CMcD8HeR3AGOVUvlKqWjgLmCdn48phG+c75CULDMFhg8x80nDgmuQkJ4aeoUpuxg0/BrktdbtwCPAX4EDwJ+01vv8eUwhfMI9X43W9kFC8lzrU0MrSZWTUjD2RvNLRHlrTRXhxu8JyrTW7wDv+Ps4QviUPV/N2p0nSThXz5dqLngG+fQrA1a0fkvJhnE3QkIGnN4FVftd1x5E2JEslEJ449Yen37M3sPGEeQtkWZM11CWVWieU7Jh9Hw4ux+q9sKFqoAWS/ieBHkhvFi1bgNp1nOcqmtmwsEKWuJjWHeuhcX5mHw1EZZAF9F3YhIht8g8mmrg7D5Tu28+H+iSCR+QIC9ER+2tDLHWOmezT1ZxLm+EGbsTQiMpWV8lDIX8eZB3DTScglM7TcAXIUuCvBAd1Z/kazNMIF+34SBDahsYcufVrgE3Qq1/fF8oZc4zbogE+RAXIvdki1C0ZNVWlqzaGuhi9J57e3y5vT3ekZRMqfCuyXcUnQCzvgE5syRdcYiSmrwQHbkF+aubGiGm46DdgyzYJWbC2C/C6Ovg3EE4XeJ5o5gIahLkhc85au/bjtV6zK9ZPidgZeqx9hZX/3gw/eMn5EKU/b9K8iBoqumKJRKGTTKPpnNwejecKTX3FIigJc01Qrg7X+7qM950CY5UeuaPHwzt8T2RkA5jb4A5/2T63A/JM6mXRdCRmrzwOUeNPaRq8A7nPndN77cP2u2RlGwQtcf3RGQ0ZE83D2sbVB80tXv3lBAioCTIC+GgtWdwKi03+eIn2nvVxCRCbGogShYaLFEmZ/3wyaaP/Zky87jUEOiSDWoS5IXfuNfgQ6JW33weWi645svK7YN2x5j55GzJ99JTcUNMf/tRX4C6cqjcA+cOgc0a6JINOhLkhXBw7zHS2g4HKuC22a5l0h7fexERkDbaPFovQluzCfoVO+FiTaBLNyiET5DXWmpZQSiketrUuTXVfH7KDNodyoOEBJvoePNIGApZ06H2KJzcJu33fhY+Qf7kdqgqg5Rc87M6JQdiUzoH/kPrTdth6kizTSjmBBe+p3UXg3bbUwpHREKi5GH3GaVg6BjzuFBlgv3ZA5IN0w/CJ8ijobHaPE7tMotiEk0gT84xz4mZJvlS60U4vsUkmUoabg/4uWabyBjTV7rsVbM8daQZAUi+DPokZHraNJ+HlkbXfFk5jMyA1EQznzwivJKSBZOkYTDxyzD6Wji5Ayp3g7U90KUKG2EU5L1oaYSzn5kHmJs53P94bFaoP2UebDUDKSQNMz0o6k64anYRFvPrIDXXHvSzJeiHG/emGpsN9h6HawtcywZTKoNAiU0xd9aOuho2/0+gSxM2wjvId9Rd7UDboKHSPNzZrG5Bf7M96GdB7lWQPtZvxQ0nQVuDd6g76Zo+VgWNl6Q9PlCi403/e7mT1ifkjte+sFlNUDhTFuiSCF/Q2owE5eBoj3e/01Vq8gMrMi7QJQgbEuT7o/Yo7H/T5NxurDbBwpu2ZtMEIJyCKkPlpXrPG3ZKyyEjBYalmvn4oaZ2KQbOYEsC50eDq7nG16xtJte2I992VKzrAm5KrrmoG2GBo5vMBd+UXPNIzYUkuZAXNOorXNOOQbun5rt6ZiVnBaRYg1pChgxF6CMS5H2p7ZK5q+/cITNviYSkLNNzo70Vao6Yh2Ndco4J+Cm5g6YHT1D2m2845Zo+XQs1F2BKvmuZ5KsZeONvgfRxpoIkN031iwR5f7K2d51329puMh6eLzfzERaY87AZpEEMqBff+oRUWy2LZ+Sy/S+7KALPi66DOb1woCgFGVfC0LHm/pfyv0kOnD6SIB8sbFY4vN7UXpKzITY50CXyi6DrN9/eQqrNNWB1xrFKWuJjiBmVYRZExpi0uiIwIiJgxFTInASVJXB8s7nPxSEmyTPfkOhEgnwwcW/fj012u5ErGxIyzR+88Jklq7aS0X6GK+pM0PjFpsPcfeQ0FdmZ7NptmnAW33CtpMsIBpZIyJkJw6fAqWI48alp3rxqhRm85MQWz+AvnCTIB6tLDXDJLehboky7fdoYGHlVYMvmAwGvwdul2c45p+MaLzKktoGyaeNcG8hF1+ASGW1ulsqaZq5vWSIhd5ap7VfsgJOfSv/6DiTIhwprm0nkdP441B4xwSc52zxLO36frFk+B8pOsfavpk/24gstALRNHMniGfYc8nLRNThFxZm89Q6R0ZA31wT/E1tNahObpEYACfKhyRHsHWJT7EHf/kgcZr4UtNW0WQrvtIaG0675snLaoyI5n21vg1dKboIKNdHxcMUC07Rz7BOo2mu6MidkwJm9gzIBmgT5cHCp3jzOHjDzKsJcMGxrNl8AKdmuJG0JGdK279DSAK1Nrlr7qjeILBjFnUV5Zj4h3byPIvTEpsCEWyC3CC7WQuZ408xzfKu5U30QBXsJ8uFI20yAB9cXgEfbfrY98Ge7Mm8ORu45ihovwZEz8PXrXcukFh/6EjPNA8xoVeO/BKPmDKpgL0F+sLG2efbPV8rU7pOzzF24ydnmNv7BUNt3vwlq73HTfOPRP16CfFhyBvurTS+dM6VhPSyhBPnBTmtoPGselJhllih7wB9hnpNGQFxqAAvpJx3a47G4DdoNEuTDXVwqXLnIBPuT201XzDC8WCtBXnRmbfPMpw8w+6HwCvTWdrhwxjVfVg7jsiE22sxHxUJ8WkCKJgZYbLLJY99+KSwzy0qQFz1Tusa03ztq9omZoZ1g7UKlq9bW2gYHK+D2q13rk3PkJqjBZsx80x351E5T0QkTEuRFz1ysNY/KUjMfYTGBPmmE6aKWNMJcwC171T5vX5aQaW5YCTbu+eMPVECbtcMgIdJUM+hEJ5hAn3uVGXM2TIK93/73KaWeBB4Equ2L/k1r/Y6/jicGmM3aeRQtpVxt/O5fBgnpri+DxOHmP1OEJbA3cTkuPINrkJCCUa5lcqfr4BUdbw/2RW7BPnTb6v1dxfqp1vppPx9DBAtvg6bYrCYvuLfc4LEprhq/43kgBotob/HMIV9aDvnDINk+MIiKMCmixeAWnQBjrner2e8KyZp9EP6OFoOGow9/9UHXsrghroCflm+6d/q6bbz2qKvLnNUK+47DDdNc6xMzzG3yQoBbsJ9t8uOcKg6p/Dj+DvKPKKW+DhQD/0drfb7jBkqpZcAygJEjR/q5OCLoNZ83j7MH4AimZp84HJKG2Z+Hmy8C98Cvde++CKr2uaY/Pw3NrWYkKAfpOim8iY6H0deaZpwTW6FiZ0h0uexXkFdKrQeGe1n1PeAF4ClA25//P+D+jhtqrVcDqwFmzpzZxSCpYtBqu+R58xaYfvwJGebCb0IGHN9iElYlZpq8PQkZ5hGT1Dn4N52DmsOu+V32kbqmjXYtkyAvLicqztTsc2ZB+Wao3BPUd872K8hrrb/Yk+2UUi8Cb/XnWEI4WdvMjUzuNzO1NpkA7kjfAKa3T3waRCean9z518KBv3heO9h1GEYPh9RE1zK56Cp6IibJ3EyVW2RGrjq73/t1qQDzZ++aEVprR9eL24G9/jqWEF61t3j2/jm737MttbUN9p6Ar7jl54+ON81BQvRUfBpM/DKMnA3HPnaN8Rwk/Nkm/xOlVCGmuaYcWO7HYwnRvY4Xy/aegLZ2mDbGtSwpS26CEn2TmAkFd5pfmKd3Q+ooKP8EmusCWiy/BXmt9d/7a99C+MSuwyYR25Q81zJHxkIh+soxrgNA5gT72LRboKUxIMWRLpRi8NpywAT4BLe++XLRVfhShAWyZ9jHpt1peuW0XRrYIgzo0YQIFhXnoPwsfGGi53K56Cr8wRJl2uxn/yPkXzOg92FITV4MTh/b+wFcPcG1LD7NXHgVwl8iYyDvC6Z2f3IbVBT7/S5aCfJi8LHZ4O1icwPUcLeeNNJUIwZKVByMvg5GzjH97Kv81/lQgrwYfIoPQ2Ut3H+D53LJPCkGWmSM6WefPcNvh5A2eTG4tFs5//xbNA1Jgmsmea6TmrwIFD+OzSA1eX+zWuH1rdDcYrrrRURAhDJDzUUot3kLWNzmIyLc5t2WRVtg0iiIifJN+cqrTO4WpcwxwEwrXP3FIxRkpkJ8zOX7kLe0mRqyUm4POk/HRUNKD9IMW61g0+b47vvsq9Y2+OmbDDl1ji1/v5Cro93+/COjTSoEIcJMeAT58+fhr5vgzP5utrtg0spau8kzcb7R9L4ATGTCM+gp+z9x0fCj+zzbdTvafRRe8EMa/YgItwBqL4/7fGIcXD8Vhthv1+94u7XWcKwK1pf07riOLwNH4HU+R5gvMlsPb+semmy+6OxF73QeGjhb5/2z6nhsxxeQTZsXOsqgtWvZpTba6pqIamnj0y9MZfvILE7tNAOHLJ6Ra2rxchOUCEPhEeQPH4Zlj/ds2ygLjOhm7E6lYEq+CeLuHIFSYy7evb8bHn7BPi6oNsu1/dkxf8l+l+Xr/2b6Y9u0ea3V/myzmWVWm3nY3JY7t3F7rqw1wRm3Y2ntOW+1wfrd8KdPun8/4qLhu3dCVKT99e7nYd9vmxVO17jK5wiejmdH2ZQyedljozz30XH6ZLX5InV/X93Pw7GsaBykJ5v9ux/Tsa1NQ2s7nKqxfyiOwE/naeBUUxvlM8dRkj608/uQktt5mRBhIDyC/KRJ8M5vTK7n7uSkuwZr7q/cDNPc4azhK2fF32PZqEzPBFj9MaGHwWjFImixp0F1VlCV57xSEGlx1agDZK17jdqP8uyPam/HS8nx67GFCJTwCPLx8TBpHMSeGtjj3nPdwB6vNywWiA/hgbYHUoRFboISYSs8grwISY4a/Km6Zo95f9foO+0/OdvckShEGJIulEIMyQt0CYTwG6nJi4Bx1KgHqgbfpbTR3W8jRIiSmrwY3GISzbixQoQpqcmLgAtYDR4gY4L0jxdhTWryIuis3XnS2YTjPu0XWYX+27cQQUCCvBgUvH5ZDJ8MCemBKZAQA0Saa0TQ6Nil8hebDtPabvNY15+mnXOWYTB0DLQ2meA+9sZ+lliI4CdBXoQ2FQG6i1xEMUk8v62Oxogk/lYzkX06n41bTUqLNcvnDGAhhQgcCfIiaHjrUumcnjkS4oaYGnhChnnEp5vRnHb/DiIi7cvTISHTTEfF8reyrQDs07WBOSkhAkyCvAheafnsjM2mLmIIi69Z1PVdqTOWdrkLR419yaqtHvNCDBYS5EVwUArih0LiMBYvWQCJwyA5ix9PlXQDQvSHBHnhf0pBXBokDTPBu7UJao9C0gj7Y5hpYvHjCPZSgxeDlQR50X8RFtMWHpNsgndChgnmjqDuNYAvCEhRhRhsJMiL3omKNYE7cRgkZprn+KEm0Gtterr4cbxKIUTvSJAXPTNuIQy9wtTWu0oDoBQoCfBCBBMJ8sK72GSTuCtxuHlOGy05XoQIQRLkhT0T4wgTzB3P0QmBLlW3pFukEN2TID/Y5c+DvLmBLkWPSWAXonckyA9m8WmQlh/oUvTa/soGlqzayrZj5i5WR+B3kC8AIVwkyIcTFWG6MiZnQUsj1Bx2rYuwmJ4wKTmQkgsp2SHRJOPgCOSOwL6/siGQxREiZEiQD2Xxaa4bipJHmCDuuPX/fLl5Tsk2gT1pRFgNVj1xRDL7KxuYOCLZuaxjzV5q9EL0M8grpRYDTwITgCKtdbHbun8FvgFYgUe11n/tz7EGvZgkc0E0Odt1gTQqtuvth+SF1QDV3nLQdGymcedo0pFALwa7/tbk9wJ3AKvcFyqlJgJ3AZOALGC9Umqc1traz+MNDpHR9tp5FiRlmVp6TFKgSxV0Ogbwnn4BCDGY9CvIa60PAKjO/ae/AvxRa90CHFNKHQaKAPmfdzn510DGeHMHqfRJ96q7mnlXF2WlRi8GK38N/5cNuI+1VmFfJryJSYSMcZB7lblwKgG+T9Ysn+PRRi+E6EFNXim1HhjuZdX3tNZvdvUyL8t0F/tfBiwDGDlyZHfFCX0qwuR8SckxzTHJ2RCbIoHdRyR/vBCeug3yWusv9mG/FYD7YJw5wOku9r8aWA0wc+ZMr18EIS8hHYZNsl80HeHXlLpCCOHOX10o1wF/UEo9g7nwOhbY7qdjBb/cIhgxNdClGFSkBi+E0d8ulLcDzwMZwNtKqRKt9Y1a631KqT8B+4F24OFB17PGEmVuOkrNNdkbhRAiAPrbu+bPwJ+7WPdD4If92X9IiYyxB/WRJrAnDocIf13XFkKInpE7XnsrItJ+F2kuRMeDzWaCekKmBHUhRNCRIN+RJcqkC7hQ5Zp35HtJHWlPDyBvmxAiNAyuaBU/FC7WeC6LjHYF8JRckzJA2+DULnvOl+EynJ0QImSFb5CPsJgA7aiFJ2eb5pXdv4fIWFdQTxzmpZnFAiOvCkixhRDCl8InyEfGQPpYE8ydWRe9nN60ewe+bEIIESDhE+SzppmHEEIIJ+kOIoQQYUyCvBBChDEJ8kIIEcYkyAshRBiTIC+EEGFMgrwQQoQxCfJCCBHGJMgLIUQYkyAvhBBhTGkdPCPuKaWqgeP92EU6cM5HxQmkcDkPkHMJRuFyHiDn4jBKa53hbUVQBfn+UkoVa61nBroc/RUu5wFyLsEoXM4D5Fx6QpprhBAijEmQF0KIMBZuQX51oAvgI+FyHiDnEozC5TxAzqVbYdUmL4QQwlO41eSFEEK4kSAvhBBhLGSDvFJqsVJqn1LKppTqstuRUqpcKVWmlCpRShUPZBl7qhfnskgpdVApdVgp9fhAlrGnlFJpSqkPlFKH7M9DutguKD+X7t5jZTxnX1+qlJoeiHL2RA/O5TqlVL39MyhRSj0RiHJ2Ryn1K6XUWaXU3i7Wh9Jn0t25+P4z0VqH5AOYAFwJbARmXma7ciA90OXt77kAFuAIMBqIBvYAEwNddi/l/AnwuH36ceC/QuVz6cl7DHwJeBdQwGxgW6DL3Y9zuQ54K9Bl7cG5zAOmA3u7WB8Sn0kPz8Xnn0nI1uS11ge01gcDXQ5f6OG5FAGHtdZHtdatwB+Br/i/dL32FeA39unfALcFrii91pP3+CvAb7XxKZCqlBox0AXtgVD5e+mW1vpjoPYym4TKZ9KTc/G5kA3yvaCB95VSO5VSywJdmH7IBk66zVfYlwWbYVrrSgD7c2YX2wXj59KT9zhUPoeelnOOUmqPUupdpdSkgSmaz4XKZ9JTPv1MIn1RIn9RSq0HhntZ9T2t9Zs93M1crfVppVQm8IFS6jP7t+mA8sG5KC/LAtL/9XLn0ovdBMXn0kFP3uOg+Ry60ZNy7sLkPGlUSn0JeAMY6++C+UGofCY94fPPJKiDvNb6iz7Yx2n781ml1J8xP2MHPJj44FwqgFy3+RzgdD/32SeXOxelVJVSaoTWutL+k/lsF/sIis+lg568x0HzOXSj23JqrRvcpt9RSv1CKZWutQ61hF+h8pl0yx+fSVg31yilEpRSSY5pYCHg9ap2CNgBjFVK5SulooG7gHUBLpM364D77NP3AZ1+pQTx59KT93gd8HV7j47ZQL2jeSrIdHsuSqnhSillny7CxIOaAS9p/4XKZ9Itv3wmgb7a3I+r1LdjvsFbgCrgr/blWcA79unRmF4Fe4B9mKaRgJe9L+din/8S8Dmm10SwnstQ4EPgkP05LZQ+F2/vMbACWGGfVsDP7evLuEzPrkA/enAuj9jf/z3Ap8DVgS5zF+fxClAJtNn/n3wjhD+T7s7F55+JpDUQQogwFtbNNUIIMdhJkBdCiDAmQV4IIcKYBHkhhAhjEuSFECKMSZAXQogwJkFeCCHC2P8PqoyQejAvpqsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(X_test, Y_mean, 'r-', label='Predictive mean');\n",
    "plt.scatter(X, Y, marker='+', label='Training data')\n",
    "plt.fill_between(X_test.ravel(), \n",
    "                 Y_mean + 2 * Y_sigma, \n",
    "                 Y_mean - 2 * Y_sigma, \n",
    "                 alpha=0.5, label='Epistemic uncertainty')\n",
    "plt.title('Prediction')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ts_env",
   "language": "python",
   "name": "ts_env"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}