inferential_statistics_1a-Q6.25.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inferential Statistics Ia - Frequentism"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning objectives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Welcome to the first Frequentist inference mini-project! Over the course of working on this mini-project and the next frequentist mini-project, you'll learn the fundamental concepts associated with frequentist inference. The following list includes the topics you will become familiar with as you work through these two mini-projects:\n",
    "* the _z_-statistic\n",
    "* the _t_-statistic\n",
    "* the difference and relationship between the two\n",
    "* the Central Limit Theorem, including its assumptions and consequences\n",
    "* how to estimate the population mean and standard deviation from a sample\n",
    "* the concept of a sampling distribution of a test statistic, particularly for the mean\n",
    "* how to combine these concepts to calculate a confidence interval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prerequisites"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For working through this notebook, you are expected to have a very basic understanding of:\n",
    "* what a random variable is\n",
    "* what a probability density function (pdf) is\n",
    "* what the cumulative density function is\n",
    "* a high-level sense of what the Normal distribution\n",
    "\n",
    "If these concepts are new to you, please take a few moments to Google these topics in order to get a sense of what they are and how you might use them.\n",
    "\n",
    "While it's great if you have previous knowledge about sampling distributions, this assignment will introduce the concept and set you up to practice working using sampling distributions. This notebook was designed to bridge the gap between having a basic understanding of probability and random variables and being able to apply these concepts in Python. The second frequentist inference mini-project focuses on a real-world application of this type of inference to give you further practice using these concepts. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this notebook, we will use data sampled from a known normal distribution. This allows us to compare our results with theoretical expectations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I An introduction to sampling from the Normal distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's explore the ways we can generate the Normal distribution. While there's a fair amount of interest in [sklearn](https://scikit-learn.org/stable/) within the machine learning community, you're likely to have heard of [scipy](https://docs.scipy.org/doc/scipy-0.15.1/reference/index.html) if you're coming from the sciences. For this assignment, you'll use [scipy.stats](https://docs.scipy.org/doc/scipy-0.15.1/reference/tutorial/stats.html) to complete your work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "from scipy.stats import t\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from numpy.random import seed\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ Call up the documentation for the `norm` function imported above. What is the second listed method?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# help ( norm )\n",
    "# 2nd method in norm -- Methods inherited from scipy.stats._distn_infrastructure.rv_continuous:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ Use the method that generates random variates to draw five samples from the standard normal distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.11348847, 0.97448309, 0.72873463, 0.35146781, 0.70760514])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seed(47)\n",
    "# draw five samples here\n",
    "samples = np.random.random ( 5 )\n",
    "samples "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ What is the mean of this sample? Is it exactly equal to the value you expected? Hint: the sample was drawn from the standard normal distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.575155828796599"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate and print the mean here, hint: use np.mean()\n",
    "mean = np.mean ( samples )\n",
    "mean "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ What is the standard deviation of these numbers? Calculate this manually here as $\\sqrt{\\frac{\\sum_i(x_i - \\bar{x})^2}{n}}$. Hint: np.sqrt() and np.sum() will be useful here and remember that numpy supports [broadcasting](https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3045536451712141"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "std = np.sqrt ( np.sum ( ( samples - mean )  ** 2 ) / len ( samples ) ) \n",
    "std "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we have calculated the actual standard deviation of a small (size 5) data set. But in this case, this small data set is actually a sample from our larger (infinite) population. In this case, the population is infinite because we could keep drawing our normal random variates until our computers die. In general, the sample mean we calculate will not be equal to the population mean (as we saw above). A consequence of this is that the sum of squares of the deviations from the _population_ mean will be bigger than the sum of squares of the deviations from the _sample_ mean. In other words, the sum of squares of the deviations from the _sample_ mean is too small to give an unbiased estimate of the _population_ variance. An example of this effect is given [here](https://en.wikipedia.org/wiki/Bessel%27s_correction#Source_of_bias). Scaling our estimate of the variance by the factor $n/(n-1)$ gives an unbiased estimator of the population variance. This factor is known as [Bessel's correction](https://en.wikipedia.org/wiki/Bessel%27s_correction). The consequence of this is that the $n$ in the denominator is replaced by $n-1$.\n",
    "\n",
    "__Q:__ If all we had to go on was our five samples, what would be our best estimate of the population standard deviation? Use Bessel's correction ($n-1$ in the denominator), thus $\\sqrt{\\frac{\\sum_i(x_i - \\bar{x})^2}{n-1}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.34050132669909267"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "std_u = np.sqrt ( np.sum ( ( samples - mean )  ** 2 ) / (len ( samples ) - 1) ) \n",
    "std_u "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ Now use numpy's std function to calculate the standard deviation of our random samples. Which of the above standard deviations did it return?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3045536451712141"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.std ( samples )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ Consult the documentation for np.std() to see how to apply the correction for estimating the population parameter and verify this produces the expected result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.34050132669909267"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.std ( samples , ddof = 1 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary of section"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, you've been introduced to the scipy.stats package and used it to draw a small sample from the standard normal distribution. You've calculated the average (the mean) of this sample and seen that this is not exactly equal to the expected population parameter (which we know because we're generating the random variates from a specific, known distribution). You've been introduced to two ways of calculating the standard deviation; one uses $n$ in the denominator and the other uses $n-1$ (Bessel's correction). You've also seen which of these calculations np.std() performs by default and how to get it to generate the other."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You use $n$ as the denominator if you want to calculate the standard deviation of a sequence of numbers. You use $n-1$ if you are using this sequence of numbers to estimate the population parameter. This brings us to some terminology that can be a little confusing.\n",
    "\n",
    "The population parameter is traditionally written as $\\sigma$ and the sample statistic as $s$. Rather unhelpfully, $s$ is also called the sample standard deviation (using $n-1$) whereas the standard deviation of the sample uses $n$. That's right, we have the sample standard deviation and the standard deviation of the sample and they're not the same thing!\n",
    "\n",
    "The sample standard deviation\n",
    "\\begin{equation}\n",
    "s = \\sqrt{\\frac{\\sum_i(x_i - \\bar{x})^2}{n-1}} \\approx \\sigma,\n",
    "\\end{equation}\n",
    "is our best (unbiased) estimate of the population parameter ($\\sigma$).\n",
    "\n",
    "If your data set _is_ your entire population, you simply want to calculate the population parameter, $\\sigma$, via\n",
    "\\begin{equation}\n",
    "\\sigma = \\sqrt{\\frac{\\sum_i(x_i - \\bar{x})^2}{n}}\n",
    "\\end{equation}\n",
    "as you have complete, full knowledge of your population. In other words, your sample _is_ your population. It's worth noting at this point if your sample is your population then you know absolutely everything about your population, there are no probabilities really to calculate and no inference to be done.\n",
    "\n",
    "If, however, you have sampled _from_ your population, you only have partial knowledge of the state of your population and the standard deviation of your sample is not an unbiased estimate of the standard deviation of the population, in which case you seek to estimate that population parameter via the sample standard deviation, which uses the $n-1$ denominator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You're now firmly in frequentist theory territory. Great work so far! Now let's dive deeper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II Sampling distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we've been dealing with the concept of taking a sample from a population to infer the population parameters. One statistic we calculated for a sample was the mean. As our samples will be expected to vary from one draw to another, so will our sample statistics. If we were to perform repeat draws of size $n$ and calculate the mean of each, we would expect to obtain a distribution of values. This is the sampling distribution of the mean. The Central Limit Theorem (CLT) tells us that such a distribution will approach a normal distribution as $n$ increases. For the sampling distribution of the mean, the standard deviation of this distribution is given by\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{mean} = \\frac{\\sigma}{\\sqrt n}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\sigma_{mean}$ is the standard deviation of the sampling distribution of the mean and $\\sigma$ is the standard deviation of the population (the population parameter)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is important because typically we are dealing with samples from populations and all we know about the population is what we see in the sample. From this sample, we want to make inferences about the population. We may do this, for example, by looking at the histogram of the values and by calculating the mean and standard deviation (as estimates of the population parameters), and so we are intrinsically interested in how these quantities vary across samples. In other words, now that we've taken one sample of size $n$ and made some claims about the general population, what if we were to take another sample of size $n$? Would we get the same result? Would we make the same claims about the general population? This brings us to a fundamental question: _when we make some inference about a population based on our sample, how confident can we be that we've got it 'right'?_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's give our normal distribution a little flavor. Also, for didactic purposes, the standard normal distribution, with its variance equal to its standard deviation of one, would not be a great illustration of a key point. Let us imagine we live in a town of 50000 people and we know the height of everyone in this town. We will have 50000 numbers that tell us everything about our population. We'll simulate these numbers now and put ourselves in one particular town, called 'town 47', where the population mean height is 172 cm and population standard deviation is 5 cm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed(47)\n",
    "pop_heights = norm.rvs(172, 5, size=50000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd7xcVbn/8c9XQMqPTgDpEUEUrhfESFGBCFw6YkHFBiga5YKKwpUiF7g0Y0EUFTQKBgEJEUFCsQAaEGkSigiIxhglEulNeuD5/bHWSebU2ZPMnL3mnO/79ZrXzOz6zDNrZs1ae8/aigjMzMyG8oq6AzAzs/K5sjAzs6ZcWZiZWVOuLMzMrClXFmZm1pQrCzMza8qVRYdI+q6k/23TttaV9G9Ji+Xn0yV9vB3bztv7uaT92rW9FvZ7oqSHJf1rgHnjJc1ZyO1uI+neissu9H4qbLttZWBR5fKzft1xjASSZkvacSHXrVw2S+PKYiHkwvKspKckPS7pekmfkjQ/nxHxqYg4oeK2hix4EfGPiFg2Il5qQ+zHSTq3z/Z3jYizF3XbLcaxDnAosHFEvKqd246I30bERu3YlqTJkk5cyDgqlYF2G+jHRC4/s9q0/f0lXdeObY10kkLSBj3P21k2h5sri4W3Z0QsB6wHTAQOB85s904kLd7ubRZiPeCRiHiw7kBGu54Wq9mQIsK3Fm/AbGDHPtO2AF4G/iM/nwycmB+PAS4DHgceBX5LqqjPyes8C/wb+AIwFgjgAOAfwLUN0xbP25sOfAm4GXgCuARYOc8bD8wZKF5gF+AF4MW8vzsatvfx/PgVwNHA34EHgR8BK+R5PXHsl2N7GPjiEHlaIa//UN7e0Xn7O+bX/HKOY/IA644H5pBaHw8Cc4GPNsxfEvhajuMB4LvA0gPlANgcuA14CvgJcEHDezPofoAJOVcv5DgvzdMPB/6Zt3cvsMMgr39ylf0Mkbsz83L/BE4EFsvz9geuy6//MeBvwK553knAS8BzOeZv5+kBbNAQ1xnAFcDT+f0YNJ994np93vZLefuPD/Ve53l/B96UH384x7Jxfv5x4Gf58XHA1Lydp4C7gHFD5CiAzwCzSGXxqw37rFKOJwD35xwfOtD7Nkh5mk3+/JM+9zeQPttzgW8Dr8zzrs37eTrn6v0DbOv1pM/f4/n1vqNPHN8BLs/5uAl4TV3fe25ZtElE3Ez6MthmgNmH5nmrAqsDR6VV4iOkD+eekboJvtKwznakgrTzILvcF/gYsCYwDzitQoy/AE4GLsj723SAxfbPt7cD6wPLkj4Ajd4GbATsABwj6fWD7PJbpC+R9fPr2Zf0BXkVsCtwf45j/0HWf1Vefy1S5fkdSSvleV8GXgtsBmyQlzmm7wYkvRK4mPTBWxk4H3hXlf1ExCTgPOArOc49JW0EHAy8OVLLcmfSl0cVQ72evs4mva8bAG8EdiJ9sfbYklRRjQG+ApwpSRHxRdKPkYNzzAcPsv0PkiqW5UgVT6V8RsQ9wKeAG/L2V8yzBnyv87xrSF+SANuSvty3a3h+TcMu3gFMAVYEptG/7PX1LmAc6QfBXqTPBFQrx28HNiTl9oiFPA7xEvA50vuwNekz8d8AEbFtXmbTnKsLGleUtARwKfArYDXg08B5uYz1+ADwf8BKwEzSe1YLVxbtdT/pC6mvF4E1gPUi4sVI/ZbNBuU6LiKejohnB5l/TkT8MSKeBv4XeF+buhM+BHw9ImZFxL+BI4F9+nSH/V9EPBsRdwB3AP0qnRzL+4EjI+KpiJgNnAJ8pIVYXgSOzzm7gvTrbCNJAj4BfC4iHo2Ip0iV4D4DbGMrYHHgtLydi0gtsqb7GSSml0i/wjeWtEREzI6Ivy7K6+m7kKTVSZXpIbkMPAic2uf1/T0ivh/pONbZpPK1esU4AC6JiN9FxMvA81TPZz8V3utrWFA5bENqFfc8347elcV1EXFFfl3nMEDZ6uPLOeZ/AN8gfblC9XL8dETcCfywYd3KImJGRNwYEfPy6/5ew2trZitSJTYxIl6IiF+TeiAa47goIm6OiHmkHy6btRpju7iyaK+1SN1MfX2V9KvgV5JmSTqiwrbua2H+34ElSL9uFtWaeXuN216c3l9EjWcvPUMq8H2NAV45wLbWaiGWR/KHpO++VgWWAWbkEwweB36Rp/e1JvDPPpVz39wOtp9+ImImcAipy+RBSVMkrbmIr6ev9Ujv59yG1/c90q/PHvPfg4h4Jj8cMOZBNOaglXwOpNl7fQ2wjaRXAYuRugHfKmksqTVye8N6fcvWUk2O2/X9HPS8F1XK8WDrVibptZIuk/QvSU+SKtmqn8M1gftyhd0YR+NnpMpnbVi4smgTSW8mvcn9zhLJv7YOjYj1gT2Bz0vaoWf2IJts1vJYp+HxuqRfrQ+T+keXaYhrMXp/6Jtt937Sl1XjtueR+rFb8XCOqe+2/tnidgbb9rPAJhGxYr6tEBEDfZDmAmvl1kiPdQZYbjD98hURP46It5FeW5C6cNrpPtKv/TENr2/5iNik4vpVhpJuXKaVfA60/SHf61zBPkM6vnBtbrn8i3TM4Lo+X5at6vs5uD8/rlKOB1u312eI1H04mDOAPwEbRsTypC5mDbF8o/uBdRrPoqR9n5G2c2WxiCQtL2kPUj/rublJ23eZPSRtkL+wniR1ZfScBvsAqU+1VR+WtLGkZYDjgQtz0/3PpF9ju+c+0aNJ3SY9HgDG9imgjc4HPifp1ZKWZcExjnmDLD+gHMtU4CRJy0laD/g8cO7Qa1ba9svA94FTJa0GIGktSQMd37mBlOuDJS0uaS/SQcmqer0/kjaStL2kJUkHep9lwXvZFhExl9SPfUouX6+Q9BpJVbs3WipTLeazZ/tr5+NBVd/ra0jHenq6nKb3eb6w/kfSSvlU7M+SWi1QrRz/r6RlJG1COr7Ss+7twG6SVs6toUOG2P9ypM/0vyW9Djiwz/yh3oubSBXTFyQtIWk86cfklOYve/i5slh4l0p6ivQr8IvA11lwQK+vDYGrSH3UNwCnR8T0PO9LwNG5+X9YC/s/h3TQ9l/AUqRfbUTEE6QDbD8g/UJ5mnRwvcdP8v0jkm4dYLtn5W1fSzrL5jnSgbeF8em8/1mkFteP8/bb4XBS196Nufl/FQP0/0fEC8C7SQeUHyedjXMZ6Zd7FWeSjk88LulnpIp3IunX9L9IXUNHLdpLGdC+pK6du0lnPF1IOi5RxTeBvSU9JqnpiQ9ZpXxmvyadufMvSQ/nac3e62tIX6zXDvJ8YV0CzCB9wV/OgtPXq5Tja0iv+WrgaxHxqzz9HNKxuNmkSvsCBncY6WSBp0gVbt9ljwPOzuXnfY0zctl8B+n41MPA6cC+EfGnJq+5Fmp+nNVsZJF0E/DdiPhh3bHYwpMUpO6fmS2uN5ZUgSzRaot5NHPLwkY8SdtJelXuhtoP+E/SAVwzq2ik/jvYrNFGpD71ZYG/Anvn4wJmVpG7oczMrCl3Q5mZWVMjshtqzJgxMXbs2LrDsLrdm0eC3qiQQT4djxVuxowZD0fEgH/GHJGVxdixY7nlllvqDmN0uv76dP+Wt9QbB8D48el++vQ6o1jA8QytpLIzSkn6+2DzRmRlYTU6Kv/loJQvIOseLjtFc2VhI9fRR9cdQdmcH2uBKwsbuXZcqCtfjh7Oj7XAZ0PZyHX77elmA3N+rAVuWdjIdUge/8194ANzfqwFblmYmVlTrizMzKwpVxZmZtaUKwszM2vKB7itvU4+ue4IFigplgZjj7i80nKzJ+7e2UBKy09p8VgvriysvUoaqqGkWEpUWn5Ki8d6cWVh7VXS+D7DHEuzFsOUWY+kB1sNQzBVlPReQXnxWC+uLKy9Shrfp6RYSlRafkqLx3pxZWHt9b3v1R2BdSuXnaK5srD28rURbGG57BTNp85ae116abqZtcplp2huWVh7nXJKut9zz3rjsO7jslM0VxY2cn3jG3VHUDbnx1rgysJGrs02qzuCsjk/1gIfs7CR66qr0s0G5vxYC9yysJHrxBPTva8INzDnx1rgloWZmTXlysLMzJpyZWFmZk25sjAzs6Z8gNvaq6TxfUqKpUSl5ae0eKwXVxbWXiWN79OGWKpeqKgrlfReQXnxWC/uhrL2Kml8n5JiKVFp+SktHuvFLQtrr5LG9ykploXQ8cuvlpaf0uKxXjpaWUiaDTwFvATMi4hxklYGLgDGArOB90XEY5IEfBPYDXgG2D8ibs3b2Q84Om/2xIg4u5Nx2yK48MK6I7Bu5bJTtOHohnp7RGwWEePy8yOAqyNiQ+Dq/BxgV2DDfJsAnAGQK5djgS2BLYBjJa00DHHbwhgzJt3MWuWyU7Q6jlnsBfS0DM4G3tkw/UeR3AisKGkNYGfgyoh4NCIeA64EdhnuoK2iyZPTzaxVLjtF63RlEcCvJM2QNCFPWz0i5gLk+9Xy9LWA+xrWnZOnDTa9F0kTJN0i6ZaHHnqozS/DKvMH3haWy07ROn2A+60Rcb+k1YArJf1piGU1wLQYYnrvCRGTgEkA48aN6zffRqFzzqk7grI5P9aCjlYWEXF/vn9Q0sWkYw4PSFojIubmbqYH8+JzgHUaVl8buD9PH99n+vROxm0jxDrrNF9mNHN+rAUd64aS9P8kLdfzGNgJ+CMwDdgvL7YfcEl+PA3YV8lWwBO5m+qXwE6SVsoHtnfK08yGdsEF6WYDc36sBZ1sWawOXJzOiGVx4McR8QtJvwemSjoA+Afw3rz8FaTTZmeSTp39KEBEPCrpBOD3ebnjI+LRDsZtI8UZZ6T797+/3jhK5fxYCzpWWUTELGDTAaY/AuwwwPQADhpkW2cBZ7U7RjMzq8bDfZiZWVOuLMzMrClXFmZm1pQHErT2Kml8n5JiKVFp+SktHuvFlYW1V0lj+5QUS4lKy09p8Vgv7oay9ippyIaSYilRafkpLR7rxZWFtVdJH/iSYilRafkpLR7rxd1Q1l7Tp9cdgXUrl52iuWVhZmZNubKw9vra19LNrFUuO0VzZWHtddll6WbWKpedovmYhY1cV1xRdwRlc36sBa4sbORaZpm6Iyib82MtcDeUjVynn55uNjDnx1rgysJGrqlT080G5vxYC1xZmJlZU64szMysKR/gtlFp7BGX1x2CWVdxy8LMzJpyy8Laq6TxfUqKpUSl5ae0eKwXVxZmXa5ql9rsibt3OBIbydwNZe1V0vg+JcVSotLyU1o81otbFtZeN9xQdwQL9IwzdNhh9cZRqtLyU1LZsX5cWVh7/fSndUdg3cplp2juhjIzs6ZcWVh7HXlkupm1ymWnaO6GsvYqqd956aXrjqBspeWnpLJj/TStLCStDpwMrBkRu0raGNg6Is7seHRmi+LnP687grI5P9aCKt1Qk4FfAmvm538GDqm6A0mLSbpN0mX5+asl3STpL5IukPTKPH3J/Hxmnj+2YRtH5un3Stq56r7NzKw9qlQWYyJiKvAyQETMA15qYR+fBe5peP5l4NSI2BB4DDggTz8AeCwiNgBOzcuRWzL7AJsAuwCnS1qshf3baHXCCelmA3N+rAVVKounJa0CBICkrYAnqmxc0trA7sAP8nMB2wMX5kXOBt6ZH++Vn5Pn75CX3wuYEhHPR8TfgJnAFlX2b6Pc1Venmw3M+bEWVDnA/XlgGvAaSb8DVgX2rrj9bwBfAJbLz1cBHs+tE4A5wFr58VrAfZBaL5KeyMuvBdzYsM3GdeaTNAGYALDuuutWDM/MzKpoWllExK2StgM2AgTcGxEvNltP0h7AgxExQ9L4nskD7aLJvKHWaYxzEjAJYNy4cf3mm5nZwhu0spD07kFmvVYSEXFRk22/FXiHpN2ApYDlSS2NFSUtnlsXawP35+XnAOsAcyQtDqwAPNowvUfjOmZmNgyGalnsOcS8AIasLCLiSOBIgNyyOCwiPiTpJ6RurCnAfsAleZVp+fkNef6vIyIkTQN+LOnrpDOyNgRubvK6rC6rrFJ3BAuUFEuJSstPafFYL4NWFhHx0Q7t83BgiqQTgduAnv9rnAmcI2kmqUWxT47jLklTgbuBecBBEdHK2Vg2nEoa36ekWEpUWn5Ki8d6qfKnvFWAY4G3kVoU1wHHR8QjVXcSEdOB6fnxLAY4mykingPeO8j6JwEnVd2fmZm1V5VTZ6cADwHvIXUPPQRc0MmgrIuVNL5PSbGUqLT8lBaP9VLl1NmVI6LxnzsnSnrnoEvb6PZI5QZn53msoaGVlp+Syo71U6Wy+I2kfYCp+fneQLXrONroM2lS3RFYt3LZKVqVbqhPAj8GXsi3KcDnJT0l6clOBmdmZmWo8qe85ZotYzbfhAnp3r8SrVUuO0WrdD0LSe8Ats1Pp0fEZZ0Lybran/9cdwQLrL123RGUrbT8lFR2rJ8qp85OBN4MnJcnfVbS2yLiiI5GZraozj237gjK5vxYC6q0LHYDNouIlwEknU36M50rCzOzUaLqZVVXJP2rGtKYTWZFGnvEghP1jrkq9X0fv+OEusIp2yH5Gmbf+Ea9cVhXqFJZfAm4TdJvSCPAbkse88msZBs/OKvuEMp2++11R2BdpMrZUOdLmk46biHg8Ij4V6cDMzOzcjT9n0W+Wt0OpOMWlwCvlOQr1ZmZjSJV/pR3OrA18IH8/CngOx2LyMzMilPlmMWWEbG5pNsAIuIxSa/scFzWrV772rojmG/Wyv2uvmuNCnqvgPLisV6qVBYvSlqMfClTSasCL3c0KuteBf379qhdPl13CGUr6L0CyovHeqnSDXUacDGwuqSTSNezOLmjUZmZWVGqnA11nqQZpIPcAO+MiHs6G5Z1rYLG9zn5F98C3MLo0fgfFBg8P7Mn7j5sMfVSUNmx/qr+KW8ZoKcraunOhWNdr6DrKK//6D/rDqFoxeWnoLJj/VUZG+oY0uVOf0r6n8UPJf0kIk7sdHDWhb70pbojsG7lslO0Ki2LDwBvzNfI7hlY8FbAlYWZ2ShR5QD3bGCphudLAn/tSDTW/d7znnQza5XLTtGqtCyeB+6SdCXpmMV/AddJOg0gIj7Twfis2xR0HeW7V1u/7hCKVlx+Cio71l+VyuLifOsxvTOhmLWXR5sdmvNjrahy6uzZwxGImZmVq8oxC7OudOqlX+PUS79WdxjFcn6sFVX/Z2HWddZ46uG6Qyia82OtGLRlIemcfP/Z4QvHzMxKNFQ31JskrQd8TNJKklZuvDXbsKSlJN0s6Q5Jd0n6vzz91ZJukvQXSRf0jGAracn8fGaeP7ZhW0fm6fdK2nnRXrKZmbVqqMriu8AvgNcBM/rcbqmw7eeB7SNiU2AzYBdJWwFfBk6NiA2Bx4AD8vIHAI9FxAbAqXk5JG0M7ANsAuwCnJ5HwTUzs2Ey6DGLiDgNOE3SGRFxYKsbjogA/p2fLpFvAWwPfDBPPxs4DjgD2Cs/BrgQ+Ha+St9ewJSIeB74m6SZwBbADa3GZMNg663rjmC+W9d6Xd0hFK24/BRUdqy/KqfOHihpU2CbPOnaiPhDlY3nFsAMYAPS1fX+CjweEfPyInOAnivUrAXcl/c5T9ITwCp5+o0Nm21cp3FfE4AJAOuuu26V8KwTChrf5yvb7V93CEUrLj8FlR3rr8o1uD8DnAeslm/nSao05nNEvBQRmwFrk1oDrx9osZ5dDTJvsOl99zUpIsZFxLhVV121SnhmZlZRlVNnP066tOrTAJK+TOoC+lbVnUTE45KmA1sBK0paPLcu1gbuz4vNAdYB5khaHFgBeLRheo/Gdaw0PWP7/PSn9cYBnHFxukbXge86quZIylRcfgoqO9ZflT/lCXip4flLDPxrv/dK0qqSVsyPlwZ2BO4BfgPsnRfbD7gkP56Wn5Pn/zof95gG7JPPlno1sCFwc4W4rQ5bb11M3/NKzz7JSs8+WXcYxSouPwWVHeuvSsvih8BNknrGh3oncGaF9dYAzs7HLV4BTI2IyyTdDUyRdCJwW8O2zgTOyQewHyWdAUVE3CVpKnA3MA84KCJewsp02GF1R2DdymWnaFUOcH89dyG9jdSi+GhE3FZhvT8Abxxg+izS8Yu+058jXWRpoG2dBJzUbJ9mZtYZlYb7iIhbSRc8Mhva+PHpfvr0OqOwbuSyUzSPDWUj1u/W27TuEIrm/FgrXFnYiPWtt36g7hCK5vxYK4Y8G0rSYpKuGq5gzMysTENWFvmso2ckrTBM8Zi1zeSpxzJ56rF1h1Es58daUaUb6jngznwN7qd7Jvra21a6peY9X3cIRXN+rBVVKovL883MzEapStfgzv/AXjci7h2GmMzMrDBVBhLcE7iddG0LJG0maVqnAzMzs3JU6YY6jvSP6+kAEXF7HqPJrL899qg7gvmufk2/gQKsQXH5KajsWH9VKot5EfFEug7RfP2GCDcDihrf5/tbvrvuEIpWXH4KKjvWX5XK4o+SPggsJmlD4DPA9Z0Ny8zMSlKlsvg08EXSNbXPB34JnNDJoKyLdWh8n7FHtH5C3pQfHwHAPh+c2NZYRori8uOxoYpW5WyoZ4Av5oseRUQ81fmwrGvtv3/dEVi3ctkpWtPKQtKbgbOA5fLzJ4CPRcSMDsdm3cgfeFtYLjtFq9INdSbw3xHxWwBJbyNdEOk/OxmYdamHH073Y8bUG4cttKpdfrMn7t7eHbvsFK1KZfFUT0UBEBHXSXJXlA1s73zFXPc7W6tcdoo2aGUhafP88GZJ3yMd3A7g/eT/XJiV7LLXbVN3CEVzfqwVQ7UsTunzvHF4Sv/Pwop37uZt7iYZYZwfa8WglUVEvH04AzFrt6VefA6A55ZYquZIyuT8WCuqnA21IrAvMLZxeQ9RbqWb/JPjgIL+R1AY58daUeUA9xXAjcCdwMudDcfMzEpUpbJYKiI+3/FIzMysWE2HKAfOkfQJSWtIWrnn1vHIzMysGFVaFi8AXyWND9VzFlQA63cqKDMzK0uVyuLzwAYR8XCng7ERoKAhGy58w451h1C04vJTUNmx/qpUFncBz3Q6EBshCvrAF/dlWJji8lNQ2bH+qlQWLwG3S/oNaZhywKfO2iAKGt9npWeeAOCxZVaoOZIyFZefgsqO9VflAPfPgJNIFzya0XAbkqR1JP1G0j2S7pL02Tx9ZUlXSvpLvl8pT5ek0yTNlPSHhuFGkLRfXv4vkvZbmBdqw2TvvReM8VOzM372Jc742ZfqDqNYxeWnoLJj/VW5nsXZC7ntecChEXGrpOWAGZKuBPYHro6IiZKOAI4ADgd2BTbMty2BM4At85lXxwLjSAfWZ0iaFhGPLWRc1kmHHlp3BNatXHaKVuUf3H9jgLGgImLIs6EiYi4wNz9+StI9wFrAXsD4vNjZpEEJD8/TfxQRAdwoaUVJa+Rlr4yIR3M8VwK7kAY2tNLsuWfdEVi3ctkpWpVjFuMaHi8FvBdo6X8WksYCbwRuAlbPFQkRMVfSanmxtYD7Glabk6cNNr3vPiYAEwDWXXfdVsKzdrr33nS/0Ub1xmHdx2WnaFW6oR7pM+kbkq4DjqmyA0nLAj8FDomIJyUNuuhAux9iet84JwGTAMaNG+dRcevyyU+me1+TwFrlslO0Kt1Qmzc8fQWppbFclY1LWoJUUZwXERflyQ9IWiO3KtYAHszT5wDrNKy+NnB/nj6+z/TpVfZvo9u5b9yt7hCK5vxYK6p0QzVe12IeMBt4X7OVlJoQZwL3RMTXG2ZNA/YDJub7SxqmHyxpCukA9xO5QvklcHLPWVPATsCRFeK2Ue6y129bdwhFc36sFVW6oRb2uhZvBT4C3Cnp9jztKFIlMVXSAcA/SMdAII1uuxswk/QnwI/m/T8q6QTg93m543sOdpsNZY0nHwJg7vKr1hxJmZwfa0WVbqglgffQ/3oWxw+1XkRcx8DHGwB2GGD5AA4aZFtnAWc1i9Ws0amXpUaxr9cwMOfHWlGlG+oS4AnSH/Geb7KsmZmNQFUqi7UjYpeOR2JmZsWqMtzH9ZLe0PFIzMysWFVaFm8D9s//5H6edBwiIuI/OxqZmZkVo0plsWvHo7CRo6Dxfb6/xbvqDqFoxeWnoLJj/VU5dfbvwxGIjRAtju8z9ojLOxQIXL3Blh3b9kiwqPlp5b2bPXH35gt5bKiiVTlmYVbdvfcuGOOnZus/Mof1H5lTdxjFKi4/BZUd669KN5RZdQWN73PyL78N+H8EgykuPwWVHevPlYW118kn1x2BdSuXnaK5srD2estb6o7AupXLTtF8zMLa6/rr082sVS47RXPLwtrrqKPSvfudrVUuO0VzZWEj1rfesk/dIRTN+bFWuLKwEet3YzerO4SiOT/WCh+zsBFr4wdmsfEDs+oOo1jOj7XClYWNWMdcPYljrp5UdxjFcn6sFa4szMysKVcWZmbWlCsLMzNrypWFmZk15VNnrb0KGt/nK9vuV3cIRSsuPwWVHevPlYW1V0Hj+9y69uvrDqFoxeWnoLJj/bkbytqroPF9Np9zD5vPuafuMIpVXH4KKjvWn1sW1l4Fje/zhWvPBgq6XkNhistPQWXH+nNlYe31ve/VHYF1K5edormysPbaaKO6I7Bu5bJTNB+zsPa69NJ0M2uVy07R3LKw9jrllHS/5571xmHdx2WnaB2rLCSdBewBPBgR/5GnrQxcAIwFZgPvi4jHJAn4JrAb8Aywf0TcmtfZDzg6b/bEiDi7UzHbyHL8DhPqDqFow5mfsUdc3nSZKbMeYav1VxmGaGxhdLIbajKwS59pRwBXR8SGwNX5OcCuwIb5NgE4A+ZXLscCWwJbAMdKWqmDMdsIcvfq63P36uvXHUaxnB9rRcdaFhFxraSxfSbvBYzPj88GpgOH5+k/iogAbpS0oqQ18rJXRsSjAJKuJFVA53cqbmuPKr8kO+2ts28HfJGfwTg/1orhPmaxekTMBYiIuZJWy9PXAu5rWG5OnjbY9H4kTSC1Slh33XXbHLZ1o09fPwXwl+FgnB9rRSlnQ2mAaTHE9P4TIyZFxLiIGLfqqqu2NTgzs9FuuCuLB3L3Evn+wTx9DrBOw3JrA/cPMd3MzIbRcFcW04CeoS73Ay5pmL6vkq2AJ3J31S+BnSStlA9s75SnmZnZMOrkqbPnkw5Qj5E0h3RW00RgqqQDgH8A782LXz69racAAAlBSURBVEE6bXYm6dTZjwJExKOSTgB+n5c7vudgt5mZDZ9Ong31gUFm7TDAsgEcNMh2zgLOamNo1kk94/v8cGa9cQBH7Xxw3SEUrbT8HLXzwfz6sPF1h2GD8D+4rb3mj+9Tf2Uxa5W16w6haKXlZ9Yqa3t8qIK5srD2mj+2T/0n2u0w8yYArt5gy5ojKVNp+dlh5k1w6cse7qNQriysvXrG99nqf+qNA/jEzRcD5XwZlqa0/Hzi5ovh4WtdWRSq/p9/NrJceGG6mbXowHce6bJTMLcsrL3GjKk7AutSjy2zgstPwVxZWHtNnpwf+F/01pq977wKJj8E++9fdyg2AHdDWXtNntxQYZhVlyqLyXWHYYNwy8JGrM/tcWjdIRTN+bFWuLKwEWvu8u4KG0qJ+blx1iPsU2F4+9kTdx+GaKyRKwtrSbPrVEyZ9Uh6sNUwBNPEHvdcC8Blr9+25kjK5PxYK1xZ2Ij14duuAPxlOBjnx1rhA9xmZtaUKwszM2vKlYWZmTXlysLMzJryAW5rqwPfeWTdIcxXUiwlKi0/pcVjvbmysLZ6bJkV6g5hvpJiKVFp+SktHuvN3VDWVnvfeVUatqEAJcVSotLyU1o81psrC2urkj7wJcVSotLyU1o81pu7oazpv7Jbsc8HJ7ZtWza6uOyUzZWFmXWdqj9wPIZU+7gbytrqEzddxCduuqjuMKwLueyUzZWFtdUOf72ZHf56c91hWBdy2Smbu6FsxNr/vcfVHULRnB9rhSsLG7GeW2KpukMomvNjrXBlYSPWh29NB0HP3dwHOQcyGvLjA+Ht48piBGvnKbHdaI8//RYY2V+Gi8L5sVb4ALeZmTXVNS0LSbsA3wQWA34QEaP2HzyjvcVg1m7urmquK1oWkhYDvgPsCmwMfEDSxvVGZWY2enRLy2ILYGZEzAKQNAXYC7i71qjazC0Gs7KN5haIIqLuGJqStDewS0R8PD//CLBlRBzcsMwEYEJ+uhFwb5PNjgEe7kC4I4XzMzTnZ2jOz9BKzc96EbHqQDO6pWWhAab1quUiYhIwqfIGpVsiYtyiBjZSOT9Dc36G5vwMrRvz0xXHLIA5wDoNz9cG7q8pFjOzUadbKovfAxtKerWkVwL7ANNqjsnMbNToim6oiJgn6WDgl6RTZ8+KiLsWcbOVu6xGKednaM7P0JyfoXVdfrriALeZmdWrW7qhzMysRq4szMysqRFbWUg6S9KDkv7YMO04Sf+UdHu+7dYw70hJMyXdK2nneqIePgPlJ0//dM7BXZK+0jB91OdH0gUNZWe2pNsb5jk/0maSbsz5uUXSFnm6JJ2W8/MHSZvXF3nnDZKbTSXdIOlOSZdKWr5hXneUnYgYkTdgW2Bz4I8N044DDhtg2Y2BO4AlgVcDfwUWq/s11JCftwNXAUvm56s5Pwvy02f+KcAxzk+v8vMrYNf8eDdgesPjn5P+L7UVcFPd8deQm98D2+XHHwNO6LayM2JbFhFxLfBoxcX3AqZExPMR8TdgJmmIkRFrkPwcCEyMiOfzMg/m6c5PA0kC3gecnyc5P3ky0POLeQUW/BdqL+BHkdwIrChpjeGJdPgNkpuNgGvz4yuB9+THXVN2RmxlMYSDc1P4LEkr5WlrAfc1LDMnTxttXgtsI+kmSddIenOe7vz0tg3wQET8JT93fpJDgK9Kug/4GnBknu78wB+Bd+TH72XBn4y7JjejrbI4A3gNsBkwl9SVABWGExklFgdWInUV/A8wNf+Kdn56+wALWhXg/PQ4EPhcRKwDfA44M093flLX00GSZgDLAS/k6V2Tm674U167RMQDPY8lfR+4LD/1cCLJHOCiSJ2pN0t6mTTgmfOTSVoceDfwpobJzk+yH/DZ/PgnwA/y41Gfn4j4E7ATgKTXAj3D0nZNbkZVy6JPP+m7SE1DSEOH7CNpSUmvBjYEbh7u+ArwM2B7mF+gX0kaGdP5WWBH4E8RMadhmvOT3A9slx9vD/R0000D9s1nRW0FPBERc+sIsC6SVsv3rwCOBr6bZ3VN2RmxLQtJ5wPjgTGS5gDHAuMlbUZq5s0GPgkQEXdJmkq6PsY84KCIeKmOuIfLIPk5Czgrn/L3ArBfbmU4P3BsRJxJGpessQvK5WdB+fkE8M3c+nqOBZcMuIJ0RtRM4Bngo8Me8DAaJDfLSjooL3IR8EPorrLj4T7MzKypUdUNZWZmC8eVhZmZNeXKwszMmnJlYWZmTbmyMDOzplxZ2KgmaWzfkXcrrPMpSfs2WWZ/Sd8eZN5RQ6wnSb9uHJV0YUm6qmFIG7NF4srCrEUR8d2I+NEibGLQyoL0f4Q7IuLJRdh+j3OA/27DdsxcWZgBi0n6fr6Gx68kLQ0g6TWSfiFphqTfSnpdnn6cpMPy4zfngSlvkPTVPq2UNfP6f1G+NoikicDS+ZoP5w0Qy4eAS3qeSNo3b/8OSefkaZMlnSHpN5JmSdouD4x5j6TJDduaRhrHymyRubIwS0MsfCciNgEeZ8Hw0ZOAT0fEm4DDgNMHWPeHwKciYmug7z9vNwPeD7wBeL+kdSLiCODZiNgsIj40wPbeCswAkLQJ8EVg+4jYlAXjLkEa8HF70oB9lwKnApsAb8ijFBARjwFLSlqleirMBjZih/swa8HfIqLnqnczgLGSlgXeAvwkDbwLpAvUzCdpRWC5iLg+T/oxsEfDIldHxBN52buB9eg9HPVAVo6Ip/Lj7YELI+JhgIhovEbCpRERku4kDZd+Z97PXcBYoOf1PAisCTzSZL9mQ3JlYQbPNzx+CVia1Op+PCI2G2K9gYaXHmq7VT5v8yS9IiJeztsfbDyenm2/3Gc/L/fZz1LAsxX2azYkd0OZDSAfYP6bpPfC/LOUNu2zzGPAU3kkVUiDDFbxoqQlBpl3L7B+fnw18L6ebiRJK7fyGvK1SF5FGjTTbJG4sjAb3IeAAyTdAdxFugRmXwcAkyTdQGoJPFFhu5OAPwxygPty0oilRMRdwEnANTmGr7cY/5uAGyNiXovrmfXjUWfNFoGkZSPi3/nxEcAaEfHZJqsNtb01SNer/q82xPZNYFpEXL2o2zLzMQuzRbO7pCNJn6W/A/svysYiYm4+jXf5NvzX4o+uKKxd3LIwM7OmfMzCzMyacmVhZmZNubIwM7OmXFmYmVlTrizMzKyp/w8j8+IuOQ8oFQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = plt.hist(pop_heights, bins=30)\n",
    "_ = plt.xlabel('height (cm)')\n",
    "_ = plt.ylabel('number of people')\n",
    "_ = plt.title('Distribution of heights in entire town population')\n",
    "_ = plt.axvline(172, color='r')\n",
    "_ = plt.axvline(172+5, color='r', linestyle='--')\n",
    "_ = plt.axvline(172-5, color='r', linestyle='--')\n",
    "_ = plt.axvline(172+10, color='r', linestyle='-.')\n",
    "_ = plt.axvline(172-10, color='r', linestyle='-.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, 50000 people is rather a lot to chase after with a tape measure. If all you want to know is the average height of the townsfolk, then can you just go out and measure a sample to get a pretty good estimate of the average height?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def townsfolk_sampler(n):\n",
    "    return np.random.choice(pop_heights, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's say you go out one day and randomly sample 10 people to measure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed(47)\n",
    "daily_sample1 = townsfolk_sampler(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEWCAYAAACe8xtsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZgcVbnH8e+PsMsWSETIQkARAZWAEVSURVDBm4heUYMgxKtGuYKCuAByAZGLiLjhAkSMURACgmiAKAIaQSJLCIGwagjxEhIh7PsSeO8f5ww0k+6ZrsnUVHfn93meerq7llPv1Jw+p2t7SxGBmZlZEStVHYCZmbUfdx5mZlaYOw8zMyvMnYeZmRXmzsPMzApz52FmZoW582gTkk6T9D/9VNZISU9IGpQ/z5D06f4oO5f3B0kH9Fd5BdZ7vKQHJP27zrRdJC3sY7nvknRnk/P2eT1NlN1vdaAKko6VdFY/l3mkpDP6s0xrjjuPFiBpgaSnJT0u6RFJMyV9TtJL/5+I+FxEfLPJsnbvaZ6I+L+IWCsiXuiH2JdpECJiz4j45fKWXTCOEcBhwFYR8Zr+LDsiroqILfqjLElTJB3fxziaqgMrkog4ISL67YcPgKRdJf1F0qOSFtSZPipPf0rSHb193zqVO4/WMS4i1gY2AU4Evgb8vL9XImnl/i6zRWwCPBgR91cdiLW9J4HJwFcaTD8HuBHYAPg6cL6koQMUW+uICA8VD8ACYPdu47YHXgTemD9PAY7P74cAFwOPAA8BV5F+CJyZl3kaeAL4KjAKCOBTwP8BV9aMWzmXNwP4FnAd8Cjwe2D9PG0XYGG9eIE9gOeA5/P6bqop79P5/UrAUcC/gPuBXwHr5mldcRyQY3sA+HoP22ndvPySXN5Rufzd89/8Yo5jSp1ldwEWkvZO7gcWA5+smb4acHKO4z7gNGCNetsA2I7UeDwO/AY4t+Z/03A9wMS8rZ7LcV6Ux38NuDeXdyewW4O/f0oz62mw7ARgfl7H3cC+efxrgT8DD+bt/2tgvW7/668AN5Ma1Z8DGwJ/yGVdDgzu9v+cCCzKMR1WU9axwFk1n98GzCTV45uAXXqIv+42qi0T+HHerl3DUuDYPG1j4AJS3bkb+EIT38vdgQXdxr0eeBZYu2bcVcDnqm5HBnrwnkeLiojrSI3Du+pMPixPG0r6Ih+ZFolPkBq/cZEOS51Us8zOwJbA+xqscn/gv0hfsqXAKU3E+EfgBODcvL5t6sw2IQ+7ApsBa5G+5LXeCWwB7AYcLWnLBqv8EakD2Sz/PfuTGszLgT2BRTmOCQ2Wf01efhipM/2JpMF52rdJDcNo4HV5nqO7FyBpVeBCUkO+PulX6IeaWU9ETCI1ziflOMdJ2gI4CHhrpD3P95Ea7Gb09PfUxvwq0v9zz7yOdwBzuiaTfjhsTKofI0gNcq0PA+8hbZ9xpI7jSNKPmJWAL3Sbf1dgc+C9wOH1DutIGgZcAhxP2o5fBi6o9wu+2W0UEQfl7boWqU49DPw+H/69iNRBDSPVs0MkNfou9GRrYH5EPF4z7qY8foXizqO1LSJ9sbp7HtgI2CQino90TL63JGXHRsSTEfF0g+lnRsQtEfEk8D/AR7tOqC+nfYHvRcT8iHgCOAIY3+3w2Tci4umIuIn0RVymE8qxfAw4IiIej4gFwHeBTxSI5XnguLzNppN+nW4hScBngEMj4qHcMJwAjK9TxtuAlYFTcjm/Je2x9bqeBjG9QNrr2UrSKhGxICLuWp6/p8G8LwJvlLRGRCyOiFsBImJeRFwWEc9GxBLge6SOudaPIuK+iLiX9Cv72oi4MSKeJXWk23ab/xu5rs0FfgHsUyee/YDpETE9Il6MiMuAWcD768xbaBvlDuh3wMERcSPwVmBoRBwXEc9FxHzgZ9T///ZmLdLeea1HgbX7UFZbc+fR2oaRDkt19x1gHvAnSfMlHd5EWfcUmP4vYBXSL8vltXEur7bslUl7TF1qr456ivQF7W4IsGqdsoYViOXBiFhaZ11DgTWBG/IFC48Af8zju9sYuLdbZ9192zZazzIiYh5wCOnX/v2SpkraeDn/nu7reJLU8X4OWCzpEklvAJD06rzOeyU9BpzFsv/3+2reP13nc/d1dq9L9f6eTYCPdG3vvM3fSfpR1D3+preRpFWA84GzI2Jqzbo27rauI3llHWzWE8A63catQzqctkJx59GiJL2V1DD+rfu0/Mv7sIjYjHQY4UuSduua3KDI3vZMRtS8H0n6VfsA6Tj3mjVxDeKVjWpv5S4ifXlry17KKxugZjyQY+pe1r0Fy2lU9tPA1hGxXh7WzYc/ulsMDMt7K11G1JmvkWW2V0ScHRHvJP1tQTqE1q8i4tKIeA+pcb6D9Msb0iGrAN4cEeuQ9ghUv5Smda9Li+rMcw9pb3e9muFVEXFig/ib3UY/IjXkR3Vb193d1rV2RNTby+nNrcBmkmr3NLbJ41co7jxajKR1JI0FppJOBM6tM89YSa/LDdhjpN36rstu7yOdEyhqP0lbSVoTOA44P9KlvP8AVpf0H/lX3VGkQwhd7gNG1V5W3M05wKGSNpW0Fi+fI1naYP66ciznAf8raW1JmwBfIv1SXi4R8SKpMf2+pFdDOibf4Jj430nb+iBJK0vai3RxQ7Ne8f+RtIWkd0taDXiG1Ikt9yXUtSRtKOkD+dzHs6Rfz13rWDt/fiSfh2h0hVER/yNpTUlbA58kXVDQ3VnAOEnvkzRI0ur5HpnhdeJvahtJ+izpkNvH8/+0y3XAY5K+JmmNvL435h9oy5C0kqTVSXvfyrGtChAR/yCdLzomj/8Q8GbSyfgVijuP1nGRpMdJv5K+Tjr2/MkG825OusrlCVJj9tOImJGnfQs4Ku+ef7nA+s8knQT+N7A6+SRoRDwK/DdwBulX/pOkk/VdfpNfH5Q0u065k3PZV5KucnkGOLhAXLUOzuufT9ojOzuX3x++RjoUeE0+fHM5dc4fRMRzwH+STlA/QvqlfjGpUW7Gz0nH7h+R9DtSR3wiae/n38CrSYdU+tNKpIssFpEOg+5M+p8CfIN09dijpBPYv+2H9f2VtC2vAE6OiD91nyEi7gH2Iv2tS0j1/ivUb5Oa3Ub7kDrmRUo3wT4h6cj8w2Mc6WKIu3M5Z5AuNqhnJ1IHNZ205/Q0UPs3jAfGkE7Inwjsnc8XrVDU+3lWM+uJpGuB0yLiF1XHUiVJo0iN8ypF9yyt/XjPw6wgSTtLek0+bHUA6bDFH6uOy2wgderdxmZl2oJ0/mUt4C7SYYvF1YZkNrB82MrMzArzYSszMyusow5bDRkyJEaNGlV1GJ3vzpydfIt+STTb+by9ivH2GlA33HDDAxFROLFjR3Ueo0aNYtasWcUXnDkzvb7jHf0bUKfaZZf0OmNGlVG0D2+vYry9ilnO9kvSv3qfa1kd1Xn02ZH5knFXVjNrNxW1X+48rLijjup9HrO+cv1qC+48rLjdV8gHp9lAcf1qC77ayoqbMycNZmVw/WoL3vOw4g45JL36HJGVwfWrLXjPw8zMCiut85A0QtJfJN0u6VZJX6wzjySdImmepJslbVcz7QBJ/8zDAWXFaWZmxZV52GopcFhEzM4PTrlB0mURcVvNPHuS0otvDuwAnArsIGl94BhS2uPIy06LiIdLjNfMzJpU2p5Hfk7y7Pz+ceB2ln1k6F7AryK5BlhP0kakB9xflp8n/TBwGbBHWbGamVkxA3LCPOf53xa4ttukYbzyeccL87hG4+uVPRGYCDBy5Mi+BXjCCX1brk2NOvyS5Vp+u+FjAZjdoJwFJ/7HcpXfjOX9G3ozEH+DNTAA38eOqj8VtV+ldx750aMXAIdExGPdJ9dZJHoYv+zIiEnAJIAxY8b0LUWw05IUMnv4llWHYJ3M38diKtpepV5tlZ95fQHw64io93jLhcCIms/DSY/KbDS+HDNnvpwfxnq13cLb2W7h7VWHYZ3K38diKtpepe15SBLpec23R8T3Gsw2DThI0lTSCfNHI2KxpEuBEyQNzvO9FziirFid26qYr175SwDGf/zEiiOxjuTvYzEdmNtqR+ATwFxJXbeLHkl6oDwRcRrpAfPvB+YBTwGfzNMekvRN4Pq83HER8VBpkZ5+emlFm5mVqqL2q7TOIyL+Rv1zF7XzBPD5BtMmA5NLCG1Zfm6AmbWritov32EOcNFFaTAzazcVtV/ObQXw3e+m13Hjqo3DzKyoitovdx5W2HG7Taw6BOtkP/hB1RFYE9x5WGG3bbhZ1SFYJxs9uuoIrAk+52GF7bhgDjsu8PMWrCSXX54Ga2ne87DCDp45FYCrR/kXopXg+OPTq58o2NK852FmZoW58zAzs8LceZiZWWHuPMzMrDCfMAfntiroyPcdVHUI1sn8fSym03JbtRXntipk/gbDqw7BOpm/j8U4t1WFnNuqkN3mXctu87o/FNKsn/j7WIxzW1XIua0K+cx1FwJwxet2qDgS60j+Phbj3FYVOv/8qiMwM+ubitovdx4AQ4ZUHYGZWd9U1H6V+RjaycBY4P6IeGOd6V8B9q2JY0tgaH6K4ALgceAFYGlEjCkrTgCmTEmvEyaUuhozs35XUftV5gnzKcAejSZGxHciYnREjCY9n/yv3R41u2ueXm7HAWnjd/0DzMzaSUXtV5mPob1S0qgmZ98HOKesWKx/HTr2sKpDsE525plVR2BNqPych6Q1SXsotXeeBfAnSQGcHhGTKgnO6lq8ztCqQ7BONmJE1RFYEyrvPIBxwNXdDlntGBGLJL0auEzSHRFxZb2FJU0EJgKMHDmy/GiNsbenf8XFW+5UcSTWkc49N71+7GPVxmE9aoWbBMfT7ZBVRCzKr/cDFwLbN1o4IiZFxJiIGDN0qH8RD4T9bpzOfjdOrzoM61SnnpoGa2mVdh6S1gV2Bn5fM+5Vktbueg+8F7ilmgjNzKyeMi/VPQfYBRgiaSFwDLAKQESclmf7EPCniHiyZtENgQsldcV3dkT8saw4zcysuDKvttqniXmmkC7prR03H9imnKjMzKw/tMI5DzMzazOtcLVV9ZzbqpADP3hE1SFYJ/P3sRjntqqQc1sV8vCa61YdgnUyfx+LqWh7+bAVOD1JQXvPvZy9515edRjWqfx9LKai7eXOA1xZC3LnYaXy97GYTstt1VZmzKg6AjOzvqmo/fKeh5mZFebOA+Dkk9NgZtZuKmq/3HkAXHxxGszM2k1F7ZfPeVhhEz5ybNUhWCeb7qSb7cCdhxX2zCqrVx2CdbI116w6AmuCD1tZYfvNvoT9Zl9SdRjWqX760zRYS3PnYYWNveMqxt5xVdVhWKc677w0WEtz52FmZoW58zAzs8LceZiZWWHuPMzMrLAyH0M7GRgL3B8Rb6wzfRfSs8vvzqN+GxHH5Wl7AD8EBgFnRMSJZcUJOLdVQeM/Xu6/w1Zw/j4W04G5raYAe/Qyz1URMToPXR3HIOAnwJ7AVsA+krYqMU4zMyuotM4jIq4EHurDotsD8yJifkQ8B0wF9urX4LpzbqtCPnPtb/nMtb+tOgzrVP4+FrOC5rZ6u6SbJP1B0tZ53DDgnpp5FuZxdUmaKGmWpFlLlizpWxR//3sarCm73XUdu911XdVhWKdyrrliKmq/qkxPMhvYJCKekPR+4HfA5oDqzBuNComIScAkgDFjxjScr0cXXNCnxczMKldR+1XZnkdEPBYRT+T304FVJA0h7WmMqJl1OLCoghDNzKyByjoPSa+RpPx++xzLg8D1wOaSNpW0KjAemFZqMEcckQYzs3ZTUftV5qW65wC7AEMkLQSOAVYBiIjTgL2BAyUtBZ4GxkdEAEslHQRcSrpUd3JE3FpWnIDPdxT0zMqrVR2CdbI11qg6gvZSUfvVa+chaUPgBGDjiNgzXzb79oj4eU/LRcQ+vUz/MfDjBtOmA07q36ImfPQbVYdgnewPf6g6AmtCM4etppD2AjbOn/8BHFJWQGZm1vqa6TyGRMR5wIsAEbEUeKHUqKylHXz1ORx89TlVh2Gd6pvfTIO1tGY6jyclbUC+XFbS24BHS43KWtqO/7qJHf91U9VhWKe64oo0WEtr5oT5l0hXO71W0tXAUNLJbjMzW0H12nlExGxJOwNbkG7guzMini89MjMza1kNOw9J/9lg0uslERFObmRmtoLqac9jXA/TAuiczmODDaqOoK08vMY6VYdgnczfx2Iq2l4NO4+I+ORABlIp57Yq5MAPHVl1CNbJ/H0splVzW0naQNIpkmZLukHSD/PVV2ZmtoJq5lLdqcAS4MOkq6yWAOeWGdSAc26rQr761yl89a9Tqg7DOpW/j8W0cG6r9SOi9o6d4yV9sKyAKvHgg1VH0Fa2u/eOqkOwTuZcc8VU1H4103n8RdJ44Lz8eW/gkvJCqsCkSVVHYGbWNxW1X80ctvoscDbwXB6mAl+S9Likx8oMzszMWlMzNwmuPRCBVGrixPTqPRAzazcVtV9NPc9D0geAnfLHGRHRWQ8Y/sc/qo6grSxee0jVIVgnGz686gjaS0XtVzPP8zgReCvw6zzqi5LeGRGHlxqZtaxDx3256hCsk511VtURWBOa2fN4PzA6Il4EkPRL4EbAnYeZ2Qqq2WeYr1fzft1mFpA0WdL9km5pMH1fSTfnYaakbWqmLZA0V9IcSbOajNEGyNGXT+Loy31+yEpyyCFpsJbWzJ7Ht4AbJf2FlFV3J6CZO1KmkB4z+6sG0+8Gdo6IhyXtCUwCdqiZvmtEPNDEemyAbXX//KpDsE42Z07VEVgTmrna6hxJM0jnPQR8LSL+3cRyV0oa1cP0mTUfrwF8lszMrE00k9tKwG6k8x6/B1aVtH0/x/EpoPap9wH8KefSmthLfBMlzZI0a8mSJf0clpmZ1dPMOY+fAm8H9smfHwd+0l8BSNqV1Hl8rWb0jhGxHbAn8HlJO9VdGIiISRExJiLGDB06tL/CMjOzHjRzzmOHiNhO0o0A+RzFqv2xcklvBs4A9oyIlxK0RMSi/Hq/pAuB7YEr+2Oddb3+9aUV3Ynmrz+s6hCsk/n7WExF26uZzuN5SYNIh5KQNBR4cXlXLGkk6YFSn4iIf9SMfxWwUkQ8nt+/FzhuedfXI99ZXsiRexxcdQjWyfx9LKai7dVM53EKcCGwoaT/JSVGPKq3hSSdA+wCDJG0EDgGWAUgIk4DjgY2AH6aTquwNCLGABsCF+ZxKwNnR8Qfi/1ZZmZWpmautvq1pBtIJ80BPhgRtzex3D69TP808Ok64+cD2yy7RImc26qQE/74I8B7IFYSfx+LaeXcVsCaQNehqzXKC6cifmZyIZs9dG/VIVgnc665YlrtGeZdJB0NfAS4gHSfxy8k/SYiji87uAHzrW9VHYGZWd9U1H41s+exD7BtRDwDLyVKnA10TudhZmaFNHOfxwJg9ZrPqwF3lRJNVT784TSYmbWbitqvZvY8ngVulXQZ6ZzHe4C/SToFICK+UGJ8A8PPMC/ktldvVnUI1slGj646gvbSws8wvzAPXWaUE4q1i+N27zFjjNny+cEPqo7AmtDMpbq/HIhAzMysfTT7PA+zl3z/opP5/kUnVx2Gdar99kuDtbRm7/Mwe8lGj/sxK1aihQurjsCa0HDPQ9KZ+fWLAxeOmZm1g54OW71F0ibAf0kaLGn92mGgAjQzs9bT02Gr04A/ApsBN5DuLu8SebyZma2AGnYeEXEKcIqkUyPiwAGMaeC9/e1VR9BWZg97Q9UhWCfz97GYirZXM5fqHihpG+BdedSVEXFzuWENMOe2KuSknSdUHYJ1Mn8fi6loezXzDPMvAL8GXp2HX0tyLm4zsxVYM5fqfpr0KNonASR9G/g78KMyAxtQXXlhLrig2jjaxKkXngDAgR86suJIrCP5+1hMRdurmc5DwAs1n1/glSfP25+PsRYy+OnHqg7BOplzzRXTquc8gF8A10rqym/1QeDnzRQuaTIwFrg/It5YZ7qAHwLvB54CJkTE7DztAF5+3O3xpaZJ+fKXSyvazKxUFbVfvZ7ziIjvAZ8EHgIeBj4ZEc1mLpsC7NHD9D2BzfMwETgVIN9HcgywA7A9cIykwU2u08zMStZUepK8NzC7aOERcaWkUT3Mshfwq4gI4BpJ60naCNgFuCwiHgLI6eD3AM4pGkNTdtklvc6YUUrxZmalqaj9qjq31TDgnprPC/O4RuOXIWkiaa+FkSNHlhOlvcLVm2zT4/RRh18yQJFYVcr8Hx/8wnAAftTG9ajs78CCE/+j1PKbUXXnUe/Ee/QwftmREZOASQBjxoypO4/1rx/tuE/VIVgHc/1qDz2e85A0SNLlJa5/ITCi5vNwYFEP483MrAX02HlExAvAU5LWLWn904D9lbwNeDQiFgOXAu/NCRkHA+/N46wFTDnvGKacd0zVYViHcv1qD80ctnoGmJtPWj/ZNbKZZ5dLOod08nuIpIWkK6hWycufBkwnXaY7j3Sp7ifztIckfRO4Phd1XNfJc6ve6kufrToE62CuX+2hmc7jkjwUFhE9HrzMV1l9vsG0ycDkvqzXzMzK1dQzzCWtAYyMiDsHICYzM2txzSRGHAfMIT3bA0mjJU0rOzAzM2tdzRy2OpZ0l/cMgIiYI2nTEmMaeGPHVh1BW7nitdtXHYJ1MNevgipqv5rpPJZGxKMpDdVLOut+Cue2KuRnO/xn1SFYB3P9Kqii9quZzuMWSR8HBknaHPgCMLPcsMzMrJX1es4DOBjYGniWlFvqMeCQMoMacLvs8nJ+GOvV1LMPZ+rZh1cdhnUo16+CKmq/mrna6ing6/khUBERj5cf1gCbMKHqCMzM+qai9qvXzkPSW0n3W6ydPz8K/FdE3FBybAPHnYeZtatW7TxID37674i4CkDSO0kPiHpzmYENqAceSK9DhlQbh5lZURW1X810Ho93dRwAEfE3SZ116GrvvdOrn+dhZu2movarYechabv89jpJp5NOlgfwMfI9H7ZiuvgN76o6BOtgrl/toac9j+92+1yb5rKz7vOwQs7arvoH0Vjncv1qDw07j4jYdSADsfax+vPPAPDMKqtXHIl1Itev9tDM1VbrAfsDo2rnbyYlu3WmKb85FoDxHz+x2kCsI7l+tYdmTphPB64B5gIvlhuOmZm1g2Y6j9Uj4kulR2JmZm2jmfQkZ0r6jKSNJK3fNZQemZmZtaxmOo/ngO8AfwduyMOsZgqXtIekOyXNk7RMshpJ35c0Jw//kPRIzbQXaqb5+SFmZi2kmcNWXwJeFxEPFClY0iDgJ8B7gIXA9ZKmRcRtXfNExKE18x8MbFtTxNMRMbrIOvvM6UkKOf9Nu1cdgnUw16+CWjg9ya3AU30oe3tgXkTMB5A0FdgLuK3B/PvwyntJBo47j0L85bYyuX4V1MKdxwvAHEl/IaVlB5q6VHcYcE/N54XADvVmlLQJsCnw55rRq0uaBSwFToyI3zVYdiIwEWDkyJG9hNSAc1sVMvipRwF4eM11K47EOpHrV0EtnNvqd3koSnXGNbozfTxwfkS8UDNuZEQskrQZ8GdJcyPirmUKjJgETAIYM2ZM3+58d26rQk793bcAX4dv5XD9KqjVclt1iYhf9rHshcCIms/DgUUN5h0PfL7behfl1/mSZpDOhyzTefSLww4rpVgzs9JV1H41c4f53dTZY4iIzXpZ9Hpgc0mbAveSOoiP1yl/C2Aw6WqurnGDgaci4llJQ4AdgZN6i7XPxo0rrWgzs1JV1H41c9hqTM371YGPAL3e5xERSyUdBFwKDAImR8Stko4DZkVE1+W3+wBTI6K2g9oSOF3Si6TLiU+svUqr3915Z3rdYovSVmFmVoqK2q9mDls92G3UDyT9DTi6iWWnk9Kb1I47utvnY+ssNxN4U2/l95vPfja9+pyHmbWbitqvZg5bbVfzcSXSnsjapUVkLe+sbd9fdQjWwVy/2kMzh61qn+uxFFgAfLSUaKwtXLzlTlWHYB3M9as9NHPYys/1sFfY6LElACxeZ2jFkVgncv1qD80ctloN+DDLPs/juPLCslb2/YvTzqivw7cyuH61h2YOW/0eeJSUEPHZXuY1M7MVQDOdx/CI2KP0SMzMrG00k5J9pqSBu2zWzMxaXjN7Hu8EJuQ7zZ8l5ayKiHhzqZGZmVnLaqbz2LP0KKrm3FaF/Gz7D1UdgnUw16+CWjW3VUT8ayACqZRzWxVyxevqZtY36xeuXwVV1H41c86j891558v5YaxXmz24kM0eXFh1GNahXL8Kqqj9auawVedzbqtCTrj0x4Cvw7dyuH4V1Kq5rVYIJ5xQdQRmZn1TUfvlzgPgHe+oOgIzs76pqP3yOQ+AmTPTYGbWbipqv7znAXDkkenV5zzMrN1U1H6587DCfvSO8VWHYB3M9as9lHrYStIeku6UNE/S4XWmT5C0RNKcPHy6ZtoBkv6ZhwPKjNOKuXrUaK4eNbrqMKxDuX61h9L2PCQNAn4CvAdYCFwvaVqdZ5GfGxEHdVt2feAY0lMLA7ghL/twWfFa87a6bz4At224WcWRWCdy/WoPZe55bA/Mi4j5EfEcMBXYq8ll3wdcFhEP5Q7jMsCZfVvE0VdM4ugrJlUdhnUo16/2UGbnMQy4p+bzwjyuuw9LulnS+ZJGFFwWSRMlzZI0a8mSJf0Rt5mZ9aLMzkN1xkW3zxcBo3KG3suBXxZYNo2MmBQRYyJizNChfmylmdlAKLPzWAiMqPk8HFhUO0NEPBgRXU8n/BnwlmaXNTOz6pTZeVwPbC5pU0mrAuOBabUzSNqo5uMHgNvz+0uB90oaLGkw8N48zszMWkBpV1tFxFJJB5Ea/UHA5Ii4VdJxwKyImAZ8QdIHgKXAQ8CEvOxDkr5J6oAAjouIh8qK1bmtijlpJ185beVx/SqoE3NbRcR0YHq3cUfXvD8COKLBspOByWXG9xLntipk9vAtqw7BOpjrV0HObVUh57YqZLuFt7Pdwtt7n9GsD1y/CnJuqwo5t1UhX70yXRTn5y1YGVy/CnJuqwqdfnrVEZiZ9U1F7Zc7D4Attqg6AjOzvqmo/fI5D4CLLkqDmVm7qaj98p4HwHe/m17Hjas2DjOzoipqv9x5WGHH7Tax6hCsg7l+tQd3HlaYU2VbmWfRpxcAAAzdSURBVFy/2oPPeVhhOy6Yw44L5lQdhnUo16/24D0PK+zgmVMB/LQ3K4XrV3vwnoeZmRXmzsPMzApz52FmZoW58zAzs8J8whyc26qgI993UNUhWAdz/SrIua0q5NxWhczfYHjVIVgHc/0qyLmtKuTcVoXsNu9adpt3bdVhWIdy/SqoE3NbSdoD+CHpMbRnRMSJ3aZ/Cfg06TG0S4D/ioh/5WkvAHPzrP8XER8oLVDntirkM9ddCMAVr9uh4kisE7l+FdRpua0kDQJ+ArwHWAhcL2laRNxWM9uNwJiIeErSgcBJwMfytKcjYmDuEjr//AFZjZlZv6uo/SrzsNX2wLyImB8RzwFTgb1qZ4iIv0TEU/njNUA1BzuHDEmDmVm7qaj9KrPzGAbcU/N5YR7XyKeAP9R8Xl3SLEnXSPpgo4UkTczzzVqyZEnfIp0yJQ1mZu2movarzHMeqjMu6s4o7QeMAXauGT0yIhZJ2gz4s6S5EXHXMgVGTAImAYwZM6Zu+b3q2vATJvRpcTOzylTUfpXZeSwERtR8Hg4s6j6TpN2BrwM7R8SzXeMjYlF+nS9pBrAtsEznYQPv0LGHVR2CdTDXr/ZQZudxPbC5pE2Be4HxwMdrZ5C0LXA6sEdE3F8zfjDwVEQ8K2kIsCPpZLq1gMXrDK06BOtgrl/tobTOIyKWSjoIuJR0qe7kiLhV0nHArIiYBnwHWAv4jSR4+ZLcLYHTJb1IOi9zYrertKxCY2+/EoCLt9yp4kisE7l+tYdS7/OIiOnA9G7jjq55v3uD5WYCbyozNuu7/W5M/1J/ua0Mrl/twXeYm5lZYe48zMysMHceZmZWmDsPMzMrzCnZwbmtCjrwg0dUHYJ1MNevgipqv9x5gPNaFfTwmutWHYJ1MNevgipqv3zYCpzbqqC9517O3nMvrzoM61CuXwVV1H658wB3HgX5y21lcv0qqAMTI7aPGTOqjsDMrG8qar+852FmZoW58wA4+eQ0mJm1m4raL3ceABdfnAYzs3ZTUfvlcx5W2ISPHFt1CNbBXL/agzsPK+yZVVavOgTrYK5f7cGHrayw/WZfwn6zL6k6DOtQrl/twZ2HFTb2jqsYe8dVVYdhHcr1qz248zAzs8JK7Twk7SHpTknzJB1eZ/pqks7N06+VNKpm2hF5/J2S3ldmnGZmVkxpnYekQcBPgD2BrYB9JG3VbbZPAQ9HxOuA7wPfzstuBYwHtgb2AH6ayzMzsxZQ5p7H9sC8iJgfEc8BU4G9us2zF/DL/P58YDdJyuOnRsSzEXE3MC+XZ2ZmLaDMS3WHAffUfF4I7NBonohYKulRYIM8/ppuyw6rtxJJE4GJ+eMTku7sc8RSoylDgAf6XG55Konr7V1vvj220Sytur2gydj07RLW3Lh+QetuswGPq4n6BSv49qpbP/tevzbpSwxldh71/pJocp5mlk0jIyYBk4qFVoykWRExpsx19IXjKq5VY3NcxTiuYsqIq8zDVguBETWfhwOLGs0jaWVgXeChJpc1M7OKlNl5XA9sLmlTSauSToBP6zbPNOCA/H5v4M8REXn8+Hw11qbA5sB1JcZqZmYFlHbYKp/DOAi4FBgETI6IWyUdB8yKiGnAz4EzJc0j7XGMz8veKuk84DZgKfD5iHihrFibUOphseXguIpr1dgcVzGOq5h+j0vph76ZmVnzfIe5mZkV5s7DzMwKW+E6D0mTJd0v6ZZu4w/OqVBulXRSHvceSTdImptf392gzGMl3StpTh7eX3JcoyQ9XbO+0xqUub6kyyT9M78OLjmufWtimiPpRUmj65RZyvbKqW66ylwgaU7NtF7T3eSLO67N2+vcfKFHaXFVXb96iKvS+tVDXANWv3qIbbSka3K5syRtn8dL0im5jt0sabsGZb4l/7/n5fl7vDmjH+LaN8dzs6SZkrZpUOYUSXfXbLNltusyImKFGoCdgO2AW2rG7QpcDqyWP786v24LbJzfvxG4t0GZxwJfHsC4RtXO10OZJwGH5/eHA98uM65uy70JmD+Q26vb9O8CR+f3WwE3AasBmwJ3AYPqLHMeMD6/Pw04sOS4Kq1fPcRVaf1qFNdA1q9GsQF/AvbM798PzKh5/wfSPWpvA65tUOZ1pPsgleffs+S43gEMzu/37CGuKcDeReJY4fY8IuJK0pVdtQ4EToyIZ/M89+fXGyOi6/6SW4HVJa1WdVwF1KZ/+SXwwQGMax/gnKLrW864gPQrEPhozfp7TXeTl3k3KU0O9O/2qhtXC9SvunEVUFb9aiauUusXNIwtgHXy+3V5+f6zvYBfRXINsJ6kjWoXzJ/XiYi/R2qxf0X/bbO6cUXEzIh4OI+/hnTPXL9Y4TqPBl4PvCsfsvirpLfWmefDwI1dDWYdB+Vdw8l92X3vQ1ybSroxj39Xg+U3jIjFAPn11QMQV5eP0fOXu4zt1eVdwH0R8c/8uV6qnO7pbjYAHomIpT3M099x1aqifvUUV5X1q6e4ulRVvw4BviPpHuBk4Ig8vpk6NiyP72me/o6r1qdIezuN/G/eZt9v5keMO49kZWAwaXfzK8B5tcciJW1Nyvj72QbLnwq8FhgNLCbtapcZ12JgZERsC3wJOFvSOo2L6Xe9ba8dgKci4pYGy5e1vbp0/1W6PKly+lPdX8sV1q9GcVVdvxrFBVRevw4EDo2IEcChpHvVoPo61iiutGJpV1Ln8bUGyx8BvAF4K7B+D/O9xJ1HshD4bd7lvA54kZRIDEnDgQuB/SPirnoLR8R9EfFCRLwI/Iz+ywBcN658+OXBvO4bSMfwX19n+fu6dp3za9HDXoXiqpk+nh5+FZa4vbrS3PwncG63eHtLd/MA6VDDyj3M099xVV2/6sbVAvWr4fbKKqtfpIwYv83vf1NTdrPpmIb3Mk9/x4WkNwNnAHt1/V+7i4jF+fv8LPALmthm7jyS35GOdyPp9cCqwAOS1gMuAY6IiKsbLdzt2OaHgEa/iPorrqHKzzeRtBkpfcv8OsvXpn85APh9mXHlzysBHyGl4K+rxO0FsDtwR0TUHh7oNd1NPgb9F1KaHOjf7VU3rhaoX43iqrp+1Y0rx1N1/VoE7JzfvxvoOqQ2DdhfyduAR7sO6XXJnx+X9La8p74//bfN6sYlaSSpU/lERPyj0cI1PwJEOg/T+zYrcna9EwbSL5bFwPOkXwKfIjV+Z+UNNht4d573KOBJYE7N0HXF0xnAmPz+TGAucDOpEm1UclwfJp1gvSmPH1dTTm1cGwBX5Ip0BbB+mXHl+XcBrqlTTunbK16+auRzdeb/OukX9J3UXOECTOflK542I3Uq80i/3lYrM66q61cPcVVav3r5Pw5I/eqh7r8TuCFvm2uBt+R5RXr43V153WNqyplT834M6XtzF/BjcpaPEuM6A3i4pn7NalD3/5zjvoX03V6rtzicnsTMzArzYSszMyvMnYeZmRXmzsPMzApz52FmZoW58zAzs8LcedgKRylrbKFr/yV9TtL+vcwzQdKPG0w7soflJOnP/XEXt6TLS0hfYrYMdx5mTYiI0yLiV8tRRMPOg5QF9aaIeGw5yu9yJvDf/VCOWY/cediKapCknyk9j+RPktYAkPRaSX9Uer7GVZLekMcfK+nL+f1bcwK5v0v6Tre9mI3z8v/Uy885ORFYQ+k5Cb+uE8u+1NxpLGn/XP5Nks7M46ZIOlXSXyTNl7RzTvp3u6QpNWVNI+WEMiuVOw9bUW0O/CQitgYeId1VDTAJODgi3gJ8GfhpnWV/Qbr7+e3AC92mjSZlfH0T8DFJIyLicODpiBgdEfvWKW9H0t3BXUkSv066a38b4Is18w0mpZ44FLgI+D6wNfAm5Yf3REq/vZqkDZrfFGbFrdz7LGYd6e6I6HrS4A3AKElrkR6e85uaJMGvSE2d81GtHREz86izgbE1s1wREY/meW8DNuGVqbrrWT8iHs/v3w2cHxEPAERE7XMbLoqIkDSXlKp8bl7PraQHOHX9PfcDGwN1k+CZ9Qd3Hraiqn1uxgvAGqQ98UcioqdHcPb22NDu5TbzHVsqaaVIWWBF4zTdXWW/2G09L3Zbz+rA002s16zPfNjKLMsnrO+W9BF46SqobbrN8zA5M2oeNb7J4p+XtEqDaXeSkjJCSjD40a7DTpLWL/I35KyorwEWFFnOrCh3HmavtC/wKUk3kTLL7lVnnk8BkyT9nbSn8GgT5U4Cbm5wwvwSUrZYIuJW4H+Bv+YYvlcw/reQss4u7XVOs+XgrLpmBUlaKyKeyO8PJ6X8/mIvi/VU3kak51+/px9i+yEwLSKuWN6yzHricx5mxf2HpCNI359/AROWp7CIWJwvG16nH+71uMUdhw0E73mYmVlhPudhZmaFufMwM7PC3HmYmVlh7jzMzKwwdx5mZlbY/wPAAfWyV7PTaQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = plt.hist(daily_sample1, bins=10)\n",
    "_ = plt.xlabel('height (cm)')\n",
    "_ = plt.ylabel('number of people')\n",
    "_ = plt.title('Distribution of heights in sample size 10')\n",
    "_ = plt.axvline(172, color='r')\n",
    "_ = plt.axvline(172+5, color='r', linestyle='--')\n",
    "_ = plt.axvline(172-5, color='r', linestyle='--')\n",
    "_ = plt.axvline(172+10, color='r', linestyle='-.')\n",
    "_ = plt.axvline(172-10, color='r', linestyle='-.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sample distribution doesn't look much like what we know (but wouldn't know in real-life) the population distribution looks like. What do we get for the mean?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "173.47911444163503"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(daily_sample1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And if we went out and repeated this experiment?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "daily_sample2 = townsfolk_sampler(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "173.7317666636263"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(daily_sample2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ Simulate performing this random trial every day for a year, calculating the mean of each daily sample of 10, and plot the resultant sampling distribution of the mean."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3debwkVXn/8c+XHQLIMgOyzDCgQASXUQeCioqAUZYJJuxBAUMyChFFwB/DokGDOCgowQhIAAcBBQSUVSMgEzUghJFBIDrKMsoyYZMdRIZ5fn+cc+9terrvrbq3q+t2z/f9etWruqurTj333KfrdG2nFBGYmZkBLFN3AGZmNn64UTAzs0FuFMzMbJAbBTMzG+RGwczMBrlRMDOzQW4UxgFJZ0r6bIfKmizpOUnL5vdzJP1jJ8rO5f1Q0gGdKq/Eek+Q9Lik/2vx2XaSHhxlue+WNL/gvKNeT4GyO5YDdZB0vKQLOlzmMZLO7mSZNjI3ChWTtEDSi5KelfSUpJskfVzSYN1HxMcj4l8LlrXjcPNExB8iYtWIeKUDsS/xRY+InSLivLGWXTKOScARwBYR8dpOlh0RP4uIzTtRlqTZkk4YZRyFcmBpEhEnRkTHftAASHqfpBslPS1pQYvPp+TPX5D0m5G+b/3IjUJ3TI+I1YCNgFnAUcA5nV6JpOU6XeY4sRHwREQ8Wncg1vOeB84FPtPm8+8CtwNrA8cCl0qa2KXYxoeI8FDhACwAdmyatjWwGHhjfj8bOCG/ngBcDTwF/BH4GanxPj8v8yLwHPD/gClAAAcBfwB+2jBtuVzeHOBLwK3A08AVwFr5s+2AB1vFC3wQ+DPwcl7fHQ3l/WN+vQxwHPB74FHg28Br8mcDcRyQY3scOHaYenpNXv6xXN5xufwd89+8OMcxu8Wy2wEPkvYmHgUWAh9t+HxF4OQcxyPAmcDKreoAeBtpo/As8D3g4ob/Tdv1ADNyXf05x3lVnn4U8FAubz6wQ5u/f3aR9bRZ9kDgvryO+4H98vTXAT8Bnsj1fyGwRtP/+jPAr0gby3OAdYEf5rKuB9Zs+n/OAB7OMR3RUNbxwAUN77cBbiLl8R3AdsPE37KOGssE/j3X68CwCDg+f7Y+cBkpd+4HPlnge7kjsKBp2mbAS8BqDdN+Bny87u1INwfvKdQgIm4lfenf3eLjI/JnE0lf0GPSIvER0kZteqTDQ19uWOa9wBuAD7RZ5f7AP5C+PIuA0wrE+CPgRODivL63tJjtwDy8D9gEWJX05W20LbA5sAPwOUlvaLPKr5Mahk3y37M/aUN4PbAT8HCO48A2y782L78BqZH8hqQ182cnkb7wU4HX53k+11yApBWA75M20GuRfjX+bZH1RMRZpI3ul3Oc0yVtDnwC2CrSnuIHSBviIob7expj/gvS/3OnvI53AvMGPib9IFiflB+TSBvaRrsD7yfVz3RSg3AM6cfJMsAnm+Z/H7Ap8NfAzFaHVyRtAFwDnECqxyOBy1r94i5aRxHxiVyvq5Jy6knginwY9ipSw7MBKc8Ok9TuuzCcLYH7IuLZhml35OlLDTcK9XmY9IVp9jKwHrBRRLwc6Zj3SB1UHR8Rz0fEi20+Pz8i7oqI54HPAnsNnIgeo/2Ar0bEfRHxHHA0sE/TYazPR8SLEXEH6Qu2ROOSY9kbODoino2IBcApwEdKxPIy8IVcZ9eSfk1uLknAPwGfjog/5i/8icA+LcrYBlgOOC2XczlpD2vE9bSJ6RXSXsoWkpaPiAURce9Y/p428y4G3ihp5YhYGBF3A0TEPRFxXUS8FBGPAV8lNbiNvh4Rj0TEQ6RfxbdExO0R8RKpgXxr0/yfz7l2J/AtYN8W8XwYuDYiro2IxRFxHXAbsHOLeUvVUW5YfgAcGhG3A1sBEyPiCxHx54i4D/gPWv9/R7IqaW+60dPAaqMoq2e5UajPBqTDQ82+AtwD/FjSfZJmFijrgRKf/x5YnvRLcKzWz+U1lr0caQ9nQOPVQi+QvnjNJgArtChrgxKxPBERi1qsayKwCjA3n+h/CvhRnt5sfeChpka4uW7brWcJEXEPcBjp1/mjki6StP4Y/57mdTxPalA/DiyUdI2kvwSQtE5e50OSngEuYMn/+yMNr19s8b55nc251Orv2QjYc6C+c51vS/qx0xx/4TqStDxwKfCdiLioYV3rN63rGF6dg0U9B6zeNG110mGtpYYbhRpI2oq0wft582f5l/IREbEJaXf+cEk7DHzcpsiR9iQmNbyeTPoV+jjpOPIqDXEty6s3liOV+zDpS9lY9iJevWEp4vEcU3NZD5Usp13ZLwJbRsQaeXhNPgzRbCGwQd67GDCpxXztLFFfEfGdiNiW9LcF6VBWR0XEf0bE+0kb3d+QfilDOnQUwJsjYnXSL3i1LqWw5lx6uMU8D5D2TtdoGP4iIma1ib9oHX2dtIE+rmld9zeta7WIaLVXMpK7gU0kNe4ZvCVPX2q4UegiSatL2hW4iHQC7c4W8+wq6fV5w/QMafd64PLSR0jH3Mv6sKQtJK0CfAG4NNIlq78FVpK0S/4VdhxpV37AI8CUxstnm3wX+LSkjSWtytA5iEVt5m8px3IJ8EVJq0naCDic9Mt2TCJiMWkj+TVJ60A65t3mmPPNpLr+hKTlJO1GuiigqFf9fyRtLml7SSsCfyI1TmO+VLiRpHUl/U0+t/AS6dfuwDpWy++fysf5211xU8ZnJa0iaUvgo6QT8c0uAKZL+oCkZSWtlO/x2LBF/IXqSNLHSIe+/j7/TwfcCjwj6ShJK+f1vTH/8FqCpGUkrUTaW1aObQWAiPgt6XzMv+Tpfwu8mXQSe6nhRqE7rpL0LOlXzbGkY7sfbTPvpqSrPp4jbaROj4g5+bMvAcfl3eQjS6z/fNLJ0/8DViKfPIyIp4FDgLNJv8qfJ53kHvC9PH5C0i9blHtuLvunpKs+/gQcWiKuRofm9d9H2oP6Ti6/E44iHZL7RT6Mcj0tjs9HxJ+BvyOd2H2K9Mv6atLGtohzSMfGn5L0A1IDO4u0t/J/wDqkQxudtAzp4oSHSYcj30v6nwJ8nnQ11dOkE7+Xd2B9/0WqyxuAkyPix80zRMQDwG6kv/UxUt5/htbbm6J1tC+pwX1Y6ebM5yQdk39QTCddRHB/Luds0kn6Vt5DaniuJe3pvAg0/g37ANNIJ7JnAXvk8zFLDY18DtNs6SXpFuDMiPhW3bHUSdIU0kZ3+bJ7gtZbvKdg1kDSeyW9Nh8+OoB0+OBHdcdl1i39eges2WhtTjq/sSpwL+nwwcJ6QzLrHh8+MjOzQZXuKeQOp54lXU2wKCKmSVqLdMXCFNKdi3tFxJNVxmFmZsVUuqeQG4VpEfF4w7QvA3+MiFn5xqw1I+Ko4cqZMGFCTJkypbI4LZufe5DevCOdhvY/11c5rq+umzt37uMRUapDvzrOKexG6vAL4DxSB2vDNgpTpkzhtttuK7+mm25K43e+s/yyS6PttkvjOXPqjKJ3uL7KcX2V04Htl6TfjzzXq1XdKASpu4YAvpk7DVt34MRdRCwcuKGomaQZpB4ZmTx58ujWfky+3NlJaGa9pqbtV9WNwrsi4uG84b9O0m+KLpgbkLMApk2b5rPh3XDccSPPYzZazq+eUGmjEBEP5/Gjkr5P6jLgEUnr5b2E9Uj9xdt4sONS95Ap6ybnV0+o7OY1SX8x0LFU7pflr4G7gCtJD14hj6+oKgYrad68NJhVwfnVE6rcU1gX+H7ucHI5Une3P5L0P8AlkgaeFrZnhTFYGYcdlsY+B2NVcH71hMoahfywiyUeqBIRT5CejmRmZuOM+z4yM7NBbhTMzGyQGwUzMxvU372knnhi3RH0FtdXJabMvKbtZwtm7dLFSGrm/Cqnpvrq70bB3VuU4/qyKjm/yqmpvvr78NFNNw31H2Ijc31ZlZxf5dRUX/29p+C+j8pxfVmVnF/l9GnfR/X65jfrjsDMbHRq2n71d6PgftvNrFfVtP3q73MKV12VBjOzXlPT9qu/9xROOSWNp0+vNw4zs7Jq2n71d6Ng5Zx6at0RWD9zfvUENwo2ZOrUuiOwfub86gn9fU7Byrn++jSYVcH51RO8p2BDTjghjf2ELKuC86sneE/BzMwGuVEwM7NBbhTMzGyQGwUzMxvU3yea3fdROa4vq5Lzqxz3fVQB931UjuvLquT8Ksd9H1XAfR+V4/qyKjm/ynHfRxVw30fluL6sSs6vctz3UQUuvbTuCMzMRqem7Vd/NwoTJtQdgZnZ6NS0/ervcwqzZ6fBzKzX1LT9cqNgZjYe1bT96u/DR1bO+efXHYH1M+dXT3CjYEMmTao7Autnzq+e0N+Hj6yciy9Og1kVnF89wXsKNuSMM9J4773rjcP6k/OrJ3hPwczMBrlRMDOzQW4UzMxsUOWNgqRlJd0u6er8fmNJt0j6naSLJa1QdQxmZlZMN040fwr4NbB6fn8S8LWIuEjSmcBBwBmVrNl9H5Xj+rIqOb/Kqam+Kt1TkLQhsAtwdn4vYHtg4K89D/hQZQFMmOD+j8pwfVmVnF/l1FRfVe8pnAr8P2C1/H5t4KmIWJTfPwhs0GpBSTOAGQCTJ08e3doHbhE/8MDRLb+0cX2NaMrMawZfX3TfEwDsk6ctmLVLLTH1DOdXOTXVV2V7CpJ2BR6NiLmNk1vMGq2Wj4izImJaREybOHHi6IJw30fluL6sSs6vcvqw76N3AX8jaWdgJdI5hVOBNSQtl/cWNgQeriyCOXMqK9rMrFI1bb8q21OIiKMjYsOImALsA/wkIvYDbgT2yLMdAFxRVQxmZlZOHfcpHAUcLuke0jmGcypb08knp8HMrNfUtP3qSqMQEXMiYtf8+r6I2DoiXh8Re0bES5Wt+Oqr02Bm1mtq2n65Qzwbcu21dUdg/cz51RPcKNiQVVapOwLrZ86vnuC+j2zI6aenwawKzq+e4EbBhlxySRrMquD86gluFMzMbJAbBTMzG+RGwczMBvnqI7MWGju+M1ua9Hej4L6PynF9WZWcX+X0W99HZmbWe/q7UXDfR+W4vqxKzq9yaqqv/j58dPPNdUfQWwb6WTnyyHrjsP7k/Cqnpu1XfzcKl11WdwRmZqNT0/arvw8fmZlZKf3dKBx9dBrMzHpNTduv/j585HMK5ay8ct0RWD9zfpUzXs8pSFoXOBFYPyJ2krQF8I6IqO6JaVaPH/6w7giWOu1uklswa5cuR9IFzq+eUOTw0WzgP4H18/vfAodVFZCZmdWnSKMwISIuARYDRMQi4JVKo7J6/Ou/psGsCs6vnlCkUXhe0tpAAEjaBni60qisHjfckAazKji/ekKRE82HA1cCr5P038BEYI9KozIzs1qM2ChExC8lvRfYHBAwPyJerjwyMzPruraNgqS/a/PRZpKIiMsrisnMzGoy3J7C9GE+C2D8Nwprr113BL3F9WVVcn6VU1N9tW0UIuKj3QykEu77qBzXl1XJ+VXOeO37SNLakk6T9EtJcyX9W74ayczM+kyRS1IvAh4DdidddfQYcHGVQXWM+z4qx/VlVXJ+lTOO+z5aKyIa7zg5QdKHqgqoo554ou4Ieov7irIqOb/KqWn7VaRRuFHSPsAl+f0eQG881fyss+qOwMxsdGrafhU5fPQx4DvAn/NwEXC4pGclPVNlcGZm1l1Fbl5brRuBVGLGjDT2HoOZ9Zqatl+Fnqcg6W+A9+S3cyLi6upC6qDf/rbuCHrLhhvWHYH1M+dXOTVtv4o8T2EWsBVwYZ70KUnbRsTMSiOz7rvggrojsH7m/OoJRfYUdgamRsRiAEnnAbcDbhSsJ7R7kA306cNszMag6DOa12h4/ZoiC0haSdKtku6QdLekz+fpG0u6RdLvJF0saYWyQVtFDjssDWZVcH71hCJ7Cl8Cbpd0I6mX1PcARe6oeAnYPiKek7Q88HNJPyR1xf21iLhI0pnAQcAZowvfOmrevLojsH7m/OoJI+4pRMR3gW1IHeBdTno+80UFlouIeC6/XT4PAWwPXJqnnwf0xo1wZmZLgSJ9HwnYgXRe4QpgBUlbFylc0rKS5gGPAtcB9wJP5Ud6AjwIbNBm2RmSbpN022OPPVZkdWZmNkZFzimcDrwD2De/fxb4RpHCI+KViJgKbAhsDbyh1Wxtlj0rIqZFxLSJEycWWZ2ZmY1RkXMKfxURb5N0O0BEPFn25HBEPCVpDukw1BqSlst7CxsCD5cNurDNNqus6L7k+rIqOb/Kqam+ijQKL0talvyLXtJEYPFIC+X5Xs4NwsrAjsBJwI2k/pMuAg4Arhhl7CPznczluL6sSs6vcmqqryKNwmnA94F1JX2RtEE/rsBy6wHn5QZlGeCSiLha0v8CF0k6gXS/wzmjC93MzDqtSN9HF0qaSzrZDPChiPh1geV+Bby1xfT7SOcXque+j8pxfY3JcDfJGc6vssZz30fAKsDAIaSVqwunw/xM2HLcV5RVyflVznh7RvMASZ8D9gQuI9289i1J34uIE6oObsy+9KW6IzAzG52atl9F9hT2Bd4aEX+CwQ7yfgmM/0bBzMxKKXKfwgJgpYb3K5JuQhv/dt89DWZmvaam7VeRPYWXgLslXUc6p/B+Uj9GpwFExCcrjG9s/IzmcqZOrTsC62fOr3LG8TOav5+HAXOqCcVqd+qpdUdg/cz51ROKXJJ6XjcCMTOz+hV9noItDT784TSYVcH51ROK3qdgS4MHH6w7Autnzq+e0HZPQdL5efyp7oVjZmZ1Gu7w0dslbQT8g6Q1Ja3VOHQrQDMz657hDh+dCfwI2ASYS7qbeUDk6WZm1kfaNgoRcRpwmqQzIuLgLsbUOe94R90R9BbXl1XJ+VVOTfVV5JLUgyW9BXh3nvTT3APq+Oe+j8pxfVmVnF/l1FRfRZ7R/EngQmCdPFwo6dCqAzMzs+4rcknqP5Ieyfk8gKSTgJuBr1cZWEcM9Bty2WX1xtErXF9WJedXOTXVV5FGQcArDe9f4dUnnccvH8Msx31FWZWcX+WM13MKwLeAWyQN9H/0IXrlEZpHHll3BGZmo1PT9qvIieavSpoDbEvaQ/hoRNxedWBmZtZ9hbq5iIhfkh6s01u22y6N58ypMwozs/Jq2n657yMbssMOdUdg/cz51RPcKNiQz3627gisnzm/esKw9ylIWlbS9d0KxszM6jVsoxARrwAvSHpNl+KxOu20UxrMquD86glFDh/9CbgzP6P5+YGJ4/rZzDY6L75YdwTWz5xfPaFIo3BNHszMrM8VekazpJWByRExvwsxmZlZTYp0iDcdmEd6tgKSpkq6surAzMys+4ocPjoe2BqYAxAR8yRtXGFMnbPrrnVH0FtcX1Yl51c5NdVXkUZhUUQ8Lb2qD7yoKJ7Oct9H5bi+rErOr3LGa99HwF2S/h5YVtKmwCeBm6oNy8zM6jDiOQXgUGBL4CXgu8AzwGFVBtUx22031H+Ijcz1ZVVyfpVTU30VufroBeDY/HCdiIhnqw+rQw48sO4IzMxGp6bt14iNgqStgHOB1fL7p4F/iIi5Fcc2dm4UbARTZvoWHBunatp+FTl8dA5wSERMiYgpwD+THrwzLEmTJN0o6deS7pb0qTx9LUnXSfpdHq85pr9gOI8/ngYzs15T0/arSKPwbET8bOBNRPwcKHIIaRFwRES8AdgG+GdJWwAzgRsiYlPghvy+GnvskQYzs15T0/ar7eEjSW/LL2+V9E3SSeYA9ibfszCciFgILMyvn5X0a2ADYDdguzzbebmso0YVvXXWXnvVHYH1M+dXTxjunMIpTe//peF1qfsUJE0B3grcAqybGwwiYqGkdcqUZRU65JC6I7B+5vzqCW0bhYh4XydWIGlV4DLgsIh4pukmuOGWmwHMAJg8eXInQrGRvPBCGq+ySr1xWH9yfvWEIlcfrQHsD0xpnL9I19mSlic1CBdGxOV58iOS1st7CesBj7ZaNiLOAs4CmDZtWm/cQd3rdt45jf1Ma6uC86snFLmj+VrgF8CdwOKiBSvtEpwD/Doivtrw0ZXAAcCsPL6icLRmZlapIo3CShFx+CjKfhfwEdIDeublaceQGoNLJB0E/AHYcxRlm5lZBYo0CudL+ifgalJXFwBExB+HWyhfutruBMIOhSM0M7OuKdIo/Bn4CnAsQ1cdBbBJVUGZmVk9ijQKhwOvj4jeuzXY3VyU4/qyKjm/yhmvfR8BdwMvVB1IJZyE5bi+rErOr3LGcaPwCjBP0o28+pzCiJek1m6g35AJE+qNo1f0QH25A7se1gP5Na7UVF9FGoUf5KH3DPQb4uuii3F9WZWcX+XUVF9FnqdwXjcCqcQRR9QdgZnZ6NS0/SpyR/P9tOjrKCLG/9VH06fXHYGZ2ejUtP0qcvhoWsPrlUg3m61VTTgdNn9+Gm++eb1xWGlL+7mD4f7+BbN26WIkVpuatl9FDh890TTpVEk/Bz5XTUgd9LGPpbGPYZpZr6lp+1Xk8NHbGt4uQ9pzWK2yiKw+Bx9cdwTWz5xfPaHI4aPG5yosAhYAflpGP9p777ojsH7m/OoJRQ4fdeS5CtYDHnggjSdNqjcO60/Or55Q5PDRisDuLPk8hS9UF5bV4iMfSWOfg7EqOL96QpHDR1cATwNzabij2czM+k+RRmHDiPhg5ZGYmVntlikwz02S3lR5JGZmVrsiewrbAgfmO5tfIj04JyLizZVGZmYttbuxzTe1WScUaRR2qjyKqrjvo3JcX1Yl51c547Xvo4j4fTcCqYT7PirH9WVVcn6VU1N9FTmn0Lvmzx/qP8RG5vqyKjm/yqmpvoocPupd7vuoHNeXVcn5Vc547fuop514Yt0RmJmNTk3br/5uFN75zrojMDMbnZq2X/19TuGmm9JgZtZratp+9feewjHHpLGPYZpZr6lp+9XfjYKVc9xxdUdg/cz51RPcKNiQHXesOwLrZ86vntDf5xSsnHnz0mBWBedXT/Cegg057LA09jkYq4Lzqyd4T8HMzAa5UTAzs0FuFMzMbJAbBTMzG9TfJ5rd91E5ri+rkvOrnH7r+0jSucCuwKMR8cY8bS3gYmAKsADYKyKerCoG931UkuvLquT8KqcP+z6aDXywadpM4IaI2BS4Ib+vjvs+Ksf1ZVVyfpXTb30fRcRPJU1pmrwbsF1+fR4wBziqqhjc91FJri+rkvOrnKWk76N1I2IhQEQslLROpWv75jcrLd7MrDI1bb/G7YlmSTOAGQCTJ08eXSGbb97BiKzTpsy8pu4Q+kq7+lwwa5cuR2IdUdP2q9uXpD4iaT2APH603YwRcVZETIuIaRMnThzd2q66Kg1mZr2mpu1Xt/cUrgQOAGbl8RWVru2UU9J4+vRKV2Nm1nE1bb+qvCT1u6STyhMkPQj8C6kxuETSQcAfgD2rWr+Nwqmn1h2B9TPnV0+o8uqjfdt8tENV67Qxmjq17gisnzm/eoK7ubAh11+fBrMqOL96wri9+shqcMIJaewnZFkVnF89wXsKZmY2yI2CmZkNcqNgZmaD3CiYmdmg/j7R7L6PynF9WZWcX+W476MKuO+jclxfViXnVzlLSd9H3eW+j8pxfVmVnF/lLCV9H3WX+z4qx/VlVXJ+ldNvfR+NC5deWncEZmajU9P2q78bhQkT6o7AzGx0atp+9XejMHt2Gh94YJ1RLPX8MB2zUahp+9XfJ5pnzx6qWDOzXlLT9qu/9xSsnPPPrzsC62fOr57gRsGGTJpUdwTWz5xfPcGNgg25+OI03nvveuOwjhrunM6CWbt0LxDnV09wo2BDzjgjjf2ltSo4v3pCf59oNjOzUtwomJnZIDcKZmY2yI2CmZkN6u8Tze77qBzXl1XJ+VWO+z6qgPs+Ksf1ZVVyfpVTU3319+Ejd3NRjuvLquT8KsfdXFTAHeKVU6C+xs2NUFapSv7P/j6WU1N99XejMGdO3RGYmY1OTduv/j58ZGZmpfR3o3DyyWkwM+s1NW2/+rtRuPrqNJiZ9Zqatl/9fU7Byrn22rojsH7m/OoJbhRsyCqr1B2B9TPnV0/o78NHVs7pp6fBrArOr57gPQUbcsklaXzIIaNafLhr22186ur/bIz5Zd3hPQUzMxtUS6Mg6YOS5ku6R9LMOmIwM7Mldb1RkLQs8A1gJ2ALYF9JW3Q7DjMzW1IdewpbA/dExH0R8WfgImC3GuIwM7MmdZxo3gB4oOH9g8BfNc8kaQYwI799TtL8Ua9RGu7TCcDjoy67OvXF5foq5B0DL07adbjZ+q6+dNIY1+z8Kmds9bVR2dXV0Si0+gtjiQkRZwFnVR6MdFtETKt6PWU5rnIcVzmOq5ylKa46Dh89CExqeL8h8HANcZiZWZM6GoX/ATaVtLGkFYB9gCtriMPMzJp0/fBRRCyS9AngP4FlgXMj4u5ux9Gg8kNUo+S4ynFc5TiucpaauBSxxOF8MzNbSvmOZjMzG+RGwczMBvVNoyDpXEmPSrqrafqhuUuNuyV9OU97v6S5ku7M4+3blHm8pIckzcvDzl2IbYqkFxvWeWabMteSdJ2k3+XxmhXGtF9DPPMkLZY0tUWZldWXpIsbyl0gaV7DZ0fnLlPmS/pAmzI3lnRLrq+L80UOlcXVzRwrGVdX8msUcXUtx9rENVXSL3KZt0naOk+XpNNyfv1K0tvalPn2/L++J88/7M0FHYhrvxzPryTdJOktbcqcLen+hvpaok6XEBF9MQDvAd4G3NUw7X3A9cCK+f06efxWYP38+o3AQ23KPB44ssuxTWmcb5gyvwzMzK9nAidVFVPTcm8C7ut2fTV9fgrwufx6C+AOYEVgY+BeYNkWy1wC7JNfnwkcXHFcXcuxknF1Jb/KxtXNHGuT+z8GdsqvdwbmNLz+Ien+qm2AW9qUeSvpfkbl+XeqOK53Amvm1zsNE9dsYI8ycfTNnkJE/BT4Y9Pkg4FZEfFSnufRPL49IgbujbgbWEnSiuMhthJ2A87Lr88DPtSlmPYFvlsu1HLaxAakX27AXg0x7AZcFBEvRcT9wD2krlSal9keuDRPKl1fZePqZo6VrK+ixpRfY4yr0hxrE1cAq+fXr2Ho3qndgG9H8gtgDUnrNS6Y368eETdH2hJ/m87VV8u4IuKmiHgyT/8F6X6vjuibRqGNzYB357UcIbUAAAYFSURBVMMG/yVpqxbz7A7cPrAhbOETeRft3NHsQo8yto0l3Z6nv7vN8utGxEKAPF6n4pgG7M3wX9iq6mvAu4FHIuJ3+X2rblM2aFpmbeCpiFg0zDydjqtRXTnWLq668mukuAbUkWOHAV+R9ABwMnB0nl4kvzbI04ebp9NxNTqItHfSzhdzfX2tyA+Tfm8UlgPWJO32fQa4pPFYn6QtgZOAj7VZ/gzgdcBUYCFpd7fq2BYCkyPircDhwHckrd6+mI4aqb7+CnghIu5qs3yV9TWg+VdkkW5TCnWtMkYtf93WnGOt4qozv4aLC6g1xw4GPh0Rk4BPA+cMhNRi3m7mV7u40oql95EahaPaLH808JfAVsBaw8w3qN8bhQeBy/Ou363AYlIHUkjaEPg+sH9E3Ntq4Yh4JCJeiYjFwH/QdFiiitjyYZAn8vrnko6Rb9Zi+UcGdmPzuOzhp8IxNXy+D8P8gqu4vpC0HPB3wMVNMY/UbcrjpN3+5YaZp9Nx1Z5jreKqOb/axtWgrhw7ALg8v/5eQ7lF8utBXn34ppP51S4uJL0ZOBvYbeB/2iwiFubv80vAtyhQX/3eKPyAdCwZSZsBKwCPS1oDuAY4OiL+u93CTccO/xZo9+ulk7FNVHrmBJI2ATYF7mux/JWkhCGPr6gqpvx+GWBPUlfnLVVcXwA7Ar+JiMZd9SuBfSStKGljUn3d2rhQPs57I7BHntSp+mob1zjJsVZx1ZlfbePK8dSZYw8D782vtwcGDmtdCeyvZBvg6YHDagPy+2clbZP3rPenc/XVMi5Jk0mNxUci4rftFm5o2EU6zzFyfZU5Kz2eB9Kvi4XAy6SW+yDSRu2CXBG/BLbP8x4HPA/MaxgGrv45G5iWX58P3An8ipQc63Uhtt1JJybvyNOnN5TTGNvawA05SW4A1qoqpjz/dsAvWpTTlfqKoSspPt5i/mNJv3jn03DVB3AtQ1cAbUJqLO4h/eJascq4upljJePqSn6N8v/YlRxrk/vbAnNzvdwCvD3PK9JDwe7N653WUM68htfTSN+be4F/J/cWUWFcZwNPNuTWbW3y/ic57rtI3+1VR4rD3VyYmdmgfj98ZGZmJbhRMDOzQW4UzMxskBsFMzMb5EbBzMwGuVGwvqHUA2ip69YlfVzS/iPMc6Ckf2/z2THDLCdJP+nEHcOSrq+o2xCzV3GjYEu1iDgzIr49hiLaNgqkXi3viIhnxlD+gPOBQzpQjtmw3ChYv1lW0n8oPQ/ix5JWBpD0Okk/Unq2wc8k/WWefrykI/PrrXLHYTdL+krTXsf6efnfaeg5E7OAlZX6qb+wRSz70XBnq6T9c/l3SDo/T5st6QxJN0q6T9J7c0dvv5Y0u6GsK0n9BZlVyo2C9ZtNgW9ExJbAU6Q7eCE94PzQiHg7cCRweotlv0W60/YdwCtNn00l9d75JmBvSZMiYibwYkRMjYj9WpT3LtLdqAMd4x1Lukv8LcCnGuZbk9SFwaeBq4CvAVsCb1J+KEqkbpJXlLR28aowK2+5kWcx6yn3R8TAU9nmAlMkrUp6KMn3Gjp9fVUXwrmvotUi4qY86TvArg2z3BART+d5/xfYiFd3qdzKWhHxbH69PXBpRDwOEBGN/eZfFREh6U5Sd9J35vXcTXoozsDf8yiwPtCy8zOzTnCjYP2m8ZkFrwArk/aIn4qI4R5FONLjE5vLLfLdWSRpmUg9eor23SkPlL24aT2Lm9azEvBigfWajZoPH1nfyyd675e0JwxeFfSWpnmeJPd0mSftU7D4lyUt3+az+aSO+CB1KrfXwOEfSWuV+RtyL5evBRaUWc6sLDcKtrTYDzhI0h2kXkJ3azHPQcBZkm4m/bJ/ukC5ZwG/anOi+RpSz59ExN3AF4H/yjF8tWT8byf1ILpoxDnNxsC9pJplklaNiOfy65mkbpk/NcJiw5W3Hun5vu/vQGz/BlwZETeMtSyz4ficgtmQXSQdTfpe/B44cCyFRcTCfHns6h24V+EuNwjWDd5TMDOzQT6nYGZmg9womJnZIDcKZmY2yI2CmZkNcqNgZmaD/j9zGdSgdu1emAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "seed(47)\n",
    "ymeans = [ np.mean ( townsfolk_sampler(10) ) for i in range ( 365 ) ]\n",
    "_ = plt.hist(ymeans, bins=20)\n",
    "_ = plt.xlabel('height (cm)')\n",
    "_ = plt.ylabel('number of people')\n",
    "_ = plt.title('Distribution of heights in sample size 10')\n",
    "_ = plt.axvline(172, color='r')\n",
    "_ = plt.axvline(172+5, color='r', linestyle='--')\n",
    "_ = plt.axvline(172-5, color='r', linestyle='--')\n",
    "_ = plt.axvline(172+10, color='r', linestyle='-.')\n",
    "_ = plt.axvline(172-10, color='r', linestyle='-.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above is the distribution of the means of samples of size 10 taken from our population. The Central Limit Theorem tells us the expected mean of this distribution will be equal to the population mean, and standard deviation will be $\\sigma / \\sqrt n$, which, in this case, should be approximately 1.58."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ Verify the above results from the CLT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "171.8660049358649"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_mean = np.sum ( ymeans )/365\n",
    "sample_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5756704135286475"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "std_ymeans = np.std ( ymeans )\n",
    "std_ymeans"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember, in this instance, we knew our population parameters, that the average height really is 172 cm and the standard deviation is 5 cm, and we see some of our daily estimates of the population mean were as low as around 168 and some as high as 176."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ Repeat the above year's worth of samples but for a sample size of 50 (perhaps you had a bigger budget for conducting surveys that year!) Would you expect your distribution of sample means to be wider (more variable) or narrower (more consistent)? Compare your resultant summary statistics to those predicted by the CLT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAdXElEQVR4nO3debwcVZn/8c+XsAQMEAIBCVuCMiioLAbEEQVZfoqIoAIGGbbBQcZhE3AIi4oOYtARHNz4IVsEJCwurKKAZIBBYRIIS0QEQlAgkoAEEsBAkmf+OOdC0+m+t+69Xd333vq+X69+3a7t1NN1u54+darqlCICMzOrjuU6HYCZmbWXE7+ZWcU48ZuZVYwTv5lZxTjxm5lVjBO/mVnFOPEPAJLOkfTlFpW1oaSFkobl4amSPteKsnN5v5J0UKvK68V6T5P0rKS/Npi2o6Qn+1juByU9XHDePq+nQNkt+w50gqRTJV3S4jJPknReK8u0xIm/ZJJmS3pF0gJJ8yXdKelwSa9v+4g4PCL+o2BZu3Q3T0T8OSJGRMSSFsS+zM4cEbtFxOT+lt3LODYAjgM2i4i3trLsiLg9IjZtRVmSLpJ0Wh/jKPQdqJKIOD0iWlZpAZD0YUm3SnpB0uwG08fm6S9L+mNP+9tg5cTfHntExKrARsAk4ATg/FavRNLyrS5zgNgIeC4i5nY6EBv0XgIuAL7UZPplwL3AmsDJwFWSRrcptvaJCL9KfAGzgV3qxm0LLAXelYcvAk7L79cCrgPmA38Dbif9QF+cl3kFWAj8OzAWCOBQ4M/AbTXjls/lTQW+CdwNvABcDYzK03YEnmwUL/BR4FXgtby++2rK+1x+vxxwCvAEMBf4CbB6ntYVx0E5tmeBk7vZTqvn5efl8k7J5e+SP/PSHMdFDZbdEXiSdFQwF5gDHFIzfSXgP3MczwDnACs32gbA1qQdfwFwJXB5zf+m6XqAw/K2ejXHeW0efwLwVC7vYWDnJp//oiLrabLswcCsvI7Hgf3z+LcBvwWey9v/UmBk3f/6S8D9pIR4PrAO8Ktc1s3AGnX/z8OAp3NMx9WUdSpwSc3wdsCdpO/xfcCO3cTfcBvVlgl8P2/Xrtdi4NQ8bQzwM9J353HgqAL75S7A7Lpx/wAsAlatGXc7cHin80irX67xd0BE3E3asT/YYPJxedpo0k54UlokDiAlrj0iNeV8q2aZHYB3Ah9pssoDgX8m7SCLgbMLxHgjcDpweV7fFg1mOzi/PgxsDIwg7aC1tgc2BXYGviLpnU1W+T1S8t84f54DScnuZmA34Okcx8FNln9rXn490g/hDyStkaedQdqptwTenuf5Sn0BklYEfkFKwqNItb9PFllPRJxLSqzfynHuIWlT4Ahgm0hHfB8hJdsiuvs8tTG/hfT/3C2v4x+BGV2TST/6Y0jfjw1IybTWp4FdSdtnD1LSP4lUAVkOOKpu/g8DmwD/D5jYqClE0nrA9cBppO14PPCzRjXnotsoIo7I23UE6Tv1PHB1bjK9lvTjsh7pe3aMpGb7Qnc2B2ZFxIKacffl8UOKE3/nPE3aKeq9BqwLbBQRr0Vqg+6pQ6VTI+KliHilyfSLI+LBiHgJ+DKwb9fJ337aHzgzImZFxELgRGBCXZPT1yLilYi4j7QTLfMDkmP5DHBiRCyIiNnAd4ADehHLa8DX8za7gVQr3FSSgH8BvhgRf8s79enAhAZlbAcsD5ydy/k56Uipx/U0iWkJ6WhjM0krRMTsiHisP5+nybxLgXdJWjki5kTETICIeDQiboqIRRExDziT9KNa63sR8UxEPEWq3d4VEfdGxCLSj+BWdfN/LX/XHgAuBPZrEM8/ATdExA0RsTQibgKmAR9rMG+vtlH+8fglcGRE3AtsA4yOiK9HxKsRMQv4MY3/vz0ZQToqrvUCsGofyhrQnPg7Zz1SU069bwOPAr+RNEvSxAJl/aUX058AViDV6PprTC6vtuzlSUcqXWqvwnmZtHPVWwtYsUFZ6/UiluciYnGDdY0GVgGm55Pr84Eb8/h6Y4Cn6n5o67dts/UsIyIeBY4h1bLnSpoiaUw/P0/9Ol4i/WgeDsyRdL2kdwBIWjuv8ylJLwKXsOz//Zma9680GK5fZ/13qdHn2QjYp2t7522+PalCUx9/4W0kaQXgKuCnETGlZl1j6tZ1Em/+Dha1EFitbtxqpCaoIcWJvwMkbUNKanfUT8s13uMiYmPSofexknbumtykyJ6OCDaoeb8hqTb5LKldd5WauIbx5oTYU7lPk3a82rIX8+bkUcSzOab6sp7qZTnNyn4F2DwiRubX6rnJoN4cYL18lNBlgwbzNbPM9oqIn0bE9qTPFqRmp5aKiF9HxK6kxPpHUo0XUjNPAO+JiNVINXE1LqWw+u/S0w3m+QvpKHNkzestETGpSfxFt9H3SEn4lLp1PV63rlUjotHRRU9mAhtLqq3hb5HHDylO/G0kaTVJHwemkE5aPdBgno9LentOPi+SDoW7Ls18htQG3lv/JGkzSasAXweuinS555+A4ZJ2z7WpU0iH3V2eAcbWXnpa5zLgi5LGSRrBG+cEFjeZv6EcyxXANyStKmkj4FhSDbVfImIpKRGeJWltSG3QTdqAf0fa1kdIWl7SnqQT8UW96f8jaVNJO0laCfg76Qeo35fZ1pK0jqRP5Lb+RaRaa9c6Vs3D83O7e7MrWXrjy5JWkbQ5cAjp5He9S4A9JH1E0jBJw/M9EOs3iL/QNpL0eVIz1Wfz/7TL3cCLkk6QtHJe37ty5WoZkpaTNJx01Ksc24oAEfEn0vmRr+bxnwTeQzpxPKQ48bfHtZIWkGonJ5PaWg9pMu8mpKspFpIS0Q8jYmqe9k3glHxIe3wv1n8x6YTlX4Hh5BN2EfEC8AXgPFLt+iXSieUuV+a/z0m6p0G5F+SybyNdTfF34MhexFXryLz+WaQjoZ/m8lvhBFLz2e9zk8fNNGgvj4hXgU+RTqbOJ9WQryMl1CLOJ7VVz5f0S9KP6CTSUcdfgbVJzRCttBzpgoCnSU2HO5D+pwBfI12l9ALpZOvPW7C+/yZty1uA/4yI39TPEBF/AfYkfdZ5pO/9l2icb4puo/1IP6pPK92guFDSSbnSsAfpxP3juZzzSCfGG/kQ6cflBtIRyytA7WeYAIwnnTyeBOydz48MKer5vKFZdUm6CzgnIi7sdCydJGksKbGu0NsjOht4XOM3qyFpB0lvzU09B5EO9W/sdFxmrTRU7/Q066tNSecbRgCPkQ7153Q2JLPWclOPmVnFuKnHzKxiBkVTz1prrRVjx47tdBhmPXs49/C8aUs6/DTrl+nTpz8bEcvcrDgoEv/YsWOZNm1ap8Mw69mOO6a/U6d2MgozACQ90Wi8m3rMzCrGid/MrGKc+M3MKsaJ38ysYpz4zcwqxonfzKxinPjNzCrGid/MrGKc+M3MKmZQ3LlrNlCNnXj9m4anzHoOgAl14xuZPWn3UmIy64lr/GZmFePEb2ZWMU78ZmYV48RvZlYxTvxmZhXjxG9mVjFO/GZmFePEb2ZWMU78ZmYV48RvZlYxTvxmZhXjxG9mVjFO/GZmFePEb2ZWMaUnfknDJN0r6bo8PE7SXZIekXS5pBXLjsHMzN7Qjhr/0cBDNcNnAGdFxCbA88ChbYjBzMyyUhO/pPWB3YHz8rCAnYCr8iyTgb3KjMHMzN6s7Br/d4F/B5bm4TWB+RGxOA8/CaxXcgxmZlajtMQv6ePA3IiYXju6wazRZPnDJE2TNG3evHmlxGhmVkVl1vg/AHxC0mxgCqmJ57vASEldz/pdH3i60cIRcW5EjI+I8aNHjy4xTDOzaikt8UfEiRGxfkSMBSYAv42I/YFbgb3zbAcBV5cVg5mZLasT1/GfABwr6VFSm//5HYjBzKyylu95lv6LiKnA1Px+FrBtO9ZrZmbL8p27ZmYV48RvZlYxTvxmZhXjxG9mVjFO/GZmFePEb2ZWMU78ZmYV48RvZlYxTvxmZhXjxG9mVjFO/GZmFePEb2ZWMU78ZmYV05beOc0GsrETr+90CGZt5Rq/mVnFOPGbmVWME7+ZWcU48ZuZVYwTv5lZxTjxm5lVjBO/mVnFOPGbmVWME7+ZWcX4zl0bEnz3rVlxrvGbmVWME7+ZWcU48ZuZVYzb+M0qqD/nRGZP2r2FkVgnuMZvZlYxTvxmZhXjxG9mVjFO/GZmFdPjyV1J6wCnA2MiYjdJmwHvj4jzS4/OzJryTWvWV0Vq/BcBvwbG5OE/AceUFZCZmZWrSOJfKyKuAJYCRMRiYEmpUZmZWWmKJP6XJK0JBICk7YAXSo3KzMxKU+QGrmOBa4C3SfofYDSwd6lRmZlZaXpM/BFxj6QdgE0BAQ9HxGulR2ZmZqVomvglfarJpH+QRET8vLuCJQ0HbgNWyuu5KiK+KmkcMAUYBdwDHBARr/YpejMz67Xuavx7dDMtgG4TP7AI2CkiFkpaAbhD0q9ITUdnRcQUSecAhwI/6k3QZmbWd00Tf0Qc0p+CIyKAhXlwhfwKYCfgs3n8ZOBUnPjNzNqmx6t6JK0p6WxJ90iaLum/8lU+PZI0TNIMYC5wE/AYMD9fEgrwJLBeX4M3M7PeK3I55xRgHvBp0tU884DLixQeEUsiYktgfWBb4J2NZmu0rKTDJE2TNG3evHlFVmdmZgUUSfyjIuI/IuLx/DoNGNmblUTEfGAqsB0wUlJXE9P6wNNNljk3IsZHxPjRo0f3ZnVmZtaNIon/VkkTJC2XX/sCPXYSImm0pJH5/crALsBDwK28cR/AQcDVfQvdzMz6osgNXJ8nXYlzSR5ejnQ377Gkc7irNVluXWCypGF5mSsi4jpJfwCmSDoNuBdwZ29mZm1U5AauVftScETcD2zVYPwsUnu/2Zu4t0mz9ij0zF1JnwA+lAenRsR15YVkZmZlKnI55yTgaOAP+XV0HmdmZoNQkRr/x4AtI2IpgKTJpLb5iWUGZmZm5Sj66MXayzdXLyMQMzNrjyI1/m8C90q6ldQ754eAE0uNygYln5w1GxyKXNVzmaSpwDakxH9CRPy17MDMzKwcRU7uCtiZ1M5/NbCiJF+OaWY2SBVp4/8h8H5gvzy8APhBaRGZmVmpirTxvy8itpZ0L0BEPC9pxZLjMjOzkhSp8b+Wu13oetj6aGBpqVGZmVlpiiT+s4FfAOtI+gZwB3B6qVGZmVlpilzVc6mk6aQTvAB7RcRD5YZlZmZlKdRXD7AK0NXcs3J54ZiZWdl6TPySvgLsA/yMdB3/hZKuzA9kMbOK6c+NerMn7d7CSKyvitT49wO2ioi/w+udtt0DOPGbmQ1CRU7uzgaG1wyvRHpoupmZDUJFavyLgJmSbiK18e8K3CHpbICIOKrE+MzMrMWKJP5f5FeXqeWEYmZm7VDkcs7J7QjEzMzao2h//GZmNkQ48ZuZVUzTxC/p4vz36PaFY2ZmZeuuxv9eSRsB/yxpDUmjal/tCtDMzFqru5O75wA3AhsD00l37XaJPN7M+siPqrROaVrjj4izI+KdwAURsXFEjKt5OembmQ1SRS7n/FdJWwAfzKNui4j7yw3LzMzKUuSZu0cBlwJr59elko4sOzAzMytHkTt3P0d6/OJLAJLOAH4HfK/MwMzMrBxFruMXsKRmeAlvPtFrZmaDSJEa/4XAXZK6+uvZCzi/vJDMzKxMRU7unilpKrA9qaZ/SETcW3ZgZmZWjkKPXoyIe0gPXzEzs0HOffWYmVWME7+ZWcV0m/glDZN0c7uCMTOz8nWb+CNiCfCypNXbFI+ZmZWsyMndvwMP5GfuvtQ10s/aNTMbnIok/uvzy8zMhoBCz9yVtDKwYUQ83IaYzMysREU6adsDmEHqmx9JW0q6psByG0i6VdJDkmZ2PckrP8jlJkmP5L9r9PdDmJlZcUUu5zwV2BaYDxARM4BxBZZbDByX+/TfDvg3SZsBE4FbImIT4JY8bGZmbVIk8S+OiBfqxkVPC0XEnHzHLxGxAHgIWA/YE5icZ5tM6vvHzMzapMjJ3QclfRYYJmkT4Cjgzt6sRNJYYCvgLmCdiJgD6cdB0tpNljkMOAxgww037M3qrB/8OECzoa9Ijf9IYHNgEXAZ8CJwTNEVSBoB/Aw4JiJeLLpcRJwbEeMjYvzo0aOLLmZmZj0oclXPy8DJ+QEskZttCpG0AinpXxoRP8+jn5G0bq7trwvM7UvgZmbWN0Wu6tlG0gPA/aQbue6T9N4Cy4nUb/9DEXFmzaRrgIPy+4OAq3sftpmZ9VWRNv7zgS9ExO0AkrYnPZzlPT0s9wHgANKPxYw87iRgEnCFpEOBPwP79CVwMzPrmyKJf0FX0geIiDsk9djcExF30PwRjTsXjM/MzFqsaeKXtHV+e7ek/086sRvAZ4Cp5YdmZmZl6K7G/5264a/WvO/xOn4zMxuYmib+iPhwOwMxM7P26LGNX9JI4EBgbO387pbZzGxwKnJy9wbg98ADwNJywzEzs7IVSfzDI+LY0iMxM7O2KNJlw8WS/kXSurlL5VGSRpUemZmZlaJIjf9V4NvAybxxNU8AG5cVlJmZladI4j8WeHtEPFt2MGZmVr4iTT0zgZfLDsTMzNqjSI1/CTBD0q2krpkBX85pZjZYFUn8v8wvMzMbAor0xz+5p3nMzGzwKHLn7uM06JsnInxVj5nZIFSkqWd8zfvhpP7zfR2/mdkg1eNVPRHxXM3rqYj4LrBTG2IzM7MSFGnq2bpmcDnSEcCqpUVkZmalKtLUU9sv/2JgNrBvKdHY68ZOvL7Py86etHsLIzFrnf58r8Hf7VYpclWP++U3MxtCijT1rAR8mmX74/96eWGZmVlZijT1XA28AEyn5s5dMzMbnIok/vUj4qOlR2JmZm1RpJO2OyW9u/RIzMysLYrU+LcHDs538C4CBEREvKfUyMzMrBRFEv9upUdhZmZtU+RyzifaEYiZmbVHkTZ+MzMbQpz4zcwqxonfzKxinPjNzCrGid/MrGKc+M3MKsaJ38ysYpz4zcwqxonfzKxinPjNzCrGid/MrGKc+M3MKqa0xC/pAklzJT1YM26UpJskPZL/rlHW+s3MrLEya/wXAfVP7poI3BIRmwC35GEzM2uj0hJ/RNwG/K1u9J7A5Px+MrBXWes3M7PG2t3Gv05EzAHIf9duNqOkwyRNkzRt3rx5bQvQzGyoG7AndyPi3IgYHxHjR48e3elwzMyGjHYn/mckrQuQ/85t8/rNzCqv3Yn/GuCg/P4g4Oo2r9/MrPKKPGy9TyRdBuwIrCXpSeCrwCTgCkmHAn8G9ilr/VU2duL1nQ7BzAaw0hJ/ROzXZNLOZa3TzMx6NmBP7pqZWTmc+M3MKsaJ38ysYpz4zcwqxonfzKxinPjNzCrGid/MrGKc+M3MKsaJ38ysYpz4zcwqxonfzKxinPjNzCrGid/MrGKc+M3MKsaJ38ysYpz4zcwqxonfzKxiSnsCl5lZq3XqsaKzJ+3ekfWWxTV+M7OKceI3M6sYJ34zs4px4jczqxif3O1Gf04kDbWTQWY2dLjGb2ZWMU78ZmYV48RvZlYxTvxmZhXjk7sl6dQdhmZmPXGN38ysYpz4zcwqxonfzKxi3MZvZlaigXgjqGv8ZmYV48RvZlYxTvxmZhXjxG9mVjFD/uSub6QyM3sz1/jNzCqmI4lf0kclPSzpUUkTOxGDmVlVtT3xSxoG/ADYDdgM2E/SZu2Ow8ysqjpR498WeDQiZkXEq8AUYM8OxGFmVkmKiPauUNob+GhEfC4PHwC8LyKOqJvvMOCwPLgp8HAPRa8FPNvicFtpoMcHjrFVHGNrOMb+2ygiRteP7MRVPWowbplfn4g4Fzi3cKHStIgY35/AyjTQ4wPH2CqOsTUcY3k60dTzJLBBzfD6wNMdiMPMrJI6kfj/F9hE0jhJKwITgGs6EIeZWSW1vaknIhZLOgL4NTAMuCAiZrag6MLNQh0y0OMDx9gqjrE1HGNJ2n5y18zMOst37pqZVYwTv5lZxQy4xC/pAklzJT1YM+5ySTPya7akGTXTTsxdPzws6SNNyhwn6S5Jj+SyVhyAMV4k6fGaMrZsV4yS1pR0q6SFkr7fTZmjJN2Ut+NNktYYgDGeKumpmjI+1sYYd5U0XdID+e9OTcrs5HYsGmMnt+O2NePvk/TJJmW2bL8uKb6W7tMtFRED6gV8CNgaeLDJ9O8AX8nvNwPuA1YCxgGPAcMaLHMFMCG/Pwf41wEY40XA3h3ajm8BtgcOB77fTZnfAibm9xOBMwZgjKcCx3doO24FjMnv3wU8NQC3Y9EYO7kdVwGWz+/XBeZ2Ddct07L9uqT4WrpPt/I14Gr8EXEb8LdG0yQJ2Be4LI/aE5gSEYsi4nHgUVKXEPXL7ARclUdNBvYaSDGWoTcxRsRLEXEH8Pceit2TtP2gzduxFzG2VC9jvDciuu5JmQkMl7RSg0U7uR2LxthSvYzx5YhYnCcPp8ENnq3er1sd30A34BJ/Dz4IPBMRj+Th9YC/1Ex/Mo+rtSYwv+Yf1WieTsfY5RuS7pd0Vsk7Y32MRa0TEXMA8t+1Wx7ZG/oaI8AReTte0N9mlB50F+OngXsjYlGDaQNlO3YXI3RwO0p6n6SZwAPA4TX7b5d27td9ia9Lu/bpXhlsiX8/3qhJQ7HuHwp1EdFCfYkR4ETgHcA2wCjghNaH9rr6GAeivsb4I+BtwJbAHNIhelkaxihpc+AM4PMlrruovsbY0e0YEXdFxOak/eFEScPrlmnnft2X+KC9+3SvDJrEL2l54FPA5TWji3T/8CwwMi/fbJ5Ox0hEzIlkEXAhJTUHNYmxqGckrZvL6WrbbLn+xBgRz0TEkohYCvyYNm9HSesDvwAOjIjHmize0e1YJMZOb8eaOB4CXiKdj6jVlv26H/G1bZ/ui0GT+IFdgD9GxJM1464BJkhaSdI4YBPg7tqFIp1luRXYO486CLh6IMUIryeArvbEvYAH6+cpMcairiFtP2j/diykaztmn6SN21HSSOB64MSI+J9ulu3YdiwaY4e347iuhC5pI1LvvLNrF2rjft2n+PK0du3Tvdfps8v1L9Ih1RzgNVJt+dB44wz54Q3mP5l0pczDwG4142/gjasXNiYl20eBK4GVBmCMvyW1Fz4IXAKMaHOMs0kntxbm+TfL488Dxuf3awK3AI/kv6MGYIwX5+14PynBrtuuGIFTSLW/GTWvtQfSduxFjJ3cjgeQTjzPAO4B9mqyz7Rsvy4pvpbu0618ucsGM7OKGUxNPWZm1gJO/GZmFePEb2ZWMU78ZmYV48RvZlYxTvw2ZEgaW9u7YsFlDpd0YA/zHKwmPYJKOqmb5STpt5JW601MTcq6ueRuE6xCnPit0iLinIj4ST+KaJr4gY8B90XEi/0ov8vFwBdaUI6ZE78NOcMk/VjSTEm/kbQygKS3SbpRqQ/62yW9I48/VdLx+f02uUOt30n6dt3Rw5i8/COSvpXnnwSsnPtav7RBLPtTczeppANz+fdJujiPu0jSj5SeNTBL0g65U7SHJF1UU9Y1pD5jzPrNid+Gmk2AH0TqQGs+qQdKSA/FPjIi3gscD/ywwbIXku7SfD+wpG7alsBngHcDn5G0QURMBF6JiC0jYv8G5X0AmA6vd4p2MrBTRGwBHF0z3xqkLoa/CFwLnAVsDrxb+eEdEfE8sJKkNYtvCrPGlu95FrNB5fGI6Hr62XRgrKQRwD8CV6ZuU4D0YJzX5T5sVo2IO/OonwIfr5nlloh4Ic/7B2Aj3tzddiOjImJBfr8TcFVEPAsQEbV9v18bESHpAVL3vw/k9cwExpK6BoDUmdsY4Lke1mvWLSd+G2pq+5ZfAqxMOrKdHxHdPfquUTe/3ZVbZN9ZLGm5SD1ciubdBneVvbRuPUvr1jMceKXAes265aYeG/LyydXHJe0Dr19ts0XdPM8DCyRtl0dNKFj8a5JWaDLtYVJHYpA6Y9u3q6lG0qjefIbcw+NbadALpFlvOfFbVewPHCrpPlLPins2mOdQ4FxJvyPV0F8oUO65wP1NTu5eD+wIEBEzgW8A/51jOLOX8b8X+H00f9KTWWHundMskzQiIhbm9xNJXREf3cNi3ZW3LvCTiNi1BbH9F3BNRNzS37LM3MZv9obdJZ1I2i+eAA7uT2ERMSdfWrpaC67lf9BJ31rFNX4zs4pxG7+ZWcU48ZuZVYwTv5lZxTjxm5lVjBO/mVnF/B/Q+tCrrl+HIwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "seed(47)\n",
    "# calculate daily means from the larger sample size here\n",
    "ymeans = [ np.mean ( townsfolk_sampler(50) ) for i in range ( 365 ) ]\n",
    "_ = plt.hist(ymeans, bins=20)\n",
    "_ = plt.xlabel('height (cm)')\n",
    "_ = plt.ylabel('number of people')\n",
    "_ = plt.title('Distribution of heights in sample size 10')\n",
    "_ = plt.axvline(172, color='r')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "171.94366080916114"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_mean = np.sum ( ymeans )/365\n",
    "sample_mean "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6736107539771146"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "std_ymeans = np.std ( ymeans )\n",
    "std_ymeans"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we've seen so far, then, is that we can estimate population parameters from a sample from the population, and that samples have their own distributions. Furthermore, the larger the sample size, the narrower are those sampling distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### III Normally testing times!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the above is well and good. We've been sampling from a population we know is normally distributed, we've come to understand when to use $n$ and when to use $n-1$ in the denominator to calculate the spread of a distribution, and we've  seen the Central Limit Theorem in action for a sampling distribution. All seems very well behaved in Frequentist land. But, well, why should we really care?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember, we rarely (if ever) actually know our population parameters but you still have to estimate them somehow. If we want to make inferences such as \"is this observation unusual?\" or \"has my population mean changed?\" then you need to have some idea of what the underlying distribution is so you can calculate relevant probabilities. In frequentist inference, you use the formulas above to deduce these population parameters. Take a moment in the next part of this assignment to refresh your understanding of how these probabilities work."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall some basic properties of the standard Normal distribution, such as about 68% of observations being within plus or minus 1 standard deviation of the mean.\n",
    "\n",
    "__Q:__ Using this fact, calculate the probability of observing the value 1 or less in a single observation from the standard normal distribution. Hint: you may find it helpful to sketch the standard normal distribution (the familiar bell shape) and mark the number of standard deviations from the mean on the x-axis and shade the regions of the curve that contain certain percentages of the population."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3hUZfbA8e9JBxJ6qCEkgYB06QKCqCigCLuurIi6WNG167rqrh111766iquoiB2xgopdQYoovbcQAoReQ0lC2vn9MTf+xpgyCZnczOR8nmeezNx73ztnksmced/7FlFVjDHGmNKEuB2AMcaY6s+ShTHGmDJZsjDGGFMmSxbGGGPKZMnCGGNMmSxZGGOMKZMlC1NuIvKAiLxVwbKXicjcUvZ/ISLjijtWRI6KSFJFnrciRORFEbm3ks4V78Qf6jyeJSJXVca5nfP9+nsLNiJysYh87XYcNZ0lixpCRNJEJMv5wNotIq+JSLTbcRWlqsNV9fUS9kWraiqAiEwRkYcr+jxev48jInJIROaLyLUi8uv/hKpeq6oP+XiuIaUdo6pbnfjzKxqz1/P9LlmX9ns7weeaIiI5zvum8HZhZT+P1/MliIiKSFjhNlV9W1XP9tdzGt9YsqhZzlPVaKAH0Bu4p+gB4lFT3hfnqWoM0Bp4FLgTeLWyn8T7gy9APe4kusLbe24HZKpeTflQMF5UdTvwBdAZfm0SeURE5gGZQJKItBCRGSJyQERSROTqIqeJEpH3nG/mS0SkW+EOEblLRDY5+9aIyB+LlBUReU5EMkRknYic6bWjxOYZ5xtnWxEZD1wM3OF80/1URP4uIh8WOf45EXnGh99HhqrOAC4ExolI4e/l19qLiDQWkc+cWsgBEZkjIiEi8iYQD3zqxHKH17fjK0VkK/B9cd+YgTYi8ovze5guIg2d5xosIulFXkuaiAwRkWHAP4ELnedbXvT35sR1j4hsEZE9IvKGiNRz9hXGMU5EtorIPhG5u6zfUXEK/x5ej71/X4NFJF1E/ubEsFNELvc6tpaIPOXEmCEic0WkFvCjc8gh5/X1k983R/YXkYVOuYUi0t9r3ywReUhE5jnvv69FpHFFXp/5LUsWNZCItALOAZZ6bb4UGA/EAFuAd4F0oAVwAfAv7w91YBTwPtAQeAf4RETCnX2bgIFAPeBB4C0Rae5Vti+QCjQG7gc+Kvyg9IWqTgLe5v+/8Z4HvAUME5H6zmsMw/Ph/2Y5zvuL85oHFrP7b86+WKApng9sVdVLga04tTZVfdyrzGlAB2BoCU/5F+AKPL/jPOC/PsT4JfAv4D3n+boVc9hlzu10IAmIBp4vcsypQHvgTOA+EelQ1nNXQDM874GWwJXARBFp4Ox7EugJ9MfzHroDKAAGOfvrO6/vJ+8TOu+Tz/H8rhoBTwOfi0gjr8PGApcDTYAI4PbKf2k1jyWLmuUTETkEzAVm4/nQKTRFVVerah6ef/JTgTtVNVtVlwGv4EkohRar6geqmovnHzYKOAVAVd9X1R2qWuA0WWwE+niV3QM8o6q5zv71wLkn8sJUdSeeb6WjnU3DgH2quricp9qB58OrqFygOdDaiXuOlj2x2gOqekxVs0rY/6aqrlLVY8C9wJ/FuQB+gi4GnlbVVFU9CvwDGFOkVvOgqmap6nJgOVBc0il0u1OjOiQi+8oRRy4wwfl9zQSOAu3F08x5BXCzqm5X1XxVna+qx30457nARlV9U1XzVPVdYB1wntcxr6nqBuf3Pg04uRwxmxJYsqhZ/qCq9VW1tapeV+RDbJvX/RbAAVU94rVtC55viL87XlUL+P9aCCLyFxFZVvgBg6e5y7spYHuRD9othWVP0OvAJc79SyhHrcJLS+BAMdufAFKAr0UkVUTu8uFc28qxfwsQzm9/TxXVwjmf97nD8NSICu3yup+Jp/ZRkied9019VS1PfPudLx9Fn6cxni8Xm8pxrkJFXxv8/r1ZntdmfGTJwhTy/vDeATQUkRivbfHAdq/HrQrvON8U44AdItIaeBm4AWikqvWBVYB4lW0pIt6P453nrGi8hT4BujrXHEbgaarymYj0xvOh87uuvap6RFX/pqpJeL7F3ubVLFdSDaOsmkcrr/vxeL6J7wOOAbW94grF0/zl63l34Llo733uPGB3GeXKKxOvOPHUSH2xD8gG2hSzr7yvDX7/3jR+YMnC/I6qbgPmA/8WkSgR6Yqnzdn7w7eniJzvNG3cAhwHFgB18PzD7wVwLmp2LvIUTYCbRCRcREbjadefWc4wd+Npj/eOOxv4AM81lF9UdasvJxKRuiIyApgKvKWqK4s5ZoRzcV2Aw0C+cys2Fh9dIiIdRaQ2MAH4wOlauwFPB4JznetA9wCRXuV2AwlScq+1d4FbRSRRPN2jC69x5JVwfEUtA8aKSKhz4f00Xwo5NdHJwNPi6UgR6lzIjsTzvimg5N/nTKCdiIwVkTDxdOPtCHx2wq/GlMqShSnJRUACnm9yHwP3q+o3Xvun47mAfBDPtYzznbbpNcBTwE94PtS6APOKnPtnIBnPN8xHgAtUdX8543sV6Og0dX3itf115zl9aYL6VESO4GkOuhvPtZfLSzg2GfgWT7v7T8ALqjrL2fdv4B4nlvJcTH0TmIKn2SQKuAk8vbOA6/BcJ9qOp6bh3TvqfefnfhFZUsx5Jzvn/hHYjOdb/I3liMtXN+OpZR3Cc53kk9IP/43bgZXAQjzNfo8BIaqaiec9Mc/5fZ7iXch5n4zA0+FgP54L4yNUtTzXUkwFiC1+ZIKJiMTjueDZTFUPux2PMcHCahYmaDjNMrcBUy1RGFO5An1kqTEAiEgdPM1eW/B0mzXGVCJrhjLGGFMma4YyxhhTpqBphmrcuLEmJCS4HYapgPX71wPQvlF7lyMxFWF/v8C2ePHifaoaW9ZxQZMsEhISWLRokdthmAoYPGUwALMum+VqHKZi7O8X2ESk6Ij4YlkzlDHGmDJZsjDGGFMmSxbGGGPKZMnCGGNMmSxZGGOMKZNfk4WIDBOR9eJZlrPE+f9F5AJnicZeXtv+4ZRbLyIlrTRmjDGmCvit66wzB/9E4Cw8M2YuFJEZzqyk3sfF4Jlt82evbR2BMUAnPIudfCsi7Zzpm40xxlQxf9Ys+gApztKOOXjWChhVzHEPAY/jmUa50Cg8k8EdV9XNeFYo61NMWWOqraycfDbvO/brbev+TPILbHodE5j8OSivJb9dNjId6Ot9gIh0B1qp6mdF1gFoiWchHe+y3ssmFpYfD4wHiI+Pr6SwjamYvPwCvlu3h5827WfJ1oOs2XGYvCLJIToyjO7x9ekR34AhHZrSJa6eS9EaUz7+TBZSzLZf/3Oc6aT/A1xW3rK/blCdBEwC6NWrl31lM644kp3Lewu38dq8NLYfyqJWeCjdWtXjmtOSaBMbTYizgmx2bj4rt2ewZOshnvt+I89+t5E+iQ256tREhnRoSkhIcW97Y6oHfyaLdH67xnAcv11nOQbPcpuznOWYmwEzRGSkD2WNcV1BgTJlfhr/+WYDR47n0SexIfef15EzTmpCWGjxLbxjnJ8ZWbm8v8iTYMa/uZg2sXX41x+70DepUdW9AGPKwZ/JYiGQLCKJeJaGHAOMLdzpLB3ZuPCxiMwCblfVRSKSBbwjIk/jucCdDPzix1iNKZfN+45xxwfLWZh2kMHtY7ntrHZ0javvc/l6tcK5amASl/VP4ItVu3jiq/VcOGkBl/VP4I5h7akdETTTtpkg4bd3pKrmicgNwFdAKDBZVVeLyARgkarOKKXsahGZBqwB8oDrrSeUqS7eXLCFhz9bQ2RYCE+N7sb5PVri1I7LLSw0hPO6teDMDk14/Mv1TJmfxvfr9vDcRd3p1sr35GOMv/n164uqzgRmFtl2XwnHDi7y+BE8C7cbUy3kFyiPfL6WyfM2M7h9LI/9qStN60ZVyrlrR4TxwMhODO/cjNumLefCST/x3EU9OKtj00o5vzEnykZwG+OD7Nx8rn97CZPnbeay/gm8Oq53pSUKb32TGvHJ9QNo3zSGa95cxBs/pVX6cxhTEZYsjCnD4excxr68gK/W7OLeER15YGQnQv3Ycyk2JpJ3x5/CGSc14b7pq3nsy3XY8sfGbZYsjClFdm4+V01ZxIr0DCaO7cGVpyZWyfPWjgjjpUt7MbZvPP+btYnnv0+pkuc1piTW5cKYEuTmF3Dd20tYuOUAz47pzjldmlfp84eGCA+P6kx2bj5PfbOBerXD+Uu/hCqNwZhCliyMKUZBgXL7+8v5ft0eHvljZ0Z2a+FKHCEhwuN/6srhrDzum76aulHh/KH77yYzMMbvrBnKmGI89uU6pi/bwR3D2nNx39auxhIWGsLzY7vTL6kRf3t/OfM37XM1HlMzWbIwpojPVuzgpR9TueSUeP56Whu3wwEgKjyUSX/pSWLjOtzwzlJ2HMpyOyRTw1iyMMbL+l1HuOODFfRs3YD7RnSq8GA7f4iJCuelS3uSk1fAtW8tJjvXxqmaqmPJwhhHRlYu17y5iDqRYbxwcQ8iwqrfv0eb2Gie/nM3VqRncN/0Vdal1lSZ6vffYIwLVJW/TVtG+sEsXri4h18G3FWWszs148Yz2jJtUTpTF24ru4AxlcCShTHAO79s5du1e/jHOR3ondDQ7XDKdMuQdgxMbsyET9ewed8xt8MxNYAlC1Pjbd53jIc/W8vA5MZc3j/B7XB8EhoiPHFBNyLCQrj1vWXk5Re4HZIJcpYsTI2Wl1/Are8tIyIshCcu6BZQCxA1qxfFw3/ozLJth3hh1ia3wzFBzpKFqdFemLWJZdsO8fAfOtOsXvW9TlGS87q1YNTJLXj2u40s33bI7XBMELNkYWqsVdszePa7jYw6uQXnuTRCuzJMGNmZJjGR3DptmXWnNX5jycLUSPkFyj8+WkmD2hFMGNnZ7XBOSL3a4Tz6p66k7j3G/6w5yviJJQtTI73xUxort2dw33kdqVc73O1wTthp7WIZ2a0F/5u1iZQ9R90OxwQhvyYLERkmIutFJEVE7ipm/7UislJElonIXBHp6GxPEJEsZ/syEXnRn3GammVnRhZPfrWeQe1iOa9r1c4k60/3juhIVHgId3+80gbrmUrnt2QhIqHARGA40BG4qDAZeHlHVbuo6snA48DTXvs2qerJzu1af8Vpap4HZ6whr0B5eFTnajWdx4mKjYnkH+d04OfNB/hwyXa3wzFBxp81iz5AiqqmqmoOMBUY5X2Aqh72elgHsK9Dxq++WbObL1fv4uYhycQ3qu12OJXuwl6t6NW6AY98voYDx3LcDscEEX8mi5aA91wE6c623xCR60VkE56axU1euxJFZKmIzBaRgcU9gYiMF5FFIrJo7969lRm7CULZuflM+Gw17ZpGc/XAJLfD8YuQEOFf53fhSHYeT3293u1wTBDxZ7Iorn7/u5qDqk5U1TbAncA9zuadQLyqdgduA94RkbrFlJ2kqr1UtVdsbGwlhm6C0eR5m9l2IIv7z+tEeGjw9u1o1zSGS05pzbu/bGXtzsNlFzDGB/78j0kHWnk9jgN2lHL8VOAPAKp6XFX3O/cXA5uAdn6K09QAew5nM/H7FIZ0aMqAto3dDsfvbhmSTN1a4Tz02Rq72G0qhT+TxUIgWUQSRSQCGAPM8D5ARJK9Hp4LbHS2xzoXyBGRJCAZSPVjrCbIPfn1enLyC7j73A5uh1Il6teO4NYh7Zi/aT/frNntdjgmCPgtWahqHnAD8BWwFpimqqtFZIKIjHQOu0FEVovIMjzNTeOc7YOAFSKyHPgAuFZVD/grVhPcVm3P4P3F6Vw+IJHExnXcDqfKXNw3nuQm0Twycy3H82xktzkxYf48uarOBGYW2Xaf1/2bSyj3IfChP2MzNYOqMuHTNTSsHcENZ7R1O5wqFRYawr0jOvKXyb8wZV4a11STJWJNYAreq3zGAN+u3cMvaQe47ex21I0K/JHa5TWoXSynt49l4g8pZGTmuh2OCWCWLEzQyi9QnvhqHUmN63Bhr1ZlFwhSdw4/iSPH83hhdorboZgAZsnCBK2PlqSzYfdRbh/anrAg7ipblpOa1eWPJ7dkyrw0dmZkuR2OCVA19z/IBLXs3Hye+XYjXePqMbxzM7fDcd2tZ7WjQJX/frfR7VBMgLJkYYLSWwu2sP1QFncOOymo5n+qqFYNa3Nx39ZMW5TOpr02K60pP0sWJugczs5l4g8pDExuXCMG4PnqhjPaEhUWwpNf2TQgpvwsWZig8+qczRzMzOXvQ9u7HUq10jg6kqsHJfHFql2sTM9wOxwTYCxZmKCSkZnL5LmbGdqpKV3j6rsdTrVz5amJ1K8dzjPfbnA7FBNgLFmYoPLq3FSOHM/jliE2lVhxYqLCuXpgEt+t28PybYfcDscEEEsWJmgcysxh8rw0hnduRofmv5uk2DjG9U+ggdUuTDlZsjBB4+U5qRw9nsfNQ5LLPrgGi44M4+pBSfywfi9Lth50OxwTICxZmKBw4FgOU+alcW7X5pzUzGoVZRnXL4GGdSJ45lsbd2F8Y8nCBIWX56SSmZvPLWdarcIXdSLDGD8oiR837GXxFpvQ2ZTNkoUJeIcyc3hjfhrndmlOctMYt8MJGH/p15qGdSJ47nubM8qUzadkISK1RMQ6rZtq6bV5aRzLya9xU5CfqNoRYVx5aiKz1u+1cRemTGUmCxE5D1gGfOk8PllEZpReypiqcSQ7lynz0zirY1O7VlEBl/ZrTUxUGBN/sNqFKZ0vNYsHgD7AIQBVXQYk+C8kY3z31oKtZGTlcsPpVquoiLpR4VzeP4EvV+9iw+4jbodjqjFfkkWeqlaojioiw0RkvYikiMhdxey/VkRWisgyEZkrIh299v3DKbdeRIZW5PlNcMvKyeeVOakMTG5Mt1Y2WruiLh+QSO2IUKtdmFL5kixWichYIFREkkXkOWB+WYVEJBSYCAwHOgIXeScDxzuq2kVVTwYeB552ynYExgCdgGHAC875jPnV1IVb2X8shxvPsB5QJ6JBnQguOaU1ny7fQdq+Y26HY6opX5LFjXg+tI8D7wKHgVt8KNcHSFHVVFXNAaYCo7wPUNXDXg/rAOrcHwVMVdXjqroZSHHOZwwAx/PyeWl2Kn0SGtInsaHb4QS8q05NJCw0hP/N2uR2KKaaKjNZqGqmqt6tqr1VtZdzP9uHc7cEtnk9Tne2/YaIXC8im/DULG4qZ9nxIrJIRBbt3bvXh5BMsPhk6XZ2Hc7meusBVSma1I3iwl6t+Ghpuq2mZ4rlS2+oH0Tk+6I3H85d3Ioz+rsNqhNVtQ1wJ3BPOctOchJYr9jYWB9CMsGgoEB56cdUOrWoy6BkW6+isowflESBwuS5m90OxVRDYT4cc7vX/SjgT0CeD+XSgVZej+OAHaUcPxX4XwXLmhrk6zW7Sd17jOcu6m6r4FWiVg1rM6Jrc975eSs3nJ5MvdrhbodkqhFfmqEWe93mqeptQF8fzr0QSBaRRBGJwHPB+jfjM0TE+8rkuUDhRDUzgDEiEikiiUAy8IsPz2mCnKryv9mbiG9Y29bW9oNrBrXhWE4+by5IczsUU8340gzV0OvW2OnGWuZ/qarmATcAXwFrgWmqulpEJojISOewG0RktYgsA24DxjllVwPTgDV4BgNer6r5FXmBJrgsSD3A8m2HGD8oibBQm62msnVsUZfB7WN5bV4a2bn2L2f+ny/NUIvxXC8QPM1Pm4ErfTm5qs4EZhbZdp/X/ZtLKfsI8Igvz2Nqjhdnb6JxdAQX9IxzO5Sgde1pbRgzaQHvL9rGpf0S3A7HVBNlJgtVTayKQIwpy5odh5m9YS9/H9qeqHAbduMvfRMbcnKr+kyak8pFfeKtBmeAUpKFiJxfWkFV/ajywzGmZJN+3ESdiFAu6dva7VCCmojw18FtuObNxcxctYuR3Vq4HZKpBkqrWZxXyj4FLFmYKrP9UBafrtjJ5f0TrJdOFTirQ1MSG9fh5R9TOa9rc+t1ZkpOFqp6eVUGYkxpXnP6/l9+qrWKVoWQEOGqgYnc/fEqFqQeoF+bRm6HZFzmywVuRORcPFN+RBVuU9UJ/grKGG8ZWbm8+8tWzuvanJb1a7kdTo3xpx5xPP31Bl6ek2rJwvjUdfZF4EI8c0QJMBqwRmNTZd79ZSvHcvK5elCS26HUKFHhofylXwLfr9vDRpu+vMbzpZtDf1X9C3BQVR8E+vHb0dXG+E1OXgGvzdvMqW0b06lFPbfDqXEu7deaqPAQXp6T6nYoxmW+JIvCWcUyRaQFkAtYw7GpEjOW72D34eNWq3BJwzoRjO7Zik+W7mDPYV/mDzXBypdk8ZmI1AeeAJYAaXimKjfGr1SVl39M5aRmMTZhoIuuGphIbkEBU+anuR2KcZEvc0M9pKqHVPVDPNcqTvIehW2Mv/y4cR/rdx/hqoFJ1nXTRa0b1WFox2a8/fNWMnN8mUPUBCNfLnAvF5F/ikgbZzGiCi2xakx5vTInlSYxkTYorBq4elAiGVm5fLA43e1QjEt8aYYaiWdOqGkislBEbheReD/HZWq4dbsOM2fjPsb1TyAizKabcFvP1g3pHl+fV+duJr/gd0vLmBrAl2aoLar6uKr2BMYCXfFMJmiM37wyZzO1wkO5uK99L6kurjo1iS37M/lmzW63QzEu8Okrm4gkiMgdeBYoOgm4w69RmRptz+Fspi/bzuhecdSvHeF2OMYxtFNT4hrU4tW51o22JvLlmsXPeOaBCgFGq2ofVX3K75GZGuuNn7aQV6BcMcB6aFcnYaEhXDEgkYVpB1m27ZDb4Zgq5kvNYpyq9lDVR1XVvlIYv8rMyeOtn7dwVoemJDSu43Y4pog/925FTFSYDdKrgXy5ZrGuKgIxBuDDJds5lJnLVQNtEF51FB0Zxtg+8XyxcifbDmS6HY6pQn7tZiIiw0RkvYikiMhdxey/TUTWiMgKEflORFp77csXkWXObUbRsib4FBQok+dupltcPXonNHA7HFOCcf0TEBFet0F6NYrfkoWIhAITgeFAR+AiEelY5LClQC9V7Qp8ADzutS9LVU92biMxQe+H9XvYvO8YV9ogvGqtRf1anNulOe8t3MaR7Fy3wzFVxJcL3A+JSJjX47oi8poP5+4DpKhqqqrm4OlJNcr7AFX9QVUL67ILAFtYuQZ7Zc5mWtSLYnjnZm6HYspw1cBEjhzPY9oiG6RXU/hSswgDfhaRriJyNrAQWOxDuZbANq/H6c62klwJfOH1OEpEFonIAhH5Q3EFRGS8c8yivXv3+hCSqa6OHc/jp9T9jOufQLit+VztdY2rT++EBrw2bzNqY/RqhDIXP1LVf4jId8DPwEFgkKqm+HDu4toRin1bicglQC/gNK/N8aq6Q0SSgO9FZKWqbioS2yRgEkCvXr3sLRvAdmVkExMRypg+NggvUFx5ahLXvrWYutE5NKxj42GCnS/NUIOAZ4EJwCzgeWeq8rKk89t1L+KAHcWcfwhwNzBSVY8XblfVHc7PVOd5u/vwnCYA5eQXsO9YDn/u1Yp6tWx97UBxVsemxDeszc6MrLIPNgHPl/r+k3gG4/1bVcfi+Sb/vQ/lFgLJIpIoIhHAGOA3vZpEpDvwEp5EscdrewMRiXTuNwYGAGt8eUEm8OzOyEaxQXiBJjREuGJAAkey8zh63GajDXa+JIt+qvrrB7WqfoTnw7tUqpoH3AB8BawFpqnqahGZICKFvZueAKKB94t0ke0ALBKR5cAPwKPeMZjgkZmTx+4jx2lYO4L4RrXdDseU0+herQgNEXZm2MJIwc6Xaxb5xWzb78vJVXUmMLPItvu87g8podx8oIsvz2EC24dLtpOXX0DzerXcDsVUQJ3IMJrWjWJnRhbbDmTSqqEl/GBl3U6MawoKlNfmbiY6MoyYqDK/t5hqqlndKMAG6QW7UpOFiISISP+qCsbULD+s30PqvmNWqwhwEWEhNKoTYYP0glypyUJVCwCbYdb4ReEgPOt2Gfia14uyQXpBzpdmqK9F5E9i8y+YSrR6R8avg/DsnRX46kSG/TpILy+/wO1wjB/4kixuA94HckTksIgcEZHDfo7LBLlX526mtg3CCypXnppE+sEsvraV9IKSL1OUx6hqiKqGq2pd53HdqgjOBKc9h7P5dPkOG4QXZAoH6b1ia10EJV+XVR0pIk86txH+DsoEt9d/SrOV8IJQ4SC9JVsPsXjLQbfDMZXMl+k+HgVuxjOCeg1ws7PNmHLLzMnjrQVbGdqxmQ3CC0Kje7WiblSYrdMdhHypWZwDnKWqk1V1MjDM2WZMuX2wOJ2MrFyuHmS1imBUJzKMsX1b8+WqXbaSXpDxdVBefa/79fwRiAl++c5KeN3j69OzdUO3wzF+cln/BEJEmDxvs9uhmErkS7L4N7BURKaIyOt41rL4l3/DMsHo27W7SdufyVWn2vrawaxZvSjO69aCaQu3kZFlg/SCRYnJQkQKJwv8CDjF+fkRnokFp1ZBbCbIvDInlbgGtRjaqanboRg/u/LURI7l5DP1l61uh2IqSWk1i/86P39S1Z2qOkNVp6vqrqoIzASXZdsOsTDtIJcPSCTMVsILep1b1qNfUiOmzE8j1wbpBYXS/mtznbW2W4rIf4veqipAExxenpNKTGQYF/ZuVfbBJihcPSiRnRnZfL5ip9uhmEpQWrIYgWctimw81ymK3ozxybYDmXyxcidjT4knOtJml60pBrdrQtsm0Uz6MRW1hboDXon/uaq6D5gqImtVdXkVxmSCzKtzNxMiwuX9rbtsTRISIlw9MJE7P1zJ/E37GdC2sdshmRPgy3QflihMhR3KzGHaom2MPLkFzepFuR2OqWKjTm5J4+hIJv1og/QCnV+vNIrIMBFZLyIpInJXMftvE5E1IrJCRL4TkdZe+8aJyEbnNs6fcRr/efvnrWTm5DN+kHWXrYmiwkO5fEACszfsZf2uI26HY06AL9N9hFbkxE65icBwoCNwkYh0LHLYUqCXqnYFPgAed8o2BO4H+gJ9gPtFpEFF4jDuOZ6Xz5T5aQxqF8tJzWzuyZrq4r7x1AoP5WWbYDCg+VKzSBGRJ4r5oC9LHyBFVVNVNQeYCozyPkBVf1DVwjkBFgBxzv2hwDeqemTecTcAACAASURBVEBVDwLf4JlmxASQ6Ut3sPfIccYPtFpFTVa/dgQX9m7F9GXb2ZWR7XY4poJ8SRZdgQ3AKyKyQETGi4gvXxNbAtu8Hqc720pyJfBFeco6sSwSkUV79+71ISRTVQoKlJfnpNKheV0GtG3kdjjGZVcMSCS/QHltvk0BEqh8ucB9RFVfVtX+wB14mod2isjrItK2lKLFrX9WbP85EbkE6AU8UZ6yqjpJVXupaq/Y2NhSX4epWt+v28PGPUe5ZlAStsiiiW9Um+FdmvPOgq0ctnW6A5JP1yyc9Sw+Bp7FsyZ3EvApMLOUoumA9wisOGBHMecfAtwNjFTV4+Upa6qvl37cRMv6tTi3a3O3QzHVxLWD2nDkeB7v/mxTgAQiX5qhNuK51vCEqnZX1adVdbeqfgB8WUq5hUCyiCSKSAQwBpjhfYCIdAdewpMo9njt+go4W0QaOBe2z3a2mQCweMsBFqYd5KqBiYTb1B7G0SWuHgPaNuLVuZs5npfvdjimnHz5T/6Lql6pqvMLNxROMqiqN5VUSFXzgBvwfMivBaap6moRmSAiI53DngCigfdFZJmIzHDKHgAewpNwFgITnG0mALw4O5X6tcNtag/zO9ee1oY9R44zfak1FAQaX+Ze+C/Qo8i254rZ9juqOpMiTVWqep/X/SGllJ0MTPYhPlONpOw5wjdrdnPTmcnUjrCpPcxvndq2MZ1a1OXFHzdxQc84QkLselagKPG/WUT6Af2BWBG5zWtXXaBCYy9M8Jv0YypR4SGM69e67INNjSMiXHNaG256dynfrt3N2Z2auR2S8VFpzVAReJqIwoAYr9th4AL/h2YCza6MbD5eup0/92pFo+hIt8Mx1dQ5nZvRqmEtXpy9ySYYDCClTSQ4G5gtIlNUdUsVxmQC1CtzUilQuNoG4ZlShIWGMH5gEvdOX83Pmw9wSpKNwwkEpa2U94xz93kRmVH0VkXxmQBx8FgO7/yylZHdWtCqYW23wzHV3OherWgcHcnEH1LcDsX4qLQrkG86P5+sikBMYJsyP43MnHz+OriN26GYABAVHsqVpyby2JfrWJmeQZe4em6HZMpQYs1CVRc7P2cXd6u6EE11d/R4HlPmp3FWx6a0axrjdjgmQFxySjwxUWG8MMtqF4GgtN5QKylheg4AZ6ZYY3jn5y1kZOVyndUqTDnERIVzWf8Env8hhZQ9R2jbxL5oVGelNUONqLIoTMDKzs3n5TmbGdC2Ed3jbRZ5Uz6X9U/g5Tmp/G9WKk/9uZvb4ZhSlNYMtaW0W1UGaaqvD5eks/fIca4fXNqcksYUr1F0JBf1iWf6su2kH8wsu4BxTWm9oeY6P4+IyOGiP6suRFNd5eYX8MIPmzi5VX36tbHuj6Zirh6YhAi8OHuT26GYUpRWszjV+RmjqnWL/qy6EE119fGS7Ww/lMVNZ7a1achNhbWoX4sLerZi2sJ0dmZkuR2OKYFPU4KKSA8RuUlEbnRmijU1XF5+ARNnpdC5ZV1Ob9/E7XBMgLtucBsKVHlpti29Wl35sp7FfcDrQCOgMTBFRO7xd2CmepuxfAdb9mdy0xnJVqswJ6xVw9qc36Ml7/6ylT2HbenV6siXmsVFQG9VvV9V7wdOAS72b1imOssvUJ7/PoUOzetyVsembodjgsT1p7clr0CZ9KPVLqojX5JFGhDl9TgSsCtRNdjnK3eSuu8YN55h1ypM5WndqA6jurXgrZ+3sO/o8bILmCpVWm+o50Tkv8BxYLWITBGR14BVwNGqCtBULwUFynPfbSS5STTDbHppU8muP6Mtx/MKeNlqF9VOaYPyFjk/FwMfe22f5bdoTLX32cqdbNxzlP9e1N0WrjGVrk1sNCO7teCNn7Zw1cAkYmNsqvvqorQpyl8/0ZOLyDDgWTyLJb2iqo8W2T8IeAboCoxx1vUu3JcPrHQeblXVkRhX5eUX8My3G2jfNIYRXZq7HY4JUjedmcyny3fw4uxN3Duio9vhGIcvvaGSReQDEVkjIqmFNx/KhQITgeFAR+AiESn6l98KXAa8U8wpslT1ZOdmiaIamL5sB6l7j3HrWclWqzB+0yY2mj92j+OtBVvYbT2jqg1fLnC/BvwPyANOB97g/6cvL00fIEVVU1U1B5gKjPI+QFXTVHUFUFCuqE2Vy80v4L/fb6RTi7oMtWsVxs9uPjOZvALlBVvvotrwJVnUUtXvAHHmhXoAOMOHci2BbV6P051tvooSkUUiskBE/lDcASIy3jlm0d69e8txalNeHy1JZ8v+TG4d0s56QBm/i29Um9E943j3l23sOGSjuqsDX5JFtoiEABtF5AYR+SPgy5Dd4j5RyrPgbryq9gLGAs+IyO/mv1bVSaraS1V7xcbGluPUpjxy8gr473cpdIurx5kdbLS2qRo3nNEWRXneahfVgi/J4hagNnAT0BO4FBjnQ7l0oJXX4zhgh6+BqeoO52cqnh5YNs2IS95buJXth7K49SyrVZiqE9egNhf2bsW0hdvYut9mpHVbmclCVReq6lHgMHCTqp6vqgt8OPdCIFlEEkUkAhgD+LR2t4g0EJFI535jYACwxpeypnJl5uTx7Hcp9EloyGntrPZmqtaNZyQTFio8/c16t0Op8XzpDdXLWTVvBbBSRJaLSM+yyqlqHnAD8BWwFpimqqtFZIKIjHTO3VtE0oHRwEsistop3gFYJCLLgR+AR1XVkoULJs/dzL6jx7lzeHurVZgq17RuFJcPSGT68h2s2WErI7iptEF5hSYD16nqHAARORVPD6kyl1VV1ZnAzCLb7vO6vxBP81TRcvOBLj7EZvzo4LEcXpqdypAOTenZuqHb4Zga6tpBbXh7wRae+Godr13ex+1waixfrlkcKUwUAKo6Fzjiv5BMdfG/2Zs4mpPHHcPaux2KqcHq1Q7nutPb8sP6vfycut/tcGqs0uaG6iEiPYBfROQlERksIqeJyAvYlB9Bb8ehLKbMT+P87nG0axrjdjimhhvXL4GmdSN57Mt1qJanU6WpLKU1Qz1V5PH9XvftrxXknvl2Ayjcelay26EYQ62IUG4Z0o5/fLSSb9bs5mwbGFrlSpsb6vSqDMRUH2t3Hub9xelcMSCRuAa13Q7HGABG94zjlTmpPPrFOk4/qQnhoT4t9GkqiS+9oeqJyNOFI6VF5CkRqVcVwZmqp6o88vla6kaFc+MZbd0Ox5hfhYWG8M9zOpC67xhvL9jidjg1ji+peTKeC9p/dm6H8fSGMkFo1vq9zE3Zx81nJlO/doTb4RjzG2ec1IT+bRrx7HcbycjKdTucGsWXZNHGWVI11bk9CCT5OzBT9fLyC3hk5loSGtXmklNaux2OMb8jItx9bgcOZeUy0aYBqVK+JIssZ2wFACIyALCZvYLQ1IXbSNlzlLuGdyAizNqDTfXUqUU9LugRx5R5aTYNSBXy5RPhWmCiiKSJSBrwPHCNX6MyVe5wdi7/+WYDfRIaMrRTU7fDMaZUtw9tT2iI8OiXa90OpcYoNVk4s822V9VueEZsd1XV7s4aFCaIPPvtRg5k5nDPiA42rYep9prWjeLa09owc+Uu5qfsczucGqHUZKGqBXjmd0JVD6uqTc4ShDbsPsKU+WmM6d2KrnH13Q7HGJ9cc1oScQ1q8cCnq8nNt/XT/M2XZqhvROR2EWklIg0Lb36PzFQJVeWBGauJjgzj70NPcjscY3wWFR7KvSM6smH3Ud78ybrS+psvEwle4fy83mubYj2igsIXq3Yxf9N+JozqRMM61lXWBJazOzZlULtY/vPNBs7r1oLYmEi3QwpavqxnkVjMzRJFEMjMyePhz9bQoXldxvaJdzscY8pNRLj/vI5k5+Xz+Jfr3A4nqPkygjtKRG4TkY9E5EMRuUVEoqoiOONfE39IYUdGNg+O7ESYTZ1gAlSb2GiuGJDI+4vTWbzlgNvhBC1fPiHeADoBz+HpNtsReNOfQRn/W7/rCC/NTuX8Hi3pk2iXoExgu/HMZFrUi+KfH62yi91+4kuyaK+qV6rqD85tPNDO34EZ/ykoUP758UpiosK459yObodjzAmLjgzjwVGdWb/7CC/PSXU7nKDkS7JYKiKnFD4Qkb7APF9OLiLDRGS9iKSIyF3F7B8kIktEJE9ELiiyb5yIbHRu43x5PuObdxduZfGWg/zznA52UdsEjbM6NmVop6Y8++1Gtuw/5nY4QceXZNEXmO81gvsn4DQRWSkiJQ7OE5FQYCIwHE/T1UUiUvRr7FbgMuCdImUb4lk/oy/QB7hfRBr49IpMqfYczubRL9bRL6kRF/T83Yq2xgS0B0d2Jjw0hHs+WWWLJFUyX7rODqvgufsAKaqaCiAiU4FRwJrCA1Q1zdlXtJFxKPCNqh5w9n/jxPFuBWMxjgc/W8Px3AIe+WNnG6ltgk6zelH8fWh77p+xmunLdvCH7i3dDilolJksVLWio11aAtu8HqfjqSlUtOzv/uoiMh4YDxAfb10/y/Llqp18vmInt53VjqTYaLfDMcYvLjmlNR8v3c4Dn66mf9tGNImxzpuVwZ/9JYv72uprvdCnsqo6SVV7qWqv2NjYcgVX0xw4lsM9n6yiU4u6/HVwG7fDMcZvQkOEJ0d3JTMnn3s+tuaoyuLPZJEOtPJ6HAfsqIKyphj3TV9FRlYuT47uZstRmqDXtkkMfzurHV+v2c2M5fbRURn8+amxEEgWkUQRiQDGADN8LPsVcLaINHAubJ/tbDMVMHPlTj5bsZObzkimQ/O6bodjTJW4amAS3ePrc9/01ew5nO12OAHPb8lCVfPwzFj7FbAWmKaqq0VkgoiMBBCR3iKSDowGXhKR1U7ZA8BDeBLOQmBC4cVuUz77jx7n3k9W0bllXa615idTg3iao7qRnZvPP6056oT50huqwlR1JjCzyLb7vO4vxNPEVFzZyXjW/zYVpKrc8cEKjmTn8bY1P5kaqE1sNH8f2p6HP1/L1IXbuMjmQKsw+/QIYm8u2MJ36/bwj3NO4qRm1vxkaqYrBiQyMLkxD366mpQ9R90OJ2BZsghS63cd4eHP1zK4fSyX9U9wOxxjXBMSIjw1uhu1I8K46d2lHM/LdzukgGTJIghl5+Zz07tLqRsVxhMXdLPBd6bGa1I3isf+1JU1Ow/zxJfr3Q4nIFmyCEKPfL6W9buP8OTobrYYjDGOszo25dJTWvPK3M38sH6P2+EEHEsWQWb6su28uWALV52ayOD2TdwOx5hq5e5zO3BSsxhufW8Z6Qcz3Q4noFiyCCIbdh/hrg9X0juhAXcOt/W0jSkqKjyU/13Sk/x85bq3l9j1i3KwZBEkjh7P49q3FlMnMoznx/awbrLGlCCxcR2e/HM3VqRn8NBna8ouYABLFkFBVbnzgxVs2Z/J82O707SuTZxmTGmGdmrGNacl8daCrXy0JN3tcAKCJYsg8Pz3KXy+cid3DG3PKUmN3A7HmIDw97Pbc0pSQ+76aCVLtx50O5xqz5JFgJu5cidPfbOB87u3ZPygJLfDMSZghIWG8MLFPWlWN4qr31jMjkNZbodUrVmyCGAr0g9x27Rl9GzdgH//qYuNpzCmnBrWieDVcb04npvPla8v4tjxPLdDqrYsWQSoXRnZXP3GIhrVieSlS3sSGRbqdkjGBKTkpjE8N7Y763cd5pb3lpFfYBMOFseSRQDKyMxl3ORfOJqdxyvjetE42gbeGXMiBrdvwn0jOvLNmt08MGO1zVBbDL/OOmsqX3ZuPle9sZDN+47x2uW9bX0KYyrJZQMS2ZmRzUs/ptI4OpKbhyS7HVK1YskigOTlF3DDO0tYtOUgz1/UgwFtG7sdkjFB5a7hJ7HvaA7/+XYDjaIjuOSU1m6HVG1YsggQBQXKXR+t5Nu1e3hoVCfO7drc7ZCMCToiwqN/6sLBzBzunb6K+rXDGdG1hdthVQt2zSIAFBQo//x4JR8sTueWIclc2i/B7ZCMCVrhoSFMHNuDXq0bcPPUZXyxcqfbIVULfk0WIjJMRNaLSIqI3FXM/kgRec/Z/7OIJDjbE0QkS0SWObcX/RlndVZQoNz9yUqmLtzGjWe05eYzrR3VGH+rFRHKa5f34eRW9bnx3aV8ucoSht+ShYiEAhOB4UBH4CIR6VjksCuBg6raFvgP8JjXvk2qerJzu9ZfcVZnBQXKPdNX8e4v27j+9DbcdlY7G0thTBWJjgxjyuW96RpXjxvesYThz5pFHyBFVVNVNQeYCowqcswo4HXn/gfAmWKfhgDk5hfwt/eX887PW7lucBtuP7u9JQpjqlhMVDivX9GHrnH1uP6dpXywuObOI+XPZNES2Ob1ON3ZVuwxqpoHZACFkxslishSEZktIgOLewIRGS8ii0Rk0d69eys3ehdl5eQz/o1FfLx0O38f2p6/D7VEYYxbYqLCeePKvpyS1JDb31/Oyz+muh2SK/yZLIr7dCs60qWkY3YC8araHbgNeEdEfjegQFUnqWovVe0VGxt7wgFXB4cyc7jk1Z+ZvWEv//pjF64/va0lCmNcFh0ZxuTLenNul+Y8MnMtj36xrsYN3PNn19l0oJXX4zhgRwnHpItIGFAPOKCev8JxAFVdLCKbgHbAIj/G67rUvUe56vVFpB/MYuLYHgzvYt1jjakuIsNC+e9F3WlQJ5wXZ28i/WAmT47uRlR4zZhqx581i4VAsogkikgEMAaYUeSYGcA45/4FwPeqqiIS61wgR0SSgGQgqOt+czfu4w8T53EoK5e3rupricKYaig0RHhoVGfuGn4Sn6/cyZ9f+ondh7PdDqtK+C1ZONcgbgC+AtYC01R1tYhMEJGRzmGvAo1EJAVPc1Nh99pBwAoRWY7nwve1qnrAX7G6SVV5fX4a4177heb1ajH9+gH0SWzodljGmBKICNee1oZJl/Zi056jjHx+Lsu2HXI7LL/z6whuVZ0JzCyy7T6v+9nA6GLKfQh86M/YqoMj2bn88+NVfLp8B2ee1IRnL+pOdKQNqjcmEJzVsSkfXtefK6csYvSL87n7nA6M658QtNcYbQS3S1bvyGDk8/P4fMUO/j60PS//pZclCmMCzEnN6vLZjacyKDmWBz5dw1/fWkJGVq7bYfmFJYsqVlCgvDp3M398YT6ZOXm8e/UpXH96W0JCgvPbiDHBrkGdCF4Z14u7z+nAt2t3c+5/57Agdb/bYVU6SxZVaMv+Y4yZtICHPlvDwLaNmXnTQPramtnGBDwR4epBSbx3TT9CQ4SLXl7Ag5+uJisn3+3QKo21e1SBvPwC3lywhce/XE9YqPDk6G78qUfLoG3bNKam6tm6AV/cPJDHvljHa/PSmLV+L/8+vwunBMGXQqtZ+NniLQcY+fw8Hvx0DX0SG/L1rYO4oGecJQpjglTtiDAeHNWZd67qS25+AWMmLeCWqUvZcySwu9hazcJPdh/O5smv1vP+4nSa1Y3ihYt7MLxzM0sSxtQQ/ds25ptbT+OFWSm8NDuV79bu4eYhyVzarzWRYYE3kM+SRSXLyMrlpdmbmDxvM/kFyjWDkrjpzGTqWE8nY2qcWhGh/O3s9pzfI477Z6zm4c/XMmV+Gn87ux2jurUMqI4t9glWSQ5n5/LmT1t4eU4qhzJzGXVyC/52VnviG9V2OzRjjMsSG9fhjSv6MGfjXh79Yh23vrecl2anctOZyQzt1IzQAEgalixO0L6jx3lt3mbemL+FI8fzGNw+ltvPbk/nlvXcDs0YU80MTI5lQJvGfLZyJ//5ZgPXvb2EpNg6/PW0Now6uSURYdX3MrIliwpakX6I1+dv4dMVO8jNL2B452ZcN7itJQljTKlCQoSR3VpwbpfmfLFqJxN/2MTfP1jBE1+tZ2zfeMb2iadJ3Si3w/wdSxblkJGVy+crdvLeom0s33aI2hGh/LlXHJf1T6Rtk2i3wzPGBJDQEGFEV0/SmLVhL6/PT+OZbzfy/PcpDO3UjAt6xTGwbWPCQqtHbcOSRRmyc/OZs3Ef05dt5+s1u8nJK6Btk2geOK8j5/eMo25UuNshGmMCmIhwevsmnN6+CWn7jvHmgi18uCSdz1fuJDYmkj+c3IIRXVvQNa6eq70pLVkUIyMrlzkb9/Llql38sG4Px3LyaVA7nIt6t+JPPePo0tLdP5oxJjglNK7DvSM6cuewk/h+3R4+WpLOlPlpvDxnMy3r12Jop2ac3akpPVs3ILyKaxyWLPCsd70iPYOfNu1j1vq9LN12iPwCpVGdCEae3JLhnZvRr02jKv/jGGNqpoiwEIZ1bsawzs04lJnDt2v38OWqnby1YAuT520mOjKMAW0bcVq7JpyS1JDExnX8/gW2RiaLg8dyWJZ+iOXbDrEo7SCLtxwkK9czh0uXlvW4bnAbTmsXS/f4BgHRpc0YE7zq147ggp5xXNAzjiPZucxL2c/sDXuZvX4PX63eDUBsTCR9ExvSI74BJ8fXp2PzupW+gl9QJ4v8AiX9YCYbdh9l7c7DrN15mDU7D7NlfyYAItC+aQwX9m5F38SG9ElsSKPoSJejNsaY4sVEhf9a41BVUvcd4+fUA/y8eT+/bD7AZyt2AhAeKrRvFkOHZnXp0Nxza9skmsbRERWugQRNsjh6PI+3f97C1v2ZbD2QSereY2zef4ycvALAkxhaN6xNx+Z1GdM7nm6t6tE1rr6tIWGMCUgiQpvYaNrERjO2bzwAuzKyWbbtEMu2HWL1jgx+WL+H9xen/1qmblQYbZpEk9ioDvGNahPf0PdBw379pBSRYcCzQCjwiqo+WmR/JPAG0BPYD1yoqmnOvn8AVwL5wE2q+lVpz7V53zHu/ngVEaEhxDWsRVLjOgxuH0ub2GjaNo2mfdMYm3LDGBPUmtWLYlg9T82j0J4j2azfdYSUPUfZtPcoKXuOsiB1Px8v246q7+f226eniIQCE4GzgHRgoYjMUNU1XoddCRxU1bYiMgZ4DLhQRDoCY4BOQAvgWxFpp6olTg6f1LgO3//jDJrGRAXUfCvGGONPTWKiaBITxcDk2N9sP56XT/rBLNo+5tt5/Nm9pw+QoqqpqpoDTAVGFTlmFPC6c/8D4EzxNKiNAqaq6nFV3QykOOcrUZ3IMJrXq2WJwhhjfBAZFkqbWN8HE/uzXaYlsM3rcTrQt6RjVDVPRDKARs72BUXKtiz6BCIyHhjvPDwuIqsqJ/RqqTGwz+0g/KixXC5B/fqwv1+gCva/XXtfDvJnsijuK37RFrKSjvGlLKo6CZgEICKLVLVXeYMMFPb6Apu9vsAVzK8NPK/Pl+P82QyVDrTyehwH7CjpGBEJA+oBB3wsa4wxpor4M1ksBJJFJFFEIvBcsJ5R5JgZwDjn/gXA96qqzvYxIhIpIolAMvCLH2M1xhhTCr81QznXIG4AvsLTdXayqq4WkQnAIlWdAbwKvCkiKXhqFGOcsqtFZBqwBsgDri+tJ5Rjkr9eSzVhry+w2esLXMH82sDH1ydano62xhhjaiSbGc8YY0yZLFkYY4wpU1AlCxF5SERWiMgyEflaRFq4HVNlEpEnRGSd8xo/FpH6bsdUmURktIisFpECEQmKrooiMkxE1otIiojc5XY8lU1EJovInmAc4yQirUTkBxFZ67wvb3Y7psokIlEi8ouILHde34OlHh9M1yxEpK6qHnbu3wR0VNVrXQ6r0ojI2Xh6jOWJyGMAqnqny2FVGhHpABQALwG3q6pP/b+rK2fKmw14TXkDXFRkypuAJiKDgKPAG6ra2e14KpOINAeaq+oSEYkBFgN/CJa/nzNbRh1VPSoi4cBc4GZVXVDc8UFVsyhMFI46FDOQL5Cp6teqmuc8XIBn/EnQUNW1qrre7TgqkS9T3gQ0Vf0RT0/GoKOqO1V1iXP/CLCWYmaSCFTqcdR5GO7cSvzMDKpkASAij4jINuBi4D634/GjK4Av3A7ClKq4KW+C5sOmJhGRBKA78LO7kVQuEQkVkWXAHuAbVS3x9QVcshCRb0VkVTG3UQCqereqtgLeBm5wN9ryK+v1OcfcjWf8ydvuRVoxvry+IOLTtDWmehORaOBD4JYirRcBT1XzVfVkPK0UfUSkxKbEgFvgQVWH+HjoO8DnwP1+DKfSlfX6RGQcMAI4UwPwglM5/n7BwKatCXBOW/6HwNuq+pHb8fiLqh4SkVnAMKDYzgoBV7MojYgkez0cCaxzKxZ/cBaTuhMYqaqZbsdjyuTLlDemmnIuAL8KrFXVp92Op7KJSGxhj0oRqQUMoZTPzGDrDfUhnul2C4AtwLWqut3dqCqPMy1KJJ5VBQEWBFlvrz8CzwGxwCFgmaoOdTeqEyMi5wDP8P9T3jzickiVSkTeBQbjmcZ7N3C/qr7qalCVREROBeYAK/F8pgD8U1VnuhdV5RGRrnjWEwrFU3GYpqoTSjw+mJKFMcYY/wiqZihjjDH+YcnCGGNMmSxZGGOMKZMlC2OMMWWyZGGMMaZMliyMcYGIHC37KGOqD0sWxhhjymTJwphSiEhvZ/2QKBGp48z737nIMY+JyHVejx8Qkb+JSLSIfCciS0RkZXHzX4nIYBH5zOvx8yJymXO/p4jMFpHFIvKVM2W2Ma6wZGFMKVR1IZ4pOh4GHgfeUtWic+dMBS70evxn4H0gG/ijqvYATgeecqaQKJMzJ9FzwAWq2hOYDATV6G8TWAJuIkFjXDABzzxP2cBNRXeq6lIRaeKszBgLHFTVrc4H/r+cBYIK8ExP3hTY5cNztgc6A984+SUU2FkZL8aYirBkYUzZGgLReBaHiQKOFXPMB8AFQDM8NQ3wrKkSC/RU1VwRSXPKe8vjtzX8wv0CrFbVfpXxAow5UdYMZUzZJgH34lk/5LESjpmKZ1bZC/AkDoB6wB4nUZwOtC6m3Bago4hEikg94Exn+3ogVkT6gadZSkQ6VcqrMaYCrGZhTClE5C9Anqq+46ypPV9EzlDV772PU9XVzjrN21W1sLnobeBTEVkELKOY6Z9VdZuITAP+r507xAEQioEof7z2WwAAAFBJREFUuD/hiJyMi5GQoLgJCsHXrMLNyLqql1T0SHIl2ef8HmOsSbYZkSXv99rzn03hm6+zAFTOUABUYgFAJRYAVGIBQCUWAFRiAUAlFgBUD1rJQvzaZNPeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(start=-3, stop=3, num=100)\n",
    "_ = plt.plot(x, norm.pdf(x))\n",
    "plt.xlim(-3,3)\n",
    "plt.ylim(0)\n",
    "plt.xlabel('x value')\n",
    "plt.ylabel('probability for x value')\n",
    "plt.title('Probability Distribution Function')\n",
    "_ = plt.axvline(-1, color='g')\n",
    "_ = plt.axvline(1, color='g')\n",
    "plt.show()\n",
    "# It shows below that the area between two green lines is approx. 68%"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating this probability involved calculating the area under the pdf from the value of 1 and below. To put it another way, we need to integrate the pdf. We could just add together the known areas of chunks (from -Inf to 0 and then 0 to $+\\sigma$ in the example above. One way to do this is using look up tables (literally). Fortunately, scipy has this functionality built in with the cdf() function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ Use the cdf() function to answer the question above again and verify you get the same answer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6826894921370859\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# This approach takes the cdf's and subtracts them\n",
    "print(norm.cdf(1)-norm.cdf(-1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ Using our knowledge of the population parameters for our townsfolk's heights, what is the probability of selecting one person at random and their height being 177 cm or less? Calculate this using both of the approaches given above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8413447460685429\n"
     ]
    }
   ],
   "source": [
    "prob_ht_177 = norm.cdf(177, 172, 5) # Note the 5 value the docs call it Scale the exercise calls it the Std Deviation\n",
    "print(prob_ht_177)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ Turning this question around. Let's say we randomly pick one person and measure their height and find they are 2.00 m tall? How surprised should we be at this result, given what we know about the population distribution? In other words, how likely would it be to obtain a value at least as extreme as this? Express this as a probability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1138987855742525e-253\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# Not sure about this one it seems low\n",
    "prob_ht_2 = norm.cdf(2, 172, 5)\n",
    "print(prob_ht_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could calculate this probability by virtue of knowing the population parameters. We were then able to use the known properties of the relevant normal distribution to calculate the probability of observing a value at least as extreme as our test value. We have essentially just performed a z-test (albeit without having prespecified a threshold for our \"level of surprise\")!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're about to come to a pinch, though here. We've said a couple of times that we rarely, if ever, know the true population parameters; we have to estimate them from our sample and we cannot even begin to estimate the standard deviation from a single observation. This is very true and usually we have sample sizes larger than one. This means we can calculate the mean of the sample as our best estimate of the population mean and the standard deviation as our best estimate of the population standard deviation. In other words, we are now coming to deal with the sampling distributions we mentioned above as we are generally concerned with the properties of the sample means we obtain.\n",
    "\n",
    "Above, we highlighted one result from the CLT, whereby the sampling distribution (of the mean) becomes narrower and narrower with the square root of the sample size. We remind ourselves that another result from the CLT is that _even if the underlying population distribution is not normal, the sampling distribution will tend to become normal with sufficiently large sample size_. This is the key driver for us 'requiring' a certain sample size, for example you may frequently see a minimum sample size of 30 stated in many places. In reality this is simply a rule of thumb; if the underlying distribution is approximately normal then your sampling distribution will already be pretty normal, but if the underlying distribution is heavily skewed then you'd want to increase your sample size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ Let's now start from the position of knowing nothing about the heights of people in our town.\n",
    "* Use our favorite random seed of 47, to randomly sample the heights of 50 townsfolk\n",
    "* Estimate the population mean using np.mean\n",
    "* Estimate the population standard deviation using np.std (remember which denominator to use!)\n",
    "* Calculate the (95%) [margin of error](https://www.statisticshowto.datasciencecentral.com/probability-and-statistics/hypothesis-testing/margin-of-error/#WhatMofE) (use the exact critial z value to 2 decimal places - [look this up](https://www.statisticshowto.datasciencecentral.com/probability-and-statistics/find-critical-values/) or use norm.ppf())\n",
    "* Calculate the 95% Confidence Interval of the mean\n",
    "* Does this interval include the true population mean?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to Calculate Margin of Error: Steps\n",
    "\n",
    "## Step 1: Find the critical value. The critical value is either a t-score or a z-score. If you aren’t sure, see: T-score vs z-score. In general, for small sample sizes (under 30) or when you don’t know the population standard deviation, use a t-score. Otherwise, use a z-score.\n",
    "Click here for a minute video that shows you how to find a critical value.\n",
    "\n",
    "The basic z score formula for a sample is:\n",
    "z = (x – μ) / σ\n",
    "For example, let’s say you have a test score of 190. The test has a mean (μ) of 150 and a standard deviation (σ) of 25. Assuming a normal distribution, your z score would be: \n",
    "z = (x – μ) / σ \n",
    "= 190 – 150 / 25 = 1.6.\n",
    "\n",
    "## Step 2: Find the Standard Deviation or the Standard Error. These are essentially the same thing, only you must know your population parameters in order to calculate standard deviation. Otherwise, calculate the standard error (see: What is the Standard Error?).\n",
    "Click here for a short video on how to calculate the standard error.\n",
    "\n",
    "## Step 3: Multiply the critical value from Step 1 by the standard deviation or standard error from Step 2. For example, if your CV is 1.95 and your SE is 0.019, then:\n",
    "1.95 * 0.019 = 0.03705\n",
    "Example question: 900 students were surveyed and had an average GPA of 2.7 with a standard deviation of 0.4. Calculate the margin of error for a 90% confidence level:\n",
    "The critical value is 1.645 (see this video for the calculation)\n",
    "The standard deviation is 0.4 (from the question), but as this is a sample, we need the standard error for the mean. The formula for the SE of the mean is standard deviation / √(sample size), so: 0.4 / √(900) = 0.013.\n",
    "1.645 * 0.013 = 0.021385\n",
    "\n",
    "## Example question: Find a critical value for a 90% confidence level (Two-Tailed Test). \n",
    "### Step 1: Subtract the confidence level from 100% to find the α level: 100% – 90% = 10%.\n",
    "### Step 2: Convert Step 1 to a decimal: 10% = 0.10.\n",
    "### Step 3: Divide Step 2 by 2 (this is called “α/2”).\n",
    "0.10 = 0.05. This is the area in each tail.\n",
    "### Step 4: Subtract Step 3 from 1 (because we want the area in the middle, not the area in the tail):\n",
    "1 – 0.05 = .95.\n",
    "\n",
    "Step 5: Look up the area from Step in the z-table. The area is at z=1.645. This is your critical value for a confidence level of 90%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "172.7815108576788"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seed(47)\n",
    "# take your sample now\n",
    "samples50 = townsfolk_sampler(50)\n",
    "mean50 = np.mean ( samples50 )\n",
    "mean50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.195424364433547"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "std50 = np.std ( samples50 , ddof=1 ) \n",
    "std50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1628909342876539\n"
     ]
    }
   ],
   "source": [
    "# margin of error = ( CV * std ) / sqrt ( sample size )\n",
    "# ppf uses mean=0 and stddev=1, which is the \"standard\" normal distribution.\n",
    "## cv = ( 1 - (( 100 - 95 )/100) ) / 2\n",
    "margin_of_error = norm.ppf(0.975)*std50 / np.sqrt(50)\n",
    "print(margin_of_error)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[171.61861992339115, 173.94440179196644]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "mean_conf_interval = [mean50 - margin_of_error, mean50 + margin_of_error]\n",
    "print(mean_conf_interval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ Above we calculated the confidence interval using the critical z value. What is the problem with this? What requirement, or requirements, are we (strictly) failing?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The problem with this approach is that in the real world we rarely (if ever) actually \n",
    "# know your population parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Q:__ Calculate the 95% confidence interval for the mean using the _t_ distribution. Is this wider or narrower than that based on the normal distribution above? If you're unsure, you may find this [resource](https://www.statisticshowto.datasciencecentral.com/probability-and-statistics/confidence-interval/) useful. For calculating the critical value, remember how you could calculate this for the normal distribution using norm.ppf()."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__A:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1829831592442817\n"
     ]
    }
   ],
   "source": [
    "\n",
    "margin_of_error = t.ppf(0.975, mean50)*std50/np.sqrt(50-1)\n",
    "print(margin_of_error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[171.5985276984345, 173.96449401692308]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "mean_conf_interval = [mean50 - margin_of_error, mean50 + margin_of_error]\n",
    "print(mean_conf_interval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "# It is wider than the previous confidence interval\n",
    "# This is because we are estimating population parameters from a sample which has greater uncertainty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is slightly wider than the previous confidence interval. This reflects the greater uncertainty given that we are estimating population parameters from a sample."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning outcomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having completed this project notebook, you now have hands-on experience:\n",
    "* sampling and calculating probabilities from a normal distribution\n",
    "* the correct way to estimate the standard deviation of a population (the population parameter) from a sample\n",
    "* what a sampling distribution is and how the Central Limit Theorem applies\n",
    "* how to calculate critical values and confidence intervals"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}