fehiepsi/changepoint.ipynb

## changepoint.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax\n",
    "import jax.config as config; config.update(\"jax_platform_name\", \"cpu\")\n",
    "import jax.numpy as np\n",
    "\n",
    "import numpyro\n",
    "import numpyro.distributions as dist\n",
    "from numpyro.mcmc import MCMC, NUTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 2000/2000 [00:06<00:00, 331.75it/s, 1023 steps of size 1.84e-03. acc. prob=0.79]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "                mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "   lambda1     17.76      0.59     17.73     16.79     18.74    195.84      1.00\n",
      "   lambda2     22.71      0.85     22.65     21.43     24.16    156.59      1.00\n",
      "       tau      0.59      0.01      0.59      0.58      0.61    149.88      1.00\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def model(data):\n",
    "    alpha = 1 / np.mean(data)\n",
    "    lambda1 = numpyro.sample('lambda1', dist.Exponential(alpha))\n",
    "    lambda2 = numpyro.sample('lambda2', dist.Exponential(alpha))\n",
    "    tau = numpyro.sample('tau', dist.Uniform(0, 1))\n",
    "    lambda12 = np.where(np.arange(len(data)) < tau * len(data), lambda1, lambda2)\n",
    "    numpyro.sample('obs', dist.Poisson(lambda12), obs=data)\n",
    "\n",
    "count_data = np.array([\n",
    "    13,  24,   8,  24,   7,  35,  14,  11,  15,  11,  22,  22,  11,  57,\n",
    "    11,  19,  29,   6,  19,  12,  22,  12,  18,  72,  32,   9,   7,  13,\n",
    "    19,  23,  27,  20,   6,  17,  13,  10,  14,   6,  16,  15,   7,   2,\n",
    "    15,  15,  19,  70,  49,   7,  53,  22,  21,  31,  19,  11,  18,  20,\n",
    "    12,  35,  17,  23,  17,   4,   2,  31,  30,  13,  27,   0,  39,  37,\n",
    "    5,  14,  13,  22,\n",
    "])\n",
    "mcmc = MCMC(NUTS(model), num_warmup=1000, num_samples=1000)\n",
    "mcmc.run(jax.random.PRNGKey(0), count_data)\n",
    "mcmc.print_summary()\n",
    "samples = mcmc.get_samples()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABdoAAAFTCAYAAAAugr49AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd5hU5d0G4OdQlCaCgkZBRSUWsGDH3ktiQZNoSLFEjSVqoqLGEgUVo7FHjbHXJPYoxq6xxBqNisbeQBRbVESUJjDfHwv7sVTh7DIs3Pd17cXOe8r85gDznnnmPe8pKpVKAAAAAACA2dOk2gUAAAAAAEBjJmgHAAAAAIASBO0AAAAAAFCCoB0AAAAAAEoQtAMAAAAAQAmCdgAAAAAAKEHQDvWkKIp+RVFUqvGcRVE0K7P9bG57eFEU/yiK4sOJNfSbnf0AwJw0P/XXRVGsUBTFH4uieLEoiq8m9tm3F0Wx+uzUAQDzs6Iodi6K4vBq1wHMvQTtwOz6ZZLFktxW7UIAgGnaJsnmSa5OsmOSXyXpmOTfRVGsVc3CAKAR2jmJoB2YrtkaVQOQpHulUpkwcXTeAdUuBgCYyvVJ/lSpVGpHwxdF8WCSwUl+k2SPKtUFAADzHCPaoYEURXFwURRPFkXxeVEUXxRF8VRRFNtPsU6XiZeSH1AUxalFUXxUFMWIoij+UhRFq6IouhZFce/Ey73fKopiz+k83cpFUTxUFMXIiZeFn1QURZMpnmuNoigeLYpidFEUQ4uiOD5JMTt1J0mlUplQ5vgAwNxgXu6vK5XKp5OH7BPbhid5I0mn2TleADA/KoriqiR7Juk08ZygUhTF4KIoWhRFcU5RFC9NPA/4aOIUqytNsf00p4EriuKqoigGz5lXATQ0I9qh4XRJcllqRo01S80l23cURfH9SqVy9xTrHpPk4dR03N2SnJ5kQpI1klya5MwkBya5siiK/1QqlZen2P62JFckOTXJtkmOn7h9vyQpiqJDkgeTfDTxOcYkOTLJ0iXrBoDGrkvmo/66KIpFkqyS5MrprQMATOXk1Ey/tk6SnSa2jUmyYJKFkvRP8mGSRVIzVdtTRVGsVKlUPqpCrUCVCNqhgVQqlSMm/T5xtNo/k6yQmmlWpvwA/HalUpk0+u3eoig2TrJ7kt0rlcpfJu7jP6np0H+UZMoP7pdWKpXTJv5+X1EUbZP0KYri3Eql8kWSw5K0TrJtpVIZMnF/9yd5t2TdANCozYf99fmpGSF/7gzWAQAmU6lU3i6K4n9JxlYqlaemWLzvpF+Komia5N4kHyf5SZJz5lyVQLWZOgYaSFEUaxVFcUdRFB8nGZfkmyRbJ1lxGqtP+YH4tYl/3jupoVKpDEvySZKlprH9jVM8vj5Jm9SMWEuS9ZM8NelD+8T9fZ3kHyXrBoBGbX7qr4uiOCbJT5McXKlU3preegDAt1cUxW5FUfy7KIovUtMnf52a/t1naJjPCNqhARRFsVRqRpYtkuSQJBuk5hKze5K0mMYmw6Z4PHYG7dPa/uPpPJ40/+oS01hnqu1mo24AaLTmp/66KIoDkvw+ye8qlcoV01oHAJg1RVHsmOSGJK+m5svs9VLTJ/8vPkPDfMfUMdAwtkuycJLdKpXK+5Mai6Jo1UDPt3iSd6Z4nCRDJ/754WRtU243uTldNwBU03zRXxdFsXuSC5OcValUTpmNugGAaeud5K1KpbLXpIaiKJqn5svwyY2euGyBSqUydrL2RRu8QmCOMaIdGsakD7rfTGooimKFJBs20PPtNsXj3km+SvLSxMdPJuk5cQTcpHpap+bGaZOb03UDQDXN8/11URS7pObGp5dNPq87ADDLxiRpOUVbq9RMFzO53ZM0naJt0v1WJk0Xl6Io2qXmqjRgHmFEOzSMB1LT2V5TFMVZqbkU/MQkQ9IwX3D9cuKN0J5Jsm1qbsbSb+KN1ZKaG7D8KjU3XuuXmhOEI5OMmt26i6JYO0mXydq7FUXxo4m/31WpVEbW26sDgIYxT/fXRVFskuS6JC8muaooip6T7WNMpVJ5vn5fHgDM015JskhRFAcm+U9qRqnfk2TnoijOSXJHkrWS/DrJF1Nse3eS4UkuLYqib5IFkxyVmi/cgXmEEe3QACqVystJfpZkmSS3p6YDPTrJvxroKXul5gZotyf5eZL+SU6erJ5Pk2yZ5NMkVyf5U2pOCOrM0TqLdR+c5KbUzEeXJLtOfHxTksXq52UBQMOZD/rrLVLzQX6NJI+nZsT8pJ9b6/WVAcC877LU3Mj890meTs3Nyi9NckqSH098vH1qrkQbPvmGE79U3yHJhNTcHP3UJOcneWgO1Q7MAUWlUql2DQAAAAAA0GgZ0Q4AAAAAACUI2gEAAAAAoARBOwAAAAAAlCBoBwAAAACAEgTtAAAAAABQQrOZLK/MkSoAYP5WlNxefw0ADa9sf12Gvh4AGl6pvt6IdgAAAAAAKEHQDgAAAAAAJQjaAQAAAACgBEE7AAAAAACUIGgHAAAAAIASBO0AAAAAAFCCoB0AAAAAAEoQtAMAAAAAQAmCdgAAAAAAKEHQDgAAAAAAJQjaAQAAAACgBEE7AAAAAACU0KzaBcDcYMDAofW6v149OtXr/gBgXlGffW599rdza10AADAj9Z1pTYvz22/HiHYAAAAAAChB0A4AAAAAACUI2gEAAAAAoARBOwAAAAAAlCBoBwAAAACAEppVuwAAAJgdAwYOrXYJAAAASYxoBwAAAACAUgTtAAAAAABQgqAdAAAAAABKELQDAAAAAEAJgnYAAAAAAChB0A4AAAAAACUI2gEAAAAAoARBOwAAAAAAlCBoBwAAAACAEppVuwAAAOrfgIFD621fvXp0qrd9AQAAzIuMaAcAAAAAgBIE7QAAAAAAUIKgHQAAAAAAShC0AwAAAABACYJ2AAAAAAAoQdAOAAAAAAAlCNoBAAAAAKCEZtUuAACAuduAgUOrXQIAAMBczYh2AAAAAAAoQdAOAAAAAAAlmDoGGkBDXGLfq0enet8nAAAAAFCeEe0AAAAAAFCCoB0AAAAAAEoQtAMAAAAAQAmCdgAAAAAAKEHQDgAAAAAAJQjaYR538sknZ6uttkrbtm1TFEUGDx481TqVSiUXXHBBunfvnlatWqVLly455JBD8sUXX8z5ggGgkXv99ddz0EEHZeWVV06rVq2y3HLL5Te/+c1U/erNN9+cDTbYIIsuumhatGiRFVdcMf3798/YsWOrVDkAwPxnzJgx6dOnTxZbbLG0bt0622+//TSzkyltttlmKYpiqp/Ro0fXWe+xxx7L+uuvnxYtWmTJJZfMcccdl3HjxjXQq6GamlW7AKBhXXzxxenatWs233zz3H777dNc5/zzz8+hhx6a448/PptttlneeOONHHvssRkyZEgGDBgwhysGgMbt/vvvz+OPP54DDzwwq622Wt5555387ne/y5NPPpmnnnoqTZrUjHX57LPPsvnmm+fII49Mu3bt8vTTT6dfv3756KOPcsEFF1T5VQAAzB9+/etf5+abb84555yTjh07pl+/ftl6663z3//+Ny1atJjhtptvvnl+//vf12lbcMEFa38fNGhQtt5662y77ba59dZb89Zbb+WYY47J119/nXPPPbdBXg/VU1QqlRktn+FCmFcMGDi02iXMVK8enWZruwkTJqRJkya54447suOOO2bQoEHp0qVLnXV69uyZTp065ZZbbqltO++883LYYYflyy+/TOvWrcuUDsxcUXJ7/TVTaQx927xqo6VaZJFFFklR/P9/7fvuuy/bbrttHn744Wy66abT3fa4447Ln/70pwwbNqzO9sBcoZr/KfX1AA3g/fffT5cuXXLFFVdkjz32SJIMHTo0yy67bC688MLsu+++0912s802S4cOHXLzzTdPd539998/999/f9544400a1Yz3vm8885Lnz59MmTIkCyxxBKlX8OcOO+f3UyqESrV15s6BpL88YTD0uen389/Hv1nDv7B5tlt/a45+ZA9MmL4sHw4ZFB+98td8+P1v5s+P/1+Br/xSu12EyZMyC1XXJADdtowP1p3ufyq18Z58Pab6uz7P4/+M30P+En23GL1/GSjlXLUHjvm+ScfqbPOdRedld03XzXvvPZSjtpjx+y2ftcc1nvbvPzcv0u/tkmj5mbkm2++ycILL1ynrV27dqlUKpnJl3EAMNvm1f530UUXnSokX2ONNZIkn3zyyUy3NXUMADC32WuvvbL22mvnzjvvTLdu3dKqVatsv/32+fzzz/PWW29l8803T+vWrbP22mvnxRdfrN1uwoQJOe2009K1a9csuOCCWWGFFXL11VfX2fedd96ZrbfeOosttljatm2bnj175r777quzTr9+/dKhQ4c8//zz6dmzZ1q1apU11lgjjz76aKnXNel5fvCDH9S2derUKRtttFHuvvvuUvtOkoEDB2azzTarDdmTZJtttsm4ceOmeo00foJ2mOjTj4bmb38+Mz876Mj86nd/yGsvPJsLT/5tzjz6oGy8ba8cdcbFGT9+XM48+qDa8PnSPxyfmy47L9v84Gc5/ryrs94W2+WCE/vkmX89ULvfj4cOyTqbbJVD+/8xvz3zkqy02to5+eDd8+rAZ+o8/5jRo/LHEw7Ntj/8eX57xiVpvsACOa3PvhkzalTtOhMmTMi4ceNm+DN+/PhZfu377rtvbrzxxtx1110ZMWJEnn/++Zx22mnZa6+90qZNm9k8ogAwc42l/x0/btyMf2bS/z7xxBNJkm7duk21bPz48Rk5cmQee+yxnHfeeTnwwAONZgcA5jpDhgzJCSeckP79++eSSy7JE088kf322y+9e/dO7969c/PNN2fcuHHp3bt37XnbIYcckv79+2e//fbLnXfemV122SV777137rjjjtr9Dho0KDvuuGOuvfba3HLLLdlggw3yve99L48//nid5x85cmT23HPP7L///rnllluy4IILZpdddsnIkSNr15nV3OS1115L586dp8o+Vl555bz22mszPSb33XdfWrVqlVatWmXbbbet8yVDkowePToLLLBAnbZJU8u8+uqrM90/jYs52mGiEV9+kdOuHpAlluqSJBn85qu57eqL8puTzs3mO/4oSc1NQ/v/es+8P+itNGvWLPfcdE0O6Xd2tthp1yTJ6j03zrD/fZIbLj4n62yyVZJk+96/qH2OCRMmZNW1N8iQd97IA7ddn5V7rFO7bOzo0dnniBOz2robJknad1w8h/feNi8/91TW3HDzJMnee+891Te/U9p0003z8MMPz9JrP/DAAzNixIjsuOOOmTBhQpJk5513zsUXXzxL+wGAWdUY+t/z+/XJQ/+oO2J+St3X6plTLpv2ZcMjR47M0UcfnU033TTdu3efannr1q0zZsyYJMkee+yRM844Y6bHDQBgTvv888/z5JNPZvnll0+SvPjiiznjjDNy9dVX1067UqlUsv322+e1115L8+bN8+c//zlXXnll9txzzyTJVlttlQ8//DAnnnhidthhhyTJwQcfXPscEyZMyOabb56XX345l19+eTbccMPaZaNGjcq5556bLbbYIkmyxBJLZI011si//vWvbLfddklmPTcZNmxY2rVrN9U67du3z7Bhw2a6nz333DNdu3bNu+++m1NOOSUbb7xxXnjhhdope1t16JT7//VEneldHrv3niTJc2++Z7rHeYygHSZabMnOtR/yk9T+vuq6G/x/29LLJkk+/+SjfPT+uymaNEnPLbbL+MnuFr3auhvl0XsHZPz48WnatGk+/fiD/PWC0/PCvx/NsE8/qf1Wd/IP+UnSrFnzrLL2+rWPl1ruu0mSzz7+sLatX79+dTqgaVlooYVm4VXXuO6663LSSSelf//+2XDDDfP222/n+OOPzz777JNrrrlmlvcHAN9WY+h/ex9weLb/8V4zfB0tW0/7CrBKpZJ99tknn3zySe68885prvPEE09k5MiRefrpp3PSSSfl4IMPzoUXXjjD5wMAmNO6dOlSG7InSdeuXZOkNvievG3o0KF5++2306RJk+yyyy4ZN9l525Zbbpnrrruu9rzt/fffz3HHHZcHHnggH374Ye152+Qhe5I0b948m222We3jSVcKvv/++7Vts5ObTOtKwkqlMtMrDE888cTa3zfeeONstdVWWWmllXLuuefW3uj0e7vukb4H/iQ3XHJuvrfrHvnwvUG55rxT06Rp0zRt0nSG+6fxEbTDRK3b1J2jvFnzmkt7Wi+08GRtzZMkY8eOyZdffJ4J48fnpxuvPM39Dfv04yzS8Tv5/aF7Z9TXX+UnBx6RJZbqkhYtW+Vvfz4zwz//rM76Ldu0qTOfevOJzz927JjatqWXXjqdO3ee4euY1UvNJ0yYkEMOOSS/+c1vcswxxyRJNtlkkyy55JLZbrvtcuihh2bNNdecpX0CwLfVGPrfjt/plA6LzeRGVdPpf3/729/m1ltvzf3335/llltumutM6mc32mijdOjQIXvuuWf69OlT54MsAEC1TTnye9KUKJO3T2obPXp0Pv3004wfP36qe8JN8uGHH2bJJZfMTjvtlBEjRuSkk05K165d07p165xwwglT3dumbdu2dc7bJn+uSWY1N2nfvn2++OKLqdb54osvpjnSfUa+853vZMMNN8xzzz1X27Z6z43zs4OOyo2X/jHX/fnMNGvWPLvtd2juvO6KLLxoh1naP3M/QTvMpjZt26Vps2Y59cpb06SY+nYHCy/SIR++NzjvvPZSTrjg2trLz5Nk7JjRU63/bTTE1DGffvppPvvss/To0aNO+6Sbtr399tuCdgDmGtXof2d36phzzjknZ555Zq6//vpsvPHG3+q5JvW5gwYNErQDAI3aIosskmbNmuXxxx+vE5BPsthii+Wtt97K888/n7vvvrt2+pekZpqY2TGruclKK62U9957L19//XVat25du85rr72WlVZaabZqmHIA5K77/jo7/HSffDx0SDosvkQmjJ+Qv114RlZcVdYyrxG0w2xabd0NM2H8+Iz8akR69NxkmuuMnfitavPJbnzxyQfv57WB/8ky3532SLwZaYipYzp27JhWrVrlueeey6677lrb/uyzzyZJ7bxiADA3qEb/OztTx/ztb39Lnz59ctZZZ2W33Xb71s816aZfyy677CzXCQAwN9liiy0yfvz4DB8+PFtvvfU015kUqE+6QWiSvPvuu3n88cez2mqrzfJzzmpuss022yRJbr311vz85z9PknzwwQd59NFHZ3kqv48//jiPP/549t5776mWtWzVOl0mnodef/HZ6bhE56y+3rcbiEHjIWiH2dSpy/LZ9ke756yjf5Vd9jwwy3dbPd+MHZMhb7+eD959Jwf3PTOdl10+iy6+RK48++T89FdHZtTXX+W6i87KIot9Z7aes0uXLrMcfD/yyCP53//+Vxuc33333enYsWO6deuWbt26pSiK7LfffjnnnHPSqlWr2jna+/btm/XWWy9rrbXWbNUKAA2hGv3v4ksulcWXXOpbr//II49kr732yjbbbJP1118/Tz31VO2yzp07117OvN1222WrrbZK9+7d07Rp0zz++OM566yz8uMf/9hodgCg0VtxxRVzwAEHpHfv3jnqqKOy9tprZ/To0Xn55Zfzxhtv5LLLLstKK62Uzp07p0+fPjn55JMzYsSI9O3bN506dZqt55zV3KRz587ZZ599cuihh6ZSqaRjx47p169flllmmdrgPUlOOumknHTSSbVzzb/44os55phjsuuuu2aZZZbJkCFDcuqpp6ZJkyY59NBDa7f7cMig/Ovu2/LdVXpk/Pjx+c+/Hsg/B9yQ351/dZo2E8vOa/yNQgn7H3NKllxmudz/97/lb38+K61at8lSy303W+38kyRJ8wUWzNFnXZpLTj0upx+5XxZdbIn8aN9f56X/PJkhb70+R2rs27dvHnnkkdrHv/rVr2rb+/XrlyQ57bTT0qFDh1x77bU59dRT07Fjx+ywww7p37//NC/vAoBqmtv734ceeijffPNN7r333tx77711lk3e/66zzjq56qqrMnjw4DRr1izLLbdcTj311BxwwAENXiMAwJzwpz/9KSussEIuvfTSnHDCCWnbtm26deuWffbZJ0nNSPa///3vOeigg/KjH/0onTt3znHHHZeHH344L7300hyp8bzzzkvr1q1z+OGHZ+TIkdl0001z3XXXpUWLFrXrTJgwIePHj699vOiii6ZSqeSYY47JZ599loUWWiibbbZZbrvttiy99NK16zVrvkBe+PdjGfCXSzJ+3Lh07b56Trz4+nRfc7058tqYs4pJd/KdjhkuhHnFgIFDq13CTPXqMXvf5gKNwqzdxXhq+mum0hj6tnmVPhvmWWX76zL09QCN0LxyTj4fnd+W6usNVQUAAAAAgBJMHQMAAAAAwDTNiZH588KoeUE7Daoh/iPOC//xAAAAAIB5h6ljAAAAAACgBEE7AAAAAACUIGgHAAAAAIASBO0AAAAAAFCCoB0AAAAAAEpoVu0CAABgXjJg4NB62U+vHp3qZT8AAEDDM6IdAAAAAABKELQDAAAAAEAJgnYAAAAAAChB0A4AAAAAACUI2gEAAAAAoARBOwAAAAAAlNCs2gUA1TFg4NB632evHp3qfZ8AAAAAMLczoh0AAAAAAEoQtAMAAAAAQAmCdgAAAAAAKEHQDgAAAAAAJQjaAQAAAACghGbVLgAAgP83YODQapcAAADALDKiHQAAAAAAShC0AwAAAABACYJ2AAAAAAAoQdAOAAAAAAAlCNoBAAAAAKCEZtUuAAAAAABgThswcGi1S2AeYkQ7AAAAAACUIGgHAAAAAIASBO0AAAAAAFCCOdqhkTBvGAAAAADMnYxoBwAAAACAEgTtAAAAAABQgqAdAAAAAABKELQDAAAAAEAJgnYAAAAAAChB0A4AAAAAACUI2gEAAAAAoARBOwAAAAAAlNCs2gXArBowcGi1SwAAAAAAqCVoBwAAAADmKgZa0tiYOgYAAAAAAEoQtAMAAAAAQAmCdgAAAAAAKEHQDgAAAAAAJQjaAQAAAACgBEE7AAAAAACUIGgHAAAAAIASBO0AAAAAAFCCoB0AAAAAAEoQtAMAAAAAQAmCdgAAAAAAKEHQDgAAAAAAJTSrdgEAAAAAML8YMHBog+6/V49ODbp/YNqMaAcAAAAAgBIE7QAAAAAAUIKgHQAAAAAAShC0AwAAAABACYJ2AAAAAAAoQdAOAAAAAAAlCNoBAAAAAKAEQTsAAAAAAJQgaAcAAAAAgBIE7QAAAAAAUIKgHQAAAAAAShC0AwAAAABACc2qXQAAADC1AQOH1tu+evXoVG/7AgAApmZEOwAAAAAAlCBoBwAAAACAEgTtAAAAAABQgqAdAAAAAABKELQDAAAAAEAJgnYAAAAAACihWbULAAAAAADqx4CBQxv8OXr16NTgzwGNjRHtAAAAAABQghHtAAAAAMC3NidGzUNjI2gHAIB5XH1+GHapOAAATM3UMQAAAAAAUIKgHQAAAAAAShC0AwAAAABACeZopw43swAAAAAAmDVGtAMAAAAAQAmCdgAAAAAAKEHQDgAAAAAAJQjaAQAAAACgBEE7AAAAAACUIGgHAAAAAIASBO0AAAAAAFCCoB0AAAAAAEoQtAMAAAAAQAmCdgAAAAAAKEHQDgAAAAAAJQjaAQAAAACgBEE7AAAAAACUIGgHAAAAAIASBO0AAAAAAFCCoB0AAAAAAEoQtAMAAAAAQAmCdgAAAAAAKEHQDgAAAAAAJQjaAQAAAACgBEE7AAAAAACUIGgHAAAAAIASBO0AAAAAAFCCoB0AAAAAAEoQtAMAAAAAQAmCdgAAAAAAKKFZtQsA5h0DBg6tdgkz1atHp2qXAAAAAMA8xoh2AAAAAAAoQdAOAAAAAAAlCNoBAAAAAKAEQTsAAAAAAJQgaAcAAAAAgBIE7QAAAAAAUIKgHQAAAAAAShC0AwAAAABACYJ2AAAAAAAoQdAOAAAAAAAlCNoBAAAAAKCEZtUuAACgsRswcGi1SwAAAKCKBO0AAAAAEAMogNln6hgAAAAAAChB0A4AAAAAACUI2gEAAAAAoARBOwAAAAAAlCBoBwAAAACAEgTtAAAAAABQQrNqFwAAADQeAwYOrXYJ09SrR6dqlwAAwHzMiHYAAAAAACjBiPZGbG4dTdQYPfXQPbnuz2dm6OB3skjHxbN971+k1+77TXf9y87omzv+dnl67b5/fnH48bXt7w96K3884dAMHfx21txg8xzU94y0bNW6dvnLzz6Vs445KH+67V912qfln7ffmPP7Hp7rHn99qnWvu+is3HXDVbn2of8mST7+4L3sv/36tctbtGqdTsssn132OjAbbbPjNNdZsEXLtOuwWFZcdY1s+6Pd033N9b7FkQKA6tJn67MBgPKcUzmnov4Z0c5879WBz+QPfX6Z73bvkeP+eGW27PXjXHPe73P7Xy+b5vrvvf1G/jnghrRqs9BUy87re1iWWKpLjvzDn/PeoDdy8+Xn1y6bMGFCLj+zX3Y/5JiZdi6za6/Djs8frh6Q3555SZZcetmc+dsD88y/HpjmOseff012++VvMmL4sBy3zw9z/cVnN0hNAFBf9Nn6bACgPOdUzqloGEa0M9+74ZJzs3KPdXNw3zOTJGusv2m+/nJ4brzknHxvtz3SvPkCdda/9PQTssNP9snDd95Sp33UyK/zxn+fz3HnXpWFF1k0X4/4Mrddc3F2P+ToJMkDt12Xps2aZbMdfthgr6VTl+Wz4mprJUlWX2/jvP3af3PPTddknU22muY6q2T9bLnTbvnbhWfk+ovOTve1embVtTdosPoAoAx9tj4bACjPOZVzKhqGEe3M9wa9/nJWX2+jOm091t8kX305PK+/8Gyd9ifuvyPvD3ozP9z7oKn2M+6bsUmSBVq0mPhny9q2kV+NyN8uPDP7HnliiqJoiJcxlSZNmmTZFbvnkw/en+m6P97/8CzScfHce9Nf5kBlADB79Nn6bACgPOdUzqloGIJ25nvfjB2TZs2b12lrvkDNt7fvD3qrtm3M6FG58uyTs/uvj0mLlq2m2s9CC7fP4p2Wzp3XXZkRw4flvo53KWEAACAASURBVL//Nct3Wy1JcuOl52b19Taq/QZ1VkwYPz7jx42r81OZMOFbbfvJB++lXYeOM12vadOmWXXdDfP6f5+b5foAYE7RZ+uzAYDynFM5p6JhmDqG+d53luqSN19+oU7bmy8NTJKM+PKL2rZbrrgg7Tsuls22n/4lT/sd3T9nHHVA/nLBaVly6WWz/9Gn5MMhg/LAbdfn3Bvvn636frZJt2m2L9Su/VRtlQkTMn7cuIz8ekTuv/W6vPnSwPzytyd/q+fpsNgSGf75/2arRgCYE/TZNfTZAEAZzqlqOKeivgname9t96Of56LfH5v7/v7XbLDV9nnzpYEZcO0lSZKmTWou+vh46JDcdu3FOfniG2Z4ydNaG22Rq//5Qj795MN8p/Myadq0aU75zS+y489+mQ6LL5m7brgqf7/qwiTJD35xUL6/254zre+Uy2/Jggu2qNN279//mqcevHuqdX9/2N61vzdr1jw7/Xy/bLfrHjM/CEkqlcq3Wg8AqkWfXUOfDQCU4ZyqhnMq6pugfQ4aMHBotUtgGrbs1TuD33glF/3+2Fx48m+zYIuW2eM3x+bSPxyfhRetudzomvNOzZobbJ5Oy3bNVyOGJ0kqlQn5ZuyYfDVieFq3aVvb8SzYsmU6LbNckuSFpx7N4DdeyRF/uDCDXn8lf7vwjJx21YAkyW/37JVuPdZJlxWm/U3tJMuttMpUd+d+5tEHprnu3kf0S7ce66Rl6zZZrNNSU93AZEY++99HWXiRmV9eBQDVos+uoc8GAMpwTlXDORX1TdDOfK9p06bZ7+hT8tNfHZlPP/4wi3daOkMnzkm24qprJkmGDn47g994ZapvT++64arcdcNVueyep9Nh8SXrLBs/blwuP7Nf9jz0uCzYomVeevaJrLrOhum8bNckyWrrbpiXnn1qph3MrFhiqS7p2n31Wd5u/Lhx+e/Tj6f7mj3rrRYAqG/6bH02AFCecyrnVDQMQTtM1KZtu7Rp2y5JcvdN12Sl1deu7QwOPuGMjBr1dZ31zzr6oHRfq2e223X3LNx+0an2d8/N16Z124Wz0bY71baNGT3q/38fNTKZSy5TuuGSc/L5/z7Otj/6ebVLAYCZ0mfrswGA8pxTOaeifgname+9/uKzeXXgM1l2he4Z+fWIPHrPgDz/5CM59Yq/164zrW9Hmy+wYDosvmRWXXuDqZZ99eUXueHic9L3wr/WtnVfs2euPvf3eeC265MkLz7zRHb/9bEN8IpmbOjgt9O2XfuM++abfDz0vTx274A898TD6X3A4Vll7fXneD0A8G3ps/XZAEB5zqmcU9EwBO3M95o2a57H7v1Hrr/o7BRNmqTbGuvm1CtvTZfvrjzb+7zuz2dlnU23zvIrr1rbttxKq2TPQ4/NX/90epJkr8N+l2VXrL/Lpb6tq86pufv2AgsumPYdFs+Kq62ZUy6/Jd3XXG+O1wIAs0Kfrc8GAMpzTuWcioZRzOQOu3PH9RzzCDdDhXlTrx6dql0CjV9Rcnv9dZXp46H69MfMAWX76zL09TCHOK+D6phLzuVK9fVN6qsKAAAAAACYH5k6BgCYLxmtBAAAQH0xoh0AAAAAAEoQtAMAAAAAQAmCdgAAAAAAKEHQDgAAAAAAJQjaAQAAAACgBEE7AAAAAACUIGgHAAAAAIASBO0AAAAAAFCCoB0AAAAAAEoQtAMAAAAAQAmCdgAAAAAAKEHQDgAAAAAAJTSrdgEAjd2AgUPrdX+9enSq1/0BAAAA0LCMaAcAAAAAgBIE7QAAAAAAUIKgHQAAAAAAShC0AwAAAABACYJ2AAAAAAAoQdAOAAAAAAAlCNoBAAAAAKAEQTsAAAAAAJTQrNoFAAAAAMDMDBg4tNolAEyXEe0AAAAAAFCCEe0AAECjV1+jHHv16FQv+wGY2zT0aHDvn8D8zoh2AAAAAAAoQdAOAAAAAAAlCNoBAAAAAKAEQTsAAAAAAJQgaAcAAAAAgBIE7QAAAAAAUIKgHQAAAAAAShC0AwAAAABACYJ2AAAAAAAoQdAOAAAAAAAlCNoBAAAAAKAEQTsAAAAAAJQgaAcAAAAAgBKaVbsAAOoaMHBove+zV49O9b5PAAAAAGoY0Q4AAAAAACUY0Q4AAABAKQ1xZS5AY2JEOwAAAAAAlCBoBwAAAACAEkwdMx0ueQIAAAAA4NsQtAMAAABUkcF+AI2fqWMAAAAAAKAEQTsAAAAAAJQgaAcAAAAAgBLM0Q4ANKj6nHO0V49O9bYvAAAAqC9GtAMAAAAAQAmCdgAAAAAAKEHQDgAAAAAAJcwTc7TX59yv0/PhkEG59ZqL8saLz2XI269n5TXWzSmX3Vy7/L//eSLH/3K3aW7bY/1N0+/CvzZ4jcC8b2bvRUny1YjhufKsk/Lvh+7JuG++Sbc11k33qy9N165dq1Q1zJ0ev/+OPHzHLXn71Rcz8qsRWbLL8tl59/2zyfd2rl3nsXtvz2P33Z7XX3wuwz79JIeceHa23Gna/T0w77rpppty7bXX5tlnn83w4cOz4oor5ogjjshPfvKTaa5/7rnn5rDDDssPf/jD3HzzzdNcB2Bu8d7bb+SS04/P6y8+m9ZtFs7Wu/TOj/c/PE2bNp3uNh9/8F723379qdo32manHPGHC2sf77xG52lu36z5Arn56XfKFw80apO//yzavn323Xff9O3bd4bvP5P8/e9/z6mnnpqXXnoprVq1yjrrrJNbbrklrVu3TpLcf//9ueKKK/Lkk0/m3XffTd++fdOvX78GfT3zRNA+Jwx5+408+9iDWXHVNTNu3DdTLV9+pVXzh6sH1Gn730cf5MzfHpg1N9x8TpUJzONm9l6UJGf+9lcZ8tZr2efIE9O6TdvcdNkfs+WWW+a///1v2rZtO4crhrnX7ddeksU6LZW9j+iXtu0WybOPPZizjz04X37xeXb4yd5JkiceuDOffPB+1t54y9x/63VVrhiolrPPPjvLLrtszjnnnHTo0CF33XVXfvrTn+bTTz/NIYccUmfdTz75JCeddFI6duxYpWoBvr2vvvwiJxzwkyy13Hdz7DlX5KP33s2VZ5+USqWSnx101Ey33+uw47Nyj7VrH7dtt0id5VPmJElyym9+kZV6rFO+eKBRm/L9Z/EMT58+fTJhwoT0799/httedtllOfjgg3PUUUfljDPOyLBhw/Lggw9m3Lhxtevcc889efHFF7Plllvm+uuvb+iXk0TQ/q2ts+nWWW/zbZMkfzhiv3z5xed1lrdqs1BWXG2tOm2vPP90mjRpko223mGO1QnM22b2XvTaC89m4JOP5KSLb8hq626YJFlh1TVy4I4b5JJLLskRRxwxx2uG+lSfV7Ed98er0rb9/38YXG3dDfP5/z7K7X+5tDZoP+IPf06TJk0yauTXgnaYj/3jH/9Ihw4dah9vscUW+eCDD3L22WdPFbQfc8wx2X777fPee+/Vaw31+f7Xq0enetsX0Ljdc9NfMnbM6Bx91qVp1WahpGcy8usRuf7is7PLngfWtM1Apy7LT5WFTG7KZW++PDBffvF5NtmuV73UDzReU77/9OrRKV9++WX69euXo446aroDBT/99NMcdthhOf/88/PLX/6ytn2XXXaps94ZZ5yRs846K0kyYMDUX/o1BHO0f0tNmsz6oXr0ngHpvlbPLLLYdxqgImB+NLP3okGvv5ymzZql+1o9a9vaLdoxq622Wu68886GLg8alclD9kmWW2mVDP/809rHs9P/A/OeyUP2SdZYY4188sknddqeeeaZ3HjjjTnttNPmVGkApTz3+ENZY/1N6wTqG2/bK2NHj85Lzz5V78/36N0D0qJlq6yzydb1vm+gcZnW+0/v3r0zatSoPPLII9Pd7sYbb0yS7LnnnjPcfzU+y/n02EA+ePedvPPaS9nYt7TAHPTN2DFp2rTpVPOZLbjggnn11VerVBU0Hq+98J8stdwK1S4DaASeeOKJdOvWrfZxpVKpvYS5UycjxoHGYejgt9Jp2eXrtHVcolMWbNEyQwe/NdPtz+97eH6w1tL5xdZr5oozT8yY0aOmu26lUsnjD9yRdTfbJgu2bFm6dqBxm9b7z9JLL51WrVrltddem+52//73v7Piiivm8ssvT+fOndO8efOst956eeKJJxq65JkydUwDefSeAWnWrHnW3/L71S4FmI98Z6kuGTtmTAa/+Wq6fHflJMmY0aPy0ksvZcSIEVWuDuZuL/z7sTz98H05uN9Z1S4FmMv985//zIABA3LFFVfUtl155ZX56KOPTNMGNCpfjRie1gstPFV7m7YL56svh093u+bNF8j3f7xnevTcNC3btMlL/3kyt151YT56/90ce+4V09zmlef+nc8+/jAbbWtAIjD995/27dtn2LBh093uo48+yuuvv57+/fvn9NNPz6KLLprTTz892223Xd58880svvjiDVn2DAnaG8ij996eHutvkoUWbl/tUoD5yBobbJrFOy2dP/c/Or8+8ey0bN0m1553aoYPH57mzZtXuzyYa338wXs5+9iDs+5m22TLnXardjnAXGzw4MH56U9/ml69emWvvfZKkgwfPjzHHntszjvvvLQ0ShNoZIpi6rZKpZJiWgsmWqTj4tnv6FNqH6+69gZpt0jHXHzqsXnn9Zez3Irdp9rm0XsGpE3bhbPGBpvWS91A4zc77z8TJkzIV199lZtuuinbbbddkmSDDTbIMssskwsuuCAnn3xyQ5U7U4L2BjDo9Vfy/qA3s+u+h8x8ZYB61Lz5Aulz2p9y1jEH56Bdak5gV15j3eyxxx558MEHq1wdzJ1GDB+Wkw/ePR2/0ymHnXJ+tcsBqmxGNx0dMXxYjvnFLmnbYYn0Pur02nWv/uMpWWjR7+SbJVbJXx99JUnyyfCRWXjhb/LFF19koYUWmmpaN4C5QZuFFs7XI76cqn3kVyPSeqFp34hwejbYavuaoP3V/04VtI8fNy5P/vOurL/l99O8+QKlagbmDdN7/xk+fHjatWs33e0WWaTmXlubbbZZbVvbtm2z1lpr5ZVXXqn3OmeFoL0BPHbvgCzQokXW3WzbapcCzIdWWGWNXHT7Y/ng3XfSpGnTLLFUl1z6u/3Ts2fPmW8M85kxo0al/6/3yjfffJP+l12dFi1bVbskYC41o/eLoYPfyVuvvJCfb1I3WHp1YM3lz48++mg22mijOV0ywEx16tI1Qwe9Xaftfx99kNGjRqZTl66ztK9JI1CnNRL1xacfy/Bhn2Xj7Xae/WKBecq03n/ee++9fP3111lppZWmu93KK6+coihSqVTqtFcqlarcAHVybobaAB677x9ZZ5Ot07JV62qXAsyniqJIpy7LZ4mluuSDd9/JAw88kH322afaZcFcZfy4cTn9qP3z4ZBBOeGCa9NukQ7VLgmYS83s/eJnBx2Zky+9sc5PlxW6ZZNNNslDDz2UVVddtUqVA8zYmhtunueffDijvv6qtu3x+27PAi1aZJW1Zm2gzhMP3JkkWX7lqd/z/nXPgLTvsFhWWXv9cgUD84xpvf/ccMMNadmyZTbddPpTTO2www6pVCp56KGHatuGDx+eZ599NquvvnqD1jwzRrR/S2NGjcqzj/0zSfL5Jx9l5Ndf5Yn770iSrLXRlrV3zH79xWfz8dAh+UWfE6pWKzDv+jbvRTdccm46L9s1bdu1z7tvvpYbL/1jevfuna233rqapcNc56JTj82zjz2YfY88MV8NH5bXX3y2dtlyK62S5gssmPfefiPvvfNGxo4dkyR5+5UX0rJlq7Rtv6gPijAfmdn7xTJdpx511XqhtvmmeZsMb/fdPDzoqyRfTbUOQLVtt+vPc+f1V+S0Pr/MD/b6VT4a+m6uv+js9Pr5fmnVZqHa9Q7YacN0X7NnDpl40/jrLjoro77+Oiv3WDutWi+Ul597Krddc1F6bvG9dFmhW53n+GbsmPz7oXuzxU67Vn20KTD3mPL95+Onv0y/fv1y+OGHp23b/5+6qmvXrtl0001z+eWXJ0nWXnvt9OrVK/vss09OO+20dOjQIaeffnqaN2+egw46qHa7d999N88880ySZOzYsXnllVdy8803p3Xr1vne977XIK+pKkH7jOY9nFt9MezTnH7UAXXaJj2++M4ns3jLpZLU3AS1VZu2WWvDzed4jcC879u8F40YPiyXn9E3X34xLB2+s0R23mP/XHbWSdUoF+ZqA5/8V5LksjP6TrXs4jufzOJLLpXH7v9Hbrj4nNr2u264OnfdcHW6r9Uzp1x28xyrFaiub/N+AdAYtWnbLidddH0u+cPvcsqhe6X1Qgtnx5/9Mr0POLzOeuPHjc+ECRNqH3fu0jW3XXNxHrjtuowdPTodllgyO+9xQHbd99dTPcezjz+UkV99mY233anBXw/QeEz5/rNI+/Y57LDD0q9fvzrrjRs3LuPHj6/T9pe//CVHHnlkDj/88IwcOTIbbrhhHnzwwbRv3752nYceeii/+MUvah/fdNNNuemmm7LMMstk8ODBDfKaiinns5nCDBfOrsYYtAM0Zr16dKp2CczY9G+p/u00SH9dX/T7AHMX5wWzrWx/XcZc3dfPyLxwHjAn/s/MC8cJoIy55PykVF9v6hiA+UBDnLjPJZ0gAAAAQNWZHAsAAAAAAEowoh0AmCaXMAMAAMC3I2gHYLbUdwhrKhoAAACgsTJ1DAAAAAAAlCBoBwAAAACAEgTtAAAAAABQgjnaAZgrNMSNN837DgAAAMwJgnYAAID5SH1+ue1LbeYHDTEgBIB5j6ljAAAAAACgBCPaAQAAYD5kpDYA1B8j2gEAAAAAoAQj2uvBbddcnOsvPjujR35d7VKA+VyLVq3Te//Ds/Me+1e7FGjU9O3ArNIHA/Mb50tAtbRp0yb9+vVLnz59ql1KHUWlUpnR8hkunF3z2uVpv9h6rQz79ONqlwGQJGnfYfFcef+z1S5jrtCIbtBWlNxef13P9O3A7NAHz75G0meX7a/L0Ncz13G+BFTTEksskQ8++KC+d1uqrzeivR702n0/3+ICc4UWrVqn1+77VbuMuUZ9f3hsJCEA9UDfDswqfTAwv3G+BFRLmzZt5rrR7IkR7QDwrTVg0G5EOwDztUbyZbYR7QDQQOaScwEj2gEAAGi86jPwnUs+qAMA85km1S4AAAAAAAAasxlOHXPiiSfek6TDnCunapZMUu+z5zcyjoFjMInj4BgkjkEyZ4/Bp3379t1udjeej/rrGfFvtv45pg3Dca1/jmnDcFynVqq/LuPEE098K8kX1XjueYB/y7PHcZs9jtvsc+xmj+M2e6Z33Mr19ZVKZb7/6devX6XaNVT7xzFwDBwHx8AxcAwa84+/L8e0sfw4ro5pY/lxXOeuH38fjp3j1jh+HDfHznFrHD8NddxMHQMAAAAAACUI2mucWO0C5gKOgWMwiePgGCSOQeIYNDb+vuqfY9owHNf655g2DMd17uLvY/Y5drPHcZs9jtvsc+xmj+M2exrkuM1wjnbg/9q7/6A7qvqO4+8PNSAQ0CIjQgsCA0EENCIVQX48YFMYkJqCBaVaYwsMtQKiEYZWBRSltApNAesULSkDAlL5UWj5FWxQIToSxwINhCAkQggC1kKBAAFP//ieB66X+/u5e/fu5vOa2QnP3r3L93zvuXvPnj171szMzMzMzMzMzKwzj2g3MzMzMzMzMzMzM5sCd7SbmZmZmZmZmZmZmU2BO9rNzMzMzMzMzMzMzKagVh3tkvaR9G+SVkpKkua02GaGpCsl/a+kZyX9RNKOHfY5kffVvLyl0MIMqFsO2pQlSTq/y353kXSrpNV535+XpEILM6AiciBp6zbvObDwAg2ohzxMl3SupIfz57pU0ok97HdfSYslPSfpAUnHFlaIKSoiBzU8Jmwmab6kR/Ix8QZJ2/ew3zrVg75zULV6UAeSTpH0Y0lPSXpc0rWSdm7a5lBJN+bXk6SJksKthG45lTRN0lmS7pT0jKRVkr4laasy4x53PdbVL0q6N+f1V5JukbRnWTFXQS95bdr+n/JxYO4o46ySHuvq/Ba/dT8sK+Y6kfRxSQ/mttRiSXt32LZru0PSnDbbvHY0JRqNYectb7expH/IbcHnJd0v6fDiSzM6BdS3hW22+e/RlGh0CqpzJ+R2wGrFuej5kqYXX5rRKaDOTVP0Q/0s7/O/NMb9MYPqJ295+3UlfSG/53lJP5d0fNM2h0lakl9fIumPii3F6A07b5J2kvSvir6OJOm0XmOpVUc7MB24GzgBWN38oqRtgNuAB4H9gZ2BzwJP97DvnYDNG5Zlwwl56DrmgN8sw+bAIXn9t9vtUNLGwM3AL4DfA44HPgN8amhRD9fQc9DgwKb3fneqwRaoWx7OBg4GPgLsCHwJ+BtJH2m3w/wd+g/gduAdwJnAuZIOG27oQzP0HDSo/DFBkoCrge2B2cRnugJYIGnDdjusUz0YNAcNqlIP6mAC+BqwJ/Eb/iLxOW3SsM2GRL0c19+ncTNB55xuAOxKHBt3Bd4PbAncIOk1I4+2OiboXleXAn8J7ALsRbRNb5C02WhDrZQJuucVAEkfINqsj4wywAqaoLecLuA3f+sOGmGMtSTpCGAe8GWi7XE7cL26X8js1u54tun1zVNKzw0x9FIVkTdJ04CbiLbg4cAOwBziuFwLBdW3Q5te2xr4P3o7p66MgurckcDfEu2rHYE/JY6r84Ydf1kKqnNnAMcS/VFvBb4OXCXpHcONvjwD5u1Sop/qGOL49cfAnQ373AO4HLgEmJn/vULS7kWUoQxF5I04D1pO9Bn393uQUqrlQnSez2la9y3gkj73MwEkYNOyyzSMHLTY5gJgaZdt/gJ4Cli/Yd1ngZWAyi7niHKwda4Hu5VdpmHlgeh4PL1p3a3AeR32cxawrGndN4BFZZdxhDmozTEBmJHL8vaGdesAjwFHrQ31YAo5qGw9qMtCXEB5CTikxWub5s9nouw4q7R0ymnDNm/Nud2l7HirsvSY141zXg8oO96qLO3yCrw5t1F3JE6Q5pYda1WWVjkF5gPXlR1b3RbgR8AFTeuWAWe22b5ru4PoHH667LJVMG/HAA8A65ZdvirlrcV7/iQfP7Ysu7zjnjvgPODWpnWnA3eXXd4xz9sjwAlN674DXFx2eUvM2x8AT3bJ2+XAzU3rFgCXll3ecc5b0/Z3A6f1Gk/dRrS3JWkdYuTyEsW0AI8rbpU8osdd3KG4ZfoWSfsVGOrI5FuTPkh0NHeyB/D9lFLjSNAbgS2IDujK6iMHk66U9Jik2/JoqSr7AXCIpC0BFLeszwRu6PCePYgRH41uBHbLo0GqZpAcTKrDMWG9/O/LI51SSr8GnidGWLZTp3owaA4m1aEeVNVGxEWRX5UdSI30ktON87/Oe+865lXSukRHz1PAT0cYV9W9Kq/5TotLgTNSSveUFViFtaure+X2732SLpD0xhJiq438nX8nr25L3UTcXdBJt3bH+pJWKKaiuK5mIz2Lytts4q73cyU9mqdVOK2CbdqWCq5vjY4Grk8pPTRgqGOnwNz9AJgp6d35/7MV8IfEHcOVV2De1qPhnC1bTW/nbGNvwLzNBn4MfCof95cppsFqnIao3bl7LaYsLDBvA1trOtqBNxKjNP6KSPgsoiF+iaT3dXjfKmJE92HE7VFLgVsk7VNsuCNxJHGw+pcu272JmDam0S8aXquyXnPwNDCXuJ3wIOAW4HJJHy42vEIdT5zQ/1zSGmIk98kppes6vKddXXgNMYK0agbJQZ2OCfcS06R8WdIminnKTgZ+l7hVr5061YNBc1CnelBV84jv76KyA6mRjjnNDdmvAtemlB4eZWAV1zKvkt4n6WnipPFEYFZKqfnYau21yuvpwC9TSv9YTkiV1yqnNxDTGrwX+DTwLuC7ktZ79dutR5sCv0XrtlS7c6te2h1LgT8jpvn6EHFsuU09PHunIorK27bElAHTiCklP0dMT3Hm0CIvV1F5e5mkGcC+9D54rSoKyV1K6TKiX+p7+Rx0BXAXcPJQoy9PUXXuRuCTknaQtI6kWbwyhVEdDJK3bYkLDW8ncvcJYjqU+Q3btDt3r3pf3qSi8jawtWl+zcmLCteklM7O//1TSbsRc2S27FhLKS0lvuSTFknamuh0/V4xoY7M0cDVKaXHe9g2Nf2tNuurpqccpJSeIDoXJt0haVPgJODiAuMr0nHAe4ir5yuAfYCvSFqeUuo0ortOdaHvHNTpmJBSWqOYV/2bwC+J2z0XANf38vamvytZDwbNQZ3qQRVJOptoHO2VUnqp7HjqoFtO80jhi4HXE8dM60GXvP4ncRfVpkR75NuS9kgprRpxmJXTKq+S9iWmzphZYmiV1a6u5g6hSXdJWky0mQ4GrhxtlLXTqi3Vsh3VS7sjpbSIhoskkm4nLpwcRwwuqYuh5o1Xpgw8Otf9xZLeAJwj6TMpzxtQA8POW6OjiY7Sf59ylONpqLnLv1efAz5OTHmxHXGh83Tg80OMu2zDrnMnEBdzluT9/Ay4EPjY0CIeDz3njTh+JeDIlNKTAJI+AdwoabOGARz97LOqisjbQNamEe1PEA/4WdK0/h6g2wMZmv2IeFhKZUmaCexGb1edH+XVV4Imbxmt7MirPnPQSmXrgaT1iVEaJ6WUrk0p3ZlSOg+4jPgxa6ddXXiR6KSsjCnkoJXK1oWU0uKU0kyi82zzlNKBwBvo/MCP2tQDGDgHrVS2HlSJpHOIkXr7p5QeKDueOuiW04bpON4GvDelVLnveRm65TWl9ExK6f6U0g9TSn8OrAGOGnWcVdMhr/sRo9pWSXpR0ovEfO1nSfIdGB30c1xNKT0CPIx/76biCeLCfqu2VD/nVh3bHbnT+I5O21RMUXlbBdzXdDH0HuJBeFW7U7OVQutbvtvto8CFKaUXBw1yTBWVuzOI+bG/kVK6K6V0FTHC/STV42HzheQtpfR4Smk2sCHxceA62wAABQdJREFU+/4WYuaBujy4eJC8rQJWTnYWZ5NT5032c7Y7d69sX16TovI2sLWmoz2l9AIxB88OTS/NIEZl9GMm8cFU2THEA6IW9LDtImBvSa9tWDeLeBjF8qFHNjr95KCVKteDaXlpHl33Ep2PC4uA329aNwu4I6W0ZnjhjcSgOWilynUBgJTSkymlx/PtxbsB13TYvE714GV95qCVyteDcSdpHjHl1/4ppXvLjqcOuuU0z1F7OdHJvl9K6dERh1hJA9bVdXjluRHWQpe8fo2opzMblkeAc4hpT6yFfutqvqPzd/Dv3cDyeeliou3UaBZwex+76tjukCTiO1GLz6rAvN0GbKd4ptukGcCzRAdOpY2gvs0mLkh8c6AAx1iBuduA1uegogaKrnMppedSSiuJGToOo/9ztrE0YN5uA7Zomlt8Rv53sp9zUZ/7rJQC8zaloGqzEHOwTzasnyVuu5kJbJVfnw28QHSwbkfc4rQGOLhhHxcBFzX8/cn8vu2BnYgRsAk4tOzyDpKDvM0GxBN2/7rNPs4Ebmn4+3XEVbDLgJ2JebCeAj5ddnlHmIOPEichOxIXa+bmunRi2eUdNA/AQuLpyRPANsTt1quB4xr20fx92AZ4Bvj7nIujch4OK7u8I8xBrY4JxJyU+xHzlL2fuPj0naZ91L0eDJKDStWDOizA+cRvz/7EiIXJZXrDNpvkz3Yifx5H5b/fVHb847h0yylxAnM1sBLYtWmb9cuOf1yXHvK6MTGabXdi1Mw7gX8mHsL8trLjH9ell2NAi/csB+aWHfu4Lj3U1enAV4gHqW2dj62LiBHtG5Udf5UX4Ijcdjoqt6XmESMz35xf77vdAZwKHJDbMzPzcWUN8K6yyzvmedsyfw/OJc7zDsh1/O/KLu84561h2wXAzWWXsUq5A07Lde6DxHnVLOB+ms4/qrwUlLfdib6obYG9iefmPQC8vuzylpi36cBDwBU5b+8h+jeuaNhmT+LO81OIuwBOIX4bdi+7vGOet3V5pQ/hfuDr+b+36xpP2QkZcnIn8pexeZnfsM0c4D6iM+1O4ENN+1gILGz4+6Sc1NXA/wDfBw4qu6xTzMHH8hdtizb7mA8sb1q3CzE31nPEVcVTAZVd3lHlgOhoX0J0Lj5F3Ib54bLLOpU8ECdSFxIdKKuJh0LObfxcm78Ped2+wE+IDoEHgWPLLusoc1C3YwIxZ+dDxA/TCuCLwLpN+6h7Peg7B1WrB3VY2nyGCTitYZs53bbx0ntOiY61dtvMKTv+cV16yOsGwFXEaOvn87/XUKMTnjLy2uY9y3FH+8A5BdYnHj73WMNv5Hxgy7Jjr8NCzM+8PB8HFgP7NLzWd7uDuHtjRd7fY/mz26Psco573vJ27yZGPq4m2rVfoKktWPWloLxtC/waOLzs8lUpd8RAhlOBZXm7h4i7sn677LKOed72JfpkniPuNrmINv05VV76yVtetwNwEzGgbCVxEX2jpm0+QPR1vEBMkVK7wWHDzhvtz4MWdotFeQdmZmZmZmZmZmZmZjaAtWaOdjMzMzMzMzMzMzOzIrij3czMzMzMzMzMzMxsCtzRbmZmZmZmZmZmZmY2Be5oNzMzMzMzMzMzMzObAne0m5mZmZmZmZmZmZlNgTvazczMzMzMzMzMzMymwB3tZmZmZmZmZmZmZmZT4I52MzMzMzMzMzMzM7MpcEe7mZmZmZmZmZmZmdkU/D8zOh6L2TQfSQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1490.4x331.2 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAbw0lEQVR4nO3deZzO5f7H8ffFjDEzGUuyzaQRUWQZUZYWlFJkSTnTKUsUnR+t2k4KJSTCqSjbsZ0S4ZgTUnHSYjmVSJT1UNFEjmKYsczM9/fHfbube8xy39Pcy1xez8djHua+vtf3uq+Pyz3v+V73wjiOIwAAbFMq1BMAACAQCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4IICMMbHGmG3GmM+NMZE52m80xmQbYwaGcn6AzQwf1QUEljEmSdJ6SRMcx3nKGFNF0mZJnzuO0zm0swPsRcABQWCMeUTSy5JulPSYpIaSGjuOcyikEwMsRsABQWCMMZKWSWonqYyk9o7jrArtrAC78RwcEASO6zfJuZKiJH1NuAGBR8ABQWCMqSZpoqSvJDU2xjwU4ikB1iPggABzb0/OlnRKUnu5gm6MMaZRSCcGWI7n4IAAM8YMlvSSpHaO43xsjCkj16sqoyQ1cxwnI6QTBCzFFRwQQO63CIySNNpxnI8lyXGcU5LulJQoaXzoZgfYjSs4AICVuIIDAFiJgAMAWImAAwBYiYADAFiJgAMAWCkiwOPzEk0AQGFMIAblCg4AYCUCDgBgJQIOAGAlAg4AYCUCDgBgJQIOAGAlAg4AYCUCDgBgJQIOAGAlAg4AYKVAf1QXAISllE37Qz0FuHVpEh+QcbmCAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADALdRj/RT16QELXt7plf78gWzNfjPt+iuaxuoR8s6evCO67V8wWw5jhOimcIXEaGeAACEg43rPtaOb77K89jxo0fUol0HXXTJZYoqG63Nn3+maS8+o1MnMtS11/1Bnil8RcABOOdlnj6tGWOH6a6BT2jS84+fdfyOex/0ut34qqv1S+o+fbR0IQEXxtiiBFCoPn36qFmzZlq2bJnq16+vmJgYdezYUYcPH9auXbvUtm1bxcbGqlmzZtq8ebPnvOzsbL344ouqU6eOoqKiVLduXc2ePdtr7GXLlql9+/aqUqWK4uLi1KJFC33wwQdefYYPH67KlStr48aNatGihWJiYpSUlKRPP/20WOpbOm+GykSV1fVd/uTzOeXKV1Tm6dPFcv8IDAIOgE9++OEHDR06VC+88IKmTp2qtWvXqn///kpOTlZycrIWLlyozMxMJScne56beuCBB/TCCy+of//+WrZsmbp166a+fftq6dKlnnH37NmjW2+9VXPnztWiRYvUqlUr3XzzzVqzZo3X/aenp6t3794aMGCAFi1apKioKHXr1k3p6emePtnZ2crMzCzwKysry2vcXw8d1IJpf1O/x4arVKmCfyRmZWYqI/24Nnz2b320dJFu7tH7j/61IoDYogTgk8OHD2vdunWqXbu2JGnz5s0aO3asZs+erV69ekmSHMdRx44dtW3bNkVGRur111/XzJkz1bu3KwhuuOEGpaam6rnnnlOnTp0kSYMGDfLcR3Z2ttq2bautW7dqxowZat26tedYRkaGJk6cqHbt2kmSqlevrqSkJH3yySfq0KGDJKlv375nXSHmdt1112n16tWe27MnjlSTltepwRUtCjzv10MHdU/7pp7bd9z7kDrd2bfAcxBaBBwAnyQmJnrCTZLq1KkjSZ7Aydm2f/9+7d69W6VKlVK3bt2UmZnp6XP99ddr3rx5ysrKUunSpbVv3z4NGTJEK1euVGpqqufqL2e4SVJkZKTatGnjuV2/fn1J0r59+zxtw4cP9wrMvJQrV87z/bavN2jtqmV6bdHqQuuPq1BJ4/6xTBkZx7Xly3VaPHOyomNidNs9Aws9F6FBwAHwSYUKFbxulylT5qz2M20nTpzQoUOHlJWVpfLly+c5XmpqqmrUqKHOnTsrLS1Nzz//vOrUqaPY2FgNHTpUBw8e9OofFxfntYWY877OqFmzphISEgqswxjj+X7GuGG6qfvdiilXTsfSjnjaT508oeNpRxVbLs7TVjoiQnUaNJYkNWzWSqVMKb09Zbw6JvdVVHR0gfeJ0CDgAAREpUqVFBERoTVr1uT53FaVKlW0a9cubdy4Ue+9955nm1FybUcWhb9blPv3/lc7t2zSu29O9+oze+JIzX31RS3+8vt8x7n4soY6dfKkDv/ys6rXrFWk+SKwCDgAAdGuXTtlZWXpyJEjat++fZ59zgRZVFSUp+3777/XmjVr1KhRI7/v098tymdemaWsrEyv48/e10Od7uyrFu1uLnCcbZu+UGSZKFW6oJrf80RwEHAAAqJevXq6//77lZycrCeeeELNmjXTiRMntHXrVu3YsUPTp0/XpZdeqoSEBA0ePFgjRoxQWlqahg0bpvj4+CLdZ2JiohITE33uXz/pyjzbq9espcubtfTcfuyujmp76+2KT6ytzNOn9fV/PtXy+bPU5e7+bE+GMQIOQMBMmjRJdevW1bRp0zR06FDFxcWpfv366tevnyTXldvixYs1cOBA3X777UpISNCQIUO0evVqbdmyJcSz/12tevW1bN7fdehgqqLKRqvGhbX0wPCXdd0tt4V6aiiACfBnqfFBbQDCUsqm/aGeAty6NIk3hffyH2/0BgBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWCki1BMAziUpm/aHegrAOYMrOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOACAlQg4AICVCDgAgJUIOKCEWf/RCj3U4wbdfuXF6t+xpVLmTi2w//Sxw9Q1KUEzx4/wat+3Z5ce79lJf77mMo178v+UkX7c6/jWDevV98YrzmrPy6p/LVDXpIQ8+85742X1bNvQc/vATz+qa1KC5yu5dT0N/vMt+uyDd/Pt86eWl2jAra01/ulB2vrVfwqdDyARcECJ8t2mLzRm8H26pEETDfnbTF3f5U+a88oo/evN6Xn2/3H3Dq1Kma+Y88qddeyVYY+o+oWJenzM6/pxzw4tnPGq51h2drZmjBuung/8VdExsQGppc8jz2rM7BQ9OW6qatSspXFP/kVffLIyzz7PvjpHPe57SGlHftWQft319pTxAZkT7BIR6gkA8N38qRN1WZMrNWjYOElSUsvrdPzoES2YOkE39+ilyMgyXv2nvTRUne7sp9XLFnm1Z6Qf145vNmrIxFkqX+l8HU87qiVzpqjnA09JklYumafSERFq06l7wGqJT6yteo2ukCQ1vuoa7d72jVa8M0fNr70hzz6Xq6Wu79xDb00eq7ffGK8GV7RQw2atAjY/lHxcwQElyJ7tW9X4qqu92pq0vFbHjh7R9q83eLWv/XCp9u3Zqe59B541TubpU5KkMmXLuv+M9rSlH0vTW5PH6d7Hn5MxJhBlnKVUqVKqVa+BDv60r9C+fxrwqCpdUFXvv/OPIMwMJRkBB5Qgp0+dVERkpFdbZBnXVdu+Pbs8bSdPZGjm+BHq+eBfVTY65qxxypWvqKrxNbVs3kylHflVHyx+U7XrN5IkLZg2UY2vutpz5eSP7KwsZWVmen052dk+nXvwpx9VofIFhfYrXbq0Gl7ZWtu/+crv+eHcwhYlUIJUuzBRO7d+7dW2c8smSVLa0d88bYv+/poqXlBFbTrmv8XY/6kXNPaJ+/WP115UjZq1NOCpkUr9YY9WLnlbExd8WKT53XVt/Tzby1WoeFabk52trMxMpR9P04f/nKedWzbpvidH5HH22SpXqa4jh38p0hxx7iDggBKkw+13641RT+uDxW+q1Q0dtXPLJs+rKEuXcm3IHNj/g5bMnaIRU+YXuMV4xdXtNHvV1zp0MFXVEi5S6dKlNfKhe3TrXfepctUaWj5/lhbPmixJuu2egbqlR+9C5zdyxiJFRZX1ant/8Zta/+/3zuo76pG+nu8jIiLV+e7+6nBHr8L/EiQ5juNTP5zbCDigBLm+S7L27vhWb4x6WpNHPKmostHq9dDTmjbmWZU/37W9N+eV0Wraqq3ia9XRsbQjkiTHydbpUyd1LO2IYs+L8wRfVHS04i+6WJL09fpPtXfHt3pszGTt2f6t3po8Vi/OSpEkPdm7i+o3aa7EunlfoZ1x8aWXn/Wqyy8+XZln376PDVf9Js0VHXueqsRfeNYLZAryv19+VvlKhW9n4txGwAElSOnSpdX/qZH68/89rkMHUlU1vqb2u597q9ewqSRp/97d2rvj27OumpbPn6Xl82dp+orPVblqDa9jWZmZmjFuuHo/PERRZaO1ZcNaNWzeWgm16kiSGl3ZWls2rC804PxR/cJE1WnQ2O/zsjIz9c3na9SgaYtimwvsRMABJdB5cRV0XlwFSdJ778zRpY2becJo0NCxysjwfsP1y08NVIMrWqjDHT1VvuL5Z423YuFcxcaV19U3dfa0nTyR8fv3GelSmGwLzp86QYd/OaCbbr871FNBmCPggBJk++YN+m7TF6pVt4HSj6fp0xUp2rjuY43++2JPn7yuiiLLRKly1Rp5vm/s2NHfNH/KBA2b/KanrUHTFpo9cZRWLnlbkrT5i7Xq+eDTAaioYPv37lZchYrKPH1aB/b/qM/eT9FXa1cr+f5HdXmzlkGfD0oWAg4oQUpHROqz99/V22+MlylVSvWTrtTomf9U4iWXFXnMea+/rObXtVfty37/OK2LL71cvR9+Wm9OekmS1OeRZ1SrXvFtT/pq1gTXqyrLREWpYuWqqteoqUbOWKQGTa8K+lxQ8pgAvxopPPY0gDCRsml/qKcAhJ0uTeID8okCvNEbAGAlAg4AYCUCDgBgJQIOAGAlAg4AYCUCDgBgJQIOAGAlAg4AYCUCDgBgJQIOAGAlAg4AYCUCDgBgJQIOAGAlAg4AYCUCDgBgJQIOAGAlAg4AYCUCDgBgJQIOAGCliFBPAIGXsml/qKcAAEHHFRwAwEoEHADASgQcAMBKBBwAwEoEHADASgQcAMBKBBwAwEoEHADASgQcAMBKBBwAwEoEHADASgQcAMBKBBwAwEoEHADASgQcAMBKBBwAwEoEHADASgQcAMBKBBwAwEoEHADASgQcAMBKBBwAwEoEHADASgQcwtL/DqYquVVddU1KUEb6cUnS4V8OaNaEEXq4R3slt6qrfh2a62/PPqzDB38O8WyBwMvrMZHb9LHD1DUpQTPHjwjy7MITAYewNGvCCyobE+vVtvu7b7T+oxW6pkMXDfnbTPV5+Bnt2LJRT/bpmu8DHrBFXo+JnH7cvUOrUuYr5rxyQZxVeCPgEHa2fvUfbVy7Wl17DfBqvyypuSYt/ljd+w5Sw+atdU2HLnp64kz9krpP61YuD81kgSDI7zGR07SXhqrTnf0UW6588CYW5gg4hJWsrCxNG/OsevR/ROUqVPI6dl658iodEeHVFn/RxYoqG60jvx4K5jSBoCnoMXHG2g+Xat+enered2CQZxfeCDiElfcXztXpUyd1S4/ePvXfu+NbnTyRoQsvrhvgmQGhUdhj4uSJDM0cP0I9H/yrykbHBHl24Y2AQ9g4+tuvemvyWN3z6FBFREYW2j87O1vTxw5TjZq11KTFtUGYIRBcvjwmFv39NVW8oIradOwe5NmFv4jCuwDB8eZrY1T38iQ1u+Z6n/rPfXW0tm/+SiOnv+NTIAIlTWGPiQP7f9CSuVM0Ysp8GWOCPLvwR8AhLPywe7tWpczXyBkLdSztiCTX1oskpR87qlKlSimqbLSn//IFs7Vk9hsaPHqS6jZsGpI5A4Hky2Niziuj1bRVW8XXquPp4zjZOn3qpI6lHVHseXHndPAZx3ECOX5AB4dvUjbtD/UUCrX+oxV68dF78z1+Q9dkDRo2TpK0duUyjXvyL+r10NPq2uv+YE0RCCpfHhO7vt2svTu+zbfP9BWfq3LVGoGYXrHq0iQ+IClMwJ0DSkLAHf31sL7fvc2rbeOa1Vo8a7KefXWOqiVcpPjE2vrmy7V6fmBP3dT9Lt37xPMhmi0QeL48JjKOH1NGhvd7QF9+aqAaXNFCHe7oqUsbXaHIMlHBnHaRBCrg2KJEWIirWEkNm7Xyajv40z5JUv2mVyk6JlY//nenRj9yr+ITa+vqGztr++YNOc4/X9UvTAzmlIGA8uUxkZfIMlGqXLXGWeeeiwg4lBg7tmxU+rGj2rvjWz11T1evY21vvUMPPT8hRDMDEI7YojwHlIQtSgDnrkBtUfI+OACAlQK6RcmVAwAgVLiCAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiY/qQomwZM4UvT1lvE6kHy+8M3COKBsTq+QBj6prrwGhnkpY4goOJULK3KmEG5DLifTjSpk7NdTTCFsEHEqELj37q2w+n54OnKvKxsSqS8/+oZ5G2Arohy2nbNrPhy0DAArEhy0DAOAHAg4AYKWAvopyU8r03ZJ+C+R9BEENST+FehJ/kA01SHbUYUMNkh11UEOY2JSiXcOGDatT7AM7jhOwr+HDhzuBHD8YX9QQPl821GFDDbbUQQ3h8xWoOtiiBABYKdAB91yAxw8GaggfNtRhQw2SHXVQQ/gISB0BfZsAAAChwhYlAMBKBBwAwEoEHADASn4HnDHmaWOMY4x5LUebMcYMN8b8ZIzJMMasNsY0KGScPu5xcn+VLUohxVDDbcaY940xv7iPtfFxrOuMMRuMMSeMMf81xtwfsIl732+x1GCMaZPPOlwa0AJ+v3+vOowxkcaYMcaYzcaY48aYVGPMW8aYmj6MFRZrUdQawm0t3G0jjDHb3HX8aoxZZYxp5cNYYbEW7ja/awjlWuRVQ67jU93HH/NhrJCsg/u+i6WOP7IWfgWcMaaFpPskbc516AlJgyU9IKm5pIOSPjTGlCtkyHRJ1XN+OY5zwp85+auAGmIlrZX0qB9j1ZK03H1ekqTRkl41xnQvntnme7/FVkMODeS9Fjv/yBx9kU8dMZKaShrp/rOLpAslrTDG5PvBBGG2FkWqIYdwWQtJ2i5poKSGkq6WtEeuOqoWMFY4rYVUhBpyCOpaFFDDmeO3y/UzttA3d4dqHdz3XWx15OD/Wvj6hjlJ5SXtltRO0mpJr7nbjaRUSUNy9I2WlCZpQAHj9ZF0LJhvJsyvhlx9KktyJLXxYbwxknbmapsuaV0JqqGNu2/lcFuLHH3ru+fYsKSthZ81lIS1iHPP8aYSvBa+1BD0tSisBkkXSdov6TJJeyU9Vsh4QV+HANVR5LXw5wpuqqSFjuP8O1d7LUnVJH1wpsFxnAxJn0gqbCsj2hjzvTFmnzFmqTEmyY/5FEV+NRRVS+Wo2+19Sc2MMZHFdB+5FXcNZ3zp3kpbZYxpW8xj58WfOuLcf/5aQJ9wXwtfajgjLNfCGFNGUn9JRyVtKqBr2K6FHzWcEcy1yLcG95X/PEkvOI7znY/jhWIdpOKv4wy/18KngDPG3CepjqRn8zhczf3ngVztB3Icy8t2SX3l2r65U9IJSWuMMZf4Mid/FVJDUVVT3nVHyHUVVawCVEOqpL9I6i7pNrnWZZUx5tpivA8v/tTh/oH0sqR3HcfZV0DXsF0LP2oIy7UwxnQyxhyT6zH6iKT2juPk/rvOKezWogg1BHUtfKjhOUn/cxzndT+GDeo6SAGro8hrUejzAcaYepJGSbrGcZxTBXTN/Y5xk0fb750dZ52kdTnuZ61cv1E9IOnBwublDz9qKIq86s6r/Q8JVA2O42yX6x/MGeuMMYmSHpPrKrxY+VOH+7e9f0iqIKmzD8OH3Vr4U0MYr8VHkprI9UPxPkkLjDEtHcdJLeCccFsLv2oI5loUVoMx5jq5ntJpUoThg7IOUuDq+ENr4cN+ah+5/jIyc3w5krLd39dz326e67xlkmb7uXc7U9J7AdgTLqyGqBx9/Xn+6hNJk3K13SHptKTIklBDPvc1TNJ3xb0O/tQh1y9f70jaJqlaSVwLf2sI17XI47ydkp4tSWvhbw3BXAsfahid4/ucx7Mk7QuHdQhkHX9kLXx5RdcSSV/mapvp/gcyStIOST9Lai/pC0kyrpf6XyPpcR/Gl/scI6mRpK99PccPhdVQ1CuidZK65mprL+lLx3FOF3HM/ASqhrw0kWtbIBAKrcP9/MDbki6XK6R/9mHcsFqLItaQl5CuRT7nlZIUVcC4YbUW+ZxXWA15CdRaFFbDL5LezHX8fbmey5pWwLjBXAcpcHXkxae1KDTgHMf5Tbn+TzdjzHFJhx3H2eK+PVHSEGPMNrkC7xlJxyS9leOcVZI+dxznr+7bwyStl6v4OLm2JRvJtddarHysoZKkmnJtJUlSHWPMb5J+PvPDyRgzxz1eL3efNyQNctc/RVJruX6LubOk1GCMeViuVzJtlVRG0t1yPSgC8lLiwupwb+m9I9dLiG+V5BhjzjyXe8RxvYAprNeiqDWE4VrEGWOekPSuXD9MLpDr5fYJkhbkOCec16JINQRzLXx5bCvXc2nGmNNyPa6352gL2Tq47zcgdfyRtSiu//D0JbneGjBJUkVJ/5F0o+M4aTn61Jb0Y47bFeR6tU01SUckbZR0reM4nxfTnPzVWa7fNs448xvFc5KGu7/3eqOu4zh7jDG3SJogVzD/JOlBx3EWBXaq+fK7Brn+wYyTFC8pQ65/RB0dx1keuGkWKEGuFx5J0oZcx+6RNMv9fTivRZFqUPitRaZc7z3qK+l8Sf+Ta5fmWsdxcr6/KZzXokg1KPzWwhfhvA7+KLa14H8TAABYic+iBABYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWImAAwBYiYADAFiJgAMAWOn/Ae7K9FodHFgsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import arviz as az\n",
    "\n",
    "infer_data = az.from_numpyro(mcmc)\n",
    "az.plot_posterior(infer_data, kind='hist');\n",
    "az.plot_posterior((infer_data.posterior.tau.values * len(count_data)).astype(int));"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}