{ "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": [ { "name": "stdout", "output_type": "stream", "text": [ "Downloaded utils.py\n" ] } ], "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": [ { "name": "stderr", "output_type": "stream", "text": [ "In /home/downey/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The text.latex.preview rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /home/downey/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The mathtext.fallback_to_cm rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /home/downey/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: Support for setting the 'mathtext.fallback_to_cm' rcParam is deprecated since 3.3 and will be removed two minor releases later; use 'mathtext.fallback : 'cm' instead.\n", "In /home/downey/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The validate_bool_maybe_none function was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /home/downey/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The savefig.jpeg_quality rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /home/downey/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The keymap.all_axes rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /home/downey/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The animation.avconv_path rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /home/downey/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The animation.avconv_args rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n" ] } ], "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": 10, "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": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
TIME(In Days)AGE OF PATIENTSCENSORAGE AT MENARACHEBREASTFEEDCONTRACEPTDETECTIONNEOADJUVANT
048441181.5212
1898430171.0112
236400152.0212
323561151.5111
4300530121.0212
\n", "
" ], "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": 10, "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": 32, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "log_duration = np.log10(duration)" ] }, { "cell_type": "code", "execution_count": 33, "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": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0u0lEQVR4nO3deXxU9b3/8dcnIWFNWELYwo7IqqAEqDtudWkVXFoXpNXWWmvpvvfWXr3e295fva2/21qvP2u92toWatXWWqrVKi61IlEWBQTZCWsIkIQEsn5+f5yTYRImC5DJzCTv5+ORR842Zz5nDsw7Z/t+zd0RERFJNmmJLkBERCQWBZSIiCQlBZSIiCQlBZSIiCQlBZSIiCQlBZSIiCQlBZS0KTN70MzubKN1DTezg2aWHo4vNrNb22Ld4fr+amafbKv1dVRmdo6ZrW2D9cwys8L2fl9JXV0SXYCkDjPbDAwEaoBaYDXwK+Ahd68DcPfbj2Fdt7r7i00t4+5bgV4nVnXk/e4CTnL3m6LWf1lbrLujc/fXgHHxfh8zc2Csu69vz/eV5KUjKDlWV7h7FjAC+E/gW8Av2/pNzKzD/PFUfwSYrDrSZy0diwJKjou7l7j7M8B1wCfNbDKAmT1qZv8eDvc3s2fN7ICZ7TOz18wszcx+DQwH/hyewvummY00MzezT5vZVuClqGnRX6BjzOwtMysxsz+ZWb/wvY46fWRmm83sIjO7FPgucF34fivC+ZFThmFd3zOzLWa2x8x+ZWa9w3n1dXzSzLaa2V4z+5emPpvwM/gfM1tkZuXA+Wb2ETNbZmalZrYtPKKrX/4xM/taOJwXvtcd4fhJ4WdnMd7nJDN7Jfws9prZwkb1dolaNnpbbzazf5jZfWa2D7gn3EeTo5bPNbNDZjYg+rM1s2+b2R8a1fHfZvbTcPgWM1tjZmVmttHMPtvU59RoHa+GgyvCfXRd430a7s9vmNlKMys3s1+a2UALTtWWmdmLZtY3avkPmdkb4batMLNZUfNuDusrM7NNZja3NXVK+1JAyQlx97eAQuCcGLO/Fs7LJTg1+N3gJT4P2EpwNNbL3X8U9ZrzgAnAJU285SeATwFDCE41/rQVNT4H/ABYGL7flBiL3Rz+nA+MJji1eH+jZc4mOOV0IfB9M5vQzNveCPwHkAW8DpSHtfcBPgJ8zszmhMu+AswKh88DNoa/Ac4FXvPYbZLdA/wN6AsMBX7WTD2NzQzfZwDwb8BTwA1R8z8OvOLuexq97nfA5WaWDZGjw48Dvw3n7wE+CmQDtwD3mdnpLRXj7ueGg1PCfbSwiUWvAS4GTgauAP5K8O+qP8H32RfDuvKAvwD/DvQDvg48GQZvT4J/N5eFZwPOBJa3VKO0PwWUtIUdBF8CjVUDg4ER7l7t7k190Ua7y93L3f1QE/N/7e7vuXs5cCfwcWubU2hzgZ+4+0Z3Pwh8B7i+0dHb3e5+yN1XACuAWEFX70/u/g93r3P3w+6+2N3fDcdXEnzR14fQK8A5ZpZGEEg/As4K550Xzo+lmuBU65DwPV4/hu3d4e4/c/ea8LP+LQ0D6kaOhE6Eu28B3gHmhJMuACrc/c1w/l/cfYMHXiEI0Fh/vByvn7n7bnffDrwGLHH3Ze5eCTwNnBYudxOwyN0XhZ/5C0ABcHk4vw6YbGbd3X2nu69qwxqljSigpC3kAftiTL8XWA/8LTyd8u1WrGvbMczfAmQQ/PV8ooaE64tedxeCI796u6KGK2j+Bo4G22FmM83sZTMrMrMS4HbCut19A3AQmErwZf4ssMPMxtF8QH0TMOAtM1tlZp9qdgubqQ94Cege1jkirOXpJl4bHWYNgszMLjOzN8PTkgcIAqEt9k+93VHDh2KM1++TEcDHwtN7B8JazgYGh3/cXEewD3aa2V/MbHwb1ihtRAElJ8TMphME1FF/vbt7mbt/zd1HE5yO+aqZXVg/u4lVtnSENSxqeDjBUcReglNoPaLqSic4tdja9e4g+FKLXncNDb8Aj0Xj9/st8AwwzN17Aw8ShEu9V4Brgczw6OAVglOCfWni9JO773L3z7j7EOCzwANmdhLBZwFRnwcwqLn6wrswf08QPDcCz7p7WRPb9gQwy8yGAleF24aZdQWeBP4LGOjufYBFjbazvWwjONruE/XT093/E8Ddn3f3iwmO8N8HfpGAGqUFCig5LmaWbWYfBRYAj7v7uzGW+Wh4Id+AUoJb02vD2bsJrvUcq5vMbKKZ9SC4dvIHd68F1gHdwpsRMoDvAV2jXrcbGBmeRovld8BXzGyUmfXiyDWrmuOoMZYsYJ+7HzazGQQhEO0VYD5Qf7PAYuALwOvh9h3FzD4WhgTAfoLQqXX3ImA7wWeVHh5ZjWlFjb8lOLKYS4zTe/XC9S8G/hfY5O5rwlmZBJ95EVBjZpcBH27F+9Y73n8TsTwOXGFml4SfQbfwpouh4Y0VV4bXoioJjl5jfsaSWAooOVZ/NrMygr9Q/wX4CcHF8FjGAi8SfAH8E3jA3ReH834IfC88/fL1Y3j/XwOPEpxu60Z4UdzdS4A7gIcJvpzLCW7QqPdE+LvYzN6Jsd5HwnW/CmwCDhMERFu5A/i38LP7PsHRSrRXCEKsPqBeJzgCepWmTQeWmNlBgqOzL7n7pnDeZ4BvAMXAJOCNlgp09yUEn9sQgpsPmvNb4CKigiw84voiwbbtJwjhZ1p63yh3AY+F/yY+fgyvO4q7bwNmE9xAUUTw7/UbBN95aQQ38OwgODV9HsH+kSRj6rBQRESSkY6gREQkKSmgREQkKSmgREQkKSmgREQkKaVcI5H9+/f3kSNHJroMERFpI2+//fZed89tPD3lAmrkyJEUFBQkugwREWkjZrYl1nSd4hMRkaSkgBIRkaSkgBIRkaSUctegYqmurqawsJDDhw8nuhSJoVu3bgwdOpSMjIxElyIiKaRDBFRhYSFZWVmMHDkSO7rjUUkgd6e4uJjCwkJGjRqV6HJEJIXE7RSfmT1iQdfZ7zUx38zsp2a2PuzCucVeN5ty+PBhcnJyFE5JyMzIycnR0a2IHLN4XoN6FLi0mfmXEbR2PRa4DfifE3kzhVPy0r4R6XjcnZ3FB1n2wa6WFz5OcTvF5+6vmtnIZhaZDfwq7AL8TTPrY2aD3X1nvGoSEZHjU11Ty8YdB3h/azHvb9nL+9uKKS2vpGtGF379L1eSnt72xzuJvIsvj4bdTheG045iZreZWYGZFRQVFbVLcceiuLiYqVOnMnXqVAYNGkReXl5kvEePoFPTzZs3Y2bceeedkdft3buXjIwM5s+fD8Bdd93V4LVTp07lwIEDACxbtoxbb70VgEcffTTymtZKT09n6tSpTJo0iSlTpvCTn/yEurq6Ntj6QFVVFeeeey41NW3Vv5+IJFJpeSVvrdnBr59/l+/+4mVu+vdn+O4vXuZXz6/krfd3UFpeCUBldQ1bdpfEpYZE3iQR67xPzM6p3P0h4CGA/Pz8pOvAKicnh+XLlwNByPTq1Yuvfz3og69Xr16R5UaPHs2zzz7LPffcA8ATTzzBpEmTGqzrK1/5SuS10X7wgx/wve9977hr7N69e6TGPXv2cOONN1JSUsLdd9993OuMlpmZyYUXXsjChQuZO3dum6xTRNpXYVEpi5dtYcnqHewoLmtx+R7dMhg3LIe6uvh8LScyoAqBYVHjQwl6uOywunfvzoQJEygoKCA/P5+FCxfy8Y9/nB07mt/ssrIyVq5cyZQpU46ad/PNN5OdnU1BQQG7du3iRz/6Eddee22z6xswYAAPPfQQ06dP56677mLLli3MmzeP8vJyAO6//37OPPNM5s2bx7XXXsvs2bMBmDt3Ltdddx1jxozhlltuoaqqirq6Op588knGjh3LnDlz+M53vqOAEkkh5Yeq+Md7hbz0zmY+KNzX7LKD+vVi3PAcxg/PYdzwHIYPyI7rNeZEBtQzwHwzWwDMBEra4vrTNXf+4YQLa8qT9zT/xd8a119/PQsWLGDQoEGkp6czZMiQBgF133338fjjjwPQt29fXn75ZQoKCpg8eXKT69y5cyevv/4677//PldeeWWLAQXB0VxdXR179uxhwIABvPDCC3Tr1o0PPviAG264gYKCAm699Vbuu+8+Zs+eTUlJCW+88QaPPfYYX/nKV/jSl77E3Llzqaqqora2FoDJkyezdOnSE/yERKQ9bN9bxu9fXs2S1Tuorqk9an56ehqjB/dhwoj+jBuew7hhOfTN6tauNcYtoMzsd8AsoL+ZFQL/CmQAuPuDwCLgcmA9UAHcEq9aksmll17KnXfeycCBA7nuuuuOmh/rFN/OnTvJzT2qod+IOXPmkJaWxsSJE9m9e3erawnuTwkedJ4/fz7Lly8nPT2ddevWAXDeeefx+c9/nj179vDUU09xzTXX0KVLF8444wz+4z/+g8LCQq6++mrGjh0LBNe5MjMzKSsrIysrq9V1iEj7Kimv5LsPvczBQ1UNpqelpTF93GBmnTaCqScNJDMjPUEVBuJ5F98NLcx34PPxev9klZmZybRp0/jxj3/MqlWr+POf/9zia7p3797sc0Rdu3aNDNeHTks2btxIeno6AwYM4O6772bgwIGsWLGCuro6unU78lfSvHnz+M1vfsOCBQt45JFHALjxxhuZOXMmf/nLX7jkkkt4+OGHueCCCwCorKxs8HoRST6L3lzfIJxGDurDBaeP5OxTh9G7Z9dmXtm+OkRLEtHa4jRcvH3ta1/jvPPOIycnp1XLT5gwgR//+Mdt9v5FRUXcfvvtzJ8/HzOjpKSEoUOHkpaWxmOPPRY5ZQfBNa4ZM2YwaNCgyA0dGzduZPTo0Xzxi19k48aNrFy5kgsuuIDi4mJyc3PVpJFIEtu17yDPLdkQGb9jzjQunJacrbx0uIBKBZMmTTrq7r160degAP74xz8yfvx4SkpKTujU2aFDh5g6dSrV1dV06dKFefPm8dWvfhWAO+64g2uuuYYnnniC888/n549e0ZeN3DgQCZMmMCcOXMi0xYuXMjjjz9ORkYGgwYN4vvf/z4AL7/8Mpdffvlx1Sci8bdtTyl3P/pa5Ogpt08PZk0dkeCqmmatPSWULPLz871xh4Vr1qxhwoQJCaqofdx3331kZWVFnoVqLxUVFZxyyim888479O7du9llr776an74wx8ybty4o+Z1hn0kksy27Snlzl++QllF8PxSRpd0vnvTWZw6ZkCCKwMze9vd8xtPV3cbKeJzn/tcg2tN7eHFF19k/PjxfOELX2gxnKqqqpgzZ07McBKRxHv0rysi4dQ1owvf+8TZSRFOzdEpvhTRrVs35s2b167vedFFF7F169ZWLZuZmcknPvGJOFckIsdj255Slq8P7vA1jH+9+RzGDW/dNfBE6jBHUKl2qrIz0b4RSay/Ld0YGZ4+fnBKhBN0kIDq1q0bxcXF+iJMQvX9QenWc5HEcHfeeK8wMn7pzDEJrObYdIhTfEOHDqWwsJBkbEhWjvSoKyLtb+WGPRw4GDxHmd2zK6eMTu7rTtE6REBlZGSot1YRkUbKKip54I9vR8ZnTsgjLS11+mfrEAElItKZuDsHD1VRUl5JaXklJeWVlByspLQi+F0/fde+gxSXHgKClsevPje17rJVQImIJJi7U3G4umHglB8dOPW/S8srqTvGa+5fuHo6A/r2bHnBJKKAEhFpJxWHq1m0ZD3bi8oaBk9FJbW1bdeBaGPXzprAjAlD4rb+eFFAiYi0k8dfeI/n39rQ8oKt0L1rBr17diW7Z9fI7z6Nxnv37Eq/7O5kJ1EDsMdCASUi0g7cnX+uKmxyfteMLkcHTq+jA6f+d0aXxHaF0R4UUCIi7WDD9v2UlgdNDWX37Mr8q/IbBE7XTH0dN6ZPRESkHRSs2xUZPv3kQUwbNziB1aSGDtGShIhIMquprWPJ6u2R8dNPVji1ho6gRETiaF/pIe5d8CZbd5cAkGbGlCRvRTxZKKBERNqYu1NYVMaqTUU8sXhNpKkhgI+eOZZe3TMTWF3qUECJiJwgd2fLrhJWbd7L6s1FrN6yN3JDRD3DmHfJKVx51tgEVZl6FFAiIsdh9/5ylqzezurNe1m9eS/lh6uaXDarR1e+dt3MlGqoNRkooEREjlHB2p3c+7s3qamtbXKZXt0zmTiiPxNH5XLulOH0TtGHZRNJASUicgw27tjPjxceHU7ZPbsyaWQuE0f2Z9LIXIYPzMYsdVoOT0YKKBGRVtpfdpgfPP4GVdVBOOX26cHV545n4sj+5PXPUiC1MQWUiEgrLfj7KvaXHem+4l/mnc2wAdkJrqrj0oO6IiKtsGVXCX9/e3Nk/MvXzlA4xZmOoEREmuHuvFiwif/960qcoA+mU8cMUFNF7UABJSLShL0lFTzw9Nus2LA7Mq1LejqfuOTUBFbVeSigRERiqKtzvvfwYooOVESmDcnJ4gvX5DNqcJ/EFdaJKKBERGLYW1IRCSfDuOKssdxw4SQyMzp+P0zJQgElIhJDbZ1Hhgf07cEnL9Vpvfamu/hERGKo8yMBlZ6mr8pEiOunbmaXmtlaM1tvZt+OMb+3mf3ZzFaY2SozuyWe9YiItFZNbV1kOC1ND+AmQtwCyszSgZ8DlwETgRvMbGKjxT4PrHb3KcAs4MdmpnboRSTholsjV0AlRjyPoGYA6919o7tXAQuA2Y2WcSDLgvZBegH7gJo41iQi0irPLdkQGR4zpG8CK+m84hlQecC2qPHCcFq0+4EJwA7gXeBL7l6HiEgC7dp3kCWrd0TGrzhTfTglQjwDKtYxsTcavwRYDgwBpgL3m9lRbYeY2W1mVmBmBUVFRW1dp4hIA6s2FTVoNWLEoN4JrqhzimdAFQLDosaHEhwpRbsFeMoD64FNwPjGK3L3h9w9393zc3Nz41awiAjAB4X7I8OTR6mTwUSJZ0AtBcaa2ajwxofrgWcaLbMVuBDAzAYC44CNcaxJRKRFm3cdiAyPydP1p0SJ24O67l5jZvOB54F04BF3X2Vmt4fzHwTuAR41s3cJTgl+y933xqsmEZGWuDu79pVHxvP6ZyWwms4tri1JuPsiYFGjaQ9GDe8APhzPGkREWmvP/nLuf6qAsorgFvP09DRysrsnuKrOS00diUin5+78/e3NPLJoBZXVR550mTF+iJ6BSiAFlIh0au7OI4tWsOjN9ZFphnHVueP4+PkTEliZKKBEpNOqqa3j508X8OqKrZFp9V1qnDwsJ4GVCSigRKQTaxxOZ0wayhevma4uNZKEAkpEOqUtu0oahNPF+aO57YrTdM0piSigRKTTcXeeeu39yPj08UP47JWnETQLKslCASUincqhymoe+OM7vPHekaZCrzpnnMIpCSmgRKTT2LanlB/99p/sKC6LTJt28mDGDdcNEclIASUinUJlVQ3/9uhr7Cs7FJn24emjueWyKQmsSpqjgBKRTuFvBZsi4ZSZkc7tV57OeVNHJLgqaU5cu3wXEUkG6wv38dSrR26K+OQlpyqcUoCOoESkw6qrc55+bS0LXlpNXV3QF2rfrO5ccPrIxBYmraKAEpEOqaq6lh88/g/e3bgnMq1bZhe+cHW+HsRNEQooEemQ/vSPdQ3C6eRhOXzp2ukM6tcrgVXJsVBAiUiHU1Nbx3NLNkTGrzhzLPM+fArp6brsnkq0t0Skw3lrzQ4OHDwMBNecblI4pSTtMRHpcDbu2B8ZnjV1OF0UTilJe01EOpxte0ojwyMG9k5gJXIiFFAi0qGs3VrM22t3RcaHK6BSlgJKRDqMkvJK7n+qAMcBOHXMAIYPzE5wVXK8dBefiHQI67YV818L3qS4NGjOqGtGF+6Yk69WylOYAkpEUt7zb23gl4tWUFtbF5l22xWnkdunRwKrkhOlgBKRlPbyO5t56M/LIuM9u2XypWunM23c4ARWJW1BASUiKWvD9v08+MyRcBo1uA/fuOEMBvbtmcCqpK0ooEQkJdXU1vFfC96kprYWgGG52fz7rbPolqmvtY5Cd/GJSErasH0/ew6UA9C9awbfmnumwqmDUUCJSErasrskMnz6yYMYnKNGYDsaBZSIpKRNOw9EhvWsU8ekgBKRlPT+1uLI8Ni8fgmsROJFASUiKaf8UBXbdgft7aWZMW54ToIrknhQQIlIylmztTjSnNGoIX10c0QHpYASkZSzZvPeyPCE4f0TWInEkwJKRFLO6i1RATVSAdVRxTWgzOxSM1trZuvN7NtNLDPLzJab2SozeyWe9YhI6jtcVcP67Uc6JJw4QgHVUcXtxK2ZpQM/By4GCoGlZvaMu6+OWqYP8ABwqbtvNbMB8apHRDqGLbtKqKsLGoUdmptNds+uCa5I4iWeR1AzgPXuvtHdq4AFwOxGy9wIPOXuWwHcfU8c6xGRDmBf2aHIsB7O7djiGVB5wLao8cJwWrSTgb5mttjM3jazT8RakZndZmYFZlZQVFQUp3JFJBWUHKyMDPfW0VOHFs+AitVLmDca7wJMAz4CXALcaWYnH/Ui94fcPd/d83Nzc9u+UhFJGfsPHo4M98nqlsBKJN7i+fBAITAsanwosCPGMnvdvRwoN7NXgSnAujjWJSIpbG9JRWS4Ty8FVEcWzyOopcBYMxtlZpnA9cAzjZb5E3COmXUxsx7ATGBNHGsSkRS2fW8Zr608cuVA16A6trgdQbl7jZnNB54H0oFH3H2Vmd0ezn/Q3deY2XPASqAOeNjd34tXTSKSmqpralmyegdPLF4T6db9pLx+nDpaN/52ZHFtH8TdFwGLGk17sNH4vcC98axDRFLTrn0HebFgE39/ZzOl5Udujkgz4/bZp5OWFutSt3QUasBKRJLOrn0HefjZ5Sz/YHekzb16aWZ86iNTGTW4T2KKk3ajgBKRpPPIX1aw7INdDablZHfnovxRXDRtFP2yuyeoMmlPCigRSTo7issiw6eOGcDlHzqJ08cOIj1dzYd2JgooEUkqtbV1lERdb/ryx2bqgdxOqtk/R8zs0ajhT8a9GhHp9J57ayMVh6sB6NEtg6zumQmuSBKlpePlKVHDX4pnISIiu/YdZMFLqyLjc84epzv1OrGWAqpx00QiInGxs/gg3//lK5Gjp0H9enHlWWMTXJUkUkvXoIaa2U8J2tWrH45w9y/GrTIR6TTqw6m+pfKMLul8/qppZHRJT3BlkkgtBdQ3ooYL4lmIiHRO5Yeq+MGv/9EgnL4z90wmjlTD0J1dswHl7o+1VyEi0vnU1Tn/9w9vRW4rz+iSzndvOotTx6gJI2lFY7Fm9kkze8fMysOfgqb6bRIRaa19pYe451ev8c66Iw/kzr8qX+EkEc0eQYVB9GXgq8A7BNeiTgfuNTPc/Vdxr1BEOpyCtTu5/6kCyiqOPO901TnjOPvUYc28Sjqblq5B3QFc5e6bo6a9ZGbXEHThroASkWPy3JIN/OLZZZFxw7jq3HHceNGkBFYlyailgMpuFE4AuPtmM8uOT0ki0lFVHK7m8ReO9KjTN6s7X7p2Oqeo2wyJoaWAOnSc80REGig/VMXv/r6aQ5VHnnP64W3nk61mjKQJLQXUBDNbGWO6AaPjUI+IdCDuzobt+3l+6UZeW7mN6prayLzZZ5+scJJmtRRQU4CBwLZG00cAO+JSkYikPHfnleVbefafH7Bp54Gj5o/J68usqSPavzBJKS0F1H3Ad919S/REM8sN510Rr8JEJHU9/dpafhN1rane8IG9uWT6aGadNoLMDLUSIc1rKaBGuvtRp/jcvcDMRsanJBFJZSXllfxh8fuR8Ywu6Zx1ylA+nD+ak4f1w0yNv0rrtBRQ3ZqZpy4tReQov3txFZXVNQAMy83mnlvPI6uHrjXJsWupJYmlZvaZxhPN7NPA2/EpSURS1evvbuOFgo2R8ZsuOUXhJMetpSOoLwNPm9lcjgRSPpAJXBXHukQkxWzaeYAHnj7yd+uM8UOYdvKgBFYkqa6lxmJ3A2ea2fnA5HDyX9z9pbhXJiIp4511u/jxwjcjp/YG9evF/Kvzdb1JTkhLR1AAuPvLwMtxrkVEUtBzSzbw8LPL8bB/026ZXfjmDWfQU121ywlqVUCJiMSyfP3uBu3q9e/dg+/edBYjBvVOYFXSUSigROS4uDu/ffHIs05j8vrynbln0TeruZt/RVqvxf6gRERieXvdLjZs3w8Ezzp9+8YzFU7SphRQInLM3J0Ff18VGf/w9FH0y9ajkdK2FFAicsyWvr8z0sZeRpd0rjpnfGILkg5JASUix8TdWfjS6sj4pTNG69SexIUCSkSOyZurt7N51wEgOHqac864xBYkHZbu4hORVqmqruWp19by9KtrI9MunzmGPr109CTxEdcjKDO71MzWmtl6M/t2M8tNN7NaM7s2nvWIyLGrra2jYO1Ovvyzv/HEy6upqQ06HezZLZPZOnqSOIrbEZSZpQM/By4GCgkann3G3VfHWO7/AM/HqxYRaR13Z8+BCj4o3Mf6wn2s376fDTv2U1Vd22C5MXl9uWP2NHqrR1yJo3ie4psBrHf3jQBmtgCYDaxutNwXgCeB6XGsRURiKC2vZP32/UEgbd/HB4X7KauobHL5nt0yuenDk7lo2ijS0tTOnsRXPAMqj4ZdxRcCM6MXMLM8glbRL6CZgDKz24DbAIYPH97mhYp0BnV1zpbdJazdWsz7W4tZt62Y3fvLW/XanOzu5I8fwnUXTNRRk7SbeAZUrD+vvNH4/wW+5e61zbV67O4PAQ8B5OfnN16HiMRQfqiKtdv2sXZbMWu3FrNu275Ia+PN6dEtg5Py+jI2rx8nDe3HSXl99RCuJEQ8A6oQGBY1PhTY0WiZfGBBGE79gcvNrMbd/xjHukQ6HHensKiMDwr3sXZrEEiFRWWRFsabkp6exujBfRgbBtFJQ/sxJKeXusmQpBDPgFoKjDWzUcB24HrgxugF3H1U/bCZPQo8q3ASaVlZRSUfFO5n7bZiPti2jw+276PicHWLr+ub1Z3xw3MYNzyHccP6MWpwHzK6pLdDxSLHLm4B5e41Zjaf4O68dOARd19lZreH8x+M13uLdDQ7iw+ycsNu1hXuY93WfewoLmvxNWlmjBzcJwikYUEo9e/dXUdHkjLi+qCuuy8CFjWaFjOY3P3meNYikqqWr9/Nvz/2eoun67J7duXkof0YO6wf44f356S8vnTL1LP4krr0r1ckyT23ZMNR4ZSWlsaowb0ZNywnEkoD+/bU0ZF0KAookSS1t6SCf67azvL1uyPTbrhoEpNH5jJ6SF8yM3TtSDo2BZRIEtmzv5w3V2/njfcK+aBwX4N5Q3Ozufa8CQmqTKT9KaBEksAry7ew6M0NrN++L+b89PQ0PjZL4SSdiwJKJMHWbSvmp08uPWp6mhmnjB7AGZPymDExTy04SKejgBJJsBcKNkWG09LSmDImDKUJQ8jqoVCSzksBJZJAlVU1vPFeYWT8nk+dy/gR/RNYkUjyUI+6Igm0ZM0ODlcF7eMNycli3PCcBFckkjwUUCIJ9PKyzZHhWaeN0HNMIlEUUCIJsmd/Oe9uKALAMM6bqq5kRKIpoEQS5O/vbI60EDHlpAH0790jwRWJJBcFlEgC1NbW8dLbmyPjF+WPanphkU5KASWSAMvW72Zf2SEgaOR1+vghCa5IJPkooEQS4MWoZ58uOG0kXdL1X1GkMf2vEGln+0oP8fbanZHxC6aNTFwxIklMASXSzl4o2ESdBzdHTByZS17/rARXJJKcFFAi7WjVpiKefPX9yPjFujlCpElq6kikHZRVVPLOul08smgFtbV1AIwY1IczJw9NcGUiyUsBJRInu/YdZOn7O3lrzQ7e37I3cloPoHfPbnxn7pm6OUKkGQookTa0Y28Zi5dt4a33d7JtT0nMZbpmdOFbN55Bbh89mCvSHAWUSBspOlDB1x/4O5XVNTHnjx3aj+njh3DulOEKJ5FWUECJtJEnFq9pEE5d0tOZMmYA08cPJn/8EPpmdUtgdSKpRwElcgzcnaIDFewsPhjzp95nrzydc6cMp1um/ouJHC/97xFpxN0pLj0UM4B27Sunpra22ddPHjWAD08f3U7VinRcCigRoKq6lj8sXsPStTvZWXyQ6prmQ6gpQ3OzufWjU9u2OJFOSgElnZ678z9/eptXV2xt1fLZPbsyOKdX1E8WeTm9GJTTS6f0RNqQ/jdJp/f8WxuPCqde3TMbhNCQnCwG5/RiUL+e9OyemaBKRToXBZR0amu27OWRv66IjM86bQQ3X3oqWT26JrAqEQEFlHRiu/eX839++89I00MjB/Xhs1ecTmZGeoIrExFQY7HSSR2qrOaHj79BWUUlEFxX+taNZyicRJKIAko6nbo6577fvxVpiig9PY1v3XAGA/r2THBlIhJNp/ik03B3NmzfzzNvfMDb6450GHjH7GmMH9E/gZWJSCxxDSgzuxT4byAdeNjd/7PR/LnAt8LRg8Dn3H0FIm1oZ/FBXl2xlddXbmNHcVmDeXPOHses00YkqDIRaU7cAsrM0oGfAxcDhcBSM3vG3VdHLbYJOM/d95vZZcBDwMx41SSdR01tHS8WbOLlZVtYv31fzGXOmDSUuRdPbufKRKS14nkENQNY7+4bAcxsATAbiASUu78RtfybgHpvkxNWWVXDj373T5av333UvK4ZXfjQpDzOnTKcKWMGYGYJqFBEWiOeAZUHbIsaL6T5o6NPA3+NYz3SCRyqrOYHj7/B6s1FkWnp6WmcPnYQ50wZTv7Jg+iq1h5EUkI8/6fG+tPUY0zDzM4nCKizm5h/G3AbwPDhw9uqPulgSssr+Y9f/6PBKb2rzhnH7LNP1oO3IikongFVCAyLGh8K7Gi8kJmdCjwMXObuxbFW5O4PEVyfIj8/P2bISee2Z3859zz2eoObID5xyanMPvvkBFYlIicingG1FBhrZqOA7cD1wI3RC5jZcOApYJ67r4tjLdIBVVXXsnZbMe9u2MPf39nMgYOHATCMWz86lUtnjklwhSJyIuIWUO5eY2bzgecJbjN/xN1Xmdnt4fwHge8DOcAD4cXqGnfPj1dNktrq6pyNO/azcuMe3t24hzVbio/qFiM9PY0vXzuDMyfrfhuRVGfuqXXGLD8/3wsKChJdhrSDqupaNu86wAeF+1m1qYh3N+2h4nB1k8v36p7JN244g8mjctuxShE5UWb2dqyDE93OJEmhrs7ZtqeU9dv38UHhftZv38eW3aXU1dU1+7ohOVmcMmYAp4zOZcqYgfToltFOFYtIvCmgJCGKSw+xZvNe1m/fx/rt+9mwYz9V1S33Yts3qzunjM7l1NEDOGXMAPr37tEO1YpIIiigpF3V1tax8KXVPPXqWjz2UwcNDMnJYkxeX04e1o9TRg9gaG6WHq4V6SQUUNJuig5UcN8TS1i7NebTBPTL6s5JeX05aWg/Tsrry5i8vvRS77UinZYCStrFm6u38/OnCxrc5DAmry9TxwyMBFK/7O4JrFBEko0CSuKqqrqWx55byXNvbYhMSzPjugsncvU540lL0+k6EYlNASVxU1hUyo8XLmHr7pLItP69e/CVj81Q/0si0iIFlLQ5d2fxsi089OyyBnfmzZyQxx1XTdN1JRFpFQWUtKlDldU89OdlvLpia2Ral/R0brnsVC6ZMVp34IlIqymgpE0UFpXyt6UbeXnZlgY3QgzJyeJr181k5OA+iStORFKSAkqOW01tHUvW7OBvb23kvU17jpp/wekj+fRHptJN/S+JyHHQN4ccl9dXbuPR51ayv+zQUfMG9evFDRdO4uxTh8V4pYhI6yig5JiUH6rioWeX8frKbQ2mG8b08YO5ZOYYdaUuIm1CASWttnpzET99cilFByoi0/r06sZF+aO4OH+U2sUTkTalgJIW1dTW8fsY7eedf9pIPv2RKXTvqhbERaTtKaCkSVXVtWzZXcIvnl3Ghu37I9N7dsvk9tmnq1NAEYkrBZRQVlHJ9qIytu8to7CoLDK8e1/5US2OTx41gC9ck6/TeSISdwqoTsLdKTpQEQTQ3iCECotK2b63jNLyyhZfn56extyLJnPlWWN1A4SItAsFVAdTVV3LzuKDbCsqbXBUtGNvGdU1LXcIGM0wBvbrycjBfbj2vPGM0sO2ItKOFFApqrS8MgifPcFRUDBcRtGBilZ1BBgto0s6ef2zyMvNYmhuFnm52Qztn8XgnF5kZqTHaQtERJqngEoBZRWVrNxYxHsb97B1dymFRaUcPFR1zOvJ7tmVobnZkTDK6x8EUm6fHjptJyJJRwGVhKpranl/azEr1u9mxYY9bNpxoNVHRfWn5RocEYXDWT26xrlyEZG2o4BKAu7O1t2lLF+/m5UbdrNq894WrxfFPC2Xm8XgfjotJyIdgwIqQYpLD7EyPEJauWEPJeWHm1w2zYyThvbj1DEDGDcsR6flRKRTUEC1k8NVNby3qYiVG3azYv0eCotKm11+cE4vTh0zkCljBjB5VC491cmfiHQyCqjj4O5UVtdSWl4Z/FRUUVYR/m40HsyvpKy8qtnrSL26Z3LK6AFMPWkgp44ZwIC+Pdtxi0REko8CiqCtudLySsoqqiiNCpqSqGllFZWUlgfBU1JeRU3tsT1T1Fh6ehoThvdnykkDmDJmIKOH9NEpOxGRKJ0qoFZtKuLv72xuEDalFVUcqqxu+cVtYPjA3pEjpIkj+tNVHfmJiDSpU31D7i2p4JXlW9pkXV3S0+ndM5OsHl3JDn/37tmVrB6ZZPfoSlbPrmTXD/fIJLtnV7qkp7XJe4uIdAadKqCaeg7IMLJ6BmGSHQZLVo/wd6Px7J7BMl0z0nVKTkQkjjpVQI0c1Js75kw7Knx6dc9U2IiIJJlOFVD9srtz4bRRiS5DRERaQRdFREQkKcU1oMzsUjNba2brzezbMeabmf00nL/SzE6PZz0iIpI64hZQZpYO/By4DJgI3GBmExstdhkwNvy5DfifeNUjIiKpJZ5HUDOA9e6+0d2rgAXA7EbLzAZ+5YE3gT5mNjiONYmISIqIZ0DlAduixgvDace6DGZ2m5kVmFlBUVFRmxcqIiLJJ54BFeu+7caN0bVmGdz9IXfPd/f83NzcNilORESSWzwDqhAYFjU+FNhxHMuIiEgnZO6t66n1mFds1gVYB1wIbAeWAje6+6qoZT4CzAcuB2YCP3X3GS2stwg40faK+gN7T3AdiaZtSA7ahuSQ6tuQ6vXDiW3DCHc/6vRY3B7UdfcaM5sPPA+kA4+4+yozuz2c/yCwiCCc1gMVwC2tWO8Jn+MzswJ3zz/R9SSStiE5aBuSQ6pvQ6rXD/HZhri2JOHuiwhCKHrag1HDDnw+njWIiEhqUksSIiKSlDprQD2U6ALagLYhOWgbkkOqb0Oq1w9x2Ia43SQhIiJyIjrrEZSIiCQ5BZSIiCSlDh1QHaE19VZswywzKzGz5eHP9xNRZ1PM7BEz22Nm7zUxPxX2QUvbkOz7YJiZvWxma8xslZl9KcYySb0fWrkNyb4fupnZW2a2ItyGu2Msk+z7oTXb0Hb7wd075A/Bs1cbgNFAJrACmNhomcuBvxI0ufQhYEmi6z6ObZgFPJvoWpvZhnOB04H3mpif1PuglduQ7PtgMHB6OJxF8AB9qv1faM02JPt+MKBXOJwBLAE+lGL7oTXb0Gb7oSMfQXWE1tRbsw1Jzd1fBfY1s0iy74PWbENSc/ed7v5OOFwGrOHoRpmTej+0chuSWvjZHgxHM8KfxnepJft+aM02tJmOHFBt1pp6ArW2vjPCQ+6/mtmk9imtzST7PmitlNgHZjYSOI3gL99oKbMfmtkGSPL9YGbpZrYc2AO84O4ptx9asQ3QRvuhIwdUm7WmnkCtqe8dgnaspgA/A/4Y76LaWLLvg9ZIiX1gZr2AJ4Evu3tp49kxXpJ0+6GFbUj6/eDute4+laBh7BlmNrnRIkm/H1qxDW22HzpyQHWE1tRbrM/dS+sPuT1oWirDzPq3X4knLNn3QYtSYR+YWQbBF/tv3P2pGIsk/X5oaRtSYT/Uc/cDwGLg0kazkn4/1GtqG9pyP3TkgFoKjDWzUWaWCVwPPNNomWeAT4R3znwIKHH3ne1daDNa3AYzG2RmFg7PINinxe1e6fFL9n3QomTfB2FtvwTWuPtPmlgsqfdDa7YhBfZDrpn1CYe7AxcB7zdaLNn3Q4vb0Jb7Ia6NxSaSx6k19fbUym24FvicmdUAh4DrPbyVJhmY2e8I7urpb2aFwL8SXFhNiX0ArdqGpN4HwFnAPODd8NoBwHeB4ZAy+6E125Ds+2Ew8JiZpRN8af/e3Z9Npe8kWrcNbbYf1NSRiIgkpY58ik9ERFKYAkpERJKSAkpERJKSAkpERJKSAkpERJKSAkqSkpn9S9ha8sqwReSZbbTeKy1Gq/CtfO1dZvb1GNPnmNnEqPF/M7OLTqTOJt7fzOwlM8sOxw+29Jpm1jXfghazPfohyvA9jqk1bTO72czuP95aYqzvFDN7tK3WJ6mrwz4HJanLzM4APkrQenVl+AWaeQyv7+LuNbHmufszHP3A9omaAzwLrA7fI17dPFwOrIjRxM/x+AdBzYsbTb8MGBv+zAT+J/zdbtz9XTMbambD3X1re763JBcdQUkyGgzsdfdKAHff6+47AMxsc/1f/GaWb2aLw+G7zOwhM/sb8CszW2JRjVSa2WIzm1b/176Z9Q7XlRbO72Fm28wsw8w+Y2ZLLWjs8kkz69FUoWZ2JnAlcG94pDfGzB41s2uj6v2Bmf3TzArM7HQze97MNtQ/3Bgu943wPVdajD52QnOBP8WowczsXjN7z8zeNbPrwulpZvZAeCT6rJktqq/L3Ze5++YY79Gq1rTN7BYzW2dmrxA8RFs//Yrws19mZi+a2cCwjg/MLDeqrvVm1t/MPhbWvcLMXo16iz8TtJwinZgCSpLR34Bh4RfgA2Z2XitfNw2Y7e43EnRN8nGA8At2iLu/Xb+gu5cQ9K9Vv+4rgOfdvRp4yt2nh41drgE+3dQbuvsbBEdk33D3qe6+IcZi29z9DOA14FGCJ+0/BPxbWN+HCY5YZgBTgWlmdm6M9ZwFvB1j+tXh66YQND1zb7jNVwMjgVOAW4EzmtqOKC22ph2u++6wnouBiVGzXyfoH+g0gn3wTXevAx4nCFjCGle4+17g+8Al4Wd9ZdR6CoBzWlGvdGAKKEk6YUOT04DbgCJgoZnd3IqXPuPuh8Lh3wMfC4c/DjwRY/mFwHXh8PXhOMBkM3vNzN4l+FI90W4b6k8pvkvQAV2ZuxcBhy1o1+zD4c8ygpagxxMEVmP9wr6QGjsb+F3YyvRu4BVgejj9CXevc/ddwMutqLU1rWnPBBa7e5EH/ZQtjJo3FHg+/Oy+wZHP7hHgE+Hwp4D/DYf/ATxqZp8haM6r3h5gSCvqlQ5MASVJKfyyXezu/wrMB64JZ9Vw5N9tt0YvK496/Xag2MxOJQihBTHe5hngMjPrRxCIL4XTHwXmu/spBEcKjd/nWFWGv+uihuvHuxCEwg/DI7Cp7n6Su/8yxnpq6k9JNhIrVJqb3pzWtqbdVBtpPwPuDz+7zxJ+du6+DdhtZhcQBNxfw+m3A98L33O5meWE6+lG0I6bdGIKKEk6ZjbOzKKPIKYCW8LhzQRhAkdCqykLgG8Cvd393cYzwyO1t4D/JuiiujaclQXstKB7h7mNXxdDWfia4/U88CkL+jrCzPLMbECM5dYCo2NMfxW4zoKO5HIJuqh/i+B02zXhNZ+BBA3etqQ1rWkvAWaZWU74GX0sal5vYHs4/MlGr3uY4FTf7+s/azMb4+5LwhtL9nIkHE8G3mtFvdKBKaAkGfUiaDF5tZmtJLjGcVc4727gv83sNaC2idfX+wPBqbvfN7PMQuAmGp6mupPgS/gFju4OIZYFwDfCGwPGtGL5Btz9b8BvgX+Gp8b+QOzA+wuxQ+ZpYCXBNbWXCK777CLoO6mQ4Iv+/xFsUwmAmX3RgpbZhwIrzezhcF2LgI0ErWn/ArgjRr07CfbHP4EXCU5L1rsLeCLcP3sbvfQZgn37v1HT7g1v7HiPIGhXhNPPD7dXOjG1Zi6SIsKbE37l7hcfw2t6ufvB8NTZW8BZYXi1OzPLB+5z92ZvfjCzrgTX0c5u6nEB6Rz0HJRIinD3nWb2CzPLPoZnoZ4Nb8TIBO5JYDh9G/gcrTtlOhz4tsJJdAQlIiJJSdegREQkKSmgREQkKSmgREQkKSmgREQkKSmgREQkKf1/7SXtV0+8s2kAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "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": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(2.0641865669356645, 0.8499523545190062)" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mu = log_duration.mean()\n", "sigma = log_duration.std()\n", "mu, sigma" ] }, { "cell_type": "code", "execution_count": 35, "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": 55, "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": 64, "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": 65, "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": 78, "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": 78, "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": 79, "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": 79, "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": 68, "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": 68, "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": 69, "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": 69, "metadata": {}, "output_type": "execute_result" } ], "source": [ "likelihood = densities.prod(axis=2)\n", "likelihood.shape" ] }, { "cell_type": "code", "execution_count": 70, "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": 70, "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": 80, "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": 80, "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": 81, "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": 81, "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": 82, "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": 82, "metadata": {}, "output_type": "execute_result" } ], "source": [ "likelihood = densities.prod(axis=2)\n", "likelihood.shape" ] }, { "cell_type": "code", "execution_count": 83, "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": 83, "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": 84, "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": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtFUlEQVR4nO3deZxcZZ3v8c+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": [ "
" ] }, "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": 75, "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": 85, "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": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA52ElEQVR4nO3deZwcZ33v+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": [ "
" ] }, "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": 86, "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": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5V0lEQVR4nO3de3xcd3nv+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": [ "
" ] }, "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": "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.8.12" } }, "nbformat": 4, "nbformat_minor": 4 }