1_Bayesian_Regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm, uniform\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from scipy.optimize import minimize\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAib0lEQVR4nO3de3TT9f3H8VdaJEFNvppKSXooWnETS0UtAlYcRyaXsq1HFDd19gyYh82egoOeTalHrZ2Xop6Dbmda58687GDH5tnQlXOsx+EocxYKdExrJwLrZh1py+hICp5GbPL7g18zctpimzb9NM3zcU7OWb75Nnlrz06ffi+f2MLhcFgAAAAGpJgeAAAAJC9CBAAAGEOIAAAAYwgRAABgDCECAACMIUQAAIAxhAgAADCGEAEAAMaMMz3AmYRCIR0+fFhOp1M2m830OAAAYADC4bA6OzuVkZGhlJQzH/MY1SFy+PBhZWZmmh4DAADEoKWlRZMnTz7jPqM6RJxOp6RT/yAul8vwNAAAYCACgYAyMzMjf8fPZFSHSM/pGJfLRYgAAJBgBnJZBRerAgAAYwgRAABgDCECAACMIUQAAIAxhAgAADCGEAEAAMYQIgAAwBhCBAAAGDOqFzQDAADx0R0Kq765Q+2dXUp3OjQ7y63UlJH/XjdCBACAJFPT6FN5dZN8/q7INq/lUFlBtvJzvCM6C6dmAABIIjWNPhVtaoiKEElq9XepaFODahp9IzoPIQIAQJLoDoVVXt2kcB+v9Wwrr25Sd6ivPeKDEAEAIEnUN3f0OhJyurAkn79L9c0dIzYTIQIAQJJo7+w/QmLZbzgQIgAAJIl0p2NY9xsOhAgAAElidpZbXsuh/m7StenU3TOzs9wjNhMhAgBAkkhNsamsIFuSesVIz/OyguwRXU+EEAEAIInk53hVWZgrjxV9+sVjOVRZmDvi64iwoBkAAEkmP8erhdkeVlYFAABmpKbYlDc1zfQYsZ+aqays1IwZM+RyueRyuZSXl6c33ngj8npXV5eKi4uVlpamc889V8uWLVNbW9uwDA0AAMaGmENk8uTJ2rBhg/bu3as9e/boq1/9qm688UZ98MEHkqR169apurpar776qmpra3X48GHdfPPNwzY4AABIfLZwODxs67i63W49+eSTuuWWWzRx4kRVVVXplltukSR9+OGHuuyyy1RXV6drrrlmQO8XCARkWZb8fr9cLtdwjQkAAOJoMH+/h+Wume7ubm3evFknTpxQXl6e9u7dq5MnT2rBggWRfaZNm6YpU6aorq6u3/cJBoMKBAJRDwAAMHYNKUTef/99nXvuubLb7brrrru0ZcsWZWdnq7W1VePHj9d5550Xtf+kSZPU2tra7/tVVFTIsqzIIzMzcyjjAQCAUW5IIXLppZdq37592rVrl4qKirR8+XI1NTXF/H6lpaXy+/2RR0tLy1DGAwAAo9yQbt8dP368LrnkEknSzJkztXv3bv3kJz/Rrbfeqs8++0zHjh2LOirS1tYmj8fT7/vZ7XbZ7fahjAQAABLIsK6sGgqFFAwGNXPmTJ111lnatm1b5LX9+/fr448/Vl5e3nB+JAAASGAxHxEpLS3VkiVLNGXKFHV2dqqqqkrbt2/Xm2++KcuydOedd6qkpERut1sul0tr1qxRXl7egO+YAQAAY1/MIdLe3q7vfOc78vl8sixLM2bM0JtvvqmFCxdKkp566imlpKRo2bJlCgaDWrx4sZ599tlhGxwAACS+YV1HZLixjggAAIlnxNcRAQAAiAUhAgAAjCFEAACAMYQIAAAwhhABAADGECIAAMAYQgQAABhDiAAAAGMIEQAAYAwhAgAAjCFEAACAMYQIAAAwhhABAADGECIAAMAYQgQAABhDiAAAAGMIEQAAYAwhAgAAjCFEAACAMYQIAAAwhhABAADGECIAAMAYQgQAABhDiAAAAGMIEQAAYAwhAgAAjCFEAACAMYQIAAAwhhABAADGECIAAMAYQgQAABhDiAAAAGMIEQAAYEzMIVJRUaFZs2bJ6XQqPT1dS5cu1f79+6P2uf7662Wz2aIed91115CHBgAAY0PMIVJbW6vi4mLt3LlTb731lk6ePKlFixbpxIkTUfutWrVKPp8v8njiiSeGPDQAABgbxsX6gzU1NVHPX3rpJaWnp2vv3r2aN29eZPvZZ58tj8cT+4QAAGDMGrZrRPx+vyTJ7XZHbX/llVd0wQUXKCcnR6Wlpfr000/7fY9gMKhAIBD1AAAAY1fMR0ROFwqFtHbtWs2dO1c5OTmR7d/+9rd14YUXKiMjQ++9957uvfde7d+/X7///e/7fJ+KigqVl5cPx0gAACAB2MLhcHiob1JUVKQ33nhD77zzjiZPntzvfm+//bZuuOEGHTx4UFOnTu31ejAYVDAYjDwPBALKzMyU3++Xy+Ua6pgAAGAEBAIBWZY1oL/fQz4isnr1am3dulU7duw4Y4RI0pw5cySp3xCx2+2y2+1DHQkAACSImEMkHA5rzZo12rJli7Zv366srKwv/Jl9+/ZJkrxeb6wfCwAAxpCYQ6S4uFhVVVV6/fXX5XQ61draKkmyLEsTJkzQoUOHVFVVpa997WtKS0vTe++9p3Xr1mnevHmaMWPGsP0DAACAxBXzNSI2m63P7S+++KJWrFihlpYWFRYWqrGxUSdOnFBmZqZuuukm3X///QO+3mMw55gAAMDoMCLXiHxRv2RmZqq2tjbWtwcAAEmA75oBAADGECIAAMAYQgQAABhDiAAAAGMIEQAAYMywfNcMAACQukNh1Td3qL2zS+lOh2ZnuZWa0vdyFziFEAEAYBjUNPpUXt0kn78rss1rOVRWkK38HFYU7w+nZgAAGKKaRp+KNjVERYgktfq7VLSpQTWNPkOTjX6ECAAgojsUVt2ho3p9379Vd+ioukND/oL2Ma87FFZ5dZP6+jfVs628uol/l/3g1AwAQBKnFmJV39zR60jI6cKSfP4u1Td3KG9q2sgNliA4IgIA4NTCELR39h8hseyXbAgRAEhynFoYmnSnY1j3SzaECAAkucGcWkBvs7Pc8loO9XeTrk2nTnHNznKP5FgJgxABgCTHqYWhSU2xqawgW5J6xUjP87KCbNYT6QchAgBJjlMLQ5ef41VlYa48VvS/I4/lUGVhLhf7ngF3zQBAkus5tdDq7+rzOhGbTv1B5dTCmeXneLUw28PKqoNEiABAkus5tVC0qUE2KSpGOLUwOKkpNm7RHSROzQAAOLUAYzgiAgCQxKkFmEGIAAAiOLWAkcapGQAAYAwhAgAAjCFEAACAMYQIAAAwhhABAADGECIAAMAYQgQAABhDiAAAAGMIEQAAYAwhAgAAjCFEAACAMYQIAAAwhhABAADGxBwiFRUVmjVrlpxOp9LT07V06VLt378/ap+uri4VFxcrLS1N5557rpYtW6a2trYhDw0AAMaGmEOktrZWxcXF2rlzp9566y2dPHlSixYt0okTJyL7rFu3TtXV1Xr11VdVW1urw4cP6+abbx6WwQEAQOKzhcPh8HC80ZEjR5Senq7a2lrNmzdPfr9fEydOVFVVlW655RZJ0ocffqjLLrtMdXV1uuaaa77wPQOBgCzLkt/vl8vlGo4xAQBAnA3m7/ewXSPi9/slSW63W5K0d+9enTx5UgsWLIjsM23aNE2ZMkV1dXXD9bEAAIO6Q2HVHTqq1/f9W3WHjqo7NCz/bYskMm443iQUCmnt2rWaO3eucnJyJEmtra0aP368zjvvvKh9J02apNbW1j7fJxgMKhgMRp4HAoHhGA8AEAc1jT6VVzfJ5++KbPNaDpUVZCs/x2twMiSSYTkiUlxcrMbGRm3evHlI71NRUSHLsiKPzMzM4RgPADDMahp9KtrUEBUhktTq71LRpgbVNPoMTYZEM+QQWb16tbZu3ao//elPmjx5cmS7x+PRZ599pmPHjkXt39bWJo/H0+d7lZaWyu/3Rx4tLS1DHQ8AMMy6Q2GVVzepr5MwPdvKq5s4TYMBiTlEwuGwVq9erS1btujtt99WVlZW1OszZ87UWWedpW3btkW27d+/Xx9//LHy8vL6fE+73S6XyxX1AACMLvXNHb2OhJwuLMnn71J9c8fIDYWEFfM1IsXFxaqqqtLrr78up9MZue7DsixNmDBBlmXpzjvvVElJidxut1wul9asWaO8vLwB3TEDABid2jv7j5BY9kNyizlEKisrJUnXX3991PYXX3xRK1askCQ99dRTSklJ0bJlyxQMBrV48WI9++yzMQ8LADAv3ekY1v2Q3IZtHZF4YB0RABh9ukNhXff422r1d/V5nYhNksdy6J17v6rUFNtIj4dRwMg6IgCA5JCaYlNZQbakU9Fxup7nZQXZRAgGhBABAAxafo5XlYW58ljRp188lkOVhbmsI4IBG5YFzQAAySc/x6uF2R7VN3eovbNL6U6HZme5ORKCQSFEAAAxS02xKW9qmukxkMA4NQMAAIwhRAAAgDGECAAAMIYQAQAAxhAiAADAGEIEAAAYQ4gAAABjCBEAAGAMIQIAAIwhRAAAgDGECAAAMIYQAQAAxhAiAADAGEIEAAAYQ4gAAABjCBEAAGAMIQIAAIwhRAAAgDGECAAAMIYQAQAAxhAiAADAGEIEAAAYQ4gAAABjxpkeAAASTXcorPrmDrV3dind6dDsLLdSU2ymxwISEiECAINQ0+hTeXWTfP6uyDav5VBZQbbyc7wGJwMSE6dmAGCAahp9KtrUEBUhktTq71LRpgbVNPoMTQYkLkIEAAagOxRWeXWTwn281rOtvLpJ3aG+9gDQH0IEAAagvrmj15GQ04Ul+fxdqm/uGLmhgDGAEAGAAWjv7D9CYtkPwCkxh8iOHTtUUFCgjIwM2Ww2vfbaa1Gvr1ixQjabLeqRn58/1HkBwIh0p2NY9wNwSswhcuLECV1xxRV65pln+t0nPz9fPp8v8vj1r38d68cBgFGzs9zyWg71d5OuTafunpmd5R7JsYCEF/Ptu0uWLNGSJUvOuI/dbpfH44n1IwBg1EhNsamsIFtFmxpkk6IuWu2Jk7KCbNYTAQYprteIbN++Xenp6br00ktVVFSko0ePxvPjACCu8nO8qizMlceKPv3isRyqLMxlHREgBnFb0Cw/P18333yzsrKydOjQId13331asmSJ6urqlJqa2ufPBINBBYPByPNAIBCv8QAgJvk5Xi3M9rCyKjBM4hYit912W+R/X3755ZoxY4amTp2q7du364YbbujzZyoqKlReXh6vkQBgWKSm2JQ3Nc30GMCYMGK371588cW64IILdPDgwX73KS0tld/vjzxaWlpGajwAAGDAiH3XzCeffKKjR4/K6+3/HKrdbpfdbh+pkQAAgGExh8jx48ejjm40Nzdr3759crvdcrvdKi8v17Jly+TxeHTo0CHdc889uuSSS7R48eJhGRwAACS+mENkz549mj9/fuR5SUmJJGn58uWqrKzUe++9p5dfflnHjh1TRkaGFi1apIcffpgjHgAAIMIWDodH7Tc0BQIBWZYlv98vl8tlehwAADAAg/n7zXfNAAAAYwgRAABgDCECAACMIUQAAIAxhAgAADCGEAEAAMaM2MqqADCWdIfCfPEdMAwIEQAYpJpGn8qrm+Tzd0W2eS2HygqylZ/T/9dYAOiNUzMAMAg1jT4VbWqIihBJavV3qWhTg2oafYYmAxITIQIAA9QdCqu8ukl9LUfds628ukndoVG7YDUw6hAiADBA9c0dvY6EnC4syefvUn1zx8gNBSQ4QgQABqi9s/8IiWU/AIQIAAxYutMxrPsBIEQAYMBmZ7nltRzq7yZdm07dPTM7yz2SYwEJjRABgAFKTbGprCBbknrFSM/zsoJs1hMBBoEQAYBByM/xqrIwVx4r+vSLx3KosjCXdUSAQWJBMwAYpPwcrxZme1hZFRgGhAgAxCA1xaa8qWmmxwASHqdmAACAMYQIAAAwhhABAADGECIAAMAYQgQAABhDiAAAAGMIEQAAYAwhAgAAjCFEAACAMYQIAAAwhhABAADGECIAAMAYQgQAABhDiAAAAGMIEQAAYEzMIbJjxw4VFBQoIyNDNptNr732WtTr4XBYDz74oLxeryZMmKAFCxbowIEDQ50XAACMITGHyIkTJ3TFFVfomWee6fP1J554Qj/96U/13HPPadeuXTrnnHO0ePFidXV1xTwsAAAYW8bF+oNLlizRkiVL+nwtHA7r6aef1v33368bb7xRkvSrX/1KkyZN0muvvabbbrst1o8FAABjSFyuEWlublZra6sWLFgQ2WZZlubMmaO6urp4fCQAAEhAMR8ROZPW1lZJ0qRJk6K2T5o0KfJaX4LBoILBYOR5IBCIx3gAAGCUGFV3zVRUVMiyrMgjMzPT9EgAACCO4hIiHo9HktTW1ha1va2tLfJaX0pLS+X3+yOPlpaWeIwHAABGibiESFZWljwej7Zt2xbZFggEtGvXLuXl5fX7c3a7XS6XK+oBAADGrpivETl+/LgOHjwYed7c3Kx9+/bJ7XZrypQpWrt2rR555BF96UtfUlZWlh544AFlZGRo6dKlwzE3AAAYA2IOkT179mj+/PmR5yUlJZKk5cuX66WXXtI999yjEydO6Hvf+56OHTum6667TjU1NXI4HEOfGkhQ3aGw6ps71N7ZpXSnQ7Oz3EpNsZkeCwCMsYXD4bDpIfoTCARkWZb8fj+naZDwahp9Kq9uks//v0X9vJZDZQXZys/xGpwMAIbXYP5+j6q7ZoCxqqbRp6JNDVERIkmt/i4VbWpQTaPP0GQAYBYhAsRZdyis8uom9XXosWdbeXWTukOj9uAkAMQNIQLEWX1zR68jIacLS/L5u1Tf3DFyQwHAKEGIAHHW3jmwL3oc6H4AMJYQIkCcpTsHdqfYQPcDgLGEEAHibHaWW17Lof5u0rXp1N0zs7PcIzkWAIwKhAgQZ6kpNpUVZEtSrxjpeV5WkM16IgCSEiECjID8HK8qC3PlsaJPv3gshyoLc1lHBEDSinllVQCDk5/j1cJsDyurAsBpCBFgBKWm2JQ3Nc30GAAwanBqBgAAGEOIAAAAYwgRAABgDCECAACMIUQAAIAxhAgAADCGEAEAAMYQIgAAwBhCBAAAGEOIAAAAYwgRAABgDCECAACMIUQAAIAxhAgAADCGEAEAAMaMMz0AYtMdCqu+uUPtnV1Kdzo0O8ut1BSb6bEAABgUQiQB1TT6VF7dJJ+/K7LNazlUVpCt/ByvwckAABgcTs0kmJpGn4o2NURFiCS1+rtUtKlBNY0+Q5MBADB4hEgC6Q6FVV7dpHAfr/VsK69uUneorz0AABh9CJEEUt/c0etIyOnCknz+LtU3d4zcUAAADAEhkkDaO/uPkFj2AwDANEIkgaQ7HcO6HwAAphEiCWR2lltey6H+btK16dTdM7Oz3CM5FgAAMSNEEkhqik1lBdmS1CtGep6XFWSznggAIGHENUQeeugh2Wy2qMe0adPi+ZFjXn6OV5WFufJY0adfPJZDlYW5rCMCAEgocV/QbPr06frjH//4vw8cxxpqQ5Wf49XCbA8rqwIAEl7cq2DcuHHyeDzx/pikk5piU97UNNNjAAAwJHG/RuTAgQPKyMjQxRdfrDvuuEMff/xxv/sGg0EFAoGoBwAAGLviGiJz5szRSy+9pJqaGlVWVqq5uVlf+cpX1NnZ2ef+FRUVsiwr8sjMzIzneAAAwDBbOBwesfXAjx07pgsvvFAbN27UnXfe2ev1YDCoYDAYeR4IBJSZmSm/3y+XyzVSYwIAgCEIBAKyLGtAf79H9MrR8847T1/+8pd18ODBPl+32+2y2+0jORIAADBoRNcROX78uA4dOiSvl1tMAQBAnEPkhz/8oWpra/XPf/5T7777rm666Salpqbq9ttvj+fHAgCABBHXUzOffPKJbr/9dh09elQTJ07Uddddp507d2rixInx/FgAAJAg4hoimzdvjufbAwCABMd3zQAAAGMIEQAAYAwhAgAAjCFEAACAMYQIAAAwhhABAADGECIAAMAYQgQAABhDiAAAAGMIEQAAYAwhAgAAjCFEAACAMYQIAAAwhhABAADGECIAAMAYQgQAABhDiAAAAGMIEQAAYAwhAgAAjCFEAACAMYQIAAAwhhABAADGECIAAMAYQgQAABhDiAAAAGMIEQAAYAwhAgAAjCFEAACAMYQIAAAwhhABAADGECIAAMAYQgQAABhDiAAAAGPiHiLPPPOMLrroIjkcDs2ZM0f19fXx/kgAAJAg4hoiv/nNb1RSUqKysjI1NDToiiuu0OLFi9Xe3h7PjwUAAAkiriGyceNGrVq1SitXrlR2draee+45nX322XrhhRfi+bEAACBBxC1EPvvsM+3du1cLFiz434elpGjBggWqq6vr82eCwaACgUDUAwAAjF1xC5H//Oc/6u7u1qRJk6K2T5o0Sa2trX3+TEVFhSzLijwyMzPjNR4AABgFRtVdM6WlpfL7/ZFHS0uL6ZEAAEAcjYvXG19wwQVKTU1VW1tb1Pa2tjZ5PJ4+f8Zut8tut8drJAAAMMrE7YjI+PHjNXPmTG3bti2yLRQKadu2bcrLy4vXxwIAgAQStyMiklRSUqLly5fr6quv1uzZs/X000/rxIkTWrlyZTw/FgAAJIi4hsitt96qI0eO6MEHH1Rra6uuvPJK1dTU9LqAFQAAJCdbOBwOmx6iP4FAQJZlye/3y+VymR4HAAAMwGD+fo+qu2YAAEByIUQAAIAxhAgAADCGEAEAAMYQIgAAwBhCBAAAGEOIAAAAYwgRAABgDCECAACMIUQAAIAxhAgAADCGEAEAAMYQIgAAwBhCBAAAGEOIAAAAYwgRAABgDCECAACMIUQAAIAxhAgAADCGEAEAAMYQIgAAwBhCBAAAGEOIAAAAYwgRAABgDCECAACMIUQAAIAxhAgAADCGEAEAAMYQIgAAwBhCBAAAGEOIAAAAY8aZHgCx6Q6FVd/cofbOLqU7HZqd5VZqis30WAAADAohkoBqGn0qr26Sz98V2ea1HCoryFZ+jtfgZAAADE7cTs1cdNFFstlsUY8NGzbE6+OSRk2jT0WbGqIiRJJa/V0q2tSgmkafockAABi8uB4R+fGPf6xVq1ZFnjudznh+3JjXHQqrvLpJ4T5eC0uySSqvbtLCbA+naQAACSGuIeJ0OuXxeOL5EUmlvrmj15GQ04Ul+fxdqm/uUN7UtJEbDACAGMX1rpkNGzYoLS1NV111lZ588kl9/vnnZ9w/GAwqEAhEPfA/7Z39R0gs+wEAYFrcjojcfffdys3Nldvt1rvvvqvS0lL5fD5t3Lix35+pqKhQeXl5vEZKeOlOx7DuBwCAabZwONzXJQd9Wr9+vR5//PEz7vP3v/9d06ZN67X9hRde0Pe//30dP35cdru9z58NBoMKBoOR54FAQJmZmfL7/XK5XAMdc8zqDoV13eNvq9Xf1ed1IjZJHsuhd+79KteIAACMCQQCsixrQH+/BxUiR44c0dGjR8+4z8UXX6zx48f32v7BBx8oJydHH374oS699NIBfd5g/kGSRc9dM5KiYqQnOyoLc7mFFwBg1GD+fg/q1MzEiRM1ceLEmIbat2+fUlJSlJ6eHtPP45T8HK8qC3N7rSPiYR0RAEACiss1InV1ddq1a5fmz58vp9Opuro6rVu3ToWFhTr//PPj8ZFJJT/Hq4XZHlZWBQAkvLiEiN1u1+bNm/XQQw8pGAwqKytL69atU0lJSTw+Limlpti4RRcAkPDiEiK5ubnauXNnPN4aAACMIXz7LgAAMIYQAQAAxhAiAADAGEIEAAAYE9cvvRutukNhbn0FAGAUSLoQqWn09VoMzMtiYAAAGJFUp2Z6lkc/PUIkqdXfpaJNDapp9BmaDACA5JQ0IdIdCqu8uqnPL4vr2VZe3aTu0IC/egcAAAxR0oRIfXNHryMhpwtL8vm7VN/cMXJDAQCQ5JImRNo7+4+QWPYDAABDlzQhku50DOt+AABg6JImRGZnueW1HOrvJl2bTt09MzvLPZJjAQCQ1JImRFJTbCoryJakXjHS87ysIJv1RAAAGEFJEyKSlJ/jVWVhrjxW9OkXj+VQZWEu64gAADDCkm5Bs/wcrxZme1hZFQCAUSDpQkQ6dZomb2qa6TEAAEh6SXVqBgAAjC6ECAAAMIYQAQAAxhAiAADAGEIEAAAYQ4gAAABjCBEAAGAMIQIAAIwhRAAAgDGjemXVcDgsSQoEAoYnAQAAA9Xzd7vn7/iZjOoQ6ezslCRlZmYangQAAAxWZ2enLMs64z628EByxZBQKKTDhw/L6XTKZuNL6UZKIBBQZmamWlpa5HK5TI8D8TsZjfidjD78TkaPcDiszs5OZWRkKCXlzFeBjOojIikpKZo8ebLpMZKWy+Xi/8yjDL+T0YffyejD72R0+KIjIT24WBUAABhDiAAAAGMIEfRit9tVVlYmu91uehT8P34now+/k9GH30liGtUXqwIAgLGNIyIAAMAYQgQAABhDiAAAAGMIEQAAYAwhgoiKigrNmjVLTqdT6enpWrp0qfbv3296LJxmw4YNstlsWrt2relRktq///1vFRYWKi0tTRMmTNDll1+uPXv2mB4raXV3d+uBBx5QVlaWJkyYoKlTp+rhhx8e0PecwLxRvbIqRlZtba2Ki4s1a9Ysff7557rvvvu0aNEiNTU16ZxzzjE9XtLbvXu3fv7zn2vGjBmmR0lq//3vfzV37lzNnz9fb7zxhiZOnKgDBw7o/PPPNz1a0nr88cdVWVmpl19+WdOnT9eePXu0cuVKWZalu+++2/R4+ALcvot+HTlyROnp6aqtrdW8efNMj5PUjh8/rtzcXD377LN65JFHdOWVV+rpp582PVZSWr9+vf7yl7/oz3/+s+lR8P++8Y1vaNKkSfrlL38Z2bZs2TJNmDBBmzZtMjgZBoJTM+iX3++XJLndbsOToLi4WF//+te1YMEC06MkvT/84Q+6+uqr9c1vflPp6em66qqr9Itf/ML0WEnt2muv1bZt2/TRRx9Jkv72t7/pnXfe0ZIlSwxPhoHg1Az6FAqFtHbtWs2dO1c5OTmmx0lqmzdvVkNDg3bv3m16FEj6xz/+ocrKSpWUlOi+++7T7t27dffdd2v8+PFavny56fGS0vr16xUIBDRt2jSlpqaqu7tbjz76qO644w7To2EACBH0qbi4WI2NjXrnnXdMj5LUWlpa9IMf/EBvvfWWHA6H6XGgU5F+9dVX67HHHpMkXXXVVWpsbNRzzz1HiBjy29/+Vq+88oqqqqo0ffp07du3T2vXrlVGRga/kwRAiKCX1atXa+vWrdqxY4cmT55sepyktnfvXrW3tys3Nzeyrbu7Wzt27NDPfvYzBYNBpaamGpww+Xi9XmVnZ0dtu+yyy/S73/3O0ET40Y9+pPXr1+u2226TJF1++eX617/+pYqKCkIkARAiiAiHw1qzZo22bNmi7du3Kysry/RISe+GG27Q+++/H7Vt5cqVmjZtmu69914ixIC5c+f2uq39o48+0oUXXmhoInz66adKSYm+5DE1NVWhUMjQRBgMQgQRxcXFqqqq0uuvvy6n06nW1lZJkmVZmjBhguHpkpPT6ex1jc4555yjtLQ0rt0xZN26dbr22mv12GOP6Vvf+pbq6+v1/PPP6/nnnzc9WtIqKCjQo48+qilTpmj69On661//qo0bN+q73/2u6dEwANy+iwibzdbn9hdffFErVqwY2WHQr+uvv57bdw3bunWrSktLdeDAAWVlZamkpESrVq0yPVbS6uzs1AMPPKAtW7aovb1dGRkZuv322/Xggw9q/PjxpsfDFyBEAACAMawjAgAAjCFEAACAMYQIAAAwhhABAADGECIAAMAYQgQAABhDiAAAAGMIEQAAYAwhAgAAjCFEAACAMYQIAAAwhhABAADG/B+Kr63E53GyIAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Bayesian Regression\n",
    "np.random.seed(42)\n",
    "x = np.random.uniform(0, 10, 10)\n",
    "y = 3 * x + np.random.normal(0, 5, 10)\n",
    "\n",
    "# Plot the data points\n",
    "plt.scatter(x, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficients:\n",
      "[2.79015326]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "lm = LinearRegression(fit_intercept=False)\n",
    "lm.fit(x.reshape(-1, 1), y)\n",
    "print(\"Coefficients:\")\n",
    "print(lm.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Maximum Likelihood Estimate\n",
    "def likelihood(beta, sd):\n",
    "    y_pred = beta * x\n",
    "    res = y - y_pred\n",
    "    l = norm.logpdf(res, loc=0, scale=sd).sum()\n",
    "    return l\n",
    "\n",
    "\n",
    "def ml_likelihood(params):\n",
    "    beta = params[0]\n",
    "    sd = params[1]\n",
    "    if sd <= 0:\n",
    "        return -np.inf\n",
    "    y_pred = beta * x\n",
    "    res = y - y_pred\n",
    "    l = norm.logpdf(res, loc=0, scale=sd).sum()\n",
    "    return -l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum Likelihood Estimates:\n",
      "[2.79016478 4.59534662]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "result = minimize(ml_likelihood, [1, 1], method=\"Nelder-Mead\")\n",
    "print(\"Maximum Likelihood Estimates:\")\n",
    "print(result.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bayesian Regression\n",
    "def posterior_estimate(beta, sd, prior_beta, prior_sd):\n",
    "    l = likelihood(beta, sd) + prior_beta, prior_sd\n",
    "    return l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Brute force result\n",
    "\n",
    "# Set a strong prior of beta = 3 with a standard deviation = 0.1\n",
    "betas = np.random.uniform(-10, 15, 10000)\n",
    "beta_p = norm.logpdf(betas, loc=5, scale=3)\n",
    "\n",
    "# SDs prior\n",
    "sds = np.random.uniform(10, 30, 10000)\n",
    "sds_p = uniform.logpdf(sds, loc=10, scale=20)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up parameters\n",
    "iterations = 20000\n",
    "burnin = 2000\n",
    "\n",
    "# Define support functions\n",
    "def jump_function(x_in):\n",
    "    b_proposal = np.random.normal(x_in[0], 1)\n",
    "    sd_proposal = np.random.uniform(x_in[1] - 0.5, x_in[1] + 0.5)\n",
    "    if b_proposal < 0:\n",
    "        b_proposal = 0.01\n",
    "    if sd_proposal < 15:\n",
    "        sd_proposal = 15\n",
    "    return {\"b_proposal\": b_proposal, \"sd_proposal\": sd_proposal}\n",
    "\n",
    "def prior_likelihood(params):\n",
    "    return np.log(norm.pdf(params[0], 5, 3)) + np.log(norm.pdf(params[1], 10, 30))\n",
    "\n",
    "def posterior_mh(beta, sd):\n",
    "    return likelihood(beta, sd) + prior_likelihood([beta, sd])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the parameter matrix\n",
    "parameters = np.empty((iterations, 5))\n",
    "parameters[:, :] = np.nan\n",
    "parameters[0, :] = [betas[0], sds[0], beta_p[0], sds_p[0], np.nan]\n",
    "\n",
    "# Run the Metropolis-Hastings algorithm\n",
    "for i in range(iterations-1):\n",
    "    # Calculate the posterior for the current iteration\n",
    "    parameters[i, 4] = posterior_mh(parameters[i, 0], parameters[i, 1])\n",
    "  \n",
    "    # Find a proposal point\n",
    "    proposals = jump_function(parameters[i, :2])\n",
    "  \n",
    "    # Calculate the posterior for the proposal\n",
    "    proposal_likelihood = posterior_mh(proposals['b_proposal'], proposals['sd_proposal'])\n",
    "  \n",
    "    acceptance_ratio = np.exp(proposal_likelihood - parameters[i, 4])\n",
    "  \n",
    "    # If the proposal likelihood is greater, always accept the proposal\n",
    "    if acceptance_ratio > 1:\n",
    "        parameters[i+1, :2] = proposals['b_proposal'], proposals['sd_proposal']\n",
    "    else:\n",
    "        # Otherwise, sometimes accept and sometimes explore\n",
    "        p = np.random.uniform(0, 1)\n",
    "        if p > (1 - acceptance_ratio):\n",
    "            parameters[i+1, :2] = proposals['b_proposal'], proposals['sd_proposal']\n",
    "        else:\n",
    "            parameters[i+1, :2] = parameters[i, :2]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb1004ec370>]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAGdCAYAAABO2DpVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB6F0lEQVR4nO2dd5wURfbAX8/mhQ3AkhaWnCQnQUQxgATRw/NOPY9TzKeHZ/YUc14MZzjPQ89T8M6A+jtRTxEkK5IRkByXnMNG2Fy/P2Z7trunuruquzpMz/vy2Q+7M91Vr7urq169evWeRAghgCAIgiAIIoCQ1wIgCIIgCBIcULFAEARBEEQYqFggCIIgCCIMVCwQBEEQBBEGKhYIgiAIgggDFQsEQRAEQYSBigWCIAiCIMJAxQJBEARBEGEkul1hbW0tHDx4EDIyMkCSJLerRxAEQRDEAoQQKCkpgdzcXAiF9O0SrisWBw8ehLy8PLerRRAEQRBEAPv27YPWrVvrfu+6YpGRkQEAYcEyMzPdrh5BEARBEAsUFxdDXl5eZBzXw3XFQl7+yMzMRMUCQRAEQWIMMzcGdN5EEARBEEQYqFggCIIgCCIMLsWipqYGHn/8cWjfvj2kpaVBx44d4dlnnwXMvI4gCIIgCACnj8WLL74IU6ZMgQ8++AB69OgBq1atghtvvBGysrLgrrvuckpGBEEQBEFiBC7FYsmSJTBu3DgYO3YsAAC0a9cOPvnkE1ixYoUjwiEIgiAIEltwLYWce+65MG/ePNi2bRsAAKxbtw4WL14MY8aM0T2noqICiouLVT8IgiAIggQTLovFww8/DMXFxdCtWzdISEiAmpoaeP7552H8+PG65+Tn58PTTz9tW1AEQRAEQfwPl8Xis88+g48++gg+/vhj+Pnnn+GDDz6AV155BT744APdcyZNmgRFRUWRn3379tkWGkEQBEEQfyIRji0deXl58PDDD8PEiRMjnz333HPw4YcfwpYtW5jKKC4uhqysLCgqKsIAWQiCIAgSI7CO31wWi9OnT0clHklISIDa2lprUiIIgiAIEii4fCwuv/xyeP7556FNmzbQo0cPWLNmDbz66qtw0003OSUfgiAIgiAxBNdSSElJCTz++OMwY8YMOHr0KOTm5sK1114LTzzxBCQnJzOVgUshCIIgCBJ7sI7fXIqFCFCxQBAEQZDYwxEfCwRBEF6ufmcp7DlR5rUYCIK4BCoWCII4yoqCk3CspMJrMRAEcQlULBAEQRAEEQYqFgiCIAiCCAMVCwRBEARBhIGKBYIgjiNJXkuAIIhboGKBIAiCIIgwULFAEARBEEQYqFggCIIgCCIMVCwQBEEQBBEGKhYIgiAIgggDFQsEQRAEQYSBigWCIC6A+00RJF5AxQJBEARBEGGgYoEgCIIgiDBQsUAQxHEw8iaCxA+oWCAI4jiEeC0BgiBugYoFgiAIgiDCQMUCQRAEQRBhoGKBIIjjoI8FgsQPqFggCIIgCCIMVCwQBEEQBBEGKhYIghiy61gpXPDyAltl4EoIgsQPqFggCGLIqdNVsOfEaVtl4G5TBIkfULFAEARBEEQYqFggCIIgCCIMVCwQBHEc9LFAkPgBFQsEQRAEQYSBigWCIIZgcCsEQXhAxQJBEMeRUDtBkLgBFQsEQRAEQYSBigWCII5DMG86gsQNXIpFu3btQJKkqJ+JEyc6JR+CIAiCIDFEIs/BK1euhJqamsjfGzZsgEsuuQSuuuoq4YIhCBIc0McCQeIHLsWiadOmqr8nT54MHTt2hAsuuECoUAiCIAiCxCZcioWSyspK+PDDD+G+++4znI1UVFRARUVF5O/i4mKrVSIIgiAI4nMsO29++eWXUFhYCDfccIPhcfn5+ZCVlRX5ycvLs1olgiAxCi6EIEj8YFmxeO+992DMmDGQm5treNykSZOgqKgo8rNv3z6rVSII4gGoFCAIwoOlpZA9e/bA3Llz4YsvvjA9NiUlBVJSUqxUgyAIgiBIjGHJYjF16lRo1qwZjB07VrQ8CIIgCILEMNyKRW1tLUydOhUmTJgAiYmWfT8RBEEQBAkg3IrF3LlzYe/evXDTTTc5IQ+CIAiCIDEMt8lh5MiRGJ4XQRAEQRAqmCsEQRAEQRBhoGKBIIjjYERvBIkfULFAEARBEEQYqFggCGKIiARiiSHsahAkXsC3HUEQx0lPTvBaBARBXAIVCwRBEARBhIGKBYIgCIIgwkDFAkEQBEEQYaBigSCIIbhTFEEQHlCxQBDEEIyziyAID6hYIAiCIAgiDFQsEAQRwtDJ870WAUEQH4CKBYIgQjhQeEb3O1xOQZD4ARULBEEMQedNBEF4QMUCQRAEQRBhoGKBIAiCIIgwULFAEARBEEQYqFggCGKIgOSmCILEEahYIAjiGLW19P0gT329EQ4a7CJBECR2QcUCQRBDiI29oq/P3Ub9fNqS3bD7RJn1ghEE8S2oWCAI4hjHSisBAIDY0U4QBIkpULFAEMQxEkMGDhqoayBIIEHFAkEQQ+w4byYYKRYIggQSVCwQBPEENFggSDBBxQJBEARBEGGgYoEgCIIgiDBQsUAQxDGM/DNwowiCBBNULBAEMUSykd8UlQcEiT9QsUAQxDX+b/V+r0VAEMRhULFAEMQQYmP/hrwUIpfwwOfr7AuEIIivQcUCQRAEQRBhoGKBIIgn0CwhVTW1HkiCIIhIuBWLAwcOwB/+8Ado0qQJpKWlQa9evWDVqlVOyIYgiA+w47wp8+4Pu0yP2XakBDo/+p3tuhAE8ZZEnoNPnToFQ4cOhYsuugi+++47aNq0KWzfvh0aNWrklHwIggSAvSdPR32m3TFSVlHtkjQIgjgJl2Lx4osvQl5eHkydOjXyWfv27YULhSBIsKil7Dt9efZWGNalqQfSIAjiJFxLIV9//TUMHDgQrrrqKmjWrBn069cP3n33XcNzKioqoLi4WPWDIEh8QYtnsf5AkfuCIAjiOFyKxa5du2DKlCnQuXNnmD17Ntxxxx1w1113wQcffKB7Tn5+PmRlZUV+8vLybAuNIIh72MluKsMSKMvtWFqnyipdrhFB4gMuxaK2thb69+8PL7zwAvTr1w9uu+02uPXWW+Htt9/WPWfSpElQVFQU+dm3b59toREEiS3sxMJwin7PzoHyqhqvxUCQwMGlWLRs2RK6d++u+uyss86CvXv36p6TkpICmZmZqh8EQeIDeUdJrf/0CgRBHIJLsRg6dChs3bpV9dm2bdugbdu2QoVCEMQ/2Mn3IVsqaM6bCIIEEy7F4t5774Vly5bBCy+8ADt27ICPP/4Y/vnPf8LEiROdkg9BEI84WHhGWFmoVyBI/MClWJx99tkwY8YM+OSTT6Bnz57w7LPPwuuvvw7jx493Sj4EQTygvKoGzp08HwDsOW/KSyEsegUqHwgSDLjiWAAAXHbZZXDZZZc5IQuCID5B+CCPWgOCxA2YKwRBEMdBtQJB4gdULBAEifCfpbvhyzUHIk6Xx0oqhJTL4rwpIl4GgiDeg4oFYsrHy/W3EyPBYub6wzB/y9HI32c/P1dIubgSgtB4bc42r0VAHAAVC8SUR2as91oExCVCobB1QbQiwBLHQq7z1e+3wutzccCJB96Yt91rERAHQMUCQeKY57/dpPo7JElClQp5eYNwFLr+QBFsOog5hRAkVkHFAkHimHd/LIj6jGgCcNvxfZD1CR5lBVdN7LP9SAnsOFrqtRhInIKKhQ2OFpdDVU2t12IgAeUfC3fAZyvdza0jUSwWIiwYTuQKKS6vgqPF5cLLDQLPz9wML83a4rUYSJyCioUNBr0wD2ZtOOy1GEhAWbDlKCwrOOFqnSFJbPht2dpxxoFkX6/N2QbXvbdCeLlBIDkhBJU46UE8AhULm8RCDoR2D3/rtQiIz1m95yQcLDwDIUmCWgJQUl4ltHw5Aqcx9e8Sy/JLDLx6npEQkjDxG+IZ3JE3EQRxEZcGh9s//Bn+MLg+mWBpeXXkd9ogX3C8DNrnNHBDNEu8s2gnnCyrhLNaZsIV/Vp5LY7roNKFeAlaLBDEp7DN8kXVpfaDyEhNMjz+olcWcpXP42MhYlDM/24LzN9yFH7ZX2S/sBgF440hXoGKBYIgIEnOznJ5yxalVDnhNBorxO+VI16DigWC+Bi3BgcJJEfr8sI0z7tNdu+J01zxNvwMhkdHvCSQisW/ftzltQgIYh8XBwdJAtXor5zpi7AeeDVg81Q77OUFUKzwLYllYkk/Ol0ZjHuO1BNIxeK5bzd7LQKCxBQSAPxt/g7HyucZ5wgALN1lf5utFYUIZ/ruU1qBikXQCKRigSAIH5LPRtSiM2zbXWliy9YR5Xc/bDsmQqyYwWeP05gYsq4gbKBi4QAYDRAJIhsO8O+w8Hp8k5WM69/HQFp+BfWK4IGKhQMMemGe1yIgiHAe+Hwd9znyoMGy5i/KL8BOObHkm4AgfgUVCwRB4EDhmcjveoNrZXU4RHSthZCOSmfQU2WVxscKHN0JAObz8TmozAUPVCwQxKeEN2pE97qr95yC2RudyVGzes8pmLv5iOExHR6ZyVyevBSivIx+z86xIBk/kiRBaXk1dH70O1fqQ6wRz7FGggoqFogwzs3HJSA3+O/P++HN+dsdKZvmNBlxBGRwmLj2n8tg/L+WRX3u1NBBm+0qP8JEXAjiPpgrBBHGwSJ0Wg0KtAF717Ey0/P0tol6FsfCk1oRHnApJHigxSKGmbvJ2GSNxDaSRB8YD5w6AxsOFLsujx08Gzxw0PI9sZAhGuEDFYsY5pZ/r/JaBMQDDhWdMT/IZ1gdOqav2Gv4vWEcC2AbtD5evhfW7iusO5lTQARBokDFAkFijJDL0Y+Mqrt52kpH6374i/WWz2VNrPb63G2wcOtRy/XEOj9sOwYnSit0v5+76QjMcdA6igaL4IGKhQu0e/hbr0WwTFCSMiHOMG+LPwbkg8rtsorPWXccuJmi3m9c//4KWGQQmXT6yn2mliM7YBcTPFCxcJCPlu+BotNsoYnjnbX7CuGSVxd5LQZCQUTH76SCSgiBcyfPj/pckgBYQm7guGaOo5lv60of/foPUFFd42BNiFugYuEgj87YANuPlngtRkxwvKQCth8t9VqMmMBuXo/Ve07CvpOnBUnDxikdBfuFmZth9/HwbhOzwWvJzuPc9bLqMzGVW8NlnL438jPacrgkEoQNiW1QseDg30t3W4o6iCBWcMo8f/f0tfAJh2nbzsBidu4/f9gFd3+6FiZ98YtpWb9/dzlTnfJAxXr/lMpHLAZr+nLNAfhpB7/SxYOTFielg622ln8v3e1YvYhzoGLBwRNfbYy7gDus/UnBcfMYBwg/yvtfW0ugsroW0pLsvba0wX5A20a2ytSDpf3sPl4G247UW6t+3M43SBrVwRrOmxYhNFaY+lMB/G/dQVtlGCmAel+9v7gAvlp7QPe8I8XlUFxOt1SVV9XAmcrwsofRLX/iq40G3yJ+BRULxDZVNbVw0SsLvRYj8Px76W64+p2lkJRg/7UlAFCi6PRzGiZHfW+X79YfUoUHf+6bTVznvzhri636JSk8gFmh11OzVX+v2n0STldW25LHKSRJsh0Lgvd0Qgj83+r9hkrgH/61HN6YS48Q+/T/NsE9n66JqjsWFTs7LNt1IpAO8qhYIIhP0c4iC89UwaGiM0IG/VW7T0Kvp74XUJI+01fug90n6n05/rW4gHqcPKvV9q9LLJr3+Zczok3xJeVqJeK3by+F5QUnLcnjNCEdJ1UnB6yC42Ww6ZB5kDY9a8fpyurIPSbqtai44nf/jA5/HwS4FIunnnoKJElS/XTr1s0p2XyJ08rl1+sOUnMteEUsv+eEECir8OcskwW9tmZ3wJBAgvIq9iU9p7di6l0OqzuT3mE8Usey82ZCyL7FwghaBNhqxoejd1RIkqCmroxY62OOFmPqAjO4LRY9evSAQ4cORX4WL17shFy+xY5zl956o5JDhWdg/f4i5jJXFJyEV2ZvtSxTkFm847hrmTTdRMj2T007dmJcsjtYW33X+M36yt9jbZgLL4XQxLa7e0iJldti6LehCF4WK86zh4vKgRACg16YB4WnK70Wx9dwKxaJiYnQokWLyE9OTo4TcvkW2gtWW0uYdov0fup72HLY2Hy4/9QZKC5nn2Wv2XsKPvDYc9qvk73TlTUxvX1Nr2MW0fVqrRDaMv0wwLKIUHymSjf2wbr9Rab3Sk6aJnIQBgA4VVYJ+d9thqPF5bByd3gJZcnO47DhAPukgZXwUkj0lfI8Q2Pnzegv7TaPkMIvRKlM+KDZ6XJO/jzYUbclXqScfr5mq3ArFtu3b4fc3Fzo0KEDjB8/HvbuNd62VlFRAcXFxaofvzLspQWmx9DawEuzt8L9n69jqsNsoDtZxqcJS5KOUC4SwPfCN0QP+PYHfVGzQmW0Sxq8Q7V2cGO5zINF5fDWgh2cNdVzvKRCE6lTDPtOnYZ3Fu2Cr9cdhJumhsOevzhrK/xn6R5BNdQTHqTr/z5dWe2bpR29Z6hUhmLJxYJo/kfocCkWgwcPhmnTpsGsWbNgypQpUFBQAOeffz6UlOgHgcrPz4esrKzIT15enm2hedhxtIR5lrCXIWgQrVM/UlwOB0w6WZlYC4PBO4it318E1Ta25JZVVMPsjYctn19VU+uqlWL/KfcCTckzRyEWC5OBh2UG/3eDAd2K8qM9hbUEHn8RO/CM1dSlCWGSaMpVDNIbDhRB9ydmOz4LlpVTo3qMfHMkqFeGMLtp8OBSLMaMGQNXXXUV9O7dG0aNGgUzZ86EwsJC+Oyzz3TPmTRpEhQVFUV+9u3bZ1toHt5asBOe4dzm5iRmL5Gf1xj1UF7S5X9fDFuPhBXNq99eyl3WxoPF8Mf/rAYAgGk/FcCnK/lyFDz59UZ4+L/mwZZEcd6L5lYuq/hl1qmHUVP+2mZchXD5+hV888tBmLZkNwAAvKez2yRciPrP/adOwxVv/aRbj6gxTi/oU2llNRxx0PlPub1W1PIOzXmTFb3+LBSqv+9u+rgs2Xkcvvh5v+Xz6/1CYq+fdhNb202zs7OhS5cusGOH/swlJSUFMjMzVT+uY7MNzFizX5gHMyEEVu85peugKbfX695bDpuZtnOJdNAisENACHL5Glbs5t+ed7SkvtOdu/ko/MAZLOlkaSWc0CwnzdpwiFsOv0DrwPzSpxlF7yw642yOnD0n2CxF2oHteGllfYp0BaKVOLlW5eBecLwMvv3lEFxlQeFmRXkdzuZnia6PB0mS6pcVXFwKmbXhMLz7o4Ei6gE+eZ2FYkuxKC0thZ07d0LLli1FySMcEf3FvZ+ug9LInmt7ZdUSgL9+vxX++eMu6vdy+T9uPw6Hi9hmNiwi7TxmnoejqobAiFd/4C5bxA6DW/69CgAA7vx4jW65LNDOuf3Dn/kFMmHXsVLPtp2Jcd7UlKldhlB8sMBiSnG7M2arJnLuXSGq38V087RZuJ6y9ft3l9laPlROLuR7LlJRopXFco+NZAhJoJis1RdW4fAypgT2FC75mmZvdC6NfBDgUiweeOABWLRoEezevRuWLFkCv/71ryEhIQGuvfZap+QTgpDOIhLzl1K+pfVk9TkbDhTBpoNqCwVL5xDetmVe//C/mmcOFdUZeb1mKl+Hk2L86aOfbTkNsqBrjXI5yiKL5cwJxGRVNT+m8HSVA7E6whWzlLpk5wmoEdRY3Q5NblaP9vvi8ip47ptNIIGkWFao//60w3Fn9Lbm8vL2op32CwkwXIrF/v374dprr4WuXbvC1VdfDU2aNIFly5ZB06ZNnZLPPsJNnPRWqa2GN5Twy7O3whvztqnKF70FzgxRtXmpVvhlmUDLlIU74dEZ67nOcdffRnxdvO1JqyDrSVRWUQ0vK2K3XNAl3P/U1hI4VGTsRG2qhAvzsRBTDo1LXl2k65Bu1GecKquEaT9ZWwaIfjbWLvBocTn8a3GBym/DzVaudHQFANh74jTXxDCiDAmU2o4FZcvhYl+mmudSLKZPnw4HDx6EiooK2L9/P0yfPh06duzolGzCEDnYsHYY3R6fxVWuHDBGKWvII+e9m6et5Dpee39rPN76IvK26Tl6bTlcAh9wbB3ceawUthy27r/y3uICmLk+7Csi393ddhK/+d0zFPQ7XK25PC0pAQAAVu89BUPy56vjIjgnniFmfgMfLrO+7XT70VI4VlJB/U5+qrTHu/VICTz1P7oju9kODt3vrPpYQP3z1bNwOjFghup8O+ZtPgKzNhyGYS8vgMLTzvoDOcno13+EFT4MNR/4XCGiTJz15nX7XRXNy5oWjoI57bNtidTM21K/pr7tCP9g6J8xy/6due8ztvgkTqB8/t/8chC2HimBmloSGbTu/MSa7whLE/aD5UcpgnLZSW8gogWpY12Wi7zfAPBInWXJjt+D3E9Q/ROAwGNfbrBcNoD+dbn17tGWMVgY+VrYh0vPsqIs7qH/+6WuDnGNccaaA1BwvAw+X7U/4nzMc8/807f5m8ArFgBiB14xjnPRrVNe+xvWpX5ZicViYWe5ZNaGwzDVxDQqpzb2CnyP1ShjdCj727X7CuHH7cdMz3/1+61woPAMbNcojNq+W5kLwqxfd8xCpShWufQRHe9CHsQl6vc6RarOUfLx8r2qY/O/28wkrl49TihprloFBb6EZsnSlPfqVJ0lYauFyY0eJ8sqo+4dTx+6l3E3khXaPfytpfP8MAnQEnjFQpSG6bRTVKjOCbNhSmJ0pSbwyPSIYp1/0bZjMGPNAfaTGfFSGdC7Fe8tLoDX5mxzvP6bpq2Ez1byxWo5U1kDX62lPwft9RCgr+9+uGwPvK6TolrJj3UZQ0+bKIwzFEtAZs1LzxRr993TS3Sl6+dkoz6j9/udRbvqymevILId07pIhujpFU4njBMJq9XDCSVK7cvGft59n611QBZxPP2/jQJLs07gFQsAUcsXctRDeln29+1LlpZCeLsReTbGig+VYV36PvM9FBwvow4Aa/cVRnI2OMnGg0Vw0MSBUEvB8TK4e/raqM/p2/zqn8g+hkixVqni6Myd2gVUpbcUobVYMA7in63aF5UO3SlUs3CNwCIG/6h7Lt8Dh/SK6O3I9uqTJLolgjbgOzWZsyJ7saCwAyJRXsfUn3Z7JoeSwCsWeqk0vlt/CIqsOO3oNKgth0sMFZjImqtOObLFQvliOe+8af52xMp6PEB422Cxw4GZzHDjXsh18CSrAwA4VlLBvLTFEw1Rf0ZprwHrKRY0Kw6A2k8i8p3ij7/83y+wVZMEUOsLoXz/bquLrWIF3mbAq2yYRvAV2A6NJDMM6S1Jum1Heb1my0ZOLfvIdYUsaBh+8rXwS/+rJPCKhR53fPQzrN1fyHy8QRiLCCw5CwhEN8rdx8vqtkFp6mRsvXa2PilroJUSi6Fr9czaMXgp1Ieidx1mz+q3by+xtTNFD6csFvzBkqLfF61k8vIKi1Pygq3mPit61M/oaTJF3y/Wd1h+xlGDrSTXV/+R2XKXHexut1RHCKX/Lh+izIPy8uwttuqlyiK8RD5isl8yIfCKxeer90N1jf0nJ5dg1ImWVerPII0az4WvLIRDReVRZbM5b5of4zRRpl4/CBVD8HTSdmwDRu+B1p9B1fGblOv2UohedSzNrqbuHsi7EyLnygneBF0KSypws+ywNOTHpHsPoD7y5uo9p7jK/mkHX/h8M4weR0ilWOgvGwHUX/O6/YUwZaG4wFRB6ab8eB2BVywAxO6Hlt+Bvs98b+n8esuHVomIjgjH2mCc1HiDpExbfQGLy/mWV+TO/cftx+Cqt5dYq9QAJ57Jom3q2TlPm6oVFIVZW6WeBdAsSB3RmwIDTYFypldmsd6fO3l+vRyc8+ZDjOH+WZEkgPH/Ws58PGv7IEDfVaO872ZLIX7ET3L6SRaZuFAsxIQGDhciF6UMqiLiucoR4VSmQIZOT1S3OOLVRdSgO16H57aC3m2zeiknSivND4JwFDwlhwrLYeVuvlmjkh+3H48y2XsxOTGrU6+NzN3sTD4FvSWupynBn6KXQujKCs0/ww7H694lnsygJ8sqmR2MX5xFXxJwKlprVD4ghnNOV1bDv5fuieyq0UOv/dQvRdX1vTHQFZ0qY+srRFDqcPhzO8SHYiGgjOMMgwvP/nmWMlgGkZnrDzOX/ft3l0XXVdej7jhaatlj3soL//f52x1JH010fndjQB79+o+6zpGEWJNh+1F18jgv+lY3lkLsPB9ZgZGzlho9d60lwSkfIqNlURoECHy34RBc845zmU83HhSX82UlQ7RHIx8Z1XMxeQS8ky09Lnl1Eb18yyWqKTxdCf2encN9nlVlsOeTsy2d5wZxoViI4KJXFpoew7VWrjPrUnWKDC+RWWryp76u39e8ZOcJVvEUgvGfwvLqv/L9Nu6onmU6Gvqb87ab5uEwdrqtgcve/JFJBqtBbMxkYC6D+M+h1kaAygg8V8Rz7BnOnD1G91Zu13/8zyr9rbCRcriqtY/sLGrw9j37DT2ctxm0fuj5mfxBw27/z2pFofW/urUUolTSCRFvXZAdaj9abj1cu7b9rd1XSI0o63cCq1hUVNdAx0dmGh4jqoPW20Jqty4R202nLdlt63xZ6gKDvBRONXtC6ju1mloCPSgaellFNfx1zjb4SI6WSAgcLKRbQvT0tIqqWthwQHwGT9YgPNqmMfr1H+gHeoXJA3ZzuSz/u82wlENB3nZEbfGZteGw6m8rM+DZG48wKBZKPw/uKgKBdhvnrI31916pAJWUK5eVaTtmxPP9piOwbn8RAAB8vfYg9/m0Ji+3pa8slCdzUOM7c8VbP8EJF5dXRBFYxaKmlkQ0SL/N8PSIct5kDr0pXpZI0XVl3/IBX2IycfUbX5xW2TheWgmbdNJ8x0gzcGRLKC/qXSH1N87JdV2W5/Ppyn2wrm7Jw0o5RgoyqwwsWJlkiqh7zia1X4sTy40iULavJ76qt6qKWE7m5RHOrMN6iFhupfV37mY5FkNgFQsWRHmEL94RnkGxPn6JsgNE/lPZsJh3hQhqeKMEzZSdcLS3W6Sd87kUUwsV8dwvo2ctsvvRu+RVlKU3SYreeTV389Go40Rgd/DV3us2jdPrPhfbaGUxWUslJJxxk6aQKHclmbXF1+aqw9b/ff4O2CTAt0K0Uq68L8qYHGbVGN3PQwYRb3ccdV5Zl5sQb0sqqzBerhv0/DwBkZ3dJXCKBV3jc5bjpfQUxnZloPV1n6/iy0NhBItD4RNf6WdhnPTFL7Bgy9Goe+4X051X1mhtKOLTnI58duoUhV4cBr372PWxWaq/f9jGF1xKArZrEG193KFxjhWFFTmX7qIv8fR+6nvLmVb/s2yPpSRqIjhgEKdDf+cWpf9mvJVD8ufrfvfrt8Rv+1ai3M2zvOCkZb8IPQXXiRTyThI4xYKKjy1JuzWmWaWoCQonizV7w9sWH6xLJewWuwxMx99tOEyN93/j1Pplk6Ml5fDWgh3wgU1fDydxbotemEVbjxkqcTxjUNh5065k5uguw1HqNvM3YEHoJRn6Oqn/ls3gLHrZSwxRH+WgVJF6JElI+6KVoE1JoBf6miUkttGs2c0tzoaS+ngt86q3l6remSe+jp6QdXnsO+q56uyuOtfo30unEh+KhUsQov8S620jjaQtpzSoRIVi8et/8Gvc2qBHdpA16cLTlZHZkwR1A53BeWUVNfDTjuOwihIF0MiHZOBzc1T1ajlQeAaOcqwfHyvxbq35jo9+9qxuGTszfbXFIrocEZFtmWVxunyDCpQOvtpcKvKM9TdTwu8pq0Mra/A12q6xOz9Rt6utOr45ZrIkhiTDZQQWBj0/V/W32c4I5bufmZYU+Z0mqjL2iFXLnCstVCHbh8uikz1WcoeoryfG9Ir4UCyMZt1abHXAQDwJKKVXpV4qbiW8+RiG/3URfPFzuFyJYzamtMywmAm1cUO+0KR3v/+ztfDIDP1lGi3LdpnvuxflxKWHV53D4u3H4fyXFlg+X/mMRUXZpNYj6AZZsRDYWUqqriVw1dvq+BM0CS7t1SLqM2XSPGPHxegv99ZltzWT3ex1a56ZapJMLPx/SXmV7rLvUU1wvUdN3k2lzFkKxYJ25/75g3GALRE8/iV7X6JHn6ejozH/e+lu0/OUV6w3kfKxsYZK4BQLPrMygdkbD8Pf52+Hi/+6EADCETW1yxNOENVRyIO0FW9yizL8cqAINuvsoNCy42gpzNpwGE5X1qh2Bpjdb0IISBLA+gNFkc++02z5Y+HDZdEzIDuRHb9aewAOa7Z28aaUN0Md2EdQmcA/cH7x837Yf8rejFRZv9c4sXQlR+20pJRQPoskIVP+TmkEu47V9zVGNRu9Z8dKKqCypla3jZkp8gkhCWoMKiAEYOexUrh52iq4YeoKw7L0aD9JHf9FKapeQjIZEanuzV6//1D6FxEod7zw8OlKdV8UaztDAqdY8LJ810n4au3ByAv+8Yq9cJPFrZWsMw5qBxAJGe5eA2JVKmS2Hi6B5MRQZE1dL4uoFu2Shxy0iGdtXuRd2X28DO6evhaW6TjLuYnWhGxGZXUtjO7ZkuscHoudjN4WU6eco3n0Lj/M3swURdb3+P9W72c6zuiatxwuidpmqnznzKyoCSHJVPk4UlwBK3afjLJY/fdni/Lrzcyp54Y/fVxnkH5vcYHuMpBRuSyczfl+svCvH3fBBsVEi/Z4Hvqv2nrqhzbPQ9wrFlo2HCgSoiHzQoA/EVHkXBdbXVKCQrGQrNUtn3PjNGdjY+h17vO3qLdCKi/hvk/X2qpTAv1Bkva51oTMcj9zs1Kpn+udajfQ2k876hUwahUCmh9zEQwH2nkdhC3HRKwUCpEpZZfXKdknyyoN1+Dltswqn2r5yuSckGTs4KnUAco5I5nqlslxrNkl58/cbBqB2Cq0/ElWUCpub8zdrjupYe1PRT0HpwicYsHTL8imSeU532047Ih2yFKm381dBAgkJ0hQGXHWM8/cRMCGwxUB+N+6g7DzWP2WwJnrD1FjKTCVx3B/lb4ce0+cph7DajbV27JZUl7F5P/CwpbDxfDK7K2GSunIHtFr+2ZYVXKtUlJexfR8nH5DtInkdOVQObSGmb/liOJ746UFGXlAN5txs+xeZN3Cqb1Gs6UQJVasXzxY6XuNTnn6fxth25ESJt83Jwfr1zXxRfRg3RTS7fFZut/5gcApFryI3P/P84B1G5AfW4mCpMRQZGYVDoxUCxsPFukeLyJOgDLHyd/n74AvbQ7KEbXI5NkPe5nu8Lh053Gm8vVYf6AI7p6+ln4u5/P/ZX8R/H3BjsjfD36+Dv71Y9jZTY4j0LhBMl+hoK+EOdU+tx0pheUMDrZOc/uH1nfwrFbsfGK9T1VydGCTVmPHKqlVSg5q4ku4rUQCALwxb7ui/nrMrpNX1qk/7Yath0vgtE5yQCVKfxczeJ+H1sfp81X74cHP14XLAmJanvJ7lu3DXhNYxYKnkYiCtbGpHZfqrCaUU9nNnu4RkqTIdUoQ9lH5/b+W6x7/R2XiIU5oA/+mQ8Xs90XrH2t1qYnjWEnSr8XOc1J1vjqfA4Tvz+4TZXC6shqGTtYPGGSGXqprWhsXtRNKz+fmu/WHVH//b51xLgaR78OeE+b9CO3y5b5fAsmwX5C3bit9F2gz54LjZVERNLXPfvTrPyq+02/rx0sq4b8K3w4CBE56GNSO6Pwe+cz/4ygVOdtuv2e+j6QZqKyuhZKKathxrBQ+VzwD7fKsFuU9eGmWeTwVrwmsYsGT08BoVs2rmZo6MgqaHKzfr28lcAp6B2p+f7Te8MozWBw4JU3lfu5n1MFu6j+XKJqBE46QYb8XdQIokfNRmnxOT6A2aCxixS76QL2/uID5WKVvFj2ZVp1CrnggtNmn0swt87d5O+CJrzYwL5caHbflcAncXzdbluFJnPXG3O3mB1nEqhIhoo1/yJGVtIRhfPllfyEAAJw6XRXJ2aJMxBaBhI9hRbsVX+Y9jrbqNIFTLJSzadZjjajWvPhPfLlBd32cEPNQxrRANwD1OyVYFY/fvl0fMEvvMkSbOMMOpuqlBJZ+wEgKK2Y9q51PvQMc0XyuLV/MSKktRXsfXrPYQU/6ot5jXKu0SSBBLalvRyt3nxSqiFEta4Jq0PcRUP7O4IchcIprdXkzYrEweUe0/YsYSfgw39WlRpuPxC5qndu7aQPPdvOS8mrYTok6rIf8mGlt80RZJfMOOwD990S7xOUlgVMsZLQd7tuLdsJTX7PuKdZ/uou2HYPNh/Qb1C3/XkUvUTV7jW4ZTuQs+G7DIfODbOD02qz+S+Zs59ObEujGiNfnbosom8Xl1ZF1cy3a4EK0joB3UNTmENA2raveXgondGY4MjPWsG0ZBKBbqESN42bJmPwKzRIxd5N+jBXl/ZL9MrT3VbscMnfzUTh1mn25wujdlOXt9eTsqFgunmPWliiXRYi5j4ITVHFEnDVKIjbmjfolLKVi9fT/6OOV+x4x/ARWsdA2tG1HSgydDLX8dsoS+HrdQVhZoHYo084weNuzkUVC9LvB4rBkF5YX2u41V1bXwjrF0o91c6lmdl8nmFa8kvJqrhDHszcegVW76x33aMs7hADsrttlEm0hYa4qirV7C6PKkiSA9xfvZi7j3k/XmR9kgKhOfZpOPhl1Cnf/orwNmw7Rw3+zQlsOOVFWyVwWi6m+pKIaylxIkMcD7fJqagl8T1tCUDAvags5gTMu9H9lFdVwgiEJpRbtOxMKyZ/Xfzb1p93Uc0M6HaqffFESvRZANIKMsrBqzylqfgsAa34SWhMfrREsLzgJV/ZvxV+4WxASNQrbvd8spk85dLFSDLay6XU99615tkejTIna+uXookYHsq+NW0d57osKBy+RO59oCteRYjF7/fVQL4UwHO+cKBGe+Goj/PXqPuH6WNujix2/ckm2rLKaOTtqaUU17DlRBj1ys5wSLYryqnrZ5HukVM7LKqrhNhMncG2+mjmbjsDDXzgdop/AE19t1N0lxvO85UlPWnKC+bG6+QH9o1kE1mKhXW4oPlMFK3efijqGPf6b4jxVGYqzDE5T7lJxwpTFu772lmKLIi8RPxbOC9lwoCjixFRfFn/9lQKyaQK4a1I8ICikth6i06brQdst8ncbbYkFL7tLvfZJizhpJxCZlWvUyyuhZdexMl1rkLbu7zcehrF/W8wlRwljIjU9VNadOmk6P0rPBMrK6coaV3a6/Pfn/XBQ4HJScoL5kKx97M9+swlmrj/kK4tFYBULLXM3G2/n4YFHY5RfGu0WOZZw0trSPl+1T1eGc+u2FsoxDIzYdLAYXp691fQ4MyQAppdXFvXKKUvCOT9svgAz1hxgSirGyksC7oUZNN8EEyOHJZyetZzwcFuiF2jv52Vv/qhzJH2gZ9U17D73sb2Mw7w7GfzpJoERdAkxzviqv5Vb65CttBbaY+6mI7CN4qhp5mNGGyfqI7Jql2XrvmeSSH3ujqOlsMUkwJrbBE6xiOy6cFB907VzUKpcUeejof2qtKIaisuruLbFsjhasZj51+yjL/GYonkZDhaVwx6d6JQ0Kqtro5cRGM7TW1MURcFx+46zVCXB5Bzq1jPeeikVaz+KBWcvkTg5c1OmTgcAaPdwfXIt5jg21AfC4KvEXaZeOZrt32ZBqUwKN4saykvvp9TO02Zb2msJQFW1cw/9ya83msZOoUETW/bR4x2famoJjHvrJ24ZvMKWYjF58mSQJAnuueceQeK4i+7auAk82T3rA2BFH3jv9LVRSYi8iILHgxXpjLynjxSXe+LRDVDfGZvdcyviLd91IsqxUosXOWliHT+tI2th3jhq8RLcuHJZtorqmoifg5vvp9WanMoVIiNqm7UyPQEPlTW1sK4u4JaZP5cfsKxYrFy5Et555x3o3bu3SHmEwbr+SHseegFIZPTWK9fVBURRIq+70oIl8foKyEXsOOq92cssG6KSxTvqnZu0L+Pwvy7iTul9SffmXMc7TXgnhrq9vTl/hyrvCKsDot2B07Bsn3U+VmBz3hR3nUz1GRwjNwuz98XKoymwmLeD5f7c++lauO+zdVBeVQPfrmfftn6kuJxpOVaPOz6MdtJUSmvWrRccL4OBz82xXL8eczfrbx/mQd6pEr0UUv83T9Znv2JJsSgtLYXx48fDu+++C40aNRItk2UOFZ1RPDh36lTtXaf4ccjOirR+hdfEL5sEF283zlXhFH9TxPfv8MhMS2VsP1LKtENGidYpjtXDnTaQnqL4CJg/BvsDld4OI7vI1pbpK/X9b+R2JtLJDGFDTn/d4ZGZUGaw/VHpwKiHBO5NTI+XVsL/1h2E9xYXcMVr2HiwiGk5VoscY4i2RZ7nmkvLq+F4Kfu2XFaKKXEo7FiXjZR8bej28PGWq/IES4rFxIkTYezYsTBixAjTYysqKqC4uFj14wQl5VUwJH8+PPn1BgBgfxB2FRCzqICykyFtlmDVi9zsPc+fyfdivzqHL5KeHcdJK7E1tM9owVbj6KZ6FJ6uUsUhsZ6mXvM3h+Lx045opVAvCyovmykdEoB7SrZTeNmpslQtP79/LzUOCX3GIGbEm/PNd9YYWWHNLLRG3xtdozYyrlGfZcciZrRrhQf5Mu3sCNl0sJg5ZcJ5nXKgQ9MGzGXv1sk9Ux95k34PWfoIP1kkuRWL6dOnw88//wz5+flMx+fn50NWVlbkJy8vj1tIFuT2f8aneeppz9ysM9CGmJXLqK0lhoPiOz/wmSKVlgh1fZSdDJIE6w+Yv3SzdKJ+0gbh6lq6BUL0evp2B6KbAtQHpZIxeqqLKYqFHZQpsEVtw411RDsT6rFwK99OMydzqvBkxzytiXBaRMlTsaKAPnkwu4Y/f7zG+AAL6PUDtIix8ntoxWoi89aCHfDK9+qdYrS+mgCB0opqzqyolsUyL9u5ornhUiz27dsHd999N3z00UeQmprKdM6kSZOgqKgo8rNv3z7zkwTw5nzzPAysfhissD5Y+eU0s1hotXi5/BpCbClQqzTxPPQwy+PRLCNF9zu99NO0F0vE7ggreDGTP1YiNpjUOs3MKhzDTOv1L7RKT2G5lFcEbh82und3fcI3iMqDO49PkRP+Ap+uUvfBRluIedoOIWC43GOVfSfr71d5VS0sqIuySYsYK6Kt1xL2KQ3vNl65XG3XatYXrTFxBBcZTkEEXIrF6tWr4ejRo9C/f39ITEyExMREWLRoEfztb3+DxMREqKmJvskpKSmQmZmp+nEEzYP6kcEPIbwrxJ3RRTljWBKJ1MZZd10Zk78zT5t7/2f6YZpnrKEnUYuqzuR7q7dO+/LrbRWTl1x41nh5OFlWyZjUyvj77UdLVdtujQ7X7gJygv8sU5vkY30pxO9YjQTLgtKRXAJ/74pxy2J2I2fcjDKOLf0AYb84Ea8Mz9LEByZLQZNUUUSjy910qNhXEwiukN7Dhw+H9evVYVJvvPFG6NatGzz00EOQkGAejjRo0NJh01B+JVsCqB2+QYvmCU5EiwzIC3X5RvU7/+tHXxIyPkfU+quWx77cAOd2bMJ9Hk1ebYIxI/QGh8LTlZCSGP0OOZmmOuZwufO0G1VSSSyYwe3sSPhkBXt2UDeZxBna++e9hXBBl6amx0kgcT9TPWVDGRV6/L+W654/+bst8MkKutXfT0onl2KRkZEBPXv2VH3WoEEDaNKkSdTn8cLd09dGfud9sLzOmx9xpPUVgRMNlVYmr/OoHbRKwetzt8OSneZRUJVUVpt3vlasBH2fmQObnxkNAOFsvDLvLd4FnZo1ZCrDP12LONRRFM2vULs8ZIdvftHfahmJUcNYlh+fzfuLC1R/ayMCHy5mX7apYHgvnMYp69wBWiZis9QPFoQpOF5mGDRR2S/4mcAkIfODtmYcQ6D+d6nOqOl0REmnOVzMv31x2a6T0K4Juxe1SGizhXkW9qcfZ7Ac8cxklMsjcjv+ULGccaaqJmZTigcZXi98R732TYrWG+S2HzF2aNabHVNFcLEL1rMQ3vLBKmF18Fo6tFh53r99e6mN+iyfKhzbisXChQsFiGEfP91UGqyKz1/+7xeHJWGH5l0vYmeFMuaCU7wyeytsOChm5hr15Gw2Nu3ptG2KSqWzqobAVkquAlbc8iNyDL+/24ztQUQfpVdGRbX3iqebk7uBz82lfn5IYKyWOZvMncr9MKGVcduibURgLBZ+gMVisefEaW4npxMc6/ci+dXfYyc2vZa5m49EdTL0LWPm/LTjOBwqKocJ769gPscq9UmKrJ7vn44OUePk4G+2K4AlPoJVYl5x1cEsAjOiT9wpFi/O2gLbbcz+jDB6SeVvdlsIwztARztHjGFJ0MUyDj/9v00wqH3jyN/VlF0qynKmLNwJffKyueSKOsb8ECq0de5Y7/Y/WLo78ruf9CanB1SesPlm6PkFibyffno2dmHZRjr2b4uhfQ7fsi7LbsUgEBjFQm7TZo17ykK184vIroHFYkEon4nCD+ZQP0ObtVmZsbGEXz5m4H/C8tytDlq0/BFm8Uj8TnmV906BNOT25NTd/acm54YkRdfFOkn6m05kz4AaGyyzo26pd7lOgDAtsf5uOUXg0qbzwjJIiCG6Af60U6z2auaIFU9YiXTqB2SxRYr6PGeIdz9T7aOOvJgzO+02zvfzeFQwtehGcclrP3CVqUWM34d/noldzlACfNF2hCDGxL1iIdI0ZfR6RdbOFZ/5OW02LVlXrFNTG22fYM4pI1waf9WHWMepcVWrRB0vrYjMqEWhFf3LtQetlxUABcNPzpixTGAUCz80atbtpk7Bm7vAiHs+XSusLD9xQuOQxbr33uzxaTukU5T8Cyxlye04FqwriLPQAsP9e2n0Z37BB12w6+BrSicwioXfkQceJ9+9V74XF2iqNsZ7Cb2Zh1PXpY1CapTLZRFDdlbsr2KHJ77a4F5lgpuviAmZXMJfXQx0JyNyMhUUrnnHeiwMUaBiIRSDXSEW319a5kHEOlafg1kGQx4T6vs/Feh+J5cS68HT4omvbCwf8CJaLY7t6QPADVP58oaYIXSXjLiiuGB1PHWSwCgWTji9ccvgQEuaPCs4jnd+wKrFgicXiAhQr0BosCYQZMUsayYTsa6d2OBQobiAXEEiMIqFzOyN/CGaRXFRt2a630UUH+5ScYQRyWVvLnakXFFKpVzOFkrUUwTxI0HSK3ivxa2MrrFGYBQLP7gEGOXA8IN8PBglwokFth0pRZUMQRDEAwKjWPgdq86bLBHgnECIiTSOKBTlCxNjCiiCBAmRuwuLzsSvfxwqFgLZclg/2NYXP4fXRmPNcoGwkZOR4rUICOIJftjq70eORgU4ix8Co1j4IbDJc9+io2W8cn6nHCHl/LDdfCsq4m+u/EfsJu+zgo+CofqKeF6KDYxiYYX9p9wP1cqrAOFsIDYQtYvjz5+sEVMQ4hk/x9kyolHMlnhGG4zPTb4UvHuIl7hWLB770sXANnWs2n3K9ToR5/l2/SGvRUAQxCanTotTBmZtPCysLF7cy4FFJziKRYxM7HmDl9iJ3R/vlFW4N5MyC6CFIIj/uWnaKq9FCATBUSwQRAPuMUcQBHGfwCgWMWKwQBAEQZBAExzFAjULBEEQBPGcwCgWCIIgCIJ4T2AUCz/EsUAQBEEQr/E6hkZgFAsEQRAEQbwHFQsEQRAECRBe2+8Do1ig8yaCIH4mr3Ga1yIgiCsER7HwWgAEQRAEQYKjWCBiSUrw2v0HQYKF5LlLHRIveN3SAqNYYLIusVTV4P1EEJHgzjUkXgiQYuG1BAiCIAiCBEaxQBAEQRDEe1CxQBAEcQG0qiLxQmAUC3xpEQRBEMR7uBSLKVOmQO/evSEzMxMyMzNhyJAh8N133zklG4IgCIIgMQaXYtG6dWuYPHkyrF69GlatWgUXX3wxjBs3DjZu3OiUfAiCIIEAraqIW6w/UORp/Yk8B19++eWqv59//nmYMmUKLFu2DHr06CFUMARBEARB+Nl0qNjT+rkUCyU1NTXw+eefQ1lZGQwZMkT3uIqKCqioqIj8XVzs7QUjCIIgSJCprfXWPMbtvLl+/Xpo2LAhpKSkwO233w4zZsyA7t276x6fn58PWVlZkZ+8vDxbAiMIgiAIok/DFMs2AyFwKxZdu3aFtWvXwvLly+GOO+6ACRMmwKZNm3SPnzRpEhQVFUV+9u3bZ0tgPTCqHYIgCIIASJK3Qb251Zrk5GTo1KkTAAAMGDAAVq5cCW+88Qa888471ONTUlIgJSXFnpQIgiAxDqYdQOIF23EsamtrVT4UCIIgCILEL1wWi0mTJsGYMWOgTZs2UFJSAh9//DEsXLgQZs+e7ZR8CIIggQDtFUi8wKVYHD16FK6//no4dOgQZGVlQe/evWH27NlwySWXOCUfgiAIgiAxBJdi8d577zklh21w+RJBEARBADz23QxOrpBdx0u9FgFBEEQXnPwgboGKhSBKK2q8FgFBEARB4p7AKBYIgiAIgngPKhYIgiAugEH8ELeQwNu1EFQsEARBXAB9LJB4ITCKBUa1QxAEQRCAxAS0WCAIggQerz31kfghp4G3aTRQsUAQBHEBNKoi8QIqFggSx6QlJXgtAoIgosE4FgiCeMXDY7p5LQKCIACQnixOyT9dWS2sLCugYoEgiKP8Y3x/r0VAkLjiTKW3ASNRsUCQOObfS3d7LULcgC4WSLyAigWCxAk9W2VGfbbzWJnj9eJmiDDovIm4heTxFiRULBDEQZTvd/ucBt4JAgApieioiSDxgNfKPCoWCOISeY3TPa2/phanzF5yvLTCaxGQOCGEFgsECS7K19vr6LC1aItHEMQFULFAEAdRrnU2Sk/2UBJULIKEyK2JSPDwOuEdKhYI4iBKi8XvBuV5JgcAQG2tp9UjCOISmN0UQeIEr1/2kT2ae1o/giDu4HVeGlQsEMRBvH7BleRmp3ktAoIgcQAqFghiwLi+ubbOV1op9p08bVecmMRPylVQwFsaPEQ+U4xjgQAAwB/OaeO1CAiFJy7rbuv8Dk3rY1eUVHgbv98r0GdUPJ2bZ0R+X/HocA8lcZb7LunitQiuIfI18VrxRMXCI5IS1I8+PTnRI0kQI+xq/tnpSYIkcZZnx/XwWgSEg39ePyDye7OMVKZzrh0Ue5OX5pkpXosQk3htJUTFgoOMFHGDf8ssXO+ON7yOY2GEfyVDaLAqE0oa4BZVX+O1lUEkqFhw8PWfz3Os7CA1KqQer3eCKPFKEq9nT0iYs9s39loExCU6Nm3oaf2oWLhA4wbRgZGiOlvsfH2J3ceSk1FvyvXaocoIHxtTEEGM6tHCaxG4EaWY92+TLaQcv9DQxHp+03ntXZKEDioWHFht4p2aeas9+pmmGbG/hjrKID7E2F4tI7/7V61wFlRaEMQcntfEbFnV62XXwCgWfu68aANKtMHC+WHnGR866I3pKW4WNe/+C4SV5QQ+bqKu89JvenstAhILxJE2znOpfu9LAqNYuIFIS3aPVlmqv/u0ztI5Uhzx8I5ue26M1yKo8NPqh5+WYq7s38prEXzFK1f1EVJORmrAdpcJGkH9PhDz4ueJNAAqFq5A6887aZxrxvRqCe9eP1C3jIu6NnVGkICRnCi2SYu8ZcG/+3Ro9zAxAbseJxjQtpHXItgmLQl3r5hxz4jOXotgCL7dLiBimUPEbNOPA5vfNe94wes1WV7aNE73WgShyO/myO728rn48R3nRdXVBeGCHMAsU7LXb3NgFAunJuNje7WEkMWyfz84HJBGhGxyEblZ/PvXY5GsNP7AUjkNvHEEZZ1hef2yW2mGn98+hPnYJpTdT07hdVpopxjbu6X5QQa4fVes+kgZBb6SACDBaqfrEV2au+ugn+bzmCRcikV+fj6cffbZkJGRAc2aNYMrrrgCtm7d6pRsXDg14XprfH8YbfHlkV8NmmJhVdlYMsl6+F4/roSIlCnLgSiXLNamczvmCK/XbfReHx4FT+9Z1vpcB+jeMtNrEWJWTVKGrOchM9W4XUma/61waa/6fttpg1zIxc61VXaaaUJBrw2QXIrFokWLYOLEibBs2TKYM2cOVFVVwciRI6GsrMwp+SzTPsdagzdCHmR41/FZl0KMzNFirB4+1Cx0qGV8M+wmCXMTr+++830N/Qr9vmZ+1/BOXotgmYEanwqzNvbd3ec7J4wgJEni6u9Y+vryqhobErnPY2PP8loEW3C5EM+aNUv197Rp06BZs2awevVqGDZsmFDB7MI6MLGgHZDPapkJ6/YVmp8nqf/3Gr/IwQLt8aUnJ8DpytjqIJSDaqzOSlnRbV8ut7tW2WlwoPCMu5Xa4Kah7eESim9FcmIIKqtrhdZ1lg+sMwDsfRGLb9k5HZpAwfHoya2y395yuIRZNisI8YFTlJFi2wk9huNYFBUVAQBA48b6oWIrKiqguLhY9eMGQvsyXygIwXTe1IOmGNLWXc9u5+8wxX5fCxWJXvtywkxsNG/g3brphtm4bRN9Z9PWjdIgkdK2bxzajrset4cT7aSrVyu2bfNG1lPe1uJkv7z5mdFMx4lwCYk1B2ojLCsWtbW1cM8998DQoUOhZ8+eusfl5+dDVlZW5CcvL89qlVw4uWffKb8iI5mVGmxLiw6csWSxqKEszGvFL8i/VGiALSoB2m7qdMel177c9sPjbec8PiBWE3kZKVeS5PX80jray7pvJFua839N0N9aDxJA37xsAAA4pwPfxEGVNVpAu2OdGMRS3+oGlhWLiRMnwoYNG2D69OmGx02aNAmKiooiP/v27bNapSFaL3GRnajVtB6yVk5TGHj9HUYpBtClNhw4jaDNmryC5YX2U8CnWMDojRDxuui1aSd8e0Q++p6t2JcH9JSQyVf2MjzPTNykBAnyGqsd8vq34Y9J4fYbYbU+s8nR57efCwAArRvxbSv+88Vs8R0u7taMq1wz2jYR79OnByEEOjVrCE//yn+RlGUsKRZ33nknfPPNN7BgwQJo3bq14bEpKSmQmZmp+nEDbQdgpyPSDmCsfbA8mxbxsicIXsOj0cjF7YJm0BK3xSJOdvS8imCyw0Gp9MRxW//jVZJ4UpDrbXP93aA2hucdK6nQ/S47PQlSEhPgx79crPpcTho2XPAgKBTNwxXxqK2UQd15Z1uSet7+Q3/D7/84rINpGSzHyBhOAiC8W2vCue30j4mlXSGEELjzzjthxowZMH/+fGjf3tsMajw8dwV9ueZZg/wZeuFxWR9a44bhwZG1Y3V+V4iabAHbM++4sKPtMlihWn58YrWQFaEVj6itSU6KN+lSPs9xI1n0Wh5PB6X3LIzqtRq3wauO02q9JRXVut9d0dc4vHlPRr8FPVpmpRrGjRDBhqdHcR1v9N5aeadzs6K3Xxop0hIAXHdOW8Mylf4iZtYjFqucfYfM2IHrSidOnAgffvghfPzxx5CRkQGHDx+Gw4cPw5kz3ntgax+sdmZh1Mj0EiLJ5xjFo2CTjfKZD8bDRQ9eZOv85IQQPDS6myBpAHIaKlKMM57jld9A1HF1/2dGxX1w7kHzXntLSufrBkb+Be1dNCHrwRNsi/XI3gy5fzLrJi5mAylvX6GNy+FkjAXZSiWn8W6QkmgpuJ1ourfMNPX3eFYx2dRu2wUAaJDiruO1WTvIr1tuU77253Wix9Dx2meHS7GYMmUKFBUVwYUXXggtW7aM/Hz66adOyWcZbZ/76370WQEBgKvPpjuU6j1n3ofm9ay6VXYaZKYmCpXj8cu6Q2KCcXmje1h3rPTyng1qz+cwFgtRIJtmpBhel56SwvMY9I41KoLXahYrS2RfTRxqesxHt5xjqw5tIjd5KyntOTilf2sndA2SE+HrO+nX/tkf2aK48r76tGtr2yQ94p+R0zC6zWhPYanTamTZyI4gk0pKDaxaAPVKnLK/+fCWwQAA0CzDm6jDenAvhdB+brjhBofEs462sSUmhKARZyd2vLQSACjbHDnfUq+NE3JiItloE3FWc3g8fPu6AVzHKwc3ViuPEwrIsM7WI2n6wRLlFbq7Qgx8QXi3CzsRbIvHudTPWwLvuyQ8Q9eKqH0umQIzoNJ9G+j306vMq4khAUsQJk3E6L3PFmzBoTVBrZXI62YamEUfllljO0qENvkBGDnWhCRJZSob1L6x6focAEDbxuH6lI1ukKIjbcqhZdoZr7rnhhWJS3uF17PdDD+rR4/cenNt60ZhE73ZE/zyT9EzIac6+i3PhvevM+8AkmOd1J0hm4bt3OqbzzP2YWIJAtci015uGbMqWilCC9MGlMcv6x6VydcOrPfTKG6EFrnvuGmouc8Y69ZUJy1uqUnhbnvvidOw6VB9XCC5Rq2ITjvtWsGpu0NbtmZ5FMq2O+IsurOsshxW64uSphTLiQh4t+S6gf9anIOc06GJ7ncPjuqq+51WaenUrKFp8JpB7RtHtkyqOhmDRm6nM1KauV/4df3Wt3VPjARCwmWnJBrP9nhqtzugK/exy+bwkElnQFMMaXRrkRH1GY8SBwCQanFmLMstQuH5k4Fj7Bu/6wu9WmWblpGYIEHbJunw4m+Mt0Na5aeHLzb8/ubz2hvuNrLa5G8w8ogHYsl/KCfDvOMXuo2d2X9HfeCWZ8cAAMCqPadUkWhpks2863y4rE+ubeOk1im5XjZ2WG+dnYmP1etU9vE0C4f2GcjKHRN11/MHymRUu8XYCme1zIxqS17PHeNKsaA5GsodBYsDlXwM0wuiOMYNR0RlbPl0ZQwIKfzSWGlo/73DWCtnKVIvSVHLrLQoJ7OuFIXACo+N7R71mdIx1Izxg82tUVqa6GRWder9Hte3FbO/QbsmDeDibvbScfuNW85n37qnR9fmYtqbHjkNU+D+S9gCRllhSIcm1ARqyjaXnZ4EIcl+O2ymY/mi+nPYVGO8tqfK8hs5oWq76i7NM+DyPvV5i2hxTSRJgh3Pj4n8nRCSVEuByuBrRkOB9qseuVlwYVe1lQWXQhzirJaZkWWHLxkcqXixa+qUwN7DH2zgiCeLdnmfXGqkwIFt2UxnAxiPM+L3Jvv7AQBuPb8D5DRMVs0KRAdV+uKOc5lNhlbieZzbqQl0btZQYYqlbY3lK5O1jfVvk81XsEWeuCxaYVNi5ZVgfQfuHs4W+MgMpX+BPICkJ5uv/f/hHPN2TMNu/940IwUkCWBsr+htuRPObQd3MdyXsMXSpiA60Nqol4Na5P1T9SXRyFaRyJIl5ag3ftevvgyT+5ecGIIBivdQFddEcUMSFctS7XMaqGLRmGUspRQHAAC/HdAaHuHceu40gVUsxvXNhScuD3eEcnhYGvIzMnvvCKk/hpazwogOijVmUe+3UeQ4eU312kF5qoYsc1udP4lfdjOM69sKBndoopJHdEeYlpxgycTKMrhLdf9Yj2euW/O3bHWRY6/I9+sBg2U8kbQwiZbo9UyTl6E6W/WM4Hm+ItvwW+OjAzSN7tkCRmtC2tOC+RFQD5yi33qtU6ZuTBTGmtNtbPOk1UB7Zo/WWXhp8Tf6tWkE947oYrgcSvXX45AJIHpyadRehnZSLuPr30fZyuK1k3FgFAsn76NWm2XRLJUv0X0Kc6jTa1/f3zsMRvVoAf+787zIbhClUuSGDFZQijS6R4vITMJSWQ6uN+abhG6O1Bmp23rl2lMXPxT2G7huSDvV5+d2NB4gWUQQEtLbhYYlsoprdLaZG8Gzo8LKPe2X18j8IB1qagnUEvqyp5OPRlk0AQJNdJwU5fvRKD1JV57rzmkLGSnOx8Ew6l+6Ns+Au0fQLEH1QtOWSey+Q0ZW2o9uOcfz5Q0eAqNY0GiemUo1ISqRHxbXfn2GY/S8x+XGk5OREtkZse6JkdRjHzcxPdPISkuCUEiCXq2zIs6aGamJqrDFevLb6XzOt7FFU0Z+Fm9fN4A5+Y/dd43XxG50iwjlezv9uYjlIK1SyX2+0opkcqyV9sOePpu/bB5YnPGGdmoCyZzRE43Epl1TG81uFrOto0oWbTsGSQkhtY8VrRCBKOVZ8MCFcFaLTMhMTYLbL4h2PJaPNbuH8nGrHxthW76C/EvZDjRpX3XedYbHyN/S/F5E8pRBjhBZWfZaBwmMYiE3RmUwlKYZKVQTot06QOLzCla+fLLp8pqBefDP68NxHrIY4mtoOxRtwzFK1jN+cFv4SidojR2UM1QWC0MuR1ZWJ5YUtP2rNsz7EMquIV4ptGIrzakNeK0wZuu6jNsIRd1Ls2LcMIR1bZ5py2dK716wxDqQQOLKhAogpoNXSvyfmwYb1vbq1X3guSvqrWoR53Sd8uzLVl9a+5wGUTFL+rfJjuz6iiwX6kjQTqNUNWF2uKYov3W/SJJk2m4fudQ8erBeu1E65/api7b66R/Vgc9YJ3Ksr+mwLk11v+vTOputEIcJjGJRT/jpGPlVTL3h7OizTJ6qtoNomZUGu17Q14b1gj3JiWNCoejtn7y+G0qMgn8lhCTL2ydZ+evVfUyPMetknZ6V0tZ4++TVh15umW2u+BhdgwTR7UgZWTIz1Z6JV6sYdWjaEDaa5GgghDClLRfhb2P0Dpltz2YlJTFk+G4DgOFDUr+XzjY4c0VMfP0ZqUlRoajtPlmj6zC7xrX7CqHwdBUAhPuhOy7sGDlHm1wtMy3J/o4SC7f0tmHm+Y5SEkPQMistaiu7crIwsF1jWPfkSMjQvOd6ofQlkAxy9Fi/Dw+O6gqdmomLHWOFwCkWkgRw74gukNdYfwbfsc6ZMspkqCA5MQTj+uZGygSI1riNIgqyNgtlvProLUNE8TtjgSxI9DJZOzpa8CEWz/pRJiG+VY5MTJLQiVqOqHuAd1FSKtvZhmkYG0MKv+DN67bpaXfnXNhVf9YRKYLhJtCsIBOG1G+XjZ6vhrGa/MsIIwWGliTKEgLHYiuDGM97KOqdZWkHN5zbDu40SBmuVPp4xLrh3HamybOMFEqthScpVO/knC8storxDWILkGVMQkiCUEiCfiY7sHjypFidSKWbTBInXtSJa3u9EwRPseA8XvuS/fumQQAA8PWdQ1XbjVR1MFSijHioffHkqI4A6gx7RhYLrRYcLZN7HplWk5eN7KE/iAuVX6eoczW7AGhV0pIR6XHnRZ1Uf2vbUpIil4r2+pIYljGs3pEr+7dWl0Mp6K3fW1kiNJaIRbm0QrOMFIVSz87se4ZFfWbfmdb51evOqtkmW31P/aqHKgqqTPPMVGEKztntot8NvbJ1LaiSZLg7i9eKE87fodxNJsAviXLPZZl5fF54YdqBJgF1p5/f8L+EjFh5eWiPUV6/6tZC3wHH7PF3ad7QMJOkiGUJvevlaed2+htCwl7oPNDyPJitP755LV2544H1ntA6NZaOI1vj5W60RMYjU/TWQf37ffXA1tTPRVq6TE37AjtZZWRN5SXwXA5PwDVW0fnSyNf//jTF4a7wTCV7YRZJSggJCFgVvpD/3Ez376A9d9okSZLCW92/+NO5tuSROadDY2hdZ5mWZZBTnat9Suw3zEQby9SsKGs4i+IAavYU05ISfJOa3R9SCMbsRVK+CE7s9xW1djr8rObw+8H0wDx2OgsR8h0qKud+2frmZUc5HMpX0blZQ+pAcAHDkoEZogY8SQoPLNcrlhqSEiS4pLvaEiMvkdnffsaOMoJfl+YZ8PYfwhYJAow+FjqyGl0DT/ba/adOMx8LoI6S61SAJyvPx+ojpWVxPVpcQT1W3T9p6zeXoGlGimqLu9X796cLO0KzzHqTekpiCGbdc37k7+2KKJI0tEvNhISdvJvrRfEEe/3axIs6RiaG/drUW1eYlkJMjglFlnAsi0eFFsdiYNtGMNgg/YQeKx4dbhjfyE0Cp1jInT/bsWytZGT35vDbAfUzQpbzeLbp6dE+pwFc3lvt5+EV2vDeIQmgaSbfOp4kSbDmiUuo3901vDM8eXn0rE7EZS+0uHSj5MnLu0fa1W2KhHVhr3PK7MykPLl1JJmknleSlBCKUmJopCUnwOie9T4URoqkds149+SxzPJos9caXcl/fz7AXO7Ato00WxKJcDdHmgWBReG2OhGhtRE5+Z6ZHLzm/ZTEBPh1P3VKdSv9x3VD2qoczCVJUllyzZbzrDwzQui5Qq4dxBdzJINzBxbrY+V5/HIeED1lSflcWVKym93PjNQkQ78/NwmeYmHgaWuVkT1awLWD2tQ7cZo8u2aZKaY+Eay4qVAY12VNELO02HqlEhI2dVq5j9qOmbb2PKBtI/hN/+jlA2XbefE3vSJlycs4yr40aguwotfp0LQhUy4Po22O2vKTEkLw7vUDdY4NH7z1udGqz2tNZqvK5/PAyOi8FspzlcXcQsm8KkkSfHvXeXT5dOo3G8wL8i8NDzbyjNHwaHPk65lgkMRMFFPG91dNSLSw5nqxQrVimVLbH4rsUkRH7yVAb6/5V/a2XCbteo2jXNJj8kgQrVhkpyfpOkqaKSFKC43edavKM/7aVwRPseB4a6y+YGZpoN+9fiDcdXG9Y5/X1gYRtG2STl33M6OP2dZAHQiwdbxWlxtys9OifF20z+mas6OXofRmj9pPv7jjXLWVS/M9i9xWlqy0W5hZzeAEwHBXAYD62rsoEnh98+ewMpGaFIIeuVlR58nl0z83W7aUmNLDs9QFEL4fz/+6Z+T3+orMy5UkvknLwHaNmfM/qOvhPsW8TMUFytcg+yNYQbnlV6iiYjOgWwSlQqy4oXpKgPKem2WBVtK1eSasshjIS/ucZUtNEMaLwCgWmalJ9QGYTDoiozVMI2T/ALMkValJCcI9d0W2NUkKxxXgMetKANBcsfRBizIpEkKI43EGrBBlptW5haEQfYlEi9HAKp9+M8U6wA7/fXzjd33rz9YTT1Fkz1ZZ0K5JukrZEI3IzlYvey3L7NvolWHZQqyErX3YR0/m//052rpES49OE/OmujZp1oWMMtgJRoNwOoTIR95yXns4v7N9fywaLM7YNIKgIFglMIrF8LOawaK/XFRnXmKE88HbaSjXMmT5tMogk+UGLVU1tVR/BjOUlx9+//VvCMtSX0pSiFsBExHm1w6WlzA156UmhZjbhJXQ7jJmSyFWL0c7iJrN8njq0Sqt4b/F9NJ23uGwKVy/d5l246CohFxm5VE/F7FlUiOnUZG/PFWfUkAvPboRRvK+c91ALr+dWsaAbgBqpeaGoe3hHKXDo8B1AyO/hRSDCMymSyHaPoHDUuJ3nNl47gGSJEFSQniWaP5AlWbB6IOVMzZRsj31q+7wyYq9NgrR/+qT2+pDyLLk2JC3idqOyGfwXUeT5SKAcARUPQcwomMSYQnza6dfNuqQ7x7eGXq3zqrPL6M47op+rQw9ubUiJSWEIP/KXvDFz/utC6sgp2EydKVYDAghwmKEHCw8E/ldmXtGxikLU21tfWItO9fy3zuGwLJdJ3W/tyq/nG1WFEZSWMnQa2SJWffESMOIsG5kyWzbJB32nDhdV599S+iPf7lIFbxOW97CBy6McpZkefaSpL6XW54dbSt0gASSqhMe1iUHMtOSoKhuG/IrV/WBBz5fpzknNgiMxUKmlhAorahmOrZpRgrVeW5c31aUo+1htdMycm6kHcMSFlwvx4QbjTYpIQRX9g/f3+z0ZN38GcSBnQB2kCSAey/pAp2a0c39PVtlMe3Y0ELrtnu2yozUycqFXZvB7HujA0KxLlnRxo/LerdUdcBG8pguI+icS40dElU23VJ07wi1s6ky2ByNAW0bw0RNUDMT8Ziw6onP+nzl45ZNGq6elbNiMFgz5SmiFVnXYPRCxvdrkw03DlUv4bFcb63OrhCqXJIsi/rzvMbppinPRSjbZkqFnnza7+v/llTZc4cJSOzoFYFTLBqmJEIZo2LxwY2D4GpF6mTWpE5K7r+kCzRIToD/3DyI+1wnMFNgdr1wacTcyfpqPTuuB2TWhaodWTd46g0kdw/vbNjBJyeG4NWr+5rWabbUwgpzdkMT3ExZfOdFYSdKERaACUPawW8MdiYY8fvBbVQWJVkaveUop9aUuzbPoFpItKmt7eTaGX9O20gSKT14trKLQmkxaJGVynyN2qBiqneJADx66VmG57/0297QKN3ceTps/IyWaUDbxvDAqK5qmRjuXffcTBjAEf3Wz7BcL+0Qvfce06Z7CMsap/zY8hqnq7JybjMJ+KJHWnKisNjsQzvRZySsgwxroBcA9qWQ64a0i1qyoPlYXdE3F+64sCN3amkavCbR5RSnMwA+5YR7WOIo+wptXIG6/41mJSIG6gnntoMLDLIhyvBsG6SF7qZ1eixLBLR6tZ/83x3nRt0/GmYzRCMeGt0NBrXX91WS3z+zonX9JjhkcWZXiBqz531lv1ZMZv5w3Akbgmm4+bz2MEmj9PQ3yc9her/qvldmIuUuow6RQdXUPmsxpDUwEDjFwklogztrc2BtuB/dco75QRbLNi+HryD5fshWgdd/109YFlUChOm6MtMS4e7hnVXR/OzcDt1OgFJoZmoi5DJkRAUA+NOFGhN8XXl6cSkUh9jCtNM3+T5dmSnTZDuc9uOzDQZqM1jrUJKdxhcXwloSMucHAI1xwTZRMps889/0b828HBF2thQz6dHjiz8NNfzerTFZuL5HuSFaS5PJ4b4lkIqFWTvTW9c3L1cvgpr4xn1prxZ1ZdNbk172Vp7GxyLzFXUZXo1wIgEaq8UiJTEB7jWYiTjJ5b1zYRYl0ZUeF3drBosevFD1Geu94/Gst1K+HqrIiwbHNUxNNFQq9Wfy4toOLWw2QDghoNYfgyqLyb1ikZTPQkY/NirypoV7FPVucxTx16v7qCybRtek52PBRJ2MZRU1Fgvgr85uV2Wlm+duNzonxJJRIzC7QmRYdoXwpLb1gpyGydBdE4xK+0Ios3DaHTyMFKPXFRleVbtpiL5FQcQLELUu7DG0zj0UkiA1xGahkaSwf0nbJg10y8tpmAzHSysV59i/fismXrle3YGP8vG/bxpkuATmZp9Ik89uhEhlhMSo7ywWrWuVEdLs64VyckCq5fCF0pPjeCk9Z4oRKYkhyExNgvIq/nNpsFyC1ffRaClE6QQbJAJnsXB7KEpMkBzIKFevHPVrkw3f3X2+8eGqMwXMbvyAjW1nVjqAtKQEauciSrdhucf1O5TCB4tZCjEu5Tyd8MVhKdgbRkZqEjXqpxfQUls7Icv029iXLeXkUKK3srNi2ahgct9qCQGDqPTG2Gjgj409C974XT9X38+B7RpBj1z26MPPjusJrbLTOHaF8MvkVwKnWADYn53wcNPQ9vD1neeZNnCe9v/gqC5wYddwR5SUEKKG0m7SMJkaOInnRWvEsNXMU1zUEt8a3w/yGqVDhxx6/A3i0gxQid1Os2FKomkG2p6UENy6yxY25OE9VZIkXYdcUXClQKdcwTkdmjAvA8jLr/J7bUS2YjeGKOuf9gNRS1C2gpfZuLbkxBBTzB4AtrbHcp8njTkragutERd1awadmxvH89FWq7RcZKbF7oJC4BQLSQLH7K60Fyg1KUF4MqFrzm4DvTRb37Q1pyQmUEM987ziiQkhyM1KtTRgOD22uqkcyvx5eGf469V9XKnL+J7LTpL2Ov8NT48y9HvQmz3r5ZBwO8Q6Lb22YSwNhzU+K3lLrPDKb9UJtyy9n0T9u+hlRbl8ngjcdneP3EmJQfLEZd1NQ97vPFYGAADndmwCt55PP9ap/oZ2f8YPDkfcpd2O7PTkSHCv9OREaJCcAA+P6RaR0U/Lw0bErkqkg1nnt14RvlYkfjFbudXwjDorESIQgbMqP8DScdXfN3ca07i+reBUWWXU5+1yGlCP532uNbXevBTae8061xDR2njukd6xzTJTYdkkcdYaWrA5noHU6JIIx64Quzwwqiv8e+luKC6vj1PUr00jgzPUDGzXGAbqpD9wsv/W3p3nf90LPlq+l9p/3j+yi+q9kSQpYnUUlqDNBQJnsQAw7kREpTPXo3MzuukrVjRNHpy8Ijse3LznPXrpWdBZJ6JmkJQbUSRwbEWUEdX+9SJnGmGl7jmUKKYANP8NvrJZD2+RJWb7NLVejgJ/M6C1KmWAlloDJ24n+OWpUY6U+5v+reHtPwxwpGweMlKTVEthsUrgFIuUpJDpurKTzLnvAkfKZe3AeK+8hhAubV0+9KqBreFBTWQ9kYSzm4aZeuPZUd/fPdw4xTcPtw7roLt9N+hwzbDr/jdLHLfhaWc6fwAwXHY0280iArdtMFaXd1TxEGwI3bhBMvSkLI3JFg8rIbhFINrCcPXZeTC6ZwuxhdZhFJPFL5Zu0QRuKeTvv+9vK7SvVbwySGg7cV45jhRb267Vu3W2pfNYUVosLqI4vMnrlG6hWrN2eHjxorNhGcBY25YczdaJ6zAzywPoPB9OYbTX2icvGw4UnqFqFlYv000rpmMJ4uzEsdDi0P1ICEmeLcsZobzcU6eroLyq1vB4EXE43ILbYvHDDz/A5ZdfDrm5uSBJEnz55ZcOiGWdrLQkVZhuoXj4UPUalPZaneysJHDvFpjOhFwSxKnbaZYwC/GOllmp0DClfsn0vE458FvGfCtmzYW3OVkdDjs1bQjz79exngr0X+rfphGM7sE20/dqdi4698jYXi2Zj01KCBla+JS35D9Ld0d/r/HCjZWlWW7FoqysDPr06QNvvfWWE/LELKap2l2Qwa4J+oAiLbYd7HQgKx6td1qz7GNhvXr9MgUX2oWS4pxW1+KHLhJbsQb5WV3QxXwbJHfZit/5IsKKH4GYnDfrZFw6abjKx8EpWG5JSXk1VFuYbYdCEnRoWu/v5ZSCfHmfXLhvpHNLoiIQfelvje/PfOxfr+4D+Vf2on4nKwmf1vmwMFlVYkOv4F8KGTNmDIwZYy1ZF+IsjllqFDg96YhksdQMLs//uqfDNTsLLWgTy7GtG7nj+6Hd3kwjVmZLWoTl0gHzZbDrzmkLf5u/Q0yFdXy8fK8lp1UZv0Z1zExLgjsv6gR/XyD2fvkJoyjPj4w9C9KSEqBrC/1JRnpKYlQCyFjAcYkrKiqguLhY9ROzePh+xlKXLmS7abikyN/jB6uXDmJ1kAMA+MM5bagOxr6/IhsCymPbi7+hz95U1Rg0IMO8FTqfv3ZNX7h6YJ5pvWbcOqwD3HxeB8Nj7jHJSaKUn/U9qaoxXns3gxZ223aIcwF9YWpSQlRqdadgyXnkNj1zMw2VCgCA2fcMg2sHhf3J/Kke0nFcscjPz4esrKzIT16e/Rc8FrE72J7VMtMX26HcwizwjlsBtGQRtAGH7PDcFb2i1l1FpJq3gt6l/PGCDlFRSLn9A1TbTcP/t8pWW2DccEY7u11j6q4f1ucoR0A8p0MTGKzJ2Nq4QTKkKp6dE9dj1xm9qqYWkhMUZXikwYp8Y2NpkLVD4wbJkb4B41gomDRpEhQVFUV+9u3b53SVzuHhU22QkujYdig/Qgvqo8RoPdJvntNTbzgbmtVF09PjodHdIp1l2yYNYIxLz1rPTD5pzFmO+Bkon80PD9L9R6ya7u0+dj1ryKtX9438rlVo599/AVzZv96x0wnnabsBqIh2OwHl9l7RNzdKaQoKfusPZHjbilHiR7/huGKRkpICmZmZqh/EnyQ7tJbXoSk9kqMRPKGC/c5F3ZoZeoYnJUiqGWX33EyYElDrlPKRtmmS7shOAdFlGoVFz0hNsmxRYF3Os5/qO1pJj0rNLklcVgCv3Tbc2l2DWCP2vEJilCBE3rR6CfPvv5D7HLOQ3sYdW2zeazek/vz2IZbPtTNzFhbqQFA5AABJifylDWjbGCbrePmzoAqAyVh9I5uRGJmV9BgafWNIVGb+EKAt6NyKRWlpKaxduxbWrl0LAAAFBQWwdu1a2Lt3r2jZEA+ZeFHHqM/cnKWYBYMxEqV36yy4S1BkziAohErO1uRKcPKRWin7twNaW77nvPV1bZ4Bv3DmDspKS4LhZzXnOmfBAxdyHa/lxqHtbJ0fzkBqDMuOFzvcPbwztHRhCy8NCSTPLSw0lM/kmz+fB3dcGN3nKgk/x9joj7j3J65atQouuqh+bfS+++4DAIAJEybAtGnThAmGuEeHpg1gV10GQJkHR3WjHuvWC6oM6c1LUkIIulNSzSvpl5fNJ49FWfyO6x2u5qHSdAi3tkdKkgSZFnIHNWmQDF9NHMp8fHudpG6siEjyZaqsSc62hXsvMd4tE+/QwqZriaXlYW7F4sILL/TtvmgvadwgGXrkGg9mftzyBKDY+SBo+BQX4ddOSfRreery7vDU/zbBTSaplvlLDj52HgdrHo9Y6DdDIQn6cCqmMsp76Na1snTXsTITjndi5SkFLleIkxh1rC2yUuHbu843PP/13/UTLJG7SBLb4CJi4DXrDK0qt+d2yoHmmcY7NJRImv+dIMpr36fwikh7RKlJ6NblNlrnzeev7GU6CTIvUzyszesP57Qxjf+gKlfy52SA930acVZzz5aTeEHFAom8dFLdSqsfsG89ob+1XZpnwPJHRtgq2cvsuaJxKx6ITN+8bFuh53k7483PjLZcl+MwXottnVNjQv9Vn2jL6b2XdDZNgqWSyaZIdnjuCuvOs7FMmybp0KZJbGRhRsUCiegS8iBjtJ/dPR8L/07iWzdKg1WP2VNOlPhDlTOH11xOS3EuSZLjoeeV9zMtWX+rKAv/vE7stl/lPXRz+cGsLt7Q8Q1cSB+AxC5xa5e8Z4SYXQOBQNPn2PE/EIWZB7Tbbj7qhFoS5DRkX07xNRz3kVfRs+u06AfsKiZ6rHx0hGEeCSV223qsKK5OIUn+zJcStB1nSuJYseD3Uk5NTIBLewUv+qU2bLUIL3QR+EEMN2SI+HH44YJdxs+X/Mbv+kIvBm99KzQ1icQqkliK2Og13987zGsRAgHaszhITgzBP8YHLyKiVpf3QzI9O9tNYw3/zaXcg2aVsno/RE9Kx/VtJbZAj6h1YFkxqIpKl+bsTqGIPj4YQhC/YXfmLKLPiZGNEsLw6lLdUmr0rk/7jIP+zCUJbDmvWq43btR0OvGsvHsBKhZIFAm6vbuLnRMl1bPqa5Oe4kxVtRAx3FiekNd/Hxt7luN1RdfNfqwT96JpRgpc0r0+kuVnq/bDmr2FwuvxE7LzKuvttJ8rJL6Jd6XKC1CxQFSvXWpSCDJSvV8hYwlDbHh+jPWmkgQwpldL1+vl2TnrRPec0zAF3r1+oNAy/eio5yVO+FjE0pKB5J9d9HGD9yMI4hsIAKx/ahQk+cDJgoSzkCEO0ywzVejWWS/p2ybbaxEcQcSuENGKRRB2/CDOgYoFolLmjZSKIR2bMM0GRUwOhnTMgY5NGxrUYVxLLK3Vj+vbyvBanYZ166yf72m/NtnQv00jr8Xgwq3bObJ7czhTWeNSbQiCioXjXBegVLg3uxjfIt8kNXW2zVTSfuK1a/pGfr9+SHDaixI3lZJ43LZrhN93t9xyXnv41+ICR+twO8JsvOO9zTvgPHtFT69FCCRmwYVi1WHrmXH+bS/27ql7zyNoPhaoJyGxBlosEOFgPxhM7GU3ZTvuzxd3grIKfbN9rCqMMqgkRBMsNRABQMUCURC0mZ5I8Na4w70mEXGDZNJmXbIJettz4/qCfg/9BioWSMxxYdempsfgzFA8bmz/DQnIHItjiLP4OmMsBfS5cR/0sUBijmk3DvJaBIQbZ4f7WBk6Yn0pB0B8YrYgWaGQMKhYIDjFQ5iw00y654pJ5kUbmJtlpGATRnTBJGzug0shSGS6h52zPtgx2bMKdGraEA4VnrEtA212G/R2G/S2h/4PwQMtFohw42zzzFRo2yRd9VlGirs6rOh1Vez8/EssBWm10izjpe05pUCF20estJBggBYLRPiM75bz28Mt56uDaa33IKMj4i/iZHxE/AjqFa6CigUivMf3gxe29xIEDzfiWMQTCSEpMHla/AwBgv2By+BSCIKjMBIsYsg0wpqnBbGOH5eSbhvWwWsRHAUVCySiV/jxBUTEs/U5a3EI/Ns+fCtYFKjDRyMH5nOqfRHiDyuqkkcuPctrERwFFQskhrpldkT3I0Haa5+SKDYOAQuiBg0mJzx/jSGICUTzv/jycSnEbdDHAgmkZoFe4PGFakYawPYcDziVUqBN43Ro4vGS0wc3xVdQP1QsEOyHEcdx2hIdK2qk30zyfsKppZAXf9PbmYI5uKCLeRqCIIGKBYLJxxhAC4h/FdDwGrrXUiB2cWq5EZU590EfC8S3A4YdsC+JH2oJUSl+QfKHiQfkeU1trbdyIOJAxQLxsbc/gmjQUxh9pEg+Nlbf499HYvoO7IaCAyoWSCAR3YHjLNi+Amq05Pbk5d1tlS2pfvd2+L7l/GDHKHAKXJINDqhYIDhoIo5jNtT3zcu2XDa23thG7n9QrwgOlhSLt956C9q1awepqakwePBgWLFihWi5EBfBtU0kliEEICShj0WsU4uaRWDgViw+/fRTuO++++DJJ5+En3/+Gfr06QOjRo2Co0ePOiEfglgCnTfF49cBmxASM887VuT0An+2LsQK3IrFq6++CrfeeivceOON0L17d3j77bchPT0d3n//fSfkQ1wgmGub2IMHEdrATHQ+R2KDyK6QQPZD8QlXHIvKykpYvXo1TJo0KfJZKBSCESNGwNKlS6nnVFRUQEVFReTv4uJii6IiTnBl/1bQulE6/G3edvjd2XleiyOM9jkN4Fd9coWV17V5hrCyYpWOTRsyH9uthfp+DWzXCNJT9Lub5pmpcFFX4yBCvVtnQXKCei40pmcLaNM4HZplpgIAwNjeLaFN43RmObUM7dTE8rlmaO8JC1cNaM1132ORwR2awIaDxXBB5xyvRUEEIRGO6erBgwehVatWsGTJEhgyZEjk87/85S+waNEiWL58edQ5Tz31FDz99NNRnxcVFUFmZqZFsREEQRAEcZPi4mLIysoyHb8d3xUyadIkKCoqivzs27fP6SoRBEEQBPEIrqWQnJwcSEhIgCNHjqg+P3LkCLRo0YJ6TkpKCqSkeJsABkEQBEEQd+CyWCQnJ8OAAQNg3rx5kc9qa2th3rx5qqURBEEQBEHiE+4kZPfddx9MmDABBg4cCIMGDYLXX38dysrK4MYbb3RCPgRBEARBYghuxeKaa66BY8eOwRNPPAGHDx+Gvn37wqxZs6B58+ZOyIcgCIIgSAzBtStEBKxepQiCIAiC+Aff7ApBEARBECR+QMUCQRAEQRBhoGKBIAiCIIgwULFAEARBEEQYqFggCIIgCCIMVCwQBEEQBBEGKhYIgiAIgggDFQsEQRAEQYTBHXnTLnI8ruLiYrerRhAEQRDEIvK4bRZX03XFoqSkBAAA8vLy3K4aQRAEQRCblJSUQFZWlu73rof0rq2thYMHD0JGRgZIkiSs3OLiYsjLy4N9+/bFXajweL32eL1uALz2eLz2eL1uALx2v1w7IQRKSkogNzcXQiF9TwrXLRahUAhat27tWPmZmZme33yviNdrj9frBsBrj8drj9frBsBr98O1G1kqZNB5E0EQBEEQYaBigSAIgiCIMAKjWKSkpMCTTz4JKSkpXoviOvF67fF63QB47fF47fF63QB47bF27a47byIIgiAIElwCY7FAEARBEMR7ULFAEARBEEQYqFggCIIgCCIMVCwQBEEQBBFGYBSLt956C9q1awepqakwePBgWLFihdciMZOfnw9nn302ZGRkQLNmzeCKK66ArVu3qo658MILQZIk1c/tt9+uOmbv3r0wduxYSE9Ph2bNmsGDDz4I1dXVqmMWLlwI/fv3h5SUFOjUqRNMmzbN6csz5Kmnnoq6rm7dukW+Ly8vh4kTJ0KTJk2gYcOG8Jvf/AaOHDmiKiMWrxsAoF27dlHXLkkSTJw4EQCC88x/+OEHuPzyyyE3NxckSYIvv/xS9T0hBJ544glo2bIlpKWlwYgRI2D79u2qY06ePAnjx4+HzMxMyM7OhptvvhlKS0tVx/zyyy9w/vnnQ2pqKuTl5cFLL70UJcvnn38O3bp1g9TUVOjVqxfMnDlT+PUqMbr2qqoqeOihh6BXr17QoEEDyM3Nheuvvx4OHjyoKoPWTiZPnqw6JtauHQDghhtuiLqu0aNHq46Jxedudt20d16SJHj55Zcjx8TqM49AAsD06dNJcnIyef/998nGjRvJrbfeSrKzs8mRI0e8Fo2JUaNGkalTp5INGzaQtWvXkksvvZS0adOGlJaWRo654IILyK233koOHToU+SkqKop8X11dTXr27ElGjBhB1qxZQ2bOnElycnLIpEmTIsfs2rWLpKenk/vuu49s2rSJvPnmmyQhIYHMmjXL1etV8uSTT5IePXqoruvYsWOR72+//XaSl5dH5s2bR1atWkXOOecccu6550a+j9XrJoSQo0ePqq57zpw5BADIggULCCHBeeYzZ84kjz76KPniiy8IAJAZM2aovp88eTLJysoiX375JVm3bh351a9+Rdq3b0/OnDkTOWb06NGkT58+ZNmyZeTHH38knTp1Itdee23k+6KiItK8eXMyfvx4smHDBvLJJ5+QtLQ08s4770SO+emnn0hCQgJ56aWXyKZNm8hjjz1GkpKSyPr16z259sLCQjJixAjy6aefki1btpClS5eSQYMGkQEDBqjKaNu2LXnmmWdU7UDZN8TitRNCyIQJE8jo0aNV13Xy5EnVMbH43M2uW3m9hw4dIu+//z6RJIns3LkzckysPnOZQCgWgwYNIhMnToz8XVNTQ3Jzc0l+fr6HUlnn6NGjBADIokWLIp9dcMEF5O6779Y9Z+bMmSQUCpHDhw9HPpsyZQrJzMwkFRUVhBBC/vKXv5AePXqozrvmmmvIqFGjxF4AB08++STp06cP9bvCwkKSlJREPv/888hnmzdvJgBAli5dSgiJ3eumcffdd5OOHTuS2tpaQkgwn7m2o62trSUtWrQgL7/8cuSzwsJCkpKSQj755BNCCCGbNm0iAEBWrlwZOea7774jkiSRAwcOEEII+cc//kEaNWoUuW5CCHnooYdI165dI39fffXVZOzYsSp5Bg8eTP74xz8KvUY9aIOMlhUrVhAAIHv27Il81rZtW/Laa6/pnhOr1z5hwgQybtw43XOC8NxZnvm4cePIxRdfrPos1p95zC+FVFZWwurVq2HEiBGRz0KhEIwYMQKWLl3qoWTWKSoqAgCAxo0bqz7/6KOPICcnB3r27AmTJk2C06dPR75bunQp9OrVC5o3bx75bNSoUVBcXAwbN26MHKO8T/IxXt+n7du3Q25uLnTo0AHGjx8Pe/fuBQCA1atXQ1VVlUrmbt26QZs2bSIyx/J1K6msrIQPP/wQbrrpJlVyvqA+c5mCggI4fPiwSsasrCwYPHiw6hlnZ2fDwIEDI8eMGDECQqEQLF++PHLMsGHDIDk5OXLMqFGjYOvWrXDq1KnIMX6+FwDhd1+SJMjOzlZ9PnnyZGjSpAn069cPXn75ZdVyVyxf+8KFC6FZs2bQtWtXuOOOO+DEiROR7+LhuR85cgS+/fZbuPnmm6O+i+Vn7noSMtEcP34campqVJ0rAEDz5s1hy5YtHkllndraWrjnnntg6NCh0LNnz8jnv//976Ft27aQm5sLv/zyCzz00EOwdetW+OKLLwAA4PDhw9R7IH9ndExxcTGcOXMG0tLSnLw0KoMHD4Zp06ZB165d4dChQ/D000/D+eefDxs2bIDDhw9DcnJyVCfbvHlz02uSvzM6xsvr1vLll19CYWEh3HDDDZHPgvrMlchy0mRUXkOzZs1U3ycmJkLjxo1Vx7Rv3z6qDPm7Ro0a6d4LuQyvKS8vh4ceegiuvfZaVbKpu+66C/r37w+NGzeGJUuWwKRJk+DQoUPw6quvAkDsXvvo0aPhyiuvhPbt28POnTvhkUcegTFjxsDSpUshISEhLp77Bx98ABkZGXDllVeqPo/1Zx7zikXQmDhxImzYsAEWL16s+vy2226L/N6rVy9o2bIlDB8+HHbu3AkdO3Z0W0xhjBkzJvJ77969YfDgwdC2bVv47LPPPB/03OS9996DMWPGQG5ubuSzoD5zJJqqqiq4+uqrgRACU6ZMUX133333RX7v3bs3JCcnwx//+EfIz8+PqTDPWn73u99Ffu/Vqxf07t0bOnbsCAsXLoThw4d7KJl7vP/++zB+/HhITU1VfR7rzzzml0JycnIgISEhaqfAkSNHoEWLFh5JZY0777wTvvnmG1iwYIFpavnBgwcDAMCOHTsAAKBFixbUeyB/Z3RMZmambwbx7Oxs6NKlC+zYsQNatGgBlZWVUFhYqDpG+WyDcN179uyBuXPnwi233GJ4XBCfuSyn0fvbokULOHr0qOr76upqOHnypJB24HU/ISsVe/bsgTlz5pimxh48eDBUV1fD7t27ASC2r11Jhw4dICcnR9W+g/zcf/zxR9i6davpew8Qe8885hWL5ORkGDBgAMybNy/yWW1tLcybNw+GDBnioWTsEELgzjvvhBkzZsD8+fOjTFw01q5dCwAALVu2BACAIUOGwPr161UvotxJde/ePXKM8j7Jx/jpPpWWlsLOnTuhZcuWMGDAAEhKSlLJvHXrVti7d29E5iBc99SpU6FZs2YwduxYw+OC+Mzbt28PLVq0UMlYXFwMy5cvVz3jwsJCWL16deSY+fPnQ21tbUTZGjJkCPzwww9QVVUVOWbOnDnQtWtXaNSoUeQYv90LWanYvn07zJ07F5o0aWJ6ztq1ayEUCkWWCWL12rXs378fTpw4oWrfQX3uAGEr5YABA6BPnz6mx8bcM3fcPdQFpk+fTlJSUsi0adPIpk2byG233Uays7NV3vJ+5o477iBZWVlk4cKFqu1Fp0+fJoQQsmPHDvLMM8+QVatWkYKCAvLVV1+RDh06kGHDhkXKkLcejhw5kqxdu5bMmjWLNG3alLr18MEHHySbN28mb731lufbLu+//36ycOFCUlBQQH766ScyYsQIkpOTQ44ePUoICW83bdOmDZk/fz5ZtWoVGTJkCBkyZEjk/Fi9bpmamhrSpk0b8tBDD6k+D9IzLykpIWvWrCFr1qwhAEBeffVVsmbNmsjOh8mTJ5Ps7Gzy1VdfkV9++YWMGzeOut20X79+ZPny5WTx4sWkc+fOqm2HhYWFpHnz5uS6664jGzZsINOnTyfp6elR2+8SExPJK6+8QjZv3kyefPJJx7ffGV17ZWUl+dWvfkVat25N1q5dq3r3ZW//JUuWkNdee42sXbuW7Ny5k3z44YekadOm5Prrr4/pay8pKSEPPPAAWbp0KSkoKCBz584l/fv3J507dybl5eWRMmLxuZu1d0LC20XT09PJlClTos6P5WcuEwjFghBC3nzzTdKmTRuSnJxMBg0aRJYtW+a1SMwAAPVn6tSphBBC9u7dS4YNG0YaN25MUlJSSKdOnciDDz6oimlACCG7d+8mY8aMIWlpaSQnJ4fcf//9pKqqSnXMggULSN++fUlycjLp0KFDpA6vuOaaa0jLli1JcnIyadWqFbnmmmvIjh07It+fOXOG/OlPfyKNGjUi6enp5Ne//jU5dOiQqoxYvG6Z2bNnEwAgW7duVX0epGe+YMECavueMGECISS85fTxxx8nzZs3JykpKWT48OFR9+PEiRPk2muvJQ0bNiSZmZnkxhtvJCUlJapj1q1bR8477zySkpJCWrVqRSZPnhwly2effUa6dOlCkpOTSY8ePci3337r2HUTYnztBQUFuu++HMtk9erVZPDgwSQrK4ukpqaSs846i7zwwguqwTcWr/306dNk5MiRpGnTpiQpKYm0bduW3HrrrVGTwVh87mbtnRBC3nnnHZKWlkYKCwujzo/lZy6DadMRBEEQBBFGzPtYIAiCIAjiH1CxQBAEQRBEGKhYIAiCIAgiDFQsEARBEAQRBioWCIIgCIIIAxULBEEQBEGEgYoFgiAIgiDCQMUCQRAEQRBhoGKBIAiCIIgwULFAEARBEEQYqFggCIIgCCIMVCwQBEEQBBHG/wNeTm57cX1dAwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the \"chain\"\n",
    "plt.plot(parameters[burnin:iterations, 0], linewidth=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([6.700e+01, 5.200e+01, 7.800e+01, 1.470e+02, 2.390e+02, 4.730e+02,\n",
       "        7.770e+02, 1.102e+03, 1.459e+03, 1.821e+03, 1.988e+03, 2.125e+03,\n",
       "        1.871e+03, 1.666e+03, 1.374e+03, 9.720e+02, 6.640e+02, 4.980e+02,\n",
       "        2.720e+02, 1.810e+02, 7.400e+01, 4.700e+01, 1.700e+01, 1.700e+01,\n",
       "        7.000e+00, 8.000e+00, 0.000e+00, 0.000e+00, 2.000e+00, 2.000e+00]),\n",
       " array([0.01      , 0.27040193, 0.53080385, 0.79120578, 1.05160771,\n",
       "        1.31200964, 1.57241156, 1.83281349, 2.09321542, 2.35361735,\n",
       "        2.61401927, 2.8744212 , 3.13482313, 3.39522506, 3.65562698,\n",
       "        3.91602891, 4.17643084, 4.43683277, 4.69723469, 4.95763662,\n",
       "        5.21803855, 5.47844047, 5.7388424 , 5.99924433, 6.25964626,\n",
       "        6.52004818, 6.78045011, 7.04085204, 7.30125397, 7.56165589,\n",
       "        7.82205782]),\n",
       " <BarContainer object of 30 artists>)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGdCAYAAAAMm0nCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAArwklEQVR4nO3df3RU9Z3/8dcQmOFXMhAgmcwaQsQVCL9E0BgRCoVNDCmtR2rld5RUXDugEHVDWgsBtybCFrFKoWwFuqewoLuAFSoSQEiVIBh2GkBJhYLRhQlWIANxDZDM948e7rdTQEmYcfIJz8c595zcz+dz731/RM2LO597xxYIBAICAAAwSItIFwAAANBQBBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHFaRrqAcKmvr9fx48cVHR0tm80W6XIAAMA1CAQCOnv2rNxut1q0uPp9lmYbYI4fP67ExMRIlwEAABrhk08+0U033XTV/mYbYKKjoyX99R9ATExMhKsBAADXwu/3KzEx0fo9fjXNNsBc+tgoJiaGAAMAgGG+bvkHi3gBAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjNMy0gUACI1uszY1+thjRVkhrAQAwo87MAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxGhRgCgsLdccddyg6OlpxcXG67777VFFRETTmyy+/lMfjUadOndS+fXuNGTNGVVVVQWMqKyuVlZWltm3bKi4uTk8//bQuXrwYNGbHjh26/fbb5XA4dMstt2jlypWNmyEAAGh2GhRgdu7cKY/Ho927d6u4uFgXLlxQenq6ampqrDEzZ87UG2+8oddee007d+7U8ePHdf/991v9dXV1ysrK0vnz57Vr1y795je/0cqVKzV79mxrzNGjR5WVlaXhw4fL6/VqxowZ+uEPf6i33norBFMGAACmswUCgUBjD/7ss88UFxennTt3aujQoaqurlaXLl20evVqff/735ckHTp0SL169VJpaanuuusuvfnmm/rOd76j48ePKz4+XpK0dOlS5eXl6bPPPpPdbldeXp42bdqkAwcOWNcaO3aszpw5o82bN19TbX6/X06nU9XV1YqJiWnsFAFjdJu1qdHHHivKCmElANB41/r7+7rWwFRXV0uSYmNjJUllZWW6cOGCRo4caY3p2bOnunbtqtLSUklSaWmp+vbta4UXScrIyJDf79fBgwetMX97jktjLp3jSmpra+X3+4M2AADQPDU6wNTX12vGjBkaPHiw+vTpI0ny+Xyy2+3q0KFD0Nj4+Hj5fD5rzN+Gl0v9l/q+aozf79f//d//XbGewsJCOZ1Oa0tMTGzs1AAAQBPX6ADj8Xh04MABrVmzJpT1NFp+fr6qq6ut7ZNPPol0SQAAIExaNuagadOmaePGjSopKdFNN91ktbtcLp0/f15nzpwJugtTVVUll8tljdmzZ0/Q+S49pfS3Y/7+yaWqqirFxMSoTZs2V6zJ4XDI4XA0ZjpAk3E961gA4EbSoDswgUBA06ZN0/r167V9+3YlJycH9Q8cOFCtWrXStm3brLaKigpVVlYqLS1NkpSWlqb9+/fr5MmT1pji4mLFxMQoJSXFGvO357g05tI5AADAja1Bd2A8Ho9Wr16t119/XdHR0daaFafTqTZt2sjpdConJ0e5ubmKjY1VTEyMpk+frrS0NN11112SpPT0dKWkpGjSpEmaP3++fD6fnnnmGXk8HusOyj//8z/r5Zdf1r/8y79oypQp2r59u1599VVt2sTfTgEAQAPvwCxZskTV1dUaNmyYEhISrG3t2rXWmBdeeEHf+c53NGbMGA0dOlQul0vr1q2z+qOiorRx40ZFRUUpLS1NEydO1OTJkzVv3jxrTHJysjZt2qTi4mL1799fP//5z/XrX/9aGRkZIZgyAAAw3XW9B6Yp4z0wMFGk1sDwHhgATcU38h4YAACASCDAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADBOy0gXACDyus3a1OhjjxVlhbASALg23IEBAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIzDi+yAELuel8IBAK4Nd2AAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIzT4ABTUlKi0aNHy+12y2azacOGDUH9NpvtituCBQusMd26dbusv6ioKOg85eXlGjJkiFq3bq3ExETNnz+/cTMEAADNToMDTE1Njfr376/Fixdfsf/EiRNB2/Lly2Wz2TRmzJigcfPmzQsaN336dKvP7/crPT1dSUlJKisr04IFC1RQUKBly5Y1tFwAANAMNfhFdpmZmcrMzLxqv8vlCtp//fXXNXz4cN18881B7dHR0ZeNvWTVqlU6f/68li9fLrvdrt69e8vr9WrhwoWaOnVqQ0sGAADNTFjXwFRVVWnTpk3Kycm5rK+oqEidOnXSgAEDtGDBAl28eNHqKy0t1dChQ2W32622jIwMVVRU6PTp01e8Vm1trfx+f9AGAACap7B+lcBvfvMbRUdH6/777w9qf/zxx3X77bcrNjZWu3btUn5+vk6cOKGFCxdKknw+n5KTk4OOiY+Pt/o6dux42bUKCws1d+7cMM0EAAA0JWENMMuXL9eECRPUunXroPbc3Fzr5379+slut+vRRx9VYWGhHA5Ho66Vn58fdF6/36/ExMTGFQ4AAJq0sAWYP/zhD6qoqNDatWu/dmxqaqouXryoY8eOqUePHnK5XKqqqgoac2n/autmHA5Ho8MPAAAwS9jWwLzyyisaOHCg+vfv/7VjvV6vWrRoobi4OElSWlqaSkpKdOHCBWtMcXGxevToccWPjwAAwI2lwQHm3Llz8nq98nq9kqSjR4/K6/WqsrLSGuP3+/Xaa6/phz/84WXHl5aWatGiRfrjH/+oP//5z1q1apVmzpypiRMnWuFk/PjxstvtysnJ0cGDB7V27Vq9+OKLQR8RAQCAG1eDP0J6//33NXz4cGv/UqjIzs7WypUrJUlr1qxRIBDQuHHjLjve4XBozZo1KigoUG1trZKTkzVz5sygcOJ0OrVlyxZ5PB4NHDhQnTt31uzZs3mEGgAASJJsgUAgEOkiwsHv98vpdKq6uloxMTGRLgc3kG6zNkW6hG/UsaKsSJcAoBm51t/ffBcSAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4Yf0yRwDN3/W894Z3yABoLO7AAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjNPgAFNSUqLRo0fL7XbLZrNpw4YNQf0PPfSQbDZb0HbvvfcGjTl16pQmTJigmJgYdejQQTk5OTp37lzQmPLycg0ZMkStW7dWYmKi5s+f3/DZAQCAZqnBAaampkb9+/fX4sWLrzrm3nvv1YkTJ6ztP//zP4P6J0yYoIMHD6q4uFgbN25USUmJpk6davX7/X6lp6crKSlJZWVlWrBggQoKCrRs2bKGlgsAAJqhlg09IDMzU5mZmV85xuFwyOVyXbHvww8/1ObNm7V3714NGjRIkvTSSy9p1KhR+rd/+ze53W6tWrVK58+f1/Lly2W329W7d295vV4tXLgwKOgAAIAbU1jWwOzYsUNxcXHq0aOHHnvsMX3++edWX2lpqTp06GCFF0kaOXKkWrRooffee88aM3ToUNntdmtMRkaGKioqdPr06Stes7a2Vn6/P2gDAADNU8gDzL333qv/+I//0LZt2/T8889r586dyszMVF1dnSTJ5/MpLi4u6JiWLVsqNjZWPp/PGhMfHx805tL+pTF/r7CwUE6n09oSExNDPTUAANBENPgjpK8zduxY6+e+ffuqX79+6t69u3bs2KERI0aE+nKW/Px85ebmWvt+v58QAwBAMxX2x6hvvvlmde7cWYcPH5YkuVwunTx5MmjMxYsXderUKWvdjMvlUlVVVdCYS/tXW1vjcDgUExMTtAEAgOYp7AHm008/1eeff66EhARJUlpams6cOaOysjJrzPbt21VfX6/U1FRrTElJiS5cuGCNKS4uVo8ePdSxY8dwlwwAAJq4BgeYc+fOyev1yuv1SpKOHj0qr9eryspKnTt3Tk8//bR2796tY8eOadu2bfre976nW265RRkZGZKkXr166d5779UjjzyiPXv26N1339W0adM0duxYud1uSdL48eNlt9uVk5OjgwcPau3atXrxxReDPiICAAA3rgYHmPfff18DBgzQgAEDJEm5ubkaMGCAZs+eraioKJWXl+u73/2ubr31VuXk5GjgwIH6wx/+IIfDYZ1j1apV6tmzp0aMGKFRo0bpnnvuCXrHi9Pp1JYtW3T06FENHDhQTz75pGbPns0j1AAAQJJkCwQCgUgXEQ5+v19Op1PV1dWsh8E3qtusTZEuwRjHirIiXQKAJuZaf3+H/CkkoDkghABA08aXOQIAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDt9GDSBirudbv48VZYWwEgCm4Q4MAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgNDjAlJSUaPXq03G63bDabNmzYYPVduHBBeXl56tu3r9q1aye3263Jkyfr+PHjQefo1q2bbDZb0FZUVBQ0pry8XEOGDFHr1q2VmJio+fPnN26GAACg2WlwgKmpqVH//v21ePHiy/q++OIL7du3Tz/96U+1b98+rVu3ThUVFfrud7972dh58+bpxIkT1jZ9+nSrz+/3Kz09XUlJSSorK9OCBQtUUFCgZcuWNbRcAADQDLVs6AGZmZnKzMy8Yp/T6VRxcXFQ28svv6w777xTlZWV6tq1q9UeHR0tl8t1xfOsWrVK58+f1/Lly2W329W7d295vV4tXLhQU6dObWjJAACgmQn7Gpjq6mrZbDZ16NAhqL2oqEidOnXSgAEDtGDBAl28eNHqKy0t1dChQ2W32622jIwMVVRU6PTp01e8Tm1trfx+f9AGAACapwbfgWmIL7/8Unl5eRo3bpxiYmKs9scff1y33367YmNjtWvXLuXn5+vEiRNauHChJMnn8yk5OTnoXPHx8VZfx44dL7tWYWGh5s6dG8bZAACApiJsAebChQv6wQ9+oEAgoCVLlgT15ebmWj/369dPdrtdjz76qAoLC+VwOBp1vfz8/KDz+v1+JSYmNq54AADQpIUlwFwKLx9//LG2b98edPflSlJTU3Xx4kUdO3ZMPXr0kMvlUlVVVdCYS/tXWzfjcDgaHX4AAIBZQr4G5lJ4+eijj7R161Z16tTpa4/xer1q0aKF4uLiJElpaWkqKSnRhQsXrDHFxcXq0aPHFT8+AgAAN5YG34E5d+6cDh8+bO0fPXpUXq9XsbGxSkhI0Pe//33t27dPGzduVF1dnXw+nyQpNjZWdrtdpaWleu+99zR8+HBFR0ertLRUM2fO1MSJE61wMn78eM2dO1c5OTnKy8vTgQMH9OKLL+qFF14I0bQBAIDJbIFAINCQA3bs2KHhw4df1p6dna2CgoLLFt9e8vbbb2vYsGHat2+ffvSjH+nQoUOqra1VcnKyJk2apNzc3KCPgMrLy+XxeLR371517txZ06dPV15e3jXX6ff75XQ6VV1d/bUfYQF/r9usTZEuAV/jWFFWpEsAEAbX+vu7wQHGFAQYXA8CTNNHgAGap2v9/c13IQEAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxWka6ACBcus3aFOkSAABhwh0YAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABinwQGmpKREo0ePltvtls1m04YNG4L6A4GAZs+erYSEBLVp00YjR47URx99FDTm1KlTmjBhgmJiYtShQwfl5OTo3LlzQWPKy8s1ZMgQtW7dWomJiZo/f37DZwcAAJqlBgeYmpoa9e/fX4sXL75i//z58/WLX/xCS5cu1Xvvvad27dopIyNDX375pTVmwoQJOnjwoIqLi7Vx40aVlJRo6tSpVr/f71d6erqSkpJUVlamBQsWqKCgQMuWLWvEFAEAQHNjCwQCgUYfbLNp/fr1uu+++yT99e6L2+3Wk08+qaeeekqSVF1drfj4eK1cuVJjx47Vhx9+qJSUFO3du1eDBg2SJG3evFmjRo3Sp59+KrfbrSVLlugnP/mJfD6f7Ha7JGnWrFnasGGDDh06dE21+f1+OZ1OVVdXKyYmprFThMH4Nurm7VhRVqRLABAG1/r7O6RrYI4ePSqfz6eRI0dabU6nU6mpqSotLZUklZaWqkOHDlZ4kaSRI0eqRYsWeu+996wxQ4cOtcKLJGVkZKiiokKnT58OZckAAMBALUN5Mp/PJ0mKj48Pao+Pj7f6fD6f4uLigoto2VKxsbFBY5KTky87x6W+jh07Xnbt2tpa1dbWWvt+v/86ZwMAAJqqZvMUUmFhoZxOp7UlJiZGuiQAABAmIQ0wLpdLklRVVRXUXlVVZfW5XC6dPHkyqP/ixYs6depU0JgrneNvr/H38vPzVV1dbW2ffPLJ9U8IAAA0SSENMMnJyXK5XNq2bZvV5vf79d577yktLU2SlJaWpjNnzqisrMwas337dtXX1ys1NdUaU1JSogsXLlhjiouL1aNHjyt+fCRJDodDMTExQRsAAGieGhxgzp07J6/XK6/XK+mvC3e9Xq8qKytls9k0Y8YM/eu//qt+97vfaf/+/Zo8ebLcbrf1pFKvXr1077336pFHHtGePXv07rvvatq0aRo7dqzcbrckafz48bLb7crJydHBgwe1du1avfjii8rNzQ3ZxAEAgLka/Bj1jh07NHz48Mvas7OztXLlSgUCAc2ZM0fLli3TmTNndM899+iXv/ylbr31VmvsqVOnNG3aNL3xxhtq0aKFxowZo1/84hdq3769Naa8vFwej0d79+5V586dNX36dOXl5V1znTxGDR6jxtXwCDbQdF3r7+/reg9MU0aAAQEGV0OAAZquiLwHBgAA4JtAgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAME7LUJ+wW7du+vjjjy9r/9GPfqTFixdr2LBh2rlzZ1Dfo48+qqVLl1r7lZWVeuyxx/T222+rffv2ys7OVmFhoVq2DHm5aOK6zdoU6RIAAE1QyBPB3r17VVdXZ+0fOHBA//RP/6QHHnjAanvkkUc0b948a79t27bWz3V1dcrKypLL5dKuXbt04sQJTZ48Wa1atdJzzz0X6nIBAICBQh5gunTpErRfVFSk7t2761vf+pbV1rZtW7lcrisev2XLFn3wwQfaunWr4uPjddttt+nZZ59VXl6eCgoKZLfbQ10yAAAwTFjXwJw/f16//e1vNWXKFNlsNqt91apV6ty5s/r06aP8/Hx98cUXVl9paan69u2r+Ph4qy0jI0N+v18HDx686rVqa2vl9/uDNgAA0DyFdVHJhg0bdObMGT300ENW2/jx45WUlCS3263y8nLl5eWpoqJC69atkyT5fL6g8CLJ2vf5fFe9VmFhoebOnRv6SQAAgCYnrAHmlVdeUWZmptxut9U2depU6+e+ffsqISFBI0aM0JEjR9S9e/dGXys/P1+5ubnWvt/vV2JiYqPPBwAAmq6wBZiPP/5YW7dute6sXE1qaqok6fDhw+revbtcLpf27NkTNKaqqkqSrrpuRpIcDoccDsd1Vg0AAEwQtjUwK1asUFxcnLKysr5ynNfrlSQlJCRIktLS0rR//36dPHnSGlNcXKyYmBilpKSEq1wAAGCQsNyBqa+v14oVK5SdnR307pYjR45o9erVGjVqlDp16qTy8nLNnDlTQ4cOVb9+/SRJ6enpSklJ0aRJkzR//nz5fD4988wz8ng83GEBAACSwhRgtm7dqsrKSk2ZMiWo3W63a+vWrVq0aJFqamqUmJioMWPG6JlnnrHGREVFaePGjXrssceUlpamdu3aKTs7O+i9MQAA4MYWlgCTnp6uQCBwWXtiYuJlb+G9kqSkJP3+978PR2kAAKAZ4LuQAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGCeuXOQJAU9Rt1qZGH3us6Ku/HgXAN4M7MAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYJeYApKCiQzWYL2nr27Gn1f/nll/J4POrUqZPat2+vMWPGqKqqKugclZWVysrKUtu2bRUXF6enn35aFy9eDHWpAADAUC3DcdLevXtr69at//8iLf//ZWbOnKlNmzbptddek9Pp1LRp03T//ffr3XfflSTV1dUpKytLLpdLu3bt0okTJzR58mS1atVKzz33XDjKBQAAhglLgGnZsqVcLtdl7dXV1XrllVe0evVqffvb35YkrVixQr169dLu3bt11113acuWLfrggw+0detWxcfH67bbbtOzzz6rvLw8FRQUyG63h6NkAABgkLCsgfnoo4/kdrt18803a8KECaqsrJQklZWV6cKFCxo5cqQ1tmfPnuratatKS0slSaWlperbt6/i4+OtMRkZGfL7/Tp48OBVr1lbWyu/3x+0AQCA5inkASY1NVUrV67U5s2btWTJEh09elRDhgzR2bNn5fP5ZLfb1aFDh6Bj4uPj5fP5JEk+ny8ovFzqv9R3NYWFhXI6ndaWmJgY2okBAIAmI+QfIWVmZlo/9+vXT6mpqUpKStKrr76qNm3ahPpylvz8fOXm5lr7fr+fEAMg5LrN2tToY48VZYWwEuDGFvbHqDt06KBbb71Vhw8flsvl0vnz53XmzJmgMVVVVdaaGZfLddlTSZf2r7Su5hKHw6GYmJigDQAANE9hDzDnzp3TkSNHlJCQoIEDB6pVq1batm2b1V9RUaHKykqlpaVJktLS0rR//36dPHnSGlNcXKyYmBilpKSEu1wAAGCAkH+E9NRTT2n06NFKSkrS8ePHNWfOHEVFRWncuHFyOp3KyclRbm6uYmNjFRMTo+nTpystLU133XWXJCk9PV0pKSmaNGmS5s+fL5/Pp2eeeUYej0cOhyPU5QIAAAOFPMB8+umnGjdunD7//HN16dJF99xzj3bv3q0uXbpIkl544QW1aNFCY8aMUW1trTIyMvTLX/7SOj4qKkobN27UY489prS0NLVr107Z2dmaN29eqEvFN+R61gwAAHAltkAgEIh0EeHg9/vldDpVXV3NepgII8AAf8UiXuDrXevvb74LCQAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgnJahPmFhYaHWrVunQ4cOqU2bNrr77rv1/PPPq0ePHtaYYcOGaefOnUHHPfroo1q6dKm1X1lZqccee0xvv/222rdvr+zsbBUWFqply5CXDADfiG6zNjX62GNFWSGsBDBfyNPAzp075fF4dMcdd+jixYv68Y9/rPT0dH3wwQdq166dNe6RRx7RvHnzrP22bdtaP9fV1SkrK0sul0u7du3SiRMnNHnyZLVq1UrPPfdcqEsGAACGCXmA2bx5c9D+ypUrFRcXp7KyMg0dOtRqb9u2rVwu1xXPsWXLFn3wwQfaunWr4uPjddttt+nZZ59VXl6eCgoKZLfbQ102AAAwSNjXwFRXV0uSYmNjg9pXrVqlzp07q0+fPsrPz9cXX3xh9ZWWlqpv376Kj4+32jIyMuT3+3Xw4MErXqe2tlZ+vz9oAwAAzVNYF5TU19drxowZGjx4sPr06WO1jx8/XklJSXK73SovL1deXp4qKiq0bt06SZLP5wsKL5KsfZ/Pd8VrFRYWau7cuWGaCQAAaErCGmA8Ho8OHDigd955J6h96tSp1s99+/ZVQkKCRowYoSNHjqh79+6NulZ+fr5yc3Otfb/fr8TExMYVDgAAmrSwBZhp06Zp48aNKikp0U033fSVY1NTUyVJhw8fVvfu3eVyubRnz56gMVVVVZJ01XUzDodDDocjBJXjSq7n6QkAAEIt5GtgAoGApk2bpvXr12v79u1KTk7+2mO8Xq8kKSEhQZKUlpam/fv36+TJk9aY4uJixcTEKCUlJdQlAwAAw4T8DozH49Hq1av1+uuvKzo62lqz4nQ61aZNGx05ckSrV6/WqFGj1KlTJ5WXl2vmzJkaOnSo+vXrJ0lKT09XSkqKJk2apPnz58vn8+mZZ56Rx+PhLgsAAAj9HZglS5aourpaw4YNU0JCgrWtXbtWkmS327V161alp6erZ8+eevLJJzVmzBi98cYb1jmioqK0ceNGRUVFKS0tTRMnTtTkyZOD3hsDAABuXCG/AxMIBL6yPzEx8bK38F5JUlKSfv/734eqLAAA0IzwXUgAAMA4BBgAAGAcAgwAADAOAQYAABgnrG/iBQCExvW8TPJYUVYIKwGaBu7AAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDk8h3SCu5wkGAACaGu7AAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBx+CoBg/B1AAAA/BUBBgCauev5y8+xoqwQVgKEDgGmEfifAQAAkcUaGAAAYBwCDAAAMA4fIX3DWIgLwCR8ZI6mijswAADAOAQYAABgHAIMAAAwTpMOMIsXL1a3bt3UunVrpaamas+ePZEuCQAANAFNdhHv2rVrlZubq6VLlyo1NVWLFi1SRkaGKioqFBcXF+nyAABfI1IPLbB4+MbQZO/ALFy4UI888ogefvhhpaSkaOnSpWrbtq2WL18e6dIAAECENck7MOfPn1dZWZny8/OtthYtWmjkyJEqLS294jG1tbWqra219qurqyVJfr8/5PXV134R8nMCAEKj68zXInLdA3MzInLd5ubS7+1AIPCV45pkgPnLX/6iuro6xcfHB7XHx8fr0KFDVzymsLBQc+fOvaw9MTExLDUCAPC3nIsiXUHzcvbsWTmdzqv2N8kA0xj5+fnKzc219uvr63Xq1Cl16tRJNpstZNfx+/1KTEzUJ598opiYmJCdt6m5EeZ5I8xRYp7NyY0wR4l5NjcNnWcgENDZs2fldru/clyTDDCdO3dWVFSUqqqqgtqrqqrkcrmueIzD4ZDD4Qhq69ChQ7hKVExMTLP+F+6SG2GeN8IcJebZnNwIc5SYZ3PTkHl+1Z2XS5rkIl673a6BAwdq27ZtVlt9fb22bdumtLS0CFYGAACagiZ5B0aScnNzlZ2drUGDBunOO+/UokWLVFNTo4cffjjSpQEAgAhrsgHmwQcf1GeffabZs2fL5/Pptttu0+bNmy9b2PtNczgcmjNnzmUfVzU3N8I8b4Q5SsyzObkR5igxz+YmXPO0Bb7uOSUAAIAmpkmugQEAAPgqBBgAAGAcAgwAADAOAQYAABiHANMAixcvVrdu3dS6dWulpqZqz549kS4p5EpKSjR69Gi53W7ZbDZt2LAh0iWFXGFhoe644w5FR0crLi5O9913nyoqKiJdVsgtWbJE/fr1s14elZaWpjfffDPSZYVVUVGRbDabZsyYEelSQqqgoEA2my1o69mzZ6TLCov//d//1cSJE9WpUye1adNGffv21fvvvx/pskKqW7dul/152mw2eTyeSJcWMnV1dfrpT3+q5ORktWnTRt27d9ezzz77td9v1BAEmGu0du1a5ebmas6cOdq3b5/69++vjIwMnTx5MtKlhVRNTY369++vxYsXR7qUsNm5c6c8Ho92796t4uJiXbhwQenp6aqpqYl0aSF10003qaioSGVlZXr//ff17W9/W9/73vd08ODBSJcWFnv37tWvfvUr9evXL9KlhEXv3r114sQJa3vnnXciXVLInT59WoMHD1arVq305ptv6oMPPtDPf/5zdezYMdKlhdTevXuD/iyLi4slSQ888ECEKwud559/XkuWLNHLL7+sDz/8UM8//7zmz5+vl156KXQXCeCa3HnnnQGPx2Pt19XVBdxud6CwsDCCVYWXpMD69esjXUbYnTx5MiApsHPnzkiXEnYdO3YM/PrXv450GSF39uzZwD/+4z8GiouLA9/61rcCTzzxRKRLCqk5c+YE+vfvH+kywi4vLy9wzz33RLqMb9wTTzwR6N69e6C+vj7SpYRMVlZWYMqUKUFt999/f2DChAkhuwZ3YK7B+fPnVVZWppEjR1ptLVq00MiRI1VaWhrByhAK1dXVkqTY2NgIVxI+dXV1WrNmjWpqaprl13F4PB5lZWUF/Tfa3Hz00Udyu926+eabNWHCBFVWVka6pJD73e9+p0GDBumBBx5QXFycBgwYoH//93+PdFlhdf78ef32t7/VlClTQvrFw5F29913a9u2bfrTn/4kSfrjH/+od955R5mZmSG7RpN9E29T8pe//EV1dXWXvQU4Pj5ehw4dilBVCIX6+nrNmDFDgwcPVp8+fSJdTsjt379faWlp+vLLL9W+fXutX79eKSkpkS4rpNasWaN9+/Zp7969kS4lbFJTU7Vy5Ur16NFDJ06c0Ny5czVkyBAdOHBA0dHRkS4vZP785z9ryZIlys3N1Y9//GPt3btXjz/+uOx2u7KzsyNdXlhs2LBBZ86c0UMPPRTpUkJq1qxZ8vv96tmzp6KiolRXV6ef/exnmjBhQsiuQYDBDc3j8ejAgQPNcj2BJPXo0UNer1fV1dX6r//6L2VnZ2vnzp3NJsR88skneuKJJ1RcXKzWrVtHupyw+du/tfbr10+pqalKSkrSq6++qpycnAhWFlr19fUaNGiQnnvuOUnSgAEDdODAAS1durTZBphXXnlFmZmZcrvdkS4lpF599VWtWrVKq1evVu/eveX1ejVjxgy53e6Q/VkSYK5B586dFRUVpaqqqqD2qqoquVyuCFWF6zVt2jRt3LhRJSUluummmyJdTljY7XbdcsstkqSBAwdq7969evHFF/WrX/0qwpWFRllZmU6ePKnbb7/daqurq1NJSYlefvll1dbWKioqKoIVhkeHDh1066236vDhw5EuJaQSEhIuC9e9evXSf//3f0eoovD6+OOPtXXrVq1bty7SpYTc008/rVmzZmns2LGSpL59++rjjz9WYWFhyAIMa2Cugd1u18CBA7Vt2zarrb6+Xtu2bWuW6wmau0AgoGnTpmn9+vXavn27kpOTI13SN6a+vl61tbWRLiNkRowYof3798vr9VrboEGDNGHCBHm93mYZXiTp3LlzOnLkiBISEiJdSkgNHjz4slca/OlPf1JSUlKEKgqvFStWKC4uTllZWZEuJeS++OILtWgRHDGioqJUX18fsmtwB+Ya5ebmKjs7W4MGDdKdd96pRYsWqaamRg8//HCkSwupc+fOBf2t7ujRo/J6vYqNjVXXrl0jWFnoeDwerV69Wq+//rqio6Pl8/kkSU6nU23atIlwdaGTn5+vzMxMde3aVWfPntXq1au1Y8cOvfXWW5EuLWSio6MvW7vUrl07derUqVmtaXrqqac0evRoJSUl6fjx45ozZ46ioqI0bty4SJcWUjNnztTdd9+t5557Tj/4wQ+0Z88eLVu2TMuWLYt0aSFXX1+vFStWKDs7Wy1bNr9fxaNHj9bPfvYzde3aVb1799b//M//aOHChZoyZUroLhKy55luAC+99FKga9euAbvdHrjzzjsDu3fvjnRJIff2228HJF22ZWdnR7q0kLnS/CQFVqxYEenSQmrKlCmBpKSkgN1uD3Tp0iUwYsSIwJYtWyJdVtg1x8eoH3zwwUBCQkLAbrcH/uEf/iHw4IMPBg4fPhzpssLijTfeCPTp0yfgcDgCPXv2DCxbtizSJYXFW2+9FZAUqKioiHQpYeH3+wNPPPFEoGvXroHWrVsHbr755sBPfvKTQG1tbciuYQsEQvhaPAAAgG8Aa2AAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMM7/A5hIKxRGqaRBAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot a histogram of the Beta posterior\n",
    "plt.hist(parameters[burnin:iterations, 0], bins=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.986173127210402\n"
     ]
    }
   ],
   "source": [
    "# Calculate and print the mean of the Beta posterior\n",
    "print(np.mean(parameters[burnin:iterations, 0]))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mcmc",
   "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.9.1"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "9276ee63f0f305db9b395444cd35bfd67dfb6213f72130623ae4c18412569e6b"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}