{
"cells": [
{
"cell_type": "markdown",
"id": "572fb94e",
"metadata": {},
"source": [
"# Permutation test for paired data\n",
"\n",
"We first look at the case of paired data - data in which we wish to compare two groups and each datapoint in one group has a counterpart in the other\n",
"\n",
"Experimental designs using paired data include matched pairs (eg brothers and sisters) and repeated measures (measurements of the same individual before- and after- an intervention, or on- and off-drug).\n"
]
},
{
"cell_type": "markdown",
"id": "2c4eabfd",
"metadata": {},
"source": [
"### Set up Python libraries\n",
"\n",
"As usual, run the code cell below to import the relevant Python libraries"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "3286b373",
"metadata": {},
"outputs": [],
"source": [
"# Set-up Python libraries - you need to run this but you don't need to change it\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy.stats as stats\n",
"import pandas \n",
"import seaborn as sns"
]
},
{
"cell_type": "markdown",
"id": "4bd9c055",
"metadata": {},
"source": [
"## Colab users\n",
"\n",
"You need to use a more recent version of scipy.stats than the default. To do this run the following code block and *after* it has run, go to the menus at the top of colab and click `runtime-->Restart Runtime`"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "c5c86d36",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Requirement already satisfied: scipy==1.10.0 in /Users/joreilly/opt/anaconda3/lib/python3.9/site-packages (1.10.0)\n",
"Requirement already satisfied: numpy<1.27.0,>=1.19.5 in /Users/joreilly/opt/anaconda3/lib/python3.9/site-packages (from scipy==1.10.0) (1.21.5)\n"
]
}
],
"source": [
"# Set-up Python libraries - you need to run this but you don't need to change it\n",
"!pip install scipy==1.10.0\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy.stats as stats\n",
"import pandas \n",
"import seaborn as sns"
]
},
{
"cell_type": "markdown",
"id": "df221c50",
"metadata": {},
"source": [
"## Toy example\n",
"\n",
"[A toy example is an example with a very small dataset, just to show how it works]\n",
"\n",
"We are interested in whether men or women own more pairs of socks. We decide on a matched pairs design in which husbands are compared to their wives, as it is hypothesised that lifestyle factors such as the size of the home in which people live and the duration of holidays taken will affect the number of pairs that can be reasonably justified, and these lifestyle factors are generally shared by both members of a married couple.\n",
"\n",
"We obtain sock-counts for the following informal sample of 10 couples:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "e0483d68",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Husband
\n",
"
Wife
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
10
\n",
"
12
\n",
"
\n",
"
\n",
"
1
\n",
"
17
\n",
"
13
\n",
"
\n",
"
\n",
"
2
\n",
"
48
\n",
"
20
\n",
"
\n",
"
\n",
"
3
\n",
"
28
\n",
"
25
\n",
"
\n",
"
\n",
"
4
\n",
"
23
\n",
"
18
\n",
"
\n",
"
\n",
"
5
\n",
"
16
\n",
"
14
\n",
"
\n",
"
\n",
"
6
\n",
"
18
\n",
"
13
\n",
"
\n",
"
\n",
"
7
\n",
"
34
\n",
"
26
\n",
"
\n",
"
\n",
"
8
\n",
"
27
\n",
"
22
\n",
"
\n",
"
\n",
"
9
\n",
"
22
\n",
"
14
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Husband Wife\n",
"0 10 12\n",
"1 17 13\n",
"2 48 20\n",
"3 28 25\n",
"4 23 18\n",
"5 16 14\n",
"6 18 13\n",
"7 34 26\n",
"8 27 22\n",
"9 22 14"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"socks = pandas.DataFrame(data=[[10,12],[17,13],[48,20],[28,25],[23,18],[16,14],[18,13],[34,26],[27,22],[22,14]], columns=['Husband','Wife'])\n",
"socks"
]
},
{
"cell_type": "markdown",
"id": "5906a5fc",
"metadata": {},
"source": [
"Let's plot those data. For paired data a scatter plot is often a good choice, but actually for this tiny dataset, I prefer showing the pairs using a plot like that shown on the right:\n",
"\n",
"
\n",
"
Why do I prefer the plot on the right? We are going to be interested in whether husbands have more socks than their wives or vice versa - I think this can be very clearly seen in the plot on the right (by inspecting whether the lines slope up or downwards) - however if there were 1000 couples in the sample rather than 10 this plot would be too crowded and hard to inspect\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Plotting code - don't get sidetracked by this, it's not that important\n",
"\n",
"plt.subplot(1,2,1)\n",
"sns.scatterplot(data=socks, x='Husband', y='Wife')\n",
"plt.plot([0,50],[0,50],'r--') # add the line x=y for comparison\n",
"\n",
"\n",
"plt.subplot(1,2,2)\n",
"sns.barplot(data=socks, ci=None, color=[0.9,0.9,0.9]) # ci=None switches off errorbars\n",
"for i in range(len(socks)):\n",
" plt.plot([0,1], [socks.Husband[i], socks.Wife[i]], '.-')\n",
" plt.xticks([0,1], labels=['Husband','Wife'])\n",
"plt.ylabel('pairs of socks owned')\n",
"plt.tight_layout\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "233aae63",
"metadata": {},
"source": [
"We can see that there is one couple in which the wife owns more socks and nine in which the husband owns more.\n",
"We also note that there is one couple in which the husband has an extreme number of socks.\n",
"\n",
"What is the mean difference in number of pairs of socks for [husband - wife]?"
]
},
{
"cell_type": "code",
"execution_count": 78,
"id": "0c6ac6e6",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6.6"
]
},
"execution_count": 78,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.mean(socks.Husband-socks.Wife)"
]
},
{
"cell_type": "markdown",
"id": "24f48262",
"metadata": {},
"source": [
"On average, the husbands own 6.6 more pairs of socks than their wives (but this will be skewed by the man with 48 pairs of socks).\n",
"\n",
"### Is this result statistically significant? \n",
"\n",
"That is, would the result (on average husbands own 6.6 more pairs of socks) be very unlikely to occur due to chance, if the null hypothesis were true, ie if there were no sex differences in the number of pairs of socks owned? \n",
"\n",
"To answer this question, we need to know what values for the mean difference in number of pairs of socks for [husband - wife] we would get due to chance - if actually all the people (male and female) were drawn from some distribution where the number of socks they owned does not depend on their sex. \n",
"\n",
"If we drew lots of samples of 10 couples from that sock-owning population, even though there is no overall difference in sock ownership between men and women (and therefore between husbands and their wives), nonetheless we would almost always get some difference between the husbands and their wives due to random chance. The distribution of these differences is called the null distribution of difference of means, that is, it is the distribution we would expect to obtain if the null hypothesis were true.\n",
"\n",
"### Obtaining the null distribution by permutation testing\n",
"\n",
"In previous weeks we have obtained simulated distributions of statistics such as the sample mean, by drawing many samples from a (known) parent population (as in the exercises on the Central Limit Theorem) or by bootstrapping. Here we will attempt to do something similar.\n",
"\n",
"We don't have access to the parent population, only the sample of 10 couples. The sample tells us several interesting things about the parent distribution, regardless of sex effects:\n",
"
\n",
"
for most couples there is a small difference in the number of pairs of socks owned\n",
"
for some couples there is a large difference, driven by one member having way more socks than one can imagine needing\n",
"
ie the distriution of differences has positive skew
\n",
"
husbands with lots of socks tend to have wives with lots of socks \n",
"
ie there is a correlation in number of socks between husbands and wives)
\n",
"
\n",
"\n",
"It also tells us about some potential sex effects:\n",
"
\n",
"
in most couples the husband has more socks\n",
"
the outlier individuals with loads of socks tend to be male\n",
"
\n",
"\n",
"What we are going to do is shuffle the data around to create many new (re)samples preserving the non-sex-related information but ignoring the sex of the sock owner. Using these simulated (shuffled) datasets we will work out how often we get a mean difference of 3.8 or more pairs of socks between husbands and wives, thus determining how likely our difference is to have occurred due to chance.\n"
]
},
{
"cell_type": "markdown",
"id": "974cdaf1",
"metadata": {},
"source": [
"### Run the simulation\n",
"\n",
"To generate new simulated datasets, we will shuffle around the datapoints in our original dataset. \n",
"\n",
"Which ones can we shuffle?\n",
"\n",
"
\n",
"
We assume there are no sex differences, so we can swap men for women\n",
"
We assume that it *does* matters which family you are in (high-sock men tend to have high-sock wives), therefore we cannot swap people between couples.\n",
"
\n",
"\n",
"Therefore, the only shuffling that we are allowed is to swap the labels 'Husband' and 'Wife' within couples. \n",
"To generate each new simulated dataset, we will randomly decide whether each couple from the original dataset gets flipped."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "2edd4003",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Husband
\n",
"
Wife
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
12
\n",
"
10
\n",
"
\n",
"
\n",
"
1
\n",
"
14
\n",
"
13
\n",
"
\n",
"
\n",
"
2
\n",
"
20
\n",
"
48
\n",
"
\n",
"
\n",
"
3
\n",
"
15
\n",
"
18
\n",
"
\n",
"
\n",
"
4
\n",
"
17
\n",
"
16
\n",
"
\n",
"
\n",
"
5
\n",
"
14
\n",
"
16
\n",
"
\n",
"
\n",
"
6
\n",
"
13
\n",
"
18
\n",
"
\n",
"
\n",
"
7
\n",
"
36
\n",
"
24
\n",
"
\n",
"
\n",
"
8
\n",
"
22
\n",
"
27
\n",
"
\n",
"
\n",
"
9
\n",
"
19
\n",
"
14
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Husband Wife\n",
"0 12 10\n",
"1 14 13\n",
"2 20 48\n",
"3 15 18\n",
"4 17 16\n",
"5 14 16\n",
"6 13 18\n",
"7 36 24\n",
"8 22 27\n",
"9 19 14"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"socks_shuffled = socks.copy() # work on a copy of the original dataframe\n",
"for i in range(len(socks_shuffled)):\n",
" if np.random.rand()>0.5: # generate a random number between 0 and 1 - if it is more than 0.5:\n",
" socks_shuffled.loc[i,'Husband'] = socks.loc[i,'Wife'] # flip number of socks for husband and wife\n",
" socks_shuffled.loc[i,'Wife'] = socks.loc[i,'Husband'] # flip number of socks for husband and wife\n",
" #else:\n",
" # don't shuffle the row!\n",
"socks_shuffled"
]
},
{
"cell_type": "markdown",
"id": "2bbe7db1",
"metadata": {},
"source": [
"#### What?\n",
"\n",
"The above might be clearer in an example where the flips are easier to see.\n",
"\n",
"Try running the code block below a few times and keep an eye on how the dataframe changes - note that in the original dataframe the man always has an odd number of pairs of socks."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "dd417139",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Husband
\n",
"
Wife
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
2
\n",
"
1
\n",
"
\n",
"
\n",
"
1
\n",
"
3
\n",
"
4
\n",
"
\n",
"
\n",
"
2
\n",
"
6
\n",
"
5
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Husband Wife\n",
"0 2 1\n",
"1 3 4\n",
"2 6 5"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df = pandas.DataFrame(data=[[1,2],[3,4],[5,6]], columns=['Husband','Wife'])\n",
"\n",
"df_shuffled = df.copy() # work on a copy of the original dataframe\n",
"for i in range(len(df)):\n",
" if np.random.rand()>0.5: # generate a random number between 0 and 1 - if it is more than 0.5:\n",
" df_shuffled.loc[i,'Husband'] = df.loc[i,'Wife'] # flip number of socks for husband and wife\n",
" df_shuffled.loc[i,'Wife'] = df.loc[i,'Husband'] # flip number of socks for husband and wife\n",
" #else:\n",
" # don't shuffle the row!\n",
"df_shuffled"
]
},