day0/IntroductiontoBayesianReasoning.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "38ef8a40",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "431bdde7",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A (brief) Introduction to Bayesian Reasoning and Statistics\n",
    "\n",
    "#### Bryan Scott, CIERA/Northwestern, \n",
    "#### DSFP Session 19: Orientation for Cohort 7 \n",
    "version 0.1 "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b570ae7f",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Part 1: Warm Up - Conditional Probability and the Law of Total Probability\n",
    "\n",
    "One of the main objects of Bayesian inference is the conditional probability distribution. This tells us what the chance of something happening given that we know something else happened are. For example,\n",
    "\n",
    "$$\n",
    "p(Good | Dog) = 1 \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "235467a4",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "A related concept is the law of total probability - this just formalizes the idea that the odds of something happening must involve a sum over all of the (conditional) ways it can happen. \n",
    "\n",
    "$$\n",
    "p(B) = \\Sigma_n P(B | A_n) P(A_n)\n",
    "$$\n",
    "\n",
    "for example,\n",
    "\n",
    "$$\n",
    "p(Pet) = p(Pet| Own) p(Own) + p(Pet| Friendly) p(Friendly)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e863d8b9",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "One more useful bit of probability theory is that conditional probabilities can be related to each other via a simple relationship:\n",
    "\n",
    "$$\n",
    "p(A| B) = \\frac{P(B|A)P(A)}{P(B)}\n",
    "$$\n",
    "\n",
    "which also holds if you swap $A \\rightarrow B$. This result is sometimes called Bayes' theorem, and the interpretation of this equation (which is simply a fact about conditional probabilities) gives rise to Bayesian statistics."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71d308dd",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Part 2: Prior Beliefs\n",
    "\n",
    "Let's consider the following problem. Suppose you're deciding whether to take the bus or walk to your destination. Unfortunately, you live in Chicago, so while there is a bus tracker app, it will tell you only a set of probable times you might expect the bus to arrive. Having previously solved this problem empirically, one learns not to rely on the bus tracker app's estimates. \n",
    "\n",
    "As such, you'll call the rate at which buses arrive $\\lambda$ and use the bus schedule to set a best guess $p_0(\\lambda)$ for the rate at which buses arrived. $\\lambda$ is a $\\textit{parameter}$ of some distribution we wish to infer. It's assumed distribution before we have observed any buses arriving, $p_0(\\lambda)$ is called the $\\textit{prior}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23aa62fc",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Part 3: Conditioning Beliefs on Evidence\n",
    "\n",
    "Now that we know how often we expect the bus to arrive, we need someway to update our knowledge of $\\lambda$ based on observing bus arrival times. This is called the $\\textbf{likelihood}$. It tells us how our observed data changes our beliefs about the value of some parameter. \n",
    "\n",
    "While choosing a distribution that captures our prior beliefs about the value of a parameter often requires care, there are more clear rules governing the proper updating functions. It turns out that the right updating function for this problem is the $\\textit{Poisson Distribution}$.\n",
    "\n",
    "For most problems, you will instead use a $\\textit{Gaussian distribution}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb144172",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Part 4: Bayes' Insight\n",
    "\n",
    "\n",
    "The deep insight is that Bayesian probabilities are statements about the probability of a parameter taking on a certain value. This differs from the classical interpretation, where probability is about realizations of the data, rather than distributions over parameters. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7ab2844",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For example, what we care about as bus riders is the rate at which buses arrive. If we believe the next bus will arrive soon, we will likely wait for it. If we believe the bus will arrive much later, perhaps we will walk. What we want to quantify is our beliefs about the value of the rate parameter of bus arrivals. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46eaab30",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "To make this quantitative, the probability distribution for the rate at which buses arrive $p(\\tau | \\lambda)$, conditioned on the time we've waited so far, is proportional to the product of our $\\textit{prior}$ belief about when a bus will arrive and the $\\textit{likelihood}$ - a function that tells us how our beliefs change as we observe (or don't observe) buses arriving. \n",
    "\n",
    "$$\n",
    "p(\\lambda | \\tau) \\propto p_0(\\lambda) p(\\tau | \\lambda )\n",
    "$$\n",
    "\n",
    "There is a missing normalization constant here, sometimes called the $\\textbf{evidence}$ or the $\\textbf{marginal likelihood}$ - we will return to this point in detail in the Bayesian inference session, but for now, you can calculate the evidence by asserting that $\\Sigma_{\\tau} p(\\lambda | \\tau) = 1$. For most problems we can neglect this constant since we only care about the maximum of the probability distribution. When it does matter, we often use the law of total probability to calculate this term. More on this in a future session."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8841241e",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "What this equation says is that my $\\textit{belief}$ about the frequency at which buses arrive depends on my prior belief about how often the bus should come and a factor that tells me how to update my belief based on how long I've been waiting at the bus stop."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "967c18e1",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Part 5: Putting this into practice\n",
    "\n",
    "Now, in the rest of this notebook, work out the posterior estimate for the rate at which buses arrive based on your wait time so far. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9b79e6d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import gamma"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d4fae89",
   "metadata": {},
   "source": [
    "The choice of prior is, in general, a subtle problem that we will return to. In the analysis of wait times, it is common to select a prior from the $\\Gamma$ distribution family described by the equation:\n",
    "\n",
    "$$\n",
    "p_0 (\\lambda) = \\frac{\\beta^a \\lambda^{a-1} e^{-\\beta \\lambda}}{\\Gamma(a)}\n",
    "$$\n",
    "\n",
    "Begin by discussing with your partner what each of these parameters means in the context of a wait time problem."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c17524ed",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Next, either use a built in function from scipy or write your own function for the $\\Gamma(a)$ distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "40281475",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fec51680d30>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAwYklEQVR4nO3deXycV3no8d8zqzTavduSZdmJHcdO4jhx4iwEQkJCQoC0BCiBQqGlKZ9CC7QUKLSXyy1w2dqGXHJJ07S0XLjkJiFAaAJpScgCzmLHNvG+W5a8SrbW0YxmO/ePM4pGmnekkTSa9fl+Pu9nNGdezZyZOM8cPe85zxFjDEoppUqfq9AdUEoplRsa0JVSqkxoQFdKqTKhAV0ppcqEBnSllCoTnkK98Lx580xbW1uhXl4ppUrSK6+80m2Mme/0WMECeltbG1u2bCnUyyulVEkSkfZMj2nKRSmlyoQGdKWUKhMa0JVSqkxMGtBF5F9F5IyI7MzwuIjIPSJyUEReFZHLct9NpZRSk8lmhP5vwC0TPH4rsDJ53AV8Z+bdUkopNVWTznIxxjwnIm0TnHI78D1jq3y9KCKNIrLYGHMyV50c7QucPQsdHdDZCa2tsG5drl9FKaVKUy6mLTYDHSn3O5NtaQFdRO7CjuJpbW2d0ots3gzf/z6Ew6NtGzdqQFdKqRG5uCgqDm2ONXmNMfcbYzYYYzbMn+84Lz6j2tqxwRzsKF0ppZSVi4DeCSxNud8CnMjB847R0pLedvIkxGK5fiWllCpNuQjojwEfSM52uQrom438eV0dNDaObUskbFBXSimVRQ5dRH4IXA/ME5FO4AuAF8AYcx/wBPAW4CAwBHxotjrb0gK9vWPbOjpg6VLH05VSqqJkM8vlzkkeN8BHc9ajCSxdCjvHzYbv6HA+VymlKk1JrRR1yqPrhVGllLJKKqA7pVY6Ouz8dKWUqnQlFdDnzwe/f2xbKATnzhWmP0opVUxKKqC7XNDcnN6uaRellCqxgA7OeXS9MKqUUiUY0J3y6MeO5b8fSilVbEouoDuVgNGArpRSJRjQm5ttLj1VTw8MDBSmP0opVSxKLqB7vbB4cXq75tGVUpWu5AI6aNpFKaWcaEBXSqkyoQFdKaXKREkG9JYWkHHbanR1wdBQYfqjlFLFoCQDelUVLFiQ3q4XRpVSlawkAzo4p13a2/PfD6WUKhYlG9Db2tLbNKArpSpZyQb0ZcvS244ezXs3lFKqaJRsQG9tTb8w2t0NwWBh+qOUUoVWsgHd74dFi9LbNe2ilKpUJRvQwTntogFdKVWpSjqg64VRpZQaVdIBXS+MKqXUqJIO6EuXOpfS7e0tSHeUUqqgSjqge73Oe4zqKF0pVYlKOqADLF+e3nbkSP77oZRShVaWAV1H6EqpSlS2AT2RyHtXlFKqoEo+oC9aZKsvpgqH4dSpwvRHKaUKpeQDuojzfHTNoyulKk3JB3TQC6NKKQVlHNAPH85/P5RSqpCyCugicouI7BORgyLyWYfHG0TkZyLyWxHZJSIfyn1XM1uxIr3txAmbS1dKqUoxaUAXETdwL3ArsAa4U0TWjDvto8BuY8w64Hrg70XEl+O+ZlRXB/Pnj20zRtMuSqnKks0I/UrgoDHmsDEmAjwI3D7uHAPUiYgAtcA5IJbTnk7ivPPS2w4dymcPlFKqsLIJ6M1A6vbLncm2VN8GLgROADuAjxtj0maCi8hdIrJFRLZ0dXVNs8vOnNIumkdXSlWSbAK6OLSZcfffDGwHlgCXAt8Wkfq0XzLmfmPMBmPMhvnjcyQz5DRCP3zYpl6UUqoSZBPQO4GlKfdbsCPxVB8CHjXWQeAIsDo3XczOkiV2F6NUoRCcPJnPXiilVOFkE9A3AytFZHnyQud7gMfGnXMMuBFARBYCFwB5TXi4XM7TFzWPrpSqFJMGdGNMDPgY8CSwB3jIGLNLRD4iIh9JnvZ3wDUisgN4CviMMaZ7tjqdiVPa5eDBfPdCKaUKw5PNScaYJ4AnxrXdl/LzCeDm3HZt6s4/P71NA7pSqlKUxUrREStW2Nouqbq7dQcjpVRlKKuAXlUFLS3p7TpKV0pVgrIK6KBpF6VU5dKArpRSZaIiAnpnp52TrpRS5azsAnpjI8ybN7bNGJ2PrpQqf2UX0AFWrUpv278///1QSql80oCulFJlomICens7DA/nvy9KKZUvZRnQ586FOXPGtiUSmkdXSpW3sgzooGkXpVTl0YCulFJlomwD+gUXpLcdOaIbRyulylfZBvS5c+2RKpHQVaNKqfJVtgFdxHmUvm9f/vuilFL5ULYBHZwD+t69+e+HUkrlQ8UF9I4OCAbz3xellJptZR3Qm5pg4cKxbcbobBelVHkq64AOzqP03bvz3w+llJptZR/QL7wwvW3Pnvz3QymlZlvZB/TVq9P3Ge3qsnuNKqVUOSn7gB4IQFtberuO0pVS5absAzo4p100j66UKjcVEdDXrElv27vXrhxVSqlyUREBffly8PvHtg0NwdGjBemOUkrNiooI6B6P8/TFXbvy3xellJotFRHQAS66KL1t587890MppWZLxQT0tWvT29rbYXAw/31RSqnZUDEBfd485zIAOn1RKVUuKiagg6ZdlFLlraICulPaZedOnb6olCoPFRXQV60Cn29s2+CgTl9USpWHrAK6iNwiIvtE5KCIfDbDOdeLyHYR2SUiz+a2m7nh9draLuO9+mr++6KUUrk2aUAXETdwL3ArsAa4U0TWjDunEfjfwNuNMWuBd+W+q7lxySXpbTt25L8fSimVa9mM0K8EDhpjDhtjIsCDwO3jznkv8Kgx5hiAMeZMbruZOxdfnN7W2QnnzuW/L0oplUvZBPRmoCPlfmeyLdUqoElEnhGRV0TkA05PJCJ3icgWEdnS1dU1vR7PUGMjLF2a3q5pF6VUqcsmoItDmxl33wNcDtwGvBn4WxFZlfZLxtxvjNlgjNkwf/78KXc2V5zSLtu3570bSimVU9kE9E4gdUzbApxwOOcXxpigMaYbeA5Yl5su5t46h57t3w+hUP77opRSuZJNQN8MrBSR5SLiA94DPDbunJ8C14mIR0QCwEagaNdgtrba1EuqeFwXGSmlStukAd0YEwM+BjyJDdIPGWN2ichHROQjyXP2AL8AXgVeBh4wxhRteBSBSy9Nb9e0i1KqlHmyOckY8wTwxLi2+8bd/wbwjdx1bXatWwfPPDO2bedOiMVsuV2llCo1FbVSNNWqVVBVNbYtHLY7GSmlVCmq2IDu8TjPdtm6Nf99UUqpXKjYgA5w2WXpbdu32wukSilVaio6oK9dm16sKxi0UxiVUqrUVHRA9/mca6Rr2kUpVYoqOqCDc9pl2zatka6UKj0VH9AvvtiW1U01MAD79hWmP0opNV0VH9CrqpzTLlu25L8vSik1ExUf0AE2bEhv27bNLjJSSqlSoQEdm3Zxmu2ii4yUUqVEAzrg9zsvMnr55fz3RSmlpksDetIVV6S3bdsGw8P574tSSk2HBvSkiy6CQGBsWySiFRiVUqVDA3qSxwOXX57ermkXpVSp0ICeYuPG9Lbdu6G/P/99UUqpqdKAnuL886GpaWxbIqGjdKVUadCAnkIErroqvX3TJjDjt8VWSqkiowF9nKuvTm87fhw6O/PfF6WUmgoN6OMsXAgrVqS3b9qU/74opdRUaEB34DRKf+klLQWglCpuGtAdbNiQXoExGNQ56Uqp4qYB3UEg4Fwn/fnn898XpZTKlgb0DF73uvS2vXuhuzv/fVFKqWxoQM9g5UpYsCC9/de/zn9flFIqGxrQMxBxHqX/+td6cVQpVZw0oE/gmmvA7R7bNjBgqzAqpVSx0YA+gbo654ujzz6b/74opdRkNKBP4g1vSG87cMCuHlVKqWKiAX0S558Pixent//qV/nvi1JKTUQD+iRE4Prr09tffNEuNlJKqWKhAT0LV18NVVVj26JR+M1vCtMfpZRyklVAF5FbRGSfiBwUkc9OcN4VIhIXkXfmrouF5/c7T2F8+mmIx/PfH6WUcjJpQBcRN3AvcCuwBrhTRNZkOO9rwJO57mQxuP56m35J1dMDW7cWpDtKKZUmmxH6lcBBY8xhY0wEeBC43eG8PwN+BJzJYf+Kxvz5cMkl6e3/9V+6+YVSqjhkE9CbgY6U+53JtteISDPwu8B9Ez2RiNwlIltEZEtXV9dU+1pwN92U3tbebqcxKqVUoWUT0MWhbfyY9G7gM8aYCTPKxpj7jTEbjDEb5s+fn2UXi8f550NbW3r7k2WZZFJKlZpsAnonsDTlfgtwYtw5G4AHReQo8E7gf4vI7+Sig8VExHmUvnMndHSktyulVD5lE9A3AytFZLmI+ID3AI+lnmCMWW6MaTPGtAGPAH9qjPlJrjtbDC67zObTX2MMrkSMX/xcE+lKqcLyTHaCMSYmIh/Dzl5xA/9qjNklIh9JPj5h3rysnD2La/du3p84zP5XjlId6cMbCyEmQWKTh+D2ADWtc6GlBZYuhdWrbQ3e8dNjlFJqFkwa0AGMMU8AT4xrcwzkxpgPzrxbRSSRsHvPPfec3eHCGM5PQE8UItHR01yJGJ27+7kg3g9Hjow+MG8eXHqpXZ3U0pLv3iulKkhWAb1iHTgADz4InZ1jmt0uaG6BI4fHnn7mDLS2QnV1SmN3N/zyl/ZYuhTe+EbYuBE8+tErpXJLo4qT4WEbyDdtynjKksXQ2WFLAIwwBo51wAWrMvxSRwd873vwk5/AjTfa4O7357TrSqnKpbVcxuvuhq9/fcJgDnbji+ZkBsXI6Md45jQMhSZ5jf5++PGP4fOftyN33QJJKZUDOkJPdegQ3HvvxGUU582DK6+E885jwbwW7vlqLQNDbjyJCP7oAI3BE4Qa2rl56R44fHjiZaQDA/Dww7YW7zvfaXPtegFVKTVNGtBHtLfDPfdAOOz8+LJl8Lu/a2euJINuFfCmW+xgO+b2E3P7CVbN49HoJVz0+29jSUMQtmyxo/2jRzO/dnc33Heffe73vhcWLsz521NKlT8xBSpEsmHDBrNly5aCvHaaEyfgm990Hpl7vfDud9tyi670DNXwsM2cDAyMbb/kEvjoR1Ma2ttt4ZdXXrEzZzLxeODNb4Zbb7WvrZRSKUTkFWPMBqfHNIfe3w933+0czOfOhc98Bl7/esdgDvaa5q23pre/+uq4Gi/LlsGHPwxf/KKd5ZIptRKLweOPw5e+ZFM2SimVpcoO6IkEPPAA9PWlP7ZkCfz1X9uphpN4/ethzpz09kcfdUihL1gAf/iH8Ld/CxdemPlJT52yF2cffnjsVBqllMqgsgP644/Dvn3p7QsWwCc+AXV1WT2N1wu3OxQUPnzYZlgcNTfDxz9u8zLz5jmfY4ydBfOlL02cg1dKKSo5oO/dawP6eHV1Npg3NEzp6a680sbo8R59dIIBtohNtn/hC/CWt9i5kE5OnYKvfQ2eeGLi/LtSqqJVZkCPRuH730/Ph4jYPPfcuVN+SpcL7rgjvf3sWTvInpDPZ4f4f/M3zvV5wQbyn/7UXrw9e3bK/VNKlb/KDOhPPglOG2zcdpudOjhNa9faY7yf/xx6e7N4giVL7EXYO+7IXBrg0CH4u7+z0yGVUipF5QX0M2dshB1v5Uob0GfoXe9KnxAzPAyPPJLlE7hccPPNdi5ka6vzOaEQ/PM/w//5PxCJzKi/SqnyUXkB/cEH05fau1zwvvdlnJo4FYsXwxvekN6+eTPs3z+FJxoZrd96a+Ypjr/+NXzlK3YevVKq4lVWQN+/H3btSm+/6SYbiXPk7W+H2tr09h/8YIplWzwe+J3fgb/8S+d5kQAnT9qg/pvf6G7VSlW4ygnoxsBjj6W3NzXlJNWSKhCAd7wjvf3UKfjP/5zGE65caeetX3aZ8+PRqK3i+N3vZi5doJQqe5UT0PfvH7d0M+kd75iVErbXXAMrVqS3P/44nD49jScMBOCuu+D3fz9zSYCXXrKj9XH125VSlaEyArox8LOfpbcvXgwbHEsizJiIc1o+FrOD6WllR0Tguuvgc5/DLF5E3B0lHBhksPEs/XPP0LPgBOfkFc59908599zf03Pul/T1vcDg4G8ZGjpAJNJNIqGlepUqV5VRbXHfPufR+VvfmpMLoZm0tNj0/JNPjm0/eBCefRauvz675zHGEI2eIRQ6SDh8jHC0ncidZ0nsPwKnTzn/0v4D0LsQzj9/3IIlweNpwOdbjM+3CL+/maqqVny+ZlyuyvjnoFS5qoz/g596Kr1tyRK4/PJZf+m3vhW2bk2f9v6jH8GaNbbKgBNj4gSDexgc3M7Q0G6iUYfFRKtWQWOD/YaIx9MfP3PaloG88EKoqRl5ZmKxXmKxXoaG9rx2qoiHqqpWqqtXJo/zcbur059TKVW0yj+gd3fDjh3p7bfdlpfNJHw++MAH4O//fmx7JAL/9m/wqU+N/SMhHG6nr28TAwNbiMcHJ3+BBQttuYLde2DIoWJkaAi2b4PzzodFCwHn92xMjFDoMKHQYeBJRFxUVS2npuYiamrW4ve3Irr5hlJFrfwD+rPPpies58zJPGNkFqxaZSsyPvfc2PZDh+yslze/OcHg4HZ6ep4iFDo49ReoDsD6S+HgIecUTCIBB/bb5aorV2auGZPCmASh0CFCoUN0d/8Uj6eJ2tp11NZeSiBwASKVcflFqVJS3gE9ErGLb8a7/vpZzZ07ueMO2L3b/sEwyvDcc6+wcOFjBAJTn/oi4sHrnYvHMwe3uxb3lW9EDrbD889DLIZxJUi4Y8TdMeLeKLFoL7HBrckUjMNE+QnEYj309j5Db+8zuN111NWtp67uCqqrV+rIXakiUd4BffNmGBoa2+bxwLXX5r0rVVXwwQ/a1IsxUFe3l9bWRwgEOti9G9avn3jgPJICqa5eRVXVMqqqWvF45qQH04XAipNw//3Qkb6CNCEJolujDL/tMiJrFxMePkY4fJR4fCDt3Ezi8QF6e5+jt/c5PJ4m6us3Ul9/FX5/7hZnKaWmrry3oPvKV+zWb6muuQb+4A9m93Un8Oij/Rw48BBz524e0z5/fnpdMBEvtbXrqKvbQCCwemoXKSMRW+bgN7/JfM6ll8IHPoAJBIjFzjE0dIBQaB9DQ/ucL8JOoqqqjYaGa6mr24DbHZjy7yulJjfRFnTlO0I/eTI9mEP2cwVzzBhDX9+vueiiHxGJhBgcd72zq8uWYF+8GPz+Zhobb0gGxqrpveDI1dgLLrA1B4aH08/Zvh3a25E/+iO8K1fS0DCXhoarXpsmGQzuZHBwB6HQPoyZvA57OHyUcPgoZ848RF3dBhobr6OqaoWmZJTKk/IN6C+/nN7W2mr39syzWGyA06e/x+Dgq4AdiW/blj7TcPv21SxZcivLll2QuyC4caOtsX7//c4rSHt6bB7ottvs4XIhIvh8C/H5FtLUdCPx+BDB4E4GBrYSDO7EmIm3xDMmSn//C/T3v4Dfv4SGhjdQX3/V9L+clFJZKc+UizG2/Oz4jSDe/W648cbZec0MgsHdnDr1XWKx/jHtZ86M7n4XDC6js/N36e+/kKYm2/Usd7/LXjRqa/g+80zmc847z+53mmlLPCAeDxMM/pb+/pcZGtqd1cgdwOXyU1+/kcbG6/H7HbZ2UkplZaKUS3kG9EOH7AbLqURsW3397LzmOMYYzp17ku7unwDOn/H+/TW8+OI76e6+mtT54atW2e1GM+1xMSPbttnaA+MvFo+oqoI777Qj+0n+SojFBhgY2Exf3yaGhzuy7kJ19Uqamm6gtvZSnf6o1BRVXg79xRfT2y68MG/BPJEY5tSpf2dgINMO0dDQcA033XQHO3bUjpvKaOuIPfigrQWT8/Tz+vU29fQv/2K/+MYLh23Vxh074L3vTVlhms7jqaOp6Qaamm4gHO6gr+83DAy8RDye4csiKRQ6QCh0AI+nicbG62lsvA63O/PrKKWyU34j9FgMPv1pCI5bNfmhD8FVV+X+9caJRns5fvzbGUesbncNixb9AbW16wCbFfryl9O7C3b3oze9aZY6mkjY0o+PP565Ulhjo72w6rSvXsanjTI4uJXe3ucJhRzq5zgQ8VJffxVNTTfg9y/J+rWUqkQTjdCz+ntXRG4RkX0iclBEPuvw+PtE5NXksUlE1s2009O2d296dPR67RS9WTY8fIqOjq9nDOaBwIW0tX3htWAOdj/qu+5yXuf0yCM2QzIrXC5429ts7YFMm2L39sI998D//b/Os2Qcn9ZLff1GWls/RVvbF2hsvB6Xa+KLocZE6et7nqNHv0hHx90MDu6kUAMNpUrZpCN0EXED+4GbgE5gM3CnMWZ3yjnXAHuMMT0icivw340xGyd63lkbof/gB+lr7K+4Aj784dy/VopQ6BDHj99LPO4w1AbmzHkz8+b9Tsac8TPPwA9/mN7u9cInPmGLJs6aoSGb43nppcznzJtn5++vWjXlp4/HwwwMvERPz6+IRE5m9Ts+30IaG2+goeFqXK7c16tXqlTNdIR+JXDQGHPYGBMBHgRuTz3BGLPJGNOTvPsi0DKTDk+bMfDqq+nt69fP6ssODe2ns/Nux2DucvlYvPiPmT//HRNeAHzDG5ynyEejcO+9s7xtaCBgZ7d8+MP2Zyfd3XZ64w9/mPVofYTbXUVj4xtoa/sCLS2fTP6FMvHFgUjkNGfO/JDDhz9LV9ejRKM9E56vlMouoDcDqTmEzmRbJn8E/NzpARG5S0S2iMiWrvH1ZHPh2DGbJkjldk8pBzxVweAeOjvvIZGIpD3mdtfR0vKX1NdPvomGCPze78HFF6c/NjQE//iPdqrjrLriCvhv/81eQM7kmWfgi1+EPXsyn5OBiFBTs5rm5j9l+fIv0dT0JlyuiVe/xuNDnDv3JEeOfI4TJx4gFDo65ddVqlJkE9CdhlKOeRoReSM2oH/G6XFjzP3GmA3GmA3z58/PvpfZ2r49vW31ajsVbxYEg7s4fvzbjgttvN55tLZ+murqtqyfz+WCP/5juw5ovP5+G9THT63PuaYmO2fyPe+xq02dnD0Ld99t6/86Xc3Ngs83jwUL3sV5532NBQvuxOfLUBg+yZgEAwObOXbsf3Ls2NcZGNia9Rx4pSpFNgG9E1iacr8FSEsAiMglwAPA7caY2Q47zpzSLZdcMisvNTS0n+PHv4Mx6Vu6+f1LaW39zKRByonfDx/7mPPGF+fOwT/8g72dVSLwxjfajanPOy/zeS+8AF/4gs29T/Mipsvlp6npetra/gfNzR8lEFg96e+EQoc4ceKfOHLkbzh37j8nnSapVKXI5qKoB3tR9EbgOPai6HuNMbtSzmkFngY+YIzZlM0L5/yiaHe3XWI53le/akedORQOt9PR8Q8kEuG0x6qqltHS8vEZz6s+exa+8Q27Mn+8uXPhL/5iwgWduZNIwK9+BT/+sU3oZ7J6tZ23vnDhjF8yHO6kt/cp+vtfdvzCHM+uQr2apqYb8Plm/vpKFbMZrxQVkbcAdwNu4F+NMV8WkY8AGGPuE5EHgDuAkWpYsUwvOCLnAf3pp+H//b+xba2tzkF+BoaHT9DR8U3HC6BVVctpafnznFUaPH3aBvUBh8q2jY3wyU/CokU5eanJdXXZFab792c+x+2Gm2+Gt7wlc7pmCmKx/mSZ3meyLu9bU7OWxsYbqKlZq0XBVFmqjKX/99wDu3aNbXvb2+ymnjkSjfZw7NhXicV60x6rqmqjpeWTOS9AdeKEnVwyvjoj2EWcf/7nzjn3WWGMLcf7yCMQCmU+b84cuypq/fqcLHVNJKIMDGymp+cphocdCow58PkW0Nj4Rurrr9GiYKqslH9Aj8XscDUybqbJ5z9vR+k5EI+H6Oj4BsPDx9Me8/uXsHTpp2Zt+frx4/aCqNNI3e+HP/mTWZ3Ik663Fx56CF7JXNoAsGmYd70LWnIzi9UYw9DQPnp7n2JwcAeZauSkGknH2KJgugGHKn3lH9APHIBvfnNsW02NHdrmZIQY4/jx/8XQ0N60x7ze+bS2/hUeT8OMX2ciJ0/aoN7Xl/6Yy2XT19ddN6tdSLdzp11FOtHUGxHbsbe9Lae1dCKRM/T0PE1//yYSiezmxQcCq2hsvD5ZFGzyfVWVKkblH9Afe8zWJEl1+eV2Tf0MGWM4derf6e9/Ie0xj6eepUs/g8+Xj6uTNo19993j9yUd9aY32b1L87pdaiQCP/+53e06NsEFTL8fbrnFdjIH+fUR8XiI/v5N9PQ8TTSa4YMZx+NpoKHhdTQ0vA6vd07O+qJUPpR/QP/619MrB77vffD618/4qc+de5KurkfT2l0uP0uXfoqqqtykdLLV1wff+pZNwzhZs8Yu+JygSOLsOH3apmF27pz4vIYGe13j2msn3kR1ioxJEAzupKfnaYaGsl30JNTWXkxDw3XU1FykpXxVSSjvgB4O2/x5Ytwiky99yW7UOQODg7/l+PHvMD5XK+Kiuflj1NTkM3E9KhSC++6zdciczJtn8+o5unyQPWNs2d2HH558WeuCBTYNc8UVOa8RPDx8it7eX9Hf/0LW6RiPpzE5ar8GrzdDsTKlikB5B/RXX7XFTlLNnWtr0s4gUITDnXR0fN0xICxc+H4aG1837efOhVjM1iHblGHWv8djSwlcd90s1FSfTCxm567/x3/YL9yJLFliA3uOZsSkisfD9Pe/QG/vM0Qip7L8LSEQWE1Dw+uorV2Hy+XNaZ+UmqnyDugPPQRPPTW27dprbR3vaYrHg7S3fxmnne+bmm5iwYJ3Tvu5c8kYm7r+8Y8zL9Rcvx7e//4CpGDATst5/HF49tn0v6DGW7LE7ml62WU5vwhgjCEU2k9v7zMMDm7PumSA2x2gru5KGhquwe9v1XntqiiUd0D/4hfTSxF++MP2T/lpMCZBZ+e3HGe01NRcTHPznxZdrnXnTrsBUaZd5RoabOXbvE5tTHXmDPzkJ5NPcwS70vSWW+DKK2dlD75otJf+/k309T1PNJp9DQWfbzENDVdTV3clXm9uVx4rNRXlG9D7++Gv/iq9/RvfmPYUuTNnHqan55dp7Xau+WeKdpFKVxf80z9BxwRbe153nZ0FUz1xgcPZc/SoDezZVGpsarIzYl73ulkprmYvou6ir+95gsEdUyj0JQQCF1Bfv5Ha2vW43YX6MFWlKt+A/vLLdmiaqrnZloCdhv7+lzl58l/S2t3uAK2tn8/b9MTpikZt9YPnn898TmOjLaQ4yyXiJ7Z/P/z0p3Dw4OTnVlXZb6IbbrArUGeBHbW/QF/fb4hGsy/rLOKhpuZi6uuvpKbmYs23q7wo34D+yit2DnTqsPTGG+Hd757yU4XDnRw79lWHUrhCS8vHqamZoEZ4kXnlFfj+9zOnYMAWofy938tTgS8nxsC+fXYNgdNm1eOJ2G+hG26w2zfNQj7b5toP0Nf362R53gmKkY3jcvmprV1HXd0GAoE1GtzVrCnfgD5iYMAGh717be51itukxeNDtLd/xXF0Nn/+HcyZc3Nu+plHPT3w7/8+cXbD47G1tG65xa77KQhj7Ij9iScyz8Mcb8kSu8XTVVfNWq17u23eFvr7NxEKZfGFk8LlqqKm5mLq6i6jpmatbqGncqr8A/oMGGM4fvxegsEdaY/V1V3B4sV/VLKzG4yx6ZdHHpl417j6enj72+3koLyuMh3v6FH4xS/sRiXZ/Lv0+WDDBpuSWb581uZnRiJd9Pe/RH//i1NKyQCIeKmpWUNt7Xpqay+ZtXo/qnJoQJ9Ad/d/cPbsz9La/f4ltLZ+tixGV2fP2pIrky3iXLjQBvbLLy/A3PVUp0/DL39pN9CYqAZ7qkWL4JprYONGe6FgFhhjCIeP0t//EgMDW7Iu6TtKqK4+n9raS6ipuQSfb2HJDhZU4WhAz2BwcCfHj3+b8StBXa5qli373LR2HCpWxtjc+sMPp2+7Ot6SJbak+eWXF3jEPjgIzz1n57FP1ukRIrbK41VXwaWXzlpKxpgEQ0P7GBh4mYGBbSQSE5QTzsDrnU9t7cXU1FxEdfUqzburrGhAdxCJdHPs2Jcdty9rbv4otbWzs3VdoYXD9jryL385cS0tsKvzb77ZxkZvIWNNPA7bttkNqg8cyP73vF676/aGDXDRRbN2oSCRiDE0tIeBgS0MDv52WsFdxEsgsIqamrUEAmvw+Rbp6F050oA+TiIR5dixrzE8nD5pe+7c25g37+0F6FV+dXfbFabZ/Ceoq7PXIK+7btayGdk7ccKO2F96aeJNNsbz+WxQX7/eBvlZmoxvg/teBge3Mjj4W+Jxh51JsuDxNBIIrCYQuJBAYDVeb2NuO6pKlgb0FLYc7nfp738p7bGamrU0N3+s6FaCzqYjR+xan2wmmLhcdmX+ddfBBRcUOM8eicDWrXYHpYm2xXPidtuZUJdeaoP73NkpxmVMglDoEIOD2xgc/G3W5X2d+HwLqK6+gEBgFdXVK3W1agXTgJ6ip+cpzpx5KK3d653LsmWfr9hZCHv3ws9+lt1aH7CFLEeuQc5SPMxeVxe8+KI9MhWLn8iSJXb0vnatneM+CyUHjDFEIicZHPwtweAOQqHDZLPjUiZe7zyqq1dSXX0e1dXna4qmgmhATxoa2kdn591py7xFvLS2fjrvtc2LjTE2Rf3EE9mtzh+xapUtnXPZZVBbO3v9m5QxcPiwXUG8ZYvzRqyT8flg5UpbWP6CC+z2ebMQKOPxIMHgToLBXQSDu6cxY2YstztAVdUKqqtXUFW1gqqqZTnbrFwVFw3oQDR6lvb2rzjmNBct+iANDVfnrS+loKMD/uu/YPPmyQsljnC5bAy87DJYt84WBSuYRMIuNtuyxV5QDQan9zw1NfYba+VKe9vcPCvVIIeHjxEM7mZoaDeh0CGMic/4eX2+RVRVtSWPZfj9LbhcudstShVGxQf0RCJKR8fXCYePpT3W2PhGFi58T176UYp6e+3Mweeft7XQpqKtzZYYuOgiu9lGwTICiYTNs2/bZhctZTsF0klVFZx3HqxYYW/b2nJ+gTWRGCYUOsjQ0F6CwT0MD3cyk/TMKMHvX4zf34rfv5SqqqX4/S0Vm2YsVRUf0Ht7n+P06R+ktQcCq2hp+YRuGJyFWMzuJfLcczbfPtV/NnV1dnr4hRfaUXxBa8gcO2bfzI4d0N4+s+cTsYua2trssWyZTdPkcJ5nPD5EKHSAoaH9hEIHkgOT3P1/6/E04ve34Pc34/MtSd4u0nnxRariA7oxhp6ep+jqeoSR/xE8niaWLfs8Hk9dXvpQTs6dG70Gefr09J5jzhybxVi50g50Fy8u0Ai+rw9277bLaPfsmX5qJpXLZYN8a6sN7iNHXW7+rcXjYcLhQ4RChwiFDhIOHyGRiOTkuUcJXu88/P7F+HyLxhw6oi+sig/oI4LBPZw8+c8kEsPJi6DL8vr65WZksLt5s12Fei77/SLSVFfbciwjA93WVjvnPa9BfuQN7d1rjwMHsi89kI26OjujZuRYtMh+k9XWzuiNGpNgeLiTUOgw4fBhwuGjRCLT/KbNgttdg8+3EK93AT7fArzeBXi98/F65+F21+hsm1mmAT1FJNJNJHKc2tp1eX/tcmaMzV5s326Pkydn/px1dbB06egAt7nZxsBZmFXoLBazb2r/fjuf8+DByfdInY5AwL6xBQtsQZ358+3P8+fbx6YhHg8SDre/dgwPt09ph6bpcrmq8HrnJY+5eL1z8XjmJn+eg8sV0IA/QxrQVd6dOQO7dtk09f79uRvoitg4t3ixjX2pcbCubpZH9ImEXal6+LCt4X7kyPRzTtkKBOwFh3nz7IT/uXNtvmruXLurUyCQ9Zu2Qf4Yw8PHGB7uJBzuIBo9PYXdmmbO5fLh8czB42nC623C42nE42nC42nA42nE7W7A46mvqMV9U6UBXRVUNGrj3549diZhe3v2UyGnwu+3wT419o3Ev6amGWc2nAWDtuxve7u97eiYWe5pqnw+++aammyOauRoaLBHfb09MtSxSSSiRCKnGB7uJBI5wfDwCYaHjxOL9eTvPaQRPJ6614K7x9OA212H212fbK/H7a5NttXicuXrT7bioAFdFZVw2A5uDxwYHehOVK89V7xeG+NGYt9IzBuJe3V19ra2doZTzQcGoLPTBvfOTjh+3Oag4jOfWz5tfv/oG6yrs2+yrm7055oaeyR/jle5iEROE4mcJBI5lTxOE42eyckc+VxyuaqSAb7mtVuXa+TnQPJ+ALc7MOZWxFOS6R8N6KqojWQyjh61x7FjNgZOVg1yNo3EtpFYV1trsxsjtyNHdbU9AgE7Rd3ny/BXQDxu81AnTtjj1Cl7e+ZMYd9oJiKjbyzlzZpAFdFAnGggTMQfJuoLEvEGiboGibr6MW5ja+WMHEWcOhFx43JV43JV43ZXJ3+uSrlNPfyv3Yr4k/f9Kfd9eUsTaUBXJScet4Pa48dHB7nHj89sTVA+uFw2Dl57LdxxRxa/kEjYFM3p0za4nz5ta9OcOWN3JinkqH6KDIaYN0LMP0zUHybqCxOtjhINRIlVx4hWRTEewO0CVzLgty6FmkLWi8gdEU8ysPvSbkW8uFzelJ99iHioq7uSqqqWKb5O5oCeVfJJRG4BvgW4gQeMMV8d97gkH38LMAR80BizdUq9VCqF2z06u2XjxtH2UMgObk+dsrFv5OjqsgUYCy2RsGn1rOOwyzV60XPt2vQn6+uzb+7sWXucOzd629OT22mVMyQI3qgfb9RP9WB92uMGQ9wTI+Ybtod3mGjb1cTqa4nFeonFeonH+xz3KCgFxsSIx2PYEJidqqplUw7oE5k0oItdRnkvcBPQCWwWkceMMbtTTrsVWJk8NgLfSd4qlVMj89WXLx/bboxNXXd12YKLZ8/a25G4d+5cfgN+TqoBuFyjFzydGGO/PXp77Zvs7bVHX9/o0d9vb2fjKvQUCYIn5sUT88JQclRe/xZYvGLMeYlElHi8n1isP3nbRzw+kLw/kPx5gHh8MFmbqTBZhlwQye1q3GxG6FcCB40xh20H5EHgdiA1oN8OfM/Y/M2LItIoIouNMTmYjazU5ERGJ3Scd17648bY0f1IzBuJeyPxrr/ffiH098NQDgaIs7R/xlgio4n+lglGecbYNzXyJgcH7c+Dg2OPYHD0mI359k4cZt+4XF5cLjt3fTLGGBKJodeCezweTN7aNvtYMHk7lHIbKoqLu4UI6M1A6tY+naSPvp3OaQY0oKuiIDJ6bW/JkonPjcdtTHOKdUNDo7eh0OjtyDFyfTMvAT1bIqOzWBYvzu534nH75kbeYDA4+obD4bFvOhSy05RG2kd+jkQm/8tghtsCikhydksNsDDr3zPGYEyMRGKIRCJMPB4ikQiRSISTtyESieHk/XDKz8MYM5y8P5y8HyGRmN40rVxXv8wmoDtdsx//N0425yAidwF3AbS2VnbtcVW83O7R0f5UxWI2pvlKvUqt2z06rXG6jLEfSCQyGuCHh8feTudDzgERSV6cbABmXud59AtiGGOirwV5Y6IkEpHk/WjysZG2KB5PbneHySagdwJLU+63ACemcQ7GmPuB+8HOcplST5UqAR5PzmpwlT4RO/nf67V/HZSx0S+IwlaozGbi5GZgpYgsFxEf8B7gsXHnPAZ8QKyrgD7NnyulVH5NOkI3xsRE5GPAk9hpi/9qjNklIh9JPn4f8AR2yuJB7JydD81el5VSSjnJah66MeYJbNBObbsv5WcDfDS3XVNKKTUVxbsuVyml1JRoQFdKqTKhAV0ppcpEwYpziUgXMN0deucB3TnsTjnSz2hi+vlMTj+jiRXq81lmjJnv9EDBAvpMiMiWTNXGlKWf0cT085mcfkYTK8bPR1MuSilVJjSgK6VUmSjVgH5/oTtQAvQzmph+PpPTz2hiRff5lGQOXSmlVLpSHaErpZQaRwO6UkqViZIK6CJyi4jsE5GDIvLZQven2IjIUhH5lYjsEZFdIvLxQvepGImIW0S2ich/FLovxSi549gjIrI3+W/p6kL3qdiIyCeT/4/tFJEfikhVofsEJRTQU/Y2vRVYA9wpImsK26uiEwP+0hhzIXAV8FH9jBx9HNhT6E4UsW8BvzDGrAbWoZ/VGCLSDPw5sMEYcxG2Cu17Ctsrq2QCOil7mxpjIsDI3qYqyRhz0hizNfnzAPZ/xObC9qq4iEgLcBvwQKH7UoxEpB54PfAvAMaYiDGmt6CdKk4eoFpEPEAAhw19CqGUAnqmfUuVAxFpA9YDLxW4K8XmbuDTwCSbXVasFUAX8N1kWuoBESnv7YamyBhzHPgmcAy7b3KfMeY/C9srq5QCelb7lioQkVrgR8AnjDH9he5PsRCRtwJnjDGvFLovRcwDXAZ8xxizHggCer0qhYg0YbMDy4ElQI2I/H5he2WVUkDPat/SSiciXmww/4Ex5tFC96fIXAu8XUSOYlN2N4jI9wvbpaLTCXQaY0b+snsEG+DVqDcBR4wxXcaYKPAocE2B+wSUVkDPZm/TiiYigs197jHG/EOh+1NsjDF/bYxpMca0Yf/9PG2MKYqRVbEwxpwCOkTkgmTTjcDuAnapGB0DrhKRQPL/uRspkgvHWW1BVwwy7W1a4G4Vm2uB9wM7RGR7su1zyS0ElcrWnwE/SA6cDqN7BI9hjHlJRB4BtmJnlm2jSMoA6NJ/pZQqE6WUclFKKTUBDehKKVUmNKArpVSZ0ICulFJlQgO6UkqVCQ3oSilVJjSgK6VUmfj/3bYbIVX+oswAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "\n",
    "a = 1\n",
    "\n",
    "x = np.linspace(gamma.ppf(0.01, a), gamma.ppf(0.99, a), 100)\n",
    "ax.plot(x, gamma.pdf(x, a),'b-', lw=5, alpha=0.6, label='a = 1')\n",
    "\n",
    "a = 2\n",
    "\n",
    "x = np.linspace(gamma.ppf(0.01, a), gamma.ppf(0.99, a), 100)\n",
    "ax.plot(x, gamma.pdf(x, a),'r-', lw=5, alpha=0.6, label='a = 2')\n",
    "\n",
    "a = 3\n",
    "\n",
    "x = np.linspace(gamma.ppf(0.01, a), gamma.ppf(0.99, a), 100)\n",
    "ax.plot(x, gamma.pdf(x, a),'y-', lw=5, alpha=0.6, label='a = 3')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95302231",
   "metadata": {},
   "source": [
    "Next, we need to pick a likelihood function that updates our estimate of the rate based on observed wait times. We assume buses arrive as a Poisson process with a certain rate, in other words, our likelihood function is\n",
    "\n",
    "p(t | $\\lambda$) = $\\lambda t e^{-\\lambda t}$\n",
    "\n",
    "plot this as a function of $\\lambda$. You can use either a built in function or write your own. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "id": "1e652bb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def poisson(rate, time):\n",
    "    return rate*time*np.exp(-rate*time)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69042209",
   "metadata": {},
   "source": [
    "Next, we use Bayes' theorem to calculate the posterior estimate of the rate at which buses arrive based on the time we have waited.\n",
    "\n",
    "p($\\lambda$ | t) $\\propto$ p(t|$\\lambda$) p($\\lambda$)\n",
    "\n",
    "(note that we are ignoring the denominator here - you can calculate it by summing over the right hand side for all values of $\\lambda$ - it is just a normalization constant.)\n",
    "\n",
    "Plot your posteriors as a function of the unknown parameter $\\lambda$ for many different values of the observed wait time. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "74b4e9bf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fec904fc970>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEICAYAAABMGMOEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAwTklEQVR4nO3dd3yW5bnA8d+VvSFkAUlYEpEpIwJCFakiQ4uVtlaPo9pS3NW2ttrjOcdje06Pp8NqrZZjGYq7bq0stQKiMhKG7CErYSWEmYTs6/zxvMEACXnf5F3Je30/n3zedd/PcxFjrtxbVBVjjDGmOWGBDsAYY0zbYAnDGGOMWyxhGGOMcYslDGOMMW6xhGGMMcYtEYEOwJdSU1O1R48egQ7DGGPajPz8/EOqmtbYZ+06YfTo0YO8vLxAh2GMMW2GiOxu6jPrkjLGGOMWSxjGGGPcYgnDGGOMW9r1GIYxxjRUXV1NYWEhFRUVgQ4l4GJiYsjKyiIyMtLtOpYwjDEho7CwkMTERHr06IGIBDqcgFFVSkpKKCwspGfPnm7Xsy4pY0zIqKioICUlJaSTBYCIkJKS4nFLyxKGMSakhHqyqNeS74MljDPU1tXy7JfP8vnezwMdijHGBBVLGGcIDwvnuQ3P8UnBJ4EOxRjTDu3atYsBAwa4XX7z5s0MHjyYIUOG8NVXX/kwsuZZwmhE1/iu7C/bH+gwjDGGd955h2uuuYbVq1dz3nnnNVteVamrq/NJLJYwGtEloQt7S/cGOgxjTDtVU1PDD37wAwYNGsR3v/tdysvLyc/PZ8yYMQwbNozx48ezf/9+5s6dyxNPPMGMGTMYO3YsAI8//jgDBgxgwIABPPHEE4DTaunbty933XUXQ4cOpaCggN///vdcdNFFDBo0iEceecQrcdu02kZkJmSy8sBKVNUGyIxppx59fwMb9x336jX7dU3ikW/1b7bcli1bmDlzJqNHj+aHP/whTz/9NG+//TbvvvsuaWlpvPbaazz88MPMmjWLO+64g4SEBB544AHy8/OZPXs2y5cvR1UZMWIEY8aMITk5mS1btjB79myeeeYZFi5cyLZt21ixYgWqyuTJk1myZAmXXnppq/59ljAa0SW+C2XVZRyvOk6H6A6BDscY085kZ2czevRoAG666SZ++9vfsn79esaNGwdAbW0tXbp0Oave0qVLufbaa4mPjwdgypQpfPrpp0yePJnu3bszcuRIABYuXMjChQsZMmQIAKWlpWzbts0Shi9kJmQCsK90nyUMY9opd1oCvnJmz0ViYiL9+/fniy++OGc9VW3ys/okUl/uV7/6FbfffnvrAj2DjWE0okuCk9n3le0LcCTGmPZoz549p5LDK6+8wsiRIykuLj71XnV1NRs2bDir3qWXXso777xDeXk5ZWVlvP3221xyySVnlRs/fjyzZs2itLQUgL1791JUVNTquK2F0Yiu8V0Bp4VhjDHe1rdvX55//nluv/12cnJyuPfeexk/fjw/+clPOHbsGDU1Ndx///307396K2jo0KHceuutDB8+HICpU6cyZMgQdu3adVq5K6+8kk2bNnHxxRcDkJCQwIsvvkh6enqr4pZzNXHautzcXG3JAUqqyoiXR/CdnO/w4PAHfRCZMSYQNm3aRN++fQMdRtBo7PshIvmqmttYeeuSaoSI2FoMY4w5gyWMJnRJ6GJdUsYY04AljCZkJmTaoLcxxjRgCaMJXeK7cKzyGGXVZYEOxRhjgoIljCZ0TbCZUsYY05AljCbUJwwb+DbGGIcljCbUr8WwTQiNMd42f/58+vTpQ+/evXnsscfO+vyll15i0KBBDBo0iFGjRrF27doARHk2W7jXhJTYFKLCothfai0MY4z31NbWcvfdd/Phhx+SlZXFRRddxOTJk+nXr9+pMj179mTx4sUkJyczb948pk2bxvLlywMYtcNaGE0IkzBnaq3NlDLGeNGKFSvo3bs3vXr1Iioqiuuvv5533333tDKjRo0iOTkZgJEjR1JYWBiIUM9iLYxz6BJvazGMabfmPQQH1nn3mp0HwsSzu5ga2rt3L9nZ2adeZ2VlnbP1MHPmTCZOnOi1EFvDEsY5ZCZksqhgUaDDMMa0I41tx9TUuTuffPIJM2fOZOnSpb4Oyy2WMM6hS3wXSipKqKipICYiJtDhGGO8qZmWgK9kZWVRUFBw6nVhYSFdu3Y9q9yXX37J1KlTmTdvHikpKf4MsUk2hnEONrXWGONtF110Edu2bWPnzp1UVVXx6quvMnny5NPK7NmzhylTpvDCCy9w/vnnByjSs/k1YYjIBBHZIiLbReShRj6/UUS+dH19LiIXulvXF04lDJspZYzxkoiICP7yl78wfvx4+vbty3XXXUf//v2ZPn0606dPB+DXv/41JSUl3HXXXQwePJjc3EY3j/U7v3VJiUg48DQwDigEVorIe6q6sUGxncAYVT0iIhOBZ4ERbtb1ulNrMcpsLYYxxnsmTZrEpEmTTnvvjjvuOPV8xowZzJgxw99hNcufLYzhwHZV3aGqVcCrwDUNC6jq56p6xPVyGZDlbl1fSItLI1zCrYVhjDH4N2FkAgUNXhe63mvKj4B5ntYVkWkikiciecXFxa0IFyLCIsiIy7AxDGOMwb8Jo7F5Y40e9yciY3ESRv1xd27XVdVnVTVXVXPT0tJaFGhDqXGpHDp5qNXXMcaYts6fCaMQyG7wOgs4a1WciAwCZgDXqGqJJ3V9ITXGEoYxxoB/E8ZKIEdEeopIFHA98F7DAiLSDXgLuFlVt3pS11dSY1MpOVnSfEFjjGnn/DZLSlVrROQeYAEQDsxS1Q0icofr8+nAfwApwDOulY81ru6lRuv6I+7U2FSOVB6huq6ayLBIf9zSGGOCkl9XeqvqXGDuGe9Nb/B8KjDV3br+kBLrrLA8fPIwGfEZ/r69MaYdOnr0KFOnTmX9+vWICLNmzeLiiy8OdFjNsq1BmpEamwrAoYpDljCMMV5x3333MWHCBN544w2qqqooLy8PdEhusYTRjPqEYeMYxhhvOH78OEuWLOG5554DICoqiqioqMAG5SZLGM041cKwmVLGtCv/u+J/2Xx4s1eveUGnC3hw+IPnLLNjxw7S0tK47bbbWLt2LcOGDePJJ58kPj7eq7H4gm0+2AxLGMYYb6qpqWHVqlXceeedrF69mvj4+EaPaQ1G1sJoRlR4FElRSRSXt27VuDEmuDTXEvCVrKwssrKyGDFiBADf/e5320zCsBaGG1JjUympsDEMY0zrde7cmezsbLZs2QLAxx9/fNp53sHMWhhuSI211d7GGO956qmnuPHGG6mqqqJXr17Mnj070CG5xRKGG1JiU1h/aH2gwzDGtBODBw8mLy8v0GF4zLqk3GAtDGOMsYThltTYVE7WnKS8um0srjHGGF+whOGGtFhnm3RrZRjT9qk2ejJCyGnJ98EShhvq95OyhGFM2xYTE0NJSUnIJw1VpaSkhJiYGI/q2aC3G+oX7xWftLUYxrRlWVlZFBYW0trTONuDmJgYsrKymi/YgCUMN9hqb2Pah8jISHr27BnoMNos65JyQ8fojoRLuG1AaIwJaZYw3BAmYaTEpFgLwxgT0ixhuCkl1hKGMSa0WcJwky3eM8aEOksYbkqLS7MxDGNMSLOE4aaUmBRKKkqo07pAh2KMMQFhCcNNqbGp1GotRyuPBjoUY4wJCEsYbjq1eM8OUjLGhChLGG6qTxg2jmGMCVWWMNx0arV3hc2UMsaEJo8ThojEi0i4L4IJZrY9iDEm1DWbMEQkTET+RUQ+EJEiYDOwX0Q2iMjvRSTH92EGXlxkHLERsZYwjDEhy50WxifAecCvgM6qmq2q6cAlwDLgMRG5yYcxBo202DQOlVvCMMaEJnd2q71CVavPfFNVDwNvAm+KSKTXIwtCqbGplFTYoLcxJjQ128KoTxYikiMis0Tk6abKtHcpsSl2JoYxJmR5Muj9AvA6TlcUIjJAROb4JKogZftJGWNCmScJI0xV5wG1AKq6Hhjgk6iCVGpsKieqTlBZWxnoUIwxxu88SRj7RKQnoAAiIkCsT6IKUrZ4zxgTyjxJGPcDfwM6i8htwKvAel8EFaxsLYYxJpS5faa3qu4SkQnAtcAgYDEwy1eBBaNT+0nZwLcxJgR5utJbVPV1Vf13VX1GVSs8qiwyQUS2iMh2EXmokc8vEJEvRKRSRB4447NdIrJORNaISJ6HcXuFdUkZY0KZ2y0MEZkBTBGRMmAf8CXwpao+5Wb9cOBpYBxQCKwUkfdUdWODYoeBnwDfbuIyY1U1YP1ByTHJCGJdUsaYkORJC+MSIENVs4EpwNtAvAf1hwPbVXWHqlbhjIFc07CAqhap6kogKNd1RIZFkhyTbAnDGBOSPEkYy4BkAFXdq6pzVfUxD+pnAgUNXhe63nOXAgtFJF9EpjVVSESmiUieiOQVF3t/rMHWYhhjQpUnCeNZYLGIPCAil4hIBw/vJY28px7UH62qQ4GJwN0icmljhVT1WVXNVdXctLQ0D0NsXmpsqo1hGGNCkicJ40Xg7zjjHncBn4vIVx7ULwSyG7zOwhkLcYuq7nM9FuF0hw334N5ekxqbarOkjDEhye1Bb6BQVR9p+IaIRHtQfyWQ41r8txe4HvgXdyqKSDzOSvMTrudXAr/24N5ekxKbwqGTh1BVnLWLxhgTGjxJGGtE5D5VfbL+DVV1e48MVa0RkXuABUA4MEtVN4jIHa7Pp4tIZyAPSALqROR+oB+QCrzt+gUdAbysqvM9iN1rUmNSqa6r5njVcTpEe9orZ4wxbZcnCSMDuEJEHgRWAWuBNar6ursXUNW5wNwz3pve4PkBnK6qMx0HLvQgVp9Ji3PGRUpOlljCMMaEFLfHMFT1OlXtC/QE/gPYCozwVWDByrYHMcaEqmZbGCLyiKo+KiKjgHWqegKnhbHK59EFoZTYFMC2BzHGhB53uqQWuh7vBwaISASwka9Xer/ho9iCkrUwjDGhqtmEoapfuB6vg1Mzo/oDA4GRQEgljMTIRKLComwthjEm5Hgy6A04M6NE5DJVfRxARPqo6hbvhxacRMRWextjQpJHCUNEOgJ/Ai4QkQqcmVJTgdu8H1rwSo2zhGGMCT3uDHpPBNaq6j5VPQrcJiLjgUM452K85dsQg09qTCp7TuwJdBjGGONX7rQwvgP8RkQygM3AGpyWxV5gjqrW+i684JQam8rqotWBDsMYY/yq2XUYqjpVVXOBP+KsvdgJjAVWALt9G15wSo1N5UjlEarrgnIXdmOM8QlPxjBuU9VTq61F5BngF94PKfjVr8U4fPIwGfEZAY7GGGP8w5Pdao+LyLD6F6qaD5zv/ZCCX1qssz3IoQob+DbGhA5PWhg/Al4QkY1APs46jJDskzm1eK/8EKQEOBhjjPETT/aS2gqMwtk8MAPYBEzyUVxBzVZ7G2NCkdstDBHJAR4CKlT1bt+FFPzqxzAsYRhjQoknYxgv4GwDcgmAiAwQkTk+iSrIRYVHkRSVZAnDGBNSPEkYYao6D6gFUNX1wACfRNUGpMWmUVJh+0kZY0KHJwljn+t4VQUQ5/i7WJ9E1QakxqZSVF4U6DCMMcZvPEkY9wN/AzqLyG3Aq8B6XwTVFqTHpVvCMMaEFLcHvVV1l4hMAL6Nc1zqYmCWj+IKehnxGRSXF1OndYSJJ3nXGGPaJk9mSX0OPOw6MCmkzsBoTEZcBjVaw+GKw6em2RpjTHvmyZ/G04B7RORjERnpq4DaivS4dAAOlh0McCTGGOMfnizcW6+q38HZP+rfROQfIjLYZ5EFufo9pA6WW8IwxoSGlnS+bwd+AxQCed4Np+3IiLOEYYwJLZ6MYfwTyAEqgI2ur1t9E1bw6xTTiYiwCOuSMsaEDLcShoiEAYuAq1T1pE8jaiPCJIz02HRrYRhjQoZbXVKqWgeMsWRxuoz4DFuLYYwJGZ6MYawRkUdcrQ2DM45hLQxjTKjw5Jd/NnA9zhYh74rIb0Tkez6Kq01Ij0vnYNlBVDXQoRhjjM95stL7OgARiQb64xygNAJ43TehBb+MuAwqais4XnWcDtEdAh2OMcb4VEvOwzipqvcAq3wWVRvRcC2GJQxjTHvXkvMwLoXQPg+j3qm1GDa11hgTAuw8jFawxXvGmFBi52G0QmpcKoLY1FpjTEhwewwD5zyMGXx9HsYEQvg8DIDIsEhSY1OthWGMCQmebD64CydJ/ATohbPy+2ZPbiYiE0Rki4hsF5GHGvn8AhH5QkQqReQBT+oGSv3UWmOMae/cThgishiIc52HsQ+IBuo8qB8OPA1MBPoBN4hIvzOKHcZJSH9oQd2AsMV7xphQ4ckYRkdVPS4iw4CpQDLOka3uGg5sV9UdqlqFc8TrNQ0LqGqRqq4Eqj2tGygZ8ZYwjDGhwZOEUS0iEcAtwP+q6iM4C/jclQkUNHhd6HrPq3VFZJqI5IlIXnFxsQfhtUx6XDonqk5QXl3u83sZY0wgeZIw/gysBa4G3ne9l+BBfWnkPXf31HC7rqo+q6q5qpqblpbmdnAtZVNrjTGhwpNB7zk4W4EMUNWTItIb+MKDexXi7EdVLwtnLMTXdX2qc3xnAJtaa4xp9zzaeVZVS+u3OFfV7ap6mwfVVwI5ItJTRKJwNjJ8zw91fcpaGMaYUOHJXlKXAv8FpOB0TT2hqivcra+qNSJyD7AACAdmqeoGEbnD9fl0EemMc+xrElAnIvcD/VyD7WfVdffevpQelw7Y9iDGmPbPk4V7s4A7gTXAMOAJEXlCVf/u7gVUdS4w94z3pjd4fgCnu8mtusEgJiKGDtEdrIVhjGn3PEkYh1T1Q9fz+SKyFFgGuJ0w2itbi2GMCQXNjmGIyBxX19BSEfkP19RagEqgwpfBtRUZcRnWJWWMaffcGfSeiTOFtRPwbWC7iHwEbCYIu4gCoWtCV/aVBcWkLWOM8Zlmu6RUdTGwuP61a5uOfsBgYJDPImtDshKyOFZ5jONVx0mKSgp0OMYY4xOezJLqBPwUSAc2AnNU9QVfBdaWZCU64/SFJwrplxIUW1wZY4zXebIO41XgBM4q7zicMY3hPomqjclOdNYUFp4oDHAkxhjjO57Mkuqiqr9zPf+HiLwGvAyM9H5YbUt9C6PgREEzJY0xpu3ypIVxWEROjVmo6g6clkbIi4+Mp1NMJwpLrYVhjGm/PGlh3A68ISKfAutwdqr9yidRtUFZCVnWwjDGtGuebD64GRgKfIIz8L0GuME3YbU9WYlZNoZhjGnXmm1hiEi3M95a5voCSBc5tfP4UVU97sXY2pSsxCwW7FpAdV01kWGRgQ7HGGO8zp0uqedxFu41PJOisdfPAXO8Flkbk52YTa3WcqD0ANlJ2c1XMMaYNsadhXtj/RFIm1VTCTs/JevoXsCZKWUJwxjTHnky6G0a2psPy6bD1vlQeZzs8HDolknhqlnQoTckpAc6QmOM8SqPEoaIZOPMjhoADAT6q2quLwILagc3wvOTISwC+k6Gvt8ira6GqBX/SsGOD2HDh3DLO9B1SKAjNcYYr3Fn0Pt24Ac4+0dFAx8A63FOvPtvn0YXjMpK4JXvQ1QCTPsEkroCznSzzC3PUpiRDpuXw5xr4OZ3IHNoQMM1xhhvcWda7a9w9pAaBvwDiME58e5NVd3qy+CCTk0V/P0WOHEQrn/5VLKol52YTWH1Mbj1A4jpCHO+7XRdGWNMO+BOwrhaVZer6leq+j3gL8D7IvJTEfHoTPA2b+G/we6lcM3TkDXsrI/rF+9ph2wnacQlwwvXwtE9AQjWGGO8q9lf+Kq6/ozX84HhOOdjfOajuILP0T2wcgZcNBUGfa/RItmJ2ZRVl3G08ih0zHa6pOrq4K1pUFvj13CNMcbb3DlxT858T1UrVfXfccY2Gi3T7nzxDIjAN37aZJGzNiHs1BOufhz2fAFLH/dHlMYY4zPudCl9IiL3nrniW0SigCwReR5X4mi3yg/Dqjkw8HvQIavJYo1ucz7oOhh4HSx6DPYs93WkxhjjM+4kjAlALfCKiOwXkY0ishPYhrOX1J9U9Tkfxhh4eTOhugxG3XvOYpkJmUAj25xf9Qcn0bw1FSpLfRWlMcb4lDtjGBWq+oyqjga6AZcDQ1S1u6r+WFXX+DrIgKqugOX/B73HQUb/cxaNiYghPTb97G3OYzrAtf/njIMs+V3jlY0xJsi5PctJRCYCnwKLgGdFJDQOTlr7CpQVw+j73CqeldjENufdL4YhN8EXT0PRJi8HaYwxvufJtNhngJ/jnLD3LPAHEWnf25urwrJnnBXbPb7hVpUmEwbAFY86C/4+eMC5tjHGtCGeJIyDqvqZqh5R1Y+A8cDDPoorOBRthENbYegtzgwpN/Tq0Iui8iKOVzWy03t8Klzxn85aji//7t1YjTHGxzxJGLtE5L9cs6MAqoETPogpeGx6HxC44Gq3q+Qk5wCw/cj2xgsM/QFkDoOFD0PFMS8EaYwx/uFJwlBgClAgIkuB7cAiEcnxSWTBYON70O1ij3aePT/5fAC2Hmli15SwMLjqj1B2CJb8wRtRGmOMX3hyROsNqtoP6A7cDzwKxAMzRKT9HWZd8hUUbYC+3/KoWkZcBolRiWw7sq3pQl2HwIU3wPLpcHhnKwM1xhj/8HgvKNc02zxVnamqP1HVMara/k4M2vSe8+hhwhARcjrmNN3CqHf5fzjbo3/0SAsDNMYY/wqtzQM9sel9pyXQ0fNceH7y+Ww/uh0910yopC4w+n7Y+C7s/rzlcRpjjJ9YwmjMsUJnW/K+k1tUPSc5h9LqUvaX7T93wVH3QmJXWPCvziaFxhgTxCxhNGbTP5zHFiaMZge+60XFwRWPwL7VsM6m2RpjgpsljMZseh/S+kJq7xZV793RqXfOge96A6+DLoPh499A9ckW3c8YY/zBrwlDRCaIyBYR2S4iDzXyuYjIn12ffykiQxt8tktE1onIGhHJ81mQVWVwcD30a1nrAiAhKoHMhMzmWxjgTLMd/99wvNDZNsQYY4JUs2d6e4uIhANPA+OAQmCliLynqhsbFJsI5Li+RgB/dT3WG6uqh3waaFQ8PLANaipadZmc5Bz3WhjgbDvS5ypY+idnVbkH6z6MMcZf/NnCGA5sV9UdqloFvApcc0aZa4A56lgGdBSRLn6M0RERBTFJrbpETsccdh3fRVVtlXsVxj3qJKlF/9Oq+xpjjK/4M2FkAg0X+BW63nO3jAILRSRfRKb5LEovOT/5fGq1lh3HdrhXITUHcn8I+c9B0WafxmaMMS3hz4TR2O59Zy5UOFeZ0ao6FKfb6m4RubTRm4hME5E8EckrLi5uebStVD9Tyu1uKYAxD0FUojPN1nazNcYEGX8mjEKg4Sq4LGCfu2VUtf6xCHgbp4vrLKr6rKrmqmpuWlqal0L3XLekbkSFRbk38F0vPgUuexC++hi2fei74IwxpgX8mTBWAjki0tO14+31wHtnlHkPuMU1W2okcExV94tIvIgkAohIPHAlsN6PsXssIiyCnh16sfbgJlbsPMy8dft5Pa+ADzceJH/3YQoOlze+EvyiH0NKb6eVUVvt/8CNMaYJfpslpao1InIPsAAIB2ap6gYRucP1+XRgLjAJZyfccuA2V/UM4G1xzqSIAF5W1fn+it0Tx8qrWbDxAIu3FrPlSAJ10Ru57tMvGi3bMS6SIdkdGdotmXH9M+iTkYhERMGV/w2vfB9WzoCRd/r5X2CMMY2Tc+531Mbl5uZqXp7vlmw0dLS8illLdzL7s12cqKwhPTGa7j3z2Vz9Iv817FXOT8kmPiqCoyerOFxWxb6jFawtOMrqgiNsKypFFXqlxXP1wC58b1gW2R/cCPtWwb2rna4qY4zxAxHJV9Xcxj7zWwujvaqrU2Ys3cFTH2/nRGUNEwd05o4x5zEoqwObD2dy3T9eJDxuJ3279AegG3Gn6v7LiG4AFJ+oZMGGA8xdt5+/fLKdv3yynZvPu4n/rFwCH/8amfxkQP5txhjTkCWMVjhUWsnP/r6WJVuLufyCdB4Y34e+Xb5ev3F+8vkkRiaSfzCfq3s1fWpfWmI0N43szk0ju7Pv6EleWr6bV1YUkFU9nh+tep4VyVcx/BvjEDePiTXGGF+whNFCy3aUcO8rqzl+spr/vnYA/zK821m/0MPDwhmSMYS8A+53i3XtGMsvxl/Avd/M4YO8LEoWLiP+w1/wrfwn+dmEvoztk26JwxgTELb5YAss2lLELbNWkBgTwTt3j+bGEd2b/CWem5HLruO7OHTSsx1NYiLD+c7Ffel07R8YELaLK8r/wQ+fy+OGvy1jbcFRL/wrjDHGM5YwPPTJ5iKmzcknJz2Bt+4cdVoXVGOGZQwDIP9gfovuFz5wCvS6jPvkNX43IYNtB0u55unPuO/V1ew/ZrvbGmP8xxKGB/65+SC3v5BPn86JvDR1BB3jopqt0zelL7ERsR51S51GBCb9Eamp4Lqip1j0i8u4Z2xv5q0/wDf/sJg/f7yNiurall3bGGM8YAnDTRv2HeOul1ZxQZdEXvyRe8kCIDIsksFpg8kvalkLA3DO5RjzS9j4Dok75vLA+D58/LMxXNYnjcc/3Mq4Py3mo40HW359Y4xxgyUMN5SUVjJtTj7JcVHM/MFFdIiL9Kh+budcth3ZxtGKoy0PYvT90HkQfPBzKCshu1Mcf71pGC//eAQxEeFMnZPHD59bye6SspbfwxhjzsESRjOqa+u4++VVHCqt5P9uHkZaYrTH16gfx1hVtKrlgYRHwrefgZNHYP6Dp94edV4qc++7hIcn9WX5jhLG/WkJT35k3VTGGO+zhNGM387dxLIdh/mfKQMZlNWxRdcYmDqQqLAo8g62ctV554FwyQOw7vWvzx0HIsPD+PGlvfjnA5cxvn9n/vTRVsY/sYRFW4padz9jjGnAEsY5LN5azOzPdnHrqB5MGZrV4utEhUcxKG1Qi2dKneaSnztdU+/dC8dP3+w3IymGp24YwktTRxAeJtw6eyV3vJDPvqM2m8oY03qWMJpwrLyaB9/4kpz0BB6aeEGrr5fbOZfNhzdzrPJY6y4UEQXfnQU1lfDmj6Hu7K6n0b1TmXffJfxifB8WbS3i8j8u5plF26mqqWvdvY0xIc0SRhMefX8DxaWVPH7dYGIiw1t9vcuyL6NO6/hwtxfOuUjNgav+ALuXwqd/bLRIdEQ4d4/tzUc/G8MlOan8bv4WJjyxhCVbA3eolDGmbbOE0Yj56w/w1uq93DO2NwOzOnjlmv069aN7Unfm7/TSruwX3gADr3POAN/9eZPFspLjePaWXGbfdhF1qtwyawW3v5BHweFy78RhjAkZljDOcKy8moffXkf/rknc883eXruuiDCx50RWHFhBcbkX/soXgasfh+Qe8PcfwLHCcxYf2yedBT+9lF+M78OSrYe4/PHF/H7BZsoqa1ofizEmJFjCOENSbAS/GN+Hx68bTGS4d789E3tORFEW7FrgnQtGJ8INr0JNBbxyA1Sdew1GfTfVJw9cxlUDu/D0J1/xzT8u4u95BdTWtd9zUYwx3mEJ4wwiwvXDu9Gnc6LXr92rQy8u6HQB83bO895F0/rAd2bCgXXwzp1Q1/zAducOMfzp+4N5885RdOkQyy/f+JKrn1rKZ9s92yDRGBNaLGH42YQeE/jy0JcUnCjw3kXPvxLG/Ro2vguLfut2tWHdk3n7rlH8+YYhHD9ZzY0zlnPzzOWs39vKmVzGmHbJEoafTew5EcB7g9/1Rt0LQ26GJb+HL552u5qIMPnCrnz88zH821V9Wbf3GFc/tZR7Xl7FV8Wl3o3RGNOmWcLws64JXRmcNph5u7zYLQXOIPi3noR+18CCf4X85z2qHhMZztRLerHkl2O595u9+efmIsY9vpifvbaGXYdsfypjjCWMgJjUaxLbjmxjbfFa7144LBymzIDeV8D798G6Nzy+RFJMJD+/sg9LfjmWqZf0Yu76/Vz++GJ++toath484d14jTFtiiWMAJh83mSSo5P565q/ev/iEVFw3QvQfRS8ORXyZrfoMqkJ0fzrpL4s+eVYfji6Bws2HODKPy3h9hfyyN99xMtBG2PaAksYARAfGc9tA27js32fsbpotfdvEBUHN74BOePgH/c74xrasmmz6YkxPHxVPz578Jv85Ju9+eKrEr7z18+Z8sxnzFu336bjGhNCRFv4i6QtyM3N1by8Vu4Q6yPl1eVMfGsiOck5zLhyhm9uUlsN794NX74Gw6fB+N8626S3QlllDa/nFTDzs50UHD5JZsdYbhzZjesv6kanePcOlTLGBC8RyVfV3MY+sxZGgMRFxvGjAT9i+f7lrDyw0jc3CY+Eb0+Hi++BFc/CnGugtHVbnsdHR3Dr6J4semAs028aSveUOH43fwsj/+dj7nt1NZ9/dYg6a3UY0y5ZCyOAKmoqmPjWRLondWf2+NmIiO9utvY1ZyA8tiN873noNsJrl9568AQvLtvN26v3cqKihu4pcUwZksW1QzLplhLntfsYY3zvXC0MSxgB9tKml3hsxWP87tLfnVqj4TMH1sNrN8HR3TD6PhjzEETGeO3yFdW1zFu/n9dWFrBsx2EAcrsnc83growf0Jn0RO/dyxjjG5Ywglh1XTW3zb+NbUe28drVr9GjQw/f3rDiGCx4GFa/AKl9nGNfsxr92WiVvUdP8u6avby9ai/bikoRgRE9OzGhf2fG9e9MZsdYr9/TGNN6ljCC3IGyA3zv/e+RFpfGy5NeJibCD3+Jb/8I3vsJHN8Lg66Hy/8dOrT8VMFz2XrwBB98uZ8P1u1ne5GzerxflySu6JvOmD5pXJjVkQgvb/RojGkZSxhtwNK9S7nzozuZkjOFR0c96p+bVhx3DmBa9ldnpfjIu5wB8vgUn91yR3EpH206yIcbD5K/+wh1Ch1iIxndO4WLz0tl1Hkp9EqN9+14jjGmSZYw2og/r/ozf1v3N3488MfcO+Re//3SPLIbPv41rH8DImJh6C0w6h7o2M2ntz1WXs3S7YdYtKWIz7YfYt+xCgDSE6O5qEcncnskk9u9Exd0SfT6VvPGmMZZwmgjautq+c2y3/DmtjeZ1HMSvxn9G6LC/bi2oWgzfP5nZ92G1jlbjAy5Gc6f4Kwg9yFVZc/hcj7bXsKKnSWs3HWEvUdPAhAdEcaAzA5cmNWRAZlJDMjsQK/UeOvGMsYHLGG0IarKzPUzeXLVkwxNH8rvx/ye9Lh0/wZxbC/kzYI1L8GJ/RCXAhdcDX2/BT0vhYhov4Sx9+hJVu0+wtqCo6wpOMr6fceoqHbO+4iOCCMnI4E+GUlc0DmR3ukJ9E5PILNjLGFh1p1lTEtZwmiD5u+cz8NLHyY8LJxb+9/Krf1vJS7Sz2saamvgq49h7auwbSFUlUJ0EvT4BvQc4ySPtAsgzD9/6dfU1rHjUBkb9h1jw97jbDl4gs0HTlB8ovJUmeiIMHqkxNMjNY4eqfFkJ8eR3SmO7ORYunaMJSYy3C+xGtNWWcJoowpPFPLEqidYsGsBabFpfL/P95nUaxLZidn+D6a6AnYuhi1zYcdiOLLTeT86CboOgcxh0HkgpPeDlPNavQWJJ46UVbG9uJTtRaV8VVTKrpIydpWUs6eknKra008gTE2IomvHWDKSYuicFEPnDjGkJUaTlhhNemI0qQnRJMdFERVh3V0mNAVNwhCRCcCTQDgwQ1UfO+NzcX0+CSgHblXVVe7UbUxbTxj11hSt4anVT7HiwAoALky7kEsyL2Fw+mAGpg70f8sD4Oge2LUUCldCYR4c3ABa63wWFgmdekKnXpDc0xk8T+rqTNtN7AzxaX7p1qqrUw6eqGBPSTkFR06y7+hJ9h87SeGRkxQdr+TA8QqOnaxutG5STASd4qPoGBflPMZGkhQbSce4SJJiIkmMiSAp1nlMiP76Kz46gtjIcOsWM21WUCQMEQkHtgLjgEJgJXCDqm5sUGYScC9OwhgBPKmqI9yp25j2kjDq7Svdx9ydc5m/cz5bjmwBIEzCyErIIjspm26J3egc35lOMZ3oFNOJDtEdSIhMID4yntiIWKLDo4kKjyJMfPDXc3UFHNoKRZugaCMc/goO73S+qhs5gCmmA8SlQlwniE12XkcnQUwSRCdCZDxExTs770bEOivSI2KdwfeIGAiPdlox4ZFOggqPcB7DIlxf4c5U4WacrKrlUGklRScqKT5RQUlZFSWlVZSUVnK4vJqj5VUcLqviaHk1x09Wc6Kyxq1vR1xUOHFR4cRGhRMb6XzFnPoKIyYynOiIMKIiwoiOCCcqIoyo8PrXYUSG13/JqecR4UJkuBARFkZEmBARHkZ4mBARJs5juBAuzvPwMCGskedhAmGu12GC6/Hr5yLYlOYQFywJ42LgP1V1vOv1rwBU9X8alPk/YJGqvuJ6vQW4DOjRXN3GtLeE0dCxymOsO7SOtcVr2XlsJ3uO76HgRAGl1c0fqxoRFkFkWCQREkF4WDhhEka4OI/1XwCCICKnPZ6p4S+Xxj4HoK4W6mqgrtr1WPv1o9a6HusQ1yN462dS6oM8/fXXwTde/qz6p7/b8H+ZhpHqGU+0wSd6ZuGzX/qdO/e31NE2xWsEr09b06K650oYEa0JykOZQEGD14U4rYjmymS6WRcAEZkGTAPo1s236wgCqUN0B76R+Q2+kfmN094vry7ncMVhSipKOFF1gtLqUsqqyjhZc5LK2koqayupqq2ipq6GGq2hpq6GOq2jTuuo1VpUFUWp0zoUPfX661+CZ/0adJ638g+Phteirg60pkFyqU8sdU5CUVdiOfWlXz9S58Ran3hUGzxyjuf1rxs8byTKsz/y9N/dRPkzLl3//dQGIdV/j05LWI39d2ksoekZ3+OGnzcS0lllm0iSTWm0jAffqsZiNe6LwTe7RfgzYTT2x8qZPxVNlXGnrvOm6rPAs+C0MDwJsD2Ii4wjLjKOrETfbPNhjAld/kwYhUDD6T1ZwD43y0S5UdcYY4wP+XPu4EogR0R6ikgUcD3w3hll3gNuEcdI4Jiq7nezrjHGGB/yWwtDVWtE5B5gAc7U2FmqukFE7nB9Ph2YizNDajvOtNrbzlXXX7EbY4yxhXvGGGMasDO9jTHGtJolDGOMMW6xhGGMMcYtljCMMca4pV0PeotIMbA70HE0kAocCnQQTQjW2Cwuz1hcngnWuCBwsXVX1bTGPmjXCSPYiEheU7MPAi1YY7O4PGNxeSZY44LgjM26pIwxxrjFEoYxxhi3WMLwr2cDHcA5BGtsFpdnLC7PBGtcEISx2RiGMcYYt1gLwxhjjFssYRhjjHGLJQw/EJFZIlIkIusDHUtDIpItIp+IyCYR2SAi9wU6JgARiRGRFSKy1hXXo4GOqSERCReR1SLyj0DH0pCI7BKRdSKyRkSCZtdNEekoIm+IyGbXz9rFQRBTH9f3qf7ruIjcH+i4AETkp66f+/Ui8oqI+Ob4vBawMQw/EJFLgVJgjqoOCHQ89USkC9BFVVeJSCKQD3xbVTcGOC4B4lW1VEQigaXAfaq6LJBx1RORnwG5QJKqXh3oeOqJyC4gV1WDaiGaiDwPfKqqM1zn2cSp6tEAh3WKiIQDe4ERqhrQhb4ikonz895PVU+KyN+Buar6XCDjqmctDD9Q1SXA4UDHcSZV3a+qq1zPTwCbcM5PDyh1lLpeRrq+guIvGxHJAq4CZgQ6lrZARJKAS4GZAKpaFUzJwuVy4KtAJ4sGIoBYEYkA4gii00UtYRgARKQHMARYHuBQgFPdPmuAIuBDVQ2KuIAngF8CdQGOozEKLBSRfBGZFuhgXHoBxcBsVzfeDBGJD3RQZ7geeCXQQQCo6l7gD8AeYD/OqaMLAxvV1yxhGEQkAXgTuF9Vjwc6HgBVrVXVwTjntw8XkYB35YnI1UCRquYHOpYmjFbVocBE4G5XV2igRQBDgb+q6hCgDHgosCF9zdVFNhl4PdCxAIhIMnAN0BPoCsSLyE2BjeprljBCnGuM4E3gJVV9K9DxnMnVfbEImBDYSAAYDUx2jRW8CnxTRF4MbEhfU9V9rsci4G1geGAjAqAQKGzQQnwDJ4EEi4nAKlU9GOhAXK4AdqpqsapWA28BowIc0ymWMEKYa3B5JrBJVR8PdDz1RCRNRDq6nsfi/E+0OaBBAar6K1XNUtUeON0Y/1TVoPjrT0TiXRMXcHX5XAkEfFaeqh4ACkSkj+uty4GATqo4ww0ESXeUyx5gpIjEuf7/vBxnbDEoWMLwAxF5BfgC6CMihSLyo0DH5DIauBnnL+X66YWTAh0U0AX4RES+BFbijGEE1RTWIJQBLBWRtcAK4ANVnR/gmOrdC7zk+u85GPhtYMNxiEgcMA7nr/ig4GqJvQGsAtbh/I4Omi1CbFqtMcYYt1gLwxhjjFssYRhjjHGLJQxjjDFusYRhjDHGLZYwjDHGuMUShjHGGLdYwjDGGOMWSxjG+JGIDBSR3SJyZ6BjMcZTljCM8SNVXYezrcgtgY7FGE9ZwjDG/4qA/oEOwhhPWcIwxv8eA6JFpHugAzHGE5YwjPEjEZkAxAMfYK0M08ZYwjDGT0QkBvgdcBfOTqQBPxTKGE9YwjDGf/4NmKOqu7CEYdogSxjG+IHrAKFxOGeCgyUM0wbZeRjGGGPcYi0MY4wxbrGEYYwxxi2WMIwxxrjFEoYxxhi3WMIwxhjjFksYxhhj3GIJwxhjjFv+H7NXzyJSqawfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "wait_time = 2\n",
    "\n",
    "plt.plot(x, gamma.pdf(x, a, scale = 1/wait_time)/np.sum(gamma.pdf(x, a, scale = 1/wait_time)), label = 'before')\n",
    "plt.plot(x, (x*wait_time*np.exp(-x*wait_time) * gamma.pdf(x, a, scale = 1/(wait_time)))/np.sum(x*wait_time*np.exp(-x*wait_time)*gamma.pdf(x, a, scale = 1/(wait_time))), label = '0.2')\n",
    "\n",
    "wait_time = 4\n",
    "\n",
    "plt.plot(x, (x*wait_time*np.exp(-x*wait_time) * gamma.pdf(x, a, scale = 1/(wait_time)))/np.sum(x*wait_time*np.exp(-x*wait_time)*gamma.pdf(x, a, scale = 1/(wait_time))), label = '6')\n",
    "\n",
    "plt.ylabel('$p(\\lambda | observed_time)$')\n",
    "plt.xlabel('$\\lambda$')\n",
    "\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01b6af59",
   "metadata": {},
   "source": [
    "Based on your plots, answer the following questions: \n",
    "- What happens to your estimate of the rate at which buses arrive as a function of how long you've waited for the first bus? What does this imply about your estimate of the time until the next bus will arrive?\n",
    "- What does the shape of these distributions tell you about the $\\Gamma$ distribution and choice of Poisson likelihood for this problem. Is this behavior generalizable? \n",
    "\n",
    "(challenge) If you've answered the above questions, now solve this problem analytically and explicitly compute the expected time until the next bus arrives using your posterior distribution for the rate at which buses arrive. "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}