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 deviation of $\\text{Bernoulli}(p)$, so the plugin estimator for the standard error is $\\sqrt{\\frac{\\hat{p}(1-\\hat{p})}{n}} = \\sqrt{\\frac{\\overline{X}_n(1 - \\overline{X}_n)}{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{\\frac{\\overline{X}_n(1-\\overline{X}_n)}{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{\\frac{\\hat{p}(1 - \\hat{p})}{n} + \\frac{\\hat{q}(1 - \\hat{q})}{n}} = \\sqrt{\\frac{\\overline{X}_n(1 - \\overline{X}_n)}{n} + \\frac{\\overline{Y}_m(1 - \\overline{Y}_m)}{n}} $$\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{\\frac{\\overline{X}_n(1 - \\overline{X}_n)}{n} + \\frac{\\overline{Y}_m(1 - \\overline{Y}_m)}{n}} $"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd1zV1R/H8dcRQVScuEUcuLeJe+TeI1ea20qzLJtm69ful78ss7IyK3/unSs3mntDLnABThREBAVkXu75/fE1fw5SknvvF7if5+PBAy73yzkfr/L2e8/3fM9RWmuEEEJkf7nMLkAIIYRtSKALIUQOIYEuhBA5hAS6EELkEBLoQgiRQ+Q2q+NixYrpChUqmNW9EEJkSwEBAVFa6+LpPWdaoFeoUAF/f3+zuhdCiGxJKXX+756TIRchhMghJNCFECKHkEAXQogcwrQx9PSkpqYSFhZGUlKS2aXYhbu7O15eXri6uppdihAiB8pSgR4WFkaBAgWoUKECSimzy7EprTXXrl0jLCyMihUrml2OECIHeuiQi1JqplIqUikV+DfPK6XUt0qpEKXUUaXUY49aTFJSEp6enjkuzAGUUnh6eubYdx9CCPNlZAx9FtDlAc93Barc+hgD/JiZgnJimP8lJ//ZhBDme+iQi9Z6h1KqwgMO6Q3M0cY6vPuUUoWVUqW11uE2qlEIkQVZrBb8L/uz4/wO4lPizS4nW2np3ZJOPp1s3q4txtDLAhfveBx263v3BbpSagzGWTze3t426FoI4WixybG8sPYF1gav5XrSdQAU8u7zn5jYYmKWDfT0/ibT3TVDaz0DmAHg6+ubJXfW8PDwID5ezjaESE9KWgr9lvRj27ltDKs7jM4+nWlfqT3F8hUzuzSBbQI9DCh3x2Mv4LIN2hVCZCEBlwP4cPuHbD6zmVm9ZzGi/gizSxL3sEWgrwZeVEotApoAN2wyfv7KK3D4cKabuUv9+jB1qk2aGjlyJAULFsTf35+IiAi++OIL+vfvb5O2hcgqElITWBS4iOn+0zl4+SD5XPMxpdMUCfNHpLUm7nIcuVxy4VHKw+btPzTQlVILgTZAMaVUGPAB4HqruOnAOqAbEAIkAKNsXmUWFR4ezq5duzh58iS9evWSQBc5htaaL/d8yb93/ZvrSdepWbwm33X9jmF1h1HIvZDZ5WV51jQrMaExRAZGEhkYSdTJKK6dusa109dIiU+h5dstaf/v9jbvNyOzXJ56yPMaGGeziv5iozNpe3riiSfIlSsXNWvW5MqVK2aXI4RNWLWV1ze+ztT9U+lepTtvtniTVt6tZNrt30iJTyH8UDjhAeFEHI4gMjCSq8evYkm0GAcoKFy+MJ7VPCnXshzFqhXDu6V9JoVkqTtFs5s8efLc/tr4f02I7O/VDa/y7YFvebnJy0zpPIVcSpZ8+ovWmqvHr3J2y1ku+18mPCCcqyeu3p4G4lHKgxJ1SuD7vC8lapegZJ2SFK9ZHNd8jlnuQwJdCHFbUGQQ3x34jud9n+frzl/LWTlgSbZwfvt5Tq85zek1p7l+1piq6VHKgzK+Zaj5ZE3K+JahTMMydhkX/yck0O+RkJCAl5fX7cevvfYar732mokVCeE4H27/kPxu+fmk7SdOHeY3I28SvC6Y07+fJnRTKCnxKeR2z02lDpVoMbEFVbpWoZB31ruWIIF+D6vVmqHjZs2adddjmbsusrv1wetZdnwZH7X5CM98nmaX41BaayKPRXLq91MErwkmbH8YaChQpgC1B9emWs9qVGxX0WFDJ49KAl0IJ6e1Zn3Iesb8PoaaxWsyscVEs0tyCEuShbNbz3J6zWmC1wRz48INAMr4lqHNh22o2qMqpRqUylbvVCTQH+Kzzz5j6dKld31vwIABvPvuuyZVJITt/HH2D9774z32hu2lYuGKzO0zlzy58zz8B7OpuPA4gtcGc3rNac74nSE1IRXXfK5U6liJ1v9qTZXuVShQuoDZZT4yCfSHePfddyW8RY6jtebJZU+y7PgyvAp6Mb37dEY1GIWbi5vZpdmU1pqIQxHGBc3fT3PZ37iJvWC5gtQbWY+qPapSsW1FcrvnjCjMGX8KIcQ/sv38dpYdX8aE5hP4uO3HuOd2N7skm0lNSOXsH2dvj4fHXY4DBV5NvGj7aVuq9axGiTolstVQSkZJoAvhhL7Z/w3F8xXnozYf5YgwT7mZQuCiQE6tPMWZzWewJFlw83DDp5MPVXtWpXLXyniUNHdKoSNIoAvhZLTW7Di/g341+pHXNa/Z5WRKXHgcB6YdIGB6AInRiRSuUJjHRj9G1R5VKf94eXLnca6Ic64/rRCCy3GXiU6Mpl7JemaX8sgijkSw7+t9HFtwDKvFSvUnqtPstWaUa1EuRw6lZJQE+j3OnTtHjx49CAz8/xaqH374IR4eHrzxxhsmViaEbQRHBwNQrVg1kyv5Z7RVE7IhhL1T9nJ2y1lc87viO9aXJi83oahPUbPLyxIk0IVwMn6hfgCUL1Te5EoyJjUxlaPzjrLv631EnYiiQJkCtJ/UnoZjGpK3SPYeMrK1LBvor2x4hcMRtl0PvX6p+kztYttVHNu0aUOTJk3YunUr169f59dff6VVq1Y27UMIW0hNS2X8+vFMD5jOE9WfoHLRymaX9EA3I29y8IeDHPzhIAlXEyjVoBR95vWh1oBauLi5mF1elpRlAz07sVgsHDhwgHXr1vHRRx+xefNms0sS4j7PrH6GuUfnMrHFRD5r91mWHWuOOhXFni/3cHTuUdKS06jasyrNXmtG+cfLZ9mas4osG+i2PpPOqL/7B/Ogf0h9+/YFoGHDhpw7d84eZQmRKYmpiSwJWsJzDZ9jUodJZpeTrtiwWLZ9tI3D/z2Mi6sL9UfVp+krTSlWTfYrzagsG+hm8fT0JCYm5q7vRUdHU7Fixb/9mb/WRXdxccFisdi1PiH+iZspN1kbvJa5R+eSnJZMr2q9zC7pPgnXEtj1+S4OTDsAGhq/2JhW77Qif4n8ZpeW7Uig38PDw4PSpUuzZcsW2rdvT3R0NBs2bODll182uzQhMuRmyk3WBa9j6fGlrDm9hkRLIiXzl+SNZm/QyaeT2eXdlhKfwt6v97L3y72kxKdQb3g9Hv/wcQqXL2x2admWBHo65syZw7hx43j99dcB+OCDD/Dx8TG5KiEeTGvNDwd/4K0tbxGfEk+J/CUYVX8UA2oNoJV3K1xyZY0LidY0K4f/e5g/3vuDm1duUv2J6rT9tC0lapUwu7RsTwI9HTVr1mTr1q0ZOnbbtm23vy5WrJiMoQtTXI67zNOrnmZj6EY6+3RmYouJtC7fOsuE+F8u7r3I+hfXE/5nOOWal2PQykF4NfV6+A+KDJFAFyKbS0hNwHeGL9eTrvN9t+953vf5LDcbJP5KPFve2sLhWYcpUKYAfRf0pfag2lmuzuxOAj2Dxo0bx+7du+/63ssvv8yoUaNMqkgIQ8DlAMLjw1k2YBn9avYzu5y7WC1WDnx/gG3vbyM1MZUWb7Wg9butcfPIWcv0ZhUS6Bn0/fffm12CEPfRWjPv6DwAWni3MLmau106cIk1z60h4nAEPp196PptVzyrOtfWdvcJCoJJk6BXLxgwwObNS6ALkU1prZm4eSIz/pzB+MbjKeVRyuySAGP2ypZ3t3DguwMUKF2AAcsGUKNvDeceXjlwAD7/HFauhHz5wNfXLt1IoAuRTc07Oo/Jeybzgu8Lpt2Id68zW86w+pnV3Dh/A98XfGn/7/a4F8r+660/kvh4I8BnzoStW6FIEXj/fXjpJShmn5ulJNCFyKbWBK+hXMFyTOs2zfSz3+TYZPze9CPgpwCKVinKqJ2j8G7pbWpNpkhNBT8/mD/fCPOEBPD2hsmT4bnnoIB99yuVQE9HREQEr7zyCgcPHiRPnjxUqFCBqVOnUq9ePapXr05SUhIFChRg3LhxjBgxAoBZs2YxYcIEypYtC0DdunWZM2eOmX8MkQNprTlw6QAzD81k9anV9K/Z3/QwP/vHWVaNWsWNizdo9noz2n7SFte8rqbW5FBWK+zZA4sWwZIlcPWqcTY+bBgMGQItWkCuXA4pRQL9Hlpr+vTpw4gRI1i0aBEAhw8f5sqVK/j4+HDo0CEAzpw5Q9++fbFarbdnugwcOJBp06aZVrvIuSJvRjLv6DxmHppJ0NUg8rnmY2CtgXza7lPTarIkWdj89mb2T92PZ1VPnt79NOWalTOtHofSGgICjBBfvBjCwsDdHXr2NEK8a1dwc/xMHgn0e2zduhVXV1fGjh17+3v169e/74ahSpUqMWXKFF5//XWZuijs6nrSdWr9UIuohCiaejVlRo8ZDKw9kIJ5CppWU8SRCJYPWc7VoKs0erERHf/TEdd8OeCsPCkJrl17+Mfx43DmDLi6Qpcu8J//GGFu5yGVh8mygb7hlQ1EHI6waZul6peiy9QuDzwmMDCQhg0bZqi9xx57jJMnT95+vHjxYnbt2gXIHHVhO9vPbScqIYqVA1fSu3pvU2vRVs2+b/ax5a0t5C2alyHrh1C5SxZcV91qhRs3MhbOd34kJPx9m3nzgqen8VGzJrzzDvTpA0Wzzm5JWTbQswOt9V2PZchF2FJiaiKzDs/iiz1f4OHmQdcqXU2tJyEqgRXDVxCyPoRqvarR69de5CuWz6RiEiAkBE6fhuBg4yMkBCIjjWCOjjZCPT1KGSH8Vzh7eUG9ev9/nN5H0aJGoGdxGQp0pVQX4BvABfhFaz3pnucLAfMA71ttfqm1/m9mCnvYmbS91KpVi2XLlmXo2EOHDlGjRg07VySczbWEa3x/8Hu+O/AdUQlRNC7bmF96/oKbi3l3V17YfYFlA5eRcDWBrtO60uiFRva/GJucbAxr3BnawcHG40uX7j62VCmoXBnq1n1wMHt6QuHCDrtI6WgPDXSllAvwPdARCAMOKqVWa62P33HYOOC41rqnUqo4cEopNV9rnWKXqu2oXbt2vPPOO/z888+MHj0agIMHD5Jwz1uxc+fO8cYbb/DSSy+ZUabIwdrNacfRK0fpXqU7b7Z4k1berUybyaKtmj1f7mHLO1soXKEwz+x7htINStuns/h42LABli+Hffvg/Pm7z7I9PaFKFWjf3vhcpQpUrWoEuclj11lFRs7QGwMhWuszAEqpRUBv4M5A10ABZfyr8wCigWy504NSihUrVvDKK68wadIk3N3db09bDA0NpUGDBrenLb700ksyTi5sKvJmJEevHOWzdp/xTqt3TK0l4VoCK4evJHhdMDUH1KTXL73IUzCPbTuJjobff4cVK2DjRuOiZLFi0K6dMe3vr+CuUiVLjVVnVRkJ9LLAxTsehwFN7jlmGrAauAwUAAZqre8bwFJKjQHGAHh7Z92bDsqUKcOSJUvu+35iYuLf/szIkSMZOXKkHasSzuDP8D8BaObVzNQ6Lh28xNL+S4mPiLfPEIvWMHYs/PorpKUZ49ijR0PfvtCyJeSWy3uPIiOvWnp/i/qex52Bw0A7wAfwU0rt1FrH3vVDWs8AZgD4+vre24YQTu/ApQMoFA3LZGymlT38+cufrBu3Do9SHjy9+2nK+JaxfSfTpsGMGfDsszBmjLG2iTOv9WIjGQn0MODOuwW8MM7E7zQKmKSNaR8hSqmzQHXggE2qFMIJRCVEMe/oPGqXqG3KHPO0lDTWv7yegOkBVOpYiX4L+tl+FsvFi0aQf/EF9OhhfC1BbjMZCfSDQBWlVEXgEjAIGHzPMReA9sBOpVRJoBpw5lEK0lqbfiuzvdw7zVGIv8Qmx9J1flcu3LjAxqEbHd7/zcibLOm3hAu7LtD8zea0/3d7crnYaCaI1Wqsb/Ljj8Z4udZGmM+cKWFuYw8NdK21RSn1IrARY9riTK11kFJq7K3npwOfALOUUscwhmgmaq2j/mkx7u7uXLt2DU9PzxwX6lprrl27hru7k648Jx7o4+0fcyj8ECsGruDxCo87tO+IIxEs6rWIm1dv0m9RP2oPrG27xk+cMMbFT56E4sXhzTeNIZaKFW3Xh7hNmXXW6Ovrq/39/e/6XmpqKmFhYSQlJZlSk725u7vj5eWFq2sOuEVa2FSN72tQoXAF1g9Z79B+T6w4wYqhK3Av4s6gVYMo09CG4+V79xpn4q6uMGUK9OsHeWw8S8YJKaUCtNbpLqiepS4lu7q6UlH+5xZO5nLcZU5GneTZBs86rE+tNTv/vZOt722lbOOyDFw5kAKlbTiXe+9eY7542bLGdMRKlWzXtvhbWSrQhXBG/peNd6rNyjlmqqIl2cLvo3/n6Nyj1Blch56/9LTtcrdXrxrbq5UuDbt3Q4kStmtbPJAEuhAm231hNy7KhXol69m9r8ToRBb3Wcz5Hedp83EbWr/X2rbXq44cMaYiXr1qnKVLmDuUBLoQJtJas+LkCtpVbEd+t/x27SvmTAzzu83n+tnr9F3QlzpP1bFh4zHG9mo//GDc0blwITz2mO3aFxkigS6EiTaf2UxwdDBvNH/Drv1c9r/Mgu4LsFqsDNs8jPKtytuu8chIqF8frlyBF16Ajz82duwRDieBLoRJki3JvLj+RXyK+DC83nC79RO6KZTFfReTv3h+hmwYQrFqNt6geMoUiIgwxsubmbtkgbOTQBfCBKlpqQxZPoTT106zfsh63HPb5/6Eo/OPsmrkKorXKs6Q9UNsO5NFa2MLtu++g4EDJcyzgJy5KLAQWVhKWgoDlg7gtxO/MaXTFLpUts/a/3un7GXF0BV4t/Rm5PaRtg3zsDDo1QsGD4batY0t2ITp5AxdCAebtGsSq06t4ruu3/Fi4xdt3r62avwm+rH3y73U7F+TPvP6kDuPDX/VDx2Ctm0hJcUYbhk/HlxcbNe+eGQS6EI4UHhcOF/s/oJ+NfrZJczTUtNY/cxqjs49SqNxjejyTRfbrckCcO4cdOsGhQrBli3G5hIiy5BAF8JBtNaMWzeOVGsq/+lg+yGK1MRUlg5YSvDaYNp+0pZW79pwp6OLF+Hnn+Gnn4wzcwnzLEkCXQgH+fXQr6w4uYIvO36JT1Efm7adHJvMwp4LOb/zPN1/7I7v2HSX+vhnrFbYtMlYJXHNGuMiaNeuxrTEmjUz376wOQl0IRzkq71f0dSrKa82e9Wm7SZEJTCvyzyuHLlC3/k2umHoxg3joueOHcbdnhMnGqskVqiQ+baF3UigC+EAWmtuJN2gRbkW5FK2G9OOj4hnToc5RIdEM3DFQKr2qJr5RiMjoUsXCAw0hlhGjgQ3t8y3K+xOAl0IO7JqKytPruSTHZ8QHh9O5aK2G3e+ceEGc9rPIS48jiHrhlCxnQ1WKtXaODM/edLYjKJz58y3KRxGAl0IO9kQsoEJfhMIjAykStEqzOo9iyF1h9ik7ZizMcxuO5uk60kM2zSMcs3LPfyHMmL1ati/39hNSMI825FAF8IOUtNS6bekH6U9SjOvzzwG1h5I7ly2+XWLDolmdrvZpMSnMHzLcNttSpGaCu+9Z8xeGTbMNm0Kh5JAF8IOTkadJCE1gY/afGSzs3KAa6evMbvtbCzJFkZsHUGpeqVs1jaTJxvj5qtWQW6JhuxI/taEsIPjV48DUKek7ZaovXriKnPazcGaZmXE1hGUrFPSNg2HhcGnn8Ivv8CTTxpj6CJbkkAXwg7CYsMAKFfQNmPbkYGRzGk/BxSM3DaS4jWLZ77RiAj4/HNjJovVCs89B599lvl2hWkk0IWwg6sJV3HN5Uph98KZbuvK0SvMaT8HFzcXhv8x3DbL3yYlQaNGEB5uTEv817+gvA3XSBemkEAXwg4SUhPI75Y/07fe/3Vm7pLHhRFbR+BZxdM2Bc6bZwy1rF9vzDkXOYIsnyuEHYTHh1M0b9FMtXH1+FVmt5uNi5uNwxyM8fJ69WRqYg4jgS6EjSWkJrD7wm6aejV95DaiTkYxu91scrnkYvgfw20b5hYLHD4MnTqBLTeIFqaTQBfCRrTWLA5cTPVp1QmPD6d3td6P1E7UqShmt50NwIitI2y/ZVxICCQnGxtTiBxFxtCFyKTY5Fj8Qv2Yun8quy7son6p+szrO4/W5Vv/47aiQ6KZ024O2qqNMK9u4zAHOHvW+CzL3+Y4EuhC/ENaa05fO83a4LWsDV7LzvM7SbWmUiJ/CWb0mMHTDZ7GJdc/38Hn+rnrzG43m7SUNEZsHWGbqYn3unHD2GVIKfCx7RK+wnwS6EJkUJIliff+eI+VJ1cSGhMKQK3itXi16at0r9qd5uWaP/Lt/bFhsbdv5x/xxwhK1C5hy9INISHQs6fx+aefoKSNbkwSWYYEuhAZFHA5gK/2fkXr8q15rdlrdKvSjQqFK2S63fiIeOa0n0PitUSGbR5Gqfo2vJ3/Tk89ZSyN6+cHbdrYpw9hKgl0ITIoLiUOgEntJ9GsXDObtJkQlcCcDnOIvRTL0I1DKduorE3avU9EBPj7G3eCSpjnWBma5aKU6qKUOqWUClFKvfU3x7RRSh1WSgUppbbbtkwhzBcUGQRA2YK2Cd2k60nM6zyPmNAYnlr9FN4tvG3S7v87SDIW2ho6FKpWNcbNu3WzbR8iS3noGbpSygX4HugIhAEHlVKrtdbH7zimMPAD0EVrfUEpZYcBQCHMcyPpBgsCF1CvZD28C2U+eFNuprCg+wKuHLvCoJWDbLM5BRgbVKxdC4sXG2EeFwdFi8LAgcaSuPXr26YfkSVlZMilMRCitT4DoJRaBPQGjt9xzGBgudb6AoDWOtLWhQphhmRLMj8c/IFPd35KdGI0M3rMyHSbliQLi3ovImxfGP2X9KdKtyo2qBS4fBlGjTI2di5a1Fg58cknoW1bcHW1TR8iS8tIoJcFLt7xOAxocs8xVQFXpdQ2oADwjdZ6zr0NKaXGAGMAvL1t/PZSCBvbF7aPQcsGcf7GeTpU6sB/OvyHx0o/lqk201LTWPrkUs5uOcsTs5+gZr+atil2+XIYPRoSE+H7742vJcSdTkYCPb17g3U67TQE2gN5gb1KqX1a69N3/ZDWM4AZAL6+vve2IUSW8umOT0m0JLJp6CY6+nTMdHvWNCsrR6zk9O+n6fZ9N+oNr2eDKoErV4wz8QYNYP58Y7xcOKWMBHoYcOeizl7A5XSOidJa3wRuKqV2APWA0wiRDcWnxLPrwi4G1hpokzDXWrNu3DoCFwbS/vP2NHqhkQ2qvOXAAUhLg6+/ljB3chmZ5XIQqKKUqqiUcgMGAavvOWYV0EoplVsplQ9jSOaEbUsVwjGuxF+hzaw2xKfEM6j2IJu0ueWdLQT8FECLt1rQ8q2WNmkTgOvX4ZNPIG9e4wxdOLWHnqFrrS1KqReBjYALMFNrHaSUGnvr+ela6xNKqQ3AUcAK/KK1DrRn4ULYQ1RCFM1nNic8LpyVg1bStmLbTLe56z+72D1pN77P+9L+3+1tUOUtN24YKyYePgy//Qb589uubZEtZejGIq31OmDdPd+bfs/jycBk25UmhOP5hfpxJuYMG4dupJNPp0y35/+TP1ve2kLtp2rTbVq3TG94cZfvvjNuFlq1yrilXzg9WT5XiFuSLcn4nfEDoKV35odFji08xtrn11KlexWemP0EKpcNw9xiMaYn1q8vYS5uk0AXTi8lLYWf/H+iyndV+O/h/zKg5gDyuebLVJvB64JZOXwl5VuXZ8DSAbi4/vPVF9NvOBjeecfY/3PnTujRwzbtihxB1nIRTu1Q+CH6LO7D+RvnaebVjJm9Z9K+YubGuS/svsCS/ksoWbckT61+Cte8NpgPvmIFTJ0KO3ZArlzQtStMmyZn5+IuEujCqb268VWSLEmsH7Kezj6dMz3GHXEkggXdF1DIuxBDNgwhT8E8mS/y7Fno29dYv/zzz2H4cChTJvPtihxHAl04rd0XdrP9/Hamdp5Kl8pdMt1edGg08zrPI0+BPAzbNIz8xW0062TNGuPz779DjRq2aVPkSDKGLpzWkqAluOd2Z3TD0ZluK+5yHHM7zsVqsTJ001AKeRfKXINaw5Yt0KULjB9vBHn16pmuU+RsEujCaW0+u5nW5Vtn+gJoYkwi8zrPI+FqAkPWD6F4jUxuHbd5M/j6QocOxhzzzz6D3buN5W+FeAAZchFOSWtN8LVgelbN3EXF1IRUFvZYyLXT1xi8bnDmN6iwWmHAAChcGGbMMJa8dXfPXJvCaUigC6eTkJrAosBFpFpT8czr+cjtpKWmsaT/ktvL4FZqXynzxYWGGrfzT54Mzz6b+faEU5FAF04j4HIAv/z5CwsCFxCbHEuVolUe+W5QbdWsGrWKkPUh9JjRw3bL4Pr7G58bNrRNe8KpSKCLHO9myk26zO/Crgu7yJs7LwNqDeDZBs/S0rvlI01T1Fqz4dUNHJt/jHb/bkfD0TYM3/37jYW26tSxXZvCaUigixxvy9kt7Lqwi4/afMT4JuMp7F44U+3t/GwnB749QNNXm9p25cSdO41x87ZtIbf8aop/Tv7ViBxNa82y48twz+3Omy3exD135i4wHvzxIFv/tZW6w+rS6ctOtlts6+hR4zZ+b2/4739t06ZwOjJtUeRYWmve3/o+c4/OZWzDsZkO86ClQawbt44q3avQ69detl1s6/PPjbNyPz8oIXusi0cjZ+gix0lNS2VT6CZmHp7J8hPLeabBM3zV+atMtRnqF8ryIcsp17wcA5bYcLEtMKYq+vlB9+5QrtzDjxfib0igixxBa83ByweZd3QeiwIXcTXhKp55PXm31bt83PZjcqlHfzN66cAlFvdZTPEaxRm8ZjCu+Wy8+fKJE3DtGrRrZ9t2hdORQBfZ3smokzy59EmORR4jj0seelXrxbC6w+hcuTNuLm6ZavvqiavM7zYfj5IeDNkwBPfCdrjJJyDA+NzIhvuMCqckgS6ytW3nttFncR/cXNz4uefP9K/ZP9OzWP5y4+IN5nWaR67cuRi6aSgFShewSbv32bLFuDNU1moRmSSBLrIdrTV7Lu5h1uFZzD4ymyqeVVg7eC0VClewWR8JUQnM6zSP5NhkRm4fSVGfojZr+y4pKbB6NfTqBS42HJcXTkkCXWQbl2IvMefIHGYdmcXpa6fJ75qfEfVGMLnTZJudlQMkxyUzv9t8rp+7ztCNQylVv5TN2r7Pt98at/oPHTcGG7oAABzmSURBVGq/PoTTkEAX2cKMgBk8v/Z5rNpK6/KteavFWwyoNQAPNw+b9mNJtrC4z2LC/wxn4IqBlG9d3qbt3+XsWfjwQ2P+eceO9utHOA0JdJEtfLv/W+qVrMeSAUuoXLSyXfqwpllZMWwFZ7ecpfes3lTrWc0u/QBGmP91R+i339qvH+FU5MYikeUFRQYRdDWI4fWG2y3MtdasfX4tx5cep9NXnag/or5d+gHg4kVo3RpiY40LohUr2q8v4VTkDF1keRP8JlAwT0EG1xlstz62vLOFP3/+k5bvtKTZa83s1g8Ar78O0dHGphX17fgfh3A6coYusqyE1ATe9HuT9SHr+eDxDyiR3z63xO/5cg+7J+2m4XMNafepnW/u2b0bli6FCRMkzIXNyRm6yJI2hmzk+bXPc/b6WZ6u/zQvNX7JLv0cmnkIvwl+1BpYi27fd7PdYlvpiY+HkSON2/snTLBfP8JpSaCLLENrjd8ZP77a+xWbQjdRzbMaW0dspU2FNnbp78TyE/w++nd8OvvQZ04fcrnY+Q3rhAnGjkTbtkH+/PbtSzglCXRhuiRLEguOLWDK3ikEXQ2ilEcpJrWfxCtNXyFP7jx26TPUL5TfnvqNsk3K8uRvT+Li5oCbejZsgL59jQuiQtiBBLowXYc5Hdh9cTd1S9ZlVu9ZDKo9yG5BDnBx70UWP7GYYtWLMXjtYNzyZ269lwwJDYWYGPB89D1MhXgYCXRhqpS0FHZf3M34xuOZ2mWqfcewgYgjESzotoACZQowdNNQ8hbJa9f+OH0aPvsM5s835pz36GHf/oRTk1kuwlRhsWEA1C5R2+5hfi34GvM6zcPNw41hm4fhUdK2d5ne5403oEYNY1bLyy8bNxP17GnfPoVTy1CgK6W6KKVOKaVClFJvPeC4RkqpNKVUf9uVKHKybee2AdDUq6ld+7lx4QZzO8xFa82wzcMoXN52a7+kKzYWvvoKeveGc+eMr0uXtm+fwuk9dMhFKeUCfA90BMKAg0qp1Vrr4+kc9x9goz0KFTlHQmoCm89sZvWp1aw4uYKyBcpSu0Rtu/V3M/ImczvOJelGEiO2jqBYtWJ26+u2o0eNz08/LVvKCYfJyBh6YyBEa30GQCm1COgNHL/nuJeA3wBZpV/cJyI+gjWn17D61Gr8zviRZEmiYJ6CdK3clVebvmq34ZbEmETmdprLjYs3GOY3jNINHHSWHBJifK5RwzH9CUHGAr0scPGOx2FAkzsPUEqVBfoA7XhAoCulxgBjALy9vf9prSIbSkxN5JMdnzB5z2QsVgvlC5Vn9GOj6VWtF63Lt870jkIPkhyXzPyu84k6EcVTvz+FdwsH/pu7dMn4XKaM4/oUTi8jgZ7eqZO+5/FUYKLWOu1BZ1pa6xnADABfX9972xA5zPZz2xn9+2iCo4MZUW8ErzV7jTol6tj94idAakIqC3su5LL/ZZ787Ul8OvnYvc+7xMRA3rzGhxAOkpFADwPu3IrcC7h8zzG+wKJbv6jFgG5KKYvWeqVNqhTZyvnr53l7y9ssDFxIxcIV2TR0Ex19HLfetyXJWNP8/I7z9J3fl+q9TdjaLTFRwlw4XEYC/SBQRSlVEbgEDALuWvZOa317/U+l1CxgjYS584lLjuPzXZ8zZe8UlFK82+pd3mn1Dvlc8zmshrSUNJY+uZTQTaH0mtmLOk/VcVjfdzl1CsrbcXMMIdLx0EDXWluUUi9izF5xAWZqrYOUUmNvPT/dzjWKLM5itTDz0Ez+tfVfRN6MZGjdofy73b8pV6jcw3/YhqwWK8uHLOf076fp9n03Goxq4ND+b9uzB/buhVGjzOlfOK0M3SmqtV4HrLvne+kGudZ6ZObLEtnFptBNvL7pdQIjA2np3ZI1T62hUVnHT3TSVs2qUas4vszYoKLRCyZMtrp4ESZOhIULoWxZGDPG8TUIpya3/otHtvP8TjrP60ylIpVYNmAZfWv0dcgFz3tpq2bN2DUcnXeUtp+2tf8GFXe6fh02boS1a2HZMtAa/vUvI9hlRUXhYBLo4pGtOrUKNxc3Dj93mAJ5CphSg9aadS+t48+f/6TVu61o/a6dVzLUGo4fNwJ87Vpjw4q0NChaFAYNgvffhwoV7FuDEH9DAl08svUh62np3dLUMN/wygb8f/Cn+YTmtP2krX07jIgwlr4NDjYe16tnnIl37w5NmoCLA5bgFeIBJNDFIzkccZjjV4/zgu8LpvSvtWbT65s48O0BmrzShA7/6WD/4Z6TJ40wf+cdeP558PKyb39C/EMS6OKhLFYLgZGB7Avbx/5L+9kftp+TUSdxc3Gjf03Hr8OmtWbjaxvZP3U/jcc3pvOUzo4Zu4+NNT737i1hLrIkCXRxn2sJ19h6biv7w/az/9J+AsIDSEhNAMAzrydNvJowqPYgelbtSUmPkg6tTWvNxlc3sv+b/TR5uQmdv3ZQmAPs2AGurlCtmmP6E+IfkkAXAKRZ0/A748evh35l1clVpFpTcXNxo0GpBjzb4FmaeDWhSdkmVCpSyZSZLHBPmL/SxHFn5lrD+vUwdy60bw+FCtm/TyEegQS6kwuLDeMn/5+YdWQWYbFheOb1ZFyjcQysPZAGpRrYdSu4f+KvC6B/jZk7LMwPHIA334Tt28HHBz7+2P59CvGIJNCd3IClAzhw6QCdfTrzdeev6Vm1Z5YJ8b9ordnw8gYOfHeApq82pdNXnRwT5tOnGxc/ixeHadNg9Ghwc8D+o0I8Igl0JxeTGEP/mv1Z3H+x2aWkS1s1a19YS8BPATR9rSmdvnRQmAPMnm1MTdy5EwqYMzVTiH9C9hR1UgmpCby1+S1CokMomd+xFzYzyppmZfUzqwn4KYCWb7d0bJgHBIC/P3TrJmEusg05Q3dCG0M2MnbtWM5dP8fT9Z/mozYfmV3SfdJS01gxbAVBi4N4/MPHefz9xx0X5hs3Qr9+xnoszz3nmD6FsAEJdCeitebTHZ/y/rb3qeZZjW0jtvF4hcfNLus+lmQLywYu49SqU3Sc3JHmbzR3XOeBgdCjB9SuDevWycbOIluRQHcSSZYknln9DAuOLWBo3aH83PNn3HO7m13WfVJuprCk7xJCN4XSdVpXGo9r7NgC9u0Di8VYaEvCXGQzEuhOIC45jt6LerP13FY+a/cZb7d827S55A+SGJPIgu4LuLT/Er3/25v6I+s7toB9++Cnn4ybh2RzCpENSaDncBduXKD/kv78Gf4nc/vMZWjdoWaXlK74iHjmdZ5H1MkoBiwdQI2+NRzX+cGD8MEHxs1DxYoZoZ5bfjVE9iP/anOghNQEVpxYwawjs9hyZgtuLm6sGLiCntV6ml1auq6fu87cjnOJC49j8NrBVOpQyXGdz54NI0cay99+/jm8+CJ4eDiufyFsSAI9BzkVdYrJeyazJGgJcSlxVCxckfcff5+R9UdSoXAFs8tLV2RQJPM6zyP1ZirDNw/Hq6kDF71KSTE2o2jcGPz8oGBBx/UthB1IoOcgg34bRPC1YJ6s9SQj6o2gVflW5FJZ91aD8zvOs6j3InLnzc3IHSMpWcfB8+GXLjW2jZsxQ8Jc5AgS6DlEbHIsx64c4+2Wb/NJu0/MLuehgpYGsWLoCopUKsKQ9UMoXKGw44tYuxZKloTOnR3ftxB2kHVP30SGXIm/wvtb36fSN5VI02m0q9jO7JIeat/UfSwbuIwyjcrw9O6nzQlzgG3bjNUTs+CMHyEehZyhZ1PB14KZvGcyc47MISUthd7Ve/NGszdo4d3C7NL+lrZqNk3YxL4p+6jRtwZ95vXBNa+rOcVYrcaWcpUrm9O/EHYggZ7NnIo6xac7P2XBsQW4ubgxsv5IXmv2GlU9q5pd2gNZki2sHLGSoMVBNH6pMZ2/7kwuFxPeIKalGbf2//yzsc55sWKOr0EIO5FAzybC48KZ4DeBhYELyeOSh1ebvsqE5hMcvmPQo0iMSWRJ3yWc23aODl90oPkbzR1/Y9P58zBzpvERFmYsifvGGzA0a87LF+JRSKBnE8+sfoat57byerPXeaP5G5TIX8LskjIkOiSaBT0WEHMmhj7z+lB3SF3HFxEWBrVqQUKCcQF06lTo2VPWNhc5jgR6FheTGMPS40tZH7KeLzp8wYQWE8wuKcPO7zzP4j6LQcPwzcMp39qk2+k3bYKbN2HPHmjWzJwahHAACfQsJjoxmp3nd7Lt3Da2nd/GkYgjaDQ1itXgpSYvmV1ehh2edZg1z62hcIXCDF47mKKVi5pTSGws/PCDMT2xaVNzahDCQSTQTRaXHMe2c9v44+wfdwW4e253mpdrzkdtPqJNhTY0Lts4y20Nlx6rxYrfm37s+3ofFdtXZMCSAeQtmtecYuLioGtXOHLEuIlIpieKHE4C3cEsVgv+l/3ZFLoJvzN+7Avbh8VqIW/uvDQv15yP235MmwptaFSmUbYI8DslXU9i2aBlhG4MpfFLjen0VSdcXF0cW4TVCrt2wbx5sGQJxMfD4sXwxBOOrUMIE0igO8iGkA3MCJjBH2f/4EbyDRSKhmUaMqH5BDpW6kjzcs2zXYDfKTIoksV9FnP93HV6zOhBw9ENHVvA6dPGQlvz5xszWvLnh759jU2eZdxcOIkMBbpSqgvwDeAC/KK1nnTP80OAibcexgPPa62P2LLQ7G7076NJtiQzoOYAOvp0pH3F9njm8zS7LJsIWhrEqlGryFMgDyP+GIF3S2/HFvDzz0Zwaw2dOsFnnxln5PnzO7YOIUz20EBXSrkA3wMdgTDgoFJqtdb6+B2HnQUe11rHKKW6AjOAJvYoODu6FHuJsNgwvunyDeObjDe7HJuxWqxseWcLeybvoVzzcgxYOoACZRy4obLFAu++C198AV26GHPMZZch4cQycobeGAjRWp8BUEotAnoDtwNda73njuP3AQ5cAzXrC4wMBKBBqQYmV2I78RHx/Db4N85tPYfv8750mdoFFzcHjZefOAH//S/MnWvcvv/88/Dtt7IphXB6GfkNKAtcvONxGA8++34GWJ/eE0qpMcAYAG9vB78tN4HWmtCYUNacXgNAuULlTK7INoLXB7Nq5CqS45LpPas39Uc4YKu4xESYM8cI8v37jfDu3h2efdb4LDNYhMhQoKf3m6LTPVCpthiB3jK957XWMzCGY/D19U23jewszZrG0StH2XlhJ7su7GLnhZ1ExEcAUKlIJcoWKGtyhZljSbaw5e0t7Pt6HyXqlGD4H8MpUctBd6y+/jr8+CPUrg1ffQVDhhhzy4UQt2Uk0MOAO08tvYDL9x6klKoL/AJ01Vpfs0152cPCYwuZfWQ2ey7uIS4lDgDvQt60r9ielt4taeXdihrFa2TpzSYeJupUFL899RsRhyJo9GIjOk3uRG53Bw1xpKXBsmUwYIAxBVHOxoVIV0Z+Iw8CVZRSFYFLwCBg8J0HKKW8geXAMK31aZtXmYV9t/87xm8YT1XPqgypM4RW5VvR0rsl3oVyxpCS1ppDMw+xYfwGcufNzaDVg6jWs5pji1i9Gq5ehX79JMyFeICHBrrW2qKUehHYiDFtcabWOkgpNfbW89OB9wFP4Idbq+hZtNa+9ivbXNGJ0Ry4dICNIRuZun8qfar3YXH/xbi6mLS2t50kXU9izXNrCFoSRMV2Fekzt49jZ7EAJCUZqyLWrGnMKxdC/K0MvWfWWq8D1t3zvel3fP0s8KxtS8saLFYLgZGB7Avbd/vj1LVTAORSuRhUexCzn5id48L8wu4LLB+ynLhLcbT/vD3NJzR3/Prl+/bBmDFw5oyxibNrznqNhbA1mef1N+Yemcsvh37B/7I/CakJABTPV5xm5Zoxot4Imno1xbeMLwXyOPiM1c7SUtLY8dkOdn66k8IVCvP07qcp29jBF3OvXYN//QumT4cyZWDlSujQwbE1CJENSaCnY8f5HYxYOYIaxWvwbINnaerVlKZeTalQuILjN2ZwoMsBl1k1ahWRxyKpO6wu3aZ1I09BBy5HEBxsrFU+a5Yx1DJ+PHzyCRTIWf9pCmEvEuj3iLwZybAVw/Ap6sP+Z/fj4eZhdkl2Z0mysO2jbeyZvAePkh6Ov/D555/w8cfGxU9XV2MXoddeMzalEEJkmAT6HS7FXqLD3A5cvXmV7SO3O0WYX9x7kdVPrybqZBT1n65P5686417Y3bFF9OgBycnw3nvwwgtQqpRj+xcih3C6QE9MTeTKzStcib9y3+c1wWuISohi49CNNCrbyOxS7SoxJpEtb28hYEYAhcoVYujGofh08nF8IfHxEB4On38Ob73l+P6FyEGyfaBrrYlPif/bkL5y8/9fR96MvH3jz70KuxemQuEKLOq3iCZeOXddMa01gQsD2fjqRhKiEmj6SlPafNSGPAVMWrr3lDFjiPImbU8nRA6S7QJ914VdTN4z+a7QTrQk3necQuGZz5OS+UtS0qMkjco2Mr6+9fjOzyXyl8jWa5Fn1JVjV1j/0nrObz9PmUZlGLJhCKUbmLg64V93gAK0bWteHULkENku0JMsSZy7fo6S+UtSxbPK34Z08fzFyZ0r2/3x7CIxOpGtH2zF/0d/3Au50+OnHjR4poHj55WDsVnz5s3GBdA1ayAyEh5/XMbNhbCBbJd4HSp14MhY2TsjIyzJFvyn+7Pj4x0kXU+i4XMNaftJW/J55nNsIampxpZwy5cbYZ6UBAULQrdu0KuXcVFUCJFp2S7QxcNpqyZwUSB/vPcH189ep1KHSnSa0omSdUxYndDf31ji9sgRY5x8zBgjxFu1Ajc3x9cjRA4mgZ7DhPqFsnniZiIORVCqfinzZq/ExcGHHxo3CpUsaZydP/GELK4lhB1JoOcQF/dcZNsH2ziz+QyFyheiz9w+1BlcB5XLwQGalmZsBffee8b4+HPPwaRJULiwY+sQwglJoGdzF3ZdYPtH2zmz+Qz5iuej05RONHqhEbnzmPBXu3mzcYfnsWPQogX8/js0buz4OoRwUhLo2ZDWmnPbzrHz052c/eMs+Uvkp+OXHfEd64tbfhPGpU+cMJa4XbcOKlaEpUtl7XIhTCCBno1Y06ycWH6CPV/s4bL/ZTxKedBpSid8n/PFNZ8JS8umphpn5D/+CPnzw+TJ8NJLkCfnz+kXIiuSQM8GkuOSOTL7CPu+3kfMmRiKVilKjxk9qDesnuO2gUvP4sUwbZoxc+XTT6F4cfNqEUJIoGdlEYcj8J/uz7H5x0iJT8GrqRcdv+xItV7VzLkp6F6bNkGxYsYZeq4sUI8QTk4CPYtJTUglcHEgAdMDuHTgErndc1NrYC18x/ri1dTL7PL+z2KB9euhUycJcyGyCAn0LCIyKJKAnwI4MucIyTeSKVa9GJ2ndqbesHrkLZrX7PL+T2sIDTXWYImKggEDzK5ICHGLBLqJYi/FErQkiMCFgVw+eBkXNxdq9KuB71hfvFt5Z43dkaKj4cAB2L/f+DhwwNgiDqBaNejSxdz6hBC3SaA7WEJUAseXHSdwUSDnd5wHDaUfK03HLztSb3g98hfPb26BiYnGXZ0bNxoBfvq08X2loGZN6N0bmjQxPmrVgtzyT0iIrEJ+Gx0gISqB02tPE7Q4iDN+Z7BarHhW8+TxDx6n9qDaFKtWzNwCtYaAAPj1V1i4EG7cMG7Xb9oURo40wtvX11hQSwiRZUmg24G2asL/DCd4fTAh60II2x8GGgp5F6Lpa02p81QdStYrmTWGVJYvN9ZcOXYM3N2hf3945hlo3VoudgqRzUig20hiTCJn/M4QvC6YkPUh3Iy8CQrKNi7L4x88TpVuVSjjWyZrhPidPvzQGCf/8UcYNEjWXBEiG5NAf0TJsclcOnCJc9vPce6Pc4TtD0OnadyLuFO5S2WqdKuCT2cf88fEHyQ11RheadUKxo41uxohRCZJoGeA1prokGjC9oZxcc9FwvaGceXYFdCgcinK+JahxcQWVO1elbKNy5IrdzYYqti717jD88IFePVVs6sRQtiABHo6kq4nEX4onEv7L90O8ISoBADyFMyDV1MvqvepTrnm5SjbpCzuhdxNrvgfSE6G11+HH36AsmVh5Upj5ooQIttz2kC3plmJvRhLdGg0MaExRIdGE306mogjEVw/e/32cZ7VPKnaoypezb0o16wcxWoUyxq33T+KyEjo0wf27IGXX4ZPPoECBcyuSghhIzk60C1JFmLOxNwV2jGhMUSHRHP93HWsqdbbx7q4uVCkUhHKNirLY6Mfo/RjpSnjW8bx+2/ay8mT0LUrREQYi2o9+aTZFQkhbCzbB3rS9STj7Drk7tCOCY0h9lIs6P8fm6dgHor4FKFUvVLU6FeDoj5FKeJThKI+RSlQtkD2PfN+GH9/447O3Llhxw5o1MjsioQQdpDtAj1sXxj7v91/O8AToxPvet6jlAdFfIpQsV1FI6wr/z+083rmzXrTBu0pJgYWLIC33wZPT/Dzg8qVza5KCGEnGQp0pVQX4BvABfhFaz3pnufVree7AQnASK31nzauFYCkG0lc2n+JIj5FqPlkzbvOsotUKoKbh5PvJG+xGMvazpoFq1ZBSopxp+dvvxkXQYUQOdZDA10p5QJ8D3QEwoCDSqnVWuvjdxzWFahy66MJ8OOtzzZXuXNlxoeOt0fT2d/PP8P77xvj5MWKwfPPw4gRUL++bAcnhBPIyBl6YyBEa30GQCm1COgN3BnovYE5WmsN7FNKFVZKldZah9u84o0bjW3PxP1OnzbWXJk+3bgA6ubk71aEcDIZCfSywMU7Hodx/9l3eseUBe4KdKXUGGAMgLe39z+t1VCwoLHqn7hf/fowaRKUK2d2JUIIE2Qk0NN7r64f4Ri01jOAGQC+vr73PZ8hzZoZu8oLIYS4S0bm6YUBd57yeQGXH+EYIYQQdpSRQD8IVFFKVVRKuQGDgNX3HLMaGK4MTYEbdhk/F0II8bceOuSitbYopV4ENmJMW5yptQ5SSo299fx0YB3GlMUQjGmLo+xXshBCiPRkaB661nodRmjf+b3pd3ytgXG2LU0IIcQ/kUPvdRdCCOcjgS6EEDmEBLoQQuQQEuhCCJFDKON6pgkdK3UVOO/gbosBUQ7uM7uR1+jB5PV5OHmNHiyzr095rXXx9J4wLdDNoJTy11r7ml1HViav0YPJ6/Nw8ho9mD1fHxlyEUKIHEICXQghcghnC/QZZheQDchr9GDy+jycvEYPZrfXx6nG0IUQIidztjN0IYTIsSTQhRAih3C6QFdKTVZKnVRKHVVKrVBKFTa7pqxGKTVAKRWklLIqpWT62S1KqS5KqVNKqRCl1Ftm15OVKKVmKqUilVKBZteSFSmlyimltiqlTtz63XrZHv04XaADfkBtrXVd4DTwtsn1ZEWBQF9gh9mFZBV3bJbeFagJPKWUkr0Q/28W0MXsIrIwC/C61roG0BQYZ49/P04X6FrrTVpry62H+zB2VxJ30Fqf0FqfMruOLOb2Zula6xTgr83SBaC13gFEm11HVqW1Dtda/3nr6zjgBMa+yzbldIF+j6eB9WYXIbKFv9sIXYh/RClVAWgA7Ld12xna4CK7UUptBkql89S7WutVt455F+Nt0HxH1pZVZOQ1EnfJ0EboQjyIUsoD+A14RWsda+v2c2Sga607POh5pdQIoAfQXjvpRPyHvUbiPrIRusgUpZQrRpjP11ovt0cfTjfkopTqAkwEemmtE8yuR2QbGdksXYh0KaUU8CtwQms9xV79OF2gA9OAAoCfUuqwUmr6w37A2Sil+iilwoBmwFql1EazazLbrQvpf22WfgJYorUOMreqrEMptRDYC1RTSoUppZ4xu6YspgUwDGh3K3cOK6W62boTufVfCCFyCGc8QxdCiBxJAl0IIXIICXQhhMghJNCFECKHkEAXQogcQgJdCCFyCAl0IYTIIf4HQ5EvU/csRNIAAAAASUVORK5CYII=\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": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1000/1000 [01:05<00:00, 15.19it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average fraction in bounds: 0.954\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# 1000 iterations with Normal distribution\n",
    "\n",
    "bounds = []\n",
    "for k in tqdm(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": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD4CAYAAAANbUbJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAajUlEQVR4nO3dfZBU9Z3v8fe3ewYGAWMCEyWC4WE1AXMTwAl58m7UvSagNxrduMG9N4q1KZYqJFLqXk1ZNw+VSpXZTVKUJQmFD1ErXkliHtabIrqr0buJqawMgoIKLhrUuQoC5mrxMMx09/f+cU4PzTBAA2f69Mz386rqmtOnf9Pn29Twmd/8zu/8jrk7IiIyvBXyLkBERAafwl5EJACFvYhIAAp7EZEAFPYiIgG05HXg8ePH++TJk/M6vIjIkLR27dqd7t5+rN+XW9hPnjyZzs7OvA4vIjIkmdkrx/N9GsYREQlAYS8iEoDCXkQkAIW9iEgACnsRkQCOGvZmdreZvWlmGw/zupnZbWa2xcyeNbPZ2ZcpIiInop6e/T3A3CO8Pg84M30sBH544mWJiEiWjjrP3t3/zcwmH6HJpcB9nqyV/EczO8XMJrj7GxnVmKl1m9bxq+/9CvbnXYmIRHX2BWfzNwv+pqHHzOKiqtOB12qed6X7Dgl7M1tI0vvnjDPOyODQx+7HS37MyY+ejKN1/EUkHxt6NwzJsLcB9g2YpO6+ElgJ0NHRkUvaFrcU2X72dn6w8Qd5HF5EJBdZzMbpAibVPJ8IvJ7B+w6KlndaKJ9SzrsMEZGGyiLsHwKuSmflfBx4u1nH6wGKe4swOu8qREQa66jDOGb2AHAeMN7MuoCvA60A7r4CWA1cBGwB9gLXDFaxJ8orTkt3CzZmoJEnEZHhq57ZOFce5XUHFmdW0SDq3dcLgI1S2ItILKGuoC11lwCwkQp7EYklZNiPaBuRcyUiIo0VMuzHjB6TcyUiIo0VKuy7u7sBaGtry7kSEZHGChX2pZ50zL5VY/YiEkuosO/en/TsW1tbc65ERKSxQoX9nn17ABjZNjLnSkREGitU2O/esxuA0aN0Ca2IxBIq7N/e+zYAY0ePzbkSEZHGChX21TH7k9pOyrkSEZHGChX2vT3Jcgk6QSsi0YQK+1IpmXrZ0pLFMv4iIkNHqLD3UnK/lEJLqI8tIhIr7CvlCqCevYjEEyrsy6XkDlXFlmLOlYiINFaosO/p7QGgbYTWxhGRWEKG/aiRo3KuRESksUKFfXXMXsM4IhJNqLD3cjobpxjqY4uIKOxFRCKIlXrJKA5W1Hr2IhJLqLD3StKzLxY1Zi8isYS6uqga9lZQz16k4dyhXE4elcqh2wPtO962zf5el18OX/pSQ//5Q4V9dTaOFkIL6Ej/GZs9GIbycWv3DQdmUChAsZg8qtsD7TvcdqEAb73V8NJDhX05/YFT2A+gXIb9+6GnJ/lauz0Y+6oB0KjAGg7MTjxkjvR6SwuMGJHNex1L2yzfa7CPa0N3VCBU2Fdn47QUc/rYlUrjwvRY92UdiK2tMHJk8hgx4tCvLS0H/6dqbR0a/9nzOm6hMKSDRvIXK+yrY/YnMhtnzx7YvBleeCF5bNoEO3fWF6rpEsuZqfbEDheq1e13vau+dvW8Xs++ESMUTCJNJlTYV4dxRrbUccPxnTsPBHrt49VXD7QpFmHaNDj1VBg9Gt7znmxD80iBXO0Ji4jUIVTYl0olWmgZeIljd/jNb+B734Nnn03CvmrUKPjgB+Hcc2H69GR7+nT4i79IgldEpMmFCvvDDuM8/TT8wz/Ab38LU6fCZZclYV4N9jPOSMZMRUSGqJhhX51nv3s3LF0Kd90F48bBbbfB3/99MlwiIjKM1NVdNbO5ZrbZzLaY2c0DvP4uM/vfZvaMmT1nZtdkX2oGkqxPwn7jRvjoR+Huu+HGG+Gll2DJEgW9iAxLR+3Zm1kRWA5cCHQBa8zsIXd/vqbZYuB5d/+cmbUDm83sfnfvGZSqj1d1bZynnoLzz4eTT4ZHH4ULLsi3LhGRQVbPMM4cYIu7vwxgZquAS4HasHdgrJkZMAZ4C8h4nuGJ84pTsQo8/DDs25f05idMyLssEZFBV88wzunAazXPu9J9tW4HpgOvAxuA69z9kKt0zGyhmXWaWeeOHTuOs+TjV6lUcHPYsQNOOUVBLyJh1BP2A10d4/2efxZYD7wPmAncbmYnH/JN7ivdvcPdO9rb24+52BPlFU8+TXd3Mi9eRCSIesK+C5hU83wiSQ++1jXALzyxBfgT8MFsSsxOpZz27IfLWikiInWqJ+zXAGea2RQzGwHMBx7q1+ZV4K8AzOxU4APAy1kWmoVKpZL07Ds7kwuiRESCOGrYu3sJuBZ4BHgB+Km7P2dmi8xsUdrsW8AnzWwD8Bhwk7vvHPgd8+MVxwsV2LABLr4473JERBqmrouq3H01sLrfvhU1268Dn8m2tOx5xQ+cbJg1K89SREQaKtQaAF7xZMwe4L3vzbcYEZEGChf2VMP+1FPzLUZEpIFihb074MmiZuPH512OiEjDxAr76jBOW5vWgheRUMKF/YCXiImIDHPxljg2h1atbCkiscTq2VfH7HV3KREJJlTYUyEZsz/ppLwrERFpqFBhb27JMI7CXkSCCRX2fZfPKuxFJJhwYe9UYNSovCsREWmoUGFfveG4evYiEk2osK94BSO9qEpEJJBQYU+F5ATtCM2zF5FYQoV9Mozj0BLqWjIRkVhhn8zGqSQ3GxcRCSRU2JsbjmvFSxEJJ1TYUwF0glZEAooV9p4uhFaI9bFFRGKlXiX9qouqRCSYUGFvni6ENnp03qWIiDRUqLDvm2ev2TgiEkyoCedWXS5BYS8iwcTq2VdXvdSYvYgEEyvsq1MvW1vzrkREpKFChX3FPVkITVMvRSSYUKnnThL2ZnmXIiLSUKHCHnfMtBCaiMQTLOwtOUmrsBeRYIKFfTpmr9k4IhJMXWFvZnPNbLOZbTGzmw/T5jwzW29mz5nZ/8m2zGy4AziMHJl3KSIiDXXU8QwzKwLLgQuBLmCNmT3k7s/XtDkF+AEw191fNbP3DlbBJ8LcknOzmo0jIsHUk3pzgC3u/rK79wCrgEv7tflb4Bfu/iqAu7+ZbZnZ8b4rq0RE4qgn7E8HXqt53pXuq3UW8G4ze8LM1prZVQO9kZktNLNOM+vcsWPH8VV8AsxNUy9FJKR6wn6gZOzfPW4BzgEuBj4L/E8zO+uQb3Jf6e4d7t7R3t5+zMWeMCdZCE1hLyLB1DMHsQuYVPN8IvD6AG12uvseYI+Z/RvwEeDFTKrMiHka8gp7EQmmnp79GuBMM5tiZiOA+cBD/dr8M/CfzazFzE4CPga8kG2pWVHPXkTiOWrP3t1LZnYt8AhQBO529+fMbFH6+gp3f8HMHgaeJVlu7E533ziYhR8Xt2RMSguhiUgwdV1K6u6rgdX99q3o9/yfgH/KrrRBUB2z1xW0IhJMqAnn5pZcWaV59iISTKzUc0sWQtOYvYgEEyzs0yWOtVyCiAQTKuyNtFOvnr2IBBMq7HGFvIjEFCrszS2ZjSMiEkyYsPdkfeNkzF5EJJg4YY+nC6GJiMQTJuyhOoyTdxUiIo0XL+xFRAIKE/bVMXv17EUkojBhX/HKgdsSiogEEybsy15W2ItIWGHC/sAwjtJeROKJE/bVqZfKehEJKEzYg66gFZG4woS9hnFEJLI4Ya9hHBEJLEzYlytlLZcgImGFCXtPF0BTz15EIooT9p4M41BQ2otIPHHCPh2zFxGJKE7Ye/UEraZeikg8YcK+4hVA66CJSEyhwl5TL0Ukqnhhn3chIiI5CBX2AKbZOCISUJiwL1VKac9eYS8i8YQL+6J69iISUJiw31faB0CLevYiElCYsN/buxdzo1VhLyIB1RX2ZjbXzDab2RYzu/kI7T5qZmUz+0J2JWajOowT5rebiEiNo2afmRWB5cA8YAZwpZnNOEy77wCPZF1kFg5cQauevYjEU09Hdw6wxd1fdvceYBVw6QDtlgA/B97MsL7M9E29VNaLSED1hP3pwGs1z7vSfX3M7HTgMmDFkd7IzBaaWaeZde7YseNYaz0hjlPwgkbsRSSkesJ+oHzsv5rYMuAmdy8f6Y3cfaW7d7h7R3t7e701ZqKn1ANAsaFHFRFpDi11tOkCJtU8nwi83q9NB7AqHQ8fD1xkZiV3/1UmVWbgQNirby8i8dQT9muAM81sCvB/gfnA39Y2cPcp1W0zuwf4dTMFPSS3JQQoaNBeRAI6ati7e8nMriWZZVME7nb358xsUfr6Ecfpm4VXdFtCEYmrnp497r4aWN1v34Ah7+4LTrys7JXKJUA9exGJKcw1RqVKEvbFlrp+v4mIDCthwr5Sri5xHOYji4j0CZN8B8bsNYwjIvGECfvqlQG6eYmIRBQm7N3Tnn1RYS8i8cQJew3jiEhg4cJeF9CKSERxwt51UZWIxBUn7MtJ2Bd0glZEAgoT9n2zcTSOIyIBhQn76s1LsP6rM4uIDH9hwr6vZ69BexEJKEzYV8fsNc9eRCKKE/auefYiEleYsC+XdfMSEYkrTNh3l7oBGFlozbkSEZHGCxP2Pb3JPWhbi7rluIjEEybsS6X05iUFhb2IxBMm7CuVZJ69xuxFJKIwYa8ljkUksjBhT7rqpXr2IhJRnLBPZl5S0D1oRSSgMMlXrlTn2edciIhIDsKE/b6efQC0WJiPLCLSJ0zy7etNwl5D9iISUZiw7+nZD4AVw3xkEZE+YZKvVOoFFPYiElOY5OtNr6DVOI6IRNSSdwGNUqokYa/zsyLDR29vL11dXXR3d+ddSuba2tqYOHEira3ZLN4YJuy9ovXsRYabrq4uxo4dy+TJk4fV/213Z9euXXR1dTFlypRM3rOufq6ZzTWzzWa2xcxuHuD1/2Zmz6aPP5jZRzKpLkPV9eyH0c+DSHjd3d2MGzduWAU9JJ3ScePGZfoXy1HD3syKwHJgHjADuNLMZvRr9ifg0+7+YeBbwMrMKsxId0/yj6YTtCLDy3AL+qqsP1c9yTcH2OLuL7t7D7AKuLS2gbv/wd3/nD79IzAx0yozsGf/HgBMl9CKSED1hP3pwGs1z7vSfYfzd8BvBnrBzBaaWaeZde7YsaP+KjPwTvc7SQ3q2YtIQPUk30BdYR+wodn5JGF/00Cvu/tKd+9w94729vb6q8zA7p60Z6+OvYhkaMyYMXmXUJd6ZuN0AZNqnk8EXu/fyMw+DNwJzHP3XdmUl52e3vSiKg3jiAxPS5fC+vXZvufMmbBsWbbvmZN6evZrgDPNbIqZjQDmAw/VNjCzM4BfAF9y9xezL/PE7OnZUzPPXmEvIo23YMECvvKVr/DJT36SqVOn8uCDDzb0+Eft2bt7ycyuBR4BisDd7v6cmS1KX18BfA0YB/wgPYNccveOwSv72Lz69quYJyGvMXuRYWoI9MDfeOMNfv/737Np0yYuueQSvvCFLzTs2HVdVOXuq4HV/fatqNn+MvDlbEvLzs69OxX2IpK7z3/+8xQKBWbMmMH27dsbeuxhn3ylSomn33haYS8iuRs5cmTfdvW+2I0ybJdLeHb7s9y7/l7u33A/2/ds52OtHwU0Zi8iMQ2rsO8p93DfM/exfM1y1m9bT2uhlYvPupgFH1nAtCfG8vN7f4eNyGZRIRERgL179zJx4oHrSK+//nquv/76HCsa2LAI+329+7hr3V1858nv0PVOFzNPm8ltc2/jyv90JeNPGg/AC4/+FoBCazHPUkVkmKlUKnW1u+eeew56vnv37kGo5vCGdNjvL+1nRecKbn3yVrbt3sanJn2KOz93J5+Z9plD1pWo7EvXxjlpVB6liojkakiGfblS5v4N9/O1x7/GK2+/wvmTz+eBv36AT7//04ddPMj3p7clHNXWyFJFJJhvf/vb/OxnPzto3xVXXMEtt9ySU0WJIRf2f3jtDyz69SI2vLmB2RNmc8fn7uDCaRce9fu8nPypZUUN44jI4LnllltyD/aBDLmwbym0sK+0j1V/vYorzr6CQp23nvJSGvYtmo0jIvEMubCfc/ocNi3eRLFwbD30SjmZ01rQPHsRCWhIJt+xBj3UDOO0DMmPLCJyQsIkX6UnWQit0DbyKC1FRIafMGHvPekSxyOG3MiViDSprVu38qEPfeigfd/4xjf47ne/m1NFhxcm+SqVdMy+RbNxRIajpQ8vZf22bNezn3naTJbNbf7VNOsRp2ffNxtHYS8i+TjvvPO46aabmDNnDmeddRa/+93vGnbsMD37vhO0GsYRGZaGSg+8VCrx1FNPsXr1ar75zW/y6KOPNuS4YXr2lbRnX2hV2ItINg53xf7h9gNcfvnlAJxzzjls3bp1MMoaUJywT3v2BU29FJGMjBs3jj//+c8H7XvrrbcYP378Yb+nuqZ9sVikVCoNan21wiRf30VVGrMXkYyMGTOGCRMm8NhjjwFJ0D/88MOce+65OVd2qDBjGpVy8rXQGub3m4g0wH333cfixYu54YYbAPj617/OtGnTcq7qUIHC3jEqulOViGRqxowZPP7443W1feKJJ/q2x48frzH7wVApOwUqcIQTJyIiw1Wonn2B+u4oIyJyIhYvXsyTTz550L7rrruOa665JqeKIoa9evYiMsiWL1+edwmHCDOMUy4p7EUkrjBhXyk7RcoKexEJKUzYl3vTsBcRCShO2JfUsxeR7G3bto358+czbdo0ZsyYwUUXXcSLL77IqFGjmDVrFtOnT2fOnDnce++9fd9zzz330N7ezsyZM5k5cyZXXXXVoNcZ5gRtqbdCC427NFlEhj9357LLLuPqq69m1apVAKxfv57t27czbdo01q1bB8DLL7/M5ZdfTqVS6ZuR88UvfpHbb7+9YbWGCfve/WnYF8L8MSMSysNLH2bb+m2ZvudpM09j7rK5h3398ccfp7W1lUWLFvXtmzlz5iEXS02dOpXvf//73HDDDblNvwyTfKWeCq30QkuY328iMsg2btzIOeecU1fb2bNns2nTpr7nP/nJT/qGcX70ox8NVol9wiTf/r1lTmY/jBqVdykiMgiO1ANvBu5+0PNGD+PU1bM3s7lmttnMtpjZzQO8bmZ2W/r6s2Y2O/tST0z33gptdMNI3XBcRLJx9tlns3bt2rrarlu3junTpw9yRYd31LA3syKwHJgHzACuNLMZ/ZrNA85MHwuBH2Zc5wlxd/budkaxD4pa4lhEsnHBBRewf/9+7rjjjr59a9as4ZVXXjmo3datW7nxxhtZsmRJo0vsU88wzhxgi7u/DGBmq4BLgedr2lwK3OfJ3yl/NLNTzGyCu7+RdcFbvv0T/uVb/37ENhU3yl6g5EXKXqC73IpT4DS26QStiGTGzPjlL3/J0qVLufXWW2lra2Py5MksW7aMl156iVmzZtHd3c3YsWNZsmRJ06+NczrwWs3zLuBjdbQ5HTgo7M1sIUnPnzPOOONYawVgZPvJtB/+JjDpcZxioZdisYeWglMsOKNHljj7s3M1Zi8imXrf+97HT3/600P279u377Dfs2DBAhYsWDCIVR2qnrAf6CokP442uPtKYCVAR0fHIa/XY9LCeUxaOO94vlVEJKx6xjS6gEk1zycCrx9HGxERyUk9Yb8GONPMppjZCGA+8FC/Ng8BV6Wzcj4OvD0Y4/UiIv31n9I4XGT9uY46jOPuJTO7FngEKAJ3u/tzZrYofX0FsBq4CNgC7AXyOwshImG0tbWxa9cuxo0bhw2jda/cnV27dtHW1pbZe1pevxU7Ojq8s7Mzl2OLyPDQ29tLV1cX3d3deZeSuba2NiZOnEhra+tB+81srbt3HOv7hbmCVkSGn9bWVqZMmZJ3GUOCJp2LiASgsBcRCUBhLyISQG4naM1sB/DKURtmbzywM4fjHg/VOniGUr2qdXAM1Vrf7+7tx/oGuYV9Xsys83jOZOdBtQ6eoVSvah0c0WrVMI6ISAAKexGRACKG/cq8CzgGqnXwDKV6VevgCFVruDF7EZGIIvbsRUTCUdiLiAQQJuyPdtP0PJjZ3Wb2ppltrNn3HjP7VzP7j/Tru2te+2pa/2Yz+2yDa51kZo+b2Qtm9pyZXdes9ZpZm5k9ZWbPpLV+s1lrrTl+0czWmdmvm7lWM9tqZhvMbL2ZdTZ5raeY2YNmtin9uf1EE9f6gfTftPp4x8yWZlqvuw/7B8nSzC8BU4ERwDPAjCao6y+B2cDGmn3/CNycbt8MfCfdnpHWPRKYkn6eYgNrnQDMTrfHAi+mNTVdvSR3ThuTbrcC/w58vBlrran5euB/Ab9u8p+DrcD4fvuatdZ7gS+n2yOAU5q11n51F4FtwPuzrLfhHySnf7xPAI/UPP8q8NW860prmczBYb8ZmJBuTwA2D1Qzyf0FPpFj3f8MXNjs9QInAU+T3De5KWslubPbY8AFNWHfrLUOFPZNVytwMvAn0kkozVzrALV/Bngy63qjDOMc7obozehUT+/ylX59b7q/aT6DmU0GZpH0mJuy3nRYZD3wJvCv7t60tQLLgP8BVGr2NWutDvyLma01s4XpvmasdSqwA/hROjx2p5mNbtJa+5sPPJBuZ1ZvlLCv64boTa4pPoOZjQF+Dix193eO1HSAfQ2r193L7j6TpNc8x8w+dITmudVqZv8VeNPd19b7LQPsa+TPwafcfTYwD1hsZn95hLZ51tpCMkT6Q3efBewhGQY5nLz/XZMiklu/XgL87GhNB9h3xHqjhP1QuiH6djObAJB+fTPdn/tnMLNWkqC/391/ke5u2noB3P3/AU8Ac2nOWj8FXGJmW4FVwAVm9uMmrRV3fz39+ibwS2BOk9baBXSlf9EBPEgS/s1Ya615wNPuvj19nlm9UcK+npumN4uHgKvT7atJxsar++eb2UgzmwKcCTzVqKLMzIC7gBfc/fvNXK+ZtZvZKen2KOC/AJuasVZ3/6q7T3T3ySQ/l7919//ejLWa2WgzG1vdJhlb3tiMtbr7NuA1M/tAuuuvgOebsdZ+ruTAEE61rmzqzeMERE4nPS4imUHyEnBL3vWkNT0AvAH0kvym/jtgHMnJuv9Iv76npv0taf2bgXkNrvVckj8TnwXWp4+LmrFe4MPAurTWjcDX0v1NV2u/us/jwAnapquVZBz8mfTxXPX/UTPWmh57JtCZ/hz8Cnh3s9aaHv8kYBfwrpp9mdWr5RJERAKIMowjIhKawl5EJACFvYhIAAp7EZEAFPYiIgEo7EVEAlDYi4gE8P8BolB2WwZg7H0AAAAASUVORK5CYII=\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": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1000/1000 [01:04<00:00, 15.42it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average fraction in bounds: 0.209\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# 1000 iterations with Cauchy distribution\n",
    "\n",
    "bounds = []\n",
    "for k in tqdm(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": {
      "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 numpy as np\n",
    "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: 0.824\n"
     ]
    }
   ],
   "source": [
    "# Estimate the mean waiting time and give a standard error for the estimate.\n",
    "theta = r.mean()\n",
    "se = np.std(r, ddof=1) / np.sqrt(len(r))\n",
    "\n",
    "## Alternatively:\n",
    "# from scipy.stats import sem\n",
    "# se = sem(r)\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: (69.541, 72.253)\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_1 = n_2 = 100$, so we can use $n$ to refer to both.\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{\\frac{\\hat{p_1}(1 - \\hat{p_1})}{n} + \\frac{\\hat{p_2}(1 - \\hat{p_2})}{n}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated mean: 0.050\n",
      "Estimated SE: 0.047\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "p_hat_1 = 0.9\n",
    "p_hat_2 = 0.85\n",
    "n = 100\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",
    "\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.010, 0.110)\n",
      "95% confidence interval: (-0.041, 0.141)\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\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": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 10000/10000 [00:05<00:00, 1863.74it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated SE: 63.420\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\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(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 53.099, 301.701\n",
      "95% confidence interval (percentile): \t 38.336, 266.201\n",
      "95% confidence interval (pivotal): \t 88.599, 316.464\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
}