All of Stats Chapter 8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.3**. (Computer Experiment)  Generate 100 observations from a $N(0, 1)$ distribution.  Compute a 95 percent confidence band for the CDF $F$.  Repeat this 1000 times and see how often the confidence band contains the true function.  Repeat using data from a Cauchy distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import norm, cauchy, uniform\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm.notebook import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1a4a5e48a60>"
      ]
     },
     "execution_count": 167,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABb5klEQVR4nO2dd3yN1x/H3ydLECP2JmLG3jP2pokRe3crbZXSqp9Wl1KqqrSqrZq1Y29iBjVTI0FCkBiRRWSv8/vjEImMexPRDOf9ej2ve+/znHOe85B87sn3fIeQUqLRaDSanI9JVk9Ao9FoNJmDFnSNRqPJJWhB12g0mlyCFnSNRqPJJWhB12g0mlyCWVbduFixYrJSpUpZdXuNRqPJkZw9ezZASlk8pWtZJuiVKlXizJkzWXV7jUajyZEIIW6ldk2bXDQajSaXoAVdo9Focgla0DUajSaXkGU29JSIiYnB19eXyMjIrJ6KJgUsLS0pV64c5ubmWT0VjUaTAtlK0H19fSlQoACVKlVCCJHV09EkQkpJYGAgvr6+2NjYZPV0NBpNCmQrk0tkZCRFixbVYp4NEUJQtGhR/deTRpONMSjoQoglQogHQohLqVwXQoj5QggvIcQFIUTDF5mQFvPsi/6/0WiyN8aYXJYCC4DlqVzvDlR9cjQDfn3yqtFoNDmWuLg4oqOjiY6OJiYmJsX3Gb3WqlUrunTpkulzNijoUsojQohKaTRxBJZLlVj9pBCisBCitJTyXmZNUqPRZD6h0aH84/sPkpxdEyE6Kprw0HDCQ8MJexymXkPDEs4l+fw4+bmY6BhiY2KJjY0lNiaWmOgY4mLjiI+Pf2lzfvvDt7NG0I2gLOCT6LPvk3PJBF0I8TbwNkCFChUy4daZj5WVFaGhoVk9DY3mpfIo8hFtlrbhgt+FrJ5K6kQDj4CQJ6+J34cA4UAUEGfEWAKwBPI8d1gBpmkcJsZdE0JgEWtBnrg85InNQ56YPFjGWCYceaPykjc6L/mj8pM3Mi/FTIu96L9OimSGoKdkWE3xK19KuRhYDNC4ceOcvSzQaHIwi88u5oLfBZb1XkZl68pZMgcpJYEPAvHx9sHH2wdfb198vH2453OPB3cf8PjR4yTthRAUKV6EEmVKUKJ6CQoXLYxVASvyWeUjv1V+8hfMr16t8pO/QP5n5wvkxyKPRZp7QPGx8cQ9jiMmJIbYkFhiH8USGxJLzKMnnxMdMQ+Tn4sLT/tbRZgKzAqZYVZQHXal7TLl3/B5MkPQfYHyiT6XA+6+8Kjjx4Ob2wsPk4T69WHevEwZatSoURQsWJAzZ85w//59vv/+e5ycnDJlbI3mZbPefT2NSjdiRL0R/9k9//33XzZt2oSHhweenp54enom+Ws4T5482NraUtO2Jl3adaFChQqUL18+4ShbtiwWFhYZureUkjC/MAKuBODv4U/AlQACPAIIvBZIRGAE0aHRafYXJgLLwpYJR8HCBbGsYImltWWS8yke1paY5zP/T5wKMkPQtwLjhBBrUJuhj14V+/m9e/c4duwYV65cwcHBQQu6Jkew5PwSTt89zbyu817qfaSU3LlzB2dnZ/766y/c3NwwMTHBxsaGatWqYW9vT7Vq1ahWrRpVq1alfPnymJqavtA94+Pieej9MIloP32NfPjM5dY8vznFahSjQqsK5C+ZP00xtixsiYVV2iv87IJBQRdCrAbaAcWEEL7AF4A5gJRyEbAT6AF4oaxaozNlZpm0kn6Z9O7dGxMTE+zs7PDz88vq6Wg0BnHxdmHszrF0tOnIuKbjMm1cKSU3btzg/PnznDt3jnPnznH+/HkePHgAQKNGjfj5558ZPHgwRYsWfeH7xUbGpijagdcCiYt+Zv7IXzI/xWsWp9agWhSvWZxiNYtRrEYxCpYrmCUCHR8PDx6AqSkUTzEB7othjJfLYAPXJTA202aUg8iTJ0/Ce/XPoNFkT64FXuPT/Z+y6comKltX5u9+f2Nq8mKrYQAvLy8mTZrEwYMHefToEQBmZmbUqlWLHj160LBhQ9q3b0/t2rVf6D5x0XH4/uPLzYM38XbxxveEb4JwCxOBdWVritUoRpXuVRJEu1iNYuS1zvvCz2gs8fHg7w8+PuDrm/LrnTsQEwOffQbffpv5c8hWof8ajSZzkVIyed9k5v0zD0szS75p/w0ftfiIfOb5MjxmWFgYV65cYefOncyYMQMLCwuGDBlCw4YNadiwIbVr106y2MkI8bHx3D17N0HAfVx9iAmPAQGl6peiybgmlGtejuI1i1OkShHMLF+ulMXHQ0CAEuXUBPvOHYh+zhRvYQHlykH58tCqlXotXx6aN38589SC/hzh4eGUK1cu4fOECROYMGFCFs5Io8k4v5z+hTkn5jC6/mi+6/gdJa1KGt03NDQUDw8P3N3duXz5Mu7u7ri7u3Pz5s2Ev0idnJz46aefKFOmzAvPNfR+KO4b3PHa7cWtI7eIfqzUsUTtEjR4owE2HWyo2KYieYtk7qpbSrWyTm1V7eurjufF2tz8mVi3aPFMrJ+eK1dOmVX+S8uOFvTnMDaYYOnSpUk+a991TXbjasBVJu2bRLcq3fjT4U+jbcbHjx/nzTffxMPDI+GchYUF1atXp2nTpowePRo7Ozvq1KlDtWrVXmiO4QHhuG905/Lay9w6fAsZLylarSh1htbBpr0NldpVIn+J/C90j4R7hYOHB7i7w+XL6tXDA27fTlmsy5ZVwtysGTg5pSzWJtkqG5YWdI0mVxIXH8foLaOxNLNkicMSo8Q8Pj6ehQsXMmHCBCpUqMA333yDnZ0ddnZ22NraYmaWOXIR+TASj00eXF57mRv7byDjJEWrF8X+f/bUHlib4nYvtlsYFgZXrjwT7cuX1XHzplqNgxLsatWUJ3OfPs+E+qlYlyiR/cTaGLSgG+Dbb79l/fr1Sc7179+fqVOnZtGMNJrUCY8J58itI6y6uIoTvidY0WcFpQuUNtjv4sWLjBkzBldXV3r16sWKFSsoXLhwps0rPjYez52enP/zPF67vYiLjqOwTWFaTmpJ7YG1KVmvZLq9TsLCkq+4UxLu6tWhSRMYNQrs7KBWLahSRV3LbWhBN8DUqVO1eGuyLfEyngt+F9h7fS97r+/l6O2jRMdFY2lmyXuN32NonaEGx9i0aRP9+/encOHCLFmyhJEjR2KSScvTx3cfc+6Pc5z7/RwhviFYlbaiybgm1B5YmzJNyqRLxCMi4PBh2LMH9u5VAv4UCwsl3E2bwujRz4Tb1jZ3CndqaEHXaHIos11nM+fEHB6EKV/vOiXq8H7T9+li2wX7CvbkNTdu83Dz5s0ULVoUd3f3TPERBwi8FojrbFcuLL9AXHQctl1s6Ta/G9V6VcPU3Dh3SSmVaO/erUT8yBGIigJLS2jTBgYNUqL9VLgzySKUo9H/BBpNDmXuybmUtirN7M6z6Vy5s1Gmlec5cuQI27Zto1GjRpki5nfP3sV1pivuG90xy2NGgzcb0OKjFhSpUsSo/sHBsG+fEvA9e5QrIKgV93vvQdeuSszz/nfu5TkKLegaTQ7EL9SP+6H3mdxycobzsWzduhUnJydsbW1ZtGhRhucipeTmoZsc++4YN/bdIE/BPLT+tDXNPmyGVUkrg/19fGDLFti8GQ4dgrg4KFwYOndWAt6li9qs1BhGC7pGkwNx9XEFoGnZphnqHxsby/jx46lZsyaHDx/O8Aaot4s3LlNd8D3pS/6S+ek0qxON3mmEZSHLVPtICZcuPRPxs2fV+Zo1YfJkeO01tYmpTSjpR/+TPcfNmzfp1asXly49q7g3ffp0rKys+Pjjj7NwZhrNMw56HyS/ef4MC7qzszPe3t5s2rQpQ2Lu7+7P7vG7ubHvBgXLFaTnrz2pP6p+mhGbAQEwaxZs2gTXr6uAm+bN1TlHR7WpqXkxtKBrNDmQ2yG3sS1ii7lpxlw4fvvtNypVqoSDg0O6+kWHRnP468OcnHsSiwIWdJnbhSZjmhgMvT98GIYOBT8/ZUp5uhIvnX6zvyYNsq2gj989Hrf7bpk6Zv1S9ZnXbV6mjtmuXTuaNWvGwYMHefjwIX/++Sf29vaZeg+N5nkeRj6kgEWBDPWdO3cuLi4uzJo1y2j3RCklVzZdYfeHuwnxDaH+6/XpNLMT+YunHcUZFwdff60OW1s4dQoaNMjQtDVGkG0FPScRGxvLqVOn2LlzJ19++SX79+/P6ilpcjGB4YGc9D3JmMZj0tXv/v37/PDDD8yZMwcnJycmTpxoVL8gryB2vb8Lr91elKxbEqe1TpRvmfYupY8PrFkDy5cre/mIEbBgARTI2HeQxkiyraBn9kraWFILdEgrAKJv376Ayvl88+bNlzEtjQaAu4/v8uHuD4mOi2Z0feNKD3h7ezN79myWLFlCTEwMI0aMYPHixQaLScRGxnJs5jGOzTyGqYUpXX/sStNxTTExS3lVHxgIGzbA338rn3FQeVBWr1Y+45qXT7YV9KyiaNGiBAcHJzkXFBSEjY1Nqn2epgo1NTUlNjb2pc5P82oSHhPOD8d/YJbrLKLjopnWZhr1StUz2M/NzY1WrVoRExPDyJEjmTx5MlWrVjXYz2u3FzvH7ST4ejC1B9Wmyw9dKFAm5eX1jh2waJEKAIqNVd4qX38NgwcrM4vmv0ML+nNYWVlRunRpDhw4QMeOHQkKCmL37t18+OGHWT01zStKUEQQzf5ohleQF/1q9mNWp1nYFjGslP7+/vTu3Rtra2uOHTtGpUqVDPaJfBTJnvF7cFvqRtHqRRm+fziVO6ZeRHr5chg5UvmJT5gAQ4ZA3br/bcpYzTO0oKfA8uXLGTt2bIKN8YsvvsBWLzU0WYCUktFbRnPr4S32DttLZ9vORvWLi4tj6NCh3L9/n6NHjxol5reO3mLTsE2E+IZgP9WeNtPaYJYndYnYtg1efx06dYLt2+EFa1poMgEt6ClgZ2fHwYMHjWp76NChhPfFihXTNnRNprLefT1br27lp24/GRTz+Ph4Ll26xNGjR9m+fTv79u1j8eLFNGnSJO1+sfEc/uowR789SmGbwrzu+jrlmpdL1k5K5T++d686du1SHivOzlrMswta0DWabMzZu2exMLVgbJPkZXtjYmI4c+YMR48e5ciRI7i6uvLw4UMAypYty+eff86bb76Z5vghviFsHLyR28duU29kPbr/3J08BZ6pc3AwuLio/Cp794K3tzpfsaJanX/9tfZcyU5oQTeSsWPH4urqmuTchx9+yOjRxnkaaDTGEhcfx5FbR1jvvp41l9ZgU9gmxYLOr732Gnv27AGgevXqODk5YW9vj729PZUqVTKYmtZrjxfOQ52Ji4qj76q+1BlSB1Abm3PmqND8U6dUPU0rK+jQASZOVLlVqlTRdvLsiBZ0I1m4cGFWT0GTi5FScuTWEdZdXsdGj434hfmRzzwfvar1YlyTccnaR0dHc+jQIYYOHcrcuXMpUaKE8feKlxz55giHph+iRO0S9F/fn2LViyVcnzJFCXrTpjB1qhLwZs1erbziORUt6BpNNmDKgSnMcp1FXrO89KzWkwF2A+hRtQf5LVKOxDx79ixRUVE4OjqmS8yjQqLYNHwTV7depe7wuvRa1AvzfM+U+u+/lZi/9x7oNYyRSAmPHyv7VEpHUFDyc8OHw/vvZ/pUtKBrNFnM+XvnmXN8DsPrDufXnr+mKuKJ2bVrFyYmJnTq1Mno+wTfCGb1a6sJuBpAt/ndaDquaRKzzM6d8MYbYG8PP/6YoUfJuUgJoaHGC/LzR1xc6mObmoK19bOjaFEoWPClPIYWdI0mC3kU+Yjhm4ZTNF9Rfur2k1FiDnDp0iVq1KiBtbW1Ue1vHb3Fur7riI+LZ/je4dh0SBoot3698iGvU0d5rVhYpPtRsidSKsH18QFfX3Ukfn/njgpxDQ5Wmwep8bwoW1tD5crJzxUpkvycldV/tuGgBT0F7t+/z/jx4zl9+jR58uShUqVKzJs3j3r16lGjRg0iIyMpUKAAY8eOZeTIkQAsXbqUSZMmUbZsWQDq1q3L8uXLs/IxNNmc6Lho+q3rx9XAq+wZtgfrvMaJM6joZWPT3rotdWPb29uwtrFm8PbBFK2atDLRjh0qNL9lS+VPXqhQep4iC5FSiXFKQp34fURE0n6mplCmDJQrB7VrQ7FiqYvx06NAgRyxC6wF/TmklPTp04eRI0eyZs0aQIVP+/n5YWtry/nz5wG4ceMGffv2JT4+PsHTZeDAgSxYsCDL5q7JWXzv+j0HvA+w1HEpHWw6GN3vzp07uLq6GoxellJy8PODHP3mKJU7VcZpnRN5rZPXbtu6VYn47t2Q37g/EP57Hj8GNzc4d+7Z4eUFkZFJ25maQtmySqwbNAAHB/W+XDkVzlquHJQqpdrlQrSgP8fBgwcxNzfn3XffTThXv379ZAFDlStXZu7cuUycOFG7LmrSjX+YP9+7fk+fGn0YWX+kUX3i4+PZuHEjX3zxBXFxcYwZk3q2xbiYOLa9tY1/l/1Lgzca0PPXnqkWZw4MVBqXbcQ8OBjOn38m3GfPgqenWpGDSqLesKGqT/dUpJ++liyZa8XaGLKtoO8ev5v7bvczdcxS9UvRbV63NNtcunSJRo0aGTVew4YNuXLlSsLntWvXcuzYMUD7qGtSZ9vVbXyw+wPCY8L5tsO3RvXZuXMnn376KRcvXqRGjRps3rw51XQUsZGxrHNah+cOT9pOb0vbz9um6pPu5QUHDkDr1hl+nBcjPh48PMDVFY4fV69eXs+uV6yoxHv4cPXaoIGuipEG2VbQcwLy6YrhCdrkokmLSw8u8dmBz9h2bRt2xe1wGelCzeI1Dfbbtm0bjo6OVKlShVWrVjFw4MBUU99Gh0azxnEN3ge96bmoJ43faZzquI8fQ+/eYGIC8+dn9KnSSXi4ilZydVXHiRPwJLqV4sWhVSvlatOokRLvYsXSHE6TFKMEXQjRDfgJMAX+kFLOfO56IWAlUOHJmHOklH+9yMQMraRfFrVq1WLDhg1GtT1//jw1axr+hdS8mkgpufjgIhvcN7DBfQMeAR7kN8/P952+Z3zz8UaVj7t27RpDhgyhYcOGHDlyhHz58qXaNvJRJH/3+Bvfk770Wd6HusPqpjn255+rxfHevZBGdujMITAQvvkGfvkFoqPVOTs76N9fiXirVirXbg7YeMzOGBR0IYQpsBDoDPgCp4UQW6WU7omajQXcpZSvCSGKA1eFEKuklNEvZdYvkQ4dOvDZZ5/x+++/89ZbbwFw+vRpwsPDk7S7efMmH3/8Me+/hOAATc7GP8yfH0/+yAb3DXgGeWIiTGhTsQ1jm4zFyc6JklYljR5r3bp1hIaGsmnTpjTFPDwwnJVdV+J3wQ+ndU7Y9bNLc9y4OFV4ok8f6NjR6Omkn6goVarom28gJETl2nVyUtWhixR5iTd+NTFmhd4U8JJS3gAQQqwBHIHEgi6BAkIZ6qyAICBHVnoQQrBp0ybGjx/PzJkzsbS0THBbvH79Og0aNEhwW3z//fe1nVyTjAl7J7D64mra27RnYouJ9K7RO10iDhAQEMCOHTtYs2YNZcuWpXz51Eu+hfqFsqLTCgI9Axm0eRBVexguYOHqqgo29++frmkZh7e3Wvbv26cM9A8fQvfu8P33yk1Q89IwRtDLAj6JPvsCzZ5rswDYCtwFCgADpZTxmTLDLKBMmTKsW7cu2fmI5/1ZEzFq1ChGjRr1EmelySmcunMKh+oOOA90Tlc/Ly8vtmzZwtatWzl27Bjx8fGULVuWSZMmpdonxDeE5R2XE+IbwpAdQ9IsRpGYHTtUbpYePdI1xVQmEaJSMj4V8aebmuXLQ79+qnTRS/0zQPMUYwQ9JaOWfO5zV8AN6ADYAvuEEEellCFJBhLibeBtgAoVKqR7shpNdudOyB08Az0ZVW9UuvrNnDmTKVOmACooberUqTg6OtKwYcNUPVSCbwSzvONyIoIiGLZ3GBVaGf87tWuX8mx54dS3R48qu01goPJ7bN8ePvhAZfSqVk3bxP9jjBF0XyDx33vlUCvxxIwGZkrl9uElhPAGagCnEjeSUi4GFgM0btz4+S8FjSbHs+bSGiQSJzsno/ucO3eOadOm0adPH3744Yc069c+JfhGMEvbLiUmPIYRB0ZQpnEZo++3YQNcvAg//2x0l+RICatWKY8UGxuVO6BVq1yUMyBnYoygnwaqCiFsgDvAIGDIc21uAx2Bo0KIkkB14EZmTlSjyQns8NxBvZL1qFrUsB0b4O7duzg5OVGiRAn+/PNPo3KzPLr9iGUdlikxdxlBqXqljJ6fv7/KpNioESSKnTPMw4fK3fDkSeVq+M8/KgCobVuV/EVvcGYLDAq6lDJWCDEO2INyW1wipbwshHj3yfVFwNfAUiHERZSJ5hMpZcBLnLdGk+2Il/GcuXuGYXWHGdU+KCiIzp074+/vj4uLi1FiHno/lOUdlxP5MJKRLiPTJeY7dsCYMfDoESxdCmZp/fZfuwaHDikBP3lS+TeCMqHUqqVs461aKfu4rj+XbTDKD11KuRPY+dy5RYne3wW6ZO7UNJqchX+YP4+jH2NXPG2Xwad8/PHHeHp6snfvXoN1PwEigiJY0WUFj+8+Zvi+4ZRuaFzEpJ8ffPghrF2rXL/XrjXgbLJsGYwercwqxYopF8OhQ9VrkyYvLfWr5sXRkaIaTSYRHqNiFSxMDduRt23bxl9//cWkSZNo166dwfZRj6NY1WMVgVcDGbJjCOVbpu7GmBh/f6hXT1lHvvoKPvnEgJl75Uol5h07wq+/6mCfHIYWdI0mE3C778bADQMxESbULpH68vfq1at8/PHHbN++nerVq/O///3P4NixkbGs7b2Wu2fuMmDjACp3Ms41EZSABwYqq0maKYri4pTif/21Kh66ZQukEcikyZ6YZPUENEmpVKkSAQF6+yEnsfzf5TT/ozmPox5zYMQBWpZvmeT648eP2bBhA8OGDaN27docPnyYWbNm4ebmRkED5ov4uHichzrj7eJN76W9qeFYw+h5ubrCX3+pws5pinlwMPTqpQR95EjYtk2LeQ5Fr9AzkdjYWMzS3GnKWUgpkVJiYqK/91Nj/439vL7lddpWasuafmsonr84oLxXtm7dypYtW3BxcSE6OpoiRYrwzjvvMG3aNEqWNBw5KqVk57ideDh70HVeV4O5WRITF6dKVpYtC9OmGWj8/fcqIOi33+Ctt7SJJQeTbdVn/PjxuLm5ZeqY9evXZ968eWm2uXnzJt27d6d169YcP36csmXLsmXLFq5evcq7775LeHg4tra2LFmyBGtra9q1a0fLli1xdXXFwcGBbdu20aBBA86ePYu/vz/Lly/nu+++4+LFiwwcOJBvvvkGgN69e+Pj40NkZCQffvghb7/9tsH537x5k169enHp0iUA5syZQ2hoKNOnT6ddu3bUr1+fU6dOERISwpIlS2jatCnTp0/n+vXr3LlzBx8fHyZPnpyQo2b27NmsW7eOqKgo+vTpw5dffpnw/O3bt+fEiRNs3ryZihUrvtg/fC7lSsAVBqwfQI1iNdg8cDMF8qgonQULFiTk+LG1tWXcuHE4OjrSsmXLdH3hu37vytlFZ2n1SSuaf9jc6H43bqgCz+fPq3wtBvOc37ypfMmN+BnUZG+yraBnJZ6enqxevZrff/+dAQMGsHHjRr7//nt+/vln2rZty+eff86XX36Z8OXw8OFDDh8+DKjNLgsLC44cOcJPP/2Eo6MjZ8+epUiRItja2vLRRx9RtGhRlixZQpEiRYiIiKBJkyb069ePokWLpjErw4SFhXH8+HGOHDnC66+/niD8Fy5c4OTJk4SFhdGgQQN69uzJpUuX8PT05NSpU0gpcXBw4MiRI1SoUIGrV6/y119/8csvv7zQfHIzF/0u0mlFJ8xNzdkyaEuCmAOsXLmSOnXqsGbNGmrWrJlqpGdaXFpziQOfHqD2oNp0nJF22Hx8PJw5o8zeW7bA5cvqfN++MHCgoQe5CHv2QNOm6Z6jJvuRbQXd0Er6ZWJjY0P9+vUBaNSoEdevX+fhw4e0bdsWgJEjR9I/UVajgc/91jg4OABQp04datWqReknCfkrV66Mj48PRYsWZf78+WzatAkAHx8fPD09X1jQBw8eDECbNm0ICQnh4ZM8046OjuTNm5e8efPSvn17Tp06xbFjx9i7dy8NGjQAIDQ0FE9PTypUqEDFihVp3tz4FeGrhJSSLVe38MbWN8hrlpcDIw5gW+RZoYnQ0FDc3NwYO3YsdnbGuS8+z60jt9g8cjMV7CvguNQRYZLyF8K5c7B4sTJ5372r8prb28PcueDoqGoYp4mnJ3TurOzlv/6aoblqshfZVtCzkjyJAiVMTU0ThDE18j/3N+3T/iYmJknGMjExITY2lkOHDrF//35OnDhBvnz5aNeuHZHP10ZMATMzM+Ljn+U8e77P8yvBp59TOi+lZMqUKbzzzjtJrt28eTPZ82gU5++dZ8LeCRy6eQi74nZsG7yNytZJVdPZ2ZmoqCh69+6doXsEXAlgTe81FLYpzKDNgzDLk/qv6MiRyrzSvbsS8B49IF1rggkTICZGBRC99ITomv8CvdtlBIUKFcLa2pqjR48CsGLFioTVekZ49OgR1tbW5MuXjytXrnDy5Emj+pUsWZIHDx4QGBhIVFQU27dvT3J97dq1ABw7doxChQpR6En59i1bthAZGUlgYCCHDh2iSZMmdO3alSVLlhAaGgqowsMPHjzI8DPlZiJjI3lvx3s0WtyIi34XWdhjIf+++28yMQdYv349lSpVonUGarqF+oWyqscqTM1NGbpzKHmLJC/o/JSYGLhyRQUMbdigKrSlS8zv3IGdO+Gdd6CG8Z4zmuyNXqEbybJlyxI2RStXrsxff2W8IFO3bt1YtGgRdevWpXr16kabN8zNzfn8889p1qwZNjY21HjuF9Ha2pqWLVsmbIo+pWnTpvTs2ZPbt28zbdo0ypQpQ5kyZfDw8KBFixYAWFlZsXLlylRLm72q+Ib40m9dP07dOcWHzT5kervpFLYsnGJbKSUnTpzAwcEh3Xbz6LBoVr+2mtD7oYw6PArrymmnAbh6FWJjVeRnhlizRhnfdT7/3MVT17T/+mjUqJF8Hnd392TnNMbRtm1befr06WTnv/jiCzl79uxMu8+r9H/kds9NlpxdUlrNsJLO7s4G2/v5+UlAzps3L133iY+Ll2v7rpXTxXTpsdnDqD6//SYlSHntWrpu9YxmzaRs0CCDnTVZCXBGpqKreoWu0aTABb8LdFzekbzmeTk54iS1StQy2CcsLAyAvHlTN5WkhMs0FzycPegyt4vRgUP//KPSrFSpkq5bKZYvVwPMnZuBzprsjBb0bEhgYCAdU6jwcuDAgVQ9YQ4dOpTi+enTp2fizF4N/EL96LS8E5ZmlhwaeSiJF0tqnDx5kqFDh2JiYkLtdJRZu/j3RY7NOEbDtxrSfLzxnkX+/ipoKN0ekUeOqLy57duryCNNrkILejakaNGimR5UpTGeX8/8in+4PxfevWCUmP/0009MnDiRcuXKcfjwYVq2bGmwD8CdU3fY8voWKrapSI8FPdJld3/82IiAocRcvw5TpqhCFBUrwt9/G8ifq8mJaC8XjSYR8TKe387+Ro+qPahTso7B9mvWrGH8+PH06tULNzc3o71bHt99zJreayhQugADNg7A1ML4zejQUGUxeRIqkTZSKiGvWVMlRJ8+HS5dglLG51HX5Bz0V7RGk4g7IXe4H3qf16q9ZrCti4sLo0aNwt7enrVr1yaJOUiL2KhY1vVbR1RIFMNODCNfMeMTYUkJ8+dDRAQMGmREh99+g5kzlV/jrFlQ2rgc6pqciRZ0jSYRF/wuAFC1SNol5DZv3szAgQOpWrUqzs7ORou5fJJwy/ekL/039KdkHcNJup7i4aEqDh0+rNKVt2ploIOnp0q12KWLKlqhk27lerTJRaN5gpSSb49+S2mr0rQo3yLFNsHBwUyePJl+/frRsGFDjhw5QrFixYy+x9nfznL+j/PYT7XHrp/xTuRff60KVVy4AL//Dnv3qlD/VImLUz7mFhawZIkW81cELejZDEP50K2srP7D2bxa7L+xnxO+J/iy3ZfkM09qBomKimLu3LnY2toyZ84cRowYwb59+yiSjuLIt11vs+uDXVTtUZV2X7Yzut+aNfD55yrZ1tWr8OabaYh5ZCTs3w9vvKESos+bp9xhNK8E2uSSieSmfOhxcXGvXNTopQcqO2U/u35Jzl++fJl+/fpx9epVunbtyqxZs6hXr166xg69H8p6p/UUrliYvqv6YmJq3Frq3j147z1o1kxVh0v24yWlSq+4d686Dh9Wom5urtLhjhiRrnlqcjbZVn3Gj4fM9tyrX18tWNIiO+dDT4yUksmTJ7Nr1y6EEPzvf/9j4MCBvPfee3Tr1g0HBwf69OmDtbU1S5Ys4c8//8Tb25tvvvmGlStXMn/+fKKjo2nWrBm//PILpqamWFlZMWHCBPbs2cMPP/yQoXwkOZmz986Szzwf1pbPwu7Xrl3LG2+8gZWVFbt27aJbt27pHjc+Np4NAzcQFRLF8H3DsSxsaXTf1atVQaGlS1MQ8x9/hNmzleqD8mR55x1lM2/TBvRfc68c2uSSAp6enowdO5bLly9TuHBhNm7cyIgRI5g1axYXLlygTp06fPnllwntn+ZDnzhxIkBCPvR3330XR0dHFi5cyKVLl1i6dCmBgYEALFmyhLNnz3LmzBnmz5+fcN5YnJ2dcXNz499//2X//v1MmjSJe/fu0aZNm4QkYnfu3MHd3R1QCbvs7e3x8PBg7dq1uLq64ubmhqmpKatWrQJUpGPt2rX5559/Xjkx3+W5i1UXV/Fuo3cT/MHnzJnDoEGDqFevHufOncuQmAPsn7KfW0du0WtxL0rULpGuvl5eUKhQKvmzZs9WGbn++ANu3QJ3d7Vi6dFDi/krSrZdoWdhOvQckQ/92LFjDB48GFNTU0qWLEnbtm05ffo09vb2zJs3D3d3d+zs7AgODubevXucOHGC+fPns2zZMs6ePUuTJk0AiIiIoEQJJTKmpqb069cvrdvmWsbuHEut4rX4tuO3AMycOZMpU6YwcOBAli9fjoWFRYbG9XD24MScEzR+rzF1hxpXQs7fX9nNV6yA06dTqT3x4IFamU+apOzlGg3ZWNCzkuyaDz0xKkdPcsqWLUtwcDC7d++mTZs2BAUFsW7dOqysrChQoABSSkaOHMl3332XrK+lpeUrZzcHFerv/dCbH7v+iKWZJceOHWPKlCkMGTKEZcuWZXhfJPBaIJtHbaZs07J0nds1zbaxseDsrER89271uX59+OEH5UKeDFdX9frki1mjAW1yMYrskg89MW3atGHt2rXExcXh7+/PkSNHaPpkKdeiRQvmzZtHmzZtsLe3Z86cOdjb2wPQsWNHNmzYkJD7PCgoiFu3bmX4WXI6QRFBjN8zHoD6peoD8OOPP1KkSBF+//33DIt5dFg06/qtw9TClP7r+6dZqAJUrYmBA1Ud0AkTVGW4p++LF0+hw6FDqtKQLh2nSYReoRtJdsiHnpg+ffpw4sQJ6tWrhxCC77//nlJPwrnt7e3Zu3cvVapUoWLFigQFBSUIup2dHd988w1dunQhPj4ec3NzFi5c+EoWgt58ZTPvbn+XwIhAvmj7BW0rtiU8PJzNmzfz0UcfkS+f8RGciZFSsuPdHTy4/IBhu4dRqEKhNNtfvKiKOr/9NvzyCxj1R9KtWyrVYgZNQZpcSmp5dV/2ofOh50xyy/+Rh7+HZDqy/qL68vy98wnn//33XwnINWvWZHjsU7+cktOZLg99dchg28hIKZs3l9LaWsqAgHTcpG1bKVu1yvAcNTkX0siHrk0umlcSn0c+ACzoviDB1AKqXB9A3brGbWA+z71z99gzfg9VulehzdQ2abaVEsaNg5MnVbFno/fEg4LgxAlo1ChDc9TkXrTJJRuSkXzomvRx+9FtAApZPjOH+Pn58f3339O3b19q1qyZ7jEjH0Wyvv968pfIT5/lfRAmaYfbr12rPA4/+wycnIy8iZ+fCtKIjoZRo9I9R03uJtsJupQy3fUYcxvZNR+6TMWzJify5/k/qVa0GrWKP6tEtGPHDkJDQ/niiy/SPZ6Ukq1vbOXR7UeMOjzKYAbF2FiYNk3lZ/nqKyNuEBGhfHlnzFCRoFOmGJk/V/MqYZTJRQjRTQhxVQjhJYT4NJU27YQQbkKIy0KIwxmZjKWlJYGBgblKOHILUkoCAwOxtDQ+yjG74hXkxQnfE7zT6J0kiwcfHx+EEMmKbxvD6V9O47HRgw4zOlC+Zfk020ZEwEcfqaCh6dON2AR9+BAaNFBL+Y4dVaj/jBk64ZYmGQZX6EIIU2Ah0BnwBU4LIbZKKd0TtSkM/AJ0k1LeFkKkLxzuCeXKlcPX1xd/f/+MdNe8ZCwtLSlXrlxWT+OFOXTzEAA9q/ZMOBccHMwff/xBnTp10h1EdO/cPfZO2EvVHlVpOTHtakWHDsFbbykxf/ddcHQ0MLiU8PrrquLQzp3QvXu65qZ5tTDG5NIU8JJS3gAQQqwBHAH3RG2GAM5SytsAUsoHGZmMubk5NjY2Gemq0RjNbq/dlLYqTbWi1RLOffDBB9y/fz8hetdYokKi2DBwA/mK56P3st4p2s0jIpQWr1oFmzZB5coqIWIK2yTJ2bhRdZo7V4u5xiDGmFzKAj6JPvs+OZeYaoC1EOKQEOKsEEKneNNkS0KiQth+bTv9avZLMLdERUWxZs0axowZQ+PGjY0eS0rJ9ne3E3wjmH6r+yWxm8fEqIjPkSOhZEm16enqqkzfFy8aKeYAZ86ozIkffJCex9S8ohizQk/JUPe8kdsMaAR0BPICJ4QQJ6WU15IMJMTbwNsAFSpUSP9sNZoXZNWFVUTFRTGkzpCEc4cPHyY2Ntbo4s5PcfvLjUurL9H+6/ZUtH8WmLV9u6otERCgEms5OcHgwdC+vZF1mePj4fhxWLdOLettbIyMNtK86hjz4+ULJN7lKQfcTaFNgJQyDAgTQhwB6gFJBF1KuRhYDNC4cWO986n5T4mX8fx48kcal2lM83LNiYqK4quvvmLWrFmUKVOGDh06GD2Wv4c/O8ftxKaDDa2nJM1MOX8+WFoqS0n37mBUdToplW/5unWwfj3cvasG6dED3n8/nU+qeVUxRtBPA1WFEDbAHWAQymaemC3AAiGEGWABNAN+zMyJajQvyp/n/sQzyJO1TmsJCwujRYsWXLp0idGjRzN37lwKFy5s1DgxETFsHLQRCysL+qzok6RYhZTw779Kh3v3Tsfkpk2Db79V6t+9OwwYAL16QYEC6XpGzauNQUGXUsYKIcYBewBTYImU8rIQ4t0n1xdJKT2EELuBC0A88IeU8tLLnLhGkx4ehD1g8v7JtKvUjv52/Tl8+DCXLl3it99+S3dxkX2T9uF3wY8hO4dQoExSwb1wQWW2TZf15sIFmDkThgyBX3+FggXTNR+N5ilGBRZJKXcCO587t+i5z7OB2Zk3NY0m89jovpGHkQ+Z320+QghCQ0MBEvLeG8uVLVc4vfA0zSc0p2r3qsmur1+vzN1Gr85DQmDYMChSBH7+WYu55oXIdpGiGs3L4FrgNfKb56d2idoA7Nq1izx58lC1anJRTo0Q3xC2vr6V0g1L03FGym4qLi7QvHkqKW+fJyYG+vcHDw/l15iOgtMaTUro5FyaVwKfEB8qFKqAEIKwsDBWrFjBwIEDsba2NtwZiI+Lx3mYM7FRsfRb3S/F/OaxsaoOrtE1J2bPVoWdFy2Czp2NfxiNJhX0Cl3zShAaHUqBPMre7erqyuPHjxky5Pm9/dQ5NvMYtw7fwnGpI0WrpZwgLTBQBRFVqWLEgIGBMGsWODjoEnKaTEOv0DWvBMGRwRTKozIrHjp0CFNTU1q1amVUX58TPhz64hC1B9em3oh6qbZ7WkXQ3NzAgLt2QbNmEBamcrJoNJmEFnRNrkdKic8jH8oWLIuUEmdnZ+zt7bGysjLYN/JRJM5DnClUvhA9f+2ZaibQa9eUl6EQYGeXymBXr0K/fsqn0cxMmVtq1UqlsUaTfrSga3I9Xx3+inuh92hVvhXnz5/n6tWrDB482GA/KSU7xuzgkc8j+v7dF8tCKWea3L4dGjeGe/dUuH/rxHFGHh7w9dcqT26NGmp1PmOGclZPRyCTRmMM2oauydU4ezgz/fB0RtYbyRsN3uD1118nX7589O/f32DfCysvJIT2l2+Rckrc8+eVo4qdnYoMrVABZR+fPx82bAB3d7Vsb9VK5TPv3x/KlMnch9RonqAFXZOr+fP8n9ha2/Jbr9/w9fXl77//5q233jLo3RJ0PYid7+2kgn2FZKH9TwkMVBaUokXVwrvE06TREybAypXQpg289x706aNFXPOfoAVdk+uIio3i/P3zHPc5juttV/rW7EseszyMHz8eU1NTPv744zT7x8XE4TzEGRMzE/qu7JsktP8pgYHQqZNKuXL4cCIxB1Uk1MFBLdk1mv8QLeiaHI9/mD/HfY6rw/c4p++cJiouCoDK1pUZWmcoR44cwdnZmRkzZlCpUqU0xzs0/RB3Tt3BaZ0ThSoUSnbd11dtgF65Alu2KIeVBO7eBU9PlTdXo/mP0YKuybHs9trNB7s+wDPIEwBzE3MalWnE2CZjaVWhFS3KtaB0gdIAfLf2OwDGjRuX5pg3D9/k2HfHqP96fWr1T+qBEhkJP/yg9jSlhK1boUuX5wZYu1Zd7Ncvcx5So0kHWtA1OZb1l9dzL/QeMzvOpFWFVjQq3Yi85nlTbOvl5UWpUqUokEb2wojgCDYN20SRKkXo/lPS6kA7d8K4ceDtrbR6zhxIcaG/fTvUrQvVq7/Ak2k0GUMLuibHcj/sPtWKVuOT1p8YbBsQEECJEqmXupVSsv3t7YTeD+WNE29gYWXx5LzyOvziC+UyfuBAGt6G8fGqwlA6IlA1msxEC7omx3In5A7lC6XsTvg8YWFh5M+fP9Xrbn+54b7BnY4zO1KmsfJIiYhQpvD162HECPjtN1VzIlX8/VX2xFQjizSal4sOLNLkSKSU3H50m/IFjRN0f39/iqeSAjHwWiC7PthFpfaVaDVJpQOIjlal4zZsUDm0li41IOagvgHAyBJFGk3mowVdkyPxDfElODKYmsVqGmwbEhLCzZs3KZOCL3hcdBwbh2zE1MKUPsv7IEwEcXGqBujOnSoR4scfq9igNLl4UYX0pxn7r9G8XLSga3IcUkqWui0FoFm5Zmm2jYuLY9iwYYSFhaWYXfHg5we5d/YeDn84ULCcKi6xeTM4O6uVuVHFjFatgqZNITgY9u17LvZfo/nv0DZ0TY4iOi6a93e+z+Jzi3Go7kCj0o3SbD9v3jy2bdvGggULsLe3T3LN28Ub1+9dafhWQ2r2fbbSX7hQhfB/9JEREzp4UBnaW7dWLoslS2bksTSaTEGv0DU5hqsBV2nzVxsWn1vMlNZTcB7gjKmJaZp9tm3bRqNGjXjvvfeSnA8PDGfTiE0UrVaUrj92TTj/+LHS6FGjVCm5NPH0VLlZqlVTTulazDVZjF6ha7I9cfFxzDs5j/8d/B95zfKyzmkd/WulnlzLx8eHo0ePcvToUU6ePMnrr7+eJO2tlJJtb20j7EEYg7cOxiK/RcI1Ly/1WreugUm5u0PHjspmvmWLrgWqyRZoQddkayJiIui3rh+7vHbhUN2BRT0XJUR/AsTHx+Ph4cGxY8cSRPz27dsAFChQgHbt2vHOO+8kGfPc7+e4sukKnWd3pnTD0kmuPU2/kqqgS6l2S0eNUpUsDh+GdNQl1WheJlrQNdmW8Jhweq/pzf4b+/mlxy+82/jdJCvtjRs38vbbbxMUFARAyZIlsbe3Z+LEidjb21O3bl1Mn7Ob+Hv4s3v8bip3rkyLCS2SXPPzg7lzYcCAVDT64kWVSXH/fhUJun27kfXmNJr/Bi3ommzLJ/s+Yf+N/fzl+Bcj6ydPdrVjxw7i4uJYsmQJ9vb22NraplpRCCA2KpaNgzdikd+C3st6I0yStt2xQ1WFmzr1uY5RUTBpktotLVQIfvoJxowxotacRvPfogVdk205cvsI3ap0S1HMAfz8/LC1tWX06NFGjbf/0/34/evH4G2DKVA6eU4XX1/1WqNGopP37qnkLSdOwNix8NVXUKRIeh9Fo/lP0IKuyXZIKVl8djGXH1ymV9Veqba7e/cuZcuWNWpMz12e/DPvH5q+35Rqvaolux4aCsuXQ82aYPF0j/TiRejaFR49gnXrlEeLRpON0YKuyVbce3yPN7a+wS6vXXSu3JmPWqTuDO7r60uzZmkHFgGE3g9l88jNlKhTgs7fd06xzaRJcOOG2uME4NIllYXLwkIVrKhTJyOPo9H8p2hB12QbPPw96LC8A48iH7Gg+wLGNBmDiUg5VOLBgwcEBARQ1YCHiYyXbB65mejQaPqt7oeZZfIf+ZgYWLIE3nwT7O2BBw+US6K5uXJKr5Z8Ra/RZEe0oGuyBRf9LtJxeUdMhAn/vPkPdUqmvSI+cOAAAE2aNEmz3Ym5J7i+9zq9futFiVopp889d04l42rT5smJRYuUqLu5aTHX5Ch0pKgmy9ngvoHWf7XG3NScw6MOGxTz1atXM2rUKKpXr56moN85dYcDUw5Qs19NGr7VMNn1+HjluNKxIxQoAG3bPjm5aBF06wb16r3oo2k0/yla0DVZxtO8LP3X96dGsRocf/041YulXuknNjaWL774giFDhtC8eXNcXV3JmzflCkWRDyPZMHADBcoW4LXfX0vmzhgVpYo8jxsHrVopk3n58sCdO8qzxcEhMx9Vo/lPMErQhRDdhBBXhRBeQohP02jXRAgRJ4RwyrwpanIry/9dzoLTCxjfbDxHRx+lYuGKqbb19vamXbt2fPXVV4wcOZK9e/dStGjRFNtKKdn65lZCfENwWuNEXuvkov/PP8o8PmMG7N6tknEByrMFtKlFkyMxaEMXQpgCC4HOgC9wWgixVUrpnkK7WcCelzFRTe7jTsgdAOZ0mZNmki1nZ+cEX/OVK1cydOjQNMc9s+gMHhs96Dy7M+Wal0uxTWioeu3QIVGucymVwpcqBc2bp+9hNJpsgDGbok0BLynlDQAhxBrAEXB/rt37wEYg7V0qjeYJtx7dolCeQmmK+ebNmxkwYACNGzdmzZo1VEqxMvMz7rvdZ89He6jSvUqy0P7EHDwIZmZgY5PopIsLuLoqG3oa5eo0muyKMYJeFvBJ9NkXSOL8K4QoC/QBOqAFXWMEUbFRbLqyideqv5ZqGxcXFwYOHEiTJk3Yu3cvBQokj+5MMubjKDYM3EC+ovlSDO1/SmSkclPs0weS1I2+cEG9OmmLoSZnYowNPaXfCvnc53nAJ1LKuDQHEuJtIcQZIcQZf39/I6eoyY2cunOKh5EPcaqZunh+//33lCpVip07dxoUcyklO8bsIMgriL5/9yV/8dRX2MePQ1CQKvychPPnIV8+HdqvybEYs0L3BRJX4i0H3H2uTWNgzRNPgmJADyFErJRyc+JGUsrFwGKAxo0bP/+loHlFkFLiEeABQKXClVJtc/HiRTp06IC1tbXBMd2WunFx1UXafdWOSm1THvMpBw+CiUkiv3OAvXthxQr48EMjCohqNNkTYwT9NFBVCGED3AEGAUmKM0opEyyRQoilwPbnxVzzavMw8iH7b+xnt9du9lzfg2+ILxamFpQrmHzTMiAggFGjRnH37t1kZeNSwt/dn51jd2LTwQb7z9JuLyVs2KAqxiWpSTFmjErk8t136X00jSbbYFDQpZSxQohxKO8VU2CJlPKyEOLdJ9cXveQ5anIoodGhzP9nPjs9d3LS9yRxMo5CeQrRqXInvmj7Bd2rdKdovqSuh0ePHmXw4MH4+/vz888/89Zbb6V5j5jwGNYPWI+FlQV9VvbBxDRtK+L583DlCowfn+jkgwcqkcvcuZCKX7tGkxMwKvRfSrkT2PncuRSFXEo56sWnpcnpxMXHMWjDIHZ47qBxmcZMaT2FblW60axcM8xMUv6xc3FxoUePHlSoUIGTJ0/SoEEDg/fZPX43/pf9GbZnWIopcZ9n3jzlwDJgwJMTwcGqaAVA/frGPZxGk03RuVw0L4XPD37ODs8dLOyxkPeavGew/YkTJ3BwcKBq1aocOnQo1aChxFxcfZFzv5+j9ZTW2HaxNdj++nVYvVqlNbe2BrZtg3feUSv0qVOhXTsjnkyjyb7o0H/NS2Hh6YX0q9mPMY3HGGy7fv16OnXqRJkyZdi3b59RYh5wJYDtb2+nfKvytP+qvcH2UirtzptXpcrl6lUV3l+smAob/eYbvRmqyfFoQddkOlJKQqNDqV60epol4eLj45k6dSoDBgygXr16HD58mFKlShkcPzosmnX91mFmaYbTGidMzAz/GO/cCQcOwMyZULYs8KSQNL/8Ao0aGftoGk22Rgu6JtP53vV74mRcmom2AFatWsWMGTN48803OXjwIKVLlzY4tpSS7e9sx9/Dn36r+1GwXEGDfQDOnFEL8Ndff3Liab25QoWM6q/R5AS0DV2TaUTHRfPHuT+YcmAKg2oPYnjd4Wm2v3nzJgC//vorZmbG/SieWXSGi6su0v7r9lTuVNnouXl5qWyKlpZPTvz5J1StCrVqGT2GRpPd0YKueWGCIoL47cxv/HzqZ+6F3qNNxTYscViSprkFICgoiDx58hgt5r7/+LL7w91U6V7FoL/58zx8qMzlgHJRdHWF2bNVhJFGk0vQgq7JMOEx4Uxzmcais4sIjwmni20X/nL8iy62XQyKOcDu3btp0SL1BFqJCXsQxnqn9RQsW5C+K/ummqclNUJCwMrqyYeDB9Vrr9QLUGs0OREt6JoMcdHvIgM3DORKwBWG1xvOxBYTqVuyrtH9r169ypUrVxg7dqzBtvGx8WwYtIHwgHBeP/46eYukP/jn7l1o3PjJhz17VIrc6mnb+DWanIYWdE26CIoI4tfTv/LN0W8olKcQe4fvpVPlTuke58aNGwA0MsLD5MBnB7h58CaOSx0p3cDwxunz7NmjrCyDBwOPHyv/8zfe0G6KmlyHFnSNUVwNuMq8k/NY9u8yImIjcKjuwOJeiylpVTJd48TFxXH8+HGWLl0KQBEDmQ0vr7vM8dnHafRuI+qPrJ/uef/7r8qGW6cOfPwx8PffKn/ukCEG+2o0OQ0t6Jo0iZfxjNw8kpUXVmJhasGwOsMY33y8wULOiQkLC2Pv3r1s3bqV7du3ExAQgLm5OY6OjtgkqTCRFL+LfmwZvYXyLcvT/afuGZr/jBmQJw/s2AEFreLhxx+hYUMw0nav0eQktKBr0sQryIuVF1Yyuv5ovuv4XbpW5Ldv3+b9999n7969REZGUrhwYXr27ImDgwPdunWjYMHUfcjDA8NZ23steQrlof+G/phapF7VKCUCA+HUKThyBDp2fBJMtGSpihBdvVqbWzS5Ei3omjTZd30fAOOajkuXmEdHR+Pk5MSVK1d4++23cXR0xN7eHnNzc4N942Li2DBgAyF3Qhh1eJTBpFuxsaq288mTz45r19Q1U1Po3Rvw91cx/23awMCBRj+HRpOT0IKuSREpJXOOz+GT/Z/QrGwz6pQw3sQC8Nlnn3H69Gk2btxI375909V378S9eLt403tZb8o1S7nIMyhX8s8+U1Gg4eHqXIkSypoyerR6bdToibvioo2qTNH8+Xp1rsm1aEHXpMjnBz/nm6Pf0N+uP0t7L8Xc1PDKOjG///47AwcOTLeYn/vzHKd+PkWLiS2oN6Jemm2XLoXTp+Gtt6B5cyXgFSumotdXr6rycnWNd63UaHIaWtA1yZBS8uf5P+lZtSdrnNZgItIfTRkeHk7lysaH5gPcPHyTHWN2YNvVlk6zDLtC3rsHNWrATz8ZMbivL1SooFfnmlyNjnvWJMMjwIN7ofd4rdprGRJzFxcXYmNjKV68uNF9gq4Hsa7fOorYFlEZFA1UHgLw8Xmy2WkMoaHP1ZzTaHIfWtA1yfj26LfkNcuLQ3WHdPe9f/8+Q4YMoUaNGgbLxz0lIjiC1a+tBgmDtw3GsrClwT5Sws2bysRiFEFBOrOiJtejBV2ThPP3zvP3xb8Z33w8pQsYH5UZGRnJggULaNSoEY8ePWLdunVYJSRPSZ246DjWO60nyCuIARsHUKRK2oFGT7l+XeVnMcokLqUyuZQpY9TYGk1ORQu6JgEpJR/v+5iieYsyudVko/pERUXxyy+/UKVKFd5//31sbW05cOAAdeoY9oqRUrJ9zHa8Xbxx+MOBSu0qGXXPiAiYPl29b9nSiA4zZqhkLkY11mhyLnpTVJPAnut7cPF2YX63+RS2LGyw/YULF+jTpw83btygVatWLFu2jA4dOhiVaRHg2HfHcFviRptpbQx6tDzl+nXo3x/On4dp06B2bQMdtmyB//0Phg1T7jAaTS5GC7oGUKvlaQenUalwJd5p/I7B9ps3b2bYsGEUKlSIPXv20LlzZ6OFHODf5f/iMtWFOkPr0O7Ldgbbx8YqF/Jp01Qo//bt0LOnETdavBhsbOCPP7SHiybXowVdA8AB7wOcuXuGP177AwtTi1TbPXjwgJkzZ/Ljjz/SpEkTNm/eTJl02qavbb/Glte3YNPBBoc/HQx+EVy4oAKFzp1TIv7LL8oDMVViYlRWruPHVfRRnz7qW0CjyeVoQdcAsNF9I/nN8zO07tAUr/v7+zN79mwWLlxIZGQkb775JvPnzydv3vTlJr997Dbr+6+nVP1SDNw8ELM8af8IHjwIDg6QPz+sXavMLcn0PzgYTpx4JuCnTj0LHa1QQWdW1LwyaEF/hZFSctL3JH+5/cXKCyvpXrU7lmbJXQZnzJjBt99+S2RkJIMHD+Z///sfNWrUSPf9/C748XevvylUoRBDdw0lT4G0V8179yoxt7VV75P5nO/dCx99BO7u6rOpKTRoAG++Ca1aqU3QcqmnDtBochta0F9B7j2+x4oLK/jL7S+uBFwhn3k+BtYeyJftvkzW9s6dO0ydOpVu3brx448/ZkjIAYJvBLOy60osrCwYtncY+YvnN9hn1izlaXjkCBQtmkKDdevg9m345hsl4E2aqKW8RvOKogX9FcI72Jvvjn3HUrelxMTH0Kp8K/50+JP+dv0pkCfljIbuT1a/kydPzrCYh/qFsqLLCmKjYhl9dDSFKxY2PFdvOHtWJUZMUcwB7t+HatVg6tQMzUujyW1oQX8F8AryYsbRGSz/dzmmJqa81fAtPmj2AdWLpV1TMyQkhK+++gozMzOj/MpTIvJRJKu6rSL0XigjDoygRK0SabaPi1O5WaZNAxMTGD48jca+vtqkotEkQgt6LueEzwk6Lu+IRDKu6TgmtZxE2YKGE6AEBQXRvXt3zp07x99//02xYsXSfe+YiBjWOK7hwaUHDN42mHLN0xbfwECVu/zYMeXN8uuvUL58Ko2lVOYWXXlIo0lAC3ou5lrgNV5b/RplC5bl0MhDRgn5U9555x3c3NzYuHEjDg7pz+kSGxXLeqf13Dpyi76r+lKlW5U029++DV27KlPL8uUqDihNb8Y7d5R3S82a6Z6bRpNbMSr0XwjRTQhxVQjhJYT4NIXrQ4UQF54cx4UQxoX9aV4aYdFhOK5xRAjBrqG70iXmUkoOHDjA8OHDMyTmcTFxbBy0Ec+dnvT6rRd1Bqdtrtm6FRo3Vulw9+xRZhaDMUArVqjXZs3SPT+NJrdiUNCFEKbAQqA7YAcMFkLYPdfMG2grpawLfA0szuyJaownPCac93a+x9WAq6x1WkuVImmvjp/n2rVrBAcH0ywDYvlUzK9svkL3n7vT6K1Gqc8zHN54AxwdlUvi8ePQtq2BG8TEwLhxqlRRjx7qm0Cj0QDGmVyaAl5SyhsAQog1gCPg/rSBlPJ4ovYnAb1T9R8SFBGE621Xjt4+ytHbRzl79ywx8TFMtZ9KB5sO6R5v8eLFmJqa0rVr13T1i4uOY+PgjXg4e9Dtp240Hdc0zfbLlsGSJfDpp/Dll2CReoCqwssLRo5Uyj9pEnz3nfI912g0gHGCXhbwSfTZF0hr6fYGsOtFJqVJGyklGz024uLtwtHbR7n04BIA5ibmNCnbhAktJtCuUju62qZPkAEePnzI4sWLGTRoEBXSjK9PSkxEDOud1uO505NuP3Wj2QeGV/f376vXb79VHi2pEh8PCxYo5c+TB1avhkGDjJ6bRvOqYIygp2TNlCk2FKI9StBbp3L9beBtIF1ioUnK8n+XM2rLKKwsrGhZviUDaw3EvoI9Tcs2Ja95+kLxn2fhwoWEhoYyadIko/tEh0azxnEN3ge96bW4V5pmlsScPau8WNIU88hIGDAAtm1TJpbFi9NRpkijebUwRtB9gcTOY+WAu883EkLUBf4AukspA1MaSEq5mCf29caNG6f4paAxzEaPjVQsVBGvD7wwM8k8R6Xw8HDmzZtH9+7dqVfPuH3tyEeR/N3zb3xP+NJneR/qDjOuCPPdu7BrF0yZktbgkSqx1u7dyjn9/fd1xkSNJg2M8XI5DVQVQtgIISyAQcDWxA2EEBUAZ2C4lPJa5k9TA+Af5s8w52Fsu7aN/nb9M1XMAZYtW0ZAQACffprMkSlFwvzDWNFpBXf+uYPTWiejxRxUBsX4eOjWLY1Gn3yi3F5+/x0++ECLuUZjAIOKIKWMFUKMA/YApsASKeVlIcS7T64vAj4HigK/PEmFGiul1O4HGUBKSXhMOMGRwQRFBBEcEUxwZDDewd58e/RbQqJC+KLtF0xpndbSNmP3XbBgAQ0bNsTe3t5g+4c3H7Ky60oe3X7EwE0DqdarmsE+YWFqP/PgQZXPHKBUqTQ6HDqknNPffNO4h9BoXnGMWuJJKXcCO587tyjR+zcB/VuXiMjYSIIjnohyZHCCMCcW6dSuxcTHpDhmi3It+P2136lVolamz/fo0aO4u7vz119/GcxP7nfRj5VdVxIbEcvw/cOp0Crl/RAp4fBhOHBAifipU8rr0MxM5dH65huVSTHFjkuWwOXLRlax0Gg0oCNF001geCDbr23n1qNbSrAjnwl0YgGPjI1MdQyBoJBlIawtrbHOa421pTXlCpZL8rlI3iIJ763zqs8VC1VMV1Wg9LB7927MzMxwcnJKs53XHi/W919PngJ5GH10NCVqp5ybRUr48EP4+WflWdioEUyYAO3bq8SIqdaPvn9flYrbvl01Hj/+xR5Mo3mF0IJuBIHhgWy+spl17us4cOMAcTIOgAIWBRLE1trSmurFqj8T4zTEuVCeQpiaZB//aT8/P7Zt20aTJk2wSlVp4cyiM+wct5MStUswZPsQCpYrmGrbmTOVmH/wAXz9NRRMvekzrl6Fjh0hIAB+/FF1TtMFRqPRJEYL+hNi4mLwD/fnQdgD/EL98Avz437ofQ7ePMj+G/uJjY+lsnVlJrWcRP9a/albsm6mb0r+13h6ejJnzhyWLVtGdHQ0v//+e4rt4uPi2TdpHyd/PEnVnlXpt7qfweIUs2Ypa8mPPxqpye7uSszj4uCff8BILxuNRvOMnK1IBgiPCU8i0Mneh/nhF6reB0ak6GmJTWEbJraYyIBaA2hQqsFLM3n8l5w5c4aZM2fi7OyMhYUFI0eOZOLEiVSrlnxjMzo0GudhzlzdcpWm7zel69yumJilrdDh4fD4sSoeZJSYb9kCo0ZB3rzK6K4Tbmk0GSLHC/rjqMesurgKt/tuyUQ7NDo0xT6F8hSiRP4SlLQqiV1xO9pVakfJ/CUpaVVSnX/yvmT+klhZWOUKEQdlWvnkk09YtmwZhQsXZsqUKbz//vuUSsXVJMgriDW91xDgEWB09Gd0NDg5KRu6UXlZPvlELeMbNVIViCpXzsCTaTQayMGCfiP4Bj//8zNL3JYQEhVCsXzFEoS4admmlMxfMkG0E78vkb9EinUzczOxsbEsWrSI//3vf4SHhzNlyhSmTJlCgQIpVykC8NzlifMQZ4SJYNieYVTulLbQPnoEf/8NixYpH/PffoNOnQxMbPlyJebjxsGcOSqsX6PRZJgcJ+hu99344tAXbLu6DVMTU/rb9efDZh/SrJxOo5oSoaGhdOzYkVOnTtGpUycWLFhA9eqpVyqKj4vn8FeHOfL1EUrWLcnATQOxtrFOtf3p06oQxdq1ytRSr54S9sGDjZjcnTvqdd48nWRLo8kEcpygB0UEcdznOJ/Zf8aYxmPSlef7VWTChAmcPn2aFStWMHTo0DTNR6H3Q3Ee6oy3izf1R9Wnx8IemOczT7Htw4cwebIK4rSygqFD4e23leXEaAvV7dvK/UWLuUaTKeQ4QW9fqT0+H/m8cmaTjLBlyxZ+//13Jk+ezLBhw9Js67Xbi80jNxP1OAqHPx1o8HqDFNs9fqxW459/Dn5+KovttGmQhvUmZaKjwdkZevVKZ0eNRpMaOU7QhRBazI0gNDSUMWPGUL9+fb7++utU20WHRrP3472c/e0sJWqXYOTBkRS3K56kjZTg6qqCN9etUyH8DRuqBIiNjEusmJxTp1QJuf79MziARqN5nhwn6BrjmDFjBvfu3UtwTUyJW0dusXnUZh7efEiLj1vQ4esOmFkm/ZGIiYE2beDkSbUKHzIEXn9dVX7LsPOPlCqICKBSpQwOotFonkcLei7i8ePH7Nmzhy1btrB27VqGDx9O8+bNk7WLiYjB5X8unPzxJNY21ow+MpoKrVPOx/Lvv0rMp05VqW7z58/g5EJCwMVFpcLdvRtu3QJzcyini1tpNJmFFvQczt27d9m6dStbtmzBxcWF6OhoihQpwtChQ5k9e3ay9jcP32THuzsIuBJA4/ca03lWZyysUlnB31KpVISAd97JgJiHhcHChbBzp7LZxMaqHdSOHZX/eY8eUKxY+h9ao9GkiBb0HEhsbCw//PADGzdu5PTp0wDY2toybtw4HB0dadmyJWZmSf9rwx6EsW/SPv5d/i+FKxVm2J5h2HZJKdWhYudOZV6Jj4cVK1RloXQRH69cX7Zsgfr1YeJElfy8ZUsjiodqNJqMoAU9B7J06VI+/fRTmjRpwrfffoujoyN2dnYpuiTGx8Vz/s/z7P90P9Gh0bT+rDVtprZJ1R3xKePHq0pv27eDjU0GJjl9uhLzn35SSbY0Gs1LRwt6DuTo0aOUKFGCf/75J02/8puHb7Jn/B7uu92nYtuK9Py1J8VrFk+1/VMuXABPT1WXOUNiDjB/viof9/77GRxAo9GkFy3oOYyAgAAOHDhAy5YtUxXz+273OfTFIa5uvUrB8gXpt6YftQbUMionTXw8jB0L1taqNnOGkBJCQ1WSrVySB0ejyQloQc9BeHl50b17dwICAhg7dmyy63fP3OXI10e4uvUqeQrmocOMDjQf3xzzvGmbVxLzyy9w7Bj89RcUN7yYT5m5c1Ua3KpVMziARqPJCFrQszlSSkJCQjh16hRDhgxBSsmBAwdo1apVQhufEz4c+foIXru8sLS2pN1X7Wj2fjMsC6cvAOv6deV80r07jByZgcnGxKhvgkmTVMrFESMyMIhGo8koWtCziJiYGB48eMC9e/e4f/8+9+/fT3j//LnISFXOztbWll27dlH1ycr31tFbHPnqCDf23yBv0bx0/K4jTd5rQp6C6c9aKCWMGaPSqixenE5LycOHKqnL/Png66s8WZYt09WGNJr/GC3omcjT1bQhgb5//z4BAQFIKZONUaRIEUqXLk2pUqVo1aoVpUqVolSpUpQuXZru3btTuFBhPHd6cnz2cW4eukn+EvnpPLszjd9tnKo/uTHs2AH79qmycUbH+kREKG+WX35RNvP27VX+3O7dtZhrNFmAFvR0EhcXx9atW7l8+XKKoh0REZGsj4WFRYIwV65cmZYtWyaI9tPXUqVKUbJkSfKkkhM8PDCc83+c58yvZ3jo/ZACZQrQ9ceuNHq7kUEXRENcv67SkZcqBe++a2Qnd3cYOBAuXVIO6xMnqgQvGo0my9CCbiRxcXGsWbOGr776imvXrgFgbW2dIMgtWrRIJtJPX62trdNd9SgmIgaf4z54H/DG+4A3d8/cRcZLKrapSMfvOlKzT01MLTKWdjY2Fk6cUMm1tm2DK1fU+VmzwMzQT8TDh8omM326ivrctUsFDGk0mixHC7oB4uLiWLduHV999RVXrlyhTp06bNy4kZ49e6a6ms4I8bHx3D1zF28XJeC3XW8TFxWHMBWUa1YO+6n21BpQixK1S2T4Hv/+q/YsV62CgACVSqVtW2U779XLQPW369dVkNCSJSqkv0cP+OMPKF06w/PRaDSZixb054iLi8PHx4dr167h4eHBb7/9hoeHB7Vr12b9+vX07dsXk0ywD0sp8b/sz40DN/A+4M2tw7eICokCoGS9kjR5rwk2HW2o2KYieQpk/ItDSvjzT2XmPn9eRd07OCgf865dVX2JNImPhzfeUJucZmaqFNFHH6lwfo1Gk614JQVdSsmDBw+4du0a165dw9PTM+G9l5cXUVFRCW1r1arFunXr6Nev3wsL+cObDxME3NvFmzC/MACsba2pNagWlTtWplL7SuQvntGUhkmJioI334SVK5X+zp+vzN1Fi6ZjkOvXYelS5cf43Xd6Ra7RZGNytaA/evQoiVgnFu+QkJCEdubm5lSpUoVq1arRo0cPqlWrRtWqValWrRqlSpVKt/37KWEPwvB28ebGgRvcdLlJ8I1gAKxKWVG5Y2VsOtpg09GGwhULZ8bjJiE8HF57TWWs/fZblfo2Q4+xf796ff99LeYaTTYnVwj6gwcPcHV1TSbcfn5+CW2EEFSsWJFq1aoxYsQIqlWrliDcFStWxDQT6lpGhURx68ithFX4g4sPAMhTKA+V2lWi2fhmVO5YmWI1i2X4SyItQkPBw0M5nixZAsePK0tJhuJ7pFQ284kToXFjqFs30+er0Wgylxwv6Ddu3KBFixY8eKDEs2TJklSrVo1evXolEW1bW1ssLTO3dF1sZCw+J3wSTCh3Tt1BxknMLM2o0LoCtQfXpnLHypRuWBoTs8z1y750SdnEL19Wx6VLcPPms+v588Py5SqDbYb46ivlydK3rxrI/MVcIzUazcsnRwt6YGAg3bt3JyYmBhcXFxo2bEihQoUy/T5SSiICIwi6HkTwjWCCPIO4ffQ2t4/dJjYyFmEqKNukLK0/bY1NRxvKtyifrJRbZrJ6tbKFg9LZGjWgeXO1d1m7NtSqpTxWMvxHh5TKNbF7d1i/XgcJaTQ5hBwr6BERETg4OHDr1i32799P69at09VfSkn042jCA8IJDwgnzD8s4X24f3jC+0e3HxF8PTjBA+UpJWqXoNE7jbDpaEOltpUyFG6fXiIjYcMGVQ6uVi2ltVWqvITF85UrcPcuTJumxVyjyUEYJehCiG7AT4Ap8IeUcuZz18WT6z2AcGCUlPJcJs81gfj4eEaMGMGJEydYu3YtrVu3Ji46jvDApGKcINKJziV+Hxcdl+L4JuYm5CuWj/zF81OgTAHKtyyPta01RWyLYG1rjbWNdYaiM6VUohwenv4jKEjViwgMVEkMf/lFZad9KcyYAZaWyr9Ro9HkGAwKuhDCFFgIdAZ8gdNCiK1SSvdEzboDVZ8czYBfn7xmOu6brvHuW59xNHAjjoUduPmeNzNGfUtMeGyqfSwL5yF/0bzkK2JJ4TJ5KV27iPpc1JJ8RSzJl/A+L/mL5cXCyjzZpqWU8MBfcPmaGR47oggKjiY8QiQ6UK/hST+HhT9rExEBUqZ/MzRfPkm+vJJ2LWIYMzqSDvYxymMlIN1DGWbvXhV5NHkylCnzEm6g0WheFsas0JsCXlLKGwBCiDWAI5BY0B2B5VJlmzophCgshCgtpbyX2ROe8/1+jgZuoADuiIdXiOYM9bhAQR6Tj3DyEfbkVR0mxBPz0JzIh5ZEXrckijxE8uzVn+TnEr/3oTyXqYU7dgSS1IHbnOgk90p8lEr0Pv9zc0rPYUkkIhz1d8/2J8fLpkUL5eeo0WhyFMYIelnAJ9FnX5KvvlNqUxZIIuhCiLeBtwEqVKiQ3rkC8OFnbbkzeRmRsR3ZfsORzfF9yWcRQ9lCYUTGmhIVa0pkjCmRsWZEx764K2LhvJHUKh1Ev9L3sSvljl2pYOxKBVGiQATmpvHpGCnvkyM9UT1ZQKFCMGiQ9mrRaHIgxgh6SjaC5/O+GtMGKeViYDFA48aNk+eONYJ6r9Viz2u1AHj8GA4dgn37zPH3L4ylpTL95smT9DWlc2lde/peHZYIUQbQ5geNRpO9MUbQfYHyiT6XA+5moE2mU6CAioZ87bWXfSeNRqPJ/hjjk3YaqCqEsBFCWACDgK3PtdkKjBCK5sCjl2E/12g0Gk3qGFyhSyljhRDjgD0ot8UlUsrLQoh3n1xfBOxEuSx6obbvRr+8KWs0Go0mJYzyQ5dS7kSJduJzixK9l0DyMvQajUaj+c/QYYAajUaTS9CCrtFoNLkELegajUaTS9CCrtFoNLkELegajUaTSxDKQSULbiyEP3ArS27+cinGy0mblR3Izc8G+vlyOrn5+RI/W0UpZfGUGmWZoOdWhBBnpJSNs3oeL4Pc/Gygny+nk5ufz9hn0yYXjUajySVoQddoNJpcghb0zGdxVk/gJZKbnw308+V0cvPzGfVs2oau0Wg0uQS9QtdoNJpcghZ0jUajySVoQc9khBCzhRBXhBAXhBCbhBCFs3pOmYkQor8Q4rIQIl4IkWtcxIQQ3YQQV4UQXkKIT7N6PpmJEGKJEOKBEOJSVs8lsxFClBdCHBRCeDz5ufwwq+eUmQghLIUQp4QQ/z55vi/Taq8FPfPZB9SWUtYFrgG5rdryJaAvcCSrJ5JZCCFMgYVAd8AOGCyEsMvaWWUqS4FuWT2Jl0QsMFFKWRNoDozNZf93UUAHKWU9oD7Q7UkRoRTRgp7JSCn3Siljn3w8iSrHl2uQUnpIKa9m9TwymaaAl5TyhpQyGlgDOGbxnDINKeURICir5/EykFLek1Kee/L+MeCBKlCfK5CK0CcfzZ8cqXqyaEF/ubwO7MrqSWgMUhbwSfTZl1wkCq8KQohKQAPgnyyeSqYihDAVQrgBD4B9UspUn8+oikWapAgh9gOlUrg0VUq55Umbqag/B1f9l3PLDIx5vlyGSOGc9ufNQQghrICNwHgpZUhWzyczkVLGAfWf7MdtEkLUllKmuB+iBT0DSCk7pXVdCDES6AV0lDnQ0d/Q8+VCfIHyiT6XA+5m0Vw06UQIYY4S81VSSuesns/LQkr5UAhxCLUfkqKga5NLJiOE6AZ8AjhIKcOzej4aozgNVBVC2AghLIBBwNYsnpPGCIQQAvgT8JBSzs3q+WQ2QojiTz3lhBB5gU7AldTaa0HPfBYABYB9Qgg3IcQiQx1yEkKIPkIIX6AFsEMIsSer5/SiPNnEHgfsQW2qrZNSXs7aWWUeQojVwAmguhDCVwjxRlbPKRNpBQwHOjz5fXMTQvTI6kllIqWBg0KIC6iFxz4p5fbUGuvQf41Go8kl6BW6RqPR5BK0oGs0Gk0uQQu6RqPR5BK0oGs0Gk0uQQu6RqPR5BK0oGs0Gk0uQQu6RqPR5BL+DxlHp/xYEyW5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# one iteration with Normal distribution\n",
    "\n",
    "n=100\n",
    "alpha=0.05\n",
    "r=norm.rvs(size=n)\n",
    "epsilon=np.sqrt((1/(2*n))*np.log(2/alpha))\n",
    "\n",
    "F=lambda x: sum(r<x)/n\n",
    "xx=sorted(r)\n",
    "F_n=np.array(list(map(F,xx)))\n",
    "L_n=np.maximum(F_n-epsilon,0)\n",
    "U_n=np.minimum(F_n+epsilon,1)\n",
    "CDF=np.array(list(map(norm.cdf,xx)))\n",
    "\n",
    "Sde=np.sqrt((F_n*(1-F_n))/n)\n",
    "NL=F_n+norm.ppf(alpha)*Sde\n",
    "NU=F_n-norm.ppf(alpha)*Sde\n",
    "\n",
    "df=pd.DataFrame({\n",
    "    'x': xx,\n",
    "    'F_n': F_n,\n",
    "    'U_n': U_n,\n",
    "    'L_n': L_n,\n",
    "     'CDF': CDF,\n",
    "    'normal_upper': NU,\n",
    "    'normal_lower': NL\n",
    "})\n",
    "df['in_bounds'] = (df['U_n'] >= df['CDF']) & (df['CDF'] >= df['L_n'])\n",
    "\n",
    "plt.plot('x','L_n',data=df, color='red')\n",
    "plt.plot('x','U_n',data=df,color='green')\n",
    "plt.plot('x','CDF',data=df, color='purple')\n",
    "plt.plot('x','normal_upper',data=df,color='black')\n",
    "plt.plot('x','normal_lower',data=df, color='blue')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "56649addc6224814b2785532d570bb9f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=1000.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Average fraction in DKW bounds: 0.958\n",
      "Average fraction in normal bounds: 0.000\n",
      "Average distance in DKW bounds: 25.3\n",
      "Average distance in normal bounds: 15.6\n"
     ]
    }
   ],
   "source": [
    "# 1000 iterations faster\n",
    "bounds=[]\n",
    "normal_bounds=[]\n",
    "dst=[]\n",
    "normal_dst=[]\n",
    "n=100\n",
    "alpha=0.05\n",
    "epsilon=np.sqrt((1/(2*n))*np.log(2/alpha))\n",
    "\n",
    "for k in tqdm(range(1000)):\n",
    "    r=norm.rvs(size=n)\n",
    "    \n",
    "    F=lambda x: sum(r<x)/n\n",
    "    xx=sorted(r)\n",
    "    F_n=np.array(list(map(F,xx)))\n",
    "    L_n=np.maximum(F_n-epsilon,0)\n",
    "    U_n=np.minimum(F_n+epsilon,1)\n",
    "    CDF=np.array(list(map(norm.cdf,xx)))\n",
    "\n",
    "    Sde=np.sqrt((F_n*(1-F_n))/n)\n",
    "    NL=F_n+norm.ppf(alpha/2)*Sde\n",
    "    NU=F_n-norm.ppf(alpha/2)*Sde\n",
    "\n",
    "    df=pd.DataFrame({\n",
    "        'x': xx,\n",
    "        'F_n': F_n,\n",
    "        'U_n': U_n,\n",
    "        'L_n': L_n,\n",
    "         'CDF': CDF,\n",
    "        'normal_upper': NU,\n",
    "        'normal_lower': NL\n",
    "    })\n",
    "    \n",
    "    all_in_bounds = ((df['U_n']>=df['CDF'])&(df['CDF']>=df['L_n'])).all()\n",
    "    bounds.append(all_in_bounds)\n",
    "    all_in_bounds = ((df['normal_upper']>=df['CDF'])&(df['CDF']>=df['normal_lower'])).all()\n",
    "    normal_bounds.append(all_in_bounds)\n",
    "    \n",
    "    d=sum(abs(CDF-L_n))+sum(abs(CDF-U_n))\n",
    "    dst.append(d)\n",
    "    d=sum(abs(CDF-NL))+sum(abs(CDF-NU))\n",
    "    normal_dst.append(d)\n",
    "\n",
    "print('Average fraction in DKW bounds: %.3f' % np.array(bounds).mean())\n",
    "print('Average fraction in normal bounds: %.3f' % np.array(normal_bounds).mean())\n",
    "\n",
    "print('Average distance in DKW bounds: %.1f' % np.array(dst).mean())\n",
    "print('Average distance in normal bounds: %.1f' % np.array(normal_dst).mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d44d1a0505504099821781df88f6893e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=1000.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Average fraction in bounds: 0.967\n"
     ]
    }
   ],
   "source": [
    "# 1000 iterations original\n",
    "bounds=[]\n",
    "n=100\n",
    "alpha=0.05\n",
    "epsilon=np.sqrt((1/(2*n))*np.log(2/alpha))\n",
    "\n",
    "for k in tqdm(range(1000)):\n",
    "    r=norm.rvs(size=n)\n",
    "    \n",
    "    F_n = lambda x : sum(r < x) / n\n",
    "    L_n = lambda x : max(F_n(x) - epsilon, 0)\n",
    "    U_n = lambda x : min(F_n(x) + epsilon, 1)\n",
    "\n",
    "    # xx = sorted(r)\n",
    "    xx = r # No need to sort without plotting\n",
    "    \n",
    "    df = pd.DataFrame({\n",
    "        'x': xx, \n",
    "        'F_n': np.array(list(map(F_n, xx))), \n",
    "        'U_n': np.array(list(map(U_n, xx))), \n",
    "        'L_n': np.array(list(map(L_n, xx))), \n",
    "        'CDF': np.array(list(map(norm.cdf, xx)))\n",
    "    })\n",
    "    \n",
    "    all_in_bounds = ((df['U_n']>=df['CDF'])&(df['CDF']>=df['L_n'])).all()\n",
    "    bounds.append(all_in_bounds)\n",
    "\n",
    "print('Average fraction in bounds: %.3f' % np.array(bounds).mean())\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1a4a4411340>"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlQElEQVR4nO3dd3Rc5Z3/8fdXo2arW5KbbCPbmCIbbGxhSiBhCWFNqCEQSHZZnJAQzgEC50eyAbwJpJ3fpizHydJCQv0lGxPYFCeBsKFvqHZwwTXIDcu2ZHVLlmS15/fHHcljWWUkz+hO+bzOmXNn7lyNPs4xnzx+5t7nmnMOERGJfyl+BxARkchQoYuIJAgVuohIglChi4gkCBW6iEiCSPXrFxcVFbnS0lK/fr2ISFz629/+VuucKx7oPd8KvbS0lNWrV/v160VE4pKZ7RrsPU25iIgkCBW6iEiCUKGLiCQI3+bQB9LZ2UllZSXt7e1+R4mKzMxMpk2bRlpamt9RRCQBxVShV1ZWkpOTQ2lpKWbmd5yIcs5RV1dHZWUlM2fO9DuOiCSgYadczOwxM9tvZhsGed/M7CdmVmFm681s4WjDtLe3U1hYmHBlDmBmFBYWJuy/PkTEf+HMoT8BLBni/YuAOcHHjcBDxxIoEcu8VyL/2UTEf8NOuTjnXjez0iEOuRx4ynnr8L5tZvlmNsU5ty9SIUVizf6D+/nfXf/L+ur3qfmwgOptU2hvzqS7K0BXRyrdXQG6uwLg9H/icrSLzs/hrqWLIv65kZhDLwF2h7yuDO47qtDN7Ea8UTwzZsyIwK8WGVsPvPsA96+6ny21W6ByMfz3f0HD7GF+qmdMskn86Ox8nbuWRv5zI1HoAw1BBrxrhnPuEeARgPLy8pi8s0Z2djYtLS1+x5AY9eN3fkxXTxf3LP4xDz/yZdKz0/m3f4czz4TiYsjIgPR0b5uaCt4sm84Olv7Oi8qnRqLQK4HpIa+nAXsj8LkiMaehvYFPn3QVa3/6Fepr4M03obzc71QinkgU+krgFjNbAZwBNEVk/vz222Ht2mP+mCMsWADLl0fko5YuXUpubi6rV6+mqqqKH/zgB1x11VUR+WyJTRv3b6S2tZaaV6/h97+H++5TmUtsGbbQzexXeP8+KDKzSuAeIA3AOfcw8BzwSaACaAU+H62wsWbfvn389a9/ZcuWLVx22WUq9AT35LonSdlXzh8f/xiXXOKNOURiSThnuXx2mPcdcHPEEvWK0Eg6mq644gpSUlIoKyujurra7zgSRbubdnP/u/eT95f3GV9sPP547/y4SOzQtzXHICMjo++59/9rkqjufvluuuuPo2H7bO64A4qK/E4kcjQVukgYXtz+IvMblwFw2WU+hxEZhAq9n9bWVqZNm9b3uO+++/yOJD5r7WylqqWKpk1ncuKJMHu4085FfBJTi3PFgp6e8C4CeeKJJ454rXPXE1ddax042LulhM9c6XcakcFphC4yjPaudmgspaVxHIsX+51GZHAaoQ/je9/7Hs8888wR+66++mqWLVvmUyIZa509nbDHa3IVusQyFfowli1bpvJOcs452LOYtIxu5s0L+B1HZFCachEZRmtnK+xZzPFlTehmUxLLVOgiw2hsb4LqUzlxnm5OIrFNhS4yjF27O+FQHiedrGVwJbap0EWG8cHfvXnzspN1rb/ENhV6Pzt37mTevHlH7Lv33nv50Y9+5FMi8duOCm/i/NSyTJ+TiAxNhS4yjD27MiG1lZNn5fodRWRIMXva4u1/vp21VWsj+pkLJi9g+ZLlEf3M8847jzPOOINXXnmFxsZGHn30Uc4999yI/g7xV+2ePJiwnfTUecMfLOIjjdAjoKuri3fffZfly5fzrW99y+84EmEN+wqwgh1+xxAZVsyO0CM9kg6XDbLI9WD7Aa680lvgY9GiRezcuTMascQnzkHDvjxyF9f7HUVkWBqh91NYWEhDQ8MR++rr6ykaYgHs3nXRA4EAXV1dUc0nY6uhAbraxjOzVGe4SOxTofeTnZ3NlClTeOmllwCvzP/85z9zzjnn+JxM/LBjZzcAs0t1iajEvpidcvHTU089xc0338wdd9wBwD333MNsLYKdlN6vaAQKmaVClzigQh9AWVkZr7zySljHvvrqq33Pi4qKNIeeYFZvrgYKOf2kKX5HERmWplxEhrChogFSOvn4KTplUWKfRuhhuvnmm3njjTeO2Hfbbbfx+c9/3qdEMhY+3N1DIHc/E7JK/I4iMiwVepgeeOABvyOID+qqM8gubAJU6NJPZye0tIzuccUVcN11EY+kQhcZQmt9ASWzG/2OIcfCOWhrG335hj4OHjz8vKMj/AwZGZCdffhRVxeVP6oKXWQQzjm6DhRTUFztd5Tk1NMDBw5AU5P3aGw8/HywfYOVsHPh/97Q4u19FBbCcccdfp2VNfBxAz2yshirO6Oo0EUGUdPUAu35TJo0gjIQj3NemQ5VvsO9PnBg+N+TkQF5eYcfOTkwffrRhRpu+Y4bBynxe66ICn0AVVVV3H777axatYqMjAxKS0tZvnw58+fP56STTqK9vZ2cnBxuvvlmrr/+egCeeOIJvva1r1FS4s21nnrqqTz11FN+/jHkGL2/rQbIoWSq7iM6IOdg/XpYsQLeeuvoQu4Z5oYgqamHizg/39sef/yRr0MfA+3L1JLGoVTo/Tjn+NSnPsX111/PihUrAFi7di3V1dXMnj2bNWvWALB9+3auvPJKenp6+s50ueaaa7j//vt9yy6R9cbmCmAWC46f7HeU2LJli1fiTz/tPQ8E4PTTYcaMgYt3sHIeNw6GWCNJRk6F3s8rr7xCWloaN910U9++BQsWHHXB0KxZs7jvvvu44447dOpignqvYg8Ap58ww+ckEdDZCa2t3nxy/+1A+wY7prISNm70ivijH4XbboNPfxqKi/3+EwoxXOi33w5r10b2MxcsgOXLhz5mw4YNLFq0KKzPW7hwIVu2bOl7/fTTT/PXv/4V0DnqiaCq2ps7L5kS5f9MnPPOmBhtyYZzTGfnyDKZeXPP48cfuS0pgS99Ca6+GqZOjc7/HjJqMVvo8cD1++ZcUy6JpbbGmzufODFCH7h3L3zxi7B799HFO9x8c3+BgFewvY/Q4i0qOrqIe7cD7Rtom5mp6ZA4FFahm9kS4MdAAPi5c+7f+72fB/wCmBH8zB855x4/lmDDjaSjZe7cuTz77LNhHbtmzRpOPvnkKCcSP3R0d7BrTxvpWa1kZIyPzId+5zvw4otwySXDl+pwhZueHplMklCGLXQzCwAPAJ8AKoFVZrbSObcp5LCbgU3OuUvNrBjYama/dM6N4Mz72HD++edz991387Of/YwvfelLAKxatYrW1tYjjtu5cydf/epXufXWW/2IKVH20vaX6GouYGrxCEfOg/nwQ3j0UbjhBnjooch8pkg/4ZxwuRiocM5tDxb0CuDyfsc4IMe82/pkA/VAXN7pwcz47W9/y1/+8hdmz57N3Llzuffee5k6dSrbtm3jtNNO4+STT+Yzn/kMt956q+bJE9SDqx8krb2E46ZGaHT+b//mnd98992R+TyRAYQz5VIC7A55XQmc0e+Y+4GVwF4gB7jGOXfU0MbMbgRuBJgxI3bPHJg6dSq//vWvj9rf1tY26M8sXbqUpUuXRjGVjJX6tnr+9Pc/UdT9MyZNjMBFJi++CL/4BXz9695FLyJREs7f1oG+Gel/6dw/AmuBqcAC4H4zyz3qh5x7xDlX7pwrL9ZpThKjmg8143B0tOQyxJ0Hw7NypTdnPncu3HlnRPKJDCacQq8EQocV0/BG4qE+D/zGeSqAHcBJkYkoMrbautrAQUtj5rGdXv3yy3DllTB/Prz2mndBjUgUhVPoq4A5ZjbTzNKBa/GmV0J9CHwcwMwmAScC20cTqP+pgIkkkf9siWRHww5oz6O7K+XYRug//CFMnuxNuUyYELF8IoMZttCdc13ALcALwGbg1865jWZ2k5n1Xk75HeBsM3sfeAn4unOudqRhMjMzqaurS8jic85RV1dHptaeiHlb67ZCq9fkoyr0/fvhu9+FF16AL3zBWzBKZAyEdR66c+454Ll++x4Oeb4XuPBYw0ybNo3KykpqamqO9aNiUmZmJtOmTfM7hgzj/er3KXAn0MAIr2hfswZ+8hP41a/g0CFYsgS+8pVoxRQ5SkxdKZqWlsbMmTP9jiFJblPtJqal/gMNeMtgh+Wmm+CnP/Uu/LnhBrjlFtBFZzLGYqrQRWLB3ua9TO8pBcIs9E2bvDL/whfgP/7DW0lQxAfxu5K7SBT0uB6qWqrI6JgChPld5re/7S0F+/3vq8zFVyp0kRBVLVV0dHeQ0TkFszDONPzNb7x1wf/1X0f5DapI5KjQRULsbNwJQKC9mIKCYe5GduiQN3e+aBEsWzYm+USGojl0kRC7m7xVLroP5g0/3bJ5M9TUwH/+55jdBFhkKBqhi4SoaqkCoK15/NBfiNbWenPmAPPmRT+YSBhU6CIhunq8RUKbGgMDj9A7O+HHP4Y5c+CZZ+BrX4OysrENKTIIFbpIiIOdBwFobLCBC/2667z7Iy5e7N3x/gc/0J19JGZoDl0kxIb9G5hVMIv6+pSBC/2dd7wFt559VkUuMUcjdJEQW2q3UFZ4Co2NA5yD3t3tzZ2XlKjMJSap0EVC1LXVkedKgQEK/aWXoKUFPvrRMc8lEg4VukiIpvYm0g5NAgYo9Ice8q4EvfTSMc8lEg4VukgIh6PzYBYABQUhb7zyCvzud94XohkZfkQTGZYKXSSoq6eLQ12H6Dro3T3xiEL/5jfhuOO8S/xFYpQKXSRoW/02ul03BTYLCCl052DdOrjsMm8RLpEYpUIXCdrWsA2ArG7vJiR9hd7QAM3NoLX6Jcap0EWC6lrrAHBt3hKLfSvhNjdz5A6R2KRCFwna27wXgK7WbDIzoe/2r+3t3lZfhkqMU6GLBK3fv57pudNpbc448gvR7m5vm6oLqyW2qdBFgirqKzip6CQaG/vNrvQW+pCLo4v4T39DRYIa2xuZMG7C0YXeO4eek+NDKpHwqdBFAOcc+5r3UTy++OhCb2rytsPej07EXyp0EeDDpg9p7mhm7sS5Rxd6S4u3zc72IZlI+FToIsCmmk0AzJs4b/ApFxW6xDgVugiwuXYzACcVnkxjY7/ZlcZGb6vz0CXGqdBF8KZcstOzGUchXV391nHpnUPPzfUlm0i4VOgiwMGOg+Sk5/QNxo8YoXd1QSCg0xYl5ulvqAjQ0N5AXmZe32D8iNmVtjZIT/cjlsiIqNBFgN0HdjM9d/rAI/R162DuXD9iiYxIWIVuZkvMbKuZVZjZnYMcc56ZrTWzjWb2WmRjikRXXWsdxVnFA4/Q162DhQv9iCUyIsMuTmFmAeAB4BNAJbDKzFY65zaFHJMPPAgscc59aGYTo5RXJCp6XA8plnL0CL29HerqYMYMv6KJhC2cEfpioMI5t9051wGsAC7vd8zngN845z4EcM7tj2xMkejqdt0ELHD0GYq7d3vbkhIfUomMTDiFXgLsDnldGdwX6gSgwMxeNbO/mdm/DPRBZnajma02s9U1NTWjSywSBZ3dnaSlpB19lf9773nb+fN9ySUyEuEUug2wz/V7nQosAi4G/hH4hpmdcNQPOfeIc67cOVdeXFw84rAi0dLZ00lawCv0tLSQO82tWePt0JeiEgfCWeC5Epge8noasHeAY2qdcweBg2b2OjAf+HtEUopEWUd3BxmBDJqavNG59Q5jtmyBOXN02qLEhXBG6KuAOWY208zSgWuBlf2O+T1wrpmlmtl44Axgc2SjikRPZ3cnqSmpfYXeZ9s2r9BF4sCwI3TnXJeZ3QK8AASAx5xzG83spuD7DzvnNpvZn4H1QA/wc+fchmgGF4mkrp4u0gJpR6/jUlMDZ5/tVyyREQnrnlrOueeA5/rte7jf6x8CP4xcNJGx0dXTRWdPJ5mpmTQ1hZzh0tMDtbVQVORnPJGw6UpRSXoHDh0AIC8j78gpl7o67/Zzkyb5F05kBFTokvR2NOwAoCS35MhC37XL2+ocdIkTKnRJeu/t8841Xzhl4ZGF/u673nbRIn+CiYyQCl2S3s7GnQQswHG5M2luDin0DRu8hdGPO87XfCLhUqFL0mtsbyQ/M5+WFsO5kEJvaIDi4pCT0kVimwpdkl7jocYj1kLvK/QdO/SFqMQVFbokvYa2BgoyC44s9Lo6WLUKzj/f12wiI6FCl6S3qWYTpfmlRxb6M89456FffLGv2URGQoUuSW1v8152Ne3iI9M/cvhe0DkO7r/fu6lFebm/AUVGQIUuSW1n404AyorLDo/QO2pg40a47jp9ISpxRYUuSa2hrQGAgnEhc+gZ7d6TI+5DJxL7VOiS1EIv+z/gPSXPNXpPsrL8CSUySip0SWrNHc0AZKdn09QEqakwbndwGX8tmytxRoUuSa2lowWAnIycwze32LzJmztXoUucUaFLUust9Ky0rMPruLz/PsyerSkXiTsqdElq9W315GbkEkgJHC70Dz6Ak07yO5rIiKnQJantP7if4vHeDcsPHIDcXGDfPpg61d9gIqOgQpek9kH9B8wsmAngjdBzur3bzk2Z4nMykZFToUvS6nE9bKrZxLzieUCw0FMPgnMwc6bP6URGToUuSau1s5XWzlam5njTK01NkOeCVxdpDXSJQyp0SVoHOw4CkJWehXPeHHpej3flKNOn+5hMZHRU6JK0Oro7AMgIZHDwoHc/6L6rRAsL/QsmMkoqdElaXT1dAH2nLALkZQbXcQkEfEolMnoqdEla3a4bgIAdLvTc9hpIT4dx43xMJjI6KnRJWs45AMzs8Ah9zyYoK/MWdRGJMyp0SVptXW0AjEsdd7jQ67bDCSf4mEpk9FTokrR613HpXWkRIL95N0yY4GMqkdFToUvS6j1tMbTQ89qqtCiXxC0VuiStvpUW07MOF7od0G3nJG6p0CVp9Z9yCQQcWW01wSUXReKPCl2SVu/dinLSvZtb5GZ1YwBz5/qaS2S0wip0M1tiZlvNrMLM7hziuNPNrNvMropcRJHoaGxvBCA/M99bxyXzkPdGaalvmUSOxbCFbmYB4AHgIqAM+KyZlQ1y3PeBFyIdUiQa6lrrGJ82nozUDBobIS8jeJWoznKROBXOCH0xUOGc2+6c6wBWAJcPcNytwH8D+yOYTyRq9rfuZ2LWRMBbaTHfDniX/E+e7HMykdEJp9BLgN0hryuD+/qYWQnwKeDhoT7IzG40s9VmtrqmpmakWUUiqqqliklZkwC8Efqh/d6t5zIy/A0mMkrhFPpA53C5fq+XA193Lrg4xiCcc48458qdc+XFxcVhRhSJjm3125hVMAsIroXeUaMbW0hcC6fQK4HQxaGnAXv7HVMOrDCzncBVwINmdkUkAopEQ3tXOzsbd3JCoXeZf1MT5HXVQX6+v8FEjkE4KxCtAuaY2UxgD3At8LnQA5xzfcMaM3sC+KNz7neRiykSWZtrNuNwzC2eS08PNDU58qmEadP8jiYyasMWunOuy8xuwTt7JQA85pzbaGY3Bd8fct5cJBZtrt0MQFlxGS0t4JyRTz2cssjnZCKjF9Yaoc6554Dn+u0bsMidc0uPPZZIdLV2tgKQl5lHY6O3L48mXVQkcU1XikpSqm+rB7zbz/UWegGNWjpX4poKXZLSs5ueZf6k+RRnFR9eOje7S3cqkrimQpek80HdB6zau4rrTr0OoG+Enj++w79QIhGgQpek8+ymZwH4zNzPACGFXqTbzkl8U6FL0nm+4nnKp5YzPc+7vKLvS9Hpuf6FEokAFboknW0N2zhl4il9r/sKfULAn0AiEaJCl6RysOMg+5r3UZpf2revsRGy7CBpgR7fcolEggpdkso7e97B4Th96ul9+xoaID/lAKRqDl3imwpdksbLO17my3/8MhmBDM6aflbf/oYGKLBGSEvzL5xIBKjQJeE1tDVw/e+u5+NPfRyA5//pefIz8/veb2yEfBohPd2XfCKRon9jSsL79mvf5pfrf8myc5ex7NxljEs78uKhhgaY4RpU6BL3NEKXhLe/dT+l+aV89/zvHlXmEJxycfWacpG4p0KXhLe9YTu5GYOfY95X6PpSVOKcCl0S2lu73+Ltyre5fv71A77f1QXNzVDg6iAra4zTiUSWCl0S2pPrniQnPYcbFt4w4PuHV1psgLy8sQsmEgUqdElom2o2MX/yfLLTswd8/4hC19K5EudU6JKwunu6eX//+5QVlQ16TEODty2gAU47bYySiUSHCl0S1pqqNTS2N/Kx0o8Neky9d58LJqQfhKKiMUomEh0qdElYb+5+E4CPHTd8oRfkdI1FJJGoUqFLwtrRsIPxaeOZmjN10GN6p1wmFLgxSiUSPSp0SUgd3R08X/E8ZcVlmNmgx/WN0CcMfoxIvNCVFJKQfvLOT9hat5U/fvaPQx5XV+vIpoX0BYN/cSoSL1ToknDWVq3lG698g0tOuISLT7h4yGPrdx6gkAY466whjxOJB5pykYTS2N7IlU9fSeG4Qh697NFhj6+r7mIC9TBjxhikE4kujdAloTy76Vl2NO7gtaWvMTFr4rDH1zUYhdRBfuEYpBOJLo3QJaHsatxFiqVw9vSzwzq+rjHVK/Rc3SBa4p8KXRLKzqadlOSUkJoS3j8+a5vTKaIWcnKinEwk+lToklC21m7lxKITwzq2uxsa2zK8EboKXRKACl0SxvaG7aypWsOiKYvCOr6hAZwzCqmHcUff+EIk3qjQJWHc8+o9pKWkcdsZt4V1fG2tty3MPgRDXHwkEi/CKnQzW2JmW82swszuHOD9fzKz9cHHm2Y2P/JRRQb32JrH+MX6X3D7mbczJWdKWD/TW+jFE7qjmExk7Az7zZGZBYAHgE8AlcAqM1vpnNsUctgO4GPOuQYzuwh4BDgjGoFFejUfamZ99XreqnyLu166i0/M+gTfOu9bYf98TY23LZqtG1tIYgjnVIDFQIVzbjuAma0ALgf6Ct0592bI8W8D0yIZUpKbc449zXtYW7W277Gueh0V9RV9x5RPLeeZq58hLRD+jZ5rqrqBAMXzJkUhtcjYC6fQS4DdIa8rGXr0fQPw/EBvmNmNwI0AM3RlngzAOceupl28U/kO7+55l7XVXoHXt9X3HXP8hONZMHkBS+cvZf7k+SyYvICSnJIhF+EaSM2eDmAcxdMyIvynEPFHOIU+0H8lA641amb/gFfo5wz0vnPuEbzpGMrLy7VeqXDg0AFW7VnFO3ve4e3Kt3lnzzvsP7gfgMzUTE6ZeAqfPvnTLJi8gAWTF3DKxFPIyYjMKYa11d3kcIDM3PSIfJ6I38Ip9EpgesjracDe/geZ2anAz4GLnHN1kYkniWjVnlU8vPph3t7zNptrNuOC44MTC09kyfFLOLPkTM6YdganTDxlRFMoI7W/xpjIfkhXoUtiCKfQVwFzzGwmsAe4Fvhc6AFmNgP4DXCdc+7vEU8pCaGivoJlLy/j1xt/TV5GHh+Z8RGumXsNZ5ScweKSxRSMKxjTPNU1RjE1EAiM6e8ViZZhC90512VmtwAvAAHgMefcRjO7Kfj+w8A3gULgweA8Zpdzrjx6sSVeHDh0gNd3vc7KrSt5fO3jpAfS+cZHv8FXz/4quRn+rp9SXRtgNtWQossxJDGEteCFc+454Ll++x4Oef5F4IuRjSbxqLWzlTd3v8nLO17m5R0vs3rvarpdN5mpmdxw2g3c87F7wj5PPNqqawOcTTWkZPsdRSQitHyuHLP6tnqeWvcUv9vyO96qfIuO7g5SU1JZXLKYu865i/Nnns9Z088iMzXT76h9uruhtiHAJKohe7LfcUQiQoUuo+Kc443db/DTv/2UZzY+w6HuQ8yfNJ+vLP4K5888n3NmnBOxs1GioabGW8dlMlWQ/RG/44hEhApdRqSzu5Mn1z3J8reXs7FmI7kZudxw2g18ufzLnDrpVL/jhW3fPm87hX2QleVvGJEIUaFLWLp7ulmxYQX3vnYvFfUVLJyykJ9f+nOunXctWenxV4hVVd52MlUqdEkYKvQ44pyjs6eTju4OOruD2wFej/a93tdHvNfjbd/b9x6bazdz6qRTWXntSi454ZIRX5kZS3pH6JOp0lrokjCSrtB7XE9faY243EZbkhH6Xd0uuqsCGkZ6IJ30QDppgTRvm+Jti8YX8fRVT3NV2VWkWPyf5rdnj7edyl4VuiSMuCv0Nz58gx+99aNRl2y0SzHFUvpKsH8phr7ufT4udRx5GXlHHxtyzFCfN5LfNdznBFKS5wKbykqYmNVCxsEOyM/3O45IRMRdoR/sPMj2hu1HlNFQpRiN4lMpxr/KSpg2vh7S8nWlqCSMuCv0C2dfyIWzL/Q7hsS5ykooTauGCVo6VxJH/E+GioxCZSVM6/kQpsTGVasikRB3I3SRY9XaCvX1MC1rM8yZ43cckYjRCF2Szu7g7VqmHdwKJ5zgbxiRCFKhS9LZscPbzmQHzJvnbxiRCFKhS9LpLfRZbNeUiyQUFbokne3bISO1y7tKNDN2VoAUOVYqdEk6O3bAzNw6UgIpkOvvTTZEIkmFLkln2zbHrIMb4KKLdNm/JBQVuiQV5+DvW3o44dB6+Od/9juOSESp0CWpfPghtLYHOCl1G1x8sd9xRCJKhS5JZdMmbzt3YQZk616iklhU6JJU1q3ztifP04Jcknh06b8kled/f4hT2ULhPK3hIolHI3RJGnV18Ma7aVzKH+Dss/2OIxJxKnRJDs3NPP+5/0d3TwqXlbwHp53mdyKRiFOhS+JbvRrKyvjD/6QzeXwT5e8/DunpfqcSiTjNoUvie/BBOhpbeX78p7nmc6mkFPgdSCQ6NEKXxNfWxut5l9Lcmsqll/odRiR6VOiS2Do74Y03WJlyBZmZcMEFfgcSiR4VuiS2FStwu3fzh44LueACGD/e70Ai0aNCl8T23HNsmPhxdlaP57LL/A4jEl0qdElsW7fy29zrAbjkEp+ziERZWIVuZkvMbKuZVZjZnQO8b2b2k+D7681sYeSjioxQSwur1mfwf3dey0UXwRRdHCoJbthCN7MA8ABwEVAGfNbMyvoddhEwJ/i4EXgowjlFwtbZCRs3wj1f3sd53S8yubCLJ5/0O5VI9IVzHvpioMI5tx3AzFYAlwObQo65HHjKOeeAt80s38ymOOf2RTrwC99bzf/5Tn6kP1ZikHN29L5Bju0hhU6XSpvLpLYrny7SMGZzacqfeOiNCygujm5WkVgQTqGXALtDXlcCZ4RxTAlwRKGb2Y14I3hmzJgx0qwA5BZnUFZUM6qflfhjdnSFH13zkGI9pFoP4wIdFKUf4OTcPZxbtJnSi+fC7HHRDyoSA8Ip9IH+++n/X1k4x+CcewR4BKC8vHywwdaQzrrxFJ65cTQ/KSKS2ML5UrQSmB7yehqwdxTHiIhIFIVT6KuAOWY208zSgWuBlf2OWQn8S/BslzOBpmjMn4uIyOCGnXJxznWZ2S3AC0AAeMw5t9HMbgq+/zDwHPBJoAJoBT4fvcgiIjKQsFZbdM49h1faofseDnnugJsjG01EREZCV4qKiCQIFbqISIJQoYuIJAgVuohIgjDv+0wffrFZDbBrjH9tEVA7xr8zEuIxdzxmBuUea/GY2+/MxznnBlzMwrdC94OZrXbOlfudY6TiMXc8ZgblHmvxmDuWM2vKRUQkQajQRUQSRLIV+iN+BxileMwdj5lBucdaPOaO2cxJNYcuIpLIkm2ELiKSsFToIiIJIuEL3cy+E7xx9Voz+x8zmxry3l3BG1tvNbN/9DNnf2b2QzPbEsz+WzPLD3kvlnNfbWYbzazHzMr7vRezuWH4m6HHCjN7zMz2m9mGkH0TzOwvZvZBcFvgZ8b+zGy6mb1iZpuDfz9uC+6P9dyZZvauma0L5v5WcH9s5nbOJfQDyA15/hXg4eDzMmAdkAHMBLYBAb/zhmS9EEgNPv8+8P04yX0ycCLwKlAesj/WcweCmWYB6cGsZX7nGiTrR4GFwIaQfT8A7gw+v7P370usPIApwMLg8xzg78G/E7Ge24Ds4PM04B3gzFjNnfAjdOfcgZCXWRy+Nd7lwArn3CHn3A68tdwXj3W+wTjn/sc51xV8+TbeXaAg9nNvds5tHeCtmM5NyM3QnXMdQO/N0GOOc+51oL7f7suBJ4PPnwSuGMtMw3HO7XPOvRd83gxsxrvvcKznds65luDLtODDEaO5E77QAczse2a2G/gn4JvB3YPd2DoWfQF4Pvg8nnKHivXcsZ5vOJNc8C5hwe1En/MMysxKgdPwRrsxn9vMAma2FtgP/MU5F7O5E6LQzexFM9swwONyAOfcMufcdOCXwC29PzbAR43pOZzD5Q4eswzowssOcZJ7oB8bYF8snTMb6/kSgpllA/8N3N7vX88xyznX7ZxbgPev5MVmNs/nSIMK645Fsc45d0GYh/4X8CfgHmLgxtbD5Taz64FLgI+74GQdcZB7EL7nHkas5xtOtZlNcc7tM7MpeKPJmGJmaXhl/kvn3G+Cu2M+dy/nXKOZvQosIUZzJ8QIfShmNifk5WXAluDzlcC1ZpZhZjOBOcC7Y51vMGa2BPg6cJlzrjXkrZjOPYRYzx3OzdBj2Urg+uDz64Hf+5jlKGZmwKPAZufcfSFvxXru4t4zzMxsHHABXofEZm6/v5WN9gNvRLABWA/8ASgJeW8Z3pkNW4GL/M7aL3cF3pzu2uDj4TjJ/Sm80e4hoBp4IR5yB/N9Eu/si23AMr/zDJHzV8A+oDP4v/UNQCHwEvBBcDvB75z9Mp+DN4W1PuTv9CfjIPepwJpg7g3AN4P7YzK3Lv0XEUkQCT/lIiKSLFToIiIJQoUuIpIgVOgiIglChS4ikiBU6CIiCUKFLiKSIP4/2meQg8g7PmUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# one iteration with Cauchy distribution\n",
    "\n",
    "n=100\n",
    "alpha=0.05\n",
    "r=cauchy.rvs(size=n)\n",
    "epsilon=np.sqrt((1/(2*n))*np.log(2/alpha))\n",
    "\n",
    "F_n=lambda x: sum(r<x)/n\n",
    "L_n=lambda x: max(F_n(x)-epsilon,0)\n",
    "U_n=lambda x: min(F_n(x)+epsilon,1)\n",
    "\n",
    "xx=sorted(r)\n",
    "\n",
    "df=pd.DataFrame({\n",
    "    'x': xx,\n",
    "    'F_n': np.array(list(map(F_n,xx))),\n",
    "    'U_n': np.array(list(map(U_n,xx))),\n",
    "    'L_n': np.array(list(map(L_n,xx))),\n",
    "     'CDF': np.array(list(map(norm.cdf,xx)))   \n",
    "})\n",
    "df['in_bounds'] = (df['U_n'] >= df['CDF']) & (df['CDF'] >= df['L_n'])\n",
    "\n",
    "plt.plot('x','L_n',data=df, color='red')\n",
    "plt.plot('x','U_n',data=df,color='green')\n",
    "plt.plot('x','CDF',data=df, color='blue')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f0a9dc9e8be04e67ac4736fb5eb348f2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=1000.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Average fraction in bounds: 0.180\n"
     ]
    }
   ],
   "source": [
    "bounds=[]\n",
    "n=100\n",
    "alpha=0.05\n",
    "epsilon=np.sqrt((1/(2*n))*np.log(2/alpha))\n",
    "\n",
    "for k in tqdm(range(1000)):\n",
    "    r=cauchy.rvs(size=n)\n",
    "    \n",
    "    F_n = lambda x : sum(r < x) / n\n",
    "    L_n = lambda x : max(F_n(x) - epsilon, 0)\n",
    "    U_n = lambda x : min(F_n(x) + epsilon, 1)\n",
    "\n",
    "    # xx = sorted(r)\n",
    "    xx = r # No need to sort without plotting\n",
    "    \n",
    "    df = pd.DataFrame({\n",
    "        'x': xx, \n",
    "        'F_n': np.array(list(map(F_n, xx))), \n",
    "        'U_n': np.array(list(map(U_n, xx))), \n",
    "        'L_n': np.array(list(map(L_n, xx))), \n",
    "        'CDF': np.array(list(map(norm.cdf, xx)))\n",
    "    })\n",
    "    \n",
    "    all_in_bounds = ((df['U_n']>=df['CDF'])&(df['CDF']>=df['L_n'])).all()\n",
    "    bounds.append(all_in_bounds)\n",
    "\n",
    "print('Average fraction in bounds: %.3f' % np.array(bounds).mean())\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.7**.  Data on the magnitudes of earthquakes near Fiji are available on the course website.  \n",
    "\n",
    "- Estimate the CDF.\n",
    "- Compute and plot a 95% confidence envelope for F.\n",
    "- Find an approximate 95% confidence interval for F(4.9) - F(4.3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "data=pd.read_csv('fijiquakes.csv', sep=',')\n",
    "r=np.array(data['mag'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABMPUlEQVR4nO3dd3yN1x/A8c/JkoVE7L33iL1H7RF7770ppRTVokapUnvW+NkrRGKFUKtmlBJ7BREqRhBJZJ3fH49drdDEzfi++3pe7r3Pufd+z33029NzznOO0lojhBAiYTIzdQBCCCFijyR5IYRIwCTJCyFEAiZJXgghEjBJ8kIIkYBZmOqLU6ZMqbNmzWqqrxdCiHjpxIkT97XWqaJb3mRJPmvWrHh7e5vq64UQIl5SSt34mPLSXSOEEAmYJHkhhEjAJMkLIUQCZrI++fcJDw/Hz8+P0NBQU4cS71lbW5MxY0YsLS1NHYoQwoTiVJL38/MjadKkZM2aFaWUqcOJt7TWPHjwAD8/P7Jly2bqcIQQJvTB7hql1GKl1D2llM8/nFdKqRlKqStKqdNKqWKfGkxoaChOTk6S4P8jpRROTk7yf0RCiGj1yS8Fav/L+TpArhdHD2DufwlIEnzMkN9RCAHR6K7RWu9XSmX9lyINgWXaWLP4iFLKQSmVTmt9J6aCFCIx01pzPfA6x24fIyQ8xNThxCqtNZERkURGRBIRHkFEeASR4a8fR0S8/TwyPPL1axF/L/9yKXWtNVoDGvS7zzVvldHaKBT1xusv34fWxvuj+IfP5fV3Rr387pfvfX2uRpWafNurX2z/nEDM9MlnAG698dzvxWt/S/JKqR4YrX0yZ84cA18tRMITERXB6b9Oc/DmwVfHnaA40maKBMJi4Ah/8VnvOxKBe0cexKsk/75+gffuRKK1XgAsAChRokSc3K3E3NycQoUKvXru5uaGLL8gYtOzsGccvX30VUI/7HeYoLAgALIkz0LVbFWpkLkCZTOWxcHaIca+N/hZMP63/fH38+fO7TuvHj96+IjgZ8E8e/aM4GfBBAe9fhwWFhbtzzczt8LCwg4zc1uUskNjS5S2Q0daoyKtMMcCM8xf/WOhzLGwMMMCMyy0GRaYY4nCEjMsMcNKa6zR2GiNDRobHYUl6o1PePsfM8wxMzMHpcBMgTIDpVBKoRVgZjw3U+rtMmYvujtfPDce8+qcmTIzzpkZn6XMjPcrs5fvUZiZGWXMzP9+XpkpctbLHWPX8UNiIsn7AZneeJ4R8I+BzzUJGxsbTp06ZeowRAL2POI5Xte82H19NwdvHuSPO38QqSNRKAqnKUzHIh2pkLkC5TOVJ1PyTB/+wH/w4MEDzp07x82bN7l16xa3bt166/HDhw/fKq+UInnytNjYpMbCIilKpQPsUdhhY22Phbk9YWH2hIYaf4IdllhhjQU2mGOLObaYYQ/Yo0lhHkYys2DsVDA2UcFYRYZgHh6CiuZudMrcDHMbSyxsLLG0syJJUitsU9hgl8oGG0drrB2tsUlhg42jjfHY0QabFK8fJ0mWxEisiVxMJHl3oJ9Sag1QGngcE/3xAwdCTOdaZ2eYNi1mPmvp0qW4u7sTHBzM1atXady4MT/99FPMfLhIcELCQ/C86smGcxvwuOTBk+dPsLawpnSG0gyrMIwKmStQJmOZT2qpR0VF4evry6lTp946bt269VY5BwdHUqbMjJ1dJrJlK0fGjJl5+jQTAQGZCArKjNbpiAgMwyzwCTYqGCebYFJYB5PcMhh7s2BsrUKwNg/G0iIY8+eBEOL/uuP6HWaWZtiltMM2pS22KW2xcUrz6rFtSltsUthgZW+Fpa0llnaWWNn9/bG5lfkn/NLiXR9M8kqp1UAVIKVSyg8YBVgCaK3nAduAusAVIBjoHFvBfg4hISE4OzsDkC1bNjZt2vSPZU+dOsXJkydJkiQJefLkoX///mTK9OktL5GwBIUFsf3ydjac38DWS1t5Fv6MFDYpaJ6/OU3zNaVqtqoksUjyUZ8ZGhrK2bNn30rmp0+f5smTJ4DR3Zg7d15y565EgQLOBAUV5P79LNy6lYnAQHsCAwE0yVQQhdLco2jyANJmDsA+ZA8EBBAZ9GLarcb4tzkYlLl6K0Hbpkz5InHbvPP668PK3kpmeMUR0Zld0/oD5zXQN8YieiGmWtwf62O6a6pVq0by5MkByJ8/Pzdu3JAkn8g9ef6ELZe2sOHcBrZf2U5oRCip7VLTvnB7muZvSuUslbE0j/5dyCEhIRw4cABPT092796Nj48PkZHG6KS9vT1FihShXr32WFo68/ixM5cvF+D8eRvOnwfQFMr6lKKpA2hQ3AfH8HtYBt7n+e0Awp+Ewl3gLtiksCFVgVSkqlmAVAVS4ZDVAbtUr1vhSZInkYQdj8WpO17jmyRJXrfCzM3NiYiIMGE0wlSCwoJwPeeK63lXPK96EhYZRvqk6elerDtN8zWlQuYKmJtFr+tBa8358+fx9PTE09OTffv2ERoaipWVFRUrVuTLL4eRJImR0C9ezI63txm//26819FBU7nwI1rUP0/KoBuEXvAlyPcJ+BrnzVPa4lQgFamqFDSSen7jsEttJ0k8AZMkL8Qn8n/qz8yjM5l3Yh6BoYFkTp6ZviX70ix/M8pkLGPMwoiGhw8fsnv37leJ3c/PD4C8efPSo0dPMmSoxa1blfH0tGX3buM95uZQpLCmQ70H5LW+gd19Xx7+4cvT/U8BCEltR9YqWclUIRNpCqcxknkqu1j5HUTcJkleiI905q8zTDk8hVVnVhGpI2mSrwkDSw+kXKZy0WoRR0ZGcvz4cbZv346npyfHjx8nKiqK5MmTU716dYYN+54kSWrx+++ZWbUK7t+HJEmgWlVN10YPyYIvZjdu4HfQl6cnn3IfCEljJPWXh1OeBLw8SFgYBAVF74iIMO5Giop6eadS3HjcrBl07PhZfi5J8u8ICgqKVrlOnTrRqVOnV8+3bNkSSxGJuEBrjdc1L34+/DM7r+7E1tKWXiV6MbDMQLI7Zv/g+588ecKuXbvYsmULW7duJSAgADMzM0qWLMnIkSMpXboWd+6UYssWC4YMgZAQcHCAenU1NQvextHPhyvu53my/QkXAfu09mStkpUsVbKQtXIcTOrPn0NgoHE8fvz68ZuvhYb++/H8+dvPg4ONxB0e/nGxvJgPb8yFjyOPXwyUfw6S5IX4F2GRYazxWcOUw1M4/ddp0tqnZULVCfQs0ZMUNin+9b3Xrl1jy5YteHh4sG/fPsLDw3F0dKROnTq4uLiQN28t9u1LgZsbjBtnNPQyZoSuXaFG4b+wueLDufU+XF8VyE0rc3LWzknFkRWNlnruWE7qoaHvT87RPT60OJ6ZGdjYgLX1+w8bG3B0/PtrSZOCvT3Y2Rl//tthZwcWkuLkF/gAT09Pvvnmm7de+9DUShH/BYYGsuDEAqYfnY7/U38KpCrAkoZLaF2w9T9Oe4yIiODw4cN4eHiwZcsWzhtTXMiXLx8DBw7ExcUFR8dybN5swaRJ8OefxvsKFYJvv4XqxR6iz/hwdq0PJ2cFoMwV2atlp9J3lcjXOB/WDtYxV8HgYLh2Da5cgatXjT9fHnfuGK3of2NhYSRhB4fXR6ZMkDz526/902Fra7RqRayTJP8BtWrVolatWqYOQ3wmgaGBjN03lgV/LCAoLIjq2auzqMEiauWo9d6Wc1RUFJ6enqxYsYLt27fz6NEjLC0tqVy5Mj179qRePRdCQnKwYQP07g3nzhnvK1cOfv4ZapZ5QtDRs/is8eG3scaN4pkrZKbu7Lrkb5Yfu9T/YbD0+XM4f/7tBP7yuH377bIpUkDOnFC2LGTI8PcE/u5hYyNJOp6QJC/ECzuu7KCbezfuBt2ldaHWDC47GOe0zu8t+/DhQ5YsWcLcuXO5evUqTk5ONGjQABcXF2rUqMnVq8nYsAHq1YNLl4zeiYoVYeZMqFMlhEe/n8VntQ8bh9wADemKp6PG5BoUaFmA5JmSf3zwWsPNm3DkyOvjjz+MQcqX0qSBHDmgWjUjob88cuQwkrxIkCTJi0TvcehjBu8czKKTi8ifKj+bWm6iZIaS7y178uRJZs+ezapVqwgJCaF8+fKMHTuWJk2acuqUFa6uMHQoXL9uTHP84gsYNAhcaoUTeOwSZ1aeYdWgy0SFR5Eyb0qqjK5CwVYFccrt9HFBBwWBt/frhH70KNy9a5yztoYSJeDLL40/8+QxEnnSpP/xlxLxkSR5kajtvLqTru5d8X/qz7DywxhVZRTWFm/3fT9//pwNGzYwe/ZsDh8+jK2tLe3ataNPnz5kyODMnDmQKxfcugWWllC9OowcCfVdonh62pczK8+wbMg5wp6GkTR9Ukp/WZpCbQuR1jlt9AdPAwJg+3Y4dMhI6mfOGCO1YHx5jRpQpoxxFCpkBCIEkuRFIvXk+RO+3vk1C/9YSN6UeTnU5RClM5Z+q8ytW7eYP38+Cxcu5N69e+TMmZOpU6fSqVMnAgIcmToV/vc/YyJJ7dowfjy4uGhCr9/l9IrTLBvpQ9CdIJIkS0L+5vkp3LYwWSpnwcw8ejdJcfUqbN4Mbm7w++9GUk+eHEqXhoYNjYReqhQ4feT/BYhERZL8O3x9fXFxccHH5/WWtqNHj8be3p6vv/7ahJGJmOJ1zYuu7l3xe+LHkHJD+OGLH95qve/du5cZM2awefNmtNa4uLjQt29fqlevwaFDZnTuDO7uYGUFHTrAV19B2iSPOLPqDKvKneH+hfuYWZqRu15uCrUtRK56ubC0iUbLWmujH93NzThe/h0sUgS++85I7EWKGB38QkSTJHmRaDx9/pShu4Yy78Q8cjvl5mDng5TNVPbV+du3b9O/f382bdqEk5MTQ4YMoVevXmTMmJVNm4wZMceOGQ3n776DPn0g7OotDnx9gMvbLgOQpVIWygwqQ/5m+bFxtPlwUOHhsG+fkdQ3bwY/PyOJV6oEv/xiJPZs2WLpFxGJQZxN8gN3DOTU3VMx+pnOaZ2ZVntajH5mlSpVKF26NL/99huBgYEsWrSIihUrxuh3iP9uz/U9dHXvyo3AGwwuO5ixX4zFxtJIwpGRkcybN4/hw4cTHh7OxIkTGTBgABER1ixebORaX19jIsqcOdChg+beUV92tN6P72++2Ka0pcoPVXDu6EzyzNGYGfPsmdG/vmkTbN1q3HRkYwO1ahl3RdWrBylTxu4PIhKNOJvk45OIiAiOHTvGtm3bGDNmDF5eXqYOSbwQFBbEN7u+YY73HHKlyMWBzgcon7n8q/NnzpyhR48eHDlyhBo1ajB37lxsbXPwww8wd65x82b58kaid3HRXN91hbU1D3Dr0C3s09lTc2pNivcojpWd1b8H8vQpbNkCrq6wbZuxboGTEzRpAo0aGaO1trax+luIxCnOJvmYbnFH1z/Ndvi3WRBNmjQBoHjx4vj6+sZGWOITePt702pDK649usZXZb5iXNVx2FoaiTQkJISxY8cyefJkHBwcWL58OS1btmXqVMWoUUYvSpMmMHgwlC6lubD5AotLH+DOH3dInjk5defUpWjnolhY/8u/QoGBRue9qyt4eho3J6VNC126QNOmxsR5ue1exDL5G/YOJycnHj169NZrDx8+JNu/9Iu+XFde1pSPG6J0FNOOTGOY1zDS2qdlb6e9VMpS6dV5Ly8vevXqxdWrV+nUqRM///wzAQFOVKxoTDdv0gQmT4asWaI4u+4sc7sdIOBsAClypqDBogYUblf4n7eme/DA6FvfsAG8vIz/WmTMaNzu2qyZcUepDJyKz0iS/Dvs7e1Jly4du3fvplq1ajx8+JAdO3YwYMAAU4cmouHes3t0cuvE9ivbaZS3EYsaLHq1kFhAQACDBw9m+fLl5MqVi927d1O5clWmTjUGUu3sYPVqaNYkkjMrTzP7x4M8vPyQVPlT0WRlEwq0KICZxXsSdGAgrF1rtNj37IHISGOwdOBAo8VesqQkdmEykuTfY9myZfTt25fBgwcDMGrUKHLkyGHiqMSH7L62m3ab2vEo5BGz686md4neKKXQWrNs2TIGDx7M48ePGTlyJN9++y03blhTsSIcPmx0i8+dCyHnrjM7jzuBvoGkLZqWFq4tyNsoL8rsPd11vr7GPpWLFhl3oObKZdzu2qwZFC0qa7uIuEFrbZKjePHi+l3nzp3722vi0yWW3zM8MlyP8Bqh1Wil887Kq/+8++erc5cvX9ZVq1bVgC5Xrpz28fHRERFaT5mitbW11o6OWq9cqXVYSLj2/NpTj1aj9czcM/XFLRd1VFTU+7/w2DGtW7bU2sxMawsLrdu109rbW+t/Ki9EDAK89UfkWmnJi3jtRuAN2mxsw6Fbh+hatCvTa0/HzsoOrTXz589n0KBBWFpaMnfuXHr06MHVq2ZUrmzcQFq/PsyfDxaBASwuu5G7p+5SvFdxav5c8++zZaKijOmOP/8M+/dDsmTGqOyXXxp97kLEUZLko6lv3778/nLH5BcGDBhA586dTRSRcD3nSjePbkRGRbK66WpaFWwFwN27d+nWrRtbt26lZs2aLF68mHTpMjBzJgwfbmylt2wZtG2r8Z57nF1f78LK3opWm1uRp0Get78kNBSWL4cpU+DiRcicGaZONXb2SJbMBLUW4uNIko+m2bNnmzoE8UJIeAiDPAcx78Q8SmUoxeqmq19twbd582a6detGUFAQM2fOpG/fvly9qqhSBQ4cMO4zWrAAkpkHsaa+O5e3XSZn7Zw0XNIQ+7T2r7/k/n3jzqdZs4zFwYoVezEq20ymPYp4Rf62injl7L2ztHJthc89H4aUG8K4quOwMrfi6dOnfPXVVyxatIiiRYuycuVK8ubNx5w5MGSIsc7MkiXG3smXt11iZefNPH/ynNozalOqX6nX90E8eADffw+LFxuteBcXo1umcmUZSBXxkiR5EW8s/3M5Pbf0JGmSpOxou4NaOY0duw4fPky7du3w9fVl+PDhjB49mqdPrWjUyLgXqXZt+PVXSO0YzrZ+O/Ge402awmnouKcjqQumfv0FmzdDz55Gou/UyVgIPl8+k9RViJgiSV7EeVE6ipF7RvLjwR+pkrUKq5uuJq19WsLDw/nhhx+YMGECmTNnZt++fVSoUIF9+6BtW7h3z5jh+OWXcPfUHRZU38j9C/cpM6gM1cZXe3236sOHRqGVK8HZGXbuhMKFTVllIWKMJHkRpz0Le0YHtw5sPL+RHsV6MKvuLCzNLbl48SLt2rXD29ubTp06MX36dOzskjF6NIwda2yEdOQIFHXWHJp8iD0j92Cb0pZ2O9uRo8Yb9zy4uxut9/v3YfRoGDFCNtwQCYok+XfY29sTFBRk6jAEcPvJbRqsacCpu6f4pdYvDCht3HU8d+5cBg8ejK2tLa6urjRp0gQ/P2NK5P790L49zJ4NKjiI5TU3cn33dfI2zkv9hfWxdXqxCNijRzBggDFzpnBhY1VIZ2fTVVaIWCJJXsRJJ/xP0GBNA548f4J7K3fq5a7HX3/9RZcuXdi2bRu1a9d+MTUyHe7u0Lmzsf7XsmVGkr9x4AYbWm4g9FEo9RfWp2jXoq8HV7duhR494K+/jEHWb781RmaFSIDibpIfOBBOnYrZz3R2NjppY0CnTp1IliwZ3t7e3L17l59++olmzZrFyGcndq7nXGm/qT2p7VJzqMshCqUphLe3N40aNeLBgwfMmjWLPn368Py54ssvYeZMY4bjmjWQM6fm98mH2D18N47ZHGm7vS1pi6Q1Pjgw0NjGaelSYx9UDw/jjUIkYHE3yccDd+7c4eDBg1y4cIEGDRpIkv+PtNb8ePBHvt3zLWUzlsWtlRup7VKzYsUKunfvTpo0aThy5AhFihTh4kVo1cpoBwwcCBMngg4JZW1jNy5uvki+pvlosKgB1slfbOu3Ywd06wZ37xot9+++M+6KEiKhi87aB0Bt4CJwBRj2nvPJAQ/gT+As0PlDnxlX166xs7OLVrmOHTvqFStWvHpub28fWyF9srjwe0ZXaHiobr+xvWY0uo1rGx0SHqIjIiL0119/rQFdpUoVfe/ePR0VpfWSJVrb2Wnt5KS1h4fxfv8T/np69un6B4sf9OFfDr9edyYwUOsuXbQGrQsU0Pr4cZPVUYiYQEyvXaOUMgdmAzUAP+C4Uspda33ujWJ9gXNa6/pKqVTARaXUSq11WEz9xyguSvJGS9D47cWnuPfsHo3XNubQrUOM/WIs31b8lsDAQBq1boSnpyf9+vVj6tSphIZa0r69MdOxShVYsQLSp9ecWPgH2/tvxy6VHZ32dSJTuUzGBx86BK1bG/umDh8Oo0ZJ610kOtHprikFXNFaXwNQSq0BGgJvJnkNJFXGyJY98BCQ3TPEB/nc86H+6vrcDbrLumbraF6gOefOnaNhw4bcuHGDhQsX0q1bN65fN2bPnD9vTJEcPhwiQ8PY3Gkbfy77kxw1c9B4RWPsUtkZ67lPmmQMqmbJYiT70qVNXVUhTCI6ST4DcOuN537Au//GzALcAX8gKdBSax317gcppXoAPQAyZ878KfHGuuDgYDK+sargoEGDGDRokAkjSri2X95Oyw0tsbeyZ3+n/ZTMUBIPDw/atm2Lra0te/fupVy5chw6ZKz3Hh5u3KdUrRrcv3if9c3Wc+/sPSqPrkylkZUwMzcDf39jes2ePUYrft48WUhMJG4f6s8BmgO/vvG8PTDznTLNgF8ABeQErgPJ/u1z42qffEISl3/PladXarMxZtp5nrO+9fiWjoqK0uPGjdNKKV2iRAl969YtrbXWy5drbWWlda5cWl+4YLzXZ62PnmA/Qf+U8id9ZeeV1x+6davWKVNqbWur9eLFsr67SJD4yD756OxJ5gdkeuN5RowW+5s6AxtfxHDlRZLP+6n/4REJ2xqfNbTf1J5KWSpxoPMBHMwcaNGiBSNHjqRt27bs37+f9OkzMnKk0SgvX964ezVntki2D9jOhpYbSF0oNT1P9jTuXg0LMxYRq1cP0qeHEyeMifOyoJgQ0equOQ7kUkplA24DrYA275S5CVQDDiil0gB5gGsxGaipjB8/nvXr17/1WvPmzfn2229NFFH8tv7setptbEeFzBXY0noL927fo1GjRvj4+PDzzz8zaNAgQkIULVsae2F362bcvRr5LITlNddyY98NynxVhuqTqmNuaQ5XrhjdMt7e0LevsamHtbWpqylE3BGd5j5QF7gEXAW+ffFaL6DXi8fpgZ3AGcAHaPehz5TumtgX135P13Ou2nyMuS6/qLx++vyp3rNnj3ZyctIODg56x44dWmutb9/WunhxrZUytuiLitL64bWHelbeWXqs1Vh9euXp1x+4YoXW9vbGHn6bNpmmUkJ8ZsTG9n9a623Atndem/fGY3+g5n/8741IwNwuuNFyQ0tKZSjF9rbb2bppK23btiV37ty4u7uTM2dO/vgDGjSAx4+NVX/r14fbx2+z2mU1keGRtN/VniyVshibZvfvb9y5WqGCMacyjg7kC2Fq0emTF+I/8bjoQYv1LSierjg72u3Aw9WDNm3aUL58eY4cOULOnDnZuBEqVgRz89f7r17YfIGllZdiaWdJ10NdjQR/6hSUKAH/+58xRfK33yTBC/EvJMmLWLX10laarmuKc1pnPNt5snndZtq3b0+lSpXYtm0bSZMmY+JEaNrUWAzy6NEXf848ytrGa0lTKA3djnQjZd6Uxs4fpUvD06fGFMkxY2QrPiE+QP4NEbFmx5UdNFnXhEJpCuHZzpNNazbRpUsXqlatiru7O+bmtnTqZKwc2bq1seOelWUUO77aydFpR8nbKC9NVjbBMomZsUvTL79ArVrGra4pU5q6ekLEC9KSj2OyZs3K/fv3//G8vb39P56LS3Ze3UmjNY3Inyo/u9rvYtNqI8FXr14dDw8PwsJsqV7dSPBjxhjd6uZR4axvvp6j045SekBpmm9ojmVkKDRsaCT4L7+ELVskwQvxEaQlH4MiIiKwkO4DvK550XBNQ/KmzItXey9cV7jSo0cPatWqxaZNm3j+3IaaNeHPP43lgVu2hGf3nrG6wWpuH7tNrWm1KDOgDNy4YXTOnzsHc+dCr16mrpoQ8U6czUg7Bu7g7qm7MfqZaZ3TUnta7X8t4+vrS506dahQoQKHDh0iQ4YMbN68mYsXL9KrVy+Cg4PJkSMHixcvxtHRkSpVqlCuXDl+//13GjRogIeHB0WLFuXEiRMEBASwbNkyfvzxR86cOUPLli0ZN24cAI0aNeLWrVuEhoYyYMAAevToEaN1NZU91/fQYHUDcqXIhVcHL9YvW0/v3r2pW7curq6uhIVZU7u2MX7q6mrk8AeXHrCyzkqe3nlKy40tydsoLxw+bKxl8Py5sUxw9eqmrpoQ8ZJ017zH5cuX6du3L2fPnsXBwQFXV1c6dOjApEmTOH36NIUKFWLMmDGvygcGBrJv3z4GDx4MgJWVFfv376dXr140bNiQ2bNn4+Pjw9KlS3nw4AEAixcv5sSJE3h7ezNjxoxXr8dn+3z34bLKhWyO2fDq4MXaJWvp3bs3Li4ubNy4kfBwI8GfOAHr1xsJ/saBGywqu4iwoDA67e1kJPiVK+GLLyBpUuNWV0nwQnyyONuS/1CLOzZly5YN5xf7fRYvXpyrV68SGBhI5cqVAejYsSPNmzd/Vb5ly5Zvvb9BgwYAFCpUiAIFCpAuXToAsmfPzq1bt3BycmLGjBls2rQJgFu3bnH58mWcnJxiu2qx5sCNA9RbVY+sDlnZ02EPaxatYcCAATRs2JB169bx/LkVderAsWOwbp3Rze6z1ge3Dm44ZHOg7ba2OGZNbkyLHDsWKlc2mvrx+DcRIi6Is0nelN5cJ97c3JzAwMB/LW9nZ/fe95uZmb31WWZmZkRERLB37168vLw4fPgwtra2VKlShdDQ0JirwGd2N+guDdc0JGOyjOzpuIdVC1cxaNAgGjduzJo1awgPt6JePaNRvmYNNGkCJxefxL2rO1kqZaHlppbYWGtjq6f166FLF6MPXvZdFeI/k+6aaEiePDmOjo4cOHAAgOXLl79q1X+Kx48f4+joiK2tLRcuXODIkSMxFepnp7Wmz9Y+BIcH49bKjRXzVjBo0CCaNWvG2rVriYiwwsXFuMFp5Upo1gxOrzyNezd3ctbOSbud7bB5HmjsArJhA0yebMyHlwQvRIyQlnw0/e9//3s18Jo9e3aWLFnyyZ9Vu3Zt5s2bR+HChcmTJw9lypSJ9nvj2nr3G85tYNOFTUysNpHNizYzbNgwWrZsyYoVKwgLs6B+fdi/H5YvN2bRnHM9h1tHN7JWyUqLjS2wOHfGWMvg0SNwczMeCyFizscsdBOThyxQFvti+/cMeBagU09OrYvPL65/nPijBnSbNm10eHi4Dg7Wunp1Y6Gx5cuN8he3XNQ/WP6gF5VfpJ8/fW4sKmZrq3WmTFqfOhWrsQqRUBAL68kL8V4DdwzkYchD+qfrz4jhI2jVqhXLli0jIsKCRo1g925YsgTatYNrXtdY13QdaYukpc3WNlitX2l0zhcsaIzGFili6uoIkSBJd00c9ODBA6pVq/a313fv3h1nZuBsubSFlWdWMqLMCMb3H0+WLFlYuHAhERHmNGlibNO3aBF07Ag39t9gdYPVpMyTknae7bDevA66doUaNYwuGhsbU1dHiARLknwc5OTkxKlTp0wdxj8KDA2k55aeFExdkGc7nnH58mV2796NpaU9zZrB9u2wYIExScbvqB+r6q3CIasD7Xe1x2b7RujUydioVRK8ELFOumvERxuycwh3g+4yMP1AZkyfQe/evalQoSotWhhLy8ydC927w52Td1hZeyV2aezo4NUBu93u0KGDMZNm82ZJ8EJ8BpLkxUfxuubFryd/ZUCxAfw09CcyZ87MxImTaN8e3N1h1ixjiZl7PvdYXmM5SZIlocPuDiQ9uN3onK9YETw8wNbW1FURIlGQ7hoRbUFhQXT36E5up9zo3zSXLl1i9+7duLsnZd06GD/e2Gb1waUHLKu+DHMrczrs6YDDcS9o08bYkXvLFnjn5jEhROyRJC+ibcTuEdwIvMHcwnPp/WVvevXqRYECVWneHMqWhW++gUfXH/G/qv9DR2k6/taRFKf3GYvFlykDW7dCPFkqWYiEQrpr4ph/W0/e19eXggULvvXa6NGj+fnnn2M9roM3DzLz2Ex6F+nNL8N/IVOmTPz000/06wfPnhkbfgT5P2ZZ1WVEhETQwasDqS4fhhYtoGRJ2LbNWHBMCPFZxdmW/MCBA2N8homzszPTpk2L0c98U0JdTz4kPISu7l3J6pAVywOWXLx4kV27drFzZ1I2bIAJEyBD8qcsrbyMkIchdNjdgTQ3jxtrGBQrZky3SZbM1NUQIlGSlvw7fH19yZcvH927d6dAgQLUrFmTkJAQTp06RZkyZShcuDCNGzfm0aNHAFSpUoURI0ZQuXJlpk+fTpUqVfjqq6+oVKkS+fLl4/jx4zRp0oRcuXIxcuTIV9/TqFEjihcvToECBViwYIGpqhsto/eO5tKDSwzKOIiZ02bSo0cPihatTt++Rg7/6stIVtVdxVP/p7Td3pb0904Zm7YWKQKenpA8uamrIETi9TG3x8bkEVeXNbh+/bo2NzfXJ0+e1Fpr3bx5c718+XJdqFAhvXfvXq211t99950eMGCA1lrrypUr6969e796f+XKlfXQoUO11lpPmzZNp0uXTvv7++vQ0FCdIUMGff/+fa211g8ePNBaax0cHKwLFCjw6vUsWbLogICAf4ytQIECb702atQoPXny5PeWj4nf85jfMW02xkx33tBZ582bV2fKlEk/fvxYt2+vtYWFsRrB3jF79WhG6/Nu57XesUPrJEm0LlZM64cP//P3CyHehixr8N9FZz35/fv3vyofnfXkkyRJ8mo9eYAZM2ZQpEgRypQp82o9+Q9RSn3U6/9VWGQYXdy7kM4+HckOJePChQv8+uuvHDyYjOXLYcQISG8ZwP5x+ynYqiB5bW8ZC8Xnywe7doGjY6zEJYSIvoTXgRwD4up68k5OTq+6iV56+PAh2bJl++B7P8WEAxPwuefDlPxTGDJkCN27d6d06ZoULAgFCsDwYVGsrOpOkqRJqN3MzlhBMk8e8PKCFCliJSYhxMeRlnw0xJX15O3t7UmXLh27d+8GjAS/Y8cOKlSo8Mmx/JPTf51m/IHxtM7bml9H/UqGDBn4+eefGToU/P2N2TR/LjyO3xE/ao8qhV3XVpAzp5Hg48j6OkIIaclHW1xZT37ZsmX07dv31X6yo0aNIkeOHJ8cy/tEREXQZXMXUtikIPXx1Jw/f54dO3bg7Z2MBQtgyBDInTqQOSN2k7NOTgrtnAJhYbBpE6RKFaOxCCH+o4/pwI/JI64OvCYkn/p7Tjo4STMaPXHNRG1mZqa7du2qg4K0zpZN61y5tH72LEovr7lcj7cbrx9NW6o1aD1lSgxHL4R4Hz5y4FVa8uItVx5eYdTeUTTM0ZBlPywjffr0TJkyhW+/hevXjV2eLm84zdWdV6k9vgIOY1tC6dIwYICpQxdCvIck+TjIVOvJa63ptaUXVuZWZDqVic3nNrNt2zZ8fJIzYwb06wfF8jxjdiNPMpXLRMmTC+DpU6OD3tw81uISQny6aCV5pVRtYDpgDvyqtZ74njJVgGmAJXBfa/1JI5Na61ibEhhfxMR68sb/1X2cZX8uY/f13QzLNoyfvv2Jzp0788UXdXB2hsyZ4ccfYXu37YQFhVG/lQ1mX75YlSx//v8UqxAi9nwwySulzIHZQA3ADziulHLXWp97o4wDMAeorbW+qZRK/SnBWFtb8+DBA5ycnBJ9ov8vtNY8ePAAa2vraL/n3rN7DNo5iHLpy7H1562kTZuWqVOnMmYMXLxo7PR0e89Fzq49S5XhpUk1oTM4OxujsEKIOCs6LflSwBWt9TUApdQaoCFw7o0ybYCNWuubAFrre58STMaMGfHz8yMgIOBT3i7eYG1tTcaMGaNd/ivPr3j6/CklfUsy/cx03NzcuHLFgcmTjZ36KpYKZU7+raQulJoKt9bA/fvGmjSWlrFYCyHEfxWdJJ8BuPXGcz+g9DtlcgOWSqm9QFJgutZ62bsfpJTqAfQAyJw589++yNLSMtZu7BH/bPvl7aw6s4p+2fsxr/s8mjZtSp06DSlRAtKkgZ9/Bq9vvAi6G0TLoZkwH/g/GDnSaMkLIeK06CT59/WbvNvhawEUB6oBNsBhpdQRrfWlt96k9QJgAUCJEiU+vtNYxLhnYc/ovbU3eZ3ycnrRaaytrZk5cyYTJ8KZM8YmToF/+nJi/gnK9CtGhp/7GH3wbyy2JoSIu6Jzx6sfkOmN5xkB//eU2aG1fqa1vg/sB4rETIgiNn3/2/fceHyDhsEN2b9vP5MnT+bBg3SMG2ds5lSrWjge3TxwyObAFyHbX9/u+sZyDUKIuCs6Sf44kEsplU0pZQW0AtzfKbMZqKiUslBK2WJ055yP2VBFTPP292ba0Wm0z96e+T/Op1KlSnTq1JWuXcHBAaZPh31j9vHwykPq982E1aK58NVXxrx4IUS88MHuGq11hFKqH+CJMYVysdb6rFKq14vz87TW55VSO4DTQBTGNEuf2Axc/DfhkeF09+hOGrs0PHV7SkhICAsWLGDuXDOOHYNVqyD85h0O/XwI5w4FyT5niLE2zQ8/mDp0IcRHiNY8ea31NmDbO6/Ne+f5ZGByzIUmYtO0I9M4dfcUw1IOY+LGiYwbNw4bmzyMGAF16kDzppEsKu2OXSo7atr9Dteuwb59YGtr6tCFEB9BVqFMhK49usaovaOom6kuy39cTsGCBfn66yH06QNaw5w5cGTqYe6eukvdL3NgM28a9OkDlSqZOnQhxEeSZQ0SGa01Pbf0xMLMglRHUuHv74+rqyubN1uxdStMnQpJwx6wd/Re8jXMTb5lwyFTJpj4t5uchRDxgCT5RGb56eV4XfNicMbBTP12Kv379yd37tI0bAglSkD/fpoVNTywsLagThYf2HzB2Kc1aVJThy6E+ASS5BORgGcBDPIcROm0pdk+dTsZM2Zk3LhxfPWVcQPrjh1wetlJbuy7Qf2RhUn6Y3Po3Blq1jR16EKITyR98onIoJ2DePL8CSWuleDcuXPMnTuXEyeSsmgRDB4MOdMFsevrXWSpmImi7mMgdWqYMsXUYQsh/gNpyScSnlc8WXF6Bb2y9GJhz4W0atWKatXqUaQIZM8Oo0bB9q6ehAeH41L0NmrGn7B5s2zGLUQ8J0k+EXgW9oxeW3uR2zE3p389jZ2dHdOmTWPcOLh06cUKk3sv47PGhyr9C5JyXhto2dLYmFsIEa9Jd00iMHrvaHwDfakfVJ9Dvx9iypQpBASkYdIk6NABKpcNY2ufraTMl5LyJ2cac+GnTzd12EKIGCBJPoE7dfcUU49MpU2WNiyctJCqVavSvn0nunc3li6YMgV++/43Ht94TH0XjcXBfcayk2nSmDp0IUQMkO6aBCxKR9Fnax+cbJx4sukJYWFhzJ8/n3nzFEeOwPLlEObrz9HpRynePh+ZF/SAKlWgSxdThy6EiCHSkk/Alv25jMN+h2lh1oIt7lsYPXo0SZLkZPhwqFULWreMwqO7B3Zp7Kj+ZCOEhsL8+SC7cgmRYEiST6AehTxi6K6hlExRko1TNuLs7MxXXw2ib1+IjIS5c+Ho9CPcPXWXOh1TYb15LXz3HeTOberQhRAxSJJ8AvXdb9/xIOQBOU/n5O7duyxcuBB3d0s8PIyFJB14xG/f/0aeujnIt/xbKFhQ9msVIgGSPvkE6OSdk8z1nkuLlC1Y98M6evbsSY4cJahfH4oVgwEDNGvrb8XM3Iw6qb1R/rfBdQNYWZk6dCFEDJOWfAITpaPou60vKaxTcHP1TRwcHBg/fjzffAP37sHChXB+3Rmuel6lao/sJP/fDOjXTzYCESKBkiSfwLwcbG0a2ZRDBw8xYcIEfHxSsHChsalT3izBeA70JEPJdJT0HAcZMsD48aYOWwgRS6S7JgF5OdhaKmUpPCZ4UKJECdq160qxYpA1K4wZA7v67SI0MJT6ZcIwm+kD7u6ywqQQCZgk+QTk5WBrnet1OOZ/jI0bNzJ1qjkXLxorTN47ep1TS09Rvkc+0izoAM2bQ/36pg5bCBGLJMknEC8HW1ulbsWqsavo2rUradKUZvx4aNYMqlYKZ17hLTjmcKTy+flgbS1LFwiRCEiffALw5mCr/3p/7O3t+fHHHxk4EMzMjN2e9o/bz8MrD3FpZI7lgT0weTKkS2fq0IUQsUySfALwv1P/M+5sVS3Yu2cv48aN4/jxVGzeDN9/D1aBf3Hop0MUaZGb7ItGGnu1du1q6rCFEJ+B0lqb5ItLlCihvb29TfLdCcmjkEfkmZWH7PbZ8f/RnxQpUnDwoDfOzhZYWMCpk1Gs+GIxj64+om+FU9hu2wB//gl585o6dCHEJ1BKndBal4hueWnJx3MvB1sLXijIrVu3mDVrFtOmWXD1KsycCWeWnuD20dvU7JAGW7dVMHKkJHghEhFJ8vHYy8HWNunasHzectq3b0/GjBVeDbaWc37GnhF7yFopE4XXjYT8+eGbb0wdthDiM5LZNfHUm4Otd13vkiRJEiZNmkTv3q8HW72+8SIsKIy6mX1QB/xg7UFZukCIREZa8vHUy8HWNlZt8PL0YsyYMZw8me7VYKu+eZNTS05RpnVWUq2cBr17Q7lypg5bCPGZycBrPPRqsDVpdv6a9Bd2dnYcPnySokUtjcHWP6JYWm4BIQ9C6Jt6HVb+vnDhAiRPburQhRD/kQy8JgIvB1uLXi2Kr68vM2fOZPp0y1eDrX/+epy//vyLWvXMsfrjiLGdnyR4IRIl6ZOPZ14OtrbL2I6lfZfSsmVLsmX7grp1jVUKyhR8yuxmv5Hji0zkWzfEmBPfpo2pwxZCmIgk+Xjk5WCrk40T9zfdx9zcnJ9//pl+/YzB1ilTYNeQXUSERlAn9QnUk8cwa5Zs5ydEIibdNfHIfO/5HPY7THub9mzbso3vvvuO06czvhpsjbzmy5mVZyjXNitO6+ZC//5QqJCpwxZCmFC0Bl6VUrWB6YA58KvWeuI/lCsJHAFaaq03/NtnysDrx7n95Db5ZuejZJqS3PzxJmZmZhw/foZixaywsICTJyJZUno+YUFh9HVag6X/DRlsFSIB+tiB1w921yilzIHZQA3ADziulHLXWp97T7lJgOfHhSyio9/2foRHheN8zZk9V/bg6enJ9OlWXL0KO3fCyblHCTgbQMs+jljOOQorVkiCF0JEq0++FHBFa30NQCm1BmgInHunXH/AFSgZoxEKNp7fiNsFN4bmH8rMjjNp3LgxuXPXpGFDY7C1dL4nzGq8l1w1spBn9WAZbBVCvBKdPvkMwK03nvu9eO0VpVQGoDEw798+SCnVQynlrZTyDggI+NhYE6XHoY/pt60fhVMX5sySM5iZmTF9+vRXywhPmQI7B+9ER2rqpDwug61CiLdEJ8m/L1u825E/DfhGax35bx+ktV6gtS6htS6RKlWqaIaYuA3zGsZfz/6itVlrtm/fztixYzlzJtOrwdbwi9c4u+4sFTpkw3GNDLYKId4Wne4aPyDTG88zAv7vlCkBrFFG6zElUFcpFaG1douJIBOrgzcPMu/EPPoU6sOMfjMoWrQo3bv3x9kZ8uSB/n0iWFxiG445HCl/YgakSQOjR5s6bCFEHBKdJH8cyKWUygbcBloBb3X4aq2zvXyslFoKbJEE/988j3hOD48eZEmehYhdEfz111+4u7vzyy8WrwZbT8w6zINLD2jTPwUWM2WwVQjxdx/srtFaRwD9MGbNnAfWaa3PKqV6KaV6xXaAidXEgxM5f/88AzMOZOH8hfTt2xcHhxJMmGAsI1wydyD7x+4nb73s5Fo5WgZbhRDvJQuUxUHnA87jPN+Zxrkbc2HCBQICAjh37jxNmiTD2xvOn4cDfddyxfMK/RrfJvnaBXDypPTFC5EIyAJl8VyUjqK7R3fsLO3IdzUff/75JzNnzmTTpmTs2QOTJsGzU5e54HaBSp2zk3z1PPjyS0nwQoj3krVr4piFJxby+63f+anET4xuMRoXFxfKl29M/vxQvjx06RDBvMLbccrjRLmj02SwVQjxr6QlH4f4P/VnqNdQvsj6Bfvn7Qdg1qxZDBqkePoUFiyAQ5MP8ujqI+rWBfMTx4xlhJMlM3HkQoi4SpJ8HNJ/e3/CIsNoSlO2bNnCDz/8wLlzWVi1CkaMgNTm9zn440EKNMpF9v+NksFWIcQHSXdNHOF2wY2N5zfyfenv+bHbjxQpUoSuXQfg7Az58sGwbzSra3pgaWNJbfsD8FjubBVCfJgk+TjgyfMnxtIFaQrzcNtD/P39cXV15YcfLLhxAw4cgDNLvLl58CYNvy+M/djhMHCgDLYKIT5IknwcMNxrOP5P/Rmfdzyd+3Smd+/emJuXZvp06NkTCmV+zJw6XmSvno0iWyfIYKsQItokyZvYoVuHmOs9l34l+jHjuxmkSZOGMWMmUKOGkcsnTtRsabMFrTX1yz1A/eANK1fKYKsQIlokyZtQWGQY3T26kzFZRtKfT88ff/zBunXrWLw4OadOgasr3Nxyhivbr1B7dGkcJjeD2rWhdWtThy6EiCckyZvQ+P3jORdwjiVVltDfpT916tShaNFmdOwIDRtCzfLPmJN/BxnLZqTksVmgNcydK4OtQohokymUJrLn+h7G7h9L+8Lt2TxtM5GRkcyaNZvevRUWFjB7NuwYsJ2woDAaNLXEbNtWGDcOsmY1dehCiHhEWvImcOfpHdq4tiFPyjzUjaxLa7fWTJo0iYMHs+HlZST4p94XOLv2LF+MKEOqnzpA8eLGWvFCCPERJMl/ZhFREbTZ2IYnz5/g3sSdppWbUqhQIdq1+4rChaFcOejUKpR5hbaRpnAayvuvhwcPwNMTLORyCSE+jmSNz2zM3jHs9d3L0oZLWf7Lcvz8/Fi3bh3ffGPJkyfG0gW7h+0i6G4Qrb7PhXmv3jB0KDg7mzp0IUQ8JH3yn9HOqzsZf2A8nZ0743jTkVmzZjFgwACePCnLihUwbBjY/nWdPxb+QdkBJUk/ZTBkzw6jRpk6dCFEPCUt+c/k9pPbtN3YlgKpCzC88HDKlChD0aJF+f77SRQvbmznN2RgOEtKepAiZwqqWP4Oly8bW0DZ2po6fCFEPCVJ/jOIiIqglWsrQsJDWNNkDd1bduf58+esWbOGCROS4OsL+/bB4Qm/8ejaIzouqohlz0HQoQPUqGHq8IUQ8Zgk+c9g5J6RHLx5kJVNVuI635V9+/axdOlSgoJy88sv0L075LC5zaJfjlC8RzGyLhgBDg4wZYqpQxdCxHOS5GPZ1ktbmfT7JHoW70mmwEy0H9Oedu3a0aRJB0qXhtSp4cdxkayr7o59Onuq57wOC47C8uWQMqWpwxdCxHOS5GPRzcc36eDWAee0znxf6nvKlChD9uzZmT17Dn36KC5cAC8v8Jl/kHtn7tFqUXWsB9SGWrWgbVtThy+ESAAkyceSsMgwWm5oSXhkOOuaraNvt77cvXuXw4cPs25dUlasgDFjoFDaAOaN3U/B1gXJ4zYJoqJk6QIhRIyRJB9LhnsN54jfEdY1W8fONTtxc3NjypQpWFoWp39/qF4dhg+LYllld5IkS0LtamHQzcPYzi9bNlOHL4RIICTJx4LNFzYz9chU+pXsR+6I3LQf3J46derQpctASpUCR0djteDjM4/gd8SPxvNrYvdtYyhWDAYMMHX4QogERJJ8DLv+6Dod3TpSIn0JxpQfQ/ky5XF0dGTJkqX06mXG1avw228Qevkmu4ftJm+jvBQ6vgTu34dt22TpAiFEjJKMEoOeRzynxYYWAKxrto6hg4dy8eJFdu3axcaNqVm7FiZMgGK5g5hfbD0OWR1o2NUJVX8hfP210ZIXQogYJEk+Bg3ZNQRvf282tdzEsZ3HWLRoESNGjMDRsRoDB0KdOvD1oEhW1txAaGAo7dxbYN2uhrF8sGznJ4SIBZLkY8im85uYeWwmg8oMokiSIjj3cKZs2bJ89dXoV/Phly2DPSO8uLH/Bo1XNCaN6xy4eNFYYdLOztRVEEIkQJLkY8DtJ7fp5tGNEulL8EPlH6hWpRpKKVauXEXPnpbcvGksW/DXnrMcmXqEUv1LUTjlHZg0Cbp2hZo1TV0FIUQCJUn+P4rSUXRw60BoRCirmqxi7OixHD16lPXr1+PhkZWNG41ZkbkcAlhYczOZymWi5sD8ULoEFCgAM2aYugpCiARMkvx/NOXQFPZc38Ov9X/l+onrTJo0iR49epA5czPatIEGDaBvt+f8WnotVvZWNF/VCPO2DSA0FDZskBUmhRCxKlrrySulaiulLiqlriilhr3nfFul1OkXxyGlVJGYDzXu+ePOH3y751ua5muKSwYXOnToQIECBfjuu19o0QLSp4clSzTuXTbz8MpDmq1tRtI5k+D3343dQfLkMXUVhBAJ3Adb8kopc2A2UAPwA44rpdy11ufeKHYdqKy1fqSUqgMsAErHRsBxxbOwZ7RxbUNqu9TMd5lPt7bdCAwMZNcuL/r2tcXfHw4ehPOLDnF+43lqTqlJ1qc+8NNP0LMntG5t6ioIIRKB6HTXlAKuaK2vASil1gANgVdJXmt96I3yR4CMMRlkXDTIcxCXHlxid4fdbN+4HTc3NyZPnszOnQVxd4fp0yFV0HWWD9tN/ub5KdMkHRSrZ2zjN22aqcMXQiQS0UnyGYBbbzz3499b6V2B7e87oZTqAfQAyJw5czRDjHs2nd/Egj8W8E35b8htlZvG/RpTvnx5Spf+ii++gKZNoWPjJywovgGnPE40mFsbVa8GRETA+vVgbW3qKgghEonoJPn3LYeo31tQqS8wknyF953XWi/A6MqhRIkS7/2MuM7/qT/dPLpRPF1xxlQZQ6P6jQgPD+eXX5bSpIk5WbPC/DkRrG+wjoiQCFpubEmS8aPg6FEjwefMaeoqCCESkegkeT8g0xvPMwL+7xZSShUGfgXqaK0fxEx4cUuUjqLDJmO65MomK1m2ZBk7duxgxoyZjB+fk3v34PBhODLGk9tHb9N8Q3NSnj8Av/wC/ftDs2amroIQIpGJTpI/DuRSSmUDbgOtgDZvFlBKZQY2Au211pdiPMo4Yurhqey+vpuF9ReSJCgJgwYNomrVqlha9mHzZmO3PnOfP/Ge4025IeXIX9QainWGkiVh8mRThy+ESIQ+mOS11hFKqX6AJ2AOLNZan1VK9Xpxfh7wPeAEzFHGZhcRWusSsRf253fyzklG7B5B47yN6VykM9WrV0cpxYgRi6lf34waNaBVlbssKb+FrFWyUm1UBahc0dj8Y+1aSJLE1FUQQiRC0boZSmu9Ddj2zmvz3njcDegWs6HFHcHhwbTZ2IZUdqlYWH8hc+bMYe/evcyb9ytDhmTB1hbmTw9lg8s6bJxsaLqmKWbDhsKJE+DmJpuACCFMRu54jYZBnoO4eP8iu9rv4uHth3zzzTfUrVuXK1e6cPIkuLlpjg534/HNx3Ta3wn7fdtg1iwYNAgaNjR1+EKIREyS/AdsvrCZ+SfmM6TcEKpkqULFihWxtramQ4eFtGql6NULUl87ws7NF6k5tSaZUoZCrW5QtixMnGjq8IUQiZwk+X/h/9Sfru5dKZauGOOqjmPq1KkcPnyYuXNXMGhQevLmhcGt/Fhd3Yu8jfJSppczlCsHlpZGP7ylpamrIIRI5CTJ/4MoHUVHt44EhwezsslKLl+4zMiRI2ncuDGenm0ICICNK0PY0mEDyTImo8HiBqivvoJTp2DrVsiU6YPfIYQQsU2S/D+YdmQaXte8mO8ynxzJc1C2dlmSJUtGhQrzGDxYMXmy5tpUN57eeUqX37tgs9Md5s+HYcOgbl1Thy+EEIAk+feafWw2Q3YNoVHeRnQv1p2xY8dy4sQJpk/fwPDhqaleHcrqw3h5XKL29NpkSBMJNXpCmTIwdqypwxdCiFckyb8hMiqSr3d+zbSj03DJ7cLyxss5deoUY8eOpVWrNixb1hRra/ip3y08mnqRr2k+SvUpDlWrQlQUrFwJFvKTCiHiDslILwSFBdF2Y1vcL7rzZakvmVprKhHhEXTs2JFUqVKROvVM1qyB9f8LZm+/DThkcaDBogaoSZPgwAFjA9fs2U1dDSGEeIskeYxZNC6rXPjzrz+ZWWcm/Ur1A+C7Md9x5swZJkzYwrffpqBHd03YWjee3XtG18Ndsb54GkaNMtaGb9fOxLUQQoi/S/RJ/s+7f+Ky2oXA0EDcW7lTL3c9AI4ePcqkSZNo27YLs2fXI3duaJnpdw4svEzd2XVJl8seiraBjBlhzhxj+QIhhIhjEnWS33Z5Gy03tCR5kuQc6HwA57TOAPj7+9OuXTsyZMjA06dTuXcPloy5ycGee8jfPD8lepeALl3g+nXYtw8cHExaDyGE+CfR2uM1IZp9bDb1V9cnV4pcHO129FWCv3r1KhUqVODu3bu0abMad/fkjBv+jDPfb8AxmyMNfm2AWr8eli6Fb7+FCu9dOl8IIeKERJfkI6MiGbhjIP2296Nernrs77yfDMkyAHDmzBkqVKjA48ePWbJkDzNnlqfqF5r0RzYR/CCY5uubkyTwL2OP1tKl4bvvTFwbIYT4d4mquyYoLIg2rm3wuOTBgNIDmFJzCuZm5gAcOXKEunXrYmNjg7v7Afr1y4+1NQwqdRDvSVepN68eaQulMqZLRkQY0yVl2QIhRByXaJL8mzNoZtWZRd9SfV+d27VrF40aNSJdunQsWOBFx45Z8fODpWN8OTHiNwq2KkjxHsWNBcf27ze6anLkMF1lhBAimhJFkj955yQN1jQgMDQQj9Ye1M31etkBV1dXWrduTb58+Rg92pNmzdJibg7bNwTxRw9XUuRMgcsCF5S3N3z/PbRsCR06mLA2QggRfQm6Tz5KRzHtyDTKLCqD1poDnQ+8leAXL15MixYtKFmyJL1776Vly7SkTg0HvJ5zc4oroY9CjX54FQ5t2kD69DBvnkyXFELEGwm2Jf9X0F902tyJHVd2UD93fRY1WEQqu1Svzk+ZMoWvv/6aWrVqUbSoK7172/HFF7BwQgA7Wq3jweUHNFzSkDSF00DXrnDtGvz2m0yXFELEKwkyyW+7vI1Obp14GvaU2XVn07tEb17sPYvWmpEjRzJhwgSaNGmOpeUKJk60onNnGFDjLOtquGNpa0kHrw5krZIVNmyAxYuN6ZKVKpm2YkII8ZESVJIPjQhl6K6hzDw2k0KpC/Fb098okLrAq/NRUVH069ePuXPn0r59d65dm8vvv5sz7ocoSgR64dbmMBnLZqT5+uYky5AMbt2C7t2hVClj+QIhhIhnEkyS97nnQ2vX1vjc82FA6QFMrD4RawvrV+fDw8Pp2LEjq1evplu3oezZM5HbtxUr5gURvnoDR/bdoGTfktSaWgtzK3OIjIT27WW6pBAiXov3SV5rzezjs/l659ckt07OtjbbqJOrzltlgoODad68Odu2baNHj4msX/8N5uawccYtLoxZT8ijEBota0SR9kVev2nyZGPJgiVLIGfOz1wrIYSIGfE6yd97do+u7l3ZcmkLdXLWYUnDJaSxT/NWmQcPHtCwYUMOHTpEp07zWbKkBzmya6a2Po53P0+SZ0pO18NdSVskrfEGrWHqVBgxApo3h44dTVAzIYSIGfE2yXte8aSjW0cehT5ieu3p9C/V/9Xg6ku+vr7Url0bX19fmjZdy9KlzaleOZyuabZwbPRpctXLRePljbFxtDHeEBEB/fsb0ySbNjVuepLpkkKIeCzeJfnnEc8Zvns4vxz5hfyp8rOz/U4Kpyn8t3InT56kbt26hISEUrr0TjZsqET3pg8pemkdF/f/RZUxVag0shLK7EUSf/IEWrQAT08YOhR+/BHMEvRtBEKIRCDeJfnlp5fzy5Ff6FuyL5NrTMbG0uZvZTw9PWnWrBm2to5YW3tx6FABfux8CTZu5ImZos3WNuSqk+v1G27eBBcXOHcOFiwwZtQIIUQCEO+SfJeiXcjjlIeKWSq+9/zSpUvp3r07KVLk59697eTP5sikL3bgu+QoaYumpYVrCxyzOb5+g7c31K8PwcGwYwdUr/6ZaiKEELEv3vVHmCmz9yZ4rTXjx4+nc+fO2NpW5sG93xhQ0pf2j2bgu/YoxXsVp8vvXd5O8G5uxg1OSZLAoUOS4IUQCU68a8m/T0REBH379mPBgvlYmLcjb/gAGqZcRvjxx2Sqk5Pqk6qTptAbs25ezqAZMgRKlgR3d0iT5p+/QAgh4ql4n+SfPXtGs2at2bHDAwd60iWJM8mCt5IyXzqqr2lA9mrZ337DuzNoli0DW1vTBC+EELEsXif5gIAAKld24fz545SiA3VJR/LUz6k2vgkFWxV8PXPmpTdn0HzzDUyYIDNohBAJWrSSvFKqNjAdMAd+1VpPfOe8enG+LhAMdNJa/xHDsb7l4sWrlCtdk8DHfrSkBYWT5eeL0ZUo2ackFkneU62XM2jOn4eFC6Fbt9gMTwgh4oQPJnmllDkwG6gB+AHHlVLuWutzbxSrA+R6cZQG5r74M1asX3GADh3ro6LC6Kg60rRvE6oOKIiNRQRcvgBBQcbx9KnxZ2Cg0WoPCYHt22WAVQiRaESnJV8KuKK1vgaglFoDNATeTPINgWVaaw0cUUo5KKXSaa3vxHTA46qN4Ic9P2OHHWPMK9MRV5LPWgizPvDG7Nlh927Inz+mQxJCiDgrOkk+A3Drjed+/L2V/r4yGYC3krxSqgfQAyBz5swfGysARcvlIcdv6VlSsRJlnB0haXGwt//nI2lS408nJ7CI10MQQgjx0aKT9d63eIv+hDJorRcACwBKlCjxt/PRUW9sR+qNlUXDhBAiOqIztcQPyPTG84yA/yeUEUII8ZlFJ8kfB3IppbIppayAVoD7O2XcgQ7KUAZ4HBv98UIIIT7OB7trtNYRSql+gCfGFMrFWuuzSqleL87PA7ZhTJ+8gjGFsnPshSyEECK6ojUSqbXehpHI33xt3huPNdA3ZkMTQgjxX8ntnkIIkYBJkhdCiARMkrwQQiRgkuSFECIBU8aYqQm+WKkA4MYnvj0lcD8Gw4lvEnP9E3PdIXHXX+puyKK1ThXdN5osyf8XSilvrXUJU8dhKom5/om57pC46y91/7S6S3eNEEIkYJLkhRAiAYuvSX6BqQMwscRc/8Rcd0jc9Ze6f4J42ScvhBAieuJrS14IIUQ0SJIXQogELM4neaWUuVLqpFJqy3vOKaXUDKXUFaXUaaVUMVPEGFs+UPcqSqnHSqlTL47vTRFjbFFK+Sqlzryom/d7zifYax+Nuif0a++glNqglLqglDqvlCr7zvmEfO0/VPePvvbxYT+8AcB5INl7zn3WDcRN4N/qDnBAa+3yGeP53L7QWv/TzS8J/dr/W90hYV/76cAOrXWzF3tY2L5zPiFf+w/VHT7y2sfplrxSKiNQD/j1H4q82kBca30EcFBKpftsAcaiaNQ9sUuw1z4xU0olAyoBiwC01mFa68B3iiXIax/Nun+0OJ3kgWnAUCDqH87/0wbiCcE0/r3uAGWVUn8qpbYrpQp8nrA+Gw3sVEqdeLEB/LsS8rX/UN0h4V777EAAsORFV+WvSim7d8ok1GsfnbrDR177OJvklVIuwD2t9Yl/K/ae1+L9nNBo1v0PjDUsigAzAbfPEdtnVF5rXQzjf837KqUqvXM+QV77Fz5U94R87S2AYsBcrXVR4Bkw7J0yCfXaR6fuH33t42ySB8oDDZRSvsAaoKpSasU7ZRLqBuIfrLvW+onWOujF422ApVIq5WePNJZorf1f/HkP2ASUeqdIQr32H6x7Ar/2foCf1vroi+cbMBLfu2US4rX/YN0/5drH2SSvtR6utc6otc6KsXn4Hq11u3eKJcgNxKNTd6VUWqWUevG4FMa1fPDZg40FSik7pVTSl4+BmoDPO8US5LWPTt0T8rXXWt8Fbiml8rx4qRpw7p1iCfLaR6fun3Lt48PsmreoRLyB+Dt1bwb0VkpFACFAK51wbl9OA2x68XfZAliltd6RSK59dOqekK89QH9g5YvZJdeAzonk2sOH6/7R116WNRBCiAQsznbXCCGE+O8kyQshRAImSV4IIRIwSfJCCJGASZIXQogETJK8EEIkYJLkhRAiAfs/e/Kp01iWb88AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n=len(r)\n",
    "alpha=0.05\n",
    "eps=np.sqrt((1/(2*n))*np.log(2/alpha))\n",
    "\n",
    "F_n=lambda x: sum(r<x)/n\n",
    "L_n=lambda x: max(F_n(x)-eps,0)\n",
    "U_n=lambda x: min(F_n(x)+eps,1)\n",
    "Sde=lambda x: np.sqrt(((sum(r<x)/n)*(1-sum(r<x)/n))/n)\n",
    "NL=lambda x: F_n(x)+norm.ppf(alpha)*Sde(x)\n",
    "NU=lambda x: F_n(x)-norm.ppf(alpha)*Sde(x)\n",
    "\n",
    "xx=sorted(r)\n",
    "\n",
    "df=pd.DataFrame({\n",
    "    'x': xx,\n",
    "    'F_n': np.array(list(map(F_n,xx))),\n",
    "    'U_n': np.array(list(map(U_n,xx))),\n",
    "    'L_n': np.array(list(map(L_n,xx))),\n",
    "    'normal_L': np.array(list(map(NL,xx))),\n",
    "    'normal_U': np.array(list(map(NU,xx)))\n",
    "}\n",
    ")\n",
    "\n",
    "plt.plot('x','F_n', data=df, color='blue')\n",
    "plt.plot('x','U_n',data=df, color='green')\n",
    "plt.plot('x','L_n',data=df, color='red')\n",
    "plt.plot('x','normal_L',data=df,color='purple')\n",
    "plt.plot('x','normal_U',data=df,color='black')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1a4a6044fd0>"
      ]
     },
     "execution_count": 180,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABc10lEQVR4nO3dd1RURxvA4d8A0gULgl2wNxArNlRiV+w99hqJxmhiTDSNL9Fo7NFYY4uJvWPHAghiF0TFLqBYwIoiRcp8f1xcGyoaFIF5ztnD3rJ3Z1h9d5iZ+46QUqIoiqJkTXoZXQBFURTl/VFBXlEUJQtTQV5RFCULU0FeURQlC1NBXlEUJQszyKg3trKykra2thn19oqiKJnSsWPHbksp86X1/AwL8ra2thw9ejSj3l5RFCVTEkKEvc35qrtGURQlC1NBXlEUJQtTQV5RFCULU0FeURQlC1NBXlEUJQt7Y5AXQiwSQkQKIU694rgQQswQQlwUQgQJIaqkfzEVRVGUd5GWlvwSoNlrjjcHSqU8BgFz/nuxFEVRlPTwxnnyUsp9Qgjb15zSBlgqtZzFB4UQuYQQBaSUN9KrkIqS3d2Ouc3R60eJTYjN6KKkKyklSYlJJCYkaj8TE3XbL+571f7EhMSn5yRo+5+c82wq9ZfSqr+0+fyO5GSZ8rqU48kpL0nZTpZPj0spdfu1bZ5/rZTP7WvUoBEj+g5611/bW0mPm6EKAVef2Q5P2fdSkBdCDEJr7VO0aNF0eGtFyXqklITcD8Hvip/uceb2mQwqDJAAPH7m5+M37HvT8YSURxKQ/OGq8jG55BeeqYK8SGVfqiuRSCnnA/MBqlWrplYrURQgMTmRoIig54L6jWitjZTLOBd1itShV6Ve1Cpci1zGuf7z+z2Of0zkzUgirkcQcSPlcT1Ctx11L4rYmFjd423o6xuQw9AEgxwm6OcwRU/fBCFMkMIEqW9CkpEJiQYmJOYwJjE5BzLJAH300EcPg2d+Gkg9DBDkQA8DCYboYSjBWEqMkZgkS0xJxgyJaXIyxsnJ6KGHPvq6n/roI1INT0+9+rjQIpueHgiBFAL0QAgBQoBeynGhp3su9PWeOZayTy9lW3e+QOgJbNuWf6vf63+RHkE+HCjyzHZh4Ho6XFdRsqRkmcyBqwfYE7IHvyt+HAg/QPTjaACKWRbjE7tPqFu0LnWL1qV8vvLoibebBJecnExoaChBQUGcOXOGq1evEh4ernvcunXrpddYWlpibV0YC4uiWBS0AswAc5KSzEhMNCMhwYz4eHPi4syIidEesbHmKedpDwNyYJ6UiGlsDKaxjzAjhjzGMeQ2eoSFfgzmIgZjGYNR4iP0H8ciEh4jkpLSXjEBRpbGGOc2wczKFNO8JpjkMcHkyc8Xn+c2wcDYAKEnEPpacNXT19O2X7dPvP6LIbNJjyDvAQwVQqwEnIAo1R+vKM9LTE7EN8yXdWfWsf7Mem5E30AgcLBxoHel3tQtWpc6RepQxLLImy/2jAcPHnDy5EmCgoKee0RHR+vOyZMnD0WKFKFAgcKULFmDHDkKk5hYmEePCnPnTmGuXy9EeHhOoqKeXlcgMRWx5LeIIZ/pIwqYxJArRww5cz7CLGcMxsm3MUy4gn58DMQ8Ijk6Bvk48eUCxoFeoh6mVqaY5jPFLJ8ZplYFMMlrgqG5ITnMcmBoboihWcpzs9fv09NXs77f1huDvBBiBdAAsBJChAM/AzkApJRzgW1AC+AiEAP0fV+FVZTMJCEpAa9QL9YGr2Xj2Y3cirmFiYEJLUq1oEO5DjQv1TzN3S9SSkJCQggICCAoKIgTJ04QFBRESEiI7pxcuXLh4OBAr159yJOnEnFxDkRGlicszJyLF+HEieevmTePpGLRBzQtfpOCZU5iERuJuB3J41tRxN+LRSZLiEJ7PMPQ3PCZoG2OqZU1pvm0bVOrlED+zHMjS6Ms1zrOTERGLeRdrVo1qbJQKllNfGI8uy/vZu2ZtWw6u4l7cfcwNzTHtbSrFthLNsfM0CxN17pz5w579uxh165d7Nq1i7AwLfmgnp4epUuXxsHBgUqVKlGwoAOxsQ5cvFiEQ4cEx45BXJx2jXz5oEwZKF3oEbamkVglR2IUFUni9VvcPRtJ/IN43fuZ5zfHuqI1uYrneilQP3luamVKDpMc6f57U9JOCHFMSlktzeerIK8o/83jpMdsu7CNdWfW4XHOgwfxD7AwsqBNmTZ0KNeBJiWaYJLD5I3XiY+P58CBA3h6erJr1y6OHTuGlBILCwsaNmxIo0aNcHCoQVxcBQICTDh4EA4dgmvXtNcbGUGVKlCzagL2llfJdT+EqNPh3DoVScztGN37mOQxwbqiNfkq5sO6grX2vEI+TPOavq9fkZKO3jbIZ1g+eUXJ7O7F3mPesXnMPDyT6w+vk8ckDx3KdaBj+Y40tGuIkYHRa18vpSQ4OFgX1H18fIiJiUFfX5+aNWvi7u5Oo0aNMTKqzpYtBixeDAEB8GSssnhxqF8fnKomUtosnBzhoVzdF0r4/HCuPE7iqr6gQJUClGlbRhfMrStaY2ZjprpPshHVkleUtxR6P5TpB6ez4PgCHiU8onHxxnzp9CVNSjQhh/7ruzKioqLw9PRk27ZteHp6cv26NhGtdOnSNGnShMaNG1O3bgOCgizYuBE2bYLQUG1GXq1a4OICNaomYWt4nfuBIYTuDeWq/1US4xIReoL8lfNj94kdti62FK1bFKOcr/+iybQSE+HRI4iOfv0jMfHp3UnJyenzPD1e7+oKXbq8U9VVS15R3pPD1w4z5cAU1gavRU/o8an9p3xV8ysq5a/0ytdIKTlz5gzbtm1j69at+Pn5kZiYSK5cuWjSpIkusOfNWxRPT1i/Hvr2hbt3te6Xxo3h+++hQeUoIr1Oc3n3ZU5Nv0LAowQAbBxsqDq4KnYudhSrVwzjXMYf6tfxZnFxcP8+REVpP599REVpx1/3iI9/fjsm5mnwfjLokN508+D1nj5/cfu/PH+y7eDwfsqfChXkFeU1kmUyW85vYbL/ZHyv+GJpZMnIWiP5wukLClsUTvU1sbGxeHt7s3XrVrZu3UpoaCgA9vb2jBw5kpYtW1KzZk3u3TNgyxb44gvw9NTiVu7c0LIltG0L9arHcmVHMEFLg1g28AoA+crnw7GPI7YuttjWt8XU6j32oz8J0u/6iI9P5aIvMDAAY+NXP0xNIU8e7bmJCeTMCebmYGam/Xzdw8wMcuR4uwCcBanuGkVJRWxCLEtPLGXqwamcv3OeYpbFGF5zOP0r9yenUc6Xzr9y5YouqO/du5fY2FhMTExo1KgRLVq0oEWLFhQtWpRr17TW+tq14Oen/RVftCi0aaMF9lrVEwnZeZ6Ty05yfut5khOSsSprhX0Pe+w/tSe3Xe70qWBysjZie/EiXLqk/XzyuHFDC9KPH7/+GjlyaN9KuXJpD0vLp8/f9LCw0AK3gWpnvi3VXaMo/0FMQgyT/Scz8/BMbsfcplrBaqzssJIO5TtgoPf8f5f4+HjWr1/PnDlz8PX1BcDOzo7+/fvTsmVLGjRogLGxMWFhWlBfuxYOHNBeW6GC1g3Trh042CdzZV8YQcuC+KPdGeIfxGOe35waQ2vg0MOB/JXzv9tAaXIyhIU9H8CfPC5der6lnSOHNpJbogTUqPF88H7Vw9g4y7Z+sxIV5BUlxf4r++m7qS8X7l6gVelWjKw9Eueizi8F2LCwMObNm8fChQuJjIykePHijB8/nrZt21KmTBmEEFy8CH/8oQX2J3+wOjrC2LHQoQOUKSOJOBGhBfYVp3h47SGGOQ0p174cDj0csHWxffu7O+/e1eZUHjqEbn7l/ftPjxsba0G8VClo3hxKlnz6KFIE9PX/0+9P+TipIK9ke7EJsfyw9wemHZxGsVzF2NtrLy52Ls+dk5yczM6dO5k9ezZbt25FCIGrqyuff/45jRs3Rk9PjzNnYNw4LbA/ubu0enX4/XctsJcoAffD7nNy+Un2LjvJrdO30DPQo2TzkjSZ0oQyrcqQwzSNNxolJsLJk1owf/I4f147JgRUrAidOkG1atrdUCVLQoECWh+0kq2oIK9kaweuHqDPpj6cv3Met2puTGw8EXNDc93x27dvs2jRIubOnUtISAg2NjaMGTOGQYMGUbRoUR4+hBkz4K+/IDhYe03t2jB1KrRvD8WKQezdWILXBrO4bxBXfLUB1CJ1itBidgsqdKqQtsHT6GjYvRv8/bWAfvQoxKZkiLS2hpo1oU8f7We1atoApaKggrySTcUmxPKT109MPTiVwhaF2d1zNw2LNwS0aY8HDx5k9uzZrFmzhvj4eOrVq8f48eNp164dhoaGXLsG334L8+ZpswFr1YKZM7U+9kKFIDEukfNbzrNq+PMDqC5jXdI+gBoRAZs3w8aNWoCPjwdDQ6hcGQYN0gK6kxPY2qq+ceWVVJBXsp2D4Qfpu6kvZ2+f5bOqnzGp8STdjBkfHx9GjBhBQEAAOXPmZMCAAbi5uVGhQgVA64aZMgVWrNDGNTt2hK+/1sYqk5OSCfMJY9NPQZxZ+44DqBcuaEF940ZtlFZKLYi7uWlTcGrV0ibQK0oaqSCvZBtxiXH87PUzkw9MplDOQnj28KRxicaAlgzsm2++YfHixdja2jJ37ly6d++Oubk5UsLOnTB5stagNjODIUPgyy/Bzg7uXb6H5zdHOPUuA6jJyVrXy5PAfiZlBajKlcHdXZtXaW+vWurKO1NBXskWjlw7Qu+NvTlz+wwDKg9gStMpWBhZIKVk6dKljBw5kvv37/Ptt9/y008/YWpqSnw8LFmitdxPndLGLSdM0HpKcueGO+fvsLG3L0HLghBCpH0ANSEBvLxgwwbw8IDr17WZLfXray321q21znxFSQcqyCtZWnxiPO7e7kz0n0jBnAXZ0X0HTUs2BeDcuXMMHjwYb29vatWqxbx587C3t+fePW0wdcYM7b4ge3st2HfrpnWJ3wq+xfqhvpxaeQp9I32chjlR+5va5CzwmsHO+Hjtz4C1a7WENPfuaX8SNGumtdZbtNDu7FSUdKaCvJJlBd8KpuvarpyMPEk/x35MbToVS2NL4uLimDBhAuPHj8fExIS5c+cycOBAQI85c2DUKG0yS5MmWnBv3FjrLYkIimDf2H0Erw0mh2kOao2sRe2va2Nm/Yr88LGxWj/P2rXaAOqDB9pdoa1ba3MqmzTRbtVXlPdIBXkly5FSsjBgIcO2D8Pc0Jytn26lRakWAHh5eTF48GDOnz9Pt27dmDp1Kvnz5yc0FPr3h717taA+efLTHFI3jt9g36/7OLvxLIY5DXEe40zN4TVTn/r46BFs26YF9q1bte08ebQR2g4doGFDNXCqfFAqyCtZSlRcFIO2DGL16dU0Kt6IpW2XUiBnAW7dusXIkSNZunQpxYsXZ+fOnTRp0gQptWmQI0dqr583DwYO1Fru4YfC2ffrPi5svYBxLmPqu9fHaZgTJrlfaH3HxGiDpuvWwfbtWgs+Xz7o0UML7A0aaGkDFCUDqCCvZBkHww/SbV03rkZdZXzD8YyqMwqBYPHixYwcOZIHDx4wZswYfvjhB0xMTLhyBQYMgF27tAb2woXaeOdV/6v4/OLDpZ2XMMljgstYF2oMrYGx5QtpfCMiYNYsmD0b7tzRRmb79dNa7c7OKk2A8nGQUmbIo2rVqlJR0kNScpKc4DtBGvxiIItNKyb9r/hLKaUMDQ2VLi4uEpB169aVp06dklJKmZws5fz5UubMKaWZmZRz5mj7Hsc8lluHbJXuuMuJ+SZKv9/9ZNyDuJff8MwZKQcOlNLISEohpGzdWsq9e6VMSvqQ1VayKeCofItYq4K8kqndeHhDNl7aWOKO7LS6k7wXe08mJyfLJUuWyJw5c8qcOXPK+fPny6SUAHzlipRNm2r/8l1cpLx8WbvOzRM35awKs6Q77nLHiB0yPjr++TdKTpbS21tKV1ftxcbGUn72mZRnz37gGivZnQrySrax8+JOaT3JWhqPNZbzjs6TycnJMjIyUrZr104C0tnZWYaEhEgptRi9YIGUFhZSmppKOWuW1vBOTk6WB/84KH81+lVOzj9ZXtx58fk3SUiQcsUKKatW1f67WFlJ6e4uZWTkh6+wokgV5JVs4HHiYznKc5TEHVlhVgV5KkLrhtm8ebO0sbGRhoaGcuLEiTIxMVFKKeXVq1I2a6b9a69fX8pLl7TrPLz5UP7b/F/pjrtc7rpcRkdGP32TBw+knDpVyqJFtReWLi3lvHlSxsR84NoqyvNUkFeytEt3L8kaf9WQuCMHbx4sYx7HyIcPH8qBAwdKQDo4OMgTJ05IKbXW+99/S2lpqbXeZ8582m1+fut5Ocl6khxrPFYe+vOQTE5O1g48eCDlt99qLwIp69WTctMm1d+ufDRUkFeyrA1nNkiL8RbScrylXHN6jZRSSj8/P1m8eHEphJCjRo2ScXHaQGlUlJSffqr9C3d2lvJiSi9MQmyC3PbFNumOu5xtP1tGnIp4+ga7d0tZrJg2mNq5s5SHDn3gGirKm6kgr2Q5ycnJcrzveIk7ssZfNWTIvRAZHx8vR48eLfX09KStra308fHRnX/kiJQlSkippyflr79KmdJrIyNORsjZ9rOlO+5y+5fbZUJsgnbg4UMp3dykrltm//4MqKWipI0K8kqWEpcQJ3tt6CVxR3Zb203GJsTKkydPSkdHRwnIfv36yaioKCml1qMyebKUOXJIWaSIlL6+2jWSk5PloZmH5K9Gv8pJ1pPk+W3nn77B3r1S2tpqrfevvlJ97spHTwV5JcuIjI6UdRfVlbgj/+f9P5mYmCinTJkijYyMZL58+eTGjRt150ZESNm8ufYvum1bKe/c0fZHR0TL5a7LpTvuclmLZfLhzYfagYcPpRwyRHtByZJS+vllQA0V5e29bZBXd7wqH6XgW8G4LnflRvQNVnVcRZNCTWjVqhXbt2+ndevW/PXXX1hbWwOwZ4+WQeDePe0GVDc3LS1BmG8Ya7usJfZuLM1mNKPG0Braoh0+PtC3L4SGwvDh2sKspmlYgk9RMiEV5JWPzo6LO+iytgumOUzx6eODeZQ5NWrUICQkhFmzZuHm5oYQgoQEbV2N8eOhbFkt4aODg/bXqf+UA+z+dje5i+emx44e2DjYaMnCxozRcgiXKKEFe2fnjK6uorxXKsgrH5U/D//Jlzu+xN7ans3dNnPc+zg9e/bExMSEvXv34pwSlMPCtPzuBw5o2SP/+ENLzx4XFYdHPw/OrD9DufblaLO4DUYWRuDrq7XeL12CYcPgt9+0FyhKFveadcmeEkI0E0KcE0JcFEJ8l8pxSyHEZiHECSHEaSFE3/QvqpKVJSQlMGTrEL7Y/gWtSrdiX599LP5jMW3btqV06dIcPXpUF+DXrQNHR221phUrYMECLV5HnIzgr+p/cXbTWZpMaUKntZ0wMkiCESO0VZekBG/vp98IipIdvKnTHtAHLgHFAUPgBFD+hXPGAL+nPM8H3AUMX3ddNfCqPHEv9p4u/8woz1HyftR92bZtWwnInj17ypiUGS8xMVIOHqyNldao8fTOVSmlDPw7UI41GSsnF5gsQ/eFajsDAqQsU0Z7wdChUkZHv/zmipLJ8B4GXmsAF6WUlwGEECuBNkDws98VQE6hLUVvnhLkE//j94+SDVy8e5FWK1px6e4lFrVehHNOZ2rXqs25c+eYNm0aX375JUIIrl+HNm20Na9HjYJff9WW4kuMS2T7l9s5Pv84tg1s6bCiA+Y2Zlq/+zffaHnd9+yBTz7J6KoqSoZIS5AvBFx9ZjsccHrhnD8BD+A6kBPoIqVMfvFCQohBwCCAokWLvkt5lSzEJ9SH9qvbA7Cr5y7izsVR/ZPq6OnpsXPnTho2bAhAQAC0agX372vLo7Zurb3+fuh9VndczY1jN6jzbR0+GfsJevfvQptPteX2XF1h8WKwssqgGipKxktLn7xIZZ98YbspEAgUBByBP4UQFi+9SMr5UspqUspq+fLle8uiKlnJprObaPxPY6zNrDnU/xCHVh+iRYsWFC1alKNHj+oC/IYNULcu6OnB/v1PA/yFbReYV2Uedy/epeumrjSa0Ai9/b5aZ/3OnVq/u4eHCvBKtpeWIB8OFHlmuzBai/1ZfYH1KV1GF4EQoGz6FFHJajaf20ynNZ2oUqAKe7rt4cehP/Ltt9/SoUMH/P39sbOzQ0r4/Xdo3x7s7eHwYahUCZKTktn7416Wt1xOrmK5GHRsEGValNDmUn7yiTbf/eBBbQaNSK19oijZzJs67dG6dC4DdjwdeK3wwjlzAPeU5zbANcDqdddVA6/Z05ZzW2SOX3LI6vOryxNnT0hHR0cphJDjx4/XZYKMi5Oyd29tvLRr16eZBqIjo+XSRkulO+5yU/9N8nHMY20VEGdn7eTevbU7WRUlCyO9B16llIlCiKHATrSZNouklKeFEINTjs8FfgWWCCFOonXvfCulvJ2eX0ZK5rf9wnbar26Pg40DEytMpGHdhiQkJLBlyxZatGgBwO3b0K4d+PlpjfOfftIa5HfO32FZ82U8uPaA1gtbU7lfZa2Dvl8/ePwY/vlHu+1VUZTnpOlmKCnlNmDbC/vmPvP8OtAkfYumZCU7L+6k3ap2VMhXgfFlx9O2RVty587Nzp07KV26NADBwdpY6Y0bsHIldOmivfaq/1VWtF6B0BP08elD4UpW8MUX8OefUKWKdnKpUhlYO0X5eKk7XpX3bvfl3bRd1ZayVmX5rfRvdGzVESsrK7y8vHSzrHbuhM6dwcREu1/JKWX+VvC6YNZ3X49lUUu6b+9OnoRIqOkKJ05oNzmNHw9GRhlXOUX5yKXpjldFeVd7Q/bSekVrSuUpxdiSY+nUuhPW1tb4+PjoAvyff0KLFmBnpw2wOjlpY0UHph1gTac1FKxakP7+/clzeCdUrQrXrsHWrTB1qgrwivIGKsgr741PqA+uy10pkacEv9r9Ste2XSlUqBA+Pj4ULlyYxEQYMkTreXF11frhixbVZtDsGL4Dz688Kde+HD09u2M6YwJ8+qkW5E+c0L4VFEV5I9Vdo7wXvmG+tFjeArvcdvxc5Ge6te+GnZ0de/bsIX/+/Dx8CB07gqendmPq+PGgrw8JMQms776esxvPUvOrmjT5X11E356wZo02yDpnjnarq6IoaaKCvJLu9l/ZT/NlzSlqWZSfCv1Ez049KVWqFLt378ba2proaGjZEvz9YeFCLXYDPLr1iBWtVnDt8DWa/dEMp05F4RMXLZfBpEnw9ddq7ruivCUV5JV0deDqAZota0Yhi0J8b/M9vbv0pmzZsuzevRsrKysePdK6Zvz9YflybbAVUqZItljGw2sP6byuM+Vs46BGDW0lkI0bn97qqijKW1F98kq6ORR+iKb/NqWAeQFGW4+mf/f+VKhQgb1792JlZUVMjJaDxtcX/v33aYC/6n+VhbUXEh8VT2+v3pQT57RcBkI8n8tAUZS3poK8ki6OXDtCk3+bYG1mzag8oxjUYxCVKlVi9+7d5MmTh9hYLYukjw8sXQpdu2qvC14XzN+f/I1JHhP6+/ejsM8yLZdBxYpPcxkoivLOVHeN8p/djrlNy+UtsTK1YmSukbj1caNatWrs2LEDS0tL4uKgbVst4+/ff0P37trrDkw7gOfXnhSuWZhua9pj+sNXsGSJdhfU4sXapHlFUf4T1ZJX/rMvd3zJ/bj7DDYazNB+Q3FycmLnzp26AN+uHezaBYsWQc+e2mv2jdunmyLZa1VLTLu10QL8zz9ryz2pAK8o6UK15JX/xOOcB8tPLqdDQge++/w7nJ2d2bJlC+bm5sTHQ4cOsGOHtkRfnz7aaw5OP4jXD1449HCgzajS6DWoo93gtGLF034cRVHShQryyju7H3efwVsGUzqxNJsnbaZevXps2bIFMzMzHj+GTp1g2zaYN09bbBvg6Lyj7Byxk3IdytGmmwl6dWtrrXYfn6e5DBRFSTequ0Z5ZyM9RxLxIAKDzQZYWFiwatUqzMzMSEjQZs5s3gyzZ8OgQdr5J/45wVa3rZRqUYoO7ZLRa+0KtrZPcxkoipLuVEteeSe7Lu1iYcBC6ofXx+eED6tXr8ba2pqEBK3HZdMmLSeNm5t2/uk1p9nUZxN2LnZ0/lQf/d7doU4d2LIFcubM2MooShamgrzy1h7GP2Tg5oHYPbbjwL8H6NixI506dSIxUZs5s369tvrekCHa+ee3nGf9p+spXKswXfsaY9D7U6hVS0syZm6esZVRlCxOddcob230ntGE3Q3DaIsRFhYWzJo1i8REbc2ONWu05JDDhmnnXt59mdUdVpPfMT+fDrbAsG93rWtm2zYV4BXlA1AteeWt+Ib5MuvILGpfrY1/kL+um6ZvX1i1SksxM2KEdm6Ybxgr26wkb5m89BieF+M+XaBaNdi+XXXRKMoHooK8kmYxCTH08+hH4djCHF1+VNdNs3atNsX9xx9h5Ejt3GuHr7G85XIsiljQc1QBTPp21lZx2rEDLCwytB6Kkp2o7holzX72+pmLty5its1M101z+7bW916tmrYeK8DNEzf5t9m/mFqZ0uv7Ipj37wIODtryT5aWGVsJRclmVEteSZPD1w4z9eBUaoTV4PDJw7pumh49tESRu3eDgQHcOnOLfxr/g6GZIb1+tsViYBeoUEFLHJ8rV0ZXQ1GyHdWSV94oPjGefpv6kS86H4ErA3XdNFu2wLJl8P33YG8Pdy/dZWnDpQg9Qa9fS5J7cFcoW1bLaZAnT0ZXQ1GyJRXklTca5zuO0zdPY7HdQtdNc/8+fPaZFtxHj4aoK1EsbbiUpMdJ9BpXhrxunaFUKa2JnzdvRldBUbIt1V2jvFbgzUDG+42nckhlAk4F6LppBg6EiAjw8AC95ESWNV9G3P04ek+sgPUXXaBkSS3tpJVVRldBUbI11ZJXXikhKYF+m/phGWXJ6dWndd00u3ZpCcdGjtTW1fb51Ydbwbfo+H0ZCozoCnZ2WoDPly+jq6Ao2Z4K8sorTfafTMC1AHLtzKXrpomOhoEDoUwZLSvwzRM38Z/oT6XmBSj5v55QtCjs3QvW1hldfEVRUN01yiucuXUGdx937C/Zc/LUSV03zbBhcOWKtoSfUY5kPPp7YGxpSJMD/4NChbQAb2OT0cVXFCWFaskrL0lKTqKfRz9M7ppwbu05XTeNry/MnAlffKHlFjs4/SA3jt2geYnzmMbe0bKSFSiQ0cVXFOUZKsgrL5l5eCYHww6S1zOvrpsmNlbLCW9nB7/9pk2X9PrJi9JVc1Lh8CJwd9emSyqK8lFR3TXKc0Lvh/L93u8pc64M506f03XTjBoFFy5oMyJNTSXrBm1BT1/QMvRPRNWqT/MZKIryUVEteUVHSonbVjdkpCRkY4ium+bIEZgyRVv8o2FDCFwcSMjeEBqXv4ZF1FVt8VYD1V5QlI9RmoK8EKKZEOKcEOKiEOK7V5zTQAgRKIQ4LYTwSd9iKh/CylMr2XF+B9a7rHXdNPHx0Lev1tU+cSI8vPEQz689KVbRnKqH52i3uzo4ZHTRFUV5hTc2v4QQ+sAsoDEQDhwRQnhIKYOfOScXMBtoJqW8IoRQ8+cymTsxd/hyx5cUO1+MsOAwli1bhrW1NT//DKdPaws4WVrC6v7bSYhNoFXkUkTFCjBmTEYXXVGU10jL39g1gItSyssAQoiVQBsg+JlzPgXWSymvAEgpI9O7oMr7NXLXSO7evEu0RzTNmjWjW7dunDihDbL27AktW8KZDWc4s+4Mn1SPIu+xs7DlIBgaZnTRFUV5jbR01xQCrj6zHZ6y71mlgdxCCG8hxDEhRK/ULiSEGCSEOCqEOHrr1q13K7GS7vZc3sOSgCXY+tmiJ/SYO3cuSUmCfv20vGLTpkHc/Ti2DdmGTQkzah/5QxtorV49o4uuKMobpKUlL1LZJ1O5TlWgIWACHBBCHJRSnn/uRVLOB+YDVKtW7cVrKBkgNiGWz7Z8hk2oDZcOX2LatGkUK1aMCRPg+HFYu1bLL+Yx0JNHEY/oln8z+qVLalMmFUX56KUlyIcDRZ7ZLgxcT+Wc21LKR8AjIcQ+oBJwHuWj9ovPL1y6dolcW3JRvXp1vvjiC86e1WJ4x47QoQOEeIUQsCCAWlViKRhwDPbtAxOTjC66oihpkJbumiNAKSGEnRDCEOgKeLxwzibAWQhhIIQwBZyAM+lbVCW9nbh5gkn+kyh5pCTRUdEsWLAAIfQZMADMzODPPyEhNoHNAzeTu5AJLsenwtChULduRhddUZQ0emNLXkqZKIQYCuwE9IFFUsrTQojBKcfnSinPCCF2AEFAMrBASnnqfRZc+W+SkpMYuHkgOcNzctHrImPGjMHBwYG5c2H/fm3NVhsb2PWtN/cu3aNX4T3ksC2kjcQqipJpCCkzpmu8WrVq8ujRoxny3gpMPzidEZtHYL3UGksTS4KCgrhzx5jy5bXx1F274GbADf6q8ReOFRJoHTRWu921YcOMLrqiZGtCiGNSymppPV/dppgNhd0P44e9P1D8RHEuh19mtfdqjI2N+eILePwY5s6F5MQkPPp7YJbLkMYnf4cBA1SAV5RMSAX5bEZKyefbPifpWhKh20MZMGAA9evXZ8MG2LABJkzQFnXym3CAm4E36Vz0ACbGuWHy5IwuuqIo70AF+Wxm1elVbDu7jUK7C5FkncTEiROJitLGUytVgq++gjvn7+Dt7k25clDuzE7YvFm73VVRlExHBfls5G7sXb7c8SWFzxQm/Hw4a9asIXfu3Hz+Ody8CRs3goGBZMtnWzAw1KP5+SnQvTu4umZ00RVFeUcqyGcjIz1Hcjv8NoZbDGnTpg0dOnRg/36YMweGD9cGXAMWBRLqHYpr0RPkNDGEP/7I6GIrivIfqCCfTewN2cvigMUU8y3G3Rx3mTVrFo8fCwYO1JZl/fVXiI6IxnOkJ0Xt9KgSshHWrNZud1UUJdNSQT4beJK6IN/FfIQdD2P27NkUKlSIX36BM2dg61YwN4d1g3aSEP0Y10fzEO3aare8KoqSqalFQ7KBX/f9ysUrF4nbGkedOnX47LPPOHMGxo2Drl2hRQu4sP0Cp1acom7hEPIZP9Rud1UUJdNTQT6LOxlxkkn+k7A7ZEd8TDzz588H9Bg0SEtdMH06PH70mK1uW7EqkIO6IUvh99+hYMGMLrqiKOlAdddkYVJKhmwbgnGIMSG+Ibi7u1O+fHnmzQM/P23VPhsb2Pm1F1FhUfTJuRaDOjW1df4URckSVJDPwpadXIbvBV/ybM9D+fLl+e6777h+HUaNAhcX6NMHrh+7zqHph6hS8j7Fws7B/EDQU3/gKUpWoYJ8FhUVF8U3u76hwPEC3Lh5g83rNmNkZMSwYRAfD/PmgUxKZvPAzZhZGtD44lz4eQyUL5/RRVcUJR2pIJ9FuXu7czPkJgZ7Dejbty+1a9dm0yZYt05LJFmqFPhPOcjNgJt0zOeFcTk7GD06o4utKEo6U0E+CzoZcZIZh2ZQeH9hHpo+ZMKECTx4AEOGgL29tnLfvZB7eP/kTeniCZS/7APrfcHIKKOLrihKOlNBPot5MthqesmU8MBwZsyYgbW1NUOGwPXrWkvewECy7fNtCJJpETIL8dlnaiEQRcmi1AhbFrPs5DJ8L/piuMsQe3t73Nzc8PfXUhd88QU4OcGplae4uOMin+QOwDK/iTZlUlGULEm15LOQJ4OtBYMKcv3mdTas2kBysgGDBkHhwjB2LMTejWXHlzsoWESP6lc3wLq1KsOkomRhKshnIe7e7ty8chPDXYZ0796devXqMWECnD6tZQvOmRM29fck9m4sPQ0Wote2DbRvn9HFVhTlPVLdNVnEyYiTzDw8k6IHimKYw5CJEydy9aqWeKxtWy1bcIhXCIGLAqld+Ar5De+p1AWKkg2olnwWoBtsDTHlypErTJo0iYIFC9KpEyQnw7RpkBiXyJbPtpA7nwH1w5bCrOlQqFBGF11RlPdMteSzgGUnl+F7yRfj3caULVuWYcOGsWsXrF0L338Ptrawb9w+7l64S8u4teSoXR0GD87oYiuK8gGolnwmFxUXxUjPkRQ6VYhr4ddYvms5UhoydKi2VuvIkRB5KpL9E/bjYPeAEuGnYX6ASl2gKNmECvKZnLu3OxHXIjDaaUTHjh1p1KgREybA+fOwbRsYGUo2D9yMkakeTULmwo/fQYUKGV1sRVE+ENWcy8SCIoKYeXgmdoft0NfTZ8qUKVy58nSwtXlzOPbXMcIPhtPE0AuzMkVgzJiMLraiKB+QaslnUrrB1iumhPiHMG7cOIoWLUqnTiClNtj66NYj9ozeQ7HCiVQK3wPrfMDYOKOLrijKB6SCfCa17OQy/C77kX93fmxK2vD111/j6akNtv76qzbYuqn/bh4/jKfF/b8QAwdCvXoZXWxFUT4wFeQzoSeDrUXOFuFq2FW2bt0KGPHFF08HW6/6XyVwUSC1CoRgnZisUhcoSjalgnwm9LP3z0TciMB0uymtW7emRYsWusHW7dvB0CCZrZ9vJWdufRrcWAGL50Hu3BldbEVRMoAaeM1kgiKC+PPwn5Q8WpKkxCSmTZv23GBrs2ZwZPYRIk5E0DRxC4a1q0GvXhldbEVRMohqyWciyTKZIduGYHbdjIv7LvLTTz9RvHhxOnZ8Otj68MZDvH70oniRx5QPPwKzjqk58YqSjakgn4ksPL4QvxA/Cu0uRO5iufn222/x9NRyxI8dqw22ru+xi8TYBFo8nI8Y8jk4OmZ0sRVFyUBpauIJIZoJIc4JIS4KIb57zXnVhRBJQoiO6VdEBeDGwxt8s+sbSl4uybVL15g+fTr6+qbPDbaGeodyctlJaltfIK+V0PpwFEXJ1t7YkhdC6AOzgMZAOHBECOEhpQxO5bzfgZ3vo6DZ3bAdw4i9F8tNj5s0bdqUNm3aPDfYaqCXxLYh27DMq4/ztZWw5C/IlSuji60oSgZLS0u+BnBRSnlZSvkYWAm0SeW8L4B1QGQ6lk8BPM55sDZ4LaWPlCbhcQJ//vknV68Kxo6Fdu20wdaD0w9yK/gWzR9v0hKQ9eyZ0cVWFOUjkJY++ULA1We2wwGnZ08QQhQC2gGfANVfdSEhxCBgEEDRokXftqzZ0oP4BwzZNoRit4pxyusUv/zyCyVLlnxusPVB+AN8/udDadt4ylw5rgZbFUXRSUskEKnsky9sTwe+lVImve5CUsr5UspqUspq+fLlS2MRs7fv93xP+O1wEjYnUKZMGUaNGqUbbP3+eyhWDHaO2IlMTKJZ2DwYMkQNtiqKopOWlnw4UOSZ7cLA9RfOqQasFEIAWAEthBCJUsqN6VHI7OrA1QPMOjKLqhercuzKMby8vHjxztZLnpcIXhuMS+Hz5H5sAL/8ktHFVhTlI5KWIH8EKCWEsAOuAV2BT589QUpp9+S5EGIJsEUF+P/mcdJjBm4eiHWMNSfWn6B37940aNCA8eOfDrbqk8i2odvIY2NA7fBVsGShGmxVFOU5b+yukVImAkPRZs2cAVZLKU8LIQYLIdTyQu/JpP2TOB1xmrx78mJhYcHkyZO5fPn5O1v9J/lz98JdWsSuw6BOTTXYqijKS9J0M5SUchuw7YV9c19xbp//Xqzs7fyd8/y671eqRVbj6LGjLFy4kLx5rfj0U9DXhxkz4F7IPXzH+VK+RBwlQk7ArONqsFVRlJeoO14/MskymUGbB2EYZ8jFlRdxdnamb9++LFsGu3bBzJlQpAisaL0DISRNL82GL4ZApUoZXXQlDRISEggPDycuLi6ji6J85IyNjSlcuDA5cuT4T9dRQf4jsyhgET5hPtQMqMmxR8eYO3cud+4IRoyAmjXBzQ3ObT7H+c3naVT0LBZxxmqwNRMJDw8nZ86c2NrakjJRQVFeIqXkzp07hIeHY2dn9+YXvIb6+/4jcjP6Jt/s+gb7R/Yc3HaQUaNGUb58eb76Cu7fh7/+guT4BHYM20G+ggbUvLIaJk1Sg62ZSFxcHHnz5lUBXnktIQR58+ZNl7/4VJD/iHy540sexT7i4fqHlChRgu+//x5PT/jnH/j2W6hYEfwm+HE/9D4totegX6eWGmzNhFSAV9Iivf6dqO6aj8SW81tYfXo1LmEueF3yYufOnUhpwuDBULo0/PAD3Aq+xf7f92NfKhbbSye1wVYVMBRFeQ3Vkv8IPIx/iNtWN0oml2T/8v107dqVJk2a4O4OISEwfz4Y5khmU79NGJro0fTCnzB0qBpsVd6Jvr4+jo6OuseECRNeee7GjRsJDn6ai/Cnn35i9+7d/7kM9+/fZ/bs2W/9Ond3dyZPnpzqsaVLl1KxYkUqVKhA+fLldef16dMHOzs7KlWqROnSpenVqxfXrl3Tvc7W1hZ7e3vd78Pf3//dKvWRUi35j8D3e78nPCqcqrurYmJiwrRp0wgIgKlTYcAAqF8fDk4/zLVD12hfIhAzIzP43/8yuthKJmViYkJgYGCazt24cSOurq6UL18egF/SaZD/SZD//PPP0+V627dvZ/r06Xh6elKwYEHi4uL4559/dMcnTZpEx44dkVIyffp0XFxcOHXqFIaGhgB4eXlhZWWVLmX52Kggn8EOhR/iz8N/0ii6Ebv9dzN79mysrPLTsiVYWcHEiXDv8j32fr+X0hVzUPHURli6VA22ZgHDh0MaY22aOTrC9Onv9trvvvsODw8PDAwMaNKkCe3bt8fDwwMfHx/Gjh3LunXr+PXXX3F1daVjx47Y2try6aef4uXlRUJCAvPnz2f06NFcvHiRb775hsGDBxMdHU2bNm24d+8eCQkJjB07ljZt2vDdd99x6dIlHB0dady4MZMmTWLSpEmsXr2a+Ph42rVrx/9SGjLjxo1j6dKlFClShHz58lG1atWXyj5+/HgmT55MwYIFAW364cCBA186TwjBiBEj2LBhA9u3b6dNm9QS6mYtKshnoISkBAZuHkh+/fwELg3EycmJzz77jOnT4fhxWLUKcuWS/NNxM3r60DLkT0TjxtCjR0YXXcnEYmNjcXwmid3o0aNp3LgxGzZs4OzZswghuH//Prly5aJ169a6oJ6aIkWKcODAAUaMGEGfPn3Yv38/cXFxVKhQgcGDB2NsbMyGDRuwsLDg9u3b1KxZk9atWzNhwgROnTql+4vC09OTCxcucPjwYaSUtG7dmn379mFmZsbKlSsJCAggMTGRKlWqpBrkT506ler+V6lSpQpnz57VBXkXFxf09fUxMjLi0KFDaf9lZgIqyGegCX4TOBl5ksanG7P37l7m7ZpHWJgeP/4Irq7QqRMELAwgZG8Irg5hWFy4D3PnqsHWLOJdW9z/VWrdNYmJiRgbGzNgwABatmyJq6trmq7VunVrAOzt7YmOjiZnzpzkzJkTY2Nj7t+/j5mZGWPGjGHfvn3o6elx7do1IiIiXrqOp6cnnp6eVK5cGYDo6GguXLjAw4cPadeuHaamps+9338l5fOJdLNyd40aeM0gPqE+uPu400i/EbvW7GL48OE4OFTCzU3LTjB7Njy8/gDPrz2xrWBKlaAlWj988eIZXXQlCzIwMODw4cN06NCBjRs30qxZszS9zsjICAA9PT3d8yfbiYmJLFu2jFu3bnHs2DECAwOxsbFJde63lJLRo0cTGBhIYGAgFy9epH///kDaphJWqFCBY8eOpanMAAEBAZQrVy7N52dmKshngIjoCLqu60oJixJcW3GNIkWK4O7uzvLlsHMnjBsHhQtLtrptJSkhiVYRCxCVHWHEiIwuupJFRUdHExUVRYsWLZg+fbqupZ8zZ04ePnz4zteNiorC2tqaHDly4OXlRVhYWKrXbdq0KYsWLSI6OhqAa9euERkZSb169diwYQOxsbE8fPiQzZs3p/o+o0ePZtSoUdy8eROA+Ph4ZsyY8dJ5UkpmzJjBjRs30vxFltmp7poPLCk5ie7ru3M/7j4dozvyZ/CfbNq0ifh4c4YPBycnbd2P06tOc37zeZrUjCLP4Yuw4zAYqI9L+e9e7JNv1qwZX375JW3atCEuLg4pJdOmTQOga9euDBw4kBkzZrB27dq3fq/u3bvTqlUrqlWrhqOjI2XLlgUgb9681KlTh4oVK9K8eXMmTZrEmTNnqFWrFgDm5ub8+++/VKlShS5duuDo6EixYsVwdnZO9X1atGhBREQEjRo1QkqJEIJ+/frpjn/zzTf8+uuvxMTEULNmTby8vHQza7I8KWWGPKpWrSqzI3cvd4k7csySMdLAwEB27txZSill795SGhhIGRQkZXRktJxoNVH+VW6KTEJI+fXXGVtoJd0EBwdndBGUTCS1fy/AUfkWsVY1DT+g3Zd38z+f/9GlRBdW/LSCQoUKMW/ePHbvhr//hjFjwN4e1nffSVxUHK3Nt6BnW0zNiVcU5Z2pIP+B3Hh4g+7ru1PWqizJW5K5cuUK+/btw9AwF599BqVKwY8/wvkt5zm5/CT16ydj7XMYduwAM7OMLr6iKJmUGnj9ABKTE+m2rhvRj6PpkdyDNavW4O7uTu3atfnlF7h8GebNA+Lj2DJ4C9alLHHe/zt07w5Nm2Z08RVFycRUS/4D+NnrZ3zCfPi98u/88ukv1K9fP2W6GEyeDP36gYsLbP5sF9E3oulS4QD6luaQMvilKIryrlSQf892XNzBb36/0adiH1a6r8TIyIh///2XuDh9evSAvHm1lPAhXiEcn3+cWg1NKLRnh9ZJny9fRhdfUZRMTgX59+hq1FV6rO+BvbU9OX1zEhAQwMaNGylcuDB9+0JwsNblntM4gX8HbCaPrQUuh36FRo1UnnhFUdKF6pN/TxKSEui6rivxSfEMzTWUmX/MZMiQIbRp04bFi2HJEi1HfJMmsPfHvdy7fI9WhY+RIylOpS5Q3qsnqYYrVKhApUqVmDp1KsnJyQB4e3s/l9Lghx9+oGnTprRr146NGzfq9pcpU4axY8fqtjt06MD69eufe58XrwVa2t93mW//XzybStje3p5Nmza99/ds0KABR48efe/vkxYqyL8n3+/9Hv+r/kysOZEfvvwBe3t7Jk+ezKlT2s1OLi7w888QfiicQ9MPUbVJHmz9loG7O5QokdHFV7KwJ7lrTp8+za5du9i2bZsu4+Ozxo0bx/79+9m4cSO1a9fW5Vm/c+cO5ubmHDhwQHfugQMHqF279gerw6skJSWlut/Ly4vAwEDWrl3LsGHDPnCpMpbqrnkPNp/bzCT/SQyuMpj1v60nOjqalStXkphoTKdOYGEBy5eDTEzEo58HOfOb0fjEZC1P7FdfZXTxlQ9k+I7hBN4MTNdrOuZ3ZHqz6Wk+39ramvnz51O9enXc3d11+6dMmcK2bdvYuXMnJiYm1KlTh1GjRgHg7++Pq6sr27dvR0pJaGgoJiYm5M+f/63KamtrS5cuXfDy8gJg+fLllCxZkj59+mBsbMzp06eJiIhg6tSpuLq6kpSUxHfffYe3tzfx8fEMGTKEzz77DG9vb/73v/9RoEABAgMDn1vk5EUPHjwgd+7cuu2pU6eyaNEiAAYMGMDw4cMJDQ3F1dWVU6dOATB58mSio6Nxd3enQYMGODk54eXlxf3791m4cCHOzs7ExsbSt29fgoODKVeuHLGxsW/1u3ifVJBPZ2H3w+i9sTeV81em8OnCzN09l/nz51OuXHl69IDz52H3bsifH/aM8eFW8C0+bX4Xo53hsGWdSl2gfHDFixcnOTmZyMhIAPbv38+5c+c4duwY5ubmAFStWpVTp07x+PFj/P39qV+/PpcvX+bMmTMEBARQp06dd3pvCwsLDh8+zNKlSxk+fDhbtmwBIDQ0FB8fHy5duoSLiwsXL15k6dKlWFpacuTIEeLj46lTpw5NmjQB4PDhw5w6dQo7O7tU38fFxQUpJZcvX2b16tUAHDt2jMWLF3Po0CGklDg5OVG/fv3nvgRSk5iYyOHDh3V/Ae3evZs5c+ZgampKUFAQQUFBVKlS5Z1+H++Diijp6HHSYzqv7UySTOKH4j/QpUUXOnbsyIABA/jrL631/uuvWlfNOY9z+I33o3LLgpTa6q4lH6tWLaOroHxAb9Pift/kM6l3S5Ysyb179/D09NTlkTcyMqJChQocP36cgwcPMmrUKC5fvoy/vz8BAQGpdtW8Knvks/u7deum+znimQR8nTt3Rk9Pj1KlSlG8eHHOnj2Lp6cnQUFBuj79qKgoLly4gKGhITVq1HhlgIenqYQvXbpEw4YNadCgAX5+frRr1w6zlJsN27dvj6+v7xvTGbdv3x7QvvhCQ0MB2Ldvn64byMHBAQcHh9de40NSffLp6Ntd33L42mH+/ORPvhn8DQULFmT+/PmcOCEYNky7r2nMGLhz4Q4bem6gQJX8tLgwHYoVg3RaVk1R3tbly5fR19fH2toaABsbG7Zt28aIESN0XSkAtWvXZt++fTx8+JDcuXNTs2ZN/P398ff3T7UlnzdvXu7du/fcvrt37z6Xt/3ZgP+q50+2pZTMnDlTl444JCRE15I3S+Nd4SVKlMDGxobg4OCXcso/YWBgoBuIBl5KjfwkpbK+vj6JiYmvLPPHQgX5dOJxzoPph6YzrMYwdvyxg9DQUJYvX46+fm46ddKW8vvnH0iMecyqdqvQy6FH59pXMTgfDHPmQMqfxYryId26dYvBgwczdOjQ54JU6dKlWb9+PT169NClHa5Tpw7z5s2jUsoC8g4ODhw8eJArV65QoUKFl65dqlQprl+/zpkzZwAICwvjxIkTz2XAXLVqle7nkwyUAGvWrCE5OZlLly5x+fJlypQpQ9OmTZkzZw4JCQkAnD9/nkePHr1VfSMjIwkJCaFYsWLUq1ePjRs3EhMTw6NHj9iwYQPOzs7Y2NgQGRnJnTt3iI+P13UhvU69evVYtmwZoK1SFRQU9Fblep9Ud006uP7wOv029aNKgSpUiqhE/+X9+eWXX6hduw5dukBICHh7g5WVZF03D26fuU33ceXINaYL9OkDzZtndBWUbORJquGEhAQMDAzo2bMnX6Uy4F+9enUWL15M69at8fLyonbt2ly+fJnRo0cDWovX2tqaIkWKoKf3cnvxyY1/ffv2JS4ujhw5crBgwQIsLS1158THx+Pk5ERycjIrVqzQ7S9Tpgz169cnIiKCuXPn6latCg0NpUqVKkgpyZcv33PTOl/nyfJ+CQkJTJgwARsbG2xsbOjTpw81atQAtIHXJytT/fTTTzg5OWFnZ6dLj/w6bm5u9O3bFwcHBxwdHXXX/Ci8TcrK9HxklVTDSclJstHSRtJ0nKncfnC7NDMzk/Xq1ZOJiYnyzz+lBCl//10713+qv3THXfqO3iplvnxSli8vZXR0xlZA+aBUquGnihUrJm/duvXS/t69e8s1a9ZkQIk+PirV8Edg+sHp7L68m9lNZ/P9599jZGTEsmXLCAjQZ8QIaNkSRo6EUO9Qdn2zi7Jty1DHdwI8egQ+PirDpKIo71WagrwQohnwB6APLJBSTnjheHfg25TNaMBNSnkiPQv6MQq8GcjoPaNpV7YdVz2ucvz4cTZs2ICZWWE6dYICBbQUNNHXH7C2y1rylMxD2+JBiI2+Wgd9NlljUlFS82RmyouWLFnyQcuR1b0xyAsh9IFZQGMgHDgihPCQUj57x0EIUF9KeU8I0RyYDzi9jwJ/LGISYui2rhtWpla4FXSj2afN6Nu3L23atKVdOwgPB19fsDRP5G/XNSTEJND7f3YYuX0BAwdCjx4ZXQVFUbKBtLTkawAXpZSXAYQQK4E2gC7ISyn9nzn/IFA4PQv5MRrpOZKzt8+yrfM2hnccTsGCBZk6dSrTp8OmTVqW4Jo1YYvbDsIPhtNpjgv5vm8HlSrBH39kdPEVRckm0hLkCwFXn9kO5/Wt9P7A9tQOCCEGAYMAihYtmsYifnw8znkw5+gcRtYaic/fPgQHB7Njxw7Ons3FqFHQrh18+SUELgnk2Nxj1P6qJuWXfgcJCbB6NZiYZHQVFEXJJtIS5FOb4Z/qXQRCCBe0IF83teNSyvloXTlUq1Yt9TsRPnI3Ht6gv0d/KuevTCvTVrhMcmHAgAFUr94UR0coUgQWLYKbATfYMngLdp/Y0VDuggMHYNUqKF06o6ugKEo2kpabocKBIs9sFwauv3iSEMIBWAC0kVLeSZ/ifVySZTK9N/bm0eNHLGqxiEEDBlGoUCEmT57CwIFw86bWUDdMimFV+1WYWZvRoY8ZetOmaKknO3fO6CooygdLNRwaGkrFihXfb2XewNvbG0tLSxwdHXFwcKBRo0a6HD3vy8dQ72elJcgfAUoJIeyEEIZAV8Dj2ROEEEWB9UBPKeX59C/mx+GPg3+w6/IupjWdxr9//Mu5c+dYtGgRa9ZYsH49/PYbVKmczLpu64i+EU3nmc6YDRsAVavClCkZXXxFAbJ2quFn0ww84ezsTGBgIEFBQVSvXp1Zs2ZlQMkyzhu7a6SUiUKIocBOtCmUi6SUp4UQg1OOzwV+AvICs1NujU6UUmapbFuBNwP5bs93tCnThorxFXGb6sbgwYMpUqQRbdpAw4ZalmCvH724vOsyreY0p9C4z0FKrXmfku9CUXSGD4eUlAHpxtERpk9P8+kZkWp4z549jBw5ksTERKpXr86cOXM4ceIEEyZMYP369WzatImuXbsSFRVFcnIy5cuX5/Lly1y6dIkhQ4Zw69YtTE1N+euvvyhbtix9+vQhT548BAQEUKVKFaa8okElpeThw4eULFkS0PLo9OvXj8uXL2Nqasr8+fNxcHDA3d0dc3NzRo4cCUDFihV1qQ2aN29O3bp18ff3p1ChQmzatAkTExOOHTtGv379MDU1pW7dVHurM0ya5slLKbcB217YN/eZ5wOAAelbtI9HTEIMn677lLwmeZnRaAaNajeiaNGijB07kaZNwdhYmw9/3uMsfr/5UXlAZaqcWQZHjsD69VC8eEZXQVFe6UOmGo6Li6NPnz7s2bOH0qVL06tXL+bMmcPQoUMJCAgAwNfXl4oVK3LkyBESExNxctLmeQwaNIi5c+dSqlQpDh06xOeff87evXsBLY/N7t270dfXf+k9fX19cXR05M6dO5iZmfHbb78B8PPPP1O5cmU2btzI3r176dWrly5Pz6tcuHCBFStW8Ndff9G5c2fWrVtHjx496Nu3LzNnzqR+/fp88803afpdfCjqjtc0+MbzG87cPoNnD0+mjZvGhQsX2LNnD5Mm5eTYMdiwAYwf3WFprw0UrFaQFi5x0H2G1lJr1y6ji698rN6ixf2+yfeQajg1586dw87OjtIpExB69+7NrFmzGD58OCVLluTMmTMcPnyYr776in379pGUlISzszPR0dH4+/vTqVMn3bXi4+N1zzt16pRqgAetu+ZJS/z3339n1KhRzJ07Fz8/P9atWwfAJ598wp07d4iKinpt+e3s7HQJ1p6kGo6KiuL+/fvUr18fgJ49e7J9e6oTDDOEykL5BpvPbWb20dl8XetrjK4Z8ccffzBkyBCE+ISJE2HQIGjZJIHVHVejb6hP56lOGLgNBCcn+P33jC6+orzR+0o1nJpnv0xe5OzszPbt28mRIweNGjXCz88PPz8/6tWrR3JyMrly5dKlGQ4MDNRlt4S0pxpu3bo1+/bte2VZhBCvTTVs9Ey365NUw1LKjzbNMKgg/1o3Ht6gn0c/HPM7MsZpDH379sXOzo5RoybQs6c2G3LqVNg2dBuRpyJpv8gVyy/7gL6+Nl3S0DCjq6Aor/U+Uw2npmzZsoSGhnLx4kUA/vnnH10LuF69ekyfPp1atWqRL18+7ty5w9mzZ6lQoQIWFhbY2dmxZs0aQAvQJ068feYUPz8/SqSsofxsemBvb2+srKywsLDA1taW48ePA3D8+HFCQkJee81cuXJhaWmJn58fgO6aHwvVXfMKyTKZPpv68OjxI5a3X477j+5cvnwZb28fvvrKnMhI8PCA86sDCFwciPMPzpTcPhMCAmDzZm0hEEX5CH2oVMOgdc8ULvz0Bvhp06axePFiOnXqpBt4HTx4MABOTk5ERERQr149QPsSsba21n35LFu2DDc3N8aOHUtCQgJdu3bVfeG8zpM+eSkllpaWLFiwAAB3d3ddemBTU1P+/vtvQJsOunTpUhwdHalevbqua+l1Fi9erBt4bdq06RvP/5DE6/58ep+qVasmjx49miHvnRbTD05nxM4RzGk5h7LRZXFxcWHYsGE4OPzBgAEwcSL0ahrBAqcFFKldhB6DjNHr2gVGjVLdNMornTlzhnIqMZ2SRqn9exFCHHub2YuqJZ+KuUfnMtJzJK1Kt6J7me44ODhQsmRJ+vb9jTp14JNPYOjAeBbUWI1xbmPaT3RCr1FNrR/+mRtEFEVRMpoK8s9ISk7i293fMuXAFFqWasmy9sv4dsS3hIWFsWfPPgYONEuZLinZMsiDe5fv0XtXD8y/6guJibBsGeTIkdHVUBRF0VFBPsWjx4/osaEHG89u5IsaXzC16VR8vHyYM2cOX331FZ6edTl6VJv2fm3DYYLXBNNwQkOKHVwF+/ZpE+VTBnQURVE+FirIo82iabWiFQE3A/ij2R8McxrGgwcP6NevH6VLl6Zx47G0aAEDBkD1guEs7uJJadfS1GlgCHV/gi5doGfPjK6GoijKS7J9kA+KCMJ1uSt3Y++yqesmXEu7IqXk66+/Jjw8nG3b/Ojf34RSpWD8jzH867yWnAVz0nZWI0TD2lCwIMydCx/xPFlFUbKvbB3kd1zcQec1nclplBPfvr5ULlAZKSWjRo1iwYIFjBr1LQsW1CIiAg74SzzdNhJ9M5q+fn0xcf8WLl8Gb2/IlSujq6IoipKqbHsz1Jwjc3Bd7kqJPCU4NOAQlQtUJikpCTc3NyZPnsyQIUMoVeo31q7VJszE7vbjwrYLNJnahEKh/rB4MYweDc7OGV0VRXkrQgi+/vpr3fbkyZOfS072ITRo0IDUplC/av+HZGtri729PY6Ojtjb27Np06b3/p7vs97ZLsgnJSfx9c6v+Xzb5zQr2Qzfvr4UtihMQkICvXr1Yt68eYwePZphw2YyfLgeLi7QqXooXj94UbFrRaq3yq/lMnBygp9/zujqKMpbMzIyYv369dy+ffudXp9aOt/MKikpKdX9Xl5eBAYGsnbtWoYNG/aBS5W+slV3zaPHj+i+vjubzm1iWI1hTG06FX09feLi4ujSpQseHh6MHz+egQO/o3FjLSvB3MnRrG+xljyl8uA6pzmibQs1XVJJFzuG7+Bm4M10vWZ+x/w0m97stecYGBgwaNAgpk2bxrhx4547FhYWRr9+/bh16xb58uVj8eLFFC1a9KV0vnfu3MHExISzZ88SFhbG4sWL+fvvvzlw4ABOTk4sWbIEADc3N44cOUJsbCwdO3ZMNW/9m7wqJbC9vT2+vr5YWlpiZWXFtGnT6NWrFz179qR37964uLjw3Xff4e3tTXx8PEOGDOGzzz7D29ub//3vfxQoUIDAwECCg4Nf+d4PHjwgd+7cuu2pU6eyaNEiAAYMGMDw4cMJDQ3F1dWVU6dOAdpfRtHR0bi7u9OgQQOcnJzw8vLi/v37LFy4EGdnZ2JjY+nbty/BwcGUK1eO2NjYt/69pFW2CfLXH16n9YrWBNwMYEazGXzh9AUA0dHRtG3blj179jBr1iwaN/6cWrUgLAzWrUnmwNfriH8QT6/dvTCaOwN8fLSuGjVdUsnEhgwZgoODgy5H/BNDhw6lV69e9O7dm0WLFjFs2DDdilDPpvPt06cP9+7dY+/evXh4eNCqVSv279/PggULqF69OoGBgTg6OjJu3Djy5MlDUlISDRs2JCgoCAcHh7cq66tSAtepU4f9+/dTrFgxihcvjq+vL7169eLgwYPMmTOHhQsXYmlpyZEjR4iPj6dOnTo0adIEgMOHD3Pq1Cns7OxSfU8XFxeklFy+fJnVq1cDcOzYMRYvXsyhQ4eQUuLk5ET9+vWf+xJITWJiIocPH9YtzrJ7927mzJmDqakpQUFBBAUFUaVKlbf6nbyNbBHkT9w8QasVrbgbexePrh60LN0SgPv379OiRQsOHTrE33//jZ1dL2rWBD092LsX4rZ5EeodSpslbbCODYMff9SW8OvdO4NrpGQFb2pxv08WFhb06tWLGTNmYPLMwvIHDhzQLePXs2fP574EXkzn26pVK4QQ2NvbY2Njg729PQAVKlQgNDQUR0dHVq9ezfz580lMTOTGjRsEBwe/dZB/VUpgZ2dn9u3bR7FixXBzc2P+/Plcu3aNPHnyYG5ujqenJ0FBQaxduxaAqKgoLly4gKGhITVq1HhlgAetu8bKyopLly7RsGFDGjRogJ+fH+3atdNlvGzfvj2+vr60bt36teVv37498DQ1McC+fft03UAODg5v/Tt5G1m6T15KyYxDM3Ba4ESyTMavn58uwEdGRuLi4sLRo0dZs2YNenq9aNQI8uWDgwfBOuqCtgBI/8o4digJ3btDgQJquqSSZQwfPpyFCxfy6NGjV57zbGbKF9P5Pkm7q6en91wKXj09PRITEwkJCWHy5Mns2bOHoKAgWrZs+Vza3rR6VUrgevXq4evri6+vLw0aNCBfvnysXbsW55TJEFJKZs6cqUtNHBISomvJpzU1cYkSJbCxsSE4OPiVaZJfl5oYnv6enqQmfrYOH0KWDfIR0RG0XN6SL3d8SeMSjQn4LADH/I4AhIeHU69ePc6dO4eHx2aCgtrTsyfUrg0HDkDCqbOs67YOm0o2NJ/ZXFv84+JF+OcfeMOfZoqSWeTJk4fOnTuzcOFC3b7atWuzcuVKQMv6+F+Wsnvw4AFmZmZYWloSERHxzgtpvColcJEiRbh9+zYXLlygePHi1K1bl8mTJ+uCfNOmTZkzZw4JCQmA1t30ui+01ERGRhISEkKxYsWoV68eGzduJCYmhkePHrFhwwacnZ2xsbEhMjKSO3fuEB8fr1ugJK11OnXqFEFBQW9VrreRJbtrtl/YTp9NfXgQ/4A/m//J59U/131rPvnz6969e2zZ4snixXX591/o0wfmzErG130P/pP8KVitIJ3XdybHNg9YuBDGjIGUvNeKklV8/fXX/Pnnn7rtGTNm0K9fPyZNmqQbeH1XlSpVonLlylSoUIHixYuneWGRli1bkiNlUkOtWrWYN29eqimBQUtP/GSGjLOzM6NHj9Z9MQ0YMIDQ0FCqVKmClJJ8+fLpxhfexMXFBX19fRISEpgwYQI2NjbY2NjQp08fatSoobt+5cqVAfjpp59wcnLCzs6OsmXLvvH6bm5uujo5OjrqrvleSCkz5FG1alWZ3mITYuWwbcMk7kj72fbyVMSp546fOnVKFihQQObNm1fu3n1UOjtLCVKOGydl1LUHcpHzIumOu9zitkUmxCVIeeWKlLlzS1m9upSPH6d7eZXsJzg4OKOLoGQiqf17AY7Kt4i1WaYlfzryNN3WdeNk5EmG1RjG741/x9jAWHf86NGjNG3aFCMjI5Yu3cfgweW5ehVWrgQnm1DmV1nL44ePafdvOxy6O0BSEvTqBY8fw/LlarqkoiiZUqYP8lJKZh+ZzchdI7EwsmDrp1tpUarFc+fs27cPV1dX8ubNy/jxe+jZszh6erBnt0Tu38/ST/eSp1Qeeu3phXUFbZ1LJk3SUhYsWgQlS374iimKoqSDTB3kbz26RT+Pfmw5v4XmJZuzuM1ibMxtnjtn8+bNdO7cGVtbW9zcdtOrVyFKlID1y2MJ/Hkj5zefp0LnCrRa0AqjnEYgJUybpvXBd+qkddYriqJkUpk2yHte8qT3xt7cjb3LH83+4IsaX7w0Jenvv/+mf//+VKlShQYNtvHll1a4uMCcn26wo8NqHoQ/oNmMZtQYWkN7bWIifPGFNk2yQwdYskRNl1QUJVPLdEE+PjGeMXvGMPXgVMrnK8/OHjtxsHn5RoIpU6YwcuRI6tdvRK5c65k0KSd9+0g+q36cNc22Y5bPjL77+lK4Zsoiww8eaHnhd+yAb7+F337T7opSFEXJxDJdkP836F+mHpzK59U+Z3KTyZjkMHnuuJSS0aNH8/vvv9OgQScuXvyHmzeNmDgugZLntrJjyAmKNy5O+2XtMcuXckPElSvg6grBwTB/PgwcmAE1UxRFSX+Zrqnat3Jf9vXZx6yWs14K8ImJiQwcOJDff/+dKlUG4+29AmNjI3b+E4n5igUE/XOC+j/Xp/v27k8D/NGjWkbJsDCtFa8CvJLFfcyphm1tbZ/Ljunt7Y2rq+uHLBru7u4UKlQIR0dHypYti5ub23N3tL4PS5YsYejQoe/l2pkuyOsJPZyLvZzDPS4uTnf3XqFCP3L8+GwGdnzAb9U34td9Lg9vPKT79u40cG+Ann5KtTduhHr1wMgI/P2hUaMPWxlFyQAq1fDzUks3PGLECF2GypMnT+Lj45MBJUsfma67JjUPHjygbdu2eHl5YWj4B+Jhb6Y32sKDjYGcN9Cj5oia1Pm2ztPW+5MZNCNHQo0asGkT2Ni8/k0UJZ0NHz6cwMDAdL2mo6Mj06dPf+05mS3V8BPu7u5cunSJa9eucfXqVUaNGsXAgQPx9vbmp59+Im/evJw7d4569eoxe/Zs9PT08PT05OeffyY+Pp4SJUqwePFizM3NsbW1pV+/fnh6ejJ06FC6du2a6ns+fvyYuLg4XabJwMBABg8eTExMDCVKlGDRokXkzp2bBg0aMHnyZKpVq8bt27epVq0aoaGhLFmyBA8PD2JiYrh06RLt2rVj4sSJACxevJjx48dToEABSpcu/Vz+n/SU6VryL4qMjKRePRe8vX0x4i96WpViUNxMHu47QdXBVRl2eRhNJjd5GuATE+Hzz+Hrr7UZNF5eKsAr2c6QIUNYtmwZUVFRz+1/kmo4KCiI7t27P7dgxpNUw1OmTAHQpRqeNm0arVq1YsSIEZw+fZqTJ0/qvrzGjRvH0aNHCQoKwsfH5z/naAkKCmLr1q0cOHCAX375hevXrwNa6uApU6Zw8uRJLl26pPtLZezYsezevZvjx49TrVo1pk6dqruWsbExfn5+qQb4adOm4ejoqAvAjo6OAPTq1Yvff/+doKAg7O3t0/SlFRgYyKpVqzh58iSrVq3i6tWr3Lhxg59//pn9+/eza9eu1+a0/68ydUs+NDQUZ+cmXLsWjqN0p02OCETEdRz7VabeD/WwLGr5/AsePNBSBe/cCd99B+PGqRk0SoZ5U4v7ffpYUw2nlpnx2X1t2rTBxMQEExMTXFxcOHz4MLly5aJGjRoUL14cgG7duuHn54exsTHBwcG6nDmPHz+mVq1aumt16dLlleUYMWIEI0eOJCEhgY4dO7Jy5UqaN2/O/fv3qZ+Sw6p379506tTpldd4omHDhlhaarGofPnyhIWFcfv2bV3mzCdlOX/+/Buv9S7SFOSFEM2APwB9YIGUcsILx0XK8RZADNBHSnk8ncv6nMDAU9RzbkJs9EN60Qs7vSQqdXOg3k/1yFMiz8sveDKD5swZWLAA+vd/n8VTlI/e8OHDqVKlCn379n3lOemRavjIkSPkzp2bPn36vDHVcN68ebl37x5WVlaAtirUk+cvlufZ7dT2Sylp3LgxK1asSPW90pJuOEeOHDRr1ox9+/bRvHnzV573bLrhV6UahufTDX80qYaFEPrALKA5UB7oJoQo/8JpzYFSKY9BwJx0LudzVv7jjVPV2sjoaAbSk6ZtXBji1422/3MkT+x1LSH87t2wYYOWHnjmTG0GzZUr2gwaFeAV5aNMNdygQQP++ecfQBsQ/ffff3FxcdEd37RpE3Fxcdy5cwdvb2+qV68OaN01ISEhJCcns2rVKurWrUvNmjXZv38/Fy9eBCAmJuatW8tSSvz9/SlRogSWlpbkzp0bX19fAP755x9dq97W1pZjx44B6BYpeR0nJye8vb25c+cOCQkJrFmz5q3K9TbS0pKvAVyUUl4GEEKsBNoAz3YitQGWpmRIOyiEyCWEKCClvJHeBR7XaAz/2zOJnFjwo149uolN2GyaA29aUL14cdizB8q/+P2kKNnXx5Zq+Mcff8TNzY1KlSohpaRZs2b06NFDd7xGjRq0bNmSK1eu8OOPP1KwYEHOnz9PrVq1+O677zh58iT16tWjXbt26OnpsWTJErp160Z8fDwAY8eOpXTp0m8sx7Rp0/j3339JSEjAwcGBzz//HNDuon8y8Fq8eHHd72fkyJF07tyZf/75h08++eSN1y9QoADu7u7UqlWLAgUKUKVKlVcuKv5fCfmK1U50JwjREWgmpRyQst0TcJJSDn3mnC3ABCmlX8r2HuBbKeXRF641CK2lT9GiRauGhYW9dYG3/PA33/z2M4vq1qdW5VyQMyeYm7/68eR43rxgkKmHIJQs4MyZM5QrVy6ji5Epubu7Y25uzsiRI5/b7+3tzeTJk9O0WEdmk9q/FyHEMSlltbReIy1RL7WOoxe/GdJyDlLK+cB8gGrVqr3+2+UVXMf2xnWsWmNVURQlLdIS5MOBIs9sFwauv8M5iqIo7+xVd+U2aNCABg0afNCyZCZpmT94BCglhLATQhgCXQGPF87xAHoJTU0g6n30xytKVvCmLlJFgfT7d/LGlryUMlEIMRTYiTaFcpGU8rQQYnDK8bnANrTpkxfRplC+ek6WomRjxsbG3Llzh7x5836wKXRK5iOl5M6dOxgbG7/55Dd448Dr+1KtWjWZWoIiRcnKEhISCA8Pf+N8cUUxNjamcOHCukXNn3gfA6+KoqSTHDlyYGdnl9HFULIRdU+/oihKFqaCvKIoShamgryiKEoWlmEDr0KIW8Db3/KqsQLebcWDrCE71z871x2yd/1V3TXFpJT50vrCDAvy/4UQ4ujbjC5nNdm5/tm57pC966/q/m51V901iqIoWZgK8oqiKFlYZg3y8zO6ABksO9c/O9cdsnf9Vd3fQabsk1cURVHSJrO25BVFUZQ0UEFeURQlC/vog7wQQl8IEZCy+tSLx4QQYoYQ4qIQIkgIUSUjyvi+vKHuDYQQUUKIwJTHTxlRxvdFCBEqhDiZUreXMtll5c8+DXXP6p99LiHEWiHEWSHEGSFErReOZ+XP/k11f+vPPjMkKPsSOANYpHLs2QXEndAWEHf6cEV7715XdwBfKaXrByzPh+YipXzVzS9Z/bN/Xd0ha3/2fwA7pJQdU9awMH3heFb+7N9Ud3jLz/6jbskLIQoDLYEFrzhFt4C4lPIgkEsIUeCDFfA9SkPds7ss+9lnZ0IIC6AesBBASvlYSnn/hdOy5Gefxrq/tY86yAPTgVFA8iuOFwKuPrMdnrIvK5jO6+sOUEsIcUIIsV0IUeHDFOuDkYCnEOJYygLwL8rKn/2b6g5Z97MvDtwCFqd0VS4QQpi9cE5W/ezTUnd4y8/+ow3yQghXIFJKeex1p6WyL9PPCU1j3Y+j5bCoBMwENn6Isn1AdaSUVdD+NB8ihKj3wvEs+dmneFPds/JnbwBUAeZIKSsDj4DvXjgnq372aan7W3/2H22QB+oArYUQocBK4BMhxL8vnJNVFxB/Y92llA+klNEpz7cBOYQQVh+8pO+JlPJ6ys9IYANQ44VTsupn/8a6Z/HPPhwIl1IeStleixb4XjwnK372b6z7u3z2H22Ql1KOllIWllLaoi0evldK2eOF07LkAuJpqbsQIr9IWSRUCFED7bO888EL+x4IIcyEEDmfPAeaAKdeOC1LfvZpqXtW/uyllDeBq0KIMim7GgLBL5yWJT/7tNT9XT77zDC75jkiGy8g/kLdOwJuQohEIBboKrPO7cs2wIaUf8sGwHIp5Y5s8tmnpe5Z+bMH+AJYljK75DLQN5t89vDmur/1Z6/SGiiKomRhH213jaIoivLfqSCvKIqShakgryiKkoWpIK8oipKFqSCvKIqShakgryiKkoWpIK8oipKF/R8icaPvRz5g0AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#faster\n",
    "\n",
    "n=len(r)\n",
    "alpha=0.05\n",
    "eps=np.sqrt((1/(2*n))*np.log(2/alpha))\n",
    "\n",
    "F=lambda x: sum(r<x)/n\n",
    "\n",
    "xx=np.array(sorted(r))\n",
    "\n",
    "F_n=np.array(list(map(F,xx)))\n",
    "L_n=np.maximum(F_n-eps,0)\n",
    "U_n=np.minimum(F_n+eps,1)\n",
    "Sde=np.sqrt((F_n*(1-F_n))/n)\n",
    "NL=F_n+norm.ppf(alpha/2)*Sde\n",
    "NU=F_n-norm.ppf(alpha/2)*Sde\n",
    "\n",
    "plt.plot(xx,F_n, color='blue', label='Estimated CDF')\n",
    "plt.plot(xx,U_n, color='green', label='DKW Upper Bound')\n",
    "plt.plot(xx,L_n, color='red',label='DKW Lower Bound')\n",
    "plt.plot(xx,NL,color='purple', label='Normal Lower Bound')\n",
    "plt.plot(xx,NU,color='black', label='Normal Upper Bound')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% confidence interval: (0.526, 0.588)\n"
     ]
    }
   ],
   "source": [
    "z_95=norm.ppf(1-alpha/2)\n",
    "theta=F(4.9)-F(4.3)\n",
    "se=np.sqrt(theta*(1-theta)/n)\n",
    "u=theta+z_95*se\n",
    "l=theta-z_95*se\n",
    "\n",
    "print('95%% confidence interval: (%.3f, %.3f)' % (l, u))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.8**.  Get the data on eruption times and waiting times between eruptions of the old faithful geyser from the course website.\n",
    "\n",
    "- Estimate the mean waiting time and give a standard error for the estimate.\n",
    "- Also, give a 90% confidence interval for the mean waiting time.\n",
    "- Now estimate the median waiting time.\n",
    "\n",
    "In the next chapter we will see how to get the standard error for the median."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated mean: 70.897\n",
      "Estimated SE: 0.823\n"
     ]
    }
   ],
   "source": [
    "data=pd.read_csv('geysers.csv', sep=',')\n",
    "r=np.array(data['waiting'])\n",
    "n=len(r)\n",
    "mu=r.mean()\n",
    "se=r.std()/np.sqrt(n)\n",
    "\n",
    "print(\"Estimated mean: %.3f\" % mu)\n",
    "print(\"Estimated SE: %.3f\" % se)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90% confidence interval: (69.544,72.250)\n"
     ]
    }
   ],
   "source": [
    "z_90=norm.ppf(.95)\n",
    "u=mu+z_90*se\n",
    "l=mu-z_90*se\n",
    "print('90%% confidence interval: (%.3f,%.3f)' %(l,u))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimated median: 76.000\n"
     ]
    }
   ],
   "source": [
    "median=np.median(r)\n",
    "print('estimated median: %.3f' %median)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.9**.  100 people are given a standard antibiotic to treat an infection and another 100 are given a new antibiotic.  In the first group, 90 people recover; in the second group, 85 people recover.  Let $p_1$ be the probability of recovery under the standard treatment, and let $p_2$ be the probability of recovery under the new treatment.  We are interested in estimating $\\theta = p_1 - p_2$.  Provide an estimate, standard error, an 80% confidence interval and a 95% confidence interval for $\\theta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimated mean: 0.050\n",
      "standard error: 0.047\n"
     ]
    }
   ],
   "source": [
    "p1=0.9\n",
    "p2=0.85\n",
    "n=100\n",
    "\n",
    "theta=p1-p2\n",
    "se=np.sqrt(p1*(1-p1)/n+p2*(1-p2)/n)\n",
    "\n",
    "print('estimated mean: %.3f' %theta)\n",
    "print('standard error: %.3f' %se)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "80% CI: (-0.010, 0.110)\n",
      "95% CI: (-0.041, 0.141)\n"
     ]
    }
   ],
   "source": [
    "z80=norm.ppf(.9)\n",
    "z95=norm.ppf(.975)\n",
    "\n",
    "print('80%% CI: (%.3f, %.3f)' %(theta-z80*se, theta+z80*se))\n",
    "print('95%% CI: (%.3f, %.3f)' %(theta-z95*se, theta+z95*se))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 8.5.10**. In 1975, an experiment was conducted to see if cloud seeding produced rainfall.  26 clouds were seeded with silver nitrate and 26 were not.  The decision to seed or not was made at random.  Get the data from the provided link.\n",
    "\n",
    "Let $\\theta$ be the difference in the median precipitation from the two groups.\n",
    "\n",
    "- Estimate $\\theta$.\n",
    "- Estimate the standard error of the estimate and produce a 95% confidence interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}