{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6d46ce96",
   "metadata": {},
   "source": [
    "### Exercise 63"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3927052b",
   "metadata": {},
   "source": [
    "#### Import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1f9caf5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import poisson, gamma\n",
    "from scipy.stats import gaussian_kde as kde"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab4b218a",
   "metadata": {},
   "source": [
    "#### Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9edf321e",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('mining.csv')\n",
    "xy = df.to_numpy()\n",
    "x, y = xy[:, 0], xy[:, 1]\n",
    "n = len(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f70b7711",
   "metadata": {},
   "source": [
    "#### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a21e25b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAAd90lEQVR4nO3de5QjZ3nn8e9P3TMYYYPtVicLMS2ZDRAuwdgtDM46YCAJ1z2BAyQkDhB7c3p3shvbsLCBsNkAORNCWC4BdoE5Bg6JBJwNl8A6BMPJwTHhlvT4bnMz9jTYgRjMLdjLbfrZP6ra1oxVUpVUpVZX/z7n1Bmp6n3fep6qmmc0r6SSIgIzM6ufxnYHYGZm1XCBNzOrKRd4M7OacoE3M6spF3gzs5pa3O4ABrVareh0OtsdhpnZjnHw4MFvRsTysG1zVeA7nQ7r6+vbHYaZ2Y4haSNrm6dozMxqygXezKymXODNzGrKBd7MrKZc4M3MaqrSAi/peEnvlfR5SZ+TdEaV+6tCv9+n0+nQaDTodDr0+/1tGcPMrKiqPyb558BHIuKZkvYCzYr3V6p+v8/a2hq33347ABsbG6ytrQFw9tlnz2wMM7NJqKrbBUu6F3AFcL/IuZNutxvz9Dn4TqfDxsZdP2Labrc5dOjQzMYwM8si6WBEdIdtq3KK5mTgG8A7JF0u6UJJ9xgS3JqkdUnr3/jGNyoMp7ivfOUrhdZXNYaZ2SSqLPCLwGnAmyPiVOA24MVHN4qIAxHRjYju8vLQb9tum5WVlULrqxrDzGwSVRb4m4CbIuKz6fP3khT8HWP//v00m0e+bdBsNtm/f/9MxzAzm0RlBT4ivg58VdID01WPB66ran9VOPvsszlw4ADtdhtJtNttDhw4UOjN0TLGMDObRGVvsgJIejhwIbAXuAE4JyK+ndV+3t5kNTObd6PeZK30Y5IRcQUwdMdmZlYtf5PVzKymXODNzGrKBd7MrKZc4M3MasoF3sysplzgzcxqygXezKymXODNzGrKBd7MrKZc4M3MasoF3sysplzgzcxqygXezKymXODNzGrKBd7MrKZc4M3MasoF3sysplzgzcxqygXezKymXODNzGrKBd7MrKZc4M3MasoF3sysplzgzcxqqtICL+mQpKslXSFpvcp9AfT7fTqdDo1Gg06nQ7/fr3qXpciKe9p8Bvu3Wi1arVausXbqcTSzo0REZQtwCGjlbb+6uhqT6vV60Ww2A7hjaTab0ev1Jh5zFrLi3rdv31T5DBs3z1g79Tia7VbAemTUVCXbqyHpENCNiG/mad/tdmN9fbIX+p1Oh42Njbusb7fbHDp0aKIxZyEr7oWFBQ4fPnyX9XnzyRp33Fg79Tia7VaSDkZEd+i2igv8jcC3SV4JvjUiDgxpswasAaysrKyOK0pZGo0Gw3KRxObm5kRjzkJW3Fny5pNn3GFj7dTjaLZbjSrwVb/JemZEnAY8CfjPkh59dIOIOBAR3YjoLi8vT7yjlZWVQuvnRVZ8CwsLhdpP0m5Ym516HM3sriot8BFxc/rnLcAHgNOr2tf+/ftpNptHrGs2m+zfv7+qXZYiK+61tbWp8hk2bp6xdupxNLMhsibnp12AewDHDTz+FPDEUX2meZM1InmDsN1uh6Rot9s75o3BrLinzWew/9LSUiwtLeUaa6ceR7PdiO14k1XS/UhetQMsAu+KiJEvA6d5k9XMbDcaNQe/WNVOI+IG4JSqxjczs9H8TVYzs5pygTczqykXeDOzmnKBNzOrKRd4M7OacoE3M6spF3gzs5pygTczqykXeDOzmnKBNzOrKRd4M7OacoE3M6upsQVe0r+VdLf08VmSzpN0fOWRmZnZVPK8gn8fcFjSzwIHgPsC76o0KjMzm1qeAr8ZET8Bng68MSJeBNy72rDMzGxaeQr8jyX9BvA84KJ03Z7qQjIzszLkKfDnAGcA+yPiRkknA39ZbVhmZjatkb/oJGkBeGlEnL21LiJuBF5VdWBmZjadka/gI+Iw0Ja0d0bxmJlZSfL8JusNwCclfQi4bWtlRLy2sqjMzGxqeQr8l9OlARxXbThmZlaWsQU+Il4OIKkZEbdXH5KZmZUhzzdZz5B0HfD59Pkpkv535ZGZmdlU8nxM8vXAE4BbASLiSuDRFcZkZmYlyHWzsYj46lGrDufdgaQFSZdLumh86+L6/T6dTodGo0Gn06Hf7w9dV1a/vDGUmU+etq1Wi1arVSiGMuIuM3ezQdNeW742gYgYuQDvBX4BuIzkG6wvBN4zrt9A/xeQ3LvmonFtV1dXo4herxfNZjOAO5Y9e/bE3r17j1jXbDaj1+tN3S9vDHn6TTvWsLZFYigj7jJzNxs07bW1m65NYD2y6m/WhrizQLeAPvAvwC1ADzhxXL+070nA3wGPq6LAt9vtzAJ39NJut6fuVySGcf2mHStPDqNiKCPuMnM3GzTttbWbrs1RBV7J9myS/l1EfHLcuoy+7wVeSfLxyhdGxFOHtFkD1gBWVlZWNzY2xg17h0ajwbj4B/bD5ubmVP2KxDCu37Rj5clhVAxlxF1m7maDpr22dtO1KelgRHSHbcszB//GnOuO3ulTgVsi4uCodhFxICK6EdFdXl7OEc6dVlZWJmo7ab8i24vsY5Kx8ow/qk0ZcZeZu9mgaa8tX5uprJf2JDcY+6/AV0nm0beWlwFXZvUb6P9K4CbgEPB14HagN6qP5+A9B28W4Tn4IphkDh54DPBHwNfSP7eWFwD3z+qXMdZZVDAHH5GcyHa7HZKi3W5Hr9cbuq6sfnljmFSRsQbbLi0txdLSUqEYyoi7zNzNBk17be2Wa3NUgc8zB9+OiI30cQM4NiK+N+Y/BkePcRYZc/CDut1urK+vFxnazGxXm3YO/pWS7inpHsA1wHWSXlQkgIi4ZFxxNzOzcuUp8A9OX7E/Dfhb4GTgOVUGZWZm08tT4PdI2kNS4D8UET8medPCzMzmWJ4C/1aST8LcA7hUUhsoNAdvZmazl+d2wW8A3jCwakPSY6sLyczMypDnBz+Q9BTgIcAxA6tfUUlEZmZWijz3g38L8OvA7wECngW0K47LzMymlGcO/hci4rnAtyP5daczgAdUG5aZmU0rT4H/f+mft0u6D/Bj4N7VhWRmZmXIMwd/kaTjgVeT3BM+gAurDMrMzKaX51M0f5w+fF/6q0zHRMR3qw3LzMymledN1mdJOi59+iLgHZJOrTYsMzObVp45+D+MiH+VdCbwS8DbgLdUG5aZmU0rT4Hf+oHtpwAHIuJvgL3VhWRmZmXIU+BvlvRWks/Cf1jS3XL2MzOzbZSnUP8acDHwhIj4DnAiyVy8mZnNscxP0Ui6Z3qb4GOAS9J1JwI/BPyrHGZmc27UxyTfBTwVOEjy2XcNbAvgfhXGZWZmU8os8Fu/wBQRJ88uHDMzK8uoKZrTRnWMiMvKD8fMzMoyaormNemfxwBd4EqSaZqHkczBn1FtaGZmNo3MT9FExGMj4rHA14DTIqIbEavAqcDNswrQzMwmk+djkg+MiKu3nkTENcCDqgvJzMzKkOdukldJuhDopc/PBq6qLiQzMytDngJ/DrAPOD99finw5soiMjOzUoydoomIH0TE6yLi6enyuoj4wbh+ko6R9I+SrpR0raSXlxNyufr9Pp1Oh0ajQafTod/v52rbarVotVq5+g0bQxKLi4tIKrTfIvuqarxJj9m0sZd9HCbZxyximFdlXP82YxFRyULyiZtj08d7gM8CjxrVZ3V1NWap1+tFs9kMki9uBRDNZjN6vV6utnn65R2jyH7z7KuMnKftX2bsZR+HSfYxixjmVRnXv1UDWI+sOpy1ocwFaJL8GtQjR7WbdYFvt9tDL9Z2u5277bh+RcYost9x+yoj52n7lxl72cdhkn3MIoZ5Vcb1b9UYVeCVbL8rSX8ZEc+RdH5E/PnQRmNIWiC51cHPAv8rIn5/SJs1YA1gZWVldWNjY5JdTaTRaDAsf0lsbm7majuuX5Exiux33L6KxpB3vDKO2SSxl30cJtnHLGKYV2Vc/1YNSQcjojts26g5+NX0R7bPlXSCpBMHlzw7jojDEfFw4CTgdEkPHdLmQCSfse8uLy/nGbY0KysruddntS3SZpLtRWLMY9rxyjhmk8Re9nGYZB+ziGFelXH92zbIemkPnAd8juTukTcANw4sN2T1GzHe/wBeOKqN5+A9B19V3GXsw3PwnoOfR0wzBw+8eVybjH7LwPHp47sDnwCeOqrPrAt8RHLhttvtkBTtdnvkRTrYdmlpKZaWlnL1GzYGEAsLC3fMXebdb5F9VTXepMds2tjLPg6T7GMWMcyrMq5/K9+oAp85Bz9I0inAL6ZPL42IsV90kvQw4J3AAslU0P+JiFeM6tPtdmN93beaNzPLa9Qc/NgvOkk6j+RN0Penq/qSDkTEG0f1S/8ROLVosGZmVo4832T9HZKPN94GIOlVwKeBkQXezMy2V56bjQk4PPD8MEf+upOZmc2hPK/g3wF8VtIH0udPA95WWURmZlaKsQU+Il4r6RLgzHTVORFxeaVRmZnZ1PK8gieSn+fzT/SZme0geebgzcxsB3KBNzOrqZEFXtKCpI/PKhgzMyvPyAIfEYeBTUn3mlE8ZmZWkjxvsn4fuFrSx4DbtlZGxHmVRWVmZlPLU+Dfz523KTAzsx0iz+fg3ynp7sBKRHxhBjGZmVkJxn6KRtK/B64APpI+f7ikD1Ucl5mZTSnPxyRfBpwOfAcgIq4A7ldZRGZmVoo8Bf7HEfHdo9b5hxfNzOZcnjdZr5X0m8CCpPuT/JTfp6oNy8zMppXnFfzvAQ8h+W3WdwPfAy6oMCYzMytBnk/R3A68NP2hj4iIf60+LDMzm1aeT9E8QtLVwFUkX3i6UtJq9aGZmdk08szBvw343Yj4BICkM0l+BORhVQZmZmbTyTMHf3iruANExD8AP6kuJDMzK0PmK3hJp6UP/17SW0neYA3g14FLqg/NzMymMWqK5jVHPf+jgcdRQSxmZlaizAIfEY+dZSBmZlauPJ+iOV7SeZJeK+kNW0uOfveV9HFJ10m6VtL55YS8M/T7fTqdDo1Gg06nQ7/f3+6QRsqKd2u9JBYXF5GUmc+4nIuMtRMM5ttqtWi1Wtt2vqu63iYdt8x+O+3v0lyJiJELybdWXwucAzxva8nR797Aaenj44AvAg8e1Wd1dTXqoNfrRbPZDJKprACi2WxGr9fb7tCGyop33759d1mflc+4nIdt3wnHJsuofGadU1XX26Tjltlvz549sXfv3h1/vVQJWI+sOpy1Ie4s1JeNa5NnAT4I/PKoNnUp8O12e+hf+na7vd2hDZUV78LCQmYBOzqfcTlnbZ/3Y5NlXD6zzKmq623SccvuV4frpUqjCryS7dkkPZ/kV50uIrldAQAR8a2RHY8cowNcCjw0Ir531LY1YA1gZWVldWNjI++wc6vRaDDsuEpic3P+7tOWFe84g/mMy3ncPub12GTJc8xmlVNV19uk45bdb5iddr1USdLBiOgO25bnc/A/Al4NfBo4mC7rBXZ+LPA+4IKjiztARByIiG5EdJeXl/MOO9dWVlYKrd9uWXEtLCzk7jcu53G5z+uxyZIn3lnlVNX1Num4Zfebtu2ulvXSfmsBbgBa49pl9N0DXAy8IE/7ukzReA7ec/Ceg/cc/Kww5Rz8R4HmuHZD+gn4C+D1efvUpcBHJBdru90OSdFut+f+gsyKd2s9A3PyWfmMy7nIWDvBYL5LS0uxtLS0bee7qutt0nHL7LfT/i7N2qgCn2cO/gMktwv+OEfOwZ83pt+ZwCeAq7nzB0L+ICI+nNWn2+3G+nru2R8zs11v1Bx8npuN/XW6FBLJPWtUtJ+ZmZUjz/3g3zmLQMzMrFxjC7ykG+Gu956JCP/wtpnZHMszRTM4t3MM8CzgxGrCMTOzsoz9HHxE3Dqw3BwRrweeUn1oZmY2jTxTNKcNPG2QvKLP88rfzMy2UZ5CPXhf+J8Ah4BfqyQaMzMrTZ5P0fi+8GZmO1CeKZq7Ac8AOoPtI+IV1YVlZmbTyjNF80HguyQ3GfvhmLZmZjYn8hT4kyLiiZVHYmZmpcpzu+BPSfr5yiMxM7NS5XkFfybw2+k3Wn9Icn+ZiIiHVRqZmZlNJU+Bf1LlUZiZWenyfExy5/+GnpnZLpRnDt7MzHYgF3gzs5pygTczqykXeDOzmnKBNzOrKRd4M7OacoE3M6spF3gzs5pygTczqykXeDOzmqqswEt6u6RbJF1T1T7y6Pf7dDodGo0GnU6Hfr+/neFMbVw+85BvkRgG27ZaLVqt1hGPJbG4uIikQmMNts3ax6TjTXIchuU2OG6e2PPGUMbxL5Lv1hhZ52rYPgbbTrrfaXMvo1+R/sPaVP73NSIqWYBHA6cB1+Tts7q6GmXq9XrRbDYDuGNpNpvR6/VK3c+sjMtnHvItEsOwtuOWImM1m83Yt2/fyH0UHS/vsSyS2549e2Lv3r25Yh8XQ5nHP0++o8bIc/wn3e+0uZfRr0j/YW2yznvR3IH1yKrDWRvKWEh+5m/bCny73R56AbXb7VL3Myvj8pmHfIvEkNV23FJkrIWFhVLHy3ssJ80tT+yjYij7+I/Ld9wYeY5/WdfspOesqnM92L/I9VA0d0YUeCXbqyGpA1wUEQ8d0WYNWANYWVlZ3dgo7+aVjUaDYflJYnNzs7T9zMq4fOYh3yIxZLUdp8yxio6X91hOE884o2Io+/iPy7eqPCe5Zic9Z1Wd68H+RY5T0dwlHYyI7tDYco9SkYg4EBHdiOguLy+XOvbKykqh9fNuXD7zkG+RGCaNq8hYCwsLpY6XN+YyjnlW7KPGLvv4j2szbnue4z/JuEX6TJrDtOd6cH2RfEr9+5r10r6MhW2eopmHOekyeQ7ec/Cegy8Wi+fga1zgI5ID2263Q1K02+0dW9y3jMtnHvItEsNg26WlpVhaWjriMdw5j1tkrMG2WfuYdLxJjsOw3AbHzRN73hjKOP5F8t0aI+tcDdvHYNtJ9ztt7mX0K9J/WJsy/r6OKvCVzcFLejdwFtAC/gX4o4h426g+3W431tfXK4nHzKyORs3B5/lN1olExG9UNbaZmY237W+ymplZNVzgzcxqygXezKymXODNzGrKBd7MrKZc4M3MasoF3sysplzgzcxqygXezKymXODNzGrKBd7MrKZc4M3MasoF3sysplzgzcxqygXezKymXODNzGrKBd7MrKZc4M3MasoF3sysplzgzcxqygXezKymXODNzGrKBd7MrKZc4M3MaqrSAi/piZK+IOl6SS+ucl959Pt9Op0OjUaDTqdDv98vpe00fWwyg8e61WrRarWOOO7jzsW4/llty7xustoOW58Vb5HYt9pKYnFxEUmZMc7rtTyL47DVdnB71r7GXXuTXJuliohKFmAB+DJwP2AvcCXw4FF9VldXoyq9Xi+azWYAdyzNZjN6vd5UbafpY5MZdqwHlz179sTevXszz8W4/uPalnHdZLXdt2/fXdYPy2fS2Mf1KZrHLOXNZ5rjMOz4F1nGnatx1+YkgPXIqsNZG6ZdgDOAiweevwR4yag+VRb4drs99IC32+2p2k7TxyaTdazHLVvnIk//cW2nvW6y2i4sLExcXCbJc1iM83otFz3vkxyHMo7/NOdsEowo8Eq2l0/SM4EnRsTvpM+fAzwyIv7LUe3WgDWAlZWV1Y2NjUriaTQaDMtVEpubmxO3naaPTSbrWI+zdS7y9B/XdtrrZtIc8iiS59F9RsW23ddy0WM2yXHYLtMcW0kHI6I7bNu2v8kaEQciohsR3eXl5cr2s7Kyknt9kbbT9LHJTHpMt/rl6T+u7bTXTVbbhYWFsbGNUyTPYfHM67VcdP+THIcyjv8kKju2WS/tp12Ysykaz8HXh+fgsxfPwU93HDwHn7/ALwI3ACdz55usDxnVp8oCH5Gc6Ha7HZKi3W6PLdh5207TxyYzeKyXlpZiaWnpiOM+7lyM65/VtszrJqvtsPVZ8RaJfast3DnXnBXjvF7LszgOW20Ht2fta9y1N8m1WdSoAl/ZHDyApCcDryf5RM3bI2L/qPbdbjfW19cri8fMrG5GzcEvVrnjiPgw8OEq92FmZsNt+5usZmZWDRd4M7OacoE3M6spF3gzs5qq9FM0RUn6BlDNV1mn1wK+ud1BVMS57UzObWcqO7d2RAz9luhcFfh5Jmk966NIO51z25mc2840y9w8RWNmVlMu8GZmNeUCn9+B7Q6gQs5tZ3JuO9PMcvMcvJlZTfkVvJlZTbnAm5nV1K4u8JLeLukWSdcMrHu4pM9IukLSuqTTB7adla6/VtLfD6yfqx8Xh2K5SbqXpP8r6co0t3MG+jxP0pfS5XnbkcvRMnI7RdKnJV2d5nLPgW0vSc/NFyQ9YWD9XJ23InlJ+mVJB9P1ByU9bqDParr+eklvkKTtyGdQ0XOWbl+R9H1JLxxYN1fnDCa6Hh+Wbrs23X5Mur7885Z1H+HdsACPBk4DrhlY91HgSenjJwOXpI+PB64DVtLnP5X+WfjHxecwtz8AXpU+Xga+leZyIsk9/U8ETkgfnzCnuf0T8Jj08bnAH6ePH5yek7uR/DbBl9NzNnfnrWBepwL3SR8/FLh5oM8/Ao8CBPzt1jnfKbkNbH8v8FfAC9Pnc3fOJjhvi8BVwCnp8yVgoarztqtfwUfEpSTF7IjVwNa/tvcC/jl9/JvA+yPiK2nfW9L1pwPXR8QNEfEj4D3Ar1YaeA4FcwvguPQVw7Fpv58ATwA+FhHfiohvAx8Dnlh17ONk5PYA4NL08ceAZ6SPfxV4T0T8MCJuBK4nOWdzd96K5BURl0fE1vm7Fri7pLtJujdwz4j4TCRV4y+Ap1Ue/BgFzxmSngbcSJLblrk7Z1A4t18BroqIK9O+t0bE4arO264u8BkuAF4t6avA/yT5qUFITtgJki5J/0v83HT9zwBfHeh/U7puHl3A8NzeBDyIpOBfDZwfEZvsrNyu5c6/7M8C7ps+zsphp+SWldegZwCXRcQPSXK4aWDbvOYFGblJOhb4feDlR7XfKecMss/bA4CQdLGkyyT9t3R9JefNBf6u9gHPj4j7As8H3pauXwRWgaeQvLL9Q0kP2J4QJ5aV2xOAK4D7AA8H3nT0fOgOcC7wu5IOAscBP9rmeMoyMi9JDwFeBfzHbYhtWlm5vQx4XUR8f7sCK0FWbovAmcDZ6Z9Pl/T4qoKo9BeddqjnAeenj/8KuDB9fBNwa0TcBtwm6VLglHT94Kuqk4CbZxRrUVm5nQP8afpfw+sl3Qj8HEkeZw30Pwm4ZCaRFhQRnyf57y/pP7xPSTfdTPb5mfvzNiIvJJ0EfAB4bkR8OV19M0kuW+YyLxiZ2yOBZ0r6M5L3vjYl/QA4yA44ZzAyt5uASyPim+m2D5PM3/eo4Lz5Ffxd/TPwmPTx44AvpY8/CJwpaVFSk+Qi/BzJmyn3l3SypL3As4EPzTjmvLJy+wrweABJPw08kOQN1YuBX5F0gqQTSC7Yi2cacU6Sfir9swH8d+At6aYPAc9O56dPBu5P8mbWjjhvWXlJOh74G+DFEfHJrfYR8TXge5Ielb6n8lySa3fuZOUWEb8YEZ2I6JD8pvOfRMSb2CHnDEZejxcDPy+pKWmR5O/jdZWdt+1+B3o7F+DdwNeAH5P8y/ofSP7bdJDkHfrPAqsD7V9E8kmaa4ALBtY/GfgiyTv8L93uvIrmRjI181GS+fdrgN8aGOdckjcmrwfO2e68RuR2fnoOvgj8Kem3tNP2L03PzRcY+GTCvJ23InmRFI3bSKbWtpatT3Z10/P4ZZL3V7RdOU16zgb6vYz0UzTzeM4mvB5/i2SO/hrgzwbWl37efKsCM7Oa8hSNmVlNucCbmdWUC7yZWU25wJuZ1ZQLvJlZTbnAm5nVlAu8WYkkLWx3DGZbXOBt15L0CkkXDDzfL+l8SS+S9E+SrpL08oHtf53eaO5aSWsD678v6TWSrgTOmG0WZtlc4G03ezvJV8K3vlL+bODrJLczOJ3kxmurkh6dtj83IlZJvnF4nqSldP09gM9GxCkR8Q8zjN9sJN9szHatiDgk6VZJpwI/DVwOPILknjuXp82OJSn4l5IU9aen6++brr8VOAy8b5axm+XhAm+73YXAbwP/huQV/eOBV0bEWwcbSToL+CXgjIi4XdIlwDHp5h9ExOEZxWuWm6dobLf7AMmvVD2C5E5/FwPnpj86gaSfSe8MeC/g22lx/zmSn1Yzm2t+BW+7WkT8SNLHge+kr8I/KulBwKfT3zz+Psnd/z4C/CdJnyO5K+Vntitms7x8N0nb1dI3Vy8DnhURXxrX3mwn8RSN7VqSHkxyn/u/c3G3OvIreDOzmvIreDOzmnKBNzOrKRd4M7OacoE3M6spF3gzs5r6/zUDFrymJb/JAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x, y, color='black')\n",
    "plt.xlabel('year')\n",
    "plt.ylabel('number of disasters')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b3bc33a",
   "metadata": {},
   "source": [
    "#### (a)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89327a0f",
   "metadata": {},
   "source": [
    "#### Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a5eaffde",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha, beta = 1, 1\n",
    "sim = 1_000\n",
    "mplot = 200"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d071c60",
   "metadata": {},
   "source": [
    "#### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8e573ba7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def logsumexp(x):\n",
    "    xstar = np.max(x)\n",
    "    return xstar + np.log(np.sum(np.exp(x - xstar)))\n",
    "\n",
    "\n",
    "def update_tau(lambda1, lambda2, y, n):\n",
    "    logprobs = np.zeros(n)\n",
    "    for i in range(n):\n",
    "        y1 = y[:i]\n",
    "        y2 = y[i:]\n",
    "        \n",
    "        logprob1 = np.sum(poisson(lambda1).logpmf(y1))\n",
    "        logprob2 = np.sum(poisson(lambda2).logpmf(y2))\n",
    "        logprob = logprob1 + logprob2\n",
    "        logprobs[i] = logprob\n",
    "    logprobs -= logsumexp(logprobs)\n",
    "    probs = np.exp(logprobs)\n",
    "    \n",
    "    return np.random.choice(np.arange(n), p=probs)\n",
    "\n",
    "\n",
    "def update_lambda1(tau, lambda2, y, n, alpha, beta):\n",
    "    y1 = y[:tau]\n",
    "    return gamma(a=alpha + np.sum(y1), scale=1 / (beta + len(y1))).rvs()\n",
    "\n",
    "\n",
    "def update_lambda2(tau, lambda1, y, n, alpha, beta):\n",
    "    y2 = y[tau:]\n",
    "    return gamma(a=alpha + np.sum(y2), scale = 1 / (beta + len(y2))).rvs()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70773fca",
   "metadata": {},
   "source": [
    "#### Gibbs sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "27d56901",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = np.zeros((sim, 3))\n",
    "sample[0] = np.array([n // 2, np.mean(y), np.mean(y)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "072cea95",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "100\n",
      "200\n",
      "300\n",
      "400\n",
      "500\n",
      "600\n",
      "700\n",
      "800\n",
      "900\n"
     ]
    }
   ],
   "source": [
    "for s in range(sim - 1):\n",
    "    if s % 100 == 0:\n",
    "        print(s)\n",
    "    tau, lambda1, lambda2 = sample[s]\n",
    "    new_tau = update_tau(lambda1, lambda2, y, n)\n",
    "    new_lambda1 = update_lambda1(new_tau, lambda2, y, n, alpha, beta)\n",
    "    new_lambda2 = update_lambda2(new_tau, new_lambda1, y, n, alpha, beta)\n",
    "    sample[s + 1] = np.array([new_tau, new_lambda1, new_lambda2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dedeed9",
   "metadata": {},
   "source": [
    "#### Posterior"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4d50d4c",
   "metadata": {},
   "source": [
    "##### Tau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "980097f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "tau_posterior = np.zeros(n)\n",
    "for i in range(n):\n",
    "    tau_posterior[i] = np.sum(sample[:, 0] == i) / sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e3e80065",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeN0lEQVR4nO3dfXBc9X3v8fdXkiXLAT8L25EwfsAkgdaE2NhmbgsEEnBaMIHCxGlzY0Iy9CHMTdLAxbnkktQkNE7btPcO3BSmYVIgQBNz27oJqWFaCHOT2Njm2SYGY8u25EekRbHDSpbk7/3jnBUr6ay0ks/Z3SN/XjM7PnsetL/fHKEP3/M7+zvm7oiIiAxUVe4GiIhIZVJAiIhIJAWEiIhEUkCIiEgkBYSIiESqKXcD4jJ9+nSfM2dOuZshIpIqW7dufcvdG6K2jZmAmDNnDlu2bCl3M0REUsXM9hTapktMIiISSQEhIiKRFBAiIhJJASEiIpEUECIiEkkBISIikRQQIiISSQEhqbdt2zaeffbZcjdDZMwZM1+Uk1PXHXfcweuvv8727dvL3RSRMUUVhKTe3r17yWaz5W6GyJijCkJSr7W1lerq6nI3Q2TMUQUhqdbV1cXhw4c5fvx4uZsiMuYoICTV9u/fDwRBISLxUkBIqrW2tgKoghBJgAJCUq2lpQUIAsLdy9wakbFFASGplqsgQFWESNwUEJJquQoCFBAicVNASKrlB4QGqkXipYCQVNMlJpHkKCAk1VpaWvq+JKcKQiReCghJrd7eXvbv38/s2bMBVRAicUs0IMxsuZntMLOdZrY6Yvufm9l2M3vZzP7DzM7K27bKzN4IX6uSbKek0+HDh+nt7WXevHmAKgiRuCUWEGZWDdwLfAw4F/ikmZ07YLcXgMXuvhBYB3w7PHYq8DVgKbAE+JqZTUmqrZJOuQHqXECoghCJV5IVxBJgp7vvcvfjwGPANfk7uPvT7v5O+HYj0BQuXwk85e7t7p4BngKWJ9hWSaHcALUqCJFkJBkQjcC+vPct4bpCPgv8dJTHyikoV0HMnz8fUECIxK0ipvs2s08Bi4FLRnjczcDNQN9ApZw6WlpaGDduHI2Nwf876BKTSLySrCBagTPz3jeF6/oxs48AdwAr3L1rJMe6+/3uvtjdFzc0NMTWcEmH1tZWGhsbGT9+PKAKQiRuSQbEZmCBmc01s1pgJbA+fwczuwC4jyAcDudt2gBcYWZTwsHpK8J1In1aWlpobGykrq4OUAUhErfEAsLde4BbCP6wvwb80N23mdkaM1sR7vZXwGnAj8zsRTNbHx7bDtxFEDKbgTXhOpE+LS0tNDU1UVtbC6iCEIlbomMQ7v4E8MSAdXfmLX9kiGMfAB5IrnWSZu5Oa2srK1asUAUhkhB9k1pSKZPJkM1mVUGIJEgBIamUu8W1qampr4JQQIjESwEhqZT7klxjY2NfBaFLTCLxUkBIKuUHhCoIkWQoICSV2traAJg+fTo1NTVUVVWpghCJmQJCUimTyTBu3DgmTJgAQG1trSoIkZgpICSVMpkMU6ZMwcwAqKurUwUhEjMFhKRSe3s7U6dO7XuvCkIkfgoISaVcBZGjCkIkfgoISaWBAaEKQiR+CghJpagKQgEhEi8FhKRSVAWhS0wi8VJASOr09vbS0dHRb5BaFYRI/BQQkjodHR24uwapRRKmgJDUyWQyABqkFkmYAkJSJyogVEGIxE8BIamTCwh9UU4kWQoISZ1CFYQCQiReCghJnfb24PHkus1VJFkKCEkdVRAipaGAkNTJZDLU1dVRX1/ft04VhEj8FBCSOplMpt8ANaiCEEmCAkJSZ+A0G6DbXEWSoICQ1Glvbx8UELlLTO5eplaJjD0KCEmdQhUEQHd3dzmaJDImKSAkdaICora2FkDjECIxUkBI6hQapAYFhEicFBCSKj09Pfz6178uWEFooFokPgoISZW3334boOAYhCoIkfgoICRVor5FDaogRJKggJBUiZrJFVRBiCRBASGpUqiCyAWEKgiR+CggJFWGu8SkCkIkPgoISRVVECKlo4CQVIl6FgSoghBJggJCUiWTyTBhwoS+iiFHg9Qi8VNASKpETbMBus1VJAmJBoSZLTezHWa208xWR2y/2MyeN7MeM7t+wLZeM3sxfK1Psp2SHoUCQhWESPxqkvrBZlYN3At8FGgBNpvZenffnrfbXuBG4NaIH5F19w8m1T5Jp+ECQhWESHySrCCWADvdfZe7HwceA67J38Hdm939ZeBEgu2QMSTqWRCgQWqRJCQZEI3Avrz3LeG6Yo03sy1mttHMPh61g5ndHO6z5ciRIyfRVEmLqJlcQRWESBIqeZD6LHdfDPwh8HdmNn/gDu5+v7svdvfFDQ0NpW+hlNxwg9SqIETik2RAtAJn5r1vCtcVxd1bw393Ac8AF8TZOEmf7u5ujh07pkFqkRJJMiA2AwvMbK6Z1QIrgaLuRjKzKWZWFy5PB/4LsH3oo2SsKzTVN0B1dTVmpktMIjFKLCDcvQe4BdgAvAb80N23mdkaM1sBYGYXmlkLcANwn5ltCw//ALDFzF4Cnga+NeDuJzkF5QJi8uTJg7aZGXV1daogRGKU2G2uAO7+BPDEgHV35i1vJrj0NPC4XwC/nWTbJH06OzsBqK+vj9xeW1urCkIkRpU8SC3STy4gBk6zkaMKQiReCghJjdwf//Hjx0dur6urUwUhEiMFhKRGroIoFBC1tbWqIERipICQ1Mj98R/qEpMqCJH4KCAkNVRBiJSWAkJSo5gKQgEhEh8FhKRGMRWELjGJxEcBIamhCkKktBQQkhqqIERKSwEhqaEKQqS0FBCSGsV8k1oVhEh8FBCSGl1dXdTW1mJmkdt1m6tIvBQQkhqdnZ0Fxx9Al5hE4qaAkNTo6uoqeHkJNEgtEjcFhKSGKgiR0lJASGqoghApLQWEpEYxFcTx48dx9xK2SmTsUkBIahRTQQB0d3eXqkkiY5oCQlKjmAoC0DiESEwUEJIaw1UQuW0ahxCJhwJCUmO4CiJ3iUkVhEg8FBCSGsVWEAoIkXgoICQ1iq0gdIlJJB4KCEkNVRAipaWAkNRQBSFSWgoISQ1VECKlVTPcDmZ2Z9R6d18Tf3NEClMFIVJawwYE8Ju85fHAVcBryTRHJJq7q4IQKbFhA8Ld/yb/vZn9NbAhsRaJRMhVBfomtUjpjGYMYgLQFHdDRIaS+6OvS0wipVPMGMQrQG56zGqgAdD4g5TUcM+jzt+mCkIkHsWMQVyVt9wDHHL3noTaIxIpFxCqIERKp5gxiD0AZnYGwSD1e80Md9+bdONEcnJVgSoIkdIZdgzCzFaY2RvAbuBnQDPw04TbJdKPKgiR0itmkPouYBnwurvPBS4HNibaKpEBVEGIlF4xAdHt7m1AlZlVufvTwOKE2yXSjyoIkdIrJiDeNrPTgGeBH5jZ/wKOFfPDzWy5me0ws51mtjpi+8Vm9ryZ9ZjZ9QO2rTKzN8LXqmI+T8auYiqImpoazEwVhEhMirmL6SXgHeBLwB8Bk4DThjvIzKqBe4GPAi3AZjNb7+7b83bbC9wI3Drg2KnA1wgqFQe2hsdmimivjEHFVBBmRl1dnQJCJCbFBMSH3f0EcAL4RwAze7mI45YAO919V3jMY8A1QF9AuHtzuO3EgGOvBJ5y9/Zw+1PAcuDRIj5XxqBiKggILjPpEpNIPApeYjKzPw2/JPd+M3s577UbKCYgGoF9ee9bwnXFKOpYM7vZzLaY2ZYjR44U+aOlUh05coS1a9fS29s7aFsxFQRAfX0977zzTiLtEznVDDUG8QhwNfCv4b+51yJ3/1QJ2jYsd7/f3Re7++KGhoZyN0dO0o9//GNWr17Nc889N2hbsRXEpEmT6OjoSKR9IqeaggHh7h3u3uzun3T3PXmv9iJ/ditwZt77pnBd0sdKSmWzWQA2bdo0aFuxFcTkyZMVECIxSfKBQZuBBWY218xqgZXA+iKP3QBcYWZTzGwKcAWaQXbMy4XAxo2Dv2Yzkgri7bffjr1tIqeixAIinK/pFoI/7K8BP3T3bWa2xsxWAJjZhWbWAtwA3Gdm28Jj2wm+oLc5fK0ZQeUiKZULCFUQIpWhmLuYRs3dnwCeGLDuzrzlzRSYOtzdHwAeSLJ9UllyIdDc3MyhQ4eYMWNG37ZcBZH7MlwhqiBE4qNnUkvFyAUEDK4iOjs7qa2tpapq6F9ZVRAi8VFASMXIZrNMmDCBmpqaQQEx3ONGcyZNmkQ2m9V3IURioICQitHZ2cmkSZNYuHBhZAUx3PgDBBUEoCpCJAYKCKkYnZ2d1NfXs3TpUp577rl+X5gbSQUBaBxCJAYKCKkYuSph2bJlHD16lF/96leDtg1HFYRIfBQQUjFyIbB06VKg/0C1KgiR0lNASMXIBcSCBQuYPHlyv4BQBSFSegoIqRi5EKiqqmLJkiWqIETKTAEhFSObzfZVCfPnz6e19d3pt1RBiJSeAkIqRu4uJoCpU6fS3t7OiRPBo0KKrSBOP/10zEwVhEgMFBBSMfKrhKlTp3LixAmOHj06aNtQqqqqmDhxoioIkRgoIKRiDAwIgPb2YI7GYisI0HxMInFRQEjFGCogiq0gQPMxicRFASEVY7gKotiAUAUhEg8FhFQEd+93F1NUBVHsJSZVECLxUEBIRejp6eHEiRP97mKC0V1iUgUhEg8FhFSEgU+MmzJlChAEhLuPaJBaFYRIPBQQUhEGBkRdXR3vec97aGtr63u2w0gqiI6ODtw9mcaKnCIUEFIRop45nfuyXO5xoyOpIE6cOMGxY8fib6jIKUQBIRVhqICI2jYUzcckEg8FhFSEuCsI0HxMIidLASEVIZvNAqogRCqJAkIqQi4Ecre5gioIkXJTQEhF0BiESOVRQEhFKBQQ3d3dtLW1AaogREpNASEVoVBAAOzfv3/QtqGoghCJhwJCKsJQAXHgwAGg+Api/Pjx1NXVqYIQOUkKCKkIhe5igpFXEPDut6lFZPQUEFIRCt3FBCOvICAYh9AlJpGTo4CQihDnGASoghCJgwJCKkIuIPKrBFUQIuWlgJCK0NnZSW1tLVVV7/5K1tfXU1dX1xcQqiBESksBIRUh6oFAZsbUqVP7pvtWBSFSWgoIqQj5jxvNl7vMBCMLCFUQIidPASEVobOzs98dTDm5gBg3bly/y0/DmTx5Mu+88w7d3d2xtVHkVKOAkIpQ6JnT06ZNA0Y2/gDvfptaVYTI6CUaEGa23Mx2mNlOM1sdsb3OzP4p3L7JzOaE6+eYWdbMXgxff59kO6X8CgVEroIYyeUleHc+Jo1DiIxeTVI/2MyqgXuBjwItwGYzW+/u2/N2+yyQcfezzWwlsBb4RLjtTXf/YFLtk8oyXECoghApvSQriCXATnff5e7HgceAawbscw3wj+HyOuByM7ME2yQVShWESOVJMiAagX1571vCdZH7uHsP0AFMC7fNNbMXzOxnZva7CbZTKkDcFUQuIDKZzEm3TeRUVamD1AeA2e5+AfDnwCNmNnHgTmZ2s5ltMbMtR44cKXkjJT7D3eY60gqioaEBgMOHD59840ROUUkGRCtwZt77pnBd5D5mVgNMAtrcvcvd2wDcfSvwJnDOwA9w9/vdfbG7L879QZB0Gu4215FWEA0NDZgZhw4diqV9IqeiJANiM7DAzOaaWS2wElg/YJ/1wKpw+XrgP93dzawhHOTGzOYBC4BdCbZVyizuS0w1NTVMnz6dgwcPxtI+kVNRYncxuXuPmd0CbACqgQfcfZuZrQG2uPt64HvAQ2a2E2gnCBGAi4E1ZtYNnAD+xN3bk2qrlF/cg9QAM2fOVAUhchISCwgAd38CeGLAujvzljuBGyKOexx4PMm2SWWJu4IAmDFjhgJC5CRU6iC1nGIKBcRpp51GTU3NqCsIXWISGT0FhJRdb28vx48fjwyI3IyuUQPYw8lVEO4eRzNFTjmJXmISKUZXVxdAwRD47ne/y7x580b8c2fOnEk2m+Xo0aNMnDjoLmkRGYYCQsou6nGj+a677rpR/dwZM2YAcOjQIQWEyCjoEpOU3XABMVozZ84E0DiEyCgpIKTskgqI/ApCREZOASFll3RAqIIQGR0FhJRdUgExffp0qqqqVEGIjJICQsoum80Che9iGq3q6moaGhpUQYiMkgJCyi6pCgKGnm7jjjvu4Omnn479M0XGCt3mKmWXZEAUmm6jq6uLu+++myNHjvDhD3849s8VGQtUQUjZJV1BRF1i2rt3L6ABbJGhKCCk7EpRQQycbmP37t2AAkJkKAoIKbukK4iuri46Ojr6rW9ubgYUECJDUUBI2eXuYkqqgoDBX5bLDwhN5icSTQEhZZerIOK+zRUKT7eRu8TU3d1NJpOJ/XNFxgIFhJRd0mMQULiCADhw4EDsnysyFiggpOw6Ozuprq6mpib+u64LTbfR3NzM2WefHblNRAIKCCm7Qk+Ti8O0adOorq7uV0Fks1kOHjzIRRddBCggRApRQEjZJRkQVVVVnHHGGf0CYs+ePQAsW7YMUECIFKKAkLJLMiBg8JflcuMPCxcupL6+XgEhUoACQsoum80mcgdTzsDpNnIBMXfuXGbOnKlBapECFBBSdqWuIHbv3s24ceOYNWtWwak4REQBIRUg6YAYON1Gc3MzZ511FlVVVcyaNUsBIVKAAkLKLumAaGpqoru7m127dgFBQMyZMwcoPJmfiCggpAIkHRArVqzAzHj44YeBICDmzp0LBAHR1tbG8ePHE/t8kbRSQEjZJR0Qs2fP5vLLL+f73/8+x44d4/Dhw/0qCBj8TWsRUUBIBUj6LiaAG2+8kebmZh588EGAQQGhy0wigykgpOySriAArr32WiZOnMjdd98NKCBEiqGAkLIrRUBMmDCBlStX0traCtA3BjFr1ixAASESRQEhZVeKgIDgMhNAXV1d3yR+Z5xxBqCAEIkS//SZIiNUqoBYtmwZ73vf+4BgjiaA2tpapk2bpoAQiaCAkLJy95IFhJnxyCOPcOzYsX7rNd2GSDQFhJRVd3c37l6SgAD40Ic+NGidviwnEk1jEFJWuedRJ32b61A03YZINAWElFWSjxstVq6CyM3VJCIBBYQkav/+/axdu5a33npr0Lbe3l5uu+02AObNm1fqpvWZOXMm2WyWo0ePlq0NIpUo0YAws+VmtsPMdprZ6ojtdWb2T+H2TWY2J2/bV8L1O8zsyiTbKcloaWnhkksuYfXq1Zx33nmsW7eub1tPTw+rVq3ioYce4q677uLKK8t3inNfltNAtUh/iQWEmVUD9wIfA84FPmlm5w7Y7bNAxt3PBv4WWBseey6wEjgPWA78n/DnJWbHjh184xvf4LrrruOee+7RH4uTtG/fPi699FIOHz7Mgw8+SFNTEzfccAMXXnghl112GYsWLeIHP/gB3/zmN/nqV79a1rbmAuK+++7rm/E1Lu7Opk2buPXWW1m5ciUPP/wwHR0dsX6GSFIsqeuuZnYR8HV3vzJ8/xUAd//LvH02hPv80sxqgINAA7A6f9/8/Qp93uLFi33Lli0jbueePXu4+uqreeWVV4BgauiWlhbMjAULFlBToxu9RuPgwYP09PTw5JNPsnTpUrq7u/nOd77DT37yk759PvGJT/D5z3++jK0MZDIZrrrqKn7xi18AweWuuMZEMpkMBw4cYNy4cUydOpVDhw5RW1vL/PnzMbNYPkNk4cKFPProo6M61sy2uvviqG1J/vVrBPblvW8Blhbax917zKwDmBau3zjg2MaBH2BmNwM3QzBj56ga2djI7Nmz+dznPsd1111HU1MT27dvZ926dX2hISN3/vnn8+Uvf5lFixYBMG7cOG6//XZuv/32MrdssClTpvDzn/+cPXv28Pjjj7Nx48bYBqzr6uq44ooruPrqq5k0aRKbNm1i3bp17N27N5afLwLvTh0TtyQriOuB5e7+ufD9fwWWuvstefu8Gu7TEr5/kyBEvg5sdPeHw/XfA37q7usoYLQVhIjIqWyoCiLJQepW4My8903hush9wktMk4C2Io8VEZEEJRkQm4EFZjbXzGoJBp3XD9hnPbAqXL4e+E8PSpr1wMrwLqe5wALguQTbKiIiAyQ2BhGOKdwCbACqgQfcfZuZrQG2uPt64HvAQ2a2E2gnCBHC/X4IbAd6gM+7e29SbRURkcESG4MoNY1BiIiMXLnGIEREJMUUECIiEkkBISIikRQQIiISacwMUpvZEWBPudsxhOnA4ClNxwb1LZ3Ut3SKu29nuXtD1IYxExCVzsy2FLpTIO3Ut3RS39KplH3TJSYREYmkgBARkUgKiNK5v9wNSJD6lk7qWzqVrG8agxARkUiqIEREJJICQkREIikgToKZPWBmh8MHH+XWfdDMNprZi2a2xcyW5G27NFy/zcx+lrd+uZntMLOdZra61P2IMpK+mdkkM/s3M3sp7Ntn8o5ZZWZvhK9VUZ9VagX6dr6Z/dLMXgn7MjFv21fCc7PDzK7MW19R520k/TKzj5rZ1nD9VjO7LO+YReH6nWb2v60Cno060nMWbp9tZsfM7Na8dRV1zmBUv48Lw23bwu3jw/Xxnzd312uUL+Bi4EPAq3nrngQ+Fi7/HvBMuDyZYPry2eH7M8J/q4E3gXlALfAScG7K+vY/gLXhcgPB1O21wFRgV/jvlHB5SoX2bTNwSbh8E3BXuHxueE7qgLnhuaquxPM2wn5dALw3XP4toDXvmOeAZYABP82d87T0LW/7OuBHwK3h+4o7Z6M4bzXAy8D54ftpQHVS500VxElw92cJ/hj2Ww3k0n4SsD9c/kPg/7r73vDYw+H6JcBOd9/l7seBx4BrEm14EUbYNwdOD/+P5bTwuB7gSuApd2939wzwFLA86bYPp0DfzgGeDZefAv4gXL4GeMzdu9x9N7CT4JxV3HkbSb/c/QV3z52/bUB9+ICuWcBEd9/owV+dB4GPJ974YYzwnGFmHwd2E/Qtp+LOGYy4b1cAL7v7S+Gxbe7em9R5U0DE74vAX5nZPuCvga+E688BppjZM2FJ/+lwfSOwL+/4lnBdJfoi0X27B/gAQWC8AnzB3U+Qrr5t490/Fjfw7iNvC/UhLX0r1K98fwA87+5dBH1oydtWqf2CAn0zs9OA24G/GLB/Ws4ZFD5v5wBuZhvM7Hkz++/h+kTOmwIifn8KfMndzwS+RPDUPAhKw0XA7xP8n/X/NLNzytPEUSvUtyuBF4H3Ah8E7hl4PTgFbgL+zMy2AqcDx8vcnrgM2S8zOw9YC/xxGdp2sgr17evA37r7sXI1LAaF+lYD/A7wR+G/15rZ5Uk1IrFHjp7CVgFfCJd/BPxDuNwCtLn7b4DfmNmzwPnh+vz/q2sCWkvU1pEq1LfPAN8KS9udZrYbeD9BPy7NO74JeKYkLR0hd/8VQflOGNy/H25qpfD5qfjzNkS/MLMm4J+BT7v7m+HqVoK+5FRkv2DIvi0FrjezbxOM/Z0ws05gKyk4ZzBk31qAZ939rXDbEwTjFw+TwHlTBRG//cAl4fJlwBvh8r8Cv2NmNWY2geCX+DWCwagFZjbXzGoJnsu9vsRtLlahvu0FLgcwsxnA+wgGpDcAV5jZFDObQvALv6GkLS6SmZ0R/lsFfBX4+3DTemBleH1+LrCAYDAwFeetUL/MbDLwE2C1u/88t7+7HwB+bWbLwjGlTxP87lacQn1z99919znuPgf4O+Bud7+HlJwzGPL3cQPw22Y2wcxqCP573J7YeSv3CH6aX8CjwAGgmyDZP0tQ9m0luENiE7Aob//bCO5kehX4Yt763wNeJ7jD4o5y92ukfSO4tPQkwfjDq8Cn8n7OTQQDuzuBz5S7X0P07QvhOXgd+BbhLAPh/neE52YHeXeGVNp5G0m/CP7o/Ibg0mDulbuzbnF4Ht8kGF+ycvVptOcs77ivE97FVInnbJS/j58iGKN4Ffh23vrYz5um2hARkUi6xCQiIpEUECIiEkkBISIikRQQIiISSQEhIiKRFBAiIhJJASFSQcysutxtEMlRQIiMkpmtMbMv5r3/ppl9wcxuM7PNZvaymf1F3vZ/CSdq3GZmN+etP2Zmf2NmLwEXlbYXIoUpIERG7wGCKQ1yUyKsBA4STMexhGDiwkVmdnG4/03uvojgG6//zcymhevfA2xy9/Pd/f+VsP0iQ9JkfSKj5O7NZtZmZhcAM4AXgAsJ5px6IdztNILAeJYgFK4N158Zrm8DeoHHS9l2kWIoIEROzj8ANwIzCSqKy4G/dPf78ncys0uBjwAXufs7ZvYMMD7c3OnuvSVqr0jRdIlJ5OT8M8FT8i4kmGlzA3BT+NAazKwxnJlzEpAJw+H9BI+GFKloqiBEToK7Hzezp4G3wyrgSTP7APDL8Jnxxwhm3/x34E/M7DWCWWE3lqvNIsXSbK4iJyEcnH4euMHd3xhuf5E00SUmkVEys3MJnnPxHwoHGYtUQYiISCRVECIiEkkBISIikRQQIiISSQEhIiKRFBAiIhLp/wOv6xv21eD2EAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x[0] + np.arange(n), tau_posterior, color='black')\n",
    "plt.xlabel('year')\n",
    "plt.ylabel('tau')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "81ddd53f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior mode: 1892\n"
     ]
    }
   ],
   "source": [
    "tauMAP = x[0] + np.argmax(tau_posterior)\n",
    "print(f'Posterior mode: {tauMAP}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97bb156b",
   "metadata": {},
   "source": [
    "##### Lambdas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "21545db7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4U0lEQVR4nO3deXhU5dn48e+dhSSAIISdBAKyKCBliQhSFRSpKIJWUNkEtaKlVG154afUClVr30qlLq+4USsILixRERWKgliVXSEIKAYECXsSViFAkvv3x8zEELJMkjlzMpn7c11zMXPOc55zT0LmnvOcZxFVxRhjTPiKcDsAY4wx7rJEYIwxYc4SgTHGhDlLBMYYE+YsERhjTJiLcjuAsqpXr54mJSW5HYYxxoSUdevWZahq/aL2hVwiSEpKYu3atW6HYYwxIUVEdha3z5qGjDEmzDmWCEQkVkRWi8gGEdkkIn8poswoETkoIuu9j984FY8xxpiiOdk0dAq4SlWPi0g08LmIfKSqKwuVe1tVxzoYhzHGmBI4lgjUM3fFce/LaO/D5rMwxpTqzJkzpKenk52d7XYoISc2NpaEhASio6P9PsbRm8UiEgmsA1oBz6vqqiKK3SwiVwBbgT+o6q4i6hkNjAZo1qyZgxEbYyqD9PR0zjvvPJKSkhARt8MJGapKZmYm6enptGjRwu/jHL1ZrKq5qtoJSAC6iUiHQkXeB5JUtSOwBJhRTD0vq2qyqibXr19k7ydjTBWSnZ1NfHy8JYEyEhHi4+PLfCUVlF5DqnoYWAZcW2h7pqqe8r6cDnQNRjzGmMrPkkD5lOfn5mSvofoicr73eRxwDfBtoTKNC7wcAGxxKp6q7PTp00yfPp3c3Fy3QzHGhCAnrwgaA8tEJBVYAyxR1YUi8qiIDPCWuc/btXQDcB8wysF4qqynnnqKu+++m9mzZ7sdijFVRs2aNQNSz+TJk/nHP/5RarlRo0Yxb968EsvMnTuX9u3bExEREdCBtU72GkoFOhex/ZECzx8CHnIqhnDx3XffAXDw4EGXIzHGOKlDhw6kpKRwzz33BLReG1lcBRw5cgSAw4cPuxuIMVXQ8ePHufrqq+nSpQsXX3wx7733HgA7duzgwgsvZNSoUbRp04Zhw4bx8ccf07NnT1q3bs3q1avz69iwYQM9evSgdevWvPLKK4Cnh8/YsWNp27Ytffr04cCBA/nlH330US655BI6dOjA6NGj8a0kedFFF9G2bduAv8eQm2vInCslJYWsrCzi4+PdDsWYgHvggQdYv359QOvs1KkTTz/9tF9lY2Njeeedd6hVqxYZGRl0796dAQM8rdtpaWnMnTuXV199lUsuuYQ33niDzz//nAULFvDEE0/w7rvvApCamsrKlSv56aef6Ny5M9dffz0rV67ku+++Y/Pmzezfv5927dpx5513AjB27FgeecTTeDJixAgWLlzIDTfcENCfQUGWCKoAX5cxY0zgqSoTJ07ks88+IyIigt27d7N//34AWrRowcUXXwxA+/btufrqqxERLr74Ynbs2JFfx8CBA4mLiyMuLo7evXuzevVqPvvsM4YMGUJkZCRNmjThqquuyi+/bNkynnzySU6cOEFWVhbt27e3RGCK95///IcZM2Zwww03MGPGDObMmcN5553ndljGBIy/39ydMnv2bA4ePMi6deuIjo4mKSkpv59+TExMfrmIiIj81xEREeTk5OTvK9yls6QuntnZ2YwZM4a1a9eSmJjI5MmTHR9hbfcIQtwnn3zC3LlzAVi0aBFbt251OSJjqpYjR47QoEEDoqOjWbZsGTt3Fjubc7Hee+89srOzyczM5NNPP+WSSy7hiiuu4O233yY3N5e9e/eybNkygPwP/Xr16nH8+PFSexIFgiWCEPfVV19x8cUX07FjR+DnHkTGmMAYNmwYa9eu5eKLL2bmzJlceOGFZa6jY8eO9O7dm+7du/PnP/+ZJk2acNNNN9G6dWvatWvH7bffTo8ePQA4//zzufvuu+nQoQO/+tWvuOSSS/Lreeedd0hISGDFihVcf/31/OpXvwrIexTf3ehQkZycrLYwzc/atm1Lp06dmDlzJtWrV+dPf/oTjz76qNthGVMhW7Zs4aKLLnI7jJBV1M9PRNapanJR5e2KIMRlZGRQv359YmJiaNmypV0RGGPKzG4WhzBVpWnTpvmzDF5++eU2P4sxpswsEYQwESE1NTX/9auvvupiNMaYUGVNQ8YYE+YsEYSwdevW8ctf/jJ/1OW7777LL37xCzIzM90NzBgTUiwRhLBdu3bxxRdfkJeXB3imo05NTWXv3r0uR2aMCSWWCEKY75t/vXr1AGjUqBEA+/btcy0mY6qKyjgN9fjx47nwwgvp2LEjN910U8AmmrREEMIyMjIA8ucZatzYs86PJQJjqqZrrrmGb775htTUVNq0acPf/va3gNRriSCEZWRkEBsbS/Xq1YGfrwisaciYwKlM01D37duXqChPZ8/u3buTnp4ekPdo3UdDWKNGjejVq1f+2IGaNWvSp08fGjZs6HJkxgRWr169ztl2yy23MGbMGE6cOMF11113zv5Ro0YxatQoMjIyGDRo0Fn7Pv30U7/PXVmnoX711Ve59dZb/X4fJbFEEMLGjRvHuHHj8l+LCEuWLHExImOqnso4DfVf//pXoqKiGDZsWEDeoyUCY0ylV9I3+OrVq5e4v169emW6Aiissk1D/dprr7Fw4UI++eSTgM0kYPcIQtjVV1/NxIkTz9r2wAMPcNlll7kUkTFVT2WahnrRokU8+eSTLFiwIP/eYCA4dkUgIrHAZ0CM9zzzVHVSoTIxwEygK5AJ3KqqO5yKqapZv379OTMM5uTksGXLFpciMqbqGTZsGDfccAMXX3wxycnJFZqGOiMj46xpqJcuXUq7du1o1qxZkdNQN2rU6KxpqMeOHcupU6e45pprAM8N4xdffLHib1JVHXkAAtT0Po8GVgHdC5UZA7zofX4b8HZp9Xbt2lWNak5OjoqITpo06aztjz/+uAJ68uRJdwIzJgA2b97sdgghraifH7BWi/lcdaxpyHvu496X0d5H4cUPBgIzvM/nAVeLTZ/pl0OHDqGq+YPJfGxQmTGmrBy9RyAikSKyHjgALFHVVYWKNAV2AahqDnAEOGcVdhEZLSJrRWTtwYMHnQw5ZPgGkxVOBDaozBhTVo4mAlXNVdVOQALQTUQ6lLOel1U1WVWT69evH9AYQ1VkZCQ33HADF1xwwVnbW7VqxZAhQwI2PN4Yt2iIrZ5YWZTn5xaU7qOqelhElgHXAt8U2LUbSATSRSQKqI3nprEpRevWrVmwYME529u0acMbb7zhQkTGBE5sbCyZmZnEx8fbYktloKpkZmYSGxtbpuOc7DVUHzjjTQJxwDXA3wsVWwCMBFYAg4Clal8DAiInJyd/KLoxoSYhIYH09HSsKbjsYmNjSUhIKNMxTn5SNAZmiEgkniaoOaq6UEQexXP3egHwL+B1EUkDsvD0HDJ+eOaZZ/j73/9OWlraOf2JO3fuTLt27Zg9e7ZL0RlTMdHR0flLsBrnOZYIVDUV6FzE9kcKPM8GBjsVQ1W2b98+MjIyiIuLO2df7dq1AzYZlTGm6rORxSEqKyuLunXrFtl+6rusNsYYf1giCFG+RFAUXyKw2y3GGH9YIghRpSWC06dP2402Y4xfrFtJiLrqqquK7RXUvXt3JkyYYN3ujDF+kVBrPkhOTta1a9e6HYYxxoQUEVmnqslF7bOmoRCVl5dX4v5jx44FbGFrY0zVZokgBJ0+fZqoqCimTJlS5P68vDzi4+N58skngxyZMSYUWSIIQb6ZR2vUqFHk/oiICJo0aWJdSI0xfrFEEIKysrIAiu01BDaWwBjjP0sEIcifRNCkSRN2794drJCMMSHMEkEI8icRNGzYkAMHDgQrJGNMCLNxBCGoWbNm3H///SQmJhZb5uabb6Z169aoqo0nMMaUyMYRGGNMGLBxBFXM8ePHyc7OLrFMdnY2mzdv5tixY0GKyhgTqiwRhKAJEybQtGnTEsusX7+e9u3b8/nnnwcpKmNMqLJEEIKysrKIj48vsUyDBg0A2L9/fzBCMsaEMEsEIaikmUd9fInAeg4ZY0pjiSAE+ZMIatasSfXq1e2KwBhTKksEIcifRACeqwJLBMaY0tg4ghB0//33+7Ww95QpU2jYsGEQIjLGhDLHxhGISCIwE2gIKPCyqj5TqEwv4D3gB++mFFV9tKR6bRyBMcaUXUnjCJy8IsgBxqnqVyJyHrBORJao6uZC5f6rqv0djKNKOXPmDOnp6TRq1Ii4uLgSy/74449s376dXr16BSc4Y0xIcuwegaruVdWvvM+PAVuAkju/m1Lt3LmTli1bMm/evFLLvvTSS/Tp06fURWyMMeEtKDeLRSQJ6AysKmJ3DxHZICIfiUj7Yo4fLSJrRWRtuC/I7s+Ecz4NGzYkNzc3/xhjjCmK44lARGoC84EHVPVood1fAc1V9RfAc8C7RdWhqi+rarKqJtevX9/ReCu7Q4cOAVCnTp1Sy9qgMmOMPxxNBCISjScJzFbVlML7VfWoqh73Pv8QiBaRek7GFOrKekUANqjMGFMyxxKBeOY+/hewRVWnFlOmkbccItLNG0+mUzFVBb5EUJYrAksExpiSlNprSETiVbU8H849gRHARhFZ7902EWgGoKovAoOA34pIDnASuE1DbV7sILv88st5+umn/boiaNGiBQsWLKBr165BiMwYE6pKHUcgIt8D64F/Ax+5/UFt4wiMMabsKroeQRvgZTzf7r8XkSdEpE0gAzT+S0tLIy0tze/yixYtYt26dQ5GZIwJdaUmAvVYoqpDgLuBkcBqEVkuIj0cj9CcZdy4cdx8881+l7/rrruYNm2agxEZY0KdX/cIgOF4rgj2A78HFgCdgLlA6ZPemIA5dOiQX/cHfOrUqZPf5dQYY4riT9PQCqAWcKOqXq+qKaqao6prgRedDc8UlpWV5VePIZ+6detaIjDGlMifRPCwqj6mqum+DSIyGEBV/+5YZKZI5bkisJHFxpiS+JMIHixi20OBDsT4p6xXBNY0ZIwpTbH3CESkH3Ad0FREni2wqxaemUVNkKkqM2bMoHXr1n4fM3HiRO6//34HozLGhLqSbhbvAdYCA4CC/Q+PAX9wMihTNBHhlltuKdMxbdpYT19jTMmKTQSqugHYICKzVdWuACqBI0eOsHbtWrp06eJ381BaWhpLly5l6NCh1KxZ0+EIjTGhqNh7BCIyx/v0axFJLfDYKCKpQYrPFJCamkqfPn3KNEBs1apV3HPPPezevdvByIwxoaykpiFfw7KtHlZJlGUKah9fWbthbIwpTrFXBKq61/s0A9ilqjuBGOAXeO4fmCAryxTUPr6ylgiMMcXxp/voZ0CsiDQF/oNnhPFrTgZlimZXBMYYJ/iTCERVTwC/Bqap6mCgyCUljbOysrKIiIigVq1afh/jSwQ2qMwYU5xS5xrCs8ZMD2AYcJd3W6RzIZnijBw5kssuu4yICP/XE4qPj2fjxo0kJiY6GJkxJpT5kwjuxzOS+B1V3SQiLYFlzoZlitKqVStatWpVpmMiIyPp0KGDQxEZY6qCUhOBqn6G5z6B7/V24D4ngzJF++STT4iNjaVnz55lOu61116jTp06DBw40KHIjDGhzJ9pqNsA/wMkFSyvqlc5F5YpysSJE6lTpw6LFi0q03H//Oc/SUpKskRgjCmSP01Dc/FMNz0dyHU2HFOSrKwsLrjggjIfZxPPGWNK4k8iyFHVF8pasYgkAjOBhoACL6vqM4XKCPAMnsntTgCjVPWrsp4rXJR15lGfOnXqlGl5S2NMePGn+8n7IjJGRBqLSF3fw4/jcoBxqtoO6A78TkTaFSrTD2jtfYwGypxwwkVeXh6HDx8u02Ayn7p161r3UWNMsfy5Ihjp/Xd8gW0KtCzpIO/I5L3e58dEZAvQFNhcoNhAYKaqKrBSRM4XkcYFRjUbr6NHj5KXl1euRGBNQ8aYkvjTa6jCaxKLSBLQGVhVaFdTYFeB1+nebZYICqlRowZr1qyhSZMmZT724YcfZuLEiQ5EZYypCvzpNVQd+CPQTFVHi0hroK2qLvTnBCJSE5gPPKCqR8sTpIiMxtN0RLNmzcpTRciLjo4mOTm5XMeef/75gQ3GGFOl+HOP4N/AaeAy7+vdwOP+VC4i0XiSwGxVTSmiyG6g4JDXBO+2s6jqy6qarKrJ9evX9+fUVc6OHTuYPn06GRkZZT7222+/5aGHHmLPHpsr0BhzLn8SwQWq+iRwBsA775CUdpC3R9C/gC2qOrWYYguA28WjO3DE7g8UbfXq1dx9993s37+/zMfu3LmT//3f/2XHjh2BD8wYE/L8uVl8WkTi8NwgRkQuAE75cVxPPDOVbhSR9d5tE4FmAKr6IvAhnq6jaXi6j95RluDDSXmmoPaxieeMMSXxJxFMAhYBiSIyG88H/KjSDlLVzynlysHbW+h3fsQQ9sozBbWPTUVtjCmJP72GlojIV3jGAghwv6qWvaHaVEhWVhZxcXHExsaW+VhLBMaYkhSbCESkS6FNvrb7ZiLSzEYAB9ehQ4fKdTUAP/caOnz4cOACMsZUGeJpnSlih4hvqulYIBnYgOeKoCOwVlV7BCXCQpKTk3Xt2rVunNpVGRkZHDp0iNatW5fr+BMnThAXF4fnHr4xJtyIyDpVLbIPerFXBKra23twCtBFVTd6X3cAJjsQpylBvXr1qFevXrmPr169egCjMcZUJf50H23rSwIAqvoNcJFzIZmivPjii3zwwQflPn7q1Kk8/fTTgQvIGFNl+JMIUkVkuoj08j5eAVKdDsyc7a9//Svz588v9/EffPAB8+bNC2BExpiqwp9EcAewCc+SlffjmTTO+vsHWUVuFoOn55CNIzDGFMWf7qPZwD+9D+OC06dP89NPP1U4EVj3UWNMUfy5IjAuq8hgMh9LBMaY4lgiCAGBSATx8fHExsaSnZ0dqLCMMVVEseMIAEQkEvi7qv5P8EIqWTiOI8jLy+PQoUPExcVZN1BjTLmUaxwBgKrmisgvnQnL+CsiIoL4+Hi3wzDGVFH+NA19LSILRGSEiPza93A8MpNv1apVTJw4sUJTRGzZsoVbb72Vb775JnCBGWOqBH8SQSyQCVwF3OB99HcyKHO2VatW8be//Y2cnJxy1/HTTz8xZ84ctm/fHsDIjDFVgT/dR23MgMt8N4srsuSkbx0D6zlkjCms1CsCEUkQkXdE5ID3MV9EEoIRnPE4dOgQtWrVIirKn+UjimZTURtjiuPvmsULgCbex/vebSZIsrKyKtR1FKB27dqIiI0uNsacw59EUF9V/62qOd7Ha0B4riDvkmPHjlU4EURERNCmTRuio6MDFJUxpqrwp60hU0SGA296Xw/Bc/PYBElKSgqnT5+ucD3ffvttAKIxxlQ1/lwR3AncAuzDs0rZIGzSuaASEWJiYtwOwxhTRZWaCFR1p6oOUNX6qtpAVW9U1R9LO05EXvXeXC6y47p3SusjIrLe+3ikPG8gHNx3333MmTOnwvVMmjSJ3/72twGIyBhTlZS0ZvEEVX1SRJ4DzpmHQlXvK6Xu14D/A2aWUOa/qmpjEkqgqrz00kvExcVxyy23VKiuLVu2kJpqS0kYY85W0j2CLd5/yzWxj6p+JiJJ5TnW/OzkyZOcPn06fxxARdgMpMaYopS0ZvH73knnLnZw0rkeIrIB2AP8j6puKqqQiIwGRgM0a9bMoVAqp0DMPOrjSwSqaovYG2PylXiPQFVzgZ4OnfsroLmq/gJ4Dni3hDheVtVkVU2uXz+8eq76+v0HIhHUrVuXM2fOcOLEiQrXZYypOvzpPrpeRBYAc4GffBtVNaUiJ1bVowWefygi00SknqpmVKTequbEiRPUqlUrIE1DzZs3p1OnTpw4cYIaNWoEIDpjTFVQ4noEACJS1ChiVdU7S63cc49goap2KGJfI2C/qqqIdAPm4blCKDGgcFyPwBhjKqrc6xFA+SedE5E3gV5APRFJByYB0d46X8QzHuG3IpIDnARuKy0JGGOMCTx/Jp1rIyKf+MYDiEhHEXm4tONUdYiqNlbVaFVNUNV/qeqL3iSAqv6fqrZX1V+oandV/bLib6fqSUlJ4bbbbgvIEpPbtm2jW7duLFmyJACRGWOqCn9GFr8CPAScAVDVVOA2J4MyP/v666+ZM2cO1apVq3BdkZGRrFmzhl27dgUgMmNMVeFPIqiuqqsLbSv/CimmTHwzj0ZE+POrKplNRR1eDh8+zPvvv8+qVavcDsVUcv58umSIyAV4RxeLyCA8cw6ZIDh06FBAuo4C1KpVi8jISEsEYWDhwoU0atSIAQMG0L17dyZPnkxeXp7bYZlKyp9E8DvgJeBCEdkNPADc62RQ5mdZWVkB6ToKnsnrzj//fEsEVdwXX3zB4MGD6dChA0uXLuX222/nL3/5C88995zboZlKyp9xBKqqfUSkBhChqsdEpIXTgRmP2rVrB+yKAOCqq66iefPmAavPVC7Hjh3j1ltvJTExkY8++oj69evTq1cv9uzZw6OPPsrIkSMrtOSpqZr8GUfwlap2KbRtnap2dTSyYtg4AmOK98ADD/Dss8+yYsUKLr300vzt69evp0uXLowbN44pU6a4GKFxS0njCIptGhKRC0XkZqC2iPy6wGMUEOtQrMaYcvrqq6947rnnuPfee89KAgCdOnVi2LBhTJs2jaNHjxZTgwlXJd0jaAv0B84Hbijw6ALc7XhkBlWlV69ezJo1K2B1TpgwgSuuuCJg9ZnKITc3l3vuuYf69evzxBNPFFlm7NixnDhxgrfeeivI0ZnKrqTZR98D3hORHqq6IogxGa9jx46xfPlybrjhhoDVefz4cTZv3hyw+kzl8MILL7B27VreeOONYu8BdOvWjQ4dOjB9+nRGjx4d3ABNpeZPr6GbRKSWiER7Rxgf9K5hbBwWyCmoferXr09WVhY5OTYUpKr4/vvvefDBB+nbty+33Vb8WE8R4Te/+Q1r1qxhw4YNQYzQVHb+JIK+3plC+wM7gFbAeCeDMh6BnILap2HDhqgqGRk2yWtVcObMGYYPH061atWYPn16qetMDB8+nOjoaF5//fUgRWhCgT+JINr77/XAXFU94mA8pgDfFUGgxhGAJxEA7N+/P2B1Gnfk5uZy5513snr1al566SUSExNLPSY+Pp5+/frx5ptvkpubG4QoTSjwJxG8LyLfAl2BT0SkPlDxGdBMqaKioujUqVP+h3cgtG7dml//+tcBmbvIuOfw4cMMHTqUWbNm8de//pXBgwf7fezQoUPZs2cP//3vfx2M0ISSUscRAIhIXeCIquaKSHWglqruczy6Itg4AhPOTp8+zcyZM3nkkUfYv38/TzzxBP/v//2/MtVx4sQJGjRowNChQ3n55ZcditRUNuUaR1Dg4GhgOPC2iMwD7gIyAxuiCTZb+iH0LF68mIsuuoi7776bhIQEVq1aVeYkAFC9enVuuukm5s2bx+nTpx2I1IQaf5qGXsDTLDTN++ji3WYc9tRTT3H55ZcHvN5mzZrx0EMPBbxe45y//OUvXHvttURHR/PBBx+watUqkpOL/HLnl6FDh3Lo0CEWL14cwChNqPJnrqFLvAvM+ywVEet7FgRbt27l+++/D3i9IsLevTaBbKh46qmnmDx5MqNGjeLFF18kJiamwnX26dOHevXq8cYbbwR0nIoJTf5cEeR6p6EGQERaAtbdIAh8axEEWoMGDazXUIhYu3Yt48ePZ9CgQUyfPj0gSQAgOjqaW265hffee4/jx48HpE4TuvxJBOOBZSLyqYgsB5YC45wNy4Cn+2ggu476NGzY0BJBCMjNzWXMmDE0bNiQ6dOnExkZGdD6hw4dysmTJ3n33XcDWq8JPaUmAlX9BGgN3Af8HmirqsucDswEdlGagiwRhIYZM2awZs0a/vGPf1C7du2A19+jRw9atWrFSy+9FPC6TWjxp9dQLJ7FaSYDk4DfereVdtyrInLAt+h9EftFRJ4VkTQRSRWRLkWVC2edOnWiW7duAa+3b9++DBkyJOD1msDJy8tjypQpdO3alaFDhzpyjoiICH7729/y+eefk5qa6sg5TGjwZz2COcAxwDcF5lDgfFUtcQSLiFwBHAdmqmqHIvZfh+cK4zrgUuAZVb20cLnCbByBCQcfffQR1113HbNmzWLYsGGOnScrK4umTZty++2325VBFVehcQRAB1W9S1WXeR93A+1LO0hVPwOySigyEE+SUFVdCZwvIo39iMcEwKlTpzhz5ozbYZhi/POf/6RJkyZlGjFcHnXr1mXIkCHMmjXLmgvDmD+J4CsR6e57ISKXAoH4St4U2FXgdbp32zlEZLSIrBWRtQcPHgzAqSu/gwcPkpiYyJtvvhnwupcvX05sbCyff/55wOs2Fbdt2zaWLFnCmDFjgjIVyIMPPsipU6eKXcfAVH3+JIKuwJciskNEdgArgEtEZKOIBKVhUVVfVtVkVU2uX79+ME7puoMHD5Kenl7qbJLlUa9ePcAmnqusZs2ahYgwcuTIoJyvTZs2+WMUdu7cGZRzmsrFn0RwLdACuNL7aOHd1h/PimXltRsoOF1ignebATIzPbN4+D60A8lmIK28VJXXX3+d3r17k5CQELTzTpo0CRFh3DjrGR6O/Ok+urOkRwXOvQC43dt7qDueSe1suKuXb70AJxJB3bp1iYmJIT09PeB1m4pZuXIl27ZtY8SIEUE9b2JiIpMmTWL+/PnMnTs3qOc27vPniqBcRORNPM1IbUUkXUTuEpF7ReReb5EPge1AGvAKMMapWEKRLxHEx8cHvO6IiAiaNWtmzQCV0OzZs4mLi+Pmm28O+rnHjx9P165dGTNmDPv2uTK5sHGJP3MNlYuqlthRXT39Vn/n1PlDXUJCAjfddJMjVwQADzzwgCOjlk355eXlkZKSQr9+/TjvvPOCfv6oqChmzpxJcnIyI0aMYPHixUREOPZd0VQi9luupPr160dKSgpxcXGO1D9mzJgS17c1wbdq1Sr27t3Lr3/9a9diaNeuHc8++ywff/wxU6ZMcS0OE1yWCMLU6dOn+f77720R+0pk/vz5REdH079/f1fjuOuuuxg0aBCPPPII33xT5MQApoqxRFBJ3XjjjVx99dWO1f/WW2/Rpk0btm/f7tg5jP9UlZSUFPr06ePIvEJlISJMmzaN2rVrM3LkSPuyEAYsEVRS+/fvD/hskwU1b94cwG4YVxLr16/nhx9+cOUmcVHq16/P888/z1dffcW///1vt8MxDrNEUEllZGQ40mPIxxJB5ZKSkkJERAQDBgxwO5R8gwYNomfPnjzyyCO2ZkEVZ4mgksrMzHSsxxBA06ZNiYiIsERQScyfP58rr7ySyjRyXkSYMmUK+/bt45///Kfb4RgHWSKohHJycjh06JCjiSA6OpqmTZtaIqgEtmzZwpYtW1ztLVScHj160L9/f5599llOnjzpdjjGIZYIKqEzZ84wduxYunfvXnrhCvj73//OnXfe6eg5TOlSUlIAuOmmm1yOpGjjxo0jIyODWbNmlV7YhKRS1yOobGw9AlPVdOnShZiYGFasWOF2KEVSVbp06cKpU6fYtGmTIxMhGudVdD0CE2RnzpwJyloBmZmZLF261LoHumjbtm18/fXXjq87UBEiwh/+8Ae2bNnC0qVL3Q7HOMASQSX0wQcfUK1aNdavX+/oed577z2uvvpqfvzxR0fPY4o3f/58gErTbbQ4gwcPpnbt2taVtIqyRFAJOTnhXEEXXHABAGlpaY6exxRv3rx5XHLJJfndeSuruLg4hgwZwvz58zly5Ijb4ZgAs0RQCQUrEbRq1QrwNE+Y4Nu5cydr1qxh0KBBbofilzvuuIPs7Gzefvttt0MxAWaJoBI6cOAA1atXp3r16o6ep3HjxsTGxloicEmoNAv5XHLJJbRr144ZM2a4HYoJMEsEldDevXtp3Lix4+eJiIigZcuW1jTkknnz5tGpU6f8JrrKTkQYPnw4X375JTt27HA7HBNAlggqoYEDB/L73/8+KOd64YUXeOyxx4JyLvOz9PR0VqxYETLNQj6+qcvfeustlyMxgWTjCIxxwbPPPsv999/Pt99+S9u2bd0Op0x69uzJsWPHSE1NdTsUUwY2jiDE7Ny5k+zs7KCca8+ePcyYMYPDhw8H5XzG4+2336ZDhw4hlwQAhgwZwsaNG9m4caPboZgAsURQyZw8eZKkpCSmTp0alPNt3LiRUaNG2R91EG3bto0vv/yS4cOHux1KuQwePJjIyEjefPNNt0MxAWKJoJLZv38/QFBuFoONJXDD66+/jogwbNgwt0Mpl4YNG9KnTx/efPNNQq1p2RTN0UQgIteKyHcikiYiDxaxf5SIHBSR9d7Hb5yMJxTs3bsXgEaNGgXlfM2bNycyMtISQZCoKrNmzaJ3794kJCS4HU65DR06lB07drBy5Uq3QzEB4FgiEJFI4HmgH9AOGCIi7Yoo+raqdvI+pjsVT6jYt28fELxEEB0dTcuWLdm6dWtQzhfuvvjiC7Zt28aIESPcDqVCbrzxRmJjY3njjTfcDsUEgJNXBN2ANFXdrqqngbeAgQ6er0rwXREEq2kIoG3btnz33XdBO184e+mll6hVq1bIdRstrFatWvTv35+3337bJi2sApxMBE2BXQVep3u3FXaziKSKyDwRSSyqIhEZLSJrRWTtwYMHnYi10ujZsyf/+Mc/grpS1TPPPMOSJUuCdr5wlZmZydy5cxk+fDg1a9Z0O5wKGzp0KAcPHuSTTz5xOxRTQW7fLH4fSFLVjsASoMix66r6sqomq2pyZVrKzwm/+MUvGDdunKML1xfWsmVLGjZsGLTzhauZM2dy6tQp7rnnHrdDCYh+/fpRu3Ztax6qApxMBLuBgt/wE7zb8qlqpqqe8r6cDnR1MJ6QsHnz5qBPC52RkcHjjz9uA4QclJuby/PPP89ll11Gx44d3Q4nIGJjY7n55ptJSUmxZSxDnJOJYA3QWkRaiEg14DZgQcECIlKwIXwAsMXBeELCyJEjg/6NMTc3lz//+c98+umnQT1vOHnnnXfYtm0bf/zjH90OJaCGDh3K8ePHeffdd90OxVSAY4lAVXOAscBiPB/wc1R1k4g8KiIDvMXuE5FNIrIBuA8Y5VQ8oSJYE84V1KBBA2rXrm03jB2iqjz55JO0atWKG2+80e1wAqp3794kJSXxyiuvuB2KqYAoJytX1Q+BDwtte6TA84eAh5yMIZTk5eWxf//+oHUd9RER6znkoGXLlrFmzRpeeOGFoN77CYaIiAjuvvtu/vSnP7F161batGnjdkimHNy+WWwK2LNnDzk5OSQmFtl5ylGWCJyhqvzpT3+iadOmjBo1yu1wHHHHHXcQFRVlVwUhzBJBJeL7IHZjIrK2bduSmZlpN/0CbOHChaxcuZJJkyYRGxvrdjiOaNy4MQMHDmT69OkcO3bM7XBMOdg01JXIgQMHWLp0Kb/61a+oU6dOUM994sQJYmJiqlzThZvOnDlD586dOXXqFJs3byY6OtrtkByzcuVKevTowVNPPVXlbohXFTYNdYho0KABt912W9CTAED16tUtCQTYtGnT2LRpE1OmTKnSSQCge/fu9OrVi6lTp3L69Gm3wzFlZImgElm8eDEbNmxw7fwTJkzg6aefdu38VcmBAweYNGkSffv2ZeDA8JhZ5cEHH2T37t28+uqrbodiysiahiqRpKQkLrvsMtdGavbs2ZOoqCiWL1/uyvmrkltvvZV3332XDRs2cOGFF7odTlCoKldeeSXffvstaWlp1KpVy+2QTAHWNBQCTp48yY8//ujqilUdO3YkNTXV5pivoPnz5zNnzhwmTZoUNkkAPN2Qp06dysGDB3niiSfcDseUgSWCSiItLQ1VdT0RHD58mPT0dNdiCHV79uzh3nvvpUuXLowfP97tcIIuOTmZ22+/nalTp7razGnKxhJBJeFm11Ef3xw4NudQ+eTm5jJ06FBOnDjB7Nmzq/wN4uJMnTqVunXrcvvtt3Pq1KnSDzCus0RQSfgSgZsjMzt06EBSUhI//fSTazGEsokTJ7J8+XJeeOGFsGoSKiw+Pp7p06eTmprKhAkT3A7H+ENVQ+rRtWtXrYpOnDihy5YtczsMU04zZsxQQO+99163Q6k0HnjgAQX0tddeczsUo6rAWi3mc9V6DZlz+P5zRETYBaM/PvroIwYOHMjll1/OokWLwrZJqLCcnBz69u3LF198waJFi+jdu7fbIYU16zVUya1atYqHH36YrKwst0Nhy5YtJCYmsmjRIrdDCQnLly/n5ptvpkOHDqSkpFgSKCAqKop58+bRunVrBg4ciH2Bq7wsEVQCc+fOZcqUKcTExLgdCi1atODQoUN88MEHbodS6f3nP/+hX79+JCUlsXjxYmrXru12SJVO3bp1WbRoEXXr1uWaa65h9erVbodkimCJoBJYvHgxV1xxBTVq1HA7FGJjY+nTpw8ffvihjScowaxZs+jfvz9t2rRh+fLlQV1jOtQkJCSwfPly6tSpQ58+fexqsxKyROCy7777jm+++YZrr73W7VDyDRgwgB07djBnzhy3Q6l0VJXHHnuMESNG0LNnT5YtW2ZJwA/Nmzfns88+o2XLlvTv35+pU6eSl5fndljGp7i7yJX1UdV6Dd19990aExOj+/fvdzuUfKdPn9ZLL71U4+Pj9dixY26HU2mcPHlSR40apYCOGDFCT5065XZIIefo0aN64403KqB9+/bVrVu3uh1S2KCEXkN2ReAi9Ta9jB49mgYNGrgczc+io6OZPXs277zzDjVr1iQ7O5udO3e6HZar0tLSuOyyy3jttdeYNGkSM2bMoFq1am6HFXLOO+88UlJSmDZtGl9++SXt27fnnnvuYdOmTW6HFt6KyxCV9VHVrghUVfPy8twOoUSPP/64xsXF6UsvveR2KK6YN2+e1qpVS+vUqaPvv/++2+FUGXv37tV77rlHY2JiFNB27drpQw89pCtXrtScnBy3w6tyKOGKwPUP9rI+qkoiWL9+va5YscLtMPyyc+dOveaaa1REdOHChW6HEzS7du3SwYMHK6DdunXTHTt2uB1SlXTw4EF95pln9KqrrtLIyEgFtE6dOnrzzTfrCy+8oN9//32l/7IUClxLBMC1wHdAGvBgEftjgLe9+1cBSaXVWRUSwbFjx7Rjx46akJCg2dnZbofjl59++km7dOmi5513nq5atcrtcBy1efNmHTt2rMbExGhMTIw+9thjdj8gSDIzM/WNN97QO++8UxMTExVQQJs3b6533XWXvvnmm3rgwAG3wwxJJSUCx0YWi0gksBW4BkgH1gBDVHVzgTJjgI6qeq+I3AbcpKq3llRvKI8szsvL4+uvv+aOO+5g06ZNvPfee/Tv39/tsPyWnp7OlVdeycmTJ9m2bRtxcXEVqi8vL4/t27fz3//+l++++46MjAwSEhK47rrrSEpKolq1asTExBATE+PIKOfc3FyysrL4/vvv2bp1Kxs3buTjjz8mNTWVqKgoRo4cycMPP0xSUlLAz21Kp6qkpaWxZMkSPv74Y5YuXcqRI0cA6NSpE3369OGyyy6jWbNmJCQk0KBBA0TE5agrr5JGFjuZCHoAk1X1V97XDwGo6t8KlFnsLbNCRKKAfUB9LSGo8iaCxYsXM3z4cI4ePVo4Tlq0aAF4VpUqvPh2ZGQkzZo1A2D//v3nTMgWFRVFYmIiqsq+ffvOWfy9WrVqNGnSBIDdu3eTnZ1NREQEDRo0oHr16sDPN40LKu+2QNZV1Lbc3FxycnKoVq0aubm5ZGZmnvPHV6NGDWJjY8nJycn/wy1YT3R0NHl5eZw5c+ac85VERBAR6tatS40aNcjOzubQoUPnxNqkSRMiIiI4cuRI/u+z4Pnj4+M5fvz4Ob9LESEuLo7HHnuMYcOGMW3aNObNm3fOe/MNiho/fjwffvjhWfsbNGjAsmXLABgzZsw5i/wkJSXlD9YbOXLkOaNt27dvn99td9CgQWzZsuWs/d26dePf//43AP369ePHH388a3+vXr14/vnnAbjyyivJyMg4a//111/Pk08+CXimjC78/3Xw4MFMnjwZVaVDhw4UNmrUKMaPH8/x48e59NJLz9n/u9/9jjFjxrB//36uuuqqc/ZPmDCBkSNH8sMPPxT5JWjy5MkMHjyYb775hltvPfc74V133UV2djbz58/n66+/Puf/Z7Vq1ahRowYREREcPXr0nP+bjRo1IjY2luPHj5/zswFo2rQpMTExHD169Jz9IkJiYiLVqlXj8OHDZGZmnnN88+bNiYqKIisr66z/mz4tWrQgIiKCjIyM/L+Ngi644AKg6M+imJgYJk2aVO41oUtKBFHlqtE/TYFdBV6nA4X/5+SXUdUcETkCxANn/QZEZDQwGsj/UC6rWrVq0bx5c/bv31+wXiIiIvKnX05LS8v/5fv+A0VFRdG5c2dEhK1bt+ZPA+HbHxMTQ+fOnQH49ttvz0k01atXz69/y5YtxMTE0LRpU2JjY88qV9Q3mfJuC2RdJW07efIkq1evJjc396w+4RdeeCGJiYkcP36cVatWnXV8REQEnTp1omXLlpw+fZqtW7dywQUXkJiYSFxcHCdPnqR69er5V09r164lNzc3/5GXl0e3bt2IjIxk+/bt7Nix46z6RYRu3boRExPDzp07z/qg9JXr3bs3derUYfv27fzwww/UrFmT8847jxo1aiAi+X9ojRs3pl27dme974JXQU2bNj1nf926dfOfJyYmnrPf96UAPEnhxIkTZ+33fSkBaNmy5Tk/++bNm+c/b9WqFTVr1jxrf8G/jzZt2pzTG61p06b5zy+66CKys7PP2t+4ceP854VjB2jYsCHg+T0Wtd83piIqKqrI/b6fT7Vq1Yrc71uvOy4ursj93bp145e//CUDBgzgkUce4ejRo5w4cYKTJ09y6tQpkpOTqV27Nrt27SoyUbRv357atWtz4MABtm7dek79F110ETVr1mTfvn1nXYX66mnbti1xcXHs2bOnyL+LCy+8kGrVqrFr164i1wC/6KKLiIyMZMeOHezdu/ec/b73HBMTw4EDB87aV7Nmzfyff6A5eUUwCLhWVX/jfT0CuFRVxxYo8423TLr39TZvmXNTtVcoNw0ZY4xb3Jp0bjeQWOB1gndbkWW8TUO1gXOvt4wxxjjGyUSwBmgtIi1EpBpwG7CgUJkFwEjv80HA0pLuDxhjjAk8x+4ReNv8xwKLgUjgVVXdJCKP4unGtAD4F/C6iKQBWXiShTHGmCBy8mYxqvoh8GGhbY8UeJ4NDHYyBmOMMSWzuYaMMSbMWSIwxpgwZ4nAGGPCnCUCY4wJc44NKHOKiBwEyjs5fj0KjVoOA/aew4O95/BQkffcXFWLXE4v5BJBRYjI2uJG1lVV9p7Dg73n8ODUe7amIWOMCXOWCIwxJsyFWyJ42e0AXGDvOTzYew4PjrznsLpHYIwx5lzhdkVgjDGmEEsExhgT5sImEYjItSLynYikiciDbsfjNBF5VUQOeBf/CQsikigiy0Rks4hsEpH73Y7JaSISKyKrRWSD9z3/xe2YgkFEIkXkaxFZ6HYswSAiO0Rko4isF5GAr8wVFvcIRCQS2Apcg2fJzDXAEFXd7GpgDhKRK4DjwExVPXfx2SpIRBoDjVX1KxE5D1gH3FjFf88C1FDV4yISDXwO3K+qK10OzVEi8kcgGailqucuflzFiMgOILmk1RsrIlyuCLoBaaq6XVVPA28BA12OyVGq+hmeNR7ChqruVdWvvM+PAVvwrItdZanHce/LaO+jSn+7E5EE4HpgutuxVBXhkgiaArsKvE6nin9AhDsRSQI6A6tcDsVx3maS9cABYImqVvX3/DQwAchzOY5gUuA/IrJOREYHuvJwSQQmjIhITWA+8ICqHnU7Hqepaq6qdsKzLng3EamyTYEi0h84oKrr3I4lyH6pql2AfsDvvE2/ARMuiWA3kFjgdYJ3m6livO3k84HZqpridjzBpKqHgWXAtS6H4qSewABvm/lbwFUiMsvdkJynqru9/x4A3sHT3B0w4ZII1gCtRaSFiFTDszbyApdjMgHmvXH6L2CLqk51O55gEJH6InK+93kcng4R37oalINU9SFVTVDVJDx/x0tVdbjLYTlKRGp4Oz8gIjWAvkBAewOGRSJQ1RxgLLAYzw3EOaq6yd2onCUibwIrgLYiki4id7kdUxD0BEbg+Za43vu4zu2gHNYYWCYiqXi+8CxR1bDoUhlGGgKfi8gGYDXwgaouCuQJwqL7qDHGmOKFxRWBMcaY4lkiMMaYMGeJwBhjwpwlAmOMCXOWCIwxJsxZIjCmgkRklIg0ceGc/xfMc5qqyxKBCQsiEuVg9aOAMiUCh+MxpkwsEZiQICJJIvKtiMwWkS0iMk9Eqnv3PSIia0TkGxF52TvCGBH5VESe9s7ffr+I3CAiq7zz2H8sIg295SaLyAwR+a+I7BSRX4vIk9753xd5p61ARLqKyHLvxF+LRaSxiAzCMx3ybO8AtriiyhUVT4H3FuGdb/78Atu+F5GGxcVc6GfzmjcO3+vjBZ6P9/5sUsNlrQJTdpYITChpC0xT1YuAo8AY7/b/U9VLvOsuxAEF56evpqrJqvoUnrn6u6tqZzzz1EwoUO4C4CpgADALWKaqFwMngeu9yeA5YJCqdgVeBf6qqvOAtcAw78RvOUWVKyYeAFQ1D3gPuAlARC4Fdqrq/lJiLpGI9AVa45mXphPQNdCTlZmqwS5PTSjZpapfeJ/PAu4D/gH0FpEJQHWgLrAJeN9b7u0CxycAb3u/oVcDfiiw7yNVPSMiG4FIwDeEfyOQhCcJdQCWeC84IoG9RcRYWrm3izjGt/0R4N945tDxlSsp5tL09T6+9r6uiScxfFaGOkwYsERgQknh+VBURGKBaXhWb9olIpOB2AJlfirw/DlgqqouEJFewOQC+06B59u5iJzRn+deycPzdyLAJlXtUUqMpZX7qZjtK4BWIlIfuBF43I+YfXLwXt2LSASehOGL5W+q+lIpMZswZ01DJpQ0ExHfB+xQPM0mvg/9DO86BIOKPNKjNj9PPz6yjOf+DqjvO7+IRItIe+++Y8B5fpQrljfxvANMxTN7amYZYt4BdPU+H4BnlTLwTLJ4p/fngog0FZEGpcViwo8lAhNKvsOzKMcWoA7wgncO/lfwTMu7GM8MnMWZDMwVkXVAmdZ+9S5xOgj4u3cWyPXAZd7drwEvimeVsMgSypXmbWA4Zzcf+RPzK8CV3vP1wHvVoar/Ad4AVnibvObxc8IyJp/NPmpCgniWnlzovSFsjAkguyIwxpgwZ1cExhgT5uyKwBhjwpwlAmOMCXOWCIwxJsxZIjDGmDBnicAYY8Lc/wfkRQgNNQwhYQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "space = np.linspace(0, 5, mplot)\n",
    "plt.plot(space, kde(sample[:, 1])(space), color='black', label='lambda1')\n",
    "plt.plot(space, kde(sample[:, 2])(space), color='black', linestyle='dashed', label='lambda2')\n",
    "plt.xlabel('parameter value')\n",
    "plt.ylabel('posterior density')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3920c01",
   "metadata": {},
   "source": [
    "#### (b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e31c518a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgyElEQVR4nO3df5RdZX3v8fcnmSETQiY/mEBCJkOghAByxR9DIBVZxtt6QQlor4qiVnutdFW9oVjaattFW9vVq65ruS1oWfFnbRFai7agKHKBJXgNhiRAIISELCBmkkBmksyQ35lJvvePs088GSczZ2b2j5lzPq+1zso5++yz9/esrJnPPM/ez/MoIjAzs/o1oegCzMysWA4CM7M65yAwM6tzDgIzszrnIDAzq3MNRRcwXC0tLTF//vyiyzAzG1dWr17dFRGzBnpv3AXB/PnzWbVqVdFlmJmNK5I2n+g9dw2ZmdU5B4GZWZ1zEJiZ1TkHgZlZnXMQmJnVOQeBmVmdcxCYmdU5B4EV7oc//CFr1qwpugyzujXuBpRZbTl69CjXXXcdkli5ciXnnHNO0SWZ1R23CKxQzzzzDN3d3ezevZurr76agwcPFl2SWd1xEFihHn30UQA+97nPsX79elauXFlwRWb1x0FghXr00UdpbW3l2muvBeC5554ruCKz+uMgsMJEBI8++ihvfvObaWtrY/LkyQ4CswI4CKwwL774Itu2bePNb34zEyZMYOHChaxfv77osszqjoPACrNixQoALrvsMgDOP/98twjMCuAgsMK8+OKLACxYsACA8847j82bN7N///4iyzKrOw4CK0xHRwctLS00NTUBpSCICDZu3FhwZWb1xUFghdm6dSutra3HXp9//vmA7xwyy5uDwArT0dFxXBAsWLCACRMmOAjMcuYgsMJ0dHQwd+7cY6+bmpo466yzHARmOXMQWCEOHjxIV1fXcS0CgDPPPJOOjo6CqjKrT5kFgaQmSSslPSVpnaS/GmCfSZL+VdImST+XND+remxs2bZtG8BxLQKA2bNn8/LLLxdRklndyrJFcAh4a0RcBLwOuELSpf32+SiwOyLOAW4BPp9hPTaGlP/q798imD17Ntu3byciiijLrC5lFgRRsjd52Zg8+v90XwP8U/L834H/KklZ1WRjx9atW4GBg2D//v3s3bt3oI+ZWQYyvUYgaaKkJ4EdwAMR8fN+u8wFtgBERB/QA5w6wHGul7RK0qrOzs4sS7aclFsEA3UNAe4eMstRpkEQEUci4nVAK7BI0oUjPM7yiGiPiPZZs2alWqMVo6Ojg6lTp9Lc3Hzc9jlz5gAOArM85XLXUER0Aw8DV/R7ayswD0BSAzAN2JlHTVasrVu3/kprANwiMCtClncNzZI0PXk+GfhNoP8N4vcAH06evxt4KHyVsC70H0xWVg6C7du3512SWd3KskUwB3hY0lrgcUrXCL4v6bOSrk72+RpwqqRNwKeAT2dYj40h/QeTlc2cOZOGhga3CMxylNni9RGxFnj9ANtvrnh+EHhPVjXY2BQRdHZ2cvrpp//KexMmTPBYArOceWSx5W7fvn0cPnyYU0/9lRvEAA8qM8ubg8By19XVBeAgMBsjHASWu507SzeGtbS0DPh+eXSxmeXDQWC5KwfBYC2CHTt2cOTIkTzLMqtbDgLL3VBBMGfOHI4ePXqsC8nMsuUgsNyVf8EP1jUEHktglhcHgeWu3CKYMWPGgO+XpxHxvFJm+XAQWO527tzJ9OnTaWgYeBhLucuoHBhmli0HgeVu586dJ7w+AA4Cs7w5CCx3XV1dJ7w+AKVpJgB27dqVV0lmdc1BYLkbqkXQ2NhIc3OzWwRmOXEQWO6GCgIodQ85CMzy4SCw3DkIzMYWB4Hl6tChQ+zdu3fQawTgIDDLk4PAcjXUqOIyB4FZfhwElqtqg2DmzJkOArOcOAgsV0PNPFp26qmn0t3dTV9fXx5lmdU1B4Hlaqi1CMrK7+/evTvzmszqnYPAcjWcawSV+5tZdhwElqvyX/jl0cMn4iAwy4+DwHLV3d3NSSedxOTJkwfdz0Fglp/MgkDSPEkPS3pW0jpJNwywz1sk9Uh6MnncnFU9NjZ0d3czffr0IfdzEJjlZ+B5gNPRB/xhRKyRNBVYLemBiHi2336PRsRVGdZhY0hPTw/Tpk0bcr9yEHjiObPsZdYiiIjtEbEmeb4HWA/Mzep8Nj5U2yKYOnUqDQ0NbhGY5SCXawSS5gOvB34+wNuLJT0l6YeSXnOCz18vaZWkVV61anyrNggkeVCZWU4yDwJJpwB3A38QEa/2e3sNcGZEXATcCvzHQMeIiOUR0R4R7eVlDG18qrZrCDzNhFleMg0CSY2UQuCOiPhu//cj4tWI2Js8vw9olDT4kFMb16ptEYCDwCwvWd41JOBrwPqI+LsT7DM72Q9Ji5J6/JNfw4YTBDNnzvTFYrMcZHnX0JuADwFPS3oy2fanQBtARNwOvBv4fUl9wAHgfRERGdZkBTp8+DAHDhyoumtoxowZPPHEExlXZWaZBUFE/BTQEPvcBtyWVQ02tvT09ABU3SKYMWOG5xoyy4FHFltuuru7geEFwd69e+nt7c2uKDNzEFh+yi2CaruGyvMRlQPEzLLhILDcjKRFAJ6K2ixrDgLLzXBbBOUg8J1DZtlyEFhu3CIwG5scBJYbB4HZ2DRkEEgafCkpsyr19PQgiVNOOaWq/csXix0EZtmqpkXwmKTvSHp7eRSw2Uh0d3czbdo0JkyoriHqFoFZPqr5iTwXWE5plPDzkv5W0rnZlmW1aDjTSwA0NjYyZcoUXyw2y9iQQRAlD0TE+4GPAR8GVkr6iaTFmVdoNWM4M4+WeXSxWfaGnGIiuUbwQUotgleA/wncA7wO+A5wVob1WQ0ZbosAHARmeahmrqEVwD8D74yIjortqyTdnk1ZVot6enqYP3/+sD7jIDDLXjXXCP48Iv66MgQkvQcgIj6fWWVWc8oXi4dj5syZDgKzjFUTBJ8eYNtn0i7Eap+7hszGphN2DUm6Eng7MFfSP1S81Qz0ZV2Y1ZajR4/y6quvjuhise8aMsvWYNcItgGrgKuB1RXb9wA3ZlmU1Z49e/YQESNqEezfv5/Dhw9z0kknZVOcWZ07YRBExFPAU5LuiAi3AGxUhju9RFnloLLTTz895arMDAbvGvq3iHgv8ISkyuUjRWl4wWszr85qxnBnHi2rnGbCQWCWjcG6hm5I/r0qj0KstqXRIjCzbJzwrqGI2J487QK2RMRmYBJwEaXrB2ZVG20Q+IKxWXaquX30EaBJ0lzgx5RGGH8zy6Ks9oy0a8gtArPsVRMEioj9wG8BX46I9wCvGfJD0jxJD0t6VtI6STcMsI8k/YOkTZLWSnrD8L+CjQfuGjIbu6oKgmRyuQ8AP0i2Tazic33AH0bEBcClwCckXdBvnyuBBcnjeuAfq6raxp2RtgjKweEgMMtONUFwA6WRxN+LiHWSzgYeHupDEbE9ItYkz/cA64G5/Xa7BvhWMsPpY8B0SXOG9Q1sXOju7ubkk0+msbFxWJ9rbGxk6tSpDgKzDA056VxEPELpOkH59QvAsuGcRNJ84PXAz/u9NRfYUvG6I9m2vXInSddTajHQ1tY2nFPbGDGS6SXKPLrYLFvVTEN9LnATML9y/4h4azUnkHQKcDfwBxHx6kiKjIjllBbHob29PYbY3cagkaxFUOb5hsyyVc001N8Bbge+ChwZzsElNVIKgTsi4rsD7LIVmFfxujXZZjVmtC0CB4FZdqoJgr6IGPZF3GR9468B6yPi706w2z3AJyXdBVwC9FSMX7Aa0t3dzaxZs0b02RkzZrBx48aUKzKzsmqC4F5JHwe+Bxwqb4yIoTpt30RpzMHTkp5Mtv0p0JZ8/nbgPkoznG4C9gO/M5zibfzo6enhnHPOGdFnvSaBWbaqCYIPJ//+UcW2AM4e7EMR8VNK8xINtk8An6iiBhvn3DVkNnZVc9eQ1yS2UYmIUQfBgQMHOHjwIE1NTekWZ2ZDjyOQdLKkP5e0PHm9QJInorOqHTx4kN7e3lHdNQQeVGaWlWoGlH0DOAz8evJ6K/A3mVVkNWek00uUOQjMslVNEPxaRHwB6AVI5h0atO/frNJIp5cocxCYZauaIDgsaTKlC8RI+jUq7h4yG8poWwSVi9OYWfqquWvoL4AfAfMk3UHpttCPZFmU1Za0uoY8zYRZNqq5a+gBSWsozSAq4IaI6Mq8MqsZ7hoyG9sGW7O4/9oA5RG/bZLayjOLmg1ltC0CT0Vtlq3BWgRfTP5tAtqBpyi1CF4LrAIWZ1ua1YrRBsHEiRNpbm52EJhlZLA1i5dExBJKLYE3RER7RLyR0nTSnhjOqtbT00NDQwOTJ08e8TE8zYRZdqq5a2hhRDxdfhERzwDnZ1eS1ZryqOLSPIQj42kmzLJTzV1DayV9FfiX5PUHgLXZlWS1ZjTTS5R5cRqz7FTTIvgdYB2lJStvAJ7Fs4TaMIxmUZoytwjMslPN7aMHgVuSh9mwpdUicBCYZaOaFoHZqLhFYDa2OQgsc2m0CGbOnMmhQ4c4cOBAOkWZ2TGDBoGkiZL+d17FWG1Kq2sIPM2EWRYGDYKIOAJcllMtVoP6+vrYt29fKl1D4NHFZlmo5vbRJyTdA3wH2FfeGBHfzawqqxnleYbSahE4CMzSV00QNAE7gbdWbAvAQWBDGu30EmUOArPsVHP7qMcM2IiNdubRMq9JYJadatYsbpX0PUk7ksfdklrzKM7GP7cIzMa+atcsvgc4I3ncm2wblKSvJ8HxzAnef4ukHklPJo+bh1O4jQ9pXSOYNm0aknzXkFkGqgmCWRHxjYjoSx7fBGZV8blvAlcMsc+jEfG65PHZKo5p40y5RTDarqEJEyYwbdo0twjMMlBNEOyU9MFkTMFESR+kdPF4UBHxCOA/3+pcWl1D4NHFZlmpJgj+B/Be4GVKaxO8m/QmnVss6SlJP5T0mhPtJOl6Saskrers7Ezp1JaHctfQ1KlTR30sr0lglo1q7hraDFydwbnXAGdGxF5Jbwf+A1hwghqWA8sB2tvbI4NaLCPd3d00NzczceLEUR/LLQKzbAy2ZvEfR8QXJN1KadzAcSJi2WhOHBGvVjy/T9KXJbVERNdojmtjSxrTS5TNmDGDLVu2pHIsM/ulwVoE65N/V2VxYkmzgVciIiQtotRNNeS1Bxtf0ph5tMwtArNsnDAIIuJeSROB/xIRNw33wJLuBN4CtEjqAP4CaEyOfTulaw2/L6kPOAC8LyLc7VNj0m4R7N69m4gY1bKXZna8Qa8RRMQRSW8ayYEj4v1DvH8bcNtIjm3jR3d3N/PmzUvlWDNmzKC3t5f9+/czZcqUVI5pZtXNNfSkJ52zkerp6eHCCy9M5ViV00w4CMzS40nnLFNpdw1BKQhaWz3LiVlaPOmcZSYi6OnpST0IPM2EWbqqmXTuXEkPlucMkvRaSX+efWk23u3Zs4ejR4+metcQeOI5s7RVM7L4K8BngF6AiFgLvC/Loqw2lH9hl/v2R+vUU08FYOdO32VslqZqguDkiFjZb1tfFsVYbSl34aQVBC0tLYCDwCxt1QRBl6RfIxldLOndlOYcMhtUOQjKXTqjdfLJJ9PU1ERXlwefm6WpmruGPkFpnp/zJG0FXgQ+kGlVVhPS7hqSREtLi4PALGXVBEFExG9ImgJMiIg9ks7KujAb/9LuGgIcBGYZqKZr6G6AiNgXEXuSbf+eXUlWK8otgrS6hqAUBJ6K3Cxdg80+eh7wGmCapN+qeKuZ0iAzs0Ht2rWLSZMmMXny5NSO2dLSwksvvZTa8cxs8K6hhcBVwHRgacX2PcDHMqzJasSuXbuYOXNmqhPEuWvILH2DzT76n8B/SlocEStyrMlqxO7du1PtFoJSEHR3d9Pb20tjY2OqxzarV9VcI3iXpGZJjckI485k3WKzQZVbBGmaNWvWsWObWTqqCYK3JauJXQW8BJwD/FGWRVlt2LVrVyYtAsDdQ2YpqiYIyu3vdwDfiYieDOuxGrJ79+7UWwQOArP0VTOO4F5Jz1FaRez3Jc0CDmZbltWCLLqGHARm6RuyRRARnwZ+HWiPiF5Ki9Nck3VhNr719vayd+9edw2ZjQNDtggkNQIfBC5PbgP8CXB7xnXZOJf29BJl5RlIPajMLD3VdA39I6XrBF9OXn8o2fa7WRVl41/aE86VTZo0ialTp7pFYJaiaoLg4oi4qOL1Q5Keyqogqw1ZtQjAg8rM0lbNXUNHkmmoAZB0NnBkqA9J+rqkHeWVzQZ4X5L+QdImSWslvaH6sm2sy2LCuTIHgVm6qmkR/BHwsKQXAAFnAtWsY/xN4DbgWyd4/0pgQfK4hFJ30yVVHNfGgay6hqAUBDt27Ej9uGb1qprF6x+UtIDS3EMAGyLiUBWfe0TS/EF2uQb4VkQE8Jik6ZLmRIQXvakBWXcNrVu3LvXjmtWrau4aagI+DlxGaZWyRyXdHhGjHUswF9hS8boj2fYrQSDpeuB6gLa2tlGe1vJQbhFMnz499WPPnj2bV155hYhIdUI7s3pVzTWCb1GajvpWSl09rwH+Ocui+ouI5RHRHhHt5blmbGzbtWsX06ZNY+LEiakfe/bs2Rw6dIieHg9yN0tDNdcILoyICypePyzp2RTOvRWYV/G6NdlmNaCzs5OsQnv27NkAvPzyy5m0OMzqTTUtgjWSLi2/kHQJsCqFc98D/HZy99ClQI+vD9SOHTt25BIEZjZ61bQI3gj8TNIvktdtwAZJT1Naz/i1A31I0p3AW4AWSR3AX5BMYBcRtwP3AW8HNgH7qe5OJBsnOjs7OfvsszM5toPALF3VBMEVIzlwRLx/iPcD+MRIjm1jX2dnJ5dcks3dwA4Cs3RVc/vo5jwKsdpx9OhRurq6MusamjFjBo2NjQ4Cs5RUc43AbFi6u7vp6+vLLAgkMXv2bAeBWUocBJa68sygp512WmbncBCYpcdBYKkrB0GWYz4cBGbpcRBY6hwEZuOLg8BSV54QLuuuoc7OTo4cGXIiXDMbgoPAUlduEZSXlczC7NmzOXr0qFcqM0uBg8BS19nZSXNzM5MmTcrsHB5LYJYeB4GlLst5hsocBGbpcRBY6nbs2JHp9QFwEJilyUFgqcujRTBnzhwAOjo6Mj2PWT1wEFjq8giCyZMnc9ppp/GLX/xi6J3NbFAOAktVRGQ6z1CltrY2B4FZChwElqqenh56e3sdBGbjiIPAUrVt2zbgl334WSoHQWlGczMbKQeBpWrr1tJqo62trZmfq62tjX379rF79+7Mz2VWyxwElqryXTxz587N/FxtbW0A7h4yGyUHgaWq3CI444wzMj+Xg8AsHQ4CS1VHRwctLS00NTVlfi4HgVk6HASWqq1bt+ZyfQBK01xPmjTJQWA2Sg4CS1VHR0cu1wcAJkyYwLx58xwEZqOUaRBIukLSBkmbJH16gPc/IqlT0pPJ43ezrMey19HRkVuLADyWwCwNmQWBpInAl4ArgQuA90u6YIBd/zUiXpc8vppVPZa9gwcP0tXVlXsQbN68ObfzmdWiLFsEi4BNEfFCRBwG7gKuyfB8VrDyYLK8uoYAzj77bLZt28b+/ftzO6dZrckyCOYCWypedyTb+vvvktZK+ndJ8zKsxzKW52CysoULFwLw/PPP53ZOs1pT9MXie4H5EfFa4AHgnwbaSdL1klZJWuWlCceuPAeTlZWDYMOGDbmd06zWZBkEW4HKv/Bbk23HRMTOiDiUvPwq8MaBDhQRyyOiPSLa85jMzEamiBbBggULkMRzzz2X2znNak2WQfA4sEDSWZJOAt4H3FO5g6TKmcmuBtZnWI9lrKOjg1NOOYXm5ubcznnyySfT1tbmFoHZKDRkdeCI6JP0SeB+YCLw9YhYJ+mzwKqIuAdYJulqoA/YBXwkq3osey+99NKx0b55WrhwoYPAbBQyCwKAiLgPuK/ftpsrnn8G+EyWNVh+Nm7cyHnnnZf7eRcuXMg3vvENIgJJuZ/fbLwr+mKx1Yi+vj42bdrEueeem/u5Fy5cyN69e9m+fXvu5zarBQ4CS8XmzZvp7e0tLAjAdw6ZjZSDwFKxceNGgEKDwHcOmY2Mg8BSUQ6C8i/lPLW2tjJjxgyeeOKJ3M9tVgscBJaKjRs3Mn36dFpaWnI/tyTa29t5/PHHcz+3WS1wEFgqNmzYwLnnnlvYXTuLFi3i6aef5sCBA4Wc32w8cxBYKjZu3FjI9YGyiy++mCNHjrh7yGwEHAQ2avv372fLli2FBsGiRYsA3D1kNgIOAhu19etLM4MUMZisbM6cOcydO9dBYDYCDgIbtcceewz45V/lRbn44otZuXJloTWYjUcOAhu1FStWMGfOnELmGaq0ePFinn/++WML5JhZdRwENmorVqxg8eLFhc/zc+WVVwJw3333DbGnmVVyENio7NixgxdeeIFLL7206FK48MILaWtr4wc/+EHRpZiNKw4CG5Xy9YHFixcXXElpYNk73vEOHnjgAQ4dOjT0B8wMcBDYKD322GM0NDTwxjcOuLhc7q666ir27dvHT37yk6JLMRs3HAQ2Kj/+8Y9pb29n8uTJRZcCwJIlS5gyZQp33XVX0aWYjRsOAhuxDRs2sHr1at773vcWXcoxkydP5kMf+hDf/va36ezsLLocs3HBQWAjdueddyKJa6+9tuhSjrNs2TIOHTrE8uXLiy7FbFxwENiIRAR33HEHS5Ys4Ywzzii6nOOcf/75vO1tb+NLX/oS+/fvL7ocszHPQWAj8tBDD7Fp0yauu+66oksZ0J/92Z+xfft2br755qF3NqtzDgIbtt7eXpYtW8b8+fPHbBBcfvnl/N7v/R633HILP/vZz4oux2xMcxDYsN1yyy08++yz/P3f//2YuVtoIF/4whdoa2tj6dKlrF69uuhyzMasTINA0hWSNkjaJOnTA7w/SdK/Ju//XNL8LOux0bv11lv5kz/5E975zneydOnSossZVHNzMw8++CDNzc0sWbKEW2+9lb6+vqLLMhtzMgsCSROBLwFXAhcA75d0Qb/dPgrsjohzgFuAz2dVj43MkSNH2Lx5M9/+9re5/PLLWbZsGe9617uO3TE01p199tk88sgjXHrppSxbtowzzzyTG2+8kbvvvpt169bR1dXF0aNHiy7TrFCKiGwOLC0G/jIi/lvy+jMAEfG/Kva5P9lnhaQG4GVgVgxSVHt7e6xatWrY9dx///3ceOONJDUc236i59Xul+fzIs69Z88ejhw5AsC8efO46aab+PjHP05DQwPjSUTw/e9/n6985Sv86Ec/ore399h7EydOZOrUqTQ0NNDY2EhDQ8Oxx4QJE3INvLzDdTyEuf3SRz/6UT71qU+N6LOSVkdE+0DvZfnTPBfYUvG6A7jkRPtERJ+kHuBUoKtyJ0nXA9cDI57quLm5mQsvvLDymEM+r3a/rJ8Xee7m5mZaW1u5+OKLueiii8ZdAJRJYunSpSxdupSDBw+ydu1aXnzxRV555RVeeeUV9uzZQ19f3688yiGYh6z+KBsr57PRO/300zM57rj4qY6I5cByKLUIRnKMxYsXj4mJ0ax4TU1NLFq0qPCFdMzGiiwvFm8F5lW8bk22DbhP0jU0DdiZYU1mZtZPlkHwOLBA0lmSTgLeB9zTb597gA8nz98NPDTY9QEzM0tfZl1DSZ//J4H7gYnA1yNinaTPAqsi4h7ga8A/S9oE7KIUFmZmlqNMrxFExH3Aff223Vzx/CDwnixrMDOzwXlksZlZnXMQmJnVOQeBmVmdcxCYmdW5zKaYyIqkTmDzCD/eQr9Ry3XA37k++DvXh9F85zMjYtZAb4y7IBgNSatONNdGrfJ3rg/+zvUhq+/sriEzszrnIDAzq3P1FgTLiy6gAP7O9cHfuT5k8p3r6hqBmZn9qnprEZiZWT8OAjOzOlc3QSDpCkkbJG2S9Omi68mapK9L2iHpmaJryYukeZIelvSspHWSbii6pqxJapK0UtJTyXf+q6JryoOkiZKekPT9omvJg6SXJD0t6UlJw1+rd6jj18M1AkkTgY3Ab1JaMvNx4P0R8WyhhWVI0uXAXuBbEXHhUPvXAklzgDkRsUbSVGA18M4a/38WMCUi9kpqBH4K3BARjxVcWqYkfQpoB5oj4qqi68mapJeA9ojIZABdvbQIFgGbIuKFiDgM3AVcU3BNmYqIRyit8VA3ImJ7RKxJnu8B1lNaF7tmRcne5GVj8qjpv+4ktQLvAL5adC21ol6CYC6wpeJ1BzX+C6LeSZoPvB74ecGlZC7pJnkS2AE8EBG1/p3/D/DHwNGC68hTAD+WtFrS9WkfvF6CwOqIpFOAu4E/iIhXi64naxFxJCJeR2ld8EWSarYrUNJVwI6IWF10LTm7LCLeAFwJfCLp+k1NvQTBVmBexevWZJvVmKSf/G7gjoj4btH15CkiuoGHgSsKLiVLbwKuTvrM7wLeKulfii0pexGxNfl3B/A9St3dqamXIHgcWCDpLEknUVob+Z6Ca7KUJRdOvwasj4i/K7qePEiaJWl68nwypRsiniu0qAxFxGciojUi5lP6OX4oIj5YcFmZkjQlufkBSVOAtwGp3g1YF0EQEX3AJ4H7KV1A/LeIWFdsVdmSdCewAlgoqUPSR4uuKQdvAj5E6a/EJ5PH24suKmNzgIclraX0B88DEVEXt1TWkdOBn0p6ClgJ/CAifpTmCeri9lEzMzuxumgRmJnZiTkIzMzqnIPAzKzOOQjMzOqcg8DMrM45CMwASXuH3quq4/ylpJuq2O+bkt6dxjnNRstBYGZW5xwEZhUknSLpQUlrkvnfr0m2z5f0XPKX/EZJd0j6DUn/T9LzkiqH/F8kaUWy/WPJ5yXptmRNjP8LnFZxzpslPS7pGUnLkxHSZrlxEJgd7yDwrmSCryXAFyt+MZ8DfBE4L3lcB1wG3AT8acUxXgu8FVgM3CzpDOBdwELgAuC3gV+v2P+2iLg4WTdiMlDz8+vb2NJQdAFmY4yAv01mdzxKabry05P3XoyIpwEkrQMejIiQ9DQwv+IY/xkRB4ADkh6mNEHY5cCdEXEE2CbpoYr9l0j6Y+BkYCawDrg3s29o1o+DwOx4HwBmAW+MiN5klsum5L1DFfsdrXh9lON/lvrP23LCeVwkNQFfprT61BZJf1lxPrNcuGvI7HjTKM133ytpCXDmCI5xTbKW8KnAWyhNBvcIcG2yiMwcSt1O8Mtf+l3JOgq+k8hy5xaB2fHuAO5NuntWMbIpnddSWhegBfjriNgm6XuUrhs8C/yC0sywRES3pK9Qmlb4ZUqhYZYrzz5qZlbn3DVkZlbnHARmZnXOQWBmVuccBGZmdc5BYGZW5xwEZmZ1zkFgZlbn/j83HLMTWS7ESQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "space = np.linspace(0, 5, mplot)\n",
    "plt.plot(space, gamma(a=(alpha + np.sum(y)), scale=1/(beta + n)).pdf(space), color='black')\n",
    "plt.xlabel('lambda')\n",
    "plt.ylabel('posterior density')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}