Chapter 08 - Estimating the CDF and Statistical Functionals.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Estimating the CDF and Statistical Functionals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.1 Empirical distribution function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **empirical distribution function** $\\hat{F_n}$ is the CDF that puts mass $1/n$ at each data point $X_i$.  Formally,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{F_n}(x) & = \\frac{\\sum_{i=1}^n I\\left(X_i \\leq x \\right)}{n} \\\\\n",
    "&= \\frac{\\text{#}|\\text{observations less than or equal to x}|}{n}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "I\\left(X_i \\leq x\\right) =\n",
    "    \\begin{cases}\n",
    "      1   & \\text{if } X_i \\leq x \\\\\n",
    "      0   & \\text{if } X_i > x\n",
    "    \\end{cases}       \n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 8.3**.  At any fixed value of $x$,\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbb{E}\\left( \\hat{F_n}(x) \\right) = F(x)\n",
    "\\quad\\mathrm{and}\\quad \n",
    "\\mathbb{V}\\left( \\hat{F_n}(x) \\right) = \\frac{F(x)(1 - F(x))}{n}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Thus,\n",
    "\n",
    "$$ \\text{MSE} = \\frac{F(x)(1 - F(x))}{n} \\rightarrow 0 $$\n",
    "\n",
    "and hence, $\\hat{F_n}(x) \\xrightarrow{\\text{P}} F(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 8.4 (Glivenko-Cantelli Theorem)**.  Let $X_1, \\dots, X_n \\sim F$.  Then\n",
    "\n",
    "$$ \\sup _x |\\hat{F_n}(x) - F(x)| \\xrightarrow{\\text{P}} 0 $$\n",
    "\n",
    "(actually, $\\sup _x |\\hat{F_n}(x) - F(x)|$ converges to 0 almost surely.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.2 Statistical Functionals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A **statistical functional** $T(F)$ is any function of $F$.  Examples are the mean $\\mu = \\int x dF(x)$, the variance $\\sigma^2 = \\int (x - \\mu)^2 dF(x)$ and the median $m = F^{-1}(1/2)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **plug-in estimator** of $\\theta = T(F)$ is defined by\n",
    "\n",
    "$$\\hat{\\theta_n} = T(\\hat{F_n}) $$\n",
    "\n",
    "In other words, just plug in $\\hat{F_n}$ for the unknown $F$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A functional of the form $\\int r(x) dF(x)$ is called a **linear functional**.  Recall that $\\int r(x) dF(x)$ is defined to be $\\int r(x) f(x) d(x)$ in the continuous case and $\\sum_j r(x_j) f(x_j)$ in the discrete."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plug-in estimator for the linear functional $T(F) = \\int r(x) dF(x)$ is:\n",
    "\n",
    "$$T(\\hat{F_n}) = \\int r(x) d\\hat{F_n}(x) = \\frac{1}{n} \\sum_{i=1}^n r(X_i)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have:\n",
    "\n",
    "$$ T(\\hat{F_n}) \\approx N\\left(T(F), \\hat{\\text{se}}\\right) $$\n",
    "\n",
    "An approximate $1 - \\alpha$ confidence interval for $T(F)$ is then\n",
    "\n",
    "$$ T(\\hat{F_n}) \\pm z_{\\alpha/2} \\hat{\\text{se}} $$\n",
    "\n",
    "We call this the **Normal-based interval**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3 Technical Appendix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 8.12 (Dvoretsky-Kiefer-Wolfowitz (DKW) inequality)**.  Let $X_1, \\dots, X_n$ be iid from $F$.  Then, for any $\\epsilon > 0$,\n",
    "\n",
    "$$\\mathbb{P}\\left( \\sup_x |F(x) - \\hat{F_n}(x) | > \\epsilon \\right) \\leq 2 e^{-2n\\epsilon^2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the DKW inequality, we can construct a confidence set.  Let $\\epsilon_n^2 = \\log(2/\\alpha) / (2n)$, $L(x) = \\max \\{ \\hat{F_n}(x) - \\epsilon_n, \\; 0 \\}$ and $U(x) = \\min \\{\\hat{F_n}(x) + \\epsilon_n, 1 \\}$.  It follows that for any $F$,\n",
    "\n",
    "$$ \\mathbb{P}(F \\in C_n) \\geq 1 - \\alpha $$\n",
    "\n",
    "To summarize:\n",
    "\n",
    "A $1 - \\alpha$ nonparametric confidence band for $F$ is $(L(x), \\; U(x))$ where \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "L(x) &= \\max \\{ \\hat{F_n}(x) - \\epsilon_n, \\; 0 \\} \\\\\n",
    "U(x) &= \\min \\{ \\hat{F_n}(x) + \\epsilon_n, \\; 1 \\} \\\\\n",
    "\\epsilon_n &= \\sqrt{\\frac{1}{2n} \\log \\left( \\frac{2}{\\alpha} \\right) }\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.5 Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.1**.  Prove Theorem 8.3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.  We have:\n",
    "\n",
    "$$\n",
    "\\hat{F_n}(x) = \\frac{\\sum_{i=1}^n I\\left(X_i \\leq x \\right)}{n}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "I\\left(X_i \\leq x\\right) =\n",
    "    \\begin{cases}\n",
    "      1   & \\text{if } X_i \\leq x \\\\\n",
    "      0   & \\text{if } X_i > x\n",
    "    \\end{cases}       \n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Thus,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{E}(\\hat{F_n}(x)) & = n^{-1} \\sum_{i = 1}^n \\mathbb{E}(I\\left(X_i \\leq x \\right)) \\\\\n",
    " & = n^{-1} \\sum_{i = 1}^n \\mathbb{P}\\left(X_i \\leq x \\right) \\\\\n",
    " & = n^{-1} \\sum_{i = 1}^n F(x) \\\\\n",
    " & = F(x)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{E}(\\hat{F_n}(x)^2) & = n^{-2} \\mathbb{E} \\left( \\sum_{i = 1}^n I\\left(X_i \\leq x \\right) \\right)^2 \\\\\n",
    "& = n^{-2} \\mathbb{E} \\left( \\sum_{i = 1}^n I\\left(X_i \\leq x \\right)^2 \n",
    "+ \\sum_{i = 1}^n \\sum_{j = 1, j \\neq i}^n I\\left(X_i \\leq x \\right) I\\left(X_j \\leq x \\right) \\right) \\\\\n",
    "& = n^{-2} \\left( \\sum_{i = 1}^n \\mathbb{E} \\left( I\\left(X_i \\leq x \\right)^2 \\right)\n",
    "+ \\sum_{i = 1}^n \\sum_{j = 1, j \\neq i}^n \\mathbb{E} \\left( I\\left(X_i \\leq x \\right) I\\left(X_j \\leq x \\right) \\right) \\right) \\\\\n",
    "& = n^{-2} \\left( \\sum_{i = 1}^n \\mathbb{E} \\left( I\\left(X_i \\leq x \\right) \\right)\n",
    "+ \\sum_{i = 1}^n \\sum_{j = 1, j \\neq i}^n \\mathbb{E} \\left( I\\left(X_i \\leq x \\right) \\right) \\mathbb{E} \\left( I\\left(X_j \\leq x \\right) \\right) \\right) \\\\\n",
    "& = n^{-2} \\left( \\sum_{i = 1}^n \\mathbb{P}\\left(X_i \\leq x \\right) \n",
    "+ \\sum_{i = 1}^n \\sum_{j = 1, j \\neq i}^n \\mathbb{P}\\left(X_i \\leq x \\right) \\mathbb{P}\\left(X_j \\leq x \\right)  \\right) \\\\\n",
    "&= n^{-2} \\left( \\sum_{i = 1}^n F(x)\n",
    "+ \\sum_{i = 1}^n \\sum_{j = 1, j \\neq i}^n F(x)^2  \\right) \\\\\n",
    "&= n^{-2} \\left( nF(x) + (n^2 - n)F(x)^2 \\right) \\\\\n",
    "&= n^{-1} ( F(x) + (n - 1) F(x)^2 )\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "\n",
    "Therefore,\n",
    "\n",
    "$$ \\mathbb{V}(\\hat{F_n}(x)) \n",
    "= \\mathbb{E}(\\hat{F_n}(x)^2) - \\mathbb{E}(\\hat{F_n}(x))^2\n",
    "= F(x) /n + (1 - 1/n)F(x)^2 - F(x)^2\n",
    "= \\frac{F(x)(1 - F(x))}{n}\n",
    "$$\n",
    "\n",
    "Finally,\n",
    "\n",
    "$$ \\text{MSE} = (\\text{bias}(F(x)))^2 + \\mathbb{V}(F(x)) = (\\mathbb{E}(\\hat{F_n}(x)) - F(x))^2 + \\mathbb{V}(F(x)) = \\mathbb{V}(F(x)) = \\frac{F(x)(1 - F(x))}{n} \\rightarrow 0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.2**.  Let $X_1, \\dots, X_n \\sim \\text{Bernoulli}(p)$ and let $Y_1, \\dots, Y_m \\sim \\text{Bernoulli}(q)$.\n",
    "\n",
    "- Find the plug-in estimator and estimated standard error for $p$.\n",
    "- Find an approximate 90-percent confidence interval for $p$.\n",
    "- Find the plug-in estimator and estimated standard error for $p - q$.\n",
    "- Find an approximate 90-percent confidence interval for $p - q$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "**(a)**\n",
    "\n",
    "$p$ is the mean of $\\text{Bernoulli}(p)$, so its plugin estimator is $\\hat{p} = \\mathbb{E}(\\hat{F_n}) = n^{-1} \\sum_{i=1}^n X_i  = \\overline{X}_n$.\n",
    "\n",
    "$\\sqrt{p(1-p)}$ is the standard error of $\\text{Bernoulli}(p)$, so its plugin estimator is $\\sqrt{\\hat{p}(1-\\hat{p})} = \\sqrt{\\overline{X}_n(1 - \\overline{X}_n)}$.\n",
    "\n",
    "**(b)**\n",
    "\n",
    "The 90-percent confidence interval for $p$ is $\\hat{p} \\pm z_{5\\%} \\hat{se}(\\hat{p})= \\overline{X}_n \\pm z_{5\\%} \\sqrt{\\overline{X}_n(1-\\overline{X}_n)}$.\n",
    "\n",
    "**(c)**\n",
    "\n",
    "The plug-in estimator for $\\theta = p - q$ is $\\hat{\\theta} = \\hat{p} - \\hat{q} = \\overline{X}_n - \\overline{Y}_m$.\n",
    "\n",
    "The standard error of $\\hat{\\theta}$ is\n",
    "\n",
    "$$\\text{se} = \\sqrt{\\mathbb{V}(\\hat{p} - \\hat{q})} = \\sqrt{\\mathbb{V}(\\hat{p}) + \\mathbb{V}(\\hat{q})} = \\sqrt{\\hat{p}(1 - \\hat{p}) + \\hat{q}(1 - \\hat{q})} = \\sqrt{\\overline{X}_n(1 - \\overline{X}_n) + \\overline{Y}_m(1 - \\overline{Y}_m)} $$\n",
    "\n",
    "**(d)**\n",
    "\n",
    "The 90-percent confidence interval for $\\theta = p - q$ is $\\hat{\\theta} \\pm z_{5\\%} \\hat{se}(\\hat{\\theta}) = \\overline{X}_n - \\overline{Y}_m \\pm z_{5\\%} \\sqrt{\\overline{X}_n(1 - \\overline{X}_n) + \\overline{Y}_m(1 - \\overline{Y}_m)} $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.3**. (Computer Experiment)  Generate 100 observations from a $N(0, 1)$ distribution.  Compute a 95 percent confidence band for the CDF $F$.  Repeat this 1000 times and see how often the confidence band contains the true function.  Repeat using data from a Cauchy distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import norm, cauchy\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from tqdm import tqdm_notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x26bbd31ea90>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3RU1drH8e9OSAIhoSXUQOgQegugFJEqPRQhSBEQ5XIFlBe7WFG89sJVr4oUQTpSpEgVpEuQXkMPIabSQkLKZPb7x4FICWSAmTnJ5PmslTWZOSdznkHyc7PPLkprjRBCiNzPzewChBBC2IcEuhBCuAgJdCGEcBES6EII4SIk0IUQwkXkM+vC/v7+ukKFCmZdXgghcqW//vorXmtdPKtjpgV6hQoV2Llzp1mXF0KIXEkpdeZOx6TLRQghXIQEuhBCuAgJdCGEcBGm9aFnJT09ncjISFJSUswuxSHy589P2bJl8fDwMLsUIYQLylGBHhkZia+vLxUqVEApZXY5dqW1JiEhgcjISCpWrGh2OUIIF5Rtl4tSaopSKlYpdeAOx5VSaqJS6rhSap9SquH9FpOSkoKfn5/LhTmAUgo/Pz+X/deHEMJ8tvShTwM63uV4J6Dqta/hwP8epCBXDPPrXPmzCSHMl22Xi9Z6o1Kqwl1OCQGma2Md3u1KqSJKqdJa67/tVKMQIgeIuRLDnANzSLiaYHYpuV6LwBZ0qNzB7u9rjz70AODsDc8jr712W6ArpYZjtOIJDAy0w6WFEI5ksVpYfWI10/ZMY/GRxaRb01HIvzQfiIZXmr6SYwM9q/+6We6aobX+AfgBIDg4OEfurOHj48OVK1fMLkMIU+2N3stPe39i1v5ZxCTFUKxAMUY2Hsm/gv9FkH+Q2eXlKFpr0hLTSIpNIikuiaTYJJLjko3nt3yfFGc8b/5qc6Oz2s7sEeiRQLkbnpcFouzwvkIIJwtPCKfv/L7sjdmLh5sHXap1YXC9wXSu2hlPd0+zyzONNcPKpTOXSDiWwPlj50kI/+fx8rnLZKRmZPlznr6eFCxekIIlClK4fGHKNC5DwRIFqdCqgkPqtEeg/wqMUkrNAZoCl+zSfz5mDOzZ88Bvc5P69eHLL+3yVkOGDKFQoULs3LmT6OhoPv74Yx5//HG7vLcQZvli2xccO3+Mrzt9TWjtUPy9/c0uyWm01iRGJd4U1tcfL5y8QEbaP6Ht6eNJsarFKNO4DEG9gvAp6YN3cW8KljDCu2DxgngX98ajgHPnnGQb6Eqp2cCjgL9SKhJ4G/AA0Fp/B6wAOgPHgWRgqKOKzWn+/vtvNm/ezJEjR+jevbsEusjVLqVc4tfwX+lYpSMjm4w0uxyHsaRaiN4dTfyR+H9C+1rLOz05PfM8dy93ilUphn8Nf6p1r4ZfNT/8qvrhV82PgiUL5shRa7aMcnkim+MasP9/fTu1pB2pR48euLm5UbNmTWJiYswuR4j79nfi33Sa2YnYpFieDX7W7HLsKvVyKme3nSViUwQRmyKI/DMys4vELZ8bRSoWwa+aHxVaV8gM7GJVi1G4XGGUW84L7bvJUTNFcxsvL6/M743/rwmR+1xJu0LzKc2JTYplef/ltK3U1uySHkhSbBJnNp3JDPDoPdFoq0a5K0o3KE3jZxsT2CKQEnVKUKRCEdw93M0u2W4k0IXI4w7EHuDUxVPM7DXTIUPpHElrzcXTF4nYFMGZjUaIJ4Qb4+Tz5c9H2YfK0nJcSwJbBlLu4XJ4+rj2jV0J9FskJydTtmzZzOdjx45l7NixJlYkhGMlpSUBULZQ2WzOzBmsGVbObjnL4UWHObr4KBdPXwQgf5H8BLYIpMGwBgS2DKRMozK4e7pO69sWEui3sFqtNp03bdq0m57L2HWRW+2JNkaTVSpayeRK7sySauHUulNGiC85SnJcMu5e7lRuX5lmLzUjsGUgJWqVyHV93vYmgS5EHnY0/ig/7f2JmsVr5rgWujXDyonVJ9g3Yx/hy8JJS0zD09eTal2qEdQziCqdquDl65X9G+UhEujZmDBhAvPnz7/ptT59+jBu3DiTKhLiwZ27fI53/3iXKbunUMCjAJO6TTK7pEwJ4QnsmbaHvdP3kngukQLFClArtBY1etagYtuK5POS2LoT+ZPJxrhx4yS8hUv55dAvDFw0kAxrBiMbj2TcI+MoUbCEqTWlJqZyaP4hdk/ZzdktZ1FuiiqdqtDxq45U61pNQtxG8qckRB6SYklhzKoxBPkHsbDvQioWNXezlbjDcez4egf7pu8j7UoaftX9aPdRO+oOqotvaV9Ta8uNJNCFyEOm7p5K5OVIpoZMNS3MtdacXn+abZ9t49iKY7h7uVO7X20a/asRZR8qmyNnYOYWEuhC5CEz98+kbsm6tK3o/MlD2qo5tOAQm/+zmeg90XgX96bVO61o/O/GFCxR0On1uCIJdCHyiOPnj7P17FbebvW2U1vBWmuOLDrChrc3EHsgFv8gf7pN6kbdgXXJl18iyJ7kT/MWp0+fpmvXrhw48M8Wqu+88w4+Pj68+OKLJlYmxP1LsaQQuiCUwvkL81SDp5xyTa014cvC2fD2BqJ3R+NXzY9eM3tRK7QWbu627H4p7pUEuhAuLtWSytAlQ9n19y6W9FtCucLlsv+hB6C15vjK42x4awNRO6MoWqkoPX7qQZ3+dXDLJ0HuSDk20MesHJM5g81e6peqz5cd7buK46OPPkrTpk1Zv349Fy9eZPLkybRs2dKu1xDifsUlxdFrXi82R2zmP23/Q/fq3R16vcg/I1n1f6uI3BZJ4fKF6T65O3UH1XWpBbByshwb6LmJxWJhx44drFixgnfffZe1a9eaXZIQpGWk0WJqCyIuRTCn9xxCa4c67FqJUYmsfXUt+2bsw6eUD13+14UGTzXIc2upmC3HBrq9W9K2utPNorvdROrVqxcAjRo14vTp044oS4h7tix8GeEJ4Szos4DeNXs75BqWVAvbv9jOpgmbyEjLoPkrzWk5rqVMyTdJjg10s/j5+XHhwoWbXjt//jwVK955zO71ddHd3d2xWCwOrU8IW2it+TbsWwJ8A+gR1MMh1zi24hi/PfcbF05coHr36nT4rAPFqhRzyLWEbeQOxS18fHwoXbo069atA4wwX7lyJS1atDC5MiFst+DQAtadWseLzV7E3c2+3R6Jfycyv+98ZnWZhbuHOwNXDaTfkn4S5jmAtNCzMH36dEaOHMkLL7wAwNtvv03lypVNrkoI22it+b9V/0fD0g0Z1WSUXd939+TdrH5xNZYUC63fb03zl5pLP3kOIoGehZo1a7J+/Xqbzt2wYUPm9/7+/tKHLkyXmpHKucRzjGw8knxu9vkVv3jmIkufWcrJNScp36o83SZ1w6+qn13eW9iPBLoQLmZn1E4ACno++HR6rTX7ZuxjxagVoKHL/7rQaHijPL+RRE4lgW6jkSNHsmXLlptee/755xk6dKhJFQlxuzUn1tBzbk8qFa3E4zUff6D3unrhKstHLOfgvIMEtgyk5/SeFKlQxE6VCkeQQLfRN998Y3YJQtzRwdiDfL7tc6bvm06t4rVYOXAlpXxK3ff7nd5wmkWDFnEl+gptJrSh+SvNZbr+/UpOhvBwOHIEjh41Hrt0gYED7X4pCXQhcimtNWtPruWzbZ+x6sQqCuQrwPCGw5nQdgJF8t9fSzojLYP1b61ny8dbKFalGE9tfYqAxgF2rtwFaQ3nzv0T2Dc+RkT8c55SUKECPPSQQ8qQQBciF9Ja89SvTzFtzzRKFizJ+63fZ0TwCPy87/9G5cUzF1nQdwHndpyj4TMNeezzx/D08bRj1blYWhpERRmhHRlpfF3//sQJowV+40bxvr4QFASPPGI8Vq9uPFapAvnzO6xMCXQhcqEpu6cwbc80Xnz4Rd5v8z5e+R5sZmb48nAWDVqEztD0md+Hmo/XtFOlucCVKzcHdFahHRt7+88VLAhlyxot7pYtbw7uUqWM1riTSaBnITo6mjFjxhAWFoaXlxcVKlTgyy+/pF69egQFBZGSkoKvry8jR45k8ODBAEybNo2XXnqJgADjn6d169Zl+vTpZn4M4aKOJRxj9G+jaVepHR+2+/CBJg5ZM6xseGcDm97fRKn6pegzv4/rTxBKT4dBg2D/fiOsL1++/ZxixYywLlsWGjX65/uAgH++L1TIlNC+Gwn0W2it6dmzJ4MHD2bOnDkA7Nmzh5iYGCpXrszu3bsBOHnyJL169cJqtWaOdAkNDeXrr782rXbh+qzayjNLn8HT3ZPpPaY/UJinXExh4cCFHFt+jPpP1afz153xKOBhx2pzqMhImDsXmjSBwYP/CekbHwsUMLvK+yKBfov169fj4eHBiBEjMl+rX7/+bROGKlWqxOeff84LL7wgQxeF08zeP5s/zvzBpG6TKO1b+r7fJ+5wHHN7zOXCyQt0/rYzwSOCXXcvz6tX4fhxOHbM+NpzbVnu556DAQPMrc3Ocmygrxyzkug90XZ9z1L1S9Hxy453PefAgQM0atTIpvdr2LAhR44cyXw+d+5cNm/eDMgYdWF/GdYM3tv4HnVL1n2gXYeOLj3KwgEL8SjgwZO/P0n5luXtWKVJUlLg5Ml/QvvGr8jIm88tUcK4WemgkSZmyrGBnhtorW96Ll0uwlEup15m9G+jOZpwlLmPz8VN3fuYcK012z7bxpqX11C6YWlCF4VSuFxhB1TrYFobQwLXrIG1a2HfPmNo4I2/j35+ULUqtG5tPF7/qlIFCufCz2wjmwJdKdUR+ApwB37UWn94y/HCwM9A4LX3/FRrPfVBCsuuJe0otWrVYsGCBTadu3v3bmrUqOHgikRet+nMJp5c/CQRlyJ4o+Ub9KnZ557fIyM9gxUjV7Br0i5q9qlJj2k98PDORf3lcXGwbh2sXm0E+fVWd5Uq0Lw5DBlyc3AXLWpquWbJNtCVUu7AN0B7IBIIU0r9qrU+dMNpI4FDWutuSqniwFGl1EytdZpDqnagNm3a8PrrrzNp0iSeeeYZAMLCwkhOTr7pvNOnT/Piiy8yevRoM8oUeUTEpQjaTG9D+cLl2Tx0Mw+Xe/ie3yP1cirzHp/HyTUnaTmuJa3Ht875a7FkZMCOHbBihfG1a5fxepEi0K4dtG9vfN1ln4K8yJYWehPguNb6JIBSag4QAtwY6BrwVcZdFR/gPJArd3pQSrFo0SLGjBnDhx9+SP78+TOHLZ44cYIGDRpkDlscPXq09JMLh9oTvQeL1cLPvX7mobL33ud7+dxlZnWeRdyhOLpP6U6DoQ0cUKWdnD9vtMCXL4eVKyE+HtzcoFkzeO896NDBGELoLsv13oktgR4AnL3heSTQ9JZzvgZ+BaIAXyBUa2299Y2UUsOB4QCBgYH3U69TlClThnnz5t32+tWrV+/4M0OGDGHIkCEOrErkRacunAKgSrEq9/yzsQdjmdlxJimXUui/vD+VO+TQNf1TU6F/f1i8GKxWo/+7UydjvZMOHYwx4cImtgR6Vv8207c8fwzYA7QBKgNrlFKbtNY3jdjXWv8A/AAQHBx863sIIW5x1WI0Inw8fe7p5yI2RzC722zy5c/H0I1DKVX//hfqciitYfhwWLgQXnwReveGxo2lFX6fbLlVHgmUu+F5WYyW+I2GAgu14ThwCgiyT4lC5F1X0q6gUPe0UUX48nBmtJ9BwRIFGbZtWM4N84gIGDoUpk+Hd9+FTz4xhhJKmN83WwI9DKiqlKqolPIE+mF0r9woAmgLoJQqCVQHTt5PQbcOBXQlrvzZhGPsOLeDOiXr2Bzo+2ftZ26PuRSvWZyhm4fmzPXLIyPh2WeNESqzZ8PLL8Obb5pdlUvI9m+J1tqilBoFrMIYtjhFa31QKTXi2vHvgPeAaUqp/RhdNK9orePvtZj8+fOTkJCAn5+fy81a01qTkJBAfgeutCZcy76YfWyK2MQzDZ+x6fyd3+9k+b+XU6FVBfot6YdXoQdbsMturFY4cAB+//2foYdaw7Bh8PrrUK5c9u8hbKLMajUGBwfrnTt33vRaeno6kZGRpKSkmFKTo+XPn5+yZcvi4ZGLxv8KU8QlxdF4UmMsVgthz4RlO81/+1fbWTVmFVW7VKXP/D7mrsmitbGk7PUAX7/eGEcOxhjxTp1g7Fgo7wIzVE2glPpLax2c1bEcNVPUw8ODijKuVORxV9Ku0G12N2KSYtg0dFO2Yb75w82se20dNXrXoPes3rh7mtQHffAg/PwzzJkD19c+CggwArxNG+NLWuMOlaMCXYi8LtWSSu95vQmLCuOXvr8QXCbLhhhgdONteGcDG8dvpE7/OvT4qQdu+Zy8TVxUlNEP/vPPxqJX7u7GUMOXXzYCvFq1HLfErCuTQBciB/l066esPrGaKd2n0COoxx3P01qz7rV1bPloC/Wfqk+3H7o5d8/PtDQYMwa++87oYmnSBCZOhNBQY/ErYQoJdCFyiCtpV/hi+xd0qdqFoQ3uPANZa83aV9ay9ZOtNBrRiC7fdHHuVP74eGO8+MaNMHo0jBpltMSF6STQhcgBtp7dyohlI0i4msC4luPueJ7WmrWvGmEe/Gwwnb/u7NwRYSkpxmJYZ87AzJnGDE+RYzi5w00IcaOE5AT+tfRfNJ/SnIspF1nSb8kdF+DSWvP7G7+z9eOtBP/bhDAH+P57Y0PkhQslzHMgaaELYYLoK9F8vu1zvg37lhRLCmMfGsu7rd+96xT/je9tZPMHm2n4TENzwjw8HN5/31hjvHNn515b2EQCXQgnik+OZ/wf45m0axJpGWk8UfsJXm/5OjWL17zrz239dCsb3t5AvSfr0fW7rs5f/vbAAWPZWqVANnHJsSTQhXCi19a+xtQ9UxlSfwivtnjVplUUw74NY81La6jVtxbdJ3d3XphrDVu3GmutzJ4NPj7GJCHZ1CXHkkAXwokOxB2gVYVW/Nj9R5vO3zt9LytGrqB69+r0/Lmnc8aZnzoFM2YYQX7iBHh7Q8+eMH48VKrk+OuL+yY3RYVwgr3Re+k8szPbI7fTsFRDm37m8KLDLBm6hIptK/L43Mdx93DCDNBly4zQfvttCAyEadMgOtqYOCRhnuNJC10IB4pNimXsqrHM2j+LwvkL81G7j3iu6XPZ/tzJdSf5pd8vBDQJoN/ifuTL76Rf1cmTjen6W7bIWiu5kAS6EA70/sb3mXtwLi83f5lXmr9C0QLZb158Luwcc0Lm4Ffdj/4r+uPp4+mESoHt22HVKmPDZQnzXEkCXQgHikqMomqxqnzY7kObzo8/Es/MTjMpWKIgA1cNpEDRAg6u8Jrly6FPHyhTBl591TnXFHYnfehCOFB4QjilfGzbMejyucvM6DADN3c3Bq0ehG9pXwdXd01kJPToYYxe2bLF6DsXuZIEuhAOcijuEPtj9xNSPSTbc1MuphgbOl9MYcDKARSr4sSNkc+cAYsF/vMfKFnSedcVdieBLoQDrD25lvYz2lMgXwH61up713MtKRbm9JhD/NF4QheFUrrB3dc/t7sLF4xHn3vbiFrkPBLoQtiRVVt5afVLtJ/RnkJehdg6bOtdN6iwZlhZNGgRZ/44Q4+felCprQlDAzdtAg8PqFvX+dcWdiU3RYWwo2Xhy/h026c83eBpvur0Fd4e3nc9f/ULqzm04BDtP21PnSfqOKnKa1JS4NtvjQW3WrSQFroLkBa6EHY0adckSvmU4tsu32Yb5tu/2s6fX/1J0+eb0uyFZk6qEMjIgKlTjTXMX3gBmjY1NqoQuZ4EuhB2sur4KpaFL+OZhs/g4X73TZoPLzrMqv9bRVDPIDp81sFJFV7z2mvw1FNQqpSxifOqVbJBhYuQLhch7CA2KZbBiwdTq3gtXmvx2l3PjdweycL+CynbtCy9fu7l3K3joqKMreIGDjTWapH9Pl2KBLoQdvDxlo+JT45nzaA1FPC482Sg8yfOM7vbbHwDfOn3az88vO/ekre7iRONLpfx4yXMXZB0uQjxgC6nXmbSrkn0qdWHOiXvfGPz6oWrzOoyC23VDFgxgILFCzqxSozlcBcsMNY1r1jRudcWTiGBLsQD+ivqLy6nXmZwvcF3PCcjLYN5vedx4eQFQheH4lfNz4kVXrNkibEcbq9ezr+2cArpchHiAVm1FYCCHlm3uLXWLH92OafXn6bH9B6Ub2nCwlenTsHQodCwIQwa5PzrC6eQFroQDyg+OR4AX6+s117Z+slWdk/ezSNvPkK9QfWcWZrh8GGjm0VrmD8f8ud3fg3CKSTQhXhAYVFheLl7Uat4rduOHV54mLWvrqVWaC0effdR5xe3Zg08/DBcuQIrV8omFS5OAl2IB7Q/dj+1StS6bez537v+ZuHAhQQ0CSBkagjK2aNK4uKge3coVw527ICHHnLu9YXTSaAL8YBOXzxNpaI3t3wT/05kdvfZePt7029JPzwKOHl4IsDMmcb0/tmzZcOKPMKmQFdKdVRKHVVKHVdKZbn6vVLqUaXUHqXUQaXUH/YtU4icKcWSwskLJ6lStErma+lX05kTMoeUiyk8sfQJfEqasEZKcrIx5rxJE6hd2/nXF6bIdpSLUsod+AZoD0QCYUqpX7XWh244pwjwLdBRax2hlCrhqIKFyElWn1iNxWqhcUBjwBjR8utTvxK1M4rQRaGUqmfb5hZ29957xsiWqVPNub4whS0t9CbAca31Sa11GjAHuHXF/v7AQq11BIDWOta+ZQqRs5w4f4J+C/oRMieE0j6laVW+FQCbPtjEgTkHaPtBW4JCgswpzmqFzz+H/v2hVStzahCmsCXQA4CzNzyPvPbajaoBRZVSG5RSfymlnszqjZRSw5VSO5VSO+Pi4u6vYiFM9t4f7xH0TRBLw5fyRss3ODLqCH7efhxZfIT1b6yn7sC6NH+luXkFpqRAWhrUM2GIpDCVLROLsro1r7N4n0ZAW6AAsE0ptV1rHX7TD2n9A/ADQHBw8K3vIUSuMH7jeFqVb8WMnjMyN6+I2ReTOaKl26Ruzh/Rcp3W8MUXxvclpOczr7El0COBcjc8LwtEZXFOvNY6CUhSSm0E6gHhCOFiMqwZNCvXLDPMk+OTmRMyB69CXoQuCiVffpMmYKelwauvGoE+YIDxJfIUW7pcwoCqSqmKSilPoB/w6y3nLAFaKqXyKaW8gabAYfuWKoT5DsUdQqMp7WOEeUZ6BvP7zCfx70T6Le6Hb5msZ4s61PnzxgbPFSoYYf7cc8bSuB4mDJUUpsq2KaG1tiilRgGrAHdgitb6oFJqxLXj32mtDyulVgL7ACvwo9b6gCMLF8IMvxz6BYUiJMgYF7Dq/1ZxesNpes7oSUCTW28tOVhUFHzwgTGSJTkZOnSAadOgfXtZGjePsunfhlrrFcCKW1777pbnnwCf2K80IXIWq7by096feKT8I5TxLcNfk/4i7JswHn7xYeoOdPIGy0lJRnAfO2Z0rYwdC3WcvCepyHFktUUhbLTq+CpOXDjB+NbjidgSwYqRK6jcoTLtPmzn/GJGjzYW3Vq92lh4Swhk6r8QNjlz8QxDlgyhctHKtPNtx7ze8yhSvgi95/R27hZyYKxrPnUqvP66hLm4ibTQhcjGobhDhC4IJdWSytq+a1nSewnpSekM/n0wBYreebs5h7h0CUaOhLp14e23nXttkeNJC12IO7iUcomxq8ZS77t6RF6OZEGfBZx48wRRO6Po+XNPitcs7vyiPvnEuBn6448yikXcRgJdiFtcTr3MF9u+oPrX1fly+5cMrT+U8FHh+CzzYe/0vbR6p5U50/oTEuC//4XevaFxY+dfX+R40uUixDURlyKY+OdEJu2axOXUyzxS/hGW9V9GcJlgTq0/xeoXVhPUI4hWb5qwPsrGjfDkk8boljfecP71Ra4ggS7yvLikOF5e+zIz9s4A4PGajzP24bE0CWgCwMUzF1nQdwH+1f3pMb0Hys2JY7y1hnHj4MMPjd2GNm2SNVrEHUmgizxLa82s/bN4fuXzXE69zOgmoxnz0BjKF/lnM4j05HTm9pxLRnoGoYtD8fL1cm6Re/YYs0AHDoT//Q98TFhbXeQaEugiz3pm6TNM3j2ZpgFNmdx9MrVK3LwnqNaapcOXEr0nmv7L+uNX1c/5RR45Yjy+8oqEuciWBLrIsxYcWkCvGr2Y9/g83N3cbzu+/cvt7J+5n9bvt6Zq56omVAhERhqPgYHmXF/kKjLKReRJaRlppGWkEVgoMMswP/X7Kda8tIYavWrQ8vWWJlR4TXq68ejl5K4ekStJoIs86aXVL3HVcpV2lW6faXnx9EXm952Pf3V/QqaFmLe2OUBMDOTPL2POhU0k0EWeM+/gPCbumMjzTZ+nS7UuNx27fhPUarGacxP0Vps3w0MPgZv8qorsyd8Skafsi9nH0CVDaVauGR+3//imY5k3QfdG03tWb3Nugt5o8mTYtQu6dMn+XCGQm6IiD0nPSKfn3J4UyV+EBX0W4OnuedPxHHET9LpNm+Df/zbWOB8zxtxaRK4hgS7yjF1/7+LkhZPM6jUrc/u4667fBA3qGUTL10y8CQqwfz907w4VK8KcOZBPfk2FbaTLReQZW85uAeDRCo/e9PrFMxdZELoAv2p+9PjJyTNBbxURAY89Bt7exlrnRYuaV4vIdeR//SLPiLgUQSGvQje1zm+cCdpvcT/zb4K+8YaxRO6OHVC+fPbnC3EDaaGLPCMuOY6i+f9p8WqtWfqMMRO018xe+FUz+Sborl0waxb8619Qq1b25wtxCwl0kSdordkSsYV6pf5Z2Gr7l9vZP2s/rce3plqXauYVl5wMr74KTZtCsWLw4ovm1SJyNQl0kSccjj/MmUtn6FLVGAJ4ct3Jf26CmjUTNCMD5s0zWuMffWQswHXoEJQpY049IteTQBd5wobTGwBoV6kdF08bN0H9g/zNuQlqscCMGVC7NoSGGjdAN2ww9gn193duLcKlSKCLPGHr2a2U8S1DWc+yzOkxB52hzbkJunUrVKtmbFbh6Wm00Pftg1YmbJohXI6MchF5wrnEc1QuUpmlTy8lZl8MA1YMoFiVYs4v5JNPIDERliyBrl1lSr+wK/nbJPKEhOQEqq+tzoE5B2j7QVuqdKxiQhEJRmv8oYeMiUMS5sLO5G+UcHmJqStkAD8AABqzSURBVImkbk0lYE4AtfrWovkrzZ1bQFISTJhgbCF3+jQ8/rhzry/yDOlyES5tb/RePpr7ET0X9KRg9YJ0n9LducvhbtkCvXsby+CGhBjBLmPMhYNIC124HK01a0+u5bGfH6PpV03xn+BPfs/8DFsxDM+Cntm/gT1NmgQpKUawL14sYS4cSlrowqVYrBb6/9Kf+YfmU9q7NOP+GAcXYNCaQRSraMJN0AsXjCn8zZo5/9oiz5EWunAZWmtGLBvB/EPzea/1e0xLmIZ1m5WOX3WkwqMVnF9QRgb8+acxTFEIJ5BAFy7jP5v/w+Tdk3nzkTfpEdWDbR9to+HwhgT/O9icgmbPNvrOQ0PNub7Ic2wKdKVUR6XUUaXUcaXUq3c5r7FSKkMpJbfxhVNtj9zOW+vfol/tfowoNIIlTy0hsEUgnf/b2fl7gqakwPPPw6BBUL++7DgknCbbQFdKuQPfAJ2AmsATSqmadzjvI2CVvYsU4m42ntlIz7k9KVe4HJ81+oy5Pebi7e9NnwV9cPd0d24xly5BkyYwcSI89xxs2wYFCji3BpFn2dJCbwIc11qf1FqnAXOAkCzOGw38AsTasT4h7ijDmsH4P8bT+qfW+Hr6srjXYn7r/xtXz1+l35J++JT0cX5R06cbOw4tXAhffQX58zu/BpFn2RLoAcDZG55HXnstk1IqAOgJfHe3N1JKDVdK7VRK7YyLi7vXWoUA4FLKJSb+OZFa39bi7Q1v079Of3Y+s5NTb54iclskIdNCKN2gdPZvZG8WC3z/PTRqBD17Ov/6Is+zZdhiVh2Q+pbnXwKvaK0z7tZfqbX+AfgBIDg4+Nb3EOKuDsYe5Juwb5i+dzpJ6Uk8VPYh5veZT+8avdn22Tb2/rSXVm+3olYfk8Z6f/ABHDwIc+eac32R59kS6JFAuRuelwWibjknGJhzLcz9gc5KKYvWerFdqhR5VoY1g2Xhy5i4YyK/n/odL3cvnqjzBCMbjyS4jDF65ejSo6x5eQ01+9Sk1VsmrVq4cye8+66xpnnfvubUIPI8WwI9DKiqlKoInAP6Af1vPEFrXfH690qpacAyCXPxIC6mXGTyrsl8HfY1py+epmyhsnzQ5gOebvg0xQsWzzwvZl8MC/svpEyjMvSYZuIGzxMmQJEi8M035lxfCGwIdK21RSk1CmP0ijswRWt9UCk14trxu/abC3GvUi2pNJnUhGPnj/FI+Uf4tP2nhASFkM/t5r+uV2KuMLvbbLwKedFvST88vD3MKfj0aWM53Ndfh0KFzKlBCGyc+q+1XgGsuOW1LINcaz3kwcsSednXO77m2PljLA5dTEhQVgOqIP1qOnN7zCUpLomhm4biW8bXyVXe4MQJ0BratzevBiGQmaIih0lOT+aDzR/wWOXH7hjm2qpZMmQJkX9G0uvnXpRpZPIenPra/X1Z31yYTBbnEjnK7P2zOX/1PK+3fP2O5/z+5u8cnHeQdh+3o0avGk6s7g5OnDAeZXNnYTIJdJEjpFhSmPjnRCZsmkCdEnVoGdgyy/N2T93N5g820/CZhjR7MYesYBgWBn5+xgYWQphIAl2YSmvN7AOzeW3da0RciqBL1S582fHLLNdfObX+FMuGL6NSu0p0/saENVruJDzcWOc8p9Qj8izp9BOmmrFvBgMWDsCvgB/rnlzHsv7LqFLs9v0+44/EM6/XPPyq+dFnfh/cPZy8RsvdnDplrHkuhMmkhS5Mo7Xmi+1fULtEbXYO34mbyrp9cSX6Cj93/Bl3T3eeWPYE+YvkoPVREhIgMhJq3rZenRBOJy10YZpVJ1axJ3oPoxqPumOYp11JY1aXWSTHJdN/eX+KVizq5CrvIjERxo83vm/a1NxahEBa6MIkqZZURv82mmp+1RhSf0iW51gtVhaELiB6TzT9fu1HmeAcMookLQ1++MEI87g4GDAAWrQwuyohJNCFOb7/63uOnz/OqoGr8MrnddtxrTXLn13OsRXH6Pp9V6p1yUHbuLVtC5s3w6OPwscfQ+PGZlckBCCBLkxgsVr4fNvntAhsQYfKHbI8Z9OETeyatIsWr7eg0fBGTq4wG1u2wLPPwtdfy8gWkaNIH7pwugWHFnDm0hleavZSlsd3T93N+jfXU3dgXdq838bJ1WXDYjFmhvr5SZiLHEcCXTjV9ZEt1fyq0bVa19uOhy8LZ+kzS6ncoTLdJ3fPOWPNr/vjD+Oxfn1z6xAiCxLowqm+/+t7dpzbwZimY24b2XJm0xnm95lP6Qal6ftLX+fvB2qLefOgYEHo1MnsSoS4jfShC6dISktiafhSRq0YRddqXRneaPhNx6P3RjO722wKBxam/4r+ePp4mlTpXVy5ArNnQ69esvGzyJEk0IXdpWeksz92P2HnwthxbgdhUWEcjDuIVVupV7Ies3rNwt3tn9b3+RPn+fmxn/Hy9WLQmkEULF7QxOrv4r//Ncae/+tfZlciRJYk0MUDS89IZ8GhBWyL3EZYVBi7/95NakYqAH4F/Ggc0JieQT1pHNCY1hVaU9Dzn8BO/DuRGe1nYLVYGbx+MIUDC5v1Me5u2TIYNw4efxya5ZBFwYS4hQS6eCAZ1gwGLx7M7AOz8fbwplHpRoxqMorGZRrTOKAxFYtUvOONzeSEZH7u8DNJsUkMXj+Y4jWKZ3meqdLSjJb5W29Bw4bw008yukXkWBLo4r5prRn922hmH5jNhDYTeLn5y7dtE3cnqZdTmdlxJgnHEhiwYgABjQMcXO19WLMGnnsOjhyBLl3gxx/B29vsqoS4IxnlIu7bR1s+4n87/8dLzV7i9Zav2xzmaUnG+izRe6LpM78PFdtUzP6HnMVqNbpXWreGDh0gPd14vmwZlCpldnVC3JUEurgv8w7O47V1r9Gvdj8+aveRzT93fS/Qs1vP0mtmL6p3q+7AKu9BcjJ89x3UqAHdusHx4/DZZ3DggNE6FyIXkC4XcU8uXL3Ay2te5sfdP/Jw2YeZGjLV5sk/llQL83rN4+S6k/SY1oNafWs5uFobLVliTOWPioJGjWDWLOPmp4eH2ZUJcU8k0IVNzl89z5IjS3ht3WvEJ8fzcrOXefvRt8mfz7a1yS0pFuY9Po/jK4/T7cdu1HuynoMrtkFMDIweDfPnQ926MHMmtGolNz1FriWBLrKUYc0gLCqMlcdXsvL4SsKiwrBqKw1LN+S3Ab/RoHQDm98rLSmNuT3mcnLtSbp814WGwxo6sPJ70LEjHDoE778PL78sLXKR60mgi0zJ6cn8cugXlh1bxpoTa7iQcgGFoklAE9585E06VulI4zKNb5oUlJ3UxFRmdZnF2S1nCZkWQv3BOWgNlJMnjUlC48aZXYkQdiGBLjgYe5Dv//qe6Xuncyn1EqV9ShMSFELHyh1pV6kdft5+9/W+Vy9cZWanmUTtjKLXrF7UDq1t58ofQESEcSNUpvALFyKBnkdprVlwaAETd0xkc8RmPN096V2jNyOCR9AysOUDr3KYFJfEzx1+Ju5QHH1/6UtQSJCdKreDtDTo29cI86efNrsaIexGAj0P+ivqL55b+Rxbz26lctHKfNzuY4bUH0LxgvaZqZn4dyIz2s3gwskL9FvSjyodq9jlfe1m7Fj480/jZmjVqmZXI4TdSKDnIfHJ8by69lWm7J5C8YLF+bHbjwypP+Se+sSzc+nsJaa3nU5iVCL9V/SnYuscNGkIYPJk+OYbeOEFY2iiEC5EAj0P0Foz7+A8Rv82mgspFxj78FjefORNCue370JYF05eYHrb6Vw9f5VBqwdRrlk5u77/A7t4EUaOhPbt4cMPza5GCLuTQHdx8cnxDPt1GL8e/ZXGZRqzrvs66pSsY//rHI1netvpWK5aePL3JynTqIzdr/HAtm6F1FR47TXIJ3/1heuxaeq/UqqjUuqoUuq4UurVLI4PUErtu/a1VSmVA2aNiDMXz9B8SnNWHV/FZx0+Y9uwbQ4J8+i90UxrNQ1rupXBGwbnzDAHCAszJg01aWJ2JUI4RLbNFKWUO/AN0B6IBMKUUr9qrQ/dcNopoJXW+oJSqhPwA9DUEQWLu7ucepkDsQfYH7Of8RvHk5yezNon19IisIVDrnfq91PM7TkXr0LG5hT+Qf4OuY5dxMRAsWLGFnJCuCBb/t3ZBDiutT4JoJSaA4QAmYGutd56w/nbgbL2LFLcLj0jnfCEcPbH7md/zH7jMXY/py+ezjynfOHybByy0SGtcoD9s/ezZMgSilUtxsCVAylUtpBDrmM3kZFQPAeuuS6EndgS6AHA2RueR3L31vcw4LesDiilhgPDAQIDA20sUVy35sQapu+bzv6Y/RyOP0xaRhoA7sqd6v7VaRrQlKcbPE2dknWoU6IO5YuUv20jZnvQWrPl4y2se3Ud5VuVJ3RRKAWK5vAJOmlpsGED9OtndiVCOIwtgZ7VDBOd5YlKtcYI9Cz/fa+1/gGjO4bg4OAs30NkbXPEZrrO7kqR/EVoVLoRj1V+LDO4g/yD8Mrn5ZQ6MtIyWPbvZeyZsofaT9QmZGoI+bxywQ3GbduM/UA7dza7EiEcxpbfxEjgxvFnZYGoW09SStUFfgQ6aa0T7FNe3mXVVs5cPMOR+CMciT/CB5s/oHzh8mx/ejvFChQzpaar568y7/F5nF5/mkfeeoRH33n0gWeUOs2mTcYN0VatzK5ECIexJdDDgKpKqYrAOaAf0P/GE5RSgcBCYJDWOtzuVbqw5PRkwhPCM4P7+tfRhKOkWFIyz6tYpCLL+i8zLcxjD8Yyt+dcLp25RM8ZPak7sK4pddy3sDAICoKiRc2uRAiHyTbQtdYWpdQoYBXgDkzRWh9USo24dvw74C3AD/j2WovNorUOdlzZudtnWz9j7am1HIk/wpmLZ9DXerAUiopFKxLkH0S7Su0I8g/K/PL3Nm/0yMH5B1kydAmePp48ue5JAlvkwvsfCQlQurTZVQjhUEprc7qyg4OD9c6dO025tpnSMtLwet+LwMKBNC/XnCD/IGr41yDIP4iqflVt3jDCGawWK2tfW8u2T7dR9uGy9F3QF98yvmaXdX8qVYLgYJg3z+xKhHggSqm/7tRgzgV3s1xLQrJxe+GlZi8xqskok6u5s6S4JBaELuD0+tMEPxtMxy864u5pvzVfnCoiAk6dgueeM7sSIRxKAt2JktOTGbBwAG7KjUalG5ldzh2dCzvHvN7zSIpNImRqCPWH5KBNKe5FQoKxENfXX4ObGzz2mNkVCeFQEuj3SWtNakYqV9KukJiaSGJaIompicbza98npiXedPzPc3/yV9RfTO85nYfLPWz2R7iN1prdk3ezYuQKfEr7MGzrMEo3zIX9zlFR8J//wJQpxiYWXbrAq69CjRpmVyaEQ+WZQL8ewLaGb2JqIlfS7x7WFqvFpmt7uHng6+VLkfxFmNFzBgPqDnDwp713KRdTWPavZRycd5BK7SvRe3ZvvP28zS7r3p07B82bG6E+YAC8+CLUqmV2VUI4hUsEeoY1gx3ndrAsfBmH4g/dMazvNYB9PX3x9fLFx9OHQl6FCCgUgK+n8fz6sRvPufE1H0+fzO+dNennfp3ZdIZFAxeRGJVImwltaP5Kc9zc7T/D1OHOnze6VRISjJUVg2Wglchbcm2gX0q5xOoTq1l2bBkrjq0gPjk+cwp8Ya/CFM5fODOA7xS6dwrinB7A9mK1WPlj/B9smrCJIhWL8NSWpwhoEmB2WfdHa3jySTh2DFaulDAXeVKuC/RNZzbxzh/vsPHMRixWC8UKFKNTlU50rdaVxyo/RtECMnHEFvFH4lk8ZDHn/jxHvcH16PTfTnj55uL/kX39NSxfDhMnQuvWZlcjhClyXaADxFyJ4YWHX6Brta48VPYh8rnlyo9hCmuGle1fbOf3N37Hs6Anvef0pnZobbPLun/p6fDGG/Dxx8bNz1E5dyioEI6W65KwRWALDjx7wOwycqW4Q3H8+vSvRG6LpHpIdbp+1xWfUj5ml3X/Dh6Ep5+G7dthxAj4/HNjvRYh8qhcF+i5ZjGoHMSSYmHjhI1s+WgLXr5e9Py5J3X618mdf5YWCyxdanSx/P47FCpkzP7s08fsyoQwXa4LdHFvTq49yYqRK0gIT6DuwLp0+LwDBYvnwh174uLgxx/hf/+Ds2ehXDljrPmwYbJphRDXSKC7qIunL7LmpTUcWnCIopWLMnD1QCq3r2x2Wfdu9274739h1ixjg+e2bY0bn127ykbPQtxCfiNcTOLfiWyasIm/fvgLt3xutH6vNc1ebEa+/LnsP/W6dfDuu8Y65t7eMHSoccNTJgkJcUe57Ldc3MnV81fZ8vEW/pz4J9Z0Kw2GNeCRNx7J+ft8ZiUy0miBlywJn31mhLmsYy5EtiTQc7m0K2ls/3I7Wz/ZSmpiKnX61+HRdx6lWBVzNsKwi7feAqvV2AO0QgWzqxEi15BAz6UsKRZ2fr+TTRM2kRyXTPWQ6rR+rzUl65Q0u7QHc/Ys/PSTsdSthLkQ90QCPZexWqzs+WkPf7z7B5fPXqZi24q0mdCGsk3Lml3ag9PaGEuuNTz/vNnVCJHrSKDnEhnpGRycd5CN4zeSEJ5AQJMAQqaGUKltJbNLs4/t242VEbdsgX79pHUuxH2QQM/hLp+7zJ5pe9j1wy4uRVyiRO0ShC4OpXr36rlzYtCtLl2C4cONyUGlSsEPPxg3QYUQ90wCPQfKSMvg6NKj7J68mxOrTqCtmgqtK9D5m85U7VwV5eYCQX7dJ5/A/Pnw9ttGC90nFy9FIITJJNBzkNiDseyevJt9M/aRHJ+Mb4AvLV5vQf0h9SlWORePWrmT1FSYNAm6dYN33jG7GiFyPQl0k6VeTuXAnAPsnrybczvO4ebhRlBIEPWfqk/lDpVz50YTtvryS4iNlRugQtiJBLoJtNZEbIpg95TdHJp/iPTkdIrXKk6HzztQd2Dd3LnWyr2KjIT33oOQEGjTxuxqhHAJEuhOkpaURvTuaI4uPcrhXw5z4cQFPH09qTOwDg2HNaRM4zKucZPTVh99BGlp8MUXZlcihMuQQLczrTWXz14mem80MXtjiNkbQ/TeaM4fPw8a3PK5UbFNRVq91Yqaj9fEw9vD7JKdLyEBpkwxNnGuWNHsaoRwGRLoD8CSYiH2YGxmaMfsjSFmXwwpF1IyzylaqSgl65WkzoA6lKpXivKtylOgaAETq84BvvsOkpONUS1CCLuRQLeB1por0VduDu69McQfjUdnaAA8vD0oWbcktfrWomTdkpSsV5KSdUriVSgX79PpCFevwjffwGOPycqJQtiZBPotMtIyiDscR8y+mJu6TJLjkjPPKRxYmJL1ShLUK4hS9UpRsl5JilUu5lrjwx3BaoXBgyE6GmbPNrsaIVxOng705Pjkm/u690QTdzgOa7oVAHcvd0rULkG1btUyg7tk3ZLSZXIvtIbTpyEsDBYvNiYRffIJtGpldmVCuJw8EehWi5WE8ISb+rlj9saQGJWYeY5PaR9K1StFlU5VKFmvJKXqlcKvmh9u+Vx4HLgjxMYa4R0WBjt2GI/x8cYxLy8YOxZeeMHcGoVwUS4X6CkXUzKD+/pj3ME4LCkWANw83CheszgV21bMDO6S9UrmjbHfjhIRAa+8Atu2wZkzxmtublCzJnTvDo0bQ5MmULs2eHqaW6sQLsymQFdKdQS+AtyBH7XWH95yXF073hlIBoZorXfZuVYA4o/GE740nOT45H++4v75/ur5q5nnehf3plS9UjQe2TgzvP2D/HH3dHdEaXlTTIyxz2dMDHTuDKNHG+HdoIGsyyKEk2Ub6Eopd+AboD0QCYQppX7VWh+64bROQNVrX02B/117tLvYA7GseWkNbh5uePt7U7B4Qbz9vSlVvxTexb0pVK5QZqvbp5RP3pqs42xhYTBsGERFwdq18PDDZlckRJ5mSwu9CXBca30SQCk1BwgBbgz0EGC61loD25VSRZRSpbXWf9u74GqeZ3i1+kI83SwoBViA6GtfwnksFggPN/b6XLxYwlyIHMCWQA8Azt7wPJLbW99ZnRMA3BToSqnhwHCAwMDAe60VgHz+RchXp9p9/ayws6eegmefBV9fsysRQmBboGfVZ6Hv4xy01j8APwAEBwffdtwmDz9sDH0TQghxE1vG5EUC5W54XhaIuo9zhBBCOJAtgR4GVFVKVVRKeQL9gF9vOedX4ElleAi45Ij+cyGEEHeWbZeL1tqilBoFrMIYtjhFa31QKTXi2vHvgBUYQxaPYwxblE0hhRDCyWwah661XoER2je+9t0N32tgpH1LE0IIcS9kXrsQQrgICXQhhHAREuhCCOEiJNCFEMJFKON+pgkXVioOOGPKxW3nD8SbXYRJ5LPnTfLZc77yWuviWR0wLdBzA6XUTq11sNl1mEE+u3z2vMYVPrt0uQghhIuQQBdCCBchgX53P5hdgInks+dN8tlzMelDF0IIFyEtdCGEcBES6EII4SIk0LOhlPpEKXVEKbVPKbVIKVXE7JqcRSnVRyl1UCllVUrl6uFctlBKdVRKHVVKHVdKvWp2Pc6klJqilIpVSh0wuxZnUkqVU0qtV0odvvZ3/Xmza3oQEujZWwPU1lrXBcKB10yux5kOAL2AjWYX4mg3bIbeCagJPKGUqmluVU41DehodhEmsAAvaK1rAA8BI3Pzf3cJ9GxorVdrrS3Xnm7H2I0pT9BaH9ZaHzW7DifJ3Axda50GXN8MPU/QWm8Ezptdh7Nprf/WWu+69n0icBhjP+RcSQL93jwF/GZ2EcIh7rTRucgjlFIVgAbAn+ZWcv9s2uDC1Sml1gKlsjg0Tmu95No54zD+eTbTmbU5mi2fPY+waaNz4ZqUUj7AL8AYrfVls+u5XxLogNa63d2OK6UGA12BttrFBu5n99nzENnoPI9SSnlghPlMrfVCs+t5ENLlkg2lVEfgFaC71jrZ7HqEw9iyGbpwMUopBUwGDmutPze7ngclgZ69rwFfYI1Sao9S6rvsfsBVKKV6KqUigYeB5UqpVWbX5CjXbnxf3wz9MDBPa33Q3KqcRyk1G9gGVFdKRSqlhpldk5M0BwYBba79fu9RSnU2u6j7JVP/hRDCRUgLXQghXIQEuhBCuAgJdCGEcBES6EII4SIk0IUQwkVIoAshhIuQQBdCCBfx//UqhQIS18OZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# One iteration wtih Normal distribution\n",
    "\n",
    "n = 100\n",
    "alpha = 0.05\n",
    "r = norm.rvs(size=n)\n",
    "epsilon = math.sqrt((1 / (2 * n)) * math.log(2 / alpha))\n",
    "\n",
    "F_n = lambda x : sum(r < x) / n\n",
    "L_n = lambda x : max(F_n(x) - epsilon, 0)\n",
    "U_n = lambda x : min(F_n(x) + epsilon, 1)\n",
    "\n",
    "xx = sorted(r)\n",
    "\n",
    "df = pd.DataFrame({\n",
    "    'x': xx, \n",
    "    'F_n': np.array(list(map(F_n, xx))), \n",
    "    'U_n': np.array(list(map(U_n, xx))), \n",
    "    'L_n': np.array(list(map(L_n, xx))), \n",
    "    'CDF': np.array(list(map(norm.cdf, xx)))\n",
    "})\n",
    "df['in_bounds'] = (df['U_n'] >= df['CDF']) & (df['CDF'] >= df['L_n'])\n",
    "\n",
    "plt.plot( 'x', 'L_n', data=df, color='red')\n",
    "plt.plot( 'x', 'U_n', data=df, color='green')\n",
    "plt.plot( 'x', 'CDF', data=df, color='purple')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "87e150145c4e4054ad60e0672ec73b97",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Average fraction in bounds: 0.963\n"
     ]
    }
   ],
   "source": [
    "# 1000 iterations with Normal distribution\n",
    "\n",
    "bounds = []\n",
    "for k in tqdm_notebook(range(1000)):\n",
    "    n = 100\n",
    "    alpha = 0.05\n",
    "    r = norm.rvs(size=n)\n",
    "    epsilon = math.sqrt((1 / (2 * n)) * math.log(2 / alpha))\n",
    "\n",
    "    F_n = lambda x : sum(r < x) / n\n",
    "    L_n = lambda x : max(F_n(x) - epsilon, 0)\n",
    "    U_n = lambda x : min(F_n(x) + epsilon, 1)\n",
    "\n",
    "    # xx = sorted(r)\n",
    "    xx = r # No need to sort without plotting\n",
    "    \n",
    "    df = pd.DataFrame({\n",
    "        'x': xx, \n",
    "        'F_n': np.array(list(map(F_n, xx))), \n",
    "        'U_n': np.array(list(map(U_n, xx))), \n",
    "        'L_n': np.array(list(map(L_n, xx))), \n",
    "        'CDF': np.array(list(map(norm.cdf, xx)))\n",
    "    })\n",
    "    all_in_bounds = ((df['U_n'] >= df['CDF']) & (df['CDF'] >= df['L_n'])).all()\n",
    "    bounds.append(all_in_bounds)\n",
    "    \n",
    "print('Average fraction in bounds: %.3f' % np.array(bounds).mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x26bbe413eb8>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZScdZ3v8fe3q6r37iSkO2sn6SRGoVEISQQdZURQZBu2UQHvVZZxMswBDnEFL+dedZRzFLgM4wFhAig3w1wREBWdDOAg6I0KSYAQSAKahEA6SSedbrL0Vuvv/vFUOtWdXiqd6nrqqfq8zqmTqud5qvqTTtcnv/7Vs5hzDhERCb4yvwOIiEhuqNBFRIqECl1EpEio0EVEioQKXUSkSIT9+sINDQ2uubnZry8vIhJIL7300l7nXONQ63wr9ObmZtauXevXlxcRCSQze3u4dZpyEREpEip0EZEioUIXESkSvs2hDyUej9Pa2kpfX5/fUcZFZWUlTU1NRCIRv6OISBEqqEJvbW2lrq6O5uZmzMzvODnlnKOjo4PW1lbmzp3rdxwRKUKjTrmY2Y/MbI+ZvT7MejOzH5jZZjNbb2aLxhqmr6+PyZMnF12ZA5gZkydPLtrfPkTEf9nMoT8EnDPC+nOBBenbUuDeYwlUjGV+SDH/3UTEf6NOuTjnfm9mzSNschGwwnnn4X3BzCaa2XTn3K4cZRSRMeqL9XHvP9/Lvr/sgzigs2UXhBPPPJHPXvXZnL9uLubQZwLbMx63ppcdUehmthRvFM/s2bNz8KVFZDg7tu3g9o/dzqR3JlGW8cu4U6v77rX4awVb6EPNIwz5E+OcWw4sB1iyZElB/lTV1tbS1dXldwyRY/b4TY9Tt6OOxlsbWfrlpYTKQ1iZpv2KWS72Q28FZmU8bgJ25uB1RWSMEtEE+1buY8OJG7j6a1cTrgyrzEtALkboTwLXm9kjwGnA/pzMny9bBuvWHfPLDLBwIdx1V05e6qqrrqK+vp61a9fS1tbGbbfdxqc//emcvLbIsdr18i7ogq7TuqiKVPkdR/Jk1EI3s58AZwANZtYKfBOIADjn7gNWAucBm4Ee4OrxCltodu3axapVq3jjjTe48MILVehSMN559R0AFn1kzHsRSwBls5fLFaOsd8B1OUt0SI5G0uPp4osvpqysjJaWFnbv3u13HJF+f/h/fyAejnPJmZf4HUXyqKCOFA2aioqK/vve/2sihWHbum1UzKjg1KZT/Y4ieaSTc4kUmUQqQWVrJfXvrdfBbCVGhT5IT08PTU1N/bc777zT70giR6Xt3TYmHJhA7dxav6NInmnKZZBUKpXVdg899NCAx9p3XQrF1s1bAZg8a7LPSSTfNEIXKTJvv+VdoWxG8wyfk0i+aYQ+iltvvZXHHntswLLPfOYz3HLLLT4lEhnZgbYDADQ2DXkdYSliKvRR3HLLLSpvCZRYRwyA2imaQy81mnIRKTJ72/YCMLNpps9JJN9U6CJFZnfbblJlKWrqa/yOInmmQhcpMl0dXaTqUtoHvQSp0EWKSCwZI7U/Rag+5HcU8YEKfZBt27bx/ve/f8Cyb33rW9xxxx0+JRLJ3sHoQSr7KglP1P4OpUiFLlJE4qk4Vb1VhOo0Qi9FBfvf+LKnlrGuLbfnQ184bSF3nZPbszieccYZnHbaaTz33HPs27ePBx98kNNPPz2nX0MkW8lUkopoBWW1GquVIv2r50AikWD16tXcddddfPvb3/Y7jpSwlEt5hV6jt3YpKtgReq5H0tkabs+AkfYYuPTSSwFYvHgx27ZtG49YIll5e//bVEQrqD+u3u8o4gP9Nz7I5MmTeffddwcs6+zspKGhYdjnHDoveigUIpFIjGs+kZFs2LWBSCLCtCnT/I4iPlChD1JbW8v06dN59tlnAa/Mn3rqKT760Y/6nExkdFt3eGdabGzQeVxKkQp9CCtWrOC73/0uCxcu5Mwzz+Sb3/wm8+fP9zuWyKh27N4BQGV9pc9JxA8FO4fup5aWFp577rmstn3++ef77zc0NGgOXXwV6/ZOzFVeW+5zEvGDRugixaTP+yNSE/E3h/hCI/QsXXfddfzhD38YsOzGG2/k6quv9imRyJFSPd4Vt8prNEIvRSr0LN1zzz1+RxAZVarXK3SN0HMkkYC+vtFvvb3ZbXfodvHF8PnP5zyuCl2kiFifd7xE0YzQUymIRo++MMdSskM9J5k8tvzhMFRWHnnr7MzN92fwlxuXVxURX/QXeiF/KLp9O9x/P7z44ugFG4sd29cyG7pQq6oO36+vH3qbkZ6T7S2c34pVoYsUkd6uXqAAp1xSKfjNb+Dee+FXvwLn4JRToK4OJk06ttIc6TmRiFfqJUKFPoS2tjaWLVvGmjVrqKiooLm5mbvuuouTTz6Z448/nr6+Purq6rjuuuu48sorAXjooYf42te+xsyZ3mW/TjrpJFasWOHnX0NKUM/BHgAi1T4VejzujcC3bRt4++MfYcsWaGyEm26CpUuhudmfjEVMhT6Ic45LLrmEK6+8kkceeQSAdevWsXv3bubPn88rr7wCwNatW7n00ktJpVL9e7pcdtll3H333b5ll9KWTCUP77ZYNU6FPlxhH7rt2OGNxg8pK4OmJjj+ePjOd+DSSyF9qgzJPRX6IM899xyRSIRrr722f9nChQuPOGBo3rx53HnnnXzlK1/RrotSELrj3UTiESgHK8vBNINz8Npr3hTJs896I+zW1qELu7kZPv5x78/MW1OTN+0heVGwhf7UsqdoW9eW09ectnAa59x1zojbvP766yxevDir11u0aBFvvPFG/+Of/vSnrFq1CtA+6pJ/e3v2Ek6EKas8huMFo1H43e+8Ev/Vr+Dtt73lixbBxz52uKjnzlVhF6CCLfQgcM4NeKwpF/FTe3c7kXiEUNVRXq0okYBnnoF/+zf49a+hq8v7oPETn4BbboHzz4cZM8YntORUVoVuZucA/wKEgAecc98btH4C8DAwO/2adzjnfnwswUYbSY+XE088kccffzyrbV955RVOOOGEcU4kkp19ffuIxCPZfSDqHLz8slfiP/kJ7NkDkyfDFVfAhRfCWWd5pS6BMurvZmYWAu4BzgVagCvMrGXQZtcBG51zJwNnAP/bzAp4R9jhnXnmmUSjUe6///7+ZWvWrOHtQ796pm3bto2vfvWr3HDDDfmOKDKkA9EDXqGP9oHoa6/BySfDkiXeboSnnw6/+AXs3AnLl8MFF6jMAyqbEfqpwGbn3FYAM3sEuAjYmLGNA+rMu6xPLdAJBPJKD2bGz3/+c5YtW8b3vvc9Kisr+3db3LJlC6ecckr/bos33HCD5smlYMSSMcKJMOHaEd7Wq1bB3/wNVFfDfffBZz/r7QcuRSGbQp8JbM943AqcNmibu4EngZ1AHXCZcy41aBvMbCmwFGD27NljyZsXM2bM4NFHHz1ieW9v77DPueqqq7jqqqvGMZXIyBKpBJF4hHDVMG/r11+HT34SZs/25sznzMlvQBl32XwcPtT+T27Q408B64AZwELgbjM74qKGzrnlzrklzrkljY26oopILqVcyhuhVw5T6KtWeYfT/8d/qMyLVDaF3grMynjchDcSz3Q18ITzbAbeAo7PTUQRyUYilSCcCA8/h37oWrnpo5ml+GRT6GuABWY2N/1B5+V40yuZ3gHOAjCzqcD7gK1jCTR4V8BiUsx/N/Hf/uh+IvEIldXDXH7uwAFvn3F94Fm0Ri1051wCuB54GtgEPOqc22Bm15rZocMpvwP8lZm9BjwL3OSc23u0YSorK+no6CjK4nPO0dHRQWWlrvUo42N/337CifDwhd7VBTU1+Q0leZXVfujOuZXAykHL7su4vxM4+1jDNDU10draSnt7+7G+VEGqrKykqanJ7xhSpFIuRTgZHn4/9J07Ydq0/IaSvCqoI0UjkQhz5871O4ZIIDkc4fgwH4rG47B+PSxYkP9gkje6SLRIkYglvP3QQxVDHPr/pS/B5s2gXWuLmgpdpEhEo1HKXNmRe7msXAn33ANf+Yp3IJEULRW6SJFI9nnXvzxiyuXll70/b701z4kk31ToIkUi2hcFhij0jg5v7xZdWKLoqdBFisSujl3AEIX+yiugs4KWBBW6SJHYs28PMKjQOzth9Wr4yEd8SiX5pEIXKRLxvjjAwL1cbroJYjG45hqfUkk+qdBFioSLekdY94/QX38dHngAvvxlOOkkH5NJvqjQRYqFN0AnXJEu9NWrvT//4R/8ySN5p0IXKRYx748BI/SKCu9izlISVOgiReKIOfRnnvE+DA0d5UWjJbBU6CJFIJaMDZxy2bwZNmyA88/3N5jklQpdpAh0xboIJ7ypllBFCO6/3xuZX3aZz8kkn1ToIkWgO9bdX+jhyjCsWAEXXKCrE5UYFbpIEeiOdxNKenPl4VQc2trg1FN9TiX5pkIXKQIDplwOdHoLp071MZH4QYUuUgT29e07PELfnb6G++zZPiYSP6jQRYrAgeiBwyP0Q4U+a5aPicQPKnSRIhBLxrwRukFZvM9bWFvrbyjJOxW6SBHoS/T1X37O+tKFHhnmYtFStFToIkXgQPQAoWSIUHnIO+R/4kSYMsXvWJJnKnSRInBoL5dwRRjWrYNTTgEzv2NJnqnQRYpAX6KPUDJ0+LD/973P70jiAxW6SBGIJWNEUhFC5WXeNUTnzPE7kvhAhS5SBJKpJJFkhHAkPc2iPVxKkgpdpAgkXdL7UDScfkuHwyM/QYqSCl2kCDjnCCfD3pQL6BzoJUqFLlIEDo3Qw+H0lItG6CVJhS5SBPoSfUSSEUKHjhJ9z3v8DSS+yKrQzewcM3vTzDab2c3DbHOGma0zsw1m9rvcxhSRkcSSMW/KJdrtLWhp8TeQ+GLU38vMLATcA3wSaAXWmNmTzrmNGdtMBH4InOOce8fMdIiaSB71fyga7YH6epg82e9I4oNsRuinApudc1udczHgEeCiQdt8DnjCOfcOgHNuT25jishIeuO9hJNhwrEeaGryO474JJtCnwlsz3jcml6W6b3AJDN73sxeMrMvDPVCZrbUzNaa2dr29vaxJRaRI3T2dnpTLiShosLvOOKTbAp9qBNCuEGPw8Bi4HzgU8D/NLP3HvEk55Y755Y455Y0NjYedVgRGVp/oVtKuyyWsGz2bWoFMs+U3wTsHGKbvc65bqDbzH4PnAz8OScpRWREnb2dlCXKCLkElJf7HUd8ks0IfQ2wwMzmmlk5cDnw5KBtfgmcbmZhM6sGTgM25TaqiAznYOxgutDjKvQSNuoI3TmXMLPrgaeBEPAj59wGM7s2vf4+59wmM3sKWA+kgAecc6+PZ3AROaw33gsJCKEpl1KW1eFkzrmVwMpBy+4b9Ph24PbcRRORbCWSCSxhhFNxXamohOlIUZEiEEp6o/Ky7v0wbZrPacQvKnSRgEumkpDw7oe7VOilTIUuEnC9id7+EXooFdO50EuYCl0k4PoSfYQT3sdhIVIq9BKmQhcJuFgydniETgKm6FRKpUqFLhJw0US0v9DDJKGhwedE4hcVukjADRyhJ2HCBJ8TiV9U6CIBt7t798BCr6vzOZH4RYUuEnBbOrcMnEPXgUUlS4UuEnDtPe0Ze7kkdeh/CVOhiwRcb7x34JSLCr1kqdBFAq4n3kMllUB6LxdNuZQsFbpIwHXHu6mmGoAykjBxos+JxC8qdJGAOxA9QK15R4eGLQVVVT4nEr+o0EUCrrO3k7oyb1fFUNngq0NKKVGhiwRcb6KXSufNoYcs5XMa8ZMKXSTgkqkkoVR6L5cK7eFSylToIgHXm+ilIlUBQCiiQi9lKnSRgOuKdVHh0oVerrd0KdO/vkjAHYwePDyHXqM9XEqZCl0k4A5ED1DpKjFSlM2d43cc8ZEKXSTgehO9RFIRbw+XqVP9jiM+UqGLBFgylew/H3qIJFRW+h1JfKRCFwmwTXs3AVBfVu+dx6W83OdE4icVukiArd6xGoCG8gZCpjMtljoVukiA7evbB0BFqoKQ05kWS50KXSTAkqkkAKlYkpCL6/JzJU6FLhJgnb2dRMoiEE0RJgE1NX5HEh+p0EUCbFfXLqbVTiPZE/P2clGhlzQVukiAbdu3jTkT55Dsi6vQRYUuEmTb9m2jeWIzyd50odfW+h1JfJRVoZvZOWb2ppltNrObR9jug2aWNLNP5y6iiAynN9FLTaTm8AhdhV7SRi10MwsB9wDnAi3AFWbWMsx23weeznVIERlaMpWkzMpIRuPeh6KTJvkdSXyUzQj9VGCzc26rcy4GPAJcNMR2NwA/A/bkMJ+IDCOaiNLR28G02mkkoknNoUtWhT4T2J7xuDW9rJ+ZzQQuAe4b6YXMbKmZrTWzte3t7UebVUQytB5oBWDOhDkk4ymv0CsqfE4lfsqm0G2IZYOvRHsXcJNzLjnSCznnljvnljjnljQ2NmabUUSGsKfb+2V4au3Uw4VeXe1zKvFTOIttWoFZGY+bgJ2DtlkCPGJmAA3AeWaWcM79IicpReQI7T3eb7kN1Q0k484r9AkTfE4lfsqm0NcAC8xsLrADuBz4XOYGzrm5h+6b2UPAr1XmIuNr18FdAEyvnU4i4QiR0Ai9xI1a6M65hJldj7f3Sgj4kXNug5ldm14/4ry5iIyPXV1eoU+tnUoy4QiTgnA2YzQpVln96zvnVgIrBy0bssidc1cdeywRGU1nbycTKiYQshDJpBEKDf5oS0qNjhQVCaiO3g4mVU0iFU8BECpL+ZxI/KZCFwmoLZ1bmDdpHoloAoBwmUbopU6FLhJQW9/dyvxJ80nGvL2FQ2iEXupU6CIBFUvGqI5Uk4ymCz2sEXqpU6GLBFTKpTCsf8olFNHbudTpJ0AkoKLJKJXhyv4Rejg81EHdUkpU6CIBFE/GiSVjVEWqDo/Q9aFoyVOhiwTQzoPe2Tdm1M04PEKv0Ym5Sp0KXSSADp1psam+6fBui/U67L/UqdBFAqj/TIs1Uw/v5TJRVysqdSp0kQDq6O0AvDMt9o/QJ6jQS50KXSSADkYPAlBfUX94hF5X5WckKQAqdJEAOhjzCr2mvObwCL1aH4qWOhW6SAB19HRQX1FPuCx8eIRerrdzqdNPgEgAtXW30VjtXcaxf4Qe8jORFAIVukgAbWzfyAmNJwCQ7IkBENJ+6CVPhS4SMNFElDf2vsH7G98PQKKrF4BwnfZDL3UqdJGAeW3PayRSCRZNXwRAojsKQLi20s9YUgBU6CIBs65tHUB/oSd7YhgpyirL/YwlBUCFLhIwG9s3UhWuYu6kuQAkemKESUCF5tBLnQpdJGBe3f0qLY0tlJn39k30qtDFo0IXCZBEKsGLrS/y4aYPH17WcYAQSWhq8jGZFAIVukiAvPXuW3THu1k8Y3H/smTnAW+EPm+ej8mkEKjQRQKks7cTgCk1U/qX9U+5TJzoVywpECp0kQB5Z/87gHfa3EMS0ZRX6JXabbHUqdBFAmTVO6uojlRz0tST+pcloknCZSkwXVO01KnQRQLk9+/8ntNmnkYkFOlfloinCId0PVFRoYsExv6+/bza9ip/PeevByxP9CUJR/RWFhW6SGC8uvtVHI4PNX1owPJE3BGu0FtZVOgigbG5czMA8yfNH7A8kYRw2I9EUmiyKnQzO8fM3jSzzWZ28xDr/5uZrU/f/mhmJ+c+qkjpenj9w9z41I3MqJvBnIlzBqyLJ8s0hy5AFoVuZiHgHuBcoAW4wsxaBm32FvAx59xJwHeA5bkOKlKKkqkk1/zyGj7/889zyrRTeOHvXqA8NPAkXImkEdGUi5DdCP1UYLNzbqtzLgY8AlyUuYFz7o/OuXfTD18AdAyySA48tfkpfrzux3z9r77Ob6/8LbMmzDpim0QqRLhSlyuS7Ap9JrA943Fretlw/g74z6FWmNlSM1trZmvb29uzTylSolasX8Hkqsl898zvEi47cqLcOUfchQhX69S5kl2hD3W0wpATdmb2cbxCv2mo9c655c65Jc65JY2NjdmnFClB9790P49ueJRrTrlmwH7nmbwLRBuRGhW6QDafjbcCmb/nNQE7B29kZicBDwDnOuc6chNPpPTsPLiTh9c/zDee/QbnLTiPW8+8ddht4z1xACKV2s1Fsiv0NcACM5sL7AAuBz6XuYGZzQaeAD7vnPtzzlOKFLlt+7bxxKYn+Nmmn/Gn7X/C4Thr7lk89pnHhh2dA8Q7DgAQqRl+Gykdoxa6cy5hZtcDTwMh4EfOuQ1mdm16/X3A/wImAz8073wSCefckvGLLRJ8KZfigZcfYPlLy3lp10sALJy2kH/6+D/xtyf8LSc0njDqa8TXeJeji7xnzihbSinI6vc059xKYOWgZfdl3P8i8MXcRhMpXm+9+xbXPHkNz297nsXTF3PbJ27j0hMuZf5x80d/cob4+k0ARFoWjEdMCRhNvInkQVesi03tm9jYvpH1u9fzry/9K2VWxoMXPsjVC6/GxnimxHh3DIDIxNpcxpWAUqGL5ND+vv1s2usVd+bt7f1v929THirnk/M+yQ/P/yGzJ8w+pq/X/6Fonc6FLip0CQjnHCmXIpFK9N/iqfiAx/3Lk0MvH+k5Iz1v1Oe4BLsO7mJj+0Z2HNzRn7kyXMnxDcfzkdkf4e8b/p4Tp5xIS2ML8ybNG3Kf8rGIt+8DaokcpxG6qNADaXCxHUuR5fQ5bnwz+KnMygiXhQmXhYmURfrvH7o1VDdw1ryzaGlooaXRuzVPbCZUNr5HcPYXerX2chEVelacc0STUbpiXf237lj3gMeDb91xb31voveoR4SjlZ8b+riuvMkstEjoyHIbqfiqwlUjPm+o52S1biyvl2X2UFmIMivMc6XE2zqAJhW6AEVY6IlUYkDZHirWYYv30LbxEdbFuki6ZNYZqsJV1JbXUlteS0W4or8cMgukPFROdaS6YEtsuOcUarGVqvi7XQAqdAECWOir3lnFHX+8Y9jC7kv0Zf1aIQtRV1FHTaSmv4Bry2uZXjt9wOPB6wffasprBmw73r9mixwSS3j/warQBQJY6D3xHt7a9xa15bVMqJzAzPqZoxfuMOvLQ+Vj3l1MpBDEEwY4FboAASz0s+efzdnzz/Y7hoj/9u8nFnWUV5gGJgLoEnQiwdXaSoxyyqsDNy6TcaJCFwmqvXtV6DKACl0kqPbsIUoFFfU6SlQ8KnSRoNq50yv042r8TiIFQoUuElRtbcSooFyFLmkqdJGg2rePaFkVFfUVfieRAqFCFwmq3buJWgXldbqeqHhU6CIB5f78F/pSFVRNqvI7ihQIFbpIQMV37CbljMpJ2stFPCp0kYDqS3iH+1dOVKGLR4UuElC9SRW6DKRCFwmo3qh3hGj15Gqfk0ihUKGLBFF3Nz0pb++W6gYVunhU6CJBtH073XgHFNVM0YFF4lGhiwTRjh10410YumqydlsUjwpdJIja2+mihupJ5YQiukKWeFToIkH07rscpJ66GXV+J5ECokIXCaKuLvYzgYlzJvqdRAqICl0kgNz+A+xjIvVzj/M7ihQQFbpIAHVv7yRKJcctUKHLYSp0kQDa80YnAI0tjT4nkUKiQhcJoB1vdgEwfdF0n5NIIcmq0M3sHDN708w2m9nNQ6w3M/tBev16M1uU+6giAkBnJ3/ZP4Wp00yH/csAoxa6mYWAe4BzgRbgCjNrGbTZucCC9G0pcG+Oc4pI2pblz7Kd2Xzggtl+R5ECE85im1OBzc65rQBm9ghwEbAxY5uLgBXOOQe8YGYTzWy6c25XrgNvvvWnPPOdF3P9siIFyzlwGA4jkQpxIFHNZOvgg9++zu9oUmCyKfSZwPaMx63AaVlsMxMYUOhmthRvBM/s2WMbXVQ01tPYMKanigSWWQozCJUlaKw/wOKrPkD5DL0RZKBsCt2GWObGsA3OueXAcoAlS5YcsT4bs5aey6yl547lqSIiRS2bD0VbgVkZj5uAnWPYRkRExlE2hb4GWGBmc82sHLgceHLQNk8CX0jv7fIhYP94zJ+LiMjwRp1ycc4lzOx64GkgBPzIObfBzK5Nr78PWAmcB2wGeoCrxy+yiIgMJZs5dJxzK/FKO3PZfRn3HaCP3EVEfKQjRUVEioQKXUSkSKjQRUSKhApdRKRImPd5pg9f2KwdeDuHL9kA7M3h6+VSoWYr1FygbGNRqLlA2cZiuFxznHNDnjfZt0LPNTNb65xb4neOoRRqtkLNBco2FoWaC5RtLMaSS1MuIiJFQoUuIlIkiqnQl/sdYASFmq1Qc4GyjUWh5gJlG4ujzlU0c+giIqWumEboIiIlTYUuIlIkAlfoZvYZM9tgZikzWzJo3TfSF6p+08w+lbF8sZm9ll73AzMb6oIcuc650MxeMLN1ZrbWzE4dLWc+mdkN6a+/wcxuK6Rs6RxfNTNnZg0Zy3zLZma3m9kb6Yug/9zMJhZCrowMI17IPY85ZpnZc2a2Kf2zdWN6+XFm9hsz+0v6z0k+ZgyZ2Stm9utCypa+dOfj6Z+zTWb24aPO5pwL1A04AXgf8DywJGN5C/AqUAHMBbYAofS61cCH8a6s9J/AuXnI+cyhr4N3auHnR8uZx+/hx4H/AirSj6cUSrZ0jll4p2t+G2gohGzA2UA4ff/7wPcLIVc6Qyj9decB5ek8Lfn+d0tnmQ4sSt+vA/6c/h7dBtycXn7zoe+fTxm/DPxf4NfpxwWRDfg/wBfT98uBiUebLXAjdOfcJufcm0Osugh4xDkXdc69hXdu9lPNbDpQ75z7k/O+KyuAi/MRFahP35/A4Ss4DZkzD3ky/SPwPedcFMA5t6eAsgH8M/B1Bl7G0NdszrlnnHOJ9MMX8K7K5XuutP4LuTvnYsChC7nnnXNul3Pu5fT9g8AmvOsLX4RXWKT/zMd78Ahm1gScDzyQsdj3bGZWD/w18CCAcy7mnNt3tNkCV+gjGO5C1TPT9wcvH2/LgNvNbDtwB/CNUXLm03uB083sRTP7nZl9sFCymdmFwA7n3KuDVvmeLcM1eL/pQWHkKoQMRzCzZuAU4EVgqktfxSz95xSfYt2FN1hIZSwrhGzzgHbgx+npoAfMrOZos2V1gYt8M7P/AqYNseoW59wvh3vaEMvcCMuP2Ug5gbOALznnfmZmn8X7n6+KyZEAAAJaSURBVPcT45nnKLKFgUnAh4APAo+a2bwCyfY/8KY3jnjaEMtymi2bnzszuwVIAP+er1xZKIQMA5hZLfAzYJlz7kAePrYalZldAOxxzr1kZmf4nWeQMLAIuME596KZ/QveFMtRv0jBcc59YgxPG+5C1a0c/vU4c/kxGymnma0Abkw/fIzDv+Ll5YLao2T7R+CJ9BTUajNL4Z0IyNdsZvYBvHnoV9MF0AS8nP5AedyzjfZzZ2ZXAhcAZ6W/d+QjVxYKIUM/M4vglfm/O+eeSC/ebWbTnXO70tOge4Z/hXHzEeBCMzsPqATqzezhAsnWCrQ6515MP34cr9CPKlsxTbk8CVxuZhVmNhdYAKxO/5py0Mw+lN675QvAcKP8XNoJfCx9/0zgLyPlzEOeTL9IZ8LM3ov3Acxev7M5515zzk1xzjU755rxfsgXOefa/M5mZucANwEXOud6MlYVwr9nNhdyz4v0e+xBYJNz7s6MVU8CV6bvX0l+3oMDOOe+4ZxrSv9sXQ781jn33wskWxuw3czel150FrDxqLP58WnuMX4SfAneGz0K7Aaezlh3C96n/W+SsScLsAR4Pb3ubtJHyI5zzo8CL+HtcfAisHi0nHn8HpYDD6e/Jy8DZxZKtkE5t5Hey8XvbHgfdm4H1qVv9xVCrowM5+HtUbIFb4rIr3+zj+JN96zP+F6dB0wGnsUb2DwLHOfzz9YZHN7LpSCyAQuBtenv3S/wpkWPKpsO/RcRKRLFNOUiIlLSVOgiIkVChS4iUiRU6CIiRUKFLiJSJFToIiJFQoUuIlIk/j+/JMfghOInLwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# One iteration wtih Cauchy distribution\n",
    "\n",
    "n = 100\n",
    "alpha = 0.05\n",
    "r = cauchy.rvs(size=n)\n",
    "epsilon = math.sqrt((1 / (2 * n)) * math.log(2 / alpha))\n",
    "\n",
    "F_n = lambda x : sum(r < x) / n\n",
    "L_n = lambda x : max(F_n(x) - epsilon, 0)\n",
    "U_n = lambda x : min(F_n(x) + epsilon, 1)\n",
    "\n",
    "xx = sorted(r)\n",
    "\n",
    "df = pd.DataFrame({\n",
    "    'x': xx, \n",
    "    'F_n': np.array(list(map(F_n, xx))), \n",
    "    'U_n': np.array(list(map(U_n, xx))), \n",
    "    'L_n': np.array(list(map(L_n, xx))), \n",
    "    'CDF': np.array(list(map(norm.cdf, xx)))\n",
    "})\n",
    "df['in_bounds'] = (df['U_n'] >= df['CDF']) & (df['CDF'] >= df['L_n'])\n",
    "\n",
    "plt.plot( 'x', 'L_n', data=df, color='red')\n",
    "plt.plot( 'x', 'U_n', data=df, color='green')\n",
    "plt.plot( 'x', 'CDF', data=df, color='purple')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "387ec21b08ed4af59f7c82c865cb8b6f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Average fraction in bounds: 0.204\n"
     ]
    }
   ],
   "source": [
    "# 1000 iterations with Cauchy distribution\n",
    "\n",
    "bounds = []\n",
    "for k in tqdm_notebook(range(1000)):\n",
    "    n = 100\n",
    "    alpha = 0.05\n",
    "    r = cauchy.rvs(size=n)\n",
    "    epsilon = math.sqrt((1 / (2 * n)) * math.log(2 / alpha))\n",
    "\n",
    "    F_n = lambda x : sum(r < x) / n\n",
    "    L_n = lambda x : max(F_n(x) - epsilon, 0)\n",
    "    U_n = lambda x : min(F_n(x) + epsilon, 1)\n",
    "\n",
    "    # xx = sorted(r)\n",
    "    xx = r # No need to sort without plotting\n",
    "\n",
    "    df = pd.DataFrame({\n",
    "        'x': xx, \n",
    "        'F_n': np.array(list(map(F_n, xx))), \n",
    "        'U_n': np.array(list(map(U_n, xx))), \n",
    "        'L_n': np.array(list(map(L_n, xx))), \n",
    "        'CDF': np.array(list(map(norm.cdf, xx)))\n",
    "    })\n",
    "    all_in_bounds = ((df['U_n'] >= df['CDF']) & (df['CDF'] >= df['L_n'])).all()\n",
    "    bounds.append(all_in_bounds)\n",
    "    \n",
    "print('Average fraction in bounds: %.3f' % np.array(bounds).mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.4**. Let $X_1, \\dots, X_n \\sim F$ and let $\\hat{F}_n(x)$ be the empirical distribution function.  For a fixed $x$, use the central limit theorem to find the limiting distribution of $\\hat{F}_n(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "We have:\n",
    "\n",
    "$$\n",
    "\\hat{F_n}(x) = \\frac{\\sum_{i=1}^n I\\left(X_i \\leq x \\right)}{n}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "I\\left(X_i \\leq x\\right) =\n",
    "    \\begin{cases}\n",
    "      1   & \\text{if } X_i \\leq x \\\\\n",
    "      0   & \\text{if } X_i > x\n",
    "    \\end{cases}       \n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Let $Y_i = I\\left(X_i \\leq x \\right)$ for some fixed $x$.  From the central limit theorem,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\sqrt{n} (\\overline{Y}_n - \\mu_Y) & \\leadsto N(0, \\sigma_Y^2) \\\\\n",
    "\\overline{Y}_n & \\leadsto N(\\mu_Y, \\sigma_Y^2 / n)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "We can estimate the mean $\\mu_Y$ as $\\mathbb{E}(\\hat{\\mu_Y}) = \\mathbb{E}(\\overline{Y}_n) = n^{-1} \\sum_{i=1}^n \\mathbb{E}(I(X_i \\leq x)) = n^{-1} \\sum_{i=1}^n F(x) = \\hat{F_n}(x)$.\n",
    "\n",
    "We can estimate the variance $\\sigma_Y^2$ as $\\mathbb{E}(\\hat{\\sigma_Y}^2) = \\mathbb{E}(\\mathbb{V}(\\overline{Y}_n)) = n^{-1} \\sum_{i=1}^n \\mathbb{E}((Y_i - \\overline{Y}_n)^2) = n^{-1} \\sum_{i=1}^n \\left( \\mathbb{E}(Y_i^2) - 2 \\mathbb{E}(Y_i \\overline{Y}_n) + \\mathbb{E}(\\overline{Y}_n^2) \\right) \\leq n^{-1} \\sum_{i=1}^n \\left( \\mathbb{E}(Y_i) + \\mathbb{E}(\\overline{Y}_n^2) \\right) \\leq 2$\n",
    "\n",
    "Therefore, for large $n$, the limiting distribution has variance that goes to 0 -- so $\\overline{Y}_n \\leadsto \\mu_Y$, or $I\\left(X_i \\leq x \\right) \\leadsto F(x)$ for every x.  Then,\n",
    "\n",
    "$\\hat{F_n}(x) \\leadsto n^{-1} \\sum_{i=1}^n F(x) = F(x)$,\n",
    "\n",
    "and, as expected, $F$ is the limiting distribution of $F_n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.5**.  Let $x$ and $y$ be two distinct points.  Find $\\text{Cov}(\\hat{F_n}(x), \\hat{F_n}(y))$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "We have:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\text{Cov}(\\hat{F_n}(x), \\hat{F_n}(y)) & = \\mathbb{E}(\\hat{F_n}(x) \\hat{F_n}(y)) - \\mathbb{E}(\\hat{F_n}(x))\\mathbb{E}(\\hat{F_n}(y)) \\\\\n",
    "&= \\mathbb{E}(\\hat{F_n}(x) \\hat{F_n}(y)) - F(x)F(y)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "But:\n",
    "\n",
    "$$\n",
    "\\hat{F_n}(x) \\hat{F_n}(y) = \\frac{1}{n^2} \\sum_{i=1}^n \\sum_{j=1}^n I(X_i \\leq x) I(X_j \\leq y)\n",
    "$$\n",
    "\n",
    "so\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{E}(\\hat{F_n}(x) \\hat{F_n}(y)) & = \\frac{1}{n^2} \\sum_{i=1}^n \\sum_{j=1}^n \\mathbb{E}(I(X_i \\leq x) I(X_j \\leq y)) \\\\\n",
    "&= \\frac{1}{n^2} \\sum_{i=1}^n \\sum_{j=1}^n \\mathbb{P}(X_i \\leq x, X_j \\leq y) \\\\\n",
    "&= \\frac{1}{n^2} \\sum_{i=1}^n \\sum_{j=1}^n \\mathbb{P}(X_i \\leq x | X_j \\leq y) \\mathbb{P}(X_j \\leq y) \\\\\n",
    "&= \\frac{1}{n^2} \\left( \\sum_{i=1}^n F(\\min\\{x, y\\}) + \\sum_{i=1}^n \\sum_{j=1, j \\neq i}^n F(x)F(y) \\right) \\\\\n",
    "&= \\frac{1}{n} F(\\min\\{x, y\\}) + \\left( 1 - \\frac{1}{n}\\right) F(x)F(y)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Therefore, assuming $x \\leq y$,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\text{Cov}(\\hat{F_n}(x), \\hat{F_n}(y)) & = \\mathbb{E}(\\hat{F_n}(x) \\hat{F_n}(y)) - \\mathbb{E}(\\hat{F_n}(x))\\mathbb{E}(\\hat{F_n}(y)) \\\\\n",
    "&= \\mathbb{E}(\\hat{F_n}(x) \\hat{F_n}(y)) - F(x)F(y) \\\\\n",
    "&= \\frac{1}{n} F(\\min\\{x, y\\}) + \\left( 1 - \\frac{1}{n}\\right) F(x)F(y) - F(x)F(y) \\\\\n",
    "&= \\frac{F(x)(1 - F(y))}{n}\n",
    "\\end{align}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.6**.  Let $X_1, \\dots, X_n \\sim F$ and let $\\hat{F}$ be the empirical distribution function.  Let $a < b$ be fixed numbers and define $\\theta = T(F) = F(b) - F(a)$. Let $\\hat{\\theta} = T(\\hat{F}_n) = \\hat{F}_n(b) - \\hat{F}_n(a)$.\n",
    "\n",
    "- Find the estimated standard error of $\\hat{\\theta}$.\n",
    "- Find an expression for an approximate $1 - \\alpha$ confidence interval for $\\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "**(a)**\n",
    "\n",
    "The estimated mean for $\\hat{\\theta}$ is $\\mathbb{E}(\\hat{\\theta}) = \\mathbb{E}(\\hat{F}_n(b) - \\hat{F}_n(a)) = \\mathbb{E}(\\hat{F}_n(b)) - \\mathbb{E}(\\hat{F}_n(a)) = F(b) - F(a) = \\theta$.\n",
    "\n",
    "The estimated variance for $\\hat{\\theta}$ is\n",
    "\n",
    "$$ \\mathbb{V}(\\hat{\\theta}) = \\mathbb{E}(\\hat{\\theta}^2) - \\mathbb{E}(\\hat{\\theta})^2 $$\n",
    "\n",
    "But\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{E}(\\hat{\\theta}^2) & = \\mathbb{E}((\\hat{F}_n(b) - \\hat{F}_n(a))^2)  \\\\\n",
    "&= \\mathbb{E}(\\hat{F}_n(a)^2 + \\hat{F}_n(b)^2 - 2 \\hat{F}_n(a)\\hat{F}_n(b)) \\\\\n",
    "&= \\mathbb{E}(\\hat{F}_n(a)^2) + \\mathbb{E}(\\hat{F}_n(b)^2) - 2 \\mathbb{E}(\\hat{F}_n(a)\\hat{F}_n(b))\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{E}(\\hat{F}_n(a)^2) & = \\mathbb{E}\\left(\\left(\\frac{1}{n} \\sum_{i=1}^n I(X_i \\leq a)\\right)^2\\right) \\\\\n",
    "& = \\frac{1}{n^2} \\left(\\sum_{i=1}^n \\mathbb{E} \\left( I(X_i \\leq a)^2 \\right) + \\sum_{i=1}^n \\sum_{j=1, j \\neq i}^n \\mathbb{E}\\left( I(X_i \\leq a) I(X_j \\leq a) \\right) \\right) \\\\\n",
    "& = \\frac{1}{n^2} \\left(\\sum_{i=1}^n \\mathbb{E} \\left( I(X_i \\leq a) \\right) + \\sum_{i=1}^n \\sum_{j=1, j \\neq i}^n \\mathbb{E}\\left( I(X_i \\leq a) I(X_j \\leq a) \\right) \\right) \\\\\n",
    "& = \\frac{1}{n^2} \\left(\\sum_{i=1}^n \\mathbb{P} \\left( X_i \\leq a \\right) + \\sum_{i=1}^n \\sum_{j=1, j \\neq i}^n \\mathbb{P}\\left( X_i \\leq a, X_j \\leq a \\right) \\right) \\\\\n",
    "& = \\frac{1}{n^2} \\left(n F(a) + n(n-1) F(a)^2 \\right) \\\\\n",
    "&= F(a) \\frac{1}{n} + F(a)^2 \\left(1 - \\frac{1}{n} \\right) \\\\\n",
    "\\mathbb{E}(\\hat{F}_n(b)^2) & = F(b) \\frac{1}{n} + F(b)^2 \\left(1 - \\frac{1}{n} \\right)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "From the previous exercise,\n",
    "\n",
    "$$ \\mathbb{E}(\\hat{F}_n(a)\\hat{F}_n(b)) = \\frac{1}{n} F(a) + \\left( 1 - \\frac{1}{n}\\right) F(a)F(b) $$\n",
    "\n",
    "Putting it together,\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "\\mathbb{E}(\\hat{\\theta}^2) \n",
    "&= \\mathbb{E}(\\hat{F}_n(a)^2) + \\mathbb{E}(\\hat{F}_n(b)^2) - 2 \\mathbb{E}(\\hat{F}_n(a)\\hat{F}_n(b)) \\\\\n",
    "&=   F(a) \\frac{1}{n} + F(a)^2 \\left(1 - \\frac{1}{n} \\right) +  F(b) \\frac{1}{n} + F(b)^2 \\left(1 - \\frac{1}{n} \\right) - 2 \\left( \\frac{1}{n} F(a) + \\left( 1 - \\frac{1}{n}\\right) F(a)F(b) \\right) \\\\\n",
    "&=  \\frac{1}{n} \\left(F(b) - F(a)\\right) + \\left( 1 - \\frac{1}{n}\\right) \\left( F(b) - F(a) \\right)^2 \\\\\n",
    "&=  \\frac{1}{n} \\theta + \\left( 1 - \\frac{1}{n} \\right) \\theta^2 \\\\\n",
    "\\mathbb{V}(\\hat{\\theta}) &=  \\mathbb{E}(\\hat{\\theta}^2) - \\mathbb{E}(\\hat{\\theta})^2 \\\\\n",
    "&=  \\frac{1}{n} \\theta + \\left( 1 - \\frac{1}{n} \\right) \\theta^2 - \\theta^2 \\\\\n",
    "&= \\frac{\\theta (1 - \\theta)}{n} \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Finally, the estimated standard error is\n",
    "\n",
    "$$\\text{se}(\\hat{\\theta}) = \\sqrt{\\mathbb{V}(\\hat{\\theta})} = \\sqrt{\\frac{\\hat{\\theta}(1 - \\hat{\\theta})}{n}}$$\n",
    "\n",
    "**(b)**\n",
    "\n",
    "An approximate $1 - \\alpha$ confidence interval is\n",
    "\n",
    "$$ \\hat{\\theta} \\pm z_{\\alpha/2}\\text{se}(\\hat{\\theta}) = \\hat{\\theta} \\pm z_{\\alpha/2} \\sqrt{\\frac{\\hat{\\theta}(1 - \\hat{\\theta})}{n}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.7**.  Data on the magnitudes of earthquakes near Fiji are available on the course website.  \n",
    "\n",
    "- Estimate the CDF.\n",
    "- Compute and plot a 95% confidence envelope for F.\n",
    "- Find an approximate 95% confidence interval for F(4.9) - F(4.3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "data = pd.read_csv('data/fijiquakes.csv', sep='\\t')\n",
    "r = np.array(data['mag'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x26bbe4be828>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3yN1x/A8c8Re8+qCqVo0aqVGkWrdkvtVat2qU0VVaWTWjVLjRg/SokdxB61R40WpTFCSok9YmSc3x8nrpvkJrmJm9zk5vt+ve4ruc9z7nPP00e/jjO+R2mtEUII4ZpSOLsCQggh4o8EeSGEcGES5IUQwoVJkBdCCBcmQV4IIVxYSmd9cc6cOXWBAgWc9fVCCJEkHT58+LrWOpe95Z0W5AsUKMChQ4ec9fVCCJEkKaX8YlNeumuEEMKFSZAXQggXJkFeCCFcmNP65G0JCgrC39+fR48eObsqDpU2bVrc3d1JlSqVs6sihEhmElWQ9/f3J1OmTBQoUACllLOr4xBaa27cuIG/vz8FCxZ0dnWEEMlMjN01SilPpdQ1pdRfUZxXSqlJSilfpdRxpVSZuFbm0aNH5MiRw2UCPIBSihw5crjcv06EEEmDPX3yc4E60Zx/HygS9uoKTHueCrlSgH/KFe9JCJE0xNhdo7XeqZQqEE2RBsB8bXIW71NKZVVK5dFaX3FQHYVI1u4/uc9+//0c+PcAD4MfOrs6iV5oiCLoUWqePExNaEgKtAatFWhl3++ADjW/E/a71gpQYb/z3L9XL/syXd+vlCD/PRzRJ58XuGT13j/sWKQgr5Tqimntkz9/fgd8tRCu5+r9q+y6uMu8Lu3iyJUjhOgQABRJ+1+FOsQNHmeGx1kgOG0MrzTh3welhycZ4UmGsJ9RvILTO/s2Y3SxxXa6vp8w3+WIIG/rT53NnUi01jOAGQAeHh6JcreSjBkzcv/+fWdXQyQTWmv+ufmPJaj/fvF3fG/6ApA2ZVrK5y3P4MqDqZK/ChXcK5AlbRan1jc4GO7cgdu37X9Zl793L+7fnS4dZMxo9coW4X2EV4YMkDIlpEgBSplXYvk9d+6qDnsmMXFEkPcH8lm9dwcuO+C6Qriky/cus+LUCrac38Kui7sICAwAIEe6HFTOX5lPyn5C5fyVKZOnDKndUsd7fe7cAV/fZ68rVyIH56evmNo/SkHWrOFfhQtHPpY5swnaadPa90qd2lxbxJ4jgvxqoKdSajFQHrjjkP74vn3h6NHnvkw4pUrBhAkOuVT79u3JnDkzhw4d4r///mP06NE0bdrUIdcWrufinYssO7mMZaeWsfvSbgBeyfYKHxT5gMr5K1M5f2Vey/FavAzSaw03b4YP5Nav69fDl8+aFbJlexaQixSJHKSjemXMaFqqIvGIMcgrpRYBVYGcSil/YDiQCkBrPR1YB3wA+AKBQIf4qmxic+XKFXbt2sXff/9N/fr1JciLcM7dOseyk8vwOuXFgX8PAFAyd0m+fe9bmhRrQrFcxRz+naGh8PffsG+feR05YgL57dvPyigF+fKZFnbjxubn09crr5huDuE67Jld81EM5zXQw2E1espBLe741LBhQ1KkSEHx4sW5evWqs6sjEoEzN87gddILr5NeHPnvCABl85RlZPWRNCnWhCI5ijj0+65fh/37zWvfPvPz7l1zLmtW8PCA1q1NAC9UyPwsWNB0gYjkIVGteE1q0qRJY/nd/F0nkiO/237MPTqXZaeW8ee1PwGo4F6BsTXH0rhYYwpmc8xK56AgOH78WSt93z7TSgdwc4MSJaBVK6hQwbyKFJGuEyFBXog4O3z5MGP3jmXpiaWE6lAq56/MxDoTaVS0Efmy5Iv5Ana4eRPWroWVK2HDBnjwwBzPk8cE8i5dzM+yZaWbRdgmQT6CwMBA3N3dLe/79+9P//79nVgjkZiE6lDW/7OesXvHsv3CdjKlzkTfCn3pXb43+bM4Zu3HxYuwapUJ7Dt2QEgIvPQStGsH771ngrq7ezKbbRIcDA8fmuk9sXkFB5uR59BQwlY8JY7fO3Y0k0sSgAT5CEJDQ+0qN3fu3HDvZW69a3sU/IgFxxcwbu84/r7+N+6Z3Rlbcyydy3R+7rnrWsOff5qgvnKlGSwFKF4cBg2Chg1NSz1JdL08eWLfRPo7d+DRo+hfjx8/+z042P46uLk9myifKlX8T3x3c4v9Z3PkiL9nEIEEeSGicT3wOtMOTmPKwSlce3CN0i+WZmHjhTQr3oxUbnFPHR0cDLt3m6C+ahWcP2/+/3/7bRgzBho0MH3qCS4kBG7ciHo1U0yvwMDor+/mFvVE+SxZop8sny4dZMoUecVTxFVQadIks3/mRE+CfAy+//57li5dGu5Ys2bNGDp0qJNqJBKC701fftr7E3OOzuFh8EM+KPIBAyoO4L0C78V5LntQEGzdCsuWwYoVZmZMmjRQsyYMHQr16kHu3A6+EVuePIELF55NlD979tnv58+bikYlZcrIk+Nfesn8zJIl5on0GTJIAE5gEuRjMHToUAnoyYjfbT8GbBzA8lPLSeWWijYl2tC/Yn9ef+H1OF3v8WPYvBm8vEyL/dYt0xitV8/MUa9TxzQ+HU5r8Pd/NlHe+uXnZ/qHn8qUycytLFnSVCpv3vCroaxf6dNLkE5iJMgLgZkCO+uPWQzYOACNZkjlIfQs15M8mfLE+loPH5qZMF5esGaNmbeeNSvUrw9Nm5qWu8PnqQcGwuHD4edXXrbKLpI9uwnkFSpAmzbhV0DlyiWB24VJkBfJ3qU7l+i8pjMbz26kWsFqzK4/mwJZC8TqGg8ewPr1JrB7e5v32bNDs2YmsFerZvKvOITW8M8/z4L5/v1w7JjpTwez6unpNBwPD3j1VVMZkSxJkBfJltYazyOe9N/Yn5DQEKZ+MJVuHt1IoeyfxnLiBIwfD4sWmRb8Cy9A27bQpAm8+66Z3OGAipog7u0Ne/aYoH7zpjmXKROULw9DhpigXq6caZkLEUaCvEiW/O/602VNF3x8fXj35XfxbODJK9leseuzWsO2bTB2rGm9p0tnAnurVlC5splA8tyCg2HXrmfzKv38TJfKG2+YfvOny1qLFnXQFwpXJUE+ggsXLlCvXj3++uvZlrYjRowgY8aMfPbZZ06smXAErTXzjs2jr09fgkKDmPz+ZD5961O7Wu9BQbB0qQnuR46YVvu330K3bpAzpwMq9+ABbNxogrq3t2mtp0kDtWrBV1+Z0doXXnDAF4nkRIK8SDYu37tM1zVdWfvPWqrkr8KcBnMolL1QjJ+7exdmzoSJE+HSJdN4njnTjF8+9wBqQIAZnV21ygT4R4/MzJZ69cwqqFq14mn6jUguEm2Q7+vTl6P/OTaffKkXSzGhjmOzW1atWpXy5cuzbds2bt++zezZs6lSpYpDv0M8H601C44voLdPbx4HP2ZC7Qn0Kt8rxtb7pUswaRLMmGECfdWqMG0avP/+c64+vXjx2WT53bvNdMb8+aFrVxPYK1d2UGe+EIk4yCclwcHBHDhwgHXr1vH111+zefNmZ1dJhLly7wrd1nZj9enVVMpXiTkN5sSY7vfoURg3DhYvNv3vzZvDgAEmtUCcnT1rAruXFxw8aI6VKAFffmkCe6lSMo1RxItEG+Qd3eK2V1SrGaNb5di4cWMAypYty4ULF+KjWiIOvM94035lex4EPWBcrXH0Kd8HtxRRD1I+eACDB8OUKaaHpFcv6NMHXn45jhU4ffpZYH+akMbDA0aNMtNvCheO44WFsF+iDfLOkiNHDm7duhXu2M2bNylYMOqc4E/zyru5uREcm0RKIl48Dn7M4M2DmbB/AqVeLMWiJosomrNotJ/ZuRM6dDCr+vv0gREjzAKmWNEaTp40Qd3LC54O3lesaEZrmzSBAgXicktCxJkE+QgyZsxInjx52LJlC9WrV+fmzZv4+PjQp08fZ1dN2OGfG//QcllL/rjyB73K9WJ0zdGkTRn16OiDB/DFF6bvvVAhk9o31kMqZ8/CnDmm1f7336bbpXJlM1LbuLHJCyyEk0iQt2H+/Pn06NGDAQMGADB8+HAKFYp5FoZwrgXHF9B9bXdSu6VmZYuVNCjaINryu3aZ1ruvr+maGTkylhtv7N1rOu+XLzeBvWpV6N0bGjWCF198rnsRwlEkyNtQvHhxtm3bZlfZ7du3W37PmTOn9Mk7wf0n9+mxrgfzj82nSv4qLGy8MNqdmQIDzXjnhAmm92TbNhOf7RISAqtXm+6XPXtMn86gQeZviZdecsTtCOFQEuRFknbkyhFaeLXg7K2zDH93OF++8yUpU0T9x3rPHmjf3qR+6dHDjIHaNQ09MBDmzoWffjJN/wIFTHdMx44yj10kahLk7dSjRw92794d7lifPn3o0KGDk2qUvGmtmXxgMgM3DSRX+lxsbbeVdwu8G2X5hw/NotFx48yU9C1bTNKwGF29ClOnws8/m800ypWDJUtMl0xK+d9HJH7yp9ROU6dOdXYVRJjrgdfpuKoja86sod6r9ZjTYA4500edV2DfPtN6P33apCAYPdrk9YrW33+bzGPz55tNNurXh88+g0qVZD67SFIkyIskZceFHbRe3pqAwAAm1plIr3K9olzDEBRkWu+jR5t9MDZuNLnco3XuHPTrZ/rd06Y1I7P9+pl0vUIkQRLkRZKgtWb07tF8sfULCmUrxL5O+yidp3SU5c+dg48+ggMHoFMn0yjPnDmaLwgNNTkLPv/cdMOMGAGffippe0WSJ0FeJHqPgx/T1bsr84/Np8XrLZj54UwypYm6v+W330waGKVM1simTWP4ggsXzADqtm1Qu7bJPpYv6tk5QiQlz5NmSYh4F/AggBr/q8H8Y/P5uurXLGqyKMoAHxgIXbpAy5bw+usmB020AV5rmD7d5JA5dMgE9/XrJcALlyIt+Qjc3NwoUaKE5f3KlSspIEvRneJkwEnq/VqPK/evsLjJYlq80SLKsn/+CS1amPHSIUPg669jSOTo5wedO5tdtmvUgNmzzbQbIVyMBPkI0qVLx9Gjjk1xLGJvg+8Gmns1J13KdOxov4NyecvZLKc1/PKLGRvNmtUMrtaoEc2FtYZZs0xayacf7tJFZswIl5Vog3zfvuaf245UqpRZ5egIc+fOZfXq1QQGBnL27FkaNWrE6NGjHXPxZG7KgSn08elDiRdKsPqj1eTPYruFfeuWic/Llpmu9PnzY9g46dIl03rfuNFMkp89WxKGCZeXaIO8szx8+JBSpUoBULBgQVasWBFl2aNHj3LkyBHSpEnDa6+9Rq9evcgn/blxFhwaTJ/1ffj50M/Uf60+CxsvJGNq26tJ9+wxs2cuX4YxY6B//2g28tAaPD1NoZAQs7Dpk0+ec+cPIZIGu4K8UqoOMBFwA2ZprUdFOJ8FWADkD7vmWK31nOepmKNa3LEVm+6a6tWrkyVLFsDku/Hz85MgH0e3H92mhVcLNp7dyGcVP2NUjVE2c7+HhMCPP5r57y+/bDZWKme7J8fw9zfNfR8fk6DG0xOiSRsthKuJMcgrpdyAqUBNwB84qJRarbU+aVWsB3BSa/2hUioXcFoptVBr/SReap1IPM0jD5JL/nmcvXmWeovq4XvTl1kfzqJTmU42y/33n9lXdcsWM4Nm+nQI+zvWtpUrzdTIx49h8mQz711a7yKZsaclXw7w1VqfA1BKLQYaANZBXgOZlFl6mBG4CUjEEzHa6beTxr81RqPZ1HYTVQtUtVnu6FH48EOTPmb2bLMQNcqx0kePTAqCqVPNTkyLFskuTCLZsqdZkxe4ZPXeP+yYtSlAMeAy8CfQR2sdGvFCSqmuSqlDSqlDAQEBcayycBVzj86lxvwa5Eifg32d9kUZ4FetMntwgOmL79gxmgB/6hSUL28C/IABpj9HArxIzrTW0b6AZph++Kfv2wKTI5RpCvwEKKAwcB7IHN11y5YtqyM6efJkpGOuwpXvLS7G7xmvGYGuPq+6vhl402aZ0FCtR4/WWimt33pL68uXo7lgaKjWs2drnT691rlyab1uXfxUXAgnAw7pGOK29cuelrw/YD2a6I5psVvrACwPq4NvWJCPflNNkWxN2j+J/hv706RYE9a3Xk+2dNkilXnyxMx2/PxzaNbMbMuXJ08UF7x7F1q3NklqKlSAY8fg/ffj9yaESCLsCfIHgSJKqYJKqdRAS2B1hDIXgeoASqncwGvAOUdW1Fk2bNhAqVKlwr0aNWrk7GolWVMPTKWPTx8aFW3EoiaLSOUWeVnqjRtQq5aZCPPVV6ZLPV26KC548CCULm1yvH//vZkDH+XfBkIkPzEOvGqtg5VSPYENmCmUnlrrE0qpbmHnpwPfAnOVUn9iumwGaa2vx6VCWusoU8c6Q+3ataldu/ZzXcP8C0v8cugXeq7vSf3X6rO46WKbAf7vv6FePTPzccEC00C3KTTUpJYcMsRsu7dzJ7z9dvzegBBJkF3z5LXW64B1EY5Nt/r9MlDreSuTNm1abty4QY4cORJVoH8eWmtu3LhB2rRpnV0Vp5p5eCbd1najbpG6LGm6hNRuqSOV2bzZJBRLk8YkhKxYMYqLXbsGH39s5r43bmzSFGSL3OUjhEhkK17d3d3x9/fH1WbepE2bFnd3d2dXw2k8j3jS1bsrdQrXwau5F2lSpolUZvp06NkTihUDb2+z0MmmzZuhbVuT02DaNLNy1UUaBELEh0QV5FOlSkVBWY3oUuYdnUfn1Z2pVagWK1qsIG3K8P+iCQ42Mx0nTYK6dU3/u82t+bSG776D4cOhaFHT926VLVQIYVuiCvLCtSw4voAOqzpQrWA1VrZYGSnA371rVq6uX2+ySI4ZA26RMxmYxU0dO5q/Adq2NS34DBkS5iaESOIkyIt4sejPRXy88mOqFqjK6o9Wky5V+Okxly+bGTSnT5tsv127RnGhq1ehYUOzG/fIkTBokHTPCBELEuSFwy05sYQ2K9pQJX8V1ny0hvSp0oc7f+UKvPeeCfQ+PlC9ehQXOn7c5DK4fh2WLweZuipErEm2JuFQy04uo9WyVryd7228W3mTIXX4bpWrV00q93//Nd00UQZ4b2+oVMl02v/+uwR4IeJIgrxwmBWnVtByWUvKu5dnXat1kXLBX7tmAvzFi7Bu3bN8NOFobea/168Pr70GBw5AmTIJcwNCuCAJ8sIh1pxeQ3Ov5pTNU5b1rddH2mw7IMC02s+fh7Vr4Z13bFwkKMhMiRwwwMx/37kT8kbMhSeEiA0J8uK5nb91no+WfUSpF0uxoc0GMqfJHO78jRtm31VfX1izxuzdEcnNm1CnDsycCV98YdIUpE9vo6AQIjZk4FU8F601XdZ0QSnFsubLyJI2/C4eN2+aAH/6tAnwNvvgz5wxuQz8/MxGrW3bJkzlhUgGJMiL5zL7yGy2nN/CtLrTIm24fesW1KwJJ0+anPA1a9q4wNatJpeBm5vZ8slmR70QIq6ku0bEmf9dfwZsHEDVAlXpWjb8RPfbt808+L/+ghUrTE9MJDNnQu3aJmvkgQMS4IWIBxLkRZxorem+tjtBIUHM+nAWKdSzP0p375qgfuwYeHnBBx/YuMCYMWYFVI0aZrsnSWchRLyQ7hoRJ7/++SveZ7wZX2s8hbIXshy/d88E+MOHTYD/8EMbHx4/3uwG0rIl/O9/kFL+GAoRX+T/LhFrV+9fpbdPbyq4V6B3+d6W4/fvm1b7gQNmckyDBjY+PGGCmSLZrJkEeCESgHTXiFjrtb4X95/cx7O+J24pTEaxBw9MFsm9e00escaNbXxwyhSTiaxJE1i4UAK8EAlA/i8TsbL81HKWnlzK99W+p1iuYoDZpKlxY9i1C3791TTSI5k2DXr1MsnGFi2CVJF3hRJCOJ4EeWG3mw9v8unaTyn1YikGvj3Qcnz6dJPefdo0aNHCxgdnzIBPPzUd9L/9JgFeiAQkQV7Yrf+G/tx4eAOfNj6W/Vn9/Ez231q1TEaCSGbPNifq1oWlSyF15G3/hBDxR/rkhV3W/7OeecfmMbjSYEq9WAowucSe5oGfMcNGmve5c6FLFzPdxsvLbN4qhEhQ0pIXMbr7+C6feH9C8VzF+fKdLy3H58413TRTptjYk/V//zO7OdWoYVZDJfONzIVwFgnyIkaDNg3i33v/sqfZHssm3JcvQ//+UKUKdO8e4QO//grt25u8wqtWSYAXwomku0ZEa/uF7Uw/PJ1+FfpR3r08YLppPv3UbL06ezaksP5TtHixSTD2zjuwejWkS2f7wkKIBCEteRGlB08e0Gl1JwpnL8w3731jOb5kiWmgjxkDRYpYfWDpUmjTxuSg8faWVMFCJAIS5EWUhm0bxrlb59j+8XbLPq0BAdCzJ5QrZ9Y1WSxfDh99BBUqmF1BMmSwfVEhRIKS7hph0z7/fUzYN4FPPT7l3QLvWo737g137oCnp8kODMDZs6YF/9ZbZuPWjBltX1QIkeAkyItIHgc/puOqjuTLko9RNUZZjq9aZbrchw2D118POxgaCp07mwVOXl6QKZPtiwohnEK6a0Qk3+38jlPXT+HT2seyV+vt22YWTcmSMHiwVeGZM2H7dvNT9mMVItGRIC/C+fPqn4zaPYp2JdtRu3Bty/EBA+DaNTOeaslKcOkSDBxo9vTr1Mk5FRZCREu6a4RFSGgIXdZ0IWvarIyrNc5yfONG0wf/+edQpkzYQa1NuoKQENOKj7TcVQiRGNgV5JVSdZRSp5VSvkqpwVGUqaqUOqqUOqGU2uHYaoqEMO3QNPb/u58JtSeQM31OwGwC0qULFC0KX31lVfh//zODrKNGya5OQiRiMXbXKKXcgKlATcAfOKiUWq21PmlVJivwM1BHa31RKfVCfFVYxI9Ldy4xZMsQaheqTasSrSzHhwwxvTK7d1stXP3vP+jbFypVgh49nFNhIYRd7GnJlwN8tdbntNZPgMVAxD1/WgHLtdYXAbTW1xxbTRGftNb0XN+TUB3KtLrTUGFdLzt3wtSp0KcPVKxoKWyWuwYG2ljuKoRIbOz5PzQvcMnqvX/YMWuvAtmUUtuVUoeVUu1sXUgp1VUpdUgpdSggICBuNRYOt+zUMlafXs03Vb+hYDbT9RIYaMZSX3kFvvvOqrCXl0k49s038NprzqmwEMJu9syusTWipm1cpyxQHUgH7FVK7dNanwn3Ia1nADMAPDw8Il5DOMHtR7fptb4XZfKUoU+FPpbjw4eDry9s3Wq1ePX6dbPctWxZk51MCJHo2RPk/YF8Vu/dgcs2ylzXWj8AHiildgIlgTOIRG3QpkFce3CNta3WkjKF+eNw4ACMH28mz7z3nlXhvn3h1i3YvFn2ZxUiibCnu+YgUEQpVVAplRpoCayOUGYVUEUplVIplR4oD5xybFWFo+3028mMP2bQr0I/yuQxcyODgswC1pdegtGjrQp7e5vNt4cOhRIlnFNhIUSsxdgc01oHK6V6AhsAN8BTa31CKdUt7Px0rfUppZQPcBwIBWZprf+Kz4qL5/M4+DFd13SlQNYCfF31a8vxMWPgzz9NluDMmcMO3r5tmvUlSpjpNkKIJMOuf3NrrdcB6yIcmx7h/RhgjOOqJuLTD7//wOkbp/Fp7UOG1KbT/Z9/zHhqs2Zmz22LgQPNtMlVq2SPViGSGJn/lgyduHaCkbtG0rpEa0vqgqcLWNOmhYkTrQpv3gyzZplA7+HhnAoLIeJMRs+SmVAdSlfvrmRKk4mfav9kOT5nDmzbZjbkzpMn7OD9+2a566uvmuk2QogkR4J8MvPLoV/Yc2kPcxvMJVeGXABcvQqffWb2aw2XZ+yLL8DPD37/XbbxEyKJku6aZOTfu/8yeMtgqhesTruSz9ar9e0LDx6YVrxlAeuuXTBlipkXX6mScyoshHhuEuSTkV7re/Ek5AnT6023pC5Yt85sBDJ0qElCBsDDh6ZJ//LL8MMPzquwEOK5SXdNMrHi1ApW/L2CUdVHUTh7YcB0uXfvDsWLR9gI5Jtv4MwZM+gqW/kJkaRJkE8G7jy6Q8/1PSmZuyT9Kz5LRzBsGFy8aDJMWmZGHjtmJst37Gg2AxFCJGkS5JOBIVuG8N/9/1jZYiWp3My2TgcPwqRJpiX/9tthBUNCzGyaHDlMoBdCJHkS5F3c3kt7mXZoGn3L9+WtvG8BJnVBly7w4oswcqRV4SlTTPRfvBiyZ3dOhYUQDiVB3oWFhIbQY10P8mbKyzfvfWM5Pn686ZVZvhyyZAk76OdnRl8/+ACaN3dOhYUQDidB3oX9cvgXjvx3hN+a/kamNJkAOHsWRoyARo3MC3i2EQjAzz/Lfq1CuBAJ8i4q4EEAQ7cOpVrBajQr3gx4lrogdWqYPNmq8JIlZi7lTz+ZaZNCCJchQd5FDd48mPtP7jPl/SmWOfHz58OWLaaxnvfp3l43b0Lv3vDWW9Crl/MqLISIFxLkXdDeS3vxPOrJwLcHUixXMQACAsxmTm+/bVrzFgMHwo0bsHEjuLk5p8JCiHgjK15djPVg67B3hlmO9+sH9+7BzJlWqQu2bQNPT5O4pmRJ51RYCBGvpCXvYp4Oti5ustgy2Lphg9nU6auvzOpWwKQu+OQTKFRIMkwK4cIkyLsQ68HW5q+baZAPHkC3bvDaaxE2dfr+e7NLyKZNkmFSCBcmQd6F2Bps/f57uHABduwwG4IAZn+/H3+Ejz+GGjWcVl8hRPyTPnkXsc9/H55HPelXoZ9lsPXMGRg7Ftq1g3feCSv4NHVB1qwwbpzzKiyESBDSkncBtgZbtTYzItOlg9GjrQpPmwb798OCBSZHjRDCpUmQdwEzDs/gjyt/hBtsXbHCzIqcOBFy5w4reOmS6ZivXRtatXJehYUQCUZprZ3yxR4eHvrQoUNO+W5XEvAggFenvEqZPGXY3HYzSikePIBixSBbNjh8GFKmxDTtGzQwq6H++gsKFnR21YUQcaCUOqy19rC3vLTkk7ghW4ZEGmz94QfTaF+4MCzAAyxbBmvWmE56CfBCJBsy8JqE7fPfx+wjsyMNto4ZA23bmo25Abh1y3TQlykDffo4r8JCiAQnLfkk6ulg60uZXop5sHXwYJPXYO1aq6a9ECI5kP/jk6joBlsnTDAbggBmb78ZM0zqgjJlnFdhIYRTyMBrEjPw9NUAABY6SURBVBTwIIDXprxG6Tylww22Fi9uNgH544+wBntwMJQta7prTp2CDBmcXXUhxHOSgddkYMiWIdx7co/J708ON9h68SLs3GnVIzNtGhw/Dl5eEuCFSKZk4DWJeTrY2rd8X4rnMtnGnq5sDTfYevUqfPkl1KoFjRs7r8JCCKeSIJ+EBIcGWwZbv3r3K8AMtvbubfLShBtsHTTIZJqcPFm28xMiGZPumiRk0v5J/HHlj3B7tq5caVIJhxts3bUL5s0zq1tffdV5FRZCOJ1dA69KqTrARMANmKW1HhVFubeAfUALrbVXdNeUgdfYuXD7Aq///DrVClZjdcvVKKUIDDQrW2WwVYjkw+EDr0opN2AqUBPwBw4qpVZrrU/aKPcjsCF2VRYx0VrTfW13UqgUTP1gqgy2CiHsZk+ffDnAV2t9Tmv9BFgMNLBRrhewDLjmwPoJYPFfi/Hx9eH7at+TP0t+wOz3MWYMtGljY7C1Zk0ZbBVCAPYF+bzAJav3/mHHLJRSeYFGwPToLqSU6qqUOqSUOhQQEBDbuiZLNwJv0MenD+XylqPHWz2AZytb06Y1gd5CBluFEBHYE+RtRYuIHfkTgEFa65DoLqS1nqG19tBae+TKlcveOiZrn236jFuPbjHzw5m4pXADng22fv11hJWt8+bBgAFmrz8hhMC+2TX+QD6r9+7A5QhlPIDFYX3FOYEPlFLBWuuVDqllMrX1/FbmHp3LkMpDeDP3mwAEBkLfvlCiBPTsGVYwOBh69IB8+Ux3jRBChLEnyB8EiiilCgL/Ai2BcDtOaK0tuWuVUnMBbwnwz+dh0EO6rulK4eyFLQnI4Nlg644dEQZbjx2DpUtlsFUIEU6MQV5rHayU6omZNeMGeGqtTyiluoWdj7YfXsTNNzu+4eyts2xpt4V0qdIBZkbkmDHQurXVnq1Xr8KwYWawtUkT51VYCJEo2bUYSmu9DlgX4ZjN4K61bv/81Urejl89zpg9Y+hQqgPVClYDIDTU7L+dIQOMH29VeNAg04cjg61CCBtkxWsiExIaQufVncmeLjtjaj6bOjNzphlb9fSEF14IO/h0sHXwYBlsFULYJEE+kZl6cCoHLx/k18a/kiN9DgAuX4bPP4dq1aB9+7CCTwdb3d1lsFUIESUJ8onIxTsX+WLLF9QpXIeWb7S0HO/VC548gV9+seqRmT5dBluFEDGSIJ9IaK35dO2naDTT6k6zpC5YuRKWLzezagoXDiv8dGVrjRoy2CqEiJYE+URi6cmlrP1nLeNrjadA1gIA3L1r5sK/+abZvc9i8GAZbBVC2EWCfCJw6+Eteq/vTdk8ZelVvpfl+BdfmP745cshVaqwg7t3w9y5ZlZN0aJOqa8QIumQIJ8IfL7pc64HXmd96/WkTGEeyd698PPPpj++XLmwgkFBMtgqhIgVCfJOtuPCDmYdmcXAtwdSOk9pwAyyduliYvl331kVnjDBDLZ6eUHGjM6psBAiSZEg70SPgh/R1bsrr2R7hRFVR1iOjx4NJ07AmjWQKVPYwXPnYPhwaNBA0ggLIewmQd6JRmwfwZkbZ9jYZiPpU6UH4PRp+PZbaN4c6tULK6g1dOtmktVMmSKDrUIIu0mQd5L1/6znx90/0rl0Z2oWqgmY1AVdu0L69DBxolXhBQtg0yYT4N3dnVNhIUSSJEHeCS7duUTbFW15M/ebTHp/kuW4p6fZym/mTKs88QEB0K8fVKwI3bs7p8JCiCTLnk1DhAMFhQTRcllLHoc8ZmmzpZYMk//9BwMHwrvvQqdOVh8YMMBMmJ8xA1LI4xJCxI605BPY0K1D2XNpD4uaLOLVHK9ajvfpY3bumzHDqst940b43//MdMk33nBOhYUQSZo0DROQ9xlvxuwZQ7ey3cLlpvH2hiVLTCx/9WncDww0g62vvgpDhzqnwkKIJE9a8gnE77Yf7Va0o/SLpfmpzk+W4/fuma72N94wmSYtRoyA8+dh+3azY7cQQsSBBPkE8CTkCc29mhOiQ1jSbAlpUz4L2l9+Cf/+a1ryqVOHHTxyxOwM0rmz6aQXQog4kiCfAAZtGsSBfw/g1cyLwtkLW47v329yjH36qZk8A5g88V26QM6cZlWUEEI8Bwny8WzFqRVM2D+BXuV60aT4s7TAd+9CmzaQN69JI2wxaRIcPgy//QbZsiV8hYUQLkWCfDw6d+scHVZ14K2X3gq3lZ/WprF+/jxs2waZM4eduHDBbMpdrx40a+aUOgshXIsE+XjyOPgxzZc2RynFb01/I03KNJZz06aZPviRI6FKlbCDWpsR2BQpYOpUSV0ghHAICfLxZMDGARy+cpiVLVZSMFtBy/HDh80C1g8+iDCbZtEi8PEx+Qzy50/4CgshXJLMk48HS04sYerBqfSv0J8GRRtYjt++bXphcueG+fOtFrDeuAF9+5rE8T16OKfSQgiXJC15B/vnxj90Xt2ZCu4VGFVjlOW41tCxI1y6ZPLT5Mhh9aHPPoNbt2DzZnBzS/hKCyFclrTkHehh0EOaezUnlVsqfmv6G6ncUlnOTZoEK1bAqFFW0yUBtm412/kNHGg2cxVCCAeSlrwD9fXpy9H/juL9kTf5szzrVz9wwMTw+vWhf3+rDzx8aHILFy5sZtUIIYSDSZB3kN/++o0Zf8xgUKVB1H21ruX4zZtmA5CXXjIN9nCTZoYNg7NnYcsWSJcuwesshHB9EuQd4OKdi3zi/Qnl85bn2/e+tRzXGtq3h8uXYdeuCGubvL1h3DiThKxatQSvsxAieZAg/5xCQkNos7wNITqEhY0XhuuHHzfO7NM6caKZOGPh5wft2kGpUvDTT5EvKoQQDiJB/jn9uPtHfr/4O/MazqNQ9kKW43v2wODB0KQJ9Opl9YEnT6BFC5OjZulSyTAphIhXds2uUUrVUUqdVkr5KqUG2zjfWil1POy1RylV0vFVTXwO/HuA4duH0+L1FrR9s63l+PXrJo6//DLMnh2hH37wYJOZzNPTDLgKIUQ8irElr5RyA6YCNQF/4KBSarXW+qRVsfPAu1rrW0qp94EZQPn4qHBicf/JfVota8VLmV5ier3pqLBIHhpqemKuXYO9eyFLFqsPrVxpumd69oSmTZ1TcSFEsmJPd005wFdrfQ5AKbUYaABYgrzWeo9V+X2AuyMrmRj1Xt+b87fPs+3jbWRNm9VyfPRoWL/epJ8pU8bqA+fPm1FYDw8YOzbB6yuESJ7s6a7JC1yyeu8fdiwqnYD1tk4opboqpQ4ppQ4FBATYX8tEZumJpcw5OochlYfwzsvvWI7v3Gl26mvRwuQas3j82MyjBJOZLE0ahBAiIdjTkreVDlHbLKjUe5ggX9nWea31DExXDh4eHjavkdhdunOJrt5dKZe3HMPfHW45HhAAH30EhQpF2IwbTNqCQ4fMkteCBSNfVAgh4ok9Qd4fyGf13h24HLGQUupNYBbwvtb6hmOql7iEhIbQdkVbgkKCwk2XfJqX5sYNWLfOKj88mBk0U6aY1JMNGzqn4kKIZMueIH8QKKKUKgj8C7QEWlkXUErlB5YDbbXWZxxey0RizJ4x7PDbgWd9z3Db+E2bZtY2TZwIJa3nFfn6QqdOUL68SVojhBAJLMYgr7UOVkr1BDYAboCn1vqEUqpb2PnpwFdADuDnsFkmwVprj/irdsI7dPkQw7YNo1nxZrQv1d5y/ORJGDAA6tSJMB/+0SOTVzhlSrOVn2WXbiGESDhKa+d0jXt4eOhDhw455btj6/6T+5T5pQwPgx9yvNtxsqUz+QkePzaN9MuX4fhxePFFqw917w7Tp5slr/XqOafiQgiXo5Q6HJtGtKx4tUM/n3743vRl68dbLQEeYMgQOHbMxPFwAX7RIhPgP/9cArwQwqkkn3wMlp1cxqwjsxhUaRBVC1S1HN+40axr6tEjQhw/fdqkD65UCb77LsHrK4QQ1qS7Jhr+d/15c9qbFMpeiN0dd5PazfSrBwSY/T2yZzczIy1ZggMDoUIFuHIFjhwBd5dfEyaESGDSXeMgoTqUj1d+zOOQx/za+FdLgNcaOnc2eeJ9fCKkge/dG/780yx5lQAvhEgEJMhHYdyecWw9v5VZH86iSI4iluO//AKrV8P48RGmS/76q8lG9sUXZqqNEEIkAtJdY8Ok/ZPot6EfjYo2YmmzpZbkY6dOQdmyUKWKaayneDqicf68yQ3/xhuwY4eZNimEEPFAumueQ3BoMP18+jHl4BQaFm3IvIbzLAH+8WNo1QoyZDDb+FkCfHAwtA1LM7xwoQR4IUSiIhEpzL3H92i5rCXr/lnHgIoD+LHGj7ilcLOc//JLOHoUVq2CPHmsPvjDD7B7NyxYAAUKJHi9hRAiOhLkMbNo6v1aj7+u/cW0utPo5tEt3PnNm0124O7doX59qxN798I330Dr1uYlhBCJTLLvk//jyh98uOhD7j2+x9JmS6lduHa48zduQIkSkDWrmS6ZPn3Yibt3TT+81qaJH253ECGEiB/SJx8La06v4aNlH5E9XXZ2d9xNidwlwp1/Ol3y+nWTXdIS4MEkqvHzM0nkJcALIRKpZLviddL+STT8rSHFchVjf+f9kQI8wKxZZse+kSNNo91i8WKYPx+GDTMrW4UQIpFKdt01EWfQLGi0gAypM0Qqd/q02b7v7bdhwwar2TR+fmaCfPHiphUvs2mEEAlIumuiEdMMmqcCAswWfunSwbx5VgE+JATatDG7dS9YIAFeCJHoJZsoFdMMmqf+/hvq1jXpg1esgJdesjo5ciTs2mW6al55JWEqLoQQzyFZBPnDlw9Tf3F97j2+h3crb+oUtp12YPt2aNTI7O+xfbvJFW+xfz+MGGE2cm3TJgFqLYQQz8+lB15DdSjj9oyj4uyKuCk3dnfcHWWAnzcPatUyLff9+yME+Hv3zHJXd3f4+ecIu3QLIUTi5bIt+Sv3rvDxyo/ZdG4TDYs2ZNaHs8iRPkekclrDV1+Z1O/Vq4OXl5kTH07v3nDhgslLE+mkEEIkXi4Z5NecXkPH1R158OQB0+tOp2vZrpYcNNYePYIOHcyMyE6dzIbcqVJFKLRkiUlWM2wYVK6cIPUXQghHcakg/zDoIQM3DWTqwamUzF2SRU0WUSxXMZtlAwKgYUPYswdGjTI79UX6e+DiRfjkE9N3M2xY/N+AEEI4mMsE+eNXj9NqWStOBJygX4V+jKw+kjQp09gse/o0fPCBmUGzdCk0bWqjUEiIyS4ZHGyyS0Zq4gshROKX5IO81prJBybz+abPyZo2Kz6tfSLln7G2fTs0bmxidqQZNNZ+/NEsdpo7FwoVioeaCyFE/EvSQf7ag2t0WNWBdf+so26Rung28OSFDC9EWX7+fJOLpnBhWLsWCha0UUhrGDfO5BZu0QLatYu/GxBCiHiWZIO8j68PH6/8mDuP7jD5/cn0eKuHzcFVsHMGDUBQEPTsCTNmQLNmMGeOTJcUQiRpSS7IPwp+xJDNQ5iwfwJvvPAGW9pt4Y0X3oiy/M2b0KULLF8ezQwagDt3TGDftAmGDDF/I6Rw6WUEQohkIMkF+YXHFzJh/wR6levFjzV+JF2qdFGW3b7djJ1evWp6YPr1i6Jh7udnchmcPm1ST3bqFG/1F0KIhJTkgnyH0h0omrMolfJHneI3KMhkIBg5EooUMRs4lS0bReGDB+HDD82keR8f058jhBAuIsn1R6RQKaIN8L6+Zs3SDz+YBvkff0QT4Jcvh3ffNekm9+yRAC+EcDlJLshHRWuTf6Z0aThzxsx/nzkTMkROFW8Kjx1rJsi/+Sbs22fywwshhItxiSB/+7ZJDtm+vdno4/jxKBY4genL6d4dBg40hbZtg9y5E7K6QgiRYJJ8kN+1y2zU5OUF338PW7dCvnxRFL5zB+rVg19+gcGDTdKadFEP3AohRFJnV5BXStVRSp1WSvkqpQbbOK+UUpPCzh9XSpVxfFXDCw6G4cNNl3rKlLB7N3zxBbhF3ujJ8PMz+7Fu3Wpm0IwcKVMkhRAuL8bZNUopN2AqUBPwBw4qpVZrrU9aFXsfKBL2Kg9MC/sZL86fh9atzayZdu00k8c8JnOK+3D+Hty/H/51757pz/nhBzODZv16qFEjvqomhBCJij1TKMsBvlrrcwBKqcVAA8A6yDcA5muzK/g+pVRWpVQerfUVR1fY59uDtBhRFLTm13T9+GjhPJgfEvMHX3nFtOJlgFUIkYzYE+TzApes3vsTuZVuq0xeIFyQV0p1BboC5M+fP7Z1BaBw0ZRUfOEs06v8SoE8GSHTYMiYMepXpkzmZ44csvG2ECLZsSfq2VojquNQBq31DGAGgIeHR6Tz9ijcrDQ+zQBKxeXjQgiRrNgz8ugPWM9XcQcux6GMEEKIBGZPkD8IFFFKFVRKpQZaAqsjlFkNtAubZVMBuBMf/fFCCCFiJ8buGq11sFKqJ7ABcAM8tdYnlFLdws5PB9YBHwC+QCDQIf6qLIQQwl52jURqrddhArn1selWv2ugh2OrJoQQ4nnJaiAhhHBhEuSFEMKFSZAXQggXJkFeCCFcmDJjpk74YqUCAL84fjwncN2B1UlqkvP9J+d7h+R9/3Lvxsta61z2ftBpQf55KKUOaa09nF0PZ0nO95+c7x2S9/3Lvcft3qW7RgghXJgEeSGEcGFJNcjPcHYFnCw5339yvndI3vcv9x4HSbJPXgghhH2SakteCCGEHSTICyGEC0v0QV4p5aaUOqKU8rZxLsE3EE9IMdx7VaXUHaXU0bDXV86oY3xRSl1QSv0Zdm+HbJx32Wdvx727+rPPqpTyUkr9rZQ6pZSqGOG8Kz/7mO491s8+KeyH1wc4BWS2cS5BNxB3gujuHeB3rXW9BKxPQntPax3V4hdXf/bR3Tu49rOfCPhorZuG7WGRPsJ5V372Md07xPLZJ+qWvFLKHagLzIqiiGUDca31PiCrUipPglUwHtlx78mdyz775EwplRl4B5gNoLV+orW+HaGYSz57O+891hJ1kAcmAJ8DoVGcj2oDcVcQ070DVFRKHVNKrVdKvZ5A9UooGtiolDoctgF8RK787GO6d3DdZ/8KEADMCeuqnKWUyhChjKs+e3vuHWL57BNtkFdK1QOuaa0PR1fMxrEkPyfUznv/A5PDoiQwGViZIJVLOJW01mUw/zTvoZR6J8J5l3z2YWK6d1d+9imBMsA0rXVp4AEwOEIZV3329tx7rJ99og3yQCWgvlLqArAYqKaUWhChjKtuIB7jvWut72qt74f9vg5IpZTKmeA1jSda68thP68BK4ByEYq46rOP8d5d/Nn7A/5a6/1h770wgS9iGVd89jHee1yefaIN8lrrIVprd611Aczm4Vu11m0iFHPJDcTtuXel1ItKKRX2eznMs7yR4JWNB0qpDEqpTE9/B2oBf0Uo5pLP3p57d+Vnr7X+D7iklHot7FB14GSEYi757O2597g8+6QwuyYclYw3EI9w702B7kqpYOAh0FK7zvLl3MCKsD/LKYFftdY+yeTZ23PvrvzsAXoBC8Nml5wDOiSTZw8x33usn72kNRBCCBeWaLtrhBBCPD8J8kII4cIkyAshhAuTIC+EEC5MgrwQQrgwCfJCCOHCJMgLIYQL+z/qJWbgg9KrdwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = len(r)\n",
    "alpha = 0.05\n",
    "epsilon = math.sqrt((1 / (2 * n)) * math.log(2 / alpha))\n",
    "\n",
    "F_n = lambda x : sum(r < x) / n\n",
    "L_n = lambda x : max(F_n(x) - epsilon, 0)\n",
    "U_n = lambda x : min(F_n(x) + epsilon, 1)\n",
    "\n",
    "xx = sorted(r)\n",
    "\n",
    "df = pd.DataFrame({\n",
    "    'x': xx, \n",
    "    'F_n': np.array(list(map(F_n, xx))), \n",
    "    'U_n': np.array(list(map(U_n, xx))), \n",
    "    'L_n': np.array(list(map(L_n, xx)))\n",
    "})\n",
    "\n",
    "plt.plot( 'x', 'L_n', data=df, color='red')\n",
    "plt.plot( 'x', 'U_n', data=df, color='green')\n",
    "plt.plot( 'x', 'F_n', data=df, color='blue')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% confidence interval: (0.526, 0.588)\n"
     ]
    }
   ],
   "source": [
    "# Now to find the confidence interval, using the result from 8.5.6:\n",
    "\n",
    "import math\n",
    "from scipy.stats import norm\n",
    "\n",
    "z_95 = norm.ppf(.975)\n",
    "theta = F_n(4.9) - F_n(4.3)\n",
    "se = math.sqrt(theta * (1 - theta) / n)\n",
    "\n",
    "print('95%% confidence interval: (%.3f, %.3f)' % ((theta - z_95 * se), (theta + z_95 * se)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.8**.  Get the data on eruption times and waiting times between eruptions of the old faithful geyser from the course website.\n",
    "\n",
    "- Estimate the mean waiting time and give a standard error for the estimate.\n",
    "- Also, give a 90% confidence interval for the mean waiting time.\n",
    "- Now estimate the median waiting time.\n",
    "\n",
    "In the next chapter we will see how to get the standard error for the median."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "data = pd.read_csv('data/geysers.csv', sep=',')\n",
    "r = np.array(data['waiting'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated mean: 70.897\n",
      "Estimated SE: 13.570\n"
     ]
    }
   ],
   "source": [
    "# Estimate the mean waiting time and give a standard error for the estimate.\n",
    "theta = r.mean()\n",
    "se = r.std()\n",
    "\n",
    "print(\"Estimated mean: %.3f\" % theta)\n",
    "print(\"Estimated SE: %.3f\" % se)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90% confidence interval: (48.576, 93.218)\n"
     ]
    }
   ],
   "source": [
    "# Also, give a 90% confidence interval for the mean waiting time.\n",
    "\n",
    "import math\n",
    "from scipy.stats import norm\n",
    "\n",
    "z_90 = norm.ppf(.95)\n",
    "\n",
    "print('90%% confidence interval: (%.3f, %.3f)' % ((theta - z_90 * se), (theta + z_90 * se)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated median time: 76.000\n"
     ]
    }
   ],
   "source": [
    "# Now estimate the median time\n",
    "median = np.median(r)\n",
    "\n",
    "print(\"Estimated median time: %.3f\" % median)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.9**.  100 people are given a standard antibiotic to treat an infection and another 100 are given a new antibiotic.  In the first group, 90 people recover; in the second group, 85 people recover.  Let $p_1$ be the probability of recovery under the standard treatment, and let $p_2$ be the probability of recovery under the new treatment.  We are interested in estimating $\\theta = p_1 - p_2$.  Provide an estimate, standard error, an 80% confidence interval and a 95% confidence interval for $\\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**. Let $X_1, \\dots, X_100$ be indicator random variables (0 or 1) determining recovery on the first group, and $Y_1, \\dots, Y_100$ indicating recovery on the second group. From the problem formulation, we can assume $X_i \\sim \\text{Bernoulli}(p_1)$ and $Y_i \\sim \\text{Bernoulli}(p_2)$.\n",
    "\n",
    "If $\\theta = p_1 - p_2$, then from exercise 8.5.2:\n",
    "\n",
    "$$\\hat{\\theta} = \\hat{p_1} - \\hat{p_2}$$\n",
    "\n",
    "$$\\text{se}(\\hat{\\theta}) = \\sqrt{\\hat{p_1}(1 - \\hat{p_1}) + \\hat{p_2}(1 - \\hat{p_2})}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated mean: 0.050\n",
      "Estimated SE: 0.466\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "p_hat_1 = 0.9\n",
    "p_hat_2 = 0.85\n",
    "\n",
    "theta_hat = p_hat_1 - p_hat_2\n",
    "se_theta_hat = math.sqrt(p_hat_1 * (1 - p_hat_1) + p_hat_2 * (1 - p_hat_2))\n",
    "\n",
    "print('Estimated mean: %.3f' % theta_hat)\n",
    "print('Estimated SE: %.3f'   % se_theta_hat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "80% confidence interval: (-0.548, 0.648)\n",
      "95% confidence interval: (-0.864, 0.964)\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "z_80 = norm.ppf(.9)\n",
    "z_95 = norm.ppf(.975)\n",
    "\n",
    "print('80%% confidence interval: (%.3f, %.3f)' % ((theta_hat - z_80 * se_theta_hat), (theta_hat + z_80 * se_theta_hat)))\n",
    "print('95%% confidence interval: (%.3f, %.3f)' % ((theta_hat - z_95 * se_theta_hat), (theta_hat + z_95 * se_theta_hat)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.10**. In 1975, an experiment was conducted to see if cloud seeding produced rainfall.  26 clouds were seeded with silver nitrate and 26 were not.  The decision to seed or not was made at random.  Get the data from the provided link.\n",
    "\n",
    "Let $\\theta$ be the difference in the median precipitation from the two groups.\n",
    "\n",
    "- Estimate $\\theta$.\n",
    "- Estimate the standard error of the estimate and produce a 95% confidence interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "data = pd.read_csv('data/cloud_seeding.csv', sep=',')\n",
    "X = data['Seeded_Clouds']\n",
    "Y = data['Unseeded_Clouds']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated mean: 177.400\n"
     ]
    }
   ],
   "source": [
    "theta_hat = X.median() - Y.median()\n",
    "\n",
    "print('Estimated mean: %.3f' % theta_hat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4e34a3bf799e4c7d9726a3008f898365",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=10000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Estimated SE: 62.363\n"
     ]
    }
   ],
   "source": [
    "# Using bootstrap (from chapter 9):\n",
    "\n",
    "nx = len(X)\n",
    "ny = len(Y)\n",
    "\n",
    "B = 10000\n",
    "t_boot = np.zeros(B)\n",
    "for i in tqdm_notebook(range(B)):\n",
    "    xx = X.sample(n=nx, replace=True)\n",
    "    yy = Y.sample(n=ny, replace=True)\n",
    "    t_boot[i] = xx.median() - yy.median()\n",
    "    \n",
    "se = np.array(t_boot).std()\n",
    "\n",
    "print('Estimated SE: %.3f' % se)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% confidence interval (Normal): \t 55.171, 299.629\n",
      "95% confidence interval (percentile): \t 37.800, 262.560\n",
      "95% confidence interval (pivotal): \t 92.240, 317.000\n"
     ]
    }
   ],
   "source": [
    "# See example 9.5, page 135\n",
    "\n",
    "from scipy.stats import norm\n",
    "\n",
    "z_95 = norm.ppf(.975)\n",
    "\n",
    "normal_conf = (theta_hat - z_95 * se, theta_hat + z_95 * se)\n",
    "percentile_conf = (np.quantile(t_boot, .025), np.quantile(t_boot, .975))\n",
    "pivotal_conf = (2*theta_hat - np.quantile(t_boot, 0.975), 2*theta_hat - np.quantile(t_boot, .025))\n",
    "\n",
    "print('95%% confidence interval (Normal): \\t %.3f, %.3f' % normal_conf)\n",
    "print('95%% confidence interval (percentile): \\t %.3f, %.3f' % percentile_conf)\n",
    "print('95%% confidence interval (pivotal): \\t %.3f, %.3f' % pivotal_conf)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}