09_eval2-resampling-and-CIs/04.2-confidence-intervals_mnist.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "STAT 451: Intro to Machine Learning (Fall 2021)  \n",
    "Instructor: Sebastian Raschka (sraschka@wisc.edu)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# L09: Model Evaluation 2 -- Confidence Intervals and Resampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Author: Sebastian Raschka\n",
      "\n",
      "Last updated: 2021-11-04\n",
      "\n",
      "Python implementation: CPython\n",
      "Python version       : 3.9.6\n",
      "IPython version      : 7.29.0\n",
      "\n",
      "numpy     : 1.21.2\n",
      "mlxtend   : 0.19.0\n",
      "matplotlib: 3.4.3\n",
      "sklearn   : 1.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -a 'Sebastian Raschka' -u -d -v -p numpy,mlxtend,matplotlib,sklearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4.  Confidence Intervals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, we are comparing the confidence intervals from using the normal approximation and using the bootstrap percentile method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.1 MNIST (5k subset) Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mlxtend.data import mnist_data\n",
    "#from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "\n",
    "X, y = mnist_data()\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y,\n",
    "                                                    test_size=0.15, \n",
    "                                                    random_state=12345,\n",
    "                                                    stratify=y)\n",
    "## Too slow\n",
    "#clf = KNeighborsClassifier(n_neighbors=3,\n",
    "#                           weights='uniform', \n",
    "#                           algorithm='kd_tree', \n",
    "#                           leaf_size=30, \n",
    "#                           p=2, \n",
    "#                           metric='minkowski', \n",
    "#                           metric_params=None, \n",
    "#                           n_jobs=1)\n",
    "\n",
    "clf = DecisionTreeClassifier(random_state=123)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.1 Out-of-Bag (OOB) Bootstrap; Bootstrapping Training Sets -- Setup Step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- If you don't tune your model on the training set, you actually don't need a test set for this approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.746178018090561"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "rng = np.random.RandomState(seed=12345)\n",
    "idx = np.arange(y_train.shape[0])\n",
    "\n",
    "bootstrap_train_accuracies = []\n",
    "\n",
    "for i in range(200):\n",
    "    \n",
    "    train_idx = rng.choice(idx, size=idx.shape[0], replace=True)\n",
    "    test_idx = np.setdiff1d(idx, train_idx, assume_unique=False)\n",
    "    \n",
    "    boot_train_X, boot_train_y = X_train[train_idx], y_train[train_idx]\n",
    "    boot_test_X, boot_test_y = X_train[test_idx], y_train[test_idx]\n",
    "    \n",
    "    clf.fit(boot_train_X, boot_train_y)\n",
    "    acc = clf.score(boot_test_X, boot_test_y)\n",
    "    bootstrap_train_accuracies.append(acc)\n",
    "\n",
    "bootstrap_train_mean = np.mean(bootstrap_train_accuracies)\n",
    "bootstrap_train_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.2 Bootstrap 1-sample Confidence Interval (Based on OOB Estimates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7213957368895877 0.7709602992915342\n"
     ]
    }
   ],
   "source": [
    "se = np.sqrt( (1. / (200-1)) * np.sum([(acc - bootstrap_train_mean)**2 \n",
    "                                       for acc in bootstrap_train_accuracies])) \n",
    "ci = 1.96 * se  # 1.96 based on T distribution\n",
    "\n",
    "bootstrap_na_lower = bootstrap_train_mean-ci\n",
    "bootstrap_na_upper = bootstrap_train_mean+ci\n",
    "\n",
    "print(bootstrap_na_lower, bootstrap_na_upper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEYCAYAAABRMYxdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmG0lEQVR4nO3deXRV9b338c8XAgQBlTAJYg3YMCWBAAFlUIKoqCggghPWBFREqy0ur5XHarHjtS2PqEutxVqgOFyEK6I+q7elCCrKFQIFZBSBFBkMAQplhsD3+SOH0wQynCQnObDzfq3FOuf8zh6++8eP7ce999nb3F0AAABBUivWBQAAAEQbAQcAAAQOAQcAAAQOAQcAAAQOAQcAAAROXHWurGnTpp6YmFidqwQAAGe5pUuX7nL3ZtFcZrUGnMTERGVnZ1fnKgEAwFnOzP4R7WVyigoAAAQOAQcAAAQOAQcAAAROtV6DU5zjx49r69atOnLkSKxLAapNfHy8WrdurTp16sS6FAAIpJgHnK1bt6pRo0ZKTEyUmcW6HKDKubt2796trVu3qk2bNrEuBwACKeanqI4cOaImTZoQblBjmJmaNGnCUUsAqEIxDziSCDeocRjzAFC1Igo4Zvaoma02s1Vm9raZxZtZgpnNNbMNodfGVV0sAABAJMoMOGZ2saQfSEp39xRJtSXdIWm8pHnuniRpXujzOScnJ0cpKSlRWdbzzz+vQ4cOlTrN1KlTtX379qis73Q5OTl66623KjRv7969o1wNAACxE+kpqjhJ9c0sTtJ5krZLGiJpWuj7aZKGRr26c0xlA86JEycqtf7SAk5+fn6p837++eeVWjcAAGeTMgOOu2+TNFHSFkk7JO1z979KauHuO0LT7JDUvLj5zWyMmWWbWXZeXl70Ko+i/Px8ZWZmqnPnzho+fHg4pMybN09du3ZVamqqRo8eraNHj5bY/uKLL2r79u3q37+/+vfvrxMnTigrK0spKSlKTU3VpEmTNGvWLGVnZ2vkyJFKS0vT4cOHlZiYqJ/97Gfq27evZs6cqddee009evRQly5ddOutt4ZrycrK0tixY3XllVeqXbt2+vDDD8/YjvHjx+vTTz9VWlqaJk2apKlTp2rEiBG6+eabdd111+nAgQMaMGCAunXrptTUVM2ZMyc8b8OGDSVJCxYsUEZGhoYPH64OHTpo5MiRcvcz1pWRkaFHH31UV111lTp27KglS5Zo2LBhSkpK0lNPPRWe7o033lDPnj2VlpamBx54IBziHnzwQaWnpys5OVkTJkwIT5+YmKgJEyaEa1y3bl1l/3oBVLOZ2d9o0tyvNDP7m1iXghqszJ+Jh66tGSKpjaS9kmaa2d2RrsDdJ0uaLEnp6eln/pfyNDOzv9GspVsjXXyZhndvrRHpl5Q6zfr16/X666+rT58+Gj16tF555RU9/PDDysrK0rx589SuXTvdc889+t3vfqexY8cW2z5u3Dg999xzmj9/vpo2baqlS5dq27ZtWrVqlSRp7969uvDCC/XSSy9p4sSJSk9PD68/Pj5eCxculCTt3r1b999/vyTpqaee0uuvv65HHnlEUsERmo8//lgbN25U//799fXXXys+Pj68nGeffVYTJ04Mh5+pU6dq0aJFWrlypRISEpSfn6/Zs2fr/PPP165du3TFFVdo8ODBZ1zw+ve//12rV69Wq1at1KdPH3322Wfq27fvGf1Wt25dffLJJ3rhhRc0ZMgQLV26VAkJCbrsssv06KOPaufOnZoxY4Y+++wz1alTRw899JDefPNN3XPPPfrlL3+phIQEnThxQgMGDNDKlSvVuXNnSVLTpk21bNkyvfLKK5o4caL+8Ic/lOvvHEBszVq6VV9s3qPL2ySUuf8Fqkok98G5RtJmd8+TJDN7V1JvSblm1tLdd5hZS0k7o1HQ1n8e1heb90RjUZKkK9o2KXOaSy65RH369JEk3X333XrxxRd17bXXqk2bNmrXrp0kKTMzUy+//LL69+9fbPu4ceOKLLNt27batGmTHnnkEQ0aNEjXXXddieu//fbbw+9XrVqlp556Snv37tWBAwc0cODA8He33XabatWqpaSkJLVt21br1q1TWlpaqdt27bXXKiEhQVLB/VeefPJJffLJJ6pVq5a2bdum3NxcXXTRRUXm6dmzp1q3bi1JSktLU05OTrEBZ/DgwZKk1NRUJScnq2XLluFt/+abb7Rw4UItXbpUPXr0kCQdPnxYzZsXHOh75513NHnyZOXn52vHjh1as2ZNOOAMGzZMktS9e3e9++67pW4fAADFiSTgbJF0hZmdJ+mwpAGSsiUdlJQp6dnQ65wSl1AOrRvX1+VtEqKxqPDyynL6EQwzK/a0jKQS20/XuHFjrVixQn/5y1/08ssv65133tEf//jHYqdt0KBB+H1WVpbee+89denSRVOnTtWCBQtKrbMshZf95ptvKi8vT0uXLlWdOnWUmJhY7L1Y6tWrF35fu3btEq/fOTVdrVq1isxTq1Yt5efny92VmZmp//zP/ywy3+bNmzVx4kQtWbJEjRs3VlZWVpE6Ti2rtHUDAFCaMgOOu39hZrMkLZOUL+nvKjjl1FDSO2Z2rwpC0IhoFDQi/ZJqP6S5ZcsWLVq0SL169dLbb7+tvn37qkOHDsrJydHXX3+t7373u5o+fbr69etXYrskNWrUSPv371fTpk21a9cu1a1bV7feeqsuu+wyZWVlFZmmJPv371fLli11/Phxvfnmm7r44ovD382cOVOZmZnavHmzNm3apPbt2xeZt6xl79u3T82bN1edOnU0f/58/eMfUX86fREDBgzQkCFD9Oijj6p58+bas2eP9u/fr3/9619q0KCBLrjgAuXm5urPf/6zMjIyqrQWAEDNEtGjGtx9gqQJpzUfVcHRnHNex44dNW3aND3wwANKSkrSgw8+qPj4eE2ZMkUjRoxQfn6+evToobFjx6pevXrFtkvSmDFjdMMNN6hly5Z6/vnnNWrUKJ08eVKSwkcxTl0sXL9+fS1atOiMWn7+85/r8ssv16WXXqrU1NQigaV9+/bq16+fcnNz9eqrrxa5/kaSOnfurLi4OHXp0kVZWVlq3LjorYlGjhypm2++Wenp6UpLS1OHDh2i2o+n69Spk37xi1/ouuuu08mTJ1WnTh29/PLLuuKKK9S1a1clJyerbdu24dODAABEi0V6yiUa0tPTPTs7u0jb2rVr1bFjx2qr4VyVlZWlm266ScOHD491KYgSxj6C6vbfLwpfZDzjgV6xLgfnADNb6u7pZU8ZubPiUQ0AAADRFPOniSMyU6dOjXUJAACcMziCAwAAAoeAAwAAAoeAAwAAAoeAAwAAAoeAI+nbb7/VHXfcocsuu0ydOnXSjTfeqK+++ko5OTlKSUmRJB07dkyjRo1SamqqunTpUuQOwxkZGWrfvr3S0tKUlpamnTvPfGrFM888o4kTJ1a61tKeGH7K3r179corr1R6XeeCBQsW6Kabbop1GQCAs0yNDzjurltuuUUZGRnauHGj1qxZo1/96lfKzc0tMt1rr70mSfryyy81d+5cPfbYY+Gb+EkFj0FYvny5li9fHn7eUlWobMA59TTvs4G7F+nDcwGPjgCAc0ONDzjz589XnTp1wncjlgoeMHnllVcWmW7NmjUaMKDgxs3NmzfXhRdeqNNvWliWFStW6Oqrr1ZSUlI4MLm7Hn/8caWkpCg1NVUzZswotX38+PH69NNPlZaWpkmTJmn16tXq2bOn0tLS1LlzZ23YsEHjx4/Xxo0blZaWpscff1wLFixQ//79dddddyk1NVWSNHToUHXv3l3JycmaPHlyuMaGDRvqscceU7du3TRgwADl5eWVuk3PPPOMRo8erYyMDLVt21Yvvvhi+LvnnntOKSkpSklJ0fPPPy+pIKB17NhRDz30kLp166ZPP/1UHTp00H333aeUlBSNHDlSf/vb39SnTx8lJSVp8eLFkqTFixerd+/e6tq1q3r37q3169eXWldOTo6uvPJKdevWTd26ddPnn38uqeCIz1VXXaVbbrlFnTp10tixY8Mhq6Rtz8jI0JNPPql+/frphRde0Lx589S1a1elpqZq9OjROnr0qJYsWaLOnTvryJEjOnjwoJKTk8NPkgcAxIC7V9uf7t27++nWrFlT5PPsDbM9689ZnvXnrDOm/Y8F/+FZf87y11a+VqR97e614XnW7l5b5Lu8Q3lnLKewF154wceNG1fsd5s3b/bk5GR3d//973/vw4cP9+PHj/umTZv8ggsu8FmzZrm7e79+/TwlJcW7dOniP/vZz/zkyZNnLGvChAneuXNnP3TokOfl5Xnr1q1927ZtPmvWLL/mmms8Pz/fv/32W7/kkkt8+/btJbbPnz/fBw0aFF7uww8/7G+88Ya7ux89etQPHTpUpG539/nz5/t5553nmzZtCrft3r3b3d0PHTrkycnJvmvXLnd3lxRe3k9/+lP//ve/X2r/TZgwwXv16uVHjhzxvLw8T0hI8GPHjnl2dranpKT4gQMHfP/+/d6pUydftmyZb9682c3MFy1aFO7j2rVr+8qVK/3EiRPerVs3HzVqlJ88edLfe+89HzJkiLu779u3z48fP+7u7nPnzvVhw4aFt61wf5xy8OBBP3z4sLu7f/XVV35q7M2fP9/r1avnGzdu9Pz8fL/mmmt85syZpW57v379/MEHH3R398OHD3vr1q19/fr17u7+ve99zydNmuTu7j/+8Y/9scce84ceesh/9atfldpv7meOfSAobnv1c7/0iQ/9tlc/j3UpOEdIyvYoZ46z7kZ/2w9sV3Zu8UdGVuat1PaD29WqYasi7fuP7Q/Ps/9Y0YdNHj1xNCp1jR49WmvXrlV6erouvfRS9e7dW3FxBd136qGY+/fv16233qrp06frnnvuOWMZQ4YMUf369VW/fn31799fixcv1sKFC3XnnXeqdu3aatGihfr166clS5aU2H7++ecXWWavXr30y1/+Ulu3btWwYcOUlJRUbP09e/ZUmzZtwp9ffPFFzZ49W5L0zTffaMOGDWrSpIlq1aql22+/XZJ09913a9iwYWX2zaBBg1SvXj3Vq1dPzZs3V25urhYuXKhbbrkl/DTzYcOG6dNPP9XgwYN16aWX6oorrgjP36ZNm/CRpeTkZA0YMEBmptTUVOXk5EgqeFBoZmamNmzYIDPT8ePHS63p+PHjevjhh7V8+XLVrl1bX331VZG+aNu2rSTpzjvv1MKFCzV8+PBSt/1U+/r169WmTRu1a9dOkpSZmamXX35Z48aN009+8hP16NFD8fHxRY5kAQCq31l3iqpVw1ZKb5Gu9BZnPpKic7POSm+RrjYXtCnS3qhuo/A8jeo2KvJdvdr1Sl1fcnKyli5dWmZdcXFxmjRpkpYvX645c+Zo79694TBx6onfjRo10l133RU+rXI6Mzvjs5fwLLCS2k9311136f3331f9+vU1cOBAffTRR8VOdypoSAWnaf72t79p0aJFWrFihbp27aojR45EVHNx6tX7dx/Xrl1b+fn5pdZfuJbT569Vq1b4c61atcLXvDz99NPq37+/Vq1apQ8++KDEek+ZNGmSWrRooRUrVig7O1vHjh0rcZtK2sbC7adqLm279uzZowMHDmj//v1l1gcAqFpnXcAZ+t2hmnL9FE25fsoZ3/2232815fopui/1viLtHRI6hOfpkFD0CdlN6zctdX1XX321jh49Gr4mRpKWLFmijz/+uMh0hw4d0sGDByVJc+fOVVxcnDp16qT8/Hzt2rVLUsFRgw8//DD8y6vTzZkzR0eOHNHu3bu1YMEC9ejRQ1dddZVmzJihEydOKC8vT5988ol69uxZYnujRo2KPGF806ZNatu2rX7wgx9o8ODBWrly5RnTnG7fvn1q3LixzjvvPK1bt07/+7//G/7u5MmTmjVrliTprbfeUt++fSVJL730kl566aVS+7Kwq666Su+9916432bPnn3GdU3lsW/fvnCQjOSxFfv27VPLli1Vq1YtTZ8+vcjF1YsXL9bmzZt18uRJzZgxI7yNJW17YR06dFBOTo6+/vprSdL06dPVr18/SQVPk//5z3+ukSNH6oknnqjwtgIAKu+sO0VV3cxMs2fP1rhx4/Tss88qPj5eiYmJ4YtiT9m5c6cGDhyoWrVq6eKLL9b06dMlSUePHtXAgQN1/PhxnThxQtdcc43uv//+YtfVs2dPDRo0SFu2bNHTTz+tVq1a6ZZbbtGiRYvUpUsXmZl+85vf6KKLLiqxvUmTJoqLi1OXLl2UlZWlI0eO6I033lCdOnV00UUX6Sc/+YkSEhLUp08fpaSk6IYbbtCgQYOK1HH99dfr1VdfVefOndW+ffsip4saNGig1atXq3v37rrgggvCFzevW7dOffr0ibhfu3XrpqysLPXs2VOSdN9996lr167hU07l9aMf/UiZmZl67rnndPXVV5c5/UMPPaRbb71VM2fOVP/+/YscNerVq5fGjx+vL7/8MnzBsVTythcWHx+vKVOmaMSIEcrPz1ePHj00duxY/elPf1JcXJzuuusunThxQr1799ZHH32kq6++WmlpaVq+fHmFthsAUDEW6amQaEhPT/fTf3m0du1adezYsdpqQOkaNmyoAwcOnNF+00036d1331XdunVjUFX0LFiwQBMnTtSHH354xnclbXtVYewjqG7//SJ9sXmPLm+ToBkP9Ip1OTgHmNlSdz/z2pRKqPFHcBCZ4gIBAABnKwIOiqjOIxixkJGRoYyMjGK/C/q2A0BNclZcZFydp8mAswFjHgCqVswDTnx8vHbv3s0OHzWGu2v37t2Kj4+PdSkAEFgxP0XVunVrbd26tcxHAgBBEh8fr9atW8e6DAAIrDIDjpm1l1T497JtJf1E0p9C7YmSciTd5u7/LG8BderUKXKHXQAAgMoq8xSVu6939zR3T5PUXdIhSbMljZc0z92TJM0LfQYAAIi58l6DM0DSRnf/h6QhkqaF2qdJGhrFugAAACqsvAHnDklvh963cPcdkhR6bV7cDGY2xsyyzSyb62wAAEB1iDjgmFldSYMlzSzPCtx9srunu3t6s2bNylsfAABAuZXnCM4Nkpa5e27oc66ZtZSk0OvOaBcHAABQEeUJOHfq36enJOl9SZmh95mS5kSrKAAAgMqIKOCY2XmSrpX0bqHmZyVda2YbQt89G/3yAAAAyi+iG/25+yFJTU5r262CX1UBAACcVWL+qAYAAIBoI+AAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAiYt1AQCA6vPB+qpfx+7DVb8OoCwcwQEAAIETUcAxswvNbJaZrTOztWbWy8wSzGyumW0IvTau6mIBAAAiEekRnBck/Y+7d5DURdJaSeMlzXP3JEnzQp8BAABirsyAY2bnS7pK0uuS5O7H3H2vpCGSpoUmmyZpaNWUCAAAUD6RHMFpKylP0hQz+7uZ/cHMGkhq4e47JCn02rwK6wQAAIhYJAEnTlI3Sb9z966SDqocp6PMbIyZZZtZdl5eXgXLBAAAiFwkAWerpK3u/kXo8ywVBJ5cM2spSaHXncXN7O6T3T3d3dObNWsWjZoBAABKVWbAcfdvJX1jZu1DTQMkrZH0vqTMUFumpDlVUiEAAEA5RXqjv0ckvWlmdSVtkjRKBeHoHTO7V9IWSSOqpkQAwLlo9+HqubGgJN3cvuxpULNEFHDcfbmk9GK+GhDVagAAAKKAOxkDAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAIeAAAIDAiYtkIjPLkbRf0glJ+e6ebmYJkmZISpSUI+k2d/9n1ZQJAAAQufIcwenv7mnunh76PF7SPHdPkjQv9BkAACDmKnOKaoikaaH30yQNrXQ1AAAAURBpwHFJfzWzpWY2JtTWwt13SFLotXlxM5rZGDPLNrPsvLy8ylcMAABQhoiuwZHUx923m1lzSXPNbF2kK3D3yZImS1J6erpXoEYAAIByiegIjrtvD73ulDRbUk9JuWbWUpJCrzurqkgAAIDyKDPgmFkDM2t06r2k6yStkvS+pMzQZJmS5lRVkQAAAOURySmqFpJmm9mp6d9y9/8xsyWS3jGzeyVtkTSi6soEAACIXJkBx903SepSTPtuSQOqoigAAIDK4E7GAAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcCIOOGZW28z+bmYfhj4nmNlcM9sQem1cdWUCAABErjxHcH4oaW2hz+MlzXP3JEnzQp8BAABiLqKAY2atJQ2S9IdCzUMkTQu9nyZpaFQrAwAAqKBIj+A8L+lHkk4Wamvh7jskKfTavLgZzWyMmWWbWXZeXl5lagUAAIhImQHHzG6StNPdl1ZkBe4+2d3T3T29WbNmFVkEAABAucRFME0fSYPN7EZJ8ZLON7M3JOWaWUt332FmLSXtrMpCAQAAIlXmERx3/z/u3trdEyXdIekjd79b0vuSMkOTZUqaU2VVAgAAlENl7oPzrKRrzWyDpGtDnwEAAGIuklNUYe6+QNKC0PvdkgZEvyQAAIDK4U7GAAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcOJiXQAAAJX1wfrqW9fN7atvXag4juAAAIDAIeAAAIDAKTPgmFm8mS02sxVmttrMfhpqTzCzuWa2IfTauOrLBQAAKFskR3COSrra3btISpN0vZldIWm8pHnuniRpXugzAABAzJUZcLzAgdDHOqE/LmmIpGmh9mmShlZFgQAAAOUV0TU4ZlbbzJZL2ilprrt/IamFu++QpNBr8xLmHWNm2WaWnZeXF6WyAQAAShZRwHH3E+6eJqm1pJ5mlhLpCtx9srunu3t6s2bNKlgmAABA5Mr1Kyp33ytpgaTrJeWaWUtJCr3ujHZxAAAAFRHJr6iamdmFoff1JV0jaZ2k9yVlhibLlDSnimoEAAAol0juZNxS0jQzq62CQPSOu39oZoskvWNm90raImlEFdYJAAAQsTIDjruvlNS1mPbdkgZURVEAAACVwZ2MAQBA4BBwAABA4BBwAABA4BBwAABA4BBwAABA4BBwAABA4BBwAABA4BBwAABA4BBwAABA4BBwAABA4BBwAABA4BBwAABA4ETyNHEAQBX6YH2sKwCChyM4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcAg4AAAgcMoMOGZ2iZnNN7O1ZrbazH4Yak8ws7lmtiH02rjqywUAAChbJEdw8iU95u4dJV0h6ftm1knSeEnz3D1J0rzQZwAAgJgrM+C4+w53XxZ6v1/SWkkXSxoiaVposmmShlZRjQAAAOVSrmtwzCxRUldJX0hq4e47pIIQJKl5CfOMMbNsM8vOy8urZLkAAABlizjgmFlDSf8taZy7/yvS+dx9srunu3t6s2bNKlIjAABAuUQUcMysjgrCzZvu/m6oOdfMWoa+bylpZ9WUCAAAUD6R/IrKJL0uaa27P1foq/clZYbeZ0qaE/3yAAAAyi8ugmn6SPqepC/NbHmo7UlJz0p6x8zulbRF0ogqqRAAAKCcygw47r5QkpXw9YDolgMAAFB53MkYAAAEDgEHAAAEDgEHAAAEDgEHAAAEDgEHAAAEDgEHAAAEDgEHAAAEDgEHAAAEDgEHAAAEDgEHAAAEDgEHAAAEDgEHAAAETiRPEwcAlMOvF/9a6/asU4eEDnqi5xOxLgcxwBiIPQIOAETZuj3rlJ2bHesyEEOMgdgj4ABAlHVI6FDkFTUPYyD2CDgAEGWckgBjIPa4yBgAAAQOAQcAAAQOp6gAIMp2Hd6loyeOql7tempav2msy0EMMAZij4ADAFH2+MePKzs3W+kt0jXl+imxLgcxwBiIPU5RAQCAwCnzCI6Z/VHSTZJ2untKqC1B0gxJiZJyJN3m7v+sujIB4Nxxb+q9GvLdIZyaqMEYA7EXySmqqZJekvSnQm3jJc1z92fNbHzoM7+JAwBJfS/uG+sSEGOMgdgr8xSVu38iac9pzUMkTQu9nyZpaHTLAgAAqLiKXmTcwt13SJK77zCz5iVNaGZjJI2RpO985zsVXB0AVK8P1se6ApytqnNs3Ny++tYVNFX+Kyp3nyxpsiSlp6d7Va8PAGJt6c739M+j29W4Xit1bz401uUgBhgDsVfRX1HlmllLSQq97oxeSQBwbluWN0cfbfudluXNiXUpiBHGQOxVNOC8Lykz9D5TEn+DAADgrBHJz8TflpQhqamZbZU0QdKzkt4xs3slbZE0oiqLBIBzyf3J3NitpmMMxF6ZAcfd7yzhqwFRrgUAACAquJMxAAAIHAIOAAAIHB62CQBR9vaGx/XN/pW6pFFn3Zn021iXgxhgDMQeR3AAIMoOHNulvce268CxXbEuBTHCGIg9juAAQJS1u7CPGtdrpWb128S6FMQIYyD2CDgAEGX9Lr4v1iUgxhgDsccpKgAAEDgEHAAAEDicogKAKNt+cJ2O5O9XfFwjtWrQIdblIAYYA7FHwAGAKPt/Ob/W5v3ZatMonVv211CMgdjjFBUAAAgcjuAAQJQNSnwifHoCNRNjIPYIOAAQZVxzAcZA7HGKCgAABA4BBwAABA6nqACcMz5YH+sKIvPxtj8o7/BmNavfhjva1lCMgdjjCA4ARNlXez/Tsl3v66u9n8W6FMQIYyD2CDgAEGUN6zbVhXVbqWHdprEuBTHCGIg9TlEBQJTdmfTbWJeAGGMMxB5HcAAAQOAQcAAAQOBUKuCY2fVmtt7Mvjaz8dEqCgAAoDIqfA2OmdWW9LKkayVtlbTEzN539zXRKg4AzkWvrR7FgxZrOMZA7FXmCE5PSV+7+yZ3PybpvyQNiU5ZAAAAFWfuXrEZzYZLut7d7wt9/p6ky9394dOmGyNpTOhjiqRVFS83EJpK2hXrImKMPqAPJPpAog8k+kCiDySpvbtH9cmklfmZuBXTdkZacvfJkiZLkpllu3t6JdZ5zqMP6AOJPpDoA4k+kOgDiT6QCvog2suszCmqrZIuKfS5taTtlSsHAACg8ioTcJZISjKzNmZWV9Idkt6PTlkAAAAVV+FTVO6eb2YPS/qLpNqS/ujuq8uYbXJF1xcg9AF9INEHEn0g0QcSfSDRB1IV9EGFLzIGAAA4W3EnYwAAEDgEHAAAEDgVDjiRPKbBzDLMbLmZrTazj8ua18wSzGyumW0IvTauaH3VoaJ9YGaXmNl8M1sbav9hoemfMbNtoXmWm9mN1bU9FVHJcZBjZl+Gvssu1F5TxkH7Qn/Py83sX2Y2LvRdoMaBmT1eaFtWmdkJM0sobd6gjYOS+qAm7Q/KGAc1Yn9QyjioSfuDC8zsAzNbERrzo8qat0LjwN3L/UcFFxVvlNRWUl1JKyR1Om2aCyWtkfSd0OfmZc0r6TeSxofej5f064rUVx1/KtkHLSV1C71vJOmrQn3wjKT/iPX2VXUfhN7nSGpazHJrxDgoZjnfSro0iOPgtOlvlvRRWfMGbRyU0gc1Zn9QUh+EPteI/UFpfXDacgK7P5D05Km/R0nNJO0JTRvV/UFFj+BE8piGuyS96+5bJMndd0Yw7xBJ00Lvp0kaWsH6qkOF+8Ddd7j7stD7/ZLWSrq42iqPnsqMg9LUiHFwmgGSNrr7P6q02qpR3se23Cnp7QjmDdo4KCzcBzVsf1BY4XFQmhoxDk4T9P2BS2pkZiapoQoCTn4Z85Z7HFQ04Fws6ZtCn7fqzH+Q7SQ1NrMFZrbUzO6JYN4W7r5DKvhHL6l5BeurDpXpgzAzS5TUVdIXhZofNrOVZvbHs/xwbGX7wCX9NdQ+plB7jRsHKriP1Ok7uiCNA0mSmZ0n6XpJ/x3BvEEbB5KK7YPC3yUq2PsDSSX2QU3ZH0gqfRwo+PuDlyR1VMHNgb+U9EN3P1nGvOUeBxUNOJE8piFOUndJgyQNlPS0mbWLcN5zQWX6oGABZg1VMLjHufu/Qs2/k3SZpDRJOyT93+iWHVWV7YM+7t5N0g2Svm9mV1VZpVUnGuOgrqTBkmYWmido4+CUmyV95u57KjDv2awyfVCwgJqxPziluD6oKfuDU0oaBzVhfzBQ0nJJrVSwTS+Z2fkRzhuxigacSB7TsFXS/7j7QXffJekTSV3KmDfXzFpKUug1ktMZsVKZPpCZ1VHBzuxNd3/31AzunuvuJ0Jp9jUVHLI7W1WqD9x9e+h1p6TZ+ve21phxEHKDpGXunnuqIYDj4JTT/8+0Ju0PTjnj/85r0P7glDP6oAbtD04p7iiNVDP2B6NUcNre3f1rSZsldShj3vKPgwpeRBQnaZOkNvr3hUDJp03TUdK80LTnqeAp4imlzSvptyp6EdFvKlJfdfypZB+YpD9Jer6Y5bYs9P5RSf8V622toj5oIKlRaJoGkj5XwdPpa8w4KPT9f0kaFeRxEJruAhWca28QybxBGwel9EGN2R+U0gc1Zn9QUh8U+i7w+wMVHJF6JvS+haRtKniielT3B5XZiBtVcLX/Rkk/DrWNlTS20DSPq+DXI6tUcNi1xHlD7U1U8B+CDaHXhFj/ZVVFH0jqq4LDbitVcJhuuaQbQ99NV8E5yZUqeLZXy+rcpmrsg7ahwbtC0uqaOA5C7edJ2i3pgtOWGcRxkKVidsw1bH9wRh/UwP1BcX1Q0/YHJf1bqBH7AxWcmvpraJtWSbq7tHkrOg54VAMAAAgc7mQMAAACh4ADAAACh4ADAAACh4ADAAACh4ADAAACh4ADAAACh4ADAAAC5/8D7sI553wyaXEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "ax.vlines( bootstrap_train_mean, [0], 80, lw=2.5, linestyle='-', label='bootstrap train mean')\n",
    "\n",
    "ax.vlines(bootstrap_na_upper, [0], 15, lw=2.5, linestyle='dotted', \n",
    "          label='CI95 bootstrap, normal approx.', color='C2')\n",
    "ax.vlines(bootstrap_na_lower, [0], 15, lw=2.5, linestyle='dotted', color='C2')\n",
    "\n",
    "ax.hist(bootstrap_train_accuracies, bins=7,\n",
    "        color='#0080ff', edgecolor=\"none\", \n",
    "        alpha=0.3)\n",
    "plt.legend(loc='upper left')\n",
    "\n",
    "plt.xlim([0.6, 0.8])\n",
    "plt.tight_layout()\n",
    "#plt.savefig('figures/bootstrap-histo-1.svg')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.3 Bootstrap Percentile Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7204597701149426 0.7696229551148738\n"
     ]
    }
   ],
   "source": [
    "bootstrap_percentile_lower = np.percentile(bootstrap_train_accuracies, 2.5)\n",
    "bootstrap_percentile_upper = np.percentile(bootstrap_train_accuracies, 97.5)\n",
    "\n",
    "print(bootstrap_percentile_lower, bootstrap_percentile_upper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEYCAYAAABRMYxdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvNUlEQVR4nO3deXRV1d3/8fc3JBCQMYwBlMEiUyZCiDIoAVRQlJmqYCVQB6QOuKwt9VGhtvpQyyPKUktxQIrYB/EnonbZp4hERaMQbEBGEUgBQQhQMMwk7N8fudwmIcPNeMnJ57VW1r1nn7PP+e67c8KXs88925xziIiIiHhJSLADEBEREaloSnBERETEc5TgiIiIiOcowRERERHPUYIjIiIinhNalQdr1qyZa9++fVUeUkRERDxs7dq1B51zzQuWV2mC0759e9LS0qrykCIiIuJhZvavwso1RCUiIiKeowRHREREPEcJjoiIiHhOld6DU5izZ8+yZ88eTp06FexQRKpMeHg4bdu2JSwsLNihiIh4UtATnD179tCgQQPat2+PmQU7HJFK55zj0KFD7Nmzhw4dOgQ7HBERTwr6ENWpU6do2rSpkhupMcyMpk2b6qqliEglCnqCAyi5kRpHv/MiIpUroATHzB4ys41mtsHM/mpm4WYWYWbLzWyb77VJZQcrIiIiEogSExwzawM8ACQ456KAWsCtwDRghXOuE7DCt1ztZGRkEBUVVSH7eu655zhx4kSx27z++uvs3bu3Qo5XUEZGBm+++WaZ6vbp06eCoxEREQmeQIeoQoG6ZhYK1AP2AsOBBb71C4ARFR5dNVPeBCcnJ6dcxy8uwcnOzi627hdffFGuY4uIiFxMSkxwnHPfA7OAXcA+4Khz7h9AS+fcPt82+4AWhdU3s7vNLM3M0jIzMysu8gqUnZ3NhAkTiImJYcyYMf4kZcWKFfTo0YPo6GgmTZrE6dOniyyfM2cOe/fuZcCAAQwYMICcnBySk5OJiooiOjqa2bNn8/bbb5OWlsb48eOJi4vj5MmTtG/fnieffJJ+/fqxZMkSXn75ZXr16kVsbCyjR4/2x5KcnMzkyZO5+uqrueKKK/jggw8uaMe0adP47LPPiIuLY/bs2bz++uuMHTuWm2++meuvv55jx44xaNAg4uPjiY6OZtmyZf669evXByAlJYWkpCTGjBlDly5dGD9+PM65C46VlJTEQw89xDXXXEPXrl1Zs2YNo0aNolOnTjz22GP+7d544w0SExOJi4vjnnvu8Sdx9957LwkJCXTv3p3p06f7t2/fvj3Tp0/3x7hly5bydq+IVLAlabuZvfxblqTtDnYoIkUq8WvivntrhgMdgCPAEjO7PdADOOfmAfMAEhISLvyXsoAlabt5e+2eQHdfojE92zI24dJit9m6dSuvvvoqffv2ZdKkSbz00kvcd999JCcns2LFCq644gruuOMO/vSnPzF58uRCy6dOncqzzz7LypUradasGWvXruX7779nw4YNABw5coTGjRvzwgsvMGvWLBISEvzHDw8PZ9WqVQAcOnSIu+66C4DHHnuMV199lfvvvx/IvULzySefsH37dgYMGMB3331HeHi4fz8zZ85k1qxZ/uTn9ddfJzU1lfXr1xMREUF2djZLly6lYcOGHDx4kKuuuophw4ZdcMPrP//5TzZu3Ejr1q3p27cvn3/+Of369bvgc6tduzaffvopzz//PMOHD2ft2rVERERw+eWX89BDD3HgwAEWL17M559/TlhYGFOmTGHRokXccccdPPXUU0RERJCTk8OgQYNYv349MTExADRr1oyvv/6al156iVmzZvHKK6+Uqs9FpHK9vXYPX+08zJUdIkr8+yoSLIE8B+daYKdzLhPAzN4B+gD7zSzSObfPzCKBAxUR0J5/n+SrnYcrYlcAXNWxaYnbXHrppfTt2xeA22+/nTlz5nDdddfRoUMHrrjiCgAmTJjAiy++yIABAwotnzp1ar59duzYkR07dnD//fczdOhQrr/++iKPf8stt/jfb9iwgccee4wjR45w7NgxBg8e7F/305/+lJCQEDp16kTHjh3ZsmULcXFxxbbtuuuuIyIiAsh9/sqjjz7Kp59+SkhICN9//z379++nVatW+eokJibStm1bAOLi4sjIyCg0wRk2bBgA0dHRdO/encjISH/bd+/ezapVq1i7di29evUC4OTJk7RokXuh76233mLevHlkZ2ezb98+Nm3a5E9wRo0aBUDPnj155513im2fiIhIYQJJcHYBV5lZPeAkMAhIA44DE4CZvtdlRe6hFNo2qcuVHSIqYlf+/ZWk4BUMMyt0WAYosrygJk2asG7dOv7v//6PF198kbfeeovXXnut0G0vueQS//vk5GTeffddYmNjef3110lJSSk2zpLk3feiRYvIzMxk7dq1hIWF0b59+0KfxVKnTh3/+1q1ahV5/8757UJCQvLVCQkJITs7G+ccEyZM4L//+7/z1du5cyezZs1izZo1NGnShOTk5HxxnN9XcccWEREpTokJjnPuKzN7G/gayAb+Se6QU33gLTP7OblJ0NiKCGhswqVVfslz165dpKam0rt3b/7617/Sr18/unTpQkZGBt999x0/+clPWLhwIf379y+yHKBBgwZkZWXRrFkzDh48SO3atRk9ejSXX345ycnJ+bYpSlZWFpGRkZw9e5ZFixbRpk0b/7olS5YwYcIEdu7cyY4dO+jcuXO+uiXt++jRo7Ro0YKwsDBWrlzJv/5V6AzzFWbQoEEMHz6chx56iBYtWnD48GGysrL48ccfueSSS2jUqBH79+/nww8/JCkpqVJjERGRmiWgqRqcc9OB6QWKT5N7Nafa69q1KwsWLOCee+6hU6dO3HvvvYSHhzN//nzGjh1LdnY2vXr1YvLkydSpU6fQcoC7776bG264gcjISJ577jkmTpzIuXPnAPxXMc7fLFy3bl1SU1MviOV3v/sdV155Je3atSM6OjpfwtK5c2f69+/P/v37mTt3br77bwBiYmIIDQ0lNjaW5ORkmjTJ/2ii8ePHc/PNN5OQkEBcXBxdunSp0M+xoG7duvH73/+e66+/nnPnzhEWFsaLL77IVVddRY8ePejevTsdO3b0Dw+KiIhUFAt0yKUiJCQkuLS0tHxlmzdvpmvXrlUWQ3WVnJzMTTfdxJgxY4IdilQQ/e5LdXXLn1P9Nxkvvqd3sMORGs7M1jrnEgqWXxRTNYiIiIhUpKDPJi6Bef3114MdgoiISLWhKzgiIiLiOUpwRERExHOU4IiIiIjnKMERERERz1GCA/zwww/ceuutXH755XTr1o0bb7yRb7/9loyMDKKiogA4c+YMEydOJDo6mtjY2HxPGE5KSqJz587ExcURFxfHgQMXzloxY8YMZs2aVe5Yi5sx/LwjR47w0ksvlftY1UFKSgo33XRTsMMQEZGLTI1PcJxzjBw5kqSkJLZv386mTZt4+umn2b9/f77tXn75ZQC++eYbli9fzsMPP+x/iB/kToOQnp5Oenq6f76lylDeBOf8bN4XA+dcvs+wOtDUESIi1UONT3BWrlxJWFiY/2nEkDvB5NVXX51vu02bNjFoUO6Dm1u0aEHjxo0p+NDCkqxbt46BAwfSqVMnf8LknOORRx4hKiqK6OhoFi9eXGz5tGnT+Oyzz4iLi2P27Nls3LiRxMRE4uLiiImJYdu2bUybNo3t27cTFxfHI488QkpKCgMGDGDcuHFER0cDMGLECHr27En37t2ZN2+eP8b69evz8MMPEx8fz6BBg8jMzCy2TTNmzGDSpEkkJSXRsWNH5syZ41/37LPPEhUVRVRUFM899xyQm6B17dqVKVOmEB8fz2effUaXLl248847iYqKYvz48Xz00Uf07duXTp06sXr1agBWr15Nnz596NGjB3369GHr1q3FxpWRkcHVV19NfHw88fHxfPHFF0DuFZ9rrrmGkSNH0q1bNyZPnuxPsopqe1JSEo8++ij9+/fn+eefZ8WKFfTo0YPo6GgmTZrE6dOnWbNmDTExMZw6dYrjx4/TvXt3/0zyIiISBM65Kvvp2bOnK2jTpk35lpduW+qSP0x2yR8mX7DtL1N+6ZI/THYvr385X/nmQ5v9dTYf2pxvXeaJzAv2k9fzzz/vpk6dWui6nTt3uu7duzvnnPvzn//sxowZ486ePet27NjhGjVq5N5++23nnHP9+/d3UVFRLjY21j355JPu3LlzF+xr+vTpLiYmxp04ccJlZma6tm3buu+//969/fbb7tprr3XZ2dnuhx9+cJdeeqnbu3dvkeUrV650Q4cO9e/3vvvuc2+88YZzzrnTp0+7EydO5IvbOedWrlzp6tWr53bs2OEvO3TokHPOuRMnTrju3bu7gwcPOuecA/z7++1vf+t+8YtfFPv5TZ8+3fXu3dudOnXKZWZmuoiICHfmzBmXlpbmoqKi3LFjx1xWVpbr1q2b+/rrr93OnTudmbnU1FT/Z1yrVi23fv16l5OT4+Lj493EiRPduXPn3LvvvuuGDx/unHPu6NGj7uzZs84555YvX+5GjRrlb1vez+O848ePu5MnTzrnnPv222/d+d+9lStXujp16rjt27e77Oxsd+2117olS5YU2/b+/fu7e++91znn3MmTJ13btm3d1q1bnXPO/exnP3OzZ892zjn3X//1X+7hhx92U6ZMcU8//XSxn5tzF/7ui1QXP537hWv36w/cT+d+EexQRByQ5grJOS66B/3tPbaXtP2FXxlZn7mevcf30rp+63zlWWey/HWyzuSfbPJ0zukKiWvSpEls3ryZhIQE2rVrR58+fQgNzf34zk+KmZWVxejRo1m4cCF33HHHBfsYPnw4devWpW7dugwYMIDVq1ezatUqbrvtNmrVqkXLli3p378/a9asKbK8YcOG+fbZu3dvnnrqKfbs2cOoUaPo1KlTofEnJibSoUMH//KcOXNYunQpALt372bbtm00bdqUkJAQbrnlFgBuv/12Ro0aVeJnM3ToUOrUqUOdOnVo0aIF+/fvZ9WqVYwcOdI/m/moUaP47LPPGDZsGO3ateOqq67y1+/QoYP/ylL37t0ZNGgQZkZ0dDQZGRlA7kShEyZMYNu2bZgZZ8+eLTams2fPct9995Genk6tWrX49ttv830WHTt2BOC2225j1apVjBkzpti2ny/funUrHTp04IorrgBgwoQJvPjii0ydOpUnnniCXr16ER4enu9KloiIVL2Lboiqdf3WJLRMIKHlBdNKENM8hoSWCXRo1CFfeYPaDfx1GtRukG9dnVp1ij1e9+7dWbt2bYlxhYaGMnv2bNLT01m2bBlHjhzxJxPnZ/xu0KAB48aN8w+rFGRmFyy7IuYCK6q8oHHjxvHee+9Rt25dBg8ezMcff1zoducTDcgdpvnoo49ITU1l3bp19OjRg1OnTgUUc2Hq1PnPZ1yrVi2ys7OLjT9vLAXrh4SE+JdDQkL897w8/vjjDBgwgA0bNvD+++8XGe95s2fPpmXLlqxbt460tDTOnDlTZJuKamPe8vMxF9euw4cPc+zYMbKyskqMT0REKtdFl+CM+MkI5g+Zz/wh8y9Y98f+f2T+kPncGX1nvvIuEV38dbpE5J8hu1ndZsUeb+DAgZw+fdp/TwzAmjVr+OSTT/Jtd+LECY4fPw7A8uXLCQ0NpVu3bmRnZ3Pw4EEg96rBBx984P/mVUHLli3j1KlTHDp0iJSUFHr16sU111zD4sWLycnJITMzk08//ZTExMQiyxs0aJBvhvEdO3bQsWNHHnjgAYYNG8b69esv2Kago0eP0qRJE+rVq8eWLVv48ssv/evOnTvH22+/DcCbb75Jv379AHjhhRd44YUXiv0s87rmmmt49913/Z/b0qVLL7ivqTSOHj3qTyQDmbbi6NGjREZGEhISwsKFC/PdXL169Wp27tzJuXPnWLx4sb+NRbU9ry5dupCRkcF3330HwMKFC+nfvz+QO5v87373O8aPH8+vf/3rMrdVRETK76IboqpqZsbSpUuZOnUqM2fOJDw8nPbt2/tvij3vwIEDDB48mJCQENq0acPChQsBOH36NIMHD+bs2bPk5ORw7bXXctdddxV6rMTERIYOHcquXbt4/PHHad26NSNHjiQ1NZXY2FjMjGeeeYZWrVoVWd60aVNCQ0OJjY0lOTmZU6dO8cYbbxAWFkarVq144okniIiIoG/fvkRFRXHDDTcwdOjQfHEMGTKEuXPnEhMTQ+fOnfMNF11yySVs3LiRnj170qhRI//NzVu2bKFv374Bf67x8fEkJyeTmJgIwJ133kmPHj38Q06l9atf/YoJEybw7LPPMnDgwBK3nzJlCqNHj2bJkiUMGDAg31Wj3r17M23aNL755hv/DcdQdNvzCg8PZ/78+YwdO5bs7Gx69erF5MmT+ctf/kJoaCjjxo0jJyeHPn368PHHHzNw4EDi4uJIT08vU7tFRKRsLNChkIqQkJDgCn7zaPPmzXTt2rXKYpDi1a9fn2PHjl1QftNNN/HOO+9Qu3btIERVcVJSUpg1axYffPDBBeuKantl0e++VFe3/DmVr3Ye5soOESy+p3eww5EazszWOucuuK+lxl/BkcAUlhCIiIhcrJTgSD5VeQUjGJKSkkhKSip0ndfbLiJSk1x0NxmLiIiIlJcSHBEREfEcJTgiIiLiOSUmOGbW2czS8/z8aGZTzSzCzJab2Tbfa5OqCFhERESkJCUmOM65rc65OOdcHNATOAEsBaYBK5xznYAVvuVq6YcffuDWW2/l8ssvp1u3btx44418++23ZGRk+B/ad+bMGSZOnEh0dDSxsbGkpKT46yclJdG5c2fi4uKIi4vjwIEDFxxjxowZzJo1q9yxlnc2ca8q+LmkpaXxwAMPALkPBrzvvvuCFZqIiARBaYeoBgHbnXP/AoYDC3zlC4ARFRhXlXHOMXLkSJKSkti+fTubNm3i6aefZv/+/fm2O/+k42+++Ybly5fz8MMP+2ehhtz5qNLT00lPT6dFixaVFm95E5y8T/QNNudcvs+wPAp+LgkJCZoPSkSkBittgnMr8Fff+5bOuX0AvtdC/1U3s7vNLM3M0jIzM8seaSVZuXIlYWFhTJ482V8WFxd3wbQCmzZtYtCgQQC0aNGCxo0bU/ChhSVZt24dAwcOpFOnTv6EyTnHI488QlRUFNHR0f6n5xZVPm3aND777DPi4uKYPXs2GzduJDExkbi4OGJiYti2bRvTpk1j+/btxMXF8cgjj5CSksKAAQMYN26cf1LLESNG0LNnT7p37868efP8MdavX5+HH36Y+Ph4Bg0aREl9NmPGDH72s59d0C6AP/7xj/Tq1YuYmBimT58O5CYiXbt2ZcqUKcTHx7N7926eeeYZ/5WxadNyLwRu376dIUOG0LNnT66++mq2bNkCQHJyMg888AB9+vShY8eO/qkVCn4uKSkp3HTTTRfEm5mZyejRo+nVqxe9evXi888/L1UfiohI9RDwc3DMrDYwDPhNaQ7gnJsHzIPcJxmXWOGfiyC9+CsUtIqGG2b+Z3nfevh7EWHFjYMe44vc1YYNG+jZs2eJYcXGxrJs2TJuvfVWdu/ezdq1a9m9e7d/KoKJEydSq1YtRo8ezWOPPVboBI7r16/nyy+/5Pjx4/To0YOhQ4eSmppKeno669at4+DBg/75qb744otCy2fOnJnvSbz3338/Dz74IOPHj+fMmTPk5OQwc+ZMNmzY4J8eICUlhdWrV7Nhwwb/jOKvvfYaERERnDx5kl69ejF69GiaNm3K8ePHiY+P53/+53948skn+e1vf1viHFSFtWvDhg1s27aN1atX45xj2LBhfPrpp1x22WVs3bqV+fPn89JLL/Hhhx/y7rvv8tVXX1GvXj0OHz4M5M7rNHfuXDp16sRXX33FlClT/BOJ7tu3j1WrVrFlyxaGDRvGmDFjLvhc8g4h5vXggw/y0EMP0a9fP3bt2sXgwYPZvHlzif0vIiLVS2ke9HcD8LVz7vzYzX4zi3TO7TOzSODCG0/K4sgu+Neq0tU5dbToOu0vnDCxLCZNmsTmzZtJSEigXbt29OnTh9DQ3I9v0aJFtGnThqysLEaPHs3ChQu54447LtjH8OHDqVu3LnXr1mXAgAGsXr2aVatWcdttt1GrVi1atmxJ//79WbNmTZHlDRs2zLfP3r1789RTT7Fnzx5GjRrln+G8oMTERH9yAzBnzhyWLl0KwO7du9m2bRtNmzYlJCSEW265BYDbb7+dUaNGlfjZFNWuf/zjH/To0QPIfYjetm3buOyyy2jXrp1//quPPvqIiRMnUq9ePQAiIiI4duwYX3zxBWPHjvUf4/Tp0/73I0aMICQkhG7dul0wlFiSjz76iE2bNvmXf/zxR7KysmjQoEExtUREpLopTYJzG/8ZngJ4D5gAzPS9LquQiBpfBu1KSEpaRedfDm9UdJ3GlxW7q+7du/uHOYoTGhrK7Nmz/ct9+vTxJxPnZ7lu0KAB48aNY/Xq1YUmOAWv6pgZRc0FFugcYePGjePKK6/kb3/7G4MHD+aVV16hY8eOF2yXd7LJlJQUPvroI1JTU6lXrx5JSUmcOnWq0P0XdiWqpG3Ot+s3v/kN99xzT751GRkZ+WJxzl1Q/9y5czRu3LjICSrr1KmTr35pnDt3jtTUVOrWrVuqeiIiUr0EdA+OmdUDrgPeyVM8E7jOzLb51s0srG6p9RgPE/9W/M8NBQ4VGVP0tsUMTwEMHDiQ06dP57t3ZM2aNXzyySf5tjtx4gTHjx8HYPny5YSGhtKtWzeys7M5ePAgAGfPnuWDDz7wf/OqoGXLlnHq1CkOHTpESkqKf9hp8eLF5OTkkJmZyaeffkpiYmKR5Q0aNCArK8u/zx07dtCxY0ceeOABhg0bxvr16y/YpqCjR4/SpEkT6tWrx5YtW/jyyy/9686dO+dP+N5880369ctNHF944YUih6oKa9fgwYN57bXX/NMffP/994V+u+z666/ntdde48SJEwAcPnyYhg0b0qFDB5YsWQLkJjHr1q0rsj1AiW3Oe7y87dAs3yIi3hTQFRzn3AmgaYGyQ+R+q6paMzOWLl3K1KlTmTlzJuHh4bRv357nnnsu33YHDhxg8ODBhISE0KZNGxYuXAjkDp0MHjyYs2fPkpOTw7XXXstdd91V6LESExMZOnQou3bt4vHHH6d169aMHDmS1NRUYmNjMTOeeeYZWrVqVWR506ZNCQ0NJTY2luTkZE6dOsUbb7xBWFgYrVq14oknniAiIoK+ffsSFRXFDTfcwNChQ/PFMWTIEObOnUtMTAydO3f2DxdB7pWejRs30rNnTxo1auS/uXnLli307ds34Ha1bt2azZs307t37kzD9evX54033qBWrVoXxJKenk5CQgK1a9fmxhtv5Omnn2bRokXce++9/P73v+fs2bPceuutxMbGFtmPMTEx+T6X80NjBc2ZM4df/OIXxMTEkJ2dzTXXXMPcuXOL3K+IiFRPVtpL/OWRkJDgCn7zaPPmzXTt2rXKYpDi1a9fv9BJJ2+66Sbeeecdateuna98xowZ1K9fn1/+8pdVFaJn6Hdfqqtb/pzKVzsPc2WHCBbf0zvY4UgNZ2ZrnXMJBcs1m7gE5Py3k0RERKoDJTiST2FXb4ozY8aMyglERESkHC6KyTarcphM5GKg33kRkcoV9AQnPDycQ4cO6Q++1BjOOQ4dOkR4eHiwQxER8aygD1G1bduWPXv2lDglgIiXhIeH07Zt22CHISLiWUFPcMLCwvI9YVdERESkvII+RCUiIiJS0ZTgiIiIiOcowRERERHPUYIjIiIinqMER0RERDxHCY6IiIh4jhIcERER8RwlOCIiIuI5SnBERETEc5TgiIiIiOcowRERERHPCfpcVCIiUnHe31r5xzh0svKPIVJeuoIjIiIinhNQgmNmjc3sbTPbYmabzay3mUWY2XIz2+Z7bVLZwYqIiIgEItArOM8Df3fOdQFigc3ANGCFc64TsMK3LCIiIhJ0JSY4ZtYQuAZ4FcA5d8Y5dwQYDizwbbYAGFE5IYqIiIiUTiBXcDoCmcB8M/unmb1iZpcALZ1z+wB8ry0qMU4RERGRgAWS4IQC8cCfnHM9gOOUYjjKzO42szQzS8vMzCxjmCIiIiKBCyTB2QPscc595Vt+m9yEZ7+ZRQL4Xg8UVtk5N885l+CcS2jevHlFxCwiIiJSrBITHOfcD8BuM+vsKxoEbALeAyb4yiYAyyolQhEREZFSCvRBf/cDi8ysNrADmEhucvSWmf0c2AWMrZwQRUTkYnToZNU8WPDmziVvI1JQQAmOcy4dSChk1aAKjUZERESkAuhJxiIiIuI5SnBERETEc5TgiIiIiOcowRERERHPUYIjIiIinqMER0RERDxHCY6IiIh4jhIcERER8RwlOCIiIuI5SnBERETEc5TgiIiIiOcowRERERHPUYIjIiIinqMER0RERDxHCY6IiIh4jhIcERER8RwlOCIiIuI5SnBERETEc5TgiIiIiOcowRERERHPUYIjIiIinhMayEZmlgFkATlAtnMuwcwigMVAeyAD+Klz7t+VE6aIiIhI4EpzBWeAcy7OOZfgW54GrHDOdQJW+JZFREREgq48Q1TDgQW+9wuAEeWORkRERKQCBJrgOOAfZrbWzO72lbV0zu0D8L22KKyimd1tZmlmlpaZmVn+iEVERERKENA9OEBf59xeM2sBLDezLYEewDk3D5gHkJCQ4MoQo4iIiEipBHQFxzm31/d6AFgKJAL7zSwSwPd6oLKCFBERESmNEhMcM7vEzBqcfw9cD2wA3gMm+DabACyrrCBFRERESiOQIaqWwFIzO7/9m865v5vZGuAtM/s5sAsYW3lhioiIiASuxATHObcDiC2k/BAwqDKCEhERESkPPclYREREPEcJjoiIiHiOEhwRERHxHCU4IiIi4jlKcERERMRzlOCIiIiI5yjBEREREc9RgiMiIiKeowRHREREPEcJjoiIiHiOEhwRERHxHCU4IiIi4jlKcERERMRzlOCIiIiI5yjBEREREc9RgiMiIiKeowRHREREPEcJjoiIiHiOEhwRERHxHCU4IiIi4jlKcERERMRzAk5wzKyWmf3TzD7wLUeY2XIz2+Z7bVJ5YYqIiIgErjRXcB4ENudZngascM51Alb4lkVERESCLqAEx8zaAkOBV/IUDwcW+N4vAEZUaGQiIiIiZRToFZzngF8B5/KUtXTO7QPwvbYorKKZ3W1maWaWlpmZWZ5YRURERAJSYoJjZjcBB5xza8tyAOfcPOdcgnMuoXnz5mXZhYiIiEiphAawTV9gmJndCIQDDc3sDWC/mUU65/aZWSRwoDIDFREREQlUiVdwnHO/cc61dc61B24FPnbO3Q68B0zwbTYBWFZpUYqIiIiUQnmegzMTuM7MtgHX+ZZFREREgi6QISo/51wKkOJ7fwgYVPEhiYiIiJSPnmQsIiIinqMER0RERDxHCY6IiIh4jhIcERER8RwlOCIiIuI5SnBERETEc5TgiIiIiOcowRERERHPUYIjIiIinqMER0RERDxHCY6IiIh4jhIcERER8RwlOCIiIuI5SnBERETEc5TgiIiIiOcowRERERHPUYIjIiIinqMER0RERDwnNNgBiIiIFOf9rVV3rJs7V92xpHLpCo6IiIh4jhIcERER8ZwSExwzCzez1Wa2zsw2mtlvfeURZrbczLb5XptUfrgiIiIiJQvkCs5pYKBzLhaIA4aY2VXANGCFc64TsMK3LCIiIhJ0JSY4Ltcx32KY78cBw4EFvvIFwIjKCFBERESktAK6B8fMaplZOnAAWO6c+wpo6ZzbB+B7bVFE3bvNLM3M0jIzMysobBEREZGiBZTgOOdynHNxQFsg0cyiAj2Ac26ecy7BOZfQvHnzMoYpIiIiErhSfYvKOXcESAGGAPvNLBLA93qgooMTERERKYtAvkXV3Mwa+97XBa4FtgDvARN8m00AllVSjCIiIiKlEsiTjCOBBWZWi9yE6C3n3Admlgq8ZWY/B3YBYysxThEREZGAlZjgOOfWAz0KKT8EDKqMoERERETKQ08yFhEREc9RgiMiIiKeowRHREREPEcJjoiIiHiOEhwRERHxHCU4IiIi4jlKcERERMRzlOCIiIiI5yjBEREREc9RgiMiIiKeowRHREREPEcJjoiIiHhOILOJi4hIOb2/NdgRiNQsuoIjIiIinqMER0RERDxHCY6IiIh4jhIcERER8RwlOCIiIuI5SnBERETEc5TgiIiIiOeUmOCY2aVmttLMNpvZRjN70FceYWbLzWyb77VJ5YcrIiIiUrJAruBkAw8757oCVwG/MLNuwDRghXOuE7DCtywiIiISdCUmOM65fc65r33vs4DNQBtgOLDAt9kCYEQlxSgiIiJSKqW6B8fM2gM9gK+Als65fZCbBAEtiqhzt5mlmVlaZmZmOcMVERERKVnACY6Z1Qf+HzDVOfdjoPWcc/OccwnOuYTmzZuXJUYRERGRUgkowTGzMHKTm0XOuXd8xfvNLNK3PhI4UDkhioiIiJROIN+iMuBVYLNz7tk8q94DJvjeTwCWVXx4IiIiIqUXGsA2fYGfAd+YWbqv7FFgJvCWmf0c2AWMrZQIRUREREqpxATHObcKsCJWD6rYcERERETKT08yFhEREc9RgiMiIiKeowRHREREPEcJjoiIiHiOEhwRERHxHCU4IiIi4jlKcERERMRzlOCIiIiI5yjBEREREc8JZKoGERE575+L4MguaHwZ9Bgf7GikKqjPqyVdwRERKY30N+GTmbmvUjOoz6slJTgiIiLiOUpwRETKYFfWLib+fSJ/WP2HYIciVegPq/+gfq8mdA+OiEgZnDh7grT9acEOQ6rYlsNb1O/VhBIcEZHSaBUNwDE7TULLS+kS0SXIAUml8/U5raLpEtEEQP1eDSjBEREpjRtmApAAzA9uJFJVfH0O8OsghiGlo3twRERExHOU4IiIiIjnaIhKRKQ09q2HU0f5tzlONLucOrXq0Kxus2BHJZXJ1+eEN+Jg49aczjmtfq8GlOCIiJTG338D/1rF/gbNGNusHgktE5g/RHfjeJqvz2nXj0ciW5C2P039Xg1oiEpEREQ8p8QrOGb2GnATcMA5F+UriwAWA+2BDOCnzrl/V16YIiIXl8j6kfyu7y81TFHD/Dz65wz/yXD1ezUQyBDV68ALwF/ylE0DVjjnZprZNN+yvj0nIjVGo9qNGPGTEcEOQ6pYvzb9gh2CBKjEISrn3KfA4QLFw4EFvvcLgBEVG5aIiIhI2ZX1JuOWzrl9AM65fWbWoqgNzexu4G6Ayy67rIyHExGpHO9vLd32vU9CM+DgSUgtZV25+BX2+1AZfX5z54rZjxSt0m8yds7Nc84lOOcSmjdvXtmHExGpEsfOHOSj3S+x9sC7wQ5FqtDaA++q36uJsiY4+80sEsD3eqDiQhIRufhlnT3Ix9//ia8zlwU7FKlCX2cuU79XE2VNcN4DJvjeTwDU0yIiInLRCORr4n8FkoBmZrYHmA7MBN4ys58Du4CxlRmkiMjFInXQ3/zvnw5iHFJ18vb5XUGMQ0qnxATHOXdbEasGVXAsIiIiIhVCTzIWERERz1GCIyIiIp6jyTZFREqh94qhNDuwik31I7mnTTsubRDDbZ3+GOywpBKd7/ODLfrxwGXd2J21Xv1eDegKjohIGeS4sxw5s5djZw4GOxSpQsfOHFS/VxO6giMiUgb1QhsR3+wGmtftEOxQpApd0bgvTeq0Vr9XA0pwRETKoFGdSMb85KlghyFVrH+bO4MdggRIQ1QiIiLiOUpwRERExHM0RCUiUgZnck6w4+gawkMb0PqSLsEOR6rI3uNbOJWdpX6vBnQFR0SkDA6d2sUrmyfxt4w/BDsUqUJ/y/iD+r2aUIIjIiIinqMhKhGRUtgY/9+EnTnK3px/c+clkYSHNgh2SFLJzvf52dqNGFq7tn+ISi5uSnBERErhxyYxANQBOgY3FKki5/scoHUQ45DS0RCViIiIeI4SHBEREfEcDVGJyEXn/a3BjqBo3ddOo+GRb9gUZrzaOpbmdTvo6bYed77Pf2wczUutfkLmyZ3q92pAV3BEREqh4ZFvaHZgFU2ObObrg+/x7ZHPgx2SVLLzfd7wyDd8e+Rz9Xs1oQRHRKQMalkYjWu3pn7tZsEORapQ/drN1O/VhIaoRETKoEW9y/lV/N+CHYZUsds6/THYIUiAdAVHREREPEcJjoiIiHhOuRIcMxtiZlvN7Dszm1ZRQYmIiIiUR5kTHDOrBbwI3AB0A24zs24VFZiIyMVs3/EtPPplNC9vnBjsUKQKvbxxovq9mijPFZxE4Dvn3A7n3Bngf4HhFROWiIiISNmZc65sFc3GAEOcc3f6ln8GXOmcu6/AdncDd/sWo4ANZQ+32msGHAx2EEFSk9sOan9Nbn9NbjvU7PbX5LZD1bW/nXOuecHC8nxN3AopuyBbcs7NA+YBmFmacy6hHMes1mpy+2ty20Htr8ntr8lth5rd/prcdgh++8szRLUHuDTPcltgb/nCERERESm/8iQ4a4BOZtbBzGoDtwLvVUxYIiIiImVX5iEq51y2md0H/B9QC3jNObexhGrzyno8j6jJ7a/JbQe1vya3vya3HWp2+2ty2yHI7S/zTcYiIiIiFys9yVhEREQ8RwmOiIiIeE55nmRc4jQNZpZkZulmttHMPimprplFmNlyM9vme21S1vgqW1nbb2aXmtlKM9vsK38wz/YzzOx7X510M7uxqtpTGuXs+wwz+8a3Li1PeU3o+855+jbdzH40s6m+ddWi76Hk9pvZI3nascHMcswsori61aX/y9p2L5z3UO6+9/y5X0z/V/tzP4C2NzKz981sne93fGJJdSu9751zpf4h96bi7UBHoDawDuhWYJvGwCbgMt9yi5LqAs8A03zvpwF/KEt8lf1TzvZHAvG+9w2Ab/O0fwbwy2C3r7La7nufATQrZL+e7/tC9vMDuQ+oqhZ9H2j7C2x/M/BxSXWrQ/+Xs+3V+rwvb/t9y54/94trf4H9VKtzP8C/e4+e7zugOXDYt23QzvuyXsEJZJqGccA7zrldAM65AwHUHQ4s8L1fAIwoY3yVrcztd87tc8597XufBWwG2lRZ5OVXnr4vjuf7voBBwHbn3L8qNdqKV9opWm4D/hpA3erQ/2VuuwfOeyhf3xenOvQ9VFz7q+O5H0jbHdDAzAyoT26Ck11C3Urt+7ImOG2A3XmW93DhyXoF0MTMUsxsrZndEUDdls65fZD7BwFoUcb4Klt52u9nZu2BHsBXeYrvM7P1ZvbaRXqptrxtd8A/fOV35ymvUX1P7nOjCv7xu9j7HgJrPwBmVg8YAvy/AOpWh/4vT9vzrmtP9TvvofztrwnnPlB8/1M9z/1A2v4C0JXcB/5+AzzonDtXQt1K7fuyJjiBTNMQCvQEhgKDgcfN7IoA617sytP+3B2Y1Sf3l3+qc+5HX/GfgMuBOGAf8D8VG3aFKG/b+zrn4smdhf4XZnZNpUVaOSqi72sDw4AleepUh76H0p2/NwOfO+cOl6Huxag8bc/dQfU976H87a8J5/55RfV/dT33A2n7YCAdaE1uW14ws4YB1q0UZU1wApmmYQ/wd+fccefcQeBTILaEuvvNLBLA9xrI0EYwlKf9mFkYuX/kFjnn3jlfwTm33zmX48t6Xyb30t7Fplxtd87t9b0eAJbynzbWiL73uQH42jm3/3xBNel7KN0ULQX/p1rdz/3ytL26n/dQzvbXkHP/vMKu0kD1PfcDaftEcofmnXPuO2An0KWEupXb92W5cYfc/6HuADrwn5uGuhfYpiuwwrdtPXJnEY8qri7wR/LfcPRMRdxoVNE/5Wy/AX8Bnitkv5F53j8E/G+w21rBbb8EaODb5hLgC3JnpK8RfZ9n/f8CE6tb3wfaft92jcgdg78kkLrVof/L2fZqfd5XQPtrxLlfVPvzrKuW536Af/f+BMzwvW8JfE/ubOJBO+/L0+Abyf0mwHbgv3xlk4HJebZ5hNxvk2wg95JskXV95U3J/Ydhm+81ItgdW9HtB/qRe3luPbmX89KBG33rFpI7drme3Hm9IquyTVXQ9o6+X+51wMaa1ve+8nrAIaBRgX1Wi74vRfuTKeQPdXU/98vadi+c9+Vsf00694v63a/W535JbSd3aOofvrZsAG4vrm5V9L2mahARERHP0ZOMRURExHOU4IiIiIjnKMERERERz1GCIyIiIp6jBEdEREQ8RwmOiIiIeI4SHBEREfGc/w88G+mLRe8vKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "ax.vlines( bootstrap_train_mean, [0], 80, lw=2.5, linestyle='-', label='bootstrap train mean')\n",
    "\n",
    "ax.vlines(bootstrap_na_upper, [0], 15, lw=2.5, linestyle='dotted', \n",
    "          label='CI95 bootstrap, normal approx.', color='C2')\n",
    "ax.vlines(bootstrap_na_lower, [0], 15, lw=2.5, linestyle='dotted', color='C2')\n",
    "\n",
    "ax.vlines(bootstrap_percentile_upper, [0], 15, lw=2.5, linestyle='dashed', \n",
    "          label='CI95 bootstrap, percentile', color='C1')\n",
    "ax.vlines(bootstrap_percentile_lower, [0], 15, lw=2.5, linestyle='dashed', color='C1')\n",
    "\n",
    "ax.hist(bootstrap_train_accuracies, bins=7,\n",
    "        color='#0080ff', edgecolor=\"none\", \n",
    "        alpha=0.3)\n",
    "plt.legend(loc='upper left')\n",
    "\n",
    "plt.xlim([0.6, 0.81])\n",
    "plt.tight_layout()\n",
    "#plt.savefig('figures/bootstrap-histo-1.svg')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.4 Normal Approximation (No Bootstrap) Based on Test Set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7405787360493972 0.8007545972839363\n"
     ]
    }
   ],
   "source": [
    "clf.fit(X_train, y_train)\n",
    "acc_test_na = clf.score(X_test, y_test)\n",
    "ci_test_na = 1.96 * np.sqrt((acc_test_na*(1-acc_test_na)) / y_test.shape[0])\n",
    "\n",
    "test_na_lower = acc_test_na-ci_test_na\n",
    "test_na_upper = acc_test_na+ci_test_na\n",
    "\n",
    "print(test_na_lower, test_na_upper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.5 Bootstrapping the Test Set predictions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Avoids retraining the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7359666666666667 0.7947333333333333\n"
     ]
    }
   ],
   "source": [
    "clf.fit(X_train, y_train)\n",
    "\n",
    "predictions_test = clf.predict(X_test)\n",
    "acc_test = np.mean(predictions_test == y_test)\n",
    "\n",
    "rng = np.random.RandomState(seed=12345)\n",
    "idx = np.arange(y_test.shape[0])\n",
    "\n",
    "test_accuracies = []\n",
    "\n",
    "for i in range(200):\n",
    "    \n",
    "    pred_idx = rng.choice(idx, size=idx.shape[0], replace=True)\n",
    "    acc_test_boot = np.mean(predictions_test[pred_idx] == y_test[pred_idx])\n",
    "    test_accuracies.append(acc_test_boot)\n",
    "    \n",
    "    \n",
    "mean_test_accuracies = np.mean(test_accuracies)\n",
    "bootstrap_lower_test = np.percentile(test_accuracies, 2.5)\n",
    "bootstrap_upper_test = np.percentile(test_accuracies, 97.5)\n",
    "\n",
    "print(bootstrap_lower_test, bootstrap_upper_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.5 OOB Bootstrap with Percentile Method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same as \"4.1.2 Normal Approximation-based Bootstrap Interval (Based on OOB Estimates)\" but using mlxtend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean Bootstrap score 0.746178018090561\n",
      "Score Std 0.012612371356695583\n"
     ]
    }
   ],
   "source": [
    "from mlxtend.evaluate import bootstrap_point632_score\n",
    "\n",
    "bootstrap_scores = bootstrap_point632_score(clf, \n",
    "                                            X_train, y_train, \n",
    "                                            n_splits=200, \n",
    "                                            method='oob',\n",
    "                                            random_seed=12345)\n",
    "\n",
    "bootstrap_oob_mean = np.mean(bootstrap_scores)\n",
    "print('Mean Bootstrap score', bootstrap_oob_mean)\n",
    "print('Score Std', np.std(bootstrap_scores))\n",
    "\n",
    "bootstrap_oob_percentile_lower = np.percentile(bootstrap_scores, 2.5)\n",
    "bootstrap_oob_percentile_upper = np.percentile(bootstrap_scores, 97.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7204597701149426 0.7696229551148738\n"
     ]
    }
   ],
   "source": [
    "print(bootstrap_percentile_lower, bootstrap_percentile_upper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.6 .632 Bootstrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The .632 Bootstrap is the default setting of `bootstrap_point632_score`; it tends to be overly optimistic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean Bootstrap score 0.8052215333155874\n",
      "Score Std 0.009690206838191797\n"
     ]
    }
   ],
   "source": [
    "bootstrap_scores = bootstrap_point632_score(clf, \n",
    "                                            X_train, y_train, \n",
    "                                            n_splits=200,\n",
    "                                            random_seed=12345)\n",
    "\n",
    "bootstrap_632_mean = np.mean(bootstrap_scores)\n",
    "print('Mean Bootstrap score', bootstrap_632_mean)\n",
    "print('Score Std', np.std(bootstrap_scores))\n",
    "\n",
    "bootstrap_632_percentile_lower = np.percentile(bootstrap_scores, 2.5)\n",
    "bootstrap_632_percentile_upper = np.percentile(bootstrap_scores, 97.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7855672805949966 0.8229009076326003\n"
     ]
    }
   ],
   "source": [
    "print(bootstrap_632_percentile_lower, bootstrap_632_percentile_upper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  4.1.7 .632+ Bootstrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The .632+ Boostrap method attempts to address the optimistic bias of the regular .632 Boostrap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean Bootstrap score 0.7971839079424146\n",
      "Score Std 0.010569723397996472\n"
     ]
    }
   ],
   "source": [
    "bootstrap_scores = bootstrap_point632_score(clf, X_train, y_train, \n",
    "                                            n_splits=200, \n",
    "                                            method='.632+',\n",
    "                                            random_seed=12345)\n",
    "\n",
    "bootstrap_632plus_mean = np.mean(bootstrap_scores)\n",
    "print('Mean Bootstrap score', bootstrap_632plus_mean)\n",
    "print('Score Std', np.std(bootstrap_scores))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "bootstrap_632plus_percentile_lower = np.percentile(bootstrap_scores, 2.5)\n",
    "bootstrap_632plus_percentile_upper = np.percentile(bootstrap_scores, 97.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7756256434146681 0.8164960245443756\n"
     ]
    }
   ],
   "source": [
    "print(bootstrap_632plus_percentile_lower, bootstrap_632plus_percentile_upper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAADQCAYAAAD4dzNkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8Y0lEQVR4nO3de7xVVbn/8c9XRN2KSl7yKKkoKqiAoFsT8YLG0cxSNI3STPSYl1KPlZhWJynrJ6csy/tdtLwlIpqk0AkveUVugvfymltLVFDRnXL5/v4YY8FisdbaFzbsveB5v177tdcac8wxnznW2uuZY8y515RtQgghhFBbVmvvAEIIIYTQcpHAQwghhBoUCTyEEEKoQZHAQwghhBoUCTyEEEKoQZHAQwghhBoUCTyEENqJpIGS/iZprqQhku6RdEyFut0lWdLqKzrOSnLcW7d3HE2RNEjS6+0dR1uLBB5CqBmStpc0UdJ7kv4u6dCiZYUEN7fo53+Klh8p6U1JL0saVFTeQ9Ijkjqt2L0B4KfAxba72B5r+0Db17dDHK2S436pOXXza7PN8o5pVdJhjuRCCKGaPPK8E7gc+E9gH+CPkvrbfqGoalfb88usOxLYGdgFuBjonRdfCHzX9oLlvAvlbAk83Q7brSmSVi99TUOMwEMItaMXsBlwge0FticCDwNHN2PdDYEG228C/wdsDSDp8Fz+WFMNSPqmpGclfSDpGUk75/LtJd0vaY6kpyUdXLTOKEmXSBqX13tcUo+87MUcxx/zbMGauZ3j8/JOks6X9Lakl4CDSuJZX9I1eVahQdLPCrMIkoZJeiivPzvPOhxYtO4Gkq6T9EZePrZo2RclTc/784ikvlX6ZNGouol9fTCv8mTe16FNbUvSK5K+L2kG8KGkH0kaXbL930q6MD8+tuj1eUnSiVXi/n7usw8kPS/pc5Xqdmi24yd+4id+OvwP0AeYC6io7M/AHflxd8BAA/A6cB2wUV62GvAC8BngS8ATQBdgOrBhM7Z9RG53V0DANqTRc2fg78APgDWA/YAPgJ55vVHAu8BupBnPG4Fbitp9BRhc9Px+4Pj8+CTgOWBzYAPgvrx/q+flY4ErgHWATwOTgBPzsmHAPOCbQCfgZOCNQt8B44BbgU/lfdgnl+8MvAV8Nq93TI5xzQr9YmCbZu7rorrN2VZ+PD3vf13u74+A9fLyTsCbwO75+UFAj/z67JPr7pyXDQJez497Av8ANit63/Ro7/d3a35iBB5CqBXPkT7wh0vqLGl/0gf12nn526QEuyVpmnxdUhLB9kJSEhsNnEFKbD8FLgL6SLpP0nhJvSnveOAXtp9w8nfbrwK7kw4ERtr+xGlW4G7ga0XrjrE9yWkK+EagXzP39yvAb2z/w/a7wHmFBZI2AQ4ETrf9oe23gAuArxat/6rtq5xODVwPbApsImnTvO5Jtmfbnmf7gbzON4ErbD/uNMtxPfBx3s/maMm+NmdbF+b9b8z9PRUYkpftB3zkPHtie5ztF/Pr8wAwAdirzHYXAGsCO0jqbPsV2y82c/86lEjgIYSaYHse6cP7IOCfwPeAP5BG29iea3uy7fm2/wWcAuwvab28/C+2d7e9D7AQqCeNGn9HGrGeC1xdYfObA+U+5DcD/pEPEApeBboVPf9n0eOPSAm/OTYjjRSL2y0ojP7fzNPPc0ij8U+X267tj/LDLqR9edf27DLb3BL4XqHN3O7mOZbmaMm+Nmdb/yhZ5yYWHxwdmZ8DIOlASY9Jeje39QVgo9KN2v47cDowAnhL0i2Smrt/HUok8BBCzbA9w/Y+tje0fQDpHPKkStXzbxUXShLpIrbTSB/wnfLo7gmg0vnef5CmZ0u9AWwuqfizdAvSdPuyepOU0IrbLY7nY9Ipgq75Zz3bOzaj3X8AG0jqWmHZz4va7Gp7bds3t3YnmoijqW2V3i7zNmCQpM8Ah5ITuKQ1gduB84FNbHcF/kTJa7+oUfsm23uSDiIM/G8b7tcKEwk8hFAzJPWVtJaktSWdQZoWHpWXfVZST0mrSdqQdHX5/bbfK2nmeGCa7enAO0CdpB2AfYFK/xJ1NXCGpF2UbCNpS+Bx4EPgzDytP4h0jv2WNtjdPwCnSfqMpE8BZxUWOF2MNwH4laT18j73kLRPU43mde8BLpX0qRz33nnxVcBJuS8laR1JB0latw3251/kiwdbuy3bs0jXCVwHvGz72bxoDdK0+Cxgfr5gb/9ybeT3yH456f8baCRNq9ecSOAhhFpyNGlk+hbwOeA/bX+cl20N3Eu6iOwp0gi1+Fw0kjYC/hv4H4B8rvYUYCLp39NOLbdR27cBPyeN+D4gXUC2ge1PgINJ55TfBi4FvmH7uTbY16uA8cCTpHO/Y0qWf4OUuJ4BZpPO72/azLaPJl3kVriu4HQA25NJ56Yvzm3+nXR6oS2MAK7P0+VfWYZt3QQMpmj63PYHpBmVP+S2jgTuqrD+mqR/KXybNOX/adJFiDWncEViCCGEEGpIjMBDCCGEGhQJPIQQQqhBkcBDCCGEGhQJPIQQQqhBcTOTEMISNtpoI3fv3r29wwghZFOmTHnb9sal5ZHAQwhL6N69O5MnT27vMEIImaRXy5XHFHoIIYRQgyKBhxBCCDUoEngr5a/++7mkF/I9aE/L5YdImpHvcTtZ0p65fPN8x6Nnle4Z/N+t3G53SY25/SfzPXR7trKtfpK+0IztHdma9mtVab9IOljSWfnxiPwVnisqlmHFN1pQukfyUjdoaEF7y7R+CCva2GkNDBw5ka3OGsfAkRMZO60tvmZ+5RAJvPWGkW400Mv29iz+7uO/ADvZ7gccx+K7G80Hvpfr7g58O3//8hIkvdKMbb9ou5/tnUi3CWzt1wD2I92xp5rupK8lXIqkDnMNhaRObdhcP4r6xfZdtke2YfstMYzm3wkqhJXK2GkNnD1mJg1zGtON3uc0cvaYmZHEsw7zAVyDTgaOLNxGMN+PF9tzi+qsQ76bTr6BwJv58QeSniXdcvCZZYxjPdJ3/yJpLeAy0m0S5wPftX1fuXLgYdL9kOvyLMF5pO8F/m1u18DepO8M3l7SdNLBwmzS7RzXAtaRdDBwJ/Ap0u0Nf2T7TkndSd9L/TjQH3iB9B3RhdsaLkXSKNLNBXYENsnx352T80hgEOl7jC+xfUW+ccQ5pH7tJ6kP6a5CB+T4r7J9kaRdgF+Tbm34NjDM9puS7s/x7Qt0Bf4rPy/tlzqg3vYpJfH2AC4BNibdOvGbpd+BLWkEsBXpO6q3y32/O+m7sxuAL9meVy5GYCDpNbtRUiMwIDd7qqQv5f4+wvZzkjYAriV9H/hHwAm2Z+SbetycY5xEhbszhVXP0Csebe8QmjTttTl8smDhEmWN8xZw5ugZ3DzptXaKamm3njig6UrLQYzAW68HMDRPk98jadvCAkmHSnoOGEcahS8hJ7f+pGTRqm3nKfQXSQnh17n82wC2+5Bu4nB9Tt5LlZNe+x8Dt+bR/K3AGcC38+zBXqS79JwF/DXXuSBvZwBwjO39SAn3UNs7kxLhryQVkkRP4ErbfYH3gW81Y9+6A/uQDhIuz/H/F/Ce7V2BXYFvStoq198N+KHtHYATSMmyf97mjZI6AxcBh9vehZTkfl60vdVt70a6mcM5+eYUpf1SyZXAqbndM0g3siinR96fQ4DfA/fl16IROKhSjLZHA5OBo3Isjbm9t3N/X5a3C/AT0h22+pJmZG7I5ecAD9nuT7q5Q/EtKReRdEJ+L0+eNWtWlV0OYcUpTd5Nla9qYgTeemsC/7ZdL+kw0ofuXgC27wDuyLfoO5d05xwAJHUh3bf2dNvv57IfAkfkKpvl0S7Aw7a/XWbbL+Yki6ShpETyeWBPUiIgj8peJY36KpWXehj4taQbgTG2X1+ci5fwZ9vvFnYJ+H95XxeSZhU2ycv+Yfvh/Pj3pLsFnV+uwSJ/yLMaf5P0EtCLdFvAvpIOz3XWB7YFPgEm2X45lw8GLs93mML2u5J6A72BP+d96USeCckKd3iaQjp4aJb8Ou4B3FbUR2tWqH5PHmXPzNu/N5fPzNvs2USMpYpjPiw/3hP4MoDtiZI2lLQ+aRblsFw+TtLscg3avpL0PqK+vj7ucLQKaK9RY0sMHDmRhjmNS5V361pXE/Evb5HAW+91UiIGuIN0f9ol2H5Q6R69G9l+O4+0bgdutD2mqN7PyaNCSa8UknMz3VW07UrTo82aNrU9UtI40vnfxyQNrlD1w6LHR5GmZ3fJSeoV0vQ65NMHxZtoThhlnos00h1fvCBPoRfHojLrC3jadqW/9sKtKBfQsr+H1YA5zXytPgawvVDSPC++BeDCvM2mYizbHkvGXO41dsnvEGrK8AN6cvaYmTTOW3y77rrOnRh+QKuu213pxBR6640F9suP9yGd40XSNoUpZEk7k+7X+04uuwZ41vavl26u1fYEXsyPHyQlVCRtR5oufb5K+QfAuoWGJPWwPdP2/5KmbnuV1iljfeCtnLz3BbYsWraFpEJS+hrwUN7OeZIOrdDeEZJWy+eXt85xjgdOzgdASNpO0jpl1p0AnFS4uC6fF34e2LgQh6TOknassj80Y5/JsycvSzoitytJOzXRbiXVYmwylqz4NR5EmmZ/v6T8QNK1CiHUhCH9u3HeYX3o1rUOkUbe5x3WhyH9u7V3aB1CjMBbQNKfgONtv0G6qOpGSd8B5gLH52pfBr4haR7pHOdQ284XRB0NzCyaIv+B7T+1IpQeuQ2RppEL276UdN54JulitWG2P5ZUqfw+4Kzc1nnAnjkJLyBdXHcPaZQ4X9KTwCjyBXNFbgT+KGkyMB0ovojrWeAYSVcAfyOdswXoQ5o5KOd54AHSNPxJtv8t6WrSVPPUfCA0CxhSZt2rSacGZuT+v8r2xXnq/cI8pbw68Bvg6QrbByjtl0qOAi6T9CPSBWW3AE9WqV+W7U+qxDiK9NoVX8RWzgjgOkkzSBexHZPLfwLcLGkqqV87zpU/ITTDkP7dImFXoMWzeSG0nXyh3t22e5dZNt72AWXKR+V1Ri//CEMl9fX1jq9SDaHjkDTFdn1peUyhhxWuXPIOIYTQMjGFHpYL26+QrqxuyTrDlkswIYSwEooReAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgS+HKi5OeSXpD0rKTTcvkhkmZImi5psqQ9c/nmku7LdZ+W9N8V2h0hqSGv/5ykyyS16nWUNEzSZk3UGSJph9a0X6tK+0XS1YU+kPSKpI1WUBxdJX2r6PkgSXcvQ3vLtH4ItWjstAYGjpzIVmeNY+DIiYyd1tDeIbWZSODLzzBgc6CX7e2BW3L5X4CdbPcDjgOuzuXzge/lursD366SOC/I6+8A9AH2WYYYqyZwYEjezlIkrd7K7ba5No5lGEX9Yvt428+0YfvN1RX4VlOVQgjljZ3WwNljZtIwpxEDDXMaOXvMzJUmiXeYD+CV0MnAkbYXAth+K/+eW1RnHcC5/E3gzfz4A0nPAt2AaoljDWAtYDaApH7A5cDawIvAcbZnlysHPgfUAzdKagQGAOcAB5MOJiYAY/LzfST9CPgycA3wCDAQuEvSC8CPcizvAEfZ/pekEUCPvA+bA7+wfVW1DpM0F7gC2Dfv01dtz5LUA7gE2Bj4CPim7eckjQLeBfoDUyVdlvdzY2ABcITtFyUNB74CrAncYfscSd2Be4CHgD2ABuAQ4KAy/XIPcIbtySXxfh04Le/748C3bC8oqfMKcFPep87ACcB5wDbAL21fnustFSMwEughaTrwZ2Ac0EXSaKA3MAX4um1L+hxwPulv+gngZNsfS/o88BvgbWBqtf4PYVkMveLR9g5hKdNem8MnCxYuUdY4bwFnjp7BzZNea6eokltPHLDMbcQIfPnpAQzN0+T3SNq2sEDSoZKeI30gH1e6Yk4u/UlJoZzv5A/1N4EXbE/P5TcA37fdF5hJSshly22PBiaTEm4/oA44FNgx1/uZ7UeAu4DhtvvZfjG319X2PrZ/RUqAu9vuT5plOLMozr6khDgA+HFT0/WkA5qptncGHiiK/0rgVNu7AGcAlxatsx0w2Pb3gBuBS2zvRErKb0raH9gW2A3oB+wiae+87ra5/o7AHODLpf1iu7FcoJK2B4YCA3P/LQCOqrBf/7A9APgrMAo4nDTL8tPcVqUYzwJezHEMz231B04nzYpsDQyUtFZud6jtPqQkfnIuvwr4ErAX8B8V4kPSCfm9OnnWrFmVqoVQU0qTd1PltSZG4MvPmsC/bddLOgy4lvQhiu07gDvyh/S5wODCSpK6ALcDp9t+v0LbF9g+X1JnYLSkr5JGiV1tP5DrXA/cJmn9cuVl2nwf+DdwtaRxQLVzpbcWPf4McKukTUkj0ZeLlt2ZE2CjpPtICWpslXYXFrX9e2BM7o898r4U6q1ZtM5tthdIWhfolvsW2/+GRclxf2Bart+FlCxfA14uOviZAnSvElupzwG7AE/kuOqAtyrUvSv/ngl0sf0B8IGkf0vqmuOrFGOpSbZfz/s2Pcf8Qd6XF3Kd64FvA/fn8r/l+r8nzQAsxfaVpAMl6uvrXXXPQyijLUaUbW3gyIk0zFn6GLxb17oOGW9LxQh8+XmdlIgB7iCNRpdg+0HSFOlGADkh3w7caHtMUxuwPQ+4F9i7qbrNaGs+KcHeTjrvfW+V6h8WPb4IuDiP/E4kTekvarZ0My0Ni/QenZNHoYWf7cvEoqVXX1R+XtG629i+Ji/7uKjeAlp2QCvg+qJ2e9oeUaFuYTsLS7a5MG+zWoyV2iqOudK+Q8v7PISVxvADelLXudMSZXWdOzH8gJ7tFFHbigS+/IwF9suP9wFeAJC0jfKQTdLO5HPHuewa4Fnbv27OBvI6e5CmWd8DZkvaKy8+GnigUnl+/AGwbm6rC7C+7T+Rpmj7ldapYH3S+WOAY0qWHSJpLUkbAoNI52bJpw/KWY00vQxwJPBQnoV4WdIRhX2WtFPpirne65KG5HprSlobGA8cl/cPSd0kfbrK/kDT+wzpYsTDC21J2kDSlk2sU0mlGJsTB8BzQHdJ2+Tnhdf4OWCrfA0BwNdaGV8INWlI/26cd1gfunWtQ6SR93mH9WFI/27tHVqbiCn0NiTpT8Dxtt8gXYB0o6TvAHOB43O1LwPfkDQPaCSdt7TSv5MdDczMU6MAP8gJtdR38gVUnYEZLD4nfAxweU5cLwHHNlE+Kpc3AgcCd+bzpgK+k+vcAlyl9G9wheRabARpersBeAzYqmjZJNJ5/i2Ac22/kWcbKo0YPwR2lDQFeI90jhnSueXL8oV0nXNMT5ZZ/2jgCkk/BeaRLmKbkM9XP5qPm+YCXyeNXisZxeJ+KTvPZvuZHM8EpX/jm0eatn61SrtlVYoxX4D3sKSnSKdIxlVY/9+SjiW9DoWL2C7PF7GdAIyT9DbpeoXeLY0vhFo2pH+3lSZhl5IdM2yh7eWr0OfaPr+k/IvA1rYvLLPOXNtdVlCIoYL6+npPnjy56YohhBVC0hTb9aXlMQIPK5Tt+CKREEJoA5HAw3JR5YKuauvE6DuEEJopLmILIYQQalAk8BBCCKEGRQIPIYQQalAk8BBCCKEGRQIPIYQQalAk8BBCCKEGRQIPIYQQalAk8BBCCKEGRQIPIYQQalAk8BBCCKEGRQIPIYQQalAk8BBCCKEGRQIPIYQQalAk8BBCCKEGRQIPIYQQalAk8BokaYGk6ZKelDRV0h6tbKe7pCObqNNV0rdaF2ltKu0XSfWSLsyPh0m6eAXGMkTSDkXP75dUvwztLdP6IYSOIxJ4bWq03c/2TsDZwHmtbKc7UDWBA12BsglcUqdWbrfNKWmr93N3ivrF9mTbp7VR2y01BNihqUph2Yyd1sDAkRPZ6qxxDBw5kbHTGto7pBCaFAm89q0HzIZFSeyXkp6SNFPS0GrlwEhgrzya/46kHSVNys9nSNo21+mRy34paZCk+yTdBMzM7Y+VNEXS05JOKAQmaa6kX+VZgr9I2rjajkgaIel3kiZK+pukbxYtGy7piRzXT3JZd0nPSroUmApsLunMvI9PShqZ6/WQdG+O8a+SeuXyUZIulPSIpJckHV6hXwZJurtMvBtLuj3H9YSkgWXqDMv980dJL0s6RdJ3JU2T9JikDSrFmGdWDgZ+mWPpkZs9Ir9OL0jaK6+/lqTr8r5Pk7RvLq+TdEvut1uBumqvwapo7LQGzh4zk4Y5jRhomNPI2WNmRhIPHd7q7R1AaJU6SdOBtYBNgf1y+WFAP2AnYCPgCUkPAntUKD8LOMP2FwEkXQT81vaNktYAOuU6vW33y3UGAbvlspfzdo+z/a6kutz27bbfAdYBptr+nqQfA+cApzSxb32B3fO60ySNA3oD2+btCrhL0t7Aa0BP4Fjb35J0IGnE+lnbHxWSI3AlcJLtv0n6LHBpUZ9tCuwJ9ALuAkaX6ZdBFWL9LXCB7YckbQGMB7YvU6830J/0ev0d+L7t/pIuAL4B/KZcjLb3k3QXcLft0TkWgNVt7ybpC7lPBwPfBrDdJx+gTJC0HXAy8JHtvpL6kg502sTQKx5tq6ba1bTX5vDJgoVLlDXOW8CZo2dw86TX2imqlrv1xAHtHUJYwSKB16bGooQ6ALhBUm9SIrrZ9gLgX5IeAHatUv5+SbuPAj+U9BlgTE4m5bY/qSh5A5wm6dD8eHNSsn0HWAjcmst/D4xpxr7dabsRaJR0Hylp7wnsD0zLdbrkbbwGvGr7sVw+GLjO9kcA+aCiC+kA5raifVmzaHtjbS8EnpG0STPiKzYY2KGo3fUkrWv7g5J69+WyDyS9B/wxl88E+jYjxlKFfpxCmu6H1EcXAdh+TtKrwHbA3sCFuXyGpBnlGswzJycAbLHFFtX2eaVTmrybKg+ho4gEXuNsPyppI2Bj0ui0nErlpW3dJOlx4CBgvKTjgZfKVP1wUcNpdDoYGJBHvfeTRpplN9GcMMo8F3Ce7SuKF0jqXhxLrle6/mrAnMIBTxkfl6zfEquR9ruxiXrF21hY9Hwh6W+wqRgrtbeAxX/D1WJvst9tX0maBaC+vr45r9NKM+IbOHIiDXOWfgm7da1bafYxrJziHHiNy9OlnUgj3geBoZI65fPNewOTqpR/AKxb1NbWwEu2LyRNJ/ctrVPG+sDsnLx7kaa/C1YDCueVjwQeyts5RVKlqfRD8vncDYFBwBOkqenj8kgVSd0kfbrMuhNyvbVzvQ1svw+8LOmIXCZJO1XZH5qxz8XbW7Qfkvo1Y52lNBFjc2N5EDgqr78dsAXwfEl5b9JrGooMP6AndZ2XvB6zrnMnhh/Qs50iCqF5YgRemwrnwCGNvI6xvUDSHcAA4EnSqOtM2/+sUv4OMF/Sk8Ao0sj565LmAf8EfpqnoR+W9BRwDzCuJJZ7gZPy1OzzwGNFyz4EdpQ0BXgPKFw81wt4uMK+Tcrb2AI41/YbwBuStgcezVPMc4Gvk0agi9i+NyfRyZI+Af4E/ICUwC6T9COgM3BL7otKZpT0y7QK9U4DLsn7vjopWZ5Upd1qKsV4C3CVpNNYfDBUzqXA5ZJmAvOBYbY/lnQZcF2OcTqpf0ORIf27AfDL8c/zxpxGNutax/ADei4qD6Gjkt2s2bIQWkzSXNtdypTfDRxm+5OS8hHAXNvnr6AQQxn19fWePHlye4cRQsgkTbG91Pc3xAg8rHCFq7tDCCG0XiTwsNyUG303UX/EcgolhBBWOnERWwghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeFhuJC2QNF3Sk5KmStqjle10l3RkE3W6SvpW6yLtuPK+P9XCdf5D0i2SXpT0jKQ/SdquNW2FEDquSOBheWq03c/2TsDZwHmtbKc7UDWBA12BsglcUqdWbrfmSBJwB3C/7R62dwB+AGzSvpF1bGOnNTBw5ES2OmscA0dOZOy0hvYOKYQmRQIPK8p6wGxISUbSLyU9JWmmpKHVyoGRwF55NP8dSTtKmpSfz5C0ba7TI5f9UtIgSfdJugmYmdsfK2mKpKclnVAITNJcSb/KswR/kbRxtR2RtE/eznRJ0yStK6lLXndqjv2QXLe7pOckXZ3360ZJgyU9LOlvknbL9UZI+p2kibn8m2W22ynv2xN5v08sE96+wDzblxcKbE+3/dfmvlCrmrHTGjh7zEwa5jRioGFOI2ePmRlJPHR4q7d3AGGlVidpOrAWsCmwXy4/DOgH7ARsBDwh6UFgjwrlZwFn2P4igKSLgN/avlHSGkCnXKe37X65ziBgt1z2ct7ucbbflVSX277d9jvAOsBU29+T9GPgHOCUKvt1BvBt2w9L6gL8O5cfavt9SRsBj0m6K5dvAxwBnAA8QZpN2BM4mDQ6HpLr9QV2z/FMkzSuZLv/Bbxne1dJawIPS5pQtH8AvYEpVWJvc0OveHRFbq7NTXttDp8sWLhEWeO8BZw5egY3T3qtnaJaNreeOKC9QwgrQIzAw/JUmELvBXweuCFP8e4J3Gx7ge1/AQ8Au1YpL/Uo8ANJ3we2tN1YYfuTSpLbaZKeBB4DNge2zeULgVvz49/nOKp5GPi1pNOArrbnAwL+n6QZwP8B3Vg8bf2y7Zm2FwJPA3+xbdLMQPeidu+03Wj7beA+0gFIsf2Bb+SDoseBDYv2YZlIOkHSZEmTZ82a1RZN1ozS5N1UeQgdRYzAwwph+9E8Mt2YlOzKqVRe2tZNkh4HDgLGSzoeeKlM1Q8XNZxG5IOBAbY/knQ/aWag7Caa2P7IPDr+AmmkPZg0ct4Y2MX2PEmvFLX/cdHqC4ueL2TJv8HS7ZY+F3Cq7fFVwnsaOLxa/OXYvhK4EqC+vr7q/peq9dHewJETaZiz9DFgt651Nb9vYeUWI/CwQkjqRZrqfgd4EBiaz+luDOwNTKpS/gGwblFbWwMv2b4QuIs09bxEnTLWB2bn5N2LlHALVmNx0jsSeChv5xRJS02lS+qRR9T/C0wGeuX238rJe19gyxZ0T8EhktaStCEwiDTdXmw8cLKkzjmO7SStU1JnIrBm8Tl0SbtK2qcV8awShh/Qk7rOS17nWNe5E8MP6NlOEYXQPDECD8tT4Rw4pNHjMbYXSLoDGAA8SRplnmn7n1XK3wHm5+nvUaSR7dclzQP+Cfw0n9t+WOnfpO4BSs8f3wuclKe4nydNoxd8COwoaQrwHlC4eK4Xabq81Ok5SS8AnsnbWxf4o6TJwHTguRb2FaSDlXHAFsC5tt+Q1L1o+dWkKfep+VTELBafPwfAtiUdCvxG0lmk8/OvAKe3Ip5VwpD+3QD45fjneWNOI5t1rWP4AT0XlYfQUSmdigth1SVpru0uZcrvBg6z/ckKiGEEMNf2+ct7W02pr6/35MmT2zuMEEImaYrt+tLyGIGHUEHhqvcQQuiIIoGHVV650Xc7xDCivWMIIdSWuIgthBBCqEGRwEMIIYQaFAk8hBBCqEGRwEMIIYQaFAk8hBBCqEGRwEMIIYQaFAk8hBBCqEGRwEMIIYQaFAk8hBBCqEGRwEMIIYQaFAk8hBBCqEGRwEMIIYQaFAk8hBBCqEGRwEMIIYQaFAk8hBBCqEGRwEMIIYQatEITuKQFkqZLelLSVEl7tLKd7pKObKJOV0nfal2kLY7nYElntVFbwyRtVvT8FUkbtbKtJvupifV/0Np1S9oZImmHtmhrRSjtN0n1ki7Mj4dJungFxrJE30m6X1L9MrS3TOuHjmfstAYGjpzIVmeNY+DIiYyd1tDeIYUVZEWPwBtt97O9E3A2cF4r2+kONJWYugJlE7ikTq3cblm277I9so2aGwZs1lSlZupO0/1UTZskcGAIsFwTuJK2ej93p6jfbE+2fVobtd1SQ1jOfRdq19hpDZw9ZiYNcxox0DCnkbPHzIwkvopYvR23vR4wG9KHL/AL4EDAwM9s31qpHBgJbC9pOnA9MAG4DliDdFDyZeBcoEeu82dgHHAO8CbQD9hB0lhgc2At4Le2r8zxzAWuAPbNMX7V9ixJ9wPTgd1y/MfZniRpGFBv+xRJo4D3gXrgP4AzbY/OyeViYB/g5RzntbZHFzpE0uF5vRslNQID8qJTJX0J6AwcYfs5SesAFwF9SK/jCNt3lvRxaT9dmMsGAWsCl9i+QtKmwK15n1YHTgYOAuryuk/bPqoozk7ANTlW5/24QFIP4BJgY+Aj4JvABsDBwD6SfgR82faLlCFpBNAD6JZfl1/YviovGw58Jcd9h+1zJHUH7gHuy301RNJQ4GhgIXCP7bPKxZX7sOxrVabfpgFn2P5iSbwbA5cDW+Si020/XFJnGCkJdwJ6A78ivU+PBj4GvmD73eb2XW72CEmXkg5S/8v2XyWtBVyW92U+8F3b90mqI/1t7AA8C9SV6/taNPSKR9s7hHY37bU5fLJg4RJljfMWcOboGdw86bV2imrZ3XrigKYrhRWewAsJYS1gU2C/XH4YKanuBGwEPCHpQWCPCuVnUfSBKukiUgK+UdIapA/Ls4DetvvlOoNIibe37Zfzdo/LH551ue3bbb8DrANMtf09ST8mJf5T8jrr2N5D0t7AtaQP5VKbAnsCvYC7gNF5H7uTEu6nSR+m1xavlBP9KXnfJue4Ad62vXM+JXAGcDzwQ2Ci7eMkdQUmSfo/2x8WNVnaTycA79neVdKawMOSJuTYxtv+eU7Oa+ekcEqh/0r0A7rZ7p3b7ZrLrwROsv03SZ8FLrW9n6S7gLuLD1aq6AvsTnoNpkkaR+rjbUmvn4C7cv+/BvQEjrX9LUkHkpLlZ21/JGmDSnGx+L1X7rUq7bdBFWL9LXCB7YckbQGMB7YvU6830J/0vv878H3b/SVdAHwD+E1z+y6/H1a3vZukL5Dem4OBbwPY7iOpFzBB0nakg7GPbPeV1BeYWm5H8nvjBIAtttiiXJXQAZUm76bKw8plRSfwxqKEOgC4QVJv0gfozbYXAP+S9ACwa5Xy90vafRT4oaTPAGPyh2C57U8qSt4Ap0k6ND/enJQk3iGN3m7N5b8HxhStczOA7QclrVeUvIqNtb0QeEbSJrlsT+C2XP5PSfdV6KNyCtufQkq2APsDB0s6Iz9fizQSfLZKO/sDffNIH2B90j4/AVwrqXOOfXoT8bwEbJ0PnMaRkkUX0gHXbUV9v2Yz9q3UnbYbgcbcR7uR+m5/0kgYoEuO+zXgVduP5fLBwHW2PwLIB2dNxVXutWquwaSZnMLz9SSta/uDknr35bIPJL0H/DGXzyS9Hi3tu+L3Q/f8eE/SjAx5duFVYDtgb9LMC7ZnSJpRrsE8+3QlQH19vavtdEcRozQYOHIiDXMalyrv1rUu+mcV0G5T6LYfzRdnbUwaVZVTqby0rZskPU6a9h0v6XhSkim1aHSaR1WDgQF5tHY/KQmW3USFx+WeQ5oaXbSpkt+tUWhvAYtfM5Gmo59vQTsCTrU9fqkFaUR7EPA7Sb+0fUOlRmzPlrQTcABp5PcV4HRgToURe0uU618B59m+oiTm7hS9prle6fqrNRFXudequVYjvX+W/gStvI2FRc8Xkl7PpmKs1F7p+6GSmkjIoeWGH9CTs8fMpHHegkVldZ07MfyAnu0YVVhR2u3fyPI0XyfSiPdBYKikTvm84t7ApCrlHwDrFrW1NfCS7QtJ06B9S+uUsT4wOyfvXqRp24LVgMIo9UjgoaJlQ/M29yRNR7/XzF1+CPiypNXySG9QhXpNxV0wnnRuXDme/s1oazxwch5pI2k7SetI2hJ4K59vvgbYOdefV6hbLB94rWb7duB/gJ1tvw+8LOmIXEc5yS8Vh6RT8qmCcg6RtJakDUl99ESO+7g8UkVSN0mfLrPuhFxv7VxvgybiqqS5r8EEFp9aQVK/ZqyzlJb0XRUPAkfl9bcjzcY8X1Lem/S3EVYSQ/p347zD+tCtax0ijbzPO6wPQ/p3a+/QwgrQXufAIY0YjrG9QNIdpIuQniSNFs60/c8q5e8A8yU9CYwijZy/Lmke8E/gp3n69GFJT5EudBpXEsu9wEl5SvF54LGiZR8CO0qaArxHTtrZbEmPkC9ia8G+3w58DngKeAF4PLddahRwuZa8iK2cc0nnTmfkJP4K8MWSOjNYsp9+S5pynZrXmUU6ZzwIGJ77by7pvCykKdUZkqYWX8RGusjsOi2+6vvs/Pso4DKlC646A7eQXrtbgKsknUY6MOoFLHGxV5FJpNdqC+Bc228Ab0jaHng0H6/MBb5OGoEuYvvenEQnS/oE+BPpSvpKcVVS2m/TKtQ7Dbgkv4dWJyXLk6q0W01z+66SS0nvm5mki9iG2f5Y0mWk12oG6QLMSa2ML3RQQ/p3i4S9ipIds2ulJM213aVM+f0UXWDWina72J6bR5eTgIG2/7ls0dYeSXcDh9n+pKR8BDDX9vntElgA0jnwyZNb9RYPISwHkqbYXur7G9rz38hWRXfni97WII0uV7nkDVD671ghhBBaLhJ4GeVG37l80DK2u0zrr+xsj2jvGEIIoVbEd6GHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDIoGHEEIINSgSeAghhFCDmkzgkizpV0XPz5A0YrlGtXQM90uqX5HbrEWS5rawfhdJV0h6UdLTkh6U9NlqbUmqk/SApE6Suks6chni/UFr1y1pZ4ikHdqirSa206L+rdLOMEmbNVGnTfdJ0saS7m2r9kLHM3ZaAwNHTmSrs8YxcORExk5raO+QwnLWnBH4x8BhkjZqzQYkrd6a9dqLpE7tHcMKdDXwLrCt7R2BYUBTr/NxwBjbC4DuQKsTONAmCRwYAiz3BN6GhgFVEzhV9qk1f1O2ZwFvShrY0nVDxzd2WgNnj5lJw5xGDDTMaeTsMTMjia/kmvNBMB+4EvgO8MPiBZK2BK4FNgZmAcfafk3SKFJi6A9MlbQh0Aj0ArYEjgWOAQYAj9seltu7DNgVqANG2z6nWmCSfgx8Kdd/BDjRtiXdD0wHdgPWA46zPSnPHPQAugGbA7+wfZWkQcA5wJtAP0k7A5cB9Xn/v2v7Pkl3ArfbvkHSicDeto+qEt+OwHXAGqSDpS/b/puksXn7awG/tX1lrj8XuAQYDMwmJbhfAFsAp9u+S9Iw4FBgTWAr4CbbPymz7eHAV3K9O0r7UlIP4LPAUbYXAth+CXipWp8DR7E4aY8Etpc0HbgeuDCXDcrbvcT2FZI2BW4lvRarAycDBwF1ed2ni/sxH0RdQ+p/A9faviDHfAnp/fYR8E1gA+BgYB9JP8p9/GK5wPPrvwWwdf79G9sX5mXfJR2cAFxt+zcV2vgVsC/p9fmq7VmS+gGXA2sDL5Leb7PLlQOfy/t1o6RG0t/AOXkf5gMTgDGl+5T74xFgIHCXpBeAH5HeW++QXsd/VXqP5/DHkl6/h8vt28po6BWPtncIK8S01+bwyYKFS5Q1zlvAmaNncPOk19opqhXn1hMHtHcI7aK558AvAY6StH5J+cXADbb7AjeSPsALtgMG2/5efv4pYD/SgcAfgQuAHYE++YMO4Ie264G+pA+vvk3EdbHtXW33JiXxLxYtW8f2HsC3SAcZBX1JyWMA8OOiqczd8vZ3AL4NYLsP8DXgeklrASfkdfYCvgec2kR8J5ESdD/Sh/brufw427vkstPyAQ7AOsD9edkHwM+A/yQl7J8Wtbsb6YO4H3BE6ekFSfsD2+Z6/YBdJO1dEtuOwPQ8km4WSWsAW9t+JRedBfzVdj/bFwD/Bbxne1fSgdg3JW1FSvjjcz/slLd7FtCY1y09COoHdLPdO78G1+XyK4FTc/+cAVxq+xHgLmB4bqts8i7SCziA1DfnSOosaRfSQeVngd1z3P3LrLsOMNX2zsADpMQLcAPw/fx3MLNaue3RwGRSwu1Het8eCuyY6/2syj51tb2P7V8BDwG72+4P3AKcWRRnpff4ZGCvcp0i6QRJkyVNnjVrVhNdGDqa0uTdVHlYOTRrKs72+5JuAE4jjaQLBgCH5ce/I40WC24rSQ5/zKPjmcC/bM8EkPQ0aSp2OvAVSSfkuDYlTSHOqBLavpLOJI1wNgCeJh0cANycY39Q0nqSuubyO203Ao2S7iN9kM8BJtl+OdfZE7gor/+cpFeB7WzPyKP++4BDbb9bJTaAR4EfSvoMadr5b7n8NEmH5sebk5LtO8AnQOE85UzgY9vzcp91L2r3z7bfAZA0Jsc7uWj5/vlnWn7eJW/jwSbibcpGpL6qZH+gr6TD8/P183afAK6V1BkYa3t6E9t5Cdha0kXAOGCCpC7AHsBtkgr11mzFPoyz/THwsaS3gE1I/XeH7Q9hUZ/uxeL+K1hImkkA+D0wJh/UdrX9QC6/PsdYtrxMPO8D/wauljQOuLtK7LcWPf4McGue3VgDeLloWbn3+FjgLSpM3edZoCsB6uvrXSWGmrKqjMwGjpxIw5zGpcq7da1bZfpgVdSSq9B/QxphrVOlTvEf/oclyz7OvxcWPS48Xz2P1M4APpdHIuNIU8xl5RHxpcDheZR2VUn90g8hN1FeHK+orA8p2TZ1DhPbN5GmQhuB8ZL2y9P1g4EBtnciJYlC3PNsF+JZ1E95irv4YKvSPhTHf14evfWzvY3ta0rqPA3sJKkl74FGqrwmebunFm13K9sTbD8I7A00AL+T9I1qG7E9mzRSv580G3I16b06p6jtfra3b0HsBcXvvQWkfq32elcNtZXrLW7Ank9KsLeTzntXu9Cs+D16EWkGqg9wIs1776/FkgfgYSUx/ICe1HVe8vKdus6dGH5Az3aKKKwIzf7wzqPNP5CSeMEjwFfz46NI03qttR7pA+o9SZsABzZRv/CB9XYenR1esnwogKQ9SdO67+XyQyStlaetB5FGh6UeJO0PkrYjnS99XtJuOa7+wBn5oKMiSVsDL+XzrHeRpjbXB2bb/khSL9KUbUv9p6QNJNWRPvRLz2mOB47L/YKkbpI+XVwhT8tOBn6iPKSVtK2kQyptNCfWTvngCdI0/7ol2z05j7SRtJ2kdfK1Em/lc7HXADvn+vMKdYvlCyZXs3078D/AzrbfB16WdESuI0k7lYtD0imSTqm0H2U8CAyRtLakdUhT2n8tU281Fr/PjgQeyu+r2fm0CsDRwAOVykvjza/R+rb/BJxOOn2w1D6VsT7pgAjS9STFKr3HtwOeqtJmqFFD+nfjvMP60K1rHSKNvM87rA9D+ndr79DCctTSq1l/BRR/MJ5GmhodTr6IrbWB2H5S0jTSyPAlmrjQxvYcSVeRpppfYelEPFvSI+SL2IrKJ5FG91sA59p+IyfpYpcCl+ep6/mkq4YhjfKPzet8j7Tv+5FGQNi+vKSdocDXJc0D/kk6j/0hcJKkGcDzwGPV9rOCh0inLLYhXcRWPH2O7QmStgcezbl5LvB10hRqseNJr+nfJX1EmlkY3sS2J5CmnP+PdHpjvqQngVHAb0lT/VPzQcEs0gHGIGB47oe5QGEEfiUwQ9LUkvPg3YDrimYHzs6/jwIuyxd2dSad+30y/75K0mmkBNuLFlyoZXuq0oWXk3LR1bZLp88hvXY7SpoCvEc+SCQl0MslrU167x7bRPmoXN5IOiC8Mx8UiXSNCGX2qdQI0lR9A+k9VHwwudR7PJfvm8vDSmhI/26RsFcxWjxju/JQugr9jNLEpnSF7lzb57dHXG1B6Sr0etstGWG25fb7k67KP7o9tt8cku4GDrP9SXvHsqJVe49LehA4JM+kVFRfX+/JkydXqxJCWIEkTckXeC8hvokttEgemd6nDvz/8ra/uCom72okbQz8uqnkHUKoHSvlCDyE0HqSZgGvNlFtI+DtFRDOsogYl11Hjw86foxtEd+WtjcuLYwEHkJoMUmTy03pdSQR47Lr6PFBx49xecYXU+ghhBBCDYoEHkIIIdSgSOAhhNa4sr0DaIaIcdl19Pig48e43OKLc+AhhBBCDYoReAghhFCDIoGHEEIINSgSeAhhCZI+L+l5SX+XdFaZ5cMlTc8/T0lakL+bf3NJ90l6VtLTkv67I8VXtLyTpGn5G/uWi2WJUVJXSaMlPZf7crncTmwZY/xOfo2fknRz0f0RVmR860v6o6QncyzHNnfd9o6xzf5WbMdP/MRP/GAboBPwIrA16TalTwI7VKn/JWBifrwp6cYzkG7E8kK1dVd0fEVl3wVuAu7uaH2Yn18PHJ8fr0G6LW2HiZF0r4KXgbr8/A/AsBUdH/AD4H/z442Bd3PdFu1bO8XYJn8rMQIPIRTbDfi77Zecvo72FqDiHeqArwE3A9h+0/bU/PgD4FnSh32HiA9A0meAg0i3qF1eWh2jpPVIt969BsD2J7bndKQYs9WBOkmrA2sDb5Rda/nGZ2DdfOOkLqTkOL+Z67ZrjG31txIJPIRQrBvwj6Lnr1PhgyXfZe3zpHuZly7rTrrt7uMdLL7fAGcCC9s4rmLLEuPWpLv4XZen+a9WusVth4nRdgNwPvAa8Cbpds0T2iG+i4HtSQcPM4H/tr2wmeu2d4yLLMvfSiTwEEIxlSmr9L+mXwIetv3uEg2ke5zfDpzudB/3DhGfpC+S7ks/pY1jKrUsfbg6sDNwme3+pFvYLo9zuMvSj58ijTS3AjYD1pH09XaI7wBgeo6hH3BxnsFoyb4ti2WJMTWwjH8rkcBDCMVeBzYvev4ZKk+PfpUlp1WR1Jn0gXSj7TEdLL6BwMGSXiFNd+4n6fcdLMbXgddtF0Zjo0kJva0tS4yDgZdtz7I9DxgD7NEO8R0LjHHyd9J5+V7NXLe9Y2ybv5W2PrEfP/ETP7X7QxoBvkQaXRUuzNmxTL31Sefz1ikqE3AD8JuOGF/J8kEsv4vYlilG4K9Az/x4BPDLjhQj8FngadK5b5Euujt1RccHXAaMyI83ARpId/5q1r61c4xt8reyepXcHkJYxdieL+kUYDzpKttrbT8t6aS8/PJc9VBggu0Pi1YfCBwNzJQ0PZf9wPafOkh8K0QbxHgqcKOkNUgJ4lja2LLEaPtxSaOBqaSLxqbRxl8X2sz4zgVGSZpJSojft/02QLl12zK+ZY1R0p60wd9KfJVqCCGEUIPiHHgIIYRQgyKBhxBCCDUoEngIIYRQgyKBhxBCCDUoEngIIYRQgyKBhxBCGfnuW4U7cd2Wv1K0tW2NknR4fny1pB2q1B0kaY+i5ydJ+kZrt12m/e9I+rek9duqzdA+IoGHEEJ5jbb72e4NfAKcVLxQUqfWNGr7eNvPVKkyiKJvNrN9ue0bWrOtCr4GPEH6H+/lQknkl+UsOjiEEJr2V2CbPDq+T9JNpC/h6CTpl5KekDRD0omwKIFdLOkZSeOATxcaknS/pPr8+POSpub7Rf8l39jiJOA7efS/l6QRks7I9ftJeixv6478veSFNv9X0iRJL0jaq9xOSOpBuivWj0iJvFDeRdJ1kmbmtr9cLr5ctiie/PwpSd3zz7OSLiV9ycvmki6TNFnpntc/KVpnV0mP5HYnSVpX0l8l9Suq87Ckvq1/yVZ+8U1sIYRQhdItMw8E7s1FuwG9bb8s6QTS3bh2lbQm8LCkCaS7S/UE+pC+QvMZ4NqSdjcGrgL2zm1tYPtdSZcDc22fn+t9rmi1G0hfW/qApJ8C5wCn52Wr295N0hdy+eAyu1O4LehfgZ6SPm37LeB/8n70ydv8VLn4mtFdPYFjbX8rt/PDvE+dgL/khPwccCsw1PYTSjf3aCTd4nUYcLqk7YA1bc9oxjZXWTECDyGE8ury11xOJt0685pcPsn2y/nx/sA3cr3HgQ2BbUn39L7Z9gLbbwATy7S/O/BgoS2X3NWtVD5n3dX2A7no+rydgsINMaYA3Ss081XgFqdbWo4Bjsjlg4FLCpVsz25pfNmrth8rev4VSVNJX7e6I7ADKcm/afuJ3O77tucDtwFfVLrJx3HAqGZsb5UWI/AQQiiv0Xa/4gJJkG7xuaiINCIeX1LvCzR9C0s1o05LfJx/L6DMZ3se/W4L/DnvR+G71i+pEEul+Oaz5OBvraLHi/pG0lbAGcCutmdLGpXrlm3X9keS/ky6VelXgPoK+xmyGIGHEELrjQdOzqNGJG0naR3gQeCr+Rz5psC+ZdZ9FNgnJzqKpqg/ANYtrWz7PWB20fnto4EHSutV8TXSnbG655/NgG6StgQmAKcUKuZz65Xie4V8i1NJO5PuxlXOeqSE/p6kTUinISBNoW8madfcxrr5NAWkafQLgSeaOeJfpcUIPIQQWu9q0nT1VKVh7SxgCHAHsB8wE3iBMonW9qx8Dn1MvmL7LeA/gT8CoyUdQrozWbFjgMuV/qWtpXcq+yqLk2jBHbn8Z8Alkp4ijeB/YntMhfhuZ/Fpgyfy/i3F9pOSppFuPfoS8HAu/0TSUOAiSXWk89+DSef9p0h6H7iuBfu1yoq7kYUQQugQJG0G3A/0yufpQxUxhR5CCKHdKX1ZzePADyN5N0+MwEMIIYQaFCPwEEIIoQZFAg8hhBBqUCTwEEIIoQZFAg8hhBBqUCTwEEIIoQb9fx5nrO3XmdBrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "labels = [\n",
    "    \"Normal approx. sample CI (test set, no bootstrap)\",\n",
    "    \"Bootstrapping the test set, percentile method\",\n",
    "    \"Bootstrap, sample CI\",\n",
    "    \"Bootstrap, percentile method\",\n",
    "    \".632 Bootstrap, percentile method\",  \n",
    "    \".632+ Bootstrap, percentile method\",    \n",
    "]\n",
    "\n",
    "mean = [\n",
    "    acc_test_na,\n",
    "    mean_test_accuracies,\n",
    "    bootstrap_train_mean,\n",
    "    bootstrap_oob_mean,\n",
    "    bootstrap_632_mean,\n",
    "    bootstrap_632plus_mean\n",
    "]\n",
    "\n",
    "lower_error = [\n",
    "    test_na_lower,\n",
    "    bootstrap_lower_test,\n",
    "    bootstrap_na_lower,\n",
    "    bootstrap_oob_percentile_lower,\n",
    "    bootstrap_632_percentile_lower,\n",
    "    bootstrap_632plus_percentile_lower\n",
    "]\n",
    "\n",
    "upper_error = [\n",
    "    test_na_upper,\n",
    "    bootstrap_upper_test,\n",
    "    bootstrap_na_upper,\n",
    "    bootstrap_oob_percentile_upper,\n",
    "    bootstrap_632_percentile_upper,\n",
    "    bootstrap_632plus_percentile_upper    \n",
    "]\n",
    "\n",
    "\n",
    "asymmetric_error = [np.array(mean) - np.array(lower_error), \n",
    "                    np.array(upper_error) - np.array(mean)]\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(7, 3))\n",
    "ax.errorbar(mean, np.arange(len(mean)), xerr=asymmetric_error, fmt='o')\n",
    "#ax.set_xlim([0.86, 1.0])\n",
    "ax.set_yticks(np.arange(len(mean)))\n",
    "ax.set_yticklabels(labels)\n",
    "ax.set_xlabel('Prediction Accuracy')\n",
    "ax.set_title('95% confidence intervals')\n",
    "plt.tight_layout()\n",
    "plt.savefig('CI_mnist.pdf')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}