probability_and_statistics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import Image, display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction \n",
    "Probability and statistics are two sides of the same coin. Probability theory gives us tools to model the world and predict the outcomes of random systems; statistics explains our observations of the world and evaluates the accuracy of our models. Oceanographers use tools from probability in fields as diverse as genomics and eddy modeling, while every experimentalist employs statistical methods in the analysis section of their papers. Probability and statistics are too broad to cover in any depth in a single 90-minute lecture, even as a review. But the basic principles of both subjects follow from a common goal: reasoning about information.\n",
    "\n",
    "## Learning Goals\n",
    "\n",
    "The learning goals for this lecture are to become familiar with some of the basic tools and ideas of probability theory and statistics: independence, Bayes' theorem, expectation, hypothesis testing, confidence intervals, and regressions. In addition, we will review some common statistical descriptors and probability distributions.\n",
    "\n",
    "# Probability Theory\n",
    "> \"Probability theory is the branch of mathematics that deals with probability\" -*en.wikipedia.org*\n",
    "\n",
    "Probability is the system by which we model random, or stochastic, phenomena- systems with a degree of uncertainty. We model such a system abstractly with a *probability space*: a set of possible outcomes for the system, each with a value (its probability) that denotes its likelihood of occurring. For this to truly be a probability space, we require that all probabilities sum to one (or integrate to one, in the case of continuous spaces):\n",
    "$$\\sum_ip_i = 1 \\qquad \\int_xp(x)dx = 1$$\n",
    "## Random Variables\n",
    "> \"You’ve gotta ask yourself a question: “Do I feel lucky?” Well, do ya, punk?\" -Harry Calahan, *Dirty Harry*\n",
    "\n",
    "Mathematically, a *random variable* associates a value with each outcome in a probability space. Random variables can be used to represent a wide range of systems. For example, imagine you are playing a game of roulette.\n",
    "\n",
    "<img src='https://www.888casino.com/blog/sites/newblog.888casino.com/files/inline-images/american-roulette-wheel_0.jpg' height=300 width=300/>\n",
    "\n",
    "The wheel has 18 black numbers, 18 red numbers, and 2 green numbers. We can model our earnings after betting \\\\$1000 on black with a random variable $X$. $X$ will be $2000 if the number is black, and zero otherwise. \n",
    "\n",
    "How much can we expect to win, if we make this same bet a large number of times? The *expectation* of a random variable lets us answer this and many similar questions:\n",
    "\n",
    "$$\\mathbb{E}\\left[X\\right] = \\sum_xxP(X=x)$$\n",
    "That is, the expectation is the weighted sum of all the values the random variable can take on, where the weight is the probability. For the roulette wheel,\n",
    "\n",
    "$$\\mathbb{E}\\left[X\\right] =  2000\\cdot\\frac{16}{38}+0 \\cdot \\frac{18}{38} \\approx \\$842.11$$\n",
    "\n",
    "Apparently, we shouldn't play roulette!\n",
    "\n",
    "One last note: when we *realize* a random variable, or observe its value, we call that *sampling*.\n",
    "\n",
    "## Independence\n",
    "> \"We're going to live on! We're going to survive! Today we celebrate our Independence Day!\" -President Whitmore, *Independence Day*\n",
    "\n",
    "Earlier, I mentioned that this subject deals with reasoning about information. Independence of random variables is best understood as a statement about mutual information:\n",
    "\n",
    "*If two random variables $X$ and $Y$ are independent, then learning something about the outcome of $Y$ doesn't give us any information about the outcome of $X$*\n",
    "\n",
    "We can write this formally using conditional probability notation. For outcomes $X = x$ and $Y = y$, the *conditional probability* $P(X=x | Y=y)$ denotes the probability of $X=x$ occurring, given that $Y=y$ has occured. For brevity, we'll write this as $P(X|Y)$, dropping the values when there's no ambiguity. So for independent random variables, $P(X|Y) = P(X)$ for every possible outcome. \n",
    "\n",
    "We can see how this works by first simulating a simple six-sided die roll and a coin flip. We'll let $D$ be a random variable representing the roll of the die, and $C$ the coin flip (where heads is 1 and tails is 0). Our probability space here has as its outcomes all the pairs $\\{(0, 0), (1, 0), \\dots, (6, 0), (0, 1), \\dots, (6, 1)\\}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Outcome (1, 0): 810 / 5096 ~ 0.15895\n",
      "Outcome (2, 0): 866 / 5096 ~ 0.16994\n",
      "Outcome (3, 0): 865 / 5096 ~ 0.16974\n",
      "Outcome (4, 0): 854 / 5096 ~ 0.16758\n",
      "Outcome (5, 0): 862 / 5096 ~ 0.16915\n",
      "Outcome (6, 0): 839 / 5096 ~ 0.16464\n",
      "Outcome (1, 1): 858 / 4904 ~ 0.17496\n",
      "Outcome (2, 1): 777 / 4904 ~ 0.15844\n",
      "Outcome (3, 1): 834 / 4904 ~ 0.17007\n",
      "Outcome (4, 1): 815 / 4904 ~ 0.16619\n",
      "Outcome (5, 1): 794 / 4904 ~ 0.16191\n",
      "Outcome (6, 1): 826 / 4904 ~ 0.16843\n",
      "Theoretical:    1   /    6 ~ 0.16667\n"
     ]
    }
   ],
   "source": [
    "# This function simulates rolling a die and flipping a coin\n",
    "def simulate_independent_random_variables():\n",
    "    \n",
    "    # Roll the die and store the result in D\n",
    "    D = np.random.randint(1, 6+1)\n",
    "    \n",
    "    # Flip the coin and store the result in C\n",
    "    C = np.random.randint(0, 1+1)\n",
    "    \n",
    "    # Return the pair (D, C)\n",
    "    return (D, C)\n",
    "\n",
    "# We will run 10,000 trials of this random system\n",
    "N = 10000\n",
    "\n",
    "simulations = [simulate_independent_random_variables() for n in range(N)]\n",
    "\n",
    "# To evaluate the independence of the two random variables, we'll take a look at\n",
    "# the fraction of the total set for each outcome \n",
    "for c in range(0, 1+1): # Python range is a weird function\n",
    "    for d in range(1, 6+1):\n",
    "        num_outcomes = len([x for x in simulations if x == (d, c)])\n",
    "        num_coin = len([x for x in simulations if x[1] == c])\n",
    "        print(f\"Outcome ({d}, {c}): {num_outcomes} / {num_coin} ~ {num_outcomes/num_coin:3.5f}\")\n",
    "print(f\"Theoretical:    {1}   /    {6} ~ {1/6:3.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about random variables that aren't independent? We can use two coin flips to help us demonstrate. Specifically, let $X$ be the random variable that takes on the sum of the two coin flips, and let $Y$ be the product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Outcome (0, 0): 2557 / 7502 ~ 0.34084\n",
      "Outcome (1, 0): 4945 / 7502 ~ 0.65916\n",
      "Outcome (2, 0): 0 / 7502 ~ 0.00000\n",
      "Outcome (0, 1): 0 / 2498 ~ 0.00000\n",
      "Outcome (1, 1): 0 / 2498 ~ 0.00000\n",
      "Outcome (2, 1): 2498 / 2498 ~ 1.00000\n"
     ]
    }
   ],
   "source": [
    "# This function simulates flipping two coins\n",
    "def simulate_nonindependent_random_variables():\n",
    "    \n",
    "    # Flip the first coin and store the result in D1\n",
    "    C1 = np.random.randint(0, 1+1)\n",
    "    \n",
    "    # Flip the second coin and store the result in D2\n",
    "    C2 = np.random.randint(0, 1+1)\n",
    "    \n",
    "    # Return the pair (C1+C2, C1*C2)\n",
    "    return (C1+C2, C1*C2)\n",
    "\n",
    "# We will run 10,000 trials of this random system\n",
    "N = 10000\n",
    "\n",
    "simulations = [simulate_nonindependent_random_variables() for n in range(N)]\n",
    "\n",
    "# To evaluate the nonindependence of the two random variables, we'll take a look at\n",
    "# the fraction of the total set for each outcome \n",
    "for y in range(0, 1+1): # Python range is a weird function\n",
    "    for x in range(0, 2+1):\n",
    "        num_outcomes = len([z for z in simulations if z == (x, y)])\n",
    "        num_y = len([z for z in simulations if z[1] == y])\n",
    "        print(f\"Outcome ({x}, {y}): {num_outcomes} / {num_y} ~ {num_outcomes/num_y:3.5f}\")\n",
    "#print(f\"Theoretical:    {1}   /    {6} ~ {1/6:3.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In theory, the probability that two coinflips sum to one is $1/2$, so clearly these events are nonindependent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayes Theorem\n",
    "\n",
    "> “Now, remember what Bayes' Theorem does: it helps us update a hypothesis based on new evidence.” - Dan Morris, *Bayes Theorem: A Visual Introduction For Beginners*\n",
    "\n",
    "Bayes Theorem gives us a way to relate conditional probabilities of events:\n",
    "\n",
    "$$P(A|B) = \\frac{P(A \\text{ and } B)}{P(B)}$$\n",
    "\n",
    "By rearranging the above expression, we can get an equivalent expression:\n",
    "$$P(A \\text{ and } B) = P(A|B)P(B)$$\n",
    "We chose the names for these events arbitrarily; therefore, we can swap them, and substitute into the original expression:\n",
    "$$P(A|B) = \\frac{P(B|A)P(A)}{P(B)}$$\n",
    "Finally, we can rewrite the denominator as a sum of two conditional probabilities (as $P(B) = P(B \\text{ and } A) + P(B \\text{ and not }A) = P(B \\text{ and } A) + P(B \\text{ and } A^C) = P(B|A)P(A)+P(B|A^C)P(A^C)$:\n",
    "$$P(A|B) = \\frac{P(B|A)P(A)}{P(B|A)P(A)+P(B|A^C)P(A^C)}$$\n",
    "\n",
    "This is an extremely powerful equation for reasoning about information, as the following example will demonstrate.\n",
    "\n",
    "Rapid point-of-care HIV testing allows for quick, cheap screening of a large number of people for HIV. When compared with viral load testing, the rapid tests have a $93.7\\%$ true positive rate and a $98.1\\%$ true negative rate. That is,\n",
    "$$P(Test +| HIV +) = .937$$\n",
    "$$P(Test -| HIV +) = .063$$\n",
    "$$P(Test +| HIV -) = .019$$\n",
    "$$P(Test -| HIV -) = .981$$\n",
    "There are approximately $1.1$ million people in the US living with HIV, out of about $327.2$ million people. What is the probability that a person who tests positive for HIV, has HIV- what is $P(HIV +|Test +)$? \n",
    "\n",
    "$$P(H+|T+) = \\frac{P(T+|H+)P(H+)}{P(T+|H+)P(H+)+P(T+|H-)P(H-)}$$\n",
    "\n",
    "$$P(H+|T+) = \\frac{.937 \\cdot 1.1 / 327.5}{.937 \\cdot 1.1 / 327.5+.019 \\cdot 326.4 / 327.5} \\approx 0.1425$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(H+ | T+) = 0.14775\n"
     ]
    }
   ],
   "source": [
    "def simulate_test():\n",
    "    status = np.random.choice(['h+', 'h-'], p = (1.1/327.2, 1 - 1.1/327.2))\n",
    "    if status == 'h+':\n",
    "        test = np.random.choice(['t+', 't-'], p = (.937, 1 - .937))\n",
    "    else:\n",
    "        test = np.random.choice(['t+', 't-'], p = (1 - .981, .981))\n",
    "    return (status, test)\n",
    "\n",
    "N = 50000\n",
    "\n",
    "simulations = [simulate_test() for x in range(N)]\n",
    "\n",
    "conditioned_simulations = [x[0] for x in simulations if x[1] == 't+']\n",
    "\n",
    "probability = len([x for x in conditioned_simulations if x == 'h+']) / len(conditioned_simulations)\n",
    "print(f\"P(H+ | T+) = {probability:3.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bayes theorem tells us that even with the high accuracy of the test, the low likelihood of a person actually having HIV means that most positive test results are in fact false positives. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linearity of Expecation\n",
    "> \"I find my life is a lot easier the lower I keep my expectations.\" -Bill Watterson, *Calvin and Hobbes*\n",
    "\n",
    "What is *expectation*? It's a generalization of the average (mean) of a set of numbers. What is the average of 3, 5, and 7? We can calculate this by dividing the sum of these numbers by three:\n",
    "$$\\text{Ave}\\left(3, 5, 7\\right) = \\frac{3 + 5 + 7}{3} = 5$$\n",
    "To extend this to expectation, first we can rewrite the fraction to give each number its own weight, $1/3$:\n",
    "$$\\text{Ave}\\left(3, 5, 7\\right) = \\frac{1}{3}3 + \\frac{1}{3}5 + \\frac{1}{3}7$$\n",
    "If we let $X$ be a random variable with a uniform distribution on $\\{3, 5, 7\\}$, then we can rewrite this average as:\n",
    "$$\\frac{1}{3}3 + \\frac{1}{3}5 + \\frac{1}{3}7 = \\sum_xP(X = x) \\cdot x$$\n",
    "\n",
    "So this average of three numbers is a weighted sum of the values taken on by the random variable $X$; each weight is the value of that number. This expression- $\\sum_xP(x)\\cdot x$- is known as the *expectation* of the random variable $X$. \n",
    "\n",
    "$$\\mathbb{E}\\left[X\\right] = \\sum_xxP(X = x)$$\n",
    "And with it written in this form, we can easily extend the idea of averaging to non-uniformly distributed random variables (even infinite ones!)\n",
    "\n",
    "Let $Y$ be a random variable with a geometric distribution on the natural numbers $\\{0, 1, 2, \\dots\\}$:\n",
    "$$P(Y = k) = p(1-p)^k$$\n",
    "This distribution models, for example, the number of flips of a weighted coin before the first heads. What is the expectation of $Y$?\n",
    "\n",
    "$$\\mathbb{E}\\left[Y\\right] = \\sum_kkP(Y = k) = \\sum_{k=0}^{\\infty}kp(1-p)^k = p(1-p)\\sum_{k=0}^{\\infty}k(1-p)^{k-1} = -p(1-p)\\frac{d}{dp}\\sum_{k=0}^{\\infty}(1-p)^k$$\n",
    "$$\\mathbb{E}\\left[Y\\right] = -p(1-p)\\frac{d}{dp}\\left[\\frac{1}{1-(1-p)}\\right] = \\frac{1-p}{p}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEYCAYAAAATRII7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcXHWZ7/HP0129pLdsXVlJJ3RHg8CQiG1kURYFDFwFZVxgRBbRiMo4Djpuc6961ZnrLI6jgyOi7I6o4BaV1YWACEIHAoSwSDYSsnX2Tu/Lc/84p5tKUVVd3amqU939fb9e9eqz1annnDpV3z5L/Y65OyIiIlEriboAERERUCCJiEiRUCCJiEhRUCCJiEhRUCCJiEhRUCCJiEhRUCCJiEhRUCCJiEhRUCCNkpk9bWanFfg1F5nZ42bWZmYfL+Rrj5SZ3WhmX426jpEayftqZhvN7IyRjovCcMtVbPUWo1x95pM/G7n8LkmcV67f00J85xUkkMIV02lmBxMeVxfitTPUk/UblWp6dz/G3e/LeXGZfRq4z91r3f1bBX7tCSGi93WImU0N/+E4NWn4LWb2MzOz0cw3cbkUPq+UzTrJ17aRzXyzfc9yWWPyaxbis1HIPaS3u3tNwuPKAr72eDEfeDrXMzWzWK7nOdYUyzpw973AtcDfDw4zs/8DHA1c5Grrq+CKZdvIZCzUmBV3z/sD2AickWJ4E7AHOD7snwPsAk5LeN7ngLXAXuAGoDLh+XOAnwKtwAbg4wnj5gE/C8ftBq4Oh98CDACdwEGCvY7PAuuAtvC13pkwn1dMn7xMwGuA+4B9BIFxbtKyfwp4EtgP/DhxGVKsk5TzAn4P9ANdYR2vTvHctMuR5j35TFhXNxAbZj1kXA7gtcBj4XN/DPwI+Opwy5Uw738I590OXAfMBO4M5/dbYGqa5fgscHvSsG8C3xpunaRZB4nva8b1SYbtk6RtngzbaoplmkuwvTUB7wY2A3NSTHcZ8KuE/heAnyT0bwaWJNZD5u05q+10uPWSNO3xwOPhtLeF8/1qNusl0zYz0u1mmNdJ9Z2QcdvI9B2TYh0M99kYmm/4mi+F0z4HvGWY9yzT9ruRzN+fDixM6L9xsK4M6+SMbN4fRvi9N/S84SbIxSN5QZLGfQh4BqgC7gb+Pel5a8I3fhrwYMIKKwFWAV8AyoFGYD3wVqAUeAL4BlANVAJvTFcPwYd+TjjP9xJs3LMz1c/LH/Aygi+Cz4d1vDncmBYlTPdIOP9p4bJekWZdDDev+4APZljPGZcjRf2rw3U7abjnZ1qOsNZNBP/VlwHvAnoT3qts1tHDBF8mc4GdBB/g1wIVBGH8xTTLMR/oAOrC/lJgG3BClsuUvA6G3usst4t022fifNJuqxney+uBuwi+7I5PM00jwZdBCTA7fA9eShi3FyhJUc9Qd9KyZLudZrWdJWwXfxduA+cDPWT3Gc64zYxku8lm/Sevkyy2jYzfMSnWQcrPRtJ3ySIS/vkAFgBNw7xnmWrcSJrtMxyfNpAyrJPBeefse++QZRpuglw8wuIOEnx4Bh8fShi/AniKIE0rkp53RUL/OcC6sPsNwItJr/M5gv8CTiT4IMcy1JMyIMPxq4HzMk2fsBG9CdhO+MEPx90KfClhuosSxv0rcE2a1x1uXveRIZCGW44U9X8g2+dnWg7gFGArYAnj/8TLXzzZrKP3JYz7KfCdhP6/BX6Roc4/AheH3WcObiNZLtMHksan3TbSbBfpts+h+WTaVjPUeSzBF8Z7hnmPNhPshVxAcKjvEeAogr2nFamWK8P2nNV2mu12Fm4XLyVtF39M2C4yfYYzbjMj2W6yWf/J62S4bYNhvmOS1kHaz0bifIGFBKF6BlA23HaZRY1pt8+w/3ACKWffe4mPQh53fIe7/zbNuO8RhNJyd+9OGrc5oXsTQeJC8J/xHDPblzC+FHiA4D+CTe7el01hZnYxcBXBfyQANUB9Ns8N69ns7gNJdc5N6N+e0N2RsAyjmVdao1iOxHWbzfPTLcccgv/MPaluEsYPt1w7Ero7U/TXZFiOHwIXAjcDfxP2Z7tMh6yDRFmuz3TbZ6JM22o65QSHYX6WYRqAlcBpBF9mKwn+2TuV4Atz5TDPTZbVdjqC7SzVdpG4vjKtl2w/C9lsN6NZ/8m1Jsv2O2a4z8YQd3/BzD4BfAk4xszuBq5y962jrDF5fLrtczRy+b03JPLLvs2sBvhPguO/XzKzaUmTzEvobiD4bwOCFb3B3ackPGrd/ZxwXEOGE31DG4eZzScIxCuB6e4+hWA311JNn8JWYJ6ZJa7LBoL/DEdq1PPKcjmSjXQ9pLMNmJt0BVhDQncu11EqtwGnmdkRwDsJA+lw3tsRrI9022eiTNtqOouBNVl84Q0G0pvC7pUEgXQq6QMp0/ac0Qi3k1TbReL6yrRecrnNZLP+U62TTOtpuO+YQcN9Ng59QfcfuvsbCULUgX8Zppbh3stM22cHwamSQbNGMO+8fKYjDySCE9Cr3P2DwG+Aa5LGf8zMjgiD6vMEJ8cgODRxwMw+Y2aTzKzUzI41s9eH47YBXzOzajOrNLOTE+a5g+A4MgTHf51g9xszu4zgcAlppk/2Z4Jj6J82s7LwOv23E5y4HKnDmVc2y5Gv5z8E9AEfN7OYmZ0PLE0Yn8t19Aru3kpwOPMGgi+eZ8JRh7NM2T433faZKNO2ms4SgkNhw1kJnE5wDmELwX/9y4DpBBcTpJJpex7OSNbpQwQX4lwZbhfnceh2kWm95HKbyWb9j3SdDPcdM2i4z8YQC35n+GYzqyC4eKmTYP2Npr5BmbbP1cDfhOtjGcE/MYkK9b03pJCB9Kuk3yH9PNxAlwFXhNNcBRxvZu9LeN4PgXsITkKuB74K4O79BCtgCcFVM7uA7wOTE8YtBF4EthCcfB30/4D/He7CnwN8nWDD2QH8FcHJP1JNb2afShzh7j3AucDZYQ3/TXA+49mRrqDDmZe7r81iOfLy/LDu84FLCU6kv5eEQ025XEcZ/JDg2PvQ4brDXKZsn5ty+0yaV9ptNUMJi8kikNz9eYLzsw+E/QfCOh4MXzeVtNtzFq+X9TpN2C4uJziUeBHwa4JDkcN9hnP5ucpm/Y9onWTxHTM4XcbPRpIK4GthfduBGQQhMuL6EmTaPv8uXIZ9wPsIzrclKsj3XiI79NBmcTGzjQQn8dOdexKRMcTM/kxwcvuGqGuR4lMMh+xEZJwys1PNbFZ4uOoS4DiCy9lFXmF8/LpXRIrVIuAnBFe7rQPe5e7boi1JilVRH7ITEZGJQ4fsRESkKIyrQ3b19fW+YMGCqMsQERkzVq1atcvd41HXAeMskBYsWEBLS0vUZYiIjBlmlrLliCjokJ2IiBQFBZKIiBQFBZKIiBQFBZKIiBQFBZKIiBQFBZKIiBQFBZKIiBSFCR9Iff0DfPsPL7Dy+daoSxERmdAmfCCVlhjX3r+eu5/ePvzEIiKSNxM+kMyMxng161sPRl2KiMiENuEDCaApXsO61vaoyxARmdAUSEBjvJrWtm7aunqjLkVEZMJSIBHsIQGs116SiEhkFEhAU7wagHU6jyQiEhkFEtAwrZrSEtMekohIhBRIQHmshIZpVdpDEhGJkAIp1FhfrT0kEZEIKZBCTTNq2LC7nf4Bj7oUEZEJSYEUaqyvpqdvgJf2dkZdiojIhKRACjXNCC79XrdL55FERKKgQAo11oeXfu9UIImIRCGWrxmb2fXA24Cd7n5sOOzHwKJwkinAPndfkuK5G4E2oB/oc/fmfNU5aFp1OZMnlbF+ly5sEBGJQt4CCbgRuBq4eXCAu793sNvMvg7sz/D80919V96qS2JmNMWrtYckIhKRvB2yc/f7gT2pxpmZAe8Bbs3X649GY7xGe0giIhGJ6hzSm4Ad7v6XNOMduMfMVpnZ8kwzMrPlZtZiZi2trYd3k72meA2tbd0cUCOrIiIFF1UgXUjmvaOT3f144GzgY2Z2SroJ3f1ad2929+Z4PH5YRTWGbdrpB7IiIoVX8EAysxhwPvDjdNO4+9bw707g58DSQtT2cqvfOo8kIlJoUewhnQE86+5bUo00s2ozqx3sBs4C1hSisIZpVZSWmNq0ExGJQN4CycxuBR4CFpnZFjO7PBx1AUmH68xsjpndEfbOBP5oZk8AjwC/cfe78lVnosFGVnXITkSk8PJ22be7X5hm+KUphm0Fzgm71wOL81XXcJriamRVRCQKaqkhSWNcjayKiERBgZSkKa5GVkVEoqBAStIYXmmnCxtERApLgZRkqJFVBZKISEEpkJJMqy5nSlUZ63Rhg4hIQSmQkphZeDtz7SGJiBSSAimFJjWyKiJScAqkFBrVyKqISMEpkFJoUiOrIiIFp0BKYejSb92sT0SkYBRIKQw2srp+lwJJRKRQFEgplMdKmK9GVkVECkqBlEZjvFo/jhURKSAFUhpN8Ro27upQI6siIgWiQEqjMV5NT/8AW/Z2RF2KiMiEoEBKo3HoduY6jyQiUggKpDSa1Oq3iEhB5fMW5teb2U4zW5Mw7Etm9pKZrQ4f56R57jIze87MXjCzz+arxkzUyKqISGHlcw/pRmBZiuHfcPcl4eOO5JFmVgp8GzgbOBq40MyOzmOdaTXFa9TIqohIgeQtkNz9fmDPKJ66FHjB3de7ew/wI+C8nBaXpcb6au0hiYgUSBTnkK40syfDQ3pTU4yfC2xO6N8SDkvJzJabWYuZtbS2tua00KYZNew62M3+TjWyKiKSb4UOpO8ATcASYBvw9RTTWIphaX8M5O7XunuzuzfH4/HcVBkavHusDtuJiORfQQPJ3Xe4e7+7DwDfIzg8l2wLMC+h/whgayHqS6ZLv0VECqeggWRmsxN63wmsSTHZo8CrzOxIMysHLgBWFKK+ZPOnVxFTI6siIgURy9eMzexW4DSg3sy2AF8ETjOzJQSH4DYCHw6nnQN8393Pcfc+M7sSuBsoBa5396fzVWcmZaUlNEyrYt1O7SGJiORb3gLJ3S9MMfi6NNNuBc5J6L8DeMUl4VFojNdoD0lEpADUUsMwmuLVamRVRKQAFEjDUCOrIiKFoUAahtq0ExEpDAXSMHTpt4hIYSiQhjGtupypamRVRCTvFEhZaIzX6JCdiEieKZCy0BSv1iE7EZE8UyBloTGuRlZFRPJNgZQFNbIqIpJ/CqQsNM3QlXYiIvmmQMpCw7SgkVVd2CAikj8KpCyUlZbQML1Ke0giInmkQMpSY70u/RYRyScFUpaa4tVs2t1BX/9A1KWIiIxLCqQsNcVrwkZWO6MuRURkXFIgZakxHl76rXsjiYjkhQIpS01qZFVEJK/yFkhmdr2Z7TSzNQnD/s3MnjWzJ83s52Y2Jc1zN5rZU2a22sxa8lXjSEwdamRVe0giIvmQzz2kG4FlScPuBY519+OA54HPZXj+6e6+xN2b81TfiDXFa9Tqt4hInuQtkNz9fmBP0rB73L0v7H0YOCJfr58PjfFqNR8kIpInUZ5D+gBwZ5pxDtxjZqvMbHmmmZjZcjNrMbOW1tbWnBeZKGhktUeNrIqI5EEkgWRm/wj0Af+TZpKT3f144GzgY2Z2Srp5ufu17t7s7s3xeDwP1b7s5QsbtJckIpJrBQ8kM7sEeBvwPnf3VNO4+9bw707g58DSwlWY3uCl3zqPJCKSewUNJDNbBnwGONfdO9JMU21mtYPdwFnAmlTTFtpgI6vaQxIRyb18XvZ9K/AQsMjMtpjZ5cDVQC1wb3hJ9zXhtHPM7I7wqTOBP5rZE8AjwG/c/a581TkSg42s6tJvEZHci+Vrxu5+YYrB16WZditwTti9Hlicr7oOV2N9jX4cKyKSB2qpYYSaZlSzcXe7GlkVEckxBdIINdXX0NvvamRVRCTHFEgj1DRDjayKiOSDAmmEGuuD3yKt26nzSCIiuaRAGqGp1eVMqy7XHpKISI4pkEahsb5ae0giIjmmQBqFxni19pBERHJMgTQKTYONrHaokVURkVxRII1CY9jI6jrtJYmI5IwCaRSawkZW1WKDiEjuKJBGYV7YyKratBMRyR0F0igMNrKqVr9FRHJHgTRKTfEa3RdJRCSHFEij1BivZpMaWRURyRkF0ig1xdXIqohILimQRqlp6HbmOo8kIpILCqRRGmxkVZd+i4jkRl4DycyuN7OdZrYmYdg0M7vXzP4S/p2a5rmXhNP8xcwuyWedozHYyKr2kEREcmPYW5ib2flZzKfL3e9IMfxG4Grg5oRhnwV+5+5fM7PPhv2fSXrNacAXgWbAgVVmtsLd92ZRS8E01ldrD0lEJEeGDSTge8AvAcswzSnAKwLJ3e83swVJg88DTgu7bwLuIymQgLcC97r7HgAzuxdYBtyaRb0F0xSv4XfP7oi6DBGRcSGbQLrT3T+QaQIz+8EIXnOmu28DcPdtZjYjxTRzgc0J/VvCYaleezmwHKChoWEEZRy+xng1P24JGlmdXFVW0NcWERlvsjmH9NHhJnD3i3JQS6JUe2Oe5rWvdfdmd2+Ox+M5LiOzJjWyKiKSM9kE0uNmdkEOX3OHmc0GCP/uTDHNFmBeQv8RwNYc1pATjYOXfu9UIImIHK5sAunNwHvDK+IW5uA1VwCDV81dQnB+KtndwFlmNjW8Cu+scFhRGWxkdf0uXdggInK4hj2H5O6bgHea2TLgQTN7FBhIGH9uuuea2a0EFzDUm9kWgivnvgb8xMwuB14E3h1O2wxc4e4fdPc9ZvYV4NFwVl8evMChmJSVljB/epX2kEREciCbixows0XAp4EHgG+TEEiZuPuFaUa9JcW0LcAHE/qvB67P5nWi1Biv0R6SiEgOZPM7pK8B5wKfdPc781/S2NIUr+G+53bS1z9ArFQNX4iIjFY2e0j9wPHu3pXvYsaixng1vf3O5r2dHFlfHXU5IiJjVjb/0v90uDAys+NzVM+YM3jpt27WJyJyeLIJpBvCq92mpXsA1+W70GKlVr9FRHIjm0N2k4FVZG46qDU35Yw9U6rKmTtlEn94tpXlpzRFXY6IyJiVzWXfCwpQx5h2yUnz+ec7nuXJLfs47ogpUZcjIjImDXvIzsw+ndD97qRx/5yPosaaC5c2UFsZ47sr10ddiojImJXNOaTEZoM+lzRuWQ5rGbNqK8u46IT53LlmGxv1myQRkVHJJpAsTXeq/gnrspMWECsp4ft/1F6SiMhoZBNInqY7Vf+ENaOukvOPn8ttLVvYdbA76nJERMacbAJpsZkdMLM24Liwe7D/r/Jc35jyoVMa6ekf4KY/bYy6FBGRMWfYQHL3Unevc/dad4+F3YP9uitdgqZ4DWcdPZObH9pEe3df1OWIiIwpanwtxz58ahP7O3v58aObh59YRESGKJBy7PiGqSxdMI3r/riB3v6sGkUXEREUSHnx4VMbeWlfJ79+suhucisiUrQUSHlw+qIZvHpmDd9duR53XYgoIpINBVIelJQYy09p4tntbax8fsI28yciMiIFDyQzW2RmqxMeB8zsE0nTnGZm+xOm+UKh6zxc5y6ew6y6SjUnJCKSpaxuYZ5L7v4csATAzEqBl4Cfp5j0AXd/WyFry6XyWAmXv/FI/umOZ3hi8z4Wz1OjqyIimUR9yO4twDp33xRxHXlx4RvCRlfvXxd1KSIiRS/qQLoAuDXNuBPN7Akzu9PMjkk3AzNbbmYtZtbS2lpc52tqKmK8/4T53LlmuxpdFREZRmSBZGblwLnAbSlGPwbMd/fFwH8Bv0g3H3e/1t2b3b05Ho/np9jDcOnJCygrKeF7D+hckohIJlHuIZ0NPObuO5JHuPsBdz8Ydt8BlJlZfaELzIUZtZX89evmctuqLbS2qdFVEZF0ogykC0lzuM7MZpmZhd1LCercXcDacupDb2qkV42uiohkFEkgmVkVcCbws4RhV5jZFWHvu4A1ZvYE8C3gAh/DvzBtjNfw1qNncfNDG9XoqohIGpEEkrt3uPt0d9+fMOwad78m7L7a3Y9x98XufoK7/ymKOnPpw6c2cqCrjx+p0VURkZSivspuwnhtw1SWHjmN6x5Yr0ZXRURSUCAV0EdObWLr/i5+9YQaXRURSaZAKqDTFsVZNLNWja6KiKSgQCogM2P5KY08t6ON+9ToqojIIRRIBfb2xXOYPbmS765Uc0IiIokUSAU22Ojqw+v3sHrzvqjLEREpGgqkCFywtIG6ypj2kkREEiiQIlBTEeP9J87nrqe3s0GNroqIAAqkyFx60pGUlarRVRGRQQqkiMRrK3jX647gdjW6KiICKJAiNdjo6o1/2hB1KSIikVMgRejI+mqWHTOLWx7axEE1uioiE5wCKWIfPrWJA119XPeA9pJEZGJTIEVsybwpvO242Xzjt8/zi8dfirocEZHIxKIuQODf372YPe09fPK2J6ibFOPNR82MuiQRkYLTHlIRqCwr5dqLmzlmTh0f+cFjPLJhT9QliYgUnAKpSNRUxLjh0tczd+okLr/pUdZuPRB1SSIiBRVZIJnZRjN7ysxWm1lLivFmZt8ysxfM7EkzOz6KOgtpek0FP7j8DdRWxLj4+kfYqFYcRGQCiXoP6XR3X+LuzSnGnQ28KnwsB75T0MoiMmfKJG6+/A0MuHPRdX9mx4GuqEsSESmIqAMpk/OAmz3wMDDFzGZHXVQhLJxRw42XvZ697T1cfN0j7OvoibokEZG8izKQHLjHzFaZ2fIU4+cCmxP6t4TDDmFmy82sxcxaWlvHz03vjjtiCt+7uJkNu9r5wI2P0tGjH86KyPgWZSCd7O7HExya+5iZnZI03lI85xX3/Xb3a9292d2b4/F4PuqMzEkL6/nWha9l9eZ9XPGDx+jpG4i6JBGRvIkskNx9a/h3J/BzYGnSJFuAeQn9RwBbC1Nd8Vh27Cy+dv5x3P98K1f9ZDX9A6/IZBGRcSGSQDKzajOrHewGzgLWJE22Arg4vNruBGC/u28rcKlF4T2vn8fnzzmKXz+5jS/8cg3uCiURGX+iaqlhJvBzMxus4YfufpeZXQHg7tcAdwDnAC8AHcBlEdVaFJaf0sSe9l6uWbmOadXlfPKsRVGXJCKSU5EEkruvBxanGH5NQrcDHytkXcXuM8sWsa+jh//6/QtMqSrn8jceGXVJIiI5o7bsxhAz45/e+Vfs7+zlK79ey5RJZfz1646IuiwRkZwo5t8hSQqlJcZ/XrCEkxdO59M/fZLfrt0RdUkiIjmhQBqDKmKlfPf9zRw7p46P/fAx/rx+d9QliYgcNgXSGFVTEeOGy5Yyb1oVH7yphcdf3Bt1SSIih0WBNIZNqy7nlsuXMrmqjHdd8xD/etezdPX2R12WiMioKJDGuNmTJ/Gbv30T5792Lv993zrO+dYDtGzU/ZREZOxRII0Dk6vK+Ld3L+bmDyylp2+Ad3/3Ib74yzUc7Fb7dyIydiiQxpFTXh3n7k+cwiUnLuDmhzfx1m/cz8rnx0+DsyIyvimQxpnqihhfOvcYbr/iRCrLSrjk+ke46ierdQsLESl6CqRx6nXzp/Gbj7+JK09fyIrVWznjP1Zyx1MTsilAERkjFEjjWGVZKZ966yJ+eeXJzJpcyUf/5zGuuGUVO3UXWhEpQgqkCeCYOZP5xUdP5rNnH8UfntvJGf+xkp+0bFar4SJSVBRIE0SstIQrTm3izr97E0fNquPTtz/Jxdc/wuY9HVGXJiICKJAmnMZ4DT9afgJfecexPLZpL2/9z/u54cENuvGfiEROgTQBlZQY7z9hPvdcdSpLj5zG//3VWs76xkpueXgTHT367ZKIRMPG03mE5uZmb2lpibqMMcXd+c1T2/juyvU89dJ+6ipjXLC0gYtPnM8RU6uiLk9E8szMVrl7c9R1gAJJQu7Oqk17ueHBjdz19HbcnbOOnsVlJy9g6ZHTCO/uKyLjTDEFUsFv0Gdm84CbgVnAAHCtu38zaZrTgF8CG8JBP3P3LxeyzonGzGheMI3mBdPYuq+TWx7exK2PvMhdT2/n6Nl1XHryAs5dPIfKstKoSxWRcarge0hmNhuY7e6PmVktsAp4h7uvTZjmNOBT7v62kcxbe0i51dnTzy9Wv8QND27g+R0HmV5dzt+8oYGLTpjPzLrKqMsTkRyY0HtI7r4N2BZ2t5nZM8BcYG3GJ0rBTSov5cKlDVzw+nn8ad1ubnhwA1f/4QW+c986/tdxs7ns5CNZMm9K1GWKyDgR6TkkM1sA3A8c6+4HEoafBvwU2AJsJdhbejrNPJYDywEaGhpet2nTpvwWPcFt2t3OTX/axG0tm2nr7uO1DVO49KQFnH3sbMpjumhTZKwppj2kyALJzGqAlcA/ufvPksbVAQPuftDMzgG+6e6vGm6eOmRXOAe7+7i9ZTM3PbSJDbvaqa2McfqiGZx59ExOWxSntrIs6hJFJAsTPpDMrAz4NXC3u/9HFtNvBJrdfVem6RRIhTcw4Nz/l1bueGobv3tmJ7vbeygrNU5squfMo2dy5mtmMmuyzjeJFKsJHUgWXD98E7DH3T+RZppZwA53dzNbCtwOzPdhilUgRat/wHn8xb3cs3YH967dwYZd7QAcd8RkznzNTM48ZiaLZtbqEnKRIjLRA+mNwAPAUwSXfQN8HmgAcPdrzOxK4CNAH9AJXOXufxpu3gqk4uHurGs9yD1rd3DP0ztYvXkfAA3TqoI9p6Nn0jx/KrFSnXcSidKEDqR8UiAVr50HuvjtMzu5d+12HnxhNz39A0ytKuP0o2Zw1tEzOXlhvc47iURAgZQnCqSx4WB3H/c/38q9a3fw+2d3sr+zlxKDo2bV0bxgavAD3flTmTNlUtSliox7CqQ8USCNPb39Azy6cQ8Pr9/Dqk17ePzFfXT09AMwZ3Jl2HrEVJrnT2PRrFpKS3T+SSSXiimQCv7DWJFEZaUlnNRUz0lN9QD09Q/wzLY2WjbtoWXjXv68YTcrntgKQG1FjCUNU3h9uAe1pGEKVeXahEXGC+0hSVFzd7bs7RwKqJaNe3l+ZxvuUFpiHDOnjtfNn8qSeVM4alYdjfFqynShhEjWimkPSYEkY87+zl4ee3EvLRuDkHpiyz66eoMLNstLS1g4o4ajZtfymll1wd/ZddTXVERctUhxUiDliQJpYurpG2D9roPuM6IyAAAKvklEQVQ8u62NZ7Yf4JltbTy77QA727qHpqmvqeA1s2s5alYtR82q4zWz62iaUU1FTK2Xy8RWTIGkA/Ay5pXHSjhqVh1HzarjHcwdGr77YDfPbW/jme1tPLPtAM9uP8BND22ipy/Ym4qVGE3xYG/qVTNqWFBfzYLp1Syor6amQh8NkULTp07Grek1FZy0sIKTFtYPDevrH2Dj7vZgL2r7AZ7d1sajG/bwy9VbD3luvLaCI6dXs6C+igX11WF3EFiTyrVXJZIPCiSZUGKlJSycUcvCGbW8ffGcoeEdPX1s3NXBxt3tbNjVzsZd7Wzc3c7vn21l18Eth8xjZl0FC6ZXc2T9yyE1f3oVc6ZMoq4ypqaRREZJgSQCVJXHOHpOHUfPqXvFuLauXjbt7hgKqg27g7/3rN3BnvaeQ6atqYgxZ0olc6dMYk74mDtlEnOnBt0zayvUXJJIGgokkWHUVpZx7NzJHDt38ivG7e/sZeOudrbs7WTrvk5e2hf83bq/k9Wb97G3o/eQ6UsMZtVVDgXU4GN2XSUz6iqYUVtJfU25QksmJAWSyGGYPKmMxfOmsDjNnXM7evrYuq8rCKkwsAZD6/EX93HHU9vo7T/0SlczmF5dwYzaijCkgqCaWVdBvLZyaFi8tkJXCcq4okASyaOq8hgLZ9SwcEZNyvEDA07rwW627+9iZ1s3O9u62HGgm9a2LnYe6GZnWzfPbDtAa1s3Ayl+oTG1qizYq6otZ3p1BdOqy5leXc60mqB/ek0506rLqa+uoG6Szm9JcVMgiUSopMSYWVfJzLrMNzHsH3B2t3ez80A3rWFw7TzQzY7w766D3Ty5dx+7D/bQ1t2Xch6xEmNqGFhBUFUE3dXlTK0uZ0pVGVOrypk8qWyou6q8VCEmBaNAEhkDSkuMGbWVzKgd/u673X397G3vZdfBbva097CnvSepu4c97d08NUyAAZSVGpMnDYZV2SHdUxLCa8qkcuomxaitLKOuMvhbHtN5MBkZBZLIOFMRK2XW5NKsbx3f3dfP/o5e9nb0sq+jh32d4d+O3rD75f6X9nXy9Nb97OvopbO3P+N8J5WVHhJSdZPKqKssSxgWdNdVllFbGaO2MkZNRRk1lTFqymNUV5Tq4o4JRoEkMsFVxEqZUVfKjGEOGybr6u3nQGcQZPs7e2nr6uVAVy8HOvs40Bl0t3X1DQ3b097Dpt0dHOgMpu9LdVIsyaSyUmoqY9RWxIKgqohRXXFo/+DfqvIYNRWlVIVhFvTHqCovpboiRkWsRIcfi1wkgWRmy4BvAqXA9939a0njK4CbgdcBu4H3uvvGQtcpIulVlpVSWTbyIIOgFfeu3oEwtHrZ39lHe3cfB7v7ONjVR1v4t72nj7auweG9HOzuY/OeDg52B9O3dfVlFWwQXHJfXR4EWlVFKdXlQVjVVMSYVF5KVXkQYpPKS6kqK2VSeenQ8EllsXB8OLwsYdryUrUwnyMFDyQzKwW+DZwJbAEeNbMV7r42YbLLgb3uvtDMLgD+BXhvoWsVkfwws6Ev/OEu6MjE3enuG+Bgdx+dPf0c7O6jo6eP9u5+Onr6OBj+fbm/j47uftp7+ujo6ae9u48dbV10dPfT0RNM09U7QE//wIjqKCu1oYCeFD4qy0qC/vLB/oTx5SVUxoJxlUPjSqiIvfy8ylgpFWXBdIPjKspKxvWeXhR7SEuBF9x9PYCZ/Qg4D0gMpPOAL4XdtwNXm5n5eGqaXEQOm9nLQZBLvf0DdPb209nzclC93N1PZ28QaJ1D/f10hY/Onn66esPn9/azt72Hrb2D0wzQFU6f7Z7dK5cZKmKHhteM2gpuu+KknK6DKEQRSHOBzQn9W4A3pJvG3fvMbD8wHdiVPDMzWw4sB2hoaMhHvSIywZSVllBWWkJdZVneXqO3fyAIsN5+unoG6O4LA6uvn+7eYFxXOGxoXG8/3b39dPeF48PpJ+U4kKMSRSCl2tdM/lchm2mCge7XAtdCcD+kwytNRKQwBkOvNo+hN9ZEcSZuCzAvof8IYGu6acwsBkwG9hSkOhERiUQUgfQo8CozO9LMyoELgBVJ06wALgm73wX8XuePRETGt4IfsgvPCV0J3E1w2ff17v60mX0ZaHH3FcB1wC1m9gLBntEFha5TREQKK5LfIbn7HcAdScO+kNDdBby70HWJiEh09GsuEREpCgokEREpCgokEREpCgokEREpCjaerqY2s1ZgU9R1RKSeFC1ZTCBafi2/ln905rt7PJfFjNa4CqSJzMxa3L056jqiouXX8mv5x/7y65CdiIgUBQWSiIgUBQXS+HFt1AVETMs/sWn5xwGdQxIRkaKgPSQRESkKCiQRESkKCqQxxsyWmdlzZvaCmX02xfirzGytmT1pZr8zs/lR1Jkvwy1/wnTvMjM3szF/KWyibJbfzN4TbgNPm9kPC11jPmWx/TeY2R/M7PHwM3BOFHXmg5ldb2Y7zWxNmvFmZt8K182TZnZ8oWs8bO6uxxh5ENyuYx3QCJQDTwBHJ01zOlAVdn8E+HHUdRdy+cPpaoH7gYeB5qjrLvD7/yrgcWBq2D8j6roLvPzXAh8Ju48GNkZddw6X/xTgeGBNmvHnAHcS3HH7BODPUdc80of2kMaWpcAL7r7e3XuAHwHnJU7g7n9w946w92GCO/KOF8Muf+grwL8CXYUsrgCyWf4PAd92970A7r6zwDXmUzbL70Bd2D2ZV96Nesxy9/vJfOfs84CbPfAwMMXMZhemutxQII0tc4HNCf1bwmHpXE7wH9N4Mezym9lrgXnu/utCFlYg2bz/rwZebWYPmtnDZrasYNXlXzbL/yXgIjPbQnDPtb8tTGlFYaTfD0Unkhv0yahZimEpr9s3s4uAZuDUvFZUWBmX38xKgG8AlxaqoALL5v2PERy2O41g7/gBMzvW3fflubZCyGb5LwRudPevm9mJBHeePtbdB/JfXuSy/n4oVtpDGlu2APMS+o8gxSEJMzsD+EfgXHfvLlBthTDc8tcCxwL3mdlGguPoK8bRhQ3ZvP9bgF+6e6+7bwCeIwio8SCb5b8c+AmAuz8EVBI0PDoRZPX9UMwUSGPLo8CrzOxIMysHLgBWJE4QHrL6LkEYjafzBzDM8rv7fnevd/cF7r6A4Bzaue7eEk25OTfs+w/8guDCFsysnuAQ3vqCVpk/2Sz/i8BbAMzsNQSB1FrQKqOzArg4vNruBGC/u2+LuqiR0CG7McTd+8zsSuBugiuOrnf3p83sy0CLu68A/g2oAW4zM4AX3f3cyIrOoSyXf9zKcvnvBs4ys7VAP/AP7r47uqpzJ8vl/yTwPTP7e4LDVZd6eAnaWGdmtxIciq0Pz5F9ESgDcPdrCM6ZnQO8AHQAl0VT6eip6SARESkKOmQnIiJFQYEkIiJFQYEkIiJFQYEkIiJFQYEkIiJFQYEkIiJFQYEkIiJFQYEkkkdmtsDMnjWzm8J71NxuZlVR1yVSjBRIIvm3CLjW3Y8DDgAfjbgekaKkQBLJv83u/mDY/QPgjVEWI1KsFEgi+ZfcPpfa6xJJQYEkkn8N4b15ILhfzx+jLEakWCmQRPLvGeASM3sSmAZ8J+J6RIqSbj8hkn8D7n5F1EWIFDvtIYmISFHQ/ZBERKQoaA9JRESKggJJRESKggJJRESKggJJRESKggJJRESKwv8HJYgwMDWdVGYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "delta_x = 0.05\n",
    "x = np.array([delta_x * (x + 1) for x in range(int(1 / delta_x))])\n",
    "y = (1 - x) / x\n",
    "\n",
    "plt.plot(x, y)\n",
    "plt.title(\"Expectation of a random variable $Y$ with a geometric distribution\")\n",
    "plt.xlabel(\"p\")\n",
    "plt.ylabel(\"E[Y]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expectation has an extremely useful property- *linearity* . \n",
    "\n",
    "$$\\text{Linearity of Expectation: For any scalars } a \\text{ and } b\\text{, and random variables } X \\text{ and } Y \\text{, } \\mathbb{E}\\left[aX + bY\\right] = a\\mathbb{E}\\left[X\\right] + b\\mathbb{E}\\left[Y\\right]$$\n",
    "\n",
    "Why is this a useful property? For one, linearity allows us to \"pass the expectation through\" an expression:\n",
    "\n",
    "$$\\mathbb{E}\\left[\\sum_{k=1}^{n}a_kX^k\\right] = \\sum_ka_k\\mathbb{E}\\left[X^k\\right]$$\n",
    "\n",
    "This gives us the ability to break down complicated functions of random variables into pieces than are much easier to calculate. The canonical expression for the variance of a random variable is derived using linearity of expectation:\n",
    "\n",
    "$$\\mathbb{E}\\left[(X-\\mu)^2\\right] = \\mathbb{E}\\left[X^2-2\\mu X + \\mu^2\\right] = \\mathbb{E}\\left[X^2\\right]-2\\mu\\mathbb{E}\\left[X\\right]+\\mu^2 = \\mathbb{E}\\left[X^2\\right] -\\mu^2$$\n",
    "where $\\mu = \\mathbb{E}\\left[X\\right]$.\n",
    "\n",
    "In the case of continuous random variables, we alter the expression for expectation somewhat:\n",
    "$$\\mathbb{E}\\left[X\\right] = \\int_xxP(X=x)dx$$\n",
    "\n",
    "You can think of $P(X=x)dx$ as the probability of $X$ taking on a value between $x$ and $x + dx$- a kind of weighted Reimann sum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to Describe a Distribution\n",
    "> \"I just wish that Statistics was as easy as arranging numbers in chronological order, finding the median, lower and upper quartiles, and placing them on a Box and whisker's chart\" -Charmaine J. Forde\n",
    "\n",
    "We've seen that random variables have probability distributions associated with them. We can, of course, describe a finite, discrete probability distribution (which is a function) by listing all the input-output pairs. But this doesn't give us much useful information for understanding what might happen if we sample from it; and we can't list all the pairs for an infinite distribution.\n",
    "\n",
    "Instead, probabilists and statisticians have developed a wide range of quantitative descriptors for the behavior of a probability distribution. Some of the more common ones are:\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\text{Mean:} &\\qquad \\mu(X) = \\mathbb{E}\\left[X\\right]\\\\\n",
    "\\text{Variance:} &\\qquad \\sigma^2(X) = \\mathbb{E}\\left[(X-\\mu(X))^2\\right] \\\\\n",
    "\\text{Standard Deviation:} &\\qquad \\sigma(X) = \\sqrt{\\mathbb{E}\\left[(X-\\mu(X))^2\\right]} \\\\\n",
    "\\text{Median:} &\\qquad P(x \\leq \\text{Median}(X)) = 0.5 \\\\\n",
    "\\text{Mode:} &\\qquad \\arg\\max_xP(X=x) \\\\\n",
    "\\text{Support:} &\\qquad \\{x | P(X = x) > 0\\} \\\\\n",
    "\\text{Range:} &\\qquad \\max(\\text{Support}(X)) - \\min(\\text{Support}(X)) \\\\\n",
    "\\text{Quartiles: } &\\qquad P(x \\leq Q_1) = 0.25; P(X \\leq Q_3) = 0.75 \\\\\n",
    "\\text{Interquartile Range:} &\\qquad Q_3 - Q_1 \\\\\n",
    "\\text{Skewness:} &\\qquad \\mathbb{E}\\left[\\left(\\frac{X - \\mu(X)}{\\sigma(X)}\\right)^3\\right] \\\\\n",
    "\\text{Kurtosis:} &\\qquad \\mathbb{E}\\left[\\left(\\frac{X - \\mu(X)}{\\sigma(X)}\\right)^4\\right]\n",
    "\\end{align*}$$\n",
    "\n",
    "You may have heard of some of these before; others (like kurtosis) are rare enough that you may never encounter them. I want to focus on three groups of these descripters for a moment:\n",
    "\n",
    "<table style=\"width:100%\">\n",
    "  <tr>\n",
    "    <th>Type</th>\n",
    "    <th>Center</th>\n",
    "    <th>Spread</th> \n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>Expectation</td>\n",
    "    <td>Mean</td>\n",
    "    <td>Variance</td> \n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>Robust</td>\n",
    "    <td>Median</td>\n",
    "    <td>IQR</td> \n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>Extreme</td>\n",
    "    <td>Mode</td>\n",
    "    <td>Range</td> \n",
    "  </tr>\n",
    "</table>\n",
    "\n",
    "Each of these descriptors represents either the central value or the spread of a distribution. But the properties of these descriptors are often quite different. \n",
    "\n",
    "The two \"Expectation\" descriptors, the mean and the variance, are derived from the expectation of (functions of) the random variable. They both fall into an entire class of mathematical descriptors known as moments (skewness and kurtosis are also types of moments). Moments of a distribution have many nice analytical properties, and expressions involving moments arise frequently in physics and engineering. The mean and variance are, however, not always defined for a distribution; in addition, moments are not particularly resistant to outliers. That means that extreme values of a distribution can have an influence on moments out of proportion to their probabilites.\n",
    "\n",
    "Consider a random variable $X$ that takes on the value $0$ with probability $0.5 * (1 - 10^{-5})$, the value $1$ with the same probability, and the value $10^{10}$ with probability $10^{-5}$. What are the mean and variance of $X$?\n",
    "\n",
    "$$\\mathbb{E}\\left[X\\right] = 0.5 * (1 - 10^{-5}) + 10^{10} * 10^{-5} \\approx 10^{5}$$\n",
    "$$\\mathbb{E}\\left[X^2\\right]-\\mu^2 \\approx 0.5 * (1 - 10^{-5}) + 10^{20} * 10^{-5} - 10^{10}\\approx 10^{15}$$\n",
    "\n",
    "The two \"Robust\" descriptors address this problem, sacrificing the analytical properties of the descriptors for resistance to outliers. You could, for example, multiply the larget value in the support of a random variable by $10^{100}$, and the median and IQR would remain unchanged. These types of descriptors are most useful when discussing samples; a box-and-whiskers plot, for example, will show the median and IQR of a dataset as the lines on the large central box:\n",
    "<img src='https://upload.wikimedia.org/wikipedia/commons/2/2f/BoxPlot.jpg' height=300 width=800/>\n",
    "Finally, the two \"Extreme\" descriptors represent extremes of the random variable. The *mode* is the value with the highest likelihood of being sampled, while the *range* is the difference between the largest and smallest possible values to sample. Modes often arise in discussions of maximum likelihood; range gives you guarantees about where your samples will lie. \n",
    "\n",
    "Some probability distributions are so common, they have names. Here are a couple important ones:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Binomial Distribution\n",
    "\n",
    "The binomial distribution represents the number of heads you get after flipping $n$ coins, where each coin has a probability of $p$ of getting heads.\n",
    "\n",
    "$$P(k; n, p) = \\binom{n}{k}p^K(1-p)^{n-k}$$\n",
    "$$\\begin{align*}\n",
    "\\text{Mean:} &\\qquad np\\\\\n",
    "\\text{Variance:} &\\qquad np(1-p) \\\\\n",
    "\\text{Median:} &\\qquad \\lfloor np \\rfloor \\text{ or } \\lceil np \\rceil \\\\\n",
    "\\text{Support:} &\\qquad \\{0, 1, \\dots, n\\}\n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8lfXd//HXJ5uQEAhJGAlJ2HuKTAdqi1iRZb0dVaij3t5VW++un9U6qtW7u96tttW7anFPloqKdSNDwt6yAiQEkjADIfvz++Nc0ENMyElykuuMz/PxOI9zzjXf1wl8rv29RFUxxhgTHiLcDmCMMab1WNE3xpgwYkXfGGPCiBV9Y4wJI1b0jTEmjFjRN8aYMGJFPwCIyN9F5L5Anq+IfCIitzRjXueLyNamjl/H9N4VkVnO5++KyGI/Tvs7IrLIX9OrNe1YEdkkIp1bYvrm7ERkiIgscTuHm6zotwIRyRWRkyJyXEQOi8g7ItLtVH9VvU1VH27tXP6ar4g8KCKVIlLivL4SkcdFpIvXvD5X1b4+TusFH7Jfpqqz/ZA9W0RURKK8pv2iqk5s7rTrcSvwmarub6Hpf42IxIjIG86/QxWRCbX6i4j8RkQOOq/fioi0Vj5fiMgwEVkpIqXO+7CzDPuCiBSIyDHn3+LpjRVVXQccEZErWiV4ALKi33quUNUEoAtwAPiLy3n87VVVTQSSgelAZ2Cld+H3B6dABfO/2/8EnndhvouB64G6Vja3AtOAocAQYDKenAFBRGKA+cALQAdgNjDf6V6X/wGyVbUdMAX4lYic49X/RQJo+VpbMP/nCUqqWga8AQw41U1E/ikiv3I+TxCRPBH5sYgUOlssN3oNmyQiz4lIkYjsFpFfnCqCzmGOL0TkTyJyRER2isg4p/teZ3qz6plvBxF525nuYedzRhOWr1JVNwJXA0XAj72Xy2ve/09E8p09g60icomITALuAa529orWOsN+IiKPiMgXQCnQo47DTSIifxGRoyKyRUQu8eqRKyLf8PruvTfxmfN+xJnn2NqHi5zfcIUz7RUiMs6r3yci8rDzu5eIyCIRSanrtxGRTKAnsNyr2z9F5Aln769ERJaLSM/G/OYNUdUKVX1MVRcD1XUMMgv4g6rmqWo+8Afgu75M2/kt3xCRV538q0RkqP/SAzABiAIeU9VyVf0zIMDFdQ2sqhtVtfzUV+fl/Zt+AlwiIrF+zhkUrOi3MhGJx1MQl51lsM5AEpAO3Aw8ISIdnH5/cfr1AC4EZgI3eo07GlgHdAReAl4BzgV64dnSe1xEEuqYZwTwLJAFZAIngccbv4QeqlqNZ+vs/Nr9RKQvcAdwrrN3cCmQq6rvAY/i2WtIUFXv4nEDni3SRGB3HbMcDewEUoAHgDkikuxD1Auc9/bOPJfWypoMvAP8Gc9v+kfgHRHp6DXYdXj+BmlADPCTeuY1GNipqlW1ul8L/BLPVux24JH6wjor8/ped/uwvHUZCKz1+r7W6earqcDrePbyXgLmiUh0PfnXnSX/X8+Sb52e2WbMurNlFJG/ikgpsAUoABae6ues2CqBBg83hiIr+q1nnogcAY4B3wR+d5ZhK4GHnK3mhcBxoK+IROJZYfxcVUtUNRfPVtkNXuPuUtVnnaL7KtDNmVa5qi4CKvCsAM6gqgdV9U1VLVXVEjyF58JmLvM+PIWgtmogFhggItGqmquqOxqY1j+dLbgqVa2so38hni3BSlV9FdgKXN6s9B6XA9tU9Xln3i/jKSTex4SfVdWvVPUk8BpQ3/Hm9kBJHd3nqOqXzsrgxbOMj6q2P8vr101aQkgAjnp9PwokiPh8XH+lqr7h/F3+CMQBY+rJP+Qs+b/vY75TGRPrC+RMKxHPRsccoLzWICV4/h5hx4p+65mmqu3xFLs7gE+l/is4DtbaGizF8w8/Bc+WpPeW7m48ewSnHPD6fBJAVWt3+9qWvojEi8iTziGjY3gOe7R3VjRNlQ4cqt1RVbcDdwEPAoUi8oqIdG1gWnsb6J9fa0twN9DQNH3Rla/vWdT+zb2Pk5/6W9XlMHUXKl/HbynHgXZe39sBx2v9nmdz+m+jqjVAHv757U+pnQ/ne10r0NNUtdo5pJUB/Fet3onAEb8lDCJW9FuZ8w9xDp6t3fMaOXoxnr2ALK9umUC+H6L9GM/u7mjnBNipwx5NuorDOc9wBfB5Xf1V9SVVPQ/Psijwm1O96plkQwUovdaWaSaePQ2AE0C8Vz/vlW1D093Hmb/3qWk35Tdfh+d8RFSDQ9bDOe9Q3+ueJk52I56TuKcMdbr56vSVaM7fPYN///ZnEJGNZ8n/97PkG1Lr7zukERmj8Dqm72xgxODZGww7VvRbmXhMxXP8dnNjxnUO2bwGPCIiiSKSBfwIz1UNzZWIZy/giHMc+4GmTEREokWkP/AynuL6xzqG6SsiFzsn0sqc+Z46wXgAyJbGX6GTBvzAmf9VQH/+fRx3DXCN028k8G2v8YqAGjznSOqyEOgjIteJSJSIXI3nJPzbjcyHquYB24BRjR3XaxoJZ3k9Wt944rk/IM75GiMicV5F9DngRyKS7hTEHwP/9Bo3V0S+e5ZY54jIDGdldheeQyl1nrNS1YFnyX9bPdP/BM+/jx84y3GH0/2jOpYzTUSuEZEEEYkUkUvxnDPxHnYC8JHXyd6wYkW/9bwlIsfxHNN/BJjlXOXSWHfi2XLdiecyvJeAZ/yQ7zGgDZ69iWXAe40c/2pn+Y4AC4CDwDmqWtcWXyzwa2de+/EU7FNbqa877wdFZFUj5r8c6O1M8xHg26p60Ol3H54tvcN4Tpi+dGokVS11hv/COZl4xrFoZxqT8RTCg8DPgMmqWtyIbN6e5MxzMK1lK56VazrwvvP51B7Mk8BbwHpgA54T10/C6cslO3L2Cw/m4znXdBjPss2o57xLk6hqBZ5LSmfi+fd1E57DpRVOxntE5N1Tg+M5lJPn5Pk9cJeqzvea5HeA+vYqQp74ftjOGNNczt7NauASVS1wO09DROQ84HZVvbae/g8CvVT1+lYN1kQiMhh4SlXHup3FLVb0jTFNFmxF39jhHWOMCSu2pW+MMWHEtvSNMSaMNPl64ZaSkpKi2dnZbscwxpigsnLlymJVTW1ouIAr+tnZ2eTk5LgdwxhjgoqI1NUm1dfY4R1jjAkjVvSNMSaMWNE3xpgwYkXfGGPCiBV9Y4wJI1b0jTEmjFjRN8aYMGJF3zTJZ18VsXL3YbdjGGMaKeBuzjKBrbpG+f2irfztkx20i4vio59MICUh1u1Yxhgf2Za+8dmxskpumb2Cv32ygyuGduVkZTWPvtOoh38ZY1xmRd/4ZHvhcaY9/gWfbyvmV9MG8Zdrh/OfF/Rkzup8luxo6kOkjDGtzYq+adBHWw4w/YkvOHqykpe+N4brx3iesnfHxb3oltyGX8zbQHlVdQNTMcYEAiv6pl6qyl8/2c7Ns3PI7BjPgjvPY1T35NP946IjeWjKIHYWneD/PtvpYlJjjK+s6Js6nayo5s6XV/Pb97YyeUhX3rhtHOnt23xtuIv6pfGtwZ35y0fb2XOw1IWkxpjGsKJvvibvcClX/m0J76wv4P9N6sefrxlGm5jIeoe/f/JAoiKE+xdswJ7EZkxgs6JvzrB850GmPP4Few+X8sysc/mvCT0RkbOO0zkpjh9N7MsnW4t4b8P+VkpqjGkKK/oG8By/f37Zbr7zj+W0j49m3u3juahfms/jzxqbxYAu7XjwrY0cL69qwaTGmOawom+oqKrhnrkbuG/eBs7vncK828fTMzWhUdOIiozgkemDKCwp54+LvmqhpMaY5vKp6IvIJBHZKiLbReTuOvr/SEQ2icg6EflQRLK8+lWLyBrntcCf4U3zFZWUc93/LePlL/fw/Qk9+cesc2kXF92kaQ3P7MB1ozL555JdbMg/6uekxhh/aLDoi0gk8ARwGTAAuFZEBtQabDUwUlWHAG8Av/Xqd1JVhzmvKX7KbfxgXd4Rpjy+mA37jvKXa4fzs0n9iIw4+/H7hvzs0n4kt43h3nkbqK6xk7rGBBpftvRHAdtVdaeqVgCvAFO9B1DVj1X11PV6y4AM/8Y0/jZvdT5X/X0pESK8cds4rhja1S/TTYqP5t7L+7N27xFe/nKPX6ZpjPEfX4p+OrDX63ue060+NwPven2PE5EcEVkmItPqGkFEbnWGySkqKvIhkmmq6hrlfxZu5q5X1zA0oz3z7xjPoPQkv85j2rB0xvboyG/f20JRSblfp22MaR5fin5d+/t17reLyPXASOB3Xp0zVXUkcB3wmIj0/NrEVJ9S1ZGqOjI1NdWHSKYpjpZWctM/V/DkZzu5fkwmL9wyukVayBQRfjV9EGWVNTy60BpkMyaQ+FL084BuXt8zgH21BxKRbwD3AlNU9fTmnaruc953Ap8Aw5uR1zTR9sISpj6xmCU7inl0+mB+NW0wMVEtd/FWz9QE/vPCHsxdnc+S7dYgmzGBwpf/9SuA3iLSXURigGuAM67CEZHhwJN4Cn6hV/cOIhLrfE4BxgOb/BXe+OZfmw4w7YklHC+v4qXvjeG60ZmtMt/bL+pFZnI8v5hvDbIZEygaLPqqWgXcAbwPbAZeU9WNIvKQiJy6Gud3QALweq1LM/sDOSKyFvgY+LWqWtFvJarK4x9t43vP55CdEs+CO87j3Ozkhkf0k7joSB6aOpCdRSd46lNrkM2YQCCB1lbKyJEjNScnx+0YQa+0ooqfvr6Od9YXMHVYV35z5RDioutvP6cl3f7iKj7YfIAP/vsCsjq2dSWDMaFORFY650/Pyu7IDUF7D5Uy469LWLihgJ9f1o/Hrh7mWsEHuG/yAGIiI7h//kZrkM0Yl1nRDzFLdxxkyuOLyT9ykme/ey7/eWHDDaa1tM5Jcfx4Yh8+/aqIheutQTZj3GRFP0SoKrOX5HL908tJbhvD/NvHM6Gv7w2mtbQbxmQxsGs7fvnWRkrKKt2OY0zYsqIfAsqrqrn7zfU8sGAjE/qkMu/28fRoZINpLc3TINtgio6X88cPrEE2Y9xiRT/IFZaUce1Ty3g1Zy93XNSL/5s5ksQmNpjW0oZ1a8/1o7OYvSTXGmQzxiVW9INYUUk5U/7yBZsLSnj8uuH85NK+RDSzwbSW9pNL+5LcNpZ75663BtmMcYEV/SD25qo89h8r48XvjWbyEP80mNbSktpEc9/k/qzNO8pL1iCbMa3Oin6QUlXmrMpjRGZ7RmR2cDtOo0wZ2pXxvTwNshWWlLkdx5iwYkU/SG3cd4yvDhxn+ojga8VaRHh46iDKK2t45B1rkM2Y1mRFP0jNWZVPdKRwxZAubkdpkh6pCdw2oSfz1+zjC2uQzZhWY0U/CFVV17BgbT4X90ujfXyM23Ga7PsTepLVMZ5fzNtAWaU1yGZMa7CiH4Q+31ZM8fEKZgThoR1vcdGRPDx1ELuKT/CkNchmTKuwoh+E5qzOp318NBcF0B23TXVBn1QmD+nCE59sJ7f4hNtxjAl5VvSDzLGyShZt3M8VQ7q26ENQWtN9kwcQGxnBffM3WINsxrSw0KgaYeS99fspr6phxoizPaY4uHRq52mQ7fNtxby9rsDtOMaENCv6QebNVXl0T2nLsG7t3Y7iVzeMzWZwehIPv72JY9YgmzEtxop+ENl7qJTluw4xY3i6680l+1tkhPDI9EGeBtkWWYNsxrQUK/pBZP6afACmDQ+dQzvehmS0Z+aYLJ5bmsv6PGuQzZiWYEU/SHiaXchnVPdkuiXHux2nxfz40r50TIjl3nnWIJsxLcGKfpBYm3eUncUnmBGiW/mntIuL5r7JA1iXd5QXl+92O44xIceKfpCYsyqPmKgIvhWkzS40xhVDunB+7xR+995Wa5DNGD+zoh8EKqpqeGvtPiYO6ES7AH1Aij+JCA9NHUR5dQ2/etsaZDPGn6zoB4FPthZyuLQypK7Nb0j3lLZ8f0JPFqzdx+fbityOY0zIsKIfBOauziclIYbze6e6HaVV3XZhT7qntOX++RutQTZj/MSKfoA7UlrBh5sLuWJoV6Ijw+vP5d0g298/3eF2HGNCQnhVkSD09roCKqpruDLIW9RsqvN6pzBlaFf++vEOdlmDbMY0mxX9ADd3dT690xIY2LWd21Fc84vJ/YmNjuB+a5DNmGazoh/AcotPsHL3YWaMyAi5ZhcaIy0xjp9e2pfPtxXzljXIZkyz+FT0RWSSiGwVke0icncd/X8kIptEZJ2IfCgiWV79ZonINuc1y5/hQ93c1fmIwLThXd2O4rrvjM5iSIY1yGZMczVY9EUkEngCuAwYAFwrIgNqDbYaGKmqQ4A3gN864yYDDwCjgVHAAyLSwX/xQ5eqMnd1PuN6dqRLUhu347guMkJ4ZNpgDh4v5w/vb3U7jjFBy5ct/VHAdlXdqaoVwCvAVO8BVPVjVS11vi4DTp11vBT4QFUPqeph4ANgkn+ih7aVuw+z51ApM4aH5wncugzOSGLm2GyeW7ab7YUlbscxJij5UvTTgb1e3/OcbvW5GXi3MeOKyK0ikiMiOUVFdiMOwJur8mkTHcmkQZ3djhJQ7ry4F9GREcxeYu3yGNMUvhT9us4g1nkJhYhcD4wEfteYcVX1KVUdqaojU1PD6wakupRVVvPOun1MGtSZtrFRbscJKB0TYrliSFfeXJVnx/aNaQJfin4e0M3rewawr/ZAIvIN4F5giqqWN2Zcc6aPthRyrKyK6SHeomZTzRqXRWlFNXNW5rkdxZig40vRXwH0FpHuIhIDXAMs8B5ARIYDT+Ip+IVevd4HJopIB+cE7kSnmzmLOavySUuMZXyvFLejBKQhGe0Z1q09zy3dTY21uW9MozRY9FW1CrgDT7HeDLymqhtF5CERmeIM9jsgAXhdRNaIyAJn3EPAw3hWHCuAh5xuph4Hj5fzydZCpg1PJzIifK/Nb8iscVnsLD7B4u3FbkcxJqj4dMBYVRcCC2t1u9/r8zfOMu4zwDNNDRhu3lq7j6oaDasWNZviW4O78Mg7m3luaS4X9LHzQMb4yu7IDTBzV+czoEs7+nUO32YXfBEbFcm1ozL5cEshew+VNjyCMQawoh9QthceZ23eUdvK99F1ozOJEOGFZXb5pjG+sqIfQOauziNCYMowa3bBF12S2nDpwE68smIvJyusvX1jfGFFP0DU1CjzVu/j/N6ppCXGuR0naMwcm83Rk5UsWJvvdhRjgoIV/QCxfNch8o+ctEM7jTS6ezL9Oicye8lua3bZGB9Y0Q8Qc1blkRAbxcQB1uxCY4gIM8dms6ngGCt3H3Y7jjEBz4p+ADhZUc27G/Zz2aDOtImJdDtO0Jk2vCuJcVHMXmondI1piBX9ALBo036Ol1cx3Q7tNEl8TBT/MbIb764voPBYmdtxjAloVvQDwJxV+XRNimNM945uRwlaN4zJolqVF5fvcTuKMQHNir7LCkvK+HxbEdNHpBNhzS40WXZKWyb0SeWlL/dQUVXjdhxjApYVfZctWLOPGoXp9rCUZps5LpuiknLe27jf7SjGBCwr+i6bsyqfoRlJ9EpLcDtK0LuwdypZHeN5bkmu21GMCVhW9F20Zf8xNhUcs3bz/SQiQrhhTBY5uw+zIf+o23GMCUhW9F00d1U+URHCFUOt2QV/uWpkN9pER/Lc0ly3oxgTkKzou6S6Rpm3Jp8JfVPpmBDrdpyQkdQmmukj0pm/Zh+HT1S4HceYgGNF3yVLdhRz4Fg5M0bYCVx/mzk2i/KqGl7L2et2FGMCjhV9l8xZlU+7uCgu7pfmdpSQ069zO0Z3T+b5ZbuptscpGnMGK/ouOFFexXsb9nP5kK7ERVuzCy1h1rhs8g6f5OMthQ0PbEwYsaLvgvc27OdkZTVXWrMLLWbigE50SYpj9tJct6MYE1Cs6Ltgzuo8MpPjOSerg9tRQlZUZATfGZ3J59uK2V543O04xgQMK/qtrODoSZbsOMj04emIWLMLLemaUZnEREbY4xSN8WJFv5XNW70PVeyGrFaQkhDL5UO68MbKPI6XV7kdx5iAYEW/Fakqc1blcU5WB7JT2rodJyzMHJvF8fIq5q7KczuKMQHBin4r2rjvGNsKj9tWfisantmBoRlJzF5qj1M0Bqzot6o5q/KJiYxg8pAubkcJKzPHZrO98DhLdhx0O4oxrrOi30oqq2tYsDafS/qn0T4+xu04YeXyIV1IbhvDbGt90xgr+q3l821FFB+vsEM7LoiLjuSac7vxr80HyDtc6nYcY1xlRb+VzFmVT4f4aCb0tWYX3PCdMVkA9jhFE/Z8KvoiMklEtorIdhG5u47+F4jIKhGpEpFv1+pXLSJrnNcCfwUPJsfKKlm06QBXDO1KTJStZ92Q3r4NEwd05pUv91BWWe12HGNc02AFEpFI4AngMmAAcK2IDKg12B7gu8BLdUzipKoOc15Tmpk3KL27voCKqhprUdNlM8dlcbi0krfW7nM7ijGu8WWzcxSwXVV3qmoF8Aow1XsAVc1V1XWAPZG6Dm+uyqdHSluGZiS5HSWsje3Rkd5pCcxemmuXb5qw5UvRTwe8GybPc7r5Kk5EckRkmYhMq2sAEbnVGSanqKioEZMOfHsPlfLlrkPMGGHNLrhNRJg5LpsN+cdYvfeI23GMcYUvRb+uStWYzaRMVR0JXAc8JiI9vzYx1adUdaSqjkxNTW3EpAPfvNX5AEyzq3YCwozh6STGRtnD003Y8qXo5wHdvL5nAD4fFFXVfc77TuATYHgj8gU1VWXu6nxGd08mo0O823EM0DY2im+PzOCd9QUUlpS5HceYVudL0V8B9BaR7iISA1wD+HQVjoh0EJFY53MKMB7Y1NSwwWbN3iPsLD7BlXYCN6DcMCaLymrllS/tcYom/DRY9FW1CrgDeB/YDLymqhtF5CERmQIgIueKSB5wFfCkiGx0Ru8P5IjIWuBj4NeqGjZFf+7qfGKjIrhscGe3oxgvPVITuKBPKi8u301ltV17YMJLlC8DqepCYGGtbvd7fV6B57BP7fGWAIObmTEoVVTVsGDtPiYO7ExiXLTbcUwts8ZmcfPsHBZtPMDl1haSCSN2p1AL+XhrIUdKK5lhJ3AD0oS+aXRLbmPt8ZiwY0W/hcxdlU9KQgzn905xO4qpQ2SEMHNMNl/mHmLTvmNuxzGm1VjRbwFHSiv4cMsBpgxNJyrSfuJAddXIDOKiI3h+Wa7bUYxpNVaRWsDb6wqorFZmjLBDO4GsfXwM04alM3d1PkdLK92OY0yrsKLfAuasyqNvp0QGdm3ndhTTgBvGZlFWWcPrK+3yTRMerOj7WW7xCVbtOcJ0a3YhKAzsmsS52R14buluqmusPR4T+qzo+9mc1fmIwLRhdmgnWMwal82eQ6V8+lWh21GMaXFW9P3I0+xCHuN7ptA5Kc7tOMZHlw7sTFpiLLOX7HY7ijEtzoq+H+XsPszeQyftBG6QiY6M4Dujs/j0qyJ2FZ9wO44xLcqKvh/NWZVHm+hILh1ozS4Em2tHdyM6Unh+qW3tm9BmRd9PyiqreXtdAZcN6kzbWJ9atzABJC0xjssGdeH1nL2cKK9yO44xLcaKvp98uLmQkrIqeyRiEJs1LpuS8irmOs9AMCYUWdH3k7mr8+jULpaxPTu6HcU00YjM9gxKb8dzS+1xiiZ0WdH3g4PHy/lkaxHThqcTGWHX5gcrEWHm2Gy+OnCcZTsPuR3HmBZhRd8P3lq7j6oaZcZwO7QT7KYM7Ur7+GieW5rrdhRjWoQV/WZSVV5cvodB6e3o2znR7TimmeKiI7n63G4s2nSAfUdOuh3HGL+zot9Mn20rZlvhcW4a393tKMZPrh+d5azM7fJNE3qs6DfT04t3kZYYy+QhXd2OYvykW3I8l/TvxMtf7qWsstrtOMb4lRX9ZvjqQAmffVXEzLFZxETZTxlKZo3N5tCJChauL3A7ijF+ZZWqGZ5ZvIu46AiuG53ldhTjZ+N7daRHaltm2x26JsRY0W+ig8fLmbM6nxkjMkhuG+N2HONnIsKssdms3XuENXuPuB3HGL+xot9ELyzbQ0VVjZ3ADWFXnpNBQmwUz9nD000IsaLfBOVV1Ty/bDcX9U2lV1qC23FMC0mIjeLKEem8va6A4uPlbscxxi+s6DfBgjX7KD5ezs3n9XA7imlhN4zNpqK6hldX2OMUTWiwot9IqsrTi3fRr3Mi43tZOzuhrldaAuf1SuGFZbupqq5xO44xzWZFv5GW7DjIlv0l3HRed3sGbpiYOTaLgqNlfLDpgNtRjGk2K/qN9PTiXaQkxDBlqN2MFS4u6d+J9PZtmL001+0oxjSbFf1G2F54nI+2FHL9mCzioiPdjmNaSWSEcMPYLJbtPMTW/SVuxzGmWXwq+iIySUS2ish2Ebm7jv4XiMgqEakSkW/X6jdLRLY5r1n+Cu6GZ7/YRUxUBNePsZuxws3VI7sRGxXB3z/d4XYUY5qlwaIvIpHAE8BlwADgWhEZUGuwPcB3gZdqjZsMPACMBkYBD4hIh+bHbn2HT1Tw5qo8pg9LJyUh1u04ppV1aBvDzed1Z+7qfJbuOOh2HGOazJct/VHAdlXdqaoVwCvAVO8BVDVXVdcBtS9vuBT4QFUPqeph4ANgkh9yt7qXvtxDWWUNN51nN2OFqzsv7k235DbcO2895VXWEJsJTr4U/XTA+yLlPKebL3waV0RuFZEcEckpKirycdKtp6KqhtlLcjm/d4q1mR/G2sRE8tDUQewsOsGTn+50O44xTeJL0a/rukRfHyDq07iq+pSqjlTVkampqT5OuvW8s34fhSXl3Gxb+WHvor5pXD64C49/vJ1dxSfcjmNMo/lS9POAbl7fM4B9Pk6/OeMGBFXlH5/voldaAhf2CbwVkml9918xgNjICO6bt8EeoG6Cji9FfwXQW0S6i0gMcA2wwMfpvw9MFJEOzgnciU63oLF81yE27jvGTePtZizj0aldHD+5tC+LtxezYG1QbcMY03DRV9Uq4A48xXoz8JqqbhSRh0RkCoCInCsiecBVwJMistEZ9xDwMJ4VxwrgIadb0Hh68S46xEczY4SvpzFMOLh+TBZDMpJ4+O1NHC2tdDuOMT7z6Tp9VV2oqn1UtaeqPuJ0u19VFzifV6iAORn1AAARZElEQVRqhqq2VdWOqjrQa9xnVLWX83q2ZRajZeQWn+Bfmw/YzVjmayIjhEenD+bQiQp+8/4Wt+MY4zO7I/csnv1iF1ERwg12M5apw6D0JL47rjsvLd/Dyt2H3Y5jjE+s6NfjaGklr6/MY8rQdNLaxbkdxwSoH03sQ5ekOO6du55Ka4XTBAEr+vV4ecUeSiuq7TJNc1YJsVE8cMVAtuwv4ZnFu9yOY0yDrOjXobLaczPW2B4dGdC1ndtxTIC7dGAnvtE/jcf+tY28w6VuxzHmrKzo1+HdDfspOFrGLefbVr5pmIjw4BTPtQsPzN9o1+6bgGZFvxZV5enPd9IjpS0X9U1zO44JEhkd4vnvb/bmwy2FvL9xv9txjKmXFf1aVu4+zNq8o9w4PpuICLsZy/juxvHd6dc5kQcXbOJ4eZXbcYypkxX9Wp5evIukNtFceU6G21FMkImOjODRGYM5UFLGHxZtdTuOMXWyou9l76FS3t+4n+tGZxIfE+V2HBOERmR24LpRmcxeksuG/KNuxzHma6zoe3n2i1wiRJg1NtvtKCaI/WxSP5LbxnLP3PVU19hJXRNYrOg7SsoqeS1nL5OHdKFzkt2MZZouqU00903uz7q8ozy/NNftOMacwYq+49UVezleXsXN5/VwO4oJAVOGduX83in8ftFX7D9a5nYcY06zog9UVdfw7Be5jMpOZnBGkttxTAgQER6eOoiK6hoeenuj23GMOc2KPrBo0wHyj5y0598av8pOacudF/Vi4fr9fLyl0O04xgBW9AHPZZqZyfF8c0Ant6OYEHPrhT3omdqW++Zv4GSFPUzduC/si/7qPYdZufswN47PJtJuxjJ+FhsVySPTB5N3+CR//mib23GMsaL/9OJdJMZGcdXIbg0PbEwTjOnRkW+fk8H/fbaTrftL3I5jwlxYF/38Iyd5d8N+rh2dSUKs3YxlWs493+pPYlwU985dT41du29cFNZFf/aSXABmjct2NYcJfcltY/j5t/qTs/swr+XsdTuOCWNhW/RPlFfx8pd7mDSoM+nt27gdx4SBq87JYFT3ZP7n3S0UHy93O44JU2Fb9F/P2UtJWRW32GWappWICI9MG0RpRRWPvrPZ7TgmTIVl0a+uUZ75IpcRme0ZntnB7TgmjPTulMitF/Rgzup8lmwvdjuOCUNhWfT/tfkAew6VWpMLxhV3XtybzOR4fjFvA+VVdu2+aV1hWfSfXryL9PZtuHSg3YxlWl9cdCQPTxvEzuIT/O2THW7HMWEm7Ir++ryjfLnrEDeOzyYqMuwW3wSIC/ukMnlIF/768Q52Fh13O44JI2FX9Z5evJO2MZH8x7l2M5Zx1/2TBxAbFcEv5m2wh6mbVhNWRX//0TLeXlfA1edm0i4u2u04JsyltYvjZ5P6smTHQeatyXc7jgkTYVX0Zy/NpUaVG8dnux3FGACuG53F0G7t+dXbmzlSWuF2HBMGfCr6IjJJRLaKyHYRubuO/rEi8qrTf7mIZDvds0XkpIiscV5/929835VWVPHS8j1MHNCZbsnxbsUw5gyREcKj0wdx5GQlv3lvi9txTBhosOiLSCTwBHAZMAC4VkQG1BrsZuCwqvYC/gT8xqvfDlUd5rxu81PuRntzVT5HT1Zyy/l2M5YJLAO7JnHjuGxe/nIvObmH3I5jQpwvW/qjgO2qulNVK4BXgKm1hpkKzHY+vwFcIiIB005xTY3yzOJdDM1I4pwsuxnLBJ7//mYfuibFce/cDVRW17gdx4QwX4p+OuDdQlSe063OYVS1CjgKdHT6dReR1SLyqYicX9cMRORWEckRkZyioqJGLYAvPt5ayK7iE9x0XncCaF1kzGltY6N4cMpAth4o4R+f73I7jglhvhT9uqpk7evL6humAMhU1eHAj4CXRKTd1wZUfUpVR6rqyNTUVB8iNc7Ti3fRJSmObw3u4vdpG+MvEwd25psDOvG/H37F3kOlbscxIcqXop8HeF/UngHsq28YEYkCkoBDqlquqgcBVHUlsAPo09zQjbFx31GW7DjIrHHZRNvNWCbA/XLKQCJEuH++XbtvWoYvVXAF0FtEuotIDHANsKDWMAuAWc7nbwMfqaqKSKpzIhgR6QH0Bnb6J7pvnlmcS5voSK49N7M1Z2tMk3Rt34YffbMPH28t4t0N+92OY0JQg0XfOUZ/B/A+sBl4TVU3ishDIjLFGexpoKOIbMdzGOfUZZ0XAOtEZC2eE7y3qWqrXZ5QWFLGW2v3cdXIDJLi7WYsExy+Oy6bAV3a8cu3NlJSVul2HBNifHpGoKouBBbW6na/1+cy4Ko6xnsTeLOZGZvshaW7qayp4cbxdpmmCR5RkRE8OmMw0//6BX9Y9BUPThnodiQTQkL2IHdZZTUvLN/DJf060T2lrdtxjGmUYd3ac/3oLGYvzWVd3hG345gQErJFf+7qfA6dqOBmezKWCVI/ndSXlIRY7pm7niq7dt/4SUgWfVXl6cW7GNi1HWN6JLsdx5gmaRcXzf2TB7Ah/xjPLd3tdhwTIkKy6H/6VRHbC49zs92MZYLc5CFduKBPKr99fwsvLd9jl3GaZgvJov/04l2kJcYyeUhXt6MY0ywiwu+vGsKIzA7cM3c9s55dQcHRk27HMkEs5Ir+1v0lfL6tmFnjsomJCrnFM2EoLTGOF24ezUNTB7Ji1yEm/ukzXs/Za1v9pklCrio+s3gXcdERXDfKbsYyoSMiQpg5Npt3f3g+/Ton8tM31vG953IoPFbmdjQTZEKq6BcfL2fumnxmjMigQ9sYt+MY43fZKW155dax/OLy/ny+rZiJj33G/DX5ttVvfBZSRf/FZXuoqKrhJrsZy4SwyAjhlvN78M4Pzie7Y1t++Moavv/iKg4eL3c7mgkCIVP0yyqreX5ZLhf1TaVXWoLbcYxpcb3SEnjjtrH8bFJfPtxcyMQ/fcZ7GwrcjmUCXMgU/UMnKujTKZGbz+vhdhRjWk1UZATfn9CLt+48j85Jcdz2wiruemW1PW/X1EsC7VjgyJEjNScnx+0YxgSdyuoanvh4O49/tJ3ktjH8+srBXNyvk9uxTCsRkZWqOrKh4UJmS9+YcBcdGcFd3+jDvNvH0yE+hpv+mcNPX1/LMWup03ixom9MiBmUnsSCO8fz/Qk9eXNVHpP+9Bmfb/P/Y0hNcLKib0wIio2K5GeT+vHmf40jLiaSG57+knvnrudEeZXb0YzLrOgbE8KGZ3Zg4Q/O55bzuvPSl3uY9L+fsWznQbdjGRdZ0TcmxMVFR/KLyQN49daxRIhwzVPL+OVbGzlZUe12NOMCK/rGhIlR3ZN594fnM3NsFs9+kcu3/vw5K3cfdjuWaWVW9I0JI/ExUTw0dRAv3jKaiqoarvr7Ev7n3c2UVdpWf7iwom9MGBrfK4X37jqf/xjZjSc/3ckVf1nM+ryjbscyrcCKvjFhKjEuml9fOYRnbzyXY2WVTPvrF/xx0VYqquzRjKHMir4xYe6ivmksuutCpg7ryp8/2s60J75gc8Ext2OZFmJF3xhDUnw0f/yPYTx1wzkUlpQz5fHFPP7RNnsgewiKcjuAMSZwTBzYmZHZydw/fwO/X/QVizYdYObYbPp3SaR3WqI9jS4EWINrxpg6vbOugPvnb+DgCU+LnVERQq+0BPp3aUf/LokM6JJE/y6JdEyIdTmpAd8bXLMtfWNMnS4f0oVJgzqzq/gEmwuOnX4t3XGQuavzTw+XlhjrrAhOrQza0T2lLVGRtlcQiKzoG2PqFels3fdKS+CKoV1Pdz90ouL0SmBTwTE2F5SwZMdOKqs9Rw5ioyLo0ymR/l0SvVYI7UhqE+3WohiHHd4xxvhFRVUNO4qOe+0VlLC54Njpw0MA6e3b0L9LOwZ4rQwyk+OJiBAXk4cGvx7eEZFJwP8CkcA/VPXXtfrHAs8B5wAHgatVNdfp93PgZqAa+IGqvt+I5TDGBImYqIjThfwUVaWopJxNXnsEmwuO8dGWA9Q425ttYyLp2/nMPYJ+nRNpG2sHIlpCg7+qiEQCTwDfBPKAFSKyQFU3eQ12M3BYVXuJyDXAb4CrRWQAcA0wEOgK/EtE+qiq3fNtTBgQEdLaxZHWLo4JfdNOdy+rrOarAyWn9wg2FRxjwdp9vLh8jzMepCbE0iYmktioCGKjnPfoCOKiIomN9uoWFUFsdCRxzvvpbt7DRUecMZ246Fr9nXFEQn+Pw5dV6Shgu6ruBBCRV4CpgHfRnwo86Hx+A3hcPL/eVOAVVS0HdonIdmd6S/0T3xgTjOKiIxmS0Z4hGe1Pd1NV8o+cPL03kH/4JOVV1ZRX1Tivasorazh6spLySq9uVTWUVXrem3u0OsYp/hEinKr/gmfl5Xk/NaSnv3c3OaPbv1ceIs4L+fd3r2memokA/bu04/HrRjRvIRrgS9FPB/Z6fc8DRtc3jKpWichRoKPTfVmtcdNrz0BEbgVuBcjMzPQ1uzEmhIgIGR3iyegQzzcHNP7ZvqpKZbWesaIoq6x2VhBeK4/Kasqcd+8VSlnlv1csp851KqAKip5eoZzqhnc3r2GUM7txRjf1muaZ3VDITI5v/A/XSL4U/br2d2qvT+sbxpdxUdWngKfAcyLXh0zGGHMGESEmSoiJiiDR7TABzJcLafOAbl7fM4B99Q0jIlFAEnDIx3GNMca0El+K/gqgt4h0F5EYPCdmF9QaZgEwy/n8beAj9ewfLQCuEZFYEekO9Aa+9E90Y4wxjdXg4R3nGP0dwPt4Ltl8RlU3ishDQI6qLgCeBp53TtQewrNiwBnuNTwnfauA2+3KHWOMcY/dnGWMMSHA15uzrHEMY4wJI1b0jTEmjFjRN8aYMGJF3xhjwkjAncgVkSJgdzMmkQIU+ylOsAi3ZQ635QVb5nDRnGXOUtXUhgYKuKLfXCKS48sZ7FASbsscbssLtszhojWW2Q7vGGNMGLGib4wxYSQUi/5TbgdwQbgtc7gtL9gyh4sWX+aQO6ZvjDGmfqG4pW+MMaYeVvSNMSaMhEzRF5FJIrJVRLaLyN1u52lpItJNRD4Wkc0islFEfuh2ptYiIpEislpE3nY7S2sQkfYi8oaIbHH+3mPdztTSROS/nX/XG0TkZRGJczuTv4nIMyJSKCIbvLoli8gHIrLNee/g7/mGRNH3enj7ZcAA4FrnoeyhrAr4sar2B8YAt4fBMp/yQ2Cz2yFa0f8C76lqP2AoIb7sIpIO/AAYqaqD8DTpfo27qVrEP4FJtbrdDXyoqr2BD53vfhUSRR+vh7eragVw6uHtIUtVC1R1lfO5BE8h+Nrzh0ONiGQAlwP/cDtLaxCRdsAFeJ5ZgapWqOoRd1O1iiigjfMkvnhC8Il7qvoZnuePeJsKzHY+zwam+Xu+oVL063p4e8gXwFNEJBsYDix3N0mreAz4GVDjdpBW0gMoAp51Dmn9Q0Tauh2qJalqPvB7YA9QABxV1UXupmo1nVS1ADwbdkCav2cQKkXfpwewhyIRSQDeBO5S1WNu52lJIjIZKFTVlW5naUVRwAjgb6o6HDhBC+zyBxLnOPZUoDvQFWgrIte7myp0hErRD8sHsItINJ6C/6KqznE7TysYD0wRkVw8h/AuFpEX3I3U4vKAPFU9tRf3Bp6VQCj7BrBLVYtUtRKYA4xzOVNrOSAiXQCc90J/zyBUir4vD28PKSIieI7zblbVP7qdpzWo6s9VNUNVs/H8jT9S1ZDeAlTV/cBeEenrdLoEzzOnQ9keYIyIxDv/zi8hxE9ee1kAzHI+zwLm+3sGDT4YPRjU9/B2l2O1tPHADcB6EVnjdLtHVRe6mMm0jDuBF50Nmp3AjS7naVGqulxE3gBW4blKbTUh2CSDiLwMTABSRCQPeAD4NfCaiNyMZ+V3ld/na80wGGNM+AiVwzvGGGN8YEXfGGPCiBV9Y4wJI1b0jTEmjFjRN8aYMGJF3xhjwogVfWOMCSP/H9YnBPUq/hNAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def binomial(k, n, p): # n and p are parameters\n",
    "    return (np.math.factorial(n)/(np.math.factorial(k) * np.math.factorial(n - k))) * p**k * (1-p)**(n-k)\n",
    "\n",
    "# CHANGE THESE TO SEE DIFFERENT BINOMIAL DISTRIBUTIONS\n",
    "n = 10\n",
    "p = 0.3\n",
    "\n",
    "x = [i for i in range(n+1)]\n",
    "y = [binomial(i, n, p) for i in x]\n",
    "plt.plot(x, y)\n",
    "plt.title(f\"Binomial Distribution (n = {n}, p = {p})\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Poisson Distribution\n",
    "\n",
    "The Poisson distribution represents the expected number of events occuring in a fixed interval of time, when the rate at which events happen is $\\lambda$, and each event is independent of any others\n",
    "\n",
    "$$P(k; \\lambda) = e^{\\lambda}\\frac{\\lambda^k}{k!}$$\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\text{Mean:} &\\qquad \\lambda\\\\\n",
    "\\text{Variance:} &\\qquad  \\lambda\\\\\n",
    "\\text{Median:} &\\qquad \\lfloor \\lambda + 1/3 - 0.02\\lambda \\rfloor\\\\\n",
    "\\text{Support:} &\\qquad \\{0, 1, 2, \\dots\\}\n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEKCAYAAAD+XoUoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt4HHd97/H3V5Il3+2VLSeOJUu2c8MhTmKvA5RyK2lwQkhoAsGhtIHmPGlOm1M4fTglLW3gSWnLpaelFw5NTkmh4ZKYUKgJzgmhAc7hgQTLjmOwnYul+iLbsdeRb5Icy5K+54+ZdTbrlTQr7X0+r+fZR7szv5n5arT67uxvfvMdc3dERCQe6sodgIiIlI6SvohIjCjpi4jEiJK+iEiMKOmLiMSIkr6ISIwo6YuIxIiSvohIjCjpy7jMbJuZvbXccURVyHgz12Vmu8zsqkKsN3vdxWBmf2VmHynW+vNlZj83s0vKHUfcKenHSJi0TppZn5kdNLN/MbOZ4y3n7pe4+49KEOK4Mn6HE2Z21Mx+amZ3mNmZ93LUeKMk8UL+7tnbK+Z+NbMW4LeBezOm7TCznmIkXjNrMrMvmdnu8G/ztJldk9Xsr4F7Cr1tyY+Sfvy8y91nAiuB1cCfljmeiXiXu88C2oFPAx8DvlTIDZhZQyHXVwYfBDa4+8mMaa8FngduKsL2GoC9wFuAOcCfAevMrCOjzXrgbWa2sAjbl4iU9GPK3fcBjxIkAszsNWb2o/DoeZuZXZ9um3mEamYfM7N94dHcc2b29ox2OedFWPdHzWyrmR0zs4fMbGrE3+GYu68H3gfcambp32XceM3sAWAx8N3wm88fhct9zMy2Av1m1pDj28BqM9tuZkfCb0pnYjUzN7PzM15/2cw+Nc72rspoX8j9dA3w46z9NQz8BLgsyv7Nh7v3u/sn3X2Xu4+4+yPAfwKrMtq8DGwCri709iU6Jf2YMrM24FrgaTObAnwX+D6wAPhvwNfM7KKsZS4C7gRWh0fa7wB2jTUv4rpvBtYAS4AVBEepkbn7z4Ee4E1R43X33wL2EH7zcffPhovdArwTmOvuQzk295vhepYBFxLxm9IY20vHWuj9dCnwXNY2pgFrw2VHZWaPhB88uR6PRPh1MbNzCPbPtqxZOyjCh45Ep6QfP98xs6MER3w/Bv4SeD0wE/i0uw+6+xPAIwQJMNMw0AQsN7Mp4VFd1zjzoqz77919v7v3EiS+yyfwe+0HmvOIdzR/7+57s7pFMv1jOL8X+AvO3kcTVej9NBc4kTXtL4B9wLKxzuW4+3XuPneUx3Xj/SLhB9jXgK+4+7NZs0+EsUmZKOnHz7vDf952d/+9MLmdB+x195GMdruBRZkLuvtO4CPAJ4FDZvagmZ03zrwo634x4/kAQfLL1yKgN2q8Y9ibx/zdBL9fIRR6Px0BZqVfmNkbCL4p3AQcI+zWK7TwhPoDwCDBt6xss4Cjxdi2RKOkLxAcJbdljoAh6H/el93Q3b/u7r9KcBLVgc+MMy/yuifKzFYTJMef5BNv+PqsRcbZXFvG88UEv1/aADA94/W5eay70PtpK0H3CmHf//3AHeG3hGcYo4vFzB4Nzzvkejw6xnJGcEL9HOAmdz+do9lrwu1LmSjpC8BTQD/wR2Y2xYKx4+8CHsxsZGYXmdmvmVkT8DJwkqALZax5kdY9EWY228yuC9f1VXf/RdR4QweBpXlu9vfNrNXMmoE/AR7KmLcFeL+Z1ZvZGoKRLJnG2l6h99OGjO3fA/wsPLmajnPUfn13vyY875DrkT0MM9MXCZL6u3J1j4V/h1XA4xP5haQwlPQFdx8EricY8XEY+F/Ab+foj20iGCJ5mKCrYQFB4ht1Xh7rzsd3zewEQVfLx4G/AT6Uo91Y8QL8FfCn4QnKj0bc9tcJTrZ2h49PZcz7MEGiPkpwwvc7WcuOur0i7Kd/Ba41s7cA7wX+e8a8LRT4ZKqZtQO/S3Ce4cWMbwa/mdHseuBH7r4/50qkJEy3SxSpTWb2l8Ahd/98uWMBMLOngNvc/ZfljiXOlPRFRGJE3TsiIjGipC8iEiNK+iIiMVJxRaXmz5/vHR0d5Q5DRKSqbNq06bC7t4zXruKSfkdHB52dneUOQ0SkqpjZ7ijt1L0jIhIjSvoiIjGipC8iEiNK+iIiMaKkLyISI0r6IiIxoqQvIhIjFTdOv1yeePYgW/ZM/IY+VyxO8LaLFxQwIhGRwlPSD/2Pb27lpf5BzPJf1h3mTJvC03/269TVTWAFIiIloqQPHB0Y5KX+Qf7k2ou5/c3L8l5+Xede/ujhrXSl+rjgnFnjLyAiUibq0we6Uv0ALGuZyP24IdmeAKBz95GCxSQiUgxK+kBXqg+ApRNM+kvmz2DejEY6dynpi0hlU9IHulP9TKk32hLTJrS8mbGyPUHn7t4CRyYiUlhK+gRH+u3zZtBQP/Hdsbojwe6XBkidOFXAyERECktJH+hO9bF0/oxJrWNVezMAm3S0LyIVLPZJ//TwCHt6B1i2YGL9+WmvXTSbxoY69euLSEWLfdLf2zvA6WGf9JF+U0M9l7XO0QgeEalosU/63enhmpM80oegi2fb/mO8fHp40usSESmG2Cf99HDNZfMnn/ST7QlODzvP7J14OQcRkWKKfdLvTvUzb0Yjc6ZPmfS6VukiLRGpcJGSvpmtMbPnzGynmd2VY/4fmtl2M9tqZv9hZu0Z84bNbEv4WF/I4Auh+3DfhK/EzZaY0cj5C2bSuUsjeESkMo2b9M2sHvgCcA2wHLjFzJZnNXsaSLr7CuBh4LMZ8066++Xh4/oCxV0wXal+lrZM7iRupmR7gk27jzAy4gVbp4hIoUQ50r8S2Onu3e4+CDwI3JDZwN1/6O4D4csngdbChlkcR/oH6e0fLNiRPgRdPMdfHmJneK5ARKSSREn6i4C9Ga97wmmjuQ14NOP1VDPrNLMnzezduRYws9vDNp2pVCpCSIXRfThdc6eAR/odwUVaGq8vIpUoStLPVSA+Z9+FmX0ASAKfy5i82N2TwPuBz5vZWbWL3f0+d0+6e7KlpSVCSIWRrq450UJruXTMmx4UX9OVuSJSgaIk/R6gLeN1K7A/u5GZXQV8HLje3c8UoHH3/eHPbuBHwBWTiLegJltoLRczY1XYry8iUmmiJP2NwAVmtsTMGoG1wKtG4ZjZFcC9BAn/UMb0hJk1hc/nA28Ethcq+MkqRKG1XJJh8bVDJ14u6HpFRCZr3Gzn7kPAncBjwA5gnbtvM7N7zCw9GudzwEzgm1lDM18DdJrZM8APgU+7e8Uk/e5UH8sK2J+fdqb4mvr1RaTCRLpdortvADZkTbs74/lVoyz3U+DSyQRYLKeHR9j90gBXX3Juwdf92kWzaWqoo3P3Ea65dGHB1y8iMlGxvSJ3b+8AQyNe0OGaaUHxtbm6MldEKk5sk/4rI3cK370DsKojwbZ9xzg5qOJrIlI5Ypv0uwtYaC2XZHuCoRHnmR4VXxORyhHjpN/P/JmFKbSWS7r4moZuikgliW3S70r1sbRIR/kAc6er+JqIVJ7YJv3uw/0sW1Cc/vw0FV8TkUoTy6SfLrRWzCN9COrwHH95iBcOqfiaiFSGWCb9YhRayyV55qYq6uIRkcoQy6SfHq5ZjDH6mdrnTWf+zEZdmSsiFSOmSb+PKfVGawELreWSLr6mi7REpFLEMul3p/rpKEKhtVyS7c3s6VXxNRGpDLFM+l2pvqL356et6gjH66uLR0QqQOyS/unhEfa8NFD0/vy0154350zxNRGRcotd0t8TFlor5N2yxtLYUKfiayJSMWKX9LuLXGgtl6SKr4lIhYhh0i9uobVckh1B8bUte1V8TUTKK3ZJvyvVV9RCa7msXJwuvqaLtESkvGKX9LtT/SXrz0+bO72RCxbMVL++iJRd7JJ+V5HuizueZEeCzSq+JiJlFqukf6R/kCMDp4teaC2XVe0qviYi5RerpJ8utFbsksq5qPiaiFSCWCX9rkPhcM0yHOmr+JqIVIJ4Jf3DfTTW1xW90FouZkayvZmNOtIXkTKKV9I/1E/7vOklKbSWS7Ijwd7ekxw6ruJrIlIesUr63Yf7SlZzJ5dVZ/r11cUjIuURm6SfLrRWyvIL2S5JF19Tv76IlElskn6pC63l0thQx2Vtc3VlroiUTWySfveZWySW70gfgqGb2/YfV/E1ESmLSEnfzNaY2XNmttPM7sox/w/NbLuZbTWz/zCz9ox5t5rZC+Hj1kIGn4+uVPpm6OU70gcVXxOR8ho36ZtZPfAF4BpgOXCLmS3PavY0kHT3FcDDwGfDZZuBTwCvA64EPmFmicKFH113qo/5M5uYM610hdZyWbW4GYDOXeriEZHSi3KkfyWw09273X0QeBC4IbOBu//Q3QfCl08CreHzdwCPu3uvux8BHgfWFCb0/HSl+st6EjdtzvQpXHiOiq+JSHlESfqLgL0Zr3vCaaO5DXg0n2XN7HYz6zSzzlQqFSGk/HWnyjtcM9Oq9mY271HxNREpvShJ33JMy5mtzOwDQBL4XD7Luvt97p5092RLS0uEkPLTGxZaK/dJ3LRke4ITLw/x/KET5Q5FRGImStLvAdoyXrcC+7MbmdlVwMeB6939VD7LFlv3mZO4FZL0O8KLtDReX0RKLErS3whcYGZLzKwRWAusz2xgZlcA9xIk/EMZsx4DrjazRHgC9+pwWkm9MlyzMrp3FjdPZ/7MJjapX19ESqxhvAbuPmRmdxIk63rgfnffZmb3AJ3uvp6gO2cm8E0zA9jj7te7e6+Z/TnBBwfAPe5e8mErXal0obXppd50TkHxtYTKLItIyY2b9AHcfQOwIWva3RnPrxpj2fuB+ycaYCF0pfrpmD+d+rpcpxjKI9mR4P9se5GDx1/mnNlTyx2OiMRELK7I7U71laWG/liSHenx+uriEZHSqfmkf3p4hD295S20lssl581m6pQ6dfGISEnVfNJPF1qrlJO4aVPq67isda5O5opISdV80u86VFnDNTMlO4LiawODQ+UORURiouaTfvfh8L64FXakD5Bsb2ZYxddEpIRqPul3HaqMQmu5rFwcXKSlm6WLSKnUfNLvPtxfMeUXsqn4moiUWu0n/VRfRXbtpKWLrw2r+JqIlEBNJ/1KK7SWy+qOsPjaQRVfE5Hiq+mkny60VmnDNTMl28OLtNTFIyIlUNNJv6vCqmvm0tY8jZZZTWzSnbREpARqOul3p/orqtBaLq8UX9ORvogUX00n/a5UX8UVWstlVXuCniMnOXj85XKHIiI1rqaTfneqv+IKreWi4msiUio1m/TThdaWLajc/vw0FV8TkVKp2aS/+6Wg0Fo1HOmr+JqIlErNJv0zwzUXVH7SB1jd0cy2/cfpP6XiayJSPDWb9LtS6UJrld+9A7CqI8HwiPOMiq+JSBHVbNLvTvXRMquJ2VMrr9BaLisXJzDTRVoiUly1m/QP97N0fnUc5QPMmTaFCxfMUtIXkaKq2aTfVeGF1nJZ1ZHg6d0qviYixVOTSb+3f5CjFV5oLZdke4ITp1R8TUSKpyaTflcVFFrLRcXXRKTYajLpV0N1zVzSxdc6VXxNRIqkJpN+V6qfxoY6FiWmlTuUvJwpvqZyDCJSJDWZ9LtTfXTMq/xCa7kkO5rZd/QkLx5T8TURKbwaTfr9Vde1k5ZsD26Wrjo8IlIMkZK+ma0xs+fMbKeZ3ZVj/pvNbLOZDZnZe7LmDZvZlvCxvlCBj2ZwaITdvQNVcyVutuXnzWbalHp18YhIUTSM18DM6oEvAL8O9AAbzWy9u2/PaLYH+CDw0RyrOOnulxcg1kj29A4wPOJVe6Q/pb6Oy9rmqPiaiBRFlCP9K4Gd7t7t7oPAg8ANmQ3cfZe7bwVGihBjXl65RWJ1Jn0Ihm5uP6DiayJSeFGS/iJgb8brnnBaVFPNrNPMnjSzd+dqYGa3h206U6lUHqs+W3eVFVrLRcXXRKRYoiT9XENg8qkTsNjdk8D7gc+b2bKzVuZ+n7sn3T3Z0tKSx6rPVm2F1nJJF1/bqH59ESmwKEm/B2jLeN0K7I+6AXffH/7sBn4EXJFHfHnrSvVVVaG1XF4pvqYRPCJSWFGS/kbgAjNbYmaNwFog0igcM0uYWVP4fD7wRmD72EtNnLvTleqvmhunjGVVR4Kn9xxV8TURKahxk767DwF3Ao8BO4B17r7NzO4xs+sBzGy1mfUA7wXuNbNt4eKvATrN7Bngh8Cns0b9FFRv/yDHTp6u+iN9gNUdCfpODfHciyq+JiKFM+6QTQB33wBsyJp2d8bzjQTdPtnL/RS4dJIxRtZ9ODiJWwtH+unia5t297L8vNlljkZEakVNXZHbdSgstFYFN0MfT2tiGgtmNanipogUVE0l/e7D1VloLRczI9mh4msiUli1lfRTfSyZN6MqC63lsqo9KL524NjJcociIjWippJ+V6q/qi/Kynam+JqO9kWkQGom6Q8OjbCnd6Bqa+7kki6+pjo8IlIoNZP0jwwMcuE5s7h44axyh1Iw6eJrukhLRAol0pDNanDO7Kk8+uE3lTuMgku2N/PFH3fRf2qIGU018+cSkTKpmSP9WpUMi69tUfE1ESkAJf0Kt7I9KL6mk7kiUghK+hVu9tQpXHSOiq+JSGEo6VeBVe0qviYihaGkXwWSYfG1Z188Xu5QRKTKKelXgVeKr6lfX0QmR0m/CpwpvqaTuSIySUr6VSBdfE1H+iIyWUr6VULF10SkEJT0q8TqDhVfE5HJU9KvEq9ZqOJrIjJ5SvpVYkp9HZe3zdVFWiIyKUr6VSTZkWD7/uP0nRoqdygiUqWU9KvIqvYEIw5b9qj4mohMjJJ+FTlTfE1dPCIyQUr6VSRdfE0nc0VkopT0q4yKr4nIZCjpVxkVXxORyVDSrzIqviYik6GkX2VaE9M4Z3YTG3VlrohMgJJ+lTEzku3NbNqlETwikr9ISd/M1pjZc2a208zuyjH/zWa22cyGzOw9WfNuNbMXwsethQo8zla1J9h/7GX2H1XxNRHJz7hJ38zqgS8A1wDLgVvMbHlWsz3AB4GvZy3bDHwCeB1wJfAJM0tMPux4S6aLr6lfX0TyFOVI/0pgp7t3u/sg8CBwQ2YDd9/l7luBkaxl3wE87u697n4EeBxYU4C4Y+1M8TV18YhInqIk/UXA3ozXPeG0KCIta2a3m1mnmXWmUqmIq46vV4qv6UhfRPITJelbjmlRrwyKtKy73+fuSXdPtrS0RFx1vCU7Euw4oOJrIpKfKEm/B2jLeN0K7I+4/sksK2NIdjSr+JqI5C1K0t8IXGBmS8ysEVgLrI+4/seAq80sEZ7AvTqcJpN0xeK5mMFG9euLSB7GTfruPgTcSZCsdwDr3H2bmd1jZtcDmNlqM+sB3gvca2bbwmV7gT8n+ODYCNwTTpNJUvE1EZmIhiiN3H0DsCFr2t0ZzzcSdN3kWvZ+4P5JxCijSHYk+PbmfQwNj9BQr+vsRGR8yhRVLNneTP/gMM++eKLcoYhIlVDSr2K/cv486uuM7/3iQLlDEZEqoaRfxRbMmspbL2zhW5t6GBrOvi5ORORsSvpV7ubVbRw6cYofP6+L2kRkfEr6Ve7XLl7A/JmNPLRx7/iNRST2lPSr3JT6Om5c2coTzx4ideJUucMRkQqnpF8Dbk62MjTifPvpnnKHIiIVTkm/Bpy/YBYrF89lXWcP7rphuoiMTkm/RtycbGPnoT42qxaPiIxBSb9GXHfZeUxvrOebnTqhKyKjU9KvETObGnjnpQv57jP76Ve5ZREZhZJ+Dbl5dRv9g8Ns0BW6IjIKJf0akmxPsHT+DNapi0dERqGkX0PMjPcm29i46wjdqb5yhyMiFUhJv8bctHIR9XXGuk6N2ReRsynp15gFs6fytota+NZmFWETkbMp6degm5NtpE6c4kfPqQibiLyakn4NetvFC5g/s4mHdEJXRLIo6degKfV13LRyEU88e4hDJ14udzgiUkGU9GvUe5NtDI843968r9yhiEgFUdKvUecvmMmq9gTrOveqCJuInKGkX8Pel2yjK9XP5j1Hyh2KiFQIJf0adu2KhUxvrGfdRo3ZF5GAkn4Nm9nUwHUrFvLIVhVhE5GAkn6NuzkZFGH7noqwiQhK+jVvVXuCpS0zWKcbp4sISvo1z8y4OdlG5+4jdKkIm0jsRUr6ZrbGzJ4zs51mdleO+U1m9lA4/ykz6wind5jZSTPbEj7+qbDhSxQ3ninCpqN9kbgbN+mbWT3wBeAaYDlwi5ktz2p2G3DE3c8H/hb4TMa8Lne/PHzcUaC4JQ8LZk3lbRct4Fub9nFaRdhEYi3Kkf6VwE5373b3QeBB4IasNjcAXwmfPwy83cyscGHKZL1vdRuH+1SETSTuoiT9RUBmv0BPOC1nG3cfAo4B88J5S8zsaTP7sZm9KdcGzOx2M+s0s85USkmpGN56UQvzZzapi0ck5qIk/VxH7NnX9Y/W5gCw2N2vAP4Q+LqZzT6roft97p5092RLS0uEkCRfU+rruGmVirCJxF2UpN8DtGW8bgX2j9bGzBqAOUCvu59y95cA3H0T0AVcONmgZWJuVhE2kdiLkvQ3AheY2RIzawTWAuuz2qwHbg2fvwd4wt3dzFrCE8GY2VLgAqC7MKFLvpa1zCTZnuAhFWETia1xk37YR38n8BiwA1jn7tvM7B4zuz5s9iVgnpntJOjGSQ/rfDOw1cyeITjBe4e79xb6l5Dobl7dRreKsInEllXaEV8ymfTOzs5yh1Gz+k8NceVf/IB3rljIZ99zWbnDEZECMbNN7p4cr52uyI2ZGU0NXLfiPB7ZekBF2ERiSEk/hm5e3crA4DDf26oibCJxo6QfQysXJ1jWMkNj9kViSEk/hjKLsO08pCJsInGipB9TN65spaHO+PSjOxgeqayT+SJSPEr6MdUyq4k/u245P9hxiE99b3u5wxGREmkodwBSPrf+Sgd7egf40k/+k8XN0/nQG5eUOyQRKTIl/Zj7k2tfQ8+RAe55ZDuL5k7j6kvOLXdIIlJE6t6Jufo64/Pvu4IVi+bwBw8+zTN7j5Y7JBEpIiV9YVpjPf9862rmz2zitq90srd3oNwhiUiRKOkLEJzY/fKHVjM4NMyHvryRYydPlzskESkCJX054/wFs7j3t5Lsfqmf//rVTQwO6daKIrVGSV9e5Q3L5vGZm1bw066X+ON/+4VKMIvUGI3ekbPcuLKVvb0n+dsfPM/i5ul8+KoLyh2SiBSIkr7k9AdvP589vQP87Q+ep615GjeubC13SCJSAEr6kpOZ8Vc3XsqBYyf52Le2snDONN6wbN74C4pIRVOfvoyqsaGOL35gFR3zZvC7D3Sy89CJcockIpOkpC9jmjNtCvd/cDWNDfV88F82kjpxqtwhicgkKOnLuNqap3P/B5O81DfIf/nKRk4ODpc7JBGZICV9iWRF61z+bu3lbN13jI889LTKMYtUKSV9iezqS87l7uuW89i2g/zlhh3lDkdEJkCjdyQvH3rjkjPlmBvqjdvftJR5M5vKHZaIRKSkL3n703cu5+jAae79cTf/8pNdXHvpuXzg9e2sak9gZuUOT0TGYJV2mX0ymfTOzs5yhyERvHDwBF97ag/f2tTDiVNDXHzuLH7rDe28+/JFzGjS8YRIKZnZJndPjttOSV8ma2BwiH/fsp8Hfrab7QeOM7OpgRtXLuIDr2/nwnNmlTs8kVhQ0peSc3ee3nuUrz65m0e2HmBwaIQrlzTzgde3s+aSc2ls0LgBkWJR0pey6u0f5OFNe/nqk3vY0zvA/JmNrF29mFtet5hFc6eVOzyRmqOkLxVhZMT5fzsP88DPdvPEswcB+LWLz+HdV5zHkvkzaJ07ndnTGnQCWGSSoib9SGfbzGwN8HdAPfDP7v7prPlNwL8Cq4CXgPe5+65w3h8DtwHDwB+4+2N5/B5S5erqjLdc2MJbLmxh39GTfOOpPTy4cQ8/2HHwTJtZTQ0sSkyjNTGN1sR0Fs3NeJ6YRmL6FH0oiBTIuEf6ZlYPPA/8OtADbARucfftGW1+D1jh7neY2VrgN9z9fWa2HPgGcCVwHvAD4EJ3H/U6fh3p177BoRGeffE4+46cpOfISfYdPUnPkQF6wtd9p4Ze1X56Y/2rPgjOnTOV6Y31NDXU09RQx9Qpwc+mKXU0NdQzNfyZa1p9nT48pDYV8kj/SmCnu3eHK34QuAHYntHmBuCT4fOHgX+04NDsBuBBdz8F/KeZ7QzX97Oov4jUnsaGOla0zmVF69yz5rk7x08O0XM0+BBIfzD0HBlg39GTbN5zdFL3722oM+rqjDqDOjPqzDDAjHB6MM/C6Zmv6+ogmBpIf/nI/BhJfyN51UeL5Xyat8l829FHXXW4eOFs/uGWK4q6jShJfxGwN+N1D/C60dq4+5CZHQPmhdOfzFp2UfYGzOx24HaAxYsXR41dapCZMWf6FOZMn8Ml583J2ebk4DAvnx7m1NAIp4aCn2den841Ld02mDbiwYfLiDsjDiPu+JlpnJnu4fTMdmnpb8iZ35PTs189LWOZyeyYSSzsk9uylFBboviDHKIk/VwHCdnvotHaRFkWd78PuA+C7p0IMUmMTWusZ1pjfbnDEKlKUQZO9wBtGa9bgf2jtTGzBmAO0BtxWRERKZEoSX8jcIGZLTGzRmAtsD6rzXrg1vD5e4AnPPheux5Ya2ZNZrYEuAD4eWFCFxGRfI3bvRP20d8JPEYwZPN+d99mZvcAne6+HvgS8EB4oraX4IOBsN06gpO+Q8DvjzVyR0REiksXZ4mI1ICoQzZVDEVEJEaU9EVEYkRJX0QkRpT0RURipOJO5JpZCtg9iVXMBw4XKJxCUlz5UVz5UVz5qcW42t29ZbxGFZf0J8vMOqOcwS41xZUfxZUfxZWfOMel7h0RkRhR0hcRiZFaTPr3lTuAUSiu/Ciu/Ciu/MQ2rprr0xcRkdHV4pG+iIiMQklfRCRGqjLpm9kaM3vOzHaa2V055jeZ2UPh/KfMrKMEMbWZ2Q/NbIeZbTOzD+do81YzO2ZmW8LH3cWOK2Pbu8zsF+F2z6poZ4G/D/fZVjNbWYKYLsrYF1vM7LiZfSSrTUn2mZndb2aHzOyXGdOazew0t2oyAAAEPUlEQVRxM3sh/JkYZdlbwzYvmNmtudoUOK7Pmdmz4d/p22Z29n0nGf9vXoS4Pmlm+zL+VteOsuyY/79FiOuhjJh2mdmWUZYt5v7KmR/K8h4LbglXPQ+C8s5dwFKgEXgGWJ7V5veAfwqfrwUeKkFcC4GV4fNZBDeTz47rrcAjZdpvu4D5Y8y/FniU4G5nrweeKsPf9UWCC0xKvs+ANwMrgV9mTPsscFf4/C7gMzmWawa6w5+J8HmiyHFdDTSEzz+TK64of/MixPVJ4KMR/s5j/v8WOq6s+f8TuLsM+ytnfijHe6waj/TP3Kjd3QeB9I3aM90AfCV8/jDwdrNJ3FU6Anc/4O6bw+cngB3kuB9wBbsB+FcPPAnMNbOFJdz+24Eud5/M1dgT5u7/l+BeEJky30dfAd6dY9F3AI+7e6+7HwEeB9YUMy53/767D4UvnyS4I11JjbK/oojy/1uUuMIccDPwjUJtL6ox8kPJ32PVmPRz3ag9O7m+6kbtQPpG7SURdiddATyVY/YbzOwZM3vUzC4pVUwE9yb+vpltsuBG9Nmi7NdiWsvo/4zl2mfnuPsBCP5pgQU52pR7v/0OwTe0XMb7mxfDnWG30/2jdFWUc3+9CTjo7i+MMr8k+ysrP5T8PVaNSX8yN2ovOjObCXwL+Ii7H8+avZmg++Iy4B+A75QiptAb3X0lcA3w+2b25qz55dxnjcD1wDdzzC7nPouinPvt4wR3pPvaKE3G+5sX2heBZcDlwAGCrpRsZdtfwC2MfZRf9P01Tn4YdbEc0ya8z6ox6U/mRu1FZWZTCP6gX3P3f8ue7+7H3b0vfL4BmGJm84sdV7i9/eHPQ8C3Cb5mZyrnTeyvATa7+8HsGeXcZ8DBdBdX+PNQjjZl2W/hybzrgN/0sOM3W4S/eUG5+0F3H3b3EeB/j7K9cu2vBuBG4KHR2hR7f42SH0r+HqvGpD+ZG7UXTdhf+CVgh7v/zShtzk2fWzCzKwn2/0vFjCvc1gwzm5V+TnAi8JdZzdYDv22B1wPH0l87S2DUI7By7bNQ5vvoVuDfc7R5DLjazBJhd8bV4bSiMbM1wMeA6919YJQ2Uf7mhY4r8xzQb4yyvSj/v8VwFfCsu/fkmlns/TVGfij9e6wYZ6qL/SAYafI8wSiAj4fT7iH4JwCYStBVsBP4ObC0BDH9KsFXrq3AlvBxLXAHcEfY5k5gG8GIhSeBXynR/loabvOZcPvpfZYZmwFfCPfpL4BkiWKbTpDE52RMK/k+I/jQOQCcJjiyuo3gPNB/AC+EP5vDtkngnzOW/Z3wvbYT+FAJ4tpJ0Mebfp+lR6qdB2wY629e5LgeCN87WwmS2cLsuMLXZ/3/FjOucPqX0++pjLal3F+j5YeSv8dUhkFEJEaqsXtHREQmSElfRCRGlPRFRGJESV9EJEaU9EVEYkRJX0QkRpT0RURi5P8D2f1xGlzPmfgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def poisson(k, mean): # mean is a parameter\n",
    "    return np.exp(-mean) * mean**k / np.math.factorial(k)\n",
    "\n",
    "# CHANGE THIS TO SEE DIFFERENT POISSON DISTRIBUTIONS\n",
    "mean = 2\n",
    "x = [i for i in range(20+1)]\n",
    "y = [poisson(i, mean) for i in x]\n",
    "plt.plot(x, y)\n",
    "plt.title(f\"Poisson Distribution ($\\lambda$ = {mean})\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometric Distribution\n",
    "\n",
    "The geometric distribution represents the number of coin flips required before flipping your first heads, when the coins have a probability $p$ of coming up heads.\n",
    "\n",
    "$$P(k; p) = p(1-p)^k$$\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\text{Mean:} &\\qquad \\frac{1-p}{p}\\\\\n",
    "\\text{Variance:} &\\qquad  \\frac{1-p}{p^2}\\\\\n",
    "\\text{Median:} &\\qquad \\left\\lceil \\frac{-1}{\\log_2(1-p)}\\right\\rceil-1\\\\\n",
    "\\text{Support:} &\\qquad \\{0, 1, 2, \\dots\\} \n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X98XHWd7/HXZyZJ0yb9adK0tIU2oUJDQX6E+gPkeqXFglhUcCmyF1R8cLnIKnd1F1y9LHB1Rd31N7qisIqi/FLc4hb5pbiiSyFFbGlLoS0tTUubQEt/t2mSz/5xzrTT6UwySWdyMmfez8djHpk55zvnfHJm8p6T7znzPebuiIhIvCSiLkBERApP4S4iEkMKdxGRGFK4i4jEkMJdRCSGFO4iIjGkcBcRiSGFu4hIDCncZVCY2TIze9dQX3bmssxsrZnNLsSysy2/kMzsS2Z2bTGWPRjM7GkzOyHqOuJC4T6Emdl8M1tkZrvMrD28f7WZWdS1pcsnAN39BHd/YoDL3mNmO8zsDTP7k5ldZWYH3rv5LruYdea7vkIuP2Nd9cBlwPcLvew+1jvOzB4I36PrzOzDvbR9wsz2mtnO8LYyo8k/AzcXt+LyoXAfoszs08A3ga8CE4AG4CrgDKAqwtL6xcwqCrCY97n7SOAY4BbgOuD2Aiz3EAWqNSofARa6+55BXu+tQCfB+/NS4Ht97H1f4+614e24jHkLgP9pZhOLVGt5cXfdhtgNGA3sAi7so91RwC+ADuBl4JNp82YATwBvAMuAeRnPXQv8HbAkXNftBH+gDwE7gMeAsXms5ydAD7AH2An8fbjs68Jl7wMqwmmz0543BfhluMzXge/k+B0PeV44bVa4zpnZ2oTr3hD+HiuBs7PVmfbc9FrbMpa1FvgssBzYCvwbUB3Oc+DYtLY/Ar6Qa7tk1trbaxS2+0xY1zbgntR6c2yn3wJ/nfb4UuBP4fM2AeuBcwv8Pq0hCPY3Z7wfbsnR/gng430s81Hg8qj/BuNwi7wA3bK8KDAX6AIqemmTABYDNxDsyTcCa4D3AJXAKuAfwnnvDoPuuLTnrwWeIgj0SUA78CxwCjAsDIt/7G09GcvKDMTnCAJ8eGYbIAn8Bfh6GBDVwJk5fs9Dlp02/RXg/2RZ9nFhkB0VPp4KNOVaVmatOX6X58P544A/cjDAc4Z7H+ub3ddrFLZ7muCDdRywAriql/dDB3B62uNbCD5Y/ipc12eAdb08/9cEHzLZbr/O8ZxTgD0Z0z4DPJij/RNhna+F2/FdWdp8C/ha1H+DcbipW2ZoqgNec/eu1ISwr/mNsP/5LOB0oN7db3b3TndfA/wAmA+8Dagl2IPqdPffEvzxXpKxnm+7+2Z33wD8AVjk7n92933AAwR/vL2tpzffcvf1nr2bYBZBaP2du+9y973u/mS/thBsJAi9TN0EH07NZlbp7mvdffUR1ArBfxXr3X0L8EUO344Dkc9r9C133xiu90Hg5F6WN4bgwyHlRODr7n6vu+8H7gSONrPqbE929/PdfUyO2/k51llL8F9Fum3AyBztryPYOZgE3AY8aGZNGW12hL+LHCGF+9D0OlCX3gfs7u9w9zHhvARB//NRYeC/YWZvEOwFNhAE53p370lb5jqCP6p0m9Pu78nyuLaP9fRmfS/zphDsRXb10qYvk4AtmRPdfRVwLXAj0G5md5vZUX0sq7daM+evI9i+Ryqf12hT2v3dBK9HLls5NFRPBO5Pezwe2OnuewdWblY7gVEZ00Zx6IfMAe6+yN13uPs+d/8xwd77eRnNRhL8tyBHSOE+NP0XQf/vBb20WQ+8nLGHNdLdzyPYq52SfkYJcDRBP3R/9baelGwXBejtQgHrCfYiB3QA08xOJwjBrHv77v4zdz+T4IPJgS/3UVNfFzWYknb/aILtC0HgjkibN6Efyy3kawRB3/ybAcxsDEHNHWnzLyI4npKVmT2UdhZL5i3X814EKsxsetq0txAcP8iHA5lnfs0g6LKTI6RwH4Lc/Q3gJuC7ZnaRmdWaWcLMTiboo4agP3a7mV1nZsPNLGlmM8PgW0RwkPTvzawyPK/6fcDdAyint/WkbCb4d7s/y3wVuMXMasys2szO6OtJZjbKzM4n+D1+6u5Ls7Q5zszebWbDgL0E/4F0D7DOlE+Y2WQzG0fwX8s94fTngA+H22Qu8D8yntfb+gr5GgEsTFv/iQS/84fNrMLM3gtcTfDfTFbufq4fPIsl83ZujufsIjgofnP4Op5BsEPyk8y2ZjbGzN4TvtYVZnYpcBbwcFqbYcBpBAdV5Qgp3Icod/8K8LcEZ5+0EwTF9wn6Lf/k7t0EYXAywRksrwE/BEa7eycwDzg3nP5d4DJ3f2EAdeRcT1qzLwGfD7ttPtOPZR5LcGC0Dbi4l6c8aGY7CPb4Pwd8DfhojrbDCA4mvkbQrTGeIJD7XWeanwGPEBxIXgN8IZz+qfD3eIPg7JRfZTwv5/oK+RqF7gTOM7PhBOF+F/B2gu6am4D3u/vyAS67N1cTHIhuB35OcJB7GRz4byC17SsJtlvqgOrfhDWln+s+D3jC3TciR8zcdZk9kTgws38iCNnjgBfd/esRl9QvZrYIuMLdn4+6ljhQuIvEjJk9SXBK5m+irkWio24ZkfiZCQy0e0diQnvuIiIxpD13EZEYimygpLq6Op86dWpUqxcRKUmLFy9+zd3r+2oXWbhPnTqV1tbWqFYvIlKSzGxdPu3ULSMiEkMKdxGRGFK4i4jEkMJdRCSGFO4iIjGUV7ib2VwzW2lmq8zs+izzP2JmHWb2XHj7eOFLFRGRfPV5KqSZJQkugjuHYPS+Z8xsQZYR5u5x92uKUKOIiPRTPnvus4BV7r4mHKb0bnq/iERRta7dwi0PvYCGTRARyS2fcJ/EoZcZa+Pwy7UBXGhmS8zsfjObkmU+ZnalmbWaWWtHR0e2Jn1aumEb//r71XTs3Deg54uIlIN8wj3zMlhw+OXDHgSmuvtJwGPAj7MtyN1vc/cWd2+pr+/z27NZNdYHl5Fc07FrQM8XESkH+YR7G4deQ3IyB68hCYC7v+7uqV3pHxBcKqsomuqDq8yt7thZrFWIiJS8fML9GWC6mU0zsypgPrAgvYGZTUx7OA9YUbgSD3XU6OFUVya05y4i0os+z5Zx9y4zu4bgQrZJ4A53X2ZmNwOt7r4A+KSZzQO6gC3AR4pVcCJhTKurZY323EVEcsprVEh3X0hwdfX0aTek3f8s8NnClpZbY30NS9u2DdbqRERKTkl+Q7Wpvpa2rbvZu7876lJERIakEg33Gnoc1r2+O+pSRESGpJIM98a61OmQ6ncXEcmmNMM9PB1yzWs6Y0ZEJJuSDPeaYRVMGFXN6nbtuYuIZFOS4Q7QNL6G1dpzFxHJqmTDvbGuljXtOzWAmIhIFqUb7vU17NjXpQHERESyKNlwbwoHEFvdrq4ZEZFMJRvuB8+Y0UFVEZFMJRvuGkBMRCS3kg331ABiGvpXRORwJRvuEHTNaM9dRORwJR3uTfW1rNcAYiIihynxcK/BNYCYiMhhSjzcNYCYiEg2JR3u0+p0PVURkWxKOtxrhlUwcXS1DqqKiGQo6XCH4IwZ7bmLiByq9MO9rpY1Hbs0gJiISJqSD/em1ABiOzSAmIhISsmHe2NqADH1u4uIHFDy4d40PjwdUgOIiYgcUPLhPnFUNdWVCQ39KyKSpuTDPTWAmPbcRUQOKvlwh+Cgqk6HFBE5KBbh3lhfS9vWPRpATEQkFItw1wBiIiKHikm4p06HVNeMiAjEJNxTA4hpdEgRkUAswj01gJi+yCQiEsgr3M1srpmtNLNVZnZ9L+0uMjM3s5bClZif4JJ72nMXEYE8wt3MksCtwLlAM3CJmTVnaTcS+CSwqNBF5qOpvpbVGkBMRATIb899FrDK3de4eydwN3BBlnb/H/gKsLeA9eWtsa6GnRpATEQEyC/cJwHr0x63hdMOMLNTgCnu/uveFmRmV5pZq5m1dnR09LvY3qTGmFG/u4hIfuFuWaYd6PswswTwdeDTfS3I3W9z9xZ3b6mvr8+/yjw06nRIEZED8gn3NmBK2uPJwMa0xyOBmcATZrYWeBuwYLAPqqYGENMl90RE8gv3Z4DpZjbNzKqA+cCC1Ex33+bude4+1d2nAk8B89y9tSgV55BIGI11tdpzFxEhj3B39y7gGuBhYAVwr7svM7ObzWxesQvsj8b6Go0OKSICVOTTyN0XAgszpt2Qo+27jrysgWmqr+U/lr7K3v3dVFcmoypDRCRysfiGakpjOIDY2tfV7y4i5S1W4Z4aQEwHVUWk3MUq3BvrgwHEVrer311Eyluswn1EVTCA2JrXtOcuIuUtVuEOqTFmtOcuIuUtduEejA6pAcREpLzFLtyb6ms1gJiIlL3YhXvqoOoqdc2ISBmLYbjrdEgRkdiF+8RR1QyvTOqgqoiUtdiFeyJhTKur0Z67iJS12IU7BBfu0ABiIlLOYhnujXU1tG3dw9793VGXIiISiViGe9P4Wg0gJiJlLZbh3liXGmNG4S4i5Sme4R6e675GZ8yISJmKZbiPqKrgqNHVOh1SRMpWLMMdgi8zaXRIESlXsQ33Jg0gJiJlLLbh3hgOINauAcREpAzFONzDM2bU7y4iZSi24Z66nupqDUMgImUotuE+IRxATKdDikg5im24JxJGY32N9txFpCzFNtwhPB1Se+4iUoZiHe5N9TVseEMDiIlI+Yl1uDfWBwOIvawvM4lImYl3uNelxphRuItIeYl3uOtcdxEpU7EO99QAYjqoKiLlJq9wN7O5ZrbSzFaZ2fVZ5l9lZkvN7Dkze9LMmgtf6sAEl9xTt4yIlJc+w93MksCtwLlAM3BJlvD+mbuf6O4nA18BvlbwSgeosa6G1e07NYCYiJSVfPbcZwGr3H2Nu3cCdwMXpDdw9+1pD2uAIZOkjfW17Ors1gBiIlJW8gn3ScD6tMdt4bRDmNknzGw1wZ77J7MtyMyuNLNWM2vt6OgYSL39dmCMmXb1u4tI+cgn3C3LtMP2zN39VndvAq4DPp9tQe5+m7u3uHtLfX19/yodoANnzKjfXUTKSD7h3gZMSXs8GdjYS/u7gfcfSVGFNGFUNSOqktpzF5Gykk+4PwNMN7NpZlYFzAcWpDcws+lpD98LvFS4Eo9MImFMq6vRGTMiUlYq+mrg7l1mdg3wMJAE7nD3ZWZ2M9Dq7guAa8xsNrAf2ApcXsyi+6upvpZnX9kadRkiIoOmz3AHcPeFwMKMaTek3f9UgesqqMb6Gh5cspG9+7uprkxGXY6ISNHF+huqKRpATETKTVmEe5PGmBGRMlMW4T5No0OKSJkpi3AfUVXBpDHDNYCYiJSNsgh3QNdTFZGyUj7hXlfDmg4NICYi5aFswr1pfDCA2ObtGkBMROKvbMK9sS4YQEz97iJSDsom3JvG63RIESkfZRPuBwYQ00FVESkDZRPuZkZjvQYQE5HyUDbhDkG/u4b+FZFyUF7hXl/Dxm172NPZHXUpIiJFVVbh3qQBxESkTJRVuM+cNBqAp19+PeJKRESKq6zCfVpdDU31NTy6YnPUpYiIFFVZhTvA7OYGFq3ZwrY9+6MuRUSkaMou3M9pbqCrx3liZXvUpYiIFE3ZhfvJU8ZSV1vFI8vVNSMi8VV24Z5MGGcf38DvV3bQ2dUTdTkiIkVRduEOMKe5gZ37unhqjc6aEZF4KstwP3N6HdWVCR5V14yIxFRZhnt1ZZJ3Tq/nsRWbdfEOEYmlsgx3CLpmXt22l+c3bI+6FBGRgivbcD/7+PEkDB5dvinqUkRECq5sw/1NtcM47ZixPLpC57uLSPyUbbhD0DWz4tXtrN+yO+pSREQKqszDfQIAj2msGRGJmbIO92l1NRw7vlanRIpI7JR1uAPMntHAope3sG23BhITkfgo+3Cf09xAd4/zOw0kJiIxkle4m9lcM1tpZqvM7Pos8//WzJab2RIze9zMjil8qcVxypQx1NUO0xjvIhIrfYa7mSWBW4FzgWbgEjNrzmj2Z6DF3U8C7ge+UuhCiyWRMGbPGM/vV3awr0vXVhWReMhnz30WsMrd17h7J3A3cEF6A3f/nbunzid8Cphc2DKL6+BAYluiLkVEpCDyCfdJwPq0x23htFyuAB7KNsPMrjSzVjNr7ejoyL/KIjvj2DqGVyb1bVURiY18wt2yTMs62paZ/TXQAnw123x3v83dW9y9pb6+Pv8qiywYSKyOx5a3ayAxEYmFfMK9DZiS9ngysDGzkZnNBj4HzHP3fYUpb/DMaW5g03YNJCYi8ZBPuD8DTDezaWZWBcwHFqQ3MLNTgO8TBHtJnlN49owGDSQmIrHRZ7i7exdwDfAwsAK4192XmdnNZjYvbPZVoBa4z8yeM7MFORY3ZI2rqaLlmHG6tqqIxEJFPo3cfSGwMGPaDWn3Zxe4rkjMaW7giwtXsH7LbqaMGxF1OSIiA1b231BNN7u5AUBjzYhIyVO4p9FAYiISFwr3DHOaG3h6rQYSE5HSpnDPoIHERCQOFO4ZTp48hvqRw9Q1IyIlTeGeITWQ2BMr2zWQmIiULIV7FnOaG9jV2c1/rX496lJERAZE4Z7FO5qCgcR0bVURKVUK9yyqK5Oc9WYNJCYipUvhnsOc5gls2r6XpRu2RV2KiEi/KdxzePfx48OBxNQ1IyKlR+Gew7iaKlqmjlO4i0hJUrj34pzmBl7YtIP1W3b33VhEZAhRuPdi9gwNJCYipUnh3oupdTVM10BiIlKCFO59SA0k9sbuzqhLERHJm8K9DxpITERKkcK9D2+ZPIbxGkhMREqMwr0PiYRx9owGfr+yQwOJiUjJULjnYU7zeA0kJiIlReGeh3c01TGiKqmuGREpGQr3PFRXJjlrej2PrdhMT48GEhORoU/hnqc5zQ1s3r5PA4mJSElQuOfp3cePpyJh/PSpdVGXIiLSJ4V7nsbWVHHFmdO4b3EbrWu3RF2OiEivFO798KnZ05k0Zjife+B59nf3RF2OiEhOCvd+GFFVwY3zTmDl5h3c8eTLUZcjIpKTwr2f5jQ3MHtGA9947CXatmooYBEZmhTuA3DjvGYAbnpwecSViIhkp3AfgMljR3Dt7Ok8unwzjyzbFHU5IiKHySvczWyuma00s1Vmdn2W+WeZ2bNm1mVmFxW+zKHnY2dO47iGkdz04HJ2d3ZFXY6IyCH6DHczSwK3AucCzcAlZtac0ewV4CPAzwpd4FBVmUzwhQ/MZMMbe/jm4y9FXY6IyCHy2XOfBaxy9zXu3gncDVyQ3sDd17r7EqCszg88feo4Lm6Zwu1/eJkXNm2PuhwRkQPyCfdJwPq0x23htH4zsyvNrNXMWjs6OgayiCHn+nOPZ2R1BZ9/4HmNOyMiQ0Y+4W5Zpg0oxdz9NndvcfeW+vr6gSxiyBlbU8Vnz5tB67qt3L+4LepyRESA/MK9DZiS9ngysLE45ZSmi06dzKyp4/inh1awZZeutSoi0csn3J8BppvZNDOrAuYDC4pbVmlJJIwvfGAmO/d28aWFK6IuR0Sk73B39y7gGuBhYAVwr7svM7ObzWwegJmdbmZtwIeA75vZsmIWPRS9uWEkH39nI/ctbuPplzWwmIhEy9yjOQjY0tLira2tkay7WHZ3djHna/9JzbAkv/6bd1JVoe+IiUhhmdlid2/pq53Sp4BGVFVw07wTeHHzTm7XwGIiEiGFe4HNbm7gnOYGvvn4i6zfooHFRCQaCvci+Md5J5Aw48YFy4iq20tEypvCvQgmjRnOtbOn8/gL7TyyfHPU5YhIGVK4F8lHz5jG8RNGctOCZezap4HFRGRwKdyLpDKZ4IsfmMnGbXs1sJiIDDqFexGddsw45p8+hduffJkVr2pgMREZPAr3Irtu7vGMHl7J53+lgcVEZPAo3ItsbE0V/3DeDBav28pPF62LuhwRKRMK90Fw4amTeOf0Om7492V847EXtQcvIkWncB8EZsYPLmvhwlMn843HXuKqny5mp86gEZEiUrgPkurKJP/8oZO44fxmHn+hnQ9+94+se31X1GWJSEwp3AeRmfGxM6dx58dm0b5jH/O+80eefOm1qMsSkRhSuEfgjGPrWPCJM5kwqprL7ljED/+wRsMUiEhBKdwjcvSbRvDLq9/BnOYGvvAfK/j0fX9h7/7uqMsSkZhQuEeoZlgF37v0NP7v7Dfzy2c3cPFtT7Fp296oyxKRGFC4RyyRMD41ezrf/1+nsWrzDt73nSdZvG5r1GWJSIlTuA8R7zlhAr+8+gyGVya55LanuLd1fdQliUgJU7gPIcdNGMmCa85g1rRx/P39S7hxwTL2d/dEXZaIlCCF+xAzZkQVP/ro6Vxx5jR+9Ke1XHb702zZ1Rl1WSJSYhTuQ1BFMsH/O7+Zf/nQW1j8ylbmfedJlm/UqJIikj+F+xB24WmTufd/v5393T2899t/4KP/9jSPLd9Mt8amEZE+WFRfnmlpafHW1tZI1l1qXtu5jzv/tJa7n1lP+459HDW6mvmzjubi06fQMKo66vJEZBCZ2WJ3b+mzncK9dOzv7uHxFZu5a9Er/OGl10gmjDkzGrj0bUdzRlMdiYRFXaKIFFm+4V4xGMVIYVQmE8ydOZG5Myey9rVd/PyZV7ivtY3fLNvEMW8awYdnHc1Fp03mTbXDoi5VRCKmPfcSt6+rm988v4m7Fr3C0y9voSqZYO7MCVz61qOZNW0cZtqbF4kTdcuUoZc27+CuRa/wi2fb2LG3i2PH13LpW4/mg6dOZvTwyqjLE5ECULiXsT2d3Ty4ZCM/W/QKz61/g8qkcdyEkZw4aQwnTR7NiZNGc9yEkVQmdbKUSKlRuAsAz2/YxsKlr7J0wzaWtG1j2579AFRVJJgxcRQnTRrNiZNHc9Lk0RxbX0uFAl9kSNMBVQFg5qTRzJw0GgB3Z/2WPSzZ8AZL24Kw/9WfN/CTp4ILd1dXJjjhqGDPPrWH31hfS1Jn4YiUHO25l7meHmft67sO7NkvbdvG8xu3sbszGFt+eGWSSWOH0zBqGA0jqxk/qjq4H/4cP7Ka8aOGMawiGfFvIlIeCrrnbmZzgW8CSeCH7n5LxvxhwJ3AacDrwMXuvra/RcvgSySMxvpaGutrueDkSQB09zhrOnaypG0byzZu59Vte9i8fS+LXt5C+4697O8+fIdg7IhKGkaF4T9yWHh/GKOHV1I7rIKaYRUHftYMS1I7rILhlUmdzSNSJH2Gu5klgVuBOUAb8IyZLXD35WnNrgC2uvuxZjYf+DJwcTEKluJLJozpDSOZ3jCSC087dJ67s3X3fjZv38vm7Xtp374vuL9jL5u376N9+15e3LSDjp37+hwmIWFQU3Vo4NekfRBUVyYZVpGgMmlUVSSoTCaoqkhQFf6sTAb3Kw9MM6qSSSqTRmVFgqQZycTBWyL12IxkMviZSEDSjIpEIrgftgtuwXVvUz9FSkk+e+6zgFXuvgbAzO4GLgDSw/0C4Mbw/v3Ad8zMXBcGjR0zY1xNFeNqqpgxcVTOdt09zuu79rF9z3527utm174udu7rYte+LnZ1Bo8Pmbav+8D9Lbt2s6uziz2dPezv7qGzK/jZFfGYOmYcDH0Ms/RphoVtzIJ5QDgtmJdaBqnnpk2zjGmQ/QMlfdLBdVj2+Yc999Aphy39CD6/juSjr1Q/OI+k6k+ePZ33veWogtWSTT7hPglIv3JEG/DWXG3cvcvMtgFvAl5Lb2RmVwJXAhx99NEDLFlKQTJhQX/8yMKNfdPd40HYd/ewvyv10+ns7qazy4PH4YdBZ3cPPT1Od4/T405X2v3unuBYQ1eP0+1+oF13+Dj1H0dPj9Pj4IQ/3XGHHj843cPpPeH01GMAh+Ax4WM/OI30aWnLSkndPXTa4Q3SP+7S96UyPwYzd7MOnz/wD84j+sgt0d0/P8LCB+N7J/mEe7YPqMzfLJ82uPttwG0QHFDNY90iBwTdK0mqK3XwVqQv+ZzU3AZMSXs8GdiYq42ZVQCjgS2FKFBERPovn3B/BphuZtPMrAqYDyzIaLMAuDy8fxHwW/W3i4hEp89umbAP/RrgYYJTIe9w92VmdjPQ6u4LgNuBn5jZKoI99vnFLFpERHqX13nu7r4QWJgx7Ya0+3uBDxW2NBERGSgNJCIiEkMKdxGRGFK4i4jEkMJdRCSGIhsV0sw6gHUDfHodGd9+HSJUV/+orv4bqrWprv45krqOcff6vhpFFu5Hwsxa8xnycrCprv5RXf03VGtTXf0zGHWpW0ZEJIYU7iIiMVSq4X5b1AXkoLr6R3X131CtTXX1T9HrKsk+dxER6V2p7rmLiEgvFO4iIjE0pMPdzOaa2UozW2Vm12eZP8zM7gnnLzKzqYNQ0xQz+52ZrTCzZWb2qSxt3mVm28zsufB2Q7ZlFaG2tWa2NFxna5b5ZmbfCrfXEjM7dRBqOi5tOzxnZtvN7NqMNoO2vczsDjNrN7Pn06aNM7NHzeyl8OfYHM+9PGzzkpldnq1NAWv6qpm9EL5OD5jZmBzP7fU1L1JtN5rZhrTX67wcz+3177cIdd2TVtNaM3sux3OLss1yZUNk76/g0mFD70YwvPBqoBGoAv4CNGe0uRr41/D+fOCeQahrInBqeH8k8GKWut4F/DqCbbYWqOtl/nnAQwRXznobsCiC13QTwZcwItlewFnAqcDzadO+Alwf3r8e+HKW540D1oQ/x4b3xxaxpnOAivD+l7PVlM9rXqTabgQ+k8dr3evfb6Hrypj/L8ANg7nNcmVDVO+vobznfuDC3O7eCaQuzJ3uAuDH4f37gbOtyFfbdfdX3f3Z8P4OYAXBNWRLwQXAnR54ChhjZhMHcf1nA6vdfaDfTD5i7v6fHH6VsPT30Y+B92d56nuAR919i7tvBR4F5harJnd/xN27wodPEVwBbdDl2F75yOfvtyh1hRnwV8DPC7W+PGvKlQ2RvL+GcrhnuzB3ZogecmFuIHVh7kERdgOdAizKMvvtZvYXM3vIzE4YpJIceMTMFltwMfJM+WzTYppP7j+4KLZXSoO7vwrBHygwPkubKLfdxwj+48qmr9e8WK4Ju4zuyNHNEOX2eiew2d1fyjG/6NssIxsieX8N5XAv2IW5i8HMaoHXB2dSAAACTUlEQVRfANe6+/aM2c8SdD28Bfg28KvBqAk4w91PBc4FPmFmZ2XMj3J7VQHzgPuyzI5qe/VHJNvOzD4HdAF35WjS12teDN8DmoCTgVcJukAyRfZeAy6h9732om6zPrIh59OyTDui7TWUw33IXpjbzCoJXry73P2XmfPdfbu77wzvLwQqzayu2HW5+8bwZzvwAMG/xuny2abFci7wrLtvzpwR1fZKsznVPRX+bM/SZtC3XXhQ7XzgUg87ZjPl8ZoXnLtvdvdud+8BfpBjnZG818Ic+CBwT642xdxmObIhkvfXUA73IXlh7rA/73Zghbt/LUebCam+fzObRbCdXy9yXTVmNjJ1n+CA3PMZzRYAl1ngbcC21L+LgyDn3lQU2ytD+vvocuDfs7R5GDjHzMaG3RDnhNOKwszmAtcB89x9d442+bzmxagt/TjNB3KsM5+/32KYDbzg7m3ZZhZzm/WSDdG8vwp9xLjAR5/PIzjivBr4XDjtZoI3PEA1wb/5q4CngcZBqOlMgn+XlgDPhbfzgKuAq8I21wDLCM4QeAp4xyDU1Riu7y/hulPbK70uA24Nt+dSoGWQXscRBGE9Om1aJNuL4APmVWA/wd7SFQTHaR4HXgp/jgvbtgA/THvux8L32irgo0WuaRVBH2zqPZY6K+woYGFvr/kgbK+fhO+fJQTBNTGztvDxYX+/xawrnP6j1Psqre2gbLNesiGS95eGHxARiaGh3C0jIiIDpHAXEYkhhbuISAwp3EVEYkjhLiISQwp3EZEYUriLiMTQfwO4Dhk5wKJNsgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def geometric(k, p): # p is a parameter\n",
    "    return p*(1-p)**k \n",
    "\n",
    "# CHANGE THIS TO SEE DIFFERENT POISSON DISTRIBUTIONS\n",
    "p = 0.5\n",
    "x = [i for i in range(20+1)]\n",
    "y = [geometric(i, p) for i in x]\n",
    "plt.plot(x, y)\n",
    "plt.title(f\"Geometric Distribution ($p$ = {p})\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normal Distribution\n",
    "\n",
    "The normal, or Gaussian Distribution (also known as the bell curve), is one of the most commonly encountered probability distributions. This is because of the central limit theorem, which we'll discuss next. The normal distribution is also the limit of Poisson and Binomial distributions, as certain parameters go to zero or infinity. The nice analytical properties of the probability distribution function mean that most people approximate noise or measurement error as Gaussian.\n",
    "\n",
    "$$P(x; \\mu, \\sigma) = \\frac{1}{\\sqrt{2\\pi\\sigma^2}}e^{-\\left[(x-\\mu)/\\sigma\\right]^2}$$\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\text{Mean:} &\\qquad \\mu\\\\\n",
    "\\text{Variance:} &\\qquad \\sigma^2\\\\\n",
    "\\text{Median:} &\\qquad \\mu \\\\\n",
    "\\text{Support:} &\\qquad \\mathbb{R}\n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEKCAYAAAD+XoUoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VfW57/HPkxlCEggZgIQwBkgABQk4tM4BQS3UTmJta3t66rXVtuf09lZ77NFT257T4d5jz+CxWrXWtopTtWhRBBUnRAkyhjCEMCQQMgAhCSE703P/2Ct2ExKyk+xk7eF5v177xd5rr+FZm+S7V37rt9ZPVBVjjDGRIcrtAowxxgwdC31jjIkgFvrGGBNBLPSNMSaCWOgbY0wEsdA3xpgIYqFvjDERxELfGGMiiIW+6TcROSAihQNYvlhErghQLWesa6C19bb+QBKRfxORfxiMdYcrEflQRGa6XUcostAPYk5wVYlIos+0vxeRdS6W5Ren9tMi0iAidSKyXkRuE5GPf+ZUdaaqrvNzXecMcH/X5Y/uthfI9XfZVjrwFeChQK+7l+2misgLInJKRA6KyBeHcvu9EZE7RKRIRDwi8ng3s/xf4L4hLissWOgHvxjguwNdiYjEBKCWvvqUqiYBE4CfA3cCjwZyAy7tVyB9FVilqqeHeLsPAC1AJnAz8GCQHTkfAX4KPNbD+yuBK0Vk7NCVFB4s9IPfr4Dvi8jIrm+ISJ6IrHOOpItFZGmX9w+IyJ0isg04JSIxzrT/IyLbnKO8R0UkU0RecY7K14rIKGf5u0RknzN9p4jc0J8dUNWTqroSuBG4RURm+dT38RG1U+thZ3u7ReRqEfkDkAO8JCKNIvKDc+xX178G5jt1nxCR34lIgrMdFZGpPtt9XER+6jw/a3tdaz3X5+7M933n8z0pIk93brcHS4C3fJa/W0Qe9Hk9SkRae1lHnzh/OX4W+GdVbVTVd/GG6Jf9XD5WRH7m7Gur83mqiGwNVI2q+mdVfRE41sP7zcAmYFGgthkpLPSDXxGwDvi+70QRiQVeAl4DMoBvA38Skeldlr8JuA4YqaptzrTPAguBacCngFeAfwLS8P5MfMeZbx9wKZAC/Bj440COrFT1Q6DCWecZnLrvAOY7fx1cAxxQ1S8Dh/D+1TBCVX95jv3q6mZnPVOcff2RHzX2tL3OOv353L8ALAYmAefhPZrvyWxgd5fXW3xezwF2OyF3FhF52fny6e7xcg/bnAa0q+oen2lbAX+P9H8KXI33/3Ek8DrwAnDWQUE/6/NXCXD+ANcRcSz0Q8M9wLed9t9OFwEjgJ+raouqvgG8jDcMff2nqpZ3aT74L1WtUtXDwDvAB6q6WVU9eH955wKo6rOqekRVO1T1aWAvsGCA+3IESO1mejsQD+SLSKyqHlDVfedYT3f71dV/O/McB37G2Z9Nf/jzuf+n87kdx/sFMecc6xsJNPi87i70ezyCVtXrVXVkD4/re1hsBHCyy7STQNI56gRARJLwHhR82flsTwHPA6mqWhag+vzVgPfzM31goR8CVHUH3mC5y2fyOKBcVTt8ph0EsrosXt7NKqt8np/u5vUIABH5iohs6TwyA2bh/WtgILKA410nqmop8A/AvwDVIrJCRMadYz3d7de55jmI9zMbKH8+96M+z5twPs8enMAJWxGJw/tXyXaf98/nzC+BQGgEkrtMS+bML5+eXAaUqepen2mjOHOfh0oSUOfCdkOahX7ouBf4Bn8LlyPAePHpDYO3Lfpwl+X6NWCCiEwAfou3yWW0qo4EdgDSn/U565yPt/53u3tfVZ9U1U/iPfGrwC863+pudj82Od7neQ7ezwy8QTzc570xfVi3v5+7v7bhbW4ByAcOq2oTgIgIcAXnONJ3zsU09vB4pYfF9gAxIpLrM+18oNiPetPxflF1bl/wNut021TTz/r8lcc5PhvTPQv9EOEcCT/N39rbPwBOAT9wTqxdgbd9fkWANpmIN/xqAETka3iP9PtMRJJF5Hqntj+q6vZu5pkuIleJSDzQjPcvjnbn7Spgcj82fbuIZItIKt5zFk8707cAXxSRaBFZDFzeZblzbS/Qn/sqn+3PBjJEZIqIDAN+gvcL8EBPC6vqEufcQ3ePJT0scwr4M3CfiCSKyCeAZcAf4OMT24/3sMkdwAUiMsep8d/w/pw83d3M/anPqSHGOXkdDUSLSIL49NRyfk7mAWt6WofpnoV+aLkPbxijqi3AUry9P2qB/wG+oqq7ArEhVd0J/D/gfbwhOBt4r4+reUlEGvA2s9wN/DvwtR7mjcfbrbMWb1NBBt6gBm+w/MhpZvp+D8t350m8J1zLnMdPnenfxRvUdXhP9r7YZbketzcIn/sTwLVOgM4GVuM9sV6K93Mvw/vZBdq3gGFANfAU8E1V7TzSH08P/9eqWoT3/Mgqp7YxwLWq2hrg+n6E94v/LuBLznPfE/FLgXWqeqSbZc05iA2XaIy7RORf8YbvNcAjqvq8i7XE4W0yOW8QgjxgROQD4OvO+S7TBxb6xgQJEakAFjl/ZRkzKCz0jQkC4r0grgpIDOYjbBP6LPSNMSaC2IlcY4yJIEF3s6q0tDSdOHGi22UYY0xI2bRpU62qpvc2X9CF/sSJEykqKnK7DGOMCSkictCf+ax5xxhjIoiFvjHGRBALfWOMiSAW+sYYE0Es9I0xJoL4Ffoisli8w9eVishd55jvc86waQU+037oLLdbRK4JRNHGGGP6p9cumyISjXcQ5YV4h7rbKCIru94fxGdEnQ98puUDy/EOwzYOWCsi01S1HWOMMUPOn376C4DSzqHQRGQF3ntvd70p1E+AX3LmWK7LgBXOMHz7RaTUWd/7Ay3cmKFUVtPIlvI6Kk82k5wQw+T0ESyYlEpstLWQmtDiT+hnceawcxXAhb4ziMhcYLyqvtzl/uNZwIYuy3Ydzg8RuRW4FSAnJ8e/yo0ZAu/ureXXa/dQdPDEWe+lJ8XzpQsn8L8un0xCbLQL1RnTd/6EfnfD4318lzZn2Lj7ga/2ddmPJ6g+DDwMUFBQYHeAM6473dLOPX/ZwbObKhifOoy7lsygMC+DrJHDqW9uZUt5Hc9sLOf+tXv4y5bD3H/jHM4fb2N0m+DnT+hXcOZYo9n8baxR8A5OPAtY5x0ukzHAShFZ6seyxgSd46da+LvHN7K1oo5vXTGF71yde8aR/LC4aK6ZOYZrZo7hnb013PX8dpY/vIEHbp7LVTMyXazcmN750yC5EcgVkUnOqDrLgZWdb6rqSVVNU9WJqjoRb3POUmdYtZXAchGJF5FJQC7wYcD3wpgAqWtq4caH3qeksp4Hb57HDxbPOGfTzaW56bx4+yfIzRzBN57YxKs7jg5htcb0Xa+hr6ptwB14x+4sAZ5R1WIRuc85mj/XssXAM3hP+r4K3G49d0ywam5t59YnNnHwWBO/++p8Fs8a49dy6UnxPPWNi5idlcJ3V2xmUzft/8YEi6AbRKWgoEDtLpvGDXc+t42ni8r5z5vmsvT8cX1e/lijh888uJ7G5jZe+e6lZCQnDEKVxnRPRDapakFv81l/M2OAlVuP8HRROXdcObVfgQ8wekQ8j95SQFNLO//4zBY6OoLrgMoYsNA3hooTTdz95+1ckDOSfyjMHdC6pmYk8S9L83mv9BgPv1MWoAqNCRwLfRPRVJV/fnEH7ar8x/K5xATgYqsvFIxn8cwx3L9mDwePnQpAlcYEjoW+iWiv7DjKm7tr+N7CaYxPHR6QdYoIP142k9joKH704g6C7byZiWwW+iZiNbW08eOXipk5LpmvXjIxoOvOTE7gB4un887eWv66vTKg6zZmICz0TcR65J39VNV7uG/ZzIA063R184UTmDEmiV++uhtPm/VUNsHBQt9EpJoGDw+9tY/FM8cwb0LqoGwjOkr44bV5HDrexB83HBqUbRjTVxb6JiL91xt78bR1cOeSGYO6ncunpXNpbhr/9cZeGppbB3VbxvjDQt9EnMqTp1nxYTlfmD+eSWmJg769H1wzg7qmVp54/+Cgb8uY3ljom4jzm3X76FDlW1dMGZLtzc5O4aoZGfz2nTIaPW1Dsk1jemKhbyJKVX0zT20s53PzsskeFZgumv749lVTnaP9A0O2TWO6Y6FvIspj7+6nrb2Db10xdUi3OzdnFJfmpvHYu/tpbrWePMY9FvomYjR62njyw0MsmT2WnNFDd5Tf6bbLp1Db2MJfthwe8m0b08lC30SMpzeW09DcxjcunezK9i+ZMpq8sck88s5+u0rXuMZC30SEtvYOHnt3P/MnjmKOS8Maigi3XjaJvdWNrNtT40oNxvgV+iKyWER2i0ipiNzVzfu3ich2EdkiIu+KSL4zfaKInHambxGR3wR6B4zxx+u7qjlcd5qvf3KSq3Vcf9440pPi+aN13zQu6TX0RSQaeABYAuQDN3WGuo8nVXW2qs4Bfgn8u897+1R1jvO4LVCFG9MXf3j/IONSEijMc3cM29joKG4sGM+bu71fQsYMNX+O9BcApapapqotwApgme8Mqlrv8zIRsAZLEzT21TTybmktX7wwZ1DusdNXyxeMR4GnP7RbM5ih589vQBZQ7vO6wpl2BhG5XUT24T3S/47PW5NEZLOIvCUil3a3ARG5VUSKRKSopsbaOk1g/XHDQWKjhRvn57hdCgDZo4Zz+bR0ni4qp629w+1yTITxJ/Slm2lnHcmr6gOqOgW4E/iRM7kSyFHVucD3gCdFJLmbZR9W1QJVLUhPT/e/emN64Wlr54XNh1k8ayzpSfFul/OxLy7Ioarew+u7qt0uxUQYf0K/Ahjv8zobOHKO+VcAnwZQVY+qHnOebwL2AdP6V6oxfbd2ZzV1Ta18oSDb7VLOcNWMDMYkJ/DkB9bEY4aWP6G/EcgVkUkiEgcsB1b6ziAivgOLXgfsdaanOyeCEZHJQC5gA4eaIfPspnLGpSRwyZQ0t0s5Q0x0FF+YP56399ZQfrzJ7XJMBOk19FW1DbgDWA2UAM+oarGI3CciS53Z7hCRYhHZgrcZ5xZn+mXANhHZCjwH3KaqxwO+F8Z04+jJZt7eU8NnLsgmOqq7Vkp3LZ8/HgFWbLSjfTN0YvyZSVVXAau6TLvH5/l3e1jueeD5gRRoTH/9eXMFHQqfmxdcTTudxo0cxpXTM3i2qILvLZwelF9MJvy433/NmEGgqjxXVMGCialMHIJ75vfXZ+dlU93gYf2+WrdLMRHCQt+EpY8OnaCs9hSfC7ITuF1dNSODpPgYXthsN2EzQ8NC34SlZ4sqGB4XzXWzx7pdyjklxEZz7eyxrN5xlKYWG2DFDD4LfRN2Tre08/K2Sq6dPZbEeL9OW7nqhguyONXSzpqdVW6XYiKAhb4JO2tLqmj0tPHZC4K7aafTgompjEtJ4EVr4jFDwELfhJ2Xth4hMzmeBZNS3S7FL1FRwrK5Wby9t5baRo/b5ZgwZ6Fvwkp9cyvrdtdw3exxIdUF8oa5WbR3KC9vPdfF7sYMnIW+CSuvFVfR0t7Bp84P7hO4XU3LTCJ/bDIvbLHQN4PLQt+ElZe2HmF86jDXRscaiBvmZrG1vI6ymka3SzFhzELfhI3jp1p4t7SWT503DpHQadrp9KnzxwHw122VLldiwpmFvgkbr+yopL1DPw7PUDMmJYF5E0bx1+0W+mbwWOibsPHS1iNMzRjBjDFJbpfSb9fOHsuuow3WxGMGjYW+CQtV9c18sP94yDbtdLp29hgAVtnRvhkkFvomLPx1WyWqcH2I9drpamzKMC7IGclftx91uxQTpiz0TVhYtb2SGWOSmJI+wu1SBuza2WMpqaxnf+0pt0sxYchC34S86vpmNh06wZJZoX2U3+la5yZx1sRjBoNfoS8ii0Vkt4iUishd3bx/m4hsF5EtIvKuiOT7vPdDZ7ndInJNIIs3BuC1nVWowhKnPTzUjRs5jLk5Iy30zaDoNfSdMW4fAJYA+cBNvqHueFJVZ6vqHOCXwL87y+bjHVN3JrAY+J/OMXONCZRXdxxlcloiuRmh37TT6dpZYyk+Us/BY9bEYwLLnyP9BUCpqpapaguwAljmO4Oq1vu8TATUeb4MWKGqHlXdD5Q66zMmIOqaWni/7BjXzBoT0r12uur8q8X67JtA8yf0s4Byn9cVzrQziMjtIrIP75H+d/q47K0iUiQiRTU1Nf7WbgxrdlbR3qEsmRUeTTudskcN5/zsFFbvsF48JrD8Cf3uDp/0rAmqD6jqFOBO4Ed9XPZhVS1Q1YL09HQ/SjLGa3XxUcalJDA7K8XtUgJu0cwxbK04ydGTzW6XYsKIP6FfAYz3eZ0NnOtWgCuAT/dzWWP81uhp4+29tWHXtNNpUX4mAGtKbEQtEzj+hP5GIFdEJolIHN4Tsyt9ZxCRXJ+X1wF7necrgeUiEi8ik4Bc4MOBl20MvLmrmpa2jrDpqtnV1IwRTEpL5LVia+IxgdPrAKKq2iYidwCrgWjgMVUtFpH7gCJVXQncISKFQCtwArjFWbZYRJ4BdgJtwO2q2j5I+2IizKvFR0kbEce8CaPcLmVQiAiL8jN57L391De3kpwQ63ZJJgz4NWq0qq4CVnWZdo/P8++eY9mfAT/rb4HGdKe5tZ03d1WzbE5WSI2Q1VcL8zN56O0y1u2uYWmI3j3UBBe7IteEpPfLjtHU0s6imZlulzKo5uaMIm1EnDXxmICx0Dchae3OKobHRXPx5NFulzKooqOEwrxM1u2uwdNmLaNm4Cz0TchRVdaWVHFZbjoJseF/gfeimZk0etrYUHbc7VJMGLDQNyFnx+F6quo9FOaHd9NOp0umpDE8LtqaeExAWOibkLNm51GiBK6akeF2KUMiITaay6els2ZnFR0dZ13baEyfWOibkLOmpJqCCamkJsa5XcqQWTQzk+oGD9sOn3S7FBPiLPRNSKk40URJZT2F+ZFxlN/pyukZREeJNfGYAbPQNyHl9ZJqAArzIqM9v9PI4XFcOCmV13baLRnMwFjom5CytqSKyemJTA6DYRH7alF+JqXVjZTVNLpdiglhFvomZNQ3t7Kh7BgLI+wov9PCmd7bR6+xo30zABb6JmS8vaeG1naNmK6aXWWNHMbMccnWxGMGxELfhIw1O6tITYzjgpzwvMGaPxbmZ/LRoRMca/S4XYoJURb6JiS0tnfw5q5qrpqREdY3WOtNYV4mqvDGrmq3SzEhykLfhISNB45T39wWcb12upo5LpmxKQmstYFVTD9Z6JuQsHZnNXExUVyam+Z2Ka4S8d6A7e09tTS32g3YTN9Z6Jugp6qsKTnKJ6aMJjHeryEgwlphfianW9t5f98xt0sxIciv0BeRxSKyW0RKReSubt7/nojsFJFtIvK6iEzwea9dRLY4j5VdlzWmN3urGyk/fjpie+10ddHkVBLjom3sXNMvvYa+iEQDDwBLgHzgJhHJ7zLbZqBAVc8DngN+6fPeaVWd4zyWBqhuE0E6+6VHent+p/iYaC6fns7rJXYDNtN3/hzpLwBKVbVMVVuAFcAy3xlU9U1VbXJebgCyA1umiWRrdlZxfnYKmckJbpcSNArzMqmq97DjiN2AzfSNP6GfBZT7vK5wpvXk68ArPq8TRKRIRDaIyKe7W0BEbnXmKaqpqfGjJBMpqhua2VJeZ0f5XVw5PYMo8Y4gZkxf+BP63XWK7vZvShH5ElAA/Mpnco6qFgBfBH4tIlPOWpnqw6paoKoF6enpfpRkIsUbnTdYs/b8M4xKjKNgYiprSqy/vukbf0K/Ahjv8zobONJ1JhEpBO4Glqrqx5cLquoR598yYB0wdwD1mgiztqSKrJHDmDEmye1Sgs7CvExKKuupONHU+8zGOPwJ/Y1ArohMEpE4YDlwRi8cEZkLPIQ38Kt9po8SkXjneRrwCWBnoIo34e10Szvv7K2lMC8Dkci9CrcnnX/9vG5H+6YPeg19VW0D7gBWAyXAM6paLCL3iUhnb5xfASOAZ7t0zcwDikRkK/Am8HNVtdA3fnm3tBZPWwcL88e4XUpQmpSWyJT0RLs61/SJX1e6qOoqYFWXaff4PC/sYbn1wOyBFGgi1+slVSTFx7BgUqrbpQStwrxMHntvP/XNrSQnxLpdjgkBdkWuCUodHcrakmoum55OXIz9mPakMD+T1nbl7T3W6834x36bTFDaWlFHbaMnYgdM8dcFOaMYNTzWum4av1nom6D0ekk10VHCFdOtC++5REcJV83I5I1d1bS2d7hdjgkBFvomKK0tqWL+xFGMHB7ndilBb2F+BvXNbRQdOOF2KSYEWOiboFN+vIldRxvsKlw/XZqbTlx0lPXiMX6x0DdBpzO8rrbQ90tifAyXTB3N2pIqVO0GbObcLPRN0Hm9pJqpGSOYlJbodikhozAvk4PHmthX0+h2KSbIWeiboFLf3MqGsmPWtNNHV+dlALBmp12da87NQt8Elbd219DWoRQ6IWb8MzZlGLOzUqxd3/TKQt8ElbUlVaQmxjE3Z5TbpYScwrxMPjp0gtpGT+8zm4hloW+CRmt7B2/uquaqGRlER9kN1vqqMD8DVXhjlzXxmJ5Z6JugUXTgBPXNbdae30/5Y5MZl5JgV+eac7LQN0FjbUkVcTFRXJqb5nYpIUlEKMzP5J29tTS3trtdjglSFvomKKgqa0uquGTKaBLj/br5q+lGYV4mp1vbWb+v1u1STJCy0DdBYV9NIwePNVnTzgBdODmVEfEx1nXT9MhC3wSFzpC62rpqDkh8TDSXT0vn9ZIqOjrs6lxzNr9CX0QWi8huESkVkbu6ef97IrJTRLaJyOsiMsHnvVtEZK/zuCWQxZvwsbakitlZKYxNGeZ2KSGvMD+D6gYP2w+fdLsUE4R6DX0RiQYeAJYA+cBNIpLfZbbNQIGqngc8B/zSWTYVuBe4EFgA3Csi1gHbnOFYo4ePDp2wo/wAuXK6t8urXahluuPPkf4CoFRVy1S1BVgBLPOdQVXfVNUm5+UGINt5fg2wRlWPq+oJYA2wODClm3Dxxq5qVLH2/AAZOTyOggmjWGNdN003/An9LKDc53WFM60nXwde6cuyInKriBSJSFFNjQ37Fmle21nFuJQEZo5LdruUsLEwP5NdRxsoP97U+8wmovgT+t1dGtntGSIR+RJQAPyqL8uq6sOqWqCqBenpNlJSJGlqaePtPTUsmjkGEbsKN1A6b0v9ujXxmC78Cf0KYLzP62zgSNeZRKQQuBtYqqqevixrItdbu2vwtHWweNYYt0sJK5PSEpmaMYK1JdZ105zJn9DfCOSKyCQRiQOWAyt9ZxCRucBDeAPf96dsNbBIREY5J3AXOdOMAeDV4qOkJsYxf2Kq26WEncK8TDaUHaO+udXtUkwQ6TX0VbUNuANvWJcAz6hqsYjcJyJLndl+BYwAnhWRLSKy0ln2OPATvF8cG4H7nGnG0NLWwRsl1SzMy7QbrA2ChfkZtHUob+2282Tmb/y63l1VVwGruky7x+d54TmWfQx4rL8FmvC1fl8tDZ42a9oZJHPGj2J0YhxrdlbxqfPHuV2OCRJ2Ra5xzerio4xwxnc1gRcdJRTmZfLmrmo8bXYDNuNloW9c0d6hvFZcxZUzMoiPiXa7nLC1ePYYGjxtrC895nYpJkhY6BtXbDp4gmOnWlg805p2BtMlU0aTFB/DKzsq3S7FBAkLfeOKV3ccJS4miium23UZgyk+Jpqr8zJYs7OKtvYOt8sxQcBC3ww5VWV18VEuy02ze+cPgcWzxnCiqZUP91vHOWOhb1xQfKSew3WnucaadobE5dMySIiN4pUdR90uxQQBC30z5F7dcfTjniVm8A2Li+aKaRmsLj5q99g3Fvpm6L1afJQLJ6UyKjHO7VIixpLZY6hu8LC5/ITbpRiXWeibIVVa3UhpdaNdkDXErpyRQWy08Ko18UQ8C30zpFYXe0NnUb6F/lBKTojlk1PTeGXHUVStiSeSWeibIbVqeyVzxo9kTEqC26VEnCWzxlJx4rQNoxjhLPTNkNlfe4riI/Vcf95Yt0uJSItmZhITJby8zS7UimQW+mbIvLzVO5TCdRb6rhg5PI5Lc9P467ZKa+KJYBb6Zsi8vK2S+RNHMTZlmNulRKzrzxvH4brTfHSozu1SjEss9M2Q2FvVwO6qBq4/z27x66aFMzOJi47i5W02gF2k8iv0RWSxiOwWkVIRuaub9y8TkY9EpE1EPtflvXZnYJWPB1cxkeelbZVEibe/uHFPckIsl09PZ9X2SrtQK0L1GvoiEg08ACwB8oGbRCS/y2yHgK8CT3azitOqOsd5LO3mfRPmVJWXtx3hwkmjyUiyXjtuu/68sVTVe9h4wO7FE4n8OdJfAJSqapmqtgArgGW+M6jqAVXdBtht/MxZSiobKKs5xfXn2wncYFCYl0lCbJT14olQ/oR+FlDu87rCmeavBBEpEpENIvLpPlVnwsLL244QHSUsmWWhHwwS42O4akYGr+yotNstRyB/Qr+7Eav70hiYo6oFwBeBX4vIlLM2IHKr88VQVFNjgziHE2/TTiWXTBlNqt1rJ2hcf944ahtb+MButxxx/An9CmC8z+tswO9T/6p6xPm3DFgHzO1mnodVtUBVC9LTbVCNcLK5vI5Dx5v4lPXaCSpXzchgRHwMf9ly2O1SzBDzJ/Q3ArkiMklE4oDlgF+9cERklIjEO8/TgE8AO/tbrAk9L3x0mPiYKOu1E2QSYqNZPGsMr2w/SnOrDZoeSXoNfVVtA+4AVgMlwDOqWiwi94nIUgARmS8iFcDngYdEpNhZPA8oEpGtwJvAz1XVQj9CtLR18NK2IyzMzyQpIdbtckwXn5mbRYOnjTU7q9wuxQwhv8aqU9VVwKou0+7xeb4Rb7NP1+XWA7MHWKMJUet2V1PX1MpnLujLeX8zVC6aPJqxKQm8sPkwnzrfmt8ihV2RawbNC5sPMzoxjktz7TxNMIqKEpbNyeKtPTXUNnrcLscMEQt9MyhONrXyekk1nzp/HLHR9mMWrD5zQRbtHcpLW+22DJHCfhvNoPjr9kpa2ju4Ya417QSzaZlJ5I9N5oXN1osnUljom0HxwuYKJqcncl52itulmF585oIstlWcpLS60e1SzBCw0DcBV368iY0HTvBPuzi4AAATzklEQVSZuVmIdHdtnwkmS88fR5R4v6hN+LPQNwH354+8TQXL5ljTTijISE7g0tx0XvjoMO12582wZ6FvAqq9Q3mmqJxPTk1jfOpwt8sxfvpCwXiOnGzmnb12G5RwZ6FvAurd0loO153mxvnje5/ZBI3C/AxSE+NY8WF57zObkGahbwJqxYeHGDU8lkUzM90uxfRBfEw0n70gi7UlVdQ0WJ/9cGahbwKmpsHDmp1VfPaCbOJjot0ux/TRjfNzaOtQnv/ITuiGMwt9EzDPf1RBW4eyfIE17YSiqRkjWDAxlac3lqNqJ3TDlYW+CQhV5emN5RRMGMXUjCS3yzH9dOP88eyvPcWGMrvPfriy0DcB8cH+4+yvPcXyBTlul2IG4NrZY0lKiGHFxkNul2IGiYW+CYgVHx4iKSGG62bbkIihbFhcNDfMzeKVHUepa2pxuxwzCCz0zYDVNnpYtf0oN8zNYlicncANdTctyKGlrYNniqz7Zjiy0DcD9tQHh2hp7+ArF090uxQTAHljk1kwKZXfrz9oV+iGIb9CX0QWi8huESkVkbu6ef8yEflIRNpE5HNd3rtFRPY6j1sCVbgJDi1tHfxhw0Eum5bO1IwRbpdjAuTvPjGRw3WnbVStMNRr6ItINPAAsATIB24Skfwusx0Cvgo82WXZVOBe4EJgAXCviIwaeNkmWLyyo5LqBg9fu2Si26WYACrMyyRr5DAeX7/f7VJMgPlzpL8AKFXVMlVtAVYAy3xnUNUDqroN6Oiy7DXAGlU9rqongDXA4gDUbYLE4+sPMCktkcun2ehY4SQmOoqvXDyBDWXHKamsd7scE0D+hH4W4HtGp8KZ5g+/lhWRW0WkSESKamrshk+hYkt5HZsP1XHLxROIirJbKIebG+ePJyE2isffO+B2KSaA/An97n6b/T2749eyqvqwqhaoakF6uh0xhorH39vPiPgYPjsv2+1SzCAYOTyOG+Zm8+KWwxw/Zd03w4U/oV8B+F5Xnw34O6DmQJY1Qay6vpm/bq/k8wXZJCXEul2OGSRf+8REPG0dPPWhXawVLvwJ/Y1ArohMEpE4YDmw0s/1rwYWicgo5wTuImeaCXGPvruf9g7lFuumGdamZSZxaW4av3vvAM2t7W6XYwKg19BX1TbgDrxhXQI8o6rFInKfiCwFEJH5IlIBfB54SESKnWWPAz/B+8WxEbjPmWZCWF1TC3/ccJDrzxvHxLREt8sxg+xbV0ylttHDs3axVliI8WcmVV0FrOoy7R6f5xvxNt10t+xjwGMDqNEEmcfXH+BUSzvfunKK26WYIXDR5FQuyBnJb94qY/mCHGKj7ZrOUGb/e6ZPGj1t/O69AxTmZTJjTLLb5ZghICLccdVUDted5sXNh90uxwyQhb7pkyc/OMjJ063ccdVUt0sxQ+jK6RnkjU3mwbf22a0ZQpyFvvFbc2s7v31nP5+cmsac8SPdLscMIRHh9iunUFZzild3HHW7HDMAFvrGb89uqqCmwWNt+RFqyayxTE5L5IE3S21krRBmoW/8crqlnQfeKGXehFFcPHm02+UYF0RHCd+8Ygo7K+vtaD+EWegbvzy+/gBH65v5wTXTEbFbLkSqG+ZmkZsxgl+t3k1be9dbbZlQYKFvelXX1ML/rCvl6hkZXGhH+REtJjqKHyyeQVntKZ4pqnC7HNMPFvqmVw+u20ejp43/s3i626WYIFCYl0HBhFH8eu0eTrfYVbqhxkLfnNORutP8bv0BPjM32/rlG8Dbk+fOJTOobvDw2Ht2v/1QY6Fvzun+NXtA4R8X5rpdigki8yemUpiXwW/W7eOE3YEzpFjomx7tOHyS5z+q4MsXTyB71HC3yzFB5geLZ3CqpY371+5xuxTTBxb6plsdHcrdL+4gNTGO71xtR/nmbNMyk/jyRRP444aDbK846XY5xk8W+qZbT208xNbyOu6+Lo+UYXa/fNO97y2aTmpiPD96cbvdniFEWOibs9Q2evjlq7u5aHIqn57j78iYJhKlDIvln6/PY2vFSRtoJURY6Juz/NuqXZzytPHTT8+yC7FMr5aeP45Lpozml6/uoqbB43Y5phcW+uYMH5Qd4/mPKvjGZZOZmpHkdjkmBIgI9y2bxenWdv51VYnb5Zhe+BX6IrJYRHaLSKmI3NXN+/Ei8rTz/gciMtGZPlFETovIFufxm8CWbwKpobmV7z+3lfGpw/i23TrZ9MHUjBF884qpvLD5sN2XJ8j1GvoiEg08ACwB8oGbRCS/y2xfB06o6lTgfuAXPu/tU9U5zuO2ANVtBsGPX9rJ4ROn+fWNcxge59egasZ87NtXTWV2Vgo//PM2quub3S7H9MCfI/0FQKmqlqlqC7ACWNZlnmXA753nzwFXizUGh5RV2yt5blMFt185lXkTUt0ux4Sg2Ogo7r9xDqdb2/n+c9vs9stByp/QzwJ8R0SucKZ1O48zkPpJoPPOXJNEZLOIvCUil3a3ARG5VUSKRKSopqamTztgBq6qvpl/emE752WnWJ98MyBTM0Zw93X5vL2nhifeP+h2OaYb/oR+d0fsXb/Ce5qnEshR1bnA94AnReSsG7io6sOqWqCqBenp6X6UZAKlvUP5389spbm1nftvnGODXpsB+9KFOVw5PZ1/XVXCrqP1bpdjuvDnN7wCGO/zOhs40tM8IhIDpADHVdWjqscAVHUTsA+YNtCiTeD8/JUS3i2t5cdLZzIlfYTb5ZgwICL84nPnMXJ4LN94oojjdm+eoOJP6G8EckVkkojEAcuBlV3mWQnc4jz/HPCGqqqIpDsnghGRyUAuUBaY0s1APb+pgt++s5+vXDyBG+fnuF2OCSMZSQk89OUCquo93P6nj2i1AVeCRq+h77TR3wGsBkqAZ1S1WETuE5GlzmyPAqNFpBRvM05nt87LgG0ishXvCd7bVPV4oHfC9N2W8jp++MJ2Lpqcyj9f37UzljEDN2f8SP7thtm8X3aMn7680+1yjMOvfnmqugpY1WXaPT7Pm4HPd7Pc88DzA6zRBFjlydP8rz8UkZEUz//cPM/a8c2g+ey8bEoq63nk3f1MG5PEzRdOcLukiGedsSNMTYOHm3/7Aac87Tx72wJSE+PcLsmEubuWzKC0ppEfvbiDpIRYlp4/zu2SIpod4kWQE6da+PKjH1B5spnffW0+eWNtJCwz+GKio3jw5nnMn5jKPz69hdeK7YpdN1noR4j65la+8tiHlNWe4pFbCpg/0S7AMkNnWFw0j311PrOzUrjjyc28tceux3GLhX4EqGnw8KVHPmDX0Xp+86UL+MTUNLdLMhFoRHwMv//aAqZmjODWJ4pYbUf8rrDQD3P7ahr5zIPvsbeqkd98aR5Xzch0uyQTwVKGx/LHv7+QvLHJ3PbHTfx+/QG3S4o4FvphrOjAcT774HqaPO2suPUirs6zwDfuS02M46lvXERhXib3rizmX1eV0GGjbg0ZC/0wpKr8fv0BvvjIB4waHsefv3UJ548f6XZZxnxsWFw0v/nSPL5y8QQefruMr/9+I8cabQCWoWChH2bqmlq49Q+buHdlMZ+YMprnv3kJE0Ynul2WMWeJjhJ+vHQm9y2byXulx1jyH++wfl+t22WFPQv9MPLO3hqu/Y93WLe7mh9dl8ejt8y3fvgmqIkIX7l4Ii/cfgkjEmK4+ZEP+MWru2hubXe7tLBloR8GKk+e5vY/fcSXH/2Q+Nhonv/mJfz9pZOJirIhDUxomDkuhZe//Uk+Py+bB9ftY+H9b/F6SZXbZYUlCbaBDgoKCrSoqMjtMkJCU0sbj68/wH+/UUp7h3L7lVO59bLJJMRGu12aMf22fl8t9/ylmNLqRgrzMrlz8XRyM2285t6IyCZVLeh1Pgv90NPoaeMP7x/kt++UcfxUC4V5Gdxz/UxyRg93uzRjAqKlrYPfvbef/3x9L02t7Vw7ayzfvnoqM8bYVeQ9sdAPQ2U1jazYWM4zReXUNbVy+bR0vnN1LvMmjHK7NGMGxfFTLTz6bhm/X3+QRk8bV83I4IsLcrhiejoxdqPAM1joh4m6phbWllTz/KYK3i87RkyUUJiXyW1XTGGOdcM0EaKuqYXfvXeAJz88RE2DhzHJCXy+IJtrZ49lxpgkbEhuC/2Qparsq2nkvdJjvLbzKBvKjtPeoYxPHcby+Tl8fl42GckJbpdpjCta2zt4vaSapz48xNt7a1CFCaOHc83MMVw+LZ15E0ZF7DktC/0QccrTxs7KerZXnGRzeR0byo5R0+C9SGVyWiKLZ43hmpljOC87xY5mjPFR0+Bhzc4qXi0+yvrSWto6lLjoKObkjGTBxFRmZaUwOzuFcSkJEfG7E9DQF5HFwH8A0cAjqvrzLu/HA08A84BjwI2qesB574fA14F24Duquvpc2wrH0D/d0k5VfTOHjjdx4Ngp9td6HwdqT3HweBOd/wVjkhO4cHIqF08ezUWTRzNh9PCI+GE1ZqAamlspOnCC98uO8f6+Y+ysrKfdubXDqOGxTEkfwcS0RCalJTJxdCIT04YzNmUYI4fFhk3X5oCFvjPG7R5gId4B0DcCN6nqTp95vgWcp6q3ichy4AZVvVFE8oGngAXAOGAtME1Ve7zyIlhCv6ND8bR10NLWgaetHc8Z/3bgae3glKeNBk8rDc1t1J92/m1uo765ldoGDzUNHqobPDR62s5Y9/C4aCaMTmRS2nCmZSYxOyuF2Vkp1mxjTIA0t7azs7Ke4sMnKT5ST5lzkFXdcOatHmKjhbQR8WQkxZOelMDI4bEkJcSQlBBLckIMyQl/e50QG0VcTBTxMdHEx0QRH+t97p0WRUyUuHqQ5m/o+zNy1gKgVFXLnBWvAJYBvoNeLgP+xXn+HPDf4t37ZcAKVfUA+50xdBcA7/u7I/46caqFzz/0Ph2qqEJ7h571vEO9bebtqnR0eN/r6HztvNfhzN8fcTFRJDs/IGkj4sgbm8xl0+LJSI4nIymB7FHDmJSWSEZSvB3BGzOIEmKjuSBnFBfknNmz7ZSnjQPHTnHwWBNV9c1UN3iorvdQ3dBMxYkmio94D966Hqj5K0q8t5eIEvn4385pf3vtPI+CaPF+UQiAQP7YZP77ixcM/AM4B39CPwso93ldAVzY0zyq2iYiJ4HRzvQNXZbN6roBEbkVuBUgJyfH39rPEBMtTM9MQoSPP+gonw89SsR5zccffI/zOdO6fqPHO9/occ5jRLw34L1HAjHEx0TmCSRjQkVifAwzx6Uwc1zKOedr71AaPW00NP/tS8DT2uWv/taOM563dngPJjsPKjtUae9wDiyd6dr53Ge64j3gVGDiENwny5/Q7+6QtOuhcE/z+LMsqvow8DB4m3f8qOksSQmxPHDz4H5DGmMiQ3SUkDIslpRhsW6XEnD+XN1QAYz3eZ0NHOlpHhGJAVKA434ua4wxZoj4E/obgVwRmSQiccByYGWXeVYCtzjPPwe8od4zxCuB5SISLyKTgFzgw8CUbowxpq96bd5x2ujvAFbj7bL5mKoWi8h9QJGqrgQeBf7gnKg9jveLAWe+Z/Ce9G0Dbj9Xzx1jjDGDyy7OMsaYMOBvl027Y5ExxkQQC31jjIkgFvrGGBNBLPSNMSaCBN2JXBGpAQ4OYBVpQG2AynGb7Utwsn0JTpG+LxNUNb23mYIu9AdKRIr8OYMdCmxfgpPtS3CyffGPNe8YY0wEsdA3xpgIEo6h/7DbBQSQ7Utwsn0JTrYvfgi7Nn1jjDE9C8cjfWOMMT2w0DfGmAgSdqEvIj8RkW0iskVEXhORcW7X1F8i8isR2eXszwsiMtLtmvpLRD4vIsUi0iEiIdmtTkQWi8huESkVkbvcrmcgROQxEakWkR1u1zIQIjJeRN4UkRLn5+u7btfUXyKSICIfishWZ19+PCjbCbc2fRFJVtV65/l3gHxVvc3lsvpFRBbhHZugTUR+AaCqd7pcVr+ISB7QATwEfF9VQ+pWqiISDewBFuIdHGgjcJOq7jzngkFKRC4DGoEnVHWW2/X0l4iMBcaq6kcikgRsAj4div8vzrjiiaraKCKxwLvAd1V1Qy+L9knYHel3Br4jkW6GZwwVqvqaqnaO0LwB78hjIUlVS1R1t9t1DMACoFRVy1S1BVgBLHO5pn5T1bfxjn0R0lS1UlU/cp43ACV0Mw53KFCvRudlrPMIeH6FXegDiMjPRKQcuBm4x+16AuTvgFfcLiKCZQHlPq8rCNFwCVciMhGYC3zgbiX9JyLRIrIFqAbWqGrA9yUkQ19E1orIjm4eywBU9W5VHQ/8CbjD3WrPrbd9cea5G+/IY39yr9Le+bMvIUy6mRayf0WGGxEZATwP/EOXv/ZDiqq2q+ocvH/VLxCRgDe99TpcYjBS1UI/Z30S+Ctw7yCWMyC97YuI3AJcD1ytQX4Cpg//L6GoAhjv8zobOOJSLcaH0/79PPAnVf2z2/UEgqrWicg6YDEQ0JPtIXmkfy4ikuvzcimwy61aBkpEFgN3AktVtcnteiLcRiBXRCaJSBzecaBXulxTxHNOfj4KlKjqv7tdz0CISHpnDz0RGQYUMgj5FY69d54HpuPtKXIQuE1VD7tbVf84A83HA8ecSRtCuCfSDcB/AelAHbBFVa9xt6q+EZFrgV8D0cBjqvozl0vqNxF5CrgC7y18q4B7VfVRV4vqBxH5JPAOsB3v7zzAP6nqKveq6h8ROQ/4Pd6fryjgGVW9L+DbCbfQN8YY07Owa94xxhjTMwt9Y4yJIBb6xhgTQSz0jTEmgljoG2NMBLHQN8aYCGKhb4wxEeT/A7X7pe+HM1uuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def normal(x, mu, sigma): # mu, sigma are parameters\n",
    "    return (1/(2*np.pi*sigma**2)**0.5)*np.exp(-1*((x - mu) / sigma)**2)\n",
    "\n",
    "# CHANGE THESE TO SEE DIFFERENT POISSON DISTRIBUTIONS\n",
    "mu = 0\n",
    "sigma = 1\n",
    "\n",
    "x = np.arange(-3*sigma, 3*sigma, sigma / 30)\n",
    "y = [normal(i, mu, sigma) for i in x]\n",
    "plt.plot(x, y)\n",
    "plt.title(f\"Normal Distribution ($\\mu$ = {mu}, $\\sigma$ = {sigma})\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Central Limit Theorem\n",
    "\n",
    "Why is the normal distribution so common? The *central limit theorem* is one reason:\n",
    "\n",
    "If $X_1, X_2, \\dots, X_N$ are independent, identically distributed (IID) random variables, with mean $\\mu$ and standard deviation $\\sigma$, then in the limit as $N \\to \\infty$, the random variable \n",
    "$$Z = \\frac{1}{N}\\sum_{i=1}^{N}X_i$$\n",
    "is normally distributed, with mean $\\mu$ and standard deviation $\\sigma / \\sqrt{n}$\n",
    "\n",
    "Note: the $X_i$ can have **any** distribution! You don't even have to know the distribution itself, just its mean and standard deviation. This means that randomly sampling from a population distribution allows us to estimate the mean of the population, with the accuracy of that estimate (i.e. the standard deviation) decreasing as the sample size increases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistics\n",
    "\n",
    "> \"There are three types of lies -- lies, damn lies, and statistics.\" -Benjamin Disraeli\n",
    "\n",
    "If probability theory is how we model random systems, statistics is how we evaluate those model. Statistical tools tell us the extent to which our system matches the model; they let us estimate the parameters of a model from observations; and they let us predict the distribution of future samples from a model. Instead of throwing a list of statistical tools at you, I want to explore the concepts underlying these uses of statistical tools."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The P-Value\n",
    "> \"When you run two experiments, one testing a clinical treatment that seems theoretically promising and the other testing whether dead salmon respond emotionally to romantic photos, and both experiments succeed with p-values of .03, you don't really want to treat the two hypotheses the same.\" -Jordan Ellenberg, *How Not to Be Wrong: The Power of Mathematical Thinking*\n",
    "\n",
    "P-values are often misunderstood by even trained scientists. P-values are, in fact, a tool from probability theory; not statistics. But in statistics, we use this tool to make a statistical conclusion about a model for a system. \n",
    "\n",
    "Consider a population of people in a country $A$, with heights that are normally distributed about an unknown mean $M$, but with a known standard deviation $\\sigma$. In neighboring country $B$, the heights are normally distributed with mean $\\mu$ and standard deviation $\\sigma$. Is this $\\mu$ the value of the mean in country $A$?\n",
    "\n",
    "First, we set out to collect a simple, random sample of $N$ members of the population in $A$. We can consider the height of each sampled person to be a random variable $X_i$, with distribution identical to the population distribution. The average of these sampled variables, $Z = \\sum_iX_i/N$, is a random variable with mean $M$ and standard deviation $\\sigma /\\sqrt{N}$. When we conducted the random sample, we essentially sampled from $Z$; the sample mean $\\mu_{sample}$ was drawn from the distribution of $Z$. And, if the mean of country $A$ was equal to that of country $B$, the probability of getting a value as large or larger than $\\mu_{sample}$ would be $$P(x \\geq \\mu_{sample}) = 1 - \\Phi\\left(\\frac{x-\\mu}{\\sigma/\\sqrt{n}}\\right) = \\int_{\\mu_{sample}}^{\\infty}\\frac{1}{\\sqrt{2\\pi\\sigma^2/n}}e^{-((x-\\mu)/(\\sigma/\\sqrt{n}))^2}dx$$\n",
    "\n",
    "The function $\\Phi$ is known as the *cumulative distribution function* for the normal distribution; it represents the probability of getting a value less than $x$, for a standard normal distribution with mean $0$ and variance $1$.\n",
    "\n",
    "The p-value is the probability of getting a value as extreme or more extreme than the observed value, given our hypothesis that $M = \\mu$. In statistical terminology, $M = \\mu$ is the *null hypothesis*; we are comparing it to the alternative hypothesis $M > \\mu$. \n",
    "\n",
    "How do we judge the validity of the null hypothesis? If this p-value is small, then the observed event is unlikely given that our model is true. We can conclude that the model must be wrong, and so we reject the null hypothesis and accept the alternative hypothesis. \n",
    "\n",
    "What if the p-value is large? Then the event that occured was highly likely, given our null hypothesis. This evidence doesn't *confirm* the null hypothesis; in fact, nothing can confirm a null hypothesis! We simply fail to reject it based on the observed evidence.\n",
    "\n",
    "What is a small p-value? You may have seen reference to a significance level of $\\alpha = 0.05$ before. For a statistical inference about a set of hypotheses to be valid, you must choose a significance level before conducting any sort of test. Researchers across many scientific communities have standardized upon $0.05$ as a reasonably small p-value; if you observe a value that should only occur one out of every twenty experiments, then your null hypothesis is likely to be wrong. $0.05$ is a somewhat arbitrary choice; many other significance levels are also used. What's important is to know what likeliness threshold is reasonable for the data you are collecting: for example, if we are running an experiment 20 times, we would expect at least one to be significant at the $\\alpha=0.05$ level about $64\\%$ of the time. \n",
    "\n",
    "Choosing a p-value to fit the observed data, or deliberately choosing too high a value for $\\alpha$, is known as p-hacking. It's unethical, and leads to improper statistical conclusions. To avoid p-hacking, choose your significance level before you run any experiments (which is why a standard of $0.05$ is useful!), and make sure that you adjust that p-value to account for multiple experiments.\n",
    "\n",
    "## Statistical Tests and Test Statistics\n",
    "> \"At school you were taught about chemicals in test tubes, equations to describe motion, and maybe something on photosynthesis—about which more later—but in all likelihood you were taught nothing about death, risk, statistics, and the science of what will kill or cure you.\" -Ben Goldacre, *Bad Science*\n",
    "\n",
    "When testing a hypothesis, you must identify the statistical test to be used, perform the experiment, run the test, and form a conclusion based on the p-value. In high-school, I learned an acronym to remember the steps for hypothesis testing: PHANTOM\n",
    "\n",
    "P: identify the **parameters** of interest for your test (e.g. population mean)\n",
    "\n",
    "H: state your null and alternative **hypotheses**- make sure these cover every possibility!\n",
    "\n",
    "A: identify the **assumptions** you are making about the population to be sampled\n",
    "\n",
    "N: choose and **name** the statistical test to be performed (1-sided Z-test, 2-sided t-test, etc.)\n",
    "\n",
    "T: calculate a **test** statistic, using data sampled from the population of interest\n",
    "\n",
    "O: **obtain** a p-value using the cumulative distribution function of the distribution of the test statistic\n",
    "\n",
    "M: **make a decision** about whether or not to reject your null hypothesis\n",
    "\n",
    "This procedure is designed for scoring points on the AP Statistics exam. You won't always explicitly state the null and alternative hypotheses in a research paper. In fact, the entire hypothesis testing procedure might be written in a single sentence! As long as the important details- the parameter of interest, the significance level, the statistical test, the value of the test statistic, and the p-value- are included in your work, your audience will understand your statistical analysis.\n",
    "\n",
    "What exactly are statistical tests, and how are they related to test statistics? *Statistical tests* are a systematic way of relating hypotheses to p-values, by linking specific hypothesis (under specific assumptions about the sample and parameter) to methods for calculating a *test statistic*. This test statistic, under those assumptions, has a fixed probability distribution, and we can easily look up or calculate the p-value for the calculated test statistic.\n",
    "\n",
    "In our example above, we performed a *one-sided Z-test*. The one-sided Z-test:\n",
    "\n",
    ">links the hypotheses \"The mean of the population is $\\mu$\" and \"The mean of the population is greater than $\\mu$\"; \n",
    "\n",
    ">under the assumption that we calculated a mean $M$ from a simple random sample of size $N$ of a population with variance $\\sigma^2$;\n",
    "\n",
    ">to the $Z$-score $(x - \\mu) / (\\sigma/\\sqrt{N})$. \n",
    "\n",
    "From the central limit theorem this $Z$-score has a standard normal distribution, and so we can compute our p-value easily. \n",
    "\n",
    "Choosing a statistical test can be difficult, but the task is made easier by identifying the population parameter of interest and the hypotheses you wish to test about that parameter. Once we have a parameter and hypotheses, the assumptions we can make about our sample typically inform our choice of test. For example, if instead of having a known standard deviation $\\sigma$, the distribution of heights in country $A$ has an unknown standard deviation, we would use a *t-test*, with a test statistic that follows *Student's t-distribution*.\n",
    "\n",
    "If for some reason you don't want to calculate a p-value, but you wish to perform a hypothesis test, you can use a table listing *critical values* for the test statistic at a series of significance levels. If your test statistic lies beyond the critical value at your significance level, its p-value is smaller than the significance level."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence Intervals\n",
    "> \"The problem with the world is that the intelligent people are full of doubts, while the stupid ones are full of confidence.\" -Charles Bukowski\n",
    "\n",
    "*Confidence intervals*, like p-values, are a tool from probability used to make statistical conclusions. Consider, again, the sample of heights from the population of country $A$. Given our sample mean $\\mu_{sample}$, what can we say about the likely value of the true population mean? Our sample mean is, by the central limit theorem, a sample from a normal distribution with mean $M$ and standard deviation $\\sigma/\\sqrt{N}$. If we repeated this \"sampling of the sample mean\" many times, we would get many different samples from this same normal distribution. For each sample, we can construct  a $(1-\\alpha)\\times 100 \\%$-confidence interval; the fraction of these intervals that contain the population mean $M$ is $1-\\alpha$, in the limit as the number of samples goes to infinity. \n",
    "\n",
    "Note: a confidence interval does **not** represent a probability that the population parameter is within some range of the estimate! The parameter is unknown, but guaranteed to be either inside or outside that range. The confidence level represents the *fraction* of confidence intervals that contain the population parameter, after sampling many times from the distribution of interest.\n",
    "\n",
    "How do we construct an $(1-\\alpha)\\times 100\\%$-confidence interval? For a cumulative distribution function $\\Phi$ of the distribution of the test statistic, we find the critical value $z^* = \\Phi^{-1}\\left(1-\\alpha/2\\right)$. If our parameter estimate is $\\bar{x}$, and the standard deviation of the distribution of that estimate (also known as the *standard error*) is $SE$, the $(1-\\alpha)\\times 100\\%$-confidence interval is:\n",
    "\n",
    "$$\\left(\\bar{x} - z^*\\cdot SE, \\bar{x} + z^* \\cdot SE\\right)$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Regression\n",
    "> \"History fancies itself linear - but yields to a cyclical temptation.\" Criss Jami, *Healology*\n",
    "\n",
    "Consider a sample of two variables, $X$ and $Y$, drawn from a population. How can we characterize the relationship between these variables, given the observed samples? *Linear regression* gives us a powerful set of tools for characterizing simple relationships, by finding the line that best fits the data. Linear regression finds the values $m$ and $b$ that minimize the squared error of approximating the relationship between $X$ and $Y$ as linear:\n",
    "$$\\text{Error }(m, b) = \\sum_i\\left[y_i-\\left(mx_i+b\\right)\\right]^2$$\n",
    "To minimize this error, we take its derivatives with respect to $b$ and $m$, and set them equal to zero:\n",
    "$$\\begin{align*}\n",
    "\\frac{dE}{db} = -2\\sum_i(\\left[y_i-\\left(mx_i+b\\right)\\right] = 0 &\\to b = \\bar{y} - m\\bar{x} \\\\\n",
    "\\frac{dE}{dm} = -2\\sum_i(x_i-\\bar{x})\\left[y_i-\\left(mx_i+b\\right)\\right] = -2\\sum_i(x_i-\\bar{x})\\left[y_i-\\left(mx_i+\\bar{y} - m\\bar{x}\\right)\\right] &\\to m = \\frac{\\sum_i(x_i-\\bar{x})(y_i-\\bar{y})}{\\sum_i(x_i-\\bar{x})^2}\n",
    "\\end{align*}$$\n",
    "Note that $\\bar{x}$ is the average of the $X$ samples.\n",
    "How strong is this relationship? The *sample correlation coefficient* can tell us. If we let $s_x$ and $s_y$ represent the standard deviations of the $N$ samples, then the correlation coefficient $r$ is:\n",
    "\n",
    "$$r = \\frac{1}{N-1}\\sum_i\\left[\\frac{x_i-\\bar{x}}{s_x}\\cdot\\frac{y_i-\\bar{y}}{s_y}\\right]$$\n",
    "\n",
    "Of course, there are statistical tests for the accuracy of all of these quantities. \n",
    "\n",
    "I want to emphasize a few points about linear regression. \n",
    "\n",
    "First, the equation $y=mx+b$ only explains the data within the range of $X$ values sampled. If a value lies outside that range, there are no guarantees that the regression will give a valid prediction of the $Y$-value. \n",
    "\n",
    "Second, you can get a better regression fit by removing outliers from your data. An outlier can be a value with too extreme an $X$-value, too extreme a $Y$-value, or both. In general, it is a good idea to inspect your data visually (i.e. with a histogram, a box-and-whiskers plot, or a scatterplot) before performing a regression or any other statistical procedure.\n",
    "\n",
    "Third, not all relationships are linear. But many relationships are linear for small variations in the independent variable. Think of this as a statistical analogue to a taylor series expansion: If $Y = f(X)$, then we can expand $f$ in a taylor series about $\\mu_X$, and truncate after the linear terms: $$Y \\approx f(\\mu_X) + (X - \\mu_X)f'(\\mu_X)$$\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}