
{
"cell_type": "markdown",
"id": "87067c37",
"metadata": {},
"source": [
"### Visualizing randoms shuffles in the sock data\n",
"\n",
"Back to our 'real' sock data\n",
"\n",
"Let's see how the distribution of differences changes over a few random shuffles.\n",
"\n",
"Below I generate 4 random shuffles of our sock data (in which some husbands and wives are randomly flipped), and plot the outcomes:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "c8a67dd9",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHVCAYAAAB8NLYkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAADHtUlEQVR4nOzdd3hb5fnw8e/RsCwveW/HI3b2ciY7zLB3ICQp0Bbetj86oLSlpVBGS1ltKdBNKS0jCYEQVpiBQNhJHGc4znISx9uWp+Qha53z/qFIsRM78dCy9HyuK3/gId1Gt47uZ5znlhRFURAEQRAEQRDGPFWgAxAEQRAEQRC8QxR2giAIgiAIIUIUdoIgCIIgCCFCFHaCIAiCIAghQhR2giAIgiAIIUIUdoIgCIIgCCFCFHaCIAiCIAghQhR2giAIgiAIIUIUdoIgCIIgCCFCFHaCIAiCIAghQjOcH37ggQd48MEH+30tLS2NxsZGABRF4cEHH+SZZ56hvb2dBQsW8Le//Y2pU6cO+TlkWaa+vp7Y2FgkSRpOeEIQURSFzs5OMjMzUalGNn7wR76ByLlQ4I18A3GNE4ZurFzjRL6FhmHlmzIM999/vzJ16lSloaHB889oNHq+/+ijjyqxsbHKa6+9ppSVlSlLlixRMjIyFLPZPOTnqKmpUQDxL0T+1dTUDCfF/J5vIudC699o8s1fOSfyLbT+Bfs1TuRbaP0bSr4Na8YOQKPRkJ6eftzXFUXhySef5J577uGaa64B4PnnnyctLY2VK1fy/e9/f0iPHxsbC0BNTQ1xcXHDDU8IEmazmZycHM/rOVK+zjcQORcKvJVvIK5xwtCMlWucyLfQMJx8G3ZhV1FRQWZmJjqdjgULFvDwww9TUFBAZWUljY2NLFq0yPOzOp2OhQsX8tVXXw2ahFarFavV6vnvzs5OAOLi4kQShoDRTv17O99A5Fwo88ZSk7jGCcMRbNc4kW+hbSj5NqyNAQsWLOCFF17ggw8+4N///jeNjY2cdtpptLa2evYEpKWl9fudvvsFBvLII49gMBg8/3JycoYTkhDCfJFvIHJOGJy4xgn+JPJN8IVhFXYXX3wx1157LdOnT+f888/nnXfeAVzTw27HVpOKopywwrz77rsxmUyefzU1NcMJSQhhvsg3EDknDE5c4wR/Evkm+MKojjuJjo5m+vTpVFRUePYIHDuSMBqNx404+tLpdJ4pYjFVLJyIN/INRM4JQyeucYI/iXwTvGFUhZ3VamXPnj1kZGSQn59Peno669ev93zfZrOxceNGTjvttFEHKggi3wR/Ezkn+JPIN8EbhnXzxM9//nMuv/xyxo0bh9Fo5KGHHsJsNnPzzTcjSRJ33HEHDz/8MEVFRRQVFfHwww8TFRXFsmXLfBW/EMJEvgn+JnJO8CeRb4IvDKuwq62tZenSpbS0tJCSksIpp5zCN998Q25uLgB33XUXFouF2267zXOY4ocffuiVIwiE8CPyTfA3kXOCP4l8E3xBUhRFCXQQfZnNZgwGAyaTSewNGMPG0us4lmIVBjaWXsOxFKswuLHyOo6VOIUTG87rKHrFCoIgCIIghAhR2AmCIAiCIIQIUdgJgiAIgiCECFHYCYIgCIIghAhR2AmCIAiCIIQIUdgJgiAIgiCECFHYCYIgCIIghAhR2IWJxu5GNjdsprG78eQ/LAijJPJN8DeRc8FBvA6BN6zOE8LYtLZiLQ989QAKCipU3H/a/VxTdE2gwxJClMg3wd+eK3uOJ0ufdOWcpOL+U0XOBcLairU8+NWDyMjivR9AYsYuxDV2N/Lg1w+i4GowIiPz4NcPitGU4BMi3wR/a+xu9BR1ALIici4Q3O99GRkQ7/1AEoVdiKs2VyMrcr+vyYpMTWdNgCISQpnIN8Hfqs3VnqLOTeSc/4n3fvAQhV2IGxc3Dgmp39dUkoqc2JwARSSEMpFvgr8NlFsi5/xvXNw4VFL/kkK8DoEhCrsQlx6dztSkqZ7/du8/SY9OD2BUQqhKj07n0oJLPf8t8k3wNYvT0u+/Rc4FRnp0Ovefen+/4u4nxT8Rr0MAiJsnQpyiKDT1NAFw17y7uCD3AvFGE3zKqTgBuHL8lfyo+Eci3wSf2tq0FYCZKTO5ffbt5MTmiJwLkGuKruG0zNP4/vrvc8h0iJSolECHFJbEjF2Iq+6sptnSjFal5boJ14kLnuBTiqKwtdH1QXvF+CtEvgk+V9JYAsCpmacyL32eyLkAS49OZ2HOQuDoayP4lyjsQpx7NDs9eTqRmsgARyOEutrOWowWIxqVhhkpMwIdjhDiFEWhpMlVPMxNmxvgaAQ392vh/vwR/EsUdiHOPWKakzYnwJEI4cD9ISsGEoI/1HbVYuwRA4lgU5xajEpSUd1ZjbHHGOhwwo4o7EKce8Q0N12MZgXfE7Mngj+5r2/Tkqah1+gDHI3gFhsRy8SEiYCYtQsEUdiFsPqueuq761FLamalzAp0OEIYcF/ExQyx4A9iRSJ4uV8Tsc/O/0RhF8LcH7JTk6YSpY0KcDRCqGvoaqCuq841kEidFehwhDAgViSCl/s1ETN2/icKuxDmXhYTo1nBH9z5NjlxMtHa6ABHI4S6xu5GartqUUkqsSIRhGanzgbgoOkgbb1tAY4mvIjCLoSJ0azgTyLfBH9y59vkxMnERMQEOBrhWAmRCRTGFwJQ2lQa4GjCiyjsQlRzTzNV5iokJLEsJviF2F8n+JNYkQh+nn12TWKfnT+Jwi5EuT9kJyVOIi4iLsDRCKGuxdLCYfNhJCSKU4sDHY4QBjwzxOIO7KAl9tkFhijsQpQYzQr+5M63CQkTMOgMAY5GCHUtlhYqTZVISMxOmx3ocIRBuIvufW37MFlNAY4mfIjCLkSJ0azgT+42YmJ/neAP7j1bRQlFYiARxJL1yeTF5aGgsN24PdDhhA1R2IWg9t52DnQcABCjWcEvxAyx4E8i38YOsc/O/0RhF4Lco9nC+EISIhMCHI0Q6jp6OzwDCfFBK/iDWJEYO9zXBLHPzn9EYReCxGhW8KetRtcFu8BQQGJkYoCjEUKdyWqior0CECsSY4G7+N7duptue3eAowkPorALQWI0K/iTyDfBn0qbSlFQyDfkk6xPDnQ4wklkxGSQFZOFU3Gyw7gj0OGEBVHYhRizzczetr2AmLET/EP06xT8SaxIjD1in51/icIuxGw3bkdBITcul5SolECHI4S4Tlsn+9r3AeKDVvAPMUM89rhfK7HPzj9EYRdixOyJ4E/bjNuQFZmc2BzSotMCHY4Q4rpsXexp2wOIa9xY4n6tylrK6HX0Bjia0CcKuxAjRrOCP4l8E/xpe/N2ZEUmOyab9Oj0QIcjDFFObA6p+lTssp2ylrJAhxPyRGEXQnrsPZS3lgNiNCv4h9jvJPiTWJEYmyRJOrrPrlHss/M1UdiFkO3N23EqTjKjM8mMyQx0OEKI67H3sLtlNyA6Tgj+4ZkhFvk25oi+sf4jCrsQIkazgj/taN6BQ3GQHp1OZrQYSAi+ZXFY2NW6CxDXuLHI/ZrtaN6B3WkPcDShTRR2IUSMZgV/6ru/TpKkAEcjhLqdzTtxyA7SotLIjskOdDjCMBUYCkjQJdDr7PVsGRJ8QxR2IaLX0evZlCpGs4I/iP11gj/1zTcxkBh7+u2zE+fZ+ZQo7EJEWUsZdtlOij6FcbHjAh2OEOKsTitlza6BhLgjVvAHsSIx9rlfO1HY+ZYo7EKEGM0K/lTWXIZNtpEUmURuXG6gwxFCnM1pY2fzTkDMEI9l7tduW9M2HLIjwNGELlHYhQhxnpjgT31nT8RAQvC1XS27sDqtJEYmkh+XH+hwhBEqii8iNiKWHkcP+9r2BTqckCUKuxBgd9o9zZXFaFbwB3dhJ/JN8Ie++SYGEmOXWqVmdupsQCzH+pIo7EJAeWs5vc5eEnQJjI8fH+hwhBBnl+1sb94OiBliwT/cRYDIt7HP/RqKws53RGEXAtyj2dlps8VoVvC5Pa17sDgsGHQGMZAQfM4hO9hm3AaIGeJQ4H4NS5tKkRU5wNGEJlHYhQAxmhX8yXOjTuocVJK4hAi+5R5IxEXEUZRQFOhwhFGanDQZvUaP2Wamor0i0OGEJHFVHuOcslOMZgW/EvvrBH/quyIhBhJjn0aloTi1GBDtxXxFvEvGuL3te+m2dxOrjWVCwoRAhyOEOKfspLSpFBDniQn+IVYkQo/YZ+dborAb47Y2ukY8xWnFqFXqAEcjhLr97fvpsncRo41hYsLEQIcjhLh+AwlR2IUM92z/1qatKIoS4GhCjyjsxjgxmhX8yZ1vxaliICH4XkVHBZ32TqK10UxMFAOJUDEteRo6tY623jYqzZWBDifkiMJuDJMVmVKjazQr9jsJ/iD21wn+5M63Wamz0Kg0AY5G8JYIdQQzUmYAYp+dL4jCbgw70HEAk9WEXqNnctLkQIcjhDhZkUW/TsGvShrFikSo8uyzaxT77LxNFHZjmGc0mzILrUob4GiEUHeo4xAd1g70Gj1TkqYEOhwhxCmKIlolhjD3rH9JU4nYZ+dlorAbwzyjWTF7IviBe3/dzJSZYiAh+Nwh0yHare1EqiOZmjQ10OEIXjYjZQYalQZjj5HartpAhxNSRGE3RvUdzYr9ToI/iHwT/MmdbzNTZqJVi4FEqNFr9ExLmgaIfXbeJgq7Meqw+TCtva1EqCKYnjw90OEIIU5RFHEHtuBX7hWJOeliIBGq3KtNYp+dd4nCboxyj3BmpMwgQh0R4GiEUFfdWU2LpcU1kEgRAwnBt8T+uvDQ9zw7wXtEYTdGeWZPxP46wQ/cI+rpKdPRqXUBjkYIdTWdNRgtRrQqrViRCGHFqcWoJBW1XbU0djcGOpyQIQq7MUhRlKPLFGK/k+AHYn+d4E/ufJuePJ1ITWSAoxF8JVobzeRE11FdYtbOe0RhNwbVddXR1NOERtIwM2VmoMMRwoDYXyf4kzvfxEAi9Im+sd4nCrsxyD2ymZo8Fb1GH+BohFBX31VPQ3eDGEgIfiP214UPsc/O+0RhNwaJ2RPBn9z5NiV5ClHaqABHI4S6hq4G6rrqUEtqZqXOCnQ4go/NTpuNhESlqZIWS0ugwwkJorAbg8R+J8GfRL4J/uQZSCSJgUQ4MOgMFCUUAVDaVBrgaEKDKOzGmKbuJmo6a1BJKopTiwMdjhAGRL9OwZ/EMmz4EfvsvEsUdmOM+6I3KXESMRExAY5GCHXGHiPVndViICH4jZghDj9in513icJujBH76wR/cl9oJyZMJDYiNsDRCKGuuaeZw+bDSEgUp4mBRLhwF3YV7RWYrKYARzP2icJujBGjWcGfRL4J/rTVeGQgkTiRuIi4AEcj+EuSPol8Qz4Kithn5wWisBtDWi2tHDIdAsQHreAfnv11osOJ4AdiP2f4EvvsvEcUdmNIqdE1kilKKMKgMwQ4GiHUtfW2cdB0EIDZqbMDHI0QDsQMcfgS++y8RxR2Y4injViquOgJvudeEimMLyQhMiHA0Qihrr23nQMdBwDX2WZCeHEXdnva9tBl6wpwNGObKOzGEM8xAGJZTPADMXsi+JN7RWK8YTyJkYkBjkbwt/TodLJjspEVme3N2wMdzpgmCrsxwmQ1sb99PyA+aAX/EHdgC/7kWZEQ17ew5X7t3bkgjIwo7MaIbcZtKCjkxeWRrE8OdDhCiDPbzOxr2weID1rBP8SKhOB+7cU+u9EZVWH3yCOPIEkSd9xxh+driqLwwAMPkJmZiV6v5+yzz6a8vHy0cYY9MZoV+eZP25pcA4ncuFxSolICHU5AiHzzn05bJ/vaw3sgIfLt6Gu/q3UXFoclwNGMXSMu7LZs2cIzzzzDjBkz+n398ccf54knnuCvf/0rW7ZsIT09nQsuuIDOzs5RBxvOwn00K/LNv8K9rZPIN//aZtyGrMiMix1HalRqoMPxO5FvLtkx2aRFpeGQHexs3hnocMasERV2XV1dLF++nH//+98kJBy9W05RFJ588knuuecerrnmGqZNm8bzzz9PT08PK1eu9FrQ4abb3s2etj1AeH7QinzzP/f+unCcPRH55n8i30S+AUiSdHSfnTjPbsRGVNj98Ic/5NJLL+X888/v9/XKykoaGxtZtGiR52s6nY6FCxfy1VdfDfhYVqsVs9nc75/Q33bjdpyKk6yYLNKj0wMdjt95M99A5NzJ9Nh72N26GwjPgYTIN/8L5xUJkW/9iX12o6cZ7i+8/PLLlJaWsmXLluO+19jYCEBaWlq/r6elpVFVVTXg4z3yyCM8+OCDww0jrITzaNbb+QYi507GPZDIjM4kIyYj0OH4lcg3/+ux97C7xTWQCLdrnMi347lzYGfzTmxOGxHqiABHNPYMa8aupqaG22+/nZdeeonIyMhBf06SpH7/rSjKcV9zu/vuuzGZTJ5/NTU1wwkpLITrfidf5BuInDsZzzEnYTZ7IvItMHY078ChOMiIziArJivQ4fiNyLeB5cflkxiZiNVpZVfLrkCHMyYNq7DbunUrRqOROXPmoNFo0Gg0bNy4kaeffhqNRuMZWbhHGm5Go/G4UYebTqcjLi6u3z/hKIvDQllLGRB+hZ0v8g1Ezp1MuB5MLPItMMJ1RULk28DEPrvRG1Zhd95551FWVsb27ds9/+bOncvy5cvZvn07BQUFpKens379es/v2Gw2Nm7cyGmnneb14MNBWXMZDtlBalQq2bHZgQ7Hr0S++V+vozdsBxIi3wIjXFckRL4Nzp0LYp/dyAxrj11sbCzTpk3r97Xo6GiSkpI8X7/jjjt4+OGHKSoqoqioiIcffpioqCiWLVvmvajDSN/R7Imm30ORyDf/K2spwy7bSdGnkBObE+hw/Erkm/9ZnVbKml0DiXCbsRP5Njh3LmwzbsMu29GqtAGOaGwZ9s0TJ3PXXXdhsVi47bbbaG9vZ8GCBXz44YfExsZ6+6nCQriOZodK5Jt39W0jFm4DiaEQ+eZdZc1l2GQbyfpkcuNyAx1O0AnXfCtKKCIuIg6zzcze1r1MT5ke6JDGFElRFCXQQfRlNpsxGAyYTKYxtzfA22xOG6etOg2r08qbV75JQXxBoEMasrH0Oo6lWH3t1g9uZVPjJu5dcC9LJi0JdDhDNpZew7EUq6/9c8c/+dv2v3Fh3oX8ceEfAx3OsIyV13GsxHmsH2/4MZ/WfMqdc+7kO9O+E+hwAm44r6PoFRvEylvLsTqtJEYmkm/ID3Q4QoizO+3saN4BhN8dsUJgiBUJYTBin93IicIuiPXtDyuWxQRfK28tp9fZS4IugQLD2JkdFsYmu3x0IBFu++uEk3MXdqVNpThlZ4CjGVtEYRfEwvXYCSEwwvlGHcH/drfuxuKwEK+LZ3z8+ECHIwSZiYkTidZG02nvpKKjItDhjCmisAtSDtnBNuM2QCxTCP4RrueJCYHhXpGYnToblSQ+ioT+NCoNs1JnAUdzRRga8W4KUnvb9tLj6CE2IpaihKJAhyOEOIfsYLtxOyD21wn+Ec79YYWhEfvsRkYUdkHKs78udY4YzQo+t69tH932bmK1sRTFi4GE4FtO2elZkRAzxMJg+hZ2QXaAR1ATFUOQEqNZwZ/cy7Cz02ajVqkDHI0Q6va176PL3kWMNoaJCRMDHY4QpKYmTSVSHUm7tZ1DpkOBDmfMEIVdEHLKTrYaxY0Tgv+I/XWCP7lXJIpTi8VAQhiUVq1lZspMQOyzGw5R2AWhAx0H6LR1EqWJYlLipECHI4Q4WZEpbSoFxI06gn+IFQlhqOakuwabYp/d0InCLgi5Z0+KU4vRqLze9U0Q+qlor8BsM6PX6JmUJAYSgm/JiixWJIQhcw82S5pKxD67IRKFXRASo1nBn9z5VpxaLJptCz53sOMgJqsJvUbPlKQpgQ5HCHLTk6ejVWlptjRT01kT6HDGBFHYBRlFUcTBxIJfif11gj+5821mykwxkBBOKlITyfTk6cDR3BFOTBR2QabSVElbbxs6tY5pSdMCHY4Q4voOJMT+OsEfRL4Jw+UedIp9dkMjCrsg0280qxajWcG3Ks19BhLJYiAh+JaiKP16YAvCUHj22Yk7Y4dEFHZBxl3YidGs4A/uEfCMlBlEqCMCHI0Q6qrMVbT2thKhimB6yvRAhyOMEbNSZ6GW1NR311PfVR/ocIKeKOyCiKIobG0U++sE/xGzJ4I/uQeu01Omo1PrAhyNMFZEaaM8N9qI5diTE4VdEKntrMVoMaJRaZiRMiPQ4QghTlEUMUMs+JXYXyeMlOgbO3SisAsintFs8nQiNZEBjkYIdbVdtRh7xEBC8I++AwkxQywMlztnxJ2xJycKuyAiZk8Ef3KPfKclTUOv0Qc4GiHU1XfX09jdiEbSeNpECcJQFacVIyFRZa6iuac50OEENVHYBRFxfp3gT2J/neBP7nybkjyFKG1UgKMRxpq4iDgmJk4ExHLsyYjCLkg0dDVQ11WHWlIzK3VWoMMRwoDocCL4k9hfJ4xW3/ZiwuBEYRck3Ik6OXEy0droAEcjhLrG7kZqu2pRSSpmpcwKdDhCGBD764TREgcVD40o7IKEmD0R/Mmdb5MTJxMTERPgaIRQ19TdRE1nDSpJRXFqcaDDEcao2WmzATjQcYD23vYARxO8RGEXJMT+OsGfxOyJ4E/u69vEhInERsQGOBphrEqMTGS8YTwApU2lAY4meInCLgi0WFo4bD6MhCRGs4JfiP1Ogj+JFQnBW9w5JPbZDU4UdkHAnaATEiZg0BkCHI0Q6losLVSaKpGQPEsbguBLYoZY8Baxz+7kRGEXBNxtxMRoVvAH9xJGUUKRGEgIPtdqaeWQ6RAAc1JFYSeMjruw29u2l05bZ4CjCU6isAsCYjQr+JPIN8GfSo2ugURhfCHxkfGBDUYY81KjUhkXOw4FhW3GbYEOJyiJwi7AOno7ONBxABAftIJ/iP11gj+JfBO8TeyzOzFR2AXYVqProldgKCAxMjHA0QihzmQ1UdFeASD21wl+4elwki4GroJ3ePbZNYp9dgMRhV2AidGs4E+lTaUoKOQb8knWJwc6HCHEmawm9rfvB8Q1TvAedy7tbt1Nj70nwNEEH1HYBZg4v07wJ5Fvgj9tN25HQSEvLk8MJASvyYzJJCM6A4fiYEfzjkCHE3REYRdAnbZO9rbtBcQHreAf7j0pYvZE8Adxo47gK6Jv7OBEYRdA243bkRWZnNgc0qLTAh2OEOK67d3sadsDiA9awT/EDLHgK+I8u8GJwi6AxOyJ4E/bjNuQFZnsmGzSo9MDHY4Q4rrt3exu3Q3AvPR5AY5GCDXuO2PLmsuwOq0Bjia4iMIugMRoVvAnkW+CP+0w7sCpOMmKyRIDCcHrxsWOI1mfjE22UdZcFuhwgooo7AKkx95DeUs5IDpOCP7hPnZC5JvgD2J/neBLkiSJfXaDEIVdgOxs2YlDcZAenU5mdGagwxFCnMVhYVfrLkB80Ar+IY5yEnxN7LMbmCjsAsQze5I2F0mSAhyNEOp2Nu/EITtIi0ojOyY70OEIIa7X0UtZi2t5TBR2gq+4c2tH8w7ssj3A0QQPUdgFiNjvJPhT33wTAwnB18payrDLdlL1qWTHioGE4BsF8QXE6+KxOCyeG3UEUdgFhNVpZWfzTsB/o1mTyURlZSUmk8kvzycEF88d2H7aXyfyLbz1bSPmr4GEyLnwo5JUnskRd875SzDnmybQAYSjXS27sMk2kiKTyI3L9fnzlZaW8vbbb6MoCpIkcfnllzN7tugTGi5sTptnIOGPGWKRb4K/99eJnAtfc9Lm8HH1x2xt2sot02/xy3MGe76JGbsA6Ht3oq9HsyaTyZOAAIqi8PbbbwflKEPwjV0tu7A6rSRGJpIfl+/T5xL5Jtiddk+bJ38UdiLnwps7x7YZt+GUnT5/vrGQb6KwCwB/7q9ra2vzJKCboii0tbX5/LmF4ODP/XUi34Ty1nJ6nb2ugYTBtwMJgKamJpFzYWxCwgRitbF02bvY177P58/X2toa9PkmCjs/s8t2tjdvB/wzmrVajz+RW5IkEhMTff7cQnDwZ4cTvV5/3NdEvoWXvufX+XogIcsymzZtOu7rIufCh1qlpjitGPDPPruqqqrjvhZs+SYKOz/b07oHi8OCQWdgfPx4nz5Xe3s7b731Vr+vufcDGAwGnz63EBwcsoNtxm2A72eIHQ4H77zzTr+viXwLP/48mPijjz7i4MGDSJLkKSJFzoUff51nt2/fPjZu3Njva8GYb+LmCT/zXPRS56CSfFdXWywWVqxYQU9PD+np6Vx77bV0dXWRmJgYVAko+JZ7IBEXEUdRQpHPnkdRFN58801qamrQ6XQsWbLEM4oV+RY+HLKDbU2ugYSvZ4i3bt3KV199BcDVV19Nbm4ubW1tIufCkDvXthq3IiuyTz5bGxoaWLNmDQBz5szhrLPOCtp8E4Wdn/ljf53T6eSVV16hpaWF2NhYli1bRlxcHCkpKT57TiE4ufNtdtpsnw4kPv30U8rKylCpVFx//fUUFBT47LmE4LWvbR89jh5iI2IpjC/02fMcPHiQdevWAXD22WczY8YMgKD7gBX8Y3LSZPQaPSariYMdB70+iDWZTKxcuRK73U5BQQGXXHIJarU6aPNNLMX6kVN2UtpUCrjOd/IFRVF45513qKysRKvVeoo6ITz5Y3/djh07PMsTl156KePH+3aLgRC83Pk2O3U2apXaJ89hNBp55ZVXUBSFGTNmsHDhQp88jzB2aFVaZqbMBLzfN9ZqtbJq1So6OztJSUnh+uuvR632TW57iyjs/Gh/+3667F1Ea6OZlDDJJ8/x1VdfUVpaiiRJLF68mIyMDJ88jxD8+g4kfFXYVVVVefZxnn766cyZIzqphDNfDyS6urpYsWIFVquVcePGccUVV4hOKgLQZznWi/vsZFnmtddeo7GxkejoaJYtW0ZkZKTXHt9XRGHnR+6LXnFqsU9Gs7t372b9+vUAXHjhhUycONHrzyGMHRUdFXTaO4nWRjMx0fu50Nrayssvv4zT6WTy5Mmcd955Xn8OYeyQFfnoioQPtprY7XZWrVqFyWQiMTGRG264AY1G7CYSXPp2oDj2OJKR+uCDD9i/fz8ajYYbbriBhIQErzyur4nCzo98eRp7XV0da9euBWDevHksWLDA688hjC3ufJuVOguNyrsfgD09PaxcuRKLxUJmZiZXX301KpW4nISzivYKzDYzeo2eyUmTvfrYsizz+uuvU1dXR2RkJMuWLSMqKsqrzyGMbdNTphOhiqC1t5Uq8/FHkgzX5s2bPUfpXH311eTk5Iz6Mf1FXIn9RFZkn9040dHRwcqVK3E4HBQWFnLRRReJ5QnhaIcTLw8kHA4Hq1evprW1FYPBwNKlS4mIiPDqcwhjT98VCW8PJD755BN2796NSqXihhtuIDk52auPL4x9OrWO6SnTgdHvs9u/fz/vvfceAOeddx5Tp04ddXz+JAo7PznUcYgOawd6jZ6pSd5Lkt7eXlauXEl3dzdpaWlcd911Qb+xU/A9RVF8MkPsbp9TVVVFREQEy5YtIzY21muPL4xdvlqR2LZtG59//jkAV1xxBXl5eV59fCF0eGOfXWNjI2vWrEFRFIqLiznjjDO8FZ7fiMLOT9wjiBkpM9CqtV55TKfTyZo1azAajcTExLBs2TJ0Op1XHlsY2w6ZDtFubSdSHenVgcTnn3/Ojh07kCSJ66+/nrS0NK89tjB29R1IeHNForKykrfffhuAs846i1mzZnntsYXQ49ln1zSyfXadnZ2sXLkSm81GXl4el1566Zhc/RKFnZ94ezSrKArvvfceBw4cQKPRsHTp0qA9U0fwP3e+zUyZ6bWBxK5du9iwYQMAl1xyCYWFvjunTBhbKs2VtPW2oVPrmJY8zSuP2dzczOrVq5FlmalTp3L22Wd75XGF0DUzZSYaSUNjdyP13fXD+l2bzcbKlSsxm80kJSWxZMmSMXtzjijs/EBRFK+32fnmm28oKXE95rXXXktWVpZXHlcIDe79dd7Kt5qaGl5//XUATjnlFObNm+eVxxVCgzvfZqTMIEI9+v2W3d3drFy5kt7eXrKzs7nqqqvEzTnCSUVpo5iSPAUYXt9YWZZZu3YtDQ0NREVFsXz58gH7Xo8V4p3iB9Wd1bRYWtCqtMxImTHqx9u7dy8ffPABAIsWLWLyZO/egSaMbf3216WPfoa4ra2NVatW4XQ6mThxIosWLRr1YwqhxZsrEna7nZdffpn29nbi4+NZunQpWq13Zp2F0DeSfXYfffQRe/fuRa1Wc8MNN5CYmOir8PxCFHZ+4B45TE+ejk49uj1w9fX1vPbaa4CrX92pp5466viE0FLTWYPRYkSr0jI9efqoHstisbBy5UpPz+FrrrlGzJwI/XhzRUJRFN566y1Pz+Hly5cTHR3tjTCFMNF3n91QlJSUeHoOX3nllYwbN85nsfmLuEL7gbdmT0wmE6tWrerXr24sbuwUfMudb9OTpxOpGfkp6QP1HBY35wjHqu2qxdhjRKPSjHpFom/P4SVLloj+1sKwFacWo5JU1HTW0NTddMKfPXDgAO+88w7Qv+fwWCcKOz/wxmjWarWycuXKMdWvTggMb+SboiisW7dO9BwWTsq9IjEtaRp6zcj3JfXtOXzZZZdRUFDglfiE8BIbEcvEBFennRMtxxqNRl599dWQ7DksCjsfq++qp6G7AY2kYVbKrBE9hrtfXVNT05jqVycEhjf2O3355Zds27YNSZK47rrrRM9hYVDeWJE4tufw7NmzvRKbEJ7cuThYYRfqPYdFYedj7tmTKUlTiNKOrAVO3351S5cuHTP96gT/a+hqoK6rDrWkZmbqzBE9xu7du/noo48AV8/hCRMmeDNEIcSMdoZY9BwWvO1E++zCoeewKOx8zHNoZ/rILnqbNm3q168uOzvba7EJoafvQCJaO/xN57W1tZ6ew/Pnz+eUU07xanxCaGnsbqSuqw6VpKI4tXjYv9/T08OKFSuwWCxkZWWJnsOCV8xJdX3eHjIdotXS6vl6357Der0+ZHsOi3eQj42mX+f+/ft5//33gbHZr07wv9Gc/t/R0cGqVatwOBwUFRVx4YUXejs8IcS4BxKTEycPeyDh7jnc1taGwWDghhtuED2HBa+Ij4ynMN51gHqpsdTz9Q0bNnh6Di9ZsiRkew6Lws6HjD1GqjurkZCGPZoNhX51gv+NdH/dsT2HFy9eLG7OEU5qpPkmeg4LvnbseXbbtm3jiy++AEK/57Ao7HzInVCTEicRGzH0i5bZbA6JfnWCfzX3NHPYfNg1kEgb+kDC6XTy6quvip7DwrCNtMPJZ599JnoOCz7l3v5U0ljCoUOHwqrncGjtGAwyI1kWs9lsrFq1CrPZTHJy8pjuVyf411ajK98mJk4kLmJoR5O4ew4fPHjQc6yJ6DksDEWLpcUzkJidNvS7WMvKyvjkk08A0XNY8B33jF19Uz2rd7t6Dk+bNo1zzjknwJH5nqgYfGi4++uO7Ve3bNmyMd2vTvCvkcye9O05fM0115CZmemT2ITQ4x64FiUUYdANbTBQXV3NG2+8AcCpp54qeg4LPpOsT6YwqpDxNeOxOqxkZ2dz5ZVXhsXqlyjsfKStt42DpoMAQx7Nrl+/PqT61Qn+Ndz9TqLnsDAaw823trY2z7EmEydO5IILLvBleEKYs9vtzKybCQ6QoqSw6jksCjsfKW1y3YlTGF9IQuTJz50rKSnh66+/BuCqq64KiX51gv+097ZzoOMAMLSBhOg5LIzWcM6v69tzOCMjg2uvvVYcayL4jKIovPnmm9ABNpWNuvy6sOo5LAo7HxnO/rq+/erOOeccpk8fXeN2Ify4b+kfbxhPYuSJZ3pNJhMrV67Ebrczfvx40XNYGDaT1URFewVw8mtc357DcXFxLF26VBxrIvjUp59+yq5du1CpVHyT+g1tlja67d0jOttzLBJDJh9xj2ZPtkzR1NTUr1/dWWed5Y/whBAz1P117p7DXV1dpKSkcN1114ljTYRhcw9c8w35JOmTBv25vj2H3ceaiJ7Dgi8d23M4IiUCp+Jku3F7YAPzo2EVdv/4xz+YMWMGcXFxxMXFceqpp/Lee+95vq8oCg888ACZmZno9XrOPvtsysvLvR50sDPbzOxr2wec+IO2q6uLlStXYrVayc3NDbl+daMl8m3ohtKv0+l0smbNGtFz+AREzg3NUPfX9e05vHjxYtLT0/0R3pgh8s27Dh8+7FqC5WjPYfdn8GB9Y0PRsAq77OxsHn30UUpKSigpKeHcc8/lyiuv9CTa448/zhNPPMFf//pXtmzZQnp6OhdccAGdnZ0+CT5YbWvahoJCblwuKVEpA/7Msf3qxLEmxxP5NjSdtk72tZ98IPHhhx9SUVEheg6fgMi5oRnK/rry8nJPz+GLLrpI9BwegMg372ltbWX1atexJlOmTPH0HHYPPgbqGxuylFFKSEhQnn32WUWWZSU9PV159NFHPd/r7e1VDAaD8s9//nPIj2cymRRAMZlMow0tYP605U/KtP9NU+7/8v4Bv+90OpXVq1cr999/v/Loo48qLS0t/g3QD3z1Ono733wZq79srNmoTPvfNOWS1y4Z9Ge++eYb5f7771fuv/9+ZdeuXX6Mzj98+RqKa1x/ndZOZcbzM5Rp/5umNHQ1DPgzNTU1yu9+9zvl/vvvV9555x0/R+gfY+UaN9bzbSi6u7uVp556Srn//vuVZ555RrHZbJ7vVZuqlWn/m6bMemGWYrFbAhjl6AzndRzxHjun08nLL79Md3c3p556KpWVlTQ2NrJo0SLPz+h0OhYuXMhXX3016ONYrVbMZnO/f2PdyUazx/arS0oafI+K4OKtfIPQy7mT5VvfnsPnn3++6Dk8ROIaN7Btxm3Iikx2TDbp0ccvrR7bc/iiiy4KQJRjj8i3kTm25/Cxx5pkx2aTqk/FITvY2bwzgJH6z7ALu7KyMmJiYtDpdPzgBz/g9ddfZ8qUKTQ2NgIc1xomLS3N872BPPLIIxgMBs+/nJyc4YYUVHrsPexu3Q0MvP+ktLQ0bPrVeYO38w1CL+dOtL/u2J7Dp59+ur/DG3PENe7ETpRvvb29rFixol/PYXGsyYmJfBs5RVF46623qKqqQqfTsWzZMmJiYvr9jCRJnvZi4bLPbtjvuIkTJ7J9+3a++eYb/u///o+bb76Z3bt3e75/7OZ/RVFOeEPA3Xffjclk8vyrqakZbkhBZbtxO07FSWZ0JhkxGf2+d+jQIdatWweER786b/B2vkFo5VyPvYfdLa7/H8fO2PXtOZyfn89ll10mbs4ZAnGNO7HBZojdPYebm5tFz+FhEPk2cp999hk7d+5EkiSuu+66QXsOh9s+u2Hv1o+IiPD09ps7dy5btmzhqaee4pe//CXgmiHIyDha0BiNxhM2eNbpdCH15vccc3LMaLa5uZlXXnklrPrVeYO38w1CK+d2NO/AoTjIiM4gKybL8/Vjew5ff/314liTIRLXuMFZHBbKW1wb+/uuSCii5/CIiXwbmb49hy+99NIT9hx25+qO5h3YnXa06tDuQDHqOXJFUbBareTn55Oens769es937PZbGzcuJHTTjtttE8zZgx0MHF3dzcrVqygt7eXnJycsOlX5wsi3/obaPZElmVee+010XPYS0TOHeUeSKRFpfUbSHz99deensPXXnut6Dk8CiLfTu7YnsNz55742J18Qz6JkYlYnVZ2te7yQ4SBNawZu1//+tdcfPHF5OTk0NnZycsvv8ynn37K+++/jyRJ3HHHHTz88MMUFRVRVFTEww8/7PlgCQe9jl7KWsqAoyMEu93Oyy+/TEdHBwkJCdxwww1h069utES+ndxA54mtX7+effv2iZ7DIyBy7sT67q9zD0737t3Lhx9+CLh6Dk+aNClg8Y01It+GbyQ9hyVJYk7aHNZXrWdr01aKU4v9EGngDKuwa2pq4sYbb6ShoQGDwcCMGTN4//33Pf9j77rrLiwWC7fddhvt7e0sWLCADz/8kNjYWJ8EH2zKWsqwy3ZS9CnkxOZ4+tXV1NQQGRnJsmXLwqpf3WiJfDsxq9NKWbNrIOGesRM9h0dH5NyJHdvhpG/P4blz54qew8Mk8m14RtNz2F3YlTSWcOv0W30caWBJiqIogQ6iL7PZjMFgwGQyjbnWM//Y8Q/+vv3vXJx3MY8vfJwNGzbw2WefoVKp+Na3vkVBQUGgQ/SbsfQ6jqVY+yppLOE7H3yHZH0yG67bwMGDB1mxYgWKonDOOeewcOHCQIfoN2PpNRxLsfZlc9o4deWp2GQbb131Fokk8u9//5uuri7Gjx/PsmXLwmof51h5HcdKnCfjcDhYsWIFlZWVxMXFceuttw7r79nXto/Fby8mShPFl0u/RKMaWw0BhvM6ivvQvWhr49H9ddu3b+ezzz4DXP3qwqmoE/yj7/46o9HIK6+8gqIozJw5U/QcFryurKUMm2wjKTKJDF2G6Dks+I2iKLzzzjuj6jlcGF9IbEQsPY4e9rbt9VGkwUEUdl5id9rZ0bwDgBxHDm+99RYAZ5xxBrNnzw5kaEKIcu93mhU3y3OsSW5uLpdffrm4OUfwOs+NYSlzeO211zw9h5cvXy56Dgs+5Y2ew2qVmjmp4XGenSjsvKS8tZxeZy+ZUiafv/O5p1/dueeeG+jQhBBkl10DCbWspv2bdtFzWPA59/66nIacfj2H4+PjAxuYENK82XPYsxe5MbTPsxOfAF5S0lRChDOC+c3zsVgsZGVlcfXVV4tT1wWf2N26G4vdwpltZ9La2Yper2f58uVERUUFOjQhBNllO9ubtzPeNJ6Otg4Arr76arKzswMbmBDSamtref311wGYP38+CxYsGNXjuc+X3WrciqzIqKTQ/HwOzb8qAEoaSjil6RTUFvWA/eoEwZtKGkuY1j6N1M5U0XNY8Lk9rXswmA3MbJsJiJ7Dgu+1t7d7vefwpMRJRGmi6LR1UtFe4YUog5Mo7LzA7rTjLHOSYk1BG6Fl+fLlx/WrEwRvKt9RzkTTRACuvPJK0XNY8Kmv933NAuMCJCTRc1jwud7eXlauXOn1nsMalcZzhl0otxcThZ0XvPHhG2R1ZiEjc91115GamhrokIQQduDAAWIOuAYOU+ZNYebMmQGOSAhlZrOZ2s9q0SgaIlMjRc9hwad83XPYvc8ulG+gEIXdKJWVlVG+ydU7saewhwlFI9/Y6Uv2xka6v9mEvbEx0KGEvdG8Fs3Nzax+ZTUqVDTENnDNRdf4IEIhlIwm32w2GytXrURtU2PWmrnwygvFsSbCCY0m3xRF4d133/Vpz2HPPrumrQTZMb5eIwq7Uejbr25/3H6mzgrOPScda9Zw4NzzqP72tzlw7nl0rFkT6JDCVseaNRw459wRvRbunsN2m50WXQvq6Wo0anH/kzC40eSbu+dwY0MjvapetmVtY3rGdB9GK4x1o/2s+frrr9m61TWT5quew1OTpqJT62jrbaPSVOn1xw8GorAbob796owxRsoSy/r16wwWPdu20fCb+0CWXV+QZRruu1/M3AWAvbGRhvvuB/cocRivRd+ew85IJ1+nfc3czODLNyF4jCbf4GjPYUkl8XXa10zKnjTmTusX/MeTb30/a35zH7ba2iH9/p49ezw9hy+88EKf9RyOUEcwM8W1fSVU99mJwm4ELBYLK1asoKenh8TURL5K+gq9Vs+kpOBofq0oCj1btlD74x9TtWz50Qu7myxjq6oOTHBhzHa46uhFz02WqfvFXfSUbht0WUCW5X49hzdlbMKmtnn2igjCQAbLt+pbbsX83nsoDsegv7tlyxZPz+GuiV20RbaJfBNOaMB8UxQOL7mB1v/+D2dn56C/W19fz9q1awFXz+FTTjnFl6EePc9OFHYCuPrVrV69mtbWVuLi4og/JR6nyklxajFaVWCPN1FsNkxvvsnhaxdTdeNNdK7/6PiiDkClIiJXNIf3t4i8XBjgzi7Lli1ULVvG4euXYFr3Dord3u/7n376Kbt27UKlUnH6JafTQAN6jZ4pSVP8FbowBg2Wb7aDB6n76Z0cuGARrc8+i9Nk6vf9AwcO8O677wJwzjnn8I3yDUBQrkgIwWOwfHO2tmJ87DEOLDybxod+j62qqt/3TSYTK1euxG63M378eC6++GKf35zjzuVQ3WcnCrthcPerO3z4sKdf3XbzdoCAjmYdbW20/OMfVJx3HvW//BW9u3cjRUYSv2QJBeveJuOh3x19w6lUZPz2QbQjaMkijI42PZ2M3z7Y77VI/tGPMCy+Fikigt6yMup//nMOnH8BLf96Bkd7e7+ew5dffjm1Ea5ljZkpMwM+kBCC20D5lvqrX5L8wx+iTkrC0dCA8Y9/ouLsc2j87W+xVlbS1NTUr+dw5vRMOqwdRKojmZoUnHuIheAwUL6l3X8f6b/7LbqiQuSeHtpfeomDF11MzW0/pHvTZs+xJl1dXaSmpvqt5/D0lOloVBqMPUZqO4e2VDyWiA0Tw/DFF1/061eXlpbG1o2ujZ6BGM327t9P+4svYnrrbRSrFQBNaioJy5cTf/11aBISANAVFhJ9xhnYqqqJyB0niroAil+8eMDXIvXOO+lYvZq2lStxNDXR/Oc/s/uVV/j0jDNAcvUcLi4uZsXGFYCYPRGGZrB8S/re/8P8zru0Pf881n37aF+5ivq1r/PxZZdh06g9PYdfO/AaADNTZ6JVi4GEcGKD5Vv84sX0fP01rc8/T/fGz+jasAHzJ5/w1UUX0WSIIzo6mmXLlvmt57Beo2d68nS2GbdR0lRCTlyOX57XX0RhN0Tl5eV8/PHHwNF+dYdMh2jrbUOn1jEteZpf4lBkme7PP6ft+efp/uprz9cjp00j8eabibvoQqQBOl5o09NFQRckBnotNImJJP/f/5F4yy10vvcelatW8UV+PrIEOdXVFK1+hU6djpKGLUBgZ4iFsWWgfFPpdMRfczWGq6+iZ9NmjM8/z3qNmm6NmhhzJwvefIuu2Fi2GTYBIt+EoRso3yRJIvq004g+7TSshyppf+lFPj5wgDpDHGqHg9M3fIJdH4XjhiVokpP9EuectDmewu7qoqv98pz+Igq7IaipqfH0q1uwYIGnX537gMMZKTOIUEf4NAa5p4eON96g/YUXsR0+7PqiSkXsBReQePNN6IuLT7gvobe3gR7LYaL0eURGZvg0VuHETvRaqCIi0Jx/Pp9VVWFrayPFbmf+5i30OBz0fP459yTBB/MjmHptYYCiF8aaE+WbJEno589jc9Vh2vbsQQcs3LIFqbmZhnt/wzVRKqKLFebNFfkmeIeuIJ/DixZRceRGi9P27iOhspKWv/6V1n/9i7jLLyfxphuJ9NFdsW5z0+bybNmzIz6oOJg/U0VhdxLt7e28/PLLOBwOJkyYwIUXXuj5Xkmj644aX45m7Q0NtK9YQfsrryKbzQCoYmOJv+46EpcvQ5uVddLHqK9/hT177wFkQMXkSb8nM/N6n8UsDO5kr4X75py2tjbi4+O5+dZbifjBD2h/6SWaX1lNdquVW96zUfP1xcQvWULCsqVo09IC9vcIwW0o7/2PP/6YPXv2oFarWXbTTWTfeScda16j+cX/EdtgZPGXwKY7qb/0UhJvvonIKeKmHWHk9u3bxwcffAC4eg6ffs89dK5fT9vzL2DZsQPT2rWY1q4lasECEm++mZizFyJ5oZ3YsWalzkItqanrqqOxu5H06KGvaAX7Z6q4eeIE+varS09P59prr/X0q1MUxXOrtC/2O1m2b6fuzjs5cP4FtD77H2SzGW3uONLuvZfCTz4h7a5fDKmo6+1t6JOAADJ79t5Db2+D12MWTuxkr4WiKLz11ltUV1ej0+lYtmwZMTExROTkkHb33bzy2Pn87zwVltQ4nB0dtP7rXxw473zqfv4LLGW7AvZ3CcFpKO/9rVu38uWXXwJwxRVXkJubizoujqTvfofd//gRf7paRU1+DDgcmN58k8prrqXqWzfS+dFHKE6n//8oYUxraGhgzZo1KIrC7NmzOf3005G0WuIuuYS81S+T9/Iq4i65GNRqejZtova22zh48cW0vbQCubvbq7FEa6OZnDgZGN6xJ2PhM1UUdoNwOp288sorNDc3Exsby9KlS/v1q6vtqsXYY0Sj0jAjZYZXnlOx2zG/+y6VS5Zw+IalmN99D5xOok45hey//53x771H4reWo46JHvJjtrV9zdEEdJOxWKoG+nHBh3oshznRa7Fx40Z27tyJJElcf/31/XoOK4rC1+YdvDtfhXXFE2T95Wmi5s4FhwPzunUcvu46Di9bjvn9D054PpkQPk6Wb4cOHeKdd94BYOHChcf1HC5p3samSSr2/v5G8l59hbjLLgONhp6SEmp/9GMOXngRbc8/j7Ory/d/jDDmmc1mz7Em+fn5XHrppcdtH9LPmkXWE09Q+NF6km69BVVcHPaqapoeeoiKs8+h6bHHsdfVeS0mz3l2jUMv7BobXyfYP1NFYTcAd7+6Q4cOodVqWbp06XH96tzr8tOSpqHX6Ef1fE6TidZnn+XABYuou/Nn9O7YiaTVYrjmGvLfeJ3c//2X2HPPGdZ0tKIo1NevYe++3wzwXRV6fe6oYhaGL0qfx0Bvua6ufezcuZNPP/0UgEsvvZTx48f3+5n67noauxvRSBpmphcTd8EF5L70InmvrcFw5RWg1WIpLaXujjs4uOhCWp/77wkPBBVC32D51ti0jqamelavXo0sy0yfPp2zzz77uJ9zX+Pmps9FP306WX/8A4Uff0TS976H2mDAXltL0yOPcmDh2TQ98gi2mhrf/kHCmGW1Wlm5ciWdnZ0kJydz/fXXn/BYE21GBqk//zlFn35C2n2/ISIvD7mzk7b//pcDFyyi9vY76CktHfUZdH37xp6MLFvZv/93HDz0pwG+G1yfqaKwG8BQ+tW5E2E0++ushyppePBBKs4+B+Mf/4SjsRF1UhLJP/oRhZ9sIPPh349oA6nD0Un57jvZs/eXKEovUVHjOfpSu/YDBNtmz3AQGZnB5Em/59i33f6K31Je/nNUKjunnXYac+cev7TvzrcpyVOI0kZ5vq6fOpXMxx6j8OOPSL7t/1AnJGCvr8f4+OODHggqhIfj8801O1Jfv4rNWxYjSS3k5ORwxRVXHDdzUt9VT313PWpJzayUWZ6va9PSSL3zpxR++gnpDz5IxPjxyN3dtD3/AgcvvIjaH/+Yni1bQvLQV2FkZFlm7dq1NDY2EhUVxfLly9HrhzYZooqKInHZMgrefYecf/2T6NNOA1mm84MPqFq23HWo+9vrjjvUfaiKU4uRkDhsPkyLpWXQn+vpqaSk5Dpqav8HQELC6QTzZ6q4eeIYQ+1X5566dVf8Q6UoCt1ffUXbCy/QvfEzz9d1kya5jiu59BJUESO/w9Zs3smu8tuxWKqRJDUF+T8lN/d7WK1GLJYq9PrcoErAcJOZeT2JiWdisVSh02VTWfkSDY3Pkpp2kMQkEwsWLB3w9zz5Nsh+Tm1qKik/+QlJ3/se5nXraHv+BawVFbS/9BLtK1YQc845JN50E1EL5vv8VHchePTNN70+l/b2Usp23UVUVBOz57zLxImnoh3geCT3QGJq0tR+Awk3lV5PwpLrib/+Orq/+NJ1Pfv8czrXf0Tn+o+InDKFxJtvIu7ii5FGcT0Txr4PP/yQffv2oVarWbp0KQlHzlcdDkmlImbhQmIWLjx6fuubb7kOdf/FLzD+4Q/Hnd86FAadgQkJE9jXvo+tTVu5MO/C436moWEt+/bfj9PZg1abwJTJj5OcfC69vQ1B+5kqZuz6qKur47XXXAdynqhfXWN3I7VdtagkVb/R7InIvb20v/oqlVdcQc0tt7qKOkki5txzGff88+S/vpb4q68acVGnKDJV1f+mZOt1WCzVREZmMXv2KvLy/g9JUhMZmUFCwilBl4DhyP1aSFISGzZEU7bzAhz2GDSaFraWXkdNzf+Om/EY6gyxKjKS+MWLyX/rTcY99x9iFi4ERaFrwwaqv/1tKq++ho61ryMfOdBaCH3ufIuISOOLL3rYWnIxnZ1pqNU2Dhz4Jbv3/Aqns6ff7ww13yRJIubMMxj372coWPc28UuWIEVG0rt7N/W//BUV551H89//jqOtzWd/nxC8tmzZwjffuFrSXX311eTkjP4g4MgJE8j43e8o/PQTUm7/CeqUZBxGI81//jMHzjmXhvvux3rw4JAfz53jxy7HOhxdlJf/jN17foHT2UN8/ALmz19HcvK5rjiC+DNVFHZHmEwmVq1ahcPhoLCw8IT96twJMDlxMjERMSd8XLvRiPGppzhwzrk0/uY+rBUHUEVFkXDjjYx//z1y/v43okc5i2KztbBjxy0cOPAoiuIgNeVi5s9bR7xBHCoarPr2HFaUIubNX0dy8nkoio39Fb9jZ9n3sdlcH4bGHiPVndWoJBXFqcVDenz3gaA5//onBe++S8KypUh6Pda9e2n49a85cO55NP/lrzhaBl9+EELLp59+Snl5OXZ7HDOmv0Be3g8BiYaGV9m85So6u/Z6ftZzx/8wViR0hYVkPPgAhZ9sIOWnP0WTmoqzuYWWp//CgbPPof7ee+ndv9/bf5YQpCoqKjw9h88991ymTfPuIf7uQ92LPv6YzMceRTdlMkpvLx2vvMKhSy+j+tb/R9fnX5x0W4A7x/veGWs2l7F5yxU0Nr0BqCjI/ymzi18kUjc2DvkXS7FAc3Nzv351ixcvPuHGzqGMZi3l5bS/8AKmd9+DI+v/2sxMEm68kfhrr0EdF+eV2NvavqR898+w2ZpRqXRMKPoNmZk3iOW2IGYymXjrrbf69RxOTEgnIf5f1Na+QMWBR2lp+ZjNWy5n6pQn2NrhKvAmJkwkNiJ22M+nK8gn/b77SLn9djpefZW2l1bgaGyk5W9/o/WZZ4i77DLX+WQ+PhBUCAyTycTmzZs9x5pcfvnlFBQUAneSkHAqu8t/Rk/PQUpKrqaw8NfoEi6kylyFhMSs1FnDfj5NQgLJ3/8eSd/9DuYPPqTt+efpLSvDtOY1TGteI/q0U0m8+WaizzzTJ+eTCYFlMpk4ePAg7733nqfn8Jlnnumz55MiIjBceSVxV1yBZetW2p5/ns6PPqb7iy/o/uILIsaPJ/HGGzFceQWqAfb2zU6dDUBFewUdvR2YjWs5cPBxFMWOTpfBtKlPEh8/tlo4hn1hV1JSwrp16zz/PWPGjJP2qxvs/DrF6aRzwwban3+BnpKj1b9+zhwSb7qJ2PPORdJ453+5LNs5VPkkVVX/AhSio4uYNvUpYmImeuXxBd8oLS3lrbfe8vx3cXEx6Ufa70iSRE7OzcTHz2NX+e309ByidNu3aNJMR4Uy7P2cx1IbDCTdeiuJN99M50cf0fa/510Hgr7+OqbXXz9yIOhNxCxciOSHRtyC75WWlvL22297Zi0KCwspLj4665uYcCrz57/N7j2/pLX1E/bvfwA56nWiVAq58ZOJixj5AFTSajFcdilxl16CZdt21wfu+vV0f/U13V99TUReHgk33Uj8VVehijp+H58w9hybb4mJiVx++eV+mWiQJImouXOJmjsXW00N7S+toGPNGmwHD9L4wAM0//nPrkPdly/rd6h7kj6JAkMBTZ0H2br929BTBkBKyiImT3oErTbe57F7W1gPl0wmU7+iDlynsJtMpkF/p9XSSqWpEgmJ2WmuSt/Z1UXr//7HwUUXUvfjn7iKOo2GuMsvJ+/VV8lb8RJxFy7yWlFnsdSytXQpVVX/BBSyMpcyb+7roqgLcu6Zur42b958XL7Fxk5h3tw3yEi/FpDJcezghylW5iTmeyUOSasl7uKL+xwIekmfA0F/yMGLL6HtxZdwdnn3QFDBv0wmU78PWYCDBw8el28REUnMnPFvioruRZK0qHp28Iu0XhYmH38awEhIkkTU7GKyn3qSwvUfkvjd76KKjcV2+DBNv/3dkVMB/oi9IXgOeBWGb6B8a29vp9vLBwsPhetQ919RuPFT0n59N9rsbNexYs880+dQ9zLPz5+bnMld6b3QU4ZKFcHECb9l+rS/j8miDsJ8xq5tgA29iqLQ1tZ23Ll1bu5l2KKEIvRGM40v/g3Ta2s9p2KrDQbib7iBhGXL0KalDvgYo9FkfJe9e3+Nw9GJRhPLpEmPkJZ6sdefR/C+4eSbRhPNlCmPExE7i337fsP4SBl13cM0x8WTknKB12LSz5pF1qxZpDb8nPaVK2l/5VXs1dU0/f73ND/9NPGLF5OwfDkR2SfvciIEl7a2tuP2Fw2Wb5IkMS7nO8THz2P9putI0NhIsKzjUOV48vNuQ5K8M4Orzcoi7a5fkHzbbZjeeIO2F1/AXlVN67P/ofW//yPuwkUk3nQT+lmzvPJ8gv8MJ9/8RR0TQ+JNN5GwfDldn3xC2/Mv0LNlC+Z16zCvW0fk7FlYvpPIVPV6JLVCuxzJhfNfIzZmbG9LCesZu8TExOOmiCVJIjExcdDf2dpYwuRqhdte7uTgogtpf+FF5O5uIgrHk/7bByn89BNSf3qH14s6p9PCnr2/ZteuH+NwdGKIK2b+vHWiqBtDRpJvh+VU/tAYSZNTj9NhZmfZD9i37wGcTu/e1arNyCD1Zz+j6JMNpN9/HxH5+Z4DQQ8uWkTtT273yoGggv+MJN+c2iweqVezpVsNKFRWPknpthu93i5JHRNN4reWM/6998j++9+JOuUUcDoxv/seh29YSuWSJZjffXfE55MJ/jeSfPMXSa0m9vzzyX3xhSOHul+JI1VNzTlbaFC/j4RC3SE1T1eqkCKyAx3uqIV1YWcwGPqt/0uSxOWXXz7g6EK22eh44w3m//pVHlzhJLO0BhSF6LPOJOfZZyl4+20Srr9+wM2Zo9XVtY8tJVdTX78akMjL/T9mz16FXj/2EzCcDCff3EqaSmh1qqiMuY5xObcAUFv3IiVbr6G7+4DXY1RFRZGwdCkF76wj55l/HT0Q9MMPXQeCXnc9prffRrHZvP7cgneNJN9Km0qxKhLfOKcwZfIfUauj6OjYxOYtl9Pc8rHXY5RUKmLPPYfc//2X/DffwHDNNUhaLb07dlJ35884cMEiWv79b5wdHV5/bsG7RpJvgaCfOhXNz86h7bca7OMVpF6JhP+omfdHNU897WD/b+7CdvhwoMMcFUkJsiG42WzGYDBgMpmI89KdoydjMploa2sjMTHxuCR0tLbSvno17atW4Wx2HQ1h1YDhqivJ/M7/Q3dM6ydvUhSFuvpVVFQ8hCxbiYhIYeqUP5GYeLrPntNbAvE6jpS/Yz1Rvh3rurevY2/bXv6w8A9clHcRLa2fsnv3L7Db21Cp9EyccD8ZGYt9ujm574Gg7oJOk5pKwrJlxC+5flgHgvqKyLfBDSffHtv8GC/teYklE5dw7yn30tNTya7y2+nsLAcgO/tmigp/iUqlO+HjjIajpYX2l49cc1tbAZD0egxXXUnijTehK/DOXtPRGis5F8z55m9OZy8VB35PXd1KAOLiZjKl6A/Y129jzz/+QHyd2fWDkkTM2WeTePNNRC1YEBSnTAzndRSF3SB69+2n7YXnMb+9zvNh5kyO5+VpZirOyGX1t9716fPb7Sb27P01zc3vA5CUtJApkx8nIiLZp8/rLcHyOg5FsMZqspo48+UzUVD45PpPSNa7Xnur1Uj57p/R3v4VAGmplzFp0kNoNMM/CmU4HG1tdKxeTdvKlZ5BjqTTYbjyShJvuhFdYaFPn/9EgvU1HEgwx3r929ezp20Pj5/1OBfnu7Z5yLKVAwf/QE3NfwGIiZnCtKlPER1d4NNYZJsN8zvv0vb881j3Hj1jL3rhWSTedBPRp50W0A/cYH4d+xorcfpaV9d+dpXfTne36yzF3HHfo6DgTlQqV+eVNypeZ/XKe1m6I5ai3UdvMNJNnEjiTTcRd9mlqHS+G9CczHBex7BeinWzNzbS/c0mbPX1dH7yCVXf+Q6VV16J6bW1KDYbkdOnk/nHP/LBH67mzVNVTBm/wKfxdHSUsHnzZTQ3v48kaSksvJuZM54dM0Wd4B3bjdtRUMiLy/MUdQA6XSrFs55nfMHPkSQ1TcZ1bNp8OSbzDp/G0+9A0McfI3LKFBSr1XUg6GWXHzkQ9HMUWfa8p+yNjT6NSfCeTlsne9tcBVTfMzpd52Pey8wZz6LVJtLVtZstJVfS0PCaT/dcqiIiiL/6KvJfX8u4558n5rzzQJLo3vgZNbfcSuUVV9D+6qvIvb0AIueEASmKQl3dy2wpuZru7v1otUnMmvk/Cgt/6SnqAOamz2NXnor7r7KSte51EpYtcx3qvm8fDffcc9yh7sGcb2F9VyxAx5o1NNx3P8hy/2+oVMQuWkTiza47tCRJYsu6F4GTt9kZKUVxcrjqn1RWPoWiONHrxzFt6lPExc3wyfMJwc19XuJA+SZJKvLy/o+EhAXsKr+D3t4atm69nvEFP2PcuFuRJN+N2aSICAxXXEHc5ZcfORD0BTo/PnogqDo52bWEpiigUpHx2weJX7zYZ/EI3rHNuA0FhXGx40iNOv7mr+Tkc1gwfx3l5XfS3vENu/fcRVvbl0yc+Fs0mhN34BkNSZKIXjCf6AXzsVVV0fbSCkyvvYa14gCNv7mP5j89QeTMmXR//rnrOi5yTjjCbjezd989GI2uFbbExDOZMuWP6AaYJMmKySItKo2mniZ2R3dwyn2/IeX2n9CxZo3rUPeGBs+h7pFTp2LZsSNor3FhPWNnb2yk4Tf3HVfUxd+whML1H5L95J+JKi5GkiS67d3sadsDDN6IfTSs1ia2bbuJQ4eeQFGcpKVdwfx5b4miLowNpcOJwTCb+fPWkZp6CYri4MDBx9i+47tYbb5vFeY+EDT7L08z/sMPSLz5ZiS9HmdLi+uCByDLNNx3f1COaoX+htJGTKdLo7j4BQoK7kSS1DQ2vcnmLVdgNpcN+jveFJGbS/o9v6Zw46ek/vKXaDMzcXZ00L1x49HruMg5ATCZtrF5y+UYje8iSRoKx9/FrJnPDVjUget65s5997VXbTCQdMstFK7/kKw/P4F+1iwUux3L9u1BfY0L68LOdrjq6IvTR9zFl6DN6n9u1zbjNmRFJismi/Ro7/aLa2n5hE2bL6O94xtUKj2TJz/G1ClP+HzPlBC8uu3d7G7dDZx8IKHVxjFt6tNMmvgQKpWOtrbP2bz5UlrbvvBHqABEZGeTdvevyPrzE8d/U5axVVX7LRZhZLY2nnwgASBJavLzfsjs4pVE6jKxWKoo2Xod1dX/QVHkE/6ut6hjY0n6zrcZ/+EHJP/oR8f/gMi5sKUoMocP/5OtpUvo7a0lMjKHObNXk5v7/ZOuZLhzv6SxpN/XJY3Gdaj7y6tI+81vjv/FIMu3sC7sIvJy4dhehSoVEbnjjvtZdwXvzdk6Wbayv+L37Nh5K3Z7GzExU5g/7y0yfXyXoxD8dhh34FScZMVkkRGTcdKflySJrKylzJv7BtHRRdhsLWzf/m0OHPwDsuy/s8AiJ00a8ntKCB499p4hDyTc4uPnMn/+OlJSLkRR7FQceJgdO/8fNlurL0PtR9JoiF98rcg5AQCrtZnt27/NwUN/QFGcpKZeyoL5b2MwzBrS77tzf2fzTmzOgY90ij3v3KDPt7Au7LTp6WT89sGjL9KRtXJt+vEzcu4K3lv763p6KinZeh01Nc8BrmME5s5Z4/M7zYSx4UT7604kJmYC8+a+QVbWMkChquqfbC1disVS44Mojzec95QQPLY3b8ehOMiIziAzZuitxLRaA9On/Y2JE3+HShVBa+unbNp8GW1tX/kw2mNiEDknAK2tG9m0+RLa2r9EpYpk8qRHmDb1qWGtfOXF5ZEUmYRNtlHWMvD2grGQb2F/80T84sVEn3EGtqpqInLHDfjiWBwWdrXuAk68/2SoGhrfYN+++3A6u9Fo4pky+TFSUs4f9eMKoWM0M8RqdSSTJv6OxITT2bP3bszmbWzafBmTJz1MWtql3g71OEN5TwnBZTT5JkkS2VnLiDfMoWzXT+jpOcC27TeRl/sD8vPvQKXy/ceMyLnwJcs2Dh76E9XVzwIQEzPpyHE8wz9+SZIk5qTN4cOqD9natHXQgXWw51tYz9i5adPTiV4wf9AXZ2fzThyyg9SoVLJjRt7tweHopnz3z9m9+2c4nd3Ex89nwfx1oqgT+ul19HpGi6OZIU5NvYj589ZhMMzG6exiV/lP2LP31zidFm+FOqiTvaeE4OKNFYmYmInMn/cGmZlLAIXDVf+gdNtSLJY6L0V5YiLnwk9PTxVbty7xFHXZWTcyd87aERV1boPtsztWMOebKOyGoO9odqR738ydu9i85QoaG18HVOTn38Hs4peIjDz5/ikhvJS1lGGX7aTqU8mJzRnVY+n1WcwuXkVe7m2ARH39araUXE1X1z7vBCuMeVan1TOQGO2KhFqtZ/Kkh5k29WnU6hhMplI2b7kMo/F9b4QqCB6NjW+57sju3IlGY2D69L8zceIDqNWjO0TY/R7Y3rwdux/3J3uTKOyGYKT7ncB1OGJ1zX8pKbkOi+UwOl06s4tXUJD/YyRJ7e1QhRDQd/bEGzfRqFQaxo//GcWzXiAiIpXu7gq2lFxNbe0Knx4wK4wNO5t3YpftJOuTGRfrnQ3gaWmXsmD+OuLiZuFwmCnb9UP27r0Xp7PXK48vhC+ns4fde35J+e6f4nR2YTDMZcH8daSmXOiVxy+ML8SgM2BxWNjTuscrj+lvorA7CZvTxs7mncDwR7M2Wxs7d36PioqHUBQbycnns2D+OhIS5vsiVCFEeGaIvbCfs6/ExNNYMH8dSUkLkWUr+/bfR9muH2K3m07+y0LI8saKxED0+hzmzH6Z3NwfAFBXv+rIbPF+rz2HEF46O3ezecuVNDSsASTy837M7OIVREYO/Yafk1FJKmanzgaOvjfGGlHYncSull1YnVYSIxPJjxt68+n29m/YvPkyWlo3oFJFMGHCA8yY/k+02sA3TBeCl91pZ0ezqzWYLzqcREQkMXPGsxQV/hpJ0tLc/AGbN19GR8eJ95MIoWs0KxIno1JpKRz/C2bNep6IiGS6u/ezpeRq6upWidliYcgURaGm9gW2lFxLT88hdBFpzC5+iYIC39yc49ln1zQ2r4uisDuJvqf/D2U0K8sODh76M6XbvoXV1kRUVAFz56wlJ/tGcTadcFLlreX0OntJ0CVQYPDN0TeSpGLcuFuYO+cV9Ppx9FrrKd22jMrKv6IoTp88pxCc7E47O4yugYQvOuq4JSWewfz575CYeCay3Mveffeyq/wn2O1mnz2nEBrs9nZ2lv2A/fsfdK18JZ3L/PnrSEg4xWfP6V4t2da0Dac89q6JorA7ieGMZnt76yndtpzDh/8KKGRkXMf8eW8SGzvZx1EKoaJvvvl6IBAXN4P5894iPe1KFMXJoco/s23bTVitTT59XiF4uAcS8bp4CuJ9e4amLiKZWTOfo3D8L5EkDUbju2zechkmU6lPn1cYu9rbN7Np82W0tHyEJEUwoeg3zJjxDBERiT593okJE4nWRtNp72R/+9jbOiAKuxNwyA62GbcBJx/NGps/YNPmyzCZSlCrY5g65c9MmfwoanWUP0IVQsRQ+nV6k0YTy5Qpf2LK5MdRq6No7/jmyIX0E788vxBYfVckVCdpt+QNkqQiN/d7zJnzCvrIcfT21rG19AYOH/6H39qRCcHPNdB8mtJty7FaG9Hr85g3dw05Od/2y8qXRqWhOLUYGJv77ML+gOIT2dO6B4vDQlxEHEUJRQP+jNPZS8WBR6irewmAuNgZTJ36JFFRuSN6zoaGhhHHG0gZGeLYltFyyA62NbkGEr7Y7zQYSZLIyLiWuLhidpXfTlfXbnbsvJWcnO9QOP4XqFSjOz5ACF6+3F93Ioa4mcyf/xZ7995Lk3EdBw/9kbb2r5g65U/odKl+jUUILr29DZTvvpOOjs0AZKRfw4QJD6DRRPs1jjlpc/ii7gtKmkr41pRv+fW5R0vM2J2Au1KfnTZ7wNFsd/cBSkqu8RR148b9P+bMWT3iok4Ib/va9tHj6CE2Ipai+IEHEr4UHV3AvLlryM6+GYCamv9SsvU6enoq/R6L4HtO2TnkFQlf0GhimTr1SSZPehSVSk97+1ds2nwpra0b/R6LEByamz9i0+bL6OjYjFod7VpNmPIHvxd1cPQ9sbVp65i70UcUdifgWRY75qKnKAp19avZvOVKurr3odUmMWvmfykq/BUqVUQgQj2pHrOd5spuesxj88DFUNLV3kvtvna62vuf6eXOt9mps1GrAnPGoUqlY+KE+5gx4xm02gQ6O8uPHC/wekDiEUZvsHzb276Xbns3sdpYJiRMCEhskiSRmXkd8+e9QUzMJOz2Nrbv+C4VFQ8jywM3YRdCj9NpZd/+B9lZ9n0cjg5iY6cxf96bZKRfFbCYpiZNJVIdSYe1g4MdB4/7/mDvq2AglmIH4ZSdlDa5NvX2Lewcjk727L0Ho/EdABITTmfKlD+h06UEJM6hOFzaQem6BlAACWZflkHe7PhAhxWWdn9Zzycv7QUFJAnO/tYkppzuOoNpsIFEIKQkn0fs/HWUl99JR8cmdu/5OW3tXzBxwoNoNDGBDk8YohPl29ZG14pEcVpxwAYSbtHRhcyds5YDBx6htu5Fqmv+Q3vHZqZNfZKoqLyAxib4Vnf3IXaV/4SuLtdhwONybmH8+J8HfJJEq9YyM3Ummxo2sbVpK4UJR9uUneh9FQzEjN0gKjoq6LR3Eq2NZmLiRABMpu1s2nw5RuM7SJKG8ePvYtas/wV1Uddjth8t6gAU2LauQczcBUBXe6/nYgCgKPDpir10tfciK7JnIOHv/U6DidSlM7v4RfLz7wBUNDa+weYtV2Lu3BXo0IQh6Grv5dNB8g0Ct79uMGq1jokTH2DG9H+g0Rjo7Cxj85YraWx8K9ChCT6gKAr19WvYvOUKurr2oNUmus7YLPp1wIs6t4HOs+tss/DJi4O/r4KBmLEbhHt/3azUWaglFVVV/+LgoSdQFAeRkdlMm/okBkNxgKM8MUVW2PNJ89Gizv11BbrbbETFaQMTWJjqMFqOfy1k6GjqoU6pwmwzo9fomZwUPMfjSJKagvwfk5BwCuXld2CxHKakZDGFhb8kJ9s/d6gJI9NhtHDs1iBFhpbaLqLiIyg1Hr8i4Q8nv0FsOvl5z1FX91t6LDso3/1TauvWk5F+BypV4E4ZEDeIeY/D0cnefffR1OQq2hMSTj1y40xagCPr79h9djaLg4/+e3ybMUUGk9FCTEKkv0MckJixG4S7X+f8lMls3/FdDhx8HEVxkJp6CfPnvR30RZ2l084XL1ZTtf34dlGSBNGJwTEiCifxqXoGqoO2vHOYzYdcH7LFqcVofHCS+mglxM9jwfx3SEm+AEWxU1HxEDt2/j9sttZAhyYMYrB827hqP1t37cZkNQXdQMJNq00jN/cpkpO/A6gwmd7lUOWt9PZWBDo0YZTM5p1s3nIFTU1vuQaOBXdSPOv5oCvqAKYnT0er0tJsaWbHrv2s/v0W6is6jvs5SQWGVL3/AxyEKOwGoCgKW5u2MlHnJM/0Am1tn6NSRTJp4u+ZNvVptNq4QId4Qg37O/n4n5U0H+5BrZXInWXwXOAlCYovyxCzdQEQkxDJ2d+ahOcGawkktUR9RQdtLxnI7pgYFPvrBqPVxjN9+j+YMOEBVKoIWls/YdPmy2hr/zrQoQkDODbfJAl00Rq62nrZ/I8mZtSfzazkWWhVwXktkCQNqSm3kDvuKTSaFGy2aioPf4+2tjVj7i5FARRFpqr6WUq2Xo/FUk2kLpPZs1eRn/dDJCmwezwHE6mJZHrSdGbVnc+Xf6+ls7WXuORI5lyUe/R9pYKzl08Kmtk6EEuxAzrYsZ8zIps4L86B4rASHT2BadOeJiba/0dQDIfTIVP+cTMHvmkDwJCuY/61WcQm65h8TgrdbTaiEyNEURdAU07PZNyURExGC4ZUPVaLgw//XU5bA1y25zaSk3Q4p8qo1cE55pIkiZzsG4k3zGVX+e309Bxk27Ybycu7jfy8n/ikb6MwcsfmmyZCzScv7eXQtmZOq7oaSemm51QbUXHBO4MfHV1MQf5/qW94hK6uL2lsepKu7hIyM+5GozEEOjxhCGy2Fnbv/gWtbZ8BkJJyEZMnPYxWG9yvX7fJyryt1yLVuW4YK5qbysLlk9DpNUxbmOV5XwVTUQdixu44Fks1B3Z9j/PiHABkZS1j3tzXg76o62q1sfG5Kk9RN35+Amffkkdssutw2ag4LSl50aKoCwIxCZFkTUwgJiGSpMwY5v4wifK0LwCo/8LK638sxdxiCXCUJxYbO5n5894gM+N6QOHw4b9Rum0Zvb31gQ5NOEbffIuM1nLh/5vK1gnv4JDsKNXRrH5oMzV72wId5glpNPHkZD9KWtrtSJKWrq4vOFT5bbp7tgc6NOEk2tq+ZNPmy2ht+8x1nNLE3zF92l+DvqirKm9l9UObkepisKusbJv8LhfcMhWd3jV47fu+CjaisMN10nVb+9fU1r7Eps2Xo7XX0iNDfcwVTJr4O9Tq4Hvh+qreaWLDM5V0NPQSoVdz6g3ZzLw4HbXm6MvbZHeypdtGk33sNTQOddtaS/m84FUOzt+ILkpDU6WZ1Q9tpqIkuHu2qtVRTJ78CFOnPolaHYPJtJVNmy/D2PwBlaZa3qveRKWpNtBhCseo6qxiS9KHvD3zaeLT9fSYbbz11Ha+eeMgTmfwtvWSJImkxOvIy/sXERE5OBzNVFX9hObm51AUJ7WWJj5v3UWtJbjfN6HO/Xna01PDgYN/ZNv2m7HZmomOLmLe3NfJzloW1DddOR0yX752gHV/2YGl005CZhRvzPgzm+I/oKH76I0/9b02vmjvpL43+M5bDPt1k/r6V9iz9x7g6AWtxq7juWaJP01bFrjAhsBhk9n+biPVO1w3SCTnRjHvmkz0x8zKvd5h4aFGMzKuSv7e9Diujg+ejZ7hzn0HduHsNK6/ch7r/7ObxkMmPny2nNo9bZyxZALaiODcgwKQnnY5hriZ7Cq/A7N5B38pe41nSUGRdEiKkd+k7+K2KRcFOkzhCHe+5eSlcv135/PFKxXs/qKere9XUbe/nQu+O5W45OC9PugjJ1CQ/x8aGp/EZHqX5pbnWNth5l+OJShSCpIi8zPDNyzPPCXQoYadgT5PATIzb2BC0b2o1cGbVwCm5h4+fLYcY1UnANMXZnHa4kLe+TCF1pYGtjZtJTMmk5X1rfx8X43nM/WPE3NYlpkU0Nj7CuvCrre3gT17f03/MygknmuW6FJ0TE+eHqjQTqq9wcKWNfV0tdlAgskLk5l0ZjKSqv9IqMnu9BR14Hq7PdRo5rToCNK0wVsshAtFUfqdJxaXpOfqnxWzeV0lW9+vYveXDTQcNHHh/5tGUpb3Dwb2Xm9iLVmZT9LGv3nWfDXKkZ3FiqTiocYU5kVuJ1vvvbvexNETI9c337QRas751iSyJyXw6Ut7aTxkZvXvt3DOtyZROCd4e7aqVFFkZf6amOi5lDU8e6SoO5pzT5jGsTChyas5J5yY6/P0+KJu4sTfkp21PDBBDcP+LY18umIf9l4nuigN5940mYJZrjNq56TNYWfLTkqaSpiTdaGnqAPXX/uLfTWcnRhLZmRw7FUN66XYHsthjjtYDIUkjcL05OlEaoJvCVZRFA5808bG/1TR1WZDH6fhrJtzmbww5biiDqDa5uTYxRUZqLGJJdlgUNtVi7HHiEalYUbKDABUahWnXDmeK2+fRZQhgvbGHl59pIRdG2uD+m5ASdJiijzL8wHrJktqKnuaAxSV0JeiKJ6jnOamH70Du2huGkvunU9afhw2i4MP/r2LT1fsxRHk1wmDYREdCT8VORcEXJ+nxy/lR0eN93ssw2G3Ovn4hT2s/89u7L1OMgoNLLl3vqeog6Pvla1NWznYYz3ur3QClRar/4I+ibAu7KL0eRz7v0BBotkhBc1p7H1Zexx8/XItOz9oQnYqZEyK4bwfFJCcO/ihnfUD7KlTATlBvLQXTtwfstOSpqHX9F+myJ6UyA33zid3WhJOh8zGVft5/1+76O0O3q4hem0Kx56Kq1Kc5EcFb3eWcFLXVUdTTxMaScOM5Bn9vheXrOfqn89m9oW5IEH55/W8+mgJrfVdAYr25GyywsqeguO+LnLO/wb6PAUVen1uAKIZmuaaTl55eAt7v2pAkmDupXlc9dNiYhP7T+rMSp2FhESVuYpXB1jlUAP5ep2foj65sC7sIiMzmDzp9xz936Dig+4UTE5V0J0n1ny4m4//WUnj/i5UaomZF6dxyvXZROgHL9D29tp5rMm1V8A9l+feYyeWYYODe79T39mTvvSxEVx62wxOX1yISi1xaHszqx/aTP2BDj9GOTSdTplHWtSuA9OOFHcqxcmdhmqxJBYk3Pk2NXkqUdrjB4RqtYpTrx7PFT+ehT4ugrb6btds8Wd1QTdbrCgKv200s8cKETiRFNc8isi5wBjo83TypN8TGRl82yYURWHnJ7WseayEjqYeouN1XPnTYhZcXoBqgKOm4iLimJQ4CUv02bxi7AGO/pVq4A8Tc4JmGRbCfI8dQGbm9SQmnonFUkWnoueOt7+FWlIzM3VmoEMDQJYV9m5sYe9nLQDEJkcw79os4tNPvExstDv5SW0HFkVhflQEv0mPpdEukxOhFkVdEBlKv05JJTHr/HFkFsXz4bPlmJotvPGnUuZfns/si/JQDbAE7292ReEXdSYO2ZykaFQ8nibTaWsmPyqFbL3YxB4shtofNmeKa7b44//tpnp3GxtX7qN2bxvnfGsSuqjgODLpmdZu3jH3ogaeykkiW9VBZY/IuUDq+3mq1+cGZVHX22Vnw4t7qNzh+kzNm5HMuTdNQh9z4sIsOfkCuiyzAfhZXhrLM5KotFjJ1+uCqqgDUdgBrpFGZGQGXxx8G4ApSVOI1kYHOCroMdnZsraO1mrXmWa5xQZmXpSOJuLEE609sszttR00O2TyI9T8MctArFpFdnDlXthr7G6krqsOlaSiOPXkLepSc+O4/p55bFy1j/2bmtj0ViW1e9s5/ztTiUkI3DKAoig82tjJph4bkRI8lR3P5EgtIGZMgo1nhngIKxJRcRFc9qOZbP+ohm/eOMjB0maMhztZdOtU0gsCewbZuyYL/2zpBuDu9FhOidYBaWKWLgi4P0+DUX1FO+uf201XuxWVRuK0awqZcU72SY9f2dfdy/u2OSCpSLbv5Od5NyJJUtAVdG5hvRR7LPdFLxj219XtMfPxPw/RWm1BE6Fi3rWZzLki86RFnVNRuLvexF6rgwS1xF+yE4gN0i4G4c49ezI5cfKQBxIRkRou+M5Uzvv2ZDQ6NXX7O1j90GYOl7X4MtQTerGth7UmCxLwSKbhSFEnBJum7iZqOmuGPJAA12xx8aJxXPOLOcQlR9LZ1svaP5ay9f3DKHJglma39dh4oNEMwE2JUVwbP/geY0EA18rX5nWVvPHENrrarcSnRbH4rrnMPDfnpEVds83Ot3YeokdWoendBw1P0mHt8E/gIyQ+8fsYzmjWV5x2me3vNLLplTrsvTIJWZGc94N8cqYNbYT8Z2MXn3XZiJDgz9nxZImbJILWaPJt0ikZLPn1PJJzYujttvPO33byxSsVOO3+PWB2Q2cvTza7Ntf/LDWGs2OD705ywcWdb5MSJxETMbyjc9Ly41hyz3yK5qWhyArfvHGIt57eTrfJv3cCVtsc3FnXgV2Bc2J03J7i/SOAhNDS1d7Lm3/expZ1lSgKTDo1nevunkvKuNiT/q7FKfPtskpqem3kRkYwy/4mEnZKm0r9EPnIiaXYI5p7mjlsPoyERHHa0Eaz3tbW0M0nzx7GbHRdLCeclsSUc1NQqYe2h+qV9h5WtLs2dv42w8BMfXBOEwsu7jtiRzpD7B51fvX6AXZuqGXHhhrqD3Sw6JapxKf5fhaj3GLnnnoTCnBdvJ5lCWLmJJgNdX/dYCL0Gi747hRyJifw2cv7qd3bzuqHNnPet6eQO9X3h7OanDI/qe2gw6kwOVLD7zMNqIK4g0Gw894Zlv41nDMsK3c08/ELe7B2O9Dq1CxcNpGJC9KH9LuyonD73mq2mnswaNS8NKOA1WVTqGovo6SphPNyzxvpn+BzYsbuiK1G12h2YuJE4iLi/PrciqKw+4t6Xn14C2ajFV20mtOX5zDtgtQhF3Vfdlk9d8D+MDmaC+PEzEkwa7G0eAYSs9Nmj/hx1FoVZ14/gUtvm0FktJbmatft+/u+8e1Fu8Hu5I7aDnoVOC06grvSYoO6TZDgnRUJSZKYfFom1/96HklZMVg67az7yw6+XFOB0+G72WK7ovDzug6qbE7SNSqeyo5HHwQ3DQnByWF38tnq/bz7jzKs3Q5SxsVy/T3zhlzUATxe2chbxg40EvxnWh5F0ZGe9477vRSsRGF3xGhnT0bKanHw4X/K+eSlvTjsMqkF0Zz3gwLSCoe+xFDRa+eX9SZk4PK4SG5JCvyNH8KJuS8MRQlFGHSj34ieNyOZJffOJ7MoHrvVyUf/28NH/92Nrdcx6sc+VpfTdXNOi1OmSKfhsUwDGlHUBbVWSyuHTIcAmJ068oGEW0J6NIt/NYfpZ2cDsP2jGtb+YSsdR46C8CZFUXio0UxJj51olcRT2fGkaMQWE2Fg7Y3drHlsK2WfuPpUzzo/h2vvmkN86tBXFF5uaOXJKlfP4T9OzOGMBNeyrbs+2Nu2F7PN7OXIvUcsxR4RiP11jZUm1v+nHHNLLyqVxIIrC0ifrh3WzEezw3WsSbesMEev5TcZcWLmZAzwRb7FJLjOYtr63mG2rKtk36ZGGitNXHjrtCHtJxkKh6Lwy3oTFVYHSWrXzEmMuDkn6JUaXXuCihKKiI+M98pjarRqzrphAtmTEtjw4h6MVa7Z4rOXTWTC/KHPjJzMc609vGXqRQU8mmlggrg5RxiAoijs/bqBz17ej8Mmo4/Vct7NU4hIsmFsbhry45R02/h5TTsAtyRFsxBbv2XrrKgs6nrq2LB3AwtSFnj973AbTdtEcUUG2nvbOdBxAGBUy2JDpcgKpR9U8fofSjG39BKbFMnVv3Cd+D6coswiK9xR20GjQyY3Qs2fsuPRiqJuTBjtfqfBqFQS8y7N56qfzSYmQYfJaGHNYyXs+Lhm1AfMKorCH5o6+ar76LEmGeJMxDHBsyKR6v0ViYJZKSy5Zz4ZhQbsvU7WP7ebj5/3zmzxB+Ze/triujnnrrRYzogJntP9heBhszhY/9xuNrywF4dNJntSAkuOdO0ZjsNWBz+r68ABXBCr47bk41e/pie4esjvbN/pjdB9QhR2HB3NjjeMJzEy0afP1W2y8vZftvP16weRZYXCuaksuXc+6fnDW46TFYV7603s7nUQr5Z4Ojseg5g5GRNMVhMV7RWA75b+MwvjWXLvfPJnJiM7Fb54tYJ3/r4TS6dtxI+5qt3CKx2uY00eyjQwVS9mTsaKk3U4Ga3YxEiu+mkx8y7NQ5Jg79eNvPpICc01nSN+zB0WG/c1mABYlhDFEnFzjjCApsNmVj+8hYotTUgqiVOuKuDyn8wi2jC8QUC7w3VzjllWmB6p5bcZA9+cMyPB1YqvrL3MK/H7wrAqgUceeYR58+YRGxtLamoqV111Ffv27ev3M4qi8MADD5CZmYler+fss8+mvLzcq0F7m7/211XvbmX1Q5up2dOORqvinBsnseiWqej0w18Rf7q5iw1dVrQSPJEVz7iI0FtVD9V8c3/I5hvySdL77m7CyGgtF/9gOmfdMAG1RkVVWSsvP7SZ2n3tw36sjZ1W/mh0fUjfnhLDeSF4rEmo5pvJamJ/+37At9c4lVrF/MsLuPKnxUTH6+ho6mHNYyXs/GT4s8V1Nic/re3ApsBZMRHcmRqax5qEas75gyIrbPuwmrWPb8XcbCE2MZKrfzabOSPoxmOTFX5W10GN3UmmVsWT2fFEDvIY7hm7/eb9WByWUf8dvjCswm7jxo388Ic/5JtvvmH9+vU4HA4WLVpEd3e352cef/xxnnjiCf7617+yZcsW0tPTueCCC+jsHPnIzde+afgGgMKEQp88vtMh89VrB3j76R1YOu0kZUVz3d3zmHJ65oj2w63t6OH5Ntcm5QfS4yiOCs1jTUI13z6r/QyAKYlTfP5ckiQx/exsFv9qLgnpUfSYbLz55DY2vXUI2Tm0uxj39tq5+8ixJtcY9NyUGJozJ6GabxuqN6CgkB2TTbI+2efPlzUhgRvunU/ejGRkh8Lnqyt49x9l9HbZh/T7nU6ZH9e20+5UmKjT8EimAXWIbjEJ1ZzztR6zjXV/28FXaw8gywrjZ6ew5N55ZIwf/o1oiqLwQKOZbRY7MSqJp7MTSNQMXhql6dNIi0xDVmTern2b5t7m0fwpPiEpo9h409zcTGpqKhs3buSss85CURQyMzO54447+OUvfwmA1WolLS2Nxx57jO9///snfUyz2YzBYMBkMhEX5/tjR1buWckjmx8BQIWK+0+7n2uKrvHa45uaLXz4n3KMh1130ExbmMXp1xaiGeTg4JOdLfRNt5Uf1XTgBL6fFM0PguSAzmM3evridfRFvvkq1sGsrVjL/V/dD4CExAOnPeDVfDsRu9XJ56/sZ8+XrhzLGG9g5uUpRBkGX1I12p18q6qNZofMgqgI/pITHPs4Rb4NzdqKtTzw1QMouC7zD572oN/yTVEUyj6t5cvXDiA7FKLjdVzw3SmoYnsH/R27ovDjmg429dhI0ah4MTcxaHpbj9WcGyzOsXqOnaNDx0f/3U2P2XbkuKcippwx+CTJyf7Of7V08c+WbtTAX3Pij7SnO7HbN93OHtMewFU33D7ldi7OvnjYf8uJjCbfRrUpy2Ry7X9ITHTtS6usrKSxsZFFixZ5fkan07Fw4UK++uqrAR/DarViNpv7/fOXxu5GHt38qOe/ZWQe/PpBGrsbvfL4FVuaeOX3mzEeNqOL0nDx96ezcOnEQYu6kzlodfCLOhNO4JK4SL4/wMbOUOaNfIPA5VxjdyMPfvWg578VFK/m28lodWrOvXEyi26dSkSkmoaDJj7+5yHq9gz89/ftOVwQoeYPWYagKOr8JSTy7esHPUUd4Nd8kySJGefksPiXc4lPi6K7w8qbf97G7k+bkQdoR9a357Becu0bDpaizl/G+meqL8lOhV0fGXnr6e30mG0kZkZz3a/mMvXMrBGfBNG35/CvPT2HT6y5t5m9pr1H40Lmqd1PBdXM3YgLO0VRuPPOOznjjDOYNm0aAI2NrgtGWlr/RsxpaWme7x3rkUcewWAweP7l5OSMNKRhqzZX97voAciKTE1nzage1251suGFPXz4n3JsvU4yxhtYcu98CopTRvyYrQ4nP6ltp0tWKNZruT89vI418Va+QeByrtpcjUz/5U9v5NtwFc1N4/p75pOaF4e9V2bTK3Vse6ehXzuyY3sOPx1mPYdDJt+UwOdbSk4s1909l0mnpqMosHdjC58/X0WPqf/S7AvH9ByeFGbHmoTCZ6qvdHfY+Ox/Vez/shUUmHpmJot/NZekrJGvWJX26Tl8c2IU1wyx53BdT93xdQMy9T31I47F20a84/5HP/oRO3fu5Isvvjjue8cWHIqiDFqE3H333dx5552e/zabzX5LxHFx41BJqn4XP5WkIid25M/fUtvFh8/uor2xBySYe3Ee8y7NQzWKD8VeWeGntSbq7TI5WjV/yoonIsxOXfdWvsHQc87bSxWRvZGoUPUr7lSo0Fl0Xn2uoZx/ZEjRc80vZrNhRRn7v2qlsqSD1moL8xdnEZei44kw7zkciHzzNl9c30YqIlLDeTdPIWdyIp+8tJfWagsb/lXJ7CsyyJwUy8edvTzVp+fwwtjwO9YkFD5TfaG23My2txuwW2W0OhXn3jSFwjmpo3rMapvrWBO7AufG6PjJMLY0ZUVlDXgdz4zKHFVM3jSiauPHP/4xb731Fp988gnZ2dmer6enuw6lPHYkYTQajxtxuOl0OuLi4vr985f06HTuP/V+VJLrf4NKUnH/qfeTHj38wzXd+0nWPFpCe2MP0YYIrryjmAVXFIyqqJMVhfsaTJT12olTuZYnEk6wsTMUeTPfIHA5lxKZwu1Tbkd15G3n3puREjnymdzRUKtVTLsgldO/lYMuWo3ZaOWTZyr5+84WVh7pOfy7MOw5HCr55s3rm7dMmJ/Oud/PJz4zEpvFyTera1n7caOn5/D1YdpzOFQ+U73JYZcpfbuBzWvqsFtlErP1nPuD/FEXdX17Dk+J1PDQMHsOB9t1fCDDmrFTFIUf//jHvP7663z66afk5+f3+35+fj7p6emsX7+e4uJiAGw2Gxs3buSxxx7zXtRedE3RNZyWeRo1nTXkxOaM6KLX221nwwt7qNzRAkDe9CTOvXky+pjRfyD+vaWb9Z1WNMCfsuLJ04XesSaDCcV8uzj7YuYmz6W+p57MqMyguBikjY/hvB8UUPJGPV9123hZ4wAk/i8+ikVh1HM4FPPNG9c3b4tJjODs7+ZR/rGRkh0dPJeoYFUkFmg1/CLMeg6HYs55g8nYy+Y1dXQ2u87dnHhGEpPPThly7/TBHNtz+MkR9hwOxut4X8OqEn74wx+ycuVK3nzzTWJjYz2jCIPBgF6vR5Ik7rjjDh5++GGKioooKiri4YcfJioqimXLlvnkD/CG9Oj0EV/w6g90sP4/5XS1W1GpJU67ppAZ52Z75eL0ZoeF/7S6Nnb+Jj2OudHhNXMSqvmWEpkSdBeCyBgNKddm8PrhNhQJZlZaydnfSdu1GhKz9YEOzy9CNd9Gc33zFZVaIv+8FO4bJ9GFQmqHg3M+a6f2AjW5swxhU9yFas6NlKIoVG7tYOcHTcgOBV2MmnlXZ5FaMPobBRVF4XcNR3sOPz3KnsPBeB13G1Zh949//AOAs88+u9/X//vf//Ltb38bgLvuuguLxcJtt91Ge3s7CxYs4MMPPyQ21ju9KoOFLCuenpyKAoZUvVd7cm7ptvHQkY2dtyRFc0V8eHy49iXyzX+aHU7uqOugV4KZajWLK2xYOuxs/O9hppyTwoTTk0L+w1bkm/+4ew4fRiFJJXHbYQdOi0LpWw0YD3VTfFk6Wl3o7+sUOXeUzeKk9O0G6ve4zudLK4xmzlWZREZ7Z5XqP63dvG0+2nO4KIRvzhnVOXa+4O9z7Eaiq93KR/8tp25/BwATF6Rz1tIJRESOPgEbGho4bHVwU1UbnbLColgdjwxzD0Ag+OOMJ18JtXOehtM8uqGhAYuscEt1G3t6HeRGqHk+N5Eou8K2dY3UlrsGF6kF0cy9OpPImODYChCK+RYu6uvreaSpk1c7LERK8Oy4RKZEatj/RSu7P2lGUSA6Qcu8a7NIzAqeAe1YzbmxcH1rrelh82t1WEwOJBVMOy+VwlMTBxxMDuf6Bq6/8wNzL7+qdx0lc3daLNePgX2co8m34LhKjyGHd7bw8fN76O22o9GpOXvpBCaeMrxEOxF3v7rOI/3qHhykX50geIO75/CeY3sOq2HetZmkFkSz471GjIe6+fifh5h7VSZphcFxKLYwNq1s7+HVIz2Hf9+n5/DEM5NJzotiy2v1dLfb2fjcYaael0rRIB/wwtinyAr7vmxlT5+Cfv7iLBIyvVfQ9+05vDwhakwUdaMlCrshctplvnr9ADs31AKQMi6WRbdMJT7Ne0nS65S5c4j96gTBG546Qc9hSZLImx1PYo6eLa/VYWqy8uWKGopOS2Tquamj3sgshJ8PWkz8yeg61uT2lBjOPabncFJOFOf+IJ9tbzdQt7uTXeuNGA91M/eq4JktFrzD0mmn5PV6mitdd+DnTI9j1qXeXYKvslj79Rz+aYj2HD6WeKcMQUdTDx88u4uWGtcFaea5OZx69XjUWu8dO6IoCnfuq2H7EPvVCcJovVjfwgtD6Dkcl6Lj7FvzKPvQyKEt7VR81UbL4R7mL84iOiG8bujxpWBaGhuOoS6N7ezs4QflVSftORwRqWb+4iwOl3aw4/0mjAePzBZfnUna+PD4YA51jRVdlLxRj63HiVorMeuSdMbN9O5NMya7g2/tPES7U2FSiPccPpYo7E5i7zcNbFy1H4fVSWSMlvNunkzedO830v7j4UbWNrWjAf6QZWB8GB1rIvjfxrZOfrXfNfv8g+RoLjGceOlDrVEx65J0UvKjKH2rgfb6Xj7+ZyXFl6eTM234jbeF8FLfa+OmnZVYZJlToiL4VfqJjzWRJIn8OQkk5USxeU0d5mYrX75Uw4TTk5hyzuiPvRACw90W7MA3bQAY0nXMvzaL2GTvHkhtlxVuLT9MRY+VVI2Kp7LjiVKFz0SJqB4GYet1sHHVPvZvagIga2I8F3xnKtHx3j8RfU1jG3867Hqeu4fYr04QRmpvt4Vbd1XiVODSuEi+lzT0owSyJseRkKlny9o6WqstbHmtHuPBbmZenI4mInwunMLQdTuc3FRWSaPNzoSoSB7Pih1yz+G4VB3n/L88dn7QROXWDvZ/2UpLVQ/zrs0kOl7MFo8lXW02Nq+po6OhF4Dx8xOYdkEqai+vTCmKwq/21/B5exdRaldRlxpmPYdFYTcAY5WZD58tx9RsQVJJzL8sn9kX5aLywX63bzq6uHOvq3fjbTmpXBMtRqKC7zTb7Hxr5yE6nTILDNHclxY97OWPKIOWM2/OZe/GFvZ+1kLVdhNtta52ZIa08DnQWDg5p6Lwg91V7OqykKTV8OKMfCI62ob1GGqtiuLLMkgtiKb07Qbaai1s+GclxZdnkD01eO9GFY6q3mli+zuNOGwyEXo1s6/MIHOib45r+Vu1kRUNbaiAf07JZZK9xyfPE8xEYdeHIivs2FDD168fRHYqxCTqWPTdqWQUxvvk+Sp7rHynrBKbonBpioF7x2fQdIJm4oIwGhanzM1lldT22snTR/DctHxsrc0jeiyVSmLKOSmk5EexZW09nS02Pvn3YaYvSqVgXoK4i1EA4IEDdaxvNROpknhhej65eh0NHSN7rKwprtnizWvraKuxsHlNHcZD3cy4KA2NF/c7C97jsMlsf7eR6h2uu1KTc6OYe00mUXG+OUPuneYOHjrk2qv6YGEWi5INNDSIwi6o+WpzcU+Dhe5D3Rzc30l9tWuaOHNyLLMvzwC9hYYGy6gef6DNxe3ujZ0OJ7Nio/jL5FxxrEmYUDod0GGD+AikWP+8BWVF4Sd7qik19xCvUbNiRgFJERpG+45KyYvmvB/ks/XNBhr3d7HjvSaMld1MOScFW7eT6KQIn13EheD2XG0z/651tVl8enIucwyj7x4QFa/lrG/nsufTZvZ93srh0g5aa1w38midCtZaC7psPVEZwXP+XbjpMdvpbrXhdCjs/KCJrlYbSDB5YTKTzkxG8tFJD9vMPfxodxUA38lK5tZs7++F7ysQ1/GhCq5oAqD53UYS9neikyQSFIUdERIJF6SRPyfeZ7MONlnmu7sqOWixkqXT8sL0fKLUYsQZDuRdZpSPja7/kEA6NxXVNN8vJz16qIG3mzvQShLPTctnfJT3lkx1URpOvSGbg5vb2bXeSMPeLhr2uu4gR4LZl2WQNzvea88nBL+PWs3cW1EHwK8LMrgiNd5rj61SSUw9N5WU/GhK1tbT2WzjwHOHmRmpRidJKIpC84RYUi4JrjZq4eBwaQel6xqgT9sDfZyGeddkkZzru/Pjanpt3FR2CIuscG5iLL8rzPLpqoG8y4yywej6O/14HR+qsK4mehosJOzv9CSAJEnM0KtJz4r0WVIoisLP99XwdUc3MWoVL80oIFUnZjTCgdLpcF0MPF8AZYPRNfLzoVUNrTxd7XreP07M4bQE7x8ZIUkShQsSOWVJdv9vKLBtXQM9ZrvXn1MITru7LHy//DAysCQ9kR+PS/XJ86Tmu2aLs3J0zIxU97uOx+/vpGeUKy3C8PSY7ccVdQCnLsv2aVHX6XBy485DNNscTI6O5F9T89D48PxXpdPhGpy7/04/XceHI6wLO2ut5bgCTiVJWOt8d0F4qqqJVxrbUUvwzNQ8JseIJYOw0WE77qKHAkqHzWdP+UV7J7/Y57o5547cNJZkJPrsuQDUmuMvqIoC3W2++xuF4NFktXPjzkN0O2VOj4/hDxOzfTpzEmFXKJYlv1/HheN1tw5wfQPsFtlnz+mQFb5Xfpi93b2kRmh4cUYBsRrf3QGr2GTkj5oG+AZgCp5rXFgvxeqy9SiK0u+ioCgKUW12FFnx+l6AN5raebTSdXPE74uyOTcpeKZuBT+IjwCJ4y5+Skk7SrIOSe/dC1JFdy+37DqMQ4ErU+O5K9/3S1PRScf/jZIE0YniaIpQ1+10cmPZIeqsdgqjdPxnWh4RPjw7TN7XibLBiNqmHHcdlxUFXRD1mQ0H/n7vK4rCPRW1fNLWiV4l8cL0ArIjfXedUYy9yO81QccAqw8SYAiea1xYz9hFZehpnxCLrLgy0X1x0OzpRH6tDsWLy0dbTN3cvrcagO9lp/DtLN9u7BSCjxSrQTo31XUR8HwRqLYgr6xB8eIMQ4vNdXOOyeFkTlwUT04a55ebc6LitMy+LAP3U0kSFF+WIW6gCHFOReFHu6vZ2WkhUavmpRkFxGt9M2+g2GXk9U0o7zeBTYHMSEwF0Z7ruKwodEyIFTdQ+Jm/3/vP1DbzfH0rEvC3KbnMivPNcq+iKMjbOpBfqXUVdTEapDnxR6/jR/bYBdMNFMETSYCkXJJOT7EBa50FXZaeSPORfVD1vcgra1Cdn4o0yqbnVRYr3y6rxCorXJgcx/2FmV6KXhhrVNPiUHKjXNP2hgiwOJHfa4QOO/JrdUjzE5HmJ4xqtrjX6bo5p6rXxrjICP43PR+9H2/OyZsdT2phNN1tNqITxV2x4eChg/W812IiQpL477R88vS+OWRdaba63i/tdtcH6rwEpAWJJKkkehosnut4iijqAsJf7/0PWkw8cKAegN+Mz+SSlHifPI9icSKvb4Ij/WwZH+2qCSLVKDPjPdfxYCrqQBR2gGvmzjO6ywAlLRL5/UZosiK/04g0w4B0ZhLSCE7I7nTK3LLzEK12B9Nj9Px9cm7Y9KsTBibFasB9IYjVoFqag/JpM8qeTpRNbSi1PaguTB/RxcLdc3izqZs4jevmnJQI/xdWUXFaUdCFiRfrW/hHjes8xCcnj2NBvPdvzlEUBWWnCeXzFnACMWrXeyT7aAHX7zouBIyv3/t9ew7fmJnE/+Wk+OR5lNoe5PeboNsJaslVA8w42s+233U8yIT1UqybqquRiLpvUHW59r9J8VpU12W7plsBZacJ+eValGFuALcrCj+vM1HRYyVDp+WFGflE+3BjpzA2HJdvESpUi9KQLkwDrQR1vcgrq1EOdQ/7sf/Z0u3qOSzBf6bmMyFadIIQfOfrbqun5/Av8tK5Ji3B68+hWJzI6xpRPj1S1OVHoVo2rl9RB8e/r4TA8OXr0GR3enoOL0yI5eEi79+co8gK8tetyK/Vu4q6BC2qJdmoZvruCDRvC85y04/0e9dg+Ow+JEVGkVSYzvotlkmLkdQS0hnJKDl65A+M0GpDXlWDtDAFaeqJG1iDa4T5cKOZzT02otQqXpyeT4YueDZXCoExWL4BqCbFoqTpXKNEoxX57QakmQakM4Y2W7zOZOGZVlcx+NiEHM5M9E3LHkEAOGB1cFedCacCi9MSuDMvzevPodRZXO+HLgeoQToj2fWeOOb6e6L3leA/vnwdup0yP6ntoNHmYEJUJP+elofWyzc4KmY78gdNUO9qVCBNjUNamIw0xjqbjK1ovUzV1YhhoysJASRFxvDZff1GGlJuNKrlOZCjB4eC8rER5f0mFKvzhI/9fFsPb5h6Pf3qpsX67hwfYWxQdTV6LnowSL4lRKC6PhupOB4AZYcJ+ZValPYTzxZv7bHxYIMZgB+OS2V5ZpJv/ghBAFodTn5S006XrHCKIZo/Tcrx6myGIivI37Qhv1bnKuritaiuz0E16/hZE9d1/DcnfF8JvjeU69tIORWFu+tN7Lc6SNZqeGlGPnFeXv1SDnQhr6xxFXUREtJFaa79dGOsqIMwL+w0psNI9D9jR1JkNOaq/l+L1qC6OhPp9CRQgbLflQBKY++Aj/uRuZenml0n7/8sNZZFyQbf/AHCmKIxHfZc9NwkRUZjOibf1BKqs5JRXZEBehU0u2aL5d1mFOX4g6KqbA5+VtuBAzgvVsc9Bce3sBMEb+mVFX5aa6LBIZOjVfPc9Hx0XjzWROl0IK+tQ9nUBgpIk2NRLc1BSj3+hgxVtxHDxnuQjjlDaKDruOBbg13ftM1lo37sPxk7+bzbhk6CF6bnM86LN+coDhn5k2bkdxrBKkOaDtWycagmjt0Vj7Au7ByGPBSp//8CBdDvXYvU297v65IkoZqbgGpxNsRpwOxAfrUWuaS934ftLoudextcDY+XxOtZlihm6gSXgfINIKbkaTTtB477upQfjWrZOMjWg11BWW9E+aAJxXr04tnhlPlJTQcmWWFqpIbfZRhEz2HBZ2RF4TcNJsp67cSpJP6SE0+iF481UQ51I6+shrpe0EpIF6a59p9GHPO+cdqI3v4fUlZfRGTtl8c/jqTCEZfrtbiEkxvs+hb/ya+I2vk8OEd2fNjL7T2sancdBfW7DAOzvdBz2E1ps7n2z+90fWZLc+Jd++sNY/vGr7Au7OSYdExn/daTjEfavhFV8SapL19EVPlKkPu3CZEyIlEty0EqigEZlC9bkd+oR+l2UG93ckdtB1YFTo+O4OdpY7fiF7zv+HyTUCQ1usatJL96JbFfPYJkNff7HSnmyGzxqYkggbKvC3lVDUpTL3ZF4We1HVTbnaRrVDyZHY/eh610BOFvLV181GlFA/wpO57cCO8UdYpDRv60GfntBuiVIVWHamkOqknHX0N11Z+R8uoVxG36Ayp7D7bUGXQW/+Do++rI3i45RvSK9afjrm+SCkdMBip7N4avHyF5zVVEDFCEn8jnXVb+0NQJwI9SYrggzjs3gymKgrzLhLyqBlptoFejuioD1RnJSOqxfw0N+5snLJMWY80+A425CkdcLhpzFXFf/h5t234MX/yWqN2rMZ9+D7bM+Z7fkXRquDgNxulRNrZAtQXz6hpuPz2GVkWmSKfhsUwDGjFzIhzj2HyTZDtxXz9G5OGPiCl7Hn3F23TO/ymWideAyrWHRFJJSPMTUbL1ro3kJjvOV2p58BwDpRqZaJXEX3ISSBZ3XAs+9GaHhedaXed53ZcRx9wo79wMprTbkN9thBbXPlKpOB7p9KTjPmDVpsPEffUIkdUbAXDqk+lccCeWCVeBpKJnyg2e95Uo6gLj2OubHJWCft9aYjc/gbbjIEnv3EJv3vmYT/0lzricEz7Wvl47v6o3IQNXGiL5rpdWvxSrE2VDM8p+13YpcvSoLkxDig6dcih0/pJRkGPSsR25ENhi0mm5di1Re14hdstTaNv2kfT2TVgKLsZ8yi+QY12HC0uShDTNgJKhx/ZeI3fnqTmgyCTJ8FSGgWg/HggrjC198w2g/cK/ElHzBYavHkbTcYj4z35D1O6XMZ9+D/b02Z6fkzL1qJblIH9s5D+ynXc0MmpF4bGkWAp14q0s+M7mbhsPNbpmk29NiuZyw+jPi1MUBWV3J8qnzeBQQK9CdUEaUn7/pTbJ1kXMtn8RvfN/SLIdRaWhe9pNdM25DSXi6Jl5x76vhMA49nWwTL6O3oILiS35K1HlK4g8/BG6ms/omnkL3bP+H4r2+ILNaHdye20HPbLCvCgt96THeeXmHKWx13XAtdkBKpBOTUKaM3aOMRkqUX0MRKWhZ+oyjDe8T/eUpSiSCv2h90h95RJitv6d/9/efce3VZ+LH/8cSbZsy5K8VzwTx4mzB2QSMhiBcJlhdEChFy7QlvKjdAAttxTaSwu9t9D2Qhm9hVso81JWS4BAEshexFnOjme8l4ZtSZbO+f0hW4ljZzhxbPn4eb9efoElWTonevTVc77rwX/MoomECP5zkZ31SSbMAY2nNraS8l41mqP/ypEJ/fNlXUD99e/jnP0gamQskQ27SXr/G9hX/ARD69Gi00qUkU8viOO50cHJwz/e42XW32vRytsG69CFzpV4/fzoSHBxzqVWM99JOvs5TppXDc4X/awumNRlRgf3pjs2qdM0ovd/QPKblxNb9CKK2oEn8wLqr/8A1+yfdEvqRHjTzDacc39Kw/Xv4R0xGyXgw/rVn0h+cwlRB/8Jx8xTb1c17qtsodavkhtp5D9HxBFxlomXpmmoW5pR364MJnU2E4brMzGcF6+7pA4ksTspLSoe57xHaLjuHbxp56H4PVi3/IHkt67AXLIcNI2/Nbfxf04PCvB4RBTjvAQrVrxWgbrfNdinENLW0kTtgb20tTQN9qEMeyd8L4wRtE66jfqbPqZt7PVoKMQc+IDkNy7HUvQiBHwUtfn4RWfPyTejzFzfqkBbAPXdKtS1jWiBnqtmhThTTf7g3mEuVWNSdASP9cPiHK3WE5wnus8dLAs2OyE4jzT2aK+zqX4Xie9/nbiVP8HYVo/flk3T4mdpXvIigfiRvT6vtHHh4WTvgz9hNE1X/IXmS/6A3zoCY2sN8Z//kIQPbsHUuJeApvHTKgd7vH7ijAp/yIzDdpajX1qrPzgPfm0jqKCMjg3Ok0/X7+btMn5zGvxJhTRd9QpRh5Zh2/AkJtcREj79Ph8W3M7v0m8BFO5LjmVRogUtzRKcB1XtQVtWi1rejnpTCobIwZv/dHD9aja99UrwqkhRmHnTLYyaNW/Qjmc4O7RhNRvfeIXgPg69vxdqTBKO+b+irfAmbGt/RWTddmwb/4v6Q19w//gn8GFifqyZH4ywYxhhQ/uyAW2XE21LM1plG/5vJWBK0G+jJQaGV9W4/0gLlR0BRkQYeWpEHOazWJyjqRruNUdQP64ElWA5vctSUTKODusa2huxbnqa6L3/h4KGaorBPe1uWifdBsYTz+nbv2YlW955Xdq4QXZow2o2vnmK7xpFwTPyUjzZFxK7/S/EFr2AuWYLSe9cx8PnPcWqmClEKPDUiDiyznJxjmd/M+obFdAWAJMS3Gx4fP8M64YzSexOl6LgyV+CN2cBlqIXOXTwS+5NuQENha979vAty9Tgw2wRGK4fEaz5uakZbbeTuj9uI+EbhUSm998y7VMJ+DuoO7if0m2bKdl4zEokTWPTm6+SPnY8MXEJA3Y8Inglu/HNzqQOQNPY+MZfMUZEkjV5OkZT949jR8pEGq95negDHxDY8hy35d5LMybGe6t4IiUOoxIX3BLiohS07BjUz+qgxkvtH74i/rrRxEw6NzUUhf5pmsYvapxsb+8g1hDsOUk4g1rZXQJuH01v7ce7v3MbqXwLhouCxdSDD+ggpvh1rFv+iMEXHOloz78S56wfoVp6r2jhbqynqngX5du3Undw37EHL23cIAi1b1r39s3jdJF7/kws8cdtmm6Kwj39u7QVXINtw295u93E/8RMAeA3/p1MiVpwxseiBVQcn5bh/iJY7o7ESAyXp6EkDo/qT5LY9ZEWEcPBKffwLftNtKsK85s28+SuBzBst3WuZlyKYjCizE4MrmL8pBZ/fTt1z2wj7oqRWGaln7OrhbaWZqr27KRq905q9u/B7/P2fg6aiqu+Xhq9Aeaqr+s2l6TLulf+jOmtV0gbM44R4yaSUTiRaHtc8E7FgHP01dwTOZcDHpV0bz2vfvV90jY5aJ34LdzTvoMWGRscXugsR6ZVe2h6bS/egy3Y/2XkoPYWi6HpTw2tfOz0YAL+c4SdkWexOMdzsJmmN/ehujrAZEC5MBFlwtFek8jKddjWPU5E516OHUnjcMx9uNvCIQA14Kf+8CGqindwpHgnztrqE76mtHED70Tt2/aP3mX7R+9iTx8RbN/GTSQpdxQGY7BdUq0ZfDT7NzxUGUz6f1z6F75Z9r90FBd07kgxs0/H4W9sp/GNfXRUBC8QlIk2lAuTTqsso15IYtdHrQGV/1fZQr2qMCrSyBNZiVCRi7HlEHFf/hxL8RuhRknJCharjlztwLO3iZb3D+E52ELC0tEYYs5+A0RVVWksO0xV8U6O7N5BS1Vlt/ujbXZSRhVQVrSl2wdOUQxYk6U3Z6BZk1NAUXo0fmZLLN5WN5U7tlG5YxsA8ZnZjBg3ifTCiTxrTmSTRyVaUfh9RizxVeNRyr8gdvv/EH3gA1wzf0j76KuCvcVLR2DZ5cP1RQWtm2rwljlJ/PpYItIGrrdYDG0fOtp5sbPm8E/TbMy0nNku/1pAxbm8HNcXFaCBKSWGxG+MpUELzhE1Ois7t/pZDkAgKj50cdy11Y/H5aRqzy6qindQvbeYDk976PkVg4HkvHyS8kZR/PnH0sYNsl7bN0UhITOH5soyHNVHcFQfofjzj4mIjiF97HhGjJtIW/44Hqj3EkDhCmskt2floNbYiWjaT+KHt/bYkeJk2rbX0/z3A2jeAEqUiYTrR9OSOPwWMkpi1wf+Y+rVJRgN/CEzHnPkbBquf4+Y3a9h3fpHIhqKSXr/G7SNvgrXzB+iWlJJvHUc7rVVOJaV4NndSG2lm4Svj8Gc2/dSY97WVqr37qKqeCfVe3fjbXUfvVNRSMzOI2PcREaMm0h8ZjaKopA2ppBNb76KpqkoioEZN90sV7KDICYugZk33dLjvRg5Yy5NleWhnoim8lKaK8tpriznf+ocfDlrMYqm8UNPHTmWUTRf/jzmslXY1v0ak7MsuLP77tdxzn2YjpSJ2C/LxZxvp+nNffhr26j97yLirhyJZUaa7ueWiLOztc3HY501h29LiOHauDPb1sTf5KHpjb34yoO9JpYZaaHeY6W8BkvRn4nd/meUgA9NMdI2/hu4zrsHNcLa+VnYSVXxThorSrslCmZLLOmFE4IXPWPHERkTvGCxJiVLGzfITtS+jZo1D2+rm+q9u4Pv655d+NpaKd+2mT17i3n1urtxW+MZ43Fyj8VM27iv4xl1OdbNfyRmzxtEH15GVPlK3FP+Dffk23t9bdUXoOWDQ7RtCe4gEJljI+FrYzDFR9FSfeKeXb1StN6KTw4ip9OJ3W7H4XBgs9m63Vc9yG/QE7VO3mhux6zAi9kJTIzu3usWnPj7FNF73zk68Xf6d7Bd8gCYzPiOuGl6fS/+hnZQwHZxDtaFWSjHTEg+/hw1TcNRfYQjxTupKt5BQ8mhbiXMjr3ySS+cQFRs79Uu2lqacNXXY01OPicNXnp69/qkJ3sfw82JjvVcxdup3ot2l5Pq4p38o6aRl8YEhyEuWvMPpu3aEOqlyBg3kYyxY8moW451258wdAS3O2kbs5SYK38DsSkE3D6a396PZ19wiCN6QiLx13XvLR7sz9SZknjrf2U+P7eWNuFQNS6ymnky4/RWwB7/XrTt6Ow18QRQooxH53tqGuz+O4FlP8XYGiwM782YRcN5P6SyzkfV7p1U7d2Fx+no9nzxmdmhi9WE7DwMJ6hLK21c78KtfesaaSrZs5tfxudRGZ9KnKOBm//+PNHedqLtcWQUTiRj3AQyk0wkbnkSc/UWAPyxGZiWPAGFVwZ7BwFfdStNr+3BXx/8XrUuzMJ2UU5og+vB/lydqbOJN0nsTtPrTW08WRe8+nwyw37S0iYRdTtDqxkBSBgJl/0GChajev20vHeItm11AJhH2km4aQxGe3C4o7q6Gr/PS83+vVQV76CqeFePZeP2tIxgQzd+Ure5CoNpqDZ6EJ5ftLvaO7ijvAmvBlcpXq7Zu7HXeUWW+ERGFOSTz27yWz4hwqCC2QYLHoQZd6IpJtxrjuD4pBQCGsY4MwlfH4s5J3iew7HRG2zhGG8tAZVbS5so7wgwIcrEC9kJp12eruu9UH0BHP84TOumYNIWmW0l4Wtjgyu0a3bCsgegbC2aBg2mHPbG/QtltT7qDx9ADQRCz2cym0PzTdMLJxLTNd90kA3VmAvHeFM1jQeqHHzm8mJT4JfNh9F2baVm/x4CPl/ocQajiZRRo8lJMTG25UOSAhXBO/Lmo132G1oPx9Hyz8Pg1zDYIkm4aQxRo+K6vdZwbOMksTsNq91e7qtsQQW+nxzLvyaexnwlTSX6wIdYN/4nxrb64G35l8Blv4ak0bR+VUvLewfRfCqGGBPmS5Mpb9nDnvVrqT24F9V/tEatMSKC1NFjO69aJ2FJSDzBiw6eodroQfg1fFUdAb5V2kRjQOUCSyRPZcaFytO5G+uDvbe7d/aME5ORTGs7+ZHljLQ2YUvPDcZb/kX4Kl00vr6XQKMHDGC7JAfr/CxqamsG5RzPlsRb//GpGt+paOar9g7STAZeyU3oU3m69PR0OmpaaXx9L/7atmCvyfwsbJdko3hbYOV/4N/0EpWtVg63pXDIl4XT5ev2HNbkFDLGTWLEuIkkjxqN0RR+RdiHasyFW7wB/KHOxUtNbZiA57Ljmd5Zni7Q0UHdof2hESp3Q323v7NbIxkVWc5ISxsxEbfhVWcDEDU2gfjrR2OM7bnqdTgmdjLH7hSOr1f37dOtV6cYaC+4Gk/uRaTtfwXWPwsHl8Ozq2DWd4ia8wPUy2PwfdKAuS2K9veqaXYcpLZpDyoBLPGJZIyfSMa4SaTmj8EUOTyWaQ93roDKvRXNNAaCNYd/c1zN4djEZMbMW8SYeYvwe73UHOjq2d1JW0szZc2RlJHP57WQWNHKyOL7GVk4hoyvPU7q96fS/N5B2ovqcX5ShvdgC9qC+G4bw4rhRdM0Hqtx8lV7xxnVHNY0DffGalo+PAx+FYM1IthrMtKGa9UzlCz7C4eboihrnYlf63peHwajkZRRBWSMm0TG+InYknvf0kToz3st7bzUdLTm8PRjag4bIyJIHzue9LHj0a69CVddLUc627e6QwdwuHyUm6eTEXMVRtWGqvlp1z7AnlOAMWbsYJ1S2JEW/ST6o16dFhkLlzwGU79F6/sPUrJzOyVvfUrpy1vwBYwYMDIpYT5j7OdTYD+P7NQJ+BfasY3OlInuw0xH5/DEIV+AJJOBP2TGnbTmsMlsJnPCZDInTO42F7Ph4F6q9u2h0Wuh0Wth85o2zOvuITc3hZEXf530fxlL2ydVeA85oMrVa31OMTy82NjKP50ejASnmPSl5rDmCaB+XkfLweAKWnNBHN5pBrZ89iyHN62mvtUEjAg9PjY+gbyp5xGXm09aQSERUbKJ9nCzsdXLf5xmzWFFUbClpmFLTaNw4aV0tLXhXn4Ya6kZBQVXRxPr6z6g2VcLzx0g5dUPGDlzPnkLriEtfzQGw+BPURosktidwLH16vLOsF6dpqo0VZZTsmYlJds2U3PIA4wJ3R9t9JGXYmDExYnYRozE/Y9yolqj4DMfWsCFUhi+3fuif2maxhO1Lta3+ohS4PeZcaRFnH7DpCgKcRmZxGVkkn7Lv9LudlG6/StK1n1Oyfav8HSY2He4iX0vPAPAyFHTmRxzIZFtkagfVKNMjUOZmxiacCz0b5mznT81BJOyB1KtzIk9/W1NtKr2YIUdlx9N0ai2lrN5zX/j+aS18xEmQCM9PYGRF1xO3vSZpOSORFGUITs0Js7OYa+fHx1x4AcWW818tw81hzW3H8MnzdgqgxcD0VOTMU1OpXCni8Orl1FT20Kd20jd52vY8PkaomNjyZt6PnnTzic6dQSRMac50qYTktj14mzq1XV42qneVxxc4bVnJx6Xs9v9KXmjGDllKiNNZaTuewFDhxu2rIfAzVju+BnV71RBZTvap3Wo5e0oC5NRIofPxorD1avNbbzT0h6sOZxhZ1zU2c0xio61Ujh3PoVz56MG/FR/9hIln/yVww1G6r2xHD60lVKliMnxCyiwn4e2rQXvwUZMV4wgIlWKq+tdUZuPX3Rua3JzfAw3xJ/eF58aUGlfVY55V0dnr0lzZ69JcK6m2dBBbmwLIydOIPeGh4lJzztn5yCGjq6aw25VY3J0BI+m2097REo73Iq6vBY8KkQoxF87Gsu04NB92tgCZt9wM21Vhyh585cc3r2H0tZ42t1uilevpHj1yqM7CRROIGP8JOxpGbofDZPErhdP17lZ5fYS2VmvLvMk9eo0TcNZV9O579IO6g4dRFO7r/DKnTSNkdPOJ2/KdGKPXfjgvBM++wXseAO2vYqx+AOsU7+DO/NKtI0taHtdaNUeDJenoqTKsIVerXR5eKouuB/hD1JiWWjt3/faYDQxYvG/MWLRzVyw7r9xrfgDJS0xHG5NYKdjBbWeMmYkLcHsiqbjtXL2RezCNCFB5j7pVIXPzw+OtODTYEGsmftSTp7Id83lrN+5h/Qj6SRHjAAUyty72dLwKXFRbs5PrGVkbBMZBeMxXPECpE8emJMRYa+r5vCRzprDvzvNmsOaX0Nb24BW1Ln9TYoZw2WpWMb1bJNiMkYx/gcvM/7IVwT++ROqDhRz2J3A4fZUmtojqDu0n7pD+yn6x9+D89c7K2Ckjh6ry/nrktgd5+3mNl5tDk7s/EW6jSkxPd/0QEcHtQf3hTbRdDd2X7ljTU4N7buUPGo0mVnZvb+YLR2uex7Ovx0++jFUF2Hf8ASWuLdpuvARfF+lgqMD9a1KlLlJKFNP/ypHDA3Fng5+WuVAA5bGRXPzafacnJGIaJj/Y6yTv8ak5T9n0u6/41cVKjoOcciskOyeTaIpncLANErW7OTj9x8lOik+OMF93ERSwnS1ojh9zkCw56QloFFoNvF4hh1jL23K8auvUyOymZF8BVERMfjVDg6birAlfMa37VuwRXjBmgGXPgUTlob2FxNC1TQeqXb0ueaw1uxDXVYL9cGymMoUe/A70HSK2BoxDeMdn5K18y2ylv+c+e4NOHxm9kfM4KB/JDWl5bQ2N3Jg7SoOrF01JHacOBOS2B1jndvLE7XBveq+k2ThctvRiZ1tLU1UFe/iSPGO3vfayS8IBYc1OaVvL5w1A/5tJRS9SuDTRzC1HCZl4620ZS2hOeEetDIVbXUDWkVbcKJ7zPCdFKonNZ2LczwazLZE8kCqdWAS97gsuOElOP92tA/uJ69pH3n8BW/mGqqjf4axJIE860QSo0awvu599n3xGfu++Cy4v1hBYSjRC5f9xcTp6dA0fnSkhVJfgFSTgacz40J71QXrsB7sTOZ24KwLDq0aMHQu7poRfI5YP/as5SwsfR5FC6CZImDuj2De/RApC3BEd881tPKJy4sJ+K8RcadVc1jd40RbWQ8dGkQb+r64y2CAyV+DsVfAl7/Ftv5ZzldXc555PS1X38wh60KO7D/QuZNAU6iDZguvheUesWdCErtOB71+flLlIABcYYvi9vho6g8fpGrPCeqwhnbHnkjamEIizGc5fGYwwLRvUZ8wk9itz2DZ9SoxFR8RZfiM5rxf0F4+CUrbUF8rx7A4FSVreE0G1ZvWzp6TBr8arDmcYe/z4pyzlnsBDUvfIWbPW1g3/x5zy35yW27FmfOvOBuWYmtN4NLM2yi3HKaobDkep4PKnUVU7iwCTr8igBh8mqbxHzVONrd1EGNQ+H1mHFZPK4e/2smR4p3U7OtZhzU7dzJTzBcS5Qm2bRHZdaQ5HsJUEizb5Mm9COesB0gZ27ci7WJ4OL7m8AzLyYc8NZ+KtrIebW+wc4XM6OB33Zlux2S2wiWPUZ95Gbb1vyGq/Avii19mavQ/yJ/5Q9qufxxHTXW3qk6OmiocNVXsWfEJEVHRpI8dHxq2PVFVp3AkGxQD2+oauL/JRwtGxgU8/Nu2z2go3omvrfXogzrrsI7o3FsufkTWafeuHL/R4Ml0naOp+SC2tY9jPrIOAI95Go2Bn6G5gyvXlPPjUWYldCtHNpiG6uadMPAbeB7x+XmwysEuT7Dm8Cs5CWRE9t+V4ZnEm+JpDtVmVDQVvzGJBvMT+Js657PkxeAc56Py4MlreAZLnY3H7/Piqq/Dmpwi5Z2OM9DxVtrYyPP1bj4mCgMa99XuJ3nraprKS7s9zhxrDU4wHzeRDEMehrWOYK9JpIrd+lesrf8HgD9uJI45P8WXdQHQt3iD4blh7GAa6Hir7Qiw3OXh6To3AYI1h/9fysmTIq3WExx6dXSAQvC77bz4Xr/fzjTezGWrsK3/NSZHGQC+lMmh+tpwTB32Pbuo3rOrlzrsuaFNtOMzs2l3NIdtGzfse+yeXreJ/43PBsUImkba+s+o2rUBCNZhzRg7nozxk0gfO35AM3Z/fD5NV/wP5tLPsa3/NVGur0jXvkFT7EN43OehbW5Gq2zHcFkqik3mPQ0V77a081jN0ZXS19qj+zWpO1NaVDzOeT+nbdyN2Nb+CnP1FlJbb8dp+Qau9q9BSRu2eiMTFl/MxMuuxONyUrVnF1XFO6jeW4y31U3plg2UbtkQnGPVlfQpCjNvuoVRs+YN7gkOU89vKeI5SxoowV63MQd2YPz8bbqKFMZnZjOic2g9MTsX/KCtqkfb0wKAyVJFkv8hTK2NqJGxuKffQ+v4b4JR2hzR07st7fyyxknXJV+h2cT3k0+8OEfTNLRtLWhrG0EFrKbgd1rGife3O1PenAXUZ87BsvOvxH71LJF120l69wbaxlyHa8b9mC1J5E6fSe70mZ31bEtCm783H6mgsayExrISdi57n4io6KM93GHYxg3rxK60sbEzqeu8KlAUVs5ZwuJEO1MmhcEYu6LgzbuY+qwLiN3xErHbnifJ/wtaI+bRHPgBVIP6WgWGi1NQ8mWLinBX2xHglzXdt795qamVG+KjSe3DnnXnkj9xLE1XvkLUoWXYNjyJvfU1ok3radQeIeBOQv37EZTz4zHPTGDkjDmMnDEnOD+r5BBVxTup3LkNV33d0SfUNDa9+SrpY8efk6tacWKljY08b0nttphh76iJ3NBQyqTCQjIKJxJ9zDxJrc6LuqwGWjoADWvkW9j8f0NRVNrGLMU14weoMUkDfyJiSKjtCPCrY5I6gH1eP/V+tdf2TWsLBLcxKQ0uVmSUJfhdFnUO20JjJK1T7qB99JVYN/6OmAPvE7Pv70SVfIp7+vc6L1oiMRgMJOeNIjlvFJOvuJa2lmaq9gQXE1UfN20hHNu4YT0pZn9jc48VXJrBAOOnkjKqIHwmTpqicE/7DnU3LaN91BIsxtWkmb5DhPEAeFXUf9agrqhD86uDfaTiJMp9AY6f96ACFb5Abw8fPIqCJ38J9Td9hGvad4mIOEKq4S5ijMtBA21TM+o7R9BcHUBw8VBq/himXnU9M268pcfTaZqKq76+x+3i3Nrf2IymdG/iNYMBy9yLGDVrXiip0zQNdVsL6lsV0NKBwdBMcuSD2A2v0JE6kYZr38ax4D8kqRMnVe4LcPw30InaN60iOF+c0jYwKigLkzFckXZuk7pjj8uSimPREzRc/Tq+pPEYfG5s658g+f+uIbJiTY/Hx8TFkz/7Qi6843vMu/07Pe4PtzZuWCd2BYnxKGr3UFRUlYLEuME5oFNQY9Npufh3NF75V7SkOFJMP8JqfBtQ0XY6Ud+oRGv0DvZhihPIjjT2+MAZgKwwGIrtjRYRg/v8e6m/8SN8eReSEPF7EiKeRKEdqjyof6tAO+Tu9jfW5JQeF0uKYsCanDyQhy44vfZNaw+gfliN9mUDBCDKsJ7UiO9iiq2jZeFvaLzm9dAcJCFO5nTaN03VUNc1ov69CloDkBCB4WuZGCYNzlZeHWlTabzubVrm/4pAVAKmlsMkfnQH8R9/F6OjvNe/iUvLCPs2blgndrmJidztqQs1foqqcrenjtzE8N7Lxpcxg4al7+Cc9zDW2PdIivg5Bpqh0Yf6egXqLgdhtiZGAKkRRh5Os4U+dAbg4TRb2AzDnkjAlknz4j/SeMVfiEg8Qmrk94lQ9gd7i/9Rg7qyPtRbHBOXwMybbkHp7ClSFAMzbro5bIYohpNTtW9aZRvq38qgpA3wEWf6EwnmJ2mfciP1Ny2jveAaUIb1V4Tog1O1b5qzA/X/jqBtbgZAmWDD8LUslKTTL2V3TigG2sdeT/3XPsY98VY0xUhU2QqS37oC66anUDpauz18KLRxw3qOHcCd0yZzaWMj+xtbKEiOIzdxiOyYbjDRNv7rtI+6DOvmP5JS/P9o9v0/vIHpaJ/XQ6kbLklDMYd30jDcXBsXzRxLJBW+AFmRxrBP6o7ly5xDw/XvEVP8OkmbH8PVdi3uwFK0HQ60SheGKzJREiIZNWse6WPH46qvx5qcHFYN3nDTW/umqRraxka0Tc2AgkmpICHiSdScTBrmvE8gbuRgH7YYok7UvmkH3Kif14FXhUgDykXJGArCa/sQzWzDNech2gtvCO1IEbvteaL3v4dz5o/w5P9LqKcu3Nu4YZ/YQfDKNtx76U6kazWjadyNWNf8B1GVRTj8t6Idaoeq/fAvuedkhZE4c6kRQyuh68YYQdvEb+HJvwLrpqdJ3P1zmjt+gNoUj/q3EpQFySgT4omJSwi7xm64OrZ901wdaB8eRqs3Agoxxk+JTViGe+5P8GYvkKoR4qwd275pfhXtiwa0XZ2LxtLMGC5LQ7GH76rq43ekMLmOEL/ix/iK38Ax92f4k8YBhHUbJ4mdTvgTx9J81V+JOvwxiV8+QYvrXwm0p6O9XYFxKlBgwVhdg5qRjpaaMdiHK4Y4NToRx/xfEjFuJ/Ff/h531WK86lS0FU0Y9lXAjESMjbUSb2FAqa3CUFVNwBODtlVFC5hRaMNufhF11iQaJ74DRv3VyxSDQ3P5ocUX3Mrky0ZoDFZpUs6LQ5mViGIcAhcPvexIEVmzlaR3ltJWeCOu8+9DCXgxOUrx23NRY9MG+4i7kcROTxQFz6jL8WbPx7LlZTq2HqA9cCGBbcC2NgLYATfm8SvwX7xosI9W6EBH8kRarnuBqP0fEvnF27jaryVwxALvtku8hQHTZyvw7s4E7KHbIpT9xIzeSuuFD6NaehZUF+JMqbucaCvq6Lb8P8aI4dJUlJwhWC2pc0eKtoJrsG34LdGHPsKy502i97+PEvCioKEpBhwXPkb72OsH+2hDZGasDmkRMbTO/i7+Wy4hJnEFwU9Z11WSAe/uESi1VYN4hEJXFAOeMVfTftv3sYxcjcRbeFBqqzqTumObeRXjoijcl/+7JHWiX2kuf8+kDlCuTh+aSd0xju5I8QodcSMxBDwonSeqaCr2L3+OwV0zyEd5lCR2Oqbas/CPn87RL9kuRgxDtKyPCF9aZCwdmdOQeAsPhqpqejbxBlS/vbeHC3F2Wnw9kjoAxaef/VV9GefjnPtwj9sVTcXkLBuEI+qdJHY6p2akQ49tIwOofay3J8TpkHgLH/JeiAEVF9nzmk4B7Pqav+mPG9lz42/FgN+WM0hH1JMkdjqnpWZgHl8JdO3+HcA8/ohMaBfnhMRb+JD3QgwkxWpCWZRyNLlTQFmUgmLV11R+NTYNx4WPhZK7rjl24bSAQl//4qJX/osXYZpYhaG6GjU9HX+qTGQX547EW/iQ90IMJMMEG1pODDh8YI/UXVLXpX3s9XgzL8DkLMNvywmrpA4ksRs2tNQMAnKlLgaIxFv4kPdCDCTFagKdJnTHUmPT8IVZQtdFhmKHiTqXj60VLupcvsE+FDEMSLyFD3kvhOh/4fy50n9aLfhgVwNPrChH1cCgwAOLsrlqQtJgH5bQKYm38CHvhRD9L9w/V9Jjp3PVTi+/+TwYgACqBk+sKA/Lqwwx9NVKvIWNOpcv9OUD8l4I0R/Km9vDvo2THjsdcnj8bCh1sq7UwZrDjh5bC6kaVDq8pFj1tQxdDA6vX2XbETfrShysPNgs8RYmKlq8oS+fLqoGP/3nYS4rTGRuno10m3lwDk6IIeSIw8u6EgdrS51srXCGfRsniZ0OaJrGwYZ21pU6WVfiYFdNa48G/VgGBTLt0qCLM1fn8rG21MG6EidbKlx4/CfehFTibXBkxZkxKPRoC3bXtrG7to3/WgV5CVHMybMxJ9fOpPRYTEOhjqcQ55g/oLG9yt3Zxjkoa/ae9PHh1sZJYjdEtXcE2FLhYm2Jk/WlDurcHd3uH5kYxdxcO3Py7JQ2efjtyu7zAcLlykIMDQFVY1dNK+tKHKwrdXKwob3b/UmWCObm2Zida6fe7eOpLyol3gZZijWSBxZld5sL9K8z0ok0KawrcbKz2k1Jk4eSJg9/21pHbKSRmTlW5uTZmZ1jIz4mYrBPQYgB09TawfoyJ2tLHGwqd9J6TMUMowKTMmKZk2dnbq6NndWtPebYhVMbJ4ndEHJsd/C2She+wNFLcbNJ4bwsK3Ny7czJtZNmOxpkU0bEMifXRqXDS6bdHFYBKMLXsUP6G8qcOD2B0H0KMD7dwtxcG3Py7IxOikZRjvb2XDgyTuItDFw1IYlZOT0/+7ecl4bT42djmZN1pcGLQ4cnwOcHWvj8QAsKMC4thjm5dubm2SlI7v7+CjHUqZrG3to21nWOPOypa+t2f3y0iVm5wd7smTlWrOaj6VJeYnSvn6twIYldGDtVd3C6LbKz4bUxNdNKlOnEa2FSrJFhF3wivJxqSN9qNjIrx8acPBuzcuzERZ+4+ZB4Cx8nei9sUSYuGZPAJWMSCKgaxbWtrCtxsrbUwYH6dnbXtLG7po0XN1STZIlgVo6NuXk2zs+2YYk0DsKZCHF23N4Am8qD7dv6MidNbf5u949NiWFO58VqYWoMhpNczIRzGyeJXZhpau0IfrGWnro7ODchSq6ixVnpy5D+hHQLJoPEmx4ZDQoT02OZmB7LXXMyqHP7WN/ZDm0ud9HQ2sE/ihv5R3EjJoMSGgWYm2dHSs+KcKVpGqVNns5eOQdFVW4Cx0wHjokwMCM7eLE6O9dOkkUf0w8ksRtkqqqx84iDFXvrWLmvjh2Vjm73n6w7WIgzUd7Yxoq9tazYV8+GQw2nPaQvho+U2EiunpDE1ROS8PlVio64Q4tlKh1etlS42FLh4g+rj5CbWMrCsSksGpvCjLwEzCbpzRODx9MRYMPhRlburWPlvnrKm7oPsWbHmZmTZ2dOro0pI2KJMOpv1zfJEgaB09PBmgMNrNhbx6p99TS4uw+x9qU7WIhT8flVtpQ2hS4eDtW3dru/L0P6YviJNBmYkWNjRo6NH8yH8mYPazsX0RQdcVPa2MZLa0t5aW0pMZFG5uYnsWhsCgvHpJBmjxrswxfDQFVLOyv31bFybx1rDzbS3nF0PnCEUWHaiFhm59qZk2cjK07/MSmJ3QDQNI1D9W5W7K1jxd46tpQ24z9m8lKs2cQFnY1hYbymm+5gMXjqXB5W7atn5d46Vh9owO09OpfEaFA4LyeeRWNTmJioyJC+6JPs+Ciy46P4+rRUWn0BDrqMrNxbz8p9ddS5vCwvrmV5cS0A49JtwSRvbApTsuIwylC+6Af+gEpRRUvoO3Vvjavb/ak2c+jiYlSsn5hhNidUErtzxOtX+arSxbpSJxsr9lDR1H17iJHJFhaOCQ5fnJ+bQGRnL0l1dfVgHK4Y4o5d4bX5yKEeQ/qJlkgWdMbbBaOTsEcHLx4k3sTZsEQauWxCOpdNSEdVNYqrnaEv2+2VLRRXOymudvLfKw+SYIlkfkEyU1MjmJFjwx4lXz/i9Dna/Wwoc/LVqhq+PFBPS9vR+cCKAlOz4kIXEePSbaGL1eHYxsknqx/VunyhpdPHb9oaaTQwc2RC6CoiN8kyiEcq9MDtDbCxa4VXqZPm9u4rvCaOsIfmPk0aYccgvSXiHDIYFCaMsDNhhJ17LxpNg9vLF/vqWbGvji/319PU6uPdbUd4l+DeXxPTLcG5nHl2RiVKr7HormuV/tqS4CKe3cet0rdHRzC/IJmFY5OZX5BCgkXmA3eRxO4s+FWN3dWtoe1IDjV6ut2fHBvBnFwbV0zNZW5+EhZZ+CDOgqZplDV7Qg3d9uNXeEUamJFlY8mUbBaMSSbFpv+5JCJ8JcWaWTo9k6XTM+kIqGwta2blvjqW76ricKOH7VWtbK9q5U/rqkizRjK7c5HYeVlWoiJknudwdKpV+qMSo7h04ggWjU1halYcJh0ufOgPkmn0UVd38NqS4KatLu/RSZoGBcanWUILH7o2bU1PTxvEIxZD2bFD+utKHFQ5uxeazo43h7YjmZxhIcJoIF32nxBhJsJoYNbIRGaNTOS2KXFUO72s67xA2VLhosbl492dDby7s4FIo8L0TGto5WJGGJVqEv2vssUb2o7kqyNuOk6xSl/at1Prc2L35Zdf8tvf/patW7dSXV3Nu+++yzXXXBO6X9M0Hn30UV544QWam5uZOXMmzzzzDOPHj+/P4x4wp+oO7tq0dW6enZk5tpNu2ir6brjFGwSH9IOrDoNfel7/0YDrWuEV/NKzkxknX3r9aTjG22BIt5lZOjmZpZOT8XSobK10hcrV1bh8rC9zsr7MyX8BuQlRoT3z9FbPdjjGW0dAZXtVa2cVJQflx228n2GLDCX1skr/zPQ5C2ltbWXy5Ml8+9vfZunSpT3uf/LJJ/nd737Hyy+/TEFBAb/61a+45JJL2LdvH1artV8Our/VuXxUtHjJiguWBmnzBbuD15UGG5r6XrqDg5sE2xkvm7aeU8Mh3k53SL9rmGq4rfAaSHqMt3AXFWFgbl6wdJmmaRxu9LC+NFg6cWeVm9ImD6VNHl77KljPdkZOsAdndq6NhGPq2R7/uRoK9Bhvvb0Pja0doQ2vN5Y7aTt2430DTM6IDW25lBMv8y3PVp8Tu8svv5zLL7+81/s0TePpp5/mZz/7Gddddx0A//u//0tqaiqvvfYad91119kd7Tnwwa6GUDFfBchNMFPp8PXoDj4/yxaaAyKbtg4cvcdbYWoMFS3eEwzpBxu6/CSp0zlQ9BZvQ42iKIxKimZUUjQ3d9az3VTuZG2Jkw1lTlra/aw40MKKznq2hakxzMmz4w+o/HVLbbei7FdNSBrs0zklvcXb8e3b3DwbDa1+9vZSh3V255Slmdk2Ys1ysdqf+nXcsKSkhJqaGi699NLQbWazmfnz57Nu3bpeA9Hr9eL1Hu2KdTqd/XlIJ1Xn8vGbz8vpSuE0oKQpeCzHdgdPy7Rilu7gsHMm8QaDF3O9xVtxbbDBO3ZIf1aODbsM6YedoRZvemCLMnFxQQIXFwTr2e6pbQv1bO+vb6e4ti30GeqiavDEinJm5diGTM9db4ZavNW5fKGkDoLt25qSo69dmBLDnLxg58hY2Xj/nOrXb4+amhoAUlNTu92emppKWVlZr3/z61//mkcffbQ/D+O0VbR40Xq5/d8vyebywkTpJQlzZxJvMHgxd6J4+8GFmVw3OVmG9MPcUIs3vTEaFCakW5iQbuGu2RnUd9az/eeeRnZUda+mompQ6fAO6cRuqMVbRYu32/zzLjdPS+Vr01JIlI33B8w56YY6PiHSNO2ESdJDDz2Ew+EI/VRUVJyLQ+pVVpyZ4w/LoMB5WTZJ6oaQvsQbDF7MZcWZOT53MyiwID9OkrohZKjEm94lx0Zy1YQkfnlZXq+fq0ydrKYdKvF2ovbthinJktQNsH5N7NLSgtt6dF1pdKmrq+tx1dHFbDZjs9m6/QyUFGskDy7KDgVj19yMoXyVN5ycSbzB4MVcijWSByTehqyhFm/DhV4/V0Mt3vT6PgxF/ToUm5eXR1paGsuXL2fq1KkA+Hw+vvjiC5544on+fKl+c9WEJGbl2Kh0eMm0D53VVELiTQysoRhvw4UeP1dDMd70+D4MRX1O7NxuNwcPHgz9XlJSQlFREQkJCWRnZ3Pffffx+OOPM3r0aEaPHs3jjz9OTEwM3/jGN/r1wPtTijVSAjBMSbyJgaTHeBsuhuLnSo/xNhTfB73pc2K3ZcsWFi5cGPr9/vvvB+DWW2/l5Zdf5ic/+Qnt7e1897vfDW2o+Omnn4btnjsivEm8iYEk8SYGksSbOBcUTdN6W6g3aJxOJ3a7HYfD0WNuQHV19SAd1dnpSwkUvZzjyd7HcHOiY9XLe3EyejlHibfB09cST3o5z6EScxJv+jjPvsSbbM4mhBBCCKETktgJIYQQQuiEJHZCCCGEEDohiZ0QQgghhE5IYieEEEIIoROS2AkhhBBC6IQkdkIIIYQQOiGJnRBCCCGETkhiJ4QQQgihE5LYCSGEEELohCR2QgghhBA6IYmdEEIIIYROSGInhBBCCKETktgJIYQQQuiEJHZCCCGEEDohiZ0QQgghhE5IYieEEEIIoROS2AkhhBBC6IQkdkIIIYQQOiGJnRBCCCGETkhiJ4QQQgihE5LYCSGEEELohCR2QgghhBA6IYmdEEIIIYROSGInhBBCCKETktgJIYQQQuiEJHZCCCGEEDohiZ0QQgghhE5IYieEEEIIoROS2AkhhBBC6IQkdkIIIYQQOiGJnRBCCCGETkhiJ4QQQgihE5LYCSGEEELohCR2QgghhBA6IYmdEEIIIYROSGInhBBCCKETktgJIYQQQuiEJHZCCCGEEDohiZ0QQgghhE5IYieEEEIIoROS2AkhhBBC6IQkdkIIIYQQOiGJnRBCCCGETkhiJ4QQQgihE5LYCSGEEELohCR2QgghhBA6IYmdEEIIIYROSGInhBBCCKETktgJIYQQQuiEJHZCCCGEEDohiZ0QQgghhE5IYieEEEIIoROS2AkhhBBC6IQkdkIIIYQQOiGJnRBCCCGETkhiJ4QQQgihE5LYCSGEEELohCR2QgghhBA6IYmdEEIIIYROSGInhBBCCKETktgJIYQQQuiEJHZCCCGEEDohiZ0QQgghhE5IYieEEEIIoROS2AkhhBBC6IQkdkIIIYQQOiGJnRBCCCGETkhiJ4QQQgihE5LYCSGEEELohCR2QgghhBA6IYmdEEIIIYROSGInhBBCCKET5yyxe/bZZ8nLyyMqKorp06ezevXqc/VSQki8iQEl8SYGksSb6Itzkti9+eab3HffffzsZz9j27ZtzJs3j8svv5zy8vJz8XJimJN4EwNJ4k0MJIk30VfnJLH73e9+x+23384dd9xBYWEhTz/9NFlZWfzpT386Fy8nhjmJNzGQJN7EQJJ4E31l6u8n9Pl8bN26lQcffLDb7Zdeeinr1q3r8Xiv14vX6w397nA4AHA6nT0e63K5+vloB4bFYjntx+rlHLveP03Tzunr9jXe4PRjTi/vxcno5Rwl3gZPX+IN9HOeAxFzEm89SbydOt76PbFraGggEAiQmpra7fbU1FRqamp6PP7Xv/41jz76aI/bs7Ky+vvQxCBwuVzY7fZz9vx9jTeQmNMziTcx0M5lzEm8ieOdTrz1e2LXRVGUbr9rmtbjNoCHHnqI+++/P/S7qqo0NTWRmJjY6+PPBafTSVZWFhUVFdhstgF5zYE20OeoaRoul4uMjIxz/lpw+vEGgx9zEm/9T+LtxIZDvIG+Y07iLfyEc7z1e2KXlJSE0WjscTVRV1fX46oDwGw2Yzabu90WFxfX34d1Wmw2m64DEQb2HM9lz0mXvsYbhE/MSbz1L4m3kxsO8Qb6ijmJt/AXjvHW74snIiMjmT59OsuXL+92+/Lly5kzZ05/v5wY5iTexECSeBMDSeJNnIlzMhR7//33c8stt3Deeecxe/ZsXnjhBcrLy7n77rvPxcuJYU7iTQwkiTcxkCTeRF+dk8TupptuorGxkccee4zq6momTJjARx99RE5Ozrl4ubNmNpt55JFHenRf64mez1HiLfzo+Rwl3sKTXs9T4i08hfN5Ktq53h9ACCGEEEIMCKkVK4QQQgihE5LYCSGEEELohCR2QgghhBA6IYmdEEIIIYRODMvE7he/+AVTpkwZlNdetWoViqLQ0tIyqK/53nvvkZ+fj9Fo5L777huwYxmuhlPMSbwNvuEUbyd6TYm5gSPxFl7xFnaJ3W233cY111zT4/bBePPC0XPPPYfVasXv94duc7vdREREMG/evG6PXb16NYqikJGRQXV1dbddq++66y6uv/56Kioq+OUvfzlgxx+OJOZOTOKt/0m8nZzEXP+SeDs5PcZb2CV24uQWLlyI2+1my5YtodtWr15NWloamzdvpq2tLXT7qlWryMjIoKCggLS0tFCdQLfbTV1dHYsXLyYjIwOr1Trg5yGGBok3MdAk5sRA0mO8DcnErrdu36effprc3NzQ76tWrWLGjBlYLBbi4uKYO3cuZWVl3f7m+eefJysri5iYGG644YZuVy6bN2/mkksuISkpCbvdzvz58/nqq6+6/b2iKPz5z3/m2muvJSYmhtGjR/PBBx90e8xHH31EQUEB0dHRLFy4kNLS0rM69zFjxpCRkcGqVau6nevVV1/NqFGjWLduXbfbFy5c2O3KbNWqVaGgW7RoEYqihJ5r3bp1XHjhhURHR5OVlcW9995La2vrWR2vXgzXmJN4GxzDNd5AYm4wSLzpK96GZGJ3Kn6/n2uuuYb58+ezY8cO1q9fz5133hnKrgEOHjzIW2+9xYcffsjHH39MUVER3/ve90L3u1wubr31VlavXs2GDRsYPXo0S5YsweVydXutRx99lBtvvJEdO3awZMkSvvnNb9LU1ARARUUF1113HUuWLKGoqIg77riDBx988KzPb8GCBaxcuTL0+8qVK1mwYAHz588P3e7z+Vi/fj0LFy7s9rdz5sxh3759ALzzzjtUV1czZ84cdu7cyeLFi7nuuuvYsWMHb775JmvWrOGee+456+MdDvQccxJv4UfP8QYSc+FG4m2IxZsWZm699VbNaDRqFoul209UVJQGaM3NzdojjzyiTZ48udvfPfXUU1pOTo6maZrW2NioAdqqVat6fY1HHnlEMxqNWkVFRei2ZcuWaQaDQauuru71b/x+v2a1WrUPP/wwdBugPfzww6Hf3W63piiKtmzZMk3TNO2hhx7SCgsLNVVVQ4954IEHQudxpl544QXNYrFoHR0dmtPp1Ewmk1ZbW6u98cYb2pw5czRN07QvvvhCA7RDhw5pK1eu7Paazc3NGqCtXLky9Jy33HKLduedd3Z7ndWrV2sGg0Frb28/42MdCiTmTk7irX9JvJ2axFz/kXg7Nb3FW1j22C1cuJCioqJuP3/+859P++8TEhK47bbbWLx4MVdeeSW///3vqa6u7vaY7OxsMjMzQ7/Pnj0bVVVDmXddXR133303BQUF2O127HY7breb8vLybs8zadKk0P9bLBasVit1dXUA7Nmzh1mzZnW7qpk9e/bp/0OcwMKFC2ltbWXz5s2sXr2agoICUlJSmD9/Pps3b6a1tZVVq1aRnZ3NyJEjT+s5t27dyssvv0xsbGzoZ/HixaiqSklJyVkfc7iTmDsxibf+J/F2chJz/Uvi7eT0Fm+mc/rsZ8hisZCfn9/ttsrKytD/GwwGtONK3HZ0dHT7/aWXXuLee+/l448/5s033+Thhx9m+fLlzJo1q9fX7AqUrv/edttt1NfX8/TTT5OTk4PZbGb27Nn4fL5ufxcREdHjeVRVBehxjP0lPz+fzMxMVq5cSXNzM/PnzwcgLS2NvLw81q5dy8qVK1m0aNFpP6eqqtx1113ce++9Pe7Lzs7ut2MPVxJzJybx1v8k3k5OYq5/SbydnN7iLSwTu1NJTk6mpqYGTdNCQVNUVNTjcVOnTmXq1Kk89NBDzJ49m9deey0UhOXl5VRVVZGRkQHA+vXrMRgMFBQUAMFVMc8++yxLliwBgmP7DQ0NfTrOcePG8d5773W7bcOGDX16jhPpmsDZ3NzMj3/849Dt8+fP55NPPmHDhg18+9vfPu3nmzZtGrt37+7x4RdBwz3mJN4G1nCPN5CYG0gSb/qKt7Acij2VBQsWUF9fz5NPPsmhQ4d45plnWLZsWej+kpISHnroIdavX09ZWRmffvop+/fvp7CwMPSYqKgobr31VrZv387q1au59957ufHGG0lLSwOCGfwrr7zCnj172LhxI9/85jeJjo7u03HefffdHDp0iPvvv599+/bx2muv8fLLL/fLv8HChQtZs2YNRUVFoasLCAbhiy++iMfj6THJ82QeeOAB1q9fz/e+9z2Kioo4cOAAH3zwAd///vf75XiHuuEecxJvA2u4xxtIzA0kiTd9xduQTOwKCwt59tlneeaZZ5g8eTKbNm3iRz/6Uej+mJgY9u7dy9KlSykoKODOO+/knnvu4a677go9Jj8/P7S65tJLL2XChAk8++yzofv/8pe/0NzczNSpU7nlllu49957SUlJ6dNxZmdn88477/Dhhx8yefJknnvuOR5//PGz/wcgGITt7e3k5+eTmpoaun3+/Pm4XC5GjRpFVlbWaT/fpEmT+OKLLzhw4ADz5s1j6tSp/Pu//zvp6en9crxD3XCPOYm3gTXc4w0k5gaSxJu+4k3RztWgtRBCCCGEGFBDssdOCCGEEEL0JImdEEIIIYROSGInhBBCCKETktgJIYQQQuiEJHZCCCGEEDohiZ0QQgghhE5IYieEEEIIoROS2AkhhBBC6IQkdkIIIYQQOiGJnRBCCCGETkhiJ4QQQgihE/8fipQ94NXf8scAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for n in range(4):\n",
" socks_shuffled = socks.copy() # work on a copy of the original dataframe\n",
" for i in range(len(socks)):\n",
" if np.random.rand()>0.5: # generate a random number between 0 and 1 - if it is more than 0.5:\n",
" socks_shuffled.loc[i,'Husband'] = socks.loc[i,'Wife'] # flip number of socks for husband and wife\n",
" socks_shuffled.loc[i,'Wife'] = socks.loc[i,'Husband'] # flip number of socks for husband and wife\n",
" #else:\n",
" # don't shuffle the row!\n",
" socks_shuffled\n",
" \n",
" plt.subplot(1,4,n+1)\n",
" sns.barplot(data=socks_shuffled, ci=None, color=[0.9,0.9,0.9]) # ci=None switches off errorbars\n",
" for i in range(len(socks)):\n",
" plt.plot([0,1], [socks_shuffled.Husband[i], socks_shuffled.Wife[i]], '.-')\n",
" plt.xticks([0,1], labels=['Husband','Wife'])\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "ed9e5ef2",
"metadata": {},
"source": [
"We note that:\n",
" \n",
"
\n",
"
Different couples are randomly flipped in different simulated datasets (shuffles - look which colour=ed lines have reversed their slope)\n",
"
It matters a lot whether the couples with a big disparity were flipped - you can get a sense of this by looking at the grey bars which show the mean number of socks for men and women\n",
"
\n",
"\n",
"### Obtain the summary statistic of interest\n",
"\n",
"We are interested in the mean difference in pairs of socks owned [husband-wife]. For each shuffle this is obtained as follows:"
]
},