Chapter 09 - The Bootstrap.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 9. The Bootstrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $X_1, \\dots, X_n \\sim F$ be random variables distributed according to $F$, and\n",
    "\n",
    "$$ T_n = g(X_1, \\dots, X_n)$$\n",
    "\n",
    "be a **statistic**, that is, any function of the data.  Suppose we want to know $\\mathbb{V}_F(T_n)$, the variance of $T_n$.\n",
    "\n",
    "For example, if $T_n = n^{-1}\\sum_{i=1}^nX_i$ then $\\mathbb{V}_F(T_n) = \\sigma^2/n$ where $\\sigma^2 = \\int (x - \\mu)^dF(x)$ and $\\mu = \\int x dF(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9.1 Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we draw iid samples $Y_1, \\dots, Y_B \\sim G$.  By the law of large numbers,\n",
    "\n",
    "$$ \\overline{Y}_n = \\frac{1}{B} \\sum_{j=1}^B Y_j \\; \\xrightarrow{\\text{P}} \\int y dG(y) = \\mathbb{E}(Y)$$\n",
    "\n",
    "as $B \\rightarrow \\infty$.  So if we draw a large sample from $G$, we can use the sample mean to approximate the mean of the distribution.\n",
    "\n",
    "More generally, if $h$ is any function with finite mean then, as $B \\rightarrow \\infty$,\n",
    "\n",
    "$$\\frac{1}{B} \\sum_{j=1}^B h(Y_j) \\; \\xrightarrow{\\text{P}} \\int h(y) dG(y) = \\mathbb{E}(h(Y))$$\n",
    "\n",
    "In particular, for functions $a(Y_j) = Y_j^2$ and $b(Y_j) = Y_j$,\n",
    "\n",
    "$$\\frac{1}{B} \\sum_{j=1}^B (Y_j - \\overline{Y}) \n",
    "= \\frac{1}{B} \\sum_{j=1}^B Y_j^2 - \\left(\\frac{1}{B} \\sum_{j=1}^B Y_j \\right)^2\n",
    "\\xrightarrow{\\text{P}} \\int y^2 dG(y) - \\left( \\int y dG(y) \\right)^2 = \\mathbb{V}(Y)\n",
    "$$\n",
    "\n",
    "So we can use the sample variance of the simulated values to approximate $\\mathbb{V}(Y)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9.2 Bootstrap Variance Estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To simulate from the distribution of a statistic $T_n$ when the data is assumed to have distribution $\\hat{F}_n$, we simulate $X_1^*, \\dots, X_n^*$ from $\\hat{F}_n$ and then compute the statistic over these values, $T_n^* = g(X_1^*, \\dots, X_n^*)$.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\text{Real world} \\quad      & F         & \\Longrightarrow\\quad & X_1, \\dots, X_n        & \\Longrightarrow & \\quad T_n = g(X_1, \\dots, X_n) \\\\\n",
    "\\text{Bootstrap world} \\quad & \\hat{F}_n & \\Longrightarrow\\quad  & X_1^*, \\dots, X_n^* & \\Longrightarrow & \\quad T_n^* = g(X_1^*, \\dots, X_n^*)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "**Drawing an observation from $\\hat{F}_n$ is equivalent to drawing one point at random from the original data set.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bootstrap Variance Estimation\n",
    "\n",
    "1. Draw $X_1^*, \\dots, X_n^* \\sim \\hat{F}_n$.\n",
    "2. Compute $T_n^* = g(X_1^*, \\dots, X_n^*)$.\n",
    "3. Repeat steps 1 and 2, $B$ times, to get $T_{n, 1}^*, \\dots, T_{n, B}^*$.\n",
    "4. Let\n",
    "\n",
    "$$ v_{\\text{boot}} = \\frac{1}{B} \\sum_{b=1}^B \\left( T_{n, b}^* - \\frac{1}{B} \\sum_{r=1}^B T_{n, r}^* \\right)^2 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9.3 Bootstrap Confidence Intervals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Normal Interval**.\n",
    "\n",
    "$$ T_n \\pm z_{\\alpha/2} \\hat{\\text{se}}_\\text{boot} $$\n",
    "\n",
    "where $\\hat{\\text{se}}_\\text{boot}$ is the bootstrap estimate of the standard error.  This is not accurate unless the distribution of $T_n$ is close to Normal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Pivotal Intervals**.\n",
    "\n",
    "Let $\\theta = T(F)$, $\\hat{\\theta}_n = T(\\hat{F}_n)$ and define the **pivot** $R_n = \\hat{\\theta}_n  - \\theta$.  Let $\\hat{\\theta}_{n, 1}^*, \\dots, \\hat{\\theta}_{n, B}^*$ define bootstrap replications of $\\hat{\\theta}_n$.  Let $H(r)$ denote the CDF of the pivot:\n",
    "\n",
    "$$ H(r) = \\mathbb{P}_F(R_n \\leq r)$$\n",
    "\n",
    "Define the interval $C_n^* = (a, b)$ where\n",
    "\n",
    "$$\n",
    "a = \\hat{\\theta}_n - H^{-1}\\left( 1 - \\frac{\\alpha}{2} \\right) \n",
    "\\quad\\text{and}\\quad\n",
    "b = \\hat{\\theta}_n - H^{-1}\\left( \\frac{\\alpha}{2} \\right) \n",
    "$$\n",
    "\n",
    "Then,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{P}(a \\leq \\theta \\leq b) &= \\mathbb{P}(a - \\hat{\\theta}_n \\leq \\theta - \\hat{\\theta}_n \\leq b - \\hat{\\theta}_n) \\\\\n",
    "&= \\mathbb{P}(\\hat{\\theta}_n - b \\leq \\hat{\\theta}_n - \\theta \\leq \\hat{\\theta}_n - a) \\\\\n",
    "&= \\mathbb{P}(\\hat{\\theta}_n - b \\leq R_n \\leq \\hat{\\theta}_n - a) \\\\\n",
    "&= H(\\hat{\\theta}_n - a) - H(\\hat{\\theta}_n - b) \\\\\n",
    "&= H\\left(H^{-1}\\left( 1 - \\frac{\\alpha}{2} \\right) \\right) - H\\left(H^{-1}\\left( \\frac{\\alpha}{2} \\right)\\right) \\\\\n",
    "&= 1 - \\frac{\\alpha}{2}  - \\frac{\\alpha}{2} = 1 - \\alpha\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Hence $C_n^*$ is an exact $1 - \\alpha$ confidence interval for $\\theta$.  \n",
    "\n",
    "Unfortunately, $a$ and $b$ depend on the unknown distribution $H$, but we can form a bootstrap estimate for it:\n",
    "\n",
    "$$\\hat{H}(r) = \\frac{1}{B} \\sum_{b=1}^B I(R_{n, b}^* \\leq r)$$\n",
    "\n",
    "where $R_{n, b}^* = \\hat{\\theta}_{n, b}^* - \\hat{\\theta}_n$.\n",
    "\n",
    "Let $r_\\beta^*$ denote the $\\beta$ sample quantile of $(R_{n, 1}^*, \\dots, R_{n, B}^*)$, and let $\\theta_\\beta^*$ denote the $\\beta$ sample quantile of $(\\theta_{n, 1}^*, \\dots, \\theta_{n, B}^*)$.  Note that $r_\\beta^* = \\theta_\\beta^* - \\hat{\\theta}_n$.  It follows an approximate $1 - \\alpha$ confidence interval $C_n = (\\hat{a}, \\hat{b})$ where\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{a} \n",
    "&= \\hat{\\theta}_n - \\hat{H}^{-1}\\left(1 - \\frac{\\alpha}{2}\\right) \n",
    "&= \\hat{\\theta}_n - r_{1 - \\alpha/2}^* \n",
    "&= 2\\hat{\\theta}_n - \\theta_{1 - \\alpha/2}^* \\\\\n",
    "\\hat{b} \n",
    "&= \\hat{\\theta}_n - \\hat{H}^{-1}\\left(\\frac{\\alpha}{2}\\right) \n",
    "&= \\hat{\\theta}_n - r_{\\alpha/2}^* \n",
    "&= 2\\hat{\\theta}_n - \\theta_{\\alpha/2}^* \n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **$1 - \\alpha$ bootstrap pivotal confidence** is\n",
    "\n",
    "$$ C_n = \\left(2 \\hat{\\theta}_n - \\hat{\\theta}_{1 - \\alpha/2}^*, \\; 2 \\hat{\\theta}_n - \\hat{\\theta}_{\\alpha/2}^* \\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 9.3**.  Under weak conditions on $T(F)$, \n",
    "\n",
    "$$ \\lim _{n \\rightarrow \\infty} \\mathbb{P}_F\\left(T(F) \\in C_n\\right) \\rightarrow 1 - \\alpha$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Percentile Intervals**.\n",
    "\n",
    "The **bootstrap percentile interval** is defined by\n",
    "\n",
    "$$ C_n = \\left( \\theta_{\\alpha/2}^*, \\; \\theta_{1 - \\alpha/2}^*\\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9.5 Technical Appendix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The Jacknife\n",
    "\n",
    "This method is less computationally expensive than bootstraping, but it is less general -- it does *not* produce consistent estimates of the standard errors of the sample quantiles.\n",
    "\n",
    "Let $T_n = T(X_1, \\dots, X_n)$ be a statistic and let $T_{(-i)}$ denote the statistic with the $i$-th observation removed.  Let $\\overline{T}_n = n^{-1} \\sum_{i=1}^n T_{(-1)}$.  The jacknife estimate of $\\mathbb{V}(T_n)$ is\n",
    "\n",
    "$$ v_\\text{jack} = \\frac{n - 1}{n} \\sum_{i=1}^n \\left(T_{(-i)} - \\overline{T}_n \\right)^2 $$\n",
    "\n",
    "and the jacknife estimate of the standard error is $\\hat{\\text{se}}_\\text{jack} = \\sqrt{v_\\text{jack}}$.\n",
    "\n",
    "Under suitable conditions on $T$ it can be shown that $v_\\text{jack} / \\mathbb{V}(T_n) \\xrightarrow{\\text{P}} 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Justification for the Percentile Interval\n",
    "\n",
    "Suppose there exists a monotone transformation $U = m(T)$ such that $U \\sim N(\\phi, c^2)$ where $\\phi = m(\\theta)$.  \n",
    "\n",
    "Let $U_t^* = m(\\theta_{n, b}^*)$.  Let $u_\\beta^*$ be the sample quantile of the $U_b^*$'s.  Since a monotone transformation preserves quantiles, we have that $u_{\\alpha/2}^* = m(\\theta_{\\alpha/2}^*)$.  \n",
    "\n",
    "Also, since $U \\sim N(\\phi, c^2)$ the $\\alpha/2$ quantile of $U$ is $\\phi - z_{\\alpha/2}c$.  Hence $u_{\\alpha/2}^* = \\phi - z_{\\alpha/2}c$.  Similarly, $u_{1 - \\alpha/2}^* = \\phi + z_{\\alpha/2}c$.  \n",
    "\n",
    "Therefore,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{P}(\\theta_{\\alpha/2}^* \\leq \\theta \\leq \\theta_{1 - \\alpha/2}^*) &=\n",
    "\\mathbb{P}(m(\\theta_{\\alpha/2}^*) \\leq m(\\theta) \\leq m(\\theta_{1 - \\alpha/2}^*)) \\\\\n",
    "&= \\mathbb{P}(u_{\\alpha/2}^* \\leq \\phi \\leq u_{1 - \\alpha/2}^*) \\\\\n",
    "&= \\mathbb{P}(U - cz_{\\alpha/2} \\leq \\phi \\leq U + cz_{1 - \\alpha/2}) \\\\\n",
    "&= \\mathbb{P}\\left(-z_{\\alpha/2} \\leq \\frac{Y - \\phi}{c} \\leq z_{1 - \\alpha/2} \\right) \\\\\n",
    "&= 1 - \\alpha\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9.6 Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 9.6.1**.  Consider the data in example 9.6.  \n",
    "\n",
    "- Find the plug-in estimate of the correlation coefficient.  \n",
    "- Estimate the standard error using the bootstrap.  \n",
    "- Find a 95% confidence interval using all three methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data from example 9.6:\n",
    "\n",
    "LSAT = [576, 635, 558, 578, 666, 580, 555, 661, 651, 605, 653, 575, 545, 572, 594]\n",
    "GPA = [3.39, 3.30, 2.81, 3.03, 3.44, 3.07, 3.00, 3.43, 3.36, 3.13, 3.12, 2.74, 2.76, 2.88, 3.96]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x16cf9d8cfd0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAWxElEQVR4nO3df5BVZ33H8fdnYbMQSQKFTUrZpKTGxNiUkLiiLaM1mElIVGKHdvxVzVBTjI1OmraC1qltxumMQeuv0YSh8QdO00YMMb8M0XRSdNABXCIsQWKlJoYNadkgaDaBzcJ++8c9K5ebuz/Zc/fc+3xeMzucc+6zd78Py70fznme+xxFBGZmlq6miS7AzMwmloPAzCxxDgIzs8Q5CMzMEucgMDNL3OSJLmC0Zs2aFXPnzp3oMszM6sq2bduejYjWao/VXRDMnTuXjo6OiS7DzKyuSPrFYI/lfmlI0iRJP5b0QJXHJOkLkvZI6pR0ad71mJnZiWoxRnAjsHuQx64CXpF9LQduq0E9ZmZWJtcgkNQGvBm4fZAm1wBfj5LNwHRJs/OsyczMTpT3GcHngBVA/yCPzwH2lu13ZcfMzKxGcgsCSW8B9kfEtqGaVTn2ksWPJC2X1CGpo7u7e9xqNDOzfM8IFgJLJD0J3AkskvRvFW26gLPL9tuAfZVPFBFrIqI9ItpbW6vOfjIzszHKLQgi4qMR0RYRc4F3AI9ExJ9XNLsPeG82e+h1wK8i4pm8ajIzs5eq+ecIJF0PEBGrgQeBq4E9wAvAslrXY43pQE8vXQcP0zZjKjOntUx0OWaFVpMgiIiNwMZse3XZ8QBuqEUNlo57tz/NyvWdNDc10dffz6ql81gy33MQzAbjtYasoRzo6WXl+k6O9PXzXO9RjvT1s2J9Jwd6eie6NLPCchBYQ+k6eJjmphP/WTc3NdF18PAEVWRWfA4CayhtM6bS13/ix1b6+vtpmzF1gioyKz4HgTWUmdNaWLV0HlOamzitZTJTmptYtXSeB4zNhlB3q4+aDWfJ/DksPG+WZw2ZjZCDwBrSzGktDgCzEfKlITOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxOUWBJKmSNoqaYekXZJurtLmDEn3l7VZllc9ZmZWXZ63quwFFkVEj6RmYJOkDRGxuazNDcBPIuKtklqBn0q6IyJezLEuMzMrk1sQREQAPdluc/YVlc2A0yQJmAb8EjiaV01mZvZSuY4RSJokaTuwH3g4IrZUNPkicCGwD9gJ3BgR/VWeZ7mkDkkd3d3deZZsZpacXIMgIo5FxHygDVgg6aKKJlcC24HfAeYDX5R0epXnWRMR7RHR3trammfJZmbJqcmsoYg4BGwEFlc8tAy4O0r2AE8Ar6xFTWZmVpLnrKFWSdOz7anA5cDjFc2eAt6UtTkLuAD4eV41mZnZS+U5a2g2sFbSJEqBsy4iHpB0PUBErAY+AXxN0k5AwMqIeDbHmszMrEKes4Y6gUuqHF9dtr0PuCKvGszMbHj+ZLGZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklLrcgkDRF0lZJOyTtknTzIO3eKGl71uZ7edVjZmbVTc7xuXuBRRHRI6kZ2CRpQ0RsHmggaTpwK7A4Ip6SdGaO9ZiZWRW5BUFEBNCT7TZnX1HR7F3A3RHxVPY9+/Oqx8zMqst1jEDSJEnbgf3AwxGxpaLJ+cAMSRslbZP03kGeZ7mkDkkd3d3deZZsZpacXIMgIo5FxHygDVgg6aKKJpOBVwNvBq4E/kHS+VWeZ01EtEdEe2tra54lm5klpyazhiLiELARWFzxUBfwUEQ8HxHPAt8HLq5FTWZmVpLnrKHWbDAYSVOBy4HHK5rdC7xe0mRJpwKvBXbnVZOZmb1UnrOGZgNrJU2iFDjrIuIBSdcDRMTqiNgt6SGgE+gHbo+Ix3KsyczMKqg0uad+tLe3R0dHx0SXYWZWVyRti4j2ao/5k8VmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZ1YEDPb3s2HuIAz294/7ced6z2MzMqjjQ00vXwcO0zZjKzGktw7a/d/vTrFzfSXNTE339/axaOo8l8+eMWz0OAjOzGhrtm/qBnl5Wru/kSF8/R+gHYMX6ThaeN2tEITISvjRkZlYj5W/qz/Ue5UhfPyvWdw55uafr4GGam058q25uaqLr4OFxq8tBYGZWI2N5U2+bMZW+/v4TjvX199M2Y+q41eUgMDOrkbG8qc+c1sKqpfOY0tzEaS2TmdLcxKql88btshDkOEYgaQrwfaAl+zl3RcQ/DtL2NcBm4O0RcVdeNZmZTaSBN/UVFWMEw72pL5k/h4XnzRrVAPNo5DlY3AssiogeSc3AJkkbImJzeSNJk4BbgO/kWIuZWSGM9U195rSWcQ+AAbkFQUQE0JPtNmdfUaXph4D1wGvyqsXGbrTT3MxseHm+qY/FmIIgu+zz1oj45jDtJgHbgPOAL0XElorH5wB/AixiiCCQtBxYDnDOOeeMpWQbg7znLptZMYx4sFjSJElXSfo68Avg7cN9T0Qci4j5QBuwQNJFFU0+B6yMiGPDPM+aiGiPiPbW1taRlmwnYSzT3MysPg17RiDpDcC7gDcDW4GFwLkR8cJIf0hEHJK0EVgMPFb2UDtwpySAWcDVko5GxD0j7oHlYmCa28AHWOD4NLcindKa2ckb8oxAUhfwSeAHwKsiYilweCQhIKlV0vRseypwOfB4eZuIODci5kbEXOAu4K8cAsVQi7nLZkWU55o+RTXcGcF64G2ULgMdk3Qv1Qd8q5kNrM3GCZqAdRHxgKTrASJi9RhrthoY6zQ3s3qW6riYSpN7hmhQum5zGfBO4GrgdOB9wIMR0TPU9+ahvb09Ojo6av1jk+VZQ5aKAz29LLzlEY70HT8TntLcxA9WLmqIf/uStkVEe7XHhh0jyKaBPgI8kn0eYDGlULiV0nV9a2BFm+ZmlpeUx8VGMlh8CfByYFdE7AbuB+7PrvubmTWElMfFhhss/jjwDWAp8KCkvxx4LCLGb+k7M7MJVos1fYpquDOCtwPzI+IFSTOBh4B/zb8sM7Pay3tNn6IaLgiODEwVjYgDkrxaqZk1tBTHxYYLgpdLui/bVtm+KI0jL8m1OjMzy91wQXBNxf6nOf45Ao1/OWZmVmvDBcF0oC0ivgQgaSvQSikMVuZcm5k1IH82pXiGC4IVwDvK9k+htD7Qy4CvAkOuPmpmVi7VT+4W3XCDv6dExN6y/U0RcSAinqIUBmZmI+IVbYtruCCYUb4TER8s2/V60GY2YmO5cbvVxnBBsKX8Q2QDJL2f0pLUZmYjkvInd4tuuDGCm4B7JL0LeDQ79mpKN6R/W56FmVlj8Yq2xTVkEETEfuCPJC0Cfj87/O2IeCT3ysys4aT6yd2iG9E9i7M3fr/5m9lJS/GTu0XnJSPMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PE5RYEkqZI2ipph6Rdkm6u0ubdkjqzrx9KujiveszMrLoRfbJ4jHqBRRHRI6kZ2CRpQ0RsLmvzBPDHEXFQ0lXAGuC1OdZkZmYVcguCiAigJ9ttzr6ios0Py3Y3A2151WNmZtXlOkYgaZKk7cB+4OGI2DJE8/cBGwZ5nuWSOiR1dHd351GqmVmycg2CiDgWEfMp/U9/gaSLqrWTdBmlIKh6H+SIWBMR7RHR3trq++GYmY2nmswaiohDwEZgceVjkuYBtwPXRMSBWtRjZmbH5TlrqFXS9Gx7KnA58HhFm3OAu4H3RMR/51WLmZkNLs9ZQ7OBtZImUQqcdRHxgKTrASJiNfBxYCZwqySAoxHRnmNNZmZWIc9ZQ53AJVWOry7bvg64Lq8azMxseP5ksZlZ4hwEZmaJcxCYnaQDPb3s2HuIAz29E12K2ZjkOVhs1vDu3f40K9d30tzURF9/P6uWzmPJ/DkTXZbZqPiMwGyMDvT0snJ9J0f6+nmu9yhH+vpZsb7TZwYNqNHP+nxGYDZGXQcP09zUxBH6f3OsuamJroOHmTmtZQIrs/GUwlmfzwjMxqhtxlT6+vtPONbX30/bjKkTVJGNt1TO+hwEZmM0c1oLq5bOY0pzE6e1TGZKcxOrls7z2UADGTjrKzdw1tdIfGnI7CQsmT+HhefNouvgYdpmTHUINJhUzvp8RmB2kmZOa+His6c7BBpQKmd9PiMwMxtCCmd9DgIzs2HMnNbSkAEwwJeGbFw0+jxrs0bmMwI7aSnMszZrZD4jsJOSyjxrs0bmILCTkso8a7NG5iCwk5LKPGuzRuYgsJOSyjxrs0bmwWIbtQM9vSfMqU5hnrVZI3MQ2KgMNkOo0edZmzUyXxqyEfMMIbPG5CCwEfMMIbPG5CCwEfMMIbPGlFsQSJoiaaukHZJ2Sbq5ShtJ+oKkPZI6JV2aVz210shLLXiGkFljynOwuBdYFBE9kpqBTZI2RMTmsjZXAa/Ivl4L3Jb9WZdSWGrBM4TMGk9uZwRR0pPtNmdfUdHsGuDrWdvNwHRJs/OqKU8pDaR6/X2zxpLrGIGkSZK2A/uBhyNiS0WTOcDesv2u7Fjl8yyX1CGpo7u7O7+CT4IHUs2sXuUaBBFxLCLmA23AAkkXVTRRtW+r8jxrIqI9ItpbW1vzKPWkeSDVzOpVTWYNRcQhYCOwuOKhLuDssv02YF8tahpvHkg1s3qV22CxpFagLyIOSZoKXA7cUtHsPuCDku6kNEj8q4h4Jq+a8uaBVDOrR3nOGpoNrJU0idKZx7qIeEDS9QARsRp4ELga2AO8ACzLsZ6a8FILZlZvcguCiOgELqlyfHXZdgA35FWDmZkNz58sNjNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDoJENfItNc1sdPJcdM4KKoVbaprZyPmMIDEp3VLTzEbGQZAY31LTzCo5CBLjW2qaWSUHQWJ8S00zq+TB4gT5lppmVs5BkCjfUtPMBvjSkJlZ4hwEZmaJcxCYmSXOQVBDXtbBzIrIg8U14mUdzKyocjsjkHS2pP+StFvSLkk3VmlzhqT7Je3I2izLq56J5GUdzKzI8rw0dBT424i4EHgdcIOkV1W0uQH4SURcDLwR+BdJp+RY04Twsg5mVmS5BUFEPBMRj2bbzwG7gcprIQGcJknANOCXlAKkoXhZBzMrspoMFkuaC1wCbKl46IvAhcA+YCdwY0T0V7RB0nJJHZI6uru7c652/HlZBzMrMkVEvj9AmgZ8D/jniLi74rE/BRYCfwO8HHgYuDgifj3Y87W3t0dHR0eOFefnQE+vl3UwswkhaVtEtFd7LNczAknNwHrgjsoQyCwD7o6SPcATwCvzrGkizZzWwsVnT3cImFmh5DlrSMCXgd0R8ZlBmj0FvClrfxZwAfDzPOrxHH4zs+ry/BzBQuA9wE5J27Njfw+cAxARq4FPAF+TtBMQsDIinh3vQjyH38xscLkFQURsovTmPlSbfcAVedUAJ87hP0JpHHrF+k4WnjfLl2jMzEhgiQnP4R8bX0ozS0fDLzHhOfyj50tpZmlp+DMCz+EfHS+HYZaehj8jAN+acTQGLqUNjKfA8Utp/nsza0xJBAH41owj5UtpZulp+EtDNjq+lGaWnmTOCGzkfCnNLC0OAqvKl9LM0uFLQ2ZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpa43O9ZPN4kdQO/yHZnAeN+I5sJ4H4UT6P0xf0ononqy+9GRGu1B+ouCMpJ6hjsZsz1xP0onkbpi/tRPEXsiy8NmZklzkFgZpa4eg+CNRNdwDhxP4qnUfrifhRP4fpS12MEZmZ28ur9jMDMzE6Sg8DMLHGFDgJJT0raKWm7pI6Kx/5OUkiaVXbso5L2SPqppCtrX3F1g/VD0oeyWndJWlV2vJD9gOp9kTRf0uaBY5IWlLUvZF8kTZd0l6THJe2W9IeSfkvSw5J+lv05o6x9PfXjU9l+p6RvSZpe1r6Q/YDqfSl7rJ5e71X7UejXe0QU9gt4EphV5fjZwHcofbBsVnbsVcAOoAU4F/gfYNJE92GwfgCXAf8JtGT7Zxa9H0P05bvAVdn21cDGovcFWAtcl22fAkwHVgEfyY59BLilTvtxBTA5O3ZLPfRjsL5k2/X2eq/2Oyn0673QZwRD+CywAigf6b4GuDMieiPiCWAPsKDaNxfEB4BPRkQvQETsz47XWz+g9Hs4Pds+A9iXbReyL5JOB94AfBkgIl6MiEOU6l2bNVsLvC3brqt+RMR3I+Jo1mwz0JZtF7IfMOTvBOro9T5EPwr9ei96EATwXUnbJC0HkLQEeDoidlS0nQPsLdvvyo4VwUv6AZwPvF7SFknfk/Sa7HiR+wHV+/LXwKck7QU+DXw0O17Uvvwe0A18VdKPJd0u6WXAWRHxDED255lZ+3rrR7m/ADZk20XtBwzSlzp8vQ/2Oyn0673oQbAwIi4FrgJukPQG4GPAx6u0VZVjRZkbW60fk4EZwOuADwPrJIli9wOq9+UDwE0RcTZwE9n/hihuXyYDlwK3RcQlwPOULgUNpi77IeljwFHgjoFDVZ6jCP2A6n35J+rv9T7Y76TQr/dCB0FE7Mv+3A98C/hjStfRdkh6ktIp76OSfptSkp5d9u1tHL9EMaGq9GMBpXrvjpKtQD+lxagK2w8YtC/XAndnTb7J8VPbovalC+iKiC3Z/l2UXrz/J2k2QPbn/rL29dQPJF0LvAV4d2QXoyluP2DwvtTb632wfhT69V7YIMhOC08b2KY0APajiDgzIuZGxFxKf4mXRsT/AvcB75DUIulc4BXA1gkq/zcG6cdjwD3Aouz4+ZQGlZ6loP2AIfuyj1JIQ6lPP8u2C9mX7N/LXkkXZIfeBPyEUr3XZseuBe7NtuuqH5IWAyuBJRHxQtm3FLIfMGhfHq231/sQ/7YK/XqfXOsfOApnAd8qnT0xGfj3iHhosMYRsUvSOkp/6UeBGyLiWE0qHVrVfkg6BfiKpMeAF4Frs/+5FbUfMHhfeoDPS5oMHAGWQ6F/JwAfAu7Ifg8/B5ZR+o/ROknvA54C/gzqsh8/ojQL5eHsd7U5Iq4veD+gel+qKnhfqvXjeQr8evcSE2ZmiSvspSEzM6sNB4GZWeIcBGZmiXMQmJklzkFgZpY4B4HZELKpsZXHLpC0UaXVVndLWlPx+OclPS2pKdtflrXdLulFHV+99ZO16ofZUDx91GwIknoiYlrFse8At0bEvdn+H0TEzmy7idIKrfsorWS6seJ7nwTaI+LZ/Ks3GxmfEZiN3mxKn3IFYCAEMpdR+rT1bcA7a1yX2Zg4CMxG77PAI5I2SLpJZTd+ofTm/x+U1mF6i6TmCanQbBQcBGajFBFfBS6ktMDeG4HN2Voxp1C6Mc89EfFrYAul9ZjMCq3Iaw2ZFVa2CutXOL5+zEWU1pE/A9iZrfFzKvAC8O2JqtNsJHxGYDZKkhYPXPLJlkSeCTxN6bLQdWWrZZ4LXCHp1Akr1mwEfEZgNrRTJXWV7X+G0prxn5d0JDv2YeDXwJXA+wcaRsTzkjYBbwW+UaN6zUbN00fNzBLnS0NmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWuP8HzSk5MVyKLGUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from tqdm import tqdm_notebook\n",
    "%matplotlib inline\n",
    "\n",
    "df = pd.DataFrame({'LSAT': LSAT, 'GPA': GPA})\n",
    "df.plot.scatter(x='LSAT', y='GPA')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated correlation coefficient: 0.5459\n"
     ]
    }
   ],
   "source": [
    "# Plug-in estimates for mean and correlation\n",
    "\n",
    "X = df['LSAT'].to_numpy()\n",
    "Y = df['GPA'].to_numpy()\n",
    "\n",
    "def corr(X, Y):\n",
    "    mu_x = X.mean()\n",
    "    mu_y = Y.mean()\n",
    "\n",
    "    return sum((X - mu_x) * (Y - mu_y)) / math.sqrt(sum((X - mu_x)**2) * sum((Y - mu_y)**2))\n",
    "  \n",
    "\n",
    "theta_hat = corr(X, Y)\n",
    "    \n",
    "print('Estimated correlation coefficient: %.4f' % corr(X, Y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "30559d3db0164669aecf51b52b625f10",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Estimated SE of correlation coefficient: 0.2674\n"
     ]
    }
   ],
   "source": [
    "# Bootstrap for SE of correlation coefficient\n",
    "\n",
    "nx = len(X)\n",
    "ny = len(Y)\n",
    "\n",
    "B = 1000000\n",
    "t_boot = np.empty(B)\n",
    "for i in tqdm_notebook(range(B)):\n",
    "    xx = np.random.choice(X, nx, replace=True)\n",
    "    yy = np.random.choice(Y, ny, replace=True)\n",
    "    t_boot[i] = corr(xx, yy)\n",
    "    \n",
    "se = t_boot.std()\n",
    "\n",
    "print('Estimated SE of correlation coefficient: %.4f' % se)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% confidence interval (Normal): \t 0.022, 1.070\n",
      "95% confidence interval (percentile): \t -0.503, 0.522\n",
      "95% confidence interval (pivotal): \t 0.569, 1.594\n"
     ]
    }
   ],
   "source": [
    "# Confidence intervals obtained from bootstrap\n",
    "\n",
    "from scipy.stats import norm\n",
    "\n",
    "z = norm.ppf(.975)\n",
    "\n",
    "normal_conf = (theta_hat - z * se, theta_hat + z * 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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 9.6.2**.  (Computer Experiment).  Conduct a simulation to compare the four bootstrap confidence interval methods.\n",
    "\n",
    "Let $n = 50$ and let $T(F) = \\int (x - \\mu)^3 dF(x) / \\sigma^3$ be the skewness.  Draw $Y_1, \\dots, Y_n \\sim N(0, 1)$ and set $X_i = e^{Y_i}$, $i = 1, \\dots, n$.  Construct the four types of bootstrap 95% intervals for $T(F)$ from the data $X_1, \\dots, X_n$.  Repeat this whole thing many times and estimate the true coverage of the four intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from tqdm import tqdm_notebook\n",
    "from scipy.stats import norm\n",
    "\n",
    "def create_data(n=50):\n",
    "    y = norm.rvs(size=n)\n",
    "    return np.exp(y)\n",
    "\n",
    "def skewness(x):\n",
    "    n = len(x)\n",
    "    mu = sum(x) / n\n",
    "    var = sum((x - mu)**2) / n\n",
    "    return sum((x - mu)**3 / (n * var**(3/2)))\n",
    "\n",
    "def bootstrap_values(x, B=10000, show_progress=True):\n",
    "    n = len(x)\n",
    "    t_boot = np.empty(B)\n",
    "    iterable = tqdm_notebook(range(B)) if show_progress else range(B)\n",
    "    for i in iterable:\n",
    "        xx = np.random.choice(x, n, replace=True)\n",
    "        t_boot[i] = skewness(xx)\n",
    "\n",
    "    return t_boot\n",
    "\n",
    "def bootstrap_intervals(theta_hat, t_boot, alpha=0.05):\n",
    "    se = t_boot.std()\n",
    "    \n",
    "    z = norm.ppf(1 - alpha/2)\n",
    "    q_half_alpha = np.quantile(t_boot, alpha/2)\n",
    "    q_c_half_alpha = np.quantile(t_boot, 1 - alpha/2)\n",
    "\n",
    "    normal_conf = (theta_hat - z * se, theta_hat + z * se)\n",
    "    percentile_conf = (q_half_alpha, q_c_half_alpha)\n",
    "    pivotal_conf = (2*theta_hat - q_c_half_alpha, 2*theta_hat - q_half_alpha)\n",
    "    \n",
    "    return normal_conf, percentile_conf, pivotal_conf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating the data\n",
    "x = create_data(n=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "654583f1040942f28c910f15e43ba0e8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=100000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "95% confidence interval (Normal): \t 1.032, 2.638\n",
      "95% confidence interval (percentile): \t 1.154, 2.757\n",
      "95% confidence interval (pivotal): \t 0.912, 2.515\n"
     ]
    }
   ],
   "source": [
    "# Nonparametric Bootstrap\n",
    "theta_hat = skewness(x)\n",
    "t_boot = bootstrap_values(x, B=100000)\n",
    "\n",
    "normal_conf, percentile_conf, pivotal_conf = bootstrap_intervals(theta_hat, t_boot, alpha=0.05)\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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: parametric bootstrap is only covered in chapter 10.  The below assumes that \"four types of bootstrap\" in the exercise refers to the 3 types of nonparametric bootstrap covered in chapter 9, plus parametric bootstrap from chapter 10.\n",
    "\n",
    "For the parametric bootstrap, assume $X = e^Y$ where $Y \\sim N(\\mu, \\sigma^2)$.  Then\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{E}(X) & = \\mathbb{E}(e^Y) = \\int_{-\\infty}^\\infty e^y \\frac{1}{\\sigma \\sqrt{2 \\pi}} e^{-\\frac{1}{2} \\left(\\frac{y - \\mu}{\\sigma} \\right)^2} dy = \\frac{1}{\\sigma \\sqrt{2 \\pi}} \\int e^{y - \\frac{1}{2}\\left(\\frac{y - \\mu}{\\sigma}\\right)^2} dy = e^{\\mu + \\sigma^2/2} \\\\\n",
    "\\mathbb{E}(X^2) &= \\mathbb{E}(e^{2Y}) = e^{2\\mu + (2\\sigma)^2/2} = e^{2\\mu + 2\\sigma^2}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Solving the two equations for the moments we get parameter estimates:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{\\mu}_Y & = 4 \\log \\mathbb{E}(X) - \\log \\mathbb{E}(X^2) \\\\\n",
    "\\hat{\\sigma}_Y & = \\sqrt{\\log \\mathbb{E}(X^2) - 2 \\log \\mathbb{E}(X)}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "From these parameter estimates, we can generate bootstrap samples:  sample values for $Y_i ~ N(\\hat{\\mu}_Y, \\hat{\\sigma}_Y^2)$, calculate $X_i = Y_i$, and calculate the statistic $T(X_i)$ on each sample generated this way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parametric bootstrap\n",
    "# Assume X = e^Y, Y ~ N(\\mu, \\sigma^2)\n",
    "\n",
    "def estimate_parameters(x):\n",
    "    log_e_x = np.log(x.mean())\n",
    "    log_e_x2 = np.log((x**2).mean())\n",
    "    \n",
    "    mu_Y_hat = 4 * log_e_x - log_e_x2\n",
    "    sigma_Y_hat = np.sqrt(log_e_x2 - 2 * log_e_x)\n",
    "    \n",
    "    return mu_Y_hat, sigma_Y_hat\n",
    "\n",
    "def bootstrap_skewness_parametric(mu, sigma, n, B = 10000, show_progress=True):\n",
    "    t_boot = np.empty(B)\n",
    "    iterable = tqdm_notebook(range(B)) if show_progress else range(B)\n",
    "    for i in iterable:\n",
    "        xx = np.exp(norm.rvs(size=n, loc=n, scale=sigma))\n",
    "        t_boot[i] = skewness(xx)\n",
    "        \n",
    "    return t_boot\n",
    "\n",
    "def bootstrap_parametric_intervals(mu, t_boot, alpha=0.05):\n",
    "    se = t_boot.std()\n",
    "    z = norm.ppf(1 - alpha/2)\n",
    "    \n",
    "    return (theta_hat - z * se, theta_hat + z * se)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dfaab0bebe074b99aec5ca4bc485566b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=100000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "95% confidence interval (Parametric Normal): \t -0.181, 3.851\n"
     ]
    }
   ],
   "source": [
    "mu_X_hat = x.mean()\n",
    "mu_Y_hat, sigma_Y_hat = estimate_parameters(x)\n",
    "t_boot = bootstrap_skewness_parametric(mu_Y_hat, sigma_Y_hat, 50, B=100000)\n",
    "\n",
    "parametric_normal_conf = bootstrap_parametric_intervals(mu_X_hat, t_boot, alpha=0.05)\n",
    "print('95%% confidence interval (Parametric Normal): \\t %.3f, %.3f' % parametric_normal_conf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "df67a6db578c4ca9980c1359eb3ea0c9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=10), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# Repeat nonparametric bootstrap many times\n",
    "n_experiments = 10\n",
    "\n",
    "normal_conf = np.empty((n_experiments, 2))\n",
    "percentile_conf = np.empty((n_experiments, 2))\n",
    "pivotal_conf = np.empty((n_experiments, 2))\n",
    "\n",
    "for i in tqdm_notebook(range(n_experiments)):\n",
    "    theta_hat = skewness(x)\n",
    "    t_boot_experiment = bootstrap_values(x, B=100000, show_progress=False)\n",
    "    normal_conf[i], percentile_conf[i], pivotal_conf[i] = bootstrap_intervals(theta_hat, t_boot_experiment, alpha=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Normal confidence lower bound: \t\t mean 1.031, SE 0.001\n",
      "Normal confidence upper bound: \t\t mean 2.639, SE 0.001\n",
      "Percentile confidence lower bound: \t mean 1.154, SE 0.002\n",
      "Percentile confidence upper bound: \t mean 2.762, SE 0.005\n",
      "Pivotal confidence lower bound: \t mean 0.907, SE 0.005\n",
      "Pivotal confidence upper bound: \t mean 2.515, SE 0.002\n"
     ]
    }
   ],
   "source": [
    "print('Normal confidence lower bound: \\t\\t mean %.3f, SE %.3f' % \n",
    "      (normal_conf[:, 0].mean(), normal_conf[:, 0].std()))\n",
    "print('Normal confidence upper bound: \\t\\t mean %.3f, SE %.3f' % \n",
    "      (normal_conf[:, 1].mean(), normal_conf[:, 1].std()))\n",
    "\n",
    "print('Percentile confidence lower bound: \\t mean %.3f, SE %.3f' % \n",
    "      (percentile_conf[:, 0].mean(), percentile_conf[:, 0].std()))\n",
    "print('Percentile confidence upper bound: \\t mean %.3f, SE %.3f' % \n",
    "      (percentile_conf[:, 1].mean(), percentile_conf[:, 1].std()))\n",
    "\n",
    "print('Pivotal confidence lower bound: \\t mean %.3f, SE %.3f' % \n",
    "      (pivotal_conf[:, 0].mean(), pivotal_conf[:, 0].std()))\n",
    "print('Pivotal confidence upper bound: \\t mean %.3f, SE %.3f' % \n",
    "      (pivotal_conf[:, 1].mean(), pivotal_conf[:, 1].std()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9277970ebba6491a927fd789915562c9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=10), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# Repeat parametric bootstrap many times\n",
    "\n",
    "n_experiments = 10\n",
    "parametric_normal_conf = np.empty((n_experiments, 2))\n",
    "\n",
    "for i in tqdm_notebook(range(n_experiments)):\n",
    "    mu_X_hat = x.mean()\n",
    "    mu_Y_hat, sigma_Y_hat = estimate_parameters(x)\n",
    "    t_boot = bootstrap_skewness_parametric(mu_Y_hat, sigma_Y_hat, 50, B=100000, show_progress=False)\n",
    "    parametric_normal_conf[i] = bootstrap_parametric_intervals(mu_X_hat, t_boot, alpha=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parametric Normal confidence lower bound: \t mean -0.185, SE 0.003\n",
      "Parametric Normal confidence upper bound: \t mean 3.854, SE 0.003\n"
     ]
    }
   ],
   "source": [
    "print('Parametric Normal confidence lower bound: \\t mean %.3f, SE %.3f' % \n",
    "      (parametric_normal_conf[:, 0].mean(), parametric_normal_conf[:, 0].std()))\n",
    "print('Parametric Normal confidence upper bound: \\t mean %.3f, SE %.3f' % \n",
    "      (parametric_normal_conf[:, 1].mean(), parametric_normal_conf[:, 1].std()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 9.6.3**.  Let  $X_1, \\dots, X_n \\sim t_3$ where $n = 25$.  Let $\\theta = T(F) = (q_{.75} - q_{.25})/1.34$ where $q_p$ denotes the $p$-th quantile.  Do a simulation to compare the coverage and length of the following confidence intervals for $\\theta$:\n",
    "\n",
    "- Normal interval with standard error from the bootstrap\n",
    "- Bootstrap percentile interval\n",
    "\n",
    "Remark: the jacknife does not give a consistent estimator of the variance of a quantile."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.  We will assume that $t_3$ represents a t-distribution with shape parameter 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import t\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "n = 25\n",
    "X = t.rvs(3, size=25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def T(x):\n",
    "    return (np.quantile(x, 0.75) - np.quantile(x, 0.25)) / 1.34"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_hat = T(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "47ce463065f9482cafbc77f38b9a08a6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "95% confidence interval (Normal): \t 0.032, 2.089\n",
      "95% confidence interval (percentile): \t 0.605, 2.378\n"
     ]
    }
   ],
   "source": [
    "# Run bootstrap\n",
    "\n",
    "B = 1000000\n",
    "t_boot = np.empty(B)\n",
    "for i in tqdm_notebook(range(B)):\n",
    "    xx = np.random.choice(X, n, replace=True)\n",
    "    t_boot[i] = T(xx)\n",
    "    \n",
    "se_boot = t_boot.std()\n",
    "\n",
    "alpha = 0.05\n",
    "z = norm.ppf(1 - alpha/2)\n",
    "q_half_alpha = np.quantile(t_boot, alpha/2)\n",
    "q_c_half_alpha = np.quantile(t_boot, 1 - alpha/2)\n",
    "\n",
    "normal_conf = (theta_hat - z * se_boot, theta_hat + z * se_boot)\n",
    "percentile_conf = (q_half_alpha, q_c_half_alpha)\n",
    "\n",
    "print('95%% confidence interval (Normal): \\t %.3f, %.3f' % normal_conf)\n",
    "print('95%% confidence interval (percentile): \\t %.3f, %.3f' % percentile_conf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 9.6.4**.  Let $X_1, \\dots, X_n$ be distinct observations (no ties).  Show that there are\n",
    "\n",
    "$$ \\binom{2n - 1}{n}$$\n",
    "\n",
    "distinct bootstrap samples.\n",
    "\n",
    "Hint: Imagine putting $n$ balls into $n$ buckets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "Each bootstrap sample (random draws with replacement) will select $a_i$ copies of $X_i$, where $0 \\leq a_i \\leq n$ and $\\sum_{i=1}^n a_i = n$, for integer $a_i$.  Each bootstrap sample is uniquely represented by this sequence of variables, and each sequence of variables uniquely determines a bootstrap sample -- so the number of distinct bootstrap samples is equal to the number of solutions to this equation, that is, the number of partitions of $n$ into $n$ buckets.\n",
    "\n",
    "Lets write $a_i$ explicitly in base 1, representing it as $a_i$ consecutive copies of the digit 1:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "0_{10} & = \\text{empty string}_1 \\\\\n",
    "1_{10} & = 1_1 \\\\\n",
    "2_{10} & = 11_1 \\\\\n",
    "3_{10} & = 111_1 \\\\\n",
    "4_{10} & = 1111_1 \\\\\n",
    "\\vdots\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "So a solution for\n",
    "\n",
    "$$a_1 + a_2 + \\dots + a_n = n$$\n",
    "\n",
    "is uniquely represented by a sequence of $2n - 1$ symbols, being $n$ digits $1$ and $n - 1$ plus signs.  For example, if $a_1 = 3$, $a_2 = 0$, $a_3 = 1$, then we write\n",
    "\n",
    "$$ 111 + + 1 + \\dots = n $$\n",
    "\n",
    "The number of such solutions is then obtained by choosing $n$ of the $2n - 1$ symbols to be digit $1$ -- that is, $\\binom{2n - 1}{n}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 9.6.5**.  Let $X_1, \\dots, X_n$ be distinct observations (no ties).  Let $X_1^*, \\dots, X_n^*$ denote a bootstrap sample and let $\\overline{X}_n^* = n^{-1}\\sum_{i=1}^nX_i^*$.  Find:\n",
    "\n",
    "- $\\mathbb{E}(\\overline{X}_n^* | X_1, \\dots, X_n)$\n",
    "- $\\mathbb{V}(\\overline{X}_n^* | X_1, \\dots, X_n)$\n",
    "- $\\mathbb{E}(\\overline{X}_n^*)$\n",
    "- $\\mathbb{V}(\\overline{X}_n^*)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "**(a)** $$\n",
    "\\begin{align}\n",
    "\\mathbb{E}(\\overline{X}_n^* | X_1, \\dots, X_n) &= \\mathbb{E}\\left(n^{-1}\\sum_{i=1}^nX_i\\right) = n^{-1}\\sum_{i=1}^n \\mathbb{E}(X_i) = \\mathbb{E}(X)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "**(b)** $$\n",
    "\\begin{align}\n",
    "\\mathbb{V}(\\overline{X}_n^* | X_1, \\dots, X_n) & =\n",
    "\\mathbb{E}((\\overline{X}_n^*)^2 | X_1, \\dots, X_n) - \\mathbb{E}(\\overline{X}_n^* | X_1, \\dots, X_n)^2 \\\\\n",
    "&= \\mathbb{E}\\left(\\left(n^{-1}\\sum_{i=1}^nX_i\\right)^2\\right) - \\mathbb{E}\\left(n^{-1}\\sum_{i=1}^nX_i\\right)^2 \\\\\n",
    "&= n^{-2} \\mathbb{E}\\left(\\sum_{i=1}^nX_i^2 + \\sum_{i=1}^n \\sum_{j=1, j \\neq i}^n X_i X_j\\right) - \\mathbb{E}(X)^2 \\\\\n",
    "&= n^{-2} \\sum_{i=1}^n \\mathbb{E}(X_i^2) + n^{-2} \\sum_{i=1}^n \\sum_{j=1, j \\neq i}^n \\mathbb{E}(X_i X_j) - \\mathbb{E}(X)^2  \\\\\n",
    "&= n^{-1} (\\mathbb{V}(X) + \\mathbb{E}(X)^2) + n^{-2} \\sum_{i=1}^n \\sum_{j=1, j \\neq i}^n \\mathbb{E}(X_i) \\mathbb{E}(X_j) - \\mathbb{E}(X)^2 \\\\\n",
    "&= n^{-1} (\\mathbb{V}(X) + \\mathbb{E}(X)^2) + n^{-1} (n - 1) \\mathbb{E}(X)^2 - \\mathbb{E}(X)^2 \\\\\n",
    "&= \\frac{1}{n}\\mathbb{V}(X)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "**(c)** $$\n",
    "\\begin{align}\n",
    "\\mathbb{E}(\\overline{X}_n^*) &= \\mathbb{E}\\left(n^{-1} \\sum_{i=1}^nX_i^*\\right) = n^{-1} \\sum_{i=1}^n \\mathbb{E}(X_i^*) = \\mathbb{E}(X)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "**(d)** $$\n",
    "\\begin{align}\n",
    "\\mathbb{V}(\\overline{X}_n^*) &= \\mathbb{E}((\\overline{X}_n^*)^2) - \\mathbb{E}(\\overline{X}_n^*)^2 \\\\\n",
    "&= \\mathbb{E}\\left(\\left(n^{-1} \\sum_{i=1}^nX_i^*\\right)^2\\right) - \\mathbb{E}(X)^2 \\\\\n",
    "&= n^{-2} \\sum_{i=1}^n\\sum_{j=1}^n \\mathbb{E}(X_i^* X_j^*) - \\mathbb{E}(X)^2\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Now, the same value $X_k$ may have been sampled twice in $\\mathbb{E}(X_i^* X_j^*)$.  This always happens when $i = j$, and this happens with probability $1 / n$ when $i \\neq j$. Thus,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{V}(\\overline{X}_n^*) \n",
    "&= n^{-2} \\left(\\sum_{i=1}^n \\mathbb{E}(X_i^2) + \\sum_{i=1}^n \\sum_{j=1, j \\neq i}^n \\left( \\frac{1}{n} \\mathbb{E}(X_i^2) + \\left(1 - \\frac{1}{n}\\right)\\mathbb{E}(X_i)\\mathbb{E}(X_j)\\right) \\right) - \\mathbb{E}(X)^2 \\\\\n",
    "&= n^{-1} \\mathbb{E}(X^2) + n^{-2} (n-1) \\mathbb{E}(X^2) + n^{-2}(n-1)^2 E(X)^2 - \\mathbb{E}(X)^2 \\\\\n",
    "&= n^{-2} (2n - 1) \\left( \\mathbb{E}(X^2) - \\mathbb{E}(X)^2 \\right) \\\\\n",
    "&= \\frac{2n - 1}{n^2} \\mathbb{V}(X)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 9.6.6**.  (Computer Experiment).  Let $X_1, \\dots, X_n \\sim N(\\mu, 1)$.  Let $\\theta = e^\\mu$and let $\\hat{\\theta} = e^{\\overline{X}}$ be the mle.  Create a dataset (using $\\mu = 5$) consisting of $n = 100$ observations.\n",
    "\n",
    "**(a)** Use the bootstrap to get the $\\text{se}$ and 95% confidence interval for $\\theta$.\n",
    "\n",
    "**(b)** Plot a histogram of the bootstrap replications for the parametric and non-parametric bootstraps.  These are estimates of the distribution of $\\hat{\\theta}$.  Compare this to the true sampling distribution of $\\hat{\\theta}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import norm\n",
    "from tqdm import tqdm_notebook\n",
    "%matplotlib inline\n",
    "\n",
    "X = norm.rvs(loc=5, scale=1, size=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1969eb4d8b3d414893d01ef9ee29750a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "95% confidence interval (Normal): \t 104.486, 160.228\n"
     ]
    }
   ],
   "source": [
    "# Get estimated value\n",
    "\n",
    "theta_hat = np.exp(X.mean())\n",
    "\n",
    "# Run nonparametric bootstrap\n",
    "\n",
    "B = 1000000\n",
    "t_boot_nonparam = np.empty(B)\n",
    "n = len(X)\n",
    "for i in tqdm_notebook(range(B)):\n",
    "    xx = np.random.choice(X, n, replace=True)\n",
    "    t_boot_nonparam[i] = np.exp(xx.mean())\n",
    "    \n",
    "se_boot = t_boot_nonparam.std()\n",
    "\n",
    "alpha = 0.05\n",
    "z = norm.ppf(1 - alpha/2)\n",
    "normal_conf = (theta_hat - z * se_boot, theta_hat + z * se_boot)\n",
    "\n",
    "print('95%% confidence interval (Normal): \\t %.3f, %.3f' % normal_conf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7be0c785288b4815b277862cb8c0a243",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "95% confidence interval (Parametric Normal): \t 106.225, 158.490\n"
     ]
    }
   ],
   "source": [
    "# Run parametric bootstrap\n",
    "\n",
    "mu_hat = X.mean()\n",
    "theta_hat = np.exp(mu_hat)\n",
    "\n",
    "B = 1000000\n",
    "t_boot_param = np.empty(B)\n",
    "n = len(X)\n",
    "for i in tqdm_notebook(range(B)):\n",
    "    xx = norm.rvs(size=n, loc=mu_hat, scale=1)\n",
    "    t_boot_param[i] = np.exp(xx.mean())\n",
    "    \n",
    "se_boot_param = t_boot_param.std()\n",
    "\n",
    "alpha = 0.05\n",
    "z = norm.ppf(1 - alpha/2)\n",
    "normal_conf_param = (theta_hat - z * se_boot_param, theta_hat + z * se_boot_param)\n",
    "\n",
    "print('95%% confidence interval (Parametric Normal): \\t %.3f, %.3f' % normal_conf_param)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the true sampling distribution,\n",
    "\n",
    "$$ \\overline{X} = \\frac{1}{n} \\sum_{i=1}^n X_i \\sim \\frac{1}{n} N(n \\mu, n) = N(\\mu, n^{-2})$$\n",
    "\n",
    "so the distribution of $\\hat{\\theta}$ is the distribution of $e^\\overline{X}$.  Its CDF is:\n",
    "\n",
    "$$\n",
    "\\mathbb{P}\\left(\\hat{\\theta} \\leq t\\right) = \n",
    "\\mathbb{P}\\left(e^\\overline{X} \\leq t\\right) = \n",
    "\\mathbb{P}\\left(\\overline{X} \\leq \\log t\\right) = F_{\\overline{X}}\\left(\\log t\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "from matplotlib import pyplot\n",
    "\n",
    "bins = np.linspace(50, 250, 500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate the CDF for theta, calculate it for each bin, and include the differences between bins\n",
    "\n",
    "def theta_cdf(x):\n",
    "    return norm.cdf(np.log(x), loc=5, scale=1/50)\n",
    "\n",
    "theta_cdf_bins = list(map(theta_cdf, bins))\n",
    "theta_cdf_bins_delta = np.empty(len(bins))\n",
    "theta_cdf_bins_delta[0] = 0\n",
    "theta_cdf_bins_delta[1:] = np.diff(theta_cdf_bins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAHSCAYAAAAezFYoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdfXRU1dn+8WsTECjgC4hCjc6ARQp5YUhCCAWZICggAoq2liiUtjqK2lZXq6XWVrTY2sJ69OeqFkesaDVVa2vlaa3VVhLfoEBoBMT6iDAjUVDAgiAGSdi/P5IZJyeTZAKTzNv3s1aWmTN7zuw5RteVO/fZ21hrBQAAAOBzXRI9AQAAACDZEJIBAAAAB0IyAAAA4EBIBgAAABwIyQAAAIADIRkAAABw6JroCURz8sknW7fbnehpAEBSeuuttyRJQ4cOTfBMACC1VVVV7bbW9o/2XFKGZLfbrXXr1iV6GgCQlEpLSyVJFRUVCZ0HAKQ6Y0ywpedotwAAAAAcCMkAAACAAyEZAAAAcEjKnuRoDh8+rJqaGtXW1iZ6KkBYjx49lJ2drW7duiV6KgAAII5SJiTX1NSoT58+crvdMsYkejqArLXas2ePampqNGjQoERPBwAAxFHKtFvU1taqX79+BGQkDWOM+vXrx183AABIQykTkiURkJF0+JkEACA9pVRITpQ9e/bI4/HI4/FowIABOu2008KPP/vss0RP76jU1dXpxBNPlCRt375dl156aYe8T3Z2tvbu3av6+nqdffbZLY47cuSI7rzzzlbPNXnyZO3fv19btmyRx+Np1zxefPFFrV69Ovz43nvv1WOPPdaucwAAgMyRMj3JidSvXz9VV1dLkhYuXKjevXvrBz/4QZMx1lpZa9WlS+r93nH66afriSee6ND3yMrK0ssvv9zi86GQvGDBgmbPha7t3//+d0nSBx980O73f/HFF3XyySerpKREknTttde2+xwAACBzpF6iSyJbtmxRbm6urr76ahUUFGj79u3h6qwkPf7447riiiskNQS7WbNmqaioSMXFxU2qmiEbN27UqFGj5PF4lJ+fr61bt0qSpk+frsLCQuXk5GjZsmWSPq8E33jjjSooKNDkyZP1r3/9S16vV4MHD9azzz4rSVq2bJkuuugiTZ48WUOHDtWiRYuifo5QZXbZsmW65JJLNHnyZA0ZMkQ/+tGPwuPuv/9+nXXWWSotLdUVV1yh66+/vtm5du3apXPPPVcFBQWaP3++rLVN5itJ7733nsaNGyePx6Pc3Fy99tprWrBggfbv3y+Px6O5c+c2u7Y7duwIV6WlhtVO5syZo7y8PH3ta1/Tp59+KklNxqxevVqTJk3SO++8o2XLlmnx4sXyeDx67bXXdMstt+juu++WJK1fv16jR49Wfn6+Lr74Yu3bt0+SNG7cOC1YsEDFxcUaOnSoXnvttRh+KgAAQDpIyUry9c9dr+qd1XE9p2eAR3dPubvdr9u8ebMeeughLV26VHV1dS2O++53v6ubbrpJJSUlCgQCuuCCC7Rp06YmY+677z794Ac/0KWXXqpDhw6FA+bDDz+svn376uDBgyoqKtLFF1+sPn36aN++fTrvvPO0ePFiTZ8+XQsXLtQ///lPvf7667rqqqt0/vnnS5LWrFmjTZs26bjjjtOoUaN0wQUXKDc3t8W5vv7661q/fr26du2qs846S9/5zndUX1+vO++8U+vXr1evXr1UWlqq4uLiZq+99dZbNWHCBN1888165plntHTp0mZjHn30UU2fPl0//OEPVV9fr08//VTFxcVatmxZuGK/ZcuWJtc22nV/8MEHVVJSorlz5+r++++PGtol6cwzz9QVV1yhk08+OTwm9EuEJF1++eXy+/0aN26cbr75Zv3sZz/TkiVLJDVUsdesWaMVK1bo9ttv13PPPdfidQMAAOkjJUNyMjnzzDM1atSoNsf94x//0FtvvRV+/N///leffvqpevbsGT72la98RYsWLVIwGNSsWbP0pS99SZJ01113acWKFZIalsJ755135PF41LNnT5177rmSpLy8PJ1wwgnq2rWr8vLyFAgEwuedPHmyTjrpJEnShRdeqFdeeaXVkDxp0iT16dNHkvTlL39Z7777rmpqanTOOeeEz3PJJZfo3Xffbfbal156KRxAZ86cGT5PpFGjRumqq65SbW2tLrzwQo0YMSLqLxitXdtBgwaFWydCIbelkNyaPXv2qLa2VuPGjZMkfeMb39CcOXPCz8+aNUuSVFhY2OSaAgCA9JaSIfloKr4dpVevXuHvu3TpEq7+SmqyNFioInnccce1eK45c+ZozJgx+utf/6pzzz1XDz/8sD777DO99NJLWr16tXr27Klx48aFzxt5ri5duqh79+7h7yNDp3MFhrZWZAidR2roJa6rq2vyudrS1vnPOeccVVRU6K9//asuu+wy/ehHP4p642DktW3rPUKPu3btqiNHjkhSTEuztfW5QtcidB0AAEBmoCc5jrp06aKTTjpJb7/9to4cOaKnn346/NykSZN07733hh+H2goibd26VV/60pf0ve99T9OmTdOGDRu0b98+9e3bVz179tQbb7yhtWvXtntezz//vPbu3auDBw/qmWee0dixY9t9jtGjR2vlypXau3evDh8+rD/96U9Rx40fPz68asT//u//av/+/c3GBINBDRgwQD6fT/PmzdO///1vde3a8PtarEF027Zt4Wvx+9//PlwJdrvdqqqqkiT98Y9/DI/v06dP1LmcfPLJ6tmzZ7jf+He/+528Xm9McwAAAOmLkBxnv/zlLzVlyhRNnDhR2dnZ4eP33nuvXn31VeXn52v48OF64IEHmr22vLxcOTk58ng82rp1qy6//HJNmzZNBw8e1IgRI3T77bdr9OjR7Z7TuHHjVFZWppEjR2r27NntXj5Nks444wzdeOONKi4u1nnnnaecnBydcMIJzcbddttt+sc//qGCggJVVFTotNNOazbmn//8p0aMGKGRI0fqmWee0Xe+8x1J0re//W3l5+dr7ty5bc4nJydHDzzwgPLz8/XJJ5/I5/NJalh95JprrtHZZ5/dpNI+c+ZMPfnkkxo5cmSzG/B+97vf6YYbblB+fr42b96sW265pV3XBgAApB/Tnj+jd5aioiK7bt26JsfefPNNDRs2LEEzSl3Lli3Tpk2bwis5HIsDBw6od+/eOnz4sGbOnKn58+dr+vTpcZhlauNnE52ttLRUklRRUZHQeQBAqjPGVFlri6I9RyUZMfvJT36ikSNHKj8/X0OHDtUFF1yQ6CkBAAB0iJS8cQ+xC63THA933XVX3M4FAACQzAjJAJBC/FV+Ve+s1qm9Tk30VAAgrRGSASCFlG8s177afYmeBgCkPXqSAQAAAAdCMgCkCH+VX5XBykRPAwAyAiG5HZ588kkNGzZMw4YN0x/+8IeoY6y1Ouecc/Txxx9r7969uu+++zp5lol3LJ970qRJ+u9//xvnGQHpoXxjefj7fbX75K/yJ3A2AJDeCMkx2r9/v375y19q7dq1WrNmjX7xi1/owIEDzcY9++yzGjFihI4//vhWw2J9fX1HT7lV1trw9s3xdjQhOTSfOXPmZOQvFkCsvC6vzup3lqSmoRkAEF+E5Bg999xzmjp1qnr37q0+ffpo0qRJ+tvf/tZs3GOPPaaZM2dKkhYsWKB33nlHHo9HN954oyoqKjRhwgSVlZUpLy9PgUBAubm54dcuWbJECxculCS98847mjJligoLC3X22WfrP//5T7P3WrhwoebMmaNzzjlHQ4YMCe/id+DAAU2cOFEFBQXKy8vTM888I0kKBAIaNmyYrrnmGhUUFGj79u2aP3++ioqKlJOTo1tvvTV8brfbrZtvvlljxoxRUVGR1q9fr8mTJ+vMM8/U0qVLw+MWL16sUaNGKT8/P/x65+duaVy0+cyYMUO///3vj/rfE5AJBvYZqBN6NN/xEgAQPym5usX111+v6urquJ7T4/G0uivdu+++q4ceekh/+ctfJEm7du3Sqac2X4Lp1Vdf1f333y9JuvPOO7Vp06bwXCsqKrRmzRpt2rRJgwYNUiAQaPH9fD6fli5dqiFDhuhf//qXrrnmGr344ovNxm3YsEGrV6/WJ598opEjR2ratGk65ZRT9PTTT+v444/X7t27VVJSohkzZkiS3nrrLT300EPhau0dd9yhvn37qr6+XhMnTtSGDRuUn58vSTr99NO1atUq3XDDDZo3b55effVV1dbWKicnR1dffbWef/55vf3221qzZo2stZoxY4ZeeumlZp+7pXFnnHFGs/lI0qFDh7Rnzx7169evxesDAADQkVIyJCeCtVbz58/XLbfcIkm6/fbbZYxpNu6jjz5Snz59WjxPcXGxBg0a1Op7HThwQK+99pq++tWvho8dOnQo6tiZM2eqZ8+e6tmzpyZMmKA1a9Zo2rRpuvnmm/XSSy+pS5cueu+99/TBBx9Iklwul0pKSsKvf/LJJ+X3+1VXV6cdO3Zo8+bN4ZAcCtZ5eXk6cOCA+vTpoz59+qhHjx7au3evnn/+eT3//PMaOXJkeN5vv/22zjjjjCZzbG2ccz6SdMopp+j9998nJAMAgIRJyZDcWsW3o7hcriaV3O3bt2vy5MnNxnXt2lVHjhxRly7RO1l69erVbGxIbW2tJOnIkSM68cQTY6qWO4O6MUaPPfaYdu3apaqqKnXr1k1utzt87sj337Ztm5YsWaK1a9fqpJNO0rx588LjJKl79+6SpC5duoS/Dz2uq6uTtVY/+tGPdNVVVzWZg7NC3tq4yPlEXoeePXu2+dkBAAA6Cj3JMZo6dapWrlypPXv2aPfu3aqoqNCUKVOajRs6dKi2bt0qSerTp4/279/f4jlPPfVUffjhh9qzZ48OHToUbuU4/vjjNWjQoPAKGtZavf7661HP8cwzz6i2tlZ79uxRRUWFRo0apX379umUU05Rt27dtHLlSgWDwaiv/fjjj9WrVy+dcMIJ+uCDD6L2WLdm8uTJ+u1vfxu+gfG9997Thx9+2OxztzQuGmutdu7cKbfb3a65AAAAxFNKVpIToXfv3lq0aJHGjBkja61+8YtfqHfv3s3GTZs2TRUVFfrSl76kfv36aezYscrNzdXUqVM1bdq0JmO7deumn/70pxo9erQGDRqkL3/5y+HnHnvsMc2fP1+LFi3S4cOH9fWvf10jRoxo9n7FxcWaNm2a3n33Xf3kJz/RF7/4RV122WWaPn26ioqK5PF4mpw30ogRIzRy5Ejl5ORo8ODBGjt2bLuuyXnnnac333xTY8aMCV+jRx99VGeeeWaTz7148eKo47Kyspqds6qqSiUlJeralR9NAACQOMZam+g5NFNUVGTXrVvX5Nibb76pYcOGJWhGsduxY4fmzp2rF154ocPfa+HCherdu7d+8IMfdPh7dZbvfe97mjFjhiZOnJjoqcQsVX42kfpKl5c2fLNcqt5ZLc8CjyrmVSRwRgCQ2owxVdbaomjP0W4RZwMHDtSVV16pjz/+ONFTSUm5ubkpFZCBzuCv8qt0eamqd8Z3VR8AQMv4m3YH+NrXvtYp7xNaUzmdXHnllYmeApB0yjeWN1SOB3hUllemcrGJCAB0NEIyAKQAz4DPWysIyQDQ8Wi3AIAUVb2zWqXLS+Wv8id6KgCQdqgkA0AKOrXXqRo4YGC4T9lX6EvwjAAgvVBJBoAUNLDPQFXMq5BngCfRUwGAtERIbgdjjL7//e+HHy9ZsiQtb567++67dfDgwRafv+KKK7R58+aYzrV8+XJdd911xzynvXv36r777mtz3M9//vNjfi8AAICUDclut2RM/L5i2eCte/fu+tOf/qTdu3d39MdrU319fYedu7WQXF9fr2XLlmn48OEd9v7RHGtIttY22QIcAACgNSkbkoNBydr4fbWwc3MTXbt2lc/n01133RVlPkFNnDhR+fn5mjhxot59911J0rx58/Td735XX/nKVzR48GA99dRTkqSKigqNHz9eF110kYYPH66rr746HOLmz5+voqIi5eTk6NZbbw2/h9vt1u23365x48bpD3/4gx544AGNGjVKI0aM0MUXXxwOtvPmzdP8+fM1YcIEDR48WJWVlfrWt76lYcOGad68eeHzPf/88xozZowKCgr01a9+VQcOHNA999yj999/XxMmTNCECRMkNeyQF9oZcNWqVSotLVVos5fnnntOBQUFGjFiRIvrG2/fvl1TpkzR0KFDddttt4WP/8///I9yc3OVm5uru+++u9XjCxYs0DvvvCOPx6Mbb7xRO3bs0Pjx4+XxeJSbm6uXX35ZCxYs0KeffiqPx6PLLrtMgUBAw4YN0zXXXKOCggJt37691Wv7wx/+UMXFxSouLtaWLVva/oEAAADpy1qbdF+FhYXWafPmzU0eS82GHJNYzterVy+7b98+63K57N69e+3ixYvtrbfeaq219oILLrDLly+31lr74IMP2pkzZ1prrf3GN75hL7nkEltfX2/feOMNe+aZZ1prrV25cqXt3r27feedd2xdXZ2dNGmS/cMf/mCttXbPnj3WWmvr6uqs1+u1r7/+urXWWpfLZX/5y1+G57N79+7w9z/+8Y/tPffcE37PSy+91B45csT++c9/tn369LEbNmyw9fX1tqCgwP773/+2u3btsmeffbY9cOCAtdbaO++80952223h99m1a1fEtZF94oknwo+9Xq9du3at/fDDD212drbdunVrk3lHeuihh+yAAQPs7t277cGDB21OTo5du3atXbdunc3NzbUHDhyw+/fvt8OHD7fr169v8fi2bdtsTk5O+LxLliyxixYtCl+njz/+OPzvKGTbtm3WGGNXrVoVPtbatQ2d7+GHH7bTpk1r9lla4vzZBOLN+5DXeh/yfv7Y67VerzfqcwCA2ElaZ1vIo6xu0U7HH3+85s6dq3vuuUc9e/YMH1+1apX+9Kc/SZLmzJmjm266KfzchRdeqC5dumj48OH64IMPwseLi4s1ePBgSdLs2bP1yiuv6JJLLtGTTz4pv9+vuro67dixQ5s3b1Z+fr4k6dJLLw2/ftOmTbrlllu0d+9eHThwQJMnTw4/N336dBljlJeXp1NPPVV5eXmSpJycHAUCAdXU1Gjz5s0aO3asJOmzzz7TmDFjon7mrKwsXXzxxc2Or169WuPHj9egQYMkSX379o36+nPPPVf9+vWTJM2aNUuvvPKKjDG66KKL1KtXr/Dxl19+WdbaqMdnzJjR5JyjRo3St771LR0+fFgXXnihPJ7oNy+5XC6VlJSEH7d2bWfPnh3+5w033BD1fAAAIDOkbLtFIl1//fV68MEH9cknn7Q4xhgT/r579+7h7xt+aWk+JvR427ZtWrJkif75z39qw4YNmjZtmmpra8NjQuFRamir+PWvf62NGzfq1ltvbTIu9J5dunRp8v5dunRRXV2drLU699xzVV1drerqam3evFkPPvhg1M/So0cPZWVlNTturW32Gdq6FqHHkdfBec5YjB8/Xi+99JJOO+00zZkzR4888kjUcZHXq61rGznPWD4XAABIX4Tko9C3b1997WtfaxIqv/KVr+jxxx+XJD322GMaN25cm+dZs2aNtm3bpiNHjuiJJ57QuHHj9PHHH6tXr1464YQT9MEHH+hvf/tbi6/fv3+/Bg4cqMOHD+uxxx5r12coKSnRq6++Gu69PXjwoP7v//5PktSnTx/t37+/zXOMGTNGlZWV2rZtmyTpo48+ijruhRde0EcffaRPP/1Uf/7znzV27FiNHz9ef/7zn3Xw4EF98sknevrpp3X22We3eNw5p2AwqFNOOUVXXnmlvv3tb2v9+vWSpG7duunw4cNR59HWtX3iiSfC/2ypqg4AADJDTO0Wxpgpkv6fpCxJy6y1dzqeN43Pny/poKR51tr1jc8FJO2XVC+pzlpbFLfZJ9D3v/99/frXvw4/vueee/Stb31LixcvVv/+/fXQQw+1eY4xY8ZowYIF2rhxY/gmvi5dumjkyJHKycnR4MGDw+0Q0fzsZz/T6NGj5XK5lJeXF1OwDenfv7+WL1+u2bNn69ChQ5KkRYsW6ayzzpLP59PUqVM1cOBArVy5stVz+P1+zZo1S0eOHNEpp5yiF154odm4cePGac6cOdqyZYvKyspUVNTwIzBv3jwVFxdLalhWbuTIka0eHzt2rHJzczV16lTl5uZq8eLF6tatm3r37h2uJPt8PuXn56ugoEB33HFHk3mMGDGi1Wt76NAhjR49WkeOHNHvf//7mK8lAABIP6atP28bY7Ik/Z+kcyXVSForaba1dnPEmPMlfUcNIXm0pP9nrR3d+FxAUpG1NuZ104qKimxo9YSQN998U8OGDQs/drtjW5EiVi6XFAjE73xtqaio0JIlS/SXv/yl894ULXK73Vq3bp1OPvnkdr/W+bMJxFvp8lJJUsW8iobHpY2PKyqaPQcAiJ0xpqqlAm4sleRiSVustVsbT/a4pJmSIneTmCnpkca7BFcbY040xgy01u44xrm3qDMDLQAAADJLLD3Jp0naHvG4pvFYrGOspOeNMVXGGN/RTjTdlJaWUkVOIoFA4KiqyAAAID3FEpKj3ebv7NFobcxYa22BpKmSrjXGjI/6Jsb4jDHrjDHrdu3aFcO0ACD9+av8qgxWJnoaAJBxYgnJNZJOj3icLen9WMdYa0P//FDS02po32jGWuu31hZZa4v69+8f2+wBIM2VbyyXJJXllSV4JgCQWWIJyWslDTHGDDLGHCfp65JWOMaskDTXNCiRtM9au8MY08sY00eSjDG9JJ0naVMc5w8Aac/r8spXSLcaAHSmNm/cs9bWGWOuk/R3NSwB91tr7RvGmKsbn18q6Vk1rGyxRQ1LwH2z8eWnSnq6cWOGrpLKrbXPxf1TAAAAAHEU02Yi1tpnrbVnWWvPtNbe0XhsaWNAVuP219c2Pp9nrV3XeHyrtXZE41dO6LWpKisrSx6PR7m5ufrqV7+qgwcPJnpKkqSf//znrT5//vnna+/evTGda+HChVqyZMkxzykQCKi8vLzVMXv37tV99913zO8FZLrKYKX8Vf5ETwMA0krq7rjndkvGxO/L7W7zLXv27Knq6mpt2rRJxx13nJYuXRrzdOvr64/+s7ahpZBsrdWRI0f07LPP6sQTT+yw94/mWENyR14vIJ2EepVDvcsAgPhI3ZAcDErWxu+rnTuTnH322eEtnS+88EIVFhYqJydHfv/n1ZzevXvrpz/9qUaPHq1Vq1bp9ttv16hRo5Sbmyufz6fQRi6lpaW64YYbNH78eA0bNkxr167VrFmzNGTIEN1yyy3h8z366KMqLi6Wx+PRVVddpfr6ei1YsECffvqpPB6PLrvsMgUCAQ0bNkzXXHONCgoKtH37drndbu3e3bCXyyOPPKL8/HyNGDFCc+bMifrZXn/9dZ1zzjkaMmSIHnjgAUkNgfvGG29Ubm6u8vLywls4t3R8wYIFevnll+XxeHTXXXfpjTfeCM89Pz9fb7/9thYsWKB33nlHHo9HN954oyoqKjRhwgSVlZUpLy+vzWv7/e9/XwUFBZo4caJYEQWZylfok9flTfQ0ACD9WGuT7quwsNA6bd68uemBcJdHnMRwvl69ellrrT18+LCdMWOGve+++6y11u7Zs8daa+3BgwdtTk6O3b17d+MpZZ944onw60PjrLX28ssvtytWrLDWWuv1eu1NN91krbX27rvvtgMHDrTvv/++ra2ttaeddprdvXu33bx5s73gggvsZ599Zq21dv78+fbhhx9uMi9rrd22bZs1xthVq1aFj7lcLrtr1y67adMme9ZZZ9ldu3Y1m0/IrbfeavPz8+3Bgwftrl27bHZ2tn3vvffsU089ZSdNmmTr6urszp077emnn27ff//9Fo+vXLnSTps2LXze6667zj766KPWWmsPHTpkDx48aLdt22ZzcnLCY1auXGm/8IUv2K1btza7ZtGubeh8t912m7322mtb+TfXsZr9bAJx5H3Ia70PeZse83qt1+ttdQwAoG2S1tkW8mgsO+6hUahiKzVUkr/97W9Lku655x49/fTTkqTt27fr7bffVr9+/ZSVlaWLL744/PqVK1fqV7/6lQ4ePKiPPvpIOTk5mj59uiRpxowZkqS8vDzl5ORo4MCBkqTBgwdr+/bteuWVV1RVVaVRo0aF53LKKadEnafL5VJJSUmz4y+++KIuueSS8KYZffv2jfr6mTNnqmfPnurZs6cmTJigNWvW6JVXXtHs2bOVlZWlU089VV6vV2vXrm3x+PHHH9/knGPGjNEdd9yhmpqacJU8muLiYg0aNCj8uKVr26VLF1166aWSpMsvv1yzZs2Kej4AAICjQUhuh1BPcqSKigr94x//0KpVq/SFL3xBpaWlqq2tlST16NFDWVlZkqTa2lpdc801WrdunU4//XQtXLgwPE6SunfvLknq0qVL+PvQ47q6Ollr9Y1vfEO/+MUv2pxnr169oh631qpxpZFWOccYY8KtIdHOGYuysjKNHj1af/3rXzV58mQtW7ZMgwcPbjYucu6tXdu25gwAAHAsUrcnOUns27dPJ510kr7whS/oP//5j1avXh11XCjcnXzyyTpw4ICeeuqpdr3PxIkT9dRTT+nDDz+UJH300UcKNvZRd+vWTYcPH47pHE8++aT27NkTPkc0zzzzjGpra7Vnzx5VVFRo1KhRGj9+vJ544gnV19dr165deumll1RcXNzi8T59+mj//v3hc27dulWDBw/Wd7/7Xc2YMUMbNmxoNsaptWt75MiR8DUsLy/XuHHj2vz8AAAAsaKSfIymTJmipUuXKj8/X0OHDo3a5iBJJ554oq688krl5eXJ7XaH2yZiNXz4cC1atEjnnXeejhw5om7duunee++Vy+WSz+dTfn6+CgoKdMcdLa+yl5OTox//+Mfyer3KysrSyJEjtXz58mbjiouLNW3aNL377rv6yU9+oi9+8Yu66KKLtGrVKo0YMULGGP3qV7/SgAEDWjzer18/de3aVSNGjNC8efNUW1urRx99VN26ddOAAQP005/+VH379tXYsWOVm5urqVOnatq0aTFf2169eumNN95QYWGhTjjhhEYxL7sAACAASURBVPANgwAAAPFgYv1zeWcqKiqy69ata3LszTff1LBhwz4/4Ha3e0WKVrlcUiAQv/OhQ/Xu3VsHDhxI9DQkRfnZBOKodHmpJKliXsXnx0obj1VUtDgGANA2Y0yVtbYo2nOpW0km0AIAAKCD0JOMlJQsVWQAAJCeCMkAAACAQ0qF5GTsn0Zm42cSAID0lDIhuUePHtqzZw+hBEnDWqs9e/aoR48eiZ4KAACIs5S5cS87O1s1NTXatWtXoqcChPXo0UPZ2dmJngYAAIizlAnJ3bp1a7JdMQAAANBRUqbdAgAAAOgshGQAAADAgZAMAAAAOBCSAQAAAAdCMgAkKX+VX5XBykRPAwAyEiEZAJJU+cZySVJZXlmCZwIAmYeQDABJzOvyylfoS/Q0ACDjEJIBAAAAB0IyAKSBymCl/FX+RE8DANIGIRkAUlyoZznUwwwAOHaEZABIcb5Cn7wub6KnAQBphZAMAAAAOBCSAQAAAAdCMgAAAOBASAYAAAAcCMkAAACAAyEZAAAAcCAkAwAAAA6EZAAAAMCBkAwAAAA4EJIBAAAAB0IyAAAA4EBIBgAAABwIyQAAAIADIRkAAABwICQDAAAADoRkAAAAwIGQDAAAADgQkgEAAAAHQjIAAADgQEgGAAAAHAjJAAAAgAMhGQCSkL/Kr8pgZaKnAQAZi5AMAEmofGO5JKksryzBMwGAzERIBoAk5XV55Sv0JXoaAJCRCMkAAACAAyEZANJEZbBS/ip/oqcBAGmBkAwAaSDUuxzqZQYAHBtCMgCkAV+hT16XN9HTAIC0QUgGAAAAHAjJAAAAgAMhGQAAAHAgJAMAAAAOhGQAAADAgZAMAAAAOBCSAQAAAAdCMgAAAOBASAYAAAAcCMkAAACAAyEZAAAAcCAkAwAAAA6EZAAAAMCBkAwAAAA4xBSSjTFTjDFvGWO2GGMWRHneGGPuaXx+gzGmwPF8ljHm38aYv8Rr4gAAAEBHaTMkG2OyJN0raaqk4ZJmG2OGO4ZNlTSk8csn6TeO578n6c1jni0AAADQCWKpJBdL2mKt3Wqt/UzS45JmOsbMlPSIbbBa0onGmIGSZIzJljRN0rI4zhsAAADoMLGE5NMkbY94XNN4LNYxd0u6SdKRo5wjAAAA0KliCckmyjEbyxhjzAWSPrTWVrX5Jsb4jDHrjDHrdu3aFcO0AAAAgI4RS0iukXR6xONsSe/HOGaspBnGmIAa2jTOMcY8Gu1NrLV+a22Rtbaof//+MU4fAAAAiL9YQvJaSUOMMYOMMcdJ+rqkFY4xKyTNbVzlokTSPmvtDmvtj6y12dZad+PrXrTWXh7PDwAAAADEW9e2Blhr64wx10n6u6QsSb+11r5hjLm68fmlkp6VdL6kLZIOSvpmx00ZAAAA6FhthmRJstY+q4YgHHlsacT3VtK1bZyjQlJFu2cIAAAAdDJ23AMAAAAcCMkAkGT8VX5VBisTPQ0AyGiEZABIMuUbyyVJZXllCZ4JAGQuQjIAJCGvyytfoS/R0wCAjEVIBgAAABwIyQCQRiqDlfJX+RM9DQBIeYRkAEgToR7mUE8zAODoEZIBIE34Cn3yuryJngYApAVCMgAAAOBASAYAAAAcCMkAAACAAyEZAAAAcCAkAwAAAA6EZAAAAMCBkAwAAAA4EJIBAAAAB0IyAAAA4EBIBgAAABwIyQAAAIADIRkAAABwICQDAAAADoRkAAAAwIGQDAAAADgQkgEAAAAHQjIAAADgQEgGAAAAHAjJAAAAgAMhGQAAAHAgJAMAAAAOhGQAAADAgZAMAAAAOBCSAQAAAAdCMgAAAOBASAYAAAAcCMkAAACAAyEZAJKEv8qv0uWlqt5ZneipAEDGIyQDQJIo31iu6p3V8gzwqCyvLNHTAYCM1jXREwAAfM4zwKOKeRWJngYAZDwqyQCQZqp3Vqt0ean8Vf5ETwUAUhaVZABII6E2jVBfs6/Ql8jpAEDKopIMAGnEV+hTxbwKeQZ4Ej0VAEhphGQAAADAgZAMAAAAOBCSAQAAAAdCMgAAAOBASAYAAAAcCMkAAACAAyEZAAAAcCAkAwAAAA6EZAAAAMCBkAwAAAA4EJIBAAAAB0IyAAAA4EBIBgAAABwIyQAAAIADIRkAAABwICQDAAAADoRkAAAAwIGQDAAAADgQkgEAAAAHQjIAAADgQEgGAAAAHAjJAAAAgAMhGQAAAHAgJAMAAAAOhGQAAADAgZAMAAAAOBCSAQAAAIeYQrIxZoox5i1jzBZjzIIozxtjzD2Nz28wxhQ0Hu9hjFljjHndGPOGMea2eH8AAAAAIN7aDMnGmCxJ90qaKmm4pNnGmOGOYVMlDWn88kn6TePxQ5LOsdaOkOSRNMUYUxKnuQNAc253w1fkY2OaHgMAoA2xVJKLJW2x1m611n4m6XFJMx1jZkp6xDZYLelEY8zAxscHGsd0a/yy8Zo8ADQTDDZ8RT62jf/bSeKg7K/yqzJYmehpAAAaxRKST5O0PeJxTeOxmMYYY7KMMdWSPpT0grX2X0c/XQCIEFk1DlWMXa7ojwOBz48nofKN5ZKksryyuJ2zMlgpf5U/bucDgEwSS0g2UY45q8EtjrHW1ltrPZKyJRUbY3KjvokxPmPMOmPMul27dsUwLQAZL7JqHKoYBwKfB+PQ45BAoGmVOcl4XV75Cn1xOVcobIfCNwCgfWIJyTWSTo94nC3p/faOsdbulVQhaUq0N7HW+q21Rdbaov79+8cwLQBQQyB2VocDgabh2DnemKZ9ys4+5jTgK/TJ6/ImehoAkLJiCclrJQ0xxgwyxhwn6euSVjjGrJA0t3GVixJJ+6y1O4wx/Y0xJ0qSMaanpEmS/hPH+QPIdO2tDgcCDRVmaxteFwrHSVxhBgB0vq5tDbDW1hljrpP0d0lZkn5rrX3DGHN14/NLJT0r6XxJWyQdlPTNxpcPlPRw4woZXSQ9aa39S/w/BgAchch+ZROtawwAkKnaDMmSZK19Vg1BOPLY0ojvraRro7xug6SRxzhHAGhdKOy2V0stGQCAjMeOewBST7SVK+IReFlTGQDQKKZKMgAklci1j+MlckUMWi8AIOMRkgFAovUCANAE7RYAUovbffQ9yAAAxIhKMoDUELlUW7xbLZxCay9TXQaAjEVIBpAa2ljH2O3+fIjLdYz5liXhACDjEZIBpI5W2iwiC8yR+TYUno85OAMAMgohGUDqiJJyW1qtLTIchxasCI2NKSzTcgEAGY0b9wAkv1Zu1otssQiJXM0tMuMGg+3YfTr0QtZMBoCMRCUZQPJz3KwXWSWWmhd7oxV/I0N0zAViepMBIGMRkgGknFBmbk9+jQzF5F4AQFtotwCQvJzbT3fQW7TaUeFysVU1AGQgKskAklcHrYkcyr2h71vtUw6VoCk/A0BGISQDSFlHW2B29iOH8m+by8W1a3kMAEAqIyQDSCmRC13EM6uG8m+Lvc6hN415eQwAQCojJANIKR3RgRFT6KbtAgAyCiEZQEroyPvmnOGYfUQAAIRkAMnJsYFIZ3Y5sDwyAIAl4AAkp2AwXMptZcO9zhcqMwMA0hqVZADJJ0oVub2bh3TElCQpQJkZADICIRlA8mlMxc7tpztb5IpvLGoBAJmFkAwgaTlXsujMsNzqim/c2QcAaY+QDCC5tNKA3JmZtNUV32i5AIC0x417AJJLxA17ySC0hXXS3DgIAOgUVJIBoBWReZ3iMQBkDirJAJJSUi37Fk2oxMxycACQlqgkA0hKHbH9dFyxTTUApDUqyQASz+1uCJuNzb9JX0WOM3+VX5XByg45d2WwUv4qf4ecGwDSGSEZQOKFysbWSoFAst271+HKN5ZLksryyuJ63tD5QucHAMSOkAwAMerIHam9Lq98hb64ntNX6JPX5Y3rOQEgUxCSASBGkW3I4bDckckZAJAwhGQAaIdA4PMbCo2R3AqwZzUApCFCMoCkkio37YXCcjAoqskAkIYIyQASqzEVhxa4kFLwpr1AgGoyAKQZQjKAxAoGG1oWFF7cIjVRTQaAtMJmIgASLuk3DolFIMDGIgCQRqgkAwAAAA6EZAAAAMCBkAwgcVJlKYtW0IoMAOmJkAwgcRpv2kvlnNxkYQsSMwCkDW7cA5AYjVXktLhpL4Sb9wAgbVBJBpAYwWAKr/fWFAVkAEg/VJIB4BhRQAaA9EMlGQDihGoyAKQPQjKAhEmDxS3CQp8jIFdDWZnEDAApjXYLAAmVJm3J4c9hTKDhRkT6LwAgpVFJBpAw6RKQAQDph5AMAB2BJS8AIKURkgF0LrdbMqahdzedNdllBACQagjJADpX4+4hgxRI9Ew6hIv79gAgLRCSAXQet1s1WS4Zkz6rWjgFAmm0gyAAZDBWtwDQeYJBnS6bESEyEJACxiV36DcC7lIEgJRCSAaADjJILAcHAKmKdgsAnSpd2ywAAOmFkAygU9F1AABIBYRkAAAAwIGQDAAdhP1EACB1EZIBdIpMDIvJsp9IZbBS/ip/oqcBACmFkAygU1QE3Zl7114CdxgpyyuTJJVvLO/09waAVEZIBtAp3Apm7l17oR1GopSV/VV+VQYrO+ytfYU+eV3eDjs/AKQrQjKAjuN2N1RQjVFNVmZWkdvqSw5VeEMVXwBAcmAzEQAdJxiUrJUxkq1L9GQSIxCI2EsklJgdFXWvyytfoa+TZwYAaA0hGQA6S5PEDABIZrRbAAAAAA6EZADoTCyeDAApgXYLAOhMtFwAQEqgkgygQxmTucsjO4UW+wAAJD8qyQA6htutgFyyNtETSbzQXiIuV8NyySIoA0DSi6mSbIyZYox5yxizxRizIMrzxhhzT+PzG4wxBY3HTzfGrDTGvGmMecMY8714fwAASSoY1CAFEj2LpBDaSyS08ltA9CUDQLJrs5JsjMmSdK+kcyXVSFprjFlhrd0cMWyqpCGNX6Ml/abxn3WSvm+tXW+M6SOpyhjzguO1ANIUbRbRDVJANkg5GQCSWSyV5GJJW6y1W621n0l6XNJMx5iZkh6xDVZLOtEYM9Bau8Nau16SrLX7Jb0p6bQ4zh9AEsvUXagBAKkvlpB8mqTtEY9r1DzotjnGGOOWNFLSv9o7SQAAAKAzxRKSo/1N0HkrTqtjjDG9Jf1R0vXW2o+jvokxPmPMOmPMul27dsUwLQAAAKBjxBKSaySdHvE4W9L7sY4xxnRTQ0B+zFr7p5bexFrrt9YWWWuL+vfvH8vcASB1sakIACS1WELyWklDjDGDjDHHSfq6pBWOMSskzW1c5aJE0j5r7Q5jjJH0oKQ3rbX/E9eZA0haNV3dqsnirr1WBQJSMJjoWQAAWtDm6hbW2jpjzHWS/i4pS9JvrbVvGGOubnx+qaRnJZ0vaYukg5K+2fjysZLmSNpojKluPHaztfbZ+H4MAMkkuz4oFkiOgcslVVZK3btL8xI9GQBApJg2E2kMtc86ji2N+N5KujbK614Ry+YDmaVxExF3oueRCgIBaXlpQ1AGACQVtqUGEF/BoEpdgUTPIqnRjgwAyY+QDCDuWB+5daHrQ1AGgOQVU7sFACC+AgHJGNGWAgBJikoygLihMgoASBeEZABxw4pm7RNa3AIAkHwIyQCQIIGA5PUmehYAgGgIyQCQYLXq3tig7E70VAAAjQjJAJBg/1JJw+Yr9KsAQNIgJANAgnXvThEZAJINIRnAMavp6paMkZVpuBsN7VJS0vFF5Oqd1SpdXip/lb9j3wgA0gTrJAM4Ztn1wYZ2ARw1l0tSBwXlsrwySQ1BWZJ8hb6OeSMASCNUkgEgCXTkLoW+Qp8q5lXIM8DTcW8CAGmGkAzg2LjdCogWi/byV/lVurw0XN2VpJosV0PrCgAg4QjJAI6a2y0pGNQgBRI8k9RTvrFc1Tur5RngCbdDZNcFVFcvloMDgCRATzKAoxa62Yx79Y6OZ4BHFfMqmhwrdQUUDEo2aBIzKQCAJEIygDjoyH7aTBO+lmRkAEgo2i0AAAAAB0IyAAAA4EBIBoAkFJCLm/cAIIHoSQZw1LbJzV17HWSQAty8BwAJREgGcNTcCkoBdtoDAKQf2i0AIAlRoAeAxCIkAzgqtMt2LJbVA4DEIiQDOCqhjUTQcbh5DwASh5AMoN3cbtoBOsMgBfhtBAAShJAMoN2CQSnAyhYAgDTG6hYAjk4wKFlWtuhoNVkuqatb2XWBRE8FADIKIRkAkpTLJY1TQAHWSwaATkdIBoAkFV7hgowMAJ2OnmQAMXO7JWNoRQYApD9CMoCYhdqQuWmv87ESHAB0LtotALQfN+11LpdLFUG3pEBi5wEAGYRKMgAku0BAbgVlDBVlAOgshGQASBHWsrcIAHQWQjIAZJDKYKX8Vf5ETwMAkh4hGUBM3G5pe5ab5S0SxeWSjNE2uY/6FGV5ZZKk8o3lcZoUAKQvQjKAmASDUnZ9aHmLQKKnk3kCAclauRXU6tVHdwpfoU9elzeu0wKAdEVIBtAmt5vicTI5dCjRMwCA9EdIBtCmYJDicTz5q/yqDFYmehoAgFYQkgHEhnJy3IR6gkM9wu3icmm0jrLfAgAQMzYTARAbNhCJK6/LK1+hr/0vDATUw5ij7ksGAMSGSjIApCD6kgGgYxGSAQAAAAdCMoA2bZObfuRk0r07fckA0MEIyQDa5BbLWySVkhL1EP0WANCRCMkAWuV2J3oGAAB0PkIygFYFg4meAVridjfsEs4vMgAQf4RkAEhBtequiqBb1vKLDAB0BEIygKhCVcrtWW5u2ktC1d1L1DVLkjENN1YCAOKKzUQARBXeO8QEpQCbiCSbkhIpuyIgSXIbk9jJAEAaopIMAGmAvmQAiC9CMoBm3O7GNgtjaLVIEfQlA0B80W4BoJlgUMpWqN8CAIDMQyUZAAAAcCAkAwAAAA6EZABIdS6XrIwCxs0NfAAQJ4RkAE243dyrl3ICAclauRWM6Qa+ymCl/FX+Dp8WAKQyQjKAJoLBhsyF9FSWVyZJKt9YnuCZAEByIyQDCKOKnOJcLm3Par3lwlfok9fl7bQpAUCqYgk4AGHhXfZIy6kpEFC2MayZDABxQCUZQHP0XKQ0l6thHxhu4gOAo0clGQDSTOj3G2MSOg0ASGlUkgGgE/mr/KoMVnbcG7hclJABIA6oJAOQFNGGTD9yhwqtKhFaZSLuAgFKyAAQB4RkAJIibtozoW/QUbwur3yFvo57g8am5G1ySQp03PsAQBqj3QLIcG53Q+GR4nEaidhcBABwdAjJQIYLVZADAdFqkYZoTwaAo0NIBvA5ln5LLy6XXqlpfXMRAEB0MYVkY8wUY8xbxpgtxpgFUZ43xph7Gp/fYIwpiHjut8aYD40xm+I5cQBAGwIBZdcH2VwEAI5CmyHZGJMl6V5JUyUNlzTbGDPcMWyqpCGNXz5Jv4l4brmkKfGYLIAORKsFAABhsVSSiyVtsdZutdZ+JulxSTMdY2ZKesQ2WC3pRGPMQEmy1r4k6aN4ThpAB6DVIj25XNomd6JnAQApJ5aQfJqk7RGPaxqPtXcMgCRD8TgDBAJyK0hfMgC0UyzrJEdbld65iGosY1p/E2N8amjV0BlnnNGelwI4SkGWRM4Y9CUDQPvEUkmukXR6xONsSe8fxZhWWWv91toia21R//792/NSAEAbGvcXoaIMADGKJSSvlTTEGDPIGHOcpK9LWuEYs0LS3MZVLkok7bPW7ojzXAHEUZNtqNlNJO017i9CRRkAYtRmu4W1ts4Yc52kv0vKkvRba+0bxpirG59fKulZSedL2iLpoKRvhl5vjPm9pFJJJxtjaiTdaq19MN4fBED7sA11BnG5Gn4Z4sZMAIhZLD3JstY+q4YgHHlsacT3VtK1Lbx29rFMEABwjAKBhr8WRKgMVspf5Zev0JeYOQFAkmPHPQDIIC6XtLG8TJJUvrE8wbMBgORFSAYyUJN+ZHqRM0ogIH30gk9elzfRUwGApBZTuwWA9EI/cgYKLW8haXuWS19a7VZJSYLnBABJjJAMAJkg4qa9bGN06JA7YVMBgFRAuwUAdBJ/lV+VwcpET0OS1L27tHp1omcBAMmLkAxkkCZLItOP3OlCN8qV5ZUldiIul0q0WocOJXYaAJDMaLcAMkiTbajpR04Ir8ub+GXXAgHpm6bNYQCQyagkA0CG6t6dbaoBoCWEZCBD0F0Bp5IStqkGgJbQbgFkiCatFiRmAABaRUgGMlGQfuSMx/IWANAq2i0AIBM17iSyPcstY+hNBgAnQjKQAZpsQx1eAw4Zr6RE2fXB8B8VCMoA8DnaLYAMwDbUaEsg0BCS3e4mm/MBQMaikgykOe7RQ6tcLoX6LULhmIoyABCSgbQXDFIZRHSVwUr5/3izIvstAgGWhQMAiZAMZA5KyogQ2ho7tFV2ZDqOKC4DQMaiJxnIFCz7hgi+Qt/nAdkh9JcHw87VADIYlWQAAADAgZAMZAJaLQAAaBdCMpCGQsshh5dE5u69hPNX+VUZrEz0NFrnaEZ2uehLBpC56EkG0lCz9mN6SxMu1P8bumEuKTmakQMB+pIBZC4qyUC6o9UiaXhdXvkKfYmeRruw0gWATEUlGUgzzTIxq1qgvULJ2OVSoLG6TEUZQKYhJANphkyMY8YacABAuwUAAADgREgG0gjtx4grlrcAkMEIyUAaabbSG6kZxyL0w+R2k5cBZBxCMpDOWB8ZbajeWa3S5aXyV/mjDwgEpGCwSZsyYRlAJiAkA+kotJsIVWS0oiyvTJ4BHlXvrA6v49yaQODzm0IJygDSHSEZSBNNOitCS1xQRUYrfIU+VcyrkGeAp/WBjsWSG4vLAJDWWAIOSBMs/YYOw5JwADIQlWQgxTXrrOBmPXQCbuQDkO4IyUCKa9ZZwc16ScVf5Vfp8lJV76xO9FSOXUTbRcTCFwCQlgjJQDqhipx0yjeWq3pntTwDPCrLK0v0dI6N4849epMBpDN6koF0QmNyUvIM8KhiXkWipxE/gQD9yQDSHpVkIIVROEai0ZsMIF0RkoEUFQom4fZjEjM6U2M6ZpMRAOmKkAykqGb353HDHo5BZbCy5V33oom4cy/Uqkx/MoB0QkgGUh276+EYhW4ojGXXvSYcd+459hwBgJRGSAZSELvrIZ58hT55Xd6jP0HjL2oBuakoA0gbrG4BpJgmvcj0ISPRQj9/1rLiBYC0QkgGUkyTVd5Y8g2JFvkXjMab+VyugNxu/rgBILXRbgGkCLafRtJr7FFmxQsA6YCQDKQItp9OPf4qvyqDlYmeRkI4NucDgJRDSAaADhJaLSLlt6NuD8fuImxdDSBVEZKBVESrRcrwurzyFfoSPY3OEyUVszQcgFRESAZSQJNM3GyrPSA+2r2hSEuiVJNpvQCQagjJQApo0n5MLzI6wFFvKBJNC3fu0XoBIJUQkoEkR2cFOsMxbyji1EL5mNYLAKmCkAwksWYbh7D9NFJNqKrc+MNM6wWAVEFIBpJYuLMilCbYfjplZPLyb81E6bNwZGcASDqEZCAJNSkac6NeSsrI5d9aE6XPgqAMIJkRkoEk06xozI16KSsVl3+L2woXTqE+ixYqysbQqwwguRCSgSTSrGjMXXvoRHFd4aIloYpyRCoO5ecoGRoAEoaQDCSRZkVjqsjoRHFf4SKayETM6hcAkhghGUgS4aJxqCGZlSyQCRw39bH6BYBkQUgGkkCTNotg8PMqG1XklMTKFu3Uyk199CoDSBRCMpBAoaKxRB5OJ6m+skX1zmqVLi/tmBv4ooksH0ck4sjODImgDKBzEZKBBAoVjdksJP2k4soWUkOw9wzwqHpndcfewBdNC2E59FToMJVlAJ2ha6InAGSqZgtXhBIzkEC+Qp98hT6VLi9N3CQiE3GUw6GnWEIcQEeikgx0sqgtFiz1BrTM7W5WOg7958KScQA6CiEZ6GThFgu5P//bsUQ5LA34q/wqXV6q6p3ViZ5KXHTYxiKxCt3QFxIKy8YooIb1laMsuwwAcUFIBjpB1FXdWMUi7ZRvLFf1zmp5BnhS9qa9kE7ZWKQtoR7lQKDpfyMRu45EW3aZsAwgHuhJBjpY5DbT4QMmSHtFmvIM8KhiXkWip3HMfIW+xAbkaCKDcmSV2eUKPxcaEio68/sngKNFJRnoIM16jyMPUD1OO+m6NnLCWy5a4lwfztFv4VwNg+oygPYiJANxFNlWITmWd2tyAOkm1ddGjiYpWi5i0UJgDsgdtRWD0AwgFoRkII6atRmzW0hGCFWRU3Vt5Jb4Cn3yurzJW02OpqUKcxuhOeKeQMIzAEkxhmRjzBRjzFvGmC3GmAVRnjfGmHsan99gjCmI9bVAqmt2U16L5WSkI3+VX1f95SpJ6VVFDgl9pqv+clXqBOWQlu7qixKaQ6JVnAnOQGZqMyQbY7Ik3StpqqThkmYbY4Y7hk2VNKTxyyfpN+14LZAyIvNvkxzscsvKKBCMCMaE47QXGZDvv+D+tKoih/gKfbr/gvslNQTlTt2uOt5aCc2BYON/v1Eqzi0F58gKNCEaSD+xrG5RLGmLtXarJBljHpc0U9LmiDEzJT1irbWSVhtjTjTGDJTkjuG1QFJyuz/fqGCb3HIrqEC0gUE1lJDZLS/t+av8TfpzQzfqpWtADgl9tvKN5aoMVqoyWNnkOpTllaXm52/pl9jINqnQ0BZOUVPj0rjshmcdL2kmYhEOACkglpB8mqTtEY9rJI2OYcxpMb42KVz/3PVpswEAGqxe3fDPkpKmxw4davh+tFarhw61+Hr3hIbf8iRpXvfuUom39TdM5Da+6BShUOx1ecP/TNmA2E6h7aqj/aLgDM0pb6Fbn//X34bVq+U+1JCOY3lF6TdjO22tukuSeuiQatVd/1KJRmt1kzHV3UtUUtL0/2sh3RterkOHGr6P/P8gkIw8Azy6e8rdiZ5GE8a2Uf0yxnxV0mRr7RWNj+dIKrbWfidizF8l/cJa+0rj439KuknS4LZeG3EOnxpaNSRpXV99jQAABcVJREFUqKS3jvGzHY2TJe1OwPumKq5X+3C92ofr1T5cr/bherUf16x9uF7tk6jr5bLW9o/2RCyV5BpJp0c8zpb0foxjjovhtZIka61fUkIb3Ywx66y1RYmcQyrherUP16t9uF7tw/VqH65X+3HN2ofr1T7JeL1iWd1iraQhxphBxpjjJH1d0grHmBWS5jauclEiaZ+1dkeMrwUAAACSSpuVZGttnTHmOkl////t3U+oVGUYx/HvjwwhKzIqKSq0yDYtNKKNKEJlKaEVFEqEWYsEg6SNWYvElVm2aJMQCi7UNEqKsNAoamV/NEFNw6xbaHKlXJgYgfa0OO/kOPe8xztQ9xw9vw9cZu4zM5f3Pvc5z7x3zp8XuARYGxH7JC1Mj68GtgKzgB+AU8CCqtf+L7+JmZmZmdl/ZDiHWxARWykmwt2x1V33A1g03Nc22AV6XaPaOF/9cb7643z1x/nqj/PVP+esP85XfxqXr/OeuGdmZmZm1jZeltrMzMzMrEdrJ8mSBiTtkbRb0jcpdrWk7ZIOptuxdY+zKSTdnnLV+TohabGkZZKOdMVn1T3WukhaK+mYpL1dsWxNSVqalmv/XtL99Yy6Ppl8vSrpQFrefoukq1J8vKQ/u+psdf4nX5wy+cpuf66v0nxt6srVgKTdKe76km6S9Jmk/ZL2SXouxd3DSlTkyz2sREW+Gt3DWnu4haQB4K6I+K0rthI4HhErJL0AjI2IJXWNsalULDd+hGJhmAXAyYh4rd5R1U/SNOAkxeqTd6RYaU2pWJ59I8WKljcAnwATI+JMTcMfcZl8zQA+TSf9vgKQ8jUe+LDzvDbK5GsZJduf66s8Xz2Pr6K4EtNy1xeoWCX3+ojYJekKYCfwEPAk7mFDVOTrRtzDhqjI12M0uIe19pPkjDnAunR/HcUf0Ia6BzgUET/XPZAmiYgvgOM94VxNzQHejoi/IuIniivD3D0iA22IsnxFxLaIOJ2+3UHxhmNk6yvH9VWRL0mieHPeOKKDarCIOBoRu9L9P4D9FKvmuoeVyOXLPaxcRX3lNKK+2jxJDmCbpJ0qVvsDGJeu70y6va620TXbXM59c3k27VpaKx+i0itXU7ml3O2sp4CPur6fIOlbSZ9LmlrXoBqobPtzfVWbCgxGxMGumOsrSZ96Tga+xD3svHry1c09rERJvhrbw9o8SZ4SEXcCM4FFadecnYeKRWFmA++k0JvArcAk4CiwqqahXWhUEmvnsU8lJL0EnAbWp9BR4OaImAw8D2yQdGVd42uQ3Pbn+qo2j3P/0Xd9JZIuB94FFkfEiaqnlsRaV2O5fLmHlSvJV6N7WGsnyRHxa7o9Bmyh+Bh/MB030zl+5lh9I2ysmcCuiBgEiIjBiDgTEX8Db9Gi3W3DlKup4Sz33kqS5gMPAo+na7CTdrn9nu7vBA4BE+sbZTNUbH+urwxJo4BHgE2dmOurIOlSignM+oh4L4XdwzIy+XIPyyjLV9N7WCsnyZLGpAPHkTQGmAHspVgye3562nzg/XpG2GjnfALTaZ7JwxR5tLNyNfUBMFfSaEkTgNuAr2oYX6NIegBYAsyOiFNd8WvTCaNIuoUiXz/WM8rmqNj+XF959wIHIuJwJ+D6+vc47TXA/oh4vesh97ASuXy5h5WryFeje9iwVty7CI0DthR/M0YBGyLiY0lfA5slPQ38Ajxa4xgbR9JlwH3AM13hlZImUewGGeh5rFUkbQSmA9dIOgy8DKygpKbS0u6bge8odsktastZ4R2ZfC0FRgPb0/a5IyIWAtOA5ZJOA2eAhREx3JPYLgqZfE0v2/5cX+X5iog1DD2nAlxfAFOAJ4A9SpfGA17EPSwnl683cA8rk8vXvCb3sNZeAs7MzMzMLKeVh1uYmZmZmVXxJNnMzMzMrIcnyWZmZmZmPTxJNjMzMzPr4UmymZmZmVkPT5LNzMzMzHp4kmxmZmZm1sOTZDMzMzOzHv8AYhd0wqME9zAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pyplot.figure(figsize=(12, 8))\n",
    "pyplot.hist(t_boot_nonparam, bins, label='Nonparametric bootstrap', color='blue', histtype='step', density=True)\n",
    "pyplot.hist(t_boot_param, bins, label='Parametric bootstrap', color='red', histtype='step', density=True)\n",
    "pyplot.step(bins, theta_cdf_bins_delta, color='green', label='True sampling distribution')\n",
    "pyplot.axvline(x=np.exp(5), color='black', label='θ (true parameter)')\n",
    "pyplot.legend(loc='upper left')\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 9.6.7**. Let $X_1, \\dots, X_n \\sim \\text{Uniform}(0, \\theta)$.  The mle is $\\hat{\\theta} = X_\\text{max} = \\max \\{ X_1, \\dots, X_n \\}$.  Generate a dataset of size 50 with $\\theta = 1$.\n",
    "\n",
    "**(a)**  Find the distribution of $\\hat{\\theta}$.  Compare the true distribution of $\\hat{\\theta}$ to the histograms from the parametric and nonparametric bootstraps.\n",
    "\n",
    "**(b)**  This is a case where the nonparametric bootstrap does very poorly.  In fact, we can prove that this is the case.  Show that, for parametric bootstrap $\\mathbb{P}(\\hat{\\theta}^* = \\hat{\\theta}) = 0$ but for the nonparametric $\\mathbb{P}(\\hat{\\theta}^* = \\hat{\\theta}) \\approx 0.632$.\n",
    "\n",
    "Hint: show that $\\mathbb{P}(\\hat{\\theta}^* = \\hat{\\theta}) = 1 - (1 - (1/n))^n$ then take the limit as $n$ gets large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "X = np.random.uniform(low=0, high=1, size=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2778c3c9192340f19f258cb3b3d22f51",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "95% confidence interval (Normal): \t 0.959, 1.037\n"
     ]
    }
   ],
   "source": [
    "# Nonparametric bootstrap\n",
    "\n",
    "theta_hat = X.max()\n",
    "\n",
    "B = 1000000\n",
    "t_boot_nonparam = np.empty(B)\n",
    "n = len(X)\n",
    "for i in tqdm_notebook(range(B)):\n",
    "    xx = np.random.choice(X, n, replace=True)\n",
    "    t_boot_nonparam[i] = xx.max()\n",
    "    \n",
    "se_boot = t_boot_nonparam.std()\n",
    "\n",
    "alpha = 0.05\n",
    "z = norm.ppf(1 - alpha/2)\n",
    "normal_conf = (theta_hat - z * se_boot, theta_hat + z * se_boot)\n",
    "\n",
    "print('95%% confidence interval (Normal): \\t %.3f, %.3f' % normal_conf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2c62801bcf6d43cea0a5b7aafdff4db0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "95% confidence interval (Parametric Normal): \t 0.960, 1.036\n"
     ]
    }
   ],
   "source": [
    "# Run parametric bootstrap\n",
    "\n",
    "theta_hat = X.max()\n",
    "\n",
    "B = 1000000\n",
    "t_boot_param = np.empty(B)\n",
    "n = len(X)\n",
    "for i in tqdm_notebook(range(B)):\n",
    "    xx = np.random.uniform(low=0, high=theta_hat, size=50)\n",
    "    t_boot_param[i] = xx.max()\n",
    "    \n",
    "se_boot_param = t_boot_param.std()\n",
    "\n",
    "alpha = 0.05\n",
    "z = norm.ppf(1 - alpha/2)\n",
    "normal_conf_param = (theta_hat - z * se_boot_param, theta_hat + z * se_boot_param)\n",
    "\n",
    "print('95%% confidence interval (Parametric Normal): \\t %.3f, %.3f' % normal_conf_param)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the true sampling distribution,\n",
    "\n",
    "$$\\hat{\\theta} = \\max \\{ X_1, \\dots X_n \\}$$\n",
    "\n",
    "Its CDF is\n",
    "\n",
    "$$\\mathbb{P}(\\hat{\\theta} \\leq x) = \\prod_{i=1}^n \\mathbb{P}(X_i \\leq x) = F_{\\text{Uniform}(0, \\theta)}(x)^n$$\n",
    "\n",
    "where\n",
    "\n",
    "$$F_{\\text{Uniform}(0, \\theta)}(x) = \\begin{cases}\n",
    "0 & \\text{if } x \\leq 0 \\\\\n",
    "\\frac{x}{\\theta} & \\text{if } 0 < x \\leq \\theta \\\\\n",
    "1 & \\text{if } \\theta < x\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.linspace(0.75, 1.05, 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate the CDF for theta, calculate it for each bin, and include the differences between bins\n",
    "\n",
    "def theta_cdf(x):\n",
    "    if x <= 0:\n",
    "        return 0\n",
    "    if x >= 1:\n",
    "        return 1\n",
    "    return x**50\n",
    "\n",
    "theta_cdf_bins = list(map(theta_cdf, bins))\n",
    "theta_cdf_bins_delta = np.empty(len(bins))\n",
    "theta_cdf_bins_delta[0] = 0\n",
    "theta_cdf_bins_delta[1:] = np.diff(theta_cdf_bins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsYAAAHSCAYAAADvxw2lAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df5xWdZ3//8eLHwsE+AMEI1EGXHJh+DHggBDKD7HASFDLUot0c6XUbv24ubbYL7WVsmTLj7dyXczUvlpKbqbf1vroGmS2FKKLCKgrCgpCMuCSEEH8eH3+mItpgBnmgpmLGeBxv92u25zzvt7nnPd13nP0yXve1zmRmUiSJElHulbN3QBJkiSpJTAYS5IkSRiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEkAtGnuBgAcd9xxWVZW1tzNkCRJh5mXXnoJgFNOOaWZW6KW4plnnlmXmd3qeq9FBOOysjIWLFjQ3M2QJEmHmbFjxwIwd+7cZm2HWo6IeK2+95xKIUmSJGEwliRJkgCDsSRJkgS0kDnGddm2bRurVq1iy5Ytzd0UHaHat29Pz549adu2bXM3RZIkHQQtNhivWrWKzp07U1ZWRkQ0d3N0hMlM1q9fz6pVq+jdu3dzN0eSJB0ELXYqxZYtW+jatauhWM0iIujatat/sZAk6QjSYoMxYChWs/L3T5KkI0uLDsbNbfbs2fTr149+/frxk5/8pM46mcmZZ57J22+/zYYNG7jtttsOciubX2M+91lnncX//u//NnGLJEmS9p/BuB4bN27km9/8Jk8//TTz58/nG9/4Bps2bdqr3qOPPsrgwYM56qij9hkQd+zYUeom71NmsnPnzpLs+0CC8a72TJ069Yj8x4QkSWp5DMb1+OUvf8nZZ59Np06d6Ny5M2eddRa/+MUv9qp33333MWXKFACmT5/OK6+8QkVFBddccw1z585l3LhxXHzxxQwcOJAVK1YwYMCAmm1nzpzJ9ddfD8Arr7zCxIkTOfXUUznjjDN48cUX9zrW9ddfz9SpUznzzDPp27cvd9xxBwCbNm1i/PjxDB06lIEDB/Lwww8DsGLFCvr168eVV17J0KFDWblyJVdccQWVlZWUl5dz3XXX1ey7rKyML37xi4wcOZLKykqeffZZJkyYwMknn8ztt99eU+/mm29m2LBhDBo0qGb7PT93ffXqas/kyZP58Y9/fMD9JEmS1FRa7F0pavvc5z7HwoULm3SfFRUV3HLLLfW+//rrr3PXXXfx85//HICqqiqOP/74ver99re/5d/+7d8AuOmmm1i8eHFNW+fOncv8+fNZvHgxvXv3ZsWKFfUeb9q0adx+++307duX3//+91x55ZX86le/2qveokWL+N3vfsef/vQnhgwZwqRJk+jevTsPPfQQRx11FOvWrWPEiBFMnjwZqH5G/F133VUzKjtjxgy6dOnCjh07GD9+PIsWLWLQoEEAnHjiicybN4/Pf/7zXHrppfz2t79ly5YtlJeX86lPfYrHHnuMl19+mfnz55OZTJ48mSeffHKvz11fvZNOOmmv9gBs3bqV9evX07Vr13rPjyRJUqkdEsG4OWQmV1xxBV/+8pcB+NrXvlbnl7HeeustOnfuXO9+hg8f3uDtvjZt2sR//dd/ccEFF9SUbd26tc66U6ZMoUOHDnTo0IFx48Yxf/58Jk2axBe/+EWefPJJWrVqxRtvvMGbb74JQK9evRgxYkTN9rNnz2bWrFls376dNWvWsHTp0ppgvCtMDxw4kE2bNtG5c2c6d+5M+/bt2bBhA4899hiPPfYYQ4YMqWn3yy+/zEknnbRbG/dVb8/2AHTv3p3Vq1cbjCVJUrM6JILxvkZ2S6VXr167jdiuXLmSCRMm7FWvTZs27Ny5k1at6p6V0rFjx73q7rLrVmA7d+7kmGOOKWpUfM9wHhHcd999VFVV8cwzz9C2bVvKyspq9l37+MuXL2fmzJk8/fTTHHvssVx66aW73Y6sXbt2ALRq1apmedf69u3byUyuvfZaPvnJT+7Whj1HwvdVr3Z7ap+HDh06NPjZJUmSSsk5xvU4++yzmTNnDuvXr2fdunXMnTuXiRMn7lXvlFNO4dVXXwWgc+fObNy4sd59Hn/88axdu5b169ezdevWmmkaRx11FL17966580Vm8txzz9W5j4cffpgtW7awfv165s6dy7Bhw/jjH/9I9+7dadu2LXPmzOG1116rc9u3336bjh07cvTRR/Pmm2/WOWd6XyZMmMAPfvCDmi8hvvHGG6xdu3avz11fvbpkJn/4wx8oKyvbr7ZIkiQ1tUNixLg5dOrUiRtvvJGRI0eSmXzjG9+gU6dOe9WbNGkSc+fO5W//9m/p2rUro0aNYsCAAZx99tlMmjRpt7pt27blq1/9Kqeddhq9e/fm7/7u72reu++++7jiiiu48cYb2bZtGxdeeCGDBw/e63jDhw9n0qRJvP7663zlK1/hXe96Fx/96Ec555xzqKyspKKiYrf91jZ48GCGDBlCeXk5ffr0YdSoUft1Tt73vvfxwgsvMHLkyJpzdO+993LyySfv9rlvvvnmOuu1bt16r30+88wzjBgxgjZt/FWUJEnNKzKzudtAZWVlLliwYLeyF154gX79+jVTi4q3Zs0aPv7xj/P444+X/FjXX389nTp14h//8R9LfqyD5bOf/SyTJ09m/Pjxzd2UOh0qv4eSpLqNHTsWqP5CvAQQEc9kZmVd7zmVopF69OjB5Zdfzttvv93cTTkkDRgwoMWGYkmSdGTx79dN4MMf/vBBOc6uex4fTi6//PLmboIkSRLgiLEkSZIEGIwlSdJhoqwMIqpf3uxIB8KpFJIk6bDw2muw654CdTyTS2qQI8aSJEkSBuN9igiuvvrqmvWZM2cell+Au+WWW9i8eXO97//DP/wDS5cuLWpfd999N5/+9Kcb3aYNGzZw2223NVjv61//eqOPJUmSBIdQMK49b6gpXsXMPWrXrh0//elPWbduXak/XoN27NhRsn3vKxjv2LGD73//+/Tv379kx69LY4NxZu72+G1JkqSGHDLBeNe8oaZ61fPU5N20adOGadOm8Z3vfKeO9rzG+PHjGTRoEOPHj+f1118H4NJLL+Uzn/kM73nPe+jTpw8PPvggUH1j8dGjR3PeeefRv39/PvWpT9UEtyuuuILKykrKy8u57rrrao5RVlbG1772NU4//XR+8pOfcMcddzBs2DAGDx7MBz/4wZowe+mll3LFFVcwbtw4+vTpw69//Ws+8YlP0K9fPy699NKa/T322GOMHDmSoUOHcsEFF7Bp0yZuvfVWVq9ezbhx4xg3bhxQ/aS6XU/omzdvHmPHjmXXA1h++ctfMnToUAYPHlzv/YdXrlzJxIkTOeWUU7jhhhtqyr/97W8zYMAABgwYwC233LLP8unTp/PKK69QUVHBNddcw5o1axg9ejQVFRUMGDCA3/zmN0yfPp0///nPVFRU8NGPfpQVK1bQr18/rrzySoYOHcrKlSv3eW7/6Z/+ieHDhzN8+HCWLVvW8C+EJEk6vGVms79OPfXU3NPSpUt3W4e9qjRKMfvr2LFj/vGPf8xevXrlhg0b8uabb87rrrsuMzM/8IEP5N13352ZmXfeeWdOmTIlMzMvueSS/NCHPpQ7duzIJUuW5Mknn5yZmXPmzMl27drlK6+8ktu3b8+zzjorf/KTn2Rm5vr16zMzc/v27TlmzJh87rnnMjOzV69e+c1vfrOmPevWratZ/tKXvpS33nprzTE/8pGP5M6dO/NnP/tZdu7cORctWpQ7duzIoUOH5n//939nVVVVnnHGGblp06bMzLzpppvyhhtuqDlOVVVVrXNDPvDAAzXrY8aMyaeffjrXrl2bPXv2zFdffXW3dtd211135Tvf+c5ct25dbt68OcvLy/Ppp5/OBQsW5IABA3LTpk25cePG7N+/fz777LP1li9fvjzLy8tr9jtz5sy88cYba87T22+/XdNHuyxfvjwjIufNm1dTtq9zu2t/99xzT06aNGmvz5K59++hJKnlqv3/9l3LY8aMyTFjxjRLe9QyAQuynkzqXSkacNRRR/Hxj3+cW2+9lQ4dOtSUz5s3j5/+9KcATJ06lS984Qs175177rm0atWK/v378+abb9aUDx8+nD59+gBw0UUX8dRTT/GhD32I2bNnM2vWLLZv386aNWtYunQpgwYNAuAjH/lIzfaLFy/my1/+Mhs2bGDTpk1MmDCh5r1zzjmHiGDgwIEcf/zxDBw4EIDy8nJWrFjBqlWrWLp0KaNGjQLgL3/5CyNHjqzzM7du3ZoPfvCDe5X/7ne/Y/To0fTu3RuALl261Ln9e9/7Xrp27QrA+eefz1NPPUVEcN5559GxY8ea8t/85jdkZp3lkydP3m2fw4YN4xOf+ATbtm3j3HPPpaKios5j9+rVixEjRtSs7+vcXnTRRTU/P//5z9e5P0mSdOQ4ZKZSNKfPfe5z3HnnnfzpT3+qt07Uui9Mu3btapZz131j9qiza3358uXMnDmTJ554gkWLFjFp0iS2bNlSU2dXYITqKRPf/e53ef7557nuuut2q7frmK1atdrt+K1atWL79u1kJu9973tZuHAhCxcuZOnSpdx55511fpb27dvTunXrvcozc6/P0NC52LVe+zzsuc9ijB49mieffJITTjiBqVOn8sMf/rDOerXPV0PntnY7i/lckiTp8GYwLkKXLl348Ic/vFuQfM973sP9998PwH333cfpp5/e4H7mz5/P8uXL2blzJw888ACnn346b7/9Nh07duToo4/mzTff5Be/+EW922/cuJEePXqwbds27rvvvv36DCNGjOC3v/1tzVzazZs38z//8z8AdO7cmY0bNza4j5EjR/LrX/+a5cuXA/DWW2/VWe/xxx/nrbfe4s9//jM/+9nPGDVqFKNHj+ZnP/sZmzdv5k9/+hMPPfQQZ5xxRr3le7bptddeo3v37lx++eVcdtllPPvsswC0bduWbdu21dmOhs7tAw88UPOzvtFzSZJ05HAqRZGuvvpqvvvd79as33rrrXziE5/g5ptvplu3btx1110N7mPkyJFMnz6d559/vuaLeK1atWLIkCGUl5fTp0+fmqkOdfnnf/5nTjvtNHr16sXAgQOLCrO7dOvWjbvvvpuLLrqIrVu3AnDjjTfy7ne/m2nTpnH22WfTo0cP5syZs899zJo1i/PPP5+dO3fSvXt3Hn/88b3qnX766UydOpVly5Zx8cUXU1lZCVSPeA8fPhyovgXckCFD9lk+atQoBgwYwNlnn82AAQO4+eabadu2LZ06daoZMZ42bRqDBg1i6NChzJgxY7d2DB48eJ/nduvWrZx22mns3LmTH//4x0WfS0mSdHiKYv+UXUqVlZW5664Hu7zwwgv069evZr2srLg7SRSrVy9YsaLp9teQuXPnMnPmTH7+858fvIOqXmVlZSxYsIDjjjtun/X2/D2UJLVcEbs/+S4Txo4dC1T/f1gCiIhnMrOyrvcOmRHjgxliJUmSdOQ5ZILxoW7s2LE1/2pV81vhv7QkSdIe/PKdJEmShMFYkiRJAgzGkiRJErAfwTgiWkfEf0fEzwvrXSLi8Yh4ufDz2Fp1r42IZRHxUkRMqH+vkiRJUsuwPyPGnwVeqLU+HXgiM/sCTxTWiYj+wIVAOTARuC0i9n6M2iGgdevWVFRUMGDAAC644AI2b97c3E0C4Otf//o+33//+9/Phg0bitrX9ddfz8yZMxvdphUrVvCjH/1on3U2bNjAbbfd1uhjSZIklUJRwTgiegKTgO/XKp4C3FNYvgc4t1b5/Zm5NTOXA8uA4Y1uaVlZ9U0Jm+pVVtbgITt06MDChQtZvHgxf/M3f8Ptt99edHN37Nhx4J+1AfUF48xk586dPProoxxzzDElO35dGhuMS3m+JEmSilHsiPEtwBeAnbXKjs/MNQCFn90L5ScAK2vVW1Uoa5zXXqu+U3dTvfbzaSFnnHFGzeOUzz33XE499VTKy8uZNWtWTZ1OnTrx1a9+ldNOO4158+bxta99jWHDhjFgwACmTZvGroepjB07ls9//vOMHj2afv368fTTT3P++efTt29fvvzlL9fs795772X48OFUVFTwyU9+kh07djB9+nT+/Oc/U1FRwUc/+lFWrFhBv379uPLKKxk6dCgrV66krKyMdevWAfDDH/6QQYMGMXjwYKZOnVrnZ3vuuec488wz6du3L3fccQdQHbKvueYaBgwYwMCBA2sen1xf+fTp0/nNb35DRUUF3/nOd1iyZElN2wcNGsTLL7/M9OnTeeWVV6ioqOCaa65h7ty5jBs3josvvpiBAwc2eG6vvvpqhg4dyvjx46mqqtqv/pMkSWpQZu7zBXwAuK2wPBb4eWF5wx71/rfw83vAx2qV3wl8sI79TgMWAAtOOumk3NPSpUt3L4C96jRKEfvr2LFjZmZu27YtJ0+enLfddltmZq5fvz4zMzdv3pzl5eW5bt26wi7JBx54oGb7XfUyMz/2sY/lI488kpmZY8aMyS984QuZmXnLLbdkjx49cvXq1blly5Y84YQTct26dbl06dL8wAc+kH/5y18yM/OKK67Ie+65Z7d2ZWYuX748IyLnzZtXU9arV6+sqqrKxYsX57vf/e6sqqraqz27XHfddTlo0KDcvHlzVlVVZc+ePfONN97IBx98MM8666zcvn17/uEPf8gTTzwxV69eXW/5nDlzctKkSTX7/fSnP5333ntvZmZu3bo1N2/enMuXL8/y8vKaOnPmzMl3vOMd+eqrr+51zuo6t7v2d8MNN+RVV121j55rOnv9HkqSWqza/2vftTxmzJgcM2ZMs7RHLROwIOvJvcU84GMUMDki3g+0B46KiHuBNyOiR2auiYgewNpC/VXAibW27wmsriOQzwJmQfUjoYtox0G3a2QWqkeML7vsMgBuvfVWHnroIQBWrlzJyy+/TNeuXWndujUf/OAHa7afM2cO3/rWt9i8eTNvvfUW5eXlnHPOOQBMnjwZgIEDB1JeXk6PHj0A6NOnDytXruSpp57imWeeYdiwYTVt6d69O3Xp1asXI0aM2Kv8V7/6FR/60IdqHnvcpUuXOrefMmUKHTp0oEOHDowbN4758+fz1FNPcdFFF9G6dWuOP/54xowZw9NPP11v+VFHHbXbPkeOHMmMGTNYtWpVzWh4XYYPH07v3r1r1us7t61ateIjH/kIAB/72Mc4//zz69yfJEnSgWowGGfmtcC1ABExFvjHzPxYRNwMXALcVPj5cGGTR4AfRcS3gXcBfYH5Td/00ts1x7i2uXPn8p//+Z/MmzePd7zjHYwdO5YtW7YA0L59e1q3rv6e4ZYtW7jyyitZsGABJ554Itdff31NPYB27doB0KpVq5rlXevbt28nM7nkkkv4xje+0WA7O3bsWGd5ZhIRDW6/Z52IqJn2Udc+i3HxxRdz2mmn8R//8R9MmDCB73//+/Tp02everXbvq9z21CbJUmSGqsx9zG+CXhvRLwMvLewTmYuAWYDS4FfAldl5mHzzao//vGPHHvssbzjHe/gxRdf5He/+12d9XYFuuOOO45Nmzbx4IMP7tdxxo8fz4MPPsjatdUD8W+99RavFeZFt23blm3bthW1j9mzZ7N+/fqafdTl4YcfZsuWLaxfv565c+cybNgwRo8ezQMPPMCOHTuoqqriySefZPjw4fWWd+7cmY0bN9bs89VXX6VPnz585jOfYfLkySxatGivOnva17nduXNnzTn80Y9+xOmnn97g55ckSdofxUylqJGZc4G5heX1wPh66s0AZjSybS3SxIkTuf322xk0aBCnnHJKnVMYAI455hguv/xyBg4cSFlZWc2UiGL179+fG2+8kfe9733s3LmTtm3b8r3vfY9evXoxbdo0Bg0axNChQ5kxo/7TXF5ezpe+9CXGjBlD69atGTJkCHffffde9YYPH86kSZN4/fXX+cpXvsK73vUuzjvvPObNm8fgwYOJCL71rW/xzne+s97yrl270qZNGwYPHsyll17Kli1buPfee2nbti3vfOc7+epXv0qXLl0YNWoUAwYM4Oyzz2bSpElFn9uOHTuyZMkSTj31VI4++uiaL/1JkiQ1lSj2T+OlVFlZmQsWLNit7IUXXqBfv35/LSgr2+87SexTr16wYkXT7U8l1alTJzZt2nTQj7vX76EkqcWKqL7xVO3lsWPHAtXT9SSAiHgmMyvrem+/RoyblSFWkiRJJdSYOcbSQdMco8WSJOnIYjCWJEmSaOHBuCXMf9aRy98/SZKOLC02GLdv357169cbTtQsMpP169fTvn375m6KJEk6SFrsl+969uzJqlWrqKqqau6m6AjVvn17evbs2dzNkCRJB0mLDcZt27bd7VHBkiRJUim12KkUkiRJ0sFkMJYkSZIwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJKCIYBwR7SNifkQ8FxFLIuKGQvn1EfFGRCwsvN5fa5trI2JZRLwUERNK+QEkSZKkptCmiDpbgTMzc1NEtAWeiohfFN77TmbOrF05IvoDFwLlwLuA/4yId2fmjqZsuCRJktSUGhwxzmqbCqttC6/cxyZTgPszc2tmLgeWAcMb3VJJkiSphIqaYxwRrSNiIbAWeDwzf19469MRsSgifhARxxbKTgBW1tp8VaFMkiRJarGKCsaZuSMzK4CewPCIGAD8K3AyUAGsAf6lUD3q2sWeBRExLSIWRMSCqqqqA2q8JEmS1FT2664UmbkBmAtMzMw3C4F5J3AHf50usQo4sdZmPYHVdexrVmZWZmZlt27dDqjxkiRJUlMp5q4U3SLimMJyB+As4MWI6FGr2nnA4sLyI8CFEdEuInoDfYH5TdtsSZIkqWkVc1eKHsA9EdGa6iA9OzN/HhH/X0RUUD1NYgXwSYDMXBIRs4GlwHbgKu9IIUmSpJauwWCcmYuAIXWUT93HNjOAGY1rmiRJknTw+OQ7SZIkCYOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSgCKCcUS0j4j5EfFcRCyJiBsK5V0i4vGIeLnw89ha21wbEcsi4qWImFDKDyBJkiQ1hWJGjLcCZ2bmYKACmBgRI4DpwBOZ2Rd4orBORPQHLgTKgYnAbRHRuhSNlyRJkppKg8E4q20qrLYtvBKYAtxTKL8HOLewPAW4PzO3ZuZyYBkwvElbLUmSJDWxouYYR0TriFgIrAUez8zfA8dn5hqAws/uheonACtrbb6qULbnPqdFxIKIWFBVVdWYzyBJkiQ1WlHBODN3ZGYF0BMYHhED9lE96tpFHfuclZmVmVnZrVu34lorSZIklch+3ZUiMzcAc6meO/xmRPQAKPxcW6i2Cjix1mY9gdWNbqkkSZJUQsXclaJbRBxTWO4AnAW8CDwCXFKodgnwcGH5EeDCiGgXEb2BvsD8pm64JEmS1JTaFFGnB3BP4c4SrYDZmfnziJgHzI6Iy4DXgQsAMnNJRMwGlgLbgasyc0dpmi9JkiQ1jQaDcWYuAobUUb4eGF/PNjOAGY1unSRJknSQ+OQ7SZIkCYOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSUEQwjogTI2JORLwQEUsi4rOF8usj4o2IWFh4vb/WNtdGxLKIeCkiJpTyA0iSJElNoU0RdbYDV2fmsxHRGXgmIh4vvPedzJxZu3JE9AcuBMqBdwH/GRHvzswdTdlwSZIkqSk1OGKcmWsy89nC8kbgBeCEfWwyBbg/M7dm5nJgGTC8KRorSZIklcp+zTGOiDJgCPD7QtGnI2JRRPwgIo4tlJ0ArKy12SrqCNIRMS0iFkTEgqqqqv1uuCRJktSUig7GEdEJ+Hfgc5n5NvCvwMlABbAG+JddVevYPPcqyJyVmZWZWdmtW7f9brgkSZLUlIoKxhHRlupQfF9m/hQgM9/MzB2ZuRO4g79Ol1gFnFhr857A6qZrsiRJktT0irkrRQB3Ai9k5rdrlfeoVe08YHFh+RHgwohoFxG9gb7A/KZrsiRJktT0irkrxShgKvB8RCwslH0RuCgiKqieJrEC+CRAZi6JiNnAUqrvaHGVd6SQJElSS9dgMM7Mp6h73vCj+9hmBjCjEe2SJEmSDiqffCdJkiRhMJYkSZIAg7EkSZIEGIwlSZIkwGAsSZIkAQZjSZIkCTAYS5IkSYDBWJIkSQIMxpIkSRJgMJYkSZIAg7EkSZIEGIwlSZIkwGAsSZIkAQZjSZIkCTAYS5IkSYDBWJIkSQIMxpIkSRJgMJYkSZIAg7EkSZIEGDFIb9wAABBVSURBVIwlSZIkwGAsSZIkAQZjSZIkCTAYS5IkSYDBWJIkSQIMxpIkSRJgMJYkSZIAg7EkSZIEGIwlSZIkwGAsSZIkAQZjSZIkCTAYS5IkSYDBWJIkSQIMxpIkSRJgMJYkSZIAg7EkSZIEGIwlSZIkwGAsSZIkAQZjSZIkCTAYS5IkSYDBWJIkSQKKCMYRcWJEzImIFyJiSUR8tlDeJSIej4iXCz+PrbXNtRGxLCJeiogJpfwAkiRJUlMoZsR4O3B1ZvYDRgBXRUR/YDrwRGb2BZ4orFN470KgHJgI3BYRrUvReEmSJKmpNBiMM3NNZj5bWN4IvACcAEwB7ilUuwc4t7A8Bbg/M7dm5nJgGTC8qRsuSZIkNaX9mmMcEWXAEOD3wPGZuQaqwzPQvVDtBGBlrc1WFcr23Ne0iFgQEQuqqqr2v+WSJElSEyo6GEdEJ+Dfgc9l5tv7qlpHWe5VkDkrMyszs7Jbt27FNkOSJEkqiaKCcUS0pToU35eZPy0UvxkRPQrv9wDWFspXASfW2rwnsLppmitJkiSVRjF3pQjgTuCFzPx2rbceAS4pLF8CPFyr/MKIaBcRvYG+wPyma7IkSZLU9NoUUWcUMBV4PiIWFsq+CNwEzI6Iy4DXgQsAMnNJRMwGllJ9R4urMnNHk7dckiRJakINBuPMfIq65w0DjK9nmxnAjEa0S5IkSTqofPKdJEmShMFYkiRJAgzGkiRJEmAwliRJkgCDsSRJkgQYjCVJkiTAYCxJkiQBBmNJkiQJMBhLkiRJgMFYkiRJAgzGkiRJEmAwliRJkgCDsSRJkgQYjCVJkiTAYCxJkiQBBmNJkiQJMBhLkiRJgMFYkiRJAgzGkiRJEmAwliRJkgCDsSRJkgQYjCVJkiTAYCxJkiQBBmNJkiQJMBhLkiRJgMFYkiRJAgzGkiRJEmAwliRJkgCDsSRJkgQYjCVJkiTAYCxJkiQBBmNJkiQJMBhLkiRJgMFYkiRJAgzGkiRJEmAwliRJkgCDsSRJkgQYjCVJkiTAYCxJ0mGjrAwiql9lZc3dGunQ06a5GyBJkprGa69BZvVyRPO2RToUNThiHBE/iIi1EbG4Vtn1EfFGRCwsvN5f671rI2JZRLwUERNK1XBJkiSpKRUzleJuYGId5d/JzIrC61GAiOgPXAiUF7a5LSJaN1VjJUmSpFJpMBhn5pPAW0Xubwpwf2ZuzczlwDJgeCPaJ0mSJB0Ujfny3acjYlFhqsWxhbITgJW16qwqlEmSJEkt2oEG438FTgYqgDXAvxTK65rqn3XtICKmRcSCiFhQVVV1gM2QJEmSmsYBBePMfDMzd2TmTuAO/jpdYhVwYq2qPYHV9exjVmZWZmZlt27dDqQZkiRJUpM5oGAcET1qrZ4H7LpjxSPAhRHRLiJ6A32B+Y1roiRJklR6Dd7HOCJ+DIwFjouIVcB1wNiIqKB6msQK4JMAmbkkImYDS4HtwFWZuaM0TZckSZKaToPBODMvqqP4zn3UnwHMaEyjJEmSpIPNR0JLkiRJGIwlSZIkwGAsSZIkAQZjSZIkCTAYS5IkSYDBWJIkSQIMxpIkSRJgMJYkSZIAg7EkSZIEGIwlSZIkwGAsSZIkAQZjSZIkCTAYS5IkSYDBWJIkSQIMxpIkSRJgMJYkSZIAg7EkSZIEGIwlSZIkwGAsSZIkAQZjSZIkCTAYS5IkSYDBWJIkSQIMxpIkSRJgMJYkSZIAg7EkSZIEGIwlSZIkwGAsSZIkAQZjSZIkCTAYS5IkSYDBWJIkSQIMxpIkSRJgMJYkSZIAg7EkSZIEGIwlSZIkwGAsSZIkAQZjSZIkCTAYS5IkSYDBWJIkSQIMxpIkSRJgMJYkSZKAIoJxRPwgItZGxOJaZV0i4vGIeLnw89ha710bEcsi4qWImFCqhkuSJElNqZgR47uBiXuUTQeeyMy+wBOFdSKiP3AhUF7Y5raIaN1krZUkSZJKpMFgnJlPAm/tUTwFuKewfA9wbq3y+zNza2YuB5YBw5uorZIkSVLJHOgc4+Mzcw1A4Wf3QvkJwMpa9VYVyiRJkqQWram/fBd1lGWdFSOmRcSCiFhQVVXVxM2QJEmS9s+BBuM3I6IHQOHn2kL5KuDEWvV6Aqvr2kFmzsrMysys7Nat2wE2Q5IkSWoaBxqMHwEuKSxfAjxcq/zCiGgXEb2BvsD8xjVRkiRJKr02DVWIiB8DY4HjImIVcB1wEzA7Ii4DXgcuAMjMJRExG1gKbAeuyswdJWq7JEmS1GQaDMaZeVE9b42vp/4MYEZjGiVJkiQdbD75TpIkScJgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJEmSAIOxJEmSBBiMJUmSJMBgLEmSJAEGY0mSJAkwGEuSJEmAwViSJEkCDMaSJEkSYDCWJOmgKSuDiOpXWVlzt0bSnto0dwMkSTpSvPYaZFYvRzRvWyTtrVHBOCJWABuBHcD2zKyMiC7AA0AZsAL4cGb+b+OaKUmSJJVWU0ylGJeZFZlZWVifDjyRmX2BJwrrkiRJUotWijnGU4B7Csv3AOeW4BiSJElSk2psME7gsYh4JiKmFcqOz8w1AIWf3Rt5DEmSJKnkGvvlu1GZuToiugOPR8SLxW5YCNLTAE466aRGNkOSJElqnEaNGGfm6sLPtcBDwHDgzYjoAVD4ubaebWdlZmVmVnbr1q0xzZAkSZIa7YCDcUR0jIjOu5aB9wGLgUeASwrVLgEebmwjJUmSpFJrzFSK44GHovpGjG2AH2XmLyPiaWB2RFwGvA5c0PhmSpIk1aGsrPoG0cByelF9p1jpwBxwMM7MV4HBdZSvB8Y3plGSJElFqfXUlDKfmqJG8pHQkiRJEgZjSZIkCTAYS5KkQ01ZGURUv3r1au7W6DDS2PsYS5IkHVy15hVLTckRY0mSJAmDsSRJOhQ4fUIHgVMpJElSy+f0CR0EjhhLkqSWyVFiHWSOGEuSpJbJUWIdZI4YS5IkSRiMJUlSS+L0CTUjp1JIkqSWw+kTakaOGEuSpOblKLFaCEeMJUlS83KUWC2EI8aSJEkSBmNJktQcnD6hFsipFJIk6eBz+oRaIEeMJUnSweEosVo4R4wlSdLB4SixWjhHjCVJUukcxFHiFfSqOdZyykp6LB2eHDGWJEmlcxBHiXuzouZQZREH5Zg6vDhiLEmSmpZziXWIMhhLkqTGqx2GoXqUOBNWrGjOVkn7xakUkiSp8fxinQ4DjhhLknQAag+QRlSvH3GcMqHDjCPGkiQdgD0HSI+Y73qVlVV/eKgOw44S6zDiiLEkSdo35w/rCOGIsSRJ2jfnD+sI4YixJEnam/OHdQQyGEuSpL2/TQhOmdARx2AsSdKRqr65w4ZhHaEMxpIkHUn8Ip1UL4OxJEmHO8OwVBSDsSRJh6HllBmGpf1kMJYk6VBWazQ48ctzUmMYjCVJOtTUMzUi+OuX53qzojlbKB2SDMaSJLVUe95CzdFgqaQMxpKkFqt2Liwra+7WlFAxAdhbqUkl5yOhJUktVu0nEe/KiYey6i/EvQZAAuz6TL16+chlqQUwGEuS1ARqh976/TUAR5iFpZamZFMpImJiRLwUEcsiYnqpjiNJUknVM81htztA7GvaQ/qFOOlQUZJgHBGtge8BZwP9gYsion8pjiVJ0n6pbz5vfS+oM+TWvgOEoVc6PJRqxHg4sCwzX83MvwD3A1NKdCxJalGOmC+MlUo99+Xda4T2QF/Q4MiuX3STjkylCsYnACtrra8qlElS0Q7VgLnrC2OZ1cslt78joC39BXWOyu45QnvAL4OupHqU6st3dX13eLevGETENGBaYXVTRLxUorY05DhgXTMdW3WzT1qmZu2X1147tO5KULutJWz34Xmt7NHZB+lcHpA921NY32e/lPrztOTzdTDs9pn/unJcRBx+18qhrTn/+9WrvjdKFYxXASfWWu8JrK5dITNnAbNKdPyiRcSCzKxs7nbor+yTlsl+aXnsk5bJfml57JOWp6X2SammUjwN9I2I3hHxN8CFwCMlOpYkSZLUaCUZMc7M7RHxaeD/Aq2BH2TmklIcS5IkSWoKJXvAR2Y+Cjxaqv03oWafzqG92Cctk/3S8tgnLZP90vLYJy1Pi+yTSB+7I0mSJJXuyXeSJEnSoeSwDcYNPZI6Iq6JiIWF1+KI2BERXQrvrYiI5wvvLTj4rT98FdEvR0fE/x8Rz0XEkoj4+2K31YFpZJ94rZRIEf1ybEQ8FBGLImJ+RAwodlsdmEb2iddKCUTEDyJibUQsruf9iIhbC322KCKG1nrP66QEGtknzX+dZOZh96L6C3+vAH2AvwGeA/rvo/45wK9qra8Ajmvuz3G4vYrpF+CLwDcLy92Atwp196tPfZW+TwrrXivN1y83A9cVlv8OeKLYbX0d3D4prHutlKZfRgNDgcX1vP9+4BdUP19hBPD7YvvT18Htk8J7zX6dHK4jxvv7SOqLgB8flJYd2YrplwQ6R0QAnagOYduL3Fb7rzF9otIppl/6A08AZOaLQFlEHF/kttp/jekTlUhmPkn1f5PqMwX4YVb7HXBMRPTA66RkGtEnLcLhGoyLfiR1RLwDmAj8e63iBB6LiGcKT+hT0yimX74L9KP6gTDPA5/NzJ1Fbqv915g+Aa+VUimmX54DzgeIiOFUP8mpZ5Hbav81pk/Aa6W51NdvXifNZ1/nvtmvk5Ldrq2ZNfhI6lrOAX6bmbX/dTMqM1dHRHfg8Yh4sfAvIDVOMf0yAVgInAmcTPX5/02R22r/HXCfZObbeK2USjH9chPwfyJiIdX/YPlvqkfyvVZKozF9Al4rzaW+fvM6aT77OvfNfp0criPGDT6SupYL2WMaRWauLvxcCzxE9Z9c1HjF9MvfAz8t/IllGbCc6rl6+9OnKl5j+sRrpXQa7JfMfDsz/z4zK4CPUz3/e3kx2+qANKZPvFaaT3395nXSfOo99y3hOjlcg3FRj6SOiKOBMcDDtco6RkTnXcvA+4A6v1mp/VZMv7wOjAcozM07BXi1yG21/w64T7xWSqrBfomIYwrvAfwD8GRhFN9rpTQOuE+8VprVI8DHC3dCGAH8MTPX4HXSnOrsk5ZynRyWUymynkdSR8SnCu/fXqh6HvBYZv6p1ubHAw9Vf8+INsCPMvOXB6/1h68i++Wfgbsj4nmq/9zyT5m5DqCubZvjcxxOGtMnEdEHr5WSKLJf+gE/jIgdwFLgsn1t2xyf43DSmD7B/6+UTET8GBgLHBcRq4DrgLZQ0yePUn0XhGXAZqr/AuZ1UkIH2ie0kOvEJ99JkiRJHL5TKSRJkqT9YjCWJEmSMBhLkiRJgMFYkiRJAgzGkiRJEmAwliRJkgCDsSRJkgQYjCVJkiQA/h/iAMmkP2xe8AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pyplot.figure(figsize=(12, 8))\n",
    "pyplot.hist(t_boot_nonparam, bins, label='Nonparametric bootstrap', color='blue', histtype='step', density=True)\n",
    "pyplot.hist(t_boot_param, bins, label='Parametric bootstrap', color='red', histtype='step', density=True)\n",
    "pyplot.axvline(x=1, color='black', label='θ (true parameter)')\n",
    "pyplot.legend(loc='upper left')\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAHWCAYAAACFXRQ+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5TdVX0//PeGgKQkAnIxlqBERK4JJyGEULBELgaMEKvWSxSKqxrAS8Xn0RZYWqmlS1ywHm3Xg8ZIvTNa2+qCVqpUcaSVIhAcCBd5DMglQjSmchO5BPbzx0zmN/lmkpyQmZy5vF5rzXLO+e7vOZ9zdg6+88k++1tqrQEAAP6P7TpdAAAAjDRCMgAANAjJAADQICQDAECDkAwAAA1CMgAANLQVkkspJ5VS7iqlrCilnDvI8QNLKf9TSnmqlPKhLTkXAABGmrK5fZJLKdsn+f+SnJhkZZIbk7yt1nrHgDF7JXlZktcn+W2t9ZJ2zwUAgJGmnU7ynCQraq331FqfTvKNJAsHDqi1/rrWemOSZ7b0XAAAGGnaCcl7J3lgwO2Vffe1Y2vOBQCAjpjQxpgyyH3tXsu67XNLKYuTLE6SnXfe+fADDzywzacAANi8u+66K0lywAEHdLgSRoply5b9pta652DH2gnJK5PsM+D21CQPtvncbZ9ba12aZGmSzJ49u950001tPgUAwObNmzcvSdLd3d3ROhg5Sin3bexYO8stbkyyfyllWillxyRvTXJlm8+9NecCAEBHbLaTXGtdW0p5X5LvJdk+yRdqrbeXUs7qO76klDIlyU1JXpjkuVLKOUkOrrU+Oti5w/ViAABgKLSz3CK11quSXNW4b8mA31eldylFW+cCAMBI1lZIHgmeeeaZrFy5Mk8++WSnS4F+O+20U6ZOnZoddtih06UAAENo1ITklStXZvLkydl3331TymCbZsC2VWvNmjVrsnLlykybNq3T5QAAQ6ity1KPBE8++WR23313AZkRo5SS3Xff3b9uAMAYNGpCchIBmRHHn0kAGJtGVUjulDVr1qTVaqXVamXKlCnZe++9+28//fTTnS7veVm7dm123XXXJMkDDzyQt7zlLcPyPFOnTs3DDz+cZ599Nq961as2Ou65557LRRddtMnHmj9/fh577LGsWLEirVZri+q45pprcv311/ffvvTSS3P55Zdv0WMAAOPHqFmT3Em77757enp6kiQXXHBBJk2alA996EPrjam1ptaa7bYbfX/v2GefffJP//RPw/oc22+/ff7rv/5ro8fXheRzzz13g2Pr3tvvfe97SZJf/epXW/z811xzTfbYY4/MnTs3SfLe9753ix8DABg/Rl+iG0FWrFiRQw89NGeddVZmzZqVBx54oL87myTf+MY38q53vStJb7B7wxvekNmzZ2fOnDnrdTXXWb58eY444oi0Wq3MmDEj99xzT5LklFNOyeGHH55DDjkkl112WZL/0wn+8Ic/nFmzZmX+/Pn5yU9+kmOPPTYvf/nLc9VVvbvuXXbZZfmTP/mTzJ8/PwcccEAuvPDCQV/Hus7sZZddlje96U2ZP39+9t9//5x33nn94z73uc/lla98ZebNm5d3vetdOeecczZ4rNWrV+fEE0/MrFmzcvbZZ6fWul69SfLLX/4yxxxzTFqtVg499NBcd911Offcc/PYY4+l1Wrl9NNP3+C9feihh/q70knvbiennXZapk+fnje/+c35/e9/nyTrjbn++utzwgkn5O67785ll12Wiy++OK1WK9ddd10+8pGP5NOf/nSS5Oabb86RRx6ZGTNm5I1vfGMeeeSRJMkxxxyTc889N3PmzMkBBxyQ6667ro0/FQDAWDAqO8nnfPec9KzqGdLHbE1p5dMnfXqLz7vjjjvyxS9+MUuWLMnatWs3Ou4v/uIv8pd/+ZeZO3du7r333rzuda/Lbbfdtt6Yz3zmM/nQhz6Ut7zlLXnqqaf6A+aXv/zlvOhFL8oTTzyR2bNn541vfGMmT56cRx55JK95zWty8cUX55RTTskFF1yQH/zgB7nlllty5pln5rWvfW2S5IYbbshtt92WHXfcMUcccURe97rX5dBDD91orbfccktuvvnmTJgwIa985Svz/ve/P88++2wuuuii3Hzzzdl5550zb968zJkzZ4NzP/axj+XVr351zj///FxxxRVZsmTJBmO+9rWv5ZRTTslf/dVf5dlnn83vf//7zJkzJ5dddll/x37FihXrvbeDve//+I//mLlz5+b000/P5z73uUFDe5Lst99+ede73pU99tijf8y6v0QkyTve8Y4sXbo0xxxzTM4///z87d/+bS655JIkvV3sG264IVdeeWU+/vGP57vf/e5G3zcAYOwYlSF5JNlvv/1yxBFHbHbc97///dx11139t3/729/m97//fSZOnNh/3x/90R/lwgsvzH333Zc3vOENecUrXpEk+dSnPpUrr+y9mvfKlStz9913p9VqZeLEiTnxxBOTJNOnT88uu+ySCRMmZPr06bn33nv7H3f+/PnZbbfdkiSvf/3r89///d+bDMknnHBCJk+enCQ58MADc//992flypU57rjj+h/nTW96U+6///4Nzr322mv7A+jChQv7H2egI444ImeeeWaefPLJvP71r89hhx026F8wNvXeTps2rX/pxLqQu7GQvClr1qzJk08+mWOOOSZJ8md/9mc57bTT+o+/4Q1vSJIcfvjh672nAMDYNipD8vPp+A6XnXfeuf/37bbbrr/7m2S9rcHWdSR33HHHjT7WaaedlqOOOirf+c53cuKJJ+bLX/5ynn766Vx77bW5/vrrM3HixBxzzDH9jzvwsbbbbru84AUv6P99YOhs7sCwuR0Z1j1O0ruWeO3ateu9rs3Z3OMfd9xx6e7uzne+8528/e1vz3nnnTfoFwcHvrebe451tydMmJDnnnsuSdramm1zr2vde7HufQAAxgdrkofQdtttl9122y0///nP89xzz+Xb3/52/7ETTjghl156af/tdcsKBrrnnnvyile8Ih/4wAeyYMGC3HrrrXnkkUfyohe9KBMnTsztt9+eG2+8cYvruvrqq/Pwww/niSeeyBVXXJGjjz56ix/jyCOPzA9/+MM8/PDDeeaZZ/Ktb31r0HF//Md/3L9rxL/927/lscce22DMfffdlylTpmTx4sU544wz8tOf/jQTJvT+fa3dIPqLX/yi/734+te/3t8J3nfffbNs2bIkyb/+67/2j588efKgteyxxx6ZOHFi/3rjr371qzn22GPbqgEAGLuE5CH2yU9+MieddFKOP/74TJ06tf/+Sy+9ND/+8Y8zY8aMHHzwwfn85z+/wbldXV055JBD0mq1cs899+Qd73hHFixYkCeeeCKHHXZYPv7xj+fII4/c4pqOOeaYLFq0KDNnzszb3va2Ld4+LUle+tKX5sMf/nDmzJmT17zmNTnkkEOyyy67bDDub/7mb/L9738/s2bNSnd3d/bee+8NxvzgBz/IYYcdlpkzZ+aKK67I+9///iTJn//5n2fGjBk5/fTTN1vPIYccks9//vOZMWNGfve732Xx4sVJencfec973pNXvepV63XaFy5cmG9+85uZOXPmBl/A++pXv5oPfvCDmTFjRu6444585CMf2aL3BgAYe8qW/DP6tjJ79ux60003rXffnXfemYMOOqhDFY1el112WW677bb+nRy2xuOPP55JkyblmWeeycKFC3P22WfnlFNOGYIqRzd/NgFGh3nz5iVJuru7O1oHI0cpZVmtdfZgx3SSadtHP/rRzJw5MzNmzMgBBxyQ173udZ0uCQBgWIzKL+7RvnX7NA+FT33qU0P2WAAAI5lOMgAANAjJAADQICQDAECDkAwAAA1C8hb45je/mYMOOigHHXRQ/vmf/3nQMbXWHHfccXn00Ufz8MMP5zOf+cw2rrLztuZ1n3DCCfntb387xBUBAGwZIblNjz32WD75yU/mxhtvzA033JBPfOITefzxxzcYd9VVV+Wwww7LC1/4wk2GxWeffXa4S96kWmv/5ZuH2vMJyevqOe2008blXywAgJFFSG7Td7/73Zx88smZNGlSJk+enBNOOCH/8R//scG4yy+/PAsXLkySnHvuubn77rvTarXy4Q9/ON3d3Xn1q1+dRYsWZfr06bn33ntz6KGH9p97ySWX5IILLkiS3H333TnppJNy+OGH51WvelV+9rOfbfBcF1xwQU477bQcd9xx2X///fuv4vf444/n+OOPz6xZszJ9+vRcccUVSZJ77703Bx10UN7znvdk1qxZeeCBB3L22Wdn9uzZOeSQQ/Kxj32s/7H33XffnH/++TnqqKMye/bs3HzzzZk/f37222+/LFmypH/cxRdfnCOOOCIzZszoP7/5ujc2brB6Tj311Hz9619/3vMEADAURuU+yeecc056enqG9DFbrdYmr0p3//3354tf/GL+/d//PUmyevXqvPjFL95g3I9//ON87nOfS5JcdNFFue222/pr7e7uzg033JDbbrst06ZNy7333rvR51u8eHGWLFmS/fffPz/5yU/ynve8J9dcc80G42699dZcf/31+d3vfpeZM2dmwYIF2WuvvfLtb387L3zhC/Ob3/wmc+fOzamnnpokueuuu/LFL36xv1v7d3/3d3nRi16UZ599Nscff3xuvfXWzJgxI0myzz775H/+53/ywQ9+MGeccUZ+/OMf58knn8whhxySs846K1dffXV+/vOf54YbbkitNaeeemquvfbaDV73xsa99KUv3aCeJHnqqaeyZs2a7L777ht9fwAAhtOoDMmdUGvN2WefnY985CNJko9//OMppWww7n//938zefLkjT7OnDlzMm3atE0+1+OPP57rrrsuf/qnf9p/31NPPTXo2IULF2bixImZOHFiXv3qV+eGG27IggULcv755+faa6/Ndtttl1/+8pf51a9+lSR52ctelrlz5/af/81vfjNLly7N2rVr89BDD+WOO+7oD8nrgvX06dPz+OOPZ/LkyZk8eXJ22mmnPPzww7n66qtz9dVXZ+bMmf11//znP89LX/rS9Wrc1LhmPUmy11575cEHHxSSAYCOGZUheVMd3+Hyspe9bL1O7gMPPJD58+dvMG7ChAl57rnnst12g69k2XnnnTcYu86TTz6ZJHnuueey6667ttUtbwb1Ukouv/zyrF69OsuWLcsOO+yQfffdt/+xBz7/L37xi1xyySW58cYbs9tuu+WMM87oH5ckL3jBC5Ik2223Xf/v626vXbs2tdacd955OfPMM9erodkh39S4gfUMfB8mTpy42dcOADBcrElu08knn5wf/vCHWbNmTX7zm9+ku7s7J5100gbjDjjggNxzzz1JksmTJ+exxx7b6GO++MUvzq9//eusWbMmTz31VP9Sjhe+8IWZNm1a/w4atdbccsstgz7GFVdckSeffDJr1qxJd3d3jjjiiDzyyCPZa6+9ssMOO+SHP/xh7rvvvkHPffTRR7Pzzjtnl112ya9+9atB11hvyvz58/OFL3yh/wuMv/zlL/PrX/96g9e9sXGDqbVm1apV2XfffbeoFgCAoTQqO8mdMGnSpFx44YU56qijUmvNJz7xiUyaNGmDcQsWLEh3d3de8YpXZPfdd8/RRx+dQw89NCeffHIWLFiw3tgddtghf/3Xf50jjzwy06ZNy4EHHth/7PLLL8/ZZ5+dCy+8MM8880ze+ta35rDDDtvg+ebMmZMFCxbk/vvvz0c/+tH84R/+Yd7+9rfnlFNOyezZs9NqtdZ73IEOO+ywzJw5M4ccckhe/vKX5+ijj96i9+Q1r3lN7rzzzhx11FH979HXvva17Lfffuu97osvvnjQcdtvv/0Gj7ls2bLMnTs3Eyb4owkAdE6ptXa6hg3Mnj273nTTTevdd+edd+aggw7qUEXte+ihh3L66afnP//zP4f9uS644IJMmjQpH/rQh4b9ubaVD3zgAzn11FNz/PHHd7qUto2WP5sA4928efOS9H6RHpKklLKs1jp7sGOWWwyxl7zkJXn3u9+dRx99tNOljEqHHnroqArIAMDY5N+0h8Gb3/zmbfI86/ZUHkve/e53d7oEAACdZAAAaBpVIXkkrp9mfPNnEgDGplETknfaaaesWbNGKGHEqLVmzZo12WmnnTpdCgAwxEbNmuSpU6dm5cqVWb16dadLgX477bRTpk6d2ukyAIAhNmpC8g477LDZyzkDAMBQGDXLLQAAYFsRkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCAhrZCcinlpFLKXaWUFaWUcwc5Xkop/9B3/NZSyqwBxz5YSrm9lHJbKeXrpZSdhvIFAADAUNtsSC6lbJ/k0iQnJzk4ydtKKQc3hp2cZP++n8VJPtt37t5J/iLJ7FrroUm2T/LWIaseAACGQTud5DlJVtRa76m1Pp3kG0kWNsYsTPKV2uv6JLuWUl7Sd2xCkomllAlJ/iDJg0NUOwAADIt2QvLeSR4YcHtl332bHVNr/WWSS5Lcn+ShJI/UWq9+/uUCAMDwayckl0Huq+2MKaXslt4u87Qkf5hk51LKOwZ9klIWl1JuKqXctHr16jbKAgCA4dFOSF6ZZJ8Bt6dmwyUTGxtzQpJf1FpX11qfSfKtJH802JPUWpfWWmfXWmfvueee7dYPAABDrp2QfGOS/Usp00opO6b3i3dXNsZcmeT0vl0u5qZ3WcVD6V1mMbeU8gellJLk+CR3DmH9AAAw5CZsbkCtdW0p5X1Jvpfe3Sm+UGu9vZRyVt/xJUmuSvLaJCuSPJHknX3HflJK+ZckNydZm+SnSZYOxwsBAIChstmQnCS11qvSG4QH3rdkwO81yXs3cu7HknxsK2oEAIBtyhX3AACgQUgGAIAGIRkAABqEZAAAaBCSAQCgQUgGAIAGIRkAABqEZAAAaBCSAQCgQUgGAIAGIRkAABqEZAAAaBCSAQCgQUgGAIAGIRkAABqEZAAAaBCSAQCgQUgGAIAGIRkAABqEZAAAaBCSAQCgQUgGAIAGIRkAABqEZAAAaBCSAQCgQUgGAIAGIRkAABqEZAAAaBCSAQCgQUgGAIAGIRkAABqEZAAAaBCSAQCgQUgGAIAGIRkAABqEZAAAaBCSAQCgQUgGAIAGIRkAABqEZAAAaBCSAQCgQUgGAIAGIRkAABqEZAAAaBCSAQCgYUKnCwAAGGpLly1N1/Ku/tuLpi/qYDWMRkIyADDmdC3vSs+qnrSmtNKzqqfT5TAKWW4BAIxJrSmtdJ/RndaUVqdLYRQSkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCAhgmdLgAAYCgsXbY0Xcu7kiQ9q3rSmtLqcEWMZjrJAMCY0LW8Kz2repIkrSmtLJq+qMMVMZrpJAMAY0ZrSivdZ3R3ugzGAJ1kAABoEJIBAKBBSAYAgAYhGQAAGoRkAABoEJIBAKBBSAYAgAYhGQAAGoRkAGDM61nV0/+zdNnSTpfDKNBWSC6lnFRKuauUsqKUcu4gx0sp5R/6jt9aSpk14NiupZR/KaX8rJRyZynlqKF8AQAAm7Jo+qK0prSSJI8//Xi6lnd1uCJGg82G5FLK9kkuTXJykoOTvK2UcnBj2MlJ9u/7WZzkswOO/X2S79ZaD0xyWJI7h6BuAIC2LD58cbrP6E5rSiuTdpzU6XIYJdrpJM9JsqLWek+t9ekk30iysDFmYZKv1F7XJ9m1lPKSUsoLk/xxkn9Mklrr07XWh4ewfgAAGHLthOS9kzww4PbKvvvaGfPyJKuTfLGU8tNSymWllJ23ol4AABh27YTkMsh9tc0xE5LMSvLZWuvMJL9LssGa5iQppSwupdxUSrlp9erVbZQFAADDo52QvDLJPgNuT03yYJtjViZZWWv9Sd/9/5Le0LyBWuvSWuvsWuvsPffcs53aAQBgWLQTkm9Msn8pZVopZcckb01yZWPMlUlO79vlYm6SR2qtD9VaVyV5oJRyQN+445PcMVTFAwDAcJiwuQG11rWllPcl+V6S7ZN8odZ6eynlrL7jS5JcleS1SVYkeSLJOwc8xPuTXN4XsO9pHAMAgBFnsyE5SWqtV6U3CA+8b8mA32uS927k3J4ks7eiRgCAQS1dtrR/3+OeVT39+yHD1nLFPQBg1Opa3pWeVT1JktaUVhZNX9Thihgr2uokAwCMVK0prXSf0d3pMhhjdJIBAKBBSAYAgAYhGQAAGoRkAABoEJIBAKBBSAYAgAYhGQAAGoRkAABoEJIBAKBBSAYAgAYhGQAAGoRkAABoEJIBAKBBSAYAgAYhGQAAGiZ0ugAAgHYtXbY0Xcu7+m/3rOpJa0qrgxUxVukkAwCjRtfyrvSs6um/3ZrSyqLpizpYEWOVTjIAMKq0prTSfUZ3p8tgjNNJBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAaJnS6AACATVm6bGm6lnclSXpW9aQ1pdXhihgPdJIBgBGta3lXelb1JElaU1pZNH1RhytiPNBJBgBGvNaUVrrP6O50GYwjOskAANAgJAMAQIOQDAAADUIyAAA0CMkAANAgJAMAQIOQDAAADUIyAAA0CMkAANAgJAMAQIOQDAAADUIyAAA0CMkAANAwodMFAAA0LV22NF3Lu5IkPat60prS6nBFjDc6yQDAiNO1vCs9q3qSJK0prSyavqjDFTHe6CQDACNSa0or3Wd0d7oMximdZAAAaBCSAQCgQUgGAIAGIRkAABqEZAAAaBCSAQCgQUgGAIAGIRkAABqEZAAAaBCSAQCgQUgGAIAGIRkAABomdLoAAIAkWbpsabqWdyVJelb1pDWl1eGKGM90kgGAEaFreVd6VvUkSVpTWlk0fVGHK2I800kGAEaM1pRWus/o7nQZoJMMAABNQjIAADQIyQAA0CAkAwBAg5AMAAANQjIAADQIyQAA0CAkAwBAg5AMAAANQjIAADS0FZJLKSeVUu4qpawopZw7yPFSSvmHvuO3llJmNY5vX0r5aSnl34eqcABg9Fu6bGnmfWle5n1pXnpW9XS6HOi32ZBcStk+yaVJTk5ycJK3lVIObgw7Ocn+fT+Lk3y2cfwDSe7c6moBgDGla3lXfzhuTWll0fRFHa4Iek1oY8ycJCtqrfckSSnlG0kWJrljwJiFSb5Sa61Jri+l7FpKeUmt9aFSytQkC5L8XZL/a2jLBwBGu9aUVrrP6O50GbCedpZb7J3kgQG3V/bd1+6YTyf5yyTPPc8aAQBgm2onJJdB7qvtjCmlvC7Jr2utyzb7JKUsLqXcVEq5afXq1W2UBQAAw6OdkLwyyT4Dbk9N8mCbY45Ocmop5d4k30hyXCnla4M9Sa11aa11dq119p577tlm+QAAMPTaCck3Jtm/lDKtlLJjkrcmubIx5sokp/ftcjE3ySO11odqrefVWqfWWvftO++aWus7hvIFAADAUNvsF/dqrWtLKe9L8r0k2yf5Qq319lLKWX3HlyS5Kslrk6xI8kSSdw5fyQAAMLza2d0itdar0huEB963ZMDvNcl7N/MY3Um6t7hCAADYxlxxDwAAGoRkAABoEJIBAKBBSAYAgIa2vrgHADBUli5bmq7lXUmSnlU9aU1pdbgi2JBOMgCwTXUt70rPqp4kSWtKK4umL+pwRbAhnWQAYJtrTWml+4zuTpcBG6WTDAAADUIyAAA0CMkAANAgJAMAQIOQDAAADUIyAAA0CMkAANAgJAMAQIOQDAAADa64BwAMu6XLlqZreVeSpGdVT1pTWh2uCDZNJxkAGHZdy7vSs6onSe8lqRdNX9ThimDTdJIBgG2iNaWV7jO6O10GtEUnGQAAGoRkAABoEJIBAKBBSAYAgAYhGQAAGoRkAABoEJIBAKBBSAYAgAYXEwEAhoVLUTOa6SQDAMPCpagZzXSSAYBh41LUjFY6yQAA0CAkAwBAg5AMAAANQjIAADQIyQAA0CAkAwBAg5AMAAAN9kkGAIaMq+wxVugkAwBDxlX2GCt0kgGAIeUqe4wFOskAANAgJAMAQIOQDAAADUIyAAA0CMkAANAgJAMAQIOQDAAADfZJBgC2iqvsMRbpJAMAW8VV9hiLdJIBgK3mKnuMNTrJAADQICQDAECDkAwAAA1CMgAANAjJAADQYHcLAGCL2RuZsU4nGQDYYvZGZqzTSQYAnhd7IzOW6SQDAECDkAwAAA1CMgAANAjJAADQICQDAECDkAwAAA22gAMA2uICIownOskAQFtcQITxRCcZAGibC4gwXugkAwBAg5AMAAANQjIAADQIyQAA0OCLewDARtn2jfFKJxkA2CjbvjFe6SQDAJtk2zfGI51kAABoEJIBAKBBSAYAgAYhGQAAGtoKyaWUk0opd5VSVpRSzh3keCml/EPf8VtLKbP67t+nlPLDUsqdpZTbSykfGOoXAAAMraXLlmbel+Zl3pfm9e9sAePNZkNyKWX7JJcmOTnJwUneVko5uDHs5CT79/0sTvLZvvvXJvm/a60HJZmb5L2DnAsAjCC2fYP2toCbk2RFrfWeJCmlfCPJwiR3DBizMMlXaq01yfWllF1LKS+ptT6U5KEkqbU+Vkq5M8nejXMBgBHGtm+Md+0st9g7yQMDbq/su2+LxpRS9k0yM8lPtrRIAADYltoJyWWQ++qWjCmlTEryr0nOqbU+OuiTlLK4lHJTKeWm1atXt1EWAAAMj3ZC8sok+wy4PTXJg6sjXHkAAAmzSURBVO2OKaXskN6AfHmt9Vsbe5Ja69Ja6+xa6+w999yzndoBAGBYtBOSb0yyfyllWillxyRvTXJlY8yVSU7v2+VibpJHaq0PlVJKkn9Mcmet9f8Z0soBAGCYbPaLe7XWtaWU9yX5XpLtk3yh1np7KeWsvuNLklyV5LVJViR5Isk7+04/OslpSZaXUtbtIXN+rfWqoX0ZAMDWWLpsabqWdyVJelb1pDWl1eGKoLPa2d0ifaH2qsZ9Swb8XpO8d5Dz/juDr1cGAEaQddu+taa0bPsGaTMkAwBjn23f4P9wWWoAAGgQkgEAoEFIBgCABmuSAWCcsqMFbJxOMgCMU+t2tEhiRwto0EkGgHHMjhYwOJ1kAABoEJIBAKDBcgsAGEd8WQ/ao5MMAOOIL+tBe3SSAWCc8WU92DydZAAAaBCSAQCgQUgGAIAGa5IBYAwbuJtFYkcLaJdOMgCMYQN3s0jsaAHt0kkGgDHObhaw5XSSAQCgQUgGAIAGyy0AYIxx6WnYejrJADDGuPQ0bD2dZAAYg3xZD7aOTjIAADToJAPAGGAdMgwtnWQAGAOsQ4ahpZMMAGOEdcgwdHSSAQCgQScZAEYp65Bh+OgkA8AoZR0yDB+dZAAYxaxDhuGhkwwAAA06yQAwiliHDNuGTjIAjCLWIcO2oZMMAKOMdcgw/HSSAQCgQScZAEY465Bh29NJBoARzjpk2PZ0kgFgFLAOGbYtIRkARiBLLKCzLLcAgBHIEgvoLJ1kABihLLGAzhGSAWCEsMQCRg7LLQBghLDEAkYOnWQAGEEssYCRQUgGgA6yxAJGJsstAKCDLLGAkUknGQA6zBILGHmEZADYxiyxgJHPcgsA2MYssYCRTycZALaBwbrHlljAyKWTDADbgO4xjC46yQAwTHSPYfTSSQaAYaJ7DKOXTjIADCPdYxidhGQAGEK2d4OxwXILABhClljA2KCTDABbyRf0YOwRkgHgeRgYjH9034+SJMe+7FjdYxgjhGQAeB7WLatoTWnl2Jcdm0XTF2Xx4Ys7XRYwRIRkAGiTZRUwfvjiHgC0yZfyYPzQSQaATdA9hvFJSAaAAQaG4sSX8mC8EpIBGPc2tlPFuv/1pTwYf4RkAMY9O1UATUIyAOOStcbApgjJAIwbLgACtEtIBmBM21gwtqwC2BQhGYAxRzAGtpaQDMCYIBgDQ0lIBmBUae5jvI5gDAwlIRmAEW9T+xivIxgDQ0lIBmDE0CUGRgohGYBtrp0wPJBgDGxrQjIAw0YYBkYrIRmArbKxIJwIw8DoJSQDsJ5Nhd7BbCwIr7tPGAZGo7ZCcinlpCR/n2T7JJfVWi9qHC99x1+b5IkkZ9Rab27nXACG11CG3sEIwsBYtNmQXErZPsmlSU5MsjLJjaWUK2utdwwYdnKS/ft+jkzy2SRHtnkuAH22NNC2Q+gF2HLtdJLnJFlRa70nSUop30iyMMnAoLswyVdqrTXJ9aWUXUspL0mybxvnAjxvwxEqO2lLA207hF6ALddOSN47yQMDbq9Mb7d4c2P2bvPcEeGc756TnlU9nS4D2ELDESo7SaCF4dezqifzvjSv02UwQGtKK58+6dOdLmM97YTkMsh9tc0x7Zzb+wClLE6y7v8VHi+l3NVGbUNtjyS/6cDzsnHmZGQacfPyo/yo0yUMmR/lRzkzZ27paSNuTkhiXkaiPXJefjOW/psxFvwoP9rj7/P3nfisvGxjB9oJySuT7DPg9tQkD7Y5Zsc2zk2S1FqXJlnaRj3DppRyU611didrYH3mZGQyLyOPORmZzMvIY05GppE4L9u1MebGJPuXUqaVUnZM8tYkVzbGXJnk9NJrbpJHaq0PtXkuAACMKJvtJNda15ZS3pfke+ndxu0LtdbbSyln9R1fkuSq9G7/tiK9W8C9c1PnDssrAQCAIdLWPsm11qvSG4QH3rdkwO81yXvbPXcE6+hyDwZlTkYm8zLymJORybyMPOZkZBpx81J68y0AALBOO2uSAQBgXBkXIbmUclIp5a5SyopSyrmDHP9wKaWn7+e2UsqzpZQX9R27t5SyvO/YTdu++rGrjXnZpZTyb6WUW0opt5dS3tnuuTw/WzknPivDpI152a2U8u1Syq2llBtKKYe2ey7Pz1bOic/KMCilfKGU8utSym0bOV5KKf/QN2e3llJmDTjmczJMtnJeOvtZqbWO6Z/0fmHw7iQvT++WdLckOXgT409Jcs2A2/cm2aPTr2Os/bQzL0nOT/LJvt/3TPK/fWO3aE79DP+c9N32WencvFyc5GN9vx+Y5Aftnutn285J322fleGZlz9OMivJbRs5/tok/5HeazjMTfKTdufTz7afl75jHf2sjIdOcv9ltWutTydZd2nsjXlbkq9vk8rGt3bmpSaZXEopSSalN5CtbfNcttzWzAnDp515OTjJD5Kk1vqzJPuWUl7c5rlsua2ZE4ZJrfXa9P43aWMWJvlK7XV9kl1LKS+Jz8mw2op56bjxEJI3dsnsDZRS/iDJSUn+dcDdNcnVpZRlpfeqgAyNdubl/01yUHovQLM8yQdqrc+1eS5bbmvmJPFZGS7tzMstSd6QJKWUOem9gtTUNs9ly23NnCQ+K52ysXnzOemsTb3/Hf2stLUF3CjX9qWx07vU4se11oF/4zm61vpgKWWvJP9ZSvlZ39+K2DrtzMv8JD1JjkuyX3rf//9q81y23POek1rro/FZGS7tzMtFSf6+lNKT3r+8/DS9HX6fleGxNXOS+Kx0ysbmzeekszb1/nf0szIeOsntXFZ7nbemsdSi1vpg3//+Osm30/vPMmy9dublnUm+1fdPMCuS/CK9a/u2ZE5p39bMic/K8NnsvNRaH621vrPW2kpyenrXi/+inXN5XrZmTnxWOmdj8+Zz0lkbff87/VkZDyG5rUtjl1J2SXJskisG3LdzKWXyut+TvCbJoN/OZIu1My/3Jzk+SfrW8h2Q5J42z2XLPe858VkZVpudl1LKrn3HkuRdSa7t6+77rAyP5z0nPisddWWS0/t2U5ib5JFa60PxOem0QedlJHxWxvxyi9reZbWT5E+SXF1r/d2A01+c5Nu931HKhCRdtdbvbrvqx6425+Vvk3yplLI8vf8c81e11t8kyWDnduJ1jCVbMyellJfHZ2VYtDkvByX5Sinl2SR3JPnzTZ3bidcxlmzNnMT/rwybUsrXk8xLskcpZWWSjyXZIemfk6vSu5PCiiRPpPdfxnxOhtnznZeMgM+KK+4BAEDDeFhuAQAAW0RIBgCABiEZAAAahGQAAGgQkgEAoEFIBgCABiEZAAAahGQAAGj4/wF0qG0LHepvNwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pyplot.figure(figsize=(12, 8))\n",
    "pyplot.step(bins, theta_cdf_bins_delta, color='green', label='True sampling distribution')\n",
    "pyplot.axvline(x=1, color='black', label='θ (true parameter)')\n",
    "pyplot.legend(loc='upper left')\n",
    "pyplot.ylim(0, 0.1)\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)**\n",
    "\n",
    "For the parametric bootstrap process, the estimated parameter $\\hat{\\theta}^*$ is used on each $k$-th bootstrap sampling $\\{ X_{k1}, X_{k2}, \\dots, X_{kn} \\}$.  But each variable $X_{kj}$ is sampled from  $\\text{Uniform}(0, \\hat{\\theta}^*)$, which is a continuous distribution -- so the probability of obtaining exactly a sample at the boundaries is 0, and $\\mathbb{P}(X_{kj} < \\hat{\\theta}^*) = 1$.  Since the bootstrap functional of each draw, $T(F_n) = \\max(F_n)$ is the largest drawn value in each sample, its values will also always be under $\\hat{\\theta}$, thus the estimated parameter via parametric bootstraping will aywals be under $[\\hat{\\theta}$, and $\\mathbb{P}(\\hat{\\theta}* = \\hat{\\theta}) = 0$.\n",
    "\n",
    "For the nonparametric bootstrap process, the estimated parameter $\\hat{\\theta}^*$ is the maximum value with a point mass in the empirical distribution function $\\hat{F}$.  Each bootstrap resample may or may not include that value when drawing from this sample -- if $\\max\\{X_1, \\dots, X_n\\} \\in \\{X_{k1}, \\dots, X_{kn} \\}$ then the estimated functional for that bootstrap sample will be the estimated parameter $\\hat{\\theta}^*$, otherwise it will necessarily be smaller.\n",
    "\n",
    "Thus, the probability that  $\\mathbb{P}(\\hat{\\theta}* = \\hat{\\theta})$ is the probability that the largest element on the original data is included in a resampling with replacement.  That turns out to be one minus the probability that it never gets included, so,  $\\mathbb{P}(\\hat{\\theta}* = \\hat{\\theta}) = 1 - (1 - (1/n))^n$.  But $\\lim_{n \\rightarrow \\infty } (1 + x/n)^n = e^x$, so the given probability goes to $1 - e^{-1} \\approx 0.632$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 9.6.8**.  Let $T_n = \\overline{X}_n^2$, $\\mu = \\mathbb{E}(X_1)$, $\\alpha_k = \\int |x - \\mu|^k dF(x)$ and $\\hat{\\alpha}_k = n^{-1} \\sum_{i=1}^n |X_i - \\overline{X}_n|^k$.  Show that\n",
    "\n",
    "$$v_\\text{boot} = \\frac{4 \\overline{X}_n^2 \\hat{\\alpha}_2}{n} + \\frac{4 \\overline{X}_n \\hat{\\alpha}_3}{n^2} + \\frac{\\hat{\\alpha}_4}{n^3}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "First, we rewrite the sample mean in terms of an expression containing the central moments.  Let $S_n = n^{-1} \\sum_{i=1}^n (X_i - \\overline{X}_n) = 0$.  Then:\n",
    "\n",
    "$$ \\overline{X}_n = S_n + \\overline{X}_n = \\frac{1}{n} \\sum_{i=1}^n (X_i - \\overline{X}_n) + \\overline{X}_n $$\n",
    "\n",
    "The bootstrap variance, $\\mathbb{V}\\left(\\overline{X}_n^2\\right)$, can be expressed as\n",
    "\n",
    "$$ \\mathbb{V}\\left(\\overline{X}_n^2\\right) = \\mathbb{E}\\left(\\overline{X}_n^4\\right) - \\mathbb{E}\\left(\\overline{X}_n^2\\right)^2 $$\n",
    "\n",
    "Note that $\\overline{X}_n$ is the mean of the distribution, and can be treated as constant when taking expectations.\n",
    "\n",
    "We now have:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{E}\\left(\\overline{X}_n^4\\right) &= \\mathbb{E}\\left( (S_n + \\overline{X}_n)^4 \\right) \\\\\n",
    "&= \\mathbb{E}\\left( S_n^4 + 4 S_n^3 \\overline{X}_n + 6 S_n^2 \\overline{X}_n^2 + 4 S_n \\overline{X}_n^3 + \\overline{X}_n^4 \\right) \\\\\n",
    "&= \\mathbb{E}(S_n^4) + 4 \\overline{X}_n \\mathbb{E}(S_n^3) + 6 \\overline{X}_n^2 \\mathbb{E}(S_n^2) + 4 \\overline{X}_n^3 \\mathbb{E}(S_n) + \\overline{X}_n^4\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Then, computing the moments of $S_n$,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{E}(S_n) &= 0 \\\\\n",
    "\\mathbb{E}(S_n^2) &= \\mathbb{E}\\left( n^{-2} \\left( \\sum_i X_i - \\overline{X}_n \\right)^2 \\right) = \\frac{\\hat{\\alpha}_2}{n} \\\\\n",
    "\\mathbb{E}(S_n^3) &= \\mathbb{E}\\left( n^{-3} \\left( \\sum_i X_i - \\overline{X}_n \\right)^3 \\right) = \\frac{\\hat{\\alpha}_3}{n^2} \\\\\n",
    "\\mathbb{E}(S_n^4) &= \\mathbb{E}\\left( n^{-4} \\left( \\sum_i X_i - \\overline{X}_n \\right)^4 + n^{-2}n^{-4} \\sum_i \\sum_{j \\neq i} (X_i - \\overline{X}_n)^2 (X_j - \\overline{X}_n)^2 \\right) = \\frac{\\hat{\\alpha}_4 + \\hat{\\alpha}_2^2}{n^3}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "and finally\n",
    "\n",
    "$$\n",
    "\\mathbb{E}\\left(\\overline{X}_n^2\\right)^2 = \\mathbb{E}\\left(\\overline{X}_n^2 + S_n\\right)^2 = \\overline{X}_n^4 + 2 \\overline{X}_n^2 \\frac{\\hat{\\alpha}_2}{n} + \\frac{\\hat{\\alpha}_2^2}{n^2}\n",
    "$$\n",
    "\n",
    "Putting everything together,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "v_\\text{boot} &= \\mathbb{E}\\left( (S_n + \\overline{X}_n)^4 \\right) - \\mathbb{E}\\left(\\overline{X}_n^2 + S_n\\right)^2 \\\\\n",
    "&= \\mathbb{E}(S_n^4) + 4 \\overline{X}_n \\mathbb{E}(S_n^3) + 6 \\overline{X}_n^2 \\mathbb{E}(S_n^2) + 4 \\overline{X}_n^3 \\mathbb{E}(S_n) + \\overline{X}_n^4 - \\left( \\overline{X}_n^4 + 2 \\overline{X}_n^2 \\frac{\\hat{\\alpha}_2}{n} + \\frac{\\hat{\\alpha}_2^2}{n^2}\\right) \\\\\n",
    "&= \\frac{\\hat{\\alpha}_4 + \\hat{\\alpha}_2^2}{n^3} + 4 \\overline{X}_n \\frac{\\hat{\\alpha}_3}{n^2} + 6 \\overline{X}_n^2 \\frac{\\hat{\\alpha}_2}{n} + 0 + \\overline{X}_n^4 - \\overline{X}_n^4 - 2 \\overline{X}_n^2 \\frac{\\hat{\\alpha}_2}{n} - \\frac{\\hat{\\alpha}_2^2}{n^2} \\\\\n",
    "&= \\frac{4 \\overline{X}_n^2 \\hat{\\alpha}_2}{n} + \\frac{4 \\overline{X}_n \\hat{\\alpha}_3}{n^2} + \\frac{\\hat{\\alpha}_4}{n^3}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "*Reference and discussion: https://stats.stackexchange.com/q/26082*"
   ]
  }
 ],
 "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
}