Chapter 11.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a2d82185",
   "metadata": {},
   "source": [
    "# Chapter 11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "36fa9e14",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.pyplot import subplots\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from ISLP.models import ModelSpec as MS\n",
    "from ISLP import load_data\n",
    "from lifelines import (KaplanMeierFitter,CoxPHFitter)\n",
    "from lifelines.statistics import (logrank_test,multivariate_logrank_test)\n",
    "from ISLP.survival import sim_time\n",
    "import scipy\n",
    "import scipy.stats as sc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa907989",
   "metadata": {},
   "source": [
    "## Exercise 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a38bf558",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['sex', 'diagnosis', 'loc', 'ki', 'gtv', 'stereo', 'status', 'time'], dtype='object')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BrainCancer=load_data('BrainCancer')\n",
    "BrainCancer.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "325121ee",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sex</th>\n",
       "      <th>diagnosis</th>\n",
       "      <th>loc</th>\n",
       "      <th>ki</th>\n",
       "      <th>gtv</th>\n",
       "      <th>stereo</th>\n",
       "      <th>status</th>\n",
       "      <th>time</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Female</td>\n",
       "      <td>Meningioma</td>\n",
       "      <td>Infratentorial</td>\n",
       "      <td>90</td>\n",
       "      <td>6.11</td>\n",
       "      <td>SRS</td>\n",
       "      <td>0</td>\n",
       "      <td>57.64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Male</td>\n",
       "      <td>HG glioma</td>\n",
       "      <td>Supratentorial</td>\n",
       "      <td>90</td>\n",
       "      <td>19.35</td>\n",
       "      <td>SRT</td>\n",
       "      <td>1</td>\n",
       "      <td>8.98</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Female</td>\n",
       "      <td>Meningioma</td>\n",
       "      <td>Infratentorial</td>\n",
       "      <td>70</td>\n",
       "      <td>7.95</td>\n",
       "      <td>SRS</td>\n",
       "      <td>0</td>\n",
       "      <td>26.46</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Female</td>\n",
       "      <td>LG glioma</td>\n",
       "      <td>Supratentorial</td>\n",
       "      <td>80</td>\n",
       "      <td>7.61</td>\n",
       "      <td>SRT</td>\n",
       "      <td>1</td>\n",
       "      <td>47.80</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Male</td>\n",
       "      <td>HG glioma</td>\n",
       "      <td>Supratentorial</td>\n",
       "      <td>90</td>\n",
       "      <td>5.06</td>\n",
       "      <td>SRT</td>\n",
       "      <td>1</td>\n",
       "      <td>6.30</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      sex   diagnosis             loc  ki    gtv stereo  status   time\n",
       "0  Female  Meningioma  Infratentorial  90   6.11    SRS       0  57.64\n",
       "1    Male   HG glioma  Supratentorial  90  19.35    SRT       1   8.98\n",
       "2  Female  Meningioma  Infratentorial  70   7.95    SRS       0  26.46\n",
       "3  Female   LG glioma  Supratentorial  80   7.61    SRT       1  47.80\n",
       "4    Male   HG glioma  Supratentorial  90   5.06    SRT       1   6.30"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BrainCancer.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5903dbdf",
   "metadata": {},
   "source": [
    "### (a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "bc7eefd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Find the confidence level of z-score=1\n",
    "alpha=1-sc.norm.cdf(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "dbf3f0d5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp8AAAKnCAYAAADAwSbIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABJEklEQVR4nO3de3yT5f3/8fedpEkaKUGpFCoFqtUKMhWLc2DxuBXRH9M5J5sOhoIP8YxMNxlTFLehTvmiU3B8PaBDhW2yzQPqqhOEwSZU2BggakWLWOjKtC1tSZrk/v3Bmi89J2165/R6Ph55PMyd+3AFb+i713Vf18cwTdMUAAAAYAFbvBsAAACA9EH4BAAAgGUInwAAALAM4RMAAACWIXwCAADAMoRPAAAAWIbwCQAAAMsQPgEAAGAZR7wbEIlQKKTPP/9cWVlZMgwj3s0BAABAK6Zpqq6uTrm5ubLZOu7fTIrw+fnnnysvLy/ezQAAAEAXdu/ercGDB3f4eVKEz6ysLEmHvkzfvn3j3BoAAAC0Vltbq7y8vHBu60hShM/mofa+ffsSPgEAABJYV49IMuEIAAAAliF8AgAAwDKETwAAAFgmKZ75BACgt5imqUAgoGAwGO+mAAnNbrfL4XD0eNlLwicAIG35/X5VVlaqoaEh3k0BkoLH49GgQYPkdDq7fQ7CJwAgLYVCIe3atUt2u125ublyOp0UMgE6YJqm/H6//v3vf2vXrl06/vjjO11IvjOETwBAWvL7/QqFQsrLy5PH44l3c4CEl5mZqYyMDH366afy+/1yu93dOg8TjgAAaa27vTdAOorF3xf+xgEAAMAyhE8AANDG0qVL1a9fv3g3I2qrV6+WYRj68ssv492UiEydOlWXXHJJvJthKcInAABJpL2w8vvf/15ut1sPPPBAfBrVA1OnTpVhGJoxY0abz66//noZhqGpU6dGfL6xY8eqsrJSXq83hq3suU8++USGYWjLli0ttj/88MNaunRpr18/kUIu4RMAgCT2xBNP6Morr9Sjjz6qH/3oR/FuTrfk5eVp+fLlamxsDG87ePCgXnjhBQ0ZMiSqczmdTg0cOLBHKxf4/f5uHxstr9eblD3MPUH4BAAgST3wwAO68cYb9fzzz2v69Onh7cuWLdPo0aOVlZWlgQMH6oorrlBVVVX48+ah6VdffVWnnHKK3G63zjjjDG3durXDa5WXl+viiy9WTk6O+vTpo9NPP11vvvlmi32GDRumX/ziF7r66quVlZWlIUOGaMmSJV1+j9NOO01DhgzRypUrw9tWrlypvLw8jRo1qsW+pmnqgQce0LHHHqvMzEydcsop+v3vf9/mux0+7L5+/XqdddZZyszMVF5enm6++WbV19e3aPfPfvYzTZ06VV6vV9dcc0277ezq2l988YWuvPJKHX300crMzNTxxx+vp59+WpKUn58vSRo1apQMw9A555wjqW2P5DnnnKObbrpJM2fO1JFHHqmcnBwtWbJE9fX1uuqqq5SVlaXjjjtOr732WviYYDCoadOmKT8/X5mZmSosLNTDDz8c/vzuu+/WM888oz/96U8yDEOGYWj16tWSpD179mjSpEk68sgj1b9/f1188cX65JNPOvm/1XOETwAAdChYNPgDcXmZphl1e++44w7de++9euWVV/Ttb3+7xWd+v1/33nuv/vGPf+iPf/yjdu3a1e7Q9e23364HH3xQGzdu1IABA/TNb35TTU1N7V7vwIEDuvDCC/Xmm29q8+bNGj9+vCZOnKiKiooW+z300EMaPXq0Nm/erOuvv17XXXed3n///S6/z1VXXRUOapL01FNP6eqrr26z309/+lM9/fTTWrx4sbZt26Zbb71V3//+97VmzZp2z7t161aNHz9el156qf75z39qxYoVWrdunW688cYW+/3yl7/UyJEjVVZWpjvvvLPdc3V17TvvvFPbt2/Xa6+9ph07dmjx4sXKzs6WJL377ruSpDfffFOVlZUtgnZrzzzzjLKzs/Xuu+/qpptu0nXXXafvfOc7Gjt2rN577z2NHz9ekydPDhdHCIVCGjx4sH77299q+/btuuuuu/STn/xEv/3tbyVJt912my6//HJdcMEFqqysVGVlpcaOHauGhgade+656tOnj9555x2tW7dOffr00QUXXNCrvb+G2Z073mK1tbXyer2qqalR3759490cAEAKOHjwoHbt2qX8/Hy53W41+AMacdcbcWnL9nnj5XFGtvT21KlT9cILL8jv9+utt97Seeed1+UxGzdu1Fe/+lXV1dWpT58+Wr16tc4991wtX75ckyZNkiT95z//0eDBg7V06VJdfvnlWrp0qWbOnNnpxJ2TTjpJ1113XTjIDRs2TOPGjdNvfvMbSYcC/cCBA3XPPfe0+0xn8/f58ssv9cQTT2jw4MF6//33ZRiGTjzxRO3evVvTp09Xv379tHTpUtXX1ys7O1t/+ctfNGbMmPA5pk+froaGBj3//PPh7/bFF1+oX79+mjJlijIzM/XrX/86vP+6det09tlnq76+Xm63W8OGDdOoUaP0hz/8ocPvGsm1v/nNbyo7O1tPPfVUm+M/+eQT5efna/PmzTr11FPbfP8//vGPkg71fAaDQa1du1bSoV5Nr9erSy+9VM8++6wkae/evRo0aJA2bNigr33ta+2294YbbtC+ffvCPbOtryMdCvgPPPCAduzYEX5Mwe/3q1+/fvrjH/+okpKSNudt/ffmcJHmNRaZBwAgyZx88smqrq7WXXfdpdNPP11ZWVktPt+8ebPuvvtubdmyRf/5z38UCoUkSRUVFRoxYkR4v8ND1FFHHaXCwkLt2LGj3WvW19frnnvu0SuvvKLPP/9cgUBAjY2NbXo+Tz755PB/G4ahgQMHthjy70h2drYuuugiPfPMMzJNUxdddFG417DZ9u3bdfDgQX3jG99osd3v97cZnm9WVlamjz76SM8991x4m2ma4QpXw4cPlySNHj260/ZFcu3rrrtO3/72t/Xee++ppKREl1xyicaOHdvld2/t8D9Du92u/v376ytf+Up4W05OjiS1+HN9/PHH9cQTT+jTTz9VY2Oj/H5/i5DbnuY/m9b3z8GDB1VeXh51uyNF+AQAQFJmhl3b542P27Wjccwxx+jFF1/UueeeqwsuuECvv/56OEDU19erpKREJSUlWrZsmY4++mhVVFRo/PjxEQ2ldjRR5/bbb9cbb7yhBx98UAUFBcrMzNRll13W5pwZGRltztccfrty9dVXh3tRH3vssTafN5/n1Vdf1THHHNPiM5fL1e45Q6GQrr32Wt18881tPjt8MtMRRxzRadsiufaECRP06aef6tVXX9Wbb76p888/XzfccIMefPDBTs/dWnt/hodva/5/1Nym3/72t7r11lv10EMPacyYMcrKytIvf/lL/f3vf+/yOxUVFbUI5s2OPvroqNocDcInAAA69AM90qHvRDBkyBCtWbNG5557rkpKSvTGG2+ob9++ev/991VdXa377rtPeXl5kqRNmza1e46//e1v4QD2xRdf6IMPPtCJJ57Y7r5r167V1KlT9a1vfUvSoWdAYz0x5fBnDcePb/uLwIgRI+RyuVRRUaGzzz47onOedtpp2rZtmwoKCnrUtkivffTRR2vq1KmaOnWqxo0bF36u1ul0Sjo0jB5ra9eu1dixY3X99deHt7XuuXQ6nW2ufdppp2nFihUaMGCApY81MuEIAIAkNXjwYK1evVr79+9XSUmJampqNGTIEDmdTv3qV7/Sxx9/rJdeekn33ntvu8fPmzdPb731lv71r39p6tSpys7O7nAtyIKCAq1cuVJbtmzRP/7xD11xxRUR92hGym63a8eOHdqxY4fs9ra9wVlZWbrtttt066236plnnlF5ebk2b96sxx57TM8880y75/zxj3+sDRs26IYbbtCWLVv04Ycf6qWXXtJNN90UVdsiufZdd92lP/3pT/roo4+0bds2vfLKK+Fh/QEDBigzM1Ovv/669u3bp5qamij/dDpWUFCgTZs26Y033tAHH3ygO++8Uxs3bmyxz7Bhw/TPf/5TO3fuVHV1tZqamnTllVcqOztbF198sdauXatdu3ZpzZo1uuWWW/TZZ5/FrH2tET4BAEhixxxzjNasWaMvv/xS3/jGN5SRkaGlS5fqd7/7nUaMGKH77ruvw2Hf++67T7fccouKiopUWVmpl156KdxD19r//M//6Mgjj9TYsWM1ceJEjR8/XqeddlrMv0/fvn077YW79957ddddd2n+/PkaPny4xo8fr5dffjm8lFFrJ598stasWaMPP/xQ48aN06hRo3TnnXdq0KBBUbetq2s7nU7Nnj1bJ598ss466yzZ7XYtX75ckuRwOPTII4/o17/+tXJzc3XxxRdHff2OzJgxQ5deeqkmTZqkM844Q/v372/RCypJ11xzjQoLCzV69GgdffTR+utf/yqPx6N33nlHQ4YM0aWXXqrhw4fr6quvVmNjY6/2hDLbHQCQljqbtZvqWs8IByIVi9nuUfd8vvPOO5o4caJyc3NlGEaLKfsdWbNmjYqKiuR2u3Xsscfq8ccfj/ayAAAASAFRh8/6+nqdcsopevTRRyPaf9euXbrwwgs1btw4bd68WT/5yU90880368UXX4y6sQAAAEhuUU/rmzBhgiZMmBDx/o8//riGDBmihQsXSpKGDx+uTZs26cEHH2xTkSERmKapxqZDs8EyM+w9qg0LAEAiOuecc7pVVQmIhV5fU2LDhg1tVsgfP368nnzySTU1NbVZy0qSfD6ffD5f+H1tbW1vNzOssSkYrnBxal4//eH6sQRQAACAGOn12e579+4Nr8TfLCcnR4FAQNXV1e0eM3/+fHm93vCreZ0yq23Z/aUqvzwYl2sDAACkIkuWWmrdc9jc1d9Rj+Ls2bNVU1MTfu3evbvX29gsM8OuTT/9evh97cEmHfAFOnwdbIr9YrEAAOsw/AxELhZ/X3p92H3gwIHau3dvi21VVVVyOBzq379/u8e4XK4Oy2T1tkMVLv5vYdute77UB/vqOtzf43KouCBb7ihLowEA4qv5sa+GhgZlZmbGuTVAcmhoaJDUtgRoNHo9fI4ZM0Yvv/xyi21//vOfNXr06B413CoOu019XO230xcIqsEXUCDEb80AkGzsdrv69eunqqoqSZLH4+EZf6ADpmmqoaFBVVVV6tevX7sVqCIVdfg8cOCAPvroo/D7Xbt2acuWLTrqqKM0ZMgQzZ49W3v27NGzzz4r6dCq+48++qhmzZqla665Rhs2bNCTTz6pF154oduNtpLbblems+M/4KZgbEuLAQCsM3DgQEkKB1AAnevXr1/47013RR0+N23apHPPPTf8ftasWZKkH/zgB1q6dKkqKytVUVER/jw/P1+rVq3Srbfeqscee0y5ubl65JFHEnKZJQBAejEMQ4MGDdKAAQPU1NQU7+YACS0jI6NHPZ7Nog6fXa0NtnTp0jbbzj77bL333nvRXippNPqjm3TksBk8IwoACcRut8fkhyqArvX6M5/JzhcIdTijPRAM6YCvSRvK218yqiNMUgIAAOmK8NmFmb/d0unnxw/oozsuODHih9SZpAQAANIZ4bMdmRl2jRrST5srvuxy3w+rDsgW5TA6k5QAAEC6Iny2wzAMPXv1V/X6vyrVx5khdzuz3X2BkGYsK4tD6wAAAJIX4bMDhmHI5bDLlWHn2UwAAIAYIXzGSesZ8syABwAA6YDwaTG7zVC9P9Bmhjwz4AEAQDogfFrM6bAp1+tR8LDZ7syABwAA6YLwGQdOh63NNmbAAwCAdNA2BQEAAAC9hJ7PGPAFouu1dDlsES9KDwAAkEoInzEQ7XqfhTlZmjtxBAEUAACkHYbdu8nlsKkwJ6tbx+7cV9dub2mjP6gDvkCHteQBAACSHT2f3WQYhuZOHBHVkHtHVZFaL7/EsksAACBVET57wDBiszD84csvsewSAABIZYTPBHH48kssuwQAAFIVz3wCAADAMoRPAAAAWIbwCQAAAMvwzGeCavRHttySwxabSU8AAABWIHzGSeslmpqrHrVedqkrLMsEAACSCeEzTlqv99lc9ejwZZe6wrJMAAAg2RA+LdRcFWnnvro2nzVXPXJn2Fssu9QVlmUCAADJhPBpofaqInVU9QgAACAVET4tFquqSAAAAMmIpZYAAABgGcJnCvAHeO4TAAAkB8JnEmtelum9ii90sCmydUEBAADiifCZxJwOm446wsVySwAAIGkQPpNcho3/hQAAIHkw2z2BNC/B1FztCAAAINUQPrvgC3T+LKXdZkS1KHxnmtf7bK52RAAFAACphvDZAYfNkMflUIMv0GkVoXp/QLleT7cDaHtVjw6vdgQAAJBKCJ8dcGfYVVyQ3elEnkZ/UBvKqyOqw96Rw6seUe0IAACkOsJnJ6zqeaTqEQAASBdMlQYAAIBlCJ8AAACwDOETAAAAliF8AgAAwDKETwAAAFiG8AkAAADLsNRSgmoutdkV0+z+GqMAAABWI3wmqEgXmy8Y0Ec3nXuc/IGQ5OrlRgEAAPQQw+4JpLnUZjQ+qjqgLxub9F7FFzrY1HkdegAAgHij5zMGfIHoQp/dZrRbC/7wUptdX/P/SnEe6XGqwRfotBQoAABAIiB89oDDZsjjcqjBF1BTMLJnNCWp3h9QrtfTYQCNttRmhs2mplDk1wcAAIgXwmcPuDPsKi7IjqrHsdEf1IbyagXppQQAAGmI8NlD0fZSAgAApDMmHAEAAMAyhE8AAABYhmH3FNLoPzTr3mGLftISAACAFQifKcBmM1R/MKAN5dWSJI/LoeKCbAIoAABIOAy7p4BfvLZDg/pmqo8rQxl2G2t+AgCAhEX4TFIuh01D+3skSZ/ub5ApKdNpl8tBbycAAEhchM8kZRiG7p54UrybAQAAEBXCJwAAACxD+AQAAIBlmO0eJ75AsMV7u81ot9Y7AABAKiF8WsxhM+RxOdTgC6gpGApvr/cHlOv1EEABAEBKI3xazJ1hV3FBdoulkBr9QW0or1aQ5ZEAAECKI3zGAYu/AwCAdMUYLwAAACxD+AQAAIBlCJ8AAACwDOETAAAAliF8AgAAwDKETwAAAFiG8Jmi/IFQ1zsBAABYjPCZYuw2Q/X+gN6r+EIHm4JdHwAAAGAhwmeKcTpsOuoIlxp8gRZVlAAAABIBFY4SiC9wqKfSbjN6VOM9w2aTT/R6AgCAxEPPZwJw2Ax5XA41BUM64GvS5zUNPLMJAABSEj2fCcCdYVdxQbYCIVON/qA2lFcryJA5AABIQfR8Jgh3hl19XA5lOu3dOv7ul7fJNAmsAAAgsXUrfC5atEj5+flyu90qKirS2rVrO93/scce0/Dhw5WZmanCwkI9++yz3WosWnI5bBra3yNJ+nR/g3wM1QMAgAQXdfhcsWKFZs6cqTlz5mjz5s0aN26cJkyYoIqKinb3X7x4sWbPnq27775b27Zt0z333KMbbrhBL7/8co8bn+4Mw9DdE0/q8PNGf5DllgAAQEKJOnwuWLBA06ZN0/Tp0zV8+HAtXLhQeXl5Wrx4cbv7/+Y3v9G1116rSZMm6dhjj9V3v/tdTZs2Tffff3+PG4/2Na/1uaG8Wus+qiaAAgCAhBFV+PT7/SorK1NJSUmL7SUlJVq/fn27x/h8Prnd7hbbMjMz9e6776qpqanDY2pra1u8EDmnw6Zcr0cZdhvrfQIAgIQSVfisrq5WMBhUTk5Oi+05OTnau3dvu8eMHz9eTzzxhMrKymSapjZt2qSnnnpKTU1Nqq6ubveY+fPny+v1hl95eXnRNBM6FEBdju5NXgIAAOgt3ZpwZBhGi/emabbZ1uzOO+/UhAkT9LWvfU0ZGRm6+OKLNXXqVEmS3d5+OJo9e7ZqamrCr927d3enmQAAAEgwUYXP7Oxs2e32Nr2cVVVVbXpDm2VmZuqpp55SQ0ODPvnkE1VUVGjYsGHKyspSdnZ2u8e4XC717du3xQsAAADJL6rw6XQ6VVRUpNLS0hbbS0tLNXbs2E6PzcjI0ODBg2W327V8+XL9v//3/2SzscwoAABAOom6wtGsWbM0efJkjR49WmPGjNGSJUtUUVGhGTNmSDo0ZL5nz57wWp4ffPCB3n33XZ1xxhn64osvtGDBAv3rX//SM888E9tvgvA6ny6HrcPHIAAAAOIp6vA5adIk7d+/X/PmzVNlZaVGjhypVatWaejQoZKkysrKFmt+BoNBPfTQQ9q5c6cyMjJ07rnnav369Ro2bFjMvkQq8gUiWx4pEPy/heVnLCuTJBXmZGnuxBG90i4AAICeMMwkqMlYW1srr9ermpqalH/+82BTUOs+qlaDLxDR/gd8TXr6r5/qw6oDLbY/PfV0meahz8cclx1R2U6HzZA7gxnyAAAgepHmtah7PtG73Bl2FRdkR7Q2Z6M/qA3l1brjghNlsxnyBULh3k+p5WLzkfC4HCouyCaAAgCAXkP4TEDRhj/DaL/Hsnmx+WAEQdYXCLIgPQAA6HWEzxTndES+okDTYc+PAgAA9AbWOgIAAIBlCJ8AAACwDOETAAAAliF8AgAAwDKETwAAAFiG2e5oodHftrISi88DAIBYIXxCUucL0rP4PAAAiBXCJyR1vCA9i88DAIBYInymgKZQSJnqea9kRwvSs/g8AACIFSYcJTGHzZDH5dB/6n3yBwiIAAAg8RE+k5g7w67ThhypI5yOiOq3AwAAxBvhM8lFU7sdAAAg3kguAAAAsAzhEwAAAJZhtjsi0t7i8xIL0AMAgOgQPtGpzhafl1iAHgAARIdh9xR198vbZJo9nwHfvPh8H1dGm1eG3cYC9AAAICr0fKYQl8Omof09+nR/gz7d3yBfIBSTHsnOZtSzAD0AAIgGPZ8pxDAM3T3xpHg3AwAAoEOETwAAAFiG8AkAAADLED4BAABgGcInAAAALEP4RI/5A8x4BwAAkSF8otuaF6B/r+ILHWxqvwISAADA4QifKcwXCMVkofmOOB02HXWEi4XmAQBAxAifKWzGsjLd8/L2Xg2gGTZuIQAAEDmSQ4pxOWwqzMkKv9+5r04+nskEAAAJgvCZYgzD0NyJI/T494vi3RQAAIA2CJ8pyDAMuTqpxw4AABAvJBQAAABYhvAJAAAAyzji3QCkhkb//63z6bAZcmfY49gaAACQqAifKaIpFFKmrA98zQvNbyivDm/zuBwqLsgmgAIAgDYIn0nOYTPkcTn077qDcjvsclo80cjpsCnX61Hwv4vM+wJBFp0HAAAdInwmOXeGXacNOVLrPvx3OAC2Fuk6ny6HTYZhRN2G1oG3Kci6ogAAoH2EzxTQVW/njGVlEZ2nMCdLcyeO6FYABQAAiASz3VNU60pHkaAaEgAA6G30fKao5kpHkYRJXyAUce8oAABATxA+U5hhsOQRAABILAy7AwAAwDKETwAAAFiGYXf0isMrHnUHVZIAAEhNhE/EVHsVj7qDKkkAAKQmwidiqnXFo+6gShIAAKmL8JlCfIFDQ912m2F5mc3DxeLaVEkCACA1ET5TQHN99wZfQE3BkOr9AeV6Pd0Kge2tC9rdspsAAACtET5TgDvDruKCbAVCphr9QW0or+72sHd7i81TdhMAAMQKSy2lCHeGXX1cDmU6o5+g01UpTspuAgCAWKHnEx2W4qTsJgAAiDXCJyRRihMAAFiDYXcAAABYhp5PJKzDqyRR8QgAgNRA+ETCaa9KEhWPAABIDYRPJJzWVZKoeAQAQOogfCIhtV4gn4pHAACkBsInItLROp9UPwIAANEgfCIiHa33SfUjAAAQDZZaQoe6qnwkUf0IAABEh55PdKijykcS1Y8AAED3ED7RKSofAQCAWGLYHQAAAJah5xNJ4/CKR12hIhIAAImJ8ImE117Fo65QEQkAgMRE+ETCa13xqCtURAIAIHERPpEUWlc86goVkQAASExMOAIAAIBlCJ8AAACwDOETPXb3y9tkmjxfCQAAukb4TFFNod595tHlsGlof48k6dP9DZTYBAAAESF8phiHzZDH5dB/6n3y92IgNAxDd088qdfODwAAUhPhM8W4M+w6bciROsLpiHhpIgAAAKt0K3wuWrRI+fn5crvdKioq0tq1azvd/7nnntMpp5wij8ejQYMG6aqrrtL+/fu71WB0LdpliVJVoz+oA76ADvgCOtgUeXUkAADQe6JOKStWrNDMmTM1Z84cbd68WePGjdOECRNUUVHR7v7r1q3TlClTNG3aNG3btk2/+93vtHHjRk2fPr3HjQfac3hFpL/s2Ke/7NindR9VE0ABAEgAUYfPBQsWaNq0aZo+fbqGDx+uhQsXKi8vT4sXL253/7/97W8aNmyYbr75ZuXn56u4uFjXXnutNm3a1OPGA+1projUx5WhPq4MZdhtVDwCACBBRBU+/X6/ysrKVFJS0mJ7SUmJ1q9f3+4xY8eO1WeffaZVq1bJNE3t27dPv//973XRRRd1eB2fz6fa2toWLyAaTodNmU67Mp12uRzUdwcAIFFEFT6rq6sVDAaVk5PTYntOTo727t3b7jFjx47Vc889p0mTJsnpdGrgwIHq16+ffvWrX3V4nfnz58vr9YZfeXl50TQTAAAACapbM1MMw2jx3jTNNtuabd++XTfffLPuuusulZWV6fXXX9euXbs0Y8aMDs8/e/Zs1dTUhF+7d+/uTjMBAACQYBzR7JydnS273d6ml7OqqqpNb2iz+fPn68wzz9Ttt98uSTr55JN1xBFHaNy4cfrZz36mQYMGtTnG5XLJ5XJF0zTE2eGLzLsctg5/GQEAAOktqvDpdDpVVFSk0tJSfetb3wpvLy0t1cUXX9zuMQ0NDXI4Wl7Gbj/0DB4lGVPHjGVl4f8uzMnS3IkjCKAAAKCNqIfdZ82apSeeeEJPPfWUduzYoVtvvVUVFRXhYfTZs2drypQp4f0nTpyolStXavHixfr444/117/+VTfffLO++tWvKjc3N3bfBJZzOWwqzMlqs33nvjrKbQIAgHZF1fMpSZMmTdL+/fs1b948VVZWauTIkVq1apWGDh0qSaqsrGyx5ufUqVNVV1enRx99VD/84Q/Vr18/nXfeebr//vtj9y3QLl8gKLvN6LVF5w3D0NyJI8JB0xcItegBBQAAaM0wk2Dsu7a2Vl6vVzU1Nerbt2+8m5PwDjYFte6jajX4Aqr3B5Tr9VhS9ehgU1BXLd0oSXp66ulyZyTGEkeHKh016bzhOerjivr3LQAAEIFI8xo/iVOQO8Ou4oJs1R08VOWHGu+HNPpjV+HIYTMSJlwDAJBMCJ8pyp1hp6LPfx1ebjNWPC6HiguyCaAAAESJ8ImU11xuM1Y9wL5AkHKdAAB0E+ETaSHWz7w2BZnNDwBAd/T+LBQAAADgvwifAAAAsAzD7ugVPV1knhKdAACkJsInekVPF5unRCcAAKmJYXfETEflNruDEp0AAKQmej4RM63LbXYHJToBAEhthM804AtEVtknFnXgDYPKPwAAoGOEzxTmsBnyuBxq8AUiWpfSyjrwqaB1uU5KbgIA0DXCZwprrvEeSSWeRn+QOvAR6qhcJyU3AQDoGuEzxRGEYq+9cp2U3AQAIDKET6Ab2ns0gZKbAAB0jYf7AAAAYBnCJwAAACxD+AQAAIBleOYTCevwxeqp9Q4AQGogfCJhHV7piFrvAACkBobdkVA6qg9PrXcAAFIDPZ9oob1SnLEouxmp1vXhqfUOAEBqIXxCUuelOK0uu5nM9eFbl9xsjRKcAIB0R/iEpI5LcVJ2MzIdldxsjRKcAIB0R/hEGIGo+9orudkaJTgBACB8AjETyWMJlOAEAKQ7ZrsDAADAMoRPAAAAWIbwCQAAAMvwzCeSRjSLzFOOEwCAxET4RNKIZrF5ynECAJCYGHZHQuuo3GZXKMcJAEBioucTCa11uc2uUI4TAIDERvhERNqr+S5ZU/c9mcttAgCAlgif6FRnNd8l6+u+p4Ku6r9Hi3rxAIBkQvhEpzqq+S5R9z1akdZ/jxb14gEAyYTwiS4RamIjkvrv0aJePAAg2RA+AQv1xuMJ1IsHACQTHtQDAACAZQifAAAAsAzhEwAAAJbhmU+krMMXpqfWOwAAiYHwiZR1eKUjar0DAJAYGHZHSumoFjy13gEASAz0fCKltK4FT613AAASC+ETPdYUCilTibMQPbXgAQBIXIRPdFtz3fd/1x2U22GnvnscxbpefDSoLQ8AiAbhE93mzrDrtCFHat2H/6a+e5z0Vr34aFBbHgAQDcIneoTezvjqjXrx0aC2PAAgWoRPIMnF+xcAassDAKJBtxUAAAAsQ/gEAACAZQifAAAAsAzPfCJtxLLCEbXiAQDoHsIn0kYsKx1RKx4AgO5h2B0praNa7z1FrXgAALqHnk+ktNa13nuKWvEAAPQM4RMx4Qv8X3lHu82I+9qTh6PWOwAAiYPwiR5pru/e4AuEFxuv9weU6/UkVAAFAACJgfCJHnFn2FVckB0ur9joD2pDeTW13tNMoz/Y5T4OGz3QAADCJ2KAQJG+7DZD9f6ANpRXd7mvx+VQcUE29wsApDnCJ4BuczpsyvV6uuzp9gWCavAFwj3kAID0RfgE0CORPtvb/EwwACC9MSMEAAAAliF8AgAAwDIMuwPd1Hrheuq9AwDQNcIn0E2tKx1R7x0AgK4x7A5EobNa8dR7BwCga/R8AlFor1Y89d4BAIgc4RO94vBa792RaPXhD0eteAAAuo/wiZhqr9Z7d1AfHgCA1ET4REy1rvXeHdSHBwAgdRE+EXMMSQMAgI4wpgkAAADLdCt8Llq0SPn5+XK73SoqKtLatWs73Hfq1KkyDKPN66STTup2owEAAJCcog6fK1as0MyZMzVnzhxt3rxZ48aN04QJE1RRUdHu/g8//LAqKyvDr927d+uoo47Sd77znR43HgAAAMkl6vC5YMECTZs2TdOnT9fw4cO1cOFC5eXlafHixe3u7/V6NXDgwPBr06ZN+uKLL3TVVVf1uPEAAABILlGFT7/fr7KyMpWUlLTYXlJSovXr10d0jieffFJf//rXNXTo0GguDQAAgBQQ1Wz36upqBYNB5eTktNiek5OjvXv3dnl8ZWWlXnvtNT3//POd7ufz+eTz+cLva2tro2kmEDddldd0OWzUfgcApLVuLbXU+oenaZoR/UBdunSp+vXrp0suuaTT/ebPn6977rmnO00D4qqrMpuFOVmaO3EEARQAkLaiGnbPzs6W3W5v08tZVVXVpje0NdM09dRTT2ny5MlyOp2d7jt79mzV1NSEX7t3746mmYClXA6bCnOyItp35766LntHAQBIZVH1fDqdThUVFam0tFTf+ta3wttLS0t18cUXd3rsmjVr9NFHH2natGldXsflcsnlckXTNCBuDMPQ3IkjOg2VvkCoy15RAADSQdTD7rNmzdLkyZM1evRojRkzRkuWLFFFRYVmzJgh6VCv5Z49e/Tss8+2OO7JJ5/UGWecoZEjR8am5Uh5vkAw/N92m5HQdd4Nw6CyEwAAEYg6fE6aNEn79+/XvHnzVFlZqZEjR2rVqlXh2euVlZVt1vysqanRiy++qIcffjg2rUZKc9gMeVwONfgCagoe6k2s9weU6/UkdAAFAABdM0zTNOPdiK7U1tbK6/WqpqZGffv2jXdzYIGDTUEFQoduzUZ/UBvKq9XHlaFMZ3L2Lh5sCuqqpRslSU9PPT3tekkb/UEd8DXpvOE56uPq1jxHAECCizSv8VMACSndwlm6aPQHu96pHQ4bjzUAQKogfALodXaboXp/QBvKq7t1vMflUHFBNgEUAFIA4RNAr3M6bMr1ehQMRf+Ujy8QVIMvEH4MAwCQ3AifACzRk8lizRPPAADJj6nDAAAAsAzhEwAAAJZh2B2wWKzLa7ocNmrFAwCSBuETsFisy2wW5mRp7sQRBFAAQFJg2B2wgMthU2FOVq+ce+e+upj3pgIA0Fvo+UTSOLzWe1cSrRa8YRiaO3FETEOiLxCKeS8qAAC9jfCJhNderfeuJGIteMOgSg8AAIRPJDx3hl3FBdkRLzLeXAu+OwuaAwCA3kX4RFKgxxAAgNSQOGOSAAAASHmETwAAAFiGYXcASaHRH/lqB+1x2JjwBQCJgPAJIKHZbYbq/QFtKK/u0Xk8LoeKC7IJoAAQZ4RPAAnN6bAp1+vp0eoFvkBQDb5AxCsmAAB6D+ETQMKLxXqtka4RCwDoXUw4AgAAgGXo+QRSQDxru7scNhmGEbfrAwCSC+ETSAHxrPFemJOluRNHEEABABFh2B0pyxcIqtF/6OWPY89gb3E5bCrMyYp3M7RzX11ce14BAMmFnk+kHIfNkMflUIMvEJ5kUu8PKNfricnElURhGIbmThwRt+DnC4Ti2uMKAEhOhE+kHHeGXcUF2eFldRr9QW0or+7RUj2JyjBYOB0AkFwIn0hJBDIAABJT6oxBAgAAIOERPgEAAGAZwicAAAAsQ/gEAACAZQifAAAAsAyz3QGkjUZ/sFfP77Cx9BUAdIXwCSDl2W2G6v0BbSiv7tXreFwOFRdkE0ABoBOETwApz+mwKdfr6dVCA75AUA2+QLi4AQCgfYRPAD0WSYlPl8MmwzAsaE37rCit2lzOFQDQMcIn0oYvELvn/ew2I6XqxPdUJDXeC3OyNHfiiLgGUABA/BE+kfIcNkMel0MNvkDMeqbq/QHlej1pHUBdDpsKc7K0c19dRPvv3FcnXyDE85AAkOYIn0h57gy7iguyY/YsXqM/qA3l1b36/GAyMAxDcyeO6HLI3RcIRdQzCgBID4RPpAV623qHYbC0EAAgOuk7ZggAAADLET4BAABgGcInAAAALEP4BAAAgGUInwAAALAM4RMAAACWIXwCAADAMoRPAAAAWIbwCQAAAMtQ4QgAYqjRH4x3E6LisFGlCoC1CJ9AN/kCLUOG3WbI6WAwIV3ZbYbq/QFtKK+Od1Oi4nE5VFyQTQAFYBnCJxAlh82Qx+VQgy+gpmAovL3eH1Cu10MATVNOh025Xo+CITPeTYmYLxBUgy+gQBK1GUDyI3wCUXJn2FVckN3iB3ajP6gN5dVJFTwQe8n4i8fhv0ABgBUIn0A3MEQJAED3JN+v6QAAAEhahE8AAABYhvAJAAAAyxA+AQAAYBnCJwAAACxD+AQAAIBlCJ8AAACwDOETAAAAliF8AgAAwDKETwAAAFiG8ppADPkCwU4/t9uMpKz/HSu+QPfqiLscNhmGEePWAADigfAJxIDDZsjjcqjBF1BTsOOAVe8PKNfrSdsAOmNZWbeOK8zJ0tyJIwigAJACCJ9ADLgz7CouyFYgZHa4T6M/qA3l1Qp2sk8qcjlsKszJ0s59dd0+x859dfIFQnJn2GPYMgBAPBA+gRghGLXPMAzNnTiiW0PuvkCo272lAIDERPgE0OsMwyCcAwAkMdsdAAAAFiJ8AgAAwDIMuwNAmmv0d75EWDOHjccnAPQc4RMA0pTdZqjeH9CG8uqI9ve4HCouyCaAAugRwicApCmnw6Zcryei5b98gaAafIFOlxMDgEgQPgEgjUVT8KCzAgoAECkmHAEAAMAy3QqfixYtUn5+vtxut4qKirR27dpO9/f5fJozZ46GDh0ql8ul4447Tk899VS3GgwAAIDkFfWw+4oVKzRz5kwtWrRIZ555pn79619rwoQJ2r59u4YMGdLuMZdffrn27dunJ598UgUFBaqqqlIgEOhx4wEAAJBcog6fCxYs0LRp0zR9+nRJ0sKFC/XGG29o8eLFmj9/fpv9X3/9da1Zs0Yff/yxjjrqKEnSsGHDetZqIIn5ApEtaxMpu82I6rk9AADiKaqfWH6/X2VlZSopKWmxvaSkROvXr2/3mJdeekmjR4/WAw88oGOOOUYnnHCCbrvtNjU2NnZ4HZ/Pp9ra2hYvINk5bIY8LoeagiEd8DXF7PV5TYP83aibDgBAPETV81ldXa1gMKicnJwW23NycrR37952j/n444+1bt06ud1u/eEPf1B1dbWuv/56/ec//+nwuc/58+frnnvuiaZpQMJzZ9hVXJAd06VqGv1BbSivjmipHAAAEkG3lloyDKPFe9M022xrFgqFZBiGnnvuOXm9XkmHhu4vu+wyPfbYY8rMzGxzzOzZszVr1qzw+9raWuXl5XWnqUBCYXHu7vP1sHfX5bB1+O8UAMA6UYXP7Oxs2e32Nr2cVVVVbXpDmw0aNEjHHHNMOHhK0vDhw2Wapj777DMdf/zxbY5xuVxyuVzRNA1AipuxrKxHxxfmZGnuxBEEUACIs6ie+XQ6nSoqKlJpaWmL7aWlpRo7dmy7x5x55pn6/PPPdeDAgfC2Dz74QDabTYMHD+5GkwGkC5fDpsKcrJica+e+uh73ngIAei7qYfdZs2Zp8uTJGj16tMaMGaMlS5aooqJCM2bMkHRoyHzPnj169tlnJUlXXHGF7r33Xl111VW65557VF1drdtvv11XX311u0PuANDMMAzNnTiiR6HRFwj1uNcUABA7UYfPSZMmaf/+/Zo3b54qKys1cuRIrVq1SkOHDpUkVVZWqqKiIrx/nz59VFpaqptuukmjR49W//79dfnll+tnP/tZ7L4FgJRlGAbPygJACunWhKPrr79e119/fbufLV26tM22E088sc1QPQAAANIPK1MDAADAMoRPAAAAWIbwCQAAAMsQPgEAAGAZwicAAAAsQ/gEAACAZbq11BKAxOILBON2bbvNkNPB77EAgMgQPoEk5rAZ8rgcavAF1BSMT+nIen9AuV4PATRNNPpj/4uOw0YhASCdED6BJObOsKu4IFuBkBmX6zf6g9pQXq1gnK4P69hthur9AW0or475uT0uh4oLsgmgQJogfAJJjh/YsILTYVOu1xPzXzR8gaAafIG4/QIFwHqETwBARHrr0Yp4PTICID54SAsAAACWoecTQNrwBXq3h83lsMkwjF69BgAkO8IngLQxY1lZr56/MCdLcyeOIIACQCcYdgeQ0lwOmwpzsiy51s59db3euwoAyY6eTwApzTAMzZ04oldDoS8Q6vVeVQBIFYRPACnPMFjEHAASBcPuAAAAsAzhEwAAAJYhfAIAAMAyhE8AAABYhglHAHrMFwh2uY/dZvRaeUYAQPIgfALoNofNkMflUIMv0GV97np/QLleDwEUANIc4RNAt7kz7CouyFYgZHa6X6M/qA3l1Qp2sR8AIPURPgH0COtnAgCiwfgXAAAALEP4BAAAgGUInwAAALAM4RMAAACWIXwCAADAMoRPAAAAWIbwCQAAAMuwzicAxJAv0Hmlp0TjcthkGEa8mwEgjRA+ASCGZiwri3cTolKYk6W5E0cQQAFYhvAJAD3kcthUmJOlnfvq4t2UqO3cVydfIBT3SlWN/mBcrw+kKofNiPvf79YInwDQQ4ZhaO7EEUk15O4LhBKil9ZuM1TvD2hDeXW8mwKkJI/LoeKC7IQKoIRPAIgBw0i83oVk4HTYlOv1KBgy490UIOX4AkE1+AIKJNjfL8InACCunA4WXgF6S1Mw8UZk+BsPAAAAyxA+AQAAYBnCJwAAACxD+AQAAIBlCJ8AAACwDOETAAAAliF8AgAAwDKETwAAAFiG8AkAAADLED4BAABgGcInAAAALEP4BAAAgGUInwAAALAM4RMAAACWccS7AQDShy8Q7NZxdpshp4PflQEgFRA+AfQ6h82Qx+VQgy+gpmAo6uPr/QHlej0EUABIAYRPAL3OnWFXcUG2AiEz6mMb/UFtKK9WsBvHAgASD+ETgCXcGfZ4NwEAkAAYwwIAAIBlCJ8AAACwDMPuAJDmfIHIJoG5HDYZhtHLrQGQ6gifAJDmZiwri2i/wpwszZ04ggAKoEcYdgeANORy2FSYkxXVMTv31UXcSwoAHaHnEwDSkGEYmjtxRERh0hcIRdw7CgBdIXwCQJoyDIMlsABYjmF3AAAAWIbwCQAAAMsQPgEAAGAZwicAAAAsQ/gEAACAZQifAAAAsAzhEwAAAJZhnU8AScEXCPboeLvNkNPB79sAEG+ETwAJzWEz5HE51OALqCnY/dKO9f6Acr0eAigAxBnhE0BCc2fYVVyQrUDI7PY5Gv1BbSivVrAH5wAAxAbhE0DCowQkAKQOxp8AAABgmW6Fz0WLFik/P19ut1tFRUVau3Zth/uuXr1ahmG0eb3//vvdbjQAAACSU9Thc8WKFZo5c6bmzJmjzZs3a9y4cZowYYIqKio6PW7nzp2qrKwMv44//vhuNxoAAADJKerwuWDBAk2bNk3Tp0/X8OHDtXDhQuXl5Wnx4sWdHjdgwAANHDgw/LLbeYYLAAAg3UQVPv1+v8rKylRSUtJie0lJidavX9/psaNGjdKgQYN0/vnn6+233+50X5/Pp9ra2hYvAAAAJL+owmd1dbWCwaBycnJabM/JydHevXvbPWbQoEFasmSJXnzxRa1cuVKFhYU6//zz9c4773R4nfnz58vr9YZfeXl50TQTAAAACapbSy0ZhtHivWmabbY1KywsVGFhYfj9mDFjtHv3bj344IM666yz2j1m9uzZmjVrVvh9bW0tARQAACAFRNXzmZ2dLbvd3qaXs6qqqk1vaGe+9rWv6cMPP+zwc5fLpb59+7Z4AQAAIPlFFT6dTqeKiopUWlraYntpaanGjh0b8Xk2b96sQYMGRXNpAAAApICoh91nzZqlyZMna/To0RozZoyWLFmiiooKzZgxQ9KhIfM9e/bo2WeflSQtXLhQw4YN00knnSS/369ly5bpxRdf1IsvvhjbbwIAAICEF3X4nDRpkvbv36958+apsrJSI0eO1KpVqzR06FBJUmVlZYs1P/1+v2677Tbt2bNHmZmZOumkk/Tqq6/qwgsvjN23AAAAQFIwTNM0492IrtTW1srr9aqmpobnPwFE7YAvoL/s2Kc+rgxlOlljOFoHm4K6aulGSdLTU0+XO4M/QyAZNPqDOuBr0nnDc9TH1a055lGJNK/1fksAACnDFwjF/Jwuh63DFVMApB7CJwAgYjOWlcX8nIU5WZo7cQQBFEgTUZfXBACkF5fDpsKcrF47/859db3SowogMdHzCQDolGEYmjtxRMwDoi8Q6pWeVACJjfAJIG34AsFePb/dZsjpSM0BJcMwmGgEICYInwBSnsNmyONyqMEXUFOw94Z36/0B5Xo9KRtAASAWCJ8AUp47w67igmwFQr23slyjP6gN5dUK9uI1ACAVED4BpAWGjAEgMTA2BAAAAMsQPgEAAGAZwicAAAAsQ/gEAACAZQifAAAAsAzhEwAAAJYhfAIAAMAyhE8AAABYhvAJAAAAyxA+AQAAYBnCJwAAACxD+AQAAIBlCJ8AAACwDOETAAAAliF8AgAAwDKETwAAAFiG8AkAAADLED4BAABgGUe8GwAAqcQXCMa7CVGx2ww5HfRDALAO4RMAYsBhM+RxOdTgC6gpGIp3cyJW7w8o1+shgAKwDOETAGLAnWFXcUG2AiEz3k2JWKM/qA3l1QomUZsBJD/CJwDEiDvDHu8mAEDCY5wFAAAAliF8AgAAwDIMuwMA4s4XSJ5JWkCy8DUF5QsEZZqJ9Vw34RMAEHczlpXFuwlAyjq7cICy3BnxbkYYw+4AgLhwOWwqzMmKdzMAWIyeTwBAXBiGobkTRzDkDvSSg/6gDviblJlgK3EQPgEAcWMYBktUAb3ENKWmUEiGYcS7KS0w7A4AAADLED4BAABgGcInAAAALEP4BAAAgGUInwAAALAM4RMAAACWIXwCAADAMoRPAAAAWIbwCQAAAMsQPgEAAGAZwicAAAAsQ/gEAACAZQifAAAAsAzhEwAAAJYhfAIAAMAyhE8AAABYhvAJAAAAyxA+AQAAYBnCJwAAACzjiHcDAADx5QsEI9rPbjPkdNBnAaBnCJ8AkKYcNkMel0MNvoCagqEu96/3B5Tr9RBAAfQI4RMA0pQ7w67igmwFQmaX+zb6g9pQXq1gBPsCQGcInwCQxtwZ9ng3AUCaYewEAAAAliF8AgAAwDKETwAAAFiG8AkAAADLED4BAABgGcInAAAALEP4BAAAgGUInwAAALAM4RMAAACWIXwCAADAMoRPAAAAWIbwCQAAAMsQPgEAAGAZwicAAAAs063wuWjRIuXn58vtdquoqEhr166N6Li//vWvcjgcOvXUU7tzWQAAACS5qMPnihUrNHPmTM2ZM0ebN2/WuHHjNGHCBFVUVHR6XE1NjaZMmaLzzz+/240FAABAcos6fC5YsEDTpk3T9OnTNXz4cC1cuFB5eXlavHhxp8dde+21uuKKKzRmzJhuNxYAAADJLarw6ff7VVZWppKSkhbbS0pKtH79+g6Pe/rpp1VeXq65c+dGdB2fz6fa2toWLwAAACS/qMJndXW1gsGgcnJyWmzPycnR3r172z3mww8/1B133KHnnntODocjouvMnz9fXq83/MrLy4ummQAAAEhQ3ZpwZBhGi/emabbZJknBYFBXXHGF7rnnHp1wwgkRn3/27NmqqakJv3bv3t2dZgIAACDBRNYV+V/Z2dmy2+1tejmrqqra9IZKUl1dnTZt2qTNmzfrxhtvlCSFQiGZpimHw6E///nPOu+889oc53K55HK5omkaAAAAkkBUPZ9Op1NFRUUqLS1tsb20tFRjx45ts3/fvn21detWbdmyJfyaMWOGCgsLtWXLFp1xxhk9az0AAACSSlQ9n5I0a9YsTZ48WaNHj9aYMWO0ZMkSVVRUaMaMGZIODZnv2bNHzz77rGw2m0aOHNni+AEDBsjtdrfZDgAAgNQXdficNGmS9u/fr3nz5qmyslIjR47UqlWrNHToUElSZWVll2t+AgAAID0Zpmma8W5EV2pra+X1elVTU6O+ffvGuzkAkHYO+AL6y4596uPKUKbTHu/mAIhAoz+oA74mnTc8R31cUfc3Ri3SvEZtdwAAAFiG8AkAAADL9H4fLAAgZfgCwZif024z5HTQFwKkC8InAKBLDpshj8uhBl9ATcFQTM9d7w8o1+shgAJpgvAJAOiSO8Ou4oJsBUKxnaPa6A9qQ3m1gjE+L4DERfgEAETEncEsdwA9xxgHAAAALEP4BAAAgGUInwAAALAM4RMAAACWIXwCAADAMoRPAAAAWIbwCQAAAMsQPgEAAGAZwicAAAAsQ/gEAACAZQifAAAAsAzhEwAAAJYhfAIAAMAyhE8AAABYhvAJAAAAyxA+AQAAYBnCJwAAACxD+AQAAIBlCJ8AAACwDOETAAAAliF8AgAAwDKETwAAAFiG8AkAAADLED4BAABgGcInAAAALEP4BAAAgGUInwAAALAM4RMAAACWIXwCAADAMoRPAAAAWIbwCQAAAMsQPgEAAGAZwicAAAAs44h3AwAA8AWC8W4CkHIS9e8V4RMAEDcOmyGPy6EGX0BNwVC8mwOkHI/LIYfNiHczWiB8AgDixp1hV3FBtgIhM95NAVKSw2bInWGPdzNaIHwCAOIq0X4wAuhdTDgCAACAZQifAAAAsAzhEwAAAJYhfAIAAMAyhE8AAABYhvAJAAAAyxA+AQAAYBnCJwAAACxD+AQAAIBlCJ8AAACwDOETAAAAliF8AgAAwDKETwAAAFiG8AkAAADLED4BAABgGcInAAAALEP4BAAAgGUInwAAALAM4RMAAACWIXwCAADAMoRPAAAAWIbwCQAAAMsQPgEAAGAZwicAAAAs44h3AyJhmqYkqba2Ns4tAQAAQHuac1pzbutIUoTPuro6SVJeXl6cWwIAAIDO1NXVyev1dvi5YXYVTxNAKBTS559/rqysLBmG0evXq62tVV5ennbv3q2+ffv2+vWQfLhH0BnuD3SFewSdSdb7wzRN1dXVKTc3VzZbx092JkXPp81m0+DBgy2/bt++fZPqfzqsxz2CznB/oCvcI+hMMt4fnfV4NmPCEQAAACxD+AQAAIBlCJ/tcLlcmjt3rlwuV7ybggTFPYLOcH+gK9wj6Eyq3x9JMeEIAAAAqYGeTwAAAFiG8AkAAADLED4BAABgGcInAAAALEP4bGXRokXKz8+X2+1WUVGR1q5dG+8mIU7mz5+v008/XVlZWRowYIAuueQS7dy5s8U+pmnq7rvvVm5urjIzM3XOOedo27ZtcWox4mn+/PkyDEMzZ84Mb+P+wJ49e/T9739f/fv3l8fj0amnnqqysrLw59wj6SsQCOinP/2p8vPzlZmZqWOPPVbz5s1TKBQK75Oq9wfh8zArVqzQzJkzNWfOHG3evFnjxo3ThAkTVFFREe+mIQ7WrFmjG264QX/7299UWlqqQCCgkpIS1dfXh/d54IEHtGDBAj366KPauHGjBg4cqG984xuqq6uLY8thtY0bN2rJkiU6+eSTW2zn/khvX3zxhc4880xlZGTotdde0/bt2/XQQw+pX79+4X24R9LX/fffr8cff1yPPvqoduzYoQceeEC//OUv9atf/Sq8T8reHybCvvrVr5ozZsxose3EE08077jjjji1CImkqqrKlGSuWbPGNE3TDIVC5sCBA8377rsvvM/BgwdNr9drPv744/FqJixWV1dnHn/88WZpaal59tlnm7fccotpmtwfMM0f//jHZnFxcYefc4+kt4suusi8+uqrW2y79NJLze9///umaab2/UHP53/5/X6VlZWppKSkxfaSkhKtX78+Tq1CIqmpqZEkHXXUUZKkXbt2ae/evS3uGZfLpbPPPpt7Jo3ccMMNuuiii/T1r3+9xXbuD7z00ksaPXq0vvOd72jAgAEaNWqU/vd//zf8OfdIeisuLtZbb72lDz74QJL0j3/8Q+vWrdOFF14oKbXvD0e8G5AoqqurFQwGlZOT02J7Tk6O9u7dG6dWIVGYpqlZs2apuLhYI0eOlKTwfdHePfPpp59a3kZYb/ny5Xrvvfe0cePGNp9xf+Djjz/W4sWLNWvWLP3kJz/Ru+++q5tvvlkul0tTpkzhHklzP/7xj1VTU6MTTzxRdrtdwWBQP//5z/W9731PUmr/G0L4bMUwjBbvTdNssw3p58Ybb9Q///lPrVu3rs1n3DPpaffu3brlllv05z//WW63u8P9uD/SVygU0ujRo/WLX/xCkjRq1Cht27ZNixcv1pQpU8L7cY+kpxUrVmjZsmV6/vnnddJJJ2nLli2aOXOmcnNz9YMf/CC8XyreHwy7/1d2drbsdnubXs6qqqo2v3Ugvdx000166aWX9Pbbb2vw4MHh7QMHDpQk7pk0VVZWpqqqKhUVFcnhcMjhcGjNmjV65JFH5HA4wvcA90f6GjRokEaMGNFi2/Dhw8OTWPk3JL3dfvvtuuOOO/Td735XX/nKVzR58mTdeuutmj9/vqTUvj8In//ldDpVVFSk0tLSFttLS0s1duzYOLUK8WSapm688UatXLlSf/nLX5Sfn9/i8/z8fA0cOLDFPeP3+7VmzRrumTRw/vnna+vWrdqyZUv4NXr0aF155ZXasmWLjj32WO6PNHfmmWe2WZ7tgw8+0NChQyXxb0i6a2hokM3WMobZ7fbwUkspfX/EcbJTwlm+fLmZkZFhPvnkk+b27dvNmTNnmkcccYT5ySefxLtpiIPrrrvO9Hq95urVq83Kysrwq6GhIbzPfffdZ3q9XnPlypXm1q1bze9973vmoEGDzNra2ji2HPFy+Gx30+T+SHfvvvuu6XA4zJ///Ofmhx9+aD733HOmx+Mxly1bFt6HeyR9/eAHPzCPOeYY85VXXjF37dplrly50szOzjZ/9KMfhfdJ1fuD8NnKY489Zg4dOtR0Op3maaedFl5WB+lHUruvp59+OrxPKBQy586daw4cONB0uVzmWWedZW7dujV+jUZctQ6f3B94+eWXzZEjR5oul8s88cQTzSVLlrT4nHskfdXW1pq33HKLOWTIENPtdpvHHnusOWfOHNPn84X3SdX7wzBN04xnzysAAADSB898AgAAwDKETwAAAFiG8AkAAADLED4BAABgGcInAAAALEP4BAAAgGUInwAAALAM4RMAOrF69WoZhqEvv/yyV6+zdOlS9evXL/z+7rvv1qmnntqr1wSAeCB8AsBhzjnnHM2cOTP8fuzYsaqsrJTX67W0HbfddpveeustS68JAFZwxLsBAJDInE6nBg4caPl1+/Tpoz59+lh+XQDobfR8AsB/TZ06VWvWrNHDDz8swzBkGIaWLl3aYti9eXj8lVdeUWFhoTwejy677DLV19frmWee0bBhw3TkkUfqpptuUjAYDJ/b7/frRz/6kY455hgdccQROuOMM7R69eoO29J62H3q1Km65JJL9OCDD2rQoEHq37+/brjhBjU1NXX7GgAQD/R8AsB/Pfzww/rggw80cuRIzZs3T5K0bdu2Nvs1NDTokUce0fLly1VXV6dLL71Ul156qfr166dVq1bp448/1re//W0VFxdr0qRJkqSrrrpKn3zyiZYvX67c3Fz94Q9/0AUXXKCtW7fq+OOPj6h9b7/9tgYNGqS3335bH330kSZNmqRTTz1V11xzTcyuAQC9jfAJAP/l9XrldDrl8XjCQ+3vv/9+m/2ampq0ePFiHXfccZKkyy67TL/5zW+0b98+9enTRyNGjNC5556rt99+W5MmTVJ5ebleeOEFffbZZ8rNzZV06JnO119/XU8//bR+8YtfRNS+I488Uo8++qjsdrtOPPFEXXTRRXrrrbd0zTXXxOwaANDbCJ8AECWPxxMOnpKUk5OjYcOGtXhGMycnR1VVVZKk9957T6Zp6oQTTmhxHp/Pp/79+0d83ZNOOkl2uz38ftCgQdq6dWtMrwEAvY3wCQBRysjIaPHeMIx2t4VCIUlSKBSS3W5XWVlZi/AoKapJRVZcAwB6G+ETAA7jdDpbTBSKhVGjRikYDKqqqkrjxo2L6bmtvAYAxAKz3QHgMMOGDdPf//53ffLJJ6qurg73LPbECSecoCuvvFJTpkzRypUrtWvXLm3cuFH333+/Vq1aFYNWW3MNAIgFwicAHOa2226T3W7XiBEjdPTRR6uioiIm53366ac1ZcoU/fCHP1RhYaG++c1v6u9//7vy8vJicn6rrgEAPWWYpmnGuxEAAABID/R8AgAAwDKETwAAAFiG8AkAAADLED4BAABgGcInAAAALEP4BAAAgGUInwAAALAM4RMAAACWIXwCAADAMoRPAAAAWIbwCQAAAMsQPgEAAGCZ/w8yvEvsOsT0PAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots(figsize=(8,8))\n",
    "km = KaplanMeierFitter()\n",
    "km_brain = km.fit(BrainCancer['time'], BrainCancer['status'],alpha=alpha)\n",
    "km_brain.plot(label='Kaplan Meier estimate', ax=ax);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "662d029b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>KM_estimate</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>timeline</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.00</th>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.07</th>\n",
       "      <td>0.988636</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.18</th>\n",
       "      <td>0.988636</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.41</th>\n",
       "      <td>0.977141</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.54</th>\n",
       "      <td>0.977141</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>65.02</th>\n",
       "      <td>0.451543</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67.38</th>\n",
       "      <td>0.451543</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>73.74</th>\n",
       "      <td>0.451543</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>78.75</th>\n",
       "      <td>0.451543</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>82.56</th>\n",
       "      <td>0.451543</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>87 rows × 1 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "          KM_estimate\n",
       "timeline             \n",
       "0.00         1.000000\n",
       "0.07         0.988636\n",
       "1.18         0.988636\n",
       "1.41         0.977141\n",
       "1.54         0.977141\n",
       "...               ...\n",
       "65.02        0.451543\n",
       "67.38        0.451543\n",
       "73.74        0.451543\n",
       "78.75        0.451543\n",
       "82.56        0.451543\n",
       "\n",
       "[87 rows x 1 columns]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "km_brain.survival_function_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe4a56bf",
   "metadata": {},
   "source": [
    "### (b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "104e1fe6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RangeIndex(start=0, stop=88, step=1)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BrainCancer.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "eb9a9677",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define indexes of the bootstrap sample\n",
    "rng=np.random.default_rng(0)\n",
    "boot_index=rng.choice(BrainCancer.index,88,replace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8853c54d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define the sample\n",
    "sample=pd.DataFrame({'obs':BrainCancer['time'].iloc[boot_index],\n",
    "                     'status':BrainCancer['status'].iloc[boot_index]})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c4076f47",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>KM_estimate</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>timeline</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.00</th>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.07</th>\n",
       "      <td>0.977273</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.18</th>\n",
       "      <td>0.977273</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.41</th>\n",
       "      <td>0.965639</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2.03</th>\n",
       "      <td>0.954004</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5.15</th>\n",
       "      <td>0.954004</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5.51</th>\n",
       "      <td>0.941928</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6.23</th>\n",
       "      <td>0.941928</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7.05</th>\n",
       "      <td>0.941928</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7.18</th>\n",
       "      <td>0.941928</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7.25</th>\n",
       "      <td>0.915764</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8.98</th>\n",
       "      <td>0.915764</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9.77</th>\n",
       "      <td>0.915764</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10.49</th>\n",
       "      <td>0.902096</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11.02</th>\n",
       "      <td>0.888217</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11.51</th>\n",
       "      <td>0.874339</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11.57</th>\n",
       "      <td>0.860460</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13.90</th>\n",
       "      <td>0.860460</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14.56</th>\n",
       "      <td>0.832249</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14.62</th>\n",
       "      <td>0.818143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16.92</th>\n",
       "      <td>0.818143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17.57</th>\n",
       "      <td>0.789436</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18.95</th>\n",
       "      <td>0.789436</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19.41</th>\n",
       "      <td>0.773957</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19.74</th>\n",
       "      <td>0.758478</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19.90</th>\n",
       "      <td>0.742340</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20.69</th>\n",
       "      <td>0.725843</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23.67</th>\n",
       "      <td>0.725843</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26.46</th>\n",
       "      <td>0.691279</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31.15</th>\n",
       "      <td>0.673997</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31.67</th>\n",
       "      <td>0.656261</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32.82</th>\n",
       "      <td>0.656261</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33.41</th>\n",
       "      <td>0.656261</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33.67</th>\n",
       "      <td>0.617657</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34.26</th>\n",
       "      <td>0.617657</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35.93</th>\n",
       "      <td>0.597733</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>36.10</th>\n",
       "      <td>0.577808</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39.25</th>\n",
       "      <td>0.577808</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39.54</th>\n",
       "      <td>0.577808</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42.07</th>\n",
       "      <td>0.577808</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42.10</th>\n",
       "      <td>0.553733</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44.39</th>\n",
       "      <td>0.553733</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>46.16</th>\n",
       "      <td>0.553733</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>47.80</th>\n",
       "      <td>0.553733</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>49.05</th>\n",
       "      <td>0.553733</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>51.02</th>\n",
       "      <td>0.522970</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54.43</th>\n",
       "      <td>0.492207</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55.80</th>\n",
       "      <td>0.492207</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>57.11</th>\n",
       "      <td>0.492207</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>57.25</th>\n",
       "      <td>0.447461</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>57.64</th>\n",
       "      <td>0.397743</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>65.02</th>\n",
       "      <td>0.397743</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67.38</th>\n",
       "      <td>0.318194</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>78.75</th>\n",
       "      <td>0.159097</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          KM_estimate\n",
       "timeline             \n",
       "0.00         1.000000\n",
       "0.07         0.977273\n",
       "1.18         0.977273\n",
       "1.41         0.965639\n",
       "2.03         0.954004\n",
       "5.15         0.954004\n",
       "5.51         0.941928\n",
       "6.23         0.941928\n",
       "7.05         0.941928\n",
       "7.18         0.941928\n",
       "7.25         0.915764\n",
       "8.98         0.915764\n",
       "9.77         0.915764\n",
       "10.49        0.902096\n",
       "11.02        0.888217\n",
       "11.51        0.874339\n",
       "11.57        0.860460\n",
       "13.90        0.860460\n",
       "14.56        0.832249\n",
       "14.62        0.818143\n",
       "16.92        0.818143\n",
       "17.57        0.789436\n",
       "18.95        0.789436\n",
       "19.41        0.773957\n",
       "19.74        0.758478\n",
       "19.90        0.742340\n",
       "20.69        0.725843\n",
       "23.67        0.725843\n",
       "26.46        0.691279\n",
       "31.15        0.673997\n",
       "31.67        0.656261\n",
       "32.82        0.656261\n",
       "33.41        0.656261\n",
       "33.67        0.617657\n",
       "34.26        0.617657\n",
       "35.93        0.597733\n",
       "36.10        0.577808\n",
       "39.25        0.577808\n",
       "39.54        0.577808\n",
       "42.07        0.577808\n",
       "42.10        0.553733\n",
       "44.39        0.553733\n",
       "46.16        0.553733\n",
       "47.80        0.553733\n",
       "49.05        0.553733\n",
       "51.02        0.522970\n",
       "54.43        0.492207\n",
       "55.80        0.492207\n",
       "57.11        0.492207\n",
       "57.25        0.447461\n",
       "57.64        0.397743\n",
       "65.02        0.397743\n",
       "67.38        0.318194\n",
       "78.75        0.159097"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "km = KaplanMeierFitter()\n",
    "km_brain = km.fit(sample['obs'], BrainCancer['status'],alpha=1)\n",
    "km_brain.survival_function_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "680c85ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now let's repeat the process 200 times\n",
    "B=200\n",
    "km_brain_list=[]\n",
    "\n",
    "for _ in range(B):\n",
    "    boot_index=rng.choice(BrainCancer.index,88,replace=True)\n",
    "    sample=pd.DataFrame({'obs':BrainCancer['time'].iloc[boot_index],\n",
    "                     'status':BrainCancer['status'].iloc[boot_index]})\n",
    "    km_brain_list.append(km.fit(sample['obs'], BrainCancer['status'],alpha=1))\n",
    "survival=km_brain_list[0].survival_function_\n",
    "for i in range(1,B):\n",
    "    survival=survival+km_brain_list[i].survival_function_\n",
    "survival=survival/B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "1ce872c7",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>KM_estimate</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>timeline</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.00</th>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.07</th>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.18</th>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.41</th>\n",
       "      <td>0.976471</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2.03</th>\n",
       "      <td>0.976471</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4.16</th>\n",
       "      <td>0.964415</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4.56</th>\n",
       "      <td>0.952208</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6.10</th>\n",
       "      <td>0.952208</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6.23</th>\n",
       "      <td>0.952208</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6.30</th>\n",
       "      <td>0.952208</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7.05</th>\n",
       "      <td>0.952208</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7.18</th>\n",
       "      <td>0.938983</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7.25</th>\n",
       "      <td>0.938983</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8.98</th>\n",
       "      <td>0.925568</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9.77</th>\n",
       "      <td>0.898740</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10.10</th>\n",
       "      <td>0.885123</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10.82</th>\n",
       "      <td>0.885123</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11.02</th>\n",
       "      <td>0.885123</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11.51</th>\n",
       "      <td>0.885123</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11.57</th>\n",
       "      <td>0.870613</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13.90</th>\n",
       "      <td>0.856103</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14.56</th>\n",
       "      <td>0.856103</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14.75</th>\n",
       "      <td>0.841342</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16.92</th>\n",
       "      <td>0.826045</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17.57</th>\n",
       "      <td>0.810160</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19.41</th>\n",
       "      <td>0.794274</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19.90</th>\n",
       "      <td>0.794274</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20.69</th>\n",
       "      <td>0.778065</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22.03</th>\n",
       "      <td>0.778065</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26.46</th>\n",
       "      <td>0.760774</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26.85</th>\n",
       "      <td>0.760774</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29.70</th>\n",
       "      <td>0.742219</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31.15</th>\n",
       "      <td>0.742219</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31.25</th>\n",
       "      <td>0.742219</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31.67</th>\n",
       "      <td>0.742219</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32.82</th>\n",
       "      <td>0.721602</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33.67</th>\n",
       "      <td>0.721602</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35.93</th>\n",
       "      <td>0.699051</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39.25</th>\n",
       "      <td>0.676501</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41.44</th>\n",
       "      <td>0.608851</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42.10</th>\n",
       "      <td>0.608851</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44.39</th>\n",
       "      <td>0.581176</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>46.16</th>\n",
       "      <td>0.581176</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>47.11</th>\n",
       "      <td>0.581176</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>47.80</th>\n",
       "      <td>0.512803</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>49.05</th>\n",
       "      <td>0.473356</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>51.02</th>\n",
       "      <td>0.473356</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>52.23</th>\n",
       "      <td>0.473356</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>52.75</th>\n",
       "      <td>0.473356</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54.43</th>\n",
       "      <td>0.414187</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55.80</th>\n",
       "      <td>0.355017</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>65.02</th>\n",
       "      <td>0.236678</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67.38</th>\n",
       "      <td>0.177509</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>73.74</th>\n",
       "      <td>0.118339</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>78.75</th>\n",
       "      <td>0.059170</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>82.56</th>\n",
       "      <td>0.059170</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          KM_estimate\n",
       "timeline             \n",
       "0.00         1.000000\n",
       "0.07         1.000000\n",
       "1.18         1.000000\n",
       "1.41         0.976471\n",
       "2.03         0.976471\n",
       "4.16         0.964415\n",
       "4.56         0.952208\n",
       "6.10         0.952208\n",
       "6.23         0.952208\n",
       "6.30         0.952208\n",
       "7.05         0.952208\n",
       "7.18         0.938983\n",
       "7.25         0.938983\n",
       "8.98         0.925568\n",
       "9.77         0.898740\n",
       "10.10        0.885123\n",
       "10.82        0.885123\n",
       "11.02        0.885123\n",
       "11.51        0.885123\n",
       "11.57        0.870613\n",
       "13.90        0.856103\n",
       "14.56        0.856103\n",
       "14.75        0.841342\n",
       "16.92        0.826045\n",
       "17.57        0.810160\n",
       "19.41        0.794274\n",
       "19.90        0.794274\n",
       "20.69        0.778065\n",
       "22.03        0.778065\n",
       "26.46        0.760774\n",
       "26.85        0.760774\n",
       "29.70        0.742219\n",
       "31.15        0.742219\n",
       "31.25        0.742219\n",
       "31.67        0.742219\n",
       "32.82        0.721602\n",
       "33.67        0.721602\n",
       "35.93        0.699051\n",
       "39.25        0.676501\n",
       "41.44        0.608851\n",
       "42.10        0.608851\n",
       "44.39        0.581176\n",
       "46.16        0.581176\n",
       "47.11        0.581176\n",
       "47.80        0.512803\n",
       "49.05        0.473356\n",
       "51.02        0.473356\n",
       "52.23        0.473356\n",
       "52.75        0.473356\n",
       "54.43        0.414187\n",
       "55.80        0.355017\n",
       "65.02        0.236678\n",
       "67.38        0.177509\n",
       "73.74        0.118339\n",
       "78.75        0.059170\n",
       "82.56        0.059170"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "survival"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "4db6f04a",
   "metadata": {},
   "outputs": [],
   "source": [
    "timeline=survival.index.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "1acfa7e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Compute the number of events at time t of the bootstrap survival function\n",
    "q=[]\n",
    "for i in range(len(timeline)):\n",
    "    q.append(np.sum(BrainCancer['status'][BrainCancer['time']<timeline[i]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "4d4c3c3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Compute the difference between the number of subjects at risk and died subjects at time t of\n",
    "#the bootstrap survival function\n",
    "r=[]\n",
    "for i in range(len(timeline)):\n",
    "    r.append(np.count_nonzero(BrainCancer['status'][BrainCancer['time']>=timeline[i]]==0)\n",
    "              +np.sum(BrainCancer['status']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "44a54ca3",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "r_q=[]\n",
    "for i in range(len(r)):\n",
    "    r_q.append(r[i]-q[i])\n",
    "rr_q=list(map(lambda x, y: x * y, r, r_q))\n",
    "inv_rr_q=list(map(lambda x, y: x / y, [1]*len(rr_q), rr_q))\n",
    "qinv_rr_q=list(map(lambda x, y: x * y, r, inv_rr_q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "45e102fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "56"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(qinv_rr_q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "id": "6c49d5ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now let's compute standard error of each each bootstraped survival observation using Greenwood formula\n",
    "ste=[]\n",
    "for i in range(len(timeline)):\n",
    "    S2=survival.iloc[i]['KM_estimate']**2\n",
    "    factor=0\n",
    "    for j in range(i+1):\n",
    "       factor+=qinv_rr_q[j] \n",
    "    ste.append(S2*factor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "id": "ddb14327",
   "metadata": {},
   "outputs": [],
   "source": [
    "ste_df=pd.DataFrame({'timeline':timeline,'Standard error':ste})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "da52d025",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>timeline</th>\n",
       "      <th>Standard error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.00</td>\n",
       "      <td>0.011364</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.07</td>\n",
       "      <td>0.022727</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.18</td>\n",
       "      <td>0.034222</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.41</td>\n",
       "      <td>0.043717</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2.03</td>\n",
       "      <td>0.055068</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>4.16</td>\n",
       "      <td>0.065060</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>4.56</td>\n",
       "      <td>0.074617</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>6.10</td>\n",
       "      <td>0.086241</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>6.23</td>\n",
       "      <td>0.098016</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>6.30</td>\n",
       "      <td>0.109947</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>7.05</td>\n",
       "      <td>0.122540</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>7.18</td>\n",
       "      <td>0.131578</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>7.25</td>\n",
       "      <td>0.144173</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>8.98</td>\n",
       "      <td>0.152499</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>9.77</td>\n",
       "      <td>0.155665</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>10.10</td>\n",
       "      <td>0.162677</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>10.82</td>\n",
       "      <td>0.174730</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>11.02</td>\n",
       "      <td>0.186971</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>11.51</td>\n",
       "      <td>0.199607</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>11.57</td>\n",
       "      <td>0.205542</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>13.90</td>\n",
       "      <td>0.210963</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>14.56</td>\n",
       "      <td>0.223599</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>14.75</td>\n",
       "      <td>0.228596</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>16.92</td>\n",
       "      <td>0.232995</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>17.57</td>\n",
       "      <td>0.236504</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>19.41</td>\n",
       "      <td>0.239690</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>19.90</td>\n",
       "      <td>0.252565</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>20.69</td>\n",
       "      <td>0.255242</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>22.03</td>\n",
       "      <td>0.268403</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>26.46</td>\n",
       "      <td>0.270387</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>26.85</td>\n",
       "      <td>0.284503</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31</th>\n",
       "      <td>29.70</td>\n",
       "      <td>0.284566</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>31.15</td>\n",
       "      <td>0.299063</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33</th>\n",
       "      <td>31.25</td>\n",
       "      <td>0.313952</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>31.67</td>\n",
       "      <td>0.329255</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35</th>\n",
       "      <td>32.82</td>\n",
       "      <td>0.326094</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>36</th>\n",
       "      <td>33.67</td>\n",
       "      <td>0.342891</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>37</th>\n",
       "      <td>35.93</td>\n",
       "      <td>0.339248</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>38</th>\n",
       "      <td>39.25</td>\n",
       "      <td>0.335316</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39</th>\n",
       "      <td>41.44</td>\n",
       "      <td>0.287724</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40</th>\n",
       "      <td>42.10</td>\n",
       "      <td>0.305376</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41</th>\n",
       "      <td>44.39</td>\n",
       "      <td>0.295134</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>46.16</td>\n",
       "      <td>0.313898</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>47.11</td>\n",
       "      <td>0.333767</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>47.80</td>\n",
       "      <td>0.276289</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>45</th>\n",
       "      <td>49.05</td>\n",
       "      <td>0.250355</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>46</th>\n",
       "      <td>51.02</td>\n",
       "      <td>0.267591</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>47</th>\n",
       "      <td>52.23</td>\n",
       "      <td>0.286263</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>48</th>\n",
       "      <td>52.75</td>\n",
       "      <td>0.306633</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>49</th>\n",
       "      <td>54.43</td>\n",
       "      <td>0.251921</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50</th>\n",
       "      <td>55.80</td>\n",
       "      <td>0.199089</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>51</th>\n",
       "      <td>65.02</td>\n",
       "      <td>0.099687</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>52</th>\n",
       "      <td>67.38</td>\n",
       "      <td>0.063951</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>53</th>\n",
       "      <td>73.74</td>\n",
       "      <td>0.033091</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54</th>\n",
       "      <td>78.75</td>\n",
       "      <td>0.010023</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55</th>\n",
       "      <td>82.56</td>\n",
       "      <td>0.013524</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    timeline  Standard error\n",
       "0       0.00        0.011364\n",
       "1       0.07        0.022727\n",
       "2       1.18        0.034222\n",
       "3       1.41        0.043717\n",
       "4       2.03        0.055068\n",
       "5       4.16        0.065060\n",
       "6       4.56        0.074617\n",
       "7       6.10        0.086241\n",
       "8       6.23        0.098016\n",
       "9       6.30        0.109947\n",
       "10      7.05        0.122540\n",
       "11      7.18        0.131578\n",
       "12      7.25        0.144173\n",
       "13      8.98        0.152499\n",
       "14      9.77        0.155665\n",
       "15     10.10        0.162677\n",
       "16     10.82        0.174730\n",
       "17     11.02        0.186971\n",
       "18     11.51        0.199607\n",
       "19     11.57        0.205542\n",
       "20     13.90        0.210963\n",
       "21     14.56        0.223599\n",
       "22     14.75        0.228596\n",
       "23     16.92        0.232995\n",
       "24     17.57        0.236504\n",
       "25     19.41        0.239690\n",
       "26     19.90        0.252565\n",
       "27     20.69        0.255242\n",
       "28     22.03        0.268403\n",
       "29     26.46        0.270387\n",
       "30     26.85        0.284503\n",
       "31     29.70        0.284566\n",
       "32     31.15        0.299063\n",
       "33     31.25        0.313952\n",
       "34     31.67        0.329255\n",
       "35     32.82        0.326094\n",
       "36     33.67        0.342891\n",
       "37     35.93        0.339248\n",
       "38     39.25        0.335316\n",
       "39     41.44        0.287724\n",
       "40     42.10        0.305376\n",
       "41     44.39        0.295134\n",
       "42     46.16        0.313898\n",
       "43     47.11        0.333767\n",
       "44     47.80        0.276289\n",
       "45     49.05        0.250355\n",
       "46     51.02        0.267591\n",
       "47     52.23        0.286263\n",
       "48     52.75        0.306633\n",
       "49     54.43        0.251921\n",
       "50     55.80        0.199089\n",
       "51     65.02        0.099687\n",
       "52     67.38        0.063951\n",
       "53     73.74        0.033091\n",
       "54     78.75        0.010023\n",
       "55     82.56        0.013524"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ste_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cae4318d",
   "metadata": {},
   "source": [
    "### (c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "id": "bf65d8da",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>se(coef)</th>\n",
       "      <th>p</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>covariate</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>sex[Male]</th>\n",
       "      <td>0.183748</td>\n",
       "      <td>0.360358</td>\n",
       "      <td>0.610119</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diagnosis[LG glioma]</th>\n",
       "      <td>-1.239530</td>\n",
       "      <td>0.579555</td>\n",
       "      <td>0.032455</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diagnosis[Meningioma]</th>\n",
       "      <td>-2.154566</td>\n",
       "      <td>0.450524</td>\n",
       "      <td>0.000002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diagnosis[Other]</th>\n",
       "      <td>-1.268870</td>\n",
       "      <td>0.617672</td>\n",
       "      <td>0.039949</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>loc[Supratentorial]</th>\n",
       "      <td>0.441195</td>\n",
       "      <td>0.703669</td>\n",
       "      <td>0.530665</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ki</th>\n",
       "      <td>-0.054955</td>\n",
       "      <td>0.018314</td>\n",
       "      <td>0.002693</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gtv</th>\n",
       "      <td>0.034293</td>\n",
       "      <td>0.022333</td>\n",
       "      <td>0.124661</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>stereo[SRT]</th>\n",
       "      <td>0.177778</td>\n",
       "      <td>0.601578</td>\n",
       "      <td>0.767597</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           coef  se(coef)         p\n",
       "covariate                                          \n",
       "sex[Male]              0.183748  0.360358  0.610119\n",
       "diagnosis[LG glioma]  -1.239530  0.579555  0.032455\n",
       "diagnosis[Meningioma] -2.154566  0.450524  0.000002\n",
       "diagnosis[Other]      -1.268870  0.617672  0.039949\n",
       "loc[Supratentorial]    0.441195  0.703669  0.530665\n",
       "ki                    -0.054955  0.018314  0.002693\n",
       "gtv                    0.034293  0.022333  0.124661\n",
       "stereo[SRT]            0.177778  0.601578  0.767597"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cleaned = BrainCancer.dropna()#delete the missing data\n",
    "coxph = CoxPHFitter # shorthand\n",
    "model_df=MS(cleaned,intercept=False).fit_transform(cleaned)\n",
    "cox_fit = coxph().fit(model_df,'time','status')\n",
    "cox_fit.summary[['coef', 'se(coef)', 'p']]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08fe14e4",
   "metadata": {},
   "source": [
    "### (d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "id": "5e790242",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 90,  70,  80, 100,  60,  40], dtype=int64)"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BrainCancer['ki'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "id": "fc830138",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "62    40\n",
       "Name: ki, dtype: int64"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cleaned['ki'][cleaned['ki']==40]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "id": "24071042",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\kabousselham\\AppData\\Local\\Temp\\ipykernel_8224\\804784288.py:1: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  cleaned['ki'] = cleaned['ki'].replace([40], [60])\n"
     ]
    }
   ],
   "source": [
    "cleaned['ki'] = cleaned['ki'].replace([40], [60])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "id": "38fad6e7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "60"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cleaned['ki'].iloc[62]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "id": "e14d3707",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp8AAAKnCAYAAADAwSbIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABjmUlEQVR4nO3deZgcVaH//8+p6q7umUkyyWTfCEEwJLJJAiSBXAFhCKDIVS9RMIACGtkMuajh4heRC4kC8gPZV1EW5QqIyA2B3CABTNhi4uWyKltYMoRANjIz3bWc3x8Dkwyzdk93dc/M+/U8/dBTfarqNKmZ/vQ5dc4x1lorAAAAIAZOqSsAAACAvoPwCQAAgNgQPgEAABAbwicAAABiQ/gEAABAbAifAAAAiA3hEwAAALEhfAIAACA2iVJXoCuiKNK7776r/v37yxhT6uoAAADgU6y12rJli0aNGiXHab99s0eEz3fffVdjx44tdTUAAADQibfeektjxoxp9/UeET779+8vqenNDBgwoMS1AQAAwKdt3rxZY8eObc5t7ekR4fOTrvYBAwYQPgEAAMpYZ7dIMuAIAAAAsSF8AgAAIDaETwAAAMSmR9zzCQAA+q4wDOX7fqmr0ee5rqtEItHtaS8JnwAAoGx99NFHevvtt2WtLXVVIKmyslIjR46U53l5H4PwCQAAylIYhnr77bdVWVmpoUOHstBMCVlrlc1m9f777+v111/XLrvs0uFE8h0hfAIAgLLk+76stRo6dKgqKipKXZ0+r6KiQslkUm+++aay2azS6XRex2HAEQAAKGu0eJaPfFs7WxyjAPUAAAAAuoTwCQAAgNgQPgEAAGJ04okn6uijjy51NUqG8AkAAFBAbYXLu+++W+l0WhdffLGuuOIK3Xrrrd06x4oVK3TwwQerqqpKAwcO1IEHHqiGhobm1zds2KDZs2erurpa1dXVmj17tjZu3NitcxYK4RMAAKCIbrrpJh133HG66qqr9KMf/UjV1dUaOHBg3sdbsWKFZs6cqdraWj399NN65plndPrpp7cYDHTsscdq9erVWrx4sRYvXqzVq1dr9uzZBXg33cdUSwAAoEew1qrBD0ty7oqkm9eo+4svvljnnXee7rzzTn3ta1+T1NQyunHjRt1333151eWss87SmWeeqfnz5zdv22WXXZqfv/jii1q8eLGefPJJ7bfffpKkG2+8UdOmTdPLL7+sCRMm5HXeQiF8AgCAHqHBDzXpvIdKcu4XLjhMlV5usWn+/Pm6+uqr9cADD+iQQw5pt9yCBQu0YMGCDo/14IMPasaMGVq3bp2eeuopHXfccZo+fbpeffVV7brrrrrooot0wAEHSGpqGa2urm4OnpI0depUVVdXa/ny5YRPAACA3ubBBx/Un/70Jy1dulQHH3xwh2XnzJmjY445psMyo0ePliS99tprkqTzzz9fl156qfbaay/99re/1Re/+EX93//9n3bZZRfV1dVp2LBhrY4xbNgw1dXV5fmOCofwCQAAeoSKpKsXLjisZOfOxR577KH169frvPPO0z777KP+/fu3W7ampkY1NTVdOm4URZKk733ve/r2t78tSfr85z+vpUuX6pZbbtHChQsltT0xv7W2LCbsZ8ARAADoEYwxqvQSJXnkGtpGjx6tZcuWae3atZo5c6a2bNnSbtkFCxaoX79+HT4ef/xxSdLIkSMlSZMmTWpxjIkTJ2rNmjWSpBEjRui9995rdZ73339fw4cPz+l9FAMtnwAAAEWwww47aNmyZTrooINUW1urhx56SAMGDGhVLpdu9x133FGjRo3Syy+/3OL1V155RYcffrgkadq0adq0aZOefvpp7bvvvpKkp556Sps2bdL06dML8da6hfAJAABQJGPGjNGjjz7aIoB+Wi7d7sYY/fCHP9RPf/pT7bnnntprr730m9/8Ri+99JLuvvtuSU2toDNnztQpp5yi66+/XpL03e9+V1/60pdKPthIyqPb/bHHHtOXv/xljRo1SsaYLk0TsGzZMk2ePFnpdFo77bSTrrvuunzqCgAA0ON80gW/ceNGHXrood2e7H3u3Lk655xzdNZZZ2nPPffU0qVLtWTJEn3mM59pLnPHHXdo9913V21trWpra7XHHnvotttu6+Y7KYycWz63bt2qPffcU9/+9reb56vqyOuvv64jjjhCp5xyim6//Xb99a9/1amnnqqhQ4d2aX8AAICepK3Vi0aOHKmXXnqpYOeYP39+i3k+P62mpka33357wc5XSDmHz8MPP7z5noKuuO6667TDDjvo8ssvl9TUFPzss8/q0ksvJXwCAAD0MUW/53PFihWqra1tse2www7TzTffLN/3lUwmi12FnAS+r/XvNo0WU7JSKsKUBJ7rtDpshZvOe/oDU1FRFlMnAEA+Ep7D3zCgDyl6+Kyrq2s1rH/48OEKgkDr169vnjJge5lMRplMpvnnzZs3F7uazda/u0b3LHwztvMBQF83fMcB+tqPJxNAgT4ilnk+P/0HxVrb5vZPLFy4UNXV1c2PsWPHFr2OAIDSeO+Nzfrow0znBQH0CkVv+RwxYkSrpZzWrVunRCKhwYMHt7nPOeeco3nz5jX/vHnz5tgC6JBRO+hrZzcodWvTGqz/O/kKRVVDC3LsTCD9+yPvS5Kuqh2qoVVJZSJfpz1/tiTpmkm/1NDKoXKMked2fjuCbcxo3Xe+JUkafucflBw9piD1BIA4BNlQd/3nM5Kk8ONVWwD0fkUPn9OmTdOf//znFtsefvhhTZkypd37PVOplFKpVLGr1qZEMqkRO4yXEg2SpC98YR+pZqeCHLs+G2jrY03ze+20934aN6RK9X69tr7mS5LWjNqkOrdBlYlKTR05VelEusPjRfX1+iDKSpIGDEopNbSyIPUEgDj4mbDUVQBQAjl3u3/00UdavXq1Vq9eLalpKqXVq1c3L+l0zjnn6Pjjj28uP2fOHL355puaN2+eXnzxRd1yyy26+eabdfbZZxfmHfQiCSchz/FUH9QrtPxRBgAAvU/OLZ/PPvusDjrooOafP+keP+GEE3Trrbdq7dq1zUFUksaPH69FixbprLPO0tVXX61Ro0bpV7/6FdMstSHlpuS5nrIft2YCAAD0NjmHzwMPPLB5wFBb2ppY9Qtf+IL+9re/5Xqq8uA3SJkt7b/uJKRkRXz1AQAAPdqJJ56ojRs3dmmVyN4oltHuPdq7q6VXHmr/8dqjTQEVAABATeHy6KOPbrHt7rvvVjqd1sUXX6wrrriizca6rnjjjTdkjGnz8Yc//KG53IYNGzR79uzmmYNmz57d7WU9C6XoA456PCcppfq3/VqQkbJbpSgo+Gkbgs4DbUToBQCg7N1000067bTTdPXVV+vkk0/u1rHGjh2rtWvXtth2ww036OKLL26xAuWxxx6rt99+W4sXL5Ykffe739Xs2bNbDQIvBcJnZ5KpppWO2hMW9v5M17hq8Bv0dN3TnZY1jRntWNCzAwCAQrr44ot13nnn6c4772we79KdbnfXdTVixIgW2/74xz9q1qxZ6tevnyTpxRdf1OLFi/Xkk09qv/32kyTdeOONmjZtml5++WVNmDChe2+qmwifZcZzPY2oGtGl0e5+tv17bwEA6HWslfz60pw7jyW358+fr6uvvloPPPCADjnkkHbLLViwQAsWLOjwWA8++KBmzJjRavvKlSu1evVqXX311c3bVqxYoerq6ubgKUlTp05VdXW1li9fTvhEa57rda2gw4ogAIA+xK+XFowqzbn/413Jq+py8QcffFB/+tOftHTpUh188MEdlp0zZ46OOeaYDsuMHj26ze0333yzJk6cqOnTpzdvq6ur07Bhw1qVHTZsWKuFf0qB8FkiDX6ojzKBGoJtLZyN2UgNpuvze2b8bSuCbN1Sr8yHmyRJCccolXBzrpNJuHLSHU9sDwAAOrfHHnto/fr1Ou+887TPPvuof/92xo9IqqmpUU1NTc7naGho0J133qn/9//+X6vX2lrC3Frb7tLmcSJ8lshz72zUK+9tUTZsbN72v+9sVFWi6+EzbNiiTxbU/Pv/LJdNNrWYppKuPju8vzw3t8kMnKpKVU2bRgAFAJSnZGVTC2Spzp2D0aNH65577tFBBx2kmTNnavHixe0G0Hy73e+++27V19e3WNxHalra/L333mt1jPfff1/Dhw/P6X0UA+GzRBKuo36ppDLhtpHyFQlXFcmu/5ME/rbWTS+dku3XX0EYqSGMZKuq5OTQ+mmzWUVb62UDVlYCAJQpY3Lq+i61HXbYQcuWLdNBBx2k2tpaPfTQQxowYECrcvl2u99888066qijNHTo0Bbbp02bpk2bNunpp5/WvvvuK0l66qmntGnTphbd86VC+CyRtOuqwnPlhNsCouMEchxfCeN1qVncTWxr2fQqK6SqSimI5GcDOekKOcmuh89ITQEUAAAUzpgxY/Too4+2CKCflk+3+z//+U899thjWrRoUavXJk6cqJkzZ+qUU07R9ddfL6lpqqUvfelLJR9sJDHJfFm5+R8X6VcvzNfvX7+yw1Wk2uJHfpFqBQAAumP06NFatmyZNm7cqEMPPbQgk73fcsstGj16tGpra9t8/Y477tDuu++u2tpa1dbWao899tBtt93W7fMWAi2fJeY5ae3Uf5Je2/JC87Z3619XYLNKmlSH+xptax19e+urGtN/gKTcBxoBAIDCaWv1opEjR+qll14q2Dk6u0+0pqZGt99+e8HOV0iEzxIzxugHky7R3z54TAmT0M3/uKjL+yadbVMyBZEvayMRPgEAQDkjfJYBY4ySjqeE4Z8DAAD0btzzCQAAgNjQ1FYIfkPXymW3TWPkBlvl+E2jy52wXiZslNm+5TNslLHRp4/QgomyLZ6bsEEmCuSEoRRkcp6TDAAAoNgIn91hXCm7VXr9sa6VDySpQpJUs/YxDUg1NTw3RFn1a3xV7nb3a/Zf/7/yOumGN5ltc4Smt7ylfglXTuQqHUZyvMHSDntLbseDlgAAAOJE+OyOREoaMFqyXZyY3bGSmsoGyX4Kk03hM4wyCv0Kme3CZ+SmFZlkh4cz7rbplSLHU5RIy0YJ2SgjBY1SFOY8/sg2NujT74ZlNwEAQKEQPrsrkUvLopW0RZIUuRWKPp4kPoocWTepyG67BTdyPEVOZ+Fz21RL1nEVuSk5JqHIyWOVItdV1FCvrU8+1eollt0EAACFQviEJMnxPCVGjJTClsGVZTcBAEAhET57ET8KJGsV2UD5zPfpeF6rbSy7CQAAConw2Uv4NtAbmXcURJIbhBofjRZDjQAAKD8nnniiNm7cqPvuu6/UVSkJ5vnsJfq5lUoZTwk58hUoUm5rwwMAgMI48cQTdfTRR7fYdvfddyudTuviiy/WFVdc0eYSnF1VV1en2bNna8SIEaqqqtLee++tu+++u0WZDRs2aPbs2aqurlZ1dbVmz55dkDXlC4Hw2UskjKukk5BrWF4TAIByctNNN+m4447TVVddpR/96Eeqrq7WwIED8z7e7Nmz9fLLL+v+++/Xc889p69+9auaNWuWVq1a1Vzm2GOP1erVq7V48WItXrxYq1ev1uzZswvwbrqP8AkAAFAkF198sU4//XTdeeedOvnkkyW13TKaixUrVuiMM87Qvvvuq5122kk/+clPNHDgQP3tb3+TJL344otavHixbrrpJk2bNk3Tpk3TjTfeqAceeEAvv/xyId5Wt3DPJwAA6BGstWoIuriqYIFVJCpkjOm84Hbmz5+vq6++Wg888IAOOeSQdsstWLBACxYs6PBYDz74oGbMmCFJOuCAA3TXXXfpyCOP1MCBA/Vf//VfymQyOvDAAyU1hdPq6mrtt99+zftPnTpV1dXVWr58uSZMmJDT+yg0wmeZ8m3QNNS8A07kd1wAAIBepCFo0H537td5wSJ46tinVJnDstUPPvig/vSnP2np0qU6+OCDOyw7Z84cHXPMMR2WGT16dPPzu+66S7NmzdLgwYOVSCRUWVmpP/7xj/rMZz4jqeme0GHDhrU6xrBhw1RXV9fl91AshM8yddl7t3daJpW1ui2GugAAgNzsscceWr9+vc477zzts88+6t+/f7tla2pqVFNT0+Vj/+QnP9GGDRv0P//zPxoyZIjuu+8+/du//Zsef/xx7b777pLUZiuttTbn1ttiIHyWkaRJaKw3XG9l3yt1VQAAKDsViQo9dWzrlfjiOncuRo8erXvuuUcHHXSQZs6cqcWLF7cbQHPpdn/11Vd11VVX6f/+7//0uc99TpK055576vHHH9fVV1+t6667TiNGjNB777XOEu+//76GDx+e0/soBsJnGTHG6ITBX27qcu+EbwNd/RbtngCAvsMYk1PXd6ntsMMOWrZsmQ466CDV1tbqoYce0oABA1qVy6Xbvb6+XpLkOC3HjLuuqyhqul9v2rRp2rRpk55++mntu+++kqSnnnpKmzZt0vTp07v9vrqL8FlmjDHyTMdrukvq9H7QbBSoMWiU4yTkua1XLgIAAMU3ZswYPfrooy0C6Kfl0u2+6667auedd9b3vvc9XXrppRo8eLDuu+8+LVmyRA888IAkaeLEiZo5c6ZOOeUUXX/99ZKk7373u/rSl75U8sFGEuGz1zHGKBsFermxTsn1/6dUeoAm1UzqdgC1jQ3qSau7m4QrJ50udTUAANDo0aObW0APPfRQjRo1Ku9jJZNJLVq0SPPnz9eXv/xlffTRR9p55531m9/8RkcccURzuTvuuENnnnmmamtrJUlHHXWUrrrqqm6/l0IgfPYyCbnq71SqwhjJSSgTNirqTmx0XUUN9dr6ZGnuscmXU1WpqmnTCKAAgNi1tXrRyJEj9dJLLxXk+LvssovuueeeDsvU1NTo9ts7H7xcCoTPzgQZyW9s//VESiqDkWPbSxhXnuPIOgkFnfXPd8LxPCVGjJTCntPuabNZRVvrZYOeU2cAAPoKwmdn7j2l49eHTZJm/qLsAmghOV7Pumc0UlMABQAA5YflNduSrJTG7NO1suteaGodBQAAQKdo+WyLMdK37pVe/LOU6icl2pjWIWiU/utb8dcNAACgByN8tseYpvs5E2kp2UMHrUS+FIVStkFK2OKcw3ElN1WcYwMAgF6H8NkLWWNkwozsh280BdBMRirWXJ+JtDRqLwIoAADoEsJnL2RNQkGyQqGy8o1VxknLdGXi+lyFgdRYryiTlZLbLiXXkTzXLfz5AABAj0f47CVMEMr4gdzIqjEbqs4PFChQNvL16kf18pwuhM9EIqdR+06YkRs2akO0UZG7bTqqlOdq99HVBFAAANAK4bOXGPebpd0+RjRqhPxZR3U5gBrjypWjiqSrKNF0KflhpEw2VBhJInsCAIBPYaqlHiyTlF4aU7jjOe/WKQobFTlBlx7WCRQ5oRw3VCrhKJVwlHS5pAAA6MiJJ56oo48+utTVKBlaPnsyY3Tet1yd//ZUVaQqJa9l13pgQ20J65U0Hf8zO36oz/36MUnSPxvfkg271mRpolBOlFG4tZ/GDthdya507QMA0MudeOKJ2rhxo+67777mbXfffbe+9a1v6YILLtAVV1wha/OfhebVV1/V2WefrSeeeEKZTEYzZ87UlVdeqeHDhzeX2bBhg84880zdf//9kprWdr/yyis1cODAvM9bKDRT9XTG6PyxT+namlWKkq5sMtH8cL2U+qUHyPMqOnwkvW1TSaVMUinjdfGRlGtcZcOsrO3eMp4A+rYwGyrbGLT7CLIsl4ue66abbtJxxx2nq666Sj/60Y9UXV2ddwjcunWramtrZYzRI488or/+9a/KZrP68pe/rCja9ll87LHHavXq1Vq8eLEWL16s1atXa/bs2QV6R91Dy2cPlTQJjfWG663se5KkNxIb5WdCeZ/6J00YV+rkFk6z3VeQpJOQdbp2WRgrhXx/AVAA697Yog/e2dru68m0q7G71ijhcTM5epaLL75Y5513nu6880597Wtfk9R2y2hX/fWvf9Ubb7yhVatWacCAAZKkX//616qpqdEjjzyiQw45RC+++KIWL16sJ598Uvvtt58k6cYbb9S0adP08ssva8KECQV7f/kgfPZQxhidMPjLamzYrEs3/lepq1OWbGOD+kJbiUm4ctI9dCEE4GOO68iraPsjKfQj+Y2hoqhIi2Wgx7DWyjY0lOTcpqJCJocZYSRp/vz5uvrqq/XAAw/okEMOabfcggULtGDBgg6P9eCDD2rGjBnKZDIyxiiV2ja/djqdluM4euKJJ3TIIYdoxYoVqq6ubg6ekjR16lRVV1dr+fLlhE/kzxjT6f2cfZLrKmqo19Ynnyp1TWLhVFWqato0Aih6NNdzlOygVTMMuLUHkm1o0Mt7Ty7JuSf8baVMZRvLbbfjwQcf1J/+9CctXbpUBx98cIdl58yZo2OOOabDMqNHj5bUFCKrqqr04x//WAsWLJC1Vj/+8Y8VRZHWrl0rSaqrq9OwYcNaHWPYsGGqq6vr8nsoFpILeh3H85QYMVIKe3+7p81mFW2tlw16/3sFgJ5kjz320Pr163Xeeedpn332Uf/+/dstW1NTo5qami4dd+jQofrDH/6g73//+/rVr34lx3H0zW9+U3vvvbfc7ebXbquV1lqbc+ttMRA+0Ss5XpGWEy0zkZoCKAD0BaaiQhP+trJk587F6NGjdc899+iggw7SzJkztXjx4nYDaC7d7pJUW1urV199VevXr1cikdDAgQM1YsQIjR8/XpI0YsQIvffee62O8f7777cYEV8qhE8AANAjGGNy6voutR122EHLli3TQQcdpNraWj300EPNg4S2l0u3+/aGDBkiSXrkkUe0bt06HXXUUZKkadOmadOmTXr66ae17777SpKeeuopbdq0SdOnT+/u2+o2wicAAECRjBkzRo8++miLAPppuXS7S02j2ydOnKihQ4dqxYoV+sEPfqCzzjqreSDRxIkTNXPmTJ1yyim6/vrrJUnf/e539aUvfankg40k5vkEAAAoqtGjR2vZsmXauHGjDj30UG3cuLFbx3v55Zd19NFHa+LEibrgggt07rnn6tJLL21R5o477tDuu++u2tpa1dbWao899tBtt93WrfMWCi2fAAAABXTrrbe22jZy5Ei99NJLBTn+z3/+c/385z/vsExNTY1uv/32gpyv0AifKAjfZqVQyoaR/ChQY9Ag4xRuMmhHrjy3bwwiAgCgNyN8oluMjPwoqzUf/UOSFIRW2TCUu76fksnCTeeQctOaVDOJAAoAQA9H+ES3JIyrfslqWScpSXJsJGsjpRMV8tzChE8/CpQJGxX1ifWKAADo3Qif6LaEk5B1mlokrY0URqE811PKLdx4tiD0C3YsAABQOox2BwAAQGwInwAAAIgN3e69SPbjeyKTcmVU+rVbAQAAPo3w2YtcWLVUkjQuHKQ5jVMJoAAAoOzQ7d7DJZXQjpmW68S+6W6Qz8hwAABQhmj57OGMMTrtg8+rscJVxjPNrZ95H89vHVptwpUMragAABTCiSeeqI0bN+q+++4rdVVKgvDZCxgZeUrIFuBYO970UKttjSNr9O7X9yeAljHb2EBbN3qcMLPtqo0aGhU6bV/FUTZUlAkVflSvMCjcymkof2EmIxtFsmEoG/acv3Lf/s53tHHjRv3x3nubt919zz2affzx+tn55+vyyy6TtTbv93TDjTfqd7/7nf62apW2bNmiD9ev18CBA1uU2bBhg34wd67u//OfJUlHHXWUrrzyyhbl1qxZo9NOO02PPPKIKioqdOyxx+rSSy+V5xV3QRfCJ2QTrhpH1ii99sM2X0+v/VAmCGWTbV8uTuQr+uR5GMkJIyloVLfu6jCu5Cbz37+vcF1FDfXa+uRTpa4JkLMgkKQKSVLj88/LtvN5F4SSHxht/SBUlk+tPiXrurLDhynculWh33Pme458XzYIFG7ZIkm6+be/1Zk//KF+dcklOun445vLffJ6rj7asEG1Bx6o2gMP1LkXXKDwo48Uui2/mB37jW/o7Xff1X/fc4+cigp9b84czZ49W3/+OIyGYagjjzxSQ4cO1RNPPKEPPvhAJ5xwgqy1uvLKK/N8513DrzEkY/Tu1/eXCVp+AzN+2GZL6Pb7mSij9OY3mzeFUSQvtDJOWurOJPOuJw3emQDaCcfzlBgxUupBLQLAJ0xgpY+/ujrplJxU238znMDK8SWnysgp4LK9KH+OY5o+a1xXxu05rd7GbKv3JZdfrvMvuki333yzvvqVr0iSvjNnjjZu2qR7f/e7vI4/9/TTJUmPPv540/kcp8X/nxdfflkPLV2qv/7PUu07eW+5/fvrxhtv1LRp0/Tyyy9rwoQJevjhh/XCCy/orbfe0qhRoyRJv/zlL3XiiSfqoosu0oABA1qfuEAIn2hiTLstm+2xjqswWS3ZbR3+kSJFNpISKSnf5TXDUAqzkg0lET474xS5ewQoFse3krZKkkzKk2knfBrXyhjJqUgQPvsYI318y5eRjCNrrYJs1MlexZHwnKZQ2SVNdT7np+fr2htv0J/+cLcOOeigVq/LNF3zCy+5RD//5aUdHvGBe+7VjP33/9RpPq6PcZqPJUlPPv2Mqqurtd+++zR37U+dOlXV1dVavny5JkyYoBUrVmi33XZrDp6SdNhhhymTyWjlypU6qEV9C4vwiW6xTstvopFsU/h0vfzDp7JN/WyfkgmzeR5PcuTKcwlpANCTBdlIN5/7bEnOfdJFU5RMdb31dfGSh3X/fz+ghx94QAd/4cAOy37vpJP0b1/9aodlRm8XEjtT9957GjZkaKvtw4YNU11dXVOZujoNHz68xeuDBg2S53nNZYqF8Imy5xhHmSCjlz58Me9jpNy0JtVMIoACAGKxx267af0HH+j8Cy/SPntPVv/+/dstW1NTo5qamoKev61WWmtti+1dKVMMhE+UvaSTVE26pqlFNQ9+FCgTNipiPDgA9GgJz9FJF00p2blzMWrkKP3X7XfokCOP0BH/+q9a9Mc/thtA8+52b8eI4cP13vvrWm1///33m1s7R4wYoaeeajlYdcOGDfJ9v1WLaKERPtEjJJ3u3fsZhD1nlCQAoG3GmJy6vktth7Fj9ciDi3XIkUdo5le+ogfvu6/NgTyF7nafuu++2rRpk55+9lnt8/nPS5Keeuopbdq0SdOnT5ckTZs2TRdddJHWrl2rkSNHSpIefvhhpVIpTZ48ucvnygfhEwAAoEjGjB6tpYse1CFHbAugn5Zrt3vde++p7r339Oqrr0mSnnv+efXv3187jBmjmpoaTdx1Vx126KGac+aZuvqyy+RWVup73/++vvSlL2nChAmSpNraWk2aNEmzZ8/WJZdcog8//FBnn322TjnllKKOdJdYXrMwgkbJ78IjaFSFmh7ZbKOymU8eGUXZQFHWV5jjI8j6CqJIQRgq2H66HVuIKecBAEB3jR41So8sXqxNmzbqsK8cpY2bNnXreNfffJOm7D9d3zujacqlg2Yepin7T9efFy1qLnPbTTdrt0mTdMTXvqbDDj9ce+yxh2677bbm113X1X//938rnU5r//331zHHHKOjjz5al17acfd/IRhryz+lbN68WdXV1dq0aVPR03izzBbplYekVH8pWdn6db9RuvPr8dSli+qN0X47jpUkLX47q9UDv9KtVYmMH2j8tU0X8uvfP6JLUzEFoZUfRhpTUykv76mWslKQlYZOkBLp/I6xnUyYUUPYoD2H7qm0W9Ht4wEoDN+3uuP2pqmWDjzIU6qdqZaCwMrPSmN3SCjJVEt9SkbS2pSnHcfuoHQqVerq9Dy2aXUot3//gs2T2tjYqNdff13jx49XOt3yM7qreY2Wz3wlUtKwSaWuRbsGBe/JVVDqagAAALTAPZ/5Mkaa+QspyOS0W70f6Kk3Nqki6SqVaMr+DVFW/1v/iiqMJ8/JbSog05hV/5feVphOakvSlfRwTvsDAADEifDZHcZIydy6hY1CeV6jPC8h7+PwGUZGTpCQY5JycxzV7URWiY+X1XKdbQ3ZDcYoq1COEjKimwoAAJQHwmcvdeC4MZIe0bhwkOY0TiWAAgCAssA9n71IUq5s/dgW2950N8hncnUAAFAmCJ+9iJFR+OZ39NQbb+nRN98udXUAAOiWT/rsyn5anj6kEJMkET57HaNKa1VR/jNoAQDQIUeSrJUfsEpduaivr5ckJZP5rzzIPZ8AAKAsuZJSYaj1H36ohJuQ0425q/sm2zTPZ2Njt+f5tNaqvr5e69at08CBA+V243iETwAAUJaMpCFhpHfr67Xm7bdKXZ2ex1rJWplUSsYpTGf3wIEDNWLEiG4dg/CJPiMTZgt2LEeuPDe3OVkBALlLStohCOUHIfd+5shmGhVtrVflvvvKrWpjtcYcJZPJbrV4foLwiV7PMY4yQUYvffhiwY6ZctOaVDOJAAoAMTCS+GubuyiyisJQ6VRKbrr7y1UXCuETRRGGkbL5jmcLrRRGskEkKSpAbVxVJQYWZISeJPmRr4/8Bm3NZhUmCrNWbl/gOpJXoLWFAQA9F+ETXWL8rs0V6kZWjdlQdevzX1feiXw5ka+PMhsU5bjc6CeiZKJpBaoiCKwvP2pU2LBRSaehKOfojVKeq91HVxNAAaCPI3yiS3a86aESnPXxvPfcOmak/nncV4sSQN0olDGu0smEPJdfoa7ww0iZbKgwUtPwVQBAn8UnJ9plE64aR9YovfbDUlclZ1Vvr5UXRbJe/vOQtccYR5ExSiUceS5T5XZVEBbiFgoAQE9H+ET7jNG7X99fJoh3eU4TBnKirOqrd5J1Ujnt6/i+Jlx+U5FqBgAAuovwiY4ZI5uM+TJxJBuFsl5S1smt5ZK2NQAAyht9hgAAAIgN4RMAAACxodu9D8iquPdsJuXKiPV2AQBA5wiffcCFVUuLevxx4SDNaZxKAAUAAJ2i272XqrBWO4YDYznXm+4G+UVuXQUAAL0DLZ+9lJF0asMUNZrifb/IKix6qyoAAOhdCJ+9mJGR14P/iZ3Iz3nqJBP52z3PSFHhJ18ykS8TZmXCBhnL5E5d4QSRnDCU/HrltcSR40pubnO+AgDKU89NJui9jJGJMkpvfjP3Xbdbg75y02uyycKv5ejbQI71lWqsV9Js+xVyjFGSX6k2+WGkdBjJcaqkRB6t8Ym0NGovAigA9AJ8UqLsWMdVmKyWrM15X+MEzc8jx5N1inGJJ9QYhno13NRiq2cS2ik1skUgRZPQRAoVSsmK3MNn6EtBoxSFrAsPAL0An5JlxLeBFEmOcfp8gLFOnilju15w6yZk3cL/f3SUUH8nIbtdOA5sqIxCBW5Srin8evI9nZVVZMOmFsx8Wj6joPMyAIAeoW8nnDLhylHapNRoMwpsqMYwq8GJ6j4fQMtZwrhqMbNUJIWWEf8AAHQmr6HQ11xzjcaPH690Oq3Jkyfr8ccf77D8HXfcoT333FOVlZUaOXKkvv3tb+uDDz7Iq8K9kecktXvlztq7aqI+V/EZpR1PEQNZAABAL5Rz+Lzrrrs0d+5cnXvuuVq1apVmzJihww8/XGvWrGmz/BNPPKHjjz9eJ510kp5//nn94Q9/0DPPPKOTTz6525XvTTwnqQonpZTjlboqAAAARZNz+Lzssst00kkn6eSTT9bEiRN1+eWXa+zYsbr22mvbLP/kk09qxx131Jlnnqnx48frgAMO0Pe+9z09++yz3a48AAAAepacwmc2m9XKlStVW1vbYnttba2WL1/e5j7Tp0/X22+/rUWLFslaq/fee0933323jjzyyHbPk8lktHnz5hYPAAAA9Hw5hc/169crDEMNHz68xfbhw4errq6uzX2mT5+uO+64Q7NmzZLneRoxYoQGDhyoK6+8st3zLFy4UNXV1c2PsWPH5lJNAAAAlKm8BhwZY1r8bK1tte0TL7zwgs4880ydd955WrlypRYvXqzXX39dc+bMaff455xzjjZt2tT8eOutt/KpJgAAAMpMTnP5DBkyRK7rtmrlXLduXavW0E8sXLhQ+++/v374wx9KkvbYYw9VVVVpxowZuvDCCzVy5MhW+6RSKaVSrGQCAADQ2+QUPj3P0+TJk7VkyRL967/+a/P2JUuW6Ctf+Uqb+9TX1yuRaHka122aQNzmsYIN0FXbL7VZ9HNFgYwNJddvOf9nHJIJqZ2eBwAAyk3Os5jPmzdPs2fP1pQpUzRt2jTdcMMNWrNmTXM3+jnnnKN33nlHv/3tbyVJX/7yl3XKKafo2muv1WGHHaa1a9dq7ty52nfffTVq1KjCvhtgOzve9FCpqxALO3qE7De/QgAFAPQIOYfPWbNm6YMPPtAFF1ygtWvXarfddtOiRYs0btw4SdLatWtbzPl54oknasuWLbrqqqv07//+7xo4cKAOPvhg/eIXvyjcuwA+ZhOuGkfWKL32w1JXJTbmnTpZP5A8lvUEAJS/vNZvPPXUU3Xqqae2+dqtt97aatsZZ5yhM844I59ToYfIqrhd3Em5Ml3pzzZG7359f5kg3qUu/ShQxvraOT1WqbjWdvcDOdf8Jp5zAQBQICwejoK4sGppUY8/LhykOY1TuxxAbTLeS9tGkrVRU+tjXOETAIAeiPBZpnwbSF1Y3j1hg+JXph1JuRoXDtKb7oain+tNd4N8hfK4ZAEA6NH4JC8zrhylTUqNNqPAdt517AdbVV3kLu/2GBnNaZwqv4jnzyoseqtqIflR0OHsua4cJYwbX4UAACgzhM8y4zlJ7V65s8IuNHtmoqxeanxJkbVt5h03hlbRUAl5hsvIGKNsFOiNzDsdlvOMp53SowmgAIA+i9RQhjynMPcMHrq++INRPkyO0IqBTPOTMK4GuFUdzl0b2FBZm1WoSAkRPgEAfVNey2uifDUopWeiz8Z2vhq/Tq5Kd99pOUkYV0kn0e6D1k4AAGj57IWM/i37U10x+lGlnC6MWMqTa4NYWlYBAEDvQvjslYxCk1Roihc+AQAA8kG3OwAAAGJD+AQAAEBsCJ8AAACIDeETAAAAsSF8AgAAIDaETwAAAMSG8AkAAIDYED4BAAAQG8InAAAAYkP4BAAAQGwInwAAAIgNa7ujx8gqzHvfpFwZmQLWJn9+FBTma5/1VfHx04z1JVuAY+bBlaOEcUtzcgBAj0P4RI9xYdXSvPcdFw7SnMapJQ2gxhhlo0BvZN4pzPH8ULt//PyfjW/JhqUJgJ7xtFN6NAEUANAlhE+UtaRcjQsH6U13Q7eO86a7Qb5CeSW85BPG1QC3StYWponSmKD5ecokZU387y2wobI2q1CREiJ8AgA6R/jsBUIbKbL5d0mXMyOjOY1T5efZ5Z5V2K0W00JLGFeFanw123XdJ52ErFOCX+dICnvptQcAKA7CZw/mylHKeAoVqiHcqmTCK3WVisLIlLTFEgAAFA6f6D2Y5yT1ucrPKPDq9Yb5QFGBunPR8xi/NK2PJgpkbCi5fsctupGVglA26yuKchxtFfiS70uZjBTRtd+bRMG2v1k2m5U1bV9ENrBSVooyoaKoPAYOAj1BlMkoymQKdrtXoRA+ezjPJCSHD+S+bsebHip1FTqU+vixrltHua0gdUH5CB1P+pf/r+mHhWcpiLLtlk1Iei+eagG9Tr9/+Repf79SV6MZ83wCPZRNuGocWVPqagAAkBNaPoGeyhi9+/X9ZYLSDfjxo0AZ62vn9FilTLKDclaZINTOw/rLS+Ta7d4o+Q3S6MlSsrKbNUY58QMr/dfHszac8/8pkWq7Sz0IrIKsNHpsQskk3e5AV0UNDYo++kgmnS51VVogfAI9mTGyydL9GttIsjaSvKTUQfhUaCXjyHhJObmGTyeUTCClUlIy1b0Ko6w4jpXUFD6N58l4bV8bxmm6X81JJeQQPoGuiyLJ92XauZ+6VAifAGLjh1HuOwWRFESK/FDqxipXKD/BdgOOskGk9tYpCAIrP5Qa/FBBmaxUhm1cR/Jcxh6g6wifAIrOGCkbRnpj/dac93XCjNywURuijYrcxiLUDqUShZI+XiT2tfe3KtFO43kUSlFotD4KVYrpbNGxlOdq99HVBFB0Gb/GAIou4RgNSCeVz2wfxrhy5agi6SpK8CerNwm3a8hOJdz2w6cjhUZKe0Yul0BZ8cNImWyoMJJY5Axdxa8xgFgknPy6Sx05cqwjL+HI5nq/KMra9jdRJFyjZDvhJZTkWCmVIHyWoyCf22nQp/GXHAAAALEhfAIAACA2hE8AAADEhvAJAACA2BA+AQAAEBvGDQLoNj8KOv0q68pRor1ZxAEAfQbhE0DejDHKRoHeyLzTaVnPeNopPZoACgB9HOGzl8paVyri1Guu3dbMlYkchWbbz56JVGbLyEqSsgVcmjEpV4Zl/pQwrga4VbKdzB4f2FBZm1WoSAlmogaAPo3w2UvNf2daUY9foUZ9JX1j87kalG5+bafUJs0b9veyC6AXVi0t2LHGhYM0p3EqAVRNAbTT/w2RFFrWZQcAMOCoV3Hka7y3sdTV0GuZamVteVxaSbkaFw4q+HHfdDfIL2BLKgAAfQUtn72AI0ee8dRoQn1vyJPaGln1dyuLem+da33pg6bnPx+9QqFJKmvdore45srIaE7j1IIFxazCgragAgDQ1xA+e4GEcbVjeqRCLyH/48EfngmVzHMt7a5w7bYbSlNOpNBERb3HtDuMjDwudQAAygKfyL1EIrBKGCvZSE42kDGOTK49346RTTAYBAAAFA/hs6dzXSmdkhozUhDKWF9uoy/XSK7JrSnSyfryB/UjgAIAgKIhfPZ0XlKauIsUNt3TGNqsNmw1ShtPKcfr8mGcbKCqV96Roo6nzAFKxQmzHd/ZYVxZt+vXPACgNAifvYGXlJRseh4ZhWFSkfEU5RA+gXJljSMTNqjqwxc6LBclUtpasxsBFADKHOETQFmzTkJBeoiM7aDdM/LlBBmJuUQBoOwRPgGUPesk1NENIY4kRUFMtQEAdEd5zAQOAACAPoGWT3Sba4OP/+uoQo0fP/dbzAWaq1AJld36nAAAoNsIn+i2Q9f/pvn5J+u9f7L6Ub4+TI7QioFfIYACANDL0O2OvIRK6MPkiKIdv8avkyvu4QMAoLeh5RP5MUYrBn6lRUDMRE7z2u4/H71CKSf3bnfXBi1aUgEAQO9C+ET+jFH4yfyikkLjqEHpj58nm9Z7BwAA2A7hE0Bs/Cgoys0+js3Ksb6yUUY2bPsExjhKOsk2XwMAxIfwCaDojDHKRoHeyLxTnONHoZwoo4YtXrsrHCWdlMb1+ywBFABKjPAJoOgSxtUAt0rWdjRVfP6MCeSYSJGblm1jWdnABvKjjGw3pv8CABQG4RNALBLGlYo0c5axkmNC+U6yzfCpSAotsycAQDkgfAJ5yqrlOuJJuTLFSlcAAPQShM9eyreBFEmOcZQ0/DMXw4VVS1v8PC4cpDmNUwmgAAB0gEnmexlXjtImpUChGmxGHwSbmoIoCiIpV+PCQW2+9qa7Qf6nWkMBAEBLNIn1Mp6T1O6VOytUpEyU1fMNryqyUdHutetrjIzmNE5tETKzClu1ggIAgLYRPnshj6lkisrIyONXBwCAvNDtDgAAgNgQPgEAABAbwicAAABiQ/gEAABAbAifAAAAiA1DdtGCEwRqtfq1Y2QTbimqAwAAehnCJyRJ1nEUppJyM75cP9viNSfryx/UjwAKAAC6jfAJSZL1Eto6YYxM1LLd08kGqnrlHSmyJaoZAADoTQifaGa9hIiYAACgmBhwBAAAgNgQPgEAABAbwicAAABiwz2fAPoM32alsHjHN8ZR0kkW7wQA0AsQPgH0esYY+WFGaz76R1HPk3RSGtfvswRQAOgA4RNAr+eapPp7A2Vt8eZzCGwgP8rI2lbLNAAAtkP47AN8G6j1skVd40RZpW2gMHLU3meqMUYJwwT0KG+uSUqmiCeIpNAGRTwBAPQOhM9ezJWjtEmp0WYU2PxudHOtrwrrK7RS1E76zEaBBrhVYvyalO3khsKkXJmiJiAAAMob4bMX85ykdq/cWWG+zZ6S5Gbkpq1sRUryvFYv+1GgNzLvFLU7sye5sGpph6+PCwdpTuNUAigAoM8ifPZyXncHPhgrmeS2x6fR2KmkXI0LB+lNd0OnZd90N8hXKI9fPQBAH8UnINBNRkZzGqfK76DLPauw01ZRAAD6AsInUABGhtZMAAC6gE5TAAAAxIbwCQAAgNgQPgEAABAbwicAAABiQ/gEAABAbAifAAAAiA3hEwAAALEhfAIAACA2hE8AAADEhvAJAACA2LAeILrGD9rebn052UDGODLbfZcx2UDGiXI+jbHbzmOygYwxOR9DkuQY2YSb374AAKBoCJ/omOtK6ZTUmJGCsNXLxvpyG325pml98+bdGn25pnX5Tk9n/W3PG7KSsXlV28n68gf1I4ACAFBmCJ/omJeUJu4ihW0HydBmtWGrUdp4kk1JTzZt3zJprLJu7sHRDbPS6k+OsYNC18v5GE42UNUr70hRfsEV6A7fZqXcv3eVjDGOkk6y1NUA0IcQPtE5LympnQ+nyCgMk4qMJ2O3BUWb9hTlET63byyN0p6iPMInUArGGPlhRms++kepq5KTpJPSuH6fJYACiA3hE2UrEWbz2s8Js03d99bI2vJo/XS3awpzbSBX5VGvniJUQsr3/t+YuCap/t7AsrnmuiKwgfwoI2tzvz8bAPKVV/i85pprdMkll2jt2rX63Oc+p8svv1wzZsxot3wmk9EFF1yg22+/XXV1dRozZozOPfdcfec738m74uj9vvLML7p3gPyya1HUG6Nz+42VJB2y/jeq7EEBpRx8mByhFQO/0iMCqMq7ii1FUmjbGUwIAEWSc/i86667NHfuXF1zzTXaf//9df311+vwww/XCy+8oB122KHNfY455hi99957uvnmm7Xzzjtr3bp1CgL+4PUmvg1ko22fupkoK5vXYCGrdf3HatiWtwpXOfR4NX6dXAUK27v9AwDQY+QcPi+77DKddNJJOvnkkyVJl19+uR566CFde+21WrhwYavyixcv1rJly/Taa6+ppqZGkrTjjjt2r9YoG64cpU1KjTaj7HateQ02qzDPrrz7Jn1Dichvsa0x8lWTGKCk6fySdRqz6v/CGoUVnqxXHneWZBVKekSStGjot+TZ3jMKPymnxUwHheTaQIeu/01Rjg0AKI2cPpmz2axWrlyp+fPnt9heW1ur5cuXt7nP/fffrylTpujiiy/WbbfdpqqqKh111FH6z//8T1VUVLS5TyaTUSaTaf558+bNuVQTMfKcpHav3FmhIjVu15i9V+UEpQuU+zJRVs83vCrfJOQ4nQ9Asq4UmqRCk5TtQliNQ7hdOPtZ1bIS1qTwxoWDNKdxatECKACgd8npk3n9+vUKw1DDhw9vsX348OGqq6trc5/XXntNTzzxhNLptP74xz9q/fr1OvXUU/Xhhx/qlltuaXOfhQsX6mc/+1kuVUMJeR+PkjXbrZdV4aSUZv2sZkm5GhcO0pvuhlJXpeDedDfIVyiP8YsAgC7I69Pi06vOWGvbXYkmiiIZY3THHXeourpaUlPX/de//nVdffXVbbZ+nnPOOZo3b17zz5s3b9bYsWPzqSpQFoyM5jROld+TJoDsRFahLqxaWupqAAB6mJzC55AhQ+S6bqtWznXr1rVqDf3EyJEjNXr06ObgKUkTJ06UtVZvv/22dtlll1b7pFIppVKpXKoGlD0jQ+sgAKDPy6lj1PM8TZ48WUuWLGmxfcmSJZo+fXqb++y///5699139dFHHzVve+WVV+Q4jsaMGZNHlQEAANBT5XxX3rx583TTTTfplltu0YsvvqizzjpLa9as0Zw5cyQ1dZkff/zxzeWPPfZYDR48WN/+9rf1wgsv6LHHHtMPf/hDfec732l3wFFZCTKSX9/+I8h0fgwAAABIyuOez1mzZumDDz7QBRdcoLVr12q33XbTokWLNG7cOEnS2rVrtWbNmuby/fr105IlS3TGGWdoypQpGjx4sI455hhdeOGFhXsXxeAkJK9Kym6VOlppJ7tVGjBaSnCbAAAAQGfyugHt1FNP1amnntrma7feemurbbvuumurrvqyl6yQdjpQijqYDN9vkF5/TLK9ZxAJAABAMTH6oSPJHnBbAAAAQA/CTIwAAACIDeETAAAAsSF8AgAAIDaETwAAAMSG8AkAAIDYED4BAAAQG6ZaQlE0FnDq08ZIyoaOXGNkrOm0vAmNPOsqjBzZiO9XxZLd7rtrJnJkO/ku65lIpvN/PgBAL0f4RFEc93Ahj5aStEeO++wqbS5kHdCKyar/roslSfPfmS5Zr8PiO6U2ad6wvxNAAaCPo1kIBZNypUmDSl0LlKvXMtXKWv7kAEBfR8snCsYY6eL9pUyBVxttiDJaXf+yKoynlNNx65okmWygyn+8IzfTwdKovYiT9eUPrJJNuLGeN6tAF338/Oejl8tr589J1rqa/860+CoGAChrhE8UlDFSusBXlY0kz43kGSvPsZ3vUOEq2HWUwigqbEXKkJMNVPXKO3IUyjrx9mcbbfv/m3IieWrn/3fv/2cAAOSA8IleyXoJdSGmApDk26zUhR4LYxwlnWTxKwSgVyN8osfwbZB3K5pjHCUNl3uxZDtILlk5kslKltBSbowx8sOM1nz0jy6VTzopjev3WQIogG7h0xhlz5WjtEmp0WYU2PxuKG0MsxqcqCaAFsmFVUs7fL3/rosV1I+TDXeLqUboCtck1d8bKGs77ycIbCA/ysha7qMA0D18EqPseU5Su1furDDPZs9MlNXzDa8qspHEND8Fk5SrceEgvelu6FL5ROWb8rdMVJpJNsqKa5Jd+72IpND2jUF8AIqL8IkewaObr+wYGc1pnCq/k5sFt1irS/otialWAIByR/gEkDcj0+4US5/w6KYFAGyH/i8AAADEhvAJAACA2BA+AQAAEBvCJwAAAGJD+AQAAEBsCJ8AAACIDeETAAAAsSF8AgAAIDZMMl8IQabrZf1QbtAgx0nIKWL2t8aVdb2iHR8AACAfhM/ucBKSVyVlt0phtmv7+JESwVYljCvXFi98OsFWZStHEUABAEBZIXx2R7JC2ulAKQq6vEuYCfRe+L76pRJKe25RquUGDRpQ96SMDWWLcgYAAID8ED67K1mR4w6BwmS9wmRSUbI44RMAAKBcMeAIAAAAsSF8AgAAIDaETwAAAMSG8AkAAIDYED4BAAAQG8InAAAAYsNUSwBikzWhsoq6XN5VqHpjmvZVqFCmzXJGgRyFysqXbWfxhqQSMqbt/QEA8SF8AojNJf2W5LzPuf3GfvzskS6UfkLtrawwViP0Hf0rARQASozwiT7Dt4FyaHTrkGMcJQ2/Pl2RlKugfpwSlW+WtB5vqU6+AnlKlrQeANDX8emJXs+Vo7RJqdFmFNiwIMdsDLManKgmgHaBkVHDm3Mk4+vno5cr5eTQ7W4DHbL+N5Kk/xlygsJ2/n+bMJATZVU/cCdZk2rxWla+LrW35l1/tOTbrFSYXyOFkSRVSpKC0G/3uFEohZFRNgzl0nBdVrJhJD8K1Bg0yDis2pcLR6481yt1NUqCT070ep6T1O6VOyssULNnJsrq+YZXFdlI7dyCiFaMZD15SsjL6Z5Pq0rb1I/uyVXYzp8sI8lRqEBJWfOpls12uuGRG2OM/DCjNR/9o2DHtKEjabIkqa7xLTlh2+nTho4UOvrwo00yboG6L1AQQWibvhSs76dkkj+IuUi5aU2qmdQnAyjhE32C59DVCnSHa5Lq7w2UtYVL89YaffTx86RJyjFtt5xZY2SNK89Jy8mh5RzF59hI1kZKJyrk0SzdZX4UKBM2KipUN0IPQ/gEegEnCAp1O2vBmWjb6HOTDWRyCA/GBi33bWewkAkDmSiQk8nKOi3LOPL1yW2eTjYrh6bQvBV6br5ou8/dZCg57WQXGzWV9XyJ7FleTCjZ0MrzrVKW8NllYaQwyipqaFRUxLsVbDZbvIN3A+ET6MGs4yhMJeVmfLl+ef6Rce22v6xuoy/XdP2bvmv9bc8bspJpOziaKJBjfbkNjbKfSieugubw6TY0Nv2MstDyi0lWTtT2tRFFjpzQldvQIIdu97IShVaJMJKpT0oJwmdXmTArJ8zIprYqcov7N8mpqpRJlNf9uIRPoAezXkJbJ4yRicr3AzkTGunJpudbJo1V1u16y6MbZqXVn+y7g8J27o1ywqycsFFbhkyUdStavJaNstLrTc837jlJntP37q8qVzaQ9GLT863jd5CT6qCcb5QdM0qfvqUXpZUNIjUEgcaMGSTXK6+AU878oFH1/keqGH2AKhNVRT2XSbhy0uminiNXhE+gh7Neoqw7km24rTXEpj1FOYTP7RtJo7SnqL0b80NJYaioIi3rtvwja6Ptzl+RUtRewkHstmvYlvU82fbCpytZI0XppAzfHcpKGEQK/UC2skIifHaZDSTrB3KqquQmixs+yxHLawIAACA2tHwCAEov/Lh7vQ02kGwoWcvsZkBvQPgEAJRc9tHOO+ISYyKJuyaAHo9udwBAaSQkZ1g537EMoBho+QQAlIQxUmqmVfCGlUlK7a1W2zQ4Q3xiAb0Ev8oAgJIxRjJuU/BsL3xKkqKmsgB6PrrdAQAAEBvCJwAAAGJD+AQAAEBsCJ8AAACIDeETAAAAsSF8AgAAIDaETwAAAMSG8AkAAIDYMMk8gD7Dj7JFPX7CeDLMhA4AHSJ8AnnybSBF2352jKNkh0u0oNSufem8oh5/VOV4fWP8GQRQAOgA3e5Ajlw5SpuUAoVqsJnmxwfBpqZAirKSMJ5GVY6P5Vzv1r+uwBa3dRUAejqaaYAceU5Su1furHC7Zs9MlNXzDa8qspFEo1dZMcboG+PPKGoo9KNs0VtVAaC3IHwCefCcZKmrgBwYY5Q0qVJXAwAgut0BAAAQI8InAAAAYkP4BAAAQGwInwAAAIgN4RMAAACxIXwCAAAgNoRPAAAAxIZ5PnsxE2aL/u3CGlfW9Yp8FgAA0FsQPnsjk1CYrJTr10tRcZf6c4KtylaOIoACAIAuIXz2QlEirc0j9peKvM64GzRoQN2TMjaULeqZAABAb0H47KWiRLrUVQAAAGiFAUcAAACIDeETAAAAsSF8AgAAIDaETwAAAMSG8AkAAIDYED4BAAAQG8InAAAAYsM8n0AB+TaQovZfd4yjpOHXDgDQd/EpCBSAK0dpk1KjzSiwYbvlGsOsBieqCaAAgD6LT0CgADwnqd0rd1bYQbNnJsrq+YZXFdlIMjFWDrHyo2ypq5CThPFkDBckgPgQPoEC8ZxkqauAMnDtS+eVugo5GVU5Xt8YfwYBFEBsGHAEAN2UMJ5GVY4vdTXy8m796wpsz2qtBdCz0fIJIDbZ0FGHI7I+xQ23tcZlQ6OwnfsVTGjkhI6ygWSt7W418/LVsaf3qBDnR1nd/M+fSpKygZV1SvP/zQZSEEomlNprfLVh0yMKJEOTSVnJBlZ+YEv2e4eeifAJIDYXPL5jTuUr1Khvprft26B0F/YqZQDsQV3Xxqj/rk1PL7nfl2wPqPv/lroCaM9ep0iVqVLXAj0F3yEBFFXSsdqxuqHU1QAAlAlaPgEUlTHS9ye/Kz/KvWXNDbPSM03Pz5vxhkLXa/scYVZO2KiPhuwh63aldRR+ZHXdK03Pf3hUUkmn7f+3xWazUvCWZDypvRnIbNBULjG2qRzKRyaI1OgH8kgTyAGXC4CiM0by3NzvCXO1bR/PtQrbOYYjK0eRsgnJuj2g+7gMmO2+DHgJo6RTmv9vNpIcVzIfP9osYyXrSolE+wEVpWFlFFrDbAnICd3uAAAAiA3hEwAAALEhfAIAACA2hE8AAADEhvAJAACA2BA+AQAAEBvCJwAAAGJD+AQAAEBsCJ8AAACITV7h85prrtH48eOVTqc1efJkPf74413a769//asSiYT22muvfE4LAACAHi7n8HnXXXdp7ty5Ovfcc7Vq1SrNmDFDhx9+uNasWdPhfps2bdLxxx+vL37xi3lXFgAAAD1bzqvkXnbZZTrppJN08sknS5Iuv/xyPfTQQ7r22mu1cOHCdvf73ve+p2OPPVau6+q+++7Lu8JAT+fbQIpKXYvCcYyjJAtuAwC6KKdPjGw2q5UrV2r+/PktttfW1mr58uXt7vfrX/9ar776qm6//XZdeOGFnZ4nk8kok8k0/7x58+ZcqgmUJVeO0ialRptRYMNSV6dgGsOsBieqCaAAgC7J6dNi/fr1CsNQw4cPb7F9+PDhqqura3Off/zjH5o/f74ef/xxJRJdO93ChQv1s5/9LJeqAWXPc5LavXJnhb2o2TMTZfV8w6uKbCSZUtcGANAT5NVUYUzLTxlrbattkhSGoY499lj97Gc/02c/+9kuH/+cc87RvHnzmn/evHmzxo4dm09VgbLiOclSVwEAgJLKKXwOGTJEruu2auVct25dq9ZQSdqyZYueffZZrVq1SqeffrokKYoiWWuVSCT08MMP6+CDD261XyqVUiqVyqVqAAAA6AFyGu3ueZ4mT56sJUuWtNi+ZMkSTZ8+vVX5AQMG6LnnntPq1aubH3PmzNGECRO0evVq7bffft2rPQAAAHqUnLvd582bp9mzZ2vKlCmaNm2abrjhBq1Zs0Zz5syR1NRl/s477+i3v/2tHMfRbrvt1mL/YcOGKZ1Ot9oOAACA3i/n8Dlr1ix98MEHuuCCC7R27VrttttuWrRokcaNGydJWrt2badzfgIAAKBvymvA0amnnqpTTz21zdduvfXWDvc9//zzdf755+dzWgAAAPRwrO0OAACA2BA+AQAAEBuWJAGAPs6Psl0qlzBem3M6A0AuCJ8A0Mdd+9J5XSo3qnK8vjH+DAIogG6h2x0A+qCE8TSqcnxO+7xb/7oC27VWUgBoDy2fANAHGWP0jfFndClM+lG2y62jANAZwicA9FHGGCUNSxkDiBfd7gAAAIgN4RMAAACxodsd3WbCbMG/xVjjyrpegY8KAABKjfCJ/JmEwmSlXL9e6uI8gV3lBFuVrRxFAAUAoJchfCJvUSKtzSP2l2xQ0OO6QYMG1D0pY0PZgh4ZAACUGuET3RIl0qWuAgAA6EEInwC6zbeBFHVcxjGOkoY/OQDQ1/FJACBvrhylTUqNNqPAhh2WbQyzGpyoJoACQB/HpwCAvHlOUrtX7qywk2bPTJTV8w2vKrKRxLLgANCnET4BdIvnJEtdBQBAD8Ik8wAAAIgN4RMAAACxIXwCAAAgNoRPAAAAxIbwCQAAgNgQPgEAABAbwicAAABiQ/gEAABAbAifAAAAiA3hEwAAALEhfAIAACA2hE8AAADEJlHqCvRVmSAs2rFdx8hL8L0CAACUH8JnzBKOUWUqofpMID+MinKOrdlAo6orCaAAAKDsED5jlk66OmDnIQoiW5TjN2RDrXh1vcIiHR8AAKA7CJ8lkE66pa4CAABASRA+AQBd5kfZgh7PRlJgJfPxo80y9uNHJJni3K2EPPlRJD8KlQkzcgIaVroqE2aUDbOytm/2UhI+AQBddu1L5xXnwJkulPlncU6N7vvt66WuQc90wOgD1E/9Sl2N2BE+AcTGt4GUQ8tVYrtWtmyUVWDaLudEvpzIVzbMyKrtQsY4SjpeLtXFxxLG06jK8Xq3noQBoPsInwCKzpWjtEmp0WYU2K5PM5aw28Jng80oaKeLylhfrs2oIayXbSfdZsIGDUoNJYDmwRijb4w/Q4EtbJe7JNmsFLwlGU8y7Xwi2aCpXGJsUzmUj0wQqcEPtdfYQapgPEOXNYQN2upvVdpNl7oqJUH4BFB0npPU7pU7K8yl2VOSE27ri92raldFbqrtckGDXH+rNg7aR1GidRdWNmrUSxtXyVpuGMyXMUZJ0/b//+6wjmTMtkebZUzTI+FIhhnkykrkRAqcQCk3pVSC8NlVkSL5kS/T3kXfyxE+AcTCc5I572O2a+lMO56s0074dCK5jq9Gt1JRorJ1gSDnUwMAioTvkAAAAIgN4RMAAACxIXwCAAAgNoRPAAAAxIbwCQAAgNgQPgEAABAbploC0CM4YabdWUJNmJEJM00LgAMAyhrhE0CPsOuy73da5pnR/6JI/WOoDQAgX3S7Ayhb1klp68DPlroaAIACouUTQPkyRm9M+alMlOmwmBM0yvU/UuRWxFQxAEC+CJ8Aypsxsm66wyLWWtnIb39xcABA2aDbHQAAALEhfAIAACA2dLujbJkwy7cjSJKscWVdr9TVAAAUAOET5cckFCYr5fr1UpQtdW1QBpxgq7KVowigANALED5RdqJEWptH7C/ZoNRVQRlwgwYNqHtSxoZiCnkA6PkInyhLUaLj0c0AAKBn4pY6AAAAxIbwCQAAgNjQ7Q6gz/CtL4XFO74xjpIOg6IAoCOETwC9nmNcpdwKZcIGBfKLdp5M2KBBqaEEUADoAOETQK/nuSlNHDhZkS1es2c2atRLG1fJ2qho5wCA3oDwCaBP8NxUcU/AzGAA0CUMOAIAAEBsCJ8AAACIDeETAAAAsSF8AgAAIDaETwAAAMSG8AkAAIDYED4BAAAQG8InAAAAYkP4BAAAQGwInwAAAIgNy2v2UpmgeGtYd8Z1jLwE32sAAEBrhM9eJuEYVaYSqs8E8sOoJHXYmg00qrqSAAoAAFohfPYy6aSrA3YeoiCyJTl/QzbUilfXKyzR+QEAQHkjfPZC6aRb6ioAAAC0iX5RAAAAxIbwCQAAgNgQPgEAABAbwicAAABiw4AjACgg3/pS6abZzZkxjpKOV+pqAOhDCJ8AUACOcZVyK5QJGxTIL3V1uiwTNmhQaigBFEBsCJ8AUACem9LEgZMV2Z7T7JmNGvXSxlWytjQLUgDomwifAFAgnpsqdRVyE5S6AgD6IgYcAQAAIDaETwAAAMSG8AkAAIDYED4BAAAQG8InAAAAYkP4BAAAQGwInwAAAIgN4RMAAACxIXwCAAAgNoRPAAAAxIbwCQAAgNgQPgEAABAbwicAAABiQ/gEAABAbAifAAAAiA3hEwAAALEhfAIAACA2eYXPa665RuPHj1c6ndbkyZP1+OOPt1v23nvv1aGHHqqhQ4dqwIABmjZtmh566KG8KwwAAICeK+fwedddd2nu3Lk699xztWrVKs2YMUOHH3641qxZ02b5xx57TIceeqgWLVqklStX6qCDDtKXv/xlrVq1qtuVBwAAQM+Sc/i87LLLdNJJJ+nkk0/WxIkTdfnll2vs2LG69tpr2yx/+eWX60c/+pH22Wcf7bLLLlqwYIF22WUX/fnPf+525QEAANCz5BQ+s9msVq5cqdra2hbba2trtXz58i4dI4oibdmyRTU1Ne2WyWQy2rx5c4sHAAAAer6cwuf69esVhqGGDx/eYvvw4cNVV1fXpWP88pe/1NatW3XMMce0W2bhwoWqrq5ufowdOzaXagIAAKBM5TXgyBjT4mdrbattbfnd736n888/X3fddZeGDRvWbrlzzjlHmzZtan689dZb+VQTAAAAZSaRS+EhQ4bIdd1WrZzr1q1r1Rr6aXfddZdOOukk/eEPf9AhhxzSYdlUKqVUKpVL1QAAANAD5NTy6XmeJk+erCVLlrTYvmTJEk2fPr3d/X73u9/pxBNP1J133qkjjzwyv5oCAIrCt76yYWOnDz/KlrqqAHqBnFo+JWnevHmaPXu2pkyZomnTpumGG27QmjVrNGfOHElNXebvvPOOfvvb30pqCp7HH3+8rrjiCk2dOrW51bSiokLV1dUFfCsAgFw4xlXKrVAmbFAgv9PymbBBg1JDlXS8GGoHoLfKOXzOmjVLH3zwgS644AKtXbtWu+22mxYtWqRx48ZJktauXdtizs/rr79eQRDotNNO02mnnda8/YQTTtCtt97a/XcAAMiL56Y0ceBkRTbstGw2atRLG1fJ2iiGmgHozXIOn5J06qmn6tRTT23ztU8HykcffTSfUwAAYuC5Xby/PihuPQD0HaztDgAAgNgQPgEAABAbwicAAABiQ/gEAABAbAifAAAAiA3hEwAAALEhfAIAACA2hE8AAADEJq9J5oHOZILOV0xB3+I6Rl6C77sA0NcRPlFQCceoMpVQfSaQH7IMH7bZmg00qrqSAAoAfRzhEwWVTro6YOchCiJb6qqgjDRkQ614db1CrgsA6PMInyi4dNItdRUAAECZov8LAAAAsSF8AgAAIDaETwAAAMSG8AkAAIDYED4BAAAQG8InAAAAYkP4BAAAQGwInwAAAIgN4RMAAACxYYUjAECX+daXwsIdz0ZSFLlSZGWi9ssoMoqiUKaA50b3ZcNI2ShUY9AoObRndVU2zJa6CiVF+AQAdMoxrlJuhTJhgwL5hTtwaGRsUoqs1E74VCTJGtnQlxxbuHOj27JRpEwYaKufVGAIn7moTFTKNX1zOWrCJwCgU56b0sSBkxXZwjY9htlIWzZ8JNdz5CTbDi+RHynMRuo/qJ9cj4BTThr8UFszvqaPHqZ+KSJFLlzjKp1Il7oaJcGVAgDoEs9NFfyYYRQq4/pyXVeO2074jCKFbqh0okJuom+2FJUrG4UKXF9VySpVJYkU6Bq+QgIAACA2hE8AAADEhvAJAACA2BA+AQAAEBvCJwAAAGJD+AQAAEBsCJ8AAACIDeETAAAAsSF8AgAAIDaETwAAAMSG8AkAAIDYED4BAAAQG8InAAAAYkP4BAAAQGwInwAAAIgN4RMAAACxIXwCAAAgNoRPAAAAxIbwCQAAgNgQPgEAABAbwicAAABiQ/gEAABAbAifAAAAiA3hEwAAALEhfAIAACA2iVJXAEDfkQnCnPdx/VDZIFJgIkWK2izjBJESQaTGbKjQ5n4OlE6Ybfr3dRzJMbbNMlEQKQoiNWRDuTHXDx3L53caIHwCKLqEY1SZSqg+E8gP2w6Q7XH9QFV+qMAGCqOg7TJBoEQQ6qNMoDDyC1FlxCTKRmr0Qxlj5dj2w6f1I5lMIMfmdv2g+CpTCSUcU+pqoAchfAIounTS1QE7D1EQtR0uOpTZItetllL9pGRl22X8lJRJaNedh0qp/t2rLGLlNwZ6O3CVTLtKeG23awbZUH5jqDETapRM87FVbhKOUTpJmzS6jt9iALHI/8MpISUdKek2PdrkSpEjpRJND/QYWSulEo68pKtkO/++vpWcwKoqlZDHvy/Q4zHgCAAAALEhfAIAACA2hE8AAADEhvAJAACA2BA+AQAAEBvCJwAAAGJD+AQAAEBsCJ8AAACIDeETAAAAsSF8AgAAIDaETwAAAMSG8AkAAIDYED4BAAAQG8InAAAAYkP4BAAAQGwInwAAAIgN4RMAAACxIXwCAAAgNoRPAAAAxIbwCQAAgNgQPgEAABAbwicAAABiQ/gEAABAbAifAAAAiA3hEwAAALEhfAIAACA2hE8AAADEhvAJAACA2BA+AQAAEBvCJwAAAGJD+AQAAEBsCJ8AAACIDeETAAAAsSF8AgAAIDaETwAAAMSG8AkAAIDYED4BAAAQG8InAAAAYkP4BAAAQGwInwAAAIgN4RMAAACxSZS6AgDQJUEmv9cAAGWF8AmgvDkJyauSslulMNt+Oa+qqSwAoKzxlxpAeUtWSDsdKEVBx+WcRFNZAEBZI3wCKH+ESgDoNRhwBAAAgNgQPgEAABAbwicAAABik1f4vOaaazR+/Hil02lNnjxZjz/+eIflly1bpsmTJyudTmunnXbSddddl1dlAQC9U+hH8rNhm4/Qj0pdPQAFlHP4vOuuuzR37lyde+65WrVqlWbMmKHDDz9ca9asabP866+/riOOOEIzZszQqlWr9B//8R8688wzdc8993S78gCAns1xjJJpV2EQKdsQtPkIg0jJtCvHMaWuLoACMNZam8sO++23n/bee29de+21zdsmTpyoo48+WgsXLmxV/sc//rHuv/9+vfjii83b5syZo7///e9asWJFl865efNmVVdXa9OmTRowYEAu1QUAlLkgGyqKOv4ochyjhOfGVCMA+ehqXsup5TObzWrlypWqra1tsb22tlbLly9vc58VK1a0Kn/YYYfp2Wefle/7be6TyWS0efPmFg8AQO+U8Fx56USHD4In0HvkFD7Xr1+vMAw1fPjwFtuHDx+uurq6Nvepq6trs3wQBFq/fn2b+yxcuFDV1dXNj7Fjx+ZSTQAAAJSpvAYcGdPyvhtrbattnZVva/snzjnnHG3atKn58dZbb+VTTQAAAJSZnFY4GjJkiFzXbdXKuW7dulatm58YMWJEm+UTiYQGDx7c5j6pVEqpVCqXqgEAAKAHyKnl0/M8TZ48WUuWLGmxfcmSJZo+fXqb+0ybNq1V+YcfflhTpkxRMpnMsboAAADoyXLudp83b55uuukm3XLLLXrxxRd11llnac2aNZozZ46kpi7z448/vrn8nDlz9Oabb2revHl68cUXdcstt+jmm2/W2WefXbh3AQAAgB4hp253SZo1a5Y++OADXXDBBVq7dq122203LVq0SOPGjZMkrV27tsWcn+PHj9eiRYt01lln6eqrr9aoUaP0q1/9Sl/72tcK9y4AAADQI+Q8z2cpMM8nAABAeSvKPJ8AAABAdxA+AQAAEBvCJwAAAGJD+AQAAEBsCJ8AAACIDeETAAAAsSF8AgAAIDaETwAAAMSG8AkAAIDYED4BAAAQG8InAAAAYkP4BAAAQGwInwAAAIgN4RMAAACxIXwCAAAgNoRPAAAAxIbwCQAAgNgkSl2BrrDWSpI2b95c4poAAACgLZ/ktE9yW3t6RPjcsmWLJGns2LElrgkAAAA6smXLFlVXV7f7urGdxdMyEEWR3n33XfXv31/GmKKfb/PmzRo7dqzeeustDRgwoOjnQ8/DNYKOcH2gM1wj6EhPvT6stdqyZYtGjRolx2n/zs4e0fLpOI7GjBkT+3kHDBjQo/7RET+uEXSE6wOd4RpBR3ri9dFRi+cnGHAEAACA2BA+AQAAEBvCZxtSqZR++tOfKpVKlboqKFNcI+gI1wc6wzWCjvT266NHDDgCAABA70DLJwAAAGJD+AQAAEBsCJ8AAACIDeETAAAAsSF8fso111yj8ePHK51Oa/LkyXr88cdLXSWUyMKFC7XPPvuof//+GjZsmI4++mi9/PLLLcpYa3X++edr1KhRqqio0IEHHqjnn3++RDVGKS1cuFDGGM2dO7d5G9cH3nnnHX3rW9/S4MGDVVlZqb322ksrV65sfp1rpO8KgkA/+clPNH78eFVUVGinnXbSBRdcoCiKmsv01uuD8Lmdu+66S3PnztW5556rVatWacaMGTr88MO1Zs2aUlcNJbBs2TKddtppevLJJ7VkyRIFQaDa2lpt3bq1uczFF1+syy67TFdddZWeeeYZjRgxQoceeqi2bNlSwpojbs8884xuuOEG7bHHHi22c330bRs2bND++++vZDKpBx98UC+88IJ++ctfauDAgc1luEb6rl/84he67rrrdNVVV+nFF1/UxRdfrEsuuURXXnllc5lee31YNNt3333tnDlzWmzbdddd7fz580tUI5STdevWWUl22bJl1lproyiyI0aMsD//+c+byzQ2Ntrq6mp73XXXlaqaiNmWLVvsLrvsYpcsWWK/8IUv2B/84AfWWq4PWPvjH//YHnDAAe2+zjXStx155JH2O9/5TottX/3qV+23vvUta23vvj5o+fxYNpvVypUrVVtb22J7bW2tli9fXqJaoZxs2rRJklRTUyNJev3111VXV9fimkmlUvrCF77ANdOHnHbaaTryyCN1yCGHtNjO9YH7779fU6ZM0b/9279p2LBh+vznP68bb7yx+XWukb7tgAMO0NKlS/XKK69Ikv7+97/riSee0BFHHCGpd18fiVJXoFysX79eYRhq+PDhLbYPHz5cdXV1JaoVyoW1VvPmzdMBBxyg3XbbTZKar4u2rpk333wz9joifr///e/1t7/9Tc8880yr17g+8Nprr+naa6/VvHnz9B//8R96+umndeaZZyqVSun444/nGunjfvzjH2vTpk3adddd5bquwjDURRddpG9+85uSevffEMLnpxhjWvxsrW21DX3P6aefrv/93//VE0880eo1rpm+6a233tIPfvADPfzww0qn0+2W4/rou6Io0pQpU7RgwQJJ0uc//3k9//zzuvbaa3X88cc3l+Ma6Zvuuusu3X777brzzjv1uc99TqtXr9bcuXM1atQonXDCCc3leuP1Qbf7x4YMGSLXdVu1cq5bt67Vtw70LWeccYbuv/9+/eUvf9GYMWOat48YMUKSuGb6qJUrV2rdunWaPHmyEomEEomEli1bpl/96ldKJBLN1wDXR981cuRITZo0qcW2iRMnNg9i5W9I3/bDH/5Q8+fP1ze+8Q3tvvvumj17ts466ywtXLhQUu++PgifH/M8T5MnT9aSJUtabF+yZImmT59eolqhlKy1Ov3003XvvffqkUce0fjx41u8Pn78eI0YMaLFNZPNZrVs2TKumT7gi1/8op577jmtXr26+TFlyhQdd9xxWr16tXbaaSeujz5u//33bzU92yuvvKJx48ZJ4m9IX1dfXy/HaRnDXNdtnmqpV18fJRzsVHZ+//vf22QyaW+++Wb7wgsv2Llz59qqqir7xhtvlLpqKIHvf//7trq62j766KN27dq1zY/6+vrmMj//+c9tdXW1vffee+1zzz1nv/nNb9qRI0fazZs3l7DmKJXtR7tby/XR1z399NM2kUjYiy66yP7jH/+wd9xxh62srLS33357cxmukb7rhBNOsKNHj7YPPPCAff311+29995rhwwZYn/0ox81l+mt1wfh81OuvvpqO27cOOt5nt17772bp9VB3yOpzcevf/3r5jJRFNmf/vSndsSIETaVStl/+Zd/sc8991zpKo2S+nT45PrAn//8Z7vbbrvZVCpld911V3vDDTe0eJ1rpO/avHmz/cEPfmB32GEHm06n7U477WTPPfdcm8lkmsv01uvDWGttKVteAQAA0HdwzycAAABiQ/gEAABAbAifAAAAiA3hEwAAALEhfAIAACA2hE8AAADEhvAJAACA2BA+AaADjz76qIwx2rhxY1HPc+utt2rgwIHNP59//vnaa6+9inpOACgFwicAbOfAAw/U3Llzm3+ePn261q5dq+rq6ljrcfbZZ2vp0qWxnhMA4pAodQUAoJx5nqcRI0bEft5+/fqpX79+sZ8XAIqNlk8A+NiJJ56oZcuW6YorrpAxRsYY3XrrrS263T/pHn/ggQc0YcIEVVZW6utf/7q2bt2q3/zmN9pxxx01aNAgnXHGGQrDsPnY2WxWP/rRjzR69GhVVVVpv/3206OPPtpuXT7d7X7iiSfq6KOP1qWXXqqRI0dq8ODBOu200+T7ft7nAIBSoOUTAD52xRVX6JVXXtFuu+2mCy64QJL0/PPPtypXX1+vX/3qV/r973+vLVu26Ktf/aq++tWvauDAgVq0aJFee+01fe1rX9MBBxygWbNmSZK+/e1v64033tDvf/97jRo1Sn/84x81c+ZMPffcc9pll126VL+//OUvGjlypP7yl7/on//8p2bNmqW99tpLp5xySsHOAQDFRvgEgI9VV1fL8zxVVlY2d7W/9NJLrcr5vq9rr71Wn/nMZyRJX//613XbbbfpvffeU79+/TRp0iQddNBB+stf/qJZs2bp1Vdf1e9+9zu9/fbbGjVqlKSmezoXL16sX//611qwYEGX6jdo0CBdddVVcl1Xu+66q4488kgtXbpUp5xySsHOAQDFRvgEgBxVVlY2B09JGj58uHbccccW92gOHz5c69atkyT97W9/k7VWn/3sZ1scJ5PJaPDgwV0+7+c+9zm5rtv888iRI/Xcc88V9BwAUGyETwDIUTKZbPGzMabNbVEUSZKiKJLrulq5cmWL8Cgpp0FFcZwDAIqN8AkA2/E8r8VAoUL4/Oc/rzAMtW7dOs2YMaOgx47zHABQCIx2B4Dt7Ljjjnrqqaf0xhtvaP369c0ti93x2c9+Vscdd5yOP/543XvvvXr99df1zDPP6Be/+IUWLVpUgFrHcw4AKATCJwBs5+yzz5brupo0aZKGDh2qNWvWFOS4v/71r3X88cfr3//93zVhwgQdddRReuqppzR27NiCHD+ucwBAdxlrrS11JQAAANA30PIJAACA2BA+AQAAEBvCJwAAAGJD+AQAAEBsCJ8AAACIDeETAAAAsSF8AgAAIDaETwAAAMSG8AkAAIDYED4BAAAQG8InAAAAYkP4BAAAQGz+f4FOXBBJHyEWAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots(figsize=(8,8))\n",
    "by_ki = {}\n",
    "for ki, df in cleaned.groupby('ki'):\n",
    "    by_ki[ki] = df\n",
    "    km_ki = km.fit(df['time'], df['status'])\n",
    "    km_ki.plot(label='Ki=%s' % ki, ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88b42065",
   "metadata": {},
   "source": [
    "## Exercise 11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "id": "6232b3b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let first create the table 11.4 in csv format and import it to do the assignments\n",
    "df=pd.read_csv('table11_4.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "id": "5c599a6c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 6 entries, 0 to 5\n",
      "Data columns (total 3 columns):\n",
      " #   Column  Non-Null Count  Dtype  \n",
      "---  ------  --------------  -----  \n",
      " 0   Y       6 non-null      float64\n",
      " 1   delta   6 non-null      int64  \n",
      " 2   X       6 non-null      float64\n",
      "dtypes: float64(2), int64(1)\n",
      "memory usage: 276.0 bytes\n"
     ]
    }
   ],
   "source": [
    "df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab3ee852",
   "metadata": {},
   "source": [
    "### (a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "id": "30402938",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['group']=np.where(df['X']<2,1,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "id": "34a7c6cc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp8AAAKnCAYAAADAwSbIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1eklEQVR4nO3dfZTXdZ3//8eHueRCIEEHSFTURBTThM0L9JtWUubR9eR+Zdc28qJOhKnIakl2Vpct0XRdLQPzmxdduOVvNz1Hiy2pVdPMUpL9uopLrSiWEF/YYkh0YGY+vz/M2UZAGS5enxm43c75nDPz/rzfn/dz5t3Ivff7c1GpVqvVAABAAf1qPQAAALsO8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFBMfa0H2BKdnZ158cUXs9tuu6VSqdR6HAAAXqdarWbt2rUZNWpU+vXb/PnNPhGfL774YkaPHl3rMQAAeBMvvPBC9tprr83e3yfic7fddkvy6g8zePDgGk8DAMDrtba2ZvTo0V3dtjl9Ij5fu9Q+ePBg8QkA0Iu92VMkveAIAIBixCcAAMWITwAAiukTz/kEAHYN1Wo17e3t6ejoqPUovE5dXV3q6+u3+W0vxScA0CusX78+y5cvz7p162o9CpsxYMCAjBw5Mo2NjVv9GOITAKi5zs7OLF26NHV1dRk1alQaGxt9sEwvUq1Ws379+vy///f/snTp0rztbW97wzeSfyPiEwCoufXr16ezszOjR4/OgAEDaj0Om9C/f/80NDTk+eefz/r169Pc3LxVj+MFRwBAr7G1Z9MoY3scH0cYAIBixCcAAMWITwCAbbRixYpceOGFOeCAA9Lc3JyWlpYce+yxuemmm3r1q/dvvvnmHH/88Rk8eHAqlUp+//vf7/B9esERAMA2ePbZZzNp0qQMHTo0V155ZQ499NC0t7dnyZIlufXWWzNq1Kiceuqpm9x2w4YNaWhoKDzx/1i3bl3e//735/3vf39mzZpVZJ/OfAIAbIPp06envr4+jz/+eM4444yMGzcuhx56aE4//fR873vfyymnnNK1bqVSyU033ZQ///M/z8CBA/O5z30uSTJv3rzsv//+aWxszNixY/ONb3yja5vnnnsulUolixYt6lr2+9//PpVKJQ888ECS5IEHHkilUsn3vve9HHbYYWlubs6RRx6ZJ5988g1nnzFjRi699NIcddRR2+8X8iac+QQAeqVqtZqXN9Tmk476N9Rt0fuMrl69Ovfdd1+uvPLKDBw4cJPrvP5xLr/88syZMyf/+I//mLq6utx999258MILc/311+e9731vvvvd7+bss8/OXnvtlRNOOKFHc19yySW54YYbMmLEiHzmM5/JqaeemiVLltT07OrriU8AoFd6eUNHDv7bH9Rk30/Pfl8GNL55Jv3qV79KtVrN2LFjuy0fPnx4XnnllSTJeeedl6uvvrrrvjPPPDPnnHNOt+/POuusTJ8+PUkyc+bMPProo7n22mt7HJ+XX355TjzxxCTJ1772tey11165++67c8YZZ/TocXYkl90BALbR689u/vznP8+iRYtyyCGHpK2trdt9EydO7Pb94sWLM2nSpG7LJk2alMWLF/d4jqOPPrrr69133z1jx47dqsfZkZz5BAB6pf4NdXl69vtqtu8tccABB6RSqeSZZ57ptny//fZ79XH6999om01dnn99vFar1a5lr72xe7Va7bp/w4YNWzTfph671pz5BAB6pUqlkgGN9TW5bWmwDRs2LCeeeGJuvPHGvPTSS1v1c44bNy4PP/xwt2WPPPJIxo0blyTZY489kiTLly/vuv9PX3z0px599NGur3/3u99lyZIlOeigg7Zqrh3FmU8AgG0wd+7cTJo0KRMnTswVV1yRt7/97enXr18ee+yxPPPMM5kwYcIbbn/JJZfkjDPOyBFHHJH3vOc9uffee3PXXXflhz/8YZJXz54eddRRueqqq7Lvvvtm1apV+exnP7vJx5o9e3aGDRuWlpaWXHbZZRk+fHhOO+20ze57xYoVWbFiRX71q18lSZ588snstttu2XvvvbP77rtv3S/kzVT7gDVr1lSTVNesWVPrUQCAHeDll1+uPv3009WXX3651qNslRdffLH6yU9+sjpmzJhqQ0NDddCgQdV3vvOd1Wuuuab60ksvda2XpHr33XdvtP3cuXOr++23X7WhoaF64IEHVr/+9a93u//pp5+uHnXUUdX+/ftXDz/88Op9991XTVK9//77q9VqtXr//fdXk1Tvvffe6iGHHFJtbGys/tmf/Vl10aJFbzj35ZdfXk2y0e22227b5PpvdJy2tNcqf/xFbLEf//jHueaaa7Jw4cIsX748d9999xsWdZI8+OCDmTlzZp566qmMGjUqn/rUpzJt2rQt3mdra2uGDBmSNWvWZPDgwT0ZFwDoA1555ZUsXbo0Y8aMSXNzc63H6XMeeOCBnHDCCfnd736XoUOH7rD9vNFx2tJe6/FzPl966aUcdthhufHGG7do/aVLl+YDH/hAjjvuuDzxxBP5zGc+kwsuuCDf+c53erprAAD6uB4/5/Okk07KSSedtMXr33TTTdl7771z/fXXJ3n1SbWPP/54rr322px++uk93T0AAH3YDn/B0U9/+tNMnjy527L3ve99ueWWW2r+eaabUu3szMvr1tZ6DNjIln7axk6vYUDi9wDQzfHHH58ePpOyZnZ4fK5YsSItLS3dlrW0tKS9vT2rVq3KyJEjN9qmra2t2xuytra27ugxu7y8bm0GXLt3sf0BPfTWiclHfyhAAfqoIu/zuak3Tt3U8tfMmTMnQ4YM6bqNHj16h88I9BG/eTxZ85taTwHAVtrhZz5HjBiRFStWdFu2cuXK1NfXZ9iwYZvcZtasWZk5c2bX962trcUCtP+A3bLu4mVp/d3qtLe9UmSf8EZebu/Mqbf/Z5JkwUfHZa/dd9FXgW54OZl3zKtfd66v7SwAbLUdHp9HH3107r333m7L7rvvvkycOHGzz/dsampKU1PTjh5tkyr9+mXAoCEZMGhITfYPr7dufXtezvNJko6h+yS7b/yxbLuE9Vv3ySEA9C49vuz+hz/8IYsWLer6WKelS5dm0aJFWbZsWZJXz1pOnTq1a/1p06bl+eefz8yZM7N48eLceuutueWWW3LxxRdvn58AAIA+o8dnPh9//PGccMIJXd+/dnn8Ix/5SG6//fYsX768K0STZMyYMZk/f34uuuiifPnLX86oUaPyxS9+0dssAQDsgnocn2/2Uv7bb799o2Xvete78otf/KKnuwIAYCdT5NXuAAA7sxUrVuTCCy/MAQcckObm5rS0tOTYY4/NTTfdlHXr1tV6vE367//+75x//vkZO3ZsBgwYkL333jsXXHBB1qxZs0P3u8NfcAQAsDN79tlnM2nSpAwdOjRXXnllDj300LS3t2fJkiW59dZbM2rUqJx66qmb3LaWH7jz4osv5sUXX8y1116bgw8+OM8//3ymTZuWF198Mf/yL/+yw/brzCcAwDaYPn166uvr8/jjj+eMM87IuHHjcuihh+b000/P9773vZxyyild61Yqldx000358z//8wwcODCf+9znkiTz5s3L/vvvn8bGxowdOzbf+MY3urZ57rnnUqlUul7snSS///3vU6lU8sADDyRJHnjggVQqlXzve9/LYYcdlubm5hx55JF58sknNzv3+PHj853vfCennHJK9t9//7z73e/O5z//+dx7771pb2/fvr+kP+HMJwDQO1WryYYaXbLewo/yXb16de67775ceeWVGThw02+F9/oP1bn88sszZ86c/OM//mPq6upy991358ILL8z111+f9773vfnud7+bs88+O3vttVe3F3lviUsuuSQ33HBDRowYkc985jM59dRTs2TJki0+u7pmzZoMHjw49fU7LhHFJwDQO21Yl1w5qjb7/syLSeObv6/yr371q1Sr1YwdO7bb8uHDh+eVV179sJrzzjsvV199ddd9Z555Zs4555xu35911lmZPn16klffSejRRx/Ntdde2+P4vPzyy3PiiScmSb72ta9lr732yt13350zzjjjTbddvXp1/v7v/z4f//jHe7TPnnLZHQBgG73+7ObPf/7zLFq0KIccckja2tq63Tdx4sRu3y9evDiTJk3qtmzSpElZvHhxj+c4+uiju77efffdM3bs2C16nNbW1px88sk5+OCDc/nll/d4vz3hzCcA0Ds1DHj1DGSt9r0FDjjggFQqlTzzzDPdlu+3335Jkv79+2+0zaYuz78+XqvVateyfv36dS17zYYNG7Zovk099uutXbs273//+zNo0KDcfffdO/wFUM58AgC9U6Xy6qXvWty24PmeSTJs2LCceOKJufHGG/PSS1v3McDjxo3Lww8/3G3ZI488knHjxiVJ9thjjyTJ8uXLu+7/0xcf/alHH3206+vf/e53WbJkSQ466KDN7ru1tTWTJ09OY2Nj7rnnnjQ3N2/Vz9ATznwCAGyDuXPnZtKkSZk4cWKuuOKKvP3tb0+/fv3y2GOP5ZlnnsmECRPecPtLLrkkZ5xxRo444oi85z3vyb333pu77rorP/zhD5O8evb0qKOOylVXXZV99903q1atymc/+9lNPtbs2bMzbNiwtLS05LLLLsvw4cNz2mmnbXLdtWvXZvLkyVm3bl2++c1vprW1Na2trUleDd66urqt/6W8AfEJALAN9t9//zzxxBO58sorM2vWrPz6179OU1NTDj744Fx88cVdLyTanNNOOy033HBDrrnmmlxwwQUZM2ZMbrvtthx//PFd69x6660555xzMnHixIwdOzZf+MIXMnny5I0e66qrrsqFF16YX/7ylznssMNyzz33pLGxcZP7XbhwYX72s58lefXpA39q6dKl2XfffXv2i9hC4hMAYBuNHDkyX/rSl/KlL33pDdfb3EeUf+ITn8gnPvGJzW43bty4/PSnP33Txzr22GPzH//xH1sw8Zt/ZPqO4jmfAAAUIz4BACjGZXcAgD6uVpfQt4YznwAAFCM+AQAoRnwCAL1GX7l0vKvaHsdHfAIANffaRzquW7euxpPwRl47PtvyEZxecAQA1FxdXV2GDh2alStXJkkGDBjwpp9JTjnVajXr1q3LypUrM3To0G369CPxCQD0CiNGjEiSrgCl9xk6dGjXcdpa4hMA6BUqlUpGjhyZPffcMxs2bKj1OLxOQ0PDdvm8d/EJAPQqdXV12yVy6J284AgAgGLEJwAAxbjsDn3Iyxs68oe29lqPURvr2zPoj1++tL491V3197ALqu9XSXODS7CwsxCf0Ic8+ZvfZ8lv19Z6jJqoa1+Xk//49U9/tTrrBg2s6TyUM6CpPsceMFyAwk5CfEIfUl/XL4Oatv6NffuyfnX/83MPaKxPv13097CraWvvyLq29rR3+tQb2FmIT+hDmuvq0r9x1zz706/yPz93c0NdKrvo72FXtKGjs9YjANuRFxwBAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABRTX+sBAHqqX8fL6bfhD7UegwLqNnSkbkN70rY2/smCrdCvPmnoX+spuvGXDPQ5A1f/Rwb8/pe1HoMC1rd3ZuCGjtTVDUkaXKyDHmscmOx3fK8KUPEJ9DnVfvXpaBhU6zEooL3SmfZqe9I0KGmoq/U40Le0tyXrX0o622s9STfiE+hzOuua01nfe/5fPDtOZzrT0dmeNAwQn7A1OtbXeoKNuIYBAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQzFbF59y5czNmzJg0NzdnwoQJeeihh95w/TvuuCOHHXZYBgwYkJEjR+bss8/O6tWrt2pgAAD6rh7H55133pkZM2bksssuyxNPPJHjjjsuJ510UpYtW7bJ9R9++OFMnTo15557bp566qn88z//cx577LF89KMf3ebhAQDoW3ocn9ddd13OPffcfPSjH824ceNy/fXXZ/To0Zk3b94m13/00Uez77775oILLsiYMWNy7LHH5uMf/3gef/zxbR4eAIC+pUfxuX79+ixcuDCTJ0/utnzy5Ml55JFHNrnNMccck1//+teZP39+qtVqfvvb3+Zf/uVfcvLJJ292P21tbWltbe12AwCg7+tRfK5atSodHR1paWnptrylpSUrVqzY5DbHHHNM7rjjjkyZMiWNjY0ZMWJEhg4dmi996Uub3c+cOXMyZMiQrtvo0aN7MiYAAL3UVr3gqFKpdPu+Wq1utOw1Tz/9dC644IL87d/+bRYuXJjvf//7Wbp0aaZNm7bZx581a1bWrFnTdXvhhRe2ZkwAAHqZ+p6sPHz48NTV1W10lnPlypUbnQ19zZw5czJp0qRccsklSZK3v/3tGThwYI477rh87nOfy8iRIzfapqmpKU1NTT0ZDQCAPqBHZz4bGxszYcKELFiwoNvyBQsW5JhjjtnkNuvWrUu/ft13U1dXl+TVM6YAAOw6enzZfebMmfnqV7+aW2+9NYsXL85FF12UZcuWdV1GnzVrVqZOndq1/imnnJK77ror8+bNy7PPPpuf/OQnueCCC/LOd74zo0aN2n4/CQAAvV6PLrsnyZQpU7J69erMnj07y5cvz/jx4zN//vzss88+SZLly5d3e8/Ps846K2vXrs2NN96Yv/mbv8nQoUPz7ne/O1dfffX2+ykAAOgTehyfSTJ9+vRMnz59k/fdfvvtGy07//zzc/7552/NrgAA2In4bHcAAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBitup9PgGgpLb2zlqPAL1WXb9KGuv6zvlE8QlAr9WvkrzS3pGnfrOm1qNAr9XcWJdD3zqkzwSo+ASg12qo65dhA5vSWa31JNA7bejoyCvrO9LRWU3qaj3NlhGfAPRqDX3kbA7USntHe61H6BF/0QAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAiqmv9QDAlmtr78wrGzpqPUZN9Gv/n5+7rT1pr6vhMBTVVJdUKrWeAthexCf0ITP+v0W1HqFm+ueVLG5+9eu//LfmvFzbcSjo4LckX5gkQGFn4bI79HL9G+ryjr2H1noMqJmnf5e07Zon/GGn5Mwn9HKVSiVfP+ed+f5/LM+gxoY0N+6a15v7ta9L/uXVr7/97lfS3tRc24HY4V7pSD50X62nALY38Ql9QKVSSVN9XZoa6tLcsIvGZ/7n526qT+r91wugT3LZHQCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFDMVsXn3LlzM2bMmDQ3N2fChAl56KGH3nD9tra2XHbZZdlnn33S1NSU/fffP7feeutWDQwAQN9V39MN7rzzzsyYMSNz587NpEmT8pWvfCUnnXRSnn766ey9996b3OaMM87Ib3/729xyyy054IADsnLlyrS3t2/z8AAA9C09js/rrrsu5557bj760Y8mSa6//vr84Ac/yLx58zJnzpyN1v/+97+fBx98MM8++2x23333JMm+++67bVMDANAn9eiy+/r167Nw4cJMnjy52/LJkyfnkUce2eQ299xzTyZOnJgvfOELeetb35oDDzwwF198cV5++eXN7qetrS2tra3dbgAA9H09OvO5atWqdHR0pKWlpdvylpaWrFixYpPbPPvss3n44YfT3Nycu+++O6tWrcr06dPz3//935t93uecOXPyd3/3dz0ZDQCAPmCrXnBUqVS6fV+tVjda9prOzs5UKpXccccdeec735kPfOADue6663L77bdv9uznrFmzsmbNmq7bCy+8sDVjAgDQy/TozOfw4cNTV1e30VnOlStXbnQ29DUjR47MW9/61gwZMqRr2bhx41KtVvPrX/86b3vb2zbapqmpKU1NTT0ZDQCAPqBHZz4bGxszYcKELFiwoNvyBQsW5JhjjtnkNpMmTcqLL76YP/zhD13LlixZkn79+mWvvfbaipEBAOirenzZfebMmfnqV7+aW2+9NYsXL85FF12UZcuWZdq0aUlevWQ+derUrvXPPPPMDBs2LGeffXaefvrp/PjHP84ll1ySc845J/37999+PwkAAL1ej99qacqUKVm9enVmz56d5cuXZ/z48Zk/f3722WefJMny5cuzbNmyrvUHDRqUBQsW5Pzzz8/EiRMzbNiwnHHGGfnc5z63/X4KAAD6hB7HZ5JMnz4906dP3+R9t99++0bLDjrooI0u1QMAsOvx2e4AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxdTXegCAnurX0ZbOjldqPQY7WL+OpP+ffF2p1HQc6JX6dXSmX0dH0v5Kkrrud7a/krS3JdVqTWbbHPEJ9DkH/uSiWo9AIYub//jFj2s6BvR+j77BfQe8N2keXGyUN+OyO9AndNb3z9rh76j1GABsI2c+gb6hUsnT77k9uy/7QToaBqWzvvnNt6FPa2tPzrzv1a//aXLS5F8s2Mj69s68vL4jh+89NP3rX3/ZfV3S9oekof+mN64Rf8pA31GppFrX9Meb+NzZdVaTl1/7ui6p1r3h6rBL6qx2prOuPalvThpe/0fSmXRs6HVPmHbZHQCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFbFV8zp07N2PGjElzc3MmTJiQhx56aIu2+8lPfpL6+vocfvjhW7NbAAD6uB7H55133pkZM2bksssuyxNPPJHjjjsuJ510UpYtW/aG261ZsyZTp07Ne97znq0eFgCAvq3H8Xndddfl3HPPzUc/+tGMGzcu119/fUaPHp158+a94XYf//jHc+aZZ+boo4/e6mEBAOjbehSf69evz8KFCzN58uRuyydPnpxHHnlks9vddttt+a//+q9cfvnlW7Sftra2tLa2drsBAND39Sg+V61alY6OjrS0tHRb3tLSkhUrVmxym1/+8pe59NJLc8cdd6S+vn6L9jNnzpwMGTKk6zZ69OiejAkAQC+1VS84qlQq3b6vVqsbLUuSjo6OnHnmmfm7v/u7HHjggVv8+LNmzcqaNWu6bi+88MLWjAkAQC+zZaci/2j48OGpq6vb6CznypUrNzobmiRr167N448/nieeeCKf/OQnkySdnZ2pVqupr6/Pfffdl3e/+90bbdfU1JSmpqaejAYAQB/QozOfjY2NmTBhQhYsWNBt+YIFC3LMMcdstP7gwYPz5JNPZtGiRV23adOmZezYsVm0aFGOPPLIbZseAIA+pUdnPpNk5syZ+fCHP5yJEyfm6KOPzs0335xly5Zl2rRpSV69ZP6b3/wmX//619OvX7+MHz++2/Z77rlnmpubN1oOAMDOr8fxOWXKlKxevTqzZ8/O8uXLM378+MyfPz/77LNPkmT58uVv+p6fAADsmirVarVa6yHeTGtra4YMGZI1a9Zk8ODBtR4HivtDW3v+bfFvM6ipIf0b62o9Ts302/CHvOXX96ejYVA66/vXehx2sFfak9P/9dWvv3NS0tzj0yWw82tr78zL69tzxD5vSf+G1/37sGFd0rY2OfB9SdNuO3yWLe01n+0OAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKqa/1AADwZl7pqPUE0Du1tb/691GtVms9yhYTnwD0eh+6r9YTQG/VL0ljFo5OBjTWepYt47I7AL1SU11y8FtqPQWwvTnzCUCvVKkkX5iUtLnkDpvV1t6Zlze0p3/9gFqPssXEJwC9VqWSNPuXCjarkqTamVQqlVqPssVcdgcAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFbFV8zp07N2PGjElzc3MmTJiQhx56aLPr3nXXXTnxxBOzxx57ZPDgwTn66KPzgx/8YKsHBgCg7+pxfN55552ZMWNGLrvssjzxxBM57rjjctJJJ2XZsmWbXP/HP/5xTjzxxMyfPz8LFy7MCSeckFNOOSVPPPHENg8PAEDfUqlWq9WebHDkkUfmiCOOyLx587qWjRs3LqeddlrmzJmzRY9xyCGHZMqUKfnbv/3bLVq/tbU1Q4YMyZo1azJ48OCejAs7hT+0teffFv82g5oa0r+xrtbj1Ey/DX/IW359fzoaBqWzvn+txwGoubb2zry8vj1H7POW9G943b8PG9YlbWuTA9+XNO22w2fZ0l7r0ZnP9evXZ+HChZk8eXK35ZMnT84jjzyyRY/R2dmZtWvXZvfdd9/sOm1tbWltbe12AwCg7+tRfK5atSodHR1paWnptrylpSUrVqzYosf4h3/4h7z00ks544wzNrvOnDlzMmTIkK7b6NGjezImAAC91Fa94KhSqXT7vlqtbrRsU771rW/liiuuyJ133pk999xzs+vNmjUra9as6bq98MILWzMmAAC9TH1PVh4+fHjq6uo2Osu5cuXKjc6Gvt6dd96Zc889N//8z/+c9773vW+4blNTU5qamnoyGgAAfUCPznw2NjZmwoQJWbBgQbflCxYsyDHHHLPZ7b71rW/lrLPOyj/90z/l5JNP3rpJAQDo83p05jNJZs6cmQ9/+MOZOHFijj766Nx8881ZtmxZpk2bluTVS+a/+c1v8vWvfz3Jq+E5derU3HDDDTnqqKO6zpr2798/Q4YM2Y4/CgAAvV2P43PKlClZvXp1Zs+eneXLl2f8+PGZP39+9tlnnyTJ8uXLu73n51e+8pW0t7fnvPPOy3nnnde1/CMf+Uhuv/32bf8JAADoM3r8Pp+14H0+2dV5n89XeZ9PgO52+vf5BACAbSE+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFBMfa0HALZcW3tHrUeoqboNHVnf3pn2Smc601nrcSikXyVpqHOuBHYW4hP6gPp+lQxoqs+6tvZs6Nh1o6tuQ3sGbuhIe7U9HZ3ttR6HQl5p78iwgU0CFHYS4hP6gOaGuhx7wPC0d1ZrPUptta1NXd2QpGlQ0jCg1tNQQFt7Z576zZrs6v/Th52J+IQ+ormhrtYj9AL1SUO/pKHu1RsAfY5rGAAAFCM+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDH1tR4AAN7Mho6OWo8AvVJf/NsQnwD0WnX9KmlurMsr6zvS3tFe63GgV2purEtdv0qtx9hiWxWfc+fOzTXXXJPly5fnkEMOyfXXX5/jjjtus+s/+OCDmTlzZp566qmMGjUqn/rUpzJt2rStHhrYxbW31XoCCmlMcuieDenodK4ENqeuXyWNna8kna+7o5f+t7LHf8133nlnZsyYkblz52bSpEn5yle+kpNOOilPP/109t57743WX7p0aT7wgQ/kYx/7WL75zW/mJz/5SaZPn5499tgjp59++nb5IYBdRL/6pHFgsv6lpGN9raehkMZaDwC9XWeSzV0YaBz46n87e5FKtVqt9mSDI488MkcccUTmzZvXtWzcuHE57bTTMmfOnI3W//SnP5177rknixcv7lo2bdq0/Pu//3t++tOfbtE+W1tbM2TIkKxZsyaDBw/uybjAzmbDy0mny68AW6RffdLQv8iutrTXepTC69evz8KFC3PppZd2Wz558uQ88sgjm9zmpz/9aSZPntxt2fve977ccsst2bBhQxoaGjbapq2tLW1t/3OquLW1tSdjAjuzQv8RBWDH6NFbLa1atSodHR1paWnptrylpSUrVqzY5DYrVqzY5Prt7e1ZtWrVJreZM2dOhgwZ0nUbPXp0T8YEAKCX2qr3+axUur+iqlqtbrTszdbf1PLXzJo1K2vWrOm6vfDCC1szJgAAvUyPLrsPHz48dXV1G53lXLly5UZnN18zYsSITa5fX1+fYcOGbXKbpqamNDU19WQ0AAD6gB6d+WxsbMyECROyYMGCbssXLFiQY445ZpPbHH300Rutf99992XixImbfL4nAAA7rx5fdp85c2a++tWv5tZbb83ixYtz0UUXZdmyZV3v2zlr1qxMnTq1a/1p06bl+eefz8yZM7N48eLceuutueWWW3LxxRdvv58CAIA+ocdv/DRlypSsXr06s2fPzvLlyzN+/PjMnz8/++yzT5Jk+fLlWbZsWdf6Y8aMyfz583PRRRfly1/+ckaNGpUvfvGL3uMTAGAX1OP3+awF7/MJANC7bWmvbdWr3QEAYGuITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx4hMAgGLEJwAAxYhPAACKEZ8AABQjPgEAKEZ8AgBQjPgEAKAY8QkAQDHiEwCAYsQnAADF1Nd6gC1RrVaTJK2trTWeBACATXmt017rts3pE/G5du3aJMno0aNrPAkAAG9k7dq1GTJkyGbvr1TfLE97gc7Ozrz44ovZbbfdUqlUdvj+WltbM3r06LzwwgsZPHjwDt8fteV471oc712L473rccxrp1qtZu3atRk1alT69dv8Mzv7xJnPfv36Za+99iq+38GDB/sf7i7E8d61ON67Fsd71+OY18YbnfF8jRccAQBQjPgEAKAY8bkJTU1Nufzyy9PU1FTrUSjA8d61ON67Fsd71+OY93594gVHAADsHJz5BACgGPEJAEAx4hMAgGLEJwAAxYjP15k7d27GjBmT5ubmTJgwIQ899FCtR2I7mDNnTv7sz/4su+22W/bcc8+cdtpp+c///M9u61Sr1VxxxRUZNWpU+vfvn+OPPz5PPfVUjSZme5ozZ04qlUpmzJjRtczx3vn85je/yV//9V9n2LBhGTBgQA4//PAsXLiw637HfOfR3t6ez372sxkzZkz69++f/fbbL7Nnz05nZ2fXOo537yU+/8Sdd96ZGTNm5LLLLssTTzyR4447LieddFKWLVtW69HYRg8++GDOO++8PProo1mwYEHa29szefLkvPTSS13rfOELX8h1112XG2+8MY899lhGjBiRE088MWvXrq3h5Gyrxx57LDfffHPe/va3d1vueO9cfve732XSpElpaGjIv/7rv+bpp5/OP/zDP2To0KFd6zjmO4+rr746N910U2688cYsXrw4X/jCF3LNNdfkS1/6Utc6jncvVqXLO9/5zuq0adO6LTvooIOql156aY0mYkdZuXJlNUn1wQcfrFar1WpnZ2d1xIgR1auuuqprnVdeeaU6ZMiQ6k033VSrMdlGa9eurb7tbW+rLliwoPqud72reuGFF1arVcd7Z/TpT3+6euyxx272fsd853LyySdXzznnnG7LPvjBD1b/+q//ulqtOt69nTOff7R+/fosXLgwkydP7rZ88uTJeeSRR2o0FTvKmjVrkiS77757kmTp0qVZsWJFt+Pf1NSUd73rXY5/H3beeefl5JNPznvf+95uyx3vnc8999yTiRMn5n//7/+dPffcM+94xzvyf/7P/+m63zHfuRx77LH50Y9+lCVLliRJ/v3f/z0PP/xwPvCBDyRxvHu7+loP0FusWrUqHR0daWlp6ba8paUlK1asqNFU7AjVajUzZ87Msccem/HjxydJ1zHe1PF//vnni8/Itvv2t7+dX/ziF3nsscc2us/x3vk8++yzmTdvXmbOnJnPfOYz+fnPf54LLrggTU1NmTp1qmO+k/n0pz+dNWvW5KCDDkpdXV06Ojry+c9/Pn/1V3+VxN94byc+X6dSqXT7vlqtbrSMvu2Tn/xk/u///b95+OGHN7rP8d85vPDCC7nwwgtz3333pbm5ebPrOd47j87OzkycODFXXnllkuQd73hHnnrqqcybNy9Tp07tWs8x3znceeed+eY3v5l/+qd/yiGHHJJFixZlxowZGTVqVD7ykY90red4904uu//R8OHDU1dXt9FZzpUrV270/5zou84///zcc889uf/++7PXXnt1LR8xYkSSOP47iYULF2blypWZMGFC6uvrU19fnwcffDBf/OIXU19f33VMHe+dx8iRI3PwwQd3WzZu3LiuF4z6G9+5XHLJJbn00kvzl3/5lzn00EPz4Q9/OBdddFHmzJmTxPHu7cTnHzU2NmbChAlZsGBBt+ULFizIMcccU6Op2F6q1Wo++clP5q677sq//du/ZcyYMd3uHzNmTEaMGNHt+K9fvz4PPvig498Hvec978mTTz6ZRYsWdd0mTpyYD33oQ1m0aFH2228/x3snM2nSpI3ePm3JkiXZZ599kvgb39msW7cu/fp1T5i6urqut1pyvHu5Gr7Yqdf59re/XW1oaKjecsst1aeffro6Y8aM6sCBA6vPPfdcrUdjG33iE5+oDhkypPrAAw9Uly9f3nVbt25d1zpXXXVVdciQIdW77rqr+uSTT1b/6q/+qjpy5Mhqa2trDSdne/nTV7tXq473zubnP/95tb6+vvr5z3+++stf/rJ6xx13VAcMGFD95je/2bWOY77z+MhHPlJ961vfWv3ud79bXbp0afWuu+6qDh8+vPqpT32qax3Hu/cSn6/z5S9/ubrPPvtUGxsbq0cccUTXW/HQtyXZ5O22227rWqezs7N6+eWXV0eMGFFtamqq/q//9b+qTz75ZO2GZrt6fXw63jufe++9tzp+/PhqU1NT9aCDDqrefPPN3e53zHcera2t1QsvvLC69957V5ubm6v77bdf9bLLLqu2tbV1reN4916VarVareWZVwAAdh2e8wkAQDHiEwCAYsQnAADFiE8AAIoRnwAAFCM+AQAoRnwCAFCM+AR4Aw888EAqlUp+//vf79D93H777Rk6dGjX91dccUUOP/zwHbpPgFoQnwB/4vjjj8+MGTO6vj/mmGOyfPnyDBkypOgcF198cX70ox8V3SdACfW1HgCgN2tsbMyIESOK73fQoEEZNGhQ8f0C7GjOfAL80VlnnZUHH3wwN9xwQyqVSiqVSm6//fZul91fuzz+3e9+N2PHjs2AAQPyF3/xF3nppZfyta99Lfvuu2/e8pa35Pzzz09HR0fXY69fvz6f+tSn8ta3vjUDBw7MkUcemQceeGCzs7z+svtZZ52V0047Lddee21GjhyZYcOG5bzzzsuGDRu2eh8AteDMJ8Af3XDDDVmyZEnGjx+f2bNnJ0meeuqpjdZbt25dvvjFL+bb3/521q5dmw9+8IP54Ac/mKFDh2b+/Pl59tlnc/rpp+fYY4/NlClTkiRnn312nnvuuXz729/OqFGjcvfdd+f9739/nnzyybztbW/bovnuv//+jBw5Mvfff39+9atfZcqUKTn88MPzsY99bLvtA2BHE58AfzRkyJA0NjZmwIABXZfan3nmmY3W27BhQ+bNm5f9998/SfIXf/EX+cY3vpHf/va3GTRoUA4++OCccMIJuf/++zNlypT813/9V771rW/l17/+dUaNGpXk1ed0fv/7389tt92WK6+8covme8tb3pIbb7wxdXV1Oeigg3LyySfnRz/6UT72sY9tt30A7GjiE6CHBgwY0BWeSdLS0pJ9992323M0W1pasnLlyiTJL37xi1Sr1Rx44IHdHqetrS3Dhg3b4v0ecsghqaur6/p+5MiRefLJJ7frPgB2NPEJ0EMNDQ3dvq9UKptc1tnZmSTp7OxMXV1dFi5c2C0ek/ToRUUl9gGwo4lPgD/R2NjY7YVC28M73vGOdHR0ZOXKlTnuuOO262OX3AfA9uDV7gB/Yt99983PfvazPPfcc1m1alXXmcVtceCBB+ZDH/pQpk6dmrvuuitLly7NY489lquvvjrz58/fDlOX2QfA9iA+Af7ExRdfnLq6uhx88MHZY489smzZsu3yuLfddlumTp2av/mbv8nYsWNz6qmn5mc/+1lGjx69XR6/1D4AtlWlWq1Waz0EAAC7Bmc+AQAoRnwCAFCM+AQAoBjxCQBAMeITAIBixCcAAMWITwAAihGfAAAUIz4BAChGfAIAUIz4BACgGPEJAEAx/z9lZX1Zksj5OgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots(figsize=(8,8))\n",
    "by_X = {}\n",
    "for group, df12 in df.groupby('group'):\n",
    "    by_X[group] = df12\n",
    "    km_group = km.fit(df12['Y'], df12['delta'])\n",
    "    km_group.plot(label='Group %s' % group, ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "441b5707",
   "metadata": {},
   "source": [
    "We see the confidence intervals of both groups are quite the same. We conjecture that there is no significant difference between the two survival curves."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "895649a9",
   "metadata": {},
   "source": [
    "### (b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "id": "3fd3e998",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>se(coef)</th>\n",
       "      <th>p</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>covariate</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>group</th>\n",
       "      <td>-0.340143</td>\n",
       "      <td>1.235876</td>\n",
       "      <td>0.783144</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               coef  se(coef)         p\n",
       "covariate                              \n",
       "group     -0.340143  1.235876  0.783144"
      ]
     },
     "execution_count": 143,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coxph = CoxPHFitter\n",
    "group_df=df[['Y','delta','group']]\n",
    "model_df=MS(['Y','delta','group'],intercept=False).fit_transform(group_df)\n",
    "cox_fit = coxph().fit(model_df,'Y','delta')\n",
    "cox_fit.summary[['coef', 'se(coef)', 'p']]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7d29536",
   "metadata": {},
   "source": [
    "The baseline group here is 1. The result shows that the risk associated with group 1\n",
    "is more than 1.40 times (i.e. $e^{0.34}=1.404$) the risk associated\n",
    "with group 2."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b82845b",
   "metadata": {},
   "source": [
    "The resulting 𝑝-value is 0.78, indicating no evidence of a difference in survival between the two groups. Therefore, the true coefficient is significantly null."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "871741fa",
   "metadata": {},
   "source": [
    "### (c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "id": "50912ecb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{1:       Y  delta    X  group\n",
       " 0  26.5      1  0.1      1\n",
       " 2  57.3      1 -0.3      1\n",
       " 4  20.2      0  1.8      1\n",
       " 5  89.8      0  0.4      1,\n",
       " 2:       Y  delta     X  group\n",
       " 1  37.2      1  11.0      2\n",
       " 3  90.8      0   2.8      2}"
      ]
     },
     "execution_count": 144,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "by_X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "id": "eff6d1bb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>t_0</th>\n",
       "      <td>-1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>null_distribution</th>\n",
       "      <td>chi squared</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>degrees_of_freedom</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>test_name</th>\n",
       "      <td>logrank_test</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div><table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>test_statistic</th>\n",
       "      <th>p</th>\n",
       "      <th>-log2(p)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.08</td>\n",
       "      <td>0.78</td>\n",
       "      <td>0.35</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>"
      ],
      "text/latex": [
       "\\begin{tabular}{lrrr}\n",
       " & test_statistic & p & -log2(p) \\\\\n",
       "0 & 0.08 & 0.78 & 0.35 \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "<lifelines.StatisticalResult: logrank_test>\n",
       "               t_0 = -1\n",
       " null_distribution = chi squared\n",
       "degrees_of_freedom = 1\n",
       "         test_name = logrank_test\n",
       "\n",
       "---\n",
       " test_statistic    p  -log2(p)\n",
       "           0.08 0.78      0.35"
      ]
     },
     "execution_count": 146,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logrank_test(by_X[1]['Y'],\n",
    "             by_X[2]['Y'],\n",
    "             by_X[1]['delta'],\n",
    "             by_X[2]['delta'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2c8ed08",
   "metadata": {},
   "source": [
    "The resulting 𝑝-value is 0.78 for the log-rank test, which is the same value found for the Cox model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d0d673f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}