Discrete Distributions.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true
   },
   "outputs": [],
   "source": [
    "%%javascript\n",
    "MathJax.Hub.Config({\n",
    "    TeX: { equationNumbers: { autoNumber: \"AMS\" } }\n",
    "});"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discrete distributions\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Bernoulli distribution\n",
    "Experiment with binary outcome. Random variable $X \\in \\{0,1\\}$ with\n",
    "    $$ P[X=1] = p   ;  P[X=0] = 1-p $$\n",
    "    \n",
    "The **probability mass function** (*pmf*) is\n",
    "\\begin{equation}\n",
    " f_X(k;p) = p^k (1-p)^{(1-k)} \\hspace{2cm} k \\in \\{0,1\\}\n",
    "\\end{equation}\n",
    "\n",
    "Expectation and variance of $X$:\n",
    "\\begin{equation}\n",
    "\\mu = p \\hspace{2cm} \\sigma^2 = p(1-p)\n",
    "\\end{equation}\n",
    "    \n",
    "\n",
    "### Binomial distribution\n",
    "$n$ repetitions of an experiment with binary outcomes, with $p$ the probability of success. Random variable $X$ is the number of successes, and the *pmf* is given by:\n",
    "\\begin{equation}\n",
    " f_X (k; n, p) = P[X=k] = \\binom{n}{k} p^k (1-p)^{(n-k)}  = \\frac{n!}{k!(n-k)!)} p^k (1-p)^{(n-k)}\n",
    "\\end{equation}\n",
    "\n",
    "Since each repetition is an independent Bernoulli trial, the Binomial random variable $X$ can be thought of as a sum of $n$ Bernoulli random variables $$ X = X_1 + X_2 + \\cdots + X_n $$\n",
    "\n",
    "This provides an elegant way to derive the expectation and  variance of $X$:\n",
    "\\begin{equation}\n",
    " \\mu_X = E(X) = \\sum_{i=1}^{n} E(X_i) = np \\hspace{2cm} \\sigma^2_X = \\sum_{i=1}^{n} \\sigma^2(X_i) = n p (1-p)\n",
    "\\end{equation}\n",
    "\n",
    "The moment generating function, with $q=1-p$, is\n",
    "$$\n",
    "M_X(t) = E(e^{tX}) = \\sum_{k=0}^n\\binom{n}{k}\\;e^{tk} p^k q^{n-k} = \\sum_{k=0}^n\\binom{n}{k}\\;(pe^t)^k q^{n-k} = (q+ pe^t)^n\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import statsmodels.api as sm\n",
    "\n",
    "import matplotlib\n",
    "matplotlib.use('nbagg')\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEICAYAAACwDehOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucFOWd7/HPV0YwBpWoGA1CkKAG8MLGETQxOicm3tagvkTFywp7vAQvcWOMq54kxph1N8Zd3WPQg2Y13uI9MZKIGnfVzSaryOgSFZCAggohERE0RlHB3/mjninbpofpudR0z8z3/XrNa6qrnqr6VXc//at66vIoIjAzMwPYqNYBmJlZ/XBSMDOznJOCmZnlnBTMzCznpGBmZjknBTMzyzkpdCNJj0o6OQ1PkfSbkmlvShrRiWX/H0n/1hVxlscj6QZJ/9CFy54u6dtdtTzrvFp9Ju1Zb2n96eC6Pi9pQUfnr7C8+yVNTsMfqs9dsOzjJf2qq5bXHg21WGlXkrQE+DiwDngP+G9gakS8XMu42isiBrY2TdKjwF5k2xfAQuAu4IqIeCfN/4/VrCct65aI2GAC2VA87SFpCnByROxTsuypXbFsq15b9aRWn0lXrVfSRcA3gTVp1HLgV8AlEbE8reu/gJ2rXNbIiDhhQ+Ui4uBOhFy6vuHAYmDjiFiblv0T4Cddsfz26i1HCl9OP2LbAX8CftiRhUiq5yR5ZkRsRraN5wCTgJmS1JUrqfP3wDqnS+pJHbsj1ZEtgSOAbYEnJW3XlStRprf8dq6nV21YRKwB7gZGt4yTNEDSP0t6SdKf0uHqR9K0JklLJZ0n6Y/Aj0vGnSPpFUnLJf1tyfK2kHSTpBWSXpT0rZYviKSLJN1SUna4pKjmhzaVG1nFNv4lIh4FJgB7A39dvm5Jm0i6RdJKSaslzZb0cUmXAJ8HpqXmoWkl6z5D0kKyo5BK8Wwt6SFJf5b0n5I+2do2thzmSxoFTAf2TutbnaZ/qDlK0imSFkl6TdIMSZ8oe1+mSlqYtuWqrk6EfU0r9ST/TDpZB6ZI+q2kK9Ln9YKkz6bxL6flTW5lvR+T9Mu03FVpePsObN97ETEXOAZYQbYTlW9XybrPk7QsfacXSNpf0kHA/wGOSd/Z36Wyj0q6RNJvgbeAEVq/OUuSpkl6XdJzkvYvmbBE0hdLXpf+Vvw6/V+d1rm31m9e/myqx6+n/58tmfaopO+l9/3Pkn4laev2vm8telVSkLQp2Rfh8ZLR3wd2AsYCI4EhwIUl07cl27P4JHBqybgtUtmTgKskfSxN+2GaNgLYDzgR+Fu6WUS8BDST/ciXm0wW41BgK2Aq8HZEfBP4L7KjjoERcWbJPIcD4yn5oShzPPA9YGtgDlUc2kbE/LTux9L6BpWXkfQF4J+Ao8n2YF8Ebi8rdiiwJ7BbKndgW+u21rVST8p1pg6MB54m++7dSvZ57klW/04g2ymp1Dy5EfBjsro4DHgbmNb+LcxExDrgXirUEUk7A2cCe6ajiwOBJRHxAPCPZEcdAyNi95LZ/obsN2Izsu9pufHA82R15DvAzyRtWUWo+6b/g9I6HyuLdUvgPuBKsvf0cuA+SVuVFDuO7DPYBugPfKOK9VbUW5LCz9Ne6OvAl4DLIEvbZB/i2RHxWkT8mewDn1Qy7/vAdyLinYh4O417D7g47XHMBN4EdpbUL817QUT8OSKWAP9C9mWphT+QJbRy75F9eUZGxLqIeDIi3mhjWf+U3qO3W5l+X0T8Op3D+CbZ3v/QjoeeOx64PiKeSsu+IC17eEmZ70fE6pQIHyFL8NZ+FetJKzpTBxZHxI/Tj/IdZDsnF6c69ivgXbIE8SERsTIifhoRb6W6eglZ0umM1urIOmAAMFrSxhGxJCKeb2NZN0TE3IhYGxHvVZj+CvCv6T27A1hAOpLvpL8GFkbEzWndtwHPAV8uKfPjiPh9qr930ok60luSwuFpL3QTsuz/n5K2BQYDm5K1K65OFeKBNL7FinQ4XWplywmf5C1gINkewMZ8eC/hRbK9qVoYArxWYfzNwIPA7ZL+IOkHkjZuY1ltnZjPp0fEm2m9n2i9eNU+Qcn7mZa9kg+/p38sGW75LKz9WqsnlXSmDvypZPhtgIgoH7feZyhpU0nXpCapN8iaVQalRNRRFetIRCwCvgZcBLwi6fbSZstWtFVHlsWHnzD6IgXUkZJlF1JHektSALLDxYj4GdlewD7Aq2RfwDERMSj9bVF2ZU17HhP7Ktke1CdLxg0DlqXhv5AloRatVbhOS3vpe5A1B31I2lP5bkSMBj5L1vxyYsvkVhbZ1vuQHxWkQ/8tyfbC/pJGt7bdbS33D5S8n5I+SnaUs6zVOaxTKtST9mirDnTGOWRXB42PiM35oFmlQ+eQ0nmOL1OhjgBExK3pqrhPkn1PL22Z1Moi2/ouDyk73zWM7PsNG/5taFcdKVl2IXWkVyUFZQ4DPgbMj4j3gR8BV0jaJpUZIqlDbdLpcPhO4BJJmyk72fp1oOWE0RxgX0nDJG1B1hTSpdLe1H5kbaVPADMrlPlfknZNe1hvkFXi99PkP5G1BbfXIZL2kdSf7NzC4xHxckSsIPtyniCpn6T/DXyqZL4/Adun+Sq5DfhbSWMlDSBr3puVmiWsAOX1pD3zVlEHOmMzsp241akd/TsdWYikBmUXOdxG9uN7eYUyO0v6QvrOrUnrLa0jw9X+K4y2Ac6StLGko4BRfFA/5wCT0rRGYGLJfCvSulurlzOBnSQdl7btGLJzf79sZ3xV6S1J4ReS3iT7AbwEmJyuPgA4D1gEPJ4OSf+dKq5V3oCvkmX9F4DfkJ1Iux4gIh4ia0N9GniSrv3Qpkn6M9kX9l+BnwIHpcRXbluyq0veIKv0/0nWpATwf4GJyq7uuLId67+VrJK+RnaEUnoN9ynAuWTNPmPIroFv8TAwF/ijpFfLFxoR/w58O23PcrKEMqm8nHWJDdWT9mi1DnTSvwIfITsaeZysqbc9jknb9zowg+z7uEdE/KFC2QFkF6G8Stb0sg0f7MTdlf6vlPRUO9Y/C9gxLfMSYGJErEzTvk323V4FfJfsPQMgIt5K5X+bmrn3Kl1oWsahZEdSK4G/Bw6NiPXqU1dQuJMdMzNLesuRgpmZdQEnBTMzyzkpmJlZzknBzMxyPe7hZ1tvvXUMHz681mFYnXnyySdfjYjBbZfs/VxHrJJq60iPSwrDhw+nubm51mFYnZFU6Vk0fZLriFVSbR1x85GZmeWcFMzMLOekYGZmOScFMzPLOSmYmVmu0KQg6SBl3dwtknR+K2WOljRP0lxJt1YqY2Zm3aOwS1LTY5uvIuvhaSkwW9KMiJhXUmZHsicTfi4iVrU83trMzGqjyCOFccCiiHghIt4l66f1sLIypwBXRcQqgIh4pcB4zMysDUUmhSF8uPu6pazfbeVOZJ1H/FbS45IOqrQgSadKapbUvGLFivWmNzU10dTU1EVhm1k9cj3vHrU+0dxA1ilFE3As8CNJg8oLRcS1EdEYEY2DB/tJBmZmRSkyKSyjpF9fYHvW71N0KTAj9Sm8GPg9WZIwM7MaKDIpzAZ2lLRD6p93ElkXeaV+TnaUgKStyZqTXigwJjMz24DCkkJErAXOBB4k6yf4zoiYK+liSRNSsQfJ+kGdBzwCnFvSp6mZmXWzQp+SGhEzgZll4y4sGQ7g6+nPzMxqrNYnms3MrI44KZgVrK07+yXtK+kpSWslTawwfXNJSyVN656IrS9zUjArUMmd/QcDo4FjJY0uK/YSMAVo7TEv3wN+XVSMZqWcFMyK1ead/RGxJCKeBt4vn1nSHsDHgV91R7BmTgpmxarmzv6KJG0E/AvwjSrKbvCuf7NqOSmY1a/TgZkRsbStgr3xrn8/1qI2Cr0k1cyqurO/NXsDn5d0OjAQ6C/pzYio+Bh6s67gpGBWrPzOfrJkMAk4rpoZI+L4lmFJU4BGJwQrmpuPzApUzZ39kvaUtBQ4CrhG0tzaRWx9nY8UzApWxZ39s8malTa0jBuAGwoIr8dqOd/w6KOP1jSO3sZHCmZmlnNSMDOznJOCmZnlnBTMzCznpGBmZjknBTMzyzkpmJlZzknBzMxyTgpmZpZzUjAzs5yTgpmZ5ZwU+qiI4KyzzmLkyJHstttuPPXUUxssP2HCBHbZZZduis6s9p544gnGjh3L2LFj2X333bnnnnvyaQ888AA777wzI0eO5Pvf/36hcTz55JPsuuuujBw5krPOOouIaLXs7NmzaWho4O677+7w+pwU+qj777+fhQsXsnDhQq699lpOO+20Vsv+7Gc/Y+DAgd0YnVnt7bLLLjQ3NzNnzhweeOABvvKVr7B27VrWrVvHGWecwf3338+8efO47bbbmDdvXmFxnHbaafzoRz/K6+sDDzxQsdy6des477zzOOCAAzq1PieFHm7JkiV8+tOf5vjjj2fUqFFMnDiRt956q8357r33Xk488UQksddee7F69WqWL1++Xrk333yTyy+/nG9961tFhG9WqCVLljBq1ChOOeUUxowZwwEHHMDbb79d1bybbropDQ3Zg6TXrFmDJCA7ghg5ciQjRoygf//+TJo0iXvvvXeDy5oyZQpTp06lsbGRnXbaiV/+8pdVxbB8+XLeeOMN9tprLyRx4okn8vOf/7xi2R/+8IcceeSRbLPNNlUtuzWFJgVJB0laIGmRpPU6B5E0RdIKSXPS38lFxtNbLViwgNNPP5358+ez+eabc/XVV3P22Wfnh76lfy2HusuWLWPo0A86BNt+++1Ztmz9DsG+/e1vc84557Dpppt22/aYdaWFCxdyxhlnMHfuXAYNGsRPf/pTAC677LKKdeSss87K5501axZjxoxh1113Zfr06TQ0NFRdd8otWbKEJ554gvvuu4+pU6eyZs0aFixYUDGGsWPHsnr1apYtW8b223/wVPXW1rVs2TLuueeeDR7xV6uw/hQk9QOuAr5E1ln5bEkzIqL8OOuOiDizqDj6gqFDh/K5z30OgBNOOIErr7yy1b2J9pgzZw7PP/88V1xxBUuWLOn08sxqYYcddmDs2LEA7LHHHvl3+dxzz+Xcc8/d4Lzjx49n7ty5zJ8/n8mTJ3PwwQd3OI6jjz6ajTbaiB133JERI0bw3HPPMXbsWObMmdPhZbb42te+xqWXXspGG3V+P7/ITnbGAYsi4gUASbcDhwHFNb71US2HtaWvzz77bB555JH1yk6aNInzzz+fIUOG8PLLL+fjly5dypAhQz5U9rHHHqO5uZnhw4ezdu1aXnnlFZqamtypifUoAwYMyIf79euXNx9ddtll/OQnP1mv/L777suVV175oXGjRo1i4MCBPPvss1XVnUoq1dMFCxZwzDHHVCz/6KOPMmTIEJYuXdrmupqbm5k0aRIAr776KjNnzqShoYHDDz+8zbjKFZkUhgAvl7xeCoyvUO5ISfsCvwfOjoiXywtIOhU4FWDYsGEFhNqzvfTSSzz22GPsvffe3Hrrreyzzz6cc845G5xnwoQJTJs2jUmTJjFr1iy22GILtttuuw+VOe200/LD0SVLlnDooYc6IViv0daRwuLFixk6dCgNDQ28+OKLPPfccwwfPpxBgwaxcOFCFi9ezJAhQ7j99tu59dZbAbjgggsYN24cRxxxxHrLu+uuu5g8eTKLFy/mhRdeYOedd2aTTTbZ4JHCoEGD2HzzzXn88ccZP348N910E1/96lcrxtpiypQpHHrooR1KCFD7E82/AIZHxG7AQ8CNlQpFxLUR0RgRjYMHD+7WAHuCnXfemauuuopRo0axatWqqtoVDznkEEaMGMHIkSM55ZRTuPrqq/NpLYfa1jWqOLe2r6SnJK2VNLFk/FhJj0maK+lpSZV3Ka0Qv/nNb9h9990ZO3YsRxxxBFdffTVbb701DQ0NTJs2jQMPPJBRo0Zx9NFHM2bMGACeeeYZtt1224rLGzZsGOPGjePggw9m+vTpbLLJJlXFcfXVV3PyySczcuRIPvWpT+VNWNOnT2f69Olds7GlIqKQP2Bv4MGS1xcAF2ygfD/g9baWu8cee0S5/fbbL/bbb7/1xvcFixcvjjFjxtQ6jJoDmqOg73Jn/tL3+nlgBNAf+B0wuqzMcGA34CZgYsn4nYAd0/AngOXAoLbWWamO9ETl9bqt1/XggAMOqDh+8uTJcdddd3VzNB9WbR0p8khhNrCjpB0k9QcmATNKC0gqba+YAMwvMB6zWsjPrUXEu0DLubVcRCyJiKeB98vG/z4iFqbhPwCvAD5UrmMPPvhgrUPotMLOKUTEWklnAg+S7S1dHxFzJV1MlrFmAGdJmgCsBV4DphQVT281fPhwnn322VqHYa2r9tzaBkkaR3ak8Xwr033erY7dcMMNtQ6hakWeaCYiZgIzy8ZdWDJ8AVmzkpm1Ih1R3wxMjoj3K5WJiGuBawEaGxtbfw6CWRtqfaLZrLdbBgwteb19GlcVSZsD9wHfjIjHuzg2s/U4KZgVq81za61J5e8BboqIjj/hzKwdnBTMChQRa4GWc2vzgTtbzq2l82lI2lPSUuAo4BpJc9PsRwP7AlNKHgXj64WtUIWeUzCzqs6tzSZrViqf7xbglsIDNCvhIwUzqwtNTU00NTXVOow+z0nBzMxyTgpmZpZzUjAzs5yTgpmZ5ZwUzMws56RgZmY5JwUz6xV8SWvXcFIwM7Ock4KZmeWcFMzMLOekYGZmOScFMzPLOSmYmVnOScHMzHJOCmZmlnNSMDOznJOCmZnlnBTMCibpIEkLJC2SdH6F6ftKekrSWkkTy6ZNlrQw/U3uvqitr3JSMCuQpH7AVcDBwGjgWEmjy4q9BEwBbi2bd0vgO8B4YBzwHUkfKzpm69sKTQpt7SGVlDtSUkhqLDIesxoYByyKiBci4l3gduCw0gIRsSQingbeL5v3QOChiHgtIlYBDwEHdUfQ1ncVlhSq3ENC0mbA3wGziorFrIaGAC+XvF6axnXpvJJOldQsqXnFihUdCtQMij1SaHMPKfkecCmwpsBYzHq1iLg2IhojonHw4MG1Dsd6sCKTQpt7OZI+AwyNiPs2tCDvBRn02OflLwOGlrzePo0rel6zDqnZiWZJGwGXA+e0VdZ7QdaDzQZ2lLSDpP7AJGBGlfM+CBwg6WPpBPMBaZxZYYpMCm3t5WwG7AI8KmkJsBcwwyebrTeJiLXAmWQ/5vOBOyNirqSLJU0AkLSnpKXAUcA1kuameV8ja16dnf4uTuPMCtNQ4LLzPSSyZDAJOK5lYkS8Dmzd8lrSo8A3IqK5wJjMul1EzARmlo27sGR4NtlOU6V5rweuLzRAsxKFHSlUs4dkZmb1pcgjhTb3kMrGNxUZi5mZtc13NJuZWc5JwczMck4KZmaWc1IwM7Ock4KZmeWcFMzMLOekYGZmOScFMzPLOSmYmVnOScHMzHJOCmZmlnNSMDOznJOCmZnlnBTMzCznpGBmZjknBbNuIOkgSQskLZJ0foXpAyTdkabPkjQ8jd9Y0o2SnpE0X9IF3R279S1OCmYFk9QPuAo4GBgNHCtpdFmxk4BVETESuAK4NI0/ChgQEbsCewBfaUkYPV1TUxNNTU21DsPKOCmYFW8csCgiXoiId4HbgcPKyhwG3JiG7wb2lyQggI9KagA+ArwLvNE9YVtf5KRgVrwhwMslr5emcRXLpP7NXwe2IksQfwGWAy8B/xwRrxUdsPVdTgpm9W0csA74BLADcI6kEeWFJJ0qqVlS84oVK7o7RutFNpgUJF2a/h/VPeGY1adO1oVlwNCS19uncRXLpKaiLYCVwHHAAxHxXkS8AvwWaCxfQURcGxGNEdE4ePDgDoRolmnrSOGQ1K7pKx6sr+tMXZgN7ChpB0n9gUnAjLIyM4DJaXgi8HBEBFmT0RcAJH0U2At4rgMx9Dk+kd0xDW1MfwBYBQyUVHpyS0BExOaFRWZWXzpcFyJiraQzgQeBfsD1ETFX0sVAc0TMAK4Dbpa0CHiNLHFAdtXSjyXNTev6cUQ83dUbZ9Zig0khIs4FzpV0b0SUXy1h1md0ti5ExExgZtm4C0uG15Bdflo+35uVxpsVpaoTzS2VQNLmkrZs+Wtrvipu2JmabsqZI+k3Fa7dNqsrHa0LZj1FW81HQHZlA3AxsIbsumnS//WugiiZp+WGnS+RXYI3W9KMiJhXUuzWiJieyk8ALgcOau9GmHWXjtQFs56kqqQAnAvsEhGvtmPZ+Q07AJJabtjJk0JElLbNfpQPKplZvepIXTDrMapNCs8Db7Vz2ZVu2BlfXkjSGcDXgf6kqywqlDkVOBVg2LBh7QzDrEt1pC6Y9RjVJoULgP+WNAt4p2VkRJzV2QAi4irgKknHAd/ig8vySstcC1wL0NjY6KMJq6XC6oJZPag2KVwDPAw8A7xf5TzV3LBT6nbg/1W5bLNa6UhdMOsxqk0KG0fE19u57PyGHbJkMIns7sycpB0jYmF6+dfAQszqW0fqglmPUW1SuD+16/+CDx8yt/pgripv2DlT0heB98huDFqv6ciszrS7Lpj1JNUmhWPJrgwqv9dgg5fhVXHDzt9VuX6zetGhumDWU1SbFEYDpwP7kFWI/wKmFxWUWR1zXbBerdqkcCNZxx5XptfHpXFHFxGUWR1zXbBerdqksEtElD6C4hFJ81otbdZ7uS5Yr1ZtJztPSdqr5YWk8UBzMSGZ1TXXBevVqj1S2IPshp2X0uthwAJJz5A9Nni3QqIzqz+uC9arVZsU/JA6s4zrgvVqVSWFiHix6EDMegLXBevtqj2nYGZmfYCTgpmZ5ZwUzApWRQ+EAyTdkabPkjS8ZNpukh6TNDf1UrhJd8ZufY+TglmBSnogPJjsbuhjK3Q7exKwKiJGAlcAl6Z5G4BbgKkRMQZoIntOmFlhnBTMipX3QBgR75I9Iv6wsjKHkd0VDXA3sL8kAQcAT0fE7wAiYmVErOumuK2PclIwK1alHgiHtFYmItYCrwNbATsBIelBSU9J+vvWViLpVEnNkppXrFjRpRtgfYuTgln9aiB78N7x6f8RkvavVDAiro2IxohoHDx4cHfGaL2Mk4JZsarpgTAvk84jbAGsJDuq+HVEvBoRb5E9hv4zhUdsfZqTglmx8h4IJfUn64FwRlmZGXzQwdRE4OGICLIOqnaVtGlKFvsBfvieFarax1yYWQdU2QPhdcDNkhYBr5ElDiJilaTLyRJLADMj4r6abIj1GU4KZgWrogfCNcBRrcx7C9llqWbdws1HZmaWc1IwM7Ock4KZmeWcFMzMLOekYGZmOScFMzPLFZoUqnhk8NclzZP0tKT/kPTJIuMxM7MNKywpVPnI4P8BGlNn53cDPygqHjOrraamJpqammodhrWhyCOFNh8ZHBGPpGe6ADxO9lwYMzOrkSKTQjWPDC51EnB/pQl+LLCZdZaPVKpTFyeaJZ0ANAKXVZruxwKbmXWPIp99VM0jg5H0ReCbwH4R8U6B8ZiZWRuKPFJo85HBkv4KuAaYEBGvFBiLmZlVobCkkLoVbHlk8HzgzpZHBkuakIpdBgwE7pI0R1L5c+bNzKwbFfro7CoeGfzFItdvZmbtUxcnms3MrD44KZiZWc5JwczMck4KZt2giueADZB0R5o+S9LwsunDJL0p6RvdFbP1TU4KZgWr8jlgJwGrImIkcAVwadn0y2nljn+zruSkYFa8Np8Dll7fmIbvBvaXJABJhwOLgbndFK/1YU4KZsWr5jlgeZl0j8/rwFaSBgLnAd/d0Ar8fDDrKk4KZvXtIuCKiHhzQ4X8fDDrKoXevGZmQHXPAWsps1RSA7AFsBIYD0yU9ANgEPC+pDURMa34sK0vclIwK17+HDCyH/9JwHFlZWYAk4HHgInAwxERwOdbCki6CHjTCcGK5KRgVrCIWCup5Tlg/YDrW54DBjRHxAzgOuBmSYuA18gSh1m3c1Iw6wZVPAdsDXBUG8u4qJDgzEr4RLOZmeWcFMzMLOekYGZmOScFMzPLOSmYmVnOScHMzHJOCmZmlnNSMDOznJOCmZnlnBTMzCznpGBmZjknBTMzyxWaFKrorHxfSU9JWitpYpGxmJmVampqoqmpqdZh1J3CkkKVnZW/BEwBbi0qDjPrfv7B7bmKfHR23lk5gKSWzsrntRSIiCVp2vsFxmFmZlUqsvmoms7Kq+JOyc3MukePONHsTsnNzLpHkUmhms7KzcysjhSZFPLOyiX1J+tzdkaB6zOrS1VchTdA0h1p+ixJw9P4L0l6UtIz6f8Xujt263sKSwoRsRZo6ax8PnBnS2flkiYASNpT0lKyvmmvkTS3qHjMaqHKq/BOAlZFxEjgCuDSNP5V4MsRsSswGbi5e6K2vqzIq4+q6ax8Nlmzkllv1eZVeOn1RWn4bmCaJEXE/5SUmQt8RNKAiHin+LCtr+oRJ5rNerBqrsLLy6Qj7NeBrcrKHAk81VpC8BV61lWcFMzqnKQxZE1KX2mtjK/Qs67ipGBWrGquwsvLSGoAtgBWptfbA/cAJ0bE84VHa32ek4JZsaq5Cm8G2YlkgInAwxERkgYB9wHnR8Rvuy1i69OcFMwKVM1VeMB1wFaSFgFfB1ouWz0TGAlcKGlO+tummzfB+phCrz4ys6quwltDdll2+Xz/APxD4QGalfCRgpmZ5ZwUzMws56RgZmY5JwUzM8s5KZiZ4d7iWjgpmFmn+Qe193BSMDOznJOCmZnlnBTMzCznpGBmZjknBTMzyzkpmJlZzknBzMxyTgpmZpZzUjAzq6Cv3pDnpGBm7dZXfzD7AicFMzPLOSmYmVmu0KQg6SBJCyQtknR+hekDJN2Rps+SNLzIeMxqpTN1QdIFafwCSQd2Z9z2gb7SZFZYUpDUD7gKOBgYDRwraXRZsZOAVRExErgCuLSoeMxqpTN1IZWbBIwBDgKuTsvrVn3lB9GKPVIYByyKiBci4l3gduCwsjKHATem4buB/SWpwJjMaqEzdeEw4PaIeCciFgOL0vKsxnprolREFLNgaSJwUEScnF7/DTA+Is4sKfNsKrM0vX4+lXm1bFmnAqcCDBs2bI8XX3yxkJit55L0ZEQ01jqOSjpTF4CLgMcj4pY0/jrg/oi4u2wdriO2QdXWkR5xojlC0vnvAAAEXklEQVQiro2IxohoHDx4cK3DMas7riPWVYpMCsuAoSWvt0/jKpaR1ABsAawsMCazWuhMXahmXrMuU2RSmA3sKGkHSf3JTpbNKCszA5ichicCD0dR7VlmtdOZujADmJSuTtoB2BF4opvitj6ooagFR8RaSWcCDwL9gOsjYq6ki4HmiJgBXAfcLGkR8BpZZTHrVTpTF1K5O4F5wFrgjIhYV5MNsT6hsBPNRWlsbIzm5uZah2F1pp5PNHc31xGrpFedaDYzs+7hpGBmZjknBTMzyzkpmJlZrsedaJa0Aqh0u+bWwKsVxvdm3uYPfDIifNcWPb6O1HuM9R4fdLKO9Lik0BpJzX3t6hNvs7VHT3jv6j3Geo8POh+jm4/MzCznpGBmZrnelBSurXUANeBttvboCe9dvcdY7/FBJ2PsNecUzMys83rTkYKZmXWSk4KZmeV6XFLoTAfoPVUV2zxF0gpJc9LfybWIs6tIul7SK6k3skrTJenK9H48Lekz3R1jT9PWd6gG8az3GUvaUtJDkham/x+rcYxDJT0iaZ6kuZL+rp7ilLSJpCck/S7F9900fof027co/Rb2b9eCI6LH/JE9dvh5YATQH/gdMLqszOnA9DQ8Cbij1nF3wzZPAabVOtYu3OZ9gc8Az7Yy/RDgfkDAXsCsWsdcz3/VfIfq4TMGfgCcn4bPBy6tcYzbAZ9Jw5sBvwdG10uc6fs/MA1vDMxK9eFOYFIaPx04rT3L7WlHCp3pAL2nqmabe5WI+DVZnwKtOQy4KTKPA4Mkbdc90fVIdfcdauUzLq27NwKHd2tQZSJieUQ8lYb/DMwHhlAncabv/5vp5cbpL4AvkP32QQfi62lJYQjwcsnrpWlcxTIRsRZ4HdiqW6IrRjXbDHBkakq5W9LQCtN7k2rfE8v0lPfr4xGxPA3/Efh4LYMplZqh/4psb7xu4pTUT9Ic4BXgIbIjwtXptw868Fn3tKRglf0CGB4Ru5F9MW5so7xZXYus7aMurpeXNBD4KfC1iHijdFqt44yIdRExlqzv7nHApzu7zJ6WFDrTAXpP1eY2R8TKiHgnvfw3YI9uiq1W3Jl9+/SU9+tPLc2A6f8rNY4HSRuTJYSfRMTP0ui6izMiVgOPAHuTNae2dLXc7s+6pyWFznSA3lO1uc1l7ekTyNo+e7MZwInpKqS9gNdLDudtfdXUm3pQWncnA/fWMBbSucjrgPkRcXnJpLqIU9JgSYPS8EeAL5HV/UfIfvs6Fl8tz+538Iz7IWRXATwPfDONuxiYkIY3Ae4CFgFPACNqHXM3bPM/AXPJrip5BPh0rWPu5PbeBiwH3iNrEz0JmApMTdMFXJXej2eAxlrHXO9/lb5DdfgZbwX8B7AQ+HdgyxrHuA9Z09DTwJz0d0i9xAnsBvxPiu9Z4MI0fkT67VuUfgsHtGe5fsyFmZnlelrzkZmZFchJwczMck4KZmaWc1IwM7Ock4KZmeWcFMzMLOekYGZmuf8P/CFi5Il0XDwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# Density functions\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.subplots_adjust(wspace=0.5,hspace=0.5)\n",
    "\n",
    "# Bernoulli\n",
    "rv = stats.bernoulli(0.4)                   # So called \"Frozen\" rv: just means parameters have been fixed.\n",
    "\n",
    "a1 = fig.add_subplot(1,2,1)\n",
    "x=[0,1]\n",
    "_ = a1.vlines(x,0,rv.pmf(x))\n",
    "_ = a1.set(ylabel=\"pmf\",title=\"Bernoulli Distribution\")\n",
    "_ = a1.text(0.5,0.8,\"p=0.4\", transform=a1.transAxes)\n",
    "\n",
    "# Binomial\n",
    "n,p = 30, 0.4 \n",
    "rv = stats.binom(n,p)\n",
    "\n",
    "a2 = fig.add_subplot(1,2,2)\n",
    "x =  np.arange(n+1)\n",
    "_ = a2.vlines(x,0,rv.pmf(x))\n",
    "\n",
    "_ = a2.set(ylabel=\"pmf\",title=\"Binomial Distribution\")\n",
    "_ = a2.text(0.5,0.8,\"n=30, p=0.4\", transform=a2.transAxes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multinomial distribution\n",
    "\n",
    "Generalization of binomial. $k$ outcomes with given individual probabilities $p_i, i=1 \\dots k$, and $n$ trials. Simplest example is throwing dice $n$ times. Random variables in the set $\\{X_i\\}$ are the number of times each outcome is obtained\n",
    "\n",
    "$$\n",
    "    f_{\\{X_i\\}} (x_1 \\dots x_k;n,p_1 \\dots p_k) = P [X_1 = x_1 \\dots X_k=x_k]=\\frac{n!}{x_1!  \\dots x_k!} p_1^{x_1} \\dots p_k^{x_k} \\hspace{2cm} \\sum_{i=1}^k x_i = n\n",
    "$$\n",
    "\n",
    "Note that this is a **multivariate and joint** probability mass function for $k$ random variables.\n",
    "\n",
    "**Example**. For six-sided dice, with equal probabilities for each face to appear, the probability of obtaining 5 sixes and 2 each of the other faces is\n",
    "$$\n",
    "f_{\\{X_i\\}} (x_1= x_2 = \\dots x_5=2,x_6=5;n=15,p_1=p_2=\\dots=p_6=1/6) = \\frac{15!}{(2!)^5 5!} (\\frac{1}{6})^{15} = 0.0072\n",
    "$$\n",
    "\n",
    "The choice of $\\{x_i\\}$ in this example is very specific. Generally, we are likely to be more interested in finding the probabilities of events such as the one in the next example.\n",
    "\n",
    "**Example**. What is the probability of getting at least one six in four throws of a fair die?\n",
    "\n",
    "We could list all possibilities (1 six, 2 sixes, etc.) and the corresponding cases for the remaining throws, but a better way to solve the problem is to consider the probability of no sixes in 4 throws:\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "P[\\textrm{at least 1 six in 4 throws}] & = 1 - P[\\textrm{no sixes in 4 throws}]\\\\\n",
    "    & = 1 - (\\frac{5}{6})^4 \\\\\n",
    "    & = 0.52\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "So, you need only 4 throws for a better than even chance of getting a six.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0007242685563578615"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "0.0007242685563578548"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Example 1\n",
    "\n",
    "x = [2,2,2,2,2,5]\n",
    "p = [1/6 for p_i in range(6)]\n",
    "\n",
    "stats.multinomial.pmf(x,n=15,p=[1/6 for p_i in range(6)])\n",
    "(math.factorial(15)/(math.factorial(5)*(2**5)))*(1/6)**15\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometric distribution\n",
    "Binary outcome, with probability of success $p$. Random variable $X$ is number of trials ***before*** first success. The *pmf* is\n",
    "\\begin{equation}\n",
    "f_X(k;p) = (1-p)^{(k-1)} p^k \\hspace{2cm} k \\in \\{1,2,3 \\dots\\}\n",
    "\\end{equation}\n",
    "\n",
    "Expectation and variance of $X$:\n",
    "\\begin{equation}\n",
    "\\mu = \\frac{1}{p} \\hspace{2cm} \\sigma^2 = \\frac{(1-p)}{p^2}\n",
    "\\end{equation}\n",
    "### Negative Binomial distribution\n",
    "Binary outcome, trials continue until $r$ failures are seen. Random variable $X$ is the number of successes before stopping.\n",
    "The *pmf* is\n",
    "\\begin{equation}\n",
    "f_X(k;p,r) = \\binom{k+r-1}{k} \\cdot p^k (1-p)^r \\hspace{2cm} k \\in {0,1,2, \\dots}\n",
    "\\end{equation}\n",
    "\n",
    "The superscript on the binomial coefficient in the *pmf* is interesting: since the $r$-th failure is always the last one, the coefficient has to count the ways of choosing $k$ successes out of $(k+r-1)$ trials.\n",
    "\n",
    "Expectation and variance of $X$:\n",
    "\\begin{equation}\n",
    "\\mu = \\frac{pr}{1-p} \\hspace{2cm} \\sigma^2 = \\frac{pr}{(1-p)^2}\n",
    "\\end{equation}\n",
    "\n",
    "The distribution is called negative binomial because\n",
    "$$ \\binom{k+r-1}{k} = (-1)^k \\binom{-r}{k} $$\n",
    "Proof:\n",
    "\n",
    "The binomial theorem can be written as:\n",
    "$$ (1+x)^a = 1 + a x + \\frac{a(a-1)x^2}{2!} + \\cdots = \\sum_{k=0}^\\infty \\binom{a}{k} x^k$$\n",
    "\n",
    "When $a$ is negative, $$ \\binom{a}{k} = \\frac{a(a-1)\\dots(a-k+1)}{k!} = \\frac{(-1)^k (-a)(-a+1)\\dots(-a+k-1)}{k!}$$\n",
    "\n",
    "Setting $r=-a$, we have $$ \\binom{-r}{k} = (-1)^k \\frac{r(r+1)\\dots(r+k-1)}{k!} = (-1)^k \\binom{k+r-1}{k}$$\n",
    "\n",
    "\n",
    "### Hypergeometric distribution\n",
    "\n",
    "$n$ draws from a population with $N$ objects, of which $K$ are of one kind. The probability that $k$ of the $n$ draws are of that kind is given by the *pmf*\n",
    "\\begin{equation}\n",
    "   f_X (k; N,K,n) = \\frac{\\binom{K}{k} \\cdot \\binom{N-K}{n-k}}{\\binom{N}{n}}\n",
    "\\end{equation}\n",
    "\n",
    "Expectation and variance:\n",
    "\\begin{equation}\n",
    "    \\mu = n\\frac{K}{N} \\hspace{2cm} \\sigma^2 = n \\frac{K}{N} \\frac{(N-K)}{N} \\frac{(N-n)}{(N-1)}\n",
    "\\end{equation}\n",
    "\n",
    "[Nice proof on StackExchange](https://math.stackexchange.com/questions/1380460/derivation-of-mean-and-variance-of-hypergeometric-distribution)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEICAYAAAD2u0vkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xm8XVV99/HP1yBQCEOQiBAyIWEIQ0O5DC0K97EIQWugLUIYJLbUPFCorZEiVgXFGS08UsMDsaTyiMwqXAVEKsRqGcwFIhAwJoRAEiIEEqQIEkJ+zx9rXdg53OHcYd99z7nf9+t1X/fsYe39O+fsdX57rz0sRQRmZmZVeEvVAZiZ2fDlJGRmZpVxEjIzs8o4CZmZWWWchMzMrDJOQmZmVhknoSFK0jhJL0oaMdSXLelSSZ/Jr1slrRiI5eblvVvSooFantVP0q2SZni9nc67TNLh/VjXSZJ+0tfynSxvoaTW/Pqzkq4cwGX/i6R/H6jlvUlENM0fMB24F/g98Ex+/feAqo6tJs5lwOElLfvDwGvAi/nvceA/gN36uKxf9LJMK7CiH/EHsGvV31GF28UzwJaFcX8HzBuEdX8WuHKQ3mfkOvoi8CxwNbBt1Z9/H76rTusw8G1gHfA/+e9h4MvANn1Yz7eBLwzWd9nf+tuXv6Y5EpL0ceAbwNeAdwA7AKcBhwCbVhhar0napJ+LuDsiRgLbAIcDLwP3Sdq738HVKONIbZgbAfxj1UEMgj/O2+guwCjSD2czuSAitgJGA38DHAz8t6QtB3IlA/BbUb2q9ygGaK9kG9Ke1V/3MN9mwNeBJ4GngUuBPypM/wiwBFgDtAE7FaYF6ahqMWnv5vPAO4G7gBeA64BNC/P/BbAAeD7Ps28e/x1gAykxvAicDUzIyz81x/ZfhXGb5HLbkY5ongLWAjd28R4/TCdHL8CPgBvy69plfxhYmt/X48BJwJ7AH3jjqOr5eGPP7P8Ct+TP/HAKe2vkPSngX0h7ucuAkwpxzAP+rrN48/su7iUfT82eWY5rXv5cFwLTCtO+DcwGbs7v5V7gnVVvn73YjpcB5+Ttb9s8bqMjIWAP4PY8zyLguMK0twE/zNvjfOALxW2BtJO2PE+/D3h3Hj+VtOf+av7cf1X8rkj15nlg78KyRpO24bd3t7138T43Otol1aufdLaNdGwfpHq7Nm+fRxXm3YlUV9eQ6u5HCtM+C1wPXJm3h4eA3YBPko44lwNHdLHedwJ3AM+RtuPvUjhao+cjoS/UjNsKWAWc2cl2L+CiHNMLOc69gZn5O1mXv5cfFtb9CeBB4BVgk2I8+X3fAFyb3/f9pKTf1ef/bdK2smX+TjfwRkvKTtQcWQHTSHXv+fyZ7VnzuZyVY/tdjmHz7rb7ZjkS+lNSRbmph/m+QtoIpwC7AmOAcwEkvYd0yHwcsCPwBHBNTfkjgf1JezVnA3OAk4GxpI3mhLys/YC5wP8m/TBcBrRJ2iwiPkRKNB+IiJERcUFh+YeRfmSP7CT27wBbAHsBbydttL3xfeDdtSPzntnFpIq9FfBnwIKIeJR0JHl3jnPbQrETgS+SKtYvOlnXO4DtSZ/vDGCOpN17CjAiDs0v/ziv89qaWN9K+pH9Cekz+AfguzXLng58jrR3vSTH2UjaSRX7rNoJ+bu6HbiK9P6nA5dImpxnmU1K4O8gfe615zfmk7b97fIyrpe0eUT8GPgScG3+3P+4WCgiXiFtPycURh8H/Cwinulue+/pzUoaBRwD3NPNbAeREu72wAXA5ZKUp11D2unZCTgW+FKuyx0+QKo7o4AHgNtI58LHAOfnWDsNjfR7sBOpTo6lH0drEfE/pO/uTXUQOAI4lPTbtA3ps30uIuaQkt8F+Xv5QKHMCcD7SYlxfSfLPJqUgDu+6xtz/ekuxt8DRwFP5fWNjIinivNI2o3UfPpPpB2RW4AfSiq2Nh1H2rGZCOxLSrhdapYktD3wbPHLkHSXpOclvSzp0LzRzgQ+FhFr8kbxJVJFhrT3Pzci7s+V7pPAn0qaUFjPBRHxQkQsJLXz/iQilkbE74Bbgf3yfDOByyLi3oh4LSKuIO2xHNzD+/hsRPw+Il4ujpS0I2njOC0i1kbEqxHxs15+Rk+RNsjObAD2lvRHEbEqv7/u3BQR/x0RGyLiD13M85mIeCXHeTNpw+yvg4GRwFciYl1E3EE6wiv+OP4gIn6Zt4Xvkn50G825wD9IGl0z/i+AZRHxHxGxPiIeAL4HfDA3i/41cF5EvBQRjwBXFAtHxJUR8Vwu+6+kHbcedw6yq3ijrkDaEbkqv+7L9n6/pOdJRxnj6DoZADwREd+KiNfye9oR2EHSWFJz+yci4g8RsQD4d+CUQtmfR8RteXu4nvTD+ZWIeJWUwCZI2pYaEbEkIm7P2/Bq4ELSTmJ/dFUHXyXt0O1BOn/9aESs6mFZF0fE8trfioL7IuKG/D4vBDan59+fehwP3Jw/m1dJR6h/RNp5Lcb2VESsIe00dlsHmyUJPQdsX2wfjYg/y3vvz5He52jSkcR9OTk9D/w4j4e0x/NEofyLueyYwnqeLrx+uZPhkfn1eODjHevJ6xqb19Gd5V2MHwusiYi1PZTvzhhSk8VG8t7P8aSjnlWSbpa0Rx/j7LA2L7fDE/T83uuxE7A8IjbULLv4Hf228Pol3vhOGkZEPExKrufUTBoPHFSzXZ1EOvIZTWqWKX43G31Pks6S9Kik3+Wy25B24OpxJ7CFpIPyjtkU4AeFuHq7vf9Jrp+bk5p3fy5p8y7mff07jYiX8suRefkdO5QdareH2jr6bE5mHcMdy9qIpB0kXSNppaQXSE169X5WXemqDt4BfJN0JPuMpDmStu5hWT3Vwden5/rScbTYX7W/kxvyuvpcB5slCd1N2vM6upt5niVtdHtFxLb5b5tIJ0ch7aWM75g5N328DVjZh3iWA18srGfbiNgiIq7O07t6dHlX45cD23W2x9YLfwn8vNOVpj3F95L2MH8NfKuPcXYYVXMCdhzp84XUXLRFYdo7elhW0VPAWEnF7XYcffuOhrrzSOcoi5V7OakJrLhdjYyI04HVwHpg58L8YzteSHo3qQn5OGBUTgC/IzU7QQ/faf7hvo501HkC8KPCj39P23t3y32VdPQykdSk3RtPkerFVoVxA7U9fIn0mewTEVuTmt3VfZGuSRpJOn/aVR28OCL2ByaTmuX+uWNSF4vsqQ4Wv/u3kLaLjjr4El3XwZ6WW/s7qbyuPn/mTZGEIuJ50nmASyQdK2krSW+RNIV0sq0jY38LuEjS2wEkjZHUcf7lauBvJE3JbdlfAu6NiGV9COlbwGl5r1GStpT0/kJleZp0VVC9728VqbnvEkmjJL1V0qE9lZM0QtJESf9GOsH/uU7m2UHS0TlpvEI6GdlxpPE0sHNNe2+9Pidp0/zj9xekphBIJ6//StIWknYlXYxR1N1ncy+pAp2dP4NWUpt/7bm7hhcRS0gndT9aGP0jYDdJH8rv/62SDpC0Z04S3wc+mz/bPdi4WWorUpJaDWwi6VyguLf9NKlpqrvfhKtIR80n8UZTHPS8vXcpNyP+DWkHcWlP8xdFxHLSRRBflrS5pH1J29NA3COzFaku/E7SGN5ICr0iaTNJ+wM3ki6s+I9O5jkgf3ZvJe2k/YGN62DdvxUF+0v6q9w69E+kut1x3m0BcGL+fZjKxs2MTwNvk7RNF8u9Dni/pD/P8X48L/uuPsQINEkSAoh0gn8WaW/v6fx3Gekqko4P6BOkk9X35EPs/yS3iUfEfwKfIbWxryJdHVNsA+9NLO2kvdhvkja8JWx8cu7LwKdz08WbTkB34UOktuNfk66i+adu5v1TSS+SrrSZR/qxOSAiHupk3reQPrenSE0FhwGn52l3kK6C+a2kZ+uME9Lh+Nq8zO+SzmX9Ok+7iHS1z9Ok9v3v1pT9LHBF/mw2Oo8UEetISeco0pHtJcAphWU3m/PJO1Hw+sntI0jb5VOkz/mrpHM7AGeSmth+SzoZfzXpBwLSCfkfA78hNaf8gY2bdDp2Ep6TdH9nwURExz14O5F2ijrG97S9d+ZXeRtdS7qA4i/zOYTeOoF0tedTpObB83Jd7q/PAX9COlq8mZTge+NsSf9DatL/f6SrEf+sppm6w9akRL6W9N08R7rVBOByYHKuDzf2Yv03kXYY1pJ+O/4qH3VCugXgA6Sr204iJUgAcl26Glia17lRE15ELCIdFf4bqQ5+gHSR1bpexLYRRbhTO7NmJOmrwDsiYtCfPmBWr6Y5EjIb7iTtIWnf3CR2IKlp6gc9lTOrUuPfbWtmHbYiNaXsRGru/Fd6vnfOrFJujjMzs8q4Oc7MzCrTNM1x22+/fUyYMKHqMGyIue+++56NiNonDwxLriPWmarrSNMkoQkTJtDe3l51GDbESHqi57mGB9cR60zVdcTNcWZmVhknITMzq4yTkJmZVcZJyMzMKuMkZGZmlSk1CUmaKmmRpCWSavtGKc7315JCUkth3CdzuUWFJ12bmVkTKe0S7fyI9tnAe0kdKs2X1JZ7fCzOtxXpqa73FsZNJj0peC/SI0j+U9Juhc6ozMysCZR5JHQgsCRS99frSH2+dNbp3OdJj6MvdhN9NHBN7lr3cdKj4Q8sMVYzM6tAmUloDBv3V7KCjXuJRNKfAGMj4ubels3lZ0pql9S+evXqfgfc2tpKa2trv5dj1oxcP6wMlV2YkHtwvJDUM1+fRMSciGiJiJbRo/1kFjOzRlPmY3tWUujnnNTHebEf8q1IfcrPS92U8w6gTdK0OsqamVkTKPNIaD4wSdJESZuSLjRo65gYEb+LiO0jYkJETCD1fz4tdxXcBkzP/bNPBCYBvywxVjMzq0BpSSgi1pP6vL8NeBS4LiIWSjo/H+10V3YhcB3wCPBj4AxfGWfNqKfbGCTNkvSIpAcl/VTS+MK01yQtyH9ttWXNGkGpT9GOiFuAW2rGndvFvK01w18EvlhacGYVq/M2hgeAloh4SdLpwAXA8XnayxExZVCDNhtgfmKCWXV6vI0hIu6MiJfy4D2k86NmTcNJyKw6dd2KUHAqcGthePN8i8I9ko7prMBA38ZgNtCaplM7s2Ym6WSgBTisMHp8RKyUtAtwh6SHIuKxYrmImAPMAWhpaYlBC9isTj4SMqtOXbciSDoc+BTp6tFXOsZHxMr8fykwD9ivzGDNyuAkZFadbm9jAJC0H3AZKQE9Uxg/StJm+fX2wCGkq0nNGoqb48wqEhHrJXXcxjACmNtxGwPQHhFtwNeAkcD1+abuJyNiGrAncJmkDaSdya/UPhzYrBE4CZlVqKfbGCLi8C7K3QXsU250ZuVzc5yZmVXGScjMzCrjJGRmZpVxEjIzs8o4CZmZWWWchMzMrDJOQmZmVhknITMzq4yTkJmZVabUJFRHr5GnSXoo9wz5C0mT8/gJkl4u9Bp5aZlxmplZNUp7bE+dvUZeFRGX5vmnARcCU/O0x9xrpJlZcyvzSKieXiNfKAxuCbi/EzOzYaTMJFRXr5GSzpD0GHAB8NHCpImSHpD0M0nv7mwF7jXSzKyxVX5hQkTMjoh3Ap8APp1HrwLGRcR+wCzgKklbd1J2TkS0RETL6NGjBy9oMzMbEGUmobp6jSy4BjgGICJeiYjn8uv7gMeA3UqK0woigo9+9KPsuuuu7Lvvvtx///3dzj9t2jT23nvvQYrOrDqXXnop++yzD1OmTOFd73oXjzxSXvdNw6kelpmE6uk1clJh8P3A4jx+dL6wAUm7AJOApSXGatmtt97K4sWLWbx4MXPmzOH000/vct7vf//7jBw5chCjMxt4EcGGDRt6nO/EE0/koYceYsGCBZx99tnMmjWrtJiGUz0sLQlFxHqgo9fIR4HrOnqNzFfCAZwpaaGkBaRmtxl5/KHAg3n8DcBpEbGmrFib0bJly9hjjz046aST2HPPPTn22GN56aWXeix30003ccoppyCJgw8+mOeff55Vq1a9ab4XX3yRCy+8kE9/+tOdLMVsaFu2bBm77747p5xyCnvvvTfLly/vsczWW79xRuD3v/89uafbbtfRlzoIw6seltqzah29Rv5jF+W+B3yvzNiGg0WLFnH55ZdzyCGH8Ld/+7dccsklrFy5kjvvvPNN806fPp1zzjmHlStXMnbsG62oO++8MytXrmTHHXfcaP7PfOYzfPzjH2eLLbYo/X2YlWHx4sVcccUVHHzwwQAcf/zxLFq06E3zzZo1i1NOOQWA2bNnc+GFF7Ju3TruuOOOHtfRWR0866yz+NjHPuZ6mLl77yY2duxYDjnkEABOPvlkLr74Ym688cZ+L3fBggU89thjXHTRRSxbtqzfyzOrwvjx419PQADXXnttj2XOOOMMzjjjDK666iq+8IUvcMUVV3Q7f2d18KyzzuKiiy7qX/A0Tz10Empitc0FknrcAxszZsxGTRMrVqxgzJiNr6y/++67aW9vZ8KECaxfv55nnnmG1tZW5s2bV8r7MCvDlltuudFwPUdCHaZPn97teZoOndVBwPWwKCKa4m///feP/jrssMPisMMO6/dyhoLHH388gLjrrrsiIuLUU0+Nr3/96z2W+9GPfhRTp06NDRs2xN133x0HHHBAj+vZa6+9BiTmMgDtMQS2z6Hw11kd6c0232z1o7fb7W9+85vXX7e1tUXH57lixYp4z3ve0+k6+lIHIwa3HlZdRyq/T8jKs/vuuzN79mz23HNP1q5dW9ee2/ve9z522WUXdt11Vz7ykY9wySWXvD5tyhQ/RcmGr29+85vstddeTJkyhQsvvPD1prhVq1axySadNyr1pQ7C8KqHbo5rYptssglXXnllr8pIYvbs2Z1OW7BgwZvGTZgwgYcffrhP8ZlVpS/b7Te+8Y1Ox99zzz2cccYZnU7rSx2E4VUPnYTMzPrhzDPPrDqEhubmuCbVyHtGZs3AdbA+TkJmZlYZJyEzM6uMk5CZlaK1tZXW1taqw7AhzknIzMwq4yRkViFJUyUtkrRE0jmdTJ8l6RFJD0r6qaTxhWkzJC3OfzNqy5o1Aichs4rk7kpmA0cBk4ETJE2ume0BoCUi9iU9Uf6CXHY74DzgIOBA4DxJowYrdrOB4iRkVp0DgSURsTQi1pE6djy6OENE3BkRHc//v4fUOSTAkcDtEbEmItYCtwNTByluswHjJGRWnTFAsSObFXlcV04Fbu1NWUkzJbVLal+9enU/wzUbeE5CZg1A0slAC/C13pSLiDkR0RIRLaNHjy4nOLN+KDUJ1XHS9TRJD0laIOkXxfZwSZ/M5RZJOrLMOM0qshIYWxjeOY/biKTDgU8B0yLild6UNRvqSktCdZ50vSoi9omIKaQTrhfmspOB6cBepHbuS/LyzJrJfGCSpImSNiVt823FGSTtB1xGSkDPFCbdBhwhaVS+IOGIPM6soZR5JFTPSdcXCoNbApFfHw1cExGvRMTjwJK8PLOmERHrgTNJyeNR4LqIWCjpfEnT8mxfA0YC1+cWg7Zcdg3weVIimw+cn8eZNZQyn6Ld2YnTg2pnknQGMAvYFHhPoew9NWU7PekKzAQYN27cgARtNpgi4hbglppx5xZeH95N2bnA3PKiMytf5RcmRMTsiHgn8Ang070s65OuZmYNrMwk1NsTp9cAx/SxrJmZNaAyk1A9J10nFQbfDyzOr9uA6ZI2kzQRmAT8ssRYzcysAqWdE4qI9ZI6TrqOAOZ2nHQF2iOiDTgzX376KrAWmJHLLpR0HfAIsB44IyJeKytWMzOrRqnde9dx0vUfuyn7ReCL5UVnZmZVq/zCBDMzG76chMzMrDJOQmZmVhknITMzq4yTkJmZVcZJyMzMKuMkZGZmlXESMjOzyjgJmZlZZZyEzMysMk5CZmZWGSchMzOrjJOQmZlVxknIzMwq4yRkZmaVKTUJSZoqaZGkJZLO6WT6LEmPSHpQ0k8ljS9Me03SgvzXVlvWzMwaX2md2kkaAcwG3gusAOZLaouIRwqzPQC0RMRLkk4HLgCOz9NejogpZcVnZmbVK/NI6EBgSUQsjYh1wDXA0cUZIuLOiHgpD94D7FxiPGZmNsSUmYTGAMsLwyvyuK6cCtxaGN5cUrukeyQd01kBSTPzPO2rV6/uf8RmZjaoSmuO6w1JJwMtwGGF0eMjYqWkXYA7JD0UEY8Vy0XEHGAOQEtLSwxawGZmNiDKPBJaCYwtDO+cx21E0uHAp4BpEfFKx/iIWJn/LwXmAfuVGKuZmVWgzCQ0H5gkaaKkTYHpwEZXuUnaD7iMlICeKYwfJWmz/Hp74BCgeEGDmZk1gdKa4yJivaQzgduAEcDciFgo6XygPSLagK8BI4HrJQE8GRHTgD2ByyRtICXKr9RcVWdmZk2g1HNCEXELcEvNuHMLrw/votxdwD5lxtZfra2tAMybN6/SOMzMGpmfmGBWoTpu6D5U0v2S1ks6tmaab+i2hjckro4zG47qvKH7SeDDwFmdLMI3dFvDcxIyq87rN3QDSOq4ofv1JBQRy/K0DVUEaFY2N8eZVae3N3TX6vGGbrOhzkdCZo2rxxu6Jc0EZgKMGzeuihjNuuUjIbPq1HVDd1fquaE7IuZEREtEtIwePbp/0ZqVwEnIrDo93tDdFd/Qbc3CScisIhGxHui4oftR4LqOG7olTQOQdICkFcAHSTdwL8zF9wTaJf0KuBPf0G0NyueEzCpUxw3d8+mki5NGuKHbrB4+EjIzs8o4CZmZWWWchMzMrDLdJiFJX83/Pzg44Zg1FteRgdHa2vr6Q4FteOnpSOh9Sn0sfHIwgjFrQK4jZv3Q09VxPwbWAiMlvVAYLyAiYuvSIjNrDK4jZv3Q7ZFQRPxzRGwL3BwRWxf+tnLlMnMdMeuvui5MiIijASRtLWm7jr+eytXRV8osSY9IelDSTyWNL0ybIWlx/pvRmzdlNtj6WkfMhru6blbND0E8H/gDEHl0ALt0U6aevlIeAFoi4iVJpwMXAMfnynse0JLXc18uu7ZX785skPSljphZ/U9M+Gdg74h4thfLrqevlDsL898DnJxfHwncHhFrctnbganA1b1Yv9lg6ksdMRv26r1P6DHgpV4uu7d9pZwK3NqbspJm5v5U2levXt3L8MwGVF/qiNmwV++R0CeBuyTdC7zSMTIiPjoQQUg6mdT0dlhvykXEHGAOQEtLS/Qwu1mZSq0jZs2q3iR0GXAH8BBQbzfDdfWVIulw4FPAYRHxSqFsa03ZeXWu16wKfakjZsNevUnorRExq5fLfr2vFFJSmQ6cWJxB0n6kyjs1Ip4pTLoN+JKkUXn4CHwzoA1tfakjZsNevUno1nz1zw/ZuKlhTVcFImK9pI6+UkYAczv6SgHaI6IN+BowErg+3XTOkxExLSLWSPo8KZEBnN/dusyGgF7XETOrPwmdQLrctPZen24vP62jr5TDuyk7F5hbZ3xmVetTHTEb7upNQpOBvwfeRapoPwcuLSsoswbkOmLWB/UmoSuAF4CL8/CJedxxZQRl1oBcR8z6oN4ktHdETC4M3ynJ/dmbvcF1xKwP6r1Z9X5JB3cMSDoIaC8nJLOG5Dpi1gf1HgntT7oR78k8PA5YJOkh0uPq9y0lOrPG4Tpi1gf1JqGppUZh1vhcR8z6oK4kFBFPlB2IWSNzHTHrm3rPCZmZmQ04JyEzM6uMk5CZmVXGScjMzCrjJGRmZpVxEjIzs8o4CZlVSNJUSYskLZFU+wRuJB0q6X5J6yUdWzNthqTF+W/G4EVtNnCchMwqImkEMBs4ivQU7hMkTa6Z7Ungw8BVNWW3A84DDgIOBM4rdAJp1jCchMyqcyCwJCKWRsQ64Brg6OIMEbEsIh7kzV2GHwncHhFrImItcDt+aoM1oFKTUD+bGl6TtCD/tZUZp1lFxgDLC8Mr8riyy5oNGfU+O67XCk0N7yVVkPmS2iKi+Hj7jqaGszpZxMsRMaWs+MyGg9zl+EyAcePGVRyN2ZuVeSTUn6YGs+FgJTC2MLxzHjdgZSNiTkS0RETL6NGj+xyoWVnKTEL9bS7YXFK7pHskHdPZDJJm5nnaV69e3Z9YzaowH5gkaaKkTYHpQL1Nz7cBR0galS9IOCKPM2soQ/nChPER0ULqJvn/SHpn7Qzey7NGFhHrgTNJyeNR4LqIWCjpfEnTACQdIGkF8EHgMkkLc9k1wOdJiWw+cH4eZ9ZQSjsnRP+aGoiIlfn/UknzgP2AxwYyQLOqRcQtwC01484tvJ5PqjudlZ0LzC01QLOSlXkk1OemhtzEsFl+vT1wCPBI96WGrtbWVlpbW6sOw8xsyCktCfWnqQHYE2iX9CvgTuArNVfVmZlZEyizOa7PTQ0RcRewT5mxmZlZ9YbyhQlmZtbknITMzKwyTkJmZlYZJyEzM6uMk5CZmVXGScjMGorvu2suTkJmZlYZJyEzM6uMk5CZmVXGScjMzCrjJGRmZpVxEjIzs8o4CZmZWWWchMzMrDJOQmZmVplSk5CkqZIWSVoi6ZxOph8q6X5J6yUdWzNthqTF+W9GmXGamVk1SktCkkYAs4GjgMnACZIm18z2JPBh4KqastsB5wEHAQcC50kaVVasZmZWjTKPhA4ElkTE0ohYB1wDHF2cISKWRcSDwIaaskcCt0fEmohYC9wOTC0xVjMzq0CZSWgMsLwwvCKPG7CykmZKapfUvnr16j4HamZm1WjoCxMiYk5EtEREy+jRo6sOx8zMeqnMJLQSGFsY3jmPK7usmZk1iDKT0HxgkqSJkjYFpgNtdZa9DThC0qh8QcIReZyZmTWR0pJQRKwHziQlj0eB6yJioaTzJU0DkHSApBXAB4HLJC3MZdcAnyclsvnA+XmcmZk1kU3KXHhE3ALcUjPu3MLr+aSmts7KzgXmlhmfmZlVq6EvTDAzs8bmJDTEtLa20traWnUYZmaDwknIzMwq4yRkVqE6nq+4maRr8/R7JU3I4ydIelnSgvx36WDHbjYQSr0wwcy6Vni+4ntJTwWZL6ktIh4pzHYqsDYidpU0HfgqcHye9lhETBnUoM0GmI+EzKrT4/MV8/AV+fUNwJ9L0iDGaFYqJyGz6tTzjMTX58n33v0OeFueNlHSA5J+Jundna3Az1e0oc5JyKwxrQIATfWyAAAFm0lEQVTGRcR+wCzgKklb187k5yvaUOckZFadep6R+Po8kjYBtgGei4hXIuI5gIi4D3gM2K30iBuMb3kY+pyEzKpTz/MV24COnoWPBe6IiJA0Ol/YgKRdgEnA0kGK22zA+Oo4s4pExHpJHc9XHAHM7Xi+ItAeEW3A5cB3JC0B1pASFcChwPmSXiV1Cnman69ojchJyKxCdTxf8Q+kB/zWlvse8L3SAzQrmZvjzMysMk5CZmZWGSchMzOrjJOQmZlVptQk5IczmplZd0q7Os4PZzQzs56UeSTkhzOWzHeDm1mjKzMJ+eGMZmbWraF6YYIfzmhmNgyUmYT8cEYzM+tWmUnID2c0M7NulXZ1nB/OaGZmPSn1AaZ+OKOZmXVnqF6YYGY2qHzLQzWchMzMrDJOQmZmVhknITMzq4yT0DDh9m4zG4qchMzMrDJOQmZmVhknITMzq4yTkJmZVcZJyMysl3yhz8BxErI3cQUzs8HiJGRmZpVxEjIzs8o4CVm/uOnOzPrDScjMrETeUeuek5CZmVWm1CQkaaqkRZKWSDqnk+mbSbo2T79X0oTCtE/m8YskHVlmnDY4vEf4Zq4jNtyVloQkjQBmA0cBk4ETJE2ume1UYG1E7ApcBHw1l51M6up7L2AqcEleng0TwyFhuY5YreGw3dcq80joQGBJRCyNiHXANcDRNfMcDVyRX98A/Lkk5fHXRMQrEfE4sCQvz+xNGrjiuo5Yn/Vmux/KdUQRUc6CpWOBqRHxd3n4Q8BBEXFmYZ6H8zwr8vBjwEHAZ4F7IuLKPP5y4NaIuKFmHTOBmQDjxo3b/4knnijlvVjjknRfRLRUHUdnXEdsKKi6jjT0hQkRMSciWiKiZfTo0VWHYzbkuI7YUFdmEloJjC0M75zHdTqPpE2AbYDn6ixr1uhcR2zYKzMJzQcmSZooaVPSSdS2mnnagBn59bHAHZHaB9uA6fnKoInAJOCXJcZqVgXXERv2NilrwRGxXtKZwG3ACGBuRCyUdD7QHhFtwOXAdyQtAdaQKiF5vuuAR4D1wBkR8VpZsZpVwXXErMQLEwZbS0tLtLe3Vx2GDTFVn3QdSlxHrDNV15GGvjDBzMwam5OQmZlVxknIzMwq4yRkZmaVaZoLEyStBrq6HXx74NlBDGcw+b11b3xE+C5Nuq0jzbwNQXO/v4avI02ThLojqb1Zr5Dye7P+avbPuZnfXzO8NzfHmZlZZZyEzMysMsMlCc2pOoAS+b1ZfzX759zM76/h39uwOCdkZmZD03A5EjIzsyHIScjMzCrT1ElI0lRJiyQtkXRO1fEMNEnLJD0kaYGkhn4ypaS5kp7JPYl2jNtO0u2SFuf/o6qMsRk1cx1ppvoBzVtHmjYJSRoBzAaOAiYDJ0iaXG1UpfhfETGl0e8VAL4NTK0Zdw7w04iYBPw0D9sAGSZ1pFnqBzRpHWnaJAQcCCyJiKURsQ64Bji64pisCxHxX6T+coqOBq7Ir68AjhnUoJqf60gDadY60sxJaAywvDC8Io9rJgH8RNJ9kmZWHUwJdoiIVfn1b4EdqgymCTV7HWn2+gFNUEdK61nVBsW7ImKlpLcDt0v6dd5bajoREZJ8P4H1xrCpH9C4daSZj4RWAmMLwzvncU0jIlbm/88APyA1rzSTpyXtCJD/P1NxPM2mqevIMKgf0AR1pJmT0HxgkqSJkjYFpgNtFcc0YCRtKWmrjtfAEcDD3ZdqOG3AjPx6BnBThbE0o6atI8OkfkAT1JGmbY6LiPWSzgRuA0YAcyNiYcVhDaQdgB9IgvQ9XhURP642pL6TdDXQCmwvaQVwHvAV4DpJp5K6IDiuugibT5PXkaaqH9C8dcSP7TEzs8o0c3OcmZkNcU5CZmZWGSchMzOrjJOQmZlVxknIzMwq4yRkZmaVcRIyM7PK/H/MW2tV/LOWPwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Density function plots\n",
    "\n",
    "# Density functions\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.subplots_adjust(wspace=0.5,hspace=0.5)\n",
    "\n",
    "k= np.arange(12)\n",
    "# Geometric\n",
    "p=0.4\n",
    "rv = stats.geom(p)\n",
    "a1 = fig.add_subplot(1,2,1)\n",
    "\n",
    "_ = a1.vlines(k,0,rv.pmf(k))\n",
    "_ = a1.set(ylabel=\"pmf\",title=\"Geometric Distribution\")\n",
    "_ = a1.text(0.5,0.8,\"p=0.4\", transform=a1.transAxes)\n",
    "\n",
    "# Negative Binomial\n",
    "r,p = 3, 0.4 \n",
    "rv = stats.nbinom(r,1-p)   # Note: scipy.stats definition is for stopping criterion\n",
    "\n",
    "a2 = fig.add_subplot(1,2,2)\n",
    "_ = a2.vlines(k,0,rv.pmf(x))\n",
    "_ = a2.set(ylabel=\"pmf\",title=\"Negative Binomial Distribution\")\n",
    "_ = a2.text(0.5,0.8,\"r=3, p=0.4\", transform=a2.transAxes)"
   ]
  }
 ],
 "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.8.2"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}