Inferential_Statistics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "MathJax.Hub.Config({\n",
       "    TeX: { equationNumbers: { autoNumber: \"AMS\" } }\n",
       "});\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "MathJax.Hub.Config({\n",
    "    TeX: { equationNumbers: { autoNumber: \"AMS\" } }\n",
    "});"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.stats as smstats\n",
    "\n",
    "import matplotlib\n",
    "#matplotlib.use('nbagg')\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inferential Statistics\n",
    "\n",
    "Inferential statistics involves drawing conclusions about a population from a sample. \n",
    "\n",
    "A sample **statistic** is a quantity calculated from the sample values. Typical notation is $\\bar{x}$ for the sample mean, $s$ for the sample standard deviation. The goal of inferential statistics is to estimate values for parameters of the population, such as its mean $\\mu$ and its standard deviation $\\sigma$. \n",
    "\n",
    "Suppose we draw a sample of size $n$. The particular set of values obtained are denoted $x_1, x_2 \\cdots x_n$. These are particular realizations of the random variables $X_1, X_2 \\cdots X_n$. We have a **random sample** if they are independent of each other and identically distributed (iid). \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Central Limit theorem\n",
    "If $X_1, X_2 \\cdots X_n$ are independent, identically distributed random variables with finite mean and variance $\\mu$ and $\\sigma^2$, and $S_n = \\sum_{i=1}^n X_i$, then the quantity $\\frac{S_n-n\\mu}{\\sigma\\sqrt{n}}$ has a normal distribution asymptotically.\n",
    "$$\n",
    "\\lim_{n \\to \\infty} \\frac{S_n-n\\mu}{\\sigma\\sqrt{n}}\\; \\sim \\; N(0,1)\n",
    "\\implies \\lim_{n \\to \\infty} P\\left(a \\le \\frac{S_n-n\\mu}{\\sigma\\sqrt{n}} \\le b \\right) = \\frac{1}{\\sqrt{2\\pi}} \\int_a^b e^{-\\frac{u^2}{2}} du\n",
    "$$\n",
    "\n",
    "**Proof**\n",
    "\n",
    "Let $Z_i = \\frac{X_i-\\mu}{\\sigma\\sqrt(n)},\\quad \\sum Z_i = \\frac{S_n-n\\mu}{\\sigma\\sqrt{n}}$ and note first that\n",
    "$$\n",
    "M_{\\sum Z_i}(t) = E(e^{t\\sum Z_i}) = E(e^{tZ_1}e^{tZ_2} \\cdots e^{tZ_n}) = E(e^{tZ_1})E(e^{tZ_2}) \\cdots E(e^{tZ_n}) = [ M_{Z_1}(t) ]^n\n",
    "$$\n",
    "because the $Z_i$ are independent. Now,\n",
    "$$\n",
    "\\begin{split}\n",
    "M_{Z_1}(t) &= E(e^{tZ_1}) = E\\left(1 + tZ_1  +\\frac{t^2Z_1^2}{2!} + \\frac{t^3Z_1^3}{3!}+\\cdots\\right)\\\\\n",
    "&= 1 + tE(Z_1) + \\frac{t^2E(Z_1^2)}{2!} + \\frac{t^3E(Z_1^3)}{3!}+\\cdots\n",
    "\\end{split}\n",
    "$$\n",
    "But $E(Z_1) = 0, \\; E(Z_1^2) = \\frac{1}{n}, \\; E(Z_1^3)=\\frac{1}{n^{\\frac{3}{2}}}E\\left(\\frac{(X_1-\\mu)^3}{\\sigma^3}\\right)$ and so on. This means that\n",
    "$$\n",
    "\\lim_{n \\to \\infty}[ M_{Z_1}(t) ]^n =  \\lim_{n \\to \\infty} \\left(1+\\frac{1}{n} \\frac{t^2}{2!} + \\frac{1}{n^{\\frac{3}{2}}}E\\left(\\frac{(X_1-\\mu)^3}{\\sigma^3}\\right) + \\cdots \\right)^n\n",
    "$$\n",
    "\n",
    "Here, we make the crucial observation that terms with powers of $n$ in the denominator greater than 1 in this expression go to zero as $n \\to \\infty$. So,\n",
    "$$\n",
    "\\lim_{n \\to \\infty} [ M_{Z_1}(t) ]^n = \\lim_{n \\to \\infty} \\left(1+\\frac{1}{n} \\frac{t^2}{2} \\right)^n = e^{\\frac{t^2}{2}}\n",
    "$$\n",
    "which is the moment generating function of a standard normal distribution. From the uniqueness of moment generating functions, this implies that \n",
    "$$\n",
    "\\frac{S_n-n\\mu}{\\sigma\\sqrt{n}} \\text{  is distributed normally as  } n \\to \\infty\n",
    "$$\n",
    "$$\\tag*{$\\blacksquare$}$$\n",
    "\n",
    "This is a powerful theorem because it says that if we have enough variables from an arbitrary distribution, their sum will be distributed normally."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimators\n",
    "A function $g(X_1, X_2 \\cdots X_n)$ designed to estimate a population parameter is called an **estimator**. For example, the function\n",
    "$$\n",
    "g(X_1, X_2 \\cdots X_n) = \\bar{X} = \\frac{1}{n}\\sum_{i=1}^n X_i\n",
    "$$\n",
    "is an estimator for $\\mu$. Note that we could choose some other function as an estimator: a trimmed mean or a median are examples. The value computed from a particular sample is the **estimate**, and is denoted by a hat over the parameter. For example, an estimate of $\\mu$ is $\\hat{\\mu}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bias and Variance of estimators\n",
    "Estimator functions are designed to estimate population parameters. Given a choice among estimators, we can develop criteria to select the best one among them. The main criteria are bias and efficiency. For an estimator function $g$ designed to estimate the population parameter $\\theta$,\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "\\text{Bias} &\\equiv E(g) - \\theta \\\\\n",
    "\\text{var}(g) &\\equiv E\\left[(g-E(g))^2\\right]\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "where the expectations are taken over all samples from the population. Clearly, we prefer estimators with low or zero bias and small variance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Mean Squared Error\n",
    "\n",
    "A metric that incorporates both estimator bias and variance is the mean squared error (MSE). For a given estimator $g(X_1,X_2,\\cdots X_n)$ and population parameter $\\theta$,\n",
    "\n",
    "$$\n",
    "MSE \\equiv E\\left[(g-\\theta)^2\\right] = E\\left[g^2 - 2 \\theta g + \\theta^2\\right] = E(g^2) -2\\theta E(g) + \\theta^2 \n",
    "$$\n",
    "\n",
    "Now, for any random variable $Y, E(Y^2) = \\text{var}(Y) + E(Y)^2$. So,\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "MSE &= E(g^2) -2\\theta E(g) + \\theta^2\\\\ \n",
    "&= [E(g)^2 + \\text{var}(g)]-2\\theta E(g)] + \\theta^2\\\\\n",
    "&= [E(g) - \\theta]^2 + \\text{var}(g)\n",
    "\\end{split}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimator for population mean $\\mu$\n",
    "\n",
    "Each of the random variables $X_i$ has an expected value $\\mu$ and variance $\\sigma^2$. We could pick one of them, say $X_1$, as our estimator for the mean. We will then have an unbiased estimator with variance $\\sigma^2$. However, it is possible to pick a better estimator, the sample mean $\\bar{X} = \\frac{1}{n}\\sum X_i$. The sample mean is unbiased, but more importantly, it is more **efficient** because it has a lower standard deviation:\n",
    "\n",
    "\\begin{equation}\n",
    "E(\\bar{X}) = E\\left[\\frac{1}{n}\\sum_{i=1}^n X_i\\right] = \\frac{1}{n}\\sum_{i=1}^n E(X_i) = \\frac{1}{n}\\; n\\mu = \\mu \\end{equation}\n",
    "\n",
    "\\begin{align}\n",
    "\\text{var}(\\bar{X}) &= \\frac{1}{n^2}\\text{var}\\sum X_i = \\frac{1}{n^2} \\sum\\text{var}( X_i) +\\sum_{i<j} 2 \\text{cov}(X_i,X_j) \\nonumber \\\\\n",
    " &= \\frac{1}{n^2} \\sum\\text{var} (X_i) = \\frac{1}{n^2}\\; n\\sigma^2 = \\frac{\\sigma^2}{n}\n",
    "\\end{align}\n",
    "\n",
    "Here, we have used the fact that $\\text{cov}(X_i,X_j) = 0\\; \\forall\\; i,j$ since the $\\{X_i\\}$ are all independent of each other.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimator for population variance $\\sigma^2$\n",
    "Let $g(X_1,X_2,\\cdots X_n) = \\sum (X_i - \\bar{X})^2$. Then\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "g &= \\sum (X_i^2 -2\\bar{X}X_i + \\bar{X}^2) = \\sum X_i^2 -2n\\bar{X}^2 + n \\bar{X}^2 = \\sum X_i^2 - n \\bar{X}^2\\\\\n",
    "  \\\\\n",
    "E(g) &= \\sum E(X_i^2) - n E(\\bar{X}^2)\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "We now use the fact that for any random variable $Y, E(Y^2) = \\text{var}(Y) + E(Y)^2$. So,\n",
    "$$\n",
    "\\begin{split}\n",
    "E(g)&= \\sum(\\sigma^2 + \\mu^2) - n\\left(\\frac{\\sigma^2}{n}+\\mu^2\\right)\\\\\n",
    "&= n(\\sigma^2 + \\mu^2) - (\\sigma^2 + n\\mu^2)\\\\\n",
    "&= (n-1)\\sigma^2\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "This implies that the sum of squared deviations from the mean is a biased estimator of the population variance. However, this is easily fixed. The **sample variance**, $s^2$, defined as\n",
    "\\begin{equation}\n",
    "s^2 \\equiv \\frac{1}{n-1} \\sum_{i=1}^n (X_i-\\bar{X})^2\n",
    "\\end{equation}\n",
    "gives us an unbiased estimator of the population variance $\\sigma^2$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true,
    "toc-nb-collapsed": true
   },
   "source": [
    "## Sampling distributions\n",
    "\n",
    "As seen above, estimators are random variables whose values vary from sample to sample. The distributions corresponding to estimators are known as **sampling distributions**. Knowing the sampling distribution of an estimator allows us to draw conclusions about population parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true,
    "toc-nb-collapsed": true
   },
   "source": [
    "### Distribution of the sample mean\n",
    "First, consider the case when the **population is normal**, which implies that $X_i \\sim N(\\mu,\\sigma^2)$ for all $i$. We proceed by considering moment generating functions. For $X \\sim N(\\mu,\\sigma^2)$, the moment generating function $M_{X}(t) \\equiv E\\left(e^{tX}\\right) = e^{\\mu t + \\sigma^2\\frac{t^2}{2}}$. This function in unique.\n",
    "\n",
    "The moment generating function for $\\bar{X}$ is given by\n",
    "\n",
    "$$\n",
    "M_{\\bar{X}}(t) = E\\left(e^{t\\sum\\frac{X_i}{n}}\\right)\n",
    "= E\\left(e^{t\\frac{X_1}{n}}\\dot  e^{t\\frac{X_2}{n}}\\cdots e^{t\\frac{X_n}{n}} \\right)\n",
    "= E\\left(e^{t\\frac{X_1}{n}}\\right)E\\left(e^{t\\frac{X_2}{n}}\\right) \\cdots E\\left(e^{t\\frac{X_n}{n}}\\right)\n",
    "$$\n",
    "\n",
    "where the last equality follows from the fact that the ${X_i}$ are independent of each other.\n",
    "$$\n",
    "\\begin{split}\n",
    "M_{\\bar{X}}(t) &= \\Pi_{i=1}^n M_X(t/n)\\\\\n",
    "&= \\Pi_{i=1}^n e^{\\mu \\frac{t}{n} + \\frac{\\sigma^2}{n^2}\\frac{t^2}{2}}\\\\\n",
    "&= e^{\\sum_{i=1}^n \\mu \\frac{t}{n} + \\frac{\\sigma^2}{n^2}\\frac{t^2}{2}}\\\\\n",
    "&= e^{\\mu t + \\frac{\\sigma^2}{n}\\frac{t^2}{2}}\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "This implies that\n",
    "$$\n",
    "\\boxed{\n",
    "\\bar{X} \\sim N(\\mu,\\frac{\\sigma^2}{n}), \\;\\text{population normal}\n",
    "}\n",
    "$$\n",
    "\n",
    "In case the **population is not normal**, we use the central limit theorem, from which we get\n",
    "$$\n",
    "\\lim_{n\\to\\infty} P\\left(a \\le \\frac{\\bar{X}-\\mu}{\\sigma/\\sqrt{n}} \\le b\\right) = \\frac{1}{\\sqrt{2\\pi}} \\int_a^b e^{-\\frac{u^2}{2}du}\n",
    "$$\n",
    "In other words,\n",
    "$$\n",
    "\\boxed{\n",
    "Z=\\frac{\\bar{X}-\\mu}{\\sigma/\\sqrt{n}} \\;\\text{ is normally distributed as } n \\to \\infty \\text{ for a non-normal population}\n",
    "}\n",
    "$$\n",
    "\n",
    "In many cases, the population variance is unknown. In that case we are forced to use a sample estimate of the population variance. To this end, note that\n",
    "$$\n",
    "\\frac{\\bar{X}-\\mu}{s/\\sqrt{n}} = \\frac{\\bar{X}-\\mu}{\\frac{s}{\\sigma}\\;\\frac{\\sigma}{\\sqrt{n}}} = \\frac{(\\bar{X}-\\mu) / (\\sigma/\\sqrt{n})}{\\sqrt{s^2/\\sigma^2}}\n",
    "$$\n",
    "\n",
    "The numerator of the last quantity is a standard normal variable. As shown below, the denominator is a $\\chi^2_{n-1}$ variable divided by its degrees of freedom, $(n-1)$. This implies that the whole quantity has a t distribution with $(n-1)$ degrees of freedom\n",
    "\n",
    "$$\n",
    "\\boxed{\n",
    "\\frac{\\bar{X}-\\mu}{s/\\sqrt{n}}\\; \\text{ is t distributed with } (n-1) \\text{ degrees of freedom, population normal}\n",
    "}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distribution of the sample variance\n",
    "Assume that the population is normal, with finite mean $\\mu$ and variance $\\sigma^2$. Now,\n",
    "$$\n",
    "\\begin{split}\n",
    "\\sum (X_i-\\mu)^2 &= \\sum\\left( X_i-\\bar{X}+\\bar{X}-\\mu \\right)^2\\\\\n",
    "&= \\sum (X_i-\\bar{X})^2 + \\sum (\\bar{X}-\\mu)^2 + 2 \\sum (X_i-\\bar{X})(\\bar{X}-\\mu)\\\\\n",
    "&= \\sum (X_i-\\bar{X})^2 + n (\\bar{X}-\\mu)^2 + 2 (\\bar{X}-\\mu)\\sum(X_i-\\bar{X})\n",
    "\\end{split}\n",
    "$$\n",
    "Now, since $\\sum(X_i-\\bar{X}) = 0$, we have\n",
    "$$\n",
    "\\sum (X_i-\\mu)^2 = \\sum (X_i-\\bar{X})^2 + n (\\bar{X}-\\mu)^2\n",
    "$$\n",
    "Dividing throughout by $\\sigma^2$,\n",
    "$$\n",
    "\\sum \\left(\\frac{X_i-\\mu}{\\sigma}\\right)^2 = \\sum \\frac{(X_i-\\bar{X})^2}{\\sigma^2} + \\left(\\frac{\\bar{X}-\\mu}{\\sigma/\\sqrt{n}}\\right)^2\n",
    "$$\n",
    "We can rewrite this as\n",
    "$$\n",
    "\\sum Z_i^2 = \\frac{(n-1)s^2}{\\sigma^2} + Z_{\\bar{X}}^2\n",
    "$$\n",
    "The $Z_i$ are all standard normal variables, and so is $Z_{\\bar{X}}$, because we are assuming that the population is normal.\n",
    "\n",
    "Let $Q = \\sum Z_i^2, R = \\frac{(n-1)s^2}{\\sigma^2}, S = Z_{\\bar{X}}^2$, and note that\n",
    "$$\n",
    "Q \\sim \\chi_n^2, \\qquad S \\sim \\chi_1^2\n",
    "$$\n",
    "since the sum of $n$ squared standard normal variables is $\\chi_n^2$ distributed.\n",
    "\n",
    "Now, for independent random variables $R, S$ and $Q=R+S$, the moment generating functions can be written as:\n",
    "$$\n",
    "M_Q(t) = E(e^{tQ}) = E(e^{tR}e^{tS}) = E(e^{tR}) E(e^{tS}) = M_R(t) M_S(t) \\implies M_R(t) = \\frac{M_Q(t)}{M_S(t)}\n",
    "$$\n",
    "Since\n",
    "$$\n",
    "M_Q(t) = \\left(1-2t\\right)^{-\\frac{n}{2}} \\text{  and  } M_S(t) = \\left(1-2t\\right)^{-\\frac{1}{2}}\n",
    "$$\n",
    "we get $M_R(t) = \\left(1-2t\\right)^{-\\frac{n-1}{2}}$.\n",
    "\n",
    "Since the moment generating function is unique,\n",
    "$$\n",
    "\\boxed{\n",
    "\\frac{(n-1)s^2}{\\sigma^2} \\sim \\chi_{n-1}^2, \\; \\text{population normal}\n",
    "}\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distribution of sample proportion\n",
    "Suppose that the population consists of individuals with a characteristic that is a Bernoulii distributed random variable $X \\in \\{0,1\\}$ with\n",
    "\n",
    "$$ \n",
    "    P[X=1] = \\pi \\qquad P[X=0] = 1-\\pi \\qquad \\mu = \\pi \\qquad \\sigma^2 = \\pi(1-\\pi) \n",
    "$$\n",
    "    \n",
    "We draw a sample of size $n$, obtaining random variables $X_1, X_2 \\cdots X_n$, and define the sample proportion $p$\n",
    "$$\n",
    "p \\equiv \\bar{X} = \\frac{X_1+X_2+\\cdots+X_n}{n}\n",
    "$$\n",
    "\n",
    "Since the $X_i$ are iid,\n",
    "$$\n",
    "\\begin{split}\n",
    "E(p) &= \\pi \\\\\n",
    "\\sigma^2_p &= \\frac{\\sigma^2}{n} =\\frac{\\pi(1-\\pi)}{n}  \\\\\n",
    "\\end{split}\n",
    "$$\n",
    "Note that $p$ can only take on discrete values $k/n$, where $0\\le k\\le n$. Now, with $c=np$, the sample count,\n",
    "$$\n",
    "P [X_1+X_2 \\cdots X_n = c] = \\binom{n}{c}\\; \\pi^c\\; (1-\\pi)^{n-c}\n",
    "$$\n",
    "\n",
    "Using the Cenrtral Limit Theorem, we know that \n",
    "$$\n",
    "\\lim_{n\\to\\infty} P\\left(a \\le \\frac{\\bar{X}-\\pi}{\\sigma/\\sqrt{n}} \\le b\\right) = \\frac{1}{\\sqrt{2\\pi}} \\int_a^b e^{-\\frac{u^2}{2}du}\n",
    "$$\n",
    "\n",
    "or, in other words, the random variable\n",
    "\n",
    "$$\n",
    "Z=\\frac{p-\\pi}{\\sqrt{\\frac{\\pi(1-\\pi)}{n}}}\n",
    "$$\n",
    "\n",
    "is asymptotically normal. This is actually a statement that the Binomial distribution is asymptotically normal. Often, textbooks suggest a rule of thumb that this approximation can be used if the conditions $n\\pi \\ge 10$ and $n(1-\\pi) \\ge 10$ are met."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interval estimation\n",
    "\n",
    "Sample statistics are point estimates of population parameters, since they give us one number as an estimate. If we know the sampling distribution of the statistic, we can calculate a range within which the population parameter falls with a pre-specified probability. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Population mean\n",
    "Suppose that we have a sample mean and we seek a range for the population mean with 95% probability. Since the population is normal, we know that the sampling distribution of the mean will be normal. This implies that\n",
    "$$\n",
    "P \\left(-z_{0.025} \\le \\frac{\\bar{X}-\\mu}{\\sigma\\sqrt{n}} \\le z_{0.025} \\right) = 0.95\n",
    "$$\n",
    "or, in general,\n",
    "$$\n",
    "P \\left(-z_{\\alpha/2} \\le \\frac{\\bar{X}-\\mu}{\\sigma\\sqrt{n}} \\le z_{\\alpha/2} \\right) = (1-\\alpha)\n",
    "$$\n",
    "Here $\\alpha$ is known as the **level of significance**, $(1-\\alpha)$ is called the **confidence level**, and $z_{\\alpha/2}$ is the point beyond which the area under the standard normal curve is $\\alpha/2$.\n",
    "\n",
    "The **confidence interval** for the population mean results from inverting the inequalities\n",
    "\n",
    "$$\n",
    "\\bar{X}-z_{\\alpha/2}\\frac{\\sigma}{\\sqrt{n}} \\le \\mu \\le \\bar{X}+z_{\\alpha/2}\\frac{\\sigma}{\\sqrt{n}}\\quad \\text{with confidence level } (1-\\alpha), \\text{ for a normal population}\n",
    "$$\n",
    "\n",
    "As is generally the case, if the population variance is unknown, we use the sample variance and critical values for the t distribution with $(n-1)$ degrees of freedom.\n",
    "\n",
    "$$\n",
    "\\boxed{\n",
    "\\bar{X}-t_c\\frac{s}{\\sqrt{n}} \\le \\mu \\le \\bar{X}+t_c\\frac{s}{\\sqrt{n}}\\quad \\text{with confidence level } (1-\\alpha), \\text{ for a normal population}\n",
    "}\n",
    "$$\n",
    "\n",
    "The critical value $t_c$ is given by\n",
    "$$\n",
    "P(T_{n-1} \\le t_c) = 1-\\frac{\\alpha}{2}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4948364398216353"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "2.2621571627409915"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "(-0.11029681585827622, 1.0999696955015468)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "(-0.11029681585827622, 1.0999696955015468)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Interval estimation of population mean assuming normally distributed population\n",
    "\n",
    "np.random.seed(12345)\n",
    "n=10\n",
    "v1 = stats.norm.rvs(size=n)  # Standard normal random variates\n",
    "\n",
    "xbar = np.mean(v1)\n",
    "s = np.std(v1, ddof=1)    # ddof=1 uses (n-1) for the denominator, so gives us sample sd\n",
    "xbar\n",
    "\n",
    "ConfLevel = 0.95\n",
    "alpha = 1 - ConfLevel     # Level of significance\n",
    "\n",
    "tc = stats.t.ppf(1-alpha/2,n-1) # critical value; P(-\\infty,tc) = prob in right tail = 1-alpha/2\n",
    "tc\n",
    "xbar-tc*s/np.sqrt(n),xbar+tc*s/np.sqrt(n)\n",
    "\n",
    "smstats.weightstats._tconfint_generic(xbar,s/np.sqrt(n),n-1,alpha,alternative=\"2-sided\") # Single function call"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Population variance\n",
    "\n",
    "For a normal population, $\\frac{(n-1)s^2}{\\sigma^2} \\sim \\chi^2_{n-1}$. The $\\chi^2$ density is not symmetric and so the left and right critical values for a given confidence level are arbitrary, to some extent. One way to pick specific critical values is to search for the values that lead to the smallest interval, but the convention is just to pick the values that lead to the same remaining probability on either side of the range. Thus, with\n",
    "$$\n",
    "P\\left(c_1 \\le \\frac{(n-1)s^2}{\\sigma^2} \\le c_2 \\right) = 1-\\alpha\n",
    "$$\n",
    "$c_1$ and $c_2$ are chosen such that\n",
    "$$\n",
    "P(\\chi^2_{n-1} \\le c_1) = \\frac{\\alpha}{2} \\quad P(\\chi^2_{n-1} \\ge c_2) = \\frac{\\alpha}{2}\n",
    "$$\n",
    "This yields\n",
    "$$\n",
    "\\boxed{\n",
    "\\frac{(n-1) s^2}{c_2} \\le \\sigma^2 \\le \\frac{(n-1) s^2}{c_1}\\quad \\text{with confidence level } (1-\\alpha), \\text{ for a normal population}\n",
    "}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample variance:  0.8459179615720821\n",
      "Critical values:  2.7003894999803584 19.02276779864163\n",
      "Probabilities up to critical values:  0.025000000000000015 0.975\n",
      "95.0 % confidence interval for the variance:  0.3385519314309392 2.3849132798950827\n"
     ]
    }
   ],
   "source": [
    "# Interval estimation of population variance assuming normally distributed population\n",
    "\n",
    "np.random.seed(12345)\n",
    "n=10\n",
    "v1 = stats.norm.rvs(size=n)  # Standard normal random variates\n",
    "s = np.std(v1, ddof=1)    # ddof=1 uses (n-1) for the denominator, so gives us sample sd\n",
    "print(\"Sample variance: \",s)\n",
    "\n",
    "ConfLevel = 0.95\n",
    "alpha = 1 - ConfLevel     # Level of significance\n",
    "\n",
    "c1 = stats.chi2.ppf(alpha/2,n-1)    # P(0,c1) = prob in left tail = alpha/2\n",
    "c2 = stats.chi2.ppf(1-alpha/2,n-1)  # P(0,c2) = 1-prob in right tail = 1-alpha/2\n",
    "\n",
    "print(\"Critical values: \", c1,c2)\n",
    "print(\"Probabilities up to critical values: \",stats.chi2.cdf(c1,n-1), stats.chi2.cdf(c2,n-1))\n",
    "\n",
    "# Interval estimate\n",
    "print(100*ConfLevel,\"% confidence interval for the variance: \", (n-1)*(s**2)/c2,(n-1)*(s**2)/c1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Population proportion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the sampling distribution of the sample proportion is known, exact calculations are possible. We can derive implicit expressions for the bounds of the confidence interval $[p_l,p_u]$. For a given confidence level $\\alpha$\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "P[k \\le c | \\pi=p_u] &= \\sum_{k=0}^c \\binom{n}{k} p_u^k (1-p_u)^{n-k} = \\frac{\\alpha}{2} \\\\\n",
    "P[k \\ge c | \\pi=p_l] &= \\sum_{k=c}^n \\binom{n}{k} p_l^k (1-p_l)^{n-k} = \\frac{\\alpha}{2}\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "These expressions can be solved for $p_l$ and $p_u$. A useful relationship between the expressions and incomplete beta functions makes this easier. If $X$ is a binomially distributed random variable, then the probability of at least $k$ successes\n",
    "$$\n",
    "P(X \\ge k ; p) = \\sum_{r=k}^n \\binom{n}{r} p^r (1-p)^{n-r}\n",
    "$$\n",
    "can be treated as a continuous function of $p$ over the interval $[0,1]$. Differentiating this with respect to $p$ gives\n",
    "$$\n",
    "\\begin{split}\n",
    "\\frac{dP(X \\ge k ; p)}{dp} &= \\sum_{r=k}^n r \\binom{n}{r} p^{r-1} (1-p)^{n-r} - \\sum_{r=k}^{n-1} (n-r) \\binom{n}{r} p^r (1-p)^{n-r-1}\\\\\n",
    "&= \\sum_{r=k}^n n\\binom{n-1}{r-1} p^{r-1} (1-p)^{n-r} - \\sum_{r=k}^{n-1} n \\binom{n-1}{r} p^r (1-p)^{n-r-1}\\\\\n",
    "&= n \\left[\\sum_{s=k-1}^{n-1} \\binom{n-1}{s} p^{s} (1-p)^{n-s-1} - \\sum_{r=k}^{n-1} \\binom{n-1}{r} p^r (1-p)^{n-r-1}\\right]\\\\\n",
    "&= n \\binom{n-1}{k-1} p^{k-1}(1-p)^{n-k}= k \\binom{n}{k} p^{k-1}(1-p)^{n-k}\\\\\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "Now, for a beta distribution, the probability density function is\n",
    "$$\n",
    "f(x;\\; \\alpha,\\beta) = \\frac{\\Gamma (\\alpha+\\beta)}{\\Gamma(\\alpha)\\Gamma(\\beta)} x^{\\alpha-1}(1-x)^{\\beta-1}\n",
    "$$\n",
    "With $\\alpha=k,\\; \\beta=n-k+1,\\; x=p$ this becomes\n",
    "$$\n",
    "\\begin{split}\n",
    "f(p;\\;k,n-k+1) &= \\frac{\\Gamma (n+1)}{\\Gamma(k)\\Gamma(n-k+1)} p^{k-1}(1-p)^{n-k}\\\\\n",
    "&= k \\frac{n!}{k! (n-k)!} p^{k-1}(1-p)^{n-k} = k \\binom{n}{k}p^{k-1}(1-p)^{n-k}\\\\\n",
    "&= \\frac{dP(X \\ge k ; p)}{dp}\n",
    "\\end{split}\n",
    "$$\n",
    "Thus,\n",
    "$$\n",
    "P(X \\ge k ; p) = \\int_0^p f(p';\\;k,n-k+1)\\;dp' \\equiv I_\\beta(p;k,n-k+1)\n",
    "$$\n",
    "\n",
    "where $I_\\beta$ is the incomplete beta function. Thus, we now have to solve\n",
    "$$\n",
    "\\begin{split}\n",
    "I_\\beta (p_l;\\; c,n-c+1) &= \\frac{\\alpha}{2}\\\\\n",
    "I_\\beta (p_u;\\; c+1,n-c) &= 1-\\frac{\\alpha}{2}\\\\\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "to obtain $p_l$ and $p_u$. This can be done with the inverse beta or percent point functions for the beta distribution.\n",
    "\n",
    "This is known as the Clopper-Pearson confidence interval. The derivation above is reproduced from [here](http://faculty.washington.edu/fscholz/DATAFILES498B2008/ConfidenceBounds.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Normality assumption interval:  0.015974234910674623 0.5840257650893254\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.015974234910674623, 0.5840257650893254)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact interval:  0.0667395111777345 0.6524528500599972\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.0667395111777345, 0.6524528500599972)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Confidence intervals for population proportions\n",
    "n = 10\n",
    "p = 0.3\n",
    "c = int(n*p)        # Count of successes\n",
    "\n",
    "ConfLevel = 0.95\n",
    "alpha = 1- ConfLevel   \n",
    "\n",
    "# Assuming normality (terrible assumption for this set of sample data)\n",
    "\n",
    "s = np.sqrt(p*(1-p)/n)\n",
    "zc = stats.norm.ppf(1-alpha/2)\n",
    "pl = p-zc*s\n",
    "pu = p+zc*s\n",
    "print(\"Normality assumption interval: \",pl,pu)\n",
    "smstats.proportion.proportion_confint(c,n,alpha,method='normal')\n",
    "    \n",
    "#  Exact interval, incomplete beta function for pl: I(pl; c, n-c+1) = alpha/2\n",
    "#                               pu: I(pu; c+1,n-c) = 1-alpha/2\n",
    "        \n",
    "pl = stats.beta.ppf(alpha/2,c,n-c+1)\n",
    "pu = stats.beta.ppf(1-alpha/2,c+1,n-c)\n",
    "print(\"Exact interval: \",pl,pu)\n",
    "smstats.proportion.proportion_confint(c,n,alpha,method='beta')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true,
    "toc-nb-collapsed": true
   },
   "source": [
    "## Two samples: difference of means\n",
    "\n",
    "In many applications (such as drug trials), comparing the statistics from two samples is quite important. The samples could be drawn from two different populations or from the same population but subject to two different treatments.\n",
    "\n",
    "Consider two independent normal random variables, $Z_1 \\sim N(\\mu_1,\\sigma^2_1)$ and $Z_2\\sim N(\\mu_2,\\sigma^2_2)$. We are interested in the distribution of $Y = Z_1 - Z_2$. Now,\n",
    "$$\n",
    "\\begin{split}\n",
    "M_Y(t) &= E(e^{tY}) = E(e^{t(Z_1-Z_2)}) = E(e^{tZ_1}\\;e^{-tZ_2}) = E(e^{tZ_1})E(e^{-tZ_2})\\\\\n",
    "&=e^{\\mu_1 t + \\sigma_1^2 t^2/2}\\; e^{-\\mu_2 t + \\sigma_2^2 t^2/2}\\\\\n",
    "&= e^{(\\mu_1-\\mu_2) t + (\\sigma_1^2+\\sigma_2^2) t^2/2}\n",
    "\\end{split}\n",
    "$$\n",
    "which means that $Y \\sim N(\\mu_1-\\mu_2, \\sigma_1^2+\\sigma_2^2)$.\n",
    "\n",
    "Now consider the situation when we have two independent samples of sizes $n$ and $m$ respectively, from normally distributed populations. Given their sample means $\\bar{X}$ and $\\bar{Y}$, what can we say about the difference of the population means?\n",
    "\n",
    "Since the populations are normal, we know that $\\bar{X} \\sim N(\\mu_X,\\frac{\\sigma_X^2}{n})$ and $\\bar{Y} \\sim N(\\mu_Y,\\frac{\\sigma_Y^2}{m})$. So,\n",
    "\n",
    "$$\n",
    "\\bar{X}-\\bar{Y} \\sim N\\left(\\mu_X-\\mu_Y, \\frac{\\sigma_X^2}{n}+\\frac{\\sigma_Y^2}{m}\\right)\n",
    "$$\n",
    "\n",
    "We often do not know the population variances, and in analogy with the one sample case, we would like to be able to set up a statistic using sample variances that is t distributed. This can be done if we assume $\\sigma_X = \\sigma_Y = \\sigma$. Under this assumption\n",
    "$$\n",
    "Z = \\frac{\\bar{X}-\\bar{Y} -(\\mu_X-\\mu_Y)}{\\sigma \\sqrt{\\frac{1}{n}+\\frac{1}{m}}} \\sim N(0,1)\n",
    "$$\n",
    "Given that $T = Z/\\sqrt{V/k}$ is t distributed if $Z \\sim N(0,1)$ and $V \\sim \\chi^2_k$, we seek a $\\chi^2$ distributed estimator of $\\sigma$. For this, we note that\n",
    "$$\n",
    "\\frac{(n-1) s_X^2}{\\sigma^2} \\sim \\chi^2_{n-1}, \\quad \\frac{(m-1) s_Y^2}{\\sigma^2} \\sim \\chi^2_{m-1},\\quad \\frac{(n-1)s^2_X + (m-1)s^2_Y}{\\sigma^2} \\sim \\chi^2_{n+m-2}\n",
    "$$\n",
    "So we define a **pooled variance** estimator\n",
    "$$\n",
    "s^2_p \\equiv \\frac{(n-1)s^2_X + (m-1)s^2_Y}{n+m-2} \\implies V = \\frac{(n+m-2) s^2_p}{\\sigma^2} \\sim \\chi^2_{n+m-2}\n",
    "$$\n",
    "This means that\n",
    "$$\n",
    "\\frac{Z}{\\sqrt{V/(n+m-2)}} = \\frac{Z}{s_p/\\sigma} = \\frac{\\bar{X}-\\bar{Y} -(\\mu_X-\\mu_Y)}{s_p \\sqrt{\\frac{1}{n}+\\frac{1}{m}}} \\sim T_{n+m-2}\n",
    "$$\n",
    "Summarizing,\n",
    "$$\n",
    "\\boxed{\n",
    "s^2_p \\equiv \\frac{(n-1)s^2_X + (m-1)s^2_Y}{n+m-2}, \\qquad \\frac{\\bar{X}-\\bar{Y} -(\\mu_X-\\mu_Y)}{s_p \\sqrt{\\frac{1}{n}+\\frac{1}{m}}} \\sim T_{n+m-2}, \\quad \\text{populations normal, pooled variance}\n",
    "}\n",
    "$$\n",
    "\n",
    "If we cannot make the assumption that the population variances are equal, we have to adopt a different procedure. We know that\n",
    "$$\n",
    "Z = \\frac{\\bar{X}-\\bar{Y} -(\\mu_X-\\mu_Y)}{\\sqrt{\\frac{\\sigma^2_X}{n} + \\frac{\\sigma^2_Y}{n}}} = \\frac{\\bar{X}-\\bar{Y} -(\\mu_X-\\mu_Y)}{\\sigma_D} \\sim N(0,1)\n",
    "$$\n",
    "and we want to use\n",
    "$$\n",
    "T = \\frac{\\bar{X}-\\bar{Y} -(\\mu_X-\\mu_Y)}{\\sqrt{\\frac{s^2_X}{n} + \\frac{s^2_Y}{n}}} = \\frac{\\bar{X}-\\bar{Y} -(\\mu_X-\\mu_Y)}{s_D} = \\frac{Z}{s_D/\\sigma_D}\n",
    "$$\n",
    "This variable $T$, as we have defined it, would have a t distribution if $s_D/\\sigma_D$ were a $\\chi^2$ distributed variable divided by its degrees of freedom. Unfortunately, that is not the case. However, it is possible to specify a value for the degrees of freedom that provides an approximation. See the [Wikipedia entry](https://en.wikipedia.org/wiki/Welch%E2%80%93Satterthwaite_equation) for details and references. \n",
    "\n",
    "We want to find $d$ such that\n",
    "$$\n",
    "d\\frac{s_D^2}{\\sigma_D^2} \\stackrel{approximately}{\\sim} \\chi^2_d\n",
    "$$\n",
    "Assuming that the distributional equivalence holds, we have\n",
    "$$\n",
    "\\text{var}(s_D^2) = \\frac{\\sigma_D^4}{d^2} \\text{var}(\\chi^2_d) = \\frac{\\sigma_D^4}{d^2} 2d = \\frac{2\\sigma_D^4}{d}\n",
    "$$\n",
    "But from the definition of $s_D$, we have\n",
    "$$\n",
    "\\begin{split}\n",
    "\\text{var}(s_D^2) &= \\frac{\\text{var}(s_X^2)}{n^2} + \\frac{\\text{var}(s_Y^2)}{m^2}\\\\\n",
    "&= \\frac{\\sigma_X^4}{n^2 (n-1)^2} 2(n-1) + \\frac{\\sigma_Y^4}{m^2 (m-1)^2} 2(m-1)\\\\\n",
    "&=  \\frac{2\\sigma_X^4}{n^2 (n-1)} + \\frac{2\\sigma_Y^4}{m^2 (m-1)}\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "where the second line results from the fact that $(n-1)s_X^2/\\sigma_X^2 \\sim \\chi^2_{n-1}$ etc.\n",
    "\n",
    "Equating the two expressions for $\\text{var}(s_D^2)$, we have\n",
    "$$\n",
    " d = \\frac{\\left(\\frac{\\sigma_X^2}{n} + \\frac{\\sigma_Y^2}{m} \\right)^2}{\\frac{\\sigma_X^4}{n^2(n-1)} + \\frac{\\sigma_Y^4}{m^2(m-1)}}\n",
    "$$\n",
    "Summarizing, we have\n",
    "$$\n",
    "\\boxed{\n",
    "d = \\frac{\\left(\\frac{\\sigma_X^2}{n} + \\frac{\\sigma_Y^2}{m} \\right)^2}{\\frac{\\sigma_X^4}{n^2(n-1)} + \\frac{\\sigma_Y^4}{m^2(m-1)}}, \\quad T = \\frac{\\bar{X}-\\bar{Y} -(\\mu_X-\\mu_Y)}{\\sqrt{\\frac{s^2_X}{n} + \\frac{s^2_Y}{n}}} \\stackrel{approximately}{\\sim} T_d, \\quad \\text{populations normal, unequal variances}\n",
    "}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hypothesis testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we know the sampling distribution of a sample statistic, we can make hypotheses about population parameters, and evaluate the hypothesis in light of sample data. The criterion for accepting or rejecting a hypothesis is a decision rule. The criterion is probabilistic, in the sense that we can specify the probability with which we are willing to risk the rejection of a true hypothesis. This will be clearer with an example.\n",
    "\n",
    "Let a sample of size $n$ drawn from a normally distributed population with mean $\\mu$ have sample mean and standard deviation are $\\bar{x}$ and $s$ respectively. Now, let's say we make a hypothesis that the population mean is $\\mu_0$. Under this hypothesis, the **test statistic**\n",
    "$$\n",
    "T = \\frac{\\bar{x}-\\mu_0}{s/\\sqrt{n}}\n",
    "$$\n",
    "can be computed. We then decide whether this value of T is unusual or extreme. A highly unusual value would be reason to reject our hypothesis. We can decide what \"highly unusual\" means before we conduct the test. This is the **level of significance** $\\alpha$ of the test. So, for instance, a 5% level of significance would mean that we reject the hypothesis if values equal to or larger than $T$ occur less than 5% of the time.\n",
    "\n",
    "The standard hypothesis testing procedure involves both a null hypothesis and an alternative hypothesis. The null hypothesis $H_0$ is stated for a specific value of the population parameter, while the alternative hypothesis $H_a$ is for a range. Three cases arise:\n",
    "\n",
    "1. $H_0: \\mu = \\mu_0 \\quad H_a: \\mu \\ne \\mu_0 \\qquad \\text{Two-tailed test}$\n",
    "2. $H_0: \\mu = \\mu_0 \\quad H_a: \\mu \\lt \\mu_0 \\qquad \\text{Lower-tailed test}$\n",
    "3. $H_0: \\mu = \\mu_0 \\quad H_a: \\mu \\gt \\mu_0 \\qquad \\text{Upper-tailed test}$\n",
    "\n",
    "The combination $H_0$ and $H_a$ together determine the **rejection region**: if $T$ falls within this region, we reject the null hypothesis. For example, if we are using a two tailed test with $\\alpha=5\\%$, we would find critical values for the t distribution such that $2.5\\%$ of the area under the curve is captured by the left and right tails. If the test statistic $T$ falls within either of the two tails, we reject th enull hypothesis.\n",
    "\n",
    "In addition to or instead of specifying a significance level, we can compute the $p$ value, which is the probability of occurrence of values more extreme than the computed test statistic $T$. Again, this depends on the form of the alternative hypothesis.\n",
    "\n",
    "1. $H_0: \\mu = \\mu_0 \\quad H_a: \\mu \\ne \\mu_0 \\quad p = P[t \\le -|T|] + P[t \\ge |T|] \\qquad \\text{Two-tailed test}$\n",
    "2. $H_0: \\mu = \\mu_0 \\quad H_a: \\mu \\lt \\mu_0 \\quad p = P[t \\le -|T|]  \\qquad \\text{Lower-tailed test}$\n",
    "3. $H_0: \\mu = \\mu_0 \\quad H_a: \\mu \\gt \\mu_0 \\quad p = P[t \\ge |T|]  \\qquad \\text{Upper-tailed test}$\n",
    "\n",
    "Similar considerations apply to hypothesis tests for other population parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAADYCAYAAACJH8VwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3gc5dXw4d+j3qslW7JkC/cquQLGxoVmCARCSSgJhCQQkrwJKV9IJYGEvGmQN5BAQgKh10CC6d0Y925JlmXLRVa3rLrqZct8f8ysWAvJars7W859Xb4saXZnzuzOM3PmaaM0TUMIIYQQQgSGELMDEEIIIYQQ7iPJnRBCCCFEAJHkTgghhBAigEhyJ4QQQggRQCS5E0IIIYQIIJLcCSGEEEIEEEnufIhSql0pNcX4+Qml1G+8tN0NSqlbRvne08Y52D4ppc5VSpWc5n2TjPeGjiau06z3SqVUpbHuhe5ctzcppVYrparcuL6blVKb3bU+bxnp5+B6rCulvqiUes+NsRxQSq02fr5bKfWMG9f9M6XUo+5an69x9/F8mu1oSqlpo3ifUko9rpRqVkrt9ERsg2x3uVLqiHG++pxS6m2l1JcHeW2OsX9h3opP+C5J7jzAKIjOfw6lVJfL718c7H2apsVpmlbqzVg9bbB90jRtk6ZpM52/K6XKlFIXuCyvMN5rd3NI9wHfNta9z83r7uOvyVIw0TTtWU3TLhrqdcO90dI0ba6maRvGGtdAiY6mab/VNG1UN2Aj3HaSUupJpVSd8e/ufsvL+p3P3nNZdr5S6rhSqlYpdV2/de5VSsV7On4PWgFcCGRpmnamF7f7a+BB43y1TtO0SzRNe9KL2xd+SjJ8D9A0Lc75s1KqDLhF07QPzIhFKaUApWmaw4zt+6DJwIHRvFEpFeqBZNPrlFJhmqbZgm3bnhJg+/RnIAbIAdKBD5VS5ZqmPe7yms8Ocj67H/gsEAp8pJR6ySgvvwN+r2lam2dD96jJQJmmaR0mbHdU5ysR3KTmzouUUmcqpbYppSxKqRNKqQeVUhEuy4fVZGDUCm0x3t+ilDqklDrfZfkGpdT/KqW2AJ3AFKXUOUqpXcbrdymlzum32qlKqZ1KqVal1KtKqRSX9b1k3I23KKU2KqXm9nvvOKXU+0qpNqXUx0qpyUPtk2vthFLqaWAS8LpRG/Cj/k0MSqlEpdS/jM+tWin1G2U02SqlphnbbVFKNSilXhxge5FKqXb0C0+BUuqY8ffZxudlUXqz2uUu73lCKfV3pdRbSqkOYM0g30Wpse/Hld7UNxt4GFhm7I/FeO2lSql9xmdc6Vor4rK/X1ZKVRj78XOX5dFGPM1KqWJgab84fqKUOmbEUayUurJfjFuUUn9WSjUCdyulUpVSrxmx7ASm9t+3EcQWqZS6XylVY/y7XykV6fo9K6V+rJSqBR5XepPlS0qpZ4x49yulZiilfqr02qJKpdRFLuv/ilLqoPHaUqXUbYPFOkDsFxrlo0Up9SCg+n0um42flfH51BmfyX6l1Dyl1NeBLwI/Mr7L143Xlxn7VAh0KKXCVL/aZyBKKfWiEfdepVSey7ZPKRfGd/sbpVQs8DaQqT6pHctU/Zp5lVKXG8erxTh+Z7ssK1NK/VApVWjs94tKqahhfmSfBf6oaVqnpmllwL+Arw7zvbGaphVpmlYA9AKpSqkzgTM0Tfv3MNfh3IdMpdR/lFL1Rrm63eXvXerU89NC45gMN37/qnG8NCul3lUu56NhbPM1pVSTUuqoUupW4+9fAx7lk/L8q0Hef6vLcVqslFpk/H2oc8xDSqk3jfftUEpNNZYdA6bwyXkxUp3arSBUKXWfse+lwKX94jndOfNmpdRm4/3Nxmd8ict7U5TeDF1jLF/nsuwypVS+sT9blVK5w/l8hZdpmib/PPgPKAMuMH5eDJyNXmOaAxwEvufyWg2YZvz8BPCbQdZ5M2ADvg+EA9cCLUCKsXwDUAHMNbY1HmgGbjR+v974PdXl9dXAPCAW+A/wjMv2vgrEA5Hod+f5LsueANqAlcbyB4DNQ+0TsBqoGuhzMn7PMd4bZvz+CvAPI750YCdwm7HseeDn6DcrUcCK03wfrvGEA0eBnwERwHnGvsx0ibcFWO5cd791xQKtLq/PAOa6fEeb+71+NTDfWFcucBL4XL/9fQSIBvKAHmC2sfz3wCYgBcgGivp9fp8HMo11Xwt0ABn9jpfvGN9/NPAC8G9jH+YZ3//mQT6zoWL7NbDd+F7SgK3APS77bAP+gH58RAN3A93AWiOep4DjxncYDtwKHHfZ/qXoyacCVqHfsCwa6DjqF/c44/u8xljv941Ybun/HRmx7AGSjO3Mdvn8nqBfWUQ/XvON7yJ6gLJ+N2B12fYPjX0M738cDlU2XNb3jPHzDOP7vdBY94/Qj+MIlzh2oh8PKejnmW8M83zVAJzp8vvPgeZ++30SqAfeA/Jclm03jo08oMaIbRswYxjb7dtn9GN4D/BL9HI5BSgF1hrL1wO3urz3XuBh4+crjM9iNvqxdSewdaDyP0AMG4G/oZ9DFhj7eN5g5bnfez+PXoaWGsfPNPRat+GcYxqBM414nwVeOM15cQOfHL/fAA6hH4MpwEcM/5x5M/rxeSv6De83je9MGcvfBF4Eko19WGX8fSFQB5xlvO/LRoyRwzm+5J/3/pkeQKD/6184+y37HvCKy+8jSe76CqLxt53AjcbPG4Bfuyy7EdjZbx3bgJtdXv97l2Vz0O+8QwfYdpIRZ6JLnK4nozjADmSfbp8YQXKHnpz2YFxEjeXXAx8ZPz8F/BO9P8xQ34drPOcCtUCIy/Lngbtd4n3qNOuKBSzA1a6xuXxHg14MjNfcD/y53/5muSzfCVxn/FwKXOyy7OsMktQYy/OBK1xiqXBZFop+Yp/l8rffDhbvMGI7BnzGZdla9CYs5/fci0tijJ6ovO/y+2eBdufxhn4joQFJg8SzDvjuQMdRv9fdBGx3+V0BVQyc3J0HHEa/+Qrpt54nGDi5++pgZd3YR9dthwAngHP7H4dDlQ2X9TmTu18A/+637mpgtUscX3JZ/keM5GcY5eMZ4L/GdzDN+G57XJYvR0/QY4CfopefJGPZAvRzyQ7gfOB24B70G5l30ZOPVYNst2+f0ROHin7Lfwo8bvx8C7De5TutBFYav78NfK3fZ9MJTB7oc3d5XTb6eSve5W+/A54YTnk29u+7A/x9OOeYR12WfQY4NNAxZfy+gU+O3/W4JO3ARQz/nHkzcNRlWYzx3gnoN6kOIHmA/fk7xo2by99KBvte5Z95/6RZ1ouU3vT0htKbOFvRL6jjRrm6as0oWYZy9Dt1p0qXnzON5a7KgYmDvL4c/W5tnFH1/3ulN/m1op9s6Bd333s1TWsHmvrFMlbOO+ATRlOABf2ONN1Y/iP0k/xOo9ljuM1ImUCldmp/xNN9LqfQ9P4316LfQZ8wmlZmDfZ6pdRZSqmPjKamFuN9/b//WpefO9GT5b5Y+8Xpuu6bXJpKLOi1cQN+R+i1a2GnW98gTheb6/v7H4v1mqZ191vXSZefu4AG7ZP+jF3G/3EASqlLlFLbjeYyC/oFcDjl5pTPzCgvA36fmqatBx4EHgLqlFL/VEolDLH+QY+N/suNY6wK95SLUz5vY92VnHrcDvZdDeV29M//CPAqeiLSN7hD07QtmqZ1aXqz7e/Qb27ONZbla5q2WtO0s4Bi9Br/36I3af4K+ArwtFJKcXqT0ZulLS7H88/QExbQWxaWKaUy0FsMHOi12s73PuDyvib0c8NETi8TaNJO7RfY/1xwOtnoifBA6x3qHDPa7+p054ShzpmnbFfTtE7jxzj0fWnSNK15gG1OBv5fv+8mG/ee74UbSHLnXX9Hr0afrmlaAvoJa6gT3WAm9jtJTkKvzXNyTfxq0Aulq0nod/tO2f2WWdGbaG5Ab+q4AEhEr8WhX9x971VKxaE3EbjGMhzaaZZVot+FjtM0Lcn4l6Bp2lwATdNqNU27VdO0TOA24G9qeNMd1ADZSinXctD/czldXGia9q6maRei3+0eQm+6HOx9zwGvoddqJqL3yxvu93+CT39HABh9ih4Bvo3e1J6E3mzrum7XeOrRmycHXN8o9D++TncsjojS++79B32U83hj395ieJ/bKZ+ZUV6yB3uxpml/0TRtMXrN9QzgDueiwd4yxPZdtx0CZPHJ59KJXlviNGEE6z3l83bZr+pB3zFMmqY1aZr2RU3TJhjlKwS9lnbQtzDwd/Fn4E5N07rQuyLs1vQ+fOHoNxenU4neLJ/k8i9e07TPGDE2ozcJX4t+fnrB5Ua3Er3p0fW90ZqmbR1imzVAijp1RG//c8FQMQ/Ub3U455jRGvScwBDnzCFUon8WSYMs+99+n2+MpmnPj3ovhEdIcudd8eh9tNqNGp5vjmFd6cDtSqlwpdTn0fuYvDXIa98CZiilblB6x+9r0S9gb7i85ktKqTlKqRj0PlQvGzUp8egniUb0i9FvB1j/Z5RSK5Q+OOQe9OaooWo1+juJ3rfmUzRNO4F+Mv+TUipBKRWilJqqlFoFoJT6vFIqy3h5M/oFZzijg3egX2R/ZHyOq9GbCF8YTsBKqfFKqSuU3gm+B71p0bndk0CWchkwg/5ZNmma1q30juY3DGc7hn8DP1VKJRv7+h2XZbHo+1xvxPUV9Jq7ARnf63/RB1bEKKXmoPedGa3ngTuVUmlKqXHofaXcNcdbBHpfvXrAZnT6HnL6EsObwFyl1FVKH5hzO6cmUX2UUkuNmtVw9P5s3Zz6XQ54bA5hscu2v4d+jGw3luUDNxg14xej9yV0Ook+GCFxkPX+G7hU6VOPhAP/z1j3UAkM0DeYY/Ugy6YqfbBNqPFZfx1wzk05SenzrkUopaKUUneg16Bu6beOC9Gb4Z3nl+PAeUofiBWJfi45nZ1Am9IHrEQbscxTSrkOInoOvdn9GuNnp4fRy8lcI5ZE4/x4Wsb5aivwO2PfcoGvMfzj+FHgh0qpxUo3zbjpGtM5Zgj/Rr8GZCmlkoGfuOzPac+Zp2O89230m+RkI+6VxuJHgG8YZUUppWKVPlDMn6e5CUiS3HnXD9Ev6G3oheRTozpHYAcwHb127X+BazRNG/Ckafz9MvSLQCN6M+ZlmqY1uLzsafT+H7XoHYpvN/7+FHp1fzV6U8t2Pu054C70JpDFwJdGsT+/Q08QLEqpHw6w/Cb0C30xegL3MnptGeidmHcofTTsa+h9X4acL1DTtF70E+0l6J/j34CbNE07NMyYQ4AfoN+dN6FfoJ0J+3r0KQxqlVLOz/lbwK+VUm3oCdBIRhD+Cv17OI5+0n7aZT+KgT+h96M8iV5TsmWAdbj6NnoTTC369/74CGLp7zfAbqAQ2A/sNf42ZkYz2e3on1Uzevl5bZjvbUDv6P579ON+OoN/LgnoZbIZ/XNuRO+oD/qI0TnGsblukPcP5FX02iXnYKarNE2zGsu+i37sWdBH4/at1zj+ngdKjW2e0uSlaVoJehn7K/px+1n06Ul6hwpIKZWNfv7ZP8hLFhvL2tDL5Bc1TXNOxRGP3vrQjH4+uBi4xPW8Y9S03mvsn9N30JOuD4BvaUNMJ2Qsvwy9D99xYx8fRW85cHoN/fus1fTRuc73voI+eOcFpXcjKUIv38NxPXrLRA36YIS7tGFOYaVp2kvo5+Hn0D+7degD3MZ6jjmdR9D7+hWgl7n/9lt+unPmUG5Eb705hD6A4nsAmqbtRh+E8aCxzqPo/feEj3GOjBF+RCl1M3qn2hVmxyKE8B9KqS+hj+j+qdmxCCE8RyYxFkKIIKFpmtseiSaE8F3SLCuEEEIIEUCkWVYIIYQQIoBIzZ0QQgghRACR5E4IIYQQIoBIcieEEEIIEUAkuRNCCCGECCCS3AkhhBBCBBBJ7oQQQgghAogkd0IIIYQQAUSSOyGEEEKIADKix4+NGzdOy8nJ8VAoQnjXnj17GjRNSzM7DicpXyLQSBkTwrMGK2MjSu5ycnLYvXu3+6ISwkRKqXKzY3Al5UsEGiljQnjWYGVMmmWFEEIIIQKIJHdCCCGEEAFEkjshhBBCiAAyoj53wv9YrVaqqqro7u42OxTTREVFkZWVRXh4uNmhCD8kZWhoUsaEp0j50420jEly50P2lDfxZmEtYaGKy/MymTcxcczrrKqqIj4+npycHJRSbojSv2iaRmNjI1VVVZxxxhlmhyNM1GOzs25fNQVVLUxOieG6pZNIjBn6RBnsZWgoUsaEU1VzJy/trqKpo5fl08axdu74MZcZKX+jK2OS3PkATdO4990S/rbhGJFhIWgaPLKplB+tncU3V08d07q7u7uDulAopUhNTaW+vt7sUISJGtp7uPnxnRRVtxIfFUZbt41HNx/n8ZuXDnkTFexlaChSxgTAB8Unuf2FfXRb7cREhPH09nIumD2eB29YSFR46KjXK+VvdGVM+tz5gEc3HedvG45x/ZnZ5P/yInb/4gIuy83kD+8c4sVdFWNefzAXCpD9D3a9NgdffWIXR+vaefhLiyi86yLe+M4KIkJDuPnxndRYuoZchxxDpyefT3ArqLTwrWf3Mi09jo0/WkPBXRdx56Wz+fDQSe54uRBN08a0fjm+Rv4ZSHJnskO1rdz7bglr547nt1fOJzoilISocP78hTyWT0vl7teKqR7GxceXxcXFmR2CCGJ/XX+EwqoW7r92ARfPy0ApxbyJiTz51aV09tr5+Sv7x3zx8aTKykrWrFnDnDlzmDt3Lg888MCAr1u3bh3FxcVDru/hhx/mqaeeAuDmm2/m5ZdfHnYsZWVlzJs3b9ivF4Gvq9fO9/+dT2pcBE999UyykmMIDVHccu4UfnjRTF4vqOGNwhNmhzkmAx33d999N/fdd59JEQ1NkjsTaZrGXa8eID4qjN9eOf+UzDwsNIQ/XJ0LwK9eO2BWiH7JZrOZHYLwEeWNHfx9wzGuWjSRi+dlnLJsWno8P7hwBh+V1PPugVqTIhxaWFgYf/rTnyguLmb79u089NBDAyZxw03uvvGNb3DTTTd5IlQRhB7bcpzS+g7+eE0uSTERpyz7xqqpzJ+YyK/fKKa9R87Lw6VpGg6HY0zrkOTORJuONLDjeBO3nz+d1LjITy3PSo7hW6un8l7xSfZXtZgQoefk5+dz9tlnk5uby5VXXklzczN1dXUsXrwYgIKCApRSVFTozdJTp06ls7OT+vp6rr76apYuXcrSpUvZsmULoN9F3XjjjSxfvpwbb7zRtP0SvuX+D44QFqr4ycWzBlx+8zk5TEuP48/vH8Hh8M3au4yMDBYtWgRAfHw8s2fPprq6+pTXbN26lddee4077riDBQsWcOzYMR555BGWLl1KXl4eV199NZ2dncDgNQ579uxh1apVLF68mLVr13LixIm+v+fl5ZGXl8dDDz3k4b0V/qSl08rDHx/jgtnpnDv900+ZCw1R/PqKudS39fD0Np96WInbrF69mu9+97ssWLCAefPmsXPnTuCTa9KyZcuYPn06jzzySN977r33XpYuXUpubi533XUXoNcOzpw5k5tuuol58+ZRWVk5prhkQIWJHlx/lMzEKK47M3vQ19y8PIdHNpXy4EdH+MeNS8a0vV+9foDimtYxraO/OZkJ3PXZuSN+30033cRf//pXVq1axS9/+Ut+9atfcf/999Pd3U1rayubNm1iyZIlbNq0iRUrVpCenk5MTAy33HIL3//+91mxYgUVFRWsXbuWgwcPAlBcXMzmzZuJjo526z4K/1TR2Mm6/Gq+fu4U0hOiBnxNWGgI314zje+9mM97xSe5eN6E067T7DJUVlbGvn37OOuss075+znnnMPll1/OZZddxjXXXANAUlISt956KwB33nkn//rXv/jOd74z4HqtVivf+c53ePXVV0lLS+PFF1/k5z//OY899hhf+cpXePDBB1m5ciV33HHHGPZUBJont5XR1m3jBxfOHPQ1Cycls3JGGo9uKuXmc3KIjhj94Aqzy99gOjs7yc/PZ+PGjXz1q1+lqKgIgMLCQrZv305HRwcLFy7k0ksvpaioiCNHjrBz5040TePyyy9n48aNTJo0iSNHjvDkk09y9tlnj3m/JLkzSXFNKzvLmvjZZ2YRGTb4wR4fFc5Ny3J4aMNRKps6yU6J8WKUntHS0oLFYmHVqlUAfPnLX+bzn/88oF+ktmzZwsaNG/nZz37GO++8g6ZpnHvuuQB88MEHpzQ9tba20t7eDsDll18uiZ3o8/T2MkKV4ivLTz91wGW5Gfzp/RIe33J8yOTOTO3t7Vx99dXcf//9JCQkDPn6oqIi7rzzTiwWC+3t7axdu3bQ15aUlFBUVMSFF14IgN1uJyMjA4vFgsViYeXKlQDceOONvP322+7ZIeHXrHYHz+4o59zp45iTefrj8Vurp3LdP7fzWkE11y6d5KUI3WewwQzOv19//fUArFy5ktbWViwWCwBXXHEF0dHRREdHs2bNGnbu3MnmzZt57733WLhwIaCX6yNHjjBp0iQmT57slsQOJLkzzVPbyogKD+ELSwavtXO6/qxJ/G3DUV7cVckP1w5+hzSUsd6deMPKlSvZtGkT5eXlXHHFFfzhD39AKcWll14KgMPhYPv27URFfbomJjY21tvhCh/V1WvnxV2VrJ03gQmJA9faOYWFhnD9mZP44zslHK1rZ1r64AOAzCpDVquVq6++mi9+8YtcddVVw3rPzTffzLp168jLy+OJJ55gw4YNg75W0zTmzp3Ltm3bTvm78yIlRH/vHTjJydYefnvl/CFfe9YZKUxPj+O5HRVjSu7MKn+pqak0Nzef8rempqa+Oef6J3/O3wf6u6Zp/PSnP+W22247ZVlZWZlbr2HS584E3VY7bxSe4LLczE91QB3IxKRoVs9M58XdldjsY+tk6QsSExNJTk5m06ZNADz99NN9tXjnnnsuzzzzDNOnTyckJISUlBTeeustVqxYAcBFF13EX//617515efne38HhM/74OBJWrttfPHM4V1IPr84m7AQxQs7xz71kLtpmsbXvvY1Zs+ezQ9+8INBXxcfH09bW1vf721tbWRkZGC1Wnn22WdPu42ZM2dSX1/fl9xZrVYOHDhAUlISSUlJbN68GWDI9Yjg8Z+9VWQkRrF6ZvqQr1VKcf2ZkyioauFAjf/1H4+LiyMjI4P169cDemL3zjvv9F2XXnzxRQA2b95MYmIiiYn63Jmvvvoq3d3dNDY2smHDBpYuXcratWt57LHH+lqcqqurqaurc3vMktyZ4KNDdbT32PjcgonDfs8XlmRR39bD9tImD0bmGZ2dnWRlZfX9+7//+z+efPJJ7rjjDnJzc8nPz+eXv/wlADk5OWia1tcMtGLFCpKSkkhOTgbgL3/5C7t37yY3N5c5c+bw8MMPm7Zfwne9ml/D+IRIzpqSOqzXp8VHct6sdF4vrPG5gRVbtmzh6aefZv369SxYsIAFCxbw1ltvfep11113Hffeey8LFy7k2LFj3HPPPZx11lksX76cWbMGHlDiFBERwcsvv8yPf/xj8vLyWLBgAVu3bgXg8ccf53/+539YsGCBT08ZI7ynqaOXjYfruTwvk9CQ4c2/duXCiYSFKF4rqPFwdJ7x1FNPcc8997BgwQLOO+887rrrLqZO1R8yEBUVxcKFC/nGN77Bv/71r7735ObmsmbNGs4++2x+8YtfkJmZyUUXXcQNN9zAsmXLmD9/Ptdcc80pN2Vuo2nasP8tXrxYE2N321O7tcX3vK/Z7I5hv6er16bN+cXb2k/+UzCibRUXF480vIA00OcA7NZGcPx7+p+UL/do7ujRpv3sTe03bxwY0fteza/WJv/4DW3n8cZT/i5laHikjAWPp7eVaZN//IZ2oLplRO/78mM7tOW//1BzOIZ/7fP18rdq1Spt165dn/r7XXfdpd17771u3dZIypjU3HlZa7eV9SV1fDYvY9h3PABR4aFcMGc87xTVYg2AplkhPOXtolqsdo0rRlAzDnD+rHQiw0J4w09rFoTwltfya5ieHsfsjPgRve/S+RlUNXdREGBTe/kiSe687KNDdfTaHFyWmzni9146P4PmTitbjzV6IDIhAsM7RbXkpMYwd4gRfP3FRoZx3qx03iqqxe5jTbNC+Ir6th52lTdxWW7miB+JddGcCYSHKt4sDJwbqA0bNrBkyaenKbv77rv54Q9/aEJEOknuvOzDg3WkxkawMDtpxO9dOSON2IhQn55NXwgzdfTY2HaskQtmjx/V8ygvmZ9BfVsP+ZXNQ79YiCD0UUkdmgYXzBl6IEV/iTHhLJ82jncPnJT+mx4myZ0XWe0ONpTUsWZWOiEjaJJ1igoPZfm0cXxcUj+ighHshSjY9z+YbD7aQK/dwXmzR37hAVg1PY3QEMWGkvpT/i7H0OnJ5xM8Pjx4kozEKOZkjKxm3On8WelUNHVyvKFj2O+R42vkn4Ekd160u6yZ1m4bF4zywgOwZlY61ZYujtS1D+v1UVFRNDY2Bm3h0DSNxsbGAefFE4Hnw4MniY8KY2lOyqjenxgTzuJJyXxU8snUBMFehoYiZSx4dFvtbDrSwHmz0kdVMw70TZ3yUb8bqMFI+RtdGZNJjL1o/aGTRISGsGKAZ/AN1+qZ+ns/OlTHjPFDd2bNysqiqqqK+vrhFaRAFBUVRVZWltlhCA9zODTWH6pn1Yw0wkNHf9+6elYaf3ynhLrWbtIToqQMDYOUseCw43gTnb12Lpg9ftTryE6JYVp6HBtK6vjaitM/PQbkGuY00jImyZ0XbTrSwJKcZOIiR/+xZyRGM2tCPB+V1HHbqqlDvj48PLxvFm0hAtmh2jYa2ntYNWP0N08Aa2am88d3SthwuJ4vLMmWMiSEYdPheiJCQzh7mPNHDmbNzDSe3FpOR4+N2CGuh1L+RkeaZb2kob2HQ7VtLJ82bszrWjUzjd1lzXT02NwQmRCBYeuxBoAxl7FZE+IZnxDJxsPBXVMgRH9bjzWyaHIS0RGDPw99OFbNSKfX7mDHcZn5wVMkufOSbcb0JedMHdsdj76OcdgcGrvLZUSfEE5bjzVyxrhYMpOix7QepRTLpqSyvbQpqPv5COGqqaOX4hOtLJ869gqKJbTFKG8AACAASURBVDnJRISG9F0XhftJcuclW481EB8ZxvyJiWNe15LJyYSFKCkYQhisdgc7ShvdcvMEsGxqKg3tPRwd5sAlIQJdXwWFG1qfosJDWTgpiW2lcg3zFEnuvGTL0UbOmpJK2Bg6ejvFRoaxIFsKhhBOhVUWOnrtbun2AHrtOCBlTAjDlmMNxEWGkZc19goK0G+gDtS00tJpdcv6xKkkufOCyqZOKpo6WT7NPbUKoBeMouoW2rqlYAix5WgjSsGyMXb0dspOiWFiUrTUjgth2Hq0gbPOSHFLBQXoZVXTkH53HiLJnRc4LxDuqlUAvWDYHRq7yprctk4h/NXWYw3MyUggOTbCbetcNjWV7aWNOORRZCLI1Vi6KGvsdEuTrNOCSUlEhYdI7biHSHLnBTvLmkiOCWd6epzb1rlost4hdetRKRgiuPXaHOyrsHDWGe6rGQf9Bqq508qh2ja3rlcIf+OsRDjrjNFNDj6QyLBQlkxOkdpxD5Hkzgv2lDezeHLKqGf0HoizQ+qO41JzJ4LbgZoWemwOluQku3W9y4zBGdJsJILd7rJmYiNCmTVh6InzR2LZ1FQO1bZh6ex163qFJHce19Dew/GGDpa6+cID+nDy4hOtdPbKfHcieO0u06cEWjLZvWUsMymazMQo9lZY3LpeIfzN7vJmFk1Odlt/O6fFRpndJ2XM7SS587C+C48HkrvFk5OxOzQKKlvcvm4h/MXu8iYmpcSQnuD+Z5sumpzMHunXKoJYa7eVQ7WtfYmYO+VlJREaothdLmXM3SS587A95U1EhIUwzw3z2/W3aJJe2PZWyGTGIjhpmsae8ma319o5LZ6cTE1LNzWWLo+sXwhft6/CgqbBksnu62/nFB0RytzMBPbIhPxuJ8mdh+0ubyZ3YiKRYWN7XMtAkmIimJYeJwVDBK2yxk4a2ntZ7IGacfjkgiY3UCJY7SlrIkTpo1s9YfHkZAoqW7DaHR5Zf7CS5M6Duq12iqpbWJLj/jsep8WTktlT3izTNYigtNtoMl3qoTI2KyOe6PDQvu4VQgSb3eXNzMlMIC4yzCPrXzw5mS6rnYMnWj2y/mAlyZ0HFVRasNo1jzUZgV4wWrqslDbIY5JE8Nld1kxCVBjT0tw3zZCr8NAQ8rITpeZOBCWrXZ9myBNNsk7OvnzSAuVektx50G7jYPVER1SnRVIwRBDbXd7E4snJhIS4b5qh/hZPTuZAjYxKF8Hn4IlWuqx2j17DMhL1UelyDXMvSe48aF+FhSlpsW6dNb+/qWmxJMWES8EQQaely8qx+g6PXnhA73dnd2gUVsmodBFcnFOUeLqMLc5JYa9cw9xKkjsPKqyykJflmU6oTkopFk9K7qslFCJY7DeSrVwPl7GFRkdyuYESwaagysK4uEgyEt0/zZCrxZOSZFS6m0ly5yEnW7upa+shN8v9U6D0t3BSEqX1HbR0WT2+LSF8RWG1Xqvg6TKWFBPBlLRY8itlolURXPZXtZCXlejWpysNZKExrVeBlDG3keTOQ5wHqTeSO2fNxX5pNhJBpLCyhUkpMSTFeK7bg1NeVpJceERQae+xcbS+nfleuIbNyognPFSRXyVlzF0kufOQ/dUthIYo5mR4I7nTt1EgBUMEkf3VLV65eQLIy0qkrq2H2pZur2xPCLMdqG5B0/B41yKAyLBQ5mQkUChPW3IbSe48pKCqhRnj44mOcP/kxf0lxUSQkxojNQsiaDS091Bt6fLKhQcgN1vfjjTNimDhHEDkjZo70Fug9le3yJytbiLJnQdomkZhlYVcDzxybDB52Ukymk8Ejf1evvDMyUggLERRKLXjIkgUVFmYmBTNuLhIr2wvLzuJ9h6bzNnqJpLceUBVcxeWTiu52V5M7rKSqG3t5mSrNBuJwFdQZUEpPPLM5oFEhYcyOyNBuj6IoOHNbg8AC4zrZb40zbqFJHce4LwA5E70TpMRQJ5RMKRpVgSDwqoWpqbFeeyRSAPJzUqksFKajUTgs3T2Ut7Y6bWacYAp4/TyLNcw95DkzgP2V7UQERrCzAnxXtvm3MxEQkOU1CyIgKd3e/BurQLozUZtPTaON3Z4dbtCeNv+ar32zFt9WgFCQhTzJyZK1wc3keTOAwqqLMzOiCcizHsfb1R4KDPHx0u/OxHwTrR009De49U+rfDJhU5qFkSgc15HvNXtwSk3O5HiE6302Oxe3W4gkuTOzRwOjaLqVo/Pmj+QvGx9Li5Nk2YjEbicFx7nCFZvmZYeR0xEqNxAiYBXUGnhjHGxJEaHe3W7C7KSsNo1Dp1o8+p2A5Ekd25W2tBBe4/N601GoM/F1dpto6yx0+vbFsJbCqsshIUo5mQkeHW7oSGKeRMTZToUEfC8PZjCyXnDJt2Lxk6SOzdz9hcwq+YOpNlIBLb91focklHhnp9Dsr8F2UkUn2il1+bw+raF8Ia6tm5OtHQz38tNsgCZiVGMi4ukQEbMjpkkd25WWNVCdHgo09LjvL7t6elxRIWHyF2PCFjOwRR5XpxmyFVuViK9NgcltdJsJAKTcw7JPC93ewBQSpGXlSjXMDeQ5M7NCqsszJuYQGiIZx+0PJCw0BDmT0yUmjsRsCqaOmnpsjLfi9MMuXIOqpBnYIpAVVDVQoiCuZne7fbglJedxLH6dtq6raZsP1BIcudGNruDAzXmDKZwys1K4kBNKza7NBuJwFPgHExhQn8ggKzkaFJiI9gvyZ0IUPurLExPjycmwntzSLrKzUpE06CoutWU7QcKSe7c6PDJdnpsDtMuPKAXjB6bg8Mn5REuIvDsr7IQEebdOSRdKeWci0v6BInA4+z24M3Ji/tzVo7IfHdjI8mdG5k5mMJJCoYIZAVVLczJSCA81LxTV15WIodPttHZazMtBiE8odrSRWNHL3kmJncpsRFkJUdTWC03UGMhyZ0bFVa3EB8VRk5qjGkx5KTGEB8VJgVDBBy7Q6OousXUCw/A/KwkHBoU10izkQgs+/u6PZhXQaFvX55UMVaS3LlRYZWF3KxElPL+YAonpZQUDBGQSuvb6ey1M9/kC48zuSyQplkRYAqqWggPVczKMKfbg1NuVhKVTV00dfSaGoc/k+TOTbqtdkpq20y/4wG9YJTUttFtlUe4iMDhTKbMrrlLT4hiQkKUDKoQAWd/tYVZExKIDPP+HJKunP3W90sL1KhJcucmh2rbsNo1rz/vciC5ExP1R7jIXFwigOyvshATEcqUNO/PIdnf/CwZVCECi8Nh/mAKJ+czbQtlWq9Rk+TOTZx38d5+3uVAnDFI06wIJAVVLcybmGjKHJL95WUlUtrQQUuXzMUlAkN5Uydt3TbTa8YBEqLCmZIWK33Hx0CSOzcpqGohNTaCzMQos0MxHuESITULImBY7Q6KT7T6RM040Nfv74BcfESAcFYGmDVBeH+5E6Xv+FhIcucmvjCYwumTubikYIjAUFLbRq/N4RM140BfkimDKkSgKKhsISo8hBnjze/2AHrf8ZOtPZxs7TY7FL8kyZ0bdPTYOFrX7hODKZxys5I4WtdOR4/MxSX8n7NjtS80GQEkx0YwKSWG/dVyAyUCw/5qC3MzEwkzcQ5JV87nR0sL1Oj4xrfo5w7UtOLQzHsk0kBysxJxaHpsQvi7wioLidHhTEoxbw7J/uZnJVJQKRce4f9sdgdF1a3M95FuDwBzMvT+tdICNTqS3LlBX18Fn0ruZFCFCByFVS0+0+3BKS8rUZ/Rv73H7FCEGJNj9R10We19tWW+IDoilOnpcVJzN0qS3LlBYVULGYlRpMebP5jCKS0+kszEKCkYwu8555D0pVoF+KTjuYzoE/6uwMcGUzg5J+TXNM3sUPyOJHdusL+6xaeaZJ3my5MqRAAoPtGKzaH5VJ9W0MuXUlAoTbPCz+2vaiEuMowp42LNDuUUuVlJNHdaqWruMjsUvyPJ3Ri1dFk53tDhcxce0AtGWWMnLZ0yF5fwX58879K3bqDiIsOYmhYngyqE3yussjBvYgIhPjCHpKu8vu5FcgM1UpLcjZGvXnjgk4Ihj3AR/qygysK4uEgyfGAOyf5yJyZSUNUizUbCb/XaHBw80dZ3vfAlMyfEExEaIi1QoyDJ3RgVGnftuT7WVwHo66NUIAVD+LH9VS3k+dhgCqfcrETq23qolbm4hJ8qqW2j1+7wydaniLAQZmfES83dKEhyN0aFlS1MTo0hMSbc7FA+JTEmnJzUmL7aRSH8TXuPjaP17T41Et3VfGk2En7OefPvi61PoPdtLapuweGQ2vGRkORujAqrLD43is/V/KwkqdIWfquougXNx+aQdDU3M0Hm4hJ+rbDKQnJMOFnJ0WaHMqDcrCTaemwcb+wwOxS/IsndGNS39VDT0s0CH3kk0kDyshKpaemmvk3m4hL+x5k0+WJ/IICo8FBmjJdmI+G/9Dkkk3yy2wN8cmMnN1AjI8ndGPRdeHw4ucvtG1QhBUP4n4LKFrKSo0mNizQ7lEHlZSWyv1oGVQj/09lr4/DJNp++hk1LiyM6PFRuoEZIkrsxKKi0EKL0phlfNTczgRCFPCZJ+KWCKovP1to55WYlYem0UtHUaXYoQoxIUbX+6ExfeWbzQMJCQ5g3MUGSuxGS5G4MCqpamDE+npiIMLNDGVRsZBjT0uNkOhThdxrbe6hq7vKpRyIN5JNmIyljwr8U9g2m8O0bqPkTkzhQ04LN7jA7FL8hyd0oaZrmF7UKoBcMeYSL8DfOZMnXy9iM8fFEhMlcXML/5FdamJgUTVq873Z7AMjLTqTb6uBIXbvZofgNSe5GqbKpC0un1af7KjjlZSfS0N7LiRaZi0v4j3yj28M8Hx6NDs65uKTZSPifwqoWn68Zh0/mbJUbqOGT5G6U8n18biBXuX1zcUnBEP6jsMrC9PR4YiN9t9uDU54xF5dd5uISfqKpo5eKpk6fb5IFyEmNJT4qTG6gRkCSu1EqrLQQGRbCzAnxZocypFkT4gkLURRIwRB+Qu/20OIXN0+g30B19NoprZdmI+EffH2aIVchIYrcrERJ7kZAkrtRKqiyMDczgfBQ3/8Io8JDmZURL0+qEH6jqrmLpo5ev+j2ADKoQvifgsoWlMJnn/7S3/yJSRyqbaXHZjc7FL/g+5mJD7LZHRRVt/rNhQdkUIXwLwV+VKsAMDUtjpiIUOn6IPxGYZWFaWlxxPlBtwfQuz5Y7RqHTrSZHYpfkORuFI7UtdNltfvNhQf0gtHabaO8UebiEr6voNJChJ90ewAIDVHMy0ykUKYcEn7AOduDP/S3c5ovT6oYEUnuRsEfnkzRn7MQF0jBEH6goKqFORkJRIT5zykqNyuR4ppWrDIXl/BxNS3dNLT3ssAPRso6TUyKJjU2Qro+DJP/nDl9SH5lCwlRYeSkxpgdyrBNHx9HZFiIFAzh8+wOjaLqFp9+ZvNA5mcl0mNzUFIrzUbCtxVU+sfkxa6UUsyXQRXDJsndKBRWWcjL9t0HLQ8kPDSEuZkJMqhC+Lyjde109tr9Yv4tV3l9z3GWMiZ8W0GVhYjQEGZl+Ee3B6fcrCSO1LXR2WszOxSfJ8ndCHVb7RyqbfObKRpc5WYlUVQjc3EJ3+aPtQoAk1NjSIgKkz5BwucVVFqYnRFPZFio2aGMSF5WIg4NDtS0mh2Kz5PkboT2GxOVLshONjuUEcvNSqSz184xmYtL+LB9lc0kRIVxRmqs2aGMiFKK3KwkaTYSPs1md1BY5X/dHsB1UIWUsaFIcjdCe8qbAVg0yf8KhrMmJL9SahaE79pT3syiycmEhPhPtwen3KxESmrb6LbKXFzCNx2qbaOz186iyf5XQZEeH0VGYlRf7b4YnCR3I7SnvJkzxsWSGufbD1oeyJRxsSRGh7PXSFCF8DUtXVYOn2xn8ST/u/AALJqUjM2hycVH+Ky9Ffr5f7EfJnegl7E9cg0bkiR3I6BpGnvLm1nkpxeekBDFoklJ7JaCIXzUPn+/8Bhx76mQMiZ8057yZsYnRDIxKdrsUEZl0eRkqi1d1LZ0mx2KT5PkbgTKGztp7Oj12wsPwJKcFI7WtWPp7DU7FCE+ZW95MyHKv+aQdJUSG8GUtFj2lElyJ3zTnvJmFk9O9qvZHlwtcd5ASSXFaUlyNwLOg8mfkztnreNeqVkQPmhPRTOzMxKI9ZNHIg1kyeRk9lQ0y6P+hM+pa+2mqrnLb1ufAOZkJhAVHsLu8iazQ/FpktyNwN6KZuIjw5ieHmd2KKO2IDuJsBDFbqlZED7G7tDIr7D49c0T6Dd/lk4rx+o7zA5FiFP4e3870OdszctKkpq7IUhyNwJ7yptZ6Kej+JyiI0KZm5kgBUP4nJLaNjp67X5dqwCweHIKAHukZkH4mD3lzUSEhTA30//maXW1JCeZAzWtMpnxaUhyN0xt3VZKTrb57Sg+V4smJ1NQZZFnYAqfsicAahUApqbFkhQTLjdQwufsKW8md2KiXz2zeSBLJqdgd2gUVMp8d4Px72/Yi/IrLWia/194QC8Y3VaHzPItfMre8mbS4iPJSvbPUXxOSikWT0qWUenCp3Rb7RRVtwbENWyhMc+s1I4PTpK7YdrTN4rPv6uz4ZMEVWoWhC/ZU97M4kn+O4rP1eKcZErrO2jqkFHpwjccqGmh1+7wy8mL+0uKiWB6epxcw05Dkrth2lXWxKwJCcRHhZsdyphNSIxiYlK03PUIn1Hb0k1FUydLcvz/wgN67TggE4YLn7HzeGB0e3BakqNPZuyQZ6UPSJK7Yei1OdhT3sxZU1LMDsVtFk9OZleZTNcgfMOO440AnHVGqsmRuEduViLhoYpdcgMlfMT20kampccxzg+frjSQRZOSae22caROnpU+EEnuhqGwykK31cHZUwLjwgNw1pQU6tt6ON4g0zUI820vbSI+Mow5mQlmh+IWUeGh5GUlsaNUkjthPpvdwe6yJs4OoAoK5/XYeWMoTiXJ3TBsL9UPnjNzAqdgLDMKxrZSKRjCfDtKG1l6RgqhfjzNUH9nT0llf3UL7T0yXYMwV1FNKx299oCqoMhKjmZiUjTbjsk1bCCS3A3DjuNNzJoQT3JshNmhuM0Z42JJj4+UgiFMV9faTWlDR0DVKgAsm5qK3aGx67jU3glz7XBWUJwROGVMKcXZU1LZXtoo/e4GIMndEHptDnaXNQfUHQ/oBWPZ1FS2lzZJvzthqu1G8hMo/e2cFk1KJiI0pK/mXwizbC9tZGpaLOnxUWaH4lbLpqbS3GnlcF2b2aH4HEnuhrC/2kKX1R5wtQqgN802tPdwrF46pArzbC9tJC4yjLkB0t/OKToilAXZSdL1QZjKZnewKwArKIC+67K0QH2aJHdD2G50iD4zwGoVQL/rASkYwlw7ShtZkpNMWGjgnY7OnppKUXULrd1Ws0MRQar4RCvtPTbOCsDkLis5huwU6Xc3kMA7m7rZ9tJGZo6PJyWA+ts5TUqJISMxqi+BFcLb6tq6OVbfEZC1CqDXjjs0pN+dMI2zW8DZAdTfztWyKansON4k/e76keTuNLqtdnYeb+qr4Qo0SimWGR1Spd+dMMOWow0AnBOgZWzhpCQiwkKkZkGYZtORBqalx5GeEFj97ZzOnpJKS5eVg7XyOE1Xktydxq6yJnpsDlbNSDM7FI9ZNjWVxo5eDp6QDqnC+zYebiAlNoJ5mf7/WL+BRIWHsnhSMpuNJFYIb3JWUKycHrjXsHOmjgM+uVEUOknuTmPj4XoiQkMC6skU/TkT148P15sciQg2DofGpiP1rJg2jpAAmt+uv1Uz0zhU28bJ1m6zQxFBZsdxvYJi5YxxZofiMRMSo5g5Pl6uYf1IcncaGw83sPSMZGIiwswOxWPSE6KYnZHAx4frzA5FBJniE600tPeyMoBrxsHlBqpELj7CuzYericiLCTgphnqb9XMNHYdb6ZDJgzvI8ndIGpbuik52RbQ1dlOq2aksbusWWbSF1616YjejLJyeuDWKgDMmhDP+IRIqVkQXrfxcD1nnZFCdESo2aF41KoZafTaHTKnpAtJ7gax8Yh+Ig70WgXQC4bNoUmfBeFVGw/XM2tCfMB29HZSSrFqRhqbjtRjszvMDkcEiRpLF0fq2jk3wG+eAJbkJBMTEcoGqR3vI8ndIDYerifd1smsE0fNDsXjFk9OJi4yTGoWhNd09NjYXd4U0IOVXK2akU5rt438SovZoYggsSmIKigiw0I5JyeJDR/lo/X0mB2OT5DkbgA2u4PNJXWce3ArqjHwq3kjwkI4Z2oqH5fUy5Qowiu2HWvEateC4sIDsGL6OEJDlNxACa/5+HA96fGRzBwfb3YoXrFqcgKV0cmU/fdts0PxCZLcDWBPeTOWHjvnH9tldihes2pmGtWWLnkUmfCKDw6eJD4yjKU5gTsS3VVidDgLs5Ok2Uh4RY/Nzscl9Zw/Ox2lAnckuqtVZyQBsOGNLSZH4hskuRvAewdqibBbWVlRYHYoXrNmZjoA7xWfNDkSEejsDo0PDp5k9ax0IsKC5xS0ZlY6+6tbqLF0mR2KCHBbjzXS0WvnojkTzA7FayYlRTGtsZL3tRQIgha3oQTPmXWYNE3jvX0VLK/cT1yAjzBylZkUTV5WIu8W1Zodighw+yqaaWjv5cI5480OxasunqdfaN87IGVMeNb7xSeJjQgN2KcrDebiozvYkTWX5qdfMDsU00ly18+h2jYqOx1cVLLV7FC8bu28CRRUSc2C8Kz3ik8SHqpYPTM4+ts5TU2LY3p6HO8ekNpx4TkOh8b7xSdZPTOdqPDgqaAAWHtsB/aQUD54a7vZoZhOkrt+3i+sRmkOzj8SfAfH2rlSsyA8S9M03jtQy7Kp40iICjc7HK9bO3cCO8uaaOroNTsUEaDyqyzUt/UEXc04wLy6Uia21PFu0jQoLjY7HFNJctfP29uPsPDkUdI7gm/KAmfNwjuS3AkPKTnZRlljZ1BeeEBvmnX2ORTCE94tqiUsRPX1ow4mCrjoyDY2Tsql/V9PmB2OqSS5c3HkZBsHu0K5rGiD2aGY5uJ5E9h5vIn6NpkrSLjfa/k1hIYoLpkXPB29Xc3NTGBiUjRv7z9hdigiADkcGq8X1LByRhqJMcFXMw5wcclWesMiWL+xCOx2s8MxjSR3Ll7bcpgQh53LDm0yOxTTfDYvE4cGbxTWmB2KCDCapvFqfg0rpo1jXFyk2eGYQinFZ/My2XikgcZ2uYES7rW7vJmalm6uWJBpdiimWVJ9kAltDbw69Wz48EOzwzGNJHcGTdN4dXcF51QdIL2j2exwTDNjfDxzMhJYt6/a7FBEgNlb0Uy1pSuoLzwAVy6ciN2h8Uah1N4J93o1v5ro8FAumB2c3R4AQjUHVxzYwMfZeTQ+8oTZ4ZhGkjtDfqWFCkcElxetNzsU0121aCIFVS0yobFwq1fza4gMC+GiucHZJOs0c4J+A/WK3EAJN+q1OXhz/wkunDOe2Mgws8Mx1ZUHPsIWGsabFV3Q2mp2OKaQ5M6w7v0CIuxW1h7eZnYopvtsXiYhCqm9E27Ta3PwRuEJLpg9nrggv/CAXnuXX2mhVG6ghJtsPFyPpdPK5XnBXTMOMKuhnFl1x/nv7FXw8stmh2MKSe6Arl47/z1s4eLD20js6TA7HNONT4hi+bRx/HdvNXaHPGtWjN17xbU0dfTy+SVZZofiEy5foN9A/WdvldmhiADx/M4K0uIjWRVk80cO5qoD68mfMJ2jjwXnhMaS3AFv7KukTYVzw763zA7FZ1y3dBLVli4+PlxndigiADy3o4KJSdGsnC4XHtBvoNbMTOfFXVX02hxmhyP8XI2li49K6rh2STbhoXJZB7iqaD3hdivPxk6F48fNDsfr5CgAnvugiCnN1ZxVWWR2KD7jornjSYuP5JntFWaHIvzc8YYOth5r5PozswkJCY6HmA/Hl5ZNpqG9h/eKZV5JMTYv7KpEA65dmm12KD5jXGcLl5Rs4eU5a+h8/Cmzw/G6oE/uDp5oZV+b4oZ9byOXnU+Eh4Zw/dJsPiqpo7Kp0+xwhB97YWcFYQq+sCDD7FB8yqrpaWSnRPP0+wfMDkX4MZvdwb93VbJyehrZKTFmh+NTbtz7Fm2Rsbz+URFowdXFKOiTu8fWHyLK2sPV+4N3PpzBXH/WJEKU4tkdUnsnRqe9x8bzO8q56PA20o9IEuMqJETxxaXZ7KjvpeTVD8wOR/ipN/efoLa1my+dPdnsUHzOkupiZtaX8dSMVWibN5sdjlcFdXJ3oqWLdfvruPbgRyR3t5kdjs/JSIzm4rkTeHZHOa3dVrPDEX7ohZ0VtPbY+fq2l4Luznk4rl08kejebv7x/Eaw2cwOR/gZTdN4+ONSpqbFcv6s4Hvc2FAUcPOe1zmQdgZbnnzV7HC8KqiTu8c2H8fhcHDLln+bHYrP+ubqqbR123h6W7nZoQg/02tz8OiHhzi7qogFbfLEk4Ekx0RwQ8E7vDppCRX3P2x2OMLPbDzSwMETrdy2cqr0Zx3EVUUfMqGtgQe706Cry+xwvCZokztLZy/PbSvjssNbyW6VEaGDmTcxkdUz0/jX5uN09krNghi+dXsrqe3W+MbWf0NoqNnh+Kyv73qFUIeDhzcchfp6s8MRfuThDccYnxDJFQtlbrvBRNpt3Lrzv2zPmsee5143OxyvCdrk7u8bjtFpdfDNnf8xOxSf9+0102jq6OWZ7VJ7J4an22rngVfzmVdXyqrje80Ox6eNb2/imv0f8PKs1VT/6JdmhyP8xKYj9WwrbeTWc6cQGSY3T6dzfcG7pHS2cP+O4BmZHpTJXbWli8c3l3LVwY+ZdeKo2eH4vCU5KayckcaD64/S3NFrdjjCDzy9/iDV9jB+uv5fMgp9GP5n279RmoN7mxNg3z6zwxE+zuHQ+P3bh8hKjubGZTKQYigx1h6+tf0lNqVM5eO/B8ekxkGZ3P3pzSKwWvnBhifNDsVv3HnpbNp7bDzw4RGzQxE+rqXTyoMfQe+wjAAAEKxJREFUHmFV2T6WlxeYHY5fmNhWzy271rFu1kryv/cLGXwiTuu1ghoO1LRyx9qZUms3TDfufYPJzTX8714LtpLDZofjcUGX3O0pa+KVwpN8Zd+bTGyT/i3DNWN8PNefOYmnt5dz5KSMLBaDu/f5rbQRxo/XP2Z2KH7lmzteZlx7M7/OOQ/Hc8+bHY7wUa3dVn739kHmTUzgs7nS1264Iu02frrhcQ6nZvP8d38HPT1mh+RRQZXc9djs/PixTWS2NXD7xmfMDsfv/ODCGSREhfHDlwux2eWRSeLTdh9v5Jkj7dy893Xm1AffI3/GIq63i598/Dh7M2by5MOvQXu72SEJH/THdw5R39bD/35uvoyQHaG1h7exvCyfP8z5DFXf+7HZ4XhUUCV3D7y4naO9Yfzm3YeItXabHY7fSY2L5FdXzKOg0sIjm+TCLU7V2WvjR09uZWJrPf9v49Nmh+OXri5az5pju/jDmV+g7O4/mB2O8DFbjzbwzPYKbj7nDPKyk8wOx+8o4Pdv/wVNgx+3jkd7ZZ3ZIXlM0CR3mw9U8/fCJr5Q+B5rSnebHY7f+mxuBhfPncCf3z9MQaXF7HCED7nrvwUc74J73/yz3DyNkgJ+986DhNutfLchle5Dgd83SAxPQ3sP330xn6lpsfxw7Qyzw/Fb2a11/Oyjx9iSncuj//ciVATmE5iCIrmrbOrke09sZ2pzDXd/8A+zw/FrSil+e9V80uIjue3pPdS1yUVcwHM7Kngpv5Zv73qFcyoKzQ7Hr01ob+Tetx+gYMJ0fnHfK2gyuCLo9docfPu5vbR0WXnwhkXERISZHZJfu6HgHS4u2cLvll3P5q98H6yB9wSmgE/uWjqt3PyXD7FabTz8n98QYw3sTpTekBIbwT9vWoylq5fbnt5DR49MbhzMPiqp4xfr9rPm+F6+KyPQ3eLiw9v4ztYXeGncPP75yNtmhyNMpGkaP/lPIdtLm/jD1fOZnZFgdkh+TwH3vXU/0xqr+Pa8z3Pkp/eYHZLbBXRy19Jp5csPb6Kyw84//3MP05qqzA4pYMzNTOT+axdQUGnha0/uoqvXbnZIwgRbjzbwP0/tYlZ9GQ++/kfCNBlo4y7f3/Qslx7axO9KNZ74v+CYm0ucyuHQ+NXrxfx3XzX/78IZXLkwy+yQAkZcbxeP/Pcewu1Wru+ZybF175kdklsFbHJX19rNtf/YSvGJNv76+r2cVXXA7JACzsXzMvjztQvYebyJmx7bQWO71IoGk3eLTnDzo9vIPlnB4y/cSWyXjO50pxA07n/9Pi46vI276+J58Ju/RevsNDss4SVWu4M7Xi7kia1l3LLiDL593jSzQwo4ky21PP/CzwCN6z6sJz//mNkhuU1AJndbjzVw6f0fU36imcde/x1rS7aaHVLAumLBRP5y/UIKq1q4/MEt7K9qMTsk4WE2u4P71uVz2zN7mV17jBef/RHpHTK4xhPCHXYefPUPfO7AR9yXmMftt9xHW9FBs8MSHlZj6eLaf2zjP3ur+MGFM/j5pbNRSqY98YRpjVU8//zPiLT18IVn9/PSum0B0c81oJK71m4rv3p5H1/653biq8p55Zk7WFGyw+ywAt5luZn8+7Zl2B0an/vbFu57t4RuqzTTBqKi6hauue99HtxezbX7P+CFZ39CUrfU2HlShMPGn9/4Ez/a8ARvZi1g7d938tE/XzI7LOEBDofGczsquOSBTZTUtvHX6xdy+/nTJbHzsOmNlbz25PdZVH2QO7Y3cevt/+BEsX/X4gXEkJuOHhvPbSnlH+8fpNERwg37P+CnHz5CXG+X2aEFjbzsJN793kp+/UYxD350lJf2VPLt86ZzzaIsoiPk8Tj+7mhdG3/fcIxX9laR0tnCAx88whUHPzY7rKChgG/teJmzK/dzx2e+x1dKx7H8B4/z3dsuZumMCXLx93MOh8Z7xbX8df1RDtS0cuYZKfz+qvlMSYszO7SgkdLVyrMv3MnjSy7n3nNvZPVjhdykvc0tt1/D+MkTzA5vxPw2ueux2dlT1sy6fdW8ta+SdofinMqD/Hj9Y+TVyvNPzZAYE86fvpDH55dkce+7JfxiXRF/fOcQVy6cyCXzMliSk0x4aEBVFge0hvYe1h+s4z+7K9hRbiHKYeUrhe9y+4anSezpMDu8oLSopoS3Hv8Ozyz8DA+dcx1feHwvM6IdfGHZFC5YlEPOuFizQxTDpGkaJSfbeKPgBK/sq6ba0kVOagwPXLeAy/MyJWE3Qajm4JZd61h7eBt/Xn4Dj85bw2MP7eT8iHY+d+Vyzp2bSXxUuNlhDosaSdvykiVLtN27vT8BcGevjarmLg6fbONwbRv7yhrZVd5Mtx1ird1ccnQ7X9z1GgtPuHnCz8REeOkluPBC9643CGiaxs7jTTy/s4K3imrptTmIjwpj8eRkcrOSmJeZwOTUWLJTok2bs0kptUfTtCWmbHwAZpUvTdNo7rRSWt/O4ZPtHDrRys5DNRyy6HM/ndFcw1WHN3PD7tdJbW8e3UYSE+Gdd+Dss90YeQCw2SAiAkbRx6czPJLXZ6/kuUWXUjBe72yfQxcLJyaQt2AqMzITmZQSQ0ZiNKEmPaZKypiu1+bgZGs3R+vaKTnZxv7qFnaUNtLQ3kuIghXT0/jCkiwumZdh2nfl95qaIDPTrc+MLf//7d1rjFx1Hcbx7+/MZWdmZy9lu9uGFgqtBWmgQqtFIYgRNUUbeUFI0WhCfKsJJIgRie9MfGGivtBEjTHeMGJENDEhUgUhTcBQuV8KbdB2qS3dZbftXuf688WZvXZ3221n55yZfT7J2dk5M5nz7GSfOf8558yZ7vX87kO7+eP2T/N+roskzvVrkmzf0sf2Leu4cm07l1+SozuXimwwvljHlrVWHRor8pvnjoA7Tvh6NDU4nL4Ocw5GDOf5rNtq86fmuTNRqjBerDBRDC/HixVGJksMnJ7g5Mgko6WZxwuqVba+38/dx9/kpv5XueXgc2TL+pRm3JgZN27u4cbNPXynUGb/oUH++dZJXjx6imfePkR11rqspz3NJe1p1uTSdOVSdGVTZFIBbckEbcnwMp0MSCWMwIzAwsefujRjZj7hdTOj2V4iB0YK/PTpmeM85q/uF1r/+7x7LXSfcsUZL5Wn+xVelhk6Pc7ASIGB8TKzKkZ7cYIb3jvEA0df5ZbDz3PdicNN91yuBrlSgb2v7GPvK/s42rWOJ7d8hP2bd7B/bCuPHZs5KWsSZ10moDubojvfRndXjo5sOuxWKuxYOhHQlgpIBEGtRxAENt2jmd7NdEsdm1GorcPGSzP9GitUGBwtcHKkwNBYcc79L+3KcMvWXm7a0sOtV/XS15m50D9JVtCmUyd48Olf8sAzv+aFDdfwj6s/xoGN23h4YJxfHDg+fb+OoMraXIruzhzd+Qzd7Wmy6bnrr7baFMxbZzFr3RXU5lO7vJh+LWtwd+zUBN/+82sXsbiFJaoVcqVJssVJsuUC2VKBjsI414wN8fHRYfpGh7h0ZJCrBo+weeQkGZvVrnQC0rm6ZwLq+g5gNcu3Jdl97Xp2XxsetzBWKPP2eyP0D0/QPzTOu8MTDI8VOTVRpH9onNcmShTKVQqlCsVKlVKl+T+5dD5OnJnku48fXJHHTpeL5EqT5EqTZEoFcqUCaybOsGV0mL6xIfpGh7li+H9sHTzKhsnTBMlZx0nm6tCvSX2TyaLcL/o5vrw0wj0Hn+Seg0/iwIlcN++099LftY6j3es5ke/hdLaD4WwHxzN5RtraKSZSFBMpJlNp3FbmcImdsHNFHvgCrVTHzCCbSpBLJ8imE7Xfk2xck2PHpjX0dbSxvjPDB/rybO3roCvXHLv2mk6pVJ/Xq3mSwK6h/7Dr2fA71UtBgnc613Ek30t/Zy/9XesZbO/mdCbPYKaDQ9kOJlNtFBIpCsk0xWS67pmmLNaxZe2Wvb6j05/YcWP4rg0PL92nN0eagTlQe+dntR8WBLXJsCCBBQaJBAQBQSpFKp/D8nlob4d8PpxyOQgiPj4rCOCuu8I8EplK1SmWqxTL1ektvlV3qrO3Ck/Pm7l+LpvWtsdql9GOfIfvv35HeGWqU/O2GtjUbbOnefMslYRUKpySKRKpJMlUMtz9l81Cby/09IS3N0oyCXv3QkZbKM7y6KNwqoGnkqlWYXgYBgbgzBm8VKJcKlMoV6mUK3ipTLUcTl716b0xVah1bmYvzaJqBdz91ONjb7jH5gV0x86dvv/ZuWdQmL995Hz2rs2/TzoR6Bi5qLnDI4/AWMTHA7vDxASMjoZZRkepjo5SHBunWCjh1dp6qlIN11mVKrhTnZoPVKde+adLtvBKzYHbF+nYsrbcJa++ir6n9y37bxW5GInAwnfDLf6p2+CDV5Pbr0+grjp33hnp4g1I1aZ6e9NsZTZFX6DATN/L2qrM4O67o06xoADI1KZ6W6xj+uiiiIiISAvR4E5ERESkhWhwJyIiItJCNLgTERERaSEa3ImIiIi0kGWdCsXMBoAjdVjuWmCwDo9TL8qztFbNs8nde+vwOHWhfjVU3DK1ah51rDGUZ2mtnGfBji1rcFcvZnYgTucXU56lKU9zidvzE7c8EL9MytNc4vb8KM/SVmMe7ZYVERERaSEa3ImIiIi0kKgGdz+LaLmLUZ6lKU9zidvzE7c8EL9MytNc4vb8KM/SVl2eSI65ExEREZGVod2yIiIiIi0k8sGdmd1vZm5mayPO8T0zO2hmr5jZY2bWHUGG3Wb2lpkdNrNvNnr5C+S5zMyeMrM3zOx1M7s3BpkSZvaimf016izNQP06K0dsOhbHfoE6tlzq2Fk51LFzaETHIh3cmdllwGeAo1HmqNkHXOvu24G3gQcbuXAzSwA/Bm4HtgFfMLNtjcywgDJwv7tvAz4KfDUGme4F3ow4Q1NQv+aKYcfi2C9Qx86bOjaXOnbeVrxjUW+5+wHwDSDyA//c/Ql3L9euPgdsbHCEXcBhd3/H3YvA74E7GpxhDnc/7u4v1H4fIfxn3BBVHjPbCHwO+HlUGZqM+jVXrDoWt36BOnYB1LG51LFzaFTHIhvcmdkdwDF3fzmqDEv4CvB4g5e5Aeifdf1dIv4nnM3MrgBuAP4VYYwfEr6QViPM0BTUrwXFtmMx6ReoY+dNHVuQOnZuDelYciUf3Mz+Dqxf4KaHgG8Rbs5umKXyuPtfavd5iHBT7sONzBZnZpYHHgXuc/czEWXYA5x093+b2SeiyBA36ldriEO/ajnUsXnUsdawGju2ooM7d//UQvPN7DrgSuBlM4Nw8/ELZrbL3U80Os+sXPcAe4DbvPHniDkGXDbr+sbavEiZWYqwFA+7+58ijHIz8Hkz+yyQATrN7Lfu/qUIM0VK/Vq22HUsRv0Cdews6tiyqWNLa1jHYnGeOzP7L/Bhd4/si33NbDfwfeBWdx+IYPlJwoNgbyMsw/PAF9399UZnmZXJgF8BQ+5+X1Q55qu94/m6u++JOkszUL+mM8SqY3HtF6hjy6WOTWdQx87TSncs6g9UxMmPgA5gn5m9ZGY/aeTCawfCfg34G+FBn3+IcmBXczPwZeCTtefkpdo7DpHlirRfEMuOqV9ST+rY2VZtx2Kx5U5ERERE6kNb7kRERERaiAZ3IiIiIi1EgzsRERGRFqLBnYiIiEgL0eBOREREpIVocCciIiLSQjS4ExEREWkhGtyJiIiItJD/A60aB2Sl8+APAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 792x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Confidence levels and tails\n",
    "\n",
    "fig = plt.figure(figsize=(11,3))     # Figure size is (width,height) in inches\n",
    "_ = fig.suptitle('Tail probabilities for standard normal distribution, 95% level of confidence')\n",
    "\n",
    "alpha = 0.05\n",
    "areas = [alpha,alpha/2,1-alpha]\n",
    "zc = stats.norm.ppf(areas)\n",
    "\n",
    "xl = 10\n",
    "x = np.linspace(-xl,xl,1000)     # Generates 100 evenly spaced numbers in the interval [-xl,xl]\n",
    "xd ={'Lower':(-xl,zc[0]),'2 tailed':(-xl,zc[1],-zc[1],xl), 'Upper':(zc[2],xl)}\n",
    "\n",
    "for i,key in enumerate(xd):\n",
    "    ax = fig.add_subplot(1,3,i+1)\n",
    "    _ = ax.axes.get_yaxis().set_visible(False)\n",
    "    _ = ax.axes.set_xlim(-5,5)\n",
    "    _ = ax.plot(x,stats.norm.pdf(x),label=key)\n",
    "    _ = ax.legend(loc='upper right')\n",
    "    \n",
    "    x1 = np.linspace(xd[key][0],xd[key][1],20)\n",
    "    _=ax.fill_between(x1,0,stats.norm.pdf(x1),color='r')\n",
    "    if (key == \"2 tailed\"):\n",
    "        x2 = np.linspace(xd[key][2],xd[key][3],20)\n",
    "        _=ax.fill_between(x2,0,stats.norm.pdf(x2),color='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def decision_rule(alpha,p):\n",
    "    if (p < alpha):\n",
    "        decision = \"Reject H0\"\n",
    "    else:\n",
    "        decision = \"Fail to reject H0\"\n",
    "    return(decision)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Population mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lower-tailed test:  1.8498368520005377 0.9513118788634078    Decision:  Fail to reject H0\n",
      "Upper-tailed test:  1.8498368520005377 0.048688121136592155     Decision:  Reject H0\n",
      "Two-tailed test:  1.8498368520005377 0.09737624227318431     Decision:  Fail to reject H0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Ttest_1sampResult(statistic=1.8498368520005373, pvalue=0.09737624227318435)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Hypothesis testing for population mean\n",
    "\n",
    "np.random.seed(12345)\n",
    "n=10\n",
    "x = stats.norm.rvs(size=n)  # Standard normal random variates\n",
    "alpha = 0.05\n",
    "\n",
    "xbar = np.mean(x)\n",
    "s = np.std(x, ddof=1)    # ddof=1 uses (n-1) for the denominator, so gives us sample sd\n",
    "\n",
    "mu = 0.0                # H0: mu = 0               \n",
    "T =  (xbar-mu)/ (s/np.sqrt(n))    \n",
    "\n",
    "# Lower-tailed test\n",
    "pl = stats.t.cdf(T,n-1)\n",
    "print(\"Lower-tailed test: \",T,pl,\"   Decision: \", decision_rule(alpha,pl))\n",
    "\n",
    "# Upper-tailed test\n",
    "pu = 1-stats.t.cdf(T,n-1)\n",
    "print(\"Upper-tailed test: \",T,pu,\"    Decision: \", decision_rule(alpha,pu))\n",
    "\n",
    "# Two-tailed test\n",
    "p2 = 2*(min(pl,pu)) # Two-tailed test, so twice the tail probability. Note that degrees of freedom = n-1\n",
    "print(\"Two-tailed test: \",T,p2,\"    Decision: \", decision_rule(alpha,p2))\n",
    "\n",
    "stats.ttest_1samp(x,0)     # Default behavior of ttest_1samp: H0: mu = 0, Ha: mu not equal to zero\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Population variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample variance =  0.7155771977102666\n",
      "H0: sigma^2 =  0.36\n",
      "Level of significance =  0.05\n",
      "Chi-squared distribution,  9 degrees of freedom\n",
      "\n",
      "Lower-tailed test:  17.889429942756667 0.9635217652248893    Decision:  Fail to reject H0\n",
      "Upper-tailed test:  17.889429942756667 0.036478234775110696    Decision:  Reject H0\n",
      "Two-tailed test:  17.889429942756667 0.07295646955022139    Decision:  Fail to reject H0\n"
     ]
    }
   ],
   "source": [
    "# Hypothesis testing for population variance\n",
    "\n",
    "np.random.seed(12345)\n",
    "n=10\n",
    "x = stats.norm.rvs(size=n)  # Standard normal random variates\n",
    "alpha = 0.05\n",
    "\n",
    "xbar = np.mean(x)\n",
    "s = np.std(x, ddof=1)    # ddof=1 uses (n-1) for the denominator, so gives us sample sd\n",
    "print(\"Sample variance = \",s**2)\n",
    "\n",
    "sigma = 0.6                        # H0\n",
    "print(\"H0: sigma^2 = \",sigma**2)\n",
    "print(\"Level of significance = \",alpha)b\n",
    "print(\"Chi-squared distribution, \", n-1, \"degrees of freedom\\n\")\n",
    "\n",
    "X2 = (n-1)*(s**2)/(sigma**2)      # Test statistic\n",
    "\n",
    "# Lower-tailed test\n",
    "pl = stats.chi2.cdf(X2,n-1)\n",
    "print(\"Lower-tailed test: \",X2,pl,\"   Decision: \", decision_rule(alpha,pl))\n",
    "\n",
    "# Upper-tailed test\n",
    "pu = 1-stats.chi2.cdf(X2,n-1)\n",
    "print(\"Upper-tailed test: \",X2,pu,\"   Decision: \", decision_rule(alpha,pu))\n",
    "\n",
    "# Two-tailed test\n",
    "p2 = 2*(min(pl,pu))\n",
    "print(\"Two-tailed test: \",X2,p2,\"   Decision: \", decision_rule(alpha,p2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Population proportion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H0: pi =  0.1\n",
      "Significance level =  0.05\n",
      "Sampling distribution: Binomial(n,pi) with n =  10 \n",
      "\n",
      "Lower-tailed test:  3 0.9872048016    Decision:  Fail to reject H0\n",
      "Upper-tailed test:  3 0.012795198399999985    Decision:  Reject H0\n"
     ]
    }
   ],
   "source": [
    "# Hypothesis testing for population proportion\n",
    "n = 10\n",
    "p = 0.3\n",
    "c = int(n*p)        # Sample count of successes\n",
    "\n",
    "pi = 0.1            # H0\n",
    "alpha = 0.05\n",
    "print(\"H0: pi = \",pi)\n",
    "print(\"Significance level = \", alpha)\n",
    "print(\"Sampling distribution: Binomial(n,pi) with n = \",n,\"\\n\")\n",
    "# Lower-tailed test\n",
    "\n",
    "pl = stats.binom.cdf(c,n,pi)   # Prob (counts <= sample counts)\n",
    "print(\"Lower-tailed test: \",c,pl,\"   Decision: \", decision_rule(alpha,pl))\n",
    "\n",
    "# Upper-tailed test\n",
    "\n",
    "pu = 1-stats.binom.cdf(c,n,pi)   # Prob (counts >= sample counts)\n",
    "print(\"Upper-tailed test: \",c,pu,\"   Decision: \", decision_rule(alpha,pu))"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}