ee046202_tutorial_03_classic_statistics_hypothesis_testing_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# <img src=\"https://img.icons8.com/dusk/64/000000/artificial-intelligence.png\" style=\"height:50px;display:inline\"> EE 046202 - Technion - Unsupervised Learning & Data Analysis\n",
    "---\n",
    "\n",
    "#### <a href=\"https://taldatech.github.io\"> Tal Daniel</a>\n",
    "\n",
    "## Tutorial 03 - Classical Methods in Statistical Inference - Hypothesis Testing 1\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/bubbles/50/000000/checklist.png\" style=\"height:50px;display:inline\"> Agenda\n",
    "---\n",
    "\n",
    "* [Hypothesis Testing](#-Hypothesis-Testing)\n",
    "* [Hypothesis Testing Steps](#-Hypothesis-Testing-Steps)\n",
    "* [Error Types](#-Error-Types)\n",
    "* [Example - Body Weight - Z-Statistic](#-Example---Body-Weight---Hypothesis-Testing-for-the-Mean)\n",
    "* [The Central Limit Theorem (CLT)](#-Reminder:-The-Central-Limit-Theorem-(CLT))\n",
    "* [Relation to Confidence Intervals](#-Relation-to-Confidence-Intervals)\n",
    "* [Hypothesis Testing for the Mean Summary](#-Hypothesis-Testing-for-the-Mean-Summary)\n",
    "* [Recommended Videos](#-Recommended-Videos)\n",
    "* [Credits](#-Credits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# imports for the tutorial\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/bubbles/50/000000/test-tube.png\" style=\"height:50px;display:inline\"> Hypothesis Testing\n",
    "---\n",
    "Let's begin with an example: consider a radar system that uses radio waves to detect aircrafts. The system receives a signal, and based on the received signal, it needs to decide whether an aircraft is present or not. We formulate two hypotheses:\n",
    "1. $H_0$: No aircraft is present.\n",
    "2. $H_1$: An aircraft is present.\n",
    "\n",
    "$H_0$ is called the **null hypothesis** (also: the *default hypothesis*) and $H_1$ is called the **alternative hypothesis**. We initially assume that $H_0$ is **true** and based on the observed data we need to decide whether or not to accept $H_1$ or reject it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/dusk/64/000000/workflow.png\" style=\"height:50px;display:inline\"> Hypothesis Testing Steps\n",
    "---\n",
    "#### Step 1\n",
    "---\n",
    "**Null & Alternative Hypotheses** -  Formulate the null hypothesis $H_0$: $\\theta \\in \\Theta_0$ (that the observations are the result of pure chance) and the alternative hypothesis $H_1$: $\\theta \\in \\Theta_1$ (that the observations show a real effect combined with a component of chance variation).\n",
    "\n",
    "#### Step 2\n",
    "---\n",
    "**Test Statistic** - Identify a test statistic that can be used to assess the truth of the null hypothesis. It is a value computed from sample data. The test statistic is used to assess the strength of evidence in support of a null hypothesis.\n",
    "   * A **statistic** is a real-valued function of the data. For example, the sample mean: $W(X_1, X_2,..., X_n)= \\frac{X_1 + X_2 + ... + X_n}{n}$ is a statistic.\n",
    "   * A **test satistic** is a statistic on we which we build our test.\n",
    "   * **Acceptence Region** $A$ - A *set* $A \\subset \\mathbb{R}$ is defined to be the set of all possible values of the test statistic for which we accept $H_0$.\n",
    "   * **Rejection Region** $R$ - A *set*  $R = \\mathbb{R} - A$ is defined to be the set of all possible values of the test statistic for which we reject $H_0$ and accept $H_1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Step 3\n",
    "---\n",
    "**P-value & Interpretation** - Compute the P-value, which is the probability that a test statistic, at least as significant as the one observed, would be obtained assuming that the null hypothesis were true. The smaller the P-value, the stronger the evidence **against** the null hypothesis.\n",
    "\n",
    "#### Step 4\n",
    "---\n",
    "**Significance Level** - Compare the p-value to an acceptable significance value $\\alpha$ (sometimes called an $\\alpha$ value, a probability threshold below which the null hypothesis will be rejected. Common values are 5% and 1%.). If $p \\leq \\alpha$ (the observed effect is statistically significant), the null hypothesis is ruled out, and the alternative hypothesis is valid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/ios/100/000000/battle.png\" style=\"height:30px;display:inline\"> The Two Competing Theories\n",
    "---\n",
    "* **NULL** Hypothesis - $H_0$ - any observed deviation from what we expect to see is due to chance variability.\n",
    "* **ALTERNATIVE** Hypothesis - $H_a$ - 'claim', or a theory you wish to test (the reason for the observed statistic).\n",
    "\n",
    "$H_0$ is assumed **true** until enough evidence goes against it (we then refute it and believe the alternative $H_a$).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/bubbles/50/000000/close-window.png\" style=\"height:50px;display:inline\"> Error Types\n",
    "---\n",
    "<img src=\"./assets/tut_01_error_types.jpg\" style=\"height:300px\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* **Type I Error (False Positive)** -  the incorrect rejection of a true null hypothesis. Usually a type I error leads one to conclude that a supposed effect or relationship exists when in fact it doesn't. \n",
    "    * For example, a test that shows a patient to have a disease when in fact the patient does not have the disease, a fire alarm going on indicating a fire when in fact there is no fire, or an experiment indicating that a medical treatment should cure a disease when in fact it does not.\n",
    "    * The chance of **rejecting the null hypthesis $H_0$, when it is TRUE**, denoted by $\\alpha$\n",
    "    * $\\rightarrow$ the chance of **accepting the null hypthesis $H_0$, when it is TRUE** is $1 - \\alpha$\n",
    "    * Formally:\n",
    "        * Denote a *test statistic* as $W$\n",
    "        * $P(\\text{Type 1 Error} | \\theta) = P(\\text{Reject } H_0|\\theta) = P(W \\in R | \\theta), \\theta \\in \\Theta_0$\n",
    "            * If $P(\\text{Type 1 Error} | \\theta) \\leq \\alpha, \\forall \\theta \\in \\Theta_0$, we say that the test has **significance level** $\\alpha$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* **Type II Error (False Negative)** - the failure to reject a false null hypothesis.\n",
    "    * For example, a blood test failing to detect the disease it was designed to detect, in a patient who really has the disease; a fire breaking out and the fire alarm does not ring; or a clinical trial of a medical treatment failing to show that the treatment works when really it does.\n",
    "    * The chance of **not rejecting the null hypothesis $H_0$, when it is FALSE**, denoted by $\\beta$\n",
    "    * $\\rightarrow$ the chance of **rejecting the null hypthesis $H_0$, when it is FALSE** is $1 - \\beta$ (also called **power**)\n",
    "    * Since the alternative hypothesis, $H_1$, usually contains more than one value of $\\theta$, the probability of type II error is usually a **function of $\\theta$**, and denoted by $\\beta$.\n",
    "    * Formally: $\\beta(\\theta) = P(\\text{Accept } H_0 | \\theta), \\text{ for } \\theta \\in \\Theta_1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/classroom.png\" style=\"height:20px;display:inline\"> Example - Error Types\n",
    "---\n",
    "* Hypothesis: \"A patient's symptoms improve after treatment A more rapidly than after a placebo treatment.\"\n",
    "\n",
    "* Null hypothesis ($H_0$): \"A patient's symptoms after treatment A are indistinguishable from a placebo.\"\n",
    "\n",
    "* A Type I error would falsely indicate that treatment A is more effective than the placebo, whereas a Type II error would be a failure to demonstrate that treatment A is more effective than placebo even though it actually is more effective."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/doodle/96/000000/body-measurements.png\" style=\"height:50px;display:inline\"> Example - Body Weight - Hypothesis Testing for the Mean\n",
    "---\n",
    "The following example will be used to demonstrate the statistic process:\n",
    "\n",
    "In the 1970s, 20–29 year old men in the U.S. had a mean $\\mu$ body weight of 170 pounds (77 kg). Standard deviation $\\sigma$ was 40 pounds (18 kg). We test whether mean body weight in the population is **bigger** now."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 1- Null & Alternative Hypotheses \n",
    "---\n",
    "* Under the **null hypothesis** there is no difference in the mean body weight between then and now, in which case $\\mu$ would still equal 170 pounds: $$H_0: \\mu=170$$\n",
    "* Under the **alternative hypothesis**, the mean weight has increased: $$H_a: \\mu > 170$$\n",
    "    * This statement of the alternative hypothesis is **one-sided**. That is, it looks only for values larger than stated under the null hypothesis.\n",
    "    * There is another way to state the alternative hypothesis. We could state it in a **“two-sided”** manner, looking for values that are either higher- or lower-than expected. For the current illustrative example, the two-sided alternative is $H_a: \\mu \\neq 170$. Although for the current illustrative example, this seems unnecessary, two-sided alternative offers several advantages and are much more common in practice."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 2- Test Statistic (TS)\n",
    "---\n",
    "* It is a measure of how far the observed data is from what is expected under the null hypothesis $H_0$.\n",
    "    * Compute the value of a test statistic (TS) from the data.\n",
    "* The particular TS computed depends on the tested parameter.\n",
    "    * For example, to test the population mean, the TS is the sample mean (or standardized sample mean), if the variance is known.\n",
    "    * It is very similar to the process we did in point estimation, for choosing the correct estimator.\n",
    "* The null hypothesis $H_0$ is rejected if the TS falls in a user-specified rejection region.\n",
    "* Different hypothesis tests use different test statistics based on the probability model assumed in the null hypothesis. Common tests and their test statistics include:\n",
    "\n",
    "| Hypothesis Test | Test Statistic |\n",
    "| ---------- | ---------- |\n",
    "|<center> Z-test </center>| <center> Z-statistic </center>|\n",
    "|<center> t-tests </center>|<center> t-statistic</center>|\n",
    "|<center> ANOVA</center>| <center> F-statistic </center>|\n",
    "|<center> Chi-square tests</center>| <center> Chi-square statistic </center>|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/classroom.png\" style=\"height:20px;display:inline\"> Example - TS - Z-statistic\n",
    "---\n",
    "The Z-statistic has the standard normal distribution under the null hypothesis. It is a **mean** test when $\\sigma$ is known. We will use this statistic to test the problem: $$z_{stat} = \\frac{\\overline{x} - \\mu_0}{\\sigma_{\\overline{x}}} $$ \n",
    "\n",
    "Assumptions:\n",
    "* $\\mu_0$ is the **population mean** assuming $H_0$ is **true**\n",
    "* $\\overline{x} = \\frac{x_1 + x_2 +.. +x_n}{n}$ is the sample mean.\n",
    "* $\\sigma_{\\overline{x}} = \\frac{\\sigma}{\\sqrt{n}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Solve for the \"Body Weight\" problem:\n",
    "* $\\mu_0 =170$\n",
    "* $\\sigma = 40$\n",
    "* We'll take sample size of $n=64 $ samples, $\\rightarrow \\sqrt{n} = 8$\n",
    "* $\\sigma_{\\overline{x}} = \\frac{40}{8} = 5 $\n",
    "    * If the **variance is unknown** (the variance of $X_i$): we use the **sample standard deviation** $S$ instead of $\\sigma$ (unbiased) $$S = \\sqrt{\\frac{1}{n-1}\\sum_{k=1}^n(X_k - \\overline{X})^2} = \\sqrt{\\frac{1}{n-1}\\big(\\sum_{k=1}^nX_k^2 -n\\overline{X}^2\\big)}$$ \n",
    "    * Note that you first estimate the variance for the **whole** population, and then still devide by $\\sqrt{n}$ to get the sample STD for $z_{stat}$. \n",
    "    * <a href=\"https://www.probabilitycourse.com/chapter8/8_4_3_hypothesis_testing_for_mean.php\">Read more (examples + exercises)</a>.\n",
    "\n",
    "Now, let's assume we found a sample mean of **173**, then:\n",
    "* $z_{stat} = \\frac{\\overline{x} - \\mu_0}{\\sigma_{\\overline{x}}} = \\frac{173 - 170}{5} = 0.6$\n",
    "\n",
    "Now, let's assume we found a sample mean of **185**, then:\n",
    "* $z_{stat} = \\frac{\\overline{x} - \\mu_0}{\\sigma_{\\overline{x}}} = \\frac{185 - 170}{5} = 3$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/dusk/64/000000/alarm.png\" style=\"height:30px;display:inline\"> Reminder: The Central Limit Theorem (CLT)\n",
    "---\n",
    "* The CLT states that given a sufficiently large sample size from a population with a finite level of variance, the mean of all samples from the same population will be approximately equal to the mean of the population. \n",
    "* When $n$ is large, the distribution of the **sample means** will approach a normal distribution. More formally:\n",
    "\n",
    "If $X_1, X_2, ..., X_n$ is a random sample of size $n$ taken from a population with mean $\\mu$ and variance $\\sigma^2$, and if $\\overline{X}$ is the sample mean, the limiting form of the distribution: $$ Z = \\frac{\\overline{X} - \\mu}{\\frac{\\sigma}{\\sqrt{n}}} $$ as $n \\rightarrow \\infty $, is the **standard normal distribution**\n",
    "\n",
    "<a href=\"http://onlinestatbook.com/stat_sim/sampling_dist/\">CLT DEMO</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Reasoning Behind $z_{stat}$\n",
    "---\n",
    "Sampling distribution of $\\overline{x}$ under $H_0$: $$ \\overline{x} \\sim N(170, 5) $$\n",
    "<img src=\"./assets/tut_01_z_stat.jpg\" style=\"height:200px\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 3 - P-value & Interpretation\n",
    "---\n",
    "* All hypothesis tests ultimately use a **p-value** to weigh the strength of the evidence (what the data are telling you about the population). The p-value is a number between 0 and 1, and is **the probability of the observed test statistic (or one more extreme) when $H_0$ is true**\n",
    "    * **P-value** is the *lowest* significance level $\\alpha$ that results in rejecting the null hypothesis. \n",
    "* It corresponds to the **Area Under the Curve (AUC)** in the tail of the Standard Normal Distribution beyond the $z_{stat}$\n",
    "* Converting Z-statistics to *P-value*:  $$\\textit{For } H_1: \\mu > \\mu_0 \\rightarrow P = Pr(Z > z_{stat}) = \\textit{right-tail beyond } z_{stat} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"./assets/tut_01_p_value.png\" style=\"height:300px\">\n",
    "\n",
    "\n",
    "* <a href=\"https://courses.lumenlearning.com/wmopen-concepts-statistics/chapter/hypothesis-test-for-a-population-proportion-2-of-3/\">Image Source</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# let's see for the body weight problem\n",
    "x = np.linspace(140, 200, 64)\n",
    "mu = 170  # H_0 is true!\n",
    "sigma = 5  # calcualted for 64 samples\n",
    "f_x = (1 / np.sqrt(2 * np.pi * sigma ** 2)) * np.exp(- (x - mu) ** 2 / (2 * sigma ** 2))\n",
    "x_normed = (x - mu) / sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p-val (AUC) = 0.298\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAGDCAYAAAAYtQWTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeVyV95n//9fFJgi4oRh33OIel+CeJphkEm3SmqTJZGmb2mlrM006mfb77Yydpk3a36Ttd2Y6bTOTNrXTLdM0SZNmbbOrRBMX0LiiqKCguCGgCCKyfX5/nANFBDnggfsczvv5ePCQc+7Pue/3+XDEy/u+7vs25xwiIiIiElxRXgcQERER6YlUZImIiIh0ARVZIiIiIl1ARZaIiIhIF1CRJSIiItIFVGSJiIiIdAEVWSI9jJn91sz+1esc7TEzZ2bjvM7RyMwSzOx1Mys3sxe8ziMi4U9Flkg3MLMCMztnZpVmdsLMfmNmSa2Mm29mZ80suZVlW83soe5JHBr8BWOdmQ1t5fl/bfFcmr9wi2n23H1mttk/78fM7E0zu6aNzd0JDAZSnHN3Bf3NdIKZZZhZUQdf02XFq5l9zcyO+wvRX5tZr0uMjTazfzWzo2ZW4f/89uuKXCKhSkWWSPf5hHMuCZgFzAYeaTnAObcBKAI+1fx5M5sKTAae7YacIcHMEvHNQznw6U68/uvAT4Dv4yueRgI/A5a28ZJRwD7nXF0nthXT/qjwZmY3AyuAG4A0YAzw3Uu85LvAAmA+0Af4LFDdtSlFQouKLJFu5pw7ArwJTG1jyO+A+1s8dz/wF+dcKYCZvdBsj8JaM5vS2orMbJmZfdDiuaY9HWbWy8z+w8wO+fewPWVmCW2sa6yZrTazUjMrMbNnmu+Z8O+t+79mtsOf63kzi2+2/Bv+vUlHzezvLj1LgK/AOg18D/hcAOObZ+3rf92DzrmXnHNnnXO1zrnXnXPfaGX8d4HvAHf793p9wcyizOwRMys0s2Ize9q/3uZ7zb5gZoeA1a2s83X/uhq/GsxsWRt5P25mu/17fI745zER3+dkaLN1DDWzOWa2wcxO++fzv80szr+etf5VbvePv7sj89aOzwG/cs7lOOdOAf8f0Nb76Q/8I/Al51yh89nlnFORJRFFRZZINzOzEcDHga1tDPlf4GNmNtI/Pgq4D3i62Zg3gfFAKvAR8Ewn4/w/4EpgBjAOGIav2Gg1OvADYCgwCRgBPNZizN8Ci4HRwFX4/xE2s8XA/wX+xp/7xgCyfQ7fnrvngIlmNiuwtwT49p7EAy8HMtg59yi+PV7PO+eSnHO/8mdfBizCt9cmCfjvFi+9Dt9c3NzKOj/hX1cSvkORx4FVbUT4FfBl51wyvuJ7tXPuLLAEONq4HufcUaAe+Bow0P8+bwC+4t/mtf71TfePf77lhszsGn+B1tZXW4dTpwDbmz3eDgw2s5RWxk4D6oA7/f8Z2GdmD7axXpEeS0WWSPd5xcxOAx8A7+P7R/0izrnD/uWf8T91A76C4S/NxvzaOVfhnDuPr9CZ3riXJVBmZsCXgK8558qccxX+TPe0kSvPOfeuc+68c+4k8J/4iozmnnDOHXXOlQGv4yvewFd8/ca/N+MsFxdnLbONxFfc/ME5dwJfcdKRvVkpQElnDv0182ngP51zB5xzlcA3gXtaHBp8zL+X7FxbKzGzK/EVyHf7f7atqQUmm1kf59wp59xHba3PObfFObfROVfnnCsAfsHFP4c2Oec+cM71u8TXB228NAnfodtGjd9f1D8IDAf64ivgR+MrMh8zs78JNKdIT6AiS6T73Ob/R2yUc+4rzrlzZvaxZoeCcpqNbX7I8LP4io1aaGoo/qGZ5ZvZGaDAP25gB/MMAnoDWxr3YgBv+Z+/iJmlmtlz/sNZZ4Dft7LN482+r8L3DzP49n41LzAK28n2WWCPc26b//EzwH1mFut/XAfEtnhNLNDg/yoFBl5mr9TQFjkLgRh8/V2N2iqagKbDlq8C33bOrfM/9y/NfuZP+Yd+Ct/ezUIze9/M5l9inVea2Z/9e4jO4CuMO/qz74xKfL1VjRq/r2hlbGPR+T3n3Dnn3A58eyQ/3oX5REKOiiwRDznn1jU7FNS8r+olYJiZLQLu4MJDhffha96+Ed/egjT/89bKJs7iK6R8A8yuaLasBN8/hlOa7cXo6z+81ZofAA64yjnXB9+etta22Zpj+A4vNhrZzvj7gTH+QuI4vr1mA/EdPgM4xF/fd6PRwGHnXAOwAV+T9W0B5mvNUXzN8M0z1wEnmj3n2nqx/zDvH4A1zrlfNL3Aue83+5k/4H8u2zm3FN/h31eAP15i/T8HcoHx/p/DvxD4z4EWhX1rXx9r46U5wPRmj6cDJxr7BFvYcYn8IhFDRZZICPIfUnsR+A1Q6Jzb3GxxMnAe396a3rRx2NFvOzDFzGb4m9Afa7aNBuCXwI/NLBXAzIaZ7yyy1iTj25tx2syGARc1kF/CH4FlZjbZzHoDj7Y10L8XZywwB9/hxhn4+pT+wF8PGf4JuMXMbvLv2RuK72zN5/zvrRxfb9mTZnabmfU2s1gzW2Jm/xZg5meBr5nZaPNdbqOxZyvQQ5CPA4nAw5caZGZxZvZpM+vr31t5Bl/fFfgKupQWh4KT/WMqzWwi8PctVnkCXw9Zq1oU9q19rWvjpU8DX/D/DPvjm+/ftrGNfGAd8C3znVwxCbgb+PMlpkKkx1GRJRK6fodvT8rTLZ5/Gt+hqyPAbmBjWytwzu3Dd5bde8B+fP1gzf0zkAds9B96eg+Y0Mbqvovv8hPl+PrDXgr0jTjn3sR3OYXV/u1ddDZeM58DXnXO7XTOHW/8An4K3GpmA5xzOcC9+PauleHbc7WJZpcUcM79J/B1fMXASXyH9h7Ct6coEL/GdxLCWuAgvj1jXw3wtfjzzQNONdtL1NalKD4LFPh/Bg/g78dzzuXiK/YO+A/pDsV3AsF9+A7T/RJo2dz+GPA7//i/7UDeS3LOvQX8G7AG3+evkGbFsvmuQfYvzV5yL77Pbym+z8u3nXNtNf6L9EjmnPbmioiIiASb9mSJiIiIdAEVWSIiIiJdQEWWiIiISBdQkSUiIiLSBVRkiYiIiHSBkLxz/MCBA11aWlqXbuPs2bMkJiZ26TYiieYz+DSnwac5DS7NZ/BpToOru+Zzy5YtJc65i+6WEZJFVlpaGps3b25/4GXIzMwkIyOjS7cRSTSfwac5DT7NaXBpPoNPcxpc3TWfZtbqrcJ0uFBERESkC6jIEhEREekCKrJEREREukBI9mSJiEh4q62tpaioiOrqaq+jhJW+ffuyZ88er2P0GMGez/j4eIYPH05sbGxA41VkiYhI0BUVFZGcnExaWhpm5nWcsFFRUUFycrLXMXqMYM6nc47S0lKKiooYPXp0QK/R4UIREQm66upqUlJSVGBJj2FmpKSkdGjvrIosERHpEiqwpKfp6Gc6oCLLzBab2V4zyzOzFa0sNzN7wr98h5nNarbsa2aWY2a7zOxZM4vvUEIREZFOePzxx5kyZQpXXXUVM2bMYNOmTQD85Cc/oaqqKmjbSUtLo6SkpNOvz8zM5NZbb231eTPj9ddfb3ru1ltvJTMzs9Pb6oy23l9aWhrTpk1j2rRpTJ48mUceeYTz588DcPToUe68884213n69Gl+9rOfXXK7CxYsANqen0t55ZVX2L17d9Pj73znO7z33nsdWkcwtFtkmVk08CSwBJgM3Gtmk1sMWwKM938tB37uf+0w4B+AdOfcVCAauCdo6UVERFqxYcMG/vznP/PRRx+xY8cO3nvvPUaMGAEEv8jqqPr6+oDHDh8+nMcff7xbttUZa9asYefOnWRlZXHgwAGWL18OwNChQ3nxxRfbfN2liqzGzOvXr+90rpZF1ve+9z1uvPHGTq+vswLZkzUHyHPOHXDO1QDPAUtbjFkKPO18NgL9zGyIf1kMkGBmMUBv4GiQsouIiLTq2LFjDBw4kF69egEwcOBAhg4dyhNPPMHRo0dZtGgRixYtAuDv//7vSU9PZ8qUKTz66KNN60hLS+PRRx9l1qxZTJs2jdzcXABKS0u56aabmDlzJl/+8pdxzjW95rbbbuPqq69mypQprFy5sun5pKQkvvOd7zB37lw2bNjAW2+9xcSJE7nmmmt46aWX2nwf06dPp2/fvrz77rsXLVu1ahUzZ85k2rRp/N3f/V3TXqS0tDS+973vcc011/DCCy+QkZHB1772Na699lomTZpEdnY2d9xxB+PHj+eRRx5pN3sgkpKSeOqpp3jllVcoKyujoKCAqVOnApCTk8OcOXOYMWMGV111Ffv372fFihXk5+czY8YMvvGNb5CZmcmiRYu47777mDZtWtM6G505c4bbb7+dyZMn88ADD9DQ0HDRmBdffJFly5axfv16XnvtNb7xjW+wcOFC8vPzWbZsWVPRd6l5a+3nfTkCObtwGHC42eMiYG4AY4Y55zab2X8Ah4BzwDvOuXda24iZLce3F4zBgwd3+e7QysrKbt/l2pNpPoNPcxp8mtPgutR89u3bl4qKCgD+3zv55J6oDOq2Jw5O4p9vGtvm8vnz5/PYY48xbtw4MjIy+NSnPsU111zD5z//eX70ox/x+uuvk5KSQkVFBStWrGDAgAHU19fziU98gsWLFzN16lSccyQlJfH+++/zy1/+kh/84Af893//N9/61reYPXs2K1as4K233mLlypVUVlbSq1cvfvrTnzJgwADOnTtHRkYGN910EykpKZw9e5axY8fy3nvvUV1dzcyZM3n99dcZO3Ysy5Yto66ujoqKCurr65vmraqqirq6Or7+9a/z3e9+l3nz5lFXV0dVVRUnT57kc5/7HK+99hrjx49n+fLl/PjHP+bBBx/EOYeZ8eabbwLw5JNPAvCXv/yFn/3sZ3zyk59k7dq19O/fn+nTp/PFL36RlJSUNrM755reX3MtnzczRo0axbZt20hNTaWhoYGKigqeeOIJli9fzt13301NTQ319fU88sgj7Nixg3Xr1gGwbt06srKy2LhxI2lpaU1zUFFRQVVVFVlZWWRlZTFy5EjuuOMOnnnmGW677bamMQDnzp2jtraWadOmsWTJEhYvXswnPvEJoqOjqa2t5dy5c+3OW2s/75aqq6sD/j0SSJHVWpeXC2SMmfXHt5drNHAaeMHMPuOc+/1Fg51bCawESE9Pd119ryHdHyq4NJ/BpzntPOcc24vKOVVVc8HzB3bsYNp4X7dDQmw0c9IGEBWl5uzOutRndM+ePU2nzsfGxRIdHR3UbcfGxV7y1Pzk5GS2bt3KunXrWLNmDZ///Of54Q9/yLJlyzAzkpKSml7/zDPPsHLlSurq6jh27BiFhYXMnz8fM+O+++4jOTmZhQsX8sYbb5CcnMzGjRt56aWXSE5O5q677qJ///5N6/vRj37Eyy+/DMCRI0c4fvw4aWlpREdH85nPfIbo6Gjy8/MZM2YMM2fOBGDZsmWsXLmS5OTkCy450Lt3b2JiYrj55pv5/ve/z7Zt24iJiaF3794cPXqUMWPGMGuWrwX6i1/8Ik8++SQrVqzAzLj//vub1hMdHc2dd95JcnIys2fPZurUqYwfPx6AsWPHcvr0adLS0trM3nK+GrX2vJmRmJhIUlISUVFRJCcnc9111/H4449TWlratAetqqqqaXnje50zZ07TXqzmP8fGZVdddRUAn/nMZ9iyZQuf/exnm8YAJCQkEBvr+1zExsaSkJBAdHT0BY/bm7fWft4txcfHN/3s2hNIkVUEjGj2eDgXH/Jra8yNwEHn3EkAM3sJWABcVGSJiATL0xsKefS1nNYXbslu+nbZgjQe++SUbkoVuR79hDdzHB0dTUZGBhkZGUybNo3f/e53LFu27IIxBw8e5D/+4z/Izs6mf//+LFu27IJT9Bv30kRHR1NXV9f0fGtnmWVmZvLee++xYcMGevfuTUZGRtO64uPjLyg0O3qW2re+9S0ef/xxYmJ8/2w3P0TZmsTExAseN76PqKioC/ZIRUVFUVdXd8nsgaqoqKCgoIArr7yS8vLypufvu+8+5s6dy1/+8hduvvlm/ud//ocxY8a0m7m5lvPV+Lj584HkbW/e2vp5d1YgPVnZwHgzG21mcfga119rMeY14H7/WYbzgHLn3DF8hwnnmVlv883EDYAuZSsiXWbb4dP86192c/3EVF7+yoILvr49L77p+8/MG8lv1xfwlx3HvI4sXWDv3r3s37+/6fG2bdsYNWoUQNMeI/D1+iQmJtK3b19OnDjRdIjtUq699lqeeeYZAN58801OnToFQHl5Of3796d3797k5uaycePGVl8/ceJEDh48SH5+PgDPPvtsu9u86aabOHXqFNu3b29aR0FBAXl5eQD87//+L9ddd12762lLoNnbUllZyVe+8hVuu+02+vfvf8GyAwcOMGbMGP7hH/6BT37yk+zYseOCn0EgsrKyOHjwIA0NDTz//PNcc801gK+9aM+ePTQ0NDTthQPaXH+w56097e7Jcs7VmdlDwNv4zg78tXMux8we8C9/CngD+DiQB1QBn/cv22RmLwIfAXXAVvyHBEVEgu3U2RoefOYjUpPj+c+/nU6/3nEXLC8/EM3Mkb5/AKYO68vuo2f4pxe3M3FIMmMHJbW2SglTlZWVfPWrX+X06dPExMQwbty4pmbu5cuXs2TJEoYMGcKaNWuYOXMmU6ZMYcyYMSxcuLDddT/66KPce++9zJo1i+uuu46RI0cCsHjxYp566imuuuoqJkyYwLx581p9fXx8PCtXruSWW25h4MCBXHPNNezatavd7X7rW99i6dKlTev4zW9+w1133UVdXR2zZ8/mgQceCHR6LhJo9pYWLVqEc46GhgZuv/12vv3tb1805vnnn+f3v/89sbGxXHHFFXznO99hwIABLFy4kKlTp7JkyRJuueWWS25n/vz5rFixgp07d3Lttddy++23A/DDH/6QW2+9lREjRjB16lQqK329f/fccw9f+tKX+MlPfnLBiQXBnrf2WHu7zryQnp7uNm/e3KXbUL9LcGk+g09z2jENDY4v/C6bD/NKeeGB+Uwf0e+iMS3n9Ojpc9zyxDpSk+N55cGFJMQFt2+op2uvJ2vSpEndG6gH0G11gqsr5rO1z7aZbXHOpbccqyu+i0iP8PP381mz9ySP3Dqp1QKrNUP7JfCTe2ayr7iCb7/a/p4EEZGOUJElImFvQ34pP3pnL5+YPpTPzhvVodded+UgvrpoHC9uKeKP2Yfbf4GISIBUZIlIWCuuqOarz24lbWAiP7hjWqful/fwjVeycFwK3351F7uPnumClCISiVRkiUjYqqtv4B+e3Url+Vp+/umrSeoVyFVpLhYdZfzk7pn0TYjlwT98REV1bZCTikgkUpElImHrx+/tY+OBMv71tmlMuOLymlsHJffiv++bxaGyKv75TzvavZ6OiEh7VGSJSFhak1vMk2vyuTt9BHdePTwo65wzegDfuHkCb+w8zm/XFwRlnSISuVRkiUjYKT5Tzdf+uI1JQ/rw3aXBvZr48o+N4cZJqXz/jT3sOlLe/gskZEVHRzNjxgymTp3KXXfdRVVV1WWv87e//S0PPfRQENJJJFCRJSJh5/ebDlF+rpb/uncm8bHBvbZVVJTxo7tmkBAbzcq1B4K67oh2xRVgFryvK65od5MJCQls27aNXbt2ERcXx1NPPdUNb1Tkr1RkiUhYqatv4PnsQ1w7fhDjUrvmKu19e8dyx6zhvLXrOGVna9p/gbTvxAlP1/exj32s6VYqjRoaGkhLS+P06dNNz40bN44TJ07w+uuvM3fuXGbOnMmNN97IiVa2t2zZMl588cWmx0lJf/08/vu//zuzZ8/mqquu4tFHH+1QVuk5VGSJSFhZnVvMiTPnuW/uyC7dzqfnjqSmvoEXt+jaWeGurq6ON998k2nTpl3wfFRUFEuXLm26592mTZtIS0tj8ODBXHPNNWzcuJGtW7dyzz338G//9m8Bb++dd95h//79ZGVlsW3bNrZs2cLatWuD+p4kPKjIEpGw8oesQwzu04sbJqZ26XbGD05mdlp/ns06rDMNw9S5c+eYMWMG6enpjBw5ki984QsXjbn77rt5/vnnAXjuuee4++67ASgqKuLmm29m2rRp/Pu//zs5OTkBb/edd97hnXfeYebMmcyaNYvc3NwLblYtkUNFloiEjcNlVby/7yR3p48gJrrrf33dN3ckB0vOsiG/tMu3JcHX2JO1bds2/uu//ou4uDiefPJJZsyYwYwZMzh69Cjz588nLy+PkydP8sorr3DHHXcA8NWvfpWHHnqInTt38otf/ILq6uqL1h8TE0NDQwMAzjlqamqavv/mN7/ZtO28vLxWCzzp+VRkiUjYeD77MAbcPadrDxU2WjJ1CP16x/JM1qFu2Z50vQcffLCp+Bk6dChmxu23387Xv/51Jk2aREpKCgDl5eUMGzYMgN/97netristLY0tW7YA8Oqrr1Jb67uI7c0338yvf/1rKisrAThy5AjFxcVd/dYkBHXu8sgiIt2str6B5zcfJmNCKsP6JXTLNuNjo/nUrOE8vaGAksrzDEzq1S3ble519913M3v2bH772982PffYY49x1113MWzYMObNm8fBgwcvet2XvvQlli5dypw5c7jhhhtITEwE4KabbmLPnj3Mnz8f8DXE//73vyc1tWsPcUvoUZElImFh1Z4TnKw4z33dtBer0b1zRvKrDw7ywuYi/j5jbLduu0cZPDi4ZxgOHtzukMY9Se1JT0+/qO9u6dKlLF269KKxy5YtY9myZf4Ig9m4cWPTsh/84AdN3z/88MM8/PDDAW1fei4dLhSRsPDMpkMM6RtPxoRB3brdcalJzB09gGezDtHQoAb4Tjt+HJwL3tfx416/I5F2qcgSkZB3qLSKdftLuHt29zS8t3Tf3JEcKqviw/ySbt+2iIQvFVkiEvKezT5ElMHds0d4sv3FU69gQGIcf9ikBngRCZyKLBEJaTV1Dbyw+TDXTxzMkL7d0/DeUq+YaO68ejjv7j5BccXFp/JL63R9MelpOvqZVpElIiHt3d0nKKms4dNdfIX39tw7ZyR1DY4XNhd5miNcxMfHU1paqkJLegznHKWlpcTHxwf8Gp1dKCIh7Q9ZhQzrl8C1V3Zvw3tLowcmsmBsCn/YdIgHrhtLdJR5mifUDR8+nKKiIk6ePOl1lLBSXV3doX/E5dKCPZ/x8fEMHz484PEqskQkZBWUnOXDvFL+z99cGRJFzX1zR/LQH7aydv9JFk3QNY8uJTY2ltGjR3sdI+xkZmYyc+ZMr2P0GF7Ppw4XikjIejbrENFRxt961PDe0k2Tr2BgkhrgRSQwKrJEJCSdr6vnhS1F3DgplcF9QuPwSVxMFHdePYLVucUcL1cDvIhcmoosEQlJb+ecoOxsDffNHeV1lAvcO2cE9Q2O57MPex1FREKciiwRCUl/2FTIiAEJfGzcQK+jXGBUSiIfGz+Q57MPUa8rwIvIJajIEpGQc+BkJRsPlHHP7JFEhUDDe0v3zRnJ0fJq3t9X7HUUEQlhKrJEJOS8uct3X7pPzQr8VOnudOPkwfRNiOWNnbp/noi0LaAiy8wWm9leM8szsxWtLDcze8K/fIeZzfI/P8HMtjX7OmNm/xjsNyEiPcua3GKmDuvDFX1Do+G9pdjoKK69chCZe4t102gRaVO7RZaZRQNPAkuAycC9Zja5xbAlwHj/13Lg5wDOub3OuRnOuRnA1UAV8HLw4otIT3PqbA0fHTrF9SF+HarrJw6ipLKGnUfKvY4iIiEqkD1Zc4A859wB51wN8BywtMWYpcDTzmcj0M/MhrQYcwOQ75wrvOzUItJjvb/vJA0OFk0M7SLruitTMYNVuerLEpHWBVJkDQOan6tc5H+uo2PuAZ7taEARiSyrc4tJSYxj+vB+Xke5pAGJccwc0Y81KrJEpA2B3FantVN7WjYhXHKMmcUBnwS+2eZGzJbjO9TI4MGDyczMDCBa51VWVnb5NiKJ5jP4InFO6xscq3ZXMWNQDGvXvh/09Qd7Tkf3quFP+2t55a3V9IuPvPOIIvEz2tU0p8Hl9XwGUmQVAc3vaTEcONrBMUuAj5xzJ9raiHNuJbASID093WVkZAQQrfMyMzPp6m1EEs1n8EXinGYXlHH2nQ3cl3EVGVe17Di4fMGe09Qrz/Cn/euoSRlPRojc+qc7ReJntKtpToPL6/kM5L9e2cB4Mxvt3yN1D/BaizGvAff7zzKcB5Q75441W34vOlQoIu1YnVtMTJTxsStD6wKkbZk0JJkhfeNZrUOGItKKdvdkOefqzOwh4G0gGvi1cy7HzB7wL38KeAP4OJCH7wzCzze+3sx6A38DfDn48UWkJ1m9p5j0tP70iY/1OkpAzIyMCam8vv0oNXUNxMVE3iFDEWlbIIcLcc69ga+Qav7cU82+d8CDbby2Cki5jIwiEgGOnD7H3hMV/MvHJ3odpUOun5jKs1mHyC4oY2GI3QJIRLyl/3aJSEhoPOR2/cTBHifpmIXjUoiLiWLVHh0yFJELqcgSkZCwJreYkQN6M3ZQotdROqR3XAzzxqSwZq+KLBG5kIosEfFcdW096/NLuH5iKmahd0Po9lw/YRAHS85ysOSs11FEJISoyBIRz23IL6W6tiHkr/LelsZDnDrLUESaU5ElIp5blXuChNho5o4e4HWUThmZ0ptxqUm6+ruIXEBFloh4yjnHmtyTLBw3kPjYaK/jdNr1E1PZdLCUyvN1XkcRkRChIktEPLXvRCVHTp/jhknheaiw0aIJqdTWOz7YX+J1FBEJESqyRMRTjX1MiyaEd5GVntaf5PgYVue2efcwEYkwKrJExFNrcouZPKQPV/SN9zrKZYmNjuLa8YNYs/ckDQ3O6zgiEgJUZImIZ05X1bC5sIzrw/SswpYWTUzlZMV5co6e8TqKiIQAFVki4pn3952kwRG2l25oKWPCIMx0KQcR8VGRJSKeWZNbzIDEOGaM6Od1lKAYmNSL6cP7sVpXfxcRVGSJiEfqGxzv7ztJxpWDiI4Kv6u8t+X6iansKDpNSeV5r6OIiMdUZImIJ7YdPsWpqtoec6iw0QzG80QAACAASURBVPUTU3EOMvee9DqKiHhMRZaIeGJ1bjHRUca1Vw7yOkpQTRnah9TkXrqUg4ioyBIRb6zaU8zVo/rTNyHW6yhBZWYsmpDKun0l1NY3eB1HRDykIktEut3R0+fIPV7RYy7d0NKiialUnK8ju6DM6ygi4iEVWSLS7dbu8/UrhftV3ttyzfiBxEYZ7/95vddRRMRDKrJEpNttPFDKwKReXDk4yesoXSKpVwwz+kaxcUcBOF39XSRSqcgSkW7lnGPjgTLmjRmAWc+5dENL81Ki2ZWSRsWuXK+jiIhHVGSJSLcqLK3i+Jlq5o1J8TpKl5o3IJr6qGg2v7vJ6ygi4hEVWSLSrTYeKAXo8UXWrP7RxNbXsnF3kddRRMQjKrJEpFs19mONHZTodZQulRBtzDi+n431yerLEolQKrJEpNtESj9Wo3nH97Jr0Gj1ZYlEKBVZItJtIqUfq9G8o3vUlyUSwVRkiUi3iZR+rEazTuxXX5ZIBFORJSLdJlL6sRol1NUw4+g+9WWJRCgVWSLSLSKtH6vRvMM71ZclEqFUZIlIt4i0fqxG8w7tVF+WSIQKqMgys8VmttfM8sxsRSvLzcye8C/fYWazmi3rZ2Yvmlmume0xs/nBfAMiEh4irR+r0awjuerLEolQ7RZZZhYNPAksASYD95rZ5BbDlgDj/V/LgZ83W/ZT4C3n3ERgOrAnCLlFJMxEWj9Wo4S68+rLEolQgezJmgPkOecOOOdqgOeApS3GLAWedj4bgX5mNsTM+gDXAr8CcM7VOOdOBzG/iIQB5xybDkZeP1Yj9WWJRKaYAMYMAw43e1wEzA1gzDCgDjgJ/MbMpgNbgIedc2dbbsTMluPbC8bgwYPJzMwM8C10TmVlZZdvI5JoPoOvJ81pcVUDx8qrGVBX6ul76tY5rayExx6D+noSanpTXx7Nbzbt5qrSE92z/W7Qkz6joUJzGlxez2cgRVZr/+1suc+7rTExwCzgq865TWb2U2AF8O2LBju3ElgJkJ6e7jIyMgKI1nmZmZl09TYiieYz+HrSnD6ffQjYyf2L5zEuNdmzHN06p2vX+oqs8nLmxvTix//4HGc37iHji490z/a7QU/6jIYKzWlweT2fgRwuLAJGNHs8HDga4JgioMg513hazYv4ii4RiSAbD5QxMCmOsYOSvI7iiYS680w/pr4skUgTSJGVDYw3s9FmFgfcA7zWYsxrwP3+swznAeXOuWPOuePAYTOb4B93A7A7WOFFJPT5ro9VytwxKRHZj9Vo3iH1ZYlEmnaLLOdcHfAQ8Da+MwP/6JzLMbMHzOwB/7A3gANAHvBL4CvNVvFV4Bkz2wHMAL4fxPwiEuIOlVVxrDzyro/Vkq6XJRJ5AunJwjn3Br5CqvlzTzX73gEPtvHabUD6ZWQUkTDWeH2s+WMGeJzEW7OO/vV6WYu8DiMi3UJXfBeRLhXp/ViNeteqL0sk0qjIEpEuo36sCzX1ZeWoL0skEqjIEpEuo36sC/21LyvL6ygi0g1UZIlIl1E/1oWa+rJydB9DkUigIktEuoz6sS70176sJPVliUQAFVki0iXUj9U69WWJRA4VWSLSJdSP1Tr1ZYlEDhVZItIl1I/VOvVliUQOFVki0iXUj9U69WWJRA4VWSISdOrHujT1ZYlEBhVZIhJ06se6NPVliUQGFVkiEnSN/VjzRqsfqzVXH1FflkgkUJElIkGXXXCKAYlxjEtVP1ZrEurOM+14HpvreqsvS6QHU5ElIkG3uaCM9FH91Y91CbOLdrNj0Biq8w54HUVEuoiKLBEJquKKagpKq5idpkOFl5JelENtdAzbV6kvS6SnUpElIkG1ueAUAOlp/T1OEtrSi/YAsHlHgbdBRKTLxHgdQER6luyCMuJjo5g6rK/XUUJa/+oKxpcUkn1e/9cV6an0t1tEgmpzwSlmjuhPbLR+vbQnvWg3W1JGU3+yxOsoItIF9FtQRIKm8nwdOUfLma1DhQGZXbSbivhE9r233usoItIFVGSJSNBsPXSKBgfpanoPyOzDOQBkb97ncRIR6QoqskQkaLILThFlMGuU9mQFYviZYq6oKCG7+LzXUUSkC6jIEpGg2VxQxuShfUjqpXNqAmH4+rKy+wzHVVV5HUdEgkxFlogERW19A1sPnSZ9lA4VdsTsot0cTx7IkbWbvI4iIkGmIktEgiLn6BnO1dbrIqQd1NiXtXl9jsdJRCTYVGSJSFBsLigD0JmFHTShpJDk82fJ8l/EVUR6DhVZIhIU2QVljErpTWqfeK+jhJVo18CsI3vYHDcQ6uu9jiMiQaQiS0Qum3OOzQWn1I/VSbOLdrNv4EhOb97udRQRCSIVWSJy2Q6UnKX0bI0OFXZSY1/Wlve3epxERIIpoCLLzBab2V4zyzOzFa0sNzN7wr98h5nNaraswMx2mtk2M9sczPAiEhoa+7F0EdLOmX58P7H1tWTtPep1FBEJonYvZmNm0cCTwN8ARUC2mb3mnNvdbNgSYLz/ay7wc/+fjRY553RzLpEeKrvgFAMS4xg7KNHrKGEpvq6Gacfz2BydBM6BmdeRRCQIAtmTNQfIc84dcM7VAM8BS1uMWQo87Xw2Av3MbEiQs4pIiMouKCN9VH9MxUGnzS7azY5BY6jOO+B1FBEJkkAuyzwMONzscREX7qVqa8ww4BjggHfMzAG/cM6tbG0jZrYcWA4wePBgMjMzA8nfaZWVlV2+jUii+Qy+cJnT09UNFJaeY97AupDP261zWlkJjz0W8BmD8eeTqD0Tw+/W72DCkcPtvyAEhMtnNJxoToPL6/kMpMhq7b+mrgNjFjrnjppZKvCumeU659ZeNNhXfK0ESE9PdxkZGQFE67zMzEy6ehuRRPMZfOEyp2/sPAZ8xD03pDNzZGg3vnfrnK5d6yuyyssDGj49PpmfPvwsdZtyyfjcN7s2W5CEy2c0nGhOg8vr+QzkcGERMKLZ4+FAy+7MNsc45xr/LAZexnf4UUR6iOyCMuJjo5gytK/XUcJa/+oKxpcUkl2hk75FeopA/jZnA+PNbLSZxQH3AK+1GPMacL//LMN5QLlz7piZJZpZMoCZJQI3AbuCmF9EPJZdUMbMEf2Ji1FxcLnSi3azJWU09cUnvY4iIkHQ7m9F51wd8BDwNrAH+KNzLsfMHjCzB/zD3gAOAHnAL4Gv+J8fDHxgZtuBLOAvzrm3gvweRMQjlefr2H30jK6PFSRzDudQEZ/IvlUbvI4iIkEQSE8Wzrk38BVSzZ97qtn3DniwldcdAKZfZkYRCVFbD52iwen6WMGSXuS7Ms7mzfuYdK/HYUTksmn/voh0WnbBKaIMZo7s53WUHmH4mWKuqCgh60S111FEJAhUZIlIp2UfLGPy0D4kx8d6HaVHMHx7s7L7jsCdPet1HBG5TCqyRKRTausb2HpYN4UOtjmHcziePJAj67K8jiIil0lFloh0Ss7RM1TXNjBb/VhB1dSXtSHH4yQicrlUZIlIp/z1ptA6szCYJpQUknz+LNkHy7yOIiKXSUWWiHRK1sEyRg7ozeA+8V5H6VGiXQOzjuwhO25QwLfkEZHQpCJLRDrMOcfmwlM6VNhF5hzOYd/AkZzevN3rKCJyGVRkiUiHHSg5S9nZGl2EtIs09mVteX+rx0lE5HKoyBKRDvtrP5b2ZHWF6cf3E1tfS1Zuy9vEikg4UZElIh226WAZAxLjGDso0esoPVJ8XQ3TjueR3ZAEznkdR0Q6SUWWiHRYdkEZs9P6Y2ZeR+mx5hzexc5BY6jel+d1FBHpJBVZItIhx8urOVx2Tk3vXWzO4Rxqo2PY+p4uSioSrlRkiUiHZPn7seaMVpHVla4+sgdzDWTvLPQ6ioh0UozXAUQkvGQdLCUxLprJQ/p4HaVH63v+LBNOFpJVp/tCioQr7ckSkQ7JPniKWaP6ExOtXx9dbe7hXXw0cDR1R495HUVEOkG/JUUkYKerath7ooI56sfqFrOLcqiKSyDnnfVeRxGRTlCRJSIB21xwCoDZ6sfqFnMO+24Snf2RzjAUCUcqskQkYNkFZcRGGzNG9PM6SkRIPXuKUaeOknVK9zAUCUcqskQkYJsOljF9eD/iY6O9jhIxZhflkD1gNA2nTnkdRUQ6SEWWiASkqqaOXUfKdaiwm805nMOphD7kv6e+LJFwoyJLRAKy7dBp6hqcmt67WWNfVlb2Xo+TiEhHqcgSkYBkFZRhBrNG9fc6SkQZdfoYgyrLyD5W5XUUEekgFVkiEpCsg2VMuqIPfRN0cczuZMCcohyykofDuXNexxGRDlCRJSLtqq1vYOuh07qVjkdmH87haJ9BFL2/0esoItIBKrJEpF27jpRzrrZeN4X2SNP1sj7c5XESEekIFVki0q5s/02hZ49WP5YXJpQUklxdSVahLuMgEk5UZIlIu7IOljF6YCKpyfFeR4lI0a6B9CN7yIofDHV1XscRkQCpyBKRS2pocGQXnGJ2mvZieWn24RzyBwyndNNHXkcRkQAFVGSZ2WIz22tmeWa2opXlZmZP+JfvMLNZLZZHm9lWM/tzsIKLSPfYX1xJ+bla9WN5bG6Rrx8re+02j5OISKDaLbLMLBp4ElgCTAbuNbPJLYYtAcb7v5YDP2+x/GFgz2WnFZFul+Xvx9KZhd6adiyPXrXnyd573OsoIhKgQPZkzQHynHMHnHM1wHPA0hZjlgJPO5+NQD8zGwJgZsOBW4D/CWJuEekmWQfLGNynFyMH9PY6SkSLa6hjxrF9ZEX1A+e8jiMiAQikyBoGHG72uMj/XKBjfgL8E9DQyYwi4hHnHNkHy5idNgAz8zpOxJtTlEPOwNFU7tztdRQRCUBMAGNa+83a8r9RrY4xs1uBYufcFjPLuORGzJbjO9TI4MGDyczMDCBa51VWVnb5NiKJ5jP4QmFOT1Y1cPxMNf1qSz3PEgzdOqeVlfDYY1BfH7RV9qrpTUN5FL/N3svUspNBW29nhcJntKfRnAaX1/MZSJFVBIxo9ng4cDTAMXcCnzSzjwPxQB8z+71z7jMtN+KcWwmsBEhPT3cZGRmBvodOyczMpKu3EUk0n8EXCnP6py1FwHY+fdNcJg3p42mWYOjWOV271ldklZcHbZXpcQn8+OHnOL9hNxlfeCRo6+2sUPiM9jSa0+Dyej4DOVyYDYw3s9FmFgfcA7zWYsxrwP3+swznAeXOuWPOuW8654Y759L8r1vdWoElIqEp62AZfeJjmDA42esoAiTVnGPKiXw21SZ4HUVEAtBukeWcqwMeAt7Gd4bgH51zOWb2gJk94B/2BnAAyAN+CXyli/KKSDfKLvD1Y0VFqR8rVMwu2s22QWM4f7DA6ygi0o6ArpPlnHvDOXelc26sc+5x/3NPOeee8n/vnHMP+pdPc85tbmUdmc65W4MbX0S6ysmK8xwoOctsXbohpMw+nENNTBw739XNokVCna74LiKt2tx4v0JdhDSkzC7y3Sw6a9tBj5OISHtUZIlIqzYdLCM+Noppw/p6HUWaSTl3hnElh8iq0K9vkVCnv6Ui0qrsgjJmjuhPXIx+TYSa2UU5bBk4mvpi7y/jICJt029PEblIRXUte46dUT9WiJpzOIeKXonkvvuh11FE5BJUZInIRbYUnqLBwRz1Y4Wkxr6s7M37PU4iIpeiIktELpJ1sIzoKGPmyH5eR5FWDD9zkmHlxWwqqfU6iohcgoosEbnI+vxSZozoR2KvQG4KIV6Yd2gHGweMpuH0aa+jiEgbVGSJyAXOVNeyo+g0C8ameB1FLmFB4Q5OJfRhz5vrvI4iIm1QkSUiF8g6UEaDgwVjB3odRS5hwaHtAGzYlOtxEhFpi4osEbnAh/kl9IqJUj9WiBtSUcqY0iI+LGvwOoqItEFFlohcYEN+KbPTBhAfG+11FGnHgsLtZA0cS+2x415HEZFWqMgSkSYllefJPV7BfPVjhYUFh3ZwNi6BHW9+4HUUEWmFiiwRabIhvxRATe9hYt6hnQCs35zncRIRaY2KLBFpsj6/lOReMbpfYZgYcO4Mk0/ks/6sLrUhEopUZIlIk/X5JcwdM4CYaP1qCBcLCnewZdBYqg8UeB1FRFrQb1IRAaDoVBWFpVW6dEOYWVi4nZqYOLa8qfsYioQaFVkiAjTrxxqnfqxwMrsoh+iGetZvL/Q6ioi0oCJLRABfP1ZKYhxXpiZ7HUU6IKnmHNOP7ePD2kRwzus4ItKMiiwRwTnH+vwS5o9NISrKvI4jHbSwcDs7Bo3hzK49XkcRkWZUZIkI+SfPcuLMefVjhan5hdtpiIoi6+1NXkcRkWZUZIkIG/JLAFiofqywNOtILr1qz7N+9xGvo4hIMyqyRIQP80oZ1i+BkQN6ex1FOiG+vpb0I3tYHzUAGnQvQ5FQoSJLJMI1NDg2HChlwdgUzNSPFa4WFG4nd+AoSrK2eh1FRPxUZIlEuN3HzlB+rlaXbghzCwq3A7Bx1WaPk4hIIxVZIhFuvb8fS03v4W3a8TySz5/lw/0nvY4iIn4qskQi3Pr8UsYOSmRwn3ivo8hliHENzD20iw3xV0BtrddxRAQVWSIRraaugayDZdqL1UPMP7Sdgn5DOPK+LuUgEgpUZIlEsB1Fp6mqqdelG3qIhQW+vqz172/zOImIgIoskYi2Pr8UM5g3RkVWT3BlySFSzp5mQ2G511FEhACLLDNbbGZ7zSzPzFa0stzM7An/8h1mNsv/fLyZZZnZdjPLMbPvBvsNiEjnfZhXwpShfejXO87rKBIEUTjmH9rBh31G4KqqvI4jEvHaLbLMLBp4ElgCTAbuNbPJLYYtAcb7v5YDP/c/fx643jk3HZgBLDazeUHKLiKX4VxNPVsPnVY/Vg+zoHA7J5JSOPDeh15HEYl4gezJmgPkOecOOOdqgOeApS3GLAWedj4bgX5mNsT/uNI/Jtb/pdvEi4SAzYVl1NQ3MH+sDhX2JAsKdwCw/oNdHicRkZgAxgwDDjd7XATMDWDMMOCYf0/YFmAc8KRzrtXTXsxsOb69YAwePJjMzMxA8ndaZWVll28jkmg+g6+r5/SFvTVEG5w/nEPmsd1dtp1Q0q2f08pKeOwxqK/vnu35OQcpZbW8mjCMEfo9GnY0p8Hl9XwGUmS1dp+Nlnuj2hzjnKsHZphZP+BlM5vqnLvov1jOuZXASoD09HSXkZERQLTOy8zMpKu3EUk0n8HX1XP645wPmTnSWHzjgi7bRqjp1s/p2rW+Iqu8+5vQF338Yd4bN49rZ8wgql+/LtuO/t4Hn+Y0uLyez0AOFxYBI5o9Hg4c7egY59xpIBNY3OGUIhJU5edq2Vl0mgXj1I/VEy0s2M7phGR2v7nO6ygiES2QIisbGG9mo80sDrgHeK3FmNeA+/1nGc4Dyp1zx8xskH8PFmaWANwI5AYxv4h0QtbBMhocLFA/Vo+04JCvL2vDJv26FfFSu4cLnXN1ZvYQ8DYQDfzaOZdjZg/4lz8FvAF8HMgDqoDP+18+BPidvy8rCvijc+7PwX8bItIRH+aVEB8bxcyRXXcoSbwzuLKMsaWHWVcDX/I6jEgEC6QnC+fcG/gKqebPPdXsewc82MrrdgAzLzOjiASRc45VuSdYOHYgvWKivY4jXSTjwBb+d+bHqSo4RO+0kV7HEYlIuuK7SITJK67kcNk5rp+U6nUU6UI35GVRExPHBy9neh1FJGKpyBKJMKtyiwG4fqKKrJ4svWg3yefPsnrnEa+jiEQsFVkiEWb1nmImD+nDkL4JXkeRLhTXUMe1Bz9ide/hNOgWOyKeUJElEkFOV9WwubCMG3SoMCLckJdFcdIAdr2+xusoIhFJRZZIBMnce5IGBzdMGux1FOkGGQe2YK6BVR9ExhX9RUKNiiyRCLIqt5iBSXFcNayv11GkGww4d4ZZR3JZXZXgu9+OiHQrFVkiEaK2voH39xazaEIqUVGt3QlLeqLr87PZOWg0JzZu9TqKSMRRkSUSIbYUnuJMdZ36sSLMDXlZAKx5c4PHSUQij4oskQixOreYuOgorhk/yOso0o0mlBQyrPwE7xWe9TqKSMRRkSUSIVbtOcHcMQNI6hXQjR6khzDghrxsPhw4luojx72OIxJRVGSJRICCkrPknzzLDboAaUS6Pj+Lc7HxbHh5tddRRCKKiiyRCPDXq7zr0g2RaN6hnfSuOcfqrYVeRxGJKCqyRCLA6twTjE9NYmRKb6+jiAfi62u5pmAbq3sNwZ0/73UckYihIkukh6uormXTgTLdEDrC3ZCfxZE+g8h9432vo4hEDBVZIj3cuv0l1DU4btChwoi2KH8zAKszd3icRCRyqMgS6eFW7SmmX+9YZo3s53UU8VDq2VNcdWwfqypidfV3kW6iIkukB6tvcKzZW0zGlYOIidZf90h3fX42WweNoXRbjtdRRCKCfuuK9GDbDp+m7GwN1+uG0ALcmJeFsyjWvP6h11FEIoKKLJEebHXuCaKjjOt0lXcBppzIZ3BFKasPnPI6ikhEUJEl0oOt2lNM+qj+9O0d63UUCQGG75Dh2v5jqTlZ4nUckR5PRZZID3Xk9Dlyj1fohtBygevzs6js1Zvsl1Z5HUWkx1ORJdJDrd5zAoAb1I8lzSws3E5cXQ2rNh/wOopIj6ciS6SHWpVbTFpKb8YMTPQ6ioSQ3rXnWVi4nVXRqbiaGq/jiPRoKrJEeqCqmjrW55dy/cTBmJnXcSTEXJ+XRWG/K8h/V2cZinQlFVkiPdCHeaXU1DWoH0tadX3j1d9Xb/U4iUjPpiJLpAd6d/dxknrFMDttgNdRJAQNqzjJxOKDvHtK/wSIdCX9DRPpYapr63lz13FumjKYuBj9FZfW3Zq7juzU8RRt3OZ1FJEeS7+BRXqY1bnFVFTXcfvMYV5HkRC2NCcTgFf/tNbbICI9WEBFlpktNrO9ZpZnZitaWW5m9oR/+Q4zm+V/foSZrTGzPWaWY2YPB/sNiMiFXt56hNTkXiwYO9DrKBLCRpwpZs7hXbxUkYCrr/c6jkiP1G6RZWbRwJPAEmAycK+ZTW4xbAkw3v+1HPi5//k64P845yYB84AHW3mtiATJqbM1ZO4tZumMoURH6axCubTbctaQ328oOX/O9DqKSI8UyJ6sOUCec+6Ac64GeA5Y2mLMUuBp57MR6GdmQ5xzx5xzHwE45yqAPYCOYYh0kb/sPEZtveM2HSqUANyS+wFxdbW8/I76skS6gjnnLj3A7E5gsXPui/7HnwXmOuceajbmz8APnXMf+B+vAv7ZObe52Zg0YC0w1Tl3ppXtLMe3F4zBgwdf/dxzz13eO2tHZWUlSUlJXbqNSKL5DL7OzOnjG89RVef414UJuj5WK7r1c1pZCXl5EOKH4v6rfBh5dQn85w3JREd3rE1Xf++DT3MaXN01n4sWLdrinEtv+XxMAK9t7Td1y8rskmPMLAn4E/CPrRVYAM65lcBKgPT0dJeRkRFAtM7LzMykq7cRSTSfwdfROT1UWsX+t9bwT4snsChjXNcFC2Pd+jlduxYeewzKy7tne51UfeV8Hrj9W8QWnObaL9zRodfq733waU6Dy+v5DOS/LUXAiGaPhwNHAx1jZrH4CqxnnHMvdT6qiFzKK9uOALB0hg4VSuAy8jfTp7qSVzbkex1FpMcJpMjKBsab2WgziwPuAV5rMeY14H7/WYbzgHLn3DHzHa/4FbDHOfefQU0uIk2cc7yy9QhzRw9gWL8Er+NIGImvr+WW3A94q89oqo4Vex1HpEdpt8hyztUBDwFv42tc/6NzLsfMHjCzB/zD3gAOAHnAL4Gv+J9fCHwWuN7Mtvm/Ph7sNyES6XYUlXOg5KyujSWdclvOGqriEnj39294HUWkRwmkJwvn3Bv4Cqnmzz3V7HsHPNjK6z6g9X4tEQmil7ceIS4miiXThngdRcLQ7KLdDCs/wcvFpy86dVxEOk9XfBcJc7X1Dby+/Sg3Tkqlb0Ks13EkDEXhuC0nk3UDx3Fy1z6v44j0GCqyRMLcB3kllJ6t4TY1vMtluG13JvVR0bz+7HteRxHpMVRkiYS5V7YeoV/vWDImpHodRcLY+NLDTD2exysnHLRz/UQRCYyKLJEwVnm+jrdzjnPLtCHExeivs1ye23LWsGNgGnmrN3odRaRH0G9lkTD29q7jVNc26KxCCYpP7llLVEM9r762wesoIj2CiiyRMPbKtiOMGJDA1aP6ex1FeoDUs6dYWLidl2v742pqvI4jEvZUZImEqRNnqvkwr4TbZwzTfQolaO7YtZqiPqlseUkN8CKXS0WWSJh6fftRGhws1aFCCaKb9m8koaaal9fkeB1FJOypyBIJUy9vPcL04X0ZO6jr7zAvkSOxtpqb92/gzwkjqDkd2je3Fgl1KrJEwtC+ExXkHD3DbdqLJV3gtpw1lMcnk6nb7IhcFhVZImHo5a1HiI4ybr1qqNdRpAe6pmAbA8+e4qWtR7yOIhLWVGSJhJlzNfU8l3WIRRNSGZTcy+s40gPFuAY+tWs176ZcSdGu/V7HEQlbKrJEwswLWw5zqqqW5deO8TqK9GCf2/I65hy/+tXbXkcRCVsqskTCSH2D43/WHWTGiH7MTtO1saTrDK0o4ZN71vK8DaH82Emv44iEJRVZImHk7ZzjHCqr4svXjtG1saTLfSnrJari4vn9z172OopIWFKRJRImnHP8Yu0B0lJ6c9OUK7yOIxFg0skCrj2whd+UJ1JdWeV1HJGwoyJLJExkHSxj++HTfOFjY4iO0l4s6R5fznqJkt79ePUX2psl0lEqskTCxMq1BxiQGMddVw/3OopEkAWF25lyPI+VeedoqKv3Oo5IWFGRJRIG9p+oYFVuMffPH0V8bLTXcSSCGLA86yXy+w5h9TO6OKlIR6jIEgkDv1x3gPjYKO6fn+Z1FIlAt+R+wLDyYlZ+eMjrKCJhRUWWSIgrPlPNK1uPctfVIxiQGOd1HIlAMa6Bv9v8ClkD0tj65gdexxEJGyqy3wefxAAAFORJREFURELcb9YXUNfQwBc/NtrrKBLB7tn+Dn2qK1n56havo4iEDRVZIiGs8nwdz2wsZPHUKxiVkuh1HIlgibXVfGbrG7zVdwwFH+3xOo5IWFCRJRLCns8+zJnqOpZfO9brKCIs2/I6sfX1/M9v3/U6ikhYUJElEqLqGhy//uAgc0YPYMaIfl7HESH17Cluz1nNC3EjKD10zOs4IiFPRZZIiMo+Xs+R0+f4sm4ELSHkS9kvcz4mjqd//qrXUURCnooskRDknOPNg7WMHZTIogmpXscRaTKutIgb92/i6ap+nK/VxUlFLkVFlkgI+jCvlEMVDSy/dgxRuoWOhJgvZb3EqYRkPth3xusoIiEtoCLLzBab2V4zyzOzFa0sNzN7wr98h5nNarbs12ZWbP9/e/ceH0V97nH88+xu7oSEa4SAQi2gAoKK1mJVihbFU8Hqq5Wec5RTrfSiPVrraWsvXipt0VJ7tFgV6/0GtWIBb4DVqBVRFAgXAxhBIARBLgmEJCTZ/fWPHWSbJmRJdjO5fN+v1zgzv/mN+/Awuzz7m9kZs9WJDFyko3LOcferH5KTZlx0Ur7f4Yj8m9NK1jCydC0vlhrV1TV+hyPSZjVZZJlZELgHGA+cAHzTzE6o1208MMibpgD3xmx7BDg/EcGKdAZzV5Ty7sbdXHRsCmkhPUJH2h4Dfvz6o+wKp/CnO5/xOxyRNiuekazTgGLn3AbnXA0wC5hYr89E4DEXtQTINbM+AM65N4DdiQxapKMqr6pl6gtFjOify9n9Q36HI9Ko0ZtXcXpaOfftyWLjOj1uR6Qh8RRZ+cCWmPUSr+1I+4hIE/6waD279h9g6sRhBEzXYknbNilrB2l1Ndx8z8s45/wOR6TNieerckOf9PXfTfH0OfyLmE0heqqRvLw8CgoKjmT3I1ZRUZH01+hMlM+W+7g8zKNvVzP26BC7ipcrp0nQqjmtqIBbboFwx/0FXqhfHyYcqODJ/flMf/QFTh3Qxe+Q2j297xPL73zGU2SVAP1j1vsBpc3oc1jOuZnATIBRo0a5MWPGHMnuR6ygoIBkv0Znony2TCTi+P97F9OjS4Q7vzWGnIwU5TQJWjWnb7wRLbLKy1vn9XxQMH06t978Y5ZP/gPPVvbge18/ly5Z6X6H1a7pfZ9YfuczntOFS4FBZjbQzFKBScC8en3mAZd7vzI8HSh3zul2wCJxmv3eFlZsKeNnFxxPTkaK3+GIxC3kIty28F4+yczl7t/N9jsckTalySLLOVcHXAMsAIqAvzjn1pjZd83su163F4ENQDHwAPD9g/ub2dPA28AQMysxsysT/GcQadd276/h9pfXctrA7nxNt2yQduiU0rVcWriAh/bnsn5lsd/hiLQZcf18yTn3ItFCKrbtvphlB1zdyL7fbEmAIh3dHS+vpaK6jtsmDsN0sbu0Uz95/VEWDP4iv3jgNWbffayOZRF0x3cRX72/aQ+zlm7hii8NZMhR2X6HI9Js3av28pPXH+XdrL4898RCv8MRaRNUZIn4pC4c4Zd/W81RXdO59pxBfocj0mKXFi5kZOlafvPebsrLKvwOR8R3KrJEfPL4kk18sG0vN114AllpuvGotH8BHFMX3svutC78/ne6E7yIiiwRH+zYW82dC9dz1uBejB92lN/hiCTMsO0fcdnyF3iitierlhb5HY6Ir1RkibQy5xw3z1vDgboIt04YqguEpcO5/s0n6F5Vzs8ee5vq2o57M1aRpqjIEmllf3y1mJdWf8L14wYzsGeW3+GIJFzOgf1MXfgnVmXlcePtc/TIHem0VGSJtKL5haXcuWg9F5+Uz3fO+pzf4Ygkzfnr3+b6N5/guYpM7nlokd/hiPhCRZZIK1m+eQ83PFPIqGO68dtLhus0oXR4P1g8i4lrCpj+YS0vLFzmdzgirU5Flkgr2FpWxVWPvU/vrmncf9kppIWCfockknQG3P7SXZy8tYjrF26icM1mv0MSaVUqskSSrOJAHVc+spQDtWEemnwqPbqk+R2SSKtJD9cyc85Ueu3fw7f/vJjSnfv8Dkmk1ajIEkmicMRx3azlrN++jxn/dTKD8nRXd+l8elaW8+Bff0WVBfn2b+ayv7rW75BEWoWKLJEkmvZSEa8U7eCWCUM5e3Avv8MR8c2QnZv449zbWRvK4dpf/5VwRL84lI5PRZZIkjz97mYeeHMjk794DJd/cYDf4Yj47ssb3uemv8/kldqu3DHjeb/DEUk6FVkiSbC4eCe//Ntqzhrci19+9QS/wxFpMyYve57Llj3P/aUBZs9Z7Hc4IkmlB6aJJNjzK0u54ZlCBvbMYsZ/nkQoqO8yIgcZcPMrM/m4W19ujIykPPAeV008Rbc0kQ5Jn/4iCRKJOKYvWMc1Ty1nWN8cnrrqdLqmp/gdlkibE3IR7n/u14xfv5jfLNnO9X9cqMfvSIekIkskAfZV1zLl8feZ8Voxk07tz1NXnU6vbN2qQaQxmbUHmDF3Gj9643GeK63j0tvm8kl5td9hiSSUiiyRFtq0az8X/2kxr63bwa0ThvLbi4eTGtJbS6QpBvzg7dncP2cqxRURJkx9nuUbPvU7LJGE0b8EIi3wjw93MmHGW3xacYDHrziNyaMH6NoSkSN03odLmPP4DaRV7OXSexfzbEGR3yGJJISKLJFmcM7x8Fsbmfzwu+R1TWPe1V9i9Od7+h2WSLs1ZOcm5j1yHadsLeJHL29g6iOvUxeO+B2WSIuoyBI5Qtv3VvPD2Su4df4HjD2uN3O+fwZH98j0OyyRdq9b9T4em/0LJr8/nz+vrWDytOf56NMKv8MSaTbdwkEkTnura7n/9Y948B8bCUcc154ziGvPGUQgoNODIomSEglz6yv3c/yOjdx2zlWMm/4a3xjak+suOom8rul+hydyRFRkiTThQF2YJ5ZsZsarH7KnspaJI/vyo68M0eiVSBJNWrmQc4vfYcYZk3gyMp7n1izk26OPZsp5Q3VrFGk3VGSJNCISccwt3Mr0BevZWlbFmYN68pPzj2NYfo7foYl0Cj0ry7ll0f18a+lcfn/2ZGZYkCeXbOKaccfz32ceS1oo6HeIIoelIkukntpwhIJ1n3LnovUUbdvL0L5dmXbJcM4cpAc8i/jhmLJPuHvu7UxZ8lduH3sFty0I8fDfi7hu4kn8x4h8MlJVbEnbpCJLBAhHHO9s3MX8wm28vHobeypr6d89g7smjeTCE/vquiuRNmDY9o94/Omf8+aAkUz78pXc8GyIm+YU8pXjevHVUwdy1uCeGt2SNkVFlnRakYhj+ZY9zC/cxgurtvHpvgNkpgY59/g8LhzRl7MH99JNRUXaoDM/XsEZD/8vS44ezvzhY3mp8gvMLdpFdiDC+cOO4qujBjD62B6k6Lmh4jMVWdJpOOfYsruKFSVlLN+8h4VrtrO1rIrUUICxQ3pz4Yi+jD2ut049iLQDARyjN69k9OaV/CoQ5K1jRjB/+Dm8XDWKZ1buoHvIMW54X045thcj++fyuV5dCGpEWlqZiizpsHbsq2bllnJWlpSxoqScVSVl7KmsBSAtFGD0sT244bzBnHt8Htn6tZJIu5USCTNm4zLGbFxGdTCFNwaezPwTz+H5/SOZtXwbAFkBx/C+2Yz4XG9G9M/lxH455Odm6AkNklRxFVlmdj5wFxAE/uycm1Zvu3nbLwAqgf9xzi2LZ1+R5qqqCbO1rJKSPVWU7Klia1kVW735lt2V7Nh3AICAweC8bMadcNRnH65DjsrWqQSRDig9XMu44ncYV/wOEYwNPfIp7DuEwv5DKSwdwMObBlITjH6p6pYCR3fPpF+vruR3zyQ/NyM6dYtOulWEtFSTRZaZBYF7gK8AJcBSM5vnnPsgptt4YJA3fQG4F/hCnPtKJ+Ccoy7iqAs7aiMR6sKO6tow1bVhqrx5dW2Eqpow1XVhKmvC7K2qpdybyiq9eVUte6tq2VNZQ5k3KnVQKGD0yU0nPzeDMwf14vg+2Yzon8vQvl3JTNWgrUhnE8Dx+V0lfH5XCZes+jsAB4Ih1vUaQGH+8XzQdxAlWT0oyu7JopzenxVfB2WHoFtGCjmZqeR2SScnK5WcjBRyM1PIyYhOGakhMlKCpKcEvHl0ykgNkh4KEAoGSAkaoUB0rpGzziWef3lOA4qdcxsAzGwWMBGILZQmAo855xywxMxyzawPMCCOfVvd9bNXULihirs+eMvPMFrEuSa2N7GDq9fsvJbP1l20j3POW3aftUW8hYhzhJ0jEoGq6mpSFr9COBLdJ+y8giocIRyJFljNYQZd0//1Q+3o7pnkZITok3PoW2e/bhn0zk7XNRciclhp4TpO/KSYEz8phvcPtUcwdmblsLVrb7bm5rG1Rz6l2b0oS8ukPDWTsrQsStKzKc/Ipjwti3CgedduBg1CB6eAETQIBIyAQdCMmnCYzH8sIBAwgmZgYES3mxle06Flr2gz7z+Gt+59FFrsSgPMPuvV8PbDb27zKiurGDPGv9ePp8jKB7bErJcQHa1qqk9+nPsCYGZTgCkAeXl5FBQUxBFa8+zZeYAUC1O7f1/SXqNVNHHg19/8b90Pvlkb6XfozRyzboe2B7z1AEY4NUJqaphATHvQIBgIRucGwUD0gyXofbCkBCE1YKQGITUAqcGDy9F5ZoqREYLAZx8Qdd5UdSjYcthfDus+hnWHT0e7U1FRkdT3QWfUqjmtqoKbbmr6G1E7VtGvHwXTp/sdRsJ1AYZ407/aB+zDOah2Afa7ADUuQI0zavDmLmaOEXZGGKPOGREgXK/NAREg4gznoCYri+D+8mib92kb/cILLnadmC/L9T7FXWPL9Q7FjntkHpKSkU7BokWQ4s+p33iKrIb+Ka//d9NYn3j2jTY6NxOYCTBq1Cg3Joml55gxUFBQQDJfo7NRPhNPOU28Vs/p+PGt91o+KCgoYMyll/odRocSPUYv8juMDsPvz9F4iqwSoH/Mej+gNM4+qXHsKyIiItLhxPPzqqXAIDMbaGapwCRgXr0+84DLLep0oNw5ty3OfUVEREQ6nCZHspxzdWZ2DbCA6G0YHnLOrTGz73rb7wNeJHr7hmKit3D41uH2TcqfRERERKQNiet37c65F4kWUrFt98UsO+DqePcVERER6eh0N0YRERGRJFCRJSIiIpIEKrJEREREkkBFloiIiEgSqMgSERERSQIVWSIiIiJJoCJLREREJAlUZImIiIgkgYosERERkSSw6M3a2xYz+xTYlOSX6QnsTPJrdCbKZ+Ipp4mnnCaW8pl4ymlitVY+j3HO9arf2CaLrNZgZu8550b5HUdHoXwmnnKaeMppYimfiaecJpbf+dTpQhEREZEkUJElIiIikgSducia6XcAHYzymXjKaeIpp4mlfCaecppYvuaz016TJSIiIpJMnXkkS0RERCRpOmyRZWYPmdkOM1vdwLYbzMyZWc+YthvNrNjM1pnZea0bbdvXUD7N7BYz22pmK7zpgphtymcTGjtGzewHXt7WmNkdMe3K6WE0cozOjjk+PzazFTHblM8mNJLTkWa2xMvpe2Z2Wsw25fQwGsnnCDN728xWmdl8M+sas035PAwz629mr5lZkfd5ea3X3t3MFpnZh968W8w+rZtT51yHnICzgJOB1fXa+wMLiN6Hq6fXdgJQCKQBA4GPgKDff4a2NDWUT+AW4IYG+iqfzc/pl4FXgDRvvbdy2vx81tv+e+Am5bNlOQUWAuO95QuAAuW0RflcCpztLV8B3KZ8xp3PPsDJ3nI2sN7L2x3AT732nwK3+5XTDjuS5Zx7A9jdwKY/AD8GYi9GmwjMcs4dcM5tBIqB0xrYt9M6TD4bonzGoZGcfg+Y5pw74PXZ4bUrp0043DFqZgZ8A3jaa1I+49BITh1wcLQlByj1lpXTJjSSzyHAG97yIuASb1n5bIJzbptzbpm3vA8oAvKJ5u5Rr9ujwEXecqvntMMWWQ0xswnAVudcYb1N+cCWmPUSr02ado2ZrfSGwQ8OySqfzTcYONPM3jGz183sVK9dOW2ZM4HtzrkPvXXls/muA35nZluA6cCNXrty2jyrgQne8teJnm0B5fOImNkA4CTgHSDPObcNooUY0Nvr1uo57TRFlpllAj8HbmpocwNt+tll0+4FjgVGAtuIno4B5bMlQkA34HTg/4C/eKMwymnLfJNDo1igfLbE94AfOuf6Az8EHvTaldPmuQK42szeJ3rKq8ZrVz7jZGZdgGeB65xzew/XtYG2pOa00xRZRIuBgUChmX0M9AOWmdlRRKvZ/jF9+3FoCFwa4Zzb7pwLO+ciwAMcGnZVPpuvBJjjot4FIkSfvaWcNpOZhYCLgdkxzcpn800G5njLz6D3fYs459Y658Y5504h+kXgI2+T8hkHM0shWmA96Zw7eFxuN7M+3vY+wMHLLlo9p52myHLOrXLO9XbODXDODSCa7JOdc58A84BJZpZmZgOBQcC7PobbLhw8iD1fIzrsDcpnS/wNGAtgZoOBVKIPN1VOm+9cYK1zriSmTflsvlLgbG95LHDwFKxy2gxm1tubB4BfAPd5m5TPJnij/A8CRc65O2M2zSP6ZQBvPjemvVVzGkrm/9xPZvY0MAboaWYlwM3OuQcb6uucW2NmfwE+AOqAq51z4VYLth1oKJ/AGDMbSXS49WPgO6B8xquRnD4EPOT9xLsGmOyiP4tRTptwmPf8JP71VKGO0Tg1coxeBdzljRBWA1NAOY1HI/nsYmZXe13mAA+D8hmnM4DLgFUxt2f5GTCN6KUWVwKbiV7r5ktOdcd3ERERkSToNKcLRURERFqTiiwRERGRJFCRJSIiIpIEKrJEREREkkBFloiIiEgSqMgSERERSQIVWSIiIiJJoCJLRDo0MzvVe4h5upllmdkaMxvmd1wi0vHpZqQi0uGZ2VQgHcgASpxzv/U5JBHpBFRkiUiHZ2apwFKij4EZrceTiEhr0OlCEekMugNdgGyiI1oiIkmnkSwR6fDMbB4wCxgI9HHOXeNzSCLSCYT8DkBEJJnM7HKgzjn3lJkFgcVmNtY596rfsYlIx6aRLBEREZEk0DVZIiIiIkmgIktEREQkCVRkiYiIiCSBiiwRERGRJFCRJSIiIpIEKrJEREREkkBFloiIiEgSqMgSERERSYJ/ApIfTzg/gLtIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(10,6))\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "ax.plot(x, f_x, label='Standard Normal Distribution')\n",
    "ax.fill_between(x[np.where(x.astype(int)==173)[0][0]:], y1=f_x[np.where(x.astype(int)==173)[0][0]:],\n",
    "                color='red', label=\"P-value\")\n",
    "ax.grid()\n",
    "ax.legend()\n",
    "ax.set_xlabel('x')\n",
    "ax.set_title('P-Value and AUC for z-stat = 0.6')\n",
    "p_val = np.sum(f_x[np.where(x.astype(int)==173)[0][0]:])\n",
    "print('p-val (AUC) = {:.3f}'.format(p_val))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# let's see for the body weight problem\n",
    "x = np.linspace(140, 200, 64)\n",
    "mu = 170  # H_0 is true!\n",
    "sigma = 5  # calcualted for 64 samples\n",
    "f_x = (1 / np.sqrt(2 * np.pi * sigma ** 2)) * np.exp(- (x - mu) ** 2 / (2 * sigma ** 2))\n",
    "x_normed = (x - mu) / sigma\n",
    "\n",
    "def plot_p_auc():\n",
    "    fig = plt.figure(figsize=(10,6))\n",
    "    ax = fig.add_subplot(1,1,1)\n",
    "    ax.plot(x, f_x, label='Standard Normal Distribution')\n",
    "    ax.fill_between(x[np.where(x.astype(int)==185)[0][0]:], y1=f_x[np.where(x.astype(int)==185)[0][0]:],\n",
    "                    color='red', label=\"P-value\")\n",
    "    ax.grid()\n",
    "    ax.legend()\n",
    "    ax.set_xlim([170, 195])\n",
    "    ax.set_ylim([0, 0.05])\n",
    "    ax.set_xlabel('x')\n",
    "    ax.set_title('P-Value and AUC for z-stat = 3.0')\n",
    "    p_val = np.sum(f_x[np.where(x.astype(int)==185)[0][0]:])\n",
    "    print('p-val (AUC) = {:.3f}'.format(p_val))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p-val (AUC) = 0.001\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAGDCAYAAABnZBdiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXxU9b3/8ddnJvtCIGFPAgmyrwmyCigoKlQUcANti1TrUrW12muvvVq1/mqvra239dZK1brWq6hVXKm4EHFBQSAgqyyyhJ2wEyDb9/fHDBhiQgZIcmaS9/PxmIcz53zPmc/kC+bN95z5fs05h4iIiIh4w+d1ASIiIiKNmcKYiIiIiIcUxkREREQ8pDAmIiIi4iGFMREREREPKYyJiIiIeEhhTKSBMbOnzey3XtdREzNzZtbR6zqOMLN4M3vTzPaY2cte1yMijYfCmEg9MLO1ZnbQzPab2VYze8rMkqpoN9jMDphZchX7FpjZzfVTcXgIBstSM2tbxfbfVtqWFQx4URW2XWlmXwZ/7pvNbLqZDa3m7S4FWgFpzrnLav3DnAQzG25mBSd4TJ2EXDObaGYrgmF1m5k9Y2ZNjtM+x8zmmVlR8L85tV2TSEOhMCZSfy50ziUBfYH+wF2VGzjnZgMFwCUVt5tZT6A78EI91BkWzCyRwM9hD/D9kzj+NuDPwO8IhKx2wN+AsdUc0h742jlXehLvFVVzq4j3KTDEOZcCdACigCpHYM0sBngd+CfQDHgGeD24XUQqURgTqWfOuY3AdKBnNU2eASZV2jYJeNs5VwhgZi+b2ZbgKMUsM+tR1YnMbLKZfVJp29GREzOLNbM/mtn64IjdFDOLr+Zcp5nZh2ZWaGY7zOx5M2taYf9aM/sPM1sUrGuqmcVV2H97cHRqk5ldffyfEhAIYruB+4CrQmhfsdaU4HE3Oededc4dcM6VOOfedM7dXkX73wB3AxOCo2jXmJnPzO4ys3XBkaBng+etOAp3jZmtBz6s4pxvBs915FFuZpOrqfd7ZrbUzPaZ2cbgzzGRwJ+TthXO0dbMBpjZbDPbHfx5/vVIyDGzWcFTLgy2n3AiP7fjcc5tcM7tqLCpDKhuBG44gbD2Z+fcYefcw4ABZ9dWPSINicKYSD0zs0zge8CCapo8Bwwzs3bB9j7gSuDZCm2mA52AlsB84PmTLOf3QGcgh8Av1nQCoaTK0oH/BtoC3YBM4N5KbS4HRgHZQG9gcvAzjAL+Azg3WPfIEGq7isBI4ItAVzPrG9pHAmAwEAe8Fkpj59w9BEbQpjrnkpxz/wjWPhkYQWAkKAn4a6VDzyLwszi/inNeGDxXEoFLoFuAD6op4R/A9c65ZAIh/UPn3AFgNLDpyHmcc5sIhKBbgebBz3kOcGPwPc8Mnq9PsP3Uym9kZkODQa66R3WXcY8cuwfYRyAs/7mapj2ARe7Y9fYWBbeLSCUKYyL1Z5qZ7QY+AT4i8Mv/O5xzG4L7fxDcdA6BYPF2hTZPOuf2OecOEwhEfY6M2oTKzAy4FrjVObfTObcvWNPEaupa5Zx7LzjSsR14iEAYqehh59wm59xO4E0CIQ8CIe0p59ziYMi4t4ba2hEIQf/nnNtKIMScyOhYGrDjZC45VvB94CHn3Brn3H7gV8DESpck7w2Ouh2s7iRm1plAkJ4Q7NuqlADdzayJc26Xc25+dedzzs1zzn3unCt1zq0F/s53+6FazrlPnHNNj/P4pIZjU4AM4EFgbTVNkwhcXq5oD/CdeyFFRGFMpD6NC/6ya++cu9E5d9DMhlW4BLWkQtuKlyp/SCCUlACYmd/MHjCz1Wa2l29/ITY/wXpaAAnAvCOjIsC/g9u/w8xamtmLwctoewncD1T5PbdUeF5E4JcyBEbTKgaRdTXU9kNgmXMuP/j6eeBKM4sOvi4FoisdEw2UBx+FQPNTvJerbaU61xG49NaqwrbqwhVw9HLp68CvnXMfB7f9V4U+nxJsegmB0dJ1ZvaRmQ0+zjk7m9lbwcvUewkE6BPt+1MSvNT+bwKjllXZD1S+ub8JgRE1EalEYUzEQ865jytcgqp4CedVIN3MRgAXc+wlyisJ3IQ+EkgBsoLbrYq3OEAgcAUamLWusG8HcBDoUWFUJCV4Wa0q/w04oLdzrgmBkbuq3rMqmwlc1jyiXQ3tJwEdgoFjC4FRuOYELtsBrOfbz31ENrDBOVcOzAYOAeNCrK8qmwjc1F+x5lJga4VtjmoELy//HzDTOff3owc497sKfX5DcNtc59xYApedpwEvHef8jwLLgU7BfvgvQu8HKv0DoKrHsBBPFQWcVs2+JUDv4OjrEb2D20WkEoUxkTAUvJT3CvAUsM4592WF3cnAYQKjPwlUc7kzaCHQwwLTDMRR4fJgMLQ8DvyPmbUEMLN0M/vO/U8V3nc/sNvM0oHv3Ah/HC8Bk82su5klAPdU1zA4KnQaMIDAZc4cAvdR/R/fXqr8F3CBmZ0XHClsS+DbqS8GP9seAve+PWJm48wswcyizWy0mf0hxJpfAG41s2wLTENy5J6yUC993g8kArccr5GZxZjZ980sJTj6uZfAfWEQCH5plS5BJwfb7DezrsBPKp1yK4F73KpU6R8AVT0+rqbO75tZOwtoH/x81d0Dlxf8DD+zwJdEjkzJ8p0vOoiIwphIOHuGwMjMs5W2P0vgktlGYCnweXUncM59TeBbhe8DKwncr1bRfwKrgM+Dl7zeB7pUc7rfEJiWYw+B+9deDfWDOOemE7jZ+8Pg+x3vl/JVwOvOua+cc1uOPIC/AGPMLNU5twS4gsBo3U4CI2FfBGs88p4PAbcRCGnbCVxSvJnAyFMoniTwZYpZwDcERtp+GuKxBOsbBOyqMOpU3RQdPwTWBvvgBoL3CzrnlhMIhWuCl5LbEvgixJUELvk9DlS+Sf9e4Jlg+8tPoN6adAc+IxDIPwVWELjnEAALzOH2X8G6iwmMSk4i8I3Yqwlcpi+uxXpEGgw79ssuIiIiIlKfNDImIiIi4qGQwpiZjbLAMhirzOyOKvabmT0c3L+o4nxAFpgI8iszyzezLysfKyIiItKY1fi1bzPzA48QmKyxAJhrZm8455ZWaDaawESOnYCBBL7tM7DC/hGVZm4WEREREUIbGRsArApOfFhM4NtKldd2Gws86wI+B5qaWZtarlVERESkwQkljKVz7MSGBcFtobZxwAwzm2dm151soSIiIiINUSizU1c1mWDlr2Aer80Q59ym4DxG75nZcufcrMqNg0HtOoC4uLjT27WraU5I8VqZgw37ykmJMZrFBf4IlJeX4/PpeyGRSv0XudR3kU39F9m+/vrrHc65KlcvCUUoYayAY2fOziAwM3VIbYIL2+Kc22ZmrxG47PmdMOacewx4DKBLly5uxYoVIX4E8dLkp+awcut+Pv7lCHw+Iy8vj+HDh3tdlpwk9V/kUt9FNvVfZDOzmpZ4O65QYvhcoFNwFuoYAosIv1GpzRvApOC3KgcBe5xzm80s0cySg4UmAucBi0+lYAkv43LS2bj7IF+u2+V1KSIiIhGpxpEx51xpcCmLdwE/8KRzbomZHVlTbQrwDoFFblcRWBz4R8HDWwGvBZcniyKw2PG/a/1TiGfO7d6K+Gg/0/I3MiA71etyREREIk4olylxzr1DIHBV3DalwnMH3FTFcWuAPqdYo4SxxNgozuvRine+2sy9F/ao+QARERE5RkhhTOR4xuWk83r+Jj76ejvRXhcjImGhpKSEgoICDh065HUpESElJYVly5Z5XYbUIC4ujoyMDKKja/e3ncKYnLKhnZqTmhjDtPyNXNbW62pEJBwUFBSQnJxMVlYWwVtV5Dj27dtHcnKy12XIcTjnKCwspKCggOzs7Fo9t75HK6cs2u9jTO82vL90KwdLtfC8iMChQ4dIS0tTEJMGw8xIS0urk9FehTGpFWNz0jlcWs68raVelyIiYUJBTBqauvozrTAmtaJvu6a0S01g9iaFMREJD/fffz89evSgd+/e5OTk8MUXXwDw5z//maKiolp7n6ysLHbsOPnll/Py8rjsssuq3G5mvPnmm0e3jRkzhry8vJN+r5NR3efLysqiV69e9OrVi+7du3PXXXdx+PBhADZt2sSll15a7Tl3797N3/72t+O+7xlnnAEEfg5jxow5oZqnTZvG0qXfLqF999138/7775/QOeqTwpjUCjNjbE5blhaWs22fbtgVEW/Nnj2bt956i/nz57No0SLef/99MjMDc5PXdhg7UWVlZSG3zcjI4P7776+X9zoZM2fO5KuvvmLOnDmsWbOG664LrHrYtm1bXnnllWqPO14YO1LzZ599dtJ1VQ5j9913HyNHjjzp89U1hTGpNWNz0nHAmws3e12KiDRymzdvpnnz5sTGxgLQvHlz2rZty8MPP8ymTZsYMWIEI0aMAOAnP/kJ/fr1o0ePHtxzzz1Hz5GVlcU999xD37596dWrF8uXLwegsLCQ8847j9zcXK6//noCszsFjBs3jtNPP50ePXrw2GOPHd2elJTE3XffzcCBA5k9ezb//ve/6dq1K0OHDuXVV1+t9nP06dOHlJQU3nvvve/s++CDD8jNzaVXr15cffXVR0elsrKyuO+++xg6dCgvv/wyw4cP59Zbb+XMM8+kW7duzJ07l4svvphOnTpx11131Vh7KJKSkpgyZQrTpk1j586drF27lp49ewKwZMkSBgwYQE5ODr1792blypXccccdrF69mpycHG6//Xby8vIYMWIEV155Jb169Tp6ziP27t3L+PHj6d69OzfccAPl5eXfafPKK68wefJkPvvsM9544w1uv/12cnJyWL16NZMnTz4aDo/3c6uqv+uDvk0ptaZjyyTaN/Hxev5Grhlau980EZHI9Zs3l7B0095aPWf3tk245zhzG5533nncd999dO7cmZEjRzJhwgTOOussfvazn/HQQw8xc+ZMmjdvDgQuZ6amplJWVsY555zDokWL6N27NxAIcfPnz+dvf/sbf/zjH3niiSf4zW9+w9ChQ7n77rt5++23jwkuTz75JKmpqRw8eJD+/ftzySWXkJaWxoEDB+jZsyf33Xcfhw4dolOnTnz44Yd07NiRCRMmHPez3nXXXdx1112ce+65R7cdOnSIyZMn88EHH9C5c2cmTZrEo48+ys9//nMgMAXDJ598AsCUKVOIiYlh1qxZ/OUvf2Hs2LHMmzeP1NRUTjvtNG699VbS0tKqrT1UTZo0ITs7m5UrV9KqVauj26dMmcItt9zC97//fYqLiykrK+OBBx5g8eLF5OfnA4FLkXPmzGHx4sVVflNxzpw5LF26lPbt2zNq1CheffXVai+DnnHGGVx00UWMGTPmO21q+rlV1d/1QSNjUqsGt4liUcEe1mzf73UpItKIJSUlMW/ePB577DFatGjBhAkTePrpp6ts+9JLL9G3b19yc3NZsmTJMZe3Lr74YgBOP/101q5dC8CsWbP4wQ9+AMAFF1xAs2bNjrZ/+OGH6dOnD4MGDWLDhg2sXLkSAL/fzyWXXALA8uXLyc7OplOnTpjZ0XNVZ9iwYQB8/PHHR7etWLGC7OxsOnfuDMBVV13FrFnfLvtcOeBddNFFAPTq1YsePXrQpk0bYmNj6dChAxs2bDhu7Sei4ijhEYMHD+Z3v/sdv//971m3bh3x8fFVHjtgwIBqp4wYMGAAHTp0wO/3c8UVVxwNmieqpp9bVf1dHzQyJrVqYBs/U7+GafmbuO3czl6XIyJh4HgjWHXJ7/czfPhwhg8fTq9evXjmmWeYPHnyMW2++eYb/vjHPzJ37lyaNWvG5MmTj5m64MhlTr/fT2npt19QqupbdXl5ebz//vvMnj2bhIQEhg8ffvRccXFx+P3+4x5/PHfeeSf3338/UVGBX9tVhZ6KEhMTj3l95HP4fL6jz4+8Li0tPW7todq3bx9r166lc+fO7Nmz5+j2K6+8koEDB/L2229z/vnn88QTT9ChQ4caa66o8s/ryOuK20Opt6afW3X9Xdc0Mia1qlmcjzNOS+P1/I01/qEXEakrK1asOGZkJz8/n/bt2wOQnJzMvn37gMC9SImJiaSkpLB161amT59e47nPPPNMnn/+eQCmT5/Orl27ANizZw/NmjUjISGB5cuX8/nnn1d5fNeuXfnmm29YvXo1AC+88EKN73neeeexa9cuFi5cePQca9euZdWqVQA899xznHXWWTWepzqh1l6d/fv3c+ONNzJu3LhjRgoB1qxZQ4cOHfjZz37GRRddxKJFi47pg1DMmTOHb775hvLycqZOncrQoUMBaNWqFcuWLaO8vJzXXnvtaPvqzl/bP7faojAmtW5sTjrrCovI37Db61JEpJHav38/V111Fd27d6d3794sXbqUe++9F4DrrruO0aNHM2LECPr06UNubi49evTg6quvZsiQITWe+5577mHWrFn07duXGTNm0K5dOwBGjRpFaWkpvXv35te//jWDBg2q8vi4uDgee+wxLrjgAoYOHXo0JNbkzjvvpKCg4Og5nnrqKS677DJ69eqFz+fjhhtuCOk8VQm19spGjBhBz549GTBgAO3atePvf//7d9pMnTqVnj17kpOTw/Lly5k0aRJpaWkMGTKEnj17cvvtt9f4PoMHD+aOO+6gZ8+eZGdnM378eAAeeOABxowZw9lnn02bNm2Otp84cSIPPvggubm5R0Mv1P7PrbZYOI5edOnSxa1YscLrMuQk5OXl0XfQEPr99n2uHNCOey/S4uGRJC8vj+HDh3tdhpyEcOu7ZcuW0a1bN6/LiBhaDilyVPVn28zmOef6new5NTImta5JXDQju7XkrUWbKC0r97ocERGRsKYwJnVibE46O/YX88mqk5+VWkREpDFQGJM6MbxLC5rERfF6/iavSxEREQlrCmNSJ2Kj/FzQuw3vLtlCUbHWqxQREamOwpjUmbE56RQVl/He0q1elyIiIhK2FMakzgzISqVtSpwuVYqIiByHwpjUGZ/PuDCnLbO+3k7h/sNelyMijYzf7ycnJ4eePXty2WWXUVRUdMrnfPrpp7n55ptroTqRbymMSZ0an5tOabnj7a82e12KiHipdWswq71H69Y1vmV8fDz5+fksXryYmJgYpkyZUg8fVOTEKYxJneraugldWyfz2oKNXpciIl7aWsv3jp7g+YYNG3Z0CZwjysvLycrKYvfub1cL6dixI1u3buXNN99k4MCB5ObmMnLkSLZW8X6TJ0/mlVdeOfo6KSnp6PMHH3yQ/v3707t3b+65554TqlUaH4UxqXPjc9NZsH433+w44HUpItIIlZaWMn36dHr16nXMdp/Px9ixY4+uafjFF1+QlZVFq1atGDp0KJ9//jkLFixg4sSJ/OEPfwj5/WbMmMHKlSuZM2cO+fn5zJs3j1mzZtXqZ5KGRWFM6tzYnHTMYJpGx0SkHh08eJCcnBz69etHu3btuOaaa77TZsKECUydOhWAF198kQkTJgBQUFDA+eefT69evXjwwQdZsmRJyO87Y8YMZsyYQW5uLn379mX58uXHLFouUlmU1wVIw9c6JY4zTktjWv5Gfj6yE2bmdUki0ggcuWesokceeYTHH38cgHfeeYfBgwezatUqtm/fzrRp07jrrrsA+OlPf8ptt93GRRddRF5e3tFFxiuKioqivDyw5JtzjuLi4qPPf/WrX3H99dfX4aeThkQjY1IvxudmsK6wiPnrd3ldiog0YjfddBP5+fnk5+fTtm1bzIzx48dz22230a1bN9LS0gDYs2cP6enpADzzzDNVnisrK4t58+YB8Prrr1NSUgLA+eefz5NPPsn+/fsB2LhxI9u2bavrjyYRTGFM6sWonq2Ji/bpRn4RCTsTJkzgn//859FLlAD33nsvl112GcOGDaN58+ZVHnfttdfy0UcfMWDAAL744gsSExMBOO+887jyyisZPHgwvXr14tJLL2Xfvn318lkkMplzzusavqNLly5uxYoVXpchJyEvL4/hw4dXue9nLyxg1srtzPmvkcRE6d8B4eh4/SfhLdz6btmyZXTr1u3bDa1b1+43Klu1gi1bau98Htu3bx/JyclelyEh+M6fbcDM5jnn+p3sOfUbUerN+Nx0dheVkLdCw/Uijc6WLeBc7T0aUBATURiTejOsU3PSEmN0qVJERKQChTGpN1F+Hxf2acsHy7ax52CJ1+WIiIiEBYUxqVcX902nuKycd7Q8kkiDF473JIucirr6M60wJvWqV3oKHVok6lKlSAMXFxdHYWGhApk0GM45CgsLiYuLq/Vza9JXqVdmxsW56fxxxtds2FlEZmqC1yWJSB3IyMigoKCA7du3e11KRDh06FCd/JKX2hUXF0dGRkatn1dhTOrd2JxAGHtj4SZuGtHR63JEpA5ER0eTnZ3tdRkRIy8vj9zcXK/LEI/oMqXUu8zUBAZkpfLq/AJdwhARkUZPYUw8MS43ndXbD7B4416vSxEREfGUwph44oJebYjx+3h1QYHXpYiIiHhKYUw8kZIQzdldW/Lmwk2UlpV7XY6IiIhnFMbEM+P7prNjfzEfr9rhdSkiIiKeURgTz4zo0pKmCdFM05xjIiLSiCmMiWdionxc0KsN7y7Zwv7DpV6XIyIi4gmFMfHUxX3TOVRSzruLt3hdioiIiCcUxsRTfds1o11qgpZHEhGRRkthTDxlZozLTefT1TvYsueQ1+WIiIjUO4Ux8dz43HScgzcWanRMREQaH4Ux8Vx280RyMpvy6nyFMRERaXwUxiQsjM9NZ/mWfSzbrOWRRESkcVEYk7AwpncbonymOcdERKTRURiTsJCWFMtZnVvwev4mysqd1+WIiIjUG4UxCRvj+6azZe8hPl9T6HUpIiIi9UZhTMLGyG6tSI6N0pxjIiLSqCiMSdiIi/Yzuldrpn+1mYPFZV6XIyIiUi8UxiSsjMtN50BxGe8t2+p1KSIiIvVCYUzCyqDsNNqkxPHa/AKvSxEREakXCmMSVnw+Y2xOOrNW7mDH/sNelyMiIlLnQgpjZjbKzFaY2Sozu6OK/WZmDwf3LzKzvpX2+81sgZm9VVuFS8N1cd90ysodby7c5HUpIiIida7GMGZmfuARYDTQHbjCzLpXajYa6BR8XAc8Wmn/LcCyU65WGoXOrZLp3qaJJoAVEZFGIZSRsQHAKufcGudcMfAiMLZSm7HAsy7gc6CpmbUBMLMM4ALgiVqsWxq4i/ums7BgD6u37/e6FBERkToVShhLBzZUeF0Q3BZqmz8DvwTKT7JGaYQu6tMWn6HRMRERafCiQmhjVWyrvF5NlW3MbAywzTk3z8yGH/dNzK4jcImTFi1akJeXF0JpEm72799fa33XPdXPC7NXkxu9CZ9V9UdMaltt9p/UL/VdZFP/NW6hhLECILPC6wyg8p3V1bW5FLjIzL4HxAFNzOyfzrkfVH4T59xjwGMAXbp0ccOHDw/1M0gYycvLo7b6bmeTAm57aSHJ2X3on5VaK+eU46vN/pP6pb6LbOq/xi2Uy5RzgU5mlm1mMcBE4I1Kbd4AJgW/VTkI2OOc2+yc+5VzLsM5lxU87sOqgphIVc7v0Zr4aD+vztelShERabhqDGPOuVLgZuBdAt+IfMk5t8TMbjCzG4LN3gHWAKuAx4Eb66heaUQSY6M4v0cr3l60iUMlWh5JREQaplAuU+Kce4dA4Kq4bUqF5w64qYZz5AF5J1yhNGrj+2YwLX8TeSu2MapnG6/LERERqXWagV/C2pDT0mieFKtLlSIi0mApjElYi/L7GJvTlpkrtrG7qNjrckRERGqdwpiEvfG56ZSUOd5atNnrUkRERGqdwpiEvR5tm9CpZZImgBURkQZJYUzCnpkxLjedL9ftYn1hkdfliIiI1CqFMYkI43IDq2tNy9fomIiINCwKYxIR0pvGMzA7ldcWbCQwk4qIiEjDoDAmEePivul8s+MACwv2eF2KiIhIrVEYk4gxqmcbYqJ8vDa/wOtSREREao3CmESMlPhozu3WijcXbaakrNzrckRERGqFwphElHG56ew8UMysr7d7XYqIiEitUBiTiHJW5xY0S4jmNc05JiIiDYTCmESUmCgfY3q35b2lW9l7qMTrckRERE6ZwphEnPF90zlcWs6/v9ridSkiIiKnTGFMIk5uZlM6NE/kpS83eF2KiIjIKVMYk4hjZkzon8mX63axats+r8sRERE5JQpjEpEu7ptBlM94cY5Gx0REJLIpjElEapEcy7ndW/Hqgo0cLi3zuhwREZGTpjAmEWtC/0x2Hijm/aXbvC5FRETkpCmMScQa1qkF6U3jeXHueq9LEREROWkKYxKx/D7jsn4ZfLxyBxt2FnldjoiIyElRGJOIdlm/TMzgZU1zISIiEUphTCJaetN4zuzUgpe+LKBUi4eLiEgEUhiTiHfFgEy27D3ErJVaPFxERCKPwphEvLO7tqJ5UozmHBMRkYikMCYRLybKxyV9M/hg+Ta27T3kdTkiIiInRGFMGoQJ/TMpK3e8Mr/A61JEREROiMKYNAgdWiQxIDuVqXM34JzzuhwREZGQKYxJgzGxfybrCouYvabQ61JERERCpjAmDcb3erUhOS6KqXN1I7+IiEQOhTFpMOKi/YzPTWf64i3sLir2uhwREZGQKIxJgzKhfybFpeW8tmCj16WIiIiERGFMGpQebVPonZHCi3N0I7+IiEQGhTFpcCb0z2TF1n0sLNjjdSkiIiI1UhiTBueiPm2Jj/Yzde56r0sRERGpkcKYNDjJcdGM6d2GN/I3ceBwqdfliIiIHJfCmDRIEwdkcqC4jLcWbfK6FBERkeNSGJMGqW+7ZnRsmcSLmnNMRETCnMKYNEhmxsT+mSxYv5sVW/Z5XY6IiEi1FMakwbq4bwbRfuNF3cgvIiJhTGFMGqzUxBjO69Ga1xZs5FBJmdfliIiIVElhTBq0K/q3Y3dRCe8u2eJ1KSIiIlVSGJMG7YzT0shoFq/Fw0VEJGwpjEmD5vMZE/pl8tnqQtYVHvC6HBERke9QGJMG77J+mfgMjY6JiEhYUhiTBq91ShwjurTk5XkFlJaVe12OiIjIMRTGpFGY0D+T7fsOM3PFdq9LEREROYbCmDQKZ3dtScvkWF6coznHREQkvEGpQWcAACAASURBVCiMSaMQ5fdx6ekZzFyxjS17DnldjoiIyFEKY9JoTOifSbmDV+bpRn4REQkfCmPSaLRPS+SM09KY+uUGysud1+WIiIgACmPSyEzon8mGnQf5bHWh16WIiIgACmPSyJzfozVNE6K1eLiIiIQNhTFpVOKi/YzPTWfGkq3sPFDsdTkiIiIKY9L4TOifSXFZOa/OL/C6FBEREYUxaXy6tm5CTmZTps7dgHO6kV9ERLylMCaN0hUDMlm5bT/z1+/2uhQREWnkQgpjZjbKzFaY2Sozu6OK/WZmDwf3LzKzvsHtcWY2x8wWmtkSM/tNbX8AkZMxpndbEmP8mpFfREQ8V2MYMzM/8AgwGugOXGFm3Ss1Gw10Cj6uAx4Nbj8MnO2c6wPkAKPMbFAt1S5y0hJjo7iwT1veWrSZfYdKvC5HREQasVBGxgYAq5xza5xzxcCLwNhKbcYCz7qAz4GmZtYm+Hp/sE108KGbdCQsTBzQjoMlZby5cLPXpYiISCMWFUKbdKDi+jEFwMAQ2qQDm4Mja/OAjsAjzrkvqnoTM7uOwKgaLVq0IC8vL5T6Jczs378/YvrOOUdGkvH4B0toe3CN1+WEhUjqPzmW+i6yqf8at1DCmFWxrfLoVrVtnHNlQI6ZNQVeM7OezrnF32ns3GPAYwBdunRxw4cPD6E0CTd5eXlEUt9dE/MNv3lzKS0659KjbYrX5Xgu0vpPvqW+i2zqv8YtlMuUBUBmhdcZwKYTbeOc2w3kAaNOuEqROjI+N52YKB8vzdXi4SIi4o1QwthcoJOZZZtZDDAReKNSmzeAScFvVQ4C9jjnNptZi+CIGGYWD4wEltdi/SKnpGlCDKN7tua1BRs5VFLmdTkiItII1RjGnHOlwM3Au8Ay4CXn3BIzu8HMbgg2ewdYA6wCHgduDG5vA8w0s0UEQt17zrm3avkziJySCf0z2XuolOmLdSO/iIjUv1DuGcM59w6BwFVx25QKzx1wUxXHLQJyT7FGkTo1uEMaWWkJvDhnA+NzM7wuR0REGhnNwC+Nnplxef9MvvhmJ2u276/5ABERkVqkMCYCXNo3A7/PmKob+UVEpJ4pjIkALZvEMbJbS176cgMHi3Ujv4iI1B+FMZGgHw/rwK6iEl6ZX+B1KSIi0ogojIkE9WvfjD6ZTfnHx2soK9eqXSIiUj8UxkSCzIxrh2WztrCI95dt9bocERFpJBTGRCoY1aM16U3jeeJjrVUpIiL1Q2FMpIIov4+rh2Yzd+0u8jfs9rocERFpBBTGRCqZ0D+T5LgoHtfomIiI1AOFMZFKkmKjuHJgO6Z/tZkNO4u8LkdERBo4hTGRKkw+IwufGU99utbrUkREpIFTGBOpQpuUeC7s05apc9ez52CJ1+WIiEgDpjAmUo0fD8vmQHEZL8xZ73UpIiLSgCmMiVSjR9sUzjgtjac/XUtxabnX5YiISAOlMCZyHNcO68CWvYd4+6tNXpciIiINlMKYyHGc1bkFHVsm8fisb3BOSySJiEjtUxgTOQ6fL7BE0tLNe5m9utDrckREpAFSGBOpwdicdJonxWgSWBERqRMKYyI1iIv2M2lwFjNXbGfl1n1elyMiIg2MwphICH4wqD1x0T6e+Pgbr0sREZEGRmFMJASpiTFc0jeD1xZsZPu+w16XIyIiDYjCmEiIrhmaTUl5Oc/NXut1KSIi0oAojImEqEOLJEZ2a8Vzn6/jYHGZ1+WIiEgDoTAmcgKuHdaBXUUl/Gt+gdeliIhIA6EwJnIC+mc1o09GCv/45BvKyzUJrIiInDqFMZETYGb8eFgHvtlxgPeXbfW6HBERaQAUxkRO0OierUlvGq9pLkREpFYojImcoCi/jx8NyWLO2p0s3LDb63JERCTCKYyJnISJA9qRHBelJZJEROSUKYyJnISk2CiuHNCO6Yu3ULCryOtyREQkgimMiZykyUOyMOCpT9d6XYqIiEQwhTGRk9QmJZ4xvdvw4pz17DlY4nU5IiISoRTGRE7Bj4d14EBxGS/OWe91KSIiEqEUxkROQc/0FAZ3SOPpz9ZSUlbudTkiIhKBFMZETtG1Z2azec8h3l602etSREQkAimMiZyi4Z1b0rFlEo9/vAbntESSiIicGIUxkVPk8xk/HprNkk17mb2m0OtyREQkwiiMidSCcbnpNE+K0RJJIiJywhTGRGpBXLSfHw7K4sPl21i1bZ/X5YiISARRGBOpJT8Y1I7YKJ9Gx0RE5IQojInUkrSkWC45PYNXF2xk+77DXpcjIiIRQmFMpBZdMzSbkrJynvt8ndeliIhIhFAYE6lFp7VI4pyurfjn5+s4VFLmdTkiIhIBFMZEatm1w7LZeaCYf80v8LoUERGJAApjIrVsQHYqvTNS+MfH31BerklgRUTk+BTGRGqZmfHjYR1Ys+MAHyzf5nU5IiIS5hTGROrA93q2Jr1pPI9/vMbrUkREJMwpjInUgSi/jx8NyWLONztZVLDb63JERCSMKYyJ1JEJ/TNJjo3icU0CKyIix6EwJlJHkuOiuXJgO975ajOrtu33uhwREQlTCmMidei6MzsQH+3nwXeXe12KiIiEKYUxkTqUlhTL9Wd24N0lW5m3bqfX5YiISBhSGBOpY9cMy6Zlciy/e2c5zmneMREROZbCmEgdS4iJ4tZzOzNv3S5mLN3qdTkiIhJmFMZE6sFlp2dwWotE/vDv5ZSWlXtdjoiIhBGFMZF6EOX38Z+jurJ6+wFe+lJrVoqIyLdCCmNmNsrMVpjZKjO7o4r9ZmYPB/cvMrO+we2ZZjbTzJaZ2RIzu6W2P4BIpDi3eyv6tW/G/7z/NUXFpV6XIyIiYaLGMGZmfuARYDTQHbjCzLpXajYa6BR8XAc8GtxeCvzCOdcNGATcVMWxIo2CmfGr73Vl+77DPKGJYEVEJCiUkbEBwCrn3BrnXDHwIjC2UpuxwLMu4HOgqZm1cc5tds7NB3DO7QOWAem1WL9IRDm9fSqjerTm7x+tZsf+w16XIyIiYSAqhDbpwIYKrwuAgSG0SQc2H9lgZllALvBFVW9iZtcRGFWjRYsW5OXlhVCahJv9+/er72pwVrNyZpSU8Z/P5vHD7rFel3MM9V/kUt9FNvVf4xZKGLMqtlWeLOm4bcwsCfgX8HPn3N6q3sQ59xjwGECXLl3c8OHDQyhNwk1eXh7qu5otLvmKqXM38OvL+5PVPNHrco5S/0Uu9V1kU/81bqFcpiwAMiu8zgA2hdrGzKIJBLHnnXOvnnypIg3HLSM7ERPl48EZK7wuRUREPBZKGJsLdDKzbDOLASYCb1Rq8wYwKfitykHAHufcZjMz4B/AMufcQ7VauUgEa5kcx4+HdeDtRZvJ37Db63JERMRDNYYx51wpcDPwLoEb8F9yzi0xsxvM7IZgs3eANcAq4HHgxuD2IcAPgbPNLD/4+F5tfwiRSHTdmR1onhTDf7+zTMskiYg0YqHcM4Zz7h0CgavitikVnjvgpiqO+4Sq7ycTafSSYqO45ZxO/Pr1JcxcsY2zu7byuiQREfGAZuAX8dDEAe3Ibp7I76evoKxco2MiIo2RwpiIh6L9Pm4/vwsrtu7jX/O1TJKISGOkMCbisdE9W5OT2ZSHZnzNoZIyr8sREZF6pjAm4jEz41eju7Jl7yGe/FTLJImINDYKYyJhYGCHNEZ2a8mjeavZdaDY63JERKQeKYyJhIn/HNWVA4dL+evMVV6XIiIi9UhhTCRMdGqVzGWnZ/Lc7HVs2FnkdTkiIlJPFMZEwsit53bG54M/aZkkEZFGQ2FMJIy0Tonj6iHZTMvfxOKNe7wuR0RE6oHCmEiYuWH4aTRLiOaB6cu9LkVEROqBwphImGkSF81Pz+7EJ6t2MOvr7V6XIyIidUxhTCQMfX9QOzJT43lg+nLKtUySiEiDpjAmEoZio/z8x3ldWLp5L68v3Oh1OSIiUocUxkTC1IW929IrPYU/vqtlkkREGjKFMZEw5fMZd4zuysbdB3lu9jqvyxERkTqiMCYSxoZ0bM5ZnVvw15mr2FNU4nU5IiJSBxTGRMLcHaO7svdQCX/7SMskiYg0RApjImGuW5smjM9N56lP17Jp90GvyxERkVqmMCYSAX5xXhcAHnrva48rERGR2qYwJhIB0pvGM/mMLP41v4DlW/Z6XY6IiNQihTGRCHHj8NNIjo3SMkkiIg2MwphIhGiaEMPNZ3ckb8V2Plu9w+tyRESkliiMiUSQSYOzSG+qZZJERBoShTGRCBIX7ee2czuzqGAPU7/c4HU5IiJSCxTGRCLM+Nx0hnRM47dvLWV9YZHX5YiIyClSGBOJMD6f8eClffD5jF+8nE+ZLleKiEQ0hTGRCNS2aTy/uagHc9fu4vGP13hdjoiInAKFMZEINT43ndE9W/PQjK9Ztllzj4mIRCqFMZEIZWbcP74XTeKjuXVqPodLy7wuSUREToLCmEgES02M4Q+X9mL5ln38z3srvS5HREROgsKYSIQ7u2srrhiQyd9nrWbu2p1elyMiIidIYUykAbjzgu5kNIvntpfy2X+41OtyRETkBCiMiTQASbFRPHR5DgW7DnL/20u9LkdERE6AwphIA9E/K5XrzuzAC3M28MGyrV6XIyIiIVIYE2lAbju3M11bJ/Of//qKnQeKvS5HRERCoDAm0oDERvn5nwk57DlYzJ2vfYVzmp1fRCTcKYyJNDDd2jThtnO7MH3xFl5bsNHrckREpAYKYyIN0HVndqBf+2bc8/oSNu0+6HU5IiJyHApjIg2Q32f86fI+lDnHf7y8kHItJi4iErYUxkQaqPZpifx6THc+W13I05+t9bocERGphsKYSAM2sX8mZ3dtye//vZxV2/Z5XY6IiFRBYUykATMzHrikFwkxfm6dupCSsnKvSxIRkUoUxkQauJbJcfxufC++2riH//1wldfliIhIJQpjIo3A6F5tGJ+bziMzV5G/YbfX5YiISAUKYyKNxL0X9aBlciy3Tc3nYHGZ1+WIiEiQwphII5ESH80fL+vDmh0HeGD6Mq/LERGRIIUxkUZkSMfmTD4ji2dmr+Pjldu9LkdERFAYE2l07hjdldNaJHL7y4vYU1TidTkiIo2ewphIIxMXHVhMfPv+w9zzxmKvyxERafQUxkQaod4ZTfnp2R2Zlr+Jtxdt9rocEZFGTWFMpJG6aURH+mSkcOe0r9i295DX5YiINFoKYyKNVLTfx58uz+FgcRm//NcinNNi4iIiXlAYE2nEOrZM4leju5K3Yjv/N2e91+WIiDRKCmMijdykwVkM7dic3761jC0HtHaliEh9UxgTaeR8PuMPl/Ym2m/8ef4hCvcf9rokEZFGRWFMRGjbNJ4nrupP4UHHj56ey/7DpV6XJCLSaIQUxsxslJmtMLNVZnZHFfvNzB4O7l9kZn0r7HvSzLaZmSY0EgljA7JTuTEnliWb9nL9c19yuFTrV4qI1Icaw5iZ+YFHgNFAd+AKM+teqdlooFPwcR3waIV9TwOjaqNYEalbuS2j+P0lvfl0VSG3TV1IWbm+YSkiUtdCGRkbAKxyzq1xzhUDLwJjK7UZCzzrAj4HmppZGwDn3CxgZ20WLSJ159LTM7jze914+6vN3P36Yk15ISJSx6JCaJMObKjwugAYGEKbdEBTe4tEoGvP7EDhgWKmfLSatMQYbjuvi9cliYg0WKGEMatiW+V/KofS5vhvYnYdgUuctGjRgry8vBM5XMLE/v371XcRrGL/DYxzLEmP4uEPV1G4eT3nto/2tjg5Lv3di2zqv8YtlDBWAGRWeJ0BbDqJNsflnHsMeAygS5cubvjw4SdyuISJvLw81HeRq3L/DTuznBufn8/zS7cyoE93xuake1ecHJf+7kU29V/jFso9Y3OBTmaWbWYxwETgjUpt3gAmBb9VOQjY45zTJUqRCBfl9/HwFbkMzE7lFy8tJG/FNq9LEhFpcGoMY865UuBm4F1gGfCSc26Jmd1gZjcEm70DrAFWAY8DNx453sxeAGYDXcyswMyuqeXPICJ1KC7az+NX9aNzq2R+8s/5zFu3y+uSREQalFAuU+Kce4dA4Kq4bUqF5w64qZpjrziVAkXEe03ionnm6gFcOuUzrn56Li/fMJjOrZK9LktEpEHQDPwiEpIWybE8d/VAYqJ8TPrHHAp2FXldkohIg6AwJiIha5eWwLNXD+BAcSmT/jFH61iKiNQChTEROSHd2jThycn92bj7oNaxFBGpBQpjInLC+mel8ugP+modSxGRWqAwJiIn5eyurfhDcB3LW6fmax1LEZGTFNK3KUVEqnLJ6RnsKirmt28vo2nCYu4f1xOzqhbkEBGR6iiMicgp+fGwwDqWj+YF1rH8hdaxFBE5IQpjInLKfnl+F3YdKOZ/P1xFamIMPxqS7XVJIiIRQ2FMRE6ZmfHbcT3ZVVTMb95cSmpijNaxFBEJkW7gF5FaEeX38ZeJuQzqEFjHcqbWsRQRCYnCmIjUmrhoP49P6keX1sn85J/ztI6liEgIFMZEpFYlx0Xz9I8G0LpJHFc/PZevt+7zuiQRkbCmMCYita5FcizPXTOQ2Cgfl/99NrO+3u51SSIiYUthTETqRGZqAi/fMJjWTeK46qk5/PXDlZRrYlgRke9QGBOROtM+LZFXbzyDi/q05Y8zvub6f85j76ESr8sSEQkrCmMiUqcSYqL484Qc7rmwOzOXb2PsXz9lxRbdRyYicoTCmIjUOTPjR0OyeeG6Qew/XMq4Rz7lzYWbvC5LRCQsKIyJSL3pn5XK2z8dSo+2TfjpCwv4f28tpaSs3OuyREQ8pTAmIvWqZZM4/u/aQUw+I4t/fPIN33/iC7bvO+x1WSIinlEYE5F6FxPl496LevDnCTksKtjNmP/9WBPEikijpTAmIp4Zl5vOazcOIS7az8THZvPc7LU4p+kvRKRxURgTEU91a9OEN24ayrBOLfj160v4xcsLOVhc5nVZIiL1RmFMRDyXkhDNE5P6cevIzry2YCMXP/oZ6wuLvC5LRKReKIyJSFjw+YxbRnbiyav6s3FXERf+9RNmrtjmdVkiInVOYUxEwsqIri1566fDaNs0nqufnstf3tcySiLSsCmMiUjYaZeWwKs/OYPxOen8z/tfc+2zX7LnoJZREpGGSWFMRMJSfIyfP13eh/vG9uCjr7dz0V8/YdnmvV6XJSJS6xTGRCRsmRmTBmcx9fpBHCwuY/zfPmXago1elyUiUqsUxkQk7J3ePpW3fjaU3ulN+fnUfO5+fTF7D+mypYg0DApjIhIRWibH8fy1A7lmaDbPzl7HWX+YyT8++YbDpZqTTEQim8KYiESMaL+PX4/pzls/HUqPtin8v7eWcs6fPmLago36xqWIRCyFMRGJOD3TU/jnjwfy3DUDSImP5udT8xnzv58w6+vtXpcmInLCFMZEJGIN69SCN28eyl8m5rD3UAmTnpzDD574gq8K9nhdmohIyBTGRCSi+XzG2Jx0PvjFWdw9pjtLNu3hwr9+ws9eWKAllUQkIiiMiUiDEBvl5+qh2Xz0yxHcPKIjM5Zu4ZyH8rj3jSUU7j/sdXkiItVSGBORBqVJXDT/cX4XPrp9BJeenslzn6/jrAfz+N8PVlJUXOp1eSIi36EwJiINUqsmcfz3xb149+dnMqRjGn9672vOejCP579YR0lZudfliYgcpTAmIg1ax5ZJ/P2H/fjXT84gKy2BO19bzPn/M4t/L96Mc5oOQ0S8pzAmIo3C6e2b8dL1g3liUj/8PuOGf87n4kc/Y843O70uTUQaOYUxEWk0zIyR3Vsx/ZZh/OGS3mzefYjL/z6ba56ey9db93ldnog0UlFeFyAiUt+i/D4u75/JRTlteerTtfwtbxXn/3kWZ5yWxticdEb1bE2TuGivyxSRRkIjYyLSaMVF+/nJ8NP4+JcjuOWcTmzcdZBfvrKI/r99n5uen897S7dSXKqb/UWkbmlkTEQavaYJMfx8ZGduOacT+Rt283r+Jt5cuIm3v9pM04RoLujVhnG56Zzerhk+n3ldrog0MApjIiJBZkZuu2bktmvGnRd045NVO5i2YCOvzt/I81+sJ6NZPGNz2jIuJ51OrZK9LldEGgiFMRGRKkT7fYzo0pIRXVpy4HApM5ZuYdqCTTyat5pHZq6mR9smjM9N58I+bWnVJM7rckUkgimMiYjUIDE2ivG5GYzPzWD7vsO8tWgT0xZs5LdvL+P+d5Yx5LTmjM1py6ierUnWjf8icoIUxkRETkCL5Fh+NCSbHw3JZs32/UzL38Tr+Ru5/ZVF3DVtMSO7t2J8Tjpndm5BTJS+IyUiNVMYExE5SR1aJHHbuZ25dWQnFmzYzesLNvLmos28vejbG//P6daS3MxmNEuM8bpcEQlTCmMiIqfIzOjbrhl92zXjrjHd+WTlDl5bsJF/zS/g+S/WA3Bai0ROb9/s6KND8yR9M1NEAIUxEZFaFe33MaJrS0Z0bcnB4jIWFuxm3rpdzF+3ixlLt/LSlwUApMRH07ddU05v34y+7ZuRk9mUhBj9L1mkMdLffBGROhIf42dQhzQGdUgDwDnHmh0Hjoazeet2MXPFdgD8PqNbm2RObxcIZ6e3b0Z603jMNHom0tApjImI1BMz47QWSZzWIonL+2UCsKeohPkbvg1nL88r4JnZ6wBo1SQ2MHLWLhDOerRN0ZcCRBoghTEREQ+lJEQfnc8MoLSsnOVb9jF/fSCczVu3i3e+2gJAbJSP3hkp9ExPoX1qAu2bJ9I+NYGMZglefgQROUUKYyIiYSTK76NneiBwTRqcBcDWvYeOjpx9uW4XU+duoKi47OgxPoPUOKPLqs9pl5pIVloC7dMSaJeaSPu0BBJj9b96kXCmv6EiImGuVZM4Rvdqw+hebYDAvWfb9x9mfWER6wqLWFd4gDnLvuHA4TL+vXgzu4pKjjm+eVIs7YMBrX0woAUeiTRLiNZ9aSIeUxgTEYkwZkbL5DhaJsfRLysVgLyYzQwfPgSAvYdKWF9YxNrCA6wrLDr6fPbqQl6dv/GYcyXHRtEuLYHMZgmkJcWQlhhDs8QYUhNjSEuMpVli9NH/xkb56/2zijQGIYUxMxsF/AXwA0845x6otN+C+78HFAGTnXPzQzlWRERqV5O46KOXOis7VFLGhp3BEbWdRawvPMDawiJWLVrJnGLHrrgknFX9JYGkKEiNjyI1KY7UJvGkJsWSFgxuzRJjjj4/8kiKjdKom0gIagxjZuYHHgHOBQqAuWb2hnNuaYVmo4FOwcdA4FFgYIjHiohIPYmL9tOpVTKdWiUfu+PsO2DmTMrMx564JHYmNKEwPoVdiU0pbJLGzsSm7IxPZmdcMjtjk9gal8yyxKYUxjeh2F/1epzRBknRRkK0n4TYqMAjLobEuCjiY6JIjPGTEBNFQoyfhFj/0XaJR7bF+EmMjSI+xh/YFmwT5dc3SqVhCWVkbACwyjm3BsDMXgTGAhUD1VjgWeecAz43s6Zm1gbICuFYEREJE35XTurBvaQe3EtHCmps74Ci6Dh2JqRQmJDCrvgmgfCW3Iyd8U0oiorjQHQsB/0xHIiKpSgqls0xcRTFxnMgOp6D0YH95dWMxlVZo0GMD6J9RozfiPb7iPH7iI7yBZ5H+YmJ9hMd7Q++Du4/0vZou8C2GL+fKL/h9xl+M3w+I8oX+K/fKjz3gc+MKJ/v2+d+w2ffHuv3BR5HjvWZYQZmgfZHXgcWXzB8we3bisrZsLMIAJ8vsN2C+6l4LIFzBZ4F9gVec3QUsmKbigOTlbcdOeZI+2/baTSzvoUSxtKBDRVeFxAY/aqpTXqIx4qISIQyILHkEIl7DpG5Z+tJncMBh6NiKIqO40B0HEUx8RRFx1IUE8+BmHgOxiZwICbu6P5iXxQl/iiKfVEU+6Mp8Qdf+6MpDj4v8UVTHBXFAd+326puG02ZL0zuhZs10+sKqmW4Cs/5zvOqtlV+fkLvZ676fSdynlpuV51W0PxUjg8ljFVVY+WfUnVtQjk2cAKz64Drgi8Pm9niEGqT8NMc2OF1EXLS1H+R65T6LhPaNoXUWqznxJQWBx4H99bqaY8MHFmFASMDosFiTv13cK3ZDrTwugg5aeuh/akcH0oYKwAyK7zOADaF2CYmhGMBcM49BjwGYGZfOuf6hVCbhBn1XWRT/0Uu9V1kM7Mv16n/IpaZfXkqx4dykX4u0MnMss0sBpgIvFGpzRvAJAsYBOxxzm0O8VgRERGRRqvGkTHnXKmZ3Qy8S2B6iiedc0vM7Ibg/inAOwSmtVhFYGqLHx3v2Dr5JCIiIiIRKKR5xpxz7xAIXBW3Tanw3AE3hXpsCB47wfYSPtR3kU39F7nUd5FN/RfZTqn/LJCjRERERMQLmjlPRERExEP1HsbM7Ekz21Zx6gozm2pm+cHHWjPLr7DvV2a2ysxWmNn59V2vHOtE+s/MsszsYIV9U6o/s9S1avoux8w+D/bPl2Y2oMI+/d0LIyfSf/q7F36q6b8+ZjbbzL4yszfNrEmFffr7FyZOpO9O+u+ec65eH8CZQF9gcTX7/wTcHXzeHVgIxALZwGrAX98163HS/ff/27t/F7mqMIzj38cEjZggiAoWiiuojY2IFkKIBBG0UFAUbRQCihILYycIQhQstPEPCHZJVEiRUkE0hYUhCKISELVwUSy0VtnwWsyd7CbOxZ3B7Dmz+X7gsmfv7MKBhxfeOffHuXXs7zz6yA74GHh4GD8CfDaMrb3Ojjnzs/Y6O0byOw3sG8YHgDeHsfXX0TFndgvV3pavjFXVKeCPWZ8NG44/BRwbTj0GHK+qv6rqJyZPa94363+1NebMTx0Zya6A6bfxa1l/D6C115k581NnRvK7Ezg1jD8BnhjG1l9H5sxuIb3dM7YX+K2qvh9+H9tmSX26OD+AlSRfJfk8yd5WE9OoV4B3kvwMvAu8Npy39pbDWH5g7S2Db4BHh/GTrL8k3frr31h2sEDt9daM5BuNbAAAAkBJREFUPcOFqyqb3k5JXbg4v1+BW6rqbuBV4OjGeyLUhZeAQ1V1M3AIODKct/aWw1h+1t5yOAAcTHIG2AP8PZy3/vo3lt1CtddNM5ZkJ/A48MGG05vZikkdmJXfsMT++zA+w+S+hzvazFAjngNODOOPWL8UYu0th5n5WXvLoarOVtVDVXUPky+yPwwfWX+dG8tu0drrphkDHgTOVtXqhnMngaeTXJVkBbgd+LLJ7PRf/pVfkhuS7BjGtzHJ78dG89NsvwD7hvF+YHqJ2dpbDjPzs/aWQ5Ibh59XAK8D0yfvrL/OjWW3aO1t6g38/6ckx4AHgOuTrAJvVNURJvtWXnDjd022XfoQ+A5YAw5W1bktnrI2mCc/Jk+gHE6yBpwDXqyqmTf/69KblR3wPPDesLL5J/ACWHs9mic/rL3ujOS3O8l095oTwPtg/fVmnuxYsPZ8A78kSVJDPV2mlCRJuuzYjEmSJDVkMyZJktSQzZgkSVJDNmOSJEkN2YxJkiQ1ZDMmSZLUkM2YpG0tyb1Jvk6yK8k1Sb5NclfreUnSlC99lbTtJXkL2AVcDaxW1duNpyRJ59mMSdr2klwJnGayZdD9bi0jqSdeppR0ObgO2A3sYbJCJkndcGVM0raX5CRwHFgBbqqqlxtPSZLO29l6ApJ0KSV5FlirqqNJdgBfJNlfVZ+2npskgStjkiRJTXnPmCRJUkM2Y5IkSQ3ZjEmSJDVkMyZJktSQzZgkSVJDNmOSJEkN2YxJkiQ1ZDMmSZLU0D+0OaEqPszmFgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_p_auc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Interpretation\n",
    "---\n",
    "* A small p-value (typically ≤ 0.05) indicates **strong evidence against the null hypothesis $H_0$**, so you reject the null hypothesis.\n",
    "\n",
    "* A large p-value (> 0.05) indicates weak evidence against the null hypothesis, so you fail to reject the null hypothesis.\n",
    "\n",
    "* p-values very close to the cutoff (0.05) are considered to be marginal (could go either way)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 4- Significance Level ($\\alpha$)\n",
    "---\n",
    "* It is the degree of certainty required in order to **reject** the null hypothesis $H_0$.\n",
    "* A test statistic, TS, with p-value **less** than some pre-determined false positive (or size) is said to be statistically significant at that level.\n",
    "* Commonly used p-values:\n",
    "\n",
    "|<center> P-Value &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; </center> |<center> Wording </center>|\n",
    "| -----------------------------------------| --------------- |\n",
    "|<center> $p > 0.05$ </center>| <center> Not Significant </center>|\n",
    "|<center> $0.01 \\leq p \\leq 0.05$ </center>|<center> Significant</center>|\n",
    "|<center> $0.001 \\leq p < 0.01$</center>| <center> Very Significant </center>|\n",
    "|<center> $p < 0.001$</center>| <center> Extremely Significant </center>|\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/cotton/64/000000/tie.png\" style=\"height:30px;display:inline\"> Formalization\n",
    "---\n",
    "Let's design a level $\\alpha$ test to choose between: $$ H_0: \\mu= \\mu_0 $$ $$ H_1: \\mu \\neq \\mu_0$$\n",
    "We initially assume $H_0$, thus $z_{stat} \\sim \\mathcal{N}(0,1)$. We will choose a threshold $c$. If $|z_{stat}| \\leq c$, we **accept** $H_0$, and if $|z_{stat}| > c$, accept $H_1$.\n",
    "\n",
    "To choose $c$: $$ P(|z_{stat}| >c | H_0) = \\alpha $$\n",
    "Since the standard normal PDF is **symmetric** around $0$, we have: $$ P(|z_{stat}| >c | H_0) = 2P(z_{stat} > c| H_0) \\rightarrow P(z_{stat} > c| H_0) = \\frac{\\alpha}{2} \\rightarrow c=z_{\\frac{\\alpha}{2}}$$\n",
    "Therefore, we accept $H_0$ if $$ \\mid\\frac{\\overline{X} - \\mu_0}{\\frac{\\sigma}{\\sqrt{n}}}\\mid \\leq z_{\\frac{\\alpha}{2}} $$ and reject it otherwise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/trust.png\" style=\"height:30px;display:inline\"> Relation to Confidence Intervals\n",
    "---\n",
    "Notice that saying we accept $H_0$ if  $$ \\mid\\frac{\\overline{X} - \\mu_0}{\\frac{\\sigma}{\\sqrt{n}}}\\mid \\leq z_{\\frac{\\alpha}{2}} $$ can be interpreted as the following **acceptance region** for $\\mu_0$: $$ \\mu_0 \\in \\big[ \\overline{X}-  z_{\\frac{\\alpha}{2}}\\cdot \\frac{\\sigma}{\\sqrt{n}}, \\overline{X}+  z_{\\frac{\\alpha}{2}}\\cdot \\frac{\\sigma}{\\sqrt{n}}\\big] $$ Which is the $(1-\\alpha)$ confidence interval for $\\mu_0$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/task.png\" style=\"height:30px;display:inline\"> Exercise - Hypothesis Testing\n",
    "---\n",
    "We continue with the radar example. Recall that the system receives a signal and based upon that signal it decides whether an aircraft is present or not. We denote:\n",
    "* $X$ - the received signal (R.V., sampled)\n",
    "* We suppose that: $$ X=W, \\text{   if no aircraft is present} $$ $$ X=1+W, \\text{   if an aircraft is present} $$ $$ W \\sim \\mathcal{N}(0,\\sigma^2=\\frac{1}{9}) $$\n",
    "* We can write instead: $$X = \\theta + W$$ where $\\theta=0$ if there is no aircraft and $\\theta=1$ otherwise.\n",
    "* The **hypotheses**: \n",
    "    * $H_0$ (null): No aircraft is present\n",
    "    * $H_1$ (alternative): An aircraft is present"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "1. Write $H_0$ and $H_1$ in terms of possible values of $\\theta$.\n",
    "2. Suggest a *simple* test statistic with level $\\alpha = 0.05$ to decide between $H_0$ and $H_1$.\n",
    "3. Find the probability of missing a present aircraft, that is, find $\\beta$ (the probability of type 2 error).\n",
    "4. If we observe $X=0.6$, is there enough evidence to reject $H_0$ at a *significance level* $\\alpha=0.01$?\n",
    "5. For a probability less than 5% to miss a present aircraft, what is the **smallest** significance level that we can achieve?\n",
    "\n",
    "\n",
    "* Reminder: $$W \\sim \\mathcal{N}(\\mu, \\sigma^2) \\rightarrow \\frac{W - \\mu}{\\sigma} \\sim \\mathcal{N}(0,1)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/idea.png\" style=\"height:30px;display:inline\"> Solution\n",
    "---\n",
    "#### Section 1 - Write $H_0$ and $H_1$ in terms of possible values of $\\theta$.\n",
    "---\n",
    "The **hypotheses**: \n",
    "   * $H_0$ (null): No aircraft is present: $\\theta=0$\n",
    "   * $H_1$ (alternative): An aircraft is present : $\\theta=1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Section 2 - Suggest a *simple* test statistic with level $\\alpha = 0.05$ to decide between $H_0$ and $H_1$.\n",
    "---\n",
    "* The *observed* data is a random variable $X$. \n",
    "* Under $H_0$ we have $X \\sim \\mathcal{N}(0, \\frac{1}{9})$ and under $H_1$ we have $X \\sim \\mathcal{N}(1, \\frac{1}{9})$. \n",
    "* We suggest the following test: set a threshold $c$. If the obsereved value of $X$ is less than $c$, choose $H_0$ ($\\theta=\\mathbb{E}[X]=0$), otherwise, choose $H_1$ ($\\theta=\\mathbb{E}[X]=1$). \n",
    "* To find the best $c$ we use the required $\\alpha$, that is, we demand: $$ P(\\text{type I error}) = P(\\text{Reject } H_0 | H_0) = P(X>c|H_0) = P(W>c) =P(\\frac{Z}{3} > c) = 1 - \\phi(3c) =\\alpha$$ (the last transition is due to the fact that we assume $H_0$ and $X \\sim \\mathcal{N}(0, \\frac{1}{9})$, which is not the standard distribution). \n",
    "* It holds that $P(\\text{type I error}) = \\alpha$, thus we get: $$ c=\\frac{1}{3}\\phi^{-1}(1- \\alpha) = \\frac{1}{3}\\phi^{-1}(1-0.05) = \\frac{1}{3}\\phi^{-1}(0.95) = 0.548 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Section 3 - Find the probability of missing a present aircraft, that is, find $\\beta$ (the probability of type 2 error).\n",
    "---\n",
    "Note that the *alternative* hypothesis is simple, that is, it contains only one value ($\\theta = 1$), so $\\beta(\\theta) = \\beta$ and we can write: $$ \\beta = P(\\text{type II error}) = P(\\text{Accept } H_0 | H_1) = P(X < c | H_1) = P(1+W < c) = P(W<c-1) = \\phi(3(c-1)) $$ Since we found out that for the given $\\alpha$, $c=0.548$ then $\\beta = 0.088$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Section 4 - If we observe $X=0.6$, is there enough evidence to reject $H_0$ at a *significance level* $\\alpha=0.01$?\n",
    "---\n",
    "* For $\\alpha=0.01$ we get $ c=\\frac{1}{3}\\phi^{-1}(1- \\alpha) = \\frac{1}{3}\\phi^{-1}(1-0.01) = \\frac{1}{3}\\phi^{-1}(0.99) = 0.775$, which is **larger** than $0.6$. \n",
    "* Thus, we **cannot** reject $H_0$ at significance level $\\alpha=0.01$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Section 5 - For a probability less than 5% to miss a present aircraft, what is the **smallest** significance level that we can achieve?\n",
    "---\n",
    "We want $\\beta = 0.05$, and from (3) we deduce that $c=1 + \\frac{1}{3}\\phi^{-1}(\\beta)=0.452$. Thus, we need $c \\leq 0.452$ to obtain $\\beta \\leq 0.05)$. Let's calculate $\\alpha$: $$ P(\\text{type I error}) = P(\\text{Reject } H_0 | H_0) = 1 - \\phi(3c) = 0.0875 $$ which means that the smallest significance level that we can achieve is 0.0875."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/nolan/64/000000/table-1.png\" style=\"height:30px;display:inline\"> Hypothesis Testing for the Mean Summary\n",
    "---\n",
    "All expansions can be found <a href=\"https://www.probabilitycourse.com/chapter8/8_4_3_hypothesis_testing_for_mean.php\">HERE</a>.\n",
    "\n",
    "* 2-sided hypothesis testing for the mean: $H_0: \\mu=\\mu_0, H_1: \\mu \\neq \\mu_0$\n",
    "\n",
    "| <center>Case </center>  | Test Statistic | Acceptance Region |\n",
    "| ------ | --- | --- |\n",
    "| $X_i \\sim \\mathcal{N}(\\mu, \\sigma)$, $\\sigma$ **known** &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;| $W = \\frac{\\overline{X}  -\\mu_0}{\\frac{\\sigma}{\\sqrt{n}}}$| $\\mid W \\mid \\leq z_{\\frac{\\alpha}{2}} $  |\n",
    "|$n$ large, $X_i$ non-normal &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;| $W = \\frac{\\overline{X}  -\\mu_0}{\\frac{S}{\\sqrt{n}}}$| $\\mid W\\mid \\leq z_{\\frac{\\alpha}{2}} $  |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "* 1-sided hypothesis testing for the mean: $H_0: \\mu \\leq \\mu_0, H_1: \\mu > \\mu_0$\n",
    "\n",
    "| <center>Case </center>  | Test Statistic | Acceptance Region |\n",
    "| ------ | --- | --- |\n",
    "| $X_i \\sim \\mathcal{N}(\\mu, \\sigma)$, $\\sigma$ **known** &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;| $W = \\frac{\\overline{X}  -\\mu_0}{\\frac{\\sigma}{\\sqrt{n}}}$| $W \\leq z_{\\alpha} $  |\n",
    "|$n$ large, $X_i$ non-normal &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;| $W = \\frac{\\overline{X}  -\\mu_0}{\\frac{S}{\\sqrt{n}}}$| $W \\leq z_{\\alpha} $  |\n",
    "\n",
    "* The only difference is the *absolute* sign on $W$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/bubbles/50/000000/video-playlist.png\" style=\"height:50px;display:inline\"> Recommended Videos\n",
    "---\n",
    "#### <img src=\"https://img.icons8.com/cute-clipart/64/000000/warning-shield.png\" style=\"height:30px;display:inline\"> Warning!\n",
    "* These videos do not replace the lectures and tutorials.\n",
    "* Please use these to get a better understanding of the material, and not as an alternative to the written material.\n",
    "\n",
    "#### Video By Subject\n",
    "\n",
    "* Hypothesis Testing - <a href=\"https://www.youtube.com/watch?v=VK-rnA3-41c\">Hypothesis Testing Statistics Problems & Examples</a>\n",
    "* p-Value - <a href=\"https://www.youtube.com/watch?v=eyknGvncKLw\">Understanding the p-value - Statistics Help</a>\n",
    "    * <a href=\"https://www.youtube.com/watch?v=-MKT3yLDkqk\">What is a P Value? What does it tell us?</a>\n",
    "* Test Statistics (t-stat is covered in the next tutorial):\n",
    "    * <a href=\"https://www.youtube.com/watch?v=QZ7kgmhdIwA\">Test Statistics: Crash Course Statistics</a>\n",
    "    * <a href=\"https://www.youtube.com/watch?v=5ABpqVSx33I\">Z-statistics vs. T-statistics</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <img src=\"https://img.icons8.com/dusk/64/000000/prize.png\" style=\"height:50px;display:inline\"> Credits\n",
    "---\n",
    "* Examples, exercises and definitions from <a href=\"https://probabilitycourse.com/\">Introduction to Probability, Statistics and Random Processes</a> - https://probabilitycourse.com\n",
    "* Icons from <a href=\"https://icons8.com/\">Icon8.com</a> - https://icons8.com\n",
    "* Datasets from <a href=\"https://www.kaggle.com/\">Kaggle</a> - https://www.kaggle.com/"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}