{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sampling and Testing the Penguins\n",
"\n",
"This notebook uses the [Palmer Penguins data](https://github.com/allisonhorst/palmerpenguins) to demonstrate confidence intervals and two-sample hypothesis tests, both using parametric methods and the bootstrap.\n",
"\n",
"Our question: What are the average values of various penguin dimensions (bill length and depth, flipper length, and body mass)? Do these dimensions differ between penguin species?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"Python modules:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import scipy.stats as sps\n",
"import statsmodels.api as sm\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Set up a random number generator:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"rng = np.random.default_rng(20200913)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Load the Penguin data:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
species
\n",
"
island
\n",
"
bill_length_mm
\n",
"
bill_depth_mm
\n",
"
flipper_length_mm
\n",
"
body_mass_g
\n",
"
sex
\n",
"
year
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
Adelie
\n",
"
Torgersen
\n",
"
39.1
\n",
"
18.7
\n",
"
181.0
\n",
"
3750.0
\n",
"
male
\n",
"
2007
\n",
"
\n",
"
\n",
"
1
\n",
"
Adelie
\n",
"
Torgersen
\n",
"
39.5
\n",
"
17.4
\n",
"
186.0
\n",
"
3800.0
\n",
"
female
\n",
"
2007
\n",
"
\n",
"
\n",
"
2
\n",
"
Adelie
\n",
"
Torgersen
\n",
"
40.3
\n",
"
18.0
\n",
"
195.0
\n",
"
3250.0
\n",
"
female
\n",
"
2007
\n",
"
\n",
"
\n",
"
3
\n",
"
Adelie
\n",
"
Torgersen
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
2007
\n",
"
\n",
"
\n",
"
4
\n",
"
Adelie
\n",
"
Torgersen
\n",
"
36.7
\n",
"
19.3
\n",
"
193.0
\n",
"
3450.0
\n",
"
female
\n",
"
2007
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" species island bill_length_mm bill_depth_mm flipper_length_mm \\\n",
"0 Adelie Torgersen 39.1 18.7 181.0 \n",
"1 Adelie Torgersen 39.5 17.4 186.0 \n",
"2 Adelie Torgersen 40.3 18.0 195.0 \n",
"3 Adelie Torgersen NaN NaN NaN \n",
"4 Adelie Torgersen 36.7 19.3 193.0 \n",
"\n",
" body_mass_g sex year \n",
"0 3750.0 male 2007 \n",
"1 3800.0 female 2007 \n",
"2 3250.0 female 2007 \n",
"3 NaN NaN 2007 \n",
"4 3450.0 female 2007 "
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"penguins = pd.read_csv('../data/penguins.csv')\n",
"penguins.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Some of these names are cumbersome to deal with. I'm going to give them shorter names:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"penguins.rename(columns={\n",
" 'bill_length_mm': 'BillLength',\n",
" 'bill_depth_mm': 'BillDepth',\n",
" 'flipper_length_mm': 'FlipperLength',\n",
" 'body_mass_g': 'Mass'\n",
"}, inplace=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A few things will be eaiser if we also split out penguins by species:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
species
\n",
"
island
\n",
"
BillLength
\n",
"
BillDepth
\n",
"
FlipperLength
\n",
"
Mass
\n",
"
sex
\n",
"
year
\n",
"
\n",
" \n",
" \n",
"
\n",
"
276
\n",
"
Chinstrap
\n",
"
Dream
\n",
"
46.5
\n",
"
17.9
\n",
"
192.0
\n",
"
3500.0
\n",
"
female
\n",
"
2007
\n",
"
\n",
"
\n",
"
277
\n",
"
Chinstrap
\n",
"
Dream
\n",
"
50.0
\n",
"
19.5
\n",
"
196.0
\n",
"
3900.0
\n",
"
male
\n",
"
2007
\n",
"
\n",
"
\n",
"
278
\n",
"
Chinstrap
\n",
"
Dream
\n",
"
51.3
\n",
"
19.2
\n",
"
193.0
\n",
"
3650.0
\n",
"
male
\n",
"
2007
\n",
"
\n",
"
\n",
"
279
\n",
"
Chinstrap
\n",
"
Dream
\n",
"
45.4
\n",
"
18.7
\n",
"
188.0
\n",
"
3525.0
\n",
"
female
\n",
"
2007
\n",
"
\n",
"
\n",
"
280
\n",
"
Chinstrap
\n",
"
Dream
\n",
"
52.7
\n",
"
19.8
\n",
"
197.0
\n",
"
3725.0
\n",
"
male
\n",
"
2007
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" species island BillLength BillDepth FlipperLength Mass sex \\\n",
"276 Chinstrap Dream 46.5 17.9 192.0 3500.0 female \n",
"277 Chinstrap Dream 50.0 19.5 196.0 3900.0 male \n",
"278 Chinstrap Dream 51.3 19.2 193.0 3650.0 male \n",
"279 Chinstrap Dream 45.4 18.7 188.0 3525.0 female \n",
"280 Chinstrap Dream 52.7 19.8 197.0 3725.0 male \n",
"\n",
" year \n",
"276 2007 \n",
"277 2007 \n",
"278 2007 \n",
"279 2007 \n",
"280 2007 "
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"chinstraps = penguins[penguins['species'] == 'Chinstrap']\n",
"chinstraps.head()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"adelies = penguins[penguins['species'] == 'Adelie']\n",
"gentoos = penguins[penguins['species'] == 'Gentoo']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That's all we need right now!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Confidence Intervals\n",
"\n",
"Remember that $s/\\sqrt{n}$ is the **standard error** of the mean. The SciPy [`sem`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.sem.html) function computes the standard error of the mean. We can multiply this by 1.96 to get the distance between the sample mean and the edge of the confidence interval.\n",
"\n",
"So let's write a function that returns the 95% CI of the mean of some data:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def mean_ci95(xs):\n",
" mean = np.mean(xs)\n",
" err = sps.sem(xs)\n",
" width = 1.96 * err\n",
" return mean - width, mean + width"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Warmup: Chinstrap Flippers\n",
"\n",
"As a warmup, let's compute confidence intervals for the chinstrap flipper length. We will do this with both the standard error and with the bootstrap.\n",
"\n",
"Let's get the mean & SD of the Chinstrap penguins:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(195.8235294117647, 7.131894258578147, 68)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p_mean = chinstraps['FlipperLength'].mean()\n",
"p_std = chinstraps['FlipperLength'].std() # defaults to sample SD\n",
"p_n = len(chinstraps)\n",
"p_mean, p_std, p_n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What's the confidence interval?"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(194.12838574870136, 197.51867307482803)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mean_ci95(chinstraps['FlipperLength'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's bootstrap the chinstraps, and compute the percentile 95% confidence interval for them:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([194.14705882, 197.5 ])"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"boot_means = [np.mean(rng.choice(chinstraps['FlipperLength'], p_n)) for i in range(10000)]\n",
"np.quantile(boot_means, [0.025, 0.975])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's break that down a bit. It is a Python *list comprehension*, a convenient way of assembling a list. It is equivalent to the following code:\n",
"\n",
"```python\n",
"boot_means = []\n",
"for i in range(10000):\n",
" # compute the bootstrap sample\n",
" boot_samp = rng.choice(chinstraps['FlipperLength'], p_n)\n",
" # compute the mean of the bootstrap sample\n",
" boot_mean = np.mean(boot_samp)\n",
" # add it to our list of bootstrap means\n",
" boot_means.append(boot_mean)\n",
"```\n",
"\n",
"It just does all of that in one convenient line, without excess variables. We then pass it to `np.quantile` to compute the confidence interval; numpy functions can take lists of numbers as well as arrays."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's see the distribution of those bootstrap sample means, with a vertical line indicating the observed sample mean:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAH9CAYAAADVvj8ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABfAUlEQVR4nO3deXxTVd4/8E/2NE2T7mkLbdn3HQXqLiCIiAiMOoiKuOCDuPLoKM9PxdFxUGbGbQZ1cBxwHRVHXFBRBARBUCgFWvaWpaX7njZt9vv7I02ktoUuae5N8nm/XnlJktvke2ybT8+5554jEwRBABEREUmSXOwCiIiIqG0MaiIiIgljUBMREUkYg5qIiEjCGNREREQSxqAmIiKSMAY1ERGRhDGoAQiCALPZDF5STkREUsOgBlBXVwej0Yi6ujqxSyEKSxaLBTKZDDKZDBaLRexyiCSFQU1ERCRhDGoiIiIJY1ATERFJGIOaiIhIwhjUREREEsagJiIikjAGNRERkYQxqImIiCSMQU1ERCRhDGoiIiIJY1ATERFJGIOaiIhIwhjUREREEsagJiIikjAGNRERkYQxqImIiCSMQU1ERCRhDGoiIiIJY1ATERFJGIOaiIhIwhjUREREEqYUuwAiorPNueUOKJW/fjQlxBixetVrIlZEJC4GNRFJyhULlkITofPd37TqjyJWQyQ+Dn0TERFJGIOaiIhIwhjUREREEsagJiIikjAGNRERkYQxqImIiCSMQU1ERCRhDGoiIiIJY1ATERFJGIOaiIhIwhjUREREEsagJiIikjAGNRERkYQxqImIiCSMQU1ERCRhDGoiIiIJY1ATERFJGIOaiIhIwhjUREREEsagJiIikjAGNRERkYQxqImIiCSMQU1ERCRhDGoiIiIJU4pdABEFjwUL70V5dW2LxxNijFi96jURKmopGGok6ggGNRG1W3l1LSYtXNbi8U2r/ihCNa0LhhqJOkLUoe+nn34aMpms2W3QoEG+561WKxYvXoy4uDjo9XrMmTMHpaWlzV4jPz8f06dPh06nQ2JiIh599FE4nc5AN4WIiKhbiN6jHjp0KL7//nvffaXy15IefvhhfPXVV1i7di2MRiPuu+8+zJ49Gzt27AAAuFwuTJ8+HUlJSfjpp59QXFyM2267DSqVCn/+858D3hYiIiJ/Ez2olUolkpKSWjxeW1uLt956Cx988AEmTpwIAFi9ejUGDx6MXbt2YcKECfjuu+9w6NAhfP/99zCZTBg1ahSeffZZPPbYY3j66aehVqsD3RwiIiK/En3W9/Hjx5GSkoI+ffpg3rx5yM/PBwBkZmbC4XBg8uTJvmMHDRqEtLQ07Ny5EwCwc+dODB8+HCaTyXfM1KlTYTabcfDgwTbf02azwWw2N7sRERFJkahBPX78eKxZswYbNmzA66+/jpMnT+LSSy9FXV0dSkpKoFarER0d3exrTCYTSkpKAAAlJSXNQtr7vPe5tixfvhxGo9F3S01N9W/DiIiI/ETUoe9p06b5/j1ixAiMHz8e6enp+PjjjxEREdFt77t06VIsWbLEd99sNjOsiYhIkkQf+j5bdHQ0BgwYgNzcXCQlJcFut6OmpqbZMaWlpb5z2klJSS1mgXvvt3be20uj0cBgMDS7ERERSZGkgrq+vh55eXlITk7G2LFjoVKpsGnTJt/zR48eRX5+PjIyMgAAGRkZyM7ORllZme+YjRs3wmAwYMiQIQGvn4iIyN9EHfp+5JFHMGPGDKSnp6OoqAjLli2DQqHA3LlzYTQaceedd2LJkiWIjY2FwWDA/fffj4yMDEyYMAEAMGXKFAwZMgS33norVqxYgZKSEjzxxBNYvHgxNBqNmE0jIiLyC1GD+syZM5g7dy4qKyuRkJCASy65BLt27UJCQgIA4KWXXoJcLsecOXNgs9kwdepUvPbar0sAKhQKrF+/HosWLUJGRgYiIyMxf/58PPPMM2I1iYiIyK9EDeoPP/zwnM9rtVqsXLkSK1eubPOY9PR0fP311/4ujYiISBJEX/CEiOhsv5ysQqNghlophzFCBadCK3ZJRKJiUBORqOxON/7z82nf/T2nqyFXnxXO6ddjyUf7sGzGUBh1KhEqJBIXg5qIRJNbVocHP9yH7FO/XrkxwKRHYqwRDpcbBVWNKDFb8WlWIX4+WYW/3zwaY9JiRKyYKPAkdXkWEYWPrw4U49q/b8fBInOznvLkwSZc2CsWF/WNx00XpqLPmW+RHqdDYU0j5r35MzJPV4tYNVHgMaiJKOD+vf0k7vvPXlgdblzaPx6fL764zWNzf9kM3e41iGwoRqPDhRv+8QMmz7sXCxbeG8CKicTDoW8iCqg3t53Ac18fBgDMz0jHUzOGwtrY0ObxDreAKXf/PzhcbqzLKkRxLVDWfyZUhz4IVMlEomKPmogC5qPd+b6QXnLVADx93VAo5LJ2fa1KIcd1I1Ng0CphtjpRnHBhd5ZKJBkMaiIKiO3HK7D002wAwD2X98EDk/pDJmtfSHtpVQpMHZoEGYCaqN74cn9RN1RKJC0MaiLqdoU1jXjgwyy4BWD2mB54/OpBnX6tlOgIXNg7FgDwzPpDqLM6/FUmkSQxqImoWzlcbtz7/l5UWewYmmLAn2cN73BP+rcu7BUDtd2M8job/rE510+VEkkTg5qIutWqbSewv6AGBq0Sb9wyFlqVosuvqZTLkVy5FwDw7x0ncaK8vsuvSSRVDGoi6ja5ZfV4ZdNxAMCyGUORGqvz22tHNRThyoEJcLgE/G3jMb+9LpHUMKiJqFsIgoClnx6A3enG5QMSMHtMD7+/x2PTPOe6v84uxrHSOr+/PpEUMKiJqFt8e7AEu09VI0KlwJ9nd/28dGsGJRkwbVgSBAF4tannThRqGNRE5HcOlxsrNhwFANx9aW/0iI7otvd6YFJ/AMBX2cU4zl41hSAGNRH53Ye7C3CiwoK4SDUWXt63W99rcLIBU4eaIAieiWVEoYZBTUR+ZXe6sbLpkqkHJ/eHXtP9KxXfdWkfAMC6rEI45epufz+iQGJQE5Fffb6vECVmK0wGDX5/YVpA3vOC9BgMTTHA6nCj2tC9PXiiQOOmHETkN263gFXbTgAA7ri4N+659z6UV9e2OC4hxojVq17z2/vKZDLcflEvPPrJAVQZBsDtFiBv5xriRFLHoCYiv9lytAzHy+oRpVFi7vg0fPmPWkxauKzFcZtW/dHv7z1jZAqe/+YIKi3AyUoL+ibo/f4eRGJgUBOR33gnc908Pg0GrarN4w7s349rb5jnu+90Orv83lqVAnPG9sSqbSdwuNjMoKaQwaAmIr84XWnBjtxKyGTArRnp5zzW4Raa9bRtjQ349rOPu1zD75qC+mSFBQ12J3RqfsRR8ONkMiLyiw93FwAALuufgJ4x/lsqtCMGmKIQYa2EWwCOlvCaagoNDGoi6jI35Fi7xxPUc8cFZqZ3W6LrPJPZDhWbRa2DyF8Y1ETUZXWRPVBRb0e8XoNJgxNFrcVYfxoKmQwV9XZU1NtErYXIHxjURNRlNVG9AQA3XNATKoW4HytKtx3pcZ6h9+Nl3P6Sgh+Dmoi6xOpwoV6XDACYPdr/O2R1Rr9Ez4zvXAY1hQAGNRF1SW5ZPQSZAoOTDehvihK7HABAn/hIyGVAlcWOSg5/U5BjUBNRl3hnV88clSJyJb/SqBRIi/UMf+eWs1dNwY1BTUSdVmd14ExNIwDPymBS4h3+5nlqCnYMaiLqNG8I6hrLunXP6c7om6CHTAZU1tthbnSIXQ5RpzGoiajTvJO1jPX5IlfSklalQLJRCwA4VWkRuRqizmNQE1GnWGxOFNdaAQAGS4HI1bSuV1wkAOBUZYPIlRB1HoOaiDrFO0kryaCFytUocjWt8wZ1QVUD3DJ+3FFw4k8uEXVKXtOwt3fSlhTF69XQa5RwugVYtOKumEbUWdxahog6rNHh8s327peoxx4RavjtVple2Tk5mNT0b5lMhvQ4HQ4WmVGvk9asdKL2YlATUYedrLBAEDw9VmNE2/tOd6ffbpXplbloVrP7veIicbDIjDoGNQUpDn0TUYedaDo/3TdBusPeXqmxEZABsKsNKK6V5rl0onNhUBNRhzhdbuRXeWZR946PFLma89MoFUg0aAAAO/MqRa6GqOMY1ETUIWeqG+FwCYjUKJAYpRG7nHbpGeNZTpRBTcGIQU1EHXKiwrN4SJ94PWQymcjVtE9qjGfVtJ0nGNQUfBjURNRuAjwTyQDPDlXBItkYAQhunKluREEVFz+h4MKgJqJ2s6pjUG9zQqWQoWeMtNb2Phe1Ug6dtQIAh78p+DCoiajdzJE9AQBpsTooFcH18RHZWAYA+CmvQuRKiDomuH7TiEhUdZE9AHjOTwebSGspAGD3qWqRKyHqGAY1EbVLUU0jrJpYAECveJ3I1XRchLUSCrkMhTWNvJ6aggpXJiOidtl02NMjTTZqoVM3/+hoz3KeYlMITgxOjkJOoRl7TlVjxsjgOcdO4Y1BTUTt8v1hzznePgktZ3u3dzlPsV2QHoucQjMyT1djxkguKUrBgUFNROdVb3P6ZksH4/lpwNPrr418D0i6BB9u2Ys97/4ZAJAQY8TqVa+JXB1R2xjURGFgwcJ7UV5d2+Lx9obU1qPlsLvcUNvrEKMTZxOOrnK4BVw15xb8e8cp2LSxuPSOJ6FWyrFp1R/FLo3onBjURGGgvLq21aHp9obUd4dKAAAGSwFksjF+rS2QorQqRGmVqLM6UWK2Ii02+CbFUfjhrG8iOie7043NRzznpw2WMyJX03XJRi0AoLiGM78pODCoieicfj5ZiTqrE/F6DSJswb+qV4rRM9u72GwVuRKi9mFQE9E5fXfQc1nWVUMSIYMgcjVdZzJ4etSlZisEIfjbQ6GPQU1EbXK7BWw85AnqKUOSRK7GP+Kj1JDLAKvDDbPVKXY5ROfFoCaiNmUX1qLEbEWkWoGMvnFil+MXSrkc8XrPPtqlHP6mIMCgJqI2eWd7XzEwEVqVQuRq/Ofs4W8iqWNQE1GbfMPeQ00iV+JfJoO3R20TuRKi82NQE1GrTlZYcKy0Hkq5DFcMTBS7HL9KaupRl9VZIUAmcjVE58agJqJWbWwa9p7QJw7GiOBcjawtMZFqqBQyOFwCbCqD2OUQnRODmoha9e3B0Bz2BgC5TIbEKE+vulEbGpPkKHQxqImohcKaRmSeroZMFjqXZf2W9zx1Y9Me20RSxaAmoha+OlAEABjXKxZJTUtuhhrvzO9GDXvUJG0MaiJqYf2BYgDAtSG8Z7M3qK2aaNicLpGrIWobg5qImjlVYcGBM7WQy4Bpw0Jz2BsADFoltCo5BJkCh4vrxC6HqE0MaiJq5qtsT2/64n7xvhW8QpFMJvP1qg+cqRG3GKJzYFATUTNf7vecn752RLLIlXQ/b1DvL6gVuRKitjGoicjneGkdjpTUQaWQYerQ0B329vLO/N7PHjVJmGSC+vnnn4dMJsNDDz3ke8xqtWLx4sWIi4uDXq/HnDlzUFpa2uzr8vPzMX36dOh0OiQmJuLRRx+F08kdcYg648umSWSX9k9AtE4tcjXdz9R0LXVeeT3qbfzcIGmSRFDv3r0b//znPzFixIhmjz/88MP48ssvsXbtWmzduhVFRUWYPXu273mXy4Xp06fDbrfjp59+wttvv401a9bgqaeeCnQTiIKeIAhY33RZ1oyRoT/sDQCRGiVUDgsEAcg+w+FvkibRg7q+vh7z5s3Dm2++iZiYGN/jtbW1eOutt/Diiy9i4sSJGDt2LFavXo2ffvoJu3btAgB89913OHToEN577z2MGjUK06ZNw7PPPouVK1fCbre3+Z42mw1ms7nZjSjcHSo240S5BWqlHJMHh95qZG2JsFUC4PA3SZfoQb148WJMnz4dkydPbvZ4ZmYmHA5Hs8cHDRqEtLQ07Ny5EwCwc+dODB8+HCbTrx8qU6dOhdlsxsGDB9t8z+XLl8NoNPpuqampfm4VUfDxXjs9cWAiorShtbb3uXiDOruQPWqSJlGD+sMPP8TevXuxfPnyFs+VlJRArVYjOjq62eMmkwklJSW+Y84Oae/z3ufasnTpUtTW1vpuBQUFXWwJUXATBOHX2d5hMuztFWGrAsChb5IupVhvXFBQgAcffBAbN26EVhvYJQo1Gg00mtC9PpSoo/YV1OBMdSN0agUmDgqtLS3PxxvU+VUNqGmwh8UkOgouovWoMzMzUVZWhjFjxkCpVEKpVGLr1q149dVXoVQqYTKZYLfbUVNT0+zrSktLkZTkuWwkKSmpxSxw733vMUR0ft5h78mDTdCpRfv7XRQKtwO94nQAOPxN0iRaUE+aNAnZ2dnYt2+f73bBBRdg3rx5vn+rVCps2rTJ9zVHjx5Ffn4+MjIyAAAZGRnIzs5GWVmZ75iNGzfCYDBgyJAhAW8TUTByu3+d7R0Oi5y0ZnjPaADAAQ5/kwSJ9qdzVFQUhg0b1uyxyMhIxMXF+R6/8847sWTJEsTGxsJgMOD+++9HRkYGJkyYAACYMmUKhgwZgltvvRUrVqxASUkJnnjiCSxevJhD20TttPtUFUrNNkRplbh8YILY5YhiRA8jvtxfxPPUJEmSHuN66aWXIJfLMWfOHNhsNkydOhWvvfaa73mFQoH169dj0aJFyMjIQGRkJObPn49nnnlGxKqJgot32Hvq0CRolAqRqxHH8J5GABz6JmmSVFD/8MMPze5rtVqsXLkSK1eubPNr0tPT8fXXX3dzZUShSYAMXzdtwhGuw94AMDTFAJkMKKxpREW9LaQ3I6HgI/p11EQkHkuECZUWO2J0KlzcL17sckQTpVWhT3wkAPaqSXoY1ERhrFafBgCYNjwZKkV4fxyM8E4o405aJDHh/ZtJFMZcbgG1kZ6gDudhb6/hPbznqWvELYToNxjURGEqv6oBboUaCVEajO8dJ3Y5ohvRNKGMl2iR1EhqMhkRBc6JinoAwJQhJijksmbPLVh4L8qrWwZWdk4OJgWkusAbkmKAXAaU1dlQarbCZAjsiolEbWFQE4UhQRBwssICAJg8pOVOWeXVtZi0cFmLxzMXzer22sSiUysxwBSFIyV1OHCmFlcNYVCTNHDomygMldXZYLG5IHc7kNGHw95evvPU3PKSJIRBTRSGTjT1pvUNxdCqwnORk9b4zlPzEi2SEAY1URg6We4J6qiGQpErkRbvmt/ZZ2ohCIK4xRA1YVAThZk6qwPl9TYAQJSlSORqpGVQUhSUchkqLXYU1VrFLocIAIOaKOzkVzUAAEwGDZRum8jVSItWpcDApCgAPE9N0sGgJgozBVWNAIC0WJ3IlUgTr6cmqWFQE4URQRBQUO3pUTOoWze8RzQArvlN0sGgJgojlRY7GuwuKOUyJBl5nXBrzu5Rc0IZSQEXPCEKIwVN56d7REdAKZfjwP79uPaGeS2OC+UVyM5ngCkKaoUctY0O5Fc1ID0uUuySKMwxqInCiHciWWrTsLfDLYTdCmS/1dofK4oeUwFtHA6cqWVQk+gY1ERhwuUWUFjDiWS/1dofK7IjZcgurEV2YS1mjEwRqTIiD56jJgoTFfU2OFwCNEo54vVqscuRNJNBAwA4wEu0SALYoyYKE0VNvelkoxYymew8R4e3xCjPRLucQjPcbgHys3YXa2tnsYQYI1avei1gNVL4YFAThYnippW2kqMjRK5E+uIi1ZC5nai3AXnl9ehvivI919bOYptW/TGQJVIY4dA3URgQ8GtQp/CyrPOSy2WIsFUCALLya8QthsIeg5ooDDiUOtTbnJDJAJOBQd0eOmsFACCroFrkSijcMaiJwkCDNgEAkKDXQKXgr3176KzsUZM08DeWKAw0aOMBACk8P91uETZPj/poaR3qbU6Rq6FwxqAmCgPeHjXPT7efymVFj+gICAJwoKBG7HIojDGoiUKc1eGCVR0NAFzfu4NGp0UDALIY1CQiBjVRiDtcbAZkckSoFNBreEVmR4xOiwEAZOVzQhmJh0FNFOJymrZrTDRouNBJB/l61Pk13EmLRMOgJgpxB854gtoUxWHvjhqaYoBaIUelxY6Cqkaxy6EwxaAmCnHZZ/WoqWM0SgWGpBgA8HpqEg+DmiiEWR0uHC+rBwAkRjGoO+Ps4W8iMTCoiULY4WIzXG4BCpeVE8k6iRPKSGwMaqIQ5p1IFmGr4kSyThqdGg0AOFhkhtXhErcYCksMaqIQ5p1IFmGrErmS4NUzJgLxeg2cbgEHi1pub0nU3RjURCHsYJEZAIO6K2QyGc9Tk6gY1EQhyuFyI7dpIpnWxvOrXeEN6r08T00iYFAThaiTFRbYXW5EqhVQOS1ilxPURqd6JpTtPV0jbiEUlhjURCHqcLFn2HtgUhQ4jaxrRqYaoZDLUGK2wq7QiV0OhRkGNVGIOlpSBwAYmGQQuZLgp1MrMSgpCgDQ2LRlKFGgMKiJQtSRpqAenBwlciWhYUzT9dQNDGoKMAY1UYjy9qgHsUftF2PTGdQkDgY1UQiqbXSgsMazicTAJPao/cHbo7ZqYuB0uUWuhsIJg5ooBHl70ylGLYwRKpGrCQ2psRGI16shyBQoq7OJXQ6FEQY1UQg6WuKZ8T0omcPe/uJZ+MTTqy6ptYpcDYUTBjVRCDrsOz/NYW9/8g5/FzOoKYAY1EQh6Hip99IsBrU/jWlaoazY3AhBEMQthsIGg5ooxAiC4NuDul+iXuRqQsuIntGA4IbF5kKd1Sl2ORQmGNREIabSYkdNgwMyGdA3gUHtTxFqhW/ddA5/U6AwqIlCzPFST286NUYHrUohcjWhR2etAMAJZRQ4SrELICL/WbDwXhxxJQAJ41BdcBzX3vAvAEB2Tg4miVxbqNDZKlCFgSg2N4pdCoUJBjVRCCmvrkXi5Tej+EwtBgwehkuuuwIAkLlolriFhRBvj7q8zganyw2lggOT1L34E0YUYqosdgBAbKRa5EpCk8ppgU6tgFsASrnwCQUAe9REIYZB7R8H9u/HtTfMa/F4Tk4O+k3QIq/cgpJaK3pER4hQHYUTBjVRCHHJVbDYXQCAmEguHdoVDreASQuXtXg8c9EsJBmbgtrMCWXU/Tj0TRRCbCojAECvUUKj5Izv7mKK0gIAShnUFAAMaqIQYlV71vbmsHf3SjRoAAB1Vica7Fz4hLoXg5oohNjUnh41g7p7aZQKxOo8/49LzZxQRt2LQU0UQrxD394Qoe5jMnp61TxPTd2NQU0UQmwc+g4YnqemQGFQE4WIBrsTDpVnbW8GdfczGX8Nau6kRd2JQU0UIk6UWwAAESoFItSc8d3d4vVqyGWA1eGGmTtpUTfqVFCfOHHC33UQURcdL/PsQc3edGAo5XIkRHnOU3P4m7pTp4K6X79+uPLKK/Hee+/BauUPKJEU5DbtQc2gDhzveWpOKKPu1Kmg3rt3L0aMGIElS5YgKSkJ99xzD3755Rd/10ZEHeDd3pJBHThnn6cm6i6dCupRo0bhlVdeQVFREf7973+juLgYl1xyCYYNG4YXX3wR5eXl/q6TiM4jt5xBHWimpqHvMrMNAmQiV0OhqkuTyZRKJWbPno21a9fihRdeQG5uLh555BGkpqbitttuQ3Fxsb/qJKJzsDldOF3ZAIBBHUgxkWqoFXI43YLv0jgif+tSUO/Zswf33nsvkpOT8eKLL+KRRx5BXl4eNm7ciKKiIsycOdNfdRLROZyqaIDLLUDusiOSM74DRi6TIbGpV92giRO5GgpVnQrqF198EcOHD8dFF12EoqIivPPOOzh9+jT+9Kc/oXfv3rj00kuxZs0a7N2795yv8/rrr2PEiBEwGAwwGAzIyMjAN99843vearVi8eLFiIuLg16vx5w5c1BaWtrsNfLz8zF9+nTodDokJibi0UcfhdPJSyUovHgnkmkcZshkHIINJO956kYGNXWTTm1z+frrr+OOO+7A7bffjuTk5FaPSUxMxFtvvXXO1+nZsyeef/559O/fH4Ig4O2338bMmTORlZWFoUOH4uGHH8ZXX32FtWvXwmg04r777sPs2bOxY8cOAIDL5cL06dORlJSEn376CcXFxbjtttugUqnw5z//uTNNIwpK3kuzNPZakSsJP97z1I1aBjV1j04F9caNG5GWlga5vHmHXBAEFBQUIC0tDWq1GvPnzz/n68yYMaPZ/eeeew6vv/46du3ahZ49e+Ktt97CBx98gIkTJwIAVq9ejcGDB2PXrl2YMGECvvvuOxw6dAjff/89TCYTRo0ahWeffRaPPfYYnn76aajVPFdH4cHbo9YyqAPO26O2qqNhdbigVfHUA/lXp4a++/bti4qKihaPV1VVoXfv3p0qxOVy4cMPP4TFYkFGRgYyMzPhcDgwefJk3zGDBg1CWloadu7cCQDYuXMnhg8fDpPJ5Dtm6tSpMJvNOHjwYJvvZbPZYDabm92IgtnZQ98UWFEaJSJUCkAmx8Ei/v8n/+tUULe1rm19fT20Wm2HXis7Oxt6vR4ajQb/8z//g3Xr1mHIkCEoKSmBWq1GdHR0s+NNJhNKSkoAACUlJc1C2vu897m2LF++HEaj0XdLTU3tUM1EUuJ0uX3Lh3LoO/BkMplvf+pDRfz/T/7XoaHvJUuWAPD8YD711FPQ6XS+51wuF37++WeMGjWqQwUMHDgQ+/btQ21tLT755BPMnz8fW7du7dBrdNTSpUt9bQEAs9nMsKagVVDdCLvLDa1KDpXTInY5YSkxSoPTlQ3IKWSPmvyvQ0GdlZUFwNOjzs7ObnYOWK1WY+TIkXjkkUc6VIBarUa/fv0AAGPHjsXu3bvxyiuv4KabboLdbkdNTU2zXnVpaSmSkpIAAElJSS1WRPPOCvce0xqNRgONRtOhOomk6nipZyJZ3wQ9ZEdELiZMedf8PljMHjX5X4eCesuWLQCABQsW4JVXXoHB4P8L/N1uN2w2G8aOHQuVSoVNmzZhzpw5AICjR48iPz8fGRkZAICMjAw899xzKCsrQ2JiIgDPRDeDwYAhQ4b4vTYiKfKuSNY/UY9ckWsJV4lNa34fK6mH3emGWsmNCcl/OjXre/Xq1X5586VLl2LatGlIS0tDXV0dPvjgA/zwww/49ttvYTQaceedd2LJkiWIjY2FwWDA/fffj4yMDEyYMAEAMGXKFAwZMgS33norVqxYgZKSEjzxxBNYvHgxe8wUNnKb1vjub4piUIvEoFVC7rLDDjWOl9VhaIpR7JIohLQ7qGfPno01a9bAYDBg9uzZ5zz2008/bddrlpWV+ZYaNRqNGDFiBL799ltcddVVAICXXnoJcrkcc+bMgc1mw9SpU/Haa6/5vl6hUGD9+vVYtGgRMjIyEBkZifnz5+OZZ55pb7OIgt7xphnffRP0IlcSvmQyGSLsVbBEJOFgkZlBTX7V7qA2Go2+FY+MRv/8EJ5vQRStVouVK1di5cqVbR6Tnp6Or7/+2i/1EAUbt1tAnnfo28SgFpPWVu0J6sJa4AJOTiX/aXdQnz3c7a+hbyLqmqLaRjTYXVApZEiP1Z3/C6jbRNiqAYDXUpPfdWrGQ2NjIxoaGnz3T58+jZdffhnfffed3wojovPzDnv3jo+EUsEJTGLS2qoAAIeKzXC7W19rgqgzOvWbPXPmTLzzzjsAgJqaGowbNw5/+9vfMHPmTLz++ut+LZCI2pZX5p3xHSVyJaRx1EGrkqPB7sLJSl7PTv7TqaDeu3cvLr30UgDAJ598gqSkJJw+fRrvvPMOXn31Vb8WSERtO94047tvIs9Pi00GAYOTPZescvib/KlTQd3Q0ICoKM9f8N999x1mz54NuVyOCRMm4PTp034tkIja5t01qz+DWhKGpjQFdSEXPiH/6VRQ9+vXD5999hkKCgrw7bffYsqUKQA8l1t1xyIoRNSSIAi+zTj6MaglYVjTZVnsUZM/dSqon3rqKTzyyCPo1asXxo8f71sp7LvvvsPo0aP9WiARta68zgaz1Qm5zDOZjMQ31BfUtW1uXkTUUZ1amex3v/sdLrnkEhQXF2PkyJG+xydNmoRZs2b5rTgiapu3N50eF8k9kCViQJIeSrkM1Q0OFNVa0SM6QuySKAR0KqgBz6YXv934Yty4cV0uiIjahyuSSY9GqUB/UxQOF5txsLCWQU1+0amgtlgseP7557Fp0yaUlZXB7XY3e/7EiRN+KY6I2uabSMYVySRlaIrBE9RFZkwZ2vYufkTt1amgvuuuu7B161bceuutSE5O9i0tSkSBk1v2665ZJL4D+/fj2hvmodIwAEi4AP/+Ygs2rfojEmKMWL3qtfO/AFEbOhXU33zzDb766itcfPHF/q6HiNqJM76lxeEWMGnhMhTWNOKTzDMQYnth0nVXYtOqP4pdGgW5Ts36jomJQWxsrL9rIaJ2qrbYUVFvB8Bz1FKToPdssVtvc6LB7hS5GgoFnQrqZ599Fk899VSz9b6JKHCOlXrOT/eIjkCkptNzQqkbqJVyREeoAHguoSPqqk79hv/tb39DXl4eTCYTevXqBZVK1ez5vXv3+qU4ImrdkRJPUA9O5hrfUpQQpUFNo4NBTX7RqaC+/vrr/VwGEXXEkRLPylcDkxjUUpQQpcHxsnqU19mgFrsYCnqdCuply5b5uw4i6gBvj3pQEpfslaLEKM956rJ6G3qKXAsFv05vYFtTU4N//etfWLp0KaqqPPuw7t27F4WFhX4rjohacrsFHPUFNXvUUpTQFNQ1DQ64ZJxDQF3TqZ+gAwcOYPLkyTAajTh16hTuvvtuxMbG4tNPP0V+fr5vr2oi8r+C6gY02F1QK+Rc41uidGol9Bol6m1OWDXRYpdDQa5TPeolS5bg9ttvx/Hjx6HVan2PX3PNNdi2bZvfiiOilrzD3v1NeigVnR4Uo27m7VVb1byUlbqmU7/lu3fvxj333NPi8R49eqCkpKTLRRFR244Ue4KaE8mkzXs9daMmRuRKKNh1Kqg1Gg3M5pb7rR47dgwJCQldLoqI2uad8T2YE8kkzdejZlBTF3XqHPV1112HZ555Bh9//DEAQCaTIT8/H4899hjmzJnj1wKJqDnvRLL3/7USn/+9+QhWdk4OJolRFLXgnfltUxthd7qhVvI0BXVOpxc8+d3vfoeEhAQ0Njbi8ssvR0lJCTIyMvDcc8/5u0YiatJgd+JkpQUAMPnGO1usSpa5iPvBS0WUVgmNUg6b07PT2dAUo9glUZDqVFAbjUZs3LgRO3bswP79+1FfX48xY8Zg8uTJ/q6PiM5yuNgMQQCUzkYuHSpxMpkMCVEanKluxMFCM4OaOq3Dv+lutxtr1qzBp59+ilOnTkEmk6F3795ISkqCIAjc8pKoG+UUes5Pa21VIldC7eEL6qJaAKlil0NBqkMnTQRBwHXXXYe77roLhYWFGD58OIYOHYrTp0/j9ttvx6xZHHYj6k7ZhbUAgAgGdVBIbJr5fbCo5eRbovbqUI96zZo12LZtGzZt2oQrr7yy2XObN2/G9ddfj3feeQe33XabX4skIo8cBnVQ8c78PlxshtstQC7niCN1XId61P/5z3/wf//3fy1CGgAmTpyIxx9/HO+//77fiiOiX1kdLhwvqwcAaG3VIldD7RGjU0PmdsJid+FU0yRAoo7qUFAfOHAAV199dZvPT5s2Dfv37+9yUUTU0pGSOrjcAuIi1VC5uBd8MJDLZdDaawBw+Js6r0NBXVVVBZPJ1ObzJpMJ1dX8S5+oO3jPTw/tYQQHUIOHd/SDQU2d1aGgdrlcUCrbPq2tUCjgdDq7XBQRtXSwKaiH9+CKZMEkwu4N6lqRK6Fg1aHJZIIg4Pbbb4dGo2n1eZvN5peiiKglb496WIoRW0WuhdrP26M+VGTmJazUKR0K6vnz55/3GM74JvI/q8PlWzp0eE8unBFMtPYaKOQyVFrsKDXbkGTUnv+LiM7SoaBevXp1d9VBROdw4EwtnG4BiVEa9IiOELsc6gC54EK/BD2OltbhYFEtg5o6jKvEEwWBvfme4dMxaTEcOg1CQ1M88wo4oYw6g0FNFAT2nm4K6vRocQuhThniC2pOKKOOY1ATSZwgCM161BR8vBtyeNdqJ+oIBjWRxBVUNaKi3g6VQoZhPTiRLBh5e9SFNY2oabCLXA0FGwY1kcR5e9NDUozQqhQiV0OdYYxQITXWMwnwEM9TUwcxqIkkzhvUYznsHdSGJntGQzihjDqKQU0kcXtOcSJZKPDO/M7hhDLqIAY1kYTVNNhxuMTTAxvXK1bkaqgrvAvVZJ9hUFPHMKiJJOznk1UQBKBvQiQSDVwoI5iN6BkNADhRYUFto0PcYiioMKiJJGzXiUoAQEbfOJEroa6KjVT7JpSxV00dwaAmkrCdeZ6gntCHQR0KvL3q/WdqRK2DgguDmkiiqix2HGnaiINBHRpGNp2nPsCgpg5gUBNJ1C8nPb3pASY94vWtby1LwcXboz7AoW/qAAY1kURx2Dv0DOthhEwGFNdaUVZnFbscChId2uaSiAJDEARsOVoOALi4X7zI1VBXHNi/H9feMM93X516DWzqaNz5+PP4cuXT4hVGQYNBTSQBCxbei/LqX4dDrSoD8tOuhUxw4RIGdVBzuAVMWrjMd991qASHi+tQ4uC+4tQ+DGoiCSivrm32Yb7ndBVycysR2ViKSA1/TUNJkkGLw8V1aNDylAa1D89RE0nQyXILACDKUihyJeRvyUZPT7pRGw+XWxC5GgoGDGoiiWl0uFBc65loFNXAoA41cZFqqBQyuOUq5JbVi10OBQEGNZHEnK6wQAAQr1dD7WwQuxzyM7lcBlOUZznYrKad0YjOhSe/iCTmRIVn2Lt3fCR+/s2MYa/snBxMCnRh5DdJRi3O1DRib341fj8uTexySOIY1EQS4nILOF3p6UX3jo/E9t/MGPbKXDQr0KWRHyUbvT3qGnELoaDAoCaSkKKaRthdbkSoFEjiblkhK6kpqI+X1WPaTbdD4f51N62EGCNWr3pNrNJIghjURBJy9rC3TCYTuRrqLjq1EkJ9JWT6OAyc9SDS4yJ9z21a9UcRKyMp4mQyIokQBAEnzwpqCm3uylMAgKIaLiVK58agJpKI6gYHahsdUMhkSIvViV0OdTN3xUkAQGFNo8iVkNQxqIkkwtub7hETAbWSv5qhTij3BHWJ2Qqnyy1yNSRl/DQgkogTFZ7FL/pw2DssCPXl0KkVcLkFlJptYpdDEsagJpIAp1yN4qZzlTw/HT56RHuWE+XwN50Lg5pIAup1KRDgWV7SEKESuxwKEAY1tQeDmkgC6nQ9ALA3HW5SmoK6uLYRbm7QQW1gUBOJzOFyo06XDADok8CgDifxejU0SjkcLgFl9TxPTa3jgidEAbRg4b0or65t9li9NhHuHpMRoVLAxNXIwopMJkNKdAROVlhQVN3I1eioVaL2qJcvX44LL7wQUVFRSExMxPXXX4+jR482O8ZqtWLx4sWIi4uDXq/HnDlzUFpa2uyY/Px8TJ8+HTqdDomJiXj00UfhdDoD2RSidimvrsWkhcua3WIm/A4A0CteBzlXIws7PE9N5yNqUG/duhWLFy/Grl27sHHjRjgcDkyZMgUWi8V3zMMPP4wvv/wSa9euxdatW1FUVITZs2f7nne5XJg+fTrsdjt++uknvP3221izZg2eeuopMZpE1CGCIPy6bGgch73D0dlBLQg8T00tiTr0vWHDhmb316xZg8TERGRmZuKyyy5DbW0t3nrrLXzwwQeYOHEiAGD16tUYPHgwdu3ahQkTJuC7777DoUOH8P3338NkMmHUqFF49tln8dhjj+Hpp5+GWq0Wo2lE7VLTtBqZ4HYiLY6rkYWjhCgNVAoZbE43Ki12scshCZLUZLLaWs+5u9jYWABAZmYmHA4HJk+e7Dtm0KBBSEtLw86dOwEAO3fuxPDhw2EymXzHTJ06FWazGQcPHmz1fWw2G8xmc7MbkRi8vWmh/AQ0SoXI1ZAYFHIZko1NvepqDn9TS5IJarfbjYceeggXX3wxhg0bBgAoKSmBWq1GdHR0s2NNJhNKSkp8x5wd0t7nvc+1Zvny5TAajb5bamqqn1tD1D7eZUPdxYdFroTExPPUdC6SCerFixcjJycHH374Ybe/19KlS1FbW+u7FRQUdPt7Ev2W1eFCUa3ng9lVdEjkakhMzc5Ti1wLSY8kgvq+++7D+vXrsWXLFvTs2dP3eFJSEux2O2pqapodX1paiqSkJN8xv50F7r3vPea3NBoNDAZDsxtRoOVXNUAQgFidGmioFrscEpHJoIFCLkOD3QW7ip9H1JyoQS0IAu677z6sW7cOmzdvRu/evZs9P3bsWKhUKmzatMn32NGjR5Gfn4+MjAwAQEZGBrKzs1FWVuY7ZuPGjTAYDBgyZEhgGkLUCacrGwAA6fGcRBbulAo5ko2ea6jrI0znOZrCjaizvhcvXowPPvgAn3/+OaKionznlI1GIyIiImA0GnHnnXdiyZIliI2NhcFgwP3334+MjAxMmDABADBlyhQMGTIEt956K1asWIGSkhI88cQTWLx4MTQajZjNI2qTIAg4XeU5P50eq8Mukesh8aXF6nCmuhH1Ea2PBFL4ErVH/frrr6O2thZXXHEFkpOTfbePPvrId8xLL72Ea6+9FnPmzMFll12GpKQkfPrpp77nFQoF1q9fD4VCgYyMDNxyyy247bbb8Mwzz4jRJKJ2qbTYYbG5oJTLfOcnKbylxnpGViwRJu5PTc2I2qNuz8X9Wq0WK1euxMqVK9s8Jj09HV9//bU/SyPqVt5h754xEVAqJDFVhESWGKWBRimHDWpkF9ZidFqM2CWRRPATgkgEpyqbhr25Ghk1kctk6BnjGV3ZkVshcjUkJQxqogBzuNwoarpeNp2rkdFZvMPf2xnUdBYGNVGAFdU0wi0AUVoloiNUYpdDEpLWFNR7T9eg0e4SuRqSCgY1UYAVNC0T2TMmAjLulkVniY5QQeWwwO5yY/epKrHLIYlgUBMF2Jlqz0Sy1BgOe1NzMpkMkY2ey1R5npq8GNREAeSSq1BmtgGAb+IQ0dn0TUHN89TkxaAmCiCLNhECgGidClFanp+mliIbPUsgHywyo4rbXhIY1EQBZWlaHpK9aWqLymXFQFMUAOCnPPaqiUFNFFDeoOb5aTqXi/vFA+B5avJgUBMFSE2DHVaNZ7UpLhtK53JJ/zgAwI/HGdTEoCYKmD2nPFtZRutUiNSIunovSdy43nFQymU4U92I002r2FH44qcFUYB4r4tlb5rO5cD+/fj9LfOhTpkMZ0Qi5v5hBWLNuUiIMWL1qtfELo9EwB41UYB4gzqFQU3n4HALmLRwGUYMHQQAiBhyJSYtXIby6lqRKyOxMKiJAsDqcCG70PNByx41tYd3OdEz1Y1wt2OnQQpdDGqiAMjKr4HDJUDpbIBByzNOdH6mKC3USjlsTrdvkRwKTwxqogDwDntHNpZzfW9qF7lchtSm6+3zqxpErobExKAmCgBvUOusZSJXQsHEO/zNoA5vDGqibuZ2C9iXXwMA0FnLxS2Ggoo3qItrG+GS8ZRJuGJQE3Wz42X1qLM5oVMroLVz5i61nzFChSitEm4BaIhIFLscEgmDmqibZeV7FjoZ0dMIGTh7l9pPJpP5etX1EUkiV0NiYVATdbO9TUE9Ji1G5EooGPmCWsegDlcMaqJutrfp/PRoBjV1QmpTUNvU0Sg1W0WuhsTAoCbqRrWNDuSW1QMARqdFi1sMBaUIlQKJURoAwHZu0hGWGNRE3WhfQQ0AID1Oh3i9RtxiKGh5h7+3c9vLsMSgJupGe097zk+PTo0WtxAKamcHtcDlRMMOg5qoG2U19ajHpPP8NHVecrQWMrcT5XU2HC2tE7scCjAGNVE3cbsF36VZnPFNXaGUyxHZtKodz1OHHwY1UTfJK69HndUJrUqOgUlRYpdDQU7fUAKA56nDEYOaqJtkNV2WNaJnNFQK/qpR1+gbPUH984kq2JwukauhQOKnB1E34UIn5E8aew3i9Ro0OlzYe7pG7HIogBjURN3ky52HAABff7ga194wD9feMA/ZOTkiV0XBSgbgkn5xAIDtudzcJZxwOxaibmC2OmBR6AEAU26Yj0iN51ctc9EsMcuiIHdJ/wR8tq8I249X4NGpYldDgcIeNVE32F9QA8hkMGiVvpAm6ooD+/fjzb88DcDz83X1TQuwYOG94hZFAcGgJuoG3nOISUatuIVQyHC4BVx95x8QG6kGZDL0mXEfyqu5bWo4YFATdQPvRLJkY4TIlVCoSYvxrFKWX9UgciUUKAxqIj87e6GTZPaoyc9S4zx//DGowweDmsjPTlTUw2x1QuZ2ciMO8rue0TrIZYDZ6oRdGSl2ORQADGoiP8ts2ogjwlYJhVwmcjUUatRKuW/uQ70uWeRqKBAY1ER+5p1IprNyqUfqHt7dtOojkkSuhAKBQU3kZ96JZAxq6i7eoLZEmOByc9vLUMegJvKj2kYHjpfVAwB0NgY1dQ9TlBZqpRwuhQY5hbxEK9QxqIn8yDvbOz1OB6XLJnI1FKrkchlSYzyzv7mbVuhjUBP50d6mHbPGciMO6mapTddTc3/q0MegJvIjb496dDqDmrpXWpwnqDNPV6PRzm0vQxmDmshP3G4B+5p61GPSokWthUJfdIQKKocFdpcbv5yqErsc6kYMaiI/OV5WjzqbEzq1AgNNUWKXQyFOJpMhsrEEALD9OLe9DGUMaiI/8S50Mio1GkoFf7Wo++m9QZ1bKXIl1J34aULkJ97rp8dwIhkFiDeoDxebUV7HqwxCFYOayE98QZ0eLW4hFDaULhuGJBsAAD/lcfZ3qGJQE/lBtcWOE+UWAMDoVPaoKXAu6R8PgJdphTIGNZEfeHvTfRIiEROpFrkaCieX9GsK6twKCAKXEw1FDGoiP/BeHnNheqzIlVC4ubBXLNQKOYprrThRYRG7HOoGSrELIAp2Cxbei1268YA2Htu/+hjXfvQCACA7JweTRK6NQl+EWoELesXgp7xKbD9egb4JerFLIj9jUBN1UWlNPWxx8YAATJk1F4YIFQAgc9EskSujcHFJ/3j8lFeJH49XYP5FvcQuh/yMQ99EXdSgjYdbAPQaJaK0/NuXAs97nnrXiUo4XW6RqyF/Y1ATdZFFmwgA6BETAZlMJnI1FI6GphgRrVOh3ubE/jM1YpdDfsagJuqihogEAECP6AiRK6FwpZDLcFHfOADA9uNcpSzUMKiJusDudKNB4xl2ZFCTmC7p5/mDcXsu1/0ONQxqoi7YV1ADQa5EhEqBGJ1K7HIojHnPU2fl18BsdYhcDfkTg5qoC7zLNqby/DSJLC1Ohz7xkXC6BezgKmUhhUFN1AU/Ne1a1DNWJ3IlRMAVAz0TGzcfKRO5EvInXktC1EkNdieyCjxLh6bG8Pw0Bd6B/ftx7Q3zfPfrI0xAyiR8vjsXL8wZAbmcozyhgEFN1Em/nKyCwyVA5aiHMYLnpynwHG4BkxYu8913ut1Yte0E7NDiULEZw3oYRayO/IVD30Sd9FOeZ9g7srGU56dJEpRyOdKaTsNw+Dt0MKiJOsk7kUzfWCJyJUS/6hUXCQDYcpRBHSoY1ESdUGWx42CRGYCnR00kFelxnh71voIaVFnsIldD/sCgJuqEbcfKIQjA4GQDVC6r2OUQ+URpVdDaqiEIwNZj7FWHAlGDetu2bZgxYwZSUlIgk8nw2WefNXteEAQ89dRTSE5ORkREBCZPnozjx483O6aqqgrz5s2DwWBAdHQ07rzzTtTX1wewFRSOfmgaVrx8QILIlRC1pG8oAgBsOcJVykKBqEFtsVgwcuRIrFy5stXnV6xYgVdffRVvvPEGfv75Z0RGRmLq1KmwWn/twcybNw8HDx7Exo0bsX79emzbtg0LFy4MVBMoDLndArY1LShxxUAGNUlPVEMhAGDrsXK43ILI1VBXiXp51rRp0zBt2rRWnxMEAS+//DKeeOIJzJw5EwDwzjvvwGQy4bPPPsPvf/97HD58GBs2bMDu3btxwQUXAAD+/ve/45prrsFf//pXpKSkBKwtFD4OFNaiymJHlEaJsekxYpdD1ILOWgljhAq1jQ5k5Vfjgl6xYpdEXSDZc9QnT55ESUkJJk+e7HvMaDRi/Pjx2LlzJwBg586diI6O9oU0AEyePBlyuRw///xzm69ts9lgNpub3YjayzvsfXG/eKgUkv0VojCWvX8fUOY5Tbho+Zu49oZ5uPaGeViw8F5xC6NOkeynTEmJ55IXk8nU7HGTyeR7rqSkBImJic2eVyqViI2N9R3TmuXLl8NoNPpuqampfq6eQtkPRz3n/TjsTVLlcAuYMH48AMDdcwwmLVyGSQuXoby6VuTKqDMkG9TdaenSpaitrfXdCgoKxC6JgkRZnRX7z9QA+HVdZSIp6hUfCbkMqLTYUd3Ay7SCmWSDOikpCQBQWtr8GtXS0lLfc0lJSSgra375gdPpRFVVle+Y1mg0GhgMhmY3ovbYdLgMggCM7GlEklErdjlEbdKqFOgZ47mmOq+cV8IEM8kGde/evZGUlIRNmzb5HjObzfj555+RkZEBAMjIyEBNTQ0yMzN9x2zevBlutxvjm4Z9iPzpu4OeUypThrb9hyCRVPRN8KxSlldmEbkS6gpRZ33X19cjNzfXd//kyZPYt28fYmNjkZaWhoceegh/+tOf0L9/f/Tu3RtPPvkkUlJScP311wMABg8ejKuvvhp333033njjDTgcDtx33334/e9/zxnf5Hd1Vgd2NG1rOXWo6TxHE4mvb4IeW46Wo8RsRZ3VIXY51EmiBvWePXtw5ZVX+u4vWbIEADB//nysWbMGf/jDH2CxWLBw4ULU1NTgkksuwYYNG6DV/jrk+P777+O+++7DpEmTIJfLMWfOHLz66qsBbwuFvq3HymF3udEnPhJ9E/Ril0N0XpEaJZKNWhTXWpFXzl51sBI1qK+44goIQtsX48tkMjzzzDN45pln2jwmNjYWH3zwQXeUR9TMM6vXA5oeMB/fjRk3vul7PDsnB5NErIvoXPol6FFca0VuWT141X9w4n7URO3QaHehQuWZ5T3xqilINs70PZe5aJZYZRGdVz+THj/mVqCwphF6RYTY5VAnSHYyGZGUbDlaBrdchSitEkkGzvam4GHQqpDSdIVCrT5d5GqoMxjURO3w5X7PJgcDTFGQyWQiV0PUMQNMUQAY1MGKQU10HnVWBzYf8VyvP7DpA48omPQ36SED0KiNw6kKTioLNgxqovP4/nApbE431PZaxOvVYpdD1GE6tRKpsZ7FTz7fVyRyNdRRDGqi8/gsy/PBZqzP57A3Ba1BSZ7RoP/uPXPOq21IehjUROdQUmvFj8c9m3BE150UuRqizuuXqIfc7UB+VQN2n6oWuxzqAAY10Tn8d+8ZuAVgXO9YaJxcL5mCl0ohh6E+HwDwSSY3IgomDGqiNgiCgE8yzwAAbhjbU+RqiLoupu4EAOCrA8VosDtFrobai0FN1IbM09U4WWGBTq3ANcOTxS6HqMt01nKkx+lgsbvw1YFiscuhdmJQE7XhP794hgenD09GpIaL+FHwkwG48YJUAMD7P+eLWwy1G4OaqBVVFju+POCZ7X3z+DSRqyHynxsvSIVKIcO+ghrkFNaKXQ61A4OaqBUf7ymA3enG8B5GjEqNFrscIr9JiNJgatN+6h/8wl51MGBQE/2Gyy3gvV2nAQC3TkjntdMUcuaN9ywl+llWIfepDgIMaqLf+OFoGc5UN8IYocKMkSlil0PkdxP6xKJvQiQa7C7flQ0kXQxqot9Ytc1zCctNF6YiQq0QuRoi/5PJZLj9ol4AgDU/nYLLzZXKpIxBTXSWfQU1+PlkFZRyGRZc3Evscoi6zZyxPWGMUOF0ZYNv0xmSJl5zQnSWVdvyAADXjUpBsjFC5GqI/OvA/v249oZ5vvuq2JFAzFA8+ta32Pe320SsjM6FQU3U5FSFBRtySgAACy/rI3I1RP7ncAuYtHCZ736d1YHVP51CjSoOB87UYETPaPGKozYxqIma/GNLLtwCEGsvwyP3L2rxfHZODiaJUBdRd4nSqjDIFIXDJXX4x+ZcrLrtArFLolYwqIng6U2vyyoEAESXZTXrdXhlLpoV6LKIut0FvWJxuNiM7w6V4kiJGYOSDGKXRL/ByWREAP6+ORcut4ArByZAZ6sUuxyigImNVMNg8Sx8snJLnsjVUGsY1BT2csvq8dk+T2/6ockDRK6GKPASqw8CANYfKEJeObdzlRoGNYW9FRuOwOUWMHmwCSO5XCiFIa29BpMHmyAIwOs/sFctNQxqCmu7T1Xhu0OlkMuAx6cNFLscItHcN7EfAGBdViEKqhpErobOxqCmsCUIAv789WEAwE0XpqFfYpTIFRGJZ1RqNC7tHw+XW8DrW9mrlhIGNYWtz/YVIiu/BhEqBR6e3F/scohE98Akz+/Bx7sLcLLCInI15MXLsygs1ducWP71EQCeIb9Eg1bkiojEc/aKZfqky1Ef2QMzn34HY+3ZWL3qNZGrIwY1haW/bzqOsjobesXpcNelvcUuh0hUZ69YVllvw/s/58OsT8OpwiMiV0YAg5rC0MGiWvxr+0kAgJDzNebM/Vez57kCGYWzOL0GQ1MMyCkyoyRuNARB4J7sImNQU1hxutz4wycH4HILMNTn4/pb7m5xDFcgo3A3oU8cjpbWoVGbgG9ySnDN8GSxSwprnExGYWXVjydwsMgMY4QKKRV7xC6HSJIiNUqMSYsBALyw4QjsTrfIFYU3BjWFjbzyerz8/XEAwFPXDoHSZRW5IiLpGpMWA6WzEacrG/DurtNilxPWGNQUFtxuAY//9wDsTjcuH5CA2WN6iF0SkaSplXIkVh0AALy88RjKzPzDViwMagp5Cxbei4vvfga7T1VD7nageOObmHHjLcjOyRG7NCJJi6k7gZGp0aizOfHsV4fFLidsMagp5BVYZCgzXQgAuGxwCqbd8SgmLVwGu8MpcmVE0iaDgOeuHwa5DPhyfxG2HSsXu6SwxKCmkNZgd6LAdDFcbgG94nQY0cModklEQWVYDyNuy+gFAFj6aTbqbfwDN9AY1BTSnvnyEGxqIyLVClw1xMTrQYk64dGpA9EzJgKFNY14/hsOgQcar6OmkPXl/iJ8uLsAEARMHZoEnZo/7kQdcfbSosoIE5AyCe/tykf2dx/j89eeFbm68MFPLgpJ+ZUN+L9PswEACTUHkRo7QOSKiILP2UuLAsCWI2U4UFiLg7oRqLLYERupFrG68MGhbwo5DXYnFr67B3U2J8akRSOxKlvskohCwiX94xGrU8Op1OEPn+yHIAhilxQW2KOmkCIIAh795ACOlNQhXq/GP24eg4Vb+GFC5A8qhRxXD0vC+ztP4PvDZZiw8Dkk1Px6zjohxsjdtroBg5pCyhtbT+CrA8VQymV4/ZaxSImOELskopCSEKWBc/8XUI2ZjbK40ciYNB3pcZEAgE2r/ihydaGJQ98UMn44WoYV33q25Xv6uqG4sFesyBURhSb3iV0YmmKAAOCbnBLUNNjFLimkMagpJOSW1eGB/2RBEIC541Ixb3ya2CURhbQrBiYgyaCFzenG+gPF3LijGzGoKeiV19lw++rdMFudGJseg6evG8rrpYm6mVIux/QRyYhUK1BpseO7QyUQwN+77sBz1BTUGuxO3PX2bpypbkSEywLLD//FnM0rmx2TnZODSSLVRxTK9Bolpo9Ixn8zC5FXbkFs/FgIgsA/lP2MQU1BZ8HCe1FeXQsBMuQnXYq6yJ5QuKxo2PIGpi5/o8XxmYtmiVAlUXhINkZgylATvskpQZVxAF7+/jgevorrFvgTh74p6JRX12LSwmVQXb7IE9JyGWaP6wd7TanYpRGFpQGmKFw5MAEA8Mqm43hn5ylxCwoxDGoKSln51dh3pgYAMHWIiZdhEYlsRM9o3/7Vy744iM/3FYpcUehgUFPQMUf2xLbjFQCAS/rFo78pSuSKiAgAEqpzcFtGOgQBePijffhv5hmxSwoJDGoKKjvzKlGQeDEAYHgPI8akRYtbEBH5yAA8PWMofn9hKtwC8Mgn+/HBz/lilxX0GNQUNLLP1OLud/ZAkCvQJz4SVwxI4OxSIomRy2X486zhmN/Us/6/ddlYs+Ok2GUFNQY1BYW88nrMX/0L6m1ORDaWYtqwJMjlDGkiKZLLZXj6uqFYeFkfAMDTXx7Ci98dhdvNdfc7g5dnkeQV1TTitrd+QZXFjmE9DHD9+DGUikvELouIfuPs/asBQACQEDMc5bHD8ermXOSVW/DXG0YiQq0Qr8ggxKAmSSuoasDcN3ehsKYRfeIjsWbBONy+reW10kQkvt/uX+31yX/eRVlyBr7KLkZ+VQPevO0CJBm1IlQYnDj0TZJ10/88gonPrceZ6kaoHXVQZr6P2xfcgeycHLFLI6IOiKk7gffvmoDYSDWyC2tx3T+2Y3vTlRt0fgxqkqTjpXXI1I+HQxWJGJ0Kt145AtPueBSTFi6D3eEUuzwi6qBxvWPx+eKLMcCkR1mdDbe89TP++OVBWB0usUuTPA59k+Rk5Vfjrrf3wKnUIU6vxuzRPaBT80eVKNilxurw2eKLsfzrI3h312ms3nEKH2zZh56lOxFhr/YdlxBjxOpVr4lYqbTw048k5fN9hXj0kwOwO93Q2qow57ILEKHixBOiUKFTK/Hs9cMwcXAi7n5zG2zqaJxInYbhPYzI6BsHrUqBTav+KHaZksKgJkmwOlz401eH8N4uz+IIkwcnovCrjxGhGi9yZUTUVb+dDe7lOHYCA29/HsdK63GgsBZHSuowoqcRTrlGhCqli0FNoss8XYWln2bjWGk9AODeK/rif6cMxMz1PBdNFAramg2euWgWpg1LxrCUBmw7Xo6Kejv2nK6GLH0m/rT+EOZf1AupsToRKpYWBjWJpqimES9/fwwf7/GsBxyvV+PFG0fhsgEJIldGRIGUGqvDzePScKLCgl9OVqGsDvjX9pP41/aTyOgTh9+N7Ylpw5Ow+L4HUF5d2+LrQ/2cNoOaAu54aR3W/HQKazPPwO50AwBuvKAnlk4bjJhItcjVEZEYZDIZ+ibo0Sc+Ep+/9yZSL7sBO/IqsPNEJXaeqMRTn+dAqRqGi2aOR6+4SKiVv160FOrntBnU1CULFt7b4i/c3/5163ILOFZah515lfjyQBGy8mt8z43rHYvHrh6IsemxgSqZiCRMJpMhqrEY7901HoU1jfhv5hl8knkG+VUNQFQ6vskpgUIuQ3qsDv0S9egdHyl2yd2OQU1dUl5d2+Lc07f/eh7fHypFVkE1svJrsL+gBhb7WddKCm4YLIWIrT0KS14Z5v7jEAYMGtLitbNzcjCpuxtARJLVIzoCD0zqj/uu7If9Z2pw93P/giNlJGobHThRYcGJCgvkMkCXfCXe3XUaVw5MQM+Y0DunzaCmLhEgQ6nZipJaK4qb/lvbew7uemdPs+P0GiWE6jO4YNRw9E/UI1IzEMBEAMDeRbPanGhCRCSXyzA6LQZJVfswcc51qKi3I7e8Hnll9ai02FGvS8aTn3lWLBxg0uPKgYm4YmAiLugVA5Ui+Nf1YlBTh5TUWpGVX42sghpk5VfjUO8bcHB3QYvj+ifqMTotGqPTYjAmLQb9EvWYedMtGDXjUhGqJqJQIZPJkBClQUKUBhl94lDdYMeaN9+Avu9YNGjjcay0HsdK6/HPbSegEBy4ZGAKLkiPwQW9YjEqNTooNwQJmaBeuXIl/vKXv6CkpAQjR47E3//+d4wbN07ssoJabaMDx0rrsL+gBln5NdibX43iWmvzg+RKaJRyJBm1SDZokWTU4ui6V/DN82tEqZmIgl9b1123djosRqeG48gW3P3gA7A6XMivasDJCgtOVzag0QFsPVaOrcfKAQByGdA7PhKDkg0YkmxA3wQ90uN0SI3VQa+RbhxKt7IO+Oijj7BkyRK88cYbGD9+PF5++WVMnToVR48eRWJiotjlSZLd6YbZ6oC50YEqix1FtVYU1zSiuNaKkxUWHCutaxnKACC4obXXQmetQIStAmeyd+Ke516DTPbr3tC5bkcAW0JEoeZc112fi1alwABTFAaYoiAIAl5+4n70GH4xLBEJaNAmwKnUIa/cgrxyC746UNzsa+Mi1UiL0yEtVtes42EyaNEnQQ9jhMqvbeyIkAjqF198EXfffTcWLFgAAHjjjTfw1Vdf4d///jcef/zxgNXx9k+ncLLC0qGvOSvf2s3lFuB0C3C63HC6mv7tdsPhanrMLTQ9/uu/HU2P11udqG10oLGdC+EnGbQY1sOI7O3f4pKrpiMxStvssogV29Y1C2mgY38NExF1B5lMBkdVIW78/VwAgCAIaLC7sPHjf2PunYtwpNiMkxUW5Fc1oLrBgUqLHZUWe7OrUrxemDMcN12YFuAW/Crog9putyMzMxNLly71PSaXyzF58mTs3Lmz1a+x2Wyw2Wy++7W1nsuLzGZzl2pZvycPP5+s6tJrBJpKKYdWKUOkRokojQpVx3bj3ltno3+CHn0To3x/Rd6wbgXi1dPgtjXA+uv/OrjdLlgt9c1e0+Z04uJ5/9vivXYvubnFsW29Rqg+LqVapPS4zdrg+7e1oR6C2y25GqX4uJRqCYbHFQDydnyFdeZ832OJAOLlShzJL8PUe56E2eqExeqCxe5Evc2JquoqGJXOLueDV1RUVIvOzXkJQa6wsFAAIPz000/NHn/00UeFcePGtfo1y5YtEwDwxhtvvPHGW0BvtbW1Hc65oO9Rd8bSpUuxZMkS3323242qqirExcV16C8ds9mM1NRUFBQUwGAwdEepognltgGh3T62LXiFcvvYNo+oqKgOv37QB3V8fDwUCgVKS0ubPV5aWoqkpKRWv0aj0UCjab47S3R0dKdrMBgMIfeD5xXKbQNCu31sW/AK5faxbR0X9FeCq9VqjB07Fps2bfI95na7sWnTJmRkZIhYGRERUdcFfY8aAJYsWYL58+fjggsuwLhx4/Dyyy/DYrH4ZoETEREFq5AI6ptuugnl5eV46qmnUFJSglGjRmHDhg0wmUzd+r4ajQbLli1rMYweCkK5bUBot49tC16h3D62rfNkgiAI3fLKRERE1GVBf46aiIgolDGoiYiIJIxBTUREJGEMaiIiIgljULdi27ZtmDFjBlJSUiCTyfDZZ581e760tBS33347UlJSoNPpcPXVV+P48eO+56uqqnD//fdj4MCBiIiIQFpaGh544AHfmuJi6mrbAOCee+5B3759ERERgYSEBMycORNHjhwJYCta54+2eQmCgGnTprX6OmLwR9uuuOIKyGSyZrf/+Z//CWAr2uav793OnTsxceJEREZGwmAw4LLLLkNjY2OAWtG6rrbt1KlTLb5v3tvatWsD3Jrm/PF9Kykpwa233oqkpCRERkZizJgx+O9//xvAVrTOH23Ly8vDrFmzkJCQAIPBgBtvvLHF4lztwaBuhcViwciRI7Fy5coWzwmCgOuvvx4nTpzA559/jqysLKSnp2Py5MmwWDw7ZxUVFaGoqAh//etfkZOTgzVr1mDDhg248847A92UFrraNgAYO3YsVq9ejcOHD+Pbb7+FIAiYMmUKXK727cjVXfzRNq+XX3654wvndyN/te3uu+9GcXGx77ZixYpANeGc/NG+nTt34uqrr8aUKVPwyy+/YPfu3bjvvvsgl4v7MdfVtqWmpjb7nhUXF+OPf/wj9Ho9pk2bFujmNOOP79ttt92Go0eP4osvvkB2djZmz56NG2+8EVlZWYFsSgtdbZvFYsGUKVMgk8mwefNm7NixA3a7HTNmzID7rE1n2qXDq4OHGQDCunXrfPePHj0qABBycnJ8j7lcLiEhIUF4880323ydjz/+WFCr1YLD4ejOcjvEX23bv3+/AEDIzc3tznI7pCtty8rKEnr06CEUFxe3eB0p6GzbLr/8cuHBBx8MYKWd09n2jR8/XnjiiScCWWqH+et3btSoUcIdd9zRnaV2WGfbFhkZKbzzzjvNXis2Nvac7Q+0zrTt22+/FeRyebNNOGpqagSZTCZs3LixQ+/PHnUHebfH1Gq1vsfkcjk0Gg22b9/e5tfV1tbCYDBAqZTuGjOdaZvFYsHq1avRu3dvpKamBqTOzmhv2xoaGnDzzTdj5cqVba4VLzUd+b69//77iI+Px7Bhw7B06VI0NDRA6trTvrKyMvz8889ITEzERRddBJPJhMsvv/ycv5NS0JnfuczMTOzbt08SI3Tn0t62XXTRRfjoo49QVVUFt9uNDz/8EFarFVdccUWgS2639rTNZrNBJpM1WwRFq9VCLpd3+OeSQd1BgwYNQlpaGpYuXYrq6mrY7Xa88MILOHPmDIqLi1v9moqKCjz77LNYuHBhgKvtmI607bXXXoNer4der8c333yDjRs3Qq1Wi1T5+bW3bQ8//DAuuugizJw5U8RqO6a9bbv55pvx3nvvYcuWLVi6dCneffdd3HLLLSJW3j7tad+JEycAAE8//TTuvvtubNiwAWPGjMGkSZPanIcgBZ35PHnrrbcwePBgXHTRRQGutmPa27aPP/4YDocDcXFx0Gg0uOeee7Bu3Tr069dPxOrPrT1tmzBhAiIjI/HYY4+hoaEBFosFjzzyCFwuV5vf27YwqDtIpVLh008/xbFjxxAbGwudToctW7Zg2rRprZ4LM5vNmD59OoYMGYKnn3468AV3QEfaNm/ePGRlZWHr1q0YMGAAbrzxRlitVpEqP7/2tO2LL77A5s2b8fLLL4tbbAe19/u2cOFCTJ06FcOHD8e8efPwzjvvYN26dcjLyxOx+vNrT/u85/zuueceLFiwAKNHj8ZLL72EgQMH4t///reY5Z9TRz9PGhsb8cEHH0i+Nw20v21PPvkkampq8P3332PPnj1YsmQJbrzxRmRnZ4tY/bm1p20JCQlYu3YtvvzyS+j1ehiNRtTU1GDMmDEdnjch3XFYCRs7diz27duH2tpa2O12JCQkYPz48bjggguaHVdXV4err74aUVFRWLduHVQqlUgVt19722Y0GmE0GtG/f39MmDABMTExWLduHebOnStS5ed3vrZt3rwZeXl5LbY8nTNnDi699FL88MMPgS+6ndr7fTvb+PHjAQC5ubno27dvoErtlPO1Lzk5GQAwZMiQZl83ePBg5OfnB7zejujI9+6TTz5BQ0MDbrvtNhEq7bjztS0vLw//+Mc/kJOTg6FDhwIARo4ciR9//BErV67EG2+8IWb559Se79uUKVOQl5eHiooKKJVKREdHIykpCX369OnQe7FH3QVGoxEJCQk4fvw49uzZ02y41Gw2Y8qUKVCr1fjiiy+ancsIBudq228JggBBEHznbaSurbY9/vjjOHDgAPbt2+e7AcBLL72E1atXi1hx+3Xk++ZtnzfkgkFb7evVqxdSUlJw9OjRZscfO3YM6enpYpTaYe353r311lu47rrrkJCQIEKFnddW27xzJH7bw1QoFB2fGS2S9nzf4uPjER0djc2bN6OsrAzXXXddh96DPepW1NfXIzc313f/5MmT2LdvH2JjY5GWloa1a9ciISEBaWlpyM7OxoMPPojrr78eU6ZMAfBrSDc0NOC9996D2WyG2WwG4BkOUSgUorQL6HrbTpw4gY8++ghTpkxBQkICzpw5g+effx4RERG45pprxGoWgK63LSkpqdUJZGlpaejdu3fA2tGarrYtLy8PH3zwAa655hrExcXhwIEDePjhh3HZZZdhxIgRYjXLp6vtk8lkePTRR7Fs2TKMHDkSo0aNwttvv40jR47gk08+EatZALreNq/c3Fxs27YNX3/9daCb0Kautm3QoEHo168f7rnnHvz1r39FXFwcPvvsM2zcuBHr168Xq1kA/PN9W716NQYPHoyEhATs3LkTDz74IB5++GEMHDiwY8V0Zcp6qNqyZYsAoMVt/vz5giAIwiuvvCL07NlTUKlUQlpamvDEE08INpvtvF8PQDh58qQ4jTpPbe1tW2FhoTBt2jQhMTFRUKlUQs+ePYWbb75ZOHLkiEgt+lVX29YaSOTyrK62LT8/X7jsssuE2NhYQaPRCP369RMeffTRZpeOiMlf37vly5cLPXv2FHQ6nZCRkSH8+OOPAW5JS/5q29KlS4XU1FTB5XIFuAVt80fbjh07JsyePVtITEwUdDqdMGLEiBaXa4nBH2177LHHBJPJJKhUKqF///7C3/72N8Htdne4Fm5zSUREJGE8R01ERCRhDGoiIiIJY1ATERFJGIOaiIhIwhjUREREEsagJiIikjAGNRERkYQxqImIiCSMQU1E5ySTyfDZZ5+JXQZR2GJQE4msvLwcixYtQlpaGjQaDZKSkjB16lTs2LFD7NL84tSpU5DJZFAoFCgsLGz2XHFxMZRKJWQyGU6dOiVOgUQSx6AmEtmcOXOQlZWFt99+G8eOHcMXX3yBK664ApWVlWKX5lc9evTAO++80+yxt99+Gz169BCpIqLgwKAmElFNTQ1+/PFHvPDCC7jyyiuRnp6OcePGYenSpc22wnvxxRcxfPhwREZGIjU1Fffeey/q6+t9z69ZswbR0dFYv349Bg4cCJ1Oh9/97ndoaGjA22+/jV69eiEmJgYPPPAAXC6X7+t69eqFZ599FnPnzkVkZCR69OiBlStXnrPmgoIC3HjjjYiOjkZsbCxmzpzZrt7w/PnzW2wXunr1asyfP7/FsTk5OZg2bRr0ej1MJhNuvfVWVFRU+J7fsGEDLrnkEkRHRyMuLg7XXnst8vLyfM97e/GffvoprrzySuh0OowcORI7d+48b51EUsOgJhKRXq+HXq/HZ599ds79vOVyOV599VUcPHgQb7/9NjZv3ow//OEPzY5paGjAq6++ig8//BAbNmzADz/8gFmzZuHrr7/G119/jXfffRf//Oc/W2z7+Je//AUjR45EVlYWHn/8cTz44IPYuHFjq3U4HA5MnToVUVFR+PHHH7Fjxw7o9XpcffXVsNvt52zrddddh+rqamzfvh0AsH37dlRXV2PGjBnNjqupqcHEiRMxevRo7NmzBxs2bEBpaSluvPFG3zEWiwVLlizBnj17sGnTJsjlcsyaNavFHsb/7//9PzzyyCPYt28fBgwYgLlz58LpdJ6zTiLJ6fpmYETUFZ988okQExMjaLVa4aKLLhKWLl0q7N+//5xfs3btWiEuLs53f/Xq1QIAITc31/fYPffcI+h0OqGurs732NSpU4V77rnHdz89PV24+uqrm732TTfdJEybNs13H2dt9fnuu+8KAwcObLZVn81mEyIiIoRvv/221VpPnjwpABCysrKEhx56SFiwYIEgCIKwYMEC4eGHHxaysrKabQH77LPPClOmTGn2GgUFBQIA4ejRo62+R3l5uQBAyM7Obvae//rXv3zHHDx4UAAgHD58uNXXIJIq9qiJRDZnzhwUFRXhiy++wNVXX40ffvgBY8aMwZo1a3zHfP/995g0aRJ69OiBqKgo3HrrraisrERDQ4PvGJ1Oh759+/rum0wm9OrVC3q9vtljZWVlzd4/IyOjxf3Dhw+3Wuv+/fuRm5uLqKgo32hAbGwsrFZrs6Hnttxxxx1Yu3YtSkpKsHbtWtxxxx2tvseWLVt8r6/X6zFo0CAA8L3H8ePHMXfuXPTp0wcGgwG9evUCAOTn5zd7rREjRvj+nZycDAAt2k8kdUqxCyAiQKvV4qqrrsJVV12FJ598EnfddReWLVuG22+/HadOncK1116LRYsW4bnnnkNsbCy2b9+OO++8E3a7HTqdDgCgUqmavaZMJmv1sd8OD3dEfX09xo4di/fff7/FcwkJCef9+uHDh2PQoEGYO3cuBg8ejGHDhmHfvn0t3mPGjBl44YUXWny9N2xnzJiB9PR0vPnmm0hJSYHb7cawYcNaDL+f3X6ZTAYAXWo/kRgY1EQSNGTIEN+1y5mZmXC73fjb3/4GudwzCPbxxx/77b127drV4v7gwYNbPXbMmDH46KOPkJiYCIPB0Kn3u+OOO3Dvvffi9ddfb/M9/vvf/6JXr15QKlt+RFVWVuLo0aN48803cemllwKA77w3USji0DeRiCorKzFx4kS89957OHDgAE6ePIm1a9dixYoVmDlzJgCgX79+cDgc+Pvf/44TJ07g3XffxRtvvOG3Gnbs2IEVK1bg2LFjWLlyJdauXYsHH3yw1WPnzZuH+Ph4zJw5Ez/++CNOnjyJH374AQ888ADOnDnTrve7++67UV5ejrvuuqvV5xcvXoyqqirMnTsXu3fvRl5eHr799lssWLAALpcLMTExiIuLw6pVq5Cbm4vNmzdjyZIlnW4/kdQxqIlEpNfrMX78eLz00ku47LLLMGzYMDz55JO4++678Y9//AMAMHLkSLz44ot44YUXMGzYMLz//vtYvny532r43//9X+zZswejR4/Gn/70J7z44ouYOnVqq8fqdDps27YNaWlpmD17NgYPHow777wTVqu13T1spVKJ+Pj4VnvLAJCSkoIdO3bA5XJhypQpGD58OB566CFER0dDLpdDLpfjww8/RGZmJoYNG4aHH34Yf/nLXzrdfiKpkwmCIIhdBBGJo1evXnjooYfw0EMPiV0KEbWBPWoiIiIJY1ATERFJGIe+iYiIJIw9aiIiIgljUBMREUkYg5qIiEjCGNREREQSxqAmIiKSMAY1ERGRhDGoiYiIJIxBTUREJGH/H91XLoAsXFSGAAAAAElFTkSuQmCC",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"sns.displot(boot_means, kde=True)\n",
"plt.axvline(np.mean(chinstraps['FlipperLength']), color='black')\n",
"plt.xlabel('Sample Mean')\n",
"plt.ylabel('Density')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can then compute the confidence interval with the quantiles of this distribution:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([194.14705882, 197.5 ])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.quantile(boot_means, [0.025, 0.975])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Penguin Statistics\n",
"\n",
"With that warmup, let's look at our various biometrics."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Parametric Estimates\n",
"\n",
"We're going to compute CIs for several measurements, but we don't want to repeat all of our code.\n",
"\n",
"Pandas `groupby` let us *apply* a function to each group, which can in turn return a series or a data frame.\n",
"\n",
"Let's write a function that, given a series of penguin data, returns our statistics:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"def mean_estimate(vals):\n",
" # vals is a series of measurements of a single variable\n",
" mean = vals.mean()\n",
" se = vals.sem() # Pandas has an SEM function too.\n",
" \n",
" ci_width = 1.96 * se\n",
" return pd.Series({\n",
" 'mean': mean,\n",
" 'std': vals.std(),\n",
" 'count': vals.count(),\n",
" 'se': vals.sem(),\n",
" 'ci_width': ci_width,\n",
" 'ci_min': mean - ci_width,\n",
" 'ci_max': mean + ci_width\n",
" })"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Flipper Length\n",
"\n",
"Now we can do this to the flipper length:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
mean
\n",
"
std
\n",
"
count
\n",
"
se
\n",
"
ci_width
\n",
"
ci_min
\n",
"
ci_max
\n",
"
\n",
"
\n",
"
species
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
Adelie
\n",
"
189.953642
\n",
"
6.539457
\n",
"
151.0
\n",
"
0.532173
\n",
"
1.043060
\n",
"
188.910582
\n",
"
190.996702
\n",
"
\n",
"
\n",
"
Chinstrap
\n",
"
195.823529
\n",
"
7.131894
\n",
"
68.0
\n",
"
0.864869
\n",
"
1.695144
\n",
"
194.128386
\n",
"
197.518673
\n",
"
\n",
"
\n",
"
Gentoo
\n",
"
217.186992
\n",
"
6.484976
\n",
"
123.0
\n",
"
0.584731
\n",
"
1.146072
\n",
"
216.040920
\n",
"
218.333064
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" mean std count se ci_width ci_min \\\n",
"species \n",
"Adelie 189.953642 6.539457 151.0 0.532173 1.043060 188.910582 \n",
"Chinstrap 195.823529 7.131894 68.0 0.864869 1.695144 194.128386 \n",
"Gentoo 217.186992 6.484976 123.0 0.584731 1.146072 216.040920 \n",
"\n",
" ci_max \n",
"species \n",
"Adelie 190.996702 \n",
"Chinstrap 197.518673 \n",
"Gentoo 218.333064 "
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"penguins.groupby('species')['FlipperLength'].apply(mean_estimate).unstack()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The confidence intervals don't overlap. Chinstrap and Adelie are the closest.\n",
"\n",
":::{note}\n",
"The `unstack` function pivots the innermost level of a hierarchical index to be column labels. Try without it to see what changes!\n",
":::"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Bill Length"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
" mean std count se ci_width ci_min \\\n",
"species \n",
"Adelie 3700.662252 458.566126 151.0 37.317582 73.142461 3627.519791 \n",
"Chinstrap 3733.088235 384.335081 68.0 46.607475 91.350650 3641.737585 \n",
"Gentoo 5076.016260 504.116237 123.0 45.454630 89.091075 4986.925185 \n",
"\n",
" ci_max \n",
"species \n",
"Adelie 3773.804713 \n",
"Chinstrap 3824.438885 \n",
"Gentoo 5165.107336 "
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"penguins.groupby('species')['Mass'].apply(mean_estimate).unstack()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Adelie and Chinstrap are very similar, with Gentoo substantially larger. Look at the CI intervals again!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Bootstrap\n",
"\n",
"Now we're going to bootstrap confidence intervals. I'm just going to do it for the flipper length — you should try it for the others!\n",
"\n",
"Let's take that bootstrap logic from beofre and put it in a function:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"def boot_mean_estimate(vals, nboot=10000):\n",
" obs = vals.dropna() # ignore missing values\n",
" mean = obs.mean()\n",
" n = obs.count()\n",
" \n",
" boot_means = [np.mean(rng.choice(obs, size=n)) for i in range(nboot)]\n",
" ci_low, ci_high = np.quantile(boot_means, [0.025, 0.975])\n",
" return pd.Series({\n",
" 'mean': mean,\n",
" 'count': n,\n",
" 'ci_low': ci_low,\n",
" 'ci_high': ci_high\n",
" })"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
mean
\n",
"
count
\n",
"
ci_low
\n",
"
ci_high
\n",
"
\n",
"
\n",
"
species
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
Adelie
\n",
"
189.953642
\n",
"
151.0
\n",
"
188.907285
\n",
"
190.980132
\n",
"
\n",
"
\n",
"
Chinstrap
\n",
"
195.823529
\n",
"
68.0
\n",
"
194.161765
\n",
"
197.500000
\n",
"
\n",
"
\n",
"
Gentoo
\n",
"
217.186992
\n",
"
123.0
\n",
"
216.032520
\n",
"
218.357724
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" mean count ci_low ci_high\n",
"species \n",
"Adelie 189.953642 151.0 188.907285 190.980132\n",
"Chinstrap 195.823529 68.0 194.161765 197.500000\n",
"Gentoo 217.186992 123.0 216.032520 218.357724"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"penguins.groupby('species')['FlipperLength'].apply(boot_mean_estimate).unstack()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
":::{tip}\n",
"The `apply` function also allows us to pass additional named parameters, which are passed directly to the function we are invoking.\n",
":::"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
mean
\n",
"
count
\n",
"
ci_low
\n",
"
ci_high
\n",
"
\n",
"
\n",
"
species
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
Adelie
\n",
"
189.953642
\n",
"
151.0
\n",
"
188.874172
\n",
"
191.033113
\n",
"
\n",
"
\n",
"
Chinstrap
\n",
"
195.823529
\n",
"
68.0
\n",
"
194.147059
\n",
"
197.500000
\n",
"
\n",
"
\n",
"
Gentoo
\n",
"
217.186992
\n",
"
123.0
\n",
"
216.048780
\n",
"
218.341463
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" mean count ci_low ci_high\n",
"species \n",
"Adelie 189.953642 151.0 188.874172 191.033113\n",
"Chinstrap 195.823529 68.0 194.147059 197.500000\n",
"Gentoo 217.186992 123.0 216.048780 218.341463"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"penguins.groupby('species')['FlipperLength'].apply(boot_mean_estimate, nboot=5000).unstack()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the interval Seaborn uses in a catplot of the statistic:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sns.catplot(x='species', y='FlipperLength', data=penguins, kind='bar')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Testing\n",
"\n",
"Let's do a *t*-test for whether the chinstrap and adelie have different flipper lengths. The `ttest_ind` function does a two-sample T-test:"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Ttest_indResult(statistic=-5.780384584564813, pvalue=6.049266635901914e-08)"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sps.ttest_ind(adelies['FlipperLength'], chinstraps['FlipperLength'], nan_policy='omit', equal_var=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$p < 0.05$, that's for sure! The probability of observing mean flipper lengths this different if there is no actual difference is quite low.\n",
"\n",
"Let's try the bill depth:"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Ttest_indResult(statistic=-0.4263555696052568, pvalue=0.6702714045724318)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sps.ttest_ind(adelies['BillDepth'], chinstraps['BillDepth'], nan_policy='omit')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$p>0.05$. Very. The evidence is not sufficient to reject the null hypothesis that Adelie and Chinstrap penguins have the same bill depth."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Bootstrapped Test\n",
"\n",
"What about the bootstrap? Remember that to do the bootstrap, we need to resample from the *pool* of measurements, because the procedure is to take bootstrap samples **from the distribution under the null hypothesis**, which is that there is no difference between penguins.\n",
"\n",
"We'll work through it a piece at a time, then define a function.\n",
"\n",
"First, we need to drop the null values from our observations:"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"obs1 = adelies['FlipperLength'].dropna()\n",
"obs2 = chinstraps['FlipperLength'].dropna()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And the lengths:"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"n1 = len(obs1)\n",
"n2 = len(obs2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And the observed difference in means (since the mean is the statistic we are bootstrapping):"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-5.869887027658734"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diff_mean = np.mean(obs1) - np.mean(obs2)\n",
"diff_mean"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we're going to pool together the observations, so we can simulate the null hypothesis (data drawn from the same distribution):"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 181.0\n",
"1 186.0\n",
"2 195.0\n",
"4 193.0\n",
"5 190.0\n",
" ... \n",
"339 207.0\n",
"340 202.0\n",
"341 193.0\n",
"342 210.0\n",
"343 198.0\n",
"Name: FlipperLength, Length: 219, dtype: float64"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pool = pd.concat([obs1, obs2])\n",
"pool"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the pool, we can now sample *bootstrap samples* of the means of samples of size $n_1$ and $n_2$ from a merged pool:"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"b1 = np.array([np.mean(rng.choice(pool, size=n1)) for i in range(10000)])\n",
"b2 = np.array([np.mean(rng.choice(pool, size=n2)) for i in range(10000)])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We use the list comprehension from before, but this time we wrap its result in `np.array` to make NumPy arrays that support vectorized computation (lists do not). That then allows us to compute the difference in means for all bootstrap samples in one go:"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-0.88206077, -0.28632645, -0.25350604, ..., 0.93523568,\n",
" -0.85196728, -0.31077133])"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"boot_diffs = b1 - b2\n",
"boot_diffs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, our goal is to compute $p$, the probability of observing a difference at least as large as `diff_mean` (in absolute value) under the null hypothesis $H_0$ (the samples are from the same distribution). So we want to compute the absolute values of the differences, both bootstrapped and observed. NumPy helpfully vectorizes this:"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5.869887027658734"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"abs_diff_mean = np.abs(diff_mean)\n",
"abs_diff_mean"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.88206077, 0.28632645, 0.25350604, ..., 0.93523568, 0.85196728,\n",
" 0.31077133])"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"abs_boot_diffs = np.abs(boot_diffs)\n",
"abs_boot_diffs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And now we can actually compare. The `>=` operator compares the bootstrapped differences with observed, giving us a logical vector that is `True` when a bootstrapped difference is at least as large as the observed difference:"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([False, False, False, ..., False, False, False])"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"boot_exceeds_observed = abs_boot_diffs >= abs_diff_mean\n",
"boot_exceeds_observed"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we want to covert this to an estimated $p$-value. The `np.mean` function, when applied to a boolean array, returns the fraction of elements that are `True`:"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.mean(boot_exceeds_observed)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"No elements were `True` (if we ran with a much larger bootstrap sample count, we might see a small number of `True` values). Since our bootstrap sample count estimates the $p$-value to approximately 5 decimal places, we conclude that the difference is statistically significant."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Bootstrap Test Function\n",
"\n",
"Let's go! We're going to put all that logic into a short function for easy reuse. Note that variables assigned within a function are *local* variables: they only have their values within the function, without affecting the global variables of the same name."
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"def boot_ind(s1, s2, nboot=10000):\n",
" ## we will ignore NAs here\n",
" obs1 = s1.dropna()\n",
" obs2 = s2.dropna()\n",
" n1 = len(obs1)\n",
" n2 = len(obs2)\n",
" \n",
" ## pool the observations together\n",
" pool = pd.concat([obs1, obs2])\n",
" ## grab the observed mean\n",
" md = np.mean(s1) - np.mean(s2)\n",
" \n",
" ## compute our bootstrap samples of the mean under H0\n",
" b1 = np.array([np.mean(rng.choice(pool, size=n1)) for i in range(nboot)])\n",
" b2 = np.array([np.mean(rng.choice(pool, size=n2)) for i in range(nboot)])\n",
" \n",
" ## the P-value is the probability that we observe a difference as large\n",
" ## as we did in the raw data, if the null hypothesis were true\n",
" return md, np.mean(np.abs(b1 - b2) >= np.abs(md))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And with this function, we can compute the bootstrap versions of the tests above:"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-5.869887027658734, 0.0)"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"boot_ind(adelies['FlipperLength'], chinstraps['FlipperLength'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The $p$-value was very small in the $t$-test; it makes sense that we would almost never actually see a mean this large under the null!\n",
"\n",
"Let's do the bill depth:"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-0.07423061940007614, 0.6724)"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"boot_ind(adelies['BillDepth'], chinstraps['BillDepth'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Yeah, those kinds of penguins seem to have the same bill depth."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conclusion\n",
"\n",
"We see extremely similar $p$-values and confidence intervals for the parameteric methods and the bootstrap. That is expected — if two procedures that claimed to do the same thing returned wildly different results, we would have a problem.\n",
"\n",
"But as we discussed in the videos, the bootstrap works in places where the parametrics won't."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.10.7 ('venv': venv)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.7"
},
"vscode": {
"interpreter": {
"hash": "83adc33feab0f9fa0b6c69e9e6552376583bf2c1c67bc73919d7687a215c1dbc"
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}