{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[](https://mybinder.org/v2/gh/kasparvonbeelen/ghi_python/main?labpath=10_-_Hypothesis_Testing.ipynb)\n",
"\n",
"\n",
"# Lecture 10: Hypothesis Testing\n",
"## Data Science for Historians (with Python)\n",
"### A very gentle and practical introduction\n",
"### Created by Kaspar Beelen and Luke Blaxill\n",
"\n",
"### For the German Historical Institute, London\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We've covered how to describe, summarize, and compare variables. However, we lack a formal procedure to assess the importance of the differences we observe. For example, we established that men are, on average, one year younger than women. But how can we establish the value or 'significance' of this gap?\n",
"\n",
"In the notebook, we move from descriptive to inferential statistics and compute the extent to which means between subgroups in our data are statistically significant."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For more background on the concepts and terminology used in this notebook, please consult the lecture by Luke Blaxill."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We repeat some of the code from the previous notebook\n",
"- load the required libraries\n",
"- load the synthetic census data"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd\n",
"import numpy as np\n",
"import seaborn as sns\n",
"import random\n",
"from scipy import stats\n",
"from tqdm.notebook import tqdm\n",
"sns.set()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"df = pd.read_csv('data/democraphic_data/london_subsample.csv',index_col=0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To make our analysis more interesting and complex, we add a dimension: place. We study how age differences between men and women vary depending on the registration district. In Pandas, adding subgroups is convenient: simply pass a list with column names `['district', 'gender']` (instead of just one column as we've done previously). `.groupby()` will read this list from left to right, i.e. groups the data first by place and then further aggregates by `gender` within each district."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"district gender\n",
"Bethnal Green F 25.920963\n",
" M 25.613501\n",
" U 20.766603\n",
"Camberwell F 27.986109\n",
" M 26.708931\n",
" ... \n",
"Whitechapel M 26.254569\n",
" U 24.051316\n",
"Woolwich F 26.653639\n",
" M 26.311769\n",
" U 21.427165\n",
"Name: age, Length: 90, dtype: float64"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"by_reg_gen = df.groupby(['district','gender'])['age'].mean()\n",
"by_reg_gen"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Working with the output of this operation requires a bit more thought. The `.groupby()` arranges data slightly differently depending on whether you group on one column or more. This becomes apparent when printing the `type()` of the `.index` attribute."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pandas.core.indexes.base.Index"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"type(df.groupby('gender')['age'].median().index)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pandas.core.indexes.multi.MultiIndex"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"type(by_reg_gen.index)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`by_reg_gen` orders the data using a `MultiIndex`, which means that the index contains multiple levels (place and gender) via which we can access our data.\n",
"\n",
"Place sits the highest level of our grouped data frame. We can access the separate districts using `.loc[]`"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"gender\n",
"F 25.920963\n",
"M 25.613501\n",
"U 20.766603\n",
"Name: age, dtype: float64"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"by_reg_gen.loc['Bethnal Green']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Similarly, we can slice the data by place."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"district gender\n",
"Bethnal Green F 25.920963\n",
" M 25.613501\n",
" U 20.766603\n",
"Camberwell F 27.986109\n",
" M 26.708931\n",
" U 23.410050\n",
"Chelsea F 31.318046\n",
" M 30.430615\n",
" U 28.800000\n",
"Name: age, dtype: float64"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"by_reg_gen.loc['Bethnal Green':'Chelsea']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But our index has two levels, place **and** gender. We can obtain the means for women for all districts using the following syntax."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"district\n",
"Bethnal Green 25.920963\n",
"Camberwell 27.986109\n",
"Chelsea 31.318046\n",
"Fulham 27.920645\n",
"Greenwich 27.593689\n",
"Hackney 28.796820\n",
"Hampstead 29.742615\n",
"Holborn 27.719117\n",
"Islington 29.047747\n",
"Kensington 30.823465\n",
"Lambeth 29.019583\n",
"Lewisham 28.879865\n",
"London City 31.323944\n",
"Marylebone 30.465945\n",
"Mile End Old Town 26.298023\n",
"Paddington 30.450756\n",
"Pancras 29.136537\n",
"Poplar 26.384635\n",
"Shoreditch 27.056387\n",
"Southwark 27.144232\n",
"St George Hanover Square 30.821499\n",
"St George In The East 24.545266\n",
"St Giles 30.129131\n",
"St Olave Southwark 26.486323\n",
"Stepney 25.985841\n",
"Strand 29.288732\n",
"Wandsworth 28.028524\n",
"Westminster 28.990447\n",
"Whitechapel 24.649451\n",
"Woolwich 26.653639\n",
"Name: age, dtype: float64"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"by_reg_gen.loc[:,'F']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice the comma in `.loc[]` (`by_reg_gen.loc[:**,**'F']`) \n",
"- the part before the comma indicates the items we want to access from the first level (place). In this case, we entered a colon, meaning from the first till the last rows. \n",
"- the part after the comma indicates the items we want to select from the second level. In this case, we want to retrieve all elements with value 'F' for the second level of the MultiIndex..\n",
"\n",
"Computing the age differences by place is then fairly straightforward."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"district\n",
"Bethnal Green 0.307461\n",
"Camberwell 1.277178\n",
"Chelsea 0.887431\n",
"Fulham 1.295697\n",
"Greenwich 0.971892\n",
"Hackney 1.357030\n",
"Hampstead 0.645606\n",
"Holborn 0.341140\n",
"Islington 1.430837\n",
"Kensington 1.410332\n",
"Lambeth 1.411260\n",
"Lewisham 1.420661\n",
"London City 3.262573\n",
"Marylebone 0.921905\n",
"Mile End Old Town 0.599087\n",
"Paddington 1.664367\n",
"Pancras 1.031302\n",
"Poplar -0.218534\n",
"Shoreditch 0.130194\n",
"Southwark 0.101532\n",
"St George Hanover Square 1.412338\n",
"St George In The East -0.011402\n",
"St Giles -0.879827\n",
"St Olave Southwark 0.088463\n",
"Stepney -0.334267\n",
"Strand -1.759558\n",
"Wandsworth 1.004169\n",
"Westminster -0.329065\n",
"Whitechapel -1.605118\n",
"Woolwich 0.341871\n",
"Name: age, dtype: float64"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f_m_diff = by_reg_gen.loc[:,'F'] - by_reg_gen.loc[:,'M']\n",
"f_m_diff"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Of course, you can be more restrictive and slice the data by place and obtain only means for women."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"district gender\n",
"Bethnal Green F 25.920963\n",
"Camberwell F 27.986109\n",
"Chelsea F 31.318046\n",
"Name: age, dtype: float64"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"by_reg_gen.loc['Bethnal Green':'Chelsea','F']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Hypothesis testing"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At this stage, we can compute and compare the distribution of variables, calculate their means or other relevant statistics. But a question then immediately appears: are the differences we observe \"significant\", and what do we mean with \"significance\" anyway. In this section, we have a closer look at hypothesis testing from a data-driven perspective.\n",
"\n",
"Traditionally, statistical methods, such as the Student's t-test arose in times of limited computing power and relied on equations and assumptions related to the distribution of the data. Explaining this requires many detours and implies a steep learning curve for the statistically uninitiated. \n",
"\n",
"In this lecture, we explore a data-driven, and hopefully intuitive procedure, for assessing the significance of observed statistics (such as the mean). First, we discuss a procedure called **bootstrapping** and then explain how to use **permutation** for significance testing. \n",
"\n",
"The question we first address regards the relation between `gender` and `age`, and we went to assess if variations we observe over the districts is significant. We focus on Whitechapel and Poplar as we noticed in the previous section that both these districts deviate from the general pattern and contain a slightly younger female population.\n",
"\n",
"The age gap is bigger in Westminster than in Poplar, but is the differences of the mean significant?"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"((18772, 4), (40838, 4))"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_whitechapel = df[df.district=='Whitechapel']\n",
"df_poplar = df[df.district=='Poplar']\n",
"df_whitechapel.shape, df_poplar.shape"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
" mean std\n",
"gender \n",
"F 26.384635 19.047276\n",
"M 26.603169 18.596804\n",
"U 22.963240 18.230697"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_poplar.groupby('gender')['age'].agg([np.mean,np.std])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bootstrap: Compute the Distribution of a Sample Statistic"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The means give us an idea of the expected age in different registration districts. However, our statistic is derived from partial data, a subsection of Londoners at a particular point in time and therefore an unreliable measure: difference samples may produce difference means, and we'd like to know how much variation we'd expect if we were to take more samples.\n",
"\n",
"Of course, we can not do this, we only have this particular data set. \n",
"But one statistic remains a weak empirical basis. To estimate the variation and produce confidence intervals for the mean we will follow a statistical procedure called bootstrapping, which is simple, intuitive and requires fewer assumptions about the distribution of the data. Also, it avoids equations and statistical theory and therefore easier to grasp. \n",
"\n",
"Will use the bootstrap method to compute the variation we can expect around the mean for women in Whitechapel.\n"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"regdist = 'Whitechapel' # change to Poplar or any other registration district\n",
"df_sub_F = df[(df.gender=='F') & (df.district==regdist)]\n",
"df_sub_M = df[(df.gender=='M') & (df.district==regdist)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Before we start, observe the different distribution of `age` variable for each value in `gender` (blue='F', orange='M') Actually the look pretty similar, no? "
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAD7CAYAAABNEGKbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU5d3//9eZLclkISSZJBD2RcIWtshmBbFiZDMpqHVp49KC+ru1lO8trRsu3CpVsWqrUm+qtVa4BZeGYlmLohawQhAIEJGwk0B2ss52Zs7vjyEDgSRMlslk+TwfDx8y55w5856TTD5zXeec61I0TdMQQgghmkkX6ABCCCE6BikoQgghWoQUFCGEEC1CCooQQogWIQVFCCFEi5CCIoQQokVIQRFCCNEiDIEO0JpKS6twuxu+7SY6Oozi4spWStR4kq95JF/zSL7maW/5dDqFrl1DfX5+pyoobrd2xYJSs11bJvmaR/I1j+Rrno6cz69dXmvXrmX69OlMnTqVFStWXLY+OzubOXPmkJKSwhNPPIGqqgBkZmYyZ84cUlNTufvuu8nNzQVg586djBs3jtTUVFJTU3nsscf8GV8IIUQj+K2g5Ofn8+qrr7Jy5UrWrFnDqlWryMnJqbXNwoULWbRoERs3bkTTNFavXu1d/vzzz7NmzRpmzZrFc889B0BWVhb33Xcfa9asYc2aNSxZssRf8YUQQjSS3wrK9u3bGT9+PJGRkZjNZlJSUtiwYYN3fW5uLjabjZEjRwIwe/ZsNmzYgMPhYP78+SQmJgIwaNAgzpw5A3gKyrZt20hLS+OBBx7wLhdCCBF4fisoBQUFWCwW7+PY2Fjy8/PrXW+xWMjPz8dkMpGamgqA2+3mjTfe4IYbbgAgPDyc9PR0MjIymDx5MgsWLPBXfCGEEI3kt5PydQ1irCiKz+sdDgePPvooqqpy//33A7B48WLv+jvuuINXXnmFiooKwsPDfcoUHR3m03YWi2/7CxTJ1zySr3kkX/N05Hx+KyhxcXHs2rXL+7igoIDY2Nha64uKiryPCwsLveurqqp48MEHiYyMZNmyZRiNRtxuN2+//Tbz5s1Dr9dfeAMG399CcXHlFa9gsFjCKSys8HmfrU3yNY/kax7J1zztLZ9Op/j8RRz82OU1ceJEduzYQUlJCVarlU2bNjFp0iTv+oSEBIKCgsjMzAQgIyPDu37hwoX07t2b119/HZPJ5Amq07F582Y2btzo3X7EiBGEhIT46y2IFpZfWs2jf9rBo2/voKC0OtBxhBAtzG8FJS4ujgULFpCenk5aWhozZ84kKSmJuXPnkpWVBcDSpUtZsmQJ06ZNw2q1kp6ezsGDB9myZQu7d+8mLS2N1NRU5s6dC8CLL77I+++/z4wZM/jkk0+8V3+J9uGjL45QVu2gotrJ2/84UGe3pxCi/VI604yN0uXlf/XlO1dp55E3t5Mytifx0Wb+su57fn1rEkn9Y9pEvrZC8jWP5GueNtvlJcTFDh4vwa1pjBsSx4Sh8USYjXy1Vy77FqIjkYIiWsWhk+cIDTbQIzYMg17H+KHx7M0porzaEehoQogWIgVFtIqc3DIGJHRBd/7S8InD4nG5NfYcLrrCM4UQ7YUUFOF3TtVNfomVnnEX+mJ7xoYRHREsBUWIDkQKivC7syXVuDWNhJgLBUVRFEYOjOHA8RLsTlcA0wkhWooUFOF3uUWe+RUSLLXnVRg5MAan6ib7eGkgYgkhWpgUFOF3BSVWAOK61r4JdVDPSEKC9Hx3uDAQsYQQLUwKivC7onIbXcJMGA36WssNeh3D+0Wz70gx7s5zO5QQHZYUFOF3xWU2YiKC61w3on8MZVUOTpxtuzd7CSF8IwVF+F1xmY3oLnUXlOH9o1EU2JsjV3sJ0d5JQRF+5dY0SipsRNfTQgkLMTIgoQt7c4pbOZkQoqVJQRF+VVbpQHVp9bZQAEYMiOFEfgWlFfZWTCaEaGlSUIRfFZfbAIi5QkEB2HtEur2EaM+koAi/KjlfUKLq6fIC6B5tJqZLMHvlrnkh2jUpKMKvyio9gz92CTXVu42iKIwYEMPBE6Vy17wQ7ZgUFOFX5dUO9DqF0BBjg9vV3DW//2jTT8673G5UlxsATdNwVxSi2auavD8hROP4bU55IQDKqhyEm43eUYbrk9grkohQE98cyGfMoNhGvUZFtYNPvjzCjgP5uFwak/sqpPIvlHOnQdFjHDKFoAm3o+jk110If5JPmPCr8ioHEQ10d9XQ63SMHRzL1u9yqbY5MQc33KKpUVRm5aWV33Gu0s41w7sRqVSSfOwdqnQQmnw7pqqzOA/8C81aRvCPH2zu2xFCNEC6vIRf+VpQACYMjUd1aew65NvYXnanizc+yaLKpvLbu0aTnjKIqa6thBndvFWRwt9OJBB87T2Yxt6GenQnzqyNzXkrQogrkIIi/KqsykEXs28FpU98OHFRZv6978pTA2uaxl83fM+pgkrmzRpC/+5dcJ3Yg+v0fkLG3cL4iSPJPFTIviNFmEZMw9BnNPadn+AsK2juWxJC1EMKivAbTdOoqL7QQtE0Dde5PDSXWuf2iqJw/agEcnLLOJJb1uC+/7XrNN8cyCft2r6MGBCDpmnY93yGEh6Dccj1pIztRUyXYNb8+xgAQRPvAnSUfP5Bi75HIcQFUlCE31TbVVSX5i0o9n+/T/Xqx6leuwRNrXsu+WtHdMMcZGDDtyfr3e/3J0pZ9XkOowbGMGNiHwBcZ3/AXXAEU9JNKDo9Br2OmRP7cOxMBdknStGFRWMaMY2qg9twFR1v6bcqhEAKivCj8ipP0YgINeEqOoEz+wt0cQNwFxzBsW99nc8JNhm4fkwCuw8VcjSv/LL1ReesLFuzn7ioEH45c4j36jFn9lYwmTEOuta77YShcYQGG/hqbx4ApuE3ophCcOxZ18LvVAgBUlCEH11cUNQj/wFFj/mmBRh6j8KRtaneVsq0cb2JCDPx1w3f17rRsdLq5NWP9uJyaTw0ezghQZ6LFDWHFfVYJsb+41AMQd7tjQY9E4bFs/uHQiqqHShBoUSMSUE9thN3uZxLEaKlSUERflNp9ZwrCQ8xop74Dn3CYJSgUIxDfwz2KtSTe+p8XkiQgXtuSuR0QSXLMvZzrtLOkdwylnyQSeE5Kw/PGU636AvTCavHdoHLgfGqay7b1+QR3VFdGt8cyAegy9UzQdHj2CdXfAnR0qSgCL+psjkBCFUcuM+dQd89EQB99yEooV1x/rCt3ueOGBDDz1IGsf9oCf/vjW08/7dMKqqd/L/bRjKoV9da2zqP/AclIhZdbP/L9pNgCaNnbBg7D3laJIbwrhj6j8V5eDuaU0Y3FqIlyY2Nwm9qCkpI1SlcgN7SDwBFp8PQbyzOA1vQHFYUU0idz58yKoHEXpHsySkiNNjI1Ymx3m6uGprDiivve4zDbkCp52785EEW/v71MUor7Fgs4RgHT0E9vB3nkW8wJU5uuTcsRCcnLRThN1VWFb1OQV96CgC9pY93naH3SHCrqLkHG9xHt+hQpo3rzaQR3S8rJgDq6f3gVjH0HlXvPpITPUO57P7Bc8OkPm4Auq4JnhP5QogWIwVF+E2VzUlosAF3+VmU0K4oJrN3nT5+IJhCUE/UfR7FV+qJPRAUij5uQL3bdIsOJSEmlF3fe7q9FEXBOPg63IXHcJWcatbrCyEukIIi/KbK6iQ0xIi7LB9dRFytdYrOgKFnEq5Te9E0d5P2r7nduE7tw9AzCUWnb3DbMYMs/HDqHGWVnvMmhv7jQNGhHt7RpNcWQlxOCorwmyqbSmiwEa0sH12XuMvWG3qNQLOW4y483qT9uwqOoNkqPN1nVzBiQAwa8N35k/O6kAj0PYbhzPmmyQVNCFGbXwvK2rVrmT59OlOnTmXFihWXrc/OzmbOnDmkpKTwxBNPoKqey0wzMzOZM2cOqamp3H333eTm5gJQXl7OvHnzmDZtGnfddReFhb4NIigCo8rmJDLIhWarqLOg6HsMAxTUU1lN2r/rxHeg6DH0HH7FbXvHhxNhNrIzO9+7zDhwAlpVCa6zh5v0+kKI2vxWUPLz83n11VdZuXIla9asYdWqVeTk5NTaZuHChSxatIiNGzeiaRqrV6/2Ln/++edZs2YNs2bN4rnnngPgtddeIzk5mfXr13Prrbfy/PPP+yu+aAFVVpU4QwUASh0FRRcSgc7SF/XU3ibtXz25B323q2qdm6mPTlEY3i+a3d8X4HJ7WiSG3qPBYELNkW4vIVqC3wrK9u3bGT9+PJGRkZjNZlJSUtiwYYN3fW5uLjabjZEjPd0Vs2fPZsOGDTgcDubPn09ioueehUGDBnHmjGf02a1btzJr1iwAZs6cyVdffYXT6fTXWxDNVGVzEqPzFJRLz6HUMPRKwl1wDLetolH7dpcX4C7N86m7q0bSgBgqrU6O5HqGdFGMQRh6JqEe/066vYRoAX4rKAUFBVgsFu/j2NhY8vPz611vsVjIz8/HZDKRmpoKgNvt5o033uCGG2647DkGg4GwsDBKSkr89RZEM6guNzaHiy5KJQC6sOg6tzP0TAI0XI3s9qq5Oqyhy4UvNbRPFHqdQtZF0wwb+oxGs5bhLjzWqNcXQlzObzc2app22bKLbzy70nqHw8Gjjz6Kqqrcf//99b6OTud7TYyODvNpO4sl3Od9BkJ7yHeuwnM1VVeDDcUUQmxC3dP6ajHDObEpAkPh91gmpvj8GnlnsjDG9CCu/+V3xzdkSN9oDhwv5YFbPMfQFXYNJ7b+GVPBfqKG+t7a8af28PNtyyRf8zQnn98KSlxcHLt27fI+LigoIDY2ttb6oqIi7+PCwkLv+qqqKh588EEiIyNZtmwZRqNnOtjY2FiKioqIj49HVVUqKyuJjIz0OVNxcSVu9+WF7GIWSziFhY3rfmlN7SXfmeIqAPS2MpSQLg1m1nUfSlXOdxTkl6H48AVBs1dhO3EQU1JKo49F8uBY/vLZQQ4dKSQqItiTsdsgyg/+B9ew1Ebtyx/ay8+3rZJ8zXNpPp1O8fmLOPixy2vixIns2LGDkpISrFYrmzZtYtKkSd71CQkJBAUFkZmZCUBGRoZ3/cKFC+nduzevv/46JtOF2f4mT55MRkYGAOvWrSM5OdlbbETbUnV+YMggtRwltGuD2xp6jUCzVeAu8q3bST21DzRXo7q7aiQP9pzL2XdJt5f7XB7usvz6niaE8IHfCkpcXBwLFiwgPT2dtLQ0Zs6cSVJSEnPnziUry9NfvnTpUpYsWcK0adOwWq2kp6dz8OBBtmzZwu7du0lLSyM1NZW5c+cCMH/+fPbs2cOMGTNYuXIlTz31lL/ii2aqPD+Ol9FRceWC0mMYKArqyX0+7Vs9/h1KSESdg0FeSc+4cGK6BLMv56KC0sNz2bGae6DR+xNCXODXwSFnzZrlvSqrxvLly73/TkxM5OOPP661fsiQIRw6dKjO/UVGRvKnP/2p5YOKFldldaKgobOVoTM33C2pBIehi+2PemofQck/aXBbzeVEPbXPM/dJI86feV9LURjRP4avs/JwOF2YjHqULnEoYdG4Th+AIdc3ep9CCA+5U174RbVNxazYUTTXFVsoAIaew3EXHsddfa7B7Vx52eC0YejT+O6uGkkDonE43Xx/0vNaiqJg6DEUNfcgmtt1hWcLIeojBUX4RbVdJVJXDeBbQel7NaChHvm2we2cR3aCMRh99yFNzpbYKxKTUcfeIxcuCtEnDAOnVS4fFqIZpKAIv7DaVaJNnkuHr9TlBaDv2h1ddG+cOd/Uu42m2lGP7cTY72oUg6ne7a7EaNAztE8U+3KKvJevGxKGAArqaTmPIkRTSUERfmFzqHQ1euaMV0IifHqOccB43IVHcZedrXO9eny3p7tr4OVT/TbWiAExFJfbyS30XN6sBIehs/TBJSfmhWgyKSjCL6x2F10M5wtKsG/XsRsGjPfM937wizrXO7//CiUsGn23q5qdb3g/z537F3d7GRKG4so/guawNnv/QnRGUlCEX1gdKhF6O+j0YKx7it9L6UK7Yuh3Nc7vv7rsj7qr8DiuvGxMw25AUZr/a9s1PIje8eHsPXLh8mF9whDQXLjO/tDs/QvRGUlBEX5hs7sI09lRgsPrneu9Lqakm8BpxbHnn95lmqZh3/UJmEIwJl7XYhlH9I/mSG4ZFdWelpQ+bgDoDVecllgIUTcpKMIvrA6VUMXuc3dXDb2lD4aBE3Hs24Cr4AgA6uHtuE5lETTmJygm31o7vhgxIAZNg/1HPQOMKgYT+riBnkuThRCNJgVF+IXNrhKCDSW48QPNBY2/HSW0K9XrXsG6ZRm2L99FH38VxqE/btGMvePD6RJqqn35cPfBuItPNno4fSGEFBThJ1a7i2CtutEtFPBMvGWesRB9bD9cuQcxDBxPSMr8K84b3+jXURSG948m62gJquv8pFsJnvtbXHnft+hrCdEZ+HXoFdE5aZqGzeEiyG1tUgsFQBcRi3n6Iy2c7HIj+sfw731nyDldRmLvrugsfcEYjCv3IMZ+V/v99YXoSKSFIlqcQ3WjaS6MrqZ1ebWmIX26YtAr3m4vRadH320Qap6cmBeisaSgiBZns6uYFQcKWpO6vFpTSJCBQb26sufi0Ye7D0Ery8ddWdzAM4UQl5KCIlqc1eEiTLEBtPkWCnguH84vqeZsiWfsMb33PIpc7SVEY0hBES3OalcJ1XnG8WoPBWXkgBgA9hz2dHvpohJQgsNRc6WgCNEYUlBEi7PZVcKUmoLStru8AGIiQ+hhCWNPzvnzKIoOfffBuPIOegePFEJcmRQU0eKsDhehuvbT5QUwamAMh0+fo9LqmWlSnzAEraoUrZ6BKoUQl5OCIlqctZ21UABGDvTcNb/v/NVehu6DAVDlPIoQPpOCIlqczeEiRLGj6Y3NmrekNfWODycyzOQ9j6JExHqmBZZxvYTwmRQU0eJsDpUQxYFiMgc6is90isLIATFkHSvBqbo90wInDEU9fQDN5Qx0PCHaBSkoosVZ7S5C9U50QaGBjtIoIwfGYHe4OHSyFABD32RwWnHJLI5C+EQKimhxVodKqM4JQe2nhQIwuHdXTEYd352/2kufMASCQnEebXieeyGEhxQU0eJsdhWzrn11eYFnrvlhfaPZc9gz17yiN2DoPRr1+HdoqiPQ8YRo86SgiBZntbs851DaWQsFPDc5llbYOZlfCYDxqongtKIe+U+AkwnR9klBES3O5lAJwYFial/nUACSBkSjgPcmR323RHRdu+PY/y+5yVGIK5CCIlqc1aZiwt4uWygRZhP9e3S5cPmwomAcnoK7+ATq8d0BTidE2yYFRbQ4l9OKDq3dnUOpMWpADCfyKygp99ztb7zqR+giu2P/5kM0hzXA6YRou6SgiBan2D2j9ra3q7xqjBzoGSxyb86FOVKCJt2DVlmM7Yv/RXOpgYwnRJslBUW0OJ3q+RbfXlso8VFm4rqGeC8fBjDEX0XQxDtRT3yH9Z8v4S4vDGBCIdomKSiiRakuNwb3+XG82tmNjTUURWHkwBi+P1GK1X6hNWIaegPBU+bhKj5F1cdPnj9R7w5gUiHaFikookVZ7SohNQNDttMuL4DRV1lQXZr35HwN48CJhN76HPr4gdi3f4Dti+VobleAUgrRtvi1oKxdu5bp06czdepUVqxYcdn67Oxs5syZQ0pKCk888QSqWrtv+vXXX+ePf/yj9/HOnTsZN24cqamppKam8thjj/kzvmiCaptKiOIZ+6q9dnkBDEjoQkyXYLYfuHz4el1YNCHT/htT8mzUnB04dn0agIRCtD1+Kyj5+fm8+uqrrFy5kjVr1rBq1SpycnJqbbNw4UIWLVrExo0b0TSN1atXA1BRUcHjjz/Ou+++W2v7rKws7rvvPtasWcOaNWtYsmSJv+KLJqq2OTErnrvK23NBURSF8UPjOXi8hHOV9jrXB42+GWPiZBx7/omr4EgAUgrRtvitoGzfvp3x48cTGRmJ2WwmJSWFDRs2eNfn5uZis9kYOXIkALNnz/au37JlC3369OHee++ttc+srCy2bdtGWloaDzzwAGfOnPFXfNFE1TaVkPPT/9KOCwrAhKFxaBr852B+vdsETbgDJSQC+7cfy42PotPzW0EpKCjAYrF4H8fGxpKfn1/veovF4l2flpbGvHnz0Ov1tfYZHh5Oeno6GRkZTJ48mQULFvgrvmgizzkUJ25DMIqufZ+i6xYdSt9uEXy1N6/eYqEYgzGNnIkrLxt34dFWTihE22Lw147r+gAqiuLz+rosXrzY++877riDV155hYqKCsLDfZtmNjrat9kDLZa2PW1tW86XfbqMEMWBLji0zeZsTK7Uyf157cPvyDtnY+RVsXVu475mGid2fYr+2DYsQ0e2ar5AkHzN05Hz+a2gxMXFsWvXLu/jgoICYmNja60vKrpwBU1hYWGt9Zdyu928/fbbl7VcDAbf30JxcSVud8PdEhZLOIWFFT7vs7W19XzVNhWzYkczhLTJnI09foN7RBBuNvLxv34goWtIvdsZ+o+l4sC/0Ub/FMUY1Gr5Wpvka572lk+nU3z+Ig5+7PKaOHEiO3bsoKSkBKvVyqZNm5g0aZJ3fUJCAkFBQWRmZgKQkZFRa/1lQXU6Nm/ezMaNG73bjxgxgpCQ+j/kovXVdHm150uGL2Y06JkyKoG9R4o5cbb+PwSG/uNBdaCezmrFdEK0LX4rKHFxcSxYsID09HTS0tKYOXMmSUlJzJ07l6wsz4du6dKlLFmyhGnTpmG1WklPT29wny+++CLvv/8+M2bM4JNPPuG5557zV3zRRJ7Lhh3og9vnTY11ufHqnpiDDPz96/rPkei7DYKgUNRjma2YTIi2xaf+oocffpg77riDiRMnNmrns2bNYtasWbWWLV++3PvvxMREPv744wZf92IDBw7kww8/bFQG0bqq7U5CdY52e5d8XczBRqaN78UnXx5l/7FihvWNvmwbRafH0HsU6vHdaG53u78gQYim8Om3/sYbb+Stt94iJSWFd955h3Pnzvk7l2inrOdbKO35HpS63Hh1T+KjzPx1/SFsjroHhzT0GAaOatxFx1s3nBBthE8FZdasWXzwwQe89dZbFBcXc+utt7Jw4UL27dvn73yinbFaHQQpThRTxzq3ZTTouXd6IiXlNt5b/32dVynquw8GQM072NrxhGgTfG6Xu91uTpw4wfHjx1FVlejoaJ555hlefvllf+YT7YzTWgW073G86jOwRyRzruvPt9kFrPvmxGXrdeYu6KJ64MqVgiI6J5/Oobz66qt8+umn9OzZkzvvvJPXX38do9FIdXU1U6ZMYeHChf7OKdoJ1/m5UDpal1eNaeN6cTK/gk++PEpcVzPJibUvddd3H4Iz+ws01YFiMAUopRCB4VNBKSkpYfny5SQmJtZabjabeeWVV/wSTLRP7prJtTpYl1cNRVG4b/pgisttLP/sIF0jgujfvYt3vaH7YJz7N+EqOo4h/qoAJhWi9fnU5eVyuS4rJjVXYP3oRz9q+VSi/appoRg7ZkEBMBn1PDwnicgwE3/8eJ93qmAAXVx/ANz5OfU9XYgOq8EWytNPP01+fj6ZmZmUlJR4l6uqytGjMm6RuJyiWkFPhzspf6kIs4lf3zqCxe/t4r3137PgthEoioIuJAIlIhZXvow+LDqfBgvKLbfcwuHDhzl06BApKSne5Xq9nlGjRvk9nGhfNE1DUW0Q1HHPoVysW3Qot1zXnxWbf2D7/rNcM7wbAPrY/rhyD3iOxxXGpxOiI2mwoAwfPpzhw4dzzTXXEBcX11qZRDtld7oIxjMXSkc9h3KpKaMT2L7/DBlfH2Xs4DiMBh36+IGoOTvQKopQIixX3okQHUSDBWX+/Pm8/vrr/PKXv6xz/dq1a/0SSrRPVruLYO9sjZ2joOgUhZ9M6sfvV+1l2/4zXDcyAX2s5zyKqyAHnRQU0Yk0WFDmzp0LwKJFi1oljGjfbA7PXfJuRd+pLpkd2ieKXrFhfJ6Zy+QR3dFF9QBDEK78HIwDJgQ6nhCtpsGrvIYNGwbA2LFj6datG2PHjuXcuXN8++23DB48uFUCivajpoWiGYIDHaVVKYrCdaMTOF1YyZHcchSdHn1sP1xypZfoZHy6bPipp55i+fLlHDlyhMWLF5Obm8sTTzzh72yinbGeb6FoHfiS4fqMHxJHSJCeL/fkAqCPG4C7+BSa8/L56IXoqHwqKPv37+eZZ55h8+bN/OQnP2HJkiXk5ub6O5toZ2x21XMOpRMWlGCTgdFXWdh9uAin6vacR9HcuAqPBTqaEK3Gp4KiaRo6nY5t27Yxfvx4AKxWq1+DifbHand5RhoO6nwFBeDqxFisdpWDx0u8Nzi6CuR+FNF5+FRQevXqxdy5czl9+jRjx47lv//7vxk0aJC/s4l2xubwtFD0HXBgSF8M6RNFSJCBXYcK0AWHo3SJkzvmRafi01heS5YsYfPmzYwZMwaj0UhycjJpaWn+zibaGavD00LRdaDJtRrDoNcxckAMew4X4XZrnhscT++XGxxFp+FTC8VsNpOcnEx5eTkHDhwgKSlJhl4Rl7HZVYJ1TvTBnbOFAjBiQDRVNpVjZ8rRxw1As5ajVRQGOpYQrcKnFsrLL7/MBx98QHT0halPFUVhy5Ytfgsm2h+rXSWIjje5VmMM6ROFokDW0WL6DL5wHkUXEXuFZwrR/vlUUNavX8+mTZtk+BXRIJetGp3See6Sr0tYiJF+3SPIOlpC6jWj5AZH0an41OXVrVs3KSbiilzeuVA6b5cXwPC+0Rw/U06FzYXe0ldGHhadhk8FZcKECbz00ktkZmZy4MAB739CXMzt8FxK3plbKADD+0ejAQePlVy4wVGVGxxFx+dTl9enn34KwIYNG7zL5ByKuIyj40+u5YveceGEBhs4cKyE5OH9QXPhKjyOoZtcai86Np8Kyueff+7vHKIjcFpBJy0UnU5hcJ8oDhwvQbl+OACu/CNSUESH51OXV1VVFYsXL+buu+/m3LlzPPXUU1RVVfk7m2hndOr5qXA7+TkUgKF9unKu0sHZagNKRByus4cCHUkIv/OpoDz33HOEh4dTXFxMUFAQlZWVPPXUU77SH9oAACAASURBVP7OJtoZvctTUDp7CwU8Q9qD5zyKIWEwrjOH0NxqgFMJ4V8+FZTs7GwWLFiAwWAgJCSEpUuXkp2d7e9soh1xqm5Mmme2xs4w/e+VxESGENc1hAPHS9AnDAWnDVeBDBQpOjafCopOV3szl8t12TLRudVMrqWhg040uVZDhvSN4tDJcxA/CFBwnd4f6EhC+JVPVeHqq6/m5Zdfxmaz8fXXX/PQQw8xbtw4f2cT7YjV4Zlcy20MlnGrzhvaJwq708XRIhc6Sx9cuQcDHUkIv/KpoDzyyCOYzWbCw8N57bXXSExM5De/+Y2/s4l2xGrrvJNr1SexV1d0isKB4yUYEobiKjiC5pBpH0THdcWCsnnzZn7+85/z5z//mdOnTxMeHs7o0aMJCgpqjXyinagZul7On1xgDjbQt3s4B46Vou+VBJob9eTeQMcSwm8aLCgZGRm8/PLL/OxnP+Ojjz7igw8+IC0tjeeff55NmzZdcedr165l+vTpTJ06lRUrVly2Pjs7mzlz5pCSksITTzyBqta+Cub111/nj3/8o/dxeXk58+bNY9q0adx1110UFsoorm1F9fnZGnWddHKt+gztE8Xxs+VYI3qjmCNRj34b6EhC+E2DBeVvf/sb7733HjNnzmTAgAH079+ftLQ03n77bf7yl780uOP8/HxeffVVVq5cyZo1a1i1ahU5ObUnG1q4cCGLFi1i48aNaJrG6tWrAaioqODxxx/n3XffrbX9a6+9RnJyMuvXr+fWW2/l+eefb8p7Fn5gs7sIVhyddnKt+gztG4WmwaGTZRj6XY16ap90e4kOq8GC4nQ66d69+2XL+/bti93e8NhE27dvZ/z48URGRmI2m0lJSak1dEtubi42m42RI0cCMHv2bO/6LVu20KdPH+69995a+9y6dSuzZs0CYObMmXz11Vc4nU4f3qbwt2q7SojixBDSOSfXqk/fbhEEm/Se8yj9xoJLRT22K9CxhPCLBguKXq+vd52maQ3uuKCgAIvF4n0cGxtLfn5+vestFot3fVpaGvPmzbvs9S9+jsFgICwsjJKSkgZziNbhOYfiwGgOC3SUNsWg15HYqysHzg8UqYvsjuPAlit+foRoj3way6sp6vrAXHw56ZXW+6ox98NER/v2x85iCW90jtbUJvMpCsHnWyjRbTHfRVr7+I0b3o09f8/CpdMTNX4mRRv+lwh7HsE9E9tEvsaSfM3TkfM1WFAOHTrE6NGjL1uuaRoOh6PBHcfFxbFr14WmfUFBAbGxsbXWFxUVeR8XFhbWWl+X2NhYioqKiI+PR1VVKisriYyMbPA5FysursTtbvibocUSTmFhhc/7bG1tNd+50nL0ioY+OLRN5qsRiOPXK8ZzXunr3ae4btgYCAolf8tKQqY/ctmXqLb6860h+ZqnveXT6RSfv4jDFQrK5s2bmxxs4sSJ/PGPf6SkpISQkBA2bdrE//zP/3jXJyQkEBQURGZmJmPGjCEjI4NJkyY1uM/JkyeTkZHBAw88wLp160hOTsZoNDY5o2g5bptnsFCdnJS/THyUmeiIILKOFDNlVAJBo1Ox71iJ6+QeDL1HBTqeEC2mwYKSkJDQ5B3HxcWxYMEC0tPTcTqd3HLLLSQlJTF37lx+9atfMXz4cJYuXcqTTz5JVVUVQ4YMIT09vcF9zp8/n0cffZQZM2YQHh7O0qVLm5xPtKya2RqloFxOURRGDbSwdU8eVrtK8NDrcWZ/ge3ffyM0biBKsJx3Eh2DonWis4PS5eU/f/nrZ9xi/5j4nz5BVZeBgY5Tr0Adv8Onz7Hkg93MnTWECUPjcRUcpfofz6PvMZyQlF+hKLqA5vOV5Gue9pavsV1eMsKjaBGa0zN0vbRQ6tY/oQtdw4PY9X0BAPrYfgSNvx3XyT04vlsb4HRCtAwpKKJF6FTPzXpSUOqmUxTGXGUh62gJVrtnRAjj0BswDJyIY9ffUY9/F+CEQjSfFBTRIpSaFkqwFJT6jBsSh+py8222534rRVEIvvYedJa+WL94G9e5vAAnFKJ5pKCIZtM0zTtbo7RQ6tevewQJllC27rlQOBSDiZCpD6EYTNg+fxvN7QpgQiGaRwqKaDaH002w4kBDQTEFBzpOm6UoCteNTODE2QqOny33LteFRRM04U7cRSeo2Pt5ABMK0TxSUESzWc8PXe/SB3mvVhJ1mzA0jiCTng3/OVlruaH/OPTxV1H65f+hqQ3fNCxEWyWfftFsVrtnci23QVonV2IONnL96AR2ZhdwprjKu1xRFExj0nBVleHM2RHAhEI0nRQU0WxWu2f6X80gc6H4IuXqXhgNOtbtOFFrub77YExxfXHuXY+muQOUToimk4Iimq2mhYJJCoovIkJNXDcqge0HznK6sNK7XFEUuoybhbvsLK4zhwKYUIimkYIims1qVzErDrnCqxFmTuyDOcjAh1sO1xp5OzRxPBhDcP7w7wCmE6JppKCIZqtpoeiCZHItX4WFGEn9UV8OHi9lT86FUbd1xiCM/a9GPbrLO/qAEO2FFBTRbFaHC7POgV4m12qU60Yl0C3azKrPc1BdF86ZGK76Eah21OO7A5hOiMaTgiKazWZzEKw4MYZIQWkMg17HT68fSEGplS++y/Uu18cNQDFHylTBot2RgiKazWn1nFiWLq/GG94visG9u/KPfx+j2uYEQFF0GPqOQT2VJd1eol2RgiKazXV+ci1FCkqjKYrCT68fQLVN5Z8XXUZs6JsMLifqqX0BTCdE40hBEc2meQuKXOXVFL3iwhk/NJ5/ZZ7mXIUdAH38IJTgcNSj0u0l2g8pKKLZ3Pbzd3xLC6XJZl3TB1V184+vjwCg6HQY+oxGPbVPhmIR7YYUFNFsmsMz/a9ikhZKU8VHmUlOjOWf2y6cSzH0TQanDdfpAwFOJ4RvpKCIZtM5zxcUaaE0y4wJvam2qWzZ7bniS999MASF4jz6bYCTCeEbKSii2Wpma5RzKM3TKy6cMYmxbNl1CqfqQtEbMPQejXpij3R7iXZBCopoFrdbw+iy4Vb0oDcFOk67N3vKAMqrnWzbfxYAY/+rwWmVbi/RLkhBEc1SfX7YFVUfjKIogY7T7g3vH0Pv+HA2fnsKt6ahTxgi3V6i3ZCCIpqlyubErDhwG6S7qyUoisK0cb3IL6lm7+EiFJ0BY5/RqCe+k24v0eZJQRHNUmVVMescIFd4tZgxgyzEdAlm/beeWR0N/caC04Z6en+AkwnRMCkoolmqbU5CFIeckG9Bep2OqVf3JOd0GTm5ZegTBntucvxhW6CjCdEgKSiiWapsMnS9P1yb1I3QYAMb/3MSRWfAcNU1qCe+w11VGuhoQtRLCopoFs85FDuGECkoLSnYZGDK6AR2/1BIfkk1psHXgebGeejrQEcTol5SUESzVFvthCgODKFdAh2lw/nx6B4YjTo++fIIui7x6HsMw3ngX2hOe6CjCVEnKSiiWZzVFegUMJgjAh2lw+kSFsT08b3ZdaiQ7BOlBI1ORbOW4zjwr0BHE6JOUlBEs6hVFQAowTK5lj/cNLYXMV2C+cu6bOyRfdD3TMKx559yLkW0SVJQRLO4bTUFJTzASTomk1HP/alDKa2w879rD2IYdzu4VGxb/4zmVgMdT4hapKCIZtHZPbM1SgvFf/p378JdU69i35Fi3vxXIcrYO3DlHsD2r2Uyo6NoU/xaUNauXcv06dOZOnUqK1asuGx9dnY2c+bMISUlhSeeeAJV9XzjysvL46677uKmm27iwQcfpKrKM9/Gzp07GTduHKmpqaSmpvLYY4/5M77wgeI4P7mWtFD86rpRCaTfNIiDx0t4+isThQNuRj2xm6rVj+M89DWa2xXoiEL4r6Dk5+fz6quvsnLlStasWcOqVavIycmptc3ChQtZtGgRGzduRNM0Vq9eDcCzzz7LnXfeyYYNGxg2bBhvvfUWAFlZWdx3332sWbOGNWvWsGTJEn/FFz4yqDUFRVoo/nbdyAQe//kYzEEGnvs2ks/CbsUVFI7ty3eo/ugJnMdkdkcRWH4rKNu3b2f8+PFERkZiNptJSUlhw4YN3vW5ubnYbDZGjhwJwOzZs9mwYQNOp5OdO3eSkpJSazl4Csq2bdtIS0vjgQce4MyZM/6KL3ygaRoGtRpVMaIYZKTh1tC3WwRP33s1t00ZwFd5ofzm+HXs73kbmk6PbfMb2Lb9Dc3tDnRM0UkZ/LXjgoICLBaL93FsbCz79u2rd73FYiE/P5/S0lLCwsIwGAy1lgOEh4czY8YMbrjhBv7v//6PBQsW8OGHH/qcKTrat2/RFkvb7r5pK/mqbU5CsOE2hdXK1Fby1acj5Pv5zC5M+1E//vyP/Szfm8fEYXfyi/4Hqdz1T0JCTETf+Au/jf7cEY5fIHXkfH4rKJqmXbbs4l/w+tY39LzFixd7l91xxx288sorVFRUEB7u2wEoLq7E7b58/xezWMIpLKzwaX+B0JbyFZRWE6bYUQ1mb6a2lK8uHS3fL6Yl0jMmlA+3HMblHsLdw1XKd63HHhLnubs+wPlam+Rrnkvz6XSKz1/EwY9dXnFxcRQVFXkfFxQUEBsbW+/6wsJCYmNjiYqKorKyEpfLVWu52+1m2bJl3uU1aloyovVVVDsJ1dlAzp8E1I1X92T2pH7852A+nzMefcJQ7DtW4i7LD3Q00cn4raBMnDiRHTt2UFJSgtVqZdOmTUyaNMm7PiEhgaCgIDIzMwHIyMhg0qRJGI1GkpOTWbduXa3lOp2OzZs3s3HjRu/yESNGEBIS4q+3IK6gotpJqGJHJ1d4BdyMCb1JTozl718fp2jIT0Gnx/bv9+ts8QvhL35toSxYsID09HTS0tKYOXMmSUlJzJ07l6ysLACWLl3KkiVLmDZtGlarlfT0dACefvppVq9ezfTp09m1axe//vWvAXjxxRd5//33mTFjBp988gnPPfecv+ILH1RUOwjT2WTYlTZAURTSUwYRZjay/PMzGEal4co9gOvEnkBHE52IonWirzByDqVlrd+ew4/2P4d+1E8wX50KtK18deno+fYdKeK1j/aRdk0vfnzmXTS3i9Bbn0PRG9tEPn+TfM3TZs+hiI7PXlkGgCFUurzaiqT+MYwfEsdn35yicuhstPJ8nPs3BzqW6CSkoIgmc50vKIpZhq5vS27/8UCCjHre/U5B33ME9t1rcVvLAx1LdAJSUESTaVZPQdGZIwOcRFwsItTEbdcP4PDpMvZ1uQ5UO47MjEDHEp2AFBTRZDpbTQtFCkpb86Ph3UjsFcnf/lOOe8AknNlbcZXmBjqW6OCkoIgm09s93ShKiHR5tTWKonD3TYk4VY2PSwaDMQj7N6sCHUt0cFJQRJO43G6C1Eoc+hAUvdxc2hbFRZlJ/VEfth2uoqjnj3Gd2od6KivQsUQHJgVFNEl5lZMIXTUuk9yD0paljO1FgiWUZd/HQrgF+zcfylD3wm+koIgmKa2wE6GzogVLd1dbZtDruGdaIsUVKjtDrsVdmovzwJZAxxIdlBQU0STnKu1EKFZ0YXJCvq3r370L14/pwYrvQ7HFJGLf+QnuyuJAxxIdkBQU0STnKqxE6KwEhUcFOorwwexJ/YiPDuWPp5PQNE3G+RJ+IQVFNEnVuVL0ikZQl+hARxE+CAky8PCcJIpdYWzVknGd3It69NtAxxIdjBQU0SRqeSEA+i6WK2wp2or4KDP33zyUtUX9KNLHYdu2Aret7Y4rJdofKSiiSbQKz1w2SlhMgJOIxkjqH83tUxN5p2Qsblsltq//Kl1fosVIQRFNoqsu8fw/TLq82psfj+nBiKtH8M/qEbiO7ULN2RHoSKKDkIIiGs3ldhPsPOe5qdEkE5y1R3Mm96eq3/UcdVqo+up93JUlgY4kOgApKKLRisvtdFUqUYO6BjqKaCKdonDvjKF8EzkdVVUp2vAn6foSzSYFRTRa4TkrUfoqOX/Szhn0OtJvuZav9BMJKfmB/O1rAh1JtHNSUESjFZZWEaWrxNQ1NtBRRDOFBBmYctudZLv7ELx/DUU5+wMdSbRjUlBEo1UWnsGkuAiJ7RnoKKIFdI0Ipmfqf3FOC8O+ZRkVJUWBjiTaKSkootGcxZ55NQxRCQFOIlpKt24WtEn3E4KN05/8HpvNFuhIoh2SgiIaTSk7A4Cua/cAJxEtqe+QYRQPuY0eWh7ff/gaqktGJRaNIwVFNIrVrhLmLMZuCEcxmQMdR7SwgdemkNv9Bvo7vue71cvlyi/RKFJQRKPkFlXRTX8OV3hcoKMIPxk04y5OdRlNYsU3ZGW8F+g4oh2RgiIaJTe/lO76Ukxx/QIdRfiJoigk3vpfHAsZRt/CL/lh3QeBjiTaCSkoolHKT+ZgUNyE9UoMdBThRzqdnsG3/5ofjIPpdvpfnNrwHprmDnQs0cZJQRGN4i7IAUAfNyDASYS/GY0GBt++gD26YUSe3ErRujfRVGegY4k2TAqK8Fl5lYMYRy5WUxS6EJlLvjMwh5hI+ulDbHGPJTg3kyPvPY3bWh7oWKKNkoIifLY/J5+rjGdQug0OdBTRirqEBzPhtnvIcE/BdTaH0g+fRM07FOhYog2SgiJ8dvbgboIVlcjEqwMdRbSyuCgzP7n7Tv4ZeSdlNo2qz35H0ed/RXNYAx1NtCFSUIRPyqocWIp349CFYOgxNNBxRACEhRiZ//+lcmL0r/jGOQjj4S/I/+tCvt/yD6qq5M56AYZABxDtw7bte5loPImr/xQUvTHQcUSA6HUKU8YNpGL4/2Pn1zvoduwzehz5lKIfNvKFMYmqhLH06tWNvt0iiO0agk5RAh1ZtCK/FpS1a9eybNkynE4n99xzD3fddVet9dnZ2Tz55JNUVlaSnJzMs88+i8FgIC8vj4ULF1JcXEzfvn1ZunQpoaGhlJeX88gjj3Dq1CmioqJ47bXXsFhkTnN/O1NUQcwPa9CMerqOuznQcUQbEG42cX3KZFzuazm95xt0BzYw2boD9eR/OHCkB586enFU6U1cfAx94sPp2y2CPvHhRHcJRpEi02Epmp/GVsjPz+eOO+7g008/xWQycfvtt/P73/+eAQMuXG46c+ZMnnvuOUaOHMnjjz/OsGHDuPPOO7n//vu5+eabmTFjBm+++SbV1dUsXLiQxYsXEx8fz7x588jIyGDr1q289tprPmcqLq7E7W747Vos4RQWVjT5fftba+c7mVvE8bX/yyjdD7jG/JTIMdMa3F6OX/O053yuklwc2V/gOPItOls5GlCsRJFji+aMqwulrlBsxi5EW6KIi4umR0IMfbpH0TU8qFXytQXtLZ9OpxAdHebz8/3WQtm+fTvjx48nMjISgJSUFDZs2MBDDz0EQG5uLjabjZEjRwIwe/Zs/vCHP3Drrbeyc+dO3nzzTe/yn/3sZyxcuJCtW7eyYsUKwFOMFi9ejNPpxGj0bxfMDydLUE/tQ3E5Liw8X4c1PP9Xztcp7aJ1F7Z1X1gHKDXP9W520fbe/Z7f9qJ1mqYRHGzAZnNe/GKXvNaF5ytotaJceHzpcy68E+X8Phx2G1rZWbo7jjNKZ8N61U1YRt+EEPXRRyUQcs3PCJ54J+7CY6insogrOEJMwVGwH76wYcX5/3LAoenJw4SqmHDpTbh1Jtx6Ey7dhX+7dSY0RY+mKKDo0NCBongfg87776AgE3anCw0FFP357XQoKKAont/vmgaS4vnn+aUoChe2qdlE8W5aa52iXNiBcum2DbTAToQYsVpb514etz6Yqq4DQdExqFckEWaT31/TbwWloKCgVndUbGws+/btq3e9xWIhPz+f0tJSwsLCMBgMtZZf+hyDwUBYWBglJSXExfk2rpSvldZiCff+22pX+ejjzfwqbL1Pz+1IKgnFFjWAnjfOod/AJJ+fd/Hxa4skX/P4lC92JAz1fFnUNA23rQq1rBC1vAi3vQpHVSXFhSWUlZRSXVaO22EF1Y7O6cDgqMSoOTHixISKSXGiR0OnyECVjbWk7GbOuiKZdW0/5qUN9+k5zfn981tBqasn7eLKXd/6Kz3vUjqd7xeqNbXL66F5aVQUj0dxOWu+y1DX15ianLW2US5+fMn6i57rfVjznemS/dfsu2tXM6Wlnks1FV2tPdV8xbrwLUy5sEcv7/G66FuZd1PvkwgONhFuCgbABT4309tbk76t6dD59DHQ9cK00V37Q9dGPF3TNE8LXHOD5kZzuz3DwbjdaJoLze0mqquZ4qLy8+s0zzq3C0270DrXuKh3QNPQLtr/xeu0S3oSNK2mba95Owe0i7ZDq/vv2sUiuoRQXtY6l1prhiD+yxwFeC779uXn1ma7vOLi4ti1a5f3cUFBAbGxsbXWFxVdmBmusLCQ2NhYoqKiqKysxOVyodfrvcvB08opKioiPj4eVVWprKz0dqn5U5ewILqE9fb76/jCYgkHU9v9gyOEvyhKzZcyz5eiur5mBkeGY3K23DmZltbWvzA0l9/uQ5k4cSI7duygpKQEq9XKpk2bmDRpknd9QkICQUFBZGZmApCRkcGkSZMwGo0kJyezbt26WssBJk+eTEZGBgDr1q0jOTnZ7+dPhBBC+MZvBSUuLo4FCxaQnp5OWloaM2fOJCkpiblz55KVlQXA0qVLWbJkCdOmTcNqtZKeng7A008/zerVq5k+fTq7du3i17/+NQDz589nz549zJgxg5UrV/LUU0/5K74QQohG8ttlw22RXDbsf5KveSRf80i+5mnuORQZekUIIUSLkIIihBCiRUhBEUII0SI61eCQOp1vYwj5ul2gSL7mkXzNI/mapz3la2zWTnVSXgghhP9Il5cQQogWIQVFCCFEi5CCIoQQokVIQRFCCNEipKAIIYRoEVJQhBBCtAgpKEIIIVqEFBQhhBAtQgqKEEKIFtGphl6pT0ZGBkuXLiU6OhqA6667jgULFpCXl8fChQspLi6mb9++LF26lNDQ0FbPl5mZyQsvvICqqkRGRvLCCy+QkJDAzp07eeihh4iPjwdgyJAhLFmypNXz1Vi7di3Lli3D6XRyzz33cNdddwUsC8Abb7zB+vXrAc/kbL/5zW947LHHyMzMJCQkBICHHnqIqVOnBixjeno6xcXFGAyej+LixYs5efJkmziOH330ER988IH38enTp0lNTcVqtQb0GFZWVnL77bfzpz/9iR49erB9+3aWLFmC3W5n2rRpLFiwAIDs7GyefPJJKisrSU5O5tlnn/Ue59bMt2rVKv72t7+hKArDhg3j2WefxWQy8cYbb/DJJ58QEREBwG233dYqP+tL89X3majvuDZIE9rixYu1tWvXXrZ83rx52meffaZpmqa98cYb2ksvvdTa0TRN07QpU6Zo2dnZmqZp2kcffaQ98MADmqZp2jvvvKP96U9/CkimS509e1abMmWKVlpaqlVVVWmzZs3SDh8+HLA827Zt0376059qdrtdczgcWnp6urZp0yZt5syZWn5+fsByXcztdmvXXHON5nQ6vcva2nGs8cMPP2hTp07ViouLA3oM9+zZo82cOVMbOnSodurUKc1qtWqTJ0/WTp48qTmdTu2+++7Ttm7dqmmaps2YMUP77rvvNE3TtMcee0xbsWJFq+c7evSoNnXqVK2iokJzu93ab37zG+0vf/mLpmmadv/992u7d+/2e6aG8mmaVufPs6Hj2hDp8gKysrLIyMjg5ptv5pFHHqGsrAyn08nOnTtJSUkBYPbs2WzYsKHVszkcDubPn09iYiIAgwYN4syZM97c27ZtIy0tjQceeMC7PBC2b9/O+PHjiYyMxGw2k5KSEpDjVcNisfDoo49iMpkwGo3079+fvLw88vLyWLRoEbNmzeIPf/gDbrc7YBmPHj2KoijMnTuXm2++mQ8++KDNHccazzzzDAsWLCA4ODigx3D16tU8/fTTxMbGArBv3z569+5Nz549MRgMzJo1iw0bNpCbm4vNZmPkyJFA631+L81nMpl45plnCAsLQ1EUrrrqKvLy8gDYv38/y5cvZ9asWSxevBi73d7q+aqrq+v8edZ3XK9ECgqePz4PP/wwa9asoVu3bixevJjS0lLCwsK8TWSLxUJ+fn6rZzOZTKSmpgLgdrt54403uOGGGwAIDw8nPT2djIwMJk+e7FuT1E8KCgqwWCzex7GxsQE5XjUGDhzo/WNy/Phx1q1bx7XXXsv48eN54YUXWL16Nbt27eLjjz8OWMby8nImTJjAm2++yXvvvceHH35IXl5emzqO4PmyYLPZmDZtGsXFxQE9hs8//zzJycnex/X93l26vLU+v5fmS0hIYOLEiQCUlJSwYsUKfvzjH1NVVcXgwYP57W9/y9///nfKy8t56623Wj1ffT/Ppn6eO1VBWb9+PZMmTar13z333MObb77JiBEjUBSFX/7yl3z11VdodQzCrCj+HXa6vnzgaak88sgjqKrK/fffD3j622uKyx133EFOTg4VFYGZXjQQx8sXhw8f5r777uO3v/0t/fr148033yQ6OpqQkBB+/vOf8+WXXwYs26hRo3jppZcwm81ERUVxyy238Ic//OGy7QJ9HD/88EPuvfdeAHr27NmmjmF9v3dt7fcxPz+fu+++mzlz5jBu3DhCQ0NZvnw5vXv3xmAwcN999wXkONb382zq8etUJ+WnTZvGtGnTai2rqKjgvffe8/7h1jQNg8FAVFQUlZWVuFwu9Ho9hYWF3mZia+YDqKqq4sEHHyQyMpJly5ZhNBpxu928/fbbzJs3D71e7922NU461iUuLo5du3Z5HxcUFPj9eF1JZmYmv/rVr3j88ceZMWMGhw4d4vjx495uzJqfdaDs2rULp9PJhAkTvHkSEhIoKirybhPo4+hwONi5cye/+93vANrcMYyLi6vzeF26vDU+v/U5cuQIc+fO5Wc/+xn33XcfAHl5eWzfvp1bbrkFCNxxrO/nWd9xvZJO1UKpi9ls5s9//jN79+4F4IMPPmDq1KkYjUaS78kAGwAAAjZJREFUk5NZt24d4LkSbNKkSQHJuHDhQnr37s3rr7+OyWQCQKfTsXnzZjZu3OjNN2LECO+VGq1t4sSJ7Nixg5KSEqxWK5s2bQrY8QI4c+YM//Vf/8XSpUuZMWMG4PmwvPDCC95zZKtWrQroFV4VFRW89NJL2O12Kisr+fvf/87LL7/cpo7joUOH6NOnD2azGWh7x3DEiBEcO3aMEydO4HK5+Oyzz5g0aRIJCQkEBQWRmZkJBO7zW1lZyS9+8Qvmz5/vLSYAwcHBvPzyy5w6dQpN01ixYkVAjmN9P8/6juuVdKoWSl30ej2vvfYazzzzDDabjT59+vDSSy8B8PTTT/Poo4+ybNkyunXrxu9///tWz3fw4EG2bNnCgAEDSEtLAzz9mcuXL+fFF19k0aJFvPnmm0RFRXlzB0JcXBwLFiwgPT0dp9PJLbfcQlJSUsDyvPPOO9jtdu83a4Dbb7+defPmcccdd6CqKjfeeCMzZ84MWMYpU6awd+9e0tLScLvd3HnnnYwZM6ZNHcdTp055L0sHSExMbFPHMCgoiN/97nc8/PDD2O12Jk+ezE033QTA0qVLefLJJ6mqqmLIkCGkp6e3er6PP/6YoqIi3n33Xd59910Arr/+eubPn8/ixYt58MEHcTqdjB492tut2Joa+nnWd1wbIjM2CiGEaBGdvstLCCFEy5CCIoQQokVIQRFCCNEipKAIIYRoEVJQhBBCtAgpKEIIIVqEFBQhhBAtQgqKEEKIFvH/Aw804Q+VKQKYAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"df_sub_F.age.plot(kind='density')\n",
"df_sub_M.age.plot(kind='density')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The histogram, however, doesn't gives us precise information about gender difference, for this we need to use the bootstrap procedure."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Bootstrapping looks as follows:\n",
"\n",
"we will draw samples of size n from our data (effectively we treat our data are the population from which we take repeated samples)\n",
"for each sample, we compute and collect the statistic of interest (the mean) and put the observation back in our data (which is called with replacement)\n",
"We repeat steps 1-2 R times. \n",
"\n",
"This procedure will generate a distribution of the sample statistic: it conveys which values are consistent with the data, and which ones are unlikely to occur. Our estimate of the mean may contain error, and we'd like to compute how much variation we could expect due to random chance.\n",
"\n",
"\n",
"In code, it is easy to implement this procedure. In the `for` loop, we randomly sample 100 observations, compute the mean, and store this statistic in `mean_sampled`. We repeat this procedure 1000 times. "
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAD7CAYAAACrOanfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU5b0/8M85Z5asEBImIUZkESUqm5jKog1SwRQkgBR7EW3K9RKrrWDzuqDUjeWn4gVa1CvQ22iv9yd4FRVC4w8j2N62/gxtSawQlVKQgkIgOyQTJjNnu38MMyRkZjIZc2bj8369eMGcZfLlZPnkeZ7zPEfQdV0HERFRCMRIF0BERLGLIUJERCFjiBARUcgYIkREFDKGCBERhYwhQkREIWOIEBFRyEyRLsAILS3t0DT/018yMlLQ1GQPY0W9E+31AdFfY7TXB0R/jdFeHxD9NUZ7fYC7xpaWdgwYkBzS+XEZIpqmBwwRzzHRLNrrA6K/xmivD4j+GqO9PiD6a4z2+oBvViO7s4iIKGQMESIiChlDhIiIQsYQISKikDFEiIgoZAwRikt8wgFReDBEKO786YszeGDthzjX7op0KURxjyFCcWfHH47hTNN5HP6qJdKlEMU9hgjFndYLLZDmVmeEKyGKfwwRiiut511wKRoAoKWNIUJkNIYIxZWGFof333aHHMFKiC4PDBGKK02tHd5/d7iUCFZCdHlgiFBc8YyDDM3uB4eTIUJkNIYIxZWm1g4kWiUMTEuEw6VGuhyiuMcQobjS3NqB9H4JSEowsSVCFAaGhkh5eTlmzpyJ6dOnY9u2bd32f/jhh5gzZw5mz56NH//4xzh37hwAoLa2Fvfeey+++93v4qGHHkJ7e7uRZVIcOWt3Ii3FiuQEM0OEKAwMC5G6ujps3LgRb7zxBnbt2oW33noLR48e9e632+1YtWoVfvWrX+E3v/kNRo4ciX//938HAKxevRoLFy5ERUUFRo0ahc2bNxtVJsUZu0NGSqIZVosEl6xFuhyiuGdYiFRWVmLixIlIS0tDUlISCgoKUFFR4d0vyzJWrVqFrKwsAMDIkSNx+vRpyLKM/fv3o6CgAAAwb968LucRBWJ3KEhJMMNqluCSVa6hRWQwwx6PW19fD5vN5n2dmZmJgwcPel8PGDAA06ZNAwB0dHTgV7/6FX7wgx+gpaUFKSkpMJncpdlsNtTV1fXqY2dkpPR4jM2W2qv3DLdorw+IvhpVVYPDqSBzYDJMkggdwID0ZJhNUqRL8yvaruGlor0+IPprjPb6gOB+ZvpjWIj4+g1QEIRu29ra2vDjH/8Yubm5uOuuu3wGhq/zAmlqsgd8ZrDNloqGhrZevWc4RXt9QHTW6FnuRNA0WBLMAIDa0+eQdOHf0SYar2Fn0V4fEP01Rnt9gLvGpiZ7yEFiWHdWVlYWGhsbva/r6+uRmZnZ5Zj6+nosXLgQubm5ePbZZwEA6enpsNvtUFX37ZkNDQ3dziPyxTNDPSXRDKvZ/aXtWQKFiIxhWIhMnjwZ+/btQ3NzMxwOB/bs2YP8/HzvflVV8eCDD2LGjBl44oknvK0Ns9mMvLw87N69GwBQVlbW5Twif9o7LoaIxezuwmKIEBnLsO6srKwslJSUoKioCLIsY/78+RgzZgyKi4uxdOlSnDlzBl988QVUVcUHH3wAABg1ahSeffZZrFy5EitWrMCWLVuQnZ2NX/ziF0aVSXHE4XS3XhOsJihw/1Iiy5xwSGQkw0IEAAoLC1FYWNhlW2lpKQBg9OjR+Nvf/ubzvJycHLz++utGlkZxyHUhMKxmCdKFEGFLhMhYhoYIUTg5vSEiQvSECFsiRIZiiFDccHZqiYgXhvtktkSIDMUQobjhDRGLBEFndxZRODBEKG44L6zaazFLgDdE2J1FZCSGCMUNl6zBYhIhCgKsF2apy1w/i8hQXAqe4kaHrHrnh5hM7i9tJcDKBUT0zTFEKG44XSqsF0LE7AkRjokQGYohQnHDJauwWi6EiORpiTBEiIzEEKG44ZRV75pZJrZEiMKCIUJxwx0i7paIJAoQAMgqx0SIjMQQobjhkjXvwLogCJAkEarKlgiRkRgiFDdkVfMOqAOA2SRAZogQGYohQnFDVlTvgDoAmCQRCruziAzFEKG4oag6TN1ChC0RIiMxRChuyIrmvSsLAEySwBAhMhhDhOKGomowSYL3NbuziIzHEKG4oaha9zERzhMhMhRDhOKCruscEyGKAIYIxQVPtxXHRIjCiyFCccETFrzFlyi8GCIUFzyTCjsPrJtN7M4iMhpDhOKCZwC9c3eWJLI7i8hoDBGKC766s8wmkQswEhmMIUJxwRMWne/OkkQuwEhkNIYIxQVvd5bEBRiJwokhQnHB251lumTGOicbEhmKIUJxQVG7t0RMkghF45gIkZEYIhQXZH8hwpYIkaEYIhQXFMXd4jBfMmNd1XToOlsjREZhiFBc8Ned5d7HECEyCkOE4oKvGesXQ4RdWkRGYYhQXPCMfVzanQWAt/kSGYghQnHB5wKMFwJFZXcWkWEYIhQXZF9LwYueEGFLhMgoDBGKC74H1t3dWZwrQmQchgjFBe+YiK+7szhXhMgwDBGKC7KqQRQEiOLFu7Mkb0uEIUJkFENDpLy8HDNnzsT06dOxbds2v8c99thj2LFjh/d1WVkZbr31VsyZMwdz5szBxo0bjSyT4oCiajB1WjcLuNgq4TwRIuOYjHrjuro6bNy4ETt27IDFYsGCBQswYcIEjBgxossxK1euxL59+zBhwgTv9pqaGqxYsQKzZs0yqjyKM4qid+nKAgBJ4sA6kdEMa4lUVlZi4sSJSEtLQ1JSEgoKClBRUdHlmPLyctx+++2YMWNGl+01NTUoKyvD7NmzsWzZMpw7d86oMilOyKrWZVAd6DSwzpYIkWEMC5H6+nrYbDbv68zMTNTV1XU5ZvHixbj77ru7nWuz2bBkyRLs2rUL2dnZWLNmjVFlUpxQfIaI+zUnGxIZx7DuLF+L3gmC4OPI7jZt2uT99+LFizFt2rRefeyMjJQej7HZUnv1nuEW7fUB0VWjZJKQYJW61GQb6P46SE62RlWtnUVrXR7RXh8Q/TVGe31AcD8z/TEsRLKyslBVVeV9XV9fj8zMzB7Pa2trw7vvvotFixYBcIeRydS7Mpua7NACzA2w2VLR0NDWq/cMp2ivD4i+Gu3tTgiAtyabLRVtrQ4AQHPL+aiq1SParuGlor0+IPprjPb6AHeNTU32kIPEsO6syZMnY9++fWhubobD4cCePXuQn5/f43lJSUl45ZVXcODAAQDA1q1bMX36dKPKpDihqHq37iyJCzASGc7QlkhJSQmKioogyzLmz5+PMWPGoLi4GEuXLsXo0aN9nidJEl544QWsWrUKHR0dGDp0KNatW2dUmRQn3Lf4XjImcmHOiMoZ60SGMSxEAKCwsBCFhYVdtpWWlnY77vnnn+/yOi8vDzt37jSyNIozsqp1u8XXO7DOGetEhuGMdYoLiuL/Fl/OEyEyDkOE4oL7Ft+ud/95185idxaRYRgiFBcUVe/yQCqATzYkCgeGCMUFRdUgiV2/nEVRgCBwxjqRkRgiFBcUVYPZ1H0yq0kS2RIhMhBDhOKCoureeSGdmSSBIUJkIIYIxQVF1byPw+3MJIl8xjqRgRgiFBcUVe/2PBGA3VlERmOIUMzTdd1vS0QS2Z1FZCSGCMU8z7Imly57AnhaIuzOIjIKQ4RinqelcelkQ882tkSIjMMQoZjnaWn4HVjnjHUiwzBEKOZ51sby353FlgiRURgiFPM8j7/1LP3emUkSoHAVXyLDMEQo5nnmgfhqiUiSyAUYiQzEEKGY522J+Jqxzlt8iQzFEKGY522J+OrOMnHGOpGRGCIU8xQOrBNFTFAhsmTJElRWVhpdC1FIlEAD6+zOIjJUUCFyxx13YPPmzSgoKMCrr76Ks2fPGl0XUdCUngbW2Z1FZJigQqSwsBBbt27F5s2b0dTUhLvvvhvLly/HwYMHja6PqEcBB9Y5Y53IUEGPiWiahhMnTuD48eNQFAUZGRlYtWoV1q9fb2R9RD3yTDaUfM4T4S2+REYyBXPQxo0bsWPHDgwePBgLFy7Eiy++CLPZjPPnz2Pq1KlYvny50XUS+eXprrr0GeuA53kibIkQGSWoEGlubkZpaSlyc3O7bE9KSsLPf/5zQwojCpanu8r/kw116LoOQejeUiGibyao7ixVVbsFyJIlSwAAt956a99XRdQLnhAx+wgRT7BwEUYiYwRsiaxcuRJ1dXWorq5Gc3Ozd7uiKDh27JjhxREFw9OdJflZCt59jOZz4J2IvpmAITJ//nwcOXIEhw8fRkFBgXe7JEm48cYbDS+OKBgX54n4HhNxH8OWCJERAobI6NGjMXr0aNxyyy3IysoKV01EveLtzvLzjHUAHFwnMkjAEHnkkUfw4osvYvHixT73l5eXG1IUUW9c7M7yvQAjcHEuCRH1rYAhUlxcDAB46qmnwlIMUSgUVYMkChB93H11sSXC7iwiIwQcaRw1ahQA4Oabb0Z2djZuvvlmnD17Fn/5y19w3XXXhaVAop4oquZzUB24uBQKZ60TGSOo21WefvpplJaW4ssvv8SaNWtw6tQpPPHEE0bXRhQURdV9DqoDF7uzOLBOZIygQuSzzz7DqlWrsHfvXtx1111Yu3YtTp06ZXRtREFRVc3n4ovAxXESRWNLhMgIQYWIrusQRREff/wxJk6cCABwOByGFkYULFnVvPNBLuWdJ8LnrBMZIqgQueqqq1BcXIyTJ0/i5ptvxr/+679i5MiRRtdGFBRV1f1OJPTOE+GMdSJDBLV21tq1a7F3717cdNNNMJvNyMvLw9y5c42ujSgocoDZ6JwnQmSsoFoiSUlJyMvLQ2trKz7//HOMGTOGy55Q1FBV3edTDYHOy56wJUJkhKBaIuvXr8fWrVuRkZHh3SYIAn77298aVhhRsJRgBtbZEiEyRFAh8v7772PPnj29XvqkvLwcW7ZsgSzLWLRoEe69916fxz322GOYMGEC5s2bBwCora3F8uXL0dTUhGHDhmHDhg1ITk7u1cemy4eiakG0RBgiREYIqjsrOzu71wFSV1eHjRs34o033sCuXbvw1ltv4ejRo92OefDBB1FRUdFl++rVq7Fw4UJUVFRg1KhR2Lx5c68+Nl1eFFX32xLxzB9hdxaRMYIKkUmTJmHdunWorq7G559/7v0TSGVlJSZOnIi0tDQkJSWhoKCgW1iUl5fj9ttvx4wZM7zbZFnG/v37vasGz5s3r9t5RJ0FHFg3cWCdyEhBdWft2LEDALr8MO9pTKS+vh42m837OjMzEwcPHuxyjGdhx+rqau+2lpYWpKSkwGRyl2az2VBXVxdMmV4ZGSk9HmOzpfbqPcMt2usDoqdGQRCQlGjuVo/NlorEZBcAICHREjX1dhaNNXUW7fUB0V9jtNcHBPcz05+gQuR3v/tdr99Y17t3HwTzeNJQz+usqckOLcC8AJstFQ0Nbb16z3CK9vqA6Kqxw6lAU7Uu9Xjqc7pUAMDZVkfU1OsRTdfQl2ivD4j+GqO9PsBdY1OTPeQgCao7q729HWvWrMEPf/hDnD17Fk8//TTa29sDnpOVlYXGxkbv6/r6emRmZvb4sdLT02G326Gq7m/+hoaGoM6jy5d7FV9/d2dxxjqRkYIKkWeeeQapqaloamqC1WqF3W7H008/HfCcyZMnY9++fWhubobD4cCePXuQn5/f48fyTGbcvXs3AKCsrCyo8+jypaiazwdSAYDEBRiJDBVUiBw6dAglJSUwmUxITEzEhg0bcOjQoYDnZGVloaSkBEVFRZg7dy5mzZqFMWPGoLi4GDU1NQHPXblyJbZv346ZM2eiqqoKP/3pT4P/H9FlR1F1nw+kAtxdoSZJ5AKMRAYJakxEvKSrQFXVbtt8KSwsRGFhYZdtpaWl3Y57/vnnu7zOycnB66+/HkxpRBfmifj/ejRJAh9KRWSQoELkW9/6FtavX4+Ojg589NFH2Lp1KyZMmGB0bURBcc8T8X/zhUkSOdmQyCBBdWctW7YMSUlJSE1NxQsvvIDc3Fw8+uijRtdG1CNd13tsiUiSwBAhMkiPLZG9e/fi1VdfxeHDh5GQkICRI0di/PjxsFqt4aiPKCD1wq3c/p4nArhnrXNgncgYAUOkrKwMmzdvxtKlS5GbmwtBEFBTU4Nnn30WTqcTd9xxR7jqJPLJ08Lwt+yJZx9bIkTGCBgir7/+Ol577TVcccUV3m1XX301xo4di8cff5whQhHnaWFwYJ0oMgKOiciy3CVAPIYNGwan02lYUUTBUoNpiYhsiRAZJWCISJLkd5+v5UmIwk32hIifpeABd0uEIUJkjKDuziKKVp5uKn+r+ALuB1NxYJ3IGAHHRA4fPozx48d3267rOlwul2FFEQVLDqI7yywJ6JDVcJVEdFkJGCJ79+4NVx1EIfG2RAJ0Z0mSCKVDCVdJRJeVgCGSk5MTrjqIQhLULb6SyIdSERmEYyIU05QgB9ZljokQGYIhQjHNM2DubxVfAJBEtkSIjMIQoZjmGVg3BxpYN/EWXyKjMEQopnlaGFJPA+vsziIyBEOEYponHAK1REyiCJUPpSIyBEOEYpqnmyrQmIhJEiArbIkQGYEhQjEtmLuzpAu3+HKpHqK+xxChmOZdxbeHGes6AI0hQtTnGCIU02Tlwt1ZAbuz3Ps4uE7U9xgiFNOCucXXM17CuSJEfY8hQjFNVjQICHyLr+fRuZy1TtT3GCIU0xRFg9kkQhAChQhbIkRGYYhQTJNVLeCzRICLLRHOWifqewwRimnyhZZIIBxYJzIOQ4RimqL2HCKSKHqPJaK+xRChmCYrvenOYkuEqK8xRCim9a47iy0Ror7GEKGYFkx3lqclwruziPoeQ4RiWnDdWRdaIhq7s4j6GkOEYpocVEvkQogobIkQ9TWGCMU0RdECrpsFAJJnYJ0tEaI+xxChmCarWsAVfAEOrBMZiSFCMU0OoiXCGetExmGIUExzj4n4XzcL6Lx2FruziPoaQ4RimntMRAp4jCdEZLZEiPocQ4RiWjB3Z3mWiWdLhKjvGRoi5eXlmDlzJqZPn45t27Z123/o0CF873vfQ0FBAZ544gkoigIAKCsrw6233oo5c+Zgzpw52Lhxo5FlUozSdd09TyTI7iyOiRD1PZNRb1xXV4eNGzdix44dsFgsWLBgASZMmIARI0Z4j1m+fDmeeeYZjBs3Do8//ji2b9+OhQsXoqamBitWrMCsWbOMKo/igKrp0PXAj8YFOLBOZCTDWiKVlZWYOHEi0tLSkJSUhIKCAlRUVHj3nzp1Ch0dHRg3bhwAYN68ed79NTU1KCsrw+zZs7Fs2TKcO3fOqDIphnlCoadbfAVBgCQKUDlPhKjPGRYi9fX1sNls3teZmZmoq6vzu99ms3n322w2LFmyBLt27UJ2djbWrFljVJkUw+QLM9B7aokA7i4tmTPWifqcYd1Zut79t77OjzANtH/Tpk3ebYsXL8a0adN69bEzMlJ6PMZmS+3Ve4ZbtNcHRL5G4awDAJA+IMlnLZ23mU0iLBZTxGu+VLTVc6lorw+I/hqjvT4guJ+Z/hgWIllZWaiqqvK+rq+vR2ZmZpf9jY2N3tcNDQ3IzMxEW1sb3n33XSxatAiAO2xMpt6V2dRkhxag68JmS0VDQ1uv3jOcor0+IDpqPNN8HgDg7JC71XJpfaIooK3dGfGaO4uGaxhItNcHRH+N0V4f4K6xqckecpAY1p01efJk7Nu3D83NzXA4HNizZw/y8/O9+3NycmC1WlFdXQ3AfUdWfn4+kpKS8Morr+DAgQMAgK1bt2L69OlGlUkxzOlSAQBWc+B5IgBglgQOrBMZwNCWSElJCYqKiiDLMubPn48xY8aguLgYS5cuxejRo7FhwwY8+eSTaG9vx/XXX4+ioiJIkoQXXngBq1atQkdHB4YOHYp169YZVSbFMJfiDhGLueffhSRJ5DwRIgMYFiIAUFhYiMLCwi7bSktLvf/Ozc3FO++80+28vLw87Ny508jSKA64ZHfLIpiWiEkSOWOdyACcsU4xyylfaImYgggRUWBLhMgADBGKWS65d91ZHBMh6nsMEYpZLiX47iwOrBMZgyFCMctzd5YliBCRJJFPNiQyAEOEYpbn7qxgB9b5jHWivscQoZjllDUIwsUFFgMxSQIUDqwT9TmGCMUsl6zCYpa6LKfjjySJUDW2RIj6GkOEYpZLVoPqygI4sE5kFIYIxSynrMHSwzLwHu5bfNmdRdTXGCIUs3rTEjFxngiRIRgiFLMcLgUJ1iBDRGR3FpERGCIUsxxOBYnW4JZ/M5v5UCoiIzBEKGadd6pICjJErGYJiqqzNULUxxgiFLN60xJJuDB24llvi4j6BkOEYlZvQsRicYeIU2ZLhKgvMUQoJimqBlnRkGgJbmDd0xLpcClGlkV02WGIUExyON1hEGxLxOrtzmJLhKgvMUQoJvU6RCxsiRAZgSFCMcnhdA+QB313FsdEiAzBEKGYdP5CSyTBaoKuuKCrCnTd/7Imnu4sJ+/OIupTwf0aRxQkvcMO5eRn0NubIST2g5R5NcS07D57f+38OSjH9qP/kRo80f8YbB++Cbvqcu8UJQgpAyH2z0TzkJFQUgZDGnQtBEuid2Dd8yArIuobDBHqE7quQ/78Qzj/8g6gOLvsE9OyYRoxEeYRkyH2s/X+vWUnlOPVkI/ug3ryc0DXYLWk4Zg6ABlXX4vE/v0BXQfkDmhtjdDOncbZj3cAugYIEqTsa2HJHoUBosqWCFEfY4hQn3D+6U3INR9AGjwG1pvmQkzLhna+BeqpL6Acq4KraidcVTshDboWpmsmwzz8WxCsyX7fT1cVqKc+g3z0z1COfwIoTggpGbCMnQnTNZPw0WEF7/z+S2y5ZYrPRRgz+ptQd6gG6tc1UL46AKH6bTzdHzh3+FPI/abBNCwPgsli5CUhuiwwROgbc33xP5BrPoD5hmmwTl4IQXAPtUmWREhpV8BywzRo9ibIR/ZBOVIJ50evwfnxVpiuGgPRNhxi/0wIJgt0Vwc0eyO0+mNQag8BLgdgSYJ5xCSYrpkEadA13ve2nz8Ki0n0u4qvaEmE6YrrYLriOlgnfB/quTq8+9obmGI5gY7/+RWEyjdgGvltWK6fCrFfZtiuFVG8YYjQN6K11sP5p/+GdOUoWCddDJBLiSkZsN44C5Zxd0JrPAH5yMdQjn/ibmVcQuifBfOwPJiG5UHKuQGC1P3LtM3hQkqSOeg6pf5Z+Ei/Ca4rZmDBKA3yF7+DXPMB5IMVkAaPhuWG70C6cgwEkfeaEPUGQ4RCpus6Oj76L0AQkZB/f1A/gAVBgGQbCsk2FJh8L3SXA5q9EVBcgCkBYnJawG4uj3aHgpTE4EMEAFKSLLB3KDDljIIp53po7S2QD/0e8t/+AEfFCxBSB8J89URIg0ZATLsCQlIau7yIesAQoZCpJ2ugnvoc1sn3QkxJD+k9BEsipPTBvT6vzeFCai9DJDXRjLbzsve1mDwA1ry7YBlfCOX4XyF/8Tu4DuwGPr04l0SwprjDJKk/hOQ0SBlXQRp0LcSBQ4N6tjtRvGOIUEh0XYPzL++6f3u/bmrYP779vIyB2Ym9Oicl0Yym1o5u2wXRBPPwb8E8/FvQZSfUxuPQW+uhtbdAb2+Bfv4sNMc5aCdPQfn7x+5z+mXCfO2tsFz/HQgJKX3yfyKKRQwRCkn7oX3Qmk4g4bZin2MWRrM75BC6s8w4UdcW8BjBbIUpeySQPdLnfq29BerJzyAf3QdX1Q64DuyG5YZpsIz5LsOELksMEeo1XVPR8oc3IQ7IgWnEpLB/fFXT0N7R+zERd3eWC7quh9wVJSYPgDjy2zCP/DbU5q/h+qQcrk//H1yff8gwocsSQ4R6TTlSCbm5Fgl3LInI3Uyt7e5xjf7JvRv0Tu+XAEXVcdbuwoBU6zeuQ0ofjMRpP4bafBKuT34D16fvecPEnDslpImVRLGGIUK9oqsKnNVlsGZfDdOQ8RGpoaXNPSM+rZdBMCg9CQBQ13y+T0LEQ0q/8kKYzIbrk13uMPn0PYiZw2HKzoU06BqIadkQUgdCEPktR/GFX9HUK/LhP0K3N2HArIfQHqG7kzwhMiAltBA53XweuUMG9Hld7jD5CbS2RshH90E58SlcBz8ADux2HyBIEPoNhDTgSogZV+H8iOuhJ18FwdR3gUYUbgwRCpquuOD6aznErBFIHD4O7Y32iNRx1h5aSyS9nxX9ki04cvIspt6YY0RpAAAxdSCsNxbCemMhdMUJtfEr6OfOQDtXB+3cGajNJ6Ec/wRnqncCkhnSFdfBNHQ8TEPHQ0zsZ1hdREZgiFDQXAffh97egoSpD0R0jsRZuxOSKCC1FzPWAfdEx+uGDMAXx1ug6TrEMPwfBJMVpkHXAIOu6bJdlzuQ0nEKTTV/hnLir+6lYP7/f7nXFhs6HqahN0FMHWh4fUTfFEOEgqK1NcD11/dgGv4tmK64LqK1NJ7rwIBUa0ghMHp4Ov78RR3+cboVV1/R34DqgiOYE5B0xTi0p14NfdI90JpPQvlHFZTjn8C577/h3PffEAcOgWQbDjFtEISkARCsye4/CckQEvtzNj1FBYYI9UjXNHT8TykgSrBOXBDpcnCm6bx3fKO3xo4YCEkU8MnhhoiGSGeCIEDKGAwpYzCseXdBO1cH5Xg1lBOfQv7yz4DrvO8TLYkQkzMgpl8JMeNKSOlXQRx4FcSktPD+B+iyxhChHrn2vwP1zN+RcFsxxJSMiNai6zrOtJzHiCtDe9BVcoIZuUMGoPrvDZh/29VRuXSJ2D8LlrEzYRk7E7quQ+9og+5ohe5sd//xvD5/DlpbA9S6I1C+/JP3fCEpDWLGVZAGDnG3ZgYOgZAyMCr/rxT7DA2R8vJybNmyBbIsY9GiRbj33nu77D906BCefPJJ2O125OXlYfXq1TCZTKitrcXy5cvR1NSEYcOGYcOGDUhO7nlRPupbuq675z8c2A1z7m0wX3tLpEvCmebzcLpUDM4MfULfTSNt+L8Vh3GqoR1XfoP3CQdBECAk9gN6GHDXne1Qm09Ca6moiGMAAA36SURBVDwBtfEEtKYTcJ38zP1gLgCwJkNMy4aUlg0xLRtif/ffQj8bBNH3cvpEwTAsROrq6rBx40bs2LEDFosFCxYswIQJEzBixAjvMcuXL8czzzyDcePG4fHHH8f27duxcOFCrF69GgsXLsSdd96JTZs2YfPmzVi+fLlRpZIPauMJOP/yNtSTn8F0zWRYby2KdEkAgEMnWgAA11wZelfUjdfY8HrFYVQdro/6EAmWYE3utlyLrrigNZ/0hop29jSUrw5AP/zRxRNFCUJKBsTkdAgp6Rf/TkmHU8mB1iFCsKYAJgtbMuSTYSFSWVmJiRMnIi3N3T9bUFCAiooKPPzwwwCAU6dOoaOjA+PGjQMAzJs3Dy+99BLuvvtu7N+/H5s2bfJuv++++wwPEV3XoZ76HLrT0/+se3b4O8P//lDO6bS99ZQVrrbuCwUG/liB37PL/i6HXnyhO1qh2Zuh1h+Ffq4OMCfCOvk+mG+4PSI/QBRVw4GjjehwqZBVDbKsYc/+r3GlLSXkMRHAPdM9d8gA/Lb6JPonW5CYYIIoCBAEAblXpSE1KT4GrAWTBVLmcEiZw7ts153t0M6dgXb2NLSWWvcjhdubodX+Dcr5s97Wy6nOJ0lmCAkp7j/mRHeomKwX/zZb3QP9ogkQBEAQAUGEIF78t/uPAEC48HfnYi++Fi7d7/13921tdYmQ2zq67vPuF7qdFm72xkTIrY4wfTQBpsGjIZgTwvTx3AwLkfr6ethsF5d9yMzMxMGDB/3ut9lsqKurQ0tLC1JSUmAymbps742MjJ5/u7TZUru8djWexMndP4ffH8Zh5CM+wkMQIaUMQELmVUiacCdSRt8GKcF/N+Kl17CvVR6sxaadn3XZ1j/FgkcW3IjMzJ7nUwSq75F7xmN16Z/w+p6/d9k+O384iueMDq3gEBh9DX1LBa4cBGBctz26pkK1n4XS1gTV3gLV0QbtfBtURxvU823QHK3QXA7oshNae6v77wt/dNl5sfssTCL2vRKkcNeXPm0R0iYU9vq8YH5m+mNYiOg+fmvu/Nusv/09nReMpiY7NM1/GNhsqWhouHQ11/5IXrgButz50979Nxmhy681vn7T8WzzV3OAcy6cl5GejKbmdr/7u7+lr48V/Dne/5Ml0ds/7gLQ3KYBbb5XvfV9DfvWNdmpeP5HEwFBgFkSYZIEJFpNMElijx+7p/osAP7P4pvR3NoBWdGg6+6vyUEZSYb/v4KtMXIsgCUbtpHXdqlPACBd+OOPrmvuINE0dyvY81rXL+679Htc19G1pdyp9ezdfOl+9+sBA5LR0my/cITe6VC922mRkJ6ehOZmP3fX9TVBgCstu9dfUzZbKpqa7CEHiWEhkpWVhaqqKu/r+vp6ZGZmdtnf2Njofd3Q0IDMzEykp6fDbrdDVVVIkuTdHg6RvvPIw9Q/FaIrvE3SaJU5IPRuq56IgoCB/Xv3TBIKTPB0W/lYl9OIXiVLRipELRqD2M1iS4WE6K2vLxi2BOvkyZOxb98+NDc3w+FwYM+ePcjPz/fuz8nJgdVqRXV1NQCgrKwM+fn5MJvNyMvLw+7du7tsJyKi6GNYiGRlZaGkpARFRUWYO3cuZs2ahTFjxqC4uBg1NTUAgA0bNmDt2rWYMWMGHA4HiorcdwCtXLkS27dvx8yZM1FVVYWf/vSnRpVJRETfgKD7GoSIcaGNiUSPaK8PiP4ao70+IPprjPb6gOivMdrrA775mEj4nyhERERxgyFCREQhY4gQEVHI4nIBRlHs+WbCYI6JpGivD4j+GqO9PiD6a4z2+oDorzHa6wO+WY1xObBOREThwe4sIiIKGUOEiIhCxhAhIqKQMUSIiChkDBEiIgoZQ4SIiELGECEiopAxRIiIKGQMESIiCllcLnviS1lZGTZs2ICMDPfTC2+77TaUlJSgtrYWy5cvR1NTE4YNG4YNGzYgOdn/c8WNVF1djeeeew6KoiAtLQ3PPfcccnJysH//fjz88MMYNGgQAOD666/H2rVrI1JjeXk5tmzZAlmWsWjRItx7770RqeNSL7/8Mt5//30AwJQpU/Doo4/iZz/7Gaqrq5GY6H564cMPP4zp06dHpL6ioiI0NTXBZHJ/y61ZswZfffVV1FzLt99+G1u3bvW+PnnyJObMmQOHwxHxa2i327FgwQL88pe/xJVXXonKykqsXbsWTqcTM2bMQElJCQDg0KFDePLJJ2G325GXl4fVq1d7r3c463vrrbfw+uuvQxAEjBo1CqtXr4bFYsHLL7+Md999F/369QMAfP/73w/b5/zSGv19b/i7tgHpl4k1a9bo5eXl3bY/8MAD+nvvvafruq6//PLL+rp168JdmtfUqVP1Q4cO6bqu62+//bb+4IMP6rqu66+++qr+y1/+MmJ1eZw5c0afOnWq3tLSore3t+uFhYX6kSNHIl2W/vHHH+v/9E//pDudTt3lculFRUX6nj179FmzZul1dXWRLk/XNE2/5ZZbdFmWvdui9Vrquq7//e9/16dPn643NTVF/Bp++umn+qxZs/QbbrhB//rrr3WHw6FPmTJF/+qrr3RZlvX7779f//3vf6/ruq7feeed+l//+ldd13X9Zz/7mb5t27aw13fs2DF9+vTpeltbm65pmv7oo4/q//mf/6nruq7/6Ec/0j/55BPDa+qpRl3XfX5eA13bQC6b7qyamhqUlZVh9uzZWLZsGc6dOwdZlrF//34UFBQAAObNm4eKioqI1OdyufDII48gNzcXADBy5EicPn3aW/vHH3+MuXPn4sEHH/RuD7fKykpMnDgRaWlpSEpKQkFBQcSuV2c2mw0rVqyAxWKB2WzG1VdfjdraWtTW1uKpp55CYWEhXnrpJWiaFpH6jh07BkEQUFxcjNmzZ2Pr1q1Rey0BYNWqVSgpKUFCQkLEr+H27duxcuVKZGZmAgAOHjyIIUOGYPDgwTCZTCgsLERFRQVOnTqFjo4OjBs3DkD4vpcvrc9isWDVqlVISUmBIAi49tprUVtbCwD47LPPUFpaisLCQqxZswZOp9Pw+nzVeP78eZ+fV3/XtieXTYjYbDYsWbIEu3btQnZ2NtasWYOWlhakpKR4m7w2mw11dXURqc9isWDOnDkAAE3T8PLLL2PatGkAgNTUVBQVFaGsrAxTpkwJrolpgPr6ethsNu/rzMzMiF2vzq655hrvD4/jx49j9+7d+Pa3v42JEyfiueeew/bt21FVVYV33nknIvW1trZi0qRJ2LRpE1577TW8+eabqK2tjcprWVlZiY6ODsyYMQNNTU0Rv4bPPvss8vLyvK/9fQ1euj1c38uX1peTk4PJkycDAJqbm7Ft2zbcfvvtaG9vx3XXXYfHHnsMO3fuRGtrKzZv3mx4fb5q9Pd5DfX7O+5C5P3330d+fn6XP4sWLcKmTZswduxYCIKAxYsX449//CN0HwsYC4Lxyzb7qxFwt0iWLVsGRVHwox/9CIC7/9wTKPfccw+OHj2KtrbwP3IzUtcrWEeOHMH999+Pxx57DMOHD8emTZuQkZGBxMRE/OAHP8Af/vCHiNR14403Yt26dUhKSkJ6ejrmz5+Pl156qdtx0XAt33zzTfzzP/8zAGDw4MFRcw09/H0NRtvXZl1dHX74wx/ie9/7HiZMmIDk5GSUlpZiyJAhMJlMuP/++yN2Lf19XkO9hnE3sD5jxgzMmDGjy7a2tja89tpr3h/Uuq7DZDIhPT0ddrsdqqpCkiQ0NDR4m3zhrhEA2tvb8dBDDyEtLQ1btmyB2WyGpmn4j//4DzzwwAOQJMl7bDgGDC+VlZWFqqoq7+v6+vqwXK9gVFdXY+nSpXj88cdx55134vDhwzh+/Li3q9LzOY+EqqoqyLKMSZMmeWvJyclBY2Oj95houJYulwv79+/H888/DwBRdQ09srKyfF63S7eH63vZly+//BLFxcW47777cP/99wMAamtrUVlZifnz5wOI7LX093n1d217EnctEV+SkpLwyiuv4MCBAwCArVu3Yvr06TCbzcjLy8Pu3bsBuO/gys/Pj1idy5cvx5AhQ/Diiy/CYrEAAERRxN69e/HBBx94axw7dqz3ropwmjx5Mvbt24fm5mY4HA7s2bMnotfL4/Tp0/jJT36CDRs24M477wTg/sZ47rnnvGNfb731VsTuzGpra8O6devgdDpht9uxc+dOrF+/Puqu5eHDhzF06FAkJSUBiK5r6DF27Fj84x//wIkTJ6CqKt577z3k5+cjJycHVqsV1dXVACL3vWy32/Ev//IveOSRR7wBAgAJCQlYv349vv76a+i6jm3btkXsWvr7vPq7tj2Ju5aIL5Ik4YUXXsCqVavQ0dGBoUOHYt26dQCAlStXYsWKFdiyZQuys7Pxi1/8IiI1fvHFF/jtb3+LESNGYO7cuQDcfZKlpaX4t3/7Nzz11FPYtGkT0tPTvbWHW1ZWFkpKSlBUVARZljF//nyMGTMmIrV09uqrr8LpdHp/gwaABQsW4IEHHsA999wDRVFwxx13YNasWRGpb+rUqThw4ADmzp0LTdOwcOFC3HTTTVF3Lb/++mvvbeQAkJubGzXX0MNqteL555/HkiVL4HQ6MWXKFHz3u98FAGzYsAFPPvkk2tvbcf3116OoqCjs9b3zzjtobGzEr3/9a/z6178GAHznO9/BI488gjVr1uChhx6CLMsYP368t9sw3AJ9Xv1d20D4ZEMiIgrZZdGdRURExmCIEBFRyBgiREQUMoYIERGFjCFCREQhY4gQEVHIGCJERBQyhggREYXsfwHQxhwD5dtFFgAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"mean_sampled = []\n",
"for _ in range(100): # repeat code in block below 1000 times\n",
" sample = df_sub_F['age'].sample(100) # randomly sample 100 observations\n",
" mean_sampled.append(sample.mean()) # append mean of subsample\n",
"pd.Series(mean_sampled).plot(kind='density') # plot the sampling distribution a statistic\n",
"df_sub_M['age'].plot(kind='density') # plot the original data distribution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this figure, we plotted the distribution of the sample statistic (blue) and the distribution of the original data (orange). Notice how different these distributions look: the sampling distribution of a statistic is narrower and centred around the mean. The phenomenon is often referred to as the **central limit theorem**: the sampling distribution of the mean will converge to a normal distribution, which is a bell-shaped curve, with most of the values within two standard deviations from the mean. \n",
"\n",
"We will not provide a formal prove, but simuation using the bootstrap technique should suffice here. Below we take 1000 samples of size=1000 and plot the mean (red), one standard deviation (blue), and two standard deviations (red) to demonstrate that the sampling distribution of the mean is close to a normal distribution."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD7CAYAAABjVUMJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAftElEQVR4nO3dfXAU9f0H8Pcld4mJOQTSO5pRpJYyPpQWLZaHWoi2TUhJTqA6MUjJ+FCMoyHASCFgUiVWZBhsSgyoZfABgQ4pQggxJHZwoNWAyEmROBEYJaCEknBIHpDc4/f3R353vYTLZe9h9+7Y92smM/ewu9/3bXbvc7ff7+1qhBACRESkOnGRDkBERJHBAkBEpFIsAEREKsUCQESkUiwAREQqxQJARKRSLABERCqljXSAQH377WW4XMr/dCE1NQUWS7fi7QZDStbU8WNhMTcplMi/QNbt+PEpMJsj938IZjsYP34sAMAcgfUd6nar5PoO5z6mxPYdK+8JcXEaDBt2vc/nYq4AuFwiIgXA3XasGDTr6dNR9XqkZjl9OvL/h0DbP336dFDzhUso7Sq9vsPWlkLbd6S3xVDxEBARkUqxABARqRQLABGRSrEAEBGpFAsAEZFKsQAQEakUCwARkUrF3O8AiKKRfkgSrkvs3Z16rA50dV6JcCKiwbEAEIXBdYlamJ7ZBQDY/fIMdEU4D5EUPARERKRSLABERCrFAkBEpFIsAEREKsVOYKIgeY/8CeeyOIqIlMICQBSk/iN/wrksjiIiJchaANauXYuGhgZoNBo8+OCDePTRR7Fs2TKYzWYkJSUBAAoLC5GRkSFnDCIi8kG2AnDo0CEcPHgQNTU1cDgcmD59OtLT09HU1ITNmzfDaDTK1TQREUkgWyfwhAkTsGnTJmi1WlgsFjidTiQmJqK1tRWlpaUwmUyoqKiAy+WSKwIREfkh6yEgnU6HiooKvPHGG8jKyoLT6cSkSZNQVlaG5ORkFBQUYPv27cjNzZW8zNTUFBkT+2cw6CPWdqCkZI2m1xNIlkjnDnbdRuo1hrosJdd3NL3uaGlDTrJ3AhcVFWHevHl48sknceDAAaxbt87z3Ny5c1FdXR1QAbBYuiNyHU6DQY/29tjompOS1QBEzesJbN1G9v/gndXfzu8ro7/c/ZcVrtcY+nar3PoO5z6mxPYdK+8JcXGaAT84y3YI6Msvv0RzczMAICkpCZmZmairq0NDQ4NnGiEEtFoORCIiigTZCsA333yDkpIS2Gw22Gw27N27Fz//+c+xcuVKdHR0wG63Y9u2bRwBREQUIbJ9/E5PT8fRo0cxc+ZMxMfHIzMzE4WFhRg2bBhmz54Nh8OBzMxM5OTkyBWBiIj8kPX4S1FREYqKivo8NmfOHMyZM0fOZomISAKeC4iISKVYAIiIVIoFgIhIpVgAiIhUigWAiEilWACIiFSKBYCISKVYAIiIVIoFgIhIpVgAiIhUigWAiEileC5mokHohyThusTeXcVmd0Y4DVH4sAAQDeK6RC1Mz+wCAOx+eUaE0xCFDw8BERGpFAsAEZFKsQAQEakUCwARkUqxABARqRQLABGRSslaANauXYvp06cjOzsbb775JgCgsbERJpMJmZmZKC8vl7N5IiLyQ7bfARw6dAgHDx5ETU0NHA4Hpk+fjsmTJ2P58uV45513kJaWhoKCAuzfvx/p6elyxSAiogHI9g1gwoQJ2LRpE7RaLSwWC5xOJzo7OzFq1CiMHDkSWq0WJpMJ9fX1ckUgIiI/ZP0lsE6nQ0VFBd544w1kZWWhra0NBoPB87zRaMT58+cDWmZqakq4Y0pmMOgj1nagpGSNptcTSJZI5w523UbqNYa6LCXXdzS97mhpQ06ynwqiqKgI8+bNw5NPPomWlparntdoNAEtz2LphsslwpROOoNBj/b2LsXbDYaUrAYgal5PYOtW+f9D/53c3b6/nd9XRn+5B2ojVKFvt8qt73DuY0ps37HynhAXpxnwg7Nsh4C+/PJLNDc3AwCSkpKQmZmJjz/+GBcuXPBM09bWBqPRKFcEIiLyQ7YC8M0336CkpAQ2mw02mw179+5FXl4eTp06hdOnT8PpdKK2thZTp06VKwIREfkh2yGg9PR0HD16FDNnzkR8fDwyMzORnZ2N4cOHY/78+bBarUhPT0dWVpZcEYiIyA9Z+wCKiopQVFTU57HJkyejpqZGzmaJiEgCXg+AKEK8LzRDFAk8FQRRhLgvNOO+2AyR0lgAiIhUigWAiEilWACIiFSKBYCISKVYAIhkpB+S5PM2UTRgASCSkfcwTw75pGjDAkBEpFIsAEREKsUCQESkUjwoSeQDT9NAasBvAEQ+8DQNpAYsAEREKsUCQESkUiwAREQqxQJARKRSHOZA1yzvkTw9Vge6Oq9EOJE0NrsTBoMeQN/csfp6KHqxANA1yz2SBwB2vzwDXRHOI1WCLt5n7lh9PRS9ZC0AlZWV2LNnD4Dei8QvWbIEy5Ytg9lsRlJS74mxCgsLkZGRIWcMIiLyQbYC0NjYiA8//BA7d+6ERqPBH/7wB/zzn/9EU1MTNm/eDKPRKFfTREQkgWydwAaDAcXFxUhISIBOp8Po0aPR2tqK1tZWlJaWwmQyoaKiAi6XS64IRETkh2zfAMaMGeO53dLSgrq6OmzduhWHDh1CWVkZkpOTUVBQgO3btyM3N1fyclNTU+SIK4m7Yy4WSMkaTa8nkCzB5h5sPpvdiQRdfFjaH2iaYLIPtizv3L5eQ6j/ZyW3k3C2pUTuaNqHgiF7J/DJkydRUFCApUuX4oc//CHWrVvneW7u3Lmorq4OqABYLN1wuYQcUf0yGPRob4+NbjcpWQ1A1LyewNat9Gn775yDrhODvk8n60Dcy/G38w80jXcGqW8egy2rf+7+bYT2f1Zuuw/nPqbE9h0r7wlxcZoBPzjL+jsAs9mMRx55BM888wxmzZqF48ePo6GhwfO8EAJaLQciERFFgmwF4Ny5c3j66aexZs0aZGdnA+h9w1+5ciU6Ojpgt9uxbds2jgAiIooQ2T5+b9y4EVarFatWrfI8lpeXhyeeeAKzZ8+Gw+FAZmYmcnJy5IpARER+SCoA77zzDmbNmoWUFOkdsCUlJSgpKfH53Jw5cyQvh4iI5CHpENCJEycwbdo0PPvsszh27JjcmYgoDPRDen9saTDoPbeJvEkqAC+88AIaGhowduxYrFixAg888AC2b98Oq9Uqdz4iCpL3RW14dTPyRXIncEpKCrKyspCTk4NLly5h69atyMrKQn19vZz5iIhIJpI+FjQ2NqKqqgoHDhzAtGnTsG7dOtx22204c+YMHn74YWRlZcmdk4iIwkxSASgrK8PDDz+MF154AXr9/36McvPNNwf0Iy4iIooekg4B1dTUYOjQodDr9Whvb8dbb73lOYdPUVGRrAGJiEgekjuB9+3b1ztDXBzMZjNWrlwpZy6imOW+oIuc54kZqA0l2qZrh6RDQEeOHEFtbS0AIDU1FWvXrsWMGQOfK4VIzfpf0EXJNpRom64dkr4B2O122Gw2z32HwyFbICIiUoakbwD33nsvHn/8ccyYMQMajQa1tbVIT0+XOxsREclIUgFYsmQJtmzZgr1790Kr1SIjIwN5eXlyZyMiIhlJKgDx8fHIz89Hfn6+3HmIiEghkgpAXV0d1qxZg46ODgjxv4uxfPrpp7IFIyIieUkqAGvXrkVxcTHuuOMOaDQauTMREZECJBWAIUOGIDMzU+4sRESkIEnDQMeNG4f9+/fLnYWIiBQk6RvA/v37sXnzZuh0Ouh0OgghoNFo2AdARBTDJBWAt956S+YYRMrRD0nynB+/x+pAV+eVCCciigxJh4BuvPFGHDt2DFVVVRg+fDiOHDmCG2+8Ue5sRLLghVKIekkqAH/729/w97//HfX19ejp6UFlZSXWrVs36HyVlZXIzs5GdnY2Vq9eDaD32gImkwmZmZkoLy8PLT0REQVNUgF47733sGHDBiQlJWHYsGGoqqrynBxuII2Njfjwww+xc+dOVFdX4/PPP0dtbS2WL1+O9evXo66uDk1NTexcJiKKEEkFQKvVIiEhwXN/yJAh0Gr9f3U2GAwoLi5GQkICdDodRo8ejZaWFowaNQojR46EVquFyWTiJSWJiCJE0gHQtLQ07Nu3DxqNBjabDRs3bhy0D2DMmDGe2y0tLairq8PcuXNhMBg8jxuNRpw/fz6gwKmpKQFNH06xdI51KVmj6fUEkiWY3O7z5IdjecGuN38Z5GgvXMtRcjsJZ1tK5I6mfSgYkgpAaWkplixZguPHj+POO+/EuHHjsGbNGkkNnDx5EgUFBVi6dCm0Wi1OnTrV5/lAf1lssXTD5RKDTxhmBoMe7e1dircbDClZDUDUvJ7A1q30ab13Tn/nyXcvT+rOHOj0vjL4yhHu9gZaTmCU2+7DuY8psX3HyntCXJxmwA/OkgrAiBEj8Pbbb+PKlStwOp1ISZH2KdxsNqOoqAjLly9HdnY2Dh06hAsXLnieb2trg9FolLQsIiIKL0kF4M033/T5+KOPPjrgPOfOncPTTz+N8vJyTJ48GUDvL4pPnTqF06dP46abbkJtbS0eeOCBIGITEVGoJBWAEydOeG7bbDaYzWZMnDjR7zwbN26E1WrFqlWrPI/l5eVh1apVmD9/PqxWK9LT05GVlRVkdCIiCoWkAvDSSy/1uX/x4kUsWbLE7zwlJSUoKSnx+VxNTY3EeEREJJegfgY5fPhwnD17NtxZiKJesCN5iKJRwH0AQgg0NTUhNTVVtlBE0crfaCKiWBNwHwDQ+7uAwQ4BERFRdAuqD4CIiGKfpAIwd+5cvz/Y2rRpU9gCERGRMiQVgLFjx+LLL79Ebm4udDoddu3aBYfDgezsbLnzERGRTCQVgE8//RRbt25FfHw8AGDKlCnIzc3FtGnTZA1HpHYcdURyknQ20IsXL8Jms3nuX758GT09PbKFIqJe7lFH3ucQIgoXSd8AcnJykJubi4yMDAghsGfPHuTn58udjYiIZCSpACxYsAB33HEHDh48iMTERJSVlWHChAlyZyMiIhlJOgQE9J4RdMyYMVi4cCF0Op2cmYiISAGSvgG8++67eOONN2C1WpGRkYGnnnoKixYtQm5urtz5SOX0Q5I8F27vsTokTdPVeUWxfLGI64vcJH0D2Lx5M7Zt24aUlBSkpqZix44dePvtt+XORoTrErWeTlD3m1Yw09D/cH2Rm6QCEBcX1+ciMGlpaZ4hoUREFJskFYChQ4eiubnZ82vgmpoa3HDDDbIGIyIieUn6/rd8+XIsWLAAZ86cwS9/+UskJiZi/fr1cmcjIiIZSSoAPT092LVrF1paWuB0OnHLLbdwJBARUYyTVAAWL16MPXv2YPTo0XLnoWtcuEageC/HG0+dEBjv9cURQeojqQDceuut2L17N8aPH4/k5GTP40OHDpUtGF2b3CNQgN4LqnSFaTluvGBLYPqvr2D/HxSbJBWAvXv3or6+vs9jGo0Gzc3Ng87b3d2NvLw8vPbaa7jpppuwbNkymM1mJCUlAQAKCwuRkZERRHQiIgqFpAJw7NixoBZ+9OhRlJSUoKWlxfNYU1MTNm/eDKPRGNQyiYgoPPwOAy0tLfXcvnjxYsALr6qqwnPPPed5s//uu+/Q2tqK0tJSmEwmVFRUwOVyBbxcIiIKnd9vAE1NTZ7bjz/+OHbu3BnQwl988cU+9y0WCyZNmoSysjIkJyejoKAA27dvD+iUEqmpKYNPJJNY6lyUkjUaXo87QyBZwplbTZ3GgWwTNrsTCbreH3veMDTZc9v7ce/bSmaMxLIi2Yac/BYAIYTP28EaOXIk1q1b57k/d+5cVFdXB1QALJZuuFyhZwmUwaBHe3tsdJFJyWoAIvJ6+u8w7e1dfvP62sHc84SDmjqN3evY37rznqZ3vcxAgi7ex+O96yuc21A49zEltu9YeU+Ii9MM+MFZ8tlA/V0TWKrjx4+joaHBc18IAa2W5yIhIooEvwXA5XKho6MDly5dgtPp9Nx2/wVKCIGVK1eio6MDdrsd27Zt4wggIqII8fvx+8SJE5g0aZLn8M/EiRM9z0kdBurttttuwxNPPIHZs2fD4XAgMzMTOTk5QcQmIqJQ+S0AX3zxRVga+eCDDzy358yZgzlz5oRluUREFDwegKeoM9BpHkheA42IkjJSiheZiU3cyyjqDHSaB5LXQCOipIyUCtcpPkhZkkcBERHRtYUFgIhIpVgAiIhUigWAiEil2AlMsuOonshT8pxHvMhM7OBeSbLjqJ7IU/KcR7zITOzgISAiIpViASAiUikWACIilWIfAAWEP/mnYHHbiT4sABQQ/uSfgsVtJ/rwEBARkUqxABARqRQLABGRSrEAEBGpFAsAxQyb3QkAip3SgOhaxwJAMSNBFw8AnpEkRBQaWQtAd3c3cnJy8M033wAAGhsbYTKZkJmZifLycjmbJiKiQchWAI4ePYrZs2ejpaUFANDT04Ply5dj/fr1qKurQ1NTE/bv3y9X80RENAjZCkBVVRWee+45GI1GAMBnn32GUaNGYeTIkdBqtTCZTKivr5ereSIiGoRsvwR+8cUX+9xva2uDwWDw3DcajTh//nzAy01NTQk5W7BiqfNRSlapr8dmd3qOvwe7jMEy3DA0ecA2KHqFsp35ejyc+5gS+2ssvSf4otipIIQQVz2m0WgCXo7F0g2X6+plyc1g0KO9PTZ+vC4lqwGQ/HoMBv2A55KXsgx/O0l7excMBr2i56un8HH//wf7H/uapv+2E859LJDtO+g2YuQ9IS5OM+AHZ8VGAY0YMQIXLlzw3G9ra/McHiIiIuUpVgDGjRuHU6dO4fTp03A6naitrcXUqVOVap6IiPpR7BBQYmIiVq1ahfnz58NqtSI9PR1ZWVlKNU9ERP3IXgA++OADz+3JkyejpqZG7iaJiEgCXg+AiIJiszvDMgrG+0IxBoOeF4tREAsAEQUlXCO3vC8U415W9I+tuTbwXEBERCrFAkBEpFIsAEREKsUCQESkUuwEpqB5jwKx2pxITOg9lw9HcZBbuEYKkTxYACho/UeBeN/mKA4CwjdSiOTBQ0BERCrFAkBEpFIsAEREKsU+ACJSnL/OYe/nOKBAXiwARKQ4f53D/Z/jgAL58BAQEZFKsQAQEakUCwARkUqxABARqRQ7gYnomuR9oRmOJvKNBYCIrkneF5rhaCLfIlIA8vPzYbFYoNX2Nl9WVoZx48ZFIgoRkWopXgCEEPjqq6+wb98+TwEgIiLlKd4J/NVXX0Gj0WDevHm4//77sXnzZqUjEBERIvANoLOzE5MnT8bzzz+Pnp4e5Ofn45ZbbsE999wjaf7U1BSZEw4s0uc1t9mdSNDFX3XbFylZ/U0z2PL9kXoOeJ4rngbjvY0MtP1b7U4kStgv5NjWYn37VbwA3HXXXbjrrrsAAMnJyXjwwQexf/9+yQXAYumGyyXkjOiTwaBHe3tku5EMBn2fTq2B8kjJagD8TtO/rUB4/5Tf3/w8VzwNpv824t5m+2+f/acx+FhWuPffaHhPkCIuTjPgB2fFDwEdPnwYBw4c8NwXQrAvgIgoAhQvAF1dXVi9ejWsViu6u7uxc+dOZGRkKB2DiEj1FP/ofd999+Ho0aOYOXMmXC4XHn74Yc8hISIiUk5Ejr0sXLgQCxcujETTRET0/3jwnYhigpRRY3KMLLuWTynBk8ERUUxwjwjyHmE20DTh5D6lhOmZXZ5CcK1gASAiUikWACIilWIBICJSKRYAIiKVurZ6NK4R3qMOrDYnEhOCOyfPQKMX9EOSAPT+lP1aG9VARNLxG0AU8h51kJgw+MgHKcvxHr3gvn0tjmogIulYAIiIVIoFgIhIpVgAiIhUigeAI0iOn5h7d/BK4f3TeXYIkxr4GxwRyP7ove94D9aIpf2IBSCC3J20QO+FLMJxaYn+yxxM/wtuRP/lLYhCM9B+F+j+2H/ficX9iIeAiIhUigWAiEilWACIiFSKBYCISKVU0wms9EUdQhlRIEWgo3fkuFAGUawIZf8Kpa2BRhlJGTXkPb2/6UKhmgIgx4ibcLbXf0TBYAIdvRPo8omuJaHuX6G0NdAoo8H2X+/p/U0XiogcAtq9ezemT5+OjIwMbNmyJRIRiIhUT/FvAOfPn0d5eTl27NiBhIQE5OXlYeLEifjRj36kdBQiIlVTvAA0NjZi0qRJGDp0KABg2rRpqK+vR2FhoaT54+I0QbdtHJYU0nICnUdKe97ThHLbe/lSpseoUZ77oSxHUlthXNaoUfCZO5pvjxo1yufj0ZLP3213dCXbDcuy/j+4HPuOHNNLeX/wN50//ubRCCFEwEsMweuvv47vvvsOixYtAgD84x//wGeffYYXXnhByRhERKqneB+Ar3qj0QT/qZ6IiIKjeAEYMWIELly44Lnf1tYGo9GodAwiItVTvAD84he/wIEDB3Dx4kVcuXIF77//PqZOnap0DCIi1VO8E3jEiBFYtGgR8vPzYbfb8eCDD+KnP/2p0jGIiFRP8U5gIiKKDjwXEBGRSrEAEBGpFAsAEZFKsQAQEamUas4G6k9lZSX27NkDAEhPT8eSJUuwdetWbNmyBUIIz2P9f7DW2tqKP/7xj7BYLLjllluwZs0aXH/99VGZtbq6GmvWrEFqaioA4N577/X8GlvJrG5btmxBfX093nnnnavm6+zsxOLFi/H1119j+PDh+Otf/wqDwSBr1lDyfvLJJygsLMT3v/99AMAdd9yBl156SfGsy5Ytg9lsRlJS7+kDCgsLkZGR0We+5uZmlJSUoLu7G3fffTdWrFgBrVb+t4Fg81ZWVuLdd9/FkCFDAAC5ubmYM2eO4lmPHDmCl156CZcvX8att96KVatWISEhoc98kXg/CJlQuY8++kg89NBDwmq1CpvNJvLz88Wbb74pMjIyxOXLl4XD4RAPPfSQ+Pe//33VvE888YSora0VQghRWVkpVq9eHbVZy8rKxO7du2XNN1jW999/XwghxMmTJ8WUKVPE73//e5/zrlixQrz++utCCCF27twpFixYENV5N27cKF577TXZMw6WNScnR5w/f97vvNnZ2eLIkSNCCCGWLVsmtmzZEtV5CwoKxKeffip7RjdfWXfs2CHuuece0dzcLIQQYtGiRT7Xm9LvB+Gg+kNABoMBxcXFSEhIgE6nw+jRo6HRaPDee+8hOTkZnZ2d6O7u9nwCcbPb7fjkk08wbdo0AMDvfvc71NfXR2VWADh27Biqq6tx//33Y/Hixejo6FA8a2trK2w2G/70pz9hwYIFA867b98+mEwmAEBOTg7+9a9/wW63R23eY8eO4aOPPsLMmTPx5JNP4ty5cxHJ2traitLSUphMJlRUVMDlcvWZ7+zZs+jp6cGdd94JQJltNpS8ANDU1IQNGzbAZDKhrKwMVqtV8axnz57FnXfeidtuuw0AUFJSctU3lUi8H4SD6gvAmDFjPDtES0sL6urqkJ6eDp1Oh6qqKvzmN7+BwWDw/PPdvv32W6SkpHi+PhsMBpw/fz4qs7rzzZ8/H7t27UJaWhrKysoikvXll1/GAw88gJtuumnAedva2jyHfLRaLVJSUnDx4sWozavX65Gfn4/q6mqkp6fLfmjNV9YpU6Zg0qRJWLlyJaqqqnD48GFs3769z3ze6xVQZpsNJe/ly5dx++23Y+nSpdi5cyc6Ozuxfv16xbMmJCQgOTkZTz/9NEwmE1555ZWrPmRF4v0gLCL9FSRanDhxQtx3331ix44dfR632+1i8eLF4uWXX+7z+H//+18xZcqUPtONHTs2KrP2d+nSJXH33XfLGdHDO+uHH34o5s+fL4QQ4uDBgwMeUvnxj38s7Ha75/6UKVNEW1tb1Obtb/z48aKzs1POmEKIgbcDIYR4//33xVNPPdXnMbPZLPLy8jz3W1paxLRp02TP6RZo3v4+//xzMWPGDLni9eGddf369WLy5MnizJkzwuFwiKVLl4qKioo+00fy/SAUqv8GAABmsxmPPPIInnnmGcyaNQvnzp2D2WwG0PsJNDs7G8ePH+8zz/Dhw9Hd3Q2n0wkAaG9vV+SkdsFk7erqwltvveW5L4RQpOOvf9ba2lqcPHkSM2bMQElJCZqamrBw4cKr5jMajZ4TBjocDnR3d3uuHxFteV0uF1599VXPduAm9/rtn/X48eNoaGjwPO/rf9z/RIxKbbPB5m1tbe3zrSBS2+33vvc9jBs3DiNHjkR8fDx++9vf4rPPPuszT6TeD0IW2foTea2trWLixImisbHR89jx48fFfffdJzo6OoTL5RLFxcWeTklv8+bNEzU1NUIIIdavXy+ef/75qMzqcDjEPffcI/7zn/8IIYR45ZVXRGlpqeJZvfn7RP3888+LV199VQghxK5du8S8efNky+kWSt5Zs2aJ9957TwjR22n92GOPyZZTCN9Zm5ubxdSpU8WlS5eEzWYTjz32mM9O/+zsbHH48GEhhBDPPvus2LBhg6xZQ8lrsVjEhAkTxJkzZ4TL5RLLli2TvbPdV9bW1lYxZcoU0draKoQQ4rnnnhPl5eVXzav0+0E4qP5cQH/+85/x7rvv4uabb/Y8lpeXB41Gg02bNiE+Ph533303li9fDp1Oh2effRa/+tWv8Otf/xpnz55FcXExLBYL0tLS8Je//AU33HBDVGY9fPgwXnzxRfT09OAHP/gBVq9eDb1er3jW2bNnAwA+/vhjVFZWeoZVrl27FkajEbNnz8alS5dQXFyMr7/+Gnq9HmvWrPF7DD7SeU+ePInS0lJ0dXVh+PDhWL16NdLS0hTP6nK5sGXLFjgcDmRmZmLx4sUAgHnz5qGoqAg/+clP8MUXX6CkpASXL1/2DFftP5wxmvI2NDTglVdegd1ux89+9jOsWLFC1rwDZU1LS0N5eTmsVituv/12rFy5EklJSRF9PwgH1RcAIiK1Yh8AEZFKsQAQEakUCwARkUqxABARqRQLABGRSrEAEBGpFAsAEZFKsQAQEanU/wFqLIdKbeN5AwAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"mean_sampled = []\n",
"for _ in range(1000): # repeat code in block below 1000 times\n",
" sample = df_sub_F['age'].sample(1000) # randomly sample 1000 observations\n",
" mean_sampled.append(sample.mean()) # append mean of subsample\n",
"ax = pd.Series(mean_sampled).plot(kind='hist',bins=100) # plot the sampling distribution a statistic\n",
"\n",
"sdist_mean = pd.Series(mean_sampled).mean()\n",
"sdist_std = pd.Series(mean_sampled).std()\n",
"\n",
"ax.axvline(x = sdist_mean, color='black', lw=2) # plot mean of distribution of the sampling statistic\n",
"ax.axvline(x = sdist_mean + sdist_std, color='blue', lw=1) # plot line one standard deviation from the mean\n",
"ax.axvline(x = sdist_mean - sdist_std, color='blue', lw=1) # plot line one standard deviation from the mean\n",
"ax.axvline(x = sdist_mean + (2*sdist_std), color='red', lw=1) # plot line two standard deviations from the mean\n",
"ax.axvline(x = sdist_mean - (2*sdist_std), color='red', lw=1) # plot line two standard deviations from the mean"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After computing the sampling distribution of the mean for women in Whitechapel, we can easily obtain the 95% confidence interval, which implies that 95% all the possible means are locature withing this interval. We can apply `.quantile()` after converting the `mean_sampled` to a `pd.Series` object (you can not apply `.quantile()` a normal Python list. \n"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.025 23.610825\n",
"0.975 25.713175\n",
"dtype: float64"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pd.Series(mean_sampled).quantile([0.025,0.975])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The mean age of women is likely to vary between 23.6 and 25.7. We'd interpret this a range of values that are produced by random chance. The probability of observing more extreme values is lower then 0.05. "
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD7CAYAAABjVUMJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAe3UlEQVR4nO3de1BU5/0G8GdhFwIFY6S7lkmMTa2TxNqaNNZLrZK0BamwUatDUCuTSw1OgqgTq6jQRNIYxzGlEjRJM+Zi1I7UKCIi2DGjrUFjJNZIhqiTiCZiBdfIxche2Pf3Bz82C7LLYdn37K7n+cw4s5dzeTies9/d877nPTohhAAREWlOWKADEBFRYLAAEBFpFAsAEZFGsQAQEWkUCwARkUaxABARaRQLABGRRukDHaCvvvnmOpxO9S9diIuLgcXSqvp6fRFKWYHQyhtKWYHQyhtKWQE5eR96KAYAUF3tv+WGhelwxx3f6/G9kCsATqcISAHoXHeoCKWsQGjlDaWsQGjlDaWsgP/znj8vZ7me8BQQEZFGsQAQEWkUCwARkUaxABARaRQLABGRRrEAEBFpFAsAEZFGhdx1AETBKHZAFG6L7Dic2qwOtDTfCHAiot6xABD5wW2Repif2w0A2PPKVLQEOA+REjwFRESkUSwAREQaxQJARKRRLABERBrFRmAiH7n3/PHnstiLiNTCAkDko+49f/y5LPYiIjVILQDr169HZWUldDodZs6ciSeeeALLly9HdXU1oqKiAABZWVlITEyUGYOIiHogrQAcO3YMR48eRWlpKRwOB6ZMmYKEhATU1NRgy5YtMJlMslZNREQKSGsEHjNmDDZv3gy9Xg+LxYL29nZERkaivr4eeXl5MJvNKCwshNPplBWBiIi8kHoKyGAwoLCwEG+99RaSk5PR3t6OcePGIT8/H9HR0cjMzMSOHTuQlpameJlxcTESE3tnNMYGbN19FUpZgdDKqyRrf/8ef26PW23bBhNZedXaDtIbgbOzszFv3jzMnz8fR44cwYYNG1zvzZ07FyUlJX0qABZLa0DuG2o0xqKxMTSa5kIpKxBaed2zejtI+/r3dF+Wv7ZHqG7bUCAnb8d+4M/lhoXpPH5xlnYK6IsvvkBtbS0AICoqCklJSSgvL0dlZaVrGiEE9Hp2RCIiCgRpBeDrr79Gbm4ubDYbbDYbDhw4gF/84hdYvXo1mpqaYLfbsX37dvYAIiIKEGlfvxMSEnDy5ElMmzYN4eHhSEpKQlZWFu644w7MmjULDocDSUlJSE1NlRWBiIi8kHr+JTs7G9nZ2V1emzNnDubMmSNztUREpADHAiIi0igWACIijWIBICLSKBYAIiKNYgEgItIoFgAiIo1iASAi0igWACIijWIBICLSKBYAIiKNYgEgItIojsVM1IvYAVG4LbLjULHZ2wOchsh/WACIenFbpB7m53YDAPa8MjXAaYj8h6eAiIg0igWAiEijWACIiDSKBYCISKNYAIiINIoFgIhIo6QWgPXr12PKlClISUnB22+/DQCoqqqC2WxGUlISCgoKZK6eiIi8kHYdwLFjx3D06FGUlpbC4XBgypQpGD9+PFasWIH33nsP8fHxyMzMxKFDh5CQkCArBhEReSDtF8CYMWOwefNm6PV6WCwWtLe3o7m5GUOHDsWQIUOg1+thNptRUVEhKwIREXkh9Upgg8GAwsJCvPXWW0hOTkZDQwOMRqPrfZPJhMuXL/dpmXFxMf6OqZjRGBuwdfdVKGUFQiuvkqz9/Xv8uT1utW0bTGTlVWs7SB8KIjs7G/PmzcP8+fNRV1d30/s6na5Py7NYWuF0Cj+lU85ojEVjY4vq6/VFKGUFgj9v94OxM6u3g7Svf4+ndfRXsG9bd6GUFZCVt2M/8Odyw8J0Hr84SzsF9MUXX6C2thYAEBUVhaSkJHz00Ue4cuWKa5qGhgaYTCZZEYiIyAtpBeDrr79Gbm4ubDYbbDYbDhw4gPT0dJw7dw7nz59He3s7ysrKMGnSJFkRiIjIC2mngBISEnDy5ElMmzYN4eHhSEpKQkpKCgYNGoQFCxbAarUiISEBycnJsiIQEZEXUtsAsrOzkZ2d3eW18ePHo7S0VOZqiYhIAd4PgChA3G80QxQIHAqCKEA6bzTTebMZIrWxABARaRQLABGRRrEAEBFpFAsAEZFGsQsCkUTuPX3arA60NN8IcCKi77AAEEnU2dMHAPa8MhWhM9INaQFPARERaRQLABGRRrEAEBFpFNsAiHrAYRpIC/gLgKgHHKaBtIAFgIhIo1gAiIg0igWAiEijWACIiDSK3RzolhWqwzDY7O0wGmMBdM0dqn8PBS8WALplheowDBGG8B5zh+rfQ8FLagEoKirCvn37AHTcJH7p0qVYvnw5qqurERUVBQDIyspCYmKizBhERNQDaQWgqqoKhw8fxq5du6DT6fDHP/4R//rXv1BTU4MtW7bAZDLJWjURESkgrRHYaDQiJycHERERMBgMGDZsGOrr61FfX4+8vDyYzWYUFhbC6XTKikBERF5I+wUwfPhw1+O6ujqUl5dj27ZtOHbsGPLz8xEdHY3MzEzs2LEDaWlpipcbFxcjI64inQ1zoSCUsgLq5O1tHTZ7OyIM4f1ejrdpfPk7e1uWe+6e/oZQ2hdCKSsgL69a20F6I/DZs2eRmZmJZcuW4Uc/+hE2bNjgem/u3LkoKSnpUwGwWFrhdAoZUb0yGmPR2BgazW6hlBWQl7f7QdTbOozG2C6NrJ50LsfbQeppGvcMSg/y3pbVPXf3dYTKvhBKWQFZeTv+j/253LAwnccvzlKvA6iursbjjz+O5557DtOnT8fp06dRWVnpel8IAb2eHZGIiAJBWgG4dOkSnn32Waxbtw4pKSkAOj7wV69ejaamJtjtdmzfvp09gIiIAkTa1+9NmzbBarVizZo1rtfS09Px9NNPY9asWXA4HEhKSkJqaqqsCERE5IWiAvDee+9h+vTpiIlR3gCbm5uL3NzcHt+bM2eO4uUQEZEcik4BnTlzBpMnT8bKlStx6tQp2ZmIyA9iB3RcbGk0xroeE7lTVABefPFFVFZWYuTIkVi1ahVmzJiBHTt2wGq1ys5HRD5yv6kN725GPVHcCBwTE4Pk5GSkpqbi2rVr2LZtG5KTk1FRUSEzHxERSaLoa0FVVRWKi4tx5MgRTJ48GRs2bMB9992HCxcuYPbs2UhOTpadk4iI/ExRAcjPz8fs2bPx4osvIjb2u4tR7r777j5dxEVERMFD0Smg0tJSDBw4ELGxsWhsbMQ777zjGsMnOztbakAiIpJDcSPwwYMHO2YIC0N1dTVWr14tMxdRyOq8oYvM8Vw8rUONddOtQ9EpoBMnTqCsrAwAEBcXh/Xr12PqVM9jpRBpWfcbuqi5DjXWTbcORb8A7HY7bDab67nD4ZAWiIiI1KHoF8DDDz+Mp556ClOnToVOp0NZWRkSEhJkZyMiIokUFYClS5di69atOHDgAPR6PRITE5Geni47GxERSaSoAISHhyMjIwMZGRmy8xARkUoUFYDy8nKsW7cOTU1NEOK7m7F88skn0oIREZFcigrA+vXrkZOTgxEjRkCn08nOREREKlBUAAYMGICkpCTZWYiISEWKuoGOGjUKhw4dkp2FiIhUpOgXwKFDh7BlyxYYDAYYDAYIIaDT6dgGQEQUwhQVgHfeeUdyDCL1xA6Ico2P32Z1oKX5RoATEQWGolNAd955J06dOoXi4mIMGjQIJ06cwJ133ik7G5EUvFEKUQdFBeDvf/87/vGPf6CiogJtbW0oKirChg0bep2vqKgIKSkpSElJwdq1awF03FvAbDYjKSkJBQUF/UtPREQ+U1QA9u7dizfffBNRUVG44447UFxc7BoczpOqqiocPnwYu3btQklJCT777DOUlZVhxYoV2LhxI8rLy1FTU8PGZSKiAFFUAPR6PSIiIlzPBwwYAL3e+09no9GInJwcREREwGAwYNiwYairq8PQoUMxZMgQ6PV6mM1m3lKSiChAFJ0AjY+Px8GDB6HT6WCz2bBp06Ze2wCGDx/uelxXV4fy8nLMnTsXRqPR9brJZMLly5f7FDguLqZP0/tTKI2xHkpZAfl5O8fJ98e6fc3qLYOM9clajkyhkNGdrLxqbQdFBSAvLw9Lly7F6dOn8cADD2DUqFFYt26dohWcPXsWmZmZWLZsGfR6Pc6dO9fl/b5eWWyxtMLpFL1P6GdGYywaG1tUX68vQikrIC+v+0HkbZz8znUrPej6On1PGXrK4e/1eVpOsOJ+CwAd/8f+XG5YmM7jF2dFBWDw4MF49913cePGDbS3tyMmRtm38OrqamRnZ2PFihVISUnBsWPHcOXKFdf7DQ0NMJlMipZFRET+pagAvP322z2+/sQTT3ic59KlS3j22WdRUFCA8ePHA+i4ovjcuXM4f/487rrrLpSVlWHGjBk+xCYiov5SVADOnDnjemyz2VBdXY2xY8d6nWfTpk2wWq1Ys2aN67X09HSsWbMGCxYsgNVqRUJCApKTk32MTkRE/aGoALz88stdnl+9ehVLly71Ok9ubi5yc3N7fK+0tFRhPCIiksWnyyAHDRqEixcv+jsLUdDztScPUTDqcxuAEAI1NTWIi4uTFoooWHnrTUQUavrcBgB0XBfQ2ykgIiIKbj61ARARUehTVADmzp3r9YKtzZs3+y0QERGpQ1EBGDlyJL744gukpaXBYDBg9+7dcDgcSElJkZ2PiIgkUVQAPvnkE2zbtg3h4eEAgIkTJyItLQ2TJ0+WGo5I69jriGRSNBro1atXYbPZXM+vX7+OtrY2aaGIqENnryP3MYSI/EXRL4DU1FSkpaUhMTERQgjs27cPGRkZsrMREZFEigrAwoULMWLECBw9ehSRkZHIz8/HmDFjZGcjIiKJFJ0CAjpGBB0+fDgWLVoEg8EgMxMREalA0S+A999/H2+99RasVisSExPxzDPPYPHixUhLS5OdjzQudkCU68btbVYHWppv+DQNfYfbizop+gWwZcsWbN++HTExMYiLi8POnTvx7rvvys5GhNsi9a5G0M4PLV+moe9we1EnRQUgLCysy01g4uPjXV1CiYgoNCkqAAMHDkRtba3rauDS0lLcfvvtUoMREZFcin7/rVixAgsXLsSFCxfwq1/9CpGRkdi4caPsbEREJJGiAtDW1obdu3ejrq4O7e3tuOeee9gTiIgoxCkqAEuWLMG+ffswbNgw2XnoFuevHijuy3HHoRP6xn17sUeQ9igqAPfeey/27NmDhx56CNHR0a7XBw4cKC0Y3Zo6e6AAHTdUafHTcjrxhi190317+fr/QaFJUQE4cOAAKioqurym0+lQW1vb67ytra1IT0/H66+/jrvuugvLly9HdXU1oqKiAABZWVlITEz0IToREfWHogJw6tQpnxZ+8uRJ5Obmoq6uzvVaTU0NtmzZApPJ5NMyiYjIP7x2A83Ly3M9vnr1ap8XXlxcjOeff971Yf/tt9+ivr4eeXl5MJvNKCwshNPp7PNyiYio/7z+AqipqXE9fuqpp7Br164+Lfyll17q8txisWDcuHHIz89HdHQ0MjMzsWPHjj4NKREXF9P7RJKEUuNiqGTtzKk0r7//Li01Giv5OzunsdnbEWEIV/xYzYzBRFZetbaD1wIghOjxsa+GDBmCDRs2uJ7PnTsXJSUlfSoAFksrnM7+Z+krozEWjY2h0UQWzFm779iNjS1e8/Y0fU+v+0pLjcZKtp37NO7bpbfX/SGY99ueyMnb8X/jz+WGhek8fnFWPBqot3sCK3X69GlUVla6ngshoNdzLBIiokDwWgCcTieamppw7do1tLe3ux53/usrIQRWr16NpqYm2O12bN++nT2AiIgCxOvX7zNnzmDcuHGu0z9jx451vae0G6i7++67D08//TRmzZoFh8OBpKQkpKam+hCbiIj6y2sB+Pzzz/2ykg8++MD1eM6cOZgzZ45flktERL7jCXgKOp6GeSC5PPWIUtJTijeZCU08yijoeBrmgeTy1CNKSU8pfw3xQepS3AuIiIhuLSwAREQaxQJARKRRLABERBrFRmCSjr16Ak/NMY94k5nQwaOSpGOvnsBTc8wj3mQmdPAUEBGRRrEAEBFpFAsAEZFGsQ2A+oSX/JOvuO8EHxYA6hNe8k++4r4TfHgKiIhIo1gAiIg0igWAiEijWACIiDSKjcAUMtQczoBIC/gLgEJG5xADnT1JiKh/pBaA1tZWpKam4uuvvwYAVFVVwWw2IykpCQUFBTJXTUREvZBWAE6ePIlZs2ahrq4OANDW1oYVK1Zg48aNKC8vR01NDQ4dOiRr9URE1AtpBaC4uBjPP/88TCYTAODTTz/F0KFDMWTIEOj1epjNZlRUVMhaPRER9UJaI/BLL73U5XlDQwOMRqPruclkwuXLl/u83Li4mH5n81UoNUD6M6vN3o4IQ7iU9XTOf/vAaI/roOCl5P/f0zQ9vR5KxxggL69a20G1XkBCiJte0+l0fV6OxdIKp/PmZclmNMaisTE0Ll73d1ajMdbjWPJK1uNtZ25sbIHRGKvqePXkP53//739H/c0Tfd9J5SOMUBW3o5t5M/lhoXpPH5xVq0X0ODBg3HlyhXX84aGBtfpISIiUp9qBWDUqFE4d+4czp8/j/b2dpSVlWHSpElqrZ6IiLpR7RRQZGQk1qxZgwULFsBqtSIhIQHJyclqrZ6IiLqRXgA++OAD1+Px48ejtLRU9iqJiEgBDgVBRD7x19Ac7jeKMRpjebMYFbEAEJFP/NVzy/1GMZ3LCp2+QKGNYwEREWkUCwARkUaxABARaRQLABGRRrERmHzm3gvEamtHZETHWD7sxUGdeBOf4MYCQD7r3gvE/TF7cRDgv55CJAdPARERaRQLABGRRrEAEBFpFNsAiEh13hqH3d9jhwK5WACISHXeGoe7v8cOBfLwFBARkUaxABARaRQLABGRRrEAEBFpFBuBieiW5H6jGfYm6hkLABHdktxvNMPeRD0LSAHIyMiAxWKBXt+x+vz8fIwaNSoQUYiINEv1AiCEwJdffomDBw+6CgAREalP9UbgL7/8EjqdDvPmzcOjjz6KLVu2qB2BiIgQgF8Azc3NGD9+PF544QW0tbUhIyMD99xzDyZMmKBo/ri4GMkJPQv0uOY2ezsiDOE3Pe5Jf7P2tvze5lWyfo4VT71x30c87f9WezsiFRwXMvY1WfuvWseF6gXgwQcfxIMPPggAiI6OxsyZM3Ho0CHFBcBiaYXTKWRG7JHRGIvGxsA2IxmNsV0atTzl8UfW7uvqC/dL+b3Nz7HiqTfd95HO/br7/ulpGnf+Pn7lfCZ0ZPbncsPCdB6/OKt+Cuj48eM4cuSI67kQgm0BREQBoHoBaGlpwdq1a2G1WtHa2opdu3YhMTFR7RhERJqn+lfvRx55BCdPnsS0adPgdDoxe/Zs1ykhIiJST0DOvSxatAiLFi0KxKqJiOj/8eQ7EYUEJb3GZPQsu5WHlGABIKKQoKTXmIyeZbfykBIcDZSISKNYAIiINIoFgIhIo1gAiIg0io3AQci914HV1o7ICN/G5PHUe+FW7tVARMqxAASh7r0OfO3V4Kn3wq3cq4GIlOMpICIijWIBICLSKBYAIiKNYhtAAMlojI0dEAVA+Q0l3C+dZ4MwaYG/Oke4HzvunTVC6ThiAQggGY2x3ZfZm+6XzrNBmG51/uoc0f3YCcXjiKeAiIg0igWAiEijWACIiDSKBYCISKM00wis9vAH/elRoERfe+/IuFEGUajoz/HVn3V56mWkpNeQ+/TepusPzRQAtYc/6G+Pgt70tfeOjBtlEIWK/h5f/VmXp15GvR2/7tN7m64/AnIKaM+ePZgyZQoSExOxdevWQEQgItI81X8BXL58GQUFBdi5cyciIiKQnp6OsWPH4sc//rHaUYiINE31AlBVVYVx48Zh4MCBAIDJkyejoqICWVlZiuYPC9P5vG7THVH9Wk5f51GyPvdp+vPYffnBsBz3x/5c1q30OFhyBNtjtdfX131exvSd0wwd2vV59+3i788tnRBC9HmJ/fDGG2/g22+/xeLFiwEA//znP/Hpp5/ixRdfVDMGEZHmqd4G0FO90el8/1ZPRES+Ub0ADB48GFeuXHE9b2hogMlkUjsGEZHmqV4AfvnLX+LIkSO4evUqbty4gf3792PSpElqxyAi0jzVG4EHDx6MxYsXIyMjA3a7HTNnzsTPfvYztWMQEWme6o3AREQUHDgWEBGRRrEAEBFpFAsAEZFGsQAQEWmUZkYD9aaoqAj79u0DACQkJGDp0qXYtm0btm7dCiGE67XuF6zV19fjT3/6EywWC+655x6sW7cO3/ve94Iya0lJCdatW4e4uDgAwMMPP+y6GlvNrJ22bt2KiooKvPfeezfN19zcjCVLluCrr77CoEGD8Le//Q1Go1Fq1v7k/fjjj5GVlYUf/OAHAIARI0bg5ZdfVj3r8uXLUV1djaiojuEDsrKykJiY2GW+2tpa5ObmorW1FaNHj8aqVaug18v/GPA1b1FREd5//30MGDAAAJCWloY5c+aonvXEiRN4+eWXcf36ddx7771Ys2YNIiIiuswXiM+DfhMa9+GHH4rHHntMWK1WYbPZREZGhnj77bdFYmKiuH79unA4HOKxxx4T//nPf26a9+mnnxZlZWVCCCGKiorE2rVrgzZrfn6+2LNnj9R8vWXdv3+/EEKIs2fPiokTJ4o//OEPPc67atUq8cYbbwghhNi1a5dYuHBhUOfdtGmTeP3116Vn7C1ramqquHz5std5U1JSxIkTJ4QQQixfvlxs3bo1qPNmZmaKTz75RHrGTj1l3blzp5gwYYKora0VQgixePHiHreb2p8H/qD5U0BGoxE5OTmIiIiAwWDAsGHDoNPpsHfvXkRHR6O5uRmtra2ubyCd7HY7Pv74Y0yePBkA8Pvf/x4VFRVBmRUATp06hZKSEjz66KNYsmQJmpqaVM9aX18Pm82GP//5z1i4cKHHeQ8ePAiz2QwASE1Nxb///W/Y7fagzXvq1Cl8+OGHmDZtGubPn49Lly4FJGt9fT3y8vJgNptRWFgIp9PZZb6LFy+ira0NDzzwAAB19tn+5AWAmpoavPnmmzCbzcjPz4fValU968WLF/HAAw/gvvvuAwDk5ube9EslEJ8H/qD5AjB8+HDXAVFXV4fy8nIkJCTAYDCguLgYv/3tb2E0Gl3/+Z2++eYbxMTEuH4+G41GXL58OSizduZbsGABdu/ejfj4eOTn5wck6yuvvIIZM2bgrrvu8jhvQ0OD65SPXq9HTEwMrl69GrR5Y2NjkZGRgZKSEiQkJEg/tdZT1okTJ2LcuHFYvXo1iouLcfz4cezYsaPLfO7bFVBnn+1P3uvXr+P+++/HsmXLsGvXLjQ3N2Pjxo2qZ42IiEB0dDSeffZZmM1mvPrqqzd9yQrE54FfBPonSLA4c+aMeOSRR8TOnTu7vG6328WSJUvEK6+80uX1//3vf2LixIldphs5cmRQZu3u2rVrYvTo0TIjurhnPXz4sFiwYIEQQoijR496PKXyk5/8RNjtdtfziRMnioaGhqDN291DDz0kmpubZcYUQnjeD4QQYv/+/eKZZ57p8lp1dbVIT093Pa+rqxOTJ0+WnrNTX/N299lnn4mpU6fKiteFe9aNGzeK8ePHiwsXLgiHwyGWLVsmCgsLu0wfyM+D/tD8LwAAqK6uxuOPP47nnnsO06dPx6VLl1BdXQ2g4xtoSkoKTp8+3WWeQYMGobW1Fe3t7QCAxsZGVQa18yVrS0sL3nnnHddzIYQqDX/ds5aVleHs2bOYOnUqcnNzUVNTg0WLFt00n8lkcg0Y6HA40Nra6rp/RLDldTqdeO2111z7QSfZ27d71tOnT6OystL1fk//x90HYlRrn/U1b319fZdfBYHab7///e9j1KhRGDJkCMLDw/G73/0On376aZd5AvV50G+BrT+BV19fL8aOHSuqqqpcr50+fVo88sgjoqmpSTidTpGTk+NqlHQ3b948UVpaKoQQYuPGjeKFF14IyqwOh0NMmDBB/Pe//xVCCPHqq6+KvLw81bO68/aN+oUXXhCvvfaaEEKI3bt3i3nz5knL2ak/eadPny727t0rhOhotH7yySel5RSi56y1tbVi0qRJ4tq1a8Jms4knn3yyx0b/lJQUcfz4cSGEECtXrhRvvvmm1Kz9yWuxWMSYMWPEhQsXhNPpFMuXL5fe2N5T1vr6ejFx4kRRX18vhBDi+eefFwUFBTfNq/bngT9ofiygv/zlL3j//fdx9913u15LT0+HTqfD5s2bER4ejtGjR2PFihUwGAxYuXIlfv3rX+M3v/kNLl68iJycHFgsFsTHx+Ovf/0rbr/99qDMevz4cbz00ktoa2vDD3/4Q6xduxaxsbGqZ501axYA4KOPPkJRUZGrW+X69ethMpkwa9YsXLt2DTk5Ofjqq68QGxuLdevWeT0HH+i8Z8+eRV5eHlpaWjBo0CCsXbsW8fHxqmd1Op3YunUrHA4HkpKSsGTJEgDAvHnzkJ2djZ/+9Kf4/PPPkZubi+vXr7u6q3bvzhhMeSsrK/Hqq6/Cbrfj5z//OVatWiU1r6es8fHxKCgogNVqxf3334/Vq1cjKioqoJ8H/qD5AkBEpFVsAyAi0igWACIijWIBICLSKBYAIiKNYgEgItIoFgAiIo1iASAi0igWACIijfo/MhVNoZSpWEsAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"ax = pd.Series(mean_sampled).plot(kind='hist',bins=100) \n",
"ax.axvline(x = df_sub_M['age'].mean(), color='blue', lw=2)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.001"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len([i for i in mean_sampled if i >= df_sub_M['age'].mean()]) / len(mean_sampled)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Please repeat this procedure for Poplar. In this case, the results should suggest that the mean age of men is well within the range of the sampling distribution of the mean computed for women."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Permutation and Hypothesis Testing"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The bootstrap is valuable for understanding the distribution of a sample statistics\n",
"\n",
"The permutation test is valuable for hypothesis testing: central to the hypothesis test is the question of whether the differences we observe are the product of random chance? The latter will be our **null hypothesis**, namely that the age of men and women are essentially the same and that random chance explains the differences between the means. \n",
"\n",
"The alternative hypothesis serves as the \"counterpoint\" to the null hypothesis: the difference between the means is **not** the result of random chance. It is by rejecting the null hypothesis that we can accept the alternative hypothesis. \n",
"\n",
"In this sense, we never prove the alternative hypothesis directly but assume it is true because the null hypothesis is inconsistent with the data."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Permutation Procedure"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The permutation follows a slightly different procedure than the bootstrap but also involves resampling from our data. \n",
"\n",
"1: We note the number of observations we have for each of our categories of interest (i.e. the number of women and men present in the data frame). Let calls these two groups `M` and `F` for convenience and `len(F)` is the number of observations for label or group `F`.\n",
"\n",
"2: We combine the results into one data frame, which embodies the null hypothesis namely that, when it comes to age, the mean age of women and men are essentially the same, i.e. we can ignore the difference.\n",
"\n",
"3: We shuffle the combined dataset (which mean randomizing the order of the rows) and take a sample of size `len(F)` (notice that this will contain data from both `M` and `F`.\n",
"\n",
"4: The rows that remain are our second is a sample of `len(M)`\n",
"\n",
"5: Now, we compute the mean for each sample, take their difference, and record this number\n",
"\n",
"6: We repeat step 2 till 5 a number times, let's say 1000, which will generate a distribution of permutation statistics. These statistics will tell what differences a random permutation of the data will generate, that ignores the difference between the categories of interest. In other words, these values are likely to occur if the null hypothesis were true."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's implement this procedure in code. First, we create three data frames: one for 'M', one for 'F' and one combined 'F' and 'M'."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"regdist = 'Poplar'\n",
"df_sub_F = df[(df.gender=='F') & (df.district==regdist)]\n",
"df_sub_M = df[(df.gender=='M') & (df.district==regdist)]\n",
"df_sub_C = df[df.gender.isin(['F','M']) & (df.district==regdist)] "
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"((19993, 4), (19376, 4), (39369, 4))"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_sub_F.shape,df_sub_M.shape, df_sub_C.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For Poplar, the difference between means is approximately 0.2 years (women being younger than men). Is this result of random chance? What is the probability of observing this difference if we assume there exists no difference for this borough?"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(26.38463462211774, 26.603168868703552)"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_sub_F['age'].mean(),df_sub_M['age'].mean()"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.21853424658581133"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_sub_F['age'].mean() - df_sub_M['age'].mean()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To answer this question, we first record the number of observations for 'F' and then implement steps 2-5 in code."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"num_F = df_sub_F.shape[0] # the number of observations for women\n",
"all_idx = set(df_sub_C.index) # get all indices"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "184e3b0f8f8a405bb95f0020cb6a2a66",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"HBox(children=(FloatProgress(value=0.0, max=5000.0), HTML(value='')))"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"permutations = [] # initiate empty list to collect permutation statistics\n",
"\n",
"for _ in tqdm(range(5000)): # we use tqmd to print the progress, running the permutation test can take a while\n",
" f_idx = set(df_sub_C.sample(num_F).index) # take sample of size num_F and get index\n",
" m_idx = all_idx - f_idx # get the remaining indices (rows not f_idx)\n",
" diff_f_m_perm = df_sub_C.loc[f_idx]['age'].mean() - df_sub_C.loc[m_idx]['age'].mean() # diff between permuted statistics\n",
" permutations.append(diff_f_m_perm) # append permutation statistic"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can plot the distribution of the permutation statistics, which indicates that if we don't assume any age differences between men and women in Poplar is likely between -0.6 and 0.6."
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD7CAYAAACBiVhwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAfSElEQVR4nO3dfVRUdf4H8PfwIOqCWXSHPGw/t+1BltVkT21ouZhtAgojG3hMLa0NV8s27UlTgajM1cpWO8c8J3dNzYc2fAiIxcFdWTlrmCW1mS1lm0/5BINPCcrTzPf3h800cIG5A3MfZub9OsdzZi6XO2/Gmfu593u/3+81CSEEiIiI3IToHYCIiIyHxYGIiGRYHIiISIbFgYiIZFgciIhIhsWBiIhkWByIiEgmTO8AvnLuXAMcDt8M2YiOjsSZM/U+2ZYWmFddzKsu5lVXZ3lDQky4+uqfdPp7AVMcHA7hs+Lg3J4/YV51Ma+6mFdd3cnLZiUiIpJhcSAiIhkWByIikmFxICIiGRYHIiKSYXEgIiIZFgciIpIJmHEORFqL6tcHvSN+/Ao1NrXi4veXdUxE5DssDkTd1DsiDJZnilzPP3g9Axd1zEPkSywORF5of7bg7e/w7IL8BYsDBRS1d8TuZwsfvJ6hKAeANr9zsYN1WDTIaFgcKKC033nr1cyjpIgYJStRR9hbiYiIZFQtDvX19UhPT8fx48fbLN+4cSOmTJnien7y5Ek88MADSE1NxWOPPYaGhgY1YxGpornFDkmKgiRF6R2FqMdUKw6ff/45Jk2ahCNHjrRZ/r///Q9vvfVWm2UvvvgiJk+eDKvVisGDB2PlypVqxSJSTa/wUFieKWrTg4nIX6lWHAoKCpCfnw+z2exa1tzcjOeffx6zZ892LWtpacEnn3yClJQUAEBmZiasVqtasYiISAHVLkgvWrRItuz1119HVlYWfvrTn7qWnTt3DpGRkQgLuxJFkiTU1NSoFYuIiBTQrLfShx9+iFOnTmH+/PnYu3eva7kQ8jsUmUwmr7cfHR3Zo3zt+Vu7MfOq+zq+zNvZtrR4DaNiXnV1J69mxaGkpATffPMNMjIycOnSJdTV1eHJJ5/Ea6+9hvr6etjtdoSGhsJms7VpilLqzJl6n926T5KiYLP5T8fCYM/b1cA0X7yOe15f7BQ625av3pNg/zyoLVDyhoSYujyo1qw4LF682PV47969WLFiBZYvXw4AuP3221FaWgqLxYLCwkIkJSVpFYsCgNKBaUSknCHGOeTn56OgoABjx47Fvn378OSTT+odiUgV7O5K/kL1M4fy8nLZssTERCQmJrqex8bGYv369WpHIQKg77QVzu6uAM9yyNg4fQYFHU5bQeQZiwP5pe7MjkpEyhnimgORt5xH/xyNTKQOFgciIpJhcSAiIhk22hL9oLNeTM7up0TBhMWB6Aed9WJi91MKRmxWIiIiGRYHIiKSYbMSBTVeTyDqGM8cKKjx7m1EHWNxICIiGRYHIiKS4TUHCgqci4nIOzxzoKDAuZiIvMNDKTK0ntx7gT2RiLqPxYEMrSf3XuDIZqLuY7MSERHJsDgQEZEMiwMREcmoXhzq6+uRnp6O48ePAwDee+89pKenw2KxYP78+WhubgYAVFdXIysrCykpKcjJyUFra6va0YiIqBOqFofPP/8ckyZNwpEjRwAAhw8fxurVq/G3v/0NxcXFcDgc2LRpEwBgzpw5yMvLQ1lZGYQQKCgoUDMaERF1QdXiUFBQgPz8fJjNZgBAr1698MILLyAyMhImkwm33HILTp48iRMnTqCxsREJCQkAgMzMTFitVjWjEXXJ2Q2WXWEpWKnalXXRokVtnsfGxiI2NhYAcPbsWWzcuBGLFy9GbW0tJElyrSdJEmpqarx6rejoyJ4HduNvO4VgyavV36lHN1hf/m3B8nnQSzDk1WWcQ01NDaZNm4asrCwkJibi008/la1jMpm82uaZM/VwOIRP8klSFGw2b3rU6yuQ87b/UDt/z9++nEr46v8wkD8PRhAoeUNCTF0eVGveW+nbb7/FpEmTcN999+Hxxx8HAMTExKCurs61js1mczVFERGR9jQtDvX19cjOzsbs2bPxyCOPuJbHxsYiIiICVVVVAIDCwkIkJSVpGY2IiNxo2qy0ZcsW1NXV4e2338bbb78NALjnnnswe/ZsLF26FLm5uWhoaEB8fDymTp2qZTQiInKjSXEoLy8HADz88MN4+OGHO1wnLi4OW7Zs0SIOERF5wBHSREQkw+JAREQyLA5ERCTD4kBERDK82Q+RAbjftc7bO94RqYHFgcgA2k/X4T/jbylQsTiQ4bjfN5qI9MFrDmQ4zvtGO4+kiUh7LA5ERCTD4kBERDJs2CW/4d6jJ5Cx5xIZAYsD+Q09bsCjB/ZcIiNgsxIREcmwOBARkQyblYgMjNcfSC8sDkQGxusPpBc2KxERkQyLAxERybA4EBGRjOrFob6+Hunp6Th+/DgAoLKyEhaLBcnJyVi2bJlrverqamRlZSElJQU5OTlobW1VOxoREXVC1eLw+eefY9KkSThy5AgAoLGxEQsWLMDKlStRWlqKAwcOoKKiAgAwZ84c5OXloaysDEIIFBQUqBmNDCaqXx9IUlRQjIAm8geqFoeCggLk5+fDbDYDAPbv34+BAwfi+uuvR1hYGCwWC6xWK06cOIHGxkYkJCQAADIzM2G1WtWMRgbDmViJjEXVrqyLFi1q87y2thaSJLmem81m1NTUyJZLkoSamhqvXis6OrJnYdvxtyNY5g0OSt83f3t/mVdd3cmr6TgHIYRsmclk6nS5N86cqYfDId9Od0hSFGw2/+lRHgh5/e3Lphcl/8+B8HkwskDJGxJi6vKgWtPiEBMTg7q6Otfz2tpamM1m2XKbzeZqiqLAxTu+eYejpUlLmnZlHTp0KA4fPoyjR4/CbrejpKQESUlJiI2NRUREBKqqqgAAhYWFSEpK0jIa6YDXGbzjHC1teaaIRZVUp+knLCIiAkuWLMETTzyBpqYmjBw5EqmpqQCApUuXIjc3Fw0NDYiPj8fUqVO1jEZERG40KQ7l5eWux8OHD0dxcbFsnbi4OGzZskWLOERE5IGiZqX169ejvr5e7SxEpJDz+oMkRSGqXx+941AAUlQcDh486Bq5/MUXX6idiYg84PUHUpui4rBw4UKUlZVh8ODBePHFF5GVlYUtW7agqalJ7XxERKQDxb2VIiMjkZqaivT0dJw/fx6bNm1CamoqRzITEQUgReejlZWVKCgowJ49e5CSkoI333wTcXFxOHbsGCZPnuzqcURERIFBUXF46aWXMHnyZCxcuBBRUT+OZP2///s/TJgwQbVwRESkD0XNSsXFxejfvz+ioqJgs9mwdu1aOBwOAMCsWbNUDUhERNpTfEF6165dV34hJARVVVX405/+pGYuIiLSkaJmpc8++wwlJSUAgOjoaLzxxhvIyMhQNRgREelH0ZlDS0sLmpubXc95lzYiosCm6Mzh7rvvRnZ2NjIyMmAymVBSUoKRI0eqnY2IvOQcLS1JUZy5lXpEUXGYO3cuNm7ciJ07dyIsLAyjR4/GxIkT1c5GRAq4T+UNwDXL7QevZ8B/7jpARqOoOISGhmLq1KmcKZXIgJxTaQBXCgKRLygqDqWlpVi6dCkuXLjQ5q5tn376qWrBiIhIP4qKwxtvvIF58+YhPj7e69t3EhGR/1FUHPr164fk5GS1sxARkUEo6so6dOhQVFRUqJ2FiIgMQtGZQ0VFBTZs2IDw8HCEh4dDCAGTycRrDkREAUpRcVi7dq3KMYiIyEgUNSvFxsbiiy++QEFBAa655hp89tlniI2N7faLFhUVIS0tDWlpaXjllVcAANXV1cjKynLdcY6jsImI9KOoOKxatQrvvvsurFYrGhsbsWLFCrz55pvdesHLly9j0aJFWL9+PYqKirBv3z5UVlZizpw5yMvLQ1lZGYQQKCgo6Nb2ydjcR/ASkXEpKg5///vf8Ze//AV9+vTB1VdfjYKCAtdEfN6y2+1wOBy4fPkyWltb0drairCwMDQ2NiIhIQEAkJmZyTvMBajeEWGuex8TkXEpuuYQFhaGXr16uZ7369cPYWHdu6l5ZGQkZs+ejTFjxqB379644447EB4eDkmSXOtIkoSamppubZ+IiHpO0R5+wIAB2LVrF0wmE5qbm7F69epuX3P46quvsHXrVvzrX/9CVFQUnn32WXz44Yey9bwdbBcdHdmtPJ3xt2YPf8tL2vCXz4W/5HQKhryKikNeXh7mzp2Lr7/+GgkJCRg6dCiWLl3q9YsBwO7duzF8+HBER0cDuNKEtHr1atTV1bnWsdlsMJvNXm33zJl6OBzC84oKSFIUbDb/mbLMn/L625fK3/nD58KfPr9A4OQNCTF1eVCtqDjExMRg3bp1uHz5Mux2OyIju3+UHhcXh9deew2XLl1Cnz59UF5ejjvuuANlZWWoqqrCbbfdhsLCQiQlJXX7NYiIqGcUFYc1a9Z0uPz3v/+91y84YsQI/Pe//0VmZibCw8MxZMgQTJ8+HaNHj0Zubi4aGhoQHx/PGWCJiHSkqDgcPHjQ9bi5uRlVVVVITEzs9otOnz4d06dPb7MsLi4OW7Zs6fY2iYjIdxQVh8WLF7d5fvbsWcydO1eVQEREpD9F4xzau+aaa3DixAlfZyEiIoPw+pqDEAIHDhxw9TYiIqLA4/U1B+DKuAc2KxEZm/u9pRubWnHx+8s6JyJ/0q1rDkRkfO3vLe0/PfPJCBQVhylTpnQ5Yvmdd97xWSAiItKfouIwePBgfPvtt5gwYQLCw8NRVFSE1tZWpKWlqZ2PiIh0oKg4fPrpp9i0aRNCQ0MBAL/5zW8wYcIEpKSkqBqOAkNUvz7oHdG9iRqJSB+KurKePXsWzc3NrucNDQ1obGxULRQFFk7TTeR/FB3OpaenY8KECRg9ejSEENi+fTunt6Au8WyByL8p+vbOnj0b8fHx+OijjxAREYGXXnoJd9xxh9rZyI85zxaAKz1liMi/KB4hHRMTg5tvvhlPPvkkwsPD1cxEREQ6U1Qctm7divnz5+Ovf/0rLl68iJkzZ/IezwTgSvORJEVBkqJc94cmIv+nqDhs2LAB7733HiIjIxEdHY1t27Zh3bp1amcjP+B+sZnXGIgCh6Jvc0hISJsb/AwYMMDVrZWIjI9TaZC3FBWH/v37o7q62jVKuri4GFdddZWqwYjIdziVBnlLUXFYsGABZs+ejWPHjmHEiBGIiIjAypUr1c5GREQ6UVQcGhsbUVRUhCNHjsBut+OGG25gjyUiogCm6IL0s88+i9DQUNx444245ZZbWBiIiAKcouIwaNAgfPDBBzh58iTOnz/v+tdd5eXlyMzMRGpqKl5++WUAQGVlJSwWC5KTk7Fs2bJub5uIiHpOUbPSzp07YbVa2ywzmUyorq72+gW/++475OfnY/PmzYiOjsZDDz2EiooK5OfnY/369RgwYABmzJiBiooKjBw50uvtExFRzykqDl988YXPXvAf//gHxo4di+uuuw4AsGzZMhw9ehQDBw7E9ddfDwCwWCywWq0sDn7GvbskEfm3LpuV8vLyXI/Pnj3rkxc8evQo7HY7srOzMW7cOGzatAm1tbWQJMm1jtlsRk1NjU9ej7Tj7C7J2VeJ/F+XZw4HDhxwPc7Ozsb777/f4xe02+3Yt28f1q9fj759+2LmzJno00c+7UJXd57rSHR0pOeVvOBvR8D+lpf0ZbTPi9HyeBIMebssDkKIDh/3xLXXXovhw4fjmmuuAQD89re/hdVqbTPiura2Fmaz2avtnjlTD4fDNxklKQo2m/8ME9Izr799SegKI32++X1TV2d5Q0JMXR5UK56V1dsj+c6MGjUKu3fvxvfffw+73Y5///vfSE1NxeHDh11NTiUlJUhKSvLJ6xERkfe6PHNwOBy4cOEChBCw2+2ux079+/f3+gWHDh2KadOmYfLkyWhpacFdd92FSZMm4ec//zmeeOIJNDU1YeTIkUhNTfX+ryEiIp/osjgcPHgQw4YNcxWExMRE18+625UVAMaPH4/x48e3WTZ8+HAUFxd3a3tE1HPud+/j5HzUZXH46quvtMpBRDpofztXTs5HToqvORBR4HG/HweROxYHIiKS4a27yGvtmyKIKPDwzIG8xqYIosDHwz+iIMM5sEgJnjkQBRnOgUVKsDgQEZEMiwMREcmwOBARkQyLAxERybA4EBGRDIsDERHJsDgQEZEMB8ERkYz7QDlO3x2cWByISMY5UA7g9N3Bis1KREQkw+JAREQyLA5ERCSja3F45ZVXMG/ePABAdXU1srKykJKSgpycHLS2tuoZjdxE9esDSYpy/SOiwKdbcdizZw/ef/991/M5c+YgLy8PZWVlEEKgoKBAr2jUjvv9GziTJ1Fw0KU4nD9/HsuWLcOjjz4KADhx4gQaGxuRkJAAAMjMzITVatUjGhERQafi8Pzzz+Opp55Cv379AAC1tbWQJMn1c0mSUFNTo0c0IiKCDuMcNm/ejAEDBmD48OHYtm0bAEAIIVvPZDJ5td3o6Eif5HPyt7Z1f8tL/kXtz5e/fX6DIa/mxaG0tBQ2mw0ZGRm4cOECLl26BJPJhLq6Otc6NpsNZrPZq+2eOVMPh0NeZLpDkqJgs/nPsB+18/rbF4F8T+3PF79v6uksb0iIqcuDas2Lw5o1a1yPt23bho8//hiLFy9Geno6qqqqcNttt6GwsBBJSUlaRyMioh8YZvqMpUuXIjc3Fw0NDYiPj8fUqVP1jkREFLR0LQ6ZmZnIzMwEAMTFxWHLli16xglKUf36oHfElY8BJ1gjIifDnDmQPpxjGABOsEZEP2JxIKIucfru4MTiQC7uO4GmZjsieoXqnIiMgNN3BycWhyDkfp3BXfudgPtjIoBnEcGExSEItb/OQKQUzyKCB6fsJiIiGZ45EJFPtW+2ZPOTf2JxICKfcm+2BNj85K9YHIJEZxehiYg6wmsOQcL9hj1ERJ6wOBARkQyLAxERybA4EBGRDIsDERHJsDgQEZEM+zYSUY+xq3Tg4ZkDEfUYu0oHHpb6AMCZMonI11gcAgBnyiQiX9OlWWnFihVIS0tDWloaXn31VQBAZWUlLBYLkpOTsWzZMj1iEZEKnGe2khSFqH599I5DCmleHCorK7F79268//77KCwsxJdffomSkhIsWLAAK1euRGlpKQ4cOICKigqtoxGRF9x3+l1xntlaniniRWs/onlxkCQJ8+bNQ69evRAeHo4bb7wRR44cwcCBA3H99dcjLCwMFosFVqtV62hE5AX3nT4FHs2Lw80334yEhAQAwJEjR1BaWgqTyQRJklzrmM1m1NTUaB2NiIh+oNs53jfffIMZM2bgueeeQ1hYGA4fPtzm5yaTyavtRUdH+jKex1NlI3Nmb26xo1d4qM5piNrq6Lvlb9+3YMirS3GoqqrCrFmzsGDBAqSlpeHjjz9GXV2d6+e1tbUwm81ebfPMmXo4HMIn+SQpCjab//T5af8f78wuSVG8VzQZTvvvlj9+3wIhb0iIqcuDas2blU6dOoXHH38cS5cuRVpaGgBg6NChOHz4MI4ePQq73Y6SkhIkJSVpHY2IiH6g+ZnD6tWr0dTUhCVLlriWTZw4EUuWLMETTzyBpqYmjBw5EqmpqVpHCwjuA+KIjIYDNv2H5sUhNzcXubm5Hf6suLhY4zSBp/2AOCIj4YBN/8G5lYiISIYjUgzOfbZL99NwzoJJRGri3sXgnLNdAm1Pw9svJyLyJTYrERGRDIsDERHJsFnJgHg9gYj0xj2QAfF6AgUD9zEP/a7qi4heV6Z64fgHY2BxICJdtB/zwPEPxsLiQESGwlHUxsDi4Ec4NQYFg85GUXc25ofUweLgRzg1BgWzzsb8kDrYlZWIiGR45mAQ7L5KpByvS6iPeyODYPdVIuU4u6v62KxEREQyPHPQmHvzUVOz3TXwh4i6h01M6mBxUMhX3ejaNx+xKYmoc0q6b7OJSR0sDgqxGx2R9th9Wz8sDj3UWTMRb8xDpL3Ompg4gM57Qb/H6umHpqtmIt6Yh0hb7mcaW5ekt2mS8jTq2v3grrMDvWBiqN5KH3zwAcaOHYvRo0dj48aNmrymc8dteaaIR/dEAcRZKJxFwcl5diFJUW2+/xG9Qjt8HKz7BcP81TU1NVi2bBm2bduGXr16YeLEiUhMTMRNN92kdzQiCiC8jqGMYYpDZWUlhg0bhv79+wMAUlJSYLVa8cc//lHR74eEmLr92uar+wBo214pSVFoampFfX1jl+u5L/d2HX99bJQcRntslBxGe2yUHD157Ny/REb2BiDfP0RG9kaEs3mq3X7DyX0doF0zVie/0xklr+euo/2jp32mSQghFCdS0VtvvYVLly7hqaeeAgBs3rwZ+/fvx8KFC3VORkQUfAxzzaGjGmUydf9sgIiIus8wxSEmJgZ1dXWu57W1tTCbzTomIiIKXoYpDnfeeSf27NmDs2fP4vLly9ixYweSkpL0jkVEFJQMc0E6JiYGTz31FKZOnYqWlhaMHz8et956q96xiIiCkmEuSBMRkXEYplmJiIiMg8WBiIhkWByIiEiGxYGIiGRYHACcPHkSDzzwAFJTU/HYY4+hoaFBtk5zczNefvll/O53v0NaWhp2796tQ9IrlOR1qq+vx7333ou9e/dqmLAtJXlra2uRnZ2NjIwM3HfffdizZ4/mOT1N/FhdXY2srCykpKQgJycHra2tmmd05ynvP//5T2RkZGDcuHGYOXMmLly4oEPKHymdWHPXrl245557NEzWMU95Dx06hClTpmDcuHHIzs7W9f31lPXLL79EVlYWxo0bhxkzZuD777/3vFFBYvr06aKkpEQIIcSKFSvEq6++KlvnzTffFE8//bRwOBzi4MGDYsSIEcLhcGgdVQihLK/T3Llzxa9//Wvx0UcfaRVPRkneZ555Rqxfv14IIcS3334r7rzzTtHa2qpZxtOnT4tRo0aJc+fOiYaGBmGxWMQ333zTZp20tDTx2WefCSGEmD9/vti4caNm+drzlPfixYvirrvuEqdPnxZCCLF8+XKxcOFCveIqen+FEMJms4nU1FQxatQoHVL+yFNeh8MhkpOTRUVFhRBCiNdee63L76GeWYUQYtKkSWLXrl1CCCEWL14s/vznP3vcbtCfObS0tOCTTz5BSkoKACAzMxNWq1W23vbt2/GHP/wBJpMJN998M9asWdPhlB9qU5oXAEpLS/GTn/wEgwYN0jJiG0rzJicnw2KxAAAGDhyIpqYmXLp0SbOc7hM/9u3b1zXxo9OJEyfQ2NiIhIQEAF2/71rwlLelpQUvvPACYmJiAACDBg3CqVOn9IrrMa9Tbm6u4sk21eQp75dffom+ffu6Buo++uijeOCBBwyZFQAcDofrjP3y5cvo3bu3x+0GfXE4d+4cIiMjERZ2ZTygJEmoqamRrXf06FF88sknyMzMxP3334+6ujqEhGj/9inNe/LkSaxbtw5z587VOmIbSvMmJyfjqquuAgCsXr0av/jFLxAV1fW9g32ptrYWkiS5npvN5jY52/+8s79DK57yXn311bj33nsBAI2NjVi1apXruR485QWAd955B/Hx8Rg6dKjW8WQ85T127BiuvfZaPPfcc7BYLMjPz0ffvn31iKrovZ03bx5ycnIwYsQIVFZWYuLEiR63a5gR0lrYvn07Fi9e3GbZz372M9l6HU34Z7fbcfr0aWzduhVff/01pk2bhu3bt6u6A+tuXofDgZycHOTl5Sk6QvCVnry/TmvXrsV7772HDRs2+Dpelzo6C3TP6ennWlOa5+LFi5g5cybi4uJw3333aRGtQ57yHjx4EDt27MDatWtx+vRpLaN1yFPe1tZWfPzxx9iwYQOGDBmC5cuXY8mSJViyZImWMQF4ztrY2IicnBysW7cOt956K9asWYPnnnsOq1at6nK7QVUcxowZgzFjxrRZ1tLSgsTERNjtdoSGhsJms3U44d+1116LtLQ0mEwmxMXF4brrrsPhw4dVneKju3kPHTqEQ4cOIScnB8CVo5zc3FwsXLgQw4YNM1xep1dffRUVFRXYuHEjrrvuOtVydiQmJgb79u1zPW8/8WP7iSG7+ju04Cmvc1l2djaGDRuGBQsWaB2xDU95rVYrbDYbsrKy0NLSgtraWkyePBmbNm3SI67HvJIkYeDAgRgyZAgAID09HbNmzdI8J+A568GDBxEREeHaV91///144403PG436JuVwsPDcfvtt6O0tBQAUFhY2OGEf6NGjXKt89133+HUqVO44YYbNM0KKMt70003oaKiAkVFRSgqKsLgwYPx8ssvq1oYepIXuHLGsHfvXrz77ruaFwbA88SPsbGxiIiIQFVVFYDO/w6teMprt9vx6KOPYsyYMcjJydF9+ntPeWfNmoWysjIUFRVh1apVMJvNuhUGJXl/9atf4ezZs/jqq68AAOXl5fjlL39pyKwDBw7E6dOncejQIQDAzp07XUWtS769bu6fjh8/Lh588EExZswY8cgjj4jz588LIYTYtGmTWL58uRDiSu+POXPmiLFjx4qxY8eK8vJyQ+d19+CDD+raW8lTXofDIW6//XZx9913i3Hjxrn+OXvaaKW4uFikpaWJ5ORksWrVKiGEENOmTRP79+8XQghRXV0tsrKyRGpqqnj66adFU1OTpvna6yrvjh07xKBBg9q8nwsWLDBsXnffffed7r2VhPCc9z//+Y/IysoSY8eOFY888oioq6szbNZdu3YJi8Ui0tPTxUMPPSSOHTvmcZuceI+IiGSCvlmJiIjkWByIiEiGxYGIiGRYHIiISIbFgYiIZFgciIhIhsWBiIhkWByIiEjm/wHd7N5tjloKUgAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"pd.Series(permutations).plot(kind='hist',bins=50)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The observed difference was -0.2 which is well the range of what we'd expect if the null hypothesis were true. In this case we don't have enough evidence to reject the null-hypothesis. "
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'Frequency')"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAVgAAAE/CAYAAAADnduiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de1zT5eIH8M9gODQwgzbyxa88HW9kmHQOJy8ZXkoZAqJQppjo8VbaTStvXCI1g4q8lOHRMjXTiqgEiaYlyUmxUuqUFmV5wTsM8MKAwdie3x8edkRMx+XZBn7er1evF3v23fP9bIyP67t9nymEEAJERNTiXBwdgIiorWLBEhFJwoIlIpKEBUtEJAkLlohIEhYsEZEkLFgiIkmUjg4g29mzFbBYmvZRX29vD5SWGlo4UfM4OtPf/+4PAMjPP2Adc3SmK2Gma3O2PEDryOTiosBNN91g023bfMFaLKLJBVt3e2fjyEyFhYVXzMDHyTbOlsnZ8gBtKxMPERARScKCJSKShAVLRCQJC5aISBIWLBGRJCxYIiJJWLBERJKwYImIJGHBEhFJ0ubP5KLrm2fH9nBXtezT3Fhdi/ILVS06J7VNLFhq09xVSoQ/m9Gic259LQLlLTojtVU8REBEJAkLlohIEhYsEZEkLFgiIklYsEREkrBgiYgkkVqwBoMBYWFhOHHiBADgww8/RFhYGMLDw7FgwQLU1NQAAAoKChAVFYXg4GDExcWhtrYWAHDq1CmMHz8eWq0WM2bMQEVFhcy4REQtSlrB/vjjjxg3bhyOHj0KADhy5AjWrl2LDz74AJmZmbBYLNi8eTMAYM6cOUhISMC2bdsghEBaWhoAYOHChYiOjoZOp4O/vz9SU1NlxSUianHSCjYtLQ2JiYnQaDQAgHbt2uGFF16Ah4cHFAoFevTogVOnTuHkyZMwGo0ICAgAAERGRkKn08FkMmHv3r0IDg6uN05E1FpIO5NryZIl9S77+vrC19cXAFBWVoZNmzYhKSkJxcXFUKvV1u3UajWKiopw9uxZeHh4QKlU1hsnImot7H6qbFFREaZOnYqoqCj07dsX33//fYNtFAoFhGj4LY4KhaLR+/P29mhSzjpqtWezbi+DM2S6PENzMtWYzGjn5trcSA3mlPk4NXVuZ/jdXcrZ8gBtK5NdC/bQoUOYNm0aHnnkEUyePBkA4OPjg5KSEus2er0eGo0GXl5eMBgMMJvNcHV1tY43VmmpoclfuatWe0Kvd66zzp0l06UZmptJrfaUsl6AXl8u7Y+1KffXWX53dZwtD9A6Mrm4KGx+4Wa3j2kZDAZMmTIFTz/9tLVcgYuHDlQqFfLz8wEAW7ZsQVBQENzc3BAYGIjs7Ox640RErYXdCjY9PR0lJSV45513EBERgYiICKxYsQIAkJKSgqSkJISEhKCqqgoxMTEAgMTERKSlpWHEiBHYt28fZs2aZa+4RETNJv0QQU5ODgBg0qRJmDRp0hW38fPzQ3p6eoNxX19fbNy4UWY8IiJpeCYXEZEkLFgiIklYsEREkrBgiYgkYcESEUnCgiUikoQFS0QkCQuWiEgSFiwRkSQsWCIiSViwRESSsGCJiCRhwRIRScKCJSKShAVLRCQJC5aISBIWLBGRJCxYIiJJWLBERJKwYImIJGHBEhFJwoIlIpKEBUtEJAkLlohIEhYsEZEkLFgiIklYsEREkrBgiYgkYcESEUnCgiUikoQFS0QkCQuWiEgSFiwRkSRSC9ZgMCAsLAwnTpwAAOTl5SE8PBzDhw/HsmXLrNsVFBQgKioKwcHBiIuLQ21tLQDg1KlTGD9+PLRaLWbMmIGKigqZcYmIWpS0gv3xxx8xbtw4HD16FABgNBoRGxuL1NRUZGdn48CBA8jNzQUAzJkzBwkJCdi2bRuEEEhLSwMALFy4ENHR0dDpdPD390dqaqqsuERELU5awaalpSExMREajQYA8NNPP6FLly649dZboVQqER4eDp1Oh5MnT8JoNCIgIAAAEBkZCZ1OB5PJhL179yI4OLjeOBFRa6GUNfGSJUvqXS4uLoZarbZe1mg0KCoqajCuVqtRVFSEs2fPwsPDA0qlst54Y3l7ezTxHtTl8WzW7WVwhkyXZ3CGTJeTmampczvb4+RseYC2lUlawV5OCNFgTKFQNHq8sUpLDbBYGs5lC7XaE3p9eZNuK4uzZLo0Q3MzyfqD0uvLpc7dWM7yu6vjbHmA1pHJxUVh8ws3u32KwMfHByUlJdbLxcXF0Gg0Dcb1ej00Gg28vLxgMBhgNpvrjRMRtRZ2K9g+ffrgyJEjKCwshNlsRlZWFoKCguDr6wuVSoX8/HwAwJYtWxAUFAQ3NzcEBgYiOzu73jgRUWtht0MEKpUKycnJePLJJ1FdXY1BgwZBq9UCAFJSUhAfH4+Kigr06tULMTExAIDExETMnz8fq1atQufOnbF06VJ7xSUiajbpBZuTk2P9uX///sjMzGywjZ+fH9LT0xuM+/r6YuPGjVLzERHJwjO5iIgkYcESEUnCgiUikoQFS0QkCQuWiEgSFiwRkSQsWCIiSViwRESSsGCJiCRhwRIRScKCJSKShAVLRCQJC5aISBIWLBGRJCxYIiJJWLBERJKwYImIJGHBEhFJwoIlIpLEbl96SK2bZ8f2cFf97+miVnvWu/7yy7YwVtei/EJVs7MROSsWLNnEXaVE+LMZ1suX/txUW1+LQHmzZyFyXjxEQEQkCQuWiEgSFiwRkSQsWCIiSViwRESSsGCJiCRhwRIRScKCJSKShAVLRCQJC5aISBIWLBGRJA4p2IyMDISGhiI0NBQvv/wyAKCgoABRUVEIDg5GXFwcamtrAQCnTp3C+PHjodVqMWPGDFRUVDgiMhFRo9m9YKuqqrBkyRJs3LgRGRkZ2LdvH/Ly8jBnzhwkJCRg27ZtEEIgLS0NALBw4UJER0dDp9PB398fqamp9o5MRNQkdi9Ys9kMi8WCqqoq1NbWora2FkqlEkajEQEBAQCAyMhI6HQ6mEwm7N27F8HBwfXGiYhaA7svV+jh4YGnn34aISEhcHd3xz333AM3Nzeo1WrrNmq1GkVFRTh79iw8PDygVCrrjTeGt7dHs/I2ZZ1T2ZwxU1PJvC/OOLez/e6cLQ/QtjLZvWB//fVXfPzxx/jqq6/g6emJ5557Drt3726wnUKhgBDiiuONUVpqgMXScB5bqNWe0Ouda8VSR2WS9aTX68tb7dyN5WzPJ2fLA7SOTC4uCptfuNn9EMGuXbvQv39/eHt7o127doiMjMS3336LkpIS6zZ6vR4ajQZeXl4wGAwwm831xomIWgObCnbjxo0wGAwtskM/Pz/k5eWhsrISQgjk5OTgnnvugUqlQn5+PgBgy5YtCAoKgpubGwIDA5GdnV1vnIioNbDpEMHBgwcRHByMwYMHY+zYsejdu3eTdzhw4ED88ssviIyMhJubG3r37o3p06dj2LBhiI+PR0VFBXr16oWYmBgAQGJiIubPn49Vq1ahc+fOWLp0aZP3TURkTzYV7OLFizFv3jxs3boVCxcuhBAC48aNQ3h4OFQqVaN3On36dEyfPr3emJ+fH9LT0xts6+vri40bNzZ6H0REjmbzMVgPDw9otVqEhYXh3Llz2Lx5M7RaLT82RUT0J2x6BZuXl4e0tDTs2bMHwcHBePPNN+Hn54djx44hOjoaWq1Wdk4iolbHpoJdtGgRoqOjsXjxYnh6/u9jL7fddhvGjBkjLRwRUWtm0yGCzMxMdOrUCZ6entDr9Vi/fj0sFgsA4KmnnpIakIiotbL5Ta6KigqMHDkSLi4uyM/Px4kTJxAfHy87H5HT8uzYHu6qxp2rc60TH4zVtSi/UNWcWOREbHp2/PDDD8jKygIAeHt7Y8WKFYiIiJAajMjZuauUCH82o0Xn3PpaBJzrPCZqDpsOEZhMJtTU1Fgv1y0lSEREf86mV7CDBw/GlClTEBERAYVCgaysLAwaNEh2NiKiVs2mgp07dy42bdqEHTt2QKlUYtiwYRg7dqzsbERErZpNBevq6oqYmBjr6atERHRtNhVsdnY2UlJScP78+XpLCH7//ffSghERtXY2FeyKFSswf/589OrVq9HrsRIRXa9sKtiOHTti+PDhsrMQEbUpNn1Mq0+fPsjNzZWdhYioTbHpFWxubi7ee+89uLm5wc3NDUIIKBQKHoMlIroKmwp2/fr1kmMQEbU9Nh0i8PX1xf79+5GWlgYvLy/88MMP8PX1lZ2NiKhVs6lg16xZg/fffx86nQ5GoxErV67Em2++KTsbEVGrZlPBfvbZZ3jrrbfQvn173HTTTUhLS7Mu/kJERFdmU8EqlUq0a9fOerljx45QKhu3TBsR0fXGppbs3Lkzdu7cCYVCgZqaGqxdu5bHYImIrsGmgk1ISMDcuXPx22+/ISAgAH369EFKSorsbERErZpNBevj44MNGzagqqoKZrMZHh4esnMREbV6NhXsunXrrjj+z3/+s0XDEBG1JTYV7MGDB60/19TUID8/H3379pUWioioLbCpYJOSkupdLisrw9y5c6UEIiJqK2z6mNblvLy8cPLkyZbOQkTUpjT6GKwQAgcOHIC3t7e0UEREbUGjj8ECFz8Xy0MERERX16RjsEREdG02FeyECROu+lUx7777bosFIiJqK2wqWH9/fxw6dAhjxoyBm5sbMjIyUFtbi9DQUNn5iIhaLZsK9vvvv8fmzZvh6uoKALjvvvswZswYBAcHN2mnOTk5WLlyJSorKzFw4EDEx8cjLy8PSUlJqK6uRkhICGbPng0AKCgoQHx8PAwGAwIDA7Fw4UIuNENErYJNH9MqKytDTU2N9XJFRQWMRmOTdnj8+HEkJiYiNTUVW7duxS+//ILc3FzExsYiNTUV2dnZOHDggPU7wObMmYOEhARs27YNQgikpaU1ab9ERPZm00vBsLAwjBkzBsOGDYMQAp9//jliYmKatMMvvvgCI0aMwC233AIAWLZsGQoLC9GlSxfceuutAIDw8HDodDp069YNRqMRAQEBAIDIyEi8/vrriI6ObtK+iYjsyaaCffrpp9GrVy988803UKlUWLRoEe65554m7bCwsBBubm6YMmUK9Ho9hgwZgu7du0OtVlu30Wg0KCoqQnFxcb1xtVqNoqKiRu3P27t5C9Oo1Z7Nur0MzpipqWTeF87tXPuyVVvKZPPBTB8fH3Tv3h2RkZH4+eefm7QzADCbzdi3bx82btyIDh06YObMmWjfvn2D7RQKBYQQVxxvjNJSAyyWhvPYQq32hF5f3qTbyuKoTLKe9Hp9Oee+wtz2wOe3bS7P5OKisPmFm03HYD/++GMsWLAAb7/9NsrLyzFz5swmHwu9+eab0b9/f3h5ecHd3R33338/du/ejZKSEus2xcXF0Gg08PHxqTeu1+uh0WiatF8iInuzqWDfe+89fPjhh/Dw8IC3tzc++eQTbNiwoUk7HDJkCHbt2oULFy7AbDbj66+/hlarxZEjR1BYWAiz2YysrCwEBQXB19cXKpUK+fn5AIAtW7YgKCioSfslIrI3mw4RuLi41Ftku3PnztaPbDVWnz59MHXqVERHR8NkMuHee+/FuHHj8Ne//hVPPvkkqqurMWjQIGi1WgBASkoK4uPjUVFRgV69ejX5zTUiInuzqWA7deqEgoIC6/HPzMxM3HjjjU3e6YMPPogHH3yw3lj//v2RmZnZYFs/Pz+kp6c3eV9ERI5iU8HGxsbi6aefxrFjxzBw4ECoVCqkpqbKzkZE1KrZVLBGoxEZGRk4evQozGYzbr/9dri5ucnORkTUqtn0Jtdzzz0HV1dXdO3aFT169GC5EhHZwKaC7dmzJ7Zu3YpTp07h3Llz1v+IiOjP2XSIYMeOHdDpdPXGFAoFCgoKpIQiImoLbCrY/fv3y85BRNTmXPUQQUJCgvXnsrIy6WGIiNqSqxbsgQMHrD9PmTJFehgiorbkqgV76WIrV1p4hYiI/pxNnyIAGr+KFRHR9e6qb3JZLBacP38eQgiYzWbrz3U6deokPSARUWt11YI9ePAg+vXrZy3Vvn37Wq/jx7SIiK7uqgX766+/2isHEVGbY/MxWCIiahwWLBGRJCxYIiJJWLBERJKwYImIJGHBEhFJwoIlIpKEBUtEJAkLlohIEhYsEZEkLFgiIklYsEREkrBgiYgkYcESEUnCgiUikoQFS0QkCQuWiEgSFiwRkSQsWCIiSRxWsC+//DLmz58PACgoKEBUVBSCg4MRFxeH2tpaAMCpU6cwfvx4aLVazJgxAxUVFY6KS0TUaA4p2D179uDTTz+1Xp4zZw4SEhKwbds2CCGQlpYGAFi4cCGio6Oh0+ng7++P1NRUR8QlImoSuxfsuXPnsGzZMjz22GMAgJMnT8JoNCIgIAAAEBkZCZ1OB5PJhL179yI4OLjeOBFRa2H3gn3++ecxe/ZsdOzYEQBQXFwMtVptvV6tVqOoqAhnz56Fh4cHlEplvXEiotZCac+dffTRR+jcuTP69++PTz75BAAghGiwnUKh+NPxxvL29mh80Euo1Z7Nur0MzpipqWTeF87tXPuyVVvKZNeCzc7Ohl6vR0REBM6fP4/KykooFAqUlJRYt9Hr9dBoNPDy8oLBYIDZbIarq6t1vLFKSw2wWBqWtS3Uak/o9eVNuq0sjsok60mv15dz7ivMbQ98ftvm8kwuLgqbX7jZ9RDBunXrkJWVhYyMDDz11FMYOnQokpKSoFKpkJ+fDwDYsmULgoKC4ObmhsDAQGRnZ9cbJyJqLZzic7ApKSlISkpCSEgIqqqqEBMTAwBITExEWloaRowYgX379mHWrFkOTkpEZDu7HiK4VGRkJCIjIwEAfn5+SE9Pb7CNr68vNm7caO9oREQtwilewRIRtUUsWCIiSRx2iIBalmfH9nBXteyv01hdi/ILVS06J9H1hAXbRrirlAh/NqNF59z6WgSc6wMzRK0LDxEQEUnCgiUikoQFS0QkCQuWiEgSFiwRkSQsWCIiSViwRESSsGCJiCRhwRIRScIzuYicFE9/bv1YsEROiqc/t348REBEJAkLlohIEhYsEZEkLFgiIklYsEREkrBgiYgkYcESEUnCgiUikoQFS0QkCQuWiEgSFiwRkSQsWCIiSViwRESSsGCJiCRhwRIRScKCJSKShAVLRCQJC5aISBKHFOzKlSsRGhqK0NBQvPLKKwCAvLw8hIeHY/jw4Vi2bJl124KCAkRFRSE4OBhxcXGora11RGQiokaze8Hm5eVh165d+PTTT7Flyxb8/PPPyMrKQmxsLFJTU5GdnY0DBw4gNzcXADBnzhwkJCRg27ZtEEIgLS3N3pGJiJrE7gWrVqsxf/58tGvXDm5ubujatSuOHj2KLl264NZbb4VSqUR4eDh0Oh1OnjwJo9GIgIAAAEBkZCR0Op29IxMRNYndv1W2e/fu1p+PHj2K7OxsTJgwAWq12jqu0WhQVFSE4uLieuNqtRpFRUWN2p+3t0ez8qrVns26vQz2zCR7XzLn59zXnv96f37bqqmZHPa13b///jseffRRzJs3D0qlEkeOHKl3vUKhgBCiwe0UCkWj9lNaaoDF0nAeW6jVntDrnetLjv8sk6wnZd2+ZM7PuRvODdjnd9pant+OdHkmFxeFzS/cHPImV35+PiZNmoRnn30Wo0ePho+PD0pKSqzXFxcXQ6PRNBjX6/XQaDSOiExE1Gh2L9jTp0/j8ccfR0pKCkJDQwEAffr0wZEjR1BYWAiz2YysrCwEBQXB19cXKpUK+fn5AIAtW7YgKCjI3pGJiJrE7ocI1q5di+rqaiQnJ1vHxo4di+TkZDz55JOorq7GoEGDoNVqAQApKSmIj49HRUUFevXqhZiYGHtHJiJqErsXbHx8POLj4694XWZmZoMxPz8/pKeny45FRNTieCYXEZEkLFgiIklYsEREkrBgiYgkYcESEUnCgiUikoQFS0QkCQuWiEgSFiwRkSQsWCIiSRy2XOH1yLNje7irmv+QX7qMnbG6FuUXqpo9JxG1PBasHbmrlAh/NqNF59z6WgSca/VMIqrDQwRERJKwYImIJGHBEhFJwoIlIpKEBUtEJAkLlohIEhYsEZEkLFgiIkl4ogHRdcizY3sA9c8KbC6eVdgQC5boOsSzCu2DhwiIiCRhwRIRScKCJSKShAVLRCQJC5aISBIWLBGRJCxYIiJJWLBERJKwYImIJOGZXJe5/IsJW+JUQp5CSNeT5n6555X+5lrr31CrKNitW7di1apVMJlMmDRpEsaPHy9tXzyFkKh5+Df0P05fsEVFRVi2bBk++eQTtGvXDmPHjkXfvn3RrVs3R0cjIroqpy/YvLw89OvXD506dQIABAcHQ6fT4YknnrDp9i4uikbvU3NT+0bfxtYcrXXuuvm7dOnSovtq7Y+LPR5zWfO31rkd4dJ9NyaHQgghZARqKatXr0ZlZSVmz54NAPjoo4/w008/YfHixQ5ORkR0dU7/KYIr9b9C4bh/yYiIbOX0Bevj44OSkhLr5eLiYmg0GgcmIiKyjdMX7IABA7Bnzx6UlZWhqqoK27dvR1BQkKNjERFdk9O/yeXj44PZs2cjJiYGJpMJDz74IO666y5HxyIiuianf5OLiKi1cvpDBERErRULlohIEhYsEZEkLFgiIklYsEREkrBgL3Hq1CmMHz8eWq0WM2bMQEVFRYNtampq8OKLL2LUqFEIDQ3Frl27HJ6pjsFgwAMPPIBvv/3W4ZmKi4sxZcoUREREYPTo0dizZ0+L59i6dStGjBiBYcOGYdOmTQ2uLygoQFRUFIKDgxEXF4fa2toWz9DYTF9++SUiIiIwcuRIzJw5E+fPn3d4pjo7d+7E0KFDpeexJdPhw4cxYcIEjBw5ElOmTHGKx+nnn39GVFQURo4ciUcffRQXLly49qSCrKZPny6ysrKEEEKsXLlSvPLKKw22efPNN8UzzzwjLBaLOHjwoBg4cKCwWCwOzVRn7ty54h//+If45ptvpOWxNdOzzz4rNm7cKIQQ4tChQ2LAgAGitra2xTKcOXNGDBkyRJw9e1ZUVFSI8PBw8fvvv9fbJjQ0VPzwww9CCCEWLFggNm3a1GL7b0qm8vJyce+994ozZ84IIYRYvny5WLx4sUMz1dHr9UKr1YohQ4ZIzWNLJovFIoYPHy5yc3OFEEK8+uqrV33e2yOTEEKMGzdO7Ny5UwghRFJSkli6dOk15+Ur2P8ymUzYu3cvgoODAQCRkZHQ6XQNtvv8888xbdo0KBQKdO/eHevWrbviegn2zAQA2dnZuOGGG9CzZ08pWRqbafjw4QgPDwcAdOnSBdXV1aisrGyxHJeustahQwfrKmt1Tp48CaPRiICAgKvmbEnXymQymfDCCy/Ax8cHANCzZ0+cPn3aoZnqxMfH27xCnexMP//8Mzp06GA9Y/Oxxx6Tuga0LZkAwGKxWP9vraqqCu7u7teclwX7X2fPnoWHhweUyosnt6nVahQVFTXYrrCwEHv37kVkZCQefvhhlJSUwMVFzsNoa6ZTp05hw4YNmDt3rpQcTck0fPhw3HjjjQCAtWvX4o477oCnZ/O/HaJOcXEx1Gq19bJGo6mX4/Lr/yxnS7pWpptuugkPPPAAAMBoNGLNmjXWy47KBADvvvsuevXqhT59+kjNYmumY8eO4eabb8a8efMQHh6OxMREdOjQwaGZAGD+/PmIi4vDwIEDkZeXh7Fjx15zXqc/VVaGzz//HElJSfXG/vKXvzTY7kqrdpnNZpw5cwYff/wxfvvtN0ydOhWff/55s8ujqZksFgvi4uKQkJBg07+o9sh0qfXr1+PDDz/Ee++916LZrvR/DZfmuNb1Mti6z/LycsycORN+fn4YPXq0QzMdPHgQ27dvx/r163HmzBmpWWzNVFtbi++++w7vvfceevfujeXLlyM5ORnJyckOy2Q0GhEXF4cNGzbgrrvuwrp16zBv3jysWbPmqvNelwUbEhKCkJCQemMmkwl9+/aF2WyGq6sr9Hr9FVftuvnmmxEaGgqFQgE/Pz/ccsstOHLkSLPXR2hqpsOHD+Pw4cOIi4sDcPFf//j4eCxevBj9+vVzSKY6r7zyCnJzc7Fp0ybccsstzcpyOR8fH+zbt896+fJV1i5fhe1qOe2VqW5sypQp6NevH2JjY6XmsSWTTqeDXq9HVFQUTCYTiouLER0djc2bNzssk1qtRpcuXdC7d28AQFhYGJ566ilpeWzJdPDgQahUKuvf+cMPP4wVK1Zcc14eIvgvNzc3BAYGIjs7GwCwZcuWK67aNWTIEOs2x48fx+nTp3H77bc7LFO3bt2Qm5uLjIwMZGRkwN/fHy+++GKzy7U5mYCLr1y//fZbvP/++y1ersC1V1nz9fWFSqVCfn7+VXPaM5PZbMZjjz2GkJAQxMXF2WVd42tleuqpp7Bt2zZkZGRgzZo10Gg0UsvVlkx33303ysrK8OuvvwIAcnJycOeddzo0U5cuXXDmzBkcPnwYALBjxw7rPwBX1bLvxbVuJ06cEI888ogICQkRkydPFufOnRNCCLF582axfPlyIcTFd4LnzJkjRowYIUaMGCFycnIcnulSjzzyiPRPEVwrk8ViEYGBgWLw4MFi5MiR1v/q3j1vKZmZmSI0NFQMHz5crFmzRgghxNSpU8VPP/0khBCioKBAREVFCa1WK5555hlRXV3dovtvbKbt27eLnj171ntMYmNjHZrpUsePH7fLpwhsyfSf//xHREVFiREjRojJkyeLkpISh2fauXOnCA8PF2FhYWLixIni2LFj15yTq2kREUnCQwRERJKwYImIJGHBEhFJwoIlIpKEBUtEJAkLlq7pxIkT6Nmz5xXPB1+wYAF69uyJsrIyBySTZ9GiRXjjjTcAANOmTcMff/wBAEhMTMTQoUOxbNky7Nq1C0OGDEFUVBSMRqMj45KTui7P5KLGU6lUOHr0KE6ePAlfX18AQGVlpfWD/G3ZW2+9Zf35ww8/xM6dO3HLLbdgwYIFeOihhzBz5kwHpiNnxoIlm7i6uiIkJARbt27FY489BgDYvn077r//frzzzjvW7XJycrBq1SqYTCa4u7tj3rx5uPvuu1FSUoLnn38epaWl0Ov18PX1xfLly+Ht7d/l+zAAAAboSURBVI2hQ4da14w9ffo0QkJCrrhwzVdffYXVq1ejpqYGZWVlGDVqFGbNmgUAWLNmDdLT03HDDTcgMDAQO3bsQE5ODmpqapCSkoK9e/fCbDajV69eiI+Ph4eHR725DQYD4uLi8Ouvv0Kj0cDV1RV///vfAQBDhw7FihUrkJSUBCEEpk2bBq1Wix07dkClUqG8vBzz5s3DqlWrsH37dlgsFvj6+iIxMRE+Pj6YMGECbrzxRhw+fBjjxo3DqFGjsGTJEhw8eBAmkwn9+/fH3LlzoVQq0bt3b0yfPh27d+9GcXExYmJiMGnSJADA6tWr8emnn0KpVKJLly5ITk6Gp6cnPvroI7z//vuwWCzo1KkTEhIS0LVrV+zbtw/JycmwWCwAgEcffdS6ChrZibzzIqitOH78uAgICBD79+8XISEh1vGJEyeK3377TfTo0UOUlpaKI0eOiLCwMFFWViaEEOLgwYPi3nvvFRUVFWL9+vVi9erVQoiL631OnTpVrF27VgghxJAhQ0RycrIQ4uK6nL17925wlozFYhGPPPKIOHLkiHW7O+64Q5SWlop///vfIjg4WJw/f15YLBaxYMEC6xlJb7zxhkhOTrau2fvaa6+JxMTEBvdxyZIlYu7cucJisYjS0lIRFBQkXn/9dWu+urN56u6rEELMmzdPvP3220IIIT799FMxa9YsYTKZhBBCfPDBB2Lq1KlCiItn1y1YsMC6r/nz54t3331XCCFEbW2teO6556xnDvXo0cO6ju7+/fuFv7+/MBqN4ssvvxTDhw+3njX30ksvidTUVPHtt9+K6OhoUVlZKYQQ4uuvv7b+jmJiYqzr9hYUFIgXXnjhar9mkoCvYMlm/v7+cHFxwYEDB+Dt7Y2Kigr06NHDen3dq666V1zAxRWJjh07hokTJ2Lfvn1Yt24djh49it9//73e8nj3338/gIuLbnh7e+P8+fO49dZb683zr3/9Czt37kRWVhYOHToEIQSqqqqQm5sLrVaLjh07AgDGjx+Pb775BsDFVfrLy8uRl5cH4OJiNd7e3g3u2549exAbGwuFQgEvLy8MGzasUY/NV199hf379yMqKgrAxVXOqqqqrNcHBgZaf965cyf279+P9PR0AGhw/LbusbjzzjtRU1ODyspK7NmzB1qt1roE5IIFCwBcXFCnsLCw3tJ558+fx7lz5xASEoJFixYhJycHAwYMwDPPPNOo+0TNx4KlRhk5ciQyMzPh5eWFiIiIetdZLBb0798fy5cvt46dPn0aGo0Gr776Kn766SdERUWhb9++qK2trbdEnEqlsv6sUCgaLB9XWVmJ0aNH44EHHkBgYCCioqLw5ZdfQggBpVJZb3tXV9d6mWJjYzFo0CAAQEVFBaqrq6943/5sDltYLBZMnToV0dHRAC5+tdClX3Ny6XqmFosFK1asQNeuXQEAFy5cqLfwS91jUTcmhICrq2u9bS5cuIALFy7AYrEgIiICc+bMsc5dXFyMG2+8EWPHjsWQIUOwe/dufP3111i5ciUyMzNbdF1eujp+ioAaJSIiAjqdDtnZ2QgLC6t3Xb9+/bB7924cOnQIAJCbm4uRI0eiuroau3btwsSJEzFq1Ch4e3sjLy8PZrPZ5v0WFhbCYDBg1qxZGDp0KL777jvU1NTAYrFg0KBB2L59O8rLywHA+soQAAYOHIhNmzZZt01ISMDSpUsbzH/fffchPT0dFosF58+fx44dOxr1uAwcOBDp6ekwGAwAgBUrVvzpAugDBw7E+vXrIYRATU0NZsyYcc31cgcMGIAvvvjCOv8bb7yB9evX495778Vnn32G4uJiAMD777+PiRMnAgDGjh2LgoICREZGYvHixbhw4YJdvtuK/oevYKlRfHx80LVrV3h6eqJTp071ruvevTsWLVqEZ555xvrKctWqVejQoQMef/xxvPLKK0hNTYWrqyv+9re/4dixYzbvt2fPnhg8eDBCQkLQsWNH3HbbbejWrRsKCwtx3333YcyYMXj44Yfh7u6O7t27o3379gCAmTNn4uWXX8bo0aNhNptxxx13YP78+Q3mf/LJJ5GYmIiQkBB4eXnVO/Rhi4ceeghFRUUYM2YMFAoFOnfu/KcLRMfFxWHJkiUIDw+HyWTCgAEDMHXq1KvOP2jQIPzxxx8YN24cgIvLVC5evBgeHh6YNm0aJk+eDIVCAQ8PD6xcuRIKhQLPPfccXnrpJSxfvhwuLi544okn8H//93+Nul/UPFxNi1q9/fv344cffkBMTAwAYN26dfjxxx/rHaogcgQWLLV6BoMBsbGxOHz4sPXV4+LFi61fLkjkKCxYIiJJ+CYXEZEkLFgiIklYsEREkrBgiYgkYcESEUny/ynntnALFWHhAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(figsize=(5, 5))\n",
"ax.hist(permutations, bins=11, rwidth=0.9)\n",
"ax.axvline(x = df_sub_F['age'].mean() - df_sub_M['age'].mean(), color='black', lw=2)\n",
"ax.set_xlabel('Mean age differences')\n",
"ax.set_ylabel('Frequency')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this figure we compare the observed difference with distribution of the permuted differences. The black line shows the difference between the means we computed earlier. You'll notice that the oberved values is well within the set of permuated values."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Oftentimes, statistical analysis requires you to set a threshold value for rejecting the null hypothesis. As we're interested in the deviation from the general pattern—women being younger than men—we calculate how likely we encounter values equal to or smaller than the observed -0.2 under the null hypothesis. \n",
"\n",
"\n",
"We can compute the probability of this scenario with values generated by the permutation test. Below we ask how likely are the compute differences equal to or smaller than -0.2."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.1282"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len([i for i in permutations if i <= (df_sub_F['age'].mean() - df_sub_M['age'].mean())]) / len(permutations)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A typical threshold value for rejecting the null hypothesis is 0.05 but this is a pure convention as nothing magically happens around this 0.05. In the social sciences, you'll often encounter 0.1 as the threshold, while in physics thresholds smaller than 0.00001 are common. Again this is a convention and changes by discipline or study.\n",
"\n",
"In any case, you'll notice that the probability of the null hypothesis being true is around 0.13 (value may differ because of randomization, substantially above 0.05 and therefore we can not accept the difference as statistically significant."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The result of the permuation test is largely similar to the Student's t-test. The p-value in this case is the probability of observed values equal to or lower than the observed difference, assuming there is no difference between men amd womenin the sample (a one-side t-test)."
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Ttest_indResult(statistic=-1.151420997318941, pvalue=0.12478303628250077)"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from scipy.stats import ttest_ind\n",
"ttest_ind(df_sub_F['age'],df_sub_M['age'],alternative='less')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can repeat the permutation for Whitechapel. In this case, the age difference should be highly significant."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Gender and Disability\n",
"### An additional example\n",
"\n",
"We can also use permutation to test the relation between categorical variables: in this case we study the relation gender and disability as reported in our census sample."
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"p = df.sample(frac=.1).groupby('gender')['disability'].agg([np.sum, len])"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
"p['no_disability'] = p['len'] - p['sum']"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"p.rename({'sum':'disability','len':'all'},axis=1,inplace=True)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
![](\"https://www.ghil.ac.uk/typo3conf/ext/wacon_ghil/Resources/Public/Images/institute_icon_small.png\"){kind=icon}
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We've covered how to describe, summarize, and compare variables. However, we lack a formal procedure to assess the importance of the differences we observe. For example, we established that men are, on average, one year younger than women. But how can we establish the value or 'significance' of this gap?\n",
"\n",
"In the notebook, we move from descriptive to inferential statistics and compute the extent to which means between subgroups in our data are statistically significant."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For more background on the concepts and terminology used in this notebook, please consult the lecture by Luke Blaxill."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We repeat some of the code from the previous notebook\n",
"- load the required libraries\n",
"- load the synthetic census data"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd\n",
"import numpy as np\n",
"import seaborn as sns\n",
"import random\n",
"from scipy import stats\n",
"from tqdm.notebook import tqdm\n",
"sns.set()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"df = pd.read_csv('data/democraphic_data/london_subsample.csv',index_col=0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To make our analysis more interesting and complex, we add a dimension: place. We study how age differences between men and women vary depending on the registration district. In Pandas, adding subgroups is convenient: simply pass a list with column names `['district', 'gender']` (instead of just one column as we've done previously). `.groupby()` will read this list from left to right, i.e. groups the data first by place and then further aggregates by `gender` within each district."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"district gender\n",
"Bethnal Green F 25.920963\n",
" M 25.613501\n",
" U 20.766603\n",
"Camberwell F 27.986109\n",
" M 26.708931\n",
" ... \n",
"Whitechapel M 26.254569\n",
" U 24.051316\n",
"Woolwich F 26.653639\n",
" M 26.311769\n",
" U 21.427165\n",
"Name: age, Length: 90, dtype: float64"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"by_reg_gen = df.groupby(['district','gender'])['age'].mean()\n",
"by_reg_gen"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Working with the output of this operation requires a bit more thought. The `.groupby()` arranges data slightly differently depending on whether you group on one column or more. This becomes apparent when printing the `type()` of the `.index` attribute."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pandas.core.indexes.base.Index"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"type(df.groupby('gender')['age'].median().index)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pandas.core.indexes.multi.MultiIndex"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"type(by_reg_gen.index)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`by_reg_gen` orders the data using a `MultiIndex`, which means that the index contains multiple levels (place and gender) via which we can access our data.\n",
"\n",
"Place sits the highest level of our grouped data frame. We can access the separate districts using `.loc[]`"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"gender\n",
"F 25.920963\n",
"M 25.613501\n",
"U 20.766603\n",
"Name: age, dtype: float64"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"by_reg_gen.loc['Bethnal Green']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Similarly, we can slice the data by place."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"district gender\n",
"Bethnal Green F 25.920963\n",
" M 25.613501\n",
" U 20.766603\n",
"Camberwell F 27.986109\n",
" M 26.708931\n",
" U 23.410050\n",
"Chelsea F 31.318046\n",
" M 30.430615\n",
" U 28.800000\n",
"Name: age, dtype: float64"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"by_reg_gen.loc['Bethnal Green':'Chelsea']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But our index has two levels, place **and** gender. We can obtain the means for women for all districts using the following syntax."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"district\n",
"Bethnal Green 25.920963\n",
"Camberwell 27.986109\n",
"Chelsea 31.318046\n",
"Fulham 27.920645\n",
"Greenwich 27.593689\n",
"Hackney 28.796820\n",
"Hampstead 29.742615\n",
"Holborn 27.719117\n",
"Islington 29.047747\n",
"Kensington 30.823465\n",
"Lambeth 29.019583\n",
"Lewisham 28.879865\n",
"London City 31.323944\n",
"Marylebone 30.465945\n",
"Mile End Old Town 26.298023\n",
"Paddington 30.450756\n",
"Pancras 29.136537\n",
"Poplar 26.384635\n",
"Shoreditch 27.056387\n",
"Southwark 27.144232\n",
"St George Hanover Square 30.821499\n",
"St George In The East 24.545266\n",
"St Giles 30.129131\n",
"St Olave Southwark 26.486323\n",
"Stepney 25.985841\n",
"Strand 29.288732\n",
"Wandsworth 28.028524\n",
"Westminster 28.990447\n",
"Whitechapel 24.649451\n",
"Woolwich 26.653639\n",
"Name: age, dtype: float64"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"by_reg_gen.loc[:,'F']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice the comma in `.loc[]` (`by_reg_gen.loc[:**,**'F']`) \n",
"- the part before the comma indicates the items we want to access from the first level (place). In this case, we entered a colon, meaning from the first till the last rows. \n",
"- the part after the comma indicates the items we want to select from the second level. In this case, we want to retrieve all elements with value 'F' for the second level of the MultiIndex..\n",
"\n",
"Computing the age differences by place is then fairly straightforward."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"district\n",
"Bethnal Green 0.307461\n",
"Camberwell 1.277178\n",
"Chelsea 0.887431\n",
"Fulham 1.295697\n",
"Greenwich 0.971892\n",
"Hackney 1.357030\n",
"Hampstead 0.645606\n",
"Holborn 0.341140\n",
"Islington 1.430837\n",
"Kensington 1.410332\n",
"Lambeth 1.411260\n",
"Lewisham 1.420661\n",
"London City 3.262573\n",
"Marylebone 0.921905\n",
"Mile End Old Town 0.599087\n",
"Paddington 1.664367\n",
"Pancras 1.031302\n",
"Poplar -0.218534\n",
"Shoreditch 0.130194\n",
"Southwark 0.101532\n",
"St George Hanover Square 1.412338\n",
"St George In The East -0.011402\n",
"St Giles -0.879827\n",
"St Olave Southwark 0.088463\n",
"Stepney -0.334267\n",
"Strand -1.759558\n",
"Wandsworth 1.004169\n",
"Westminster -0.329065\n",
"Whitechapel -1.605118\n",
"Woolwich 0.341871\n",
"Name: age, dtype: float64"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f_m_diff = by_reg_gen.loc[:,'F'] - by_reg_gen.loc[:,'M']\n",
"f_m_diff"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Of course, you can be more restrictive and slice the data by place and obtain only means for women."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"district gender\n",
"Bethnal Green F 25.920963\n",
"Camberwell F 27.986109\n",
"Chelsea F 31.318046\n",
"Name: age, dtype: float64"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"by_reg_gen.loc['Bethnal Green':'Chelsea','F']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Hypothesis testing"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At this stage, we can compute and compare the distribution of variables, calculate their means or other relevant statistics. But a question then immediately appears: are the differences we observe \"significant\", and what do we mean with \"significance\" anyway. In this section, we have a closer look at hypothesis testing from a data-driven perspective.\n",
"\n",
"Traditionally, statistical methods, such as the Student's t-test arose in times of limited computing power and relied on equations and assumptions related to the distribution of the data. Explaining this requires many detours and implies a steep learning curve for the statistically uninitiated. \n",
"\n",
"In this lecture, we explore a data-driven, and hopefully intuitive procedure, for assessing the significance of observed statistics (such as the mean). First, we discuss a procedure called **bootstrapping** and then explain how to use **permutation** for significance testing. \n",
"\n",
"The question we first address regards the relation between `gender` and `age`, and we went to assess if variations we observe over the districts is significant. We focus on Whitechapel and Poplar as we noticed in the previous section that both these districts deviate from the general pattern and contain a slightly younger female population.\n",
"\n",
"The age gap is bigger in Westminster than in Poplar, but is the differences of the mean significant?"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"((18772, 4), (40838, 4))"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_whitechapel = df[df.district=='Whitechapel']\n",
"df_poplar = df[df.district=='Poplar']\n",
"df_whitechapel.shape, df_poplar.shape"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"