{
"cells": [
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"tags": [
"remove_input"
]
},
"outputs": [],
"source": [
"path_data = '../../data/'\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"import matplotlib\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('fivethirtyeight')\n",
"\n",
"import warnings\n",
"warnings.filterwarnings('ignore')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# A/B Testing\n",
"In modern data analytics, deciding whether two numerical samples come from the same underlying distribution is called *A/B testing*. The name refers to the labels of the two samples, A and B.\n",
"\n",
"We will develop the method in the context of an example. The data come from a sample of newborns in a large hospital system. We will treat it as if it were a simple random sample though the sampling was done in multiple stages. [Stat Labs](https://www.stat.berkeley.edu/~statlabs/) by Deborah Nolan and Terry Speed has details about a larger dataset from which this set is drawn. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Smokers and Nonsmokers ###\n",
"The table `births` contains the following variables for 1,174 mother-baby pairs: the baby's birth weight in ounces, the number of gestational days, the mother's age in completed years, the mother's height in inches, pregnancy weight in pounds, and whether or not the mother smoked during pregnancy."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Birth Weight
\n",
"
Gestational Days
\n",
"
Maternal Age
\n",
"
Maternal Height
\n",
"
Maternal Pregnancy Weight
\n",
"
Maternal Smoker
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
120
\n",
"
284
\n",
"
27
\n",
"
62
\n",
"
100
\n",
"
False
\n",
"
\n",
"
\n",
"
1
\n",
"
113
\n",
"
282
\n",
"
33
\n",
"
64
\n",
"
135
\n",
"
False
\n",
"
\n",
"
\n",
"
2
\n",
"
128
\n",
"
279
\n",
"
28
\n",
"
64
\n",
"
115
\n",
"
True
\n",
"
\n",
"
\n",
"
3
\n",
"
108
\n",
"
282
\n",
"
23
\n",
"
67
\n",
"
125
\n",
"
True
\n",
"
\n",
"
\n",
"
4
\n",
"
136
\n",
"
286
\n",
"
25
\n",
"
62
\n",
"
93
\n",
"
False
\n",
"
\n",
"
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
\n",
"
\n",
"
1169
\n",
"
113
\n",
"
275
\n",
"
27
\n",
"
60
\n",
"
100
\n",
"
False
\n",
"
\n",
"
\n",
"
1170
\n",
"
128
\n",
"
265
\n",
"
24
\n",
"
67
\n",
"
120
\n",
"
False
\n",
"
\n",
"
\n",
"
1171
\n",
"
130
\n",
"
291
\n",
"
30
\n",
"
65
\n",
"
150
\n",
"
True
\n",
"
\n",
"
\n",
"
1172
\n",
"
125
\n",
"
281
\n",
"
21
\n",
"
65
\n",
"
110
\n",
"
False
\n",
"
\n",
"
\n",
"
1173
\n",
"
117
\n",
"
297
\n",
"
38
\n",
"
65
\n",
"
129
\n",
"
False
\n",
"
\n",
" \n",
"
\n",
"
1174 rows × 6 columns
\n",
"
"
],
"text/plain": [
" Birth Weight Gestational Days Maternal Age Maternal Height \\\n",
"0 120 284 27 62 \n",
"1 113 282 33 64 \n",
"2 128 279 28 64 \n",
"3 108 282 23 67 \n",
"4 136 286 25 62 \n",
"... ... ... ... ... \n",
"1169 113 275 27 60 \n",
"1170 128 265 24 67 \n",
"1171 130 291 30 65 \n",
"1172 125 281 21 65 \n",
"1173 117 297 38 65 \n",
"\n",
" Maternal Pregnancy Weight Maternal Smoker \n",
"0 100 False \n",
"1 135 False \n",
"2 115 True \n",
"3 125 True \n",
"4 93 False \n",
"... ... ... \n",
"1169 100 False \n",
"1170 120 False \n",
"1171 150 True \n",
"1172 110 False \n",
"1173 129 False \n",
"\n",
"[1174 rows x 6 columns]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"births = pd.read_csv(path_data + 'baby.csv')\n",
"births"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One of the aims of the study was to see whether maternal smoking was associated with birth weight. Let's see what we can say about the two variables.\n",
"\n",
"We'll start by selecting just `Birth Weight` and `Maternal Smoker`. There are 715 non-smokers among the women in the sample, and 459 smokers."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"smoking_and_birthweight = births[['Maternal Smoker', 'Birth Weight']]"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Maternal Smoker
\n",
"
count
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
False
\n",
"
715
\n",
"
\n",
"
\n",
"
1
\n",
"
True
\n",
"
459
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Maternal Smoker count\n",
"0 False 715\n",
"1 True 459"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#smoking_birthweight = smoking_and_birthweight.groupby('Maternal Smoker').count()\n",
"\n",
"smoking_birthweight1 = smoking_and_birthweight.groupby([\"Maternal Smoker\"]).agg(\n",
" count=pd.NamedAgg(column=\"Maternal Smoker\", aggfunc=\"count\")\n",
")\n",
"\n",
"smoking_birthweight1.reset_index()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"smoker = births[births['Maternal Smoker'] == False]\n",
"\n",
"non_smoker = births[births['Maternal Smoker'] == True]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's look at the distribution of the birth weights of the babies of the non-smoking mothers compared to those of the smoking mothers. To generate two overlaid histograms, we will use `hist` with the optional `group` argument which is a column label or index. The rows of the table are first grouped by this column and then a histogram is drawn for each one."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA60AAAFDCAYAAADRdOgEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABA2UlEQVR4nO3deVyU9f7//+cMm4omZgbIkgKmoagdPWlalnpKo3A3XOt4ci/OyTTX1NzChTxpollU3xY1LTWXXEqPmVLoUTPtmGapiSa4hQvJNjO/P/o5n0a2QWeYAR73243bzbmu93XN63qBA8+5rus9hoyMDIsAAAAAAHBDRlcXAAAAAABAYQitAAAAAAC3RWgFAAAAALgtQisAAAAAwG0RWgEAAAAAbovQCgAAAABwW4RWAAAAAIDbIrQCAAAAANwWobUMOXr0qKtLqJDou2vQd9eg765B312H3rsGfQdQEi4LrXPnzlXbtm0VEhKi8PBwxcbG6tChQ0Vu88svv8jPzy/f15YtW0qpagAAAABAafJ01RPv3LlTzzzzjP7yl7/IYrHolVdeUZcuXbRr1y7VqFGjyG1XrlypRo0aWR8XNx4AAAAAUDa5LLSuWrXK5vHixYsVGhqqlJQUPfbYY0Vue/vtt8vf39+Z5QEAAAAA3IDb3NN69epVmc1m+fn5FTu2f//+ioiIUIcOHbRmzRrnFwcAAAAAcAlDRkaGxdVFSNLf//53/fzzz/ryyy/l4eFR4JgLFy5o6dKlatmypTw9PbVhwwa9+uqrWrRokWJjYwvdNzf7AwAAlB316tVzdQmSpLy8PGVmZrq6DKDc8/X1ladn4RcBu0VoHT9+vFatWqVNmzapTp06Jdp25MiR+uabb/T11187pzg3cvToUbd5Ea9I6Ltr0HfXoO+uQd9dh967Rlnoe15enq5cuSI/Pz8ZDAZXlwOUWxaLRRkZGapWrVqhwdXllwePGzdOK1eu1Nq1a0scWCWpWbNmOnbsmOMLAwAAQIWVmZlJYAVKgcFgkJ+fX5FXNbhsIiZJGjNmjFatWqX169fr7rvvvql9HDx4kEmZAAAA4HAEVqB0FPd/zWWhddSoUVq+fLk+/PBD+fn5KT09XdIf1zNXrVpVkjRlyhTt3btXa9eulSQtXbpUXl5eaty4sYxGozZt2qSkpCS9/PLLrjoMAAAAAIATuSy0JiUlSZI6d+5ss3zMmDEaN26cJCktLU3Hjx+3WZ+QkKDU1FR5eHgoPDxcCxYsKHISJgAAAABA2eWy0JqRkVHsmEWLFtk87tOnj/r06eOkigAAAICinTnjrdOnS29amKAgswIDc0rt+cqKxx9/XJGRkZozZ46rS8lnyZIlGj16tE6fPu3qUsoNl97TCgC4ec74w+ny5Tq6dKmSQ/d5K/hjDYC7OX3aqFGjCv54RmdISJACA+0fP2zYMC1btkz9+/fX66+/brNu0qRJmj9/vjp06KDly5fbvU8/Pz+99957+a6QdHfvv/++3nrrLR07dkweHh4KDg5WdHS0XnrpJVeX5haioqKUmppa6PrWrVvrs88+K8WKCkdoBYAyyhl/OGVnV5aPT+n9MVackv6xBgCQgoODtXr1as2cOVO+vr6S/vgIn+XLlys4ONhldeXk5Mjb27tUnuuDDz7QmDFj9Morr+ihhx5STk6ODh8+rN27d5fK8ztabm6uvLy8HLrPbdu2yWQySZK+//57de/eXf/5z38UFBQkSfm+V6X5/buRyz/yBgAAAIDjNGzYUGFhYVq9erV12ebNm+Xj46MHHnjAZuy+ffvUtWtXhYWFKSQkRB07drQJdlFRUZKkp59+Wn5+ftbHkrRx40Y99NBD8vf3V+PGjTVt2jTl5OTYbBsfH69nn31WoaGhGjRokJYsWaKgoCBt375d999/v2rXrq0nnnhCJ06csG53/Phx9e7dW3fffbdq166tNm3aaNOmTSXqwcaNGxUTE6MBAwYoLCxMDRo0UJcuXfTKK69Yx8THx+v+++/X0qVLFRUVpaCgIA0fPlw5OTlKSkpSw4YNVbduXY0fP15ms9m6XUZGhoYOHaq77rpLAQEB6ty5s3744YdCa8nIyFCHDh3UrVs3ZWZmymKxaN68eWratKkCAgLUqlUrmzPfv/zyi/z8/PTJJ58oJiZGAQEBevfdd0t0/Pa444475O/vL39/f91+++2SpJo1a1qX1a1bV2+99Zb69eun2rVra+rUqdqxY4f8/Px04cKFfPV+++231mWHDx/Wk08+qeDgYEVEROiZZ56xTrx7MwitAAAAQDnTv39/LVmyxPr4ww8/VN++ffN9tMiVK1cUGxurjRs3auvWrYqKilLPnj2toWTbtm2SpPnz5+vIkSPWx1u3btXgwYM1aNAgpaSkaMGCBVqzZo2mTp1qs/+FCxfq7rvv1pdffqlJkyZJkrKzszV37lwtWLBAn3/+uS5duqQXXnjBus3Vq1f1yCOPaPXq1dq5c6c6deqk/v3768cff7T7+P39/bV3716bMFyQkydPasOGDVq+fLnef/99rVmzRn369NG+ffu0atUqzZ8/X2+++abWrVtn3WbYsGHau3evli5dqq1bt6py5crq0aOHrl27lm//aWlpio6OVmBgoD766CP5+vpq+vTp+uCDD5SQkKCUlBSNGDFCI0aM0ObNm222nTJligYOHKiUlBQ9/vjjBdY/YsQIBQUFFflV1CXAxZk1a5YeffRRff311xo4cKBd21w/5nvuuUdbt27Vp59+qqtXr6p379424b8kuDwYAAAAKGd69uypiRMn6ueff1bVqlW1detWzZ492+ZMoyQ99NBDNo9nz56ttWvXasuWLYqNjdUdd9whSapevbr8/f2t4xISEhQXF6d+/fpJkurWrauXX35ZQ4YM0bRp06zhuFWrVvrXv/5l3S4lJUV5eXlKSEhQvXr1JElxcXF69tlnZTabZTQaFRUVZXNGd9SoUdq0aZPWrFmjF1980a7jHzNmjL7//ns1bdpUYWFhat68udq2basePXrYXGZrMpmUmJio6tWrKzIyUu3bt1dycrJ++OEHeXt7q379+mrRooV27typzp076+eff9bGjRv12WefqXXr1pKkxYsXKyoqSh9//LGeeuop676PHTumrl27qn379kpISJDRaFRmZqYSExO1atUqtWrVSpJUp04d7d27V0lJSerQoYN1+8GDBxd7H/H48eMVFxdX5JjAW7jPpmvXrjbHZE8Afvvtt9WoUSNNmTLFumzx4sWqU6eOvv32WzVr1qzEdRBaAQAAgHLGz89PTzzxhD788ENVr15dDzzwgEJCQvKNO3funGbMmKEdO3bo3LlzMplMunbtmk6dOlXk/r/77jvt27dP8+bNsy4zm826du2a0tPTFRAQIEm69957823r4+NjDaySFBAQoNzcXF26dEk1atRQZmamZs2apc2bNystLU15eXnKyspSw4YN7T7+gIAAffHFFzp06JCSk5O1e/dujRgxQgsXLtTmzZtVpUoVSX/c/1u9enXrdnfeeaciIiJs7t288847de7cOUnSkSNHZDQadd9991nXXw+8hw8fti7LyclRx44d1alTJyUkJFiXHzlyRFlZWerRo4fNWe/c3FyFhobaHENBvbtRrVq1VKtWLXvbUmL21HCj7777Tl9//bX13tg/O378OKEVAAAAwB/69eunYcOGydfXV+PHjy9wzLBhw3T27Fm98sorCg0NlY+Pjzp16mRzb2pBzGazxowZoy5duuRbd/3srCTrRFB/5ulpG0Guh7frl45OnDhRW7Zs0bRp0xQeHq4qVapo6NChxdZUkMjISEVGRmrQoEH65ptv9Nhjj2n16tXq27evJOWb3MhgMBRY3/UJiywWS6HP9ecQ6uXlpbZt2+rzzz/XyZMnrYH0+jEuW7Ys35sINz5vQb270YgRI7RixYoix6SkpBT4hoU9bqzBaPzj7tI/9yEvL89mjNls1qOPPqrp06fn29/NBmxCKwAAAFAOPfTQQ/Ly8tKFCxcKvScyJSVFM2fOtF6Wevbs2XwT5nh5eVlD23VNmjTRjz/+qLCwMIfXnZKSol69elkvjc3KytLx48cVHh5+S/tt0KCBJCkzM/OW9mE2m7V7927r5cGXL1/WoUOH1KdPH+s4g8GgRYsWaejQoYqJidH69esVEhKi+vXry8fHR6mpqfkuzb4Zzr48+EbX35BIS0uz/vvgwYM2Y5o0aaLVq1crJCTEYTMeE1oBAACAcshgMCg5OVkWi0U+Pj4FjgkPD9eKFSvUvHlz/f7775o0aVK+jzUJDQ3V9u3b1bp1a/n4+MjPz0+jR49WbGysQkJC1LVrV3l6euqHH37Q3r17803GVFLh4eFav369oqOj5eXlpVmzZik7O7tE+3jhhRcUEBCgNm3aqHbt2kpPT1dCQoKqVKmidu3a3VJt0dHRGjFihF577TVVr15d06ZNU7Vq1dSzZ0+bsUajUW+88YaGDh2qJ554whpc4+LiNHHiRFksFrVu3VpXr17Vnj17ZDQa9fe//71E9Tj78uAbhYWFKTg4WDNnztTLL7+skydPas6cOTZjBg4cqPfee08DBgzQ888/rzvuuEMnTpzQ6tWrNX36dFWrVq3Ez0toBQAAAOwUFGTWn25RLJXnuxXFBYQFCxbo+eef18MPP6yAgACNHTvW5uNMJGn69OmaMGGCGjZsqMDAQB08eFDt27fXihUrNGfOHC1YsECenp4KDw+3Odt4s2bMmKG4uDhFR0fLz89Pw4YNK3Foffjhh7VkyRK9++67unDhgmrUqKGmTZtq9erVioiIuKX6Fi5cqLFjx6p3797Kzs5WixYt9Mknn6hy5cr5xhqNRi1atEjDhg1TTEyM1q1bpwkTJqhWrVpasGCBRo4cqWrVqikqKspmwip35eXlpbffflsjR47UAw88oKioKE2aNEmxsbHWMYGBgdq8ebOmTJmi7t27Kzs7W8HBwWrbtm2hb54Ux5CRkVH4hdlwK0ePHrW5aR2lg767Bn0v3p49lTRqlIdD95mdnX3Tv1CcISHBpObNs1xdhtPx8+469N41ykLfL126ZDNBDwDnKur/HJ/TCgAAAABwW4RWAAAAAIDbIrQCAAAAANwWoRUAAAAA4LYIrQAAAAAAt0VoBQAAAAC4LUIrAAAAAMBtEVoBAAAAAG6L0AoAAAAAcFueri4AAAAAKCu8jWdkNJ0uteczewQpxxxYas9XVjz++OOKjIzUnDlzXF1KPkuWLNHo0aN1+nTp/ZyUd4RWAAAAwE5G02l5XBhVek9YM0Ey2B9ahw0bpmXLlql///56/fXXbdZNmjRJ8+fPV4cOHbR8+XK79+nn56f33ntPnTt3tnsbd/D+++/rrbfe0rFjx+Th4aHg4GBFR0frpZdecnVpKCEuDwYAAADKkeDgYK1evVqZmZnWZXl5eVq+fLmCg4NdVldOTk6pPdcHH3ygMWPG6B//+Id27Nihzz//XKNHj9bvv/9eajU4Um5urqtLcClCKwAAAFCONGzYUGFhYVq9erV12ebNm+Xj46MHHnjAZuy+ffvUtWtXhYWFKSQkRB07dtTu3but66OioiRJTz/9tPz8/KyPJWnjxo166KGH5O/vr8aNG2vatGk2wTQqKkrx8fF69tlnFRoaqkGDBmnJkiUKCgrS9u3bdf/996t27dp64okndOLECet2x48fV+/evXX33Xerdu3aatOmjTZt2lSiHmzcuFExMTEaMGCAwsLC1KBBA3Xp0kWvvPKKdUx8fLzuv/9+LV26VFFRUQoKCtLw4cOVk5OjpKQkNWzYUHXr1tX48eNlNput22VkZGjo0KG66667FBAQoM6dO+uHH34otJaMjAx16NBB3bp1U2ZmpiwWi+bNm6emTZsqICBArVq1sjnz/csvv8jPz0+ffPKJYmJiFBAQoHfffbdEx1/eEFoBAACAcqZ///5asmSJ9fGHH36ovn37ymAw2Iy7cuWKYmNjtXHjRm3dulVRUVHq2bOnLly4IEnatm2bJGn+/Pk6cuSI9fHWrVs1ePBgDRo0SCkpKVqwYIHWrFmjqVOn2ux/4cKFuvvuu/Xll19q0qRJkqTs7GzNnTtXCxYs0Oeff65Lly7phRdesG5z9epVPfLII1q9erV27typTp06qX///vrxxx/tPn5/f3/t3bvXJgwX5OTJk9qwYYOWL1+u999/X2vWrFGfPn20b98+rVq1SvPnz9ebb76pdevWWbcZNmyY9u7dq6VLl2rr1q2qXLmyevTooWvXruXbf1pamqKjoxUYGKiPPvpIvr6+mj59uj744AMlJCQoJSVFI0aM0IgRI7R582abbadMmaKBAwcqJSVFjz/+eIH1jxgxQkFBQUV+paam2t03d8U9rQAAAEA507NnT02cOFE///yzqlatqq1bt2r27Nk2Zxol6aGHHrJ5PHv2bK1du1ZbtmxRbGys7rjjDklS9erV5e/vbx2XkJCguLg49evXT5JUt25dvfzyyxoyZIimTZtmDcetWrXSv/71L+t2KSkpysvLU0JCgurVqydJiouL07PPPiuz2Syj0aioqCibM7qjRo3Spk2btGbNGr344ot2Hf+YMWP0/fffq2nTpgoLC1Pz5s3Vtm1b9ejRQ15eXtZxJpNJiYmJql69uiIjI9W+fXslJyfrhx9+kLe3t+rXr68WLVpo586d6ty5s37++Wdt3LhRn332mVq3bi1JWrx4saKiovTxxx/rqaeesu772LFj6tq1q9q3b6+EhAQZjUZlZmYqMTFRq1atUqtWrSRJderU0d69e5WUlKQOHTpYtx88eHCx9xGPHz9ecXFxRY4JDCz7E3kRWgEAAIByxs/PT0888YQ+/PBDVa9eXQ888IBCQkLyjTt37pxmzJihHTt26Ny5czKZTLp27ZpOnTpV5P6/++477du3T/PmzbMuM5vNunbtmtLT0xUQECBJuvfee/Nt6+PjYw2skhQQEKDc3FxdunRJNWrUUGZmpmbNmqXNmzcrLS1NeXl5ysrKUsOGDe0+/oCAAH3xxRc6dOiQkpOTtXv3bo0YMUILFy7U5s2bVaVKFUl/3P9bvXp163Z33nmnIiIi5O3tbbPs3LlzkqQjR47IaDTqvvvus66/HngPHz5sXZaTk6OOHTuqU6dOSkhIsC4/cuSIsrKy1KNHD5uz3rm5uQoNDbU5hoJ6d6NatWqpVq1a9ralzCK0AgAAAOVQv379NGzYMPn6+mr8+PEFjhk2bJjOnj2rV155RaGhofLx8VGnTp2KnTTJbDZrzJgx6tKlS75118/OSpKvr2++9Z6ethHkeni7ft/oxIkTtWXLFk2bNk3h4eGqUqWKhg4delMTOUVGRioyMlKDBg3SN998o8cee0yrV69W3759JcnmrOv1Wgqqz2QySZIsFkuhz/XnEOrl5aW2bdvq888/18mTJ62B9PoxLlu2LN+bCDc+b0G9u9GIESO0YsWKIsekpKQU+IZFWUJoBQAAAMqhhx56SF5eXrpw4UKh90SmpKRo5syZ1stSz549q/T0dJsxXl5e1tB2XZMmTfTjjz8qLCzM4XWnpKSoV69e1ktjs7KydPz4cYWHh9/Sfhs0aCBJNrMq38w+zGazdu/ebb08+PLlyzp06JD69OljHWcwGLRo0SINHTpUMTExWr9+vUJCQlS/fn35+PgoNTU136XZN4PLgwEAAACUWQaDQcnJybJYLPLx8SlwTHh4uFasWKHmzZvr999/16RJk2wujZWk0NBQbd++Xa1bt5aPj4/8/Pw0evRoxcbGKiQkRF27dpWnp6d++OEH7d27N99kTCUVHh6u9evXKzo6Wl5eXpo1a5ays7NLtI8XXnhBAQEBatOmjWrXrq309HQlJCSoSpUqateu3S3VFh0drREjRui1115T9erVNW3aNFWrVk09e/a0GWs0GvXGG29o6NCheuKJJ6zBNS4uThMnTpTFYlHr1q119epV7dmzR0ajUX//+99LVA+XBwMAAACwYfYIkmomFD/Qkc9nLn5cYapVq1bk+gULFuj555/Xww8/rICAAI0dO9Y6c/B106dP14QJE9SwYUMFBgbq4MGDat++vVasWKE5c+ZowYIF8vT0VHh4uM3Zxps1Y8YMxcXFKTo6Wn5+fho2bFiJQ+vDDz+sJUuW6N1339WFCxdUo0YNNW3aVKtXr1ZERMQt1bdw4UKNHTtWvXv3VnZ2tlq0aKFPPvlElStXzjfWaDRq0aJFGjZsmGJiYrRu3TpNmDBBtWrV0oIFCzRy5EhVq1ZNUVFRNhNWwZYhIyOj8Auz4VaOHj1qc9M6Sgd9dw36Xrw9eypp1CgPh+4zOzu70HfjXSEhwaTmzbNcXYbT8fPuOvTeNcpC3y9dumQzQQ8A5yrq/xyf0woAAAAAcFuEVgAAAACA2yK0AgAAAADcFqEVAAAAAOC2CK0AAAAAALdFaAUAAAAKYLHwIRtAaSju/xqhFQAAALiBr6+vMjIyCK6Ak1ksFmVkZMjX17fQMZ6lWA8AVDjexjMymk47Zd+RdYyaNt7g0H2aTWYZPWzfzzyVFqLF79zl0OcBAHfn6empatWq6fLly64uBSj3qlWrJk/PwqMpoRUAnMhoOi2PC6Ocsm/fTKPCCn9T8qaYzRYZjTcE4YBXJRFaAVQ8np6eql69uqvLACo8Lg8GAAAAALgtQisAAAAAwG0RWgEAAAAAbot7WgEAbstg8NCePZVcXYZTBQWZXV0CAABujdAKAHBb589L8fEeri7DqRISJOZ5AQCgcGXu8uC5c+eqbdu2CgkJUXh4uGJjY3Xo0CFXlwUAAAAAcIIyF1p37typZ555Rps3b9batWvl6empLl266LfffnN1aQAAAAAABytzlwevWrXK5vHixYsVGhqqlJQUPfbYYy6qCoCjBd6RrUqWPa4u45Z5GK65ugQAAIAyrcyF1htdvXpVZrNZfn5+ri4FgAP5eJyVx4Wpri7j1tUc5+oKAAAAyrQyH1rHjh2rqKgo3XfffYWOOXr0aClW5Fzl6VjKEvpe+ur4S9lZ2a4u45Z55ObK5KTjMJm8ZTY7/i4Ps9li+9hkVna2a74Xubkeys42ueS5S8vly9dUvTqvM65E712jqL7Xq1evFCsB4O7KdGgdP368UlJStGnTJnl4FD67ZHl54Tt69Gi5OZayhL67Ru7lVPlU8nF1GbfOy0ueTjoOk8koo4Mzq9lskdFosFlm9DDKx8c13wsvL8nHp0z/qirWbbf9cXy8zrgGr/GuQd8BlESZ/Utg3LhxWrVqldatW6c6deq4uhwAAAAAgBOUydA6ZswYrVq1SuvXr9fdd9/t6nIAAAAAAE5S5kLrqFGjtHz5cn344Yfy8/NTenq6JMnX11dVq1Z1cXUAAAAAAEcqc5/TmpSUpCtXrqhz586qX7++9ev11193dWkAAAAAAAcrc2daMzIyXF0CAAAAAKCUlLkzrQAAAACAiqPMnWkFAJSuehEWTRu/0yXP3aiRNG38re/nVFqIFr9z163vCAAAlDpCKwCgSL6VzivMd6ZLnvu2a1KYrwN2FPCqJEIrAABlEZcHAwAAAADcFqEVAAAAAOC2CK0AAAAAALdFaAUAAAAAuC1CKwAAAADAbdkdWpOTk3X+/PlC11+4cEHJyckOKQoAAAAAAKkEoTUmJkbbtm0rdP327dsVExPjkKIAAAAAAJBKEFotFkuR63NycmQ0crUxAAAAAMBxPItaefnyZV26dMn6+OLFi0pNTc03LiMjQytXrlRgYKDjKwQAAAAAVFhFhtaFCxdq9uzZkiSDwaBx48Zp3LhxBY61WCyaOHGi4ysEAAAAAFRYRYbWhx9+WJUqVZLFYtHUqVPVrVs3RUVF2YwxGAyqUqWK7r33XjVv3typxQIAAAAAKpYiQ2vLli3VsmVLSVJ2drZiYmLUsGHDUikMAAAAAIAiQ+ufjR071pl1AAAAAACQT6GhddmyZZKkXr16yWAwWB8Xp3fv3o6pDAAAAABQ4RUaWocPHy6DwaDu3bvL29tbw4cPL3ZnBoOB0AoAAAAAcJhCQ+t3330nSfL29rZ5DAAAAABAaSk0tIaGhhb5GAAAAAAAZzO6ugAAAAAAAApj9+zBkvTll1/qvffe04kTJ/Tbb7/JYrHYrDcYDNq/f78j6wMAAAAAVGB2h9ZFixZpwoQJuuOOO9S8eXPdc889zqwLAAAAAAD7Q2tiYqJat26tlStXWidnAgAAAADAmey+p/XChQvq1q0bgRUAAAAAUGrsDq1NmzbVyZMnnVkLAAAAAAA27A6tM2bM0NKlS/XVV185sx4AAAAAAKzsvqc1Pj5et912m7p06aLw8HCFhITIw8PDZozBYNCKFSscXiQAAAAAoGKyO7QePnxYBoNBwcHBys7O1k8//ZRvjMFgcGhxAAAAAICKze7QevDgQWfWAQAAAABAPnbf0woAAAAAQGmz+0xramqqXeNCQkJuuhgAcJTcXINyc1x/y4JnVYPyMp3z/qDZ7JTdAgAAuBW7Q2vjxo3tumf14sWLt1QQADhCbo5Bp065ugrpNi/pspPqCAhwzn4BAADcid2hdcGCBflCq8lk0i+//KKPPvpId955pwYOHOjwAgEAAAAAFZfdobVv376Frnv++efVrl07Xb161SFFAQAAAAAgOWgipqpVq6pv375auHChI3YHAAAAAIAkB84e7OXlpTNnzjhqdwAAAAAAOCa0Hjx4UG+88Ybq16/viN0BAAAAACDJAbMHX7p0SZcvX1bVqlWVmJjo0OIAAAAAABWb3aG1devW+UKrwWCQn5+fwsLC1L17d/n5+Tm6PgAAAABABWZ3aF20aJEz6wAAAAAAIB+HTcQEAAAAAICjEVoBAAAAAG6L0AoAAAAAcFuEVgAAAACA23JpaE1OTlavXr10zz33yM/PT0uWLCly/C+//CI/P798X1u2bCmligEAAAAApcmu2YOzsrI0b948/fWvf1W7du0c9uSZmZmKjIxU7969NXToULu3W7lypRo1amR9XKNGDYfVBAAAAABwH3adaa1UqZL+/e9/69SpUw598kcffVSTJk1S586dZTTaf9L39ttvl7+/v/XL29vboXUBAAAAANyD3UkxKipKx44dc2Ytduvfv78iIiLUoUMHrVmzxtXlAAAAAACcxK7LgyVp0qRJevrpp3X//ferQ4cOzqypUFWrVtW0adPUsmVLeXp6asOGDRowYIAWLVqk2NjYQrc7evRoKVbpXOXpWMoS+l766vhL2VnZN729yeQts9n1c81ZLBaZzRYn7dvglH3fuE9nHkNxHHWMZpNZ2dk3//PkTJcvX1P16rzOuBK9d42i+l6vXr1SrASAu7M7tM6fP19+fn7q3bu3ateurTp16qhy5co2YwwGg1asWOHwIq+rWbOm4uLirI/vvfdeXbx4UfPmzSsytJaXF76jR4+Wm2MpS+i7a+ReTpVPJZ+b3t5kMqoEdx04jcFgkNFocNK+5fB9m82WfPt05jEUx1HHaPQwysfn5n+enOm22/74VczrjGvwGu8a9B1ASdgdWg8fPiyDwaDg4GBJ0smTJ/ONMRhK/4+aZs2aFTvrMAAAAACgbLI7tB48eNCZddy0gwcPyt/f39VlAAAAAACcwO7Q6gxXr161Tu5kNpt16tQpHThwQDVq1FBISIimTJmivXv3au3atZKkpUuXysvLS40bN5bRaNSmTZuUlJSkl19+2YVHAQAAAABwlhKFVpPJpJUrV+qrr77SuXPnNHHiRDVq1EgZGRnatm2b7r//fgUEBNi9v2+//VYxMTHWx/Hx8YqPj1fv3r21aNEipaWl6fjx4zbbJCQkKDU1VR4eHgoPD9eCBQuKvJ8VAAAAAFB22R1aL126pG7dumnfvn2qWrWqMjMzNXz4cElStWrVNGHCBPXq1UuTJk2y+8kffPBBZWRkFLp+0aJFNo/79OmjPn362L1/AAAAAEDZZvfcmlOmTNHhw4f18ccfa//+/bJY/u8jCDw8PBQTE6MvvvjCKUUCAAAAAComu0PrZ599psGDB+tvf/tbgbMEh4eHKzU11aHFAQAAAAAqNrtDa0ZGhurWrVvoeovFopycHIcUBQAAAACAVILQGhoaqkOHDhW6Pjk5WREREQ4pCgAAAAAAqQShtWfPnnr//feVnJxsXXb9MuHFixdr/fr1TJIEAAAAAHAou2cPHjFihPbs2aNOnTopIiJCBoNBY8eO1cWLF5Wenq7HH39cQ4YMcWatAAAAAIAKxu7Q6uXlpRUrVujjjz/Wp59+KoPBoLy8PDVp0kTdunXTk08+WeAETQAAAAAA3Cy7Q+t1PXv2VM+ePZ1RCwAAAAAANkocWiXp+++/t368TUhIiBo2bMhZVgAAAACAw5UotK5cuVKTJ0/Wr7/+KovFIumPyZhq166tyZMncwYWAAAAAOBQdofWJUuW6LnnnlO9evU0ZcoURUREyGKx6Oeff9b777+vIUOGKCcnR3379nVmvQAAAACACsTu0Dp37lw1a9ZM69evV6VKlWzWDRo0SNHR0Zo7dy6hFQAAAADgMHZ/Tuvp06fVs2fPfIFVkipVqqTY2Fj9+uuvDi0OAAAAAFCx2R1aGzRooDNnzhS6/tdff1X9+vUdUhQAAAAAAFIJQuvUqVP13nvvafXq1fnWrVy5Uu+//76mTZvm0OIAAAAAABWb3fe0vv7666pZs6aeeeYZjR07VnXr1pXBYNCxY8d07tw5hYeHa/78+Zo/f751G4PBoBUrVjilcAAAAABA+Wd3aD18+LAMBoOCg4MlyXr/qo+Pj4KDg5Wdna0jR47YbMNntwIAAAAAboXdofXgwYPOrAMAAAAAgHzsvqcVAAAAAIDSRmgFAAAAALgtQisAAAAAwG0RWgEAAAAAbovQCgAAAABwW4RWAAAAAIDbsju0NmnSRBs2bCh0/aZNm9SkSROHFAUAAAAAgFSCz2k9efKkMjMzC12fmZmp1NRUhxQF4NZ4G8/IaDrt6jJuiaePWcpxdRUAAABwNbtDqyQZDIZC1/3000+qVq3aLRcE4NYZTaflcWGUq8u4JZbbynb9AAAAcIwiQ+vSpUu1bNky6+OEhAS99957+cZlZGTo0KFD6tChg+MrBAAAAABUWEWG1szMTKWnp1sfX7p0SWaz2WaMwWBQlSpV9PTTT2vs2LHOqRIAAAAAUCEVGVoHDRqkQYMGSZIaN26smTNnKjo6ulQKAwAAAADA7ntaDxw44Mw6AAAAAADIp0QTMUnSlStXdOrUKf3222+yWCz51rdu3dohhQEAAAAAYHdo/e233zRmzBitXr1aJpMp33qLxSKDwaCLFy86tEAAAAAAQMVld2gdMWKE1q9fr0GDBql169by8/NzYlkAAAAAAJQgtG7ZskVDhgzRjBkznFkPAAAAAABWRnsHent7Kzw83Jm1AAAAAABgw+7Q2rlzZ33xxRfOrAUAAAAAABt2h9a4uDilpaVp6NCh+u9//6u0tDSdO3cu3xcAAAAAAI5i9z2tzZo1k8Fg0P79+7VixYpCxzF7MAAAAADAUewOraNHj5bBYHBmLQAAAAAA2LA7tI4bN86ZdQAAAAAAkI/d97T+mclk0sWLF5WXl+foegAAAAAAsCpRaN23b5+6dOmi2rVrKyIiQsnJyZKkCxcu6Mknn9T27dudUiQAAAAAoGKyO7Tu3r1b0dHROn78uHr16iWLxWJdV7NmTV29elUffPCBU4oEAAAAAFRMdofWadOmKTw8XLt27dKkSZPyrX/wwQe1Z88ehxYHAAAAAKjY7A6t+/btU79+/VSpUqUCZxEOCgpSenq6Q4sDAAAAAFRsdodWo9Eoo7Hw4enp6apcubJDiipKcnKyevXqpXvuuUd+fn5asmSJ058TAAAAAOAadofWpk2batOmTQWuy8nJ0ccff6z77rvPYYUVJjMzU5GRkZo5c2aphGQAAAAAgOvYHVpfeOEFffXVV3ruued08OBBSVJaWpq2bNmiTp066fjx4xo5cqTTCr3u0Ucf1aRJk9S5c+ciz/wCAAAAAMo+T3sHtm3bVosXL9aLL76opUuXSpKGDRsmi8Wi6tWrKykpSX/961+dVigAAAAAoOKxO7RKUo8ePRQdHa1t27bp559/ltlsVt26ddW+fXtVrVrVWTXesqNHj7q6BIcpT8dSlpS1vtfxv6zKWdmuLuOWeNwmZd/CMZhM3jKbXX81hsVikdlsKX7gTe3b4JR937hPZx5DcRx1jGaTWdnZ7vl/4vLla6pevey9zpQn9N41iup7vXr1SrESAO6uRKFVkqpUqaLHH3/cGbU4TXl54Tt69Gi5OZaypCz2vZLlkjxyfFxdxi3Jk+RT6eaPwWQyyh3uIDAYDDIa88+47ph9y+H7Npst+fbpzGMojqOO0ehhlI+Pe/6fuO22P34Vl7XXmfKiLL7Glwf0HUBJ2P0n3YYNG/Tiiy8Wuv7FF18sdKImAAAAAABuht2h9fXXX9fvv/9e6PqsrCzNmzfPIUUBAAAAACCVILQeOnRITZs2LXR9kyZNdPjwYUfUVKSrV6/qwIEDOnDggMxms06dOqUDBw4oNTXV6c8NAAAAAChddofWvLw8Xbt2rdD1165dK5VJLr799lu1adNGbdq00bVr1xQfH682bdrolVdecfpzAwAAAABKl90TMUVGRmrt2rV67rnn8n0+qtls1tq1a9WgQQOHF3ijBx98UBkZGU5/HgAAAACA69l9pnXo0KHau3evevfurf379ys7O1vZ2dnav3+/+vTpo71792rIkCHOrBUAAAAAUMHYfaa1e/fuOn78uOLj4/XFF19I+uNjECwWiwwGg8aMGaPY2FinFQoAAAAAqHhK9Dmto0aNUo8ePbRu3TqdOHFCFotFdevWVUxMjOrUqeOkEgEAAAAAFZVdofXatWt68sknFRsbq379+ikuLs7ZdQEAAAAAYN89rZUrV9Z3330nk8nk7HoAAAAAALCyeyKmBx54QF9//bUzawEAAAAAwIbdoXXWrFnat2+fJk6cqBMnTshsNjuzLgAAAAAA7J+I6a9//assFosSExOVmJgoo9EoLy8vmzEGg0G//vqrw4sEAAAAAFRMdofWrl27ymAwOLMWAAAAAABs2B1aFy1a5Mw6AAAAAADIx+57WgEAAAAAKG0lCq0nT57UP//5TzVt2lQhISHauXOnJOnChQsaOXKk9u/f74waAQAAAAAVlN2XBx85ckQdO3aU2WxW8+bNdfLkSevnttasWVP//e9/lZ2drQULFjitWAAAAABAxWJ3aJ08ebKqVaumLVu2yMPDQxERETbrH330UX366aeOrg8AAAAAUIHZfXnw119/rYEDB+rOO+8scBbhkJAQnTlzxqHFAQAAAAAqNrtDa15ennx9fQtd/9tvv8nDw8MhRQEAAAAAIJUgtEZGRmrHjh0FrrNYLFq3bp2aNm3qqLoAAAAAALA/tA4bNkxr1qzR7NmzdfHiRUmS2WzWjz/+qH/84x/69ttvFRcX57RCAQAAAAAVj90TMXXv3l2pqamaMWOGZs6caV0mSR4eHpo+fboeeeQR51QJAAAAAKiQ7A6tkvT888+rR48eWrt2rY4dOyaz2ay6deuqU6dOuuuuu5xVIwAA5ZbB4KHU1Dq6dKmSq0txmqAgswIDc1xdBgCgjCo2tGZnZ2vDhg06ceKEbr/9dnXo0EHDhw8vjdoAACj3zp+XXn65snx8yu9khgkJUmCgq6sAAJRVRYbW9PR0RUdH6/jx47JYLJIkX19fLV++XK1bty6VAgEAAAAAFVeREzFNnz5dJ06c0PDhw7V8+XLFx8fLx8dHo0ePLq36AAAAAAAVWJFnWv/zn/+od+/emj59unXZnXfeqYEDB+r06dMKCgpyeoEAAAAAgIqryDOt6enpatGihc2yli1bymKx6NSpU04tDAAAAACAIkOryWRSpUq2sxlef5yVleW8qgAAAAAAkB2zB584cUJ79+61Pr58+bIk6ejRo6patWq+8c2aNXNgeQAAAACAiqzY0BofH6/4+Ph8y2+cjMlischgMOjixYuOqw4AAAAAUKEVGVoTExNLqw4AAAAAAPIpMrT26dOntOoAAMBp6kVYNG38TleXUaBGjaSZk8wyehQ5zYQk6VRaiBa/c1cpVAUAgPso9vJgAADKOt9K5xXmO9PVZRTotmtSWFWLjEZD8YMDXpVEaAUAVCzFv60LAAAAAICLEFoBAAAAAG6L0AoAAAAAcFuEVgAAAACA2yK0AgAAAADcFqEVAAAAAOC2CK0AAAAAALdFaAUAAAAAuC1CKwAAAADAbXm6ugAApS8316DcHIOryyiSsYqncjJv/n01s9mBxQAAAMBlCK1ABZSbY9CpU66uomjVPA26cgs1BgQ4rhYAAAC4DpcHAwAAAADcFqEVAAAAAOC2CK0AAAAAALdFaAUAAAAAuC2XT8SUlJSk+fPnKz09XQ0aNFB8fLxatWpV4NhffvlFTZo0ybf8k08+0d/+9jdnlwoAgEvVi7Bo2vidri6jxCLrWFTJ8seU3maPIOWYA11cEQCgLHFpaF21apXGjh2rV199VS1btlRSUpJ69uyplJQUhYSEFLrdypUr1ahRI+vjGjVqlEa5AAC4lG+l8wrznenqMkrMN1Py0P//OVQ1EyQDoRUAYD+XhtbExET16dNHTz/9tCRpzpw52rp1q9555x1Nnjy50O1uv/12+fv7l1aZqEC8jWdkNJ22WVbH/7IqWS65qKKb42G45uoSAAAAAIdwWWjNycnR/v37FRcXZ7O8Xbt22rVrV5Hb9u/fX1lZWQoPD9fw4cPVuXNnZ5aKCsRoOi2PC6NsllXOypZHjo+LKrpJNce5ugIAAADAIVwWWi9cuCCTyaRatWrZLK9Vq5bOnj1b4DZVq1bVtGnT1LJlS3l6emrDhg0aMGCAFi1apNjY2EKf6+jRow6t3ZXK07G4ozr+l1U5Kzvf8uwClrkzj9xcmYqo2WTyltns/vOwmc2Wm97WYjHc0vaOYrFYnFaHs47xxn068xiK46hjdOUxFMdiMUiy7+fdnY+jKCaTWdlZOZKka1cu60S6e/0u43eraxTV93r16pViJQDcncsnYjIYDDaPLRZLvmXX1axZ0+bM7L333quLFy9q3rx5RYbW8vLCd/To0XJzLO6qkuVSvrOq2VnZ8qlUxs60ennJs4iaTSajjO6fWWU0FvxaYA+D4da2dxSDweC0OpxxjGazJd8+nXkMxXHUMbryGIpz/VeePfW583EUxcPDw/o66lntNtW7zX1+l/G71TXoO4CScNmfrTVr1pSHh0e+s6rnz5/Pd/a1KM2aNdOxY8ccXR4AAAAAwA24LLR6e3uradOm2rZtm83ybdu2qUWLFnbv5+DBg0zKBAAAAADllEsvD3722Wc1ZMgQNWvWTC1atNA777yjtLQ0DRgwQJI0ZcoU7d27V2vXrpUkLV26VF5eXmrcuLGMRqM2bdqkpKQkvfzyyy48CgAAAACAs7g0tHbr1k0XL17UnDlzlJ6ernvuuUcrVqxQaGioJCktLU3Hjx+32SYhIUGpqany8PBQeHi4FixYUOT9rAAAAACAssvlEzENHDhQAwcOLHDdokWLbB736dNHffr0KY2yAAAAAABuoAzMHwoAAAAAqKgIrQAAAAAAt0VoBQAAAAC4LUIrAAAAAMBtEVoBAAAAAG6L0AoAAAAAcFuEVgAAAACA2yK0AgAAAADcFqEVAAAAAOC2CK0AAAAAALdFaAUAAAAAuC1CKwAAAADAbRFaAQAAAABui9AKAAAAAHBbhFYAAAAAgNsitAIAAAAA3BahFQAAAADgtjxdXQAAACj/fs/8433yTBl16EQlF1fzfy5frqNLl269nqAgswIDcxxQEQDgRoRWAADgVHl5UlraH/8+lmnQxFc8XFvQn2RnV5aPz63Xk5AgBQY6oCAAQD5cHgwAAAAAcFuEVgAAAACA2yK0AgAAAADcFve0wmG8jWdkNJ12dRm3xMNwzdUlAAAAAPgTQiscxmg6LY8Lo1xdxq2pOc7VFQAAAAD4Ey4PBgAAAAC4LUIrAAAAAMBtEVoBAAAAAG6Le1qBG+TmGpSbY7A+Npm8ZTKVrfd3PKsalJdZeM1mcykWAwB/Ui/Comnjd7q6DCuzySyjR8lf40+lhWjxO3c5oSIAwI0IrcANcnMMOnXq/x6bzUYZy1Zm1W1e0uVTha8PCCi9WgDgz3wrnVeY70xXl2FlNltkNBqKH3ijgFclEVoBoDSUsT/FAQAAAAAVCaEVAAAAAOC2CK0AAAAAALdFaAUAAAAAuC1CKwAAAADAbRFaAQAAAABui9AKAAAAAHBbfE4rAADALTIYPLRnTyVXl+FUQUFmBQbmuLoMABUQoRUAAOAWnT8vxcd7uLoMp0pIkAIDXV0FgIqIy4MBAAAAAG6L0AoAAAAAcFuEVgAAAACA2+KeVpTImTPeOn264Pc6IusY5ZtZtt8H8axqkNns6ioAAAAAXEdoRYmcPm3UqFEFTzQxbbxBYb6lXJCD3eYlVSG0AgAAAG6D0OomvI1nZDSdLnJMHf/LqmS5VEoVFSyyjlHTxhsKXHd3+O/KSyvlggAAAACUa4RWN2E0nZbHhVFFjqmclS2PHJ9SqqhgvpnGQs+mVqk8VpdLtxwAAAAA5Ryh1YGKut+zOPbcD2oyectkcu09o9zvCQAAAKA0lcnQmpSUpPnz5ys9PV0NGjRQfHy8WrVq5eqyirzfszj23A9qNhtldPE8RwEBrn1+AAAAABVLmZvqddWqVRo7dqxGjhypr776Svfdd5969uyp1NRUV5cGAAAAAHCwMhdaExMT1adPHz399NOqX7++5syZI39/f73zzjuuLg0AAAAA4GBlKrTm5ORo//79ateunc3ydu3aadeuXS6qqvQYjQXP2gvnou+uYTDQd1fg59016Lvr0HvXqFevnqtLAFCGGDIyMiyuLsJeZ86c0T333KPPPvtMrVu3ti6fNWuWPv74Y+3Zs8eF1QEAAAAAHK1MnWm97sYzMBaLhbMyAAAAAFAOlanQWrNmTXl4eOjs2bM2y8+fP69atWq5qCoAAAAAgLOUqdDq7e2tpk2batu2bTbLt23bphYtWrioKgAAAACAs5S5z2l99tlnNWTIEDVr1kwtWrTQO++8o7S0NA0YMMDVpQEAAAAAHKxMnWmVpG7duik+Pl5z5szRgw8+qJSUFK1YsUKhoaGuLs0h0tLSNHToUIWHh8vf318tWrTQzp07restFovi4+PVoEEDBQQE6PHHH9cPP/zgworLPpPJpOnTp6tx48by9/dX48aNNX36dOXl5VnH0Pdbl5ycrF69eumee+6Rn5+flixZYrPenh5nZ2frxRdfVFhYmGrXrq1evXrp9OnTpXkYZVJRvc/NzdXkyZPVqlUr1a5dW/Xr19fAgQPzffY1vS+54n7m/+xf//qX/Pz89Prrr9ssp+8lZ0/ff/rpJ/Xr10+hoaEKDAxUmzZtdOTIEet6+l5yxfX96tWrevHFFxUZGamAgAA1b95ciYmJNmPoO4DClLnQKkkDBw7UwYMHdfbsWW3fvt1mJuGyLCMjQx06dJDFYtGKFSu0a9cuzZ492+Z+3Xnz5ikxMVGzZs3Sf/7zH9WqVUtdu3bVlStXXFh52fbaa68pKSlJs2bN0u7duzVz5ky99dZbmjt3rnUMfb91mZmZioyM1MyZM1W5cuV86+3p8bhx47Ru3Tq9/fbb2rBhg65cuaLY2FiZTKbSPJQyp6je//777/ruu+80atQobd++XUuXLtXp06fVo0cPmzdu6H3JFfczf92aNWu0b98+BQYG5ltH30uuuL6fOHFCHTp00F133aW1a9fqm2++0UsvvSRfX1/rGPpecsX1fcKECfr888/1xhtvaNeuXRo5cqSmTJmijz76yDqGvgMoTJn6yJvyburUqUpOTtbmzZsLXG+xWNSgQQMNGjRIo0aNkiRdu3ZN9erV07Rp07hE+ibFxsaqRo0aeuONN6zLhg4dqt9++03Lly+n704QFBSk2bNnq2/fvpLs+9m+dOmSIiIilJiYqCeffFKSdOrUKUVFRemTTz5R+/btXXY8ZcmNvS/I4cOH1bJlSyUnJ6thw4b03gEK6/vJkyfVoUMHffrpp+rRo4cGDx6suLg4SaLvDlBQ3wcOHCiDwaC33nqrwG3o+60rqO/333+/YmJiNH78eOuy6OhoNWzYUHPmzKHvAIpUJs+0llefffaZmjVrpgEDBigiIkIPPPCA3nzzTVksf7yv8Msvvyg9PV3t2rWzblO5cmW1atVKu3btclXZZV7Lli21c+dO/fjjj5L++IN9x44deuSRRyTR99JgT4/379+v3NxcmzHBwcGqX78+3wcHu35228/PTxK9d5a8vDwNHDhQo0aNUv369fOtp++OZzabtWnTJtWvX1/du3dXeHi42rZtq1WrVlnH0HfnaNmypTZt2qRTp05Jknbt2qXvv//eGkbpO4CilLmJmMqzEydO6O2339bw4cP1/PPP6+DBgxozZowkafDgwUpPT5ekfB/vU6tWLZ05c6bU6y0vnn/+eV29elUtWrSQh4eH8vLyNGrUKA0cOFCS6HspsKfHZ8+elYeHh2rWrJlvzI0fg4Wbl5OTo5deekkdO3ZUUFCQJHrvLPHx8apRo4aeeeaZAtfTd8c7d+6crl69qrlz52r8+PGaPHmyvvrqKw0aNEhVqlRRx44d6buTzJo1SyNGjFCjRo3k6fnHn5+zZ89Wx44dJfHzDqBohFY3Yjabde+992ry5MmSpCZNmujYsWNKSkrS4MGDreMMBoPNdhaLJd8y2G/VqlX66KOPlJSUpAYNGujgwYMaO3asQkND9dRTT1nH0Xfnu5ke831wnLy8PA0ePFiXLl3SsmXLih1P72/ezp07tXTpUu3YsaPE29L3m2c2myX9cVnqc889J0lq3Lix9u/fr6SkJGuAKgh9vzWLFy/Wrl27tGzZMoWEhOjrr7/WxIkTFRoaqr/97W+FbkffAUhcHuxW/P39810idvfdd1svpfH395ekfO84nj9/Pt8ZKthv0qRJeu6559S9e3c1bNhQvXr10rPPPqt///vfkuh7abCnx3feeadMJpMuXLhQ6BjcvLy8PD3zzDP63//+pzVr1uj222+3rqP3jrdjxw6lpaWpfv36qlmzpmrWrKnU1FRNnjxZkZGRkui7M9SsWVOenp5F/q6l74537do1TZ06VVOmTNFjjz2mRo0aafDgwerWrZt1xmz6DqAohFY30rJlS/300082y3766SeFhIRIku666y75+/tr27Zt1vVZWVn65ptv1KJFi1KttTz5/fff5eHhYbPMw8PD+o48fXc+e3rctGlTeXl52Yw5ffq0jhw5wvfhFuXm5mrAgAH63//+p3Xr1lnfRLiO3jvewIEDlZycrB07dli/AgMDNXz4cK1Zs0YSfXcGb29v/eUvf9HRo0dtlv/5dy19d7zc3Fzl5uYW+buWvgMoCpcHu5Hhw4fr0UcfVUJCgrp166YDBw7ozTff1MSJEyX9cenksGHD9Oqrr6pevXqKiIhQQkKCfH191aNHDxdXX3Z17NhRr732mu666y41aNBABw4cUGJionr16iWJvjvK1atXdezYMUl/XKJ36tQpHThwQDVq1FBISEixPa5evbr69++vSZMmqVatWqpRo4YmTJighg0b6uGHH3bhkbm/onofGBiop59+Wt9++62WLVsmg8Fgvcf4tttuU+XKlen9TSruZ/7Gs0eenp7y9/dXvXr1JPEzf7OK6/s///lPDRgwQK1atVKbNm20Y8cOrVq1yvq5ovT95hTX99atW2vKlCny9fVVSEiIkpOT9dFHH2nKlCmS6DuAovGRN25m8+bNmjp1qn766ScFBwdr0KBBGjJkiPV+DovFopkzZ+r//b//p4yMDDVr1kwJCQnWy8lQcleuXNGMGTO0fv16nT9/Xv7+/urevbtGjx6tSpUqSaLvjrBjxw7FxMTkW967d28tWrTIrh5nZWVp4sSJ+uSTT5SVlaU2bdro1VdfVXBwcGkeSplTVO/Hjh2rJk2aFLhdYmKi9SMr6H3JFfczf6OoqCibj7yR6PvNsKfvS5Ys0dy5c3X69GmFhYXphRdesHkTkr6XXHF9T09P15QpU7Rt2zb99ttvCgkJ0VNPPaXnnnvO+jcOfQdQGEIrAAAAAMBtcU8rAAAAAMBtEVoBAAAAAG6L0AoAAAAAcFuEVgAAAACA2yK0AgAAAADcFqEVAAAAAOC2CK0AAAAAALdFaAUAAAAAuC1CKwAAAADAbf1/paZhZxPRabIAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"unit = ''\n",
"\n",
"fig, ax = plt.subplots(figsize=(10,5))\n",
"\n",
"ax.hist(smoker['Birth Weight'], density=True, \n",
" label='Maternal Smoker = True', \n",
" color='blue', \n",
" alpha=0.8, \n",
" ec='white', \n",
" zorder=5)\n",
"\n",
"ax.hist(non_smoker['Birth Weight'], density=True, \n",
" label='Maternal Smoker = ', \n",
" color='gold', \n",
" alpha=0.8, \n",
" ec='white', \n",
" zorder=10)\n",
"\n",
"y_vals = ax.get_yticks()\n",
"\n",
"y_label = 'Percent per ' + (unit if unit else 'unit')\n",
"\n",
"x_label = ''\n",
"\n",
"ax.set_yticklabels(['{:g}'.format(x * 100) for x in y_vals])\n",
"\n",
"plt.ylabel(y_label)\n",
"\n",
"plt.xlabel(x_label)\n",
"\n",
"plt.title('')\n",
"\n",
"ax.legend(bbox_to_anchor=(1.04,1), loc=\"upper left\")\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The distribution of the weights of the babies born to mothers who smoked appears to be based slightly to the left of the distribution corresponding to non-smoking mothers. The weights of the babies of the mothers who smoked seem lower on average than the weights of the babies of the non-smokers. \n",
"\n",
"This raises the question of whether the difference reflects just chance variation or a difference in the distributions in the larger population. Could it be that there is no difference between the two distributions in the population, but we are seeing a difference in the samples just because of the mothers who happened to be selected?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Hypotheses ###\n",
"We can try to answer this question by a test of hypotheses. The chance model that we will test says that there is no underlying difference in the popuations; the distributions in the samples are different just due to chance. \n",
"\n",
"Formally, this is the null hypothesis. We are going to have to figure out how to simulate a useful statistic under this hypothesis. But as a start, let's just state the two natural hypotheses.\n",
"\n",
"**Null hypothesis:** In the population, the distribution of birth weights of babies is the same for mothers who don't smoke as for mothers who do. The difference in the sample is due to chance.\n",
"\n",
"**Alternative hypothesis:** In the population, the babies of the mothers who smoke have a lower birth weight, on average, than the babies of the non-smokers."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Test Statistic ###\n",
"The alternative hypothesis compares the average birth weights of the two groups and says that the average for the mothers who smoke is smaller. Therefore it is reasonable for us to use the difference between the two group means as our statistic. \n",
"\n",
"We will do the subtraction in the order \"average weight of the smoking group $-$ average weight of the non-smoking group\". Small values (that is, large negative values) of this statistic will favor the alternative hypothesis. \n",
"\n",
"The observed value of the test statistic is about $-9.27$ ounces."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Maternal Smoker
\n",
"
Birth Weight
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
False
\n",
"
123.085315
\n",
"
\n",
"
\n",
"
1
\n",
"
True
\n",
"
113.819172
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Maternal Smoker Birth Weight\n",
"0 False 123.085315\n",
"1 True 113.819172"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"means_table = smoking_and_birthweight.groupby('Maternal Smoker').mean()\n",
"\n",
"means_table = means_table.reset_index()\n",
"\n",
"means_table"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-9.266142572024918"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"means = means_table['Birth Weight']\n",
"observed_difference = means[1] - means[0]\n",
"observed_difference"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are going compute such differences repeatedly in our simulations below, so we will define a function to do the job. The function takes three arguments:\n",
"\n",
"- the name of the table of data\n",
"- the label of the column that contains the numerical variable whose average is of interest\n",
"- the label of the column that contains the Boolean variable for grouping\n",
"\n",
"It returns the difference between the means of the `True` group and the `False` group."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"def difference_of_means(table, label, group_label):\n",
" reduced = table[[label, group_label]]\n",
" means_table = reduced.groupby(group_label).mean()\n",
" means = means_table[label]\n",
" return means[1] - means[0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To check that the function is working, let's use it to calculate the observed difference between the means of the two groups in the sample."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-9.266142572024918"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"difference_of_means(births, 'Birth Weight', 'Maternal Smoker')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That's the same as the value of `observed_difference` calculated earlier."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Predicting the Statistic Under the Null Hypothesis\n",
"\n",
"To see how the statistic should vary under the null hypothesis, we have to figure out how to simulate the statistic under that hypothesis. A clever method based on *random permutations* does just that.\n",
"\n",
"If there were no difference between the two distributions in the underlying population, then whether a birth weight has the label `True` or `False` with respect to maternal smoking should make no difference to the average. The idea, then, is to shuffle all the labels randomly among the mothers. This is called *random permutation*. \n",
"\n",
"Take the difference of the two new group means: the mean weight of the babies whose mothers have been randomly labeled smokers and the mean weight of the babies of the remaining mothers who have all been randomly labeled non-smokers. This is a simulated value of the test statistic under the null hypothesis.\n",
"\n",
"Let's see how to do this. It's always a good idea to start with the data."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Maternal Smoker
\n",
"
Birth Weight
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
False
\n",
"
120
\n",
"
\n",
"
\n",
"
1
\n",
"
False
\n",
"
113
\n",
"
\n",
"
\n",
"
2
\n",
"
True
\n",
"
128
\n",
"
\n",
"
\n",
"
3
\n",
"
True
\n",
"
108
\n",
"
\n",
"
\n",
"
4
\n",
"
False
\n",
"
136
\n",
"
\n",
"
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
\n",
"
\n",
"
1169
\n",
"
False
\n",
"
113
\n",
"
\n",
"
\n",
"
1170
\n",
"
False
\n",
"
128
\n",
"
\n",
"
\n",
"
1171
\n",
"
True
\n",
"
130
\n",
"
\n",
"
\n",
"
1172
\n",
"
False
\n",
"
125
\n",
"
\n",
"
\n",
"
1173
\n",
"
False
\n",
"
117
\n",
"
\n",
" \n",
"
\n",
"
1174 rows × 2 columns
\n",
"
"
],
"text/plain": [
" Maternal Smoker Birth Weight\n",
"0 False 120\n",
"1 False 113\n",
"2 True 128\n",
"3 True 108\n",
"4 False 136\n",
"... ... ...\n",
"1169 False 113\n",
"1170 False 128\n",
"1171 True 130\n",
"1172 False 125\n",
"1173 False 117\n",
"\n",
"[1174 rows x 2 columns]"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"smoking_and_birthweight"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are 1,174 rows in the table. To shuffle all the labels, we will draw a random sample of 1,174 rows without replacement. Then the sample will include all the rows of the table, in random order. \n",
"\n",
"We can use the Table method `sample` with the optional `with_replacement=False` argument. We don't have to specify a sample size, because by default, `sample` draws as many times as there are rows in the table."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"smoking_and_birthweight2 = smoking_and_birthweight.copy()\n",
"\n",
"shuffled_labels = smoking_and_birthweight2.sample(len(smoking_and_birthweight2), replace = False)\n",
"\n",
"shuffled_labels = np.array(shuffled_labels['Maternal Smoker'])\n",
"\n",
"smoking_and_birthweight2['Shuffled Label'] = shuffled_labels\n",
"\n",
"original_and_shuffled = smoking_and_birthweight2"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Maternal Smoker
\n",
"
Birth Weight
\n",
"
Shuffled Label
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
False
\n",
"
120
\n",
"
False
\n",
"
\n",
"
\n",
"
1
\n",
"
False
\n",
"
113
\n",
"
False
\n",
"
\n",
"
\n",
"
2
\n",
"
True
\n",
"
128
\n",
"
False
\n",
"
\n",
"
\n",
"
3
\n",
"
True
\n",
"
108
\n",
"
False
\n",
"
\n",
"
\n",
"
4
\n",
"
False
\n",
"
136
\n",
"
False
\n",
"
\n",
"
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
\n",
"
\n",
"
1169
\n",
"
False
\n",
"
113
\n",
"
False
\n",
"
\n",
"
\n",
"
1170
\n",
"
False
\n",
"
128
\n",
"
False
\n",
"
\n",
"
\n",
"
1171
\n",
"
True
\n",
"
130
\n",
"
False
\n",
"
\n",
"
\n",
"
1172
\n",
"
False
\n",
"
125
\n",
"
False
\n",
"
\n",
"
\n",
"
1173
\n",
"
False
\n",
"
117
\n",
"
False
\n",
"
\n",
" \n",
"
\n",
"
1174 rows × 3 columns
\n",
"
"
],
"text/plain": [
" Maternal Smoker Birth Weight Shuffled Label\n",
"0 False 120 False\n",
"1 False 113 False\n",
"2 True 128 False\n",
"3 True 108 False\n",
"4 False 136 False\n",
"... ... ... ...\n",
"1169 False 113 False\n",
"1170 False 128 False\n",
"1171 True 130 False\n",
"1172 False 125 False\n",
"1173 False 117 False\n",
"\n",
"[1174 rows x 3 columns]"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"original_and_shuffled"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Each baby's mother now has a random smoker/non-smoker label in the column `Shuffled Label`, while her original label is in `Maternal Smoker`. If the null hypothesis is true, all the random re-arrangements of the labels should be equally likely.\n",
"\n",
"Let's see how different the average weights are in the two randomly labeled groups."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Shuffled Label
\n",
"
Birth Weight
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
False
\n",
"
119.019580
\n",
"
\n",
"
\n",
"
1
\n",
"
True
\n",
"
120.152505
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Shuffled Label Birth Weight\n",
"0 False 119.019580\n",
"1 True 120.152505"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"shuffled_only = original_and_shuffled.drop(columns=['Maternal Smoker'])\n",
"\n",
"shuffled_group_means = shuffled_only.groupby('Shuffled Label').mean()\n",
"\n",
"shuffled = shuffled_group_means.reset_index()\n",
"\n",
"shuffled"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The averages of the two randomly selected groups are quite a bit closer than the averages of the two original groups. We can use our function `difference_of_means` to find the two differences."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.132925027042674"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"difference_of_means(original_and_shuffled, 'Birth Weight', 'Shuffled Label')"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-9.266142572024918"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"difference_of_means(original_and_shuffled, 'Birth Weight', 'Maternal Smoker')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But could a different shuffle have resulted in a larger difference between the group averages? To get a sense of the variability, we must simulate the difference many times. \n",
"\n",
"As always, we will start by defining a function that simulates one value of the test statistic under the null hypothesis. This is just a matter of collecting the code that we wrote above. But because we will later want to use the same process for comparing means of other variables, we will define a function that takes three arguments:\n",
"\n",
"- the name of the table of data\n",
"- the label of the column that contains the numerical variable\n",
"- the label of the column that contains the Boolean variable for grouping\n",
"\n",
"It returns the difference between the means of two groups formed by randomly shuffling all the labels."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"def one_simulated_difference(table, label, group_label):\n",
" \n",
" births1 = table.copy()\n",
"\n",
" shuffled_labels = births1.sample(len(births1), replace = False)\n",
"\n",
" shuffled_labels = np.array(shuffled_labels[group_label])\n",
"\n",
" births1['Shuffled Label'] = shuffled_labels\n",
"\n",
" original_and_shuffled = births1\n",
"\n",
" shuffled_only = original_and_shuffled.drop(columns=['Maternal Smoker'])\n",
"\n",
" shuffled_group_means = shuffled_only.groupby('Shuffled Label').mean()\n",
"\n",
" table1 = shuffled_group_means.reset_index()\n",
" \n",
" return difference_of_means(table1, label, 'Shuffled Label') "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Run the cell below a few times to see how the output changes."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.5407315995551158"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"one_simulated_difference(births, 'Birth Weight', 'Maternal Smoker')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Permutation Test\n",
"Tests based on random permutations of the data are called *permutation tests*. We are performing one in this example. In the cell below, we will simulate our test statistic – the difference between the averages of the two groups – many times and collect the differences in an array. "
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0.18495361, 0.00251383, 1.04707101, ..., 0.77519996,\n",
" 1.81260265, -0.38025199])"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"differences = np.array([])\n",
"\n",
"repetitions = 5000\n",
"for i in np.arange(repetitions):\n",
" new_difference = one_simulated_difference(births, 'Birth Weight', 'Maternal Smoker')\n",
" differences = np.append(differences, new_difference)\n",
"\n",
"differences"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The array `differences` contains 5,000 simulated values of our test statistic: the difference between the mean weight in the smoking group and the mean weight in the non-smoking group, when the labels have been assigned at random. \n",
"\n",
"### Conclusion of the Test\n",
"The histogram below shows the distribution of these 5,000 values. It is the empirical distribution of the test statistic simulated under the null hypothesis. This is a prediction about the test statistic, based on the null hypothesis."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Observed Difference: -9.266142572024918\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAqwAAAFuCAYAAABECkoSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABSaklEQVR4nO3deVxN+f8H8NfVLnIxlUoLlaEYDVkmy1iztRAJGctYs42toWHGPoUwlhhjm4ydMrKMPbTQMGM3yJpIaBW5Uff3h1/367q33MvNPen1fDw8Zu45n3PO+3zuvfXqcz/nXFFmZqYUREREREQCVU7bBRARERERFYeBlYiIiIgEjYGViIiIiASNgZWIiIiIBI2BlYiIiIgEjYGViIiIiASNgZXoIwkICIBYLEZMTIzccrFYjC5dupTYcYODg5Uel14r7J+NGzdqu5QidenSBWKxGHfv3tV2KR+FsvcEX8fF27hxI8RiMYKDg7VdCgA+X6R5DKz0SRGLxXL/qlSpAltbW3Ts2BHr1q1Dfn6+tkvUuNIQuN5WWHNxv1xjYmJKPMwLhdADab169SAWi2FlZYWUlBSlbb799ltBBZRP7TUmtEBK9LHparsAopIwadIkAEB+fj5u376NPXv24NSpUzh27BjCw8O1XJ28v//+G0ZGRiW2/6FDh6J79+6oXr16iR2DyoZnz55h9uzZCAsL03YpJHD8uUOaxsBKn6SgoCC5x5cvX0a7du2wa9cuxMfHw83NTUuVKapVq1aJ7r9q1aqoWrVqiR6DygZ7e3ts3rwZw4YNwxdffKHtckjA+HOHNI1TAqhMcHZ2RrNmzQAA//zzD4D/fRwYEBCAq1evom/fvqhZsybEYjEuXLgg23bXrl3w9vaGnZ0dzMzM0KBBA0yfPh3Z2dlKj3Xs2DF06tQJlpaWsLOzQ58+fXDt2rUiayvqI8n8/HysX78enTp1gq2tLczNzfHFF19g8ODBOHv2LIDXHyXPnTsXADBy5Ei56RCFHy8XN5fsxIkT8PX1RY0aNWBmZob69etj0qRJePz4sULbN+fg7tq1C23atIGFhQXs7OwwcOBA3L9/v8hz1KT3rePcuXOyER9ra2t4e3sjISGh2GOlpqZi8uTJaNCgAczNzWFra4tu3brh+PHjCm3f/Mg2ISEBPj4+sLW1hVgsRmZmZpHHEIvFiIuLAwDUr19f9vzVq1dPaft169bBzc0N5ubmcHR0xJgxY4rcvzr1q2L69OkoKCjA1KlTVd6muI/chTjPcdWqVRCLxQgJCVG6Pjs7G5aWlnB2dpZNMXpzWs7+/fvRvn172ft/wIABuH37ttJ9paamIjAwEPXr14eZmRlq1KiBnj17IjY2Vq5dQEAARo4cCQCYO3eu3PtcWd9duHABPXv2hI2NDSwsLNCpUyecOnVKaQ0FBQVYv349OnToABsbG5ibm+Orr77CwoULkZeXp9A+JiYGfn5+cHZ2hpmZGRwcHNCqVStMmTIFUun/vum9qOdW1e2J3sYRVirzbt++DXd3d3z++efo1asXsrKyUL58eQDAhAkTsGbNGlhZWcHDwwNisRhnzpzBL7/8goMHD+LAgQOoWLGibF+7du3CwIEDoaenh65du8LS0hKnTp1C+/btUbduXZVrysvLQ58+fXD48GFUq1YN3bp1Q+XKlZGcnIyYmBjY29vjyy+/RJ8+fQAAcXFx6Ny5s1zIqVSpUrHHWLduHcaPHw8jIyN4e3ujWrVqSEhIwMqVK7F371789ddfsLa2VthuzZo1+Ouvv9C5c2c0a9YMZ86cwc6dO3Hx4kXExcXBwMBA5fP8EOrUkZCQgK5du0IikcDT0xP29va4fPkyPD090bJlS6X7v3z5Mrp164bHjx+jTZs26Ny5M9LT07F371507doVS5YswTfffKOw3d9//42FCxfCzc0N/fr1Q0pKCnR0dIo8j0mTJmHTpk24d+8ehg8fLnvelD1/06ZNw9GjR9GxY0e0bt0aMTExWL9+PW7cuIF9+/ZppP7itG/fHq1bt0Z0dDT279+Pjh07qrV9adCrVy/MnDkTf/zxBwIDAxWeuy1btuD58+cYM2aMwrrdu3fj8OHD8PT0RIsWLXDhwgX8+eefiImJwcGDB2Fvby9re/fuXXTq1AkPHjxAs2bN4OPjg4cPH+LPP//E4cOH8csvv6Bfv34AXv9hmpWVhX379qFZs2Zo3ry5bD82NjZyNZw7dw5LlixBkyZN0K9fPyQnJyMqKgre3t44ceIEPv/8c1nbV69eoW/fvti/fz8cHBzQvXt3GBgYIC4uDjNnzsTx48cREREBXd3XUeHgwYPw8/NDxYoV0alTJ1hZWSEzMxM3b97EypUrMWPGDFlbZT50eyrb+MqgMuG///6TjWI1aNBAbt2pU6cwfvx4/PTTT3LLt27dijVr1sDDwwOrVq2Sm2c6f/58zJkzB8HBwfj5558BADk5ORg7dixEIhH27t0LV1dXWfsff/wRS5cuVbneuXPn4vDhw2jVqhU2bdokC9DA65HXwhFQf39/JCUlIS4uDl26dIG/v79K+09KSsKkSZNQvnx5HD58GHXq1JGtmz17NkJDQzFhwgRs27ZNYdujR4/i+PHjqF27tmzZ4MGDsWPHDuzduxc+Pj4qn+eHULUOqVSKUaNGITc3F+Hh4fD29pa1X7VqFQIDAxX2nZ+fj/79+yMrKwu7d++WCwgPHz5E27ZtERgYiA4dOsDMzExu2+joaPzyyy8YMGCASucRFBSE2NhY3Lt3DwEBAbC1tS2y7T///IOTJ0/CysoKwOvA4enpifj4eJw5c0b2mvuQ+t9l1qxZaNmyJX766Se0a9dO8AEjNja2yAuVkpKSFJZVrFgRfn5+WL16Nfbv368wOvz7779DV1dXFibftH//fmzduhUdOnSQLVu6dCl+/PFHBAYGIjIyUrZ83LhxePDgASZPnozJkyfLlo8aNQrt2rVDYGAg2rRpg+rVq8PDw0MWWJs3b64w5elNBw4cwMqVK+Hn5ydbtm7dOowbNw4rV67EwoULZcsXLVqE/fv3Y8iQIQgJCZEF8IKCAowbNw7h4eFYvXo1hg8fDgBYv349pFIpdu/ejfr168sdNz09/Z2vhQ/dnso2TgmgT1JwcDCCg4Mxe/ZsDBkyBK1bt0Zubi48PDxkUwMKmZmZyS7SetPy5cuho6ODpUuXKlwUNX78eFStWlUu0O3btw8ZGRnw8fGRC6sA8P3338PExESl2vPz87F69WoYGBhg8eLFcmEVAHR0dFCtWjWV9lWUbdu2IS8vD4MGDZILqwAQGBgICwsLHDx4EA8ePFDYdtiwYXIhEQD69+8PAPj3338/qC51qFpHQkICEhMT0aRJE7mwCgCDBg1CzZo1FfZ98OBB3LhxA4MGDZILewBQrVo1jB49Gi9evMCuXbsUtq1bt67KYVVd33//vSysAoCuri769u0LQP6cP6T+d6lbty769u2L69evY926de95Jh9PXFwc5s6dq/Tf5s2blW4zePBgAFA4v1OnTuHKlSvo2LEjLC0tFbZr2bKlXFgFXn+cX716dRw9elT2frp//z6OHj0KS0tLjB8/Xq69s7Mzvv32W0gkEmzdulXt8/3qq6/kwioA9O3bF7q6unKvkYKCAvz6668wNTVFcHCw3GhxuXLlMHPmTIhEIrkaypV7HRne/pkEAFWqVHlnbR+6PZVt/HOGPkmF8zpFIhEqVqyI+vXrw9fXV2mQqFu3rsLH2Lm5ubhw4QIqV66MX3/9Vekx9PX1kZKSgvT0dFSpUgXnz58HAIVADLwetfniiy8U5qYpc/36dWRlZaF+/frFjrZ9iMJalX0cbmBggKZNm2Lnzp24cOGCwi9mFxcXhW0KQ1RxczU1TdU6inteypUrh6ZNm+LWrVtyywvntiYnJysdnStsf/36dYV1b/+xokmqnvOH1K+KKVOmIDIyEiEhIejZs+c7p59o06RJk4ockYyJiYGnp6fC8tq1a6N58+Y4evQo7ty5Azs7OwD/C7CDBg1Suj9lrzFdXV00adIEycnJsvdT4Rz5pk2bQl9fX2GbVq1aISwsTPbaVYey14ienh7MzMzkXiM3btxAWloaatSogfnz5yvdl5GRERITE2WPe/bsiaioKLRt2xbdunVDixYt0KhRI5V/Tn3o9lS2MbDSJ0md4KTsI9GMjAxIpVKkp6fLwm9RcnJyUKVKFdlFWKampiofR5msrCwAUDqCoymFtRZVk7m5uVy7NykbKS4cnVH1PreFIy0FBQVFtilcV9j2fet4n+clPT0dABAVFYWoqKgia3z27JlK+9MUVc/5Q+pXhbm5OcaMGYOff/4ZCxYswMyZM99rP0I2ZMgQxMbGIjw8HNOmTUNGRgZ27dqFmjVrolWrVkq3Keq5L3ztFb4WP+T99y5FfZKjo6Oj9DVy+/btd/6MK+Th4YGIiAgsXboUmzdvlt0i0MnJCZMmTVL4BEPT21PZxsBKZZ5IJFJYVvhD38nJCfHx8Srtp3AbZVfYA8CjR49U2k/haFVRN2jXhMJai6opNTVVrl1JHT8jI6PINoW/UD909O59npfCbdavXw8vLy+1jqfs9fSxfUj9qho9ejTCw8OxcuXKIkccgdf9UdQfMoV/nAlRly5dYGlpiQ0bNiAoKAibNm3CixcvMGDAgCKf46LeT4WvvcLnRdvvvzf33bFjR2zZskXl7dq2bYu2bdsiNzcX//zzDw4fPow1a9ZgwIABCvOlS2J7Krs4h5VIiQoVKsDJyQmJiYlIS0tTaZvCiwgKL+5609OnT+VulVWcWrVqoVKlSvjvv/9w7969d7ZXd3TzzVqV3RJHIpHIPlJ++8IITSm8Y0JRt9p5c506d1dQprjnpaCgQGkNjRo1AgCcPHnyg46tqjcvdtGEj1G/kZERpk6dColEghkzZhTZTiwWIzk5Wem6wtuzCZGuri769++Px48fY8+ePQgPD4eBgUGxFzYqe429evVK9n4qvHdt4X8TEhKU3jqq8LZjb368/z7v8+IU/pz5559/lNbwLkZGRmjevDmmT5+OWbNmQSqVKtypoiS3p7KHgZWoCCNHjsTLly8xYsQIpSOBT58+xZkzZ2SPO3fuDLFYjMjISLnlADBv3jyVP97T0dHBkCFDIJFIMHbsWOTm5sqtz8/Px8OHD2WPC2/OXVQoUKZnz57Q19fHmjVrFOYxLly4EA8ePIC7uzssLCxU3qc63NzcUKNGDVy6dAnr169XWH/hwgVs2LABurq6CheQqKtJkyZwdHREQkKCwkVGa9asUZi/Crx+LmvWrIl169YV+Uv0/PnzslHgD1X4HKryB4oqPlb9vXr1gouLCyIjI3Hu3DmlbRo1aoTk5GQcPHhQbnl4ePg774OrbQMGDICenh5++OEHXL9+Hd7e3sXeDP/EiRM4cOCA3LIVK1YgOTkZrVu3lk3zsbKyQtu2bXH//n0sXrxYrv1///2HtWvXwsDAAD179pQtf5/3eXF0dXUxfPhwPH78GBMnTsTz588V2qSlpcn9oX3s2DGl7QpHhA0NDYs95oduT2UbpwQQFcHf3x/nz5/Hb7/9BhcXF7Rt2xY2NjbIyspCUlIS4uPj0bp1a2zatAnA61HZxYsXY+DAgejSpQu6desGS0tLnDx5EleuXIGbm5vK0wu+//57nD17FkeOHEGDBg3QsWNHVK5cGQ8ePEBMTAz69u0ru5Dk66+/Rrly5fDrr78iIyNDNi9u6NChRX6cbmNjg7lz52L8+PFo3bo1unbtCnNzcyQkJCAuLg5WVlZYsGCBBnpROR0dHfz222/o3r07xowZg82bN8PV1RV6enq4fv06Dhw4gPz8fMybNw81atT4oGOJRCIsXboU3bp1w8CBA+XuwxodHY127drh8OHDctvo6elhw4YN8PHxQZ8+feDq6or69evD2NgY9+/fx4ULF5CYmIgTJ05o5Orm1q1bY+fOnfjuu+/g7e0NY2NjVKpUCUOHDn2v/X2s+kUiEWbPng0PDw+lwR8AxowZg8OHD6Nv377o2rUrTE1Nce7cOZw7dw4dOnRQCHhCYm5uDg8PD+zcuRMA8O233xbbvlOnTvD394eXlxfs7Oxw4cIFHD58GFWqVEFoaKhc24ULF6Jjx46YM2cOTpw4gUaNGsnuw5qbm4vFixfLfa1p48aNUaFCBURGRkJfXx/Vq1eHSCSCn5+fwr1YVRUYGIgrV65g/fr1OHjwIFq2bAkrKys8efIEt2/fxqlTpzB48GDZiPDUqVORlJSEZs2awcbGBoaGhrh8+TKOHDmCKlWqyO7SUZQP3Z7KNgZWomLMmzcP7u7uWLNmDWJjY5GRkYFKlSrB0tISgwYNgq+vr1x7b29vREREYO7cudi1axf09fXh5uaGQ4cOYdGiRSoHVn19fWzbtg3h4eHYvHkztm/fjlevXsHc3BzNmjVDp06dZG0dHBywZs0aLF68GBs2bJCNyL7r6u2BAweiZs2aWLp0Kfbu3Ytnz57BwsICQ4cOxcSJE0v04iHg9chbbGwsli1bhujoaKxevRr5+fkwNTWFl5cXhg0bhsaNG2vkWE2bNsVff/2FWbNm4ciRIzhy5AgaNmyIPXv24MiRIwqBFXg9fzkuLg4rVqzAvn37sHnzZkilUpibm6N27doYPXo0HB0dNVJf3759cf/+fWzbtg1hYWF4+fIlrK2t3zuwfsz6mzdvjs6dOxc5ktu8eXNs3boVISEhiIqKkntP7Nq1S9CBFXj93OzcuRNOTk5o2rRpsW09PDwwYMAAhIaGYv/+/dDT04O3tzemTZumcPs0W1tbHDt2TNb21KlTMDY2RrNmzTBmzBi0aNFCrn2lSpWwceNGBAcHIzIyEjk5OQBev7bfN7Dq6upi/fr1iIiIwMaNG3Ho0CHZRaTW1tYYN24cevXqJWs/YcIE7N27F2fPnpVNJ7K0tERAQABGjBghF7CV+dDtqWwTZWZm8rvQiIiIlFiwYAFmzZqF0NBQ2f1Z3xYcHIy5c+ciLCxM5S/vICL1cA4rERGREjk5OVi1ahVMTEw+eC41EX0YTgkgIiJ6w19//YWzZ8/i0KFDePjwIaZNm4aKFStquyyiMo2BlYiI6A1RUVHYvHkzzMzMMHbsWIwZM0bbJRGVeZzDSkRERESCxjmsRERERCRoDKxEREREJGgMrEREREQkaGUisCYmJmq7hFKHfaYe9pf62GfqYX+pj32mPvaZethf6nvfPisTgZWIiIiISi8GViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNF1tF0BEpAkpKfq4f//T/Rs8O9sOFSrow8IiT9ulEBF9dAysRPRJuH+/HCZO1NF2GSVGIjHC0qXlYGGh7UqIiD6+T3c4goiIiIg+CVoLrKtWrYKbmxusra1hbW2N9u3b48CBA7L1AQEBEIvFcv/atWunrXKJiIiISEu0NiXA0tISM2bMgL29PQoKCrB582b4+/vj2LFjqFu3LgCgVatWWLlypWwbfX19bZVLRERERFqitcDapUsXucc//vgj1qxZg9OnT8sCq4GBAczNzbVRHhEREREJhCDmsObn5yMiIgLPnj1D48aNZctPnjwJBwcHNGzYEGPGjMHjx4+1WCURERERaYMoMzNTqq2DX758Ge7u7njx4gWMjY2xatUqdOjQAQAQEREBIyMj2NraIikpCbNnz0ZBQQGOHTsGAwODIveZmJj4sconIgG5d88OU6YYabuMEjVnTi6sre9ouwwiIo1zdHQsdr1WA2teXh6Sk5ORlZWFqKgohIeHY8+ePXByclJom5KSgnr16mHt2rXw8vJS6ziJiYnv7AiSxz5TD/tLfZruszNnDD/x21pJsHSpLlxdX2i7lFKD70v1sc/Uw/5S3/v2mVbvw6qvr4+aNWsCAL788kv8+++/WL58OZYtW6bQ1sLCApaWlrh169bHLpOIiIiItEgQc1gLFRQUIC9P+be4pKWlISUlhRdhEREREZUxWhthnT59Otzd3WFlZYWcnBzs2LEDsbGx2LZtG3JychASEgIvLy+Ym5sjKSkJM2fOhKmpKTw8PLRVMhERERFpgdYCa2pqKoYOHYpHjx7BxMQEzs7O2LFjB9q2bYvc3FxcuXIFW7ZsQVZWFszNzdGiRQusW7cOFStW1FbJRERERKQFWgusK1asKHKdkZERIiMjP2I1RERERCRUgprDSkRERET0NgZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0rX1xABERqUck0sGZM4baLqNEWVkVwMIiT9tlEJHAMLASEZUST54AwcE62i6jRIWGAhYW2q6CiISGUwKIiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNAYWImIiIhI0LQWWFetWgU3NzdYW1vD2toa7du3x4EDB2TrpVIpgoODUbt2bVSrVg1dunTBf//9p61yiYiIiEhLtBZYLS0tMWPGDBw/fhzR0dFo2bIl/P39cenSJQDA4sWLERYWhrlz5+Lo0aMwNTVFt27d8PTpU22VTERERERaoLXA2qVLF7Rv3x41a9aEg4MDfvzxR1SoUAGnT5+GVCrFihUrMHbsWHh7e8PJyQkrVqxATk4OduzYoa2SiYiIiEgLBDGHNT8/HxEREXj27BkaN26Mu3fvIjU1FW3atJG1MTIygpubGxISErRYKRERERF9bLraPPjly5fh7u6OFy9ewNjYGBs2bICzs7MslJqamsq1NzU1RUpKSrH7TExMVGs5FY19ph72l/o02WfZ2XaQSIw0tj8hevnyFSSSfG2XUaKys3ORmHhHY/vj+1J97DP1sL/Up6zPHB0di91Gq4HV0dERMTExyMrKQlRUFAICArBnzx7ZepFIJNdeKpUqLFO2z7clJia+syNIHvtMPewv9Wm6z7KyDGFgoKOx/QmNRCKBnp4uDAy0+mO7xJmY6GrsdcH3pfrYZ+phf6nvfftMqz/59PX1UbNmTQDAl19+iX///RfLly/HxIkTAQCPHj1C9erVZe2fPHmiMOpKRERERJ82QcxhLVRQUIC8vDzY2trC3Nwc0dHRsnUvXrzAyZMn0aRJEy1WSEREREQfm9ZGWKdPnw53d3dYWVnJrv6PjY3Ftm3bIBKJEBAQgAULFsDR0REODg4IDQ2FsbExevTooa2SiYiIiEgLtBZYU1NTMXToUDx69AgmJiZwdnbGjh070LZtWwDAd999h9zcXAQGBiIzMxMNGzZEZGQkKlasqK2SiYiIiEgLtBZYV6xYUex6kUiEoKAgBAUFfaSKiIiIiEiIBDWHlYiIiIjobQysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaCoH1ri4ODx58qTI9WlpaYiLi9NIUUREREREhVQOrJ6enoiOji5y/fHjx+Hp6amRooiIiIiICqkcWKVSabHr8/LyUK4cZxgQERERkWbpFrcyOzsbWVlZssfp6em4d++eQrvMzExERETAwsJC8xUSERERUZlWbGBdvnw55s2bBwAQiUQICgpCUFCQ0rZSqRQ//vij5iskIiIiojKt2MDaqlUrGBoaQiqVYubMmfDx8UG9evXk2ohEIpQvXx5ffvklXF1dS7RYIiIiIip7ig2sTZs2RdOmTQEAEokEnp6ecHZ2/iiFEREREREB7wisb5o8eXJJ1kFEREREpFSRgXXz5s0AgF69ekEkEskev0vv3r01UxkREREREYoJrCNGjIBIJEL37t2hr6+PESNGvHNnIpGIgZWIiIiINKrIwHr+/HkAgL6+vtxjIiIiIqKPqcjAamNjU+xjIiIiIqKPgV9NRURERESCpvJdAgDg2LFjCA8Px507d5CRkaHwda0ikQjnzp3TZH1EREREVMapHFhXrFiBKVOm4LPPPoOrqyvq1KlTknUREREREQFQI7CGhYWhWbNmiIiIkF2I9SEWLlyI3bt348aNG9DX14erqyumTZsGJycnWZuAgACF22m5urri8OHDH3x8IiIiIiodVJ7DmpaWBh8fH42EVQCIjY3FoEGDcODAAURFRUFXVxddu3ZFRkaGXLtWrVrh2rVrsn/bt2/XyPGJiIiIqHRQeYTVxcUFSUlJGjtwZGSk3OOVK1fCxsYGp06dQqdOnWTLDQwMYG5urrHjEhEREVHpovII65w5c7Bp0yacOHGiRArJyclBQUEBxGKx3PKTJ0/CwcEBDRs2xJgxY/D48eMSOT4RERERCZMoMzNT+u5mgK+vL+7cuYObN2/C3t4e1tbW0NHRkd+ZSIRt27a9VyEDBgzAzZs3cezYMdl+IyIiYGRkBFtbWyQlJWH27NkoKCjAsWPHYGBgoHQ/iYmJ73V8Iird7t2zw5QpRtouo0RNmaKDOXPytV1GiZozJxfW1ne0XQYRfWSOjo7Frld5SsDVq1chEolQvXp1SCQS3LhxQ6GNSCRSv0IAP/zwA06dOoX9+/fLheDu3bvL/t/Z2RkuLi6oV68eDhw4AC8vL6X7UnbCiYmJ7+wIksc+Uw/7S32a7rOsLEMYGOi8u2EpJZFIoKenCwMDte5GWOqYmOhq7HXB96X62GfqYX+p7337TOWffBcvXlR756oICgpCZGQkdu/eDTs7u2LbWlhYwNLSErdu3SqRWoiIiIhIeLT6p/qkSZMQGRmJPXv2oFatWu9sn5aWhpSUFF6ERURERFSGqBxY7927p1I7a2trldpNnDgRW7duxYYNGyAWi5GamgoAMDY2RoUKFZCTk4OQkBB4eXnB3NwcSUlJmDlzJkxNTeHh4aFq2URERERUyqkcWL/44guV5qimp6ertL/Vq1cDALy9veWWT5o0CUFBQdDR0cGVK1ewZcsWZGVlwdzcHC1atMC6detQsWJFVcsmIiIiolJO5cC6bNkyhcCan5+Pu3fvYsuWLTAzM8PgwYNVPnBmZmax642MjBTu1UpEREREZY/KgdXf37/IdWPHjkWbNm2Qk5OjkaKIiIiIiAqp/MUBxalQoQL8/f2xfPlyTeyOiIiIiEhGI4EVAPT09JCSkqKp3RERERERAdBQYL148SJ+/fVXfP7555rYHRERERGRzAffJSArKwvZ2dmoUKECwsLCNFocEREREZHKgbVZs2YKgVUkEkEsFqNmzZro3r07xGKxpusjIiIiojJO5cC6YsWKkqyDiIiIiEgpjV10RURERERUEhhYiYiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQVAqsL168wNy5c3H06NGSroeIiIiISI5KgdXQ0BCLFi1CcnJySddDRERERCRH5SkB9erVw61bt0qyFiIiIiIiBSoH1p9++gnr16/HgQMHSrIeIiIiIiI5Kn/T1ZIlSyAWi9G7d29YWlrCzs4ORkZGcm1EIhG2bdum8SKJiIiIqOxSObBevXoVIpEI1atXBwAkJSUptBGJRJqrjIiIiIgIagTWixcvlmQdRERERERK8T6sRERERCRoagXW/Px8bNu2DaNGjYKfnx8uXboEAMjMzMTOnTvx8OHDEimSiIiIiMoulQNrVlYW3N3dMWzYMOzatQuHDh1CWloaAKBixYqYMmUKfvvttxIrlIiIiIjKJpUD64wZM3D16lVs374d586dg1Qqla3T0dGBp6cnDh06VCJFEhEREVHZpXJg3bt3L4YOHYp27dopvRuAvb097t27p9HiiIiIiIhUDqyZmZmoUaNGkeulUiny8vI0UhQRERERUSGVb2tlY2ODK1euFLk+Li4ODg4OGimKiDQrJUUf9+8L66Yg2dl2yMoy1Nj+cnN1NLYvIiISFpUDq6+vL3755Rd4enqiTp06AP73RQErV67Enj178PPPP5dMlUT0Qe7fL4eJE4UV6CQSIxgYaK6moCCN7YqIiARG5cA6btw4nDlzBl5eXnBwcIBIJMLkyZORnp6O1NRUdOnSBcOGDSvJWomI6BMnEungzBnNjLxrehRfU6ysCmBhwSl0ROpQObDq6elh27Zt2L59O/7880+IRCK8evUK9evXh4+PD3r27KnWV7MuXLgQu3fvxo0bN6Cvrw9XV1dMmzYNTk5OsjZSqRQhISEIDw9HZmYmGjZsiNDQUNkILxERfVqePAGCgzUz8q7pUXxNCQ0FLCy0XQVR6aJyYC3k6+sLX1/fDz5wbGwsBg0ahAYNGkAqleLnn39G165dkZCQgMqVKwMAFi9ejLCwMISFhcHR0RHz5s1Dt27dcPr0aVSsWPGDayAiIiIi4VM7sALApUuXZLewsra2hrOzs1qjqwAQGRkp93jlypWwsbHBqVOn0KlTJ0ilUqxYsQJjx46Ft7c3AGDFihVwdHTEjh07MHDgwPcpnYiIiIhKGbUuG46IiEDdunXRsmVL9OnTB3369EHLli1Rr149bN++/YMKycnJQUFBAcRiMQDg7t27SE1NRZs2bWRtjIyM4ObmhoSEhA86FhERERGVHiqPsG7cuBGjRo2Co6MjZsyYAQcHB0ilUty8eRPr16/HsGHDkJeXB39///cqZPLkyahXrx4aN24MAEhNTQUAmJqayrUzNTVFSkpKkftJTExUazkVjX2mHiH3V3a2HSQSI22XoUAikWhsXy9f6kAiydfY/oTo5ctXZeAcNfs8avI1pinZ2blITLyj7TKKJOSfZULE/lKfsj5zdHQsdhuVA+vChQvRsGFD7NmzB4aG8lddDhkyBJ07d8bChQvfK7D+8MMPOHXqFPbv3w8dHfkJ8m9PNZBKpcVOP1B2womJie/sCJLHPlOP0PsrK8tQcBefSCQSGBgYaGx/enqAgcF7zXIqFSQSCfT0dD/pcwQ0+zxq+jWmKSYmuoL9eSH0n2VCw/5S3/v2mcpTAu7fvw9fX1+FsAoAhoaG8PPzw4MHD9QuICgoCBEREYiKioKdnZ1subm5OQDg0aNHcu2fPHmiMOpKRERERJ8ulQNr7dq1i/0o/sGDB/j888/VOvikSZOwY8cOREVFoVatWnLrbG1tYW5ujujoaNmyFy9e4OTJk2jSpIlaxyEiIiKi0kvlz11mzpyJ/v37o379+ujWrZvcuoiICKxfvx7r169X+cATJ07E1q1bsWHDBojFYtmcVWNjY1SoUAEikQgBAQFYsGABHB0d4eDggNDQUBgbG6NHjx4qH4eIiIiISjeVA+vSpUtRtWpVDBo0CJMnT0aNGjUgEolw69YtPH78GPb29liyZAmWLFki20YkEmHbtm1K97d69WoAkN2yqtCkSZMQ9P/fsfjdd98hNzcXgYGBsi8OiIyM5D1YiYiIiMoQlQPr1atXIRKJUL16dQCQzVc1MDBA9erVIZFIcO3aNbltirs4KjMz853HFIlECAoKkgVYIiIiIip7VA6sFy9eLMk6iIiIiIiUUuuLA4iIiIiIPjYGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0lQNr/fr1sW/fviLX79+/H/Xr19dIUUREREREhVQOrElJSXj27FmR6589e4Z79+5ppCgiIiIiokJqTQko7purbty4wa9MJSIiIiKNK/abrjZt2oTNmzfLHoeGhiI8PFyhXWZmJq5cuYIOHTpovkIiIiIiKtOKDazPnj1Damqq7HFWVhYKCgrk2ohEIpQvXx79+/fH5MmTS6ZKIiIiIiqzig2sQ4YMwZAhQwAAX3zxBUJCQtC5c+ePUhgREREREfCOwPqmCxculGQdRERERERKqRxYCz19+hTJycnIyMiAVCpVWN+sWTONFEZEREREBKgRWDMyMjBp0iTs3LkT+fn5CuulUilEIhHS09M1WiARERERlW0qB9Zx48Zhz549GDJkCJo1awaxWFyCZRERERERvaZyYD18+DCGDRuGOXPmlGQ9RERERERyVP7iAH19fdjb25dkLUREREREClQOrN7e3jh06FBJ1kJEREREpEDlwDp69Gg8fPgQw4cPx+nTp/Hw4UM8fvxY4R8RERERkSapPIe1YcOGEIlEOHfuHLZt21ZkO94lgIiIiIg0SeXA+v3330MkEpVkLUREREREClQOrEFBQSVZBxERERGRUirPYX1Tfn4+0tPT8erVK03XQ0REREQkR63A+u+//6Jr166wtLSEg4MD4uLiAABpaWno2bMnjh8/XiJFEhEREVHZpXJg/fvvv9G5c2fcvn0bvXr1glQqla2rWrUqcnJy8Mcff5RIkURERERUdqkcWGfNmgV7e3skJCTgp59+UljfokULnDlzRq2Dx8XFoVevXqhTpw7EYjE2btwotz4gIABisVjuX7t27dQ6BhERERGVbioH1n///Rd9+/aFoaGh0rsFWFlZITU1Va2DP3v2DE5OTggJCYGRkZHSNq1atcK1a9dk/7Zv367WMYiIiIiodFP5LgHlypVDuXJF59vU1NQiQ2dR3N3d4e7uDgAYMWKE0jYGBgYwNzdXa79ERERE9OlQeYTVxcUF+/fvV7ouLy8P27dvR+PGjTVWWKGTJ0/CwcEBDRs2xJgxY/htWkRERERljMojrOPHj0ePHj0watQo+Pr6AgAePnyIw4cPIzQ0FLdv30ZYWJhGi2vXrh08PT1ha2uLpKQkzJ49G15eXjh27BgMDAyUbpOYmKjWcioa+0w9Qu6v7Gw7SCTqfQLyMUgkEo3t6+VLHUgk+RrbnxC9fPmqDJyjZp9HTb7GNCU7OxeJiXe0XUaRhPyzTIjYX+pT1meOjo7FbqNyYG3dujVWrlyJwMBAbNq0CcDri6KkUikqVaqE1atXo1GjRmqWXLzu3bvL/t/Z2RkuLi6oV68eDhw4AC8vL6XbKDvhxMTEd3YEyWOfqUfo/ZWVZQgDAx1tlyFHIpEU+Yfn+9DTAwwMVP6RVupIJBLo6el+0ucIaPZ51PRrTFNMTHQF+/NC6D/LhIb9pb737TO1fir06NEDnTt3RnR0NG7evImCggLUqFEDbdu2RYUKFdQ+uLosLCxgaWmJW7dulfixiIiIiEgY1P4ztnz58ujSpUtJ1PJOaWlpSElJ4UVYRERERGWIyhdd7du3D4GBgUWuDwwMLPKirKLk5OTgwoULuHDhAgoKCpCcnIwLFy7g3r17yMnJwdSpU/H333/j7t27iImJQa9evWBqagoPDw+1jkNEREREpZfKgXXp0qV4/vx5ketfvHiBxYsXq3Xws2fPomXLlmjZsiVyc3MRHByMli1b4ueff4aOjg6uXLmCPn36wNXVFQEBAXBwcMDBgwdRsWJFtY5DRERERKWXylMCrly5Ah8fnyLX169fH3v27FHr4C1atEBmZmaR6yMjI9XaHxERERF9elQeYX316hVyc3OLXJ+bmyvI24cQERERUemmcmB1cnJCVFQUCgoKFNYVFBQgKioKtWvX1mhxREREREQqB9bhw4fjn3/+Qe/evXHu3DlIJBJIJBKcO3cOffr0wT///INhw4aVZK1EREREVAapPIe1e/fuuH37NoKDg3Ho0CEAgEgkglQqhUgkwqRJk+Dn51dihRIRERFR2aTWfVgnTpyIHj16YPfu3bhz5w6kUilq1KgBT09P2NnZlVCJRERERFSWqRRYc3Nz0bNnT/j5+aFv374YPXp0SddFRERERARAxTmsRkZGOH/+PPLz80u6HiIiIiIiOSpfdNW8eXPEx8eXZC1ERERERApUDqxz587Fv//+ix9//BF37txRensrIiIiIiJNU/miq0aNGkEqlSIsLAxhYWEoV64c9PT05NqIRCI8ePBA40USERERUdmlcmDt1q0bRCJRSdZCRERERKRA5cC6YsWKkqyDiIiIiEgpleewEhERERFpg1qBNSkpCWPGjIGLiwusra0RGxsLAEhLS8OECRNw7ty5kqiRiIiIiMowlacEXLt2DR07dkRBQQFcXV2RlJQkuy9r1apVcfr0aUgkEixbtqzEiiUiIiKiskflwDpt2jRUrFgRhw8fho6ODhwcHOTWu7u7488//9R0fURERERUxqk8JSA+Ph6DBw+GmZmZ0rsFWFtbIyUlRaPFERERERGpHFhfvXoFY2PjItdnZGRAR0dHI0URERERERVSObA6OTkhJiZG6TqpVIrdu3fDxcVFU3UREREREQFQI7AGBARg165dmDdvHtLT0wEABQUFuH79Or799lucPXsWo0ePLrFCiYiIiKhsUvmiq+7du+PevXuYM2cOQkJCZMsAQEdHB7Nnz0b79u1LpkoiIiIiKrNUDqwAMHbsWPTo0QNRUVG4desWCgoKUKNGDXh5ecHW1rakaiQiIiKiMuydgVUikWDfvn24c+cOqlSpgg4dOmDEiBEfozYiIiIiouIDa2pqKjp37ozbt29DKpUCAIyNjbF161Y0a9bsoxRIRERERGVbsRddzZ49G3fu3MGIESOwdetWBAcHw8DAAN9///3Hqo+IiIiIyrhiR1iPHj2K3r17Y/bs2bJlZmZmGDx4MO7fvw8rK6sSL5CIiIiIyrZiR1hTU1PRpEkTuWVNmzaFVCpFcnJyiRZGRERERAS8I7Dm5+fD0NBQblnh4xcvXpRcVURERERE/++ddwm4c+cO/vnnH9nj7OxsAEBiYiIqVKig0L5hw4YqHzwuLg5Lly7F+fPnkZKSgrCwMPj7+8vWS6VShISEIDw8HJmZmWjYsCFCQ0NRp04dlY9BRERERKXbOwNrcHAwgoODFZa/feGVVCqFSCSSfQuWKp49ewYnJyf07t0bw4cPV1i/ePFihIWFISwsDI6Ojpg3bx66deuG06dPo2LFiiofh4iIiIhKr2IDa1hYWIke3N3dHe7u7gCgcG9XqVSKFStWYOzYsfD29gYArFixAo6OjtixYwcGDhxYorURERERkTAUG1j79OnzsepQcPfuXaSmpqJNmzayZUZGRnBzc0NCQgIDKxEREVEZodZXs35MqampAABTU1O55aampkhJSSlyu8TERLWWU9HYZ+oRcn9lZ9tBIjHSdhkKJBKJxvb18qUOJJJ8je1PiF6+fFUGzlGzz6MmX2Oakp2di8TEO9ouo0hC/lkmROwv9SnrM0dHx2K3EWxgLSQSieQeF86VLYqyE05MTHxnR5A89pl6hN5fWVmGMDDQ0XYZciQSCQwMDDS2Pz09wMBA8D/S3ptEIoGenu4nfY6AZp9HTb/GNMXERFewPy+E/rNMaNhf6nvfPiv2tlbaZG5uDgB49OiR3PInT54ojLoSERER0adLsIHV1tYW5ubmiI6Oli178eIFTp48qfBlBkRERET06dLqZ0s5OTm4desWAKCgoADJycm4cOECKleuDGtrawQEBGDBggVwdHSEg4MDQkNDYWxsjB49emizbCIiIiL6iLQaWM+ePQtPT0/Z48J7vvbu3RsrVqzAd999h9zcXAQGBsq+OCAyMpL3YCWNSknRx/37H/ZhQ3a2HbKyDN/dUEtyc4U1f5WIiEgdWg2sLVq0QGZmZpHrRSIRgoKCEBQU9PGKojLn/v1ymDjxwwKdRGIkuIua3sS3EBERlWaCncNKRERERAQwsBIRERGRwH3aN/QjIiISGJFIB2fOCHPOuybm41tZFcDCIk9DFRG9xsBKRET0ET15AgQHC3POuybm44eGAhYWGiqI6P9xSgARERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQmaoANrcHAwxGKx3L9atWppuywiIiIi+oh0tV3Auzg6OmLPnj2yxzo6OlqshoiIiIg+NsEHVl1dXZibm2u7DCIiIiLSEkFPCQCAO3fuoE6dOvjiiy/w7bff4s6dO9ouiYiIiIg+IkGPsLq6umL58uVwdHTEkydPMH/+fLi7u+PUqVOoUqWK0m0SExPVWk5FKyt9lp1tB4nE6IP3I5FINFBNyXj5UgcSSb62y1CgyT4T6jlq0suXr8rAOWr2eRTi+1Lor9UP7bPs7FwkJt7RTDGlQFn5XalJyvrM0dGx2G0EHVjbt28v99jV1RUuLi7YtGkTRo0apXQbZSecmJj4zo4geWWpz7KyDGFg8GFzoyUSCQwMDDRUkebp6QEGBsJ6u2u6z4R4jpokkUigp6f7SZ8joNnnUajvSyG/VjXRZyYmumXm90dZ+l2pKe/bZ4KfEvCmChUqoHbt2rh165a2SyEiIiKij6RUBdYXL14gMTGRF2ERERERlSHC/Ezi/02dOhUdO3ZE9erVZXNYnz9/jt69e2u7NCIiIiL6SAQdWB88eIDBgwcjLS0Nn332GVxdXXHo0CHY2NhouzQiIiIi+kgEHVjXrl2r7RKIiIiISMtK1RxWIiIiIip7GFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNB0tV0AERERfTpEIh2cOWOo7TJKlJVVASws8rRdRpnCwErFSknRx/37n/ZAfG6ujrZLICL6ZDx5AgQHf9o/V0NDAQsLbVdRtjCwUrHu3y+HiRM/7R88QUHaroCIiIiK82kPnRERERFRqcfASkRERESCxsBKRERERILGwEpEREREgsbASkRERESCxsBKRERERILGwEpEREREgsbASkRERESCxsBKRERERILGwEpEREREgsbASkRERESCxsBKRERERILGwEpEREREglYqAuvq1avxxRdfwNzcHF9//TXi4+O1XRIRERERfSS62i7gXSIjIzF58mQsWLAATZs2xerVq+Hr64tTp07B2tpaq7WlpOjj/v1SkfnVlp1th6wsQ+Tm6mi7FCIiIkERiXRw5oyh7Hflp8jKqgAWFnnaLkNG8IE1LCwMffr0Qf/+/QEA8+fPx5EjR7B27VpMmzZNq7Xdv18OEyd+moFOIjGCgYEOgoK0XQkREZGwPHkCBAfryH5XfopCQwELC21X8T+izMxMqbaLKEpeXh4sLCywZs0adO3aVbZ84sSJuHLlCvbt26e94oiIiIjooxD059lpaWnIz8+Hqamp3HJTU1M8evRIS1URERER0cck6MBaSCQSyT2WSqUKy4iIiIjo0yTowFq1alXo6OgojKY+efJEYdSViIiIiD5Ngg6s+vr6cHFxQXR0tNzy6OhoNGnSREtVEREREdHHJPi7BIwcORLDhg1Dw4YN0aRJE6xduxYPHz7EwIEDtV0aEREREX0Egh5hBQAfHx8EBwdj/vz5aNGiBU6dOoVt27bBxsbmvfcplUrRvXt3iMVi7Nq1S4PVfnrGjBkDFxcXVKtWDfb29ujduzeuXbum7bIEKSMjA4GBgWjUqBGqVasGZ2dnjB8/Hunp6douTdB+//13eHh4wMbGBmKxGHfv3tV2SYLDL09RXVxcHHr16oU6depALBZj48aN2i5J0BYuXIjWrVvD2toa9vb28PPzw5UrV7RdlqCtWrUKbm5usLa2hrW1Ndq3b48DBw5ou6xSY8GCBRCLxQgMDFRrO8EHVgAYPHgwLl68iEePHuH48eNo1qzZB+1v2bJl0NH5NO+bpmlffvklli9fjoSEBEREREAqlaJr1654+fKltksTnJSUFKSkpGDGjBmIj4/HypUrER8fj0GDBmm7NEF7/vw52rRpg8mTJ2u7FEEq/PKUCRMm4MSJE2jcuDF8fX1x7949bZcmSM+ePYOTkxNCQkJgZGSk7XIELzY2FoMGDcKBAwcQFRUFXV1ddO3aFRkZGdouTbAsLS0xY8YMHD9+HNHR0WjZsiX8/f1x6dIlbZcmeKdPn0Z4eDicnZ3V3lbQ92EtCWfPnkXfvn1x7NgxODo6Ijw8HN7e3touq9S4dOkSmjdvjtOnT8PR0VHb5QjewYMH4efnh7t378LExETb5Qja2bNn0bp1a5w/fx62trbaLkcw2rZtC2dnZyxZskS2rEGDBvD29tb6l6cInZWVFebNmwd/f39tl1Jq5OTkwMbGBhs3bkSnTp20XU6pYWdnh2nTpnG6YjGysrLw9ddfY/HixZg3bx6cnJwwf/58lbcvFSOsmvL06VMMGjQIixYt4l0G3sOzZ8+wceNGVK9e/YOmZJQlT58+hYGBAcqXL6/tUqgUysvLw7lz59CmTRu55W3atEFCQoKWqqJPWU5ODgoKCiAWi7VdSqmQn5+PiIgIPHv2DI0bN9Z2OYI2duxYeHt74+uvv36v7QV/0ZUmjR8/Hm3btoW7u7u2SylVVq9ejWnTpuHZs2dwdHREVFQUDAwMtF2W4GVmZmLOnDno168fdHXL1FuNNIRfnkIf2+TJk1GvXj2Gr3e4fPky3N3d8eLFCxgbG2PDhg3v9TF3WREeHo5bt25h5cqV772PUj/COnv2bIjF4mL/xcTEYMuWLbh06RJmzZql7ZK1TtU+K+Tr64sTJ05g7969sLe3R//+/fH8+XMtnsHHpW5/Aa9Ho3v37g0LCwvMnDlTS5Vrz/v0GRWNX55CH8MPP/yAU6dO4Y8//uB1Hu/g6OiImJgYHD58GIMGDUJAQAAvVitCYmIiZs6ciVWrVkFfX/+991Pq57CmpaUhLS2t2DbVq1fHhAkTsGXLFpQr97+Mnp+fj3LlyqFx48bYv39/SZcqGKr2mbKPsfPy8mBnZ4eFCxeiV69eJVWioKjbXzk5OfD19QUAbN++HRUqVCjxGoXmfV5jnMOqKC8vDxYWFlizZg26du0qWz5x4kRcuXIF+/bt015xpQDnsKouKCgIkZGR2L17N2rVqqXtckodb29vWFtbY9myZdouRXA2btyIkSNHyv0RlJ+fD5FIhHLlyuHBgwcqfWpb6j+nrFq1KqpWrfrOdj/++CNGjx4tt8zNzQ2zZs1Cly5dSqo8QVK1z5SRSqWQSqXIy8vTcFXCpU5/PX36FL6+vpBKpdixY0eZDKvAh73G6H/e/PKUNwNrdHQ0vLy8tFcYfVImTZqEyMhI7Nmzh2H1PRUUFJSp34vq6NKlC7788ku5ZSNHjoS9vT3Gjx+v8qhrqQ+sqrK0tISlpaXC8urVq8POzu7jF1QK3Lp1C1FRUWjVqhWqVq2KBw8eYNGiRdDX10eHDh20XZ7gPH36FD4+Pnj69Ck2btyI58+fy6ZOVK5c+YM+CvmUpaamIjU1FTdu3AAAXLt2DVlZWbC2tkblypW1XJ328ctT1JOTk4Nbt24BeB0ikpOTceHCBVSuXBnW1tZark54Jk6ciK1bt2LDhg0Qi8VITU0FABgbG5fZP7jfZfr06XB3d4eVlRVycnKwY8cOxMbGYtu2bdouTZAKp4G9qXz58qhcuTKcnJxU3k+ZCaykPn19fcTGxmLZsmXIysqCmZkZ3NzccOjQIZibm2u7PME5d+4cTp8+DQBo2LCh3Lrdu3ejRYsW2ihL8NauXYu5c+fKHvfs2RMAEBYWxo9y8frLU9LT0zF//nykpqaiTp06H/zlKZ+ys2fPwtPTU/Y4ODgYwcHB6N27N1asWKHFyoRp9erVAKBwe8dJkyYhKChIGyUJXmpqKoYOHYpHjx7BxMQEzs7O2LFjB9q2bavt0j5ppX4OKxERERF92kr9XQKIiIiI6NPGwEpEREREgsbASkRERESCxsBKRERERILGwEpEREREgsbASkRERESCxsBKVEoFBwcr3IwZAJYvX44vv/wSVatWRb169T5+YURERBrGwEokABs3bpR9G4hYLIa5uTlq164NHx8f/Prrr3j69KlK+4mNjcUPP/yA+vXrY+nSpQgODi7hyoWvXr16cn1rZmYGFxcXBAYGyr7VR105OTkIDg5GTEyMhqsVpqNHj6Jv376oXbs2TE1NYWNjg7Zt22LOnDl48OCBtst7bzExMbLXxYYNG5S2GTBggOw9SUTaw2+6IhKQyZMno0aNGnj58iUePXqE2NhYBAUFISwsDJs3b0bdunVlbQMDAzFu3Di57QsD1C+//KJ09LWscnZ2xpgxYwAAEokEFy9eRHh4OE6cOIH4+Hjo6Oiotb9nz57Jvp3rU/4GM6lUigkTJmDt2rWoXbs2+vXrB2tra+Tm5uLcuXNYuXIl1q1bJ/ta3dLK0NAQ27dvR9++feWWZ2dnY//+/TA0NIRUyu/YIdImBlYiAWnbti0aNWokezx+/HgcP34cvXr1Qu/evfH333/DyMgIAKCrqwtdXfm38JMnTwBAo2H1+fPnKF++vMb2pw3VqlWDn5+f3LJKlSohNDQUFy9ehIuLi3YKE7iwsDCsXbsWQ4cORUhICMqVk/9QLiQkBEuWLHnnfnJzc2WvWyFyd3fHnj17kJKSAgsLC9nyqKgoFBQUoE2bNoiOjtZihUTEKQFEAvf1118jMDAQ9+7dw7Zt22TL357DKhaLsWbNGtn/i8ViuSkB0dHR8PDwQPXq1WFpaQkPDw8kJCTIHatwn1evXsXw4cNRo0YNNG3a9L32cfPmTYwbNw41atSAlZUV+vfvj/T0dIXzi46OhqenJ6ytrVG9enV8/fXXWL9+vVybs2fPws/PDzY2NqhWrRratGmD/fv3q9+Zbyj8iFdPT09u+cOHD/Hdd9+hdu3aMDMzQ4MGDbB48WLZCNvdu3fx+eefAwDmzp0r6+uAgABcunQJYrEYO3fulO3v9u3bEIvFcqPjADBu3DjUqlXrvc4zOzsbU6dORb169WBmZoa6deti+vTpkEgkcu3EYjHGjRuHQ4cOoUWLFjA3N0eDBg2wY8eOd/ZPbm4uFixYgM8//xzBwcEKYRUATExMMHXqVLll9erVQ/fu3XHixAm0a9cO5ubm+OWXXwAAGRkZGD9+PD7//HOYmZmhcePGWLZsmdzo5d27dyEWi7Fx40aF49WrVw8BAQGyx4VTaU6cOIHAwEDUrFkTVlZW6NevHx4+fPjOcyzUoUMHVKhQQaFftm/fjnbt2qFy5cpKt1Pl/ZCUlIQJEyagUaNGsLCwgI2NDfz8/PDff//JtSucnrBjxw4sW7YM9erVg7m5Odq3b4/z58/LtX306BFGjx4NZ2dnmJmZoXbt2vDz88Ply5dVPmei0oaBlagUKBwdPHr0aJFtVq5ciZYtW8r+f+XKlfD09AQA7NixA927d4eOjg6mTJmCKVOmID09HV5eXjhz5ozCvgYOHIiMjAxMmTIFw4cPf699DBo0CA8ePMCUKVPQr18/7NmzB99//71cmy1btsDHxwcPHz7E6NGjMWPGDDRs2BAHDhyQtYmNjUXHjh3x6NEjBAYGYsaMGdDX10fv3r0RFRWlUv+9fPkSaWlpSEtLw4MHD3Dw4EEsXrwYLi4ucHJykrV7/Pgx2rVrhwMHDqB///6YO3cuXF1dMW3aNAQFBQEAPvvsM8yfPx8A4OHhIevrgQMHwtnZGWKxGHFxcbJ9xsXFoVy5ckhOTsbdu3dly+Pj4/HVV1+pfZ65ubnw8PDAH3/8AR8fH8ybNw8dOnTAsmXLMHDgQIVzP336NEaOHInOnTtj1qxZKF++PIYOHYpr164V22cJCQnIyMhAjx491J4ycevWLfTr1w9ubm6YO3cuGjVqBIlEAk9PT4SHh8PLywtz5syBra0tpk6dih9++EGt/b9t8uTJOHfuHL7//nsMGDAAf/31F3x8fJCXl6fS9oaGhvD09MT27dtlyx4+fIiYmBj07NlT6Taqvh/Onj2LuLg4eHp6Ijg4GAEBATh79iw6d+6sdA71smXLsHnzZgwdOhSTJ0/GjRs34O/vj5cvX8ra9O/fH7t27ULv3r0RGhqKYcOGoaCgoNRPzSAqDqcEEJUCVlZWMDExwe3bt4ts4+fnh1OnTuHEiRNyH38/e/YMEydOhJ+fH1asWCFbPnDgQDRt2hQzZ85UCH4ODg74448/PmgftWrVwm+//SZ7LJVKsWrVKixYsACVKlVCdnY2vv/+ezg7O+PAgQMwNjaWa1v433HjxqFx48bYtWuXbJRvyJAh6NChA3766Sd4eXm9s/9OnDgBe3t7uWUNGzbEli1bIBKJZMtmz54NiUSCuLg4mJmZyc6xWrVqWLZsGQICAmBrawsvLy8EBgbC2dlZYapB06ZNER8fL3t88uRJtGnTBgkJCTh58iRsbW2RlpaG69ev49tvv1X7PJcvX47ExEQcO3ZMNtILAHXq1MHEiRMRHx8PNzc32fKrV68iLi5O1rZr166oW7cuNmzYgFmzZhXZZ1evXpXt901SqVRhpNzExERupPr27dvYtGkTOnfuLFv222+/4dKlS1iyZAn69esHABg8eDC++eYb/Prrrxg8eLDCc6SOPXv2wMDAAABQu3ZtjB49Gps2bcKAAQNU2r5nz57YuHEjrl69itq1a2P79u2oUKECOnbsKPcHFKDe+6F9+/bw9vaW297Pzw9fffUV/vjjD0ycOFFuXXZ2NuLj42FoaAgAcHR0RN++fXH06FF06NABWVlZOHnyJGbNmoXRo0fLtnt7PjvRp4YjrESlRIUKFZCTk6P2dtHR0cjMzETPnj1lo4xpaWnIzc1Fq1atcPLkSbnRG+D16Kim99GsWTPk5+cjOTlZts/s7GxMmDBBLqwCkIXIixcvIjExET179kRGRobsuBkZGWjXrh3u3LmDpKSkd/bBl19+iT///BN//vkntm7dip9++gk3b95E3759kZubC+B1ENu1axc6dOgAHR0dufNs27YtCgoK5EZOi+Lm5ob//vsPGRkZAF6PpLZo0QKNGzeWBdm4uDhIpVLZCKs657lz5040adIEn332mVyNrVq1AvA6nL+pRYsWcsHWzMwMjo6OuHPnTrHnUXhniooVK8otT09Ph729vdy/U6dOybWxsrKSC6sAcODAAVStWhX+/v6yZSKRCGPGjIFUKsXBgweLrac4AwcOlIVVAOjduzcqVaqk1j5btGgBCwsL2Sjr9u3b4eHhIQuOb1Ln/fDm/O/nz58jPT0dlSpVgr29Pc6dO6ewb39/f7ljNm/eHABkz5ehoSH09PQQGxsre40RlQUcYSUqJXJycvDZZ5+pvd3NmzcBAN26dSuyTVZWlty+7ezsPngf1tbWcusL59sW/pItHC1+8yP5omofPXq03GjSm548eQIbG5si9wEAVapUkQU64PWcRUdHR3zzzTdYv349hg0bhidPniAzMxMbNmwo8hZHhRe1FcfNzQ1SqRTx8fFo2LAhbt++DTc3N7x69QpbtmwB8DrEmpiYyOa1qnOeN2/exKVLl4ocjXy7xrefB+D1c/GusFMYVN++pZqJiQn+/PNPAK+nMYSGhipsa2trq7AsKSkJ9vb2CtMLCsO0Kn94FOXtvtDV1YWtrS3u3bun8j7KlSsHHx8fbN++Hb6+vrhw4QJmzpyptK0674cXL17g559/xrZt2xTm1VatWlVhu3e9bwwMDDBt2jRMmzYNjo6OcHV1Rfv27dGzZ0+lzzXRp4KBlagUuH//PrKzs1GzZk21ty0oKADw+qNkS0tLpW1MTEzkHr99Rff77KOoeY9vftwPQO4j+aJqnz59epFX8js4OBS5fXEK5/vGx8fL5gACQI8ePRRub1RIlf53cXGBsbEx4uPjIZFIUL58ebi4uODly5eYNWsWHj9+LJu/WvjRvzrnWVBQgJYtW2L8+PFK2739/LzreShK7dq1AQBXrlyBh4eHbLmenp4s/KelpSnd9kPuCKDK60GVbd7nNlS+vr4ICwtDYGAgqlWrVuQty9R5P0yePBnr16/H0KFD0bRpU5iYmKBcuXIICgpSej6qPF+jRo2Ch4cH9u3bh2PHjmH+/PlYuHAhNm3ahK+//lqtcyYqLRhYiUqBrVu3AgDatGmj9rY1atQA8PpioTdHGT/2Pt5WGP6uXLmicLX828etUKGCxo5b6NWrVwBez0cEXp+biYkJXr169c5jFReqdHV14erqivj4eOTl5aFRo0bQ09NDw4YNYWBggP379+Py5cvw8fGRbaPOedaoUQM5OTka74+3NWnSBJUrV8aOHTswYcIEtS+8epuNjQ3Onz+P/Px8uX1dv35dth6A7Ir8rKwsue0lEkmRV/7fuHEDrVu3lj1+9eoVkpKS0KxZM7VqdHFxQa1atRATE4MRI0YUec7qvB8iIyPRq1cvhISEyC3PzMxElSpV1KrvTXZ2dhgxYgRGjBiB5ORktGzZEosWLWJgpU8W57ASCdzx48cxf/582NraFnnFcnHatm0ru+fo27c9AlT7mFsT+3hb69atYWJigoULF+L58+dy6wpHk1xcXGBvb4+lS5cqBJj3PW6hwvmNhR/L6+jowMvLC3v27FE6tzArK0thXmJmZqbSfbu5ueHChQs4fPiw7AIoAwMDNGjQAEuWLEF+fr7chVHqnKePjw/+/fdf7Nu3T6Fdbm7ue81zVsbIyAgTJkzA9evXixwNVGcUs0OHDnjy5Ak2b94st/3SpUshEong7u4O4PVUhM8++0zhW8TWrl2L/Px8pftet26d3Oty8+bNyMrKQvv27VWur9CcOXMwadIkhTnYb1Ln/aCjo6PQTzt27EBKSoratQGv58EWzrsuVL16dZiamhb5eiT6FHCElUhAjhw5glu3buHVq1d4/PgxTpw4gejoaFhbW2Pz5s1KLwB5l4oVK2Lx4sUYNGgQmjdvDl9fX5ibm+P+/fuIiYmBsbHxO+/LqYl9vM3ExATBwcEYNWoUWrduDV9fX1SpUgX//fcfUlJSsGHDBpQrVw7Lli1D9+7d0bRpU/j7+8PGxgYPHz7E6dOnce/ePYULfpR5+PChbJQ6Ly8Ply9fxu+//46qVati6NChsnbTp09HXFwcOnbsiG+++QZOTk54+vQprly5gt27d+Pff/+Fubk5KlSoAEdHR0RGRsLBwQFVqlSBra0tXF1dAQBfffUV8vPzZfNXCzVr1gyhoaEwMjKS++hfnfMcPXo0Dh48iG+++QY9e/ZEw4YNIZFIcOPGDezcuRPbt2+X+/KJDzFy5EjcvHkTv/32G06cOAEvLy9YW1vj2bNnuHbtGiIiImBsbKx0Lubb+vXrh/Xr12Ps2LG4ePEiHBwccOjQIRw8eBDDhw+Xm4c6YMAAhIaGYsSIEWjUqBHOnj2L48ePF3scT09PdO/eHUlJSfjtt99Qu3Zt9OnTR+1zbt++/TuDrjrvh06dOmHLli2oWLEinJyccPHiRURGRirME1fVjRs34OXlha5du6J27dowMDDAwYMHce3atWLv+kBU2jGwEglI4ceG+vr6qFy5MpycnBAcHAx/f3+Fq7XV0bVrV1hYWGDhwoVYvnw5cnNzYW5uDldXV9kthj7GPt7m7+8PU1NTLFq0CAsXLoSOjg7s7e0xePBgWZuvvvoKR44cwbx58/D7778jOzsbpqamqFu3ruzeqO9y+fJlDBs2DMDrcFi1alV4eHhgypQpcnMQP/vsMxw5cgTz58/H3r178fvvv6NSpUpwcHDA5MmT5W4gHxYWhqCgIEydOhUSiQS9e/eWBdZGjRpBX18fAGTLCs+lcFnhenXP08jICFFRUVi8eDEiIyNlodHOzg4BAQFwdHRUqU9UIRKJsGjRInTp0gVr165FeHg40tLSUL58eTg4OGD48OEYMGBAkfM432RoaIioqCjMmjULO3fuREZGBmxtbTFr1iyMGjVKru3EiRORnp6OyMhI/Pnnn2jevDl27dolu6/w20JCQhAVFYW5c+dCIpGgQ4cOmD9/vtydAzRN1fdDSEgI9PT0sHPnTmzYsAEuLi6IiIjAjz/++F7HrV69Onx9fXHixAns2LEDIpFINjr/zTffaOr0iARHlJmZyS9IJiKiUmfjxo0YOXIkDh06pLFRZSISJs5hJSIiIiJBY2AlIiIiIkFjYCUiIiIiQeMcViIiIiISNI6wEhEREZGgMbASERERkaAxsBIRERGRoDGwEhEREZGgMbASERERkaAxsBIRERGRoP0fLt5M2pxOg+QAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"test_conclusion = pd.DataFrame({'Difference Between Group Means':differences})\n",
"\n",
"print('Observed Difference:', observed_difference)\n",
"\n",
"unit = ''\n",
"\n",
"fig, ax = plt.subplots(figsize=(10,5))\n",
"\n",
"ax.hist(test_conclusion, density=True, color='blue', alpha=0.8, ec='white')\n",
"\n",
"y_vals = ax.get_yticks()\n",
"\n",
"y_label = 'Percent per ' + (unit if unit else 'unit')\n",
"\n",
"x_label = 'Diference Between Group Means'\n",
"\n",
"ax.set_yticklabels(['{:g}'.format(x * 100) for x in y_vals])\n",
"\n",
"plt.ylabel(y_label)\n",
"\n",
"plt.xlabel(x_label)\n",
"\n",
"plt.title('Prediction Under the Null Hypothesis');\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice how the distribution is centered around 0. This makes sense, because under the null hypothesis the two groups should have roughly the same average. Therefore the difference between the group averages should be around 0.\n",
"\n",
"The observed difference in the original sample is about $-9.27$ ounces, which doesn't make an appearance on the horizontal scale of the histogram. The observed value of the statistic and the predicted behavior of the statistic under the null hypothesis are inconsistent. \n",
"\n",
"The conclusion of the test is that the data favour the alternative over the null. The average birth weight of babies born to mothers who smoke is less than the average birth weight of babies born to non-smokers.\n",
"\n",
"If you want to compute an empirical P-value, remember that low values of the statistic favor the alternative hypothesis. "
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"empirical_P = np.count_nonzero(differences <= observed_difference) / repetitions\n",
"empirical_P"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The empirical P-value is 0, meaning that none of the 5,000 permuted samples resulted in a difference of -9.27 or lower. This is only an approximation. The exact chance of getting a difference in that range is not 0 but it is vanishingly small."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Another Permutation Test\n",
"We can use the same method to compare other attributes of the smokers and the non-smokers, such as their ages. Histograms of the ages of the two groups show that in the sample, the mothers who smoked tended to be younger."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA6UAAAFDCAYAAADCo6jwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABD30lEQVR4nO3dfVxUdf7//+eZ4UoRxTUFBBRB0iDUVjdLSzNTy8tMzdRat01NLD+bZV61Wl6FJrllXmRZfTW1vEjLLLMyMy9iTe3CMo1C8iJFo4AkAWHm94c/Z5sQHGBmDsLjfrtxuznnnPd5v86bw8En58rIysqyCwAAAAAAE1jMLgAAAAAAUH0RSgEAAAAApiGUAgAAAABMQygFAAAAAJiGUAoAAAAAMA2hFAAAAABgGkIpAAAAAMA0hFIAAAAAgGlMC6VFRUWaMWOGWrRooZCQELVo0UIzZsxQYWGhWSU5SU1NNbuEaouxNw9jbx7G3jyMvXkYe/Mw9gAqEx+zOn7mmWe0ZMkSLVq0SHFxcfrmm2+UmJgoPz8/jRs3zqyyAAAAAABeZFoo3b17t2699VbddtttkqTGjRvrtttu0969e80qCQAAAADgZaZdvnvddddpx44d+u677yRJBw8e1Pbt29WlSxezSgIAAAAAeJmRlZVlN6Nju92uGTNmaO7cubJarSosLNTYsWP173//u9R23AMBAABw+YiNjTW7BACVnGmX765bt06vv/66lixZoubNm2v//v2aMGGCGjVqpL///e8ltvPWgS01NZWDqEkYe/Mw9uZh7M3D2JuHsTdPdRz7wsJC5ebmml0GUC0FBgbKx6fk6GlaKJ0yZYoefPBB9evXT5IUHx+vo0eP6j//+U+poRQAAAAoi8LCQv32228KDg6WYRhmlwNUK3a7XVlZWQoKCioxmJp2T+nvv/8uq9XqNM1qtcpms5lUEQAAAKqi3NxcAilgEsMwFBwcXOqVCqadKb311lv1zDPPqHHjxmrevLm++uorLViwQHfddZdZJQEAAKCKIpAC5rnUz59pofSpp57SzJkz9cgjj+jnn39WSEiIhg4dyjtKAQAAAKAaMS2UBgUFadasWZo1a5ZZJQAAAAAATGZaKAUAAADMdOKEn44f994jVsLDbQoLK/Baf5eLHj16KC4uTnPmzDG7lGJWrFihcePG6fjx42aXUqURSnFZ87OckKWoah0kbNZwFdjCzC4DAIAq7/hxi8aOtV56QTdJTpbCyvArPjExUa+99pruuecePffcc07zpkyZonnz5qlbt25atWqVy+sMDg7W0qVL1adPH9cLqQSWLVumF198UWlpabJarYqIiFD37t3173//2+zSKoWEhAQdPXq0xPnt27fXO++848WKyoZQisuapei4rJljzS7DveolSwahFAAASBEREVq/fr1mzZqlwMBASedfcbNq1SpFRESYVldBQYH8/Py80terr76q8ePH68knn1THjh1VUFCggwcPavfu3V7p393OnTsnX19ft65z69atKioqkiR9/fXX6tevnz766COFh4dLUrHvlTe/f64w7ZUwAAAAAEoXHx+v6OhorV+/3jFt8+bN8vf31w033OC07L59+9S3b19FR0crMjJSt956q1NwS0hIkCQNHTpUwcHBjs+StGnTJnXs2FEhISFq0aKFpk+froKCAqe2SUlJeuCBB9SoUSMNHz5cK1asUHh4uLZt26brr79eDRs2VM+ePZWenu5od/jwYQ0aNEhXXnmlGjZsqA4dOui9994r0xhs2rRJvXr10r333qvo6Gg1b95ct99+u5588knHMklJSbr++uu1cuVKJSQkKDw8XKNGjVJBQYGWLFmi+Ph4NWnSRJMmTXJ6BWVWVpZGjhypxo0bKzQ0VH369NG3335bYi1ZWVnq1q2b7rjjDuXm5sput+vZZ59Vq1atFBoaqnbt2jmduf7xxx8VHBystWvXqlevXgoNDdUrr7xSpu13xRVXXKGQkBCFhIToL3/5iySpXr16jmlNmjTRiy++qLvvvlsNGzbUtGnTtH37dgUHByszM7NYvZ9//rlj2sGDB3XnnXcqIiJCTZs21X333aeMjAy31k8oBQAAACqxe+65RytWrHB8Xr58uYYMGVLsNRu//fabBg4cqE2bNmnLli1KSEjQgAEDHKFj69atkqR58+bp0KFDjs9btmzRiBEjNHz4cKWkpGj+/Pl66623NG3aNKf1L1y4UFdeeaU+/vhjTZkyRZKUn5+vuXPnav78+Xr//feVnZ2thx9+2NHmzJkz6tKli9avX68dO3aod+/euueee/Tdd9+5vP0hISHau3evU9i9mCNHjujdd9/VqlWrtGzZMr311lsaPHiw9u3bp3Xr1mnevHl64YUX9PbbbzvaJCYmau/evVq5cqW2bNmiGjVqqH///jp79myx9Z88eVLdu3dXWFiYXn/9dQUGBmrGjBl69dVXlZycrJSUFI0ZM0ZjxozR5s2bndpOnTpVw4YNU0pKinr06HHR+seMGaPw8PBSv0q7RPdSZs+era5du2rXrl0aNmyYS20ubPNVV12lLVu26M0339SZM2c0aNAgp3BfUVy+CwAAAFRiAwYM0OTJk/XDDz+oVq1a2rJli5566imnM4WS1LFjR6fPTz31lDZs2KAPP/xQAwcO1BVXXCFJqlOnjkJCQhzLJScna/To0br77rslSU2aNNETTzyh+++/X9OnT3eE33bt2ulf//qXo11KSooKCwuVnJys2NhYSdLo0aP1wAMPyGazyWKxKCEhwemM7NixY/Xee+/prbfe0qOPPurS9o8fP15ff/21WrVqpejoaLVp00adOnVS//79nS6DLSoq0oIFC1SnTh3FxcWpc+fO2rlzp7799lv5+fmpWbNmatu2rXbs2KE+ffrohx9+0KZNm/TOO++offv2kqTFixcrISFBa9as0d///nfHutPS0tS3b1917txZycnJslgsys3N1YIFC7Ru3Tq1a9dOkhQVFaW9e/dqyZIl6tatm6P9iBEjLnkf76RJkzR69OhSlwkry03Jf9K3b1+nbXIl4L700ku6+uqrNXXqVMe0xYsXKyoqSp9//rlat25d7nr+iFAKAAAAVGLBwcHq2bOnli9frjp16uiGG25QZGRkseVOnz6tmTNnavv27Tp9+rSKiop09uxZHTt2rNT1f/nll9q3b5+effZZxzSbzaazZ88qIyNDoaGhkqRrrrmmWFt/f39HIJWk0NBQnTt3TtnZ2apbt65yc3M1e/Zsbd68WSdPnlRhYaHy8vIUHx/v8vaHhobqgw8+0IEDB7Rz507t3r1bY8aM0cKFC7V582bVrFlT0vn7b+vUqeNo16BBAzVt2tTp3skGDRro9OnTkqRDhw7JYrHo2muvdcy/EGgPHjzomFZQUKBbb71VvXv3VnJysmP6oUOHlJeXp/79+zudtT537pwaNWrktA0XG7s/q1+/vurXr+/qsJSZKzX82Zdffqldu3Y57k39o8OHDxNKAQAAgOri7rvvVmJiogIDAzVp0qSLLpOYmKhTp07pySefVKNGjeTv76/evXs73Rt6MTabTePHj9ftt99ebN6Fs6uSHA9a+iMfH+c4cSGcXbi0c/Lkyfrwww81ffp0xcTEqGbNmho5cuQla7qYuLg4xcXFafjw4fr000912223af369RoyZIgkFXt4kGEYF63vwgOB7HZ7iX39MWT6+vqqU6dOev/993XkyBFH4Lywja+99lqxPxL8ud+Ljd2fjRkzRqtXry51mZSUlIv+QcIVf67BYjl/J+cfx6GwsNBpGZvNpq5du2rGjBnF1ufOAE0oBQAAACq5jh07ytfXV5mZmSXek5iSkqJZs2Y5Lhs9depUsQfS+Pr6OkLZBS1bttR3332n6Ohot9edkpKiu+66y3Hpal5eng4fPqyYmJgKrbd58+aSpNzc3Aqtw2azaffu3Y7Ld3NycnTgwAENHjzYsZxhGFq0aJFGjhypXr16aePGjYqMjFSzZs3k7++vo0ePFrt0ujw8ffnun134g8PJkycd/96/f7/TMi1bttT69esVGRnp9icG/xGhFAAAAKjkDMPQzp07Zbfb5e/vf9FlYmJitHr1arVp00a///67pkyZUuy1H40aNdK2bdvUvn17+fv7Kzg4WOPGjdPAgQMVGRmpvn37ysfHR99++6327t1b7GFHZRUTE6ONGzeqe/fu8vX11ezZs5Wfn1+mdTz88MMKDQ1Vhw4d1LBhQ2VkZCg5OVk1a9bUzTffXKHaunfvrjFjxuiZZ55RnTp1NH36dAUFBWnAgAFOy1osFj3//PMaOXKkevbs6Qimo0eP1uTJk2W329W+fXudOXNGe/bskcVi0T/+8Y8y1ePpy3f/LDo6WhEREZo1a5aeeOIJHTlyRHPmzHFaZtiwYVq6dKnuvfdePfTQQ7riiiuUnp6u9evXa8aMGQoKCnJLLYRSAAAAVEvh4Tb94RZBr/RXEZcKAPPnz9dDDz2km266SaGhoZowYYLT6z4kacaMGXrssccUHx+vsLAw7d+/X507d9bq1as1Z84czZ8/Xz4+PoqJiXE6W1heM2fO1OjRo9W9e3cFBwcrMTGxzKH0pptu0ooVK/TKK68oMzNTdevWVatWrbR+/Xo1bdq0QvUtXLhQEyZM0KBBg5Sfn6+2bdtq7dq1qlGjRrFlLRaLFi1apMTERPXq1Utvv/22HnvsMdWvX1/z58/XI488oqCgICUkJDg9EKqy8vX11UsvvaRHHnlEN9xwgxISEjRlyhQNHDjQsUxYWJg2b96sqVOnql+/fsrPz1dERIQ6depU4h9HysPIysoq+WLqaiw1NdXppm14T1nGPsC+R9bMsR6uqOzOnTN0rsC49IIXkRuYrAPp7dxckWtycn5T7doX/4UXHm5TWFjZ7/+AazjmmIexNw9jb57qNvbZ2dlOD8AB4H2l/RxyphTwgHMFhi7xoLsSpeUamvyk1b0FuSg/v4b8/S/ed3Ky5MbbGAAAAABJksXsAgAAAAAA1RehFAAAAABgGkIpAAAAAMA0hFIAAAAAgGkIpQAAAAAA0xBKAQAAAACmIZQCAAAAAEzDe0oBuMQwrNqzJ8DsMtwuPNymsLACs8sAAACotgilAFzy889SUpLV7DLcLjlZCgszuwoAgBn8LCdkKTrutf5s1nAV2Pil82c9evRQXFyc5syZY3YpxaxYsULjxo3T8ePe20/M8OOPP6ply5baunWrrrnmGq/3TygFAABAtWQpOi5r5ljvdVgvWTJcD6WJiYl67bXXdM899+i5555zmjdlyhTNmzdP3bp106pVq1xeZ3BwsJYuXao+ffq43KYyWLZsmV588UWlpaXJarUqIiJC3bt317///W+zS6s0kpKSNHv27GLTly9frp49e5pQkesIpQAAAEAlFRERofXr12vWrFkKDAyUJBUWFmrVqlWKiIgwra6CggL5+fl5pa9XX31V48eP15NPPqmOHTuqoKBABw8e1O7du73Sv7udO3dOvr6+Hll3bGysNm7c6DQtODjYI325k2kPOkpISFBwcHCxrzvvvNOskgAAAIBKJT4+XtHR0Vq/fr1j2ubNm+Xv768bbrjBadl9+/apb9++io6OVmRkpG699Van4JaQkCBJGjp0qIKDgx2fJWnTpk3q2LGjQkJC1KJFC02fPl0FBQVObZOSkvTAAw+oUaNGGj58uFasWKHw8HBt27ZN119/vRo2bKiePXsqPT3d0e7w4cMaNGiQrrzySjVs2FAdOnTQe++9V6Yx2LRpk3r16qV7771X0dHRat68uW6//XY9+eSTjmWSkpJ0/fXXa+XKlUpISFB4eLhGjRqlgoICLVmyRPHx8WrSpIkmTZokm83maJeVlaWRI0eqcePGCg0NVZ8+ffTtt9+WWEtWVpa6deumO+64Q7m5ubLb7Xr22WfVqlUrhYaGql27dk5nrn/88UcFBwdr7dq16tWrl0JDQ/XKK6+UafvLwsfHRyEhIU5f/v7+WrVqlTp16qSIiAg1bdpUQ4cO1U8//VTies6dO6dx48apefPmatCggeLj4/XEE0845hcUFOjxxx9XXFycGjZsqE6dOmnLli3lrtu0ULp161YdOnTI8bVt2zYZhqHbb7/drJIAAACASueee+7RihUrHJ+XL1+uIUOGyDAMp+V+++03DRw4UJs2bdKWLVuUkJCgAQMGKDMzU9L5/39L0rx583To0CHH5y1btmjEiBEaPny4UlJSNH/+fL311luaNm2a0/oXLlyoK6+8Uh9//LGmTJkiScrPz9fcuXM1f/58vf/++8rOztbDDz/saHPmzBl16dJF69ev144dO9S7d2/dc889+u6771ze/pCQEO3du9cp7F7MkSNH9O6772rVqlVatmyZ3nrrLQ0ePFj79u3TunXrNG/ePL3wwgt6++23HW0SExO1d+9erVy5Ulu2bFGNGjXUv39/nT17ttj6T548qe7duyssLEyvv/66AgMDNWPGDL366qtKTk5WSkqKxowZozFjxmjz5s1ObadOnaphw4YpJSVFPXr0uGj9Y8aMUXh4eKlfR48edXnc/qigoEATJ07Ujh07tGrVKmVmZuq+++4rcfnnn39e77zzjl566SXt3btXL7/8spo2beqY/8ADD2jnzp168cUXtWvXLg0aNEh33XWX9u/fX676TLt894orrnD6/OqrryooKIhQCgAAAPzBgAEDNHnyZP3www+qVauWtmzZoqeeesrpTKEkdezY0enzU089pQ0bNujDDz/UwIEDHf//rlOnjkJCQhzLJScna/To0br77rslSU2aNNETTzyh+++/X9OnT3eE33bt2ulf//qXo11KSooKCwuVnJys2NhYSdLo0aP1wAMPyGazyWKxKCEhwemM7NixY/Xee+/prbfe0qOPPurS9o8fP15ff/21WrVqpejoaLVp00adOnVS//79nS6DLSoq0oIFC1SnTh3FxcWpc+fO2rlzp7799lv5+fmpWbNmatu2rXbs2KE+ffrohx9+0KZNm/TOO++offv2kqTFixcrISFBa9as0d///nfHutPS0tS3b1917txZycnJslgsys3N1YIFC7Ru3Tq1a9dOkhQVFaW9e/dqyZIl6tatm6P9iBEjLnkf76RJkzR69OhSlwm7xNMZDx06pPDwcMfnyMhIpaSk6J577nFMi4qK0ty5c3Xttdfq+PHjTstfcPToUcXExKhdu3YyDEORkZFq27atpPNnv9euXauvvvpKkZGRju37+OOP9f/+3//T008/XWqNF1Mp7im12+169dVXNXDgQNWsWdPscgAAAIBKIzg4WD179tTy5ctVp04d3XDDDY4w8EenT5/WzJkztX37dp0+fVpFRUU6e/asjh07Vur6v/zyS+3bt0/PPvusY5rNZtPZs2eVkZGh0NBQSbroU1n9/f0dgVSSQkNDde7cOWVnZ6tu3brKzc3V7NmztXnzZp08eVKFhYXKy8tTfHy8y9sfGhqqDz74QAcOHNDOnTu1e/dujRkzRgsXLtTmzZsd+SEiIkJ16tRxtGvQoIGaNm3qdO9rgwYNdPr0aUnnA5zFYtG1117rmH8h0B48eNAxraCgQLfeeqt69+6t5ORkx/RDhw4pLy9P/fv3dzprfe7cOTVq1MhpG1x5om39+vVVv359V4flopo0aaI1a9Y4Pvv4nI97X3zxhWbPnq39+/crKytLdrtdknTs2LGLhtLBgwerb9++at26tW6++WZ16dJFXbp0kcVi0Zdffim73a7rrrvOqU1+fr46dOhQrrorRSjdunWrfvzxR6cEX5LU1FQvVOT9vuDM1bGPCslRjbx8D1dTdkVFfrLZynd1vK3Ipvx887appL7PnbMqP7/Iy9V4Xk7OWaWmpptdhiSOOWZi7M3D2JvHW2P/x8CC8rv77ruVmJiowMBATZo06aLLJCYm6tSpU3ryySfVqFEj+fv7q3fv3k73hl6MzWbT+PHjL3rF4h+vbrzwoKU/uhB6LrgQzi7ctzl58mR9+OGHmj59umJiYlSzZk2NHDnykjVdTFxcnOLi4jR8+HB9+umnuu2227R+/XoNGTJEkoo9PMgwjIvWV1R0/v8zF4LZxfwxZPr6+qpTp056//33deTIEUfgvLCNr732WrE/Evy534uN3Z+NGTNGq1evLnWZlJSUi/5B4gI/Pz9FR0c7TcvNzVW/fv100003afHixapfv74yMzN12223lfh9aNWqlb766itt2bJFn3zyiRITE3X11VfrzTfflM1mk2EY+uijj4qNeUBA+d5pXylC6dKlS/XXv/5VLVq0uOSy3jqwpaamchA1SVnGPsCeLWuBv4crKruiIoss5bxj22K1yN/fnG3Kz88vsW9fX8nfv1IcMtyqdm2fSvGzzjHHPIy9eRh78zD2l5+OHTvK19dXmZmZJd6TmJKSolmzZjkuGz116pQyMjKclvH19XWEsgtatmyp7777rliYcYeUlBTdddddjktX8/LydPjwYcXExFRovc2bN5d0PnBVZB02m027d+92XL6bk5OjAwcOaPDgwY7lDMPQokWLNHLkSPXq1UsbN25UZGSkmjVrJn9/fx09erTYpdPl4Y7Ldy8mNTVVmZmZmjx5sqKioiRJGzZsuGS7C7dW3n777Ro8eLBuueUWpaWlqUWLFrLb7crIyCj3mdE/M/1/mKdPn9a7777rdCocAAAAwP8YhqGdO3fKbreX+AfkmJgYrV69Wm3atNHvv/+uKVOmFHttS6NGjbRt2za1b99e/v7+Cg4O1rhx4zRw4EBFRkaqb9++8vHx0bfffqu9e/cWe9hRWcXExGjjxo3q3r27fH19NXv27DJfEfbwww8rNDRUHTp0UMOGDZWRkaHk5GTVrFlTN998c4Vq6969u8aMGaNnnnlGderU0fTp0xUUFKQBAwY4LWuxWPT8889r5MiR6tmzpyOYjh49WpMnT5bdblf79u115swZ7dmzRxaLRf/4xz/KVI87Lt+9mIiICPn7++vFF1/U8OHDdejQoWL3I//Z/PnzFRoaqoSEBPn6+mrNmjWqXbu2GjZsqJo1a+rOO+/UqFGjNHPmTLVs2VK//vqrduzYocaNG6t3795lrtH0ULpy5Ur5+/vrjjvuMLsUAAAAVCM2a7hUz3snRmzWcMl26eVKEhQUVOr8+fPn66GHHtJNN92k0NBQTZgwwfHk3QtmzJihxx57TPHx8QoLC9P+/fvVuXNnrV69WnPmzNH8+fPl4+OjmJgYp7OF5TVz5kyNHj1a3bt3V3BwsBITE8scSm+66SatWLFCr7zyijIzM1W3bl21atVK69evd3oibHksXLhQEyZM0KBBg5Sfn6+2bdtq7dq1qlGjRrFlLRaLFi1apMTERPXq1Utvv/22HnvsMdWvX1/z58/XI488oqCgICUkJDg9EMpsV1xxhRYtWqRp06Y5Xo8zc+ZM9evXr8Q2QUFBmjdvntLS0mQYhuPhTxfu312wYIGSk5M1ZcoU/fTTT6pbt67++te/6sYbbyxXjUZWVlbJF1N7mN1uV5s2bdS+fXvNmzfPrDIuistazFO2y3f3yJo51sMVld3vuRZd4pkCJUrLfVqTn7zh0gt6QGmX706cKCUlebkgL0hOLlKbNnlml8Exx0SMvXkYe/NUt7HPzs52egAOAO8r7efQtPeUStL27dv1ww8/aOjQoWaWAQAAAAAwiamX73bo0EFZWVlmlgAAAAAAMJGpZ0oBAAAAANUboRQAAAAAYBpCKQAAAADANIRSAAAAVHl2u2kvnACqvUv9/BFKAQAAUKUFBgYqKyuLYAqYwG63KysrS4GBgSUuY+rTdwEAAABP8/HxUVBQkHJycswuBaiWgoKC5ONTcvQklAIAAKDK8/HxUZ06dcwuA8BFcPkuAAAAAMA0hFIAAAAAgGkIpQAAAAAA0xBKAQAAAACmIZQCAAAAAExDKAUAAAAAmIZXwgCVTGxTu6ZP2mFK37YimyzWi/+t6uqrpemTyr/uYycjtfjlxuVfAQAAAKokQilQyQQG/KzowFmm9G2z2WWxGBedV/usFB1YgZWHPi2JUAoAAABnXL4LAAAAADANoRQAAAAAYBpCKQAAAADANIRSAAAAAIBpCKUAAAAAANMQSgEAAAAApiGUAgAAAABMQygFAAAAAJiGUAoAAAAAMA2hFAAAAABgGlND6cmTJzVy5EjFxMQoJCREbdu21Y4dO8wsCQAAAADgRT5mdZyVlaVu3brpuuuu0+rVq1WvXj39+OOPql+/vlklAQAAAAC8zLRQOm/ePIWGhmrx4sWOaVFRUWaVAwAAAAAwgWmX777zzjtq3bq17r33XjVt2lQ33HCDXnjhBdntdrNKAgAAAAB4mWlnStPT0/XSSy9p1KhReuihh7R//36NHz9ekjRixIgS26WmpnqrRK/2BWeujn1USI5q5OV7uJqyKyryk81Wvr/52O122Wzm/XGmpL7tdqNCddmKbMrPr3zfq5ycs0pNTTe7DEkcc8zE2JuHsTePt8Y+NjbWK/0AuHyZFkptNpuuueYaPf7445Kkli1bKi0tTUuWLCk1lHrrwJaamspB1CRlGfsAe7asBf4erqjsioosspTzOgTDMGSxGO4tyEU2m73Evg1DFarLYrXI37/yfa9q1/apFD/rHHPMw9ibh7E3D2MPoDJx+b/NO3fu1M8//1zi/MzMTO3cudPljkNCQtSsWTOnaVdeeaWOHTvm8joAAAAAAJc3l0Npr169tHXr1hLnb9u2Tb169XK54+uuu07ff/+907Tvv/9ekZGRLq8DAAAAAHB5czmUXuoBRAUFBbKU4XrFUaNG6bPPPlNycrLS0tL05ptv6oUXXtCwYcNcXgcAAAAA4PJW6j2lOTk5ys7Odnz+5ZdfdPTo0WLLZWVl6Y033lBYWJjLHf/1r3/VihUrNG3aNM2ZM0cRERGaNGkSoRQAAAAAqpFSQ+nChQv11FNPSTr/8JWJEydq4sSJF13Wbrdr8uTJZeq8W7du6tatW5naAAAAAACqjlJD6U033aSAgADZ7XZNmzZNd9xxhxISEpyWMQxDNWvW1DXXXKM2bdp4tFgAAAAAQNVSaii97rrrdN1110mS8vPz1atXL8XHx3ulMAAAAABA1efye0onTJjgyToAAAAAANVQiaH0tddekyTdddddMgzD8flSBg0a5J7KAAAAAABVXomhdNSoUTIMQ/369ZOfn59GjRp1yZUZhkEoBQAAAAC4rMRQ+uWXX0qS/Pz8nD4DAAAAAOAuJYbSRo0alfoZAAAAAICKsphdAAAAAACg+nL56buS9PHHH2vp0qVKT0/Xr7/+Krvd7jTfMAx98cUX7qwPAAAAAFCFuRxKFy1apMcee0xXXHGF2rRpo6uuusqTdQEAAAAAqgGXQ+mCBQvUvn17vfHGG46HHwEAAAAAUBEuh9LMzEw98sgjBFK43YkTfjp+/H+3N+fkRCk7O8CltnFRFgXmVr5bo202sysAAAAALg8uh9JWrVrpyJEjnqwF1dTx4xaNHWt1fM7PryF/f2spLf5n+iRD0YGeqqz8QkPNrgAAAAC4PLh8imnmzJlauXKlPvnkE0/WAwAAAACoRlw+U5qUlKTatWvr9ttvV0xMjCIjI2W1Op/NMgxDq1evdnuRAAAAAICqyeVQevDgQRmGoYiICOXn5+v7778vtoxhGG4tDgAAAABQtbkcSvfv3+/JOgDgsuNnOSFL0XG3rCsqJEcB9my3rKsibNZwFdjCzC4DAABUIy6HUgCAM0vRcVkzx7plXTXy8mUt8HfLuiqkXrJkEEoBAID3uBxKjx496tJykZGR5S4GAAAAAFC9uBxKW7Ro4dI9o7/88kuFCgIAAAAAVB8uh9L58+cXC6VFRUX68ccf9frrr6tBgwYaNmyY2wsEAAAAAFRdLofSIUOGlDjvoYce0s0336wzZ864pSgAAAAAQPVgccdKatWqpSFDhmjhwoXuWB0AAAAAoJpwSyiVJF9fX504ccJdqwMAAAAAVANuCaX79+/X888/r2bNmrljdQAAAACAaqLCT9/Nzs5WTk6OatWqpQULFrjccVJSkmbPnu00rUGDBvruu+9cXgcAAAAA4PLmciht3759sVBqGIaCg4MVHR2tfv36KTg4uEydx8bGauPGjY7PVqu1TO0BAAAAAJc3l0PpokWL3N+5j49CQkLcvl4AAAAAwOXBbQ86Ko/09HRdddVVatGihf75z38qPT3dzHIAAAAAAF7m8plSd2vTpo0WLlyo2NhY/fzzz5ozZ466du2qlJQU/eUvfymxXWpqqtdq9GZf1VlOTpTy82s4TcvPz3epra3IJpvN7omyKsRuN8pdl91uN3WbSuq7Itsknf9eufp99aacnLNKTU0vV9uokBzVyHPfNuW7cV3ldfa3HKVnVL9jH8d78zD25vHW2MfGxnqlHwCXL9NCaZcuXZw+t2nTRq1atdLKlSv14IMPltjOWwe21NRUDqJekp0dIH///91PnJ+fL39/f5faWqwWWSzFH8BlNsNQuesyDMO0bbLZ7CX2XZFtks5/r1z9vnpT7do+5f5ZD7Bny1rgnm3Kz8uXf4D54+MTVFuxtavXsY/jvXkYe/Mw9gAqE1Mv3/2jWrVqqXnz5kpLSzO7FAAAAACAl1SaUJqXl6fU1FQefAQAAAAA1YhLoTQvL0+zZ8/WRx995LaO//3vf2vHjh1KT0/Xnj17NHToUP3+++8aNGiQ2/oAAAAAAFRuLoXSgIAA/ec//9GxY8fc1vFPP/2kYcOG6W9/+5vuuece+fn56YMPPlCjRo3c1gcAAAAAoHJz+UFHCQkJbr3f8+WXX3bbugAAAAAAlyeX7ymdMmWKli1bps2bN3uyHgAAAABANeLymdJ58+YpODhYgwYNUsOGDRUVFaUaNZzfLWkYhlavXu32IgEAAAAAVZPLofTgwYMyDEMRERGSpCNHjhRbxjAq3/siAQAAAACVl8uhdP/+/Z6sAwAAAABQDVWa95QCAAAAAKqfMoXSoqIirV69Wg8++KAGDhyor7/+WpKUlZWl9evX6+TJkx4pEgAAAABQNbkcSrOzs9W1a1fdf//9euutt/TBBx8oMzNTkhQUFKTHHntML7zwgscKBQAAAABUPS6H0qlTp+rgwYNas2aNvvjiC9ntdsc8q9WqXr166YMPPvBIkQAAAACAqsnlUPrOO+9oxIgRuuWWWy76lN2YmBgdPXrUrcUBAAAAAKo2l0NpVlaWmjRpUuJ8u92ugoICtxQFAAAAAKgeXA6ljRo10oEDB0qcv3PnTjVt2tQtRQEAAAAAqgeXQ+mAAQO0bNky7dy50zHtwmW8ixcv1saNGzV48GD3VwgAAAAAqLJ8XF1wzJgx2rNnj3r37q2mTZvKMAxNmDBBv/zyizIyMtSjRw/df//9nqwVAAAAAFDFuBxKfX19tXr1aq1Zs0ZvvvmmDMNQYWGhWrZsqTvuuEN33nnnRR+AhMrDz3JClqLjZpdRTFyURdMn/W/fsRXZZLG6dhL/ypjfVcjrcQEAAIDLlsuh9IIBAwZowIABnqgFHmYpOi5r5lizyygmMNei6MD/fbbZ7LJYXPsDR80aE5TjoboAAAAAeF6ZQ6kkff31147Xv0RGRio+Pp6zpAAAAACAMitTKH3jjTf0+OOP66effpLdbpd0/mFHDRs21OOPP84ZVAAAAABAmbgcSlesWKEHH3xQsbGxmjp1qpo2bSq73a4ffvhBy5Yt0/3336+CggINGTLEk/UCAAAAAKoQl0Pp3Llz1bp1a23cuFEBAQFO84YPH67u3btr7ty5hFIAAAAAgMtcDqXHjx/XiBEjigVSSQoICNDAgQP1xBNPuLM2AFVIbFO7pk/aYXYZxcRF2RVgt5WrrdU46+ZqAAAAqh+XQ2nz5s114sSJEuf/9NNPatasmVuKAlD1BAb8rOjAWWaXUUxgrmRV+UKp6k10bzEAAADVkGsvg5Q0bdo0LV26VOvXry8274033tCyZcs0ffp0txYHAAAAAKjaXD5T+txzz6levXq67777NGHCBDVp0kSGYSgtLU2nT59WTEyM5s2bp3nz5jnaGIah1atXe6RwAAAAAMDlz+VQevDgQRmGoYiICEnnL9eVJH9/f0VERCg/P1+HDh1yasO7SwEAAAAApXE5lO7fv9+TdQAAAAAAqiGX7yn1tKefflrBwcF69NFHzS4FAAAAAOAllSKUfvbZZ1q6dKni4+PNLgUAAAAA4EWmh9Ls7GwNHz5czz33nIKDg80uBwAAAADgRaaH0oceekh9+vRRx44dzS4FAAAAAOBlLj/oyBOWLl2qtLQ0LV682OU2qampHqzIvL68ISokRzXy8s0uo5iiIj/ZbM5/H7HZ7C61tdvtLi/rTXa7Ue66zN6mkvquyDadb185v1dFRTbl5xWUq6313DkVufFnKr8S/Hye/S1H6RlV69jniqp2vL+cMPbm8dbYx8bGeqUfAJcv00Jpamqqpk2bpk2bNsnPz8/ldt46sKWmpla5g2iAPVvWAn+zyyimqMgiyx8yqc1ml8Xi2uuEDMNweVlvMgyVuy4zt6m0sa/INp1vXzm/V1arVf4B5fy58PWVT3nb/kl+Xn7563Ajn6Daiq1dtY59l1IVj/eXC8bePIw9gMrE5ct3W7ZsqXfffbfE+e+9955atmzpcse7d+9WZmamrr/+etWrV0/16tXTzp07tWTJEtWrV0/5+eafMQAAAAAAeJbLZ0qPHDmi3NzcEufn5ubq6NGjLnfco0cPXXPNNU7THnjgAcXExOjhhx8u09lTAAAAAMDlqUyX7xpGyZfeff/99woKCnJ5XcHBwcWetluzZk3VrVtXcXFxZSkLAAAAAHCZKjWUrly5Uq+99prjc3JyspYuXVpsuaysLB04cEDdunVzf4UAAAAAgCqr1FCam5urjIwMx+fs7GzZbDanZQzDUM2aNTV06FBNmDChQsW88847FWoPAAAAALi8lBpKhw8fruHDh0uSWrRooVmzZql79+5eKQwAAAAAUPW5fE/pV1995ck6AAAAAADVUJnfU/rbb7/p2LFj+vXXX2W324vNb9++vVsKAwAAAABUfS6H0l9//VXjx4/X+vXrVVRUVGy+3W6XYRj65Zdf3FogAKBszp0zdK6g5KellyZXFh1ID3BzRRUXHm5TWFiB2WUAAAAPcDmUjhkzRhs3btTw4cPVvn37Yq9zAQBUDucKDB07Vr62abmGJj9pdW9BbpCcLIWFmV0FAADwBJdD6Ycffqj7779fM2fO9GQ9AAAAAIBqxOLqgn5+foqJifFkLQAAAACAasblUNqnTx998MEHnqwFAAAAAFDNuBxKR48erZMnT2rkyJH67LPPdPLkSZ0+fbrYFwAAAAAArnL5ntLWrVvLMAx98cUXWr16dYnL8fRdAAAAAICrXA6l48aNk2GU7xUDAAAAAABcjMuhdOLEiZ6sAwAAAABQDbkcSv+oqKhI2dnZql27tnx8yrUKAKg0fs91+fZ6Jz61DBWWs+2fFRX5qajIPeuy2dyyGgAAAK8oU6Lct2+fpk2bpk8//VTnzp3T+vXr1bFjR2VmZioxMVEPPPCAOnbs6KlaAcDtCgulkyfL17a2r5RzzD112GwWWdyTSRUa6p71AAAAeIPL/wXavXu3unfvrsOHD+uuu+6S3W53zKtXr57OnDmjV1991SNFAgAAAACqJpdD6fTp0xUTE6P//ve/mjJlSrH5N954o/bs2ePW4gAAAAAAVZvLoXTfvn26++67FRAQcNGn8IaHhysjI8OtxQEAAAAAqjaXQ6nFYpGllBueMjIyVKNGDbcUBQAAAACoHlwOpa1atdJ777130XkFBQVas2aNrr32WrcVBgAAAACo+lwOpQ8//LA++eQTPfjgg9q/f78k6eTJk/rwww/Vu3dvHT58WI888ojHCgUAAAAAVD0uvxKmU6dOWrx4sR599FGtXLlSkpSYmCi73a46depoyZIl+tvf/uaxQgEAAAAAVU+Z3lPav39/de/eXVu3btUPP/wgm82mJk2aqHPnzqpVq5anagQAAAAAVFFlCqWSVLNmTfXo0cMTtQAAAAAAqhmX7yl999139eijj5Y4/9FHHy3xQUgAAAAAAFyMy6H0ueee0++//17i/Ly8PD377LMud/ziiy+qXbt2ioyMVGRkpLp06aLNmze73B4AAAAAcPlzOZQeOHBArVq1KnF+y5YtdfDgQZc7btiwoaZOnapt27Zp69at6tChg4YMGaKvv/7a5XUAAAAAAC5vLt9TWlhYqLNnz5Y4/+zZs8rPz3e54z/flzp58mS99NJL+uyzz3T11Ve7vB4AAAAAwOXL5TOlcXFx2rBhg2w2W7F5NptNGzZsUPPmzctVRFFRkd544w3l5ubq2muvLdc6AAAAAACXH5fPlI4cOVLDhg3ToEGDNHHiRF111VWSpG+//VazZs3S3r17tWjRojJ1/s0336hr167Ky8tTYGCgli9frvj4+FLbpKamlqmPivBmX94QFZKjGnmun832lqIiP9lszn8fsdnsLrW12+0uL+tNdrtR7rrM3qaS+q7INp1vz/fqUty1ropsk63IVqarXrwlJ+esUlPTPbb+qna8v5ww9ubx1tjHxsZ6pR8Aly+XQ2m/fv10+PBhJSUl6YMPPpAkGYYhu90uwzA0fvx4DRw4sEydx8bGavv27crOztaGDRuUmJiojRs3Ki4urtQ23pCamlrlDqIB9mxZC/zNLqOYoiKLLH/IpDabXRaL4VJbwzBcXtabDEPlrsvMbSpt7CuyTefb870qTVn2+0upyDZZrBb5+1e+40Tt2j4eOyZXxeP95YKxNw9jD6AyKdN7SseOHav+/fvr7bffVnp6uux2u5o0aaJevXopKiqqzJ37+fkpOjpaknTNNddo3759WrhwoebPn1/mdQEAKi62qV3TJ+0wu4xi4qLsCrAXv33EFTZruApsYW6uCAAAuItLofTs2bO68847NXDgQN19990aPXq0R4qx2WwqKCjwyLoBAJcWGPCzogNnmV1GMYG5klXlC6WqlywZhFIAACorlx50VKNGDX355ZcqKipyW8dPPPGEdu3apR9//FHffPONpk6dqh07dmjAgAFu6wMAAAAAULm5fPnuDTfcoF27dmno0KFu6TgjI0MjRozQqVOnVLt2bcXHx2vt2rXq3LmzW9YPAAAAAKj8XA6ls2fP1h133KHJkyfrvvvuU6NGjWSxuPxGmWLK+qReAAAAAEDV43Io/dvf/ia73a4FCxZowYIFslgs8vX1dVrGMAz99NNPbi8SAAAAAFA1uRxK+/btK8OofK9zAAAAAABcvlwOpVxuCwAAAABwt/LfFAoAAAAAQAWVKZQeOXJE//d//6dWrVopMjJSO3acf8F6ZmamHnnkEX3xxReeqBEAAAAAUEW5fPnuoUOHdOutt8pms6lNmzY6cuSI472l9erV02effab8/HzNnz/fY8UCAAAAAKoWl0Pp448/rqCgIH344YeyWq1q2rSp0/yuXbvqzTffdHd9AABIkn7PLd8dJ7my6EB6QInzc3KilJ1d8nxPCQ+3KSyswOv9AgBQ2bgcSnft2qWxY8eqQYMG+uWXX4rNj4yM1IkTJ9xaHAAAklRYKJ08Wb62abmGJj9pLXF+fn4N+fuXPN9TkpOlsDCvdwsAQKXj8p+dCwsLFRgYWOL8X3/9VVar93+pAwAAAAAuXy6H0ri4OG3fvv2i8+x2u95++221atXKXXUBAAAAAKoBl0NpYmKi3nrrLT311FOOy3dtNpu+++47/fOf/9Tnn3+u0aNHe6xQAAAAAEDV4/I9pf369dPRo0c1c+ZMzZo1yzFNkqxWq2bMmKEuXbp4pkoAAAAAQJXkciiVpIceekj9+/fXhg0blJaWJpvNpiZNmqh3795q3Lixp2oEAAAAAFRRlwyl+fn5evfdd5Wenq6//OUv6tatm0aNGuWN2gAAAAAAVVypoTQjI0Pdu3fX4cOHZbfbJUmBgYFatWqV2rdv75UCAQCoiNimdk2ftKPE+bYimyzW8r0DtSLiouwKsNvK1dZmDVeBjffJAACqhlJD6YwZM5Senq5Ro0apQ4cOSktL05w5czRu3Djt3LnTWzUCAFBugQE/KzpwVonzbTa7LBbDixWdF5gr5eeWr21uYLIOpDdxb0FuEB5uU1hYgdllAAAuM6WG0o8++kiDBg3SjBkzHNMaNGigYcOG6fjx4woPD/d4gQAAVEWFhdLJk+Vrm5ZraPKTle/d4MnJUhgncAEAZVTq9UoZGRlq27at07TrrrtOdrtdx44d82hhAAAAAICqr9RQWlRUpICAAKdpFz7n5eV5rioAAAAAQLVwyafvpqena+/evY7POTk5kqTU1FTVqlWr2PKtW7d2Y3kAAAAAgKrskqE0KSlJSUlJxaaPGzfO6bPdbpdhGPrll1/cVx0AAAAAoEorNZQuWLDAW3UAAAAAAKqhUkPp4MGDvVUHAAAAAKAa8v7bwgEAAAAA+P8RSgEAAAAApjEtlM6dO1edOnVSZGSkYmJiNHDgQB04cMCscgAAAAAAJjAtlO7YsUP33XefNm/erA0bNsjHx0e33367fv31V7NKAgAAAAB42SVfCeMp69atc/q8ePFiNWrUSCkpKbrttttMqqryO3HCT8ePl+9vCXFRFgXmVr4rtm02sysAAAAAYBbTQumfnTlzRjabTcHBwWaXUqkdP27R2LHWcrWdPslQdKCbC3KD0FCzKwAAAABglkoTSidMmKCEhARde+21pS6XmprqpYq825ercnKilJ9fo1xtbUU22Wx2N1dUcXa7UawuV+u02+2XzTa53tbcbSqp74ps0/n2fK8uxV3rqkzb5C6e3iYztrki22Qrsik/P9/NFVVcTs5Zpaaml6lNZfxdW114a+xjY2O90g+Ay1elCKWTJk1SSkqK3nvvPVmtpZ8F9NaBLTU1tVIeRLOzA+TvX74zpRarRRaL4eaKKs4w5FSXzWZ3uU7DMC6LbSpbW/O2qbSxr8g2nW/P96o0ZdnvL6WybJM7eXKb3Dn2ZVGRbbJYLfL393dzRRVXu7ZPmX53VtbftdUBYw+gMjE9lE6cOFHr1q3T22+/raioKLPLAQAAAAB4kamhdPz48Vq3bp02btyoK6+80sxSAAAAAAAmMC2Ujh07VqtWrdLy5csVHBysjIwMSVJgYKBq1aplVlkAAAAAAC8y7f0gS5Ys0W+//aY+ffqoWbNmjq/nnnvOrJIAAAAAAF5m2pnSrKwss7oGAAAAAFQSpp0pBQAAAACAUAoAAAAAMA2hFAAAAABgGkIpAAAAAMA0hFIAAAAAgGkIpQAAAAAA0xBKAQAAAACmIZQCAAAAAExDKAUAAAAAmIZQCgAAAAAwDaEUAAAAAGAaQikAAAAAwDSEUgAAAACAaQilAAAAAADTEEoBAAAAAKYhlAIAAAAATEMoBQAAAACYhlAKAAAAADANoRQAAAAAYBpCKQAAAADANIRSAAAAAIBpCKUAAAAAANMQSgEAAAAApiGUAgAAAABMQygFAAAAAJjG1FC6c+dO3XXXXbrqqqsUHBysFStWmFkOAAAAAMDLTA2lubm5iouL06xZs1SjRg0zSwEAAAAAmMDHzM67du2qrl27SpJGjRplZikAAAAAABNwTykAAAAAwDSmniktj9TU1CrZl6tycqKUn1++S51tRTbZbHY3V1RxdrtRrC5X67Tb7ZfNNrne1txtKqnvimzT+fZ8ry7FXeuqTNvkLp7eJjO2uSLbZCuyKT8/380VVVxOzlmlpqaXqU1l/F1bXXhr7GNjY73SD4DL12UXSr11YEtNTa2UB9Hs7AD5+1vL1dZitchiMdxcUcUZhpzqstnsLtdpGMZlsU1la2veNpU29hXZpvPt+V6Vpiz7/aVUlm1yJ09ukzvHviwqsk0Wq0X+/v5urqjiatf2KdPvzsr6u7Y6YOwBVCZcvgsAAAAAMM1ld6YUAABUToZh1Z49AS4vn5MTpexs15c3Q3i4TWFhBWaXAQBVmqmh9MyZM0pLS5Mk2Ww2HTt2TF999ZXq1q2ryMhIM0sDAABl9PPPUlKS67eY5OfXKPctKd6SnCyFhZldBQBUbaaG0s8//1y9evVyfE5KSlJSUpIGDRqkRYsWmVgZAACVV2xTu6ZP2mF2GcVcfbU0fZLry9uKbLJYz99JFBJWWxkncjxUWfnFRdkVYLeVu73NGq4CG6kWAEpjaii98cYblZWVZWYJAABcdgIDflZ04Cyzyyim9lkpOtD15f/4kKna9SYoMKfybVNgrmRV+UOp6iVLBqEUAErDg44AAAAAAKYhlAIAAAAATEMoBQAAAACYhlAKAAAAADANoRQAAAAAYBpCKQAAAADANIRSAAAAAIBpCKUAAAAAANMQSgEAAAAApiGUAgAAAABMQygFAAAAAJiGUAoAAAAAMA2hFAAAAABgGkIpAAAAAMA0hFIAAAAAgGkIpQAAAAAA0/iYXQAAAEBl9ntu+f+GnyuLDqQHuLEa9wgICDO7BABwIJSWIOyKfAXY95hdRjFxURZNn2SUq+2VMb+r8KSbCwIAoAorLJROVuB3Z1quoclPWt1XkJtMmeJvdgkA4EAoLYG/9ZSsmdPMLqOYwFyLogPL17ZmjQnKcW85AAAAAFAh3FMKAAAAADANoRQAAAAAYBpCKQAAAADANIRSAAAAAIBpCKUAAAAAANOYHkqXLFmiFi1aKCQkRB07dtSuXbvMLgkAAAAA4CWmhtJ169ZpwoQJeuSRR/TJJ5/o2muv1YABA3T06FEzywIAAAAAeImpoXTBggUaPHiwhg4dqmbNmmnOnDkKCQnRyy+/bGZZAAAAAAAvMS2UFhQU6IsvvtDNN9/sNP3mm2/Wf//7X5Oq+p+goNpml1BtWSyG2SVUW4y9eRh78zD25mHszVO7dpDZJQCAg5GVlWU3o+MTJ07oqquu0jvvvKP27ds7ps+ePVtr1qzRnj17zCgLAAAAAOBFpj/oyDCc/0pqt9uLTQMAAAAAVE2mhdJ69erJarXq1KlTTtN//vln1a9f36SqAAAAAADeZFoo9fPzU6tWrbR161an6Vu3blXbtm1NqgoAAAAA4E0+Znb+wAMP6P7771fr1q3Vtm1bvfzyyzp58qTuvfdeM8sCAAAAAHiJqfeU3nHHHUpKStKcOXN04403KiUlRatXr1ajRo280v/OnTt111136aqrrlJwcLBWrFjhND8xMVHBwcFOX7fccotXaqvK5s6dq06dOikyMlIxMTEaOHCgDhw44LSM3W5XUlKSmjdvrtDQUPXo0UPffvutSRVXHa6MPfu9Z7z44otq166dIiMjFRkZqS5dumjz5s2O+ezznnOpsWef946nn35awcHBevTRRx3T2O+942Jjz34PoDIx/UFHw4YN0/79+3Xq1Clt27bN6Um8npabm6u4uDjNmjVLNWrUuOgyN910kw4dOuT4WrNmjdfqq6p27Nih++67T5s3b9aGDRvk4+Oj22+/Xb/++qtjmWeffVYLFizQ7Nmz9dFHH6l+/frq27evfvvtNxMrv/y5MvYS+70nNGzYUFOnTtW2bdu0detWdejQQUOGDNHXX38tiX3eky419hL7vKd99tlnWrp0qeLj452ms997XkljL7HfA6g8TL1812xdu3ZV165dJUmjRo266DL+/v4KCQnxZllV3rp165w+L168WI0aNVJKSopuu+022e12LVq0SA899JD69OkjSVq0aJFiY2O1du1aLu+ugEuN/QXs9+7Xo0cPp8+TJ0/WSy+9pM8++0zx8fHs8x5U2thfffXVktjnPSk7O1vDhw/Xc889p6eeesoxnWO955U09hew3wOoLEw/U1rZffrpp2ratKlat26t//u//9Pp06fNLqnKOXPmjGw2m4KDgyVJP/74ozIyMnTzzTc7lqlRo4batWun//73vyZVWTX9eewvYL/3rKKiIr3xxhvKzc3Vtddeyz7vRX8e+wvY5z3nQujs2LGj03T2e88raewvYL8HUFlU6zOll3LLLbeoV69eaty4sY4cOaIZM2aod+/e+vjjj+Xv7292eVXGhAkTlJCQ4PgPYkZGhiQVezVQ/fr1deLECa/XV5X9eewl9ntP+uabb9S1a1fl5eUpMDBQy5cvV3x8vOM/4OzznlPS2Evs8560dOlSpaWlafHixcXmcaz3rNLGXmK/B1C5EEpL0a9fP8e/4+Pj1apVKyUkJGjz5s3q3bu3iZVVHZMmTVJKSoree+89Wa1Wp3mGYTh9ttvtxaah/Eoae/Z7z4mNjdX27duVnZ2tDRs2KDExURs3bnTMZ5/3nJLGPi4ujn3eQ1JTUzVt2jRt2rRJfn5+JS7Hfu9+row9+z2AyoRQWgZhYWFq2LCh0tLSzC6lSpg4caLWrVunt99+W1FRUY7pF+5vOXXqlCIiIhzTf/7552J/UUf5lDT2F8N+7z5+fn6Kjo6WJF1zzTXat2+fFi5cqLFjx0pin/ekksZ+/vz5xZZln3eP3bt3KzMzU9dff71jWlFRkXbt2qWXX35ZKSkpktjvPeFSY//TTz8VOxvKfg/ATNxTWgaZmZk6ceIEDwVwg/Hjx2vt2rXasGGDrrzySqd5jRs3VkhIiLZu3eqYlpeXp08//VRt27b1dqlVTmljfzHs955js9lUUFDAPm+CC2N/Mezz7tGjRw/t2rVL27dvd3xdc8016tevn7Zv366mTZuy33vIpcb+YmdP2e8BmKlanyk9c+aM4y+CNptNx44d01dffaW6deuqbt26mjVrlnr37q2QkBAdOXJE06ZNU/369dWzZ0+TK7+8jR07VqtWrdLy5csVHBzsuK8oMDBQtWrVkmEYSkxM1NNPP63Y2Fg1bdpUycnJCgwMVP/+/U2u/vJ2qbE/c+YM+72HPPHEE+ratavCw8N15swZrV27Vjt27NDq1avZ5z2stLFnn/ecC+++/KOaNWuqbt26iouLkyT2ew+51Niz3wOobKp1KP3888/Vq1cvx+ekpCQlJSVp0KBBmjt3rg4cOKDXX39d2dnZCgkJ0Y033qhXXnlFQUFBJlZ9+VuyZIkkOV4BcMH48eM1ceJESdK//vUvnT17Vo8++qiysrLUunVrrVu3jrGvoEuNvdVqZb/3kIyMDI0YMUKnTp1S7dq1FR8fr7Vr16pz586S2Oc9qbSxP3v2LPu8idjvzcGxHkBlY2RlZdnNLgIAAAAAUD1xTykAAAAAwDSEUgAAAACAaQilAAAAAADTEEoBAAAAAKYhlAIAAAAATEMoBQAAAACYhlAKAAAAADANoRQAAAAAYBpCKQAAAADANP8fh7CaYuPyRHEAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"smoking_and_age = births[['Maternal Smoker', 'Maternal Age']]\n",
"\n",
"unit = ''\n",
"\n",
"fig, ax = plt.subplots(figsize=(10,5))\n",
"\n",
"ax.hist(smoker['Maternal Age'], density=True, label='Maternal Smoker = True', color='blue', alpha=0.8, ec='white', zorder=5)\n",
"\n",
"ax.hist(non_smoker['Maternal Age'], density=True, label='Maternal Smoker = False', color='gold', alpha=0.8, ec='white', zorder=10)\n",
"\n",
"y_vals = ax.get_yticks()\n",
"\n",
"y_label = 'Percent per ' + (unit if unit else 'unit')\n",
"\n",
"x_label = ''\n",
"\n",
"ax.set_yticklabels(['{:g}'.format(x * 100) for x in y_vals])\n",
"\n",
"plt.ylabel(y_label)\n",
"\n",
"plt.xlabel(x_label)\n",
"\n",
"plt.title('')\n",
"\n",
"ax.legend(bbox_to_anchor=(1.04,1), loc=\"upper left\")\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The observed difference between the average ages is about $-0.8$ years."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.8076725017901509"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"observed_age_difference = difference_of_means(births, 'Maternal Age', 'Maternal Smoker')\n",
"observed_age_difference"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Remember that the difference is calculated as the mean age of the smokers minus the mean age of the non-smokers. The negative sign shows that the smokers are younger on average.\n",
"\n",
"Is this difference due to chance, or does it reflect an underlying difference in the population?\n",
"\n",
"As before, we can use a permutation test to answer this question. If the underlying distributions of ages in the two groups are the same, then the empirical distribution of the difference based on permuted samples will predict how the statistic should vary due to chance."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"age_differences = np.array([])\n",
"\n",
"repetitions = 5000\n",
"for i in np.arange(repetitions):\n",
" new_difference = one_simulated_difference(births, 'Maternal Age', 'Maternal Smoker')\n",
" age_differences = np.append(age_differences, new_difference)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The observed difference is in the tail of the empirical distribution of the differences simulated under the null hypothesis. "
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Observed Difference: -0.8076725017901509\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAFuCAYAAAB9QTkMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABPRklEQVR4nO3dd1wT9/8H8FfYqGAUAQdDBSqiViqu4qg46gJx41bqxFUXCrVfrauAonVRqnVU6qgLFUfdqAiKoyLWiXWigsp0QBTI7w8f5GdMAgESAvH1fDx8tLn73N37PrnAi7vPXQTp6eliEBEREWkJHU0XQERERKRKDDdERESkVRhuiIiISKsw3BAREZFWYbghIiIircJwQ0RERFqF4YbKHB8fHwiFQkRFRUlNFwqF6N69u9q2GxAQIHe79EF+/2zZskXTpSjUvXt3CIVCPHz4UNOllAp5nwkexwXbsmULhEIhAgICNF0KAL5f6sJw85kSCoVS/6pWrQpbW1t06dIFGzduRG5urqZLVLny8Mv5U/k1F/SDOCoqSu3Br6wo6+GlUaNGEAqFqFWrFp49eya3zXfffVemfplp2zFW1sILaYaepgsgzZo1axYAIDc3F/fv38eBAwdw/vx5nDp1Cps2bdJwddIuXLgAY2Njta1/zJgx6NOnD6ysrNS2Dfo8vHnzBgsXLkRISIimS6Eyjj931IPh5jPn7+8v9fr69evo2LEj9u3bh5iYGLi6umqoMllffPGFWtdvZmYGMzMztW6DPg92dnbYtm0bxo4diy+//FLT5VAZxp876sHLUiSlQYMGaNWqFQDg8uXLAP7/lLSPjw9u3bqFIUOGoG7duhAKhYiPj5csu2/fPnh6eqJ27dqwsLBAkyZN8NNPPyEzM1Putk6dOoWuXbuiZs2aqF27NgYNGoTbt28rrE3RafHc3FyEhYWha9eusLW1haWlJb788kuMGjUKV65cAfDhckZQUBAAYMKECVKX5PIvcRR07fvMmTPo168f6tSpAwsLCzRu3BizZs3CixcvZNp+PGZo3759aN++PWrUqIHatWvD29sbT548UbiPqlTcOuLi4iR/SVpbW8PT0xOxsbEFbis5ORl+fn5o0qQJLC0tYWtri169euH06dMybT++bBAbG4vevXvD1tYWQqEQ6enpCrchFAoRHR0NAGjcuLHk/WvUqJHc9hs3boSrqyssLS3h4OCAyZMnK1x/UepXxk8//YS8vDz8+OOPSi9T0GWfsjgu4/fff4dQKERgYKDc+ZmZmahZsyYaNGggucz98aXhw4cPo1OnTpLP/4gRI3D//n2560pOToavry8aN24MCwsL1KlTB/3798fZs2el2vn4+GDChAkAgKCgIKnPuby+i4+PR//+/WFjY4MaNWqga9euOH/+vNwa8vLyEBYWhs6dO8PGxgaWlpb4+uuvsWzZMrx7906mfVRUFLy8vNCgQQNYWFjA3t4e7dq1w+zZsyEW//+3Hil6b5VdnuTjmRtS2v379/Htt9+iXr16GDBgADIyMlChQgUAwPTp07F+/XrUqlUL7u7uEAqFuHTpEpYvX46jR4/iyJEjMDExkaxr37598Pb2hr6+Pnr27ImaNWvi/Pnz6NSpExo2bKh0Te/evcOgQYNw/PhxVK9eHb169UKVKlWQmJiIqKgo2NnZ4auvvsKgQYMAANHR0ejWrZvUL8TKlSsXuI2NGzdi2rRpMDY2hqenJ6pXr47Y2FisWbMGBw8exN9//w1ra2uZ5davX4+///4b3bp1Q6tWrXDp0iXs2bMH165dQ3R0NAwNDZXez5IoSh2xsbHo2bMnRCIRPDw8YGdnh+vXr8PDwwNt27aVu/7r16+jV69eePHiBdq3b49u3bohNTUVBw8eRM+ePbFy5UoMHTpUZrkLFy5g2bJlcHV1xbBhw/Ds2TPo6uoq3I9Zs2Zh69atePz4McaNGyd53+S9f3PnzsXJkyfRpUsXuLm5ISoqCmFhYbh79y4OHTqkkvoL0qlTJ7i5uSEyMhKHDx9Gly5dirR8eTBgwADMnz8ff/75J3x9fWXeu7/++gtv377F5MmTZebt378fx48fh4eHB9q0aYP4+Hjs3bsXUVFROHr0KOzs7CRtHz58iK5du+Lp06do1aoVevfujaSkJOzduxfHjx/H8uXLMWzYMAAf/ojJyMjAoUOH0KpVK7Ru3VqyHhsbG6ka4uLisHLlSrRo0QLDhg1DYmIiIiIi4OnpiTNnzqBevXqStjk5ORgyZAgOHz4Me3t79OnTB4aGhoiOjsb8+fNx+vRp7N69G3p6H36lHj16FF5eXjAxMUHXrl1Rq1YtpKen47///sOaNWswb948SVt5Sro8MdzQJ27evCn567hJkyZS886fP49p06Zhzpw5UtO3b9+O9evXw93dHb///rvUuJglS5Zg0aJFCAgIwM8//wwAeP36NaZMmQKBQICDBw+iadOmkvb/+9//sGrVKqXrDQoKwvHjx9GuXTts3bpVEraAD2d08s+sDB48GI8ePUJ0dDS6d++OwYMHK7X+R48eYdasWahQoQKOHz+O+vXrS+YtXLgQwcHBmD59Onbs2CGz7MmTJ3H69Gk4OjpKpo0aNQq7du3CwYMH0bt3b6X3sySUrUMsFmPixInIysrCpk2b4OnpKWn/+++/w9fXV2bdubm5GD58ODIyMrB//36pXyZJSUno0KEDfH190blzZ1hYWEgtGxkZieXLl2PEiBFK7Ye/vz/Onj2Lx48fw8fHB7a2tgrbXr58GefOnUOtWrUAfPjl5OHhgZiYGFy6dElyzJWk/sIsWLAAbdu2xZw5c9CxY8cy/8vo7NmzCgfhPnr0SGaaiYkJvLy8sG7dOhw+fFjmrNMff/wBPT09SfD42OHDh7F9+3Z07txZMm3VqlX43//+B19fX4SHh0umT506FU+fPoWfnx/8/Pwk0ydOnIiOHTvC19cX7du3h5WVFdzd3SXhpnXr1jKX3T925MgRrFmzBl5eXpJpGzduxNSpU7FmzRosW7ZMMv2XX37B4cOHMXr0aAQGBkrCWl5eHqZOnYpNmzZh3bp1GDduHAAgLCwMYrEY+/fvR+PGjaW2m5qaWuixUNLliZelPnsBAQEICAjAwoULMXr0aLi5uSErKwvu7u6Sy1P5LCwsJAOQP/brr79CV1cXq1atkhnwO23aNJiZmUn98j906BDS0tLQu3dvqWADADNnzoSpqalStefm5mLdunUwNDTEihUrpIINAOjq6qJ69epKrUuRHTt24N27dxg5cqRUsAEAX19f1KhRA0ePHsXTp09llh07dqxUoACA4cOHAwD++eefEtVVFMrWERsbi4SEBLRo0UIq2ADAyJEjUbduXZl1Hz16FHfv3sXIkSOlggEAVK9eHZMmTUJ2djb27dsns2zDhg2VDjZFNXPmTEmwAQA9PT0MGTIEgPQ+l6T+wjRs2BBDhgzBnTt3sHHjxmLuSemJjo5GUFCQ3H/btm2Tu8yoUaMAQGb/zp8/jxs3bqBLly6oWbOmzHJt27aVCjbAh0tKVlZWOHnypOTz9OTJE5w8eRI1a9bEtGnTpNo3aNAA3333HUQiEbZv317k/f3666+lgg0ADBkyBHp6elLHSF5eHn777TeYm5sjICBA6iyUjo4O5s+fD4FAIFWDjs6HX62f/kwCgKpVqxZaW0mXJ565+ezlj0MRCAQwMTFB48aN0a9fP7m/dBo2bChzKSUrKwvx8fGoUqUKfvvtN7nbMDAwwLNnz5CamoqqVavi6tWrACATnoAPfw1++eWXMtfS5blz5w4yMjLQuHHjAv+KL4n8WuVdkjE0NETLli2xZ88exMfHy/wQd3Z2llkm/xduQWNLVE3ZOgp6X3R0dNCyZUvcu3dPanr+WJzExES5f/Xnt79z547MvE+DrSopu88lqV8Zs2fPRnh4OAIDA9G/f/9CL4Fq0qxZsxSe6YiKioKHh4fMdEdHR7Ru3RonT57EgwcPULt2bQD/H3ZGjhwpd33yjjE9PT20aNECiYmJks9T/pi+li1bwsDAQGaZdu3aISQkRHLsFoW8Y0RfXx8WFhZSx8jdu3eRkpKCOnXqYMmSJXLXZWxsjISEBMnr/v37IyIiAh06dECvXr3Qpk0bNGvWTOmfUyVdnhhuPntF+SUr77R8WloaxGIxUlNTJUFJkdevX6Nq1aqSAcbm5uZKb0eejIwMAJD7l6Gq5NeqqCZLS0updh+TdwYq/68+ZZ8jlP8XXF5ensI2+fPy2xa3juK8L6mpqQCAiIgIREREKKzxzZs3Sq1PVZTd55LUrwxLS0tMnjwZP//8M5YuXYr58+cXaz1l2ejRo3H27Fls2rQJc+fORVpaGvbt24e6deuiXbt2cpdR9N7nH3v5x2JJPn+FUXSGWFdXV+4xcv/+/UJ/xuVzd3fH7t27sWrVKmzbtk3yWA0nJyfMmjVL5syoqpcnhhsqAoFAIDMt/weEk5MTYmJilFpP/jLy7jQCgOfPnyu1nvy/ghU9LE0V8mtVVFNycrJUO3VtPy0tTWGb/B++JT0rUJz3JX+ZsLAw9OjRo0jbk3c8lbaS1K+sSZMmYdOmTVizZo3CMxnAh/5QFHrzg3xZ1L17d9SsWRObN2+Gv78/tm7diuzsbIwYMULhe6zo85R/7OW/L5r+/H287i5duuCvv/5SerkOHTqgQ4cOyMrKwuXLl3H8+HGsX78eI0aMkBnfpY7lP3ccc0MlUqlSJTg5OSEhIQEpKSlKLZM/QC5/4PLHXr16JXV7eUG++OILVK5cGTdv3sTjx48LbV/UsyYf1yrvNlKRSCS5rPHpoD9Vyb9zTNHtqR/PK8pdZvIU9L7k5eXJraFZs2YAgHPnzpVo28r6eCCnKpRG/cbGxvjxxx8hEokwb948he2EQiESExPlzst/pEFZpKenh+HDh+PFixc4cOAANm3aBENDwwIH7cs7xnJyciSfp/xnA+X/NzY2Vu7t1vm36n98iak4n/OC5P+cuXz5stwaCmNsbIzWrVvjp59+woIFCyAWi2Xu2FPn8p8rhhsqsQkTJuD9+/cYP3683DMMr169wqVLlySvu3XrBqFQiPDwcKnpALB48WKlTzHr6upi9OjREIlEmDJlCrKysqTm5+bmIikpSfI6/0FZin6ByNO/f38YGBhg/fr1MuMuli1bhqdPn+Lbb79FjRo1lF5nUbi6uqJOnTr4999/ERYWJjM/Pj4emzdvhp6enszgyKJq0aIFHBwcEBsbKzOAdv369TLjbYAP72XdunWxceNGhT9wr169Kjm7VFL576EyYVYZpVX/gAED4OzsjPDwcMTFxclt06xZMyQmJuLo0aNS0zdt2lToc4Y0bcSIEdDX18cPP/yAO3fuwNPTs8AH0505cwZHjhyRmhYaGorExES4ublJLjXXqlULHTp0wJMnT7BixQqp9jdv3sSGDRtgaGiI/v37S6YX53NeED09PYwbNw4vXrzAjBkz8PbtW5k2KSkpUn+UnTp1Sm67/DNNRkZGBW6zpMsTL0uRCgwePBhXr17F2rVr4ezsjA4dOsDGxgYZGRl49OgRYmJi4Obmhq1btwL4cLZnxYoV8Pb2Rvfu3dGrVy/UrFkT586dw40bN+Dq6qr0Ja6ZM2fiypUrOHHiBJo0aYIuXbqgSpUqePr0KaKiojBkyBDJIMlvvvkGOjo6+O2335CWlia5jj9mzBiFl3RsbGwQFBSEadOmwc3NDT179oSlpSViY2MRHR2NWrVqYenSpSroRfl0dXWxdu1a9OnTB5MnT8a2bdvQtGlT6Ovr486dOzhy5Ahyc3OxePFi1KlTp0TbEggEWLVqFXr16gVvb2+p59xERkaiY8eOOH78uNQy+vr62Lx5M3r37o1BgwahadOmaNy4MSpWrIgnT54gPj4eCQkJOHPmjEru8nBzc8OePXvw/fffw9PTExUrVkTlypUxZsyYYq2vtOoXCARYuHAh3N3d5YZEAJg8eTKOHz+OIUOGoGfPnjA3N0dcXBzi4uLQuXNnmTBQllhaWsLd3R179uwB8OH7swrStWtXDB48GD169EDt2rURHx+P48ePo2rVqggODpZqu2zZMnTp0gWLFi3CmTNn0KxZM8lzbrKysrBixQqpry5o3rw5KlWqhPDwcBgYGMDKygoCgQBeXl4yz7pRlq+vL27cuIGwsDAcPXoUbdu2Ra1atfDy5Uvcv38f58+fx6hRoyRnmn788Uc8evQIrVq1go2NDYyMjHD9+nWcOHECVatWldytqEhJlyeGG1KRxYsX49tvv8X69etx9uxZpKWloXLlyqhZsyZGjhyJfv36SbX39PTE7t27ERQUhH379sHAwACurq44duwYfvnlF6XDjYGBAXbs2IFNmzZh27Zt2LlzJ3JycmBpaYlWrVqha9eukrb29vZYv349VqxYgc2bN0vO9BR2F4u3tzfq1q2LVatW4eDBg3jz5g1q1KiBMWPGYMaMGWodGAt8+Iv+7NmzWL16NSIjI7Fu3Trk5ubC3NwcPXr0wNixY9G8eXOVbKtly5b4+++/sWDBApw4cQInTpyAi4sLDhw4gBMnTsiEG+DDeKvo6GiEhobi0KFD2LZtG8RiMSwtLeHo6IhJkybBwcFBJfUNGTIET548wY4dOxASEoL379/D2tq62OGmNOtv3bo1unXrpvAMUevWrbF9+3YEBgYiIiJC6jOxb9++Mh1ugA/vzZ49e+Dk5ISWLVsW2Nbd3R0jRoxAcHAwDh8+DH19fXh6emLu3LkyjxywtbXFqVOnJG3Pnz+PihUrolWrVpg8eTLatGkj1b5y5crYsmULAgICEB4ejtevXwP4cGwXN9zo6ekhLCwMu3fvxpYtW3Ds2DHJDRLW1taYOnUqBgwYIGk/ffp0HDx4EFeuXJFc0q5ZsyZ8fHwwfvz4Qr9HqqTLEyBIT0/nc5yJiKhEli5digULFiA4OFjy/JtPBQQEICgoCCEhIUo/SJOoODjmhoiISuT169f4/fffYWpqWuKxX0SqwMtSRERULH///TeuXLmCY8eOISkpCXPnzpX6DjkiTWG4ISKiYomIiMC2bdtgYWGBKVOmYPLkyZouiQgAx9wQERGRluGYGyIiItIqDDdERESkVRhuiIiISKsw3Cjp46+zp5Jjf6oO+1K12J+qw75ULfan8hhuiIiISKsw3BAREZFWYbghIiIircJwQ0RERFqF4YaIiIi0CsMNERERaRWGGyIiItIqDDdERESkVRhuiIiISKsw3BAREZFWYbghIiIiraKn6QKIiFTtmUiAJ2/yirVspoEFMlLFKq6ocLUq6qCGYelvl0gbMdwQkdZ58iYPM86+KNayIlE2DA1FKq6ocMGtzVHDUFDq2yXSRrwsRURERFqF4YaIiIi0CsMNERERaRWGGyIiItIqDDdERESkVRhuiIiISKsw3BAREZFWYbghIiIircJwQ0RERFqF4YaIiIi0CsMNERERaRWGGyIiItIqDDdERESkVRhuiIiISKsw3BAREZFWYbghIiIircJwQ0RERFqF4YaIiIi0CsMNERERaRWGGyIiItIqDDdERESkVTQabqKjozFgwADUr18fQqEQW7ZskZovFosREBAAR0dHVK9eHd27d8fNmzel2ohEIvj6+qJu3bqoWbMmBgwYgCdPnpTmbhAREVEZotFw8+bNGzg5OSEwMBDGxsYy81esWIGQkBAEBQXh5MmTMDc3R69evfDq1StJG39/f+zfvx/r16/HoUOH8OrVK3h5eSE3N7c0d4WIiIjKCI2Gm2+//RZz5syBp6cndHSkSxGLxQgNDcWUKVPg6ekJJycnhIaG4vXr19i1axcAICMjA3/++Sfmz58PNzc3ODs7Y82aNbh+/TpOnTqlgT0iIiIiTSuzY24ePnyI5ORktG/fXjLN2NgYrq6uiI2NBQDExcXh/fv3Um2srKxQr149SRsiIiL6vOhpugBFkpOTAQDm5uZS083NzfHs2TMAwPPnz6GrqwszMzOZNs+fP1e47oSEhGLVVNzlSD72p+qwL6VlGlhAJMou9vIlWba4Ml9lIiFF8c+t8orHpmqxPz9wcHAocH6ZDTf5BAKB1GuxWCwz7VOFtSmsU+RJSEgo1nIkH/tTddiXsjJSxTA0FBVrWZEoG4aGRiquqHCmJqZwqFq51LerTjw2VYv9qbwye1nK0tISAGTOwLx8+VJyNsfCwgK5ublISUlR2IaIiIg+L2U23Nja2sLS0hKRkZGSadnZ2Th37hxatGgBAHB2doa+vr5UmydPnuD27duSNkRERPR50ehlqdevX+PevXsAgLy8PCQmJiI+Ph5VqlSBtbU1fHx8sHTpUjg4OMDe3h7BwcGoWLEi+vbtCwCoXLkyhg4dijlz5sDc3BxVqlTB7Nmz0aBBA7Rr106De0ZERESaotFwc+XKFXh4eEheBwQEICAgAAMHDkRoaCi+//57ZGVlwdfXF+np6XBxcUF4eDhMTEwky/z888/Q1dWFt7c3srOz0bZtW/z222/Q1dXVxC4RERGRhmk03LRp0wbp6ekK5wsEAvj7+8Pf319hGyMjIyxZsgRLlixRQ4VERERU3pTZMTdERERExcFwQ0RERFqF4YaIiIi0CsMNERERaRWGGyIiItIqDDdERESkVRhuiIiISKsw3BAREZFWYbghIiIircJwQ0RERFqF4YaIiIi0CsMNERERaRWGGyIiItIqGv1WcCIq256JBHjyJk/TZRRZVvkrmYhUiOGGiBR68iYPM86+0HQZRebfrJqmSyAiDeJlKSIiItIqDDdERESkVRhuiIiISKsw3BAREZFWYbghIiIircJwQ0RERFqF4YaIiIi0CsMNERERaRWGGyIiItIqDDdERESkVRhuiIiISKsw3BAREZFWYbghIiIircJwQ0RERFqF4YaIiIi0CsMNERERaRWlw010dDRevnypcH5KSgqio6NVUhQRERFRcSkdbjw8PBAZGalw/unTp+Hh4aGSooiIiIiKS+lwIxaLC5z/7t076OjwKhcRERFpll5BMzMzM5GRkSF5nZqaisePH8u0S09Px+7du1GjRg3VV0hERERUBAWGm19//RWLFy8GAAgEAvj7+8Pf319uW7FYjP/973+qr5CIiIioCAoMN+3atYORkRHEYjHmz5+P3r17o1GjRlJtBAIBKlSogK+++gpNmzZVa7FEREREhSkw3LRs2RItW7YEAIhEInh4eKBBgwalUhgRERFRcRQYbj7m5+enzjqIiIiIVEJhuNm2bRsAYMCAARAIBJLXhRk4cKBqKgOQm5uLgIAA7NixA8nJybC0tET//v3h5+cHPb0PpYvFYgQGBmLTpk1IT0+Hi4sLgoODUb9+fZXVQUREROWHwnAzfvx4CAQC9OnTBwYGBhg/fnyhKxMIBCoNN8uXL8e6desQGhoKJycnXL9+HT4+PjAwMMDMmTMBACtWrEBISAhCQkLg4OCAxYsXo1evXrh48SJMTExUVgsRERGVDwrDzdWrVwEABgYGUq9L04ULF9ClSxd07doVAGBra4uuXbvi8uXLAD6ctQkNDcWUKVPg6ekJAAgNDYWDgwN27doFb2/vUq+ZiIiINEthuLGxsSnwdWlo2bIl1q9fjzt37uCLL77ArVu3EBUVhalTpwIAHj58iOTkZLRv316yjLGxMVxdXREbG8twQ0RE9BlSekCxJkyZMgWvX79GixYtoKuri5ycHMyYMQOjRo0CACQnJwMAzM3NpZYzNzfHs2fPFK43ISGhWPUUdzmSj/2pOurqy0wDC4hE2WpZtzq9z3lforo1sc+ZrzKRkPK81Lerbvycqxb78wMHB4cC5xcp3Jw6dQqbNm3CgwcPkJaWJvOVDAKBAHFxcUUuUpHw8HD89ddfWLduHRwdHXHt2jX4+fnBxsYGw4YNk9rux8Riscy0jxXWKfIkJCQUazmSj/2pOursy4xUMQwNRWpZtzrp6+nD0NCoWMuKRNnFXrYkTE1M4VC1cqlvV534OVct9qfylA43oaGhmD17NqpVq4amTZuWyt1Ic+bMwcSJE9GnTx8AQIMGDfD48WP88ssvGDZsGCwtLQEAz58/h5WVlWS5ly9fypzNISIios+D0uEmJCQErVq1wu7duyWDjNXt7du30NXVlZqmq6uLvLw8AB8GGFtaWiIyMhJNmjQBAGRnZ+PcuXOYP39+qdRIREREZYvS4SYlJQXTp08vtWADAF26dMHy5ctha2sLR0dHxMfHIyQkBAMGDADw4XKUj48Pli5dCgcHB9jb2yM4OBgVK1ZE3759S61OIiIiKjuUDjfOzs549OiROmuRsXjxYixatAjTp0/Hy5cvYWlpieHDh0uecQMA33//PbKysuDr6yt5iF94eDifcUNE5YpARweXUvM0XUaR1aqogxqG4sIbEpUipcPNokWLMHDgQLi5uaFt27bqrEnCxMQEgYGBCAwMVNimsG8rJyIqD15m5SLg4ktNl1Fkwa3NUcNQ8Q0cRJqgdLgJCAiAqakpevbsCTs7O1hbW8uMhxEIBNixY4fKiyQiIiJSltLh5tatWxAIBLCysoJIJMLdu3dl2hR0+zURERFRaVA63Fy7dk2ddRARERGphI6mCyAiIiJSJaXP3Dx+/FipdtbW1sUuhoiIiKiklA43X375pVJjalJTU0tUEBEREVFJKB1uVq9eLRNucnNz8fDhQ/z111+wsLCQfKElERERkaYoHW4GDx6scN6UKVPQvn17vH79WiVFERERERWXSgYUV6pUCYMHD8avv/6qitURERERFZvK7pbS19fHs2fPVLU6IiIiomJRSbi5du0afvvtN9SrV08VqyMiIiIqthLfLZWRkYHMzExUqlQJISEhKi2OiIiIqKiUDjetWrWSCTcCgQBCoRB169ZFnz59IBQKVV0fERERUZEoHW5CQ0PVWQcRERGRSvDrF4iIiEirMNwQERGRVmG4ISIiIq3CcENERERaheGGiIiItIpS4SY7OxtBQUE4efKkuushIiIiKhGlwo2RkRF++eUXJCYmqrseIiIiohJR+rJUo0aNcO/ePXXWQkRERFRiSoebOXPmICwsDEeOHFFnPUREREQlovQTileuXAmhUIiBAweiZs2aqF27NoyNjaXaCAQC7NixQ+VFEhERESlL6XBz69YtCAQCWFlZAQAePXok00beF2sSERERlSalw821a9fUWQcRERGRSvA5N0RERKRVihRucnNzsWPHDkycOBFeXl74999/AQDp6enYs2cPkpKS1FIkERERkbKUDjcZGRn49ttvMXbsWOzbtw/Hjh1DSkoKAMDExASzZ8/G2rVr1VYoERERkTKUDjfz5s3DrVu3sHPnTsTFxUEsFkvm6erqwsPDA8eOHVNLkURERETKUjrcHDx4EGPGjEHHjh3l3hVlZ2eHx48fq7Q4IiIioqJSOtykp6ejTp06CueLxWK8e/dOJUURERERFZfS4cbGxgY3btxQOD86Ohr29vYqKYqIiIiouJQON/369UNYWBiio6Ml0/IvT61ZswYHDhzAoEGDVF8hERERUREo/RC/qVOn4tKlS+jRowfs7e0hEAjg5+eH1NRUJCcno3v37hg7dqw6ayUiIiIqlNLhRl9fHzt27MDOnTuxd+9eCAQC5OTkoHHjxujduzf69+/Pr18gIiIijVM63OTr168f+vXrp45aiIiIiEqsyOEGAP7991/Jbd/W1tZo0KABz9oQERFRmVCkcLN7927MnTsXT58+lTzETyAQoGbNmpg7dy7P6BAREZHGKX231JYtWzBq1ChUqFAB8+bNw9atW7FlyxbMmzcPxsbGGDt2LLZs2aLyApOSkjBu3DjY2dnB0tISLVq0wNmzZyXzxWIxAgIC4OjoiOrVq6N79+64efOmyusgIiKi8kHpMzfLli2Di4sLDhw4ACMjI6l5o0ePRrdu3bBs2TIMHjxYZcWlp6ejc+fOaNmyJXbs2AEzMzM8fPgQ5ubmkjYrVqxASEgIQkJC4ODggMWLF6NXr164ePEiTExMVFYLERERlQ9Kn7l58uQJ+vXrJxNsAMDIyAheXl54+vSpSotbuXIlqlevjjVr1sDFxQW1a9fGN998g3r16gH4cNYmNDQUU6ZMgaenJ5ycnBAaGorXr19j165dKq2FiIiIygelw42joyOePXumcP7Tp08loUNVDh48CBcXF3h7e8Pe3h6tW7fG2rVrJeN9Hj58iOTkZLRv316yjLGxMVxdXREbG6vSWoiIiKh8UPqy1Pz58zF8+HA0btwYvXr1kpq3e/duhIWFISwsTKXFPXjwAOvXr8f48eMxZcoUXLt2DbNmzQIAjBkzBsnJyQAgdZkq/3VBQSwhIaFY9RR3OZKP/ak66urLTAMLiETZalm3Or3PeV+iujWxzyWtWVMyX2UiIeW5wvn8nKsW+/MDBweHAucrHW5WrVoFMzMzjBw5En5+fqhTpw4EAgHu3buHFy9ewM7ODitXrsTKlSslywgEAuzYsaPYxefl5eGrr77C3LlzAQCNGzfGvXv3sG7dOowZM0ZqOx8Ti8UF3ppeWKfIk5CQUKzlSD72p+qosy8zUsUwNBSpZd3qpK+nD0ND2UvoyhCJsou9bEmUpGZNMjUxhUPVynLn8XOuWuxP5Skdbm7dugWBQAArKysAkIyvMTQ0hJWVFUQiEW7fvi21TEmffWNpaSlzqeuLL75AYmKiZD4APH/+XFIXALx8+VLmbA4RERF9HpQON9euXVNnHXK1bNkSd+/elZp29+5dWFtbAwBsbW1haWmJyMhINGnSBACQnZ2Nc+fOYf78+aVeLxEREWme0gOKNWH8+PG4ePEigoODce/ePezduxdr167FqFGjAHw4M+Tj44Ply5cjIiICN27cwPjx41GxYkX07dtXw9UTERGRJhTr6xdKS5MmTbBlyxbMnz8fS5YsgZWVFX744QdJuAGA77//HllZWfD19UV6ejpcXFwQHh7OZ9wQERF9psp0uAGAzp07o3PnzgrnCwQC+Pv7w9/fvxSrIiIiorKqTF+WIiIiIioqhhsiIiLSKgw3REREpFWUDjeNGzfGoUOHFM4/fPgwGjdurJKiiIiIiIpL6XDz6NEjvHnzRuH8N2/e4PHjxyopioiIiKi4inRZqqAnDt+9e5e3XxMREZHGFXgr+NatW7Ft2zbJ6+DgYGzatEmmXXp6Om7cuFHgLdtEREREpaHAcPPmzRvJN28DQEZGBvLy8qTaCAQCVKhQAcOHD4efn596qiQiIiJSUoHhZvTo0Rg9ejQA4Msvv0RgYCC6detWKoURERERFYfSTyiOj49XZx1EREREKlHkr1949eoVEhMTkZaWBrFYLDO/VatWKimMiIiIqDiUDjdpaWmYNWsW9uzZg9zcXJn5YrEYAoEAqampKi2QiIiIqCiUDjdTp07FgQMHMHr0aLRq1QpCoVCNZREREREVj9Lh5vjx4xg7diwWLVqkznqIiIiISkTph/gZGBjAzs5OnbUQERERlZjS4cbT0xPHjh1TZy1EREREJaZ0uJk0aRKSkpIwbtw4XLx4EUlJSXjx4oXMPyIiIiJNUnrMjYuLCwQCAeLi4rBjxw6F7Xi3FBEREWmS0uFm5syZBX5xJhEREVFZoHS48ff3V2cdRERERCqh9Jibj+Xm5iI1NRU5OTmqroeIiIioRIoUbv755x/07NkTNWvWhL29PaKjowEAKSkp6N+/P06fPq2WIomIiIiUpXS4uXDhArp164b79+9jwIABUt8rZWZmhtevX+PPP/9US5FEREREylI63CxYsAB2dnaIjY3FnDlzZOa3adMGly5dUmlxREREREWldLj5559/MGTIEBgZGcm9a6pWrVpITk5WaXFERERERaV0uNHR0YGOjuLmycnJMDY2VklRRERERMWldLhxdnbG4cOH5c579+4ddu7ciebNm6usMCIiIqLiUDrcTJs2DWfOnMHEiRNx7do1AEBSUhKOHz+OHj164P79+5g+fbraCiUiIiJShtIP8XNzc8OaNWvg6+uLrVu3AgB8fHwgFotRuXJlrFu3Ds2aNVNboURERETKUDrcAEDfvn3RrVs3REZG4r///kNeXh7q1KmDDh06oFKlSuqqkYiIiEhpRQo3AFChQgV0795dHbUQERERlZjSY24OHToEX19fhfN9fX0VDjgmIiIiKi1Kh5tVq1bh7du3CudnZ2djxYoVKimKiIiIqLiUDjc3btyAs7OzwvmNGzfGrVu3VFETERERUbEpHW5ycnKQlZWlcH5WVhZEIpFKiiIiIiIqLqXDjZOTEyIiIpCXlyczLy8vDxEREXB0dFRpcURERERFpXS4GTduHC5fvoyBAwciLi4OIpEIIpEIcXFxGDRoEC5fvoyxY8eqs1YiIiKiQil9K3ifPn1w//59BAQE4NixYwAAgUAAsVgMgUCAWbNmwcvLS22FEhERESmjSM+5mTFjBvr27Yv9+/fjwYMHEIvFqFOnDjw8PFC7dm01lUhERESkPKUuS2VlZcHDwwObN29G7dq1MWnSJCxduhTLli3DpEmTSi3YLF26FEKhUOp5O2KxGAEBAXB0dET16tXRvXt33Lx5s1TqISIiorJHqTM3xsbGuHr1Kvr27avuehS6ePEiNm3ahAYNGkhNX7FiBUJCQhASEgIHBwcsXrwYvXr1wsWLF2FiYqKhaolkPRMJ8OSN7ID8kso0sEBGqljl6wWALNWXS0SkdkpflmrdujViYmIwfPhwddYjV0ZGBkaPHo1Vq1Zh8eLFkulisRihoaGYMmUKPD09AQChoaFwcHDArl274O3tXeq1Einy5E0eZpx9ofL1ikTZMDRUz2MY/JtVU8t6iYjUSem7pYKCgvDPP//gf//7Hx48eCD3lnB1yQ8v33zzjdT0hw8fIjk5Ge3bt5dMMzY2hqurK2JjY0utPiIiIio7lD5z06xZM4jFYsklIB0dHejr60u1EQgEePr0qUoL3LRpE+7du4c1a9bIzEtOTgYAmJubS003NzfHs2fPFK4zISGhWLUUdzmS73Prz0wDC4hE2WpZt7rW+z7nvdrWrU4lrVsT+1xe+zrzVSYSUp4rnP+5fc7Vjf35gYODQ4HzlQ43vXr1gkAgKHFBRZGQkID58+fj77//hoGBgcJ2n9aVf3u6IoV1iqJairMcyfc59mdGqlgtl48+XJYyUvl6AUBfT19t61anktStzv4sSHnta1MTUzhUrSx33uf4OVcn9qfylA43oaGh6qxDrgsXLiAlJQVff/21ZFpubi5iYmKwYcMGnD9/HgDw/PlzWFlZSdq8fPlS5mwOERERfR6K9Jyb0ta9e3d89dVXUtMmTJgAOzs7TJs2Dfb29rC0tERkZCSaNGkC4MO3k587dw7z58/XRMlERESkYUoPKAaAR48eYfLkyXB2doa1tTXOnj0LAEhJScH06dMRFxen0uKEQiGcnJyk/lWoUAFVqlSBk5MTBAIBfHx8sHz5ckRERODGjRsYP348KlasqNHb1omIiEhzlD5zc/v2bXTp0gV5eXlo2rQpHj16hNzcXACAmZkZLl68CJFIhNWrV6utWHm+//57ZGVlwdfXF+np6XBxcUF4eDifcUNERPSZUjrczJ07FyYmJjh+/Dh0dXVhb28vNf/bb7/F3r17VV2fjIMHD0q9FggE8Pf3h7+/v9q3TURERGWf0pelYmJiMGrUKFhYWMi9E8na2rrA26+JiIiISoPS4SYnJwcVK1ZUOD8tLQ26uroqKYqIiIiouJQON05OToiKipI7TywWY//+/XB2dlZVXURERETFonS48fHxwb59+7B48WKkpqYCAPLy8nDnzh189913uHLlCiZNmqS2QomIiIiUofSA4j59+uDx48dYtGgRAgMDJdMAQFdXFwsXLkSnTp3UUyURERGRkor0EL8pU6agb9++iIiIwL1795CXl4c6deqgR48esLW1VVeNREREREorNNyIRCIcOnQIDx48QNWqVdG5c2eMHz++NGojIiIiKrICw01ycjK6deuG+/fvQywWAwAqVqyI7du3o1WrVqVSIBEREVFRFDigeOHChXjw4AHGjx+P7du3IyAgAIaGhpg5c2Zp1UdERERUJAWeuTl58iQGDhyIhQsXSqZZWFhg1KhRePLkCWrVqqX2AomIiIiKosAzN8nJyWjRooXUtJYtW0IsFiMxMVGthREREREVR4HhJjc3F0ZGRlLT8l9nZ2erryoiIiKiYir0bqkHDx7g8uXLkteZmZkAgISEBFSqVEmmvYuLiwrLIyIiIiqaQsNNQEAAAgICZKZ/OqhYLBZDIBBInl5MREREpAkFhpuQkJDSqoOIiIhIJQoMN4MGDSqtOoiIiIhUQukvziQiIiIqDxhuiIiISKsw3BAREZFWYbghIiIircJwQ0RERFqF4YaIiIi0CsMNERERaRWGGyIiItIqDDdERESkVRhuiIiISKsw3BAREZFWYbghIiIircJwQ0RERFqF4YaIiIi0CsMNERERaRWGGyIiItIqepougIiIyi+Bjg4upebJnZdpYIGMVHEpV6ScWhV1UMOwbNZGJcdwQ0RExfYyKxcBF1/KnScSZcPQUFTKFSknuLU5ahgKNF0GqQkvSxEREZFWYbghIiIircJwQ0RERFqF4YaIiIi0CsMNERERaZUyHW6WLVsGNzc3WFtbw87ODl5eXrhx44ZUG7FYjICAADg6OqJ69ero3r07bt68qaGKiYiISNPKdLg5e/YsRo4ciSNHjiAiIgJ6enro2bMn0tLSJG1WrFiBkJAQBAUF4eTJkzA3N0evXr3w6tUrDVZOREREmlKmn3MTHh4u9XrNmjWwsbHB+fPn0bVrV4jFYoSGhmLKlCnw9PQEAISGhsLBwQG7du2Ct7e3JsomIiIiDSrTZ24+9fr1a+Tl5UEoFAIAHj58iOTkZLRv317SxtjYGK6uroiNjdVQlURERKRJZfrMzaf8/PzQqFEjNG/eHACQnJwMADA3N5dqZ25ujmfPnilcT0JCQrG2X9zlSL7PrT8zDSwgEmWrZd3qWu/7nPdqW7c6lbRuTeyztvZ1Wd2nzFeZSEh5rukyiuxz+7mpiIODQ4Hzy024+eGHH3D+/HkcPnwYurq6UvMEAulHaIvFYplpHyusU+RJSEgo1nIk3+fYnxmpYrU8iv7DI+6NVL5eANDX01fbutWpJHWrsz8Loo19ram+VIapiSkcqlbWdBlF8jn+3CyucnFZyt/fH7t370ZERARq164tmW5paQkAeP5cOn2/fPlS5mwOERERfR7KfLiZNWsWdu3ahYiICHzxxRdS82xtbWFpaYnIyEjJtOzsbJw7dw4tWrQo7VKJiIioDCjTl6VmzJiB7du3Y/PmzRAKhZIxNhUrVkSlSpUgEAjg4+ODpUuXwsHBAfb29ggODkbFihXRt29fDVdPREREmlCmw826desAQHKbd75Zs2bB398fAPD9998jKysLvr6+SE9Ph4uLC8LDw2FiYlLq9RIREZHmlelwk56eXmgbgUAAf39/SdghIiKiz1uZH3NDREREVBQMN0RERKRVGG6IiIhIqzDcEBERkVZhuCEiIiKtwnBDREREWoXhhoiIiLQKww0RERFpFYYbIiIi0ioMN0RERKRVGG6IiIhIqzDcEBERkVZhuCEiIiKtUqa/FZxInmciAZ68ydN0GUWWVf5KJiIqlxhuqNx58iYPM86+0HQZRebfrJqmSyAi+izwshQRERFpFYYbIiIi0ioMN0RERKRVGG6IiIhIqzDcEBERkVZhuCEiIiKtwnBDREREWoXhhoiIiLQKww0RERFpFYYbIiIi0ioMN0RERKRVGG6IiIhIqzDcEBERkVZhuCEiIiKtoqfpAoiIiEqbQEcHl1LzNF1GkRhVqqbpEsoNhhsiIvrsvMzKRcDFl5ouo0jmNDbUdAnlBi9LERERkVZhuCEiIiKtwnBDREREWoXhhoiIiLQKww0RERFpFYYbIiIi0ioMN0RERKRVGG6IiIhIq2hNuFm3bh2+/PJLWFpa4ptvvkFMTIymSyIiIiIN0IpwEx4eDj8/P0yfPh1nzpxB8+bN0a9fPzx+/FjTpVEZI8jLQ4eLR/DDxjnocPEIBHnl6/HrRERUOK34+oWQkBAMGjQIw4cPBwAsWbIEJ06cwIYNGzB37lwNV1e2PRMJ8ORN6f+CzzSwQEaquFjLZhWzXEFeHlYuG4OmN2Nh9F4Ej7N7cKl+C0yethZiHa3I+USkxYyMjHGpmD83NalWRR3UMCzdust9uHn37h3i4uIwadIkqent27dHbGysyrbj4OCgsnWVJU/e5GHG2Rca2nrxtuvfrHhfHtf+8jFJsAEAo/ciuNy8ALfLx3CyWedirbMsMDQ00nQJWoX9qTrsS9V6laeHAI39vC6+4NbmqGEoKNVtCtLT08tfDPzIs2fPUL9+fRw8eBCtWrWSTA8KCsLOnTtx6dIlDVZHZYnR1Kkw3LhRZrrou++QvWyZBioiIiJ10Jpz8QKBdCoUi8Uy0+jzltOuHcRG0n9Jio2MkNOunWYKIiIitSj34cbMzAy6urp4/vy51PSXL1/C3NxcQ1VRWZTj4YGc1q0lAUdsZISc1q2R4+6u4cqIiEiVyn24MTAwgLOzMyIjI6WmR0ZGokWLFhqqisokHR283bEDb9euhei77/B27Vq83bED4GBiIiKtUu4HFAPAhAkTMHbsWLi4uKBFixbYsGEDkpKS4O3trenSqKzR0UFOjx7I6dFD05UQEZGaaMWfrL1790ZAQACWLFmCNm3a4Pz589ixYwdsbGyKtb4//vgD7u7usLGxgVAoxMOHDwtdZsuWLRAKhTL/srOzi1WDNilOfwLAvn370KJFC1hYWKBFixbYv3+/mistH0QiEXx9fVG3bl3UrFkTAwYMwJMnTwpchsfnB0V92Of169fRrVs3VK9eHfXr10dQUBDE4nJ9D4ZKFaU/Hz58KPcYPH78eClWXDZFR0djwIABqF+/PoRCIbZs2VLoMjw2C6YV4QYARo0ahWvXruH58+c4ffq01J1TRfX27Vu0b98efn5+RVquQoUKuH37ttQ/IyPeClmc/rxw4QK+++479OvXD1FRUejXrx9GjBjBu98A+Pv7Y//+/Vi/fj0OHTqEV69ewcvLC7m5uQUu97kfn0V92GdmZiZ69eoFCwsLnDx5EoGBgVi1ahVWr15dypWXTcV9eOru3buljsG2bduWUsVl15s3b+Dk5ITAwEAYGxsX2p7HZuHK/a3g6nTlyhW4ubnh6tWrsLW1LbDtli1bMHPmzEL/gv6cFaU/vb29kZaWhr1790qmeXp6olq1ali/fr2aKy27MjIyYG9vj5CQEPTv3x8AkJiYiEaNGmHXrl3o0KGD3OV4fAIdOnRAgwYNsHLlSsm0Jk2awNPTU+7DPtevX4+ffvoJd+7ckfzCWbJkCTZs2IAbN2589ndjFrU/Hz58iMaNGyMyMhJfffVVaZZartSqVQuLFy/G4MGDFbbhsVk4rTlzUxZkZWWhYcOGcHJygpeXF65evarpksqtixcvon379lLTOnTooNIHM5ZHcXFxeP/+vVTfWFlZoV69eoX2zed8fOY/7PPTY6qgh31euHABX3/9tdRf0h06dMCzZ8+UvrSqrYrTn/mGDh0Ke3t7dO7cGfv27VNnmVqLx2bhGG5UxMHBAatXr8bWrVuxbt06GBoaokuXLvjvv/80XVq5lJycLHMrv7m5ucwt/5+b58+fQ1dXF2ZmZlLTC+ubz/34TElJQW5ubpGOqefPn8ttnz/vc1ac/qxUqRIWLFiAjRs3YufOnWjbti28vb2xffv20ihZq/DYLJxW3C2ljIULFyI4OLjANvv370ebNm2Ktf7mzZujefPmktctWrRAmzZtsGbNGixevLhY6yzL1N2fwOf1YEZl+1ORwvrmczs+FSnqMSWvvbzpn6ui9KeZmZnU1+R89dVXSE1NxYoVK+Dl5aXWOrURj82CfTbhxsfHRzJGQRErKyuVbU9XVxfOzs64d++eytZZlqi7Py0tLT+rBzMq258XL15Ebm4uUlJSUK3a/3/H1suXL+Hq6qr09rT9+PxUcR72aWFhIbc9AK09DpWlqoenuri4KHVnEEnjsVm4zybcmJmZyZzKVyexWIzr16+jYcOGpbbN0qTu/mzWrBkiIyMxefJkyTRtfjCjsv3p7OwMfX19REZGol+/fgCAJ0+e4Pbt20XqG20/Pj/18cM+e/bsKZkeGRmJHgqeedS8eXP89NNPyM7OltxVFhkZiRo1ahQ6IF7bFac/5bl27RosLS3VUKF247FZOI65kSM5ORnx8fG4e/cuAOD27duIj49HWlqapE2PHj0wb948yevAwECcOHECDx48QHx8PCZOnIjr16/ju+++K/X6y5ri9Oe4ceNw5swZLFu2DHfu3MGyZcsQFRUFHx+fUq+/LKlcuTKGDh2KOXPm4NSpU7h69SrGjh2LBg0aoN1H35HF41PWhAkTsHXrVoSFheH27duYNWuW1MM+582bJ/WLuW/fvjA2Nsb48eNx48YNREREYPny5Rg/fjxP/aPo/bl161bs3LkTt2/fRkJCAlatWoV169ZhzJgxmtqFMuP169eIj49HfHw88vLykJiYiPj4eMlt9Tw2i+6zOXNTFBs2bEBQUJDkdf7lgpCQEMnteffv30etWrUkbTIyMvD999/j+fPnMDU1xZdffolDhw7BxcWldIsvg4rTn/lPml64cCECAgJQp04dbNiwAU2bNi3d4sugn3/+Gbq6uvD29kZ2djbatm2L3377Dbq6upI2PD5l9e7dG6mpqViyZAmSk5NRv359qYd9JiUl4f79+5L2lStXxp49ezBjxgy4ublBKBRiwoQJmDhxoqZ2oUwpan8CQHBwMB4/fgxdXV3Y2dlh9erVHG+DD4/J8PDwkLwOCAhAQEAABg4ciNDQUB6bxcDn3BAREZFW4WUpIiIi0ioMN0RERKRVGG6IiIhIqzDcEBERkVZhuCEiIiKtwnBDREREWoXhhrRWQEAAhEKhzPRff/0VX331FczMzNCoUaNCpxMRUfnCcEPlwpYtWyAUCiX/LC0t4ejoiN69e+O3337Dq1evlFrP2bNn8cMPP6Bx48ZYtWoVAgICCpz+OWvUqJFUn1tYWMDZ2Rm+vr5ITk4u1jpfv36NgIAAREVFqbjasunkyZMYMmQIHB0dYW5uDhsbG3To0AGLFi3C06dPNV1esUVFRUmOi82bN8ttM2LECMlnlai08QnFVK74+fmhTp06eP/+PZ4/f46zZ8/C398fISEh2LZtm9R3Jfn6+mLq1KlSy+f/Ul2+fLnUWR1F0z93DRo0kHy/l0gkwrVr17Bp0yacOXMGMTExUk9FVsabN28kT6suyTfGl3VisRjTp0/Hhg0b4OjoiGHDhsHa2hpZWVmIi4vDmjVrsHHjRslXkpRXRkZG2LlzJ4YMGSI1PTMzE4cPH4aRkZHk26qJShPDDZUrHTp0QLNmzSSvp02bhtOnT2PAgAEYOHAgLly4AGNjYwCAnp4e9PSkD/H8b879NMAoml4Sb9++RYUKFVS2Pk2oXr26zOPxK1eujODgYFy7dg3Ozs6aKayMCwkJwYYNGzBmzBgEBgZCR0f6JHlgYCBWrlxZ6HqysrIkx3NZ9O233+LAgQN49uwZatSoIZkeERGBvLw8tG/fHpGRkRqskD5XvCxF5d4333wDX19fPH78GDt27JBM/3TMjVAoxPr16yX/LxQKJW3kTc8XGRkJd3d3WFlZoWbNmnB3d0dsbKxUDfnruXXrFsaNG4c6deqgZcuWxVrHf//9h6lTp6JOnTqoVasWhg8fjtTUVJn9joyMhIeHB6ytrWFlZYVvvvkGYWFhUm2uXLkCLy8v2NjYoHr16mjfvj0OHz5cxB6Wln+ZQV9fX2p6UlISvv/+ezg6OsLCwgJNmjTBihUrJH+5P3z4EPXq1QMABAUFSfrax8cH//77L4RCIfbs2SNZ3/379yEUCmW+uXzq1Kn44osvirWfmZmZ+PHHH9GoUSNYWFigYcOG+OmnnyASiaTaCYVCTJ06FceOHUObNm1gaWmJJk2aYNeuXYX2T1ZWFpYuXYp69eohICBAJtgAgKmpKX788UepaY0aNUKfPn1w5swZdOzYEZaWlli+fDkAIC0tDdOmTUO9evVgYWGB5s2bY/Xq1VJnRR4+fAihUIgtW7bIbK9Ro0ZSXzqbf5n3zJkz8PX1Rd26dVGrVi0MGzYMSUlJhe5jvs6dO6NSpUoy/bJz50507NgRVapUkbucMp+HR48eYfr06WjWrBlq1KgBGxsbeHl54ebNm1Lt8i+R7dq1C6tXr0ajRo1gaWmJTp064erVq1Jtnz9/jkmTJqFBgwawsLCAo6MjvLy8cP36daX3mcoHhhvSCvlnF06ePKmwzZo1a9C2bVvJ/69ZswYeHh4KpwPArl270KdPH+jq6mL27NmYPXs2UlNT0aNHD1y6dElmG97e3khLS8Ps2bMxbty4Yq1j5MiRePr0KWbPno1hw4bhwIEDmDlzplSbv/76C71790ZSUhImTZqEefPmwcXFBUeOHJG0OXv2LLp06YLnz5/D19cX8+bNg4GBAQYOHIiIiAil+vX9+/dISUlBSkoKnj59iqNHj2LFihVwdnaGk5OTpN2LFy/QsWNHHDlyBMOHD0dQUBCaNm2KuXPnwt/fHwBQrVo1LFmyBADg7u4u6Wtvb280aNAAQqEQ0dHRknVGR0dDR0cHiYmJePjwoWR6TEwMvv766yLvZ1ZWFtzd3fHnn3+id+/eWLx4MTp37ozVq1dLvsn6YxcvXsSECRPQrVs3LFiwABUqVMCYMWNw+/btAvssNjYWaWlp6Nu3b5Ev2927dw/Dhg2Dq6srgoKC0KxZM4hEInh4eGDTpk3o0aMHFi1aBFtbW/z444/44YcfirT+T/n5+SEuLg4zZ87EiBEj8Pfff6N379549+6dUssbGRnBw8MDO3fulExLSkpCVFSU5AtyP6Xs5+HKlSuIjo6Gh4cHAgIC4OPjgytXrqBbt25yx3ytXr0a27Ztw5gxY+Dn54e7d+9i8ODBeP/+vaTN8OHDsW/fPgwcOBDBwcEYO3Ys8vLyyv3lQZLFy1KkFWrVqgVTU1OZbyH+mJeXF86fP48zZ85IXWpp2LCh3Olv3rzBjBkz4OXlhdDQUMl0b29vtGzZEvPnz5cJCfb29vjzzz9LtI4vvvgCa9eulbwWi8X4/fffsXTpUlSuXBmZmZmYOXMmGjRogCNHjqBixYpSbfP/O3XqVDRv3hz79u2TnD0YPXo0OnfujDlz5qBHjx4FdyqAM2fOwM7OTmqai4sL/vrrLwgEAsm0hQsXQiQSITo6GhYWFpJ9rF69OlavXg0fHx/Y2tqiR48e8PX1RYMGDWQud7Vs2RIxMTGS1+fOnUP79u0RGxuLc+fOwdbWFikpKbhz5w6+++67Iu/nr7/+ioSEBJw6dUpyBgkA6tevjxkzZiAmJgaurq6S6bdu3UJ0dLSkbc+ePdGwYUNs3rwZCxYsUNhnt27dkqz3Y2KxWOYMnKmpqdQZsPv372Pr1q3o1q2bZNratWvx77//YuXKlRg2bBgAYNSoURg6dCh+++03jBo1SuY9KooDBw7A0NAQAODo6IhJkyZh69atGDFihFLL9+/fH1u2bMGtW7fg6OiInTt3olKlSujSpYtU2AaK9nno1KkTPD09pZb38vLC119/jT///BMzZsyQmpeZmYmYmBgYGRkBABwcHDBkyBCcPHkSnTt3RkZGBs6dO4cFCxZg0qRJkuU+HZdH2oFnbkhrVKpUCa9fv1bZ+iIjI5Geno7+/ftLzl6kpKQgKysL7dq1w7lz56T+KgQ+nHVR9TpatWqF3NxcJCYmStaZmZmJ6dOnSwUbAJLAce3aNSQkJKB///5IS0uTbDctLQ0dO3bEgwcP8OjRo0L74KuvvsLevXuxd+9ebN++HXPmzMF///2HIUOGICsrC8CHX9r79u1D586doaurK7WfHTp0QF5entQZGUVcXV1x8+ZNpKWlAfhwhqZNmzZo3ry5JPRER0dDLBZLztwUZT/37NmDFi1aoFq1alI1tmvXDsCHIPexNm3aSIUgCwsLODg44MGDBwXuR/6deyYmJlLTU1NTYWdnJ/Xv/PnzUm1q1aolFWwA4MiRIzAzM8PgwYMl0wQCASZPngyxWIyjR48WWE9BvL29JcEGAAYOHIjKlSsXaZ1t2rRBjRo1JGdvdu7cCXd3d0nI+FhRPg8fj1d7+/YtUlNTUblyZdjZ2SEuLk5m3YMHD5baZuvWrQFA8n4ZGRlBX18fZ8+elRxjpL145oa0xuvXr1GtWjWVre+///4DAPTq1Uthm4yMDKlt1q5du8TrsLa2lpqfP24o/wdy/tmpjy8LKap90qRJUn+lfuzly5ewsbFRuA4AqFq1quSXP/BhjIWDgwOGDh2KsLAwjB07Fi9fvkR6ejo2b96s8Lbg/AHbBXF1dYVYLEZMTAxcXFxw//59uLq6IicnB3/99ReAD4HH1NRUMg6nKPv533//4d9//1V4luPTGj99H4AP70VhvxjzQ82njycwNTXF3r17AXy4lBYcHCyzrK2trcy0R48ewc7OTuYSV37wUiakKvJpX+jp6cHW1haPHz9Weh06Ojro3bs3du7ciX79+iE+Ph7z58+X27Yon4fs7Gz8/PPP2LFjh8w4IDMzM5nlCvvcGBoaYu7cuZg7dy4cHBzQtGlTdOrUCf3795f7XlP5xnBDWuHJkyfIzMxE3bp1VbbOvLw8AB8uZ9SsWVNuG1NTU6nXn97ZUpx1KBqn8fElJwBSl4UU1f7TTz8pvKPJ3t5e4fIFyR+fFBMTIxmzAAB9+/aVuSU4nzLvi7OzMypWrIiYmBiIRCJUqFABzs7OeP/+PRYsWIAXL15IxtvkX34qyn7m5eWhbdu2mDZtmtx2n74/hb0Pijg6OgIAbty4AXd3d8l0fX19SVBMSUmRu2xJ7oxS5nhQZpni3Lrdr18/hISEwNfXF9WrV1d4m39RPg9+fn4ICwvDmDFj0LJlS5iamkJHRwf+/v5y90eZ92vixIlwd3fHoUOHcOrUKSxZsgTLli3D1q1b8c033xRpn6lsY7ghrbB9+3YAQPv27VW2zjp16gD4MBD247MXpb2OT+UHhRs3bsjcNfTpditVqqSy7ebLyckB8GH8BPBh30xNTZGTk1Potgr6Baynp4emTZsiJiYG7969Q7NmzaCvrw8XFxcYGhri8OHDuH79Onr37i1Zpij7WadOHbx+/Vrl/fGpFi1aoEqVKti1axemT59e5EHFn7KxscHVq1eRm5srta47d+5I5gOQ3JmUkZEhtbxIJFJ4B9Tdu3fh5uYmeZ2Tk4NHjx6hVatWRarR2dkZX3zxBaKiojB+/HiF+1yUz0N4eDgGDBiAwMBAqenp6emoWrVqker7WO3atTF+/HiMHz8eiYmJaNu2LX755ReGGy3DMTdU7p0+fRpLliyBra2twjs0iqNDhw6SZ7p8eqswoNylFlWs41Nubm4wNTXFsmXL8PbtW6l5+X+lOjs7w87ODqtWrZL5ZVfc7ebLH4+Rf2lIV1cXPXr0wIEDB+SOhcjIyJAZR5Geni533a6uroiPj8fx48clg3sNDQ3RpEkTrFy5Erm5uVKDfouyn71798Y///yDQ4cOybTLyspS2XgtY2NjTJ8+HXfu3FF4lqEoZ0c6d+6Mly9fYtu2bVLLr1q1CgKBAN9++y2AD5fDqlWrJvP05w0bNiA3N1fuujdu3Ch1XG7btg0ZGRno1KmT0vXlW7RoEWbNmiUzZuxjRfk86OrqyvTTrl278OzZsyLXBnwYt5M/TiyflZUVzM3NFR6PVH7xzA2VKydOnMC9e/eQk5ODFy9e4MyZM4iMjIS1tTW2bdsmdxBjcZmYmGDFihUYOXIkWrdujX79+sHS0hJPnjxBVFQUKlasWOhzT1Sxjk+ZmpoiICAAEydOhJubG/r164eqVavi5s2bePbsGTZv3gwdHR2sXr0affr0QcuWLTF48GDY2NggKSkJFy9exOPHj2UGs8qTlJQkOSv27t07XL9+HX/88QfMzMwwZswYSbuffvoJ0dHR6NKlC4YOHQonJye8evUKN27cwP79+/HPP//A0tISlSpVgoODA8LDw2Fvb4+qVavC1tYWTZs2BQB8/fXXyM3NlYy3ydeqVSsEBwfD2NhY6vJTUfZz0qRJOHr0KIYOHYr+/fvDxcUFIpEId+/exZ49e7Bz506pB0SWxIQJE/Dff/9h7dq1OHPmDHr06AFra2u8efMGt2/fxu7du1GxYkW5Y0c+NWzYMISFhWHKlCm4du0a7O3tcezYMRw9ehTjxo2TGjczYsQIBAcHY/z48WjWrBmuXLmC06dPF7gdDw8P9OnTB48ePcLatWvh6OiIQYMGFXmfO3XqVGgoKsrnoWvXrvjrr79gYmICJycnXLt2DeHh4TLj2pR19+5d9OjRAz179oSjoyMMDQ1x9OhR3L59u8C736h8YrihciX/FLWBgQGqVKkCJycnBAQEYPDgwTJ3p6hCz549UaNGDSxbtgy//vorsrKyYGlpiaZNm0puyy2NdXxq8ODBMDc3xy+//IJly5ZBV1cXdnZ2GDVqlKTN119/jRMnTmDx4sX4448/kJmZCXNzczRs2FDy7JnCXL9+HWPHjgXwIUiYmZnB3d0ds2fPlhozUa1aNZw4cQJLlizBwYMH8ccff6By5cqwt7eHn5+f1MPcQkJC4O/vjx9//BEikQgDBw6UhJtmzZrBwMAAACTT8vclf1r+/KLup7GxMSIiIrBixQqEh4dLAkbt2rXh4+MDBwcHpfpEGQKBAL/88gu6d++ODRs2YNOmTUhJSUGFChVgb2+PcePGYcSIEQrHnXzMyMgIERERWLBgAfbs2YO0tDTY2tpiwYIFmDhxolTbGTNmIDU1FeHh4di7dy9at26Nffv2SZ7b9KnAwEBEREQgKCgIIpEInTt3xpIlS6TuoFI1ZT8PgYGB0NfXx549e7B582Y4Oztj9+7d+N///les7VpZWaFfv344c+YMdu3aBYFAIDnrN3ToUFXtHpURgvT0dH7xBxHRZ2TLli2YMGECjh07prKzVURlCcfcEBERkVZhuCEiIiKtwnBDREREWoVjboiIiEir8MwNERERaRWGGyIiItIqDDdERESkVRhuiIiISKsw3BAREZFWYbghIiIirfJ/O2fa1qTKfUcAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"\n",
"mean_differences = pd.DataFrame({'Difference Between Group Means':age_differences})\n",
"\n",
"unit = ''\n",
"\n",
"print('Observed Difference:', observed_age_difference)\n",
"\n",
"fig, ax = plt.subplots(figsize=(8,5))\n",
"\n",
"ax.hist(mean_differences, density=True, alpha=0.8, ec='white', zorder=5)\n",
"\n",
"ax.scatter(observed_age_difference, 0, color='red', s=30, zorder=10).set_clip_on(False)\n",
"\n",
"y_vals = ax.get_yticks()\n",
"\n",
"y_label = 'Percent per ' + (unit if unit else 'unit')\n",
"\n",
"x_label = 'Difference Between Group Means'\n",
"\n",
"ax.set_yticklabels(['{:g}'.format(x * 100) for x in y_vals])\n",
"\n",
"plt.ylabel(y_label)\n",
"\n",
"plt.xlabel(x_label)\n",
"\n",
"plt.title('Prediction Under the Null Hypothesis');\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The empirical P-value of the test is the proportion of simulated differences that were equal to or less than the observed difference. This is because low values of the difference favor the alternative hypothesis that the smokers were younger on average."
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.008"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"empirical_P = np.count_nonzero(age_differences <= observed_age_difference) / 5000\n",
"empirical_P"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The empirical P-value is around 1% and therefore the result is statistically significant. The test supports the hypothesis that the smokers were younger on average."
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 1
}