{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Survival Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:27.970806Z",
     "iopub.status.busy": "2021-04-16T19:37:27.970215Z",
     "iopub.status.idle": "2021-04-16T19:37:27.972416Z",
     "shell.execute_reply": "2021-04-16T19:37:27.971942Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# If we're running on Colab, install empiricaldist\n",
    "# https://pypi.org/project/empiricaldist/\n",
    "\n",
    "import sys\n",
    "IN_COLAB = 'google.colab' in sys.modules\n",
    "\n",
    "if IN_COLAB:\n",
    "    !pip install empiricaldist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:27.976032Z",
     "iopub.status.busy": "2021-04-16T19:37:27.975593Z",
     "iopub.status.idle": "2021-04-16T19:37:27.977681Z",
     "shell.execute_reply": "2021-04-16T19:37:27.977307Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Get utils.py\n",
    "\n",
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://github.com/AllenDowney/ThinkBayes2/raw/master/soln/utils.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:27.980924Z",
     "iopub.status.busy": "2021-04-16T19:37:27.980273Z",
     "iopub.status.idle": "2021-04-16T19:37:28.658584Z",
     "shell.execute_reply": "2021-04-16T19:37:28.659033Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from utils import set_pyplot_params\n",
    "set_pyplot_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The following cell downloads the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.662755Z",
     "iopub.status.busy": "2021-04-16T19:37:28.662317Z",
     "iopub.status.idle": "2021-04-16T19:37:28.664128Z",
     "shell.execute_reply": "2021-04-16T19:37:28.664481Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "download('https://github.com/AllenDowney/ThinkBayes2/raw/master/examples/Osogbo2021.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll use Pandas to load the data into a `DataFrame`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.669139Z",
     "iopub.status.busy": "2021-04-16T19:37:28.668327Z",
     "iopub.status.idle": "2021-04-16T19:37:28.679256Z",
     "shell.execute_reply": "2021-04-16T19:37:28.678840Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TIME(In Days)</th>\n",
       "      <th>AGE OF PATIENTS</th>\n",
       "      <th>CENSOR</th>\n",
       "      <th>AGE AT MENARACHE</th>\n",
       "      <th>BREASTFEED</th>\n",
       "      <th>CONTRACEPT</th>\n",
       "      <th>DETECTION</th>\n",
       "      <th>NEOADJUVANT</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>48</td>\n",
       "      <td>44</td>\n",
       "      <td>1</td>\n",
       "      <td>18</td>\n",
       "      <td>1.5</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>898</td>\n",
       "      <td>43</td>\n",
       "      <td>0</td>\n",
       "      <td>17</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>36</td>\n",
       "      <td>40</td>\n",
       "      <td>0</td>\n",
       "      <td>15</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>23</td>\n",
       "      <td>56</td>\n",
       "      <td>1</td>\n",
       "      <td>15</td>\n",
       "      <td>1.5</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>300</td>\n",
       "      <td>53</td>\n",
       "      <td>0</td>\n",
       "      <td>12</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   TIME(In Days)  AGE OF PATIENTS  CENSOR  AGE AT MENARACHE  BREASTFEED  \\\n",
       "0             48               44       1                18         1.5   \n",
       "1            898               43       0                17         1.0   \n",
       "2             36               40       0                15         2.0   \n",
       "3             23               56       1                15         1.5   \n",
       "4            300               53       0                12         1.0   \n",
       "\n",
       "   CONTRACEPT  DETECTION  NEOADJUVANT  \n",
       "0           2          1            2  \n",
       "1           1          1            2  \n",
       "2           2          1            2  \n",
       "3           1          1            1  \n",
       "4           2          1            2  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.read_csv('Osogbo2021.csv')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "duration = df['TIME(In Days)']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "log_duration = np.log10(duration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.691538Z",
     "iopub.status.busy": "2021-04-16T19:37:28.691003Z",
     "iopub.status.idle": "2021-04-16T19:37:28.894479Z",
     "shell.execute_reply": "2021-04-16T19:37:28.894857Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvsElEQVR4nO3dd3xU95nv8c8jIVElqmgC0YsxBgwC4k7c7cQtceK4pDjF8Xqdspt1kk3blL3Zvevdzd0Ux9dxvE7ia+M4sTeOg0tccDfNFBswWFRJNCGBJCQQKs/94xwNIzGSBtBoZqTv+/XSS6fNmefMgfnqd8rvmLsjIiKSajKSXYCIiEgsCigREUlJCigREUlJCigREUlJCigREUlJCigREUlJCijpVGZ2r5l9t5PWVWBmh8wsMxxfamaf74x1h+t72sw+3Vnr667M7Dwz29QJ61lkZiVd/b6SvnoluwBJH2a2HRgBNACNwAbgt8B97t4E4O63n8C6Pu/uz7e1jLvvBAacWtWR9/s+MNndb4la/xWdse7uzt1fBaYl+n3MzIEp7l7Ule8rqUstKDlRV7l7DjAO+FfgG8CvO/tNzKzb/PHU3AJMVd3ps5buRQElJ8XdK939SeAG4NNmNhPAzB40s38Oh4eZ2VNmdtDMKszsVTPLMLPfAQXAn8NDeF83s/Fm5mb2OTPbCbwYNS36C3SSmS03s0oz+5OZDQnf67jDR2a23cwuNrPLgW8BN4TvtzacHzlkGNb1HTPbYWb7zOy3ZjYwnNdcx6fNbKeZ7Tezb7f12YSfwS/NbImZ1QAfNLMPmdlqM6sys+KwRde8/G/M7GvhcH74XneE45PDz85ivM9kM3s5/Cz2m9mjrertFbVs9LZ+xsxeN7OfmFkF8KNwH82MWj7PzA6b2fDoz9bMvmlmf2hVx3+Z2U/D4VvNbKOZVZvZVjP7YlufU6t1vBIOrg330Q2t92m4P+8ys3VmVmNmvzazERYcqq02s+fNbHDU8h8wszfCbVtrZoui5n0mrK/azLaZ2c3x1CldSwElp8TdlwMlwHkxZn8tnJdHcGjwW8FL/JPAToLW2AB3/7eo11wAnAZc1sZbfgr4LDCa4FDjT+Oo8Rngx8Cj4fvNjrHYZ8KfDwITCQ4t/rzVMucSHHK6CPiemZ3WztveBPwvIAd4DagJax8EfAj4GzO7Nlz2ZWBROHwBsDX8DXA+8KrH7pPsR8BzwGBgDPCzduppbWH4PsOBHwKPAzdGzf848LK772v1ukeAK80sFyKtw48DD4fz9wEfBnKBW4GfmNncjopx9/PDwdnhPnq0jUU/ClwCTAWuAp4m+Hc1jOD77MthXfnAX4B/BoYA/wD8MQze/gT/bq4IjwacDazpqEbpegoo6Qy7CL4EWqsHRgHj3L3e3dv6oo32fXevcffDbcz/nbu/6+41wHeBj1vnHEK7GfhPd9/q7oeAfwQ+0ar19gN3P+zua4G1QKyga/Ynd3/d3Zvc/Yi7L3X3d8LxdQRf9M0h9DJwnpllEATSvwHnhPMuCOfHUk9wqHV0+B6vncD27nL3n7l7Q/hZP0zLgLqJY6ET4e47gLeBa8NJFwK17v5WOP8v7r7FAy8TBGisP15O1s/cfa+7lwKvAsvcfbW71wFPAGeGy90CLHH3JeFn/ldgJXBlOL8JmGlmfd19t7uv78QapZMooKQz5AMVMabfDRQBz4WHU74Zx7qKT2D+DiCL4K/nUzU6XF/0unsRtPya7YkarqX9CzhabIeZLTSzl8yszMwqgdsJ63b3LcAhYA7Bl/lTwC4zm0b7AfV1wIDlZrbezD7b7ha2Ux/wItA3rHNcWMsTbbw2OsxaBJmZXWFmb4WHJQ8SBEJn7J9me6OGD8cYb94n44CPhYf3Doa1nAuMCv+4uYFgH+w2s7+Y2fROrFE6iQJKTomZzScIqOP+enf3anf/mrtPJDgc8/dmdlHz7DZW2VELa2zUcAFBK2I/wSG0flF1ZRIcWox3vbsIvtSi191Ayy/AE9H6/R4GngTGuvtA4F6CcGn2MnA9kB22Dl4mOCQ4mDYOP7n7Hnf/gruPBr4I3GNmkwk+C4j6PICR7dUXXoX5e4LguQl4yt2r29i2x4BFZjYGuC7cNsysN/BH4N+BEe4+CFjSaju7SjFBa3tQ1E9/d/9XAHd/1t0vIWjhvwf8Kgk1SgcUUHJSzCzXzD4MLAYecvd3Yizz4fBEvgFVBJemN4az9xKc6zlRt5jZDDPrR3Du5A/u3ghsBvqEFyNkAd8Beke9bi8wPjyMFssjwN+Z2QQzG8Cxc1YNJ1FjLDlAhbsfMbMFBCEQ7WXgTqD5YoGlwJeA18LtO46ZfSwMCYADBKHT6O5lQCnBZ5UZtqwmxVHjwwQti5uJcXivWbj+pcB/A9vcfWM4K5vgMy8DGszsCuDSON632cn+m4jlIeAqM7ss/Az6hBddjAkvrLg6PBdVR9B6jfkZS3IpoORE/dnMqgn+Qv028J8EJ8NjmQI8T/AF8CZwj7svDef9C/Cd8PDLP5zA+/8OeJDgcFsfwpPi7l4J3AHcT/DlXENwgUazx8Lf5Wb2doz1PhCu+xVgG3CEICA6yx3AD8PP7nsErZVoLxOEWHNAvUbQAnqFts0HlpnZIYLW2VfcfVs47wvAXUA5cDrwRkcFuvsygs9tNMHFB+15GLiYqCALW1xfJti2AwQh/GRH7xvl+8Bvwn8THz+B1x3H3YuBawguoCgj+Pd6F8F3XgbBBTy7CA5NX0CwfyTFmB5YKCIiqUgtKBERSUkKKBERSUkKKBERSUkKKBERSUlp10nksGHDfPz48ckuQ0REOsmqVav2u3te6+lpF1Djx49n5cqVyS5DREQ6iZntiDVdh/hERCQlKaBERCQlKaBERCQlKaBERCQlKaBERCQlJSygzOwBCx6d/W4b883MfmpmReEjnDt86qaIiPQciWxBPQhc3s78Kwh6u54C3Ab8MoG1iIhIJ3J3dpcfYvX7ezpe+CQl7D4od3/FzMa3s8g1wG/DR4C/ZWaDzGyUu+9OVE0iInJy6hsa2brrIO/tLOe9Hft5r7icqpo6emf14nffvprMzM5v7yTzRt18Wj52uiScdlxAmdltBK0sCgoKuqQ4EZGerKqmjvd2lrNpZzkbd+5nS+lBGhqPf65jXX0DO/ZWMnH04E6vIZkBFesx0DEfTuXu9wH3ARQWFuoBViIiCVBSVsXS1TtYtmEXu8qrO1y+X58spo0dSlNTYr6WkxlQJcDYqPExBE+4FBGRLlJz+Civv1vCi29v5/2SinaXHTlkANMKhjK9YCjTCoZSMDwXs1htjc6RzIB6ErjTzBYDC4FKnX8SEekapfur+f1LG1i2YRf1DccfusvMzGDiqEGcNm4Y0wqGMm3sUAbn9OnSGhMWUGb2CLAIGGZmJcA/AVkA7n4vsAS4EigCaoFbE1WLiIgcU1lTx7fue4lDh4+2mJ6RkcH8aaNYdOY45kweQXZWZpIqDCTyKr4bO5jvwN8m6v1FRCS2JW8VtQin8SMHceHc8Zw7aywD+/dOYmUtpd3jNkRE5OTtqTjEM8u2RMbvuHYeF82bkMSK2qaAEhHpIYr3VfGDB1+NtJ7yBvVj0ZxxSa6qbQooEZEeoHhfFd/99ctU19YBkNUrkzuuLUzIDbadRQElItIDPPj02kg49c7qxbc+eQ4zJxz3lPWUooASEenmivdVsaZoLwCG8U+fOY9pBUOTXFXHUrdtJyIineK5FVsjw/Onj0qLcAIFlIhIt+buvPFuSWT88oWTkljNiVFAiYh0Y+u27OPgoSMA5PbvzRkThye5ovgpoEREuqnq2jru+Z9VkfGFp+WTkZG4vvM6my6SEBFJM+7OocNHqaypo6qmjsqaOioP1VFVG/xunr6n4hDlVYeBoOfxj5w/LcmVnxgFlIhIkrk7tUfqWwZOzfGB0/y7qqaOJj+xR1x86SPzGT64f4K2IDEUUCIiXaT2SD1LlhVRWlbdMnhq62hsbErY+16/6DQWnDY6YetPFAWUiEgXeeiv7/Ls8i0dLxiHvr2zGNi/N7n9e0d+D2o1PrB/b4bk9iU3hTqAPREKKBGRLuDuvLm+pM35vbN6HR84A44PnObfWb2S+yiMrqCAEhHpAltKD1BVE3Q1lNu/N3deV9gicHpn6+u4NX0iIiJdYOXmPZHhuVNHMm/aqCRWkx50H5SISII1NDaxbENpZHzuVIVTPNSCEhFJoIqqw9y9+C127q0EIMOM2ZPSpzeHZFJAiYh0MnenpKya9dvKeGzpxkhXQwAfPnsKA/pmJ7G69KGAEhE5Re7Ojj2VrN++nw3by9iwY3/kgohmhvHJy87g6nOmJKnK9KOAEhE5CXsP1LBsQykbtu9nw/b91Bw52uayOf1687UbFqZVR62pQAElInKCVm7azd2PvEVDY2Obywzom82MccOYMSGP82cXMDBNb5ZNJgWUiMgJ2LrrAP/x6PHhlNu/N6ePz2PG+GGcPj6PghG5mKVPz+GpSAElIhKnA9VH+PFDb3C0PginvEH9+Mj505kxfhj5w3IUSJ1MASUiEqfFL6znQPWxx1d8+5PnMnZ4bpKr6r50o66ISBx27KnkhVXbI+NfvX6BwinB1IISEWmHu/P8ym3899PrcIJnMM2aNFxdFXUBBZSISBv2V9ZyzxOrWLtlb2Rar8xMPnXZrCRW1XMooEREYmhqcr5z/1LKDtZGpo0emsOXPlrIhFGDkldYD6KAEhGJYX9lbSScDOOqc6Zw40Wnk53V/Z/DlCoUUCIiMTQ2eWR4+OB+fPpyHdbrarqKT0QkhiY/FlCZGfqqTIaEfupmdrmZbTKzIjP7Zoz5A83sz2a21szWm9mtiaxHRCReDY1NkeGMDN2AmwwJCygzywR+AVwBzABuNLMZrRb7W2CDu88GFgH/YWbqh15Eki66N3IFVHIksgW1AChy963ufhRYDFzTahkHcizoH2QAUAE0JLAmEZG4PLNsS2R40ujBSayk50pkQOUDxVHjJeG0aD8HTgN2Ae8AX3H3JkREkmhPxSGWbdgVGb/qbD3DKRkSGVCx2sTeavwyYA0wGpgD/NzMjus7xMxuM7OVZrayrKyss+sUEWlh/bayFr1GjBs5MMkV9UyJDKgSYGzU+BiCllK0W4HHPVAEbAOmt16Ru9/n7oXuXpiXl5ewgkVEAN4vORAZnjlBDxlMlkQG1ApgiplNCC98+ATwZKtldgIXAZjZCGAasDWBNYmIdGj7noOR4Un5Ov+ULAm7UdfdG8zsTuBZIBN4wN3Xm9nt4fx7gR8BD5rZOwSHBL/h7vsTVZOISEfcnT0VNZHx/GE5SaymZ0toTxLuvgRY0mravVHDu4BLE1mDiEi89h2o4eePr6S6NrjEPDMzg6G5fZNcVc+lro5EpMdzd15YtZ0Hlqylrv7YnS4Lpo/WPVBJpIASkR7N3XlgyVqWvFUUmWYY150/jY9/8LQkViYKKBHpsRoam/jFEyt5Ze3OyLTmR2pMHTs0iZUJKKBEpAdrHU5nnT6GL390vh6pkSIUUCLSI+3YU9kinC4pnMhtV52pc04pRAElIj2Ou/P4q+9FxudPH80Xrz6ToFtQSRUKKBHpUQ7X1XPP/7zNG+8e6yr0uvOmKZxSkAJKRHqM4n1V/NvDb7KrvDoybd7UUUwr0AURqUgBJSI9Qt3RBn744KtUVB+OTLt0/kRuvWJ2EquS9iigRKRHeG7ltkg4ZWdlcvvVc7lgzrgkVyXtSegj30VEUkFRSQWPv3LsoohPXzZL4ZQG1IISkW6rqcl54tVNLH5xA01NwbNQB+f05cK545NbmMRFASUi3dLR+kZ+/NDrvLN1X2Ran+xefOkjhboRN00ooESkW/rT65tbhNPUsUP5yvXzGTlkQBKrkhOhgBKRbqehsYlnlm2JjF919hQ+eekZZGbqtHs60d4SkW5n+cZdHDx0BAjOOd2icEpL2mMi0u1s3XUgMrxoTgG9FE5pSXtNRLqd4n1VkeFxIwYmsRI5FQooEelWNu0sZ9WmPZHxAgVU2lJAiUi3UVlTx88fX4njAMyaNJyCEblJrkpOlq7iE5FuYXNxOf+++C3Kq4LujHpn9eKOawvVS3kaU0CJSNp7dvkWfr1kLY2NTZFpt111JnmD+iWxKjlVCigRSWsvvb2d+/68OjLev082X7l+PvOmjUpiVdIZFFAikra2lB7g3iePhdOEUYO468azGDG4fxKrks6igBKRtNTQ2MS/L36LhsZGAMbm5fLPn19En2x9rXUXuopPRNLSltID7DtYA0Df3ll84+azFU7djAJKRNLSjr2VkeG5U0cyaqg6ge1uFFAikpa27T4YGda9Tt2TAkpE0tJ7O8sjw1PyhySxEkkUBZSIpJ2aw0cp3hv0t5dhxrSCoUmuSBJBASUiaWfjzvJId0YTRg/SxRHdlAJKRNLOxu37I8OnFQxLYiWSSAooEUk7G3ZEBdR4BVR3ldCAMrPLzWyTmRWZ2TfbWGaRma0xs/Vm9nIi6xGR9HfkaANFpcceSDhjnAKqu0rYgVszywR+AVwClAArzOxJd98Qtcwg4B7gcnffaWbDE1WPiHQPO/ZU0tQUdAo7Ji+X3P69k1yRJEoiW1ALgCJ33+ruR4HFwDWtlrkJeNzddwK4+74E1iMi3UBF9eHIsG7O7d4SGVD5QHHUeEk4LdpUYLCZLTWzVWb2qVgrMrPbzGylma0sKytLULkikg4qD9VFhgeq9dStJTKgYj0lzFuN9wLmAR8CLgO+a2ZTj3uR+33uXujuhXl5eZ1fqYikjQOHjkSGB+X0SWIlkmiJvHmgBBgbNT4G2BVjmf3uXgPUmNkrwGxgcwLrEpE0tr+yNjI8aIACqjtLZAtqBTDFzCaYWTbwCeDJVsv8CTjPzHqZWT9gIbAxgTWJSBor3V/Nq+uOnTnQOajuLWEtKHdvMLM7gWeBTOABd19vZreH8+91941m9gywDmgC7nf3dxNVk4ikp/qGRpZt2MVjSzdGHus+OX8Isybqwt/uLKH9g7j7EmBJq2n3thq/G7g7kXWISHraU3GI51du44W3t1NVc+ziiAwzbr9mLhkZsU51S3ehDqxEJOXsqTjE/U+tYc37eyN97jXLMOOzH5rDhFGDklOcdBkFlIiknAf+spbV7+9pMW1obl8uLpzAxfMmMCS3b5Iqk66kgBKRlLOrvDoyPGvScK78wGTmThlJZqa6D+1JFFAiklIaG5uojDrf9NWPLdQNuT1Uu3+OmNmDUcOfTng1ItLjPbN8K7VH6gHo1yeLnL7ZSa5IkqWj9vLsqOGvJLIQEZE9FYdY/OL6yPi1507TlXo9WEcB1bprIhGRhNhdfojv/frlSOtp5JABXH3OlCRXJcnU0TmoMWb2U4J+9ZqHI9z9ywmrTER6jOZwau6pPKtXJn973TyyemUmuTJJpo4C6q6o4ZWJLEREeqaaw0f58e9ebxFO/3jz2cwYr46he7p2A8rdf9NVhYhIz9PU5PyfPyyPXFae1SuTb91yDrMmqQsjiaOzWDP7tJm9bWY14c/Ktp7bJCISr4qqw/zot6/y9uZjN+TeeV2hwkki2m1BhUH0VeDvgbcJzkXNBe42M9z9twmvUES6nZWbdvPzx1dSXXvsfqfrzpvGubPGtvMq6Wk6Ogd1B3Cdu2+PmvaimX2U4BHuCigROSHPLNvCr55aHRk3jOvOn8ZNF5+exKokFXUUULmtwgkAd99uZrmJKUlEuqvaI/U89NdjT9QZnNOXr1w/nzP02AyJoaOAOnyS80REWqg5fJRHXtjA4bpj9zn9y20fJFfdGEkbOgqo08xsXYzpBkxMQD0i0o24O1tKD/Dsiq28uq6Y+obGyLxrzp2qcJJ2dRRQs4ERQHGr6eOAXQmpSETSnrvz8pqdPPXm+2zbffC4+ZPyB7NozriuL0zSSkcB9RPgW+6+I3qimeWF865KVGEikr6eeHUT/y/qXFOzghEDuWz+RBadOY7sLPUSIe3rKKDGu/txh/jcfaWZjU9MSSKSzipr6vjD0vci41m9MjnnjDFcWjiRqWOHYKbOXyU+HQVUn3bm6ZGWInKcR55fT119AwBj83L50ecvIKefzjXJieuoJ4kVZvaF1hPN7HPAqsSUJCLp6rV3ivnryq2R8VsuO0PhJCetoxbUV4EnzOxmjgVSIZANXJfAukQkzWzbfZB7njj2d+uC6aOZN3VkEiuSdNdRZ7F7gbPN7IPAzHDyX9z9xYRXJiJp4+3Ne/iPR9+KHNobOWQAd36kUOeb5JR01IICwN1fAl5KcC0ikoaeWbaF+59ag4fPN+2T3Yuv33gW/fWodjlFcQWUiEgsa4r2tuhXb9jAfnzrlnMYN3JgEquS7kIBJSInxd15+Plj9zpNyh/MP958DoNz2rv4VyR+HT4PSkQkllWb97Cl9AAQ3Ov0zZvOVjhJp1JAicgJc3cWv7A+Mn7p/AkMydWtkdK5FFAicsJWvLc70sdeVq9MrjtvenILkm5JASUiJ8TdefTFDZHxyxdM1KE9SQgFlIickLc2lLJ9z0EgaD1de9605BYk3Zau4hORuBytb+TxVzfxxCubItOuXDiJQQPUepLESGgLyswuN7NNZlZkZt9sZ7n5ZtZoZtcnsh4ROXGNjU2s3LSbr/7sOR57aQMNjcFDB/v3yeYatZ4kgRLWgjKzTOAXwCVACUHHs0+6+4YYy/1v4NlE1SIi8XF39h2s5f2SCopKKigqPcCWXQc4Wt/YYrlJ+YO545p5DNQTcSWBEnmIbwFQ5O5bAcxsMXANsKHVcl8C/gjMT2AtIhJDVU0dRaUHgkAqreD9kgNU19a1uXz/PtncculMLp43gYwM9bMniZXIgMqn5aPiS4CF0QuYWT5Br+gX0k5AmdltwG0ABQUFnV6oSE/Q1OTs2FvJpp3lvLeznM3F5ew9UBPXa4fm9qVw+mhuuHCGWk3SZRIZULH+vPJW4/8H+Ia7N7bX67G73wfcB1BYWNh6HSISQ83ho2wqrmBTcTmbdpazubgi0tt4e/r1yWJy/mCm5A9h8pghTM4frJtwJSkSGVAlwNio8THArlbLFAKLw3AaBlxpZg3u/j8JrEuk23F3Ssqqeb+kgk07g0AqKauO9DDelszMDCaOGsSUMIgmjxnC6KED9JgMSQmJDKgVwBQzmwCUAp8AbopewN0nNA+b2YPAUwonkY5V19bxfskBNhWX835xBe+XVlB7pL7D1w3O6cv0gqFMKxjKtLFDmDBqEFm9MrugYpETl7CAcvcGM7uT4Oq8TOABd19vZreH8+9N1HuLdDe7yw+xbsteNpdUsHlnBbvKqzt8TYYZ40cNCgJpbBBKwwb2VetI0kZCb9R19yXAklbTYgaTu38mkbWIpKs1RXv559+81uHhutz+vZk6ZghTxg5hesEwJucPpk+27sWX9KV/vSIp7pllW44Lp4yMDCaMGsi0sUMjoTRicH+1jqRbUUCJpKj9lbW8ub6UNUV7I9NuvPh0Zo7PY+LowWRn6dyRdG8KKJEUsu9ADW9tKOWNd0t4v6Sixbwxeblcf8FpSapMpOspoERSwMtrdrDkrS0UlVbEnJ+ZmcHHFimcpGdRQIkk2ebicn76xxXHTc8w44yJwznr9HwWzMhXDw7S4yigRJLsryu3RYYzMjKYPSkMpdNGk9NPoSQ9lwJKJInqjjbwxrslkfEfffZ8po8blsSKRFKHnqgrkkTLNu7iyNGgf7zRQ3OYVjA0yRWJpA4FlEgSvbR6e2R40ZnjdB+TSBQFlEiS7DtQwztbygAwjAvm6FEyItEUUCJJ8sLb2yM9RMyePJxhA/sluSKR1KKAEkmCxsYmXly1PTJ+ceGEthcW6aEUUCJJsLpoLxXVh4Ggk9f500cnuSKR1KOAEkmC56PufbrwzPH0ytR/RZHW9L9CpItVVB1m1abdkfEL541PXjEiKUwBJdLF/rpyG00eXBwxY3we+cNyklyRSGpSQIl0ofXbyvjjK+9Fxi/RxREibVJXRyJdoLq2jrc37+GBJWtpbGwCYNzIQZw9c0ySKxNJXQookQTZU3GIFe/tZvnGXby3Y3/ksB7AwP59+Mebz9bFESLtUECJdKJd+6tZunoHy9/bTfG+ypjL9M7qxTduOou8QboxV6Q9CiiRTlJ2sJZ/uOcF6uobYs6fMmYI86eP5vzZBQonkTgooEQ6yWNLN7YIp16ZmcyeNJz500dROH00g3P6JLE6kfSjgBI5Ae5O2cFadpcfivnT7ItXz+X82QX0ydZ/MZGTpf89Iq24O+VVh2MG0J6KGhoaG9t9/cwJw7l0/sQuqlak+1JAiQBH6xv5w9KNrNi0m93lh6hvaD+E2jImL5fPf3hO5xYn0kMpoKTHc3d++adVvLJ2Z1zL5/bvzaihA6J+csgfOoCRQwfokJ5IJ9L/Junxnl2+9bhwGtA3u0UIjR6aw6ihAxg5pD/9+2YnqVKRnkUBJT3axh37eeDptZHxRWeO4zOXzyKnX+8kViUioICSHmzvgRr+98NvRroeGj9yEF+8ai7ZWZlJrkxEQJ3FSg91uK6ef3noDapr64DgvNI3bjpL4SSSQhRQ0uM0NTk/+f3ySFdEmZkZfOPGsxg+uH+SKxORaDrEJz2Gu7Ol9ABPvvE+qzYfe2DgHdfMY/q4YUmsTERiSWhAmdnlwH8BmcD97v6vrebfDHwjHD0E/I27r0WkE+0uP8Qra3fy2rpidpVXt5h37bnTWHTmuCRVJiLtSVhAmVkm8AvgEqAEWGFmT7r7hqjFtgEXuPsBM7sCuA9YmKiapOdoaGzi+ZXbeGn1DopKK2Iuc9bpY7j5kpldXJmIxCuRLagFQJG7bwUws8XANUAkoNz9jajl3wL09DY5ZXVHG/i3R95kTdHe4+b1zurFB07P5/zZBcyeNBwzS0KFIhKPRAZUPlAcNV5C+62jzwFPJ7Ae6QEO19Xz44feYMP2ssi0zMwM5k4ZyXmzCyicOpLe6u1BJC0k8n9qrD9NPcY0zOyDBAF1bhvzbwNuAygoKOis+qSbqaqp43/97vUWh/SuO28a15w7VTfeiqShRAZUCTA2anwMsKv1QmY2C7gfuMLdy2OtyN3vIzg/RWFhYcyQk55t34EafvSb11pcBPGpy2ZxzblTk1iViJyKRAbUCmCKmU0ASoFPADdFL2BmBcDjwCfdfXMCa5Fu6Gh9I5uKy3lnyz5eeHs7Bw8dAcAwPv/hOVy+cFKSKxSRU5GwgHL3BjO7E3iW4DLzB9x9vZndHs6/F/geMBS4JzxZ3eDuhYmqSdJbU5OzddcB1m3dxztb97FxR/lxj8XIzMzgq9cv4OyZut5GJN2Ze3odMSssLPSVK1cmuwzpAkfrG9m+5yDvlxxg/bYy3tm2j9oj9W0uP6BvNnfdeBYzJ+R1YZUicqrMbFWsxokuZ5KU0NTkFO+roqi0gvdLDlBUWsGOvVU0NTW1+7rRQ3M4Y9JwzpiYx+xJI+jXJ6uLKhaRRFNASVKUVx1m4/b9FJVWUFR6gC27DnC0vuOn2A7O6csZE/OYNXE4Z0wazrCB/bqgWhFJBgWUdKnGxiYefXEDj7+yCY9910ELo4fmMCl/MFPHDuGMicMZk5ejm2tFeggFlHSZsoO1/OSxZWzaGfNuAobk9GVy/mAmjxnC5PzBTMofzAA9vVakx1JASZd4a0Mpv3hiZYuLHCblD2bOpBGRQBqS2zeJFYpIqlFASUIdrW/kN8+s45nlWyLTMsy44aIZfOS86WRk6HCdiMSmgJKEKSmr4j8eXcbOvZWRacMG9uPvPrZAz18SkQ4poKTTuTtLV+/gvqdWt7gyb+Fp+dxx3TydVxKRuCigpFMdrqvnvj+v5pW1OyPTemVmcusVs7hswURdgScicVNASacoKaviuRVbeWn1jhYXQowemsPXbljI+FGDkleciKQlBZSctIbGJpZt3MVzy7fy7rZ9x82/cO54PvehOfTR85dE5CTom0NOymvrinnwmXUcqD583LyRQwZw40Wnc+6ssTFeKSISHwWUnJCaw0e576nVvLauuMV0w5g/fRSXLZykR6mLSKdQQEncNmwv46d/XEHZwdrItEED+nBx4QQuKZygfvFEpFMpoKRDDY1N/D5G/3kfPHM8n/vQbPr2Vg/iItL5FFDSpqP1jezYW8mvnlrNltIDken9+2Rz+zVz9VBAEUkoBZRQXVtHaVk1pfurKSmrjgzvrag5rsfxmROG86WPFupwnogknAKqh3B3yg7WBgG0PwihkrIqSvdXU1VT1+HrMzMzuPnimVx9zhRdACEiXUIB1c0crW9kd/khisuqWrSKdu2vpr6h4wcCRjOMEUP6M37UIK6/YDoTdLOtiHQhBVSaqqqpC8JnX9AKCoarKTtYG9eDAKNl9cokf1gO+Xk5jMnLIT8vlzHDchg1dADZWZkJ2gIRkfYpoNJAdW0d67aW8e7WfezcW0VJWRWHDh894fXk9u/NmLzcSBjlDwsCKW9QPx22E5GUo4BKQfUNjby3s5y1RXtZu2Uf23YdjLtV1HxYrkWLKBzO6dc7wZWLiHQeBVQKcHd27q1iTdFe1m3Zy/rt+zs8XxTzsFxeDqOG6LCciHQPCqgkKa86zLqwhbRuyz4qa460uWyGGZPHDGHWpOFMGztUh+VEpEdQQHWRI0cbeHdbGeu27GVt0T5KyqraXX7U0AHMmjSC2ZOGM3NCHv31kD8R6WEUUCfB3amrb6Sqpi74qT1KdW34u9V4ML+O6pqj7Z5HGtA3mzMmDmfO5BHMmjSc4YP7d+EWiYikHgUUQV9zVTV1VNcepSoqaCqjplXX1lFVEwRPZc1RGhpP7J6i1jIzMzitYBizJw9n9qQRTBw9SIfsRESi9KiAWr+tjBfe3t4ibKpqj3K4rr7jF3eCghEDIy2kGeOG0VsP8hMRaVOP+obcX1nLy2t2dMq6emVmMrB/Njn9epMb/h7Yvzc5/bLJ7debnP69yW0e7pdNbv/e9MrM6JT3FhHpCXpUQLV1H5Bh5PQPwiQ3DJacfuHvVuO5/YNlemdl6pCciEgC9aiAGj9yIHdcO++48BnQN1thIyKSYnpUQA3J7ctF8yYkuwwREYmDToqIiEhKSmhAmdnlZrbJzIrM7Jsx5puZ/TScv87M5iayHhERSR8JCygzywR+AVwBzABuNLMZrRa7ApgS/twG/DJR9YiISHpJZAtqAVDk7lvd/SiwGLim1TLXAL/1wFvAIDMblcCaREQkTSQyoPKB4qjxknDaiS6Dmd1mZivNbGVZWVmnFyoiIqknkQEV67rt1p3RxbMM7n6fuxe6e2FeXl6nFCciIqktkQFVAoyNGh8D7DqJZUREpAcy9/ie1HrCKzbrBWwGLgJKgRXATe6+PmqZDwF3AlcCC4GfuvuCDtZbBpxqf0XDgP2nuI5k0zakBm1Dakj3bUj3+uHUtmGcux93eCxhN+q6e4OZ3Qk8C2QCD7j7ejO7PZx/L7CEIJyKgFrg1jjWe8rH+MxspbsXnup6kknbkBq0Dakh3bch3euHxGxDQnuScPclBCEUPe3eqGEH/jaRNYiISHpSTxIiIpKSempA3ZfsAjqBtiE1aBtSQ7pvQ7rXDwnYhoRdJCEiInIqemoLSkREUpwCSkREUlK3Dqju0Jt6HNuwyMwqzWxN+PO9ZNTZFjN7wMz2mdm7bcxPh33Q0Tak+j4Ya2YvmdlGM1tvZl+JsUxK74c4tyHV90MfM1tuZmvDbfhBjGVSfT/Esw2dtx/cvVv+ENx7tQWYCGQDa4EZrZa5EniaoMulDwDLkl33SWzDIuCpZNfazjacD8wF3m1jfkrvgzi3IdX3wShgbjicQ3ADfbr9X4hnG1J9PxgwIBzOApYBH0iz/RDPNnTafujOLaju0Jt6PNuQ0tz9FaCinUVSfR/Esw0pzd13u/vb4XA1sJHjO2VO6f0Q5zaktPCzPRSOZoU/ra9SS/X9EM82dJruHFCd1pt6EsVb31lhk/tpMzu9a0rrNKm+D+KVFvvAzMYDZxL85RstbfZDO9sAKb4fzCzTzNYA+4C/unva7Yc4tgE6aT9054DqtN7Ukyie+t4m6MdqNvAz4H8SXVQnS/V9EI+02AdmNgD4I/BVd69qPTvGS1JuP3SwDSm/H9y90d3nEHSMvcDMZrZaJOX3Qxzb0Gn7oTsHVHfoTb3D+ty9qrnJ7UHXUllmNqzrSjxlqb4POpQO+8DMsgi+2P+fuz8eY5GU3w8dbUM67Idm7n4QWApc3mpWyu+HZm1tQ2fuh+4cUCuAKWY2wcyygU8AT7Za5kngU+GVMx8AKt19d1cX2o4Ot8HMRpqZhcMLCPZpeZdXevJSfR90KNX3QVjbr4GN7v6fbSyW0vshnm1Ig/2QZ2aDwuG+wMXAe60WS/X90OE2dOZ+SGhnscnkCepNvSvFuQ3XA39jZg3AYeATHl5KkwrM7BGCq3qGmVkJ8E8EJ1bTYh9AXNuQ0vsAOAf4JPBOeO4A4FtAAaTNfohnG1J9P4wCfmNmmQRf2r9396fS6TuJ+Lah0/aDujoSEZGU1J0P8YmISBpTQImISEpSQImISEpSQImISEpSQImISEpSQElKMrNvh70lrwt7RF7YSeu92mL0Ch/na79vZv8QY/q1ZjYjavyHZnbxqdTZxvubmb1oZrnh+KGOXtPOuu60oMdsj76JMnyPE+pN28w+Y2Y/P9laYqzvDDN7sLPWJ+mr294HJenLzM4CPkzQe3Vd+AWafQKv7+XuDbHmufuTHH/D9qm6FngK2BC+R6Ie83AlsDZGFz8n43WCmpe2mn4FMCX8WQj8MvzdZdz9HTMbY2YF7r6zK99bUotaUJKKRgH73b0OwN33u/suADPb3vwXv5kVmtnScPj7ZnafmT0H/NbMlllUJ5VmttTM5jX/tW9mA8N1ZYTz+5lZsZllmdkXzGyFBZ1d/tHM+rVVqJmdDVwN3B229CaZ2YNmdn1UvT82szfNbKWZzTWzZ81sS/PNjeFyd4Xvuc5iPGMndDPwpxg1mJndbWbvmtk7ZnZDOD3DzO4JW6JPmdmS5rrcfbW7b4/xHnH1pm1mt5rZZjN7meAm2ubpV4Wf/Woze97MRoR1vG9meVF1FZnZMDP7WFj3WjN7Jeot/kzQc4r0YAooSUXPAWPDL8B7zOyCOF83D7jG3W8ieDTJxwHCL9jR7r6qeUF3ryR4vlbzuq8CnnX3euBxd58fdna5EfhcW2/o7m8QtMjucvc57r4lxmLF7n4W8CrwIMGd9h8AfhjWdylBi2UBMAeYZ2bnx1jPOcCqGNM/Er5uNkHXM3eH2/wRYDxwBvB54Ky2tiNKh71ph+v+QVjPJcCMqNmvETwf6EyCffB1d28CHiIIWMIa17r7fuB7wGXhZ3111HpWAufFUa90YwooSTlhR5PzgNuAMuBRM/tMHC990t0Ph8O/Bz4WDn8ceCzG8o8CN4TDnwjHAWaa2atm9g7Bl+qpPrah+ZDiOwQPoKt29zLgiAX9ml0a/qwm6Al6OkFgtTYkfBZSa+cCj4S9TO8FXgbmh9Mfc/cmd98DvBRHrfH0pr0QWOruZR48p+zRqHljgGfDz+4ujn12DwCfCoc/C/x3OPw68KCZfYGgO69m+4DRcdQr3ZgCSlJS+GW71N3/CbgT+Gg4q4Fj/277tHpZTdTrS4FyM5tFEEKLY7zNk8AVZjaEIBBfDKc/CNzp7mcQtBRav8+Jqgt/N0UNN4/3IgiFfwlbYHPcfbK7/zrGehqaD0m2EitU2pvennh7026rj7SfAT8PP7svEn527l4M7DWzCwkC7ulw+u3Ad8L3XGNmQ8P19CHox016MAWUpBwzm2Zm0S2IOcCOcHg7QZjAsdBqy2Lg68BAd3+n9cywpbYc+C+CR1Q3hrNygN0WPN7h5tavi6E6fM3Jehb4rAXPOsLM8s1seIzlNgETY0x/BbjBggfJ5RE8on45weG2j4bnfEYQdHjbkXh6014GLDKzoeFn9LGoeQOB0nD4061edz/Bob7fN3/WZjbJ3ZeFF5bs51g4TgXejaNe6cYUUJKKBhD0mLzBzNYRnOP4fjjvB8B/mdmrQGMbr2/2B4JDd79vZ5lHgVtoeZjquwRfwn/l+MchxLIYuCu8MGBSHMu34O7PAQ8Db4aHxv5A7MD7C7FD5glgHcE5tRcJzvvsIXh2UgnBF/3/JdimSgAz+7IFPbOPAdaZ2f3hupYAWwl60/4VcEeMencT7I83gecJDks2+z7wWLh/9rd66ZME+/a/o6bdHV7Y8S5B0K4Np38w3F7pwdSbuUiaCC9O+K27X3ICrxng7ofCQ2fLgXPC8OpyZlYI/MTd2734wcx6E5xHO7et2wWkZ9B9UCJpwt13m9mvzCz3BO6Feiq8ECMb+FESw+mbwN8Q3yHTAuCbCidRC0pERFKSzkGJiEhKUkCJiEhKUkCJiEhKUkCJiEhKUkCJiEhK+v9HiCi4ynXb0wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from empiricaldist import Cdf\n",
    "from utils import decorate\n",
    "\n",
    "cdf = Cdf.from_seq(log_duration)\n",
    "cdf.plot()\n",
    "    \n",
    "decorate(xlabel='Survival time (log10 days)', \n",
    "         ylabel='CDF',\n",
    "         title='Distribution raw survival times')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.064186566935665, 0.849952354519006)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu = log_duration.mean()\n",
    "sigma = log_duration.std()\n",
    "mu, sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.898796Z",
     "iopub.status.busy": "2021-04-16T19:37:28.898324Z",
     "iopub.status.idle": "2021-04-16T19:37:28.900180Z",
     "shell.execute_reply": "2021-04-16T19:37:28.900554Z"
    }
   },
   "outputs": [],
   "source": [
    "from empiricaldist import Pmf\n",
    "\n",
    "def make_uniform(qs, name=None, **options):\n",
    "    \"\"\"Make a Pmf that represents a uniform distribution.\"\"\"\n",
    "    pmf = Pmf(1.0, qs, **options)\n",
    "    pmf.normalize()\n",
    "    if name:\n",
    "        pmf.index.name = name\n",
    "    return pmf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.905031Z",
     "iopub.status.busy": "2021-04-16T19:37:28.904601Z",
     "iopub.status.idle": "2021-04-16T19:37:28.906705Z",
     "shell.execute_reply": "2021-04-16T19:37:28.907185Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "qs = np.linspace(2, 4, num=101)\n",
    "prior_mu = make_uniform(qs, name='mean')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I chose the lower and upper bounds by trial and error.\n",
    "I'll explain how when we look at the posterior distribution.\n",
    "\n",
    "Here's the prior distribution for `sigma`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.912722Z",
     "iopub.status.busy": "2021-04-16T19:37:28.912090Z",
     "iopub.status.idle": "2021-04-16T19:37:28.914846Z",
     "shell.execute_reply": "2021-04-16T19:37:28.914286Z"
    }
   },
   "outputs": [],
   "source": [
    "qs = np.linspace(0.01, 1, num=90)\n",
    "prior_sigma = make_uniform(qs, name='std')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use `make_joint` to make the joint prior distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.918892Z",
     "iopub.status.busy": "2021-04-16T19:37:28.918336Z",
     "iopub.status.idle": "2021-04-16T19:37:28.920584Z",
     "shell.execute_reply": "2021-04-16T19:37:28.921030Z"
    }
   },
   "outputs": [],
   "source": [
    "from utils import make_joint\n",
    "\n",
    "prior = make_joint(prior_mu, prior_sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we'll start by working with the data from the control group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.924782Z",
     "iopub.status.busy": "2021-04-16T19:37:28.924226Z",
     "iopub.status.idle": "2021-04-16T19:37:28.926918Z",
     "shell.execute_reply": "2021-04-16T19:37:28.927340Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(38, 51)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "complete = (df['CENSOR'] == 0)\n",
    "data_complete = log_duration[complete]\n",
    "data_ongoing = log_duration[~complete]\n",
    "\n",
    "len(data_complete), len(data_ongoing)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Likelihood\n",
    "\n",
    "First update"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.931426Z",
     "iopub.status.busy": "2021-04-16T19:37:28.930801Z",
     "iopub.status.idle": "2021-04-16T19:37:28.936844Z",
     "shell.execute_reply": "2021-04-16T19:37:28.937332Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(90, 101, 38)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu_mesh, sigma_mesh, data_mesh = np.meshgrid(\n",
    "    prior.columns, prior.index, data_complete)\n",
    "\n",
    "mu_mesh.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.942522Z",
     "iopub.status.busy": "2021-04-16T19:37:28.941934Z",
     "iopub.status.idle": "2021-04-16T19:37:28.966243Z",
     "shell.execute_reply": "2021-04-16T19:37:28.965778Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(90, 101, 38)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "densities = norm(mu_mesh, sigma_mesh).pdf(data_mesh)\n",
    "densities.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.969673Z",
     "iopub.status.busy": "2021-04-16T19:37:28.968985Z",
     "iopub.status.idle": "2021-04-16T19:37:28.973058Z",
     "shell.execute_reply": "2021-04-16T19:37:28.972560Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(90, 101)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "likelihood = densities.prod(axis=2)\n",
    "likelihood.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.978511Z",
     "iopub.status.busy": "2021-04-16T19:37:28.977443Z",
     "iopub.status.idle": "2021-04-16T19:37:28.981763Z",
     "shell.execute_reply": "2021-04-16T19:37:28.981310Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(90, 101)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from utils import normalize\n",
    "\n",
    "posterior = prior * likelihood\n",
    "normalize(posterior)\n",
    "posterior.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second update"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.931426Z",
     "iopub.status.busy": "2021-04-16T19:37:28.930801Z",
     "iopub.status.idle": "2021-04-16T19:37:28.936844Z",
     "shell.execute_reply": "2021-04-16T19:37:28.937332Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(90, 101, 51)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu_mesh, sigma_mesh, data_mesh = np.meshgrid(\n",
    "    prior.columns, prior.index, data_ongoing)\n",
    "\n",
    "mu_mesh.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.942522Z",
     "iopub.status.busy": "2021-04-16T19:37:28.941934Z",
     "iopub.status.idle": "2021-04-16T19:37:28.966243Z",
     "shell.execute_reply": "2021-04-16T19:37:28.965778Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(90, 101, 51)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "densities = norm(mu_mesh, sigma_mesh).sf(data_mesh)\n",
    "densities.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.969673Z",
     "iopub.status.busy": "2021-04-16T19:37:28.968985Z",
     "iopub.status.idle": "2021-04-16T19:37:28.973058Z",
     "shell.execute_reply": "2021-04-16T19:37:28.972560Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(90, 101)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "likelihood = densities.prod(axis=2)\n",
    "likelihood.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:28.978511Z",
     "iopub.status.busy": "2021-04-16T19:37:28.977443Z",
     "iopub.status.idle": "2021-04-16T19:37:28.981763Z",
     "shell.execute_reply": "2021-04-16T19:37:28.981310Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(90, 101)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from utils import normalize\n",
    "\n",
    "posterior2 = posterior * likelihood\n",
    "normalize(posterior2)\n",
    "posterior2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:29.059387Z",
     "iopub.status.busy": "2021-04-16T19:37:29.048016Z",
     "iopub.status.idle": "2021-04-16T19:37:29.235472Z",
     "shell.execute_reply": "2021-04-16T19:37:29.235823Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtFUlEQVR4nO3deZxcZZ3v8c+3ujvppDv7AlkIhJ2gCBgRNxZFBJRhvIMDuI2MirvO6DjD3OsVxX28egWXYRiGCyiKK4IMCCjDACJIghgIa4iQhEBC9q076a763T/OqaS6UtV9OunqPkl/369XvbrO/junT9Wvnuc85zmKCMzMzPKmMNQBmJmZ1eIEZWZmueQEZWZmueQEZWZmueQEZWZmueQEZWZmueQElVOSNkk6cKjjaLRG7qekkyQtqxheKOmkAVr3OyTdVjEckg4eiHWn68vV/1/ShyStSOOaNNTx7CpJn5P0g11YbsDOHcvOCWqQSbpT0vv6mi8i2iNiccZ1DuiXY8ZtHpBut3l31tOf/dxdEXFkRNzZ2zxZ9ysiro2IUwcirlrnxGAel75IagG+CZyaxrV6qGMabFnOHRt4TlA2JHY3se3u8nvqtofIPkArsHCoA7HhxQlqCEl6v6RFktZIulHS9Ipp20tFkq6S9F1J/ylpo6T7JR2UTrsrXeRPafXLOTW28x5Jv5P0bUnrJT0u6Q0V06en21+TxvP+imnHSZonaUNaxfPNdFJ5u+vS7b4qnf9vJT0maa2kWyXtX7VPH5H0FPBUjf0cJ+kaSS9KelbSZyQVqvbh/0paA3yuxn6OSo/VWkmPAq+omv6MpFP6u1+1tp2Ou6cqhDMkLZa0StLXK2LvUa1UWUqT9CXgdcB30u19ZxeOyz2S/k+633+WdHrV/35xet78WdI7qo9bOt9ISd+StDx9fSsddyjwRMUxuaPGsuX9OV/S0jSOD0p6haQFktaV96uv41EntgslPZ3uw6OS3lq1f73t/2xJ/50uezswudY20nknS7opjXeNpLsrjnPluTNK0tXp9h6T9I/qWZX8jKRPp/u+WdJ/SNpH0i1pHL+RNKFi/p9KekHJZ/MuSUfWi3HYiQi/BvEF3Am8D3g9sAo4FhgJfBu4q2K+AA5O318FrAGOA5qBa4Hras1bZ5vvAbqBvwdagHOA9cDEdPp/A98j+ZV8NPAi8IZ02u+Bd6Xv24Hj0/cHpNttrtjOXwKLgCPSOD8D3FsV5+3ARGBUjf28BrgBGJOu/0ngvVX78LF03aNq7OdXgbvT9e8HPAIsq5j+DHDKLuzXTttOx91TtW//lW57Vhr7+9JpnwN+UDFvj22QnhNV+9Kf49IFvB9oAj4ELAcEtAEbgMPSeacBR9Y5Ry4G7gOmAlOAe4Ev1DsmVcuWp19Gcg6dCnQCv0zXNwNYCZyY5XjUWP/bgOkkP6jPATYD0/ra/4r/8zdJPmMnABsrt121na+k+9CSvl5XsZ5n2HHufJXkMzMBmAksYOfz7D6Skmd53x8EjknjuAO4qGL+v03/tyOBbwEPDfX3VF5eQx7AcHuxI0H9B/AvFePb0w/aAelwdYK6omLeM4DHK4azJKjtH9p03B+Ad5F8kReBMRXTvgJclb6/C/g8MLlqnTt9qQC3kH5xpsMFYAuwf0Wcr69aTwAHp18uW4E5FdM+ANxZsQ9L+ji2i4HTKoYvqPHFccou7NdO26Z2gqrc9oeB36bvP8cuJqiMx2VRxbTR6bL7kiSodcBfUSOhV23vaeCMiuE3Ac/UOyZ1jtmMinGrgXMqhn8O/F2W45HhM/QQcFaG/Z9F8sOirWL6D6mfoC4m+SGw02ep6txZDLypYtr7apxn76ja93+tGP4Y8Ms6MYxP4x+X5Vjs7S9X8Q2d6cCz5YGI2ETyoZ5RZ/4XKt5vIUlo/fFcpJ+A1LNpDNOBNRGxsWpaOY73AocCj0t6QNJbetnG/sAlaRXJOpJSn+i5T0vrLDsZGEHFMamKo7dly6ZXzfNsvRnp335l2Xb1POXju7uyHJft50ZEbEnftkfEZpISxweB55VUER9eZzs9zkd2Lf4VFe87agz395wFQNK7JT1UcV69hJ5VdTX3nyT+telxKOvtnPg6SQ3AbWm16IV15qs+z2qdG5mOhaQmSV9NqzA3kCQ36KUqcjhxgho6y0m+0AGQ1AZMAp5r0PZmSFLF8Kw0huXAREljqqY9BxART0XEeSRVNV8DfpbGWqsb/KXAByJifMVrVETcWzFPreUgqe7souKYVMbRx7Jlz5OUCCuXr6mf+5Vl29TY9vL0/WaSX/Zl+/Zj3VmOS10RcWtEvJGkeu9x4N/rzNrjfKRn/AOtr+OxnZJrmP8OfBSYFBHjSapuVW+ZCs8DE9L/a1lv58TGiPhURBwInAl8UhXXaqvWO7NieL8a82T1duAs4BRgHElpErLt317PCWro/BA4X9LRkkYCXwbuj4hndmFdK4C+7pmZCnxcUoukt5FcJ7o5IpaSXG/4iqRWSUeRlC6uBZD0TklTIqJEUl0ESZXgi0CparuXAf9cvsir5OL+27LsQEQUgZ8AX5I0Jv1i+iTQn3tWfpJuf4KkmSRVKTX1c7+y+nS67f2ATwA/Tsc/BJwgaZakccA/Vy1X9/+3O8clvTD/F+kX9FZgE8k+1vIj4DOSpkiaDHw2yzZ20UP0fjwqlX80vAgg6XySElSfIuJZYB7weUkjJL2WJPHUJOktkg5Of8htIDlWtY5X5Xk2gyR57qoxJP+b1SRJ+8u7sa69jhPU0IiI+C3wv0nqp58HDgLO3cX1fQ64Oq0C+es689wPHELyi/xLwNmx436W80h+uS0Hrie5gHt7Ou00YKGkTcAlwLkR0ZlWpXwJ+F263eMj4nqS0sh1aXXFI8D2FlUZfIzk1/Vi4B6SJH5lP5b/PEkVzp+B24Dv9zJv5v3qx/ZvAOaTfAH/J8l1RtJj+WOSi+nzgZuqlrsEODttFXZpjfXu6nEpAJ8i+b+uAU4kuTZWyxdJvswXAA+TXNT/YoZt9FuG41E576PAN0gaO6wAXgr8rh+bezvwSpL9v4ikwUk9hwC/IUnkvwe+F7XvfboYWEZynv0G+BlJktkV15Ccs88Bj5I0rrBUuYWKDRJJDwIXR8QvB3Gb7yG5CP/awdqm2XAh6UMkP3BOHOpY9jYuQQ2itOrrCOCPQx2Lme0aSdMkvUZSQdJhJKXU64c6rr1RwxKUpCslrZT0SJ3pknSpkhtDF0g6tlGx5IGkr5FUO/1TWjduZnumEcC/kdxTdQdJ1e73hjSivVTDqvgknUBSl3tNROx0UVPSGSR162eQ1BFfEhGvbEgwZma2x2lYCSoi7iK5MFnPWSTJKyLiPmC8pGmNisfMzPYsQ9np5Qx63uC2LB33fPWMki4g6RWAtra2lx9+eL17Dc3MbKjNnz9/VURM2d31DGWCqnUjWs36xoi4HLgcYO7cuTFv3rxGxmVmZrtB0oBcZx/KVnzL6HkH9kwad+e6mZntYYYyQd0IvDttzXc8sD4idqreMzOz4alhVXySfgScBExW8qyUi0i6sCciLgNuJmnBt4ik89PzGxWLmZnteRqWoNKOOHubHsBHGrV9MzPbs7knCTMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzyyUnKDMzy6WGJihJp0l6QtIiSRfWmD5O0q8k/UnSQknnNzIeMzPbczQsQUlqAr4LnA7MAc6TNKdqto8Aj0bEy4CTgG9IGtGomMzMbM/R3NcMkqYCrwGmAx3AI8C8iCj1sehxwKKIWJyu5zrgLODRinkCGCNJQDuwBuju706Ymdnep26CknQycCEwEfgjsBJoBf4SOEjSz4BvRMSGOquYASytGF4GvLJqnu8ANwLLgTHAObUSn6QLgAsAZs2a1edOmZnZnq+3EtQZwPsjYkn1BEnNwFuANwI/r7O8aoyLquE3AQ8BrwcOAm6XdHd10ouIy4HLAebOnVu9DjMz2wvVTVAR8elepnUDv+xj3cuA/SqGZ5KUlCqdD3w1IgJYJOnPwOHAH/pYt5mZ7eX6vAYFIOnNwJEkVXwARMTFfSz2AHCIpNnAc8C5wNur5lkCvAG4W9I+wGHA4myhm5nZ3ixLI4nLgNHAycAVwNlkKOFERLekjwK3Ak3AlRGxUNIH0+mXAV8ArpL0MEmV4D9FxKpd3RkzM9t7KKld62UGaUFEHFXxtx34RUScOjgh9jR37tyYN2/eUGzazMwykDQ/Iubu7nqy3AfVkf7dImk60AXM3t0Nm5mZ9SbLNaibJI0Hvg48SNIS74pGBmVmZtZngoqIL6Rvfy7pJqA1ItY3NiwzMxvusjSSaALeDBxQnl8SEfHNxoZmZmbDWZYqvl8BncDDQF/dG5mZmQ2ILAlqZkQc1fBIzMzMKmRpxXeLpCFpUm5mZsNXlhLUfcD1kgokTcwFRESMbWhkZmY2rGVJUN8AXgU8HH3d1WtmZjZAslTxPQU84uRkZmaDKUsJ6nngTkm3AFvLI93M3MzMGilLgvpz+hqRvszMzBouS08Snx+MQMzMzCpl6UniV+z8JNz1wDzg3yKisxGBmZnZ8JalkcRiYBPw7+lrA7ACODQdNjMzG3BZrkEdExEnVAz/StJdEXGCpIWNCszMzIa3LCWoKZJmlQfS95PTwW0NicrMzIa9LCWoTwH3SHqapBeJ2cCHJbUBVzcyODMzG76ytOK7WdIhwOEkCerxioYR32pgbGZmNozVTVCSXh8Rd0j6H1WTDkyfB/WLBsdmZmbDWG8lqBOBO4Aza0wLwAnKzMwapm6CioiL0r/nD144ZmZmiT5b8Un6hKSxSlwh6UE/H8rMzBotSzPzv42IDcCpwFTgfOCrDY3KzMyGvSwJSunfM4D/FxF/qhhnZmbWEFkS1HxJt5EkqFsljQFKjQ3LzMyGuyw36r4XOBpYHBFbJE0iqeYzMzNrmLolKEkHAEREKSIejIh16fDqiFiQNpqYOThhmpnZcNNbCerrkgrADcB84EWgFTgYOBl4A3ARsKzRQZqZ2fDT231Qb5M0B3gH8LfANGAL8BhwM/AlPwvKzMwapddrUBHxKPC/BikWMzOz7bK04jMzMxt0TlBmZpZLTlBmZpZLWe6DQtIMYP/K+SPirkYFZWZm1meCkvQ14BzgUaCYjg6gzwQl6TTgEqAJuCIidurDT9JJJA8+bAFWRcSJ2UI3M7O9WZYS1F8Ch0XE1v6sWFIT8F3gjST3Sj0g6ca0ZWB5nvHA94DTImKJpKn92YaZme29slyDWkxSuumv44BFEbE4IrYB1wFnVc3zduAXEbEEICJW7sJ2zMxsL5SlBLUFeEjSb4HtpaiI+Hgfy80AllYMLwNeWTXPoUCLpDuBMcAlEXFN9YokXQBcADBr1qwMIZuZ2Z4uS4K6MX31V61HckSN7b+cpNukUcDvJd0XEU/2WCjicuBygLlz51avw8zM9kJ9JqiIuFrSCJLSDsATEdGVYd3LgP0qhmcCy2vMsyoiNgObJd0FvAx4EjMzG9ayPPL9JOApkgYP3wOelHRChnU/ABwiaXaa4M5l55LYDcDrJDVLGk1SBfhY9vDNzGxvlaWK7xvAqRHxBICkQ4EfkVTN1RUR3ZI+CtxK0sz8yohYKOmD6fTLIuIxSb8GFpA8BPGKiHhk13fHzMz2FlkSVEs5OQFExJOSMrXqi4ibSXo+rxx3WdXw14GvZ1mfmZkNH1kS1DxJ/wF8Px1+B8nzoczMzBomS4L6EPAR4OMkLfPuIrkWZbbHiTptQFWrzamZDaksrfi2At9MX2Z7jAgopa9Ih+vdoyCSJFWoeJnZ0KqboCT9JCL+WtLD1PhcR8RRDY3MbBdEQLEiMQloKiTNVVXYkYiqlyknsFLAttKO5Zrk0pXZUOmtBPWJ9O9bBiMQs11VTjDFUpKcCkoSS0shW3KR0rvKlTQ3ba5Y39YSNDtRmQ2JuvdBRcTz6dsPR8SzlS/gw4MTnlnvSgHbitBVTBLIyCYY0ZSUfnY1oZSr+lrSdRVLyTZK7sPEbFBl6Sz2jTXGnT7QgZj1R6SJaVsxSUYjmpKSzkCXcgrakfC2FZNkZWaDo7drUB8iKSkdKGlBxaQxwO8aHZhZPcUSdJWSareRTY2vepOgOS1VbUufiNbkZ1GbNVxv16B+CNwCfAW4sGL8xohY09CozGqISBJTRFKqydLSrrOryMbOIl3dJbYVg67uEgG0NIkRTQVamsW4US2MbOk74xTS61pdpR3VgGbWOHUTVESsB9YD5wGkDxNsBdoltZef4WQ2GMpVekqr3OqVmiKCdVu6WbO5i/Ud3ZRKwZjWZkY2J8mobUQLEmwrlujqDjZvKfLsqk5GjSgwqX0Ek9pbGNFcP1mVS05dxd7jMLPdl+WR72eS3AM1HVgJ7E/SoeuRjQ3NLFFuCNFba7pSBKs2bmP5uq0UJKaMaWHfcW2MHlFAfWSRUgTrt3SzetM2lq3pZMaEkUwbP7LucuWSUzGSqj8za4wsPUl8ETge+E1EHCPpZNJSlVmjlZNTS6H+dZ9VG7fx7OoORrU0ccDkUYwb1dxnUqpUkJjQ1sKEthY6thV5euUW1nd0c+i+bTTVqMdT2sKvq+jm52aNlOVSb1dErAYKkgoR8V/A0Y0Nyyy95lQuOdU4U7uKJZ58YTPL1nRy6D5tzJnRzvjRLXWTU0TQVSxRqtffETBqRBNzZrTT3CQeW76JrjrN9spbcMtzs8bJUoJaJ6mdpA++ayWtBLobG5YNd+VrTgUlCaraui1dLFqxhcljRnDQ1NE7lXTWdnTx2IrNPLZyE0+v7qCjq0hXMWgqiOaC2H9CK4dOaWPOPm3sN661R1IrSBw8dTTPrOrgyRe2MGd6205JT0qSZrEEhaaGHAKzYS9LgjoL6AT+nqQn83HAxY0Myqw7bSlXKzm9uHEbz67q4JB9RjNudM8nvzz54mZ+9eiLrNi0jcOntDFnn3bOOnIq7SOaaGkq0FQQm7cVWbx6C0+8uJkr7n+OWeNbOe+YabSN2JFpJHHA5FEsWLqJNZu7mNQ+Yqc4Ckpa9JlZYyh6qe7Io7lz58a8efOGOgxroPJ9TrXucVqxfivL1nRyxIx2RlcklGXrO7lh4UpWbtrGmUdM4diZYylkuDjUVSxxw8KVPLR8I+9++XQOndLWY/razV08u7qDl+03ZqdSVARsLQ7OvVhmexJJ8yNi7u6up7cbde+JiNdK2kjPqnYBERFjd3fjZtXK9zrVasK9ZnMXy9Z0MmdGO6PS5BQR3P7Uau5YtIbTDpvMB46fQHM/blBqaSpw9lH7csQ+7Vz5wHO897iZHDJ59Pbp40c3s3SN2NDRvVNpzUnJrLF6uw/qtenfMYMXjg135c5eq3NMuXXd4dPatienYin46YIXWLy6gwtPns34UTsSSESw+MXN3PnoSuY/s5a2kc1MaBvBxLYRHL3/eI4/eFKP9R+5Tzt/8/LpXPXAc3z2jQcxMq1blMS4Uc1s6CzulKDMrLGy3Ad1CXBdRPx+EOKxYSwiufY0oql6fPDkC1vYb2IrY1qTU7ZYCq584Dk6u0r8/Qn7M6olWai7WOLf7ljMT+5fSlexxElHTOVNL92Xrd0l1mzaxprN27jwxws4cEobl777WNpbd3wEjtinnYMmj+a/n17DqYdN3j5+zKhmXli/tfEHwMx6yNJI4kHgf0s6FLge+HFE+CKQDbhSndLTig3baG4S+4zd0VDh+kdW0NlV4kOv3m97ld6K9Z184vt/pLmpwL+e/3KOmL7zdSOAT7/5MC76+ULOv/wPXPWB42gbueNj8KpZ4/j1E6t7JKiWJlEs7nytdg+7fGu2x+nzPqiIuDoizgCOA54EvibpqYZHZsNOuXqvUnexxLI1nRwwedT2ZPPQ8g08/MIm3nvcjO3J6e4nXuTMb9zDqw+dzNUfOI45M8bWvR+qpanAF89+CQdObed9VzxAZ7kHWOCAiaNYsq6DYtWzNWrloqD2AxDNbGD0p0/mg4HDgQOAxxsSjQ1b5afZNlV92b+4cRvjRjXTNjKtwisFP1uwgncdO317K74HFq/hk9f+iUvffQwfP/WQmr0/VCsUxFf++qV0dQd3PfHi9vGjWppobS6wYeuOW/2KpdgpLthR4jOzxugzQUkql5guBh4BXh4RZzY8MhtWyl/21aWRVRu7mDJmR9XevKXrmdo+goPTlnZbu4v8848X8MWzX7JTw4e+FAripfuNY8nqLT3Gd3aXtl/TAtiytcjokTvfjVss+bEbZo2U5RrUn4FXRcSqRgdjw1et0sjWrhKdXSXGjt5xmt6xaA1vfcnU7cPX3P0ss6e286aj9t2l7c6YOIplazq2D3d2lygFjKwoMm3aWmT86J4flVLsqOIzs8bI8vvvcuA0SZ8FkDRL0nGNDcuGm1pf9pu3FWlvbdp+w+3mbUVWbenisKk7bqa9dcELvPu1++/ydp9esYn9J+247+mxFZs4aNKO613dpWDdli7GVzUxL5bcUaxZo2VJUN8FXsWOHsw3puPMBkzEzl/2nduKtFY8SHDJ2g5mjW/dnrA2dHTx+PMbeMXsibu4zeD3i1bzyoqqwT8u38gxM3bcg7564zbGjur5jKiIpEGHq/fMGivLR+yVEfERkv74iIi1wM4dk5nthlolqK5i9EgMazu6mVRRklm5oZMpY0bSWn3jVEa/mPccrS1NHDE9uRd9xcatPL5yM0enw8VS8NzarUwb1/N0r3czsZkNrEyP25DURNrSVtIUwF1k2oASOzflLghKFc29W1sKbO3ecerNnDCa5es6d2oSnsXS1Vv48g2P8X/feTSSiAh+umAFpx46iTHpfVHL122lbWRTjx4kyjcT9/LQXTMbIFk+ZpeS3KA7VdKXgHuALzc0KjOSVnaVj2Ma1dKz+XfriCYmt4/g0ec29Gu9azdv4yNXPciHTjmII6Yn1Xn3LVnP2o4uTj4oqS7s2FbkhXVbOWDyqO3LlfsJbHLpyWxQZLlR91rgH4GvAM8DfxkRP210YDa8FJS0jKvUNrKJTRUJ6cCJo1m+YSvrOrq2j/voqYfwP3+yoO6DBas9u2ozf/3t3/PqQyfx3hNnA/Doik3csHAl7ztuJk0FUSwFT76wmf0mtTKy4hpYMZIk5dKT2eCo+1GTNLH8AlYCPwJ+CKxIx5kNmFoJakxrM5u3FrdX4Y1sLnDM9LHcv2T99nnOPX4/Jo8ZyRd/+Sjrt3RRz+pNW/mXmx7nrd/6He98zf5ceOYRSOKZNR1cPW8573/lTKaNHZl0MrtyC20jm3p0rVRKq/Za/GgNs0HT231Q89lx7XoWsDZ9Px5YAsxudHA2fJQf/lfZmq+pIMaNbmHF+q1Mn9AKwEkHTeDSe5bw0mljmD52JJL4l3OP4qJfLOR1X7iD4w6ayCkv2Yfp40exauNWXty4lSef38hvFq7g9JdN45ZPn8A+41opRXDn02v59ROreNex0zho0ui0B/QOOrtKzJnRvr2peSl9um9LwVV7ZoOpzwcWSroMuDEibk6HTwdOiYhPDUJ8O/EDC/detR7xvmVbkYXLNnH0/mNoSdt1379kHbc8vopPvG5/JlQ8YmNjZxe3P7yC3y5cybot25gydiRTxoxkxoRRnHnsdCa1jwRg9ZYuvj9/OcVS8K6XT2dq+wiKpWDRii10l4LDp7Vt7y6p/Oj5poKr9syyGqgHFmZJUPMj4uVV4+YNxMZ3hRPU3qtY2vG4jcpqtMUrt1AsBQfvM3p7qeb2J1fz20WredNhk3ntAeO3J6/evLBxK/cvWc+9z6zjDYdM5JRDJlGQ6Owq8tQLWxjZUuDgfUZvv8/Kycls1zT8iboVVkn6DPADkiq/dwKrd3fDZtXKffEVA5orEtT+k0fx6HOb+POqDmanvZq/8dBJHD61jZsee5HfPrWa4/cfz+yJo9i3fQQTRrcQAes6u1mzZRvPrd/KH5auZ11HN3NnjuVTJx7A1PYRRAQr1m9lyepOZk5sZd9xI3aq1msu7NyBrZkNjiwlqInARcAJJAnqLuDiiFjT58ql04BLgCbgioj4ap35XgHcB5wTET/rbZ0uQe3dyolhRFPP6z3dxaRlXaEAB00d3aPE9Oc1Hfxp+UaWrutkxaatbNpWJALGjGxi4ugWpraP4JjpYzl8alJ1FxGs7+hmyepOBBw4dfT23tIhKcl1lZJrTu4twqz/Bq0ElSaiT/R3xenNvd8F3ggsAx6QdGNEPFpjvq8Bt/Z3G7b3KV+D6ir2rOprbhKHT2/j2VUd/GnJRvab2MrkMSNoKojZE0cxe+KO+5W2dpdoKmj7s6LKukvBqo3beGH9VrqLwX6TWpnY1rK91FS+CbcUOydIMxt8War4dtVxwKKIWAwg6TrgLODRqvk+BvwceEUDY7E9SPl6z7Zi0qy7nCgKErOnjGbKmKT08+zqDsaNamFiWwtjRzXT3CQKSpqjRwRbu0ps7S6xZVuRtZu72NjRzZhRzew7biST2nsmplLFTbjV18DMbGg0MkHNAJZWDC8DXlk5g6QZwFuB19NLgpJ0AXABwKxZswY8UMufHkmq0PNZUe2tzcyZ0U5XscTazV2s3tzFs6s76E6fgdFUEMUImguitaVAa0uBKWNGcMi+bTuVqkqRlNbApSazvGlkgqr1Ua++4PUt4J8ioljv8dwAEXE5yWM/mDt3bv87XrM9UnOamLqKSXJqrroPqaWpwNSxI5k6duT2caVSUCwFhYLqPlm33Bt5ufOJpoIfnWGWR3UTlKRvs3NC2S4iPt7HupcB+1UMzwSWV80zF7guTU6TgTMkdUfEL/tYtw0ThbTKrRg77pMq94VXK6EUCqJQIzGVq/GKsePhiM2F+usxs6HXWwmq3FTuNcAc4Mfp8NtIepnoywPAIZJmA88B5wJvr5whIrb3RiHpKuAmJyerJiXNzpu0I8l0ldLkQu0kE+kTb8t/y0mpoKTK0EnJLP/qJqiIuBpA0nuAkyOiKx2+DLitrxVHRLekj5K0zmsCroyIhZI+mE6/bPfDt+FEaZJqYkeJKEg7cS3tPK/SvwWclMz2RFmuQU0HxgDl+57a03F9SrtHurlqXM3EFBHvybJOM9iRrMxs75UlQX0V+KOk/0qHTwQ+17CIzMzM6CNBSSoAT5A0Dy83Eb8wIl5odGBmZja89ZqgIqIk6RsR8SrghkGKyczMLNMj32+T9Ffq7UYlMzOzAZblGtQngTagW1InSeOoiIixDY3MzMyGtSydxY4ZjEDMzMwqZerqSNIE4BCgtTwuIu5qVFBmZmZ9JihJ7yN53MZM4CHgeOD3JB28mpmZNUSWRhKfIOlp/NmIOBk4BnixoVGZmdmwlyVBdUZEJ4CkkRHxOHBYY8MyM7PhLss1qGWSxgO/BG6XtJadeyU3MzMbUFla8b01ffu5tLujccCvGxqVmZkNe709D2pijdEPp3/b2dF5rJmZ2YDrrQQ1n+RpBgJmAWvT9+OBJcDsukuamZntprqNJCJidkQcSPI8pzMjYnJETALeAvxisAI0M7PhKUsrvlekz3UCICJuIXnkhpmZWcNkacW3StJngB+QVPm9E1jd0KjMzGzYy1KCOg+YAlxP0tR8ajrOzMysYbI0M19D0puEmZnZoMnSF9+hwD8AB1TOHxHui8/MzBomyzWonwKXAVcAxcaGY2ZmlsiSoLoj4l8bHomZmVmFLI0kfiXpw5KmSZpYfjU8MjMzG9aylKD+Jv376YpxARw48OGYmZklsrTic5dGZmY26LI+8v0lwBx6PvL9mkYFZWZmlqWZ+UXASSQJ6mbgdOAewAnKzMwaJksjibOBNwAvRMT5wMuAkQ2NyszMhr0sCaojIkpAt6SxwErcQMLMzBosyzWoeekj3/+d5BlRm4A/NDIoMzOzLK34Ppy+vUzSr4GxEbGgsWGZmdlw12cVn6Tflt9HxDMRsaBynJmZWSPULUFJagVGA5MlTSB53DvAWGD6IMRmZmbDWG9VfB8A/o4kGc1nR4LaAHy3sWGZmdlwVzdBRcQlwCWSPhYR3x7EmMzMzOpfg5L0Ckn7lpOTpHdLukHSpVk7i5V0mqQnJC2SdGGN6e+QtCB93SvpZbu+K2ZmtjfprZHEvwHbACSdAHyVpPeI9cDlfa1YUhNJVeDpJL1QnCdpTtVsfwZOjIijgC9kWa+ZmQ0PvV2Dakof9w5wDnB5RPwc+LmkhzKs+zhgUUQsBpB0HXAW8Gh5hoi4t2L++4CZ/YjdzMz2Yr2VoJoklRPYG4A7KqZlucF3BrC0YnhZOq6e9wK3ZFivmZkNA70lmh8B/y1pFdAB3A0g6WCSar6+qMa4qDmjdDJJgnptnekXABcAzJo1K8OmzcxsT9dbK74vpTfkTgNui4hycikAH8uw7mXAfhXDM4Hl1TNJOgq4Ajg9IlbXieVy0utTc+fOrZnkzMxs79JrVV1E3Fdj3JMZ1/0AcIik2cBzwLnA2ytnkDQL+AXwrn6s18zMhoFMDyzcFRHRLemjwK1AE3BlRCyU9MF0+mXAZ4FJwPckAXRHxNxGxWRmZnsO7ai52zPMnTs35s2bN9RhmJlZHZLmD0RhI8vzoMzMzAadE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSE5SZmeWSImKoY+gXSRuBJ4Y6jn6aDKwa6iB2wZ4Yt2MePHti3I55cBwWEWN2dyXNAxHJIHsiIuYOdRD9IWnenhYz7JlxO+bBsyfG7ZgHh6R5A7EeV/GZmVkuOUGZmVku7YkJ6vKhDmAX7Ikxw54Zt2MePHti3I55cAxIzHtcIwkzMxse9sQSlJmZDQNOUGZmlku5SVCS9pP0X5Iek7RQ0idqzCNJl0paJGmBpGMrpp0m6Yl02oU5ivkdaawLJN0r6WUV056R9LCkhwaqWeYAxXySpPVpXA9J+mzFtEE/zv2I+9MVMT8iqShpYjptKI51q6Q/SPpTGvPna8yTt3M6S8x5O6ezxJzHczpL3Lk6pyviapL0R0k31Zg2cOd0ROTiBUwDjk3fjwGeBOZUzXMGcAsg4Hjg/nR8E/A0cCAwAvhT9bJDGPOrgQnp+9PLMafDzwCTc3icTwJuqrHskBznrHFXzX8mcMcQH2sB7en7FuB+4PiqefJ2TmeJOW/ndJaY83hO9xl31fxDfk5XbPuTwA/rHNMBO6dzU4KKiOcj4sH0/UbgMWBG1WxnAddE4j5gvKRpwHHAoohYHBHbgOvSeYc85oi4NyLWpoP3ATMbHVdvMh7neobkOMMuxX0e8KPBiK2e9DzdlA62pK/qVkl5O6f7jDmH53SW41zPUJ7T/Y17yM9pAEkzgTcDV9SZZcDO6dwkqEqSDgCOIflFUWkGsLRieFk6rt74QdNLzJXeS/LLoiyA2yTNl3RBA8OrqY+YX5VWPdwi6ch03JAfZ+j7WEsaDZwG/Lxi9JAc67Qq5CFgJXB7ROT+nM4Qc6VcnNMZY87dOZ31WOfpnAa+BfwjUKozfcDO6dx1dSSpneSf8HcRsaF6co1Fopfxg6KPmMvznEzyYX5txejXRMRySVOB2yU9HhF3NT7iPmN+ENg/IjZJOgP4JXAIQ3ycIduxJqkK+V1ErKkYNyTHOiKKwNGSxgPXS3pJRDxSMUvuzukMMQP5OqczxJzLczrrsSYn57SktwArI2K+pJPqzVZj3C6d07kqQUlqIfnyuTYiflFjlmXAfhXDM4HlvYxvuAwxI+kokuLwWRGxujw+Ipanf1cC15MUgRuur5gjYkO56iEibgZaJE1mCI8zZDvWqXOpqgoZqmNdsf11wJ0kv4Ir5e6cLusl5tyd0xXbX0eNmPN6TlfEt446xzqVl3P6NcBfSHqGpIru9ZJ+UDXPwJ3TvV2gGswXSXa9BvhWL/O8mZ4X3/6Qjm8GFgOz2XHx7cicxDwLWAS8ump8GzCm4v29wGk5iXlfdtzEfRywJF1uSI5z1rjT+cYBa4C2HBzrKcD49P0o4G7gLVXz5O2czhJz3s7pLDHn8ZzuM+68ndNVcZ1E7UYSA3ZO56mK7zXAu4CH0zpZgP9J8mEgIi4DbiZpIbII2AKcn07rlvRR4FaSliJXRsTCnMT8WWAS8D1JAN2R9Ey8D0mRHpJ/3A8j4tc5ifls4EOSuoEO4NxIzrChOs5Z4wZ4K3BbRGyuWHaojvU04GpJTSS1FT+JiJskfbAi5ryd01lizts5nSXmPJ7TWeKGfJ3TNTXqnHZXR2Zmlku5ugZlZmZW5gRlZma55ARlZma55ARlZma55ARlZma55ARlVoekkPT9iuFmSS/W6sG5Adv+lqQTBmhd10k6ZCDWZTaYnKDM6tsMvETSqHT4jcBzjd6okscpHB8D123Nv5L0nWa2R3GCMuvdLSR3xkNVb9KS2iRdKemB9Nk4Z6XjD5B0t6QH09er0/EnSbpT0s8kPS7pWqV3WlY5G/h1xXaekfRlSb+XNE/SsZJulfR0+QbJdN03VSzzHUnvSQfvBk6RlKcb88365ARl1rvrgHMltQJH0bMH9f9F8nyeVwAnA1+X1EbSM/UbI+JY4Bzg0opljgH+DphD8lyc19TY5muA+VXjlkbEq0iSzVUkSex44OK+diAiSiR39b+sr3nN8sS/qMx6ERELlDze4zySLlwqnUrSceY/pMOtJF0vLQe+I+looAgcWrHMHyJiGUDaZdMBwD1V650GvFg17sb078MkD7nbCGyU1Jn2hN2XlcB0dk58ZrnlBGXWtxuB/0PSOeakivEC/ioinqicWdLngBUkJZYC0FkxeWvF+yK1P4MdJMmuUnm5UtU6Suk6uulZI1K9fGu6XrM9hqv4zPp2JXBxRDxcNf5W4GPl60iSjknHjwOeT6vW3kXSMWZ/PAYc3M9lngXmSBopaRzwhqrphwKD1Qmq2YBwgjLrQ0Qsi4hLakz6AsljuhdIeiQdBvge8DeS7iNJDJtrLNub/yQprfUnxqXAT4AFwLXAH8vTJO0DdETE8/2Mw2xIuTdzsxySdA/Js4HWDcC6/h7YEBH/sduBmQ0il6DM8ulTpM+6GgDrgKsHaF1mg8YlKDMzyyWXoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJf+Pw9hgF2DBhreAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from utils import plot_contour\n",
    "\n",
    "plot_contour(posterior2, cmap='Blues')\n",
    "\n",
    "decorate(xlabel='Mean (mu)', \n",
    "         ylabel='Standard deviation (sigma)',\n",
    "         title='Joint posterior distributions of mu and sigma')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior Marginal Distributions\n",
    "\n",
    "I'll use `marginal`, which we saw in <<_MarginalDistributions>>, to extract the posterior marginal distributions for the population means."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:29.239347Z",
     "iopub.status.busy": "2021-04-16T19:37:29.238897Z",
     "iopub.status.idle": "2021-04-16T19:37:29.243057Z",
     "shell.execute_reply": "2021-04-16T19:37:29.242677Z"
    }
   },
   "outputs": [],
   "source": [
    "from utils import marginal\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what they look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:29.283479Z",
     "iopub.status.busy": "2021-04-16T19:37:29.263758Z",
     "iopub.status.idle": "2021-04-16T19:37:29.391453Z",
     "shell.execute_reply": "2021-04-16T19:37:29.390972Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA52ElEQVR4nO3deZwcZ33v+8+ve/bRMlpGu2TJtmws70YYY7PYrDabw4Fc7AQceJ0TX19wCFkul+RwWU6S++Lcy0kICcHH7A6LAxwCBgw2wZjN2Ja8W9iSZVm2lpE0WmbRLJqZ7t/9o2q6q1s9M9U9XdM9mu/79ZqXq6ueqnq6XJrfPE/96nnM3REREak3qVpXQEREpBQFKBERqUsKUCIiUpcUoEREpC4pQImISF1SgBIRkbqkACVzhpn9oZndneDxv2Jmfxsuv8LMtlfx2D82sz8Kl99jZr+u4rETvS7lMrNWM/uBmfWa2bdrXR+pHQUoSYSZ7TazITM7bmYHzezLZjZvGsf7uJl9bTp1cvevu/vrp3OMMs71K3c/e6pycb+Xu1/j7l+dbr3MbL2ZuZk1RI49Y9clpncAy4El7v77ta6M1I4ClCTpLe4+D7gEeAnwkVpVJPoLuYJ9zcxq8m+llueuodOAHe4+VuuKSG3NtRtfasDd9wE/Bs4DMLO3mtk2M+sxs3vN7Jzxsmb2f5nZPjPrN7PtZvYaM7sa+GvgnWGL7LGw7EIz+6KZdYX7/K2ZpcNt7zGz35jZP5jZUeDjxV1jZna5mW0Ju5K2mNnlkW33mtnfmdlvgEHg9OLvZWYXm9nDYV3/DWiJbLvSzPZW+L1OOne47r8Unt7+Kaz702b2msiG3Wb22sjnaCvtl+F/e8JzvqyC6/I34bXtN7O7zWxpuK3FzL5mZkfC/7dbzGx5qXvCzM4Jj9UT3gtvDdd/Avho5Jr85xL7ftzMvh2eq9/MnjCzs8zsr8zskJntMbPXR8pPdj2kjilASeLMbC3wRuARMzsL+CbwQaATuBP4gZk1mdnZwM3AS9x9PvAGYLe7/wT4f4B/c/d57n5heOivAmPAmcDFwOuB6C/xlwK7gGXA3xXVaTHwI+AzwBLg74EfmdmSSLF3AzcC84Hni/ZvAr4H/CuwGPg28PYJvn+532vScxd9t6XAx4Dvht9pKq8M/9sRnvO3RXWNc13+AHgvwXVtAv4yXP9HwEJgbbjvTcBQcQXMrBH4AXB3eIw/Ab5uZme7+8covCZfnOB7vIXg2i8CHgHuIvh9thr4b8D/jHEtpM4pQEmSvmdmPcCvgV8Q/OJ5J/Ajd/+pu48CnwJagcuBDNAMbDKzRnff7e7Pljpw+Jf5NcAH3X3A3Q8B/wBcFym2393/yd3H3L34F+WbgGfc/V/D7d8Enib4xTfuK+6+Ldw+WrT/ZUAj8Gl3H3X37wBbJrgOsb9XzHMDHIqc+9+A7eF3mq441+XL7r4jvKbfAi4K148SBKYz3T3j7g+5e1+Jc1wGzAM+6e4j7n4P8EPg+jLq+St3vyvsBvw2wR87nwyv1e3AejPrKON4UocUoCRJv+fuHe5+mru/L/yFtopIi8Dds8AeYLW77yRoWX0cOGRmt5vZqgmOfRpBgOgKu4l6CP5qXhYps2eSuhXUI/Q8wV/gcfff54WjLZdq6VDm94pzbiY491THjCPOdTkQWR4kCDYQtGjuAm43s/1m9v+GraVS59gT/r+f6BxTORhZHgIOu3sm8plIvWSWUoCSmbafILgAwYMUgi6hfQDu/g13f3lYxoH/HhYtHnZ/D3ACWBoGwQ53X+Du50bKTDZUf0E9QuvG6xFj/y5gdVj/6P4llfG94pybCc69P1weANoi21aUcdw416WksDX3CXffRNAifjNwwwTnWGuFyR+xzlGhya6H1DEFKJlp3wLeFCYJNAJ/QRBo7jOzs83s1WbWDAwT/CU8/lfxQYJumxSAu3cRPMP4H2a2wMxSZnaGmb0qZj3uBM4ysz8wswYzeyewiaCrKY7fEjz/+kC4/38CLi1VsJzvVYZl4bkbzez3gXPC7wTwKHBduG0zQdr2uG4gS4mkj1DF18XMrjKz8y1IVOkj6PLLlCj6AEHQ+FBYxysJuhBvn+ocFXqUia+H1DEFKJlR7r4deBfwT8Bhgl9Mb3H3EYLnNJ8M1x8g+CX81+Gu4y9sHjGzh8PlGwge0v8OOAZ8B1gZsx5HCP7C/wvgCPAh4M3ufjjm/iPAfwLeE577ncB3Jyhe7veK4wFgY3jMvwPeEX4ngP8bOCOs1yeAb0TqPRiW/03YNXpZ0feaznVZQfD/oA94iuC540nZcuG1eyvBM8TDwL8AN7j707G+efkmvB5S30wTFoqISD1SC0pEROqSApSIiNQlBSgREalLClAiIlKXKh5Asx4tXbrU169fX+tqiIhIGR566KHD7t5ZvP6UClDr169n69atta6GiIiUwcxKjsKiLj4REalLClAiIlKXFKBERKQuKUCJiEhdUoASEZG6pAAlIiJ1SQFKJMLd8bGRWldDRDjF3oMSmQ7PjDF27y1k9z1Jw8W/R/r8q2tdJZE5TS0okVDmyZ+Q3fcE4Iw9/kO1pERqTAFKBMj27CfzxJ35FZlRvOup2lVIRBSgRDybJXPfbZAtnJ08+8KjtamQiAAKUCJkn/4Z2cPPhZ8sv37PY3hR0BKRmaMAJXOaHz/C2CPfz31OX/hmrHVhsG1kAD/0bK2qJjLnKUDJnJZ9bgtkRgGwRWtIn3c1qbUX5bfvebQ2FRMRBSiZ27I9+3LL6bNeiaUbSK27KL99z6O4ew1qJiIKUDKneU9Xbtk6VgX/XX4W1tgSbD9+BI8EMRGZOQpQMme5O953IPfZOlYG/003YKvPz61XNp9IbShAydx1/Ej++VPLfKx5Xm5Tau2FuWUFKJHaUICSOct79+eWbeHKgm2p1edBKh2UO7YHHzg6o3UTEQUomcMKnj8VBShraiXVeUa+7NE9M1YvEQkoQMmc5b3RBImVJ223hSvyZfu7Z6ROIpKnACVzVkGAWlgiQM3vzJftOzQjdRKRPAUomZPcnWyJFPMoW7A8X75fAUpkpilAydw0eAzGTgBgTe3QMv/kMvOX5RbVghKZeQpQMicVJkiswMxOKmPzljA+eKwPHMXDlHQRmRkKUDInTZUgAWANTVj74vE9gvemRGTGKEDJnDRVgkRuW0GixMFE6yQihRINUGZ2tZltN7OdZvbhEttfZGa/NbMTZvaX5ewrMh2lxuArxRZEnkMp1VxkRiUWoMwsDXwWuAbYBFxvZpuKih0FPgB8qoJ9RSri7kUtqBUTli0IUEqUEJlRSbagLgV2uvsudx8BbgeujRZw90PuvgUofvo85b4iFRvuw0cGAbCGFmhbNGFRi2byKdVcZEYlGaBWA9HxYfaG66q6r5ndaGZbzWxrd7e6YGRqhd17K0tm8OUo1VykZpIMUKX+1ced+S32vu5+q7tvdvfNnZ2dpYqIFIibIAFKNReppSQD1F5gbeTzGmD/BGWrua/IpIpbUJNRqrlI7SQZoLYAG81sg5k1AdcBd8zAviJ09wzyxR89yq0/eIThkbGCbeW0oKA4UUKp5iIzpSGpA7v7mJndDNwFpIEvufs2M7sp3H6Lma0AtgILgKyZfRDY5O59pfZNqq5y6hgeGeO7v3iaO+57htGxDABtzQ286/X5GXJ9IN8Kir7nNBGb3wldTwX7KlFCZMYkFqAA3P1O4M6idbdElg8QdN/F2ldkMs919fC3t/2anuPDBet/9vBurnvNuTSkU0GK+WBvfuMkGXzjCltQSsQRmSkaSUJOGZ/7/kMnBSeAvoETPLzjQPDhxABkgy4/a2zBGpunPK5SzUVqQwFKTgk9x4d5dt8xAFKpFDe/bTNve8XZue0/e+g5AHzwWH6nGK0nQKnmIjWiACWnhMd25pMXXrRuCVddsp7XvHhDbt3DOw5wtG8omGYjZG0dsY5t85eiVHORmacAJaeER57JB6iLzgwmGly5ZB6b1gdJEFl37n30eXywJ1cudoBKNxakmnv/4WpUWUSmoAAls56789iz+QB18cb8TLivffH63PLPH34eH6igi4/CRAn0HEpkRihAyay3+0AvfQPB7LgL2pvZsLIjt+2yTatpbW4EYP+Rfg7s25fbFrcFBUXTbihAicwIBSiZ9R555kBu+YIzlhWMrdfc1MDLz88PSvL8c8/nlisPUOriE5kJClAy6z0aef508ZknT53xmkg3X9/hQ/lRHcvo4iP3DIrCbkIRSYwClMxqwyNjPL0nPzLEBWcsO6nMmasX0dIUvJPePNafG2HC2st4BhUJZj7UU2FtRaQcClAyqz2x6xCZTBaAdcsXsnhB60llzIz1Kzpo8DGafYTBE6OQSkPzvNjnKQhmg2pBicwEBSiZ1R7bmU9YuPjM5ROW27Cyg3kEkxQODo9hrR2TzwNVrGUBuXehhvrxzNjk5UVk2hSgZFZ7dGc+QeLCSQPUwnyAOjFaVoIEgKUbsNYF4SeHod5Jy4vI9ClAyax1uHeQriPHAWhsSHPOaUsnLHv6qkXMZwCAweFRKDNAQWE3n6ubTyRxClAya71wsC+3fMaqRTQ1picsu6ZzPgssGEj2xGiGkcb4z59yookSyuQTSZwClMxaew7lA9TaZQsmKRm0sNbO99znQ8NNZZ+voFtQLSiRxClAyaxVToACWN2eT2zYNzBxa2siBanmkTH9RCQZClAya0UD1LrlUweoZc35Ucif65uk4ETa9AxKZCYpQMms5O4FAWpN59QBanE6P5nhjiM+ScnSCrr49AxKJHEKUDIrHe4d4sRo0GU3r7WJjnmTz4zr2QzzLB+gnjmcyY0oEVdhFl9PWfuKSPkUoGRWKn7+NOVLt8N9pFPQ3JhmyFoYdWNvd395J21dmFv0oV48W16AE5HyKEDJrFRu9954i6etpZHjBMMh7dpfXjedNTRhLfPDA2ZhuJIHWSISlwKUzErlZvDlAlRzI/20AcE8UmVr7cgfU8+hRBKlACWz0t7u8gIUuRZUA8fDALVrf0/Z59VzKJGZowAls05xBl+8FlTQ2mlrbswFqOe6enAvL5sv+i6UXtYVSZYClMw6h3uHGB6Jn8EH+dZOY0OabJjscGJ0jANHB8o7eSTVXF18IslSgJJZpzhBIta0GePdcQbzF+dHPd93uLxMPnXxicycRAOUmV1tZtvNbKeZfbjEdjOzz4TbHzezSyLb/szMtpnZk2b2TTNrSbKuMnuU270Hha2d+UvzAWp8NPS41MUnMnMSC1BmlgY+C1wDbAKuN7NNRcWuATaGPzcCnwv3XQ18ANjs7ucBaeC6pOoqs0u5CRLuXtDaWbRsRW65/ADVkT+uuvhEEpVkC+pSYKe773L3EeB24NqiMtcCt3ngfqDDzFaG2xqAVjNrANqA/QnWVWaRsltQo0OQGQmW000sW7Ykt+lAmQGqYDy+ofKTLEQkviQD1GpgT+Tz3nDdlGXcfR/wKeAFoAvodfe7E6yrzBKVZfDl33eytoWsWDo/97nsFlRjM9YUZAGSzcBwmaNRiEhsSQaoUk+ui//cLFnGzBYRtK42AKuAdjN7V8mTmN1oZlvNbGt3d/e0Kiz170hf+Rl80enZrXUhKxa1Y+Gt190zyFgmW14lNKq5yIxIMkDtBdZGPq/h5G66icq8FnjO3bvdfRT4LnB5qZO4+63uvtndN3d2dlat8lKforPoxs3g82grp3UhTY1pFi8Icm4c5+Cx8lLNCycu7ClrXxGJL8kAtQXYaGYbzKyJIMnhjqIydwA3hNl8lxF05XURdO1dZmZtFvwGeg3wVIJ1lVmi7BEkoKgFFeyzcsk0uvna1YISmQmJBSh3HwNuBu4iCC7fcvdtZnaTmd0UFrsT2AXsBD4PvC/c9wHgO8DDwBNhPW9Nqq4ye0RfrF21dF6sfXwoMqhry3iAyu9bboAq6OJTJp9IYhqSPLi730kQhKLrboksO/D+Cfb9GPCxJOsns080627F4ngBqvgZFBQGqHIz+UzPoERmhEaSkFkl+rxo+eL2WPtEW1DjXXwrIvtO510oPYMSSY4ClMwamUyWQz2Duc8rFsUNUJFpNcZbUJFU8wNHy+ziiz6DGjha3r4iEpsClMwah/uGyGaDlPCOeS00N8XsoR6OtqCCwBQNboeOlZdqXjCaxKBe1hVJigKUzBoHj1bQvZfN4MORFlJzEKCaGtMsWRDMrFt2qnljK6SbguXMaDBShYhUnQKUzBrRrrjYCRInjjP+frg1z8PS+VZXpanmZnZSK0pEqk8BSmaNghZURc+fCt+bml4m38L8BwUokUQoQMmsUdiCihegKMjgW1iwaTqZfIXDHfWUt6+IxKIAJbNG9CXduF18BS/pTtaCKjOTT118IslTgJJZwb0wkWHFkrgv6UZaUC2FAWrFNEaTKGiNKUCJJEIBSmaF/sERhk6MAtDc2MCCtqZY+/nwyePwjVsZaYWVm2qOWlAiiVOAklmheASJOKOYQ3EXX+EzqOmkmhd08Q31xN5PROJTgJJZIZrBtzJuijmUHMk8KppqXk4mX3Q8PnXxiSRDAUpmha5IEkPcl3Sh9EjmURVn8kWCnQ/14dlM/H1FJBYFKJkVoll2cd+BAiZNM4fiTL4yuvjSjVjz+L6uqd9FEqAAJbPCoWORQWJjZvD52Ag+PgyRpaD55MC2rGBMvvJm1lWihEiyFKBkVog+H4rdgoq0aqy19PTw0S4+Tf0uUl8UoKTujYxmONoftIRSZnR2tMXar9Q0G8WiLaiDxwbKGpk82mWoTD6R6lOAkroXbdksXdhGQzrmbRttQZVIkACY19pEW0sjAKNjGXqOn4hfMQ13JJIoBSipe4UjSJSTwTd5ivm4ZR2VPYfSgLEiyVKAkrpX+Pwp/jtQk43DFxV9DlVegOrIn0sBSqTqFKCk7hWPIhFbQQuq9DMoKHwOdaCcRAkFKJFEKUBJ3TtYMIp5GV18w9GXdOdPWG55hanmGk1CJFkKUFL3Cp5BlTXM0eQv6Y4rfBdqcMJyJ2meB6k0AD46hI+NxN9XRKakACV1rXiajWUxU8whfpLE8oJ3ocqc+r1g2o1jsfcVkakpQEldO9Y/zOhYMM5de0sT7a0xp9lwn3Qk86jOhfmgd7hnqKxpNwoTJXonLigiZVOAkrp2qNIEidFhyATzR5FugobmCYs2NaZZND8/7cbh3jK6+ZQoIZIYBSipawXde2UNElvYvTfV/FHRRImD5Qwaqy4+kcQkGqDM7Goz225mO83swyW2m5l9Jtz+uJldEtnWYWbfMbOnzewpM3tZknWV+lSQYl7W86d4CRLjli3KH7vSTL6CoZVEZNoSC1BmlgY+C1wDbAKuN7NNRcWuATaGPzcCn4ts+0fgJ+7+IuBC4Kmk6ir1q9J3oApSzCdJkBhXkMnXoy4+kXqQZAvqUmCnu+9y9xHgduDaojLXArd54H6gw8xWmtkC4JXAFwHcfcTdexKsq9SpaNp3eV18kRbUBOPwRa2IvqxbcRdfT+z9RGRqSQao1cCeyOe94bo4ZU4HuoEvm9kjZvYFMyv528nMbjSzrWa2tbu7u3q1l7pwsMKJCgtGFy+3BaXhjkTqwqQByszujiz/VZnHLvVUungug4nKNACXAJ9z94uBAeCkZ1gA7n6ru292982dnZ1lVlHq2Vgmy9G+YQCM+NNsAIUtqOi8TROIdh8e6ilnuKPoM6iesqbrEJHJTdWCiv7G//0yj70XWBv5vAbYH7PMXmCvuz8Qrv8OQcCSOaS7ZxAP/6ZZvKCFxoZ07H2j7yRN9pLuuMXzW0mH03j0DZxgeGQs1nmssRlrbAk+ZDNwosxZeUVkQlMFqOn8ObgF2GhmG8ysCbgOuKOozB3ADWE232VAr7t3ufsBYI+ZnR2Wew3wu2nURWahwgSJMoY4oqiLL0YLKpWyglEqykk1L5wXSqnmItXSMMX2083sDoKuuPHlHHd/60Q7uvuYmd0M3AWkgS+5+zYzuyncfgtwJ/BGYCcwCLw3cog/Ab4eBrddRdtkDjhU8A5UGd17AIPxRjKPWtbRTlc4tcfBYwOctiLeftbWgfd2hec9BovXTr6DiMQyVYCKZt19qtyDu/udBEEouu6WyLID759g30eBzeWeU04d0VZMWQkSmVF8ZHxfg+aJRzKPWr64HZ4Nz13pu1ADR2PvJyKTmzRAufsvxpfNrDNcp1Q5mREFXXyVppi3LsBS8ZJVK83ko31xbtEH1MUnUi1TZfGZmX3MzA4DTwM7zKzbzD46M9WTuSyaTVfOO1AFIzrEeP5U6hxlpZpHA5SeQYlUzVR/Wn4QeDnwEndf4u6LgJcCV5jZnyVdOZnbKu7ii47DF+Ml3XErKh2Prz0ycaG6+ESqZqoAdQNwvbs/N77C3XcB7wq3iSRiYGiE40PBBIAN6TSL5rfE3zkaoNriJTpA0bxQPQPx32nSMyiRREwVoBrd/XDxyvA5VGMyVRIpHA9vWUfblKORRxWM6NDaEXu/ea1NtLUEt/XIaIZj/cOx9ivs4tPLuiLVMlWAmmwOa81vLYmpdJBYoGgUifgtKCicUj7ucyhrbMaawjT4bAaiA9WKSMWmClAXmlmfmfWHP33jn4HzZ6KCMjcdqnQeKApbUHHfgRq3vMJBYwu7+ZQoIVINU6WZxx9bRqSKKk2QgKIsvjID1IpIa60rMlDtVKx9Md6zLzj/wFFYur6s84rIySYNUGbWAtwEnAk8TjAaRLxBykSmoTDFvMxRJCpMkoDCLr7yMvnyz6GUySdSHVN18X2VYDSHJwiGJPofiddIhMLgsKKMcfg8m8GHx1s+BmWkmUPh864DZbWgNB6fSLVNNdTRJnc/H8DMvgg8mHyVZK5z95Oy+GIb7mN8jGNrmYelyuulLmhBlTOahFLNRapuqhbU6PiCuvZkphzrH2Z0LANAe0sT7a1NsfeNTrNRTor5uCULCqfdGDoxOsUegcIuPrWgRKohbhbfeObeBZGsPuXSSiKiXWvlp5iXNw9UsVTKWN5R/ogSBe9CqQUlUhWTBih3T7v7gvBnvrs3RJbL/9cvEkM0vXvlknLngYomSHRUdP7llWTytXUwPkG0D/XhGXU4iExXvGGeRWbQgSP5oLCy3IkKC0aRKC+Db1zhy7qDk5TMs3RDpMXmBS05EamMApTUnf2RALViOl18ZaaYj6tKJp+6+USmTQFK6k40KJTdxTdY+Uu64wpHk4gfoDSahEh1KUBJXXH33LTrACuXxpsNN7f/UPlTvRdbsUQv64rUAwUoqSt9gyMMjwQJBi1NDSxoi59iDhRl8XVUVIfovFDdPYOMZbKx9lMmn0h1KUBJXek63J9bXrF4XnnTbGSzeGQkcypIMwdoakyzaH4rAFl3DvfGS5RAo0mIVJUClNSV6aSYc6Kf3CgSTe1YeqqBUia2vILZdU3PoESqSgFK6kr0+dOq6SRIVPgO1LgVS8pPlNAzKJHqUoCSuhJ9MXZFuS2ooZ7cYqUJEuMqmheqdSGEY//5yAA+pjk9RaZDAUrqSrS1Uu47UD6NmXSLraxkZl2zwsQMtaJEpkUBSupGcYp5OdNsQHVGkRhXONxROanmellXpFoUoKRuHB8aYXA4GD28ubGBRfNbyjtAFd6BGrc8EhwPHDmOu8fbsSDVXIkSItORaIAys6vNbLuZ7TSzD5fYbmb2mXD742Z2SdH2tJk9YmY/TLKeUh8KWk9Lyksxh+IkiekFqAVtTbQ0BVmAJ0bH6BuM9zxJExeKVE9iAcrM0sBngWuATcD1ZrapqNg1wMbw50bgc0Xb/xR4Kqk6Sn0pGEGi3DH4AK9ikoSZFXQxRt/PmnS/aCbf8cPTqoPIXJdkC+pSYKe773L3EeB24NqiMtcCt3ngfqDDzFYCmNka4E3AFxKso9SRrum8AwUFSQnWtniSgvFE6xAdwHYyNn9Zbtn7u6ddB5G5LMkAtRrYE/m8N1wXt8yngQ8B8caZkVnvwHQSJMZG8OGwlWOpaXfxAaxZlh+JYl93zBZUQYA6NO06iMxlSQaoUg8Qip80lyxjZm8GDrn7Q1OexOxGM9tqZlu7u/UX62w2nVHMC1tPi7DwfaTpWBMZqHZfzC4+2hfl34Ua6sNHT0y7HiJzVZIBai+wNvJ5DbA/ZpkrgLea2W6CrsFXm9nXSp3E3W91983uvrmzs7NadZcaKE6SKEc0pdvmTb97D2B1ZyRAxW1BpdLYvKX5eh3XH00ilUoyQG0BNprZBjNrAq4D7igqcwdwQ5jNdxnQ6+5d7v5X7r7G3deH+93j7u9KsK5SY/2DJzg+FGTKNTakWVxmirkfP5L/0F6dABUdaunA0ePxRzWfn/9DyfvUzSdSqcQClLuPATcDdxFk4n3L3beZ2U1mdlNY7E5gF7AT+DzwvqTqI/UtOiBruaOYA/hAPkBZ+5Kq1Km5qYHOjjYgGNW8K3aiRKQlr0QJkYpVPtxzDO5+J0EQiq67JbLswPunOMa9wL0JVE/qSMEo5pWkmBd08VUnQAGsXjqf7p5guo293X2sXTb1FB5KlBCpDo0kIXUhmoRQ9iCxAMejLajqdPFB0XOouO9CRbv41IISqZgClNSFPYfyA72u6SxvmncoakFVM0AtrSBRQs+gRKpCAUrqQjRAxelGi/JspnCg2CoGqDWdkXeh4qaaz1vK+BsUPnhM026IVEgBSmpuLJMtGKmh3ADFYA94kGFnrQuwhqaq1a041TzOoLGWbihIdS/IMBSR2BSgpOb2HzlONhsEmM6ONlqbG8va349XP4Nv3ML2ZtpagvoMj4xxrH841n7RRAmUKCFSEQUoqbnpdO9B0bxLVezeg2DQ2NUVjCihMflEpk8BSmruhYP5aTLWdlYSoCItqCqmmI+raESJBUqUEJkuBSipub3TbEEllWI+LtqC2lvRoLFqQYlUQgFKau6Fg9Xr4kuiBVWYydc3SckIvawrMm0KUFJTo2OZglHM11QSoBIYhy+qoi6+SKD040fwzFjV6yVyqlOAkpraf/g42TB1e1lHe26a9bjcvegl3eq3oJYvaiedDv6pHOkbYujE6JT7WEMT1jY+/bvDgFLNRcqlACU1Nd0MPob7IBu0TqyxFWtqrVbVchrSKVZGJlDcf7j8QWP1HEqkfApQUlPRDL51y6fZvZfA86dxlaWaK5NPZDoUoKSmqvkOVBLde+Oiz6GiWYeT0ajmItOjACU1tae7iinmVZpJt5R1kbrtPtA7ScmIBdEAdbjaVRI55SlASc2MjGY4cCSYB8qwgnTuuJIaxbzYhlUdueVn9x+LtU/hMyi1oETKpQAlNbPvcD9OkMG3fHE7TY3pso8xk8+gmhuDDMOe48Mc7Ruacp+giy8c1bzvkEY1FymTApTUzLQz+Ji5Z1BmxoaVHbnPu7p6pt6nsRnLdfM53rM/kbqJnKoUoKRmogFqXSUJEl74flESo0hEnR7p5tsVt5tv0Zrcsh/dU+0qiZzSFKCkZqbdghoZxEfD6S/SjdBcwVTxZThj1aLc8q79PbH2SS1em1tWgBIpjwKU1MxzkW6yilLMI4kHNm8pZlaNak1oQyUtqGiAOra32lUSOaUpQElNHO0b4nDvIACNDenKAlTkmY51rKxa3SayZul8GhuCRI4jfUP0DpyYcp9oF1/22N5YM/KKSEABSmpix958csPGNYtzY92VoyBALVxVlXpNJp1OsX7FwtznZ/fFaEW1LsRawpd8x06AhjwSiU0BSmpiZzRArV40ScmJFQSoRckHKCh6DtU1dYAyM2xRvpsve0zPoUTiUoCSmti+JxKg1lb2gq33dOWWZ6IFBYWZfM/FTJQwJUqIVEQBSmZcJpNl5758gDp7bfnp4T4yiA+GLZhUumDUhiRFW1CxR5QoCFBKlBCJSwFKZtyeQ32MjGYAWLKglcULyp8io6B7b8EKLF3ePFKVWrNsAQ3pIFGiu2eQ/sHyEiVcXXwisSUaoMzsajPbbmY7zezDJbabmX0m3P64mV0Srl9rZj83s6fMbJuZ/WmS9ZSZVZAgUUHrCYoz+Gamew+CuaFOW5HPOIzzPpQtWB68pwX4YA8+HG+6DpG5LrEAZWZp4LPANcAm4Hoz21RU7BpgY/hzI/C5cP0Y8Bfufg5wGfD+EvvKLLV9T370h7PWVPr8qTYBCuD0leV181kqVVBHPYcSiSfJFtSlwE533+XuI8DtwLVFZa4FbvPA/UCHma109y53fxjA3fuBp4DVCdZVZtDOvflf6mdVnCCRD1CpRTN7a5xR8MJuT6x9UnphV6RsSQao1UD0T8W9nBxkpixjZuuBi4EHSp3EzG40s61mtrW7W++Y1LuBoRH2hnNApVIpTo8MwFqOwnegkn9JN+r0aKJEnHehKEyUyCpAicSSZIAqNe5M8Wv0k5Yxs3nA/wI+6O4lpzF191vdfbO7b+7snJlMLqncM5Ff6OtXLKS5qfzkBh8+nn+Ok26EeUurVb1YTluxMDc1yKGeAQ4eG5hyHw0aK1K+JAPUXmBt5PMaoHi+gQnLmFkjQXD6urt/N8F6ygzaEXn+tLEaz58WrsRSM5uM2pBOce76/B9DTzw79WSEQYAK54bqPaC5oURiSPJf9hZgo5ltMLMm4DrgjqIydwA3hNl8lwG97t5lwaifXwSecve/T7COMsOe2Rt9/6nSALUvt5ya4QSJcReckZ/O/dGdB6csb40t+Xe1PFvwHUSktMQClLuPATcDdxEkOXzL3beZ2U1mdlNY7E5gF7AT+DzwvnD9FcC7gVeb2aPhzxuTqqvMDHdnx57CMfgqOk50BIkaBagLz1ieW35i16FYg8Ba54bcsndtT6ReIqeSRN9udPc7CYJQdN0tkWUH3l9iv19T+vmUzGJ7u/s5PhR0bc1rbWLlksrmb/Le2qWYj1u3fAEL2pvpGzjB8aERdh/oLZhxt5TUiheR3RXk+mS7fkf6/KtnoKYis5dGkpAZ89D2fMvn3PWdFc3f5O74sdoHKDMru5svteqc3HL20E49hxKZggKUzJgtT+cD1EteVGFq+HAfPhJmzTU0Q3tl3YTVEO3me/zZGM+h2hblU+KzGfzgM0lVTeSUoAAlM6J/8ATbXwgy+AzjkrMrC1B+rDBBIulZdCdz/un5FtRTzx/JjS84mYJWVNdTidRL5FShACUz4uEdB/DwFbeNaxezsL25ouNEs99q1b03rrOjjVVLgskIR8cyPP3CkSn2gNTKSIDa/7vE6iZyKlCAkhkR7d7bXGHrCSAbyX6zJadNq07VEH0OFaubb/nZkApe8vWeffhQyffPRQQFKJkBY5ksjzxzIPf5JedU1vLxzBjZgztyn6OtkVopCFC7Yryw29hMqvP03Gd184lMTAFKErftuW6GR8YAWNbRztrO+RUdxw/vhrFg/iVrXwIzNEnhZM7b0ImFb0Ts2tcTb36oSGB1BSiRCSlASeK2RtLLN79oZcWJDdmu/DOb1KpzapogMa69tYkz1wSDxzrOb7dNPUJE8XOoOC/5isxFClCSKHdny9P595am8/zJu57OLVsddO+Ne8UF63LL//HQc1OWtyXrscZgFmEf6sV7D0yxh8jcpAAliXrhUB/dPYMAtDY3cu6GyrrlfGSQ7OHxX/5GasWLqlTD6XvlhWtJp4N/Ss/uO8bzB3onLW+pFLYyX3/fvy3R+onMVgpQkqj7nsjPfXTRmctpSFd2y/nBZ8CzQDC3krVUNkxSEua3NXPpi/KJH/c8vHvKfVKr8hNEZ3bep24+kRIUoCQxI6MZ7tqyK/f58vPWTFJ6ctn9+WSC6Muu9eK1m/MDwd776POMZbKTlk+t3xzMZUWYbt79bKL1E5mNFKAkMb949PlcVltnRxsvrTC9HIoSJOqoe2/cBacvY8mC4LnS8aERHnyqeOqzQtbURnrDpbnPmR2/TLR+IrORApQkwt354W935j6/6WUbc89pyj7WwFG8L3wJNt2ILTuzGlWsqlTKuOqS9bnPsbr5zr4yt5zdvTU/S7CIAApQkpBHnjnI3u5glISWpgZeE/nlXa7oy6ypZWdiDU3TrV4iot/x0WcOcrh3cNLyqSXrSC0NuwazGbLP/CbB2onMPgpQkogf3Jcf8eF1mzfQ1tJY8bEKnj/VUXp5sWWL2nMDyDrOTx6Y+rlS6qxX5pYzO36JZyd/diUylyhASdU9f6CXx58Nhv0xjDdeVnmXnA/1kn3hkdxnq8MEiajXvSQ/jNEd9z1D15Hjk5ZPrd+MNbUD4ANHlHIuEqEAJVX33V/lX6i97NzVLFvUXvGxMk/9DLLBMEm25DRs0dpp1y9Jl5+7OjeVfSaT5as/eXzS8tbQROrMy3OfM0/9TCnnIiEFKKmqB5/az68f35P7/JbLN1Z8LB8ZJLv9F7nP6fOurovhjSZjZvyXN12UG59vy9P7eXjH5CNFpM96JYTls11Pkd29NelqiswKClBSNb0DJ/jc9x/KfX7ZuWs4e92Sio+X3f4LfHQYAFuwgtS6i6ddx5lw5prFXHVJfiqQL9/52KTvRdmCZaTPekXuc+aBb+JDk49GITIXKEBJVbg7n/veQ/QNBO89LZrfyv/+1soDio+NBN17ofR5b6j71lPUu15/fi4xZP+Rfr7/6x2Tlk+/+O3BCO2Ajwwwdv831NUnc54ClFTFzx95vmBQ2Pe/7cXMb6ts1lyA7LP35d4LsrZFpCIvtc4GC9ubeedV+eGMvvkf2/jNk3snLG+NLTRc/u7c5+yeR8k+tyXROorUOwUombat27v4/A/zmXZvuPQMLt64ouLj+cggmSfvzn1On/s6LN0wrTrWwtUvPYMzVuen4vj0tx8smHqkWGrlOeHzqEDmwW+SPbw76WqK1C0FKJmWu7fs4pNfu4+R0QwAK5fM44Y3nF/x8XxshLF7PosPHAHAmtpJnfnyqtR1pjWkU/zXd7+c1UuDCRqz2Sz/3zfv54lJZt4t7OobZPTuvye778kZqa9IvVGAkopkMlm+8dMn+Z93PIwTPCvp7Gjjr951BS1NlbV2PDPG2C8/T/ZQfoik9KXXYY2VdxXW2sL2Zj723leyrCNItR/LZPib237Nt37+O0bHMieVt8YWGl51I9bUFqwYO8HoPZ8lo1EmZA6yU+lB7ObNm33rVqXoJimbdX79xB7+7Z7fceBo/iXU01ct4q/fdQWL5rdUdFzPjDL2238lu+uB3LqGze8gvel1065zPThw9Dgf+cIvONY/lFu3eul8bnzrJZxXYo4s7z3A6H98JteShGCKjvQFbya17IwZqbPITDGzh9x980nrFaBkKu7Onu5+Htrexb2PPJ8bY2/cJWet4C/eeVlFLSfPjJLd+RsyT/wEHzyWW58+72oaLnnbtOteT/Yf7ufT33mQZ/cdK1i/bvlCXnHBWl554TqWLmzLrffBHkZ/9s/4sT0F5VOrNpE64/Jg2vvm+pkXS6RSNQlQZnY18I9AGviCu3+yaLuF298IDALvcfeH4+xbigLU9IyMZug5PkzfwAkOHBtgX3c/+w73s/2FIyUHPm1raeRtrziba684K/ZI5T42gh8/jHc/hx9+juy+JwsCE0B648tJX/auWZVWHlc26/zkwWf5+k+fZHhk7KTtKxbPY8PKDk5f1cHKJfNY2mZ0PvsjWvY9xMmXw7Al60gtPR1bsCz4mbcUmudBUxuWUg++zA4zHqDMLA3sAF4H7AW2ANe7++8iZd4I/AlBgHop8I/u/tI4+5ZSaYDa9+wOtv/wK2XvV0te4kN03fj/Vw+XnWBCWncn6042myXjwbOksUyW0UyWzBST7I1rSKd40WlLOWfdYpoa0+MnCU+czf9kxyAzhmfGYGQQP3Ecxk5MeFxrmU/6/GtInX3VKf/L9XDvIN/4j23c9+Teks+iii3yXi63JznXdtOQCkasMIOUWbAMYORGsMCM0VQzGWskm2ogYw24pcmSwi3/AxbeN4YHB8jx3IdT7w8Fqa4VL3kDmy5/VcX7TxSgkszdvRTY6e67wgrcDlwLRIPMtcBtHvw2vd/MOsxsJbA+xr5VM9BzjOaDypSaTDplLGhvpmNeMwvntdCQ3gdd+6jG2NvWMp/0ua8nddarZnVCRDmWLmzjA29/CX/85ot44Hf7+cVjz7Nt9+EJ/0g4Zgv5EVdwX/Y8zsk8x2l0sdK7pzjL8ElrDIUbqb7+7sozdyeTZIBaDUQ7z/cStJKmKrM65r4AmNmNwI0A69atm16N5zgzaEynaGhI0dSQpqWpgZamBlqbG2hraaxOl1sqjbUswBavI9W5AVu6Hus8o27neEpaa3MjV158GldefBqjYxleONjHrv3HeO5AL4d7BjnSN8TR/iH6B0ZwnGO2kPvsIu7jIpp9hDUcZDG9dHgfi+hnHoO0Mkyzj9b6q4lMW5IBqtRvs+L+xInKxNk3WOl+K3ArBF185VRw3PL1ZzD0uvdVsmviYocEK/hP0O2TCj6lMFJht1DKjHQq+GlIp4Ng1JimOfxviQcd8WpiFnQvmYGlgp9UGtKN0NCINbRAyzxobD0lny1VQ2NDmjNWL8q93Bvl7pwYzTA4PMrQyBhjY1lGxzKMZbJksh7+ZMlmw6T/bAZGhrDsKDY2AtkRLJvBPYNlM5g7kA3/VWXDz/mu4fEuWyv9z06kwPLTz07kuEkGqL1AdG6ENcD+mGWaYuxbNQuXLOXCq65O6vAi02ZmuRatyFyR5JPoLcBGM9tgZk3AdcAdRWXuAG6wwGVAr7t3xdxXREROYYn9OebuY2Z2M3AXQar4l9x9m5ndFG6/BbiTIINvJ0Ga+Xsn2zepuoqISP3Ri7oiIlJTE6WZn9ovm4iIyKylACUiInVJAUpEROqSApSIiNSlUypJwsy6geencYilwOEqVSdJqmf1zIY6gupZbbOhnrOhjlCdep7m7ifNO3NKBajpMrOtpTJJ6o3qWT2zoY6gelbbbKjnbKgjJFtPdfGJiEhdUoASEZG6pABV6NZaVyAm1bN6ZkMdQfWsttlQz9lQR0iwnnoGJSIidUktKBERqUsKUCIiUpfmRIAys7Vm9nMze8rMtpnZn5YoY2b2GTPbaWaPm9klkW1Xm9n2cNuHa1zPPwzr97iZ3WdmF0a27TazJ8zsUTNLZNTcmHW80sx6w3o8amYfjWyrp2v5f0bq+KSZZcxscbgt8WsZnqfFzB40s8fCen6iRJma3psx61jT+7KMetbDvRmnnjW/N8Nzpc3sETP7YYltyd+X7n7K/wArgUvC5fnADmBTUZk3Aj8mmDr2MuCBcH0aeBY4nWAixceK953hel4OLAqXrxmvZ/h5N7C0Dq7llcAPS+xbV9eyqPxbgHtm8lqG5zFgXrjcCDwAXFZP92bMOtb0viyjnvVwb05Zz3q4N8Nz/TnwjQmuWeL35ZxoQbl7l7s/HC73A08Bq4uKXQvc5oH7gQ4zWwlcCux0913uPgLcHpatST3d/T53PxZ+vJ9gtuEZE/NaTqSurmWR64FvJlGXyYT32/HwY2P4U5y5VNN7M04da31fhnWIcy0nMpP3Zrn1rMm9aWZrgDcBX5igSOL35ZwIUFFmth64mOCvlqjVwJ7I573huonWJ2qSekb9Z4K/YMY5cLeZPWRmNyZYPWDKOr4s7ML4sZmdG66ry2tpZm3A1cD/iqyesWsZdqM8ChwCfurudXdvxqhjVM3uy5j1rPm9Gfd61vje/DTwISA7wfbE78vEZtStR2Y2j+B/9Afdva94c4ldfJL1iZminuNlriL4RfDyyOor3H2/mS0DfmpmT7v7L2tQx4cJxtY6bmZvBL4HbKROryVBF8pv3P1oZN2MXUt3zwAXmVkH8O9mdp67Pxn9GqV2m2R91cWoI1D7+zJGPevi3ox7PanRvWlmbwYOuftDZnblRMVKrKvqfTlnWlBm1kjwi+rr7v7dEkX2Amsjn9cA+ydZX6t6YmYXEDS7r3X3I+Pr3X1/+N9DwL8TNLVnvI7u3jfeheHudwKNZraUOryWoeso6kKZqWtZdM4e4F6Cv5ij6uLehEnrWPP7Mk496+XenKqeEbW6N68A3mpmuwm66F5tZl8rKpP8fRn3YdVs/iGI6LcBn56kzJsofOD3YLi+AdgFbCD/wO/cGtZzHbATuLxofTswP7J8H3B1jeq4gvxL4JcCL4T71dW1DMstBI4C7TN9LcPjdwId4XIr8CvgzfV0b8asY03vyzLqWQ/35pT1rId7M3LOKymdJJH4fTlXuviuAN4NPBH2+wL8NcE/Ktz9FuBOgqyUncAg8N5w25iZ3QzcRZCd8iV331bDen4UWAL8i5kBjHkwkvBygq4CCG6Qb7j7T2pUx3cA/4eZjQFDwHUe3Ln1di0B3gbc7e4DkX1n6lpCkG34VTNLE/RofMvdf2hmN0XqWet7M04da31fxq1nPdybceoJtb83TzLT96WGOhIRkbo0Z55BiYjI7KIAJSIidUkBSkRE6pIClIiI1CUFKBERqUsKUHJKCUd9Hh8B+tvhUDHVPP69ZrZ5ijIfjJ7XzO4MRwyY1cxsZalRrSs81vlm9pVqHEtOXQpQcqoZcveL3P08YAS4qQZ1+CCQC1Du/kYPRgyY7f4c+Hw1DuTuTwBrzGxdNY4npyYFKDmV/Qo408wWm9n3wjlr7g+H5MHMPm5m/2pm95jZM2b2x+H6K6MtBTP7ZzN7T/HBzexzZrbVInP6mNkHgFXAz83s5+G63eFwOpjZn4etuyfN7IPhuvUWzFv1+fBYd5tZa4nzfSU858/NbJeZvcrMvhTu+5VIudeb2W/N7OGwFTkvXP9RM9sSnvtWC9/2DFuF/92COYp2mNkrJriebwd+Eu7znvCa/sDMnjOzm8Pv9kh4jRdHjr05XF4aDp0z7gcEQ/mIlKQAJackM2sgmJfoCeATwCPufgHBaBK3RYpeQDBky8uAj5rZqjJO81/D0RIuAF5lZhe4+2cIxh27yt2vKqrTiwnetn8pwdAwf2xmF4ebNwKfdfdzgR6CYFDKIuDVwJ8R/IL/B+Bc4HwzuygMhB8BXuvulwBbCVo+AP/s7i8JW5etwJsjx21w90sJWn8fKz6pmW0Ajrn7icjq84A/IBgy6O+AQXe/GPgtcMME9Y/aCkwUDEXmzFBHMne0RoY2+hXwRYJpNt4O4O73mNkSM1sYlvm+uw8BQ2GL51KCABHH/2bBdAcNBMPXbAIen6T8y4F/Hx+6xsy+S/AL+g7gOXcfr/dDwPoJjvEDd3czewI4GHaVYWbbwn3WhPX4TdhAaiIIGABXmdmHCLofFwPbCIIcwPhguhOdeyXQXbTu5x7MtdVvZr2RYz1BELSncoigtSlSkgKUnGqG3P2i6IrxrqwiXvTf6PoxCnsXWop3DlsUfwm8xN2PhV1sJ5Ur3m2SbdGWSYaghTNZuWzRPlmCf88ZgvmFri+qbwvwL8Bmd99jZh8vqu/4sTKU/r0wxMnfr/j80bqNHyN6LYv3bwmPK1KSuvhkLvgl8IcQPF8CDnt+bqhrzazFzJYQjNq8BXge2GRmzWFL6zUljrkAGAB6zWw5QXfiuH6CaeZL1eP3zKzNzNoJBgP91TS/W7H7gSvM7EwIJrwzs7PIB4fD4TOpd5R53B1M3KqbzG7gxeFy8TnPAkrNgSQCqAUlc8PHgS+b2eMEoy7/UWTbg8CPCEY5/xsP59oxs28RdNc9AzxSfEB3f8zMHiHoJtsF/Cay+Vbgx2bWFX0O5e4Phy2tB8NVX3D3RyyY8bcq3L07TOj4ppk1h6s/4u47zOzzBN1vuwkCcTnHHTCzZ83sTHffWcaunwK+ZWbvBu4p2nYVwbUXKUmjmcucFXZzHXf3T9W6LrOBmb0NeLG7f6QKx2oGfgG83N3Hpl05OSWpBSUisbj7v4ddodWwDviwgpNMRi0oERGpS0qSEBGRuqQAJSIidUkBSkRE6pIClIiI1CUFKBERqUv/P8kyhtaR+4V2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pmf_mean = marginal(posterior, 0)\n",
    "pmf_mean.plot()\n",
    "\n",
    "pmf_mean2 = marginal(posterior2, 0)\n",
    "pmf_mean2.plot()\n",
    "\n",
    "decorate(xlabel='Population mean (mu)', \n",
    "         ylabel='PDF', \n",
    "         title='Posterior distributions of mu')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:37:29.283479Z",
     "iopub.status.busy": "2021-04-16T19:37:29.263758Z",
     "iopub.status.idle": "2021-04-16T19:37:29.391453Z",
     "shell.execute_reply": "2021-04-16T19:37:29.390972Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5V0lEQVR4nO3de3xcd3nv+88zo4slW1dLtmXJji+xndhxEhLnAgl3UpK0u+nZpy1waLPb00LZlO7SfWV3n9Lu9vXa3ZeevSktGwqUlpQeaGmhDRCuhQCBBOLc7Pga3++2bEuWLVmWZuY5f6zRzJrRSLNGnpFGo+/79dIra81aa+anFUuPfr/fs56fuTsiIiLVJjbXDRARESlEAUpERKqSApSIiFQlBSgREalKClAiIlKV6ua6AeXU1dXla9asmetmiIhICZ577rnz7t6d/3pNBag1a9awffv2uW6GiIiUwMyOFnpdQ3wiIlKVFKBERKQqKUCJiEhVUoASEZGqpAAlIiJVSQFKRESqkgKUiIhUJQUokVni7qRO7yF17sBcN0VkXqipB3VFqpVfuUDih4+ROrMXgPo3vpfYqtvmuFUi1U0BSqSC3J3UvidJPP9FSFzLvJ544R+o77sVM5vD1olUNwUokQrxVIrEd/43qZM7Jx8bPEXq2AvEb7hjDlomMj9oDkqkQlInXsoJTtbWQ+yGOzP7yR1P4O5z0TSReUE9KJEK8ZMvZ7Zja+6i7r5/AeOjjJ3YAclxfOA4fnIn1nfrHLZSpHqpByVSAe5O6tTuzH580xuweD22qIX4xtdlXk/u+Ip6USJTUIASqQC/dAYfvgiA1S/CutZkjsW3/ATEgsGL1Pkj+Ok9c9FEkaqnACVSAX5qV2bbem7C4tnRdGtuJ77h/sx+8qUvqxclUoAClEgFhIf3Yiu3TDoev+WtEIsH5/YfxAdOzFrbROYLBSiRMvPEGKmz+zP7hQKULe4ktur27DWh80UkoAAlUmZ+7gAkxwGw1uXYkqUFz4st35DZTp19ZVbaJjKfKECJlFkqNP8UW7l5yvMsFKD83AHNQ4nkUYASKbPUyXCAmjy8N8Hae7H6JgB89DIMna1420TmEwUokTLy4Yv4pdPBTiyOLd845blmhi27MbOfOnew0s0TmVcUoETKKHUq+0xTbPlGrL5x2vNjoQDl5zQPJRJW0QBlZg+a2T4zO2BmHyhw3Mzsw+njO8zsjtCx3zKzXWb2spl91swWVbKtIuUQdf5pgi0P9aCUKCGSo2IBysziwEeAh4DNwDvMLP8n9iFgQ/rr3cBH09f2Av8K2ObutwBx4O2VaqtIObh7TlUIixKglq6BeH1w/ZXz+MhghVonMv9Usgd1N3DA3Q+5+xjwOeCRvHMeAR7zwDNAu5n1pI/VAU1mVgc0A6cq2FaR6zd8ER8bAcAaFmPtvUUvsXgdsVAZJPWiRLIqGaB6geOh/RPp14qe4+4ngT8CjgGngUvu/o1CH2Jm7zaz7Wa2vb+/v2yNFylVJjkCsI6VkRcjDCdSuJaDF8moZIAq9NOZ/6BHwXPMrIOgd7UWWAksNrNfKPQh7v5xd9/m7tu6u7uvq8Ei1yMnQLWuiHxdOFFCPSiRrEoGqBPAqtB+H5OH6aY65y3AYXfvd/dx4AvAayrYVpHr5pfOZLatLXqAsu71YMGPog+eygwTiix0lQxQzwIbzGytmTUQJDk8nnfO48Cj6Wy+ewmG8k4TDO3da2bNFoyTvBnQmgRS1XIDVM80Z+ay+kasc+LvNNcwn0haxQKUuyeA9wFfJwguf+vuu8zsPWb2nvRpTwCHgAPAJ4D3pq/9EfB3wPPAznQ7P16ptoqUgw+GhvhKCFCQN8ynACUCVHjJd3d/giAIhV/7WGjbgV+f4trfBX63ku0TKRcfvYKPDQc78QZY3FHS9bZ8A+z5p+C9zipAiYAqSYiURU6CRNuKyBl8E3J6UBeO4MlE2domMl8pQImUQX6AKpUtasGa072uVBK/fK5cTROZtxSgRMpgphl8YdbRl32/gZPX3SaR+U4BSqQMZprBF2YdK7PvpyXgRSqbJCGyUJQjQF1rXsGVy6M4zrUDexlsPUtnaxN93S0lz2mJ1AIFKJHr5IkxfPhCsGMxrKX0iiaJZIr/8dVjPNA/AMDlUzv4xM7vA/CLP7GVn3ntprK1V2S+0BCfyHXyoVDvqaUbi5f+d983tx9m1wUjlf6RbPFhGn0MgM99ezfnBobL01iReUQBSuQ6+WAoQJVQg2/C6FiCz39nDymLc8HaWNJUT9viRtYtvgrAeCLJZ775ctnaKzJfKECJXKecHlR76QHqy0+/wqXhUQBGGpexcdVSNqzq5Ffu68qc84Odx9l79Pz1N1ZkHlGAErlOM61iDnB55Br/+NT+zP5Nt24lFgsSIlY1XOHVW7Kp55/66ksExVdEFgYFKJHrlFODr720DL5/fGo/I6PjAKxc2sLWO27Pvu/ASR5961bq4nEADp4c4MkXjl5/g0XmCQUokevgeVUfSulBXRy6ypefztbde/tbtlDXtTr73oMn6W5v5pH7NmRe++tv7WJ0TGWQZGFQgBK5HlfOQyoJgDW1YQ1NkS/94vf3MZ4Irl2zop3XbOmFpjasYTEAPj4Kwxf556+/iY6W4H0HLl/lW9sPl/mbEKlOClAi12GmD+i6Oz98OVst4p0P3IKZBV8dvdnzBk6wqKGOn339TZnXnt6lMkiyMChAiVyHmQaoE/2XGbwSZO4tXtTA7Tcuz75PKECl0jX5Xn1LH0aQPLHv2AUGLo9eV7tF5gMFKJHrkFvFfPk0Z+bacTA7b7V1XXcmcw/yisYOBgGqbXEjN68J0s4d59m9p2bcZpH5QgFK5Dr40NnMdilVzHceCgeoZTnHrD1cNDY7nHfv5mzP6hkN88kCoAAlch38yoXMtrUsm+bMrGQyxcuH+zP7W9fnB6heSA/n+dBZPBmkod+zORu4dh7u58rVsZk2W2ReUIASmSFPjuNXL6X3DJraIl138NQAV68FQWdpaxMrly7JOW71jVhLuoqEpzLDiF1tzazvDRY1TKVSPLfvNCK1TAFKZKaGBzKb1tweuUjsjvDw3vplBZfSyM3k0zCfLEwKUCIzlFliA7AlSyNft/Ngdnjv1nWFhwWDYb7054QC1D2hAPXCgbN6aFdqmgKUyAz5lYvZncWdka65NpZgz7Fs0df8BIkJhTL5AHq7Wli1LBhKHE8keX7/mUnXitQKBSiRGfLhbICK2oPae+wCyWQKgL7uVjpbC1eeCD9T5ZfO5hwLJ0s8s1vDfFK7FKBEZio8xBexBzVdenmYtXSDBT+ePnwBH7+WOfbq0DDfc/tOZ8olidQaBSiRGcrpQUUMUOEEiVvXTxOg4nU5S8f75Wwv6oYVbSzvCOr1jY4l2HvswqTrRWqBApTIDIWfgWJx8SG+yyPXOHRyEADD2LKma9rzw5XRw8N8ZsZtodJI4WeqRGqJApTIDLg7PhJKM4/Qg3r5cD9OsODg+t4OFjc1THt+eHXecEklgC1rsr2r3Ue00q7UJgUokZm4OphdZqNxCVbfWPSSvUezPa7phvcm5PagcrP1Nod6X/uPX2RsXPNQUnsUoERmIGd4L2IG35EzlzLbN6YrQkwnXNsvvwfV2dpET7oCRSKZ5JUTFxGpNQpQIjPgw6UN77k7R84MZvbX9LQXvSYnQA2dxVOpnOObw8N8RzXMJ7VHAUpkJq6UlmJ+YehqprhrU2M9y9qbi15jDc3YRH2/VDJYvTcknGSx+4gSJaT2KECJzECpZY7Cw3s3LG8rWH+vkNxeVP48VLYHtffYBRLJ3B6WyHynACUyAznLbERIMc8Z3lsRreo5TJ8o0d3eTHe6JzY2nuTgyQFEaokClMgMhB/SjVKH7+iZocx2SQGqbeoABbBlbbYXtUvDfFJjFKBESuTueXX4ogSowcx2lASJzHsXC1B6HkpqmAKUSKnGhiGRro1X1wgNi6c/fTzJqfNXgKCCxOplrZE/Kj/V3N1zjoefh9pz9HymEK1ILVCAEilRfg2+YgkPx85eylSQ6Fm6hMaGaAsbAtDcEQRBwMdG4NqVnMPLOxazNF0RfXQsweHTg9HfW6TKKUCJlOpKaUViczL4Sph/gqDuXk4vavD0pOPhbL5dGuaTGqIAJVKi60kxX9NTWoCC6VPNIXeYT89DSS1RgBIp0fWlmLeX/HnTpZpDbibf7qPnJ81TicxXFQ1QZvagme0zswNm9oECx83MPpw+vsPM7ggdazezvzOzvWa2x8xeXcm2ikSVk2JepAfl7hw9G+pBlTjEB8Uz+VYuXULr4mCeamR0nFMXrkw6R2Q+qliAMrM48BHgIWAz8A4z25x32kPAhvTXu4GPho79MfA1d78JuA3YU6m2ipQiN0li+qKv/YMjjIyOA7CkqSGT0FCK3OXfJwcoM2NDb3YubL8WMJQaUcke1N3AAXc/5O5jwOeAR/LOeQR4zAPPAO1m1mNmrcDrgD8HcPcxdx+sYFtFoiuhDt9MSxyFWUtX7vLvibFJ52xYlW3HK6ooITWikgGqFzge2j+Rfi3KOeuAfuAvzOwFM/ukmRV82MTM3m1m281se3+/Joilsnz8Gj6R6m0xaGqf9vzcCualD+8BWLweW5JNhCiUKLGxL9SDOq4elNSGSgaoQn8q5s/eTnVOHXAH8FF3fxUwDEyawwJw94+7+zZ339bd3V3oFJHyyVlFtwOLTf8jdDScwTeDBInMZ+UM852ddPzG3g4s/eN09Mwlro0lZvxZItWikgHqBLAqtN8HnIp4zgnghLv/KP363xEELJE5VXoG3/UlSGQ+a5rFCwEWNzXQ290CQMqdQ3pgV2pAJQPUs8AGM1trZg3A24HH8855HHg0nc13L3DJ3U+7+xnguJltSp/3ZmB3BdsqEkkpz0CNjiU4e3EYgJgZfd3RSxzls7bl2TYUSJQAuLEvm7Cx/7hW2JX5r4SaK6Vx94SZvQ/4OhAHPuXuu8zsPenjHwOeAB4GDgAjwC+H3uI3gL9OB7dDecdE5kQpVcyPnsmWOOrtaqGhPj7jz815Fmpo8hAfwKZVS3nyhaMA7NcS8FIDKhagANz9CYIgFH7tY6FtB359imtfBLZVsn0iJSthqfcT/dklNlYtn/nwHuT1oIbO4u6TMgI3hBIlXlGihNQAVZIQKYGHkySa26c9d6KCOZCZH5opa1yCLUq/R3Ichif3kFYva8300i4MXeXi0NXr+kyRuaYAJVICHxnM7jRP/5DuqfOXM9u9S5dc92dba3geanKiRDwe48bwA7sa5pN5TgFKJKJgocLcNPPphEsOrey6vh4U5BeNLTwPtaE326ZXlCgh85wClEhU41chma7iUNcI9VOXLUomU5y+GA5Q5e5BFc7kC1eUUA9K5jsFKJGIcnpPze3Tli06NziSWd22o6WJpsb66/78YkVjATauyqa+Hzw5oBV2ZV5TgBKJKidBYvrhvZPh+acyDO9BtFTzpa1NdLYEPbtr4wmOnxsqeJ7IfKAAJRJROEGi6PxTKED1lCFBAoAlXRALsvT86qVgCfgCNmqYT2qEApRIROEhPmYxxXyCxWJYy7Jse6ZKlMgpHKsAJfOXApRIVCU9A5XtQa0sVw+K/Hmo4gHqgJbekHlMAUokIi9hDqrcKeaZz81JNS+cKLE+VNn8xLkhRlXZXOYpBSiRiHKG+KaZgxoZHWfgclDFIR6Psbyj4FJmM5KTKDFFJt+ihjr6lgWFaR3noHpRMk8pQIlEFU6SmKYHlZMg0bmEWKz0VXSnEiXVHIL1oSZomE/mKwUokQh8fBQfT9e2i8Whcep5pdPh4b0yzj9B3sO6l8/hqWTB8xSgpBYoQIlEkdd7mu4h3RPhZ6DKlMGX+eyGJqwpXRk9lYQrhauW5wQopZrLPKUAJRKBj5RQg+98ZRIkMp8fIVFiTU878Xjw431ucJih4Wtlb4dIpSlAiUSQ+wxU9Id0yz3EB9ESJeriMdauaM/sa5hP5iMFKJEIPGKChLtz6kIoQFW6BxU5UULDfDL/KECJRBHxId0LQ1cZGw8SF5Y0NdC6uLHsTYmy7Abkz0OpByXzjwKUSARRn4Gq9PwT5GXyTReg8ipKuHtF2iNSKQpQIhFEHeLLmX8qwxpQBS3uhHiwfIePXsavXSl4Wl93C4sa6gC4NDzKBS0BL/OMApRIFBGH+CqxzEY+M4uUKGFmrF8ZWmFX6eYyzyhAiRThibFsL8VisKh1ynMrssxGATNJlFDJI5lvFKBEigkP7zW1YbGpf2xyelDdUwey62XtPZltHzw15XnheahXlCgh84wClEgRURcqHBtPcn4wmOcxjJ7O8hWJzWdtoQB16fSU5+X0oE4pUULmFwUokSLCVSSme0j37MAwThAAutqbqK+LV6xNsfaV2fYNTh2gutubM6nuV6+N5/TwRKqdApRIMTkZfO1TnhYuElvJ+Scgd/n3kQF8rHCGnplpHkrmLQUokSKi1uE7czEcoCqTwZdpR7wub/n36RIlNA8l85MClEgRUYf4zlwczmyvqOD80wSLOMy3ISdRQqnmMn8oQIkUEa4iMf0Q3+ykmGfaEjFRYkNfNqgePjPIeKLwGlIi1WbaAGVm3wht/8fKN0ekCuXMQXVOeVp4DmpF5ywEqJxU86kDVEtzY6Y9yWSKI2cuVbxtIuVQrAfVHdr+uUo2RKQaeTKBXx1K7xk0FX62aTyRm2K+vGMWhvhyelBTPwsFcGOoF7X/uIb5ZH4oFqD00IQsbFeHmPgxsKZWLF5X8LT8FPOG+sqlmE+w1mVBZQvAr1zEx6delHBDb7hwrAKUzA+Ff9qy1pnZ44CFtjPc/acr1jKRKuAjoV/m08w/5SZIVH54D8Di9VhLd7qiueNDZ7Glqwueu3GVEiVk/ikWoB4Jbf9RJRsiUo1yEySmzuCb7fmnCdbek1lywwdPwRQBas2KYAn4ZDLF6QtXuDxyjZbm8q9VJVJO0wYod//uxLaZdadf6690o0SqRsRnoGb1Id2QYB7qRWD6TL6G+jhrVrRlHtQ9eGqQ229cPuX5ItWgWBafmdnvmtl5YC+w38z6zeyDs9M8kbnlw9nhsOl6UGfmLEBFexYKch/Y3X/8QsXaJFIuxZIk3g/cD9zl7kvdvQO4B7jPzH6r0o0TmWs5D+kunjrFPFxFYsVsBqj28LIb0weojeEVdlVRQuaBYgHqUeAd7n544gV3PwT8QvqYSE3LmYOaYogvkUxxbmAks79iFlLMJwQLFxoAfrkfT4xNeW441fyVkxdV2VyqXrEAVe/u5/NfTM9D1VemSSJVJGcl3cI9qHCK+dLW2Ukxz7SprgFr6UrvOX753JTn9na10Lwo+LEdGr5G/+DIlOeKVINiAWrqP8emPyYy7wUP6U6ULzJobit4Xu78U2WLxBaS88DuNPNQ+ZXN9yvdXKpcsQB1m5kNmdnl9NfQxD6wtdibm9mDZrbPzA6Y2QcKHDcz+3D6+A4zuyPveNzMXjCzL5f2bYmUwcggmYd0m9uwWOGe0WwXic0XdXVdyHtgVwFKqlyxNPMZj1WYWRz4CPAAcAJ41swed/fdodMeAjakv+4BPpr+74TfBPYAlVs7W2QKUTP45irFfELUorGQuwS8elBS7YqlmS8ys/eb2Z+a2bvNrNiDvWF3Awfc/ZC7jwGfI/fBX9L7j3ngGaDdzHrSn90H/CTwyRI+U6Rsqj2Db0LUIT7IXXrj0KlBEslUxdolcr2KDfF9GtgG7AQeBv7fEt67Fzge2j+Rfi3qOR8C/j2gnyCZG5GX2Qj1oGaxisSEnAB1+dy0mXwdLYvoamsGggK3x86qsrlUr2IBarO7/4K7/xnws8BrS3hvK/Bafl5rwXPM7KeAc+7+XNEPCXp2281se3+/ilxI+eQM8U3RgwpSzOd4Dqq+EWtJLzzgKfzS1KvrQm4vap8qm0sVKxagxic23D1R4nufAFaF9vuA/Bncqc65D/hpMztCMDT4JjP7TKEPcfePu/s2d9/W3d1d6BSRGckd4is8B9U/OEIq/TxRZ0sTjQ2ljIKXj3X0ZbZ94Pg0Z8LNN3Rltvcem/QUiUjViJrFN5G5d2soq2+oyLXPAhvMbK2ZNQBvBx7PO+dx4NF0Nt+9wCV3P+3u/9Hd+9x9Tfq6b7v7L5T+7YnMXJRCsTlFYudg/mlCboA6Me25m0KVzfcdU8kjqV4Vy+Jz94SZvQ/4OhAHPuXuu8zsPenjHwOeIJjbOgCMAL88088TKbsIQ3zhBIm5yOCbEOtcxcRC7n5x+gC1pqed+ro444kk/YMjXBi6ytLWpso3UqREFR2PcPcnCIJQ+LWPhbYd+PUi7/Ek8GQFmicyJU+M4WPpuSWLwaLCTzrkLrMx+/NPE/J7UO6OWaEpXqiLx9jQ18nuI8Gc7b5jF3jNLX0FzxWZS8WG+EQWprxnoCxW+EelWnpQLO7E6oNekI+N5JRoKuTm1Usz2xrmk2qlACVSgI8MZrYjrwM1BynmE8wM68g+xVFsHmpjKEDtVYCSKqUAJVJAOMWcKRIkEskUZ0Nljua0BwVYRzYhNlVkHiqcKHHo9CBj48lpzhaZGwpQIgV4zkq6U1cxn0gxX9o6dynmE6wzeiZfS3Mjfd3BvFoqleIVlT2SKqQAJVJIhHWgTs9xFfN8paSaA2wMp5trhV2pQgpQIgVEqSIx10Vi81n7SjKLFw5NX/II4CYlSkiVU4ASKSCnisQUc1Cnzl/ObK/sqoIAVdeAtS5P7zk+eHLa8zflJUpohV2pNgpQIoWUOMS3smvuh/iA3Ey+IokSvV0tLGlqAODK1TFOhb4fkWqgACWSx8dH8fGrwU6sDhoL947CPahqGOKD0uahzIxNqzTMJ9VLAUokj+f1ngpVZLg2luDCUBDEYmYs75i7KhJhsRIy+WDyMJ9INVGAEskXYSXd8DLvyzoWUxevjh+lQiWPpqNECalm1fFTJVJFoiyzcepCdSVIZDR3YA1Bb87HR+HK9EHnxt4O4ungeqJ/iEvD1yreRJGoFKBE8kRJMT91vrpSzCeYWUkP7DY21HFjbzYIv3xYi35K9VCAEskXrsM3VYp5uAdVBQ/phoWH+VIR5qG2rl2W2d6lACVVRAFKJE9uDypKFYnq6UFB3jzUxelX1wW4ZV12JeqXD52rSJtEZkIBSiRPOIuPCFUkquUZqAm2dHVm2y8cKXr+plVLM/NQJ89f5mI6O1FkrilAiYS4e26h2AJDfFeujjGUTiaor4vT1VZdq9Fa20qoawSCZUO8yNpQDfXxnOehdh85X9H2iUSlACUSNjYMiXQmW10jNDRPOiV3Fd0lU65cO1csFiMW6kWl+g8XvWbL2tAwn+ahpEooQImE+OVs78GWdBUMPuEA1VtNKeYh1rUusx1lmG9rToDSPJRUBwUokRC/nO09WEt3wXNOVmGJo3zWtSaz7RF6UBtXdVJfFweCAHxB81BSBRSgRMKuhHpQLV0FT6nmDL4JsVCASl04iqdS055fX5c7D6V0c6kGClAiITk9qCWFA1TuMhvVlcE3wRZ3Yk1twU7iGn7pVNFrtqzNfr87lW4uVUABSiTEc3pQk4f43H1e9KAArGttZtvPHyl6/tZ1emBXqosClEhIfpJEvsEr1xgdSwDQ1FhP2+LGWWtbqWI5Aar4PNSGvk4a6oN5qLMDw/QPjlSsbSJRKECJpHkyEaoiYQUf0s1fRbfaUszDwokSUVLN6+KxnOrmqiohc00BSmTC8EUgWJ7CmtuxuoZJp8yX4T2YCFBBAPXBU/h48Urlt4Tq8u3UMJ/MMQUokTSPkMGXm2JenQkSE6x+Eda2Ir3n+IWjRa/ZGqrL99KBs0XXkxKpJAUokbQoGXzHzw1ltlcta614m65XrDv0wG6EeagbeztpaQ7m1QavjHL49GClmiZSlAKUSFq4B0WNBKiceagIASoWM161YXlmf/u+05VolkgkClAiacWqSIyMjnP+UpDZFo/HWFnlc1BQeqo5wJ0bezLbz+8/U+4miUSmACWSlpNiXiBAnejP9p5Wdi6hLl79Pz7WvhLiQbKHjwzgocUYp3LbjcuwdHLFgRMDmcrtIrOt+n/CRGaBu8OV6XtQ4eG9vnkwvAdgsXhuZfMIw3wtzY1sXB2k2DvOiwfOVqx9ItNRgBIBGBvGx0eD7bpGaJw8fHfsbDZArV4+PwIUgIUTJc4diHTNHRtXZLY1DyVzRQFKhGjLbBw7dymzPR8SJCbElm/IbKdO74t0TXge6sVXzpJKKd1cZp8ClAhRU8yzz0CtXt5W8TaViy3bABb8qPvAcfzalSJXwJoVbXS0BCsFD4+Osf/4hYq2UaQQBSgRKLrMxpWrYwxcDtZIqovH6ems/gy+CdbQlLv8xpn9xa+x3HRzZfPJXFCAEqF4ink4QaK3u4VYrHpr8BViKzZltv303kjXhOehnlOAkjmgACVCXpmjAkN8x85m559Wz6P5pwmxFTdltlNnogWo29YvJ55OpT9yZpCLWmVXZpkClAjFn4EKzz/NpwSJCda9DmJ1APjQWXxkoOg1zYvquXl1Nlg/p2w+mWUKULLg5S6zQcFlNnJ6UPMoxXyC1TUQW7Y+sx81m2/bTdlsvh/uOlH2dolMRwFKJMIyG8f751cNvkIsNMznEYf5Xr2lN7P98qF+LqmqhMyiigYoM3vQzPaZ2QEz+0CB42ZmH04f32Fmd6RfX2Vm3zGzPWa2y8x+s5LtlIWt2DLvl4avZcr9NNTHWd6xeNbaVk6xUKJE6sy+SEtpdLU1sym9iGHKnadfVi9KZk/FApSZxYGPAA8Bm4F3mNnmvNMeAjakv94NfDT9egL4N+5+M3Av8OsFrhUpi2LPQIWH91Yta63qVXSnY0tvCKpkQDCkeTnagoT33bIqs/0DBSiZRZXsQd0NHHD3Q+4+BnwOeCTvnEeAxzzwDNBuZj3uftrdnwdw98vAHqAXkQootszGfFtiYyoWryO2fGNmP3Um2jzUq7f0ZorH7jlyXtl8MmsqGaB6geOh/RNMDjJFzzGzNcCrgB8V+hAze7eZbTez7f39WqJaSlfsGaicGnzL5k8FiUJiPeFhvmjzUJ2tTWxeEwRux3l618mKtE0kXyUDVKFxkPxB72nPMbMlwN8D73f3oQLn4u4fd/dt7r6tu3vyLxeRYnzoXGa70BBfrfSgIO+B3YjzUAD33xoe5js+zZki5VPJAHUCWBXa7wNORT3HzOoJgtNfu/sXKthOWcA8mcAvZZ/vsfae3OPu87ZIbCHWsQprCJI8fPQyPhAt2Ny7uZdYeu5t37ELmYUbRSqpkgHqWWCDma01swbg7cDjeec8Djyazua7F7jk7qctmIX+c2CPu//PCrZRFjgfOgOeAsAWL8UamnOOD1weZWR0HIBFDXV0tTXNehvLycywldl8o9SxlyJd17q4ka3rlmX2f7BTyRJSeRULUO6eAN4HfJ0gyeFv3X2Xmb3HzN6TPu0J4BBwAPgE8N706/cBvwi8ycxeTH89XKm2ysLlA9n5FOvom3T8aI1k8IXFVt+e2U4dez7ydfdtzd4fDfPJbKir5Ju7+xMEQSj82sdC2w78eoHrnqLw/JRIWeUGqMmJoq+cyFaYWLeyY1baVGmx3q0Qr4fkOD54Cr90BmtbUfS6ezb38mdfeoFkMsXBkwOcOn+ZlV0ts9BiWahUSUIWNB/IDlUVC1Ab+yaXQJqPrL6RWE94mO/FSNctaWrgjg3ZQPbN7cWXjxe5HgpQsqCFA1Qsb4jP3dl/PBugNqyqjQAFELvhjsx2soRhvrdsW5vZ/vbzRxgbT5a1XSJhClCyYPm1K/jV9BxTvB7ynoE6c3GYK1fHAFi8qIGVS+fPIoXFxPq2ZlfZvXAUv3KxyBWBOzasoLs9SCS5cnWMp1VAVipIAUoWrJz5p/aVWCyeczw8vLehr6MmEiQmWONiYj2hNaIi9qJiMcvpRX392UNlb5vIBAUoWbDCASrWPv380401Mv8UFlv9qsx26tgLka97y51ricWCXx37jl3g6JlLRa4QmRkFKFmwiiVI7Dt+IbO9adXSWWnTbAoCVNArTJ07mB3uLKJ9ySLuuXllZl+9KKkUBShZsFI5ASo3QWJsPMmRUM9gQ19tpJiH2aIWYstvTO85qePRHtoFePCe7OKHT75wlKvXxsvcOhEFKFmgPJXCB0MljvJ6UEfODJJMBhUmepYuoaW5cVbbN1tiq7PZfKmj0bP5tqzpojf9DNS18QTf36EHd6X8FKBkYbrSD8kgQ8+aWrFFuQ+c7gunl9fg/NOEnKoSp/cG60RFYGb8xF3rMvtf/dHByIVnRaJSgJIFKZWTwTe5xFFuBl/tBihb3BlaaddJ7v9+5Gvf+KobaKgPMh+Pnb3Es3tPF7lCpDQKULIgFa0gcbz2KkhMJbbpDZnt1IGn8GQi0nWLmxp48O7sXNTffHu3elFSVgpQsiBNVyT20vA1zg0OA1AXj7Omp302mzbrYn23Yk3BQox+dYjU8RcjX/vI/Ruprwt6UUfODKoXJWWlACULkg9OXSR2fyi9fN3Kduritf1jYvE6Yhtem9lP7ftu5GvblyzioXvUi5LKqO2fPJECfPxadpl3i02q5L2QhvcmxDfenyl9lDq7n9Rg9J6QelFSKQpQsuDk9J5al2Px+pzj+0/UZoHY6VhzB7FVt2b2U/tL60U9eHc2o0+9KCkXBShZcHyaB3STyRQHTg5k9hdKDwogvvH1me3Uwafx8WuRr33k/k3qRUnZKUDJguP92dI8+Uts7D12IVMVoaOlKVO5eyGwnpuxlmBZdx8fJXX4R5Gv7WjJ7UX91dd3Mp7QUhxyfRSgZEFxd1Kndmf2bcXGnOPP7j2V2d62qaemKpgXY2bEN2V7Ucld34iccg7wM6/dRFNjMFx66sJlvvC9fWVvoywsClCyoPjgyUxRVGtoxpauyR5zzxmaujtUEHWhiN34Gqwh6DX65X5Sr0R/cLd9ySLe+cAtmf2//95eTp2/XPY2ysKhACULiod7Tz03Y7Hsj8Dx/sucuXgFgMb6Orau6550fa2zhmbiWx/K7Cd3fAUfH418/VvvWsf63qCwbjKZ4uNfekEJEzJjClCyoISH92Irt+Qce3ZPdnjvVRtXZCb9F5rYpjdgzUGQ8dHLJHd/K/q1MeNfPnInll7GY+ehcyokKzOmACULhifGSJ19JbMfW3lzzvHw/NPdN/XMWruqjdU1EL/9pzP7yV3fwK8ORb5+bU87P/nqGzP7f/HVl7g8Ej0jUGSCApQsGH52P6SCSX9r68EWZ1PILw5dzRSIjZlxx8YVBd9joYituxdrSwfpxDWSO54o6fp3vGULS1ubABgavsaf/P12DfVJyRSgZMHIHd7bnHNs+75scsTNN3TV7PpPUVksRt2d/zyzn9z/XXzobOTrFzXU8a5/ll1S/rn9p/n8k3vK2kapfQpQsmBMF6DCw3t3LcDsvUKsdyuxZemhOk+ReOovSko7v+umlTxyXzaN/2+/vYfn958pdzOlhilAyYLgwxfxS+leUqwOW579xTk6lmDHwf7M/kJMLy/EzIjf9bZsjb7zh0nuLG2o750P3MKWNUE2pON86PM/5uzAcNnbKrVJAUoWhNSp7PBSbPkGrK4hs//CK2dIJIOqB6uXt7G8Y/Gst69axZaupu72RzL7yR1PkDqzP/L18XiMf/22e+hsCeajhkfH+O//39MMXx0re1ul9ihAyYKQOrUrsx3ryc3e+3EovfyuTQs3e28qsVvemrPqbuKpT+HXoveC2pcs4t++/V7i6WVLjpwZ5A8ee0pBSopSgJKa56kUfjrbg7Le7PNP/YMj/ODlbPFYDe9NZmbU3fdLWEPQs/SRARLPfKakrLxNq5fya6GkiVdOXOQPHnuKkdHxsrdXaocClNQ8P/cKPjYCgDW1Yu3ZBQq/+P19JJMpIPglOlEFQXLZ4k7ir/nFzH7q6PMkn/2bkoLUm+9cy7t+KjdI/f6nv68gJVNSgJKal9z5tcx2rO+2TAHYC0NX+dZzhzPH3vbGzQuqOGyp4qtflbMkR3Lvd0g+/4WSgtSD96yfFKT+nz//bqbElEiYApTUtFT/IVKnJ9LLjfiWBzLH/iHUe9rQ18mt65fNQQvnl/jdbyN2w52Z/eSub5B88fGS3uPBe9bzq6EgdfTMIP/uo/+U8yyaCChASY1LvvTlzHZs3d1Y63IgqBzxze3Z3tPPq/cUicXi1N3/fxPruy3zWnLnEySe/wKeSkV+n4fuWc97f+bOTOLEyOg4f/iZH/DZb+0ilVLFCQkoQEnNSp0/HMreM+q2Ppw59o9P7c8sqLe+t4NXbVg+By2cnyxeR93r30WsN7u0RvLlr5P41ocyS5lE8eY71/JffvUNmZJIAH/33T38u4/+E3uPni9rm2V+UoCSmpV86SuZ7djau7C2oL7e4JVRvv5sdlXdn3vDzeo9lcji9dS9/tdyKsKnzuxj/Et/QOp09JJGN/Z18kfvfUvO8OqRM4P8p08+yZ984VkGr0Rf6kNqjwKU1KTUhaOkTu5M7xnxdO8pkUzxoc//ONN7WrOinW169mlGrK6Buje9j/itPwnp5TV89DLj3/xjEj98DB8ZiPQ+rYsb+Z1HX8s7H7glZ4mTJ184ynv/59f4iydeon9wpBLfglQ5q6UKw9u2bfPt27fPdTNkjnkqReLbf5KpvRdbs436170Ld+d/ff7H/GBndn2i3/6F+7hTAeq6pU7tDh7gHQ2toBuvJ37zm4jf8mBmld5izg0M8+mv7eCZ3SdzXo+Z8ZqtffzkvTeyoa9TPd4aY2bPufu2Sa8rQEkt8WSCxA8/TerwjzOv1f+zD2LtK/nLr+7gy09n14N6+5u38HNvuLnQ28gM+MggiR8+llO1A8DqFhFbe1ewnHzX2kjB5cUDZ/nLJ17ieP/kdai625u5f+sq7tu6ijUr2hSsaoAClNQ8T46T+N4nSR1/MfNa/KY3Er/rbfz99/by2W9lf3G+9e71vOunbtcvtwpIndpN4vkv4hePTTpmbT3E1mwj1nMz1rUGi029arG78/z+Mzz+g1d4+fC5gud0tDSxdV03W9ctY8vabpa1N+v/6TykACU1zcdGSHzvE7lLamx8HS93vom//tYuDp8ezLx+7+Ze/s3b7iUW0y+ySnF3Uke2k3zpS1OuI2X1TVjPJmJL12KdfVjnKqypreC5B08O8LUfH+SZ3SenrTyxpKmBtT3trO1p54blbfQsXULP0iW0NDcocFWxOQlQZvYg8MdAHPiku//XvOOWPv4wMAL8krs/H+XaQhSgFpaJGnvJg0+TOvZCsFquB8tnHO+8i78ZuIndeenKm9d08zuP3k9D/dR/uUv5uDt+7iCpAz8gefQ5SEy/9Ls1LsFaumFJF9bShTV3QFMrtqgFW9RKoq6Jl44N8dTOkzy3/3TkMknNi+rpamtmaWsTS1ub6Gxtom1xIy2LG2ltbmBJUwPNi+ppbqyjubE+83yWzI5ZD1BmFgf2Aw8AJ4BngXe4++7QOQ8Dv0EQoO4B/tjd74lybSEzDVAXz57h4LPfK/k6qQB3IPg3ae64p7BUEvMElkpg41epuzZI3eggDdcGscQoiaSTSKYYT6S4em2cp/xWnrZbIfQXc31dnIfvWc/Pv2kzixrq5uibW9h8fJTU8R346d2kTu/BRwZn/F5Wtwivb2JozDl/JcnZy+OcuzzOaAKSxNNfMVIYKWKZL8dIYTn/dbOJf3GZ/8ZicerrYtTXxamvixOPx6iri1EXixGPx4jFYsRjMWIxIx4zYrEYMTNiMSNmhsWMGME/QbMYWJDokdkn+8/TzDK9u4l/seHOXrjnl98JNKL3CqfrQF5v73LFpttYtWFT8ROn/vyCAaqSP6l3Awfc/VC6AZ8DHgHCQeYR4DEPouQzZtZuZj3AmgjXls35E0cZ++FnKvHWMov6rYPttoU9sXWZ12KxGA9sW8vPvv4mOkMPhMrss/pFxNfdDevuDnpWl84EhXwvnsAHTpAaOFG0hzXBE6OQGKUVaG2CdU1AN1xLJBkZHWNkNMHoWIJr4wlGx5LTV6codCgFRF88uOx8iu1qdSqZvK4ANZVKBqhe4Hho/wRBL6nYOb0RrwXAzN4NvBtg9erV19dimXdGbBF7WMcuW895CyqRty1exIa+Dm7s6+S1t65iReeSOW6l5DMzrL0H2rMp/u4OwwP4lfOZL0YG8atD+OgQjF6GsRF8fIqHdw0a6+M01jfR0RJ63WE8kWQskWQ8kcr8N5HM/UqmnGQqRSrp8yIoLASVDFCF+oz5/9+nOifKtcGL7h8HPg7BEF8pDZywpLOLkzfcO5NLpSzyxy3SD31igJGyOB6L47E6krEGxhvbSTS2kWxsx5pb2NTcyJ3NjSxZVM/Krha62po0IT4PmRks6cSWdAIbpzzPUylIjMLYVTwxFvS6EmOQHINkAk+OQ3IcUknwJKRS1KUSNKVSuKfAU9mh5PB+errDUxNBK8l4IvhKTQSw9Osp9+ArGQQ1d4L99LZ7EOTcU5m3dveJj8Dx4ONDv9ZySxnmBUmf/AswvG95UzWTQmyFI25375qKvG8lA9QJYFVovw84FfGchgjXls3KtetZ+Wu/Xam3F5EyslgMGpqhobmEGZjS1FfofaU0lUxVeRbYYGZrzawBeDuQX5f/ceBRC9wLXHL30xGvFRGRGlaxHpS7J8zsfcDXCVLFP+Xuu8zsPenjHwOeIMjgO0CQZv7L011bqbaKiEj10YO6IiIyp6ZKM9fTaCIiUpUUoEREpCopQImISFVSgBIRkaqkACUiIlWpprL4zKwfOFriZV3A+aJnLQy6F7l0P7J0L3LpfmSV417c4O7d+S/WVICaCTPbXii9cSHSvcil+5Gle5FL9yOrkvdCQ3wiIlKVFKBERKQqKUClK6ELoHuRT/cjS/cil+5HVsXuxYKfgxIRkeqkHpSIiFQlBSgREalKCyZAmdmDZrbPzA6Y2QcKHDcz+3D6+A4zu2Mu2jkbItyLd6bvwQ4z+6GZ3TYX7ZwNxe5F6Ly7zCxpZj87m+2bbVHuh5m9wcxeNLNdZvbd2W7jbInwc9JmZl8ys5fS9+KX56Kds8HMPmVm58zs5SmOV+b3Z7AMcW1/EawpdRBYR7Ba70vA5rxzHga+SrD++L3Aj+a63XN4L14DdKS3H1rI9yJ03rcJ1i/72blu9xz/22gHdgOr0/vL5rrdc3gvfhv4b+ntbuAi0DDXba/Q/XgdcAfw8hTHK/L7c6H0oO4GDrj7IXcfAz4HPJJ3ziPAYx54Bmg3s57ZbugsKHov3P2H7j6Q3n0G6JvlNs6WKP8uAH4D+Hvg3Gw2bg5EuR//F/AFdz8G4O61ek+i3AsHWszMgCUEASoxu82cHe7+PYLvbyoV+f25UAJUL3A8tH8i/Vqp59SCUr/PXyH4y6gWFb0XZtYL/B/Ax2axXXMlyr+NjUCHmT1pZs+Z2aOz1rrZFeVe/ClwM3AK2An8prunZqd5Vacivz8rtuR7lbECr+Xn10c5pxZE/j7N7I0EAer+irZo7kS5Fx8C/oO7J4M/lGtalPtRB9wJvBloAp42s2fcfX+lGzfLotyLtwIvAm8C1gPfNLPvu/tQhdtWjSry+3OhBKgTwKrQfh/BXz2lnlMLIn2fZnYr8EngIXe/MEttm21R7sU24HPp4NQFPGxmCXf/h1lp4eyK+nNy3t2HgWEz+x5wG1BrASrKvfhl4L96MAlzwMwOAzcBP56dJlaVivz+XChDfM8CG8xsrZk1AG8HHs8753Hg0XQ2yr3AJXc/PdsNnQVF74WZrQa+APxiDf5lHFb0Xrj7Wndf4+5rgL8D3lujwQmi/Zz8I/BaM6szs2bgHmDPLLdzNkS5F8cIepKY2XJgE3BoVltZPSry+3NB9KDcPWFm7wO+TpCd8yl332Vm70kf/xhBhtbDwAFghOCvo5oT8V58EFgK/O90zyHhNVi5OeK9WDCi3A9332NmXwN2ACngk+5eMPV4Pov4b+MPgL80s50EQ1z/wd1rcgkOM/ss8Aagy8xOAL8L1ENlf3+q1JGIiFSlhTLEJyIi84wClIiIVCUFKBERqUoKUCIiUpUUoEREpCopQEnNSFcbf9HMXjazz6ef0ynn+z9pZtOm25vZ+8Ofa2ZPmFl7OdsReu/bzezhaY4fMbOuAq+bmX3bzFqnufaTZra5XG2d4jM+Z2YbKvkZMr8pQEktuerut7v7LcAY8J45aMP7gUyAcveH3X2wQp91O8GzJ6V6GHhpupI87v6r7r57pg2L6KPAv6/wZ8g8pgAlter7wI1m1mlm/5Beo+aZdAknzOz3zOyv0j2JV8zsXenX32BmX554EzP7UzP7pfw3N7OPmtn29DpA/zn92r8CVgLfMbPvpF/L9GLM7F+ne3cvm9n706+tMbM9ZvaJ9Ht9w8yaCnzez6Wve8nMvpeubvD7wNvSvca3mdnS9PUvmNmfUbg+GsA7CSpCYGaLzewr6fd92czeln4901s0s18xs/3p1z5hZn+afv0v0/fhO2Z2yMxeb8G6QXvM7C+nu1eh/0dvMbMFUTBASqcAJTUn/QvvIYIK0/8ZeMHdbyVYv+ex0Km3Aj8JvBr4oJmtLOFj/lO6usatwOvN7FZ3/zBB/bE3uvsb89p0J8HT9fcQrJfzLjN7VfrwBuAj7r4FGAT+zwKf90Hgre5+G/DT6SUgPgj8TbrX+DcET/c/5e6vIig9s3qKtt8HPJfefhA45e63pXueX8tr90rgd9JtfoCg1lxYB0Gx1N8CvgT8L2ALsNXMbp/qXgGkK38fIKjlJzKJApTUkiYzexHYTlAn7c8JKrH/FYC7fxtYamZt6fP/0d2vpsvTfIdgDaCoft7MngdeIPiFXGy+5n7gi+4+7O5XCGodvjZ97LC7v5jefg5YU+D6HxCU1XkXQemdQl4HfAbA3b8CDExxXqe7X05v7yToxfw3M3utu1/KO/du4LvuftHdx4HP5x3/UrpY6k7grLvvTAeeXaHvY7p7dY6g1ykyibrWUkuuuvvt4RfMCq6R4Xn/Db+eIPcPt0X5F5vZWuDfAne5+0B6OGvSefmXTXPsWmg7SbCMRW7D3N9jZvcQ9PheDPVOJp1apB0ACTOLuXvK3fene3cPA39oZt9w99+P2O5w21N530cKqItwrxYBVyO0WRYg9aCk1n2PYM4FM3sDwVIRE8kBj5jZIjNbSlAI81ngKLDZzBrTPa03F3jPVmAYuGRBFeuHQscuAy1TtONnzKzZzBYTLIL4/ajfhJmtd/cfufsHgfMESxvkf1b4e32IYPitkH0ES5lPDOGNuPtngD8iWNY77McEw3Id6aHTQsOP05nuXkGwAOKuEt9TFgj1oKTW/R7wF2a2g6DK8r8IHfsx8BWCuZo/cPdTAGb2twTVul8hGJbK4e4vmdkLBL9YDxEMv034OPBVMzsdnody9+fTvYeJtYI+6e4vmNmaiN/H/0inZBvwT8BLBMOYH0gPa/4hwXzbZ9PDad9NHy/kKwQB+QCwNf3eKWAc+Jd53+tJM/svwI8I5td2A/nDgFOa7l6lA9bVGl3WRspA1cxlQTKz3wOuuPsfzXVbZpuZ9QCPufsDEc9f4u5X0j2oLxIsPfHFMrTjt4Ahd//z630vqU0a4hNZYNI9lk/YNA/q5vm9dC/tZeAw8A9lasog8OkyvZfUIPWgRESkKqkHJSIiVUkBSkREqpIClIiIVCUFKBERqUoKUCIiUpX+f9npvKLAhc+VAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pmf_std = marginal(posterior, 1)\n",
    "pmf_std.plot()\n",
    "\n",
    "pmf_std2 = marginal(posterior2, 1)\n",
    "pmf_std2.plot()\n",
    "\n",
    "decorate(xlabel='Population std (sigma)', \n",
    "         ylabel='PDF')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.7800000000000002, 0.4660674157303371)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pmf_mean2.idxmax(), pmf_std2.idxmax()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'lifelines'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_195977/4180200288.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mimport\u001b[0m \u001b[0mlifelines\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'lifelines'"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Think Bayes, Second Edition\n",
    "\n",
    "Copyright 2020 Allen B. Downey\n",
    "\n",
    "License: [Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}