{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "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": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Birth Weight</th>\n",
       "      <th>Gestational Days</th>\n",
       "      <th>Maternal Age</th>\n",
       "      <th>Maternal Height</th>\n",
       "      <th>Maternal Pregnancy Weight</th>\n",
       "      <th>Maternal Smoker</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>120</td>\n",
       "      <td>284</td>\n",
       "      <td>27</td>\n",
       "      <td>62</td>\n",
       "      <td>100</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>113</td>\n",
       "      <td>282</td>\n",
       "      <td>33</td>\n",
       "      <td>64</td>\n",
       "      <td>135</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>128</td>\n",
       "      <td>279</td>\n",
       "      <td>28</td>\n",
       "      <td>64</td>\n",
       "      <td>115</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>108</td>\n",
       "      <td>282</td>\n",
       "      <td>23</td>\n",
       "      <td>67</td>\n",
       "      <td>125</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>136</td>\n",
       "      <td>286</td>\n",
       "      <td>25</td>\n",
       "      <td>62</td>\n",
       "      <td>93</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1169</th>\n",
       "      <td>113</td>\n",
       "      <td>275</td>\n",
       "      <td>27</td>\n",
       "      <td>60</td>\n",
       "      <td>100</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1170</th>\n",
       "      <td>128</td>\n",
       "      <td>265</td>\n",
       "      <td>24</td>\n",
       "      <td>67</td>\n",
       "      <td>120</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1171</th>\n",
       "      <td>130</td>\n",
       "      <td>291</td>\n",
       "      <td>30</td>\n",
       "      <td>65</td>\n",
       "      <td>150</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1172</th>\n",
       "      <td>125</td>\n",
       "      <td>281</td>\n",
       "      <td>21</td>\n",
       "      <td>65</td>\n",
       "      <td>110</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1173</th>\n",
       "      <td>117</td>\n",
       "      <td>297</td>\n",
       "      <td>38</td>\n",
       "      <td>65</td>\n",
       "      <td>129</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1174 rows × 6 columns</p>\n",
       "</div>"
      ],
      "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": 2,
     "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": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "smoking_and_birthweight = births[['Maternal Smoker', 'Birth Weight']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Maternal Smoker</th>\n",
       "      <th>count</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>False</td>\n",
       "      <td>715</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>True</td>\n",
       "      <td>459</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Maternal Smoker  count\n",
       "0            False    715\n",
       "1             True    459"
      ]
     },
     "execution_count": 4,
     "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": 5,
   "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": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA60AAAFDCAYAAADRdOgEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABA2UlEQVR4nO3deVyU9f7//+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": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "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": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Maternal Smoker</th>\n",
       "      <th>Birth Weight</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>False</td>\n",
       "      <td>123.085315</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>True</td>\n",
       "      <td>113.819172</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Maternal Smoker  Birth Weight\n",
       "0            False    123.085315\n",
       "1             True    113.819172"
      ]
     },
     "execution_count": 7,
     "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": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-9.266142572024918"
      ]
     },
     "execution_count": 8,
     "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": 9,
   "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": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-9.266142572024918"
      ]
     },
     "execution_count": 10,
     "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": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Maternal Smoker</th>\n",
       "      <th>Birth Weight</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>False</td>\n",
       "      <td>120</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>False</td>\n",
       "      <td>113</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>True</td>\n",
       "      <td>128</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>True</td>\n",
       "      <td>108</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>False</td>\n",
       "      <td>136</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1169</th>\n",
       "      <td>False</td>\n",
       "      <td>113</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1170</th>\n",
       "      <td>False</td>\n",
       "      <td>128</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1171</th>\n",
       "      <td>True</td>\n",
       "      <td>130</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1172</th>\n",
       "      <td>False</td>\n",
       "      <td>125</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1173</th>\n",
       "      <td>False</td>\n",
       "      <td>117</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1174 rows × 2 columns</p>\n",
       "</div>"
      ],
      "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": 11,
     "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": 12,
   "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": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Maternal Smoker</th>\n",
       "      <th>Birth Weight</th>\n",
       "      <th>Shuffled Label</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>False</td>\n",
       "      <td>120</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>False</td>\n",
       "      <td>113</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>True</td>\n",
       "      <td>128</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>True</td>\n",
       "      <td>108</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>False</td>\n",
       "      <td>136</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1169</th>\n",
       "      <td>False</td>\n",
       "      <td>113</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1170</th>\n",
       "      <td>False</td>\n",
       "      <td>128</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1171</th>\n",
       "      <td>True</td>\n",
       "      <td>130</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1172</th>\n",
       "      <td>False</td>\n",
       "      <td>125</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1173</th>\n",
       "      <td>False</td>\n",
       "      <td>117</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1174 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      Maternal Smoker  Birth Weight  Shuffled Label\n",
       "0               False           120           False\n",
       "1               False           113            True\n",
       "2                True           128            True\n",
       "3                True           108           False\n",
       "4               False           136           False\n",
       "...               ...           ...             ...\n",
       "1169            False           113            True\n",
       "1170            False           128           False\n",
       "1171             True           130           False\n",
       "1172            False           125           False\n",
       "1173            False           117            True\n",
       "\n",
       "[1174 rows x 3 columns]"
      ]
     },
     "execution_count": 13,
     "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": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Shuffled Label</th>\n",
       "      <th>Birth Weight</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>False</td>\n",
       "      <td>118.955245</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>True</td>\n",
       "      <td>120.252723</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Shuffled Label  Birth Weight\n",
       "0           False    118.955245\n",
       "1            True    120.252723"
      ]
     },
     "execution_count": 14,
     "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": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2974785563020959"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "difference_of_means(original_and_shuffled, 'Birth Weight', 'Shuffled Label')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-9.266142572024918"
      ]
     },
     "execution_count": 16,
     "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": 17,
   "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": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.4498499322028806"
      ]
     },
     "execution_count": 18,
     "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": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2.32985359, -0.41244725,  1.42625958, ...,  0.49617441,\n",
       "       -1.12432012, -0.10480369])"
      ]
     },
     "execution_count": 19,
     "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": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observed Difference: -9.266142572024918\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqwAAAFuCAYAAABECkoSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABO70lEQVR4nO3dd1gU1/s28HtpCyK4agABKQoYBI1EsQTU2CtFsWCL0VixxQJRook9oKKJBQ2J5QtRsaJi7xWEaKLRaIwYC4qIShXFRWHfP3zZn+sCLrKwA9yf6+JKdubMzDPs7Hpz5syMKCMjQwYiIiIiIoHS0nQBRERERETFYWAlIiIiIkFjYCUiIiIiQWNgJSIiIiJBY2AlIiIiIkFjYCUiIiIiQWNgJSonfn5+kEgkOHv2rMJ0iUSCnj17ltl2g4KCCt0uvVHw+9m0aZOmSylSz549IZFIcO/ePU2XUi4K+0zwOC7epk2bIJFIEBQUpOlSAPD9IvVjYKVKRSKRKPzUqlULNjY26NatGzZs2IC8vDxNl6h2FSFwvaug5uL+cT179myZh3mhEHogbdy4MSQSCSwtLZGcnFxom6+++kpQAaWyHWNCC6RE5U1H0wUQlYXp06cDAPLy8nDnzh3s27cPcXFxOHXqFMLDwzVcnaLff/8dBgYGZbb+0aNHo0+fPqhbt26ZbYOqhufPn2PBggUIDQ3VdCkkcPzeIXVjYKVKKTAwUOH1tWvX0KlTJ+zZswexsbFwc3PTUGXKGjRoUKbrr127NmrXrl2m26Cqwc7ODpGRkRgzZgw++eQTTZdDAsbvHVI3DgmgKsHZ2Rnu7u4AgD/++APA/50O9PPzw40bNzBkyBDUr18fEokEV65ckS+7Z88eeHt7w9bWFqampmjatCnmzJmDrKysQrd16tQpdO/eHRYWFrC1tcWgQYPw77//FllbUack8/LyEBERge7du8PGxgZmZmb45JNPMHLkSFy6dAnAm1PJixYtAgCMHz9eYThEwenl4saSnTlzBv369UO9evVgamqKJk2aYPr06Xjy5IlS27fH4O7ZswcdOnSAubk5bG1tMXz4cCQlJRW5j+r0oXVcvnxZ3uNjZWUFb29vxMfHF7utlJQUzJgxA02bNoWZmRlsbGzQu3dvnD59Wqnt26ds4+Pj4ePjAxsbG0gkEmRkZBS5DYlEgpiYGABAkyZN5O9f48aNC22/YcMGuLm5wczMDA4ODpg0aVKR6y9J/aqYM2cO8vPzMWvWLJWXKe6UuxDHOf7666+QSCQIDg4udH5WVhYsLCzg7OwsH2L09rCcQ4cOoXPnzvLP/7Bhw3Dnzp1C15WSkoKAgAA0adIEpqamqFevHvr3749z584ptPPz88P48eMBAIsWLVL4nBf2u7ty5Qr69+8Pa2trmJubo3v37oiLiyu0hvz8fERERKBr166wtraGmZkZPvvsMyxbtgy5ublK7c+ePQtfX184OzvD1NQU9vb2aNeuHWbOnAmZ7P+e9F7Ue6vq8kTvYg8rVXl37txBly5d8PHHH2PAgAHIzMxEtWrVAADTpk3DunXrYGlpCQ8PD0gkEly8eBE//fQTjhw5gsOHD8PIyEi+rj179mD48OHQ1dVFr169YGFhgbi4OHTu3BmNGjVSuabc3FwMGjQIx44dQ506ddC7d2/UrFkTDx48wNmzZ2FnZ4dPP/0UgwYNAgDExMSgR48eCiGnRo0axW5jw4YNmDp1KgwMDODt7Y06deogPj4eYWFh2L9/Pw4ePAgrKyul5datW4eDBw+iR48ecHd3x8WLF7Fr1y5cvXoVMTExEIvFKu9naZSkjvj4ePTq1QtSqRSenp6ws7PDtWvX4OnpibZt2xa6/mvXrqF379548uQJOnTogB49eiAtLQ379+9Hr169sGLFCnzxxRdKy/3+++9YtmwZ3NzcMHToUCQnJ0NbW7vI/Zg+fTo2b96M+/fvY+zYsfL3rbD3b/bs2Thx4gS6deuG9u3b4+zZs4iIiMCtW7dw4MABtdRfnM6dO6N9+/Y4efIkDh06hG7dupVo+YpgwIABmDdvHn777TcEBAQovXdbtmzBixcvMGnSJKV5e/fuxbFjx+Dp6Yk2bdrgypUr2L17N86ePYsjR47Azs5O3vbevXvo3r07Hj58CHd3d/j4+ODRo0fYvXs3jh07hp9++glDhw4F8OYP08zMTBw4cADu7u5o3bq1fD3W1tYKNVy+fBkrVqxAy5YtMXToUDx48ADR0dHw9vbGmTNn8PHHH8vbvn79GkOGDMGhQ4dgb2+PPn36QCwWIyYmBvPmzcPp06exc+dO6Oi8iQpHjhyBr68vjIyM0L17d1haWiIjIwP//fcfwsLCMHfuXHnbwpR2earaeGRQlfDPP//Ie7GaNm2qMC8uLg5Tp07F999/rzB969atWLduHTw8PPDrr78qjDNdsmQJFi5ciKCgIPzwww8AgOzsbEyePBkikQj79++Hq6urvP13332HlStXqlzvokWLcOzYMbRr1w6bN2+WB2jgTc9rQQ/o4MGDkZiYiJiYGPTs2RODBw9Waf2JiYmYPn06qlWrhmPHjqFhw4byeQsWLEBISAimTZuGbdu2KS174sQJnD59Go6OjvJpI0eOxI4dO7B//374+PiovJ+loWodMpkMEyZMQE5ODsLDw+Ht7S1v/+uvvyIgIEBp3Xl5efjyyy+RmZmJvXv3KgSER48eoWPHjggICEDXrl1hamqqsOzJkyfx008/YdiwYSrtR2BgIM6dO4f79+/Dz88PNjY2Rbb9448/cP78eVhaWgJ4Ezg8PT0RGxuLixcvyo+50tT/PvPnz0fbtm3x/fffo1OnToIPGOfOnSvyQqXExESlaUZGRvD19cXatWtx6NAhpd7h//3vf9DR0ZGHybcdOnQIW7duRdeuXeXTVq5cie+++w4BAQGIioqST58yZQoePnyIGTNmYMaMGfLpEyZMQKdOnRAQEIAOHTqgbt268PDwkAfW1q1bKw15etvhw4cRFhYGX19f+bQNGzZgypQpCAsLw7Jly+TTf/zxRxw6dAijRo1CcHCwPIDn5+djypQpCA8Px9q1azF27FgAQEREBGQyGfbu3YsmTZoobDctLe29x0Jpl6eqjUMCqFIKCgpCUFAQFixYgFGjRqF9+/bIycmBh4eHfGhAAVNTU/lFWm9bvXo1tLW1sXLlSqWLoqZOnYratWsrBLoDBw4gPT0dPj4+CmEVAL755hsYGxurVHteXh7Wrl0LsViM5cuXK4RVANDW1kadOnVUWldRtm3bhtzcXIwYMUIhrAJAQEAAzM3NceTIETx8+FBp2TFjxiiERAD48ssvAQB//vlnqeoqCVXriI+PR0JCAlq2bKkQVgFgxIgRqF+/vtK6jxw5glu3bmHEiBEKYQ8A6tSpg4kTJ+Lly5fYs2eP0rKNGjVSOayW1DfffCMPqwCgo6ODIUOGAFDc59LU/z6NGjXCkCFDcPPmTWzYsOED96T8xMTEYNGiRYX+REZGFrrMyJEjAUBp/+Li4nD9+nV069YNFhYWSsu1bdtWIawCb07n161bFydOnJB/npKSknDixAlYWFhg6tSpCu2dnZ3x1VdfQSqVYuvWrSXe388++0whrALAkCFDoKOjo3CM5Ofn4+eff4aJiQmCgoIUeou1tLQwb948iEQihRq0tN5Ehne/kwCgVq1a762ttMtT1cY/Z6hSKhjXKRKJYGRkhCZNmqBfv36FBolGjRopncbOycnBlStXULNmTfz888+FbkNPTw/JyclIS0tDrVq18NdffwGAUiAG3vTafPLJJ0pj0wpz8+ZNZGZmokmTJsX2tpVGQa2FnQ4Xi8Vo1aoVdu3ahStXrij9w+zi4qK0TEGIKm6sprqpWkdx74uWlhZatWqF27dvK0wvGNv64MGDQnvnCtrfvHlTad67f6yok6r7XJr6VTFz5kxERUUhODgY/fv3f+/wE02aPn16kT2SZ8+ehaenp9J0R0dHtG7dGidOnMDdu3dha2sL4P8C7IgRIwpdX2HHmI6ODlq2bIkHDx7IP08FY+RbtWoFPT09pWXatWuH0NBQ+bFbEoUdI7q6ujA1NVU4Rm7duoXU1FTUq1cPS5YsKXRdBgYGSEhIkL/u378/oqOj0bFjR/Tu3Rtt2rRB8+bNVf6eKu3yVLUxsFKlVJLgVNgp0fT0dMhkMqSlpcnDb1Gys7NRq1Yt+UVYJiYmKm+nMJmZmQBQaA+OuhTUWlRNZmZmCu3eVlhPcUHvjKr3uS3oacnPzy+yTcG8grYfWseHvC9paWkAgOjoaERHRxdZ4/Pnz1Van7qous+lqV8VZmZmmDRpEn744QcsXboU8+bN+6D1CNmoUaNw7tw5hIeHY/bs2UhPT8eePXtQv359tGvXrtBlinrvC469gmOxNJ+/9ynqTI62tnahx8idO3fe+x1XwMPDAzt37sTKlSsRGRkpv0Wgk5MTpk+frnQGQ93LU9XGwEpVnkgkUppW8KXv5OSE2NhYldZTsExhV9gDwOPHj1VaT0FvVVE3aFeHglqLqiklJUWhXVltPz09vcg2Bf+glrb37kPel4JlIiIi4OXlVaLtFXY8lbfS1K+qiRMnIjw8HGFhYUX2OAJvfh9F/SFT8MeZEPXs2RMWFhbYuHEjAgMDsXnzZrx8+RLDhg0r8j0u6vNUcOwVvC+a/vy9ve5u3bphy5YtKi/XsWNHdOzYETk5Ofjjjz9w7NgxrFu3DsOGDVMaL10Wy1PVxTGsRIWoXr06nJyckJCQgNTUVJWWKbiIoODirrc9e/ZM4VZZxWnQoAFq1KiBf/75B/fv339v+5L2br5da2G3xJFKpfJTyu9eGKEuBXdMKOpWO2/PK8ndFQpT3PuSn59faA3NmzcHAJw/f75U21bV2xe7qEN51G9gYIBZs2ZBKpVi7ty5RbaTSCR48OBBofMKbs8mRDo6Ovjyyy/x5MkT7Nu3D+Hh4RCLxcVe2FjYMfb69Wv556ng3rUF/42Pjy/01lEFtx17+/T+h3zOi1PwPfPHH38UWsP7GBgYoHXr1pgzZw7mz58PmUymdKeKslyeqh4GVqIijB8/Hq9evcK4ceMK7Ql89uwZLl68KH/do0cPSCQSREVFKUwHgMWLF6t8ek9bWxujRo2CVCrF5MmTkZOTozA/Ly8Pjx49kr8uuDl3UaGgMP3794eenh7WrVunNI5x2bJlePjwIbp06QJzc3OV11kSbm5uqFevHv7++29EREQozb9y5Qo2btwIHR0dpQtISqply5ZwcHBAfHy80kVG69atUxq/Crx5L+vXr48NGzYU+Y/oX3/9Je8FLq2C91CVP1BUUV71DxgwAC4uLoiKisLly5cLbdO8eXM8ePAAR44cUZgeHh7+3vvgatqwYcOgq6uLb7/9Fjdv3oS3t3exN8M/c+YMDh8+rDBtzZo1ePDgAdq3by8f5mNpaYmOHTsiKSkJy5cvV2j/zz//YP369RCLxejfv798+od8zoujo6ODsWPH4smTJ/D398eLFy+U2qSmpir8oX3q1KlC2xX0COvr6xe7zdIuT1UbhwQQFWHw4MH466+/8Msvv8DFxQUdO3aEtbU1MjMzkZiYiNjYWLRv3x6bN28G8KZXdvny5Rg+fDh69uyJ3r17w8LCAufPn8f169fh5uam8vCCb775BpcuXcLx48fRtGlTdOvWDTVr1sTDhw9x9uxZDBkyRH4hyeeffw4tLS38/PPPSE9Pl4+LGz16dJGn062trbFo0SJMnToV7du3R69evWBmZob4+HjExMTA0tISS5cuVcNvsXDa2tr45Zdf0KdPH0yaNAmRkZFwdXWFrq4ubt68icOHDyMvLw+LFy9GvXr1SrUtkUiElStXonfv3hg+fLjCfVhPnjyJTp064dixYwrL6OrqYuPGjfDx8cGgQYPg6uqKJk2awNDQEElJSbhy5QoSEhJw5swZtVzd3L59e+zatQtff/01vL29YWhoiBo1amD06NEftL7yql8kEmHBggXw8PAoNPgDwKRJk3Ds2DEMGTIEvXr1gomJCS5fvozLly+ja9euSgFPSMzMzODh4YFdu3YBAL766qti23fv3h2DBw+Gl5cXbG1tceXKFRw7dgy1atVCSEiIQttly5ahW7duWLhwIc6cOYPmzZvL78Oak5OD5cuXKzzWtEWLFqhevTqioqKgp6eHunXrQiQSwdfXV+lerKoKCAjA9evXERERgSNHjqBt27awtLTE06dPcefOHcTFxWHkyJHyHuFZs2YhMTER7u7usLa2hr6+Pq5du4bjx4+jVq1a8rt0FKW0y1PVxsBKVIzFixejS5cuWLduHc6dO4f09HTUqFEDFhYWGDFiBPr166fQ3tvbGzt37sSiRYuwZ88e6Onpwc3NDUePHsWPP/6ocmDV09PDtm3bEB4ejsjISGzfvh2vX7+GmZkZ3N3d0b17d3lbe3t7rFu3DsuXL8fGjRvlPbLvu3p7+PDhqF+/PlauXIn9+/fj+fPnMDc3x+jRo+Hv71+mFw8Bb3rezp07h1WrVuHkyZNYu3Yt8vLyYGJiAi8vL4wZMwYtWrRQy7ZatWqFgwcPYv78+Th+/DiOHz+OZs2aYd++fTh+/LhSYAXejF+OiYnBmjVrcODAAURGRkImk8HMzAyOjo6YOHEiHBwc1FLfkCFDkJSUhG3btiE0NBSvXr2ClZXVBwfW8qy/devW6NGjR5E9ua1bt8bWrVsRHByM6Ohohc/Enj17BB1YgTfvza5du+Dk5IRWrVoV29bDwwPDhg1DSEgIDh06BF1dXXh7e2P27NlKt0+zsbHBqVOn5G3j4uJgaGgId3d3TJo0CW3atFFoX6NGDWzatAlBQUGIiopCdnY2gDfH9ocGVh0dHURERGDnzp3YtGkTjh49Kr+I1MrKClOmTMGAAQPk7adNm4b9+/fj0qVL8uFEFhYW8PPzw7hx4xQCdmFKuzxVbaKMjAw+C42IiKgQS5cuxfz58xESEiK/P+u7goKCsGjRIoSGhqr88A4iKhmOYSUiIipEdnY2fv31VxgbG5d6LDURlQ6HBBAREb3l4MGDuHTpEo4ePYpHjx5h9uzZMDIy0nRZRFUaAysREdFboqOjERkZCVNTU0yePBmTJk3SdElEVR7HsBIRERGRoHEMKxEREREJGgMrEREREQkaAysRERERCRoDq8AlJCRougSqIHiskKp4rJAqeJyQqsrjWGFgJSIiIiJBY2AlIiIiIkFjYCUiIiIiQWNgJSIiIiJBY2AlIiIiIkFjYCUiIiIiQWNgJSIiIiJBY2AlIiIiIkFjYCUiIiIiQWNgJSIiIiJBY2AlIiIiIkHT0XQBRESkmuRkPSQllb6fISvLFpmZ+mqoSP0sLfNhbp6r6TKISGAYWImIKoikJC34+2uXej1SqQHE4tKvpyyEhADm5pqugoiEhkMCiIiIiEjQGFiJiIiISNAYWImIiIhI0DQWWH/99Ve4ubnBysoKVlZW6Ny5Mw4fPiyf7+fnB4lEovDTqVMnTZVLRERERBqisYuuLCwsMHfuXNjZ2SE/Px+RkZEYPHgwTp06hUaNGgEA2rVrh7CwMPkyenp6miqXiIiIiDREY4G1Z8+eCq+/++47rFu3DhcuXJAHVrFYDDMzM02UR0REREQCIYgxrHl5edi5cyeeP3+OFi1ayKefP38e9vb2aNasGSZNmoQnT55osEoiIiIi0gRRRkaGTFMbv3btGrp06YKXL1/C0NAQv/76K7p27QoA2LlzJwwMDGBjY4PExEQsWLAA+fn5OHXqFMRicZHrTEhIKK/yiYjK1f37tpg500DTZZSphQtzYGV1V9NlEFE5c3BwKHa+RgNrbm4uHjx4gMzMTERHRyM8PBz79u2Dk5OTUtvk5GQ0btwY69evh5eXlwaq1YyEhIT3volEAI+VquDiRX01PThAWuwf/poUEpIHV9eXmi6DwO8UUl15HCsafdKVnp4e6tevDwD49NNP8eeff2L16tVYtWqVUltzc3NYWFjg9u3b5V0mEREREWmQIMawFsjPz0dubuHPkE5NTUVycjIvwiIiIiKqYjTWwzpnzhx06dIFlpaWyM7Oxo4dO3Du3Dls27YN2dnZCA4OhpeXF8zMzJCYmIh58+bBxMQEHh4emiqZiIiIiDRAY4E1JSUFo0ePxuPHj2FsbAxnZ2fs2LEDHTt2RE5ODq5fv44tW7YgMzMTZmZmaNOmDTZs2AAjIyNNlUxEREREGqCxwLpmzZoi5xkYGCAqKqocqyEiIiIioRLUGFYiIiIioncxsBIRERGRoDGwEhEREZGgMbASERERkaAxsBIRERGRoDGwEhEREZGgMbASERERkaAxsBIRERGRoDGwEhEREZGgMbASERERkaAxsBIRERGRoDGwEhEREZGgMbASERERkaAxsBIRERGRoDGwEhEREZGgMbASERERkaAxsBIRERGRoDGwEhEREZGgMbASERERkaAxsBIRERGRoDGwEhEREZGgMbASERERkaAxsBIRERGRoDGwEhEREZGgMbASERERkaAxsBIRERGRoDGwEhEREZGg6Whqw7/++is2bNiA+/fvAwAcHR3h7++Prl27AgBkMhmCg4MRHh6OjIwMNGvWDCEhIWjYsKGmSiYiAUtO1kNSUuX+GzwnR1vTJRARaYTGAquFhQXmzp0LOzs75OfnIzIyEoMHD8apU6fQqFEjLF++HKGhoQgNDYWDgwMWL16M3r1748KFCzAyMtJU2UQkUElJWvD3r9yBLjBQ0xUQEWmGxrojevbsic6dO6N+/fqwt7fHd999h+rVq+PChQuQyWRYs2YNJk+eDG9vbzg5OWHNmjXIzs7Gjh07NFUyEREREWmAIM6f5eXlYefOnXj+/DlatGiBe/fuISUlBR06dJC3MTAwgJubG+Lj4zVYKRERERGVN40NCQCAa9euoUuXLnj58iUMDQ2xceNGODs7y0OpiYmJQnsTExMkJycXu86EhIQyq1dTKuM+UdmoysdKVpYtpFIDTZdRpl690oZUmqeWdUmlUrWsR92ysnKQkHBX02XQ/1eVv1OoZEp7rDg4OBQ7X6OB1cHBAWfPnkVmZiaio6Ph5+eHffv2yeeLRCKF9jKZTGlaYeusTBISEirdPlHZqOrHSmamPsTiyj2GVVcXEItL/7UtlUohFovVUJH6GRvrVOnjWEiq+ncKqa48jhWNBlY9PT3Ur18fAPDpp5/izz//xOrVq+Hv7w8AePz4MerWrStv//TpU6VeVyIiIiKq3AQxhrVAfn4+cnNzYWNjAzMzM5w8eVI+7+XLlzh//jxatmypwQqJiIiIqLxprId1zpw56NKlCywtLeVX/587dw7btm2DSCSCn58fli5dCgcHB9jb2yMkJASGhobo27evpkomIiIiIg3QWGBNSUnB6NGj8fjxYxgbG8PZ2Rk7duxAx44dAQBff/01cnJyEBAQIH9wQFRUFO/BSkRERFTFaCywrlmzptj5IpEIgYGBCOSdsomIiIiqNEGNYSUiIiIiehcDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJmsqBNSYmBk+fPi1yfmpqKmJiYtRSFBERERFRAZUDq6enJ06ePFnk/NOnT8PT01MtRRERERERFVA5sMpksmLn5+bmQkuLIwyIiIiISL10ipuZlZWFzMxM+eu0tDTcv39fqV1GRgZ27twJc3Nz9VdIRERERFVasYF19erVWLx4MQBAJBIhMDAQgYGBhbaVyWT47rvv1F8hEREREVVpxQbWdu3aQV9fHzKZDPPmzYOPjw8aN26s0EYkEqFatWr49NNP4erqWqbFEhEREVHVU2xgbdWqFVq1agUAkEql8PT0hLOzc7kURkREREQEvCewvm3GjBllWQcRERERUaGKDKyRkZEAgAEDBkAkEslfv8/AgQPVUxkREREREYoJrOPGjYNIJEKfPn2gp6eHcePGvXdlIpGIgZWIiIiI1KrIwPrXX38BAPT09BReExERERGVpyIDq7W1dbGviYiIiIjKAx9NRURERESCpvJdAgDg1KlTCA8Px927d5Genq70uFaRSITLly+rsz4iIiIiquJUDqxr1qzBzJkz8dFHH8HV1RUNGzYsy7qIiIiIiACUILCGhobC3d0dO3fulF+IRURERERU1lQew5qamgofHx+1hdVly5ahffv2sLKygp2dHXx9fXH9+nWFNn5+fpBIJAo/nTp1Usv2iYiIiKhiULmH1cXFBYmJiWrb8Llz5zBixAg0bdoUMpkMP/zwA3r16oX4+HjUrFlT3q5du3YICwuTv2bvLhEREVHVonJgXbhwIQYOHIj27dujbdu2pd5wVFSUwuuwsDBYW1sjLi4O3bt3l08Xi8UwMzMr9faIiIiIqGJSObAGBQXB2NgYvXr1gp2dHaysrKCtra3QRiQSYdu2bR9USHZ2NvLz8yGRSBSmnz9/Hvb29qhRowbc3d3x3XffwcTE5IO2QUREREQVj8qB9caNGxCJRKhbty6kUilu3bql1EYkEn1wITNmzEDjxo3RokUL+bROnTrB09MTNjY2SExMxIIFC+Dl5YVTp05BLBYXup6EhIQPrkGoKuM+UdmoysdKVpYtpFIDTZdRpl690oZUmqeWdUmlUrWsR92ysnKQkHBX02XQ/1eVv1OoZEp7rDg4OBQ7X+XAevXq1VIVUpxvv/0WcXFxOHTokEKvbZ8+feT/7+zsDBcXFzRu3BiHDx+Gl5dXoet63w5XNAkJCZVun6hsVPVjJTNTH2Kx9vsbVmC6uoBYXKLbZxdKKpUW+Ue/phkb61Tp41hIqvp3CqmuPI6V0n/zlVJgYCCioqKwd+9e2NraFtvW3NwcFhYWuH37dvkUR0REREQap3JgvX//vkrtrKysVN749OnTERUVhX379qFBgwbvbZ+amork5GRehEVERERUhagcWD/55BOVxqimpaWptD5/f39s3boVGzduhEQiQUpKCgDA0NAQ1atXR3Z2NoKDg+Hl5QUzMzMkJiZi3rx5MDExgYeHh6plExEREVEFp3JgXbVqlVJgzcvLw71797BlyxaYmppi5MiRKm947dq1AABvb2+F6dOnT0dgYCC0tbVx/fp1bNmyBZmZmTAzM0ObNm2wYcMGGBkZqbwdIiIiIqrYVA6sgwcPLnLe5MmT0aFDB2RnZ6u84YyMjGLnGxgYKN2rlYiIiIiqHpUfzVqc6tWrY/DgwVi9erU6VkdEREREJKeWwAoAurq6SE5OVtfqiIiIiIgAqCmwXr16FT///DM+/vhjdayOiIiIiEiu1HcJyMzMRFZWFqpXr47Q0FC1FkdEREREpHJgdXd3VwqsIpEIEokE9evXR58+fSCRSNRdHxERVSEikTYuXtTXdBllytIyH+bmuZoug6hCUTmwrlmzpizrICIiwtOnQFBQ5X7EbkgIYG6u6SqIKha1XXRFRERERFQWGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNBUCqwvX77EokWLcOLEibKuh4iIiIhIgUqBVV9fHz/++CMePHhQ1vUQERERESlQeUhA48aNcfv27bKshYiIiIhIicqB9fvvv0dERAQOHz5clvUQERERESlQ+UlXK1asgEQiwcCBA2FhYQFbW1sYGBgotBGJRNi2bZvaiyQiIiKiqkvlwHrjxg2IRCLUrVsXAJCYmKjURiQSqa8yIiIiIiKUILBevXq1LOsgIiIiIioU78NKRERERIJWosCal5eHbdu2YcKECfD19cXff/8NAMjIyMCuXbvw6NGjMimSiIiIiKoulQNrZmYmunTpgjFjxmDPnj04evQoUlNTAQBGRkaYOXMmfvnllzIrlIiIiIiqJpUD69y5c3Hjxg1s374dly9fhkwmk8/T1taGp6cnjh49WiZFEhEREVHVpXJg3b9/P0aPHo1OnToVejcAOzs73L9/X63FERERERGpHFgzMjJQr169IufLZDLk5uaqpSgiIiIiogIqB1Zra2tcv369yPkxMTGwt7dXS1FERERERAVUDqz9+vVDREQEYmJi5NMKhgaEhYVh3759GDRokPorJCIiIqIqTeUHB0yZMgUXL16El5cX7O3tIRKJMGPGDKSlpSElJQU9e/bEmDFjyrJWIiIiIqqCVA6surq62LZtG7Zv347du3dDJBLh9evXaNKkCXx8fNC/f38+mpWIiIiI1E7lwFqgX79+6NevX6k3vGzZMuzduxe3bt2Cnp4eXF1dMXv2bDg5OcnbyGQyBAcHIzw8HBkZGWjWrBlCQkLQsGHDUm+fiIiIiCqGD3o0699//42DBw/i4MGD+PvvvxXuyaqqc+fOYcSIETh8+DCio6Oho6ODXr16IT09Xd5m+fLlCA0NxaJFi3DixAmYmJigd+/eePbs2YeUTUREREQVUIl6WHfu3InZs2fj4cOH8pAqEolgYWGB2bNnl6jnNSoqSuF1WFgYrK2tERcXh+7du0Mmk2HNmjWYPHkyvL29AQBr1qyBg4MDduzYgeHDh5ekdCIiIiKqoFQOrJs2bcKECRPg4OCAuXPnwt7eHjKZDP/99x8iIiIwZswY5ObmYvDgwR9USHZ2NvLz8yGRSAAA9+7dQ0pKCjp06CBvY2BgADc3N8THxzOwEhEREVURKgfWZcuWoVmzZti3bx/09fUV5o0aNQo9evTAsmXLPjiwzpgxA40bN0aLFi0AACkpKQAAExMThXYmJiZITk4ucj0JCQkftH0hq4z7RGWjKh8rWVm2kEoNNF1GmXr1ShtSaZ5a1iWVStWyHnVT5z4KVVZWDhIS7mq6DJVU5e8UKpnSHisODg7Fzlc5sCYlJWH06NFKYRUA9PX14evrizlz5pS4QAD49ttvERcXh0OHDkFbW1th3rt3HpDJZMXejeB9O1zRJCQkVLp9orJR1Y+VzEx9iMXa729YgenqAmJxia+VVSKVSiEWi9VQkfqpax+FzNhYp0J8Vqv6dwqprjyOFZUvunJ0dCy2Z/Phw4f4+OOPS1xAYGAgdu7ciejoaNja2sqnm5mZAQAeP36s0P7p06dKva5EREREVHmpHFjnzZuH8PBw7Nq1S2nezp07ERERgfnz55do49OnT8eOHTsQHR2NBg0aKMyzsbGBmZkZTp48KZ/28uVLnD9/Hi1btizRdoiIiIio4lL5vMvKlStRu3ZtjBgxAjNmzEC9evUgEolw+/ZtPHnyBHZ2dlixYgVWrFghX0YkEmHbtm2Frs/f3x9bt27Fxo0bIZFI5GNWDQ0NUb16dYhEIvj5+WHp0qVwcHCAvb09QkJCYGhoiL59+5Zyt4mIiIioolA5sN64cQMikQh169YF8GYIAACIxWLUrVsXUqkU//77r8IyxY01Xbt2LQDIb1lVYPr06QgMDAQAfP3118jJyUFAQID8wQFRUVEwMjJStWwiIiIiquBUDqxXr15V64YzMjLe20YkEiEwMFAeYImIiIio6vmgJ10REREREZUXBlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNJUDa5MmTXDgwIEi5x86dAhNmjRRS1FERERERAVUDqyJiYl4/vx5kfOfP3+O+/fvq6UoIiIiIqICJRoSUNyTq27dusUnUBERERGR2hX7pKvNmzcjMjJS/jokJATh4eFK7TIyMnD9+nV07dpV/RUSERERUZVWbGB9/vw5UlJS5K8zMzORn5+v0EYkEqFatWr48ssvMWPGjLKpkoiIiIiqrGID66hRozBq1CgAwCeffILg4GD06NGjXAojIiIiIgLeE1jfduXKlbKsg4iIiIioUCoH1gLPnj3DgwcPkJ6eDplMpjTf3d1dLYUREREREQElCKzp6emYPn06du3ahby8PKX5MpkMIpEIaWlpai2QiIiIiKo2lQPrlClTsG/fPowaNQru7u6QSCRlWBYRERER0RsqB9Zjx45hzJgxWLhwYVnWQ0RERESkQOUHB+jp6cHOzq4sayEiIiIiUqJyYPX29sbRo0fLshYiIiIiIiUqB9aJEyfi0aNHGDt2LC5cuIBHjx7hyZMnSj9EREREROqk8hjWZs2aQSQS4fLly9i2bVuR7XiXACIiIiJSJ5UD6zfffAORSFSWtRARERERKVE5sAYGBpZlHUREREREhVJ5DOvb8vLykJaWhtevX6u7HiIiIiIiBSUKrH/++Sd69eoFCwsL2NvbIyYmBgCQmpqK/v374/Tp02VSJBERERFVXSoH1t9//x09evTAnTt3MGDAAMhkMvm82rVrIzs7G7/99luZFElEREREVZfKgXX+/Pmws7NDfHw8vv/+e6X5bdq0wcWLF9VaHBERERGRyoH1zz//xJAhQ6Cvr1/o3QIsLS2RkpJSoo3HxMRgwIABaNiwISQSCTZt2qQw38/PDxKJROGnU6dOJdoGEREREVVsKt8lQEtLC1paRefblJQUGBgYlGjjz58/h5OTEwYOHIixY8cW2qZdu3YICwuTv9bT0yvRNoiIiIioYlO5h9XFxQWHDh0qdF5ubi62b9+OFi1alGjjXbp0wffffw9vb+8iw7BYLIaZmZn8p2bNmiXaBhERERFVbCoH1qlTp+LMmTOYMGECrl69CgB49OgRjh07Bi8vL9y5cwfTpk1Te4Hnz5+Hvb09mjVrhkmTJvHxr0RERERVjMpDAtq3b4+wsDAEBARg8+bNAN6MMZXJZKhRowbWrl2L5s2bq7W4Tp06wdPTEzY2NkhMTMSCBQvg5eWFU6dOQSwWF7pMQkKCWmsQgsq4T1Q2qvKxkpVlC6m0ZMOSKppXr7QhleapZV1SqVQt61E3de6jUGVl5SAh4a6my1BJVf5OoZIp7bHi4OBQ7HyVAysA9O3bFz169MDJkyfx33//IT8/H/Xq1UPHjh1RvXr1UhVamD59+sj/39nZGS4uLmjcuDEOHz4MLy+vQpd53w5XNAkJCZVun6hsVPVjJTNTH2KxtqbLKFO6uoBYXKKv7UJJpdIi/+jXNHXto5AZG+tUiM9qVf9OIdWVx7FS4m+FatWqoWfPnmVRy3uZm5vDwsICt2/f1sj2iYiIiKj8qTyG9cCBAwgICChyfkBAQJEXZalLamoqkpOTYWZmVqbbISIiIiLhUDmwrly5Ei9evChy/suXL7F8+fISbTw7OxtXrlzBlStXkJ+fjwcPHuDKlSu4f/8+srOzMWvWLPz++++4d+8ezp49iwEDBsDExAQeHh4l2g4RERERVVwqB9br16/DxcWlyPlNmjTBjRs3SrTxS5cuoW3btmjbti1ycnIQFBSEtm3b4ocffoC2tjauX7+OQYMGwdXVFX5+frC3t8eRI0dgZGRUou0QERERUcWl8hjW169fIycnp8j5OTk5Jb7qtE2bNsjIyChyflRUVInWR0RERESVj8qB1cnJCdHR0ZgwYYLSTf7z8/MRHR0NR0dHtRdIRKWXnKyHpCSVT6hUSDk5lfsOAUREVZnKgXXs2LEYOXIkBg4ciMDAQDRs2BAA8M8//yA4OBh//PEH1qxZU2aFEtGHS0rSgr9/5Q50gYGaroCIiMqKyoG1T58+uHPnDoKCgnD06FEAgEgkgkwmg0gkwvTp0+Hr61tmhRIRERFR1VSi+7D6+/ujb9++2Lt3L+7evQuZTIZ69erB09MTtra2ZVQiEREREVVlKgXWnJwc9O/fH76+vhgyZAgmTpxY1nUREREREQFQ8bZWBgYG+Ouvv5CXV7mf70xEREREwqPyZcOtW7dGbGxsWdZCRERERKRE5cC6aNEi/Pnnn/juu+9w9+5d5Ofnl2VdREREREQASnDRVfPmzSGTyRAaGorQ0FBoaWlBV1dXoY1IJMLDhw/VXiQRERERVV0qB9bevXtDJBKVZS1EREREREpUDqx8KAARERERaULlflYjEREREVV4JQqsiYmJmDRpElxcXGBlZYVz584BAFJTUzFt2jRcvny5LGokIiIioipM5SEB//77L7p164b8/Hy4uroiMTFRfl/W2rVr48KFC5BKpVi1alWZFUtEREREVY/KgXX27NkwMjLCsWPHoK2tDXt7e4X5Xbp0we7du9VdHxERERFVcSoPCYiNjcXIkSNhampa6N0CrKyskJycrNbiiIiIiIhUDqyvX7+GoaFhkfPT09Ohra2tlqKIiIiIiAqoHFidnJxw9uzZQufJZDLs3bsXLi4u6qqLiIiIiAhACQKrn58f9uzZg8WLFyMtLQ0AkJ+fj5s3b+Krr77CpUuXMHHixDIrlIiIiIiqJpUvuurTpw/u37+PhQsXIjg4WD4NALS1tbFgwQJ07ty5bKokIiIioipL5cAKAJMnT0bfvn0RHR2N27dvIz8/H/Xq1YOXlxdsbGzKqkYiIiIiqsLeG1ilUikOHDiAu3fvolatWujatSvGjRtXHrURERERERUfWFNSUtCjRw/cuXMHMpkMAGBoaIitW7fC3d29XAokIiIioqqt2IuuFixYgLt372LcuHHYunUrgoKCIBaL8c0335RXfURERERUxRXbw3rixAkMHDgQCxYskE8zNTXFyJEjkZSUBEtLyzIvkIiIiIiqtmJ7WFNSUtCyZUuFaa1atYJMJsODBw/KtDAiIiIiIuA9gTUvLw/6+voK0wpev3z5suyqIiIiIiL6/957l4C7d+/ijz/+kL/OysoCACQkJKB69epK7Zs1a6bG8oiIiIioqntvYA0KCkJQUJDS9HcvvJLJZBCJRPKnYKkiJiYGK1euxF9//YXk5GSEhoZi8ODBCusMDg5GeHg4MjIy0KxZM4SEhKBhw4Yqb4OIiIiIKrZiA2toaGiZbvz58+dwcnLCwIEDMXbsWKX5y5cvR2hoKEJDQ+Hg4IDFixejd+/euHDhAoyMjMq0NiIiIiIShmID66BBg8p04126dEGXLl0AQOlhBDKZDGvWrMHkyZPh7e0NAFizZg0cHBywY8cODB8+vExrIyIiKgsikTYuXtR/f0MNy8qyRWZmyeu0tMyHuXluGVREVVmJHs1anu7du4eUlBR06NBBPs3AwABubm6Ij49nYCUiogrp6VMgKEhb02W8l1RqALG45HWGhADm5mVQEFVpgg2sKSkpAAATExOF6SYmJkhOTi5yuYSEhDKtSxMq4z5R2SjqWMnKsoVUalDO1ZSvV6+0IZXmabqMMqXOfZRKpWpZj7rxfRSWDzlOsrJykJBwV/3FkKCVNqs4ODgUO1+wgbWASCRSeF1wcVdR3rfDFU1CQkKl2ycqG8UdK5mZ+h/UU1KR6OoCYrHgv9JKRV37KJVKIRaL1VCR+vF9FI4PPU6MjXX471YVUx5Zpdj7sGqSmZkZAODx48cK058+farU60pERERElZdgA6uNjQ3MzMxw8uRJ+bSXL1/i/PnzSk/fIiIiIqLKS6PnJLKzs3H79m0AQH5+Ph48eIArV66gZs2asLKygp+fH5YuXQoHBwfY29sjJCQEhoaG6Nu3rybLJiIiIqJypNHAeunSJXh6espfFzykYODAgVizZg2+/vpr5OTkICAgQP7ggKioKN6DlYiIiKgK0WhgbdOmDTIyMoqcLxKJEBgYiMDAwPIrioiIiIgERbBjWImIiIiIAAZWIiIiIhI4BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjRBB9agoCBIJBKFnwYNGmi6LCIiIiIqRzqaLuB9HBwcsG/fPvlrbW1tDVZDREREROVN8IFVR0cHZmZmmi6DiIiIiDRE0EMCAODu3bto2LAhPvnkE3z11Ve4e/eupksiIiIionIk6B5WV1dXrF69Gg4ODnj69CmWLFmCLl26IC4uDrVq1Sp0mYSEhHKusuxVxn2islHUsZKVZQup1KCcqylfr15pQyrN03QZZUqd+yiVStWyHnXj+ygsH3KcZGXlICHhrvqLIUErbVZxcHAodr6gA2vnzp0VXru6usLFxQWbN2/GhAkTCl3mfTtc0SQkJFS6faKyUdyxkpmpD7G4co//1tUFxGJBf6WVmrr2USqVQiwWq6Ei9eP7KBwfepwYG+vw360qpjyyiuCHBLytevXqcHR0xO3btzVdChERERGVkwoVWF++fImEhARehEVERERUhQj6nMSsWbPQrVs31K1bVz6G9cWLFxg4cKCmSyMiIiKiciLowPrw4UOMHDkSqamp+Oijj+Dq6oqjR4/C2tpa06VRJZKcrIekpAp1sqFQWVm2yMzUL3ReTk7lHr9KRESVm6AD6/r16zVdAlUBSUla8Pev+IFOKjUo8sKqwMByLoaIiEiNKn63EhERERFVagysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaAysRERERCRoDKxEREREJGgMrEREREQkaDqaLoCIiIgqD5FIGxcv6mu6jDJlaZkPc/NcTZdRpTCwEhERkdo8fQoEBWlruowyFRICmJtruoqqhUMCiIiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNAYWImIiIhI0CpEYF27di0++eQTmJmZ4fPPP0dsbKymSyIiIiKiciL4R7NGRUVhxowZWLp0KVq1aoW1a9eiX79+iIuLg5WVlabLq/SSk/WQlFQh/q75YDk5lfsRgkRERBWd4ANraGgoBg0ahC+//BIAsGTJEhw/fhzr16/H7NmzNVxd5ZeUpAV//8od6AIDNV0BERERFUeUkZEh03QRRcnNzYW5uTnWrVuHXr16yaf7+/vj+vXrOHDggOaKIyIiIqJyIehzvampqcjLy4OJiYnCdBMTEzx+/FhDVRERERFReRJ0YC0gEokUXstkMqVpRERERFQ5CTqw1q5dG9ra2kq9qU+fPlXqdSUiIiKiyknQgVVPTw8uLi44efKkwvSTJ0+iZcuWGqqKiIiIiMqT4O8SMH78eIwZMwbNmjVDy5YtsX79ejx69AjDhw/XdGlEREREVA4E3cMKAD4+PggKCsKSJUvQpk0bxMXFYdu2bbC2ttZ0aRohk8nQp08fSCQS7NmzR9PlkMCkp6cjICAAzZs3R506deDs7IypU6ciLS1N06WRAPAhLPQ+y5YtQ/v27WFlZQU7Ozv4+vri+vXrmi6LBG7p0qWQSCQICAgos20IPrACwMiRI3H16lU8fvwYp0+fhru7u6ZL0phVq1ZBW7ty3xeVPlxycjKSk5Mxd+5cxMbGIiwsDLGxsRgxYoSmSyMNK3gIy7Rp03DmzBm0aNEC/fr1w/379zVdGgnIuXPnMGLECBw+fBjR0dHQ0dFBr169kJ6erunSSKAuXLiA8PBwODs7l+l2BH0fVlJ06dIlDBkyBKdOnYKDgwPCw8Ph7e2t6bJI4I4cOQJfX1/cu3cPxsbGmi6HNKRjx45wdnbGihUr5NOaNm0Kb29vPoSFipSdnQ1ra2ts2rQJ3bt313Q5JDCZmZn4/PPPsXz5cixevBhOTk5YsmRJmWyrQvSwEvDs2TOMGDECP/74I++QQCXy7NkziMViVKtWTdOlkIbk5ubi8uXL6NChg8L0Dh06ID4+XkNVUUWQnZ2N/Px8SCQSTZdCAjR58mR4e3vj888/L/NtCf6iK3pj6tSp6NixI7p06aLpUqgCycjIwMKFCzF06FDo6PDjXlXxISz0oWbMmIHGjRujRYsWmi6FBCY8PBy3b99GWFhYuWyP/4Jp0IIFCxASElJsm7179yIpKQl///230u29qOpQ9Vhp06aN/PXz588xcOBAmJubY968eWVdIlUAfAgLlcS3336LuLg4HDp0iNdOkIKEhATMmzcPBw8ehJ6eXrlsk2NYNSg1NRWpqanFtqlbty6mTZuGLVu2QEvr/0Zw5OXlQUtLCy1atMChQ4fKulTSMFWPlYLT/tnZ2ejXrx8AYPv27ahevXqZ10jClZubC3Nzc6xbtw69evWST/f398f169dx4MABzRVHghQYGIioqCjs3bsXDRo00HQ5JDCbNm3C+PHjFf6QycvLg0gkgpaWFh4+fAixWKzWbTKwVgAPHz5ERkaGwjQ3NzcsXLgQPXv2hK2trUbqImF69uwZ+vXrB5lMhh07dsDIyEjTJZEAdOzYEY0aNcLy5cvl05o1awYvLy9edEUKpk+fjqioKOzbtw8ff/yxpsshAcrIyMDDhw8Vpo0fPx52dnaYOnUqGjZsqPazNxwSUAFYWFjAwsJCaXrdunUZVknBs2fP4OPjg2fPnmHTpk148eIFXrx4AQCoWbNmuZ26IeHhQ1hIFf7+/ti6dSs2btwIiUSClJQUAIChoSHP1JCcRCJRuhCvWrVqqFmzJpycnMpkmwysRJXI5cuXceHCBQBves/e9u4YV6pafHx8kJaWhiVLliAlJQUNGzas0g9hocKtXbsWAJRumTh9+nQEBgZqoiQiABwSQEREREQCx/uwEhEREZGgMbASERERkaAxsBIRERGRoDGwEhEREZGgMbASERERkaAxsBIRERGRoDGwElVQQUFBSjduBoDVq1fj008/Re3atdG4cePyL4yIiEjNGFiJBGDTpk3yJ4dIJBKYmZnB0dERPj4++Pnnn/Hs2TOV1nPu3Dl8++23aNKkCVauXImgoKAyrlz4GjdurPC7NTU1hYuLCwICAuRP8Smp7OxsBAUF4ezZs2quVphOnDiBIUOGwNHRESYmJrC2tkbHjh2xcOFCpcczViRnz56VHxcbN24stM2wYcPkn0ki0hw+6YpIQGbMmIF69erh1atXePz4Mc6dO4fAwECEhoYiMjISjRo1krcNCAjAlClTFJYvCFA//fRTob2vVZWzszMmTZoEAJBKpbh69SrCw8Nx5swZxMbGQltbu0Tre/78ORYtWgQAlfrpYTKZDNOmTcP69evh6OiIoUOHwsrKCjk5Obh8+TLCwsKwYcMG3Lp1S9Olloq+vj62b9+OIUOGKEzPysrCoUOHoK+vD5mMz9gh0iQGViIB6dixI5o3by5/PXXqVJw+fRoDBgzAwIED8fvvv8PAwAAAoKOjAx0dxY/w06dPAUCtYfXFixeoVq2a2tanCXXq1IGvr6/CtBo1aiAkJARXr16Fi4uLZgoTuNDQUKxfvx6jR49GcHAwtLQUT8oFBwdjxYoV711PTk6O/LgVoi5dumDfvn1ITk6Gubm5fHp0dDTy8/PRoUMHnDx5UoMVEhGHBBAJ3Oeff46AgADcv38f27Ztk09/dwyrRCLBunXr5P8vkUgUhgScPHkSHh4eqFu3LiwsLODh4YH4+HiFbRWs88aNGxg7dizq1auHVq1afdA6/vvvP0yZMgX16tWDpaUlvvzyS6SlpSnt38mTJ+Hp6QkrKyvUrVsXn3/+OSIiIhTaXLp0Cb6+vrC2tkadOnXQoUMHHDp0qOS/zLcUnOLV1dVVmP7o0SN8/fXXcHR0hKmpKZo2bYrly5fLe9ju3buHjz/+GACwaNEi+e/az88Pf//9NyQSCXbt2iVf3507dyCRSBR6xwFgypQpaNCgwQftZ1ZWFmbNmoXGjRvD1NQUjRo1wpw5cyCVShXaSSQSTJkyBUePHkWbNm1gZmaGpk2bYseOHe/9/eTk5GDp0qX4+OOPERQUpBRWAcDY2BizZs1SmNa4cWP06dMHZ86cQadOnWBmZoaffvoJAJCeno6pU6fi448/hqmpKVq0aIFVq1Yp9F7eu3cPEokEmzZtUtpe48aN4efnJ39dMJTmzJkzCAgIQP369WFpaYmhQ4fi0aNH793HAl27dkX16tWVfi/bt29Hp06dULNmzUKXU+XzkJiYiGnTpqF58+YwNzeHtbU1fH198c8//yi0KxiesGPHDqxatQqNGzeGmZkZOnfujL/++kuh7ePHjzFx4kQ4OzvD1NQUjo6O8PX1xbVr11TeZ6KKhoGVqAIo6B08ceJEkW3CwsLQtm1b+f+HhYXB09MTALBjxw706dMH2tramDlzJmbOnIm0tDR4eXnh4sWLSusaPnw40tPTMXPmTIwdO/aD1jFixAg8fPgQM2fOxNChQ7Fv3z588803Cm22bNkCHx8fPHr0CBMnTsTcuXPRrFkzHD58WN7m3Llz6NatGx4/foyAgADMnTsXenp6GDhwIKKjo1X6/b169QqpqalITU3Fw4cPceTIESxfvhwuLi5wcnKSt3vy5Ak6deqEw4cP48svv8SiRYvg6uqK2bNnIzAwEADw0UcfYcmSJQAADw8P+e96+PDhcHZ2hkQiQUxMjHydMTEx0NLSwoMHD3Dv3j359NjYWHz22Wcl3s+cnBx4eHjgt99+g4+PDxYvXoyuXbti1apVGD58uNK+X7hwAePHj0ePHj0wf/58VKtWDaNHj8a///5b7O8sPj4e6enp6Nu3b4mHTNy+fRtDhw6Fm5sbFi1ahObNm0MqlcLT0xPh4eHw8vLCwoULYWNjg1mzZuHbb78t0frfNWPGDFy+fBnffPMNhg0bhoMHD8LHxwe5ubkqLa+vrw9PT09s375dPu3Ro0c4e/Ys+vfvX+gyqn4eLl26hJiYGHh6eiIoKAh+fn64dOkSevToUegY6lWrViEyMhKjR4/GjBkzcOvWLQwePBivXr2St/nyyy+xZ88eDBw4ECEhIRgzZgzy8/Mr/NAMouJwSABRBWBpaQljY2PcuXOnyDa+vr6Ii4vDmTNnFE5/P3/+HP7+/vD19cWaNWvk04cPH45WrVph3rx5SsHP3t4ev/32W6nW0aBBA/zyyy/y1zKZDL/++iuWLl2KGjVqICsrC9988w2cnZ1x+PBhGBoaKrQt+O+UKVPQokUL7NmzR97LN2rUKHTt2hXff/89vLy83vv7O3PmDOzs7BSmNWvWDFu2bIFIJJJPW7BgAaRSKWJiYmBqairfxzp16mDVqlXw8/ODjY0NvLy8EBAQAGdnZ6WhBq1atUJsbKz89fnz59GhQwfEx8fj/PnzsLGxQWpqKm7evImvvvqqxPu5evVqJCQk4NSpU/KeXgBo2LAh/P39ERsbCzc3N/n0GzduICYmRt62V69eaNSoETZu3Ij58+cX+Tu7ceOGfL1vk8lkSj3lxsbGCj3Vd+7cwebNm9GjRw/5tF9++QV///03VqxYgaFDhwIARo4ciS+++AI///wzRo4cqfQelcS+ffsgFosBAI6Ojpg4cSI2b96MYcOGqbR8//79sWnTJty4cQOOjo7Yvn07qlevjm7duin8AQWU7PPQuXNneHt7Kyzv6+uLzz77DL/99hv8/f0V5mVlZSE2Nhb6+voAAAcHBwwZMgQnTpxA165dkZmZifPnz2P+/PmYOHGifLl3x7MTVTbsYSWqIKpXr47s7OwSL3fy5ElkZGSgf//+8l7G1NRU5OTkoF27djh//rxC7w3wpndU3etwd3dHXl4eHjx4IF9nVlYWpk2bphBWAchD5NWrV5GQkID+/fsjPT1dvt309HR06tQJd+/eRWJi4nt/B59++il2796N3bt3Y+vWrfj+++/x33//YciQIcjJyQHwJojt2bMHXbt2hba2tsJ+duzYEfn5+Qo9p0Vxc3PDP//8g/T0dABvelLbtGmDFi1ayINsTEwMZDKZvIe1JPu5a9cutGzZEh999JFCje3atQPwJpy/rU2bNgrB1tTUFA4ODrh7926x+1FwZwojIyOF6WlpabCzs1P4iYuLU2hjaWmpEFYB4PDhw6hduzYGDx4snyYSiTBp0iTIZDIcOXKk2HqKM3z4cHlYBYCBAweiRo0aJVpnmzZtYG5uLu9l3b59Ozw8POTB8W0l+Ty8Pf77xYsXSEtLQ40aNWBnZ4fLly8rrXvw4MEK22zdujUAyN8vfX196Orq4ty5c/JjjKgqYA8rUQWRnZ2Njz76qMTL/ffffwCA3r17F9kmMzNTYd22tralXoeVlZXC/ILxtgX/yBb0Fr99Sr6o2idOnKjQm/S2p0+fwtraush1AECtWrXkgQ54M2bRwcEBX3zxBSIiIjBmzBg8ffoUGRkZ2LhxY5G3OCq4qK04bm5ukMlkiI2NRbNmzXDnzh24ubnh9evX2LJlC4A3IdbY2Fg+rrUk+/nff//h77//LrI38t0a330fgDfvxfvCTkFQffeWasbGxti9ezeAN8MYQkJClJa1sbFRmpaYmAg7Ozul4QUFYVqVPzyK8u7vQkdHBzY2Nrh//77K69DS0oKPjw+2b9+Ofv364cqVK5g3b16hbUvyeXj58iV++OEHbNu2TWlcbe3atZWWe9/nRiwWY/bs2Zg9ezYcHBzg6uqKzp07o3///oW+10SVBQMrUQWQlJSErKws1K9fv8TL5ufnA3hzKtnCwqLQNsbGxgqv372i+0PWUdS4x7dP9wNQOCVfVO1z5swp8kp+e3v7IpcvTsF439jYWPkYQADo27ev0u2NCqjy+3dxcYGhoSFiY2MhlUpRrVo1uLi44NWrV5g/fz6ePHkiH79acOq/JPuZn5+Ptm3bYurUqYW2e/f9ed/7UBRHR0cAwPXr1+Hh4SGfrqurKw//qamphS5bmjsCqHI8qLLMh9yGql+/fggNDUVAQADq1KlT5C3LSvJ5mDFjBiIiIjB69Gi0atUKxsbG0NLSQmBgYKH7o8r7NWHCBHh4eODAgQM4deoUlixZgmXLlmHz5s34/PPPS7TPRBUFAytRBbB161YAQIcOHUq8bL169QC8uVjo7V7G8l7HuwrC3/Xr15Wuln93u9WrV1fbdgu8fv0awJvxiMCbfTM2Nsbr16/fu63iQpWOjg5cXV0RGxuL3NxcNG/eHLq6umjWrBnEYjEOHTqEa9euwcfHR75MSfazXr16yM7OVvvv410tW7ZEzZo1sWPHDkybNq3EF169y9raGn/99Rfy8vIU1nXz5k35fADyK/IzMzMVlpdKpUVe+X/r1i20b99e/vr169dITEyEu7t7iWp0cXFBgwYNcPbsWYwbN67IfS7J5yEqKgoDBgxAcHCwwvSMjAzUqlWrRPW9zdbWFuPGjcO4cePw4MEDtG3bFj/++CMDK1VaHMNKJHCnT5/GkiVLYGNjU+QVy8Xp2LGj/J6j7972CFDtNLc61vGu9u3bw9jYGMuWLcOLFy8U5hX0Jrm4uMDOzg4rV65UCjAfut0CBeMbC07La2trw8vLC/v27St0bGFmZqbSuMSMjIxC1+3m5oYrV67g2LFj8gugxGIxmjZtihUrViAvL0/hwqiS7KePjw/+/PNPHDhwQKldTk7OB41zLoyBgQGmTZuGmzdvFtkbWJJezK5du+Lp06eIjIxUWH7lypUQiUTo0qULgDdDET766COlp4itX78eeXl5ha57w4YNCsdlZGQkMjMz0blzZ5XrK7Bw4UJMnz5daQz220ryedDW1lb6Pe3YsQPJycklrg14Mw62YNx1gbp168LExKTI45GoMmAPK5GAHD9+HLdv38br16/x5MkTnDlzBidPnoSVlRUiIyMLvQDkfYyMjLB8+XKMGDECrVu3Rr9+/WBmZoakpCScPXsWhoaG770vpzrW8S5jY2MEBQVhwoQJaN++Pfr164datWrhn3/+QXJyMjZu3AgtLS2sWrUKffr0QatWrTB48GBYW1vj0aNHuHDhAu7fv690wU9hHj16JO+lzs3NxbVr1/C///0PtWvXxujRo+Xt5syZg5iYGHTr1g1ffPEFnJyc8OzZM1y/fh179+7Fn3/+CTMzM1SvXh0ODg6IioqCvb09atWqBRsbG7i6ugIAPvvsM+Tl5cnHrxZwd3dHSEgIDAwMFE79l2Q/J06ciCNHjuCLL75A//790axZM0ilUty6dQu7du3C9u3bFR4+URrjx4/Hf//9h19++QVnzpyBl5cXrKys8Pz5c/z777/YuXMnDA0NCx2L+a6hQ4ciIiICkydPxtWrV2Fvb4+jR4/iyJEjGDt2rMI41GHDhiEkJATjxo1D8+bNcenSJZw+fbrY7Xh6eqJPnz5ITEzEL7/8AkdHRwwaNKjE+9y5c+f3Bt2SfB66d++OLVu2wMjICE5OTrh69SqioqKUxomr6tatW/Dy8kKvXr3g6OgIsViMI0eO4N9//y32rg9EFR0DK5GAFJw21NPTQ82aNeHk5ISgoCAMHjxY6WrtkujVqxfMzc2xbNkyrF69Gjk5OTAzM4Orq6v8FkPlsY53DR48GCYmJvjxxx+xbNkyaGtrw87ODiNHjpS3+eyzz3D8+HEsXrwY//vf/5CVlQUTExM0atRIfm/U97l27RrGjBkD4E04rF27Njw8PDBz5kyFMYgfffQRjh8/jiVLlmD//v343//+hxo1asDe3h4zZsxQuIF8aGgoAgMDMWvWLEilUgwcOFAeWJs3bw49PT0AkE8r2JeCaQXzS7qfBgYGiI6OxvLlyxEVFSUPjba2tvDz84ODg4NKvxNViEQi/Pjjj+jZsyfWr1+P8PBwpKamolq1arC3t8fYsWMxbNiwIsdxvk1fXx/R0dGYP38+du3ahfT0dNjY2GD+/PmYMGGCQlt/f3+kpaUhKioKu3fvRuvWrbFnzx75fYXfFRwcjOjoaCxatAhSqRRdu3bFkiVLFO4coG6qfh6Cg4Ohq6uLXbt2YePGjXBxccHOnTvx3XfffdB269ati379+uHMmTPYsWMHRCKRvHf+iy++UNfuEQmOKCMjgw9IJiKiCmfTpk0YP348jh49qrZeZSISJo5hJSIiIiJBY2AlIiIiIkFjYCUiIiIiQeMYViIiIiISNPawEhEREZGgMbASERERkaAxsBIRERGRoDGwEhEREZGgMbASERERkaAxsBIRERGRoP0/9VSKAxoRwNUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "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": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 21,
     "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": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6UAAAFDCAYAAADCo6jwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABD30lEQVR4nO3dfVxUdf7//+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": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "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": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.8076725017901509"
      ]
     },
     "execution_count": 23,
     "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": 24,
   "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": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observed Difference: -0.8076725017901509\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAFuCAYAAAB9QTkMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABOmUlEQVR4nO3dd1xT1/8/8FcYAg6MIsQFqEClTioqiqNuHCCKA7dSJ666UKmt2yKKtg6KWsdH66gLK466URAUR0WsE+tEBZXpgCiQ3x9+uT9jGAkGQuLr+Xj4aHPvuTfv3NwkL849915RSkqKDEREREQ6Qk/TBRARERGpE8MNERER6RSGGyIiItIpDDdERESkUxhuiIiISKcw3BAREZFOYbihEsfb2xtisRjh4eFy08ViMbp161Zkz+vn55fr89IHOdtn27Ztmi4lT926dYNYLMbDhw81XUqxyO0zwf04f9u2bYNYLIafn5+mSwHA96uoMNx8ocRisdy/ihUrwtraGp07d8amTZuQlZWl6RLVTht+nD+VU3N+X8Th4eFFHvxKipIeXurXrw+xWIxq1arh2bNnubb57rvvStSPma7tYyUtvJBmGGi6ANKsGTNmAACysrJw//59HDx4EOfPn8fp06exefNmDVcn78KFCzAxMSmy9Y8aNQq9evVC9erVi+w56Mvw5s0bLFy4EIGBgZouhUo4fu8UDYabL5yvr6/c4+vXr6NDhw7Yv38/IiMj4ezsrKHKFH311VdFun4zMzOYmZkV6XPQl8HGxgY7duzA6NGj0aBBA02XQyUYv3eKBg9LkZy6deuiRYsWAIDLly8D+P9d0t7e3rh16xYGDRqEWrVqQSwWIyYmRlh2//79cHd3R40aNWBhYYFGjRph7ty5SEtLy/W5Tp8+jS5duqBq1aqoUaMGBgwYgNu3b+dZW17d4llZWdiyZQu6dOkCa2trSCQSNGjQACNGjMCVK1cAfDic4e/vDwAYN26c3CG5nEMc+R37DgsLQ58+fVCzZk1YWFigYcOGmDFjBl68eKHQ9uMxQ/v370e7du1QpUoV1KhRA15eXnjy5Emer1GdCltHdHS08JekpaUl3N3dERUVle9zJSQkYObMmWjUqBEkEgmsra3Rs2dPnDlzRqHtx4cNoqKi4OHhAWtra4jFYqSkpOT5HGKxGBEREQCAhg0bCu9f/fr1c22/adMmODs7QyKRwM7ODhMnTsxz/arUr4y5c+ciOzsbP/74o9LL5HfYpySOy/j9998hFouxePHiXOenpaWhatWqqFu3rnCY++NDw0eOHEHHjh2Fz/+wYcNw//79XNeVkJAAHx8fNGzYEBYWFqhZsyb69u2Ls2fPyrXz9vbGuHHjAAD+/v5yn/Pctl1MTAz69u0LKysrVKlSBV26dMH58+dzrSE7OxtbtmyBi4sLrKysIJFI0Lx5cyxfvhzv3r1TaB8eHg5PT0/UrVsXFhYWsLW1RZs2bTBr1izIZP//rkd5vbfKLk+5Y88NKe3+/fvo1KkTateujX79+iE1NRWlS5cGAEydOhUbNmxAtWrV4OrqCrFYjEuXLuHXX3/FsWPHcPToUZQrV05Y1/79++Hl5QVDQ0P06NEDVatWxfnz59GxY0fUq1dP6ZrevXuHAQMG4MSJE6hcuTJ69uyJChUqIC4uDuHh4bCxscE333yDAQMGAAAiIiLQtWtXuR/E8uXL5/scmzZtwpQpU2BiYgJ3d3dUrlwZUVFRWLt2LQ4dOoS///4blpaWCstt2LABf//9N7p27YoWLVrg0qVL2LdvH65du4aIiAgYGRkp/To/hyp1REVFoUePHpBKpXBzc4ONjQ2uX78ONzc3tG7dOtf1X79+HT179sSLFy/Qrl07dO3aFUlJSTh06BB69OiBlStXYvDgwQrLXbhwAcuXL4ezszOGDBmCZ8+eQV9fP8/XMWPGDGzfvh2PHz/GmDFjhPctt/dvzpw5OHXqFDp37oy2bdsiPDwcW7Zswd27d3H48GG11J+fjh07om3btggNDcWRI0fQuXNnlZbXBv369cP8+fPxxx9/wMfHR+G9+/PPP/H27VtMnDhRYd6BAwdw4sQJuLm5oVWrVoiJicFff/2F8PBwHDt2DDY2NkLbhw8fokuXLnj69ClatGgBDw8PxMfH46+//sKJEyfw66+/YsiQIQA+/BGTmpqKw4cPo0WLFmjZsqWwHisrK7kaoqOjsXLlSjg5OWHIkCGIi4tDSEgI3N3dERYWhtq1awttMzMzMWjQIBw5cgS2trbo1asXjIyMEBERgfnz5+PMmTPYu3cvDAw+/KQeO3YMnp6eKFeuHLp06YJq1aohJSUF//33H9auXYt58+YJbXPzucsTww194ubNm8Jfx40aNZKbd/78eUyZMgWzZ8+Wm75z505s2LABrq6u+P333+XGxSxduhSLFi2Cn58ffv75ZwDA69evMWnSJIhEIhw6dAiNGzcW2v/0009YtWqV0vX6+/vjxIkTaNOmDbZv3y6ELeBDj05Oz8rAgQPx6NEjREREoFu3bhg4cKBS63/06BFmzJiB0qVL48SJE/j666+FeQsXLkRAQACmTp2KXbt2KSx76tQpnDlzBvb29sK0ESNGYM+ePTh06BA8PDyUfp2fQ9k6ZDIZxo8fj/T0dGzevBnu7u5C+99//x0+Pj4K687KysLQoUORmpqKAwcOyP2YxMfHo3379vDx8YGLiwssLCzklg0NDcWvv/6KYcOGKfU6fH19cfbsWTx+/Bje3t6wtrbOs+3ly5dx7tw5VKtWDcCHHyc3NzdERkbi0qVLwj73OfUXZMGCBWjdujVmz56NDh06lPgfo7Nnz+Y5CPfRo0cK08qVKwdPT0+sX78eR44cUeh1+t///gcDAwMheHzsyJEj2LlzJ1xcXIRpq1atwk8//QQfHx8EBwcL0ydPnoynT59i5syZmDlzpjB9/Pjx6NChA3x8fNCuXTtUr14drq6uQrhp2bKlwmH3jx09ehRr166Fp6enMG3Tpk2YPHky1q5di+XLlwvTf/nlFxw5cgQjR47E4sWLhbCWnZ2NyZMnY/PmzVi/fj3GjBkDANiyZQtkMhkOHDiAhg0byj1vUlJSgfvC5y5PPCz1xfPz84Ofnx8WLlyIkSNHom3btkhPT4erq6tweCqHhYWFMAD5Y7/99hv09fWxatUqhQG/U6ZMgZmZmdyP/+HDh5GcnAwPDw+5YAMA06dPh6mpqVK1Z2VlYf369TAyMsKKFSvkgg0A6Ovro3LlykqtKy+7du3Cu3fvMHz4cLlgAwA+Pj6oUqUKjh07hqdPnyosO3r0aLlAAQBDhw4FAPzzzz+fVZcqlK0jKioKsbGxcHJykgs2ADB8+HDUqlVLYd3Hjh3D3bt3MXz4cLlgAACVK1fGhAkTkJGRgf379yssW69ePaWDjaqmT58uBBsAMDAwwKBBgwDIv+bPqb8g9erVw6BBg3Dnzh1s2rSpkK+k+ERERMDf3z/Xfzt27Mh1mREjRgCAwus7f/48bty4gc6dO6Nq1aoKy7Vu3Vou2AAfDilVr14dp06dEj5PT548walTp1C1alVMmTJFrn3dunXx3XffQSqVYufOnSq/3ubNm8sFGwAYNGgQDAwM5PaR7OxsrFmzBubm5vDz85PrhdLT08P8+fMhEonkatDT+/DT+ul3EgBUrFixwNo+d3liz80XL2ccikgkQrly5dCwYUP06dMn1x+devXqKRxKSU9PR0xMDCpUqIA1a9bk+hylSpXCs2fPkJSUhIoVK+Lq1asAoBCegA9/DTZo0EDhWHpu7ty5g9TUVDRs2DDfv+I/R06tuR2SMTIyQrNmzbBv3z7ExMQofIk7ODgoLJPzg5vf2BJ1U7aO/N4XPT09NGvWDPfu3ZObnjMWJy4uLte/+nPa37lzR2Hep8FWnZR9zZ9TvzJmzZqF4OBgLF68GH379i3wEKgmzZgxI8+ejvDwcLi5uSlMt7e3R8uWLXHq1Ck8ePAANWrUAPD/w87w4cNzXV9u+5iBgQGcnJwQFxcnfJ5yxvQ1a9YMpUqVUlimTZs2CAwMFPZdVeS2jxgaGsLCwkJuH7l79y4SExNRs2ZNLF26NNd1mZiYIDY2Vnjct29fhISEoH379ujZsydatWqFJk2aKP099bnLE8PNF0+VH9ncuuWTk5Mhk8mQlJQkBKW8vH79GhUrVhQGGJubmyv9PLlJTU0FgFz/MlSXnFrzqkkikci1+1huPVA5f/Upex2hnL/gsrOz82yTMy+nbWHrKMz7kpSUBAAICQlBSEhInjW+efNGqfWpi7Kv+XPqV4ZEIsHEiRPx888/Y9myZZg/f36h1lOSjRw5EmfPnsXmzZsxZ84cJCcnY//+/ahVqxbatGmT6zJ5vfc5+17Ovvg5n7+C5NVDrK+vn+s+cv/+/QK/43K4urpi7969WLVqFXbs2CFcVqNOnTqYMWOGQs+oupcnhhtSgUgkUpiW8wVRp04dREZGKrWenGVyO9MIAJ4/f67UenL+Cs7rYmnqkFNrXjUlJCTItSuq509OTs6zTc6X7+f2ChTmfclZZsuWLejevbtKz5fb/lTcPqd+ZU2YMAGbN2/G2rVr8+zJAD5sj7xCb06QL4m6deuGqlWrYuvWrfD19cX27duRkZGBYcOG5fke5/V5ytn3ct4XTX/+Pl53586d8eeffyq9XPv27dG+fXukp6fj8uXLOHHiBDZs2IBhw4YpjO8qiuW/dBxzQ5+lbNmyqFOnDmJjY5GYmKjUMjkD5HIGLn/s1atXcqeX5+err75C+fLlcfPmTTx+/LjA9qr2mnxca26nkUqlUuGwxqeD/tQl58yxvE5P/XieKmeZ5Sa/9yU7OzvXGpo0aQIAOHfu3Gc9t7I+HsipDsVRv4mJCX788UdIpVLMmzcvz3ZisRhxcXG5zsu5pEFJZGBggKFDh+LFixc4ePAgNm/eDCMjo3wH7ee2j2VmZgqfp5xrA+X8NyoqKtfTrXNO1f/4EFNhPuf5yfmeuXz5cq41FMTExAQtW7bE3LlzsWDBAshkMoUz9opy+S8Vww19tnHjxuH9+/cYO3Zsrj0Mr169wqVLl4THXbt2hVgsRnBwsNx0AFiyZInSXcz6+voYOXIkpFIpJk2ahPT0dLn5WVlZiI+PFx7nXCgrrx+Q3PTt2xelSpXChg0bFMZdLF++HE+fPkWnTp1QpUoVpdepCmdnZ9SsWRP//vsvtmzZojA/JiYGW7duhYGBgcLgSFU5OTnBzs4OUVFRCgNoN2zYoDDeBvjwXtaqVQubNm3K8wv36tWrQu/S58p5D5UJs8oorvr79esHBwcHBAcHIzo6Otc2TZo0QVxcHI4dOyY3ffPmzQVeZ0jThg0bBkNDQ/zwww+4c+cO3N3d870wXVhYGI4ePSo3LSgoCHFxcWjbtq1wqLlatWpo3749njx5ghUrVsi1v3nzJjZu3AgjIyP07dtXmF6Yz3l+DAwMMGbMGLx48QLTpk3D27dvFdokJibK/VF2+vTpXNvl9DQZGxvn+5yfuzzxsBSpwcCBA3H16lWsW7cODg4OaN++PaysrJCamopHjx4hMjISbdu2xfbt2wF86O1ZsWIFvLy80K1bN/Ts2RNVq1bFuXPncOPGDTg7Oyt9iGv69Om4cuUKTp48iUaNGqFz586oUKECnj59ivDwcAwaNEgYJPntt99CT08Pa9asQXJysnAcf9SoUXke0rGysoK/vz+mTJmCtm3bokePHpBIJIiKikJERASqVauGZcuWqWEr5k5fXx/r1q1Dr169MHHiROzYsQONGzeGoaEh7ty5g6NHjyIrKwtLlixBzZo1P+u5RCIRVq1ahZ49e8LLy0vuOjehoaHo0KEDTpw4IbeMoaEhtm7dCg8PDwwYMACNGzdGw4YNUaZMGTx58gQxMTGIjY1FWFiYWs7yaNu2Lfbt24fvv/8e7u7uKFOmDMqXL49Ro0YVan3FVb9IJMLChQvh6uqaa0gEgIkTJ+LEiRMYNGgQevToAXNzc0RHRyM6OhouLi4KYaAkkUgkcHV1xb59+wB8uH9Wfrp06YKBAweie/fuqFGjBmJiYnDixAlUrFgRAQEBcm2XL1+Ozp07Y9GiRQgLC0OTJk2E69ykp6djxYoVcrcuaNq0KcqWLYvg4GCUKlUK1atXh0gkgqenp8K1bpTl4+ODGzduYMuWLTh27Bhat26NatWq4eXLl7h//z7Onz+PESNGCD1NP/74Ix49eoQWLVrAysoKxsbGuH79Ok6ePImKFSsKZyvm5XOXJ4YbUpMlS5agU6dO2LBhA86ePYvk5GSUL18eVatWxfDhw9GnTx+59u7u7ti7dy/8/f2xf/9+lCpVCs7Ozjh+/Dh++eUXpcNNqVKlsGvXLmzevBk7duzA7t27kZmZCYlEghYtWqBLly5CW1tbW2zYsAErVqzA1q1bhZ6egs5i8fLyQq1atbBq1SocOnQIb968QZUqVTBq1ChMmzatSAfGAh/+oj979ixWr16N0NBQrF+/HllZWTA3N0f37t0xevRoNG3aVC3P1axZM/z9999YsGABTp48iZMnT8LR0REHDx7EyZMnFcIN8GG8VUREBIKCgnD48GHs2LEDMpkMEokE9vb2mDBhAuzs7NRS36BBg/DkyRPs2rULgYGBeP/+PSwtLQsdboqz/pYtW6Jr16559hC1bNkSO3fuxOLFixESEiL3mdi/f3+JDjfAh/dm3759qFOnDpo1a5ZvW1dXVwwbNgwBAQE4cuQIDA0N4e7ujjlz5ihccsDa2hqnT58W2p4/fx5lypRBixYtMHHiRLRq1Uquffny5bFt2zb4+fkhODgYr1+/BvBh3y5suDEwMMCWLVuwd+9ebNu2DcePHxdOkLC0tMTkyZPRr18/of3UqVNx6NAhXLlyRTikXbVqVXh7e2Ps2LEF3kfqc5cnQJSSksLrOBMR0WdZtmwZFixYgICAAOH6N5/y8/ODv78/AgMDlb6QJlFhcMwNERF9ltevX+P333+HqanpZ4/9IlIHHpYiIqJC+fvvv3HlyhUcP34c8fHxmDNnjtw95Ig0heGGiIgKJSQkBDt27ICFhQUmTZqEiRMnarokIgAcc0NEREQ6hmNuiIiISKcw3BAREZFOYbghIiIincJw84WKjY3VdAlfFG7v4sNtXby4vYsPt7XyGG6IiIhIpzDcEBERkU5huCEiIiKdwnBDREREOoXhhoiIiHQKww0RERHpFIYbIiIi0ikMN0RERKRTGG6IiIhIpzDcEBERkU5huCEiIiKdYqDpAoio5HomFeHJm2xNl6GStFIWwFs9pGZkaboUlVQro4cqRjJNl0GkExhuiChPT95kY9rZF5ouQyVSaQbmtiwHv4svNV2KSgJamqOKkUjTZRDpBB6WIiIiIp3CcENEREQ6heGGiIiIdArDDREREekUhhsiIiLSKQw3REREpFMYboiIiEinMNwQERGRTmG4ISIiIp2i0XATERGBfv364euvv4ZYLMa2bdvk5stkMvj5+cHe3h6VK1dGt27dcPPmTbk2UqkUPj4+qFWrFqpWrYp+/frhyZMnxfkyiIiIqATRaLh58+YN6tSpg8WLF8PExERh/ooVKxAYGAh/f3+cOnUK5ubm6NmzJ169eiW08fX1xYEDB7BhwwYcPnwYr169gqenJ7KytOu+MkRERKQeGg03nTp1wuzZs+Hu7g49PflSZDIZgoKCMGnSJLi7u6NOnToICgrC69evsWfPHgBAamoq/vjjD8yfPx9t27aFg4MD1q5di+vXr+P06dMaeEVERESkaSV2zM3Dhw+RkJCAdu3aCdNMTEzg7OyMqKgoAEB0dDTev38v16Z69eqoXbu20IaIiIi+LCX2ruAJCQkAAHNzc7np5ubmePbsGQDg+fPn0NfXh5mZmUKb58+f57nu2NhYNVernbgdipc2bu+0UhaQSjM0XYbK3me+17q6016lITYx7++tkkwb921txW39gZ2dXb7zS2y4ySESieQey2QyhWmfKqhNQRvlSxAbG8vtUIy0dXunJslgZCTVdBkqkUozYGhgCCMjY02XohLTcqawq1he02WoTFv3bW3Eba28EntYSiKRAIBCD8zLly+F3hwLCwtkZWUhMTExzzZERET0ZSmx4cba2hoSiQShoaHCtIyMDJw7dw5OTk4AAAcHBxgaGsq1efLkCW7fvi20ISIioi+LRg9LvX79Gvfu3QMAZGdnIy4uDjExMahQoQIsLS3h7e2NZcuWwc7ODra2tggICECZMmXQu3dvAED58uUxePBgzJ49G+bm5qhQoQJmzZqFunXrok2bNhp8ZURERKQpGg03V65cgZubm/DYz88Pfn5+6N+/P4KCgvD9998jPT0dPj4+SElJgaOjI4KDg1GuXDlhmZ9//hn6+vrw8vJCRkYGWrdujTVr1kBfX18TL4mIiIg0TKPhplWrVkhJSclzvkgkgq+vL3x9ffNsY2xsjKVLl2Lp0qVFUCERERFpmxI75oaIiIioMBhuiIiISKcw3BAREZFOYbghIiIincJwQ0RERDqF4YaIiIh0CsMNERER6RSGGyIiItIpDDdERESkUxhuiIiISKcw3BAREZFOYbghIiIincJwQ0RERDqF4YaIiIh0CsMNERER6RSGGyIiItIpDDdERESkUxhuiIiISKcw3BAREZFOYbghIiIincJwQ0RERDqF4YaIiIh0CsMNERER6RSGGyIiItIpDDdERESkUxhuiIiISKcw3BAREZFOYbghIiIincJwQ0RERDqF4YaIiIh0CsMNERER6RSGGyIiItIpDDdERESkUxhuiIiISKcw3BAREZFOYbghIiIincJwQ0RERDqF4YaIiIh0CsMNERER6RSGGyIiItIpDDdERESkUxhuiIiISKeU6HCTlZWFhQsXokGDBpBIJGjQoAEWLlyIzMxMoY1MJoOfnx/s7e1RuXJldOvWDTdv3tRg1URERKRJSoebiIgIvHz5Ms/5iYmJiIiIUEtROX799VesX78e/v7+uHDhAhYvXozff/8dy5cvF9qsWLECgYGB8Pf3x6lTp2Bubo6ePXvi1atXaq2FiIiItIPS4cbNzQ2hoaF5zj9z5gzc3NzUUlSOCxcuoHPnzujSpQusra3RtWtXdOnSBZcvXwbwodcmKCgIkyZNgru7O+rUqYOgoCC8fv0ae/bsUWstREREpB2UDjcymSzf+e/evYOennqPcjVr1gxnz57FnTt3AAC3bt1CeHg4OnbsCAB4+PAhEhIS0K5dO2EZExMTODs7IyoqSq21EBERkXYwyG9mWloaUlNThcdJSUl4/PixQruUlBTs3bsXVapUUWtxkyZNwuvXr+Hk5AR9fX1kZmZi2rRpGDFiBAAgISEBAGBubi63nLm5OZ49e5bnemNjY9Vap7bidihe2ri900pZQCrN0HQZKnuf+V7r6k57lYbYxOeaLqNQtHHf1lbc1h/Y2dnlOz/fcPPbb79hyZIlAACRSARfX1/4+vrm2lYmk+Gnn34qZJm5Cw4Oxp9//on169fD3t4e165dw8yZM2FlZYUhQ4YI7UQikUItn077WEEb5UsQGxvL7VCMtHV7pybJYGQk1XQZKpFKM2BoYAgjI2NNl6IS03KmsKtYXtNlqExb921txG2tvHzDTZs2bWBsbAyZTIb58+fDw8MD9evXl2sjEolQunRpfPPNN2jcuLFai5s9ezbGjx+PXr16AQDq1q2Lx48f45dffsGQIUMgkUgAAM+fP0f16tWF5V6+fKnQm0NERERfhnzDTbNmzdCsWTMAgFQqhZubG+rWrVsshQHA27dvoa+vLzdNX18f2dnZAABra2tIJBKEhoaiUaNGAICMjAycO3cO8+fPL7Y6iYiIqOTIN9x8bObMmUVZR646d+6MX3/9FdbW1rC3t0dMTAwCAwPRr18/AB96jby9vbFs2TLY2dnB1tYWAQEBKFOmDHr37l3s9RIREZHm5RluduzYAQDo168fRCKR8Lgg/fv3V09lAJYsWYJFixZh6tSpePnyJSQSCYYOHYrp06cLbb7//nukp6fDx8cHKSkpcHR0RHBwMMqVK6e2OoiIiEh7iFJSUnI9x7tChQoQiUSIj49HqVKlUKFChYJXJhIhKSlJ7UWS+nFgWvHS1u19KUmGaWdfaLoMlUilGZjbsjr8LuZ90dGSKKClORpXzPtEiJJKW/dtbcRtrbw8e26uXr0KAChVqpTcYyIiIqKSLM9wY2Vlle9jIiIiopKoRN84k4iIiEhVSp8tBQCnT5/G5s2b8eDBAyQnJyvckkEkEiE6Olqd9RERERGpROlwExQUhFmzZqFSpUpo3Lgxvv7666Ksi4iIiKhQlA43gYGBaNGiBfbu3SsMMiYiIiIqaZQec5OYmAgPDw8GGyIiIirRlA43Dg4OePToUVHWQkRERPTZlA43ixYtwvbt2xEWFlaU9RARERF9FqXH3Pj5+cHU1BQ9evSAjY0NLC0tFW5qKRKJsGvXLrUXSURERKQspcPNrVu3IBKJUL16dUilUty9e1ehjUikfZcOJyIiIt2idLi5du1aUdZBREREpBa8QjERERHpFKV7bh4/fqxUO0tLy0IXQ0RERPS5lA43DRo0UGpMTVJS0mcVRERERPQ5lA43q1evVgg3WVlZePjwIf78809YWFhgxIgRai+QiIiISBVKh5uBAwfmOW/SpElo164dXr9+rZaiiIiIiApLLQOKy5Yti4EDB+K3335Tx+qIiIiICk1tZ0sZGhri2bNn6lodERERUaGoJdxcu3YNa9asQe3atdWxOiIiIqJC++yzpVJTU5GWloayZcsiMDBQrcURERERqUrpcNOiRQuFcCMSiSAWi1GrVi306tULYrFY3fURERERqUTpcBMUFFSUdRARERGpBW+/QERERDqF4YaIiIh0CsMNERER6RSGGyIiItIpDDdERESkU5QKNxkZGfD398epU6eKuh4iIiKiz6JUuDE2NsYvv/yCuLi4oq6HiIiI6LMofZ2b+vXr4969e0VZC5FOeyYV4cmbbE2XoZJ07SqXiAiACuFm9uzZGDp0KJo3bw4XF5eirIlIJz15k41pZ19ougyV+DappOkSiIhUpnS4WblyJcRiMfr374+qVauiRo0aMDExkWsjEomwa9cutRdJREREpCylw82tW7cgEolQvXp1AMCjR48U2uR2Y00iIiqYSE8Pl5K07zigcVn27lHJo3S4uXbtWlHWQUT0RXuZngW/iy81XYbKZjc00nQJRAp4nRsiIiLSKSqFm6ysLOzatQvjx4+Hp6cn/v33XwBASkoK9u3bh/j4+CIpkoiIiEhZSoeb1NRUdOrUCaNHj8b+/ftx/PhxJCYmAgDKlSuHWbNmYd26dUVWKBEREZEylA438+bNw61bt7B7925ER0dDJpMJ8/T19eHm5objx48XSZFEREREylI63Bw6dAijRo1Chw4dcj0rysbGBo8fP1ZrcURERESqUjrcpKSkoGbNmnnOl8lkePfunVqKIiIiIiospcONlZUVbty4kef8iIgI2NraqqWoj8XHx2PMmDGwsbGBRCKBk5MTzp49K8yXyWTw8/ODvb09KleujG7duuHmzZtqr4OIiIi0g9Lhpk+fPtiyZQsiIiKEaTmHp9auXYuDBw9iwIABai0uJSUFLi4ukMlk2LVrF6KiorBkyRKYm5sLbVasWIHAwEDhruXm5ubo2bMnXr16pdZaiIiISDsofRG/yZMn49KlS+jevTtsbW0hEokwc+ZMJCUlISEhAd26dcPo0aPVWtzKlStRuXJlrF27VphWo0YN4f9lMhmCgoIwadIkuLu7AwCCgoJgZ2eHPXv2wMvLS631EBERUcmndM+NoaEhdu3ahTVr1sDW1hZfffUVMjMz0bBhQ6xZswZ//PGH2m+/cOjQITg6OsLLywu2trZo2bIl1q1bJ5yp9fDhQyQkJKBdu3bCMiYmJnB2dkZUVJRaayEiIiLtoHTPTY4+ffqgT58+RVGLggcPHmDDhg0YO3YsJk2ahGvXrmHGjBkAgFGjRiEhIQEA5A5T5Tx+9uxZnuuNjY0tuqK1CLdD8Up7lQapNEPTZajkfeZ7rasZ0M66tbHmD4z4XVKMuK0/sLOzy3e+yuEGAP7991/htG9LS0vUrVu3SG6amZ2djW+++QZz5swBADRs2BD37t3D+vXrMWrUKKHdp88tk8nyraegjfIliI2N5XYoRrGxsTAtZwojI6mmS1GJoYEhjIyMNV2GSqTSDK2sWxtrzsHvkuLB723lqRRu9u7dizlz5uDp06fCoSGRSISqVatizpw5au/RkUgkqF27tty0r776CnFxccJ8AHj+/Llwt3IAePnypUJvDhEREX0ZlB5zs23bNowYMQKlS5fGvHnzsH37dmzbtg3z5s2DiYkJRo8ejW3btqm1uGbNmuHu3bty0+7evQtLS0sAgLW1NSQSCUJDQ4X5GRkZOHfuHJycnNRaCxEREWkHpXtuli9fDkdHRxw8eBDGxvJdpyNHjkTXrl2xfPlyDBw4UG3FjR07Fp06dUJAQAA8PDwQExODdevW4aeffgLwodfI29sby5Ytg52dHWxtbREQEIAyZcqgd+/eaquDiIiItIfSPTdPnjxBnz59FIINABgbG8PT0xNPnz5Va3GNGjXCtm3bsG/fPjRv3hwLFizADz/8gBEjRghtvv/+e4wdOxY+Pj5o27Yt4uPjERwcjHLlyqm1FiIiItIOSvfc2Nvb53sG0tOnTxXGx6iDi4sLXFxc8pwvEong6+sLX19ftT83ERERaR+le27mz5+PzZs3Y9++fQrz9u7diy1btmDBggVqLY6IiIhIVUr33KxatQpmZmYYPnw4Zs6ciZo1a0IkEuHevXt48eIFbGxssHLlSqxcuVJYRiQSYdeuXUVSOBEREVFulA43t27dgkgkEk65zhlfY2RkhOrVq0MqleL27dtyyxTFtW+IiIiI8qN0uLl27VpR1kFERESkFkqPuSEiIiLSBgw3REREpFMYboiIiEinMNwQERGRTmG4ISIiIp3CcENEREQ6Relw07BhQxw+fDjP+UeOHEHDhg3VUhQRERFRYSkdbh49eoQ3b97kOf/Nmzd4/PixWooiIiIiKiyVDkvld8Xhu3fv8k7cREREpHH5XqF4+/bt2LFjh/A4ICAAmzdvVmiXkpKCGzdu5Hv3biIiIqLikG+4efPmDRISEoTHqampyM7OlmsjEolQunRpDB06FDNnziyaKomIiIiUlG+4GTlyJEaOHAkAaNCgARYvXoyuXbsWS2FEREREhaH0jTNjYmKKsg4iIiIitVA63OR49eoV4uLikJycDJlMpjC/RYsWaimMiIiIqDCUDjfJycmYMWMG9u3bh6ysLIX5MpkMIpEISUlJai2QiIiISBVKh5vJkyfj4MGDGDlyJFq0aAGxWFyEZREREREVjtLh5sSJExg9ejQWLVpUlPUQERERfRalL+JXqlQp2NjYFGUtRERERJ9N6XDj7u6O48ePF2UtRERERJ9N6XAzYcIExMfHY8yYMbh48SLi4+Px4sULhX9EREREmqT0mBtHR0eIRCJER0dj165debbj2VJERESkSUqHm+nTp+d740wiIiKikkDpcOPr61uUdRARERGphdJjbj6WlZWFpKQkZGZmqrseIiIios+iUrj5559/0KNHD1StWhW2traIiIgAACQmJqJv3744c+ZMkRRJREREpCylw82FCxfQtWtX3L9/H/369ZO7r5SZmRlev36NP/74o0iKJCIiIlKW0uFmwYIFsLGxQVRUFGbPnq0wv1WrVrh06ZJaiyMiIiJSldLh5p9//sGgQYNgbGyc61lT1apVQ0JCglqLIyIiIlKV0uFGT08Penp5N09ISICJiYlaiiIiIiIqLKXDjYODA44cOZLrvHfv3mH37t1o2rSp2gojIiIiKgylw82UKVMQFhaG8ePH49q1awCA+Ph4nDhxAt27d8f9+/cxderUIiuUiIiISBlKX8Svbdu2WLt2LXx8fLB9+3YAgLe3N2QyGcqXL4/169ejSZMmRVYoERERkTKUDjcA0Lt3b3Tt2hWhoaH477//kJ2djZo1a6J9+/YoW7ZsUdVIREREpDSVwg0AlC5dGt26dSuKWoiIiIg+m9Jjbg4fPgwfH5885/v4+OQ54JiIiIiouCgdblatWoW3b9/mOT8jIwMrVqxQS1FEREREhaV0uLlx4wYcHBzynN+wYUPcunVLHTURERERFZrS4SYzMxPp6el5zk9PT4dUKlVLUURERESFpXS4qVOnDkJCQpCdna0wLzs7GyEhIbC3t1drcZ9atmwZxGKx3NgfmUwGPz8/2Nvbo3LlyujWrRtu3rxZpHUQERFRyaV0uBkzZgwuX76M/v37Izo6GlKpFFKpFNHR0RgwYAAuX76M0aNHF1mhFy9exObNm1G3bl256StWrEBgYCD8/f1x6tQpmJubo2fPnnj16lWR1UJEREQll9Kngvfq1Qv379+Hn58fjh8/DgAQiUSQyWQQiUSYMWMGPD09i6TI1NRUjBw5EqtWrcKSJUuE6TKZDEFBQZg0aRLc3d0BAEFBQbCzs8OePXvg5eVVJPUQERFRyaXSdW6mTZuG3r1748CBA3jw4AFkMhlq1qwJNzc31KhRo4hKhBBevv32W7lw8/DhQyQkJKBdu3bCNBMTEzg7OyMqKorhhoiI6AukVLhJT09H37594enpiUGDBmHChAlFXZdg8+bNuHfvHtauXaswLyEhAQBgbm4uN93c3BzPnj3Lc52xsbHqLVJLcTsUr7RXaZBKMzRdhkreZ77XupoB7axbG2v+wIjfJcWI2/oDOzu7fOcrFW5MTExw9epV9O7dWy1FKSs2Nhbz58/H33//jVKlSuXZTiQSyT3OOVSWl4I2ypcgNjaW26EYxcbGwrScKYyMtOuMQkMDQxgZGWu6DJVIpRlaWbc21pyD3yXFg9/bylN6QHHLli0RGRlZlLUouHDhAhITE9G8eXOYmZnBzMwMERERWL9+PczMzFCxYkUAwPPnz+WWe/nypUJvDhEREX0ZlA43/v7++Oeff/DTTz/hwYMHuZ4Srm7dunVDZGQkwsPDhX/ffPMNevXqhfDwcNja2kIikSA0NFRYJiMjA+fOnYOTk1OR10dEREQlj9IDips0aQKZTIbAwEAEBgZCT08PhoaGcm1EIhGePn2qtuLEYjHEYrHctNKlS6NChQqoU6cOAMDb2xvLli2DnZ0dbG1tERAQgDJlyhT7ITQiIiIqGZQONz179sx3HIumfP/990hPT4ePjw9SUlLg6OiI4OBglCtXTtOlERERkQYoHW6CgoKKsg6lHTp0SO6xSCSCr68vfH19NVQRERERlSRKj7khIiIi0gYqhZtHjx5h4sSJcHBwgKWlJc6ePQsASExMxNSpUxEdHV0UNRIREREpTenDUrdv30bnzp2RnZ2Nxo0b49GjR8jKygIAmJmZ4eLFi5BKpVi9enWRFUtERERUEKXDzZw5c1CuXDmcOHEC+vr6sLW1lZvfqVMn/PXXX+quj4iIiEglSh+WioyMxIgRI2BhYZHrWVOWlpb53vKAiIiIqDgoHW4yMzNRpkyZPOcnJydDX19fLUURERERFZbS4aZOnToIDw/PdZ5MJsOBAwfg4OCgrrqIiIiICkXpcOPt7Y39+/djyZIlSEpKAgBkZ2fjzp07+O6773DlypVivVs4ERERUW6UHlDcq1cvPH78GIsWLcLixYuFaQCgr6+PhQsXomPHjkVTJREREZGSlA43ADBp0iT07t0bISEhuHfvHrKzs1GzZk10794d1tbWRVUjERERkdIKDDdSqRSHDx/GgwcPULFiRbi4uGDs2LHFURsRERGRyvINNwkJCejatSvu378PmUwGAChTpgx27tyJFi1aFEuBRERERKrId0DxwoUL8eDBA4wdOxY7d+6En58fjIyMMH369OKqj4iIiEgl+fbcnDp1Cv3798fChQuFaRYWFhgxYgSePHmCatWqFXmBRERERKrIt+cmISEBTk5OctOaNWsGmUyGuLi4Ii2MiIiIqDDyDTdZWVkwNjaWm5bzOCMjo+iqIiIiIiqkAs+WevDgAS5fviw8TktLAwDExsaibNmyCu0dHR3VWB4RERGRagoMN35+fvDz81OY/umgYplMBpFIJFy9mIiIiEgT8g03gYGBxVUHERERkVrkG24GDBhQXHUQERERqYXSN84kIiIi0gYMN0RERKRTGG6IiIhIpzDcEBERkU5huCEiIiKdwnBDREREOoXhhoiIiHQKww0RERHpFIYbIiIi0ikMN0RERKRTGG6IiIhIpzDcEBERkU5huCEiIiKdwnBDREREOoXhhoiIiHQKww0RERHpFANNF0CkqmdSEZ68ydZ0GSpJK2UBQ+0qmUgpxsYmuJQk03QZKqtWRg9VjLSvblIOww1pnSdvsjHt7AtNl6ESqTQDc1uW03QZRGqX9E6GgCjt+jwCQEBLc1QxEmm6DCoiPCxFREREOoXhhoiIiHQKww0RERHplBIdbpYvX462bdvC0tISNjY28PT0xI0bN+TayGQy+Pn5wd7eHpUrV0a3bt1w8+ZNDVVMREREmlaiw83Zs2cxfPhwHD16FCEhITAwMECPHj2QnJwstFmxYgUCAwPh7++PU6dOwdzcHD179sSrV680WDkRERFpSok+Wyo4OFju8dq1a2FlZYXz58+jS5cukMlkCAoKwqRJk+Du7g4ACAoKgp2dHfbs2QMvLy9NlE1EREQaVKJ7bj71+vVrZGdnQywWAwAePnyIhIQEtGvXTmhjYmICZ2dnREVFaahKIiIi0qQS3XPzqZkzZ6J+/fpo2rQpACAhIQEAYG5uLtfO3Nwcz549y3M9sbGxRVekFtHW7ZBWygJSaYamy1DZ+8z3Wle3NtYMaGfd2ljzB+W1su60V2mITXyu6TJUpq3f2+pmZ2eX73ytCTc//PADzp8/jyNHjkBfX19unkgkfyEmmUymMO1jBW2UL0FsbKzWbofUJBmMjKSaLkMlUmkGDA0MYWRkrOlSVKKNNXNbFz9trNu0nCnsKpbXdBkq0ebv7eKmFYelfH19sXfvXoSEhKBGjRrCdIlEAgB4/lw+fb98+VKhN4eIiIi+DCU+3MyYMQN79uxBSEgIvvrqK7l51tbWkEgkCA0NFaZlZGTg3LlzcHJyKu5SiYiIqAQo0Yelpk2bhp07d2Lr1q0Qi8XCGJsyZcqgbNmyEIlE8Pb2xrJly2BnZwdbW1sEBASgTJky6N27t4arJyIiIk0o0eFm/fr1ACCc5p1jxowZ8PX1BQB8//33SE9Ph4+PD1JSUuDo6Ijg4GCUK8ebFBIREX2JSnS4SUlJKbCNSCSCr6+vEHaIiIjoy1bix9wQERERqYLhhoiIiHQKww0RERHpFIYbIiIi0ikMN0RERKRTGG6IiIhIpzDcEBERkU5huCEiIiKdwnBDREREOoXhhoiIiHQKww0RERHpFIYbIiIi0ikMN0RERKRTGG6IiIhIpzDcEBERkU5huCEiIiKdwnBDREREOoXhhoiIiHQKww0RERHpFIYbIiIi0ikMN0RERKRTGG6IiIhIpzDcEBERkU5huCEiIiKdwnBDREREOoXhhoiIiHQKww0RERHpFIYbIiIi0ikMN0RERKRTDDRdABERUXET6enhUlK2pstQiXHZSpouQWsw3BAR0RfnZXoW/C6+1HQZKpnd0EjTJWgNHpYiIiIincKemy/cM6kIT95oV9dsunaVS0RExYzh5gv35E02pp19oekyVOLbhMediYgobzwsRURERDqF4YaIiIh0CsMNERER6RSGGyIiItIpDDdERESkUxhuiIiISKcw3BAREZFO0Zlws379ejRo0AASiQTffvstIiMjNV0SERERaYBOhJvg4GDMnDkTU6dORVhYGJo2bYo+ffrg8ePHmi6N8iHKzkb7i0fxw6bZaH/xKETZvPQwERF9Pp0IN4GBgRgwYACGDh2K2rVrY+nSpZBIJNi4caOmS6M8iLKzsXL5KCxcMw19Qv/EwjXTsHL5KAYcIiL6bFp/+4V3794hOjoaEyZMkJverl07REVFFVsdWnePJjNbXEqSaew+Te0uH0fjm1Ewfi8FABi/l8Lx5gW0vXwcp5q4aKaoImRkZKzpEr4Y3NbFy9DAUNMlfDHKlxfjUpIW/c78n2pl9FDFSFasz6n14SYxMRFZWVkwNzeXm25ubo7nz58XWx1VjGSoYiQqtudTp9PdLYr9OY1D/4HR/wWbHCbvM/DzmyvI6D64wOVdNFCzOmhj3dpYM6CddWtjzQDrLj4yANr4O1O8wQbQkcNSACASyb/hMplMYRqVHJlt2kBmLP8XtszYGJlt2mimICIi0hlaH27MzMygr6+v0Evz8uVLhd4cKjky3dyQ2bKlEHBkxsbIbNkSma6uGq6MiIi0ndaHm1KlSsHBwQGhoaFy00NDQ+Hk5KShqqhAenp4u2sX3q5bB+l33+HtunV4u2sXoKf1uyQREWmY1o+5AYBx48Zh9OjRcHR0hJOTEzZu3Ij4+Hh4eXlpujTKj54eMrt3R2b37pquhIiIdIhO/Jns4eEBPz8/LF26FK1atcL58+exa9cuWFlZabq0EuN///sfXF1dYWVlBbFYjIcPHyq13P79++Hk5AQLCws4OTnhwIEDRVyp9pNKpfDx8UGtWrVQtWpV9OvXD0+ePMl3mW3btkEsFiv8y8jIKKaqtYeqF+y8fv06unbtisqVK+Prr7+Gv78/ZLLiH+CojVTZ1g8fPsx1Hz5x4kQxVqydIiIi0K9fP3z99dcQi8XYtm1bgctwv86fToQbABgxYgSuXbuG58+f48yZM2jRooWmSypR3r59i3bt2mHmzJlKL3PhwgV899136NOnD8LDw9GnTx8MGzYMly5dKsJKtZ+vry8OHDiADRs24PDhw3j16hU8PT2RlZWV73KlS5fG7du35f4ZG/O05o+pesHOtLQ09OzZExYWFjh16hQWL16MVatWYfXq1cVcufYp7MVR9+7dK7cPt27dupgq1l5v3rxBnTp1sHjxYpiYmBTYnvt1wUQpKSmMel+QK1euoG3btrh69Sqsra3zbevl5YXk5GT89ddfwjR3d3dUqlQJGzZsKOJKtVNqaipsbW0RGBiIvn37AgDi4uJQv3597NmzB+3bt891uW3btmH69OkF9vB86dq3b4+6deti5cqVwrRGjRrB3d0dc+bMUWi/YcMGzJ07F3fu3BF+NJYuXYqNGzfixo0bPKMyH6pu64cPH6Jhw4YIDQ3FN998U5yl6pRq1aphyZIlGDhwYJ5tuF8XTGd6bkj9Ll68iHbt2slNa9++fbFeHFHbREdH4/3793LbrXr16qhdu3aB2y09PR316tVDnTp14OnpiatXrxZ1uVol54Kdn+6T+V2w88KFC2jevLncX8Pt27fHs2fPlD40+yUqzLbOMXjwYNja2sLFxQX79+8vyjK/WNyvC8ZwQ3lKSEjQ+MURtc3z58+hr68PMzMzuekFbTc7OzusXr0a27dvx/r162FkZITOnTvjv//+K+qStUZhLtj5/PnzXNvnzKPcFWZbly1bFgsWLMCmTZuwe/dutG7dGl5eXti5c2dxlPxF4X5dMJ04W+pLtXDhQgQEBOTb5sCBA2jVqlWhn4MXR/xA2W2dl4K2W9OmTdG0aVPhsZOTE1q1aoW1a9diyZIlqhesw1TdJ3Nrn9t0UqTKtjYzM5O7Dc4333yDpKQkrFixAp6enkVa55eI+3X+GG60mLe3tzCuIy/Vq1cv9PolEgkvjvh/lN3WFy9eRFZWFhITE1GpUiVh3suXL+Hs7Kz08+nr68PBwQH37t0rdM26pjAX7LSwsMi1PYAvcj9Wlroujuro6KjUmT+kGu7XBWO40WJmZmYKhz/UqUmTJggNDcXEiROFaV/qxRGV3dYODg4wNDREaGgo+vTpAwB48uQJbt++rdJ2k8lkuH79OurVq1fomnXNxxfs7NGjhzA9NDQU3fO4VlLTpk0xd+5cZGRkCGeehYaGokqVKgUOqP+SFWZb5+batWuQSCRFUOGXjft1wTjm5guRkJCAmJgY3L17FwBw+/ZtxMTEIDk5WWjTvXt3zJs3T3g8ZswYhIWFYfny5bhz5w6WL1+O8PBweHt7F3v92qJ8+fIYPHgwZs+ejdOnT+Pq1asYPXo06tatizYf3Tfr0229ePFinDx5Eg8ePEBMTAzGjx+P69ev47vvvtPAqyi5xo0bh+3bt2PLli24ffs2ZsyYIXfBznnz5sn9+Pbu3RsmJiYYO3Ysbty4gZCQEPz6668YO3Ysu+8LoOq23r59O3bv3o3bt28jNjYWq1atwvr16zFq1ChNvQSt8fr1a8TExCAmJgbZ2dmIi4tDTEyMcNo992vVsefmC7Fx40b4+/sLj3MOsQQGBgqnHN6/fx/VqlUT2uRc7XnhwoXw8/NDzZo1sXHjRjRu3Lh4i9cyP//8M/T19eHl5YWMjAy0bt0aa9asgb6+vtDm022dmpqK77//Hs+fP4epqSkaNGiAw4cPw9HRURMvocTy8PBAUlISli5dioSEBHz99ddyF+yMj4/H/fv3hfbly5fHvn37MG3aNLRt2xZisRjjxo3D+PHjNfUStIaq2xoAAgIC8PjxY+jr68PGxgarV6/meBslXLlyBW5ubsJjPz8/+Pn5oX///ggKCuJ+XQi8zg0RERHpFB6WIiIiIp3CcENEREQ6heGGiIiIdArDDREREekUhhsiIiLSKQw3REREpFMYbkhn+fn5QSwWK0z/7bff8M0338DMzAz169cvcDoREWkXhhvSCtu2bYNYLBb+SSQS2Nvbw8PDA2vWrMGrV6+UWs/Zs2fxww8/oGHDhli1ahX8/Pzynf4lq1+/vtw2t7CwgIODA3x8fJCQkFCodb5+/Rp+fn4IDw9Xc7Ul06lTpzBo0CDY29vD3NwcVlZWaN++PRYtWoSnT59qurxCCw8PF/aLrVu35tpm2LBhwmeVqLjxCsWkVWbOnImaNWvi/fv3eP78Oc6ePQtfX18EBgZix44dcvdi8vHxweTJk+WWz/lR/fXXX+V6dfKa/qWrW7eucG8xqVSKa9euYfPmzQgLC0NkZKTcVZeV8ebNG+FK2Z9zt/qSTiaTYerUqdi4cSPs7e0xZMgQWFpaIj09HdHR0Vi7di02bdok3A5FWxkbG2P37t0YNGiQ3PS0tDQcOXIExsbGwt2qiYoTww1plfbt26NJkybC4ylTpuDMmTPo168f+vfvjwsXLsDExAQAYGBgAAMD+V085865nwaYvKZ/jrdv36J06dJqW58mVK5cWeHy+eXLl0dAQACuXbsGBwcHzRRWwgUGBmLjxo0YNWoUFi9eDD09+U7yxYsXY+XKlQWuJz09XdifS6JOnTrh4MGDePbsGapUqSJMDwkJQXZ2Ntq1a4fQ0FANVkhfKh6WIq337bffwsfHB48fP8auXbuE6Z+OuRGLxdiwYYPw/2KxWGiT2/QcoaGhcHV1RfXq1VG1alW4uroiKipKroac9dy6dQtjxoxBzZo10axZs0Kt47///sPkyZNRs2ZNVKtWDUOHDkVSUpLC6w4NDYWbmxssLS1RvXp1fPvtt9iyZYtcmytXrsDT0xNWVlaoXLky2rVrhyNHjqi4heXlHGYwNDSUmx4fH4/vv/8e9vb2sLCwQKNGjbBixQrhL/eHDx+idu3aAAB/f39hW3t7e+Pff/+FWCzGvn37hPXdv38fYrFY4c7okydPxldffVWo15mWloYff/wR9evXh4WFBerVq4e5c+dCKpXKtROLxZg8eTKOHz+OVq1aQSKRoFGjRtizZ0+B2yc9PR3Lli1D7dq14efnpxBsAMDU1BQ//vij3LT69eujV69eCAsLQ4cOHSCRSPDrr78CAJKTkzFlyhTUrl0bFhYWaNq0KVavXi3XK/Lw4UOIxWJs27ZN4fnq168vd8PbnMO8YWFh8PHxQa1atVCtWjUMGTIE8fHxBb7GHC4uLihbtqzCdtm9ezc6dOiAChUq5LqcMp+HR48eYerUqWjSpAmqVKkCKysreHp64ubNm3Ltcg6R7dmzB6tXr0b9+vUhkUjQsWNHXL16Va7t8+fPMWHCBNStWxcWFhawt7eHp6cnrl+/rvRrJu3AcEM6Iad34dSpU3m2Wbt2LVq3bi38/9q1a+Hm5pbndADYs2cPevXqBX19fcyaNQuzZs1CUlISunfvjkuXLik8h5eXF5KTkzFr1iyMGTOmUOsYPnw4nj59ilmzZmHIkCE4ePAgpk+fLtfmzz//hIeHB+Lj4zFhwgTMmzcPjo6OOHr0qNDm7Nmz6Ny5M54/fw4fHx/MmzcPpUqVQv/+/RESEqLUdn3//j0SExORmJiIp0+f4tixY1ixYgUcHBxQp04dod2LFy/QoUMHHD16FEOHDoW/vz8aN26MOXPmwNfXFwBQqVIlLF26FADg6uoqbGsvLy/UrVsXYrEYERERwjojIiKgp6eHuLg4PHz4UJgeGRmJ5s2bq/w609PT4erqij/++AMeHh5YsmQJXFxcsHr1auFO1x+7ePEixo0bh65du2LBggUoXbo0Ro0ahdu3b+e7zaKiopCcnIzevXurfNju3r17GDJkCJydneHv748mTZpAKpXCzc0NmzdvRvfu3bFo0SJYW1vjxx9/xA8//KDS+j81c+ZMREdHY/r06Rg2bBj+/vtveHh44N27d0otb2xsDDc3N+zevVuYFh8fj/DwcOHmvJ9S9vNw5coVREREwM3NDX5+fvD29saVK1fQtWvXXMd8rV69Gjt27MCoUaMwc+ZM3L17FwMHDsT79++FNkOHDsX+/fvRv39/BAQEYPTo0cjOztb6w4OkiIelSCdUq1YNpqamCncp/pinpyfOnz+PsLAwuUMt9erVy3X6mzdvMG3aNHh6eiIoKEiY7uXlhWbNmmH+/PkKIcHW1hZ//PHHZ63jq6++wrp164THMpkMv//+O5YtW4by5csjLS0N06dPR926dXH06FGUKVNGrm3OfydPnoymTZti//79Qu/ByJEj4eLigtmzZ6N79+75b1QAYWFhsLGxkZvm6OiIP//8EyKRSJi2cOFCSKVSREREwMLCQniNlStXxurVq+Ht7Q1ra2t0794dPj4+qFu3rsLhrmbNmiEyMlJ4fO7cObRr1w5RUVE4d+4crK2tkZiYiDt37uC7775T+XX+9ttviI2NxenTp4UeJAD4+uuvMW3aNERGRsLZ2VmYfuvWLURERAhte/TogXr16mHr1q1YsGBBntvs1q1bwno/JpPJFHrgTE1N5XrA7t+/j+3bt6Nr167CtHXr1uHff//FypUrMWTIEADAiBEjMHjwYKxZswYjRoxQeI9UcfDgQRgZGQEA7O3tMWHCBGzfvh3Dhg1Tavm+ffti27ZtuHXrFuzt7bF7926ULVsWnTt3lgvbgGqfh44dO8Ld3V1ueU9PTzRv3hx//PEHpk2bJjcvLS0NkZGRMDY2BgDY2dlh0KBBOHXqFFxcXJCamopz585hwYIFmDBhgrDcp+PySDew54Z0RtmyZfH69Wu1rS80NBQpKSno27ev0HuRmJiI9PR0tGnTBufOnZP7qxD40Oui7nW0aNECWVlZiIuLE9aZlpaGqVOnygUbAELguHbtGmJjY9G3b18kJycLz5ucnIwOHTrgwYMHePToUYHb4JtvvsFff/2Fv/76Czt37sTs2bPx33//YdCgQUhPTwfw4Ud7//79cHFxgb6+vtzrbN++PbKzs+V6ZPLi7OyMmzdvIjk5GcCHHppWrVqhadOmQuiJiIiATCYTem5UeZ379u2Dk5MTKlWqJFdjmzZtAHwIch9r1aqVXAiysLCAnZ0dHjx4kO/ryDlzr1y5cnLTk5KSYGNjI/fv/Pnzcm2qVasmF2wA4OjRozAzM8PAgQOFaSKRCBMnToRMJsOxY8fyrSc/Xl5eQrABgP79+6N8+fIqrbNVq1aoUqWK0Huze/duuLq6CiHjY6p8Hj4er/b27VskJSWhfPnysLGxQXR0tMK6Bw4cKPecLVu2BADh/TI2NoahoSHOnj0r7GOku9hzQzrj9evXqFSpktrW999//wEAevbsmWeb1NRUueesUaPGZ6/D0tJSbn7OuKGcL+Sc3qmPDwvlVfuECRPk/kr92MuXL2FlZZXnOgCgYsWKwo8/8GGMhZ2dHQYPHowtW7Zg9OjRePnyJVJSUrB169Y8TwvOGbCdH2dnZ8hkMkRGRsLR0RH379+Hs7MzMjMz8eeffwL4EHhMTU2FcTiqvM7//vsP//77b569HJ/W+On7AHx4Lwr6YcwJNZ9ensDU1BR//fUXgA+H0gICAhSWtba2Vpj26NEj2NjYKBziygleyoTUvHy6LQwMDGBtbY3Hjx8rvQ49PT14eHhg9+7d6NOnD2JiYjB//vxc26ryecjIyMDPP/+MXbt2KYwDMjMzU1iuoM+NkZER5syZgzlz5sDOzg6NGzdGx44d0bdv31zfa9JuDDekE548eYK0tDTUqlVLbevMzs4G8OFwRtWqVXNtY2pqKvf40zNbCrOOvMZpfHzICYDcYaG8ap87d26eZzTZ2trmuXx+csYnRUZGCmMWAKB3794KpwTnUOZ9cXBwQJkyZRAZGQmpVIrSpUvDwcEB79+/x4IFC/DixQthvE3O4SdVXmd2djZat26NKVOm5Nru0/enoPchL/b29gCAGzduwNXVVZhuaGgoBMXExMRcl/2cM6OU2R+UWaYwp2736dMHgYGB8PHxQeXKlfM8zV+Vz8PMmTOxZcsWjBo1Cs2aNYOpqSn09PTg6+ub6+tR5v0aP348XF1dcfjwYZw+fRpLly7F8uXLsX37dnz77bcqvWYq2RhuSCfs3LkTANCuXTu1rbNmzZoAPgyE/bj3orjX8amcoHDjxg2Fs4Y+fd6yZcuq7XlzZGZmAvgwfgL48NpMTU2RmZlZ4HPl9wNsYGCAxo0bIzIyEu/evUOTJk1gaGgIR0dHGBkZ4ciRI7h+/To8PDyEZVR5nTVr1sTr16/Vvj0+5eTkhAoVKmDPnj2YOnWqyoOKP2VlZYWrV68iKytLbl137twR5gMQzkxKTU2VW14qleZ5BtTdu3fRtm1b4XFmZiYePXqEFi1aqFSjg4MDvvrqK4SHh2Ps2LF5vmZVPg/BwcHo168fFi9eLDc9JSUFFStWVKm+j9WoUQNjx47F2LFjERcXh9atW+OXX35huNExHHNDWu/MmTNYunQprK2t8zxDozDat28vXNPl01OFAeUOtahjHZ9q27YtTE1NsXz5crx9+1ZuXs5fqQ4ODrCxscGqVasUfuwK+7w5csZj5Bwa0tfXR/fu3XHw4MFcx0KkpqYqjKNISUnJdd3Ozs6IiYnBiRMnhMG9RkZGaNSoEVauXImsrCy5Qb+qvE4PDw/8888/OHz4sEK79PR0tY3XMjExwdSpU3Hnzp08exlU6R1xcXHBy5cvsWPHDrnlV61aBZFIhE6dOgH4cDisUqVKCld/3rhxI7KysnJd96ZNm+T2yx07diA1NRUdO3ZUur4cixYtwowZMxTGjH1Mlc+Dvr6+wnbas2cPnj17pnJtwIdxOznjxHJUr14d5ubmee6PpL3Yc0Na5eTJk7h37x4yMzPx4sULhIWFITQ0FJaWltixY0eugxgLq1y5clixYgWGDx+Oli1bok+fPpBIJHjy5AnCw8NRpkyZAq97oo51fMrU1BR+fn4YP3482rZtiz59+qBixYq4efMmnj17hq1bt0JPTw+rV69Gr1690KxZMwwcOBBWVlaIj4/HxYsX8fjxY4XBrLmJj48XesXevXuH69ev43//+x/MzMwwatQood3cuXMRERGBzp07Y/DgwahTpw5evXqFGzdu4MCBA/jnn38gkUhQtmxZ2NnZITg4GLa2tqhYsSKsra3RuHFjAEDz5s2RlZUljLfJ0aJFCwQEBMDExETu8JMqr3PChAk4duwYBg8ejL59+8LR0RFSqRR3797Fvn37sHv3brkLRH6OcePG4b///sO6desQFhaG7t27w9LSEm/evMHt27exd+9elClTJtexI58aMmQItmzZgkmTJuHatWuwtbXF8ePHcezYMYwZM0Zu3MywYcMQEBCAsWPHokmTJrhy5QrOnDmT7/O4ubmhV69eePToEdatWwd7e3sMGDBA5dfcsWPHAkORKp+HLl264M8//0S5cuVQp04dXLt2DcHBwQrj2pR19+5ddO/eHT169IC9vT2MjIxw7Ngx3L59O9+z30g7MdyQVsnpoi5VqhQqVKiAOnXqwM/PDwMHDlQ4O0UdevTogSpVqmD58uX47bffkJ6eDolEgsaNGwun5RbHOj41cOBAmJub45dffsHy5cuhr68PGxsbjBgxQmjTvHlznDx5EkuWLMH//vc/pKWlwdzcHPXq1ROuPVOQ69evY/To0QA+BAkzMzO4urpi1qxZcmMmKlWqhJMnT2Lp0qU4dOgQ/ve//6F8+fKwtbXFzJkz5S7mFhgYCF9fX/z444+QSqXo37+/EG6aNGmCUqVKAYAwLee15EzLma/q6zQxMUFISAhWrFiB4OBgIWDUqFED3t7esLOzU2qbKEMkEuGXX35Bt27dsHHjRmzevBmJiYkoXbo0bG1tMWbMGAwbNizPcScfMzY2RkhICBYsWIB9+/YhOTkZ1tbWWLBgAcaPHy/Xdtq0aUhKSkJwcDD++usvtGzZEvv37xeu2/SpxYsXIyQkBP7+/pBKpXBxccHSpUvlzqBSN2U/D4sXL4ahoSH27duHrVu3wsHBAXv37sVPP/1UqOetXr06+vTpg7CwMOzZswcikUjo9Rs8eLC6Xh6VEKKUlBTe+IOI6Auybds2jBs3DsePH1dbbxVRScIxN0RERKRTGG6IiIhIpzDcEBERkU7hmBsiIiLSKey5ISIiIp3CcENEREQ6heGGiIiIdArDDREREekUhhsiIiLSKQw3REREpFP+H6JeCWNGJPemAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "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": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0102"
      ]
     },
     "execution_count": 26,
     "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.6.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}