
{
"cell_type": "code",
"execution_count": 81,
"id": "7544e3ef",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"mean difference for the last shuffle = -3.4\n"
]
}
],
"source": [
"mDiff = np.mean(socks_shuffled.Husband - socks_shuffled.Wife)\n",
"print('mean difference for the last shuffle = ' + str(mDiff))"
]
},
{
"cell_type": "markdown",
"id": "801dc763",
"metadata": {},
"source": [
"### Plot the null distribution for a large number of shuffles\n",
"\n",
"Now we can repeat the process for a large number of shuffles and get the mean difference in pairs of socks owned [husband-wife] for each shuffle. The distribution of these difference is the null distribution to which our observed difference (husbands own 6.6 more pairs) is to be compared."
]
},
{
"cell_type": "code",
"execution_count": 83,
"id": "4445af3d",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGdCAYAAAD0e7I1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvBklEQVR4nO3df3RU9Z3/8ddIQn5gEklGZpKamKiRgMGSBjYr6oJLCGVF6uIxdRHWH+iBgtEIHGxKdxk8NFnTFbIbKl1cJAiHxtPT0mLrD4JW1EYsRq0kBlpXcPiRmI6mmSDjJIT7/cNlvjuGoEwmuZOb5+Oce473cz/3k/cn7SGv+cz9YTMMwxAAAIBFXWR2AQAAAAOJsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACwtyuwCIsGZM2d04sQJJSQkyGazmV0OAAD4GgzDUGdnp9LS0nTRRX2v3xB2JJ04cULp6elmlwEAAEJw9OhRXXbZZX0eJ+xISkhIkPTFLysxMdHkagAAwNfh9XqVnp4e+DveF8KOFPjqKjExkbADAMAQ81WXoHCBMgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTeeg6YwO12y+Px9GsMu92ujIyMMFUEANZF2AEGmdvtVk7OOPl8p/o1TlxcvA4ebCbwAMBXIOwAg8zj8cjnO6WCe1crMTUzpDG8LUf05lNr5PF4CDsA8BUIO4BJElMzlZwx1uwyAMDyuEAZAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYmqlhJzMzUzabrde2dOlSSZJhGHK5XEpLS1NcXJymTZumpqamoDH8fr9KSkpkt9s1atQozZkzR8eOHTNjOgAAIAKZGnb279+vlpaWwFZXVydJuv322yVJlZWVWrdunTZs2KD9+/fL6XRqxowZ6uzsDIxRWlqqnTt3qra2Vq+//rpOnjyp2bNnq6enx5Q5AQCAyGJq2Ln00kvldDoD229+8xtdeeWVmjp1qgzDUFVVlVatWqW5c+cqNzdXW7du1alTp7Rjxw5JUkdHhzZv3qzHH39chYWFysvL0/bt23XgwAHt2bPHzKkBAIAIETHX7HR1dWn79u269957ZbPZdPjwYbW2tqqoqCjQJyYmRlOnTlV9fb0kqaGhQd3d3UF90tLSlJubG+hzLn6/X16vN2gDAADWFDFh51e/+pX++te/6u6775Yktba2SpIcDkdQP4fDETjW2tqqkSNHavTo0X32OZeKigolJSUFtvT09DDOBAAARJKICTubN2/WrFmzlJaWFtRus9mC9g3D6NX2ZV/Vp6ysTB0dHYHt6NGjoRcOAAAiWkSEnY8++kh79uzRfffdF2hzOp2S1GuFpq2tLbDa43Q61dXVpfb29j77nEtMTIwSExODNgAAYE0REXa2bNmiMWPG6Oabbw60ZWVlyel0Bu7Qkr64rmfv3r2aMmWKJCk/P1/R0dFBfVpaWtTY2BjoAwAAhjfT33p+5swZbdmyRXfddZeiov5/OTabTaWlpSovL1d2drays7NVXl6u+Ph4zZs3T5KUlJSkhQsXavny5UpJSVFycrJWrFihCRMmqLCw0KwpAQCACGJ62NmzZ4/cbrfuvffeXsdWrlwpn8+nJUuWqL29XQUFBdq9e7cSEhICfdavX6+oqCgVFxfL5/Np+vTpqqmp0YgRIwZzGgAAIEKZHnaKiopkGMY5j9lsNrlcLrlcrj7Pj42NVXV1taqrqweoQgAAMJRFxDU7AAAAA8X0lR0AQORwu93yeDz9GsNutysjIyNMFQH9R9gBAEj6Iujk5IyTz3eqX+PExcXr4MFmAg8iBmEHACBJ8ng88vlOqeDe1UpMzQxpDG/LEb351Bp5PB7CDiIGYQcAECQxNVPJGWPNLgMIGy5QBgAAlkbYAQAAlsbXWACAiMIdYQg3wg4AIGJwRxgGAmEHw05/PzXyiREYONwRhoFA2MGwEo5PjXxiBAYed4QhnAg7GFb6+6mRT4wAMPQQdjAs8akRAIYPbj0HAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWZnrYOX78uObPn6+UlBTFx8dr4sSJamhoCBw3DEMul0tpaWmKi4vTtGnT1NTUFDSG3+9XSUmJ7Ha7Ro0apTlz5ujYsWODPRUAABCBTA077e3tuv766xUdHa3nn39e77//vh5//HFdcsklgT6VlZVat26dNmzYoP3798vpdGrGjBnq7OwM9CktLdXOnTtVW1ur119/XSdPntTs2bPV09NjwqwAAEAkiTLzhz/22GNKT0/Xli1bAm2ZmZmB/zYMQ1VVVVq1apXmzp0rSdq6dascDod27NihRYsWqaOjQ5s3b9a2bdtUWFgoSdq+fbvS09O1Z88ezZw5c1DnBAAAIoupKzu7du3SpEmTdPvtt2vMmDHKy8vTk08+GTh++PBhtba2qqioKNAWExOjqVOnqr6+XpLU0NCg7u7uoD5paWnKzc0N9AEAAMOXqWHnww8/1MaNG5Wdna0XX3xRixcv1oMPPqinn35aktTa2ipJcjgcQec5HI7AsdbWVo0cOVKjR4/us8+X+f1+eb3eoA0AAFiTqV9jnTlzRpMmTVJ5ebkkKS8vT01NTdq4caP++Z//OdDPZrMFnWcYRq+2Lztfn4qKCq1Zs6af1QMAgKHA1JWd1NRUjR8/Pqht3LhxcrvdkiSn0ylJvVZo2traAqs9TqdTXV1dam9v77PPl5WVlamjoyOwHT16NCzzAQAAkcfUsHP99dfr0KFDQW1/+tOfdPnll0uSsrKy5HQ6VVdXFzje1dWlvXv3asqUKZKk/Px8RUdHB/VpaWlRY2NjoM+XxcTEKDExMWgDAADWZOrXWA8//LCmTJmi8vJyFRcX6w9/+IM2bdqkTZs2Sfri66vS0lKVl5crOztb2dnZKi8vV3x8vObNmydJSkpK0sKFC7V8+XKlpKQoOTlZK1as0IQJEwJ3ZwEAgOHL1LAzefJk7dy5U2VlZXr00UeVlZWlqqoq3XnnnYE+K1eulM/n05IlS9Te3q6CggLt3r1bCQkJgT7r169XVFSUiouL5fP5NH36dNXU1GjEiBFmTAsAAEQQU8OOJM2ePVuzZ8/u87jNZpPL5ZLL5eqzT2xsrKqrq1VdXT0AFQIAgKHM9LADAPiC2+2Wx+MJ+Xy73a6MjIwwVgRYA2EHACKA2+1WTs44+XynQh4jLi5eBw82E3iALyHsAEAE8Hg88vlOqeDe1UpMzbzg870tR/TmU2vk8XgIO8CXEHYAIIIkpmYqOWOs2WUAlmLqc3YAAAAGGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYGmEHAABYWpTZBQAXwu12y+PxhHx+c3NzGKsBAAwFhB0MGW63Wzk54+Tzner3WN3+rjBUNLT1NzhKkt1uV0ZGRpgqAoCBQdjBkOHxeOTznVLBvauVmJoZ0hgtB95Q465NOn36dHiLG2LCFRzj4uJ18GAzgQdARCPsYMhJTM1UcsbYkM71thwJbzFDVDiCo7fliN58ao08Hg9hB0BEI+wAw1h/giMADBXcjQUAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNu7EAAJbU3yem89BM6yDsAAAsxdfxiSSb5s+f369xeGimdRB2AACW0n2qU5KhifMe0aVZOSGNwUMzrYWwAwCwpIvHZPDQTEjiAmUAAGBxpoYdl8slm80WtDmdzsBxwzDkcrmUlpamuLg4TZs2TU1NTUFj+P1+lZSUyG63a9SoUZozZ46OHTs22FMBAAARyvSVnWuuuUYtLS2B7cCBA4FjlZWVWrdunTZs2KD9+/fL6XRqxowZ6uzsDPQpLS3Vzp07VVtbq9dff10nT57U7Nmz1dPTY8Z0AABAhDH9mp2oqKig1ZyzDMNQVVWVVq1apblz50qStm7dKofDoR07dmjRokXq6OjQ5s2btW3bNhUWFkqStm/frvT0dO3Zs0czZ84c1LkAAIDIY/rKzp///GelpaUpKytLd9xxhz788ENJ0uHDh9Xa2qqioqJA35iYGE2dOlX19fWSpIaGBnV3dwf1SUtLU25ubqDPufj9fnm93qANAABYk6lhp6CgQE8//bRefPFFPfnkk2ptbdWUKVP0ySefqLW1VZLkcDiCznE4HIFjra2tGjlypEaPHt1nn3OpqKhQUlJSYEtPTw/zzAAAQKQwNezMmjVLt912myZMmKDCwkL99re/lfTF11Vn2Wy2oHMMw+jV9mVf1aesrEwdHR2B7ejRo/2YBQAAiGSmX7Pzf40aNUoTJkzQn//8Z916662Svli9SU1NDfRpa2sLrPY4nU51dXWpvb09aHWnra1NU6ZM6fPnxMTEKCYmZmAmgWGhP4+h7+8j7AEAFyaiwo7f71dzc7NuvPFGZWVlyel0qq6uTnl5eZKkrq4u7d27V4899pgkKT8/X9HR0aqrq1NxcbEkqaWlRY2NjaqsrDRtHrCucD2GXpK6/V39LwgA8JVMDTsrVqzQLbfcooyMDLW1tWnt2rXyer266667ZLPZVFpaqvLycmVnZys7O1vl5eWKj4/XvHnzJElJSUlauHChli9frpSUFCUnJ2vFihWBr8WAcAvHY+hbDryhxl2bdPr06fAWBwA4J1PDzrFjx/RP//RP8ng8uvTSS/W3f/u32rdvny6//HJJ0sqVK+Xz+bRkyRK1t7eroKBAu3fvVkJCQmCM9evXKyoqSsXFxfL5fJo+fbpqamo0YsQIs6aFYaA/j6H3thwJbzEAgPMyNezU1tae97jNZpPL5ZLL5eqzT2xsrKqrq1VdXR3m6gAAgBVE1DU7AABrCPVCfC7gx0Ag7AAAwiZcF/FzAT/CibADDGF8ekak6e9F/FzAj4FA2AGGID49I9KFehE/F/BjIBB2gCGIT8/oCw+8BHoj7ABDGJ+ecRYPvAT6RtgBAAvggZdA3wg7AGAhPPAS6M3Ut54DAAAMNMIOAACwNMIOAACwNMIOAACwNMIOAACwNMIOAACwNMIOAACwNMIOAACwNMIOAACwNMIOAACwNMIOAACwNMIOAACwNMIOAACwNMIOAACwNMIOAACwNMIOAACwNMIOAACwNMIOAACwtCizCwAwtDU3N/frfLvdroyMjDBVAwC9EXYAhMTX8Ykkm+bPn9+vceLi4nXwYDOBB8CAIewACEn3qU5JhibOe0SXZuWENIa35YjefGqNPB4PYQfAgCHsAOiXi8dkKDljrNllAECfuEAZAABYGmEHAABYGl9jAUAYuN1ueTyekM/v711tAPoWUti54oortH//fqWkpAS1//Wvf9W3vvUtffjhh2EpDgCGArfbrZyccfL5TvV7rG5/VxgqAvB/hRR2jhw5op6enl7tfr9fx48f73dRADCUeDwe+XynVHDvaiWmZoY0RsuBN9S4a5NOnz4d3uIAXFjY2bVrV+C/X3zxRSUlJQX2e3p69NJLLykzMzNsxQHAUJKYmhnynWneliPhLQZAwAWFnVtvvVWSZLPZdNdddwUdi46OVmZmph5//PGwFQcAANBfFxR2zpw5I0nKysrS/v37ZbfbB6QoAACAcAnp1vPDhw+HPehUVFTIZrOptLQ00GYYhlwul9LS0hQXF6dp06apqakp6Dy/36+SkhLZ7XaNGjVKc+bM0bFjx8JaGwAAGLpCvvX8pZde0ksvvaS2trbAis9ZTz311AWNtX//fm3atEnXXnttUHtlZaXWrVunmpoaXX311Vq7dq1mzJihQ4cOKSEhQZJUWlqqZ599VrW1tUpJSdHy5cs1e/ZsNTQ0aMSIEaFODwAAWERIKztr1qxRUVGRXnrpJXk8HrW3twdtF+LkyZO688479eSTT2r06NGBdsMwVFVVpVWrVmnu3LnKzc3V1q1bderUKe3YsUOS1NHRoc2bN+vxxx9XYWGh8vLytH37dh04cEB79uwJZWoAAMBiQlrZ+elPf6qamhotWLCg3wUsXbpUN998swoLC7V27dpA++HDh9Xa2qqioqJAW0xMjKZOnar6+notWrRIDQ0N6u7uDuqTlpam3Nxc1dfXa+bMmef8mX6/X36/P7Dv9Xr7PQ8AABCZQgo7XV1dmjJlSr9/eG1trd5++23t37+/17HW1lZJksPhCGp3OBz66KOPAn1GjhwZtCJ0ts/Z88+loqJCa9as6W/5AABgCAjpa6z77rsv8FVSqI4ePaqHHnpI27dvV2xsbJ/9bDZb0L5hGL3avuyr+pSVlamjoyOwHT169MKKBwAAQ0ZIKzuff/65Nm3apD179ujaa69VdHR00PF169Z95RgNDQ1qa2tTfn5+oK2np0evvvqqNmzYoEOHDkn6YvUmNTU10KetrS2w2uN0OtXV1aX29vag1Z22trbzrjzFxMQoJibm600WAAAMaSGFnffee08TJ06UJDU2NgYd+6pVl7OmT5+uAwcOBLXdc889ysnJ0SOPPKIrrrhCTqdTdXV1ysvLk/TF12d79+7VY489JknKz89XdHS06urqVFxcLElqaWlRY2OjKisrQ5kaAACwmJDCzu9+97t+/+CEhATl5uYGtY0aNUopKSmB9tLSUpWXlys7O1vZ2dkqLy9XfHy85s2bJ0lKSkrSwoULtXz5cqWkpCg5OVkrVqzQhAkTVFhY2O8aAQDA0Bfyc3YGw8qVK+Xz+bRkyRK1t7eroKBAu3fvDjxjR5LWr1+vqKgoFRcXy+fzafr06aqpqeEZOwAAQFKIYeemm24679dVL7/8ckjFvPLKK0H7NptNLpdLLperz3NiY2NVXV2t6urqkH4mAACwtpDCztnrdc7q7u7Wu+++q8bGxl4vCAUAADBTSGFn/fr152x3uVw6efJkvwoCAAAIp5Ces9OX+fPnX/B7sQAAAAZSWMPOG2+8cd4HBAIAAAy2kL7Gmjt3btC+YRhqaWnRW2+9pX/5l38JS2EAAADhEFLYSUpKCtq/6KKLNHbsWD366KNBL+UEAAAwW0hhZ8uWLeGuAwAAYED066GCDQ0Nam5uls1m0/jx4wOvdQAAAIgUIYWdtrY23XHHHXrllVd0ySWXyDAMdXR06KabblJtba0uvfTScNcJAAAQkpDuxiopKZHX61VTU5M+/fRTtbe3q7GxUV6vVw8++GC4awQAAAhZSCs7L7zwgvbs2aNx48YF2saPH6+f/OQnXKAMAAAiSkgrO2fOnFF0dHSv9ujoaJ05c6bfRQEAAIRLSGHn7//+7/XQQw/pxIkTgbbjx4/r4Ycf1vTp08NWHAAAQH+FFHY2bNigzs5OZWZm6sorr9RVV12lrKwsdXZ28vZxAAAQUUK6Zic9PV1vv/226urqdPDgQRmGofHjx6uwsDDc9QEAAPTLBa3svPzyyxo/fry8Xq8kacaMGSopKdGDDz6oyZMn65prrtFrr702IIUCAACE4oLCTlVVle6//34lJib2OpaUlKRFixZp3bp1YSsOAACgvy4o7Pzxj3/Ut7/97T6PFxUVqaGhod9FAQAAhMsFhZ2PP/74nLecnxUVFaW//OUv/S4KAAAgXC4o7HzjG9/QgQMH+jz+3nvvKTU1td9FAQAAhMsFhZ1/+Id/0L/+67/q888/73XM5/Np9erVmj17dtiKAwAA6K8LuvX8hz/8oX75y1/q6quv1gMPPKCxY8fKZrOpublZP/nJT9TT06NVq1YNVK0AAAAX7ILCjsPhUH19vb73ve+prKxMhmFIkmw2m2bOnKknnnhCDodjQAoFAAAIxQU/VPDyyy/Xc889p/b2dn3wwQcyDEPZ2dkaPXr0QNQHAADQLyE9QVmSRo8ercmTJ4ezFgAAgLAL6d1YAAAAQwVhBwAAWBphBwAAWBphBwAAWBphBwAAWBphBwAAWBphBwAAWBphBwAAWBphBwAAWBphBwAAWBphBwAAWBphBwAAWJqpYWfjxo269tprlZiYqMTERF133XV6/vnnA8cNw5DL5VJaWpri4uI0bdo0NTU1BY3h9/tVUlIiu92uUaNGac6cOTp27NhgTwUAAESokN96Hg6XXXaZ/u3f/k1XXXWVJGnr1q36zne+o3feeUfXXHONKisrtW7dOtXU1Ojqq6/W2rVrNWPGDB06dEgJCQmSpNLSUj377LOqra1VSkqKli9frtmzZ6uhoUEjRowwc3oABonb7ZbH4wn5fLvdroyMjDBWBCCSmBp2brnllqD9H/3oR9q4caP27dun8ePHq6qqSqtWrdLcuXMlfRGGHA6HduzYoUWLFqmjo0ObN2/Wtm3bVFhYKEnavn270tPTtWfPHs2cOXPQ5wRgcLndbuXkjJPPdyrkMeLi4nXwYDOBB7AoU8PO/9XT06Of//zn+uyzz3Tdddfp8OHDam1tVVFRUaBPTEyMpk6dqvr6ei1atEgNDQ3q7u4O6pOWlqbc3FzV19f3GXb8fr/8fn9g3+v1DtzEAAwoj8cjn++UCu5drcTUzAs+39tyRG8+tUYej4ewA1iU6WHnwIEDuu666/T555/r4osv1s6dOzV+/HjV19dLkhwOR1B/h8Ohjz76SJLU2tqqkSNHavTo0b36tLa29vkzKyoqtGbNmjDPBICZElMzlZwx1uwyAEQg0+/GGjt2rN59913t27dP3/ve93TXXXfp/fffDxy32WxB/Q3D6NX2ZV/Vp6ysTB0dHYHt6NGj/ZsEAACIWKaHnZEjR+qqq67SpEmTVFFRoW9+85v6j//4DzmdTknqtULT1tYWWO1xOp3q6upSe3t7n33OJSYmJnAH2NkNAABYk+lfY32ZYRjy+/3KysqS0+lUXV2d8vLyJEldXV3au3evHnvsMUlSfn6+oqOjVVdXp+LiYklSS0uLGhsbVVlZadocAAw9zc3NppwLYOCZGnZ+8IMfaNasWUpPT1dnZ6dqa2v1yiuv6IUXXpDNZlNpaanKy8uVnZ2t7OxslZeXKz4+XvPmzZMkJSUlaeHChVq+fLlSUlKUnJysFStWaMKECYG7swDgfHwdn0iyaf78+f0eq9vf1f+CAISdqWHn448/1oIFC9TS0qKkpCRde+21euGFFzRjxgxJ0sqVK+Xz+bRkyRK1t7eroKBAu3fvDjxjR5LWr1+vqKgoFRcXy+fzafr06aqpqeEZOwC+lu5TnZIMTZz3iC7NyglpjJYDb6hx1yadPn06vMUBCAtTw87mzZvPe9xms8nlcsnlcvXZJzY2VtXV1aqurg5zdQCGk4vHZIR8N5e35Uh4iwEQVqZfoAwAADCQCDsAAMDSCDsAAMDSCDsAAMDSCDsAAMDSIu6hggAAWIXb7ZbH4wn5fLvdzgtqw4CwAwDAAHC73crJGSef71TIY8TFxevgwWYCTz8RdgAAGAAej0c+3ykV3LtaiamZF3y+t+WI3nxqjTweD2Gnnwg7AAAMoMTUzJAfWInw4AJlAABgaYQdAABgaYQdAABgaVyzg6+NWygBDDfNzc2mnIvwIuzga+EWSgwk/qAg0vg6PpFk0/z58/s9Vre/q/8FoV8IO/hauIUSA4E/KIhU3ac6JRmaOO8RXZqVE9IYLQfeUOOuTTp9+nR4i8MFI+zggnALJcKJPyiIdBePyQj53zxvy5HwFoOQEXYAmI4/KAAGEndjAQAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAAS+NFoBhUzc3NppwLABi+CDsYFL6OTyTZNH/+/H6P1e3v6n9BAIBhg7CDQdF9qlOSoYnzHtGlWTkhjdFy4A017tqk06dPh7c4AIClEXYwqC4ek6HkjLEhnettORLeYgAAwwIXKAMAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEszNexUVFRo8uTJSkhI0JgxY3Trrbfq0KFDQX0Mw5DL5VJaWpri4uI0bdo0NTU1BfXx+/0qKSmR3W7XqFGjNGfOHB07dmwwpwIAACKUqWFn7969Wrp0qfbt26e6ujqdPn1aRUVF+uyzzwJ9KisrtW7dOm3YsEH79++X0+nUjBkz1NnZGehTWlqqnTt3qra2Vq+//rpOnjyp2bNnq6enx4xpAQCACGLqc3ZeeOGFoP0tW7ZozJgxamho0N/93d/JMAxVVVVp1apVmjt3riRp69atcjgc2rFjhxYtWqSOjg5t3rxZ27ZtU2FhoSRp+/btSk9P1549ezRz5sxBnxcAAIgcEXXNTkdHhyQpOTlZknT48GG1traqqKgo0CcmJkZTp05VfX29JKmhoUHd3d1BfdLS0pSbmxvo82V+v19erzdoAwAA1hQxYccwDC1btkw33HCDcnNzJUmtra2SJIfDEdTX4XAEjrW2tmrkyJEaPXp0n32+rKKiQklJSYEtPT093NMBAAARImLCzgMPPKD33ntPP/vZz3ods9lsQfuGYfRq+7Lz9SkrK1NHR0dgO3r0aOiFAwCAiBYRYaekpES7du3S7373O1122WWBdqfTKUm9Vmja2toCqz1Op1NdXV1qb2/vs8+XxcTEKDExMWgDAADWZGrYMQxDDzzwgH75y1/q5ZdfVlZWVtDxrKwsOZ1O1dXVBdq6urq0d+9eTZkyRZKUn5+v6OjooD4tLS1qbGwM9AEAAMOXqXdjLV26VDt27NCvf/1rJSQkBFZwkpKSFBcXJ5vNptLSUpWXlys7O1vZ2dkqLy9XfHy85s2bF+i7cOFCLV++XCkpKUpOTtaKFSs0YcKEwN1ZAABg+DI17GzcuFGSNG3atKD2LVu26O6775YkrVy5Uj6fT0uWLFF7e7sKCgq0e/duJSQkBPqvX79eUVFRKi4uls/n0/Tp01VTU6MRI0YM1lQAAECEMjXsGIbxlX1sNptcLpdcLleffWJjY1VdXa3q6uowVgcAAKwgIi5QBgAAGCiEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGmmhp1XX31Vt9xyi9LS0mSz2fSrX/0q6LhhGHK5XEpLS1NcXJymTZumpqamoD5+v18lJSWy2+0aNWqU5syZo2PHjg3iLAAAQCSLMvOHf/bZZ/rmN7+pe+65R7fddluv45WVlVq3bp1qamp09dVXa+3atZoxY4YOHTqkhIQESVJpaameffZZ1dbWKiUlRcuXL9fs2bPV0NCgESNGDPaUIpbb7ZbH4wn5/Obm5jBWAwDA4DE17MyaNUuzZs065zHDMFRVVaVVq1Zp7ty5kqStW7fK4XBox44dWrRokTo6OrR582Zt27ZNhYWFkqTt27crPT1de/bs0cyZMwdtLpHM7XYrJ2ecfL5T/R6r298VhooAABg8poad8zl8+LBaW1tVVFQUaIuJidHUqVNVX1+vRYsWqaGhQd3d3UF90tLSlJubq/r6esLO//J4PPL5Tqng3tVKTM0MaYyWA2+ocdcmnT59OrzFAQAwwCI27LS2tkqSHA5HULvD4dBHH30U6DNy5EiNHj26V5+z55+L3++X3+8P7Hu93nCVHdESUzOVnDE2pHO9LUfCWwwAAIMk4u/GstlsQfuGYfRq+7Kv6lNRUaGkpKTAlp6eHpZaAQBA5InYlR2n0ynpi9Wb1NTUQHtbW1tgtcfpdKqrq0vt7e1BqzttbW2aMmVKn2OXlZVp2bJlgX2v10vgAQBEpP7eIGK325WRkRGmaoamiA07WVlZcjqdqqurU15eniSpq6tLe/fu1WOPPSZJys/PV3R0tOrq6lRcXCxJamlpUWNjoyorK/scOyYmRjExMQM/CQAAQuTr+ESSTfPnz+/XOHFx8Tp4sHlYBx5Tw87Jkyf1wQcfBPYPHz6sd999V8nJycrIyFBpaanKy8uVnZ2t7OxslZeXKz4+XvPmzZMkJSUlaeHChVq+fLlSUlKUnJysFStWaMKECYG7swAAGIq6T3VKMjRx3iO6NCsnpDG8LUf05lNr5PF4CDtmeeutt3TTTTcF9s9+tXTXXXeppqZGK1eulM/n05IlS9Te3q6CggLt3r078IwdSVq/fr2ioqJUXFwsn8+n6dOnq6amhmfsAAAs4eIxGSHfXIIvmBp2pk2bJsMw+jxus9nkcrnkcrn67BMbG6vq6mpVV1cPQIUAAGCoi9hrdgAAQHgM94ucCTtDAK96AACEgoucv0DYiXC86gEAECoucv4CYSfC8aoHAEB/DfeLnAk7QwSvegAAIDQR/7oIAACA/iDsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAAS4syuwAAABD5mpubQz7XbrcrIyMjjNVcGMIOAADok6/jE0k2zZ8/P+Qx4uLidfBgs2mBh7AzwNxutzweT8jn9ydJAwDQX92nOiUZmjjvEV2alXPB53tbjujNp9bI4/EQdqzI7XYrJ2ecfL5T/R6r298VhooAAAjNxWMylJwx1uwyQkLYGUAej0c+3ykV3LtaiamZIY3RcuANNe7apNOnT4e3OAAAhgnCziBITM0MOQ17W46EtxgAAIYZbj0HAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWRtgBAACWZpmw88QTTygrK0uxsbHKz8/Xa6+9ZnZJAAAgAlgi7DzzzDMqLS3VqlWr9M477+jGG2/UrFmz5Ha7zS4NAACYzBJhZ926dVq4cKHuu+8+jRs3TlVVVUpPT9fGjRvNLg0AAJhsyD9BuaurSw0NDfr+978f1F5UVKT6+vpznuP3++X3+wP7HR0dkiSv1xvW2k6ePClJ+vSjQzrt94U0hrflI0lSx/E/KzrKNmTHiIQawjFGJNQQjjEioYZIGSMSagjHGJFQQzjGiIQaImWMSKghHGN4W7/4luXkyZNh/zt7djzDMM7f0Rjijh8/bkgyfv/73we1/+hHPzKuvvrqc56zevVqQxIbGxsbGxubBbajR4+eNysM+ZWds2y24LRpGEavtrPKysq0bNmywP6ZM2f06aefKiUlpc9zhiqv16v09HQdPXpUiYmJZpcz6Jj/8J6/xO9guM9f4ndg5fkbhqHOzk6lpaWdt9+QDzt2u10jRoxQa2trUHtbW5scDsc5z4mJiVFMTExQ2yWXXDJQJUaExMREy/2f/EIw/+E9f4nfwXCfv8TvwKrzT0pK+so+Q/4C5ZEjRyo/P191dXVB7XV1dZoyZYpJVQEAgEgx5Fd2JGnZsmVasGCBJk2apOuuu06bNm2S2+3W4sWLzS4NAACYzBJh57vf/a4++eQTPfroo2ppaVFubq6ee+45XX755WaXZrqYmBitXr2619d2wwXzH97zl/gdDPf5S/wOhvv8JclmGF91vxYAAMDQNeSv2QEAADgfwg4AALA0wg4AALA0wg4AALA0ws4w89vf/lYFBQWKi4uT3W7X3LlzzS5p0Pn9fk2cOFE2m03vvvuu2eUMmiNHjmjhwoXKyspSXFycrrzySq1evVpdXV1mlzZgnnjiCWVlZSk2Nlb5+fl67bXXzC5p0FRUVGjy5MlKSEjQmDFjdOutt+rQoUNml2WaiooK2Ww2lZaWml3KoDp+/Ljmz5+vlJQUxcfHa+LEiWpoaDC7rEFH2BlGfvGLX2jBggW655579Mc//lG///3vNW/ePLPLGnQrV678ykeLW9HBgwd15swZ/dd//Zeampq0fv16/fSnP9UPfvADs0sbEM8884xKS0u1atUqvfPOO7rxxhs1a9Ysud1us0sbFHv37tXSpUu1b98+1dXV6fTp0yoqKtJnn31mdmmDbv/+/dq0aZOuvfZas0sZVO3t7br++usVHR2t559/Xu+//74ef/xxy78x4JzC8jZORLzu7m7jG9/4hvHf//3fZpdiqueee87IyckxmpqaDEnGO++8Y3ZJpqqsrDSysrLMLmNA/M3f/I2xePHioLacnBzj+9//vkkVmautrc2QZOzdu9fsUgZVZ2enkZ2dbdTV1RlTp041HnroIbNLGjSPPPKIccMNN5hdRkRgZWeYePvtt3X8+HFddNFFysvLU2pqqmbNmqWmpiazSxs0H3/8se6//35t27ZN8fHxZpcTETo6OpScnGx2GWHX1dWlhoYGFRUVBbUXFRWpvr7epKrM1dHRIUmW/N/7fJYuXaqbb75ZhYWFZpcy6Hbt2qVJkybp9ttv15gxY5SXl6cnn3zS7LJMQdgZJj788ENJksvl0g9/+EP95je/0ejRozV16lR9+umnJlc38AzD0N13363Fixdr0qRJZpcTEf7nf/5H1dXVlnytisfjUU9PT6+XATscjl4vDR4ODMPQsmXLdMMNNyg3N9fscgZNbW2t3n77bVVUVJhdiik+/PBDbdy4UdnZ2XrxxRe1ePFiPfjgg3r66afNLm3QEXaGOJfLJZvNdt7trbfe0pkzZyRJq1at0m233ab8/Hxt2bJFNptNP//5z02eRei+7vyrq6vl9XpVVlZmdslh93V/B//XiRMn9O1vf1u333677rvvPpMqH3g2my1o3zCMXm3DwQMPPKD33ntPP/vZz8wuZdAcPXpUDz30kLZv367Y2FizyzHFmTNn9K1vfUvl5eXKy8vTokWLdP/992vjxo1mlzboLPFurOHsgQce0B133HHePpmZmers7JQkjR8/PtAeExOjK664YkhfsPl157927Vrt27ev17thJk2apDvvvFNbt24dyDIH1Nf9HZx14sQJ3XTTTYGX5lqR3W7XiBEjeq3itLW19VrtsbqSkhLt2rVLr776qi677DKzyxk0DQ0NamtrU35+fqCtp6dHr776qjZs2CC/368RI0aYWOHAS01NDfo3X5LGjRunX/ziFyZVZB7CzhBnt9tlt9u/sl9+fr5iYmJ06NAh3XDDDZKk7u5uHTlyZEi/MPXrzv8///M/tXbt2sD+iRMnNHPmTD3zzDMqKCgYyBIH3Nf9HUhf3IZ60003BVb2LrrImou7I0eOVH5+vurq6vSP//iPgfa6ujp95zvfMbGywWMYhkpKSrRz50698sorysrKMrukQTV9+nQdOHAgqO2ee+5RTk6OHnnkEcsHHUm6/vrrez1u4E9/+tOQ/jc/VISdYSIxMVGLFy/W6tWrlZ6erssvv1w//vGPJUm33367ydUNvIyMjKD9iy++WJJ05ZVXDptPuydOnNC0adOUkZGhf//3f9df/vKXwDGn02liZQNj2bJlWrBggSZNmhRYxXK73Za8Rulcli5dqh07dujXv/61EhISAqtcSUlJiouLM7m6gZeQkNDr+qRRo0YpJSVl2Fy39PDDD2vKlCkqLy9XcXGx/vCHP2jTpk2WXdE9H8LOMPLjH/9YUVFRWrBggXw+nwoKCvTyyy9r9OjRZpeGQbB792598MEH+uCDD3oFPMMwTKpq4Hz3u9/VJ598okcffVQtLS3Kzc3Vc889N2w+1Z69LmPatGlB7Vu2bNHdd989+AVh0E2ePFk7d+5UWVmZHn30UWVlZamqqkp33nmn2aUNOpthxX/lAAAA/pc1v7AHAAD4X4QdAABgaYQdAABgaYQdAABgaYQdAABgaYQdAABgaYQdAABgaYQdAABgaYQdAABgaYQdAABgaYQdAABgaYQdAABgaf8PmHOcCrzuNVwAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"proportion >6.6 = 0.06%\n"
]
}
],
"source": [
"nReps = 10000 # (number of shuffles)\n",
"mDiff = np.empty(nReps) # array to store mean difference for each shuffle\n",
"\n",
"for j in range(nReps):\n",
" socks_shuffled = socks.copy() # work on a copy of the original dataframe\n",
" for i in range(len(socks)):\n",
" if np.random.rand()>0.5: # generate a random number between 0 and 1 - if it is more than 0.5:\n",
" socks_shuffled.loc[i,'Husband'] = socks.loc[i,'Wife'] # flip number of socks for husband and wife\n",
" socks_shuffled.loc[i,'Wife'] = socks.loc[i,'Husband'] # flip number of socks for husband and wife\n",
" #else:\n",
" # don't shuffle the row!\n",
" mDiff[j] = np.mean(socks_shuffled.Husband - socks_shuffled.Wife)\n",
" \n",
"sns.histplot(mDiff)\n",
"plt.show()\n",
"\n",
"print('proportion >6.6 = ' + str(100*np.mean(mDiff>6.6)) + '%')"
]
},