Introduction+to+Bootstrap.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:3595dcf94f9e25241e065eadbb79541d7a4c601c2c99b1e27e61d58950062e1f"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "from scipy import stats"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## We are going to be exploring bootstrapping using the empirical distribution. \n",
      "\n",
      "As a reminder:\n",
      "$$ \\textrm{Var}_F(T_n) \\approx \\textrm{Var}_{\\hat{F_n}}(T_n) \\approx v_{boot}$$\n",
      "There are two approximations going on here.  To compute $v_{boot}$:\n",
      "\n",
      "For i in B:\n",
      "1. Draw $X_1^*, X_2^*, ..., X_n^* \\sim \\hat{F_n}$\n",
      "2. Compute $T_{n,i}^* = g(X_1^*, X_2^*, ..., X_n^*)$\n",
      "\n",
      "Now that you have $T_{n,1}^*,T_{n,2}^*,...,T_{n,B}^*$, compute\n",
      "$$v_{boot} = \\frac{1}{B} \\sum_{j=1}^B \\left( T_{n,j}^* - \\frac{1}{B} \\sum_{k=1}^B T_{n,k}^* \\right)^2$$\n",
      "\n",
      "We have 3 methods for constructing confidence intervals\n",
      "1. *The Normal Interval* is simple, but only correct to use if $T_n$ is approximately normal. Then $\\hat{se} = \\sqrt{v_{boot}}$ and the confidence interval is $$ C_n = (T_n - z_{\\alpha/2}\\hat{se}, T_n + z_{\\alpha/2}\\hat{se})$$\n",
      "2. *Pivot Interval*.\n",
      "   \n",
      "   Let $\\theta = T(F)$, $\\hat{\\theta} = T(\\hat{F_n})$\n",
      "   \n",
      "   Let $\\hat{\\theta}_{n,1}^*, \\hat{\\theta}_{n,2}^*,..., \\hat{\\theta}_{n,B}^*$ be the bootstrap replications of $\\hat{\\theta}$.  \n",
      "   \n",
      "   Let $\\hat{\\theta}^*_{\\beta}$ denote the $\\beta^{\\textrm{th}}$ sample quartile. \n",
      "   $$ C_n = (2\\hat{\\theta_n} - \\hat{\\theta}^*_{1-\\alpha/2}, 2\\hat{\\theta_n} - \\hat{\\theta}^*_{\\alpha/2}) $$\n",
      "3. *Percentile Interval* $$ C_n = (\\hat{\\theta}^*_{\\alpha/2}, \\hat{\\theta}^*_{1-\\alpha/2} )$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N = 300\n",
      "B = 1000\n",
      "\n",
      "x = normal(size=N)\n",
      "boot_means = [random.choice(x, size=N, replace=True).mean() for i in range(B)]\n",
      "\n",
      "sampled_means = [normal(size=N).mean() for i in range(B)]\n",
      "\n",
      "boot_density = stats.kde.gaussian_kde(boot_means)\n",
      "sampled_density = stats.kde.gaussian_kde(sampled_means)\n",
      "\n",
      "figsize(18, 6)\n",
      "y = arange(-0.8, 0.8, .01)\n",
      "plot(y, boot_density(y), 'b', label=\"Bootstrap\")\n",
      "plot(y, sampled_density(y), 'r', label=\"Sampling\")\n",
      "plot(y, stats.norm.pdf(y,0,1.0/sqrt(N)), 'k', label=\"Theoretical\")\n",
      "legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "<matplotlib.legend.Legend at 0x10af5ff90>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAABA0AAAFwCAYAAADT+YJVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl01dW99/HPyQABMhBICIQpBGSGhJCcDBAIg0UUUMSx\niFZtvb1tbbW9t8+tE7FXq7UqKtqW+tyl1Vq8rY8DUkVFCVPmgYRAlEFCgEQQgQyEEJKc549DIoEE\nM/x2pvN+rZXVkHPy3ZucriXnk/39bpvD4XAIAAAAAADgIm6dvQEAAAAAANA1ERoAAAAAAIAmERoA\nAAAAAIAmERoAAAAAAIAmERoAAAAAAIAmERoAAAAAAIAmtTg0eOKJJzRp0iRNmTJF3//+93X27FmT\n+wIAAAAAAJ2sRaFBYWGhXn75ZWVnZ2vnzp2qra3Vm2++aXpvAAAAAACgE3m05Em+vr7y9PRUZWWl\n3N3dVVlZqaFDh5reGwAAAAAA6EQtOmkwYMAA/epXv9KIESMUHBys/v37a/78+ab3BgAAAAAAOlGL\nQoP9+/frueeeU2FhoYqLi1VRUaE33njD9N4AAAAAAEAnalF7QmZmpuLi4jRw4EBJ0vXXX6/k5GQt\nX7684TljxozR/v37zewSAAAAAAC02ejRo7Vv375Wf1+LQoPx48frv//7v3XmzBl5eXlp48aNstvt\njZ6zf/9+ORyOVm8A3VtiYqISExM7exvoYLzurofX3DXxurseXnPXxOvuenjNXZPNZmvT97WoPSEs\nLEy33367IiMjNXXqVEnSPffc06YFAQAAAABA99CikwaS9Otf/1q//vWvTe4FAAAAAAB0IS06aQA0\nJyEhobO3gE7A6+56eM1dE6+76+E1d0287q6H1xytYXNYNIjAZrMx0wAAAAAAgC6ore/ZW9yeAAAA\nAABwLQMGDNDJkyc7extoBX9/f504ccKyepw0AAAAAAA0ifd53U9zr1lbX0tmGgAAAAAAgCYRGgAA\nAAAAgCYRGgAAAAAAgCYRGgAAAAAAgCYRGgAAAAAAupWQkBD17dtXPj4+GjBggBYtWqTDhw+3q2Zh\nYaHc3NxUV1fXoue7ubnpyy+/bNea3QGhAQAAAACgW7HZbFq/fr3Ky8tVUlKioKAg3XvvvZbUbs0N\nA5d7bk1NjRXb6XSEBgAAAACAbqt3795atmyZdu/eLUkqLS3V7bffrkGDBikkJESPP/54w5t7h8Oh\nxx57TCEhIQoKCtIdd9yhsrIySdKsWbMkSf3795ePj4/S0tK0b98+zZ49W/3791dgYKBuvfXWRs8N\nCwuTj4+P/vnPfyopKUnDhg3TU089pSFDhujuu+/WqVOntGjRIg0aNEgDBgzQ4sWLdeTIkYa9JyQk\n6De/+Y2io6Pl5+en6667TidPnuywn11LEBoAAAAAALqd+iCgsrJS//u//6vY2FhJ0r333qvy8nId\nOHBAmzdv1muvvaZXXnlFkvTKK6/or3/9q5KSkvTll1+qoqJCP/vZzyRJW7duleQMHcrLyxUdHa2H\nH35YV111lU6dOqUjR440nGbYsmWLJCkvL0/l5eW68cYbJUlHjx7VyZMnVVRUpDVr1qiurk533323\nioqKVFRUpD59+jSsV+/111/XK6+8opKSEnl4eOjnP/+54Z9c69gcrTl7cblCNlurjnEAAAAAALq2\nrvo+LyQkRN988408PDx0+vRpDRo0SBs2bNCECRPUt29f5ebmavz48ZKkv/zlL1q7dq02bdqkefPm\n6cYbb9SPf/xjSdKePXs0efJkVVVVqaioSKGhoaqpqZGbm/P363fccYe8vLz0yCOPaOjQoY324Obm\npn379ik0NFSSlJSUpAULFqi8vFy9evVqct87duzQ3LlzdeLECUnSnDlzFBsbq9/97neSpIKCAoWH\nh6uqqko2m61NP5vmXrO2vpacNAAAAAAAtInNZs1H69e16b333tPJkyd19uxZrV69WrNnz9aRI0d0\n7tw5jRw5suG5I0aMaGgJKCkpueSxmpoaHT16tMl1nnrqKTkcDtntdk2ePLnhxEJzAgMDGwUGlZWV\n+rd/+zeFhITIz89Ps2fPVmlpaaM378OHD2+0n3Pnzun48eOt+4EYRGgAAAAAAGgTh8Oaj/aw2Wxa\nunSp3N3dlZqaKk9PTxUWFjY8XlRUpGHDhkmSgoODL3nMw8NDQUFBTf5mPygoSH/5y1905MgRrVmz\nRj/5yU8ue2PCxTWeeeYZ7dmzR+np6SotLdXmzZvlcDgahQZFRUWNPvf09FRAQEBrfwzGEBoAAAAA\nALqdC4cbvvfeezp16pSmTJmim266SQ8++KAqKip08OBBrVq1Srfddpsk6dZbb9WqVatUWFioiooK\nPfDAA7rlllvk5uamwMBAubm5af/+/Q1r/POf/2y4yrF///6y2WwNrQtBQUGNntuUiooK9enTR35+\nfjpx4oQeffTRS/4Of/vb31RQUKDKyko98sgjuvHGG9vcmmACoQEAAAAAoNtZvHixfHx85Ofnp4cf\nflh//etfNWHCBK1evVr9+vVTaGio4uPjtXz5ct15552SpLvuuksrVqzQrFmzFBoaqr59+2r16tWS\npL59++rBBx/UjBkzNGDAAKWlpSkzM1MxMTHy8fHRtddeqxdeeEEhISGSpMTERN1xxx3y9/fXW2+9\nJZvNdsmb/fvuu09nzpxRQECA4uLitHDhwkbPsdlsWrFihX7wgx9oyJAhqq6u1gsvvNAxP8AWYhAi\nAAAAAKBJvM8za86cOVqxYoXuuusuy2oyCBEAAAAAgB6iq4cyhAYAAAAAAHSSrjS/oCm0JwAAAAAA\nmsT7vO6H9gQAAAAAANAhCA0AAAAAAECTCA0AAAAAAECTCA0AAAAAAECTCA0AAAAAAECTCA0AAAAA\nAGiDxMRErVixQpJUVFQkHx+fHnfbBKEBAAAAAKDb2bZtm+Li4tS/f38NHDhQM2fOVGZmZofuwWaz\nNXw+YsQIlZeXN/paT+DR2RsAAAAAAKA1ysrKtGjRIq1Zs0Y33XSTzp49q61bt6p3794duo+edqqg\nKZw0AAAAAAB0K3v27JHNZtPNN98sm80mLy8vXXnllZoyZYr279+vuXPnKiAgQIGBgbrttttUWlra\n8L0hISF6+umnNXXqVPn4+Ojuu+/W0aNHtXDhQvn5+enKK6/UqVOnJEmFhYVyc3PTyy+/rKFDhyo4\nOFjPPPNMk3uqf25dXZ0kKSEhQY888ohmzpwpX19fLViwQN98803D81977TWNHDlSAQEBeuyxxxQS\nEqJPP/3U4E+tbQgNAAAAAADdyrhx4+Tu7q4f/OAH2rBhg06ePNno8QcffFAlJSUqKCjQoUOHlJiY\n2PCYzWbT22+/rU8//VRffPGF1q9fr4ULF+rJJ5/UsWPHVFdXpxdeeKFRvaSkJO3bt08ff/yxfv/7\n37f4zf3atWv16quv6tixY6qurtbTTz8tSdq9e7d++tOfau3atSopKVFpaamKi4u7ZGsDoQEAAAAA\noFvx8fHRtm3bZLPZ9KMf/UiDBg3Stddeq2PHjmn06NGaN2+ePD09FRAQoPvvv1+bN29u9P333nuv\nAgMDFRwcrPj4eMXGxiosLEy9e/fW0qVLlZOT0+j5K1euVJ8+fTR58mTdeeedWrt27Xfu0Waz6c47\n79SYMWPk5eWlm266STt27JAkvfXWW1qyZIni4uLk6emp3/72t10yMJAIDQAAAAAAbWWzWfPRBuPH\nj9crr7yiQ4cOKT8/X8XFxbrvvvt07Ngx3XLLLRo2bJj8/Py0YsWKRm0BkhQUFNTweZ8+fRr92cvL\nSxUVFY2eP3z48IbPR4wYoeLi4hbtcfDgwY3Wqa9bXFysYcOGNXps4MCBLarZ0VocGnzxxReaNm1a\nw4efn98lRzYAAAAAAC7E4bDmo53GjRunO+64Q/n5+XrggQfk5uam/Px8lZaW6vXXX2+YM9D8X+Py\neygqKmr0+dChQ9u13+DgYB0+fLjhz2fOnLkk2OgqWhwajBs3Tjk5OcrJyVFWVpb69u2rpUuXmtwb\nAAAAAACX+OKLL/Tss8/qyJEjkqRDhw5p7dq1io2NVXl5uby9veXr66sjR47oD3/4Q7vXe+yxx3Tm\nzBnt2rVLr776qm6++eYWfV9zYcSyZcv0/vvvKyUlRdXV1UpMTOyyNzG0qT1h48aNGj16dKMjGgAA\nAAAAdAQfHx+lpaUpOjpa3t7eio2N1dSpU/XMM89o5cqVys7Olp+fnxYvXqxly5Z957yACx+32WyX\nPH/27NkaM2aM5s+fr//8z//U/Pnzm3zuxd/XXN1JkyZp9erVuuWWWxQcHCwfHx8NGjSow6+MbAmb\now1xxl133aXIyEj95Cc/+baQzdZlkxEAANBK+/ZJISGSh0dn7wQA0Ilc/X1eYWGhQkNDVVNTIzc3\ncyMBKyoq5O/vr3379mnkyJHtqtXca9bW17LVf+vq6mq9//77uvHGG1u9GAAA6AZOnZLCw6X586WS\nks7eDQAAPdL777+vyspKnT59Wv/xH/+hqVOntjswMKHVvz748MMPNX36dAUGBl7y2IV3XyYkJCgh\nIaE9ewMAAJ3h5ZelxYul8eOl6dOlN96Q5szp7F0BANApTF2FuG7dOt1+++1yOByKiorSm2++aWn9\npKQkJSUltbtOq9sTbrnlFi1cuFB33HFH40IufmwFAIAeoaZGCg2V3nnHGRh88ol0++3ST38qPfCA\nZPBoJgCg6+F9XvdjdXtCq0KD06dPa+TIkTpw4IB8fHws2QAAAOhC/vEPafVqaevWb7925Ih0yy2S\nt7f0+utSQEDn7Q8A0KF4n9f9dOpMg379+un48eOXBAYAAKBjVVVJn38uffih9NJL0n/8h3TzzdK2\nbe0svGqV0hct0ty5c/XnP//Z+Y+LoUOlzz6Tpk6VIiKk/HxL/g4AAKDra9PtCU0WIoECAMC4nBzn\nL/0LC6URI6RRo5zdBKNGSb17S08+KaWmOi8+aK1z27bpvxcv1ppevfTII4/o1Vdf1cCBA/U///M/\nGjp0qPNJv/+9lJsr/f3vVv61AABdFO/zuh+rTxpwjxIAAN3EJ59Iy5c7uwduuEFyd7/0OTabdO21\nUnKy1K9fy2vv2rVLty9ZoiHBwdqxcaOGDBmie+65R0888YQiIiK0atUq3XrrrbLdfbc0ZoxUXi5x\n8hAAgB6PkwYAAHQDr7/ubEF46y0pPr755zkc0p13SpWV0v/+rzNEuJza2lo999xzevJ3v9MTVVW6\nu7hYNj+/Rs/JysrS7bffrkmTJumPf/yjAu66y5la3H67BX8zAEBXxvu87qdTZxoAAICO5XBITzwh\nPfSQtGnT5QMDyRkS/PnP0sGDzlaFy6moqNDcuXP13nvvKf366/XDf/u3SwIDSZo+fbqysrI0YsQI\nhYWFKT0mxpliAADQRSUmJmrFihWdvY0mvfHGG1qwYEG767i5uenLL7+0YEeXR3sCAABdVG2t9POf\nO4cbJic75xG2hJeX9Pbbkt3unF14zTVNP+/FF1/UwIED9c9XXpF7aKiUlXWZml56+umnNWXKFP38\npZeUsmePbMXFUnBwG/5mAAC0j7e3t2znj9OdPn1aXl5ecj/ft7dmzZqGxzpbYWGhQkNDVVNTI7fz\n1xYvX75cy5cv7+SdtRwnDQAA6ILOnnV2AHzxhbRlS8sDg3pDhzpbGe6803nLwsXKy8v17LPP6vHH\nH5f7a69Jc+a0aHriihUrVH76tD6y26W1a1u3KQAALFJRUaHy8nKVl5dr5MiRWr9+fcOfv//973dY\nS0VNTU2LntedWzwIDQAA6IJWr3bOJfjgA6mJjoEWiY11tjZce6106lTjx1566SXNnz9fE8aOlZ5/\nXrr//hbVdHNz08qVK5V4+LActCgAALoom82m6upq3XHHHfL19dXkyZOVdcGJuuLiYi1btkyDBg1S\naGioVq9e3fDY2bNndd9992no0KEaOnSo7r//flVXV0uSkpKSNGzYMD311FMaMmSI7r77bjkcDj35\n5JMaM2aMAgICdPPNN+vkyZOSpFmzZkmS+vfvL19fX6WmpurVV19V/AX9hrt27dKVV16pgQMHavDg\nwXriiSckSenp6YqNjZW/v7+Cg4N177336ty5c8Z/dhcjNAAAoIupqJCeflp65hmpV6/21br7bmnu\nXOk3v/n2a/WnDB5++GHpX/+S/P2luLgW17zhhhtUbrPpoyNHpPz89m0QAAADHA6H1q1bp1tvvVWl\npaVasmSJfvazn0mS6urqtHjxYk2bNk3FxcX69NNP9dxzz+njjz+WJD3++ONKT09Xbm6ucnNzlZ6e\nrscee6yh9tGjR3Xy5EkVFRVpzZo1euGFF7Ru3Tpt2bJFJSUl8vf3109/+lNJ0tatWyVJpaWlKisr\nU0xMTKN9lpeXa/78+br66qtVUlKiffv2ad68eZIkDw8PPf/88/rmm2+UkpKiTz/9VH/84x+N/+wu\nRmgAAEAX89JLUkKCNHmyNfX+8z+dMw7q6px/fvHFF52nDCZMkFatcp4yaEXvZ/1pg5Wenpw2AAB0\nWfHx8brqqqtks9l02223KTc3V5KUkZGh48eP66GHHpKHh4dGjRqlH/7wh3rzzTclOQcVPvLIIwoI\nCFBAQIBWrlyp1y/4752bm5seffRReXp6ysvLS2vWrNFjjz2m4OBgeXp6auXKlXrrrbdUV1f3nW0J\n69evV3BwsO6//3716tVL3t7estvtkqSIiAjZ7Xa5ublp5MiRuueee7R582ZDP63mMQgRAIAupLxc\nevZZ500JVgkNlQIDpfR0adKkcq1atcr5j44vv5QKCqQbb2x1zRtuuEG/feABbfif/9HCJ56Q3Pg9\nBAC4IqsGDpro+Q8KCmr4vG/fvqqqqlJdXZ0OHjyo4uJi+fv7NzxeW1vb0EpQUlKikSNHNjw2YsQI\nFRcXN/w5MDBQvS44ClhYWKilS5c2DDqUnKcEjh49+p17PHTokEJDQ5t8bM+ePfrlL3+prKwsVVZW\nqqamRpGRkS34m1uL/8IDANCFrF4tzZ8vTZxobd0lS6R165ynDK688krnKYOUFOcdjp6era7n5uam\nlU88ocQzZ+RISrJ2swCAbsPhcFjyYbXLhRnDhw/XqFGjdPLkyYaPsrIyrV+/XpIUHByswsLChucX\nFRUp+ILbgi6uPWLECG3YsKFRvcrKSg0ZMuQ7Q5URI0Y0e23iv//7v2vixInat2+fSktL9fjjj6uu\n/thgByI0AACgiygrk557TnrkEetrL1kivfOO85TBww8/7PxiWprzXsY2WrZsmU77+GjDk09atEsA\nAKxxuSDCbrfLx8dHTz31lM6cOaPa2lrl5+crMzNTknTrrbfqscce0/Hjx3X8+HH99re/1YoVK5qt\n9+Mf/1gPPPCAioqKJElff/211q1bJ8l5KsHNzU379+9v8nuvueYalZSU6Pnnn9fZs2dVXl6u9PR0\nSc4bInx8fNS3b199/vnn+tOf/tSmn0V7ERoAANBFPP+8dNVV0rhx1te226XDh19UbOyVGj9+vPOL\n6elSdHSba7q5uWnlo48q8bPP5KistGinAAC0n81mu+S3/PV/dnd31/r167Vjxw6FhoYqMDBQ99xz\nj8rKyiRJDz30kCIjIzV16lRNnTpVkZGReuihhy6pU+8Xv/iFlixZou9973vy9fVVbGxswxv/vn37\n6sEHH9SMGTM0YMAApaWlNdqbj4+PPvnkE73//vsaMmSIxo4dq6TzJ/iefvpp/f3vf5evr6/uuece\n3XLLLY3Wtqo15LvYHBadBbHZbN367kkAADrTqVPSFVdIycnO/7VaWVmZgoLG6L77tuiJJ8ZLZ89K\nAwZIR49K3t5trltXV6epvr566t57dfX5K6IAAD0H7/O6n+Zes7a+lpw0AACgC3juOemaa8wEBpJz\nloHdfqXS08+fMsjNlcaMaVdgIJ0/bXDHHUr885/5RyUAAD0QJw0AAOhkJ086w4K0NGn0aOvrl5WV\nacyYMdqwYYsSEsbr4EHJ/40XncHByy+3u37dqVMKGzhQv//b33T1rbdasGMAQFfB+7zuh5MGAAD0\nMM8+K117rZnAQJJeffVVzZ07VxER4zV7trRhg5wJRTvmGVzIrX9/PRQTo6cTEy2pBwAAug5CAwAA\nOtE330h//KN0wXwly7311lsNU5/rr15s7xDEi137y18qe/9+ff3115bVBAAAnY/QAACATvSnP0nX\nXy+NGmWm/tGjR5WXl6d58+ZJkhYtklI/PClHcbE0caJl63gtWKDvSVr3j39YVhMAAHQ+QgMAADrR\ne+9Jy5ebq79u3TpdddVV8vLykiQNGSItCsrQqdHTJXd36xby9tb1V1yhd/76V+tqAgCATkdoAABA\nJzl6VNq7V5oxw9wab7/9tq6//vpGX1s2LE3ZHnbL17p62TJtyc1tuOcaAAB0f4QGAAB0kg0bpPnz\nJU9PM/VLS0u1fft2LVy4sNHXI2rS9f8ORcvqYdi+ixYpvndvffDBB9YWBgB0Gn9/f9lsNj660Ye/\nv7+l/x8gNAAAoJN88IF09dXm6v/rX//S7Nmz5ePj8+0XHQ75FKQpx9OuXbssXjAyUtfX1OidtWst\nLgwA6CwnTpyQw+Hgoxt9nDhxwtL/DxAaAADQCWpqpE8+kS46BGCpploTVFgom4eHopYO0/vvW7yg\nh4eWzJ6tjzZuVFVVlcXFAQBAZyA0AACgEyQnO29MGDLETP0zZ87ok08+0eLFixs/cP6qxcVLbM6r\nFy0WeM01Cvf11SeffGJ9cQAA0OEIDQAA6ASmWxM+/vhjTZ8+XQEBAY0fSEuT7HbNni0VFEhffWXx\nwvPn6/qqKr399tsWFwYAAJ2B0AAAgE5gOjRosjVBajhp0KuXtGCB9K9/WbzwuHG6rndvvf/uu6qp\nqbG4OAAA6GiEBgAAdLBDh6TiYslu/a2HkqRz585p/fr1uu666y5+QNqxQ4qMlCQtWSLrWxRsNo1Y\nuFCjfH21ZcsWi4sDAICORmgAAEAH+/BD52/53d3N1N+8ebPGjBmjYcOGNX4gP18aOVLy9ZXk3MOm\nTVJtrcUbmD9fS729aVEAAKAHIDQAAKCDdVprwvl5BvUCApyDGHfvtngD8+bp+kOH9O6776qurs7i\n4gAAoCMRGgAA0IHOnnX+dn/BAjP16+rq9O6772rp0qWXPpiWJkVHN/pSTIyUmmrxJgYP1viRI+Xr\n6amMjAyLiwMAgI5EaAAAQAfaulWaONH5W34T0tLSNGDAAI0dO/bSB88PQbyQkdBAcrYojBxJiwIA\nAN0coQEAAB2o01oTysqkgwelyZMbfdlkaHB9RYXefvttORwOAwsAAICO0OLQ4NSpU7rhhhs0YcIE\nTZw4UalG/oUBAEDP9sEH0jXXmKntcDj0zjvvNN2akJkphYdLnp6NvjxlijNLOHXK4s3MmqWIggJV\nnz2rXbt2WVwcAAB0lBaHBr/4xS909dVXq6CgQHl5eZowYYLJfQEA0OPs3y+Vljrfu5uwc+dO1dbW\nKrypBS4agljPw0OaPl2yfPSAj49s06dr6fTptCgAANCNtSg0KC0t1datW3XXXXdJkjw8POTn52d0\nYwAA9DQffCAtXCi5GWoOrG9NsNlslz7YxDyDekZbFHr31jvvvGOgOAAA6Agt+mfLgQMHFBgYqDvv\nvFMRERH60Y9+pMrKStN7AwCgRzE9z2DDhg1atGjRpQ84HM2eNJDMhgYzCgp08OBBlZSUGFgAAACY\n1qLQoKamRtnZ2frJT36i7Oxs9evXT08++aTpvQEA0GNUVkrbtklXXmmmfnl5ufLz8xUbG3vpg4cP\nSzU1UkhIk98bHe0MDSyfVxgVJffCQsXb7dq8ebPFxQEAQEfwaMmThg0bpmHDhikqKkqSdMMNNzQZ\nGiQmJjZ8npCQoISEBEs2CQBAd7dpk3N2gKnuvu3btysyMlJeXl6XPljfmtBU24Kk4GDJ21vat0+6\n4goLN+XpKc2erYT+/ZWUlKRbbrnFwuIAAOBykpKSlJSU1O46LQoNBg8erOHDh2vPnj0aO3asNm7c\nqEmTJl3yvAtDAwAA8C3TrQlJSUnNh/Xp6dL54L859S0KloYGkjR/vhI++0xr0tMtLgwAAC7n4l/k\nP/roo22q0+JRTKtXr9by5csVFhamvLw8PfDAA21aEAAAV7Rxo7Rggbn6mzdvbj40yMxscWhgufnz\nNTU7W0ePHtVXX31lYAEAAGCSzeGwpoPRZrPJolIAAPQox49Lo0dLJ05I7u7W16+oqNDgwYN1/Pjx\nS9sT6uqkAQOkvXulwMBma6SkSD/7mZSVZfHmHA5pyBBdO2WKvv/DH+rmm2+2eAEAANASbX3PbujS\nJwAAUC811XlxgYnAQHLOM5g+fXrT8wz27pX6979sYCBJ06ZJn3/uHNhoKZtNio9XQmCgJX2VAACg\nYxEaAABgWEqK1NSlBla57DyDjIzvbE2QJC8vafJkAycNJCk+XrPPnCE0AACgGyI0AADAsE4NDVow\nz6CesbkG8fEK271bX331FXMNAADoZggNAAAwqKbG+cv+mBgz9SsqKrRz507FNLdAC08aSAZDg6lT\n5X70qOLtdm3ZssXAAgAAwBRCAwAADNq5Uxo+XPL3N1M/OTlZ06dPV58+fS59sKZGys2VIiJaVCsm\nxnkqwvK5xu7uUmysEoYMoUUBAIBuhtAAAACDOqI1Yfbs2U0/uGuXNGyY5OfXolohIVJtrXT4sHX7\nazBrlhLOnSM0AACgmyE0AADAoO4yz0ByXnRgdK7B55+rpKRER48eNbAAAAAwgdAAAACDkpOluDgz\ntU+fPq28vDxL5hnUMxYaREXJ/YsvNDMmRps3bzawAAAAMIHQAAAAQ44dk06ckMaPN1M/OTlZERER\n6tu3b9NPyMiQIiNbVdNYaNC7txQRoYThwwkNAADoRggNAAAwJCVFio6W3Az91/ay8wyqqqSCAik8\nvFU1IyOlHTuk6moLNnix+Hgl1NYy1wAAgG6E0AAAAEM6dZ5BXp40dqzU3CmEZvj4SGPGOC9dsNys\nWQrfu1fC6XmbAAAgAElEQVRHjhzRsWPHDCwAAACsRmgAAIAhKSlm5xnk5uYqtrlUog3zDOoZa1GI\njZV7drbiZ8ygRQEAgG6C0AAAAAPOnZOyspztCSYkJydr2rRpls4zqGcsNPD1lcaNU8KoUbQoAADQ\nTRAaAABgQG6uNGqU832yCZs3b26+NUHqmicNJCk+XrMlThoAANBNEBoAAGBAR8wzaHYIYkWFVFgo\nTZnSptrjxknffOO8/cFy8fEK37dPhw8fZq4BAADdAKEBAAAGmJ5nsGPHjubnGWRnOwMDT8821Xdz\nc7ZVpKW1Y5PNiY+XR2qqZs6YoS1bthhYAAAAWInQAAAAA5KTzZ00SElJUXh4uPr169f0E9oxz6Ce\nsRaFQYOkwYOVMG4ccw0AAOgGCA0AALBYSYlUXu688dAEk/MM6pmea5Dg7s5cAwAAugFCAwAALJaS\n4nzTbbOZqZ+UlHT50CAzs92hgd3uzB5qa9tVpmnx8Qo/cECHDh3S119/bWABAABgFUIDAAAsZnII\nYlVVlbKzs5ufZ3DihHOC4bhx7Vpn4EBpyBBp9+52lWlafLw8tm9XXFyctm3bZmABAABgFUIDAAAs\nZnIIYk5OjsaPH9/8PIPMTGnaNMndvd1rGWtRCAmRPDwUO3asUo31QAAAACsQGgAAYKHqaiknx3m8\n34TU1FTFxMQ0/wQL5hnUMxYa2GxSfLxi3N0JDQAA6OIIDQAAsFBOjnTFFZK3t5n63xkaWDDPoJ7p\nYYj2r75SVlaWampqDC0CAADai9AAAAALmZxnIHXsSYMpU6SDB6VTpywp11h8vPzS0jRy5Ejt3LnT\nwAIAAMAKhAYAAFjI5DyDkpISVVRUaMyYMc09QTpzRho1ypL1PDyk6dOdOYTlJk6UTp5UzNSptCgA\nANCFERoAAGCh5GRzJw3S0tIUExMjW3N3OWZkSJGRlt71aKxFwc1NmjFD0T4+SktLM7AAAACwAqEB\nAAAWKS52/qJ/9Ggz9TuyNaGe0bkGM2YoprKSkwYAAHRhhAYAAFikfgahhb/ob+Q7Q4Pt253v8i0U\nHe0MDRwOS8s6xcVp0hdfqLi4WCdOnDCwAAAAaC9CAwAALJKZ6ewOMKGmpkZZWVmyN3eXY2mpcwNz\n5li6bnCw8yaIffssLesUGSn33bsVOW2a0tPTDSwAAADai9AAAACLmAwN8vPzNWLECPn5+TX9hI8+\nkmbOlPr1s3xtYy0KffpIU6YoZvhwWhQAAOiiCA0AALCAw2E2NEhNTVV0dHTzT1i/Xlq0yMjaRuca\nxMUpxs2N0AAAgC6K0AAAAAscOiS5uzuP85tw2XkGtbXShx9229Ag+quvlJaWprq6OkOLAACAtiI0\nAADAAvWnDDplCGJqqjOtGDHCyNrTpkmffy5VVhooHhenoKws+fv7a8+ePQYWAAAA7dGq0CAkJERT\np07VtGnTmh/EBACACzLZmnDixAkVFxdr0qRJTT/BYGuCJHl5SZMnS1lZBooHB0u+voqZOJEWBQAA\nuqBWhQY2m01JSUnKyclhyjEAABcwGRqkp6crMjJS7u7uTT9h/Xpp8WIzi59nfK6Bry+hAQAAXVCr\n2xMcRi5qBgCg+6ofgjh9upn6l21NKCyUjh6VoqLMLH6e0dBgxgzFnD5NaAAAQBfU6pMG8+fPV2Rk\npF5++WVTewIAoFs5cMB50+HgwWbqXzY0WL9euvpq5xRGg2JipJQUZ0Biubg4hX/xhfbu3auKigoD\nCwAAgLZqVWiwfft25eTk6MMPP9RLL72krVu3mtoXAADdhsnWhLq6OqWlpTV/3WIHtCZIUkiI85KG\nw4cNFJ8yRb1KShQ2aZIyMzMNLAAAANrKozVPHjJkiCQpMDBQS5cuVXp6uuLj4xseT0xMbPg8ISFB\nCQkJlmwSAICuzGRosGfPHvn7+ysoKOjSBysqpO3bpX/8w8ziF7DZvm1RGD7c4uLu7pLdruh+/ZSa\nmsq/HwAAsEBSUpKSkpLaXafFoUFlZaVqa2vl4+Oj06dP6+OPP9bKlSsbPefC0AAAAFeRmSn9n/9j\npvZlWxM++cT5Tt7X18ziF6kPDW680UDxuDjF5OZqLXMNAACwxMW/yH/00UfbVKfF7QlHjx5VfHy8\nwsPDFR0drUWLFul73/temxYFAKCnqKtzXkXYKUMQO6g1oZ7xYYhffaXU1FSGLgMA0IXYHBb9l9lm\ns/EfeQCAy9mzR1qwwDkM0YTw8HD95S9/kd1ub/xAXZ0UHOxsTxg92sziFykvdw57PHlS6tXL4uKl\npXIMHapgHx+lpKQoJCTE4gUAAHBtbX3P3uorFwEAwLdMzjOoqKjQ3r17FR4e3vTCAwZ0WGAgST4+\n0pgxUm6ugeJ+frKFhipm3DiuXgQAoAshNAAAoB1MhgaZmZkKCwtTr6Z+rd/BrQn16q9eNGLGDMX4\n+BAaAADQhRAaAADQDiZDg++cZ7BokZmFL8PoXIO4OMVUVBAaAADQhRAaAADQRrW1Uk6OFBFhpn6z\nocGRI9LBg1JsrJmFLyM21uBJg7g4RX7xhfLy8nT27FlDiwAAgNYgNAAAoI2++MI5GNDf3/raDoej\n+dBg/Xpp4ULJo8U3J1tm7FiptFQqKTFQPDRU/erqNG7UKOXk5BhYAAAAtBahAQAAbWSyNeHgwYNy\nc3PT8OHDL32wk1oTJMnNzeBcA5vN2aIQHKwUY8cZAABAaxAaAADQRiZDg/T0dNntdtlstsYPHD0q\nbdvmvOexk8TFmW1RiHY4lJGRYWgBAADQGoQGAAC0UWamNH26mdr1ocElnnlGuu02Mz0RLWR0rsGM\nGYoqLlZ6erqhBQAAQGvYHA6Hw5JCNpssKgUAQJdXUyP5+Tl7+319ra8/e/ZsPfzww5o/f/63X/z6\na2ncOCkvTxo2zPpFW6i8XBoyRDpxQmrqNsh2OXtWtf7+8vfwUGFhoQYMGGDxAgAAuKa2vmfnpAEA\nAG2we7c0fLiZwKC2tlbZ2dmKvLj3YdUq6eabOzUwkCQfH2n0aGnHDgPFe/eW+7RpiggNpUUBAIAu\ngNAAAIA2MDnPoKCgQMHBwerfv/+3XzxxQlqzRvqv/zKzaCuZblGw+/gQGgAA0AUQGgAA0AamhyBG\nRUU1/uJzz0nXXy+NHGlm0VaKjZWSk80Vjzp9mrkGAAB0AYQGAAC0gcnQICMjo/EQxFOnpD/+UfrN\nb8ws2AZGTxrExsq+f7/S09OZlwQAQCcjNAAAoJWqq6X8fCk83Ez9S04arF4tLVokhYaaWbANrrhC\nqqyUjhwxUHzwYI3w91fduXM6fPiwgQUAAEBLERoAANBK+fnO9+/e3tbXrqqq0ueff67w+kSirEx6\n4QXpgQesX6wdbDYpJsbcaQNbXJzsw4Yx1wAAgE5GaAAAQCtlZ0sREWZq79ixQ+PGjVOfPn2cX3jp\nJWnBAmnsWDMLtkNcnNkWhSgPD+YaAADQyQgNAABopdxcc60JjeYZVFQ4ByA++KCZxdrJ+FyDb77h\npAEAAJ2M0AAAgFbKzZWmTjVTu9E8gz/9SZozR5owwcxi7RQV5fxZnD1roHhYmKK+/lqZmZmqq6sz\nsAAAAGgJQgMAAFrB4ZDy8qSwMDP1G04aVFZKzzwjPfSQmYUs4O3t7JrIzjZQ3NNTAdOna2Dfvtqz\nZ4+BBQAAQEsQGgAA0ApFRVKfPlJgoPW1T506pSNHjmjChAnSmjXSjBnS5MnWL2Qh0y0KUQMHMtcA\nAIBORGgAAEArmDxlkJWVpfDwcHmcOyf94Q/Sww+bWchCRkODuDjZa2qYawAAQCciNAAAoBVyc82F\nBunp6c7WhP/7f50DA0xNW7SQ8ZMGR45w0gAAgE5EaAAAQCvk5ZkbgpiRkaGo8HDpqae6xSkDSRo9\n2jkI8dAhA8UHDVJEQIDyd+5UdXW1gQUAAMB3ITQAAKAVjJ80KCx0phKRkWYWsZjNJsXFScnJZup7\nz5ih0AEDlJeXZ2YBAABwWYQGAAC0UGWlcxDiuHHW1y4uLlZVVZVGvfxytzllUM90i4K9Xz/mGgAA\n0EkIDQAAaKH8fGn8eMnT0/raGRkZihoyRLYJE6SYGOsXMMj4XIPycuYaAADQSQgNAABoIaPzDFJT\nFXX4cLc7ZSA5Oyny86WqKgPFp06V/dQpZaSmGigOAAC+C6EBAAAtZHSewfvvyz58uDRzppkFDOrX\nz3kCIyvLQHEPD02JjNSBAwdUXl5uYAEAAHA5hAYAALSQqZMGjnPnlFlQoKjEROuLdxCTLQqeM2Zo\namCgsrOzzSwAAACaRWgAAEALOBzO0MDESYN9zz8vH09PBS1dan3xDmJ8roGbG3MNAADoBIQGAAC0\nwKFDkpeXFBhoceHaWmU8+6yi7Hbn/YXdVHS0ZOw9fUyM7F9/rQxCAwAAOhyhAQAALWBsnsFbbynd\n4ZD9mmsMFO84o0dLp09LxcUGig8apKiAAKUnJxsoDgAALofQAACAFjAWGrzxhjJ8fJwnDboxm02y\n26WMDDP1r5g1S6dOntSxY8fMLAAAAJrUqtCgtrZW06ZN0+LFi03tBwCALsnIEMTTp3Vu0yblFhdr\n+vTpFhfveCZDA7e4OEX2768MUwsAAIAmtSo0eP755zVx4kTZunHPJQAAbWHkpMHGjdo1YYJGjBgh\nX19fi4t3vKgog3MNYmNlr64mNAAAoIO1ODQ4fPiwPvjgA/3whz+Uw+EwuScAALqUykqpqEgaN87i\nwuvWKX3UKNm7eWtCvago50mDujoDxadMUdTp00rfvt1AcQAA0JwWhwb333+//vCHP8jNjTEIAADX\nkp8vjR8veXpaWLS2Vlq/Xhk2m6Kioiws3HkGD5Z8faV9+wwU9/CQfdo0ZaSn88sLAAA6kEdLnrR+\n/XoNGjRI06ZNU1JSUrPPS0xMbPg8ISFBCQkJ7dweAACdz8g8g7Q0KShI6QUFuudXv7K4eOepn2sw\ndqz1tYfOmSPPnTtVWFioUaNGWb8AAAA9SFJS0mXfv7eUzdGCuP6BBx7Q66+/Lg8PD1VVVamsrEzL\nli3Ta6+99m0hm43kHwDQI917rzRqlPTLX1pY9L/+S6fr6hT44os6efKkevfubWHxzvPUU9KRI9Lz\nzxsovn69rrvrLn3/xRd10003GVgAAICeq63v2VvUa/C73/1Ohw4d0oEDB/Tmm29q7ty5jQIDAAB6\nMiMnDd57TzlXXKHJkyf3mMBAcp40MDkMMaqsTOmpqYYWAAAAF2vTgAJuTwAAuAqHwxkaWHpzwp49\nUmmpMsrLe8w8g3rTpzt/XtXVBooPHCj7kCHKsOCoJQAAaJlWhwazZ8/WunXrTOwFAIAu59AhyctL\nCgy0sOj770uLFys9I6PH3JxQz8fH2cqRn2+mfuTcucrevVu1tbVmFgAAAI1wFQIAAJeRm2vxKQNJ\neu89ackSZWRk9LiTBpLz6kVTLQr+V16pIR4eKigoMLMAAABohNAAAIDLsLw14fhxKTdX34SH69ix\nYxo3bpyFxbsGo3MN4uMVde4ccw0AAOgghAYAAFxGbq7FQxA/+ECaN0+Z+fmaPn263N3dLSzeNdRf\nu2jE0KGy+/go45NPDC0AAAAuRGgAAMBlWH7S4HxrQnp6eo+bZ1BvyhTpyy+l8nIz9aOio5WekmKm\nOAAAaITQAACAZlRWSkVFkmUdBFVV0saN0jXX9Nh5BpLUq5czOMjONlN/2pIlKiguVlVVlZkFAABA\nA0IDAACakZ/vDAw8PS0quGmTNHWqHAEBPfqkgWR2rkGfefM0zs1NO3JyzCwAAAAaEBoAANCMvDyL\n5xmsWyctWaLDhw/L4XBo+PDhFhbvWozONRg9WnZPT2Vs2GBoAQAAUI/QAACAZlgaGtTVOUODa69V\nRkaG7Ha7bDabRcW7HqM3KNhsipo4URkbNxpaAAAA1CM0AACgGTt3WhgaZGdLPj7S2LFKT0/vsfMM\n6o0ZI506JR09aqa+/XvfU/quXWaKAwCABoQGAAA0weFwhgaTJ1tU8HxrgqSGkwY9mZubFBVlrkVh\n4g036HBZmUpLS80sAAAAJBEaAADQpGPHnMHB4MEWFfzgA2nRItXV1SkzM1ORkZEWFe66TM418AgL\n0zR3d2V+9JGZBQAAgCRCAwAAmpSf77w20JKxA6dPSwUFUnS09uzZo4CAAAUEBFhQuGszOtfAzU1R\nI0Yo4913DS0AAAAkQgMAAJpkaWtCRoZzOELv3i4xz6BeVJQzNHA4zNS3x8YqPS3NTHEAACCJ0AAA\ngCbl51sYGqSkSLGxklxjnkG94GCpTx/pwAEz9aOuv14Zhw6ZKQ4AACQRGgAA0CRToYErnTSQzLYo\nhC5apMpz51Sye7eZBQAAAKEBAAAXq6uTdu2SJk2yoJjD0RAaVFdXKz8/XxERERYU7h5Mhga2Xr0U\nNXCgMv7+dzMLAAAAQgMAAC5WVCT5+Un+/hYU279f8vKShg1TXl6eRo8erX79+llQuHuon2tgin3y\nZKV/+qm5BQAAcHGEBgAAXIR5BtaJjJR27JDOnTNT375ggTJoTwAAwBhCAwAALsI8A+v4+UnDhzvb\nPUyIuvVWZZSVyVFebmYBAABcHKEBAAAX4aSBtaKjJVM3IwaFhMind2/tf/ttMwsAAODiCA0AALjI\nzp3SlCkWFKqokPbskaZNU3l5uQ4cOKDJlqUR3YfJ0ECSokJClP7uu+YWAADAhREaAABwgXPnnO/z\nJ0ywoFhGhhQWJvXurezsbIWFhcnT09OCwt2L6dDAPmOGMjIyzC0AAIALIzQAAOAC+/ZJw4ZJffta\nUMzF5xnUmzJFOnhQKiszUz9q6VKll5RIZ8+aWQAAABdGaAAAwAWYZ2A9T08pPNx58MKE6bNmKVfS\nuc2bzSwAAIALIzQAAOACloUGDoeUmtoQGqSmprpsaCBJdru5FgVfX1+NHDhQO//2NzMLAADgwggN\nAAC4gGWhwb59Up8+0tChOnz4sM6cOaMxY8ZYULh7io6W0tPN1Y+LjVXKJ5+YWwAAABdFaAAAwAUs\nCw0uaE1ISUlRbGysbDabBYW7p/phiA6Hmfqxixcr+cQJ6dAhMwsAAOCiCA0AADjvzBmpqEgaO9aC\nYheFBnFxcRYU7b5GjpTq6sy9p4+bOVMpnp7Shg1mFgAAwEURGgAAcF5BgXTFFc7Bfe12QWiQnJys\n2POfuyqbzezVi2PHjtUpSV+9846ZBQAAcFGEBgAAnGdZa0J5ubR3rzRtmqqqqrRz506XvW7xQiZD\nAzc3N8XGxipl82aputrMIgAAuCBCAwAAzrMsNMjIcN4x2KuXsrKyNGHCBPXt29eCwt2bydBAkmJn\nz1ayj4+UnGxuEQAAXAyhAQAA55kagujq8wzqRUVJOTlSTY2Z+nFxcc65Bh9+aGYBAABcUItDg6qq\nKkVHRys8PFwTJ07Ub37zG5P7AgCgw5kIDZhn8C0/P2n4cOfP2QS73a6c48dV/cEHZhYAAMAFtTg0\n8PLy0qZNm7Rjxw7l5eVp06ZN2rZtm8m9AQDQYUpLpRMnpJCQdhZyOKTUVCk2Vg6Hg5MGFzHZouDt\n7a2x48crp6hIOnLEzCIAALiYVrUn1PdjVldXq7a2VgMGDDCyKQAAOtquXdLEiZJbexv39u6V+vWT\ngoNVWFgom82mESNGWLLHnsD4XIPYWCWHhnL1IgAAFmnVP43q6uoUHh6uoKAgzZkzRxMnTjS1LwAA\nOpTJeQY2m82Cwj2D6dCAuQYAAFirVaGBm5ubduzYocOHD2vLli1KSkoytC0AADrWzp3SlCkWFGKe\nwWVNmSIdPCiVlZmpHxsbq+TDh6VPP5XOnTOzCAAALsSjLd/k5+ena665RpmZmUpISGj4emJiYsPn\nCQkJjR4DAKAry8+XliyxoFBKinTXXec/TdHy5cstKNpzeHo6b6PMzJTmzrW+fmhoqM7V1urQ0KEa\nnpIizZpl/SIAAHQDSUlJlvyi3+ZwOBwteeLx48fl4eGh/v3768yZM1qwYIFWrlypefPmOQvZbGph\nKQAAuhSHQwoMdJ42GDKkHYXKy50FTpzQ6XPnNGjQIJ04cUK9e/e2bK89wa9+JQUESKYuYrruuut0\nq82mmydMkH73OzOLAADQzbT1PXuL2xNKSko0d+5chYeHKzo6WosXL24IDAAA6M6OHXMGB4MHt7NQ\nZqYUFib16qWMjAyFhYURGDTBbmeuAQAA3UWL2xOmTJmi7Oxsk3sBAKBT1A9BbPe8wvR05ztiMc/g\ncqKjpfvucwY1JmZExsbG6lf//KdzeEJJSTuPjwAA4Nrae7EUAADdnmU3J6SlOd8R69ubE3CpkSOl\nujrp0CEz9SMjI7Vr926dmTOHqxcBAGgnQgMAgMuzLDQ4f9LA4XAoJSWFkwbNsNnMXr3Yp08fTZo0\nSZlXXEGLAgAA7URoAABweZaEBkeOSGfPSqNGae/evfL29lZwcLAl++uJoqOdGYspcXFxSvHwkDZu\nlGpqzC0EAEAPR2gAAHBpdXUWhQZpac55BjYb8wxawORJA8k51yA5P18aMcLsQgAA9HCEBgAAl7Z3\nr/O6RX//dhZKT2eeQStERUnZ2eYOAcTFxSklJUWOq66iRQEAgHYgNAAAuLTsbGnaNAsK1Z80EDcn\ntISfnzR8uPOUhwnDhw+Xp6envgwLIzQAAKAdCA0AAC4tJ0eKiGhnkdpaKStLsttVWlqqwsJChYWF\nWbK/niw+Xtq82Vz9uLg4pZw9K+3fL331lbmFAADowQgNAAAuLTvbgtCgoEAKCpIGDFBaWpoiIiLk\n6elpyf56snnzpE8/NVc/NjZWyenpzoU++sjcQgAA9GCEBgAAl+VwOE8atLs9gXkGbTJnjrRli/m5\nBlq4UNqwwcwiAAD0cIQGAACXVVQk9e4tDR7czkLMM2iTQYOclxtkZZmpP23aNO3Zs0flM2dKH3/s\nbCMBAACtQmgAAHBZlg1BPH/SoK6uTmlpaYqJibGgqGsw2aLQq1cvhYeHK/3IEWnoUOfrBAAAWoXQ\nAADgsiwZglhZKe3ZI4WFKS8vT0FBQRo0aJAl+3MFc+dKn31mrv7MmTO1detWZ4sCtygAANBqhAYA\nAJdlyUmD7Gxp0iTJy0tJSUmaM2eOJXtzFbNnO7s7qqrM1J8zZ46SkpIIDQAAaCNCAwCAy7LkpMEF\n8wySkpKUkJDQ7n25El9fafJkKTnZTP0ZM2YoMzNTVRERzhMhx46ZWQgAgB6K0AAA4JKOHpXOnJFG\njmxnofPzDGpra7VlyxZCgzYw2aLg4+OjKVOmKCUry7nQxx+bWQgAgB6K0AAA4JLqr1q02dpZ6PxJ\ng/p5BoPbfRWD6zE5DFGSEhISaFEAAKCNCA0AAC4pO9uC1oSjR6XSUumKK7Rp0ybmGbRRXJyUny+V\nlZmpP2fOHG3atMkZGnD1IgAArUJoAABwSfUnDdolPV2KipLc3Jhn0A5eXs6xEFu2mKkfFxen7Oxs\nVQ4cKAUFSZmZZhYCAKAHIjQAALgkS04aXDDPYOvWrYQG7WCyRcHb21thYWFKSUmhRQEAgFYiNAAA\nuJxTp5ydBVdc0c5C5+cZ7NixQ8HBwRo0aJAl+3NFHTbX4KqrpA0bzC0EAEAPQ2gAAHA5O3ZIYWGS\nu3s7itTVSRkZkt2uTZs2ccqgnaZPl4qKzN2ImJCQ4JxrMHOmVFAgHT9uZiEAAHoYQgMAgMvJzrZg\nnsHevZKfnxQUpKSkJIYgtpOHhzRrlrRpk5n6cXFx2rFjh07X1EgJCVy9CABACxEaAABcTk6OdfMM\nampqtG3bNs2ePduSvbkyky0K/fr1U3h4uJKTk5lrAABAKxAaAABcjiUnDc7PM8jJydHw4cMVGBho\nyd5c2bx50mefmas/Z84c51yDhQuljz5ytpgAAIDLIjQAALiUykrpwAFp0qR2FkpLk6KjmWdgoUmT\npPJy6eBBM/Ub5hqMHCkFBEhZWWYWAgCgByE0AAC4lLw8acIEqVevdhSpqpJ27ZIiIphnYCGbTZo7\n11yLQmxsrPLy8lRRUcEtCgAAtBChAQDApeTkWNCasGOHNG6cznl6avv27Zo1a5Yle4PZuQZ9+/ZV\nRESEtm/f7gwNPvrIzEIAAPQghAYAAJeSnW3BEMRt26TYWGVnZyskJEQBAQGW7A3OkwaffSY5HGbq\nN8w1mDnTGf5UVJhZCACAHoLQAADgUiwZgvjpp9K8ecwzMCA0VPLykgoKzNRvmGvQt68UGSlt2WJm\nIQAAeghCAwCAy6iudr4ZnTq1nUW2b5cSEphnYIjJFoXY2Fjl5+ervLzc7EIAAPQQhAYAAJexe7cU\nEiL169eOIunp0pgxOufrq+TkZOYZGHDNNdLbb5up7eXlpcjISG3bto3QAACAFiA0AAC4jJwcC+YZ\nnG9NyMzM1OjRozVgwABL9oZvXX2185aLoiIz9RvmGkRFOe/f/PprMwsBANADEBoAAFwG8wy6h969\npWXLpLVrzdRvmGvg6SnFx0ubNplZCACAHqDFocGhQ4c0Z84cTZo0SZMnT9YLL7xgcl8AAFiu3ScN\nTp92Jg/x8cwzMOy226TXXzdzi0J0dLR2796tsrIyWhQAAPgOLQ4NPD09tWrVKu3atUupqal66aWX\nVGBqtDEAABarrZVyc9t50mDrVikiQtWenkpJSVF8fLxl+0NjM2dK5eXONgWreXl5yW63a+vWrYQG\nAAB8hxaHBoMHD1Z4eLgkydvbWxMmTFBxcbGxjQEAYKU9e6TAQKl//3YUOd+akJGRobFjx8rf39+y\n/aExNzdp+XLpjTfM1J8zZ46zRWHyZKmsTDp40MxCAAB0c22aaVBYWKicnBxFR0dbvR8AAIxISpLa\nfdHBZ59J8+Zpw4YNmj9/vhXbwmXcdpv09787T4lYbf78+froo4+c6QSnDQAAaFarQ4OKigrdcMMN\nev755+Xt7W1iTwAAWO78IYG2O3FC2rtXstv1r3/9S4sWLbJsb2jaxInSoEHS5s3W17bb7frqq690\n8A1Wy9wAACAASURBVOBBQgMAAC7DozVPPnfunJYtW6bbbrtN11133SWPJyYmNnyekJDAVGkAQJdQ\nV+cckP/cc+0o8v/bu/O4Ksv8/+Ovw6LpuOaCC66EK6KmSG4pbqTilqZmqKk5mZZjLtmU85upKZPB\nqXSsrPlabmiZmQsouZLivpNaYC6FlOaGGyrb/fvjTq3R6njOfQ4cfD978FA6N5/rAxdHuD/nuj7X\nhg3QsiVpp0/z3Xff0bx5c8vyk98WGQnz50O7dtbG9fb2pnPnzsTFxTGyc2f429/Mros2m7UDiYiI\n5JGEhATziGEn2QzDvr7EhmEwePBgypQpw1tvvXV7IJsNO0OJiIi41b590K8fJCc7EWTkSKhZk/+W\nLMmGDRtYsGCBZfnJb/vhB6hf3/yzSBFrYy9atIg5c+YQFxcHNWvCihXmYCIiIgWQo/fsdm9P2Lx5\nM/Pnz2fDhg00btyYxo0bEx8ff9cDioiIuJvTWxN+ESQ2NpauXbtakpf8sUqVICTEvJ+3WqdOndi0\naRMZGRnmN8jatdYPIiIi4uHsXmnwh4G00kBERPKprl1hyBDo08fBACdOQKNGXPv+e/wqVuTo0aOU\nKVPG0hzlt82dC4sXw/Ll1scOCwtj3LhxRFy+bHZddMUgIiIi+YDLVxqIiIh4oqwsSEyEsDAngqxf\nD2FhJGzcSHBwsAoGbtarF2zcCGfOWB87IiKC2NhYs2nCxo2QnW39ICIiIh5MRQMRESnQdu40t6s7\ndZ//89aEuLg4bU3IA8WLQ5cusGiR9bG7du1KXFwcRrlyULUq7Npl/SAiIiIeTEUDEREp0JzuZ2AY\nsG4dRrt2xMbG6qjFPHLjFAWr1a5dm8KFC5OUlKSjF0VERO5ARQMRESnQ1q938ri+lBSw2TiUmUlu\nbi711V0/T3TsCEeOmG9WstlsdO3a1dyioKKBiIjIbVQ0EBGRAisjw9ye0Lq1E0FubE1YuZKIiAhs\nNptl+Yn9fH3NYzNjYqyPHRERYR67+PDDsGMHXL1q/SAiIiIeSkUDEREpsLZsgYYNzT3xDtNRi/nG\nwIEwezbk5Fgb9+GHH+bgwYOcvn4dgoNh82ZrBxAREfFgKhqIiEiBtW6dk1sTcnMhIYFzDz7Ivn37\nCHPqCAZxVkgIVKgAS5daG7dw4cK0b9+e+Ph4bVEQERH5HyoaiIhIgbV+vZNNEPftg/Ll+SIpiTZt\n2lCkSBHLchPHjB8P0dFmf0orqa+BiIjInaloICIiBVJ6Ohw6BA895ESQXxy1qFMT8ocePeDsWet3\nEHTp0oXVq1eT1bQpJCfD6dPWDiAiIuKhVDQQEZECaeNGs2Bw331OBFm6lOxOnVi1ahVdunSxLDdx\nnLc3jB0LU6daG7dixYoEBASwedcuc7VBfLy1A4iIiHgoFQ1ERKRAcrqfQUoKHDnCthIlqFKlClWq\nVLEsN3HO4MGwdau5IMBKERER5haFiAiIjbU2uIiIiIdS0UBERAokp/sZzJ0LAwYQ98UXOjUhnyla\nFJ55Bt5809q4Xbt2NY9e7NIFVq+GrCxrBxAREfFAKhqIiEiBc+oUpKbCgw86GCA31ywaPPkksbGx\n6meQD40cCYsWwU8/WRezSZMmnD9/niNXrsADD0BionXBRUREPJSKBiIiUuBs2ABt2oCPjxMBypTh\nu5IlOXnyJM2aNbM0P3Fe+fLQrx/MmGFdTC8vr1urDbRFQUREBFDRQERECqD1653sZzBnDgweTFxc\nHJ07d8bb29uy3MQ6Y8fCzJmQkWFdzJtHL0ZEQFycdYFFREQ8lIoGIiJS4DjVBPHSJVi+HAYMYMmS\nJXTv3t3S3MQ6tWpBy5Ywe7Z1MTt27Mj27ds5V60aXLwIhw9bF1xERMQDqWggIiIFyvHjcPkyBAU5\nGGDxYmjThrSsLHbv3q0miPnc+PFmQ8ScHGviFS9enEceeYRFixdD165abSAiIvc8FQ1ERKRAWbsW\nwsLAZnMwwJw58OSTLFy4kEcffZQiRYpYmp9Yq2VLs7/B0qXWxYyMjGT+/Plm0UB9DURE5B6nooGI\niBQoc+ZA374OfvCxY3DwIHTtyvz584mMjLQ0N3GN8eMhOhoMw5p44eHhJCcnczQgALZvN7cpiIiI\n3KNUNBARkQLjwAE4ehS6dXMwwNy50L8/XyUnc/bsWdq0aWNpfuIaPXrA2bOwebM18QoVKkS/fv1Y\nsGyZuZRhzRprAouIiHggFQ1ERKTAeP99eOop8PV14INzc29uTYiJieGJJ57Ay0s/Jj2Bt7d5ksLU\nqdbFvLFFwdAWBRERucfptyERESkQrlyBBQvMooFDEhOhaFFyGzUiJiZGWxM8zODBsGULJCdbEy80\nNJTs7Gx2V60KK1eaRSUREZF7kIoGIiJSIHz8MbRqBVWqOBhgzhwYPJgvN26kTJkyBDl8/ILkhaJF\n4ZlnzJMUrGCz2czVBhs2QNmysHOnNYFFREQ8jM0wrGkbZLPZsCiUiIjIXWvaFP75T+jc2YEPvnIF\n/P3h0CGGTZpE3bp1GT9+vOU5imv99BPUrm2uNihf3vl4hw8fpnXr1pwYOBCf++4zv8FEREQ8lKP3\n7FppICIiHm/XLrMRXni4gwE+/xyaN+dqqVJ8/vnnPP7445bmJ+5Rvjz06wczZlgTLzAwkOrVq7PW\nzw/i4qwJKiIi4mFUNBAREY83cyY8/TQ43Lfw5waIsbGxNGnShMqVK1uan7jP2LHm90NGhjXxBg4c\nyLzdu+G77yAtzZqgIiIiHkRFAxER8Wjp6fDZZzBkiIMBDhyAr76C7t2ZP3++GiB6uFq1zFMSZ8+2\nJl7fvn2JW7WKS+3ba7WBiIjck1Q0EBERjzZ/vrktwc/PwQDR0fCXv3Dm8mUSEhLo1auXpfmJ+40f\nbzZEzMlxPla5cuVo3bo1S8uW1dGLIiJyT1LRQEREPJZhmEvRR4xwMMD335s3gs88w6JFi+jSpQsl\nSpSwNEdxv5Ytzf4GS5daEy8yMpL533wDCQlw9ao1QUVERDyEigYiIuKxNm+G7Gxo08bBAG++CUOH\nQqlS2ppQwIwfby4iseJgp27durFj715+DA62rhIhIiLiIXTkooiIeKzISPOoxTFjHPjgs2chMBAO\nHODI1as0b96ctLQ0fH19Lc9T3C8nB+rUgY8+glatnI83ZMgQgr28eD4pCXbsAJvN+aAiIiJupCMX\nRUTknnL6tNmXbvBgBwPMmAGPPgqVKhETE0P//v1VMChAvL3NkxSmTrUm3sCBA5m/b5/ZeTMx0Zqg\nIiIiHkArDURExCNFR8OhQ+YryXftyhWoUQM2bSI3MJBatWoRExNDaGio5XlK3snIgOrVYdMmqF3b\nuVg5OTlUq1aNVYMH0+DQIfj8c0tyFBERcReXrzQYOnQofn5+NGjQ4K4HERERsdKpU2Y7gueeczDA\nhx9C69ZQuzYrVqygVKlSNGvWzNIcJe8VLQojR8K//+18LG9vb0aMGMFbqanmSoPDh50PKiIi4gHs\nXmmwadMmihUrxqBBg/jqq69uD6SVBiIi4gaGAd26QXAwTJ7sQICsLHjgAfj0U2jWjFatWvHcc8/R\nr18/y3OVvHfmDNSqBQcPQsWKzsU6e/YsgYGBHIiMpFJurrnFRURExEO4fKVB69atKV269F0PICIi\nYqUPPoAff4R//MPBAJ98AjVrQrNmbN26lbS0NHr37m1lipKPlC0LTzwB06Y5H6tMmTJERkYyPTcX\nYmLg3Dnng4qIiORzaoQoIiIeIyUFXn4Z5s+HQoUcCGAYEBUFL74IQHR0NGPHjsXHx8faRCVfGTsW\n/vtfuHjR+VjPP/88//fxx1zq0gXef9/5gCIiIvmcigYiIuIRsrJg4ED4+9+hbl0Hg6xcCT4+0KkT\nKSkpbNq0iaFDh1qap+Q/NWpAeLg19/g1atSgQ4cO/LdiRXN7Qmam80FFRETyMUtfWvnHL9aKtm3b\nlrZt21oZXkRE7mGvvw6lSsGoUU4EiYqCF14Am40333yTESNG8Kc//cmyHCX/euEF6NoVRo+GwoWd\nizVhwgR69erFc3Xq4PvxxzBokDVJioiIWCghIYGEhASn49zVkYvHjx+nW7duaoQoIiJutX07dO8O\ne/dCpUoOBklMNG/uUlL46dw5ateuzTfffIOfn5+luUr+FR4O/fqBFYtL2rVrx9AmTYhcs8b8xrTZ\nnA8qIiLiQi5vhPj444/TokULUlJSqFKlCh85dDC2iIjI3bl8GSIj4Z13nCgYALzyitkQwceHd955\nh759+6pgcI+ZOBH+9S/IzXU+1oQJE4hevRojMxPWr3c+oIiISD51VysNfjeQVhqIiIgLPP00XLsG\nc+Y4EWTzZrPykJJCRlYW1atXJzExkVq1almWp+R/hgHNmsGkSdCjh7OxDIKDg/l3hw50SkmBuDhr\nkhQREXERl680EBERcafcXBg3ztxVMH26k8FurDLw9eWjjz6iZcuWKhjcg2w2s7dBVJRZQHAulo3x\n48cTnZQEu3fD119bk6SIiEg+o5UGIiKS71y7ZrYfOHUKli6F0qWdCHZjlUFyMjne3tSqVYt58+bR\nokULy/IVz5GTA3XqwEcfQatWzsXKzMykZs2arOjVi8YpKbBqFXjp9RgREcmftNJAREQKhHPnoGNH\n81XhL75wsmAAt1YZFCrEkiVLqFChggoG9zBvbxg/3lxt4KxChQoxevRopp45A1evwuTJzgcVERHJ\nZ7TSQERE8o3jx6FzZ4iIMG/qnH7RdssWGDAAUlIwfH0JDQ3lpZdeomfPnlakKx7q6lWoWRPWrIGg\nIOdiXbhwgZo1a7Jn5Uqq9ewJCxZAWJg1iYqIiFhIKw1ERMSj7dkDLVvCyJEQHW3RKu9frDJYtWoV\nly5donv37hYEFk9WpAg895w1CwNKlizJU089xWv/938wd665FebkSecDi4iI5BNaaSAiInkuNhaG\nDoX334devSwK+otVBtdycwkKCmLGjBk88sgjFg0gnuzyZahbFz7+2CxWOSM9PZ169eqxdOlSmsXG\nmt0716wx90KIiIjkE1ppICIiHscwYMoU81jF5cstLBiAucrgpZegUCGio6Np2LChCgZyU7Fi8K9/\nwbPPms0RnVGqVCmioqIYNWoUOZMmmf/z1VedT1JERCQf0EoDERHJE1evwrBhcPiweUJC5coWBt+6\nFfr3h8OHOZaWRkhICHv27KFq1aoWDiKezjCgbVt4/HEYMcLZWAYPP/wwkZGRPN2jBzRpArNnm109\nRURE8gFH79lVNBAREbdLS4OePSEwEGbNMveYW+qRR+DRR+HPf6Z79+40b96cv/71rxYPIgVBUhJ0\n6ABffw1lyjgbK4kOHTpw8OBByh04YG6P2b0bKlWyJlkREREnaHuCiIh4hO3bITQUeveGmBgXFAy2\nbjXvAJ98khUrVpCcnMzYsWMtHkQKiuBg6NcPbuwqcC5WME888YRZoAoLg1GjzOLVxYvOBxcREckj\nWmkgIiJu88knZtf6WbOgWzcXDJCbCy1awNNPc7V/f+rXr8/MmTPp1KmTCwaTguL8ebMp4qpV0Lix\nc7EuXLhA3bp1WbJkCQ81a2Y2Tdi9G+LjoXRpaxIWERFxgFYaiIhIvvbxxzBmDKxb56KCAZjVCC8v\nGDyYqKgomjRpooKB/KHSpeG118z7e2df/yhZsiTR0dGMHDmSHMOAd94xj2do3x7OnLEmYRERETfS\nSgMREXG5xYvNG7I1a6BBAxcNcuYM1K8PX3zBkeLFCQ0NZe/evVSpUsVFA0pBkptrbpsZPRoGDnQu\nlmEYtG3bln79+jFy5EizEvHyy7BiBaxdC35+1iQtIiJyF9QIUURE8qVly8wjFePjoVEjFw40fDgU\nLYrx9ttERETw8MMPM3HiRBcOKAXN9u1mC4Kvv4YSJZyLdeDAAcLCwjh48CDly5c3Cwf//CcsWGAu\nt7H0uBAREZE/pqKBiIjkO7Gx5rGKK1eaJ9C5zLZtN+/2Fq9Zw6RJk0hKSqJQoUIuHFQKoqFD4f77\nYepU52ONGzeOU6dOMW/ePGw2m/k/o6Lgv/81CwfVqjk/iIiIiJ1UNBARkXwlPh4GDTILB82auXCg\nnBwICYFx4zjYqBFhYWGsWLGC0NBQFw4qBdWpU2YzxBkzzDqUMy5fvkzz5s15+umnefbZZ289MG0a\nvPWWWTgICHBuEBERETs5es/u44JcRETkHrdmjVkwWLbMxQUDgPfeg5IlORseTo+HHmLq1KkqGIjD\n/PwgLg7Cw6FkSbN/oaOKFSvG8uXLad68OXXr1qX9jWB/+QsULgxt25o9DmrXtiR3ERERV9BKAxER\nsdTmzdCzJ3z+ObRq5eLBTp2CoCCy1q4l/PnnadKkCdHR0S4eVO4FGzdCnz7WrJRJSEigX79+JCYm\nEhgYeOuBjz6CSZNg9WqziaeIiIgLaXuCiIjkuaQk6NgR5s41X6l1uUGDoEIFns3I4NixYyxfvhxv\nb283DCz3ghUrzP6a69dDvXrOxZo5cybTp09n69atlCxZ8tYDCxbAuHGwapWLO4WKiMi9ztF7di8X\n5CIiIvegI0egc2eYPt1NBYONGyEhgff9/Vm3bh0LFixQwUAs1a2b2RDxkUfgu++cizVixAjCwsIY\nMGAAOTk5tx4YMMBsoBAeDjt3OjeIiIiIC2ilgYiIOO2HH6B1a3jhBfN4RZc7dw5CQ/ly4ED6vvPO\n7cu+RSw0fbp5X79pk9nzwFFZWVmEh4cTEhJCVFTUrx9cscI8auTzz6FlS+cSFhERuQOtNBARkTxx\n7pz5IumwYW4qGGRlQZ8+HGvThn7vvktMTIwKBuJSo0fD44+bKw7OnHE8jq+vL59++imLFy9m3rx5\nv36wWzeYN89sCJKQ4FS+IiIiVlLRQEREHHblCkREQKdO8Ne/umFAw4BRo/gWCN+4kZdeeokOHTq4\nYWC51/3jH9C1q3kc45dfOh6nTJkyLF++nHHjxrFo0aJfPxgeDosWwWOPmc0RRURE8gFtTxAREYdc\nvw49ekDFivDhh2CzuWHQt95i84wZ9L58mX+88gojRoxww6Ait8THw5Ah5qqaSZPAx8HDq/ft20f3\n7t0ZOXIkEydOxPbLJ9DmzdCrl/nEioiwJnEREbnn6fQEERFxmytXzHuakiVh4ULHb5zuSmwsCwcO\n5C/e3syLiSHcLd0WRW73448QGQnZ2RATA/7+jsVJS0sjIiKCJk2a8N577+Hr63vrwR07zC0L774L\nvXtbk7iIiNzT1NNARETcIj3dXEVdubL7CgbG/v283q8fL953H+s2bFDBQPJUxYrm7oHwcGja1Oxh\n6IjKlSuzadMmTp48SefOnUlPT7/1YLNm5rKGUaPMYxlFRETyiIoGIiJit9OnoV07ePBBmDXLPQWD\nzNRUhrZsyed+fmzbs4cGDRq4flCRP+DtDS+9BJ99Bs8+C888Axcv3n2cYsWKsWzZMurVq0fLli05\nfvz4rQcbN4a1a2HCBPjoI8tyFxERuRsqGoiIiF3S0qBNG+jcGaZNAy83/AQ5tn8/4UFBnK9cmS+/\n+oqKFSu6flCRu9CyJezfb25VCAqC2Ni7j+Ht7c306dN5+umnadGiBfHx8bceDAqC9evh//0/8+xH\nbQUVERE3U08DERH5Q0ePQseO8Oc/w8SJrh8vIyODKRMm8M7MmYxr1oyJGzfi/cv93iL50Pr15nMk\nJMQsrJUvf/cxVq9ezciRI6lfvz5vvvkmAQEB5gNHj0L37tC8OcyYAYULW5u8iIgUeOppICIiLrF/\nv7nCYNw41xcMDMNg0aJF1KlenW9nzWL/G2/w0tatKhiIR2jXDpKSoGpVaNAA5s69+4UBnTp14uDB\ngzRv3pzQ0FBefvllLl++DDVrwtatcPYshIWZ3RhFRETcQEUDERG5o5wciIqCDh1g6lQYOdK14yUl\nJREWFsbk559nfk4OC9aswf+FF1w7qIjFihY1nzerVsFbb5nPn1277i5G4cKFefHFF9m/fz/fffcd\ndevWZeHChRjFisHixdCli7mcYft213wSIiIiv6DtCSIicpsjR2DwYPD1hdmzoVo114yTk5PDF198\nwaxZs0hMTOSVmjUZnpGB9/LlrhtUxE2yssyGoa+9Zh6G8OqrZouCu5WYmMjo0aMxDIPhw4czYMAA\nSm3aBMOGmRWKIUOsT15ERAocbU8QERGnGQZ88AE89BD06QPr1rnm3v3o0aNMmjSJatWq8corrxDe\nqBGHq1VjROXKeG/erIKBFAi+vjBiBBw+DK1aQfv2EBkJ3357d3FatWrFzp07mTJlCgkJCVSvXp2B\nixaREB2N8cYb5iCnT7vmkxARkXue3UWD+Ph46tSpQ2BgIFFRUa7MSURE8sAPP0BEBLz/Pnz5JYwZ\nY90JCYZhkJyczAcffEC7du146KGHyMjIIH7aNLbXqMGfp0+nxGOPwaJFUKyYNYOK5BNFisDYsWax\noHZtsyj31FPmaYrXrtkXw9vbm/DwcBYtWsS3335LkyZNeG7qVAKzs3ll9242BwSQOW4cnDrl2k9G\nRETuOXb9OpiTk8Ozzz5LfHw8hw4dYuHChXz99deuzk08QEJCQl6nIHlA815wZGWZR8T16QN168KD\nD5q91urV+/V1dzvnhmFw8OBB3n33Xfr160elSpXo1KkTiYmJjBw5khOLF/Pm4cMEPfccNG1qdoaf\nMME95ziK3fRct1bx4vC3v0FKCtSoYZ6iWL48dO1qnqaYkmJf48SyZcsyZswYkpKSWPjJJ1xo3Zrn\nqlWjzH/+Q3t/f15t0YIvlyzhmr0ViV/QnN+bNO/3Hs253A0fey7asWMHDzzwANWrVwegf//+LFu2\njLp167oyN/EACQkJtG3bNq/TEDfTvHu+r76COXNg/nwICIAnnzT3Xpcseefrf2vOr1+/zrfffkty\ncvLNt5SUFL755htKly5NmzZt6NKlC1GvvUb169dh3z545x04ftw8iuHTT+G++1z5qYoT9Fx3jfvv\nh5dfNt/OnzdXHMTHw7/+ZW5pCAqCChVuvVWsaP5ZrhyUKmU+TwsVMvemhoSEEBISAkB6ejqJy5bx\n5fTpTOjbl4M2GzX9/andoIH5Vrv2zbdSpUrdMTfN+b1J837v0ZzL3bCraJCWlkaVKlVuvu/v7892\ndewVEcmfsrPhyhVyL1zidPI5Th46x7nDZ7lw7BxX085yLe0cPlcvMaD6JSbVv0Rh4xKZH1wi883L\nHM/N5ZKPDxdtNi55eXHRy4tLwOaTJxm/dCkns7M5df06p65e5eSlS1zIyKBa+fLU9venduXKtPT3\nZ2jPntSuUAG/06fN8xrfegueecbsUxAcDEOHQv/+5t2RyD2udGl47DHzzTDg66/NbQwnT5qnKh44\nYBYVTp6En36CCxcgPR0KFzYLCKVKmTHKloWyZUtRrtxgKj4xmAlPnqT2xrdI/2otqavXkLwxkZUl\nivN2bg4p587h7eODX/nyVKhUCT8/P/z8/KhQoQK7du1iwYIFFC9enBIlStz8s1ixYhQuXJhChQpR\nqFAhfHx8sNlsef3lExERN7CraGDvD4VWhf3uanCDu+/caMn5DHYGcWQsRz4nx/JwbBzHPqfflpaT\nQdzkGXfxMS74+jgQ0p3fe+6aKavn9veczLnKZ6+/7eBH28+hz8nBU1zc+fVzxI3vWRsGYGAzDGw/\nP2IzcrGRi5eRixe52DDIxYssvLiGN1leXmTZIMsG2TaDHAyyjBwyD+Vg+9pGIR8fCvn64uvry5+K\nFKFE0aIUL1KEEvfdR/HChSlRuDCXr1yhfI0aBPv44Af45ebil5lJucxMfG5sKbh61bzb+fZbc5tB\nQAC0aQOjR5v7HYoWdeNXTMTz2GzmU+V/twf9L8OAK1fM4kF6Opw7B2fOmL0Qz5yB77+H3acrcOZi\nFOkl4WLVHMqdS6bWyZ0Mz95BM9sOil/7nsuX0jh/5Bjf2+7jhFdhfvDy5WjOZeav3s4VI4fLRg5X\nyOWKkc0VI4csDLKNXLIw/x3xwYYvXvjYvMy/27x+ft+GD15422yY/5n7YW0//6tlA7xsN/5+63Fu\n/v3GlebvoJ6waelWxvnTH2X3ffZl1r7xnltyuZP8/dUrmL7Lvsy6PJxzca2Y3TupGlTVsnh2FQ0q\nV65MamrqzfdTU1Px9/f/1TUBAQFsPnLEssTEc/yYezWvU5A8cDr7el6nIL8r9+e3bPOP32AYBtcy\nM7mWmQnA2fT037x22922fJcC4ZVXXsnrFMQC3wCbgFlwe8XTyICcDMgx3/06N8OumNkYZJMDxs8f\nqJO3PVpqpn3zLgWH5rzgqtbgzqdQBQQEOBTPZtjxslx2dja1a9dm3bp1VKpUiWbNmrFw4UL1NBAR\nEREREREpwOxaaeDj48OMGTMIDw8nJyeHYcOGqWAgIiIiIiIiUsDZtdJARERERERERO49DveWOXfu\nHB07dqRWrVp06tSJ9N/YB/vGG29Qv359GjRowIABA7h+XfugPZm9856enk6fPn2oW7cu9erVY9u2\nbW7OVKxi75wD5OTk0LhxY7p16+bGDMUV7Jn31NRUwsLCqF+/PkFBQUyfPj0PMhVnxcfHU6dOHQID\nA4mKirrjNaNHjyYwMJCGDRuyd+9eN2corvBH8x4TE0PDhg0JDg6mZcuWJCUl5UGWYiV7nusAO3fu\nxMfHhyVLlrgxO3EVe+Y9ISGBxo0bExQUpGMYC4g/mvczZ87wyCOP0KhRI4KCgpg9e/bvBzQcNGHC\nBCMqKsowDMOYMmWKMXHixNuuOXbsmFGjRg3j2rVrhmEYRt++fY3Zs2c7OqTkA/bMu2EYxqBBg4xZ\ns2YZhmEYWVlZRnp6uttyFGvZO+eGYRj//ve/jQEDBhjdunVzV3riIvbM+48//mjs3bvXMAzDuHTp\nklGrVi3j0KFDbs1TnJOdnW0EBAQYx44dMzIzM42GDRveNodxcXFG586dDcMwjG3bthmhoaF5Tepc\nzwAABoBJREFUkapYyJ5537Jly82f3atWrdK8ezh75vzGdWFhYUbXrl2NxYsX50GmYiV75v38+fNG\nvXr1jNTUVMMwDOP06dN5kapYyJ55//vf/268+OKLhmGYc37//fcbWVlZvxnT4ZUGy5cvZ/DgwQAM\nHjyYpUuX3nZNiRIl8PX1JSMjg+zsbDIyMqhcubKjQ0o+YM+8X7hwgU2bNjF06FDA7IlRsmRJt+Yp\n1rFnzgFOnDjBypUreeqppxw+9lDyD3vmvUKFCjRq1AiAYsWKUbduXX744Qe35inO2bFjBw888ADV\nq1fH19eX/v37s2zZsl9d88vvhdDQUNLT0zl16lRepCsWsWfemzdvfvNnd2hoKCdOnMiLVMUi9sw5\nwH/+8x/69OlDuXLl8iBLsZo9875gwQJ69+5982S8smXL5kWqYiF75r1ixYpcvHgRgIsXL1KmTBl8\nfH673aHDRYNTp07h5+cHgJ+f3x1/gbj//vsZN24cVatWpVKlSpQqVYoOHTo4OqTkA/bM+7FjxyhX\nrhxDhgzhwQcfZPjw4WRk6EgXT2XPnAM8//zzREdH4+XlCSdqyx+xd95vOH78OHv37iU0NNQd6YlF\n0tLSqFKlys33/f39SUtL+8NrdAPp2eyZ91+aNWsWXbp0cUdq4iL2PteXLVvGM888A4DNZnNrjmI9\ne+b98OHDnDt3jrCwMJo2bcq8efPcnaZYzJ55Hz58OAcPHqRSpUo0bNiQadOm/W7M3z09oWPHjpw8\nefK2///666//6n2bzXbHf1iOHDnC22+/zfHjxylZsiSPPfYYMTExPPHEE7+blOQtZ+c9OzubPXv2\nMGPGDEJCQhgzZgxTpkzh1VdfdVnO4hxn5zw2Npby5cvTuHFjEhISXJWmWMzZeb/h8uXL9OnTh2nT\nplGsWDHL8xTXsfem4H9XD+lmwrPdzfxt2LCBDz/8kM2bN7swI3E1e+b8xu9rNpsNwzC0arAAsGfe\ns7Ky2LNnD+vWrSMjI4PmzZvz0EMPERgY6IYMxRXsmffJkyfTqFEjEhISOHLkCB07dmT//v0UL178\njtf/btFgzZo1v/mYn58fJ0+epEKFCvz444+UL1/+tmt27dpFixYtKFOmDACPPvooW7ZsUdEgn3N2\n3v39/fH39yckJASAPn36MGXKFJflK85zds63bNnC8uXLWblyJdeuXePixYsMGjSIuXPnujJtcZKz\n8w7mLxu9e/cmMjKSnj17uipVcZHKlSuTmpp68/3U1NSbS1R/65oTJ05oq6GHs2feAZKSkhg+fDjx\n8fGULl3anSmKxeyZ8927d9O/f3/AbJK2atUqfH196d69u1tzFevYM+9VqlShbNmyFClShCJFivDw\nww+zf/9+FQ08mD3zvmXLFl5++WUAAgICqFGjBsnJyTRt2vSOMR1eR9y9e3fmzJkDwJw5c+74y2Kd\nOnXYtm0bV69exTAM1q5dS7169RwdUvIBe+a9QoUKVKlShZSUFADWrl1L/fr13ZqnWMeeOZ88eTKp\nqakcO3aMjz/+mHbt2qlg4OHsmXfDMBg2bBj16tVjzJgx7k5RLNC0aVMOHz7M8ePHyczM5JNPPrnt\nBqF79+43n8/btm2jVKlSN7euiGeyZ96///57Hn30UebPn88DDzyQR5mKVeyZ86NHj3Ls2DGOHTtG\nnz59eO+991Qw8HD2zHuPHj1ITEwkJyeHjIwMtm/frvs1D2fPvNepU4e1a9cC5pbU5ORkatas+dtB\nHe3KePbsWaN9+/ZGYGCg0bFjR+P8+fOGYRhGWlqa0aVLl5vXRUVFGfXq1TOCgoKMQYMGGZmZmY4O\nKfmAvfO+b98+o2nTpkZwcLDRq1cvnZ7gweyd8xsSEhJ0ekIBYM+8b9q0ybDZbEbDhg2NRo0aGY0a\nNTJWrVqVl2mLA1auXGnUqlXLCAgIMCZPnmwYhmHMnDnTmDlz5s1rRo0aZQQEBBjBwcHG7t278ypV\nsdAfzfuwYcOM+++//+ZzOyQkJC/TFQvY81y/4cknnzQ+++wzd6coLmDPvEdHR9+8X5s2bVpepSoW\n+qN5P336tBEREWEEBwcbQUFBRkxMzO/GsxmGNiyJiIiIiIiIyO3U5lxERERERERE7khFAxERERER\nERG5IxUNREREREREROSOVDQQERERERERkTtS0UBERERERERE7khFAxERERERERG5IxUNRERERERE\nROSOVDQQERERERERkTv6/9c86yG+HiwiAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10a99aa10>"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "Conduct a simulation to compare the various bootstrap confidence interval methods. Let $n=50$ and let $T(F)=\\int(x-\\mu)^3dF(x) / \\sigma^3$ be the skewness. Draw $Y_1, Y_2, ... Y_n \\sim N(0,1)$ and set $X_i = e^{Y_i}$.  Construct the 95 percent intervals for $T(F)$ from the data $X_1, X_2, ... X_n$. Repeat this whole thing many times and estimate the true coverage.\n",
      "\n",
      "First, let's write a function that let us estimate the skewedness."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Of course something has been written already, so we could use:\n",
      "# from scipy.stats import skew as skewness\n",
      "def skewness(sample):\n",
      "    mu = mean(sample)\n",
      "    third_moment_about_mean = sum((sample-mu)**3)/float(len(sample))\n",
      "    sigma3 = std(sample)**3\n",
      "    return third_moment_about_mean/sigma3\n",
      "\n",
      "def boot(sample_data, statistic_func, N=None, B=1000):\n",
      "    N = N or len(sample_data)\n",
      "    bootstrap_data = []\n",
      "    for i in range(B):\n",
      "        resampled_data = random.choice(sample_data, size=N, replace=True)\n",
      "        resampled_stat = statistic_func(resampled_data)\n",
      "        bootstrap_data.append(resampled_stat)\n",
      "    return bootstrap_data\n",
      "\n",
      "def normal_interval(alpha, theta, boot_data):\n",
      "    se = std(boot_data)\n",
      "    z = stats.norm.ppf(1-alpha/2)\n",
      "    return (theta - se*z, theta + se*z)\n",
      "\n",
      "def pivot_interval(alpha, theta, bootdata):\n",
      "    return (2*theta - percentile(bootdata, 100*(1-alpha/2.0)), 2*theta - percentile(bootdata, 100*alpha/2.0))\n",
      "\n",
      "def percentile_interval(alpha, bootdata):\n",
      "    return (percentile(bootdata, 100*alpha/2.0), percentile(bootdata, 100.0*(1-alpha/2.0)))\n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "random.seed(0)\n",
      "\n",
      "runs = 500\n",
      "alpha = 0.05\n",
      "true_value = 6.1848771386325572\n",
      "func = skewness # variables and arguments can be set to functions\n",
      "\n",
      "N = 50\n",
      "B = 10000\n",
      "\n",
      "normal_capture_count = 0\n",
      "pivot_capture_count = 0\n",
      "percentile_capture_count = 0\n",
      "\n",
      "for i in range(runs):\n",
      "    data = exp(normal(size=N)) # the default values for mean and variance are the standard normal\n",
      "    theta = func(data)\n",
      "    bootdata = boot(data, func, N, B)\n",
      "    \n",
      "    a,b = normal_interval(alpha, theta, bootdata)\n",
      "    if a <= true_value <= b:\n",
      "        normal_capture_count += 1\n",
      "        \n",
      "    a,b = pivot_interval(alpha, theta, bootdata)\n",
      "    if a <= true_value <= b:\n",
      "        pivot_capture_count += 1\n",
      "        \n",
      "    a,b = percentile_interval(alpha, bootdata)\n",
      "    if a <= true_value <= b:\n",
      "        percentile_capture_count += 1\n",
      "        \n",
      "    if (i+1)%50 == 0:\n",
      "        print \"Finished run %d\" % (i+1)\n",
      "        "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Finished run 50\n",
        "Finished run 100"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Finished run 150"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Finished run 200"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Finished run 250"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Finished run 300"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Finished run 350"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Finished run 400"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Finished run 450"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Finished run 500"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "normal_capture = float(normal_capture_count)/runs\n",
      "pivot_capture = float(pivot_capture_count)/runs\n",
      "percentile_capture = float(percentile_capture_count)/runs\n",
      "        \n",
      "print \"Supposed: %.3f\\nNormal: %.3f\\nPivot: %.3f\\nPercentile: %.3f\" % (1-alpha, normal_capture, pivot_capture, percentile_capture)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Supposed: 0.950\n",
        "Normal: 0.400\n",
        "Pivot: 0.368\n",
        "Percentile: 0.300\n"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 279
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 279
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "_mu = sqrt(e)\n",
      "_sigma = sqrt(e**2-e)\n",
      "\n",
      "_num = _mu**9 - 3*_mu**5 + 2*_mu**3\n",
      "_den = _sigma**3\n",
      "\n",
      "_num/_den"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "6.1848771386325572"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in range(1,9):\n",
      "    print i, skewness(exp(normal(size=10**i)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 1.63820639099\n",
        "2 2.70077782089\n",
        "3 5.70857376582\n",
        "4 5.82184964254\n",
        "5 6.94182273591\n",
        "6 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "6.20437371151\n",
        "7 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "6.50476339665\n",
        "8 "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "6.18860954964\n"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}