soln/chap05.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimating Counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "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/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:19.989007Z",
     "iopub.status.busy": "2021-04-16T19:35:19.988518Z",
     "iopub.status.idle": "2021-04-16T19:35:19.991079Z",
     "shell.execute_reply": "2021-04-16T19:35:19.990545Z"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "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:35:19.994541Z",
     "iopub.status.busy": "2021-04-16T19:35:19.994130Z",
     "iopub.status.idle": "2021-04-16T19:35:19.995904Z",
     "shell.execute_reply": "2021-04-16T19:35:19.996257Z"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "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:35:19.998961Z",
     "iopub.status.busy": "2021-04-16T19:35:19.998546Z",
     "iopub.status.idle": "2021-04-16T19:35:20.677236Z",
     "shell.execute_reply": "2021-04-16T19:35:20.677597Z"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "from utils import set_pyplot_params\n",
    "set_pyplot_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous chapter we solved problems that involve estimating proportions.\n",
    "In the Euro problem, we estimated the probability that a coin lands heads up, and in the exercises, you estimated a batting average, the fraction of people who cheat on their taxes, and the chance of shooting down an invading alien.\n",
    "\n",
    "Clearly, some of these problems are more realistic than others, and some are more useful than others.\n",
    "\n",
    "In this chapter, we'll work on problems related to counting, or estimating the size of a population.\n",
    "Again, some of the examples will seem silly, but some of them, like the German Tank problem, have real applications, sometimes in life and death situations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Train Problem\n",
    "\n",
    "I found the train problem \n",
    "in Frederick Mosteller's, [*Fifty Challenging Problems in\n",
    "  Probability with Solutions*](https://store.doverpublications.com/0486653552.html):\n",
    "\n",
    "> \"A railroad numbers its locomotives in order 1..N.  One day you see a locomotive with the number 60.  Estimate how many locomotives the railroad has.\"\n",
    "\n",
    "Based on this observation, we know the railroad has 60 or more\n",
    "locomotives.  But how many more?  To apply Bayesian reasoning, we\n",
    "can break this problem into two steps:\n",
    "\n",
    "* What did we know about $N$ before we saw the data?\n",
    "\n",
    "* For any given value of $N$, what is the likelihood of seeing the data (a locomotive with number 60)?\n",
    "\n",
    "The answer to the first question is the prior.  The answer to the\n",
    "second is the likelihood.\n",
    "\n",
    "We don't have much basis to choose a prior, so we'll start with\n",
    "something simple and then consider alternatives.\n",
    "Let's assume that $N$ is equally likely to be any value from 1 to 1000.\n",
    "\n",
    "Here's the prior distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:20.681706Z",
     "iopub.status.busy": "2021-04-16T19:35:20.681246Z",
     "iopub.status.idle": "2021-04-16T19:35:20.682954Z",
     "shell.execute_reply": "2021-04-16T19:35:20.683296Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from empiricaldist import Pmf\n",
    "\n",
    "hypos = np.arange(1, 1001)\n",
    "prior = Pmf(1, hypos)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now let's figure out the likelihood of the data.\n",
    "In a hypothetical fleet of $N$ locomotives, what is the probability that we would see number 60?\n",
    "If we assume that we are equally likely to see any locomotive, the chance of seeing any particular one is $1/N$.\n",
    "\n",
    "Here's the function that does the update:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:20.687062Z",
     "iopub.status.busy": "2021-04-16T19:35:20.686637Z",
     "iopub.status.idle": "2021-04-16T19:35:20.688814Z",
     "shell.execute_reply": "2021-04-16T19:35:20.688368Z"
    }
   },
   "outputs": [],
   "source": [
    "def update_train(pmf, data):\n",
    "    \"\"\"Update pmf based on new data.\"\"\"\n",
    "    hypos = pmf.qs\n",
    "    likelihood = 1 / hypos\n",
    "    impossible = (data > hypos)\n",
    "    likelihood[impossible] = 0\n",
    "    pmf *= likelihood\n",
    "    pmf.normalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function might look familiar; it is the same as the update function for the dice problem in the previous chapter.\n",
    "In terms of likelihood, the train problem is the same as the dice problem.\n",
    "\n",
    "Here's the update:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:20.692485Z",
     "iopub.status.busy": "2021-04-16T19:35:20.692048Z",
     "iopub.status.idle": "2021-04-16T19:35:20.698353Z",
     "shell.execute_reply": "2021-04-16T19:35:20.698751Z"
    }
   },
   "outputs": [],
   "source": [
    "data = 60\n",
    "posterior = prior.copy()\n",
    "update_train(posterior, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the posterior looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:20.702138Z",
     "iopub.status.busy": "2021-04-16T19:35:20.701696Z",
     "iopub.status.idle": "2021-04-16T19:35:20.889486Z",
     "shell.execute_reply": "2021-04-16T19:35:20.889838Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA21klEQVR4nO3deZwV1Z338c+PXukNmn3VBkSQzUYbJBIN7qIxaKKJGkfNJHFM4kzyzBNniHk0OhMnTDQTR+Oo6BhN4hpNlBgS0SiuiKAism8iNPvaLL03v+ePqr7cbnq5DX373u7+vl+v++q6VXWqThXQX07VqTrm7oiIiCSbLomugIiISEMUUCIikpQUUCIikpQUUCIikpQUUCIikpQUUCIikpQUUCKtyMy+bmZz4rj9x8zsp+H0GWa2shW3/Rczuy6cvt7M3m7Fbcf1vEjHpICSdsvM1ptZmZkdMLNtZvZrM8s5hu3dbma/O5Y6ufsT7n7+sWyjBft6y91HNLderMfl7lPd/fFjrZeZFZiZm1lq1Lbb7LxIx6GAkvbuEnfPAU4BJgD/L1EVif6FfBRlzcwS8u8xkfsWaYr+UkqH4O6bgL8AYwDM7EtmttTM9prZXDM7qXZdM/tXM9tkZvvNbKWZnWNmFwK3AF8LW2Qfh+t2M7P/NbMtYZmfmllKuOx6M3vHzH5pZruB2+tfGjOz081sgZmVhD9Pj1o218zuNLN3gFJgaP3jMrPxZvZhWNdngMyoZVPMrPgoj+uIfYfzvlV393ZfWPcVZnZO1IL1ZnZu1PfoVtqb4c+94T4/dxTn5d/Dc7vfzOaYWa9G/uilA1NASYdgZoOBi4CPzOxE4CngB0BvYDbwJzNLN7MRwE3ABHfPBS4A1rv7X4H/AJ5x9xx3Pznc9ONANXACMB44H4j+JX4asA7oA9xZr049gD8D9wI9gf8C/mxmPaNW+zvgBiAX+Kxe+XTgBeC3QA/g98BXGjn+lh5Xk/uud2y9gJ8AfwiPqTlnhj+7h/ucV6+usZyXq4FvEJzXdOCHMexXOhgFlLR3L5jZXuBt4A2CX8ZfA/7s7q+4exVwN9AVOB2oATKAUWaW5u7r3X1tQxs2s77AVOAH7n7Q3bcDvwSujFpts7vf5+7V7l5WbxMXA6vd/bfh8qeAFcAlUes85u5Lw+VV9cpPAtKAe9y9yt2fAxY0ch5iPq4Y9w2wPWrfzwArw2M6VrGcl1+7+6rwnD4LFLbCfqWdUUBJe3epu3d39+Pd/bvhL7QBRLUI3P0QsBEY6O5rCFpWtwPbzexpMxvQyLaPJwiILeGlwr3AQwT/q6+1sYm61alH6DNgYAvKb/K6b3RuqKVDC48rln3TyL6b22YsYjkvW6OmS4Gj7vwi7ZcCSjqizQThAgQ3UoDBwCYAd3/S3T8fruPAf4ar1n+1/0agAugVhmB3d89z99FR6zQ1HECdeoSOq61HDOW3AAPD+keXb1ALjiuWfdPIvjeH0weBrKhl/Vqw3VjOi4gCSjqkZ4GLw04CacD/JQiad81shJmdbWYZQDlQRnB5DGAbUFDbo83dtwBzgF+YWZ6ZdTGzYWb2hRjrMRs40cyuNrNUM/saMAp4Kcby8wjuf/1TWP7LwMSGVmzJcbVAn3DfaWZ2BXBSeEwAi4Arw2VFwOVR5XYAh2ig00foWM+LdBIKKOlw3H0lcA1wH7CT4N7GJe5eSXCfZkY4fyvBL+FbwqK/D3/uMrMPw+lrCW7SLwP2AM8B/WOsxy7giwQBuQv4F+CL7r4zxvKVwJeB68N9fw34QyOrt/S4YjEfGB5u807g8vCYAG4FhoX1ugN4MqrepeH674SXRifVO65jOi/SeZgGLBQRkWSkFpSIiCQlBZSIiCQlBZSIiCQlBZSIiCSlo365ZXvSq1cvLygoSHQ1RESkAR988MFOd+9df36nCKiCggIWLlyY6GqIiEgDzKzBN6ToEp+IiCQlBZSIiCSluAaUmV0YjkuzxsymN7DczOzecPliMzsllrJm9o/hsqVm9vN4HoOIiCRG3O5BWTCo2/3AeUAxsMDMZrn7sqjVphK8SmU4wdgzDwCnNVXWzM4CpgHj3L3CzKLfLC0iraSqqori4mLKy8sTXRXpIDIzMxk0aBBpaWkxrR/PThITgTXuvg7AzJ4mCJbogJoG/CZ8pf97ZtbdzPoDBU2U/Q4ww90rAMIxekSklRUXF5Obm0tBQQF1X2ou0nLuzq5duyguLmbIkCExlYnnJb6B1B1vppi64700tU5TZU8EzjCz+Wb2hplNaGjnZnaDmS00s4U7duw4hsMQ6ZzKy8vp2bOnwklahZnRs2fPFrXI4xlQDf2trv9m2sbWaapsKpBPMNrozcCz1sC/IHef6e5F7l7Uu/cR3esTpqK8oYFLRZKTwklaU0v/PsUzoIoJBomrNYjDg501t05TZYuBP3jgfYJxZ3q1Yr3j5oO/reXfr3qOmT96Bb1FXkSkafEMqAXAcDMbYmbpwJXArHrrzAKuDXvzTQJKwkHimir7AnA2gJmdSDBWT7sYR+a5e+ZRU13Dp0u2sfitBp9LE5EoKSkpFBYWMmbMGK644gpKS0tbVH79+vU8+eSTza/YgNNPP/2oyjXnrbfeYvTo0RQWFjJv3jxmz57dfKEm7N27l//5n/85qrIXXXQRe/fubVGZZ599llGjRjF69GiuvvrqyPzHH3+c4cOHM3z4cB5//PGjqk99cQsod68GbgJeBpYDz7r7UjO70cxuDFebDawD1gAPA99tqmxY5lFgqJktAZ4GrvN22BzZvHZ3oqsgkvS6du3KokWLWLJkCenp6Tz44IMtKn80AVVTEwxE/O6777a4TCyeeOIJfvjDH7Jo0SJWrlzZ4oCqrq6u872pgGquXrNnz6Z79+4x73v16tX87Gc/45133mHp0qXcc889AOzevZs77riD+fPn8/7773PHHXewZ8+emLfbmLg+B+Xus939RHcf5u53hvMedPcHw2l39++Fy8e6+8KmyobzK939Gncf4+6nuPtr8TyGeCnZ1bL/CYp0dmeccQZr1qxh9+7dXHrppYwbN45JkyaxePFiAN544w0KCwspLCxk/Pjx7N+/n+nTp/PWW29RWFjIL3/5S2pqarj55puZMGEC48aN46GHHgJg7ty5nHXWWVx99dWMHTsWgJycHCDofXbzzTczZswYxo4dyzPPPNNomWjf+c53KCoqYvTo0fzkJz8B4JFHHuHZZ5/l3/7t37jqqqu47bbbeOaZZygsLOSZZ57h4MGD/P3f/z0TJkxg/PjxvPjiiwA89thjXHHFFVxyySWcf/75dfYzffp01q5dS2FhITfffHOD9br00ks59dRTGT16NDNnzoyULSgoYOfOnaxfv56TTjqJb3/724wePZrzzz+fsrKyI47p4Ycf5nvf+x75+fkA9OkTPOXz8ssvc95559GjRw/y8/M577zz+Otf/3o0f8x1dIp38SWjfbuO/MMXSVY/uuR3cdv2z/50TbPrVFdX85e//IULL7yQn/zkJ4wfP54XXniB1157jWuvvZZFixZx9913c//99zN58mQOHDhAZmYmM2bM4O677+all14CYObMmXTr1o0FCxZQUVHB5MmTI7/w33//fZYsWXJEF+g//OEPLFq0iI8//pidO3cyYcIEzjzzzCbLANx555306NGDmpoazjnnHBYvXsy3vvUt3n77bb74xS9y+eWX89hjj7Fw4UJ+9atfAXDLLbdw9tln8+ijj7J3714mTpzIueeeC8C8efNYvHgxPXr0qLOfGTNmsGTJEhYtWgQEwVm/Xo8++ig9evSgrKyMCRMm8JWvfIWePXvW2c7q1at56qmnePjhh/nqV7/K888/zzXX1P2zWbVqFQCTJ0+mpqaG22+/nQsvvJBNmzYxePDhbgODBg1i06ZNzf65NkcBlSD71IISaVZZWRmFhYVA0IL65je/yWmnncbzzz8PwNlnn82uXbsoKSlh8uTJ/PM//zNf//rX+fKXv8ygQYOO2N6cOXNYvHgxzz33HAAlJSWsXr2a9PR0Jk6c2GDQvP3221x11VWkpKTQt29fvvCFL7BgwQLy8vIaLQPBvZqZM2dSXV3Nli1bWLZsGePGjWvyeOfMmcOsWbO4++67gaCr/4YNGwAiLZRY1K/Xvffeyx//+EcANm7cyOrVq48IqCFDhkTO9amnnsr69euP2G51dTWrV69m7ty5FBcXc8YZZ7BkyZIGO321Rg9QBVSCqAUl0rzae1DRGvtlOH36dC6++GJmz57NpEmTePXVV49Yz9257777uOCCC+rMnzt3LtnZ2Q3Woalb3I2V+fTTT7n77rtZsGAB+fn5XH/99TE9/+PuPP/884wYMaLO/Pnz5ze6r+bqNXfuXF599VXmzZtHVlYWU6ZMabAuGRkZkemUlJQGL/ENGjSISZMmkZaWxpAhQxgxYgSrV69m0KBBzJ07N7JecXExU6ZMibm+jVFAtaG0jFSqKoIbnFWV1c2sLZI8YrkM11bOPPNMnnjiCW699Vbmzp1Lr169yMvLY+3atYwdO5axY8cyb948VqxYweDBg9m/f3+k7AUXXMADDzzA2WefTVpaGqtWrWLgwPrvDzhyfw899BDXXXcdu3fv5s033+Suu+5ixYoVjZbZt28f2dnZdOvWjW3btvGXv/ylwV/Yubm5R9Tvvvvu47777sPM+Oijjxg/fnyT9au/jfpKSkrIz88nKyuLFStW8N577zW5vaZceumlPPXUU1x//fXs3LmTVatWMXToUIYNG8Ytt9wS6RgxZ84cfvaznx31fmopoNpQt55Z7Ny8L/Ld3fUgpEgL3X777XzjG99g3LhxZGVlRbo033PPPbz++uukpKQwatQopk6dSpcuXUhNTeXkk0/m+uuv5/vf/z7r16/nlFNOwd3p3bs3L7zwQpP7u+yyy5g3bx4nn3wyZsbPf/5z+vXr12RAnXzyyYwfP57Ro0czdOhQJk+e3OB6Z511FjNmzKCwsJAf/ehH3HrrrfzgBz9g3LhxuDsFBQWR+2eN6dmzJ5MnT2bMmDFMnTqViy++uM7yCy+8kAcffJBx48YxYsQIJk2a1OT2mnLBBRcwZ84cRo0aRUpKCnfddVfkUuGtt97KhAnBi31uu+22mC9HNsXaYQ/tFisqKvJkGLDwvu/PZvO6w93Lf/y7y8nplpnAGok0bvny5Zx00kmJroZ0MA39vTKzD9y9qP66Gg+qDdX/z4A6SoiINE4BlUDqKCEi0jgFVALpYV1Jdp3hFoC0nZb+fVJAtaH6fzYlOxVQkrwyMzPZtWuXQkpaRe14UJmZsd93Vy++BNI9KElmgwYNori4GI2nJq2ldkTdWCmgEkgBJcms9mFMkUTRJb4E0iU+EZHGKaDaUP1r+Xt3lOr6vohIIxRQCVRZXkXp/opEV0NEJCkpoBJsz7aDia6CiEhSUkC1pQau5u3ZfqDt6yEi0g4ooBJs73a1oEREGqKAakMNdYjYrUt8IiINUkAl2F5d4hMRaZACKsHUSUJEpGEKqATbs/2gnoUSEWmAAirB9CyUiEjDFFBtKLqh1KXL4VOvy3wiIkdSQCVIft/syLSehRIROVJcA8rMLjSzlWa2xsymN7DczOzecPliMzulubJmdruZbTKzReHnongeQ2uKvtfUvU9UQKkFJSJyhLgFlJmlAPcDU4FRwFVmNqrealOB4eHnBuCBGMv+0t0Lw8/seB1DPPXsnxuZ3rVlfwJrIiKSnOLZgpoIrHH3de5eCTwNTKu3zjTgNx54D+huZv1jLNuu9RoQFVCbFVAiIvXFM6AGAhujvheH82JZp7myN4WXBB81s/yGdm5mN5jZQjNbmIwjgvaMCqidCigRkSPEM6CsgXn1H/hpbJ2myj4ADAMKgS3ALxraubvPdPcidy/q3bt3TBVuS/l9czALDrNk50GqKqsTXCMRkeQSz4AqBgZHfR8EbI5xnUbLuvs2d69x90PAwwSXA9uHqHhOTUshv29O5PuuLerJJyISLZ4BtQAYbmZDzCwduBKYVW+dWcC1YW++SUCJu29pqmx4j6rWZcCSOB5DXNW9D7UvgTUREUk+qfHasLtXm9lNwMtACvCouy81sxvD5Q8Cs4GLgDVAKfCNpsqGm/65mRUStEfWA/8Qr2NobfVfadRzQC58GEzrPpSISF1xCyiAsAv47HrzHoyaduB7sZYN5/9dK1czIczUk09EpCl6k0QC9RyQF5nWs1AiInUpoBKoZ//DnSR0iU9EpC4FVIKYGd37ZEdeGrtvVymVFepqLiJSSwHVhuoP+5SamkKPflGtqE3qySciUksBlWC9Bx++D7V9Y0kCayIiklwUUG2ooZFz+x7XPTK9fYMCSkSklgIqwfqoBSUi0iAFVKKEbxvsE9WC2qYWlIhIhAKqLR15hY/eg/IiL43dveUA1VU1bVwpEZHkpIBKkNpQSs9IjYyue+jQIT0PJSISUkAlgb7HdYtM6z6UiEhAAZUgFjXiVZ/BUQGl+1AiIoACqk011M0coE9UC2rbhr1tVBsRkeSmgEoC0c9CbV2/N2H1EBFJJgqoJND3+G6RThO7Nu+norwqwTUSEUk8BVQbauQKH2npqfQelBeu42z7TPehREQUUAli0b0kgP5D8iPTWz7d09bVERFJOgqoJNF/aFRArVNAiYgooBKlbgOqTgtq87rdbVwZEZHko4BqQ411M4e6Lahtn5Vw6NChtqiSiEjSUkAlidzuXcnN7wpAZXkVu7YcSHCNREQSSwGVROreh9JlPhHp3BRQbanxK3wADBjaIzJdvFoBJSKdmwIqQep3MwcYdGLPyHTx6l1tWR0RkaSjgEoig+sFVE2NOkqISOcV14AyswvNbKWZrTGz6Q0sNzO7N1y+2MxOaUHZH5qZm1mveB5DW8rrkUVezywAqiqqNfSGiHRqcQsoM0sB7gemAqOAq8xsVL3VpgLDw88NwAOxlDWzwcB5wIZ41T/eGrjCB8BxIw7nbfEqXeYTkc4rni2oicAad1/n7pXA08C0eutMA37jgfeA7mbWP4ayvwT+hWa7HSSXpp6DqjVw+OHLfBsVUCLSicUzoAYCG6O+F4fzYlmn0bJm9iVgk7t/3NTOzewGM1toZgt37NhxdEeQANH3oTau3JnAmoiIJFY8A6qhi1j1mxCNrdPgfDPLAn4M3Nbczt19prsXuXtR7969m61sW4ihAcWg4T0jPfy2fbZXQ2+ISKcVz4AqBgZHfR8EbI5xncbmDwOGAB+b2fpw/odm1q9Va94GGupmDpDRNS0yBLy7s2mNnocSkc4pngG1ABhuZkPMLB24EphVb51ZwLVhb75JQIm7b2msrLt/4u593L3A3QsIguwUd98ax+NoNbHcg4K6z0NtWKHLfCLSOcUtoNy9GrgJeBlYDjzr7kvN7EYzuzFcbTawDlgDPAx8t6my8aprsikY3ScyvX7p9gTWREQkcVLjuXF3n00QQtHzHoyaduB7sZZtYJ2CY69lgjTSzRxgSFRAfbZsB4cOHaJLFz1TLSKdi37rJaEe/XLI6xE8sFteWsnW9XsTWyERkQRQQCVIY50kapcNGXO4FfXpEl3mE5HORwHVllrwWLHuQ4lIZ6eASlL1W1Cx9gAUEekoFFBtqCUh02dwN7JyMwA4uK+cHcX74lUtEZGkpIBKUmZW5zLf2sXt4lEvEZFWo4BKkCb6SEScUHj4BRmrP9oSx9qIiCQfBVQSGz6+f2R63eJtVFfXJLA2IiJtSwGVIE11M6/Vs38u+X1yAKgoq2LjSg2/ISKdhwKqDbW0I56ZMfyUw60oXeYTkc5EAZXkoi/zrVFAiUgnooBqQ0fzLNPQcX0jlwOLV++idH9Fa1dLRCQpKaCSXFZORmT4DXfXZT4R6TQUUIkSQzfzWiOKBkaml88vjkNlRESSjwKqHRg1aVBkeuXCzepuLiKdggIqQWLpZl6r3/HdI93Ny0sr9XZzEekUFFBt6Sjf92pmnHTa4VbUsvd0mU9EOj4FVDsRfZlv+fxivd1cRDo8BVQbOpZQOX5UbzKz0wEo2XmQzev2tFa1RESSkgKqnUhNTeGkiYd7833y9mcJrI2ISPwpoBKkBX0kIsadURCZXvzmZ7rMJyIdmgKqDR1rnpwwvl/kMt+e7QfYuEovjxWRjqvJgDKzx6Kmr4t7bTqRlnQzr5WamsKY04+LfF/85vpWrJGISHJprgV1ctT09+NZEYnNuDOPj0x/8vYGDh06lMDaiIjET3MBdUwXpczsQjNbaWZrzGx6A8vNzO4Nly82s1OaK2tm/x6uu8jM5pjZgGOpY3szdGxfsrtlArBvdylrF29LcI1EROKjuYAaFAbIfVHTkU9TBc0sBbgfmAqMAq4ys1H1VpsKDA8/NwAPxFD2Lncf5+6FwEvAbTEea8LV6dRwFJ0kAFJSunDymQWR7x+8svbYKiUikqRSm1l+c9T0whZueyKwxt3XAZjZ08A0YFnUOtOA33jwm/s9M+tuZv2BgsbKuvu+qPLZHGMrrz0qOm8Y7/5pBQBL522k9EAFWTkZCa6ViEjrajKg3P3xY9j2QGBj1Pdi4LQY1hnYXFkzuxO4FigBzjqGOrZL/YfkM2BYDzav3U11VQ0fv7Gez108ItHVEhFpVc314pvV1KeZbTd0Eat+a6exdZos6+4/dvfBwBPATY3U/QYzW2hmC3fs2NFMVdtIK7b1Jpx/QmR6oS7ziUgH1Nwlvs8RtGSeAubTsjsnxcDgqO+DgM0xrpMeQ1mAJ4E/Az+pv8DdZwIzAYqKijrcZcBxZx7Pnx/5gOqqGjav3c3mdbsZMLRHoqslItJqmusk0Q+4BRgD/DdwHrDT3d9w9zeaKbsAGG5mQ8wsHbgSqN/qmgVcG/bmmwSUuPuWpsqa2fCo8l8CVjR7lEnoaJ6DipaVk8Hozx3O8Pf/uuZYqyQiklSaDCh3r3H3v7r7dcAkYA0w18z+sbkNu3s1weW3l4HlwLPuvtTMbjSzG8PVZgPrwu0+DHy3qbJhmRlmtsTMFgPn04mfz5p44eGs/vC1dZQeqEhgbUREWldzl/gwswzgYuAqgt519wJ/iGXj7j6bIISi5z0YNe3A92ItG87/Siz7TnbH2IACYMiYPvQryGfr+j1UVVTzwStrOeOy+j35RUTap+Y6STwOvAucAtzh7hPc/d/dfVOb1K6Dae2Xu5oZp19yuPfevJdW6c0SItJhNHcP6u+AEwkuo80zs33hZ7+Z7WumrLSBwikFZOUGz0Dt2X6A5e/r/w4i0jE0dw+qi7vnRn3ywk+uu+e1VSU7iniMjpGWnlrnXtQ7L7bLPiMiIkdo7hJfppn9wMx+FT5X1Ow9K2l7p100nC4pwR/lp0u28dnyJHnuS0TkGDR3ie9xoAj4BLgI+EXca9RJHGs382jde2VTOKUg8v21Zz5ptW2LiCRKcwE1yt2vcfeHgMuBM9qgTnIUplwxJhJ6qz7YzKa1uxNcIxGRY9NcQFXVToTPJklrab0GFAC9B+YxZvLhwQznPrukdXcgItLGmh2wMLrnHjBOvfiOXmt3M6/vrK+OiUwveXcDWz7dE9f9iYjEU3O9+FLq9dxLVS++5NV/SD4nnTYo8n3ObxclrjIiIseouRaUtKY2eGXtedecHLkXtWLBJj5duj3+OxURiQMFVAfTvyC/To++lx//KO6XFkVE4kEB1Yaig6I1u5nXd941J5OSmgLAZ8t3sOy9jc2UEBFJPgqoDii/Tw6nTT38dok/P/IhVZXqhCki7YsCKkHi2IAC4JyrxtZ5R9+bf1ge3x2KiLQyBVQHlZWbwfnXFka+v/H7pezZfiBxFRIRaSEFVAc24fxhDBgWDANfVVnNSw9/kOAaiYjETgHVhtqqk0StLl268KUbJ0S+L3tvI0ve3RD3/YqItAYFVAd3/MjeFJ13QuT7iw8soHS/hoYXkeSngOoELvrmKeT1yALgwN4yXnp4YYJrJCLSPAVUJ9A1O53Lbjot8v2j1z9lxQKNvCsiyU0BlSBtcQ8q2sgJAymcMiTy/bn/nse+3aVtWgcRkZZQQLWRZHjd0CU3FEUu9R0sKef3v5yXFPUSEWmIAqoTycrN4Kv/9/RI623Noi288fyyBNdKRKRhCqhOZti4fky5YnTk+yu/+1hvPBeRpKSA6oTOvmosx43sDcChmkM88bM32bvzYIJrJSJSV1wDyswuNLOVZrbGzKY3sNzM7N5w+WIzO6W5smZ2l5mtCNf/o5l1j+cxtJa2fki3KampKVz1L58nu1smENyP+t2db+qFsiKSVOIWUGaWAtwPTAVGAVeZ2ah6q00FhoefG4AHYij7CjDG3ccBq4AfxesYOrLuvbO5evoZdEkJ/gpsWrOLP/5qvjpNiEjSiGcLaiKwxt3XuXsl8DQwrd4604DfeOA9oLuZ9W+qrLvPcffa/+q/BwxCjsrQMX354rdPjXz/6PVPee2ZJQmskYjIYfEMqIFA9Eh5xeG8WNaJpSzA3wN/aWjnZnaDmS00s4U7duxoYdVbX7I2TCZddGKdVyG9+sTHLHxlTQJrJCISiGdANXSjpf6v6cbWabasmf0YqAaeaGjn7j7T3Yvcvah3794xVLdzMjOmfXcCJxT2j8z746/e15smRCTh4hlQxcDgqO+DgM0xrtNkWTO7Dvgi8HVvJzdNkqmTRH2pqSlcc8uZDBgaDM1x6NAhnpzxFus+2ZbgmolIZxbPgFoADDezIWaWDlwJzKq3zizg2rA33ySgxN23NFXWzC4E/hX4krvrXT2tJKNrGtfffhb5fXKAYPyox+54nXVLFFIikhhxC6iwI8NNwMvAcuBZd19qZjea2Y3harOBdcAa4GHgu02VDcv8CsgFXjGzRWb2YLyOIW6SqwEVkZvflW/+9JzI65CqKqp5/I65rF+mB3lFpO1ZO7lCdkyKiop84cLEDjFRXV3DrZc9BUCXlC7c+cLVCa1PU3Zs2scjt7waeZlsemYa1902haFj+ya4ZiLSEZnZB+5eVH++3iQhR+g9MI9v3nkOOd27AlBZXsWvf/KaRuMVkTalgGorUQ3VZOsk0ZA+g7rx7f84N3K5r7qqhidnvMX7L69OcM1EpLNQQEmj+gzuxo13nU+vAXlA0BPxj7+azytPfKw3TohI3Cmg2kh7/YWe3yeHf/j5+Qw8oWdk3mtPf8KTM96iorwqgTUTkY5OASXNyumWybf/41yGjz/8MO+Sdzfw0L/MYc/2AwmsmYh0ZAqoBGgHt6COkNE1jet+chanXzIyMm/Lp3u4/5//yupFWxJYMxHpqBRQErOUlC5cckMRX/7HSaSkpgDBUB2/vu015vx2ETU1hxJcQxHpSBRQ0mITzj+Bb/70cDd0d+f1Z5fwyC2vauBDEWk1Cqg24u2sm3lzhozuwz/dd1Gdl8yuX7ad/77pz3z42rp22ylERJKHAkqOWm73rvz9v53N+X9XGAnd8oOV/P6X7/Lbn77B/j1lCa6hiLRnCqg20lFbFGbGWV8dww0zzqNHv9zI/OXvF3PP915Sa0pEjpoCSlpFwajgkt+ki0dE5pXur+D3v3yXmT96ha2f7U1c5USkXVJAJUL7vwXVoIzMNKbdOIFv3XluZNgOgPVLt3PfP83mz//7ARVlerhXRGKjgJJWN2xcP75//8Wc+ZXRdEkJ/oodOnSIt19Yzi/+YRbz/7paXdJFpFkKKImLjMw0pl4/nu/fd3GdYTr27ynjhfvn8983/Zll8zfq/pSINEoB1VY6WDfzWPUZ3I1v3XkuX/vhZPJ6ZkXm7ygu4bc/fYOH/nUOaz7eqqASkSOkJroC0vGZGYVfGMKoSYN5908reeP3SykvrQTgs+U7+N//9yrHn9Sbc64axwmF/TpVgItI49SCaiNqIUB6RipTLh/NDx/+EpO/NDLyuiQIgurR2/7GAze/rEt/IgKoBSUJkJ2XyRe/XcTkaSOZ+/ulfPDqOmqqawDYuHInv/3pG/QakMfkaSMZf/YQMjLTElxjEUkEtaAkYfL75HDZ907jhw9/iUkXnVinRbVz8z5efOB9/vMbL/DybxbpHX8inZBaUG2k7rv4ElePZNS9VzbTvjORKV8dw7uzVvD+y2soPxjcoyo7UMHc3y/hjeeWMqJoAKdNPZETT+1Ply76v5VIR6eAkqTRrWcWU79xCmdfOZYPXl3LO7NWsnvrfiC4h7diwSZWLNhEt17ZTDh/GKeeN4zuvbITXGsRiRcFVAKol1rTMrqmcfolI5l08Ykse6+Y92avYu3HWyPLS3Ye5NUnF/O3pz5hyJg+jD9rCKNPP46u2ekJrLWItDYFlCStLl26MOb04xhz+nHs3LyPBXPW8sGrazlYUg4Erap1n2xj3SfbePGBBYycMJDCs4Yw4tQBpKalNLN1EUl2Cqg2UqfbtBpQLdZrQB5Trx/PeV8fx7L5xSx4eQ1rox7wra6qYcm7G1jy7gYys9IZMWEAY04/juGn9FcvQJF2Kq4BZWYXAv8NpACPuPuMesstXH4RUApc7+4fNlXWzK4AbgdOAia6+8J4HoMkl9S0FMZ9/njGff549u0uZfFbn/HR65+yee3uyDrlpZV8/MZ6Pn5jPWnpqZwwvh9jTj+OkRMGkpWbkcDai0hLxC2gzCwFuB84DygGFpjZLHdfFrXaVGB4+DkNeAA4rZmyS4AvAw/Fq+7SPuT1yOLz007i89NOYvvGEhbN/ZSP3/ws0rECoKqymuXzi1k+vxgz47iRvRhRNJATTx3AgKH5uh8oksTi2YKaCKxx93UAZvY0MA2IDqhpwG88uE7znpl1N7P+QEFjZd19eTgvjlWPA70YIa76DO7G+X9XyHnXnMzW9XtZ8u4Gls7byLaocajcnc+W7+Cz5TuY89tF5HTvyomn9mfEqQM4obC/WlciSSaeATUQ2Bj1vZigldTcOgNjLNskM7sBuAHguOOOa0lRacfMjP5D8uk/JJ/zvn4yOzbtY+m8jSx7byPFq3bVuRd4YG8ZH/5tHR/+bV2k3JCxfRk2ri8Fo/uoV6BIgsUzoBpq4tRvRzS2Tixlm+TuM4GZAEVFRUnVfml3rb92rPfAPKZcPpopl4/mQEk5qz/awsqFm1j94RZK91dE1nN3Nq/bzeZ1u3nnxeWYGQOG9WDo2L4MHdeX40b2IitHLSyRthTPgCoGBkd9HwRsjnGd9BjKirRITrdMxk8ZwvgpQzh06BDFq3ez6oPNrPpwM5tW7+bQocODKLo7m9bsYtOaXbz1x+CqdJ/B3ThuZG+OP6kXx43sTe9BefrPhkgcxTOgFgDDzWwIsAm4Eri63jqzgJvCe0ynASXuvsXMdsRQtt3S77TE69KlC8eN6MVxI3px7tXjKC+tZP2yHaxbHDxXtXnt7iPeqL59YwnbN5aw8JU1AGRmp3PcyCCsBg3vycATepDTLTMRhyPSIcUtoNy92sxuAl4m6Cr+qLsvNbMbw+UPArMJupivIehm/o2mygKY2WXAfUBv4M9mtsjdL4jXcbQWDR+R3DKz0hlZNJCRRQMBKDtYyadLtrHuk+18tmw7m9fuqdPCAig/WBm0wD443Ljv1iubAcPyGXhCTwYO68HAE3qQm9+1TY9FpKOwzvCLs6ioyBcuTOzjUgf3lfPTrz8HQFZuBrc+eUVC6yMtU1lRTfHqXWxYsZMNK3awYcXOyBstmpPXI4sBw4KOG32P707f47vTa2Auqal624UIgJl94O5F9efrTRJtpBP8P6BDS89IZeiYvgwd0xcIWsS7tx4IwmplcK9q66d7qaqsPqLsvt2l7NtdyooFmyLzUlJT6DUwl35hYPUrCH7m98nWfS2RkAJK5CiYGT3759Kzfy7jzxoKQE3NIbZvLGHz2t1sXruHTWt2sXndHqoqjgytmuoatn22t85zWgDpmWn0GpBL70F59BqYF/nZa2CuXtkknY4CKgH0P+SOKSWlC/0L8ulfkM+p5wTzDh06xM5N+9m0dnckkLZ9VsKe7Qca3EZleVWku3t9eT2z6F0ntPLo0S+H/D7ZejmudEgKKJE46tKlC30Gd6PP4G515peXVrJtQwlb1+9l+4YStq7fw9b1e+s8m1Xfvl2l7NtVytrFW+vMNzPyembRo19OvU8uPfrlkJ2Xof8USbukgBJJgMysdI4f2ZvjR/aOzHN3Du6rYOemfezctI8dxfvYEf7cvfUAh2oONbgtd6dk50FKdh7k0yXbjlienpkWaWl165VFt97ZkenuvbPJ7dGVlBSNUCzJRwHVRjTchjTHzMjplklOt0wKRvWps6y6uoY92w7WCa/dWw+we+sB9u442ORjDJXlVWELbU+j+83rmUX33lmR0OrWK4vufbLJ65lFXo+u5HTPpEsXhZi0LQWUSDuQmpoS3H8amHfEsurqGkp2lIaBtT8SXLWf8tLKJrcd3QJrjJmRm9+V3B5dyc3PJK9nVuR7Xo/DPxVk0poUUG1F3cwlTlJTUyI9CqF/nWXuTtmBykhLK/iUUhJOl+wsZf+esmb34e6R7vJNqQ2ynPxMcvO7kt0tg5zuXcnplkFO90yyu2VGfmZ3y9CzYNIkBZRIB2ZmZOVmkJWbwaDhPRtcp7qqhpKdpZTsLI2EWO30vl1l7NtV2mTnjWixBlmtzOx0croHlzWjg6t2Oqh7Oll5wTGkZ+hXVmeiP+02En2PQD2qJJmkpkW3wBpWXVXD/j1l7N9Tzr5dQatr/+4ySnaVcmBPOft2B98P7ovt7Rq1yg9WUn6wkp2b9sW0flp6Kl1zagMrnazcDLqGIZadd3i6NpSz8jLompOuTiDtlAJKRJqVmpZCfp8c8vvkNLledXVNJLAO7C3jYEkFB/aWc7CknAN7yzlQUl5nXktftVZVWU3V7uqYW2i1MrqmkZmdTtecdDKz0+iaHf7MySAzO1yWHb1O3fV0Xy0xFFAi0mpSU1Po3jub7r2zm13X3SndXxEJrANhiNWGWem+Cg7uq6DsQCWl+yoo3V9BdVXNUdWroqyKirKqJjuCNKV+wGVmpZGRlRbMz0ojPfyZ0TU1Mj8jq3Ze7bqpCroWUkAlgK7wiQSXurPzMsnOyzziQeaGuDuVFdWU7a+kdH8FpbU/w/Aq3V9J2f5g+uC+ish6ZQcqj3k0gWMNuFppGalRoZZaL8AOB1l6ZrBeemYqaRkppGcG0/U/aekpHfqWgQJKRNoFMyMjM42MzLSYWmi13J3y0qrwflcVZeF9r7IDwffyg5WUHayi7EBF5Hv5wSrKDlRSdrCSitKqVhsup6qimqqKag7sbb7nZCzMjLSMI4OrNtgOh1wqGV1TG103PSOV1PSUIAwzgrJpGakJv3engGojdZ7T7cD/4xFJNmZG1+x0umanH1X5xgKuorSKirJqykurqCyrna6kMpxXUVYVWaci/N7awxu5O5XlVVSWV7XqdmulpKaQFgZXWhhctQGWmn54uv+QfM64bFSr718BJSLShGMNuFq1lyjrh1Z0kJWXBtOV5dXhJ3r6yM/R3pOLVU11DTXVNZQ30yflYEmFAqo96wwDQ4pI46IvUbaWmppDVFU0HmC1n6qKairKoqerqKqoCZZXVFNVURNefqyp8z3W31up6fF54FoBJSLSTqWkdCElK53MrGNr3TXE3ampPhQJrOrKmkbCrJq8nlmtvn9QQImISAPMjNS0lGCssaYff4sbdcpPBPWREBFplgJKRESSkgKqraibuYhIiyigREQkKSmg2oi6mYuItIwCSkREklJcA8rMLjSzlWa2xsymN7DczOzecPliMzulubJm1sPMXjGz1eHP/Hgeg4iIJEbcnoMysxTgfuA8oBhYYGaz3H1Z1GpTgeHh5zTgAeC0ZspOB/7m7jPC4JoO/Gs8jqFkVykvzVzYKtuqrKiOTKuPhIhI8+L5oO5EYI27rwMws6eBaUB0QE0DfuPBDZr3zKy7mfUHCpooOw2YEpZ/HJhLnAKqoqyKJe9uaPXtqhefiEjz4nmJbyCwMep7cTgvlnWaKtvX3bcAhD/7NLRzM7vBzBaa2cIdO3Yc9UHEw4mnDkh0FUREkl48W1ANNRPqd2VrbJ1YyjbJ3WcCMwGKioqOqgtdXo+uXD39jKMp2qjsvEyGjGkwU0VEJEo8A6oYGBz1fRCwOcZ10psou83M+rv7lvBy4PZWrXWUzKx0xk4+Pl6bFxGRJsTzEt8CYLiZDTGzdOBKYFa9dWYB14a9+SYBJeFlu6bKzgKuC6evA16M4zGIiEiCxK0F5e7VZnYT8DKQAjzq7kvN7MZw+YPAbOAiYA1QCnyjqbLhpmcAz5rZN4ENwBXxOgYREUkc6wxvOCgqKvKFC1unu7iIiLQuM/vA3Yvqz9ebJEREJCkpoEREJCkpoEREJCkpoEREJCl1ik4SZrYD+OwYNtEL2NlK1WnPdB50DmrpPAR0HlrnHBzv7r3rz+wUAXWszGxhQz1MOhudB52DWjoPAZ2H+J4DXeITEZGkpIASEZGkpICKzcxEVyBJ6DzoHNTSeQjoPMTxHOgelIiIJCW1oEREJCkpoEREJCkpoJpgZhea2UozW2Nm0xNdn3gys8Fm9rqZLTezpWb2/XB+DzN7xcxWhz/zo8r8KDw3K83sgsTVvnWZWYqZfWRmL4XfO+M56G5mz5nZivDvxOc66Xn4P+G/hyVm9pSZZXaG82Bmj5rZdjNbEjWvxcdtZqea2SfhsnvNrKHBaBvn7vo08CEY5mMtMJRgAMWPgVGJrlccj7c/cEo4nQusAkYBPwemh/OnA/8ZTo8Kz0kGMCQ8VymJPo5WOhf/DDwJvBR+74zn4HHgW+F0OtC9s50HYCDwKdA1/P4scH1nOA/AmcApwJKoeS0+buB94HMEo6T/BZjaknqoBdW4icAad1/n7pXA08C0BNcpbtx9i7t/GE7vB5YT/AOdRvDLivDnpeH0NOBpd69w908JxvSa2KaVjgMzGwRcDDwSNbuznYM8gl9Q/wvg7pXuvpdOdh5CqUBXM0sFsghG9u7w58Hd3wR215vdouMORzzPc/d5HqTVb6LKxEQB1biBwMao78XhvA7PzAqA8cB8oK8HoxwT/uwTrtZRz889wL8Ah6LmdbZzMBTYAfw6vNT5iJll08nOg7tvAu4mGBh1C8GI33PoZOchSkuPe2A4XX9+zBRQjWvoWmmH75NvZjnA88AP3H1fU6s2MK9dnx8z+yKw3d0/iLVIA/Pa9TkIpRJc3nnA3ccDBwku6TSmQ56H8B7LNILLVgOAbDO7pqkiDcxr9+chBo0d9zGfDwVU44qBwVHfBxE07zssM0sjCKcn3P0P4extYVOd8Of2cH5HPD+TgS+Z2XqCS7pnm9nv6FznAILjKnb3+eH35wgCq7Odh3OBT919h7tXAX8ATqfznYdaLT3u4nC6/vyYKaAatwAYbmZDzCwduBKYleA6xU3Yu+Z/geXu/l9Ri2YB14XT1wEvRs2/0swyzGwIMJzghmi75e4/cvdB7l5A8Of9mrtfQyc6BwDuvhXYaGYjwlnnAMvoZOeB4NLeJDPLCv99nENwb7aznYdaLTru8DLgfjObFJ6/a6PKxCbRvUWS+QNcRNCbbS3w40TXJ87H+nmC5vdiYFH4uQjoCfwNWB3+7BFV5sfhuVlJC3vnJPsHmMLhXnyd7hwAhcDC8O/DC0B+Jz0PdwArgCXAbwl6qnX48wA8RXDfrYqgJfTNozluoCg8d2uBXxG+vSjWj151JCIiSUmX+EREJCkpoEREJCkpoEREJCkpoEREJCkpoEREJCkpoESimJmb2S+ivv/QzG5vpW0/ZmaXt8a2mtnPFeEbyF+vN7/AzK4+ym2+2zq1E4mdAkqkrgrgy2bWK9EViWZmKS1Y/ZvAd939rHrzC4AGAyp8GWqj3P30FuxfpFUooETqqgZmAv+n/oL6LSAzOxD+nGJmb5jZs2a2ysxmmNnXzez9cCycYVGbOdfM3grX+2JYPsXM7jKzBWa22Mz+IWq7r5vZk8AnDdTnqnD7S8zsP8N5txE8dP2gmd1Vr8gM4AwzWxSOc3S9mf3ezP4EzDGzHDP7m5l9GG53WtS+oo91rh0eK+qJ2jF+wuNeFh7D3S0+8yL1NPm/JpFO6n5gsZn9vAVlTgZOIhiiYB3wiLtPtGDgx38EfhCuVwB8ARgGvG5mJxC8AqbE3SeYWQbwjpnNCdefCIzxYBiDCDMbAPwncCqwhyBgLnX3fzOzs4EfuvvCenWcHs6vDcbrCcbqGefuu8NW1GXuvi9sQb5nZrP8yKf5xwOjCd6r9g4w2cyWAZcBI93dzax7C86dSIPUghKpx4O3uP8G+KcWFFvgwZhaFQSvdakNmE8IQqnWs+5+yN1XEwTZSOB84FozW0QwxElPgveZQfBOszrhFJoAzPXgRabVwBMEYzi11CvuXjvujwH/YWaLgVcJhkbo20CZ99292N0PEbwSqwDYB5QDj5jZl4HSo6iLSB0KKJGG3UNwLyc7al414b+Z8LJWetSyiqjpQ1HfD1H3SkX91kjtsAT/6O6F4WeIB+MOQTDURUNaNnR246K3/3WgN3CquxcC24DMBspEH2sNkBqG5ESCt+FfCvy1leonnZgCSqQBYaviWYKQqrWe4JIaBOMEpR3Fpq8wsy7hfamhBC/XfBn4TjjcCWZ2ogUDBDZlPvAFM+sVdqC4CnijmTL7gdwmlncjGA+ryszOAo6P4XgI65wDdHP32QSXMwtjLSvSGN2DEmncL4Cbor4/DLxoZu8TvM25sdZNU1YSBElf4EZ3LzezRwguk30Ytsx20MzQ2O6+xcx+BLxO0Jqa7e7NDWWwGKg2s4+BxwjuXUV7AviTmS0kuHS3IvbDIpfg3GSG9Tmik4lIS+lt5iIikpR0iU9ERJKSAkpERJKSAkpERJKSAkpERJKSAkpERJKSAkpERJKSAkpERJLS/wdrsGppCB6+egAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from utils import decorate\n",
    "\n",
    "posterior.plot(label='Posterior after train 60', color='C4')\n",
    "decorate(xlabel='Number of trains',\n",
    "         ylabel='PMF',\n",
    "         title='Posterior distribution')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not surprisingly, all values of $N$ below 60 have been eliminated.\n",
    "\n",
    "The most likely value, if you had to guess, is 60."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:20.894441Z",
     "iopub.status.busy": "2021-04-16T19:35:20.893628Z",
     "iopub.status.idle": "2021-04-16T19:35:20.896691Z",
     "shell.execute_reply": "2021-04-16T19:35:20.896316Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "60"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "posterior.max_prob()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That might not seem like a very good guess; after all, what are the chances that you just happened to see the train with the highest number?\n",
    "Nevertheless, if you want to maximize the chance of getting\n",
    "the answer exactly right, you should guess 60.\n",
    "\n",
    "But maybe that's not the right goal.\n",
    "An alternative is to compute the mean of the posterior distribution.\n",
    "Given a set of possible quantities, $q_i$, and their probabilities, $p_i$, the mean of the distribution is:\n",
    "\n",
    "$$\\mathrm{mean} = \\sum_i p_i q_i$$\n",
    "\n",
    "Which we can compute like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:20.900289Z",
     "iopub.status.busy": "2021-04-16T19:35:20.899837Z",
     "iopub.status.idle": "2021-04-16T19:35:20.902096Z",
     "shell.execute_reply": "2021-04-16T19:35:20.902445Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "333.41989326370776"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(posterior.ps * posterior.qs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can use the method provided by `Pmf`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:20.906600Z",
     "iopub.status.busy": "2021-04-16T19:35:20.905906Z",
     "iopub.status.idle": "2021-04-16T19:35:20.909638Z",
     "shell.execute_reply": "2021-04-16T19:35:20.909046Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "333.41989326370776"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "posterior.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean of the posterior is 333, so that might be a good guess if you want to minimize error.\n",
    "If you played this guessing game over and over, using the mean of the posterior as your estimate would minimize the [mean squared error](http://en.wikipedia.org/wiki/Minimum_mean_square_error) over the long run."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sensitivity to the Prior\n",
    "\n",
    "The prior I used in the previous section is uniform from 1 to 1000, but I offered no justification for choosing a uniform distribution or that particular upper bound.\n",
    "We might wonder whether the posterior distribution is sensitive to the prior.\n",
    "With so little data---only one observation---it is.\n",
    "\n",
    "This table shows what happens as we vary the upper bound:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:20.921293Z",
     "iopub.status.busy": "2021-04-16T19:35:20.920199Z",
     "iopub.status.idle": "2021-04-16T19:35:20.935921Z",
     "shell.execute_reply": "2021-04-16T19:35:20.935395Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "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>Posterior mean</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Upper bound</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>500</th>\n",
       "      <td>207.079228</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1000</th>\n",
       "      <td>333.419893</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000</th>\n",
       "      <td>552.179017</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             Posterior mean\n",
       "Upper bound                \n",
       "500              207.079228\n",
       "1000             333.419893\n",
       "2000             552.179017"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.DataFrame(columns=['Posterior mean'])\n",
    "df.index.name = 'Upper bound'\n",
    "\n",
    "for high in [500, 1000, 2000]:\n",
    "    hypos = np.arange(1, high+1)\n",
    "    pmf = Pmf(1, hypos)\n",
    "    update_train(pmf, data=60)\n",
    "    df.loc[high] = pmf.mean()\n",
    "    \n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we vary the upper bound, the posterior mean changes substantially.\n",
    "So that's bad.  \n",
    "\n",
    "When the posterior is sensitive to the prior, there are two ways to proceed:\n",
    "\n",
    "* Get more data.\n",
    "\n",
    "* Get more background information and choose a better prior.\n",
    "\n",
    "With more data, posterior distributions based on different priors tend to converge.  \n",
    "For example, suppose that in addition to train 60 we also see trains 30 and 90.\n",
    "\n",
    "Here's how the posterior means depend on the upper bound of the prior, when we observe three trains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:20.946439Z",
     "iopub.status.busy": "2021-04-16T19:35:20.945991Z",
     "iopub.status.idle": "2021-04-16T19:35:20.958263Z",
     "shell.execute_reply": "2021-04-16T19:35:20.958659Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "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>Posterior mean</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Upper bound</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>500</th>\n",
       "      <td>151.849588</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1000</th>\n",
       "      <td>164.305586</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000</th>\n",
       "      <td>171.338181</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             Posterior mean\n",
       "Upper bound                \n",
       "500              151.849588\n",
       "1000             164.305586\n",
       "2000             171.338181"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame(columns=['Posterior mean'])\n",
    "df.index.name = 'Upper bound'\n",
    "\n",
    "dataset = [30, 60, 90]\n",
    "\n",
    "for high in [500, 1000, 2000]:\n",
    "    hypos = np.arange(1, high+1)\n",
    "    pmf = Pmf(1, hypos)\n",
    "    for data in dataset:\n",
    "        update_train(pmf, data)\n",
    "    df.loc[high] = pmf.mean()\n",
    "    \n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The differences are smaller, but apparently three trains are not enough for the posteriors to converge."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power Law Prior\n",
    "\n",
    "If more data are not available, another option is to improve the\n",
    "priors by gathering more background information.\n",
    "It is probably not reasonable to assume that a train-operating company with 1000 locomotives is just as likely as a company with only 1.\n",
    "\n",
    "With some effort, we could probably find a list of companies that\n",
    "operate locomotives in the area of observation.\n",
    "Or we could interview an expert in rail shipping to gather information about the typical size of companies.\n",
    "\n",
    "But even without getting into the specifics of railroad economics, we\n",
    "can make some educated guesses.\n",
    "In most fields, there are many small companies, fewer medium-sized companies, and only one or two very large companies.\n",
    "\n",
    "In fact, the distribution of company sizes tends to follow a power law, as Robert Axtell reports in *Science* (<http://www.sciencemag.org/content/293/5536/1818.full.pdf>).\n",
    "\n",
    "This law suggests that if there are 1000 companies with fewer than\n",
    "10 locomotives, there might be 100 companies with 100 locomotives,\n",
    "10 companies with 1000, and possibly one company with 10,000 locomotives.\n",
    "\n",
    "Mathematically, a power law means that the number of companies with a given size, $N$, is proportional to $(1/N)^{\\alpha}$, where $\\alpha$ is a parameter that is often near 1.\n",
    "\n",
    "We can construct a power law prior like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:20.963600Z",
     "iopub.status.busy": "2021-04-16T19:35:20.963012Z",
     "iopub.status.idle": "2021-04-16T19:35:20.965471Z",
     "shell.execute_reply": "2021-04-16T19:35:20.965831Z"
    },
    "tags": [
     "hide-output"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.178368103610282"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha = 1.0\n",
    "ps = hypos**(-alpha)\n",
    "power = Pmf(ps, hypos, name='power law')\n",
    "power.normalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For comparison, here's the uniform prior again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:20.970595Z",
     "iopub.status.busy": "2021-04-16T19:35:20.969882Z",
     "iopub.status.idle": "2021-04-16T19:35:20.972668Z",
     "shell.execute_reply": "2021-04-16T19:35:20.972223Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1000"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hypos = np.arange(1, 1001)\n",
    "uniform = Pmf(1, hypos, name='uniform')\n",
    "uniform.normalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what a power law prior looks like, compared to the uniform prior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:20.993580Z",
     "iopub.status.busy": "2021-04-16T19:35:20.987921Z",
     "iopub.status.idle": "2021-04-16T19:35:21.145573Z",
     "shell.execute_reply": "2021-04-16T19:35:21.145133Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAf9UlEQVR4nO3df5QdZZ3n8fenu9OBJECQtBCTQIJEmMiqYIw4jj9m1BnCKlHXmQ2jIo57sqyA4OjswjhndOYcZ3QURxk55ERgBEGRUVzjGAV1wVldgTQRAyFE2gimTYBGlAABku7+7h9Vt1N9c7vv7U5X15PO53XOPalbVU/V91Z37qefus+tUkRgZmaWmraqCzAzM2vEAWVmZklyQJmZWZIcUGZmliQHlJmZJckBZWZmSXJAmQGSnpJ0fEnbfr2k3sLzTZJeP0HbfqekWwrPQ9IJE7HtfHulHRezZhxQNiVJelDSM/kb7COS/lXSrJHWj4hZEbF1MmqLiBdHxG2jrSNpYR42HU22dX1E/PFE1CXpNkn/rW77k3ZczOo5oGwqe0tEzAJOBV4B/E39Cs0CoJn9bX+g7ttsMjigbMqLiF8D3wFOhqHTYOdJegB4oDDvhHz6CEnXSuqT9JCkv5HUli87R9KPJf2zpMeBj9XvT9Khkr4o6beS7iMLx+LyByW9MZ9eJqlb0s68p/eZfLX/yP/9Xd4LfFWjfefzflRXwhmStkp6TNKnCrV/TNJ1hTqGemmSPg68Bvh8vr/Pj+O4/EjSp/PX/UtJywv7Oiev6cl82TvH8jO0g5P/ArMpT9IC4AzgpsLstwKvBJ5p0ORfgCOA44GjgFuAHcBV+fJXAjcAzwemNWj/UeCF+WMmWTiO5HPA5yLiS/kpyJPz+a8FfgnMjoj+/HWc2GDf/7XBNt8GLAVmAd8HtgBXjlIDEfERSa8GrouIkdZt5bhcA8wBVgFXSZoHzAAuA14REVskzQWeN1o9ZuAelE1t/1vS74AfAT8E/qGw7B8j4vGIGBZQktrJ3vQviYgnI+JB4FLg3YXVtkfEv0REf3373J8BH8+3v43szXkke4ATJM2JiKci4vYmr6nZvgE+me/7V8BngbOabLOpFo/LQxHxhYgYIAuqucDR+bJB4GRJh0bEjojYtL812dTngLKp7K0RMTsijouI99e9oW8boc0coBN4qDDvIWBeC21rXlC3zkMjrQi8D3gRcL+k9ZLe3GTbzfZdv85DeT37q5Xj8nBtIiJ25ZOzIuJpsnA7F9gh6duSTpqAmmyKc0DZwWqky/g/RtarOa4w71jg1y20rdkBLKhr37iIiAci4iyyU3afBL4maeYo+2jl9gP1+96eTz9Ndrqt5pgxbLuV4zKiiLg5It5E1qu6H/hCK+3s4OaAMivIT0/dCHxc0mGSjgP+Erhu9JbD3AhcIulISfOBC0ZaUdK7JHVFxCDwu3z2ANBHdlpsPN9B+qt83wuAC4Gv5vPvBl4r6VhJRwCX1LV7ZKT97c9xkXS0pDPz4H0OeIrsNZqNygFltq8LyHobW8k+v/oycPUY2v8d2emvX5INJPjSKOueDmyS9BTZgImVEfFsfors48CPJf1O0mlj2P83gbvIAunb5IMYIuJ7ZGG1MV/+73XtPge8Ix+F1+hzs/EelzbgQ2Q9uceB1wHvH8PrsYOUfMNCMzNLkXtQZmaWJAeUmZklyQFlZmZJckCZmVmSptSljubMmRMLFy6sugwzMxuDu+6667GI6KqfP6UCauHChXR3d1ddhpmZjYGkhldb8Sk+MzNLkgPKzMyS5IAyM7MkOaDMzCxJDigzM0uSA8rMzJI0pYaZj1f074ZndmZP2jvQjNmV1mNmZg4oAOKxB9lzy6UAtB29mGl/8uGKKzIzM5/iMzOzJJUaUJJOl7RFUo+kixssP0nSTyQ9J+nDhfkLJN0qabOkTZIuLLPOYXx/LDOzJJR2ik9SO3A58CagF1gvaW1E3FdY7XHgA8Bb65r3Ax+KiA2SDgPukvS9urZmZjaFldmDWgb0RMTWiNgN3ACsKK4QEY9GxHpgT938HRGxIZ9+EtgMzCuxVjMzS0yZATUP2FZ43ss4QkbSQuAU4I4Rlq+S1C2pu6+vbzx1gjS+dmZmVpoyA6rRu/6YPuCRNAv4OnBRROxstE5ErImIpRGxtKtrn6u1j50/gzIzS0KZAdULLCg8nw9sb7WxpGlk4XR9RNw0wbXV763czZuZ2ZiVGVDrgcWSFknqBFYCa1tpKEnAVcDmiPhMiTWamVmiShvFFxH9ks4HbgbagasjYpOkc/PlqyUdA3QDhwODki4ClgAvAd4N3CPp7nyTfx0R68qq18zM0lLqlSTyQFlXN291YfphslN/9X7EZJ538xk+M7Pk+EoS+/AgCTOzFDigzMwsSQ6oeh5mbmaWBAcU4A+hzMzS44AyM7MkOaDAlzoyM0uQA8rMzJLkgNqHB0mYmaXAAWVmZklyQNUJDzM3M0uCAwo8SMLMLEEOKDMzS5IDCvAXdc3M0uOAMjOzJDmgzMwsSQ4oMzNLkgOqnoeZm5klwQEFHmZuZpYgB1Q996DMzJLggAI8zNzMLD0OKDMzS5IDyszMkuSAAp/hMzNLkANqHx4kYWaWglIDStLpkrZI6pF0cYPlJ0n6iaTnJH14LG3NzGxqKy2gJLUDlwPLgSXAWZKW1K32OPAB4NPjaFsODzM3M0tCmT2oZUBPRGyNiN3ADcCK4goR8WhErAf2jLXtxPKHUGZmqSkzoOYB2wrPe/N5E9pW0ipJ3ZK6+/r6xlWomZmlp8yAatQtafX8WcttI2JNRCyNiKVdXV0tF2dmZmkrM6B6gQWF5/OB7ZPQdux8LT4zs+SUGVDrgcWSFknqBFYCayeh7X7yIAkzsxR0lLXhiOiXdD5wM9AOXB0RmySdmy9fLekYoBs4HBiUdBGwJCJ2NmpbVq1mZpae0gIKICLWAevq5q0uTD9MdvqupbaTwsPMzcyS4CtJAB5mbmaWHgeUmZklyQEFHsVnZpYgB5SZmSXJAWVmZklyQAEeJGFmlh4HVD0PMzczS4IDyszMkuSA2od7UGZmKXBAAfIwczOz5DigzMwsSQ4oMzNLkgPKzMyS5ICq52HmZmZJcECZmVmSHFD7cA/KzCwFDijw1czNzBLkgDIzsyQ5oMzMLEkOKMBXMzczS48Dqp6HmZuZJcEBZWZmSXJA7cM9KDOzFJQaUJJOl7RFUo+kixssl6TL8uUbJZ1aWPZBSZsk3SvpK5IOKbHQ0jZtZmbjU1pASWoHLgeWA0uAsyQtqVttObA4f6wCrsjbzgM+ACyNiJOBdmBlWbUO4w6UmVkSyuxBLQN6ImJrROwGbgBW1K2zArg2MrcDsyXNzZd1AIdK6gBmANvLK9U9KDOz1JQZUPOAbYXnvfm8putExK+BTwO/AnYAT0TELY12ImmVpG5J3X19fRNWvJmZVavMgGrULak/gdZwHUlHkvWuFgEvAGZKelejnUTEmohYGhFLu7q69qtgMzNLR5kB1QssKDyfz76n6UZa543ALyOiLyL2ADcBv19apR4kYWaWnDIDaj2wWNIiSZ1kgxzW1q2zFjg7H813GtmpvB1kp/ZOkzRDkoA3AJtLrLXAoyTMzFLQUdaGI6Jf0vnAzWSj8K6OiE2Szs2XrwbWAWcAPcAu4L35sjskfQ3YAPQDPwXWlFWrmZmlp7SAAoiIdWQhVJy3ujAdwHkjtP0o8NEy62vIlzoyM0uCryQBeJi5mVl6HFBmZpYkB5SZmSXJAQU+w2dmliAH1D48SMLMLAUOKDMzS5IDqp6HmZuZJcEBBb7UkZlZghxQZmaWJAeUmZklyQEFeJy5mVl6HFD78CAJM7MUOKDMzCxJDqh6HmZuZpYEBxR4mLmZWYIcUGZmliQHFOBRfGZm6XFAmZlZkkYNKElfLEy/p/RqzMzMcs16UC8tTF9YZiGV8iAJM7PkNAuog2/MtYeZm5kloaPJ8vmSLiMbRVCbHhIRHyitMjMzO6g1C6i/Kkx3l1lIMtyDMjNLwqgBFRHXTFYh1fJnUGZmqRk1oCStHW15RJzZpP3pwOeAduDKiPhE3XLly88AdgHnRMSGfNls4ErgZLLPwv4iIn4y2v7MzGzqaHaK71XANuArwB2MoashqR24HHgT0Ausl7Q2Iu4rrLYcWJw/Xglckf8LWXB9NyLeIakTmNHqvsds2KvyKT4zsxQ0C6hjyALmLODPgW8DX4mITS1sexnQExFbASTdAKwAigG1Arg2IgK4XdJsSXOBp4HXAucARMRuYHerL2rsCgnlz6DMzJIw6jDziBiIiO9GxHuA04Ae4DZJF7Sw7Xlkva+a3nxeK+scD/QB/yrpp5KulDSz0U4krZLULam7r6+vhbIabqTwxAFlZpaCppc6kjRd0tuB64DzgMuAm1rYdqPTgfXv/iOt0wGcClwREaeQ9agubrSTiFgTEUsjYmlXV1cLZTUp1flkZpaEZoMkriEbpPAd4O8i4t4xbLsXWFB4Ph/Y3uI6AfRGxB35/K8xQkBNCO3N6XBCmZkloVkP6t3Ai8guc/QTSTvzx5OSdjZpux5YLGlRPshhJVA/KnAtcLYypwFPRMSOiHgY2CbpxHy9NzD8s6vyxOCk7MbMzEbX7HtQ477aeUT0SzofuJlsmPnVEbFJ0rn58tXAOrIh5j1kw8zfW9jEBcD1ebhtrVs2sXwtPjOz5DQ7xXcIcC5wArCRLGT6W914RKwjC6HivNWF6SD7XKtR27uBpa3ua78UA8o9KDOzJDTrIV1DFhL3kPV0Li29okp4mLmZWWqafQ9qSUT8JwBJVwF3ll9SBeSAMjNLTbMe1J7axFhO7R1w5BsLm5mlplkP6qWF0XoCDs2fi+wjpMNLra4SQUQgD5wwM6tUs1F87ZNVSJUcRmZm6fG5rSH+HMrMLCUOqBoPNTczS4oDqsYXjDUzS4oDaohP8ZmZpcQBVeOBEmZmSXFA1RS/C+UelJlZ5RxQjXiQhJlZ5RxQOXmQhJlZUhxQQzxIwswsJQ6oGvegzMyS4oAaUuxBVVeFmZllHFA1bcVD4YQyM6uaA6oRj+IzM6ucA6rGNy00M0uKA2qIB0mYmaXEAVXjK0mYmSXFAdWIA8rMrHIOqBr5UJiZpaTUd2VJp0vaIqlH0sUNlkvSZfnyjZJOrVveLumnkv69zDqznRWmPYrPzKxypQWUpHbgcmA5sAQ4S9KSutWWA4vzxyrgirrlFwKby6pxOI/iMzNLSZk9qGVAT0RsjYjdwA3Airp1VgDXRuZ2YLakuQCS5gP/GbiyxBr38qWOzMySUmZAzQO2FZ735vNaXeezwP8EJud8m0fxmZklpcyAanSL2vp3/obrSHoz8GhE3NV0J9IqSd2Suvv6+sZT5z6lhAPKzKxyZQZUL7Cg8Hw+sL3FdV4NnCnpQbJTg38k6bpGO4mINRGxNCKWdnV1jb9an+IzM0tKmQG1HlgsaZGkTmAlsLZunbXA2flovtOAJyJiR0RcEhHzI2Jh3u7/RMS7Sqx1OPegzMwq11HWhiOiX9L5wM1AO3B1RGySdG6+fDWwDjgD6AF2Ae8tq56m5KuZm5mlpLSAAoiIdWQhVJy3ujAdwHlNtnEbcFsJ5Y2yU38Pysysar58Qo2vJGFmlhS/K9f4dhtmZklxQNU4oMzMkuKAasSfQZmZVc4BlZM/gzIzS4rflWuGneJzD8rMrGoOqCHFgKquCjMzyziganypIzOzpDigajyKz8wsKQ6oRvwZlJlZ5RxQNcPuB+WAMjOrmgOqpq1wWcLBgerqMDMzwAG1V1v73mn3oMzMKueAqikEVAz2V1iImZmBA2qvYg9q0D0oM7OqOaBqhgWUe1BmZlVzQOUk96DMzFLigKoZ1oPaU10dZmYGOKD2GjbM3D0oM7OqOaBq2gqHwp9BmZlVzgFV4x6UmVlSHFA1hR6UvwdlZlY9B1SNL3VkZpYUB1TNsIvFOqDMzKpWakBJOl3SFkk9ki5usFySLsuXb5R0aj5/gaRbJW2WtEnShWXWCaB296DMzFJSWkAp++br5cByYAlwlqQldastBxbnj1XAFfn8fuBDEfF7wGnAeQ3aTqxh34NyQJmZVa3MHtQyoCcitkbEbuAGYEXdOiuAayNzOzBb0tyI2BERGwAi4klgMzCvxFqHB9SAB0mYmVWtzICaB2wrPO9l35Bpuo6khcApwB2NdiJplaRuSd19fX3jr7Z92t7pAV9JwsysamUGlBrMi7GsI2kW8HXgoojY2WgnEbEmIpZGxNKurq5xF0vH9L3bHNg9/u2YmdmEKDOgeoEFhefzge2triNpGlk4XR8RN5VYJwDq6Nz7pN8BZWZWtTIDaj2wWNIiSZ3ASmBt3TprgbPz0XynAU9ExA5JAq4CNkfEZ0qsca92B5SZWUo6mq8yPhHRL+l84GagHbg6IjZJOjdfvhpYB5wB9AC7gPfmzV8NvBu4R9Ld+by/joh1ZdVLoQcV/c+VthszM2tNaQEFkAfKurp5qwvTAZzXoN2PaPz5VHl8is/MLCm+kkROxVN8HiRhZlY5B1TNtL2j+NyDMjOrngOqpuOQocnY/QzZ2UczM6uKA6qmo3PvFc0H+92LMjOrmAMqJwlNn7l3xnNPV1eMmZk5oIaZPmtoMnY/VWEhZmbmgCrQ9Bl7nzzrgDIzq5IDqqjYg/IpPjOzSjmgCoZ/BuUelJlZlRxQRTOOHJqMXb+tsBAzM3NAFWjWUUPT8eRjFVZiZmYOqALNmjM0HU85oMzMquSAKij2oHiyz1eTMDOrkAOqaMaRqDMbah67d4F7UWZmlXFAFUhCRy0cej742IOV1WJmdrBzQNXRnIVD0/FoT3WFmJkd5BxQddqOXjw0Pbjtbn8OZWZWEQdUHR1zIurMvrAbu35HPHx/xRWZmR2cHFB11NZO28KlQ88HNn7bvSgzswo4oBpof/GbQNmhGXzkAQY3/6DiiszMDj4OqAZ0WBftS9449Ly/+2sMbP6Be1JmZpPIATWC9pe+hbahEX1B//ob6b/lUgZ33O+gMjObBB1VF5AqdXTS8YYL2PP9y4jfPATkp/u+989o1lG0zX8JOubELMQOnY2kags2M5tiNJV6A0uXLo3u7u4J3WYM7GFgwzcYuP9WiMGG66hzBpr9Aph1FJpxJJp5JJpxJEyfCZ0zs9t4dM5A7f57wMysnqS7ImJp/fxS3zElnQ58DmgHroyIT9QtV778DGAXcE5EbGil7US75C3Xjbp8VsdLOOnw+zl25q/obNvdYI2fN91Hf3SwZ3Aa/YMd9EcHA9HOQLTTHx30D+59PhDtDEYbg7QRoWHTA9FOkM/L5w9GtiwoPKLBdN08YGhbtelG7WrLaoamg4bLwb3Jsvzjt95VdQlmk6a0gJLUDlwOvAnoBdZLWhsR9xVWWw4szh+vBK4AXtli20n1VP9hdD/+Cjb89lSeP/1Rnn9I9jhi2hNMa9vT0jY61E9He38WuQeBofCKBuHWYHnDEBxteey7bqPno86PJstb2G7D5Q1OTOzv9gF2f/tXAIVTyiO0aXl5k3nFbYzpNHaTdffrlHizbe9n+1GbNmu7n3+clbn9Ztvej+PWvvgPaJt38phLaqbMHtQyoCcitgJIugFYARRDZgVwbWTnGW+XNFvSXGBhC20rMRjtPPzsXB5+dm4+Jzi0/RkOn7aTmR1Pc2j7M8zo2MWM9l10tu1mWtseprc/R2fbbtToXWsKq71eqcnrdoerZbXPQw+u3yRLXdsxJ5ay3TIDah6wrfC8l6yX1GydeS22BUDSKmAVwLHHHjvuYss+dRIRsOdZ2PMM9O8mBnZD/95H9O+G2rzBfhgcIAYH8unBoXnDlw3AQD/EQPb5WET+b3E6iBjMtkE0WG8we7cb9ryubfYC2NtlqZse+hzTb5tmNnHKDKhGfxfXv4ONtE4rbbOZEWuANZANkhhLgZNJEnQemj2Y2p2GoYE3tSAb9i+FZTRYPnTejn0CsX47+6xTeD5s3rDqmixv8Dpaat/onF6z5WOpq8ny8exrxNc/hloaNi/xv2HTbe/nvkfdfpnbngCjbD9K/pm1HTl/v9qPpMyA6gUWFJ7PB7a3uE5nC20tUUOfjxzgQ+8P7OrNDnxlflF3PbBY0iJJncBKYG3dOmuBs5U5DXgiIna02NbMzKaw0npQEdEv6XzgZrJxa1dHxCZJ5+bLVwPryIaY95ANM3/vaG3LqtXMzNLjL+qamVmlRvqirq/FZ2ZmSXJAmZlZkhxQZmaWJAeUmZklaUoNkpDUBzy0H5uYAzw2QeWUyXVOnAOhRnCdE+lAqBEOrjqPi4iu+plTKqD2l6TuRiNJUuM6J86BUCO4zol0INQIrhN8is/MzBLlgDIzsyQ5oIZbU3UBLXKdE+dAqBFc50Q6EGoE1+nPoMzMLE3uQZmZWZIcUGZmliQHVE7S6ZK2SOqRdHGFdSyQdKukzZI2Sbown/8xSb+WdHf+OKPQ5pK87i2S/mQSa31Q0j15Pd35vOdJ+p6kB/J/j6yyTkknFo7Z3ZJ2Srqo6uMp6WpJj0q6tzBvzMdO0svzn0GPpMukib0J1wh1fkrS/ZI2SvqGpNn5/IWSnikc09UV1znmn3FFdX61UOODku7O51dyPEd5D5r838+IOOgfZLf0+AVwPNnNEn8GLKmolrnAqfn0YcDPgSXAx4APN1h/SV7vdGBR/jraJ6nWB4E5dfP+Cbg4n74Y+GTVddb9nB8Gjqv6eAKvBU4F7t2fYwfcCbyK7P6K3wGWT0Kdfwx05NOfLNS5sLhe3XaqqHPMP+Mq6qxbfinwt1UeT0Z+D5r030/3oDLLgJ6I2BoRu4EbgBVVFBIROyJiQz79JLAZmDdKkxXADRHxXET8kuzeWsvKr3TUeq7Jp68B3lqYX3WdbwB+ERGjXW1kUuqMiP8AHm+w75aPnaS5wOER8ZPI3g2uLbQprc6IuCUi+vOnt5Pd8XpEVdU5iqSOZ03eu/gz4CujbaPsOkd5D5r0308HVGYesK3wvJfRQ2FSSFoInALckc86Pz+tcnWhe11l7QHcIukuSavyeUdHdldk8n+fn0CdNSsZ/p8/teM51mM3L5+unz+Z/oLsL+OaRZJ+KumHkl6Tz6uyzrH8jKs+nq8BHomIBwrzKj2ede9Bk/776YDKNDovWun4e0mzgK8DF0XETuAK4IXAy4AdZKcCoNraXx0RpwLLgfMkvXaUdSs9xpI6gTOBf8tnpXg8RzJSTVUf048A/cD1+awdwLERcQrwl8CXJR1OdXWO9Wdc9c/+LIb/AVXp8WzwHjTiqiPUs991OqAyvcCCwvP5wPaKakHSNLJfjOsj4iaAiHgkIgYiYhD4AntPO1VWe0Rsz/99FPhGXtMjede+diri0arrzC0HNkTEI5Dm8WTsx66X4afXJq1WSe8B3gy8Mz99Q36K5zf59F1kn0W8qKo6x/EzrvJ4dgBvB75am1fl8Wz0HkQFv58OqMx6YLGkRflf2iuBtVUUkp+HvgrYHBGfKcyfW1jtbUBtFNBaYKWk6ZIWAYvJPpgsu86Zkg6rTZN9cH5vXs978tXeA3yzyjoLhv11mtrxLOy75WOXn2Z5UtJp+e/N2YU2pZF0OvC/gDMjYldhfpek9nz6+LzOrRXWOaafcVV15t4I3B8RQ6fEqjqeI70HUcXv50SN/DjQH8AZZKNVfgF8pMI6/oCsG7wRuDt/nAF8Cbgnn78WmFto85G87i1M8KijUeo8nmzkzs+ATbVjBhwF/AB4IP/3eVXWme93BvAb4IjCvEqPJ1lY7gD2kP2l+b7xHDtgKdkb7y+Az5NfHabkOnvIPnOo/X6uztf9L/nvws+ADcBbKq5zzD/jKurM538ROLdu3UqOJyO/B03676cvdWRmZknyKT4zM0uSA8rMzJLkgDIzsyQ5oMzMLEkOKDMzS5IDyg56kkLSpYXnH5b0sQna9hclvWMittVkP3+aX3361rr5CyX9+Ti3+f8mpjqz8XFAmcFzwNslzam6kKLalzRb9D7g/RHxh3XzFwINAyq/esGIIuL3x7B/swnngDLLrie3Bvhg/YL6HpCkp/J/X59fwPNGST+X9AlJ75R0Z37/mxcWNvNGSf83X+/Neft2ZfdVWp9fzPS/F7Z7q6Qvk33JtL6es/Lt3yvpk/m8vyX7cuVqSZ+qa/IJ4DXK7if0QUnnSPo3Sd8iu9DvLEk/kLQh3+6Kwr6Kr/U2SV9Tdh+o6/MrA5C/7vvy1/DpMR95s1GM+heU2UHkcmCjpH8aQ5uXAr9HdvuErcCVEbFM2Q3eLgAuytdbCLyO7MKlt0o6geyyL09ExCskTQd+LOmWfP1lwMmR3bpgiKQXkN1/6eXAb8kC5q0R8feS/ojs3kfddTVenM+vBeM5ZPfneUlEPJ73ot4WETvzHuTtktbGvt/gPwV4Mdm11H4MvFrSfWSXEDopIkL5jQvNJop7UGZAZFdrvhb4wBiarY/s3jnPkV3KpRYw95CFUs2NETEY2W0UtgInkV278Gxld0+9g+wyMovz9e+sD6fcK4DbIqIvsvsxXU92A7yx+l5E1O5JJOAfJG0Evk92O4SjG7S5MyJ6I7vw6t3569sJPAtcKentwK4G7czGzQFlttdnyT7LmVmY10/+/yQ/rdVZWPZcYXqw8HyQ4Wcn6nsjtVsRXBARL8sfiyKiFnBPj1DfRN3Wu7j9dwJdwMsj4mXAI8AhDdoUX+sA2R11+8l6e18nuxHddyeoPjPAAWU2JO9V3EgWUjUPkp1Sg+zOodPGsek/ldSWfy51PNkFNW8G/oey2xog6UX5VeFHcwfwOklz8gEUZwE/bNLmSbLbdo/kCODRiNgj6Q+B41p4PeQ1zyK7AO86stOZL2u1rVkr/BmU2XCXAucXnn8B+KakO8mu4DxS72Y0W8iC5GiyK1Y/K+lKstNkG/KeWR9NbocdETskXQLcStabWhcRzW5fsBHol/Qzsitm/7Zu+fXAtyR1k526u7/1l8VhZMfmkLyefQaZmO0PX83czMyS5FN8ZmaWJAeUmZklyQFlZmZJckCZmVmSHFBmZpYkB5SZmSXJAWVmZkn6/60E3md6XhntAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "uniform.plot(color='C4')\n",
    "power.plot(color='C1')\n",
    "\n",
    "decorate(xlabel='Number of trains',\n",
    "         ylabel='PMF',\n",
    "         title='Prior distributions')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the update for both priors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.149556Z",
     "iopub.status.busy": "2021-04-16T19:35:21.148612Z",
     "iopub.status.idle": "2021-04-16T19:35:21.152684Z",
     "shell.execute_reply": "2021-04-16T19:35:21.152290Z"
    }
   },
   "outputs": [],
   "source": [
    "dataset = [60]\n",
    "update_train(uniform, dataset)\n",
    "update_train(power, dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here are the posterior distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.168591Z",
     "iopub.status.busy": "2021-04-16T19:35:21.168151Z",
     "iopub.status.idle": "2021-04-16T19:35:21.299407Z",
     "shell.execute_reply": "2021-04-16T19:35:21.299044Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAy6ElEQVR4nO3de5xddX3v/9d77rmHkJCEJJAAQQkgFyNSFYpileAlaqtCPQWp50E5hdrW9nfE9vRobT0/tbX2R7WkqFRoVcSqNVYsUH7iHU24CAQIJCHAkMmFBHKbzGXP/pw/1pqZNTt7ZvZMZu29k7yfD/Zj1vqu73etz1oz7E++a3/3dykiMDMzqzcNtQ7AzMysHCcoMzOrS05QZmZWl5ygzMysLjlBmZlZXXKCMjOzuuQEZTaBJL1P0l057v/Lkv46Xb5A0voJ3Pf3JV2ZLr9f0k8mcN+5Xhc7MjlB2WFL0mZJByTtk7RN0j9LmnoI+/uYpH89lJgi4isR8aZD2ccYjvXjiHjZaPUqPa+IWBERtxxqXJIWSwpJTZl9V+262JHDCcoOd2+LiKnAucCrgP9Vq0Cyb8jjaCtJNfn/sZbHNhuJ/yjtiBARzwPfB84AkPR2SeskvSTpXkmn9deV9GFJz0vaK2m9pIslXQL8GfDetEf2q7TuDElfktSRtvlrSY3ptvdL+qmkz0raBXys9NaYpNdIWiNpd/rzNZlt90r6hKSfAp3ASaXnJekcSQ+ksX4daMtsu0hS+zjP66Bjp2X/fejh9Q9p7E9IujizYbOkN2bWs720H6U/X0qP+WvjuC5/lV7bvZLukjQ73dYm6V8l7Ux/t2skzS37R2GHPScoOyJIWgRcCjwo6VTga8AfAXOAO4DvSmqR9DLgOuBVETENeDOwOSL+E/g/wNcjYmpEnJXu+hagAJwCnAO8Cci+ib8a2AQcB3yiJKZZwPeAG4Bjgb8Dvifp2Ey13wGuBqYBz5S0bwH+HfgXYBbwDeA3hzn/sZ7XiMcuObfZwEeBb6XnNJoL058z02P+vCTWSq7LbwNXkVzXFuBP0/IrgRnAorTtNcCBCmKyw5ATlB3u/l3SS8BPgB+SvBm/F/heRNwdEb3A3wKTgNcAfUArsExSc0RsjoiN5Xac/st8BfBHEbE/IrYDnwUuy1TbEhH/EBGFiCh9o3wL8FRE/Eu6/WvAE8DbMnW+HBHr0u29Je3PB5qBv4+I3oj4N2DNMNeh4vOq8NgA2zPH/jqwPj2nQ1XJdfnniHgyvaa3A2en5b0kiemUiOiLiPsjYs8ExGR1yAnKDnfviIiZEXFiRPx++oZ2PJkeQUQUgeeABRGxgaRn9TFgu6TbJB0/zL5PJEkQHentpJeAfyL5V32/50aIbUgcqWeABWNo/3wMndG5XE+HMZ5XJcdmmGOPts9KVHJdtmaWO4H+wS//AtwJ3CZpi6RPS2qegJisDjlB2ZFoC0lyAZIPUkhuCT0PEBFfjYjXpXUC+FRatXRq/+eAbmB2mgRnRsT0iDg9U2ekxwEMiSN1Qn8cFbTvABak8WfblzWG86rk2Axz7C3p8n5gcmbbvDHst5LrUlbam/vLiFhG0iN+K3DFaO3s8OQEZUei24G3pIMEmoE/IUk0P5P0MklvkNQKdJF8ftGXttsGLO4f0RYRHcBdwGckTZfUIOlkSb9eYRx3AKdK+m1JTZLeCywD/qPC9j8n+fzrg2n7dwHnlas4lvMag+PSYzdLejdwWnpOAA8Bl6XblgO/lWm3AyhSZtBHatzXRdLrJZ2ZDlTZQ3LLr2+UZnaYcoKyI05ErAf+G/APwAskn228LSJ6SD6n+WRavpXkTfjP0qbfSH/ulPRAunwFyYf0jwEvAv8GzK8wjp0k/8L/E2An8D+Bt0bECxW27wHeBbw/PfZ7gW8NU32s51WJXwBL031+Avit9JwA/gI4OY3rL4GvZuLuTOv/NL01en7JeR3KdZlH8jvYAzxO8rnjIX13zeqX/MBCMzOrR+5BmZlZXXKCMjOzuuQEZWZmdckJyszM6tK4J7c8EsyePTsWL15c6zDMzI5q999//wsRMae0/KhOUIsXL2bt2rW1DsPM7KgmqewMKbne4pN0STqr8gZJ15fZLkk3pNsflnRuZtvNkrZLerSkzdclPZS+Nkt6KC1frOTZQP3bVuV5bmZmlq/celDpN70/D/wG0A6skbQ6Ih7LVFtB8kXApSQzJ9+Y/gT4MvA54NbsfiPivZljfAbYndm8MSLOntATMTOzmsizB3UesCEiNqXfiL8NWFlSZyVwayTuA2ZKmg8QET8Cdg2383SOsPeQPFbBzMyOMHkmqAUMnS25naGzFVdaZzgXANsi4qlM2RJJD0r6oaQLyjWSdLWktZLW7tixo8JDmZlZteWZoFSmrHRepUrqDOdyhvaeOoATIuIc4EPAVyVNP2jnETdFxPKIWD5nzkGDRszMrE7kmaDaSR5x0G8hg1P1j6XOQSQ1kUyi+fX+sojo7p/IMiLuBzYCp44rcjMzq7k8E9QaYKmkJemjqy8DVpfUWQ1ckY7mOx/YnT7iYDRvBJ6IiPb+Aklz0oEZSDqJZODFpok4kUMRhR6iWKx1GGZmh53cRvFFREHSdSRPv2wEbo6IdZKuSbevInkuzKXABpKnZl7V317S14CLgNmS2oGPRsSX0s2XcfDgiAuBj0sqkDwf5pqIGHaQRTUUdz1L4a6/h6YWmt/yETRpRi3DMTM7rBzVj9tYvnx55PlF3Z5v/hmxP3l8TsOis2l+/f/I7VhmZocrSfdHxPLScs/Fl6P+5AQQLzxdw0jMzA4/TlDV0leodQRmZocVJ6gqiaITlJnZWDhBVUtfb60jMDM7rDhBVUt4qLmZ2Vg4QeWpobHWEZiZHbacoHKkhuZah2BmdthygspT41H9PEgzs0PiBJUnJygzs3FzgspTgxOUmdl4OUHlyQnKzGzcnKDy1OhBEmZm4+UElSP5Mygzs3FzgsqTE5SZ2bg5QeXJn0GZmY2bE1SenKDMzMbNCSpPvsVnZjZuTlB58ig+M7Nxc4LKkXyLz8xs3Jyg8lQym3lE1CgQM7PDT64JStIlktZL2iDp+jLbJemGdPvDks7NbLtZ0nZJj5a0+Zik5yU9lL4uzWz7SLqv9ZLenOe5VUQll7fYV5s4zMwOQ7klKEmNwOeBFcAy4HJJy0qqrQCWpq+rgRsz274MXDLM7j8bEWenrzvS4y0DLgNOT9v9YxpD/fBj383MKpZnD+o8YENEbIqIHuA2YGVJnZXArZG4D5gpaT5ARPwI2DWG460EbouI7oh4GtiQxlA/+pygzMwqlWeCWgA8l1lvT8vGWqec69JbgjdLOmYs+5J0taS1ktbu2LGjgkMdipLPnIq9OR/PzOzIkWeCUpmy0lECldQpdSNwMnA20AF8Ziz7ioibImJ5RCyfM2fOKIeaYAUnKDOzSuWZoNqBRZn1hcCWcdQZIiK2RURfRBSBLzB4G2/M+6q28CAJM7OK5Zmg1gBLJS2R1EIygGF1SZ3VwBXpaL7zgd0R0THSTvs/o0q9E+gf5bcauExSq6QlJAMvfjkRJzJhfIvPzKxiuX2TNCIKkq4D7gQagZsjYp2ka9Ltq4A7gEtJBjR0Alf1t5f0NeAiYLakduCjEfEl4NOSzia5fbcZ+L10f+sk3Q48BhSAayOivrosHiRhZlaxXKc6SIeA31FStiqzHMC1w7S9fJjy3xnheJ8APjGuYPNQ+sVcJygzs4p5Jolq8vegzMwq5gRVTX3+DMrMrFJOUFUUTlBmZhVzgqomfwZlZlYxJ6hq8vegzMwq5gRVTX09tY7AzOyw4QRVTf4MysysYk5QVRQF96DMzCrlBJWr0i/qOkGZmVXKCaqa3IMyM6uYE1Q1+TMoM7OKOUFVkT+DMjOrnBNUNbkHZWZWMSeoPJXOZu4elJlZxZygqsmj+MzMKuYEVU2+xWdmVjEnqCryIAkzs8o5QVWTe1BmZhVzgqom96DMzCrmBFVNTlBmZhXLNUFJukTSekkbJF1fZrsk3ZBuf1jSuZltN0vaLunRkjZ/I+mJtP63Jc1MyxdLOiDpofS1Ks9zGw8/UdfMrHK5JShJjcDngRXAMuBySctKqq0Alqavq4EbM9u+DFxSZtd3A2dExCuAJ4GPZLZtjIiz09c1E3IiE6lYIIrFWkdhZnZYyLMHdR6wISI2RUQPcBuwsqTOSuDWSNwHzJQ0HyAifgTsKt1pRNwVEf3PTr8PWJjbGeTBvSgzs4rkmaAWAM9l1tvTsrHWGcnvAt/PrC+R9KCkH0q6oFwDSVdLWitp7Y4dO8ZwqAniL+uamVUkzwSlMmUxjjrldy79OVAAvpIWdQAnRMQ5wIeAr0qaftDOI26KiOURsXzOnDmVHGpiFdyDMjOrRJ4Jqh1YlFlfCGwZR52DSLoSeCvwvohkwruI6I6Ineny/cBG4NRxR58X96DMzCqSZ4JaAyyVtERSC3AZsLqkzmrginQ03/nA7ojoGGmnki4BPgy8PSI6M+Vz0oEZSDqJZODFpok7nYkRTlBmZhVpymvHEVGQdB1wJ9AI3BwR6yRdk25fBdwBXApsADqBq/rbS/oacBEwW1I78NGI+BLwOaAVuFsSwH3piL0LgY9LKgB9wDURcdAgi5rzd6HMzCqSW4ICiIg7SJJQtmxVZjmAa4dpe/kw5acMU/5N4JvjDjYPpY/bACcoM7MKeSaJavMwczOzijhBVZlnNDczq4wTVLV5kISZWUWcoKqtt7vWEZiZHRacoHJVbpCEE5SZWSWcoKosnKDMzCriBFVtvV21jsDM7LDgBFVtHsVnZlYRJ6gqC/egzMwq4gRVbf4MysysIk5Q1eYEZWZWESeoKgt/D8rMrCJOUNXmHpSZWUWcoPJUbjZzD5IwM6uIE1S1uQdlZlYRJ6gqi95uolzPyszMhnCCqrrwM6HMzCrgBFULvs1nZjYqJ6gcRbnZzMEDJczMKuAEVQOe0dzMbHROULXgHpSZ2ahyTVCSLpG0XtIGSdeX2S5JN6TbH5Z0bmbbzZK2S3q0pM0sSXdLeir9eUxm20fSfa2X9OY8z+2QeEZzM7NR5ZagJDUCnwdWAMuAyyUtK6m2Aliavq4Gbsxs+zJwSZldXw/cExFLgXvSddJ9Xwacnrb7xzSGuuMZzc3MRpdnD+o8YENEbIqIHuA2YGVJnZXArZG4D5gpaT5ARPwI2FVmvyuBW9LlW4B3ZMpvi4juiHga2JDGUH96D9Q6AjOzupdngloAPJdZb0/Lxlqn1NyI6ABIfx43ln1JulrSWklrd+zYMepJTBQ1Txpc6XGCMjMbTZ4JSmXKSsddV1JnIo9HRNwUEcsjYvmcOXPGeahxaJk8GINv8ZmZjSrPBNUOLMqsLwS2jKNOqW39twHTn9sPYV/V09w2uNzTWbs4zMwOE3kmqDXAUklLJLWQDGBYXVJnNXBFOprvfGB3/+27EawGrkyXrwS+kym/TFKrpCUkAy9+OREnMhHU4lt8ZmZj0ZTXjiOiIOk64E6gEbg5ItZJuibdvgq4A7iUZEBDJ3BVf3tJXwMuAmZLagc+GhFfAj4J3C7pA8CzwLvT/a2TdDvwGFAAro2IvrzOryLZSWGH3OJzgjIzG01uCQogIu4gSULZslWZ5QCuHabt5cOU7wQuHmbbJ4BPjDfeXLkHZWY2Jp5JokrU7B6UmdlYOEFVS0t2kIQTlJnZaJygqkSZz6D8RV0zs9E5QVVLdpCEh5mbmY1qxAQl6cuZ5StHqGqjaWiChnRqwGIf4QljzcxGNFoP6qzM8h/mGcgRTxo63ZFnkzAzG9FoCWq80w5ZOUOGmvs2n5nZSEb7HtRCSTeQzHPXvzwgIj6YW2RHgijJ7yVf1i03eaCZmSVGS1D/T2Z5bZ6BHA3U3DbYJfVQczOzEY2YoCLilpG221jII/nMzMZgxAQlqXRy1yEi4u0TG86RTa1TB1e699UuEDOzw8Bot/h+jeQhgF8DfkH5Zy5ZpVqnDCxGlxOUmdlIRktQ84DfAC4Hfhv4HvC1iFiXd2BHImUSFD37axeImdlhYMRh5hHRFxH/GRFXAueTPBbjXkl/UJXojjSZW3zR7QRlZjaSUR+3IakVeAtJL2oxcAPwrXzDOgJJQ3tQvsVnZjai0QZJ3AKcAXwf+MuIeLQqUR2pfIvPzKxio/WgfgfYD5wK/KGk/q/xiOR5g9PzDO5II9/iMzOr2Gjfg/Js54ekZCaJtswwc9/iMzMb0Wi3+NqAa4BTgIeBmyOiUI3AjkjNk0k7n0TvAaLYh/pnODczsyFG6yHdAiwHHgEuBT6Te0RHMDU0DH1woWeTMDMb1mgJallE/LeI+Cfgt4ALxrJzSZdIWi9pg6Try2yXpBvS7Q9LOne0tpK+Lumh9LVZ0kNp+WJJBzLbVo0l1qrxl3XNzCoy2iCJ3v6FiChIlU8kIakR+DzJF33bgTWSVkfEY5lqK4Cl6evVwI3Aq0dqGxHvzRzjM8DuzP42RsTZFQdZA2qdQuxNVzxQwsxsWKMlqLMk7UmXBUxK1ysZxXcesCEiNgFIug1YCWQT1Erg1ogI4D5JMyXNJ/m+1YhtlWTL9wBvqOhMa6H0cRsw5Mu6HmpuZja80WaSaIyI6elrWkQ0ZZZHG2K+gGQev37taVkldSppewGwLSKeypQtkfSgpB9KKns7UtLVktZKWrtjx45RTmHiKTOSz7f4zMyGl+cw8nL3A0u7FMPVqaTt5SST2PbrAE6IiHOADwFflXRQEo2ImyJieUQsnzNnzrDB5ybbg+raO3w9M7Oj3KhTHR2CdmBRZn0hsKXCOi0jtZXUBLwLeGV/WUR0A93p8v2SNpJ8wbiuHrSoSYM5M7p2j1DTzOzolmcPag2wVNISSS3AZUDp86VWA1eko/nOB3ZHREcFbd8IPBER7f0FkuakgyuQdBLJwItNeZ3cuE2aMbAYB/aMUNHM7OiWWw8qHfV3HXAn0EjyJd91kq5Jt68C7iD5ftUGoBO4aqS2md1fxtDbewAXAh+XVAD6gGsiYlde5zde2R4UB9yDMjMbTp63+IiIO0iSULZsVWY5gGsrbZvZ9v4yZd8EvnkI4VaHe1BmZhXxXHvVkn6HzD0oM7PKOEFVW8sUSOffi94uotBT44DMzOqTE1SVSUJt2V6Ub/OZmZXjBJWrMjNJAGSHmvs2n5lZWU5QNaBJMweWo8s9KDOzcpygasADJczMRucEVQvZW3ydTlBmZuU4QeWp3GzmgCbPHFxxD8rMrCwnqBrQlFkDy7G/7ia7MDOrC05Q1ZJ92KMTlJnZqJygakCTjxlYjv27iGFuBZqZHc2coHLU01Vgy8ZdbN38In2F4kC5Wiah5knJSl8vdPvBhWZmpZygcrT5se3s2dnJS9v388A9JU/+8G0+M7MROUHl6OlHtw8sP3DPxiHbPFDCzGxkTlA5mjSlZWA5Sp5in01QOEGZmR3ECSpHU2a0DbvNPSgzs5E5QeVo0tTmIevdB3oHV7IJat/OaoVkZnbYcILKkYbe1WNH++DEsJp67MBy7N1RrZDMzA4bTlBVtKN9cFojTZ87sBx7d/i7UGZmJZygqijbg6J1yuB3oQrdnpPPzKxErglK0iWS1kvaIOn6Mtsl6YZ0+8OSzh2traSPSXpe0kPp69LMto+k9ddLenOe5zYeHU+/OLAsCaYfN7Aee7aXa2JmdtTKLUFJagQ+D6wAlgGXS1pWUm0FsDR9XQ3cWGHbz0bE2enrjrTNMuAy4HTgEuAf0/3UTOlduy0bXxyyrmmZBOXPoczMhsizB3UesCEiNkVED3AbsLKkzkrg1kjcB8yUNL/CtqVWArdFRHdEPA1sSPdTQ0Mz1J5dnex96cDAurI9qL3bqhaVmdnhIM8EtQB4LrPenpZVUme0tteltwRvltQ/82olx0PS1ZLWSlq7Y0f1ey3ZXtSQHpRv8ZmZDZFnglKZstKhasPVGantjcDJwNlAB/CZMRyPiLgpIpZHxPI5c+aUaZKvjk2DX8od0oPa4x6UmVlWngmqHViUWV8IbKmwzrBtI2JbRPRFRBH4AoO38So5Xs09vyGToGbMH1iOPduIvkItQjIzq0t5Jqg1wFJJSyS1kAxgWF1SZzVwRTqa73xgd0R0jNQ2/Yyq3zuBRzP7ukxSq6QlJAMvfpnXyVWkzFebnnl88DtPapk0OOVRsY/Y69t8Zmb9mvLacUQUJF0H3Ak0AjdHxDpJ16TbVwF3AJeSDGjoBK4aqW26609LOpvk7X8z8Htpm3WSbgceAwrAtRHRl9f5VaTMl2/3vniAXVv3cez8aQBo5oKBufjixedh5vFVDdHMrF7llqAA0iHgd5SUrcosB3BtpW3T8t8Z4XifAD4x3nirZfNj2wcT1DHHw/OPABAvPQ+8qoaRmZnVD88kkaNs/2nqMZMHljevGxw9qJmDAw3jxbr7yMzMrGacoKrk+JOOGVje/NjgZ006JpOgXmqvakxmZvXMCapK5iyaQVNzMrHFC8/v4aUd+wHQ9Hmg5NcQ+3YSvV01i9HMrJ44QeUpM0iiubmRJWcMfu9p/f3J7Tw1Ng0dbr7zmerFZ2ZWx5ygquhlywdv561f+/zAcsPsxQPL8cLmKkZkZla/nKByVDrI/NRXDg4h3/irbRR6k1Hwmr1koLz4wtPVCM3MrO45QVWLxOzjpzFrXjK8vKerl02PJNMbac5gggonKDMzwAkqVyrpQ0nitFcvHFh/5CfPJuUz5kNjCwDR+RLROfSxHGZmRyMnqCpROpXtKy44caDs0Z89S6G3DzU00jB7sLy4w70oMzMnqDyVmYtv0anHcsxxUwHo2t/DUw92AKA5Jw022/ZkVcIzM6tnTlA5KpOfkMQrLhzsLT1wzyYAGua9fKCsuHV93qGZmdU9J6hq0eDjqs69eLC39Nh97eze2YnmnAwNyRd546UtRNfeqodoZlZPnKByVa4PBcctnMFJZ84FoFgssvbujai5lYbscHP3oszsKOcEVSUqed7vq1csHVj+5X8+RaHQh+a9bKAsOh6vVmhmZnXJCSpPQzpQQzPUsl9bxNSZkwDYs7OTh37wNA3zTxvYXmx/ZODBhmZmRyMnqByNlF6amhp57crBgRH3fmMdxVmLUcuUpO2B3Z6Xz8yOak5Q1aKDi86/dCmTprYCsLNjLw//+DkaFp45sL343ENVCs7MrP44QeVplDt0bZNbeO3bBz93uvPWhyjMPX1gvfjsQzkFZmZW/5ygqqRMBwqA173jtCGfRf3slw3Q2AxA7O6g+KIfYmhmRycnqGopHcaXap3UzJt+56yB9Xu/tYGuGYO9quKmX+QemplZPco1QUm6RNJ6SRskXV9muyTdkG5/WNK5o7WV9DeSnkjrf1vSzLR8saQDkh5KX6vyPLfKZO7xDZOgAF75xpM4/qRZAPT2FLjrR20DDzssbvoFUSzmGqWZWT3KLUFJagQ+D6wAlgGXS1pWUm0FsDR9XQ3cWEHbu4EzIuIVwJPARzL72xgRZ6eva/I5s4nX0NDAuz54Pg2Nya/jwcdb2bEjSVBxYDfR8VgtwzMzq4k8e1DnARsiYlNE9AC3AStL6qwEbo3EfcBMSfNHahsRd0VEIW1/H7CQejWG7zEtOHkWr3/PGUkzGvj5o9M5sLcbgL4n7s0jOjOzupZngloAPJdZb0/LKqlTSVuA3wW+n1lfIulBST+UdMF4A8/D8Df4Bl30ntNZuPRYAJ7cfQrtG16k0NNH8flHiN1b8w3QzKzO5Jmgyr0nl3YphqszaltJfw4UgK+kRR3ACRFxDvAh4KuSph8UlHS1pLWS1u7YsWOUU6iupqZG3veRC5kyvY19hWk8u2cuzz35AsVCkb7H//9ah2dmVlV5Jqh2YFFmfSGwpcI6I7aVdCXwVuB9kc4HFBHdEbEzXb4f2AicWhpURNwUEcsjYvmcOXPGeWqVKcmoFbWZOWcKl/3P19HQ0MD6PS+ju7OX557aSe+TP/GTds3sqJJngloDLJW0RFILcBmwuqTOauCKdDTf+cDuiOgYqa2kS4APA2+PiM7+HUmakw6uQNJJJAMvNuV4frk55ax5/OYfns/2ruPY1TOLA3u7ee7xrXSt+W6tQzMzq5qmvHYcEQVJ1wF3Ao3AzRGxTtI16fZVwB3ApcAGoBO4aqS26a4/B7QCdyvpldyXjti7EPi4pALQB1wTEbvyOr9K6BAmez33DSfRubeb+/91KxfNvZcDe7vZ/L1/Z8GSX2f6CSeOvgMzs8NcbgkKICLuIElC2bJVmeUArq20bVp+yjD1vwl881DizVOFd/iGeN3K0yCC7f/xGMe1baens5v7P/1JTrr2LznxtOMmPkgzszrimSTq3OvesYx5K6+if9zILJ7ljr/6Z3723Sf8OA4zO6I5QVWJxtOFSp35tguZc9FbaGxKfl3nzFjDXV/8CV/6i3t4acf+iQrRzKyuOEHlaCI7OPPffiUnnnsybZNbmNR4gFfP/gUbf9XBZ3//P7j339ZR6O2buIOZmdUBJ6jDhFomM+Xi/86Jy+Zw7PxpHD+pg9OmP05PVy933vIgf3/tf/Doz571bT8zO2LkOkjCMsZ/h29Aw4IzaDr9jcxpuIepx0yi5enH2FuYRnvnInZ27OUr/++PmLf4GF7/3jM487UnHNJtRTOzWnMPKkdDezMTkywaz30XDXOXMmlqC0vOOI53veoJFs4cHE2/dfOLfO1TP+Yzv7ean3zncTr3dU/Icc3Mqs0JqlomqDejxiaaLroGTZuDJGYd18ZVFz/OpW+dSktb80C9nR17+d4X7+eTV36bb95wHxsf3urbf2Z2WPEtviqZyLttap1K88UfpPfOvyUO7Kaxoch5U+7m3I9ezk8fnMrPv/ckXft7gOT5Umvv3sDauzcw/djJnHXhYs547QksOvVY3wI0s7rmBHWY0vTjaH7zn9B7198RnS9BsUDTA//CG85cwYW/tZKHf/Qs933vSTqeHpy/b8/OTn787cf48bcfY8qMNl62/HhOO28hp5w9j7bJLbU7GTOzMpygDmOaPpfmN/8pvfd8jtiTPI6j75Hv07j1SV71uqt41Zsu5dn1L/DQvZt55MfPsH9P10Db/bu7eOCeTTxwzyYkseCUWSw5Yy4nnTmXE5fNYdIUJywzqy0dzZ9LLF++PNauXZvb/tf+xR/T/OKTAEx524c45c0X53Kc6Omk8KMvUNySefJuUytNr3gLDaddjBqbKBT62PirbTz8482sX7uF/bu7ht2fJOaeOJMFp8xi4dJjWbj0WOYtnklTc2Mu8ZvZ0U3S/RGxvLTcPagqyfPzHrVMpukN19H36H/S96v/gChCoZvCA99CT/2ExnPfSeMJ5/CyVx7Py155PBFB+1M7eWLN8zyx5nk6Nr04ZABFRLB184ts3fwi9//XRgAamxqZe8IM5p44g+NOmMm8E2cy98QZzJwzxZ9lmVkunKCOEGpoTHpMx59O4adfJnZ3ABB7t1P44T+hGfNpPHMFDSe+EjU2sejU2Sw6dTa/8b6z6NzXzTOP7eDpR7fz9KPbeH7DroNG/PUV+tiyaRdbNg2dIL6lrZnjFk3n2OOncez8acyaN3Xg57RjJjl5mdm4OUFVSbXephtmL6b5rf+L4pM/pO+h7xK9BwCI3R0UfnIzWvsNGk55DY2nvA5NT2ZEnzy1ldPOW8hp5y0EoLurl+c37EpeT+2k/amd7OzYW/Z4PV29tKd1SrW0NTNr3lSOOW4KM2ZPHnhNP3bwZ0ur/wTNrDy/OxyB1NhE42kX07DkPPoeu5u+J+6FQvKF3ejaS9+jd9L36J3o2BNpPOFcGk48B02fO9C+ta2Zk86Yy0lnDJZ17utm6+aX2P7sbrY98xLbnt3N9md3Dxl4Uaqnq3fgVuFwJk9rZfqxk5l2TBtTZ7YxdeYkps5sY8qM1oHl/vWmJn8GZnY0cYLK0ZDbZDW41aW2aTSd+y4aT38TfU/cS/GpHydD0vvj2/kMhZ3PwIPfRtPn0jDv5Wj+y2mYdypqnTpkX5Onth6UtAD27e5i+7O72bV1Lzs79rGzYy+7tiY/+7+LNZLOvd107u1m6+bRz2fS1FamTG9l0tQWJk1rYdKUFiZNa2XytJakLH1NntbK5GmttE1pZtLUFpqaG32r0eww5AR1FFDrVJrOeitx5gri+Ufpe/LHFLesSwZTpGLPNvr2bIMnfwgIzZhHw7EnomNPRLMXo2MWoqaDh55PndHG1DPbOOnMuQdt69zbzc6Ovex+oTPz2s/uFzrZs/MAu3d2UuwrHtRuOAf2dXNgHFM3NTQ20Da5mdZJzbRObqZ1UlP6szkpz2xrG9jeQktbE82tjbS0NdHSOrjshGdWHU5QOaq3tzA1NKJFZ9Gw6Cyiez/F9kcoPvMAxY7HoK83UzOI3R307e6ATff1t0bTZqMZ89HM+Wj6PDRjHpo2B1qnln3D7u/JLDq1fDzFYpF9L3WxZ+cB9r3Uxb6XDrBvd/fA8v6B5S727+4a91RNxb7iQE9tIkiiqeXgxNXcv5z+bG5tSspbGmlqaaS5pZHG5kaamhsGlvu3NTY1DCw3pXX6l5vT7U6KdrRxgqqWOntzUesUGk8+n8aTzycKPcSOTRS3PkF0PEFx5zNDeleJIPbuIPbugPaHh25qbEZTZ6MpswZ+MnkGmjQD2qYlP1unooahUz82NDQwfdZkps+aPGq8xWKRzr09dO7t5sDeHg7s6+HAvm46B5YH1zv3dg+Ude3vpa8wsc/Kigh6uwv0dheo5uMim5rLJ7CpM9v4wF/l8x07s1pygqqWOktQWWpqST57mv9yOIckYe16lnhhM8WdzxAvbE4SE8P0YPp6id0dA0PbhzkKmjQN2qYnn2+1TkYtU5KfzZOhdUpJ2SRoboOmVmhspqGhIbmdOKNtzOdX6O2jq7OX7s5eug/00tXZS8+BQmY5+dl9oEB3Zy9dnT10d/bS01Wgt7uPnu7kZ293gZ6uQs0eDlno7St77Bmzp9QgGrP8OUHl6HCdpUNNLei4U+C4U+gfNxeFHmLPVuKlrQPJKPZsh/07id7hR/INCuLAHjiwZ7g0N1JEqLk1SVZp0krW26C5FTW1QmMLNDZBY3PSoxtYbkGNzUxqbGJyYzNMaoapTdDUAg2tqHEqNDRAQ1PmZyM0NA17S61YLNLb0zeYwLoKaY8qWe7pHlzv7S7Q091HX28fvT19FHr6KPQWk2TT05cmneJBy73pel9vcqyReoGNTX4ogR2Zck1Qki4B/j+gEfhiRHyyZLvS7ZcCncD7I+KBkdpKmgV8HVgMbAbeExEvpts+AnwA6AM+GBF35nl+Y1G//afKqKkFzToBZp1w0Lbo6ST2vQD7dhH7XiD27yK69sKB3YNJqedQboZFkgR7u+DA7rSkGpQmq0aU/uxPZGpoolUNtDY0ocZGUONAXdSQvNoaYLKSbWoACfVvkwbrNTQMLg/Z1pRZFn1FUeyDvr6grwiFQlDsEzSIvg0/TeKV0p8Ng712ZcvTGPrLIak70JYhdQf3xWBZdvuQfQzuS+X21V8HDf0for/OkPXM8bK/jyH/I5Xbnml70P5GOp4OquvP/GovtwQlqRH4PPAbQDuwRtLqiMhMGMcKYGn6ejVwI/DqUdpeD9wTEZ+UdH26/mFJy4DLgNOB44H/knRqREz4/ZiuPXt55O8+MWq9hn0j3fI6cqhl8rDJq1/0FaBrT5KwuvcTPZ3Qs5/o7oSe5BXd+wfLerug0J0kpmKhimczJOrk2MXCiAmxmv3khvTVXFJe2FTFII46wyS7YZPjwesaklxLk2HJPitJjKXHGuv2iusMs1ImxpZ3/nX5fRyCPHtQ5wEbImITgKTbgJVANkGtBG6N5F7YfZJmSppP0jsaru1K4KK0/S3AvcCH0/LbIqIbeFrShjSGn0/0iRV7eml54ZExtVHD0f2vMTU2wZRZyQCKMYq+AvT1pEmrJ0lahe40gXUn5X29yWdhfb1JUin0QF8hU9YLhd5MvcLAMsVCMiikr0AU+yD6oFibz5msHmX+CZK9bT+Gf5kcnjf7ay/PBLUAeC6z3k7SSxqtzoJR2s6NiA6AiOiQdFxmX/eVtFlQGpSkq4GrAU44Yfh/8U+kvoY25r3ynKoc60iUfJ7UBC3JaL9qpPqISJJWsS/zKkCxOJDQolgo2Z6+og8iiCim+ygmb2z96xFpnZLytN5AuyGvbL04eHsSNEGkdQMo87OYvlUOtBlsW7Z+mX0NfLbav62SfUXJcZOVoW/4B62TaTdSvUj/K5dIStuXqVO6L6sbeSaocu8jpb/94epU0nY8xyMibgJuguRxG6Pss6zmaVOY/Obfr6yyGljwmvNpmz5tPIeyGpE0+LnScHWqGI9V35BBTkOS7kBhZhsldTN1+tuW7qtkNwe9XQ2XOIcdfDVMPGV3P1ydCT7mIcozQbUDizLrC4EtFdZpGaHtNknz097TfGD7GI43IZpbW1n6trfksWszqxMa9vMmq5Y8x6euAZZKWiKphWQAw+qSOquBK5Q4H9id3r4bqe1q4Mp0+UrgO5nyyyS1SlpCMvDil3mdnJmZ5Su3HlREFCRdB9xJMlT85ohYJ+madPsq4A6SIeYbSIaZXzVS23TXnwRul/QB4Fng3WmbdZJuJxlIUQCuzWMEn5mZVYcf+Z7jI9/NzGx0wz3y3V9BNzOzuuQEZWZmdckJyszM6pITlJmZ1aWjepCEpB3AM4ewi9nACxMUTl7qPUbHd2jqPT6o/xgd36GZiPhOjIg5pYVHdYI6VJLWlht5Uk/qPUbHd2jqPT6o/xgd36HJMz7f4jMzs7rkBGVmZnXJCerQ3FTrACpQ7zE6vkNT7/FB/cfo+A5NbvH5MygzM6tL7kGZmVldcoIyM7O65AQ1TpIukbRe0gZJ19cohkWSfiDpcUnrJP1hWv4xSc9Leih9XZpp85E05vWS3lyFGDdLeiSNY21aNkvS3ZKeSn8eU4v4JL0sc40ekrRH0h/V+vpJulnSdkmPZsrGfM0kvTK99hsk3SBNzEONhonvbyQ9IelhSd+WNDMtXyzpQOZarqpRfGP+neYV3wgxfj0T32ZJD6XlVb2GI7yvVP9vMCL8GuOL5BEgG4GTSB6u+CtgWQ3imA+cmy5PA54ElgEfA/60TP1laaytwJL0HBpzjnEzMLuk7NPA9eny9cCnahVfye90K3Bira8fcCFwLvDooVwzkueh/RrJw3+/D6zIMb43AU3p8qcy8S3O1ivZTzXjG/PvNK/4houxZPtngP9di2vI8O8rVf8bdA9qfM4DNkTEpojoAW4DVlY7iIjoiIgH0uW9wOPAghGarARui4juiHia5Dlc5+Ufadk4bkmXbwHekSmvVXwXAxsjYqSZRaoSX0T8CNhV5tgVXzMlT5ueHhE/j+Sd4tZMmwmPLyLuiohCunofyROth1Xt+EZQ9es3WoxpL+M9wNdG2kdeMY7wvlL1v0EnqPFZADyXWW9n5MSQO0mLgXOAX6RF16W3W27OdMVrEXcAd0m6X9LVadncSJ6cTPrzuBrG1+8yhr4h1Mv16zfWa7YgXS4tr4bfJfnXcr8lkh6U9ENJF6RltYhvLL/TWl6/C4BtEfFUpqwm17DkfaXqf4NOUONT7j5qzcbrS5oKfBP4o4jYA9wInAycDXSQ3C6A2sT92og4F1gBXCvpwhHq1uS6SmoB3g58Iy2qp+s3muFiqtW1/HOSJ1p/JS3qAE6IiHOADwFflTS9BvGN9Xday9/15Qz9x1JNrmGZ95Vhqw4TxyHH5wQ1Pu3Aosz6QmBLLQKR1EzyR/SViPgWQERsi4i+iCgCX2DwNlTV446ILenP7cC301i2pd3//tsU22sVX2oF8EBEbEtjrZvrlzHWa9bO0Ntsuccq6UrgrcD70ls6pLd9dqbL95N8PnFqteMbx++06tcPQFIT8C7g6/1ltbiG5d5XqMHfoBPU+KwBlkpakv7r+zJgdbWDSO9Vfwl4PCL+LlM+P1PtnUD/SKHVwGWSWiUtAZaSfIiZV3xTJE3rXyb5IP3RNI4r02pXAt+pRXwZQ/7FWi/Xr8SYrll6C2avpPPTv5MrMm0mnKRLgA8Db4+Izkz5HEmN6fJJaXybahDfmH6n1Y4v443AExExcGus2tdwuPcVavE3eKgjPo7WF3ApyeiWjcCf1yiG15F0mR8GHkpflwL/AjySlq8G5mfa/Hka83omcFTSMPGdRDK651fAuv7rBBwL3AM8lf6cVYv40uNNBnYCMzJlNb1+JMmyA+gl+VfoB8ZzzYDlJG/EG4HPkc4ck1N8G0g+h+j/O1yV1v3N9Hf/K+AB4G01im/Mv9O84hsuxrT8y8A1JXWreg0Z/n2l6n+DnurIzMzqkm/xmZlZXXKCMjOzuuQEZWZmdckJyszM6pITlJmZ1SUnKLNhSApJn8ms/6mkj03Qvr8s6bcmYl+jHOfd6azUPygpXyzpt8e5z59NTHRmI3OCMhteN/AuSbNrHUhW/5c2K/QB4Pcj4vUl5YuBsgkqnc1gWBHxmjEc32zcnKDMhlcAbgL+uHRDaQ9I0r7050XphJ63S3pS0iclvU/SL9Pn4pyc2c0bJf04rffWtH2jkmcrrUknNv29zH5/IOmrJF84LY3n8nT/j0r6VFr2v0m+dLlK0t+UNPkkcIGS5wv9saT3S/qGpO+STO47VdI9kh5I97syc6zsud4r6d+UPAvqK+mMAaTn/Vh6Dn875itvBoz4LyUz4/PAw5I+PYY2ZwGnkTxOYRPwxYg4T8mD3/4A+KO03mLg10kmMf2BpFNIpoPZHRGvktQK/FTSXWn984AzInmkwQBJx5M8g+mVwIskCeYdEfFxSW8geQ7S2pIYr0/L+xPj+0me2/OKiNiV9qLeGRF70h7kfZJWx8Hf7D8HOJ1kjrWfAq+V9BjJdEIvj4hQ+vBCs7FyD8psBJHM4nwr8MExNFsTyTN1ukmmeOlPMI+QJKV+t0dEMZLHKmwCXk4yX+EVSp6m+guS6WWWpvV/WZqcUq8C7o2IHZE8k+krJA/EG6u7I6L/GUUC/o+kh4H/InlMwtwybX4ZEe2RTML6UHp+e4Au4IuS3gV0lmlnNionKLPR/T3JZzlTMmUF0v9/0ttaLZlt3ZnlYma9yNC7FqW9kf5HFPxBRJydvpZERH+C2z9MfBP1KPLs/t8HzAFeGRFnA9uAtjJtsufaR/JU3QJJb++bJA+o+88Jis+OMk5QZqNIexW3kySpfptJbqlB8kTR5nHs+t2SGtLPpU4imWjzTuB/KHncAZJOTWeCH8kvgF+XNDsdQHE58MNR2uwleZz3cGYA2yOiV9LrgRMrOB/SmKeSTL57B8ntzLMrbWuW5c+gzCrzGeC6zPoXgO9I+iXJzM7D9W5Gsp4kkcwlmcG6S9IXSW6TPZD2zHYwymOyI6JD0keAH5D0pu6IiNEea/AwUJD0K5IZtF8s2f4V4LuS1pLcunui8tNiGsm1aUvjOWiQiVklPJu5mZnVJd/iMzOzuuQEZWZmdckJyszM6pITlJmZ1SUnKDMzq0tOUGZmVpecoMzMrC79X2a0grRAEchsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "uniform.plot(color='C4')\n",
    "power.plot(color='C1')\n",
    "\n",
    "decorate(xlabel='Number of trains',\n",
    "         ylabel='PMF',\n",
    "         title='Posterior distributions')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The power law gives less prior probability to high values, which yields lower posterior means, and less sensitivity to the upper bound.\n",
    "\n",
    "Here's how the posterior means depend on the upper bound when we use a power law prior and observe three trains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.310454Z",
     "iopub.status.busy": "2021-04-16T19:35:21.308708Z",
     "iopub.status.idle": "2021-04-16T19:35:21.322722Z",
     "shell.execute_reply": "2021-04-16T19:35:21.323129Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "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>Posterior mean</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Upper bound</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>500</th>\n",
       "      <td>130.708470</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1000</th>\n",
       "      <td>133.275231</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000</th>\n",
       "      <td>133.997463</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             Posterior mean\n",
       "Upper bound                \n",
       "500              130.708470\n",
       "1000             133.275231\n",
       "2000             133.997463"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame(columns=['Posterior mean'])\n",
    "df.index.name = 'Upper bound'\n",
    "\n",
    "alpha = 1.0\n",
    "dataset = [30, 60, 90]\n",
    "\n",
    "for high in [500, 1000, 2000]:\n",
    "    hypos = np.arange(1, high+1)\n",
    "    ps = hypos**(-alpha)\n",
    "    power = Pmf(ps, hypos)\n",
    "    for data in dataset:\n",
    "        update_train(power, data)\n",
    "    df.loc[high] = power.mean()\n",
    "    \n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the differences are much smaller.  In fact,\n",
    "with an arbitrarily large upper bound, the mean converges on 134.\n",
    "\n",
    "So the power law prior is more realistic, because it is based on\n",
    "general information about the size of companies, and it behaves better in practice."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Credible Intervals\n",
    "\n",
    "So far we have seen two ways to summarize a posterior distribution: the value with the highest posterior probability (the MAP) and the posterior mean.\n",
    "These are both **point estimates**, that is, single values that estimate the quantity we are interested in.\n",
    "\n",
    "Another way to summarize a posterior distribution is with percentiles.\n",
    "If you have taken a standardized test, you might be familiar with percentiles.\n",
    "For example, if your score is the 90th percentile, that means you did as well as or better than 90\\% of the people who took the test.\n",
    "\n",
    "If we are given a value, `x`, we can compute its **percentile rank** by finding all values less than or equal to `x` and adding up their probabilities.\n",
    "\n",
    "`Pmf` provides a method that does this computation.\n",
    "So, for example, we can compute the probability that the company has less than or equal to 100 trains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.326884Z",
     "iopub.status.busy": "2021-04-16T19:35:21.326371Z",
     "iopub.status.idle": "2021-04-16T19:35:21.329117Z",
     "shell.execute_reply": "2021-04-16T19:35:21.328729Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2937469222495771"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "power.prob_le(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a power law prior and a dataset of three trains, the result is about 29%.\n",
    "So 100 trains is the 29th percentile.\n",
    "\n",
    "Going the other way, suppose we want to compute a particular percentile; for example, the median of a distribution is the 50th percentile.\n",
    "We can compute it by adding up probabilities until the total exceeds 0.5.\n",
    "Here's a function that does it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.332553Z",
     "iopub.status.busy": "2021-04-16T19:35:21.332041Z",
     "iopub.status.idle": "2021-04-16T19:35:21.334402Z",
     "shell.execute_reply": "2021-04-16T19:35:21.333955Z"
    }
   },
   "outputs": [],
   "source": [
    "def quantile(pmf, prob):\n",
    "    \"\"\"Compute a quantile with the given prob.\"\"\"\n",
    "    total = 0\n",
    "    for q, p in pmf.items():\n",
    "        total += p\n",
    "        if total >= prob:\n",
    "            return q\n",
    "    return np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The loop uses `items`, which iterates the quantities and probabilities in the distribution.\n",
    "Inside the loop we add up the probabilities of the quantities in order.\n",
    "When the total equals or exceeds `prob`, we return the corresponding quantity.\n",
    "\n",
    "This function is called `quantile` because it computes a quantile rather than a percentile.\n",
    "The difference is the way we specify `prob`.\n",
    "If `prob` is a percentage between 0 and 100, we call the corresponding quantity a percentile.\n",
    "If `prob` is a probability between 0 and 1, we call the corresponding quantity a **quantile**.\n",
    "\n",
    "Here's how we can use this function to compute the 50th percentile of the posterior distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.337871Z",
     "iopub.status.busy": "2021-04-16T19:35:21.337255Z",
     "iopub.status.idle": "2021-04-16T19:35:21.340128Z",
     "shell.execute_reply": "2021-04-16T19:35:21.339696Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "113"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "quantile(power, 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result, 113 trains, is the median of the posterior distribution.\n",
    "\n",
    "`Pmf` provides a method called `quantile` that does the same thing.\n",
    "We can call it like this to compute the 5th and 95th percentiles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.344566Z",
     "iopub.status.busy": "2021-04-16T19:35:21.343907Z",
     "iopub.status.idle": "2021-04-16T19:35:21.346263Z",
     "shell.execute_reply": "2021-04-16T19:35:21.346610Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 91., 243.])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "power.quantile([0.05, 0.95])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is the interval from 91 to 243 trains, which implies:\n",
    "\n",
    "* The probability is 5% that the number of trains is less than or equal to 91.\n",
    "\n",
    "* The probability is 5% that the number of trains is greater than 243.\n",
    "\n",
    "Therefore the probability is 90% that the number of trains falls between 91 and 243 (excluding 91 and including 243).\n",
    "For this reason, this interval is called a 90% **credible interval**.\n",
    "\n",
    "`Pmf` also provides `credible_interval`, which computes an interval that contains the given probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.350591Z",
     "iopub.status.busy": "2021-04-16T19:35:21.349787Z",
     "iopub.status.idle": "2021-04-16T19:35:21.353876Z",
     "shell.execute_reply": "2021-04-16T19:35:21.353340Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 91., 243.])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "power.credible_interval(0.9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The German Tank Problem\n",
    "\n",
    "During World War II, the Economic Warfare Division of the American\n",
    "Embassy in London used statistical analysis to estimate German\n",
    "production of tanks and other equipment.\n",
    "\n",
    "The Western Allies had captured log books, inventories, and repair\n",
    "records that included chassis and engine serial numbers for individual\n",
    "tanks.\n",
    "\n",
    "Analysis of these records indicated that serial numbers were allocated\n",
    "by manufacturer and tank type in blocks of 100 numbers, that numbers\n",
    "in each block were used sequentially, and that not all numbers in each\n",
    "block were used.  So the problem of estimating German tank production\n",
    "could be reduced, within each block of 100 numbers, to a form of the\n",
    "train problem.\n",
    "\n",
    "Based on this insight, American and British analysts produced\n",
    "estimates substantially lower than estimates from other forms\n",
    "of intelligence.  And after the war, records indicated that they were\n",
    "substantially more accurate.\n",
    "\n",
    "They performed similar analyses for tires, trucks, rockets, and other\n",
    "equipment, yielding accurate and actionable economic intelligence.\n",
    "\n",
    "The German tank problem is historically interesting; it is also a nice\n",
    "example of real-world application of statistical estimation.\n",
    "\n",
    "For more on this problem, see [this Wikipedia page](https://en.wikipedia.org/wiki/German_tank_problem) and Ruggles and Brodie, \"An Empirical Approach to Economic Intelligence in World War II\", *Journal of the American Statistical Association*, March 1947, [available here](https://web.archive.org/web/20170123132042/https://www.cia.gov/library/readingroom/docs/CIA-RDP79R01001A001300010013-3.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Informative Priors\n",
    "\n",
    "Among Bayesians, there are two approaches to choosing prior\n",
    "distributions.  Some recommend choosing the prior that best represents\n",
    "background information about the problem; in that case the prior\n",
    "is said to be **informative**.  The problem with using an informative\n",
    "prior is that people might have different information or\n",
    "interpret it differently.  So informative priors might seem arbitrary.\n",
    "\n",
    "The alternative is a so-called **uninformative prior**, which is\n",
    "intended to be as unrestricted as possible, in order to let the data\n",
    "speak for itself.  In some cases you can identify a unique prior\n",
    "that has some desirable property, like representing minimal prior\n",
    "information about the estimated quantity.\n",
    "\n",
    "Uninformative priors are appealing because they seem more\n",
    "objective.  But I am generally in favor of using informative priors.\n",
    "Why?  First, Bayesian analysis is always based on\n",
    "modeling decisions.  Choosing the prior is one of those decisions, but\n",
    "it is not the only one, and it might not even be the most subjective.\n",
    "So even if an uninformative prior is more objective, the entire analysis is still subjective.\n",
    "\n",
    "Also, for most practical problems, you are likely to be in one of two\n",
    "situations: either you have a lot of data or not very much.  If you have a lot of data, the choice of the prior doesn't matter;\n",
    "informative and uninformative priors yield almost the same results.\n",
    "If you don't have much data, using relevant background information (like the power law distribution) makes a big difference.\n",
    "\n",
    "And if, as in the German tank problem, you have to make life and death\n",
    "decisions based on your results, you should probably use all of the\n",
    "information at your disposal, rather than maintaining the illusion of\n",
    "objectivity by pretending to know less than you do."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "This chapter introduces the train problem, which turns out to have the same likelihood function as the dice problem, and which can be applied to the German Tank problem.\n",
    "In all of these examples, the goal is to estimate a count, or the size of a population.\n",
    "\n",
    "In the next chapter, I'll introduce \"odds\" as an alternative to probabilities, and Bayes's Rule as an alternative form of Bayes's Theorem.\n",
    "We'll compute distributions of sums and products, and use them to estimate the number of Members of Congress who are corrupt, among other problems.\n",
    "\n",
    "But first, you might want to work on these exercises."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Suppose you are giving a talk in a large lecture hall and the fire marshal interrupts because they think the audience exceeds 1200 people, which is the safe capacity of the room. \n",
    "\n",
    "You think there are fewer then 1200 people, and you offer to prove it.\n",
    "It would take too long to count, so you try an experiment:\n",
    "\n",
    "* You ask how many people were born on May 11 and two people raise their hands.  \n",
    "\n",
    "* You ask how many were born on May 23 and 1 person raises their hand.  \n",
    "* Finally, you ask how many were born on August 1, and no one raises their hand.\n",
    "\n",
    "How many people are in the audience?  What is the probability that there are more than 1200 people.\n",
    "Hint: Remember the binomial distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.359329Z",
     "iopub.status.busy": "2021-04-16T19:35:21.358751Z",
     "iopub.status.idle": "2021-04-16T19:35:21.363183Z",
     "shell.execute_reply": "2021-04-16T19:35:21.363696Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "200"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# I'll use a uniform prior from 1 to 2000\n",
    "# (we'll see that the probability is small that there are\n",
    "# more than 2000 people in the room)\n",
    "\n",
    "hypos = np.arange(1, 2000, 10)\n",
    "prior = Pmf(1, hypos)\n",
    "prior.normalize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.370055Z",
     "iopub.status.busy": "2021-04-16T19:35:21.369426Z",
     "iopub.status.idle": "2021-04-16T19:35:21.371847Z",
     "shell.execute_reply": "2021-04-16T19:35:21.372439Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# We can use the binomial distribution to compute the probability\n",
    "# of the data for each hypothetical audience size\n",
    "\n",
    "from scipy.stats import binom\n",
    "\n",
    "likelihood1 = binom.pmf(2, hypos, 1/365)\n",
    "likelihood2 = binom.pmf(1, hypos, 1/365)\n",
    "likelihood3 = binom.pmf(0, hypos, 1/365)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.380246Z",
     "iopub.status.busy": "2021-04-16T19:35:21.379424Z",
     "iopub.status.idle": "2021-04-16T19:35:21.382846Z",
     "shell.execute_reply": "2021-04-16T19:35:21.383434Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.006758799800451805"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# Here's the update\n",
    "\n",
    "posterior = prior * likelihood1 * likelihood2 * likelihood3\n",
    "posterior.normalize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.406019Z",
     "iopub.status.busy": "2021-04-16T19:35:21.400248Z",
     "iopub.status.idle": "2021-04-16T19:35:21.545515Z",
     "shell.execute_reply": "2021-04-16T19:35:21.545140Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+zUlEQVR4nO3deXxV9bno/8+TGQghBMIYIAHCDAIGcKiKihYsStsjFmuFVnssVU+nV++t3t7T6z3n9Bx7vHbwZ49WW6ceK2qtihZHnAVknscwJhAgDAFCyPz8/lgrOyubzGTvtffO83699itr+q717JWd/eS71nd9v6KqGGOMMZEmzu8AjDHGmMZYgjLGGBORLEEZY4yJSJagjDHGRCRLUMYYYyJSgt8B+Kl3796anZ3tdxjGGNOprVmz5piqZgYv79QJKjs7m9WrV/sdhjHGdGoisr+x5XaJzxhjTESyBGWMMSYiWYIyxhgTkTr1PShjjGmLqqoqCgsLKS8v9zuUqJSSkkJWVhaJiYmt2t4SlDHGtFJhYSHdu3cnOzsbEfE7nKiiqhw/fpzCwkJycnJaVcYu8RljTCuVl5fTq1cvS07tICL06tWrTbVPS1DGGNMGlpzar63nzi7xmWYdLzrDge3F1NYqXVKTyBnXly7dkvwOyxjTCViCMo3aseYQbz+zjsP7TjZYHhcfx8i8AcxcMIk+g3r4FJ0x5kK89tprjBgxgjFjxrSp3OLFi9m6dSv33XdfiCJryBKUaaC6uoa/PbKCdR/ubXR9bU0t274oZMfqQ1x9yziuvXW8XfIwJsq89tprzJ49u00Jqrq6mptuuombbrqpTWUSEtqfZixBmYCK8ir+8h+fsnPtocCyxKQEcsb3IbVHCkcLTlG46zjgJKqlL2zkyIES5v74MpKS7aNkTDjs27ePmTNnMm3aNNatW8eIESN47rnnWL58OT/96U+prq5mypQpPPbYYyQnJ3PfffexePFiEhISuP766/n617/O4sWL+fjjj/m3f/s3XnnlFQDuueceiouL6dq1K08++SSjRo3i29/+NhkZGaxbt47Jkyczfvx4Vq9ezaOPPsr+/fu54447KC4uJjMzk6effprBgwefV+bhhx9u93u1bxUDOE1AX/71sgbJaeL0HG64czLd07sElhXtPcnrj61k/7ZiADZ/foDK8mpu/99XkZAQH/a4jfHL/Tf+d8j2/R9vfKvZ9Tt27OBPf/oTl19+OXfccQe//vWv+cMf/sDSpUsZMWIE8+fP57HHHmP+/Pm8+uqrbN++HRGhpKSE9PR0brrpJmbPns3NN98MwLXXXsvjjz9Obm4uX3zxBXfffTcffPABADt37uT9998nPj6eZ555JhDDvffey/z581mwYAFPPfUUP/jBD3jttdfOK3MhrBWfAeDjV7ayZXlBYH7GNydwy08ua5CcAPrn9OQf/+M6LrtxVGDZzjWH+OtvlqOqYYvXmM5s0KBBXH755QB861vfYunSpeTk5DBixAgAFixYwCeffEJaWhopKSl897vf5W9/+xtdu3Y9b1+lpaUsW7aMuXPnMnHiRL73ve9RVFQUWD937txGE83y5cv55je/CcDtt9/OZ5991mKZtrIalKFg5zHefW59YP6yG0dx7a0Tmtw+Pj6OG+/KI6VbIh8s2gTAhk/20S87nelzx4U6XGM6vdbe901ISGDlypUsXbqURYsW8eijjwZqRnVqa2tJT09n/fr1je6jW7dubY6ptWVaYgmqk6upqeW1368M1H6yx/Rh1h2TWlV2xjcnUHa6ghVLdgLw7p83MGhkb4ZN6BeyeI2JFC1dhgulAwcOsHz5ci699FJeeOEFZsyYwR/+8Afy8/MZPnw4f/7zn7nqqqsoLS2lrKyMG264gUsuuYThw4cD0L17d86cOQNAWloaOTk5vPzyy8ydOxdVZePGjVx00UXNxnDZZZexaNEibr/9dp5//nm+9KUvdfj7tEt8ndwXS3ZyaM8JwGkQcfOPL231vSQR4Sv/eDHZY/oAzn2slx5eRllpRcjiNcbA6NGjefbZZ5kwYQInTpzgxz/+MU8//TRz585l/PjxxMXFsXDhQs6cOcPs2bOZMGECV111Fb/5zW8AmDdvHg899BCTJk1i9+7dPP/88/zpT3/ioosuYuzYsbz++ustxvDII4/w9NNPM2HCBP785z/zu9/9rsPfp3Tm+wZ5eXnamQcsLC+r5KHvvk7ZGSehXH/7RK6+pe2X6E6fKOORHyzh7CmnC5PJ1w5l7o8u69BYjYkE27ZtY/To0b7GsG/fPmbPns3mzZt9jaO9GjuHIrJGVfOCtw1pDUpEZorIDhHJF5HznuwSxyPu+o0iMtmz7ikROSoim4PKvCgi693XPhFZ7y7PFpFznnWPh/K9xYLlb+4MJKeefVO54mvt+8NLy+jKV++eGphfu3RPg9aAxhjTHiFLUCISD/wemAWMAW4VkeCnwmYBue7rLuAxz7pngJnB+1XVb6jqRFWdCLwC/M2zenfdOlVd2FHvJRZVnKvis9e2BeavmTeehMT2t7oZd9lgJlwxJDD/xh9WU11Vc0ExGmPOl52dHbW1p7YKZQ1qKpCvqntUtRJYBMwJ2mYO8Jw6VgDpItIfQFU/AU40tXNxmozcArwQkuhj3Bdv7QrUnjL6dWfi9OwL3ueN35tCittP37FDp/ns9e0XvE9jIk1nvi1yodp67kKZoAYCBZ75QndZW7dpyhXAEVXd5VmWIyLrRORjEbmisUIicpeIrBaR1cXFxa08VGypra1lxd93Buanzx3bIQ/ZpvZI4brb6punf/jiZk6fKLvg/RoTKVJSUjh+/LglqXaoGw8qJSWl1WVC2cy8sYb6wb/V1mzTlFtpWHsqAgar6nERuRh4TUTGqurpBjtXfQJ4ApxGEq08VkzZsfoQJ4+WAtC1e3KH1J7qTLthBCvfyefI/hIqy6v48MXNzPn+1JYLGhMFsrKyKCwspLP+c3uh6kbUba1QJqhCYJBnPgsIvnPemm3OIyIJwNeBi+uWqWoFUOFOrxGR3cAIoPM202uCt/aUd/1wEpM67mMQHx/HrO9M5pkHnIcBV76dz2U3jSJzYFqHHcMYvyQmJrZ6NFhz4UJ5iW8VkCsiOSKSBMwDFgdtsxiY77bmuwQ4papFwTtqxAxgu6oW1i0QkUy3YQYiMhSn4cWejngjseTEkdJACzsRYdqs3A4/xojJ/QMP69bW1vL+8xs6/BjGmNgXsgSlqtXAvcA7wDbgJVXdIiILRaSuhd0SnCSSDzwJ3F1XXkReAJYDI0WkUETu9Ox+Huc3jrgS2CgiG4C/AgtVtclGFp3V+o/qh9HIndyfjL6pHX4MEeHLCyYG5jd9doDig6ebLmCMMY0IaVdHqroEJwl5lz3umVbgnibK3trMfr/dyLJXcJqdmyaoaoNxniZfMzRkxxo0ojcjLh7AzjWHUFU+enmzPbxrjGkT6+qoEzmYf4Jjbk0muUsio6e1/mZle1zzjfGB6fUf7Qs0zDDGmNawBNWJeGtPYy8bFPJBBoeMziRnXF/AGeDwk1e2hvR4xpjYYgmqk1BVNn9+IDA/6erQXd7zuvob9X37rX5vtz0XZYxpNUtQnUTBzuOB5NAtLYWh4/uE5bjDL+rHoBG9AaiuquGz16x3CWNM61iC6iS8o+WOnpZFXFx4fvUiwvRbxgbmv3hrF+fOVobl2MaY6GYJqhNQVbYsq7+8N/bSQc1s3fFGT82i7+B0ACrLq1jz/u6wHt8YE50sQXUCRw6c4niRM3pmUkoiwy7qG9bjiwiX3jgyML/sjR3U1taGNQZjTPSxBNUJbF91MDA9Mm9Ah3Zt1FqTrs6ha/dkAE4eKW0QkzHGNMYSVCewc01994ajprS2s/iOlZScwJQvDw/ML1u8w5c4jDHRwxJUjCsvq2T/tvqel3Mn9fctlktuGBFonLF742EO7y/xLRZjTOSzBBXjdm84TG2Nc79nwNAMuvfs4lss6ZndGHNpfe8Vy96wWpQxpmmWoGLczrX1ncOPzBvgYySOy28aFZhe/+HewKi+xhgTzBJUDFPVBvefcif7n6CGjM5kwNAMAKoqq1n7gY2IYoxpnCWoGHa86AwlxWcBp3PYQSN7+RyR0+T8kq+MCMyvfDvfhs82xjTKElQM273xSGA6Z1wfEhLifYym3oQrh5CUkghAceEp9m214bONMeezBBXD9m6qT1BD3RFuI0FySiITp2cH5le+vcu/YIwxEcsSVIxS1QY1qGETwtt7REu8Q81v/vyANZYwxpzHElSMKi48TWnJOQC6pCbTP6enzxE1NGBoBlm5zj2x6qoaayxhjDmPJagYtXvj4cD00PF9EBEfo2nc1Jn1tShrLGGMCWYJKkbt3Xw0MD10fGRd3qsT3Fhi75ajLZQwxnQmIU1QIjJTRHaISL6I3NfIehGRR9z1G0VksmfdUyJyVEQ2B5V5QEQOish693WDZ9397r52iMiXQ/neIpmqst/TMq5u2PVIE9xYwobhMMZ4hSxBiUg88HtgFjAGuFVExgRtNgvIdV93AY951j0DzGxi979R1Ynua4l7vDHAPGCsW+6/3Bg6nRNHSgOj56Z0TaLvkB4+R9S0KdfXdyC76bMDVJyr8jEaY0wkCWUNaiqQr6p7VLUSWATMCdpmDvCcOlYA6SLSH0BVPwFOtOF4c4BFqlqhqnuBfDeGTueAp3PYwaN7h2303PYYODwjMJhhVUU1mz4/0HwBY0ynEcpvroFAgWe+0F3W1m0ac697SfApEalrntaqfYnIXSKyWkRWFxfH5gOi3gdfh4zO9DGSlokIF88YGpi3y3zGmDqhTFCNNRsLbqbVmm2CPQYMAyYCRcDDbdmXqj6hqnmqmpeZGdlf3u21z9PYIHtMHx8jaZ2J03MCtbx9W44GRv81xnRuoUxQhcAgz3wWcKgd2zSgqkdUtUZVa4Enqb+M1+Z9xaKyMxUcLTgFQFx8HFkj/O9/ryXde3ZhxMX141StWWrPRBljQpugVgG5IpIjIkk4DRgWB22zGJjvtua7BDilqkXBO/Kqu0fl+hpQ18pvMTBPRJJFJAen4cXKjngj0eTAjmOB6YHDMkhKDv/w7u1x8Yxhgem1S/fYM1HGmNAlKFWtBu4F3gG2AS+p6hYRWSgiC93NlgB7cBo0PAncXVdeRF4AlgMjRaRQRO50V/2niGwSkY3A1cCP3eNtAV4CtgJvA/eoak2o3l+kKvAkqMGjevsYSduMmjqQrt2TATh17Cz5Gw63UMIYE+tC+u+12wR8SdCyxz3TCtzTRNlbm1h+ezPH+yXwy3YFGyO8CWrQyOhJUAkJ8UycnsOyN7YDTi0qd6J/w9MbY/wXue2PTZupKgU7jwfmo+H+k5e3Nd+WZQWcO1vpYzTGGL9Zgoohxw6dodz9Uu+WlkJG31SfI2qbAUMzAp3aVlVWs+mz/T5HZIzxkyWoGNLw8l6viOwgtiWTr/U+E2Wt+YzpzCxBxRBvgsoaET33n7wmTs8hLt75WB7YXkzxwdM+R2SM8YslqBjivf8UTS34vFJ7pDBqSn0HIOtsnChjOi1LUDGiurqGw/tKAvMDh2f4F8wF8l7mW/vBXnsmyphOyhJUjDh64BQ11c5jXz37ptI1NdnniNpvZN6ABs9E7dl0pIUSxphYZAkqRhzaXd/x+8Bh0Vt7AueZqIuuyg7M23DwxnROlqBixMHdJwPTA6I8QQFMvqb+Mt/mzwuoKLdxoozpbCxBxQhvDSoWEtTA4RlkZjkDLVaWV7F1eUELJYwxscYSVAyora2laK+3BtWzma2jg4gw+ZqcwPzaD/b6GI0xxg+WoGLAsYNnqKqoBiAtoyvd07v4HFHHmHh1TuBh490bDlNy7KzPERljwskSVAw42ODyXvTXnuqk9+7G0Al9AaefwfUf7fM3IGNMWFmCigGxdv/Jy9tYYt0HNk6UMZ2JJagY0KCJeRQ/oNuYsZcNIiklEYCjBaca1BaNMbHNElSUU1UOxVgTc6/klETGXTYoML/OGksY02lYgopyJw6XUl7mDLHRtXsyPXp19TmijjfJc5lvw8f7qK7udAMlG9MpWYKKcof2NLz/FI1DbLRk6Pg+9OjdDYCzp8vZueaQzxEZY8LBElSU817ei7X7T3Xi4uKYOD07MG+X+YzpHCxBRblYbsHn5W3Nt23lQcrOVPgYjTEmHEKaoERkpojsEJF8EbmvkfUiIo+46zeKyGTPuqdE5KiIbA4q85CIbHe3f1VE0t3l2SJyTkTWu6/HQ/neIoGqNmjVFu2dxDanz6AeZOX2AqCmuoaNn9pw8MbEupAlKBGJB34PzALGALeKyJigzWYBue7rLuAxz7pngJmN7Po9YJyqTgB2Avd71u1W1Ynua2GHvJEIdvrEOc6eKgcguUsiGf1SfY4otLyNJayHc2NiXyhrUFOBfFXdo6qVwCJgTtA2c4Dn1LECSBeR/gCq+glw3kMvqvquqla7syuArJC9gwjnHaCwX3Z6TDaQ8LroyiHEJ8QDzvD2Nhy8MbEtlAlqIODtgrrQXdbWbZpzB/CWZz5HRNaJyMcickVjBUTkLhFZLSKri4uL23CoyHNkf0lguu+QdN/iCJduaSmMzBsQmF/3oTWWMCaWhTJBNfbvfHA/Na3ZpvGdi/wcqAaedxcVAYNVdRLwE+AvIpJ23s5Vn1DVPFXNy8zMbM2hIpY3QfXLTvctjnDy9nC+/kMbDt6YWBbKBFUIDPLMZwHBD7C0ZpvziMgCYDZwm7rfUKpaoarH3ek1wG5gRLujjwLeS3ydoQYFMDJvYGA4+JNHS9m7+ajPERljQiWUCWoVkCsiOSKSBMwDFgdtsxiY77bmuwQ4papFze1URGYCPwNuUtUyz/JMt2EGIjIUp+FFzN5Jr6mp5WjBqcB8Z6lBJSTGM+GKIYF5ayxhTOwKWYJyGzLcC7wDbANeUtUtIrJQROpa2C3BSSL5wJPA3XXlReQFYDkwUkQKReROd9WjQHfgvaDm5FcCG0VkA/BXYKGqxmzPoseLzlBd5XT5k9arK11Tk32OKHy8rfk2fXaAyorqZrY2xkSrhFDuXFWX4CQh77LHPdMK3NNE2VubWD68ieWvAK+0O9goE9yCrzMZNKIXvQemcezgaWc4+BUFTLwqp+WCxpioYj1JRKkGDSQ6yf2nOs5w8J5nopbaZT5jYpElqCjVmWtQAJOurq8x5a8/zKnjZc1sbYyJRpagolRnewYqWHpmN4ZN6Ac4XT5t+HifvwEZYzqcJagoVFFexYnDpYDT03dm1nmPe3UKkzzPRK1dasPBGxNrLEFFoaMHTgW+jHsN6E5iUkjbukSscZcNJjHZee9HDpRwaM/JFkoYY6KJJago1Bl7kGhMcpfg4eCtsYQxscQSVBQ6vL/+Ad3OeP/Jy/tM1HobDt6YmGIJKgod3ld/Kasz16AAhk3oS1qvrgCcPVVO/rrDPkdkjOkolqCiUIMm5p28BhUXF9egybl1fWRM7LAEFWVKT5UHBilMSon9QQpbw5ugtq4opKzUhoM3JhZYgooyDXowH9wj5gcpbI2+g9MZOLx+OPhNnx3wOSJjTEewBBVljh4oCUz3GdzDv0AijHecKGvNZ0xssAQVZYoL64c572sJKmDCldnExTsf5/3bijledMbniIwxF8oSVJTxjgGVmWUJqk5qj4bDwVtjCWOinyWoKHO0oL4GlTmoc3Zx1BRvD+frPrDh4I2JdpagokhZaQWlJecASExKoGefbj5HFFlGTRlIl1TPcPBbbDh4Y6KZJagoUuypPfUe2J24OPv1eSUkxnPRlfXDwa/7YK+P0RhjLpR9w0WR4sL6+099Btn9p8Y0HA5+P1WVNhy8MdHKElQUaXj/yRJUYwaN6EXvAc69uYpzVWxdUehzRMaY9rIEFUUatuCzBhKNEZGG40RZaz5jolZIE5SIzBSRHSKSLyL3NbJeROQRd/1GEZnsWfeUiBwVkc1BZTJE5D0R2eX+7OlZd7+7rx0i8uVQvjc/eBOUXeJrmrfro11rizh9woaDNyYahSxBiUg88HtgFjAGuFVExgRtNgvIdV93AY951j0DzGxk1/cBS1U1F1jqzuPuex4w1i33X24MMaGqspqSo2cBp5bQe0B3nyOKXD37pDJ0fF/AhoM3JpqFsgY1FchX1T2qWgksAuYEbTMHeE4dK4B0EekPoKqfACca2e8c4Fl3+lngq57li1S1QlX3AvluDDGhuPB0/Si6/buTkBgzuTckvI0l1tozUcZEpVAmqIFAgWe+0F3W1m2C9VXVIgD3Z5+27EtE7hKR1SKyuri4uMU3ESm8XRzZA7otG3fZIBKTnOHgD+87ycH8xv7XMcZEslAmqMa62Q7+N7Y123Tk8VDVJ1Q1T1XzMjMz23mo8LP7T22T0jWJcZfXDwe/6t18H6MxxrRHswlKRJ7xTC9o474LgUGe+SzgUDu2CXak7jKg+7Ouu4D27CtqHD1gLfjaaurM3MD0+o/2UXGuysdojDFt1VIN6iLP9A/buO9VQK6I5IhIEk4DhsVB2ywG5rut+S4BTtVdvmvGYqAuWS4AXvcsnyciySKSg9PwYmUbY45YVoNquyGjMwPnqrK8ig2f7PM3IGNMm7SUoNp9Z1lVq4F7gXeAbcBLqrpFRBaKyEJ3syXAHpwGDU8Cd9eVF5EXgOXASBEpFJE73VUPAteJyC7gOnceVd0CvARsBd4G7lHVmvbGH0lqamo5XlQamLcaVOuICFO+PDwwv+odu8xnTDRJaGF9log8gnN/p246QFV/0FxhVV2Ck4S8yx73TCtwTxNlb21i+XHg2ibW/RL4ZXMxRaMTh0upqXZybVqvrqR0TfI5ougx+ZqhvPPseqqraijcdZxDe04wYGiG32EZY1qhpRrU/wDWAKs9096XCQO7vNd+XbsnM/YyT2MJq0UZEzWarUGp6rPNrTfhUWwJ6oJMnZkbeFh33Yd7mfmdSSSnJPoblDGmRc0mKBEJbtTQgKre1LHhmMZYH3wXJmdsH3oPTOPYwdNUnKti06f7ybtueMsFjTG+auke1KU4D7++AHxB488amRDzPqRrNai2ExGmfnk4S55aC8DKd/ItQRkTBVq6B9UP+F/AOOB3OK3mjqnqx6r6caiDM05fcjbM+4WbfO1Q4hOc7qEKdhyjaO9JnyMyxrSk2QSlqjWq+raqLgAuwWkO/pGI/FNYojOcOl5GZbnzgGmX1GRSe6T4HFF06paWwjhPY4mV1ljCmIjXYldH7oOvXwf+G6dJ+CPA30IdmHF4h3nvMygNEbvK2l7eZ6LWf7iXinLrWcKYSNZSV0fPAsuAycD/VdUpqvqvqnowLNEZa2LegYaO7xsYbbe8rJL1H+71OSJjTHNaqkHdDozA6eZouYicdl9nROR0C2VNB2jQgs8S1AURES75yojA/PI3d9owHMZEsJbuQcWpanfPK819dVdVu1sfBtaCr2NNvnYoSe4zUEcOlLBn0xGfIzLGNKWlS3wpIvIjEXnUHUeppWbppoM1vMRn/xNcqC7dkph8Tf2Q8Cv+vtPHaIwxzWnpEt+zQB6wCbgBeDjkEZmAsjMVnD1VDkBicgLpmd18jig2XDp7ZGB6y/ICSorP+hiNMaYpLSWoMar6LVX9A3AzcEUYYjKuBvefBloLvo7SZ1APhk3oBzjPmX3x1i6fIzLGNKalBBVoh+sOn2HCqGEDCbu815EuvbG+FrXqnXyqKu3jbUykaXHAQm/LPWCCteILn6MF1kAiVEZNGRi4ZHr2dDkbP93vc0TGmGAtteKLD2q5l2Ct+MKnuNCegQqV+Pg4pt1gTc6NiWQt9iRh/FNsNaiQmnL9MBISnf75DuYfp2DncZ8jMsZ4WYKKUBXlVZw86gzzHhcXR0b/VJ8jij3d0lK46KrswPyyxdv9C8YYcx5LUBHq2MEzgemM/qkkuD1xm47lbXK+6fMDlByzJufGRApLUBHK+uALj4HDMsgZ1xeA2ppali3e4XNExpg6IU1QIjJTRHaISL6I3NfIehGRR9z1G0VkcktlReRFEVnvvvaJyHp3ebaInPOsezyU7y3UbJj38Lnia6MD06veyae8rNLHaIwxdUKWoEQkHvg9MAsYA9wqImOCNpsF5Lqvu4DHWiqrqt9Q1YmqOhF4hYZDf+yuW6eqC0P13sLBujgKn1FTBtJ7YH0v56ve3e1zRMYYCG0NaiqQr6p7VLUSWATMCdpmDvCcOlYA6SLSvzVlxelW4Rac4ehjToNRdLOsBhVKIsKXvlpfi1q2eDvV1TU+RmSMgdAmqIFAgWe+0F3Wmm1aU/YK4IiqevupyRGRdSLysYg02i2T2+ntahFZXVxc3Pp3E0bV1TUcL6pvJGG9SITe5Gty6OaOVlxSfJYtywpaKGGMCbVQJqjGOo4LfhKyqW1aU/ZWGtaeioDBqjoJ+AnwFxE575tdVZ9Q1TxVzcvMzGwyeD+dKCqltqYWgB69u5HsDg9hQicxKYFLPWNFffrqNntw1xifhTJBFQKDPPNZwKFWbtNsWXfYj68DL9YtU9UKVT3uTq8BduMMthh17P6TP6bdMILEJGdEmYP5x8nfcNjniIzp3EKZoFYBuSKSIyJJwDxgcdA2i4H5bmu+S4BTqlrUirIzgO2qWli3QEQy3cYViMhQnIYXe0L15kLJBin0R2qPFC6eMTQw//HLW3yMxhgTsgTl9n5+L/AOsA14SVW3iMhCEalrYbcEJ4nkA08CdzdX1rP7eZzfOOJKYKOIbAD+CixU1RMheXMhZsO8++fKfxhDXJzzZ7F742H2b4/M+5TGdAYhHSFXVZfgJCHvssc90wrc09qynnXfbmTZKzjNzqNegwSVZZf4wqlnn1QmTs9m7QdO5fvjl7cw/5+n+xuUMZ2U9SQRYVTVLvH57KqbxwYGh9y2spCifSd9jsiYzskSVIQpKT5LVYUzeF7X7smkuk2fTfj0GdSDsZfWt9Gxe1HG+MMSVISxQQojw1U3jw1Mb/x0P8cO2ficxoSbJagI4x2k0BpI+Ccrtxe5k/oDzmXXD1/c7HNExnQ+lqAizNED9gxUpLj21gmB6XUf7qX4oNWijAknS1ARxtuCr+/gdP8CMQwZndmgFrX0hY0+R2RM52IJKoKoasNOYq0G5bsZt10UmN74yf4G/0AYY0LLElQEKS0p51xpBQDJXRLp0aurzxGZwSN7M+LiAYDVoowJN0tQESR4FN26Z3GMv67z1KI2fXaAw/tL/AvGmE7EElQEadjFkV3eixRZub0YPTULcGpR7/55vb8BGdNJWIKKIN4WfNZAIrLMuK2+Rd+2LwrZv8366DMm1CxBRZDgS3wmcgwYmsFFV2UH5t96eq2NF2VMiFmCiiDWgi+yXfeti4hPiAdg/7Zitq866HNExsQ2S1ARouxMBaUl5wBndNeefbr5HJEJ1qtfd6bNyg3Mv/3MOmpra32MyJjYZgkqQgQPsVE3JpGJLFd/YxzJXRIB53e2dmlUjolpTFSwb8EIYS34okNqjxSu+PqYwPx7/72RivIqHyMyJnZZgooQ3vtP1oIvsn3pq6NIy3Aeoj59osyG4zAmRCxBRYijB0oC09ZJbGRLTknk+vn1D+9++uo2Thwp9TEiY2KTJagI0bAFnzUxj3STrxlKVm4vAKqranjrqbU+R2RM7LEEFQHKyyo5dewsAPEJ8WT0S/U5ItMSEWH2XXmB+c3LDrBn0xEfIzIm9oQ0QYnITBHZISL5InJfI+tFRB5x128UkcktlRWRB0TkoIisd183eNbd726/Q0S+HMr31pGKC+trT70HdCfBfdbGRLYhozKZOD0nMP/GE6upqbFm58Z0lJAlKBGJB34PzALGALeKyJigzWYBue7rLuCxVpb9japOdF9L3DJjgHnAWGAm8F/ufiJegx4kBtvlvWjy5QUTSUxOAODwvpMsf3OHzxEZEztCWYOaCuSr6h5VrQQWAXOCtpkDPKeOFUC6iPRvZdlgc4BFqlqhqnuBfHc/Ec/bB19mljWQiCbpvbtxzbzxgfn3/nsjJe7lWmPMhQllghoIFHjmC91lrdmmpbL3upcEnxKRnm04HiJyl4isFpHVxcWR0eHnkQadxFoNKtp86aujAo8GVJZX8eYTq/0NyJgYEcoE1dhgRsG9aza1TXNlHwOGAROBIuDhNhwPVX1CVfNUNS8zM7ORIuF3eF9JYLpfds+mNzQRKSEhnq/eU19Z37K8gG0rC32MyJjYEMoEVQgM8sxnAYdauU2TZVX1iKrWqGot8CT1l/Fac7yIU1Za0bAFX39rwReNssf0Ie+64YH5xY+vsh4mjLlAoUxQq4BcEckRkSScBgyLg7ZZDMx3W/NdApxS1aLmyrr3qOp8Ddjs2dc8EUkWkRychhcrQ/XmOsqR/Q0v71kLvug16zuT6No9GYCS4rO88+x6fwMyJsqFLEGpajVwL/AOsA14SVW3iMhCEVnobrYE2IPToOFJ4O7myrpl/lNENonIRuBq4MdumS3AS8BW4G3gHlWtCdX76yiH950MTPcdYvefolnX7snM/seLA/PL39xhz0YZcwESQrlztwn4kqBlj3umFbintWXd5bc3c7xfAr9sb7x+8N5/6p9j95+i3cTpOWz8dH9grKi//m45P3z0KySnJPocmTHRx3qS8Jk3QfUdku5bHKZjiAhfu3caKd2SADh5pNQu9RnTTpagfKSqHNlfEpjvl53uWyym46RldOVGTzdIy9/cwe6Nh32MyJjoZAnKRyePnqXinNPSq2v3ZLr37OJzRKajTLo6h1FT6h/De+nXyyg7U+FjRMZEH0tQPvI2kOif0xORxh7lMtGo7lJfXau+08fLePXRL3BuuxpjWsMSlI/s/lNsS8voyj/88JLA/OZlB1j1br6PERkTXSxB+ahhC7503+IwoTNm2iAuuWFEYP7NJ9ZwtPBUMyWMMXUsQfnIalCdww13TqaPOwhlVWU1f3nwU+tlwphWsATlk6rKao4ddMaBEhEbZiOGJSYlMO9/fomERKeXkCP7S3jt9yvtfpQxLbAE5ZOjBacDX1AZ/VLtQc4Y1z+7J3O+PyUwv/6jvaxYstPHiIyJfJagfBLcgs/EvrzrhjPl+voOZf/+x7Uc2HHMx4iMiWyWoHxi9586pxu/l8fA4b0AqKmu4fl//8QGODSmCZagfNJwDKh03+Iw4ZWYlMBt919R/3zUiTL+/K8fW6MJYxphCcoHqkrR3vpLfJagOpeefVL55n1XEBfv/Pkd2nOClx5eZo0mjAliCcoHp46XcfZUOQDJXRLp1b+7zxGZcBs2oR9fvbt+FN6tKwp42zqVNaYBS1A+OLT7RGB6wLAM6+Kok5py/XCu+NqYwPwnr2xh2RvbfYzImMhiCcoHB/MbJijTec389kRGT8sKzL/xxGrWfbTXx4iMiRyWoHzgrUENHGZNzDuzuLg4vvHTyxk8KjOw7K+/Xc6ONYd8jMqYyGAJygeHdtc3kKhrcmw6r+SURBb8Yjp9B6cDUFtTy/P//gl7Nttw8aZzswQVZqdPlHH6RBkASSmJ9B5oDSSMMx7Yd/7lGnr2SQWcrrCeeeBDS1KmU7MEFWbe2lP/nHTi4uxXYBw9enXljn+9hrSMrgBUVViSMp1bSL8dRWSmiOwQkXwRua+R9SIij7jrN4rI5JbKishDIrLd3f5VEUl3l2eLyDkRWe++Hg/le2uvg/nHA9N2ec8E6z0gje/++4zzkpQNGW86o5AlKBGJB34PzALGALeKyJigzWYBue7rLuCxVpR9DxinqhOAncD9nv3tVtWJ7mthaN7ZhfH2vZaVay34zPkyB56fpJ7+Px+yedkBnyMzJrxCWYOaCuSr6h5VrQQWAXOCtpkDPKeOFUC6iPRvrqyqvquq1W75FUAWUUJVKdhRX4Pyttwyxis4SdVU1/CXBz9l5Tu7fI7MmPAJZYIaCBR45gvdZa3ZpjVlAe4A3vLM54jIOhH5WESuaCwoEblLRFaLyOri4uLWvZMOcuzQGc6VVgDQLS2FjH6pYT2+iS6ZA9NY+ND19B6QBjj/4Lz66BcsXbTJukUynUIoE1Rj3SME/1U1tU2LZUXk50A18Ly7qAgYrKqTgJ8AfxGRtPN2ovqEquapal5mZnhrMAe21yfEQSN7WQ8SpkU9+6Tyvf+8vsH9yvef38CL/+9zqiqrmylpTPQLZYIqBAZ55rOA4KcPm9qm2bIisgCYDdym7r+Sqlqhqsfd6TXAbmBEh7yTDmKX90x7pPZI4R//fQbDJ/YPLNvwyT6e/F/vc+bkOR8jMya0QpmgVgG5IpIjIknAPGBx0DaLgflua75LgFOqWtRcWRGZCfwMuElVy+p2JCKZbuMKRGQoTsOLPSF8f21W4GkgMWhkbx8jMdEmuUsiC/7PdKbNqv+fq2DHMR798Vvs3xbeS9XGhEvIEpTbkOFe4B1gG/CSqm4RkYUiUtfCbglOEskHngTubq6sW+ZRoDvwXlBz8iuBjSKyAfgrsFBV6/sU8llFeVVgiA0RsRZ8ps0SEuKZ8/0p3Pi9KYHLw6ePl/HEfe/x6atb7b6UiTnSmT/UeXl5unr16rAca9f6Ip7656WAM4Lujx6dHZbjmti0a10Rix76jLIzFYFlo6dlcfMPLw0MhmhMtBCRNaqaF7zcujEIk72b6nsDGDq+r4+RmFiQO6k///S7GxpcKt72RSG/u/fv1tGsiRmWoMJkjyUo08HSM7tx14PXcflNowLLTp8o45kHPuD1x1baMPIm6lmCCoOK8ioKdta34Mse28fHaEwsSUiIZ/Y/5rHgF1eTmt4lsHzFkp389u432b7qoI/RGXNhLEGFwYHtx6itqQWc+0+pPVJ8jsjEmlFTBvLDR7/C2Evrn84oKT7Ls//yIc8/+Amnjpc1U9qYyGQJKgzs/pMJh9QeKdx2/5Xc8pPL6JZW/0/Q5s8P8PD3FvPBi5vs4V4TVSxBhUH++vqeqC1BmVASESZdPZQfPzabydcODSyvqqjmvf/ewK8XvsG6j/ZSW1vrY5TGtI4lqBA7U3KOwl3O/ae4uDiGTrAEZUKvW1oKc390Gd/95Qz6ZfcMLC8pPstLD3/O7+79Oxs+2WfPTpmIZgkqxHauORT4EhgyJpOuqfaMigmfYRP68U+/m8VX75nW4LLf0YJTLHroM357z5ts/Gy/JSoTkRL8DiDWeVtRjcwb4GMkprOKi4tj2sxcJlwxhM9e28ayxTsoL6sEnET1wq8+5b0BaVx640gunjGU5JREnyM2xmE1qBCqrq5h19qiwPzoqVEzdJWJQV26JXHdbRfxP/44h6tvGUeSJxEdO3SaN/6wigcXvMrf/7SGE0dKfYzUGIfVoEJo7+ajVJxzHpbs2SeVzKzzRv8wJuy6dk/m+tsncvmcUXz66ja+eGsX5WedGlV5WSWfvbaNz1/fzvCJ/bh4xjDGXJJFYpJ9VZjws09dCG34eF9gevS0LBv/yUSUbmkpzFwwiau/MY61S/ew7I0dHDt4GnAGR9y1rohd64pI6ZrE+C8NZvK1QxkyOtM+xyZsLEGFSGVFNZs+OxCYnzg9279gjGlGckoil35lJJfcMIKda4tYtng7u9YVBRpOlJdVsurdfFa9m09aRldGT8ti7KWDyBnfh4SEeJ+jN7HMElSIbF1RQKXbF1rvgWlk5fZqoYQx/hIRRl48gJEXD6Ck+CzrPtzL2qV7OHbodGCb0yfK+OKtnXzx1k5SuiUxMm8AIyYPYPjEfqRldPUxehOLLEGFyLoP9gamJ18z1C6LmKiSntmNq28Zx/S5Yzmw/RhrP9jDlmUFnD1dHtim/GwlGz7eF7iU3XdwOsMn9Wf4Rf0YPKq3DfthLpglqBA4WnCKXevqW+/Z5T0TrUSEIaMzGTI6kznfn8K+rcVsXV7AluUFlBSfbbDtkQMlHDlQwuevbwOgz6AeTtkxmQwa0ZveA7sTF2cNh03rWYIKgQ9f2hy4fj9qykB69kn1OSJjLlxcXBxDx/Vl6Li+fOW7F1O09yQ71xaRv66IfVuLqamuabD90YJTHC04xap38wFITEqgX046/XN6MmBYBgOG9qTvkHSSku1ryDTOPhkd7HjRmQat967+xnj/gjEmRESEAUMzGDA0g+k3j6WivIp9W4rJX1/E/q3FHNx9ItCDf52qymoKdhyjYMexBvtJz+xG74HdyczqQe+B3ek9MI3MrDR69Opql8Y7OUtQHUhVeeOJ1YHa0/CJ/RnsGfHUmFiVnJIYaGABTivWwl3H2b+1mAPbizm0+ySnT5w/5IeqcvJoKSePlja4LA4QnxBPemZX0jO7kd6nGz37OD979O5GWkYXUtNT6JKaZEkshlmC6kArluxkx+r6ro1mfHOCj9EY45+k5ITA5cA6pafKKdpzkoO7T3B4r/Pz+KEzTfYDWFNdw/GiMxwvOtPkceIT4klNT6F7zxT3Z33i6pKaRJfuyfXT7ispOcGSWpQIaYISkZnA74B44I+q+mDQenHX3wCUAd9W1bXNlRWRDOBFIBvYB9yiqifddfcDdwI1wA9U9Z1Qvj+vNUt38/c/rg3MXz5nNENGZ4br8MZEvNQeKeRO6k/upP6BZdVVNZw4XEpx4SmKD56huPAUxw6e4djB05SdqWhxnzXVNZw6dpZTx862uG2d+IR4uqQmkdItkaSUhMAruUvz8wmJccQnxpOQGE9CYhwJifHEuz8Tk+qn69ZZErxwIUtQIhIP/B64DigEVonIYlXd6tlsFpDrvqYBjwHTWih7H7BUVR8Ukfvc+Z+JyBhgHjAWGAC8LyIjVLXhndsOUHGuiu2rDnKutJJzpZXsWlfE3s31gxL2y+7Jl+df1NGHNSbmJCTG02dQD/oM6nHeuoryKkqOnuXk0bOUHD1LSbHz89TxMkpPllNaUh7o9LYtaqprKC05R2nJuY54C02KT4gnPiGOhMQ4JE6Ij3d+xsWfPx0XL858nBDnznunA9vHCQjOT0Dcn3WtIyVOEHHu7YkA0vx8ffng/UmDbXGK4p2om6/b1xVfG9Ph5zCUNaipQL6q7gEQkUXAHMCboOYAz6lTx18hIuki0h+ndtRU2TnAdLf8s8BHwM/c5YtUtQLYKyL5bgzLO/qNnTtbyaKHPmt0Xd8h6cz/56us7zJjLlBySiJ9B6fTd3B6k9tUVlRTWlLOmZPnKD15jjMnyzl7qjzwz+O50krKzlRwrrSS8rNVlJ2poLqqw/9nbVRNdQ011TVUlre8bSyItgQ1ECjwzBfi1JJa2mZgC2X7qmoRgKoWiUgfz75WNLKvBkTkLuAugMGDB7fh7dTr0i3pvGUiwoQrh/C1e6aR3MWGKzAmHJKSE8jom0pG39Y/ylFVWU3ZmUoqzlVRea6ayvK6VxWV5dVUNFhW7WxXXk11VQ011bXUVNVSVVnjzFfVUl1VQ7X7s34+PEkw1oUyQTV2ATb4bmhT27SmbHuOh6o+ATwBkJeX165R2pJSEphwxRBSuiWR0jWRnn1TGXvpILr37NKe3RljwigxKYEevUJ7hUNVqamppbqyhtoapbZWqa2pPW+6pqYWrdUWlrvlamqpqVFwG5XU1jo/tbZ+XhVQ56dq8/ONla+bb6xs/XsjEEOox7kM5W+pEBjkmc8CDrVym6Rmyh4Rkf5u7ak/cLQNx+sQIsKt//OKUOzaGBMDRISEhHjrTPcChbLfkVVArojkiEgSTgOGxUHbLAbmi+MS4JR7+a65souBBe70AuB1z/J5IpIsIjk4DS9WhurNGWOMCa2Q1aBUtVpE7gXewWkq/pSqbhGRhe76x4ElOE3M83GamX+nubLurh8EXhKRO4EDwFy3zBYReQmnIUU1cE8oWvAZY4wJD2nqIbnOIC8vT1evXu13GMYY06mJyBpVzQtebl0LG2OMiUiWoIwxxkQkS1DGGGMikiUoY4wxEckSlDHGmIjUqVvxiUgxsP8CdtEbONbiVv6K9BgjPT6I/BgjPT6I/Bgtvgt3ITEOUdXzhn/o1AnqQonI6saaRkaSSI8x0uODyI8x0uODyI/R4rtwoYjRLvEZY4yJSJagjDHGRCRLUBfmCb8DaIVIjzHS44PIjzHS44PIj9Hiu3AdHqPdgzLGGBORrAZljDEmIlmCMsYYE5EsQbWTiMwUkR0iki8i9/kUwyAR+VBEtonIFhH5obv8ARE5KCLr3dcNnjL3uzHvEJEvhyHGfSKyyY1jtbssQ0TeE5Fd7s+ePsY30nOe1ovIaRH5kZ/nUESeEpGjIrLZs6zN50xELnbPfb6IPCIijY063ZExPiQi20Vko4i8KiLp7vJsETnnOZePhzrGJuJr8+/Uh3P4oie+fSKy3l3uxzls6vslfJ9FZ2hfe7XlhTNG1W5gKM7ovxuAMT7E0R+Y7E53B3YCY4AHgJ82sv0YN9ZkIMd9D/EhjnEf0Dto2X8C97nT9wG/8iu+Rn6vh4Ehfp5D4EpgMrD5Qs4ZzoCdlwICvAXMCnGM1wMJ7vSvPDFme7cL2k9IYmwivjb/TsN9DoPWPwz8wsdz2NT3S9g+i1aDap+pQL6q7lHVSmARMCfcQahqkaqudafPANuAgc0UmQMsUtUKVd2LM1Dk1NBH2mgcz7rTzwJf9Sz3M75rgd2q2lzvIiGPUVU/AU40ctxWnzMR6Q+kqepydb4hnvOUCUmMqvquqla7syuArOb2EcoYmziHTYmYc1jHrWHcArzQ3D5CfA6b+n4J22fRElT7DAQKPPOFNJ8YQk5EsoFJwBfuonvdSy1PeargfsStwLsiskZE7nKX9VXVInD+CIA+PsbnNY+GXwiRcg6h7edsoDsd7jjr3IHzn3KdHBFZJyIfi8gV7jI/YmzL79TPc3gFcERVd3mW+XYOg75fwvZZtATVPo1dP/Wtvb6IpAKvAD9S1dPAY8AwYCJQhHOpAPyJ+3JVnQzMAu4RkSub2da38yoiScBNwMvuokg6h81pKh4/z+XPgWrgeXdRETBYVScBPwH+IiJpPsTY1t+pn7/rW2n4z5Jv57CR75cmN20ilnbHaAmqfQqBQZ75LOCQH4GISCLOh+d5Vf0bgKoeUdUaVa0FnqT+ElTY41bVQ+7Po8CrbixH3Gp/3SWKo37F5zELWKuqR9x4I+Ycutp6zgppeIktLHGKyAJgNnCbezkH95LPcXd6Dc69iRHhjrEdv1O/zmEC8HXgxbplfp3Dxr5fCONn0RJU+6wCckUkx/3Pex6wONxBuNep/wRsU9Vfe5b392z2NaCuldBiYJ6IJItIDpCLc/MyVPF1E5HuddM4N9E3u3EscDdbALzuR3xBGvzHGinn0KNN58y99HJGRC5xPyfzPWVCQkRmAj8DblLVMs/yTBGJd6eHujHuCXeMbf2d+nEOXTOA7aoauCzmxzls6vuFcH4WO6K1R2d8ATfgtGrZDfzcpxi+hFNV3gisd183AH8GNrnLFwP9PWV+7sa8gw5skdREfENxWvVsALbUnSegF7AU2OX+zPAjPs8xuwLHgR6eZb6dQ5xEWQRU4fz3eWd7zhmQh/MlvBt4FLfnmBDGmI9zD6Lus/i4u+0/uL//DcBa4MZQx9hEfG3+nYb7HLrLnwEWBm3rxzls6vslbJ9F6+rIGGNMRLJLfMYYYyKSJShjjDERyRKUMcaYiGQJyhhjTESyBGWMMSYiWYIyYSUiKiIPe+Z/KiIPdNC+nxGRmztiXy0cZ67bw/OHoT5WI8d+QER+2obtb5I29LYvIukicrdnfrqIvNnWODuK9/2KyL+IyAy/YjHhZwnKhFsF8HUR6e13IF51D0G20p3A3ap6daji6SiqulhVH2xDkXTg7pY28oOq/kJV3/c7DhM+lqBMuFUDTwA/Dl4RXAMSkVL353S3g8yXRGSniDwoIreJyEpxxpgZ5tnNDBH51N1utls+Xpyxila5HYV+z7PfD0XkLzgPcAbHc6u7/80i8it32S9wHmB8XEQeCtp+uoh8Is5YSFtF5HERiXPXXS8iy0VkrYi87PZvhohcK04HoJvE6cA02V2+T0R+5b7HlSIyvJH4honI2+J0xPupiIxqZJtvi8ijnvP7iIgsE5E9TdQ2HwSGiTPmUN37SxWRv4oz1tPzbm8AdWP8fOwe/52gnhrqjn+jiHzhvsf3RaSvu7xBTdA9x9nu9M/FGU/ofWCkZ5vA56OpY4vIR57ztlPcTlXdz8D/c8/zRhH5p9a+B+Ojjn4K3l72au4FlAJpOONE9QB+CjzgrnsGuNm7rftzOlCCMz5NMnAQ+L/uuh8Cv/WUfxvnH69cnKfzU4C7gP/tbpMMrMYZr2Y6cBbIaSTOAcABIBNIAD4Avuqu+wjIa6TMdKAcpweNeOA94GagN/AJ0M3d7mfAL9zYCoAR7vLncDrkxD0/dT1vzAfedKcfwB3TCOcp/lx3ehrwQSMxfRt41HN+XnbPzxicIWOCt8+m4RhK04FTOP2nxQHLcRJ0IrAMyHS3+wbwVCP76wmBDgG+Czwc/D7c+c3usS/G+WehK87nJN/zfp9xz2eTx3Z/N3XHuAF4353+Pk6fcnXjVWW09j3Yy79XAsaEmaqeFpHngB8A51pZbJW6XfyLyG7gXXf5JsB7qe0ldToD3SUie4BROH0ATvDUGHrgJLBKnL7C9jZyvCnAR6pa7B7zeZwB5l5rIc6VqrrHLfMCzpd5OU5C+NytfCThfNGPBPaq6k637LPAPcBv3fkXPD9/4z2IWwO7DHhZ6gcnTW4hNoDX3POzta420wor1e0XTpwRXrNx/mEYB7znHj8ep9ueYFnAi27NJAlo7Fx7XQG8qm5ffiLSWB+XI1s4dl2npmvcWMHp3+5xdcerUtUTIjKule/B+MQSlPHLb3H6FHvas6wa97KzexkpybOuwjNd65mvpeHnOLjvrrru/v9JVd/xrhCR6Tg1qMa0d9jspo7/nqreGnT8iW3YV/B+44ASVW1pH8G857G179FbpgbnfAuwRVUvbaHs/wf8WlUXu+f7AXd54HftSvFMt9T/WkvHrou3Lta6MsH7be17MD6xe1DGF6p6AngJp8FBnX04l3jAGZ0zsR27nisice59qaE4nVa+A3xfnKEDEJER4vSu3pwvgKtEpLc4DShuBT5uxfGnitPLfRzOJaPPcEaXvbzuPpKIdBWREcB2INtzf+n2oGN8w/Nzufcg6ozLs1dE5rr7FBG5qBXxteQMzvDeLdkBZIrIpe7xE0VkbCPb9cC5JAv1PWCD87ue7JadjHPJFZxLoV8TkS7i9IR/4wUc2+tdYKE4Q1kgIhnt3I8JI0tQxk8P49yfqfMkTlJYiXNPpanaTXN24HzJv4XTI3Q58EdgK7BWRDYDf6CFqwfu5cT7gQ9xe5BW1dYMEbAcp6HBZpzLWa+6lwm/DbwgIhtxEtYoN7bv4Fym24RTG3zcs69kEfkC5z7beY1KgNuAO0Wkrrf4Oa2Ir1nqjDn0udto4aFmtqvEuR/0K/f463EuOQZ7AOf9fQoc8yx/BchwLxl+H2dkANQZYvxFd3+vAJ9ewLG9/ohzT3GjW+ab7dyPCSPrzdyYDuJewvqpqs7ugH3tw2mIcaylbY2JVVaDMsYYE5GsBmWMMSYiWQ3KGGNMRLIEZYwxJiJZgjLGGBORLEEZY4yJSJagjDHGRKT/H4U779bM9V8QAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# And here's the posterior distribution\n",
    "\n",
    "posterior.plot(color='C4', label='posterior')\n",
    "decorate(xlabel='Number of people in the audience',\n",
    "         ylabel='PMF')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.548960Z",
     "iopub.status.busy": "2021-04-16T19:35:21.548544Z",
     "iopub.status.idle": "2021-04-16T19:35:21.553168Z",
     "shell.execute_reply": "2021-04-16T19:35:21.552768Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "486.2255161687084"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# If we have to guess the audience size,\n",
    "# we might use the posterior mean\n",
    "\n",
    "posterior.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.557382Z",
     "iopub.status.busy": "2021-04-16T19:35:21.556868Z",
     "iopub.status.idle": "2021-04-16T19:35:21.560089Z",
     "shell.execute_reply": "2021-04-16T19:35:21.560441Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.011543092507699223"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# And we can use prob_gt to compute the probability\n",
    "# of exceeding the capacity of the room.\n",
    "\n",
    "# It's about 1%, which may or may not satisfy the fire marshal\n",
    "\n",
    "posterior.prob_gt(1200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** I often see [rabbits](https://en.wikipedia.org/wiki/Eastern_cottontail) in the garden behind my house, but it's not easy to tell them apart, so I don't really know how many there are.\n",
    "\n",
    "Suppose I deploy a motion-sensing [camera trap](https://en.wikipedia.org/wiki/Camera_trap) that takes a picture of the first rabbit it sees each day.  After three days, I compare the pictures and conclude that two of them are the same rabbit and the other is different.\n",
    "\n",
    "How many rabbits visit my garden?\n",
    "\n",
    "To answer this question, we have to think about the prior distribution and the likelihood of the data:\n",
    "\n",
    "* I have sometimes seen four rabbits at the same time, so I know there are at least that many.  I would be surprised if there were more than 10.  So, at least as a starting place, I think a uniform prior from 4 to 10 is reasonable.\n",
    "\n",
    "* To keep things simple, let's assume that all rabbits who visit my garden are equally likely to be caught by the camera trap in a given day.  Let's also assume it is guaranteed that the camera trap gets a picture every day."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.563963Z",
     "iopub.status.busy": "2021-04-16T19:35:21.563510Z",
     "iopub.status.idle": "2021-04-16T19:35:21.565271Z",
     "shell.execute_reply": "2021-04-16T19:35:21.565631Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "hypos = np.arange(4, 11)\n",
    "prior = Pmf(1, hypos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.568713Z",
     "iopub.status.busy": "2021-04-16T19:35:21.568291Z",
     "iopub.status.idle": "2021-04-16T19:35:21.570022Z",
     "shell.execute_reply": "2021-04-16T19:35:21.570340Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# The probability that the second rabbit is the same as the first is 1/N\n",
    "# The probability that the third rabbit is different is (N-1)/N\n",
    "\n",
    "N = hypos\n",
    "likelihood = (N-1) / N**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.594482Z",
     "iopub.status.busy": "2021-04-16T19:35:21.590147Z",
     "iopub.status.idle": "2021-04-16T19:35:21.735651Z",
     "shell.execute_reply": "2021-04-16T19:35:21.735296Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAenklEQVR4nO3dfZQddZ3n8ffHACpP8hQgJMEwTJSJrEQMAR8OiohLopKgokEFBpkJKBGYIzMb3dVhdM4OIvjAEckJEIVVQAQZsk4EOVl0l1EwDyIQmAxtCNAkJJEgICgQ8tk/qhrLy026b9KVW+l8Xufc01W/+v2qvkWO/bHqVv9KtomIiGiaV3S7gIiIiHYSUBER0UgJqIiIaKQEVERENFICKiIiGikBFRERjZSAim2CpPMkfbcBdYyRZEnbbWD7RuuUtETSO+uqb2PKuv9yA9v+WtLtW7qmGNoSUDEkSPp95bNe0h8q6x8b5GN9R9Lz5b7XSrpV0kGDeYwNsf0G2z8t6+g3dCUtr/y3WCXp25J23hK1RmyuBFQMCbZ37vsADwPvr7R9r4ZDXlAeayTwKHBFDccYLO8vaz0UOAz4H60dNnRFF9FNCajYluwg6SpJT5e3yib0bZC0n6QbJK2R9KCkswayQ9t/AK4Dxlf29V5Jv5L0lKRHJJ3XZugnJK2QtFLSZ1q2vUrS98s6F0s6pLLv5ZLeLelY4HPAR8qro18PoNZHgR8DB5f7sqQzJT0APFC2/a2knvLKcK6k/Vp2M1nSMkm/lfQVSW1/h0g6qLyyXCtpqaQPV7Z9R9K3JP24rP3fJe0r6euSnpD0H5Le1N/5xNCXgIptyXHAtcBuwFzgmwDlL9n/Dfya4oroaOAcSf+1vx1K2gk4EeipND8DnFwe573AJyVNbRl6FDAWeA8wU9K7K9umAD8A9gCuBv5V0vbVwbZvBv4n8P3yKvEQ+iFpNDAZ+FWleSpwODBO0ruAfwE+DIwAHqL471V1PDCB4mpsCvCJNsfZCbi1rH1viv8+35L0hkq3D1Ncye0FPAf8Alhcrl8PfLW/84mhLwEV25Lbbc+z/SLwv4C+X+qHAcNtf9H287aXAZcB0zayr3Ml/Q54Gng7cFLfBts/tX2P7fW27wauAd7RMv6fbD9j+x7g2xS/xPsssn297RcoflG/CjhiU0+aIuB+B9wO/Iwi2Pr8i+215ZXgx4A5thfbfg74LPAWSWMq/b9c9n8Y+HpL3X3eByy3/W3b62wvBm4APlTpc6PtRbb/CNwI/NH2VeW/zfeBXEFFAiq2KY9Vlp+luJW2HfBaYD9Jv+v7UNw+22cj+7rQ9m7AGOAPwOv7Nkg6XNJt5e3CJ4EzKK4Mqh6pLD8E7Ndum+31QG/L9k5Ntb2b7dfa/lQZRu3q2K+spe/Yvwcep7iqHEjdfV4LHN7y3/NjwL6VPqsqy39os54HOSIBFUHxS/fB8pd432cX25P7G1heSZwNfEPSq8vmqyluIY62/RpgFqCWoaMry/sDK9ptK28/jmrZ/tLh+6tvAKr7WEERLn3H3gnYk+IhkJfVxsvr7vMI8LOW/5472/7kINQb25AEVAT8EnhK0n+T9GpJwyQdLOmwgQy2fSvFL+rpZdMuwFrbf5Q0Efhom2Gfl7Rj+b3MqRS3tfq8WdIHyqu7cyi+o7mjzT5WAWM29KDCJrgaOFXSeEmvpLgVeKft5ZU+fy9p9/L7rLNb6u7zI+B1kk6StH35OUzSXw1SnbGNSEDFNq/83uP9FE/iPQj8FrgceE0Hu/kK8A/lL/ZPAV+U9DTwBYqn/Fr9jOLBivkUtwt/Utl2E/AR4AmK77Y+UH4f1eoH5c/HJS3uoNa2bM8HPk/xfdFK4EBe/j3cTcAi4C7g32jzeL3tpyke/phGEdyPAV8GXrm5Nca2RXlhYURENFGuoCIiopESUBER0UgJqIiIaKQEVERENNI2MUHkXnvt5TFjxnS7jIiIaGPRokW/tT28tX2bCKgxY8awcOHCbpcRERFtSHqoXXtu8UVERCMloCIiopESUBER0UgJqIiIaKQEVERENFICKiIiGikBFRERjZSAioiIRkpARUREI20TM0lsig9+/votcpwbvvShLXKciIitTa6gIiKikRJQERHRSLUGlKRjJS2V1CNpZpvtH5N0d/n5uaRD+hsraQ9Jt0p6oPy5e53nEBER3VFbQEkaBlwCTALGASdKGtfS7UHgHbbfCHwJmD2AsTOB+bbHAvPL9YiIGGLqvIKaCPTYXmb7eeBaYEq1g+2f236iXL0DGDWAsVOAK8vlK4Gp9Z1CRER0S50BNRJ4pLLeW7ZtyGnAjwcwdh/bKwHKn3sPSrUREdEodT5mrjZtbttROooioN7e6dgNHlyaDkwH2H///TsZGhERDVDnFVQvMLqyPgpY0dpJ0huBy4Epth8fwNhVkkaUY0cAq9sd3PZs2xNsTxg+/GVvEo6IiIarM6AWAGMlHSBpB2AaMLfaQdL+wA+Bk2z/5wDHzgVOKZdPAW6q8RwiIqJLarvFZ3udpBnALcAwYI7tJZLOKLfPAr4A7Al8SxLAuvKqp+3YctfnA9dJOg14GDihrnOIiIjuqXWqI9vzgHktbbMqy38D/M1Ax5btjwNHD26lERHRNJlJIiIiGikBFRERjZSAioiIRkpARUREIyWgIiKikRJQERHRSAmoiIhopARUREQ0UgIqIiIaKQEVERGNlICKiIhGSkBFREQjJaAiIqKRElAREdFICaiIiGikBFRERDRSAioiIhqp1oCSdKykpZJ6JM1ss/0gSb+Q9Jykcyvtr5d0V+XzlKRzym3nSXq0sm1ynecQERHdUdsr3yUNAy4BjgF6gQWS5tq+r9JtLXAWMLU61vZSYHxlP48CN1a6fM32hXXVPtR88PPXb5Hj3PClD22R40TEtqHOK6iJQI/tZbafB64FplQ72F5tewHwwkb2czTwG9sP1VdqREQ0TZ0BNRJ4pLLeW7Z1ahpwTUvbDEl3S5ojafd2gyRNl7RQ0sI1a9ZswmEjIqKb6gwotWlzRzuQdgCOA35Qab4UOJDiFuBK4KJ2Y23Ptj3B9oThw4d3ctiIiGiAOgOqFxhdWR8FrOhwH5OAxbZX9TXYXmX7RdvrgcsobiVGRMQQU2dALQDGSjqgvBKaBsztcB8n0nJ7T9KIyurxwL2bVWVERDRSbU/x2V4naQZwCzAMmGN7iaQzyu2zJO0LLAR2BdaXj5KPs/2UpB0pngA8vWXXF0gaT3G7cHmb7RERMQTUFlAAtucB81raZlWWH6O49ddu7LPAnm3aTxrkMiMiooEyk0RERDRSAioiIhopARUREY2UgIqIiEZKQEVERCMloCIiopESUBER0UgJqIiIaKQEVERENFICKiIiGikBFRERjZSAioiIRkpARUREIyWgIiKikRJQERHRSAmoiIhopFpfWCjpWOAbFG/Uvdz2+S3bDwK+DRwK/HfbF1a2LQeeBl4E1tmeULbvAXwfGEPxRt0P236izvOIZvng56/fIse54Usf2iLHiYj2aruCkjQMuASYBIwDTpQ0rqXbWuAs4ELaO8r2+L5wKs0E5tseC8wv1yMiYoip8xbfRKDH9jLbzwPXAlOqHWyvtr0AeKGD/U4BriyXrwSmDkKtERHRMHUG1Ejgkcp6b9k2UAZ+ImmRpOmV9n1srwQof+7dbrCk6ZIWSlq4Zs2aDkuPiIhuqzOg1KbNHYx/m+1DKW4RninpyE4Obnu27Qm2JwwfPryToRER0QB1BlQvMLqyPgpYMdDBtleUP1cDN1LcMgRYJWkEQPlz9aBUGxERjVJnQC0Axko6QNIOwDRg7kAGStpJ0i59y8B7gHvLzXOBU8rlU4CbBrXqiIhohNoeM7e9TtIM4BaKx8zn2F4i6Yxy+yxJ+wILgV2B9ZLOoXjiby/gRkl9NV5t++Zy1+cD10k6DXgYOKGuc4iIiO6p9e+gbM8D5rW0zaosP0Zx66/VU8AhG9jn48DRg1hmREQ0UGaSiIiIRkpARUREIyWgIiKikRJQERHRSAmoiIhopARUREQ0UgIqIiIaKQEVERGNlICKiIhGqnUmiYjYuLwdOGLDcgUVERGNlICKiIhGSkBFREQjJaAiIqKRElAREdFICaiIiGikWgNK0rGSlkrqkTSzzfaDJP1C0nOSzq20j5Z0m6T7JS2RdHZl23mSHpV0V/mZXOc5REREd9T2d1CShgGXAMcAvcACSXNt31fpthY4C5jaMnwd8BnbiyXtAiySdGtl7NdsX1hX7RER0X11XkFNBHpsL7P9PHAtMKXawfZq2wuAF1raV9peXC4/DdwPjKyx1oiIaJg6A2ok8EhlvZdNCBlJY4A3AXdWmmdIulvSHEm7b2DcdEkLJS1cs2ZNp4eNiIguqzOg1KbNHe1A2hm4ATjH9lNl86XAgcB4YCVwUbuxtmfbnmB7wvDhwzs5bERENECdAdULjK6sjwJWDHSwpO0pwul7tn/Y1257le0Xba8HLqO4lRgREUNMnQG1ABgr6QBJOwDTgLkDGShJwBXA/ba/2rJtRGX1eODeQao3IiIapLan+GyvkzQDuAUYBsyxvUTSGeX2WZL2BRYCuwLrJZ0DjAPeCJwE3CPprnKXn7M9D7hA0niK24XLgdPrOoeIiOieWl+3UQbKvJa2WZXlxyhu/bW6nfbfYWH7pMGsMSIGT14fEoMpM0lEREQjJaAiIqKRElAREdFICaiIiGikjQaUpO9Ulk+pvZqIiIhSf1dQh1SWz95gr4iIiEHWX0B1NDVRRETEYOnv76BGSbqY4m+S+pZfYvus2iqLiIhtWn8B9feV5YV1FhIREVG10YCyfeWWKiQiIqJqowElaaOTu9o+bnDLiYiIKPR3i+8tFC8dvIbihYFt58eLiBhqMq9g9/UXUPsCxwAnAh8F/g24xvaSuguLiIht20YfMy9fDHiz7VOAI4Ae4KeSPr1FqouIiG1Wv6/bkPRK4L0UV1FjgIuBH25sTERExObq7yGJK4GDgR8D/2Q7b6+NiIgtor8rqJOAZ4DXAWdL6ptZQoBt71pncRERse3q7zuoV9jepfLZtfzsMpBwknSspKWSeiTNbLP9IEm/kPScpHMHMlbSHpJulfRA+XP3Tk44IiK2Dv3NZv4qSedI+qak6ZIG/Ip4ScOAS4BJwDjgREnjWrqtBc4CLuxg7Exgvu2xwPxyPSIihpj+Jou9EpgA3ANMBi7qYN8TgR7by2w/D1wLTKl2sL3a9gLghQ7GTinr6qtvagc1RUTEVqK/K6Jxtv8LgKQrgF92sO+RFH/k26cXOHwQxu5jeyWA7ZWS9m63A0nTgekA+++/fwdlR0REE/QXUC9d2dheJ3U0kUS7zgN9fcfmjC0627OB2QATJkzIa0MiYpu2Nc6M0V9AHSLpqXJZwKvL9YE8xdcLjK6sjwJWDLCujY1dJWlEefU0Alg9wH1GRMRWpL+n+Ia1PLm3XQdP8S0Axko6QNIOwDRgo5PPDnDsXKDv9fOnADcNcJ8REbEVGfBTeZ0qbwnOAG4BhgFzbC+RdEa5fZakfSneM7UrsF7SORTfez3Vbmy56/OB6ySdBjwMnFDXOURERPfUFlAAtucB81raZlWWH6O4fTegsWX748DRg1tpREQ0TX+PmUdERHRFAioiIhopARUREY2UgIqIiEZKQEVERCMloCIiopESUBER0UgJqIiIaKQEVERENFICKiIiGikBFRERjZSAioiIRkpARUREIyWgIiKikRJQERHRSAmoiIhopFoDStKxkpZK6pE0s812Sbq43H63pEPL9tdLuqvyeap82y6SzpP0aGXb5DrPISIiuqO2N+pKGgZcAhwD9AILJM21fV+l2yRgbPk5HLgUONz2UmB8ZT+PAjdWxn3N9oV11R4REd1X5xXURKDH9jLbzwPXAlNa+kwBrnLhDmA3SSNa+hwN/Mb2QzXWGhERDVNnQI0EHqms95ZtnfaZBlzT0jajvCU4R9Lu7Q4uabqkhZIWrlmzpvPqIyKiq+oMKLVpcyd9JO0AHAf8oLL9UuBAiluAK4GL2h3c9mzbE2xPGD58eAdlR0REE9QZUL3A6Mr6KGBFh30mAYttr+prsL3K9ou21wOXUdxKjIiIIabOgFoAjJV0QHklNA2Y29JnLnBy+TTfEcCTtldWtp9Iy+29lu+ojgfuHfzSIyKi22p7is/2OkkzgFuAYcAc20sknVFunwXMAyYDPcCzwKl94yXtSPEE4Oktu75A0niKW4HL22yPiIghoLaAArA9jyKEqm2zKssGztzA2GeBPdu0nzTIZUZERANlJomIiGikBFRERDRSAioiIhopARUREY2UgIqIiEZKQEVERCMloCIiopESUBER0UgJqIiIaKQEVERENFICKiIiGikBFRERjZSAioiIRkpARUREIyWgIiKikRJQERHRSLUGlKRjJS2V1CNpZpvtknRxuf1uSYdWti2XdI+kuyQtrLTvIelWSQ+UP3ev8xwiIqI7agsoScOAS4BJwDjgREnjWrpNAsaWn+nApS3bj7I93vaESttMYL7tscD8cj0iIoaYOq+gJgI9tpfZfh64FpjS0mcKcJULdwC7SRrRz36nAFeWy1cCUwex5oiIaIg6A2ok8EhlvbdsG2gfAz+RtEjS9EqffWyvBCh/7t3u4JKmS1ooaeGaNWs24zQiIqIb6gwotWlzB33eZvtQituAZ0o6spOD255te4LtCcOHD+9kaERENECdAdULjK6sjwJWDLSP7b6fq4EbKW4ZAqzquw1Y/lw96JVHRETX1RlQC4Cxkg6QtAMwDZjb0mcucHL5NN8RwJO2V0raSdIuAJJ2At4D3FsZc0q5fApwU43nEBERXbJdXTu2vU7SDOAWYBgwx/YSSWeU22cB84DJQA/wLHBqOXwf4EZJfTVebfvmctv5wHWSTgMeBk6o6xwiIqJ7agsoANvzKEKo2jarsmzgzDbjlgGHbGCfjwNHD26lERHRNJlJIiIiGikBFRERjZSAioiIRkpARUREIyWgIiKikRJQERHRSAmoiIhopARUREQ0UgIqIiIaKQEVERGNlICKiIhGSkBFREQjJaAiIqKRElAREdFICaiIiGikBFRERDRSrQEl6VhJSyX1SJrZZrskXVxuv1vSoWX7aEm3Sbpf0hJJZ1fGnCfpUUl3lZ/JdZ5DRER0R21v1JU0DLgEOAboBRZImmv7vkq3ScDY8nM4cGn5cx3wGduLJe0CLJJ0a2Xs12xfWFftERHRfXVeQU0Eemwvs/08cC0wpaXPFOAqF+4AdpM0wvZK24sBbD8N3A+MrLHWiIhomDoDaiTwSGW9l5eHTL99JI0B3gTcWWmeUd4SnCNp90GrOCIiGqPOgFKbNnfSR9LOwA3AObafKpsvBQ4ExgMrgYvaHlyaLmmhpIVr1qzpsPSIiOi2OgOqFxhdWR8FrBhoH0nbU4TT92z/sK+D7VW2X7S9HriM4lbiy9iebXuC7QnDhw/f7JOJiIgtq86AWgCMlXSApB2AacDclj5zgZPLp/mOAJ60vVKSgCuA+21/tTpA0ojK6vHAvfWdQkREdEttT/HZXidpBnALMAyYY3uJpDPK7bOAecBkoAd4Fji1HP424CTgHkl3lW2fsz0PuEDSeIpbgcuB0+s6h4iI6J7aAgqgDJR5LW2zKssGzmwz7nbafz+F7ZMGucyIiGigzCQRERGNlICKiIhGSkBFREQjJaAiIqKRElAREdFICaiIiGikBFRERDRSAioiIhopARUREY2UgIqIiEZKQEVERCMloCIiopESUBER0UgJqIiIaKQEVERENFICKiIiGikBFRERjVRrQEk6VtJSST2SZrbZLkkXl9vvlnRof2Ml7SHpVkkPlD93r/McIiKiO2oLKEnDgEuAScA44ERJ41q6TQLGlp/pwKUDGDsTmG97LDC/XI+IiCGmziuoiUCP7WW2nweuBaa09JkCXOXCHcBukkb0M3YKcGW5fCUwtcZziIiILtmuxn2PBB6prPcChw+gz8h+xu5jeyWA7ZWS9m53cEnTKa7KAH4vaemmnESH9gJ+28kA/XNNlWy+oXQuMLTOZyidCwyt8xlK5wJb7nxe266xzoBSmzYPsM9Axm6U7dnA7E7GbC5JC21P2JLHrMtQOhcYWuczlM4Fhtb5DKVzge6fT523+HqB0ZX1UcCKAfbZ2NhV5W1Ayp+rB7HmiIhoiDoDagEwVtIBknYApgFzW/rMBU4un+Y7AniyvH23sbFzgVPK5VOAm2o8h4iI6JLabvHZXidpBnALMAyYY3uJpDPK7bOAecBkoAd4Fjh1Y2PLXZ8PXCfpNOBh4IS6zmETbNFbijUbSucCQ+t8htK5wNA6n6F0LtDl85Hd0Vc7ERERW0RmkoiIiEZKQEVERCMloAaJpGGSfiXpR92uZXNJWi7pHkl3SVrY7Xo2h6TdJF0v6T8k3S/pLd2uaVNJen35b9L3eUrSOd2ua1NJ+jtJSyTdK+kaSa/qdk2bQ9LZ5bks2Rr/XSTNkbRa0r2Vtq5OLZeAGjxnA/d3u4hBdJTt8UPgbzq+Adxs+yDgELbifyPbS8t/k/HAmykeLLqxu1VtGkkjgbOACbYPpngYalp3q9p0kg4G/pZiFpxDgPdJGtvdqjr2HeDYlrauTi2XgBoEkkYB7wUu73Yt8SeSdgWOBK4AsP287d91tajBczTwG9sPdbuQzbAd8GpJ2wE78vK/k9ya/BVwh+1nba8DfgYc3+WaOmL7/wJrW5q7OrVcAmpwfB34B2B9l+sYLAZ+ImlROWXU1uovgDXAt8vbr5dL2qnbRQ2SacA13S5iU9l+FLiQ4k9FVlL8DeRPulvVZrkXOFLSnpJ2pPjzmdH9jNka/NnUckDbqeXqkoDaTJLeB6y2vajbtQyit9k+lGI2+TMlHdntgjbRdsChwKW23wQ8wxCY/b784/XjgB90u5ZNVX6XMQU4ANgP2EnSx7tb1aazfT/wZeBW4Gbg18C6rhY1BCSgNt/bgOMkLaeYdf1dkr7b3ZI2j+0V5c/VFN9xTOxuRZusF+i1fWe5fj1FYG3tJgGLba/qdiGb4d3Ag7bX2H4B+CHw1i7XtFlsX2H7UNtHUtwqe6DbNQ2Crk4tl4DaTLY/a3uU7TEUt13+j+2t9v8JStpJ0i59y8B7KG5fbHVsPwY8Iun1ZdPRwH1dLGmwnMhWfHuv9DBwhKQdJYni32arfYAFoO/NCpL2Bz7A1v9vBF2eWq7O2cxj67QPcGPxO4PtgKtt39zdkjbLp4HvlbfFllFOp7W1Kr/fOAY4vdu1bA7bd0q6HlhMcSvsV2z90wTdIGlP4AXgTNtPdLugTki6BngnsJekXuAf6fLUcpnqKCIiGim3+CIiopESUBER0UgJqIiIaKQEVERENFICKiIiGikBFds8SZZ0UWX9XEnnDdK+vyPpQ4Oxr36Oc0I5W/ttm7GP32+g/aeSXjZpsKTjJM0sl6dKGrepx45oJwEVAc8BH5C0V7cLqZI0rIPupwGfsn3UIO2vX7bn2j6/XJ0KJKBiUCWgIoo/FJ0N/F3rhtYroL6rDEnvlPQzSddJ+k9J50v6mKRflu/SOrCym3dL+n9lv/eV44dJ+oqkBZLulnR6Zb+3SboauKdNPSeW+79X0pfLti8AbwdmSfpKS/+X7U/Sv5YTAS9pnQxY0kWSFkuaL2l4ZdPHJf28PO7Esu9fS/qmpLdSzA34lfI9VQdKOkvSfeW5XTuwf4aIP5eZJCIKlwB3S7qggzGHULxmYS3FLBWX254o6WyKGSzOKfuNAd4BHAjcJukvgZMpZvA+TNIrgX+X1Deb90TgYNsPVg8maT+KCUnfDDxBMeP8VNtflPQu4Fzb7V4w2bq/T9heK+nVwAJJN9h+HNiJYo6/z5Sh94/AjHLMTrbfWk4cPAc4uG/ntn8uaS7wI9vXl7XOBA6w/Zyk3Tr4bxrxklxBRQC2nwKuoniJ3kAtsL3S9nPAb4C+gLmHIpT6XGd7ve0HKILsIIo5Dk+WdBdwJ7An0PeCu1+2hlPpMOCn5QSr64DvUbzvqj+t+ztL0q+BOyheCdF33PXA98vl71JclfW5Bl56Z9CuAwiduymmmPo4mdU7NlECKuJPvk7xXU71nVHrKP93Uk5qukNl23OV5fWV9fX8+d2J1vnEDAj4dN8bcm0fUHkf0jMbqE8DPI9WL+1P0jspZhJ/i+1DKObA29Cr1r2B5Xbrrd5LcVX6ZmBR+VLCiI4koCJKttcC11GEVJ/lFL9koXh/0fabsOsTJL2i/F7qL4ClwC3AJyVtDyDpdQN4meKdwDsk7VU+8HAixZtbO/Ea4Anbz0o6CDiisu0VQN/3bR8Fbq9s+0hZ59spbk0+2bLfp4G+WfBfAYy2fRvFizx3A3busM6IfAcV0eIi/vS9C8BlwE2SfgnMZ8NXNxuzlCJI9gHOsP1HSZdT3AZcXF6ZraGf12nbXinps8BtFFdT82x3+vqDm4EzJN1d1nVHZdszwBskLQKepAyl0hOSfg7sCnyizX6vBS6TdBbFa2eukPSass6v2f5dh3VGZDbziIhoptzii4iIRkpARUREIyWgIiKikRJQERHRSAmoiIhopARUREQ0UgIqIiIa6f8DDwIC3a/1LQoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "posterior = prior * likelihood\n",
    "posterior.normalize()\n",
    "\n",
    "posterior.bar(alpha=0.7)\n",
    "decorate(xlabel='Number of rabbits',\n",
    "         ylabel='PMF',\n",
    "         title='The Rabbit Problem')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Suppose that in the criminal justice system, all prison sentences are either 1, 2, or 3 years, with an equal number of each.  One day, you visit a prison and choose a prisoner at random.  What is the probability that they are serving a 3-year sentence?  What is the average remaining sentence of the prisoners you observe?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.741573Z",
     "iopub.status.busy": "2021-04-16T19:35:21.741016Z",
     "iopub.status.idle": "2021-04-16T19:35:21.743433Z",
     "shell.execute_reply": "2021-04-16T19:35:21.743783Z"
    }
   },
   "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>probs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.333333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.333333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.333333</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "1    0.333333\n",
       "2    0.333333\n",
       "3    0.333333\n",
       "Name: , dtype: float64"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# Here's the prior distribution of sentences\n",
    "\n",
    "hypos = np.arange(1, 4)\n",
    "prior = Pmf(1/3, hypos)\n",
    "prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.750272Z",
     "iopub.status.busy": "2021-04-16T19:35:21.749630Z",
     "iopub.status.idle": "2021-04-16T19:35:21.752134Z",
     "shell.execute_reply": "2021-04-16T19:35:21.752482Z"
    }
   },
   "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>probs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.166667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.333333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.500000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "1    0.166667\n",
       "2    0.333333\n",
       "3    0.500000\n",
       "Name: , dtype: float64"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# If you visit a prison at a random point in time,\n",
    "# the probability of observing any given prisoner\n",
    "# is proportional to the duration of their sentence.\n",
    "\n",
    "likelihood = hypos\n",
    "posterior = prior * likelihood\n",
    "posterior.normalize()\n",
    "posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.756049Z",
     "iopub.status.busy": "2021-04-16T19:35:21.755543Z",
     "iopub.status.idle": "2021-04-16T19:35:21.758197Z",
     "shell.execute_reply": "2021-04-16T19:35:21.757838Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1666666666666665"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# The mean of the posterior is the average sentence.\n",
    "# We can divide by 2 to get the average remaining sentence.\n",
    "\n",
    "posterior.mean() / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** If I chose a random adult in the U.S., what is the probability that they have a sibling?  To be precise, what is the probability that their mother has had at least one other child.\n",
    "\n",
    "[This article from the Pew Research Center](https://www.pewsocialtrends.org/2015/05/07/family-size-among-mothers/) provides some relevant data.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "source": [
    "From it, I extracted the following distribution of family size for mothers in the U.S. who were 40-44 years old in 2014:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.774819Z",
     "iopub.status.busy": "2021-04-16T19:35:21.773579Z",
     "iopub.status.idle": "2021-04-16T19:35:21.902764Z",
     "shell.execute_reply": "2021-04-16T19:35:21.903248Z"
    },
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYeUlEQVR4nO3de7RkZX3m8e9jQ4IKCMiRtIC2GlDRDM2kRSdogggOIhFd3mCEoMHB5USDE4zBCwnEmCGJoGG8jC0QWhASEQ3IRLFF28sMAzbI1VZxISrS0g2GaxxHmt/8Uftg0ZzTp8/ps+u8p/v7WatW7dr7rff91W6op95du/ZJVSFJUmseNdcFSJI0EQNKktQkA0qS1CQDSpLUJANKktQkA0qS1CQDSrMqyf9IcuIs9fWkJPclWdA9XpHkjbPRd9ff55McPVv9TWPcv0pyR5KfTrL9zUlu717742d57Idec5LXJ/nGLPc/a//+UvwdlDZWkluAXYAHgHXAt4FPAEur6sEZ9PXGqvrSNJ6zAji3qs6Yzljdc08CfrOqjpzuc2dTkt2B7wFPrqo1E2zfGrgHeF5VXdtzLa9n8G/w/D7HkWbKGZSm6/erajvgycApwJ8BZ872IEm2mu0+G/Fk4M6JwqmzC7ANcOPoSpLaZEBpRqrq7qq6GHgtcHSSZwMkOTvJX3XLOye5JMldSX6W5OtJHpXkHOBJwOe6w1jvSLIoSSU5JsmPgC8PrRsOq6cluTLJ3UkuSrJTN9b+SW4drjHJLUkOTHIw8C7gtd1413bbHzpk2NX1niQ/TLImySeSPK7bNl7H0Ul+1B2ee/dk+ybJ47rnr+36e0/X/4HAcuCJXR1nr/e8PYHvdg/vSvLlbv3fJ/lxknuSXJXkBUPPOSnJBUnOTXJvkuuT7Jnknd3r+HGSFw+1n/AwaZIPJzl1vXWfS/K2CdomyQe6/u9Oct0k//7j/77jtwe7WRtJnpFkefffxXeTvGay/aktlwGlTVJVVwK3Ai+YYPPx3bYxBjODdw2eUkcBP2IwG9u2qv526Dm/BzwT+I+TDPkHwB8CT2RwqPH0jajxC8BfA//Ujbf3BM1e391eCDwV2Bb40Hptng88HXgR8OdJnjnJkP8deFzXz+91Nb+hO5z5EuC2ro7Xr1fn94BndQ93qKoDuuVvAouBnYDzgAuSbDP01N8HzgF2BL4FXMrg/+1dgb8EPjZJncOWAUckeRQMPlx0r/P8Cdq+GPhdYE9gBwYfUu5cv1FVjf/7bgu8CvgpcFmSxzII6vOAJwBHAB9J8qz1+9CWzYDSbLiNwZvn+n4JLGTwfcsvq+rrNfWXnidV1f1V9fNJtp9TVTdU1f3AicBr0p1EsYleB5xWVTdX1X3AO4HD15u9nVxVP+++G7oWeETQdbW8FnhnVd1bVbcApwJHzbSwqjq3qu6sqgeq6lTg1xkE5bivV9WlVfUAcAGDDwSnVNUvgX8EFiXZYYoxrgTuZhBKAIcDK6rq9gma/xLYDngGg++xV1XV6sn67maGnwBeW1U/Bg4Fbqmqf+he09XAhQxCTHqIAaXZsCvwswnW/x3wfeCLSW5OcsJG9PXjaWz/IbA1sPNGVblhT+z6G+57KwYzv3HDZ939G4NZ1vp2Bn5tgr52nWlhSY5Psqo7nHYXg9nZ8GseDpGfA3dU1bqhx0xS6/qWAeMnkRzJYFb2CFX1ZQazyw8DtydZmmT7SWp/HHARcGJVfb1b/WTgud2h37u61/Q64Dc2okZtQQwobZIkz2Hw5vuI05W7GcTxVfVUBoeh/iTJ+Cf0yWZSU82wdh9afhKDT/N3APcDjxmqawGDmcTG9nsbgzfO4b4f4OFv/hvjjq6m9fv6yTT7AaD7vunPgNcAO1bVDgxmOplJf1M4Fzgsyd4MDrP+82QNq+r0qvptBock9wT+dP023eHC84CvVNXwYcYfA1+tqh2GbttW1Ztn8bVoM2BAaUaSbJ/kUAaHkM6tqusnaHNokt9MEganTq/rbjB443/qDIY+MsleSR7D4PuVT3ezhe8B2yR5aQanar+HwaGwcbczONQ12X/z5wP/NclTkmzLr76zemA6xXW1fAp4X5LtkjwZ+BMGb/4zsR2DoFwLbJXkz4EJZyubqqpuZfB91znAhZMdZk3ynCTP7fbz/cD/5Vf/rsPeBzwWOG699ZcAeyY5KsnW3e05G/hOT1soA0rT9bkk9zL4FPxu4DTgDZO03QP4EnAfcDnwkapa0W37b8B7ukM8b5/G+OcAZzM43LYN8McwOKsQ+C/AGQxmK/czOEFj3AXd/Z1Jrp6g37O6vr8G/IDBm+5bp1HXsLd249/MYGZ5Xtf/TFwKfJ5BAP+wq2uqw6CbYhnwW0xyeK+zPfBx4F+7mu4E3j9BuyOA5wH/OnQm3+uq6l4GJ1oczmDm+lPgb3j4BwrJH+pK+pUkv8tgtrdouj++lmabMyhJwENXsTgOOMNwUgsMKEl03//cxeBnAR+c02Kkjof4JElNcgYlSWrSvLgg584771yLFi2a6zIkabN21VVX3VFVY1O3HI15EVCLFi1i5cqVc12GJG3Wkvxw6laj4yE+SVKTDChJUpMMKElSkwwoSVKTDChJUpMMKElSkwwoSVKTDChJUpMMKElSkwwoSVKT5sWljjS3Xnnip+e6hEld+N5XzXUJknriDEqS1KTeAyrJgiTfSnJJ93inJMuT3NTd79h3DZKk+WcUM6jjgFVDj08ALquqPYDLuseSJD1MrwGVZDfgpcAZQ6sPA5Z1y8uAl/dZgyRpfup7BvVB4B3Ag0Prdqmq1QDd/RMmemKSY5OsTLJy7dq1PZcpSWpNbwGV5FBgTVVdNZPnV9XSqlpSVUvGxpr5A4+SpBHp8zTz/YCXJTkE2AbYPsm5wO1JFlbV6iQLgTU91iBJmqd6m0FV1TurareqWgQcDny5qo4ELgaO7podDVzUVw2SpPlrLn4HdQpwUJKbgIO6x5IkPcxIriRRVSuAFd3yncCLRjGuJGn+8koSkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJvUWUEm2SXJlkmuT3Jjk5G79SUl+kuSa7nZIXzVIkuavrXrs+xfAAVV1X5KtgW8k+Xy37QNV9f4ex5YkzXO9BVRVFXBf93Dr7lZ9jSdJ2rz0+h1UkgVJrgHWAMur6opu01uSXJfkrCQ7TvLcY5OsTLJy7dq1fZYpSWpQrwFVVeuqajGwG7BvkmcDHwWeBiwGVgOnTvLcpVW1pKqWjI2N9VmmJKlBIzmLr6ruAlYAB1fV7V1wPQh8HNh3FDVIkuaXPs/iG0uyQ7f8aOBA4DtJFg41ewVwQ181SJLmrz7P4lsILEuygEEQfqqqLklyTpLFDE6YuAV4U481SJLmqT7P4rsO2GeC9Uf1NaYkafPhlSQkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElN6i2gkmyT5Mok1ya5McnJ3fqdkixPclN3v2NfNUiS5q8+Z1C/AA6oqr2BxcDBSZ4HnABcVlV7AJd1jyVJepjeAqoG7usebt3dCjgMWNatXwa8vK8aJEnzV6/fQSVZkOQaYA2wvKquAHapqtUA3f0TJnnusUlWJlm5du3aPsuUJDWo14CqqnVVtRjYDdg3ybOn8dylVbWkqpaMjY31VqMkqU0jOYuvqu4CVgAHA7cnWQjQ3a8ZRQ2SpPmlz7P4xpLs0C0/GjgQ+A5wMXB01+xo4KK+apAkzV9b9dj3QmBZkgUMgvBTVXVJksuBTyU5BvgR8Ooea5AkzVO9BVRVXQfsM8H6O4EX9TWuJGnz4JUkJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTeotoJLsnuQrSVYluTHJcd36k5L8JMk13e2QvmqQJM1fW/XY9wPA8VV1dZLtgKuSLO+2faCq3t/j2JKkea63gKqq1cDqbvneJKuAXfsaT5K0eRnJd1BJFgH7AFd0q96S5LokZyXZcZLnHJtkZZKVa9euHUWZkqSG9B5QSbYFLgTeVlX3AB8FngYsZjDDOnWi51XV0qpaUlVLxsbG+i5TktSYXgMqydYMwumTVfUZgKq6varWVdWDwMeBffusQZI0P/V5Fl+AM4FVVXXa0PqFQ81eAdzQVw2SpPmrz7P49gOOAq5Pck237l3AEUkWAwXcArypxxokSfNUn2fxfQPIBJv+pa8xJUmbjz5nUNIW7ZUnfnquS5jQhe991VyXIG0UL3UkSWqSASVJapIBJUlqkgElSWqSASVJatIGAyrJ2UPLR/dejSRJnalmUHsPLR/XZyGSJA2bKqBqJFVIkrSeqX6ou1uS0xlcEWJ8+SFV9ce9VSZJ2qJNFVB/OrS8ss9CJEkatsGAqqploypkVFq9/Ax4CRpJGrbBgEpy8Ya2V9XLZrccSZIGpjrE9x+AHwPnM/hz7RNdnVySpFk3VUD9BnAQcATwn4D/CZxfVTf2XZgkacu2wdPMuz/N/oWqOhp4HvB9YEWSt46kOknSFmvKvweV5NeBlzKYRS0CTgc+029ZkqQt3VQnSSwDng18Hji5qm4YSVWSpC3eVFeSOArYk8Flji5Pck93uzfJPRt6YpLdk3wlyaokNyY5rlu/U5LlSW7q7necnZciSdqcTPUd1KOqaruh2/bdbbuq2n6Kvh8Ajq+qZzL4/uqPkuwFnABcVlV7AJd1jyVJepiprma+TZK3JflQkmOTTPmd1biqWl1VV3fL9wKrgF2Bw4DxHwAvA14+o8olSZu1qQ7xLQOWANcDhwCnzmSQJIuAfRj8lmqXqloNgxADnjDJc45NsjLJyrVr185kWEnSPDbVjGivqvotgCRnAldOd4Ak2wIXAm+rqnuSjfutb1UtBZYCLFmyxKuqS9IWZqoZ1C/HF6rqgel2nmRrBuH0yaoaPzX99iQLu+0LgTXT7VeStPmb8g8WDp+5B/y7aZzFF+BMYFVVnTa06WJg/K/zHg1cNNPiJUmbr6muZr5gE/rej8Fp6tcnuaZb9y7gFOBTSY4BfgS8ehPGkCRtpjb6rLzpqqpvMPnFZV/U17iSpM3DVIf4JEmaEwaUJKlJBpQkqUkGlCSpSQaUJKlJBpQkqUkGlCSpSQaUJKlJBpQkqUkGlCSpSQaUJKlJBpQkqUkGlCSpSQaUJKlJBpQkqUkGlCSpSQaUJKlJBpQkqUm9BVSSs5KsSXLD0LqTkvwkyTXd7ZC+xpckzW99zqDOBg6eYP0Hqmpxd/uXHseXJM1jvQVUVX0N+Flf/UuSNm9z8R3UW5Jc1x0C3HGyRkmOTbIyycq1a9eOsj5JUgNGHVAfBZ4GLAZWA6dO1rCqllbVkqpaMjY2NqLyJEmtGGlAVdXtVbWuqh4EPg7sO8rxJUnzx0gDKsnCoYevAG6YrK0kacu2VV8dJzkf2B/YOcmtwF8A+ydZDBRwC/CmvsaXJM1vvQVUVR0xweoz+xpPkrR56S2gJGm6Xnnip+e6hEld+N5XzXUJWxwvdSRJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqUm8BleSsJGuS3DC0bqcky5Pc1N3v2Nf4kqT5rc8Z1NnAweutOwG4rKr2AC7rHkuS9Ai9BVRVfQ342XqrDwOWdcvLgJf3Nb4kaX4b9XdQu1TVaoDu/gkjHl+SNE80e5JEkmOTrEyycu3atXNdjiRpxEYdULcnWQjQ3a+ZrGFVLa2qJVW1ZGxsbGQFSpLaMOqAuhg4uls+GrhoxONLkuaJPk8zPx+4HHh6kluTHAOcAhyU5CbgoO6xJEmPsFVfHVfVEZNselFfY0qSNh/NniQhSdqyGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJvX2J983JMktwL3AOuCBqloyF3VIkto1JwHVeWFV3TGH40uSGuYhPklSk+YqoAr4YpKrkhw7RzVIkho2V4f49quq25I8AVie5DtV9bXhBl1wHQvwpCc9aS5qlCTNoTmZQVXVbd39GuCzwL4TtFlaVUuqasnY2NioS5QkzbGRB1SSxybZbnwZeDFww6jrkCS1bS4O8e0CfDbJ+PjnVdUX5qAOSVLDRh5QVXUzsPeox5UkzS+eZi5JapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJatKcBFSSg5N8N8n3k5wwFzVIkto28oBKsgD4MPASYC/giCR7jboOSVLb5mIGtS/w/aq6uar+H/CPwGFzUIckqWGpqtEOmLwKOLiq3tg9Pgp4blW9Zb12xwLHdg+fDnx3pIWqTzsDd8x1EfOM+2z63GfT9/Sq2m6uixi31RyMmQnWPSIlq2opsLT/cjRqSVZW1ZK5rmM+cZ9Nn/ts+pKsnOsahs3FIb5bgd2HHu8G3DYHdUiSGjYXAfVNYI8kT0nya8DhwMVzUIckqWEjP8RXVQ8keQtwKbAAOKuqbhx1HZpTHrqdPvfZ9LnPpq+pfTbykyQkSdoYXklCktQkA0qS1CQDStOW5Kwka5LcMM3nLZrsOUn+MsmB3fKKJI84PTjJ65N8aGZV9yfJ7km+kmRVkhuTHDeN526W+2QqSbZJcmWSa7t9dvI0nrtF7bMkC5J8K8kls9zvZPtnSZLTu+VJ91WS+2aznonMxe+gNP+dDXwI+MRsdVhVfz5bfc2BB4Djq+rqJNsBVyVZXlXf3pRON2WfJNmqqh7YlPF79gvggKq6L8nWwDeSfL6q/s+mdLqZ7rPjgFXA9hNtTHJLVS2arcGqaiUwo99DJVlQVetmqxZnUJq2qvoa8LMNtUmyS5LPdp+Qr03yO92mBUk+3n1q/mKSR3ftz+6uMrJ+P29I8r0kXwX2m/UXMwuqanVVXd0t38vgzWTX9dv1vU+69qcl+QrwN0meluQLSa5K8vUkzxhqd3qS/53k5onG6FsNjH8C37q7PeKMrS19nyXZDXgpcMYm9vOOJNd3+/CUoU2v7may30vygq7t/hPN1jL4adDlSb6Z5L1D6/fP4AjCecD13Yzv77p21yV501C7FUk+neQ7ST6ZZKILNzzEgFJfTge+WlV7A/8eGP8pwR7Ah6vqWcBdwCsn6yDJQuBkBm8oBzG4uHDTkiwC9gGumGDzKPbJnsCBVXU8g1OG31pVvw28HfjIULuFwPOBQ4FTmAPdG9k1wBpgeVW5zx7pg8A7gAdn2kGSlwAvZ3BJub2Bvx3avFVV7Qu8DfiLKbr6e+CjVfUc4KfrbdsXeHdV7QUcA9zdtXsO8J+TPKVrt0831l7AU5niQ6cBpb4cAHwUoKrWVdXd3fofVNU13fJVwKIN9PFcYEVVre0uLPxPPdU6K5JsC1wIvK2q7pmgySj2yQVVta6r5XeAC7oQ+BiDN9hx/1xVD3aHIXfZ2Nc4m7p9sJjB1WT2TfLsCZptsfssyaHAmqq6aoJt705yTVfnE8eXk3x4gq4OBP6hqv4NoKqGj358prufah/CIEzO75bPWW/blVX1g275xcAfdLVdATyewQeK8Xa3VtWDwDVTjel3UBq1XwwtrwMePUX7efFDvQy+R7kQ+GRVfWaq9uuZzX1yf3f/KOCuLgCmGnODh1n6VlV3JVkBHAxs7Ik3W8I+2w94WZJDgG2A7ZOcW1VHVtX7gPfBQ99BLd5AP2Hy1z/+mtaxcXkwWT/3Dy2HwSz00ocVkezPI//dNjimMyj15TLgzfDQoZwJv+CdwhXA/kke3wXAq2ezwNnSHUc/E1hVVadtoOnI9kk3g/tBkleP15hk7xmM14skY0l26JYfzeBT/ncmaLrF7rOqemdV7dadAHE48OWqOnIGXX0R+MMkjwFIstMMS/pfXR0Ar9tAu0uBN3f7miR7JnnsTAY0oDRtSc4HLgeenuTWJMdM0Ow44IVJrmdw+OBZ0x2nqlYDJ3VjfQm4esZF92s/4CjggKFDLYdM0G7U++R1wDFJrmXw3U1Lf3dtIfCVJNcxuD7n8qqa6DRq99kmqqovMLje6crusNvbZ9jVccAfJfkm8LgNtDsD+DZwdQY/B/gYMzxa56WOJElNcgYlSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWrS/wfBzZCMr1CZuwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "qs = [1, 2, 3, 4]\n",
    "ps = [22, 41, 24, 14]\n",
    "prior = Pmf(ps, qs)\n",
    "prior.bar(alpha=0.7)\n",
    "\n",
    "plt.xticks(qs, ['1 child', '2 children', '3 children', '4+ children'])\n",
    "decorate(ylabel='PMF',\n",
    "         title='Distribution of family size')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "source": [
    "For simplicity, let's assume that all families in the 4+ category have exactly 4 children."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.911850Z",
     "iopub.status.busy": "2021-04-16T19:35:21.911271Z",
     "iopub.status.idle": "2021-04-16T19:35:21.916004Z",
     "shell.execute_reply": "2021-04-16T19:35:21.916486Z"
    }
   },
   "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>probs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.094828</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.353448</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.310345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.241379</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "1    0.094828\n",
       "2    0.353448\n",
       "3    0.310345\n",
       "4    0.241379\n",
       "Name: , dtype: float64"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# When you choose a person a random, you are more likely to get someone\n",
    "# from a bigger family; in fact, the chance of choosing someone from\n",
    "# any given family is proportional to the number of children\n",
    "\n",
    "likelihood = qs\n",
    "posterior = prior * likelihood\n",
    "posterior.normalize()\n",
    "posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.920242Z",
     "iopub.status.busy": "2021-04-16T19:35:21.919748Z",
     "iopub.status.idle": "2021-04-16T19:35:21.922344Z",
     "shell.execute_reply": "2021-04-16T19:35:21.921990Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9051724137931034"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# The probability that they have a sibling is the probability\n",
    "# that they do not come from a family of 1 \n",
    "\n",
    "1 - posterior[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.926250Z",
     "iopub.status.busy": "2021-04-16T19:35:21.925631Z",
     "iopub.status.idle": "2021-04-16T19:35:21.928050Z",
     "shell.execute_reply": "2021-04-16T19:35:21.928401Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9051724137931034"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# Or we could use prob_gt again\n",
    "\n",
    "posterior.prob_gt(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** The [Doomsday argument](https://en.wikipedia.org/wiki/Doomsday_argument) is \"a probabilistic argument that claims to predict the number of future members of the human species given an estimate of the total number of humans born so far.\"\n",
    "\n",
    "Suppose there are only two kinds of intelligent civilizations that can happen in the universe.  The \"short-lived\" kind go exinct after only 200 billion individuals are born.  The \"long-lived\" kind survive until 2,000 billion individuals are born.\n",
    "And suppose that the two kinds of civilization are equally likely.\n",
    "Which kind of civilization do you think we live in?  \n",
    "\n",
    "The Doomsday argument says we can use the total number of humans born so far as data.\n",
    "According to the [Population Reference Bureau](https://www.prb.org/howmanypeoplehaveeverlivedonearth/), the total number of people who have ever lived is about 108 billion.\n",
    "\n",
    "Since you were born quite recently, let's assume that you are, in fact, human being number 108 billion.\n",
    "If $N$ is the total number who will ever live and we consider you to be a randomly-chosen person, it is equally likely that you could have been person 1, or $N$, or any number in between.\n",
    "So what is the probability that you would be number 108 billion?\n",
    "\n",
    "Given this data and dubious prior, what is the probability that our civilization will be short-lived?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.931981Z",
     "iopub.status.busy": "2021-04-16T19:35:21.931530Z",
     "iopub.status.idle": "2021-04-16T19:35:21.933115Z",
     "shell.execute_reply": "2021-04-16T19:35:21.933463Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "hypos = [200, 2000]\n",
    "prior = Pmf(1, hypos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.940200Z",
     "iopub.status.busy": "2021-04-16T19:35:21.939749Z",
     "iopub.status.idle": "2021-04-16T19:35:21.942386Z",
     "shell.execute_reply": "2021-04-16T19:35:21.942831Z"
    }
   },
   "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>probs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>200</th>\n",
       "      <td>0.909091</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000</th>\n",
       "      <td>0.090909</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "200     0.909091\n",
       "2000    0.090909\n",
       "Name: , dtype: float64"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "likelihood = 1/prior.qs\n",
    "posterior = prior * likelihood\n",
    "posterior.normalize()\n",
    "posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-16T19:35:21.946847Z",
     "iopub.status.busy": "2021-04-16T19:35:21.946057Z",
     "iopub.status.idle": "2021-04-16T19:35:21.948588Z",
     "shell.execute_reply": "2021-04-16T19:35:21.948054Z"
    }
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# According to this analysis, the probability is about 91% that our \n",
    "# civilization will be short-lived.  \n",
    "# But this conclusion is based on a dubious prior.\n",
    "\n",
    "# And with so little data, the posterior depends strongly on the prior.  \n",
    "# To see that, run this analysis again with a different prior, \n",
    "# and see what the results look like.\n",
    "\n",
    "# What do you think of the Doomsday argument?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}