
{
"cell_type": "markdown",
"id": "0e88511f",
"metadata": {},
"source": [
"We can see that the null distribution for the mean difference in socks owned between husbands and wives is a bit bimodal. This is probably due to the large influence of outliers (points to the right of the plot above are probably those where the two high-sock individuals were assigned to be husbands; points to the left are probably cases where the two high-sock individuals were assiged to be wives)\n",
"\n",
"### The $p$ value\n",
"\n",
"We can also calculate the proportion of cases in which the mean difference in socks owned for [Husband-Wife] exceeds the value we observed in our original sample, 6.6. This proportion is about 0.06% (it will actually vary on each run of the permutation test as the permutations are random - but hopefully not much). It tells us that if we simulate a situation in which sex does not determine the number of socks owned (but preserving some other important features of the dataset like the high skew, and the correlation between husabnds and their wives), there is only a 0.06% chance that we would get an apparent sex difference as large as the one we observed in our 'real' data.\n",
"\n",
"The probability that the test statistic (in this case, the mean difference in pairs of socks owned) would be observed if the null hypothesis were true, is sometimes called the $p$-value. \n",
"\n",
"Our permutation test shows that the $p$-value associated with the observed difference of means is 0.0006.\n",
"\n",
"The result is considered statistically significant if $p$ is smaller than some predetermined level, known as $\\alpha$. Usually $\\alpha = 0.05$ or $\\alpha = 0.01$ is used, so the result is significant if $p<0.05$ or $p<0.01$. Our result is therefore statistically significant."
]
},