week06/2021F_Week06.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Week06\n",
    "\n",
    "Goals:\n",
    "- Write down the variance for $\\hat \\beta$\n",
    "- Write down the distribution for $\\hat \\beta$\n",
    "- Use the $\\hat \\beta$ distribution to compute confidence intervals\n",
    "- Use our estimates of $\\beta$ to perform hypothesis tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variance, $\\mathcal{I}^{-1}$, for $\\hat \\beta_{0}$ and $\\hat \\beta_{1}$\n",
    "\n",
    "We found that given a dataset of $N$ pairs of points ($x,y$) the MLE for $\\hat \\beta_{1}$ was \n",
    "\n",
    "\\begin{align}\n",
    "    \\hat \\beta_{1} = \\frac{\\sum_{i=1}^{N} (x_{i} - \\bar x) (y_{i} - \\bar y)}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2}  }\n",
    "\\end{align}\n",
    "\n",
    "To compute the variance for \\hat \\beta_{1} we need to consider repeated sampels of our datset from random variables ($X_{i},Y_{i}$).\n",
    "If we shift our focus from a fixed dataset to random variables that can generate our fixed data then the MLE is a funciton of random variables.\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat \\beta_{1} = \\frac{\\sum_{i=1}^{N} (X_{i} - \\bar X) (Y_{i} - \\bar Y)}{\\sum_{i=1}^{N} (X_{i} - \\bar X)^{2}  }\n",
    "\\end{align}\n",
    "\n",
    "Above we use capital $X$ and capital $Y$ to represent random variables. \n",
    "To be clear, $\\text{Var}(y_{i}) = 0$ and $\\text{Var}(Y_{i}) \\ge 0$.\n",
    "\n",
    "#### Variance of $\\hat \\beta_{1}$\n",
    "\n",
    "Before we compute the variance of $\\hat \\beta_{1}$ we can simplify the numerator of the MLE\n",
    "\\begin{align}\n",
    "    &\\sum_{i=1}^{N} (X_{i} - \\bar X) (Y_{i} - \\bar Y) \\\\ \n",
    "    &\\sum_{i=1}^{N} (X_{i} - \\bar X) Y_{i} -  (X_{i} - \\bar X) \\bar Y \\\\ \n",
    "    &\\sum_{i=1}^{N} (X_{i} - \\bar X) Y_{i} - \\sum_{i=1}^{N} (X_{i} - \\bar X) \\bar Y \\\\\n",
    "    &\\sum_{i=1}^{N} (X_{i} - \\bar X) Y_{i} -  \\bar Y \\sum_{i=1}^{N} (X_{i} - \\bar X) & \\bar Y\\text{ is a constant } \\\\\n",
    "    &\\sum_{i=1}^{N} (X_{i} - \\bar X) Y_{i} & \\text{the second sum is zero}\n",
    "\\end{align}\n",
    "\n",
    "Our MLE for $\\hat \\beta_{1}$ is now simpler\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat \\beta_{1} = \\frac{\\sum_{i=1}^{N} (X_{i} - \\bar X) Y_{i}}{\\sum_{i=1}^{N} (X_{i} - \\bar X)^{2}  }\n",
    "\\end{align}\n",
    "\n",
    "We can think of the above as \n",
    "\n",
    "\\begin{align}\n",
    "    \\hat \\beta_{1} = \\sum_{i=1}^{N} k_{i}Y_{i}\n",
    "\\end{align}\n",
    "\n",
    "where we now consider\n",
    "\n",
    "\\begin{align}\n",
    "    k_{i} = \\frac{(X_{i} - \\bar X)}{\\sum_{i=1}^{N} (X_{i} - \\bar X)^{2} }\n",
    "\\end{align}\n",
    "\n",
    "a constant. That is, we assume we are given $x_{i}$ values\n",
    "\n",
    "\\begin{align}\n",
    "    k_{i} = \\frac{(x_{i} - \\bar x)}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} }\n",
    "\\end{align}\n",
    "\n",
    "Our goal is to find the variance of $\\hat \\beta_{1}$ or \n",
    "\\begin{align}\n",
    "    Var \\left(\\hat \\beta_{1}\\right) = Var \\left( \\sum_{i=1}^{N} k_{i}Y_{i} \\right)\n",
    "\\end{align}\n",
    "\n",
    "#### The variance of the sum (or differnece) of random variables\n",
    "We will need a fact about the variance of the sum of independent random variables.\n",
    "\n",
    "For two independent random variables $X$ and $Y$ and constants $a$ and $b$ the following is true\n",
    "\n",
    "\\begin{align}\n",
    "    Var(aX+bY) = a^{2}Var(X) + b^{2}Var(Y)\n",
    "\\end{align}.\n",
    "\n",
    "This also means that \n",
    "\n",
    "\\begin{align}\n",
    "    Var(X-Y) = Var((1)X + (-1)Y) =1^{2} Var(X) + (-1)^{2} Var(Y) = Var(X) +  Var(Y) \n",
    "\\end{align}.\n",
    "\n",
    "\n",
    "The same relationship holds for $N$ random variables $Y_{1},Y_{2},\\cdots,Y_{N}$ which are all independent from one another and for constants $k_{1},k_{2},\\cdots,k_{N}$\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "    Var(a_{1}Y_{1}+a_{2}Y_{2}+\\cdots + a_{N}Y_{N}) = Var\\left(\\sum_{i=1}^{N}a_{i}Y_{i}\\right) = \\sum_{i=1}^{N} a_{i}^{2} Var\\left(Y_{i}\\right)  \n",
    "\\end{align}.\n",
    "\n",
    "This mean then \n",
    "\n",
    "\\begin{align}\n",
    "    Var \\left(\\hat \\beta_{1}\\right) &= Var \\left( \\sum_{i=1}^{N} k_{i}Y_{i} \\right) \\\\ \n",
    "    &= \\sum_{i=1}^{N} k^{2}_{i} Var(Y_{i})\n",
    "\\end{align}\n",
    "\n",
    "However, our linear regression model states that every $Y_{i}$ has the following same distribution\n",
    "\\begin{align}\n",
    "    Y_{i} | x_{i} \\sim \\mathcal{N}\\left( \\beta_{0} + \\beta_{1}x_{i}, \\sigma^{2} \\right)\n",
    "\\end{align}\n",
    "\n",
    "That is, the variance for each $Y_{i}$ is $\\sigma^{2}$.\n",
    "Then the above equation simplifies to \n",
    "\n",
    "\\begin{align}\n",
    "    Var \\left(\\hat \\beta_{1}\\right) &= \\sum_{i=1}^{N} k^{2}_{i} \\sigma^{2}\\\\\n",
    "    &=\\sigma^{2} \\sum_{i=1}^{N} k^{2}_{i} \\\\\n",
    "\\end{align}\n",
    "\n",
    "Our final step is to figure out what $\\sum_{i=1}^{N} k^{2}_{i}$ equals. \n",
    "\n",
    "\\begin{align}\n",
    "    k_{i} = \\frac{(x_{i} - \\bar x)}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} }\n",
    "\\end{align}\n",
    "\n",
    "and if we square $k_{i}$ we are left with\n",
    "\n",
    "\\begin{align}\n",
    "    k_{i}^{2} = \\frac{ (x_{i} - \\bar x)^{2} }{\\left(\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2}\\right)^{2} }\n",
    "\\end{align}\n",
    "\n",
    "Summing up the $k_{i}^{2}$s we find \n",
    "\n",
    "\\begin{align}\n",
    "    \\sum_{i=1}^{N} k_{i}^{2} &= \\sum_{i=1}^{N} \\frac{ (x_{i} - \\bar x)^{2} }{\\left(\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2}\\right)^{2} }\\\\\n",
    "     &=\\frac{1}{\\left(\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2}\\right)^{2} } \\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} \\\\\n",
    "     &= \\frac{1}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} }\n",
    "\\end{align}\n",
    "\n",
    "#### Variance of $\\hat \\beta_{1}$\n",
    "The final result is the variance for $\\hat \\beta_{1}$\n",
    "\\begin{align}\n",
    "    Var(\\hat \\beta_{1}) = \\sigma^{2} \\frac{1}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} }\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Variance of $\\hat \\beta_{0}$\n",
    "\n",
    "To compute the variance for $\\hat beta_{0}$ we can proceed the same way as we did for $\\hat \\beta_{1}$: (i) consider the random variables that generated the data, (ii) take the variance of the MLE (now a random variable), and (iii) simplify. \n",
    "\n",
    "The MLE for $\\hat \\beta_{0}$ is \n",
    "\n",
    "\\begin{align}\n",
    "    \\hat \\beta_{0} = \\bar Y - \\hat \\beta_{1} \\bar x\n",
    "\\end{align}\n",
    "\n",
    "The variance is then \n",
    "\n",
    "\\begin{align}\n",
    "    Var (\\hat \\beta_{0}) &= Var(\\bar Y) + \\bar{x}^{2} Var(\\hat \\beta_{1}) &\\text{ } \\bar x \\text{ is a constant} \\\\ \n",
    "    &= Var\\left(\\frac{\\sum_{i=1}^{N}Y_{i}}{N}\\right) + \\bar{x}^{2} Var(\\hat \\beta_{1}) &\\text{ }  \\text{ Definition of sample average} \\\\ \n",
    "    &= \\frac{1}{N^{2}}\\sum_{i=1}^{N}Var(Y_{i}) + \\bar{x}^{2} Var(\\hat \\beta_{1}) & \\frac{1}{N}  \\text{is a constant} \\\\ \n",
    "    &= \\frac{1}{N^{2}}\\sum_{i=1}^{N}\\sigma^{2} + \\bar{x}^{2} Var(\\hat \\beta_{1}) & \\text{The variance of every } Y_{i} \\text{ is } \\sigma^{2} \\\\ \n",
    "    &= \\frac{1}{N^{2}}N\\sigma^{2} + \\bar{x}^{2} Var(\\hat \\beta_{1})  \\\\\n",
    "    &= \\frac{1}{N^{2}}N\\sigma^{2} + \\bar{x}^{2} \\sigma^{2} \\frac{1}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} }  \\\\\n",
    "    &= \\sigma^{2}\\frac{1}{N} +  \\sigma^{2} \\frac{\\bar{x}^{2}}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} }  \\\\\n",
    "    &= \\sigma^{2}\\left(\\frac{1}{N}  + \\frac{\\bar{x}^{2}}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} }\\right)  \\\\\n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The distribution for $\\hat \\beta$\n",
    "\n",
    "We know that as our dataset grows that the MLE for $\\hat \\beta_{0}$ and for $\\hat \\beta_{1}$ approaches a normal distribution centered at our true parameter values with mean parameter equal to the correpsonding MLEs and variance parameter equal to the corresponding variances we found above.\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat \\beta_{0} \\sim \\mathcal{N}\\left( \\bar Y - \\hat \\beta_{1}\\bar x , \\sigma^{2}\\left(\\frac{1}{N}  + \\frac{\\bar{x}^{2}}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} }\\right)  \\right)\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat \\beta_{1} \\sim \\mathcal{N}\\left( \\sum_{i=1}^{N} \\frac{(x_{i}-\\bar x)y_{i} }{ (x_{i}-\\bar x)^{2} }, \\sigma^{2} \\frac{1}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} } \\right)   \n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing confidence intervals for $\\hat \\beta$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equation to compute a $1-\\alpha$ confidence interval for a normally distributed random variable with mean equal to $\\mu$  and variance equal to $\\sigma^{2}$ is \n",
    "\n",
    "\\begin{align}\n",
    "    \\text{CI}_{1-\\alpha} = ( \\mu - z_{1-\\alpha/2}\\sigma,\\mu + z_{\\alpha/2}\\sigma  )\n",
    "\\end{align}\n",
    "\n",
    "where $z_{\\alpha/2}$ is the $\\alpha/2$ quantile value from a random variable $Z\\sim\\mathcal{N(0,1)}$. \n",
    "The $p$ quantile value $q$ is the value $q$ such that $P(Z<p) = q$\n",
    "\n",
    "\n",
    "| $\\alpha$    | $z_{1-\\alpha/2}$ |\n",
    "| ----------- | -----------------|\n",
    "| 0.01        | 2.58             |\n",
    "| 0.05        | 1.96             |\n",
    "| 0.10        | 1.64             |\n",
    "| 0.20        | 1.28             |\n",
    "\n",
    "\n",
    "### Why this formula works (Shifting and scaling Normal distributions)\n",
    "\n",
    "For a random variable $X\\sim \\mathcal{N}\\left(\\mu, \\sigma^{2}\\right)$\n",
    "We can derive the distribution of a second random variable $Y = X + c$ where $c$ is a constant. \n",
    "Adding or subtracting a constant from a normally distributed random variable,X,---shifting---creates a new random variable $Y$ that is normally distributed \n",
    "\n",
    "\\begin{align}\n",
    "    X &\\sim \\mathcal{N}\\left(\\mu,\\sigma^{2}\\right) \\\\ \n",
    "    \\text{if  }\n",
    "    Y &= X + c\\\\\n",
    "    \\text{then  }\n",
    "    Y &\\sim \\mathcal{N}\\left(\\mu+c,\\sigma^{2}\\right)\n",
    "\\end{align}\n",
    "\n",
    "Multiplying a random variable $X$ by a constant $c$---scaling---creates a new random variable $Y$ that is normally distributed\n",
    "\n",
    "\\begin{align}\n",
    "    X &\\sim \\mathcal{N}\\left(\\mu,\\sigma^{2}\\right) \\\\ \n",
    "    \\text{if  }\n",
    "    Y &= cX\\\\\n",
    "    \\text{then  }\n",
    "    Y &\\sim \\mathcal{N}\\left(\\mu,c^{2}\\sigma^{2}\\right)\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "### From a $N(0,1)$ to any $N(\\mu,\\sigma^{2})$\n",
    "  \n",
    "Let $Z \\sim \\mathcal{N}\\left(0,1\\right)$. \n",
    "Then we can transform $Z$ into $Y$, a $\\mathcal{N}\\left(0,\\sigma^{2}\\right)$, by mutliplying by the constant $\\sigma$\n",
    "\n",
    "\\begin{align}\n",
    "  \\text{if  } Y = \\sigma Z \\text{ then } Y\\sim \\mathcal{N}\\left(0,\\sigma^{2}\\times 1\\right)\\\\\n",
    "  Y \\sim \\mathcal{N}\\left(0,\\sigma^{2}\\right)\n",
    "\\end{align}\n",
    "\n",
    "Then we can create a normal distribution, $X$, centered at $\\mu$ by shifting $Y$ by the constant $\\mu$\n",
    "\n",
    "\\begin{align}\n",
    "  \\text{if  } X = Y+\\mu \\text{ then } Y\\sim \\mathcal{N}\\left(\\mu,\\sigma^{2}\\right)\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "We transformed $Z$, a $\\mathcal{N}(0,1)$ random variable, into $X$, a $\\mathcal{N}(\\mu,\\sigma^{2})$ random variable by first multiplying $Z$ by $\\sigma$ and then adding $\\mu$ \n",
    "\n",
    "\\begin{align}\n",
    "    X = \\mu + \\sigma Z\n",
    "\\end{align}\n",
    "\n",
    "This also means that any $X \\sim \\mathcal{N}\\left(\\mu,\\sigma^{2}\\right)$ can be converted to a $Z\\sim \\mathcal{N}(0,1)$ by running the above process in reverse\n",
    "\n",
    "\\begin{align}\n",
    "    Z = \\frac{X - \\mu}{\\sigma}\n",
    "\\end{align}\n",
    "\n",
    "### Confidence intervals\n",
    "\n",
    "We will compute a $1-\\alpha$ confidence interval for a random variable $X$ by computing two quantiles values: (i) the quantile value $l$ such that $P(X<l) = \\alpha/2$ and the quantile value $u$ such that $P(X<u) = 1-\\alpha/2$. \n",
    "\n",
    "We can leverage the above transformation to make computing these two quantiles easier **as long as our random variable is normally distributed**. The idea is to compute the quantile values for $Z$ and transform those values so that they are confidence intervals for an arbitrary $\\mathcal{N}(\\mu,\\sigma^{2})$ random variable.\n",
    "\n",
    "For the lower quantile value:\n",
    "\\begin{align}\n",
    "    P(Z < l) &= \\alpha/2\\\\\n",
    "    P\\left(\\frac{X-\\mu}{\\sigma} < l\\right) &= \\alpha/2\\\\\n",
    "    P(X-\\mu < l\\sigma) &= \\alpha/2\\\\\n",
    "    P(X < \\mu + l\\sigma) &= \\alpha/2\\\\\n",
    "\\end{align}\n",
    "\n",
    "For the upper quantile value\n",
    "\\begin{align}\n",
    "    P(Z < u) &= 1 - \\alpha/2\\\\\n",
    "    P\\left(\\frac{X-\\mu}{\\sigma} < u\\right) &= 1-\\alpha/2\\\\\n",
    "    P(X-\\mu < u\\sigma) &= 1-\\alpha/2\\\\\n",
    "    P(X < \\mu + u\\sigma) &= 1-\\alpha/2\\\\\n",
    "\\end{align}\n",
    "\n",
    "One final step is recognizing that because the Normal distribution is symmetric the quantile value $l$ and $u$ are the same value but opposite sign."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1b60eb310>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABExElEQVR4nO3deVwW1f7A8c+XXQFxAxdQFnNHwcQ1U9zKLTPbrEy9Wdliy+1X3Uqrm2a3tLre1puVaabXFlu0LLVMrVwSd8UdUXBBREFQ2c/vj3k0VJAHeGDmeTjv1+u8Hp55zsx8B5gvh5kz54hSCk3TNM11uZkdgKZpmla5dKLXNE1zcTrRa5qmuTid6DVN01ycTvSapmkuzsPsAC5Vv359FRYWZnYYmqZZ2e7dxmvLlubGYSEbNmw4oZQKLO4zyyX6sLAw4uLizA7DGpKSjNcmTcyNQ9OsJjbWeF2xwswoLEVEDpb0meUSvVbE3Xcbr/qXWdO0CtCJ3somTjQ7Ak2zJn1ulIlO9FbWr5/ZEWiaNelzo0x0rxsrS0gwiqZpF9u82SiaXXSL3sruucd41dfoNe1ijz9uvOpzwy52tehFZICI7BaRfSLyzBXq3SwiSkRiiix71rbebhG53hFBVxsvvWQUTdO0Cii1RS8i7sC7QH8gGVgvIguVUvGX1PMHHgPWFVnWBhgBtAUaAz+LSAulVIHjDsGF9epldgSaprkAey7ddAb2KaUSAERkPnAjEH9JvcnAa8BTRZbdCMxXSuUAB0Rkn217ayoaeLWgHwox3ckzuWxJSifp1Fkyzubh5iYE+nkTEehLZHAAPp7uZoeoaaWyJ9EHA0lF3icDXYpWEJGrgSZKqR9E5KlL1l17ybrBl+5ARO4H7gdo2rSpfZFXB+PGGa/6OmSVyjiXx9cbk/l642G2H8mgpCkbfDzduKZZfW7v1IS+rRvg7iZVG6im2anCN2NFxA14ExhT3m0opWYAMwBiYmL0TCjnvfKK2RFUK2dy8pmxKoEPf0vgbG4B7UMCeKJfC7pE1COsXk1q1/SiUCmOn85hT0omv+87wZIdx7h/zgZC69XkqetbMrhdI0R0wq90+twoE3sS/WGg6DP4IbZl5/kDkcAK2y94Q2ChiAy1Y13tSrp3NzuCamPF7uP8Y8FWUk7nMKhdQx6KvYrI4IBi6zatV5Om9WrSr00DJg5uzdL4FN76ZS/j521iTvhBXr81iiZ1a1bxEVQz+twoEyltKkER8QD2AH0xkvR64E6l1I4S6q8AnlRKxYlIW2AexnX5xsAvQPMr3YyNiYlReqwbm+3bjdfISHPjcGF5BYW8/H08s9ccpHmQH6/e3J6OoXXKvJ2CQsXn65P41+KdKOCV4e0YGtXY8QFrhtWrjVed8C8QkQ1KqZjiPiu1Ra+UyheR8cASwB2YqZTaISKTgDil1MIrrLtDRL7AuHGbDzyse9yUwfjxxqu+Rl8p0s/m8tDcjazen8Y914Tz9ICW5b656u4m3NmlKdc2r8/fP9/Mo//bxL6UTP7ev4W+lFMZnnvOeNXnhl1KbdFXNd2iL2L9euO1Uydz43BBRzPOcdeH60g+dY5/DW/HzR1DHLbt3PxCJn67jS/ikhkW3Zhpt0bh6a4fQncoPXrlZSrUotdMpBN8pUg+dZY7P1zHyTO5zL2vC53C6jp0+14ebrx2c3tC6/kybclu8goV/7k9Gg+d7DWT6ERvZefH8oiONjMKl5JyOpsRM9aScS6Pz+7tQnST2pWyHxHh4d5X4eXuxpTFOxHgPyM66C6Ymil0orcyPZ6HQ2Vm5zHmk/WcOpPLvPu6ElVJSb6o+3pGUKgU//pxF4H+3rx4Q9tK36emXUoneiubPt3sCFxGXkEhD83dyJ6UTGaO6VQlSf68cb2akXI6h5l/HKBJnZrc0yO8yvbtsvS5USY60VuZvmTjMP9avIvf9p5g6s3t6dWi2Gk1K9WEwa05nH6WyT/EExHoS2zLoCqPwaXoc6NM9N0hK1u//q+eN1q5LdxyhJl/HGBM9zBu62TO/LvubsL02zvQsoE/j3++mcPp50yJw2X8/LNRNLvoRG9lTz1lFK3c9qRk8o+vthITWofnBrU2NZYaXu68P7Ij+QWKh+ZuJCdfP1JSbi+/bBTNLjrRW9k77xhFK5fsvALGz9uIr7cH7951NV4e5v+6h9f35fVb27MlKZ1Xf9xldjhaNWH+b75WsshIPfxBBUxbsps9KVm8fmt7GtTyMTucCwZENmJM9zA++SORP/adMDscrRrQid7KVq/+a0wPrUz+2HeCj38/wKhuoZa88fmPAa2ICPTlyS+3kHEuz+xwNBenE72VPffcX2N6aHY7nZ3Hk19uISLQl2cHmntdviQ1vNz5923RHM/M4aWFxY4PqGkOo7tXWtkHH5gdgVOa9tNujp3O5puHrqGGl3VngIpqUpuHe1/FW7/s5YaoxvRuZb3/PCxLnxtlolv0VtaypZ5GsIw2HDzFZ+sOMrpbWKUNb+BI43tfxVVBfjz/3XbO5epeOHbT50aZ6ERvZStXGkWzS15BIc99vY2GtXx48nrnSAJeHm5MGRZJ8qlz/OeXvWaH4zwWLTKKZhd96cbKXnzReNVj3dhlxqoEdqdk8uGoGPy8nedXu0tEPW6LCeGj3xIY1qExrRrWMjsk63vjDeP1hhvMjcNJ6Ba9lc2caRStVEfSz/H28r1c37YB/ds0MDucMnt2YGtq1fBkwjfbsdocEZrzsyvRi8gAEdktIvtE5JliPn9ARLaJyGYR+V1E2tiWh4nIOdvyzSLyX0cfgEuLiDCKVqqpP+2iUMHEwW3MDqVc6vh68Y8BLdlw8BSLth41OxzNxZSa6EXEHXgXGAi0Ae44n8iLmKeUaqeUigamAm8W+Wy/UiraVh5wUNzVgx7Pwy4bD53i281HuO/acKeelPuWjk1o06gWry7eSXaevjGrOY49LfrOwD6lVIJSKheYD9xYtIJS6nSRt76A/t/TEfR4HqUqLFS8tCieIH9vHoq9yuxwKsTdTXh+SBuOZGTz0W8JZoejuRB77lgFA0lF3icDXS6tJCIPA08AXkCfIh+Fi8gm4DQwUSn1WzHr3g/cD9C0aVO7g3d5c+aYHYHlfbv5MFuS0nnj1ih8negGbEm6NavHgLYNeW/Ffm6NaWKpoRssRZ8bZeKwm7FKqXeVUs2AfwATbYuPAk2VUh0w/gjME5HLuhQopWYopWKUUjGBgVU/VrhlNWliFK1Y2XkFTP1pN1EhAdzUIdjscBzm2UGtyC9QTFuy2+xQrEufG2ViT6I/DBT9jobYlpVkPjAMQCmVo5RKs329AdgPtChXpNXRTz8ZRSvWp2sSOXY6m2cHtcbNheZiDa3ny5hrwliwMZk9KZlmh2NNn39uFM0u9iT69UBzEQkXES9gBLCwaAURaV7k7WBgr215oO1mLiISATQH9MVHe736qlG0y5zOzuO9Ffvp2SKQrhH1zA7H4R7s1Qw/Lw/eWKpb9cV6/32jaHYp9aKmUipfRMYDSwB3YKZSaoeITALilFILgfEi0g/IA04Bo22r9wQmiUgeUAg8oJQ6WRkH4pLmzzc7Asv6aFUC6WfzeNpJnoAtqzq+XtzXM4I3l+1hS1J6lc5xq7keu+5eKaUWA4svWfZCka8fK2G9BcCCigRYrTVsaHYElnQiK4ePfj/A4HaNiAwOMDucSnNPj3BmrU7k9aW7mTP2sv4PmmY3/WSslenxPIr17q/7yMkv5InrXPt2j5+3Bw/FNuO3vSdYsz/N7HA0J6YTvZW98cZfY3pogDHUwdy1h7jl6hCaBfqZHU6lG9k1lEYBPry+dLceGkErN+fveOzKvvrK7Ags5/0V+1EoHu3XvPTKLsDH051H+zbn2a+3sWJPKr0tOFuWKfS5USa6RW9l9esbRQPgWEY2n69P4paOIQTXrmF2OFXm/PG+/cte3ao/T58bZaITvZV9/bVRNAA+WLWfAqWcfqiDsvJ0d+OB2GZsPJTOan2t3jBrllE0u+hEb2VvvWUUjdTMHOatO8RNHYKdeuCy8rq1YwgNannz9nI9OQmgE30Z6URvZd99ZxSNj35LIK+gkId7V6/W/Hk+nu7c37MZaxNOsj5RP4qilY1O9FYWEGCUau7kmVzmrD3I0KjGhNf3NTsc09zZuSn1/bx4S085qJWRTvRWpsfzAODj3xM4l1fA+D7VszV/Xg0vd+69NoLf9p5gc1K62eFoTkQneivT43mQmZ3Hp6sPMiiyEVcF+ZsdjulGdg2ldk1P3tateq0MdD96K1u8uPQ6Lu5/fx4iMyefB3o1MzsUS/Dz9mBM9zCm/7yXPSmZtGhQTf/46XOjTHSL3spq1jRKNZWbX8jHvx+ge7N6tAvR9yrOG9UtDB9PN2asqsYDwVbzc6OsdKK3ss8+M0o19d3mw6SczmGcbs1fpK6vF7fFNOG7zYc5lpFtdjjmeO89o2h20Yneyj76yCjVUGGh4sPfEmjV0J+ezfUTkJe6t0cEBYWKT1YfMDsUc3zxhVE0u+hEb2XLlhmlGlqx5zh7UrJ4oFczRFxn9ihHaVqvJgPbNWLe2kNkZueZHY5mcTrRW5mnp1Gqof+uTKBxgA+D2zcyOxTLGtczgsycfOb/mWR2KJrF2ZXoRWSAiOwWkX0i8kwxnz8gIttEZLOI/C4ibYp89qxtvd0icr0jg3d51fQx702HTvHngZOMvTYCT3fdFilJ+5DadI2oy8w/DpBXUGh2OJqFlXoW2eZ8fRcYCLQB7iiayG3mKaXaKaWiganAm7Z122DMMdsWGAC8d34OWc0O1TTRz1iVQEANT0Z0alJ65WpuXM9mHM3IZtGWI2aHolmYPc2lzsA+pVSCUioXmA/cWLSCUup0kbe+wPmxVG8E5iulcpRSB4B9tu1p9lixwijVyMG0M/y04xgjuzbF11s/5lGa2JaBtGzgz4xVCdVrCONqeG5UhD2JPhgoehEw2bbsIiLysIjsx2jRP1rGde8XkTgRiUtNTbU3ds0FfbrmIO4ijOoWZnYoTkFEuK9nBLuOZfL7vhNmh6NZlMMugCql3lVKNQP+AUws47ozlFIxSqmYwMBAR4Xk/D780CjVxJmcfL6IS2Jgu0Y0qOVjdjhO44aoRtT38+aTPxLNDqXqvP66UTS72JPoDwNFL5aG2JaVZD4wrJzrakVVs0HNvtl0mMzsfMZ0DzU7FKfi7eHOyK5NWb7rOAmpWWaHUzW+/94oml3sSfTrgeYiEi4iXhg3VxcWrSAiRSfwHAycH3FpITBCRLxFJBxoDvxZ8bCriZ9/Nko1oJTi0zWJRAbX4uqmdcwOx+nc1SUUL3c3Zq9ONDsUzYJKTfRKqXxgPLAE2Al8oZTaISKTRGSordp4EdkhIpuBJ4DRtnV3AF8A8cBPwMNKqQLHH4bm7NYkpLEnJYtR3cL0A1LlEOjvzZCoRny5IZmMc/oBKu1idnVrUEotBhZfsuyFIl8/doV1pwBTyhtgtXZ+LI+HHjI3jiowe3UidWp6MjSqsdmhOK17rgnn642H+TIuiXuvjTA7HM1C9NMoVrZokVFcXPKpsyyLT2FE56b4eOrHLMorMjiAzmF1mbU6kYJCF+9qWaOGUTS76ERvZT/+aBQX99naQ4AxqYZWMff0CCP51DmWxaeYHUrlqibnhqPoRK+ZKjuvgPnrD3Fdm4YE1y5fC+3ZZ59l+vTpjg0MCAsL42fbzfDOnTuzY8cOh+/D0frbvo+f/FFNR7XUiqUTvZX95z9GcWELNx8h/Wweo7uHlVgnLCyMGjVq4OfnR4MGDRgzZgxZWUY3wtTUVD799FPGjRvH3Llz8fPzu6yICJMmTapQnE8++SQvvPBC6RXL6ciRIzRq1IixY8cSGhqKv78/0dHR/FjGVqu7mzC6eyjrDpxkx5GMSorWAiZPNopmF53oreyXX4ziopRSzFqdSMsG/nSNqHvFuosWLSIrK4uNGzcSFxfHyy+/DMCsWbMYNGgQNWrU4K677iIrK+uiMn36dBo0aMB9991nV0wpKcVf8hg6dCi//vorx44dK9tB2mnx4sUMGDCAJk2asHLlSjIyMnj55Ze57bbbSExMLNO2bo9pSk0vd9d+gMrFzw1H04neyhYuNIqL2nDwFPFHTzOqe6jdXSqDg4MZOHAg27dvB+DHH3+kV69exdbdtGkTjz/+OPPnz6dRo5KHO05PT+f999+nc+fOjBkzptg6Pj4+dOzYkSVLltgVZ1H5+flMnjyZsLAw6tWrx7x585g6dSpTpvzVGW3x4sXccMMN/POf/yQsLAw3NzeGDBlCeHg4GzZsKNP+Amp6cvPVISzcfITUzJwyx6u5Hp3oNdPMWp1ILR8Pbupw2fBHJUpKSmLx4sV06NABgG3bttGyZcvL6qWnp3PLLbfw/PPPExsbe9nnhYWFLF26lDvuuIPQ0FCWLl3KhAkTWHiFP6ytW7dmy5YtF94PGTKE2rVrF1uGDBlyod7EiRNZuXIlW7ZsYdasWUyePJnZs2fz6KPGkFB5eXmsWrWK/v37X7S/lJQU9uzZQ9u2be3+/pw35powcgsKmf/noTKvq7kePTyglZ0fy+PJJ82NoxKknM7mp+3HGNM9jJpepf8aDhs2DA8PDwICAhg8eDDPPfccYCR0f3//i+oqpRg1ahSRkZE8/fTTl23rnXfeYerUqdSvX58xY8bw9ttvU79+6dMV+vv7c/To0Qvvv7fjEfzTp08zffp04uPjCQgIoEuXLuzatYspU6ZciHvVqlVERUVddBx5eXncddddjB49mlatWpW6n0s1C/Tj2ub1mbvuEA/GNsNDj+tfrelEb2Vr1pgdQaWZu/YgBUrZPUrlt99+S79+/S5bXqdOHTIzMy9a9tprr7Fjxw42bNhQ7CWhAwcOcOrUKfr160dUVBT16tWzK4bMzExq165tV93zli9fTosWLYiIMB5gys3NJSAggEceeeRCncWLFzNo0KAL7wsLC7n77rvx8vLinXfeKdP+ihrVLYz7Po1jWXwKA9u52Exddv7MNIP+M29lCxYYxcXk5Bcw789D9GkZRNN6NSu0rfbt27Nnz54L71esWMGUKVP46quvSkzKb7zxBvv37ycyMpJHHnmE8PBwnn/+efbu3Vts/fN27txJVFTUhfcDBw4stpePn58fAwcOBIzeNI0b//W074wZMwgODr6o9V400SulGDt2LCkpKSxYsADPCkwl2adVEMG1azB7TWK5t2FZLnpuVBad6LUqt3jbUU5k5V6xS6W9Bg0axMqVKwE4evQoI0aMYPr06Reu4ZckKCiIJ554gq1bt7JgwQLS09Pp1q0b99xzT7H1s7Oz2bBhw0XX0X/88cfLevmcL+e7RYaEhLB582aOHj3KunXrmDNnDsePHyc3Nxcw/rvIycmhdevWADz44IPs3LmTRYsWUaOCT366uwl3dwtlbcJJdh/LLH0FzWXpRG9lr75qFBcze/VBIgJ96XFV6dfFSzNq1CgWL17MuXPn+PDDD0lJSeGxxx67rIX9wAMPlLiNjh078vbbb3PkyJES6y1atIjY2NiLWuf2GDBgANdddx2tW7fmjjvu4OuvvyY6Opo+ffoA8MMPP1xozR88eJAPPviAzZs307Bhwwuxz507t0z7LOq2mCZ4ebgxZ21iubdhSc8+axTNLmK16cdiYmJUXFyc2WFYw4gRxuv8+ebG4UBbktK58d0/+OcNbRhzTbhDtvncc88RFBTE448/7pDtFadLly58/PHHREZGOnS7gwYNYvz48Rddo3e0J7/cwuJtR1n7XF9q+ZT/UpClnO9JpacTvEBENiilYor7TN+MtTIXSvDnzV6diK+XOzd3DHHYNl955RWHbask69atq5TtxsbG0rt370rZ9nmju4Xx1YZkFmxI5m8O+uOqORd96UarMieycvh+61Fu6RiCv6u0LCvo6aefrvC1+NK0CwmgQ9PazFlzkEJXH9VSK5ZO9FbmYuN5/G/dIXILChnlgJuwWtmM6hZKwokz/LFfTyBeHdmV6EVkgIjsFpF9IvJMMZ8/ISLxIrJVRH4RkdAinxWIyGZbcd3n+SvD7t1GcQF5BYV8tu4g1zavT7NAP7PDqXYGtWtEPV8vZq8+aHYojhESYhTNLqVeoxcRd+BdoD+QDKwXkYVKqfgi1TYBMUqpsyLyIDAVuN322TmlVLRjw64mPvvM7AgcZsmOY6SczuGVm9qZHUq15O3hzh2dm/Luin0knTxLk7oVe37BdC50blQFe1r0nYF9SqkEpVQuMB+4sWgFpdSvSqmztrdrAf2nVrvI7NWJNK1bk9iWQWaHUm3d2aUpAny2zkVa9Zrd7En0wUBSkffJtmUlGQsUHUTbR0TiRGStiAwrbgURud9WJy41NdWOkKqJF14wipPbcSSD9YmnuLtrKO5ueuJvszSuXYPr2jTki/VJZOcVmB1OxTz+uFE0uzj0ZqyIjARigGlFFofa+nbeCUwXkWaXrqeUmqGUilFKxQQGBjoyJOeWlGQUJ/fp6oPU8HTntpgmZodS7Y3qHsqps3ks2nLE7FAqZvNmo2h2sacf/WGg6BkaYlt2ERHpB0wAeimlLgyCrZQ6bHtNEJEVQAdgfwVirj4++cTsCCrs1Jlcvt18mOFXhxBQU3epNFu3iHo0D/Lj0zUHuaVjiN3zAGjOzZ4W/XqguYiEi4gXMAK4qPeMiHQAPgCGKqWOF1leR0S8bV/XB64Bit7E1Vzc53FJ5OQXMrq7nvjbCkSEUd1C2XY4g81J6WaHo1WRUhO9UiofGA8sAXYCXyildojIJBEZaqs2DfADvrykG2VrIE5EtgC/Aq9e0ltHuxInH8+joFAxZ81BukbUpVXDWmaHo9ncdHUIft4efLpG35StLuwaAkEptRhYfMmyF4p8fflA4cby1YDuT1deaWlmR1AhP+9M4XD6OZ4f0trudZ599lkaNGhQqePWXKpz58588skn5ZrJyRn5eXtwS8cQ5q07xHODWhPo7212SGXXooXZETgV/WSslc2YYRQnNXt1Io0DfOjXusGFZWFhYdSoUQM/Pz8aNGjAmDFjyMrKAiA1NZVPP/2UcePGMXfu3GLHeRcRJk2a5NA4n3zySV4oZ++md955h5iYGLy9vUucb7ayHDlyhJCQEHJychg7diyhoaH4+/sTHR19YZjkkozsGkpuQSGfr3fSqQad/NyoajrRa5Vib0omq/encVfX0MumsVu0aBFZWVls3LiRuLg4Xn75ZQBmzZrFoEGDqFGjBnfddddlY7xPnz6dBg0acN9999kVQ0pKil31hg4dyq+//sqxY8fKdpBA48aNmThxYonj2FemxYsXM2DAAPLz82nSpAkrV64kIyODl19+mdtuu43ExMQS170qyI8eVxlTDeYXFFZd0JopdKK3siefdNr5YmevScTLw407OjctsU5wcDADBw5k+/btgDGRR69evYqtu2nTJh5//HHmz59Po0YlT4uXnp7O+++/T+fOnS9qYYsI+/btu/B+zJgxTJw4EQAfHx86duzIkiVLynKIAAwfPpxhw4bZPR2hvfLz85k8eTJhYWHUq1ePefPmMXXqVKZMmXKhzvmZqXx9ffnnP/9JWFgYbm5uDBkyhPDwcDZs2HDFfYzqFsrRjGyWxtv3B9FS7r/fKJpd9DDFVnbunNkRlMvp7Dy+3niYoVGNqevrVWK9pKQkFi9ezPDhwwHYtm0bLVu2vKxeeno6t9xyC88//zyx58chL6KwsJCff/6ZTz75hMWLF9OnTx8mTJhQpjHeW7duzZYtWy68HzJkCL///nuxdXv06GHXxOAVMXHiROLi4tiyZQurVq3i6aefxs3NjbVr1wLG5OGrVq1i9uzZl62bkpLCnj17Sr3n0Ld1A0Lq1GD26kQGOducskWmj9RKpxO9lb37rtkRlMuXccmczS1gTAmjVA4bNgwPDw8CAgIYPHgwzz33HGAk9KJzqYIxh+qoUaOIjIzk6aefvmxb77zzDlOnTqV+/fqMGTOGt99+m/r1yz5zlb+/P0ePHr3wvrIT+ZWcPn2a6dOnEx8fT0BAAF26dGHXrl1MmTLlwvdn1apVREVFXfb9ysvL46677mL06NG0atXqivtxdxPu7hrKv37cxc6jp2ndSPeMclX60o3mUIWFijlrEukYWofI4IBi63z77bekp6dz8OBB3nvvvQvjsdepU4fMzIvnNn3ttdfYsWMHs2fPLvbhngMHDnDq1Cmio6OJiooq9yWUzMzMEicTd4TY2FhEpNjSo0ePi+ouX76cFi1aEBERAUBubi4BAQE88sgjF+oUnVD8vMLCQu6++268vLx455137IrrtpgmeHu48akrTiCuXaATvZU54XgeK/ekkph2tlwTf7dv3549Rf4lX7FiBVOmTOGrr74qMQm/8cYb7N+/n8jISB555BHCw8N5/vnn2bt370X1atasydmzZy+8v/TG686dO4mKirrwfuDAgcX2+vHz82PgwIFlPrYVK1aglCq2XHqJ6MiRIxfNTTtjxgyCg4Mvar1fmuiVUowdO5aUlBQWLFiAp6d9TyHX8fViWHQw32w6TMbZvDIfl+YcdKLXHGrW6kSC/L0Z0LZhmdcdNGgQK1euBODo0aOMGDGC6dOn06FDhyuuFxQUxBNPPMHWrVtZsGAB6enpdOvW7aKeMNHR0cybN4+CggJ++umnC/sByM7OZsOGDfTv3//Csh9//PGyXj/nS9Gui/n5+WRnZ1NQUEBBQQHZ2dnk5+eX+diLCgkJYfPmzRw9epR169YxZ84cjh8/Tm5uLmD8F5OTk0Pr1n89n/Dggw+yc+dOFi1aVOYZq0Z3DyM7r5Av4pxoXKXoaKNo9implWFW6dixo9KcU0Jqlgr9x/fq38t2l1gnNDRULVu2rNjPUlNTVXBwsDp79qx66aWXFKB8fX0vK+PGjSs1lpycHLVu3boL79evX6/atGmj/Pz81MiRI9WIESPUhAkTlFJKffHFF+qmm24q49EaXnzxRQVcVF588cVybato7HfffbcKCAhQ4eHhauPGjapfv37qmmuuUUop9fbbb6uHH374Qv3ExEQFKG9v74u+T5999pnd+7z1/dWqx2u/qPyCwgrFrpkHiFMl5FUxPreOmJgYFRcXZ3YYWjm8tGgHn609yB/P9CHI36dc23juuecICgqq0idju3Tpwscff0xkZGSV7bMiBg0axPjx48vUq6g03289wvh5m/h4dAx9izzgpjkPEdmgjJGCL6N73VjZww8br07Q+yYrJ5+v4pIZ1K5RuZM8wCuvvOLAqOyzbt26Kt9nRcTGxtK7d2+HbvP6tg1pWMuHWasTnSPRjxxpvOqZpuyir9FbWY0aRnECX29MJjMnv8QulZrjPP3002W+Dl8aT3c37urSlN/2nmB/apZDt10pkpONotlFJ3ore/11o1hcYaFi1upEoprUpkPTOmaHo5XTiM5N8XJ3Y44e1dLl6ESvVdhv+06QkHqGv+nWvFML9PdmcPtGfLUhmaycivUc0qxFJ3orc5LxPGb9cYBAf2/ne4xeu8zo7mFk5eSzYIO+LOJKdKK3snr1jGJhB06c4dfdqdzVpSleHvrXydlFN6lNVJPazF6TiNV65F2kWzejaHax68wUkQEisltE9onIM8V8/oSIxIvIVhH5RURCi3w2WkT22spoRwbv8v71L6NY2OzViXi6C3d2KXmUSs25jO4WSkLqGX7fd8LsUErmBOeGlZSa6EXEHXgXGAi0Ae4QkTaXVNsExCil2gNfAVNt69YFXgS6AJ2BF0VE361zEZnZeXy1IZkh7RtXqEulZi2D2zeivp8Xs1cnmh2K5iD2tOg7A/uUUglKqVxgPnBj0QpKqV+VUucHElkLhNi+vh5YppQ6qZQ6BSwDBjgm9Grgb38zikUtsN20010qXYu3hzt3dG7KL7uOk3TybOkrmOHmm42i2cWeRB8MFB0EI9m2rCRjgfODgdi1rojcLyJxIhKXmppqR0jVRJMmRrGgwkLF7DUH6dDUuKaruZY7uzTFTYQ5ay3a1TItzennVK5KDr17JiIjgRhgWlnWU0rNUErFKKViAgMDHRmSc5s0ySgWtHJvKgdOnNGteRfVKKAGA9o25PP1SZzLLTA7HK2C7En0h4GizcoQ27KLiEg/YAIwVCmVU5Z1Necz2zZK5cBI3aXSVY3uHkbGuTy+2aRPWWdnT6JfDzQXkXAR8QJGAAuLVhCRDsAHGEn+eJGPlgDXiUgd203Y62zLNHuMHPnXmB4WkpCaxYrdqYzsGqq7VLqwTmF1iAyuxcw/DlBYaOGullqpSj1LlVL5wHiMBL0T+EIptUNEJonIUFu1aYAf8KWIbBaRhbZ1TwKTMf5YrAcm2ZZp9mjZ0igW8+mag3i5X3nib835iQhje4Sz73gWK/da7N5Z375G0eyihynWyiTjXB7d//UL17dtyJu3R5sdjlbJcvMLuXbqclo08GfO2C5mh6NdwZWGKdb/d2tlMv/PQ5zJLWDsteFmh6JVAS8PN0Z3D+O3vSfYdey02eFo5aQTvZWNGGEUi8grKGTW6kS6N6tH28bFT/ytuZ47Ozelhqc7M38/YHYofxk40CiaXXSitzKLzYu5eNtRjmZkc9+1EWaHolWh2jW9uKVjCN9uOkJqZk7pK1SFc+eMotlFJ3ore+YZo1iAUooPf0ugWaAvvVroZx2qm79dE0ZeYaF1H6DSrkgnes0uaxNOsv3wae69NgI3NzE7HK2KRQT60bdVEHPXHiQ7Tz9A5Wx0orcyC43n8dFvCdTz9eKmDlca/UJzZWN7RJB2Jpdv9QNUTkdPDm5lFhlve39qFr/sOs5jfZvj4+ludjiaSbpG1KVt41p8/PsBbu/UBBET/7MbMsS8fTshneit7MknzY4AgI9/P4CXhxt3dwstvbLmss4/QPXEF1tYtfeEufdqLHJuOAt96Ua7orSsHBZsSGZ4h2Dq+3mbHY5mMmPuAW8++i3B7FC0MtCJ3sqGDjWKieauO0ROfiH36gekNIwHqMZcYzxAtf1whnmBxMYaRbOLTvRWZvJ4HudyC5i9OpHeLQO5KsjftDg0axnZNRR/bw/+u3K/2aFodtLX6K3sscdM3f2XG5JIO5PLg7FXmRqHZi21fDy5q2soM1bt52DaGULr+ZodklYK3aLXipVXUMgHKxPoGFqHTmF6ml/tYvdcE4aHuxszVulr9c5AJ3orM3E8j++3HuFw+jke7NXM3G50miUF1fLhlo4hfLkhmeOZ2WaHo5VCJ3oru+EGo1SxwkLF+yv207KBP31aBVX5/jXncP+1EeQXFPLJH4lVv/PbbjOKZhd9jd7KHnrIlN0u33WcPSlZ/Pv2KD3cgVaisPq+DGzXiM/WHOTB2GbU8vGsup2bdG44K7ta9CIyQER2i8g+EblslC0R6SkiG0UkX0RuueSzAtusUxdmntKsSynFeyv2EVKnBje0b2x2OJrFPdirGZk5+cxde6hqd3z2rFE0u5Sa6EXEHXgXGAi0Ae4QkTaXVDsEjAHmFbOJc0qpaFsxt1O4s+nXzyhV6M8DJ9l4KJ37e0bg4a6v7GlXFhkcwLXN6/Px7weqdrCzQYOMotnFnjO5M7BPKZWglMoF5gM3Fq2glEpUSm0FCishxurr9tuNUoXeW7Gfer5e3NqxSZXuV3NeD8Y240RWDl9tSDY7FK0E9iT6YCCpyPtk2zJ7+YhInIisFZFhxVUQkfttdeJSUy02CbGZ7rvPKFVkW3IGK/ekck+PcGp46cHLNPt0i6jH1U1r8/6K/eTm67aeFVXF/+ahtglr7wSmi0izSysopWYopWKUUjGBgXpSC7P855c9BNTw1IOXaWUiIjzWrwWH08/pVr1F2ZPoDwNF/48PsS2zi1LqsO01AVgBdChDfNVbFY7nsf1wBj/vPM7YHuFV23tCcwk9m9cnuklt3v11n27VW5A9iX490FxEwkXECxgB2NV7RkTqiIi37ev6wDVAfHmDrXbGjDFKFZj+815q+Xgw5pqwKtmf5lqMVn1zDqef4+uNVdCqr8JzwxWU2o9eKZUvIuOBJYA7MFMptUNEJgFxSqmFItIJ+AaoA9wgIi8ppdoCrYEPRKQQ44/Kq0opnejtVUW/yEZrPoW/92uhW/NaucW2CCQqJIB3ft3HzR1D8KzMXls6yZeJXQ9MKaUWA4svWfZCka/XY1zSuXS91UC7CsZYfeXlGa+elZt83/pFt+a1ijvfqr9nVhzfbDzMbZ0qsefWiRPGa/36lbcPF6I7SltZ//5GqUQ7jmSwND6Fe3qEE1BDt+a1iundMoj2tlZ9XkElXqu/5RajaHbRid7K7r3XKJXorV/24u/jwd+u0ROLaBUnIjzWtzmHTp7Vk4hbiB7rxspGjqzUzW9LzmDJjhQe69tct+Y1h+nTKoh2wQG8tXwvN0YH4+Wh25Nm0z8BK6vk8TymLtlFnZqejNXTBGoOJCI8eX1Lkk6eY/76Kh4DRyuWTvRWVonjeazef4Lf9p7godirdE8bzeF6Nq9P14i6vPXLPs7k5JsdTrWnL91Y2YMPVspmlVJM/Wk3jQJ89FOwWqUQEZ4e0Irh763mkz8OML5Pc8fuoJLODVelE72VVdKAZsviU9iclM6rw9vh46nHtNEqx9VN69C/TQM+WJnAXV1CqePr5biNV/Fgf85OX7qxsowMozhQQaFi2pLdRNT35ZaOlz36oGkO9dT1LcnKzef9lfsdu+GkJKNodtGJ3spuvNEoDvTtpsPsPZ7F/13XUo83r1W6Fg38ualDMLNXJ3I045zjNnz33UbR7KLPdCt79FGjOMi53ALeWLqbdsEBDIxs6LDtatqV/L1fCwqVYvqyvWaHUm3pRG9lw4cbxUE+/j2BIxnZTBjcWs8Fq1WZJnVrMqpbGF9sSGLHEcdeitTsoxO9lZ048deYHhV0/HQ2763Yz/VtG9A1op5Dtqlp9nq0T3Nq1/Dk5e93opQyO5xqRyd6K3PgeB5vLN1DXkEhzw5s7ZDtaVpZBNT05PF+LViTkMay+BSzw6l2dPdKK/u//3PIZuKPnOaLDUmMvSacsPq+DtmmppXVnV2aMmftQV5ZvJPYlkEVGxrBQedGdaFb9FZ2ww1GqQClFC//EE/tGp484uiHVjStDDzd3ZgwuDWJaWf5dE1ixTbmgHOjOtGJ3sqOHTNKBSyNT2H1/jQe79eCgJp6qAPNXL1bBtGzRSD/+WUvJ8/kln9Du3cbRbOLXYleRAaIyG4R2ScizxTzeU8R2Sgi+SJyyyWfjRaRvbYy2lGBVwsjRhilnM7lFjBpUTwtGvhxZ5emDgxM08rv+cGtOZtbwNSfdpV/I+PGGUWzS6nX6EXEHXgX6A8kA+tFZOElUwIeAsYAT16ybl3gRSAGUMAG27qnHBO+i3vmsr+pZfLOr3s5nH6Oz+/vWrnTumlaGTRv4M/YHuHMWJXArTFN6Bhax+yQXJ49Z39nYJ9SKkEplQvMBy56XFMplaiU2gpcOqXM9cAypdRJW3JfBgxwQNzVw4ABRimH/alZzFiVwPAOwXTR3Sk1i3msb3Ma1vJh4rfbya/Mmag0wL5EHwwUHVQi2bbMHnatKyL3i0iciMSlpqbauelqoJzjeSilePG7Hfh4uPPMoFaVEJimVYyvtwcv3tCGnUdP8+mag2aH4/Is8f+8UmqGUipGKRUTGBhodjjWUc7xPH7YdpTf953g/65rQZC/TyUEpmkVNyCyIb1aBPLmsj2knM42OxyXZk+iPwwUnc49xLbMHhVZV5s40ShlkHE2j5cWxdO2cS1GdtVjzWvWJSK8NLQtuQWFTP4+vvQViirHuVGd2ZPo1wPNRSRcRLyAEcBCO7e/BLhOROqISB3gOtsyzR79+hmlDF7+IZ6TZ3J57eb2enRKzfLC6vsyvvdVfL/1KD+X5YnZcpwb1VmpmUAplQ+Mx0jQO4EvlFI7RGSSiAwFEJFOIpIM3Ap8ICI7bOueBCZj/LFYD0yyLdPskZBgFDv9tjeVLzckM65nBJHBAZUYmKY5zgO9mtGqoT8Tvt1Gxrk8+1bavNkoml3EagMMxcTEqLi4OLPDsIbYWON1xYpSq57Jyee6f6/C29ONxY9eq2eO0pzKtuQMhr33B8M7BDPt1qjSVyjDuVFdiMgGpVRMcZ/psW6s7KWX7K46bclujmSc48tx3XSS15xOu5AAxvWM4L0V+xncvhGxLYPMDsml6Iu4Vtarl1FKsXr/CWavSWRU11BiwupWQWCa5niP9m3OVUF+PPv1NjKz7byEo9lFJ3ors2M8j4yzefzfF1sIr+fLPwbqPvOa8/LxdGfaLe1JOZ3Niwt3mB2OS9GJ3srsGM/j+e+2k5qZw/QR0dT00lfiNOfWoWkdHunTnK83HmbhliNmh+MydGawsldeueLH324yToYnr2tB+5DaVROTplWyR/pcxW97U5nwzTaublqbkDo1L69UyrmhXUy36K2se3ejFCP51Fme/3Y7MaF1eDD2qioOTNMqj4e7G/8Z0QGl4O+fb6agsJiegVc4N7TL6URvZdu3G+USOfkFPDxvEwD/vj0adz3Rt+ZimtStyeRhbVmfeIp3lu+7vMLq1UbR7KIv3VjZ+PHG6yV9hV/5YSdbktL578iraVK3mH9rNc0F3NQhhFV7TjD9lz1cHVqba5sXGQfrueeMV92P3i66RW9l06YZpYiFW44we81B7u0RzoDIRiYFpmlVY8pNkbQI8ufR/23icPo5s8NxWjrRW1mnTkax2Xc8k2cWbKVjaB3dlVKrFmp6efD+yKvJL1A89NkGcvILzA7JKelEb2VFxvPIOJfHuDkbqOHpzrt3Xq1njNKqjYhAP6bdGsWW5AxeWlTGUS41QF+jt7bHHwcg/5fljJ+3kUMnzzJnbBcaBugx5rXqZUBkQx7o1Yz/rtxPq4b+jDI7ICejE72VTZ8OwMs/7OS3vSd47eZ2dNXTAmrV1FPXt2Tf8Uz+uXAHbR5/QQ/3UQb6/38ri45mTnYdZq1O5L5rw7m9U1OzI9I007i7Cf8Z0YFWDWsxZlMeuxpGmB2S09DDFFvYn5//xCs/7qJu7DV8OCqm2vWXz8vLIzk5mexsPc1ccXx8fAgJCcHT09PsUKrU0YxzTPn727i7CRPefISgWvpSJlx5mGKd6C3qzwMnUbGx1PRyJ3z7evy8q99VtgMHDuDv70+9evUQqV5/5EqjlCItLY3MzEzCw8PNDqfKnenWgx1HT/P8Y+/w+biu1K7pZXZIprtSorfr0o2IDBCR3SKyT0SeKeZzbxH53Pb5OhEJsy0PE5FzIrLZVv5boSOpJnYePc3Y2ev54Na/02Tux9UyyQNkZ2frJF8CEaFevXrV9r8dX28PWjbw58CJM4z5ZD1ncvLNDsnSSk30IuIOvAsMBNoAd4hIm0uqjQVOKaWuAv4NvFbks/1KqWhbecBBcbushNQsRs38E18vDyZPHEHtzlebHZKpdJIvWXX/3gTU8OStOzqw7XAG98+JIztP97EviT0t+s7APqVUglIqF5gP3HhJnRuB2bavvwL6SnX/LSyH/alZjJixlsJCxZyxnQmO36TH89C0KxgQ2ZCpN7fnj31pPDx3o072JbAn0QcDSUXeJ9uWFVvHNpl4BnC+H2C4iGwSkZUicm1xOxCR+0UkTkTiUlNTy3QArmJ/ahZ3zFhLoVL87/6uNG/gb4zncX5MD80Ufn5+5V53wIAB1K5dmyFDhpRY5+DBg/Tt25f27dsTGxtLcnLyhc8OHTrEddddR+vWrWnTpg2JiYnljsWV3dwxhJeHRfLLruPc92kc53J1sr9UZV/8PQo0VUqliUhH4FsRaauUOl20klJqBjADjJuxlRyT5ew6dpq7P/4TpRT/u8+W5AE++MDcwLQKeeqppzh79iwfXOHn+OSTTzJq1ChGjx7N8uXLefbZZ5kzZw4Ao0aNYsKECfTv35+srCzc3HRv6Asu+Z6O7BqKl4cb/1iwlb/N+pOPRneqtve2imPPb85hoEmR9yG2ZcXWEREPIABIU0rlKKXSAJRSG4D9QIuKBu1K/jxwklv/uwY34eIkD9CypVE0Q2wszJplfJ2XZ7z/7DPj/dmzxvvPPzfeZ2QY77/+2nh/4oTxftEi4/2xY5Uebt++ffH3979infj4ePr06QNA7969+e677y4sz8/Pp3///oDxn0XNmnqk0guKOTdui2nC9NujWZ94irs/XsepM7kmBWc99iT69UBzEQkXES9gBLDwkjoLgdG2r28BliullIgE2m7mIiIRQHMgwTGhO78lO44x8uN1BPp7s+DB7hcneYCVK42imS4zM5Po6OhiS3x8+cdfiYqK4mvbH6NvvvmGzMxM0tLS2LNnD7Vr12b48OF06NCBp556ioICfUnigkWL/vqjXcSN0cG8e+fV7DhymuHvrybxxBkTgrMgpVSpBRgE7MFokU+wLZsEDLV97QN8CewD/gQibMtvBnYAm4GNwA2l7atjx47K1RUWFqoPV+1X4c98r4a9+7s6mZVTfMVevYxSTcXHx5sdgvL19a3Q+r/++qsaPHhwiZ8fPnxY3XTTTSo6Olo9+uijKjg4WJ06dUp9+eWXqlatWmr//v0qLy9PDR8+XH300UeXrW+F75EpSjk31h9IU9EvLVEdJi1VcYknqywsMwFxqoS8atdFLKXUYmDxJcteKPJ1NnBrMestABbY/VenGsjOK+C5r7fx9abDDIxsyBu3RZU8qffMmVUbnFaizMxMrr222L4EzJs3j8zMTMbZJnKfNGkSQ4cOtWu7jRs3vtCiz8rKYsGCBdSuXZuQkBCio6OJiDAe8x82bBhr165l7NixDjga1xcTVpdvHrqGMZ/8yR0fruW1m9txU4cQs8Myjb5bUYWOZWQzbk4cW5IzeKJ/C8b3vgq3Kw1rEKHH8rAKf39/NtuGjC5JaZ8X58SJE9StWxc3Nzf+9a9/cc899wDQqVMn0tPTSU1NJTAwkOXLlxMTU+xDj1oJwur78vVD1/DgZxv4++db2HgwnYlDWuPt4W52aFVO38avIj/HpzDwP6vYezyL/47syKN9m185yQP8/LNRNKd07bXXcuutt/LLL78QEhLCkiVLAHjhhRdYuNC4zbVixQpatmxJixYtSElJYcKECQC4u7vz+uuv07dvX9q1a4dSivvuu8+0Y3FWdX29mHtvF8b1jGDO2oPc9sHaajlTlR7rppJl5xXw6o+7mLU6kTaNavHWHR24KsjOvtmxscZrNZ0Xc+fOnbRu3drsMCyt2n6PynFu/LT9GE99uQURmDwskqFRjV3q6eIrjXWjL91Uoi1J6fxjwVZ2HcvknmvC+cfAlmX7t9HWn1rTtEuU49wYENmQ1o38eeKLLTw2fzNLd6Tw8rBI6vi6/oBoOtFXgrO5+by5dA8z/zhAoL83n4zpRO9WQWXfUJMmpdfRtOqonOdGaD1fvhjXjQ9W7effy/aw7sBJXrihDTe0b+RSrftL6UTvQEoplsWnMPmHeJJOnuPOLk15ZmAravmUc7zwn34yXgcMcFyQmuYKzj8Yd/vtZV7V3U14KPYqYlsE8czXW3n0f5uY/+chJt0Yaf9lVSejE72DxB85zeTv41mTkMZVQX7Mv79rxaf9e/VV41Unek272PvvG6/lSPTntWlci28euoZ5fx5i2k+7GPifVYzuFsbDva9yucs5OtFX0IETZ3hn+T6+3pRM7RqeTLqxLXd0boqnuwM6NM2fX/FtaJpWInc34e6uoQyMbMjUn3Yx848DfL4+iXG9IrinR3jJz7g4Gdc4ChPsT83ineX7+G7zYbw83Li3RzjjezcnoKYDp3Vr2NBx29I0rUT1/byZeksU914bwbQlu3l96R4++SOR0d3DGNUt1OlnsNL96MtAKcXve09w7+w4+r25kp+2H2Nsj3B+e7oPEwa3cWyShxLH89CsLy0tjd69e+Pn58f48eNLrLdlyxa6detGu3btuOGGGzh9+q+BXbdu3Uq3bt1o27Yt7dq1q7azSVWlFg38+XBUDAse7E5Uk9q8uWwP3V9dzqRF8SSdPGt2eOWmW/R2OHkml0VbjjBn7UH2Hc+inq8XD8dexZhrwqjv5115O37jDeP1hhsqbx9apfDx8WHy5Mls376d7du3l1jv3nvv5fXXX6dXr17MnDmTadOmMXnyZPLz8xk5ciRz5swhKiqKtLS0ajcJuJk6htZh5phO7Dp2mhkrE5i9JpFPVh/g2uaB3Nm5CX1bN3DM5dkqohN9CbLzCli+6zhfbzzMit3HyS9UtA8J4I1boxjcvhE+nlXwGPVXX1X+PpzJ+YdkirrtNnjoIWOY4kGDLv98zBijnDgBt9xy8Wd2PGyTmJjIkCFDLiTr119/naysLP75z39ecT1fX1969OjBvn37rlhvz5499OzZE4D+/ftz/fXXM3nyZJYuXUr79u2JiooCoF69Ct7YdzVVdG60aliLN2+P5qkBLfl8fRKfr0/igc82Ut/Pm8HtGjIkqjEdm9Yp/Sl3k+lEX8TJM7n8uus4P+9MYdWeVM7kFhDk7809PcK5qUMwrRvVqtqA6tev2v1pdps2bRpz5869bHnPnj1566237N5O27Zt+e677xg2bBhffvklSUnGZG579uxBRLj++utJTU1lxIgRPP300w6L3+lV8bnRKKAGj/czxqdauSeVL+OSmb8+idlrDtIowIeBkY3o0yqITuF1LDmWTrVO9OdyC9h46BTrEtJYvT+NjYdOUaggyN+bodHBDGrXkO7N6uNu1l/r85NmDB9uzv6t5kot8Jo1r/x5/foOHUriqaee4qmnnqrwdmbOnMmjjz7K5MmTGTp0KF5exk2//Px8fv/9d9avX0/NmjXp27cvHTt2pG/fvhXep0s4PwHNmDFVulsPdzf6tm5A39YNyMrJ5+f4FL7feoTP1h5k5h8HqOHpTrdm9YhtGUiX8Ho0D/KzRGu/2iT6wkJFYtoZth85zfbDGWw4eIqtyenkFSjcBCKDAxjfpzn9WgcR2TjAEj8czrcMdaI3jYeHB4WFhRfen78h6qgWfatWrVi6dClgtOJ/+OEHAEJCQujZsyf1bS3XQYMGsXHjRp3ozzMp0Rfl5+3BsA7BDOsQzNncfNYmpLFidyordqeyfNdxAPx9PLi6aR1iQutwdWgdWjeqRV0T+ui7XKIvLFQcTj9HwokzJKRmkZB6ht0pmcQfOU1WTj4AXu5utA2uxdgeEXSJqEtMaB38y/v0amWyTSunmadBgwYcP36ctLQ0/Pz8+P777xkwYIDDWvTHjx8nKCiIwsJCXn75ZR544AEArr/+eqZOncrZs2fx8vJi5cqV/P3vf6/w/rTKUdPLgz6tGtCnVQMADqWdZX3iSeIOnmLDwZO8sSz1Qt0Gtbxp1bAWrRvVomVDP8Lq+RJe37dSu3DalehFZADwH8Ad+Egp9eoln3sDnwIdgTTgdqVUou2zZ4GxQAHwqFJqicOiLyLldDajZ/7JgRNnyMn/qwXm7+1BsyA/buoQTGRwLSKDA2jRwN857pgHBJgdQbXn6enJCy+8QOfOnQkODqZVq1Z2rxsWFsbp06fJzc3l22+/ZenSpbRp04Z7772XBx54gJiYGP73v//x7rvvAjB8+HD+9re/AVCnTh2eeOIJOnXqhIgwaNAgBg8eXCnHqDle03o1aVqvJjd3NCY7yTibx7bDGew8epqdx06z82gmq/cnkFfw1+jBdWp60qN5IG/f0cHh8ZQ6TLFtztc9QH8gGWMO2TuUUvFF6jwEtFdKPSAiI4CblFK3i0gb4H9AZ6Ax8DPQQilV4uSX5R2mODe/kIfmbiC8vi8RgX62V18C/bydd7CiCozn4Qqq7RC8ZVBtv0cuMIR3bn4hh06e5cCJMySeOMOBtDPUqenJU9fb35goqqLDFHcG9imlEmwbmw/cCBSdEflG4J+2r78C3hEju94IzFdK5QAHRGSfbXtrynMgV+Ll4cZHozs5erPmcsB4HpqmWZOXhxtXBflVyUBq9iT6YCCpyPtkoEtJdZRS+SKSAdSzLV97ybrBl+5ARO4H7gdo2rSpvbG7vsWLS6+jadWRPjfKxBIXqpVSM5RSMUqpmMDAQLPDsY6aNY1SjVltBjQrqdbfG31ulIk9if4wUHSU/xDbsmLriIgHEIBxU9aedbWSfPaZUaopHx8f0tLSqndCK4FSirS0NHx8fMwOxRzvvWcUzS72XLpZDzQXkXCMJD0CuPOSOguB0RjX3m8BliullIgsBOaJyJsYN2ObA386KniX99FHxuvIkebGYZKQkBCSk5NJTU0tvXI15OPjQ0hIiNlhmOOLL4zXhx4yNw4nUWqit11zHw8sweheOVMptUNEJgFxSqmFwMfAHNvN1pMYfwyw1fsC48ZtPvDwlXrcaJdYtszsCEzl6elJeHi42WFomtOzqx+9UmoxsPiSZS8U+TobuLWEdacAUyoQY/WlRyvUNM0BLHEzVivBrFl/PeqtaZpWTjrRW5lO9JqmOUCpT8ZWNRFJBQ5WYBP1gRMOCsdMrnIcoI/FqlzlWFzlOKBixxKqlCq2f7rlEn1FiUhcSY8BOxNXOQ7Qx2JVrnIsrnIcUHnHoi/daJqmuTid6DVN01ycKyb6GWYH4CCuchygj8WqXOVYXOU4oJKOxeWu0WuapmkXc8UWvaZpmlaETvSapmkuzuUSvYhMFpGtIrJZRJaKSGOzYyovEZkmIrtsx/ONiNQ2O6byEpFbRWSHiBSKiNN1hRORASKyW0T2icgzZsdTESIyU0SOi8h2s2OpCBFpIiK/iki87XfrMbNjKi8R8RGRP0Vki+1YXnLo9l3tGr2I1FJKnbZ9/SjQRin1gMlhlYuIXIcxEmi+iLwGoJT6h8lhlYuItAYKgQ+AJ5VSZZ8v0iT2TKfpTESkJ5AFfKqUijQ7nvISkUZAI6XURhHxBzYAw5zx52Kbkc9XKZUlIp7A78BjSqm1paxqF5dr0Z9P8ja+gNP+JVNKLVVK5dversUYz98pKaV2KqV2mx1HOV2YTlMplQucn07TKSmlVmGMMuvUlFJHlVIbbV9nAjspZgY7Z6AMWba3nrbisNzlcokeQESmiEgScBfwQmn1ncQ9wI9mB1FNFTedplMmFFclImFAB2CdyaGUm4i4i8hm4DiwTCnlsGNxykQvIj+LyPZiyo0ASqkJSqkmwFxgvLnRXllpx2KrMwFjPP+55kVaOnuORdMcTUT8gAXA45f8R+9UlFIFSqlojP/cO4uIwy6r2TUevdUopfrZWXUuxjj6L1ZiOBVS2rGIyBhgCNBXWfyGShl+Ls5GT4lpUbbr2QuAuUqpr82OxxGUUuki8iswAHDIDXOnbNFfiYg0L/L2RmCXWbFUlIgMAJ4GhiqlzpodTzV2YTpNEfHCmEFtockxVXu2G5gfAzuVUm+aHU9FiEjg+V51IlID48a/w3KXK/a6WQC0xOjhcRB4QCnllK0v29SM3hgTrQOsdeIeRDcBbwOBQDqwWSl1valBlYGIDAKm89d0mk47a5qI/A+IxRgSNwV4USn1salBlYOI9AB+A7ZhnO8Az9lmxHMqItIemI3x++UGfKGUmuSw7btaotc0TdMu5nKXbjRN07SL6USvaZrm4nSi1zRNc3E60Wuaprk4neg1TdNcnE70mqZpLk4nek3TNBf3/3ch3wbbTGTWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.stats\n",
    "\n",
    "dom = np.linspace(-3,3,10**3)\n",
    "\n",
    "normal = scipy.stats.norm(0,1)\n",
    "\n",
    "fig,ax = plt.subplots()\n",
    "ax.plot( dom,normal.pdf(dom) )\n",
    "\n",
    "alpha = 0.05\n",
    "l = normal.ppf(alpha/2)\n",
    "u = normal.ppf(1-alpha/2)\n",
    "\n",
    "ax.axvline(l,color=\"red\",linestyle=\":\",label=\"l=-1.96\")\n",
    "ax.axvline(u,color=\"red\",linestyle=\"--\",label=\"u=1.96\")\n",
    "\n",
    "ax.text( 0.50,0.70,r\"P(Z<l)=$\\alpha/2$  \",transform=ax.transAxes,ha=\"center\",va=\"top\",fontsize=12  )\n",
    "ax.text( 0.50,0.60,r\"P(Z<u)=$1-\\alpha/2$  \",transform=ax.transAxes,ha=\"center\",va=\"top\",fontsize=12  )\n",
    "\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means that for a normal distribution $u = -l$, and so \n",
    "\n",
    "\\begin{align}\n",
    "    P(X < \\mu + l\\sigma) &= \\alpha/2\\\\\n",
    "    P(X < \\mu + u\\sigma) &= 1-\\alpha/2\\\\\n",
    "\\end{align}\n",
    "\n",
    "becomes\n",
    "\\begin{align}\n",
    "    P(X < \\mu - u\\sigma) &= \\alpha/2\\\\\n",
    "    P(X < \\mu + u\\sigma) &= 1-\\alpha/2\\\\\n",
    "\\end{align}\n",
    "\n",
    "or\n",
    "\n",
    "\\begin{align}\n",
    "    P(\\mu - u\\sigma < X < \\mu + u\\sigma ) &= 1 - \\alpha/2 - \\alpha/2 = 1 - \\alpha\\\\\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Application to $\\hat \\beta$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we are given the following data set:\n",
    "\n",
    "|x|y|\n",
    "|-|-|\n",
    "|1|0|\n",
    "|2|3|\n",
    "|3|6|\n",
    "|4|7|\n",
    "|5|2|\n",
    "\n",
    "\n",
    "We can compute $\\hat \\beta_{1}$ as\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat \\beta_{1} &= \\frac{\\sum_{i=1}^{N}(x_{i}-\\bar x)y_{i}}{\\sum_{i=1}^{N}(x_{i}-\\bar x)^{2}}\\\\\n",
    "        &= \\frac{\\sum_{i=1}^{N}(x_{i}-3)y_{i}}{\\sum_{i=1}^{N}(x_{i}-3)^{2}}\\\\\n",
    "    &=0.8\n",
    "\\end{align}\n",
    "\n",
    "and we can compute the variance of $\\hat \\beta_{1}$ by first finding$\\hat{\\sigma}^{2}$ and then computing the MLE of the variance: \n",
    "\\begin{align}\n",
    "    Var(\\hat \\beta_{1}) = \\bar{\\sigma}^{2} \\frac{1}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} }\n",
    "\\end{align}\n",
    "\n",
    "From the above data you can find $\\hat{\\sigma}^{2}$ equals $6.7$ and so \n",
    "\n",
    "\\begin{align}\n",
    "    Var(\\hat \\beta_{1}) &= 5.36 \\frac{1}{\\sum_{i=1}^{N} (x_{i} - 3)^{2} } \\\\ \n",
    "                        &= 0.536 \\\\ \n",
    "\\end{align}\n",
    "\n",
    "or \n",
    "\n",
    "\\begin{align}\n",
    "    sd(\\hat \\beta_{1}) &=  \\sqrt{0.536} = 0.732 \\\\ \n",
    "\\end{align}\n",
    "\n",
    "\n",
    "Then our distribution for $\\hat \\beta_{1}$ is \n",
    "\n",
    "\\begin{align}\n",
    "    \\hat \\beta_{1} \\sim \\mathcal{N}(0.8,0.732^{2})\n",
    "\\end{align}\n",
    "\n",
    "We can compute a $95\\%$ confidence interval (or 1-0.05 confidence interval) using this information and our table of quantile values above. \n",
    "\n",
    "$CI(\\hat \\beta_{1})_{1-0.05} = (0.8-1.96*0.732, 0.8+1.96*0.732)$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hypothesis Testing\n",
    "\n",
    "A **hypothesis** is a statement about a parameter, not about a random variable.\n",
    "Even more, we often divide our hypothesis into two parts: a null hypotheis ($H_{0}$) and an alternative hypotehsis $(H_{1})$. \n",
    "Our goal is often to collect data to prove our null hypothesis false, concluding our alternative hypothesis is more likely instead. \n",
    "\n",
    "Our simple linear regression model states \n",
    "\n",
    "\\begin{align}\n",
    "    Y_{i} \\sim \\mathcal{N}(\\beta_{0} + \\beta_{1}x_{i},\\sigma^{2})\n",
    "\\end{align}\n",
    "\n",
    "Two hypotheses of interest come to mind: a hypothesis about $\\beta_{0}$ and a separate hypothesis about $\\beta_{1}$. \n",
    "For $\\beta_{0}$ we may be interested to know whether or not this true parameter equals 0. \n",
    "Our null hypothesis could then be \n",
    "\n",
    "\\begin{align}\n",
    "    H_{0}: \\beta_{0} = 0\n",
    "\\end{align}\n",
    "\n",
    "and so our alternative hypothesis would be \n",
    "\n",
    "\\begin{align}\n",
    "    H_{1}: \\beta_{0} \\neq 0\n",
    "\\end{align}\n",
    "\n",
    "Now that we have written down our null and alternative hypotheses, the next question we may ask is \"how do we determine if the data we collected support our null or alternative hypothesis?\". \n",
    "Our startegy will to build a distribution over a test statistic assuming our null hypothesis is true and compare the distribution of test statistics assuming is the null is true to a test statistic estimated from our data, called the observed test statistic.   \n",
    "\n",
    "From maximum likelihood theory we know that \n",
    "\n",
    "\\begin{align}\n",
    "\\hat \\beta_{0} \\sim \\mathcal{N}\\left(\\beta_{0},  V(\\beta_{0})\\right)\n",
    "\\end{align}\n",
    "\n",
    "where the variance of $\\beta_{0}$ is \n",
    "\n",
    "\\begin{align}\n",
    "    V(\\beta_{0}) = \\sigma^{2}\\left( \\frac{1}{N}  + \\frac{\\bar{x}^{2}}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} }  \\right)\n",
    "\\end{align}\n",
    "\n",
    "Because $\\hat \\beta_{0}$ has a Normal distribution we can rewrite the above distribution for $\\hat \\beta_{0}$\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat \\beta_{0} &\\sim \\mathcal{N}\\left(\\beta_{0},  V(\\beta_{0})\\right)\\\\\n",
    "    \\hat \\beta_{0} - \\beta_{0} &\\sim \\mathcal{N}\\left(0,  V(\\beta_{0})\\right) & \\text{Shift}\\\\\n",
    "    \\frac{ \\hat \\beta_{0} - \\beta_{0}}{\\sqrt{V(\\beta_{0})}} &\\sim \\mathcal{N}\\left(0,1 \\right) & \\text{Scale}\\\\\n",
    "\\end{align}\n",
    "\n",
    "#### Note:\n",
    "The quantity $\\sqrt{V(\\beta_{0})}$ is often called the standard deviation of $\\beta_{0}$ or $sd(\\beta_{0})$ and the maximum likelihood estimate of $sd(\\beta_{0})$ is often called the standard error of $\\beta_{0}$ or $se(\\beta_{0})$.\n",
    "The above says that, even though we do not know the true parameter value $\\beta_{0}$ we do know that as we collect more and more data points $\\frac{ \\hat \\beta_{0} - \\beta_{0}}{\\sqrt{V(\\beta_{0})}}$ will converge to something we do know---a normal distribution with mean 0 and variance 1 (called a standard normal distribution).\n",
    "\n",
    "Because $\\frac{ \\hat \\beta_{0} - \\beta_{0}}{\\sqrt{V(\\beta_{0})}}$ converges to a distribution we know let us use this quantity as our test statistic. \n",
    "If our Null hypothesis was true than $\\beta_{0}$ would equal 0 and so\n",
    "\n",
    "\\begin{align}\n",
    "    \\frac{ \\hat \\beta_{0} - \\beta_{0}}{\\sqrt{V(\\beta_{0})}} &\\sim \\mathcal{N}\\left(0,1 \\right)\\\\\n",
    "    \\frac{ \\hat \\beta_{0} - 0}{\\sqrt{V(\\beta_{0})}} &\\sim \\mathcal{N}\\left(0,1 \\right) & \\text{Assumed our Null hypothesis}\\\\\n",
    "    \\frac{ \\hat \\beta_{0}}{\\sqrt{V(\\beta_{0})}} &\\sim \\mathcal{N}\\left(0,1 \\right)\n",
    "\\end{align}\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above says that **if** $\\beta_{0}$ was truly 0 then $\\frac{ \\hat \\beta_{0}}{\\sqrt{V(\\beta_{0})}}$ should follow a standard normal distribution.\n",
    "Even more, if we estimate  $\\frac{ \\hat \\beta_{0}}{\\sqrt{V(\\beta_{0})}}$ using our MLE estimate for $\\beta_{0}$ than if  the true $\\beta_{0}$ was zero we should expect this fraction to be close to 0. \n",
    "\n",
    "Lets approximate. \n",
    "\n",
    "\\begin{align}\n",
    "    \\sqrt{V(\\beta_{0})} &= \\sqrt{\\sigma^{2}\\left(\\frac{1}{N}  + \\frac{\\bar{x}^{2}}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} }\\right)} & \\text{Exact}\\\\\n",
    "    \\sqrt{ \\hat V(\\beta_{0})} &= \\sqrt{\\hat \\sigma^{2}\\left(\\frac{1}{N}  + \\frac{\\bar{x}^{2}}{\\sum_{i=1}^{N} (x_{i} - \\bar x)^{2} }\\right)} & \\text{Estimated}\n",
    "\\end{align}\n",
    "\n",
    "The expected value of $\\hat \\beta_{0}$ is $\\beta_{0}$, $E(\\hat \\beta_{0}) = \\beta_{0}$. \n",
    "However, we know our MLE for $\\beta_{0}$ will get closer to the true $\\beta_{0}$ as our data increase.\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat \\beta_{0} = \\bar Y - \\hat \\beta_{1} \\bar x\n",
    "\\end{align}\n",
    "\n",
    "If the null hypothesis was true than\n",
    "\n",
    "\\begin{align}\n",
    "     \\frac{ \\hat \\beta_{0}}{\\sqrt{V(\\beta_{0})}} &\\sim \\mathcal{N}\\left(0,1 \\right)\n",
    "\\end{align}\n",
    "\n",
    "and we can compare to this null distribution the following observed test statistic\n",
    "\n",
    "\\begin{align}\n",
    "    \\frac{\\bar Y - \\hat \\beta_{1} \\bar x}{ \\sqrt{ \\hat V(\\beta_{0})} }\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "For example, suppose we observed the following data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.9/site-packages/seaborn/_decorators.py:36: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x1c0026bb0>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEHCAYAAABFroqmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhGUlEQVR4nO3df5RcdZnn8ffTnYTOb/KjDTGQNC1ZAUlEphdhNI6TjC5mcRhUZHAHWQdt5ywaHJ0dEWfPOMdxRvyBS1ZFm0VFVydkRcR1M1EE5xD3ANpBCD8CkxATSE4InQRI0kknndSzf3RVT6erqrt+3Hu/91Z9Xufk0F1V3fWlknOf+32e7/f5mrsjIiIyUkvoAYiISPooOIiISBEFBxERKaLgICIiRRQcRESkiIKDiIgUmRDqjc2sDXgAOCU/jh+6+9+a2ZnAGmAOsBG42t2PjfW75s6d6x0dHTGPWCRZ2/r6Aehsnxp4JNKoNm7cuNfd20s9Fyw4AEeB5e5+yMwmAr8ys38GPg58xd3XmNk3gGuBW8f6RR0dHfT29sY/YpEEXfnNBwG488MXBx6JNCoz21HuuWBpJR9yKP/txPwfB5YDP8w/fgfwJ8mPTkSkuQWtOZhZq5k9CrwI3As8C7zs7sfzL9kJLAg0PBGRphU0OLj7CXc/HzgduBA4u9KfNbNuM+s1s96+vr64higi0pRSsVrJ3V8GfglcDJxqZoVayOnArjI/0+PuXe7e1d5esp4iIiI1ChYczKzdzE7Nfz0ZeBuwmaEg8Z78y64B7gkyQBGRJhZytdJ84A4za2UoSK1195+a2VPAGjP7e+C3wO0BxygiUrNcztm+r589BwaYN6ONjjlTaWmx0MOqSLDg4O6bgDeUeHwbQ/UHEZHMyuWc9U++wMfXPsrAYI62iS3c/N7zueR1p2UiQKSi5iAi0mi27+sfDgwAA4M5Pr72Ubbv6w88ssooOIiIxGDPgYHhwFAwMJjjxYMDgUZUHQUHEZEYzJvRRtvEky+xbRNbeNX0tkAjqo6Cg4hIDDrmTOXm954/HCAKNYeOOdnolRVytZKISMNqaTEued1pnL1qGS8eHOBV07VaSUREGAoQne3T6GyfFnooVVNwEBGJUJb3Noyk4CAiEpE49jaECjYqSIuIRCTqvQ2FYLNy9Qauuu1hVq7ewPonXyCX8yiHXZKCg4hIRKLe2xByI53SSiIiESnsbRgZIKrd2zAyjXRk8ETZYBN3kVszBxGRiNS7t2F0Gumx518OtpFOMwcRkYjUu7dhdBppbe9Orl+xmFvu23JSgTuJjXQKDiIiEapnb8PomsXuVwb47oM7uOMDF+J4ohvpFBxERFKiVM3ipcPHaJ9+SuIb6VRzEBFJiTT1Y9LMQUQkRV47bzpf/08XMPWUCcybfgoLZ4fZYR0sOJjZGcB3gXmAAz3ufouZfQb4ENCXf+mN7r4uzChFRJJRbnf1wtmlZw1x75wOOXM4DnzC3R8xs+nARjO7N//cV9z9SwHHJiKSqHIb3s5etayo3pDEEaTBag7uvtvdH8l/fRDYDCwINR4RkZCq2V2dxM7pVBSkzawDeAPwcP6hj5jZJjP7lpnNKvMz3WbWa2a9fX19pV4iIpIZ1Zwcl8QRpMGDg5lNA+4CPubuB4BbgdcA5wO7gS+X+jl373H3Lnfvam9vT2q4IhJQLuds6zvEg8/uZVvfoUQa0FWjnvFVs1IpiSNIg65WMrOJDAWG77v7jwDcfc+I528DfhpoeCKSIknk2UOOr5rd1YVAMvq9olzyGnK1kgG3A5vd/eYRj8939935by8HnggxPhGJVr2ra6op2IYQxfgq3V2dxBGkIWcObwKuBh43s0fzj90IXGVm5zO0vHU78OEQgxOR6ERx1z9Wnj0NwSHp8cV9BGmw4ODuvwJK/avQngaRBhPFXXUU7bDjlPbxVSt4QVpEGl8Uq2vS1FqilLSPr1pqnyEisYvirjqJPHs90j6+amnmICKxi+quupBnv6hzLp3t01J34U37+KqhmYOIxK7R7qqbgYKDiCQi7tU1Ei0FBxHJtLi7kzYrBQcRyay075rOMhWkRSSzkuhO2qwUHEQks5LoTtqslFYSkcyKe1dyM9czNHMQkcyKc1dyoZ6xcvUGrrrtYVau3sD6J19IXZvwuGjmICKZFef+ibR3gY2bgoOIVCStKZa49k+kvQts3BQcRBpQ1BfyEEtGQwejRuuyWi3VHEQaTBy58qSXjKYh399oXVarpZmDSIOJI1eedIolDfn+Zu8HpZmDSIOJY+1/Egfaj5SW/QtRd1nN5ZxtfYd48Nm9bOs7lOqVT8GCg5mdYWa/NLOnzOxJM7s+//hsM7vXzLbk/zsr1BhFsiiOC3nSKZakg1ES0pAqq0bImcNx4BPufi5wEXCdmZ0L3ADc5+6Lgfvy34tIheK4kBdSLOtWLWNN9xtZt2pZrMXouINRiDv4rLX6CHmG9G5gd/7rg2a2GVgAXAa8Nf+yO4B/AT4ZYIgimRRXrjzJlttx5vtDNevL2tLYVBSkzawDeAPwMDAvHzgAXgDmlfmZbqAbYOHChQmMUiQ7GuHshLj+H0IVu+tZGhtiWW/wgrSZTQPuAj7m7gdGPufuDpSc77l7j7t3uXtXe3t7AiMVkbFkpdgaqthda6osVK0i6MzBzCYyFBi+7+4/yj+8x8zmu/tuM5sPvBhuhCJSiSydqxBqc1utqbJQM52Qq5UMuB3Y7O43j3jqJ8A1+a+vAe5JemwiUp0sFVtDbm6rZWlsqJlOyJnDm4CrgcfN7NH8YzcCnwfWmtm1wA7gvWGGJyKVylKxNWub20LNdEKuVvoVUO5vY0WSYxGR+mStD1GWCvaFmc7olF3cM51UrFYSkWwLdQFrBqFmOgoOIlK3ai5gobutZlGImY6Cg4hEopILWJZWNTW74PscRKR5ZGlVU7NTcBCRxKSl26qMT8FBRBLTiN1WG5WCg4jEolQ7jWY/XS1LVJAWkbqNXoG0cNYUfr55T8nCc5Y2oDUzBQeRJhbFstJSK5B6ru4asx9QJcsyteQ1LAUHkSYV1bLSUiuQenfsr6udhpa8hqeag0iTimpZaakVSDmnqPC8aM5kJk9srailt5a8hqfgIJJxtZ6jENWy0lIrkP7PY7u46d1Lhx9fNGcyH12+mCt7HqroTAIteQ1PaSWRDKs1/ZLLOVMmTYikWV6pvkqfvOQc3n7OPJYsmMmLBweYPLGVK3seqvhMgqw18mtEmjmIZFgt6ZdCQFm15hFWLV9c97LSQl+ldauWsab7jaxbtYxLXncaEya0DJ9dcPjYiapmAlryGp5mDiIZVss5CiMDyvce2sG1b+6ktQVWnP0qliw4taaC73h9laqdCWTtzIVGpJmDSIbVsuN4ZEDZ/coAX/vlVlbft5Ujgydiu/jWMhOo5dQ0iU7oM6S/BVwKvOju5+Uf+wzwIaAv/7Ib3X1dmBGKpNvofP+iOZP57GVL2HNgYPj50RfVEPl8zQSyJ3Ra6TvAV4Hvjnr8K+7+peSHI5ItIy+6+/uPsuvlAbq/1ztmcTqpg3lKbWLLyulrEjg4uPsDZtYRcgwiWVdIvwD82e2/HndFUBJ38drEln2hZw7lfMTM3g/0Ap9w95dCD0gk7aopTsd9stj2ff3ctH4z1765E8vHgpvWb+bs06Zr5pARaQwOtwKfBTz/3y8Dfz76RWbWDXQDLFy4MMnxiaRSnLWEavsc7es/ypVdC1l9/5bhmcOq5YvZ339UwSEjUrdayd33uPsJd88BtwEXlnldj7t3uXtXe3t7soMUSaG49gYUUkQrV2+oaHczwKTWluHAAEMzmNX3b2Fia/Elp9Yd3mONN8rf16xSN3Mws/nuvjv/7eXAEyHHI5IVcdUSym20K7e7GSi76e3wsRMnPRZ1bUK1jugEnTmY2T8BDwKvNbOdZnYt8AUze9zMNgF/CPxlyDGKZEmhlnBhxxwAHv7dvrrvnmvpc1Ru/8W8GSenuKJusKeGfdEJvVrpqhIP3574QEQaSNR3z7XUMipdLlvLDu+xRP37mlnq0koiUr2RBeMpk1qrTgONpZZ9EZWmuOopopcqkqthX3QUHEQybvRMYdWKsyK9e661llHJctlaN+SVmx29/Zx5iWzwawYKDiIZNzrPXjhop96753I7nKN27vzp3PGBCzl87DgLZ0/lzLnjB57n9vfz9AsH+OCyTgDu2riTj699dLgjrNp01E/BQSTjRufZ79q4k1XLF5+0x6Dau+ckVv2Ue48z544/a3jkuZfpeWDbSXsovvfQjuHZkdp01E/BQSTj5s9sY9WKsygsSLpr407u7H2OO7sv4sjgiZrunmtZvlqtWt9j+75+brz78aI9FN1v6VRtIUIKDiIZlcs5v9vbz+bdBzCGgsJLh49x/YrFLJ43reazGSCZVT+1vke5n/t386arthCh1O2QFmkEce/SzeWc+5/Zwz2P7eLpPQcB+Iu3dDJryiRuuW8LZ86p7/yDWs6JSOo9yv3cOafNUG0hQgoOIhGrpd1EtZ7b38+WPYfoeWAbX71/K998YBuHB0/w/osXMTCYo+9Q6Q1qlQatqFpxjPV+tb5HuZ8br1Yh1VFaSSRiSeTr9xw4yi33ndy76Jb7tvDF97y+7N13NUXmKFpxjPd+9SyR1Yqk+GnmIBKxWtpNjJbLOQODJzgwMFjyDr//2PGS73Hk2PGyd9/Vtpao95jOSt6v1vfQEaLxU3AQiVi9+frCHfemXa+weffBkmmpRbOnlnyPs0+bXna5aRRBqxpJv19U1NV1iIKDSMTqzdcX7rg9f00qdcd95tzS73HeGCuUSgWtRXMmM3liaywXwiSK2lFLol6UFeae/f/prq4u7+3tDT0MkWGF3cW15MQffHYvV932cNHja7rfyEWdc2t+j9E1gEVzJvPR5Yv5mx8/UVQTAKo63KeS98tC++xtfYdYuXpD0e7ydRHWi9LEzDa6e1ep51SQFolBPcdwVto8rtr3GF3InTyxlSt7HiqqCbz2o8t4Zs/Bui/qWSwcq6vrv1FaSSRlCmmpwtnLUTaPG1nILXcgz3P7SxeSf7N9f9Wpp6wVjrOYCouLgoNIAqopchbuuJcumMk586cPN5OL+sJa7kI4ZdKEkkFjw9a9DZ+Dj+uo1SxSWkkkZrXk3ltajLaJrbRNbI0tnVGuXfa8GaeUTGu5x7NnI02ymAqLS9DgYGbfAi4FXnT38/KPzQbuBDqA7cB73f2lUGMUqVcSm+LKKdV2u3ChK3chBIqCRqHraWH8jZyDr6de1EhCzxy+A3wV+O6Ix24A7nP3z5vZDfnvPxlgbCJFxrrYlnvtv+45yAeXdXLXxp3sfmVojX8SF9hKZizlLoSFoLFjXz+/ff5lvvfQjuGxN2sOvtmEPkP6ATPrGPXwZcBb81/fAfwLCg6SAtWkh0q9tnD3vfuVgUQusCNnLPNntvGuC07n6RcOsODUySxZMHPMVEkhaHTMmcqRwRwvHT4GNHcOvtmEnjmUMs/dd+e/fgGYF3IwIgXVpIdKvXb1/Vu49s2d3P6rbYlcYAvLMufPbOPqixYNH/7T88C2ipemKgffvMquVjKzdSXu6hPlQzv0Si6LMLNuM+s1s96+vr6ERybNaLx2ECNXJPUdPMqsKZOKXrt0wYzYVh+NNm9GG4vmTOZTK89h4PgJPrisk/kz28btqTRa1pajSjTGmjl8G/i5md0BfMHdBxMa0x4zm+/uu81sPvBiqRe5ew/QA0M7pBMamzSxsTanlUojXb9iMd998ORc/eJ50xMrdC6cNYWPLl/MX//wsZKprUYuKkv9ys4c3P1/AxcAM4BeM/srM/t44U+MY/oJcE3+62uAe2J8L5GKjbUGvlQa6Zb7tnBF1+lFr03Kcy8dHm6NURjT6vu38K4LTldRWcY1Xs3hGNAPnAJMB3Jjv7w6ZvZPDBWf55rZTuBvgc8Da83sWmAH8N4o31OkVmPl38ulnN5wxqms6X5jkFx9uTG1tqCisoyrbHAws0uAmxm6k7/A3Q9H/ebuflWZp1ZE/V4iUSi39LNcymnRnKmxpG4qWVJbbkwrzn5VXedLS3MYq33Gp4Er3P2GOAKDSCMZr+1ClGcEVNpWutyYFBikEmVnDu6+LMmBiNQrl3Oe29/PngNH6T92nEWzp3Lm3GRSOWOlnKJuXV3pklotQ5V6pHGfg0jVcjnn/mf2sGXPoeGzlZM+P6Bcyinq9hm1tJVugGNbJGHqyioNYfu+fjbtfGU4MMD4ZyQnJcrjMnM5Z8qkCRW1ldapZlIPBQdpCHsODJDLdw0dKQ1nFkd1RkDhYr9qzSOsWr543LbS5WYsoYOlZIPSStIQ5s1oo9Wo6AS1pJVrjV3NUtJcznl818s8/cIB3vn6Bax/YjfXvrmT1hbKrj4qNWOZNWUSfQeP1nX8pzQHBQdpCB1zprLk9Jlcv2JxUc0h9Hr+egvD4zXx+/3XzCn5u0YvZZ0/s433X7yIa77968yc6SzhKDhIQ2hpMZa/dh5ntU/jgoWzOHzsOAsTXK1UyfhqPSNgvCZ+pWZGuZzTYvAPly/hxrsfZ2AwxxVdp5esyTTqwT1SHwUHaRgtLUbH3Gl0zG2sC121O50LM42b1m/mit87g7/749cxf2YbBqyucpWTNC8FB5E61HL4T7X5/mp3Om/f189N6zdzZdfC4TbdbRNb+MfLl7JozmR27Dty0u8JXZORdNJqJZEaVbNUtJ5lpdXudN5zYIBLly4YDgwwNEP41N2b+OxlS8Zd5SQCmjmI1Kzew38qzfdXW9CeN6ON1pbSy3onthrrtGNaKqDgIFKFkamhFjNmTZk0fF4DlM/h17KreaRqCtodc6by7xfNLpmKmjejrebCuDQXpZVEKjQ6NXTNt3/N+y9exPyZ/5azL5fDj2ojXCVaWoyLO+dw07uXKoUkNdPMQaRC5Q706X5LJ6vv2zrmBTiKjXDVmDChhXcufTVLFsxUCklqouAgUqF6DvQJ0SG1nr0VIgoOIhWq90AfXawlS1JbczCz7Wb2uJk9ama9occjMt6BPiKNJO0zhz90972hByHVqXWzV2jjjVuH50gzSXtwkIyp5dSzNASTSsedRGqo8HkcGBhkUmvLUJ8kBSBJWGrTSoADPzezjWbWPfpJM+s2s14z6+3r6wswPCml2jME0nIgTVrOPhj5eWzefZBNu17RAT0SRJqDw5vd/QLgHcB1ZvaWkU+6e4+7d7l7V3t7e5gRSpFqTz1Ly0U5ytPa6jH683BHB/RIEKkNDu6+K//fF4G7gQvDjkgqUe1mr3IX5T0HBtjWd4gHn93Ltr5Dsd85J7lJbSxpCVIiqQwOZjbVzKYXvgbeDjwRdlRSiWpX9JS7KA+e8ERTTQtnTUnFjuK0BCkRc09fLtPMOhmaLcBQ0fwH7v65cq/v6ury3l6tdk2LQkG1khU9pQrBN717KTff+0xRa+l1MR1KM/L8g0uXLqC1BboWzeb3O+cwYUKy90+jPw8z+Nr7LtBpbRILM9vo7l0ln0tjcKiWgkO2jQ4m+/qPcsU3Hip63ZruN3JR59zI339b3yFWrt5QtLktrmA0nsLncd0PHmFSawt3/5c3KTBILMYKDlrKKsGVWh5aaidyXKmVejumRq3wecxomzj8fSlpWAIsjSuVNQdpbqXqFje9eyn7+o/GUpzOYp6/kH76wHd+za+27uPHj+7i/23dy/HjufF/WKQCmjlI6ozcibznwACDJ5z/ds/j7Nh3pKJNddVKumNqFModBXrTu5fyzqWv1gxC6qbgIKlUSK0AJ9UDqjlBrZr3ylpbjHJHgX7yrk0sWTBTzf2kbkorSaolte6/EIwu6pxLZ/u0VAcGGPsoUO2JkCho5iBjCl30LNcmO831gCSMdRRos382Eg3NHKSsNPQ9Upvs0nQUqMRNMwcpq1zfoyjz/ePJYj0gKdUcBRp6BijZo+AgZSW9/r/cBUwnqJVXyWdTSxt1EaWVpKwk1/+nIYXVqNLS+VayRcFBykoy368LWHzU6VVqobSSlJVkvj9tLSwaiVZ8SS00c5AxJbX+P4stLLJCK76kFpo5SCpksYVFVmjFl9RCwUFSQReweGnFl1RLwUHGldQaeV3AoqW9DVIPBQcZk9bIZ5P+3qReqS1Im9klZvaMmW01sxtCj6dZaYnp2HI5Z1vfIR58dm8sZ03USn9vUq9UzhzMrBX4GvA2YCfwGzP7ibs/FXZkzUdLTMtL8925/t6kXqkMDsCFwFZ33wZgZmuAy4CSwWFbXz9XfvPBBIfXPAYGT2AGI48aN4N/XPc0bRNbww2sDgODJzh2Isek1pa6/h8GBk+wadcrw5/NwGCO637wCEsXzIzks3lq9wGAmv5tN+LfmyQrrWmlBcDzI77fmX9smJl1m1mvmfUODg4mOrhm0jaxlbPap2H5G2EzOKt9WmYvMPv7j7Fp1yts3n2QTbteYX//sZp/17ETuZMuvjB0MT52IvxRnY329ybJS+vMYVzu3gP0AHR1dfmdH7448IgaV2HVS9aXmG7rO8TK1RuGL+ju8PxLh/nm1b9XU6ql8PtG7zz+2vsuiCR1U5gx1Ppvu1H+3iQ+a/+i/HNpnTnsAs4Y8f3p+cckgKydklZO1D2G0r7zuFH+3iSMtM4cfgMsNrMzGQoKfwq8L+yQJOui7jGkjXvSyFI5c3D348BHgJ8Bm4G17v5k2FFJHJJcChrHnb7uzqVRpXXmgLuvA9aFHofEJ+mloLrTF6lcKmcO0hxCbNTSnb5IZRQcJBgdQiOSXqlNK0njq7dAXGljOTWgE6megoMEU88ZDpXWK9Lc4kIkzRQcJJh6CsTl6hVnr1p20ga0Sl8nIidTcEi5Rk+J1HqGQ6WN5dSATqQ2Cg4plsaUSFqCVaX1iqg3vok0C61WSrG09eQvBKuVqzdw1W0Ps3L1BtY/+ULiZxjkco47fOk9r+f6FWcxf2Zb2XpF2ltciKSVZg4plraUSBry96VmU/9w+RIuWHgqC2cXz2K08U2kNpo5pFghJTLSeCmRONtRpGFfQqkAdePdj5Nzyl7wtfFNpHoKDilWbUok7rRPLcEqamkIUCLNQGmlFKs2JVJv2me8YnM9+xKiogKzSDIUHFKumqWe9dQoKlkZlYb8fRoClEgzUHBISBJLQOu5q6501lHrvoSopCFAiTQDBYcEJLVfoZ676rStjBpL6AAl0gwUHBKQ1BLQeu6qlcsXkZG0WikBSa6wqXXZpjaLichIqZs5mNlngA8BffmHbsyfCpdZWbgrVy5fREZK68zhK+5+fv5PpgMDZOeuXJvFRKQgdTOHRqS7chHJmrTOHD5iZpvM7FtmNqvUC8ys28x6zay3r6+v1EtSRXflIpIlQYKDmf3CzJ4o8ecy4FbgNcD5wG7gy6V+h7v3uHuXu3e1t7cnN/iMibPXkog0riBpJXf/o0peZ2a3AT+NeTgNK43nQYhINqQurWRm80d8eznwRKixZF3azoMQkexIY0H6C2Z2PuDAduDDQUeTYVna9Swi6ZK64ODuV4ceQ9zG67NUTx+mkT87ZdKE1O+vEJF0Sl1waHTj1QHqqROM/tlFcybz939yHn/z4yfUwVREqqLgkLDx+izV04dp9M/u2HeE/3H/Fu7svogjgye0v0JEKpa6gnSjG6/PUj19mEr97I59RzgyeEL7K0SkKgoOCRvvqM16juJMwzGeItIYFBwqEOVGsvH6LNXThykrPZxEJP1UcxhH1BvJxuuzVE8fJvVwEpGoKDiMI46DesY7yayek850SpqIREFppXEkeVCPiEhaKDiMQ0VeEWlGCg7jUJFXRJqRag7jUJFXRJqRgkMFVOQVkWajtJKIiBRRcBARkSIKDiIiUqRpaw71nJkgItLogswczOwKM3vSzHJm1jXquU+Z2VYze8bM/kMc719oibFy9Qauuu1hVq7ewPonX6irZ5KISCMJlVZ6AngX8MDIB83sXOBPgdcBlwBfN7PWqN9cZyuLiIwtSHBw983u/kyJpy4D1rj7UXf/HbAVuDDq91dLDBGRsaWtIL0AeH7E9zvzj0VKLTFERMYWW3Aws1+Y2RMl/lwW0e/vNrNeM+vt6+ur6mfVEkNEZGyxrVZy9z+q4cd2AWeM+P70/GOlfn8P0APQ1dVVVSVZLTFERMaWtqWsPwF+YGY3A68GFgO/juON1BJDRKS8UEtZLzezncDFwP81s58BuPuTwFrgKWA9cJ27nwgxRhGRZhZk5uDudwN3l3nuc8Dnkh2RiIiMlLbVSiIikgIKDiIiUkTBQUREiph79vsJmVkfsCP0OGowF9gbehApo8/kZPo8iukzOVk9n8cid28v9URDBIesMrNed+8a/5XNQ5/JyfR5FNNncrK4Pg+llUREpIiCg4iIFFFwCKsn9ABSSJ/JyfR5FNNncrJYPg/VHEREpIhmDiIiUkTBISXM7BNm5mY2N/RYQjKzL5rZ02a2yczuNrNTQ48pFDO7JH9c7lYzuyH0eEIyszPM7Jdm9lT+iOHrQ48pDcys1cx+a2Y/jfp3KzikgJmdAbwdeC70WFLgXuA8d18K/CvwqcDjCSJ/PO7XgHcA5wJX5Y/RbVbHgU+4+7nARcB1Tf55FFwPbI7jFys4pMNXgL8Gmr4A5O4/d/fj+W8fYuhMj2Z0IbDV3be5+zFgDUPH6DYld9/t7o/kvz7I0AUx8lMis8TMTgf+I/A/4/j9Cg6B5U/G2+Xuj4UeSwr9OfDPoQcRSCJH5maRmXUAbwAeDjyU0P47QzeVuTh+edoO+2lIZvYL4LQST30auJGhlFLTGOvzcPd78q/5NEOphO8nOTZJNzObBtwFfMzdD4QeTyhmdinwortvNLO3xvEeCg4JKHdkqpktAc4EHjMzGEqhPGJmF7r7CwkOMVHjHSFrZv8ZuBRY4c271rriI3ObhZlNZCgwfN/dfxR6PIG9CfhjM1sJtAEzzOx/ufufRfUG2ueQIma2Hehy96ZtKmZmlwA3A3/g7n2hxxOKmU1gqCC/gqGg8BvgffnTEpuODd093QHsd/ePBR5OquRnDn/l7pdG+XtVc5C0+SowHbjXzB41s2+EHlAI+aL8R4CfMVR8XdusgSHvTcDVwPL8v4tH83fNEhPNHEREpIhmDiIiUkTBQUREiig4iIhIEQUHEREpouAgIiJFFBxEIpbvIPo7M5ud/35W/vuOwEMTqZiCg0jE3P154Fbg8/mHPg/0uPv2YIMSqZL2OYjEIN/qYSPwLeBDwPnuPhh2VCKVU28lkRi4+6CZ/VdgPfB2BQbJGqWVROLzDmA3cF7ogYhUS8FBJAZmdj7wNoZOLftLM5sfdkQi1VFwEIlYvoPorQydOfAc8EXgS2FHJVIdBQeR6H0IeM7d781//3XgHDP7g4BjEqmKViuJiEgRzRxERKSIgoOIiBRRcBARkSIKDiIiUkTBQUREiig4iIhIEQUHEREpouAgIiJF/j8zqxzLOu+rPwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np \n",
    "x = np.random.normal(0,2,100)\n",
    "e = np.random.normal(0,4,100)\n",
    "y = 10+ 4*x + e\n",
    "\n",
    "fig,ax = plt.subplots()\n",
    "sns.scatterplot(x,y,ax=ax)\n",
    "ax.set(xlabel=\"X\",ylabel=\"Y\")\n",
    "\n",
    "ax.axhline(0)\n",
    "ax.axvline(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we can compute the above observed test statistic, we would find it equals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MLE for Beta 1 is 3.98\n",
      "MLE for Beta 0 is 9.34\n",
      "MLE for sigma2 is 15.40\n",
      "MLE for V(b0) is 0.15\n",
      "MLE for square root of V(b0) is 0.39\n",
      "Test statistic (W) is 23.75\n"
     ]
    }
   ],
   "source": [
    "avgy = np.mean(y)\n",
    "avgx = np.mean(x)\n",
    "N = len(y)\n",
    "\n",
    "mlebeta1 = sum((y-avgy)*(x-avgx))/sum((x-avgx)**2)\n",
    "print(\"MLE for Beta 1 is {:.2f}\".format(mlebeta1))\n",
    "\n",
    "mleBeta0 = avgy - mlebeta1*avgx\n",
    "print(\"MLE for Beta 0 is {:.2f}\".format(mleBeta0))\n",
    "\n",
    "mlesigma2 = np.mean( (y - (mleBeta0+mlebeta1*x))**2 )\n",
    "print(\"MLE for sigma2 is {:.2f}\".format(mlesigma2))\n",
    "\n",
    "mleV = mlesigma2*((1/N) + (avgx**2)/sum((x-avgx)**2))\n",
    "print(\"MLE for V(b0) is {:.2f}\".format(mleV))\n",
    "\n",
    "print(\"MLE for square root of V(b0) is {:.2f}\".format(mleV**0.5))\n",
    "\n",
    "W = mleBeta0/np.sqrt(mleV)\n",
    "print(\"Test statistic (W) is {:.2f}\".format(W))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can form our statistic W and compare this statistic to a standard normal distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1c42c3130>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsV0lEQVR4nO3de3RV1bn38e9D7uTCXVsJCGhAATFoxAtqVWzFimBra7V1KKc9MuyBoZ5qrX21avXYYuvRvih9qefIsKeXl2q1r7Si1l5otdRKBIoGi3IrBFC5h9xvz/vHXjtswk6yk+wkK8nvM0YGe88111pzhZ0nM8+aa05zd0REpO8a0NMNEBGRrqVALyLSxynQi4j0cQr0IiJ9nAK9iEgfl9rTDWhu+PDhPmbMmJ5uhohIcm3cGPl3woQuOfxbb721191HxNsWukA/ZswYiouLe7oZIiLJddFFkX9XruySw5vZP1vaptSNiEgfp0AvItLHKdCLiPRxocvRx1NXV0dpaSnV1dU93RTpYzIzM8nPzyctLa2nmyLSZRIK9GY2E/jfQArw3+6+sIV6VwO/BM5y9+Kg7JvAV4AG4BZ3f6W9jSwtLSU3N5cxY8ZgZu3dXSQud2ffvn2UlpYyduzYnm6OSJdpM3VjZinAYuByYCJwnZlNjFMvF7gV+FtM2UTgWmASMBP4YXC8dqmurmbYsGEK8pJUZsawYcP0l6L0eYnk6KcBm9x9i7vXAsuAOXHqPQg8DMT+1MwBlrl7jbtvBTYFx2s3BXnpCvpcSX+QSKAfCeyIeV8alDUxszOAUe7+Ynv3FekN3J1fvlXKmu0HeropIu3W6VE3ZjYAeBS4vRPHmGdmxWZWvGfPns42SSTp/rjxI+549u9c/X9WUV5T39PNEWmXRAL9TmBUzPv8oCwqF5gMrDSzbcA5wHIzK0pgXwDc/Ul3L3L3ohEj4j7B2+PMjNtvP/K77JFHHuH+++9vdZ9t27YxefJkAFauXMmsWbOOqbNy5UoGDRrE1KlTmTBhAhdeeCG/+c1vmrYvWbKE//mf/2nxHCtXrmTVqlUtbl++fDkLF0bunc+dO5df/vKXrba5ue985ztHvT/vvPPatX9L5s6dy8iRI6mpqQFg7969xE59sXv37qO+X9/97nc5+eSTmTBhAq+8Ev9+/t13382oUaPIyck5qvyJJ55g6dKlnWrvS29/AIA7vP7+3k4dS6S7JRLoVwMFZjbWzNKJ3FxdHt3o7ofcfbi7j3H3McAbwOxg1M1y4FozyzCzsUAB8GbSr6IbZGRk8Pzzz7N3b/J/yC+44ALWrl3Lxo0bWbRoEQsWLOD3v/89ADfffDM33HBDi/u2Fujr6+uZPXs2d911V4fb1jzQt/ZLpb1SUlJaDMCPPvooN910EwAbNmxg2bJllJSU8PLLL/Nv//ZvNDQ0HLPPlVdeyZtvHvvx+vKXv8zjjz/eqba+s6uMc8cNI2WAUbLrUKeOJdLd2hxe6e71ZrYAeIXI8Mql7l5iZg8Axe6+vJV9S8zsGWADUA/Md/djf0Lb4du/LmHDrrLOHOIYE0/I474rJ7VaJzU1lXnz5vHYY4/x0EMPHbVt7ty5zJo1i8997nMA5OTkUF5e3qG2FBYWcu+99/LEE08wY8YM7r//fnJycrjjjjtYtGgRS5YsITU1lYkTJ7Jw4UKWLFlCSkoKP/3pT3n88cd56qmnyMzMZO3atUyfPp0pU6ZQXFzME088AcDvfvc7Fi5cSFlZGY8++iizZs3i6aefPqrOrFmzuOOOO3j55ZepqqqisLCQSZMm8bOf/azp2tydO++8k5deegkz45577uELX/gCK1eu5P7772f48OG88847nHnmmfz0pz+Ne9Pztttu47HHHmsK6LGee+45/uM//gOAF154gWuvvZaMjAzGjh3LySefzJtvvsm555571D7nnHNO3O/pwIEDGTNmDG+++SbTprV/LEBdQyObPjrMV84fx97yGt7dfbjdxxDpSQmNo3f3FcCKZmX3tlD3ombvHwIeile3t5k/fz5Tpkzhzjvv7NLznHHGGXz/+98/pnzhwoVs3bqVjIwMDh48yODBg7n55pubfhEAPPXUU5SWlrJq1SpSUlJ4+umnjzrGtm3bePPNN9m8eTMXX3wxmzZtarEdCxcu5IknnmDdunXHbHv++edZt24df//739m7dy9nnXUWF154IQBr166lpKSEE044genTp/OXv/yF888//5hjjB49mvPPP5+f/OQnXHnllU3lW7duZciQIWRkZACwc+fOo4J4fn4+O3cekwFsVVFREa+99lqHAv2WPRXUNTinfCyX0gOVrNtxsN3HEOlJveLJ2Fht9by7Ul5eHjfccAOLFi0iKyury87T0oLtU6ZM4Utf+hJXXXUVV111VYv7f/7znyclJf7jCtdccw0DBgygoKCAcePG8Y9//KNDbXz99de57rrrSElJ4fjjj+cTn/gEq1evJi8vj2nTppGfnw9E/kLZtm1b3EAP8M1vfpM5c+ZwxRVXNJXt3r2bZN+rOe644zp8rVv3Rv46O/m4HLbtq+A363dTU99ARmq7HwkR6RGa66adbrvtNp566ikqKiqaylJTU2lsbASgsbGR2traTp1j7dq1nHrqqceUv/jii8yfP581a9Zw1llnUV8ff/RHdnZ2i8dunkIxs6PaD3T6AaJoTxwiefiW2glQUFBAYWEhzzzzTFNZVlbWUW0YOXIkO3YcGaVbWlrKyJHtG6VbXV3d4V/OOw9G2jJycBYnDI4c44NDeshKeg8F+nYaOnQo11xzDU899VRT2ZgxY3jrrbeAyCiXurq6Dh9//fr1PPjgg8yfP/+o8sbGRnbs2MHFF1/Mww8/zKFDhygvLyc3N5fDhxPPGT/77LM0NjayefNmtmzZwoQJExgzZgzr1q1rOkfsDc20tLS413PBBRfwi1/8goaGBvbs2cOf//znDqVFIDJa5pFHHml6P378eLZt29b0fvbs2Sxbtoyamhq2bt3K+++/33SuGTNmJJTGee+995pGQLXX7oNVZKWlMHhgGiODQL/roAK99B4K9B1w++23HzX65qabbuJPf/oTp59+On/9619b7VHH89prrzUNr5w/fz6LFi1ixowZR9VpaGjg+uuv57TTTmPq1KnccsstDB48mCuvvJJf/epXFBYW8tprr7V5rtGjRzNt2jQuv/xylixZQmZmJtOnT2fs2LFMnDiRW265hTPOOKOp/rx585pSRrE+85nPMGXKFE4//XQuueQSvve97/Gxj32sXdcdNWnSpKPOmZ2dzUknndR0/2DSpElcc801TJw4kZkzZ7J48WJSUlJobGxk06ZNDB06FIA777yT/Px8Kisryc/PP2r461/+8hc++clPdqh9uw5V8fHBmZhZU49+18GqDh1LpEe4e6i+zjzzTG9uw4YNx5RJ3/b888/73Xff3Wqdt99+2//93/+9zWOtWbPGr7/++ha3t/X5mv3E6379f7/h7u5VtfV+4jd+44t+916b5xU5yic+EfnqIkRGQcaNq73uZqz0D5/5zGfYt29fq3UmT57Mo48+2uax9u7dy4MPPtjhtnxwqIrxx0VuDmempTBkYBq7y5S6kd5DgV5C61//9V+TcpyOpmwg8hfv/opahuceucE8PCeD/eWdu+Eu0p2UoxdpRVl1PXUNzrDs9Kayodnp7Kuo6cFWibSPAr1IK/ZXRHruw3KOBPrhORnsq1CPXnoPBXqRVuwrj/Tch2YfSd0MzU5v+gUg0hso0CeotLSUOXPmUFBQwEknncStt97aqQej5s+fT2FhIRMnTiQrK4vCwkIKCwsTnl3y4MGD/PCHP2zXOX/wgx9QWVnZ4XovvPDCUU/kRmeUjPr1r3/N7Nmz29WmsIv23Junbg5W1lHX0NjSbiKhokCfAHfns5/9LFdddRXvv/8+7733HuXl5dx9993tOk7sjIuLFy9m3bp1rFixgpNOOol169axbt26ponR2tITgf68887jjTfeaHr/17/+lby8PD766CMgMrNlsqYxDot95fFSN5HXByrVq5feQYE+AX/4wx/IzMzkX/7lX4DIY/2PPfYYS5cupbKykqeffpoFCxY01Z81axYrV64EIjNZ3n777U0PU7WmoqKCL3/5y0ybNo2pU6fywgsvAFBSUsK0adMoLCxkypQpvP/++9x1111s3ryZwsJCvv71rx9znCuuuILTTz+dyZMn84tf/IJFixaxa9cuLr74Yi6++GIAvvrVr1JUVMSkSZO47777AOLWixoxYgR5eXlNDzLt3LmTq6++umnq4lWrVjF9+vSOfItDa39FNHUT26OPpHH2aeSN9BK9b3jlbbdBnNkUO6WwEH7wgxY3l5SUcOaZZx5VlpeXx+jRo1ud/REiQffss8/mP//zP9tsxkMPPcQll1zC0qVLOXjwINOmTePSSy9lyZIl3HrrrXzpS1+itraWhoYGFi5cyDvvvBN3ZsmXX36ZE044gRdfjKzseOjQIQYNGsSjjz7KH//4R4YPH950vqFDh9LQ0MCMGTNYv349t9xyyzH1Yk2fPp1Vq1bR0NBAQUEB55xzDq+88gqzZs3i73//O2eddVab19mb7C2vJTcj9agJzKK9e+XppbdQj76LpaSkcPXVVydU97e//S0LFy6ksLCQiy66iOrqarZv3865557Ld77zHR5++GH++c9/tjk512mnncarr77KN77xDV577TUGDRoUt94zzzzDGWecwdSpUykpKWHDhg1ttvG8885j1apVrFq1inPPPZdp06bxt7/9jbVr13LKKaeQmZmZ0LX2Fvsrahkak7aBI/l6jbyR3qL39ehb6Xl3lYkTJx5zk7SsrIzt27dz8skns379+hZnf8zMzGxxyuDm3J3nnnuOCRMmHFV+6qmncvbZZ/Piiy/y6U9/mh/96EeMGzeuxeOMHz+eNWvWsGLFCu655x5mzJjBvfcevXzA1q1beeSRR1i9ejVDhgxh7ty5Cc1aOX36dB5//HEaGhq46aabyM3Npbq6mpUrV/a5/DxEAn3sjViAYTnR1I3G0kvvkFCP3sxmmtlGM9tkZsesS2dmN5vZ22a2zsxeN7OJQfkYM6sKyteZ2ZJkX0B3mDFjBpWVlU1rtzY0NHD77bczd+7cptWLWpr9sT0uu+wyHn/88ab56NeuXQvAli1bGDduHLfccgtz5sxh/fr1rc5auWvXLgYOHMj111/P17/+ddasWQNw1D5lZWVkZ2czaNAgPvzwQ1566aWm/Vs79qmnnsquXbt4/fXXmTp1KhCZc37JkiV9Lj8PQY++WaAfnJWGGRxQj156iTYDvZmlAIuBy4GJwHXRQB7j5+5+mrsXAt8DYicg2ezuhcHXzUlqd7cyM371q1/x7LPPUlBQwPjx48nMzGxaT7W12R/b41vf+hZ1dXVMmTKFSZMm8a1vfQuIpFgmT55MYWEh77zzDjfccAPDhg1j+vTpTJ48+ZibsW+//XbTzdtvf/vb3HPPPUBkJsqZM2dy8cUXc/rppzN16lROOeUUvvjFLx4VpGPrxftenH322QwbNoy0tDQAzj33XLZs2dIne/SHqurIy0o7qmzAACM3I5VDVR2fjlqkO1m099hiBbNzgfvd/bLg/TcB3P27LdS/DrjB3S83szHAb9w94YnAi4qKvLi4+Kiyd999N+5CHCLJ0Nrn67T7X+HqM/K5f/bRK5td8L0/UHTiUB77QmE3tFD6hIsuivwbjMhLNjN7y92L4m1LJHUzEtgR8740KGt+kvlmtplIj/6WmE1jzWytmf3JzC5ooYHzzKzYzIr37NmTQJNEul5jo1NeU39Mjx4gLzNNPXrpNZI26sbdF7v7ScA3gHuC4t3AaHefCnwN+LmZ5cXZ90l3L3L3omSvFSrSUYdr6nGHvMxjxywMylKgl94jkUC/ExgV8z4/KGvJMuAqAHevcfd9weu3gM3A+I40tK0Uk0hHtPa5KgsC+aA4PfpBWWlN20XCLpFAvxooMLOxZpYOXAssj61gZgUxb68A3g/KRwQ3czGzcUABsKW9jczMzGTfvn0K9pJU7s6+fftaHPsf7bHHS92oRy+9SZvj6N293swWAK8AKcBSdy8xsweILF21HFhgZpcCdcAB4MZg9wuBB8ysDmgEbnb3/e1tZH5+PqWlpSh/L8mWmZlJfn5+3G1l1UGgz4yTo1egl14koQem3H0FsKJZ2b0xr29tYb/ngOc600CAtLQ0xo4d29nDiLRLWVU9AHlZ8XP0NfWNVNc1kJmW2ANxIj1FUyCItCDao4+Xo4+mc6J1RMJMgV6kBWWt5OijI3F0Q1Z6AwV6kRaUVdVhBjnp8VM3AIeC9I5ImCnQi7SgrLqe3IxUBgywY7ZFA7169NIbKNCLtKCsqo5BA49N20Bsj16BXsJPgV6kBYeq6uIOrQTdjJXeRYFepAVl1S0H+pyMSN7+cLVy9BJ+CvQiLSirqo87hh4gI3UAaSlGeY0CvYSfAr1IC8qq6+KOoYfIvPw5GamUq0cvvYACvUgLWsvRA+RkpqpHL72CAr1IHPUNjVTWNpDbWqDPSFOOXnoFBXqROCpqGgDIjTMXfVRuRirlNRp1I+GnQC8Sx+EggEdH18Sj1I30Fgr0InFEe/Q5rfTodTNWegsFepE4oimZbPXopQ9QoBeJozzao89oea753IxU3YyVXiGhQG9mM81so5ltMrO74my/2czeNrN1Zva6mU2M2fbNYL+NZnZZMhsv0lWiKZmcjNZG3aRSU99IbX1jdzVLpEPaDPTBmq+LgcuBicB1sYE88HN3P83dC4HvAY8G+04kssbsJGAm8MPoGrIiYVYRpGSyW+nRR/P3FUrfSMgl0qOfBmxy9y3uXgssA+bEVnD3spi32UB0Fe85wDJ3r3H3rcCm4HgioXY4CN65bfToAeXpJfQSWTN2JLAj5n0pcHbzSmY2H/gakA5cErPvG832HRln33nAPIDRo0cn0m6RLpVIjz46xl55egm7pN2MdffF7n4S8A3gnnbu+6S7F7l70YgRI5LVJJEOq6ipJzNtAKkpLf+IRPP36tFL2CUS6HcCo2Le5wdlLVkGXNXBfUVC4XBNfasPS8GRHL2ejpWwSyTQrwYKzGysmaUTubm6PLaCmRXEvL0CeD94vRy41swyzGwsUAC82flmi3StikQCveakl16izRy9u9eb2QLgFSAFWOruJWb2AFDs7suBBWZ2KVAHHABuDPYtMbNngA1APTDf3Ru66FpEkqa8ur7Vh6XgSI5eqRsJu0RuxuLuK4AVzcrujXl9ayv7PgQ81NEGivSE8nb06DUNgoSdnowViSORQD8wPQUz9egl/BToReKoqGk7dRNdZUo5egk7BXqROMpr6luduTIqMie9Ar2EmwK9SByJpG4gmMFSPXoJOQV6kWbqGxqprmtMLNCrRy+9gAK9SDPRRUfaytED5GSmNc2LIxJWCvQizZTXRic0SzBHX60nYyXcFOhFmonm3BPq0St1I72AAr1IM9HAncioG92Mld5AgV6kmaZA38oUxVE5GalU1DbQ2Oht1hXpKQr0Is1U1LS9jGBUdGRORa169RJeCvQizRzJ0SfQo29aTlBz9Ul4KdCLNFOewDKCUdkZmpNewk+BXqSZ8gSWEYzKbQr06tFLeCnQizRTUVNPRmrrywhGNa0ypZE3EmIK9CLNHK6pb1pUpC3Z6UrdSPglFOjNbKaZbTSzTWZ2V5ztXzOzDWa23sx+b2YnxmxrMLN1wdfy5vuKhE0iUxRHHVllSqkbCa82P81mlgIsBj4JlAKrzWy5u2+IqbYWKHL3SjP7KvA94AvBtip3L0xus0W6TiLrxUY13YzVNAgSYon06KcBm9x9i7vXAsuAObEV3P2P7l4ZvH0DyE9uM0W6z+EE1ouNit6wrahVj17CK5FAPxLYEfO+NChryVeAl2LeZ5pZsZm9YWZXxdvBzOYFdYr37NmTQJNEuk5FbX1CE5oBZKSmkJ4yQKtMSagl9mlOkJldDxQBn4gpPtHdd5rZOOAPZva2u2+O3c/dnwSeBCgqKtKz5NKjyqvrGTc88R+NnMzUpqdpRcIokR79TmBUzPv8oOwoZnYpcDcw291rouXuvjP4dwuwEpjaifaKdLnymoaEJjSLys5I0QyWEmqJBPrVQIGZjTWzdOBa4KjRM2Y2FfgRkSD/UUz5EDPLCF4PB6YDsTdxRUKnvKYu4ZuxEJkTR6kbCbM2P83uXm9mC4BXgBRgqbuXmNkDQLG7Lwe+D+QAz5oZwHZ3nw2cCvzIzBqJ/FJZ2Gy0jkiotGcZwajcDKVuJNwS+jS7+wpgRbOye2NeX9rCfquA0zrTQJHu1J5lBKOyM1LYW17bVU0S6TQ9GSsSoz3LCEblZKapRy+hpkAvEqM9ywhG5WSkaIFwCTUFepEY7Zm5MipHOXoJOQV6kRhNc9G3a3hlKpW1DTRoOUEJKQV6kRgVNR1J3UQnNlOvXsJJgV4kxpGFwdsxvLJpOUEFegknBXqRGNGbse0J9Nnq0UvIKdCLxFDqRvoiBXqRGOXBMoJpCSwjGNUU6DUNgoSUAr1IjPJ2LCMYlaMcvYScAr1IjPJ2LCMYFV03Vg9NSVgp0IvEaM8yglFN68YqdSMhpUAvEqM9ywhGResrdSNhpUAvEqM9ywhGpaUMICN1gEbdSGgp0IvEKO9Ajx4i6RsFegkrBXqRGO1dRjAqO0OBXsIroUBvZjPNbKOZbTKzu+Js/5qZbTCz9Wb2ezM7MWbbjWb2fvB1YzIbL5Js7V1GMEozWEqYtRnozSwFWAxcDkwErjOzic2qrQWK3H0K8Evge8G+Q4H7gLOBacB9ZjYkec0XSZ7oMoLR4ZLtkZ2RqnVjJbQS6dFPAza5+xZ3rwWWAXNiK7j7H929Mnj7BpAfvL4MeNXd97v7AeBVYGZymi6SXNFlBDuSusnNSKWiVoFewimRQD8S2BHzvjQoa8lXgJfas6+ZzTOzYjMr3rNnTwJNEkm+6DKCOe1YdCQqOyNV4+gltJJ6M9bMrgeKgO+3Zz93f9Ldi9y9aMSIEclskkjCKpqmKE5r9745GnUjIZZIoN8JjIp5nx+UHcXMLgXuBma7e0179hUJg8PV7V9GMCpXo24kxBIJ9KuBAjMba2bpwLXA8tgKZjYV+BGRIP9RzKZXgE+Z2ZDgJuyngjKR0KnowDKCUdkZqVTXNVLf0JjsZol0WpufaHevN7MFRAJ0CrDU3UvM7AGg2N2XE0nV5ADPmhnAdnef7e77zexBIr8sAB5w9/1dciUinVTegbnoo3KapkFoYNBAPZ4i4ZLQJ9rdVwArmpXdG/P60lb2XQos7WgDRbpLR5YRjIruc7imjkED25/jF+lK6nqIBDqyjGDUkTnpG5LaJpFkUKAXCXRkGcGoI+vG1iW1TSLJoEAvEujIMoJRTakbjaWXEFKgFwl0ZBnBqFylbiTEFOhFAh1ZRjBKqRsJMwV6kUBFTX2HJjSDI6mbcvXoJYQU6EUCh6vrOzShGUB2euRpWs13I2GkQC8SqKht/8LgUakpA8hKS9EMlhJKCvQigfLqjgd60Jz0El4K9CKB8pqGDt+MBa0bK+GlQC8SqKip79Bc9FHZGSlaTlBCSYFeBKhraKSqroHczI7PU5OjxUckpBToRTjyRGteB0fdQGTBEqVuJIwU6EWAw9WRB50616NPUaCXUFKgF+FIj76jUyBAZAZL5egljBIK9GY208w2mtkmM7srzvYLzWyNmdWb2eeabWsws3XB1/Lm+4qEQVkSevTZGakcVqCXEGqz+2JmKcBi4JNAKbDazJa7+4aYatuBucAdcQ5R5e6FnW+qSNcpqwpy9FmdGF6ZkUptfSO19Y2kp+qPZQmPRD6N04BN7r7F3WuBZcCc2Aruvs3d1wNaMFN6pWiOPq+TPXpA6RsJnUQC/UhgR8z70qAsUZlmVmxmb5jZVe1pnEh3SUqOvmliMwV6CZeOf6oTd6K77zSzccAfzOxtd98cW8HM5gHzAEaPHt0NTRI5WjRH35kpEKK/JBToJWwS6dHvBEbFvM8PyhLi7juDf7cAK4Gpceo86e5F7l40YsSIRA8tkjSHq+vJTk8htQOrS0Vlq0cvIZXIp3o1UGBmY80sHbgWSGj0jJkNMbOM4PVwYDqwofW9RLrf4eq6To24gZjUjZ6OlZBpM9C7ez2wAHgFeBd4xt1LzOwBM5sNYGZnmVkp8HngR2ZWEux+KlBsZn8H/ggsbDZaRyQUyqrqOzXiBiAvK/KLIpoGEgmLhD7Z7r4CWNGs7N6Y16uJpHSa77cKOK2TbRTpcodrOt+jHxQE+kNVCvQSLhrsK0IkR9+ZETdwZGjmoUoFegkXBXoRoKyq8z369NTIKlNK3UjYKNCLEOnRd2bmyqhBWWlK3UjoKNBLv+fuQeqmcz16iEyhEJ1OQSQsFOil36upb6S2obHTOXpQj17CSYFe+r1oTj06PLIz8jLTlKOX0FGgl34vGatLRalHL2GkQC/9XllVdC76zgf6PAV6CSEFeun3DgaBefDA9E4fKy8rsm5sY6N3+lgiyaJAL/3ewcpaAAYnIUc/KCsN9yPpIJEwUKCXfu9g8CTrkGT06IP0j27ISpgo0Eu/Fw30yRh1o/luJIwU6KXfO1RVR15mKikDrNPHylOglxBSoJd+72BlbVJuxMKRHn2ZAr2EiAK99HsHKusYPLDzaRtQj17CSYFe+r2DVXXJ79HrZqyESEKB3sxmmtlGM9tkZnfF2X6hma0xs3oz+1yzbTea2fvB143JarhIshyqrE3K0Eogsu7sAOOA5qSXEGkz0JtZCrAYuByYCFxnZhObVdsOzAV+3mzfocB9wNnANOA+MxvS+WaLJE+kR5+cQG9mDMlObxqbLxIGifTopwGb3H2Lu9cCy4A5sRXcfZu7rwcam+17GfCqu+939wPAq8DMJLRbJCkaG51DVXVJ69EDDB2Yzv4KBXoJj0QC/UhgR8z70qAsEQnta2bzzKzYzIr37NmT4KFFOq+sug53GJSkHD3AkOw0BXoJlVDcjHX3J929yN2LRowY0dPNkX7kyFOxSezRZ6tHL+GSSKDfCYyKeZ8flCWiM/uKdLkjE5olN9DrZqyESSKBfjVQYGZjzSwduBZYnuDxXwE+ZWZDgpuwnwrKREIhetN0UFbyUjdDB6ZzoLKWBs1gKSHRZqB393pgAZEA/S7wjLuXmNkDZjYbwMzOMrNS4PPAj8ysJNh3P/AgkV8Wq4EHgjKRUIimbpLZox+SnY67HpqS8EhopQV3XwGsaFZ2b8zr1UTSMvH2XQos7UQbRbrMgSROURw1NDvy18H+itqm1yI9KRQ3Y0V6yv6KWgZYcqYojooG9wMaSy8hoUAv/dre8hqGZmcwIAkzV0ZFf2nsK1egl3BQoJd+bW95LcNzkpteUY9ewkaBXvq1feU1DOuiQK+x9BIWCvTSr+2rqGVYdkZSj5mZlsLA9BQOKNBLSCjQS7+2r7w26T16iOTp1aOXsFCgl36ruq6B8pp6huckt0cPkfTNPgV6CQkFeum3ooE42TdjAUbkZrDncE3SjyvSEQr00m/tK48E4mTn6AGOz8vgIwV6CQkFeum3ouPcuyJHPyI3k30VNdQ3NF+iQaT7KdBLv/VhWTUQSbMk23G5GbhHxumL9DQFeum3PggC/XG5mUk/9vF5kWN+dLg66ccWaS8Feum3PiyrZnhOBumpyf8xOC74K+GjMuXppecp0Eu/tftQNR8blPy0DcBxeUGg1w1ZCQEFeum3PjhUzcfykp+2ARiek4HZkfsAIj1JgV76rQ/Lqpty6cmWljKAYdnp6tFLKCQU6M1sppltNLNNZnZXnO0ZZvaLYPvfzGxMUD7GzKrMbF3wtSTJ7RfpkOq6Bg5U1nVZjx7gY4My2X2oqsuOL5KoNleYMrMUYDHwSaAUWG1my919Q0y1rwAH3P1kM7sWeBj4QrBts7sXJrfZIp0TvUl6/KCuC/T5gweyaU95lx1fJFGJ9OinAZvcfYu71wLLgDnN6swBfhy8/iUww8ySt5KDSJJFe9of78JAP3JIFjsPVOGuRcKlZyUS6EcCO2LelwZlcesEi4kfAoYF28aa2Voz+5OZXRDvBGY2z8yKzax4z5497boAkY7YcSAS6EcNGdhl58gfkkVVXYNmsZQe19U3Y3cDo919KvA14Odmlte8krs/6e5F7l40YsSILm6SCGzfX8kAgxMGZ3XZOUYGxy49oDy99KxEAv1OYFTM+/ygLG4dM0sFBgH73L3G3fcBuPtbwGZgfGcbLdJZO/ZX8vFBWV3ysFTUyCGRQL/zoAK99KxEPuWrgQIzG2tm6cC1wPJmdZYDNwavPwf8wd3dzEYEN3Mxs3FAAbAlOU0X6bjt+ysZPbTr0jYA+UFaqPRAZZeeR6QtbY66cfd6M1sAvAKkAEvdvcTMHgCK3X058BTwEzPbBOwn8ssA4ELgATOrAxqBm919f1dciEh7bN9fySUTjuvScwzKSiMvM5V/7lOgl57VZqAHcPcVwIpmZffGvK4GPh9nv+eA5zrZRpGkqqytZ8/hGkYP69oePcC4ETls2VPR5ecRaY2ejJV+J9rDHtXFqRuAk4/LYbPG0ksPU6CXfue9Dw8DMOH43C4/10kjcvjocA1l1XVdfi6RlijQS7/z3oeHSR1gjB2e3eXnOmlE5BybP1KvXnqOAr30Oxs/OMy4EdldOrQy6uTjcgDYpEAvPUiBXvqdjR8eZnw3pG0ARg8dSEbqAP7xweFuOZ9IPAr00q+U19SzY39Vt+TnAVJTBjDphDzWlx7slvOJxKNAL/3Kuu0HAZgyanC3nXNK/mDe2VlGfUNjt51TJJYCvfQrb/3zAGYwdfTgbjtn4ajBVNU18L7y9NJDFOilXyn+534mHJ9LXmZat53z9OCvhzXbD3TbOUViKdBLv1HX0Mja7Qc588Qh3XreMcMGMnJwFq+9t7dbzysSpUAv/cabW/dTXlPPJ8Z371TYZsaF40fwl017qVOeXnqAAr30G78t+YDMtAFcUND9ax58YvwIDtfUs3qb5vST7qdAL/1CfUMjL5d8wAUFI8hKT+n28184fjg5Gak891bzpRxEup4CvfQLL5d8wIdlNXyhaFTblbvAwPRUrjz9BF58exeHqjTvjXQvBXrp8xoanSV/2syJwwZyySldOwd9a64/ZzTVdY08+efNPdYG6Z8U6KXP+/Gqbbyzs4yvfXI8AwZYj7Vj0gmDmFN4Av/92tamGTRFukNCgd7MZprZRjPbZGZ3xdmeYWa/CLb/zczGxGz7ZlC+0cwuS2LbRdr0m/W7eGjFu1xyynHMPv2Enm4O/+vTp5KXlcZXfryarXu1IIl0jzYDfbDm62LgcmAicJ2ZTWxW7SvAAXc/GXgMeDjYdyKRZQUnATOBH0bXkBXpCo2Nzodl1fy25AO++tO3WPDztUzJH8Si66Zi1nO9+ajj8zL5rxuKKK+u54pFr/GdFe9SvG2/5quXLpXIUoLTgE3uvgXAzJYBc4ANMXXmAPcHr38JPGGRn6o5wDJ3rwG2BmvKTgP+mpzmH3GwspbPLTn6sO5+TL1jSo6tEq/omGPFr9P8OHHOH+988Q6WhPPHa0P8Ookcq+1rSeT7Fq9e/O9J29cbr7C6voG6hsiG3MxUbplRwPyLTyIjNTz9i8JRg1m+4HwWvvwPnnp9K0/+eQsAGakDyMlIJSs9hQFmmIEBA4IXRmRMfs//upKOeCSY2O6OR//UYp1TPp7H49dNTfq5Ewn0I4EdMe9LgbNbqhMsJn4IGBaUv9Fs35HNT2Bm84B5AKNHj0607UdJGWDxZySM81PRvCheTy/eD1PzavHrWJt14rep2X4Jtbvt48SrF79jm8B+8fY6pk7bx4l3rER7222dLz11ACOHZFFwXA5TRw8OVYCPNWroQBZ/8QwOVtby5tb9bNlbwf6KWspr6qmubcCBRnfcaXqNx+88SO+QlRb5LBYcn9NinVFDsrrk3AktDt7V3P1J4EmAoqKiDn2SczPTWPylM5LaLpGuNnhgOp+a9LGeboZ0h/+KdER/+KUzu/3UidyM3QnEDj7OD8ri1jGzVGAQsC/BfUVEpAslEuhXAwVmNtbM0oncXF3erM5y4Mbg9eeAP3gkybocuDYYlTMWKADeTE7TRUQkEW2mboKc+wLgFSAFWOruJWb2AFDs7suBp4CfBDdb9xP5ZUBQ7xkiN27rgfnu3tBF1yIiInEklKN39xXAimZl98a8rgY+38K+DwEPdaKNIiLSCXoyVkSkj1OgFxHp4xToRUT6OAV6EZE+zuI9ot6TzGwP8M+ebkcHDAf626Kguub+QdfcO5zo7nGXTwtdoO+tzKzY3Yt6uh3dSdfcP+iaez+lbkRE+jgFehGRPk6BPnme7OkG9ABdc/+ga+7llKMXEenj1KMXEenjFOhFRPo4BfpOMrPvm9k/zGy9mf3KzAbHbOuTC6Ob2efNrMTMGs2sqNm2vnrNM4Nr2mRmd/V0e7qKmS01s4/M7J2YsqFm9qqZvR/8O6Qn25hMZjbKzP5oZhuCz/StQXmfumYF+s57FZjs7lOA94BvQp9fGP0d4LPAn2ML++o1B9ewGLgcmAhcF1xrX/Q0kf+7WHcBv3f3AuD3wfu+oh643d0nAucA84P/2z51zQr0neTuv3X3+uDtG0RW0YKYhdHdfSsQXRi913P3d919Y5xNffWapwGb3H2Lu9cCy4hca5/j7n8msqZErDnAj4PXPwau6s42dSV33+3ua4LXh4F3iaxr3aeuWYE+ub4MvBS8jreo+jELo/cxffWa++p1Jep4d98dvP4AOL4nG9NVzGwMMBX4G33smkOxOHjYmdnvgHgrON/t7i8Ede4m8mfgz7qzbV0lkWuW/sfd3cz63JhsM8sBngNuc/cyM2va1heuWYE+Ae5+aWvbzWwuMAuY4UceTOjVC6O3dc0t6NXX3Iq+el2J+tDMPu7uu83s48BHPd2gZDKzNCJB/mfu/nxQ3KeuWambTjKzmcCdwGx3r4zZ1B8XRu+r17waKDCzsWaWTuSG8/IeblN3Wg7cGLy+Eegzf9FZpOv+FPCuuz8as6lPXbOejO2kYEH0DGBfUPSGu98cbLubSN6+nsifhC/FP0rvYmafAR4HRgAHgXXuflmwra9e86eBHwApwNJgLeQ+x8z+L3ARkWl6PwTuA/4f8AwwmsgU4te4e/Mbtr2SmZ0PvAa8DTQGxf+LSJ6+z1yzAr2ISB+n1I2ISB+nQC8i0scp0IuI9HEK9CIifZwCvYhIH6dALyLSxynQi4j0cf8fIEdxBx7de2wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.stats\n",
    "\n",
    "dom = np.linspace(-25,25,10**3)\n",
    "standardnormpdf = scipy.stats.norm(0,1).pdf\n",
    "\n",
    "fig,ax = plt.subplots()\n",
    "ax.plot(dom,standardnormpdf(dom),label = \"Null Distribution N(0,1)\")\n",
    "\n",
    "ax.axvline(W,color=\"red\",label=\"Our Test stat W\")\n",
    "\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whoa, if the null hypothesis was true our test statistic W would be **incredably** unlikely. \n",
    "It must be that the null hypothesis is false.\n",
    "But what are the chances that we are wrong to reject the null hypothesis? \n",
    "We can compute the probability that the null distribution would generate our test stat $W$ or value more extreme and use this probablity as measure of the probability that the null hypothesis is true and the null distribution generated  our test stat $W$. \n",
    "\n",
    "We want to compute $p( [Z > W] \\cup p(Z < -W) ) = p( Z > W ) + p(Z < -W)$, and because the Normal distribution is symmetric $ p( Z > W ) + p(Z < -W) = 2*p( Z > W )$. This is the **pvalue**.\n",
    "\n",
    "Small pvalues indicate that the probability that the null distribution generated our test statistic, or test statistic more extreme than what we observed, is small indicating the null hypothesis is likely false. Large pvalues suggest it is likely the null distribution generated our test stat, providing no evidence that our null hypothesis is false."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MLE for Beta 1 is 3.81\n",
      "MLE for Beta 0 is 0.29\n",
      "MLE for sigma2 is 16.04\n",
      "MLE for V(b0) is 0.16\n",
      "MLE for square root of V(b0) is 0.40\n",
      "Test statistic (W) is 0.73\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.9/site-packages/seaborn/_decorators.py:36: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1c4345880>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAex0lEQVR4nO3dfXBV9ZkH8O9zQzASXvNiwIQQI6mIkioTFVvpdqAvlNKyagV11rFbp7GzWnBtd9U6O/tHR0fGDrum2q64dlddW6FLHa1SfMHuilO0DYpURCXGgLABYqS8JBMScp/9I/eGm5x7c9/OOb/fuef7mek0yQ0nT26c85zf8/xeRFVBRESUKGI6ACIisg+TAxEROTA5EBGRA5MDERE5MDkQEZHDONMBuKGiokLr6upMh0EB197VAwCoryw1HAmRP7Zv3/6JqlYme60gkkNdXR1aW1tNh0EBt/LhbQCA9TdfbjgSIn+IyN5Ur7GsREREDkwORETkwORAREQOTA5EROTA5EBERA4FMVuJiChsolFFR3cPDh3rQ9XkEtSVlyISEdeuz+RARBQw0ahi866DuH3DDvQNRFFSHMHaFRdhyQXTXUsQLCsREQVMR3fPcGIAgL6BKG7fsAMd3T2u/QwmByKigDl0rG84McT1DURx+Hifaz+DyYGIKGCqJpegpHjk7bukOIKzJpW49jOYHIiIAqauvBRrV1w0nCDiPYe6cvf2BWNDmogoYCIRwZILpmPOqoU4fLwPZ00qoNlKIjITwOMAqgAogHWq+oCIlAFYD6AOQAeAFap6xFScREQ2ikQE9ZUTUV850Zvre3LVzJwC8ANVnQtgAYBbRGQugDsBbFHVBgBbYp8TEZGPjCUHVe1U1TdjHx8HsBtANYDlAB6LfdtjAP7aSIBERCFmRUNaROoAXAzgDQBVqtoZe+kghspOyf5Ns4i0ikhrV1eXP4ESEYWE8eQgIhMBbARwm6oeS3xNVRVD/QgHVV2nqk2q2lRZmfQgIyIiypHR5CAixRhKDE+q6m9iXz4kIjNir88AcNhUfEREYWUsOYiIAHgUwG5VXZvw0rMAbox9fCOAZ/yOjYgo7Eyuc/g8gBsA/FlEdsS+9iMA9wHYICI3AdgLYIWZ8IiIwstYclDV1wCkWrGx2M9YiIhoJK6QJiIreH0+AWWHyYGIjPPjfALKjvGprEREfpxPYItoVNHedQLbPvwE7V0nEI0mna1vHEcORGTcWOcTeLV3kAlBGiFx5EBExvlxPoENgjRCYnIgIuP8OJ/ABn6c4OYWlpWIyDg/ziewQXyElJggbB0hceRARFaIn0+woL4C9ZUTjSYGr5rGQRohceRARJTAy6ZxkEZIHDkQkVG2Te30umls0whpLBw5EJExNk7tDMu02nQ4ciAiY2yc2hmWabXpMDkQkTE2Tu30s2lsW0ktEctKRGSMjVM7/Woa21hSS8SRA1HA2fz0mY7XT+m5vjd+NI1tLKkl4siBKMBsf/pMx8undNvfG9sb3xw5EAWYiadPt0cqXj2l2/5kbnvjm8mBKMD8bujGn8aXtmzFdY+8gaUtW7F510ErS1k2NrsT2b5ammUlogDzuqE7+nQ2VSR9Gp+zaqFvpZBMT4yzsdmdyPbV0kwORAEWf/ocXVd34+kzWc3+J9/6rNE6eTZ9BC/fG7fES2o29BhGY3IgCjAvnz6T1ez3HD5u9Gk8VR8h2cjF9idz27HnQBRw+TR0x2ouJ6vZb2jdj3uvnGesTp5tHyEo+xjZiCMHopBKV6JJVrM/0tuP+bVTscnQ07jtfYRCwpEDUUilm+qZajZNbVmpsadx22f4FBKOHIhCKt0iLBtr9jbGVKiYHIhCKpMSjY2zaWyMqRCxrEQUUizROAV5nyq3ceRAFFIs0Yxk+15MfuPIgSjEONXzNNv3YvIbkwMRAWBJxfa9mPzGshIRsaQCrqEYjcmBiNDR3YM1m3fjpivqIbFcsGbzbsyZPinnWUGZbpBniyDsxeQnJgciQnfPSaxsqkXLK3uGb4yrFjXg056TOSWHII5E2KAfiT0HIsL4oshwYgCGau0tr+xBcVFut4igNnfZoD+NyYGI0Ns/mLQZ29s/mNP12NwNPpaViChlM7Zqcm7NWLeau0HrWxQSjhyIyPXV0m5cL0hHkhYijhyIAsSrJ2m3m7FuXC+bg33IfUaTg4j8AsAyAIdV9cLY18oArAdQB6ADwApVPWIqRiJbeD0DyO0N7fK9XrpdY8lbpstK/wlgyaiv3Qlgi6o2ANgS+5wo9GyfAeT2Cut43yJRmBel+c3oyEFVXxWRulFfXg7gi7GPHwPwPwDu8C8qIjt195wcsUht4/b96DzaZ8WTtBejGrcXpbG5nR0bew5VqtoZ+/gggKpk3yQizQCaAaC2ttan0IjMiEYV//eXPjz6WvuIRWrrW/dZ8STtRX/AzT5IEBflmWa6rDQmVVUAScemqrpOVZtUtamystLnyIj81dHdgzs27nQsUvvx8nlWbO/g1boGtxal2V6Ss5GNyeGQiMwAgNj/HzYcD5FxqW6+xUVixZOv7f0BLsrLno3J4VkAN8Y+vhHAMwZjIfJcJo3cVDffXBepuc32U+VsT142Mj2V9VcYaj5XiMh+AP8M4D4AG0TkJgB7AawwFyGRtzKthdu+Y6jtm9bZ/v7ZSIbK+sHW1NSkra2tpsOggFv58DYAwPqbL/ftZ7Z3ncDSlq2ObSbWNy/AvOqpI26u8dk2Nt58g4Dvn5OIbFfVpmSv2ThbiSg0UtXCt7x3GAf+0jdiBOH2IrWw4fuXHRt7DkShkaoWPhhF0tk0YTrKM0y/q404ciAyKFktfNWiBjzx+l7HVhFhmqsfpt/VVhw5EBkUb+Sub16AVYtn46Yr6vHE63vRebTPMZsmTHP1w/S72orJgciwSEQwr3oq5kyfjEdfax9ODKNn06TqT3xw6HjBlV24LsE8lpWILJDJVNBUB+gMDCqe2XEA82qmYNF5Vb6WXbzar8itw4Iodxw5EFki3VYRyRaarV7cgHs37cbDr7Zjz6ET2Pepf2UXLw/jsX1RXRhw5EAUEImjiz2HTmBX51E8vm2oPwEAD2zZg/m101BX4c9UTS8P47F9UV0YMDkQBUh8dNHR3YOWLW0jXusbiKK3/5RvsXh9GA/XJZjFshJRAM0qK026PqK2zL+yC/crKmxMDkQBdE5F8pr8ORX+JQf2BQoby0pELvLrtDEbavI2xEDeYXIgcolXq3pTJRwbavI2xEDeYHIgckl89s60CeNx1fwaTCopQjSq+N8PDmNWeWlOT9UmtpHgWcsEMDkQuebQsT5MmzAeNyyYhfWt+7CyqRY//O+387qpezldNJlTp6J4/p3O4SNJuadReLEhTeSSqskluKapBi2v7MGyxmq0vLInp72BEncj7Tp+EtMmjB/xulfbSESjij+0dzvOquaeRuHEkQP5ptDLFXXlpfjMWZPQNxCFCHJaA5CsjLR6ccOIxW5eTRft6O5B695PPV27QMHBkQP5wsutFmwRiQjOnzF5xNTORJnc1JOVkR7YsgfXNNUMX8Or6aKHjvUhqrnFTYWHyYF8EZYtmOPrD3779gGsWtSQ9RqAVKuOL545FU81X4ZNqxZ6Vv+vmlySNO41Vzdy7UIIsaxEvvB6qwVbDM/9nz4Jn/acxPrmBejtH8y4jJZqN9JZ5aWev0915aW4Y8n5WLN5N266oh5FEaBpVhk+V19eUOU/ygyTA/kiTFsw5zP3P9nJcH6tOk5MbJksaiv0HlLYMTmQL0ze9ILE9KrjTBMbj/EsfEwO5AvTN70gCcKqY7/XX5D/mBzIN0G46VFmwtJDCjMmB/JVIdapR/9OtdMmYN+R3oL6HUdL7CHNmFKCq+bXoCgCnFk8bnh6cqH9ncOGyYF8E5Q6dTYJbPTvNKv8THx/UQN+GlslXRQBLplVhsvryzFuXOHMHI/3kNZs3o2VTbXDq8HXvdqOtSsuwvhxglt/+ZbVf2caW+H810rWC8Jah76BQWz6cyee2XEAt294O+1ivdG/07LGavz0lT1Y2VSLR19rR8uWNnz3iVY8/05nwS34W3LBdLRce3HSbUJ27j9q9d+Z0uPIgXzjR50617JVNKr4tKcfbV0ncOuv3kJJcQR3LZmD4ydP4b2Dx1A99UzMq57iuNbo30kESfdVumPjTsyrnlJQ9fhIRNDbP5j0bzo6D7IfETxMDuQbt9Y6pEoA+ZStOrp70NZ1Ahq7qU2bMB69A4N48PdtI8olo6+V7HcqiuS2r5Ib74HfUv1NR4dSqGtaChnLSuQbN46VHGuPpnzKVoeO9Q0nBgC4an4NHtiSflfV0b/Tb98+gItnTvV0fyKb9qlK9TdtrJnC40MDjiMH8o0bax3Gml+fT9mqanIJRDCcIDLdVTXZ71Qz5UysubpxxJkI9145DxEZurHn+4Rv0xqDVH9TANjENS2BxuRAvsp3rcNYCSCfslVdeSlmV04cLi0VCTK+VrLf6RuNZ+PCs6dg98Fj+ODQcdz/wvs40tvvyqwd29YYpPqbck1LsLGsRIESTwCJ4jftfMpWkYigrHQ8Gqun4Knmy3DlxdV5lcAiEYEI8MNfv42WLW3oPNrn2qydsd4DIrdw5ECBMtYeTW6UrUqKi7CgvgIAUFtWmte1vHrC5z5V5AcmB/KFW7Nr0iUAN7foyPdaXu1Ey32qyA9MDuQ5t1dGB2WPJi+f8IPyHlBwMTmQ52yaXeMnPuFTkKVsSIvIJhGp8zEWKlBj1d4LXfwJf0F9BeorJzIxUGCMNVvpPwC8KCJ3i0ixXwHFicgSEXlfRNpE5E6/fz65h7Nr8heNKtq7TmDbh5+gvetEQe3TRHZKWVZS1V+LyO8A/BOAVhF5AkA04fW1XgUlIkUAHgLwZQD7AfxJRJ5V1Xe9+pmUuWyby5nW3m3ZEsI2QdnNlgpLup5DP4AeAGcAmISE5OCxSwG0qWo7AIjIUwCWA0iaHNq7erDy4W0+hUbxDepUh1YSz66ciLLS8Wn/3XlVk9A/GMX4ogge+0MHHvtDhyvXdcu7nccAwLr/lvoGBrHzwNHh1dt9A1Hc8ss30Vg9BSXFRWaDo4KVMjmIyBIAawE8C2C+qvb6FhVQDeDjhM/3A7gs8RtEpBlAMwBMnHGuf5GFXN/A4IgN6lSBtq4TaByf/kZVUlyU8nvyuW6h6x+Mjtj3CRh6f/oHo6F/b8g7Y40c7gZwjaru8iuYbKjqOgDrAKCpqUnX33y54YjCYduHn+C6R94Y8TVV4K6lc4YXj2UqsYwUEcFt63eg8+jpJnWu181VfMRg239L7V0nsLRlq2O9xEPXzy/o2V7kvQ3fS/3aWD2HhV4Ek6EDAGYmfF4T+xoZlsvCrmS9BACOOvrqxQ14fNve4QTBpvUQrogmE2xd5/AnAA0icg6GksK1AK43GxIB2d+oUjVT586Y5Fj78MCWPWj+Qj1atrTxBpiA6yXIBCuTg6qeEpFbAbwAoAjAL2wtb4VNtjeqVAvgHvvbS5Oufbh45lQ81XxZ4G+Abs+84opo8puVyQEAVHUTgE2m4yCnbG5UqRbA9fafSlqemlVeGvgbIKeeUiHglt2Ut7EWaKVaAFdblv+pcLbK50Q6IltYO3KgYEj3lJyqR3FORSnOqchvS2xb2XYYD1EumBwoL+k21UvXo8injm7rimqvtuom8hOTA+Ulk6fkxB6FWzf0VCOWr5xfhX1Heo0mDE49pULA5EB5yeYp2c1GbaoRy7obmtD8RKvRRjCnnlIhYEOa8pLNuc1uNmpTjVha936a0fXHaqK7sftpsq26ubMqBQlHDpSXbJ6S3WzUphqxDI7aGjLZ9VONYFSBI739w1tVZDPySFcuy3XUZGtfhQofRw6Ut0wPtHHzXIdkI5Y1VzfiuZ0jd1lJdv1UI5je/lNo6zqR9cgmfuNf2rIV1z3yBpa2bMXmXQdHjAxyGTVlcl0irzA5kG+yKUGlK8HERyybVi3EU82XYdOqhfj6hTNwx5Lz014/1Qjm5Cnn7qeZnFiXyY0/l9PwuF6CTGJZiXyTaQkq0xJMspXamVw/VUlqXJFABCMSRCYjm0zKZblMb+V6CTKJIwfyVSYlqGyemEePMACkvX6yEczqxQ3Y192LmqlnZr1qO5NyWTajpmyuS+QVjhzIOmM9MdeVlw43aGdMKcG7ncezbvImjmD2dvfgrY//gse37UVP/yAGBvvwvS+ei5qpZ+K86ZMwr3pq2gZwJusacpneyvUSZJLo6CJrADU1NWlra6vpMALBxtkvo2OKCLDkAefhNptXLxyRDFYtno11r7Y7vm9TbHV2JpIdXgQAty6ajStml2d80FD8d3B7XYNX1yUCABHZrqpNyV7jyCFE4rX8NZt3Y1ljNYoiwCWzynB5fTnGjTNTYUzVX3jw+ovx4+feHRHnqUEdUW6KKvKuyafqBUQEWZVvvNpSm1t1kylMDiHS0d2DNZt3Y2VTLVpe2TN8M15zdSO+0Xi2r0+k8Sfiju4evH/wGKZNGI/Oo33D/YXnbr0Ct3/5PNyxcedwnD/51mcdySDfU+lqp00YUboRAVYvbkBD1USWbyjUmBxC5NCxPixrrB5ODMDQk/YdG3diXvUU355Ok40WVi1qwBOv7x1OEB8f6R1ODPE49xw+PiIZbNy+H6sXN+CBLXsyqsmPtR/T899fiL97cjuKIoKvXTgdtWUs31C4cbZSiFRNLkFRJHUpxi/JZiO1vLIHV82vATD09D9h/DhHnBta9+PeK+cNz+A50tuPhqqJeP77p9c6jNWMTjULat+RXpx71kRMnTAek0qKUVeReiEfUVhw5BAideWluGRWmfHtpFPNRhI5PcWzavIZjjiP9PZjfu1UbEoy4+fcs3I/lY7rBoicOHIIkUhEcHl9OdZc3Wj0BLZU8/cXzq4YfvpPdVJcbVlpRlt1ZPNzuW6AyIkjh5AZNy6CbzSejXnVU4xNj0w1f/+SurIRcbi97TXXDRBljskhhNyaHpnrmolMF4S5PY2T5ywQZY7JgXKS78E9pubvc90AUWbYc6Cc5LtjKA++IbIbRw6Uk3xm/iSOOqZNGI9rmmrwmbMm4fwZk3FOBcs8RDZgcgiwXGr+bu2tlMsW1HHxUce0CeNxw4JZI1ZrJytN2bgfFFGhY3KwQK43+Wxr/vn2CRLlM/MnPuq4an6NY7X27Rt2YE7CxnluxkxEmWNyMCzXm1+qmv+cMXYkzeXfpJLPzJ/4qOOMcZGkpalDx06XptyMmYgyx4a0Ybk2dnM5djKXfzOWTM+OHi0+6qivLE26KG3C+CLPYiaizDA5GJbrzS+X1b62rBCOjzqmnDkOqxY1jFgFvWpRAwYGT78ftsRMFDZMDoblevPL5djJXP6NVyIRQfXUCVjfug83XVGPWxfNxk1X1GN96z6UlZ5hZcxEYcKT4AzLp+Gayylhfp4slq7RfupUFM+/0znizIaxZit5HfPKh7cBANbffLnr1yayEU+Cs1g+jd1cVvv6tUI4XdKLRhUv7j6EtS+9j5uuqEdRBGiaVYbP1Zd7vo0GEaXH5GCBQrz5pZtllPj6Q79vA5D9+c9E5B32HGK4nYO70jXaOQuJyG4cOYALrbyQagV15cSSMV/nLCQiO3DkgPw3kSOnZLOMVi9uwEfdQ6MyzkIishtHDuDxkV6IRARzZ0xC8xfqEVVAFXh8214c6e0f7ivwbAUiezE5gCUOr3Qe7UPLljbH1+NJtxAb8USFwkhZSUSuEZFdIhIVkaZRr90lIm0i8r6IfNWPeFjiyE26Jj5XNxMFl6mRwzsArgLwcOIXRWQugGsBXADgbAAvi8hnVHXQy2B4fGT2Mmni88xmouAykhxUdTcAiDhuvssBPKWqJwF8JCJtAC4FsM3rmPwucQT9jIJMdktl0iUKLtt6DtUAXk/4fH/saw4i0gygGQBqa2u9jyxLY938C2HqbKZNfPYViILJs+QgIi8DmJ7kpbtV9Zl8r6+q6wCsA4b2Vsr3em5Kd/P38owCv0YkbOITFTbPGtKq+iVVvTDJ/8ZKDAcAzEz4vCb2tUBJt27Cq9XB8aS0tGUrrnvkDSxt2YrNuw56stqbTXyiwmZbWelZAL8UkbUYakg3APij2ZCyl67k4tVTt5+nprGfQFTYTE1lvVJE9gO4HMDzIvICAKjqLgAbALwLYDOAW7yeqeSFdFM4vXrq9nu/olxPgiMi+5marfQ0gKdTvHYPgHv8jchd6aZwevXUzT4AEbnFtrJSQRjr5j+6YXxpnfP8gly5ua4g6FNtiSg/TA4eSTaF0+sprG6NSAphqi0R5Ye7svrIj91f3egDcJdaImJy8FFQDrgxEScPWyKyC8tKPgpKw9jvOFnGIrIPRw4+CsrCMb/jZBmLyD4cOfgoKAvH/I6Thy0R2YfJwWdB2YjOzziDUm4jChOWlci4oJTbiMKEIwcyLijlNqIwYXIgKwSl3EYUFiwrERGRA5MDERE5sKxESXHjPaJwY3IgB65YJiKWlVxQaPsCccUyEYV25OBW2aQQn7K5YpmIQjlyiN/Ql7ZsxXWPvIGlLVuxedfBnJ74C/EpO90xp0RU+EKZHNy8oQdlG+5scMUyEYWyrORm2aQQ9wXiimUiCuXIwc2ySaE+ZbtxohwRBVcoRw7xG/roJnIuN3Q+ZRNRIQplcnD7hs59gYio0IQyOQC8oRMRjSWUPQciIhobkwMRETkwORARkQOTAxEROTA5EBGRA5MDERE5MDkQEZEDkwMRETkwORARkQOTAxEROTA5EBGRA5MDERE5MDkQEZEDkwMRETkYSQ4icr+IvCciO0XkaRGZmvDaXSLSJiLvi8hXTcRHRBR2pkYOLwG4UFUbAXwA4C4AEJG5AK4FcAGAJQB+JiJFhmIkIgotI8lBVV9U1VOxT18HUBP7eDmAp1T1pKp+BKANwKUmYiQiCjMbeg7fAfC72MfVAD5OeG1/7GsOItIsIq0i0trV1ZV3ENGoor3rBLZ9+Anau04gGtW8r0lEFFSeHRMqIi8DmJ7kpbtV9ZnY99wN4BSAJ7O9vqquA7AOAJqamvK6k0ejis27DuL2DTvQNxBFSXEEa1dchCUXTM/5XGkioiDzLDmo6pfGel1Evg1gGYDFqhq/uR8AMDPh22piX/NUR3fPcGIAgL6BKG7fsANzVi3kGdNEFEqmZistAfCPAL6pqr0JLz0L4FoROUNEzgHQAOCPXsdz6FjfcGKI6xuI4vDxPq9/NBGRlTwbOaTxIIAzALwkIgDwuqp+T1V3icgGAO9iqNx0i6oOeh1M1eQSlBRHRiSIkuIIzppU4vWPJiKykpHkoKqzx3jtHgD3+BgO6spLsXbFRY6eQ115qZ9hEBFZw9TIwSqRiGDJBdMxZ9VCHD7eh7MmlaCuvJTNaCIKLSaHmEhEUF85kQ1oIiLYsc6BiIgsw+RAREQOTA5EROTA5EBERA5MDkRE5CCnd64ILhHpArA3z8tUAPjEhXCCjO8B34M4vg/heA9mqWplshcKIjm4QURaVbXJdBwm8T3gexDH94HvActKRETkwORAREQOTA6nrTMdgAX4HvA9iOP7EPL3gD0HIiJy4MiBiIgcmByIiMiByWEUEfmBiKiIVJiOxQQRuV9E3hORnSLytIhMNR2TX0RkiYi8LyJtInKn6Xj8JiIzReT3IvKuiOwSkdWmYzJJRIpE5C0Rec50LCYwOSQQkZkAvgJgn+lYDHoJwIWq2gjgAwB3GY7HFyJSBOAhAF8DMBfAdSIy12xUvjsF4AeqOhfAAgC3hPA9SLQawG7TQZjC5DDSv2DobOvQdulV9UVVPRX79HUANSbj8dGlANpUtV1V+wE8BWC54Zh8paqdqvpm7OPjGLoxVpuNygwRqQHwdQD/bjoWU5gcYkRkOYADqvq26Vgs8h0AvzMdhE+qAXyc8Pl+hPTGCAAiUgfgYgBvGA7FlH/F0INiNM33FaxQnQQnIi8DmJ7kpbsB/AhDJaWCN9b7oKrPxL7nbgyVGZ70MzYyT0QmAtgI4DZVPWY6Hr+JyDIAh1V1u4h80XA4xoQqOajql5J9XUTmATgHwNsiAgyVUt4UkUtV9aCPIfoi1fsQJyLfBrAMwGINz0KYAwBmJnxeE/taqIhIMYYSw5Oq+hvT8RjyeQDfFJGlAEoATBaR/1LVvzEcl6+4CC4JEekA0KSqhb4jo4OILAGwFsBfqWqX6Xj8IiLjMNSAX4yhpPAnANer6i6jgflIhp6MHgPwqareZjgcK8RGDj9U1WWGQ/Edew402oMAJgF4SUR2iMi/mQ7ID7Em/K0AXsBQI3ZDmBJDzOcB3ABgUexvvyP29EwhxJEDERE5cORAREQOTA5EROTA5EBERA5MDkRE5MDkQEREDkwORC6L7W76kYiUxT6fFvu8znBoRBljciBymap+DODnAO6Lfek+AOtUtcNYUERZ4joHIg/EtqHYDuAXAL4L4CJVHTAbFVHmQrW3EpFfVHVARP4BwGYAX2FioKBhWYnIO18D0AngQtOBEGWLyYHIAyJyEYAvY+hEtb8XkRlmIyLKDpMDkctiu5v+HEPnIewDcD+An5iNiig7TA5E7vsugH2q+lLs858BOF9E/spgTERZ4WwlIiJy4MiBiIgcmByIiMiByYGIiByYHIiIyIHJgYiIHJgciIjIgcmBiIgc/h+KXlnrKJq/rwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2d0lEQVR4nO3deXxU1d348c83k5WQBLKwJpAQwhIQEwkBwqKIC1YErUvR+iit1ccW61K11WqtS7XYWrUqfai/Rx/bqqW4VaoIrqgIypaAhDUJW9iyQBayL+f3x0ziEBIyIZPcmcn3/XrNizvnnnPnOyH5zp1zzz1HjDEopZTyXX5WB6CUUqpraaJXSikfp4leKaV8nCZ6pZTycZrolVLKx/lbHUBL0dHRJj4+3uowlOo+O3fa/x050to4lFfbuHFjkTEmprV9Hpfo4+Pj2bBhg9VhKNV9zjvP/u+qVVZGobyciOxra5923SillI/TRK+UUj5OE71SSvk4j+ujb01dXR35+flUV1dbHYryMcHBwcTGxhIQEGB1KEp1Ga9I9Pn5+YSFhREfH4+IWB2O8hHGGIqLi8nPzychIcHqcJTqMi513YjILBHZKSI5InLfaepdKSJGRNKcyu53tNspIhefSZDV1dVERUVpklduJSJERUXpN0Xl89o9oxcRG7AIuBDIB9aLyDJjzLYW9cKAO4BvnMqSgXnAGGAQ8LGIjDDGNHQ0UE3yqivo75XqCVzpukkHcowxeQAisgSYC2xrUe8x4EngXqeyucASY0wNsEdEchzHW9vZwJXqLtV1DazaWUBOwQkC/f1IietL2tC++Pnph4TyDq4k+sHAAafn+cBE5woicg4QZ4x5X0TubdH26xZtB7d8ARG5BbgFYMiQIa5FrlQXa2w0vLJmL899upuSyrqT9o3o35tH5oxlcmKURdEp5bpOD68UET/gaeDuMz2GMeZFY0yaMSYtJqbVO3gtJyLcffd3b/Gpp57i4YcfPm2bvXv3MnbsWABWrVrF7NmzT6mzatUqIiIiSE1NZeTIkUyfPp333nuvef/ixYv5+9//3uZrrFq1ijVr1rS5f9myZSxcuBCA+fPn8+abb5425paeeOKJk55nZGR0qH1b5s+fz+DBg6mpqQGgqKgI56kvDh8+fNLP6/e//z3Dhw9n5MiRrFy5stVjPvDAA8TFxdG7d++Tyl944QVefvnlDsVXWVvPj/+2nkff28a42D68/pOJ7HhsFlkPXcgzPzibqroGrv1/X/PiF7no4j3K07mS6A8CcU7PYx1lTcKAscAqEdkLTAKWOS7IttfWawQFBfH2229TVFTk9mNPmzaNzMxMdu7cyXPPPcdtt93GJ598AsCtt97KDTfc0Gbb0yX6+vp65syZw333tXn9vF0tE/3pPlQ6ymaztZmAn376aW6++WYAtm3bxpIlS8jOzmbFihX87Gc/o6Hh1Ms8l112GevWrTul/Mc//jHPP/+8y3FV1zUw/+X1fLGrkMcuH8vffjSBjOHRBAfY6NMrkCtSY/nwznO59KyBPLF8By98muPysZWygitdN+uBJBFJwJ6k5wHXNe00xpQC0U3PRWQVcI8xZoOIVAGvi8jT2C/GJgGn/iV2wCP/yWbbobLOHOIUyYPC+e1lY05bx9/fn1tuuYVnnnmGxx9//KR98+fPZ/bs2Vx11VUA9O7dmxMnTpxRLCkpKTz00EO88MILzJw5k4cffpjevXtzzz338Nxzz7F48WL8/f1JTk5m4cKFLF68GJvNxquvvsrzzz/PSy+9RHBwMJmZmUyZMoVx48axYcMGXnjhBQA+/vhjFi5cSFlZGU8//TSzZ8/mlVdeOanO7Nmzueeee1ixYgVVVVWkpKQwZswYXnvtteb3Zozhl7/8JR988AEiwoMPPsgPfvADVq1axcMPP0x0dDRbt25l/PjxvPrqq61e9Lzzzjt55plnmhO6s7feeovf/e53ALz77rvMmzePoKAgEhISGD58OOvWrWPy5MkntZk0aVKrP9NevXoRHx/PunXrSE9PP+3P3xjDvW9uYd3eYzx3bSpzzh7Uar2QQBvPX5tKoL8ff/poF7GRIVyRGnvaYytllXbP6I0x9cBtwEpgO7DUGJMtIo+KyJx22mYDS7FfuF0BLDiTETeeYsGCBbz22muUlpZ26eucc8457Nix45TyhQsXkpmZyZYtW1i8eDHx8fHceuut3HXXXWRlZTFt2jTAft/BmjVrePrpp085xt69e1m3bh3vv/8+t95662mHFi5cuJCQkBCysrJ47bXXTtr39ttvk5WVxebNm/n444+59957OXz4MACZmZk8++yzbNu2jby8PL766qtWjz9kyBCmTp3KP/7xj5PK9+zZQ9++fQkKCgLg4MGDxMV998UwNjaWgwc79sUwLS2NL7/8st16r32zn/9sPsS9F49sM8k38fMTnrxyHBMTIvn121vJLTyzD3eluppLN0wZY5YDy1uUPdRG3fNaPH8ceLy1umeivTPvrhQeHs4NN9zAc889R0hISJe9Tlt9vuPGjeOHP/whl19+OZdffnmb7a+++mpsNlur+6655hr8/PxISkpi2LBhrX6guGL16tVce+212Gw2+vfvz7nnnsv69esJDw8nPT2d2Fj72W1KSgp79+5l6tSprR7n/vvvZ+7cuVx66aXNZYcPH8bd12r69evX7nvdW1TB4+9vZ1pSND89N9Gl4wb6+/Hnealc8ucvuOtfWbzzsynYdDSO8jA6100H3Xnnnbz00ktUVFQ0l/n7+9PY2AhAY2MjtbW1nXqNzMxMRo8efUr5+++/z4IFC9i0aRMTJkygvr6+1fahoaFtHrtlF4qInBQ/0OkbiJrOxMHeD99WnABJSUmkpKSwdOnS5rKQkJCTYhg8eDAHDnw38Cs/P5/Bg08ZvHVa1dXVp/1wNsbw63e+JcAm/OGqcR0aOjkgIphH5o5lS34pr33T5kyxSllGE30HRUZGcs011/DSSy81l8XHx7Nx40bAPsqlrq6urebt2rJlC4899hgLFiw4qbyxsZEDBw4wY8YMnnzySUpLSzlx4gRhYWGUl5e7fPw33niDxsZGcnNzycvLY+TIkcTHx5OVldX8Gs4XNAMCAlp9P9OmTeNf//oXDQ0NFBYW8sUXX7Tb/92WBx54gKeeeqr5+YgRI9i7d2/z8zlz5rBkyRJqamrYs2cPu3fvbn6tmTNnutSNs2vXruYRUK35cNtR1uQWc+/FIxkY0fFva5eNG8jU4dH8ccVOCsr1TlvlWTTRn4G77777pNE3N998M59//jlnn302a9euPe0ZdWu+/PLL5uGVCxYs4LnnnmPmzJkn1WloaOD666/nrLPOIjU1ldtvv50+ffpw2WWX8c4775CSkuJSH/SQIUNIT0/nkksuYfHixQQHBzNlyhQSEhJITk7m9ttv55xzzmmuf8sttzR3GTm74oorGDduHGeffTbnn38+f/jDHxgwYECH3neTMWPGnPSaoaGhJCYmkpOT07z/mmuuITk5mVmzZrFo0SJsNhuNjY3k5OQQGRkJwC9/+UtiY2OprKwkNjb2pOGvX331FRdeeGGrr2+M4Ynl20nq15tr08/sPg4R4ZG5Y6isa2CRjsJRnsYY41GP8ePHm5a2bdt2SpnybW+//bZ54IEHTlvn22+/NXfddVe7x9q0aZO5/vrr29y/PnOLGfqr98ynO452OM6W7ntrixn+6/fNgWMVrjc691z7Q6lOADaYNvKqntErj3TFFVfQ3trBY8eObXVkUUtFRUU89thjre5rNIby6npS4vpw3ojOXwC+feZwRITnPtnd6WMp5S6a6JXH+slPfuKW41x44YVtfmiUVNZR32iaE3RnDYwI4br0Iby96SBHSrWvXnkGTfSqxzLGUFheQ6BNmDGyn9uOe9PUBBqNfZ4cpTyBJnrVY5VX11NT30Dv4AC3TlccF9mLS8YO5PVv9nGipu2hpUp1F030qscqrqglwOZHSID7/wx+Mi2Bsup6lq4/0H5lpbqYJnoX5efnM3fuXJKSkkhMTOSOO+7o1I1RCxYsICUlheTkZEJCQkhJSSElJcXl2SVLSkr4y1/+0qHXfPbZZ6msrDzjeu++++5Jd+Q2zSjZ5D//+Q9z5px2VgyPUVPXQHl1HZGhgV2y+EjqkL6MH9qXV7/Zp7NbKstponeBMYbvf//7XH755ezevZtdu3Zx4sQJHnjggQ4dx3nGxUWLFpGVlcXy5ctJTEwkKyuLrKys5onR2mNFos/IyODrr79bXmDt2rWEh4dTUFAA2Ge2dNc0xl2tuKIWESEyNLDLXuPa9CHkFVawbs+xLnsNpVyhid4Fn376KcHBwfzoRz8C7Lf1P/PMM7z88stUVlbyyiuvcNtttzXXnz17NqtWrQLsM1nefffdzTdTnU5FRQU//vGPSU9PJzU1lXfffReA7Oxs0tPTSUlJYdy4cezevZv77ruP3NxcUlJSuPfee085zqWXXsrZZ5/N2LFj+de//sVzzz3HoUOHmDFjBjNmzADgpz/9KWlpaYwZM4bf/va3AK3WaxITE0N4eHjzjUwHDx7kyiuvbJ66eM2aNUyZMuVMfsTdqrHRcLyylojgAAJsXfcncOlZAwkL9uef6/Z32Wso5QqXJjXzKHfeCVlZ7j1mSgo8+2ybu7Ozsxk/fvxJZeHh4QwZMqQ56bWloqKCiRMn8qc//andMB5//HHOP/98Xn75ZUpKSkhPT+eCCy5g8eLF3HHHHfzwhz+ktraWhoYGFi5cyNatW8lq5WexYsUKBg0axPvvvw9AaWkpERERPP3003z22WdER0c3v15kZCQNDQ3MnDmTLVu2cPvtt59Sz9mUKVNYs2YNDQ0NJCUlMWnSJFauXMns2bPZvHkzEyZMaPd9Wq2suo6GRkPf0IAufZ2QQBtXpA5myfoDPFxZS59eXfftQanT0TP6Lmaz2bjyyitdqvvhhx+ycOFCUlJSOO+886iurmb//v1MnjyZJ554gieffJJ9+/a1O3PmWWedxUcffcSvfvUrvvzySyIiIlqtt3TpUs455xxSU1PJzs5m27aWywCfKiMjgzVr1rBmzRomT55Meno633zzDZmZmYwaNYrg4GCX3quVSirrCLD50Tuo689z5k0YQm19I//O9Mr1dpSP8L4z+tOceXeV5OTkUy6SlpWVsX//foYPH86WLVvanP0xODi4zSmDWzLG8NZbbzFy5MiTykePHs3EiRN5//33+d73vsdf//pXhg0b1uZxRowYwaZNm1i+fDkPPvggM2fO5KGHTp5Ves+ePTz11FOsX7+evn37Mn/+fJdmrZwyZQrPP/88DQ0N3HzzzYSFhVFdXc2qVau8on++rqGR8up6osO65iJsS8mDwkkeGM47WYeYPyWhy19Pqda4dEYvIrNEZKeI5IjIKevSicitIvKtiGSJyGoRSXaUx4tIlaM8S0QWu/sNdIeZM2dSWVnZvHZrQ0MDd999N/Pnz29evait2R874uKLL+b5559vHqWRmZkJQF5eHsOGDeP2229n7ty5bNmy5bSzVh46dIhevXpx/fXXc++997Jp0yaAk9qUlZURGhpKREQER48e5YMPPmhuf7pjjx49mkOHDrF69WpSU1MB+5zzixcv9or++ZLKOgyGvt3YjXJ56iA2Hyhhb1FF+5WV6gLtJnoRsQGLgEuAZODapkTu5HVjzFnGmBTgD9gXC2+Sa4xJcTxudVPc3UpEeOedd3jjjTdISkpixIgRBAcHN6+nerrZHzviN7/5DXV1dYwbN44xY8bwm9/8BrB3sYwdO5aUlBS2bt3KDTfcQFRUFFOmTGHs2LGnXIz99ttvmy/ePvLIIzz44IOAfSbKWbNmMWPGDM4++2xSU1MZNWoU11133UlJ2rleaz+LiRMnEhUVRUCAvY978uTJ5OXlecUZ/fHKWnoF+hMc4Nq3LHe47OxBiMC7WYe67TWVcibtjfEVkcnAw8aYix3P7wcwxvy+jfrXAjcYYy4RkXjgPWNM2xOBt5CWlmY2bNhwUtn27dtbXYhDqY6ormtg19FyBvUJIbr3d4ujdMfv17wX11JQVsMnd597apfReefZ/3WM1FLqTIjIRmNMWmv7XOm6GQw4396X7yhr+SILRCQX+xn97U67EkQkU0Q+F5FpbQR4i4hsEJENhYWFLoSkVMeVVtkXUIkI6drRNq25PGUweUUVfHuwa9cbVqo1bht1Y4xZZIxJBH4FPOgoPgwMMcakAr8AXheR8FbavmiMSTPGpLl7rVClmpRW1REa5N+lY+fbcsnYgQTa/LT7RlnCld/4g0Cc0/NYR1lblgCXAxhjaowxxY7tjUAuMOJMAtXbyFVnVNc1UF3XcMrZfHf9XkX0CmBqUjQrth7R32XV7VxJ9OuBJBFJEJFAYB6wzLmCiCQ5Pb0U2O0oj3FczEVEhgFJQF5HgwwODqa4uFj/QNQZa63bxhhDcXFxt439nzVmAAdLqsg+VNYtr6dUk3bH0Rtj6kXkNmAlYANeNsZki8ij2JeuWgbcJiIXAHXAceBGR/PpwKMiUgc0ArcaYzo88UdsbCz5+flo/706U0fLqvETIac86KTy4OBgYmNjuyWGC5L74/c2rNh6hLGDW7+JTamu4NINU8aY5cDyFmUPOW3f0Ua7t4C3OhMgQEBAAAkJerOJOjM5BSeY/7fPefiyZKaPtu73KDI0kIkJUazMPsI9F49sv4FSbqJTICiftzL7CACzxg60OBKYNXYAuwtOkFNwwupQVA+iiV75vI+2HeXs2AgGRFg/D89FY/oD3334KNUdNNErn1ZYXsPm/BJmju5vdSiAffHwlLg+muhVt9JEr3zaZzsKMAZmjnbf4t+ddfGYAWzJL+VwaZXVoageQhO98mkfbz/KoIhgkgeecp+eZZo+dD7boaPIVPfQRK98VnVdA1/uLmLm6P7dMiWxq5L69WZwnxA+3VFgdSiqh9BEr3zW2txiquoaPKrbBuwzgM4YFcNXOUVU1zW030CpTtJEr3zWx9uP0ivQxqRhUVaHcorzR/Wjqq5BFw5X3UITvfJJxhg+3VHAtKTobp173lWTh0UT5O+n3TeqW2iiVz5p++FyDpdWM3OUZwyrbCkk0EZGYhSf7SxAZ3BSXU0TvfJJn++yj2g5d6TnTnt9/qh+7Cuu1H561eU00Suf9MWuQkYNCKN/uPV3w7Zlxij7ReKSyjqLI1G+ThO98jkVNfVs2HeM6SM892weILZvL0b0783xylqrQ1E+ThO98jlf5xVT12CYnuTZiR7g3BExlFfX06BrLagupIle+ZwvdhUSHOBHWnxfq0Np17SkGIwxlFfVWx2K8mGa6JXP+WJ3EZOGRXnksMqW0hMiERFKqrSfXnUdlxK9iMwSkZ0ikiMi97Wy/1YR+VZEskRktYgkO+2739Fup4hc7M7glWrpwLFK9hRVeEW3DUBwgI3wYH9Kq7SfXnWddhO9Y83XRcAlQDJwrXMid3jdGHOWMSYF+APwtKNtMvY1ZscAs4C/NK0hq1RX+GK3fVilp1+IdRYREkBVbQNHy6qtDkX5KFfO6NOBHGNMnjGmFlgCzHWuYIxxXu04FJrvAZkLLDHG1Bhj9gA5juMp1SW+2FXI4D4hJMaEWh2KyyJ6BQLw5e4iiyNRvsqVRD8YOOD0PN9RdhIRWSAiudjP6G/vYNtbRGSDiGzQBcDVmapraGRNTjHTR0R71GyV7ekVaCPA5seXu/V3X3UNt12MNcYsMsYkAr8CHuxg2xeNMWnGmLSYGO/5yq08S9aBEspr6r2mf76JYO+++SqniMZGHWap3M+VRH8QiHN6Husoa8sS4PIzbKvUGftyVyF+AhnDo60OpcMiQgIoOlHL9iNl7VdWqoNcSfTrgSQRSRCRQOwXV5c5VxCRJKenlwK7HdvLgHkiEiQiCUASsK7zYSt1qq9yixkX24eIkACrQ+mwiF72mLWfXnWFdhO9MaYeuA1YCWwHlhpjskXkURGZ46h2m4hki0gW8AvgRkfbbGApsA1YASwwxugMTsrtTtTUs/lACVOGe97c864ItPkxsn8YqzXRqy7g70olY8xyYHmLsoectu84TdvHgcfPNEClXLF+zzHqGw0Zid7XbdNkWlI0f/96H1W1DYQE6ihk5T56Z6zyCWtyiwj092P8UM+f9qAtU5Oiqa1vZOO+41aHonyMJnrlE77KKWb8kL5eMe1BWybER+LvJ6zJ1e4b5V6a6JXXO15Ry7bDZV7bP98kNMiflLg+rMkttjoU5WM00SuvtzbPnhgne3H/fJOMxCi25JdQVq2TnCn30USvvN6a3CJCA22Mi42wOpROm5wYTaOxX1xWyl000SuvtyanmInDogiwef+vc+qQPgT5+2n3jXIr7//LUD3a4dIq8ooqyEj07v75JsEBNtLi+2qiV26liV55tbWOhOjN4+dbykiMZvvhMo5V6Bz1yj000Suv9lVOMX17BTBqQJjVobjNZMe3k7V6Vq/cRBO98lrGGNbmFjE5MQo/P++Zlrg94wZH0DvIX8fTK7fRRK+81t7iSg6VVvtUtw2Av82P9IRIPaNXbqOJXnmtpjNeX7kQ6ywjMYq8ogoOl1ZZHYryAZrolddak1PMwIhgEqK9Z9lAV2k/vXInTfTKKzU2GtbmFTM5Mcqrlg101egB4fTpFaDDLJVbaKJXXmnHkXKOVdQyxcf655v4+QmTh0WxNrcYY3R5QdU5muiVV2rqn5/sg/3zTTISozhYUsX+Y5VWh6K8nEuJXkRmichOEckRkfta2f8LEdkmIltE5BMRGeq0r0FEshyPZS3bKnUm1uQWkxAdyqA+IVaH0mWaJmnT7hvVWe0mehGxAYuAS4Bk4FoRSW5RLRNIM8aMA94E/uC0r8oYk+J4zEGpTqpvaGTdnmM+fTYPkBgTSr+wIE30qtNcOaNPB3KMMXnGmFpgCTDXuYIx5jNjTNP3y6+BWPeGqdR3thws5URNvc/2zzcRETISo1ibW6T99KpTXEn0g4EDTs/zHWVtuQn4wOl5sIhsEJGvReTy1hqIyC2OOhsKCwtdCEn1ZE1DDicNi7Q4kq6XkRhN0YladhecsDoU5cXcejFWRK4H0oA/OhUPNcakAdcBz4pIYst2xpgXjTFpxpi0mJgYd4akfNDa3GJG9g8jqneQ1aF0uabuqTU5Oh2COnOuJPqDQJzT81hH2UlE5ALgAWCOMaamqdwYc9Dxbx6wCkjtRLyqh6upb2DDPt/vn28SF9mLuMgQ7adXneJKol8PJIlIgogEAvOAk0bPiEgq8FfsSb7AqbyviAQ5tqOBKcA2dwWvep7NB0qprmvsMYkeIGNYNN/sOUZDo/bTqzPTbqI3xtQDtwErge3AUmNMtog8KiJNo2j+CPQG3mgxjHI0sEFENgOfAQuNMZro1Rlbk1uECExK6EGJfngUpVV1bD9cZnUoykv5u1LJGLMcWN6i7CGn7QvaaLcGOKszASrlbG1uMWMGhRPRK8DqULrN5GGOfvrcIsYO9v51cVX30ztjldeormsgc3+Jz01L3J5+4cEkxoRqP706Y5roldfYuO84tQ2NzWe4PUlGYjTr9hyjrqHR6lCUF9JEr7zGmtwibH7ChATfHz/fUkZiFJW1DWzJL7E6FOWFNNErr7E2t5hxsfZl9nqaSU399DnafaM6ThO98gonaurZnF/aI7ttAPqGBpI8MFz76dUZ0USvvML6vfZx5D3tQqyzjMQoNu4/TnVdg9WhKC+jiV55hbW5xQTYhPFD+1odimUyhkdRW9/Ipn3HrQ5FeRlN9MorrM0tJnVIX0ICbVaHYpkJ8ZHY/ES7b1SHaaJXHq+0so6th3pu/3yTsOAAxsVGNK+upZSrNNErj/fNnmKMsfdR93STh0WxJd8+H79SrtJErzze2rxigvz9SBnSx+pQLJeRGE19o2H93mNWh6K8iCZ65fHW5hYzIT6SIP+e2z/fZPzQvgTa/JoXX1HKFZrolUcrPlHDjiPlPWpa4tMJCbSROqSP9tOrDtFErzza13n2LgpN9N/JSIwm+1AZJZW1VoeivIQmeuXR1uYVERpo4yydnrdZxvAojPnuQ1Cp9miiVx5tTW4xExIiCbDpr2qTs2P7EBJgY6123ygXufTXIyKzRGSniOSIyH2t7P+FiGwTkS0i8omIDHXad6OI7HY8bnRn8Mq3HS2rJq+wQodVthDo78eEhEi9cUq5rN1ELyI2YBFwCZAMXCsiyS2qZQJpxphxwJvAHxxtI4HfAhOBdOC3ItJz72FXHdI0smTysJ47v01bMhKj2F1wgoLyaqtDUV7AlTP6dCDHGJNnjKkFlgBznSsYYz4zxlQ6nn4NxDq2LwY+MsYcM8YcBz4CZrkndOXrVucUERESQPKgcKtD8ThN33J0mKVyhSuJfjBwwOl5vqOsLTcBH3SkrYjcIiIbRGRDYWGhCyEpX2eM4aucIqYMj8LmJ1aH43HGDIogLNhfE71yiVuvcInI9UAa8MeOtDPGvGiMSTPGpMXExLgzJOWl8ooqOFxazZTh2m3TGpufMDEhSvvplUtcSfQHgTin57GOspOIyAXAA8AcY0xNR9oq1dLq3fYRJdOG6wd/WzISo9h/rJL845XtV1Y9miuJfj2QJCIJIhIIzAOWOVcQkVTgr9iTfIHTrpXARSLS13ER9iJHmVKntTqniLjIEIZE9bI6FI+VMVz76ZVr2k30xph64DbsCXo7sNQYky0ij4rIHEe1PwK9gTdEJEtEljnaHgMew/5hsR541FGmVJvqGxr5OreYqXo2f1oj+oURFRqoiV61y6VVlo0xy4HlLcoectq+4DRtXwZePtMAVc+zOb+U8pp6pmr//Gn5+QmTEqP4KrcIYwwietFatU5vN1Qe56ucIkR0/nlXTB0ezdGyGnYXnLA6FOXBNNErj7N6dxFjB0XQNzTQ6lA83rQk+7eeL3bpsGTVNk30yqNU1NSzaf9xHVbpoti+vRgWE8qXu3XeG9U2TfTKo6zbc4z6RtN8pqraNz0phm/2FFNd12B1KMpDaaJXHuXL3UUE+fsxfqhOieSq6SOiqa5rZMPe41aHojyUJnrlUb7KKWJCfCTBAbpsoKsmJkQRYBO+3K399Kp1muiVxygoq2bn0XKmardNh4QG+TN+aF8+1wuyqg2a6JXHWJ1jv6Co4+c7bvqIGHYcKaegTKctVqfSRK88xue7ConuHUjyQJ2WuKOmJ9nvIm76sFTKmSZ65REaGg1f7CpkelIMfjotcYclDwwnKjRQx9OrVmmiVx5hS34JxyvrOHekzm9zJvz8hKlJ0azOKaKx0VgdjvIwmuiVR/h8VyF+8l0XhOq4aUkxFJ2oZfuRMqtDUR5GE73yCKt2FnJ2XB+d9qATvpsOQfvp1ck00SvLHauoZXN+CeeO0LP5zugfHsyoAWHaT69OoYleWe7L3YUYA+eN7Gd1KF5v+ogYNuw7xomaeqtDUR7EpUQvIrNEZKeI5IjIfa3sny4im0SkXkSuarGvwbEYSfOCJEo5+3xnIZGhgYwbHGF1KF7v/FH9qGswrNa7ZJWTdhO9iNiARcAlQDJwrYgkt6i2H5gPvN7KIaqMMSmOx5xW9qserLHR8PmuQqYnReuwSjcYP7QvYcH+fLK9oP3KqsdwZYWpdCDHGJMHICJLgLnAtqYKxpi9jn2NXRCj8mFbD5VSXFGrwyrdJMDmx7kjYvhsZyGNjUY/PBXgWtfNYOCA0/N8R5mrgkVkg4h8LSKXdyQ45fs+31mI6LBKt5o5uh9FJ2r49mCp1aEoD9EdF2OHGmPSgOuAZ0UksWUFEbnF8WGwobBQ+xZ7kk93FjBucARRvYOsDsVnnDuiHyLw6Q7tvlF2riT6g0Cc0/NYR5lLjDEHHf/mAauA1FbqvGiMSTPGpMXE6JldT1FQXk3WgRIuGN3f6lB8SmRoIOcM6auJXjVzJdGvB5JEJEFEAoF5gEujZ0Skr4gEObajgSk49e2rnu2zHQUYAxcka6J3t/NH9ePbg6U6m6UCXEj0xph64DZgJbAdWGqMyRaRR0VkDoCITBCRfOBq4K8iku1oPhrYICKbgc+AhcYYTfQKgI+2FTC4TwijBoRZHYrPOX+U/Z6Ez3bqWb1ybdQNxpjlwPIWZQ85ba/H3qXTst0a4KxOxqh8UFVtA6tzCpk3YQgiOjLE3UYNCGNQRDCfbC/gBxOGWB2OspjeGasssTqniOq6Ru2f7yIiwoxR/VidU0RNvS4a3tNpoleW+HjbUcKC/ElPiLQ6FJ81c3Q/KmsbWJtbbHUoymKa6FW3a2w0fLLjKOeN6kegv/4KdpWMxGhCA22szD5idSjKYvpXprpd5oESik7UcsFoncSsKwUH2Jgxqh8fZh+lQRcj6dE00atu9/H2o/j7CeeN0ETf1WaNHUBxRS0b9h6zOhRlIU30qlsZY1iZfYT0hEgiegVYHY7PmzHS3j22QrtvejRN9Kpb7Tp6grzCCi45a6DVofQIoUH+TE+KYeXWIxij3Tc9lSZ61a3e//YwInDxGB1W2V1mjR3AodJqneSsB9NEr7rVB98eJj0+kn5hwVaH0mNcMLofNj9hxVbtvumpNNGrbrP7aDm7C05w6TjttulOfXoFMnlYFCu0+6bH0kSvuk1Tt82sMQOsDqXHuXjsAPKKKth19ITVoSgLaKJX3eaDb48wYWgk/cK126a7XTymP34C7205ZHUoygKa6FW3yCk4wc6j5XzvLD2bt0K/sGAyEqNZtvmQdt/0QJroVbd4f8thAGaN1f55q8xJGcS+4ko25+vom55GE73qcsYY3s06yKRhkQyI0G4bq8waO4BAfz/ezXJ5gTjlIzTRqy63Ob+UvKIKrkjtyJryyt3CgwOYMTKG97Yc1rlvehiXEr2IzBKRnSKSIyL3tbJ/uohsEpF6Ebmqxb4bRWS343GjuwJX3uPfmQcJ9PfTbhsPMDdlMIXlNXydp1MX9yTtJnoRsQGLgEuAZOBaEUluUW0/MB94vUXbSOC3wEQgHfitiPTtfNjKW9Q1NPKfzYe4cHR/IkJ0bhurnT+qH72D/LX7podx5Yw+HcgxxuQZY2qBJcBc5wrGmL3GmC1AY4u2FwMfGWOOGWOOAx8Bs9wQt/ISX+4upLiilsu128YjBAfYuGhMfz7YeoTqOl15qqdwJdEPBg44Pc93lLnCpbYicouIbBCRDYWFhS4eWnmDdzIP0bdXAOeOiLE6FOVw5TmxlFfX64IkPYhHXIw1xrxojEkzxqTFxGhC8BXl1XV8mH2E2eMG6UpSHmTysChi+4bwxoZ8q0NR3cSVv76DQJzT81hHmSs601Z5ufe2HKamvpErztFuG0/i5ydcPT6O1TlFHDhWaXU4qhu4kujXA0kikiAigcA8YJmLx18JXCQifR0XYS9ylKkeYMm6/YzsH0ZqXB+rQ1EtXJUWiwi8uVHP6nuCdhO9MaYeuA17gt4OLDXGZIvIoyIyB0BEJohIPnA18FcRyXa0PQY8hv3DYj3wqKNM+bjsQ6Vszi9lXnocImJ1OKqFwX1CmDo8mjc35uuY+h7A35VKxpjlwPIWZQ85ba/H3i3TWtuXgZc7EaPyQkvWHSDQ309vkvJg16TF8fN/ZrImt4hpSXptzJfpFTLldlW1Dfw78yCXnjWQPr0CrQ5HteGiMf3p0yuAJesPtF9ZeTVN9Mrt3ttyiPKaeuZNiGu/srJMkL+NK8+JZeXWIxSUVVsdjupCmuiVWxljePWb/QyLCSU9IdLqcFQ7/mvSUBqM4bVv9lsdiupCmuiVW23aX8LmAyXMz4jXi7BeID46lPNGxPD6uv3U1re8sV35Ck30yq3+76s9hAX7c+U5rV6bVx7oxox4Cstr+GDrYatDUV1EE71ym0MlVXyw9QjzJsQRGuTSgC7lAaYnxZAQHcora/ZaHYrqIproldv84+t9GGO4YXK81aGoDvDzE/5r0lAy95eQdaDE6nBUF9BEr9yisraef67bz8VjBhAX2cvqcFQHXTMhjvBgf/5nVY7VoaguoIleucU/1x2gpLKOn0xLsDoUdQZ6B/lzY0Y8K7OPklNQbnU4ys000atOq65r4MUvcpk0LJLxQ3VIpbeanxFPcIAfiz/PszoU5Waa6FWnvbkxn6NlNfz8/CSrQ1GdENU7iHkThvDvzIMcLKmyOhzlRproVafUNTTyP6tySR3Sh4zEKKvDUZ108/RhAPy/L/Ss3pdooled0nT29/Pzh+sNUj5gcJ8Qrhofy+vf7Cf/uM5V7ys00aszVlPfwLMf72bs4HBmjOxndTjKTW6fmQQCf/54t9WhKDfRRK/O2Ktf7+dgSRX3zRqtZ/M+ZFCfEP5r0lDe2pRPTsEJq8NRbqCJXp2Rsuo6Xvh0N1OHRzM1KdrqcJSb/ey8REICbDz90U6rQ1Fu4FKiF5FZIrJTRHJE5L5W9geJyL8c+78RkXhHebyIVIlIluOx2M3xK4u8+Hkexyvr+NWsUVaHorpAVO8gbpo2jOXfHmHjPl0Uztu1m+hFxAYsAi4BkoFrRSS5RbWbgOPGmOHAM8CTTvtyjTEpjsetbopbWWhfcQUvfpnH3JRBnBUbYXU4qov89/RhDAgP5qF3s3W5QS/nyhl9OpBjjMkzxtQCS4C5LerMBf7m2H4TmCnaaeuzHvnPNgL8hF9/b7TVoaguFBrkzwOXjib7UBn/XKfz1XszVxL9YMB5rbF8R1mrdRyLiZcCTYOqE0QkU0Q+F5Fprb2AiNwiIhtEZENhYWGH3oDqXh9vO8qnOwq444Ik+ocHWx2O6mKzxw1k0rBInvpwJ8craq0OR52hrr4YexgYYoxJBX4BvC4i4S0rGWNeNMakGWPSYmJ0kWJPVVFTz8P/yWZ4v978aIrOadMTiAiPzBlLeXU9TyzfbnU46gy5kugPAs6Lf8Y6ylqtIyL+QARQbIypMcYUAxhjNgK5wIjOBq2s8fsPtnOwpIqF3z+LAJsO2OopRg4I47+nD+ONjfl8trPA6nDUGXDlr3U9kCQiCSISCMwDlrWoswy40bF9FfCpMcaISIzjYi4iMgxIAvTeai+0JqeIV7/ez4+nJJAWrxOX9TR3XJBEUr/e3P/Wt5RW1VkdjuqgdhO9o8/9NmAlsB1YaozJFpFHRWSOo9pLQJSI5GDvomkagjkd2CIiWdgv0t5qjNGxWl6mtKqOe9/cQkJ0KPdcNNLqcJQFgvxtPHX12RSeqOGRZdlWh6M6yKX13owxy4HlLcoectquBq5upd1bwFudjFFZyBjDL9/czNGyat64dTIhgTarQ1IWOTuuDz8/fzjPfrybSYlRXJMW134j5RG0o1Wd1itr9rIy+yj3XTKK1CF9rQ5HWezn5yeRkRjFQ+9uZecRXaDEW2iiV236Jq+YJ5Zv54LR/bhpqo6yUWDzE56dl0LvoABufXUjJZU65NIbaKJXrdpTVMF/v7qRuMhe/OnqFJ20TDXrFxbMX354DvnHK/npq5uorW+0OiTVDk306hTHK2q56ZX1CPB/8ycQ0SvA6pCUh0lPiGTh98exNq+YB975FmN0igRP5tLFWNVzlFbVcf1L35BfUsWrN01kaFSo1SEpD3Xl+Fj2FVfw3Kc5hAUH8JvZOl21p9JEr5qVV9cx///WsetoOS/ekEZ6go6XV6d314UjKKuu5+Wv9hAU4McvLx6pyd4DaaJXABSW1/CjV9ax/XA5i647R1eMUi4REX57WTK1jrWDT1TX8/CcMdj8NNl7Ek30ir1FFdz4f+soKKvhf29M0ySvOkRE+N3csYQF+/PXz/MoLK/h2XkpBAfoPReeQi/G9nAfbzvKZS+spqyqjtdunqhJXp0RPz/h/ktG85vZyazIPsJVi9dw4JguLu4pNNH3ULX1jTy5Ygc/+fsG4qNCWXbbVM7RG6JUJ900NYH/d0Ma+4srufS5L/kw+4jVISk00fdIW/JLuOz51fzPqlzmTYjjjVsnExfZy+qwlI+4MLk/798+jSFRvbjlHxu5c0kmx3Que0tpH30PUnSihj9/vJvX1+0nuncgL92YxszR/a0OS/mguMhevPXTDP7yWS5/WZXDF7uL+MWFI/jBhDid4toCmuh7gOMVtfzj6328+EUeVXUNXJc+hHsuHklEiN4IpbpOkL+Nuy4cwffOGsiD//6WB/+9lZdX7+GuC0dwydgB+GvC7zaa6H1YTsEJ/rF2L0s35FNV18CFyf351axRDO/X2+rQVA8yckAYS/97Mp9sL+DJFTv4+T8zGdwnhB9Niefq8XF653U30ETvYw6VVPHRtqO8vSmfzfmlBNiEuSmDuXnaMEYOCLM6PNVDiQgXJPdnxqh+fLL9KP+7eg+/e387f1ixk/NH9ePy1EFMHxFDr0BNSV1Bf6perqSylsz9JXyz5xirdhawwzF17OiB4Tx46WjmpAyiX5gu4q08g81PuGjMAC4aM4CtB0t5e9NBlm0+xIrsIwTa/EhPiGT6iGjS4iMZMyicIH8di+8OLiV6EZkF/BmwAf9rjFnYYn8Q8HdgPFAM/MAYs9ex737gJqABuN0Ys9Jt0fcglbX17D9WSW5BBTkFJ8gpPEH2wVLyiioA8PcT0uL7cv8lozh/VD+S+uvZu/JsYwdHMHZwBL/+3qjmE5XPdxXyxPIdAATYhOSB4SQPiiAxJpTEmN4kxvRmUJ9g7d/voHYTvWPN10XAhUA+sF5ElhljtjlVuwk4bowZLiLzgCeBH4hIMvY1ZscAg4CPRWSEMabB3W/EUzU2GmobGu2PevujzrFdU99IZW0DZVV1lFXXUV5dT1lVHeU19RyrqOVoWTVHSqs5UlZNeXV98zFFILZvCCP7h3NVWiypcX0ZFxtBaJB+QVPex9/mx5Th0UwZHs0Dl8LRsmoy9x8n80AJWftL+GDrYUoqv1un1k8guncQ/cODHY8gIkICCA8JIDw4gPAQf8KDA+gd7E+Qvx9B/jb7vwFO2/5+PWpOHlcyQzqQY4zJAxCRJcBcwDnRzwUedmy/Cbwg9p/iXGCJMaYG2ONYUzYdWOue8L9TUlnLVYvXYozBABgw0Dx9qn0bDMb+r9Osqk1tnPc7t6GpTWvHc3q978oNjQbqGhqpb+z49K3BAX5EhAQwIDyYYTGhZCRG0T8imMF9QhjerzfDonvrkn7KZ/UPD2bW2IHMGjuwuexYRS15hSfILTzBweNVHCmr5mhZDfnHK9m0/zhlVXUd/luz+Ql+An4ijm378+Ztx36bY1sEhO8+HJw/J5w/Mpw/QE75KGmnzeiB4Tx/bWqH3ocrXEn0g4EDTs/zgYlt1THG1ItIKRDlKP+6RdvBLV9ARG4BbgEYMmSIq7GfxOYnjGzqrhD7D1FEHP82FzeX4fhPO3nfd2Xf/V81/QefvN/5eDjta9r2Ewiw+RHo73jY7GcRzmUBNj96B/kTFmw/AwkL9icsOIBAf/1aqpSzyNBAIkMjSYtvfUZVYwxVdQ2UVdU7vh3bvyHXOL4519Q1fLdd30B1XSMNjY00Gvu37kZjaGiERtO0bf+3sREaHGUtTw6bt0+Kg1bLT9fG+Ulc3xCXfyYd4RHf9Y0xLwIvAqSlpZ3RCgZhwQEs+uE5bo1LKeUdRIRegf70CvRnQIQOPmjJlVPHg4Dzcu+xjrJW64iIPxCB/aKsK22VUkp1IVcS/XogSUQSRCQQ+8XVZS3qLANudGxfBXxq7N9TlgHzRCRIRBKAJGCde0JXSinlina7bhx97rcBK7EPr3zZGJMtIo8CG4wxy4CXgH84LrYew/5hgKPeUuwXbuuBBT1pxI1SSnkCl/rojTHLgeUtyh5y2q4Grm6j7ePA452IUSmlVCfo8A6llPJxmuiVUsrHaaJXSikfp4leKaV8nDjfreUJRKQQ2Gd1HGcgGiiyOohupu+5Z9D37B2GGmNiWtvhcYneW4nIBmNMmtVxdCd9zz2Dvmfvp103Sinl4zTRK6WUj9NE7z4vWh2ABfQ99wz6nr2c9tErpZSP0zN6pZTycZrolVLKx2mi7wIicreIGBGJtjqWriYifxSRHSKyRUTeEZE+VsfUFURklojsFJEcEbnP6ni6mojEichnIrJNRLJF5A6rY+ouImITkUwRec/qWNxFE72biUgccBGw3+pYuslHwFhjzDhgF3C/xfG4nYjYgEXAJUAycK1j4XtfVg/cbYxJBiYBC3rAe25yB7Dd6iDcSRO9+z0D/JJTl4z0ScaYD40x9Y6nX2NfRczXpAM5xpg8Y0wtsAT7wvc+yxhz2BizybFdjj3xnbLes68RkVjgUuB/rY7FnTTRu5GIzAUOGmM2Wx2LRX4MfGB1EF1gMHDA6Xmri9z7KhGJB1KBbywOpTs8i/1ErdHiONzKIxYH9yYi8jEwoJVdDwC/xt5t41NO956NMe866jyA/ev+a90Zm+paItIbeAu40xhTZnU8XUlEZgMFxpiNInKexeG4lSb6DjLGXNBauYicBSQAm0UE7F0Ym0Qk3RhzpBtDdLu23nMTEZkPzAZmGt+8MaNHLnIvIgHYk/xrxpi3rY6nG0wB5ojI94BgIFxEXjXGXG9xXJ2mN0x1ERHZC6QZY7xtBrwOEZFZwNPAucaYQqvj6Qoi4o/9QvNM7Al+PXCdMSbb0sC6kNjPVv4GHDPG3GlxON3OcUZ/jzFmtsWhuIX20avOegEIAz4SkSwRWWx1QO7muNh8G7AS+0XJpb6c5B2mAP8FnO/4f81ynOkqL6Rn9Eop5eP0jF4ppXycJnqllPJxmuiVUsrHaaJXSikfp4leKaV8nCZ6pZTycZrolVLKx/1/UR4x0Vk+apEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np \n",
    "x = np.random.normal(0,2,100)\n",
    "e = np.random.normal(0,4,100)\n",
    "y = 0 + 4*x + e\n",
    "\n",
    "fig,ax = plt.subplots()\n",
    "sns.scatterplot(x,y,ax=ax)\n",
    "ax.set(xlabel=\"X\",ylabel=\"Y\")\n",
    "\n",
    "ax.axhline(0)\n",
    "ax.axvline(0)\n",
    "\n",
    "\n",
    "avgy = np.mean(y)\n",
    "avgx = np.mean(x)\n",
    "N = len(y)\n",
    "\n",
    "mlebeta1 = sum((y-avgy)*(x-avgx))/sum((x-avgx)**2)\n",
    "print(\"MLE for Beta 1 is {:.2f}\".format(mlebeta1))\n",
    "\n",
    "mleBeta0 = avgy - mlebeta1*avgx\n",
    "print(\"MLE for Beta 0 is {:.2f}\".format(mleBeta0))\n",
    "\n",
    "mlesigma2 = np.mean( (y - (mleBeta0+mlebeta1*x))**2 )\n",
    "print(\"MLE for sigma2 is {:.2f}\".format(mlesigma2))\n",
    "\n",
    "mleV = mlesigma2*((1/N) + (avgx**2)/sum((x-avgx)**2))\n",
    "print(\"MLE for V(b0) is {:.2f}\".format(mleV))\n",
    "\n",
    "print(\"MLE for square root of V(b0) is {:.2f}\".format(mleV**0.5))\n",
    "\n",
    "W = mleBeta0/np.sqrt(mleV)\n",
    "print(\"Test statistic (W) is {:.2f}\".format(W))\n",
    "\n",
    "import scipy.stats\n",
    "\n",
    "dom = np.linspace(-5,5,10**3)\n",
    "standardnormpdf = scipy.stats.norm(0,1).pdf\n",
    "\n",
    "fig,ax = plt.subplots()\n",
    "ax.plot(dom,standardnormpdf(dom),label = \"Null Distribution N(0,1)\")\n",
    "\n",
    "ax.axvline(W,color=\"red\",label=\"Our Test stat W\")\n",
    "\n",
    "ax.legend()"
   ]
  }
 ],
 "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}