PPC Difference

ppc_difference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Posterior Predictive from Inference Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc3 as pm\n",
    "from scipy import stats\n",
    "import arviz as az\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'3.9.3'"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pm.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.10.0'"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = stats.norm().rvs(10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [mu, sigma]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:01<00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds.\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [mu, sigma]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:01<00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model:\n",
    "    sigma = pm.HalfNormal(\"sigma\", 10)\n",
    "    mu = pm.Normal(\"mu\", 1, 100)\n",
    "    \n",
    "    trace = pm.Normal(\"y\", mu, sigma, observed=y)\n",
    "    \n",
    "    trace = pm.sample(random_seed=0)\n",
    "    inf_data = pm.sample(random_seed=0, return_inferencedata=True)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [4000/4000 00:03<00:00]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [4000/4000 00:25<00:00]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with model:\n",
    "    trace_ppc = pm.sample_posterior_predictive(trace, var_names=[\"y\", \"mu\"], random_seed=0)\n",
    "    inf_ppc = pm.sample_posterior_predictive(inf_data, var_names=[\"y\", \"mu\"], random_seed=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f7ba36163d0>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD8CAYAAACINTRsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3xb1f3/8deRZFnee29n7xAcIDsQwt6jjJZZSIGWQksLpYPRlkL58it0UGhoC3QEKKOEACkhhCSQhCwynTiJY8fx3ntqnN8fko3jOLETpEiWP8/HQw/bd0gfXdtvH5977rlKa40QQgj/YPB2AUIIIdxHQl0IIfyIhLoQQvgRCXUhhPAjEupCCOFHTN4uIDY2VmdmZnq7DCGEGFK2bt1ao7WO67vc66GemZnJli1bvF2GEEIMKUqpov6WS/eLEEL4EQl1IYTwIxLqQgjhRyTUhRDCj0ioCyGEH5FQF0IIPzJgqCulHlNK6T6Pil7rlWubMqVUu1JqtVJqgmfLFkII0Z/BttT3AUm9HpN6rXsQeAC4F5gOVAEfK6XC3FinEENOU4eVl9YWkF/V7O1SxDAy2FC3aa0rej2qwdlKB+4HntJav6213g3cAoQBN3qmZCGGhkeX5vLEh3u5599fYnfIfQvEqTHYUM92da8UKqVeV0plu5ZnAYnAiu4NtdbtwFpgpntLFWLoaO208f7OMoLNRvZXtrCnrMnbJYlhYjChvhG4FbgAuBNniK9XSsW4Pgeo7LNPZa91R1FKLVJKbVFKbamurj7hooXwdZsO1WG1a564ciIAn+fXeLkiMVwMOPeL1np576+VUl8ABTi7Wb44mRfVWi8GFgPk5OTI/6XC72w73IBBwfkTEsmKzWdnSYO3SxLDxAkPadRatwC5wCigexRMQp/NEnqtE2LY2VPWRHZcKMFmE+OTwsmV7hdxipxwqCulLMBYoBwoxBneC/usnwOsd1ONQgw5e8ubGJ8UDsC4pDAO17XR2mnzclViOBjMOPVnlFLzlFJZSqkzgbeAEOBVrbUGngMeUkpdpZSaCLwCtABLPFi3ED6rw2qntKGdEXGhAGTEhABQXN/mzbLEMDGY+dRTgdeAWKAaZz/6WVrr7rl8nwaCgOeBKJwnVs/TWsvgXDEsFdU6wzszNhiA9Gjnx+K6dsYmhnutLjE8DOZE6fUDrNfAY66HEMPeodpWALJinS30NFeoH66TlrrwPJn7RQg3O1TjDPXubpeo4ABCA00US6iLU0BCXQg3O1TbSnSImYigAACUUqRFB0uoi1NCQl0INyusaSUzJviIZWlRQdL9Ik4JCXUh3OxwbRuZrq6XbunRwRTXt+E8BSWE50ioC+FGNruDiqYOkiODjlieEhVEh9VBfZvVS5WJ4UJCXQg3qmruxKEhKdJyxPKkCOfX5Y3t3ihLDCMS6kK4UXdoJ0cc2VJPdH1d0dhxymsSw4uEuhBuVNbgDO2+3S+J4d0tdQl14VkS6kK4UXdLvW/3S1xYIEaDkpa68DgJdSHcqKyhg9BAE+GWgCOWGw2K+LBAKpok1IVnSagL4Ublje09J0X7SoywSEtdeJyEuhBuVN7YQVKf/vRuSREWGf0iPE5CXQg3KmvoIPlYLfXwIMobO+QCJOFREupCuEmnzU5NSydJEcduqbd12WmWm2UID5JQF8JNuvvL+4586ZboasFLv7rwJAl1Idyke4x6ynH61EHGqgvPklAXwk16xqgfZ/QLQIWcLBUeJKEuhJt0t8CP1aceHyYtdeF5EupCuElZQztRwQEEmY39rjebDMSGBkqfuvAoCXUh3KS8seOYrfRuzrHqEurCcyTUhXCTsoZ2ko8x8qWbXFUqPE1CXQg3KW1oP2p2xr6SIiwy/4vwKAl1IdygqcNKc4ftmMMZuyVGWGhst9LWJRcgCc+QUBfCDcoanMMUU6IGbqmDXIAkPEdCXQg3KK13hfpALfXw498B6Z9fFHHva9uolC4acZJM3i5ACH9Q2jC4UD/eVaVbi+r4xbu7Aei02ll8c46bqxTDgbTUhXCD0vp2zEbnOPTj6bmqtJ+W+D83FBFuMfHt2Vl8vLeyp0tHiBMhoS6EG5S6hjMaDOqodb2n2rUEGIkKDjhqXnWHQ/PZgRrOHZfAjWemozWsyqvyeN3C/5xwqCulHlZKaaXUn3otU0qpx5RSZUqpdqXUaqXUBPeWKoRn/XdbCb/5cC+HalpPeN/ShvYjTpIerG7hZ//dxaynVpH90w+Z/sRKfvL2Tg5Wt5AQfvRY9X2VzdS2djFjRAzZsSEkRVhYf7Dma78nMfycUJ+6UuosYBGws8+qB4EHgFuBfcAjwMdKqTFa62Y31CmER/1jwyEeWZoLwNtbS1j5w3lEhZgHvX9pfTvzRsfhcGheWHOQ51bux6AUC8bFc8VpyRTXtfPejjLe+bKUjJjgo/rU1+U7A3zWyFiUUuRkRvNlUb3b3p8YPgYd6kqpCODfwO3Ao72WK+B+4Cmt9duuZbcAVcCNwF/cWbAQ7tbUYeXJD/M4e0wcD5w3hsufX8f/rdjHb66cNKj9O212qpo7iQ8P5LtLvmT57gounpzEY5dOIC7sqz726uZOHn5nFyv3VmIJOPKf5PUHa8mKDem5eGlCcjjLdpTR0NZFZPDg/7gIcSLdL4uBt7TWn/ZZngUkAiu6F2it24G1wMyvXaEQHvb+jnLarXbuO3c0E1Mi+NaZ6by+6fCgT1SWu+ZRf39HGf/LreDnF4/jTzecdkSgA8SFBfKXm05nSmoEHVYHf/usEACr3cHGglpmjojp2XZ8UjgAe8qb3PEWxTAyqFBXSt0JjAR+3s/qRNfHyj7LK3ut6/t8i5RSW5RSW6qrqwdbqxAe8ebWYkYnhDIlNQKAO+Zko4HXNx0e1P4F1c4++OL6dl745jTumJON8x/YoxkNirvnjwDgNx/uYUdxAztLGmntsjNrZGzPduOTXaFeJqEuTsyAoa6UGgP8BrhRa211x4tqrRdrrXO01jlxcXHueEohTkpVUwfbDjdw+dSUniBOiw7m7DHxvLa5GKvdcdz9O6x2fvPhXgAeuWQ8F0xMGvA1RyeEARBqMfGD/2xnzT5nw+as7K9a6rGhgSSEB5IroS5O0GBa6jOAWCBXKWVTStmAecA9rs9rXdsl9NkvAahwW6VCeMDaA84TlGePiT9i+U1nZVDd3MnHe/r+A/qVLpuDe/79JfnVLZgMiptnZA7qNdOjg7EEGJgxIpaC6lbe+bKE8UnhRPc5MTsmMZwDVTLOQJyYwYT6u8AkYGqvxxbgddfn+3GG98LuHZRSFmAOsN7N9QrhVmv2VxMXFsi4pLAjls8bHUd8WCAf7Czvdz+r3cF9r29jVV4VYxPDGBEX2u8Y9f6YjAYmJkdQ3dzJpZOTKGloZ1JK+FHbZceGUFjdesQ4dyEGMmCoa60btNa7ez+AVqDO9bUGngMeUkpdpZSaCLwCtABLPFm8EF+H3aH57EA180bHHdUHbjAoFo5P4NN9VXRY7Ues6w705bsr+MUl47HaHWTFhpzQa0/LiGJXSSNT0iIB2F/ZctQ2WbEhtHbZqW7uPMF3JoYzd11R+jTwLPA8zlZ8EnCejFEXvmxnSQMNbVbmje7/vM4FExNp67Lz2YGvLgKy2R3c//p2PtzlHOVyy4wMDte1kRV3YqF+7rgEuuwOfv3BXkIDTWwrbmDDwdojtun+Q1FwEhdDieHrpEJdaz1fa/29Xl9rrfVjWuskrbVFaz3P1aIXwmet2V+NQcHsXqNOejsrO4Zwi4mPcp2nhmx2B/e9sZ0PdpXzs4vGccecbEob2rHa9Qm31HMyopiW7mylP3zRWBLDLTyzYt8RXS3dz1kooS5OgMzSKIatNfurmZIWecwrRwOMBs4dl8DKvZW0d9n48Vs7+WBnOT+9aCx3zs0GnNMBgLP/+0QYDIpXbz+DyqZORsaHolD89L+7+HRfFeeMdY45SI4MwmwynNS0BWL4kgm9xLBU39rFjuKGY3a9dDtvQiINbVZue2UL7+8s5ycXjmXR3BE963NLnUMOxySGHespjinMEsDI+FAArs1JJSMmmP+3Yn9Pa91oUGTGBEv3izghEupiWFp7oBqHhrkDhPqsETEYleKLgloeumAsd80bccT63LImMmOCCbMEfK16AowGvnv2SHLLmnqGWYKzC0a6X8SJkFAXw9L7O8tJCA9kSmrkMbexOzQ/X7obu9aEBpr4jqvLpbfc8kYmJEe4paYrpqaQGG7hhdX5PcuyYkMpqm3F7pBhjWJwJNTFsNPYbmXNvmoumZyM8Rhjy+0OzY/e3MHS7WVcNCmRlk4bO0sbj3qe4rr2nkv6vy6zycAdc7L4oqCOLw87Z2jMjg3BateU1Le55TWE/5NQF36vrKGdfRVfja5dur2ULruDy6Yk97u9w6H56Tu7+O+2Un503mievHIyJoPqGQXTbbcr5Ce4KdQBbjgjncjgAF5YfRCAjJhgAIpqJdTF4MjoF+HXPtlbyaJ/bsXu0Jw7Lp6bZ2Tyx1X5TEuPZHLq0d0mWmt+9cEe3thSzPfPGcn3zhkFOIc3Lt9VzoPnj+m5UGntgWoCjM65z90lJNDELTMy+f0nBzhQ2UxGjHNUTVGdhLoYHGmpC79ltTt4ZGkuo+JD+eHC0azLr+Xmv2+irdPGr6+Y1O9Mis+uPMDL6w5x26xMfrBwdM/yK09L4VBtG5sK63qWrdlXTU5GNKGB7m0b3TIzk6AAIy+sOUh8WCBmk4FiCXUxSBLqwm+t2VdNaUM7P1w4mu8vGMWaH8/nz9+cxkc/mNtvP/iSjYf5wycH+EZOKr+4ePwRoX/RpCTCAk28vrkYgMqmDvIqmgccPXMyokPMXDc9jfe2l1HZ3EF6dDBFtTICRgyOhLrwW//LrSDcYmK+awbG+HALF01KIjUq+Khtvyio5ZGlu5k3Oo4nr5p81ORcQWYjV5+eyrIdZRTVtvKaa671heP7Tk7qHt+enYUGXl53iPToYA7XDe6GHUJIqAu/tT6/hjmj4jCbjv9jXlzXxt3/2kp6TDB/uOG0Y46IuWveCIICjFy/+Av+/OlBzhuf0HPxkLulRQdz8aQklmw8TGKEhcO1MlujGBwJdeGXSurbKGvs4Iys45/E7J4T3e7Q/PXmHCKCjn0RUWKEhb/ekkN8uIVZI2N4YpD3MD1Zi+Zm09Jpo6qpg9YuO3WtXR59PeEfZPSL8EubDzlPaE4fYGTK/32Ux67SRv5y0+lkxw3c6j4zO4al353llhoHMjElglkjY9jiei9FdW3EhAYOsJcY7qSlLvzSpsJ6wiym487JsmZ/NS99VshNZ2Vw/oR+b6frdYvmjqCh3QYgI2DEoEioC7+UW9bIlNTIY/aPVzd38sB/tjMmIYyfXTzuFFc3eHNHxTI6wfkfhMzWKAZDQl34HYdDc6CyhVEJ/XenOByaB97cQXOHjT/eeBqWAOMprnDwlFLcPd85iVjvMfJCHIuEuvA7pQ3ttFvtjE7ov+vlb58XsnZ/Nb+4ZPwxt/Ell0xOxmxU7CxpHHhjMexJqAu/c6DKOc/L6H5a6jtLGnj6ozzOn5DAN89MP9WlnZQAo4GxSeE0d9rYXtzg7XKEj5NQF36n+ybOI+OPbIU3d1j5/mvbiA0N5LdXT+53mgBf1X0zj97T8grRHwl14Xf2VzaTEB54xJhzrTUPv7OLw3Vt/P7604gM7v8Wdr6q+yKnFbmVMmWAOC4JdeF3DlS2HNVXvmTTYd7fWc4D540Z8IIkX5Qe7ZzawGhQ/PWzQi9XI3yZhLrwKw6HJr+q5YjL93PLGnl82R7mjo7j7j63oxsqukN9Yko4b24tlqtLxTFJqAu/0nfkS0unje8t2UZUcAC/+8aUoybqGiqiQ8yEBprIig2hw+rgHxsOebsk4aMk1IVf2V/pHPkyKj4UrZ13MCqqbeUP159G7BC+xF4pRVp0MA1tVhaMjefV9Ydo77J7uyzhgyTUhV/pHvkyKiGM1zYV896OMn64cDRnZsd4ubKvLzMmmKK6NhbNzaa+zcpbW4u9XZLwQRLqwq8cqHKOfCmtb+exZbnMGRXLPfNHersst0iPCaakrp3TM6KYmhbJXz8vxO6Q6XjFkSTUhV85UNlCVmwI31vyJZFBATx73dQh24/eV0Z0CF12B5XNnXxnbjZFtW1H3QxbCAl14Te6R75UNHVwqLaVP9wwtPvR+8qMcY6AKapp5bwJiWTGBPOXNQfl5hniCBLqwm+U1DtHvhyqaeNH54/hLD/oR+8tvTvU69owGhTfnpPNjpJGNspEX6KXAUNdKfVdpdROpVST67FBKXVxr/VKKfWYUqpMKdWulFqtlJrg2bKFONqHu8sByMmI4q65Q3M8+vEkRQRhNho45Lqi9NrTU4kOMbN4bYGXKxO+ZDAt9RLgIWAakAOsAt5VSk12rX8QeAC4F5gOVAEfK6V8f/o74TeaO6z8+VPnvCh/vOE0v+lH781oUKRGB3G41nmzDEuAkVtmZLIqr6pnKKcQA4a61nqp1nq51jpfa71fa/0zoBmYoZwzIt0PPKW1fltrvRu4BQgDbvRo5UL08stle2jqsJEUYSEpMsjb5XhMRnQwRbVf3QHpphkZWAIM0loXPU6oT10pZVRKXQ+EAuuBLCARWNG9jda6HVgLzDzO8yxSSm1RSm2prq4+qcKF6LYit4I3t5YQbjFxWnqkt8vxqIyYEIpqW3tOjkaHmPlGThpLt5dS3tju5eqELxhUqCulJimlWoBO4EXgSq31LpyBDlDZZ5fKXuuOorVerLXO0VrnxMXFnUTZQjjVtnTy8Du7GJsYRlOHjXGJ4d4uyaMyYoJp7bJT22vulzvnZKM1vLD6oBcrE75isC31fcBU4EzgBeBVpdREj1UlxCA9uTyPpg4rN8/IBGBSaoR3C/KwjO4RML2m302LDubanFRe31Tc01p3ODQOuTBpWBpUqGutu1x96lu11g8D24EfAN1XPiT02SWh1zohPGLLoTre2lrCHXOyKW1wDvPLyRx60+qeiIyYEIAj+tUB7pk/EofWvLD6IIU1rcz9v0+Z+ssVbDhY640yhRed7Dh1AxAIFOIM74XdK5RSFmAOzj53ITzC7tD8/N3dJEdYuPeckWwsqGNSSgShgSZvl+ZRqVFBKHV0qPdurT/41g6a2q0EmY08+PYOmUpgmBnwN0Ap9RTwAVDMV6Na5gMXa621Uuo54KdKqTxgP/BzoAVY4qmixfB0oLKZv31eSGxoIIkRgeRVNPP8jdMA2FHSwLdnZ3u5Qs8LNBlJjgiisOboux/dM38k/9lczOZD9Tx0wVgyY4K5+99f8mleFeeO7/vPtPBXg2nWJAL/cn1sBHYCF2qtP3KtfxoIAp4HooCNwHlaaxk4K9ymqrmDa17cQJfNQbvVjsmgmJAUxkWTElm2sxyrXffcx9PfjUoIJb+q5ajladHBjEsKZ3dZE3NGxTI2MYyE8ED+8UWRhPowMphx6rdqrTO01oFa63it9bm9Ah3t9JjWOklrbdFaz3ONVxfCbZ74YC/tXXaW3Tubb52Vjs2h6bJr2q12Xl5XSEpkEGcOwdvUnYyRcaEcrG7pv1vFdc3Vf7YUYzIauHpaKuvya2hss57aIoXXyNwvwucV1rSydHsZd8zJIjHCwvJdFYxJDCW/uoWcX69k2+EGfrBwtF9eRdqfUQmhdNoclNQf2a/e0mkjr7yZCcnhvLbpMIdr2zh3fAJ2h2bNAbkeZLiQUBc+779flmBQcOvMTF5df4ja1i5+e/UUFt+Uw+yRsTx51SSuOT3V22WeMiPjnTNw9O2C2VxYh82huWf+CAxK8ezK/UxJjSQmxMwne/teSiL8lX8PFRB+4eO9VeRkRhMeFMDL6wqZPyaOqWnOK0cXDsO+4u6bah+oamHBuK/e/7r8GswmAwvGJXDbrCz+svYgi+ZmM39MPCv3VmKzOzAZpR3n7+Q7LHxacV0be8ubWDgugbe2llDT0sV3/HAGxhMRERRAUoSFveVNRyxff7CW09OjsAQYuXveCMICTTzz0T7OGRtPY7uVbcUNXqpYnEoS6sKndXcbnD02npc+K2BKWiRnZQ+PE6LHMyU1km2Hvwrp+tYu9pQ3MWukcw75iOAA7po/gk/yqggJNGI0KFblVXmrXHEKSagLn7ZybxUj40PZV9FMUW0bd83Nxjk56PA2LSOSw3Vt1LR0ArChwHnl6IwRsT3b3DYzi/iwQP60Kp+cjCg+lVAfFiTUhc9q6rDyRUEt546LZ/Hag2TFhnDehGPOEzesnJYeBcB2V2v907wqwgJNTO41902Q2ch9545iS1E96dHB5FU0U9YgMzn6Owl14bNW76vG5tCkRwezo6SR22dlYhwmwxYHMiklApNBsaWoni6bg49yK1g4IYGAPidCv5GTRmZMMJsPOW959+k+aa37Owl14bNW7qkkJsTMhoO1hAWauGra8Bm2OBBLgJHpmdF8vKeCdfk1NHXYuHhS0lHbBRgNPHDeGA7VthEdYmbVXgl1fyehLnyS1e7g031VzBwRw/LdFVyTk0qIn0/WdaIum5rMwepWbntlMzEhZmaPiu13u4snJTEhOZwum53P8mto6pCrS/2ZhLrwSZsL62jusKEBu9bc4povXXzl6mmpTE2LxKDgZxePI9Bk7Hc7g0Hx4AVjaem002Vz8HGuXIjkz6TpI3zSij2VBJoMbDhYy/zRcWTGhni7JJ9jNhl4864Z1LZ0kRhhOe62c0fFcmZWFJsP1fPOthKuHkZX4A430lIXPkdrzcq9lYyMD6W2tYtbZmZ6uySfFWA0DBjoAEopfnLhOBwa1uXXcrjPfOzCf0ioC5+zr7KZkvp2mtqtZMWGMHfU8JhS19NOS49i3mhnv/tLnxd4uRrhKRLqwud8uKsCBRTXt3PzjIxhM/viqfDzi8cD8J/NxbR32b1cjfAECXXhU7TWLNtRRmxoICFm47CaffFUGJUQxvzRcXTaHLy8vtDb5QgPkFAXPmV3aROFNa3UtXVx1bRUwiwB3i7J7/z6yoko4IXVB3HI/Uv9joS68Cn/3liEyaCwOzQ3z8jwdjl+KTUqmLmj42jusPH65sPeLke4mYS68BkNbV28u70Us9HArJExjEoI83ZJfuvpayajgGc/3u/tUoSbSagLn/GnVfl0Wh20We3cLBcbeVRCuIXpWVFUt3TJXZH8jIS68An7Kpp5ef0h4sICSYkMYsHYeG+X5PeeuGISAE8uz/NyJcKdJNSF12mteWTpboLNRqqaO/nWWRly27VTYFRCGBkxweRXtbC/omngHcSQIL85wuve21HGxsI6xiWFYTYZuG56mrdLGjbuPWcUAI++l+vlSoS7SKgLr2rptPGbD/cyPimcXSWNXDYlmegQs7fLGjYun5qMJcDAhoI6DtW0ersc4QYS6sKrfrdiP1XNnZyeEUW71cHts7K8XdKwEmA0cMlk5zzsz62UkTD+QEJdeM3u0kZeWV/I9dPTWLGnglkjYxifHO7tsoad7puPvLejTG535wck1IVXaK15fFku0SFmJqaEU9nUyR1zsr1d1rB0RmY0EUEBaA1//1ymDhjqJNSFV6zZX83mQ/Xct2AU//6imJHxocyT2Ri9wmQ0cOHERIwGxZKNRTS2yZ2RhrIBQ10p9bBSarNSqkkpVa2UWqaUmthnG6WUekwpVaaUaldKrVZKTfBc2WKo+9OqfFKjgsiICWZPeRN3zM6S2Ri96IKJidgcmjarg39tLPJ2OeJrGExLfT7wZ2AmcA5gA1YqpaJ7bfMg8ABwLzAdqAI+VkrJdd7iKLtLG9lSVM9ts7J4ZX0RMSFmrjgtxdtlDWszR8QSFmgiMdzCy+sK6bDKtLxD1YChrrU+X2v9stZ6t9Z6F3ATEAfMAmcrHbgfeEpr/bbWejdwCxAG3Oi50sVQ9cr6QwSbjUxNi2BVXhU3zcjAEtD//TXFqWE2GVgwLp6WThs1LV28/WWJt0sSJ+lk+tTDXPvVu77OAhKBFd0baK3bgbU4W/dHUUotUkptUUptqa6uPokSxFBV39rFezvKuGpaCq+sLyLEbORWuV2dTzh/QiItnTayY0N4aW0BdpmWd0g6mVD/PbAd2OD6OtH1se+sQJW91h1Ba71Ya52jtc6Ji5OTY8PJh7vL6bI5mDMqjvd3lnHTjEwig+ViI18wb0wcgSYD6THBHKptY0VuhbdLEifhhEJdKfU7YDZwtdZaOt3ECVu2o4wRcSGsyK0g0GTgjjlysZGvCDabmDs6jryKZtKjg/jL2gK0ltb6UDPoUFdKPQvcAJyjte5919ruP+cJfXZJ6LVOCKqaO9hYWMfcUXEs3V7G9dPTiQ0N9HZZopcLJyZS0djBwvEJbC9uYPOh+oF3Ej5lUKGulPo9XwV633k6C3GG98Je21uAOcB6N9Up/MBn+2vQGiqbO1AKvjNPLjbyNRdMTCTEbKS+1Up0iJnFawsG3kn4lMGMU38euA3nSJZ6pVSi6xEKoJ3/nz0HPKSUuso1hv0VoAVY4rHKxZCz9kA1UcEBfJxbyTWnp5IUEeTtkkQfwWYTl05J5n+5FVyXk8bKvZXkV7V4uyxxAgbTUr8H54iXT4DyXo8f9drmaeBZ4HlgC5AEnKe1bnZrtWLIcjg0nx+oITwoAKUU33NN+Sp8z7U5abR12YkLMxNoMvDXz6S1PpQMZpy6OsbjsV7baK31Y1rrJK21RWs9zzVeXQgA9pQ3UdvaxeG6Nm48M52USGml+6pp6ZGMSQjjjc0lXD0thXe+LKWqucPbZYlBkrlfxCmxZr/zeoRAo4Hvnj3Sy9WI41FKcff8EeyrbGZccgRWh4N/rJepA4YKCXVxSnTf3Pi22VnEhcmIF193yeQk0qODeWtLMQvHxfPPL4po7bR5uywxCBLqwuPsDs2OkkYCDIrvzJURL0OByWjgu2ePYEdJIxOSI2lst/KfLcXeLksMgoS68Li3thRjd2gWjk+Qq0eHkGtOT2NiSjivbTrMaWmR/O3zQmx2h7fLEgOQUBce5XBonvvkAAA/PG+Ml6sRJ8JoUDx+2UQqmjpICLdQUt/O8t1yPaGvk1AXHvXfbaWUN3YQGRTAiGtwSGwAABpGSURBVLgQb5cjTtDpGVFcNS2FT/IqSYsKYrFMHeDzJNSFx7R12fi/j/YRYFTMHhmLc5ZmMdT85MKxBJqMBJtN7CptZENBrbdLEschoS485qW1hVQ0dWC1a6ZnRQ+8g/BJ8WEW7j93FPsqmwm3mHhJpg7waRLqwiMqGjt4cc1BpqRFAM5/48XQdcvMTFKjggizBPDpvmr2V8rF4r5KQl14xDMr9mF3aDJjQggxGxmbKHc2HMoCjAZunZlJaUM7gSaDTPTlwyTUhdvtLm3k7S9LuHVWJvsqmpmWEYXJKD9qQ903pqdhCTCQGRPM0u2llDe2e7sk0Q/5TRNu5XBoHl+WS1SwmZtmZLCvspnpmdKf7g/CLQEsHJ9IRVMnWmteWH3Q2yWJfkioC7d6a2sJmw/V89AFY8ivakFryMmU/nR/ceVpyTS2W5k5IpbXNxVT1iCtdV8joS7cpralk98s38v0zCiuPT2NLYfqMBkUU9MivV2acJM5o+IIt5gIDTShkda6L5JQF27z5PI8WjpsPHHlJAwGxeZD9UxIiSDYbPJ2acJNAowGzhkbz4aCWq4+PZU3Nktr3ddIqAu3+KKglre2lnDn3GxGJ4TRYbWzo7iBHBnK6HfOm5BIfZuV2SNj0Wj+vDrf2yWJXiTUxdfWZXPw83d3kxoVxPdddzTacqieTpuDWSNjvFydcLe5o+MIMCp2lTZybU6atNZ9jIS6+Npe+qyA/KoWfnX5RILMRgDW7K/CbDRwVraEur8JDTSRkxHNmn3V3DN/BAB/+lRa675CQl18LUW1rfzhkwNcNCmRs8fG9yxfs7+aM7KipT/dT80fE0deRTMmg4Hrp6fzxuZiCmtavV2WQEJdfA1aax5ZmkuA0cAjl0zoWZ5f1cL+yhbmj4nzYnXCk+a5vrdr91dz74KRBJoMPLNin5erEiChLr6GD3dVsGZ/NQ+cN5rECEvP8n9vLMJoUFw2NdmL1QlPGpMQRmK4hdX7q4gPs3DHnGw+2FnOjuIGb5c27Emoi5PS1GHl8WW5TEwJ5+YZmT3Lq5o6WLLxMFeelkJ8mOXYTyCGNKUU80bH8dmBGmx2B3fOySImxMxTy/NkvnUvk1AXJ+V3K/ZT3dLJnXOyeXdbKavyKtlaVMeDb+/E5tDce85Ib5coPGzemDiaO2xsL24gzBLAveeMZENBLWsP1Hi7tGFNzmKJE7azpIFX1h9iQnI4972+/Yh1RoPil5dPICNG7nLk72aNjMVoUKzeV01OZjQ3npnB39YV8tvlecwZGYvBIDdF8QYJdXFC7A7NT/+7i3CLidyyJm44I5075mTR1G6ltqWLkfGhZMZKoA8HEUEBTEuPZM3+an50/hjMJgM/Om8M972+nf9uK+Xq01O9XeKwJN0v4oT8Y8Mhdpc2EWA0MCE5nF9fMZERcaGclh7FueMTJNCHmXmj49hV2khNSycAl05OZkpaJL/9Xx6tnTYvVzc8SaiLQStraOeZj/YxLimM2tYu7lswCqP8iz2szRvtvDZh7f5qAAwGxSOXjKequZMX18hkX94goS4GpXtMukODQpEdG8K54xK8XZbwsgnJ4cSGmlnjCnVw3rrw8qnJLF5bQEl9mxerG54GFepKqblKqfeUUqVKKa2UurXPeqWUekwpVaaUaldKrVZKTTjG04kh6KPcClbureSWmRnsKW/i+jPS5ESYwGBQzB8Tz6q8Kjqs9p7lD10wFqWcM3eKU2uwLfVQYDdwH9DfzD0PAg8A9wLTgSrgY6WU3JjSDzR1WHlkaS7jk8IJCnDO7XLJZLmwSDhdPjWZ5g4bn+ZV9SxLjgzirnkj+GBnOZsK67xY3fAzqFDXWn+otf6p1votwNF7nVJKAfcDT2mt39Za7wZuAcKAG91dsDj1nv5fHjUtnTx19SQ+2FXO9MwokiODvF2W8BEzR8QSFxbIO9tKj1j+nbkjSIqw8Mv3c3E45IKkU8UdfepZQCKwonuB1rodWAvMdMPzCy/aWlTHv744zK0zszCbDOyvbOGyKdJKF18xGhSXT0lm9b4q6lu7epYHmY385MKx7C5t4q2tJV6scHhxR6gnuj5W9lle2WvdEZRSi5RSW5RSW6qrq/vbRPiALpuDh9/ZRUpkEA+cN5plO8owGhQXTkrydmnCx1w5LQWrXR8V3pdNSWZaeiRPf7SP5g6rl6obXrwy+kVrvVhrnaO1zomLk5n8fNXitQfZX9nCLy+fQLDZyLId5cwcEUNsaKC3SxM+ZkJyBGdkRfPK+kPY7F/10CqlePTSCdS0dPL8pzLE8VRwR6hXuD72Hd+W0GudGGIKqlv4w6p8Lp6UxIJxCewoaeRwXRuXSteLOIY7ZmdR2tDO/3KP/LWfkhbJVdNS+PvnhRySOdc9zh2hXogzvBd2L1BKWYA5wHo3PL84xewOzU/e3kWgycCjl44H4L3tZZiNBs6f0G+PmhAsGJdAZkwwL31WeNRMjQ9dMBazycDP390tszh62GDHqYcqpaYqpaa69kl3fZ2und+h54CHlFJXKaUmAq8ALcASTxUuPOfldYVsOlTHo5dOID7cgt2heX9nGfPGxBERFODt8oSPMhoU356dxY7ihqOGMSaEW3jwgjF8nl/D0u1lXqpweBhsSz0H2OZ6BAGPuz7/pWv908CzwPPAFiAJOE9r3ezWaoXH5Ve18PRH+zh3XDxXT0sBYP3BGqqaO6XrRQzomtPTiA4xs3htwVHrvnlmBlPTIvnV+3toaOvqZ2/hDoMdp75aa636edzqWq+11o9prZO01hat9TzXeHUxhNjsDh54cwfBZiO/uWoSzksQ4O+fFxIbGsj5E2RaAHF8QWYjN8/I4JO8Kg5UHtmmMxoUT141iYZ2K09+KFeaeorM/SJ6/Hn1QXYUN/DrKyb23LVoT1kTn+6r5uYZGQSajF6uUAwFN8/IxBJg4KXPjm6tj0sK5445WbyxpZiNBbVeqM7/SagLwNnF8tzK/VwxNblnCgCHQ/PI0t1EBQdw84wML1cohoroEDPXnp7Gf7eVUtXUcdT6+xaMIjUqiJ+8s4v2Lns/zyC+Dgl1QXVzJ/e9vp3M2BCeuHISAJ02O48vy2VLUT0PXzSOyGCzl6sUQ8nts7Ow2jVv9nMlabDZxNPXTKawppUnl+/1QnX+TUJ9mLPaHdz3+jaa2q38+ZvTMBkVf/2sgLlPf8qrG4r49uwsrpU72IgTlBUbwplZ0fxnS3G/877MHBHL7bOy+MeGop652IV7SKgPY1prHn0vl/UHa/n1FRNJjQrmmy9t5Ncf7CU7NpR/fftMfnHJ+J4TpkKciOump1FU28bGY8zS+OAFYxgZH8qP39pBXauMhnEXCfVh7G+fF7Jk42Hunj+Ca3PS+PGbO9hW3MDvr5/Ka4vOYvaoWG+XKIawCycmEWYx8cbmw/2utwQYee66qdS3Wrnv9W3YZSZHt5BQH6aWbi/liQ/3ctGkRH583hj+t7uc5bsr+NF5Y7h8aoq3yxN+IMhs5PKpySzfXUFjW/+TeU1MieDxyyfw2QHniXrx9UmoD0Pv7SjjB29sZ3pmNP/v2qkoBc+tPMDohFDunJPl7fKEH7l+ejqdNgdLd5QeZ5s0vpGTyh9X5fPxnr6TvYoTJaE+zLy7rZT7X99GTmY0r9w2nSCzkY2FdeRVNHP7rCxMRvmREO4zMSWC8Unh/GdL8TG3UUrxy8snMjk1gu+/to0dxQ2nsEL/I7/Bw4TWmj+tOsD9b2x3TpF623SCzSYAXl1/iMjgAOl2ER5xbU4qu0ub2FvedMxtLAFG/nbLdGLDzNz+ymaKamU2x5Mloe6HOqx21uyvZvOhOrTWdNkc/PitnTyzwnlx0au3n9ET6KUN7azYU8l109MIMssVo8L9Lp+aQoBRDXj3o7iwQF657QzsWnPL3zdR0Xj0hUtiYCZvFyDca1VeJT99ZzcVriv5JiSHE2AwsL2kgfsWjOL+c0cdMUTxX18UobXmprPkilHhGdEhZhaMTeDdbaX85MKxBByni29EXCh/v3U6N/11Ize+9AWvLTqLhHDLKax26JOWuh95Y/Nh7nh1C5HBAbx823TunJ3FnrImdpQ08ORVE/nBwtFHBHqH1c7rmw6zcHwCqVHBXqxc+Ltrc1Kpbe1iVV7VgNtOS4/i1dvPoLKpgxsWfyEt9hMkoe4nXlpbwENv72LOqDjevmsmu0sa+du6QpIiLBgNipV7qo66su+97WXUt1m5ZWamd4oWw8a80XHEhQXyry+KBrV9TmY0r7iC/eoX1lNQ3eLhCv2HhPoQ53Bofv3+Hp74cC8XT0ri6Wsmc/eSL/l/H+/n0inJfPzDefzikvF8klfFC2u+ukekze7g+dX5jE8KZ0Z2jBffgRgOTEYDt87M5LMDNeSWNQ5qn+mZ0by+aAYdVjvXvLiBnSUyKmYwJNSHsOK6Nu74xxb++nkht8zI4OYZGVz+p3V84brs/7nrphISaOLmGRlcOiWZ/7diH+vyawB4Y0sxRbVtR3XJCOEp3zozgxCzsd8baBzLpNQI3rp7JsFmIzcs/oLPD9R4sEL/IKE+BHVY7fx+5QHO/d0aNhys5ZFLxpEcaeHGv27EbDLwzj0z+dZZGT1hrZTiqasmkR0XyveWfMkzH+3jV+/v4azsaM4dF+/ldyOGi4jgAG48M533d5aTV3Hs4Y19ZcWG8PbdM0mNCubWlzcdc9oB4aS8fRPYnJwcvWXLFq/WMJR8uq+KR5bupriunYsnJ3Hv2SN5ZsU+Vu6t4sKJifz2msmEW/q/j2hhTSvf+ecW9le2kJMRxYs3nU5saOApfgdiOKtv7WLB79aQHh3Mm3fNOO5ImL6aOqx8999f8tmBGhbNzebmGRn8/fNDrD9YQ0pkEN+akcHZY4ZPI0UptVVrnXPUcgn1oaGutYtfLsvl3e1ljIwP5fHLJhBsNvK9Jduoau7gZxeN45aZmQN2pWitaWizEhkcIN0uwiuW7Sjj3te2cdVpKTxx5SSCzEYcDs3B6hZ2ljTS1GElPszC2KQwsmNDjvg5tdkdPLI0lyWbDmNUCoMBzsyKoaC6hbLGDh69dDy3zRoeU10cK9RlnLqP01rz7vZSfvX+Xpo7rNy3YBR3z8/mtU3F/ObDvcSHWXjzrplMTYsc1PMppYgKkRteCO+5dEoyBdWtPLtyP5/kVZEcGURJXRvNnbajtk2LDuLsMfHMHRVHbFggO0sa2OC6DZ5da1Ijg3n00vFkxIRw72tf8viyPaRHB7Ng3PC9n6601H2U1e5gw8FaXlh9kA0FtZyWHslTV00mKdLCQ2/tZPnuCs4dF88z106RuxKJIWljQS3vfFlKdUsnyZEWpqRGMjUtkpjQQMob29l2uIHV+6pYl19Lu/Wr296NSwrnoQvGYDYa+P7r22nttPHElRO5aFISV/15PWWN7Sy/bw5JEUFefHeeJ90vQ8TS7aW8sPogh2pb6bA6iA01c9+CUXzzzAz2lDfx3SVfUlLfzk8uGMsdc7KkC0X4vQ6rnZ0ljTS2W8mKDWFE3FddMlVNHXzvtW1sKqzjiqnJ3DYrkxte2siklAiW3HkWRoP//n5I98sQYQkwkhwZxOyRseRkRjN/TByBJgNLNh3m8WV7iA4288ais8jJjPZ2qUKcEpYAI2dk9f/zHh9uYckdZ/KnT/P506p8Ps+v5ZrTU/nHhiJeXHOQ75498hRX633SUvdxNS2dPPzOLj7eU8nc0XE8+40pxMiIFSGOklvWyI/f3Mme8iaiQ8w0tHXx5l0zOD3DPxtAx2qpyzh1H/a/3RWc/+xa1uyv5ucXj+OVW6dLoAtxDBOSI3jve7N44sqJgMah4bq/fMEbmw/j7cbrqSQtdR9UXNfG48v2sHJvJROSw3n2uqmMTgjzdllCDBktnTae+3g/f1tXiNYQF2bmW2dmcMVpKWTEhBx33+rmTgprWokLCyQzJthnz1vJidIhoK61i8VrC3h5XSFGg+Lec0bx7dlZmE3yD5UQJ+NgdQu3/n0TxfXtPctGxYdy7vgE5o+O47T0qJ7frx3FDfx9XSEf7CzH5pr8bmJKOI9eOoHpPngOS0Ldhx2ubePfm4r414Yi2qx2LpuSzE8uHOv3Q7KEOBUcDs37u8p5cfVB9vS5+5LJoEiKsNBld1DZ1EmI2cj1Z6Qzf0wcBdWtLF5bQFljO4vmZvPDhaMJNH39G8lUNXXw3o4ydpU28svLJhIR3P8V4AM5JaGulLoH+DGQBOQC92utPzvePsM11CubOliVV8WHu8r5PL8GBVw0KYn7FoxilHS1COER+yqa+XBXOSv3VpJb1v/8M9mxIYxOCGN0YhgZ0UGsyqvig10VjE0M47nrpzI2MfykXruty8aLawpYvPYgHVYHKZFBvHr7dEbGn9zvu8dDXSl1HfAv4B7gc9fH24DxWutjzsDj76He1mWjrKGDsoZ28qta2F3ayM7SRvKrnPNDp0QGcW1OKtdNT5OWuRCnUGOblZrWTqKDzVgCjGwvbmBrUR07Sxo5UNVCUW0rvW9BoBSgITU6iOSIIExGhc2uae20YbVrwoNMZMaEMD45nInJEcSGBWI2GbDaHCzfXc7fPi+kpqWLhePjWTR3BFmxIUQFm096LP2pCPWNwE6t9Z29lh0A3tJaP3ys/U421PeUNdHYbiXAqDAZDZgMigCjAZNRYXZ9NBkMR60f6AA6HBqbQ9Nld9DWZaOt005bl935uetja6edNqudtk4bLZ02mjtsNHVYae6w0dzz0UZju5XGdusRzx8XFsjklAimZUSxYFw8YxLCfPZEjBDDWYfVTn5VCweqmjlc205BdQsbCmqpbunEXR0cG3+64KRv1+fRi4+UUmbgdOCZPqtWADPd8Rp9Pf1RHqv3VZ/wfgblnLA/wOAMe4fW2Owau0NjdThO+JulFIQGmgi3BBBmMRFmMZEQbmFkvHNZYoSFlMggkiODyIgJlvstCjFEWAKMTEyJYGJKxBHLbXYH9W1WAk0G7A7d03hr67LT2GElv7KZ/KoWiura6LQ6iAoJYFxiOJmxwYBCA91BE2Zx//Wf7nrGWMAIVPZZXgmc23djpdQiYBFAenr6Sb3gwxeOY9HcbGx2jc3hwGp3hrPV7sBqd2BzaGx21/Je620OB112h/NzuwODQWEyfNWaN7pa9AFGRZDZRIjZSLDZSLDZdOTHQNfnAUYMfnwpshDiSCajgbiwr64XOWqCvAmJp7iiI3llmgCt9WJgMTi7X07mOcYkhgFyQlEIIXpz1wDoGsAO9J3vMgGocNNrCCGEGIBbQl1r3QVsBRb2WbUQWO+O1xBCCDEwd3a//A74p1JqE7AOuAtIBl5042sIIYQ4DreFutb6DaVUDPBznBcf7QYu0loXues1hBBCHJ9bT5Rqrf8M/NmdzymEEGLwZKYoIYTwIxLqQgjhRyTUhRDCj3h96l2lVDVwMidTY3GOjxfuJcfV/eSYup8cU8jQWsf1Xej1UD9ZSqkt/U1mI74eOa7uJ8fU/eSYHpt0vwghhB+RUBdCCD8ylEN9sbcL8FNyXN1Pjqn7yTE9hiHbpy6EEOJoQ7mlLoQQog8JdSGE8CMS6kII4Ud8MtSVUoFKqT8qpWqUUq1KqfeUUqmD2O8epVShUqpDKbVVKTXnGNsppdRypZRWSl3j/nfgmzxxXJVS0a7nzFNKtSulipVSL7hm7PQ7g/0Z67X9PNd2HUqpAqXUXV/3Of2Ru4+rUuphpdRmpVSTUqpaKbVMKTXRs+/CR2itfe4BvACU4bzJxjRgNbAdMB5nn+sAK3AnMA74I9ACpPez7Y+ADwANXOPt9zuUjyswEXgHuAwYCcwDcoEV3n6/Hjh+g/4Zc22fBbS6thvn2s8KXH2yz+mPDw8d14+A21w/n5OA/+K8C1u0t9+vx4+ntwvo5xsWAXQB3+y1LA1wAOcfZ7+NwEt9lh0AnuyzbDpQDMQPp1D39HHts/4i1/OGe/t9u/kYntCxAH4LHOiz7K/Ahq9zfP3t4Ynj2s8+oThvuXmpt9+vpx++2P1yOhAArOheoLUuBvYCM/vbQSlldu23os+qFb33UUqFAUuARVrrKveW7fM8dlz7EQ50Am1fo16fcpLHYkY/238E5CilAr7G8fUbnjiux9gnDGd3c/1Jljpk+GKoJ+L8i9p3sp5K17r+xAJG1zbH2+dF4H9a6+VuqHOo8eRx7aGUigR+hbPlZTvpan3PCR8L1/L+tje5nu9kntPfeOK49uf3OLsaN5xcmUPHKQt1pdSvXScmj/eY78HXvwmYAvzYU6/hDd4+rn1qCQWWAaXAg6fiNYUYiFLqd8BsnH3udm/X42luvZ3dAJ4D/jXANoeBs3D+5Y4FqnutSwA+O8Z+NThboQl9lifgPDkCsAAYD7QopXpv84ZSaoPWevZAb8BHefu4Aj2B/qHry0u01h0DVj60DPpY9FJxjO1trudTJ/Gc/sYTx7WHUupZ4HrgbK11wdeudgg4ZS11rXWN1jpvgEcbsBXnmeyF3fu6ht2NA9Yf47m7XPst7LNqYa99fgZMBqb2eoBzJMzN7nmXp54PHNfucxX/w/lH4yKtdYs736MvGOyx6GPDMbbforW2nuRz+hVPHNfuBUqp3wM3AOdorfPcU/EQ4O0ztf09cA69KwHOBU4DPqXP0DsgD/her6+vwzm64w6cQfV7nMOiMo7zOsNm9IunjivOE1AbcA5jHIWzv7P7Yfb2e3bz8RvoWPwD+Eev7buH3j3n2v4O1/59hzSe0M+tvz08dFyfB5qAc/r8TIZ6+/16/Hh6u4BjfJMDcY5BrcU5gmIZkNZnGw081mfZPcAhnCMvtgJzB3id4Rbqbj+uwHzXPv095nv7PXvgGB7vWKwGVvfZfh7wpWv7QuCuE3nO4fJw93E9zs/kY55+L95+yCyNQgjhR3xxSKMQQoiTJKEuhBB+REJdCCH8iIS6EEL4EQl1IYTwIxLqQgjhRyTUhRDCj0ioCyGEH/n/8SGlZYNwWU4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "az.plot_dist(trace_ppc[\"mu\"], ax=ax)\n",
    "az.plot_dist(inf_ppc[\"mu\"], ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f7ba28e08d0>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxb9Znv8c8jyZZ3O4kdOwnZ9wQIIQkQyEYhLUtpoS2lCzC0Fxia0pZb2plpS++UOx2my5SW6ZSh0AFKuW0zhbK2QGjYCQnE2UhCEgjOnjh2Ysf7Ium5fxwZgrBj2ZZ8pKPn/Xr5lVjnHOl7XpYeH//ObxFVxRhjjHf53A5gjDEmuazQG2OMx1mhN8YYj7NCb4wxHmeF3hhjPC7gdoBYpaWlOm7cOLdjGGNMWqmsrKxV1bLutqVcoR83bhxr1651O4YxxqQVEdnd0zZrujHGGI+zQm+MMR4Xd6EXkWUiUiUibSJSKSIL4zxugYiERGRz/2MaY4zpr7gKvYhcAdwB3AbMBlYBT4nImF6OGwI8AKwcYE5jjDH9FO8V/TeB+1X1HlV9S1W/BhwEvtLLcf8N/BZ4bQAZjTHGDECvhV5EsoE5wIqYTSuAs09w3DKgHPjhQAIaY4wZmHiu6EsBP1Ad83g1UNHdASJyCvDPwJWqGu7tBUTkehFZKyJra2pq4ohkjDEmXgnvRy8iQWA58C1VrYrnGFW9G7gbYO7cuTZvsklZf3urmhWbD1HX2kFedoDL54xmweRSt2MZc0LxFPpaIIzTDHO8cuBQN/uPAKYD94nIfdHHfICISAi4SFVjm4GMSWn76lq4+t7Xebem+QOPP7bhAHPGlPCnG+bj81lvZZOaen1nqmoHUAksjdm0FKf3Taz9wCnAacd93QW8E/1/d8cYk7Je2H6YxT99gXdrmjl1VDG/+/IZvPwP5/LgtWcwblgelXvque6BSrdjGtOjeJtubgd+JyKvA68CNwAjcQo4IvIAgKperaqdwAf6zIvIYaBdVa0vvUkrVbVNXPvbtURU+bdPncLnz3i/R/HooXms/OZizv7xc6zcdpg7n3+bZedOdjGtMd2L629NVV0O3ATcAmwAFuA0wXTNrTAm+mWMZ3SGIlz2q1WEIsptl36wyHfx+308eeNCsgM+fvrMDlbtrHUhqTEnFnejoqreqarjVDWoqnNU9aXjti1R1SUnOPYHqnryALMaM6hueLCS+tZOLps9ks+f2fN1TFlRkPu/NA8FvnTfGzS2dQ5eSGPiYHePjOnG61VHWLntMGUFQX52+axe9z97YinLlkygPRThHx/eNAgJjYmfFXpjuvG/l28E4O6r58Tdm+bbH5tGQTDA05sP0dQWSmY8Y/rECr0xMZ7dcoj99a2cPqaE2WOGxH2ciHDD4glEFG59cksSExrTN1bojYnxw7++BcC/x9FkE2vZkokEAz4eW3+ASCSS6GjG9IsVemOOs/1gA7uPtDClvIAJZQV9Pt7n83HZ7FF0hCPc+cLOJCQ0pu+s0BtznB88sRWA7144vd/PccvFM/AJ/ObluGYAMSbprNAbE9XY1smaqiMMzc9mybTh/X6egpwACyaVUt/aycq3YucCNGbwWaE3JuonT28novDlc8YN+Lluudj5i+Dnz+4Y8HMZM1BW6I2JemLjAfw+4YZFEwf8XFMqiqgozmHLgQaOtXYkIJ0x/WeF3hhgw5466ls7mTN2CIFAYj4WV505FgV++dw7CXk+Y/rLCr0xwM//5jSxfOO8SQl7zmsXjscn8Oj6/Ql7TmP6wwq9yXiRSIRVO49QEPRzzqSyhD1vMMvPzJHF1DZ1sK+uJWHPa0xfWaE3Ge+RDQfoDCvnT49dW2fgrjrLmQzNmm+Mm6zQm4z3m5feBeDmj05N+HN/Zs5J+H3Cs1utm6VxjxV6k9FCoQjbqxspLwwyemhewp/f5/MxY0QRR5s72F/fmvDnNyYeVuhNRnts434iCh+Z3v8BUr353LzRANxlUyIYl1ihNxlt+Rv7ALhu4YSkvcYVc0fjE6z5xrjGCr3JaJv21VMQDPRrArN4BQI+JpYVcKihjbqW9qS9jjE9sUJvMtb6PXW0hSLMHRv/nPP9ddnsUQDc+8qupL+WMbGs0JuMde+ruwC4an7y17W/av5YAJ5682DSX8uYWFboTcZ69Z1aAj7h3KnJuxHbpTAni/KiIO/WNhMO24IkZnBZoTcZqaaxjaPNHUytKIx7TdiBWjK1jIjCk3ZVbwaZFXqTkf472lb+ydNGDtprfvmc8QA8tHbfoL2mMWCF3mSoZ7Y4V9VfPHPsoL3m1IoicrJ8bNhXP2ivaQxYoTcZKByOsPtICyOKc8gPBgb1taePKKKxLcTBYzZK1gweK/Qm4zy1+RARhYWTSwf9tS85dQQAD7y2e9Bf22QuK/Qm4/yp0mkjv3r+4DXbdLkiOh3C32yUrBlEVuhNxtmwt57sgI+TR5UM+mvnB7MoLcimqrYZVR301zeZyQq9yShHm9s51trJlOHJm/KgN/PGDSUUUd7YddS1DCazWKE3GeWPr+8FYOmMxC8yEq/L5zjNN39Ys9e1DCazWKE3GeWZLYcA+OJZg98+32XJ1FJ8Aq+9e8S1DCazWKE3GWV7dSNFOQFKC4KuZfD5fIwemkd1QxudIZsOwSSfFXqTMd6ubqStM8KpJxW7HYVFk0tR4IlNB9yOYjKAFXqTMR5c7fRd/+Rpo1xO8n7T0aPr97ucxGQCK/QmY7z0di2QGoV+WkURwYCPDXttOgSTfFboTUZQVfYcbWF4YZDsQGq87aeUF9DQFuJoU4fbUYzHpcY73pgke33XUcIRZc4grCYVr4/NrADg96/bdAgmuazQm4zQNTXwYE5L3JvPn+GsbGWLhptks0JvMsIbu44iwNLp7g2UijWsIEhBMMDbh5vcjmI8Lu5CLyLLRKRKRNpEpFJEFp5g38UiskpEjohIq4hsE5FvJSayMX23r66V0oIgfn9qXdtMKS+gpSNMTWOb21GMh8X1rheRK4A7gNuA2cAq4CkR6WlV5SbgP4BFwAzgh8CtIrJswImN6aPN++sJRTQl+s/HWhJdr/bP66ybpUmeeC9vvgncr6r3qOpbqvo14CDwle52VtVKVf2jqm5R1SpVfRB4BujxrwBjkuXhSqeIXnByhctJPuzyOScBsPIta6c3ydNroReRbGAOsCJm0wrg7HheRERmR/d9sa8BjRmoVTud/vMXRxf9SCUjSnLJyfKx7VCj21GMh8VzRV8K+IHYS45q4ISXSCKyT0TagbXAnap6Vw/7XS8ia0VkbU1NTRyRjInfriMtlORmkZc9uMsGxmtimdOfvrGt0+0oxqOSfWdqITAXuAG4SUSu6m4nVb1bVeeq6tyysrIkRzKZZO/RFtpDEaaNKHQ7So8WRJc0fHyDzXtjkiOeQl8LhIHYfmnlwKETHRhtn39TVe8Bbgd+0J+QxvTXoxuc9vmum56p6NOnO+30z1o7vUmSXgu9qnYAlcDSmE1LcXrf9OW13Jsb1mSkl3Y4TYGXptBAqVhTygvJ8gub9x9zO4rxqHgbLW8HficirwOv4jTFjATuAhCRBwBU9ero918DqoDt0eMXAd8C7kxYcmPisL26kZyAj4riXLejnNBJQ/LYVdtMOBxJub7+Jv3FVehVdbmIDANuAUYAm4GLVLVrko7Y/vR+4MfAOCAE7AT+iegvBmMGQ2tHiIbWENMqUrd9vsu8cUOoqm1m5bbDfHRm6nUDNekt7ksHVb1TVcepalBV56jqS8dtW6KqS477/heqOlNV81W1WFVPjx5vy+mYQfP0ZucW0pnjh7mcpHeXzHK6fv71zYMuJzFeZH8jGs96ZqtT6FNpIrOenDOxFBGo3F3ndhTjQVbojWdt3HsMn8DsMSVuR+mVz+djeGGQA8faUFW34xiPsUJvPElVqW5oo6IoBxFxO05cTh1VQjii1vvGJJwVeuNJ6/bUEVGYNTr1r+a7fOxkZ6jKn20dWZNgVuiNJ3WNMv3ozNSZf743F53s3JBdvfOIy0mM11ihN560+t2jAFw4M/UmMutJXjBAUU6AqiPNbkcxHmOF3njS7iPNFOdmkZPtdztKn0yrKKKtM8LB+la3oxgPsUJvPOdQQxttoQhTy1N/oFSsxVOdSf0esXZ6k0BW6I3nPLLOWQh80ZRSl5P03WWzRwHwwg6brtskjhV64zldE5l1Fc10MrIkl5yAj20HG9yOYjzECr3xnLcONRIM+Bg1JM/tKP0yrjSfhrYQLe0ht6MYj7BCbzylozNMfUsnY4amZ5EHOHP8UACe2XLC5R6MiZsVeuMpK7Y6i3ecES2W6eiSWc7cPF3nYsxAWaE3nvJ09Cr4E7NSfyKznswZOwSfwKZ99W5HMR5hhd54yoa99fjEmd89XYkIwwuDHDrWbhOcmYSwQm885eCxNoYX5uDzpfdb+5STSgirTXBmEiO9Pw3GHOfN/fWEI8opJxW7HWXAls5w5uh5LDpnjzEDYYXeeMaj0dGkS6enz0RmPbnoZGc5wdVVNsGZGTgr9MYzXtvpTGR28anpM5FZTwpyssgP+qmqsQnOzMBZoTeeUVXbTGFOgPxgXGvep7yJpQU0d4RpbOt0O4pJc1bojSfUNrXR2hlm8vACt6MkzFkTnUXN/2ILhpsBskJvPOHR9c5Ny3ScyKwnXYuar3zLBk6ZgbFCbzzhhe3ORGafmn2Sy0kSZ+bIYvw+4c39NsGZGRgr9MYTth5sINvvY8ywfLejJFR5YZDDDW02cMoMiBV6k/ZCoQhHmzsYPTTX7SgJd8pJxUQU3txnA6dM/1mhN2nvue2HAZg7Ln0nMuvJ+dExAY9vtBWnTP9ZoTdp76/RXikf90D/+VgXntI1cOqoy0lMOrNCb9Leuj31CHBOtDuilxQEsyiwgVNmgKzQm7R3oL6VssJg2k9k1pMJXQOnWm3glOkfb34yTMbYfqiBUESZObLI7ShJMz/6l8qTm2yCM9M/VuhNWnskOlDqfA9MZNaTroFTz0XHChjTV1boTVpbtbMWeH/5PS+a8d7AKetiafrHCr1JaztrmsgP+inKzXI7SlJ1DZwypj+s0Ju01dDaSXN7mIll3pnIrCddA6c27bV1ZE3fWaE3aeuJ6M3Jsz3YrTJW1z2IxzbaDVnTd1boTdpa+ZYzIvbS00a5nCT5LjrFGQy2eqetOGX6zgq9SVub9x8j4BOmjfBu18ou+cEAhcEAVUds4JTpOyv0Ji1FIhFqmtoZWZLjdpRBM2l4AS0dYY42t7sdxaQZK/QmLa2pOooqzB5d4naUQbNgsrOoStciK8bEK+5CLyLLRKRKRNpEpFJEFp5g30+JyAoRqRGRRhFZIyKfSExkY+Dx6E3JC0723kRmPfn06c6iKiu32YpTpm/iKvQicgVwB3AbMBtYBTwlImN6OGQx8BxwcXT/vwKPnOiXgzF98cauOsDbI2JjjSvNJzvgY+sBW3HK9E28V/TfBO5X1XtU9S1V/RpwEPhKdzur6jdU9Ueq+rqqvqOqtwKVwKWJiW0y3d6jLQzJyyIrkFmtj2OH5lHX0kl7Z9jtKCaN9PopEZFsYA6wImbTCuDsPrxWIVDXh/2N6db+uhbaQxGmVXi/t02sMyc4i6s8s+WQy0lMOonncqgU8AOxDYPVQEU8LyIiXwVOAn7Xw/brRWStiKytqbGJm8yJPbrBWW1p8dQyl5MMvktnO2MGntpshd7EL+l/94rIp4GfAl9Q1d3d7aOqd6vqXFWdW1aWeR9e0zcv7nAmMrvsNO9OZNaTuWOH4hdhg02FYPognkJfC4SB2Lte5cAJLytE5DM4V/FXq+oT/UpoTIzthxoJBnyUF3tvMfB4VBTnUN3QRiQScTuKSRO9FnpV7cC5kbo0ZtNSnN433RKRz+IU+WtU9aGBhDSmS2tHiGOtnYwvzXc7imtmjykhorCmym55mfjE23RzO3CNiFwrItNF5A5gJHAXgIg8ICIPdO0sIp8D/h/wT8BLIlIR/Rqa4Pwmwzy5yVkIfH4GTGTWk65F0B+P3qswpjeBeHZS1eUiMgy4BRgBbAYuOq7NPbY//Q3R5/5F9KvLi8CSgQQ2ma3rJuRn5452OYl7usYOrNl11OUkJl3EVegBVPVO4M4eti050ffGJMrGffVk+YXpGTCRWU8Cfh/D8rPZe7TF7SgmTWTWaBOT1jo6wxxp6mDM0Dy3o7hu5qgiOsPKO4cb3Y5i0oAVepM2nooOEjprQua2z3dZOsNpvnm40trpTe+s0Ju08dfojdjL557kchL3fTK6GPor79gAQ9M7K/QmbazfW4/fJ5w2eojbUVxXlJtNYU6Adw7bQiSmd1boTVoIhyPUNLYzekhmDpLqztSKQlo7wxw61uZ2FJPirNCbtPDs1moUOHO8DcXoct604QA8VLnX5SQm1VmhN2nhiU3OQiOfOt3a57tcPscZS/Dc9sMuJzGpzgq9SQvr9tTjEzjDrujfU1oYJC/bz45D1sXSnJgVepPyIpEIhxraGFmSi4i4HSelTB5eQFN7mPrmDrejmBRmhd6kvBd31KIKZ4yzq/lYS6Y67fQPr9/nchKTyqzQm5TXtdDIpbMzb/753nSNKXh2iy0Ybnpmhd6kvLW76/AJLJhU6naUlHPSkDxysnxsPWgLhpueWaE3KS0SiXCwvpWKohx8Pnu7dmdSWQENbSEaWjvdjmJSlH1yTEpb/e5RIgpzxtlo2J4snuIsv/nwOmunN92zQm9S2p+jxesTs0a5nCR1fXae059+hbXTmx5YoTcpbU3VUUTgvGm2aHxPxg7LJyfLx+YDx9yOYlKUFXqTslSV/fWtjLD2+V5NGV5IY1uIOutPb7phnx6Tsl59p9ba5+N03nSnP/0f39jjchKTiqzQm5T18Dqn//xls619vjefP9NZtvnpzdZObz7MCr1JWa9H2+cXT7b2+d4ML8yhMCfA9mrrT28+zAq9SUnhcIQDx1oZWZyL329v03icOqqYts6IrSNrPsQ+QSYlPbf9MKpw1gSb3yZel5zmTBHx4GprpzcfZIXepKRH1jvt85fPHe1ykvRx2WnOvYyXdtg6suaDrNCblLR2Vx1+n9iKUn0QzPIzvDDI7iMtqKrbcUwKsUJvUk5HKMzhxnbGDs2z+ef76IzxQwmr8so7dlVv3meF3qScJzcdBOAcm62yzz4bbepa/rrNe2PeZ4XepJzHNjjrw3b1DTfxWzi5FL9PWLPrqNtRTAqxQm9Szvo9dQQDPmaMKHI7StoREcYPy6OmsZ2W9pDbcUyKsEJvUsqB+lYa2kJMq7Ai318fmV4OwJ8qrfnGOKzQm5Ry36tVgE17MBBXzx8LwBMbD7icxKQKK/QmpTy71Zmr5XPzTnI5Sfo6aUge+UE/Ww7YdAjGYYXepIxQKMLuoy2MKsklJzvgdpy0NmtUCa2dYXYebnI7ikkBVuhNynhkw35UYclUm8RsoC493Wn6euC13S4nManACr1JGV03D//XgvEuJ0l/l84aiQAv7DjsdhSTAqzQm5SxaV89+UE/E8oK3I6S9rKz/FQU57D3aAvhcMTtOMZlVuhNSti8/xhtnRFOH2OrSSXK4illRBSeePOg21GMy6zQm5Rwb7Rb5efm2WyVidLVBPaHNTZtcaazQm9Switv1+ITuPDkCrejeMbk8kIKggE27K13O4pxmRV647pjrR0cbmxnYlkBPp+9JRNp/oShtIcivLaz1u0oxkVxf6pEZJmIVIlIm4hUisjCE+w7QkR+LyLbRCQsIvcnJK3xpPte2QXAJbNGuhvEg65bOBGAe16ucjmJcVNchV5ErgDuAG4DZgOrgKdEpKfpBYNALfAjYE0CchoPezw6VP+as8e6nMR7zpgwlJyAjzXvHnE7inFRvFf03wTuV9V7VPUtVf0acBD4Snc7q+ouVf26qt4P2HyppkdtnWGqapsZVZJLUW6223E8ae64oTR3hNlobfUZq9dCLyLZwBxgRcymFcDZyQhlMsfvXtuNAhedMsLtKJ517UKn982dL7zjchLjlniu6EsBP1Ad83g1kJAuEiJyvYisFZG1NTW2BFomeahyLwA3LJ7gchLvWjJ1ONkBH6veseabTJUSXRxU9W5Vnauqc8vKbJ6TTBEKR3j7cBPlRUGGFQTdjuNpp48uobE9xNYDx9yOYlwQT6GvBcJAeczj5cChhCcyGeOPb+wlorB0RuxbyyTalxd0Nd/sdDmJcUOvhV5VO4BKYGnMpqU4vW+M6ZcHVzszK3713EkuJ/G+pTPKyfILL+2wptFMFG/Tze3ANSJyrYhMF5E7gJHAXQAi8oCIPHD8ASJymoicBhQBQ6Pfz0hkeJO+mttDbD/USEVxDiOKc92O43kiwmmjS2hoC/GWNd9knLgKvaouB24CbgE2AAuAi1S1a7LrMdGv462Pfi0ELon+/68JyGw84L9eeAcFLp9jK0kNlusXOje8f7HybZeTmMEW9zI+qnoncGcP25Z085j0P5bxuj9V7kOAZUsmuh0lYyydWUFOwMeL1nyTcVKi143JLLtqm6luaGfGiCJybcnAQbV4ahltnRGe3mxTF2cSK/Rm0P1sxXYA/n6xXc0Ptm9/bBoA//m8DZ7KJFbozaBbue0wwYCPS2bZaNjBNml4AeVFQbbsb6C5PeR2HDNIrNCbQfX8tmpaOsIsmFSKiN3GccMVc0ejwJ12VZ8xrNCbQfXL55zi8q2PTnU5Seb6ypKJiMD/RBdjN95nhd4MmvbOMBv21jM0P5vpI4vcjpOxcrMDnDKqmJrGdrYdanA7jhkEVujNoPnV8+8QUfiM9Z133U3nTwbgx09tczmJGQxW6M2g+d3q3fgEvnHeZLejZLyPTCunIBjg5bdrCYUibscxSWaF3gyKl9+uoa6lk3njhpIftL7zqeCy2SMJRZTfvGLLDHqdFXozKH4UbSL4/sdtuqNU8a2PTkWA+1ZZofc6K/Qm6Wob29lyoIGRxTmcPKrY7Tgmqjgvm1mjS6huaOeNXbbip5dZoTdJ971H3wRsOuJU9H+if2H94PEtLicxyWSF3iRVU1uIZ7dWU5gT4ItnjXU7jolx+tghjB6ay5YDDeyubXY7jkkSK/QmqW59cgsRhesW2Zqwqep7F00H4LuPvOlyEpMsVuhN0nSEIjy6fj85AR832nTEKeuCk0cwLD+bVTuPUNfc4XYckwRW6E3S/Otf3qIzrFx51lh8PnurpbKvnzcZBf7pz5vcjmKSwD59Jik6Q2H+8Poesv0+/vGCaW7HMb24ev5YinICrNhaTW1ju9txTIJZoTdJceuTW+kIR/jCmWPICtjbLNWJCP9wwTRU4Zv/s9HtOCbB7BNoEq6pvZM/vL6XYMDHLRdPdzuOidOVZ41lWH42L71dQ1WN9cDxEiv0JuFuXr6RcET56rmTCPjtLZZO/vWykwH4+wfXupzEJJJ9Ck1C7ahu4Jmt1QzJy+JrH7EBUunmgpNHMH1EITuqm3hknc1X7xVW6E1CXfvbSgD+/fJZtoJUmvrN381DBL736GY6bWZLT7BCbxLm1y/uZM/RFk4fU8J508vdjmP6aVRJLp+fN4aWjjA3/8luzHqBFXqTELWN7fz0me0EfMJ915zhdhwzQD+8dCbFuQEe33iATfvq3Y5jBsgKvUmIq+9dQyiifO/i6RTnZbkdxwyQz+fjV188HYAv3fcGkYg14aQzK/RmwO5/tYqtBxuZWlHIl84Z73YckyALJpVxwcxyjjR38PU/bnA7jhkAK/RmQLYfauD/PrmVgE944MvWZOM1v/rC6QzJy+LJTQd5fON+t+OYfrJCb/qtvTPM5b9+jYjCL644jfKiHLcjmQTz+30s//v5+AT+9/KN7KhudDuS6Qcr9Kbfrrh7NQ2tIa6YN5qPzxrpdhyTJFPKC/nJZ04lHFEu+9Wr1LfYDJfpxgq96Zev/2EdG/bWM2l4AT/+9KluxzFJ9pk5o7l+0XiaO8JceMfLhMJ2czadWKE3ffbDv2zl8Y0HGZafzRM3nuN2HDNIvnvRDJZMKePgsTY+++vX3I5j+sAKvemTf3lyK795uYr8bD/P3LSI3OyA25HMILrvS/OYVJbPuj31fPGe1aiq25FMHKzQm7j940Mb+e9XqigIBlh582JKC4NuRzKDTET4y9cXMnpILq/uPMInf/WqNeOkASv0pleNbZ189OcvsnztPopzAzx382IqinPdjmVcEszy89zNS5haXsimfceY/28rOVDf6nYscwJW6M0Jrdtdx5m3rWRHdROnnlTM6u+cz3DrRpnxsgI+nr5pIR+bWU5NUweLfvI8j26wfvapygq96VYkEuE7f97Ep/5rFS0dYa5fNJ7Hb1xAbrbf7WgmRYgIv75qLv98yQwiqtz0xw1cfe8aOjrDbkczMSTVbqbMnTtX1661RQ/c9JdNB7nl0Tepa+kkP+jnritPZ+Hk4W7HMimsqraZz9+zmkPH2sjL9vOdC6dx1fxxbsfKKCJSqapzu91mhd50ebhyLz95ZjvVDc7i0B+bWc4vP3862bbmq4mDqvL9xzbz+zV7iCiU5GVx3cIJ3LBoAn5baSzprNCbHtU0tvGzFTt4YtMBmtudP7nnjRvC7ZefxuhheS6nM+motqmdm/9nIy+9XYMqZPmF86YN56bzpzBtRJHb8TwrIYVeRJYB3wZGAFuAm1T15RPsvxi4HZgJHAB+oqp39fY6VuiTKxKJ8Pz2Gh6q3Mfa3XXUNDpX7wGf8JFpw7n1EzMZUWI9aszANbZ18uOnt/Hndftp6XAuIgqCAeZPHMYXzxzDosml+Hx2pZ8oAy70InIF8CCwDHgl+u+XgBmquqeb/ccDm4F7gTuBBdF/P6eqD5/otazQJ05rR4jV7x5h1c6jbN5fT1VtC4cb24hEf+R+nzCxLJ8rzxrLlWeOsQ+dSZpnthzivld3sX5PHe3R5QkFGFaQzfjSfGaMKGLe+KGcM6mUIXnZ7oZNU4ko9GuATap63XGPvQ08pKrf6Wb/HwOfUtXJxz32G2Cmqs4/0WtZoX9fKBShqSNEY3uIprZOmjvCtLSHaOkI09weorapg5rGdmqb2qlr6aCuuZOjLR0ca+mkpTNMOPLBn/dI4fkAAAfiSURBVK3fJ5QXBTlrwjA+O3c0Z44fauu6mkG3emcty9fuY92eOg7Ut9IZjnmfipAX9FOck0VRbhZD87MpLQgyvDCbkrxsinIDFEW3DcnLYkheNsW5WeQHA/h9krHv6RMV+l7Hr4tINjAH+PeYTSuAs3s4bH50+/GeAf5ORLJUtbO31+2rbYcauPH362lpD1HT1P7hHfQD/3Srx20a+218zV0f2iuOwwZ6x8TvE/Ky/IwozqG8KIcp5QXMGzeURZPLbCSrSQlnTSzlrIml731/rLWDVe8c4fVdR9l6oIED9a3UtXRwsKGNfQkYiNVV9o+v/4K8v6GH/Xt8viT8HgkGfFQU57JkShm3fHxGwp8/nolKSgE/UB3zeDVwfg/HVAB/62b/QPT5Dh6/QUSuB64HGDNmTByRPiwn4GdqeSH1rR0fKpbv/6A//BP64A///f/IBx95bz/50L4f3DN2P+c/sc8Wu59w/Psu4BOyAz6yAz6CAX/0Xx85AT85WT5ysv3kZfkZkpdNRUkOo0pyGT0kj5K8rIy9mjHpqzg3mwtPGcGFp4zodvuxlg52H2lhz9EW6ls7aWjtjP6VG6KpPURze4jWzjCdoQhhVcIRJaxKJKKEIxA57jFVp3eQ6gcvqhTt9WLw+MYPjdl5oBdoudlO/UrW/bGUmJFKVe8G7gan6aY/zzGuNP+9NS6NMd5RnJfNqXnZnDq6xO0oaSueu2+1QBgoj3m8HDjUwzGHetg/FH0+Y4wxg6TXQq+qHUAlsDRm01JgVQ+HvdbD/muT0T5vjDGmZ/H2p7sduEZErhWR6SJyBzASuAtARB4QkQeO2/8uYJSI/CK6/7XANXz4hq4xxpgki6uNXlWXi8gw4BacAVObgYtUdXd0lzEx+1eJyEXAz4Gv4AyY+npvfeiNMcYkXtw3Y1X1TpxBT91tW9LNYy8CdnfUGGNcZkMhjTHG46zQG2OMx1mhN8YYj0u5aYpFpAbY3ctupXizP74Xz8uL5wTePC8vnhN487y6O6exqlrW3c4pV+jjISJre5q8J5158by8eE7gzfPy4jmBN8+rr+dkTTfGGONxVuiNMcbj0rXQ3+12gCTx4nl58ZzAm+flxXMCb55Xn84pLdvojTHGxC9dr+iNMcbEyQq9McZ4nBV6Y4zxOE8UehE5Q0SeFZEmEWkUkVUiUtr7kalNHE+JiIrIZ9zOMxAiMlREfiki20SkVUT2ish/RWdFTRsiskxEqkSkTUQqRWSh25kGQkS+IyJviEiDiNSIyBMicrLbuRIpeo4qIv/pdpaBEpERIvLb6M+qTUS2isji3o5L+0IvImfiLET+AnAW7y9k7oUFTm4GIm6HSJCRwCjgH4BTgCuBRcAf3AzVFyJyBXAHcBswG2fhnadEpH8LHaeGJTiz0p4NfARnFbi/ichQN0MlioichbMe9Sa3swyUiJQAr+IsL30xMB34GnC412PTvdeNiKwCnlfV77mdJZFEZB7wZ5xfXNXA5ar6kLupEiu6ZsGTQImqNridpzcisgbYpKrXHffY28BDqvod95IljogUAMeAS1X1CbfzDISIFAPrgGuBfwY2q+qN7qbqPxG5DVisquf09di0vqIXkeHAfOCgiLwiIodF5GUROc/tbAMhIoXA74HrVbXX39ZprAhoB1rcDtIbEcnG+aW7ImbTCpyrYa8oxKkLdW4HSYC7cX4JP+92kAS5FFgjIsujtW6DiNwoItLbgWld6IEJ0X9vBe4FPga8DDwjIrNcSzVwdwFPq+pTbgdJluifof8C3KOqIbfzxKEU8OP8dXW8aqBi8OMkzR3ABpx1n9OWiFwHTMJZFc8rJgDLgHdxat0dwI+Ar/Z2YEoWehH5YfTmyYm+lvB+/l+r6r2qul5Vvwu8Adzg2gl0I95zEpGrgFnAt93OHI8+/KyOP6YAeALYj9Nmb1KAiNwOLAA+rapht/P0l4hMxbmP8gVV9cK9ui4+YJ2qfida6+4D/oM4Cn3cSwkOsl8AD/ayzx6gPPr/rTHbthKzjm0KiPecrgFmAE0xf5EtF5HXVHVBcuL1W7znBbxX5P8a/fbjqtqWrGAJVguEef8916UcODT4cRJLRH4OfA44V1XfdTvPAM3H+Qtsy3GfIT+wSERuAPJVtd2tcANwkA/XureAb/R2YEoWelWtJY75o0VkF87C41NjNk0B3kx8sv7rwzl9D6fX0PHeBL4FPJaEaAMS73nBe/censLpNXCBqjYlM1siqWqHiFQCS4E/HbdpKZDWi96LyB3AFThFfpvbeRLgUWBtzGP3AW/jXOl3DHqixHiV7mtdb+t3pGahj5eqqoj8FLhVRDYB64HP4nSzTMu766q6H6dJ4z3Rq5K96XylFS3yK3BuwF4K5ItIfnTzUVVNhw/f7cDvROR1nA/dDTjdRu9yNdUAiMivgKtwfiZ1ItJ1v6EpnX4RH09V64H64x8TkWac99lmd1IlxM+BVdGLweU4XXy/Dny3twPTutADqOovRCQI/AwYBmwBLlTVje4mMzHm4PwCBtgRs+1cnHEQKU1Vl0cHeN0CjAA2Axepaq9XVClsWfTflTGP3wr8YHCjmBNR1TdE5FKcv0q+j9Mk+n2ccRAnlPb96I0xxpxYSva6McYYkzhW6I0xxuOs0BtjjMdZoTfGGI+zQm+MMR5nhd4YYzzOCr0xxnicFXpjjPG4/w/M0khCWG+D1wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "az.plot_dist(trace_ppc[\"y\"], ax=ax)\n",
    "az.plot_dist(inf_ppc[\"y\"], ax=ax)"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}