
{
"cell_type": "markdown",
"id": "ae1e1fe3",
"metadata": {},
"source": [
"## Use a built in function\n",
"\n",
"Now you have seen how the permutation test works, we can learn how to run it more easily using the built in function scipy.stats.permutation_test\n",
"\n",
"Note- For those NOT using colab - You need scipy stats version > 1.8.0 to run this. You may need to check your version by running the following code block."
]
},
{
"cell_type": "code",
"execution_count": 102,
"id": "58c0b46d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'1.10.0'"
]
},
"execution_count": 102,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import scipy as scipy\n",
"scipy.version.version"
]
},
{
"cell_type": "markdown",
"id": "b4a2e8ea",
"metadata": {},
"source": [
"If this is less than 1.8.0 you need to update it - see the technical note in the first page of this chapter\n",
"\n",
"For those who are using Colab - check you followed the instructions at the top of this page"
]
},
{
"cell_type": "markdown",
"id": "4876cd34",
"metadata": {},
"source": [
"### Syntax of stats.permutation_test\n",
"\n",
"Here is how we run the permutation test (same as the one we did with our own code above, although note how much more quickly this one runs!)"
]
},
{
"cell_type": "code",
"execution_count": 88,
"id": "200125b4",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"PermutationTestResult(statistic=6.6, pvalue=0.005859375, null_distribution=array([ 6.6, 5. , 5.6, ..., -5.6, -5. , -6.6]))"
]
},
"execution_count": 88,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def mDiff(x, y):\n",
" return np.mean(x-y)\n",
"\n",
"stats.permutation_test((socks.Husband, socks.Wife), mDiff, permutation_type='samples', alternative='two-sided', n_resamples=10000)"
]
},