LR_Bayes_vs_ML.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tensorflow version:  2.0.0-rc0  running in colab?:  False\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "try: #If running in colab \n",
    "    import google.colab\n",
    "    IN_COLAB = True \n",
    "except:\n",
    "    IN_COLAB = False\n",
    "if (not tf.__version__.startswith('2')): #Checking if tf 2.0 is installed\n",
    "    if IN_COLAB: #If running in colab install tf 2.0\n",
    "        !pip install tensorflow-gpu==2.0.0-rc0\n",
    "        !pip install tensorflow-probability==0.8.0-rc0\n",
    "    print('Please install tensorflow 2.0 to run this notebook')\n",
    "print('Tensorflow version: ',tf.__version__, ' running in colab?: ', IN_COLAB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TFP Version 0.8.0-rc0\n",
      "TF  Version 2.0.0-rc0\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.model_selection import train_test_split\n",
    "import tensorflow_probability as tfp\n",
    "\n",
    "\n",
    "%matplotlib inline\n",
    "plt.style.use('default')\n",
    "\n",
    "tfd = tfp.distributions\n",
    "tfb = tfp.bijectors\n",
    "print(\"TFP Version\", tfp.__version__)\n",
    "print(\"TF  Version\",tf.__version__)\n",
    "np.random.seed(1)\n",
    "tf.random.set_seed(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6cAAAF1CAYAAAAUdZx9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df5Td510f+PcHRSHDjzJhI0I0jrGbE7QNqFhUJwsr6AYTUPhRrKgLhKUcoJwa9pCz7S4VK5NdSH9Q6VSwaffAQk2bDZQUCBtnYnBAJIjTnGbJEoVJYptEbQiBeOwQkyAIZOr6x7N/zIw9Gt8ZzY/7vd/74/U6R0dzv/fOvZ974+iZ9zzP83mqtRYAAADo06f1XQAAAAAIpwAAAPROOAUAAKB3wikAAAC9E04BAADonXAKAABA74RTAAAAeiecQs+q6qer6n8f8nN+V1X9h2E+JwCwO8Zj2B3hFPapqj5cVS/d6/e31r6vtfZPhlnTblTVq6vq5/t6fQDoyn7H6LXnGEnANB6DcAqdqqpn9F0DAABMAuEU9qGq/m2SG5P8SlX9RVX9YFW1qvqeqvqjJJfWHvfLVfXRqvqzqnp7VX3Rhud4XVX907WvX1JVD1TVD1TVx6rqoar67h3U8V9V1d1V9edV9TtJXrDp/n9ZVR9Zu//dVfWVa9dfluSHknzrWv3vXbv+3VX1/qr6ZFV9qKq+dzifGACMxhZj9JdV1f9bVVer6r1V9ZINj/+utTHvk1X1B1X17VX115L8dJIvX3uOq9d5TeMx7INwCvvQWvuOJH+U5G+11j4ryRvW7vrvkvy1JCfXbv9akhcm+bwkv5vk9ds87ecn+ZwkC0m+J8lPVtWzr1PKTyb5z0mel+Tvrv3Z6F1JbknyuUn+XZJfrqpntdZ+Pck/S/JLrbXPaq19ydrjP5bkG5P8lSTfneQ1VfWl16kBAMbGgDH69UnuSfJPszoe/sMkb6yqQ1X1mUn+zyRf11r77CT/bZL3tNben+T7kvz22jg5f52XNR7DPgin0I1Xt9b+srW2kiSttde21j7ZWnskyauTfElVfc4W3/tokn/cWnu0tfaWJH+R5MhWL1RVB5L87SQ/vPaa9yX52Y2Paa39fGvt4621x1prP57k07d7ztbaPa2132+r/n2S30jylTt98wAwhv5Okre01t7SWnuitfbWJJeTfP3a/U8k+eKqmmutPdRau383T248hv0TTqEbH1n/oqoOVNX5qvr9qvrzJB9eu+s5W3zvx1trj224/akkn7XNax1K8oyNr5nkDzc+oKr+4dqyoD9bW5L0Odu8fqrq66rqnVX1ibXHf/12jweACfAFSb55bUnv1bXx7SuSPK+19pdJvjWrs6QPVdU9VfVf7/L5jcewT8Ip7F+7zrX/IcltSV6a1UHoprXrNaTXfzjJY0mev+HajetfrO1n+cEk35Lk2WtLkv5sw+tfU39VfXqSNyb5sSTPXXv8W4ZYLwCMysYx7iNJ/m1rbX7Dn89srZ1Pktbaxdba12R1Se4HkvzMgOfYjvEY9kk4hf374yR/dZv7PzvJI0k+nuQzsrqnZGhaa48nuSvJq6vqM6rqRUm+c9PrP5bVQfMZVfXDWd27su6Pk9xUVev/Hjwzq8uMHk7yWFV9XZKvHWbNADAiG8fon0/yt6rq5NqqpmetNSK8oaqeW1W3re09fSSrW2qe2PAcN1TVM7d7IeMx7J9wCvt3Lsn/trbc5r8fcP/PZXVZz3KS30vyzg5qeGVWl/5+NMnrkvzfG+67mOTXk/zHtTr+c65dcvTLa39/vKp+t7X2yST/U1abO/1pVmd+7+6gZgDo2sYx+luzupLph7Ia+D6S5ExWfx7+tCT/S5IHk3wiq40N/8e157iU5P4kH62qP7nO6xmPYR+qtZ2uVAAAAIBumDkFAACgd8IpTIiqun/tYO7Nf76979oAYFYYj6E7lvUCAADQOzOnAAAA9O4ZfRew0XOe85x200039V0GAFPi3e9+95+01g71XcckMzYDMEzbjc1jFU5vuummXL58ue8yAJgSVfWHfdcw6YzNAAzTdmOzZb0AAAD0TjgFAACgd8IpAAAAvRNOAQAA6J1wCgAAQO+EUwAAAHonnAIAANA74RQAAIDeCacAAAD0TjgFgClUVa+tqo9V1X0brr26qpar6j1rf76+zxoBYKNn9F0AACTJ4tJyLly8kgevruTw/FzOnDySU8cW+i5rkr0uyU8k+blN11/TWvux0ZcDQF8mZYwVTgHo3eLScu64696sPPp4kmT56kruuOveJBnLwXMStNbeXlU39V0HAP2apDHWsl4Aenfh4pUnB811K48+ngsXr/RU0VR7ZVW9b23Z77MHPaCqbq+qy1V1+eGHHx51fQAM0SSNscIpAL178OrKrq6zZz+V5AVJbknyUJIfH/Sg1tqdrbXjrbXjhw4dGmV9AAzZJI2xwikAvTs8P7er6+xNa+2PW2uPt9aeSPIzSV7cd00AdGuSxljhFIDenTl5JHMHD1xzbe7ggZw5eaSniqZTVT1vw82XJ7lvq8cCMB0maYzVEAmA3q03ZJiEToKToqp+IclLkjynqh5I8iNJXlJVtyRpST6c5Ht7KxBgBvXRNXeSxthqrfVdw5OOHz/eLl++3HcZAEyJqnp3a+1433VMMmMzwHBs7pqbrM5gnjt9dCyDYle2G5vNnAIAAHRsu6654xJO+z4PVTgFAADo2Lh0zd0qgI7DeagaIgEAAHRsHLrmrgfQ5asraXkqgK4H1r7PQxVOAQCAmbe4tJwT5y/l5rP35MT5S1lcWh7q849D19ztAug4zOwKpwAAwEzbbkZxWE4dW8i500ezMD+XSrIwPzfyZkjbBdCtZnA/raqzwL6ZPacAAMBMG1WzolPHFnptfnR4fi7LAwLq+t7Tzd2Ek+TxtdNdRrEH1cwpAAAw08ZhSesobLe0ePPM7oGqp31/13tQzZwCAAAzbbsZxWmyPuO51XExG2d2bz57z8Dn6DKwC6cAAMBMG7SkddTNikZlp0uL+wjslvUCAAAzbRTNirruBjxsfXQXNnMKAADMvC6bFa13A16fmR1Fc6H9ut4S4C4IpwAAAB0aVTfgYRt1d2HhFAAAoEOT1g14cWl5pDOm6+w5BQAA6NBWTYTGsRvw+hLk5asraXlqCfIo9sgKpwAAAB3qo7nQVq7XmGm7Jchds6wXAACgQ300FxpkUGOmM7/83vyjX7k/Vz/16JbHxySjWYIsnAIAAHRs1M2FBhk0K/roEy1/+qlHk6yG1UrSBnzvKJYgW9YLAAAwA3Yy+9mS1KZro1qCLJwCAADMgJ3OfrYkC/NzqbW/z50+qlsvAAAAwzGoMdN2XvOtt+QdZ28d2XJke04BAABmwObGTJ8zdzB/+V8ey6OPP32X6foRMhu/r2vCKQAAwIzY3JhpcWk5Fy5eGdild/0ImVGFU8t6AQAAZtSpYwt5x9lbn9YEad0ojpBZJ5wCAADMuK2aJY3iCJl1wikAAMCMG9QsaVRHyKyz5xQAAGDGbW6WdHh+LmdOHhnZftNEOAUAACBPb5Y0apb1AgAA0LuhhNOqem1Vfayq7ttw7XOr6q1V9Z/W/n72MF4LAACA6TOsmdPXJXnZpmtnk/xma+2FSX5z7TYAAMDMW1xazonzl3Lz2Xty4vylLC4t911S74YSTltrb0/yiU2Xb0vys2tf/2ySU8N4LQAAgEm2uLScO+66N8tXV9KSLF9dyR133TvzAbXLhkjPba09tPb1R5M8t8PXAgAAmAgXLl7JyqOPX3Nt5dHHc+Hilc4aEi0uLffaiXcnRtKtt7XWqqoNuq+qbk9ye5LceOONoygHAACgNw9eXdnV9f1an6ldD8TrM7VJxiqgdtmt94+r6nlJsvb3xwY9qLV2Z2vteGvt+KFDhzosBwAAoH+H5+d2dX2/tpupHSddhtO7k3zn2tffmeTNHb4WAAAwBjT6ub4zJ49k7uCBa67NHTyQMyePdPJ6o56p3athHSXzC0l+O8mRqnqgqr4nyfkkX1NV/ynJS9duAwAAU0qjn505dWwh504fzfzcwSevPetgd/OGo56p3athdev9ttba81prB1trN7TW/k1r7eOtta9urb2wtfbS1trmbr4AAMAUmZTlo+PikceeePLrP/3Uo50F+VHP1O7VSBoiAQAA02/cl4+OU8faUXbsXX++cXnvWxFOAQCAoTg8P5flAUF0HJaPjlvH2lEH+VPHFsYujG7WZUMkAABghozz8tFxW3I8KftAR0k4BQAAhmK90c/C/FwqycL8XM6dPjoWM3bjtuR4nIN8XyzrBYApVFWvTfKNST7WWvvitWufm+SXktyU5MNJvqW19qd91QhMp3FdPjpuS44nZR/oKAmnADCdXpfkJ5L83IZrZ5P8ZmvtfFWdXbv9v/ZQG8DInTl55Jo9p0n/M5XjGuT7YlkvAEyh1trbk2w+xu22JD+79vXPJjk10qIAejTOS45ZZeYUAGbHc1trD619/dEkzx30oKq6PcntSXLjjTeOqDSA7pmpHG9mTgFgBrXWWpK2xX13ttaOt9aOHzp0aMSVATCrhFMAmB1/XFXPS5K1vz/Wcz0A8CTLegFgdtyd5DuTnF/7+839lgMwXRaXlnXf3QfhFACmUFX9QpKXJHlOVT2Q5EeyGkrfUFXfk+QPk3xLfxUCTJfFpeVrugEvX13JHXfdmyQC6g4JpwAwhVpr37bFXV890kIAZsSFi1euOaYmSVYefTwXLl4RTnfInlMAAIB9evDqyq6u83RmTgEAAPbp8PxclgcE0cPzc52/9rTsdTVzCgAAsE9nTh7J3MED11ybO3ggZ04e6fR11/e6Ll9dSctTe10Xl5Y7fd0umDkFAADYp/WZylHNYK7Plg6arZ3Uva7CKQAAwBCcOrbwtEDYxZLbzZ2BB5nEva7CKQAAQAe6Ol5mUGfgzUax13XYhFMAAIAB9jvr2dXxMtebFR3FXtcuCKcAAACbDGPWs6vjZbbqDJwkC7r1AgAATI/tZj13aqultftdcrtVZ+B/8a235B1nb53IYJoIpwAAAE8zjFnPro6XOXVsIedOH83C/Fwqq7Ol504fndhQus6yXgAAgE22Wjq7m1nPLo+XGdQZeNIJpwAAAJucOXnkace17GXWcxpDZFeEUwAAgE26nPVkMOEUAABgALOeoyWcAgAATJn9ntHaB+EUAABgigzjjNY+OEoGAABgigzjjNY+CKcAAABTZBhntPZBOAUAAJgiW53FupszWvsgnAIAAEyRMyePZO7ggWuu7eWM1lHTEAkAAJhJk9jRdicm9YxW4RQAAJg61wuek9rRdqcm8YxW4RQAABiqvmckdxI8t+toO2mhblrYcwoAAAzNejBcvrqSlqeC4eLS8shq2MlRKpPa0XaaCacAAMDQjMMZmzsJnpPa0XaaCacAAMDQ7GZGcnFpOSfOX8rNZ+/JifOXhja7upPgOakdbaeZcAoAAAzNTmcku1z+u5PgeerYQs6dPpqF+blUkoX5uZw7fdR+0x5piAQAAAzNmZNHrmlGlAyekeyyIdFOj1KZxI6200w4BQAAhmanwbDrhkSC5+QRTgEAgH3b7fExh+fnsjwgiGpINLvsOQUAAPZlL/tHNSRiM+EUAADYl70cH6MhEZtZ1gsAAOzLXveP2hfKRsIpAADMoN3uEd2O/aMMg2W9AAAwY4Z9xqj9owyDcAoAADNmL3tEt2P/KMNgWS8AAMyYYZwxOmhZ8DvO3jqsEplBZk4BAGDGbLUXdKd7RIe9LBgS4RSYUotLyzlx/lJuPntPTpy/ZLAEgA32u0d02MuCjdsklvUCU2j9t7nrg+b6b3OT2PsCAHlqPNxrt95hLAteZ9xmXefhtKo+nOSTSR5P8lhr7XjXrwnMtu1+m2uQA4BV+zljdJhHxxi3WTeqZb1f1Vq7RTAFRmGYv82FaVRVH66qe6vqPVV1ue96gMkzzKNjjNuss6wXmDoOAocd+arW2p/0XQQwmdZnNF999/25uvJokuRZB/c279XluD2oo3Bfs7HjVMu4GsXMaUvyG1X17qq6ffOdVXV7VV2uqssPP/zwCMoBpp2DwAFgNB557Iknv/7TTz26p469XY3b49RReJxqGWejCKdf0Vr70iRfl+T7q+pvbryztXZna+14a+34oUOHRlAOMO0cBA7X5RfHwL4Nq2NvV+P2sDsKT0st46zzZb2tteW1vz9WVW9K8uIkb+/6dYHZtp8mDzADvqK1tlxVn5fkrVX1gdbak2Nza+3OJHcmyfHjx1tfRQK7N8qlo8PcK9rFuD1Oe1nHqZZx1mk4rarPTPJprbVPrn39tUn+cZevCQBszy+OYfINCqFJRnoky7j3eBin+saplnHW9bLe5yb5D1X13iS/k+Se1tqvd/yaAMAWquozq+qz17/O6i+O7+u3KmA3ttq/+I9+5f6RLh0d9x4P41TfONUyzjqdOW2tfSjJl3T5GgDArjw3yZuqKln9OeDf+cUxTJat9i9uvrauq6Wj67Ox49qBdpzqG6daxpmjZABghvjFMYyXvewR3W3Y7HLp6Lj3eBin+saplnE1im69AADAJns9XmSrsDk/d9DSUSaacAoAAD3Y6/EiW+1ffPU3fZGj1JholvUCAEAP9nq8yPX2LwqjTCrhFAAAOnC9/aT7OV7E/kWmkWW9AAAwZDvZT+p4EbiWcAoAAEO2k/2kp44t2CMKG1jWCwAAQ7bVvtHlqytZXFq+Zn+oMAqrhFMAABiyrfaTJskdd92bZHwaF+3lrFXogmW9AAAwZIP2k67byXExo7LXs1ahC8IpAAAM2fp+0q1c77iYUdnrWavQBeEUAAA6cOrYQha2OBZmJ8fFjMJ2e2Nh1IRTAABmxuLSck6cv5Sbz96TE+cvdb58ddDy3oMHKn/5yGMjq2E7W4XkSiztZeSEUwAAJs5eQmYf+ys3Hxfz7M84mLTk6sqjY7HH88zJI6kB11tiaS8jJ5wCADBR9hoy+9pfeerYQt5x9tb8wflvyGc88xl59Ik28hq2q61tcd+47ItldginAABMlK1C5g+84b3bzqRuFbZGGcLGoYbNxn1fLLNDOAUAYKJsFeQeb23bmdStwtYoQ9g41LDZoH2xcwcP5MzJIz1VxKwSTgEAmCg7CXKDlsqOQwgbhxo227gvNkkOVD35+WmKxCgJpwAATJRBAW+QzTOsm5sTLczP5dzpozl1bKGjSp9uHGrYqq71z/XxtroLte9mTcyeZ/RdAAAA7MZ6kLtw8UoevLqST6t6MlBtNGiG9dSxhbEIgn3XMMh2DaPGsV6mj3AKAMDE2Rjw1rv3bgxWfS+VnUTj2KyJ2WJZLwAAE21cl8pOmnFs1sRsMXMKAMDEG9elspPkzMkjZqDplXAKAAA8bS/v4fm5nDl5ROhnZIRTAAAgiRlo+mXPKQAAAL0zcwoAABNkcWnZ0lumknAKAAC71FdA3HxszvLVldxx1725/IefyG994OGR1LPf9y5csxXhFAAAdmGrgJik85B14eKVa7rpJsnKo4/n9e/8o7S129vVM4xguZ/33udnx/iz5xQAAHZhq4B44eKVzl/7wasrA6+3TbcH1bMeDJevrqTlqWC4uLS849ff73vv87Nj/AmnAACwC1sFxK2uD9Ph+bkdP3ZzPcMIhvt9731+dow/4RQAAHZhq4C4m+C4V2dOHsncwQPXXKstHru5nmEEw/2+9z4/O8afcAoAALswKCDOHTyQMyePdP7ap44t5Nzpo1mYn0slWZify7d/2Y07qmcYwXC/773Pz47xpyESAADswnrjnr46zp46tvC01zr+BZ973XrOnDxyTTOiZPfBcL/vve/PjvFWrW3ePt2f48ePt8uXL/ddBgBToqre3Vo73ncdk8zYzKwYp+NNuqxlp889Tp8H02W7sdnMKQAAM22cjjfpupZBs66jrgG2Ys8pAAAzbZyONxmHWsahBmaTcAoAwEwbp+NNxqGWcaiB2SScAgAw08bpeJNxqGUcamA2CacAAMy0cTreZBxqGYcamE0aIgEAMNPG6XiTcahlHGpgNjlKBoCp5SiZ/TM2AzBM243NlvUCwIypqpdV1ZWq+mBVne27HgBILOsFgJlSVQeS/GSSr0nyQJJ3VdXdrbXf67cyJtni0rIloLvg84LBhFMAmC0vTvLB1tqHkqSqfjHJbUmEU/ZkcWk5d9x175PnYi5fXckdd92bJALXAD4v2JplvQAwWxaSfGTD7QfWrj2pqm6vqstVdfnhhx8eaXFMngsXrzwZtNatPPp4Lly80lNF483nNdji0nJOnL+Um8/ekxPnL2VxabnvkuiBcAoAXKO1dmdr7Xhr7fihQ4f6Locx9+DVlV1dn3U+r6dbn01evrqSlqdmkwXU2SOcAsBsWU7y/A23b1i7BntyeH5u4PWWmAEbYKvPa6vrs8BsMuuEUwCYLe9K8sKqurmqnpnkFUnu7rkmJtiZk0cyd/DAwPvMgD3doM9r7uCBnDl5pKeK+mc2mXXCKQDMkNbaY0lemeRikvcneUNr7f5+q2KSnTq2kHOnj2Zhi5k/M2DX2vh5VZKF+bmcO310ppshmU1mnW69ADBjWmtvSfKWvutgepw6tpBTxxZy89l70gbcbwbsWuufF6vOnDxyTQfjxGzyrDJzCgDAUJgBYy/MJrOu85nTqnpZkn+Z5ECSf91aO9/1awIAMHpmwNgrs8kkHYfTqjqQ5CeTfE1Wz1F7V1Xd3Vpz0DcAwJRZDxcXLl7Jg1dXcnh+LmdOHtlR6FhcWt7T9wHTo+uZ0xcn+WBr7UNJUlW/mOS2JMIpAMAU2ssM2Po5l+szrutdftefj+nkFxJs1nU4XUjykQ23H0jy32x8QFXdnuT2JLnxxhs7LgcAgFHZafjY6pzLH3jDe5OMJqAKSqPlFxIM0ntDpNbana21462144cOHeq7HAAAhmA9fCxfXUnL9meebtXN9/HWRnJO6m5qZTi2+oWEY4dmW9fhdDnJ8zfcvmHtGgAAU2w34WO7br4rjz6ef/BL78mJ85c6C4uC0uht9QsJxw7Ntq7D6buSvLCqbq6qZyZ5RZK7O35NAAB6tpvwcebkkcwdPLDt83U5mykojZ5jhxik03DaWnssySuTXEzy/iRvaK3d3+VrAgDQv92Ej/VzLg9UbfucXc1mCkqjN+gXEo4dovM9p621t7TWvrC19oLW2o92/XoAAPTveuFjcWk5J85fys1n78mJ85eSJD/+LV9y3RnULmYzBaXRW/+FxML8XCrJwvxczp0+qhnSjOu6Wy8AADNouzNPt+rUeu700Zw7fTQXLl7J8hYhtIvZzP2cz8re7eXYIaabcAoAQCe2Ch/bNSB6x9lbBwbYpNvZTEEJ+iecAgAwUjtpQGQ2E2aPcAoAwEgdnp8buGx385Jds5kwWzpviAQAABtpQAQMYuYUAICRsmR3dxaXln1WzAThFACAkbNkd2e26mycxOfH1LGsFwAAxtR2nY1h2ginAAAwpnbS2RimhXAKAABjanMH4+tdh0kmnAIAwJjS2ZhZoiESAACMKZ2NmSXCKQAA7EBfR7robMysEE4BAOA6HOkC3bPnFAAArsORLtA9M6cAAHSqr+Www+RIF+iemVMAADqzvhx2+epKWp5aDru4tNx3abviSBfonnAKAEBnpmU5rCNdoHuW9QIA0JlpWQ7rSBfonnAKAEBnDs/PZXlAEJ3E5bCOdIFuWdYLAEBnLIcFdsrMKQAAnbEcFtgp4RQAgE5ZDgvshGW9AAAA9E44BQAAoHdTt6x3cWnZngYAYCL5OQaYZVMVTheXlnPHXfc+edDz8tWV3HHXvUniH3YAYKz5OQaYdVO1rPfCxStP/oO+buXRx3Ph4pWeKgIA2Bk/xwCzbqrC6YMDDnje7joAwLjwcwww66YqnB6en9vVdQCAceHnGGDWTVU4PXPySOYOHrjm2tzBAzlz8khPFQEA7Mws/hyzuLScE+cv5eaz9+TE+UtZXFruuySgR1PVEGm9WYAudwDwdFX16iR/L8nDa5d+qLX2lv4qYqNZ+zlGAyhgs6kKp8nqP2b+QQOALb2mtfZjfRfBYLP0c8x2DaBm5TMYFkcQMS2mLpwCADD+NIAaDjPQTJOp2nMKAFzXK6vqfVX12qp69qAHVNXtVXW5qi4//PDDgx4C+6YB1HA4gohpIpwCwBSpqrdV1X0D/tyW5KeSvCDJLUkeSvLjg56jtXZna+14a+34oUOHRlg9s2QWG0B1wQw008SyXgCYIq21l+7kcVX1M0l+teNyGLFx33u4ub6//TcW8lsfeHhs650Eh+fnsjwgiJqBZhIJpwAwI6rqea21h9ZuvjzJfX3Ww3DtZu9hHyF2UH1vfPdyzp0+KpDuw5mTR675XBMz0Ewuy3oBYHb886q6t6rel+SrkvzPfRfE8Ox07+F6SFy+upKWp0Js12eM2hvZjVPHFnLu9NEszM+lkizMzwn8TCwzpwAwI1pr39F3DXRnp3sP+zrCxd7I7szSEURMNzOnAABTYKfdb/sKibrzAtcjnAIATIGddr/tKyTqzgtcj3AKADAFdrr3sK+QaG8kcD32nAIATImd7D1cv7+PI2fsjQS2I5wCAMwYIREYR8IpAMCI9XHO6CTWBMwW4RQAYITWzxldP85l/ZzRJL2FwXGsCZg9GiIBAIzQdueM9mUcawJmj3AKADBCfZ0zup1xrAmYPcIpAMAI9XXO6HbGsSZg9ginAAAj1Nc5o9sZx5qA2aMhEgDACA3rnNFhdtft8+xTgHWdhdOqenWSv5fk4bVLP9Rae0tXrwcAMI62CpH7CX5ddNd19inQt66X9b6mtXbL2h/BFACYKeshcvnqSlqeCpGLS8v7el7ddYFpZM8pAEBHugqRuusC06jrcPrKqnpfVb22qp7d8WsBAIyVrkKk7rrANNpXOK2qt1XVfQP+3Jbkp5K8IMktSR5K8uNbPMftVXW5qi4//PDDgx4CADCRugqRO+muu7i0nBPnL+Xms/fkxPlL+15KDNC1fTVEaq29dCePq6qfSfKrWzzHnUnuTJLjx4+3/dQDADBOzpw8ck3jomQ4R7Rcr7tuFzdfEcoAAA6ASURBVA2TALrWZbfe57XWHlq7+fIk93X1WgAA46jLI1q266673V5X4RQYV12ec/rPq+qWJC3Jh5N8b4evBQAwlvo4okXDJGASdRZOW2vf0dVzAwCwtcPzc1keEEQ1TALGmaNkAACmzE4aJgGMmy6X9QIA0IMu97oCdEU4BQCYQn3sdQXYD8t6AQAA6J2ZUwBgZiwuLVvqCjCmhFMAYCYsLi3njrvuffL8z+WrK7njrnuTREAFGAOW9QIAM+HCxStPBtN1K48+ngsXr/RUEQAbCacAwEx4cMC5n9tdB2C0hFMAYCYcnp/b1XUARks4BQBmwpmTRzJ38MA11+YOHsiZk0d6qgiAjTREAgBmwnrTo0nr1qvDMDArhFMAYOLtNMCdOrbwtOvjHP50GAZmiXAKAEy0/QS43XxvHyF2uw7Dwikwbew5BQAm2n6OiNnp966H2OWrK2l5KsQuLi3vu/7t6DAMzBLhFACYaPsJcDv93r7OSNVhGJglwikAMNF2EuAWl5Zz4vyl3Hz2npw4f+nJGc+dhr++ZjB1GAZmiXAKAEy06wW47Zbk7jT89TWDeerYQs6dPpqF+blUkoX5uZw7fdR+U2AqaYgEAEy0rY6ISZIT5y9lecDs5vqS3HecvXXg924Of2dOHrmmcVIyuhnMQR2GAaaRcAoATLzNAW5zF95B1pfk7iT8TeoZqQCTRDgFAKbOoAZGm+12Sa4ZTIBu2XMKAFOkqr65qu6vqieq6vim++6oqg9W1ZWqOtlXjaNwvUZFmgoBjB/hFACmy31JTid5+8aLVfWiJK9I8kVJXpbk/6qqA0//9umw3ayopkIA48myXgCYIq219ydJVW2+67Ykv9haeyTJH1TVB5O8OMlvj7K+xaXlkezb3KqBkVAKML6EUwCYDQtJ3rnh9gNr156mqm5PcnuS3HjjjUMrYHOTovUjXZIMPTBqYAQweYRTAJgwVfW2JJ8/4K5XtdbevN/nb63dmeTOJDl+/Hjb7/OtG9SkaP1Ily5CY1cNjPYy+zuqGePrGZc6AAYRTgFgwrTWXrqHb1tO8vwNt29YuzYyWzUpul7zonGyl9nfUc4Yb2dc6gDYioZIADAb7k7yiqr69Kq6OckLk/zOKAvYqknRbo906dN2s7/D/J4ujEsdAFsRTgFgilTVy6vqgSRfnuSeqrqYJK21+5O8IcnvJfn1JN/fWtv+INAhO3PySOYOXtsgeBRHuiwuLefE+Uu5+ew9OXH+UhaX9j5hvJfZ33GZMR6XOgC2YlkvAEyR1tqbkrxpi/t+NMmPjraip/TRpGjYS1kPz89leUCY2272dy/fs5297hsddh0Aw2bmFAAYmVPHFvKOs7fmD85/Q95x9tbO9zoOeynrXmZ/hzljvB62l6+upOWpsL2T2eC+Zq4Bdko4BQCm1rCXsp46tpBzp49mYX4ulWRhfu66Z6fu5Xu2sp+wPcw6ALpgWS8AMLW6WMq6lyNqhnWszX7DdlfH6wAMg5lTAGBqTdtS1mnoeAywFeEUAJha07aUddrCNsBGlvUCAFNt81LW9aNlRtUxeJj66HgMMCrCKQAwM4Z9tEwf7BsFppVlvQDAzBj20TIADI9wCgDMjGEfLQPA8AinAMDM0O0WYHwJpwDAzNDtFmB8aYgEAMwM3W4BxpdwCgDMlGnodru4tCxgA1NHOAUAmCDTcBwOwCD2nAIATBDH4QDTSjgFAJggjsMBppVwCgAwQRyHA0wr4RQAYII4DgeYVhoiAQBjaxhdaaets63jcIBpJZwCAGNpGF1pp7Wz7TQchwOwmWW9AMBYGkZXWp1tASaHcAoAjKVhdKXV2RZgcginAMBYGkZXWp1tASaHcAoAjKVhdKXV2RZgcuwrnFbVN1fV/VX1RFUd33TfHVX1waq6UlUn91cmADBrTh1byLnTR7MwP5dKsjA/l3Onj+6qEdAwngOA0dhvt977kpxO8q82XqyqFyV5RZIvSnI4yduq6gtba48//SkAAAYbRldanW0BJsO+Zk5ba+9vrQ1qd3dbkl9srT3SWvuDJB9M8uL9vBYAAADTq6tzTheSvHPD7QfWrj1NVd2e5PYkufHGGzsqBwBg5xaXlnPh4pU8eHUlh+fncubkEbOvAB27bjitqrcl+fwBd72qtfbm/RbQWrszyZ1Jcvz48bbf5wMA2I/FpeXccde9T56Punx1JXfcdW+SCKgAHbpuOG2tvXQPz7uc5Pkbbt+wdg0AYKxduHjlyWC6buXRx3Ph4hXhFKBDXR0lc3eSV1TVp1fVzUlemOR3OnotAIChefDqyq6uAzAc+z1K5uVV9UCSL09yT1VdTJLW2v1J3pDk95L8epLv16kXAJgEh+fndnUdgOHYb7feN7XWbmitfXpr7bmttZMb7vvR1toLWmtHWmu/tv9SAQC6d+bkkcwdPHDNtbmDB3Lm5JGeKgKYDV116wUAmEjr+0p16wUYLeEUAOjUJB7LcurYwtjXCDBthFMAoDOOZQFgp7rq1gsAsO2xLACwkZlTAKAzs3YsyyQuYQYYF1MdTg0QANCvw/NzWR4QRNePZZmmsdoSZoD9mdplvesDxPLVlbQ8NUAsLi33XRoAzIztjmWZtrHaEmaA/ZnacGqAAGAWVdU3V9X9VfVEVR3fcP2mqlqpqves/fnpUdRz6thCzp0+moX5uVSShfm5nDt9NKeOLUzdWD1rS5gBhm1ql/UaIACYUfclOZ3kXw247/dba7eMuJ4tj2WZtrH6ekuYAdje1M6cbjUQGCAAmGattfe31iZi6nHaxurtljADcH1TG04NEADwNDdX1VJV/fuq+sq+i5m2sXq7JcwAXN/ULutdHwimpQMgAKyrqrcl+fwBd72qtfbmLb7toSQ3ttY+XlV/I8liVX1Ra+3PBzz/7UluT5Ibb7xxWGU/zTSO1VstYQbg+qY2nCYGCACmU2vtpXv4nkeSPLL29bur6veTfGGSywMee2eSO5Pk+PHjbX/Vbs9YDcC6qV3WCwA8paoOVdWBta//apIXJvlQv1UBwFOEUwCYIlX18qp6IMmXJ7mnqi6u3fU3k7yvqt6T5P9J8n2ttU/0VScAbDbVy3oBYNa01t6U5E0Drr8xyRtHXxEA7IyZUwAAAHonnAIAANA74RQAAIDeCacAAAD0TjgFAACgd8IpAAAAvRNOAQAA6F211vqu4UlV9XCSP9zi7uck+ZMRljMuZvV9J7P73mf1fSez+9697+58QWvtUMevMdWuMzZPoln9/9te+Kx2x+e1cz6rnZvGz2rLsXmswul2qupya+1433WM2qy+72R23/usvu9kdt+79w2j47+7nfNZ7Y7Pa+d8Vjs3a5+VZb0AAAD0TjgFAACgd5MUTu/su4CezOr7Tmb3vc/q+05m97173zA6/rvbOZ/V7vi8ds5ntXMz9VlNzJ5TAAAAptckzZwCAAAwpSYmnFbVP6mq91XVe6rqN6rqcN81jUpVXaiqD6y9/zdV1XzfNY1CVX1zVd1fVU9U1Ux0Kauql1XVlar6YFWd7bueUaiq11bVx6rqvr5rGbWqen5V/VZV/d7af+t/v++aRqGqnlVVv1NV71173/+o75qYTVX1A1XVquo5fdcyrmb1Z5DdmMWxey9mdczbr6o6UFVLVfWrfdcyChMTTpNcaK399dbaLUl+NckP913QCL01yRe31v56kv+Y5I6e6xmV+5KcTvL2vgsZhao6kOQnk3xdkhcl+baqelG/VY3E65K8rO8ievJYkh9orb0oyZcl+f4Z+d/8kSS3tta+JMktSV5WVV/Wc03MmKp6fpKvTfJHfdcy5mb1Z5AdmeGxey9mdczbr7+f5P19FzEqExNOW2t/vuHmZyaZmc2yrbXfaK09tnbznUlu6LOeUWmtvb+1dqXvOkboxUk+2Fr7UGvtvyT5xSS39VxT51prb0/yib7r6ENr7aHW2u+uff3JrA4+C/1W1b226i/Wbh5c+zMz/6YzNl6T5Afjv71tzerPILswk2P3XszqmLcfVXVDkm9I8q/7rmVUJiacJklV/WhVfSTJt2e2Zk43+rtJfq3vIujEQpKPbLj9QPyjPTOq6qYkx5L8f/1WMhpry5Tek+RjSd7aWpuJ9814qKrbkiy31t7bdy0Txs8gT2fs3oNZG/P24V9k9ZdoT/RdyKg8o+8CNqqqtyX5/AF3vaq19ubW2quSvKqq7kjyyiQ/MtICO3S99772mFdldUnE60dZW5d28r5h2lXVZyV5Y5J/sGmVyNRqrT2e5Ja1/Wtvqqovbq3N3L5jurPd+JLkh7K6pJfM7s8g9GMWx7y9qKpvTPKx1tq7q+olfdczKmMVTltrL93hQ1+f5C2ZonB6vfdeVd+V5BuTfHWbovN/dvG/+SxYTvL8DbdvWLvGFKuqg1kdpF/fWrur73pGrbV2tap+K6v7joVThmar8aWqjia5Ocl7qypZ/bf2d6vqxa21j46wxLExqz+DDImxexdmfczbpRNJvqmqvj7Js5L8lar6+dba3+m5rk5NzLLeqnrhhpu3JflAX7WMWlW9LKtT+t/UWvtU3/XQmXcleWFV3VxVz0zyiiR391wTHarVn4z/TZL3t9b+j77rGZWqOrTe8bOq5pJ8TWbo33T61Vq7t7X2ea21m1prN2V1GeaXzmowvR4/g1yXsXuHZnXM26vW2h2ttRvW/p16RZJL0x5MkwkKp0nOV9V9VfW+rC7FmaX20z+R5LOTvHXtKJ2f7rugUaiql1fVA0m+PMk9VXWx75q6tNZw4pVJLma1ScAbWmv391tV96rqF5L8dpIjVfVAVX1P3zWN0Ikk35Hk1rX/b79n7Tek0+55SX5r7d/zd2V1z+lMtMiHCTSTP4Ps1KyO3Xs0q2Meu1BWZwAAANC3SZo5BQAAYEoJpwAAAPROOAUAAKB3wikAAAC9E04BAADonXAKAABA74RTAAAAeiecAgAA0Lv/Hx4Yr+SXAEhUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sigma = 3\n",
    "N_train = 4\n",
    "N_test  = 100\n",
    "def get_data():\n",
    "    x=np.arange(-5,5,0.01)#x range\n",
    "    noise=np.random.normal(loc=0,scale=sigma,size=len(x))# const var\n",
    "    y=2.7*x+noise #y linear model with const var\n",
    "    y=y.reshape((len(y),1))\n",
    "    x=x.reshape((len(x),1))\n",
    "    train_idx=np.random.choice(range(0,len(x)),N_train,replace=False)# sample 4 trainpoints\n",
    "    test_range=np.delete(range(0,len(x)),train_idx)\n",
    "    test_idx=np.random.choice(test_range,N_test,replace=False)# sample 100 testpoints, without the 4 from train\n",
    "    return x[train_idx], y[train_idx], x[test_idx], y[test_idx]\n",
    "    \n",
    "x_train, y_train, x_test, y_test = get_data()\n",
    "\n",
    "plt.figure(figsize=(16,6))\n",
    "plt.subplot(1,2,1)\n",
    "plt.scatter(x_train, y_train)\n",
    "plt.title(\"train_data\")\n",
    "plt.subplot(1,2,2)\n",
    "plt.scatter(x_test, y_test)\n",
    "plt.title(\"test_data\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "model_skl = LinearRegression()\n",
    "\n",
    "def get_nll_lin(x_train, y_train,x_test,y_test):\n",
    "    res = model_skl.fit(x_train, y_train)\n",
    "    r2 = (y_test - model_skl.predict(x_test))**2\n",
    "    nll = N_test*np.log(np.sqrt(2*np.pi*sigma**2)) + np.sum((r2/(2*sigma**2)))\n",
    "    return nll,res.coef_,res.intercept_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def posterior(Phi, t, alpha, beta, return_inverse=False):\n",
    "    \"\"\"Computes mean and covariance matrix of the posterior distribution.\"\"\"\n",
    "    S_N_inv = alpha * np.eye(Phi.shape[1]) + beta * Phi.T.dot(Phi)\n",
    "    S_N = np.linalg.inv(S_N_inv)\n",
    "    m_N = beta * S_N.dot(Phi.T).dot(t)\n",
    "\n",
    "    if return_inverse:\n",
    "        return m_N, S_N, S_N_inv\n",
    "    else:\n",
    "        return m_N, S_N\n",
    "\n",
    "\n",
    "def posterior_predictive(Phi_test, m_N, S_N, beta):\n",
    "    \"\"\"Computes mean and variances of the posterior predictive distribution.\"\"\"\n",
    "    y = Phi_test.dot(m_N)\n",
    "    # Only compute variances (diagonal elements of covariance matrix)\n",
    "    y_var = 1 / beta + np.sum(Phi_test.dot(S_N) * Phi_test, axis=1)\n",
    "    \n",
    "    return y, y_var\n",
    "\n",
    "def expand(x, bf, bf_args=None):\n",
    "    if bf_args is None:\n",
    "        return np.concatenate([np.ones(x.shape), bf(x)], axis=1)\n",
    "    else:\n",
    "        return np.concatenate([np.ones(x.shape)] + [bf(x, bf_arg) for bf_arg in bf_args], axis=1)\n",
    "\n",
    "def identity_basis_function(x):\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_nll_bayes(x_train, y_train,x_test, y_test):\n",
    "    beta =  1.0 / sigma**2\n",
    "    alpha = 1.0 #Prior\n",
    "\n",
    "    # Mean and covariance matrix of posterior\n",
    "    Phi = expand(x_train, identity_basis_function) #Designmatrix\n",
    "    m, S = posterior(Phi, y_train, alpha, beta)\n",
    "    Phi_test = expand(x_test, identity_basis_function)\n",
    "    # Mean and variances of posterior predictive \n",
    "    yt, yt_var = posterior_predictive(Phi_test, m, S, beta)\n",
    "    nll = 0\n",
    "    for i in range(N_test):\n",
    "        n = tfd.Normal(loc=yt[i], scale=np.sqrt(yt_var[i])) #TODO make for all N_train at once\n",
    "        nll += n.log_prob(y_test[i])\n",
    "    return -nll,m[1],m[0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficients a_ml:  [2.54714062] Coefficients b_ml [0.53985247]  nll (per data point) 2.4814315448019544\n",
      "Coefficients a_bayes:  2.5171243658426117 Coefficients b_ml [0.49328651]  nll (per data point) tf.Tensor([2.48321624], shape=(1,), dtype=float64)\n"
     ]
    }
   ],
   "source": [
    "N_train = 100 #Large number Bayes should be similar to MaxLike\n",
    "x_train, y_train, x_test, y_test = get_data()\n",
    "\n",
    "nll,a_ml,b_ml = get_nll_lin(x_train, y_train, x_test, y_test)\n",
    "print('Coefficients a_ml: ', a_ml[0], 'Coefficients b_ml', b_ml, ' nll (per data point)', nll/N_test)\n",
    "\n",
    "nll,a_ml,b_ml = get_nll_bayes(x_train, y_train, x_test, y_test)\n",
    "print('Coefficients a_bayes: ', a_ml[0], 'Coefficients b_ml', b_ml, ' nll (per data point)', nll/N_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficients a_ml:  [3.31185201] Coefficients b_ml [-0.48170139]  nll (per data point) 2.731037956960999\n",
      "Coefficients a_bayes:  1.8095932661150993 Coefficients b_ml [-0.34235386]  nll (per data point) tf.Tensor([2.99381802], shape=(1,), dtype=float64)\n"
     ]
    }
   ],
   "source": [
    "N_train = 4 #Small Bayes\n",
    "x_train, y_train, x_test, y_test = get_data()\n",
    "\n",
    "nll,a_ml,b_ml = get_nll_lin(x_train, y_train, x_test, y_test)\n",
    "print('Coefficients a_ml: ', a_ml[0], 'Coefficients b_ml', b_ml, ' nll (per data point)', nll/N_test)\n",
    "\n",
    "nll,a_ml,b_ml = get_nll_bayes(x_train, y_train, x_test, y_test)\n",
    "print('Coefficients a_bayes: ', a_ml[0], 'Coefficients b_ml', b_ml, ' nll (per data point)', nll/N_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_train = 4\n",
    "\n",
    "runs = 100\n",
    "nll_bayes = np.zeros(runs)\n",
    "nll_max_like = np.zeros(runs)\n",
    "for i in range(runs):\n",
    "    x_train, y_train, x_test, y_test = get_data()\n",
    "    m,_,_= get_nll_lin(x_train, y_train, x_test, y_test)\n",
    "    nll_max_like[i] = m\n",
    "    b ,_,_ = get_nll_bayes(x_train, y_train, x_test, y_test)\n",
    "    nll_bayes[i] = b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Means  277.52757920461846 369.22127207394334\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(269.0087899282315, 273.9339284936441)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA54AAAF1CAYAAACAkzWPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df5Sld10n+PfHNAH5IUlIG7MJPR2HAJNlJDBlBgZGIQEMhDWZM6hhAVvE07uICKiDjZ4jjmc9G9QRcXZGtxcCETGEH0EyBB0yMch6FgKdECAhQAI0kJCkGyH8UoHAZ/+4T4eboqp/VNVT91b363VOnXp+3brv+va9/Tzv+zz3VnV3AAAAYCzfN+sAAAAAHN4UTwAAAEaleAIAADAqxRMAAIBRKZ4AAACMSvEEAABgVIonAACsUFU9q6reNTXfVfWQWWaCeaR4whqrqt1V9Y9V9bWq+lJVXV5VD551LgA4nA37329W1fGLln9wKINbV/nzlyyU3f2G7n7Kan42HAkUTxjH/9Ld909yYpI7kvznGecBgCPBp5M8c99MVf3LJPedXRxgH8UTRtTd/5TkLUlOS5KqOmd45fUrVfW5qvrtfdsOZ0ZfOH37qvpwVf27YfrhVXVFVX2xqj5eVT89td3TquqjVfXVqrq1qn5tXX5BAJgvr0/ys1Pz25L82b6ZA+yHf6aqPl1VPzDMP7Wqbq+qzfu7w6r6uar6u2XWPX64nycM88vuy+Fwp3jCiKrqvkl+Jsn7hkVfz2SHeEySc5I8v6rOG9ZdlOTZU7d9ZJKTklxeVfdLckWSv0jyg0nOT/Jfq+q0YfPXJPnfuvsBSR6R5G/G/L0AYE69L8kPVNW/qKqjMtlf/vnU+mX3w919SZL/L8kfV9WDMtm3/kJ3711JkKo6O8nFSf59d7/7IPblcFhTPGEcf1lVdyb5cpInJ/n9JOnud3f3R7r7O9394Ux2SD8+3OayJA+tqlOH+eckuaS7v5nk6Ul2d/dru/uu7v5gkrcm+alh228lOa2qfqC7v9Td167LbwkA82ffWc8nJ7kxya37VhxgP5wkL0hyZpJ3J/lv3f2OFWb4qST/d5Kndvf7h2UH2pfDYU3xhHGc193HJLlPkl9K8rdV9UNV9a+r6qqq2ltVX07yvyc5Prn7stxLkjy7qr4vk/eovH74ef8syb+uqjv3fSV5VpIfGtb/+yRPS/KZqvrbqnrsev2iADBnXp/kf03yc5m6zDZJ9rcfTpLuvjPJmzO5eug/rSLDi5O8qbuvn1p2oH05HNYUTxhRd3+7uy9N8u0kj8/k8prLkjy4ux+Y5E+T1NRNLspkJ3RWkn/o7vcOyz+X5G+7+5ipr/t39/OH+/lAd5+byaU7f5nkTevx+wHAvOnuz2TyIUNPS3LpotX73Q9X1elJfj6TM6F/vIoYP5XkvKp60dSy/e7L4XCneMKIauLcJMdmcrnPA5J8sbv/qarOyOQV2bsNRfM7mbzK+vqpVe/I5DLc51TVvYavHx3ew3L08DfEHtjd30ryleFnAMCR6nlJzuzury9avux+uKruk8n7QX8jyXOTnFRVv7jo9kdX1X2mvo5a5v4/n8mLyC+qqn3Fctl9+ep+VdgYFE8Yx3+rqq9lUgJ/N8m27r4hyS8m+Z2q+mqS38rSZyb/LMm/zNSHIXT3V5M8JZMPIvh8ktuTvCLJvYdNnpNkd1V9JZPLhp41xi8FABtBd3+yu3ctsWp/++H/M8nnuvtPuvsbmXzg3/8x9dkLSXJDkn+c+nrufjJ8NpPyuaOqfuEg9uVwWKvunnUGYEpV/WyS7d39+FlnAQCAteCMJ8yR4c+v/GKSnbPOAgAAa0XxhDlRVT+RZG+SOzL58AMAADgsuNQWAACAUTnjCQAAwKgUTwAAAEa1aT3v7Pjjj++tW7eu510CcBi75pprvtDdm2edYyOzbwZgLS23b17X4rl169bs2rXUn1QCgENXVZ+ZdYaNzr4ZgLW03L7ZpbYAAACMSvEEAABgVIonAAAAo1I8AQAAGJXiCQAAwKgUTwAAAEaleAIAADAqxRMAAIBRKZ4AAACMSvEEAABgVIonAAAAo1I8AQAAGJXiCQAAwKg2zTrAetm64/JD2n73BeeMlAQADj+Hup9N7GsBjiTOeAIAADAqxRMAAIBRKZ4AAACMSvEEAABgVIonAAAAo1I8AQAAGJXiCQAAwKgUTwAAAEaleAIAADAqxRMAAIBRKZ4AAACMSvEEgA2oql5SVTdU1fVVdXFV3aeqTqmqq6vq5qq6pKqOnnVOAEgUTwDYcKrqpCS/nGShux+R5Kgk5yd5RZJXdvdDknwpyfNmlxIAvkvxBICNaVOS76+qTUnum+S2JGcmecuw/qIk580oGwDcwwGLZ1VdWFV7qur6RctfWFUfGy7z+b3xIgIA07r71iR/kOSzmRTOLye5Jsmd3X3XsNktSU5a6vZVtb2qdlXVrr17965HZACOcAdzxvN1Sc6eXlBVT0xybpJHdvf/nMnODwBYB1V1bCb74VOS/E9J7pdF++r96e6d3b3Q3QubN28eKSUAfNcBi2d3vyfJFxctfn6SC7r7G8M2e0bIBgAs7UlJPt3de7v7W0kuTfK4JMcMl94myclJbp1VQACYtunAmyzpoUn+bVX9bpJ/SvJr3f2BpTasqu1JtifJli1bVnh362/rjssPafvdF5wzUhIA+B6fTfKYqrpvkn9MclaSXUmuSvKMJG9Msi3J22eWEACmrPTDhTYlOS7JY5L8hyRvqqpaakOX8wDA2uruqzP5EKFrk3wkk/35ziS/nuRXqurmJA9K8pqZhQSAKSs943lLkku7u5O8v6q+k+T4JD6hAADWQXe/PMnLFy3+VJIzZhAHAPZrpWc8/zLJE5Okqh6a5OgkX1irUAAAABw+DnjGs6ouTvKEJMdX1S2ZvLp6YZILhz+x8s0k24aznwAAAHAPByye3f3MZVY9e42zAAAAcBha6aW2AAAAcFAUTwAAAEaleAIAADAqxRMAAIBRKZ4AAACMSvEEAABgVIonAAAAo1I8AQAAGJXiCQAAwKgUTwAAAEaleAIAADAqxRMAAIBRKZ4AAACMSvEEAABgVIonAAAAo1I8AQAAGJXiCQAAwKgUTwAAAEaleAIAADAqxRMAAIBRKZ4AAACMSvEEAABgVIonAAAAo1I8AQAAGJXiCQAAwKgUTwDYYKrqYVV13dTXV6rqxVV1XFVdUVU3Dd+PnXVWAEgUTwDYcLr74919enefnuRfJfmHJG9LsiPJld19apIrh3kAmDnFEwA2trOSfLK7P5Pk3CQXDcsvSnLezFIBwJQDFs+qurCq9lTV9Uus+9Wq6qo6fpx4AMABnJ/k4mH6hO6+bZi+PckJS92gqrZX1a6q2rV37971yAjAEe5gzni+LsnZixdW1YOTPCXJZ9c4EwBwEKrq6CQ/meTNi9d1dyfppW7X3Tu7e6G7FzZv3jxySgA4iOLZ3e9J8sUlVr0yyUuzzE4NABjdU5Nc2913DPN3VNWJSTJ83zOzZAAwZUXv8ayqc5Pc2t0fOohtXc4DAON4Zr57mW2SXJZk2zC9Lcnb1z0RACzhkItnVd03yW8k+a2D2d7lPACw9qrqfkmenOTSqcUXJHlyVd2U5EnDPADM3KYV3OafJzklyYeqKklOTnJtVZ3R3bevZTgAYGnd/fUkD1q07O8z+ZRbAJgrh1w8u/sjSX5w33xV7U6y0N1fWMNcAAAAHCYO5s+pXJzkvUkeVlW3VNXzxo8FAADA4eKAZzy7+5kHWL91zdIAAABw2FnRp9oCAADAwVI8AQAAGJXiCQAAwKgUTwAAAEaleAIAADAqxRMAAIBRKZ4AAACMSvEEAABgVIonAAAAo1I8AQAAGJXiCQAAwKgUTwAAAEaleAIAADAqxRMAAIBRKZ4AAACMSvEEAABgVIonAAAAo1I8AQAAGJXiCQAAwKgUTwAAAEaleAIAADAqxRMAAIBRKZ4AAACMSvEEAABgVIonAGxAVXVMVb2lqj5WVTdW1WOr6riquqKqbhq+HzvrnACQKJ4AsFG9Kslfd/fDkzwyyY1JdiS5srtPTXLlMA8AM6d4AsAGU1UPTPJjSV6TJN39ze6+M8m5SS4aNrsoyXmzSQgA93TA4llVF1bVnqq6fmrZ7w+X9ny4qt5WVceMGxMAmHJKkr1JXltVH6yqV1fV/ZKc0N23DdvcnuSEmSUEgCkHc8bzdUnOXrTsiiSP6O4fSfKJJC9b41wAwPI2JXl0kj/p7kcl+XoWXVbb3Z2kl7pxVW2vql1VtWvv3r2jhwWAAxbP7n5Pki8uWvau7r5rmH1fkpNHyAYALO2WJLd099XD/FsyKaJ3VNWJSTJ837PUjbt7Z3cvdPfC5s2b1yUwAEe2tXiP588n+as1+DkAwEHo7tuTfK6qHjYsOivJR5NclmTbsGxbkrfPIB4AfI9Nq7lxVf1mkruSvGE/22xPsj1JtmzZspq7AwC+64VJ3lBVRyf5VJLnZvKC8puq6nlJPpPkp2eYDwDutuLiWVU/l+TpSc4a3keypO7emWRnkiwsLCy7HQBw8Lr7uiQLS6w6a72zAMCBrKh4VtXZSV6a5Me7+x/WNhIAAACHk4P5cyoXJ3lvkodV1S3D5Tv/V5IHJLmiqq6rqj8dOScAAAAb1AHPeHb3M5dY/JoRsgAAAHAYWotPtQUAAIBlKZ4AAACMSvEEAABgVIonAAAAo1I8AQAAGJXiCQAAwKgUTwAAAEaleAIAADAqxRMAAIBRKZ4AAACMSvEEAABgVIonAAAAo1I8AQAAGNWmWQc4XGzdcfkh32b3BeeMkAQAAGC+OOMJAADAqBRPAAAARqV4AgAAMCrFEwAAgFEpngAAAIxK8QQAAGBUiicAAACjUjwBAAAYleIJAADAqBRPAAAARqV4AgAAMKpNsw4AABy6qtqd5KtJvp3kru5eqKrjklySZGuS3Ul+uru/NKuMALCPM54AsHE9sbtP7+6FYX5Hkiu7+9QkVw7zADBziicAHD7OTXLRMH1RkvNmmAUA7nbA4llVF1bVnqq6fmrZcVV1RVXdNHw/dtyYAMAineRdVXVNVW0flp3Q3bcN07cnOWE20QDgng7mjOfrkpy9aJlLeQBgth7f3Y9O8tQkL6iqH5te2d2dSTn9HlW1vap2VdWuvXv3rkNUAI50Byye3f2eJF9ctNilPAAwQ9196/B9T5K3JTkjyR1VdWKSDN/3LHPbnd290N0LmzdvXq/IABzBVvqptgd9Kc9w+c/2JNmyZcsK746V2Lrj8kPafvcF54yUBIC1VFX3S/J93f3VYfopSX4nyWVJtiW5YPj+9tmlBIDvWvWfU+nurqolL+UZ1u9MsjNJFhYWlt0OADhoJyR5W1Ulk335X3T3X1fVB5K8qaqel+QzSX56hhkB4G4rLZ53VNWJ3X3b/i7lAQDWXnd/Kskjl1j+90nOWv9EALB/K/1zKvsu5UlcygMAAMB+HMyfU7k4yXuTPKyqbhku37kgyZOr6qYkTxrmAQAA4Hsc8FLb7n7mMqtcygMAAMABrfRSWwAAADgoiicAAACjUjwBAAAYleIJAADAqBRPAAAARqV4AgAAMCrFEwAAgFEpngAAAIxK8QQAAGBUiicAAACjUjwBAAAYleIJAADAqBRPAAAARqV4AgAAMCrFEwAAgFEpngAAAIxK8QQAAGBUiicAAACjUjwBAAAYleIJAADAqBRPAAAARqV4AgAAMCrFEwAAgFEpngAAAIxK8QQAAGBUiicAbFBVdVRVfbCq3jHMn1JVV1fVzVV1SVUdPeuMAJAongCwkb0oyY1T869I8srufkiSLyV53kxSAcAiqyqeVfWSqrqhqq6vqour6j5rFQwAWF5VnZzknCSvHuYryZlJ3jJsclGS82aTDgDuacXFs6pOSvLLSRa6+xFJjkpy/loFAwD264+SvDTJd4b5ByW5s7vvGuZvSXLSLIIBwGKrvdR2U5Lvr6pNSe6b5POrjwQA7E9VPT3Jnu6+ZoW3315Vu6pq1969e9c4HQB8rxUXz+6+NckfJPlsktuSfLm737VWwQCAZT0uyU9W1e4kb8zkEttXJTlmeDE4SU5OcutSN+7und290N0LmzdvXo+8ABzhNh14k6VV1bFJzk1ySpI7k7y5qp7d3X++aLvtSbYnyZYtW1YRla07Lp91BADmQHe/LMnLkqSqnpDk17r7WVX15iTPyKSMbkvy9pmFBIApq7nU9klJPt3de7v7W0kuTfJvFm/kVVUAWDe/nuRXqurmTN7z+ZoZ5wGAJKs445nJJbaPqar7JvnHJGcl2bUmqQCAg9Ld707y7mH6U0nOmGUeAFjKat7jeXUmH9l+bZKPDD9r5xrlAgAA4DCxmjOe6e6XJ3n5GmUBAADgMLTaP6cCAAAA+6V4AgAAMCrFEwAAgFEpngAAAIxK8QQAAGBUiicAAACjUjwBAAAYleIJAADAqBRPAAAARqV4AgAAMCrFEwAAgFEpngAAAIxK8QQAAGBUiicAAACjUjwBAAAYleIJAADAqBRPAAAARqV4AgAAMCrFEwAAgFEpngAAAIxK8QQAAGBUiicAAACjUjwBAAAYleIJAADAqBRPAAAARqV4AgAAMCrFEwA2mKq6T1W9v6o+VFU3VNV/HJafUlVXV9XNVXVJVR0966wAkCieALARfSPJmd39yCSnJzm7qh6T5BVJXtndD0nypSTPm2FGALjbqopnVR1TVW+pqo9V1Y1V9di1CgYALK0nvjbM3mv46iRnJnnLsPyiJOfNIB4AfI/VnvF8VZK/7u6HJ3lkkhtXHwkAOJCqOqqqrkuyJ8kVST6Z5M7uvmvY5JYkJ80qHwBMW3HxrKoHJvmxJK9Jku7+ZnffuVbBAIDldfe3u/v0JCcnOSPJww/2tlW1vap2VdWuvXv3jpYRAPbZtIrbnpJkb5LXVtUjk1yT5EXd/fXpjapqe5LtSbJly5ZV3N3hZ+uOy2cdgRz6v8PuC84ZKQnAoevuO6vqqiSPTXJMVW0aznqenOTWZW6zM8nOJFlYWOh1CwvAEWs1l9puSvLoJH/S3Y9K8vUkOxZv1N07u3uhuxc2b968irsDAJKkqjZX1THD9PcneXImb3e5Kskzhs22JXn7bBICwD2tpnjekuSW7r56mH9LJkUUABjXiUmuqqoPJ/lAkiu6+x1Jfj3Jr1TVzUkelOHtMAAwayu+1La7b6+qz1XVw7r740nOSvLRtYsGACyluz+c5FFLLP9UJu/3BIC5spr3eCbJC5O8YfgD1Z9K8tzVRwIAAOBwsqri2d3XJVlYoywAAAAchlb7dzwBAABgvxRPAAAARqV4AgAAMCrFEwAAgFEpngAAAIxK8QQAAGBUiicAAACjUjwBAAAYleIJAADAqBRPAAAARqV4AgAAMCrFEwAAgFFtmnWAldi64/JZR2BwqP8Wuy84Z6QkAADAvHLGEwAAgFEpngAAAIxK8QQAAGBUiicAAACjUjwBAAAYleIJAADAqBRPAAAARqV4AgAAMCrFEwAAgFEpngAAAIxK8QQAAGBUiicAAACjUjwBAAAYleIJABtMVT24qq6qqo9W1Q1V9aJh+XFVdUVV3TR8P3bWWQEgWYPiWVVHVdUHq+odaxEIADigu5L8anefluQxSV5QVacl2ZHkyu4+NcmVwzwAzNxanPF8UZIb1+DnAAAHobtv6+5rh+mvZrIfPinJuUkuGja7KMl5s0kIAPe0quJZVScnOSfJq9cmDgBwKKpqa5JHJbk6yQndfduw6vYkJ8woFgDcw2rPeP5Rkpcm+c4aZAEADkFV3T/JW5O8uLu/Mr2uuztJL3O77VW1q6p27d27dx2SAnCk27TSG1bV05Ps6e5rquoJ+9lue5LtSbJly5aV3h0ctK07Lp91BIDRVdW9Mimdb+juS4fFd1TVid19W1WdmGTPUrft7p1JdibJwsLCkuUUANbSas54Pi7JT1bV7iRvTHJmVf354o26e2d3L3T3wubNm1dxdwBAklRVJXlNkhu7+w+nVl2WZNswvS3J29c7GwAsZcXFs7tf1t0nd/fWJOcn+ZvufvaaJQMAlvO4JM/J5EXf64avpyW5IMmTq+qmJE8a5gFg5lZ8qS0AMBvd/XdJapnVZ61nFgA4GGtSPLv73UnevRY/CwAAgMPLWvwdTwAAAFiW4gkAAMCoFE8AAABGpXgCAAAwKsUTAACAUSmeAAAAjErxBAAAYFSKJwAAAKNSPAEAABiV4gkAAMCoFE8AAABGpXgCAAAwKsUTAACAUW2adQDmx9Ydl886AgAAcBhyxhMAAIBRKZ4AAACMSvEEAABgVIonAAAAo1I8AQAAGJXiCQAAwKgUTwAAAEaleAIAADAqxRMAAIBRKZ4AAACMSvEEAABgVIonAAAAo1I8AQAAGJXiCQAbTFVdWFV7qur6qWXHVdUVVXXT8P3YWWYEgGkrLp5V9eCquqqqPlpVN1TVi9YyGACwrNclOXvRsh1JruzuU5NcOcwDwFxYzRnPu5L8anefluQxSV5QVaetTSwAYDnd/Z4kX1y0+NwkFw3TFyU5b11DAcB+rLh4dvdt3X3tMP3VJDcmOWmtggEAh+SE7r5tmL49yQmzDAMA0zatxQ+pqq1JHpXk6iXWbU+yPUm2bNmyFncHM7V1x+WHfJvdF5wzQhKApXV3V1Uvt35e9s2H+v+p/0sBNq5Vf7hQVd0/yVuTvLi7v7J4fXfv7O6F7l7YvHnzau8OAFjaHVV1YpIM3/cst6F9MwDrbVXFs6rulUnpfEN3X7o2kQCAFbgsybZheluSt88wCwDcw2o+1baSvCbJjd39h2sXCQDYn6q6OMl7kzysqm6pqucluSDJk6vqpiRPGuYBYC6s5j2ej0vynCQfqarrhmW/0d3vXH0sAGA53f3MZVadta5BAOAgrbh4dvffJak1zAIAAMBhaNUfLgQAAAD7o3gCAAAwKsUTAACAUSmeAAAAjErxBAAAYFSKJwAAAKNSPAEAABiV4gkAAMCoFE8AAABGpXgCAAAwKsUTAACAUSmeAAAAjErxBAAAYFSbZh2AI8vWHZcf0va7LzhnpCTr61B/77Gtx7geqf/W82Ylj71D/bfwbw0AHIgzngAAAIxK8QQAAGBUiicAAACjUjwBAAAYleIJAADAqBRPAAAARqV4AgAAMCrFEwAAgFEpngAAAIxK8QQAAGBUiicAAACjUjwBAAAYleIJAADAqDat5sZVdXaSVyU5Ksmru/uCNUkFAKyIffPqbN1x+awjrLvdF5wz6wjAGjvU/8vW4/+BFZ/xrKqjkvyXJE9NclqSZ1bVaWsVDAA4NPbNAMyr1Vxqe0aSm7v7U939zSRvTHLu2sQCAFbAvhmAubSa4nlSks9Nzd8yLAMAZsO+GYC5tKr3eB6MqtqeZPsw+7Wq+viiTY5P8oWxc4xgI+becJnrFUk2YO7B3OYexnU5M8l9gEwHY27H+wDmLvdB/FusKvMa/FtP+2dr+tOOEAexb55L9Yr5e74cgrnNfpDPybnNfxA2cvZkY+eXfXYOKf967JtXUzxvTfLgqfmTh2X30N07k+xc7odU1a7uXlhFjpnYiLk3YuZE7vUm9/raiLk3YuYjyJrsm+fVRn7sbeTsycbOv5GzJxs7v+yzM4/5V3Op7QeSnFpVp1TV0UnOT3LZ2sQCAFbAvhmAubTiM57dfVdV/VKS/57JR7Zf2N03rFkyAOCQ2DcDMK9W9R7P7n5nkneuMsOGu9RnsBFzb8TMidzrTe71tRFzb8TMR4w12jfPq4382NvI2ZONnX8jZ082dn7ZZ2fu8ld3zzoDAAAAh7HVvMcTAAAADmjU4llVD66qq6rqo1V1Q1W9aNH6X62qrqrjh/mqqj+uqpur6sNV9egx860kd1W9sKo+Niz/vanlLxtyf7yqfmKeclfV6VX1vqq6rqp2VdUZw/J5Ge/7VNX7q+pDQ+7/OCw/paquHvJdMnxQRqrq3sP8zcP6rXOU+Q3DY+D6qrqwqu41LJ/rsZ5a/8dV9bWp+ZmP9ZBjufGuqvrdqvpEVd1YVb88tXxux7uqzqqqa4fn5N9V1UOG5XMx3lP5j6qqD1bVO4b5uX1Ocnioqt1V9ZF9+6th2XFVdUVV3TR8P3ZYPvPn+fD//J6qun5q2SHnraptw/Y3VdW2GWb/7aq6dRj/66rqaVPrljzOqaqzh2U3V9WOdcq+3PHORhn75fLP/fjvZ792yPuH5X6nGeV/XVV9emrsTx+Wz9VjZ7jfVe+bZzH2SZLuHu0ryYlJHj1MPyDJJ5KcNsw/OJMPP/hMkuOHZU9L8ldJKsljklw9Zr5DzZ3kiUn+R5J7D+t+cPh+WpIPJbl3klOSfDLJUXOU+11Jnjo1xu+es/GuJPcfpu+V5Oohz5uSnD8s/9Mkzx+mfzHJnw7T5ye5ZI4yP21YV0kunso812M9zC8keX2Sr01tP/OxPsB4PzfJnyX5vmHdvufkXI/38Nz8F1Nj/Lp5Gu+p/L+S5C+SvGOYn9vnpK/D4yvJ7gzHBFPLfi/JjmF6R5JXDNMzf54n+bEkj05y/UrzJjkuyaeG78cO08fOKPtvJ/m1JbZd8jhn+Ppkkh9OcvSwzWnrkH25452NMvbL5Z/78c8aHbMt9zutw9gvl/91SZ6xxPZz9dgZ7ntV++ZZjX13j3vGs7tv6+5rh+mvJrkxyUnD6lcmeWmS6TeZnpvkz3rifUmOqaoTx8y4lP3kfn6SC7r7G8O6PVO539jd3+juTye5OckZc5S7k/zAsNkDk3x+mJ6X8e7u3neW7V7DVyc5M8lbhuUXJTlvmD53mM+w/qyqqnWKm2T5zN39zmFdJ3l/Jn9Db1/muR3rqjoqye9n8pycNvOxTvb7GHl+kt/p7u8M200/J+d2vLP/5+TMxztJqurkJOckefUwX5nj5ySHtenH1+LH3Uyf5939niRfXLT4UPP+RJIruvuL3f2lJFckOXtG2Zez3HHOGUlu7u5Pdfc3k7xx2HZU+zne2Shjv7/j46XMzfiv4THbTI6d95N/OXP12FmjffPMesu6vcdzOL37qCRXV9W5SW7t7g8t2uykJJ+bmr8l+38ijjXTuD8AAAWzSURBVG46d5KHJvm3w+nqv62qHx02m/fcL07y+1X1uSR/kORlw2Zzk3u4bOC6JHsyefJ+Msmd3X3XEtnuzj2s/3KSB61v4u/N3N1XT627V5LnJPnrYdHcjvWQ+5eSXNbdty3afC7GOlk29z9P8jM1uYT8r6rq1MW5B/M23r+Q5J1VdUsmj5MLhs3nZryT/FEmL0R8Z5h/UOb8OclhoZO8q6quqartw7ITpv5vuj3JCcP03DzPFznUvPP2e/zScEnhhTVcqpo5zr7oeGfDjf2i/MkGGP81Omab2djv5/jtd4exf2VV3Xtx/kU5Z5V/LfbNMxv7dSmeVXX/JG/NpADdleQ3kvzWetz3akzn7u6vZPLnZ47L5FT7f0jypnl8VX+J3M9P8pLufnCSlyR5zSzzLaW7v93dp2dyhvCMJA+fcaQDWpy5qh4xtfq/JnlPd/+/s0m3vCVy/1iSn0ryn2ebbP+WGe97J/mn7l5I8v8kuXCWGZeyTO6XJHlad5+c5LVJ/nCWGRerqqcn2dPd18w6C0ecx3f3o5M8NckLhv+f7jZcTbJhPo5/o+VN8ieZvKB3epLbkvyn2cbZvyWOd+62EcZ+ifwbYvw34jHbtGX2yy/L5Pf40UyO9X99hhGXdDjsm0cvnsOZn7cmeUN3X5rJE+qUJB+qqt2Z/KNfW1U/lOTWTN77uc/Jw7J1t0TuZPKKwKXD6fb3Z/Jqw/GZ/9zbkuybfnO+ezp9bnLv0913JrkqyWMzuZxh39+anc52d+5h/QOT/P06R73bVOazh0wvT7I5k2vw95nnsX5ikockuXl4Tt63qm4eNpursU6+Z7xvyXcf229L8iPD9DyP91OTPHLqFdZLkvybYXpexvtxSX5yeDy8MZPLeF6VDfKcZOPq7luH73syeU6fkeSOfZfQDt/3XVI/d8/zwaHmnZvfo7vvGA7Kv5PJi3kHOl6YWfZljnc2zNgvlX8jjX+y6mO2mT/up48nhsufuydvp3tt5nPs12rfPLOxH/tTbSuTs2s3dvcfJkl3f6S7f7C7t3b31kwOHB/d3bcnuSzJz9bEY5J8eYlL/0a3VO7BX2ZykJ6qemgmb+T+Qia5z6/Jp0edkuTUTN7ft672k/vzSX58mD4zyU3D9LyM9+aqOmaY/v4kT87k/Q5XJXnGsNm2JG8fpi8b5jOs/5vhlc11s0zmj1XVL2Ry3f8zhx3HPvM81td09w9NPSf/obsfMpV7pmO9n9wfy9RzMpPH+CeG6Xke7xuTPHD4PyRTy/blnvl4d/fLuvvk4fFw/pDjWZnj5yQbX1Xdr6oesG86yVOSXJ97Pr4WP+5m/jxfwqHm/e9JnlJVxw6XVj5lWLbu6p7vkf13mYx/svxxzgeSnFqTT9U8OpP/Ly5bh5zLHe9siLFfLv9GGP81PGabybHzfo7f9r1gUZm8R3J67OfisbOG++bZ9ZYe91OXHp/JZQ4fTnLd8PW0Rdvsznc/1baS/JdMrhX/SJKFMfMdau5MiuafZ/JgvDbJmVO3+c0h98czfILsHOV+fJJrMvkEq6uT/Ks5G+8fSfLBIff1SX5rWP7DmTwRbs7kTO2+TxO+zzB/87D+h+co813DeO4b/33L53qsF20z/am2Mx/rA4z3MUkuH8b0vZmcSZz78c7kgOIjw3Py3fvGdV7Ge9Hv8IR895Pz5vY56Wvjfw2Prw8NXzck+c1h+YOSXJnJi6b/I8lxw/KZP88z+fTy25J8K5MX0p+3krxJfn54/tyc5LkzzP76IduHMzk4PXFq+yWPczI5zvjEsO431yn7csc7G2Xsl8s/9+OfNTxmW+53mlH+vxnG/vpMjvX3ffLtXD12pu77CVnFvnkWY9/dqeHOAQAAYBTr9qm2AAAAHJkUTwAAAEaleAIAADAqxRMAAIBRKZ4AAACMSvEEAABgVIonAAAAo1I8AQAAGNX/D8QOULFdc4psAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(16,6))\n",
    "plt.subplot(1,2,1)\n",
    "plt.hist(nll_bayes,30)\n",
    "plt.title('Bayes')\n",
    "plt.subplot(1,2,2)\n",
    "plt.hist(nll_max_like,30)\n",
    "plt.title('MaxLike')\n",
    "print('Means ',np.mean(nll_bayes),np.mean(nll_max_like))\n",
    "np.median(nll_bayes), np.median(nll_max_like)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}