week01-probability-and-distributions.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import math\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Focus questions\n",
    "\n",
    "Goals for student knowledge at the end of the lesson:\n",
    "\n",
    "1. Give examples of situations where the following distributions are an appropriate model for observed data:\n",
    "    - Poission\n",
    "    - Weibull\n",
    "    - log-normal\n",
    "2. Describe a situation where the central limit theorem can be applied to use statistics based on normal distributions, when population being sampled is non-normal.\n",
    "\n",
    "## Review\n",
    "\n",
    "#### Measurement types\n",
    "* Ordinal: Hurricane scale\n",
    "* Nominal: species classification\n",
    "* Scale\n",
    "  * Interval: Zero point is arbitrarily defined\n",
    "  * Ratio: Variables with a natural zero point, Kelvin scale, length, mass\n",
    "\n",
    "#### Error types\n",
    "\n",
    "Accuracy\n",
    "* closeness to a known value\n",
    "* close to other values (repeatability)\n",
    "    \n",
    "Measurement Errors\n",
    "* Systematic - bias, drift\n",
    "* Random - noise\n",
    "  \n",
    "Computational errors:\n",
    "* Numerical errors: approximations with limited computational resources\n",
    "* Round-off, truncation\n",
    "* Data errors: garbage in - garbage out\n",
    "* Errors from previous computations\n",
    "* Modeling errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probability Density functions\n",
    "\n",
    "__Probability Functions__\n",
    "\n",
    "$P(x)$ = probability mass function (for discrete data)\n",
    "\n",
    "_Example:_ The probability an number $x$ coming up on a dice role. The PMF is plotted as vertical lines because the probability of rolling a non-integer (e.g. 3.5) is equal to zero.\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmYHGW59/HvjwQS2QLCiJAEgoAcgyBCE1QUWZQTXIi+Rkj0KFGOyIW4vOKCvigQccENXDgqsoWwBIxyTjyigQOCwBHIJGyGGBliIENQBhOWoBhC7vePekYrne6u7slUemby+1zXXNNV9TxV91NV3XfXU9VVigjMzMwa2azdAZiZ2cDnZGFmZoWcLMzMrJCThZmZFXKyMDOzQk4WZmZWyMmin0kKSXv2se5SSW+qM+0NkhbXKivp85Iu7FvEA5syl0haKemujbzsX0o6vqR5HyXpP8uYd53l3SVpn4IyZ0t6QtKfNlZcfSVpmqTbcsN9et+1c/8qkm+TpB9K+kI743GyGCQi4taI2LvOtK9ExL8DSBqXdrLhGzfC0rweeDMwJiImlLUQSWdKujw/LiKOjogZJS3yK8DXcssPSc9KWiXpUUnfljQsTVsqabWkHativifVG5eGL03lVuX+jkvFvwlMrxeMpLHAqcD4iHhpv7Z0YNso+9eGioiTIuJL7YzByaIFQ+gDeDDZDVgaEc+2O5D+IukgYFRE3FE16VURsTVwJPAe4EO5aX8EpubmsS/wohqz/3pEbJ37uzqNnwMcLmnnOmHtBvwlIh7vQ3v6/X3Rmyg3giG3f5Vlk08W6Vvb5yQ9kA5FL5E0Mk07TFK3pM+mQ/NL0vgPSeqStELSHEm7VM32LZKWpEP6b0jaLNXbQ9JNkv6Spl0habuqugc1iqVOG/Lfin+T/j+Zvlm+McW5b678SyT9TVJHjXlNk3S7pHMlPZna8bo0fpmkx/NdM5LeKuluSU+n6Wfmpo2UdHlq75OS5knaKbecJZKekfRHSe+tEcsJwIXAa1Nbzqrufkjl8ofrl0o6X9Iv0rzvlLRHruw+km5I6+TPyrrwJgKfB45Ly7k3lb1ZUu8R22aSTpf0cFoHl0kalab1Hs0dL+mRtG3/X61tlRwN3FJvYkT8HrgVeGVu9Ezg/bnh44HLGiyjep7PAfOBo6qnKevOvAHYJbX/0jT+GEkL07a7WdIrcnWWpvfFfcCztRKGpO+kfeJpSfMlvaFefGm7/UDSdZKeJUtso9J67knr/fTe91IrJO2S3qcr0vv2Q2n8evtXjbp93u6p7mmSHkrvgWskvbhBnJ+W9Jik5ZI+WGP9nJ0bfpuyI8snJf2vpP1aXS8ti4hN+g9YCvwOGAu8GLgdODtNOwxYA5wDjCD7JncE8ARwQBr3PeA3ufkF8Os0r12BPwD/nqbtSXbIOwLoIPtgP6+FWLqryr4pvT4TuDy9HpdiGJ4r+x/AObnhjwM/r7M+pqU2fwAYBpwNPAKcn+I+CngG2DoX175kXzz2A/4MvCNN+zDwc2DLNK8DgW2BrYCngb1TuZ2BfRrEc1u94dw63zO9vhRYAUwAhgNXALPStG2Ax8i6W0am4YOr12Fuvjfntt0HgS7gZcDWwM+AmVXr/Mdk+8irgL8Dr6jTpp8An27QhvHAn4AT8tsaWAy8Iq3LZWTfigMYl2v72Q329e8C364z7TDW3b9eDjxLtr9uDnwmtX+LXEz3kO2rL6ozz38Ddkjb4dTUppF1yl4KPAUckvalkWTJ8L/SdhpH9l7qXSfV+8U/1l+Ned9C9h4YCewP9ABH1tufqur2ebsDnwDuAMaQvXd+BFxVZzkTyd47ryR7f1zJ+vt172fBAcDjwMFpXzg+bY8RpX5WljnzwfCXVvJJueG3AA/l3kCr8zs4cBHZoX7v8NbA87k3bAATc9NPBm6ss+x3AHe3EEtfk8XBZB8um6XhTuDYOjFNAx7MDe+b5rdTbtxfgP3r1D8PODe9/iDwv8B+VWW2Ap4E3kWdD5qqeFpNFhdWrcPfp9dT8+u7ah7/WIe5cTfzz2RxI3BybtreabsPz63zMbnpdwFT6izrhvx2zrXhaWAl8BBZku7dXkvJksXpwFfJPlhuSMuuThbPpXX7JPBE1TK+DFxcJ6bq/esLwDW54c2AR4HDcjF9sMX32kqyrrZa0y4FLssNDyP74B2fG/dh4OY6+0XNZEGWzF4AtsmN+ypwab39qap+n7c7sIiUlNLwzr11ayznYuBrueGXUz9Z/AD4UlX9xcAbW9kerf5t8t1QybLc64eBfLdST2SH8L12SWUAiIhVZB+eo4vmp6z7Z5ayE5hPA5cD65y0LIilTyLiTrJviW+U9C9kRzhzGlT5c+7139I8qsdtDSDpYEm/Tl0FTwEn8c82zQTmArPSofXXJW0eWf/wcansY6nL6F82tJ05+at5/tobK9kHx0N9nOc62z29Hg7s1MRyq60k+7Zc7YCI2D4i9oiI0yNibdX0mWTnMqZRvwvqmxGxXfqr3re2IUsizajez9eS7Zv19vP1SDpV0iJJT0l6EhjF+vt7Xn5+OwJbsP46H01rdgFWRMQzfZzPhmz33YBrU1fRk2TJ44WquvnlVL/369kNOLV3vmneY+mHz4pGnCwyY3OvdwWW54ajquxyso0FgKStyA61H21ifl9N89svIrYlO0xXC7E0ozreXjPS8t4HzK5KgBviSrLEMzYiRgE/JLUpIp6PiLMiYjzwOuBtpH73iJgbEW8m+7b1e7JD+WY8S9atBYCkVq7cWQbsUWdavfXWa53tTrZt1rBuYm3WfWTfHFsSEQ+Tneh+C1l3SKteAdzbZNnq/Vxk+2Z+P6+7ztL5ic8CxwLbR8R2ZN1M1ft7Xn5+T5B9C69e54/SmuXAiyXlk3Mr89mQ7b4MODqXvLeLiJERUWvZj7H+e7/RfL9cNd8tI+KqJmLqMyeLzEckjUknnz4PXN2g7JXAByTtL2kE2SWQd0bE0lyZT0vaXtnliB/PzW8bYBXZyefRwKc3MJZaeoC1ZH2seTOBd5IljKZPjDZhG7Jvbs9JmkD2zRcASYdL2lfZlS1Pk735X5C0Uzp5uhVZV8Mqsm9czbgX2Cet/5Fk3UfN+m/gpZI+IWmEpG0kHZym/RkY1+AE6lXA/5W0u6Stybb71RGxpoXl97oOeGMf6gGcABwRLV69k/bVA8m6r5pxDfBWSUdK2pzsnMPfyboVm7EN2YdqDzBc0hfJzlc1JSJeSDF8OW2n3YBPkh2NNy0ilqWYv6rsgov9yNbhFU3OYkO2+w9T/LsBSOqQNKlO2WuAaZLGS9oSOKPBfH8MnJSO6iVpK2UXmtQ6Wu03ThaZK4HrgSXp7+x6BSPiRrL+3J+SfRvYA5hSVey/yK48uQf4Bdl5DoCzyE5OPZXG1/p22HQsdeL7K1nf9O3pEPU1aXw3sIDs29utrcyzwMnAdEnPAF8k2+l7vRSYTZYoFpGdaLycbL87lexb2wqyD86Tm1lYRPyB7PcC/wM8CNzWuMY6dZ8hO2H7drKugweBw9Pkn6T/f5G0oEb1i8kS7m/Ivt0/B3y02WVXxbEAeCqXqFqp+1BEdPZhsceQ9fc3daQaEYvJvlh8j+xb/tuBt0fE6iaXNxf4JdlJ6YfJ1lfDbqsaPkp2JLmEbDtfSbYdWjWV7PzCcuBa4IyIaDZpbsh2/w7ZUff16f1xB9n5w/VExC/JzvfdRHZC/aZ6M03b/0PA98m6NLvIuiZLpXRyZJMlaSnZScz/aXcsZZN0MbA8Ik5vdyybOklHkZ04fcdGWt6dZFcS/W5jLM+GHv/IbBOh7Fe+/wd4dXsjMYCIuJ7sCHJjLa/loxizPHdDbQIkfYns9xvfiIg/tjseMxt8NvluKDMzK+YjCzMzKzRkzlnsuOOOMW7cuHaHYWY2qMyfP/+JiFjvPnHVhkyyGDduHJ2dfbmi0Mxs0yWp0a/F/8HdUGZmVsjJwszMCjlZmJlZIScLMzMr5GRhZmaFnCzMzKyQk4WZmRVysjAzs0JOFmZmVmjI/IJ7YzrzzNqvhwq3b3Bz+wavgdy2IXPX2UqlEhvrdh/KPUV4iKy+dbh9g5vbN3i1o22S5kdEpaicu6HMzKyQk4WZmRVysjAzs0KlJgtJEyUtltQl6bQa0w+VtEDSGkmTq6btKul6SYskPZCeIW1mZm1QWrKQNAw4HzgaGA9MlTS+qtgjwDTgyhqzuIzsmdGvACYAj5cVq5mZNVbmpbMTgK6IWAIgaRYwCXigt0BELE3T1uYrpqQyPCJuSOVWlRinmZkVKLMbajSwLDfcncY14+XAk5J+JuluSd9IRyrrkHSipE5JnT09Pf0QspmZ1VJmslCNcc1eOTwceAPwKeAg4GVk3VXrzizigoioRESlo6PwEbJmZtZHZSaLbmBsbngMsLyFundHxJKIWAP8J3BAP8dnZmZNKjNZzAP2krS7pC2AKcCcFupuL6n3cOEIcuc6zMxs4yotWaQjglOAucAi4JqIWChpuqRjACQdJKkbeDfwI0kLU90XyLqgbpR0P1mX1o/LitXMzBrzvaH6YCjfmwbcvsHO7Ru8fG8oMzMb1JwszMyskJOFmZkVcrIwM7NCThZmZlbIycLMzAo5WZiZWSEnCzMzK+RkYWZmhZwszMyskJOFmZkVcrIwM7NCThZmZlbIycLMzAo5WZiZWSEnCzMzK1RqspA0UdJiSV2STqsx/VBJCyStkTS5xvRtJT0q6ftlxmlmZo2VliwkDQPOB44GxgNTJY2vKvYIMA24ss5svgTcUlaMZmbWnDKPLCYAXRGxJCJWA7OASfkCEbE0Iu4D1lZXlnQgsBNwfYkxmplZE8pMFqOBZbnh7jSukKTNgG8Bny4od6KkTkmdPT09fQ7UzMwaKzNZqMa4Zh9BfjJwXUQsa1QoIi6IiEpEVDo6OloO0MzMmjO8xHl3A2Nzw2OA5U3WfS3wBkknA1sDW0haFRHrnSQ3M7PylZks5gF7SdodeBSYArynmYoR8d7e15KmARUnCjOz9imtGyoi1gCnAHOBRcA1EbFQ0nRJxwBIOkhSN/Bu4EeSFpYVj5mZ9Z0imj2NMLBVKpXo7OzcKMtS7mzMEFl963D7Bje3b/BqR9skzY+ISlE5/4LbzMwKOVmYmVkhJwszMyvkZGFmZoWcLMzMrJCThZmZFXKyMDOzQk4WZmZWyMnCzMwKOVmYmVkhJwszMyvkZGFmZoWcLMzMrJCThZmZFXKyMDOzQk4WZmZWqNRkIWmipMWSuiSt91hUSYdKWiBpjaTJufH7S/qtpIWS7pN0XJlxmplZY6UlC0nDgPOBo4HxwFRJ46uKPQJMA66sGv9X4P0RsQ8wEThP0nZlxWpmZo0NL3HeE4CuiFgCIGkWMAl4oLdARCxN09bmK0bEH3Kvl0t6HOgAniwxXjMzq6PMbqjRwLLccHca1xJJE4AtgIdqTDtRUqekzp6enj4HamZmjZWZLFRjXEuPIJe0MzAT+EBErK2eHhEXREQlIiodHR19DNPMzIqUmSy6gbG54THA8mYrS9oW+AVwekTc0c+xmZlZC8pMFvOAvSTtLmkLYAowp5mKqfy1wGUR8ZMSYzQzsyaUliwiYg1wCjAXWARcExELJU2XdAyApIMkdQPvBn4kaWGqfixwKDBN0j3pb/+yYjUzs8YU0dJphAGrUqlEZ2fnRlmWcmdjhsjqW4fbN7i5fYNXO9omaX5EVIrK+RfcZmZWyMnCzMwKOVmYmVkhJwszMyvkZGFmZoWcLMzMrJCThZmZFXKyMDOzQk4WZmZWyMnCzMwKNZUs0lPvzMxsE9XskUWXpG/UeCyqmZltAppNFvsBfwAulHRHekLdtiXGZWZmA0hTySIinomIH0fE64DPAGcAj0maIWnPUiM0M7O2a/qchaRjJF0LfAf4FvAy4OfAdSXGZ2ZmA0Cz3VAPApOAb0TEqyPi2xHx54iYDfyqXiVJEyUtltQl6bQa0w+VtEDSGkmTq6YdL+nB9Hd8K40yM7P+NbzJcu+PiNvyIyQdEhG3R8THalVIV1CdD7yZ7Hnc8yTNiYgHcsUeAaYBn6qq+2Kyrq4KEMD8VHdlk/GamVk/avbI4rs1xn2voM4EoCsilkTEamAW2dHJP0TE0oi4D1hbVfdfgRsiYkVKEDcAE5uM1czM+lnDIwtJrwVeB3RI+mRu0rZA0W8vRgPLcsPdwMFNxlWr7ugm65qZWT8r6obaAtg6ldsmN/5pYHLNGv+kGuOafapsU3UlnQicCLDrrrs2OWszM2tVw2QREbcAt0i6NCIebnHe3cDY3PAYYHkLdQ+rqntzjfguAC4AqFQqQ+zR7WZmA0dRN9R5EfEJ4PuS1vswjohjGlSfB+wlaXfgUWAK8J4m45oLfEXS9mn4KOBzTdY1M7N+VtQNNTP9/2arM46INZJOIfvgHwZcHBELJU0HOiNijqSDgGuB7YG3SzorIvaJiBWSvkSWcACmR8SKVmMwM7P+oYih0XtTqVSis7NzoyxLuTMqQ2T1rcPtG9zcvsGrHW2TND8iKkXlirqh7qfBSemI2K8PsZmZ2SBT1A31to0ShZmZDWhFV0O1egWUmZkNQQ1/wS3ptvT/GUlPV//fOCGamVm7FR1ZvD7936ZROTMzG9qavZEgkg4AXk92wvu2iLi7tKjMzGxAafZ5Fl8EZgA7ADsCl0o6vczAzMxs4Gj2yGIq8OqIeA5A0teABcDZZQVmZmYDR7O3KF8KjMwNjwAe6vdozMxsQCr6Ud73yM5R/B1YKOmGNPxm4LZGdc3MbOgo6obqvX/GfLJ7OPW6uZRozMxsQCq6dHbGxgrEzMwGrqZOcEvaC/gqMJ7cuYuIeFlJcZmZ2QDS7AnuS4AfAGuAw4HL+Ofty83MbIhrNlm8KCJuJLul+cMRcSZwRHlhmZnZQNLs7yyek7QZ8GB6oNGjwEvKC8vMzAaSZo8sPgFsCXwMOBB4H3B8USVJEyUtltQl6bQa00dIujpNv1PSuDR+c0kzJN0vaZEkP1LVzKyNmjqyiIh5AOno4mMR8UxRHUnDgPPJfpPRDcyTNCciHsgVOwFYGRF7SpoCnAMcB7wbGBER+0raEnhA0lURsbSFtpmZWT9p9t5QlfTUvPuA+yXdK+nAgmoTgK6IWBIRq4FZwKSqMpPI7jkFMBs4UpLIfvi3laThwIuA1YBviW5m1ibNdkNdDJwcEeMiYhzwEbIrpBoZDSzLDXencTXLRMQa4CmymxXOBp4FHgMeAb4ZESuqFyDpREmdkjp7enqabIqZmbWq2WTxTETc2jsQEbcBRV1RqjGu+nne9cpMAF4AdgF2B06VtN5vOiLigoioRESlo6OjIBwzM+urontDHZBe3iXpR8BVZB/mx1F8y49uYGxueAywvE6Z7tTlNApYAbwH+FVEPA88Lul2oAIsKWqQmZn1v6IT3N+qGj4j97r6KKHaPGAvSbuTXWo7hSwJ5M0hu6rqt8Bk4KaICEmPAEdIupzsKqzXAOcVLM/MzEpSdG+ow/s644hYk36TMRcYBlwcEQslTQc6I2IOcBEwU1IX2RHFlFT9fLJzIr8j66q6JCLu62ssZma2YRRRdIAAkkaRHVUcmkbdAkyPiKdKjK0llUolOjs7iwv2A+XOtDSx+gYdt29wc/sGr3a0TdL8iKgUlWvlaqhngGPT39MUXw1lZmZDRLO3+9gjIt6VGz5L0j1lBGRmZgNPs0cWf5P0+t4BSYcAfysnJDMzG2iaPbI4CbgsnbsAWEkT94YyM7OhoTBZpPtB7R0Rr5K0LUBE+NYbZmabkMJuqIhYC5ySXj/tRGFmtulp9pzFDZI+JWmspBf3/pUamZmZDRjNnrP4INkvtk+uGu9ncJuZbQKaTRbjyRLF68mSxq3AD8sKyszMBpZmk8UMsh/ifTcNT03jji0jKDMzG1iaTRZ7R8SrcsO/lnRvGQGZmdnA0+wJ7rslvaZ3QNLBwO3lhGRmZgNNs0cWBwPvT7cOB9gVWJQetRoRsV8p0ZmZ2YDQbLKYWGoUZmY2oDWVLCLi4bIDMTOzgavZcxZmZrYJKzVZSJooabGkLkmn1Zg+QtLVafqdksblpu0n6beSFkq6X9LIMmM1M7P6SksWkoaRPR71aLIf9U2VNL6q2AnAyojYEzgXOCfVHQ5cDpwUEfsAhwHPlxWrmZk1VuaRxQSgKyKWRMRqYBYwqarMJLIf9wHMBo6UJOAo4L6IuBcgIv4SES+UGKuZmTVQZrIYDSzLDXencTXLRMQa4ClgB+DlQEiaK2mBpM/UWoCkEyV1Surs6enp9waYmVmmzGShGuOqH0Fer8xwsvtQvTf9f6ekI9crGHFBRFQiotLR0bGh8ZqZWR1lJotuYGxueAywvF6ZdJ5iFLAijb8lIp6IiL8C1wEHlBirmZk1UGaymAfsJWl3SVsAU4A5VWXm8M/Hs04GboqIAOYC+0naMiWRNwIPlBirmZk10OwvuFsWEWsknUL2wT8MuDgiFkqaDnRGxBzgImCmpC6yI4opqe5KSd8mSzgBXBcRvygrVjMza0zZF/nBr1KpRGdn50ZZlnJnWobI6luH2ze4uX2DVzvaJml+RFSKyvkX3GZmVsjJwszMCjlZmJlZIScLMzMr5GRhZmaFnCzMzKyQk4WZmRVysjAzs0JOFmZmVsjJwszMCjlZmJlZIScLMzMr5GRhZmaFnCzMzKyQk4WZmRVysjAzs0KlJgtJEyUtltQl6bQa00dIujpNv1PSuKrpu0paJelTZcZpZmaNlZYsJA0DzgeOBsYDUyWNryp2ArAyIvYEzgXOqZp+LvDLsmI0M7PmlHlkMQHoioglEbEamAVMqiozCZiRXs8GjpSyBwtKegewBFhYYoxmZtaEMpPFaGBZbrg7jatZJiLWAE8BO0jaCvgscFajBUg6UVKnpM6enp5+C9zMzNZVZrJQjXHVjyCvV+Ys4NyIWNVoARFxQURUIqLS0dHRxzDNzKzI8BLn3Q2MzQ2PAZbXKdMtaTgwClgBHAxMlvR1YDtgraTnIuL7JcZrZmZ1lJks5gF7SdodeBSYArynqswc4Hjgt8Bk4KaICOANvQUknQmscqIwM2uf0pJFRKyRdAowFxgGXBwRCyVNBzojYg5wETBTUhfZEcWUsuIxM7O+U/ZFfvCrVCrR2dm5UZal3JmWIbL61uH2DW5u3+DVjrZJmh8RlaJy/gW3mZkVcrIwM7NCThZmZlbIycLMzAo5WZiZWSEnCzMzK+RkYWZmhZwszMyskJOFmZkVcrIwM7NCThZmZlbIycLMzAo5WZiZWSEnCzMzK+RkYWZmhZwszMysUKnJQtJESYsldUk6rcb0EZKuTtPvlDQujX+zpPmS7k//jygzTjMza6y0ZCFpGHA+cDQwHpgqaXxVsROAlRGxJ3AucE4a/wTw9ojYl+wZ3TPLitPMzIqVeWQxAeiKiCURsRqYBUyqKjMJmJFezwaOlKSIuDsilqfxC4GRkkaUGKuZmTVQZrIYDSzLDXencTXLRMQa4Clgh6oy7wLujoi/Vy9A0omSOiV19vT09FvgZma2rjKThWqMq34EecMykvYh65r6cK0FRMQFEVGJiEpHR0efAzUzs8bKTBbdwNjc8Bhgeb0ykoYDo4AVaXgMcC3w/oh4qMQ4zcysQJnJYh6wl6TdJW0BTAHmVJWZQ3YCG2AycFNEhKTtgF8An4uI20uM0czMmlBaskjnIE4B5gKLgGsiYqGk6ZKOScUuAnaQ1AV8Eui9vPYUYE/gC5LuSX8vKStWMzNrTBHVpxEGp0qlEp2dnRtlWcqdaRkiq28dbt/g5vYNXu1om6T5EVEpKudfcJuZWSEnCzMzK+RkYWZmhZwszMyskJOFmZkVcrIwM7NCThZmZlbIycLMzAo5WZiZWSEnCzMzK+RkYWZmhZwszMyskJOFmZkVcrIwM7NCThZmZlbIycLMzAqVmiwkTZS0WFKXpNNqTB8h6eo0/U5J43LTPpfGL5b0r2XGaWZmjZWWLCQNA84HjgbGA1Mlja8qdgKwMiL2BM4Fzkl1x5M9s3sfYCLwH2l+ZmbWBmUeWUwAuiJiSUSsBmYBk6rKTAJmpNezgSMlKY2fFRF/j4g/Al1pfmZm1gbDS5z3aGBZbrgbOLhemYhYI+kpYIc0/o6quqOrFyDpRODENLhK0uL+Cb0pOwJP5J+ZO8S4fYOb2zd4bey27dZMoTKTRa2mVj+CvF6ZZuoSERcAF7Qe2oaT1NnMQ84HK7dvcHP7Bq+B2rYyu6G6gbG54THA8nplJA0HRgErmqxrZmYbSZnJYh6wl6TdJW1BdsJ6TlWZOcDx6fVk4KaIiDR+SrpaandgL+CuEmM1M7MGSuuGSucgTgHmAsOAiyNioaTpQGdEzAEuAmZK6iI7opiS6i6UdA3wALAG+EhEvFBWrH3Ulu6vjcjtG9zcvsFrQLZN2Rd5MzOz+vwLbjMzK+RkYWZmhZwsWiTpYkmPS/pdu2Mpg6Sxkn4taZGkhZI+3u6Y+pOkkZLuknRvat9Z7Y6pv0kaJuluSf/d7lj6m6Slku6XdI+kznbH098kbSdptqTfp/fga9sdUy+fs2iRpEOBVcBlEfHKdsfT3yTtDOwcEQskbQPMB94REQ+0ObR+ke4QsFVErJK0OXAb8PGIuKOg6qAh6ZNABdg2It7W7nj6k6SlQCUinmh3LGWQNAO4NSIuTFeRbhkRT7Y7LvCRRcsi4jdkV24NSRHxWEQsSK+fARZR49fzg1VkVqXBzdPfkPnGJGkM8FbgwnbHYq2RtC1wKNlVokTE6oGSKMDJwhpIdwF+NXBneyPpX6mb5h7gceCGiBhK7TsP+Aywtt2BlCSA6yXNT7f7GUpeBvQAl6RuxAslbdXuoHo5WVhNkrYGfgp8IiKebnc8/SkiXoiI/cnuDDBB0pDoTpT0NuDxiJjf7lhKdEhEHEB2N+uPpG7hoWI4cADwg4h4NfAssN6jHdrFycLWk/ryfwpcERE/a3c8ZUmH+DeT3QZ/KDgEOCb1688CjpB0eXtD6l8RsTz9fxy4lqF1N+puoDt3pDubLHkMCE4Wto50AvgiYFFEfLvd8fQ3SR2StkuvXwS8Cfh9e6PqHxG/crEwAAACGElEQVTxuYgYExHjyO6GcFNE/Fubw+o3krZKF12QumeOAobMVYkR8SdgmaS906gjye5iMSCUedfZIUnSVcBhwI6SuoEzIuKi9kbVrw4B3gfcn/r1AT4fEde1Mab+tDMwIz1MazPgmogYcpeYDlE7Addm32cYDlwZEb9qb0j97qPAFelKqCXAB9oczz/40lkzMyvkbigzMyvkZGFmZoWcLMzMrJCThZmZFXKyMDOzQk4WZiWStKq4lNnA52RhZmaFnCzMWiDpHEkn54bPlHSGpBslLUjPWphUo95h+edLSPq+pGnp9YGSbkk3x5ubbhNvNqA4WZi1ZhZwXG74WOAS4J3pBneHA99Kt00plO7D9T1gckQcCFwMfLl/QzbbcL7dh1kLIuJuSS+RtAvQAawEHgPOTXdAXUv2/I+dgD81Mcu9gVcCN6T8MizNz2xAcbIwa91sYDLwUrIjjfeSJY4DI+L5dNfXkVV11rDukXzvdAELI2LAPD7TrBZ3Q5m1bhbZXV0nkyWOUWTPkXhe0uHAbjXqPAyMlzRC0iiyO4oCLAY6ep+1LGlzSfuU3gKzFvnIwqxFEbEw3Sr70Yh4TNIVwM8ldQL3UOOW5xGxTNI1wH3Ag8DdafxqSZOB76YkMpzsaXcLN1JzzJriu86amVkhd0OZmVkhJwszMyvkZGFmZoWcLMzMrJCThZmZFXKyMDOzQk4WZmZW6P8DdvZjnz+JN3AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# generate a uniform distribution\n",
    "# see https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.randint.html#scipy.stats.randint\n",
    "low = 1\n",
    "high = 7\n",
    "rv = stats.randint(low, high)\n",
    "\n",
    "# plot the probability mass function\n",
    "x = np.arange(low,high)\n",
    "plt.vlines(x, 0, rv.pmf(x), colors='b', linestyles='-', lw=3)\n",
    "plt.xlim(0.5,6.5)\n",
    "plt.ylim(0,0.175)\n",
    "plt.title('probability mass function (PMF) for a roll of one die')\n",
    "plt.xlabel('value')\n",
    "plt.ylabel('probability');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$F(x)$ = cumulative distribution function\n",
    "\n",
    "The cumulative distribution is the probability of occurrence of a value of $x$ or lower. For the example of one die roll, the CDF is a step function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHkdJREFUeJzt3XmYXGWd9vHvTQJESSCvJKiQsCioRF4UaHEBFUQdkCXqIBDAbVTGcRB5xf11pAcdddSggzAqg+wCBtyi4iAqoCgiHQI4ISIRwTRBCMgWBALkN388T8tJ5amu6k6drkrn/lxXX13LWX7n1Klzn/OcpRQRmJmZNdqg2wWYmVlvckCYmVmRA8LMzIocEGZmVuSAMDOzIgeEmZkVOSCGIeltkq5ci/5/JOmtnaypyXgul/TO/PgIST/u4LAXSdorP+6XdG4Hh/0xSad1angNw/6UpLsl/bmO4Q8z3q9K+peahj1L0kAdw24yvmMkfbZFN3tIulnSCkmvH6vaRktSSNo+Pz5T0qdGOZw3SFqap3uXzlY5etVpkvRySTetzfAcEB1SWnlGxH4RcdZY1hER34iI17bqrt0vR0Q8PyIuX9u6JO0labBh2J+OiHeu7bAL45oJHAfMiohndHr4lfGssQEREe+OiE/WNMpPAl9oqOFwSQN5RXVH3ijZM7/XL+kxSQ/mv99LOlnSMyv97yVpVe5/6O/7+e1TgSMlbTFMTScAJ0fE5Ij4boent5d9ATg6T/fCbhdTEhG/iIjnrs0wHBBWJGlit2tYC9sA90TEXd0upFPySn1v4LuV194PfAn4NPB0YGvgP4HZlV6/GRFTgKcBbwCeASyohgSwLK/ohv4OBIiIR4AfAW8ZprRtgEWjnKaOLmNjvMyOerrXKRHRk3/ATODbwHLgHtJWCkA/cG6lu22BACbm55cDnwJ+BawAvg9sDnwDeAC4Bti21G+l/3fmx28Drqy89x/A0jycBcDL8+v7AiuBx/I4r68OC9gYuA/YqTKs6cDDwBb5+QHAdbm7XwE7DzNvXgP8DrgfOBm4olQzIOCLwF252xuAnYCjcq0rh+ZR7v5W4MO5u0eBifm1V1fm/UXAN4EHgWuBF1TqCmD7yvMz82exSZ7WVXl8K4AtC5/lQaQv3X153u1Yee9W4AO5tvtzDZMK8+bVDeM6E9gLGGzornG65gFn5+laBPQNtywCOwKPAE/k8dxXneZKv+8ClgB/AeYDWzbMr3cDNwP3AqcAavKZvwX4SeX5Znm8bxpmOVlt/ubXJgDXA1/Iz9eYNw3dHwFc1uS9P+T5/HCuZeP8uc7P07sEeFdDPRcB55K+Q+8sDHN/YGF+fynQP0xtewGDpGX2z8A5bc7z7UufVcOwNwA+DtxG+v6cnef5xnlaA3gI+EOT/l9GWtfcn/+/rGEd80ngl3l5+zEwrfL+S0jrgPvyZ7XXMPNgF9L38EHSd+KCoWlq/GzzZ/Mt0nL8R+CYVuvhntyDkDQB+AHpw9kW2Io04e06DHhz7u/ZwFXAGaStqMXA8aMs7RrghXk45wEXSpoUEf9N2or7ZqQtsBdUe4qIR0krmDmVlw8BroiIuyTtCpwO/CMpzL4GzJe0cWMBkqaRPuSPA9NIX9I9mtT7WuAVwHOAqcChpC3rU0mB+bmobDFmc0hf0qkR8XhhmLOBCyvz4LuSNmwy/qHpfwjYj9W3VJc1TNdzgPOBY0nheTHwfUkbVTo7hBTG2wE7k8KwcVw/aRjXGt00cRBpGZtKWqmcnOsqLosRsZi0cr8qj2dq4wAlvQr4TK77mXkYjcvxAcCLgBfk7v6uSX3/F6i2J78UmAR8p83pAyAingC+B7y8zV4W59pKw3o28CfgwDwPHiV9hoOkldHBwKcl7VPpbTYpJKaSlsFGD5HCcCppOfynFsc2nkFaFrcBjmpznrfjbflvb+BZwGTSRuqjETE5d/OCPA9WI+lpwA+Bk0jf5xOBH0ravNLZ4cDbgS2AjUgbP0jaKvf7qTxdHwC+JWl6YTwbkfYoz8ndXgj8fWliJG1A2li+nrQM7wMcK6nZ8gb0bhPT7qQF7IMR8VBEPBIRIzlYfEZE/CEi7iftIv8hIn6SV3gXklJ3xCLi3Ii4JyIej4i5pK2Jdtv4zmP1gDg8vwZpi+drEXF1RDwR6bjFo6QtiUavA26MiIsi4jFSE0OzA7GPAVOA55G2TBdHxB0t6jwpIpZGxMNN3l9QGfeJpJVUqc6ROhT4YURcmof9BeAppC2xam3LIuIvpIX9hR0Y75ArI+LivAI9hydXimuzLB4BnB4R1+aV50eBl0rattLNZyPivoj4E3AZzadpKmkrccjmwN1NQryVZaQVypAtJd1X+Tuk8t6DpC3nlvKxnz2BD+f5dB1wGmljbchVEfHdiFhVWsYi4vKI+G1+/wZS4LxymNGuAo7PK+6HaW+et+MI4MSIuCUiVuThHNZmM9b+wM0RcU5eV5xP2uOvboidERG/zzXP48nP/Ujg4rwsroqIS4EB0ve+0UuADYEvRcRjEXERaSO25EXA9Ig4ISJWRsQtwH+RNqab6tWAmAncNsqFH+DOyuOHC88nMwqSjpO0WNL9ku4jfXGmtdn7z4CnSHqxpG1IC8TQ1t82wHHVLylpHmxZGM6WpF1vACLtOy4tdEdE/Iy0JXwKcKekUyVt2qLO4rBK70fEKp7cWlxbW5K29qrDXkra2hlSDcK/MsrPsYnGYU/KK4O1WRYbp2kFqYlqNNN0Lynsh9wDTBtlu/tWpOaXIcsiYmrlb17lvSmkZpJ2bAn8JSKqQXYbq0/vsMtX/n5cJmm5pPtJe2nDfceWRzpWUq2h1Txvx2rDyY8nko71jLTfof7b+dy3Ad7UsC7Yk7Q3VBrP7XkdUB1PyTY0bAgAH2s1Pb0aEEuBrZss/A8BT608X5uzVB7K/1sOT9LLSW2dhwD/Jzcp3E9q54fUJtlUXuHNI+1FHA78oPJFWgr8W8OX9Kl5y6PRHaSV1lBdqj4vjPekiNgNeD6pqemDLeptdXvf6rg3AGaQtkghLejN5mWr4S4jLcRDwx6arttb9NeO1ZaZ3Gy0xi57E8MtiyOdpk1IW/6jmaYbSJ/fkKtIx0BGdGpp/swOBH7RZi87kpol2rEMeJqkapBtzerT22qenUdq4psZEZsBX+XJ71hJ4/A6Nc9XGw5pOh5n9Y3Ndvsd6r+dGpaSjqVU1wWbRETpdOM7gK3yd6U6nmbD/WPDcKdERGnP5G96NSB+Q5r4z0raRNIkSUPt7NcBr5C0taTNSLt+oxIRy0kf2pGSJkj6B9Ixi5IppAVkOTBR0ieA6tb4ncC2+QvYzHmkppQjeLJ5CdKu3rvz1pPyNO/f8EUb8kPg+ZLemFdax9A81F6Uh7khaSU5dFB1qN5nDVNrM7tVxn0sqSns1/m964DD87zcl9WbBu4ENs+fWck8YH9J++R6j8vD/tUoamz0e9Iewf552B8nNQ+2Y7hl8U5gRsNxkqrzgLdLemE+nvRp4OqIuHUU03ApsKukSQC5+fQTwCmSXi/pqZI2lLSfpM819pzf25HUZPMMUvNgO15JaqZtKSKWkj6vz+T5tDPwDsrHGpqZQtoLeUTS7qSNqZHo1Dw/H/h/kraTNJknjzG2syd5MfAcpVOQJ0o6FJhFOpbVyrnAgZL+Ln+PJimdijyj0O1VpHXSMXk8byQ1iZb8BnhA0oclPSUPeydJLxqumJ4MiNwOfCCwPekg2CBpxUpuk/smaYtqAe3N9OG8i7RVfQ9pK7vZCukS0hfl96TduEdYfXf5wvz/HknXlgYQEVeTVtRbUvnSRcRAruNkUlPCEgoHYHO3dwNvAj6ba96BdDZEyaak8Lk313wPT55H/3VgVt7dHMn5698jfRb3ktqW35iPGQC8j/S53UcKwb8NNyJ+R/rS3ZLHuVqzVETcRGp//TJwdx7OgRGxcgS1FeWV6XtI7eG3kz6DwWF7erLfpssiqdlwEfBnSXcX+v0p8C+kkwruIG18DNvmO0wdd+bxza68diLwflLgLSctj0dTme/AoZJWkD6T+aRlYLdoOEmgJIfR64CRXMszh3QwfxmpCfX4/J1t13uAEyQ9SArAeS26X00H5/nppGNRPyed8fMI8N42a7iHdPLBcaT5/SHggPzdbdXvUtJn/DGe/Ew/SGFdnb8bbyStK+4lLZffbjLcoeX4hXl67iZ9H4Y9vqTVm6/MrFdJmkVaWe8eY/DFlfReUlPPh+oel/UmB4SZmRXV1sQk6XRJd0n6nybvS9JJkpZIukHpWgAzM+sRdR6DOJN0UVMz+5Haz3cgXdn7lRprMTOzEaotICLi56x+rnWj2cDZkfwamKrV7w9jZmZd1M0bsm3F6mcBDebX1rjSV9JRpL0MNtlkk92e97znjUmBZmZ1WLgQVq0a67EuuDsi2r3+B+huQJQufikeMY9076BTAfr6+mJgYMxuiW9m1nEa7tK/+sba7CrrproZEIOsfgVw9YpcM7P1wlidSDqaUOrmhXLzgbfks5leAtwfrW8kZ2ZmY6S2PQhJ55PuRz5N6ZfEjifdeZCI+CrpcvTXka4a/ivp1rdmZtYjaguIiJjT4v0A/rmu8ZuZ2drpyXsxmZlZ9zkgzMysyAFhZmZFDggzMytyQJiZWZEDwszMihwQZmZW5IAwM7MiB4SZmRU5IMzMrMgBYWZmRQ4IMzMrckCYmVmRA8LMzIocEGZmVuSAMDOzIgeEmZkVOSDMzKzIAWFmZkUOCDMzK3JAmJlZkQPCzMyKHBBmZlbkgDAzsyIHhJmZFTkgzMysyAFhZmZFDggzMytyQJiZWZEDwszMihwQZtZz5s6FKVNAGp9/6woHhJn1nP5+WLGi21XUb/LkblcwPAeEmfWc9SUc+vu7XcXwJna7ADOz4UR0u4L1l/cgzMysqNaAkLSvpJskLZH0kcL7W0u6TNJCSTdIel2d9ZiZWftqCwhJE4BTgP2AWcAcSbMaOvs4MC8idgEOA/6zrnrMzGxk6tyD2B1YEhG3RMRK4AJgdkM3AWyaH28GLKuxHjMzG4E6A2IrYGnl+WB+raofOFLSIHAx8N7SgCQdJWlA0sDy5cvrqNXMzBrUGRCly0Eaz0eYA5wZETOA1wHnSFqjpog4NSL6IqJv+vTpNZRqZmaN6gyIQWBm5fkM1mxCegcwDyAirgImAdNqrMnMzNpUZ0BcA+wgaTtJG5EOQs9v6OZPwD4AknYkBYTbkMzMekBtARERjwNHA5cAi0lnKy2SdIKkg3JnxwHvknQ9cD7wtghfFmNm1gtqvZI6Ii4mHXyuvvaJyuMbgT3qrMHMzEbHV1KbmVmRA8LMzIocEGZmVuSAMDOzIgeEmZkVOSDMzKzIAWFmZkUOCDMzK3JAmJlZkQPCzMyKHBBmZlbkgDAzsyIHhJmZFTkgzMysyAFhZmZFDggzMytyQJiZWZEDwszMihwQZmZW5IAwM7MiB4SZmRU5IMzMrMgBYbYOmjsXpkwBaXz+WW9wQJitg/r7YcWKbldRv8mTu13B+s0BYbYOWl/Cob+/21Ws3yZ2uwAzWzsR3a7AxivvQZiZWZEDwszMihwQZmZW5IAwM7MiB4SZmRU5IMzMrMgBYWZmRQ4IMzMrqjUgJO0r6SZJSyR9pEk3h0i6UdIiSefVWY+ZmbWvrSupJU2IiCdGMmBJE4BTgNcAg8A1kuZHxI2VbnYAPgrsERH3StpiJOMwM7P6tLsHsUTS5yXNGsGwdweWRMQtEbESuACY3dDNu4BTIuJegIi4awTDNzOzGrUbEDsDvwdOk/RrSUdJ2rRFP1sBSyvPB/NrVc8BniPpl3m4+5YGlMc3IGlg+fLlbZZsZmZro62AiIgHI+K/IuJlwIeA44E7JJ0lafsmvZXu6t54W7GJwA7AXsAcUgBNLYz/1Ijoi4i+6dOnt1OymZmtpbYCQtIESQdJ+g7wH8Bc4FnA94GLm/Q2CMysPJ8BLCt0872IeCwi/gjcRAoMMzPrsnZv930zcBnw+Yj4VeX1iyS9okk/1wA7SNoOuB04DDi8oZvvkvYczpQ0jdTkdEu7xZuZWX3aDYi3RMSV1Rck7RERv4yIY0o9RMTjko4GLgEmAKdHxCJJJwADETE/v/daSTcCTwAfjIh7Rj01ZmbWMYo2fm1E0rURsWur18ZCX19fDAwMjPVozXpK9Xeb/YNB1g5JCyKibyT9DLsHIemlwMuA6ZLeX3lrU9JegZmZjVOtmpg2Aibn7qZUXn8AOLiuoszMrPuGDYiIuAK4QtKZEXHbGNVkZmY9oFUT05ci4ljgZElrtHRGxEG1VWZmZl3VqonpnPz/C3UXYmZmvaVVE9OC/P+KsSnHzMx6Rasmpt+y5u0x/iYidu54RWZm1hNaNTEdMCZVmJlZz2nVxOQzl8zM1lPD3qxP0pX5/4OSHmj8PzYlmplZN7Tag9gz/58yXHdmZjb+tHuzPiTtCuxJOmh9ZUQsrK0qMzPrunZ/D+ITwFnA5sA00u25P15nYWZm1l3t7kHMAXaJiEcAJH0WuBb4VF2FmZlZd7X7m9S3ApMqzzcG/tDxaszMrGe0ulDuy6RjDo8CiyRdmp+/BrhyuH7Num3uXOjvhxUrul2J2bqpVRPT0C/zLAC+U3n98lqqMeug9SEcJk/udgU2nrU6zfWssSrErNPWh3Do7+92FTaetXWQWtIOwGeAWVSORUTEs2qqy6yj/LOcZiPX7kHqM4CvAI8DewNn8+StwM3MbBxqNyCeEhE/BRQRt0VEP/Cq+soyM7Nua/c6iEckbQDcLOlo4HZgi/rKMjOzbmt3D+JY4KnAMcBuwJuBt9ZVlJmZdV9bexARcQ1A3os4JiIerLUqMzPrunbvxdSXf13uBuC3kq6XtFu9pZmZWTe1ewzidOA9EfELAEl7ks5s8k+OmpmNU+0eg3hwKBwAIuJKwM1MZmbjWKt7Me2aH/5G0teA80n3YjoU327DzGxca9XENLfh+fGVx7421cxsHGt1L6a9x6oQMzPrLe2exbSZpBMlDeS/uZI2q7s4MzPrnnYPUp9OOih9SP57gHQWk5mZjVPtnub67Ij4+8rzf5V0XR0FmZlZb2h3D+LhfO0DAJL2AB6upyQzM+sF7e5BvBs4u3Lc4V58LyYzs3Gt5R5Evv/ScyPiBaQrp3eOiF0i4oY2+t1X0k2Slkj6yDDdHSwpJPWNqHozM6tNy4CIiFXA0fnxAxHxQDsDljQBOAXYj/RLdHMkzSp0N4V0l9irR1C3mZnVrN1jEJdK+oCkmZKeNvTXop/dgSURcUtErAQuAGYXuvsk8DngkfbLNjOzurV7DOIfSFdOv6fh9eF+k3orYGnl+SDw4moHknYBZkbEDyR9oNmAJB0FHAWw9dZbt1mymZmtjXb3IGaRmouuB64Dvgw8v0U/Krz2t9tz5GMbXwSOazXyiDg1Ivoiom/69OltlmxmZmuj3YA4C9gROIkUDjvm14YzCMysPJ8BLKs8nwLsBFwu6VbgJcB8H6g2M+sN7TYxDZ3FNOQySde36OcaYAdJ25F+w/ow4PChNyPifmDa0HNJlwMfiIiBNmsyM7MatbsHsVDSS4aeSHox8MvheoiIx0lnP10CLAbmRcQiSSdIOmi0BZuZ2dhQROu7dktaDDwX+FN+aWvSSn8VEBExZr8s19fXFwMD3smw1lQ5CtbGYm42rklaEBEjasJvt4lp31HUY2Zm67C2AiIibqu7EDMz6y3tHoMwM7P1jAPCzMyKHBBmZlbkgDAzsyIHhJmZFTkg1mNz58KUKel6gfH4Z2ZrxwGxHuvvhxUrul1F/SZP7nYFZusmB8R6bH0Jh/7+bldhtm5q90pqG+d8Kwoza+Q9CDMzK3JAmJlZkQPCzMyKHBBmZlbkgDAzsyIHhJmZFTkgzMysyAFhZmZFDggzMytyQJiZWZEDwszMihwQZmZW5IAwM7MiB4SZmRU5IMzMrMgBYWZmRQ4IMzMrckCYmVmRA8LMzIocEGZmVuSAMDOzIgeEmZkVOSDMzKyo1oCQtK+kmyQtkfSRwvvvl3SjpBsk/VTSNnXWY2Zm7astICRNAE4B9gNmAXMkzWrobCHQFxE7AxcBn6urHjMzG5k69yB2B5ZExC0RsRK4AJhd7SAiLouIv+anvwZm1FiPmZmNQJ0BsRWwtPJ8ML/WzDuAH5XekHSUpAFJA8uXL+9giWZm1kydAaHCa1HsUDoS6AM+X3o/Ik6NiL6I6Js+fXoHSzQzs2Ym1jjsQWBm5fkMYFljR5JeDfx/4JUR8WiN9ZiZ2QjUuQdxDbCDpO0kbQQcBsyvdiBpF+BrwEERcVeNtZiZ2QjVFhAR8ThwNHAJsBiYFxGLJJ0g6aDc2eeBycCFkq6TNL/J4MzMbIzV2cRERFwMXNzw2icqj19d5/jNzGz0fCW1mZkVOSCGMXcuTJkC0vj8MzMbjgNiGP39sGJFt6uo3+TJ3a7AzHqRA2IY60s49Pd3uwoz60W1HqQeT6J4iZ+Z2fjlPQgzMytyQJiZWZEDwszMihwQZmZW5IAwM7MiB4SZmRU5IMzMrMgBYWZmRQ4IMzMrckCYmVmRA8LMzIocEGZmVuSAMDOzIgeEmZkVOSDMzKzIAWFmZkUOCDMzK3JAmJlZkQPCzMyKHBBmZlbkgDAzsyIHhJmZFTkgzMysyAFhZmZFDggzMytyQJiZWZEDwszMihwQZmZW5IAwM7MiB4SZmRXVGhCS9pV0k6Qlkj5SeH9jSd/M718tads66zEzs/bVFhCSJgCnAPsBs4A5kmY1dPYO4N6I2B74IvDvddVjZmYjU+cexO7Akoi4JSJWAhcAsxu6mQ2clR9fBOwjScMNdMECkMbmz8xsfTaxxmFvBSytPB8EXtysm4h4XNL9wObA3dWOJB0FHJWfrgDdVEvFZdOAu8dxYEyjYX6PM+N5+sbztIGnr9O2GWkPdQZEaZUao+iGiDgVOLUTRY2UpIGI6OvGuMeCp2/dNZ6nDTx9vaDOJqZBYGbl+QxgWbNuJE0ENgP+UmNNZmbWpjoD4hpgB0nbSdoIOAyY39DNfOCt+fHBwM8iYo09CDMzG3u1NTHlYwpHA5cAE4DTI2KRpBOAgYiYD3wdOEfSEtKew2F11bMWutK0NYY8feuu8Txt4OnrOnmD3czMSnwltZmZFTkgzMysyAHRhKTTJd0l6X+6XUunSZop6TJJiyUtkvS+btfUSZImSfqNpOvz9P1rt2uqg6QJkhZK+kG3a+k0SbdK+q2k6yQNdLueTpM0VdJFkn6Xv4cv7XZNJT4G0YSkVwArgLMjYqdu19NJkp4JPDMirpU0BVgAvD4ibuxyaR2Rr8bfJCJWSNoQuBJ4X0T8usuldZSk9wN9wKYRcUC36+kkSbcCfRExLi+Uk3QW8IuIOC2f5fnUiLiv23U18h5EExHxc8bpNRkRcUdEXJsfPwgsJl3VPi5EsiI/3TD/jastIUkzgP2B07pdi42MpE2BV5DO4iQiVvZiOIADYr2X76C7C3B1dyvprNz8ch1wF3BpRIyr6QO+BHwIWNXtQmoSwI8lLci32hlPngUsB87ITYSnSdqk20WVOCDWY5ImA98Cjo2IB7pdTydFxBMR8ULSFfy7Sxo3zYSSDgDuiogF3a6lRntExK6ku0H/c27yHS8mArsCX4mIXYCHgDV+DqEXOCDWU7lt/lvANyLi292upy551/1yYN8ul9JJewAH5Xb6C4BXSTq3uyV1VkQsy//vAr5Dujv0eDEIDFb2ai8iBUbPcUCsh/JB3K8DiyPixG7X02mSpkuamh8/BXg18LvuVtU5EfHRiJgREduS7j7ws4g4sstldYykTfLJE+Sml9cC4+Zswoj4M7BU0nPzS/sAPXmCSJ13c12nSTof2AuYJmkQOD4ivt7dqjpmD+DNwG9zOz3AxyLi4i7W1EnPBM7KP1q1ATAvIsbdqaDj2NOB7+SfhpkInBcR/93dkjruvcA38hlMtwBv73I9RT7N1czMitzEZGZmRQ4IMzMrckCYmVmRA8LMzIocEGZmVuSAMOsASStad2W2bnFAmJlZkQPCrEDSv0t6T+V5v6TjJf1U0rX5twpmF/rbq/r7DJJOlvS2/Hg3SVfkG9Bdkm+7btazHBBmZRcAh1aeHwKcAbwh30Rub2Buvm1JS/neV18GDo6I3YDTgX/rbMlmneVbbZgVRMRCSVtI2hKYDtwL3AF8Md9ZdBXpNzSeDvy5jUE+F9gJuDRnyoQ8PLOe5YAwa+4i4GDgGaQ9iiNIYbFbRDyW76Y6qaGfx1l9z3zofQGLIqInf1rSrMRNTGbNXUC6W+rBpLDYjPQ7DI9J2hvYptDPbcAsSRtL2ox0p06Am4DpQ789LGlDSc+vfQrM1oL3IMyaiIhF+bbTt0fEHZK+AXxf0gBwHYVbiEfEUknzgBuAm4GF+fWVkg4GTsrBMZH0q3CLxmhyzEbMd3M1M7MiNzGZmVmRA8LMzIocEGZmVuSAMDOzIgeEmZkVOSDMzKzIAWFmZkX/C62yUSe1gS+SAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make a step function from points, function adapted from http://stackoverflow.com/a/41185239\n",
    "def get_x_y_steps(x, y, where=\"post\"):\n",
    "    zipxx = list(zip(x, x))\n",
    "    zipyy = list(zip(y, y))\n",
    "    if where == \"post\":\n",
    "        x_step = [x[0]] + [_x for tup in zipxx[1:] for _x in tup]\n",
    "        y_step = [_y for tup in zipyy[:-1] for _y in tup] + [y[-1]]\n",
    "    elif where == \"pre\":\n",
    "        x_step = [_x for tup in zipxx[:-1] for _x in tup] + [x[-1]]\n",
    "        y_step = [y[0]] + [_y for tup in zipyy[1:] for _y in tup]\n",
    "    return x_step, y_step\n",
    "\n",
    "# create step function for CDF\n",
    "xcdf = np.hstack([0,x,x[-1]+1])\n",
    "ycdf = np.hstack([0,rv.cdf(x),1])\n",
    "(x_step,y_step) = get_x_y_steps(xcdf,ycdf)\n",
    "\n",
    "# plot CDF\n",
    "plt.plot(x_step, y_step,lw=3,color='b')\n",
    "plt.xlim(0.5,6.5)\n",
    "plt.ylim(0,1)\n",
    "plt.title('cumulative distribution function (CDF) for a roll of one die')\n",
    "plt.xlabel('value')\n",
    "plt.ylabel('probability');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Mutiplication Principle__ - If two events $A$ and $B$ are unrelated, the probability of joint occurrence $P(A,B)$, i.e. both events happening, is multiplication of two probabilities,\n",
    "\n",
    "$P(A,B)  = p(B)p(A)$.\n",
    "\n",
    "This is a simplification of __Bayes' theorem__ for the special case where events are independent. In the more general case where the probability of events are not independent, \n",
    "\n",
    "$P(A,B) = p(B|A)p(A) = p(A|B)p(B)$,\n",
    "\n",
    "where $p(A|B)$ is the probability of $A$ occuring, given that $B$ has already occurred. In this case, the two events depend on each other, which is not the case for dice.\n",
    "\n",
    "__Additive principle__ - Probability of one event or another mutually exclusive event is the sum of the probabilities. Mutually exclusive means that the events cannot occur at the same time.\n",
    "\n",
    "$P(C|D) = P(C) + P(D)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sum of two dice\n",
    "\n",
    "From the rules above, we can calculate the probability mass function and cumulative distribution function for the sum of two dice. From the multiplcation rule, the probability of rolling a particular sequence of dice is (1/6) $\\times$ (1/6) = 1/36. This is the probability of a joint occurrence of two independent events. For example, the combinations (1 and 1), (1 and 2), (2 and 1) and (1 and 3) are all joint occurrences of two independent events.\n",
    "\n",
    "The probability of obtaining a certain sum, which may be arrived at through several possible sequences of dice rolls, can be calculated from the addition rule. There are three ways of obtaining a sum of four: (1 and 3) or (2 and 2) or (3 and 1). These combinations of rolls are mutually exlusive. Therefore the probability of rolling a four is (1/36) + (1/36) + (1/36) = 1/12.\n",
    "\n",
    "__Exercises__:\n",
    "\n",
    "* Graph the probability mass function and cumulative distribution function for the sum of two dice.\n",
    "\n",
    "* Roll actual dice and calculate the sum of two dice from multiple trials. Plot the observed probability mass function from this limited set of trials. Compare the observed distribution with the theoretical distribution derived from probability theory above and discuss."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distributions\n",
    "\n",
    "### Normal Distribution\n",
    "\n",
    "<img src=\"images/norm_dist_rule.png\" width=\"450\">\n",
    "[Wikipedia](https://en.wikipedia.org/wiki/Normal_distribution#/media/File:Empirical_Rule.PNG)\n",
    "\n",
    "\n",
    "* __Mean__ - 1st moment\n",
    "* __Variance__ - 2nd moment\n",
    "* __Skew__ - 3rd moment, Describes the spread (Normal Distribution, skew = 0)\n",
    "* __Kurtosis__ - 4th moment, Describes how pointy the distribution is. Low value would be flat (compared to a Normal distribution (k=3))\n",
    "\n",
    "### Non-normal distribution \n",
    "__Example: Life expectancy__<br>\n",
    "\n",
    "Negative skew towards later life (spike at birth). The cumulative probability function is useful because we are not necessarily interested in the probability of dying at a very specific age, like exactly 67 years. We are more interested in the probability of living until a certain age or less. The cumulative probability distribution tapers off to one for very high ages. The probability of living to 120 years or less is equal to one.\n",
    "\n",
    "<div style=\"float: left;\" width =\"1000\">\n",
    "    <img src=\"images/life_expectancy_1.png\" style=\"float: left;\" width=\"50%\"><img src=\"images/life_expectancy_2.png\" width=\"50%\" style=\"float: right;\">\n",
    "</div>\n",
    "source: http://www.cureffi.org/2013/06/19/how-censoring-by-life-expectancy-affects-age-of-onset-distributions/\n",
    "\n",
    "\n",
    "### Poisson Distribution\n",
    "  * Often used as a model for count data\n",
    "  * Assumptions:\n",
    "    * One event does not affect the probability of the next (they are independent)\n",
    "    * Events cannot occur in the same time and place in the interval\n",
    "    * Events occur at a constant rate<br>\n",
    "\n",
    "    $k$ = # of events counted per interval (ie floods per century)<br>\n",
    "    $\\lambda$ = expected values, true mean\n",
    "    \n",
    "     $P(k$ events per interval$)$ = $\\frac{\\lambda e^{-\\lambda}}{k!} $\n",
    "\n",
    "    Count uncertainty  = $\\sqrt{k}$ - More counts => less uncertainty\n",
    "\n",
    "__Exercise__: if on average 10 flood occur per century. What is probability of counting 15 floods occuring in a centry?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Weibull Distribution \n",
    "\n",
    "A theoretical distribution that is often fit to wind speed data. Wind speed has only positive values, so the probability of a negative wind speed is zero. The distribution of wind speed is generally skewed, since very high values do occur, but are rare. This theoretical distribution is useful for engineering wind power projects, since it can be defined by a few parameters. The parameters obtained from a fit to data at a given location can be used as inputs to models for designing equipment with maximum efficiency.\n",
    "\n",
    "<img src=\"images/weibull.png\" width=\"500\" align='center'>\n",
    "Source: http://www.wind-powemr-program.com/wind_statistics.htm\n",
    "\n",
    "\n",
    "### Log-Normal Distribution \n",
    "\n",
    "Many biological data, such as growth rates, follow a log-normal distribution. The distributions are skewed, with a peak close to zero but a long tail of rare high values. These data can be transformed by taking the logarithm, giving a distribution that is more symmetric and closer to a normal distribution.\n",
    "\n",
    "<img src=\"images/log_normal.png\" width=\"600\" align='center'>\n",
    "source: Landry, M.R. and B.M. Hickey (eds.) (1989) Coastal Oceanography of Washington and\n",
    "Oregon, 607 pp., Elsevier Science, Amsterdam, The Netherlands."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Back to Normal Distributions\n",
    "\n",
    "Random instrument error is often normally distributed. The mean of errors tends to be zero if error is randomly distributed and _accurate_, i.e. not biased.\n",
    "\n",
    "<img src='images/instr_error_dist.png' width='500'>\n",
    "source: T. Connolly\n",
    "\n",
    "In this example, error is defined as: $T_{sat} - T_{mooring}$ units of [$^o$C]. There is a fairly large spread in the data. We therefore cannot be extremely confident in any particular satellite image. However, we can be confident that the mean of a large number of errors is close to zero.\n",
    "\n",
    "Normal distributions can be \"standardized\", where the\n",
    "* mean = 0\n",
    "* standard deviation = 1\n",
    "* Unitless\n",
    "\n",
    "The __z-score__ transforms the data so that each data point is unitless and described in terms of standard deviations away from the mean.\n",
    "\n",
    "## $ z_i = \\frac{x_i - \\bar{x}}{ s } $ ##\n",
    "\n",
    "The z-score can be helpful for identifying outliers (large Z-scores). Outliers are often classified as having a z-score of $\\pm$ 3. However, the basic rules of scientific integrity require that you have a good reason for excluding outliers from your analysis, and report that they have been excluded. Never delete outliers from the original data set.\n",
    "\n",
    "### Describing Uncertainty of Mean Values\n",
    "- Random instruemtn error is often normally distributed\n",
    "- We have better confidence in the mean than any one value\n",
    "- We need to be able to quantify our uncertainty in our estimate of the mean.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Student's $t$ distribution:\n",
    "\n",
    "If samples come from a parent population that is normally distributed with true mean $\\mu$, the $t$ statistic of a group of random samples can be defined as \n",
    "\n",
    "## $ t = \\frac{\\bar{x} - \\mu} { \\frac{s}{\\sqrt{N}}}$ ##\n",
    "\n",
    "where $N$ is the number of samples. The probability distribution of this $t$ statistic (i.e. from multiple groups of random samples) forms a $t$ distribution. The denominator in this equation is equal to the standard error. The shape of the $t$ distribution depends on the degrees of freedom $\\nu$. The t distribution generally has less kurtosis than the normal distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x112e556d8>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGHCAYAAACNjTnqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xl4VNX5wPHvmw0ChCCENQTcAEFAIICABSxUEPeqVVBq\nXepuVbSbttXu7a8qttbWpbUuValWrVA3rLLLvqNsLqwJ+xICSSDL+f1x7iT33plJMslkMknez/Pk\ngXvm3HvP7O/c855zxBiDUkoppVR9SKjvBiillFKq6dJARCmllFL1RgMRpZRSStUbDUSUUkopVW80\nEFFKKaVUvdFARCmllFL1RgMRpZRSStUbDUSUUkopVW80EFFKKaVUvdFApIkQkTkiMsu13V1EykTk\nuhic+3rnXN1cZVtFZEZdn9s512jn/KNicb4Q5/+BiHwpIiUisrI+2hBv6vs5qQuhXueqaiJyvois\nEpFCESkVkdaV1B0sIp+IyFGnbv9YtjWcxvh6jiUNROKciAwXkYcre3NWU6i5/COe319EbheR79Tg\n3P5zRX1tgSraVi9rGYjIOOD/gPnA9cCDdXy+y0TkAxHJEZEiEdkhIv8WkTPr8rw1VOVzIiJ9ROR1\nJ5A7JiL7RGSuiFwUpv4Zzv3PF5EDIvKSiGSEqXuTiKx3vgA3i8hdUbg/umZGBESkLfAaUADcAXwb\nOBambhLwBnAScK9Td1tsWlot+tzXUFJ9N0BVaQTwEPA8cCRaBzXGbBORVKA4wl3vAPYBL0awz0vA\nNGPMiQjPFamQbTPGzBWR1BicP5SvA6XATcaY0hicrx9wEPgjsB/oBNwILBWRYcaYdTFoQzR1B1oB\nLwC5QAvgCmCGiNxijPl7oKKIZGIDvkPAj4E04AdAXxEZaowpcdW9FXgK+DfwGDASeMJ5nTwSizum\nABiCfX5/aoyZXUXd04Bu2PfS83XeMhUzGojEP6mrA9f1F7OItDDGFBi7smJ9BAHl6ikIAegIFEYz\nCBGR5saYolC3GWN+FaL+c8BO4HZssNZgGGPeB953l4nIk8BK4D7g766bfgKkAgOMMTlO3WXA/7BX\no/7ulDUHfg381xhztbPvcyKSCPxMRJ41xuTV2Z2KMqfdCcaYSH9UxIOOzr/VebyrXTfw2VPjVqmY\n0q6ZOCYiDwN/cDa3On2QpVX1QYvILSLyhYgUiMhiEflaiDpBOSIi0lFEnncu5xeJSK6IvB04n4hs\nAc4EznX2LQvknbj6x0eJyF9FZA+ww3dbULtF5DxX//BnIvJN3+0/F5GyEPt5jllF20L234rIt0Rk\nufM47RORf4pIF1+dF5zL/F2cxyJfRPaKyCMiUmmQ6LT7O0BL13N3nXNbooj8zHmeikRki4j8RkRS\nfMfYKiIzRGSciCwTkULglsrOG8I+7KXvNlVVFJFLROQdV9fOFyLyUxFJ8NWbIyJrRaS3iMx2uk12\nisgPQhwz03nsjorIHhGZCjSjhkG2E9juCHF/LgfeCQQhTt2Pgc3AVa56XwfaAn/17f8X7K/zC6tq\ng9guo1nOa2eHiPyEMJ+nIjJBROY59/+I8/j2CVHvW857oNB5bC9zXn9bXHUC79v7ROQeEfkCKAJ6\nO7eniMgvRORz5/nbLiL/539dOXUnu17/B0Rkmoh09dU5XUTeFJFdTrt2OPXSqvEYVfr+EpHZ2Ctd\nAMud+/WPMMd6HpiD7f54w/f+DrxHTxWR90TkCPCya9+zxXbXHXZep3NEZESIc3QRkX+IyG7nsftU\nRG4IUa/ar+cIP2OynNdGvvNeusO5vZ+IfOycb6uITKr8kW949IpIfHsT6AlMBO4BDjjl+8LtICI3\nAU8DC4DHgVOBGdjL9durON9b2A+0J7B9rx2A87CXQ7c7bXgSyMf+ohRgj7NvoH/0r8Be4BdAS9dt\nofpPewL/ctr7AnAD8G8RGe98gVS2r7+8srbhq4uIXA/8A1iCvYzfEdvvPEJEBhpjjrj2SwBmAouB\n+4FvYH+NfwE8E6JtAZOBW7GXn29y2rTQue054DrgdeBR4GzgAeAMbNeDu91nAK8653oW2FTJOQP3\nLx1IxnbNTMF2U3xU1X7YKwf52O6Ko8AY4JfO/j/ytast9mrFW9jn8Urg9yKy1hgz02lHc2AW0BX4\nE7AL27c/hgj61EWkBfZqRzpwKTABmOa6vQv29bo8xO5LnfoBA51/V/jqrQDKnNtfraQtHbFfignA\nb7FB3i3YgMBf99vY1/YHwA+xXUu3A/Od19l2p96F2MdwDfb1eBL2NZJD6MfpRuyX3zPAceCgiAjw\nX2x37jPARmxX3RSgBzZQC7TrJ9jn9V/A34D2wN3A3MDrX0SSgQ+xr6MngN1AJnARNgjMr+Qxup6q\n31+/xr6WbwZ+CmwFvgxzyKexV/V+gn0dLcP72ZOEfY/Ox75HC5x2jAHew74ufo59fm8AZonI14wx\ny516HZy2ljr3dT/2NfOciKQZY55w6lX79VyDz5j3gbnY7sRrgT+LyDHgN9jA6k3gNuBFEVlojImn\n/JjaMcboXxz/Yd9UpUC3atRNwn5YLAeSXOU3Yd+As1xl3Z2y65ztdGf7virOsc59HFf5d5z95wAS\n4jbPfQC2OGWXusrSsB+8y11lDwOlYc7nP2a4to126o7yPU6rgRRXvQuc+/Cwq+x5Z98HfcdcASyt\nxnPyPHDEV9bfOc/TvvI/OOcaHeJx+kaEr5sNzjnKsJeyf1HN/ZqFKHsK+6WT7Cqb7bTrGldZMjaP\n43VX2T1OvctdZc2xVynKn5NqtOsp1/0pwSY4prtuz3ZuuzbEvv/nnCvZ2f4zcCLMefYAr1TRlsed\n42W7ytphc1PKX5PYQPwg8JRv//ZO3addZWuxwX+qq2ykc5++cpUF3reHgLa+407G5nwN95Xf4rRr\nmLPdzan3I1+9Ptgu1B8722c55/pmhK+9SN5fgffxoGocd7Sz/+W+8sB79Nch9tkEvOt/jWMDng9c\nZX/HBjptfHVfdZ7DZpG8niN8DALt/6GrLB2btFsCXOkq7+ns/1Akz0m8/2nXTOMyGPur8GnjSszD\nJm9W1a9aiP0QOldEqryEH4YB/macd0w15BpjppfvbEw+NrF1oPMLpa4EHqe/GlfuiDHmPeyvyFCX\n5v1XPuZjrzbVxAXYx+pxX/lj2Ksm/vNvMcZU52qG2/XAeOyv7w1AqthRB5UyxhwP/F9EWolIO+zV\ntRbYKzNuR40xr7r2LcZefXA/LhOAXcaYt1z1irBXdiLxOPZK1HXYX7iJ2C+UgFTn3+MEK/LVSSV8\nzlKRq144E4DFxpjyKyrGmAPAK75647BfKP8SkXaBP+xzvwTbRYSIdAb6Ai8aYwpdx5yPDa5DecMY\nc9BXdiX2ud7sO99s7Ovq6069K5ztf/vq7QU+d9ULfGacLzaxvbpq8v6KhqfdGyIyAHslaJrvfqYB\nHwPurtrLsVeTEn11P8Re/Rnk1Kvu67kmj8Fzrnp52CDqmDHmDVf5ZuAwNf/siUvaNdO4dMd+yH3h\nLjTGlIjIV5XtaIw5ISI/wnYT7BGRxcA7wEvGmD2V7euzNYK6X4Qo2+z8ezL2g7EuBB6nzSFu2wic\n4ysrcr5o3A5hL5/X9PxlBD9Pe0TksHO72xYiZIxZEvi/iLyG/YIC2z0Qltjchd9gv4zcQ8YN9kvV\nbWeIQxzCdgcEdCf081xl95Kb8wEceL5eFpGZ2C+Os52ywBd4M/++2F+s7jqFQFDOhKtuYZjbArpj\nu+n8/PfpdOwXfqjRIIaKL/rA8x2qW+ILKrqS3LaGKOuBDRZDdd0a7BdjoF0JhH5eyhPLjTFbReQx\nbDfkZBGZj+3mfdlUdCuEEun7KxpKjDH+12MP59+XwuxT5nRhpmCDjVuwXal+7seuuq/naHzG5BH6\nPZZHzT974pIGIqqcMeZPYicZuwz7a/qXwAMi8nVjzJpqHqaqD/GImxWmPDHK56lMXQ27re6Vo1o9\npsaYw05i37VUEog4H8rzsL+4fgp8hb1CkA38nuBkzHCPS52N9HJ5A3haRHoYYz7H9tUDdA5RtzNw\n0FSMKtmF/eWbYYzZH6jk5ES0w3YvRUMC9jmejDdfKaAkRFl1hXpNJGCvoEwh9HOww1WvDDjf+dfv\naOA/xpgfiMgL2Lyccdj8iR+LHQoerccpGkJdCQu8Xu/H5t6EchQIzDPzMuGnJVhb86ZVS7j3Un2+\nx2JGA5H4V+2EPmwfs2B/CcwJFDqX5E/B9ldWfjJjtmAvgz8uIqdh38D3Yy+JR9qeqpweoqyX8+9W\n599DACLS2vcr7OQQ+1a3bYHHqReux8l1/rpOAtuG/ZDsgeuXlNMd1aaOzh9I9KzMudhfWpcaYz5x\nteu0Wpx3G3Y0k5+/mydSga6CdABjTK6I7MNeEvcbive1vxr7/A/GJpEGDME+L1W9T7ZR8WvbzX+f\nvnTOs88YMytEfffxIPT7IVRZOF8C/U3V83EE2rXVGBPq172HMeYz4DPgtyIyDJtwfRt2fqNQ6vv9\nFRC4wpRf2ePvvG7ygcQqnieo/us5Xh6DBkFzROJfYJbB6uRtLMdelr3Nlw9wQ1X7i0iqiPgva2/B\nvkHd5ceq2Zbq6CKu4bpiZ4/9NrDKGBPolgl8aI5y1WtJRWDkVt22Lcd2+9zm/AoOHHcCdtTQOxHe\nj0i9h71P9/rK78cGU+/W9MAi0j5E2cnAWOxIg8qUOu0q/1wQO+yzNnOPvId9nstHAjkjYG6uzs5h\n7k8SNsmxEFjvuulN4CKxE5sF6o7FJvi97qo3C5uAeLvv0LdjX0NVPf7vAcNEpDzocdp5ja/eTOwk\nhA+Gys8RZ8ZXY8wu4FPgOuexCdw+Gm83V1VeB7qKSNBjKyLNXcd+CydhMtRBxM52ioikiZ2jxO0z\nZ99QXWAB9f3+CliB/fz4vvOZ4eF6/Muwr50rJMQMxOKdmbe6r+d4eQwaBL0iEv9WYL8cfisi/8Jm\nu89wJ7UFOLkgP8Umbc12cgNOwQYi4YbFBfQEPhaR17Ef7iXYBK4OuIZJOu25zRn+9wWw1/ULLNLL\nhZuBv4vIEOyl65uc833HVedD7NDhf4jII1QMv9sLZPmOV622OY/Tj7BD6+aJyDTsMNe7sd0Rf4zw\nfkTEGLNWRF4EbhGRk7BD9s7GBldvGWPm1uLw60TkY+yv+kPY5/VG7Hv9x1Xsu9DZ5yURecIpm0zt\nroL9DbgL+KfzxR0Y7hhyGu8QnnEC1HnYEVWdsF1MvbAjvNyTVv0Wm7A5R0T+hE1K/D72qt4LgUrG\nmCIR+RnwpPN6n4kNdK/Bjo46XEWb/uDch5nOeQqwX0RbsSOiAufJF5HbsTkKK5337z7sqJULsUnA\ndzvVHwTeBhaKnTOjLXAntqulVbUeKfgndr6Up0Tk68An2C7M3sC3sF0rK40xXzmfE78VkVOc8+Zj\nEyAvwyZmT8UOSX1SRP6Nfa8mYV+jJdgv7pBq8P6qk24GY4wRke9ig4fPnMc1BzsE+evYXItLneo/\nxl4RXCIif8N+BrbFdkuOoaL7plqv5/r+jGlw6nvYjv5V/Yf9kNqODUKqHMqLTbj6AvsBuQSbGDUL\n+NhVp7tzrMDw3bbY/t/PsL/iDmK/mPxD5TpgE9YOO/vPcsrDDsMj9FDbr4Dp2JEQq522fkaIoYLA\nAKcthdirNHeHOWa4tnmG77rqX4n95VJAxdTwnX11ngfyQrTpYWyCXFXPXbj9E7B5GIEJqbYCv8I1\nRNb9OEXwWnnIec73Y/vNd2D7vs+s5v7DsF9gR519f+s8R57HD5uAuSbM/f3SV9YV+A/2y24PdnTQ\neaGekxDHuwobKOQ692e/s31hmPq9sfMx5GPn3XkRaB+m7k3YL5xC7Bft9yJ4nM/EvqeOYd+bD2AD\n5KD3JzbIeQ/7njrmnOs5YKCv3rec90AhNgC5BDsF/Wch3rdTwrQrERt8rXVe1/uxI5l+ArTy1b0M\nGwQfcf4+w86Ncbpz+8nYL97NTrv3YeeiObeaj1F13l+RDt/1DJ2t7D3mur2/8zjuddryFfbH1bm+\nehnYz8Ct2PdkDvaH0I01fT1X8zEI9xkR7j0W0WdCQ/gT544ppZSKMyKyCntlb3x9t0WpuhI3OSIi\ncqfYaa4LxU5LPqSa+50jIsUSYnl1sdPrbnCOucbpn1NKqbgiIkn+fAwRORc7qVhVyadKNWhxEYiI\nyNXYy1sPY8fMr8H2v4Zcvtu1Xzr2UlfQZE9i1xJ4FXtpcQC2G+BtCbHGg1JK1bNMYKOIPCwiN4td\nv+RdbJdUZcsIKNXgxUXXjDN51hJjzD3OtmD7p58wxvyhkv2mYfsvy7BDDge5bvsX0MIYc4mrbBF2\nREaDWoFUKdW4OQm5z2DzudpjczI+Ah4wdki9Uo1WvV8RcYY2ZWOn3AXKV9f8CBheyX43YEeE/CJM\nleEEXymZWdkxlVKqPhhjjhhjJhljuhljUo0xGcaYiRqEqKYgHobvZmAzvf0zD+6hYnIrDxHpgc3m\n/5oxpkxCr8beKcwxO9WqtUoppZSKmngIRCIiIgnYxaUeNsYE5saIyjh0sYscjadi+JZSSimlqqc5\ndtj3TBO8dk5Y8RCI7MeOv+7oK++IXUbZLw07NfMAEfmLU5aATS05AYwzxsxx9q3uMQPGE7yCplJK\nKaWq71rsYJFqqfdAxBhTLCIrsFNQz4DyZNWx2Mll/I5gl8x2uxM7U94VVKxRsijEMc5zysPZCvDy\nyy/Tu3fvSO5GkzZlyhQef9y/or2qij5uoc3bNo8f/u+HFJcWB9/4AXapNpekxCT+8I0/MPrk0TFp\nX0Okr7XI6WMWuQ0bNjB58mSIbBX2+g9EHFOBF5yAZCl29cgWONMyi8jvgC7GmO84iazu9SUQkb3Y\nZZQ3uIr/hJ3q+T7sMLhJ2KTYyta4KALo3bs3gwYNqqSacktPT9fHqwb0cQu2YPsCfvjuDynu6A1C\nkhOSOafbOaxbuo4jXY9QXFZxewklPLD+Af6X/T9Gdh8Z6yY3CPpai5w+ZrUSUWpDvY+aATDGvI6d\nlviXwCrslLzjjTH7nCqdCF5XpKpjLsKuHXELdgrxy7FDfNdXuqNSql7syNvB5a9d7gkyAMafNp4v\n7/6S2d+ZzYisEXx595eMP8070ejx0uNc/vrl5BzJiWWTlVJREBeBCIAx5q/GmJOdoWvDjTHLXbfd\nYIwZU8m+v3DPIeIqf9MYc4ZzzP7GmJl11X6lVM2VmTK+/Z9vs69gn6f8hgE38O4175KVXvE7JCs9\ni3eueYcbBtzgqbu/YD+T/zOZ0rLSmLRZKRUdcROIKKWarqeWPcXcbd5FhyecPoG/Xfw3EhP8K9FD\nUkISz178LBf1vMhTPmfrHJ5a/lSdtlUpFV0aiKhamzRpUn03oUHSx83KOZLDjz76kaese3p3Xr3i\n1aAgxP2YJSUk8dJlL9EtvZunzoMfP0hufm7dNbgB0tda5PQxi524mOI9XojIIGDFihUrNElJqRi5\ncfqNPL/6eU/ZR9/+iLGnjq3W/nO2zuHrL37dUza5/2T++c1/Rq2NKj5t376d/fv313czmoyMjAy6\ndesW9vaVK1eSnZ0NkG2MCVqINpx4GTWjlGqC1uxewwurX/CU3TTwpmoHIQDnnnwu3znrO7y45sXy\nslfWvsIPR/yQfh37RaupKs5s376d3r17U1BQUN9NaTJatGjBhg0bKg1GakIDEaVUvXlw1oMYKq7K\ntkppxa/H/Dri4zxy3iP8Z+N/OHL8CAAGw89m/4y3J74dtbaq+LJ//34KCgp03qcYCcwRsn//fg1E\nlFKNw6pdq3jv8/c8ZT8650d0ahX5clDtW7bn/uH38/Cch8vLpm+azrKcZQzJHFLrtqr4pfM+NXya\nrKqUqhe//+T3nu2MFhlMGTalxsebMmwKGS0yPGWPLnq0xsdTSsWGBiJKqZj7/MDnvLH+DU/ZvWff\nS8uUljU+ZlqzNL4//PuesjfWv8GWQ1tqfEylVN3TQEQpFXN/WvInykxZ+XZaShp3Dr2z1se9dfCt\ntEppVb5dZsr44+I/1vq4Sqm6o4GIUiqmjp44yktrXvKU3Tb4Nto0b1PrY7dp3oabB3mXk3pu1XPl\nSaxKqfijgYhSKqZeXfcq+Sfyy7cTJIG7ht4VtePfc/Y9JErFRGjHio8xbd20qB1fKRVdGogopWLG\nGBM0BfsFPS4Imh21Nrq36c7FvS72lD278tmoHV+pWFi/fj1XXXUVp512Gi1btqR9+/aMHj2ad955\np76bFnUaiCilYmZJzhJW717tKbt98O1RP88tg27xbK/ctZIVuSuifh6l6sq2bds4evQo119/PU88\n8QQPPfQQIsIll1zC3//+9/puXlTpPCJKqZjxz6J6cpuTGX/a+KifZ9xp4+iW3o3tedvLy55d8SzP\ndHkm6udSqi5MmDCBCRMmeMruuusuBg0axNSpU/nud79bTy2LPr0iopSKieMlx3n9s9c9Zd8d+N2Q\nq+vWVmJCIjcNvMlT9tpnr1FUUhT1cykVKyJCVlYWhw8fru+mRJVeEVFKxcR7n7/HoaJDnrJr+19b\nZ+e7YcAN/HzOz8unkM87nsd7n7/H5b0vr7NzqvhTVgYHDsTufO3aQUIUf+IXFBRQWFhIXl4e06dP\n5/333290KwNrIKKUiomX173s2R7VfRQntzm5zs6XlZ7F6JNHM2frnPKyV9e9qoFIE3PgAHToELvz\n7d0L7dtH73j3338/zzxjuxQTEhK44oor+POf/xy9E8QBDUSUUnXuUOEh3tnszfaf3G9ytfYtLYWl\nS+HTT0EEzjwThgyBpGp8el3T9xpPIPLO5nfIK8ojvXl6JM1Xqt5MmTKFb33rW+Tm5vL6669TWlrK\n8ePH67tZUaU5IkqpOvfG+jc4UXqifDslMYUr+1xZ6T7GwIsvwumnw4gRcMstcPPN9v+nnWZvM6bS\nQ3BFnytITkgu3z5eepz/bPxPre6LUrHUs2dPxowZw+TJk5kxYwb5+flccskl9d2sqNJARClV597c\n8KZn+8IeF3JS6klh6xcVwcSJcP31sHVr8O3bt9vbJk6EwsLw522b2pYLelzgKXt13avVb7hScebK\nK69k2bJlfP755/XdlKjRrhmlVJ06VHiIj7d87Cmb2Hdi2PpFRXDppfDhh1Uf+/XXIS8Ppk+HZs1C\n17mm3zVM3zS9fHvWllkcLDxI29S21Wq/atjatbN5G7E8X10qdCLvvLy8uj1RDGkgopSqU+9sfoeS\nspLy7WaJzZhw+oSw9e+6K3QQ0quXzRHZuNFbPnMm3H47/OMfoY93YY8LaZ7UvHzobqkp5Z3N73Dd\nWddFfF9Uw5OQEN3k0VjZt28f7X0NLykp4cUXXyQ1NZU+ffrUU8uiTwMRpVSdemvjW57tcaeNI61Z\nWsi6zz8Pzz3nLWvd2uaDXHqp3Z4+3XbLuH8QPv88jBwJN9wQfMyWKS0Zd9o4ZmyaUV72n43/0UBE\nxbVbb72VI0eOMGrUKDIzM9m9ezevvPIKmzZtYurUqbRo0aK+mxg1miOilKozx04c44MvPvCUhRs+\nm5sL997rLWvZEj76CC67zF4NEbH///BDe5vbvffaY4TyzTO+6dme+cVMCooLIrovSsXSxIkTSUxM\n5Omnn+aOO+7g8ccfJysrixkzZnDPPffUd/OiSgMRpVSd+eCLDzyzmSZKIhf3vDhk3XvvhSNHvGXP\nPWeH6voNHQr+5TaOHIG77w7djot7XuxZkbewpJCZX8ys1n1Qqj5cddVVzJw5k9zcXI4fP87+/fuZ\nOXMmF154YX03Leo0EFFK1Rl/t8y5J59LuxbB2Xzz5sG//+0tu+46uPrq8MeeONHWcXvzTZg/P7hu\nuxbtGNV9lKdMh/EqFR80EFFK1YmSshLe//x9T1mobhlj4Cc/8Za1awePPVb1OR57LHiUwoMPhp5f\nxN8989/N//Uk0Sql6ocGIkqpOrE0Z2nQ2jKhumU+/BAWLPCW/eIXkJFR9TkyMuDnP/eWLVgAH3wQ\nXPeyMy7zbB8uOszinYurPolSqk7FTSAiIneKyBYRKRSRxSISome4vO45IrJARPaLSIGIbBCRe311\nviMiZSJS6vxbJiKanaZUjLz3+Xue7X4d+pGVnhVU75e/9G537w6RrHB+yy1w8snest//PrheVnoW\n/Tr085T5r9gopWIvLgIREbkaeAx4GBgIrAFmiki430THgD8DI4EzgF8BvxYR/8dXHtDJ9dc9+q1X\nSoXiD0RCzR2ydCksXOgte+ih8JOThZKSAg8/7C2bNw+WLQuu659l9f0vNBBRqr7FRSACTAGeMca8\nZIzZCNwGFAA3hqpsjFltjHnNGLPBGLPdGPMqMBMbmPiqmn3GmL3O3746vRdKKQBy83NZtXuVp8wf\nBAD86U/e7a5d4dvfjvx811wDXbp4y0LlmPiDoVW7V7Erf1fkJ1RKRU29ByIikgxkA+VzQBtjDPAR\nMLyaxxjo1J3ju6mViGwVke0i8raINJ6p6JSKY/65Q1o3a82IrBGespwcO0W72113QXIyEUtJge99\nz1v2xhv2HG4jskbQulnrStuqlIqteg9EgAwgEdjjK9+D7U4JS0R2iEgRsBT4izHmedfNm7BXVC4B\nrsXe14Ui0iX4SEqpaPJ3y4w7bRzJid4I429/gxLXoJXUVLu6bk3deqt3krPSUnjhBW+d5MRkzjv1\nPG9bv/C2VSkVWw19ivevAa2AYcD/icgXxpjXAIwxi4HylHgRWQRsAG7F5qKENWXKFNLT0z1lkyZN\nYtKkSdFtvVKNUHFpMR9+6V0s5oLTvd0yZWXBQcJ110HbWqxDd9JJMGmSd6Kz556DBx6w642Ut6XH\nBZ7VgP/35f8oLi0OCpSUUuFNmzaNadOmecpquhBfPAQi+4FSoKOvvCOwu7IdjTHbnP9+JiKdgJ8D\nr4WpWyIiq4DTq2rQ448/zqBBg6qqppQKYdHOReSfyPeUnX/6+Z7tOXNg2zZPUa2uhriP4Q5EtmyB\nWbPgG98I35a843ks3rmYkd39KWZKqXBC/ThfuXIl2dnZER+r3rtmjDHFwApgbKBMRMTZXhhuvxAS\ngbC59iIlM04ZAAAgAElEQVSSAPQDNDNNqTr00VcfebYHdBpA57TOnrLnn/ds0q8fRCP2HzLEHsvN\nPxV8l7QuDOg0wFPmb7NSKnbqPRBxTAVuFpHrROQM4GmgBfACgIj8TkReDFQWkTtE5CIROd35uwm4\nH/inq87PROQ8ETnFSWZ9BegG+D6WlFLR9PGWjz3b3zjlG57tvDw7FbvbDTfYBe1qSyR4DpLp0yHf\ne4EmqE3+NiulYicuAhFjzOvA94FfAquA/sB413DbToB7JqQE4HdO3WXA7cAPjDHu3I+TgGeB9cC7\n2FyS4c7wYKVUHThy/AhLdi7xlI09daxn+623oLCwYjspCSZPjl4brr3WHjOgqMgGI5W1aUnOEvKP\n+6IVperR+vXrueqqqzjttNNo2bIl7du3Z/To0bzzzjtBdTdu3Mj5559PWloa7dq147rrrmP//v0h\nj/vcc8/Rp08fUlNT6dmzJ08++WRd35UqxUUgAmCM+asx5mRjTKoxZrgxZrnrthuMMWNc208aY/oZ\nY9KMMScZYwYbY571He8+Y8wpzvG6GGMuNsasjeV9UqqpmbdtHqWmtHw7OSGZkd28uRf+xe0uuADa\nt49eG9q1g3HjvGX/+pd3e2S3kSQnVCSnlpSVMG/bvOg1Qqla2rZtG0ePHuX666/niSee4KGHHkJE\nuOSSS/i7q78xJyeHkSNH8tVXX/H73/+eH/zgB7z77ruMGzeOkhLvWkrPPPMMN998M/369ePJJ59k\nxIgR3H333TzyyCOxvnse8ZCsqpRqJD7+ytvFMTxrOC1TKsbUHjoE//ufd5+JE6PfjokT4T3XqNyZ\nM+HAgYoF8lqmtGR41nBP8PHxlo+5sGfjW2JdNUwTJkxgwgTvBHx33XUXgwYNYurUqXzX6YP8zW9+\nQ2FhIatXryYzMxOAIUOGcN555/HCCy+U1ysqKuKnP/0pF198Ma+9Zsd03HTTTZSWlvKrX/2KW265\nJWi0aKzEzRURpVTD58+1GHuKtwvk7be9c4c0awYXXRT9dlx2GTRvXrFdUhKcl+JvmyasqngnImRl\nZXH48OHysrfeeouLLrqoPAgBGDt2LD179uR114yBs2fP5uDBg9xxxx2eY955550cPXqUd999t+7v\nQBh6RUQpFRV7ju5h3d51njL/l72/W2bCBEhLi35b0tJsgPPGG95z33JLxfY3Tv0GD8+pSCtbt3cd\ne47uoWMr/0wCqiErM2UcKDgQs/O1a9GOBIneb/yCggIKCwvJy8tj+vTpvP/+++XDZnNzc9m7dy+D\nBw8O2m/o0KG8/37FWkqrVtklF/zDa7Ozs0lISGDVqlVcc801UWt3JDQQUUpFxawtszzbrVJaMTRz\naPl2qG6Zq66qu/ZcfbU3EJkzBw4fhjZt7PaQLkNoldKKoyeOlteZtWUWk/rpxIWNyYGCA3R4tEPM\nzrf3+3tp3zJ6SU/3338/zzzzDAAJCQlcccUV/PnPfwZg1y47G0Xnzp2D9uvcuTMHDx6kuLiY5ORk\ndu3aRWJiIhkZ3rVkk5OTadeuHbm5uVFrc6S0a0YpFRX+bpnR3Ud7ZiudMSM23TIB55/vXcW3pMSb\nN5KcmMy5J5/r2UeH8ap4M2XKFD766CNeeuklLrjgAkpLSzl+/DgAhc7ws2Yhlqtu7vRNBuoUFhaS\nkpIS8hzNmzcvr1cfNBBRSkVFVfkhM2Z469dVt0xAq1Yw1tuE4GG8vjZqIKLiTc+ePRkzZgyTJ09m\nxowZ5Ofnc8kllwCQmpoKUB6YuBUVFXnqpKamcuLEiZDnKCoqKq9XHzQQUUrV2rbD29h6eKunzD1X\nx/Hj8KF3+RkuvbTu23XZZd7t99+3bQkYc8oYz+1bD29le972um+YUjV05ZVXsmzZMj7//PPyLplA\nF43brl27aNu2LcnOctadO3emtLQ0aH6R4uJiDhw4QJcu9bcerOaIKKVqzT8HR7vUdvTt0Lfi9nlw\ntCIVAxF7RaSuXXyxPZcxdjs/3+aKjB9vt/t26Evb1LYcLDxYvs/8bfO5tv+1dd84FRPtWrRj7/f3\nxvR8dSnQhZKXl0ePHj1o3749y5cvD6q3dOlSBgyoWMpgwIABGGNYvnw5559fsd7SsmXLKCsr89SN\nNQ1ElFK15g9ERnYf6Rk54J8McuhQ6BiDwSmdOsHZZ8PixRVlb79dEYgkSAIju41k+qaKPpu52+Zq\nINKIJEhCVJNHY2Xfvn209830V1JSwosvvkhqaip9+vQB4IorruCll14iJyenfAjvxx9/zObNm7n/\n/vvL9x0zZgxt27blqaee8gQiTz31FC1btuTCC+tvDh0NRJRStTZvuzcQGdVtVPn/jYH//tdbvy6T\nVP0uu8wbiLz7rm1TYG2bUd1HeQIRnWFVxYNbb72VI0eOMGrUKDIzM9m9ezevvPIKmzZtYurUqbRo\n0QKABx98kDfeeINzzz2Xe+65h/z8fB599FHOOussrr/++vLjNW/enF/96lfcddddXHXVVYwfP555\n8+bx6quv8tvf/pY2geFk9UADEaVUrew+upvNBzZ7ykZ1rwhENmyALVu8+8QyELn4Yvjxjyu2d+yw\nbXJ+UHraCrDpwCadT0TVu4kTJ/Lcc8/x9NNPc+DAAdLS0sjOzuaRRx7xXL3o2rUrc+fO5b777uOB\nBx4gJSWFiy66iEcffbQ8PyTg9ttvJyUlhccee4z//ve/ZGVl8cc//pHvfe97sb57HhqIKKVqZcH2\nBZ7ttJQ0zup0Vvm2v1uma1c46yxipndvyMqyAUjABx9UBCIDOg0gLSWN/BMVi97N3z6fK/tcGbtG\nKuVz1VVXcVU1J9rp3bu3Z/Kyytx0003cdNNNtWla1OmoGaVUrfi7Ms7pdg5JCRW/cfwzR190UUW3\nSCyI2DlF3D74oOL/SQlJnNPtHM/tc7fOjUHLlFKggYhSqpb8gYg7PyQ/HxYu9Navj5w4fyAydy4c\nO1ax7W4zBOe8KKXqjgYiSqkaO1R4iLV71nrK3DkXc+d6Z1NNToZzz41R41zGjoXExIrtEyds2wL8\neSLr9qzzDOlVStUdDUSUUjX2yY5PMJjy7eZJzRncpWIBLv/aMsOH2xlPYy09HUaM8Ja5u2eGZA6h\neVLFcr0GE5T7opSqGxqIKKVqzN8tM6zrMJolVax74Q9EzjsvFq0KrbI8kZTEFIZ3He65XYfxKhUb\nGogopWps/vb5nm13rsXOnXaYrFs8BSKffw5ffVWx7e+e0UBEqdjQQEQpVSPHThxjea53amn3l/lH\nH3nrt2kDgwdTbwYMAN9ElcyaVfF/fyCyavcqCooLYtAypZo2DUSUUjWyJGcJJWUVmahJCUkM6zqs\nfNvfLTNmjDdhNNYSEmwb3NyByNDMoSRKRQNLykpYlrMsRq1TqunSCc2UUjWycId3XO7ATgNpmdIS\ngLKy4Csi9dktEzBmDLz2WsX2rFkV0723SmnFWZ3OYuWuleW3L9yxkNEnj66Hlqrq2uDv/1N1oi4f\nZw1ElFI1smjnIs/2iKyKYSnr1sFe34Kn8RCIjB3r3d6zB9avhzPPtNsjuo7wBiI7fZOgqLiRkZFB\nixYtmDx5cn03pclo0aIFGRkZUT+uBiJKqYiVmTIW7QgfiPi7ZU45BU47LRYtq9ypp0K3brB9e0XZ\nrFmuQCRrBE8ue7L8toU7FmKMQWI5Fayqlm7durFhwwb2799f301pMjIyMujWrVvUj6uBiFIqYpsP\nbOZQ0SFPmXv46+zZ3vrxcDUEbBfMmDHwwgsVZbNmQWDNL/9U7wcLD7L5wGZ6ZfSKXSNVtXXr1q1O\nvhhVbGmyqlIqYv78kK6tu5KVngXYmVQX+OYC8yeJ1id/W+bMgdJS+/+s1llkpmV6bvffV6VUdGkg\nopSKmP/L2d0ts3o1HDnirT86jvI9/YHI4cOwapX9v4h47gvY2WOVUnVHAxGlVMT8iarubpk5c7x1\ne/WCTp1i0Khqysy0bXJzD+P1ByJ6RUSpuhU3gYiI3CkiW0SkUEQWi8iQSuqeIyILRGS/iBSIyAYR\nuTdEvW85txWKyBoRmVC390Kpxu9Q4SHW71vvKXN/ebsXk4P4uhoSUNl8Iv5AZMP+DboAnlJ1KC4C\nERG5GngMeBgYCKwBZopIuHFCx4A/AyOBM4BfAb8Wke+6jjkCeBX4GzAAmA68LSJ96up+KNUULN65\n2LPdPKk5AzoNAGyuxTzfzOj1sdpuVfyByPz5UFxs/z+w00BSk1I9t/vvs1IqeuIiEAGmAM8YY14y\nxmwEbgMKgBtDVTbGrDbGvGaM2WCM2W6MeRWYiQ1MAu4G3jfGTDXGbDLGPASsBO6q27uiVOPm75YZ\n3GUwKYkpAKxZE9/5IQH+4KigAFY604ckJyYzJNN7QVa7Z5SqO/UeiIhIMpANfBwoM8YY4CNgeLj9\nfMcY6NSd4yoe7hzDbWZ1j6mUCi0oUbVrRVeGPz+kRw/o0iUGjYpQRgb08V0bdV/Jcd8n0IRVpepS\nvQciQAaQCOzxle8BKk1xE5EdIlIELAX+Yox53nVzp5ocUykVXmlZKUtylnjKhmdVxPb+/JB47JYJ\nGOVd484biPjyRJbmLKW4tDgGrVKq6YmHQKQ2voa9mnIbMMXJNVFK1ZFP937K0RNHPWWBETOh8kPi\nsVsmwB+ILFhg18gBb3AFUFBcwNo9a2PUMqWalniYWXU/UAp09JV3BHZXtqMxZpvz389EpBPwcyCw\npNXumhwTYMqUKaSnp3vKJk2axKRJk6raValGzd8tc+pJp9KxlX2brVtn5+Rwi+dAZORI7/bhw/Dp\np9C/P2S0yKBnu55sPrC5/PYlOUvI7pId41YqFZ+mTZvGtGnTPGV5eXk1Ola9ByLGmGIRWQGMBWYA\niF3YYSzwRASHSgSaubYXhTjGeU55pR5//HEGDRoUwamVahoqW+jOnx9y+unQtWsMGlVDXbvatWe+\n+qqibN48G4gAnJ15ticQWbxzMXcMuSPGrVQqPoX6cb5y5UqysyMP1uOla2YqcLOIXCciZwBPAy2A\nFwBE5Hci8mKgsojcISIXicjpzt9NwP3AP13H/BNwvojcJyK9ROTn2G6cJ1FK1Yj/ioh7IrOGMH+I\nX2V5Imdnnu25zZ8bo5SKjrgIRIwxrwPfB34JrAL6A+ONMfucKp2ALNcuCcDvnLrLgNuBHxhjHnYd\ncxFwDXALsBq4HLjUGOOdiUkpVS17j+3ly0NfesoCV0TKyhrG/CF+oQIRY+z/z+7qDUQ2H9jMoULv\nQn9Kqdqr966ZAGPMX4G/hrntBt/2k1TjyoYx5k3gzag0UKkmblnOMs92i+QW9O3QF4ANG+Cgb/LR\nhnhFZM8e+Pxz6NkT+nfsT7PEZhwvPV5++9KcpYw/fXyMW6lU4xYXV0SUUvFvac5Sz/bgLoNJSrC/\nZT7xTbPRrRtkZRH3Tj01eJ6TwJWdlMSUoORU7Z5RKvo0EFFKVcuyXO8VkSFdKmYf9Qci55wTixbV\nnkjw6BnNE1EqtjQQUUpVyRgTdEVkaObQ8v831EAEIkxY3bkEE0giUUpFhQYiSqkqbTm8hQOFBzxl\ngUBk92740pvD2qADkW3b7B8EJ6weKDwQlLCrlKodDUSUUlXyXw1p36I93dO7A8FXQ9LSoF+/WLWs\n9vr0gbZtvWXz59t/u6d3p0PLDp7bluzU7hmlokkDEaVUlfwjZoZkDsHOOxgciAwbBomJsWpZ7SUk\nBOeJBO6TiGieiFJ1TAMRpVSVlub68kO6hM8P+drXYtGi6PJ3JbnvkwYiStUtDUSUUpUqKSthRe4K\nT1kgP6SgAFau9NZvSPkhASO8i+3y6acV6+b480RW717N8ZLjKKWiQwMRpVSl1u9bT2FJoadsSKYd\nurtsGZSUVJQnJsLZ3u/tBiE7G1JSKraNgcWL7f+HdBmCIOW3nSg9wardq2LcQqUaLw1ElFKV8ieq\nntLmFDJaZACwYIG37llnQatWsWpZ9DRvDoMHe8sC3TPpzdPp3b635zZNWFUqejQQUUpVqrHOH+Kn\neSJK1Q8NRJRSlQoXiJSVwaJF3rqNKRBZsgSKi+3/NRBRqu5oIKKUCquguIBP937qKQsEIuvXVyR0\nBjTkQMSfsFpQAGvW2P/7E1a/OvQV+47tQylVexqIKKXCWrVrFaWmtHw7QRIY2GkgENwt0707dO0a\ny9ZFV/v2dtVdt8B97NuhLy2SW3hu818pUkrVjAYiSqmw/F+2fTv0pWVKS6Bx5YcEhMsTSUpIIruz\ndyXexTsXx6hVSjVuGogopcLyr7jrnsjMP2KmsQYigTXu/Hki/sdGKVUzGogopcIKl6i6axds2eKt\n2xgDkdzcigXwAnOnBCzLXaYr8SoVBRqIKKVCOlAQvNJs4MvY3y3TujX07RurltWdXr2gXTtvWeC+\nuoctAxwsPMiWw75oTCkVMQ1ElFIh+bseUpNSObP9mUDDX+guHJHg0TOB+9o9vXv5RG4B/sUAlVKR\n00BEKRWS/0t2UOdBJCcmA40zUTXAf18WLrT/ighDuni7Z3TkjFK1p4GIUiok/4q7gS/hY8dglW+p\nlcYciKxbB0eO2P/7AxFNWFWq9jQQUUoFMcaETVRdurRxLHQXTnY2JCdXbJeVuRbA8yWsrty1ktKy\nUpRSNaeBiFIqyPa87ew9ttdTFghE/N0yAwY0zIXuwklNtcGIW+A++6+IHCs+xob9G2LUMqUaJw1E\nlFJB/F0ObVPbcupJpwKNOz8kINzEZh1bdaRbejfPbZonolTtaCCilAri/3Id0mUIItLoFroLx3+f\nFi+u6I4KyhPRkTNK1YoGIkqpIOHyQz77DPLyvHUbYyDiH8J77BisXWv/rwmrSkWXBiJKKY/SslJW\n7FrhKQuXH9K9O2RmxqplsdOxI5x+urcs3MRma/es5XjJ8Ri1TKnGRwMRpZTHxv0bOXriqKcscBWg\nMa4vE064PJHsLtkIUl5eXFbMmj1rYtgypRqXuAlEROROEdkiIoUislhEhlRS95si8qGI7BWRPBFZ\nKCLjfHW+IyJlIlLq/FsmIgV1f0+Uatj83TLd0rvRsVVHIPiKyNe+FqtWxV64QKR1s9b0yujluU0T\nVpWqubgIRETkauAx4GFgILAGmCkiGWF2GQV8CEwABgGzgf+KyFm+enlAJ9df9+i3XqnGJVx+SG4u\nbN3qrduUrojs3Anbt9v/a56IUtETF4EIMAV4xhjzkjFmI3AbUADcGKqyMWaKMeZRY8wKY8yXxpif\nAJ8DFwdXNfuMMXudv311ei+UagT8X6pDu4TOD2ndGs48M1atir0zzoCTTvKWBbqm/HkiOnJGqZqr\n90BERJKBbODjQJmxa2t/BAyv5jEESAMO+m5qJSJbRWS7iLwtIn2i1GylGqWikqKgfIdwiarDhzeO\nhe7CSUgIvwCe/4rIxv0bOXL8SIxaplTjUu+BCJABJAJ7fOV7sN0p1fEDoCXwuqtsE/aKyiXAtdj7\nulBEutSqtUo1Yqt3r6akrGL+dkEY1HkQ0DQmMvPz58AEHoOzOp1FUkJSebnBsCLXO9JIKVU98RCI\n1IqIXAP8DPiWMWZ/oNwYs9gY87IxZq0xZj5wObAPuLWemqpU3PPnh/Rp34e0ZmmNfqG7cMItgNc8\nqTn9O/b33KZ5IkrVTFLVVercfqAU6Ogr7wjsrmxHEZkIPAtcaYyZXVldY0yJiKwCTq+sHsCUKVNI\nT0/3lE2aNIlJkyZVtatSDVpQfojTLbNkCZS61nZLTISh3jSJRmnwYLsAXnGx3Q4sgDdunO2eWblr\nZXldDURUUzJt2jSmTZvmKcvzz3ZYTfUeiBhjikVkBTAWmAHlOR9jgSfC7Scik4C/A1cbYz6o6jwi\nkgD0A96tqu7jjz/OoEGDqncHlGpEwo2YaewL3YUTWAAvsPou2Mdi3Dj72Dyz4pnyck1YVU1JqB/n\nK1euJNu/YmQ1xEvXzFTgZhG5TkTOAJ4GWgAvAIjI70TkxUBlpzvmReB+YJmIdHT+Wrvq/ExEzhOR\nU0RkIPAK0A0bvCilfA4XHWbzgc2eskBSZlOaP8Qv3Hwi/oTVbXnbglYsVkpVLS4CEWPM68D3gV8C\nq4D+wHjXcNtOQJZrl5uxCa5/AXJdf3901TkJ222zHnsVpBUw3BkerJTyWZ673LPdLLEZ/Tr2o7S0\naSx0F064BfB6t+9Ni+QWntv0qohSkYuLQATAGPNXY8zJxphUY8xwY8xy1203GGPGuLa/boxJDPF3\no6vOfcaYU5zjdTHGXGyMWRvr+6VUQ+HvlhnYeSApiSl89plN0HRrSoFIqAXw1qyBpISk8hFFAZon\nolTk4iYQUUrVL38gEq5b5uSToUsTGgTfsSP06OEtK18Ar4tvYjMNRJSKmAYiSimg+omqTelqSEDY\nPJFM31TvOcuw8zEqpapLAxGlFDlHcth1dJenTAORCqECEWOCE1b3Fexje972GLZMqYZPAxGlVNDV\nkDbN23B629Ob3EJ34fjvc06OXQDv1JNOpW1qW89tuhKvUpHRQEQpFfTlObjLYBIkocktdBdOr17Q\n1htv8MknICK6Eq9StaSBiFKq2ivuNvaF7sKJZAE8DUSUiowGIko1cWWmLOzU7k15IjO/6iasrshd\nQWlZKUqp6tFARKkmbvOBzUFL2A/NHNpkF7oLx3/f166FvLzgKyL5J/LZdGBTDFumVMOmgYhSTZw/\nPyQzLZPOaZ2b7EJ34QwZAikpFdvG2FlWO6d1JjMt01NXZ1hVqvo0EFGqifN/aYbrlhk4EFq2jFWr\n4k/z5nYBPLfyic0ydWIzpWpKAxGlmriluTqRWXVVdwE8DUSUqj4NRJRqwo6XHGf17tWesqGZQ5v8\nQnfh+B+DJUvsAnj+hNXVu1dzovREDFumVMOlgYhSTdjaPWuDvjCzO2c3+YXuwgm3AN7gLoM95SdK\nT7B2j66xqVR1aCCiVBPmT1Q9I+MM0punN/mF7sLp0CH0AnhtmrehZ7uennJNWFWqejQQUaoJ0/yQ\nyFU3T8T/2CqlQtNARKkmLGjETJgZVTUQqeCf1G3BgtAL4OkVEaWqRwMRpZqovKI8Nu7f6Ckbmjk0\n5EJ3TXlGVT9/UJabC9u2wdldz/aUr9+3nvzj+TFsmVINkwYiSjVRK3atwGDKt1MSU+jfsX/Q1ZD0\n9Ka50F04vXpBu3besk8+gQGdBpCUkFReZjAsz10e49Yp1fBoIKJUE+VPVB3QaQDNkpqxYIG33vDh\ndtE3ZYmEXgCveVJzzup4lqfc/xgrpYLpx4tSTZT/SzKQHzJ/vree5ocEC5ew6p9hVRNWlaqaBiJK\nNVFBgUjmUPLyYLV3fjNGj45hoxoIfyCybp1dAO/sTG+eyJKdS2LYKqUaJg1ElGqCco7kkJOf4ykb\nmjmUTz6xI0ACmjWzi70pr8GDQy+A578ikpOfQ86RHJRS4WkgolQT5F8LJb1ZOj3a9WDePG+9oUPt\nYm/Kq3lzG4y4LVgAvTJ60bpZa0+5rjujVOU0EFGqCfJ3ywzJHEKCJAQFIqNGxbBRDYx/SPO8eZAg\nCUHziWj3jFKViygQEZExIpJUdU2lVDwLlahaUADLfaNNNRAJz//YLFkCRUWasKpUpCK9IvI/oG1g\nQ0QWi0hmdJuklKpLZaYsqLtgaOZQliyB4uKKssREO3RXhXbOOXYob8Dx47B0aXAgsixnGWWmLMat\nU6rhiDQQEd/2mUCzKLVFKRUDmw9s5shx79K6QzOHBnXLDBoEaWkxbFgD06YNDBjgLZs7N3jkTP6J\n/KAZbJVSFTRHRKkmxt8t07V1VzqndQ4KREaOjGGjGij/0Oa5c6FzWme6tu7qKdeJzZQKL9JAxDh/\n4bZrTETuFJEtIlLodPmEHTQoIt8UkQ9FZK+I5InIQhEZF6Let0Rkg3PMNSIyIRptVaohC0pU7TKE\nEydg0SJvPc0PqZr/MVq4EE6cCJEnooGIUmHVpGvmYxFZKSIrgRbAfwPbrvLIDipyNfAY8DAwEFgD\nzBSRjDC7jAI+BCYAg4DZTjvK51cWkRHAq8DfgAHAdOBtEekTafuUakxCTWS2ciUUFnrr6UJ3VfNf\nNSoshBUrKmapDViSoyNnlAon0hEwv/BtT49SO6YAzxhjXgIQkduAC4EbgT/4KxtjpviKfiIilwIX\nY4MYgLuB940xU53th0TkPOAu4I4otVupBuV4yXFW7/ZOnTo0cyjz/u2t17dv8MJuKlhGhn2sPv20\nomzePDj7am+eyNo9ayksLiQ1OTXGLVQq/kUUiBhj/IFIrYlIMpAN/NZ1HiMiHwHVytkXEQHSgIOu\n4uHYqyxuM4FLa9VgpRqwNXvWUFxWMTRGELI7ZzNV80NqbNQobyAydy7ccW82gpSvblxSVsLq3asZ\nnqXDkJTyq3GyqohkiMhgEckWkdr8dsoAEoE9vvI9QKdqHuMHQEvgdVdZp1oeU6lGx98tc0bGGbRK\nTg9acVfzQ6rPn7C6YAGkJqbRp723F1i7Z5QKLeJARETOFJF52C/1JcBSYK+IzBKRXtFuYDXacw3w\nM+Bbxpj9sT6/Ug1JqPyQTz+1C7a56RWR6vMHbfn5sGZN8DBeTVhVKrSIumZEpBMwF9gH3AdsxCaw\n9gFuBuaLSF9jzN4IDrsfKAU6+so7AruraM9E4FngSmPMbN/Nu2tyTIApU6aQnp7uKZs0aRKTJk2q\nalel4lqoQMQ/bPe00yBTpymstk6doGdP2Ly5omzuXBg6aij/WP2P8jINRFRjMm3aNKZNm+Ypy/P/\noqmmSJNVpwDbgHOMMUWu8g9E5ClggVPngeoe0BhTLCIrgLHADCjP+RgLPBFuPxGZBPwduNoY80GI\nKotCHOM8p7xSjz/+OIMGDaruXVCqQThcdJhNBzZ5yoZmDuX3c7319GpI5EaP9gYi8+bBw5O8I2e+\nPPQl+wv2k9Ei3GBApRqOUD/OV65cSXZ2dsTHirRr5jzg/3xBCADGmELgEWB8xK2AqcDNInKdiJwB\nPNwHnGYAACAASURBVI0dGvwCgIj8TkReDFR2umNeBO4HlolIR+fPvezln4DzReQ+EeklIj/HJsU+\nWYP2KdXgrchd4dlOSUyhb/v+zJnjrefPeVBV83fPzJ8PfTL6kprkHSWzLEdX4lXKL9JA5FSgsnlC\nljt1ImKMeR34PvBLYBXQHxhvjNnnVOkEZLl2uRmb4PoXINf190fXMRcB1wC3AKuBy4FLjTHrI22f\nUo2BP1lyYKeBbFqfwoED3npf/3oMG9VI+IO3gwdh04ZkBnX2XlnV7hmlgkXaNZMGHKnk9nygVU0a\nYoz5K/DXMLfd4Nuu1kelMeZN4M2atEepxmbxzsWe7SFdhjDbl1l16qnQvXsMG9VIZGXBKafAli0V\nZfPmwdAeQ/lkxyflZboSr1LBajJ8N01EWof6A1oTvDCeUqqeGWOCApHhWcOZNctbT6+G1Jy/e2bu\n3OCp3pfsXIIxUVkVQ6lGoyZTvG8GDoX52xR+V6VUffnq0FfsK9jnKRvSeVjQiBkNRGrO3z0zezYM\n6ewdwnug8ABbDm9BKVUh0q4Z/ZhSqgHyXw3p0LIDh7ecEjR/iAYiNTdmjHf7wAHI33EyGS0y2F9Q\nMcXR0pylnHpSxKl0SjVakV4RWQCcjZ2O/VHsonNLjTFz3X/RbqRSqnb8gciwrsOYPdvbi9qrF3Tp\nEstWNS7du9scG7dZsyRk94xSqkKkgciD2CAkH8gB7sGOXFFKxbFFO73T5wzvOjwoUVWvhtTe2LHe\n7Y8/Dp5hdXGONyhUqqmLNBC5DrjDGHO+MeYy7Gq314pIjdesUUrVrYLiAtbsWeMpG9xpGPPne+v5\nuxZU5PyByLx5kN3JG4is3LWS4yXHY9gqpeJbpAFEN+D9wIYx5iPAAHpBV6k4tSJ3BSVlJeXbCZIA\nuYM5dsxb79xzY9uuxsh/VenoUUja7Q1ETpSeYOWuyqZjUqppiTQQSQL8s6oWA8nRaY5SKtr8+SH9\nO/ZnyXzvdD99+0L79rFsVePUoQP06+ctWzqvDWe2P9NT5u8qU6opi3TUjAAviIj7umJz4GkRKf99\nZYy5PBqNU0rVnv9Lb1jmMGa97q2j3TLRM3YsrFtXsT1rFgyfMpzP9n1WXqaBiFIVIr0i8iKwF8hz\n/b2MnV7dXaaUigPGmKAvveyOw1i40FtPE1Wjxx/ULVwIgzoO95btWKgTmynliOiKiH+qdaVUfNtx\nZAe7j+72lKXsHU6Rq4NVRBe6i6bRoyExEUpL7faJE5Cye4SnTm5+LjuO7KBberd6aKFS8UVHuyjV\niC3a4b0aclLzk9jwSQ9P2aBBcNJJsWxV49a6NQwZ4i3bvKgnJzX3Psj+50appkoDEaUasVATmf3v\nQ+9EZuedF8sWNQ3+7pnZsxIY1nWYp0zzRJSyNBBRqhHzf9md1XY4K30jR8eNi2GDmgj/fCIrVsDA\nDG+eiAYiSlkaiCjVSB0vOc6q3as8ZQm7huHOkWzRAkaMQEXZiBHQrFnFdlkZNNvrfaBX7lpJYXFh\njFumVPzRQESpRmrlrpWcKD1Rvi0IOxZ51z0ZPdr7hamio3lzOOccb9n2xUPtZHKOkrISVuxaEeOW\nKRV/NBBRqpHy54f0ad+HOTPTPWXaLVN3/I/trPfT6Nuhr6dME1aV0kBEqUbLn4PQO20YO3Z462gg\nUnfOP9+7vWULnJnm7Z7RPBGlNBBRqlEyxvDJjk88ZUm7vaM2unSB3r1j2aqmpX9/6NTJW5a4Wyc2\nU8pPAxGlGqFtedvIzc/1lO1e5k1aGDfOTmam6oZI8BWnnQu9gcieY3vYenhr7BqlVBzSQESpRuiT\n7d6rIe1S27Hs/TM8ZdotU/f83TNLPjidjNQMT5l2z6imTgMRpRqhBdsXeLbPaHkOx456L3/457pQ\n0Xfeed6rToUFQo9U33wimrCqmjgNRJRqhPz5Ic32eLtlBg60S9arupWRAYMHe8uSduvEZkq5aSCi\nVCNzuOgwn+791FO245Oveba1WyZ2xo/3bucs9gYiq3ev5tiJYzFskVLxRQMRpRqZRTsWYagYidEs\nsRmfz8321Lnggli3quny54l8NX8ISVKx8HmpKQ2a80WppkQDEaUaGX9+SLfEIVBaMX1qejoMH+7f\nS9WVs8+2j3m54pZ0Sx7kqTN/+/zYNkqpOKKBiFKNjD8/hJ3e/JDx4yE5OYYNauKSkuAb3/CWJe8a\n6dnWQEQ1ZRqIKNWInCg9wZKcJZ6yHQu8+SHaLRN7/jyRnZ94A5HFOxdTXFocwxYpFT/iJhARkTtF\nZIuIFIrIYhEZUkndTiLyiohsEpFSEZkaos53RKTMub3M+Suo23uhVP1atWsVRSVFnrKiz73Tik+Y\nEMsWKQh+zI9t8AaHBcUFrNy1MoYtUip+xEUgIiJXA48BDwMDgTXATBHJCLNLM2Av8CtgdSWHzgM6\nuf66R6vNSsUjf35I29I+UNi2fHvIEB22Wx+6doUBA1wFhe1oV9rHU0e7Z1RTFReBCDAFeMYY85Ix\nZiNwG1AA3BiqsjFmmzFmijHmZeBIJcc1xph9xpi9zt++6Dddqfjhzw8p3eL95X3hhbFsjXK7+GLv\ndslX3udGAxHVVNV7ICIiyUA28HGgzNhVoD4Capvb30pEtorIdhF5W0T6VL2LUg2TMSboikjeOm+i\nquaH1B9/IJK3zpsnsmD7AspMWQxbpFR8qPdABMgAEoE9vvI92O6UmtqEvaJyCXAt9r4uFJEutTim\nUnHri4NfsK/Ad9Fve8Wv7o4dITsbVU+ys32r8W7zBiIHCw+ycf/G2DZKqTiQVHWVhskYsxgonyVI\nRBYBG4BbsbkoYU2ZMoV0z8B/mDRpEpMmTaqDlioVHf6rISknOnHi0Cnl2xMmQEI8/PRoohISbNfY\nc885BXndaVaUxfHmO8rrzN82nz7t9cKtin/Tpk1j2rRp/9/efcdHVaWPH/88KZAAEkCEgFKkSocE\nV8GCDSlSBJQmgqgsKq6Kurqu7rJf+89VRFR00UXFgiC4giBFUVyxoEsookgNvRcDBEJCcn5/nEmY\nOzOpTObOJM/b17yYe+bcO0/GZOaZUx1laWlpJbpWOCQiB4BsoLZPeW1gT7CexBhzSkRWAE0Kq/vi\niy+SlJRUWDWlwsrXW792HGdtvBQ4veOadsu4r3dvr0QEyNxwGbT5IO/4m23fMLrjaBciU6p4An05\nT0lJIbkEza6ufz8yxmQBy4G8vUBFRDzH3wXreUQkCmgD7A7WNZUKJ0u2LHEcmy1d8u5XqOC/1LgK\nvWuugYqnF7nFbNWFzZRyPRHxGA+MEpHhInIB8DpQCXgbQESeEZF3vE8QkXYi0h6oApzjOW7h9fjf\nRKSriJwvIh2A94H6wJuh+ZGUCp0tv29ha9pWn8Ir8u5ecw2cdVZoY1L+KleGq67yKvBJRLalbWNb\n2rbQBqWUy8IiETHGzAAeBB4HVgBtgW5e020TgXo+p63AtqQkAUOBFGCe1+PVgcnAr57yKkAnz/Rg\npcqUr7c4u2U4fjbsPz3W4PrrQxyQypdj9syBFojXOi9gx4koVZ6ERSICYIyZZIxpaIyJN8Z0Msb8\nz+uxkcaYq3zqRxljon1ujbwev98Yc77nenWNMb2NMatD+TMpFSq+40PYejkY++ctAn36uBCUCqhX\nL68DE4XZ5pxird0zqrwJm0REKVVyvuNDvLtlOne2U3dVeKhXz2catU/3jF9SqVQZp4mIUhFue9p2\nUn9PdRZ6DVTVbpnwM2CA18HWLo7HfjvwG7uP6ph6VX5oIqJUhPP7Bn2iOuxrk3fYt2+IA1KFciQi\nu5Mgo6rjcb8WLqXKME1ElIpwfh9aXuNDWrWCpk1DH5MqWLNm0Lq15yAnxv4/8/Jl6pehD0opl2gi\nolSE82sR0W6ZiOBoFdlypeOxr7Z8FdpglHKRJiJKRbCdR3ay8dBGZ6HXQNV+/UIbjyo6RyKS6kxE\nNh3exPa07ShVHmgiolQE82sNyUiAvW0BaNwYdKeC8NW6tVe32d52dmyPF20VUeWFJiJKRbDA40Oi\nARg40K4hosKTiFeriIlydKmBJiKq/NBERKkI5vdh5fVhNmhQiINRxebsnnGs2ciXqV9ijAltQEq5\nQBMRpSLUlt+3BBgfYscaNGsGbdu6EJQqluRkaNDAc+AzYHVb2jb/9WGUKoM0EVEqQi3evNhZcPxs\n2NMesK0h2i0T/kRg8GDPwb5WkH6O4/GvUrV7RpV9mogoFaG+SP3CWbD56rz1QwYOdCEgVSJDh+be\nE8eMJ9BxIqp80EREqQiUY3L8W0RSrwagZUuvxbJU2GvTxv4/A/ym8X615SsdJ6LKPE1ElIpAa/at\nYf/x/c7CzdcA2hoSaUS8WkV8xonsOrqLdQfXhT4opUJIExGlItAXm326ZQ6fD4cbAV5jDlTEGDLE\nc+dAczhS1/HY55s+D31ASoWQJiJKRSC/RMTTGnLhhdC8uQsBqTPSqBFcdBGAwOaujscWbV7kSkxK\nhYomIkpFmMzsTP8VVT2JyPDhLgSkgiKve2ZjN0f5V6lfkZmdGfqAlAoRTUSUijA/7PiB41nHnYWp\nVxETo90ykWzgQIiKIi+pzJWelc732793JyilQkATEaUijF+3zO4OcLwmPXpAzZruxKTOXGIidO0K\nHD8Hdjk3CVq4aaE7QSkVApqIKBVh8hsfot0yke+WWzx3Njm7ZxZt0nEiquzSRESpCHLoxCGW7Vzm\nLNx8NQkJ0KuXOzGp4Ln+eqhWDdh0raM8ZXcK+9P3Bz5JqQiniYhSEeTzTZ+TY3JOF2TFwdbLGTQI\n4uLci0sFR1ycZ9Dq9k6QWTmv3GBYnLo4/xOVimCaiCgVQeZvnO8s2HIlnIpnxAh34lHBN3IkkF3R\nb7l37Z5RZZUmIkpFiByTw4KNC5yFG3rQsiV06uROTCr4kpM9S/Rv9B8nosu9q7JIExGlIsTKPSvZ\nm77XWbixB6NG6U67ZYmIp1XEZ5zIzqM7+XX/r+4EpVQp0kREqQgxf4NPt8zBJlQ41oSbb3YnHlV6\nhg2D6LRm8Ht9R7lf15xSZYAmIkpFCL8PoY09uOEGOPtsd+JRpadWLbi+r8DGHo7yuevnuhSRUqUn\nbBIRERkjIqkickJEfhCRCwuomygi74vIOhHJFpHx+dS7UUTWeq65SkR6BKqnVLg7fOKw/+qaG2y3\njCqb7roLWO+ck/3N1qUcPnHYnYCUKiVhkYiIyCDgBWAc0AFYBSwUkfzWiawI7AOeAFbmc83OwAfA\nG0B7YDbwiYi0DG70SpW+zzd/Tg7OabuNY66gSxf3YlKl68oroWnMVXaKtkcO2brKqipzwiIRAcYC\n/zLGTDXG/AbcARwHbg1U2Riz1Rgz1hjzHnAkn2veA8w3xow3xqwzxvwdSAHuLoX4lSpVs1b7Ttu9\ngtG3xusg1TJMBMb8sRKkXuUon7lqnksRKVU6XE9ERCQWSAbyVusxdo7aF8CZTErs5LmGt4VneE2l\nQi47J5t56z5zlMVu7cFtt7kUkAqZESMgNtXZPfPZxs/Izsl2KSKlgs/1RASoCUQDPvMS2QsknsF1\nE0vhmkqF3Depy0iXfY6y/q2vo0YNlwJSIVOtGvRvc52j7ASHWLrlB5ciUir4wiERUUoV4Pm5s50F\n+1rx6F2N3QlGhdxDo+vDnraOsvFzdfaMKjti3A4AOABkA7V9ymsDe87guntKes2xY8eSkJDgKBsy\nZAhDhgw5g3CUKj5jYPH2OXDW6bL6J/rSpo17ManQSkqCc9N7sZPVeWWLts7FmGd0jJByzbRp05g2\nbZqjLC0trUTXcj0RMcZkichy4GpgDoCIiOd44hlc+vsA1+jqKS/Qiy++SFJS0hk8tVLBMWPxejLO\n+s1RNubqvi5Fo9xy5zW9eCz16bzjjKprmL5wC4O7N3QvKFWuBfpynpKSQnJycrGvFS5dM+OBUSIy\nXEQuAF4HKgFvA4jIMyLyjvcJItJORNoDVYBzPMctvKq8BHQXkftFpLmI/AM7KPaV0v9xlAqOJz5y\ndstEH6/DfQM7uhSNcsufh/yBqAznagZPzpydT22lIktYJCLGmBnAg8DjwAqgLdDNGLPfUyURqOdz\n2gpgOZAEDMVOzc2b12aM+d5T/kfsWiP9gb7GGN2sQUWE336DX045P2w6Vu1Nhdiw+LNVIVQhNpqO\nVZ2zZ37J/ph161wKSKkgCpt3NGPMJGNMQ2NMvDGmkzHmf16PjTTGXOVTP8oYE+1za+RTZ5Yx5gLP\nNdsaY3QlIBUxHn9hH9T7zlH2YC/tlimvHujZ31nQ4Bueesl3YqBSkSdsEhGl1Gk7dsCMFfNATm/7\nHmsq06vVVQWcpcqyPq26UsFUOV0ghmkps9m9272YlAoGTUSUCkMTJkB2s48dZdc27k5cTFw+Z6iy\nLi4mju6NnWuKnGo6i+efdykgpYJEExGlwsyhQ/D6279D40WO8sHtrncpIhUubkry6Z45/0tee+sw\n+/cHrq9UJNBERKkwM348pJ83B2Iy88oqRFWgd7PeLkalwkHPpj2pEFXxdEH0KU7Um8v4gPuPKxUZ\nNBFRKowcOAAvvQS0muEo7960OwlxCYFPUuVGlQpV6N60m7Ow5SxeeQUOHnQnJqXOlCYiSoWR55+H\nY6f8u2VubHmjSxGpcNP/Ap/umcYLOZZ5zCawSkUgTUSUChP79sHLLwPNZ0N0Vl55xeiK9Gnex73A\nVFjp3bw3MVFei2LHZkDzOUyYgI4VURFJExGlwsRzz8Hx40Crjxzl3Zp0o2rFqu4EpcJOjfgaXNPo\nGmdhmw84ehSefjrwOUqFM01ElAoDO3bAq68Ccdotowo3tPVQZ0HjhRB/kEmTYOtWd2JSqqQ0EVEq\nDDz2GGRkAC0+1m4ZVajrL7jeuaZM9CloOZPMTBg3zr24lCoJTUSUctnKlTB1queg7buOx7o36a7d\nMsrPWRXP8p/O3eYDwP4urVnjQlBKlZAmIkq5yBh48EH7Lwnb4PwljseHtR3mSlwq/A1t49M90+Ab\nqLodY+Chh9yJSamS0EREKRctWACLF3sO2rzveCyhYgK9mvXyP0kpoEeTHiRU9FpbRgy0ng7A/Pkw\nb14+JyoVZjQRUcolmZnwwAO5RwbaObtlBrYaqHvLqHxVjKnIDS1vcBZ6JbNjx8LJkyEOSqkS0ERE\nKZdMmABr13oO6qTAOWsdj9/c9ubQB6Uiil/3TJ2VkLgSgA0b0EXOVETQREQpF2zfDo8/7lXg0xrS\nsFpDLql/SWiDUhGnS4Mu1Ktaz1nYYUre3SeegF27QhyUUsWkiYhSLrj/fkhP9xxEZ+bNeMg1rM0w\nokT/PFXBoqOiuaX9Lc7Ctu9BTAYAx455d/8pFZ70nU6pEFu0CGbO9CpoPgcqO9fmvrmddsuoovFL\nROIP298pjw8/1IGrKrxpIqJUCB07BqNHO8tiL57sOL6s/mU0O7tZCKNSkaxR9UZc2fBKR1nMhVMc\nx3feCUePhjIqpYpOExGlQuiRR2DLFq+Caqlk1f/cUWdU0qiQxqQi360dbnUcZzdcZNel8di+HR59\nNNRRKVU0mogoFSJffw2vvOIsq3f9m47janHV/KdkKlWI/i36O1bgNRjO7/eWo84rr8B334U6MqUK\np4mIUiGQng633eYsi6ucRcYFzib04W2HEx8bH8LIVFlQKbYSQ1oPcZQdbzGZipVO71tkDIwYYbsH\nlQonmogoFQIPPQSbNjnLho6by/6MPY6yUcnaLaNKZnSyc/DR3hO7uPGxTxxlGzfaGVtKhRNNRJQq\nZZ98ApMmOcsuuQS21HL203Q6rxOta7UOYWSqLOlQpwOd63V2lG1NfIWLL3bWe+MNmD07hIEpVQhN\nRJQqRdu3w63OcYRUqgQPvbCaL7d86Sj3/UarVHHdfeHdjuNvtv2Xx15eTeXKznq33w67d4cwMKUK\noImIUqUkOxtuvhkOH3aWT5wIc/ZOdJTVqlyLwa0HhzA6VRYNaDmA2pVrO8rm7H6Vic5fNw4cgKFD\n4dSpEAanVD40EVGqlPzjH3amjLcbb4Q+gw/w/s/OnXbv7HgnFWMqhi44VSZViK7g17L23s/vcf3g\nw/Tr56y7ZAn87W+hi02p/GgiolQp+M9/4MknnWUNGsDkyfBGymQyTmXklcdGxXJHxztCHKEqq0Z3\nHE1MVEze8fGs47y+/DUmT4bzznPWffZZmDMHpVwVNomIiIwRkVQROSEiP4jIhYXUv0JElotIhois\nF5ERPo+PEJEcEcn2/JsjIsdL96dQyu6oO3y4syw6Gt5/H+KrnOTVn151PDa49WASqySGMEJVltU9\nq67fWjQvLXuJygkn+OgjiI111h8+3O7Uq5RbwiIREZFBwAvAOKADsApYKCI186nfEJgLLAbaAS8B\nb4pIV5+qaUCi161BKYSvVJ60NLj+ev+1Gl580c6UeXf1u+w66twO9d6L7g1hhKo8+HPnPzuO96Xv\nY+qqqVx8MYwf76yblga9evmPZVIqVMIiEQHGAv8yxkw1xvwG3AEcB27Np/6dwGZjzEPGmHXGmFeB\nmZ7reDPGmP3GmH2e237/SykVHJmZMGAArF/vLB8+HO6+G07lnOLZpc86HuvSoAvJdZNDGKUqD5Lq\nJHFNo2scZc9//zzZOdmMGQNDnGufsX69/d3NzAxhkEp5uJ6IiEgskIxt3QBs9gB8AXTK57SLPY97\nWxigfhUR2SIi20TkExFpGaSwlXIwBkaNgsWLneVJSfD66yACH/3yEZsOO1c1e/Qy3QBElY6HL3nY\ncbzx0EY+XvsxInasUvv2zvpffQV33WV/l5UKJdcTEaAmEA3s9Snfi+1OCSQxn/pVRSR36sE6bItK\nH+Am7M/6nYjUDUbQSnkbNw6mTnWW1aoFH38M8fGQY3J4eunTjsc71u3o961VqWC5+vyrSaqT5Ch7\n6punyDE5VKkCn34Kdeo4z/n3v+GJJ0IYpFKERyJSKowxPxhj3jPGrDbGfAP0B/YDumqUCqpJk/zf\nvCtVgnnz7EwZgE/XfcqafWscdR697FFEJERRqvJGRPxaRVbtXcXHaz8G7AyaTz+1ibK3cePwW3dE\nqdIUU3iVUncAyAZq+5TXBvb4VwdPeaD6R4wxJwOdYIw5JSIrgCaFBTR27FgSEhIcZUOGDGGIb8eq\nKvemTIExY5xlUVEwfTp07GiPs3Oyeeyrxxx1Wp7Tkj7N+4QoSlVeDWgxgBY1W7D2wNq8snFLxtHv\ngn5ER0WTnGxncw0Y4OySufdeSEiwm+QpFci0adOYNm2aoywtLa1E1xITBh2CIvIDsMwYc6/nWIBt\nwERjzD8D1H8W6GGMaedV9gFQzRjTM5/niAJ+AeYZYx7Mp04SsHz58uUkJSUFqqJUng8+gGHD/PvU\nX3sN7vBaFuS91e9x839udtR5v//7DG0zNARRqvLuo18+YuDMgY4y39+/SZMCJ9TTpsFA56lK5Ssl\nJYXk5GSAZGNMSlHPC5eumfHAKBEZLiIXAK8DlYC3AUTkGRF5x6v+60AjEfl/ItJcRO4CbvBcB885\nfxORriJyvoh0AN4H6gNvhuZHUmXZhx/a2TC+ScjjjzuTkMzsTP7+1d8dddrUaqPLuauQGdByAG1r\nt3WU/WPJPziVc3p997vugqedQ5jIybGza959NxRRqvIsLBIRY8wM4EHgcWAF0Bbo5jXdNhGo51V/\nC3AdcA2wEjtt9zZjjPdMmurAZOBXYB5QBejkmR6sVIn96192n47sbGf5I4/AY84eGN5MeZPU31Md\nZU9f/TRREhZ/eqociJIonrjSOYhpw6ENvLH8DUfZX/4CDz3kPDcnx3bPvOGsqlRQhUXXTLjQrhlV\nmOeeg4cf9i+/7z67UJT32NO0jDSavdKMfen78so61+vM0pFLdZCqCiljDBe9eRE/7fopr6xmpZps\n+NMGqsVV86pnx4e8/LL/NV54Ae6/PxTRqkgV6V0zSoW17Gx44IHAScg99/gnIQCPf/24IwkBeObq\nZzQJUSEnIjzX9TlH2YHjB3jyv0/61IOXXoIHA4yie+ABm6T4tgQqdaY0EVGqEMeOQf/+/ktjA/z9\n7zBhgn8S8tuB35j4o3MOZJ/mfbi8weWlGKlS+bui4RX0b9HfUTZx2UQ2HHRuNCNiW/4C7cw7caKd\nYXNcd+1SQaSJiFIF2LEDLrss8A6lL7wA//d//kmIMYZ7F9zrGAxYIboC468NkMkoFULPXfMcFaIr\n5B1n5WRx74J78e2iF7EDr5991vcKMHs2XH45bNtW2tGq8kITEaXysXgxJCfDypXO8thYeOut/PvL\nZ62dxaJNixxlD3R6gMY1GpdSpEoVTeMajbnvovscZfM3zmfGLzMC1n/4YbvOiO+OvcuX2+0LvvDd\naEOpEtBERCkfOTl2pdSuXWGfc4gHNWrA55/DLbcEPvfQiUOM+cy5IEPds+ry18v+WjrBKlVMj17+\nKHWqONd2v2fBPRw6cShg/aFDYdEiqFbNWX7wIHTrZltNcnJKK1pVHmgiopSX3buhZ0879sN3Qlmz\nZvDDD9ClS/7n37/wfr8Bqi9c+wJVKlQphWiVKr6qFavycg/ntJh96ft46POH8jkDrrgCvvsOmjZ1\nlufk2Gnr114LO3eWQrCqXNBERCmP6dOhdWtYuND/seuug++/938j9rZg4wLeWfWOo6x3s94MajUo\nyJEqdWb6t+jvt8XAv1f8m/kb5ud7TosW8NNP0Lev/2OLF0ObNjBzZrAjVeWBJiKq3Dt4EAYPtrdD\nPq3TUVHwzDN2sGqNGvlfY3/6fkbOHukoq1qxKq9d95pO11VhR0R4teerfi11I2eP9GvR85aQYHeU\nfuop/0Hahw/DjTfCTTf5d2kqVRBNRFS5lZNjN61r3ty2hviqUwe+/NKuOBlVwF+KMYZbZt/CnmPO\nPRqf7/o851Y9N8hRKxUc51U9j392dW7ltTd9L7fNuc1vFo23qCj4619tK8h55/k//sEHtvVkrndS\nFwAAFWhJREFUyhT/7k2lAtFERJVLq1fbKYi33WZbRHwNHAg//1zweJBcE5dN5LMNnznKujfpzu1J\ntwcpWqVKx+jk0VzX9DpH2dz1c3nlx1cKPffKK+3fUaBN8Q4dsn9buXWUKogmIqpc2b3bbkqXlATf\nfuv/eI0adkO76dPh7LMLv963277lz5//2VFWu3Jt3u77tnbJqLAnIkzpO4ValWs5yu9fdD9Lty0t\n9Pzq1e3fy9SpttvG19dfQ4cOMGoU7Nnj/7hSoImIKieOHLErRTZpYjetC7RM9YABthVkUBHHlm5P\n286AGQPIyslylE/tN5XaVWoHIWqlSl+tyrV4q+9bjrJTOae4YcYN7DxS+FQYEbj5Zvjtt8B/Ozk5\n8OabdqD3E0/Yv0WlvGkiosq0tDQ72LRxY3jyycBLUzdqBJ99Zkf8161btOsezzpOv+n92Ju+11H+\n8CUPc23ja4MQuVKh07NpTx659BFH2d70vfSb3o/0zPQiXSMx0baOfPYZNGjg//ixY3ZafMOGdrCr\nJiQqlyYiqkw6cAAee8y+If71r/bYV3y8fWNcswZ69Cj6tbOysxg0cxDLdy93lHdr3I2nrnrqDCNX\nyh1PXPkE3Rp3c5T9tOsnBs0c5NiuoDA9esCvv9ol4itV8n/88GH7t9mwoa2zf/8ZBq4iniYiqkxZ\nvRpGj7YJyFNP2RYRX1FRcPvtsGGD3SsmPr7o188xOdz+6e3MXT/XUd60RlOmDZhGdFT0Gf4ESrkj\nOiqaDwZ8QKPqjRzl8zbMY/SnowucSeOrUiXbFbphA9x6q/9UX7AJybhxUL++HUOyZs2Z/gQqUmki\noiJeVhZ89JGd4dKuHUyenP/uoH362HEgb7wB5xZzZq0xhrELxjJ11VRHeULFBGYPnk31+Ool/AmU\nCg814mswd8hcqsc5f5enrJwScHO8wtStC//+N6xYAf36Ba6TkWHHkLRpY7dVmDEDTp4s6U+gIpEm\nIioiGQMpKXDffTahGDgQ/vvfwHVF4IYb7EZds2dDy5bFf74ck8Od8+5k4o8THeVxMXHMHTqXFue0\nKMFPoVT4aXFOC+YOnUtcTJyj/OUfX+aueXeRY4q/sUy7dnYhtBUroH///Ot98YUd8HruufZv++ef\ni/1UKgJpIqIiSmoq/POf0Lat3Rn3pZfy72OOjobhw+GXX2yLSVJSyZ4zKzuLkbNH8q/l/3JeX6KZ\neeNMLq1/ackurFSY6lyvM9NvmE6UOD8iXl/+OiNnjyQzO7NE123fHmbNglWrYMQI/119cx08aP+2\n27a103+feQY2bizRU6oIIMVtaivLRCQJWL58+XKSSvqppYLKGDvw7T//Of2NqjC1a9txIn/8Y/G7\nX3wdPnGYGz66gS9Tv3SUR0s07/Z7lyFthpzZEygVxqavmc5NH99EtnHOd7/q/KuYNXAW1eKq5XNm\n0ezZA6+9Zm9FGbSalGSXke/bFy64IPDYE+WelJQUkpOTAZKNMSlFPU8TES+aiISH9HS7ENKiRTB/\nPqxfX7TzLrkExoyx64FUqHDmcaw7sI4+H/Zh/UFnALFRsXx4w4f0b1FAG7NSZcSsX2cxeNZgv5kz\nLWq24ONBH3NBzQvO+DkyMmyr5ZQpsGRJ0c45/3y7GWXPnnZ34OIMOlelQxORINBExB1ZWbapdvFi\nm3wsXQqZRWz5rV/fdr8MH17wzrjFNXXVVO6adxfpWc41FOJi4pg1cBY9m/YM3pMpFebmrp/LwI8G\ncuLUCUd55djKvN7rdYa1HRa059q4Ed56C95+G3btKto58fF2y4YuXWxS0rFj/t0+qvRoIhIEmoiE\nxtGj8MMPNuFYutTez2+WSyBnnw3XXw/Dhtk3n4I2pCuuQycOce+Ce3lv9Xt+jyVWSWTO4DlceO6F\nwXtCpSLE/3b9j97Tevtt7ghwc9ubmdB9AjXiC9iiupiys+Gbb+wsmpkzi7feSOXKtoW0Sxe4+GKb\nmFStGrTQVD40EQkCTUSCLz3dtnYsX25nuaSk2PUCcoo58P7cc+1o+/794dJLISYmuHEaY5i2Zhr3\nLbiP/cf93/Ha1W7Hp0M+pV5CveA+sVIRZOvvW+nzYR9W7/Xfya5W5VpM7D6Rga0GBn2fpVOn7Ky4\nGTNg7lzYWfjK8w4idkzJH/5w+ta6NcTFFX6uKjpNRIJAE5GSy86GzZvtwNK1a22ykZJi958oya9Y\nTAx06gTXXgvdutkZMsFs+fD2484fefiLh1myZUnAx0cljWJC9wlUig2wTKRS5cyJrBPcu+Be3kh5\nI+DjXRp04Z9d/1lqLYfG2IUL582zy8l//33xv9iAnVXXtKmdmdO2rV3HpG1buxiiDoItGU1EgkAT\nkYIZY6fVbd5sb+vXn0481q0780WImjeHq6+2yceVV5Z+U+qqPat48psnmfnrzICPV61Ylcm9JjOo\ndRF3wVOqHJn28zTumHcHR04G3jTmxpY38silj9ChTodSjePgQTvANfd2piu0VqkCzZrZJKVZM+f9\n6rpmYYE0EQmC8p6IGAOHDtlmzx07YOtWm3Bs2nQ6+Th6NDjPFRtrWzkuvdTeOneGc84JzrULkmNy\nWLhxIS98/wKLUxfnW69v87683ONl7YpRqgA7j+xkzGdjmL1udr51ujbqyv2d7qdro64h2QLhwAHb\njbNkiR1/tnKlHRAfDGefbTfJrF/ftpzk3nKPq1cv360pmogEQVlNRHJy7L4O+/bZAV/79tlkIzfh\nyP131y47jS7YoqLsaqbJyXYdgNxboA2xSsv6g+t5Z+U7vLv6XbYf2Z5vvXpV6/FS95fo1yKf9aiV\nUg7GGGatncUDix5gW9q2fOvVq1qPEe1GMLzdcJqeHcQpboU4edKOU/vxR1i2zN42bCid56pSBerV\ngzp17C0xMfD9hISymbBoIhIEkZCInDxpk4rff3feDh+2SUbuzTvpOHjQjuEIhZo1bdLRsqXtc01O\ntv+GMukAyM7JZtnOZXy67lPmbpjLmn0Ft9dWj6vOo5c9ypg/jPFb2lopVbiMUxm88uMrPPXNU/ye\n8XuBdVud04q+zfvSp3kfkusmExMV5NHnhTh61K64/PPPdrzJ6tX2/uHDoXn+uDj7Xlmzpm1l8b75\nltWsCdWq2a7qYA/SD7aIT0REZAzwIJAIrAL+ZIz5qYD6VwAvAK2AbcBTxph3fOrcCDwONATWA38x\nxswv4JpBT0SMscnDsWP539LT/cvS0pxJRu790mixKK7oaNsM2agRtGgB6enTGDFiCC1ahKZ7JZC0\njDRW7FnBt9u+Zen2pXy3/bt8+6691apci7svvJs/XfSnM14lsrimTZvGkCG6Mmtx6GtWMqF83Q6f\nOMyknyYx8ceJ7EvfV2j9syqcxSX1L6FLgy5cfN7FtE9sH/K/RbDv1Tt32vFu69fDnDnTiI4ewvr1\ndmuJU6cKv0Zpq1TJtqYkJNjEpKD7lSvnf6tUyd6C3SoT0YmIiAwC3gH+CPwIjAVuBJoZYw4EqN8Q\nWANMAv4NXANMAHoaYz731OkMfA08DMwDbvLc72CM+TWfOJKA5c88s5w6dZLIyLAf/CdO4Hc/UJnv\n47kJRqhaI4JFxC6Tft55NtnwvdWr58zM+/Tpw5w5c0o9rlM5p9h1dBdbft/Clt+38NuB31i9dzU/\n7/u5wCbhQNrWbsuf/vAnhrUd5loLSKhet7JEX7OSceN1O5F1gqmrpjI5ZTIpu4v8mQRAw2oNaZ/Y\nnhY1W9C4emOa1GhC4xqNqXtWXb/9b0qL92uWlQVbttguna1bYds2+2/ubdeuks0OdFulSoETlfh4\n22pT3Nv27SmMGVP8RCRcGnrGAv8yxkwFEJE7gOuAW4HnAtS/E9hsjHnIc7xORC71XOdzT9k9wHxj\nzHjP8d9FpCtwN3BXQcE88siZ/CjhrWpV22pRt65dm+Pcc23C4X2/Tp3SXZUwOyeb9Kx0jmUe87sd\nOnGI/en72X98PweOH2D/8f3sT9/PrqO72H5ku98y08VRq3IthrYeyoj2I2if2D6IP5FSyld8bDyj\nO45mdMfRrNyzkikrpjDjlxnsTd9b6Lm5XzY+4RNHeYXoCtSpUofEKokkVknMu18jvgYJcQkkVExw\n/Fu1YlXiYuKIj4k/o4GysbF25kx+qzdnZtrWlK1b7b+7d9t9dHbvdt7/veAeq5A7ftzeirNYXGlw\nPRERkVggGXg6t8wYY0TkC6BTPqddDHzhU7YQeNHruBO268a3Tt9Cg7r0aahW06vAK9UVU0CZV3mg\nshLXDXx+dIyhQizEVoDYWEPFilAxzv4bVxEqVDRUjPPcr2DLJer0+ZkYUoFUwOQY2A5sB+N5Pu/W\nMuMVQ265wZCdk81PO37iirevIDM7k6ycLLKys/L9N+NUht8y0aWp1Tmt6N2sN72b9+aicy8Kyah9\npZRT+8T2TOwxkQndJ7BsxzJmr5vNnHVzWHtgbbGuk5mdyda0rWxN21rsGGKiYoiLicu7xcfEExcT\nR0xUDDFRMURHRRMt0URHRRMTFcPqHavp8X4PR1nufcH2aYiI//1KIE0EmkAsQgMRGgCCcCobTp4Q\nTpwQMjLg5Enh5Ek4mSFkZAgnM/D8K55WdinRGimuSdsPS4t/muuJCFATiAZ80+S9QPN8zknMp35V\nEalojDlZQJ3EQiNqOQvqFlrLddnACc/NzynPLT3Qg0GWDl9v/ToET1SwaImmXWI7Lq13KZfWv5RL\n6l9C3bMi4H+kUuVElETRqV4nOtXrxLPXPMveY3v579b/8vXWr/lu+3es2beGrJwgzbX1cSrnVF7L\na5Gkw4KNC0olFj/RQGXPLZLtImITkXBiBwv4jUpRBcrA/gKGSFxMHHXOqsO5Z51Lo+qNaHp2U5pU\nb0LDag2pEOPZdvck7Nmwhz3474sRLtLS0khJKV7feXmnr1nJhPPr1pjGNE5szK2Jt5KVncWW37ew\n7uA6NhzcwPYj29lxZAc7juzg5KkzXDGxuEL8vlYmnP7sLNbAu3BIRA5gv9zX9imvDfl+iuzJp/4R\nT2tIQXUK+mRqCMDHBcarApkcuqfKIINUz39LS5J+hxHPCHNVDPqalYy+biUQwve1MqYh8F1RK7ue\niBhjskRkOXA1MAdA7I5JVwMT8znte6CHT9m1nnLvOr7X6OpTx9dC7OyaLdh8WCmllFJFE4dNQhYW\n56Rwmb47EHgbuIPT03dvAC4wxuwXkWeAusaYEZ76DYGfsdN3p2ATjtzpu1946nQClgCPYKfvDgH+\nAiTlN31XKaWUUqHleosIgDFmhojUxC4+VhtYCXQzxuROKkoE6nnV3yIi12FnydwD7ABuy01CPHW+\nF5GhwFOe2wagryYhSimlVPgIixYRpZRSSpVPoVmiTimllFIqAE1ElFJKKeUaTUQKICLXicgPInJc\nRA6JiE7sLQIRqSAiK0UkR0Tauh1POBORBiLypohs9vyebRCRf3hWHFYeIjJGRFJF5ITnb/JCt2MK\nZyLyiIj8KCJHRGSviPxHRJq5HVckEZG/eN7Dxhdeu3wTkboi8q6IHPC8j63y7N1WJJqI5ENEBgBT\nsZvqtQE6Ax+4GlTkeA47gFgHIBXuAkCAUUBL7IyxO7ADrBV5m2K+AIwDOmB3517oGeCuArsMeBm4\nCLspaCywSETiXY0qQngS3T9if9dUAUSkGvAtcBLoBrQAHgAOF/kaOljVn4hEY9cS+Zsx5m13o4ks\nItIDeB4YAPwKtDfGrHY3qsgiIg8CdxhjmrgdSzgQkR+AZcaYez3Hgt0ZaaIxJtCmmMqHJ2nbB1xu\njInsVQBLmYhUAZZjN1f9G7DCGHO/u1GFLxF5FuhkjOlS0mtoi0hgSXh2mxGRFBHZJSKfiUgrl+MK\nayJSG7sW4TDy2QJHFUk14JDbQYQDr00xF+eWGfvtqaBNMZW/atgWSv29KtyrwKfGmC/dDiRC9Ab+\nJyIzPN2AKSJye3EuoIlIYI2wzeXjsGubXIdtZlriaYZSgb0FTDLGrHA7kEglIk2Au4HX3Y4lTBS0\nKWbhG1iq3BakCcBSXUepYCIyGGiPXQhTFU0jbOvROuwK568BE0Xk5qJeoFwlIiLyjGfwUX63bM+A\nrtzX5UljzCeeD9aR2G8UN7r2A7igqK+ZiNwDVAH+X+6pLobtumL8rnmfcy4wH5hujJniTuSqDJqE\nHX802O1AwpmInIdN2G4yxpTOFsBlUxSw3BjzN2PMKmPMG8Ab2LFuRRIWK6uG0PPYb+0F2YynWwZY\nm1tojMkUkc1A/VKKLVwV5TVLBa7ENpWftF/A8vxPRN43xowspfjCVVF/1wA76hz4EvutdXRpBhZh\nSrIppvIQkVeAnsBlxpjdbscT5pKBc4AUOf0mFg1cLiJ3AxWNDqoMZDden5Uea4H+Rb1AuUpEjDEH\ngYOF1fNswncSaI5nB0FPX3VDYGsphhh2ivGa/Ql41KuoLnbjo4HY/YPKlaK+bpDXEvIl8BNwa2nG\nFWlKuCmmIi8J6Qt0McZsczueCPAFdoakt7exH6rPahKSr2+xn5XemlOMz8pylYgUlTHmqIi8Dvyf\niOzAvqAPYbtmPnI1uDBljNnhfSwi6djumc3GmF3uRBX+PC0hS7CtSg8BtXK/jBljfMdFlFfjgbc9\nCUnuppiVsB8SKgARmYTd6LMPkO4ZSA6QZozRncUDMMakY2f65fG8jx00xvh+41envQh8KyKPADOw\nU8Zvxy5JUCSaiOTvQSALu5ZIPLAMuMoYk+ZqVJFFv0EUrit2sFcj7JRUsAmcwTYLl3tF2BRT+bsD\n+zu0xKd8JPY9TRWNvocVwhjzPxHpBzyLne6cCtxrjPmwqNfQdUSUUkop5ZpyNWtGKaWUUuFFExGl\nlFJKuUYTEaWUUkq5RhMRpZRSSrlGExGllFJKuUYTEaWUUkq5RhMRpZRSSrlGExGllFJKuUYTEaWU\nUkq5RhMRpZRSSrlGExGlVMQQkaYi0sntOJRSwaOJiFIqkjyE/5bjSqkIppveKaUihoikApcYY3a5\nHYtSKji0RUQpFfZEpJeIvAoIcIuIdHY7JqVUcMS4HYBSShXGGDNXRGrZu+Zpt+NRSgWPtogopSLF\n5cBXbgehlAouHSOilIoIIrIJ+ANwGKhmjDnkckhKqSDQFhGlVNgTkbOBk8aYg8BwINblkJRSQaKJ\niFIqEhwGVonILcA2Y8xel+NRSgWJds0opZRSyjXaIqKUUkop12giopRSSinXaCKilFJKKddoIqKU\nUkop12giopRSSinXaCKilFJKKddoIqKUUkop12giopRSSinXaCKilFJKKddoIqKUUkop12giopRS\nSinXaCKilFJKKdf8fyCF7TibTfF9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10febdda0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(-6,6,0.001)\n",
    "plt.figure()\n",
    "plt.plot(x,stats.t.pdf(x,3),lw=3)\n",
    "plt.plot(x,stats.t.pdf(x,300),lw=3)\n",
    "plt.legend(['3','300'])\n",
    "plt.xlabel('$t$')\n",
    "plt.ylabel('PDF')\n",
    "plt.title('t distribution for 3 and 300 degrees of freedom')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The t distribution has a narrower distribution for many degrees of freedom. This is because with more samples, the sample mean gets closer to the true mean. There will be less spread in $\\bar{x}-\\mu$, and therefore also less spread in the value of $t$.\n",
    "\n",
    "We can define $t_{1-\\alpha/2, \\nu}$ as the value of $t$ where a fraction $1-\\alpha/2$ of the entire distribution is less than $t_{1-\\alpha/2, \\nu}$. Again, $\\nu$ represents the degrees of freedom in the mean calculation, which is equal to the sample size $N$. This parameter determines the shape of the distribution, as shown in the figure above. For $\\alpha$ = 0.05, 95% of the distribution falls between $t_{0.975,\\nu}$ and -$t_{0.975,\\nu}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Values for $t_{1-\\alpha/2, \\nu}$ can be obtained in statistical tables for different values of $\\alpha$ and $\\nu$. The `scipy.stats` library in Python can also be used. For example, for $\\alpha$ = 0.05 and $\\nu$ = 20, the value of $t_{1-\\alpha/2, \\nu}$ is 2.086. This can be verified in Appendix D.3 in Emery and Thomson."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0859634472658364"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "stats.t.ppf(0.975,20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `t` in this command refers to the $t$ distribution and `ppf` is short for \"percent point function.\" This is the inverse of the cumulative distribution function, which can also be found using Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.97499997769549873"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.t.cdf(2.085963,20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Confidence intervals for the mean\n",
    "\n",
    "If many groups of $N$ samples are taken from the same normally distributed parent population, a fraction 1-$\\alpha$ of the sample means will fall in within the confidence intervals\n",
    "\n",
    "$ \\bar{x} - t_{1-\\alpha/2} \\frac{S} {\\sqrt{N}} < \\mu <  \\bar{x} + t_{1-\\alpha/2} {\\frac{S} {\\sqrt{N}}} $.\n",
    "\n",
    "For a given group of $N$ samples and $\\alpha$ = 0.05, this interval means that there is a 95% chance that the true mean falls within these values, assuming that the parent population is normally distributed.\n",
    "\n",
    "This interval does _not_ describe the spread of the data, or where 95% of the samples fall. The spread is determined by the standard deviation. \n",
    "\n",
    "For $\\alpha$ = 0.05, these values represent the 95% confidence intervals. The use of 95% confidence intervals is common, but there is no reason that this should necessarily be a cut-off. We can never be 100% certain that the true mean falls within a certain interval, but we want to be pretty certain.\n",
    "\n",
    "#### Exercise\n",
    "\n",
    "You are calibrating a nitrate sensor with bottle samples. The sensor has less accuracy than the bottle samples, but you are only able to collect a handful of discrete samples. You obtain differences $\\Delta$NO$_3^-$ of 3.25, 0.43, -0.17 and 3.89 $\\mu$M. Calculate the sample mean, and the 95% confidence intervals for the mean, assuming that the instrument error is normally distributed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}