
{
"cell_type": "markdown",
"id": "a13b89d1",
"metadata": {},
"source": [
"Firstly, to reassure you this is doing a very similar job to our home-made code, check the p-value (should be about 0.06). \n",
"\n",
"We can also plot the null distribution, which hopefully looks similar to what we got from the home-made code:"
]
},
{
"cell_type": "code",
"execution_count": 89,
"id": "ea965aee",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGdCAYAAAD0e7I1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAl3klEQVR4nO3df1TVdZ7H8df11xUKSCDvhQkCGyZ1qHTB4Sw5K65K49TMtM7ktKjb76PrjyJrNZYab54JNpuIszDa2Ezm5CE7e2bccWdtkmrSZqhZxCxl0MYNg1UYBmMBg7n8+u4frXe7i2jihe+XD8/HOd9z+v64l/f3lvjse79wXZZlWQIAADDUGLsHAAAAGErEDgAAMBqxAwAAjEbsAAAAoxE7AADAaMQOAAAwGrEDAACMRuwAAACjjbN7ACfo6+vTqVOnFBERIZfLZfc4AADgc7AsS+3t7YqPj9eYMQNfvyF2JJ06dUoJCQl2jwEAAAahvr5eV1111YD7iR1JERERkj59sSIjI22eBgAAfB5tbW1KSEgI/D0+EGJHCrx1FRkZSewAADDCXOgWFG5QBgAARiN2AACA0YgdAABgNGIHAAAYjdgBAABGI3YAAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARrP1U8/379+vp556SlVVVWpoaNCuXbt06623SpK6u7v16KOPas+ePfrwww8VFRWl+fPn65/+6Z8UHx8feA6/36+HH35YL730kjo7OzVv3jxt3rxZV111lU1nhZGmrq5Ozc3Ndo+h2NhYJSYm2j0GcFH484ORwNbY+eSTT3TDDTforrvu0re//e2gfR0dHTp48KAee+wx3XDDDWppaVFubq6++c1v6sCBA4HjcnNz9W//9m/auXOnYmJi9NBDD+mWW25RVVWVxo4dO9ynhBGmrq5OU6dOU2dnh92jKCwsXEeP1vANGyMGf34wUtgaOwsXLtTChQvPuS8qKkrl5eVB20pKSvSVr3xFdXV1SkxMVGtrq37yk5/oxRdf1Pz58yVJO3bsUEJCgl577TXddNNNQ34OGNmam5vV2dmhjLs3KDIuybY52hpO6HfPP67m5ma+WWPE4M8PRgpbY+ditba2yuVy6YorrpAkVVVVqbu7W9nZ2YFj4uPjlZqaqoqKigFjx+/3y+/3B9bb2tqGdG44X2RckqITr7V7DGBE4s8PnG7E3KD85z//WY888ohycnIUGRkpSWpsbNSECRM0adKkoGM9Ho8aGxsHfK7CwkJFRUUFloSEhCGdHQAA2GdExE53d7duv/129fX1afPmzRc83rIsuVyuAffn5eWptbU1sNTX14dyXAAA4CCOj53u7m4tXrxYtbW1Ki8vD1zVkSSv16uuri61tLQEPaapqUkej2fA53S73YqMjAxaAACAmRwdO2dD5w9/+INee+01xcTEBO1PS0vT+PHjg25kbmho0JEjR5SZmTnc4wIAAAey9QblM2fO6Pjx44H12tpaHTp0SNHR0YqPj9d3vvMdHTx4UL/85S/V29sbuA8nOjpaEyZMUFRUlO655x499NBDiomJUXR0tB5++GFdd911gZ/OAgAAo5utsXPgwAHNnTs3sL527VpJ0h133CGfz6fdu3dLkmbMmBH0uF//+tfKysqSJD3zzDMaN26cFi9eHPilgi+88AK/YwcAAEiyOXaysrJkWdaA+8+376yJEyeqpKREJSUloRwNAAAYwtH37AAAAFyqEfVLBQHT1dTU2D2CJD5nyOmc8nlUTvnvFbgQYgdwgM7W05JcWrp0qd2jSOJzhpzMSZ9HdVa3v8vuEYDzInYAB+juaJdkaUbOel2ZPNXWWficIWdzyudRSVLD4bd1ZPdW9fT02DoHcCHEDuAgl09O5DOG8Lk44fOo2hpO2Pr1gc+LG5QBAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0YgdAABgNGIHAAAYjdgBAABGI3YAAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0YgdAABgNGIHAAAYjdgBAABGI3YAAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0YgdAABgNGIHAAAYjdgBAABGI3YAAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRbY2f//v36xje+ofj4eLlcLv3rv/5r0H7LsuTz+RQfH6+wsDBlZWWpuro66Bi/3681a9YoNjZWl112mb75zW/qv/7rv4bxLAAAgJPZGjuffPKJbrjhBpWWlp5z/6ZNm1RUVKTS0lJVVlbK6/VqwYIFam9vDxyTm5urXbt2aefOnfrNb36jM2fO6JZbblFvb+9wnQYAAHCwcXZ+8YULF2rhwoXn3GdZloqLi5Wfn69FixZJkrZv3y6Px6OysjItX75cra2t+slPfqIXX3xR8+fPlyTt2LFDCQkJeu2113TTTTcN27kAAABncuw9O7W1tWpsbFR2dnZgm9vt1pw5c1RRUSFJqqqqUnd3d9Ax8fHxSk1NDRxzLn6/X21tbUELAAAwk2Njp7GxUZLk8XiCtns8nsC+xsZGTZgwQZMmTRrwmHMpLCxUVFRUYElISAjx9AAAwCkcGztnuVyuoHXLsvpt+/8udExeXp5aW1sDS319fUhmBQAAzmPrPTvn4/V6JX169SYuLi6wvampKXC1x+v1qqurSy0tLUFXd5qampSZmTngc7vdbrnd7iGaHDBDTU2N3SNIkmJjY5WYmGj3GKqrq1Nzc7PdYzjm3wswkjg2dpKTk+X1elVeXq6ZM2dKkrq6urRv3z49+eSTkqS0tDSNHz9e5eXlWrx4sSSpoaFBR44c0aZNm2ybHRjJOltPS3Jp6dKldo8iSQoLC9fRozW2Bk9dXZ2mTp2mzs4O22b4/7r9XXaPAIwYtsbOmTNndPz48cB6bW2tDh06pOjoaCUmJio3N1cFBQVKSUlRSkqKCgoKFB4erpycHElSVFSU7rnnHj300EOKiYlRdHS0Hn74YV133XWBn84CcHG6O9olWZqRs15XJk+1dZa2hhP63fOPq7m52dbYaW5uVmdnhzLu3qDIuCTb5pCkhsNv68jurerp6bF1DmAksTV2Dhw4oLlz5wbW165dK0m644479MILL2jdunXq7OzUypUr1dLSooyMDO3du1cRERGBxzzzzDMaN26cFi9erM7OTs2bN08vvPCCxo4dO+zng8+PtwSc7/LJiYpOvNbuMRwlMi7J9tekreGErV/fyZzy59kpb73i/9gaO1lZWbIsa8D9LpdLPp9PPp9vwGMmTpyokpISlZSUDMGEGAq8JQAglHjrFRfi2Ht2YC7eEgAQSrz1igshdmAb3hIAEEq89YqBOP737AAAAFwKYgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0YgdAABgNGIHAAAYjdgBAABGI3YAAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0YgdAABgNGIHAAAYjdgBAABGI3YAAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0YgdAABgNGIHAAAYjdgBAABGI3YAAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0YgdAABgNEfHTk9Pjx599FElJycrLCxMU6ZM0caNG9XX1xc4xrIs+Xw+xcfHKywsTFlZWaqurrZxagAA4CSOjp0nn3xSzz77rEpLS1VTU6NNmzbpqaeeUklJSeCYTZs2qaioSKWlpaqsrJTX69WCBQvU3t5u4+QAAMApHB07b7/9tr71rW/p5ptvVlJSkr7zne8oOztbBw4ckPTpVZ3i4mLl5+dr0aJFSk1N1fbt29XR0aGysjKbpwcAAE7g6NiZPXu2Xn/9dX3wwQeSpPfee0+/+c1v9PWvf12SVFtbq8bGRmVnZwce43a7NWfOHFVUVNgyMwAAcJZxdg9wPuvXr1dra6umTp2qsWPHqre3V0888YT+9m//VpLU2NgoSfJ4PEGP83g8+uijjwZ8Xr/fL7/fH1hva2sbgukBAIATOPrKzssvv6wdO3aorKxMBw8e1Pbt2/WDH/xA27dvDzrO5XIFrVuW1W/bZxUWFioqKiqwJCQkDMn8AADAfo6OnX/4h3/QI488ottvv13XXXedli1bpgcffFCFhYWSJK/XK+n/rvCc1dTU1O9qz2fl5eWptbU1sNTX1w/dSQAAAFs5OnY6Ojo0ZkzwiGPHjg386HlycrK8Xq/Ky8sD+7u6urRv3z5lZmYO+Lxut1uRkZFBCwAAMJOj79n5xje+oSeeeEKJiYn68pe/rHfffVdFRUW6++67JX369lVubq4KCgqUkpKilJQUFRQUKDw8XDk5OTZPDwAAnMDRsVNSUqLHHntMK1euVFNTk+Lj47V8+XJ973vfCxyzbt06dXZ2auXKlWppaVFGRob27t2riIgIGycHAABO4ejYiYiIUHFxsYqLiwc8xuVyyefzyefzDdtcAABg5HD0PTsAAACXitgBAABGI3YAAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0YgdAABgNGIHAAAYjdgBAABGI3YAAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0QYVO1OmTNHp06f7bf/v//5vTZky5ZKHAgAACJVBxc6JEyfU29vbb7vf79fJkycveSgAAIBQGXcxB+/evTvwz6+++qqioqIC6729vXr99deVlJQUsuEAAAAu1UXFzq233ipJcrlcuuOOO4L2jR8/XklJSXr66adDNhwAAMCluqjY6evrkyQlJyersrJSsbGxQzIUAABAqFxU7JxVW1sb6jkAAACGxKBiR5Jef/11vf7662pqagpc8Tnr+eefv+TBAAAAQmFQsfP4449r48aNSk9PV1xcnFwuV6jnAgAACIlBxc6zzz6rF154QcuWLQv1PAAAACE1qN+z09XVpczMzFDPAgAAEHKDip17771XZWVloZ4FAAAg5Ab1Ntaf//xnbd26Va+99pquv/56jR8/Pmh/UVFRSIYDAAC4VIOKnffff18zZsyQJB05ciRoHzcrAwAAJxlU7Pz6178O9RwAAABDYlD37AAAAIwUg7qyM3fu3PO+XfXGG28MeiAAAIBQGlTsnL1f56zu7m4dOnRIR44c6fcBoQAAAHYaVOw888wz59zu8/l05syZSxoIAAAglEJ6z87SpUv5XCwAAOAoIY2dt99+WxMnTgzlUwIAAFySQb2NtWjRoqB1y7LU0NCgAwcO6LHHHgvJYAAAAKEwqNiJiooKWh8zZoyuvfZabdy4UdnZ2SEZDAAAIBQGFTvbtm0L9RwAAABDYlCxc1ZVVZVqamrkcrk0ffp0zZw5M1RzAQAAhMSgYqepqUm333673nzzTV1xxRWyLEutra2aO3eudu7cqSuvvDLUcwIAAAzKoH4aa82aNWpra1N1dbU+/vhjtbS06MiRI2pra9P9998f0gFPnjyppUuXKiYmRuHh4ZoxY4aqqqoC+y3Lks/nU3x8vMLCwpSVlaXq6uqQzgAAAEauQcXOr371K23ZskXTpk0LbJs+fbp++MMf6pVXXgnZcC0tLbrxxhs1fvx4vfLKK/r973+vp59+WldccUXgmE2bNqmoqEilpaWqrKyU1+vVggUL1N7eHrI5AADAyDWot7H6+vo0fvz4ftvHjx+vvr6+Sx7qrCeffFIJCQlBN0QnJSUF/tmyLBUXFys/Pz/w4/Dbt2+Xx+NRWVmZli9fHrJZAADAyDSoKzt//dd/rQceeECnTp0KbDt58qQefPBBzZs3L2TD7d69W+np6brttts0efJkzZw5U88991xgf21trRobG4N+3N3tdmvOnDmqqKgY8Hn9fr/a2tqCFgAAYKZBxU5paana29uVlJSka665Rl/84heVnJys9vZ2lZSUhGy4Dz/8UFu2bFFKSopeffVVrVixQvfff79++tOfSpIaGxslSR6PJ+hxHo8nsO9cCgsLFRUVFVgSEhJCNjMAAHCWQb2NlZCQoIMHD6q8vFxHjx6VZVmaPn265s+fH9Lh+vr6lJ6eroKCAknSzJkzVV1drS1btujv/u7vAse5XK6gx1mW1W/bZ+Xl5Wnt2rWB9ba2NoIHAABDXdSVnTfeeEPTp08PvO2zYMECrVmzRvfff79mzZqlL3/5y3rrrbdCNlxcXJymT58etG3atGmqq6uTJHm9XknqdxWnqamp39Wez3K73YqMjAxaAACAmS4qdoqLi3XfffedMw6ioqK0fPlyFRUVhWy4G2+8UceOHQva9sEHH+jqq6+WJCUnJ8vr9aq8vDywv6urS/v27VNmZmbI5gAAACPXRcXOe++9p6997WsD7s/Ozg76HTiX6sEHH9Q777yjgoICHT9+XGVlZdq6datWrVol6dO3r3Jzc1VQUKBdu3bpyJEjuvPOOxUeHq6cnJyQzQEAAEaui7pn549//OM5f+Q88GTjxulPf/rTJQ911qxZs7Rr1y7l5eVp48aNSk5OVnFxsZYsWRI4Zt26ders7NTKlSvV0tKijIwM7d27VxERESGbA4C9ampqRvXXB3BpLip2vvCFL+jw4cP64he/eM7977//vuLi4kIy2Fm33HKLbrnllgH3u1wu+Xw++Xy+kH5dAPbrbD0tyaWlS5faPYokqdvfZfcIAAbhomLn61//ur73ve9p4cKFmjhxYtC+zs5Obdiw4bxhAgAXo7ujXZKlGTnrdWXyVNvmaDj8to7s3qqenh7bZgAweBcVO48++qh+/vOf60tf+pJWr16ta6+9Vi6XSzU1NfrhD3+o3t5e5efnD9WsAEapyycnKjrxWtu+flvDCdu+NoBLd1Gx4/F4VFFRob//+79XXl6eLMuS9OlbSTfddJM2b9583h/5BgAAGG4X/UsFr776au3Zs0ctLS06fvy4LMtSSkqKJk2aNBTzAQAAXJJB/QZlSZo0aZJmzZoVylkAAABCblCfjQUAADBSEDsAAMBoxA4AADAasQMAAIxG7AAAAKMROwAAwGjEDgAAMBqxAwAAjEbsAAAAoxE7AADAaMQOAAAwGrEDAACMRuwAAACjETsAAMBoxA4AADAasQMAAIxG7AAAAKMROwAAwGjEDgAAMBqxAwAAjEbsAAAAoxE7AADAaMQOAAAwGrEDAACMRuwAAACjETsAAMBoxA4AADAasQMAAIxG7AAAAKMROwAAwGjEDgAAMBqxAwAAjEbsAAAAoxE7AADAaMQOAAAwGrEDAACMRuwAAACjETsAAMBoxA4AADAasQMAAIxG7AAAAKMROwAAwGgjKnYKCwvlcrmUm5sb2GZZlnw+n+Lj4xUWFqasrCxVV1fbNyQAAHCUERM7lZWV2rp1q66//vqg7Zs2bVJRUZFKS0tVWVkpr9erBQsWqL293aZJAQCAk4yI2Dlz5oyWLFmi5557TpMmTQpstyxLxcXFys/P16JFi5Samqrt27ero6NDZWVlNk4MAACcYkTEzqpVq3TzzTdr/vz5Qdtra2vV2Nio7OzswDa32605c+aooqJiwOfz+/1qa2sLWgAAgJnG2T3AhezcuVMHDx5UZWVlv32NjY2SJI/HE7Td4/Hoo48+GvA5CwsL9fjjj4d2UAAA4EiOvrJTX1+vBx54QDt27NDEiRMHPM7lcgWtW5bVb9tn5eXlqbW1NbDU19eHbGYAAOAsjr6yU1VVpaamJqWlpQW29fb2av/+/SotLdWxY8ckfXqFJy4uLnBMU1NTv6s9n+V2u+V2u4ducAAA4BiOvrIzb948HT58WIcOHQos6enpWrJkiQ4dOqQpU6bI6/WqvLw88Jiuri7t27dPmZmZNk4OAACcwtFXdiIiIpSamhq07bLLLlNMTExge25urgoKCpSSkqKUlBQVFBQoPDxcOTk5dowMAAAcxtGx83msW7dOnZ2dWrlypVpaWpSRkaG9e/cqIiLC7tEAAIADjLjYefPNN4PWXS6XfD6ffD6fLfMAAABnc/Q9OwAAAJdqxF3ZweDV1dWpubnZ7jFUU1Nj9wgAMKSc8n0uNjZWiYmJdo9hO2JnlKirq9PUqdPU2dlh9ygB3f4uu0cAgJDqbD0tyaWlS5faPYokKSwsXEeP1oz64CF2Ronm5mZ1dnYo4+4NioxLsnWWhsNv68jurerp6bF1DgAIte6OdkmWZuSs15XJU22dpa3hhH73/ONqbm4mduweAMMrMi5J0YnX2jpDW8MJW78+AAy1yycn2v69Fv+HG5QBAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0YgdAABgNGIHAAAYjdgBAABGI3YAAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0YgdAABgNGIHAAAYjdgBAABGI3YAAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0YgdAABgNGIHAAAYjdgBAABGI3YAAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRHx05hYaFmzZqliIgITZ48WbfeequOHTsWdIxlWfL5fIqPj1dYWJiysrJUXV1t08QAAMBpHB07+/bt06pVq/TOO++ovLxcPT09ys7O1ieffBI4ZtOmTSoqKlJpaakqKyvl9Xq1YMECtbe32zg5AABwinF2D3A+v/rVr4LWt23bpsmTJ6uqqkp/9Vd/JcuyVFxcrPz8fC1atEiStH37dnk8HpWVlWn58uV2jA0AABzE0Vd2/r/W1lZJUnR0tCSptrZWjY2Nys7ODhzjdrs1Z84cVVRUDPg8fr9fbW1tQQsAADDTiIkdy7K0du1azZ49W6mpqZKkxsZGSZLH4wk61uPxBPadS2FhoaKiogJLQkLC0A0OAABsNWJiZ/Xq1Xr//ff10ksv9dvncrmC1i3L6rfts/Ly8tTa2hpY6uvrQz4vAABwBkffs3PWmjVrtHv3bu3fv19XXXVVYLvX65X06RWeuLi4wPampqZ+V3s+y+12y+12D93AAADAMRx9ZceyLK1evVo///nP9cYbbyg5OTlof3Jysrxer8rLywPburq6tG/fPmVmZg73uAAAwIEcfWVn1apVKisr0y9+8QtFREQE7sOJiopSWFiYXC6XcnNzVVBQoJSUFKWkpKigoEDh4eHKycmxeXoAAOAEjo6dLVu2SJKysrKCtm/btk133nmnJGndunXq7OzUypUr1dLSooyMDO3du1cRERHDPC0AAHAiR8eOZVkXPMblcsnn88nn8w39QAAAYMRx9D07AAAAl4rYAQAARiN2AACA0YgdAABgNGIHAAAYjdgBAABGI3YAAIDRiB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0YgdAABgNGIHAAAYjdgBAABGG2f3AKarq6tTc3Oz3WOopqbG7hEAADZwwvf/2NhYJSYm2vb1iZ0hVFdXp6lTp6mzs8PuUQK6/V12jwAAGAadraclubR06VK7R1FYWLiOHq2xLXiInSHU3Nyszs4OZdy9QZFxSbbO0nD4bR3ZvVU9PT22zgEAGB7dHe2SLM3IWa8rk6faNkdbwwn97vnH1dzcTOyYLDIuSdGJ19o6Q1vDCVu/PgDAHpdPTrT97yC7cYMyAAAwGrEDAACMRuwAAACjETsAAMBoxA4AADAasQMAAIxG7AAAAKMROwAAwGjEDgAAMBqxAwAAjEbsAAAAoxE7AADAaMQOAAAwGrEDAACMRuwAAACjETsAAMBoxA4AADAasQMAAIxG7AAAAKMROwAAwGjEDgAAMBqxAwAAjEbsAAAAoxE7AADAaMQOAAAwmjGxs3nzZiUnJ2vixIlKS0vTW2+9ZfdIAADAAYyInZdfflm5ubnKz8/Xu+++q69+9atauHCh6urq7B4NAADYzIjYKSoq0j333KN7771X06ZNU3FxsRISErRlyxa7RwMAADYbZ/cAl6qrq0tVVVV65JFHgrZnZ2eroqLinI/x+/3y+/2B9dbWVklSW1tbSGc7c+aMJOnjj46px98Z0ue+WG0NH0mSWk/+QePHuZjFYbM4ZQ5mcfYczOLsOZhlgDkaP32X5cyZMyH/e/bs81mWdf4DrRHu5MmTliTrt7/9bdD2J554wvrSl750zsds2LDBksTCwsLCwsJiwFJfX3/eVhjxV3bOcrmCq9WyrH7bzsrLy9PatWsD6319ffr4448VExMz4GNGqra2NiUkJKi+vl6RkZF2jzPsOP/Rff4Sr8FoP3+J18Dk87csS+3t7YqPjz/vcSM+dmJjYzV27Fg1NjYGbW9qapLH4znnY9xut9xud9C2K664YqhGdITIyEjj/iO/GJz/6D5/iddgtJ+/xGtg6vlHRUVd8JgRf4PyhAkTlJaWpvLy8qDt5eXlyszMtGkqAADgFCP+yo4krV27VsuWLVN6err+8i//Ulu3blVdXZ1WrFhh92gAAMBmRsTOd7/7XZ0+fVobN25UQ0ODUlNTtWfPHl199dV2j2Y7t9utDRs29HvbbrTg/Ef3+Uu8BqP9/CVeg9F+/pLksqwL/bwWAADAyDXi79kBAAA4H2IHAAAYjdgBAABGI3YAAIDRiJ1R5t///d+VkZGhsLAwxcbGatGiRXaPNOz8fr9mzJghl8ulQ4cO2T3OsDlx4oTuueceJScnKywsTNdcc402bNigrq4uu0cbMps3b1ZycrImTpyotLQ0vfXWW3aPNGwKCws1a9YsRUREaPLkybr11lt17Ngxu8eyTWFhoVwul3Jzc+0eZVidPHlSS5cuVUxMjMLDwzVjxgxVVVXZPdawI3ZGkZ/97GdatmyZ7rrrLr333nv67W9/q5ycHLvHGnbr1q274K8WN9HRo0fV19enH/3oR6qurtYzzzyjZ599Vv/4j/9o92hD4uWXX1Zubq7y8/P17rvv6qtf/aoWLlyouro6u0cbFvv27dOqVav0zjvvqLy8XD09PcrOztYnn3xi92jDrrKyUlu3btX1119v9yjDqqWlRTfeeKPGjx+vV155Rb///e/19NNPG/+JAecUkk/jhON1d3dbX/jCF6wf//jHdo9iqz179lhTp061qqurLUnWu+++a/dIttq0aZOVnJxs9xhD4itf+Yq1YsWKoG1Tp061HnnkEZsmsldTU5Mlydq3b5/dowyr9vZ2KyUlxSovL7fmzJljPfDAA3aPNGzWr19vzZ492+4xHIErO6PEwYMHdfLkSY0ZM0YzZ85UXFycFi5cqOrqartHGzZ//OMfdd999+nFF19UeHi43eM4Qmtrq6Kjo+0eI+S6urpUVVWl7OzsoO3Z2dmqqKiwaSp7tba2SpKR/77PZ9WqVbr55ps1f/58u0cZdrt371Z6erpuu+02TZ48WTNnztRzzz1n91i2IHZGiQ8//FCS5PP59Oijj+qXv/ylJk2apDlz5ujjjz+2ebqhZ1mW7rzzTq1YsULp6el2j+MI//mf/6mSkhIjP1alublZvb29/T4M2OPx9PvQ4NHAsiytXbtWs2fPVmpqqt3jDJudO3fq4MGDKiwstHsUW3z44YfasmWLUlJS9Oqrr2rFihW6//779dOf/tTu0YYdsTPC+Xw+uVyu8y4HDhxQX1+fJCk/P1/f/va3lZaWpm3btsnlculf/uVfbD6Lwfu8519SUqK2tjbl5eXZPXLIfd7X4LNOnTqlr33ta7rtttt077332jT50HO5XEHrlmX12zYarF69Wu+//75eeuklu0cZNvX19XrggQe0Y8cOTZw40e5xbNHX16e/+Iu/UEFBgWbOnKnly5frvvvu05YtW+webdgZ8dlYo9nq1at1++23n/eYpKQktbe3S5KmT58e2O52uzVlypQRfcPm5z3/73//+3rnnXf6fTZMenq6lixZou3btw/lmEPq874GZ506dUpz584NfGiuiWJjYzV27Nh+V3Gampr6Xe0x3Zo1a7R7927t379fV111ld3jDJuqqio1NTUpLS0tsK23t1f79+9XaWmp/H6/xo4da+OEQy8uLi7oe74kTZs2TT/72c9smsg+xM4IFxsbq9jY2Asel5aWJrfbrWPHjmn27NmSpO7ubp04cWJEf2Dq5z3/f/7nf9b3v//9wPqpU6d000036eWXX1ZGRsZQjjjkPu9rIH36Y6hz584NXNkbM8bMi7sTJkxQWlqaysvL9Td/8zeB7eXl5frWt75l42TDx7IsrVmzRrt27dKbb76p5ORku0caVvPmzdPhw4eDtt11112aOnWq1q9fb3zoSNKNN97Y79cNfPDBByP6e/5gETujRGRkpFasWKENGzYoISFBV199tZ566ilJ0m233WbzdEMvMTExaP3yyy+XJF1zzTWj5v92T506paysLCUmJuoHP/iB/vSnPwX2eb1eGycbGmvXrtWyZcuUnp4euIpVV1dn5D1K57Jq1SqVlZXpF7/4hSIiIgJXuaKiohQWFmbzdEMvIiKi3/1Jl112mWJiYkbNfUsPPvigMjMzVVBQoMWLF+s//uM/tHXrVmOv6J4PsTOKPPXUUxo3bpyWLVumzs5OZWRk6I033tCkSZPsHg3DYO/evTp+/LiOHz/eL/Asy7JpqqHz3e9+V6dPn9bGjRvV0NCg1NRU7dmzZ9T8X+3Z+zKysrKCtm/btk133nnn8A+EYTdr1izt2rVLeXl52rhxo5KTk1VcXKwlS5bYPdqwc1kmfpcDAAD4X2a+YQ8AAPC/iB0AAGA0YgcAABiN2AEAAEYjdgAAgNGIHQAAYDRiBwAAGI3YAQAARiN2AACA0YgdAABgNGIHAAAYjdgBAABG+x873T5RmkNBTwAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"res = stats.permutation_test((socks.Husband, socks.Wife), mDiff, permutation_type='samples', alternative='two-sided', n_resamples=10000)\n",
"sns.histplot(res.null_distribution)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "1fdec81e",
"metadata": {},
"source": [
"However, the syntax may be a bit unfamiliar.\n",
"\n",
"Firstly, we had to give the function stats.permutation_test our two samples (socks.Husband, socks.Wife) as a pair of series (individual columns from the dataframe), rather than giving it the whole pandas dataframe as we do for many other stats functions.\n",
"\n",
"Secondly, to tell stats.permutation_test the test statistic we want to get the null distribution of, we had to pass it a function called mDiff, and this function had to have the property that it takes in two series (socks.Husband, socks.Wife) and returns a single number mean(socks.Husband, socks.Wife)"
]
},