{
 "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>5</th>\n",
       "      <td>138</td>\n",
       "      <td>244</td>\n",
       "      <td>33</td>\n",
       "      <td>62</td>\n",
       "      <td>178</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>132</td>\n",
       "      <td>245</td>\n",
       "      <td>23</td>\n",
       "      <td>65</td>\n",
       "      <td>140</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>120</td>\n",
       "      <td>289</td>\n",
       "      <td>25</td>\n",
       "      <td>62</td>\n",
       "      <td>125</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>143</td>\n",
       "      <td>299</td>\n",
       "      <td>30</td>\n",
       "      <td>66</td>\n",
       "      <td>136</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>140</td>\n",
       "      <td>351</td>\n",
       "      <td>27</td>\n",
       "      <td>68</td>\n",
       "      <td>120</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\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",
       "5           138               244            33               62   \n",
       "6           132               245            23               65   \n",
       "7           120               289            25               62   \n",
       "8           143               299            30               66   \n",
       "9           140               351            27               68   \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",
       "5                        178            False  \n",
       "6                        140            False  \n",
       "7                        125            False  \n",
       "8                        136             True  \n",
       "9                        120            False  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "births = pd.read_csv(path_data + 'baby.csv')\n",
    "births.head(10)"
   ]
  },
  {
   "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>5</th>\n",
       "      <td>False</td>\n",
       "      <td>138</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>False</td>\n",
       "      <td>132</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>False</td>\n",
       "      <td>120</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>True</td>\n",
       "      <td>143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>False</td>\n",
       "      <td>140</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\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",
       "5            False           138\n",
       "6            False           132\n",
       "7            False           120\n",
       "8             True           143\n",
       "9            False           140"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smoking_and_birthweight.head(10)"
   ]
  },
  {
   "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>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>False</td>\n",
       "      <td>113</td>\n",
       "      <td>False</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>True</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>5</th>\n",
       "      <td>False</td>\n",
       "      <td>138</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>False</td>\n",
       "      <td>132</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>False</td>\n",
       "      <td>120</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>True</td>\n",
       "      <td>143</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>False</td>\n",
       "      <td>140</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Maternal Smoker  Birth Weight  Shuffled Label\n",
       "0            False           120            True\n",
       "1            False           113           False\n",
       "2             True           128            True\n",
       "3             True           108            True\n",
       "4            False           136           False\n",
       "5            False           138           False\n",
       "6            False           132           False\n",
       "7            False           120            True\n",
       "8             True           143            True\n",
       "9            False           140           False"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "original_and_shuffled.head(10)"
   ]
  },
  {
   "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>119.581818</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>True</td>\n",
       "      <td>119.276688</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Shuffled Label  Birth Weight\n",
       "0           False    119.581818\n",
       "1            True    119.276688"
      ]
     },
     "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.head(10)"
   ]
  },
  {
   "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": [
       "-0.30512972865913923"
      ]
     },
     "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.3690235690235681"
      ]
     },
     "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([ 0.78235446, -1.2423694 , -0.54122827, ..., -1.41050018,\n",
       "       -2.52660237,  2.47797127])"
      ]
     },
     "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/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABR/ElEQVR4nO3dd1gU1/s28HtpC4K4agDpKGAQYsRusMSKjaIYxRZLrNhiI0qMsQfssaAxtmDsIkYssWMBlK+JGrtixYKoSBHERWHfP3zZn+susIvADnJ/rssr2ZkzM8+c3YWbmTMzotTUVBmIiIiIiARKR9sFEBEREREVhIGViIiIiASNgZWIiIiIBI2BlYiIiIgEjYGViIiIiASNgZWIiIiIBI2BlaiUBAQEQCKR4NSpUwrTJRIJOnfuXGLbDQ4OVrldeievfzZt2qTtUvLVuXNnSCQS3L9/X9ullApV3wl+jgu2adMmSCQSBAcHa7sUAHy/qPgxsNInRSKRKPyrUqUK7O3t0aFDB6xfvx45OTnaLrHYlYXA9aG8mgv65Xrq1KkSD/NCIfRAWrt2bUgkElhbWyMxMVFlm++++05QAeVT+4wJLZASlTY9bRdAVBImTZoEAMjJycHdu3exd+9enDlzBsePH0dYWJiWq1P0v//9D0ZGRiW2/qFDh6Jbt26wsbEpsW1Q+ZCZmYnZs2cjNDRU26WQwPHnDhU3Blb6JAUFBSm8vnLlCtq2bYvdu3cjNjYWHh4eWqpMWc2aNUt0/VWrVkXVqlVLdBtUPjg6OmLLli0YNmwYvvzyS22XQwLGnztU3DgkgMoFNzc3NG3aFADw77//Avi/04EBAQG4fv06+vbtixo1akAikeDixYvyZXfv3g1fX184ODjA3Nwc9erVw/Tp05Genq5yW8ePH0fHjh1hZWUFBwcH9O7dGzdu3Mi3tvxOSebk5GDDhg3o2LEj7O3tYWFhgS+//BKDBw/G+fPnAbw7lTx37lwAwMiRIxWGQ+SdXi5oLNnJkyfRvXt3VK9eHebm5qhTpw4mTZqEZ8+eKbV9fwzu7t270bp1a1haWsLBwQEDBw7Eo0eP8t3H4lTUOi5cuCA/4mNrawtfX1/ExcUVuK2kpCRMnjwZ9erVg4WFBezt7dG1a1ecOHFCqe37p2zj4uLg5+cHe3t7SCQSpKam5rsNiUSCmJgYAECdOnXk71/t2rVVtl+/fj08PDxgYWEBZ2dnjBkzJt/1a1K/OqZPn47c3Fz89NNPai9T0Cl3IY5zXL16NSQSCUJCQlTOT09Ph5WVFdzc3ORDjN4flnPgwAG0a9dO/v0fMGAA7t69q3JdSUlJCAwMRJ06dWBubo7q1aujR48eiI6OVmgXEBCAkSNHAgDmzp2r8D1X1XcXL15Ejx49YGdnB0tLS3Ts2BFnzpxRWUNubi42bNiA9u3bw87ODhYWFvjqq6+waNEiZGdnK7U/deoU/P394ebmBnNzczg5OaFly5aYMmUKZLL/e9J7fu+tussTfYhHWKncu3v3Ljw9PfH555+jZ8+eSEtLQ4UKFQAAEyZMwNq1a2FtbQ0vLy9IJBL8888/+PXXX3Ho0CEcPHgQFStWlK9r9+7dGDhwIPT19dGlSxdYWVnhzJkzaNeuHb744gu1a8rOzkbv3r1x5MgRVKtWDV27dkXlypXx8OFDnDp1Co6Ojqhbty569+4NAIiJiUGnTp0UQk6lSpUK3Mb69esxfvx4GBkZwdfXF9WqVUNcXBxWrVqFffv24e+//4atra3ScmvXrsXff/+NTp06oWnTpvjnn3+wa9cuXLp0CTExMRCLxWrv58fQpI64uDh06dIFUqkU3t7ecHR0xJUrV+Dt7Y0WLVqoXP+VK1fQtWtXPHv2DK1bt0anTp3w4sUL7Nu3D126dMHSpUvx7bffKi33v//9D4sWLYKHhwf69euHxMRE6Orq5rsfkyZNwubNm/HgwQMMHz5c/r6pev+mTZuGY8eOoUOHDmjVqhVOnTqFDRs24NatW9i/f3+x1F+Qdu3aoVWrVoiKisKBAwfQoUMHjZYvC3r27ImZM2fizz//RGBgoNJ7t3XrVrx69QpjxoxRmrdnzx4cOXIE3t7eaN68OS5evIi//voLp06dwqFDh+Do6Chve//+fXTs2BGPHz9G06ZN4efnhydPnuCvv/7CkSNH8Ouvv6Jfv34A3v1hmpaWhv3796Np06Zo1qyZfD12dnYKNVy4cAFLly5F48aN0a9fPzx8+BCRkZHw9fXFyZMn8fnnn8vbvn37Fn379sWBAwfg5OSEbt26QSwWIyYmBjNnzsSJEyewc+dO6Om9iwqHDh2Cv78/KlasiI4dO8La2hqpqam4ffs2Vq1ahRkzZsjbqvKxy1P5xk8GlQvXrl2TH8WqV6+ewrwzZ85g/Pjx+PnnnxWmb9u2DWvXroWXlxdWr16tMM50/vz5mDNnDoKDg/HLL78AADIyMjB27FiIRCLs27cPDRo0kLefOnUqli1bpna9c+fOxZEjR9CyZUts3rxZHqCBd0de846A9unTBwkJCYiJiUHnzp3Rp08ftdafkJCASZMmoUKFCjhy5Ahq1aolnzd79mwsWLAAEyZMwPbt25WWPXbsGE6cOAEXFxf5tMGDByM8PBz79u2Dn5+f2vv5MdStQyaTYdSoUcjKykJYWBh8fX3l7VevXo3AwECldefk5KB///5IS0vDnj17FALCkydP0KZNGwQGBqJ9+/YwNzdXWDYqKgq//vorBgwYoNZ+BAUFITo6Gg8ePEBAQADs7e3zbfvvv//i9OnTsLa2BvAucHh7eyM2Nhb//POP/DP3MfUXZtasWWjRogV+/vlntG3bVvABIzo6Ot8LlRISEpSmVaxYEf7+/lizZg0OHDigdHT4jz/+gJ6enjxMvu/AgQPYtm0b2rdvL5+2bNkyTJ06FYGBgYiIiJBPHzduHB4/fozJkydj8uTJ8umjRo1C27ZtERgYiNatW8PGxgZeXl7ywNqsWTOlIU/vO3jwIFatWgV/f3/5tPXr12PcuHFYtWoVFi1aJJ++ePFiHDhwAEOGDEFISIg8gOfm5mLcuHEICwvDmjVrMHz4cADAhg0bIJPJsGfPHtSpU0dhuy9evCj0s/Cxy1P5xiEB9EkKDg5GcHAwZs+ejSFDhqBVq1bIysqCl5eXfGhAHnNzc/lFWu9bsWIFdHV1sWzZMqWLosaPH4+qVasqBLr9+/cjJSUFfn5+CmEVAH744QeYmpqqVXtOTg7WrFkDsViMJUuWKIRVANDV1UW1atXUWld+tm/fjuzsbAwaNEghrAJAYGAgLC0tcejQITx+/Fhp2WHDhimERADo378/AODcuXMfVZcm1K0jLi4O8fHxaNy4sUJYBYBBgwahRo0aSus+dOgQbt26hUGDBimEPQCoVq0aRo8ejdevX2P37t1Ky37xxRdqh1VN/fDDD/KwCgB6enro27cvAMV9/pj6C/PFF1+gb9++uHnzJtavX1/EPSk9MTExmDt3rsp/W7ZsUbnM4MGDAUBp/86cOYOrV6+iQ4cOsLKyUlquRYsWCmEVeHc638bGBseOHZN/nx49eoRjx47BysoK48ePV2jv5uaG7777DlKpFNu2bdN4f7/66iuFsAoAffv2hZ6ensJnJDc3F7/99hvMzMwQHByscLRYR0cHM2fOhEgkUqhBR+ddZPjwZxIAVKlSpdDaPnZ5Kt/45wx9kvLGdYpEIlSsWBF16tRB9+7dVQaJL774Quk0dlZWFi5evIjKlSvjt99+U7kNAwMDJCYm4sWLF6hSpQr+++8/AFAKxMC7ozZffvml0tg0VW7evIm0tDTUqVOnwKNtHyOvVlWnw8ViMZo0aYJdu3bh4sWLSr+Y3d3dlZbJC1EFjdUsburWUdD7oqOjgyZNmuDOnTsK0/PGtj58+FDl0bm89jdv3lSa9+EfK8VJ3X3+mPrVMWXKFERERCAkJAQ9evQodPiJNk2aNCnfI5KnTp2Ct7e30nQXFxc0a9YMx44dw7179+Dg4ADg/wLsoEGDVK5P1WdMT08PjRs3xsOHD+Xfp7wx8k2aNIGBgYHSMi1btkRoaKj8s6sJVZ8RfX19mJubK3xGbt26heTkZFSvXh3z589XuS4jIyPEx8fLX/fo0QORkZFo06YNunbtiubNm6Nhw4Zq/5z62OWpfGNgpU+SJsFJ1SnRlJQUyGQyvHjxQh5+85ORkYEqVarIL8IyMzNTezuqpKWlAYDKIzjFJa/W/GqysLBQaPc+VUeK847OqHuf27wjLbm5ufm2yZuX17aodRTlfXnx4gUAIDIyEpGRkfnWmJmZqdb6iou6+/wx9avDwsICY8aMwS+//IKFCxdi5syZRVqPkA0ZMgTR0dEICwvDtGnTkJKSgt27d6NGjRpo2bKlymXye+/zPnt5n8WP+f4VJr8zObq6uio/I3fv3i30Z1weLy8v7Ny5E8uWLcOWLVvktwh0dXXFpEmTlM5gFPfyVL4xsFK5JxKJlKbl/dB3dXVFbGysWuvJW0bVFfYA8PTpU7XWk3e0Kr8btBeHvFrzqykpKUmhXUltPyUlJd82eb9QP/boXVHel7xlNmzYAB8fH422p+rzVNo+pn51jR49GmFhYVi1alW+RxyBd/2R3x8yeX+cCVHnzp1hZWWFjRs3IigoCJs3b8br168xYMCAfN/j/L5PeZ+9vPdF29+/99fdoUMHbN26Ve3l2rRpgzZt2iArKwv//vsvjhw5grVr12LAgAFK46VLYnkqvziGlUgFExMTuLq6Ij4+HsnJyWotk3cRQd7FXe97+fKlwq2yClKzZk1UqlQJ165dw4MHDwptr+nRzfdrVXVLHKlUKj+l/OGFEcUl744J+d1q5/15mtxdQZWC3pfc3FyVNTRs2BAAcPr06Y/atrrev9ilOJRG/UZGRvjpp58glUoxY8aMfNtJJBI8fPhQ5by827MJkZ6eHvr3749nz55h7969CAsLg1gsLvDCRlWfsbdv38q/T3n3rs37b1xcnMpbR+Xdduz90/tF+Z4XJO/nzL///quyhsIYGRmhWbNmmD59OmbNmgWZTKZ0p4qSXJ7KHwZWonyMHDkSb968wYgRI1QeCXz58iX++ecf+etOnTpBIpEgIiJCYToAzJs3T+3Te7q6uhgyZAikUinGjh2LrKwshfk5OTl48uSJ/HXezbnzCwWq9OjRAwYGBli7dq3SOMZFixbh8ePH8PT0hKWlpdrr1ISHhweqV6+Oy5cvY8OGDUrzL168iI0bN0JPT0/pAhJNNW7cGM7OzoiLi1O6yGjt2rVK41eBd+9ljRo1sH79+nx/if7333/yo8AfK+89VOcPFHWUVv09e/aEu7s7IiIicOHCBZVtGjZsiIcPH+LQoUMK08PCwgq9D662DRgwAPr6+vjxxx9x8+ZN+Pr6Fngz/JMnT+LgwYMK01auXImHDx+iVatW8mE+1tbWaNOmDR49eoQlS5YotL927RrWrVsHsViMHj16yKcX5XteED09PQwfPhzPnj3DxIkT8erVK6U2ycnJCn9oHz9+XGW7vCPChoaGBW7zY5en8o1DAojy0adPH/z333/4/fff4e7ujjZt2sDOzg5paWlISEhAbGwsWrVqhc2bNwN4d1R2yZIlGDhwIDp37oyuXbvCysoKp0+fxtWrV+Hh4aH28IIffvgB58+fx9GjR1GvXj106NABlStXxuPHj3Hq1Cn07dtXfiHJ119/DR0dHfz2229ISUmRj4sbOnRovqfT7ezsMHfuXIwfPx6tWrVCly5dYGFhgbi4OMTExMDa2hoLFy4shl5UTVdXF7///ju6deuGMWPGYMuWLWjQoAH09fVx8+ZNHDx4EDk5OZg3bx6qV6/+UdsSiURYtmwZunbtioEDByrchzUqKgpt27bFkSNHFJbR19fHxo0b4efnh969e6NBgwaoU6cOjI2N8ejRI1y8eBHx8fE4efJksVzd3KpVK+zatQvff/89fH19YWxsjEqVKmHo0KFFWl9p1S8SiTB79mx4eXmpDP4AMGbMGBw5cgR9+/ZFly5dYGZmhgsXLuDChQto3769UsATEgsLC3h5eWHXrl0AgO+++67A9h07dkSfPn3g4+MDBwcHXLx4EUeOHEGVKlWwYMEChbaLFi1Chw4dMGfOHJw8eRINGzaU34c1KysLS5YsUXisaaNGjWBiYoKIiAgYGBjAxsYGIpEI/v7+SvdiVVdgYCCuXr2KDRs24NChQ2jRogWsra3x/Plz3L17F2fOnMHgwYPlR4R/+uknJCQkoGnTprCzs4OhoSGuXLmCo0ePokqVKvK7dOTnY5en8o2BlagA8+bNg6enJ9auXYvo6GikpKSgUqVKsLKywqBBg9C9e3eF9r6+vti5cyfmzp2L3bt3w8DAAB4eHjh8+DAWL16sdmA1MDDA9u3bERYWhi1btmDHjh14+/YtLCws0LRpU3Ts2FHe1snJCWvXrsWSJUuwceNG+RHZwq7eHjhwIGrUqIFly5Zh3759yMzMhKWlJYYOHYqJEyeW6MVDwLsjb9HR0Vi+fDmioqKwZs0a5OTkwMzMDD4+Phg2bBgaNWpULNtq0qQJ/v77b8yaNQtHjx7F0aNHUb9+fezduxdHjx5VCqzAu/HLMTExWLlyJfbv348tW7ZAJpPBwsICLi4uGD16NJydnYulvr59++LRo0fYvn07QkND8ebNG9ja2hY5sJZm/c2aNUOnTp3yPZLbrFkzbNu2DSEhIYiMjFT4TuzevVvQgRV4997s2rULrq6uaNKkSYFtvby8MGDAACxYsAAHDhyAvr4+fH19MW3aNKXbp9nb2+P48ePytmfOnIGxsTGaNm2KMWPGoHnz5grtK1WqhE2bNiE4OBgRERHIyMgA8O6zXdTAqqenhw0bNmDnzp3YtGkTDh8+LL+I1NbWFuPGjUPPnj3l7SdMmIB9+/bh/Pnz8uFEVlZWCAgIwIgRIxQCtiofuzyVb6LU1FQ+C42IiEiFhQsXYtasWViwYIH8/qwfCg4Oxty5cxEaGqr2wzuISDMcw0pERKRCRkYGVq9eDVNT048eS01EH4dDAoiIiN7z999/4/z58zh8+DCePHmCadOmoWLFitoui6hcY2AlIiJ6T2RkJLZs2QJzc3OMHTsWY8aM0XZJROWe1sawrl69GuvXr5ffxsXFxQUTJ06UP4c5ICBA6TnPDRo0UHlxBBERERF9urR2hNXKygozZsyAo6MjcnNzsWXLFvTp0wfHjx+X3yi8ZcuWWLVqlXwZVc9cJiIiIqJPm9YCa+fOnRVeT506FWvXrsXZs2flgVUsFsufqUxERERE5ZMg7hKQk5ODnTt3IjMzU+G+i6dPn4aTkxPq16+PMWPG5PsscCIiIiL6dGn1PqxXrlyBp6cnXr9+DWNjY6xevVo+hnXnzp0wMjKCvb09EhISMHv2bOTm5uL48eMQi8X5rjM+Pr60yiciIiKiYlDYg0y0Glizs7Px8OFDpKWlITIyEmFhYdi7dy9cXV2V2iYmJqJ27dpYt24dfHx8Sqym+Pj4Ynt6zaeOfaUZ9pf62FeaYX+pj32lPvaVZthf6itKX2n1tlYGBgbyx9XVrVsX586dw4oVK7B8+XKltpaWlrCyssr3edVERERE9GkSxBjWPLm5ucjOzlY5Lzk5GYmJibwIi4iIiKic0doR1unTp8PT0xPW1tbIyMhAeHg4oqOjsX37dmRkZCAkJAQ+Pj6wsLBAQkICZs6cCTMzM3h5eWmrZCIiIiLSAq0F1qSkJAwdOhRPnz6Fqakp3NzcEB4ejjZt2iArKwtXr17F1q1bkZaWBgsLCzRv3hzr16/n4/GIiIiIyhmtBdaVK1fmO8/IyAgRERGlWA0RERERCZWgxrASEREREX2IgZWIiIiIBI2BlYiIiIgEjYGViIiIiASNgZWIiIiIBE2rT7oiIiL1JSYa4NEj5eMM6ekOSEsz1EJFxc/aOheWlqofIENE5RcDKxFRGfHokQ4mTtRVmi6VGkEsVp5eFi1YAFhaarsKIhIaDgkgIiIiIkFjYCUiIiIiQWNgJSIiIiJBY2AlIiIiIkHjRVdE9EnI7wr6T0lW1qdxYRURkaYYWInok5DfFfSfkqAgbVdARKQdn/bhCCIiIiIq8xhYiYiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQGFiJiIiISNAYWImIiIhI0BhYiYiIiEjQtBZYV69eDQ8PD9ja2sLW1hbt2rXDwYMH5fNlMhmCg4Ph4uKCatWqoXPnzrh27Zq2yiUiIiIiLdFaYLWyssKMGTNw4sQJREVFoUWLFujTpw8uX74MAFiyZAlCQ0Mxd+5cHDt2DGZmZujatStevnyprZKJiIiISAu0Flg7d+6Mdu3aoUaNGnBycsLUqVNhYmKCs2fPQiaTYeXKlRg7dix8fX3h6uqKlStXIiMjA+Hh4doqmYiIiIi0QBBjWHNycrBz505kZmaiUaNGuH//PpKSktC6dWt5GyMjI3h4eCAuLk6LlRIRERFRadPT5savXLkCT09PvH79GsbGxti4cSPc3NzkodTMzEyhvZmZGRITEwtcZ3x8/EfXVRzrKC/YV5phf6lP075KT3eAVGpUQtUIw5s3upBKc1TOk0qlpVxNyUhPz0J8/L0S3Qa/h+pjX2mG/aW+D/vK2dm5wPZaDazOzs44deoU0tLSEBkZiYCAAOzdu1c+XyQSKbSXyWRK01St82PEx8d/9DrKC/aVZthf6itKX6WlGUIs1i2hioRBXx8Qi5V/bEulUojFYi1UVPxMTfVK9HvC76H62FeaYX+pryh9pdXAamBggBo1agAA6tati3PnzmHFihWYOHEiAODp06ewsbGRt3/+/LnSUVciIiIi+rQJYgxrntzcXGRnZ8Pe3h4WFhaIioqSz3v9+jVOnz6Nxo0ba7FCIiIiIiptWjvCOn36dHh6esLa2lp+9X90dDS2b98OkUiEgIAALFy4EM7OznBycsKCBQtgbGyMb775RlslExEREZEWaC2wJiUlYejQoXj69ClMTU3h5uaG8PBwtGnTBgDw/fffIysrC4GBgUhNTUX9+vURERGBihUraqtkIiIiItICrQXWlStXFjhfJBIhKCgIQUFBpVQREREREQmRoMawEhERERF9iIGViIiIiASNgZWIiIiIBI2BlYiIiIgEjYGViIiIiASNgZWIiIiIBI2BlYiIiIgEjYGViIiIiASNgZWIiIiIBI2BlYiIiIgEjYGViIiIiASNgZWIiIiIBI2BlYiIiIgEjYGViIiIiASNgZWIiIiIBI2BlYiIiIgEjYGViIiIiASNgZWIiIiIBI2BlYiIiIgEjYGViIiIiASNgZWIiIiIBI2BlYiIiIgEjYGViIiIiASNgZWIiIiIBI2BlYiIiIgEjYGViIiIiASNgZWIiIiIBI2BlYiIiIgEjYGViIiIiARN7cAaExOD58+f5zs/OTkZMTExxVIUEREREVEetQOrt7c3oqKi8p1/4sQJeHt7F0tRRERERER51A6sMpmswPnZ2dnQ0VF/hMGiRYvQqlUr2NrawtHREf7+/rh69apCm4CAAEgkEoV/bdu2VXsbRERERFT26RU0Mz09HWlpafLXL168wIMHD5TapaamYufOnbC0tFR7w9HR0Rg0aBDq1asHmUyGX375BV26dEFcXBwqV64sb9eyZUusWrVK/trAwEDtbRARERFR2VdgYF2xYgXmzZsHABCJRAgKCkJQUJDKtjKZDFOnTlV7wxEREQqvV61aBTs7O5w5cwYdO3aUTxeLxbCwsFB7vURERET0aSkwsLZs2RKGhoaQyWSYOXMm/Pz8ULt2bYU2IpEIFSpUQN26ddGgQYMiF5KRkYHc3FxIJBKF6adPn4aTkxMqVaqEpk2bYurUqTAzMyvydoiIiIiobBGlpqYWPDj1/wsJCYG3tzfc3NxKpJABAwbg9u3bOH78OHR1dQEAO3fuhJGREezt7ZGQkIDZs2cjNzcXx48fh1gsVrme+Pj4EqmPiITtwQMHTJlipO0yStSUKbqYMydH22WUqDlzsmBre0/bZRBRKXN2di5wfoFHWN83efLkjy4mPz/++CPOnDmDAwcOyMMqAHTr1k3+/25ubnB3d0ft2rVx8OBB+Pj4qFxXYTtcmPj4+I9eR3nBvtIM+0t9RemrtDRDiMW6hTcsw/T1AbFY+ce2VCrN94/4ssbUVK9Evyf8HqqPfaUZ9pf6itJX+QbWLVu2AAB69uwJkUgkf12YXr16aVRAUFAQIiIisGfPHjg4OBTY1tLSElZWVrhz545G2yAiIiKisivfwDpixAiIRCJ069YNBgYGGDFiRKErE4lEGgXWSZMmISIiAnv37kXNmjULbZ+cnIzExERehEVERERUjuQbWP/77z8A/3cbqbzXxWXixInYtm0bNm7cCIlEgqSkJACAsbExTExMkJGRgZCQEPj4+MDCwgIJCQmYOXMmzMzM4OXlVay1EBEREZFw5RtY7ezsCnz9sdasWQMA8PX1VZg+adIkBAUFQVdXF1evXsXWrVuRlpYGCwsLNG/eHOvXr0fFihWLtRYiIiIiEi61L7oqbqmpqQXONzIyUrpXKxERERGVPxoF1uPHjyMsLAz37t1DSkqK0uNaRSIRLly4UJz1EREREVE5p3ZgXblyJaZMmYLPPvsMDRo0QK1atUqyLiIiIiIiABoE1tDQUDRt2hQ7d+6UX4hFRERERFTSdNRtmJycDD8/P4ZVIiIiIipVagdWd3d3JCQklGQtRERERERK1A6sc+bMwebNm3Hy5MmSrIeIiIiISIHaY1iDg4NhamqKLl26wNHREba2ttDVVXxut0gkwvbt24u9SCIiIiIqv9QOrNevX4dIJIKNjQ2kUilu3bql1EYkEhVrcUREREREagfWS5culWQdREREREQqqT2GlYiIiIhIG9Q+wvrgwQO12tna2ha5GCIiIiKiD6kdWL/88ku1xqi+ePHiowoiIiIiInqf2oF1+fLlSoE1JycH9+/fx9atW2Fubo7BgwcXe4FEREREVL6pHVj79OmT77yxY8eidevWyMjIKJaiiIiIiIjyqB1YC2JiYoI+ffpgxYoVGDp0aHGskoiIyiGRSBf//GNYYutPT3dAWlrJrV8d1ta5sLTM1moNRGVNsQRWANDX10diYmJxrY6IiMqh58+B4GDdwhsWkVRqBLG45NavjgULAEtLrZZAVOYUy22tLl26hN9++w2ff/55cayOiIiIiEjuo+8SkJaWhvT0dJiYmCA0NLRYiyMiIiIiUjuwNm3aVCmwikQiSCQS1KhRA926dYNEIinu+oiIiIionFM7sK5cubIk6yAiIiIiUomPZiUiIiIiQWNgJSIiIiJBY2AlIiIiIkFjYCUiIiIiQWNgJSIiIiJBUyuwvn79GnPnzsWxY8dKuh4iIiIiIgVq3dbK0NAQixcvxrx580q6HiIqAYmJBnj0qOycUCnK896zsrT7uE0iIio5at+HtXbt2rhz505J1kJEJeTRIx1MnFh2Al1RnvceFFRCxRARkdapfcjl559/xoYNG3Dw4MGSrIeIiIiISIHaR1iXLl0KiUSCXr16wcrKCg4ODjAyMlJoIxKJsH379mIvkoiIiIjKL7WPsF6/fh1v376FjY0NdHR0kJCQgBs3bij9U9eiRYvQqlUr2NrawtHREf7+/rh69apCG5lMhuDgYLi4uKBatWro3Lkzrl27pv7eEREREVGZp/YR1kuXLhXrhqOjozFo0CDUq1cPMpkMv/zyC7p06YK4uDhUrlwZALBkyRKEhoYiNDQUzs7OmDdvHrp27YqzZ8+iYsWKxVoPEREREQmT2oG1uEVERCi8XrVqFezs7HDmzBl07NgRMpkMK1euxNixY+Hr6wsAWLlyJZydnREeHo6BAwdqo2wiIiIiKmUa3ecmJycH27dvx6hRo+Dv74/Lly8DAFJTU7Fr1y48efKkyIVkZGQgNzcXEokEAHD//n0kJSWhdevW8jZGRkbw8PBAXFxckbdDRERERGWL2kdY09LS4Ofnh3PnzsHExASZmZkYMWIEAKBixYqYMmUKevbsiZ9//rlIhUyePBm1a9dGo0aNAABJSUkAADMzM4V2ZmZmSExMzHc98fHxRdp+ca+jvGBfaUZb/ZWe7gCp1KjwhgIilUo1av/mjS6k0pwSqkYYCtpHTftLqErjfdR2X6WnZyE+/p5Wa1AXf8Zrhv2lvg/7ytnZucD2agfWGTNm4Pr169ixYwfq1q0LJycn+TxdXV14e3vj8OHDRQqsP/74I86cOYMDBw5AV1fx3osikUjhtUwmU5r2vsJ2uDDx8fEfvY7ygn2lGW32V1qaocb3NdUmqVQKsVis0TL6+oBYrLVRTqUiv30sSn8JVUm/j0LoK1NTvTLxs5M/4zXD/lJfUfpK7SEB+/btw9ChQ9G2bVuVgdHR0REPHjzQaOMAEBQUhJ07dyIyMhIODg7y6RYWFgCAp0+fKrR//vy50lFXIiIiIvp0qR1YU1NTUb169Xzny2QyZGdna7TxSZMmITw8HJGRkahZs6bCPHt7e1hYWCAqKko+7fXr1zh9+jQaN26s0XaIiIiIqOxS+7yLnZ2d0n1S3xcTE6MwTKAwEydOxLZt27Bx40ZIJBL5mFVjY2OYmJhAJBIhICAACxcuhLOzM5ycnLBgwQIYGxvjm2++UXs7RERERFS2qR1Yu3fvjl9//RXe3t6oVasWgP8bX7pq1Srs3bsXv/zyi9obXrNmDQDIb1mVZ9KkSQj6/w8F//7775GVlYXAwECkpqaifv36iIiI4D1YiYiIiMoRtQPruHHj8M8//8DHxwdOTk4QiUSYPHkyXrx4gaSkJHTu3BnDhg1Te8OpqamFthGJRAgKCpIHWCIiIiIqf9QOrPr6+ti+fTt27NiBv/76CyKRCG/fvkWdOnXg5+eHHj16FHj1PhERERFRUWh875Du3buje/fuJVELEREREZGSIt3s7vLly/JbWNna2sLNzY1HV4mIiIioRGgUWHfu3Ilp06bh8ePHkMlkAN6NM7WyssK0adN45JWIiIiIip3agXXTpk0YNWoUnJ2dMWPGDDg5OUEmk+H27dvYsGEDhg0bhuzsbPTp06ck6yUiIiKickbtwLpo0SLUr18fe/fuhaGhocK8IUOGoFOnTli0aBEDKxEREREVK7WfdPXo0SN0795dKawCgKGhIfz9/fH48eNiLY6IiIiISO3A6uLigsTExHznP378GJ9//nmxFEVERERElEftwDpz5kyEhYVh165dSvN27tyJDRs2YNasWcVaHBERERGR2mNYly1bhqpVq2LQoEGYPHkyqlevDpFIhDt37uDZs2dwdHTE0qVLsXTpUvkyIpEI27dvL5HCiYiIiKh8UDuwXr9+HSKRCDY2NgAgH68qFothY2MDqVSKGzduKCzDe7MSERER0cdSO7BeunSpJOsgIiIiIlJJ7TGsRERERETawMBKRERERILGwEpEREREgsbASkRERESCxsBKRERERILGwEpEREREgqZ2YK1Tpw7279+f7/wDBw6gTp06xVIUEREREVEetQNrQkICMjMz852fmZmJBw8eFEtRRERERER5NBoSUNCTq27duoWKFSt+dEFERERERO8r8ElXmzdvxpYtW+SvFyxYgLCwMKV2qampuHr1Ktq3b1/8FRIRERFRuVZgYM3MzERSUpL8dVpaGnJzcxXaiEQiVKhQAf3798fkyZNLpkoiIiIiKrcKDKxDhgzBkCFDAABffvklQkJC0KlTp1IpjIiIiIgIKCSwvu/ixYslWQcRERERkUpqB9Y8L1++xMOHD5GSkgKZTKY0v2nTpsVSGBERERERoEFgTUlJwaRJk7Br1y7k5OQozZfJZBCJRHjx4kWxFkhERERE5ZvagXXcuHHYu3cvhgwZgqZNm0IikZRgWURERERE76gdWI8cOYJhw4Zhzpw5JVkPEREREZECtR8cYGBgAEdHx5KshYiIiIhIidqB1dfXF4cPHy7WjcfExKBnz56oVasWJBIJNm3apDA/ICAAEolE4V/btm2LtQYiIiIiEja1A+vo0aPx5MkTDB8+HGfPnsWTJ0/w7NkzpX+ayMzMhKurK0JCQmBkZKSyTcuWLXHjxg35vx07dmi0DSIiIiIq29Qew1q/fn2IRCJcuHAB27dvz7edJncJ8PT0hKenJwBgxIgRKtuIxWJYWFiovU4iIiIi+rSoHVh/+OEHiESikqxFpdOnT8PJyQmVKlVC06ZNMXXqVJiZmZV6HURERESkHWoH1qCgoJKsQ6W2bdvC29sb9vb2SEhIwOzZs+Hj44Pjx49DLBarXCY+Pv6jt1sc6ygv2Fea0VZ/pac7QCpVPexGqKRSqUbt37zRhVSqfI/oT0lB+6hpfwlVabyP2u6r9PQsxMff02oN6uLPeM2wv9T3YV85OzsX2F7jJ10BQE5ODtLS0mBqago9vSKtQi3dunWT/7+bmxvc3d1Ru3ZtHDx4ED4+PiqXKWyHCxMfH//R6ygv2Fea0WZ/paUZQizW1cq2i0Iqleb7R2l+9PUBsbjkfh4JQX77WJT+EqqSfh+F0Fempnpl4mcnf8Zrhv2lvqL0ldoXXQHAuXPn0KVLF1hZWcHJyQkxMTEAgOTkZPTo0QMnTpzQaOOasrS0hJWVFe7cuVOi2yEiIiIi4VA7sP7vf/9Dp06dcPfuXfTs2RMymUw+r2rVqsjIyMCff/5ZIkXmSU5ORmJiIi/CIiIiIipH1D7vMmvWLDg6OuLo0aPIzMzEhg0bFOY3b94c27Zt02jjGRkZ8qOlubm5ePjwIS5evIjKlSujcuXKCAkJgY+PDywsLJCQkICZM2fCzMwMXl5eGm2HiIiIiMoutY+wnjt3Dn379oWhoaHKuwVYW1sjKSlJo42fP38eLVq0QIsWLZCVlYXg4GC0aNECv/zyC3R1dXH16lX07t0bDRo0QEBAAJycnHDo0CFUrFhRo+0QERERUdml9hFWHR0d6Ojkn2+TkpLyvfl/fpo3b47U1NR850dERGi0PiIiIiL69Kh9hNXd3R0HDhxQOS87Oxs7duxAo0aNiq0wIiIiIiJAg8A6fvx4nDx5EqNGjcKlS5cAAE+ePMGRI0fg4+ODu3fvYsKECSVWKBERERGVT2oPCWjVqhVWrVqFwMBAbN68GQAQEBAAmUyGSpUqYc2aNWjYsGGJFUpERERE5ZNGd2f+5ptv0KlTJ0RFReH27dvIzc1F9erV0aZNG5iYmJRUjURERERUjmn8OJEKFSqgc+fOJVELEREREZEStcew7t+/H4GBgfnODwwMzPeiLCIiIiKiolI7sC5btgyvXr3Kd/7r16+xZMmSYimKiIiIiCiP2oH16tWrcHd3z3d+nTp1cP369eKoiYiIiIhITu3A+vbtW2RlZeU7PysrC1KptFiKIiIiIiLKo3ZgdXV1RWRkJHJzc5Xm5ebmIjIyEi4uLsVaHBERERGR2oF1+PDh+Pfff9GrVy9cuHABUqkUUqkUFy5cQO/evfHvv/9i2LBhJVkrEREREZVDat/Wqlu3brh79y6Cg4Nx+PBhAIBIJIJMJoNIJMKkSZPg7+9fYoUSERERUfmk0X1YJ06ciG+++QZ79uzBvXv3IJPJUL16dXh7e8PBwaGESiQiIiKi8kytwJqVlYUePXrA398fffv2xejRo0u6LiIiIiIiAGqOYTUyMsJ///2HnJyckq6HiIiIiEiB2hddNWvWDLGxsSVZCxERERGRErUD69y5c3Hu3DlMnToV9+7dU3l7KyIiIiKi4qb2RVcNGzaETCZDaGgoQkNDoaOjA319fYU2IpEIjx8/LvYiiYiIiKj8Ujuwdu3aFSKRqCRrISIiIiJSonZgXblyZUnWQURERESkktpjWImIiIiItEGjwJqQkIAxY8bA3d0dtra2iI6OBgAkJydjwoQJuHDhQknUSERERETlmNpDAm7cuIEOHTogNzcXDRo0QEJCgvy+rFWrVsXZs2chlUqxfPnyEiuWiIiIiMoftQPrtGnTULFiRRw5cgS6urpwcnJSmO/p6Ym//vqruOsjIiIionJO7SEBsbGxGDx4MMzNzVXeLcDW1haJiYnFWhwRERERkdqB9e3btzA2Ns53fkpKCnR1dYulKCIiIiKiPGoHVldXV5w6dUrlPJlMhj179sDd3b246iIiIiIiAqBBYA0ICMDu3bsxb948vHjxAgCQm5uLmzdv4rvvvsP58+cxevToEiuUiIiIiMontS+66tatGx48eIA5c+YgJCREPg0AdHV1MXv2bLRr165kqiQiIiKickvtwAoAY8eOxTfffIPIyEjcuXMHubm5qF69Onx8fGBvb19SNRIRERFROVZoYJVKpdi/fz/u3buHKlWqoH379hgxYkRp1EZEREREVPAY1qSkJHh4eGDQoEGYOXMmxo4diwYNGiAmJqZYNh4TE4OePXuiVq1akEgk2LRpk8J8mUyG4OBguLi4oFq1aujcuTOuXbtWLNsmIiIiorKhwMA6e/Zs3Lt3DyNGjMC2bdsQHBwMsViMH374oVg2npmZCVdXV4SEhMDIyEhp/pIlSxAaGoq5c+fi2LFjMDMzQ9euXfHy5cti2T4RERERCV+BQwKOHTuGXr16Yfbs2fJp5ubmGDx4MB49egRra+uP2rinpyc8PT0BQGmYgUwmw8qVKzF27Fj4+voCAFauXAlnZ2eEh4dj4MCBH7VtIiIiIiobCh0S0LhxY4VpTZo0gUwmw8OHD0u0sPv37yMpKQmtW7eWTzMyMoKHhwfi4uJKdNtEREREJBwFHmHNycmBoaGhwrS8169fvy65qvAuLAOAmZmZwnQzM7MCHwEbHx//0dsujnWUF+wrzWirv9LTHSCVKg+7ETKpVKpR+zdvdCGV5pRQNcJQ0D5q2l9CVRrvo7b7Kj09C/Hx97Rag7r4M14z7C/1fdhXzs7OBbYv9C4B9+7dw7///it/nZ6eLt+QiYmJUvv69eurVai6RCKRwmuZTKY07X2F7XBh4uPjP3od5QX7SjPa7K+0NEOIxWXn0clSqRRisVijZfT1AbFYozv1lTn57WNR+kuoSvp9FEJfVaokRlpaba3WoI709JcwNa2o8XLW1rmwtMwugYqEjb8T1VeUvir0p0JwcDCCg4OVpn944VVekMx7CtbHsrCwAAA8ffoUNjY28unPnz9XOupKRERUVjx/DgQHC/8PSKnUqEh/6C5YAFhalkBBVK4VGFhDQ0NLqw4l9vb2sLCwQFRUFOrVqwfg3TCE06dPY+bMmVqri4iIiIhKV4GBtXfv3iW68YyMDNy5cwcAkJubi4cPH+LixYuoXLkybG1tERAQgIULF8LZ2RlOTk5YsGABjI2N8c0335RoXUREREQkHFod8HX+/Hl4e3vLX+cNP+jVqxdWrlyJ77//HllZWQgMDERqairq16+PiIgIVKyo+ZgaIiIiIiqbtBpYmzdvjtTU1Hzni0QiBAUFISgoqPSKIiIiIiJBKfA+rERERERE2sbASkRERESCxsBKRERERILGwEpEREREgsbASkRERESCxsBKRERERILGwEpEREREgsbASkRERESCxsBKRERERILGwEpEREREgsbASkRERESCpqftAoi0LTHRAI8elfzfbunpDkhLMyzx7aiSlaWrle0SEREVBwZWKvcePdLBxIklH+ikUiOIxdoJjkFBWtksERFRseCQACIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0QQfW4OBgSCQShX81a9bUdllEREREVIr0tF1AYZydnbF37175a11dXS1WQ0RERESlTfCBVU9PDxYWFtoug4iIiIi0RNBDAgDg3r17qFWrFr788kt89913uHfvnrZLIiIiIqJSJOgjrA0aNMCKFSvg7OyM58+fY/78+fD09MSZM2dQpUoVlcvEx8d/9HaLYx3lxafQV+npDpBKjUplW1KptFS286E3b3QhleZoZdtFpWlflcV91FRB+6itz1ZxK433Udt9VZY+q0Xpq/T0LMTH3yv+YsqAT+F3Ymn5sK+cnZ0LbC/owNquXTuF1w0aNIC7uzs2b96MUaNGqVymsB0uTHx8/Eevo7z4VPoqLc0QYnHJj42WSqUQi8Ulvh1V9PUBsVjQX3cFRemrsraPRZHfPmrzs1XcSvp9FEJflZXPalH7ytRU75P43aCpT+V3YmkoSl8JfkjA+0xMTODi4oI7d+5ouxQiIiIiKiVlKrC+fv0a8fHxvAiLiIiIqBwR9DmJn376CR06dICNjY18DOurV6/Qq1cvbZdGRERERKVE0IH18ePHGDx4MJKTk/HZZ5+hQYMGOHz4MOzs7LRdGhERERGVEkEH1nXr1mm7BCIiIiLSsjI1hpWIiIiIyh8GViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISNAZWIiIiIhI0BlYiIiIiEjQGViIiIiISND1tF0BERESfDpFIF//8Y6jtMkqUtXUuLC2ztV1GucLASkRERMXm+XMgOFhX22WUqAULAEtLbVdRvnBIABEREREJGgMrEREREQkahwRQgRITDfDokeq/a9LTHZCWVvbHKWVlfdqnroiIiMo6BlYq0KNHOpg4UXWgk0qNIBaX/bAXFKTtCoiIiKggHBJARERERILGwEpEREREgsbASkRERESCxsBKRERERILGwEpEREREgsbASkRERESCxsBKRERERILGwEpEREREgsbASkRERESCxiddEREREWlAJNLFP/8oPpr8U3lceR5r61xYWmZruwy5MhFY16xZg6VLlyIpKQkuLi4IDg6Gh4eHtssiIiKicuj5cyA4WPHR5J/K48rzLFgAWFpqu4r/I/ghAREREZg8eTImTJiAkydPolGjRujevTsePHig7dKIiIiIqBQIPrCGhoaid+/e6N+/Pz7//HPMnz8fFhYWWLdunbZLIyIiIqJSIEpNTZVpu4j8ZGdnw9LSEmvXrkWXLl3k0ydOnIirV69i//792iuOiIiIiEqFoI+wJicnIycnB2ZmZgrTzczM8PTpUy1VRURERESlSdCBNY9IJFJ4LZPJlKYRERER0adJ0IG1atWq0NXVVTqa+vz5c6WjrkRERET0aRJ0YDUwMIC7uzuioqIUpkdFRaFx48ZaqoqIiIiISpPg78M6cuRIDBs2DPXr10fjxo2xbt06PHnyBAMHDtR2aURERERUCgR9hBUA/Pz8EBwcjPnz56N58+Y4c+YMtm/fDjs7u1LZvkwmQ7du3SCRSLB79+5S2WZZM2bMGLi7u6NatWpwdHREr169cOPGDW2XJUgpKSkIDAxEw4YNUa1aNbi5uWH8+PF48eKFtksTpD/++ANeXl6ws7ODRCLB/fv3tV2SoKxZswZffvklLCws8PXXXyM2NlbbJQlSTEwMevbsiVq1akEikWDTpk3aLkmwFi1ahFatWsHW1haOjo7w9/fH1atXtV2WYK1evRoeHh6wtbWFra0t2rVrh4MHD2q7rDJh4cKFkEgkCAwMVKu94AMrAAwePBiXLl3C06dPceLECTRt2rTUtr18+XLo6n46T64oCXXr1sWKFSsQFxeHnTt3QiaToUuXLnjz5o22SxOcxMREJCYmYsaMGYiNjcWqVasQGxuLQYMGabs0QXr16hVat26NyZMna7sUweFDVdSXmZkJV1dXhISEwMjISNvlCFp0dDQGDRqEgwcPIjIyEnp6eujSpQtSUlK0XZogWVlZYcaMGThx4gSioqLQokUL9OnTB5cvX9Z2aYJ29uxZhIWFwc3NTe1lBH0fVm07f/48+vbti+PHj8PZ2RlhYWHw9fXVdlmCd/nyZTRr1gxnz56Fs7OztssRvEOHDsHf3x/379+HqamptssRpPPnz6NVq1b477//YG9vr+1yBKFNmzZwc3PD0qVL5dPq1asHX19fTJs2TYuVCZu1tTXmzZuHPn36aLuUMiEjIwN2dnbYtGkTOnbsqO1yygQHBwdMmzaNQxfzkZaWhq+//hpLlizBvHnz4Orqivnz5xe6XJk4wqoNL1++xKBBg7B48WLekUADmZmZ2LRpE2xsbEpt2EZZ9/LlS4jFYlSoUEHbpVAZkZ2djQsXLqB169YK01u3bo24uDgtVUWfooyMDOTm5kIikWi7FMHLycnBzp07kZmZiUaNGmm7HMEaO3YsfH198fXXX2u0nOAvutKW8ePHo02bNvD09NR2KWXCmjVrMG3aNGRmZsLZ2RmRkZEQi8XaLkvwUlNTMWfOHPTr1w96evw6knr4UBUqLZMnT0bt2rUZwApw5coVeHp64vXr1zA2NsbGjRs1OtVdnoSFheHOnTtYtWqVxsuWqyOss2fPhkQiKfDfqVOnsHXrVly+fBmzZs3Sdslao25f5enevTtOnjyJffv2wdHREf3798erV6+0uAelS9P+At4dje7VqxcsLS0xc+ZMLVVe+orSV6QaH6pCJenHH3/EmTNn8Oeff/JajgI4Ozvj1KlTOHLkCAYNGoSAgABeqKZCfHw8Zs6cidWrV8PAwEDj5cvVGNbk5GQkJycX2MbGxgYTJkzA1q1boaPzf3k+JycHOjo6aNSoEQ4cOFDSpWqdun2l6jR2dnY2HBwcsGjRIvTs2bOkShQUTfsrIyMD3bt3BwDs2LEDJiYmJV6jUBTls8UxrIqys7NhaWmJtWvXokuXLvLpEydOxNWrV7F//37tFSdwHMOqnqCgIERERGDPnj2oWbOmtsspU3x9fWFra4vly5druxRB2bRpE0aOHKnwx09OTg5EIhF0dHTw+PHjAs/MlqtzkFWrVkXVqlULbTd16lSMHj1aYZqHhwdmzZqFzp07l1R5gqJuX6kik8kgk8mQnZ1dzFUJlyb99fLlS3Tv3h0ymQzh4eHlKqwCH/fZonfef6jK+4E1KioKPj4+2iuMPgmTJk1CREQE9u7dy7BaBLm5ueXq95+6OnfujLp16ypMGzlyJBwdHTF+/PhCj7qWq8CqLisrK1hZWSlNt7GxgYODQ+kXJGB37txBZGQkWrZsiapVq+Lx48dYvHgxDAwM0L59e22XJzgvX76En58fXr58iU2bNuHVq1fyoROVK1cu0mmST1lSUhKSkpJw69YtAMCNGzeQlpYGW1tbVK5cWcvVaRcfqqK+jIwM3LlzB8C7MPHw4UNcvHgRlStXhq2trZarE5aJEydi27Zt2LhxIyQSCZKSkgAAxsbG5e6Pa3VMnz4dnp6esLa2RkZGBsLDwxEdHY3t27druzTByRvy9b4KFSqgcuXKcHV1LXR5Blb6KAYGBoiOjsby5cuRlpYGc3NzeHh44PDhw7CwsNB2eYJz4cIFnD17FgBQv359hXl79uxB8+bNtVGWYK1btw5z586Vv+7RowcAIDQ0tNyf0vXz88OLFy8wf/58JCUloVatWqX6UJWy5Pz58/D29pa/Dg4ORnBwMHr16oWVK1dqsTLhWbNmDQAo3cJx0qRJCAoK0kZJgpaUlIShQ4fi6dOnMDU1hZubG8LDw9GmTRttl/bJKVdjWImIiIio7ClXdwkgIiIiorKHgZWIiIiIBI2BlYiIiIgEjYGViIiIiASNgZWIiIiIBI2BlYiIiIgEjYGVqIwKDg5WugkzAKxYsQJ169ZF1apVUbt27dIvjIiIqJgxsBIJwKZNm+RPAZFIJLCwsICLiwv8/Pzw22+/4eXLl2qtJzo6Gj/++CPq1KmDZcuWITg4uIQrF77atWsr9K25uTnc3d0RGBgof4qPpjIyMhAcHIxTp04Vc7XCdOzYMfTt2xcuLi4wMzODnZ0d2rRpgzlz5uDx48faLq/ITp06Jf9cbNy4UWWbAQMGyL+TRKQ9fNIVkYBMnjwZ1atXx5s3b/D06VNER0cjKCgIoaGh2LJlC7744gt528DAQIwbN05h+bwA9euvv6o8+lpeubm5YcyYMQAAqVSKS5cuISwsDCdPnkRsbCx0dXU1Wl9mZqb8CVyf8tPJZDIZJkyYgHXr1sHFxQX9+vWDra0tsrKycOHCBaxatQrr16+XPzq3rDI0NMSOHTvQt29fhenp6ek4cOAADA0NIZPxGTtE2sTASiQgbdq0QcOGDeWvx48fjxMnTqBnz57o1asX/ve//8HIyAgAoKenBz09xa/w8+fPAaBYw+qrV69QoUKFYlufNlSrVg3+/v4K0ypVqoQFCxbg0qVLcHd3105hAhcaGop169Zh6NChCAkJgY6O4km5kJAQLF26tND1ZGVlyT+3QuTp6Ym9e/ciMTERlpaW8umRkZHIzc1F69atERUVpcUKiYhDAogE7uuvv0ZgYCAePHiA7du3y6d/OIZVIpFg7dq18v+XSCQKQwKioqLg5eUFGxsbWFlZwcvLC3FxcQrbylvn9evXMXz4cFSvXh1NmjQp0jpu376NcePGoXr16rC2tkb//v3x4sULpf2LioqCt7c3bG1tYWNjg6+//hobNmxQaHP+/Hn4+/vDzs4O1apVQ+vWrXHgwAHNO/M9ead49fX1FaY/efIE33//PVxcXGBubo569ephyZIl8iNs9+/fx+effw4AmDt3rryvAwICcPnyZUgkEuzatUu+vrt370IikSgcHQeAcePGoWbNmkXaz/T0dPz000+oXbs2zM3N8cUXX2D69OmQSqUK7SQSCcaNG4fDhw+jefPmsLCwQL169RAeHl5o/2RlZWHhwoX4/PPPERwcrBRWAcDU1BQ//fSTwrTatWujW7duOHnyJNq2bQsLCwv8+uuvAICUlBSMHz8en3/+OczNzdGoUSMsX75c4ejl/fv3IZFIsGnTJqXt1a5dGwEBAfLXeUNpTp48icDAQNSoUQPW1tbo168fnjx5Uug+5mnfvj1MTEyU+mXHjh1o27YtKleurHI5db4PCQkJmDBhAho2bAhLS0vY2dnB398f165dU2iXNzwhPDwcy5cvR+3atWFhYYF27drhv//+U2j79OlTjB49Gm5ubjA3N4eLiwv8/f1x5coVtfeZqKxhYCUqA/KODh47dizfNqtWrUKLFi3k/79q1Sp4e3sDAMLDw9GtWzfo6upiypQpmDJlCl68eAEfHx/8888/SusaOHAgUlJSMGXKFAwfPrxI6xg0aBAeP36MKVOmoF+/fti7dy9++OEHhTZbt26Fn58fnjx5gtGjR2PGjBmoX78+Dh48KG8THR2NDh064OnTpwgMDMSMGTNgYGCAXr16ITIyUq3+e/PmDZKTk5GcnIzHjx/j0KFDWLJkCdzd3eHq6ipv9+zZM7Rt2xYHDx5E//79MXfuXDRo0ADTpk1DUFAQAOCzzz7D/PnzAQBeXl7yvh44cCDc3NwgkUgQExMjX2dMTAx0dHTw8OFD3L9/Xz49NjYWX331lcb7mZWVBS8vL/z555/w8/PDvHnz0L59eyxfvhwDBw5U2vezZ89i5MiR6NSpE2bNmoUKFSpg6NChuHHjRoF9FhcXh5SUFHzzzTcaD5m4c+cO+vXrBw8PD8ydOxcNGzaEVCqFt7c3wsLC4OPjgzlz5sDe3h4//fQTfvzxR43W/6HJkyfjwoUL+OGHHzBgwAD8/fff8PPzQ3Z2tlrLGxoawtvbGzt27JBPe/LkCU6dOoUePXqoXEbd78P58+cRExMDb29vBAcHIyAgAOfPn0enTp1UjqFevnw5tmzZgqFDh2Ly5Mm4desW+vTpgzdv3sjb9O/fH7t370avXr2wYMECDBs2DLm5uWV+aAZRQTgkgKgMsLa2hqmpKe7evZtvG39/f5w5cwYnT55UOP2dmZmJiRMnwt/fHytXrpRPHzhwIJo0aYKZM2cqBT8nJyf8+eefH7WOmjVr4vfff5e/lslkWL16NRYuXIhKlSohPT0dP/zwA9zc3HDw4EEYGxsrtM3777hx49CoUSPs3r1bfpRvyJAhaN++PX7++Wf4+PgU2n8nT56Eo6OjwrT69etj69atEIlE8mmzZ8+GVCpFTEwMzM3N5ftYrVo1LF++HAEBAbC3t4ePjw8CAwPh5uamNNSgSZMmiI2Nlb8+ffo0Wrdujbi4OJw+fRr29vZITk7GzZs38d1332m8nytWrEB8fDyOHz8uP9ILALVq1cLEiRMRGxsLDw8P+fTr168jJiZG3rZLly744osvsHHjRsyaNSvfPrt+/bp8ve+TyWRKR8pNTU0VjlTfvXsXmzdvRqdOneTTfv/9d1y+fBlLly5Fv379AACDBw/Gt99+i99++w2DBw9Weo80sXfvXojFYgCAi4sLRo8ejc2bN2PAgAFqLd+jRw9s2rQJ169fh4uLC3bs2AETExN06NBB4Q8oQLPvQ7t27eDr66uwvL+/P7766iv8+eefmDhxosK89PR0xMbGwtDQEADg7OyMvn374tixY2jfvj3S0tJw+vRpzJo1C6NHj5Yv9+F4dqJPDY+wEpURJiYmyMjI0Hi5qKgopKamokePHvKjjMnJycjKykLLli1x+vRphaM3wLujo8W9jqZNmyInJwcPHz6UrzM9PR0TJkxQCKsA5CHy0qVLiI+PR48ePZCSkiLfbkpKCtq2bYt79+4hISGh0D6oW7cu/vrrL/z111/Ytm0bfv75Z9y+fRt9+/ZFVlYWgHdBbPfu3Wjfvj10dXUV9rNNmzbIzc1VOHKaHw8PD1y7dg0pKSkA3h1Jbd68ORo1aiQPsjExMZDJZPIjrJrs565du9C4cWN89tlnCjW2bNkSwLtw/r7mzZsrBFtzc3M4Ozvj3r17Be5H3p0pKlasqDD9xYsXcHR0VPh35swZhTbW1tYKYRUADh48iKpVq6JPnz7yaSKRCGPGjIFMJsOhQ4cKrKcgAwcOlIdVAOjVqxcqVaqk0TqbN28OS0tL+VHWHTt2wMvLSx4c36fJ9+H98d+vXr3CixcvUKlSJTg6OuLChQtK6+7Tp4/CNps1awYA8vfL0NAQ+vr6iI6Oln/GiMoDHmElKiMyMjLw2Wefabzc7du3AQBdu3bNt01aWprCuh0cHD56Hba2tgrz88bb5v2SzTta/P4p+fxqHz16tMLRpPc9f/4cdnZ2+a4DAKpUqSIPdMC7MYvOzs749ttvsWHDBgwbNgzPnz9HamoqNm7cmO8tjvIuaiuIh4cHZDIZYmNjUb9+fdy9exceHh54+/Yttm7dCuBdiDU1NZWPa9VkP2/fvo3Lly/nezTywxo/fB+Ad+9FYWEnL6h+eEs1U1NT/PXXXwDeDWNYsGCB0rL29vZK0xISEuDo6Kg0vCAvTKvzh0d+PuwLPT092Nvb48GDB2qvQ0dHB35+ftixYwe6d++OixcvYubMmSrbavJ9eP36NX755Rds375daVxt1apVlZYr7HsjFosxbdo0TJs2Dc7OzmjQoAHatWuHHj16qHyviT4VDKxEZcCjR4+Qnp6OGjVqaLxsbm4ugHenkq2srFS2MTU1VXj94RXdRVlHfuMe3z/dD0DhlHx+tU+fPj3fK/mdnJzyXb4geeN9Y2Nj5WMAAeCbb75Rur1RHnX6393dHcbGxoiNjYVUKkWFChXg7u6ON2/eYNasWXj27Jl8/GreqX9N9jM3NxctWrTA+PHjVbb78P0p7H3Ij4uLCwDg6tWr8PLykk/X19eXh//k5GSVy37MHQHU+Tyos0xRbkPVvXt3hIaGIjAwENWqVcv3lmWafB8mT56MDRs2YOjQoWjSpAlMTU2ho6ODoKAglfujzvs1atQoeHl5Yf/+/Th+/Djmz5+PRYsWYfPmzfj666812meisoKBlagM2LZtGwCgdevWGi9bvXp1AO8uFnr/KGNpr+NDeeHv6tWrSlfLf7hdExOTYttunrdv3wJ4Nx4ReLdvpqamePv2baHbKihU6enpoUGDBoiNjUV2djYaNmwIfX191K9fH2KxGAcOHMCVK1fg5+cnX0aT/axevToyMjKKvT8+1LhxY1SuXBnh4eGYMGGCxhdefcjOzg7//fcfcnJyFNZ18+ZN+XwA8ivy09LSFJaXSqX5Xvl/69YttGrVSv767du3SEhIQNOmTTWq0d3dHTVr1sSpU6cwYsSIfPdZk+9DREQEevbsiZCQEIXpqampqFKlikb1vc/BwQEjRozAiBEj8PDhQ7Ro0QKLFy9mYKVPFsewEgnciRMnMH/+fNjb2+d7xXJB2rRpI7/n6Ie3PQLUO81dHOv4UKtWrWBqaopFixbh1atXCvPyjia5u7vD0dERy5YtUwowRd1unrzxjXmn5XV1deHj44O9e/eqHFuYlpamNC4xNTVV5bo9PDxw8eJFHDlyRH4BlFgsRr169bB06VLk5OQoXBilyX76+fnh3Llz2L9/v1K7rKysIo1zVsXIyAgTJkzAzZs38z0aqMlRzPbt2+P58+fYsmWLwvLLli2DSCSCp6cngHdDET777DOlp4itW7cOOTk5Kte9fv16hc/lli1bkJaWhnbt2qldX545c+Zg0qRJSmOw36fJ90FXV1epn8LDw5GYmKhxbcC7cbB5467z2NjYwMzMLN/PI9GngEdYiQTk6NGjuHPnDt6+fYtnz57h5MmTiIqKgq2tLbZs2aLyApDCVKxYEUuWLMGgQYPQrFkzdO/eHRYWFnj06BFOnToFY2PjQu/LWRzr+JCpqSmCg4MxatQotGrVCt27d0eVKlVw7do1JCYmYuPGjdDR0cHy5cvRrVs3NGnSBH369IGdnR2ePHmCs2fP4sGDB0oX/Kjy5MkT+VHq7OxsXLlyBX/88QeqVq2KoUOHyttNnz4dMTEx6NChA7799lu4urri5cuXuHr1Kvbs2YNz587BwsICJiYmcHZ2RkREBJycnFClShXY29ujQYMGAICvvvoKOTk58vGreZo2bYoFCxbAyMhI4dS/Jvs5evRoHDp0CN9++y169OiB+vXrQyqV4tatW9i1axd27Nih8PCJjzFy5Ejcvn0bv//+O06ePAkfHx/Y2toiMzMTN27cwM6dO2FsbKxyLOaH+vXrhw0bNmDs2LG4dOkSnJyccPjwYRw6dAjDhw9XGIc6YMAALFiwACNGjEDDhg1x/vx5nDhxosDteHt7o1u3bkhISMDvv/8OFxcX9O7dW+N9bteuXaFBV5PvQ8eOHbF161ZUrFgRrq6uuHTpEiIiIpTGiavr1q1b8PHxQZcuXeDi4gKxWIxDhw7hxo0bBd71gaisY2AlEpC804YGBgaoXLkyXF1dERwcjD59+ihdra2JLl26wNLSEosWLcKKFSuQlZUFCwsLNGjQQH6LodJYx4f69OkDMzMzLF68GIsWLYKuri4cHR0xePBgeZuvvvoKR48exbx58/DHH38gPT0dZmZm+OKLL+T3Ri3MlStXMGzYMADvwmHVqlXh5eWFKVOmKIxB/Oyzz3D06FHMnz8f+/btwx9//IFKlSrByckJkydPVriBfGhoKIKCgvDTTz9BKpWiV69e8sDasGFDGBgYAIB8Wt6+5E3Lm6/pfhoZGSEyMhJLlixBRESEPDQ6ODggICAAzs7OavWJOkQiERYvXozOnTtj3bp1CAsLQ3JyMipUqAAnJycMHz4cAwYMyHcc5/sMDQ0RGRmJWbNmYdeuXUhJSYG9vT1mzZqFUaNGKbSdOHEiXrx4gYiICPz1119o1qwZdu/eLb+v8IdCQkIQGRmJuXPnQiqVon379pg/f77CnQOKm7rfh5CQEOjr62PXrl3YuHEj3N3dsXPnTkydOrVI27WxsUH37t1x8uRJhIeHQyQSyY/Of/vtt8W1e0SCI0pNTeUDkomIqMzZtGkTRo4cicOHDxfbUWUiEiaOYSUiIiIiQWNgJSIiIiJBY2AlIiIiIkHjGFYiIiIiEjQeYSUiIiIiQWNgJSIiIiJBY2AlIiIiIkFjYCUiIiIiQWNgJSIiIiJBY2AlIiIiIkH7fwx7Ez5pO1mkAAAAAElFTkSuQmCC\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": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAFuCAYAAACIv7K8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABPdklEQVR4nO3dd1wT9/8H8FdABFQgihAXoAKV4qKiYnHUjQNEceAedeKqC4XaautoQNHWQVHr+ErdAyuOulEQFK0VsU7cE1SmKESB/P7wQX7GEAgQIMTX8/Hw0ebuc3fvfHIhr9x97iJISUmRgoiIiEgL6ZR1AUREREQlhUGHiIiItBaDDhEREWktBh0iIiLSWgw6REREpLUYdIiIiEhrMeiQxvHy8oJQKERERITcdKFQiJ49e5bYdsVicZ7bpQ9y+2fr1q1lXYpSPXv2hFAoxMOHD8u6lFKR13uC+3H+tm7dCqFQCLFYXNalAODrVRoYdD5TQqFQ7l+1atVgZWWFbt26YdOmTcjOzi7rEtWuPHxQfyq35vz+KEdERJR4CNQUmh5kGjduDKFQiNq1a+P58+d5tvn222816oNN2/YxTQsyVPYqlHUBVLbmzJkDAMjOzsb9+/dx8OBBnD9/HqdPn8bmzZvLuDp5Fy5cgKGhYYmtf9y4cejbty/q1KlTYtugz8ObN2+waNEiBAYGlnUppOH4d6fkMeh85nx9feUeX7t2DZ07d8b+/fsRFRUFZ2fnMqpM0RdffFGi6zc1NYWpqWmJboM+D9bW1ti+fTvGjx+PJk2alHU5pMH4d6fk8dQVyWnYsCFat24NALh06RKA/z9s7eXlhZs3b2Lo0KGoX78+hEIhYmNjZcvu378f7u7uqFu3LszNzdGsWTP89NNPSEtLy3Nbp0+fRvfu3VGrVi3UrVsXgwcPxq1bt5TWpuzQeXZ2NoKDg9G9e3dYWVlBJBKhSZMmGDNmDC5fvgzgwykPf39/AMCkSZPkTtvlngbJ71x5eHg4+vfvj3r16sHc3BxNmzbFnDlz8PLlS4W2H48x2r9/Pzp27IiaNWuibt26GDVqFJ4+far0OapTUeuIiYmRfcO0sLCAu7s7oqOj891WQkICfHx80KxZM4hEIlhZWaFPnz44c+aMQtuPTy1ER0fDw8MDVlZWEAqFSElJUboNoVCIyMhIAEDTpk1lr1/jxo3zbL9p0yY4OztDJBLB1tYWU6dOVbr+wtSvip9++gk5OTn44YcfVF4mv1NDmjiO448//oBQKISfn1+e89PS0lCrVi00bNhQdir849PHR44cQZcuXWTv/5EjR+L+/ft5rishIQHe3t5o2rQpzM3NUa9ePQwYMABnz56Va+fl5YVJkyYBAPz9/eXe53n1XWxsLAYMGABLS0vUrFkT3bt3x/nz5/OsIScnB8HBwXBxcYGlpSVEIhG+/vprLF++HO/evVNoHxERAU9PTzRs2BDm5uawsbFB+/btMXfuXEil///LS8peW1WXp4LxiA6p7P79++jatSsaNGiAgQMHIjU1FZUqVQIAzJw5Exs2bEDt2rXh6uoKoVCIf/75B7/99huOHTuGo0ePwsjISLau/fv3Y9SoUdDT00Pv3r1Rq1YtnD9/Hl26dEGjRo1Urundu3cYPHgwTpw4gRo1aqBPnz6oWrUqnjx5goiICFhbW+Orr77C4MGDAQCRkZHo0aOH3IejiYlJvtvYtGkTZsyYAUNDQ7i7u6NGjRqIjo7G2rVrcejQIfz999+wsLBQWG7Dhg34+++/0aNHD7Ru3Rr//PMP9u3bh6tXryIyMhL6+voqP8/iKEwd0dHR6N27NyQSCdzc3GBtbY1r167Bzc0N7dq1y3P9165dQ58+ffDy5Ut07NgRPXr0QFJSEg4dOoTevXtj5cqVGDZsmMJyFy5cwPLly+Hs7Izhw4fj+fPn0NXVVfo85syZg23btuHx48eYMGGC7HXL6/WbP38+Tp06hW7duqFDhw6IiIhAcHAw7ty5g8OHD6ul/vx06dIFHTp0QFhYGI4cOYJu3boVavnyYODAgViwYAH+/PNPeHt7K7x2O3bswNu3bzF16lSFeQcOHMCJEyfg5uaGtm3bIjY2Fn/99RciIiJw7NgxWFtby9o+fPgQ3bt3x7Nnz9C6dWt4eHggPj4ef/31F06cOIHffvsNw4cPB/DhC01qaioOHz6M1q1bo02bNrL1WFpaytUQExODlStXwsnJCcOHD8eTJ08QGhoKd3d3hIeHo0GDBrK2WVlZGDp0KI4cOQIbGxv07dsX+vr6iIyMxIIFC3DmzBns3bsXFSp8+Eg9duwYPD09YWRkhO7du6N27dpISUnB3bt3sXbtWvz888+ytnkp7vIkjz1Fcm7cuCH71tysWTO5eefPn8eMGTMwb948uek7d+7Ehg0b4Orqij/++ENuHM3SpUuxePFiiMVi/PLLLwCA9PR0TJs2DQKBAIcOHULz5s1l7X/88UesWrVK5Xr9/f1x4sQJtG/fHtu2bZMFL+DDkZ7cIy5DhgzBo0ePEBkZiZ49e2LIkCEqrf/Ro0eYM2cOKlWqhBMnTuDLL7+UzVu0aBECAgIwc+ZM7Nq1S2HZU6dO4cyZM7Czs5NNGzNmDPbs2YNDhw7Bw8ND5edZHKrWIZVKMXnyZGRkZGDz5s1wd3eXtf/jjz/g7e2tsO7s7GyMGDECqampOHDggNwHS3x8PDp16gRvb2+4uLjA3NxcbtmwsDD89ttvGDlypErPw9fXF2fPnsXjx4/h5eUFKysrpW0vXbqEc+fOoXbt2gA+fFC5ubkhKioK//zzj2yfK079BVm4cCHatWuHefPmoXPnzhr/wXT27FmlA3gfPXqkMM3IyAienp5Yv349jhw5onA06n//+x8qVKggCyEfO3LkCHbu3AkXFxfZtFWrVuHHH3+Et7c3QkJCZNOnT5+OZ8+ewcfHBz4+PrLpkydPRufOneHt7Y2OHTuiTp06cHV1lQWdNm3aKJya/9jRo0exdu1aeHp6yqZt2rQJ06dPx9q1a7F8+XLZ9F9//RVHjhzB2LFj4efnJwtuOTk5mD59OjZv3oz169djwoQJAIDg4GBIpVIcOHAATZs2ldtuUlJSgftCcZcneTx19ZkTi8UQi8VYtGgRxo4diw4dOiAjIwOurq6yU1i5zM3NZYOXP/b7779DV1cXq1atUhgsPGPGDJiamsoFgcOHDyM5ORkeHh5yIQcAZs+eDWNjY5Vqz87Oxvr166Gvr48VK1bIhRwA0NXVRY0aNVRalzK7du3Cu3fvMHr0aLmQAwDe3t6oWbMmjh07hmfPniksO378eLlwAQAjRowAAPz777/FqqswVK0jOjoacXFxcHJykgs5ADB69GjUr19fYd3Hjh3DnTt3MHr0aLmQAAA1atTAlClTkJmZif379yss26hRI5VDTmHNnj1bFnIAoEKFChg6dCgA+edcnPoL0qhRIwwdOhS3b9/Gpk2bivhMSk9kZCT8/f3z/Ld9+/Y8lxkzZgwAKDy/8+fP4/r16+jWrRtq1aqlsFy7du3kQg7w4bRTnTp1cOrUKdn76enTpzh16hRq1aqFGTNmyLVv2LAhvv32W0gkEuzcubPQz/frr7+WCzkAMHToUFSoUEFuH8nJycGaNWtgZmYGsVgsd3RKR0cHCxYsgEAgkKtBR+fDR+unf5MAoFq1agXWVtzlSR5j4Wcud9yKQCCAkZERmjZtiv79++f5AdSoUSOF0y0ZGRmIjY1F1apVsWbNmjy3UbFiRTx//hxJSUmoVq0arly5AgAKQQr48C2xSZMmCufe83L79m2kpqaiadOm+X67L47cWvM6baOvr49WrVph3759iI2NVfiD7uDgoLBM7odvfmNR1E3VOvJ7XXR0dNCqVSvcu3dPbnru2J0nT57keTQgt/3t27cV5n0actVJ1edcnPpVMXfuXISEhMDPzw8DBgwo8DRpWZozZ47SIyARERFwc3NTmG5nZ4c2bdrg1KlTePDgAerWrQvg/4PP6NGj81xfXvtYhQoV4OTkhCdPnsjeT7ljAFu1aoWKFSsqLNO+fXsEBgbK9t3CyGsf0dPTg7m5udw+cufOHSQmJqJevXpYunRpnusyNDREXFyc7PGAAQMQGhqKTp06oU+fPmjbti1atGih8t+p4i5P8hh0PnOF+cDN69B9cnIypFIpkpKSZKFJmfT0dFSrVk02ONnMzEzl7eQlNTUVAPL8xqguubUqq0kkEsm1+1heR6Zyvw2qep+i3G92OTk5StvkzsttW9Q6ivK6JCUlAQBCQ0MRGhqqtMY3b96otD51UfU5F6d+VYhEIkydOhW//PILli1bhgULFhRpPZps7NixOHv2LDZv3oz58+cjOTkZ+/fvR/369dG+ffs8l1H22ufue7n7YnHefwVRduRYV1c3z33k/v37Bf6Ny+Xq6oq9e/di1apV2L59u+xWHfb29pgzZ47CEVN1L0/yGHRIZQKBQGFa7h8Le3t7REVFqbSe3GXyumIJAF68eKHSenK/HSu7MZs65NaqrKaEhAS5diW1/eTkZKVtcv8QF/doQVFel9xlgoOD0atXr0JtL6/9qbQVp35VTZkyBZs3b8batWuVHuEAPvSHsgCcG+o1Uc+ePVGrVi1s2bIFvr6+2LZtGzIzMzFy5Eilr7Gy91Puvpf7upT1++/jdXfr1g07duxQeblOnTqhU6dOyMjIwKVLl3DixAls2LABI0eOVBgPVhLL0//jGB0qlipVqsDe3h5xcXFITExUaZncwXW5g54/9vr1a7lL1vPzxRdfwMTEBDdu3MDjx48LbF/Yoykf15rXpakSiUR26uPTAYPqknsFmrJLXj+eV5ir1fKS3+uSk5OTZw0tWrQAAJw7d65Y21bVx4NA1aE06jc0NMQPP/wAiUSCn3/+WWk7oVCIJ0+e5Dkv9zYJmqhChQoYMWIEXr58iYMHD2Lz5s3Q19fPd8B/XvtYVlaW7P2Ue++h3P9GR0fneQl37uX/H5+GKsr7PD+5f2cuXbqUZw0FMTQ0RJs2bfDTTz9h4cKFkEqlClf+leTyxKBDajBp0iS8f/8eEydOzPPIw+vXr/HPP//IHvfo0QNCoRAhISFy0wFgyZIlKh+G1tXVxdixYyGRSDBt2jRkZGTIzc/OzkZ8fLzsce5NuZR9mORlwIABqFixIjZs2KAwTmP58uV49uwZunbtipo1a6q8zsJwdnZGvXr18N9//yE4OFhhfmxsLLZs2YIKFSooDKwsLCcnJ9ja2iI6Olph8O2GDRsUxucAH17L+vXrY9OmTUr/+F65ckV21Km4cl9DVYKtKkqr/oEDB8LBwQEhISGIiYnJs02LFi3w5MkTHDt2TG765s2bC7yPUVkbOXIk9PT08P333+P27dtwd3fP9yZ44eHhOHr0qNy0oKAgPHnyBB06dJCdjq5duzY6deqEp0+fYsWKFXLtb9y4gY0bN0JfXx8DBgyQTS/K+zw/FSpUwIQJE/Dy5UvMmjULb9++VWiTmJgo9wXt9OnTebbLPQJlYGCQ7zaLuzzJ46krKrYhQ4bgypUrWLduHRwcHNCpUydYWloiNTUVjx49QlRUFDp06IBt27YB+HAUaMWKFRg1ahR69uyJPn36oFatWjh37hyuX78OZ2dnlU+DzZ49G5cvX8bJkyfRrFkzdOvWDVWrVsWzZ88QERGBoUOHygZYfvPNN9DR0cGaNWuQnJwsO+8/btw4pad9LC0t4e/vjxkzZqBDhw7o3bs3RCIRoqOjERkZidq1a2PZsmVq6MW86erqYt26dejbty+mTp2K7du3o3nz5tDT08Pt27dx9OhRZGdnY8mSJahXr16xtiUQCLBq1Sr06dMHo0aNkruPTlhYGDp37owTJ07ILaOnp4ctW7bAw8MDgwcPRvPmzdG0aVNUrlwZT58+RWxsLOLi4hAeHq6Wq0U6dOiAffv24bvvvoO7uzsqV64MExMTjBs3rkjrK636BQIBFi1aBFdX1zwDIwBMnToVJ06cwNChQ9G7d2+YmZkhJiYGMTExcHFxUQgGmkQkEsHV1RX79u0D8OH3vPLTvXt3DBkyBL169ULdunURGxuLEydOoFq1aggICJBru3z5cnTr1g2LFy9GeHg4WrRoIbuPTkZGBlasWCH38wktW7ZElSpVEBISgooVK6JOnToQCATw9PRUuJeOqry9vXH9+nUEBwfj2LFjaNeuHWrXro1Xr17h/v37OH/+PMaMGSM7AvXDDz/g0aNHaN26NSwtLWFgYIBr167h5MmTqFatmuyqR2WKuzzJY9AhtViyZAm6du2KDRs24OzZs0hOToaJiQlq1aqF0aNHo3///nLt3d3dsXfvXvj7+2P//v2oWLEinJ2dcfz4cfz6668qB52KFSti165d2Lx5M7Zv347du3cjKysLIpEIrVu3Rvfu3WVtbWxssGHDBqxYsQJbtmyRHQEq6GqYUaNGoX79+li1ahUOHTqEN2/eoGbNmhg3bhxmzZpVooNqgQ/f9M+ePYvVq1cjLCwM69evR3Z2NszMzNCrVy+MHz8eLVu2VMu2WrVqhb///hsLFy7EyZMncfLkSTg6OuLgwYM4efKkQtABPozPioyMRFBQEA4fPozt27dDKpVCJBLBzs4OU6ZMga2trVrqGzp0KJ4+fYpdu3YhMDAQ79+/h4WFRZGDTmnW36ZNG/To0UPpkaM2bdpg586d8PPzQ2hoqNx7Yv/+/RoddIAPr82+fftgb2+PVq1a5dvW1dUVI0eOREBAAI4cOQI9PT24u7tj/vz5CrcxsLKywunTp2Vtz58/j8qVK6N169aYOnUq2rZtK9fexMQEW7duhVgsRkhICNLT0wF82LeLGnQqVKiA4OBg7N27F1u3bsXx48dlF1dYWFhg+vTpGDhwoKz9zJkzcejQIVy+fFl22rtWrVrw8vLCxIkTC/xdq+IuT/IEKSkpvJc0EREVy7Jly7Bw4UIEBATI7q/zKbFYDH9/fwQGBqp8006i4uIYHSIiKpb09HT88ccfMDY2LvZYMSJ146krIiIqkr///huXL1/G8ePHER8fj/nz58v9ph2RJmDQISKiIgkNDcX27dthbm6OadOmYerUqWVdEpECjtEhIiIircUxOkRERKS1GHSIiIhIazHoEBERkdZi0MlDXFxcWZdQrrH/io99WDzsv+JjHxYf+7B41NV/DDpERESktRh0iIiISGsx6BAREZHWYtAhIiIircWgQ0RERFqLQYeIiIi0FoMOERERaS0GHSIiItJaDDpERESktRh0iIiISGsx6BAREZHWqlDWBRCR5nouEeDpm5yyLqPQDKpUL+sSiEhDMOgQkVJP3+Rg1tmXZV1Goc1rql/WJRCRhuCpKyIiItJaDDpERESktRh0iIiISGsx6BAREZHWYtAhIiIircWgQ0RERFqLl5cTkdYxMDDEP0nSsi6jUGpX1kFN/fJVM1F5wKBDRFon6Z0UAdHl6/4/AW3MUFNfUNZlEGmdMj11FRkZiYEDB+LLL7+EUCjE1q1b5eZLpVKIxWLY2dmhRo0a6NmzJ27cuCHXRiKRwNvbG/Xr10etWrUwcOBAPH36tDSfBhEREWmoMg06b968gb29Pfz8/GBoaKgwf8WKFQgMDIS/vz9OnToFMzMz9OnTB69fv5a18fX1xYEDB7BhwwYcPnwYr1+/hqenJ7Kzs0vzqRAREZEGKtOg07VrV8ybNw/u7u7Q0ZEvRSqVIigoCNOmTYO7uzvs7e0RFBSE9PR07NmzBwCQmpqKP//8EwsWLECHDh3g4OCAtWvX4tq1azh9+nQZPCMiIiLSJBp71dXDhw+RkJCAjh07yqYZGhrC2dkZ0dHRAICYmBi8f/9erk2dOnXQoEEDWRsiIiL6fGnsYOSEhAQAgJmZmdx0MzMzPH/+HADw4sUL6OrqwtTUVKHNixcvlK47Li6uwO2r0oaUY/8Vnyb0YVpFc0gkmWVdRhGYlLu6016nIS5R+d+tsqAJ+2B5xz4sHmX9Z2trq/I6NDbo5BII5K9CkEqlCtM+VVCbgjooLi6uUJ1I8th/xacpfZiaJIW+vqSsyygSfX2Dsi6hUIyNjGFbzaSsy5DRlH2wPGMfFo+6+k9jT12JRCIAUDgy8+rVK9lRHnNzc2RnZyMxMVFpGyIiIvp8aWzQsbKygkgkQlhYmGxaZmYmzp07BycnJwCAg4MD9PT05No8ffoUt27dkrUhIiKiz1eZnrpKT0/HvXv3AAA5OTl48uQJYmNjUbVqVVhYWMDLywvLli2Dra0tbGxsEBAQgMqVK6Nfv34AABMTEwwbNgzz5s2DmZkZqlatirlz56Jhw4Zo3759GT4zIiIi0gRlGnQuX74MNzc32WOxWAyxWIxBgwYhKCgI3333HTIyMuDt7Y2UlBQ4OjoiJCQERkZGsmV++eUX6OrqYtSoUcjMzES7du2wZs0a6OrqlsVTIiIiIg1SpkGnbdu2SElJUTpfIBDA19cXvr6+StsYGBhg6dKlWLp0aQlUSEREROWZxo7RISIiIiouBh0iIiLSWgw6REREpLUYdIiIiEhrMegQERGR1mLQISIiIq3FoENERERai0GHiIiItBaDDhEREWktBh0iIiLSWgw6REREpLUYdIiIiEhrMegQERGR1mLQISIiIq3FoENERERai0GHiIiItBaDDhEREWktBh0iIiLSWgw6REREpLUYdIiIiEhrMegQERGR1mLQISIiIq3FoENERERai0GHiIiItBaDDhEREWktBh0iIiLSWgw6REREpLUYdIiIiEhrMegQERGR1mLQISIiIq3FoENERERai0GHiIiItBaDDhEREWktBh0iIiLSWgw6REREpLUYdIiIiEhrMegQERGR1mLQISIiIq3FoENERERai0GHiIiItBaDDhEREWktBh0iIiLSWhoddLKzs7Fo0SI0adIEIpEITZo0waJFi5CVlSVrI5VKIRaLYWdnhxo1aqBnz564ceNGGVZNREREmkKjg85vv/2G9evXw9/fHxcuXICfnx/++OMPLF++XNZmxYoVCAwMhL+/P06dOgUzMzP06dMHr1+/LsPKiYiISBOoHHQiIyPx6tUrpfMTExMRGRmplqJyXbhwAd26dUP37t1hZWWFHj16oHv37rh06RKAD0dzgoKCMG3aNLi7u8Pe3h5BQUFIT0/Hnj171FoLERERlT8qBx03NzeEhYUpnX/mzBm4ubmppahcrVq1wtmzZ3H79m0AwM2bNxEREYEuXboAAB4+fIiEhAR07NhRtoyhoSGcnZ0RHR2t1lqIiIio/KmgakOpVJrv/Hfv3kFHR71nwqZNm4b09HQ4OTlBV1cXWVlZmDVrFsaMGQMASEhIAACYmZnJLWdmZobnz58rXW9cXFyB21alDSnH/is+TejDtIrmkEgyy7qMIjApd3WnvU5DXOKLsi5Djibsg+Ud+7B4lPWfra2tyuvIN+ikpaUhNTVV9jgpKQmPHz9WaJeSkoK9e/eiZs2aKm9YFSEhIdixYwfWr18POzs7XL16FT4+PrC0tMTw4cNl7QQCgdxyUqlUYdrHCuqguLi4QnUiyWP/FZ+m9GFqkhT6+pKyLqNI9PUNyrqEQjE2MoZtNZOyLkNGU/bB8ox9WDzq6r98g87vv/+OJUuWAPgQJnx9feHr65tnW6lUih9//LHYBX1s3rx5mDx5Mvr27QsAaNiwIR4/foxff/0Vw4cPh0gkAgC8ePECderUkS336tUrhaM8RERE9PnJN+i0b98eBgYGkEqlWLBgATw8PNC4cWO5NgKBAJUqVcJXX32F5s2bq7W4t2/fQldXV26arq4ucnJyAABWVlYQiUQICwtDs2bNAACZmZk4d+4cFixYoNZaiIiIqPzJN+i0atUKrVq1AgBIJBK4ubmhYcOGpVIYAHTr1g2//fYbrKysYGdnh9jYWAQGBmLgwIEAPoQsLy8vLFu2DLa2trCxsUFAQAAqV66Mfv36lVqdREREpJlUHozs4+NTknXkacmSJVi8eDFmzpyJV69eQSQSYcSIEZg9e7aszXfffYeMjAx4e3sjJSUFjo6OCAkJgZGRUanXS0RERJpFadDZvn07AGDgwIEQCASyxwUZNGiQeioDYGRkBD8/P/j5+SltU9DYISIiIvp8KQ06EydOhEAgQN++fVGxYkVMnDixwJUJBAK1Bh0iIiKi4lAadK5cuQIAqFixotxjIiIiovJCadCxtLTM9zERERGRptPoH/UkIiIiKg6Vr7oCgNOnT2Pz5s148OABkpOTFX4WQiAQICYmRp31ERERERWZykEnKCgIc+fORfXq1dG8eXN8+eWXJVkXERERUbGpHHQCAwPRunVr7N27VzZAmYiIiEiTqTxGJzExER4eHgw5REREVG6oHHQcHBzw6NGjkqyFiIiISK1UDjqLFy/Gtm3bEB4eXpL1EBEREamNymN0xGIxjI2N0bt3b1hbW8PCwkLhl8UFAgF27dql9iKJiIiIikLloHPz5k0IBALUqVMHEokEd+7cUWgjEAjUWhwRERFRcagcdK5evVqSdRARERGpHe+MTERERFpL5SM6jx8/VqmdhYVFkYshIiIiUieVg06TJk1UGoOTlJRUrIKIiIiI1EXloLN69WqFoJOdnY2HDx9ix44dMDc3x5gxY9ReIBEREVFRqRx0hgwZonTetGnT0LFjR6Snp6ulKCIiIiJ1UMtg5CpVqmDIkCH4/fff1bE6IiIiIrVQ21VXenp6eP78ubpWR0RERFRsagk6V69exZo1a9CgQQN1rI6IiIhILYp91VVqairS0tJQpUoVBAYGqrU4IiIiouJQOei0bt1aIegIBAIIhULUr18fffv2hVAoVHd9REREREWmctAJCgoqyTqIiIiI1I4/AUFERERai0GHiIiItBaDDhEREWktBh0iIiLSWgw6REREpLVUCjqZmZnw9/fHqVOnSroeIiIiIrVRKegYGBjg119/xZMnT0q6HiIiIiK1UfnUVePGjXHv3r2SrIWIiIhIrVQOOvPmzUNwcDCOHj1akvUQERERqY3Kd0ZeuXIlhEIhBg0ahFq1aqFu3bowNDSUayMQCLBr1y61F0lERERUFCoHnZs3b0IgEKBOnToAgEePHim0yetHP4mIiIjKispB5+rVqyVZBxEREZHa8T46REREpLUKFXSys7Oxa9cuTJ48GZ6envjvv/8AACkpKdi3bx/i4+NLpEgiIiKiolA56KSmpqJr164YP3489u/fj+PHjyMxMREAYGRkhLlz52LdunUlVigRERFRYakcdH7++WfcvHkTu3fvRkxMDKRSqWyerq4u3NzccPz48RIpkoiIiKgoVA46hw4dwrhx49C5c+c8r66ytrbG48eP1VocERERUXGoHHRSUlJQr149pfOlUinevXunlqI+Fh8fjwkTJsDa2hoikQhOTk44e/as3HbFYjHs7OxQo0YN9OzZEzdu3FB7HURERFT+qBx0LC0tcf36daXzIyMjYWNjo5aicqWkpMDFxQVSqRS7du1CdHQ0lixZAjMzM1mbFStWIDAwUPajo2ZmZujTpw9ev36t1lqIiIio/FE56PTv3x/BwcGIjIyUTcs9hbV27VocPHgQgwcPVmtxK1euRI0aNbB27Vo4Ojqibt26+Oabb9CgQQMAH47mBAUFYdq0aXB3d4e9vT2CgoKQnp6OPXv2qLUWIiIiKn9UDjrTp09HmzZt0KtXL3Tv3h0CgQA+Pj6ws7ODj48PevTogfHjx6u1uEOHDsHR0RGjRo2CjY0N2rRpg3Xr1skGQj98+BAJCQno2LGjbBlDQ0M4OzsjOjparbUQERFR+aPynZH19PSwa9cu7N69G3/99RcEAgGysrLQtGlTeHh4YMCAAWr/CYgHDx5gw4YNmDhxIqZNm4arV69izpw5AIBx48YhISEBAOROZeU+fv78udL1xsXFFbhtVdqQcuy/4tOEPkyraA6JJLOsyygCk3JXd9rrNMQlvijrMuRowj5Y3rEPi0dZ/9na2qq8DpWDTq7+/fujf//+hV2sSHJycvDVV19h/vz5AICmTZvi3r17WL9+PcaNGydr92nAkkql+YaugjooLi6uUJ1I8th/xacpfZiaJIW+vqSsyygSfX2Dsi6hUIyNjGFbzaSsy5DRlH2wPGMfFo+6+q/QQQcA/vvvP9ml5BYWFmjYsGGJ/KCnSCSSjcfJ9cUXX+DJkyey+QDw4sUL2Y+NAsCrV68UjvIQERHR56dQPwGxd+9eNGrUCO3atcPgwYMxePBgtGvXDo0bN8bu3bvVXlyrVq1w584duWl37tyBhYUFAMDKygoikQhhYWGy+ZmZmTh37hycnJzUXg8RERGVLyof0dm6dSsmT54MW1tb/Pzzz7CxsYFUKsXdu3cRHByM8ePH4927dxgyZIjaips4cSK6du2KgIAAeHh4IDY2FuvWrcOPP/4I4MMpKy8vLyxbtgy2trawsbFBQEAAKleujH79+qmtDiIiIiqfVA46y5cvh6OjIw4ePAgDA/lz32PHjkWPHj2wfPlytQadZs2aYevWrViwYAGWLl2KOnXq4Pvvv8eYMWNkbb777jtkZGTA29sbKSkpcHR0REhICIyMjNRWBxEREZVPKgedp0+fYty4cQohBwAMDAzg6emJn376SZ21AQBcXFzg4uKidL5AIICvry98fX3Vvm0iIiIq31Qeo2NnZ5fvJdvPnj1TGDhMREREVJZUDjoLFizA5s2bsW/fPoV5e/fuRXBwMBYuXKjW4oiIiIiKQ+VTV6tWrYKpqSlGjx4NHx8f1KtXDwKBAPfu3cPLly9hbW2NlStXYuXKlbJlBAIBdu3aVSKFExERERVE5aBz8+ZNCAQC2f1qnj17BgDQ19dHnTp1IJFIcOvWLbllSuLeOkRERESqUjnoXL16tSTrICIiIlK7Qt0wkIiIiKg8YdAhIiIircWgQ0RERFqrSD/qSUSF91wiwNM3OSq1TatojtQkaQlXVLAM1colNRDo6OCfJM3pcFX3wdqVdVBTv+z3VSJlGHSISsnTNzmYdfalSm0lkkzo60tKuKKC+baoXtYlfDZeZWRDfPFVWZcho+o+GNDGDDX1eYUtaS6euiIiIiKtpXLQadq0KQ4fPqx0/pEjR9C0aVO1FEVERESkDioHnUePHuHNmzdK57958waPHz9WS1FERERE6lCoU1f53en4zp07MDIyKnZBREREROqS72Dkbdu2Yfv27bLHAQEB2Lx5s0K7lJQUXL9+HS4uLuqvkIiIiKiI8g06b968QUJCguxxamoqcnLkL38UCASoVKkSRowYAR8fn5KpkoiIiKgI8g06Y8eOxdixYwEATZo0gZ+fH3r06FEqhREREREVl8r30YmNjS3JOoiIiIjUrtA3DHz9+jWePHmC5ORkSKWKd8Ns3bq1WgojIiIiKi6Vg05ycjLmzJmDffv2ITs7W2G+VCqFQCBAUlKSWgskIiIiKiqVg8706dNx8OBBjB07Fq1bt4ZQKCzBsoiIiIiKT+Wgc+LECYwfPx6LFy8uyXqIiIiI1EblGwZWrFgR1tbWJVkLERERkVqpHHTc3d1x/PjxkqyFiIiISK1UDjpTpkxBfHw8JkyYgIsXLyI+Ph4vX75U+EdERESkKVQeo+Po6AiBQICYmBjs2rVLaTtedUVERESaQuWgM3v27Hx/1JOIiIhI06gcdHx9fUuyDiIiIiK1U3mMzseys7ORlJSErKwsdddDREREpDaFCjr//vsvevfujVq1asHGxgaRkZEAgMTERAwYMABnzpwpkSKJiIiIikLloHPhwgX06NED9+/fx8CBA+V+58rU1BTp6en4888/S6RIIiIioqJQOegsXLgQ1tbWiI6Oxrx58xTmt23bFv/8849aiyMiIiIqDpWDzr///ouhQ4fCwMAgz6uvateujYSEBLUWR0RERFQcKgcdHR0d6Ogob56QkABDQ0O1FEVERESkDioHHQcHBxw5ciTPee/evcPu3bvRsmVLtRVGREREVFwqB50ZM2YgPDwckydPxtWrVwEA8fHxOHHiBHr16oX79+9j5syZJVYoERERUWGpfMPADh06YO3atfD29sa2bdsAAF5eXpBKpTAxMcH69evRokWLEiuUiIiIqLBUDjoA0K9fP/To0QNhYWG4e/cucnJyUK9ePXTq1AlVqlQpqRqJiIiIiqRQQQcAKlWqhJ49e5ZELURERERqpfIYncOHD8Pb21vpfG9vb6WDlYmIiIjKgspBZ9WqVXj79q3S+ZmZmVixYoVaiiIiIiJSB5WDzvXr1+Hg4KB0ftOmTXHz5k111ERERESkFioHnaysLGRkZCidn5GRAYlEopailFm2bBmEQqHcKTSpVAqxWAw7OzvUqFEDPXv2xI0bN0q0DiIiIiofVA469vb2CA0NRU5OjsK8nJwchIaGws7OTq3FfezixYvYvHkzGjZsKDd9xYoVCAwMhL+/P06dOgUzMzP06dMHr1+/LrFaiIiIqHxQOehMmDABly5dwqBBgxATEwOJRAKJRIKYmBgMHjwYly5dwvjx40ukyNTUVIwdOxarVq2CUCiUTZdKpQgKCsK0adPg7u4Oe3t7BAUFIT09HXv27CmRWoiIiKj8UPny8r59++L+/fsQi8U4fvw4AEAgEEAqlUIgEGDOnDnw9PQskSJzg8w333yDJUuWyKY/fPgQCQkJ6Nixo2yaoaEhnJ2dER0djVGjRpVIPURERFQ+FOo+OrNmzUK/fv1w4MABPHjwAFKpFPXq1YObmxvq1q1bIgVu3rwZ9+7dw9q1axXm5f5aupmZmdx0MzMzPH/+XOk64+LiCtyuKm1IOfaforSK5pBIMlVuX5i2JeV91nuNqKPwTMpd3ZrY16rUk/Y6DXGJL0qhmvKJfwuLR1n/2draqrwOlYJORkYGBgwYAE9PTwwdOhRTpkxReQPFERcXhwULFuDvv/9GxYoVlbYTCARyj3OPMilTUAfFxcUVqhNJHvsvb6lJUujrqzZgXyLJhL6+QQlXVDC9CnoaUUdRlLe6Na2vVd0HjY2MYVvNpBQqKn/4t7B41NV/Ko3RMTQ0xJUrV5CdnV3sDRbGhQsXkJiYiK+//hqmpqYwNTVFZGQk1q9fD1NTU1SrVg0A8OKF/LeJV69eKRzlISIios+PyoOR27Rpg6ioqJKsRUHPnj0RFRWFiIgI2b+vvvoKffv2RUREBGxsbCASiRAWFiZbJjMzE+fOnYOTk1Op1kpERESaR+UxOv7+/vDw8MCPP/6I0aNHw9LSEjo6KuekIhEKhXJXWQEffmuratWqsLe3B/DhF9SXLVsGW1tb2NjYICAgAJUrV0a/fv1KtDYiIiLSfCoHnRYtWkAqlSIwMBCBgYHQ0dGBnp6eXBuBQIBnz56pvcj8fPfdd8jIyIC3tzdSUlLg6OiIkJAQGBkZlWodREREpHlUDjp9+vTJd4BvaTl06JDcY4FAAF9fX/j6+pZRRURERKSpVA46QUFBJVkHERERkdqV7CAbIiIiojJUqKDz6NEjTJ06FQ4ODrCwsMDZs2cBAImJiZg5cyZiYmJKokYiIiKiIlH51NWtW7fQrVs35OTkoHnz5nj06JHsvjqmpqa4ePEiJBIJVq9eXWLFEhERERWGykFn/vz5MDIywokTJ6CrqwsbGxu5+V27dsVff/2l7vqIiIiIikzlU1dRUVEYM2YMzM3N87z6ysLCIt/flyIiIiIqbSoHnaysLFSuXFnp/OTkZOjq6qqlKCIiIiJ1UDno2NvbIyIiIs95UqkUBw4cgIODg7rqIiIiIio2lYOOl5cX9u/fjyVLliApKQkAkJOTg9u3b+Pbb7/F5cuXS+1XzYmIiIhUofJg5L59++Lx48dYvHgx/Pz8ZNMAQFdXF4sWLUKXLl1KpkoiIiKiIlA56ADAtGnT0K9fP4SGhuLevXvIyclBvXr10KtXL1hZWZVUjURERERFUmDQkUgkOHz4MB48eIBq1arBxcUFEydOLI3aiIiIiIol36CTkJCAHj164P79+5BKpQCAypUrY+fOnWjdunWpFEhERERUVPkORl60aBEePHiAiRMnYufOnRCLxdDX18fs2bNLqz4iIiKiIsv3iM6pU6cwaNAgLFq0SDbN3NwcY8aMwdOnT1G7du0SL5CIiIioqPI9opOQkAAnJye5aa1atYJUKsWTJ09KtDAiIiKi4so36GRnZ8PAwEBuWu7jzMzMkquKiIiISA0KvOrqwYMHuHTpkuxxWloaACAuLg5VqlRRaO/o6KjG8oiIiIiKrsCgIxaLIRaLFaZ/OiBZKpVCIBDI7ppMREREVNbyDTqBgYGlVQcRERGR2uUbdAYPHlxadRARERGpnco/6klERERU3jDoEBERkdZi0CEiIiKtxaBDREREWotBh4iIiLQWgw4RERFpLQYdIiIi0loMOkRERKS1GHSIiIhIazHoEBERkdZi0CEiIiKtxaBDREREWotBh4iIiLQWgw4RERFpLQYdIiIi0loMOkRERKS1GHSIiIhIazHoEBERkdZi0CEiIiKtxaBDREREWotBh4iIiLSWRged5cuXo0OHDrCwsIC1tTU8PT1x/fp1uTZSqRRisRh2dnaoUaMGevbsiRs3bpRRxURERKRJNDronD17FqNHj8bRo0cRGhqKChUqoHfv3khOTpa1WbFiBQIDA+Hv749Tp07BzMwMffr0wevXr8uwciIiItIEFcq6gPyEhITIPV67di0sLS1x/vx5dO/eHVKpFEFBQZg2bRrc3d0BAEFBQbC1tcWePXswatSosiibiIiINIRGH9H5VHp6OnJyciAUCgEADx8+REJCAjp27ChrY2hoCGdnZ0RHR5dRlURERKQpNPqIzqd8fHzQuHFjtGzZEgCQkJAAADAzM5NrZ2ZmhufPnytdT1xcXIHbUqUNKcf+U5RW0RwSSabK7QvTtqS8z3qvEXUUnkm5q1sT+1qVetJepyEu8UUpVFM+8W9h8SjrP1tbW5XXUW6Czvfff4/z58/jyJEj0NXVlZsnEAjkHkulUoVpHyuog+Li4grViSSP/Ze31CQp9PUlKrWVSDKhr29QwhUVTK+CnkbUURTlrW5N62tV90FjI2PYVjMphYrKH/4tLB519V+5OHXl6+uLvXv3IjQ0FHXr1pVNF4lEAIAXL+S/Tbx69UrhKA8RERF9fjQ+6MyZMwd79uxBaGgovvjiC7l5VlZWEIlECAsLk03LzMzEuXPn4OTkVNqlEhERkYbR6FNXs2bNws6dO7FlyxYIhULZmJzKlSujSpUqEAgE8PLywrJly2BrawsbGxsEBASgcuXK6NevXxlXT0RERGVNo4PO+vXrAUB26XiuOXPmwNfXFwDw3XffISMjA97e3khJSYGjoyNCQkJgZGRU6vUSERGRZtHooJOSklJgG4FAAF9fX1nwISIiIsql8WN0iIiIiIqKQYeIiIi0lkafuiIiIs0m0NHBP0k5ZV1GodWurIOa+tKyLoNKAYMOEREV2auMbIgvvirrMgotoI0Zauorv7EsaQ+euiIiIiKtxSM6VO48lwjw9E35O1SeUf5KJiIq9xh0qNx5+iYHs86+LOsyCs23RfWyLoGI6LPDU1dERESktRh0iIiISGsx6BAREZHWYtAhIiIircWgQ0RERFqLQYeIiIi0FoMOERERaS0GHSIiItJaDDpERESktRh0iIiISGsx6BAREZHWYtAhIiIircWgQ0RERFqLQYeIiIi0FoMOERERaS0GHSIiItJaDDpERESktRh0iIiISGsx6BAREZHWYtAhIiIircWgQ0RERFqLQYeIiIi0FoMOERERaS0GHSIiItJaDDpERESktRh0iIiISGsx6BAREZHWYtAhIiIircWgQ0RERFqLQYeIiIi0FoMOERERaa0KZV0AERFRaRPo6OCfpJwS3UZaRXOkJknVtr7alXVQU1996/tcMOgQEdFn51VGNsQXX5XoNiSSTOjrS9S2voA2ZqipL1Db+j4XDDqfuecSAZ6+Ue+3GnV/i/lURsl+CSMiIi2iNUFn/fr1WLlyJRISEmBnZwexWAxnZ+eyLkvjPX2Tg1lnX6p1ner+FvMp3xbVS2zdRESkXbRiMHJISAh8fHwwc+ZMhIeHo2XLlujfvz8eP35c1qVRKRHk5KDTxaP4ftM8dLp4FIIcHvYhIiItCTqBgYEYPHgwRowYgQYNGmDp0qUQiUTYuHFjWZdGpUCQk4OVy8dh0ZpZ6B+2A4vWzMLK5eMYdoiIqPyfunr37h1iYmIwZcoUuekdO3ZEdHR0kdZpa2tb6GVKYqxLaSiJ8S76+gbqX2k+Ol46juY3omHw/sPpMoP3EjjeuIAOl47jVAuXUq1FXUq7D7WNXgW9si6h3OM+WHzsw+IpymdxXsp90ElMTER2djbMzMzkppuZmeHFixelVkdNfWm5HQ1/upd5WZdQaC4f1WwQ9i/038uPCTJ8n4lf3lxGZq9hpV1avlzKeV+XJ+Wx7vJYM8C6SbNpxakrABAI5EOGVCpVmEbaKat9e0gN5L85SQ0MkNW+fdkUREREGqPcBx1TU1Po6uoqHL159eqVwlEe0k5Zbm7IatNGFnakBgbIatMGWa6uZVwZERGVtXIfdCpWrAgHBweEhYXJTQ8LC4OTk1MZVUWlSkcHb3ftwtt16yD59lu8XbcOb3ftAnTK/e5NRETFVO7H6ADApEmTMH78eDg6OsLJyQkbN25EfHw8Ro0aVdalUWnR0UFWr17I6tWrrCshIiINohVfeT08PCAWi7F06VK0bdsW58+fx65du2BpaVngsv/73//g6uoKS0tLCIVCPHz4sMBltm7dCqFQqPAvMzNTHU+n3ClKHwLA/v374eTkBHNzczg5OeHAgQMlXKlmkkgk8Pb2Rv369VGrVi0MHDgQT58+zXeZz30fXL9+PZo0aQKRSIRvvvkGUVFR+ba/du0aevTogRo1auDLL7+Ev78/pNLP+zeDCtOHDx8+zHN/O3HiRClWrDkiIyMxcOBAfPnllxAKhdi6dWuBy3AflFfYPizOPqgVQQcAxowZg6tXr+LFixc4c+YMWrdurdJyb9++RceOHeHj41Oo7VWqVAm3bt2S+2dg8HleSliUPrxw4QK+/fZb9O/fHxEREejfvz9GjhyJf/75pwQr1Uy+vr44cOAANmzYgMOHD+P169fw9PREdnZ2vst9rvtgYW8QmpaWhj59+sDc3BynTp2Cn58fVq1ahdWrV5dy5ZqjqDdZ3bt3r9z+1q5du1KqWLO8efMG9vb28PPzg6GhYYHtuQ8qKmwf5irKPqgVp66KY+LEiQCAy5cvF2o5gUAAkUhUEiWVO0Xpw6CgILRt2xazZs0CADRo0AAREREICgrChg0bSqROTZSamoo///wTgYGB6NChAwBg7dq1aNy4MU6fPo1OnTopXfZz3Qc/vkEoACxduhQnT57Exo0bMX/+fIX2u3fvRkZGBoKCgmBoaAh7e3vcvn0bv//+OyZPnvxZXp1Z2D7MVa1atc9yn/tU165d0bVrVwD///cvP9wHFRW2D3MVZR/UmiM6pS0jIwONGjWCvb09PD09ceXKlbIuqVy5ePEiOnbsKDetU6dORb7JY3kVExOD9+/fy/VFnTp10KBBgwL74nPcB3NvEPrpvpPfDUIvXLiAr7/+Wu5bY6dOnfD8+XOVT7Nqk6L0Ya5hw4bBxsYGLi4u2L9/f0mWqVW4D6pPUfZBBp0isLW1xerVq7Ft2zasX78e+vr66NatG+7evVvWpZUbCQkJZX6TR03w4sUL6OrqwtTUVG56QX3xue6DRblB6IsXL/Jsnzvvc1OUPqxSpQoWLlyITZs2Yffu3WjXrh1GjRqFnTt3lkbJ5R73weIrzj6olaeuFi1ahICAgHzbHDhwAG3bti3S+lu2bImWLVvKHjs5OaFt27ZYu3YtlixZUqR1apqS7kNAu2/yqGr/KVNQX3wO+2B+Crvv5NU+r+mfk8L0oampqdzP7Hz11VdISkrCihUr4OnpWaJ1agvug8VTnH1QK4OOl5cXBgwYkG+bOnXqqG17urq6cHBwwL1799S2zrJW0n0oEom0+iaPqvbfxYsXkZ2djcTERFSvXl0279WrV3B2dlZ5e9q4D+alKDcINTc3z7M9AK3Z3wpDXTdZdXR0VOlqI+I+WFJU3Qe1MuiYmpoqnAooSVKpFNeuXUOjRo1KbZslraT7sEWLFggLC8PUqVNl07TpJo+q9p+DgwP09PQQFhaG/v37AwCePn2KW7duFaovtHEfzMvHNwjt3bu3bHpYWBh6KbmHUsuWLfHTTz8hMzNTdlVaWFgYatasCSsrq9IoW6MUpQ/zcvXqVQ5MVhH3wZKh6j742Y/RSUhIQGxsLO7cuQMAuHXrFmJjY5GcnCxr06tXL/z888+yx35+fjh58iQePHiA2NhYTJ48GdeuXcO3335b6vVrgqL04YQJExAeHo7ly5fj9u3bWL58OSIiIuDl5VXq9ZclExMTDBs2DPPmzcPp06dx5coVjB8/Hg0bNkT7j36ri/vg/5s0aRK2bduG4OBg3Lp1C3PmzJG7QejPP/8s94Hdr18/GBoaYuLEibh+/TpCQ0Px22+/YeLEiZ/taYPC9uG2bduwe/du3Lp1C3FxcVi1ahXWr1+PcePGldVTKFPp6emIjY1FbGwscnJy8OTJE8TGxsouz+c+WLDC9mFx9kGtPKJTGBs3boS/v7/sce7phsDAQAwZMgQAcP/+fdSuXVvWJjU1Fd999x1evHgBY2NjNGnSBIcPH4ajo2PpFq8hitKHuXewXrRoEcRiMerVq4eNGzeiefPmpVu8Bvjll1+gq6uLUaNGITMzE+3atcOaNWugq6sra8N98P95eHggKSkJS5cuRUJCAr788ku5G4TGx8fj/v37svYmJibYt28fZs2ahQ4dOkAoFGLSpEmYPHlyWT2FMlfYPgSAgIAAPH78GLq6urC2tsbq1as/2/E5ly9fhpubm+yxWCyGWCzGoEGDEBQUxH1QBYXtQ6Do+6AgJSXl8701IxEREWm1z/7UFREREWkvBh0iIiLSWgw6REREpLUYdIiIiEhrMegQERGR1mLQISIiIq3FoENaSywWQygUKkz//fff8dVXX8HU1BSNGzcucDoREZVfDDpULmzduhVCoVD2TyQSwc7ODh4eHlizZg1ev36t0nrOnj2L77//Hk2bNsWqVasgFovznf45a9y4sVyfm5ubw8HBAd7e3khISCjSOtPT0yEWixEREaHmajXTqVOnMHToUNjZ2cHMzAyWlpbo1KkTFi9ejGfPnpV1eUUWEREh2y+2bNmSZ5uRI0fK3qtEZemzvzMylS8+Pj6oV68e3r9/jxcvXuDs2bPw9fVFYGAgtm/fLvdbT97e3pg+fbrc8rkfsL/99pvc0R5l0z93DRs2lP0emUQiwdWrV7F582aEh4cjKipK7u7Nqnjz5o3sLtpt27ZVe72aQiqVYubMmdi4cSPs7OwwfPhwWFhYICMjAzExMVi7di02bdok+9mU8srAwAC7d+/G0KFD5aanpaXhyJEjMDAwkP1KN1FZYdChcqVTp05o0aKF7PGMGTNw5swZDBw4EIMGDcKFCxdgaGgIAKhQoQIqVJDfxXN/MfjTMKNsenG8ffsWlSpVUtv6ykKNGjUUbrFuYmKCgIAAXL16FQ4ODmVTmIYLDAzExo0bMW7cOPj5+UFHR/7guZ+fH1auXFngejIyMmT7sybq2rUrDh48iOfPn6NmzZqy6aGhocjJyUHHjh0RFhZWhhUS8dQVaYFvvvkG3t7eePz4MXbt2iWb/ukYHaFQiA0bNsj+XygUytrkNT1XWFgYXF1dUadOHdSqVQuurq6Ijo6WqyF3PTdv3sSECRNQr149tGrVqkjruHv3LqZPn4569eqhdu3aGDFiBJKSkhSed1hYGNzc3GBhYYE6dergm2++QXBwsFyby5cvw9PTE5aWlqhRowY6duyII0eOFLKH5eWeitDT05ObHh8fj++++w52dnYwNzdHs2bNsGLFCtk3+ocPH6JBgwYAAH9/f1lfe3l54b///oNQKMS+fftk67t//z6EQqHCL7JPnz4dX3zxRZGeZ1paGn744Qc0btwY5ubmaNSoEX766SdIJBK5dkKhENOnT8fx48fRtm1biEQiNGvWDHv27CmwfzIyMrBs2TI0aNAAYrFYIeQAgLGxMX744Qe5aY0bN0bfvn0RHh6Ozp07QyQS4bfffgMAJCcnY8aMGWjQoAHMzc3RsmVLrF69Wu5oycOHDyEUCrF161aF7TVu3FjuB3NzTwWHh4fD29sb9evXR+3atTF8+HDEx8cX+Bxzubi4oEqVKgr9snv3bnTu3BlVq1bNczlV3g+PHj3CzJkz0aJFC9SsWROWlpbw9PTEjRs35Nrlnkbbs2cPVq9ejcaNG0MkEqFLly64cuWKXNsXL15gypQpaNiwIczNzWFnZwdPT09cu3ZN5edM5Q+DDmmF3KMOp06dUtpm7dq1aNeunez/165dCzc3N6XTAWDPnj3o27cvdHV1MXfuXMydOxdJSUno1asX/vnnH4VtjBo1CsnJyZg7dy4mTJhQpHWMHj0az549w9y5czF8+HAcPHgQs2fPlmuzY8cOeHh4ID4+HlOmTMHPP/8MR0dHHD16VNbm7Nmz6NatG168eAFvb2/8/PPPqFixIgYNGoTQ0FCV+vX9+/dITExEYmIinj17hmPHjmHFihVwcHCAvb29rN3Lly/RuXNnHD16FCNGjIC/vz+aN2+O+fPnw9fXFwBQvXp1LF26FADg6uoq6+tRo0ahYcOGEAqFiIyMlK0zMjISOjo6ePLkCR4+fCibHhUVha+//rrQzzMjIwOurq74888/4eHhgSVLlsDFxQWrV6+W/Wr3xy5evIhJkyahR48eWLhwISpVqoRx48bh1q1b+fZZdHQ0kpOT0a9fv0Kf2rt37x6GDx8OZ2dn+Pv7o0WLFpBIJHBzc8PmzZvRq1cvLF68GFZWVvjhhx/w/fffF2r9n/Lx8UFMTAxmz56NkSNH4u+//4aHhwfevXun0vIGBgZwc3PD7t27ZdPi4+MREREh+3HfT6n6frh8+TIiIyPh5uYGsVgMLy8vXL58GT169MhzjNjq1auxfft2jBs3Dj4+Prhz5w6GDBmC9+/fy9qMGDEC+/fvx6BBgxAQEIDx48cjJyen3J9CpPzx1BVphdq1a8PY2Fjh124/5unpifPnzyM8PFzudEyjRo3ynP7mzRvMmjULnp6eCAoKkk0fNWoUWrVqhQULFigEBhsbG/z555/FWscXX3yBdevWyR5LpVL88ccfWLZsGUxMTJCWlobZs2ejYcOGOHr0KCpXrizXNve/06dPR8uWLbF//37ZUYWxY8fCxcUF8+bNQ69evfLvVADh4eGwtraWm+bo6IgdO3ZAIBDIpi1atAgSiQSRkZEwNzeXPccaNWpg9erV8PLygpWVFXr16gVvb280bNhQ4ZRYq1atEBUVJXt87tw5dOzYEdHR0Th37hysrKyQmJiI27dv49tvvy308/z9998RFxeH06dPy44sAcCXX36JWbNmISoqCs7OzrLpN2/eRGRkpKxt79690ahRI2zZsgULFy5U2mc3b96UrfdjUqlU4cicsbGx3JGx+/fvY9u2bejRo4ds2rp16/Dff/9h5cqVGD58OABgzJgxGDZsGNasWYMxY8YovEaFcfDgQejr6wMA7OzsMGXKFGzbtg0jR45UafkBAwZg69atuHnzJuzs7LB7925UqVIF3bp1kwveQOHeD126dIG7u7vc8p6envj666/x559/YtasWXLz0tLSEBUVBQMDAwCAra0thg4dilOnTsHFxQWpqak4d+4cFi5ciClTpsiW+3QcH2kfHtEhrVGlShWkp6erbX1hYWFISUnBgAEDZEc1EhMTkZGRgfbt2+PcuXNy3xaBD0dj1L2O1q1bIzs7G0+ePJGtMy0tDTNnzpQLOQBk4ePq1auIi4vDgAEDkJycLNtucnIyOnfujAcPHuDRo0cF9sFXX32Fv/76C3/99Rd27tyJefPm4e7duxg6dCgyMjIAfPgA379/P1xcXKCrqyv3PDt16oScnBy5IzXKODs748aNG0hOTgbw4chN27Zt0bJlS1kAioyMhFQqlR3RKczz3LdvH5ycnFC9enW5Gtu3bw/gQ6j7WNu2beUCkbm5OWxtbfHgwYN8n0fuFYBGRkZy05OSkmBtbS337/z583JtateuLRdyAODo0aMwNTXFkCFDZNMEAgGmTp0KqVSKY8eO5VtPfkaNGiULOQAwaNAgmJiYFGqdbdu2Rc2aNWVHdXbv3g1XV1dZ4PhYYd4PH49ve/v2LZKSkmBiYgJra2vExMQorHvIkCFy22zTpg0AyF4vAwMD6Onp4ezZs7J9jD4PPKJDWiM9PR3Vq1dX2/ru3r0LAOjTp4/SNqmpqXLbrFu3brHXYWFhITc/d5xR7h/n3KNWH586Ulb7lClT5L69fuzVq1ewtLRUug4AqFatmiwIAB/GZNja2mLYsGEIDg7G+PHj8erVK6SkpGDLli1KLzXOHeydH2dnZ0ilUkRFRcHR0RH379+Hs7MzsrKysGPHDgAfwo+xsbFs3E5hnufdu3fx33//KT368WmNn74OwIfXoqAPydyA8+ktD4yNjfHXX38B+HC6LSAgQGFZKysrhWmPHj2CtbW1wmmw3BCmSmBV5tO+qFChAqysrPD48WOV16GjowMPDw/s3r0b/fv3R2xsLBYsWJBn28K8HzIzM/HLL79g165dCuOGTE1NFZYr6H2jr6+P+fPnY/78+bC1tUXz5s3RpUsXDBgwIM/XmrQHgw5phadPnyItLQ3169dX2zpzcnIAfDjlUatWrTzbGBsbyz3+9AqZoqxD2biOj09LAZA7daSs9p9++knplVE2NjZKl89P7nimqKgo2RgHAOjXr5/CZca5VHldHBwcULlyZURFRUEikaBSpUpwcHDA+/fvsXDhQrx8+VI2Pif3FFVhnmdOTg7atWuHGTNm5Nnu09enoNdBGTs7OwDA9evX4erqKpuup6cnC42JiYl5LlucK6xU2R9UWaYol4P3798fgYGB8Pb2Ro0aNZTeOqAw7wcfHx8EBwdj3LhxaNWqFYyNjaGjowNfX988n48qr9fkyZPh6uqKw4cP4/Tp01i6dCmWL1+Obdu24ZtvvinUc6byg0GHtMLOnTsBAB07dlTbOuvVqwfgwyDaj49qlPY6PpUbGq5fv65w9dGn261SpYratpsrKysLwIfxFsCH52ZsbIysrKwCt5Xfh3GFChXQvHlzREVF4d27d2jRogX09PTg6OgIfX19HDlyBNeuXYOHh4dsmcI8z3r16iE9PV3t/fEpJycnVK1aFXv27MHMmTMLPSD5U5aWlrhy5Qqys7Pl1nX79m3ZfACyK5xSU1PllpdIJEqvpLpz5w46dOgge5yVlYVHjx6hdevWharRwcEBX3zxBSIiIjBx4kSlz7kw74eQkBAMHDgQfn5+ctNTUlJQrVq1QtX3sbp162LixImYOHEinjx5gnbt2uHXX39l0NFiHKND5d6ZM2ewdOlSWFlZKb3Soyg6deoku2fMp5cfA6qdjlHHOj7VoUMHGBsbY/ny5Xj79q3cvNxvrw4ODrC2tsaqVasUPviKut1cueM3ck8f6erqolevXjh48GCeYydSU1MVxl2kpKTkuW5nZ2fExsbixIkTsoHB+vr6aNasGVauXIns7Gy5AcOFeZ4eHh74999/cfjwYYV2GRkZahvfZWhoiJkzZ+L27dtKjz4U5qiJi4sLXr16he3bt8stv2rVKggEAnTt2hXAh1Nm1atXV7jr9MaNG5GdnZ3nujdt2iS3X27fvh2pqano0qWLyvXlWrx4MebMmaMwxuxjhXk/6OrqKvTTnj178Pz580LXBnwY55M7rixXnTp1YGZmpnR/JO3AIzpUrpw8eRL37t1DVlYWXr58ifDwcISFhcHCwgLbt2/PcwBkURkZGWHFihUYPXo02rRpg/79+0MkEuHp06eIiIhA5cqVC7yvijrW8SljY2OIxWJMnjwZHTp0QP/+/VGtWjXcuHEDz58/x5YtW6Cjo4PVq1ejb9++aNWqFYYMGQJLS0vEx8fj4sWLePz4scJA2LzEx8fLjpa9e/cO165dw//+9z+Ymppi3LhxsnY//fQTIiMj0a1bNwwbNgz29vZ4/fo1rl+/jgMHDuDff/+FSCRClSpVYGtri5CQENjY2KBatWqwsrJC8+bNAQBff/01srOzZeNzcrVu3RoBAQEwNDSUO0VVmOc5ZcoUHDt2DMOGDcOAAQPg6OgIiUSCO3fuYN++fdi9e7fczSiLY9KkSbh79y7WrVuH8PBw9OrVCxYWFnjz5g1u3bqFvXv3onLlynmONfnU8OHDERwcjGnTpuHq1auwsbHB8ePHcezYMUyYMEFunM3IkSMREBCAiRMnokWLFrh8+TLOnDmT73bc3NzQt29fPHr0COvWrYOdnR0GDx5c6OfcpUuXAgNSYd4P3bt3x44dO2BkZAR7e3tcvXoVISEhCuPgVHXnzh306tULvXv3hp2dHfT19XHs2DHcunUr36voqPxj0KFyJfcwdsWKFVG1alXY29tDLBZjyJAhCle5qEPv3r1Rs2ZNLF++HL///jsyMjIgEonQvHlz2aW+pbGOTw0ZMgRmZmb49ddfsXz5cujq6sLa2hpjxoyRtfn6669x8uRJLFmyBP/73/+QlpYGMzMzNGrUSHZvm4Jcu3YN48ePB/AhVJiamsLV1RVz586VG2NRvXp1nDx5EkuXLsWhQ4fwv//9DyYmJrCxsYGPj4/cjeMCAwPh6+uLH374ARKJBIMGDZIFnRYtWqBixYoAIJuW+1xyp+XOL+zzNDQ0RGhoKFasWIGQkBBZ2Khbty68vLxga2urUp+oQiAQ4Ndff0XPnj2xceNGbN68GYmJiahUqRJsbGwwYcIEjBw5Uuk4lY8ZGBggNDQUCxcuxL59+5CcnAwrKyssXLgQkydPlms7a9YsJCUlISQkBH/99RfatGmD/fv3y+4L9Sk/Pz+EhobC398fEokELi4uWLp0qdyVWOqm6vvBz88Penp62LdvH7Zs2QIHBwfs3bsXP/74Y5G2W6dOHfTv3x/h4eHYs2cPBAKB7GjgsGHD1PX0SAMJUlJS+EMkRESfka1bt2LSpEk4fvy42o5iEWkqjtEhIiIircWgQ0RERFqLQYeIiIi0FsfoEBERkdbiER0iIiLSWgw6REREpLUYdIiIiEhrMegQERGR1mLQISIiIq3FoENERERa6/8A70kgfTPeCBoAAAAASUVORK5CYII=\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.012"
      ]
     },
     "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
}