{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "uoP34MT3umZC"
   },
   "source": [
    "# Stochastic simulation helps you grasp concepts of statistics\n",
    "## Dr. Tirthajyoti Sarkar ([LinkedIn](https://www.linkedin.com/in/tirthajyoti-sarkar-2127aa7/), [Github](https://github.com/tirthajyoti)), Fremont, CA, July 2020\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation helps distilling concepts\n",
    "### Grasping statistic concepts can be hard\n",
    "Do you find grasping the concepts of statistical analysis - law of large numbers, expectation value, confidence interval, p-value - somewhat difficult and troublesome?\n",
    "\n",
    "You are not alone.\n",
    "\n",
    "Our human brain and psyche have not evolved to deal with rigorous statistical methods. In fact, [a study](https://www.sciencedaily.com/releases/2018/10/181012082713.htm) of why people struggle to solve statistical problems reveals a preference for complicated rather than simpler, more intuitive solutions - which often leads to failure in solving the problem altogether.\n",
    "\n",
    "We are good with a small set of numbers. The short-term working memory of the human brain is around [7–8 items/numbers](https://human-memory.net/short-term-working-memory/).\n",
    "\n",
    "Therefore, whenever a process presents itself with a scale of thousands or millions, we tend to lose our grasp on the 'inherent nature' of that process. The laws and patterns, which are only manifested at the limit of large numbers, seem random and meaningless to us.\n",
    "\n",
    "Statistics deals with large numbers and almost all theories and results in the statistical modeling and analysis are valid at the limit of large numbers only.\n",
    "\n",
    "![frustrated-at-stat](https://raw.githubusercontent.com/tirthajyoti/Stats-Maths-with-Python/master/images/Frustrated-at-stat.png)\n",
    "\n",
    "### Data science/Machine learning is rooted in statistics - what to do?\n",
    "In this era of data science and machine learning, where knowledge of the core statistical concepts are considered essential for success in those fields, this can be worrisome for data science practitioners and folks who are on their path to learn the trade.\n",
    "\n",
    "But do not despair. There is a surprisingly easy way to tackle this. And it is called 'simulation'. In particular, discrete, stochastic, event-based simulation.\n",
    "\n",
    "## Let me show you the simplest possible example\n",
    "\n",
    "![dice](https://cdn.pixabay.com/photo/2016/09/08/18/45/cube-1655118_1280.jpg)\n",
    "\n",
    "Suppose we are throwing a (fair) dice with 6 possible faces - 1 to 6. This event of the dice face taking up a value from the set {1,2,3,4,5,6} is represented by a random variable. In a formal setting, the so-called 'expectation value' (denoted by $\\text{E}[X]$) of any  random variable $X$ is given by,\n",
    "\n",
    "where $f(x)$ is the probability distribution function (PDF) or probability mass function (PMF) for $X$ i.e. the mathematical function that describes the distribution of the possible values that $X$ can assume.\n",
    "\n",
    "For a dice throwing situation, the random variable $X$ is of discrete nature i.e. it can assume discrete values only, so it has a PMF (and not a PDF). And it is a very simple PMF,\n",
    "\n",
    "$$f(x) = \\frac{1}{6}$$\n",
    "\n",
    "This is because the random variable has a 'uniform probability distribution' over the sample space {1,2,3,4,5,6} i.e. any dice throw can result in any one of these values, completely randomly, and without any bias towards any particular value. Therefore, the expected value is,\n",
    "\n",
    "$$\\text{E}[X] = \\sum_{x}x.f(x) = \\frac{1}{6}.(1+2+3+4+5+6)=\\frac{21}{6}=3.5$$\n",
    "\n",
    "So, as per theory, this is the expected value of the dice throwing process.\n",
    "\n",
    "**Is it the most probable value?** No. Because a dice does not even have a face with 3.5! So, what's the meaning of this quantity?\n",
    "\n",
    "**Is it some kind of probability**? No. Because the value is clearly greater than 1 and probability values are always between 0 and 1.\n",
    "\n",
    "**Does it mean we can expect the face to turn up either 3 or 4 most times (3.5 is the average of 3 and 4)?** No. Because the PMF tells us that all the faces are equally likely to turn up.\n",
    "\n",
    "Fortunately, the answer is provided by a fundamental tenet of statistics - The Law of large numbers - which says that, in the long run, the expected value is simply the average of all the values that the random variable will take.\n",
    "\n",
    "Notice the phrase \"_in the long run_\". How do we verify this? Can we simulate such a scenario?\n",
    "\n",
    "Sure we can. Simple Python code can help us simulate the scenario and verify the Law of Large Numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "zvw2bqmQumZD"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as stats\n",
    "from scipy import mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The `dice` array and a simple dice-throwing function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "dice = np.array([1,2,3,4,5,6])\n",
    "\n",
    "def dice_throw(dice):\n",
    "    \"\"\"\n",
    "    Simulates a single dice throw\n",
    "    \"\"\"\n",
    "    return np.random.choice(dice)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Throw the dice a few times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Here are 10 throws...\n",
      "3,3,5,1,3,2,3,1,2,6,\n",
      "\n",
      "Here are 10 more throws...\n",
      "3,6,4,2,3,3,5,2,6,3,"
     ]
    }
   ],
   "source": [
    "print(\"Here are 10 throws...\")\n",
    "for i in range(10):\n",
    "    print(dice_throw(dice),end=',')\n",
    "    \n",
    "print(\"\\n\\nHere are 10 more throws...\")\n",
    "for i in range(10):\n",
    "    print(dice_throw(dice),end=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate for a long time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "av = []\n",
    "n_throws = []\n",
    "for i in [5,10,15,20,25,50,75,100,150,200,250,500,750,1000]:\n",
    "    throws = []\n",
    "    for j in range(i):\n",
    "        throws.append(dice_throw(dice)) \n",
    "    mean = np.array(throws).mean()\n",
    "    av.append(mean)\n",
    "    n_throws.append(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAFbCAYAAAA0vux4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd5hU5fn/8fe9sAgIK9jQqAGCvSAK9gJYoliS6A+NMajEQhRjokGTGKPYjTFoYjCKX4wFa8QexRJhsRdsGLtr16hREVmaK9y/P54z7jDM7JydndmZs/N5Xde5ZufMKc88OzP3eeoxd0dERESqQ025EyAiIiLtR4FfRESkiijwi4iIVBEFfhERkSqiwC8iIlJFFPhFRESqiAK/xGZmo83MzWxYudMi7cPMupvZxWb2npktMbN3CjzOVWbm+dZVCjN7x8zqy52O9lAp3+tSpsPM+kXHPr3Yx04iBf4yMLNh0YcwfWk0s2fN7AQz61zuNIpEfgscB9wEjAaOL2tqisjMjjez0eVOhxRHFNxPN7NB5U5LpVOAKa8bgHsAA9YADgUuBDYCxpQxXblMAW4Evi53QqTd7A686O4nleDYRwFHl+C4cR0PvANcVcY0VIKO8r3uB4wn/E+fz3jtXaAb8E37JqkyKfCX17Pufm3qiZn9HXgVONLMTnH3/5Uvactz9yXAknKno1KZWU93n1fudBTZGsB7pTiwuzcBTaU4drmZWS3Qyd0XlTst+VTD99rDFLUV/79oL6rqryDuPh94glADMCC13szqs7WtZmu3SmtGGG1mPzOzl8xssZm9a2a/yXKMd6Ljb2hmd5vZPDOba2ZTzWyNjG2Xa4NLW7eLmZ1oZg3R+V43s8OynK+TmZ0apWeRmc02sx9HVXRuZv3y5ZOZfd/MbjKzt8xsoZl9aWb3m9nQjO1uMrOvzWzVLMfYIDrfXzLW/9jMHonyYYGZPWlmI7Ps71Eb9a7R9o3AXdFr3zGzCWb2vJnNid7ny2b2WzPrlOVY/czsFjP7Ksr7O8ysf652ZjPbLXq/X6blYeySs5l1jtLycrT/52Z2m5ltlrbNaAvt7/2BoWlNUqfnOXZXM7vAzD6K/jdPmdn3c2ybtY3fzNaw0K/greiz9KmZPWBmu2dst56ZTTGz/0b/53eic68YIw8c6Jvx3pb7/MX8XqQ+u5uY2YVm9gEhyGwbvb6qmV1iZu9H6Xw/er5K2jH6Zsvf6P/sZnZ8xvonzezlGO9zezObZmYfR//rD83sHjPbNm2blr7Xu5rZaRa+rwuj86be19Dosz8/+h+cmi2fzeyqLOtjteebWU8zOzs672fR5+FNM/ujmXVPPx4wI3p6Zdr/sz56PWsbf5zvQub+ZraPmT0dbf/f6DOXqEJ0ohJbJVIB/4s2HudooA9wBfAlMAo438w+cPfrM7ZdC6gHbgNOAjYHfg7UAVl/tLM4l1CVNglYDBwDXGVmb7r7o2nbTYzSNgP4M7Aa8Hfg7Va8t9HAysA1wAdR+o8EHjSz4e7+cLTd1cCBwEHRedMdmrYNAGZ2NnAKcC9wKrAU2A+42cx+4e6XZBxjCPD/gP9LPw4wENifkJ8NQC0wAvgj8D1C3qbOuQrwMOF/dRnwCrATIX+WC2BmNiba7gngHGA+oTr+UjMbELNK/jpCvjwAXEoo1R8LPG5mO7n7c8BDwCHARcBn0bkAZuc59g3AjwgXQfcRPs+3EvP/GwXeRwn5cQ0wi5AP2wK7RWnGzAYD0wmf7UnAh4TP7S+BHcxsaFSjkEu29waQXsvW2u/FdcBCYALgwH/NbCXgMWBd4B/As8AWhO/HLma2tbvPc/d3zextYFfg9Og9dgF2IHwOdwX+Eq2vAwYTPgc5mdkGhPz6GPgr8Anhf71D9F6eaGn/yB+BTtH+XYBxwH0WLuqvAC6n+fN0ppm9nV6LWQSp7/YtwPWEqvqhwG8I+bhHtN1DhN+g30dpSv0GfJLn+HG+C+n2AsYS8v4fwA+BE4E50fmTwd21tPMCDCP8MJwGrEoIfpsBl0Trn8rYvh54J8tx+kXbn57l2B8BvdLWdyf8qD2ecYx3ou0PzFifSsuGaetGR+uGZVn3HNAlbf1ahAuAG9LWbRJtey9Qk7Z+M0JVowP9YuTfilnW9SH8iN+Ttq4T8N8s+WmENr/Zaeu2jM5/bpZj3w58BfRMW+fRsluW7bsBlmX9lOh9rpm27k/RcX6asW1qfX3aujUJJcnrsxz7r9GxB+TJu92j496UnkbCxco3wMNZPh/1LR0zbdvvR8e+KmP9j1L5lbH+qizr7om23SPL8dM/My8QmsV6ZmyzX7T/6BjpzfneWvm9OD31vwI6Z2x/TvTa2Iz1x0brz0pbN5nQzr5i9HznaJsp0eevc7R+32j9/nne3y+j7bbOs91ocn+vn2XZ7/UPovXfAFulre9C+K5l/r4s93mIcc5hGcetzbL/WZnvjebfvuX+92T/rYz9XUjbfz5pv1GE35L/AP+N8x2plEVV/eV1BiEYf0ooSY0llI5+UIRjX+nuX6aeuPsCwhX+elm2/cjd/5mxbnr0uG7M8/3d3b/tHOTuHwKvZ5xvn+jxr+6+NG3bFwmlw1g8NIkAYGY9olLzEuBJYJu07ZYQrui3MrMN0w4xDPguy5bSf0r4Yl8dVc1+uwB3Aj2B7TKS8oK7/ztL+hZ66lfBrIuZrRwd5z5C89qQtM33Jfxg3pBxmD9neesjgRWAK7Kk8a7o2Ltm2S/dftHjOak0RmmeDfwL2NHMVstzjFx+FD1ekL7S3W8HXsu3s5mtDOwJ3Ovuy30eUp+ZqBp2IKEEuEJGPjxC+HGOW1PVktZ+L/7i7pmdx/YjfMcvz1g/iXChul/auumE2qEdo+e7EH4b/kr4/G0VrR9OqAWoz5P+udHjD82sa55tc7k0/XtNc0n6CXd/OrUy2uYpsv++FMzdv/ao5iaqlu8d/Z9T37ttcu+dVyHfhdvd/Z20bZ1QO7eGmfVoQ1ralQJ/eV1OuOrcizBs6gtgbYrTCeWtLOs+B1bJsj7XtuTYvtDz9Y8eswWBvIEhxcwGmNmNZjYHmEf4Af0fIR97Z2yeCu6Hpq07lHChcF3auo0IV++vRsdKX66ItumTcezXc6Svs5n9wcxeJ/wvP4+OMyXaJD2N/YE30y+EANz9U0I1drqNosd/Z0njAznSmKk/IWi8kuW1/6RtU4jvRcfOli/ZzpdpXcL/ILN6NVMqH1IXzunLp4SmgXz5EEdrvxfZ3nd/4LXMC4Lo+WuEPEtJXVTskvY4g1DqnpOx/gV3z9cceCPhs/J74Aszmx61Z/fNs1+6ZfLA3edEf2ZruplD/N+L2MxsrJnNJtQgfkH4P9dHL2d+31ujkO9CMX4ry05t/OX1RlqJcZqZPUIosVxGaJdOyTXJSUv/v9b00m1pW2vjMSzH3wWJrqofIvy4/wV4kRD8lwIn0/zjCITaBDN7HhhlZqcQquH/H3C/u3+ckTYntMXnei8vZTxfkGO7C2ke+34OIRg1EZoTzqfwC+5U/h1KqCXIJtsPU7ZjlEJLx45z3tQ2+Sb1SW03gdBslM2cHOtbo7Xfi1yfh1jc/WMze4XQ9t+dUJo9zt2XmtlMYFczu4xQ23FhjOMtBnY3s60JbeE7A2cCp5vZwe5+W4xk5cqDto4CiBV7zOzXhP/z/cDFhCbMrwlNiVfRtsJrId+FYvxWlp0CfwVx98fMbApwqJld7O6PRS99QejMk+l7WdZVslQpYQOWD1AbxDzGrsB3gMPd/cr0F6LOedlcTejINZzQTt6TZav5Ad4gVDO/5+5xSqctOQR4yN3TL94ws2zVw+8A65pZTXqp38xWB3plSSPAZ9maGGJqIASBjVi+o97G0WNrOlpmHvv7wPosf5G04fKbL+cNQtDfIsZ2AEvakA+Q/wKjGN4CNjCzzuml/qgX+Pos/z2YTuj4ty+hffvBaP2DhOafEYQAM52Y3P0pQjU8ZrYOoUblbEKnxVL7gtARN1Pc365DCN+RERnfjz2zbNva/2cpvwsVTVX9lecswlXlmWnrXgd6RlfuAJhZDXBCO6etre6KHn8VpR/4ts12j+y7LCd1xb3M1bWFIWO52vtSvYEPjZa5wB0Z26Sq4c+17EPuVo+ZvlQaM9O3Itn/X3cRLkZ+krH+xCzb/pNQ3XmGmXXLksaVzGyFPGm7PXo82cy+TaOZbUroW/KIFz5/RCpPlxlZYGY/IsaFXVR1PQ0YYWa7Zb6elt7nCFWxR5vZcgEkamrJFmwyNZI9KBXT7YTOu0dmrD8qWp8ZfKcTfpfHEy5CG9LWr0Co1fqG5rb2nCzLMFbCKJj/Ufr3nfI6sF3G0LvewM9i7p/q9Jv+We0M/C7Lto3RY9z3VsrvQkVTib/CuPubZnYj8NNoOMnDhL4A44DbzOyvhKqukSTs/+fuL5nZ5YRZCf9tZrcRfvyOJfyYDyb/VfsjhOFJE6KhXx8AgwglgxcJIwQyz/upmU0j5FlX4ArPmFjF3Z82s/GEduPnzexmQrXimlG69iKUwOKYCvzczG4itLH2AQ6nuS0w3fnAwYSxx1sT+hjsSBhy9Rlp+eHuH5jZMYTe369EtUPv0jwq5EeEkso7uRLm7g+Y2T8JTUm9zexfNA9hWkToCV4Qd7/PzO4CDosC772E4Xw/JwTqTWMc5heE4W/TzOxq4BlC88w20fv6rbu7mR1CCIazzewfhBqG7oR+AvsTAuRVec71BHCEmZ1FaOddCtyV3nm0CP4EHABcYmZbEj7nWwBHENr4/5Sx/YwoHRuRln53f9nMPib8fx/3eBNF/SG6IP4XoeRqhJqEDbOct1QmAtcC06PPay/CRc+7hM9dPlOB8wifh1sJQykPJvvETy8Tmv3GmtkCQh+ZT909a+1IKb8LFa/cwwqqcaF52MmJOV7fiHClOyNt3V6EaSgXEwLS+YRSVK7hfKOzHPcqlh8+9Q5ZhjRlOw4xh+CkvVZPxjBEwhC78YTZ4BYTqtgOJFRjOrB6jPwbSAgqqc599YSx78u9v7R9/h/NQ/B2aOHYexN6338Rpe99Qin0mIztsg5Til7rTujZ/i7hB+QNQgll12z/G0IHoluj9/IVoeTcn4zhiWnb70AoKX5KuAj8iBAwxgFdY+RfZ0Jn0ldo7jB1O7BZlm2zfj5aOHY3Qpvsx4Qx7U8TanOyffay/r8I7beXRZ+Rrwljse8Hds3Yrm+03TvRdp8TLhTOA9aJkdbVCePDvyAEWycaqpXrfZP9e3F6+r5Z9knNVfEBIWB9QBgWuGqO7Z+JjndIxvrrovVnx/xfDCP0M3kn+l98QRj5ciTLDl8bTeu+17mG6OX6f55E+C4sjj5zh8c9J+H34mTgzWj/dwkXLRuR8dsXbb8XoTPkItKGw5JlOF9rvgu59o/z/6/ExaKEi5RVVFLcBajzMAyvqkVDFD8DJrl7OeezF5EORm380q5ytE0PJHRaml6NQT9bnhBKIdA8TE9EpChU4pd2ZWFO+UOBuwmdjDYktPnXEKrg843h7nAszCf+LmF62k6EJoF9CG3dO1fjxZCIlI4Cv7SrqAPbWYQOeSsT2rUfAc5w92fKmbZyMbNxhIuhfoQ28g8Ibf5neMe725+IlJkCv4iISBVRG7+IiEgVSdQ48EKtuuqq3q9fv6Idb/78+ay4Yt5bflc95VM8yqd4lE/5KY/iqYZ8euaZZz5z96w33KqKwN+vXz9mzZpVtOPV19czbNiwoh2vo1I+xaN8ikf5lJ/yKJ5qyCczezfXa6rqFxERqSIK/CIiIlVEgV9ERKSKKPCLiIhUEQV+ERGRKlLWwG9ma5lZo5m5mfVoYbutzOxKM3vTzBaY2WtmNt7MurZnegEaGuCii9ajrg5qaqCuDsaODetFREQqXblL/BcAjTG2+zHhvt7nE267eAnwa8JtKtvNtGkwcCDcffeazJsH7jBvHkyeHNZPm9aeqREREWm9so3jN7OdgD2BcwkXAC05393/l/a83swWAZPMrK+75xyvWCwNDTByJCxYAJnXS01NYRk5EmbPhgEDSp0aERGRwpSlxG9mnYC/AWcS7jneooygn5K6i9vqRUxaThMmhODekqYmuOii9kiNiIhIYcpV1X800JVQZV+o7YGlwGtFSVEe114bL/BPmdIeqRERESlMu9+dz8xWAd4ARrn7PWY2GrgS6Onucdr7MbM1gNnAPe4+Osc2Ywj3eadPnz6Db7zxxjale5ddhuJuMdLmTJ8+s03n6igaGxvp0SNnn02JKJ/iUT7lpzyKpxryafjw4c+4+5Bsr5Uj8F8G9HX3EdHz0bQi8JtZF+DfwNrAYHefk2+fIUOGeFvn6q+rCx354mw3d26bTtVhVMN82MWgfIpH+ZSf8iieasgnM8sZ+Nu1qt/MNgEOB84ws15m1gvoHr28kpl1y7O/AdcAmwB7xQn6xTJqFNTWtrxNbS0cckj7pEdERKQQ7d3Gvx5QCzwOzImWVDv/B4QOfy25CPgh8EN3f7VUicxm3Lh4gf+EE9onPSIiIoVo7+F8jwDDM9btCfyWMD7/rVw7mtnJwHHAge7+SMlSmMOAATB1ahiyt3jxUpYsab5mqq0Ny9SpGsonIiKVrV1L/O7+mbvXpy9AquT+sLu/BhDN0HdFaj8zO5gw3v8a4EMz2zZtWa290j9iRBinv88+H31b+q+rgzFjwvoRI9orJSIiIoUp2wQ+eXQGOqU9/370ODpa0v0MuKrkKYoMGADHH/8mvXuvzfTp8G7Jpw4SEREpnrIHfne/iozA7e79Mp6PZvmAX1Y1NbB0ablTISIi0jrlnqs/scwU+EVEJHkU+AtUUxNu0iMiIpIkCvwFUlW/iIgkkQJ/gRT4RUQkiRT4C6TALyIiSaTAXyAFfhERSSIF/gIp8IuISBIp8BdIgV9ERJJIgb9ACvwiIpJECvwFUuAXEZEkUuAvkAK/iIgkkQJ/gTRlr4iIJJECf4FU4hcRkSRS4C+Q5uoXEZEkUuAvkEr8IiKSRAr8BaqJck6lfhERSRIF/gKlAr9K/SIikiQK/AVS4BcRkSRS4C+QAr+IiCSRAn+BFPhFRCSJFPgLpMAvIiJJpMBfIAV+ERFJIgX+Ainwi4hIEinwF8gsPCrwi4hIkijwF0gT+IiISBIp8BdIVf0iIpJECvwFUuAXEZEkUuAvkAK/iIgkkQJ/gRT4RUQkiRT4C6TALyIiSaTAXyAFfhERSaKyBn4zW8vMGs3MzaxHnm1XMrMrzWyOmc01s+vMbJX2SmsmBX4REUmicpf4LwAaY257EzAMOBIYDWwF3F6SVMWQCvxvvw1jx0JdXVhXVxeeNzSUK2UiIiK5lS3wm9lOwJ7An2Nsux2wB3CYu9/i7rcBo4AdzWy30qY0u1Tg33tvmDwZ5s0Lk/nMmxeeDxwI06aVI2UiIiK5lSXwm1kn4G/AmcBnMXYZAXzi7g+lVrj7U8Db0Wvt7tNPw+PChdDUtOxrTU2wYAGMHKmSv4iIVJZylfiPBroCl8TcfkPg1SzrX4lea3dxSvNNTXDRRaVPi4iISFztHvijDnlnAb9296Z820d6A19mWT8neq3dPfpo/m2ammDKlNKnRUREJK7OZTjnOcCT7n5PK/fLdjscy7EeMxsDjAHo06cP9fX1rTxdbo2NjSxc6NHpWzZvnlNfP7No506SxsbGouZ7R6V8ikf5lJ/yKJ5qz6d2DfxmtglwOLCzmfWKVnePHlcysyXuvjDLrnOA1bKs70X2mgDc/XLgcoAhQ4b4sGHD2pL0ZdTX19Otm7EwW0oz9OxpFPPcSVJfX1+17701lE/xKJ/yUx7FU+351N5V/esBtcDjhGA+h+Z2/g8IHf6yeZXsbfm52v5Lbued829TWwuHHFL6tIiIiMTV3oH/EWB4xnJ+9NpehHH92UwD1jCzHVMrzGwI8L3otXb3ox/l36a2Fk44ofRpERERiatdq/rd/TOgPn2dmfWL/nzY3RujdW8CM939iGi/x83sPuAaMzsRWEq4YHjE3f/dPqlf1lprhceuXWHJkmWH9NXWhmXqVBgwoBypExERya7cM/fl0hnolLHuIGAm8A/gGuAZYL92Tte3UhP4XH89HH548/q6OhgzBmbPhhFlmWFAREQkt7IHfne/yt0tVdqP1vVz99EZ233p7j9z917uXufuB0c1CGWRCvzf+U7zWP1jjoG5c2HiRJX0RUSkMpU98CdV+k16liwJf6ceRUREKpUCf4HSA3/qDn26U5+IiFQ6Bf4CWTR3T3rgV4lfREQqnQJ/gVIlfndV9YuISHIo8BdIVf0iIpJECvwFyhb4VeIXEZFKp8BfIPXqFxGRJFLgL5Cq+kVEJIkU+Aukqn4REUmi2IHfzFYws2PM7Aozu9/M1ovW/9jMNipdEiuTqvpFRCSJYt2kx8zWBx4AViLMkT8M6Bm9vBOwN3BoCdJXsVTVLyIiSRS3xH8x8B7QD9gDsLTXZgI7ZtmnQ1NVv4iIJFHc2/LuBBzg7l+aWeZd8z4B1ixusipftqp+lfhFRKTSxS3xLwK65XhtLeDL4iQnOTRlr4iIJFHcwP8A8HszWyltnZvZCsBxwD1FT1mFS5+yV4FfRESSIm5V/0nAo8CbhIsAB04DNgG6APuXJHUVTFX9IiKSRLFK/O7+PrA5cBmhg18DoV3/ZmCwu39cqgRWKnXuExGRJIpb4sfd5wCnRkvVU+AXEZEk0sx9BVJVv4iIJFHcCXzeJrTrZ7MU+Ap4AZjo7s8UKW0VTSV+ERFJorgl/lsIFwk9gSeBf0WPdUAtMAvYFnjCzPYoQTorjmbuExGRJIrbxv8p8Dqwj7svSq00s27AXYRZ/TYF7gTOAO4rcjorjubqFxGRJIpb4v8lcGF60Adw94XARcCx7r4E+D9gs+ImsTKpql9ERJIobuDvBfTJ8VofoEf091ygKsKfqvpFRCSJ4gb+fwF/MrP9zawLgJl1MbORwJ+i1yGU9huKn8zKo6p+ERFJorht/EcDVwNTCVP1ziN09DNCG/8x0XYfAb8vdiIrUWqufk3ZKyIiSRIr8Lv7l8APzWwTYAiwBvAxMMvdX0rbbmpJUlmBVNUvIiJJFHvmPoAoyL+Ud8MqoKp+ERFJolYFfjNbG1gf6Jr5mrtX1R361KtfRESSKO7MfT2BfwLfT62KHtNn8+tUxHRVPFX1i4hIEsXt1X8e8F1gJ0LQ3w8YBlwBvE2Yta+qqKpfRESSKG7g3ws4hzBNL8BH7v6Qu48B7gBOKkXiKplK/CIikkRxA38f4P1odr75wMppr91DcxNAi8xspJk9Zmafm9kiM3vNzP6Qmhughf2GmNn90X5fmNm/zWybmGkvCbXxi4hIEsUN/O8Dq0Z/vwHsk/baNsCi5fbIbhVgBnAkMAL4B3AKcGGuHcxsHeDfhP4IhwKHRH/fb2Z9Y5636FTVLyIiSRS3V/8DwG7AbYS5+a82s8HAYmBnYEKcg7j7pIxVM8ysDjjWzI5z92y3/t2bMFnQ/tF8ApjZY8BnhCaIS2O+h6JSVb+IiCRR3MD/W6A7gLtPMbNGYCTQDfgFkBnQW+NzoKWq/lrgG6AxbV1jtM6y7tEOVNUvIiJJlLeq38xWIPTi751a5+63uftP3X1/d7/U3VtV1jWzTmbW3cx2JNz579IcpX2AW4AFwAQzW93MVifUOswBbm7NeYspNWWvqvpFRCRJ8gZ+d18MTAa+U8Tzzo+Wh4GZtDAqwN0/AoYD/w/4JFr2B/Zw9/8VMU2tkirxp8/Vr6p+ERGpdHGr+l8kzNg3s0jn3Z7QdLA1cBowERibbUMzW5Nwc6BnCJ0CAY4F7jaz7d39vRz7jQHGAPTp04f6+voiJR0aGxuZObMeGMZbb71DY+NiYAOWLHHq64uVRcnX2NhY1HzvqJRP8Sif8lMexVP1+eTueRdgB5p783eOs0/chdBT34EBOV6/EHgHqE1b1wV4F7g4zjkGDx7sxTRjxgx3d6+pcf/DH9wvu8w9lP2LeprES+WTtEz5FI/yKT/lUTzVkE+Em+hljYlxS/y3E0rodxBuyzuHZafrxd1XL/Da49nosT/QkOX1DYGX3L0p7Vxfm9lLwIACz1kUNTXLdu6D8HdN3EGSIiIi7Sxu4L+EjEBfRDtEj2/neP1dYC8z6+LuX8O3HQ43Be4qUZpiyRb4lyxR4BcRkcoVK/C7++nFOJmZ3UuYjOclYAkh6I8DbnL3hmibN4GZ7n5EtNtkQtv+bWb2d8IQvmOBNYHLi5GuQqUCf3pvfnXwExGRStba2/L2JpS01wGmufscM+sKfO3xhvQ9DYwG+hHG4b8FnAxclpGmb+/05+7PmNmewHhgSrT6RWB3d3+hNekvtlwlfhERkUoV97a8nYFzCSXtboRq/60IY+lvAWYRAnOL3P1U4NQ82/TLsu5B4ME4aW1PCvwiIpI0cVujzwGOIszS9z2WnTHvDmDfIqcrEVTVLyIiSRO3qv9Q4HfufqWZdcp4rYFwMVB1VOIXEZGkiVvi70X2oXYQxtRnXgxUBTMFfhERSZa4gf8/wA9zvDaC5rH4VaWmJkzbo6p+ERFJirhV/WcDt5hZN8KNcRwYZGb7AT8HflCi9FU0VfWLiEjSxCrxu/sdwMHAbsA0Que+yYSheYe4+32lSmAlyzVzn4iISKWKPY7f3f8J/NPM1gdWBb4AXovmBK5K2Xr1q8QvIiKVLO44/uFAfTT3/+vA66VNVjKoql9ERJImbue+B4GPzOxiM9u+lAlKElX1i4hI0sQN/JsB/wd8H3jEzN4zswvMbHDpklb5VNUvIiJJE7dz30vufpq7bwhsCVwH7Ac8bWZvmtnZpUxkpVJVv4iIJE2rbyDr7s+7+8nuvi5hGF83wo12qo6q+kVEJGladXc+ADNbGdgf+DEwFFgIXF/kdCWCqvpFRCRp4vbqryNU7f8Y2JVwS927gYOAu919cclSWMFU4hcRkaSJW+L/lDBb332ESXvudPf5pUpUUmiufhERSZq4gf9o4FZ3/6qUiUmabHP1K/CLiEglixX43f2qEqcjkVTVLyIiSRO7c5+Z9QNGAesDXTNfd/cDi5aqhNBwPhERSZq4nfsGAzOB9wmBfzawEqfNCcEAACAASURBVNAP+AB4s0Tpq2jq1S8iIkkTdxz/BcAtwKaEO/Md4e7fA3YkdPr7U2mSV9lU1S8iIkkTN/APIozVT4W1rgDu/hhwBvDH4iet8qmqX0REkiZu4Hfg6+gWvJ8CfdNeex9Yr9gJSwJV9YuISNLEDfwvAwOivx8HTjCz9cysL/AboKEUiat0quoXEZGkidur/3KaS/m/B+4HXo2ezwdGFjldiaCqfhERSZq44/inpP39ipltBGxHuEHPE+7+aYnSV9GyVfWrxC8iIpWs1TfpAXD3RuCBIqclcTRlr4iIJE2rb8srzVJT9i5dCp06hXUK/CIiUskU+Nsgvaq/tjasU1W/JFVDA4wdC3V14bNdVxeeN1Rl112RjkuBvw3SO/elAr9K/JJE06bBwIEweTLMmxdqsubNC88HDgyvi0jHoMDfBgr80hE0NMDIkbBgATQ1LftaU1NYP3KkSv4iHYUCfxuoql86ggkTlg/4mZqa4KKL2ic9IlJacW/S848WXl4KfAU8D9wa9fivCirxS0dw7bXxAv+UKTBxYvukSURKJ26JfzNgL2A0MAIYEj2OBvYB9gYmA6+Y2fq5DmJmI83sMTP73MwWmdlrZvYHM+uSLwFmtr+ZPW1mC6P97zWzFWOmvySyBX6V+CUpliyB554LbflxNFbNJb1IxxY38J8GfAls4+5ruvtAd18T2BaYC5wEbADMI9zJL5dVgBnAkYQLh38ApwAXtnRyMzuScJOgadF+RwJvUOA8BMWSrapfJX6pVEuWwDPPhKr9H/wAVlkFttwy/v49epQubSLSfuIGzj8B49396fSV7v6UmZ0OnO/uG5nZH4G/5jqIu0/KWDXDzOqAY83suOgmQMsws1WBi4Dj3P3/0l66LWbaSya9xN85ykkFfqkU33wTSvQzZ0J9PTz8MHz1VXhtvfXgwANh6FC47z648cb81f1Dh5Y8ySLSDuIG/nWBhTleWwD0i/5+F1ihlWn4HGipqv/A6PHqVh635NIDf5foHaiqX8rlm2/g2WdDkJ85MwT6VDX+BhvAQQfBsGEhgH/nO837bbst3HJLy4HfDO66C/bdF/7853A8EUmmuIH/OWC8mT3l7h+nVprZmsB44JloVV/go3wHM7NOhAuELYFfApdmK+1HtgFeA44ws1OAPsCzwAnu/ljM9JeEqvqlnJqaQtV9qkT/yCPN7fAbbQQ//WkI9DvvDGuumfs4AwbA1KlhyF5T07IXALW1YbnhBnj9dTjrLNh0U/jFL+C006B371K+QxEphbiB/2jgPuAdM3sG+B+wGqGT3+fAHtF23wH+L+sRljWf5pqBawh9BHJZg9B/4A+EWwB/Hj3ea2brufsnMd9D0Zk1T9mrwC+l1tQEs2Y1l+gfeQTmzw+vbbwxHHpoKM0PHQp9+rTu2CNGwOzZYcjelCnhAqJHDzjkEDjhhHBxAOEcp50GF18M11wDZ54JP/95c1OXiFQ+y13QztjQrBtwOCHYrwF8DDwNXOnuuZoBch1rS6A7sDWh4+D17j42x7YPALsBI9z93mhdHaFZYaK7n5pjvzHAGIA+ffoMvvHGG1uTxBY1NjbSo0cPzjhjY956a0VWWGEp3bsv4YUXenHMMW9y4IEfFO1cSZbKJ2lZrnxqajJefbUnL7zQixde6MV//rMSixaFm0L079/I5pvPZdCgLxk48Et6987TQF9kDQ0rcskl6/Lcc73p23c+Y8c2sPXWX5T0nPo85ac8iqca8mn48OHPuPuQrC+6e1kX4FDAgQE5Xr8per1rxvp/A7fEOcfgwYO9mGbMmOHu7gcd5L7++u6bb+6+227u4H7++UU9VaKl8klalsqnRYvcH3rI/ayz3Hfd1b1bt/CZAvfNNnM/7jj3W25x//TT8qY3ZelS9zvucF933ZDGESPcX365dOfT5yk/5VE81ZBPwCzPERNbXUFnZp3J0hnP3Re09liRZ6PH/kC2SUFfIQR+y0wKYfKgslGvfmmLxYvhySfh6qv7ctZZ8NhjsGhRaEIaOBDGjAnV9jvvHIbeVRqzMCxwzz3DxD5nngmbbRZu7DN+fGWmWURijuM3szozm2hmHwGLCOP1M5dC7RA9vp3j9X8RgvzwtPSsBAwGXmjDedtME/hIayxaFNrnzzgDhg+HXr1CYL/66n7MmQNHHw233w6ffQbPPw9/+Qvst1/lB9AuXeDXv4Y33gjt/ZdcEoYLXnxx/iGCItL+4pb4JxFm6JsMvAx8XcjJzOxeQhX9S8ASQtAfB9zk7g3RNm8CM939CAB3n2VmdwBXmNnvgM8InfuagEsKSUexqFe/tGThQnjiiRDs6+tD6X7x4vC5GTQolIyHDQN4lH333bGsaS2G1VYLQf+YY8KFwK9+BX//O1x4Yeg8aJl1diJSFnED/x6E4XOT23i+pwnT/PYDvgHeAk4GLstIU6eM/UYRZgS8kNAp8FFgF3ef08b0tEm2qn6V+KvXggXw+OPNw+uefBK+/jp8TrbcMgyBGzYMdtwxlPZT6uu/KVeSS2LTTcOkQPfcEy4A9t4b9tgjzBi4ySblTp2IxA3884E2d1X30AM/ay/8tG36ZVnXCBwTLRUjPfB36hSeV3uJv6Eh/MBfey00Ng6lRw8YNQrGjWseEtZRzJ8f2uVTgf6pp0LVdqdOIdD/6lch0O+wA6y0UrlT277MQsDffXe49FI4/XTYfPPQFHDGGbDqquVOoUj1ihv4JwBjzex+d1eZNpJe1a/AD9OmZU4CY8ybB5Mnw9VXh0liRowodyoL19gYAn1qHP1TT4XZ8jp1giFDQul26NAQ6Ovqyp3aytClS7gAGjUqBP9LL4Xrrgud/449tnnGSxFpP3ED/1rA5sBrZjaDcMOedO7uvy1qyhIgvcRfUxMCQLVW9Tc0hKC/IMvYjtSFwMiRYZKYpJT8582DRx9tLtHPmhUCfefOsNVWcOKJoUS//fbQs2e5U1vZVlkF/va30P4/bly4SLr00lA7tM8+av8XaU9x7843kjB0rjOwO3BAlqXqZAv8cUr8DQ2hY1ddXdivri48b8g2mDEhJkyId0/3iy5qn/QU4quvQq3Fb38L22wTpqMdMSK8t5oa+M1v4P77Yc6cUPI/77zQdq2gH9/GG4c8vuee8H35wQ/g+9+HF18sd8pEqkesEr+79y91QpLIrPVV/ctXh9MhqsOvvTZe4J8yJYz5rgRz54Zpb1NV98880zw0c5tt4OSTQ9X9dtvBiiuWO7Udy4gRsNtuMGlSqPYfNAiOOirMBbD66uVOnUjHphm226Cmpnmu/jhV/R2xOjwldXOYYm1XCl9+2Rzo6+vDLWtTd1bcdls45ZRQdb/tttC9e/nSWS1qa8NIh4MPDgH/kkvCzYBOPRWOOw5WaO19PkUklpxV/Wa2VzQnfurvFpf2S3LlaG1Vf0eoDs8l7rTXNTVw1VXw+ectb1eM5pA5c+DOO0N78uDBsPLK4bayEyeG9J56KkyfHi4IZs4MwWeXXRT029vKK4fJil58EXbaCU46KQz7u/32cGEtIsXVUon/X8C2wFPR39mmzU1xlh973+Fl69XfUok/idXhcY0aFZorWnp/NTWhFPezn4X82nnnMDPdj34E66zTvF2hzSFffAEPPdTcGe+FF0Lg6No1VNePHx9K9NtsE9ZJZdlwQ/jXv8IcAL/+dfhsDB8eLoQ337zcqRPpOFoK/P2B/6b9LRlaW+JPQnV4ocaNC0G5pcDftWuYinbuXLjttrD88pdhGTIk/NAPHhy/OaRXrxDoU230s2c3B/rttw/jxYcNg623VrVxkuyxR7hou/zycAvgLbaAI4+EPfesLXfSRDqEnIHf3d/N9rc0a+0EPj16hJJrPkm8W+SAAaEkni1o19aGZepUWHfdsG7wYDj7bHj99XABcOutoY09jkWLQjD/IroLbLduYez8WWeFznhbbaVAn3SdO4emnZ/8JHxOLr4YrrtuG954I8wLoBobkcLFHc4HgJmtYGbfM7ONM5dSJbCSpVf1x+ncN2pU85z+udTWwiGHFDed7WXEiDBlLYTAa+bU1YW7zM2enb16fv31w/C5J5+EDz6I94O+dGkYenfOOWGc/ZdfwgMPhAuHHXdU0O9IevcOfWNeegkGDfqS3/0uDAm85Ra1/4sUKu7d+b5jZv8CFgBvAC+mLf+JHqtOa6v6x42LF/hPOKG46WxP/40ah+68E6ZPn8ncuaG/QpxRCmutFW5iE8fSpfD734cqfc3+1vGtvz6cc85/eOCBMLRy5MjQjPPcc+VOmUjyxC3xTwaGAL8G9gR2SVuGR49Vp7VV/anq8M5ZGlhqa0Nv8qlTkzeUL119fXh/229f2P5xmzmS2BwibbfbbiHYX3YZvPxyaDI64ojmC04RyS9u4N8B+KW7/9XdH3D3mZlLKRNZqVpb1Q+hurt//1C6Tc3n3qVLy9XhSTJzZmhjLzQwd/TmEGm7zp3DzX7efDPUok2ZEmoEzjsv9P8QkZbFDfyfAgtLmZAkKmTK3vfegzfeCMOV5s4Nvdl32il+dXgla2yEp58OHewKVQ3NIVIcK60EF1wQSv677RaafjbcEG6+We3/Ii2JG/hPA36bmtBHgvSZ++JO2Xv33eFx773D46BBYYhbR/iheuyxcBObYcMKP0aqOaR79+UvADpKc4gU17rrhpEh06eHIZ4HHhjmiJg1q9wpE6lQ7p53AW4G3gPmAPcD/8xYbopznHItgwcP9mI67LDDHHA4zUPIdofxDrMdpkavheWoo45aZt+993avq/skbZtjo/3XWma/1DJ+/Pjlzr/PPvtk3TbbMmnSpOX233LLLWPvf+eddy63/5prrplj+7MdmhxW/HbdrFmzlts/3rm/5/A379lzidfUuNfVuR97rPujj34cO+3h472sWbNmxd53zTXXXG7/O++8M/b+W2655XL7T5o0Kfb+++yzz3L7jx8/Pvb+mZ89d/ejjjoq9v6V9tmbMWNGC5+91FLjcITDxw7uhx3m/uGHrf3sheXD9B3d/cMPP9Rnr0o/e+4t/e4tvxTyu5ftPbcFMMtzxMS4Jf5VgQbgeaAWWC1jqdLbaizN+HspLU1guHBhKJWss076IIjno8dBxU9euxsGPA3ML8Kx3gKO49VXP2bJEr4dHdCvX4zbH0oVWwpcAazHYYd9zA03wHrrhbkAFqqxUgSIWdXv7sPzLaVOaGVKD/xLoiV3ls6YEX58vvvd9MA/O3pM+pyk3YGtgars5ykVZx7HHfchr7wSOsyeeipssAHAj8udMJHyy1UV0JGWYlf1z5gxw93dzz03Vc0f/h4yxH2vvXLvN3as+4orui9atOz6AQPcR44sahLb3f33h3yYNq15XSqfpGXKp3jakk/19e5bbBE+o9tt5/7kk8VLVyXRZymeasgnWqjqj31bXjP7DrAPsDaw3Pxq7v6b4l2OJENNzbJ/t9S5zz3cgGS33ZafWS7VwS/JZs4MHRx32KHcKRFZ3tChYcTJ1VeH3v/bbBOGjp53Hqy9drlTJ9K+4s7ctx+h0fUS4AjggIxlZKkSWMnSA3+nTi2P43/ppTCUL9WbP93mm4cxyXHm8a9U9fVhMpWePcudEpHsOnWCww8Pw2lPPjkM+1t//XAzp2w3hRLpqOJ27juX0Ju/j7uv5e79M5bvlTCNFSuzxN/SOP7UML699lr+tUFRv77Zs5d/LQkWLICnnmrbMD6R9tKzJ5x7Lrz6Kuy7L5x+emj/v+66/BNwiXQEcQP/OsDF7v5FKROTNK2p6r/77hDg11pr+ddSgf+FF4qfxvbw+OPhdrkK/JIk/frBTTfBww9Dnz6h6n/77eGJJ8qdMpHSihv4HwM2KGVCkihuVf+cOWFym2zV/BDaGFdeObnt/PX1at+X5Npxx1BjddVVoTluu+3gpz+F998vd8pESiNu4P81MMbMDovu1Nc9cyllIitV3BL/ffeF9bkCv1lo509y4N9yy+Z7D4gkTU0NHHYYvP46/OEPcOutofr/tNPCVNQiHUncwD8b2Ay4EngfmJdlqTrZ2vizlfjvvhtWXRW23jr3sQYNghdfDFPeJona96Uj6dEDzjortP//8Ifh7w02gGuuUfu/dBxxh/MdTphWUNKYNf+dqurPLPEvWQLTpoVJRDrlntSPQYPCncXeeAM22qg06S2FJ56Ar79W4JeOpW9fuOEGOO64cEOoww6Dv/0N/vIXNWlJ8sUK/O5+VYnTkUhxqvqfego+/zx3NX9KqoPf888nK/DX14f3veOO5U6JSPFtv33ovHr99fC734XP+Y9/DOefHy4ORJIoblW/ZBGnqv/uu8P6PfZo+VgbbghduiSvnV/t+9LR1dSEHv+vvQbjx8Odd4bq/1NOSfbcG1K94k7g8z8z+7SlpdQJrUTZevVnlvjvvjtUDfbu3fKxunSBjTdOVuBfuBCefFLV/FIdVlwxjPl/7TUYOTLMBbD++nDllWr/l2SJW+K/JMtyA/AJsBCYWJLUVbh8Vf0ffhgCeb5q/pSkTd2bat8fOrTcKRFpP+usA9deG5oA+vYNswFutRU89FC5UyYST9w2/tOzrTczA/4JJKwvenHkq+q/557w2JrAf9VV8PHHsMYaRUtmyah9X6rZttuG4H/DDaH9f+jQUBPwpz9B//7lTp1Ibm1q44/uADQZ+EWc7c1spJk9Zmafm9kiM3vNzP5gZl1i7l9jZs+YmZvZPm1JezFkVvVnlvjvvjuUCDbeON7x0jv4JcHMmbDFFtCrV7lTIlIeZnDwwWH435lnhov9DTcM9wL46qtyp04ku2J07vseECtwA6sAM4AjgRHAP4BTgAtj7n8kkGXS2/LILPEvWACffBI6utXUwB13hHHBb70V73gDB4bHuIG/oQHGjm0+X11deN7Q0Lr3UYhFi0JVv9r3RaB7dzj11DAB0EEHwR//GNr/J0/OPY23SLnE7dw3NstyvJlNAi4A7opzHHef5O6nuPtt7j7D3c8nBP1RUbNBS2noDZxDuFCoCOmB/4UXQgn/q69CT1+PZj149dUQ0KdNy3+83r1DDUGcwD9tWjju5MnN55s3LzyPe762eOIJWLxY7fsi6dZaK9z696mnYMAAOOqocNfK+vpyp0zSC0q77DK0XQtKlSZuiX9iluU8YBfg78Ss6s/hc+LVGJwFPAo82IZzFVV64J8wIfuV/ZIloSZg5Mh4H7BBg/LfrKehIRxvwYJwc5x0TU2tO1+hZs4M1Zw77VS6c4gk1VZbwSOPwI03hnt1DB8O+++/7HeynDV21Wb5gpK1a0Gp0sQK/O5ek2Xp5u7ruftv3H1+a05qZp2iOf53BH4JXBr1F8i1/UDgZ8CJrTlPqaUH/nxT7TY1wUUX5T/moEFhuND8FnJ0woTlA36h5ytUfb3a90VaYhYm+3n1VTj7bLj//tDf5ze/galTy1tjV00qoaBUadrcxm9mw82stR/T+dHyMDATOCnP9n8DLnH3NwtIYsmkB/587XhNTTBlSv5jDhoUfgT+85/c21x7bbzAH+d8hVi0KPRmVjW/SH7duoXJft54I9z174IL4IADFIjaSyUUlCpNi8P5zKwXsCewDvA2cIe7N0WvHQD8FtgSeL2V590e6A5sDZxGaDoYmyMNBxFuCbxva05gZmOAMQB9+vShvoiNbI2NjdTX1/PSS6sCm8beb948p75+ZovbLFzYFdiWm256jYUL/5vj/EOBFrtExD5fIZ5/fiUWL96CVVZ5kfr6z3Nul8onaZnyKZ6OkE+HHgoNDRvz0EOr0dJ3ePHipYwb9xHHH9+6sk56Hrmnhhcb7kRL+t/LPk9t19p9li5t3f6F7JPvnC3tc8UVG9HU1MKNUgiB/8orv2HkyEdald+J5e5ZF8Ld+P4LLE1bZgF9CW3tS4AXgZ8ANbmOk28BDiXcAGhAltdqCXcDPAHoFS0Do+1/DPSMc47Bgwd7Mc2YMcPd3W+9ddmPWb6lri7/sZcuDdsdfXTubXr2LN75CnHGGe5m7l980fJ2qXySlimf4uko+RT3+wvuvXq5r7RS+C737Om+4oru3bu7d+vm3rWre5cu7rW17p07u3fq5G62tFW/SVqal5qacn8yiguY5TliYksl/nOBr4AfAS9EAf9vwNPACsBh7n5tEa49no0e+wOZlVsrAmsTev5nDvm7Mdp+3SKkoSCZ4/hbqu6vrYVDDsl/TLP8HfxGjQptgS1VX8U9XyHq60Ma801DLCLLa2yMv+2oUeF3xqx5aen5e++9R//+fWNvH3ebJD/fYYeW+0yl9OhR+P80aVoK/EOAX7n7k9Hz18zsGOANYEyRgj5A6iaXb2d5rREYnrFuDcJ0wb8HphcpDQVJD/ydO+cP/CecEO+4gwbBFVeE42W7le+4cWHIUL7AH/d8rbF4cWjfP/ro4h9bpBr06BHv5j51deFWwK1RX/82w4b1LSxhHdShh5a3oFSJWurc1wd4J2Nd6nmeAWfZmdm9ZnaimY0ws++b2RnABOAmd2+ItnnTzK4AcPdv3L0+fQGeiA73YtpFSVmkB/7TTw/BP1NtbZjcY+rUMK43jkGDwhVqrs49AwaEqX0hXNFmau35WuOpp0LnPk3cI1KYUaPC70JLqi0QldK4cfHyuxQFpUqVr1e/51hf6Nz8TwOjgZsJc/zvC5wMpH/EOwMt98SoEOmBf7vtYPToUEJfYYWwrq4OxoyB2bNhxIj4x9188/DY0kQ+L74YHg84oHkc8IorhnXHH9+687VGfb3G74u0hQJR+xowIBSEundfPt8LKZh1BPkC/30Zt95NdTN/sJDb8rr7qe6+qbv3cPde7r6lu//No5EC0Tb93H10C8d4x93N3f8V55yllNnGv+qq4XHs2FCdN3cuTJzY+g/UxhuH2oNc7fyffw5/+UsY8nPTTeE8S5aEtsOddw7D+PINXylUfX0YZ7zyyqU5vkhHp0DU/kaMCAWwMWNCQcnMCy6YdQQttfGf0W6pSKhsd+dLzdTXvXvhx+3aFTbaKHeJ/4ILQpA//fTlXzvpJNh3X7j55nDzkGJKte+PGVPc44pUm1QguuiicKHe2BgKC4ccEkr6CvrFN2BAKIhNnAj19TMZVsXtlTkDv7sr8OeRGfhTd+ebP7+52r1QgwbBg1kmJ/7kk9Dh5yc/gU02Wf71vfYKdwe74IKwTct3QGidp5+GhQvVvi9SDOmBSKQ9FePufFUrs6o/1QN//vy2lfghtPN/9BF8mtGI8sc/hpL3+PG503TiiaG2INuFQ1uk5k5R+76ISHIp8LdBtqp+CEN1ilHih2Xb+T/8EC69NAxPWX/93PuOGgVrrBFK/cWUat9fZZXiHldERNqPAn8bpFejp6r6Idyatxglflg28J9zTmhKOPXUlvddYQX45S/DTUHy3ekvrq+/hsceUzW/iEjSKfC3Qa6q/mKU+FddFdZeu7mD37vvhkkojjwS+vfPv//RR4c0TJjQtnSkqH1fRKRjUOBvg2yd+yAE/raW+BsawoXEDTeE4663Xijt//Sn8fbv3TtcJNxwA7z/ftvSAmrfFxHpKBT426BUbfzTpoW29Pff59u7UKXG5e+xR/x7dZ9wQtj3r38tPC0pM2fCZpuFmggREUkuBf42KEVVf0NDmJhnwYLUrSqbLV3aunt19+0LBx4Il18eJvkp1Ndfw6OPqppfRKQjUOBvg1xV/UuXFl7VP2FC/ln3mprCxB9xnHRSuBCZNKmw9ADMmhUuOIYOLfwYIiJSGRT42yBXVT8UXuK/9tp4gX/KlHjH22IL2HXXUN3/9deFpWnmzPC4886F7S8iIpVDgb8NclX1Q+El/rj36m7NPb1PPDFMBnTDDYWlqb4eNt0UVlutsP1FRKRyKPC3Qa6qfii8xN+jR3G3g9AhcLPN4M9/Dp39WqOpCR55RO37IiIdhQJ/G7RU1V9oib8U9+o2C6X+//wH7r23delR+76ISMeiwN8GmVX9xSjxl+pe3QcdBGutFUr9raH2fRGRjkWBvw0yp+wtRom/VPfq7tIFjj8epk+HZ5+Nv199fbgL4Oqrt+58IiJSmRT426AUvfqh+V7dY8ZAXV04dl1deD57dni9EEcdBT17xr95T6p9X9X8IiIdhwJ/G7RU1d/WKXtT9+qeOzdM1Tt3bnje2pJ+upVWgp//HG6+Gd55J//2zz4bbjGsjn0iIh2HAn8blKrEX0q/+lVooogzAVBqfn6V+EVEOg4F/jYoxTj+Ult7bTj4YLjiCvjii5a3ra+HjTdW+76ISEeiwN8GpRjH3x5OPDFU4V92We5tvvlG7fsiIh2RAn8blGIcf3vYbLMwqc/FF8OiRdm3efbZMDug2vdFRDoWBf42KGXnvlI76ST45JNwb4Bs1L4vItIxKfC3Qa4Sf9euy75WiXbZJdzA589/Xv72vxAm7tlwQ+jTp/3TJiIipVPh4amy5Qr8ldy+n2IWSv2vvQZ3373sa998Aw8/rGp+EZGOSIG/DXJ17ktC4Ac44AD47neXn9Dnuedg3jwFfhGRjkiBvw1SgT71mCrxV3r7fkrnzmHO/4cfhiefbF6v9n0RkY5Lgb8NUnP1Zwb+pJT4AY48Enr1WrbUP3MmbLABrLFG+dIlIiKlocDfBpkBP/U8KSV+gB494Jhj4JZbwi2B6+pCm/9bb8HYsdDQUO4UiohIMSnwt0Guqv4klfgh3H0P4PrrQ9s+hBv0TJ4MAwfCtGnlS5uIiBSXAn8bZAb+JJb4GxrCXf8A3Jd9rakJFiyAkSNV8hcR6SgU+Nsgs6SfxBL/hAkhwLekqSneTX1ERKTytWvgN7ORZvaYmX1uZovM7DUz+4OZdWlhn63M7Eoze9PMFkT7jDezru2Z9myS3qsfwsx9cQL/lCntkx4RESmt9i7xrwLMAI4ERgD/AE4BLmxhnx8DA4Dzgb2AS4BfA9eVNKUxpAf+hgY477zwfNKk0Ekuq2bdEgAAEUxJREFUCZ3jGhuLu52IiFS2zu15MneflLFqhpnVAcea2XHuma3MAJzv7v9Le15vZouASWbW193fLVmC80gF/m++CZ3gvv66+bV580LnuKuvhqlTYcSI8qQxnx49mjv05dtORESSrxLa+D8Hclb1ZwT9lOeix7LeKT4V+L/6KnSC++abZV9PQue4UaOgtrblbWpr4ZBD2ic9IiJSWmUJ/GbWycy6m9mOwC+BS3OU9nPZHlgKvFaSBMaUmsAnn0ruHDduXLzAf8IJ7ZMeEREprXKV+OdHy8PATOCkuDua2RqEfgFT3P2r0iQvblribVfJneMGDAhNEd27L38BUFsb1k+dGrYTEZHks9YVtIt0UrMtge7A1sBpwPXuPjbGfl2AfwNrA4PdfU4L244BxgD06dNn8I033liMpAPQ2NhIj6jRe/jwoUD+KwAzZ/r0mUVLQ7F9+GFXbr55bR54YA0WLuxEt25L2H33jznggA9Ya61FBR0zPZ8kN+VTPMqn/JRH8VRDPg0fPvwZdx+S7bWyBP5lEmB2KHA1sK6752wJNzMDbgB2B3Zw91fjnmPIkCE+a9asNqc1pb6+nmHRrevilvrr6mDu3KIlIRHS80lyUz7Fo3zKT3kUTzXkk5nlDPyV0Lnv2eixf57tLgJ+CPywNUG/1Gpi5KA6x4mISKWohMC/Q/T4dq4NzOxk4DhglLs/0i6piqlz5/ylfnWOExGRStGu4/jN7F5CG/1LwBJC0B8H3JSq5jezN4GZ7n5E9Pxg4FzgKuBDM9s27ZANOYb7tZvOnWG11WDOnNCJL30WvNrasKhznIiIVIr2LvE/DYwGbgb+CewLnAykV4R3BjqlPf9+9DgaeDxj2bukqY2hpgZWWglmzw43u6mrC+vq6sLz2bMrd/IeERGpPu09c9+pwKl5tumX8Xw0IehXpJqasAwYABMnhkVERKRSVUIbf6KlAr+IiEgSKGS1UU1N8135REREKp0CfxupxC8iIkmikNVGCvwiIpIkClltpKp+ERFJEgX+NjJTiV9ERJJDIauNVNUvIiJJopDVRqrqFxGRJFHgbyOV+EVEJEkUstpIgV9ERJJEIatADQ0wdiy88w48+GCYm3/s2LBeRESkUinwF+DJJ1dm4ECYPBncw7p588LzgQNh2rTypk9ERCQXBf5WamiA00/fhAULlr0FL4TnCxbAyJEq+YuISGVS4G+lCROgqcla3KapCS66qJ0SJCIi0goK/K107bWwZEnL2dbUBFOmtFOCREREWkGBv5UaG4u7nYiISHtS4G+lHj2Ku52IiEh7UuBvpVGjoFOnpS1uU1sLhxzSTgkSERFpBQX+Vho3DmprvcVtamvhhBPaKUEiIiKtoMDfSgMGwOmnv0T37iHAp6uthe7dYerUsJ2IiEilUeAvwDbbfMHs2TBmTJixr6YmPI4ZA7Nnw4gR5U6hiIhIdp3LnYCkGjAAJk4Mi4iISFKoxC8iIlJFFPhFRESqiAK/iIhIFVHgFxERqSIK/CIiIlVEgV9ERKSKmHvLs9B1BGb2P+DdIh5yVeCzIh6vo1I+xaN8ikf5lJ/yKJ5qyKe+7r5atheqIvAXm5nNcvch5U5HpVM+xaN8ikf5lJ/yKJ5qzydV9YuIiFQRBX4REZEqosBfmMvLnYCEUD7Fo3yKR/mUn/IonqrOJ7Xxi4iIVBGV+EVERKqIAn9MZraxmT1oZgvM7CMzO9PMOpU7Xe3FzA4wszvN7EMzazSzZ8zsJxnbmJn93szeN7OFZvaQmQ3KcqyqyEszWyvKKzezHmnrqz6fzKyzmf3OzN4ws8Vm9oGZXZSxjfLJ7CAzezb6HH1oZteY2XcytqmqfDKzdc1skpm9YGZLzKw+yzZFy5O4x0oUd9eSZwF6Ax8B/wZ2B44G5gNnlztt7ZgHjwPXAwcCuwB/Bhw4Lm2bk4GFwC+A3YB7CGNl16jGvIzy6+Mon3oon5bJmynR+/s5MBQYBZybsU1V5xPwg+izMxHYNcqjd4BngZpqzSfgh8D7wM3AK0B9lm2KlidxjpW0pewJSMIS/ePnAHVp634DLEhf15EXYNUs664H3o7+7grMBU5Le31F4H/pX6RqyUtgJ+AL4ETSAr/yyQH2BJqAjVvYRvkENwLPZKxLXQxsVK35xLIXPVMzA38x8yTusZK2qKo/nhHAfe7+Vdq6G4FuhNJKh+fu2Wa5eg5YPfp7e6AO+GfaPvOBuwj5l9Lh8zKqKvwbcCbLzw6mfILDgenu/nIL2yifoJYQdNJ9GT1a9Fh1+eTuS/NsUsw8iXusRFHgj2dD4NX0Fe7+HuHKcMOypKgybA+kfrw3BJYAb2Rs8wrL5lE15OXRhJLCJVleUz7BNsDrZjbRzL6K2ldvzWi7Vj7BP4CdzOxQM6szs/Xh/7d39kFeVlUc/3xbUkR8QVTUtLDS8qUyHZkoJbDUwRdMsjT1D2wcsTQtHTCjRsZ3USxnUosmdabSSlMkUZFBgaUQX9JKHFAnV9IAERZFeVHx9Me5v+Hh2ee3++yyy7r7nM/Mnd/e89x7n3PPb/Z37vvlSuCxTKMp7NSSzrRJ2bJ6FOH4yzGATS3tLM3pWeWQ9DV8rq3m3AYAb5vZxlzSZqCfpG0y6XqtLSUNBK4ALjKz9wqShJ1gD2AMcAhwGnAWcBhwn6RaT7bydjKz6bidpuA9/8VAAzA6k6zydiqgM21StqweRZ/uVqAHUXTggerIezWSBuPz+/eb2R2ZR/VslH/Wm215FbDAzB5sJU3V7aQUTjKzlQCSlgJz8IWjs1K6SttJ0gjgV8BNwEPAIGAi3kD6esYZVdpOdehMm5Qtq8cQjr8czcDOBfKdKG4x9lok7YL/CC3BVxnXaAZ2kNSQax3vDKzN9H57rS0lHYTPXw+TVKtjv/S5k6SNhJ3A6/afmtNPzAPeBQ7EHX/YCSYD08zskppA0rP48PRJwL2EnYroTJuULatHEUP95VhEbj5H0j746s5FhTl6IZL6AQ8A2wDHp0UuNRbhw5CfzmXLz6P1Zlvuhy/Imo//YDSzaSrkVXzBX9jJ50eLEFBbuBV28no9mxWY2WJ8a9mnkijs1JLOtEnZsnoU4fjL8RBwrKQdMrJT8X/AOd2j0tZFUh983+x+wEgzez2X5O/AW8C3Mnn6ASfi9qvRm205DxiRC9elZ8cB1xN2Am88fl7SrhnZMLzR9M8UDzvBK8ChWYGkA/BV501JFHZqSWfapGxZPYvu3k/YEwK+wGMpMBM/wOEc4G168D7ODthgCj6fdQHwpVzYNqW5FF8Rex5+4Mh0fDvboKraEl+cVXSAT2XthG+PWoKPjJwInI4fyDIzl67qdroQHwGZnOp2Br7A72Vg+6raCZ8+OyWF+cDCTLxfZ9ukTFk9LXS7Aj0l4HOPj+KtwaX4yu2G7tZrK9a/KTmwojA4pREwAR/WXgc0Al+ssi0pdvyVtxM+dPogflJaM3AHMCCXptJ2SvX/HvCvZKfXgD8Bn6yynYDBW/O3qGxZPSnE7XxBEARBUCFijj8IgiAIKkQ4/iAIgiCoEOH4gyAIgqBChOMPgiAIggoRjj8IgiAIKkQ4/iAIgiCoEOH4g6AOkiZKMkkzCp7dI2n2VtRleNLl4K31zvYg6QBJjZLeSXoObkfeGyQ1ZeJjUhn9u0DV/LvHSxpeIDdJ53f1+4OgOwjHHwRtc4ykw7tbiQ851+MXl4wChuKHoXSU6amMtZ2gV1uMB4ZvhfcEwYeGuJ0vCFpnFX5i1wTgG92sS5chqa+Zrd+CIj6L3yQ3q82UbWBmK4AVW1pOd9AJdgyCLid6/EHQOgZcDYyS9Ll6idK0wBsF8s2GjCU1paHtH0taKulNSZPlHCdpoaQ1kqZKGlDwqr0kPZCG1JdIOrfgnUdImiNpraSVkn6TvYgkM5Q+RNJsSeuAca3U7RBJs1J5zZL+IGlQejZYkuG3xf0olTu7lbJ2lnRn0n+ppAkFaVoM9UvaTtIkSa9I2iDpZUnX5PKdney3IaUbX0+PlL4JGAhclt5nuWH/BklXS1oh6XVJN0vatowdJR0laYGk9ZKWS7olV59XJF2aiY9NZV2QkV0s6bXW6hAEHSEcfxC0zd3AC3ivvzM4DRgCnAVMAi4CbsTPCf8ZcC7wVeCagry/xc9uH43fDnarpBNqDyV9Bb/Pfhl+ackP8ZsBby8o6y78przj0mcLJO0GzMYvRjkd+EHSbaakbfAh/aHpfXemv7/fSt1vB0Ymvc4Bjkn2qIskAffj59bfnPS9DNg1k2YccCswFTgh/X1FG/P0JwNv4jYdmsI/Ms8vBvYCzsSnMsbiF+fk2cyOkg4EHsYvcvlm0vV04J5Mnkb8RsIaw4D1wJE5WWMr+gdBx+juywIiRPiwBmAi8Eb6ewywEdg/xe8BZhelzZVhwPmZeBPwEpmLQIAngPeBfTOyScDyTHx4KmtKrvyZwOOZeCPwWC7NUSnvwZm6GHBhCRtcC6wGdszIhqT838nV64Y2yjoo5Ts1I+uPT6c0ZWQ1/fqn+LEpPqpOuTvit6pdlpNfjjdI6l5EgzvniXW+t7k52dScrQvtCPwReDH3HX87pR2a4mPxRsdHUnwJ8EtgWYor6XZed/8fROh9IXr8QVCO3+M/zpe2lbAEs81sYyb+Eu74Xs7Jdku96iz35eL3AodJapDfEz4U+LOkPrUAzAPeAw7L5Z1eQtchwCNm9lZNYGZP4I7+iBL5s9QWSE7LlPU23nhpjaOAVWY2rc7zocD2wN25ej8KDAL2bqeeNR7JxZ+vU1bejkOA+3Lf8V/wxl3NZo14g+ULaQfE3nhjb1dJ++GNpIFEjz/oAsLxB0EJzOx9/If5TEmf2MLiVufi79aRCcg7/tcL4n3wYe8BQANwC+7oa2ED8FFgn1ze5SV03bNOuuXALiXyZ9kDWGNm63LyfJ3yDKT1XQK1If+FbF7vx5I8X++yFH0nfQvS5e3TwmapEbCSZDMzex7v0R+ZwnNmtgR4NiNbDTzXQd2DoC6xqj8IynMb8FPgkoJn68k56TqL87aU3Qvi7+NOpC8+nDwRv+s+z/9y8TJ3ci8teCd4T/rpEvmzLAN2kLRdzvkXlZ9lJe5M67EqfZ5AcSNlcXkVO0Teji1sJqkBb8CsyojnscnBz02yxiTrC/zNzD7oCoWDahM9/iAoiZltAG4AvktLR/Qq7tQ+lpEd0wVqnFwQf9rMNprZO8DjwGfM7KmCkHf8ZVgAHJvbFXA4MBh3XO3hyfQ5KlNWf+DoNvLNAnbJLmLMMR9YB+xVp95rWim7Xi9+S1gAnJycfY3ReEcra7Oakx/GJsc/l009/hjmD7qE6PEHQfv4NfAT4MvAnIz8Ydz53CZpMrAvvjq/sxkp6ar07tG40zwp83w8MEvSB/gCxDXAx4HjgQlm9kI733cjvpp+hqTr8MV41wL/xuetS2NmCyVNw3ci7Ij3jMfR9kE9M4EZwJ2SLsdX3u8JDDOzsWa2WtJE4KY0DTMX79TsD4wws3xjKcsi4HhJD+MLBBe30VAow5XAM8BUSbfi8/fXATPMbH4m3VxgMj56UnP88/CtkRCOP+gioscfBO3AzNYCPy+Q17Zu7Y2v/j4T38LV2ZwNHMqmbWvnZRe9mdk8vAe5G/A74K94Y+C/lJvT3wzzw3RG4FMZd+Hb6RqBo83s3Q7oPwZfNPcLfBvdLHwVfGs6GD6yMQXfBvgQ7lzfyKSZhG8PHIlv/bsLOIO2nec44B18gd6TtFwA2W7MbGHSY3d88eWVSZ9TckmfwRsbL5rZspR3Bd4YWQ88taW6BEER8v+pIAiCIAiqQPT4gyAIgqBChOMPgiAIggoRjj8IgiAIKkQ4/iAIgiCoEOH4gyAIgqBChOMPgiAIggoRjj8IgiAIKkQ4/iAIgiCoEOH4gyAIgqBC/B/igs4upYN56AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,5))\n",
    "plt.title(\"Running average of dice throw simulation\",fontsize=18)\n",
    "plt.plot(n_throws,av,color='blue',marker='o',markersize=10)\n",
    "plt.hlines(y=3.5,xmin=0,xmax=1100,linestyle='--',lw=3)\n",
    "plt.xticks(fontsize=15)\n",
    "plt.yticks(fontsize=15)\n",
    "plt.xlabel(\"Number of dice throw\",fontsize=15)\n",
    "plt.ylabel(\"Running average\",fontsize=15)\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "2TTVgVHpumZG"
   },
   "source": [
    "## Confidence interval\n",
    "\n",
    "### Some essential definitions\n",
    "\n",
    "***Population***: The whole collection of which we want to measure some property. We can (almost) never get enough data about the whole population. Therefore, we can never know the true values of population properties.\n",
    "\n",
    "***Sample***: A fraction (subset) of data from the population, which we can gather, and which helps us estimate the properties of the population. Because we cannot measure the true values of the population properties, we can only estimate them. This is the central job of statisticians.\n",
    "\n",
    "***Statistic***: A statistic is a function of a sample. It is a random variable because every time you take a new sample (from the same population) you will get a new value for the statistic. Examples are the sample mean or the sample variance. These are good (unbiased) estimates of the population.\n",
    "\n",
    "$$\\bar{X_n} = \\frac{1}{n}\\sum_{i=1}^{i=n}X_i, \\ \\ S_n=\\frac{1}{n-1}\\sum_{i=1}^{i=n}(X_i-\\bar{X_n})^2$$\n",
    "\n",
    "### Confidence interval\n",
    "\n",
    "A range/bound around the statistic (of our choice). We need this min/max bound to quantify the uncertainty of the random nature of our sampling. Let's clarify this further with the example of the confidence interval for the mean.\n",
    "\n",
    "Depending on where and how we are drawing the sample, we may get a good representation of the population or not. So, if we **repeat the process of drawing the sample many times**, in some cases the sample will contain the true mean of the population, and in other cases, it will miss it. \n",
    "\n",
    "**Can we say anything about the proportion of our success in drawing a sample which contains the true mean**? \n",
    "\n",
    "The answer to this question is found in the confidence interval. If some assumptions are met, then we can calculate the confidence interval that will contain the true mean (when we sample a large number of times) with a certain fraction.\n",
    "\n",
    "The necessary formulas are given below. We won't get into details about this formula or why the particular t-distribution is used in this equation. Readers can refer to any undergraduate level stats text or excellent online resources to understand the rationale.\n",
    "\n",
    "<br>\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/tirthajyoti/Stats-Maths-with-Python/master/images/Conf_interval.png\" width=600 height=420></img>\n",
    "\n",
    "**Source**: https://psu.instructure.com/courses/1844486/pages/chapter-3-confidence-intervals\n",
    "\n",
    "<br>\n",
    "Confidence intervals are a calculated range or boundary around a parameter or a statistic that is supported mathematically with a certain level of confidence. \n",
    "\n",
    "This is *__different__* than having a 95% probability that the true population proportion is within our confidence interval. Essentially, if we were to repeat this process, 95% of our calculated confidence intervals would contain the true proportion.\n",
    "\n",
    "The equation to create a confidence interval can also be shown as:\n",
    "\n",
    "$$Population\\ Proportion\\ or\\ Mean\\ \\pm (t-multiplier *\\ Standard\\ Error)$$\n",
    "\n",
    "The _Standard Error_ is calculated differenly for population proportion and mean:\n",
    "\n",
    "$$Standard\\ Error \\ for\\ Population\\ Proportion = \\sqrt{\\frac{Population\\ Proportion * (1 - Population\\ Proportion)}{Number\\ Of\\ Observations}}$$\n",
    "\n",
    "$$Standard\\ Error \\ for\\ Mean = \\frac{Standard\\ Deviation}{\\sqrt{Number\\ Of\\ Observations}}$$\n",
    "\n",
    "Therefore, the $(1-\\alpha)$% C.I. (**C**onfidence **I**nterval) is given by,\n",
    "\n",
    "$$ C.I. = \\bar{X_n} \\pm t_{\\alpha,n-1}*\\frac{S_n}{\\sqrt{n}} $$\n",
    "\n",
    "where,\n",
    "\n",
    "$\\bar{X_n} = \\text{sample mean}$, $S_n = \\text{sample standard dev}$, $n = \\text{number of samples}$, $t_{\\alpha,n-1}$ is the t-statistic for parameter $\\alpha$ and degrees of freedom $(n-1)$.\n",
    "\n",
    "### What is the practical utility?\n",
    "Be careful about the definition and the process to understand the true practical utility of the confidence interval.\n",
    "\n",
    "When you are calculating a 95% confidence interval of mean, you are not calculating any probability (0.95 or otherwise). You are calculating two specific numbers (min and max bounds around the sample mean) which creates a range of values that will contain the true population mean (unknown) if we were to repeat the process.\n",
    "Here lies the practical utility. We are not repeating the process. We are just drawing the sample once and constructing this range.\n",
    "\n",
    "If we could repeat the process a million times, we would be able to verify the claim that the true mean lies inside this range in 95% cases.\n",
    "\n",
    "But sampling a million times can be quite expensive and downright impossible in real life. So, the theoretical calculation of the confidence interval provides us with the min/max range, just from one draw of the sample. This is amazing, isn't it?\n",
    "\n",
    "### But in the simulation, we can experiment a million times!\n",
    "Yes, simulation is fantastic. We can repeat the sampling process a million times and verify the claim that our theoretical confidence interval truly contains the population mean, approximate 95% of the time.\n",
    "\n",
    "Let's verify it using a real-life example of factory production. Let's say in a factory, a certain machine produces 20 tons of product on average, with a standard deviation of 5 tons. These are the true population mean and standard deviation. So, we can write simple Python code to generate a typical production run over a year (52 weeks) and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "40xYmAb_umZG"
   },
   "outputs": [],
   "source": [
    "num_weeks = 52\n",
    "production = np.random.normal(loc=20,scale=5,size=num_weeks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 288
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1635,
     "status": "ok",
     "timestamp": 1567211402670,
     "user": {
      "displayName": "Tirthajyoti Sarkar",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AAuE7mD6d7dlMqdpzL4sermIF1ujmpSRxY2WnE4tuB-UsQ=s64",
      "userId": "01914075970409030121"
     },
     "user_tz": 420
    },
    "id": "BGgac_VsumZI",
    "outputId": "63758a0d-a59f-403b-d367-2e54187ed75c"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAEdCAYAAABACGBRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydd7gdZbX/P+vknCTnpJIeAiEkgUCAkNBDKAGpNooBRYLgFeIloIjozy7ovagXBRRRKUFAsNIEhYBCckIVCDVAAuGkQArpvZ7y/v5YM9lz5uwys/fM7PZ+nmee2Xvqu2dPWbPW911LjDFYLBaLxWKxWIpDTbEbYLFYLBaLxVLNWGPMYrFYLBaLpYhYY8xisVgsFouliFhjzGKxWCwWi6WIWGPMYrFYLBaLpYhYY8xisVgsFouliFhjzFJWiIgJMCyKcH8XOdscFtU2Pdue6Gx7YoBlvysiH4hIi4i8HkM7rhERez/wICKLROSumLad9piLyDDnnLgojv1aLJbSpLbYDbBYQjLe9/0h4A3gGs+0HRHu71Fnn8sj3GYoROQI4Frg58DfgU0R72IicDXwv0BbxNu2pGci6Y/5cvR8aypCmywWS5GwxpilrDDG/Mf7XUR2AKv90yPc3ypgVRzbDsH+zvgWY8yCorYkICLSxRgTpVFcVvvPF6fNsZzL5UC5/m9BqfTfZ8kfG5awVAwi0kVEVonIjWnmueHG/Zzvd4nIEhE5WkReFpHtTljqKxnWG+abfomIvCoi20RknYjMEpGjPfN/5MzfICKrRWSGiByVx29qBO5yvjY5bbnGmXe5iLwgImtFZL2I/EdEPpFmG91E5Gci0iQiO0TkIxF5QEQGOtu62lm02Q31etYdLCJ/cH7DDhF5U0QmZzhGx4nIfSKyHnhRRL7hrNPft7yIyAIR+XOO325E5FoR+Z7zX20TkadFZKz/GInIsyLyKRF5zTHQpzrzjhCRJ0Vks4hsEZGnHE+jf19XOP//dhGZLSLHplnmGu+x8Uy/yx8az/eYZwpTishkEXnDad9qEblHRAb7llkkIveKyOdEZK7ze2eLyDHZjrNn/dOc82mbc97+XURGeeb/VkRWiEitb70uzjXwS8+0fiLyOxFZ6vz+eSIyxbde2vMmS/tGOr97odPGBc4+dsvxuyY5+zk4zbxGEXnB871WRL7jtHeHiCwTketFpKtvvR9JjutbUjKEs0XkdhFZBazI1lZL9WKNMUvF4Lxx3glc6L95Al8GZhlj5nmm9QT+CtwNnAk0Ajf5H4R+ROQXwG3Aq8C5wGTgaWCoZ7EhwI3Odi8CVgJPi8iYkD9rKvBT5/PZaAhrmvN9mPP5HOCzwGzgnyJyuqetnYF/A19FjbpPApcDa4HdnPXvcBY/xtn+eGfdbsAs4HTgu85vmQPc43+wOvwRWAhMAr4N/B4NwX3Rt9wpwN7ArQF+/xeAjzttvggYCDwlIn18y+0L3AT8GjjVWWaM0/7dnHW/gP7ns7wPZhH5EvBLYKbzG+8C/uysF5pCjnmG7U0B7gHmoufAt53fOEtEuvsWPxa4CvgBek50Qs+J3jnafBoakt/srHcpcCDwrIgMcRb7AzAA/f+8fBLo7bQREekJPAd8ApUPfAL4B/A78b3sOPjPm0zsDiwBvob+/h8DHwMey/bb0ND+MvQesAvH0Dye9ufhvcD3gT857f4p8CWnjV7CXN+/BgS4wFnWYumIMcYOdijbAVgE3Ov5vjfQClzgmTYGMMDnPNPu8k9zpv8bWAyI8/0iZ7lhzveRzvZvCNHGTqgk4F3gV57pE51tT8yx/sXeNmRYpsbZx7+Ahz3T/8tZ99NZ1r3GWabWN/3ydO0DnkQfPp18x+jGNNu+C3jfPZ7OtAeBeQGOmwFWA90804YBzcD/eKY1okbfWN/69wPrgd6eaT1Ro+hBz3H7EHjct+5nnf3f5T9OGX7jooiO+TBn+kWec2cFMNO33DHOcl/1XQvrgN080w5zlvt8jmM9G5jvbQ96LTV7z3XgPeDPvnX/Drzj+f4DYDuwj2+5253/szbXeRPwuqr1HIdxOZa9BtjgO5ducI5XvfP9WGdbX/Cte74zfWyGbee6vh/K5/fZoboG6xmzVBTGmIXAE7R/C/4yqvt60Ld4K/CAb9pfUA/XENJzEvoAvy1bO0TkJBGZKSJrgBb0obYvMCrbemEQkUNF5J8issKzj5N9+zgF+MgY80geuzgOWGqMafRNvxfoD4z2TX8ozTZ+C4xAPRg4obVPEcwrBvCYMWaL+8UYswjVVPk9SYuMMf5epscB/zTGrPesvxF4BPWIAOzhDH/zrfsAekzzoZBj7mcU6o1q55kxxjyLvjQc71v+BWPMOs/3Oc54KBlwPKCHAH81xuz6zc619JxvH/cCZ4hID2fdPqjn9A+eZU5Dw40LnbBfrRPafALoS7DzJl07O4v2Kp4nItvQ8/0ZZ3au6+o2oAE4z9lWV+BC4A/GmG2edu8EHvC1+1/O/OM8bQlzfQf6fZbqxhpjlkrkt8AEETnQedBMBu40xuz0LbfOGNPsm+ZqOjIZY32d8ZJMOxeRQ9DQyWY0xHEUcDja69MfPs0LEdkTeAroA3wFONrZx+O+ffQFlua5mz6k70X6kWe+lw7LGmNeQr0u/+1Muhh9eN0dsA3pNDYr6Pj/pGtntva7IUhXd9VuP45RsiZgG/0Ucsz9uMc40+/w/wdrvV9MSiye7bzbDQ2jBdnHPc62JjnfPwfU0d5YHIAaLs2+4T5nfl/aE7Sn8k9RD9e9aAjxCDRsCzmuK2PMMuBhUufhOejv8r4UDAA6o9ett90rve3O4/ouWk9sS/lge1NaKpHH0JDNl9EbZA/Se7J2E5E6n0E20BlnepiudsZD0LBEOj6DGhxne7ftCI3XZ1gnLKcBvYBzjTG7DEMRaUjT3gPz3Mda0r/pD3LGfmOlg7jd4XfArY726GLgPmPM2gzL+hmYYZr//0m377WetnoZRMpocR+U7fbjeET8RsN2Z15nn2HvX66QY+7HbWem3zE7gn2sQ49fpn3s+p+NMQtF5DmcFxxn3GiM+dCzzhrUgLkiw/78102m88bP51BP1v+6E9Jo5rLxW1RLeCh6b3jGGPOOZ/4a9D/u0HnDYZkzDnt9B/19lirGesYsFYcxpg19470A1T09aYxJl7epE3pj9fI54AMyG2NPovqkdAJ2lwY0BOrtlXgiWUJFeeAaXd6Hwb7ABN9y/wIGicinsmzL9Z7U+6bPAvYQEf82P48+bOcGbOuf0dxof0KPwS0B1wP4uOPdBLS3IeqJeCHTCh5mAZ9wQ2rO+j3QMOksZ9ISVDN2rm/dz9DxZXWxM95laDnC+KN9yxVyzP28i3rtPuedKNpzdy9SvyNvnDDwK8A5ItLJs4+90N/m38c9wETRZMXjaR+iBPXO7gd8YIyZnWbIN09eA57z3cHfOSQjxpgZ6Dl7A3qd+M9D16vcK0O7XWMsievbUm0UW7RmBzsUMuAT8Hum90ffcg36BuuffxewETW8Lkd7Z93lLH+hZ7mL8InngV+gBtltaE+y09FUBZ915p/qrHMvqpW6FH2rXoJ6EdztTCRPAT9wAPpgegLVKF3oHIsFtBeT1wHPoyGV76Gat7PQB9F+zjJnONu/BjgSOMyZ3g0VbK9w2nAa+iA2wJQ0x2hklt9wg7PMmyH+W4MaSs+hvdY+C8xDPRh9PMs1As+mWX8MsA14CTWuzka1TNuAgz3LfcnZ153Of3eZ819toL2Avzfq+XjF+d8/42xvcYTHfBgeAb8zbYrnfDrNae9Hzn/TPcC1YIBrchzr01ADYzpqrJ7nbH8VsLtv2d7otbUE2Ar08M3vhRo976JhwROc4/UN2ncuyXne+Lb7Z2d/U9Fz/ha0c0i745VjG191ll8FdEkz/0+op/AHzrlwMnAJqvvaN8/r+6Qo73l2qMyh6A2wgx0KGTI9gJx5Tzg3ydo08+5ybp5HAy87D5fFeHqnOcu5D4xhvun/DbyJejjWogbBeM/8r6Dd9bc52z/JWabRs4x7s56Y4zem7U2JenPmOW1/G/We3IXHMHCW645m71+MCpSXoz0NBzjzOwG/Qb1dbXh6DKKaqnvQ0NsO5zdPznCMshlj451lLgvx3xq08sB3nf9qOyrY9veabCSNMebMOxL1Zm4GtqA6uyPSLHeFc3y2o6G/Y5xz6y7fcsc4/+dW1FiZHOUxJ40x5kyfjIbcd6DG6D3A4CDXAgGMMWe501CP4zbUEH0YGJVh2fuc7f4pw/zd0NQPC53fv9L5774W5rzxbbMf2sFmnTP8EdVqhTHGBjvL/zzD/BrnXHjDORc2OJ+vQz1m+Vzf1hizQ87B7b5vsVQUjn7jA+CXxpgfpJl/F3qT3CPptlUjInIt+pDb3WiPxiDrGOBaY8z3Y22cpWoQkUtQCcO+xpj3i90ei8XFCvgtFYVotvdR6IO/BhXtWoqEiIwj9X/cFtQQs1iiRERGoylWfgT83RpillLDGmOWSuMTqPbnA1T7ZbuVF5eH0J6KT5AqAWSxJM1vUUnC86hG1GIpKWyY0mKxWCwWi6WI2NQWFovFYrFYLEWkbMOU/fr1M8OGDSt2MwKxZcsWunXrlntBS8HYY50c9lgngz3OyWGPdTJU63F+5ZVXVhtj+qebV7bG2LBhw5g9O4rk0/HT2NjIxIkTi92MqsAe6+SwxzoZ7HFODnusk6Faj7OILM40z4YpLRaLxWKxWIqINcYsFovFYrFYiog1xiwWi8VisViKiDXGLLHQ1ARTp0LPnlBTo+OpU3W6xWKxWCyWFNYYs0TO9OkwZgxMmwabNoExOp42TadPn17sFlosFovFUjpYY8wSKU1NMGkSbN0Kzc3t5zU36/RJk6yHzGKxWCwWF2uMWSLl+us7GmF+mpvhxhuTaY/FYrFYLKWONcYskXLvvcGMsXvuSaY9FovFYrGUOtYYs0TK5s3RLmexWCwWS6VjjTFLpHTvHu1yFovFYrFUOtYYs0TK5MlQV5d9mbo6uOCCZNpjsVgsFkupk6gxJiJdReQlEXlDRN4WkR850/uIyL9FZL4z3i3Jdlmi46qrghljV16ZTHssFovFYil1kvaM7QBONMYcDIwFThORo4BvA08ZY/YBnnK+W8qQESPg/vuhoQE6dWo/r65Op99/vy5nsVgsFoslYWPMKK50u84ZDHAGcLcz/W7gzCTbZYmW00+HN9+EY45JTauvhylTdPrppxevbRaLxWKxlBpijEl2hyKdgFeAkcBvjDHfEpH1xpjenmXWGWM6hCpFZAowBWDgwIGH/uUvf0mq2QWxefNmulehYv2vf92TW25RF9hXv/oeZ521LPZ9VuuxLgb2WHdk6dKu/O1ve/LkkwPZtq0T9fWtnHTSCs4990OGDNme1zbtcU4Oe6yToVqP8wknnPCKMeawdPNqk26MMaYVGCsivYGHROTAEOveBtwGcNhhh5mJEyfG08iIaWxspFzaGiWNjanPvXrty8SJ+yawz+o81sXAHuv2TJ+u3t/m5lSuva1ba5k+fQhPPjmE++/Pzytsj3Ny2GOdDPY4d6RovSmNMeuBRuA0YIWIDAZwxiuL1S5LdHhzia1YUbx2WCxxY8uAWSyWQki6N2V/xyOGiNQDJwHzgEeAC53FLgQeTrJdlnjwGmMffVS8dlgscWPLgFkslkJI2jM2GJgpIm8CLwP/Nsb8E/gZcLKIzAdOdr5bypwtW1KfrWfMUsnYMmAWi6UQku5N+aYxZpwxZowx5kBjzI+d6WuMMR8zxuzjjNcm2S5LPFjPmKVasGXALFHR1ARTp0LPnlBTo+OpU22Iu9KxGfgtseHXjCXccddiSQxbBswSBdOnw5gxMG0abNqk98xNm/T7mDE631KZWGPMEhteY2z7dti4sXhtsVjixJYBsxSK7QRS3VhjzBIb/pCM1Y1ZKhVbBsxSKLYTSHVjjTFLbLgC/kGDdGyNMUulkq0MWG2tLQNmyY3tBFLdWGPMEhuuZ8x9AFkRv6WSSVcGDODww20ZMEtubCeQ6sYaY5bYcG8aw4fr2HrGLJXOiBFw4on6ea+9dNy3r/WIWXJjO4FUN9YYs8RCayts2wYisPfeOs16xizVwKpVOv7MZ3T8zDN6PVgs2bCdQKoba4xZYsHVi3XrBoMH62frGbNUAyudYm6HHqresQ0b4K23itsmS+ljO4FUN9YYs8SCG6Ls3j0l4LeeMUs14HrGBgyA447Tz08/Xbz2WMoDbycQP3V1thNIpWONMUssuJ6x7t1h4ED9bD1jlmrANcb694fjj9fPs2YVrz2W8uH00zsa7nV1MGWK7QRS6dQWuwGWysT1jHXrZlNbWKoLN0zZv397z5gxqqG0WLKxaZOORfScOf54uPnm4rbJEj/WM2aJBW+Y0vWMffSRLYlkqWza2mD1av3crx+MHKkvI6tWwbvvFrdtlvLgzTd1fOSROl66tHhtsSSHNcYsseA1xhoadLxzp4qZLZZKZe1aNch694bOndW7YXVjljC88YaO3ZCkNcaqA2uMWWLBa4yBFfFbqgOvXszFGmOWMLiesYkTob5ea/q6oUtL5WKNMUsseAX8YEX8lurA25PSxSvit2F6SzZaWlJpUMaMgSFD9LP1jlU+1hizxIJXwA/WM2apDrzifZfRo6FPH1iyBBYtKkqzLGXC++/D9u0wdKiGuq0xVj1YY8wSC/4wpfWMWaqBdJ6xmho49lj9bEOVlmy4erGDD9bxHnvoeMmS4rTHkhzWGLPEQibNmDXGLJVMOs8YWN2YJRiuXmzMGB1bz1j1YI0xSyxk8ozZMKWlkkkn4AdrjFmC4feMWWOserDGmCUWbJjSUo2kC1MCjB0LPXqoJmjZsuTbZSkP/J4xG6asHqwxZokFf29KK+C3VAOZwpS1tTBhgn623jFLOtauhQ8/1HQWI0fqNOsZqx6sMWaJBX9vSusZs1QDmTxjYEOVluzMmaPjAw+ETp30szXGqgdrjFliIVuY0uZaslQqmTxjYI0xS3b8ejHQiEJNjd43m5uL0y5LMlhjzBILfmOsvh569tQbyrp1xWuXxRIXbW2wZo1+7tev4/zDDoOuXeHtt1P1Ky0WF79eDDS8PWiQvsAuX16cdlmSwRpjlljwG2Ng01tYKhu3LuVuu0FdXcf5XbrAUUfp52efTbZtltInnWcMbKiyWrDGmCUW/AJ+sOktLJVNthClixuqnDUr/vZYyofW1lQZpIMOaj/PGmPVQW2QhUSkBjgVOAk4AhgEdAXWAu8BzwEPGmM+iKmdljLDL+AHK+K3VDbZxPsubp1KqxuzeJk/P1UGabfd2s+z6S2qg6yeMRHpISI/BJYAfwdOAN53Pt8NPAPUA98GFojIv0Tk2HibbCkHsoUprWfMUolkSvjq5aijVAf0+uuwYUMy7bKUPun0Yi7F9ow1NcHUqar5ranR8dSpOt0SHbk8YwuB14BvAo8YYzZlWlBExgHnAg+KyI+MMTdH10xLObFzpwr1a2uhc+fUdOsZs1QyQcKUDQ1w+OHwwgvw/PNw+unJtM1S2mTSi0FxjbHp02HSJL2fu705N22CadPg7rvh/vvtORwVuTRjpxpjTjbG/DGbIQZgjHnNGPMdYC/gychaaCk7vF4xkdR06xmzVDJBwpRgU1xYOpLNM1asMGVTkxpiW7d2TKvR3KzTJ02yHrKoyGqMGWNeCbtBY8xWY8y8/JtkKXfShSjBesYslU0QzxhYY8zSkVL0jF1/fe7cZs3NcOONybSn0gncm1JEBojI3p7vIiJTROSXIvKpeJpnKUfS9aQE6xmzlC5R6GKCesZ2313Hzz9vNTjlQpy6qXXrOpZB8uI1xpJMmH3vvcGMsXvuSaY9lU6Y1BZ3AVd6vv8I+C1wGvCQiFyUawMisqeIzBSRuSLytohc4Uy/RkSWisjrzvDxEO2ylBjpelKC9YxZSpPp0zU8NG2a6mGMSelixozR+UEIIuCfPj1VoxLy35clOaI6PzLhhii9ZZC8dOsGvXrBjh2ppMJJ4N7Ho1rOkp0wxtghwAzYleriUuC7xpj9gGuBrwXYRgtwlTFmf+Ao4DIRGe3Mu9EYM9YZHgvRrrKlUnup5ApTrlypyTEtlmITpS4mV5jSuy8/VoNTmiShm8qmF3NxdWNJhir99+9Cl7NkJ4wx1gtw7fJDgT7AH53vM4A0Dtb2GGOWG2NedT5vAuYCQ0K0oWKI+22rmGQyxrp0gd69oaVFs5VbLMUmSl1MrjBlmH15X9ROPPH4inlRKzeS0E1l04u5FEM3Nnly+koSXurq4IILkmlPpRMo6avDEmA0mlvsE8A8Y4x7avQCtofZsYgMA8YBLwITgMtF5AvAbNR71qGCoYhMAaYADBw4kMbGxjC7LBqbN29u19alS7ty8cWHs317R5+024X47LNbmTbtZYYMCXVYS4KXXhoAjGbr1pU0Nr7Tbl6PHkewfn0D//jHS+y9dxoXQYH4j7UlPirhWN999zE0N2e/DTY3w513tjBpUuYaRq2tsHr18YDw1luzmDevo7gn6L6mTWvljjuguVloba0BhE2b4Lbb2rjzTsM117zNkUfat5k48J/TUZ0f2XjuuUOAnrS2vkZjY/rkczU1o4DBzJjxLg0NyRSpnDChK3feeTjNzWlipw6dOrUyfvzLNDaGe05Vwr0jcowxgQbgO8AG4D5gK3CFZ95PgGdCbKs78ApwtvN9INAJ9dRdC/w+1zYOPfRQUy7MnDmz3fdLLzWmrs4Y9YelH+rqjLnssuK0t1Buv11/w5e+1HHe8cfrvCefjGff/mNtiY9KONYi2a9Dd6ipyb6dlSt1uT59Ct9XrqGhwZj334/2OFgU/zkd1fmRiZYWY+rrdRtr12Ze7gc/0GV++MP89pMvjz2m51ttbcfnU0ODzs+HSrh35AMw22SwaQKHKY0xPwW+AnzkjG/yzO4DTAuyHRGpAx4A/miMedDZ9gpjTKsxpg24HS25VLFUei+VTAJ+KH8Rf6Xq/KqVqHQxQcT7UWlrbDqB5IhbN/X++7BtG+y5Z8cySF6Kld7i9NNV0/aJT6SmicCUKTrdJnyNjlCFwo0xfzDGfMUYc4dj5bnT/9sYc3eu9UVEgDuAucaYGzzTB3sWOwt4K0y7yo1K76WSSTMG5Z3eopJ1ftVKVLqYIDnGguwrCFG/qNkXjMzErZsKoheD4mbhHzECPvvZ1HdjVEs3YkTybalkQhljACLSRUSGi8ho/xBg9QnABcCJvjQW14nIHBF5E61/eWXWrZQ5ld5LJZsxVq6eMZuNujK56qpgD9src9yRguQYC7KvoET1omZfMLIT1fmRiSA9KaH4xcL99+sPPyxOOyqZMElfdxeRf6J6sfnAHM/wljPOijHmWWOMGGPGGE8aC2PMBcaYg5zpnzbGJKNQLBKV3ksliGes3Iwxm426MhkxQuvrNTS0L90Feg02NOj8XF6AIJ4x777817+7r/r6YO2O4kXNvmDkZsQI+P3v088Lc35kohw8Y9Dxfr14cXHaUcmE8YxNAw4Dvo4mej3RM5zgjC0BiPttq9gE8YyVW5iy0nV+1Yyri+nRIzWttjacLiZo9n13X1OmtA8Luvu66KLkXtTsC0Yw3P926NCUsRz2/MiE6xnLZYz16wedO2u2/nR56uLGNcZqHIvBGmPRE8YYmwB81RjzK2PMv40xs/xDXI2sNKJ6Gy9VMpVDgvINU1a6zq/aGTo0dd6Cho1uvjn4NRhEwO8yYoRue8MGTYmxYUNqX0m+qNkXjGDcdZeOr7sOZjlPudGjw50f6Vi3Dj74IHMZJC8iqTJaxfCOuffrAw/UsTXGoieMMbYS2BZXQ6oN9w25T5/20//rv8q/l0q23pTlKuCvdJ1flJSjIHzxYjWM3P9v0aJw6wctEp6LbKHM2tpoX9TsC0Zu5syBV17RckRnnAH77KPT33+/8DqRucog+SlGFn4X1xg7wslzYI2x6AljjP0Q+JaI9IyrMdXGgAH6dlRbq25ogO98p3w9Yi7ZwpRuGGfVKn34lQuVrvOLinIVhL//vo4PP1w9FWvXwsaNwdcPGqYMgj+UCfrUP+qoaF/U7AtGbu52cgScdx507aoVRPr101DhsmWFbTuoeN+lmLox1xg78kgdh31ZseQmjDF2NjAUWCwi/xKRv/mGv8bUxorlxRe1RuO4cTBsmE4rN49ROrIZY507qzewtTXZoreFUuk6vygoZ0G4a4zts0/qWgzzwInKM+biDWVeccV8QNsV5YuafcHITnOzhnJBtXwuXu9YIQQV77sUyxgzJnV+H364jq1nLHrCGGP9gCbgdaAO6O8bIngnrC6ef17HRx8Ng51Ma8sroB9pNmMMylM3NmIE3HJL+nki5a/zi4JyFoS7D9aRI1PG2MKFwdeP0jPmZ//9NwHw8svRbte+YGTniSf0HrXffqnwHKSMsfnzC9t+WM9YsdJbrF8PO3fq/XzUqFQbyimyUQ6EycB/Qq4hzoZWIq4xNmFC+Wqp0pFNwA/lm95ijpO8ZcQIDR+5nS969NC33HLW+UVBOQvCvcbY3nvr56CeMa+Xt2/fyJvG8OGb6dwZ3n1XPWVRESTVRjW/YLjC/Ysuat/RyhXbF2KMtbbCW05q81IPU7r36UGDNFQ7aBC0tBQeprW0J3TSVxenrJElT9ra4IUX9PP48SkDpZI8Y+kE/FCe6S02boRbb9XPf/6zPhSbm/Xhu3Gj/p/VTjkLwgvxjK1Zo6GcPn1U/xk1dXVmVyhr9uxot+3q0y65pP30ai93s2YNPPKIdkCZPLn9vCg8Y0HLIHkptjHm3rf32kvHNlQZLaGMMRE5WkSmi8gmYLuIbBKRx0RkfEztq1jeeUcf4kOHqvvZDVOWk4GSDmOCG2Pl5Bm77Tb9v44/PqWb6NQJTjlFP5eqMD1JylUQ3toKCxbo5+HDw3vG4gxRurjnXNShSlDP1y9+0X7az35WvR4x0Beu5ma9vl0jyKUQY8ztaTxunH5fvjx4T+NihSndZ5J733ZfVqwxFi1hMvCfDDQCewA/B6Y64z2ARhE5KY4GVipevRhUTphy2zY1yLp2zewlKLffunMn/PKX+vmb38Ei618AACAASURBVGw/z/UcPP54sm0qRYIIwmtrS08Q/uGH+uDdfXd9gQjrGYtavJ+OOI0xUF2Ql2ovd+MNUfrxCvjDeMS9PY23OUmiWlqC9zT2vrAnqdfK5BmzPSqjJYxn7FrgEWCMMebHxphbnfEY4J/AT2JpYYXy3HM6do2xShHw5xLvQ/l5xv7yFw0NjB7dMWzjesYaG1M32GoliCC8pQVOOKG0cpF5Q5TQ3jMWJJdUEp4xV0BujbH48ecW89Ozp/7X27cH101F0dO4c2fdb2trsvdOG6ZMhjDG2EHA7cakvT3d5sy3BMQr3ofy8xZlIogxVk6/1ZhUCOcb30iVA3EZOBAOPVRvzI2NiTevpHAF4ek8onV1qcSWF1ygiS5LJReZ3xjr00fP340bNQ9gLsJk38+XUaO0TR9+GM+D2P87i1WQuhTw5xZLR9hQZVQ9jd2QaZL/jzXGkiGMMbYeyKQiGOnMtwRg5Up9ADQ0pHrSeA2UQjM7F5NcPSmhvDxjTzyhb8qDB8PnP59+GddbVqq6sSS9UKefrqkAQM9vb+3Ft9+GT3xCPYjbt5dOLjK/MSYSTjeWRJiyUyc1+iEe75j1jCktLZI2t5ifsD0qo+ppXIws/NYYS4Ywxth9wE9FZLKIdAUQka4iMhkNYf4tjgZWIm4vyiOPTHkR3OzOzc2a/btcySXeh/JKbfHzn+v4iiugS5f0y5SybizpjPibN8PcuWo8rFjRvvbiqFHae8xfj9VP0rnI/MYYhNONJRGmhJRu7KWXot+2NcaUl1/ukza3mJ+wnrGoehoXo0dlNmOsnB0HpUYYY+xbqDbsbmCLiGwAtjjf/+nMtwTAL953KafwXSaChCldD0Kpl0R69VWYMUN/y5e/nHm5I45QQ3r+/NLKLl+MjPj/+Y/+p+PGpT8H/vjH3DfwpHORpTPGwnjGkghTQrwiftcY23NPHVeTMeb1HH/3u1oJu3fvVA/bdIQ1xqLqaVwKYcoePTSUv317yitsKZwwSV+3GWPOBw4ALkK9YRcBBxhjJhtjtsfSwgoklzFWziL+IMZYXZ3Wd2trg9Wrk2lXPrhesSlT9Oacidra0kxxUYyM+M88o+Njj00/v9RykbW1pYxRbyqHMCWR3AdSUp6xl1+O3iPhasYOcpS/1aIZ83uOQd22s2dn9xyHLYkUVemppMOUxnQ0xsCGKuMgTGqL40SkuzFmnjHmHmPMdc54noh0E5Hj4mxopbBjR+rN9qij2s+rhFxjQYwxKM3Er35t1V/+oiG1M8/MvW4phiqLkRE/lzFWarnIli3TN/wBA9yi3Eo+Ycq4PWPDhulLzJo10acVcD1jrjFWDZ6xbJ7jlpbsnmPXi9rUFCy9RVSlp5IOU27cqM+shob216Q1xqInTJhyJjA6w7z9nPmWHLz2mp7co0erq9dLJXjGggj4ofRE/Om0VS6nnZbb43XqqTqeMUMf7qVA0l6onTs1TAlwzDHplym14tTpQpRQegJ+0BeDuEKVrjE2bJg+eDdtirb0UilSiOe4Rw+9h23fHsyL6PY0rq/vOC9M6amkjbF0XjGwxlgchDHGssluuwNbC2xLVZApRAmV5RnLJuCH0tLHZXtDNiaYtmrwYBg7VnsKPv10vO0NStJeqFdf1d+/336ZDZNSK06dyRjzesayhQRbW7XDjUg8dSn9xGWMuWHK3XZL6cYqPVRZqOc4rG7s9NPhd7/Tz506te9pHLT0lDcLfxLieWuMJUdWY8wJTf5QRH7oTLrY/e4ZfgLcBMyJvbUVQDZjrJQMlHwJG6YsBc9YVNqqUktxEcQLBWo8udqpQlJguCHKTF4xyF6culOn5ItTuw9SvzHWu7cOW7dm1zXGXZfST1w9Kl3PmNcYq/RQZaGe43zKIrlRj0svbd/TOOj53rOn3lu3bk3Gc2mNseTI5Rk7EviKMxjgHM93d7gQWA1cHl8zKwNjOmbe91IJYcqgxlgppbeISltVarqxIF4o0Af7sccWngLj2Wd1nEkv5uIWp54yRR8ubqqLvn2TL06dyTMGwXRjSYn3XVxj7JVXou2J7BpjvXtXjzFWqOc4H2PMLfR+2GHB1/GTZKgykzEWpoOLJRhZjTFjzM+NMf2NMf2BD4CJ7nfPMMQY8zFjzKvJNLl8WbxYvV59+sC++3acX0lhynIS8EelrRo/XkuozJtXGjepESPgrLPSz3N1Ktdco//V888XlgKjrS24Mea27eab9e1+yxY9bitXJteL0sU1xtwHq5cgurGkxPsuAwfC0KF6zObNi267XmPMDYVVujFWqH4xbI9KUCMaUgl88yHJ9BbWM5YcYVJb7G2MeSPOxlQ6Xq9YusSX1jNWHKLSVtXWwkkn6edSCFUuXAgPPKCfzzqrffjR1alcfTV86lO5t5UrTDt3rmqnhgxJvTUHpb5eH4wAd9wRbt1CMCb1IE0XJgrjGUvKGIN4dGOuZszrGat0zVih+sWwWfjdXrANDakqFfmQZHqLTMZYnz6qC964sWPCYEt+5NKMXSAincJsUERGikiAd+Pqw1+P0k+fPnrxr19fOj3ywhK2N2UpeMai7OFXSrqxK6/U8+j88+HBB9ULlU6n8s9/5t5WrjCtN6VFrgz76bj4Yh3fc09yBdc/+ki9fn37qlbKTxjPWFJhSojeGDOmOsOUrn4xXf3JID0cvektgoSMXa/Y2LGF6QtLIUwpYr1jUZPLM3YV0CQi/yMiB2daSET6isj5IvIP4DVgcJSNrBSyifdBvRalJGzPh6C9KUvpd0bZw++003Q8Y4amMCkW06fDww+rUewmr81EFGHaXPnFcjF2rIZu1q+Hhx7KbxthyaYXg2CesaTDlBC9iH/zZg0zNzRA587VE6YEfXn6xjf0c10diJjAPRy7d1dpyc6dwY5VFHoxKI0wJVhjLGpyacbGomWOTgBeE5GNIvKiiDwqIg+KyAwRWQisBH4FNAH7GWNsnUofmzbpxV1bm/1iLPdQZdAwZf/++na1erUmWCwm3h5+fsLkAAK9UY4Zox5CV0OVNDt2wFe/qp+vuSalRcxEFGHaMHqxTLjesWnT8t9GGIIaY9k8Y0kL+CGlN3rjjWgMfq9XDNp7xqqh9qBrbP/61zBjxqxQPRzDiPij0ItBccKU7nPJizXGoiWnZswY81djzDHAPsA3gdeBFqAbsAKtTXkaMNgY8zVjTIIlTMuHl17St89x49I/9F3KXcQf1BirrdVs4sakvAvF5PTTNSFvJycoLxI+B5CL6x0rVqjyhhvU0Nh//5RRlo1Cw7QffKBD795wwAHh2+ty3nmqH5s5M5woOl/CGGOZsqwXwzPWq5dqjpqb9dwsFG+OMXf7PXpoCLca9EBuuNf1OIYhjIjf9YwVaoyVQpgSbI/KqAkj4G8yxtxqjPmyMeYMY8ypxpjzjDHXGGP+bYzJkRyguskVonSpFs8YlJaIH9QwbG1Vz0BbW/gcQC5J6sb8ucF69IAf/EDn/frXwVJbFBqmdUOUEyZoG/KlVy8491z9/Pvf57+doOQyxnr0UD3Zjh2Zz9FiCPghWt2Y3zMG1ROqXL8e3nsPunRJlYIKQ1DP2OrV+sJSqHgfkjPGNm9Wg7xrV70W/FjPWLQUcOu0hMHtSZlJvO9S7p6xoAJ+KC0RP2jYB+DgjOrIYAwcqMbLO+/klzw1KOlKOG3erAZlp06qZQlCtkSsIrnDtIXqxby4oco774w/fJ3LGIPcIv5iCPghfmMsThF/ocmFo8T1Vo0bF+zFxU/QHpVuiHLcuJT3PV8GDNDIwqpV2cPUhR5nr1csXacca4xFizXGEqCtDV54QT+PH5992XLPwh9UwA+l5xmLwhibPl01ga4hkU/y1CBkK+EEapDlyg3mJVMi1poamDUre5g2SmNswgQYNUrP/8ceK3x7mfCmtchmjOUS8RcjTAnJGWNRi8TTvUDEcX0EpZAQJQT3jEUl3gc15tyX9mXL0i8TxXHOFqIEa4xFjTXGEmDx4m5s3KjJGl33fybiDFPG/Uba2qppCUTSF8T1U6qesTFj8lvfayD5hc9Bk6cGJaoSTl68iVjb2rT4eWsrPPlk5nXWrFEPYNeu0TxoRJIR8q9erTmSevXKXlMym2espUV/f1J1Kb246RHeeUcfsoXg14xBPJ6xbC8QUV8fQXF7pOZrjLmG/IIF2T25UYn3XbKFKqM6zrmMsUGDtPftqlW6TUthJGqMicieIjJTROaKyNsicoUzvY+I/FtE5jvjNFl/ype33uoJ5NaLQXxhyiTeSN0QZbduwbRDpZTeAgr3jMVhIGUiqhJO2XA1YjffnHlfbvj9yCP1xhwFX/iCGhqPPhqfLsbrFcuWFy2bZ2zNGh337Vt46CksXbuqxskYLdBeCElpxpK8PoLiesaOOCK/9Rsa1DBqblZNWCai9IxB9vQWUR3nXMZYTU3KaM/22y3BSNoz1gJcZYzZHzgKuExERgPfBp4yxuwDPOV8L2u8XqgbbtDaR0uW5H4bicMzltQbaRjxPpRWSHbtWv1/6uuzh62ykYSB5BJVCadsnHKK9shcuhTuuy/9MlGGKF0GDIAzzlDv3N13R7ddL0FClJDdM1Ys8b6La0AUGqpMSjOW5PURhGXL9Nzu2TN9Oayg5OpRuXKlHsdu3dKXwcuHbOktojrO7n05kzEGtkdllCRqjBljlrs1LI0xm4C5wBDgDDRFBs74zCTbFTV+LxToq/d//pPbC+XVUWXqTh+WpN5IwxpjpeQZc1MEHHRQ/l6OJAwkl6hKOGVDBL72Nf18443pc07FYYxBKlR5xx3RXQdeghpj2TxjxRLvu0SlG3ONsXRhyig1Y0leH0Fwj9thhxXWCziXiD9K8b5LtjBlVMc5l2cMrG4sSgIXZRCROuAK4GxgD6BDEQljTODbkogMA8YBLwIDjTHLnW0sF5G02xGRKcAUgIEDB9LY2Bh0d4mxdGlXLr74cLZv73jVtbTocPbZrUyb9jJDhqSvedS9+wQ2b67jH/94ll69Cu9Sdvfdx9DcnP2vbm6GO+9sYdKk/DOVzp/fHTgM2ERj4ys5l//ww27A4TQ1baGxMZraLps3b87rvLj//iHAPvTvv4zGxvfy2nd9/TFs3Zr7kuratYXGxsIywp5wwj48+uhgWlszP0U6dWrjhBOW0diYf9KuoUNr6NlzPLNn13Hzza9y0EEbd81bvXobs2e3UVMjNDc/S2NjgJowAamrgz59xrNgQRcaGlrZubOG+vpWTjppBeee+2HGaycozz23PzCQlpZ5NDZmds3u2FEDHMfixW089dTT7R6mjY39gQOAlTQ2vlNQe7KR+ZzW6+eZZ7bR2Phi3tufP/9AoB9LlrxFY+NqALZu7QQcy+LFrcyc+UxeJa78JHl9BOH++/cG9mLQoMU0Nqq1nc/9o6ZmT2AEM2cu4YADOl5r99+/F7A3gwZ9SGNjNIK4DRsGAKN59dWO515Ux/nttw8A+rN69ds0NqZPBtnWpr/t6acXM2pUllIVPvK9T1c0xphAA3Az0Aw8BFwLXO0fQmyrO/AKcLbzfb1v/rpc2zj00ENNKXLppcbU1RmjfoT0Q12dMZddlnkb++2ny82ZE02bRLK3xx1qagrbz9NP63YmTAi2/Ecf6fJ9+xa2Xy8zZ87Ma70vflHb8utf57/vKP77oLz/vjENDdn31dCgyxXKd7+r25s0qf30669/zYAxhxxS+D78PPZY+mNZV6e/67HHCtv+EUfo9p55JveyAwfqsh980H76TTfp9KlTC2tLLjKd083NxtTXaxtWrsx/+8cfr9t46qn203v10umrVuW/bS+XXpr7XhTV9RGEU07RfT7wQGpaPvePBx/U7Xz84+nnn3GGzr/nnvzamY5Zs3SbRx/dcV5U96Gjj9ZlGxszL3PXXbrM5z8frv353qfLHWC2yWDThHHOngN82xhzljHme8aYH/mHIBtxPGwPAH80xjzoTF4hIoOd+YPR8kplSRTx+qhF/EmEtCB8mLJfPw0PrFmT+5jFTRRpLaKscZkLb24wf4glbAmnXEydqoL6Bx9srw2ZM6cXEH2I0tU4pjsnotI4Bg1TQmZdTLHSWrgsXqy9QUFDSfn2jk6nGYPodWOuKZCNqK6PIG0pNK2FS670Fm6YMirxPmQPU151VfZOKRDsONswZbKEMcYEKKj4hogIcAcw1xhzg2fWI8CFzucLgYcL2U8xiSJeH7WIv9ByN0EJa4x16pR6kK0sovnd0gJvv62f801rAdmTp0ZtIEEqN9jQofq9kBJO2RgyBD77WdVu3XxzanpcxljcGse1a3Xo1i37g8bFFfH7dWPFFPC7ulS3DYX0jk6X2gKi1Y3deivccou+OHTpksz1kY2mJv3dAwfmTjeUC7e9Cxd2TG+xYoUev+7doxPvA+y+u46XLeuoqezVK9Wz2a9RC3OcwxhjVsBfOGGMsduB8wrc3wTgAuBEEXndGT4O/Aw4WUTmAyc738uSKLxQUXvGkvLYhDXGoDQSv773nmayHjYs5WnIF3/yVJcLLojWQHIZMSL1AHjnnfxLOOXCPTduv10f+mrA6sE65pho9xV3rzvXc5QrrYVLLs9Y0gJ+b+9o/4M4H89hJs9YVOktHn1UPXagRtnbb+v14SaGrqmJ/gUiF96UFoXq4err1XBtaenoIXK9YoccUlgngXT77NtX/29/bd/vflfvxePHw5e/3P4+9MUvBjvOW7fqNjp37nheeNljD/1dy5YFr/hhyUCm+KV/AL4KLAJmAt8FpvqGS4NuK4ohEc3Y1VcHE1uBMZdcYowJGK9nh7mMX7efePXVu3Z73XU66evDHwy+/1tv7dj+Qw7ZNf8xTjP1bEm7ale2mseu/k/H9QcPDr7/2bPNr52ftEuLEGC9U3jcgDGP/mF1+30vXRp837BrtV1ahNmzA6/7p96XGjDm05/27P+RR4LvO51o6tZbjQFzGC8ZMOZ5jsq8/ic/mfe5t4L+Bozp3t2Y1lbP+pdcErz9nnNvF5/8ZIfljuFpA8bcxOXmJQ4zYMy+zMt57uUcHnmk3aqBNY607Dr3OpDt/+ZzBoz5DPfptKVL26/rO/du5RIDxlzE79tNP5ZZBozpIH8Jce6ZwYM7tj3HuXcpN5s6dmTdbNp7TJqh5ROf3vW1paX9ufdjvm/AmG/zk8zbcO57xqg+8dJLjelRt9UIraYH680k/mq6Oved7/PjduvupNZ0q9vR8S9Ic+5lHPI8977GDQaM+fHn57ZbdebMmaHve8YYc8IJ+nX69PbnnnsMr+T69OvnOPeyDWN43YAxr7ySWv0/d71jhFZTxw4zl1G7lt2PdwwY8zpjAp17CxhmwJg9+CBzG5z73p576temJrPrvhdoKOC+5z/3dlHofS9miEgz9ktgKHA88L+ooN8/VD2BvFA0cyWZYyy7wpTb+0TWrtN5nOu5CoBamqmhlS5oj7SRzOfUcYXHCcOUQnIZiLrEVqzO/7WxieG7crqdeOLxqp352Z40MTzQ+m80jwYKr0mZjpGoOGk+BSQyysKrHAJoRvYo37zT4Z6z3+d/ORbNabGQvZn6pwmRZk0P7F0mvxwI7n/h/je5GMYiABYxrN30lahLLGnP2L1MppnsGXab6cw95NYdbGzWUhk9e3YMae2JusSWkDuO1y6dT3M9hho20Yv7OYftNPAxnuTH/LDdOnW0cMxQzRY6a1bOXUTKy6hQ7PB91keyvUy6sdmoUOwwZkeyHy9DUMGYqxtrbYXL/m8ohhq+zg3sx7u7lh3lfH6XUYG2vQKNTbr352wE1Y01MZyp3ExP1lNDKz0f/2vRapKWIoFv38aYmhxDwjmoS5N2uiFpLyCoYycNbOF+JjGCBRm3sStMuSOLfzgP3Avs69xAK7WsoS+7s5S3GMMfZxUonCB8mLKJ4czhQAD+6+u98xIgT+c0xvCmp7KAqHbm7/0Yw5tM57Sc23ijJX5j7H3yzCSbA9cYi6rMSjY6swOhjY30YoeT2aaZzkx7dr9I6wpOnqzXSjbq2MkF5BendP+LoMbY3iwEOhpjq1CxWNKasc30CLhc7gtxXbMu49eLAeyBisU+ZM+s28heJ1VjgC8wngVpXo5OGKbHNsksBy102nXdxG2MvYJemIeSO9VPWFxjzNX03XYbvDK3G3vwId/nf9stux/zAJjHfoG2HbUxtus+zSVsopca6y0NRatJWpJkcpmV+lCqqS1c3n9fw3U9exoj0mZ69tTvQVINzJmjXtT99ou2TWedpdu9997UNLdr8pAhxmzZUtj2r7hCt3XDDbmXfewxTVFQU9PecxwmdUFU6R3cqMT8+cF+Zxjc43veedFv2xhjzj5bt3/33fFs3yXJVBpx72v8eN1G0N7127dr6LRTJ00nYYyOQafvCu/FhD8NQI8ewaIwPXvm3vYrr+iyBx/ccd68eTpv+PDs2ygklcJ//qPz9903d1uj4vXXM/+ufFMu/P3vus3TTktNW75cp/Xo4ZMQRMQ11+j2v/c9TW3Su7d+v+++jsv+3omwB01B4UYbv/jF3Mu6qW+uuSb9/CTvHaUOEYUpEZHeIvItEfmHiDznjP+fiETrwqkAvEWXZ8yYFUpYHVd9ynTpGy64QMWlS5fCL35R2PaDesaiEiBH0etu1SrttdqtGwwPFtUMhZs6IVOplEJx6xLG7RlLsq5gtl6pLuPG5f9/hUlrAdr7b/fdNQzkeiFWa25U+vVLvi5llL2jM4n3ISXgX7IkexWEQjpcHHKI3i/ee09F4EkQVUoLL+lKIsUl3ge9N86cqZ+vvVZ7PK9fDxMmwGc+03H5UU508t13O85LR5CelC65elSWYk3SUiTwKSIiI4A5wI/R1M8fOOMfA2868y0RsNtuejNdvx62F5ZofBebNsGCBdo7ZpRHNlBTAzc4SUb+7/8KuyG6hcJzGWNRXZxR9LrzlkGKQ3MVpzG2Zo3eAOvr2/+ncZB0XUF/r9SaGh2ff76eX889B7/7XfjtbtigBnjXrqn0AEHwl0UqZo6xKHtHpyuF5NKtm07fuTNlfKajkHQ+dXWp1ChJhSoLLQ6ejuHDtVfmwoWp68QtDh71i5Krz3vWkzzf3eerr8Ljj3dcx2uMGZN7H64x5uqXs+FeG5nClKVWk7RUCfP4uRFYDww3xpxojDnPGHMiMMKZfkPWtS2BqalJvZFE5R2bM0fHo0d3vJEffzycdZZ6pL73vfz3EVTAH9XFGUVOtyiSvWZjwAA1Htat09xWUfLaazoeO1aTssZJMeoKer3Lra06vvderVcJamy8/nq4bXrTWoQxvv0Fw938XsWoSxllPjs3x1im9AVBEr8Wms5n4kQdJ2WMvfSSjqP0jHXtqseqtTV1jsSR7NUbVWhNU31s27b0UYW+fXXYvDlY/sp8PGOZjLFSq0laqoQxxiYCPzTGtMv563z/EXBChO2qeqIOVboeoExJTa+7Tm/kd90F55zT3hsRVFQfNEwZ1cUZRU63uI0xkfi8Y94wSNwkVcUhCOeeq/mTdu6EM8/UouJBz9ewIUqXUvKMQXvPYUODTqutDZ+vK1uYEoIZY4WGTV1jzA27xcm2bfpiWlMT/XXjF/HH4RkrJKqwn6Pdnzcv937CGGNu0ukPP0wfzi6le0cpE8YYM0AmdUSNM98SEVFn4c9ldIwcCZ/8pH6+/363Z2K4rN5BjbGoLs4otDNxG2MQvzGWRE/KpKo4BOXGG9VAWrwYfv/74OdrocaY6/UotjEGKc/hjBn6fezY8Al/gxpj2bLwFxo2PeQQ6NFD/5sosv1n4/XX1aN0wAHhUvAEwWuMLVum9+6ePcOfa9koJKoQRjcWxhirr1cPcXNz+udVqd07SpUwxthM4H9EZC/vROf7j4GnomxYtZO0Z6ypCZ54Iv28oKL6oMZYVBdnoQ+B5mbNWg+qGYuLuIwxV7yfhGcsybqbQVi2LPXA8Gtgsp2v+Rpj/pJIxQxT+nHbkE9JsWyaMQiWhd8bNvVnsw8SNq2tTenG4s43FkeI0sVrjMUl3i8kqhCXMQbZQ5Wldu8oVcKcJl8DugDzReQ/IvKwiLwAzAc6A1+Po4HViusZi8IYa2tLGWOZPEBRiOqDGmNRXZwjRsDdd6efJ5L7ITBvnv6m4cP1zTwu4jDG1q9XQ6NLF9UBxk3SdTdzcf31HesA+kl3vlaSZ8zFNcZWrAgmzvYShWYMUmFTN2Qatk5qUqHKOHpSuniv87i81oVEFYIaY9u3qz6ztjazke4nW4/KESPgb3/LXHaqri7Ze0epEibp6yJgP7Qs0ttAHfAOcDmwvzPfEhFRhikXLVJDadCgzA+PKET1QXtT5kpdEObB7vb+HDJEb/4i+jQyRnsVZXsIJBGihPTd3gvFFe8ffHBuwzYq/D0cRUwshcmDkO/5mq8xtuee6uFYulTrmJaSZ6xbNx127NAwbRii0Iy5dO6s94BevdRQDpPOJx8Rf1MTuypvBNW3xtGT0iWdZyxK8T4UFlUIaox5z+2gXr1cPSrdjGLdu6fuHa7h3qcPnHxysP1UMqEcqMaYncaYW4wxXzLGfNwZ32aMsSVCIybKMGUQoyMKUX2YckjpUhe4F/611wZ7sBsDv/2tfv7Vr1I53c45R6e98EL29ZMyxuLwjCUp3vdSSP68KMnnfN2yRV9uOndOhd+CUlen6xijhkkpecagvXcsDEHDlEG0XM8/r+Px48OH5saN0/tAU1Mww69d+aWAesH16zWfWZcu8cgShg/X371oEbz4ok6L2jNWSFRhxAj1di1apB0ZMhE2RAm5e1T+6lc6vvrq1L1j82Y1EFesgIcfDr6vSiXmanaWfInSM5ZLLwaFi+p37lRPRG2tPuyC4E9d4OaN+uMfg4VbnnpK3/KGDIEzzkhNnzw5tZ1sJGWMDR6sItdVq/S3RkFSyV5LlXzOV9djMnx4folavbqxUjXGwurGcoUpXWNsmfI2mgAAIABJREFU6dLsiV9Bc7+BJh4Ni1c3lss7lq38Uja9oNu7cezYeLzJXbpoz8K2Nj0/evWK/iWlELlAXZ2e+8ZkfzGM2hh7+2148klt25e+lJouop5MgN/8Jvi+KpWsxpiIrBSRcc7nVc73jEMyTa4OkvaMFSqq9+rFMmkDcnHBBZoLZ/bs9gkNM+FewF/+cvs8W6edpq7vN99MGaLpyKWjiwpveouoiuIWyzNWKgQ5Xzt1an++5huidPHqxkopTAmpB2dYYyxXmLK+XqsMNDfn9roVYoxB8FBlvvrWOEOUoNe2V8e4eTNcfnn0hbAzJUQOIhcIEqqM2hi76SYdX3hhRw/shRdqJGXmzFRnqmoll2fsN7CrUuhvAgyWiHAvhBUrcr+R5iKIZ6xQUX3YIuHpqK9PvSndkCOF8AcfwCOPaJsuuaT9vM6dNRcVZPaOrVihQ48eqRtJnLgGgL+QcD5s3Kjhlro6OPDAwrdXjgQ5X1tbVcfjaos+/3md/q9/hS9IDynP2Pz56lGqqVGjvxQoNEyZyRiDYLqxzZv1pa9Tp/yNnROcTJW5jLF89YJxivfdsOlSTxbO1tbgaYHCki4hchC5QNzG2KJF7aMaa9em/oevfrXjer16pV6Yqt07ltUYM8b8yBizzPl8jfM945BMk6uDrl31BtncXFjm9s2b9aFTV5dK+peOQnvLBRXv52LqVDWmHn44uyv9ttvUSP3MZ9KX7HBDlX/6U3pj1vUWjhkTTxkkP1Hqxtys8wcdFDwkXGnkOl9dT+nXvqa9TadNU4E7aEg9n4ek6xlzw139+iVz7gQhnzDlzp0a0uvUKft1G0Q39uKLahSMHZt//q6xY9XDs2CBvmxlIl99a1xpLbxh0zBpVopBEGPMjcaEMcZ69dJh27b2pbNuv12nnXpq5ufPZZfp+A9/0BfNaiVMbcoZIpL2cIrIviIyI7pmWSCaUGW2Mkh+vO5vV1PTrVsw93cY8X42Bg3S+oPGpESffnbs0IscUheyn6OP1ofnkiXw9NMd5yelF3OJ0hirdr2YS7Zwzdy58HUn2Y6rZ/SSz0PS9Yy5xlip6MUgvzCl1yuWTVoQxDPmivfzDVGC3nOOO04/Z/OO5aMXXL5cvVY9e8K+++bdxLSUUyFs1xjLloU/H88YdOxR2dKS8nZdcUXm9Q48UP/3zZuruz5l2HJIPTPM6wkcV3BrLO2IItdYWF2U6/4+80z9fscdwdzfUYQpXdxQ6J13pgTGXh54QB86Bx2U+eYvokYdaFjDTzkbY9WuF/OSKVwzcqS+kecS6od5SLoPGzd9RCkZY/mEKYOEKCGYMVaoXswliG5s8uTctVhdvaAbonbvX1u2RK/jKqdC2K53KlvB8HyNMb9u7KGH9JzZd1/1jGXj8st1/JvfhM+VVymEdbJ3OEwi0hk4EYgoV7zFJYoeld5wXBhcr5ybxysXURpjBx2keWe2bEl5wLy4b1uXXZb9jd41xu67TxMZeklKvO9iPWPJc++96YspewnzkBwypL0RUCrifSjMM5YrsWeuMGVrayqNTKHGWBDd2JVX5tbRtrZqLz43/YWbyiEOHVc5FcLu10//740bMxvuURljbmTjK1/JHc4/80zYfXf1aCdRo7QUydWb8moRaRWRVtQQ+4/73TN9G/BTII3/wVIIUYQp8zU6dt9dx0ENwSiNMUiFmG66qf1b5+uva0ikZ8+UsZWJ/fdXg2XjRnj00dT0HTv0ohdJTgC/xx7a9f2jjwq7KW/ZoiGG2tp4SzhVAlE/JDt1ShVFhtL0jOUbpsxGLs/Y22/rNbbXXmqwFsLBB6v2aOHCzDmrGhtTxlg6vWDnzvpfNTaGT3+RD+VUCFskt24sCmPslVfUW9qzp/aYzEVdncoLoHqF/Lk8Y4+hGfevAAS4wfnuHf4bmGiM+XaM7axKCvWMecsg5esZK5YxduqpqnNbulQ9Wy7uhXrRRcH2lS5UOXeu6hlGjoy+WHAmamo0xw8U9gB44w39Xw84QDt5WDITx0PSDVVCaXnG8glT5sox5pLLGHNDlEcfHXzfmcilG1u5Er71Lf38y1+m1wu+8077vIOZiErHVW6FsLMZYzt3pnoK9+0bbrveHpVuOosvfSl4qbkpU/Ql8+GH4y8YX4rk6k35sjHmN8aYm4EvAj9zvnuH240xzyTT3OqiUM/Y4sWqbxk0KPyDI6wxFlVvSheRlHbs+utVR7B+fSpVhZsCIxef+5zeWB59NNUrNWm9mEsUoUqrFwtO1A/Jpqb2Yfuf/jS/FBlx0Levnufr1ukDNQhBPWOut2vZsvRh3yjE+16yhSqvukp/4ymnaKqETOkd/v3v3PuJSsdVboWws4n4Xc9q//7hEiM3NcFf/qKfH3lEe0YCfOpTwbcxeLC+hLe26n8YtMRVpRBGM/Y6cGS6GSLycREJ6Xux5KJQAX++ejHI3zMWpafp/PM1j9Orr6qRt9tuqv3YY4/cAl6XwYPhpJP0xnv//TqtEowxqxfLTZQPSTeP1HvvpaZt3x5fHqmw1NSkwqZudYBcBNWMdemiL3OtrenvRVGJ910yFQ1/6in1cHftqmXQsulFk9RxFZoWKGm8In4/+YQo3WvDX9JIBD75yeDXxvTp+h+DvlAEKXFVSYQxxm4kgzEGHO7Mt0RIoWHKQkTqxQ5Tgr4Zuz3Xtm5NTV++PNzF6eYcc0OV5WyMWfF+cKJ6SHrzSPmF46WURyqsbiyoZwwyhyqXL1d9V48e0WkYx4zRNi1erCEvUMP30kv18/e/n/s/S1rHVUhW/KTJFqZ0jbF0uRvT4b02vNUHQI2poNeGux1/RysorWssTsIYY4cAz2WY9wIwrvDmWLwUGqYsxDPWt696n9aty15U1iVqY8y9ONN1GW9tDXdxnnmmZvd/5hm9uRdyXAqhUGNs2zbVw9TUJN/2ciWKh2S55JEK26MyqGYMMhtjrlfsqKPyq/eZjk6d4LDD9PP+++t/1quXVj4YMQK++c3c2yiGjivfrPhJ44YAFy5MJUJ2CesZi+raKJdrLE7CGGOdgExBqG5AleYCj48+ffSGsX59MIPITyGesZqacMZg1MZYlBdnjx7wsY/p5333hTVr9PP//V+yb1qFGmNvvqk3+dGj1atjCUahD8lyySMVVsQfxjPmprfIZIxFFaIE9Xi7iZq3b1cPi6uDW7o0FcrKRrnpuJKkSxdNYNzW1vH+F9YYi+raKJdrLE7CGGMvA1MyzJsCzC68ORYvIil3cdiac94ySK5bOixhQpVRC/ijvDinT08Jer3bTFqLMHSo/h9LluRnXFvxfnEolzxSYT1jQTVjkPKM+Xu5ueL9KHpSQsojnqkTwvbtwTzi5abjSppMIv6wxlhU10a5XGNxEsYYuwb4mIi8KCJTReRsEblMRF4ETgB+EEsLq5x8RfxvvaVvlKNH51+/MIwxFrWAP6qL0725+93xkLwWobY2VVJnwYLw61u9WHEolzxSYTVjhYYpt27Vc7KmRsOUURClR7ycdFxJk0nEH9YYi+raKJdrLE4CG2PGmKeBU4A24NfA/cCvgBbgZJveIh7yFfFHoYvKxxiL6mKJ6uIsNS1CIaFK6xkrDuWSRyrpMOXLL6toe8yY4LmkchF1uKpcdFxJk0nEH9YYi+raKJdrLE5ClUMyxjQaY8YDPYA9gZ7GmAnWEIuPfEX8UZT7KaYxFtXFWWpaBNcYmz8/3Ho7dqi3UwTGjo2+XZbMlIv+KOkwZRx6MRuuSoaojLGoro1yucbiJGxtSgCMMVuNMUuNMVtzL20phGr1jEV1cZbazT1fz9icOeqFGDWqsl31pUi56I/CeMbcJMoQzDM2ZIi+CCxfnkphEIcxZsNVyeA1xryFucMaY1FdG+VyjcVJYGNMRP6Wa4izodVKPp4xY6L1jAUpFh61MRbVxVlqN/d8jTGrFysu5aA/CqMZ27ZNRfJdugQrq1VXpy+GbW16P2hri644uBcbrkqGgQP1/F23LpUkuKVFe5qLaEHxoER1bZTDNRYnYTxj/dMMo4BPAxOAnH+fiPxeRFaKyFueadeIyFIRed0ZPh7qF1Q4+Qj4Fy3SZKkDBxZWPy9MsfCoe1NCNBdnqd3c8zXGrF6s+JS6/shrjHm9HekI4xVz8erG5s3TB/mQIakQZhTYcFUypCsYvmqVnjf9+gWvcOIS1bVR6tdYnIQR8J+QZjgY2AdYTrAM/HcBp6WZfqMxZqwzPBa0TdVAPmHKfIuD+wkapjQmnnJIUPjFWWo397320qSWH3yQvoenl6YmrcvWsyfcdptOe/rpys5Cbcmf+noV0jc3p4ytTITRi7l4dWPeEGW2skRhseGq5PD3qMynFJIlOvLSjHkxxnwI/BS4LsCyTwNrC91nNZFPmDKqcj8DBqg3atWq7CL4bdvUIOvaNbos3FFRajf3zp3VIDNGM2Bnwq33Nm1aqiQUwGOPVUedNkt+BA1V5uMZ86a3iEMv5lLt4aqk8HvGrDFWXEI6IzPSCuxRwPqXi8gX0MSxVxlj1qVbSESm4CSeHThwII2NjQXsMjk2b96cd1t37qwBjmP58jZmzHiamgDm81NPHQD0p65uLo2NIbPF+ujdezxr13bh739/gf7907ty1q2rAybQpctOGhufL2h/hZLuWNfXw223deW++/bg3/8exLZtnaivb+Xkkz/inHOWUF+/nSRPpT59xrBgQR8eemgO48ev6TB/6dKuXHzx4Wzf3tGybW7W4eyzW5k27WWGDElTzC0hCjmvLcEJc5y7dh0H9OKJJ15j+fINGZd74YU+wBja2tbQ2Dgn0LZ37NgTGMELLyzhxRf7AA106TKbxsZ4er9MmqSDlw8/7FgFIEqq6ZxuaekPHMDzz6+msfEtGhsHAvsDK2hsnBvrvqvpOAfGGBNoAEanGcYC5wHvAzMCbmcY8Jbn+0C01FINcC3w+yDbOfTQQ025MHPmzILW3203Y8CYVauCLT9ypC7/xhsF7dYYY8y4cbqtl17KvExTky4zbFjh+yuUQo91EkydqsfrxhvTz7/0UmPq6nSZTENdnTGXXZZsu/2Uw7GuBMIc5zPP1PPjvvuyL3fvvbrceecFb8df/qLrTJig44YGY5qbg69fDlTTOf3mm/o/7rOPfr/uOv1+5ZXx77uajrMXYLbJYNOECVO+BczxDa8AfwTWABfnaQyuMMa0GmPagNuBI/LZTiUTRsTvLYPkagIKIUiPyjjE+5VMLhF/qeVGs5QPQXONFRKmdEOURx4ZXuhtKR1GjlS934IF2rPWhimLS5hL6YQ007YDS4wxS/NtgIgMNsa4EvGzUKPP4mHQIJg7V4X0Bx6YfVm3DNL+++dfBslLEBF/XOL9SiWXMVZqudEs5UNcmrGmJvjd79pP27hRp1sxfXlSXw/Dhql2dcECa4wVm8DGmDFmVqE7E5E/AxOBfiKyBLgamCgiYwEDLAK+XOh+Ko0wIv6oelK6BElvEXWOsUonlzHWvXt70X4m7PG2+Ama+DVMXcrp01W75ffWvv663mfuv9+K6suVUaPUGHv33dTzxRpjxSGrMSYiQ8NszBjzQY7556WZfEeYfVQjQdJbNDVpHcZp0/T7ffdpN/errirszTWMZ8waB8HYe28NDyxapA84fy/PyZP1f8wWqrSJLy3pCBumzJXaoqlJDbGtaWqttLbq9EmT9CXQesjKj1Gj4PHHNW+c9YwVl1yasUXAwhCDJQZyeca8aRDcB/iOHfq90DQI1hiLnq5dYehQfZgtXtxxfqnlRrOUD0E9Y0HDlNdfH0y/eGOQLJOWksOb3sIaY8UllzH2KTTD/qeBycAy4CngMuAcZzzDmX5+fM2sbrIJ+L1vrv6bZnNz6s0130Sh1hiLh2yhyhEjVMSfDpv40pKNoJqxoGFK25mksnGNsXfegdWr9XMhVVss+ZPVGDPGPOoOwKnAP40xpxhjbjHGPOiMTwYeBT6RRIOrkWxhyrjfXG1vynhwjbH589PPX7JExwMG2MSXluBEHaa0nUkqG7fH/auvar3Rvn1ze+Ut8RAmtcXZwIMZ5j2Aes8sMZAtTBn3m6trCK5YoWG1dNjelOHJ5hnbuRN+/nP9fMst1VmnzZIfvXtruokNG2B7lnzAQcOUQV+w7ItYeTJ4sP537jPEhiiLRxhjbBtwTIZ5x6JpLiwxkM0zFveba+fOWji2rU3LImXbtr0hByebMfanP2mW8f33hzPOSLZdlvKmpgb699fPma5XCG6MTZ4cTL9oO5OUJ96C4WCNsWISxhj7HfB9EblZRE4RkbHO+DfAd4Fb4mmipU8fveFt2KB1IL0E9UYVYijl0o1ZYyw8mYyx1lb42c/083e+Q6DyVxaLl1yhyra2lDHWq1f2bdnOJJWPNcZKg8C3emPMNcDX0cSsj6PZ9x93vn/DGHN1HA20aEI+N4Frt26qHZo6VZMw7tyZe/1C31ytMRY9w4freOFCaGlJTX/oIe3ZtNde8LnPFadtlvImV4/KzZvVIOvePbehNWKEdhZpaOi4rO1MUv40NbXXrT7wgD5b8u3wZcmfUO/dxphfAXsCw4Gjgb2BPYwxv4yhbRZSaStckbwxmhD01lv1otm5Ezp1rCfdjkLfXHMZY1bAH56GBhgyRLUabuFjY+CnP9XP/+//WSGtJT9yecbCZt8//XTtNDJliu1MUkm4z5ZXX01Na26OJiWSJTyhgyDGmDZjzCLgVWPMYqempCUGsiVcbHOOeufOapjF+eYa1DNmBfzh8Icq//UvvTEOHAhf/GLx2mUpb3Klt8inLuWIEdp5xHYmqQy8zxZ/x6woUiJZwhPKGBORo0VkuohsAraLyCYReUxExsfUvqomSNoKY+C11+J9c82V3sKGKfPDb4z95Cc6vvJKrRtnseRDrjBlmFJIlsrEJvMtPQIbYyJyMtAI7AH8HJjqjPcAGkXkpDgaWM2ESVsR55ur1YzFg9cYe+45ePppFVRfemlx22Upb4KGKXPlGLNULjaZb+kRuFA4cC3wCHCOMcZ4pv9YRB4AfgI8GWXjqp1SSbiYq1i4NcbC09QEM2fq5xtugJtu0s/nn69eTYslX3J5xvIJU1oqi1J5tlhShAlTHgTc7jPEXG5z5lsipFQSLlrPWLS4wtkZM1LT3B6Vd95phbOWwohDM2apLErl2WJJEcYYWw9kCnqNdOZbIqRUEi56KwCkM8Vtb8rgeIWz3pQWLtu2WeGspTByhSldzZgNU1YvpfJssaQIY4zdB/xURCaLSFcAEekqIpPREObf4mhgNVMqCRfr61XLtHMnrF3bfl5rqxoQIlZ0HgQrnLXEjZuBf+XKVK9rL9YzZimVZ4slRRhj7FvAP4G7gS0isgHY4nz/pzPfEiGllHAxU49K1yvWrZsaZJbsWOGsJW66dNGXp9bWlBfMizXGLKX0bLEoYTLwbzPGnA8cAFyEesMuAg4wxkw2xtjalDFQKgkXM+nGrF4sHFY4a0mCbKFKm9rCAqXzbLEogXpTOmHJDcBnjTF/B+bF2ipLO9y0FTffXLw2WGMsGrp31woKQZazWPJlwAB47z3tUbn//u3n2dQWFpdSeLZYlECeMcfrtRJIIzm2VAOZ0ltY8X44rHDWkgTZPGM2TGmxlB5hNGO3Al8VEVsxrwrJ5RmzpZCCYYWzliTIlmvMGmMWS+kRJulrb+BAYJGIPAWsALyJDowxxor4KxQbpowGVzg7aZIK9b1i/ro6Haxw1lIo2XKN2dQWFkvpEcYY+wyww/l8bJr5BtujsmKxxlh0uMLZG2/UXpObN+vxu+AC9YhZQ8xSKJnClC0tqlkUgR49km+XxWJJT2BjzBizd5wNsZQ2mVJbWGMsP6xw1hInmcKUGzfquFcv7T1nsVhKg5zGmIjUAx8HhgHLgaeMMRmqnlkqFa9nzJhUTjFrjFkspUemMKXVi1kspUlWY0xEhqPFv4d5Jm8UkXONMf+Ks2GW0qJHDxXpb9mib9e9eul025vSYik9MoUprV7MYilNcjmqrwPaUI1YA5rw9TW0Z6WlihBJrxuzvSktltIjU5jSesYsltIklzE2Hvi+MeY5Y8x2Y8xc4MvAUBEZHH/zLKVENmPMesYsltKhVy/o3Fmvz61bU9OtMWaxlCa5jLHBwALftCb+f3v3HyxXeddx/P1Nclt+VkBCKIQaBWRsgcKYVgQGU1pa0iJ0pkh/SAfHTjMOaimiDFZHrU47Vm1xxuAfFrEISPkhtMiYaWlqSiFpadA0xfIjhJ+BmABJIGlouIGvfzxn524uS7ib7O45u3m/Zu6c3WdP7n32O5fwyfd59hwI4NC+zEiNZRiThkPERHfsmWcmxg1jUjNN5fM0+fqnaE/Q6ROVhjGpmTotVbpnTGqmqVza4hsR0ek2SIsnj2fmIb2ZlpqoU2fMDfxSM3XaxG9nTGqm1wtjnx3ILDQU3MAvDY9OnTHDmNRMOw1jmdnTMBYRVwFnAesz89hq7CDgBsrlMx4DzsvMjb38ueqNTjcLd5lSaqZO1xpzmVJqpkFfg/krwJmTxi6jXEj2aGBx9VwN5AZ+aXi4TCkNj4GGscy8E9gwafgc4Orq8dXABwc5J02dYUwaHi5TSsOjmxuF98uszFwLkJlrI+I1PwQQEQuABQCzZs1iyZIlg5nhbtqyZcvQzHVnMmFs7DReeGEaixbdyd57v8KmTacAY6xYcRePPNLpcx6DNSq1HgbWejB2tc5r1x4IvJ0HH9zIkiU/BOCpp94B7MuqVT9g+/af9HSeo8Df6cGwzq8WmYO9ckVEzAFub9sztikzD2h7fWNmvu6Ohrlz5+by5cv7Ns9eWrJkCfPmzat7Gj0xZw48/jisWgVHHVUuLDk+Dtu2lcd1G6VaN521HoxdrfOKFXDiiXDccbByZRk77LDS2V6zBg4/vLfzHAX+Tg/GnlrniLg3M+d2em3Qe8Y6Wde6mn91XP8656tG7UuVL71UgtiMGc0IYpImuEwpDY8mhLHbgAuqxxcAX69xLnod7Z+odL+Y1FwzZ5bjs8/Cyy+X7vWLL5Z/PO2zT71zk7SjgYaxiLgeWAYcExFrIuITwF8DZ0TEKuCM6rkaqr0zZhiTmmtsDA46CF55BTZsgOefL+MHHlhulySpOQa6gT8zP/oaL717kPPQrjOMScPjkENKEFu3roQzcIlSaqImLFNqiLSHMW+FJDVb+7XG3C8mNZdhTF1pv1m4t0KSmq19E79hTGquJlxnTEPEZUppeLTfEmla9U9vb4UkNY9hTF0xjEnDo32Z8o1vLI/tjEnNYxhTV2bOhOnTy6bg554rY4YxqZnalyn33788NoxJzeOeMXVl2jQ49NDyePXqcjSMSc3UvkzZ2jPmMqXUPIYxda21VLlqVTkaxqRmal+m3LixPLYzJjWPYUxda4Wxhx4qRz9NKTWTn6aUhoNhTF1rhbFHHy1HO2NSM3mdMWk4GMbUtVYY2769HA1jUjPttx/stRds3Qpr1pQx94xJzWMYU9daYazFMCY1U8TEUuXDD5ejnTGpeQxj6tphh+343DAmNVdrqfKll8rRMCY1j2FMXZvcGXMDv9Rcrc5Yi2FMah7DmLrmMqU0PNrD2N57T1yJX1JzGMbUtVmzyl6UFsOY1FytZUqwKyY1lWFMXZsxo9wWqcUwJjVXe2fMMCY1k2FMu6R9qdIwJjVXexjzshZSMxnG1LXVq8uNwltmz4YLL5y4V6Wk5nCZUmo+w5i6smgRHH/8xAUkATZvhiuvLOOLFtU3N0mv5jKl1HyGMU3Z6tVw7rnlat6ZO742Pl7Gzz3XDpnUFKtXw+WXTzy/6Sa72FITGcY0ZV/8YgldOzM+vuNf/pLq0epiX3fdxNj4uF1sqYkMY5qya6+dWhi75prBzEdSZ+1d7Mn/zdrFlprHMKYp27Klt+dJ6g+72NJwMYxpyqZ6CQsvdSHVyy62NFwMY5qy88+HsbGdnzM2Bh//+GDmI6kzu9jScDGMacouuWRqYeziiwczH0md2cWWhothTFN25JFw882wzz6vDmVjY2X85pvLeZLqYxdbGi6GMXVl/nxYuRIWLIA3vQmmTSvHBQvK+Pz5dc9Qkl1sabgYxtS1I4+EhQvh+efh5ZfLceFCO2JSU9jFloaLYUySRpBdbGl4zKh7ApKk/mh1sRcurHsmknbGzpgkSVKNDGOSJEk1MoxJkiTVKDKz7jnskoh4Bni87nlM0cHAs3VPYg9hrQfHWg+GdR4caz0Ye2qdfy4zZ3Z6YWjD2DCJiOWZObfueewJrPXgWOvBsM6DY60Hwzq/msuUkiRJNTKMSZIk1cgwNhj/VPcE9iDWenCs9WBY58Gx1oNhnSdxz5gkSVKN7IxJkiTVyDAmSZJUI8NYn0XEmRHxYEQ8HBGX1T2fURIRV0XE+oi4r23soIi4IyJWVccD65zjKIiIIyLivyLi/oj434i4qBq31j0WEXtFxD0R8cOq1p+txq11H0TE9Ij4n4i4vXpunfsgIh6LiB9FxIqIWF6NWes2hrE+iojpwBXAfOCtwEcj4q31zmqkfAU4c9LYZcDizDwaWFw91+7ZDlySmb8EnAT8bvV7bK17bxtwema+HTgBODMiTsJa98tFwP1tz61z/7wrM09ou76YtW5jGOuvdwIPZ+YjmfkS8FXgnJrnNDIy805gw6Thc4Crq8dXAx8c6KRGUGauzcz/rh5vpvzP63Csdc9lsaV6OlZ9Jda65yJiNvAB4Mq2Yes8ONa6jWGsvw4Hnmx7vqYaU//Mysy1UEIEcEjN8xkpETEHOBH4Pta6L6qlsxXAeuCOzLTW/fH3wKXAK21j1rk/EvhmRNwbEQuqMWvdZkbdExhx0WHMa4loKEXEfsC/A5/OzBciOv16a3dl5svACRFxAHBrRBxb95xGTUScBazPzHsjYl7d89kDnJKZT0fEIcAdEfFA3RNqGjtj/bUGOKLt+Wzg6ZoMsXMcAAAFCElEQVTmsqdYFxFvBqiO62uez0iIiDFKELsuM2+phq11H2XmJmAJZV+kte6tU4CzI+IxyvaR0yPiWqxzX2Tm09VxPXArZQuPtW5jGOuvHwBHR8TPR8QbgI8At9U8p1F3G3BB9fgC4Os1zmUkRGmB/TNwf2Z+qe0la91jETGz6ogREXsD7wEewFr3VGb+cWbOzsw5lL+Xv52Z52Odey4i9o2I/VuPgfcC92Gtd+AV+PssIt5P2ZswHbgqMz9X85RGRkRcD8wDDgbWAX8OfA24EXgL8ATwG5k5eZO/uhARpwLfBX7ExP6az1D2jVnrHoqI4ymbmadT/rF8Y2b+ZUT8LNa6L6plyj/MzLOsc+9FxC9QumFQtkb9W2Z+zlrvyDAmSZJUI5cpJUmSamQYkyRJqpFhTJIkqUaGMUmSpBoZxiRJkmpkGJM0VCLiExGR1b0F28e/UI2fP2n8vdX4yT36+XOq73dWL76fJBnGJA2bpdVxcrg6GdjaYfxXgW3AvX2elyTtEsOYpGHzALCBttBV3a7pl4F/pXNIW56Z2wY2Q0nqgmFM0lDJcqXqZewYuk4EAvhH4Li2269MA36FqpsWEadGxHciYmtEPBcRX26d2xIRb4mIr0bEhuq8b0TEMTubU0TMi4jNEfH56vkBEXFlRDwdET+NiCci4su9qoGk0WIYkzSMlgInVPdvhLIUeS/lnnebKAEM4G3AzwB3R8QpwGLg/4BzgU8D7wf+pfVNI+Ig4C7gGOB3gPOAfYFvtf2sHUTE+4D/BP42Mz9TDX8JOBW4GHgf5fZR3u5EUkcz6p6AJO2CpcAY8A7gTkqXbFlmZkR8r3r+LSa6Z0uBW4Clmfnh1jeJiKeAxRFxbGbeRwlP+wIntO6TFxF3A48Bvw1c0T6JiDibcn+9P83Mv2t76Z3AFZl5Q9vYtb1445JGj50xScPoHmA7E2HrZMrSJcD3Jo2vAn5C6Z7dGBEzWl+ULtg4Zb8ZwHuAO4AX2s7ZTOm6zZ00hw8BNwGXTApiACuAP4qICyPiF3f73UoaaYYxSUMnM7dSAs/J1SUuZjMRxpYBJ1X7xU4G7gYOBKZT9pSNt31to3TYjqj+7MHAhyedMw68q+2clrMpHyS4tcMUfw/4GvBnwIMRsSoiPrJ771rSqHKZUtKwWgp8jBK4HsvMtdX494H9gV8DjgK+QNlHlsBfUPZ3TfZ0ddwA3Ab8VYdzNk96/vvAHwB3RMRpmflc64XM3AR8CvhURBwPXApcFxErM/PHXb5PSSMuygeTJGm4RMR5wA2UcPV8Zn6s7bWVwJOUDfpvy8wfR8RS4NHM/M2dfM/PUzbtH5eZL77GOXOAR4Ffp3ThvgP8FDg9M194jT/zZkrg+1Bm3tLlW5U04uyMSRpWd1fH+cBFk15bBnwS2AjcX41dStms/wpwM6XT9RbgA8CfZOZDlE9Bng98OyL+AXgKmEXpst2Vmde3/5DMfC4izgC+C9weEWdm5taIuIuyfHkfpSP3Scq+tXt69eYljQ73jEkaSpn5FPAE5fpiyya9vKw1Xl2XjMy8CzgNmAlcA/wHJaA9CayrznkWOIlyYdnLgW8Cf0O5PMbK15jHWuDdwBzgloh4Q/Xzf4sS+m6k7EWbn5lrdvuNSxo5LlNKkiTVyM6YJElSjQxjkiRJNTKMSZIk1cgwJkmSVCPDmCRJUo0MY5IkSTUyjEmSJNXIMCZJklSj/wdcqS8lAB4KggAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,4))\n",
    "plt.title(\"Typical factory production over a year\",fontsize=16)\n",
    "plt.plot(production,c='blue',lw=2,marker='o',markersize=10)\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"Weeks\",fontsize=15)\n",
    "plt.ylabel(\"Production (tons)\",fontsize=15)\n",
    "plt.hlines(y=production.mean(),xmin=-2,xmax=54,color='red',linestyle='--',lw=3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "hMisk5v3umZM"
   },
   "outputs": [],
   "source": [
    "n = len(production)\n",
    "m = production.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample standard deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1506,
     "status": "ok",
     "timestamp": 1567211402671,
     "user": {
      "displayName": "Tirthajyoti Sarkar",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AAuE7mD6d7dlMqdpzL4sermIF1ujmpSRxY2WnE4tuB-UsQ=s64",
      "userId": "01914075970409030121"
     },
     "user_tz": 420
    },
    "id": "uXeImgLCumZO",
    "outputId": "e4c0ae30-25ca-4497-92c1-f4d42e6a4f1f"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.919252299560373"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "production.std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample standard error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1496,
     "status": "ok",
     "timestamp": 1567211402671,
     "user": {
      "displayName": "Tirthajyoti Sarkar",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AAuE7mD6d7dlMqdpzL4sermIF1ujmpSRxY2WnE4tuB-UsQ=s64",
      "userId": "01914075970409030121"
     },
     "user_tz": 420
    },
    "id": "rNII4IARumZQ",
    "outputId": "2daf79cc-ee3e-476d-af24-2f0b93011c2e"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6821775539618872\n"
     ]
    }
   ],
   "source": [
    "std_err=production.std()/np.sqrt(n)\n",
    "print(std_err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 90% confidence interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1491,
     "status": "ok",
     "timestamp": 1567211402672,
     "user": {
      "displayName": "Tirthajyoti Sarkar",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AAuE7mD6d7dlMqdpzL4sermIF1ujmpSRxY2WnE4tuB-UsQ=s64",
      "userId": "01914075970409030121"
     },
     "user_tz": 420
    },
    "id": "tjCO4IQBumZU",
    "outputId": "5a599449-94f2-4646-e55e-6be4bae00b0e"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90% confidence interval of mean from  18.19130937896271  to  20.47618007932221\n"
     ]
    }
   ],
   "source": [
    "confidence = 0.9\n",
    "h = std_err * stats.t.ppf((1 + confidence) / 2, n)\n",
    "i90 =[m-h,m+h]\n",
    "print(\"90% confidence interval of mean from \",m-h,\" to \",m+h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 99% confidence interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1484,
     "status": "ok",
     "timestamp": 1567211402672,
     "user": {
      "displayName": "Tirthajyoti Sarkar",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AAuE7mD6d7dlMqdpzL4sermIF1ujmpSRxY2WnE4tuB-UsQ=s64",
      "userId": "01914075970409030121"
     },
     "user_tz": 420
    },
    "id": "Ww8G6CVmumZX",
    "outputId": "4b4a7bdf-33c7-4a7a-8a78-00f811ae6f72"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "99% confidence interval of mean from  17.509783661041904  to  21.157705797243015\n"
     ]
    }
   ],
   "source": [
    "confidence = 0.99\n",
    "h = std_err * stats.t.ppf((1 + confidence) / 2, n)\n",
    "i99 =[m-h,m+h]\n",
    "print(\"99% confidence interval of mean from \",m-h,\" to \",m+h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1477,
     "status": "ok",
     "timestamp": 1567211402672,
     "user": {
      "displayName": "Tirthajyoti Sarkar",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AAuE7mD6d7dlMqdpzL4sermIF1ujmpSRxY2WnE4tuB-UsQ=s64",
      "userId": "01914075970409030121"
     },
     "user_tz": 420
    },
    "id": "TUkK7WPfumZZ",
    "outputId": "75861ce1-7970-419a-e3cc-d404daddd87c"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(18.19130937896271, 20.47618007932221)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i90[0],i90[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1471,
     "status": "ok",
     "timestamp": 1567211402673,
     "user": {
      "displayName": "Tirthajyoti Sarkar",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AAuE7mD6d7dlMqdpzL4sermIF1ujmpSRxY2WnE4tuB-UsQ=s64",
      "userId": "01914075970409030121"
     },
     "user_tz": 420
    },
    "id": "mWojRg6PumZb",
    "outputId": "2721d35c-0322-4180-c406-0aa4f712eec2"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(17.509783661041904, 21.157705797243015)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i99[0],i99[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "yyikuQkG6LnN"
   },
   "source": [
    "### Repeat the random process many times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "pQyXDKYhumZd"
   },
   "outputs": [],
   "source": [
    "def repeat(n):\n",
    "    \"\"\"\n",
    "    Simulates the factory run `n` number of times\n",
    "    Counts the frequency where population mean (i.e. 20) is contained in the C.I.\n",
    "    \"\"\"\n",
    "    interval_90_count = 0\n",
    "    interval_99_count = 0\n",
    "    num_weeks = 52\n",
    "    \n",
    "    for i in range(n):    \n",
    "        production = np.random.normal(loc=20,scale=5,size=num_weeks)\n",
    "        m = production.mean()\n",
    "        std_err=production.std()/np.sqrt(num_weeks)\n",
    "        # For 90% C.I. \n",
    "        confidence = 0.9\n",
    "        h = std_err * stats.t.ppf((1 + confidence) / 2, num_weeks)\n",
    "        if m-h <= 20 <= m+h:\n",
    "            interval_90_count+=1\n",
    "        # For 99% C.I.\n",
    "        confidence = 0.99\n",
    "        h = std_err * stats.t.ppf((1 + confidence) / 2, num_weeks)\n",
    "        if m-h <= 20 <= m+h:\n",
    "            interval_99_count+=1\n",
    "    return (interval_90_count,interval_99_count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "HTb5JIMwumZf"
   },
   "outputs": [],
   "source": [
    "repeatations = 10000\n",
    "int_90,int_99 = repeat(repeatations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1458,
     "status": "ok",
     "timestamp": 1567211402674,
     "user": {
      "displayName": "Tirthajyoti Sarkar",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AAuE7mD6d7dlMqdpzL4sermIF1ujmpSRxY2WnE4tuB-UsQ=s64",
      "userId": "01914075970409030121"
     },
     "user_tz": 420
    },
    "id": "DtUlX-6fumZh",
    "outputId": "0bb67e2a-9903-493a-a2a4-3677fe94b8c3"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.898"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "round(int_90/repeatations,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1452,
     "status": "ok",
     "timestamp": 1567211402674,
     "user": {
      "displayName": "Tirthajyoti Sarkar",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AAuE7mD6d7dlMqdpzL4sermIF1ujmpSRxY2WnE4tuB-UsQ=s64",
      "userId": "01914075970409030121"
     },
     "user_tz": 420
    },
    "id": "VrJpX6L8umZk",
    "outputId": "a6ddbf00-a448-46bd-b2b4-630e90193254"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.99"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "round(int_99/repeatations,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary and thoughts for simulation\n",
    "In this notebook, we demonstrated the power of simulation to understand concepts of statistical estimation like expected value and confidence interval. In reality, we do not get the chance to repeat a statistical experiment thousands of times, but we can simulate the process on a computer, which helps us to distill down these concepts in a clear and intuitive manner.\n",
    "\n",
    "Once you master the art of simulating a stochastic event, you can investigate the properties of the random variables and the esoteric statistical theory behind them, with a new weapon of analysis.\n",
    "\n",
    "For example, you can investigate, using stochastic simulation,\n",
    "\n",
    "- The convergence of the mean of many stochastic events to a Normal distribution (verifying the Central Limit Theorem by numerical experiment)\n",
    "- Check what happens when you mix or transform many statistical distributions together in this way or that? what kind of resulting distributions do you get?\n",
    "- If a stochastic event does not follow the theoretical assumptions, what kind of aberrant behavior you can get in the result? In this case, the simulation could be your only friend because the standard theory fails if the assumptions are not met.\n",
    "- What kind of statistical properties emerges from the operation of a Deep Learning network?\n",
    "\n",
    "For learning the foundational principles of data science and machine learning, the importance of these kinds of exercise cannot be emphasized enough.\n",
    "\n",
    "![simulation-ds](https://raw.githubusercontent.com/tirthajyoti/Stats-Maths-with-Python/master/images/Simulation-problems-DS.png)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Conf_inv_sampling_hypothesis.ipynb",
   "provenance": [],
   "version": "0.3.2"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}