
{
"cell_type": "markdown",
"id": "a7afbd3b",
"metadata": {},
"source": [
"### Defining a function\n",
"\n",
"You will have come across this in datacamp but we haven't used it since. Don't be scared! It's unfamiliar but quite handy. On the other hand for a pairwise permutation test, the function I have given you for mDiff is always going to work, so if in doubt you can just copy it :-)\n",
"\n",
"A function is a little computer programme that takes in some information (in this case, it takes in two series, (socks.Husband, socks.Wife) and returns some value (in this case the mean difference mean(socks.Husband, socks.Wife)"
]
},
{
"cell_type": "code",
"execution_count": 92,
"id": "ab69c8ea",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6.6"
]
},
"execution_count": 92,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# define a function\n",
"def mDiff(x, y):\n",
" return np.mean(x-y)\n",
"\n",
"# run the function for some inputs\n",
"mDiff(socks.Husband, socks.Wife)"
]
},
{
"cell_type": "markdown",
"id": "4be2278c",
"metadata": {},
"source": [
"Here's another example:"
]
},
{
"cell_type": "code",
"execution_count": 98,
"id": "a59f97f8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.30000000000000004"
]
},
"execution_count": 98,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# definte a new function that divides one element of each pair by the other, and then adds up the result across pairs\n",
"def bananas(x,y):\n",
" return sum(x/y)\n",
"\n",
"cats = np.array([1,2,3]) # one input array - have given it an arbitrary name\n",
"dogs = np.array([10,20,30]) # another input array - have given it an arbitrary name\n",
"\n",
"bananas(cats,dogs)"
]
},