29-Solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Moldy Gelatin\n",
    "\n",
    "- Airborne spores produce tiny mold colonie on gelatin plates in a laboratory\n",
    "- The many plates average 3 colonies per plate\n",
    "- *What fraction of plates has exactly 3 colonies?*\n",
    "    - *If the average is a large integer $m$, what fraction of plates has exactly $m$ colonies?*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_____\n",
    "\n",
    "- At first, thinking about this problem, I was puzzled\n",
    "    - We know that the expected number of colonies is 3, but the question is like one of those classic Binomial distribution problems (e.g. \"What's the probability of rolling exactly one 6 in six rolls of a die?\")\n",
    "        - To calculate this, we need some value for $p$\n",
    "            - But we can't make $p=3$ since that wouldn't make sense\n",
    "            \n",
    "- The key idea is to split the plate into sections\n",
    "    - This eventually leads to using the Poisson distribution to solve for our answer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- If we divide our plate into 4 sections, we can treat each section having a colony as a Bernoulli trial\n",
    "    - Since the expected number of colonies on a single plate is equal to 3, the probability of any of the four sections having a colony is equal to 3/4\n",
    "        - So, following this logic, the probability of having exactly 3 colonies is equal to:\n",
    "        \n",
    "$$\n",
    "P(\\text{3 colonies}) = \\binom{4}{3}(3/4)^{3}(1/4)^{1} = 4\\cdot\\frac{3^{3}}{4^{4}} = \\left (\\frac{3}{4}\\right )^{3} = 0.421875\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- But this can't be the complete answer\n",
    "    - For one thing, modeling each section of the plate as a Bernoulli trial assume the area can only have 0 or 1 colony per section\n",
    "        - So, under this model, the plate can only have 0, 1, 2, 3, or 4 colonies\n",
    "            - We know that it's not super likely, but we should be able to have way more colonies (maybe even infinitely many)\n",
    "                - This idea will lead to us eventually using the Poisson distribution (described later)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- As a test, let's try dividing the plate into 100 evenly sized sections\n",
    "    - Here, for our Binomial distribution, $p=3/100$ so:\n",
    "    \n",
    "$$\n",
    "P(\\text{3 colonies}) = \\binom{100}{3}\\cdot0.03^{3}\\cdot0.97^{97}\n",
    "$$\n",
    "- We won't calculate this by hand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import binom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2274741274821556"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "binom.pmf(3, 100, 0.03)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\implies P(\\text{3 colonies}) \\approx 0.2275\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Let's plot these probabilities for increasing $n$ to see what it converges to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = 10, 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "array_n = np.arange(4, 1000)\n",
    "list_P = []\n",
    "\n",
    "for n in array_n:\n",
    "    p = 3/n\n",
    "    P = binom.pmf(3, n, p)\n",
    "    list_P.append(P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Probability')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAJQCAYAAADhfTKrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xu8XWV97/vPb851yZ3EEOQSYmLNRkOQACm4d7cibgQsrbCtWhSqVizbfaR66rYH3K16SvXUbXfRbTenVlus9QIqVo0WD60WFW+UpKRKoEhEJIsESLiTkMta63f+mGOtzKwsslbCGGOumXzer9d6zTme8YxnPnPO1wxfnmc8Y0RmIkmSpO7S6HQHJEmStP8McZIkSV3IECdJktSFDHGSJEldyBAnSZLUhQxxkiRJXcgQJ0mS1IUMcZIkSV3IECdJktSFejrdgTocfvjhuXjx4k53Q5IkaUJr1qzZkpkLJqp3SIS4xYsXs3r16k53Q5IkaUIR8YvJ1HM6VZIkqQsZ4iRJkrqQIU6SJKkLHRLnxEmSdKjbtWsXAwMDbN++vdNdUWHatGksXLiQ3t7eAzreECdJ0iFgYGCA2bNns3jxYiKi09055GUmDz30EAMDAyxZsuSA2nA6VZKkQ8D27duZP3++AW6KiAjmz5//jEZGDXGSJB0iDHBTyzP9PgxxkiRJXcgQJ0mSatFsNlmxYgXLly/nNa95Ddu2bQPgqaee4vTTT2doaIhf/OIXnHLKKaxYsYLjjz+ej33sY+O2tWvXLi6//HKWLl3K8uXLOfXUU/nGN76xz9dfvHgxW7ZsOaC+v/e97+Wb3/zmfh/39a9/nfe9730H9JoTMcRJkqRaTJ8+nbVr13LbbbfR19c3GtCuvvpqXvWqV9FsNjnqqKP4wQ9+wNq1a7n55pv54Ac/yMaNG/dq6z3veQ+bNm3itttu47bbbuNrX/saTzzxRGV9v+KKKzjzzDP3+7hzzz2XVatWjQbWMhniJElS7V784hezfv16AD772c9y3nnnAdDX10d/fz8AO3bsYHh4eK9jt23bxic+8Qn+/M//fLTus5/9bF772tcCcM0113DCCSewfPlyLrvssnFf/8orr2T58uUsX76cj3zkIwDcc889vOAFL+B3fud3OP744znrrLN46qmnAHjTm97EddddB8CaNWs4/fTTOeWUUzj77LPZtGkTAB/96EdZtmwZL3zhC7nggguA1nlvL33pS/n617/+zD+0MbzEiCRJh5g/+to6bt/4eKltLjt6Du/79eMnVXdwcJBvfOMbnHPOOezcuZO7776bxYsXj+7fsGED5557LuvXr+dP//RPOfroo/c4fv369SxatIg5c+bs1fbGjRu57LLLWLNmDfPmzeOss87iK1/5Cueff/5onTVr1vDJT36Sm2++mczktNNO4/TTT2fevHncddddXHPNNXziE5/gta99LV/60pe46KKLRo/dtWsXv/u7v8tXv/pVFixYwOc//3n+4A/+gKuvvpoPfvCD/PznP6e/v59HH3109JiVK1dy0003jYbMslQ6EhcR50TEnRGxPiIu30e9V0dERsTKYvvlEbEmIn5SPL6sre63izbXFn9HVPkeJElSOZ566ilWrFjBypUrWbRoERdffDFbtmxh7ty5e9Q79thj+fGPf8z69ev51Kc+xQMPPDDp17jlllt46UtfyoIFC+jp6eHCCy/ku9/97h51vve97/Gf//N/ZubMmcyaNYtXvepV3HTTTQAsWbKEFStWAHDKKadwzz337HHsnXfeyW233cbLX/5yVqxYwfvf/34GBgYAeOELX8iFF17IZz7zGXp6do+THXHEEeNOCT9TlY3ERUQTuAp4OTAA3BIRqzLz9jH1ZgNvB25uK94C/HpmboyI5cANwDFt+y/MzNVV9V2SpIPZZEfMyjZyTtzYsqe7VtrRRx/N8ccfz0033cSrX/3q0fLnPe953HvvvTzxxBPMnj17j2Myc8J+7KvOyPQstBZijEynth97/PHH88Mf/nCvY//+7/+e7373u6xatYo//uM/Zt26dfT09LB9+3amT58+Yb/2V5UjcacC6zPz7szcCVwLnDdOvT8GPgSMfoOZeWtmjkTWdcC0iOgf51hJktTF5s2bx9DQ0GiQGxgYGA1OjzzyCN///vc57rjj9jhmxowZXHzxxbz97W9n586dAGzatInPfOYznHbaaXznO99hy5YtDA0Ncc0113D66afvcfxLXvISvvKVr7Bt2za2bt3Kl7/8ZV784hdPqr/HHXccmzdvHg1xu3btYt26dQwPD7NhwwbOOOMMPvShD/Hoo4/y5JNPAvDTn/6U5cuXH/iH9DSqDHHHABvatgfYczSNiDgJODYz93W2328At2bmjrayTxZTqe+Jp7lSXkRcEhGrI2L15s2bD/AtSJKkqp111ll873vfA+COO+7gtNNO48QTT+T000/nXe96FyeccMJex7z//e9nwYIFLFu2jOXLl3P++eezYMECjjrqKP7kT/6EM844gxNPPJGTTz55dNHEiJNPPpk3velNnHrqqZx22mm85S1v4aSTTppUX/v6+rjuuuu47LLLOPHEE1mxYgU/+MEPGBoa4qKLLuKEE07gpJNO4vd+7/dGp4lvvPFGzj333Gf4Ke0tJjPseEANR7wGODsz31Js/xZwamb+brHdAP4JeFNm3hMR3wbe1T5NGhHHA6uAszLzZ0XZMZl5XzEN+yXgM5n5t/vqy8qVK3P1amdfJUmHrjvuuIMXvOAFne7GuG699VauvPJKPv3pT3e6K6V74IEHeP3rX8+3vvWtcfeP971ExJrMXDlR21WOxA0Ax7ZtLwTaz+qbDSwHvh0R9wAvAla1LW5YCHwZeMNIgAPIzPuKxyeAz9GatpUkSV3qpJNO4owzzmBoaKjTXSndvffey5/92Z9V0naVlxi5BVgaEUuA+4ALgNeP7MzMx4DDR7bbR+IiYi7w98C7M/P7bXV6gLmZuSUieoFfA/b/8smSJGlKefOb39zpLlTil3/5lytru7KRuMwcBC6ltbL0DuALmbkuIq6IiFdOcPilwPOA94y5lEg/cENE/BhYSyscfqKq9yBJ0sGkqlOodGCe6fdR2TlxU4nnxEmSDnU///nPmT17NvPnz+dp1gSqRpnJQw89xBNPPMGSJUv22DfZc+K8Y4MkSYeAhQsXMjAwgFdsmDqmTZvGwoULD/h4Q1wJrrpxPTesu59Vl/7HTndFkqRx9fb27jXio+5W6W23DhVbntzBzzdv7XQ3JEnSIcQQV4LAcwskSVK9DHElOfiXh0iSpKnEEFeCCJdtS5KkehniSuBkqiRJqpshriSOw0mSpDoZ4krQmk7tdC8kSdKhxBBXgoggHYuTJEk1MsSVwHPiJElS3QxxJXE6VZIk1ckQV4ZwYYMkSaqXIa4E3rFBkiTVzRBXFofiJElSjQxxJYjA1amSJKlWhrgSOJkqSZLqZogriatTJUlSnQxxJQhXp0qSpJoZ4krg6lRJklQ3Q1xJ0vlUSZJUI0NcCZxOlSRJdTPElcDJVEmSVDdDXEmcTZUkSXUyxJUhHIuTJEn1MsSVwAgnSZLqZogrkStUJUlSXQxxJRiZTTXDSZKkuhjiSjBysV8znCRJqoshrgSua5AkSXUzxJXIc+IkSVJdDHElGBmIM8JJkqS6GOJK4HSqJEmqmyGuRM6mSpKkuhjiShAxsjrVFCdJkuphiJMkSepChrgSOZ0qSZLqYogrgQsbJElS3QxxJQhMcZIkqV6GuBI5nSpJkupiiCvByHSqq1MlSVJdDHElcDJVkiTVzRBXIqdTJUlSXQxxJdg9nSpJklQPQ1wJXJ0qSZLqZogrUTqfKkmSamKIK4HTqZIkqW6GuBI5ECdJkupiiCtBeN8tSZJUM0NcmRyJkyRJNTHElWBkHM47NkiSpLoY4krgbKokSaqbIa5ELmyQJEl1qTTERcQ5EXFnRKyPiMv3Ue/VEZERsbKt7N3FcXdGxNn722addk+nSpIk1aOnqoYjoglcBbwcGABuiYhVmXn7mHqzgbcDN7eVLQMuAI4Hjga+GRH/rtg9YZt1c3WqJEmqW5UjcacC6zPz7szcCVwLnDdOvT8GPgRsbys7D7g2M3dk5s+B9UV7k22zI7xjgyRJqkuVIe4YYEPb9kBRNioiTgKOzcyvT/LYCdvsBO/YIEmS6lZliBtvjnE050REA/gw8N/249h9trlHAxGXRMTqiFi9efPmSXT3wDmZKkmS6lZliBsAjm3bXghsbNueDSwHvh0R9wAvAlYVixue7tiJ2hyVmR/PzJWZuXLBggXP8K1MjrOpkiSpLlWGuFuApRGxJCL6aC1UWDWyMzMfy8zDM3NxZi4GfgS8MjNXF/UuiIj+iFgCLAX+eaI2O6aYT/Viv5IkqS6VrU7NzMGIuBS4AWgCV2fmuoi4AlidmU8bvop6XwBuBwaBt2XmEMB4bVb1HibL6VRJklS3ykIcQGZeD1w/puy9T1P3pWO2PwB8YDJtThkOxEmSpJp4x4YSuDpVkiTVzRBXgnBCVZIk1cwQVyJXp0qSpLoY4kqwezrVFCdJkuphiCvByGSqI3GSJKkuhrgShKfESZKkmhniSuRAnCRJqoshrgQjq1PT+VRJklQTQ1wZnE6VJEk1M8SVyIE4SZJUF0NcCRyIkyRJdTPElSBcnipJkmpmiCuR06mSJKkuhrgSjF7s14uMSJKkmhjiSuBsqiRJqpshrkROp0qSpLoY4kowMhJnhpMkSXUxxJUgvMiIJEmqmSGuRN52S5Ik1cUQVwKnUyVJUt0McZIkSV3IEFciZ1MlSVJdDHEl2H3bLVOcJEmqhyGuBKMRzgwnSZJqYogrgXdskCRJdTPElciBOEmSVBdDXAlGLvbrdKokSaqLIa4ETqdKkqS6GeJKlE6oSpKkmhjiSuDqVEmSVDdDXAmcTpUkSXUzxJXIkThJklQXQ1wpitWpnhMnSZJqYogrgdOpkiSpboa4EjmdKkmS6mKIK4EDcZIkqW6GuBKE86mSJKlmhrgSOZ0qSZLqYogrwejFfl2dKkmSamKIK4GzqZIkqW6GuBI5nSpJkupiiCvByEicGU6SJNXFEFeC8CIjkiSpZoa4EqXzqZIkqSaGuDI4nSpJkmpmiCvB6CVGTHGSJKkmhrgSeMcGSZJUN0NcqRyKkyRJ9TDElcDpVEmSVDdDXAmcTZUkSXUzxJXIgThJklQXQ1wJRi7263SqJEmqiyGuBE6nSpKkuhniSuQdGyRJUl0qDXERcU5E3BkR6yPi8nH2vzUifhIRayPiexGxrCi/sCgb+RuOiBXFvm8XbY7sO6LK9zAZo6tTO9oLSZJ0KOmpquGIaAJXAS8HBoBbImJVZt7eVu1zmfmxov4rgSuBczLzs8Bni/ITgK9m5tq24y7MzNVV9X2/OZ0qSZJqVuVI3KnA+sy8OzN3AtcC57VXyMzH2zZnMv5g1uuAayrrZYmcTZUkSXWpbCQOOAbY0LY9AJw2tlJEvA14J9AHvGycdn6TMeEP+GREDAFfAt6f45yMFhGXAJcALFq06ED6P2mjq1OdUJUkSTWpciRuvEnGvVJOZl6Vmb8EXAb84R4NRJwGbMvM29qKL8zME4AXF3+/Nd6LZ+bHM3NlZq5csGDBgb6HSXF1qiRJqluVIW4AOLZteyGwcR/1rwXOH1N2AWOmUjPzvuLxCeBztKZtpwYH4iRJUk2qDHG3AEsjYklE9NEKZKvaK0TE0rbNc4G72vY1gNfQCncjZT0RcXjxvBf4NaB9lK4jXJ0qSZLqVtk5cZk5GBGXAjcATeDqzFwXEVcAqzNzFXBpRJwJ7AIeAd7Y1sRLgIHMvLutrB+4oQhwTeCbwCeqeg+TFc6nSpKkmlW5sIHMvB64fkzZe9uev2Mfx34beNGYsq3AKeX2sjyuTpUkSXXxjg0lGBmIc3WqJEmqiyGuBKPnxJnhJElSTQxxJfCUOEmSVDdDXIkciJMkSXUxxJWiuGOD86mSJKkmhrgSOJ0qSZLqZogrkeNwkiSpLoa4EowOxJniJElSTQxxJfCODZIkqW6GuBJ5sV9JklQXQ1wJvNivJEmqmyGuBM6mSpKkuhniSuRInCRJqoshrgQxcrHfDvdDkiQdOgxxJXA6VZIk1c0QVyJvuyVJkupiiCuREU6SJNXFEFcCp1MlSVLdDHElcjZVkiTVxRBXgth9ud+O9kOSJB06DHElGJlOdSROkiTVxRBXAs+JkyRJdTPElciBOEmSVBdDXAlG79hgipMkSTUxxJXA6VRJklQ3Q1yJ0glVSZJUE0NcCUYvMGKGkyRJNTHElcDpVEmSVDdDXIkciJMkSXUxxJViZHWqMU6SJNXDEFcCp1MlSVLdDHGSJEldyBBXAlenSpKkuhniShDOp0qSpJoZ4krkxX4lSVJdDHElcDpVkiTVzRBXAmdTJUlS3QxxJXIkTpIk1cUQV4IYudhvh/shSZIOHYa4EoxMp3rHBkmSVBdDnCRJUhcyxJXIcThJklQXQ1wJRlenmuIkSVJNDHEl8I4NkiSpboa4EnnHBkmSVBdDXAm8Y4MkSaqbIa4EzqZKkqS6GeJK5ECcJEmqiyGuBKN3bDDFSZKkmhjiSuB0qiRJqpshrkSuTpUkSXUxxJXA1amSJKluhrgyOJ0qSZJqVmmIi4hzIuLOiFgfEZePs/+tEfGTiFgbEd+LiGVF+eKIeKooXxsRH2s75pTimPUR8dGYQrdLcCBOkiTVpbIQFxFN4CrgFcAy4HUjIa3N5zLzhMxcAXwIuLJt388yc0Xx99a28r8ALgGWFn/nVPUeJmtkdarzqZIkqS5VjsSdCqzPzLszcydwLXBee4XMfLxtcyYTDGZFxFHAnMz8YWYm8LfA+eV2e/9NnbFASZJ0qKgyxB0DbGjbHijK9hARb4uIn9EaiXt7264lEXFrRHwnIl7c1ubARG12iuNwkiSpLlWGuPHGp/bKOZl5VWb+EnAZ8IdF8SZgUWaeBLwT+FxEzJlsmwARcUlErI6I1Zs3bz6gNzBZrk6VJEl1qzLEDQDHtm0vBDbuo/61FFOjmbkjMx8qnq8Bfgb8u6LNhZNpMzM/npkrM3PlggULDvhNTMbI2oo0xUmSpJpUGeJuAZZGxJKI6AMuAFa1V4iIpW2b5wJ3FeULioURRMRzaS1guDszNwFPRMSLilWpbwC+WuF7mBRPiZMkSXXrqarhzByMiEuBG4AmcHVmrouIK4DVmbkKuDQizgR2AY8AbywOfwlwRUQMAkPAWzPz4WLffwX+BpgOfKP4mxIch5MkSXWpLMQBZOb1wPVjyt7b9vwdT3Pcl4AvPc2+1cDyErv5jIVXGJEkSTXzjg0lCCdUJUlSzQxxJXIgTpIk1cUQV4bR6VRjnCRJqochrgTesUGSJNXNECdJktSFDHEl8I4NkiSpboa4EozescGlDZIkqSaGuBI0vE6cJEmqmSGuBI1iJG7YECdJkmpiiCvRsENxkiSpJoa4EjS8xogkSaqZIa4EIxlu2PlUSZJUE0NcCRqjq1MlSZLqYYgrwchkqufESZKkuhjiShBeYkSSJNXMEFeCiCAC0hQnSZJqYogrSeB14iRJUn0McSVpRHjbLUmSVBtDXEkiHImTJEn1McSVJCJc2CBJkmpjiCtJ4MIGSZJUH0NcSVrnxEmSJNXDEFeSRnjbLUmSVB9DXEkiwoUNkiSpNoa4kkTgJUYkSVJtDHElaS1s6HQvJEnSocIQV5JGI1ydKkmSamOIK4m33ZIkSXUyxJXE225JkqQ6GeJK4upUSZJUJ0NcSSK8Y4MkSaqPIa4kjXB1qiRJqo8hriRBMGyKkyRJNTHElcSROEmSVCdDXElc2CBJkupkiCuJt92SJEl1MsSVpBHhdKokSaqNIa4kEbiwQZIk1WZSIS4ivhQR50aEoe9pOBInSZLqNNlQ9hfA64G7IuKDEfH8CvvUlVr3TjXFSZKkekwqxGXmNzPzQuBk4B7gHyPiBxHx2xHRW2UHu0VrYYMkSVI9Jj09GhHzgTcBbwFuBf4XrVD3j5X0rMtEhLfdkiRJtemZTKWI+Dvg+cCngV/PzE3Frs9HxOqqOtdNvNivJEmq06RCHPBXmXl9e0FE9GfmjsxcWUG/uk4jvO2WJEmqz2SnU98/TtkPy+zIwcA7NkiSpLrscyQuIo4EjgGmR8RJtBZhAswBZlTct67iJUYkSVKdJppOPZvWYoaFwJVt5U8A/72iPnWlCFzYIEmSarPPEJeZnwI+FRG/kZlfqqlPXakR4SVGJElSbSaaTr0oMz8DLI6Id47dn5lXjnPYIcnbbkmSpDpNNJ06s3icVXVHul14TpwkSarRRNOpf1k8/lE93ele3nZLkiTVaaLp1I/ua39mvr3c7nQvL/YrSZLqNNF06ppaenEQaC1sMMVJkqR6TGZ1qiYhAoaHO90LSZJ0qNjnHRsi4iPF49ciYtXYv4kaj4hzIuLOiFgfEZePs/+tEfGTiFgbEd+LiGVF+csjYk2xb01EvKztmG8Xba4t/o7Y/7ddvnAkTpIk1Wii6dRPF4//c38bjogmcBXwcmAAuCUiVmXm7W3VPpeZHyvqv5LWBYXPAbYAv56ZGyNiOXADrTtHjLgwM1fvb5+q1FrY0OleSJKkQ8VE06lrisfvREQf8HwggTszc+cEbZ8KrM/MuwEi4lrgPGA0xGXm4231ZxZtk5m3tpWvA6ZFRH9m7pjUu+qARgRDpjhJklSTiUbiAIiIc4GPAT+jNei0JCL+S2Z+Yx+HHQNsaNseAE4bp+23Ae8E+oCXjd0P/AZw65gA98mIGAK+BLw/x7nfVURcAlwCsGjRon10sxxe7FeSJNVpn+fEtfkz4IzMfGlmng6cAXx4gmNinLK9Uk5mXpWZvwRcBvzhHg1EHA/8D+C/tBVfmJknAC8u/n5rvBfPzI9n5srMXLlgwYIJuvrMedstSZJUp8mGuAczc33b9t3AgxMcMwAc27a9ENi4j/rXAuePbETEQuDLwBsy82cj5Zl5X/H4BPA5WtO2HedInCRJqtNEF/t9VfF0XURcD3yB1mjaa4BbJmj7FmBpRCwB7gMuAF4/pv2lmXlXsXkucFdRPhf4e+Ddmfn9tvo9wNzM3BIRvcCvAd+c8F3WICJc2CBJkmoz0Tlxv972/AHg9OL5ZmDevg7MzMGIuJTWytImcHVmrouIK4DVmbkKuDQizgR2AY8AbywOvxR4HvCeiHhPUXYWsBW4oQhwTVoB7hMTv83qNQJv2SBJkmoz0erU334mjWfm9cD1Y8re2/b8HU9z3PuB9z9Ns6c8kz5VxUuMSJKkOk12deo04GLgeGDaSHlmvrmifnUdb7slSZLqNNmFDZ8GjgTOBr5Da5HCE1V1qht52y1JklSnyYa452Xme4Ctxf1UzwVOqK5b3Se8xIgkSarRZEPcruLx0eI2WIcBiyvpUZdqBIxzzWFJkqRKTOqcOODjETEPeA+wCphVPFchCK8TJ0mSajOpEJeZf1U8/Q7w3Oq6070aDa8wIkmS6jOp6dSImB8Rfx4R/xIRayLiIxExv+rOdRNH4iRJUp0me07ctbRus/UbwKuBLcDnq+pUN4oY58awkiRJFZnsOXHPysw/btt+f0Sc/7S1D0ER4XSqJEmqzWRH4m6MiAsiolH8vZbWvU1VcHWqJEmq0z5H4iLiCVqzhAG8E/hMsasBPAm8r9LedZFGhLfdkiRJtZno3qmz6+pIt2vdO9UUJ0mS6jHZc+KIiFcCLyk2v52ZX6+mS93Jc+IkSVKdJnuJkQ8C7wBuL/7eUZSpEJ4TJ0mSajTZkbhfBVZk5jBARHwKuBW4vKqOdZuGlxiRJEk1muzqVIC5bc8PK7sj3c6L/UqSpDpNdiTuT4BbI+JGWufwvwR4d2W96kLedkuSJNVpwhAXEQF8D3gR8Mu0QtxlmXl/xX3rKuElRiRJUo0mDHGZmRHxlcw8BVhVQ5+6UuDCBkmSVJ/JnhP3o4j45Up70uUaES5skCRJtZnsOXFnAG+NiHuArYwOPOULq+pYt4nwYr+SJKk+kw1xr6i0FweBhhf7lSRJNZro3qnTgLcCzwN+Avx1Zg7W0bFu5EicJEmqy0TnxH0KWEkrwL0C+LPKe9SlGuHVfiVJUn0mmk5dlpknAETEXwP/XH2XulPDc+IkSVKNJhqJ2zXyxGnUfWstbOh0LyRJ0qFiopG4EyPi8eJ5ANOL7ZHVqXMq7V0XaTS87ZYkSarPPkNcZjbr6ki3a4QhTpIk1WeyF/vVBJoRDDmfKkmSamKIK0lrOtVbb0mSpHoY4krSjABc3CBJkuphiCtJs/gkPS9OkiTVwRBXkihG4jwvTpIk1cEQV5JmY2Q61RAnSZKqZ4grSdOROEmSVCNDXEkaIyNxwx3uiCRJOiQY4krSbGU4p1MlSVItDHElGRmJGzLESZKkGhjiStIYuU6c58RJkqQaGOJK0nQkTpIk1cgQVxJXp0qSpDoZ4koyck6cA3GSJKkOhriSFBnOkThJklQLQ1xJPCdOkiTVyRBXElenSpKkOhniSuJInCRJqpMhriQj58R52y1JklQHQ1xJRqdTHYmTJEk1MMSVZHQ61XPiJElSDQxxJfHeqZIkqU6GuJKM3LEhDXGSJKkGhriSNEZvu9XhjkiSpEOCIa4kjeKT9Jw4SZJUh0pDXEScExF3RsT6iLh8nP1vjYifRMTaiPheRCxr2/fu4rg7I+LsybbZKU1Xp0qSpBpVFuIioglcBbwCWAa8rj2kFT6XmSdk5grgQ8CVxbHLgAuA44FzgP83IpqTbLMjXJ0qSZLqVOVI3KnA+sy8OzN3AtcC57VXyMzH2zZnAiMJ6Dzg2szckZk/B9YX7U3YZqeMrE51JE6SJNWhp8K2jwE2tG0PAKeNrRQRbwPeCfQBL2s79kdjjj2meD5hm53gxX4lSVKdqhyJi3HK9ko4mXlVZv4ScBnwhxMcO6k2ASLikohYHRGrN2/ePMkuH7imq1MlSVKNqgxxA8CxbdsLgY37qH8tcP4Ex066zcz8eGauzMyVCxYs2M+u7z9Xp0qSpDpVGeJuAZZGxJKI6KO1UGFVe4WIWNq2eS5wV/E4ON5DAAAcTUlEQVR8FXBBRPRHxBJgKfDPk2mzU5qeEydJkmpU2TlxmTkYEZcCNwBN4OrMXBcRVwCrM3MVcGlEnAnsAh4B3lgcuy4ivgDcDgwCb8vMIYDx2qzqPewPLzEiSZLqVOXCBjLzeuD6MWXvbXv+jn0c+wHgA5NpcyqI8BIjkiSpPt6xoSROp0qSpDoZ4kri6lRJklQnQ1xJRlanDjudKkmSamCIK4nTqZIkqU6GuJKM3LFhyBAnSZJqYIgryehtt5xOlSRJNTDElWRkOtVLjEiSpDoY4koyujrVDCdJkmpgiCvJyOrU9Jw4SZJUA0NcSRresUGSJNXIEFeSkXPiBg1xkiSpBoa4kvQ0XJ0qSZLqY4griSNxkiSpToa4kkQEzUYwOOzNUyVJUvUMcSVqhThH4iRJUvUMcSXqbQRDXihOkiTVwBBXIkfiJElSXQxxJeppNjwnTpIk1cIQV6KeRnixX0mSVAtDXIl6GsGg58RJkqQaGOJK1Gx6TpwkSaqHIa5EPY2GIU6SJNXCEFei1jlxLmyQJEnVM8SVqOk5cZIkqSaGuBL1eE6cJEmqiSGuRJ4TJ0mS6mKIK5HnxEmSpLoY4krUbAS7PCdOkiTVwBBXop6md2yQJEn1MMSVyHPiJElSXQxxJfKcOEmSVBdDXIm8TpwkSaqLIa5EvU2nUyVJUj0McSVqNlzYIEmS6mGIK1FPI9g15DlxkiSpeoa4EjkSJ0mS6mKIK1GP58RJkqSaGOJK1ONInCRJqokhrkRNz4mTJEk1McSVqNfbbkmSpJoY4krUbDS82K8kSaqFIa5Evc1g1/AwmQY5SZJULUNcifqaDTJxSlWSJFXOEFei3p7Wx7nLKVVJklQxQ1yJeputj3OnK1QlSVLFDHEl6msGgJcZkSRJlTPElWh0JG7QECdJkqpliCvRSIhzJE6SJFXNEFeivh5DnCRJqochrkS7p1NdnSpJkqpliCtRX48LGyRJUj0McSXynDhJklQXQ1yJXJ0qSZLqYogrkRf7lSRJdak0xEXEORFxZ0Ssj4jLx9n/zoi4PSJ+HBHfiojnFOVnRMTatr/tEXF+se9vIuLnbftWVPke9kdf09tuSZKkevRU1XBENIGrgJcDA8AtEbEqM29vq3YrsDIzt0XEfwU+BPxmZt4IrCjaeRawHviHtuN+PzOvq6rvB8pLjEiSpLpUORJ3KrA+M+/OzJ3AtcB57RUy88bM3FZs/ghYOE47rwa+0VZvyur1tluSJKkmVYa4Y4ANbdsDRdnTuRj4xjjlFwDXjCn7QDEF++GI6B+vsYi4JCJWR8TqzZs370+/D5gLGyRJUl2qDHExTtm4J4tFxEXASuBPx5QfBZwA3NBW/G7g+cAvA88CLhuvzcz8eGauzMyVCxYs2P/eH4Dd06meEydJkqpVZYgbAI5t214IbBxbKSLOBP4AeGVm7hiz+7XAlzNz10hBZm7Klh3AJ2lN204Ju0fihjrcE0mSdLCrMsTdAiyNiCUR0UdrWnRVe4WIOAn4S1oB7sFx2ngdY6ZSi9E5IiKA84HbKuj7Adl9TpwjcZIkqVqVrU7NzMGIuJTWVGgTuDoz10XEFcDqzFxFa/p0FvDFVibj3sx8JUBELKY1kvedMU1/NiIW0JquXQu8tar3sL9GplO9TpwkSapaZSEOIDOvB64fU/betudn7uPYexhnIURmvqzELpaqt9EKcTtc2CBJkirmHRtK1GgEfT0NV6dKkqTKGeJK1t/TYPsuFzZIkqRqGeJK1t/TdDpVkiRVzhBXsmm9DXZ4iRFJklQxQ1zJ+nsajsRJkqTKGeJK1t/TZIfnxEmSpIoZ4krWmk51JE6SJFXLEFey1kicIU6SJFXLEFeyfhc2SJKkGhjiSta6TpwjcZIkqVqGuJJN6206EidJkipniCuZlxiRJEl1MMSVzDs2SJKkOhjiSjat13unSpKk6hniStbf0zTESZKkyhniSja9r8lwwk6nVCVJUoUMcSWb2dcEYNvOwQ73RJIkHcwMcSWb0d8DwNadTqlKkqTqGOJKNmNkJG6HI3GSJKk6hriSzexzJE6SJFXPEFcyR+IkSVIdDHElm+k5cZIkqQaGuJLNcHWqJEmqgSGuZKMjcTsciZMkSdUxxJXMkThJklQHQ1zJZvQ5EidJkqpniCtZsxFM6204EidJkipliKvAzL4ethriJElShQxxFZjR32Sb06mSJKlChrgKOBInSZKqZoirwIy+Jtu82K8kSaqQIa4CM/t72OpttyRJUoUMcRVwJE6SJFXNEFcBz4mTJElVM8RVwNWpkiSpaoa4CjgSJ0mSqmaIq8CMvh627xpmcGi4012RJEkHKUNcBeZMb90/9YntjsZJkqRqGOIqMHdGLwCPPrWrwz2RJEkHK0NcBeZO7wPg0W07O9wTSZJ0sDLEVeAwR+IkSVLFDHEVmDu9FeIe22aIkyRJ1TDEVWDeDKdTJUlStQxxFZhTjMQ94kicJEmqiCGuAs1GMGdaD495TpwkSaqIIa4ic2f0OZ0qSZIqY4iryNwZva5OlSRJlTHEVeSw6b086jlxkiSpIoa4isyd0ec5cZIkqTKGuIrMnd7rOXGSJKkyhriKzJ3Ry2NP7WJ4ODvdFUmSdBAyxFVk3ow+hhMecTROkiRVwBBXkWfPmQbAg0/s6HBPJEnSwajSEBcR50TEnRGxPiIuH2f/OyPi9oj4cUR8KyKe07ZvKCLWFn+r2sqXRMTNEXFXRHw+IvqqfA8H6tlz+gFDnCRJqkZlIS4imsBVwCuAZcDrImLZmGq3Aisz84XAdcCH2vY9lZkrir9XtpX/D+DDmbkUeAS4uKr38EwcMbsYiXt8e4d7IkmSDkZVjsSdCqzPzLszcydwLXBee4XMvDEztxWbPwIW7qvBiAjgZbQCH8CngPNL7XVJjnAkTpIkVajKEHcMsKFte6AoezoXA99o254WEasj4kcRMRLU5gOPZubgJNvsmGm9TeZM63EkTpIkVaKnwrZjnLJxr7cRERcBK4HT24oXZebGiHgu8E8R8RPg8f1o8xLgEoBFixbtT79L8+w503jgcUfiJElS+aociRsAjm3bXghsHFspIs4E/gB4ZWaOJp7M3Fg83g18GzgJ2ALMjYiR8Dlum8VxH8/MlZm5csGCBc/83RyAI+b08+ATjsRJkqTyVRnibgGWFqtJ+4ALgFXtFSLiJOAvaQW4B9vK50VEf/H8cOBXgNszM4EbgVcXVd8IfLXC9/CMPHu2I3GSJKkalYW44ry1S4EbgDuAL2Tmuoi4IiJGVpv+KTAL+OKYS4m8AFgdEf9KK7R9MDNvL/ZdBrwzItbTOkfur6t6D8/Ugjn9bH5iB63sKUmSVJ4qz4kjM68Hrh9T9t6252c+zXE/AE54mn1301r5OuUdOWcaO4eGeWjrTg6f1d/p7kiSpIOId2yo0KJnzQBgw8PbJqgpSZK0fwxxFRoJcfca4iRJUskMcRVaOM+ROEmSVA1DXIWm9zU5Yna/I3GSJKl0hriKLXrWDEOcJEkqnSGuYoueNYMNDz/V6W5IkqSDjCGuYovmz2DjY0+xY3Co012RJEkHEUNcxZYcPpNMuGeLU6qSJKk8hriKHXfkbADufOCJDvdEkiQdTAxxFXvu4bPoaQQ/vd8QJ0mSymOIq1hfT4Mlh890JE6SJJXKEFeDf3fkbH5qiJMkSSUyxNXguGfP5t6Ht7Ft52CnuyJJkg4ShrgaLDtqDplw232Pd7orkiTpIGGIq8GKRXMBWLvhkQ73RJIkHSwMcTU4fFY/C+dNZ+2GRzvdFUmSdJAwxNXkpEXzWHuvIU6SJJXDEFeTFcfOZeNj27n/se2d7ookSToIGOJqctqSZwHww7u3dLgnkiTpYGCIq8myo+Ywb0Yv37vroU53RZIkHQQMcTVpNIL/8EuH8/31W8jMTndHkiR1OUNcjX7leYdz/+PbuevBJzvdFUmS1OUMcTV62fOPAOAf1t3f4Z5IkqRuZ4ir0ZGHTWPFsXO5Yd0Dne6KJEnqcoa4mp19/JH85L7H2PDwtk53RZIkdTFDXM1+/cSjiIAvrhnodFckSVIXM8TVbOG8Gbx46QK+uHoDQ8OuUpUkSQfGENcBv7nyWDY9tp2b7trc6a5IkqQuZYjrgDOXHcGzZvbxmR/d2+muSJKkLmWI64D+nia/9aLn8M07HuCOTY93ujuSJKkLGeI65M2/soRZ/T387xvXd7orkiSpCxniOuSwGb284d8/h+t/sol/u9/ROEmStH8McR30Oy9+LodN7+X/XrXO+6lKkqT9YojroHkz+3jXWcfxo7sf5ms/3tTp7kiSpC5iiOuw1526iBOOOYw/WrWOB5/Y3unuSJKkLmGI67BmI7jytSfy5I5Bfv+LP2bYCwBLkqRJMMRNAUufPZs//LVlfOenm/noP93V6e5IkqQu0NPpDqjlotMWsfbeR/nIN+9iyeEzOW/FMZ3ukiRJmsIMcVNERPD/vGo5Gx7Zxn/7wr/S12zwihOO6nS3JEnSFOV06hTS39Pkr9+4khOPncul19zK1/51Y6e7JEmSpihD3BQze1ovn3rzqZy8aC6/e82t/O9/ustryEmSpL0Y4qagWf09fPri0zh/xdH8z3/4Kf/HZ/+FR7ft7HS3JEnSFGKIm6Km9Tb58G+u4L//6vP5x9sf4OyPfJdv3/lgp7slSZKmCEPcFBYRXPKSX+Irb/sV5kzr5U2fvIX/8unV3PvQtk53TZIkdZghrgssP+Ywvva7/5HfP/s4brprC2de+R3+6Gvr2PTYU53umiRJ6pA4FE6aX7lyZa5evbrT3SjFA49v58/+4U7+7l/uoxHBq04+hjf8+8UsO3pOp7smSZJKEBFrMnPlhPUMcd1pw8Pb+Nh3fsZ1awbYMTjMSYvm8rpTF3H28Udy2PTeTndPkiQdIENcm4MxxI14bNsuvvQvA3z25l/ws81b6W0GL1m6gF894Sj+0wuOYO6Mvk53UZIk7QdDXJuDOcSNyEzWbniU63+yiet/cj/3PfoUjYATFs7lJUsP58VLF7Di2Ln09XgapCRJU5khrs2hEOLaZSb/OvAY377zQW66awtrNzzK0HDS39PghQsP46RF8zh50VxOWjSPZ8+Z1unuSpKkNoa4NodaiBvrsad28cOfPcQt9zzMv9z7COvue5ydQ8MAzJ/Zx/OPms1xz57D84+azfOPnM3SI2Yzva/Z4V5LknRommyI66mjM+qsw6b3cs7yIzln+ZEA7BgcYt3Gx1l776P82/2P82/3P8Hn/vkXbN81PHrMkXOm8Zz5M1g8fybPObx4nD+DhfNmMGdaDxHRqbcjSZIwxB2S+nuanLxoHicvmjdaNjSc/OKhrdx5/xPc9eCT3PPQVu59aBvf+rcH2fLkjj2On9HX5MjDpnH0YdM58rBpHHXYtNHHw2f1M39WP/Nn9jGt19E8SZKqYogTAM1G8NwFs3juglm8Ysy+J3cM8ouHtnLPlm1seuwpNj66nfsff4pNj23n++u38MDj2xkeZ1Z+Vn8P82f1MX9mH/Nn9XP4rD6eNbOPudP7OGx6L3Om9xaPPRxWPJ/V7yifJEmTUWmIi4hzgP8FNIG/yswPjtn/TuAtwCCwGXhzZv4iIlYAfwHMAYaAD2Tm54tj/gY4HXisaOZNmbm2yvdxqJvV38PxRx/G8UcfNu7+waFhNj+5g02PbefhJ3fy0NYdbHlyJw8Vzx96cicbHt7G2g2P8vDWnQyNl/gKzUYwZ1pbqJvWw8y+Hmb29zCzvzn6fEZfk1n9Pczo72FWf5MZfT3M6i/q9TWZ3tdkWm+T3qarcSVJB6fKQlxENIGrgJcDA8AtEbEqM29vq3YrsDIzt0XEfwU+BPwmsA14Q2beFRFHA2si4obMfLQ47vcz87qq+q7909NscNRh0znqsOkT1h0eTp7cOcjjT+3iseJvz+eDo88fe2oXT+4Y5KEnt7F15yDbdgzx5I5BdgwOT/g6o31rBNN6m0zrbdDf03psbTeZPlLe22TaHvsaxb4m/T0NepsN+toe+8Zs9zaD/p4Gfc0mvT1BX7NB70i9ZoNGw5FFSVL5qhyJOxVYn5l3A0TEtcB5wGiIy8wb2+r/CLioKP9pW52NEfEgsAB4FHW1RiOYM62XOdN6WThv4vrjGRwaZuvOIbbtHGTrjkG27hhqPe4ceWwFvu27htg+OMT2XcNs3zXEU7uG2FE83z7Y2n546062D7aVF/X2MVi433oasUcQ7O9pPe9pBM1iX+txz+2eRqtOTzOKx/bt4vhm0NtoP74x2s5o/bHHFH+NRtCM4nmMlDP6fHdZ2/MIGg3anrc9tu8v6js1LknVqTLEHQNsaNseAE7bR/2LgW+MLYyIU4E+4GdtxR+IiPcC3wIuz8wdY4/Twaun2eCw6Y3Kbi+WmewaSrYPDrFzcHj0b9fQMDuKx9Z2snNoiJ2Dyc6hYXYNDrcei/171i3aGcricZjBoWEGh5Oh4WTX0DBDw8ngULJ1cJDB4vng8PDo8/Z6o4/F8fuaou6kCPYIe83G2GDIOCEwaEQrTEa0P2eP7cZe2606jdjz+FbZvuuM7m+02gzGvkZRp7Hv44NJ1AkIGO1b0F62+/V3PzLaHmPLGekjwJ5tj5QHrcI9ytvqskc7e/aF9vK214x4+ueNcfoy8nmOZPr299n+mbBHO3u+Jm2fw8j/GozuG+3v+OW764851v/JUJerMsSN9+sY9780EXERsJLWuW7t5UcBnwbemJkjc2jvBu6nFew+DlwGXDFOm5cAlwAsWrTowN6BDkkRQV9PdNXdLYaHk6FsC35D2Qp/bc+HhlvBcyT0DWW2jht9zjhlrceh4WQ4k6Fhxikbs3+vsvbXgqHh4eL43W2NbTOz1c5wtkL18Jjt3ftbZUPDw3vVGS7e03BxLczhtnb2OH74aV5j5JjhPfuz1/FTMz9rP40GTHaHu9ijfHcaHFv+dMcyNkBO8Brs1da+X2OPfk8iuO7VH8Y/lqer8zSvwV5t7fs19uzHmDbHec/t9vh89/iM9r2/PZLsfczu19xje0z5yEG/f9ZxLD585t6d64AqQ9wAcGzb9kJg49hKEXEm8AfA6e0jahExB/h74A8z80cj5Zm5qXi6IyI+CbxrvBfPzI/TCnmsXLnSf2Z1UGs0ggZB66ouXtqlbjlOwBsb9MbWSYrH9ufsDqnjlbcCY3tZq87w8O66tJe3BdJsO3ZsedJqcHi8vuzVp3HK2bO/I8F57PscLo5jzOsMt7VDW7+Gc8++wNi293y/jL7P3d/L6GfStm+Pz6rtoH3VaS9nj/Kn78/Yvkz2NWgvH/PexnvP7FGeY+rsXc7YNifZn8m8BnuUT+I1hkdLx7z27s+rvc09ysbU2eM7at9u/x7G6es+j92jTmvrqV1DTBVVhrhbgKURsQS4D7gAeH17hYg4CfhL4JzMfLCtvA/4MvC3mfnFMccclZmbohW3zwduq/A9SNKEIoJmQHPcCQhJqkZlIS4zByPiUuAGWkMDV2fmuoi4AlidmauAPwVmAV8shkDvzcxXAq8FXgLMj4g3FU2OXErksxGxgNZI6FrgrVW9B0mSpKnKe6dKkiRNIZO9d2r3nLktSZKkUYY4SZKkLmSIkyRJ6kKGOEmSpC5kiJMkSepChjhJkqQuZIiTJEnqQoY4SZKkLmSIkyRJ6kKGOEmSpC5kiJMkSepChjhJkqQuZIiTJEnqQoY4SZKkLmSIkyRJ6kKGOEmSpC5kiJMkSepChjhJkqQuZIiTJEnqQoY4SZKkLmSIkyRJ6kKRmZ3uQ+UiYjPwi4qaPxzYUlHbOnB+L1OT38vU5PcyNfm9TE11fC/PycwFE1U6JEJclSJidWau7HQ/tCe/l6nJ72Vq8nuZmvxepqap9L04nSpJktSFDHGSJEldyBD3zH280x3QuPxepia/l6nJ72Vq8nuZmqbM9+I5cZIkSV3IkThJkqQuZIh7BiLinIi4MyLWR8Tlne7PoSIijo2IGyPijohYFxHvKMqfFRH/GBF3FY/zivKIiI8W39OPI+Lkzr6Dg1tENCPi1oj4erG9JCJuLr6Xz0dEX1HeX2yvL/Yv7mS/D2YRMTcirouIfyt+N//e30vnRcTvFf+G3RYR10TENH8vnRERV0fEgxFxW1vZfv9GIuKNRf27IuKNVffbEHeAIqIJXAW8AlgGvC4ilnW2V4eMQeC/ZeYLgBcBbys++8uBb2XmUuBbxTa0vqOlxd8lwF/U3+VDyjuAO9q2/wfw4eJ7eQS4uCi/GHgkM58HfLiop2r8L+D/y8znAyfS+n78vXRQRBwDvB1YmZnLgSZwAf5eOuVvgHPGlO3XbyQingW8DzgNOBV430jwq4oh7sCdCqzPzLszcydwLXBeh/t0SMjMTZn5L8XzJ2j9B+kYWp//p4pqnwLOL56fB/xttvwImBsRR9Xc7UNCRCwEzgX+qtgO4GXAdUWVsd/LyPd1HfCfivoqUUTMAV4C/DVAZu7MzEfx9zIV9ADTI6IHmAFswt9LR2Tmd4GHxxTv72/kbOAfM/PhzHwE+Ef2DoalMsQduGOADW3bA0WZalRMKZwE3Aw8OzM3QSvoAUcU1fyu6vMR4P8Chovt+cCjmTlYbLd/9qPfS7H/saK+yvVcYDPwyWKa+68iYib+XjoqM+8D/idwL63w9hiwBn8vU8n+/kZq/+0Y4g7ceP8H5FLfGkXELOBLwP+ZmY/vq+o4ZX5XJYuIXwMezMw17cXjVM1J7FN5eoCTgb/IzJOAreyeFhqP30sNimm284AlwNHATFrTdGP5e5l6nu67qP07MsQduAHg2LbthcDGDvXlkBMRvbQC3Gcz8++K4gdGpn2KxweLcr+revwK8MqIuIfW6QUvozUyN7eYLoI9P/vR76XYfxh7T2fomRsABjLz5mL7Olqhzt9LZ50J/DwzN2fmLuDvgP+Av5epZH9/I7X/dgxxB+4WYGmxkqiP1gmpqzrcp0NCcR7IXwN3ZOaVbbtWASOrgd4IfLWt/A3FiqIXAY+NDJGrPJn57sxcmJmLaf0e/ikzLwRuBF5dVBv7vYx8X68u6juyULLMvB/YEBHHFUX/Cbgdfy+ddi/wooiYUfybNvK9+HuZOvb3N3IDcFZEzCtGWs8qyirjxX6fgYj4VVojDU3g6sz8QIe7dEiIiP8I3AT8hN3nXv13WufFfQFYROsfyNdk5sPFP5D/m9YJptuA387M1bV3/BASES8F3pWZvxYRz6U1Mvcs4FbgoszcERHTgE/TOqfxYeCCzLy7U30+mEXEClqLTfqAu4HfpvU/8f5eOigi/gj4TVor7m8F3kLrHCp/LzWLiGuAlwKHAw/QWmX6FfbzNxIRb6b13yOAD2TmJyvttyFOkiSp+zidKkmS1IUMcZIkSV3IECdJktSFDHGSJEldyBAnSZLUhQxxkqT/vx07tG0ojKEw6iuVZIuQoo6T2YKySXeoVFIVZoXQgMghBaV55Jf1zkGGF34yMJCIAwAYSMQBvCjJMclvkkuSnySfSQ6rdwH7IuIAtnmvqnN3f1TVrapOi/cAOyPiALa5dvf33/1VVceFW4AdEnEA29z/3Y+qels1BNgnEQcAMJCIAwAYKN29egMAAC/yiQMAGEjEAQAMJOIAAAYScQAAA4k4AICBRBwAwEAiDgBgIBEHADDQE62ODEkz10VEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(array_n, list_P, label = 'P(3 Colonies)')\n",
    "plt.legend()\n",
    "plt.xlabel('n')\n",
    "plt.ylabel('Probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- As we can see, the probability converges pretty quickly to a value around 0.225"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Now, we introduce the Poisson distribution\n",
    "    - Consider a situation where we're using the Binomial distribution, but:\n",
    "        1. $p$ is very small\n",
    "        2. $n$ is very big\n",
    "        \n",
    "- In this scenario, the $\\binom{n}{k}$ value will be massive, and the $p^{k}$ value will be tiny\n",
    "    - We can simplify this calcuation by approximating it using the Poisson distribution\n",
    "    \n",
    "$$\n",
    "P(\\text{3 colonies}) = \\frac{\\lambda^{k}e^{-\\lambda}}{k!}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Here, $\\lambda$ is the expected number of colonies i.e. $\\lambda = 3$\n",
    "    - **Note**: $p$ gives us a probability where $\\lambda$ gives us an expected value\n",
    "        - So $p=1/\\lambda$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- So, using the Poisson distribution for our calculation:\n",
    "\n",
    "$$\n",
    "P(\\text{3 colonies}) = \\frac{3^{3}e^{-3}}{3!} \\approx 0.224042\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- And this value fits in perfectly with our plot above!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Now, expanding this solution where $m$ is very large, we can just replace the $lambda = 3$ with $\\lambda = m$"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}