
{
"cell_type": "markdown",
"id": "6bc1764a",
"metadata": {},
"source": [
"Now we can see how we could run stats.permutation_test on our function bananas and our data cats and dogs"
]
},
{
"cell_type": "code",
"execution_count": 99,
"id": "51006dfc",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"PermutationTestResult(statistic=0.30000000000000004, pvalue=0.25, null_distribution=array([ 0.3, 10.2, 10.2, 20.1, 10.2, 20.1, 20.1, 30. ]))"
]
},
"execution_count": 99,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"stats.permutation_test((cats, dogs), bananas, permutation_type='samples', alternative='two-sided', n_resamples=10000)"
]
},
{
"cell_type": "markdown",
"id": "416776bf",
"metadata": {},
"source": [
"## Recap\n",
"\n",
"To run a permutation test on paired data, we randomly flipped some of the pairs so that the husband's sock count was assigned to the wife and vice versa. We did NOT move people between couples, as we want to retain the characteristic of the original dataset that high-sock husbands tend to have high-sock wives\n",
"\n",
"For each shuffle we calculated the mean (pairwise) difference in the number of socks - husband-wife. \n",
"\n",
"Permutation testing in this way gives us a null distribution for the mean difference. Values of mean difference that occur rarely in the null distribution are considered statistically significant.\n",
" \n",
"To run the permutation test with scipy.stats we need the option `permutation_type='samples'`"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7513a6b4",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.9.13"
}
},
"nbformat": 4,
"nbformat_minor": 5
}