{ "cells": [ { "cell_type": "markdown", "id": "572fb94e", "metadata": {}, "source": [ "# Permutation test for paired data\n", "\n", "Turning to the permutation test itself, 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": 3, "id": "3286b373", "metadata": { "tags": [] }, "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 as pd\n", "import seaborn as sns\n", "sns.set_theme(style='white')\n", "import statsmodels.api as sm\n", "import statsmodels.formula.api as smf" ] }, { "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", "#### Question & design\n", "\n", "We suspect men own ore socks than women. \n", "\n", "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", "\n", "#### Hypotheses\n", "\n", "We can state our hypotheses as follows:\n", "\n", "$\\mathcal{H_o}:$ The mean sex difference quantity of socks owned is zero \n", "* on average a wife and her husband have the same number of pairs of socks\n", "\n", "$\\mathcal{H_a}:$ The mean sex difference [husband-wife] quantity of socks owned greater than zero\n", "* on average a husband has more socks than his wife\n", "\n", "#### Data\n", "\n", "We obtain sock-counts for the following informal sample of 14 couples:" ] }, { "cell_type": "code", "execution_count": 4, "id": "e0483d68", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
HusbandWife
01012
11713
24820
32825
42318
51614
61813
73426
82722
92214
101210
111317
122221
131516
\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\n", "10 12 10\n", "11 13 17\n", "12 22 21\n", "13 15 16" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "socks = pd.read_csv('https://raw.githubusercontent.com/jillxoreilly/StatsCourseBook_2024/main/data/socks.csv')\n", "socks" ] }, { "cell_type": "markdown", "id": "5906a5fc", "metadata": {}, "source": [ "Let's plot those data. \n", "\n", "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", "* You won't be required to reproduce this type plot for the course, however you can find instructions on how to make it in section 2.5 of this website\n", "\n", "\n", "Why do I prefer the plot on the right? \n", "\n", "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)\n", "\n", "*however* if there were 1000 couples in the sample rather than 10 this plot would be too crowded and hard to inspect\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "7588c35f", "metadata": { "jupyter": { "source_hidden": true }, "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxAAAAF/CAYAAADZxC9bAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAC/tklEQVR4nOzdd3xV9fnA8c+5eyY3exHITthDloCKgFZRVNxarbvuaqt1tgIOrHXUgf4cqLXaqlSts1pXFRRUQFFZYYWZvW9y973n98cNF1KirCQnIc/79eKl96w8N4Tc85zv832+iqqqKkIIIYQQQgixD3RaByCEEEIIIYToPSSBEEIIIYQQQuwzSSCEEEIIIYQQ+0wSCCGEEEIIIcQ+kwRCCCGEEEIIsc8kgRBCCCGEEELsM0kghBBCCCGEEPtMEgghhBBCCCHEPpMEQgghhBBCCLHPJIEQQgjR4+zYsYPi4uI9/vzzn/8EYM2aNZx33nmMGDGCyZMn8+yzz2ocsRBC9B0GrQPobqNHjyYQCJCSkqJ1KEIIoYmamhpMJhPLli3TOpSfVFpaitls5uOPP0ZRlNh2p9NJQ0MDF110EdOmTWPOnDmsWLGCOXPm4HK5OO200w7o68lngxCir9ufz4Y+l0D4/X7C4bDWYQghhGZCoRCqqmodxs9at24dubm5pKam7rHvhRdewGQyMXv2bAwGA/n5+WzZsoVnnnnmgBMI+WwQQvR1+/PZ0OcSiJ0fRp988onGkQghhDamTp2qdQh7VVpaSkFBQYf7li1bxpgxYzAYdn2EjR8/nqeeeoq6ujqSkpL2++vJZ4MQoq/bn88GmQMhhBCix1m3bh11dXWce+65TJgwgXPOOYdFixYBUFlZSXp6ervjdyYA5eXl3R6rEEL0NX1uBEIIIUTPFggE2Lx5M1arlZtuugmbzcbbb7/NZZddxvPPP4/P58NkMrU7x2w2A9FSJCGEEF1LEgghhBA9islkYunSpRgMhliiMGTIEDZu3Mizzz6LxWIhEAi0O2dn4mCz2bo9XiGE6GukhEkIIUSPY7PZ9hhlKCoqoqqqivT0dKqrq9vt2/k6LS2t22IUQoi+ShIIIYQQPcratWsZOXLkHq0EV65cSUFBAWPGjGH58uXtuiYtWbKE3NzcA5pALYQQYv9IAiGEEKJHKSoqorCwkDlz5rBs2TI2btzIvffey4oVK7jiiis47bTTaGlp4fbbb2fDhg288cYbvPDCC1x++eVahy6EEH2CzIEQQgjRo+h0Op588kkeeOABrr/+epqbmxk0aBDPP/88xcXFAMyfP5977rmHmTNnkpKSwk033cTMmTM1jlwIIfoGSSCEEEL0OImJicydO/cn9w8bNoxXX321GyMSQgixk+YlTDt27KC4uHiPP//85z8BWLNmDeeddx4jRoxg8uTJPPvssxpHLIQQQgghRN+l+QhEaWkpZrOZjz/+GEVRYtudTicNDQ1cdNFFTJs2jTlz5rBixQrmzJmDy+XitNNO0zBqIYQQQoSDXhSdgXDAjd7kRI2E0ButWoclhOhimicQ69atIzc3N7aK6O5eeOEFTCYTs2fPxmAwkJ+fz5YtW3jmmWckgRBCCCE0FAn52fHtC1T++Aphvxu92UnGsLPJGnUROoNZ6/CEEF1I8xKm0tJSCgoKOty3bNkyxowZg8GwK88ZP348ZWVl1NXVdVeIQgghhNhNOOhl+/Ln2bHsGcJ+d3Sb3832pc+w49vnCQe9GkcohOhKmicQ69ato66ujnPPPZcJEyZwzjnnsGjRIgAqKytJT09vd/zOkYry8vJuj1UIIYQQoOgMVP74Sof7Kn54BUWneYGDEKILafovPBAIsHnzZqxWKzfddBM2m423336byy67jOeffx6fz7fHSqRmc3RY1O/3axGyEEJoKhIMojMatQ5D9HHhgDs28rDHPr+bcKAFnTWhm6MSQnQXTRMIk8nE0qVLMRgMsURhyJAhbNy4kWeffRaLxUIgEGh3zs7EwWazdXu8QgihpfJ3/031x58y+K5ZGJ1OrcMRfZje5ERvdnaYROjNTvQmhwZRCSG6i+YlTDabbY9RhqKiIqqqqkhPT6e6urrdvp2v09LSui1GIYTQ2o5/vUXZM8/SWlZGzeeLtA5H9HFqJETGsLM73Jc+9GwCrTXdHJEQojtpmkCsXbuWkSNHsmzZsnbbV65cSUFBAWPGjGH58uWEw+HYviVLlpCbm0tSUlJ3hyuEEJrYtuA1Nv/1bwD0O+M0Mk44XuOIRF+nN1rJGnUR/cZcht4cHQ3Tm51kjb6MjGFnsfa93+Gu+F7jKIUQXUXTEqaioiIKCwuZM2cOs2bNIiEhgQULFrBixQpee+01kpOTmT9/PrfffjuXXnopP/zwAy+88AJz5szRMmwhhOgWqqqy9R+vsH3BawD0P/dsss86Q+OohIjSGcxkjryArMMuIRxoQW9yoIaDbP3qcYyWOKxJ+VqHKIToIpomEDqdjieffJIHHniA66+/nubmZgYNGsTzzz9PcXExAPPnz+eee+5h5syZpKSkcNNNNzFz5kwtwxZCiC6nqipb/vYSO954E4ABF5xPv1NP0TQmIf7XzkXjYhOm9UYGTLwOVY2gN1ho2vYNRlsitqSO27ULIXonzfusJSYmMnfu3J/cP2zYMF599dVujEgIIbQXcrtjcx1yL72IzBknahyREPtGp4/Oa2za9g1r3v0NBks8Q0//K2ZnhsaRCSE6i+YJhBBCiD0Z4+IYcvdsmlevJW3aFK3DEWK/2VNKsCYMwOIagNEm8xaFOJRIAiGEED2EGg7TWrYZR0G0dtyamYk1M1PjqIQ4MAZLHINOeQqDyYGiM1Bf9jnOjOEYLS6tQxNCHCTN27gKIYSIJg/rH5nHDzffRsPyb7UOR4hOYbS4UHQGata+S+m/b2DNO78hHPRqHZYQ4iDJCIQQQmgsEgqx7qGHqftyCYpeT9jn0zokITqVPXUwBkscjtRB6AxmrcMRQhwkSSCEEEJDkWCQ0vsfpP7rpSgGA8U33UDSuLFahyVEp7Il5jLszH9gcqShKAr1m/6La8AkdHqj1qEJIQ6AlDAJIYRGIoEAa+/9czR5MBopufUmSR7EIcvsTEdRFLYve5bS929k46dzUNWI1mEJIQ6AjEAIIYQGIoEAa+75E40rvkdnMjHw9ltwjRiudVhCdDl7chGKTo81IRdQtA5HCHEAJIEQQggNKAYDxoQEdBYLg/54G/FDBmsdkhDdIiHnCIaf8zpWVzaqqtK4dTGu/hO0DksIsR8kgRBCCA0oOh2F115Fv1NPwdY/W+twhOhWO5OHTZ/dQ/XqfzFgwnVkjvyV1mEJIfaRzIEQQohuEmppYesrC1DDYQAUvV6SB9FnKYqCJS4LFB16c5zW4Qgh9oOMQAghRDcINjezatadtG4qI+R2k3fZJVqHJITmsg67CFfOJOxJhUTCQVpr1uBMH6Z1WEKIvZARCCGE6GKBxkZW/mEWrZvKMMbHk3bMNK1DEqLHiCYPAdZ9cDOr/nUZDVsWax2SEGIvZARCCCG6kL+unlV3zMa7fQfGhASG3D0bW79+WoclRI+iKProAnOKDjUS0jocIcReSAIhhBBdxF9Ty8o/zsJXUYkpOZkhd8/GmpGhdVhC9DiKTk/BtDvx1G/AkTKQcNBLoKWyrdWrEKKnkRImIYToAmo4zKrZd+GrqMScmsrQuXdJ8iDEz9DpjThSBhLyu1nzztWs/NdleBs2ax2WEKIDkkAIIUQXUPR6ci++AFv/bIbOvQtLWqrWIQnROygKkZAfNRIi5HdrHY0QogNSwiSEEJ1IVVUUJbq6bsJho3CNGI6i12sclRC9h8HkYOCMeQQ8tdiTCgn6GlEjYUy2JK1DE0K0kREIIYToJK1btvL9b3+Pt7w8tk2SByH2n9GagD2pEH9LNaveuIw1b18joxFC9CCSQAghRCdo2bSJlbffQWtZGZuf/5vW4QhxSIiEfIR8TYR8jQS99VqHI4RoIyVMQghxkNzr1rNq9l2EW1txFBZSeN01WockxCHB6urPwJMfR2+0Y4nLxN9SjdGagE5v1Do0Ifo0GYEQQoiD0LxmLavumEO4tRXnwBIG33kHBodD67CEOGTYkwqxxGXibShj5esXsuGjP6BGwlqHJUSfJgmEEEIcoKYfV0ZHHrxe4oYMZvCsP2Cw2bQOS4hDkr+lmqCnHk/9RkL+Zq3DEaJPkxImIYQ4AKqqsvXlV4n4fLhGDKfktpvRm81ahyXEIcuVPY6SE/6CPaUEozWBQGsNRltyrOuZEKL7yAiEEEIcAEVRKLn1JjJPnsHA22+R5EGIbuDqfzhGawLuyh9Y8fIZ7Fj+nNYhCdEnSQIhhBD7wVdZGft/o9NJ7sUXojOZNIxIiL6npXo1Yb+bxi1fEgkHtQ5HiD5HSpiEEGIf1X7xJeseeoTcSy4k44TpWocjRJ+VMexsDOY4EvOORqc3EvQ1YrS4tA5LiD5DRiCEEGIfVH/2OaUPPowaDuNetwFVVbUOSYg+LaV4OnqjlYbNi/j2bzOo3/SZ1iEJ0WdIAiGEEHtR9fEnrH/4MYhESJ02hcLfXC0TN4XoIeo3/ZdI0EPtuve1DkWIPkNKmIQQ4mdUvP8Bm558BoD0439B3q8vRdHJsxcheoq8ybdhSy4ibfBpAIQCLRhMshaLEF1JPgWFEOInlL/zbix5yJhxAnmXXybJgxA9jKIzkDHsbHR6I9Vr3ua7l07BU7dB67CEOKTJJ6EQQvyEiD8AQNapp5B7yUVStiRED6ZGwlSt/hchbwM1pf/WOhwhDmlSwiSEED+h3+mn4igqJH7oEEkehOjhFJ2ekhMepmbN22SMOA+AcMiH3mDRODIhDj0yAiGEEG1UVaXyPx8S8nhj21zDhkryIEQvYbTEkznyfAC2L3uOH//5K4K+Rm2DEuIQJAmEEEIQTR42//VvbHziKdbcPRc1HNY6JCHEAQr5m6lcuQBv/UbqN36qdThCHHKkhEkI0eepqkrZM89R8V60bjp54uEoer3GUQkhDpTREs+gk56gaftS0gafiqqqqJEQOr1R69CEOCRIAiGE6NPUSISNTz5N1X8+AkUh/8pfk/6LY7UOSwhxkGyJedgS81DVCJsXPYDfXUnRcfdJEiFEJ5AEQgjRZ6nhMBvmPUH1p5+BTkfhtVeROuVorcMSQnQib30ZVav/hRoO0lz+La7scVqHJESvJwmEEKLPKnv2+VjyUPTb35By5BFahySE6GS2pHyKj7uPkL8FV/Y41EgYFJ00RxDiIMgkaiFEn5V27DEYE1wU//53kjz0YGVlZYwcOZI33ngjtm3NmjWcd955jBgxgsmTJ/Pss89qGKHo6RJyjiSleDqRcIB1H97Ktm/+T+uQhOjVJIEQQvRZ9pwBHPbk4yRPOFzrUMRPCAaD3HjjjXg8nti2hoYGLrroInJycnj99de59tpreeSRR3j99dc1jFT0Bo1bv6J+4yeUf/sC3satWocjRK8lJUxCiD4j7Pez/i+PkHHiCcQPGQyA3iKLTPVkjz32GHa7vd22BQsWYDKZmD17NgaDgfz8fLZs2cIzzzzDaaedplGkojdIzD2SAROux5ZUgNXVH1VVpZRJiAMgIxBCiD4h7POx5u57qVvyNaV/fpCw3691SGIvli5dyquvvsp9993XbvuyZcsYM2YMBsOuZ2Djx4+nrKyMurq67g5T9DKZI8/H1f9wQoEW1rx9FXUbPtI6JCF6HUkghBCHvJDHy+o5d9P0w4/oLBaKb74RvdmsdVjiZzQ3N3PTTTfxhz/8gYyMjHb7KisrSU9Pb7ctNTUVgPLy8m6LUfRuVStfo2n7N2z6/F7CgVatwxGiV5ESJiHEIS3U0srqO+/GXboOvc3GoFl/IK6kWOuwxF7Mnj2bESNGMGPGjD32+Xw+TCZTu23mtoTQLyNLYh9ljjgfv7uC1IGnoDfZpZxJiP0gCYQQ4pAVdLtZNesuWjduxOBwMGj2H3EWFmgdltiLN998k2XLlvHOO+90uN9isRAIBNpt25k42Gy2Lo9PHBoUnZ68o24FIOCpY91/biF30o3YU+QBgxB706NKmKRVnxCiM+14481o8hAXx5C750jy0Eu8/vrr1NXVMXnyZEaOHMnIkSMBmDVrFieccALp6elUV1e3O2fn67S0tG6PV/R+W5c8irv8WzZ8MgtVjWgdjhA9Xo8Zgfi5Vn3Tpk1jzpw5rFixgjlz5uByuaTThhBir/qfezbB5mayTp6BrX9/rcMR++iBBx7A5/O123bsscfym9/8hunTp/Pee+/xyiuvEA6H0ev1ACxZsoTc3FySkpK0CFn0cjmTbiQcaKX/4dehKD3q2aoQPVKPSSCkVZ8QojMEm90YHHYUnQ6d0UjhtVdrHZLYTz81ipCUlERWVhannXYa8+fP5/bbb+fSSy/lhx9+4IUXXmDOnDndHKk4VBjMToqPfwAAb+M2yhb+icJpd2G0JWocmRA9U49Is6VVnxCiM/iqq/nh9zez6elnUVVV63BEF0lKSmL+/PmUlZUxc+ZM5s2bx0033cTMmTO1Dk30cqqqsuHjP9C07Ss2LfyT1uEI0WNpPgKxt1Z9RUVF7bbt3qpPhqqFEDv5KitZ+YdZ+GtqafzuO0LuFoxxTq3DEp2ktLS03ethw4bx6quvahSNOFQpikLBtLsoW3gfuUferHU4QvRYmo9ASKs+IcTB8u4o58fb/oi/phZLZiZD7rlLkgchxAGxuvoz6KTHMdmSaKlew/qP/kAkHNQ6LCF6FE1HIKRVnxDiYHm2bmPlH2cTbGzEmt2PIXfOxpSYoHVYQoheLhzysfa96wh66jDHZdF/3JVahyREj6FpArF7q77dzZo1i2effZbMzExp1SeE+EmtZZtZecccQs3N2HIGMOTOWRjj47UOSwhxCNAbLBRMnUP5938na+SvtA5HiB5F0wRCWvUJIQ6Gr7KKUEsL9vw8Bs++Q8qWhBCdytX/cOKzx6MoCo1bl9Bc/h3Z466UFatFn6dpAiGt+oQQByPp8HEMvP0W4kpKMDjsez9BCCH2k6Io+N0VrP3371DDAWyJeSQXHad1WEJoSvMuTD9nZ6u+e+65h5kzZ5KSkiKt+oTo45rXrMWcnIQ5JQWAxNGHaRyREOJQZ3Zm0H/8NbRU/kBi/lStwxFCcz0ugZBWfUKIn9L4w4+sufteTAkJDL33bpksLYToNpkjfomqnouiKNSU/huAlOLpGkclhDZ6XAIhhBAdafj2O9be+2cigQCWjHT0dunEJoToXoqi0Fz+HRs+vgMUHdaEHBypg7QOS4huJwmEEKLHq/9mKWvvewA1FCJhzGhKbroB3f+sESOEEN3BmTGclOLp6E0O7CklWocjhCYkgRBC9Gi1i5ew7oG/oIbDJB0+nqIbrkdnNGodlgDKy8v36/jMzMwuikSI7qMoOvKnzAJFh6IoVHz/Mo70oTjThmgdmhDdRhIIIUSPVf/NUkrvfwgiEZKPmEjRb69DaWvpLLQ3ZcqU/WpnuWbNmi6MRojuo+iiv4eq17zN5i8ewGCOZ/g5CzDZkzWOTIjuIQmEEKLHchQVYs1Ix1FUROG1V0ny0MPMnTs3lkA0NTXxwAMPcPjhh3P88ceTkpJCY2Mjn376KZ999hm33HKLxtEK0fmS8qdStep1EnOPwmiT9alE3yEJhBCixzK5XAz901wMDjuKTqd1OOJ/nHrqqbH/v/rqq5k5cyZ33XVXu2NmzJjBPffcw/vvv89ZZ53V3SEK0aX0JjuDZ85Hpzeiqio7lj9PSsmJmOwpWocmRJeST2QhRI9S8d77VH/639hrY5xTkode4Msvv+S44zpeXGvy5Ml899133RyREN1Dp4/Oydr2zZNs/Woeq9++mkg4qHFUQnQt+VQWQvQYO956m01Pz2f9Y0/QsnGT1uGI/ZCQkMCKFSs63PfVV1+RlpbWvQEJ0c1SB56EyZFG5ojzYkmFEIcqKWESQvQI2/75Oltf+gcA/U49BXtersYRif1xxhln8MQTT+D1epkyZQqJiYnU1tbywQcf8PLLL3PbbbdpHaIQXcoSl8WIc19Hb7QSCQep+P7vZAw7B53BrHVoQnQ6SSCEEJpSVZVtryxg2ysLAMg+5yyyzzpjv7r7CO1deeWVuN1u/vrXv/Lss88C0b9bi8XCddddxy9/+UuNIxSi6+mNVgA2fjKb2vUf0FK9muLj/qxxVEJ0PkkghBCaUVWVLS/+nR2v/wuAAef/kn6nn7qXs0RPpCgKN998M1dddRUrVqygqamJhIQERo4cic0mq4aLviV18Ewati4mtWSG1qEI0SVkDoQQQjMNy7+NJQ85F18oycMhwG63k5KSQlxcHMOHDycQCGgdkhDdLj5rNKPOf5uEnCMIB71UrvwnqqpqHZYQnUZGIIQQmkk4bBSZp5yEJS2NjOkdd/ARvcdbb73Fgw8+SE1NDYqi8M9//pPHHnsMo9HIgw8+iMlk0jpEIbqNwexEjYRZ887VuCu+J+htIHvMr7UOS4hOISMQQohupYbDRNqeSiuKQu5FF0jycAj497//zc0338z48eN56KGHiEQiABx77LEsXLiQJ554QuMIheh+ik5PcuFx6M1xuLIP1zocITqNjEAIIbqNGg6z/tHHCTY3M/C2m9EZpdXhoeLJJ5/k7LPPZvbs2YTD4dj2U089lbq6OhYsWMD111+vXYBCaCR96JkkFRyD0ZpA0NtAc8V3JOVN0TosIQ6KjEAIIbpFJBRi3UOPUPPZ5zSu+B732lKtQxKdqKysjGOOOabDfcOHD6eqqqqbIxKi5zBaEwj5mln1r8tY98HN1G38VOuQhDgokkAIIbpcJBik9P6HqP3iSxSDgZKbbiR+6BCtwxKdKCkpiY0bN3a4b+PGjSQlJXVzREL0LHqzE2fGCEz2FGyJss6N6N2khEkI0aUigQBr73uAhmXLUYxGSm75PYmjD9M6LNHJpk+fzqOPPkpqaipHHXUUEJ3jsnLlSp544glOPPFEjSMUQluKopB31K0EvQ2Y7Mn4mssJeupwpg/VOjQh9pskEEKILhP2+1k79z4aV3yPzmSi5LabSRg5QuuwRBe4/vrrWbduHddffz06XXRw+/zzz8fj8TB69Giuu+46jSMUQnuKTh9NHpq2sepfvyYc9DB45jPYk4u0Dk2I/SIJhBCiy/gqKnGXrkNnNjPwD7fiGiZP2g5VJpOJ+fPn8+WXX7JkyRKamppwOp2MHTuWo446SlYWP0R5QyEMOoWWYBCH0UgoomI1yK3F3hhtyZjjsgj5mjBaE7QOR4j9Jv/KhRBdxp4zgEGz/gCqStyggVqHI7rBxIkTmThxotZhiG7gD4d5cd1qFmwsxR0M4DSaODO/mAuKB2PW67UOr0fTG62UnPAX1EgIozWB1tr1GMwOzM4MrUMTYp9IAiGE6FShlhZ81TU48qKTBOMGlmgckeguX375Jf/973/xer2xdSB2UhSFuXPnahSZ6GzeUIgX163m2bU/xra5g4HY6/OLBslIxF4YzE4A3FUrWfPONRitiQyZOR+jLVHjyITYO/nXLYToNMFmN6tmzcFXVc2Qu2fjyMvTOiTRTebPn88DDzyA2WwmMTFxj5IlKWE6tBh0Cgs2dtyKecHGUi4qGdzNEfVeJlsyeqMdo9WFope1cUTvIAmEEKJTBBqbWHXHbDxbtmKMj0PRSZfovuTvf/87M2bM4J577sFkMmkdjuhiLcEg7mCgw33uYIBan4/3Nm9kVEoao1LSujm63sXsTGfwzKcxWhPRG624q1ZiS8xHb7RqHZoQP0k+4YUQBy1Q38DK2++IJg8JLobcfSf2nBytwxLdqK6ujtNPP12Shz7CYTTiNHb8d+00mnCZzPxz0zr+W74ttv2rqgoa/L7uCrFXscRloTdaadi8sG2xuZuIhINahyXET5IEQghxUPy1dfx4+x/xbt+OKSmJoXPvwtY/W+uwRDcbNGgQ69ev1zoM0U1CEZUz84s73HdmfjFNAT9T+/XnmH4DAGgNBrlxyWdMf+8Nfv35h7y4bjVb3M3dGXKvoDfHoShKtJRJjez9BCE0IiVMQogD5q+t48fb/oC/qhpzagpD7p6DJU3KFfqi2267jeuvvx6bzcbw4cOxWvcsv8jMzNQgMtEVrAYDFxRH5zn8VBemm0aMjR1f7fWQ54yntKmB7+tq+L6uhnkrv6O/w8njR0wj1WrT6q30KHEZIxhy6vPYkvJRdAZaqldjTxkoc4hEjyMJhBDigBmcDixpaSiKjiF3z8ackqJ1SEIj55xzDpFIhNtuu+0nb3bWrFnTzVGJrmTW6zm/aBAXlQxutw5ERy1cc+Pi+dvU6VR6WvmiYgeLKrazrKaK1mCQZEs02fxkx1a+rNjBERlZTEjP6rOtYO0p0ZGd6tVvsfGzu8kY/ksGTLhOkgjRo0gCIYQ4YHqzmYG330LY48WUKIsh9WV333231iEIDexs1Zpgjt7sG/dSGJ1us3N6fhGn5xfREgyyraUZXduN8YfbNvNZ+TY+2FbGByechlmvp8brASClr45QqBHCgVZABSSBED2HJBBCiP3SumUr9d8spd/pp6IoCnqLBb3FonVYQmNZWVmMHj0anXTfEvvIYTQyMCEp9vqcghIybHaaAwHiTGYAXt1YyovrVjMoIYkjMrI4IqMfBXGuPvE0PnXQyVhc2TgzRqIoCp66DdiSCrQOSwhAEgghxH5o2bSJVXfcScjtxuBwkHH8L7QOSfQQv/rVr4iPj2fSpElMnTqVI444AqfTqXVYohcZkZzKiOTUdttqvB4UYHVDHasb6nhq9Q9k2OxcVDyEk3MP/ZvpuMxRqKrK9mXz2fb1kxRMu5OU4ulahyWEJBBCiH3jXr+BVbPuJNzaiqMgn+RJE7QOSfQg7733HgsXLmTRokXccsstRCIRRo0axZQpUzj66KMZMGCA1iGKXmjOmIn8ZugovqzcwcKKHXxTXUGFpzW2X1VV/rxiKSOSU5mQlonzEG0jHPTUAyq+5u1ahyIEIAmEEGIfNK9Zy+o77yHs8eAsLmbQrNsx2O1ahyV6kPz8fPLz87nooovwer0sWbKEhQsX8uKLL3LfffeRk5PD+++/r3WYohdKslg5KaeAk3IK8IVCLK2pZGhiMgCljQ28UbaeN8rWo1cURiancmRGP47I6Eem3aFx5J1DURRyjrgR14BJJAyYgKqq+Ju3Y4mXdtlCO1KsKoT4WU2rVrFq9l2EPR7ihgxm0Ow/SvIgflZNTQ11dXW0trYSDodRVZX6+nqtwxKHAIvBwBEZ/XCZo/Ou4s0mflU0iFxnPGFVZVlNFQ/9sJwzPnyH1mB0ITZvKEREVbUM+6Apiq4teYhQtvBPfP/qubRUrdI6LNGHyQiEEOInBRobWX3nXCI+H/HDhzHw9lvQm81ahyV6oAULFvDNN9+wdOlSqqqqsNvtjB49mgsuuIBx48YxcOBArUMUh6AMm4Orh4zk6iEj2dbi5ouK7Sys2IFZr8NuNALw9Jof+HDbZialZ/HbYYdhMfTeWx81EsLbsIVI0IunbgOOtMFahyT6qN77r0gI0eVMLhe5F11A/dJllNx8I7pDtL5YHLw77rgDRVEYPHgwf/zjH5k6dWqf6JQjeo5sh5NzCgdyTuFAwrut4ry8upJan5evqytia0t8VVVOjdfLxPQsEntRFzmd3kTJ9AdpLv+OhJxJqJEQQW8jJnuy1qGJPkYSCCHEHtRwGKXtgzb9uGNJO3YairTnFD/j/vvv5+uvv2bJkiVce+21DBgwgLFjxzJ+/HjGjh1LcrLc4Ijuo1d2/b6aP/kXfFtbTWswEEtqX9mwliVVFSjA0KQUjmxrEZvjjNco4n2nN9lJyJlEJBxg/Ye301pbypBTn8Vkl4U8RfeRBEII0U7tl4vZ/s83GHznHRjj4gAkeRB7NWPGDGbMmAHAtm3b+Oqrr/jqq6+47777qK6uJi8vj/fee0/jKEVfZNLrGZ+W0W7biORU6n0+Spsa+KGuhh/qapi3cgXXDBnJ+UWDAIioamyRu54o5HfTWltKoKUaT916SSBEt5IEQggRU/3ZQtY/8hhEIlS8+2/6n3u21iGJXig9PZ0BAwawfft2tm3bRlVVFTU1NVqHJUTMhcVDuLB4CFWeVhZV7GBRxXaW1VQxqm0dika/j3M+fo/D0zI5IiOLcWkZ2AxGjaNuz2RLYtDJ/4e/eQfx/cYSDvkgEkFv6qOrdotuJQmEEAKAqo8/ZcO8J0BVSZ06heyzztA6JNGLrF27lsWLF7N48WKWL1+Oz+cjNzeXo446ihtuuIHDDjtM6xCF2EOazc7p+UWcnl9ESzCIrW2C9eLKcur9Pt7buon3tm7CqNMxOiWdIzOyOCmnAEMPGZW1xGVhicsiFGhh7Xu/RVH0DDzxEXQGaXYhupYkEEIIKj/4kI3/9xQQnfOQd/llUrYk9sspp5yCyWRi7Nix3HDDDRx11FFkZ0ufetF7OIy7RhiOzc4h3WZnUcV2FlZsZ3trC0uqytnW0szM3EIAVtfXYdDpKIx3ad4wwN+0g9aatSiKgrdxK/bkQk3jEYc+SSCE6OPK33mPsvnPAZAx4wRyL7lI8w9D0fs8/vjjTJgwAavVqnUoQhw0g07HqJQ0RqWk8ZuhoyhzN7OoYjs2gyH2+/HJ1d/zdXUFOc44Xp52oqbzJewpxZSc8Bf0Jjv25EJCgRb0RhuKIg+CRNeQBEKIPizs9VL+9rsAZJ16CgN+dZ4kD+KATJ06lfr6ep544gm+/vprmpubSUhIYPTo0Vx44YUkJSVpHaIQB0RRFPLi4smL29WhSVVVHEYjZr2ebEdcLHl4fOV3lLe2cGRmPw5PyyTO1H2lRPFZowEItNay5p2rics8jJwjfi+/00WXkARCiD5Mb7Uy5K5Z1H31DZknz5APGnHAKisrOeuss6ivr2fEiBEMGjSImpoann/+ed58801ee+010tLS9vl6dXV1/OlPf2LRokX4/X7GjBnDTTfdREFBAQBr1qzhnnvuYeXKlbhcLs4//3wuueSSrnp7QrSjKApzxx2BLxSiKeAHoknF+1vLqPF5+XjHVvSKwojkVI7M6MeUrP6kWrtncrO78gc8dRsJehvJOuxiWSNCdAlJIIToY1RVxbttG7b+/QGwpKeTdcpJGkclerv7778fg8HAv//973ZzH7Zt28bFF1/MX/7yF/70pz/t8/WuvPJKdDodzzzzDDabjUceeYQLL7yQjz76CJ/Px0UXXcS0adOYM2cOK1asYM6cObhcLk477bSueHtCdMhiMLRb2fpP449kYfl2FlVuZ1NzE8trqlheU0WSxcIx/XIIRSKsbaxnUEJSl5U8JeVPoWDqbJwZwzHZkwkHveiNUlooOpfmCYQ8ZRKi+6iqypYXXqT87Xcpvvn3JI0bo3VI4hDxxRdfcNttt+0xcTo7O5urr76aP//5z/t8rYaGBvr168eVV15JYWF0MuhVV13FySefzPr161myZAkmk4nZs2djMBjIz89ny5YtPPPMM5JACM0oisKQxGSGJCZz1ZARbG9xs6hyB19W7ODwtEwAfqyv4YqFH5NktjApox9HZmQxOjUdi75zb8dSSk4EwNu4lTVvX0W/0ZeROujkTv0aom/TfHbNlVdeybZt23jmmWd47bXXsFgsXHjhhXi9XhoaGrjooovIycnh9ddf59prr+WRRx7h9ddf1zpsIXodVVUpm/8cO/71Fmo4TKBW+vKLzhMOh0lISOhwX2JiIi0tLft8rYSEBB566KFY8lBbW8uzzz5Leno6BQUFLFu2jDFjxmDY7cnv+PHjKSsro66u7uDeiBCdpJ/DyTkFJcw7YioOowmACk8rdoOROr+PtzZv4IYln3Psu69xx9IvuySGuvX/we+uoHzFi0TCwS75GqJv0nQEQp4yCdE91EiEjU8+Q9V/PgQg/8rLST/uWI2jEoeS4uJi3nrrLY488sg99r355psUFRUd0HX/+Mc/smDBAkwmE//3f/+HzWajsrJyj+ulpkYXACsvL5cJ26LHmt4/j2P6DeDbmupYi9gqrwdfKBQ75h/r1xCKRDgiox+5u03cPhBZoy8FRU/aoFPQ6Y1EwkF0+p61IJ7onTRNIHY+Zdrpf58yPfbYYx0+ZXrqqaeoq6uTDwkh9oEaDrPh8Sep/uRTUBQKrr2KtKlTtA5LHGKuuuoqLrnkEhobG5kxYwbJycnU1tbyzjvvsHjxYh599NEDuu4FF1zAWWedxcsvv8zVV1/NP/7xD3w+HyaTqd1xZnO0243f7z/o9yJEVzLq9IxLy2BcWgY3DB/N+qYG1LZ9EVXlH+vXUOPz8sG2Mv4xLVqKVOP1kGC27PcCdoqi0G/0xQC0VK2i9IPfU3jsXOIyRnTiOxJ9keZzIHaSp0xCdD41HGb9I/Oo+Xwh6HQUXf8bUo46QuuwxCFo4sSJ3Hfffdx///18+eWucozk5GTmzp3LMcccc0DX3Tkf7q677mLFihW89NJLWCwWAoFAu+N2Jg42W/d0uhGiMyiKQpErMfY6rEa4qGQIiyq2MyQxJbb9D998wSZ3ExPTszgiox/jUzOwG/dvJKH8u78RaKli+9JnGHTS4532HkTf1GMSCHnKJEQXUBR0JhOKXk/RDb8leeLhWkckDmEnn3wyJ510Eps2baKpqYn4+Hjy8vL2uz1wXV0dS5Ys4fjjj0ev1wOg0+nIz8+nurqa9PR0qqur252z8/X+tIoVoqcx6vScllfEaXm7Hp76wiG2tLhpDgR4f2sZ728tw6jTcVhyGlcOHkFJQuLPXHGX/KmzMTkzyB5zGQBqJISi6zG3gaKX0XwS9U4FBQUMGTKEu+66i379+slTJiE6gaLTkX/V5Qy9b64kD6JbKIpCfn4+o0aNIj8//4DWFqmuruaGG27gm2++iW0LBoOsXr2a/Px8xowZw/LlywmHw7H9S5YsITc3V0amxSHHojfw7vEzefLIaZxbOJB+difBSISvqiswtSXYVR4P89f8wLrGelRV7fA6eqOVnInXozfZadiymO9fORu/u7I734o4hGiaQNTV1fHuu++2+xCQp0xCHJyw38/21/+F2vbvStHpcBYWaByVEPuupKSESZMmMWfOHJYtW8a6deu4+eabaW5u5sILL+S0006jpaWF22+/nQ0bNvDGG2/wwgsvcPnll2sduhBdwqDTMTI5jeuGjuK1Y2fwyrQTuXH4aHKdcQAsrNjGM2t+5PxP3+fkD97k/hVL+bqqgmAkvMe11EiYrUsewdtQxo7vXujutyIOEZomEPKUSYjOFfb5WHP3vWz520tseOIprcMR4oAoisLDDz/M+PHjuf766znjjDNoamri73//O5mZmSQlJTF//nzKysqYOXMm8+bN46abbmLmzJlahy5El1MUhdy4eM7IL46N8GU74jgqox8WvZ4qr4fXNq3jN19+ylULP46dF1Yj0fN1ekpOeJj0YWeTM/F3AKht+4TYV5oWv+3+lOnuu+8mLi6OJ598MvaUyWw2M3/+fG6//XYuvfRSfvjhB1544QXmzJmjZdhC9Eghj5c1d8+ledVqdBYLqVOO0jokIQ6Y0+lk9uzZzJ49u8P9w4YN49VXX+3eoIToocanZTA+LQNfOMSy6koWVuzgi4rtjE5Njx1z1cKP0Sk6fpGdwym5BeQe8XsAakrfo3rN25Sc8LCsWC32maYJxM6nTA8++CDXX389breb0aNHx54yAcyfP5977rmHmTNnkpKSIk+ZhOhAqLWV1XPuwV1ait5mY9CsPxBXUqx1WEIIIbqRRW9gUkY/JmX0I6KOJdBWwdHo9/F9XQ0qkGyxcEputKz1w7K1NCx+gX6ejVSv/hcZw8/VMHrRm2g+/V6eMglxcIJuN6tn30XLho0YHA4Gzf6jzHkQmnjzzTdJSEjgqKOOYs2aNdx4441UVFRw3HHHMXv27D266oneLxT0otMZCATcmExOIpEQBnmK3SPoFAVL2zpaLrOF139xEosqdlAQ7wIgFIlw/6ofaY4/mfiECEeF8jiyYjujk9Ow7meLWNH39JguTEKI/aeqKmvuvjeaPMTFMfiu2ZI8CE08//zz3HrrraxevRqAOXPm0NTUxBlnnMHHH398wAvJiZ4rHPazZuUL/OvVabz56jT+9eo01qx8gXBY2qz3RFl2J2cXlDA6JVrW1BoMMjY1A7vBSFNEx9tbNnLjks859p1XWFVfq3G0oqeTBEKIXkxRFPqfcxbm1BSG3jMHR16u1iGJPmrBggVceumlXHnllZSXl7NixQquuuoqbr31Vm644Qbee+89rUMUnSgU9LL6h+dZ9f3TBANuAIIBN6u+f5rVPzxPKOjVOEKxN/FmM/eMncR/TjyNxyZNYWZWGgnhZoiESWpaA8CSynIu/ew/vFC6ik3NTT/ZIlb0PZqXMAkh9p+qqrHuG64Rwxn1xGPoZMhZaGj79u0ceeSRAHz++ecoisKUKVMAyMvLo66uTsvwRCfT6QysW/Nyh/vWrXmZQcMuoq5mJTZ7GlZbSofHiZ7BqNMzNjWDsakZXOhwU9HaRHruEQB8Vr6NH+tr+bG+lidWraCf3cERGf2YmtWfoUny99qXSQIhRC/jr6mh9MGHKbjmSmz9+gFI8iA0l5iYSG1ttOzhv//9L3l5eaSnR0slSktLSU5O1jI80ckCAXds5OF/BQNufN46vvnyTowmB9OmPwfA6h+ep7FhHTZ7Gv36H01y6vDotfzNGE3OA1p0UHSu9MEzSSfa1nXLlw9zkr0fxSPGsLBiB8tqKtne2sLLG9bS4PfFEoiV9bXkOONxyOdQnyIJhBC9iK+qipV/mI2/upqNjz/JkLl3yYeu6BGmTJnCgw8+yJIlS1i4cCG//e1vgejciMcff5xTTz1V4whFZzKZnBhNToIBN2azC4s1GZ+3Fr+/EaPJidmSiKLocTj7xc6pqviaqorouk92RybJqcNRVZU3Xz0WRdFhs6dis6djs6dhs6eTW3ASDmcWkUiYcNiH0WjX6u32OfVln1Hx/d9B0XH8ua9xat7RtAaDfF1dwaKK7UzN6g9EJ2Jf/+V/8YVDnJFXxHXDDtM4ctFdJIEQopfwlpez8g+zCdTVYcnMoOiG30ryIHqMW2+9lXA4zNKlSzn77LO5+OKLAXjllVc46qijuP7667UNUHSqSCTE0BFXYnJkk5E5hmaflziLlfLypQRbtoGqctzJ7UucSgafT0bWRDytlSQlDwEg4G8kEgkA4G7eirt5a+z4zH5H4HBm0dxUxgdvnYnR5MTpzOaYE19EURQa6tbS1LgxlnRYbWno9fIUvDMk5h5N+rBzsKcUY3UNAMBuNDIlqz9T2pIHgCpvKy6zmW0tAeLN5tj2WUu/JNvh5IiMfhTFJ8hn1SFIEgghegHPtu2s/OMsgg2NWLP7MeTO2ZgSE7QOS4gYk8nEnXfeucf2t99+G7PZTENDAzabTYPIRFcwGK0MKD6DF9aU8c9vF+MOhnAaDZxZ0I9fDRyPQa/f45yMfhPJ6Dex3TazJYEzzv8KT2tV9I+nEk9LJZ7WKhzOLAC8nmogWhrlDzTHbka3bfmY1T88t9vVFCzWJDL7TWLsxDsAqCz/imCgBZs9DWdcDiazswu+G4ceRVHIPeJGACLhIJs+m0tKyQnEZ41ud1yW3clrx57EFndzrISp0tPKB9s2A/DMmh9Js9qYlJHFkRn9OCwlDaNuz58N0fsccALx+eefs3jxYqqrq/nd737HmjVrGDx4MFlZWZ0ZnxB9Xuvmzay6Yw7BpmZsOQMYPGcWJle81mEJ0c68efO45ppr9thuNpt55513mDt3LkuWLNEgMtEVvKEwL5Vu5rk1m2Pb3MEQz67ZDIrCecU5WA37dqOo15twxmXjjMvucH9G1gROO3chntYqQiFPbLvdkUlqxhi8bclHOOzH560lGGyNHbN25d+oLP8KgFHjbqJo4NmoqspnH16FxZqIzZaGzbGzbCqDuLgB6A3mPWLoyypWvEjN2rdpKPuckb96G4PJsccxA5xxsf93GE388bDxLCzfztfVFVR5Pby+aT2vb1rPm8edTIbNQXPAjwrEm+R73VvtdwLh9Xq5+uqrWbx4MQ6Hg9bWVi699FJefvllVq9ezUsvvURhYWFXxCpEn7TlpX8QbGrGnp/H4Nl3YIyTJ2ii53n88cdRVZVrr702tq2yspI77riDhQsXMmzYMA2jE53NoFNYsH5bh/sWrN/GhQM7t6W00eQg/n9uXPOLTiW/KDq3RlVV/P5GPK2V6HS7yphciUWEQj48rZXYHRkABALNVFV83eHXmTb9eZJTh+Nu3sayJXOxOdJxOLIYPPxSAIKBFhRF36cWy8sY/kuaK74nfeiZHSYP/8thNHLigHxOHJCPLxxiWXUViyq2U+lpJcMWPf/Nsg08ufp7hielcGRGP47I6Ec/h3y29Sb7nUA89NBDrFq1ir/+9a+MHj2aIUOidYx//vOfueSSS3jkkUeYN29epwcqRF9VdP11bP7bi+T86nwMDplEKHqmu+66izvuuINIJMJ1113HSy+9xF/+8hf0ej2zZs3i7LPP1jpE0YlaAiHcwVCH+9zBEHVeP69v3I7LbKI4wUmhy0m8qevmJyiKgsWSgMXSvrRzxOjr9zhWrzcx4ah78bRW0dpa2VY+VYmntRKbPdo5rMW9PZZk2OwZsQRi7aoXWfX9M5jM8dHRC3t0BCMpZSi5+ScC0a5SBqO1XSLTm+kMZkpOeBhFUQgHPGxa+Ceyx16OJW7vFScWvYFJGVlMymh/bJm7ibCq8m1tNd/WVvPwj9+S64zn7IISTsmVxVB7g/1OIN5//31+97vfMX78eMLhcGx7SkoKV155ZYc1sEKI/eOvq8OclASAwWGn4KorNI5IiJ93+umnY7PZuOmmm3jnnXfYsWMHJ5xwArfeeitJbT/L4tDhMBlwGg0dJhFOo4F4s4l3ysppDAQBGJQYx3NTx9LgD/D6hu0MTIxjYoY2rX0NBiv9c3/xs8e4EgoYN2kOntYqdLpdt0o+bz0AAX8TAX8TjQ3rYtt3JhBLFv2Biu1fYrUlM2L0bxmQdxwA69a8gtWaEi2ZsqVhsSaiKL1jPd+d807KFv2Z2tL38NSuZ9hZfz/g+GeNnsBlA4exqGI7iyp28G1tFWXuJtzBQOyYp1Z/z0BXEmNT07EYZMpuT7PffyPNzc0/Oc8hPj4ej8fT4T4hxL5p+nElq+++l/7nnk3WyTO0Dmev3J4ATS1+Wr1B7FYj8Q4zTptJ67CEBqZPn47NZuO6665j8uTJPPDAA1qHJLpIKKJyZmE2z64u22PfmYXZNPoDnFmYzbpGN+saWyh2RctT1tQ3M3/1JnKcdiZmJNMaDHHL4h8ocjm5cmg+Bp2OiKqi07hrj9WWQm7Bnr9/x0y4neGjf9Nu1MLTWkWcKy92jM9bC6h4PTWxUYhAwM23X/+53bV0OiNWexpHHP0ArsQifL4Gtm/+pG1UIwNXQs97Et9//NV46jeSe+TNB538ZNodnFVQwlkFJbgDARZXlTMsKZpUlre28NzalQCcmlvIzSPHAtAaDGKX9SZ6hP1OIAoLC3nnnXeYNGnSHvs+/fRTmf8gxEFo+G4Fa+feRyQQoPG7FWSeOB2lg24mPUVNo5fHFnzHd6U1sW0ji1O49syRpLj6To1wX3Xrrbd2uL1///589tln/PrXv46NPiiKwty5c7szPNGFrAY9vyrJAaJzHmJdmAqz+VVJDma9nosH7bqpDqsqAPFmI9MHZJBijU6eXd/UwtLqeja7W7l2eCGhSITpby+kv9NOUYKTIpeTYpeTvHg75h7yu9BkcmIyOX/yBv/YE1/C72vA01qFva2TVCQcIDvnmFjC4fXUEIkEaXVvx9g2r6CpYT3Lvor+G7Ha0jj5zPcB2FD6Gtu3/DfarrZt9CLaVSobuyOzG97xLiZ7CkNP/xuKohD0NrD1q3nkTLwBvengOqw5TSZ+kZ0Te60ocGZ+MYsqtjMxPfo99IfDnPj+G+TFuWLzJnKdcbHREW8ohEGn0BIM4jAaCUVUrDJy0WX2+zt75ZVXcs0119DY2MjRRx+NoigsXbqUN954g1deeYUHH3ywK+IU4pBXv2w5a+/9M2ooRMKYwyi56cYenTy4PYE9kgeA70preGzBd/z+vNEyEnGI+/rrjieiAmRkZLBhwwY2bNgAIH3gD0FmvZ7zinO4cGAuLcEQDqOBUETt8EZf3/b3PzgxnsFjd3WRy3ZYufWwgYTaEoyy5laagyFW1jexsr6p3fk5ThtFCU4uG5xPpt1KMBLBqOt5JUCKosNiTcJi3VW6Z7EmMXHyfbHXkUgQr6cGT2sVVlsqEC2tysw+sm1bSuzY+to1VJbv2cEsK3syR0x9CIDvlj5ES/N2bPY0cgpOJCl5MACe1ios1qR2ZVgH//4UVFWl9P0bcFd8TzjgoegX93ba9QEybA5uGD6a3w07DLVt2+qGWjyhECvra1lZX8sTq1bQz+7gpJwCzi4o5sV1q1mwsRR3MIDTaOLM/GIuKB7cYxLPQ81+/0RNmzaN+++/nwcffJDPP/8cgD/96U8kJSUxe/ZsjjvuuE4PUohDXd2Sryl94CHUUIikw8dRdMNv0fXwYdqmFv8eycNO35XW0NTilwTiEPfpp59qHYLQ2M5WrQnm6L91437ezydZzJyct6ssOi/ewSu/OJx1jW5KG92sa3CzrtFNUyDIxuZWNja3cuWQ6JP/qz5bTp0vwB/HDGJkSgJb3K1YDXpSLOYen7DqdEbsjsx2IwhJKUM5curDexxbUHI6yalDY2tltLZW4m2txBk/IHZMVcVSGutLAUhJG0FS8mBCQS9v//P4toQmBftuoxdFg36JzZ5KOBwgGGzFbHbt1/dMURQGTPwdGz+ZQ/bYyw/8G7EPX2dnVCOT03jn+Jl8UbGDRRXbWVZTyfbWFgriXLxQuppn1/4YO88dDMRen180SEYiusA+fUeffPJJZs6cSVpaGgAzZsxgxowZbNq0icbGRuLi4sjLy0PXA58ECNHT1Sz6knUPPQyRCMlHTKTot9f16JGHnVq9wYPaLw49S5cu5auvvoq1cl25ciVPPvkkv/71r6WNq9gnekUhJ85OTpydY/tHOyKpqkqN109po5uy5lZSrGYiqsqGxha84TAuc/Rhy1++W8dXVXXcPnogM3KzWF5dT70/QJHLSbbDpvm8igOVmDSQxKSBP3vMiNHX4W7eiqe1CldiCQBeby2KYkBVQ3g9VXg9VbHj84tPB6Cu5kc+/eAy9HozCUkDmTY9ujBfdeVy3M1b2hKOaOmU0di+C6AzbQjDz34FRafH17yD6lX/Inv8VV06MTzVauPUvEJOzSvEEwryXU01o1PTmbVscYfHL9hYykUlg7ssnr5snxOIMWPGkJaWxsCBA3n11VcZNmwYeXl5ez9ZCPGzAvV1EImQcvRkCq+9qlckDwB268+PkOxtvzi0/Pe//+Waa65hxIgRsQTCYDBQXl7OL3/5S5577jnGjBmjcZSiN1IUhVSbhVSbhSMyo6U9CvDmiZNY3+imvzN6YxuMRNARHcUAeG3Ddv67I7qKtVWvp8DloNjljM2tyItzYNIfGg8+0zPHk545vt02Z1w2Z5y/BL+vPjb3ojXWrjb6QNjni3aVCof9RCK7HvpsLfsPG0pfa3c9o9HBgPzjGT0+Ovdp2+aPCYV9WMyJbPv4LsKt1eiMFvqNvrQr32qMzWBkYkYWDX5fu+5Nu3MHA7QEgySYe8fnam+yTwmEw+HgueeeY+vWrdEVHD/7jE2bNv3k8aecckpnxSfEIS/r5JOwZmWRMGokyj6M4vWUrkfxDjMji1M6LGMaWZxCvENWGO1L5s2bx0knncS99+6qhS4pKeGNN97g1ltv5aGHHuLll1/WMEJxqIk3GRmdmhh7/fjkw/CFwxjaRhry4u1UeePY2DZS8WNdEz/W7ZpXMX1ABneMHUxZcwvfVNUzNCmeQYnxe3yd3kyn02O1pWC1pZCUMnSP/f1zjiHr/KPweqoJh/yx7fEJBWT2mxRbJyMYcBMMtoCqxo5Z8+Nfqa9bDcDggjOIbF9BUtHxfPju+dFOUruNXthsabgSizAYLJ3+Hh1GI06jqcMkwmk04ejh5cC91T4lEJdddhn33Xcfn3zyCYqi8MQTT/zksYqiSAIhxF7UfvElrpEjMNijT84SRx+2T+f1pK5HTpuJa88c2WE8vzlzpMx/6GM2bdrE73//+w73nXTSSVx11VXdHJHoiyy7jeBeNjifywbnE4pE2NbiobRtPsW6trkVRW2tZZdU1vHo9+s5KiuF+yYMZ0NTC8+t3sTIlATOKMhGVdUeP6fiYOj1JhzOfu22FZacSWHJmbHXwWArntYq9PpdD4ZS0kZiNDnwtFaSUTSdxMNvpLW1gvraVdTXrtrj60yf+QZx8TnU167m228ewO7IIM6Vx+BhlwDg8zWgoGAyx+/X9zsUUTkzv7jdHIidzswvJhRR93tujti7fUogLrjgAk477TSampqYOnUq8+bNY+DAn6/HE0J0bMdb77D5ub/iLC5myN2z0Zn27Ua7J3Y9SnFZ+f15o3vEiIjQVlxcHJs2bWL8+PF77NuyZQt2u6yiLrRh0OnIjXOQG+fguAEZQHRexc7Wspl2K0dkJjOmbTRjZV0jn26vpiUY4oyCbLa4PVz12XKKEqItZYvayqCy7NZeO69ifxmNduJd7cvWR469YY/jlFCItJCVQNhLfP5kIkbLbqt8R8um3M3bqK1eQW31ChISS2IJxKoVT7F+7QL0ekt05KJtle+09DHk5E8HwOupwWC0YzTuahtrNRi4oHgQoLJg47rdujAVcUHxIMx6mUDdFfb5u+pwOHA4HNx7770cdthhJCQk7P0kIUQ72197gy0v/h2A+KGDUfZjaLWndj1y2kySMAiOO+44HnnkETIzM5k8eXJs++eff86jjz4qHfpEj6IoSqzUaXJWKpOzUmP7hie7uGZYAWnWaLnNukY39f4AX1XW8VVlXew4m0FPYVtCcdXQAqwGPWFVjbWs7YvsCbnkDjydlqpVlBwxB4PZuccxKWkjmXDUn/C0VmE07XqwEAi0ABAO+3A3b8bdvBmItsXdmUB88ekN1NWuxGSKY9ykOWT1P4pg0MOW9W9xjCObC44/iWa/hzizjfLyb9hS+hp5hSdjMMq6RJ1tnxKIN998k6OOOoqEhAQURYm1b/0pUsIkRHuqqrLt1X+y7eVXAcg+5yyyzzpjv4ZppeuR6Mmuu+46fvjhB6644gqMRiMul4vGxkZCoRDDhw/nd7/7ndYhCrFPdo5W7HRkZgrzp4xp11Z2Y1MLnlCY72sbWd/o5rcjilBVlZPeXUSSxcz9E4eTZrNQ1txKqtWM3dg3noIrisKACdejhgPoDGaatn+DqkZwZe8ambTZU+mfe+we5x5+5N2MnfBHPJ7q2KRvT2slCYnFsWOCwWiSEQg0Y2gbhdDpDPy44kmCATdmswuLNRmftxa/vxGjyUlB8Wld/K77pn36ib7llltYsGABCQkJ3HLLLT97rMyBEKI9VVXZ+tI/2P7aGwAMOP+X9Dv91P2+jnQ9Ej2ZzWbjH//4B59//jnLli2jqakJp9PJ6NGjmTx5srT5Fr2WxaBnSFI8Q5J2TbAORSJscXtY1+imORBEpyhUtHqp8wVo8gdJtERHZa9f+C3VXj/PTBnDkKR4FpXXYNTpKHQ5SLIcmo0mFEVBMZhxV/7ImnevQ1F0DDn1OewpxXs9V28w44zLxhmX3eH+6TPfIBhowdNahc0RLUcLBtwEA24A/P5G/P7G2PHBgJtAsAWLXqpmOts+JRCzZ88mLi4OgE8++aRLAxLiULN9wWux5CHn4gvJOnnGAV1Huh6Jnk5RFCZPnsy4ceNoaWnB5XJhlA4o4hBk0OnIj3eQH79rpCLdZuGN6RPZ3uLBqNPhDYWJACqQExct1fm/HzewqbkVgGSLKTafosgVnV+RabceMhO27SklxGUdht5gwZqY22nXNZocxJt2fd9NpjiMJmcsiWh/rBOT0bHHdnHw9jmBUBSF/v37M2nSJCZOnMi4ceNkUpwQ+yDp8PFUvPc+2WedQcYJxx/wdaTrkejpli1bxv3338+PP/6I2jZBddiwYfz2t7/tcHK1EIcSRVHItFvJtEfr7a0GPe+ceARN/gAOowFVVSl0OQlFVLa1eKj1BaitrGPxbvMqrhlWwHnFOayobWRHi4cRyS6yHLaf+pI9mk5vpPi4+9HpjSg6A3UbP8GWVIjV1b9Tv04kEqJo4Dms+v7pPfYVDTyHSCSETi8PMjrbPiUQ77zzDt988w3Lly/no48+4u9//zsGg4Hhw4czceJEJk2axNChQw+ZrFmIzmTrn82oJx7D4Dj4hFu6Home6ttvv+XCCy8kOzubq666iuTkZKqrq3nvvfe49NJLefHFFxk5cqTWYQrR7eLN0d/PiqIwZ9wQADyhEBsaW1jX6I7NrdjU3EJBfHTS8Xuby3mnrJwLB+ZwxZACFlfU8tmOaiZlJHPkbhO+ezp92+Tlug0fs+4/t2BypDHsjBcx2hL3cua+MxitDBp2EQDr1rxMMODGaHJSNPAcBg27qF3rWdF59imBKCwspLCwkF/+8pcAbN26lW+++YZly5bxxhtv8OijjxIfH8/48eOZNGkSZ5xxRpcGLURPpobDbHxqPskTD8c1fBhApyQPO0nXI9ETPfzww4wePZpnn30W/W69+K+55houueQSHnvsMZ577jkNIxSi57AZDAxLdjEs2RXbFoxE2PkYNi/OzshkF0PaFrb7pqqet8vKMev1HJmVysLyGh7/YX20BKqtDKrY5cRl7pmfDXGZo7C4+hOfNRqDpfMX69PrzQwccgGDh11CINiCyeggEglJ8tCFDqgtQP/+/enfvz+nn346AF999RWvvPIKH3/8MR9++KEkEKLPioRCrH/4UWoXfUntoi8Y/cz/YXBI/aU49P344488+OCD7ZIHAJ1Ox3nnncfNN9+sUWRC9A7G3RoNnFM0gHOKBsReH5WVglmvY1RqdDJwaUMzW9wetrg9fLStKnZcqtVMkcvJkKR4LhwYnXfQExbCM9oSGXLa8xjMcSiKQs2690kYMKnDNq8Hamer1p0TpqVsqWsdUAJRWVnJokWLWLx4Md988w319fXExcUxbdo0Jk2a1NkxCtErRIJB1j34F+qWfI1iMFD4m2skeRB9ht1uJxQKdbgvGAzG5kQIIfbfyJQERqbs6iR0RkE2Q5LiY21lSxvdbG/xUu31U+31U+8LcOHAXOp8fs75YAnFCXE8fORI9IrCVncrmXYrhm7ujGZsG3koX/ESW778C86MkQw6+f/kRr+X2qcEIhAIsHTpUhYtWsQXX3zBxo0b0ev1DBs2jF/+8pcyB0L0eZFAgLV/foCGpctRDAZKbvk9iWNGax2WEN1m1KhRPPnkk0ycOLFdg42WlhaefvppRo+Wfw9CdBaX2cTh6ckcnp4c29YaDLG+MZpQ7Fx3Yl2jm+ZgiGqvD72i4AmFOOuDJZh0Ot6dcQROk5FPtlWRZrNQEO/AYtD/1JfsNPFZo9GbHLiyx6Ho+sb6GIeiffqbGzduHD6fj+zsbCZMmMD111/P+PHjccjTVSEI+/2snXsfjSu+R2cyUXLbzSSMHKF1WEJ0qxtuuIFTTz2VadOmMXnyZFJSUqipqeGzzz7D7/czd+5crUMUXSAQ8qLXGfAF3FhMTsKRECaDrPqrBbvRwIiUBEbsNlIxOjWRF6aNpSUYHR0sb/VhNeixGww4TUZCkQizv1lJMKKiA/o77bH5FDvnVsSbOneEwJ5SwohzX8NkT0FVVWrX/4ek/Gkouq5PXkTn2acEwuv14nK5GD9+POPHj2fs2LGSPAjRpuKd96LJg9nMwD/cimvYUK1DEqLbDRgwgFdffZV58+axcOFCmpqaiI+PZ9y4cVxzzTUUFBRoHaLoZKGwn8VrX2Dp+pfxBd1YjE7GFp7DxIEXYZDJqz2CUaejOCEu9rog3sHHp0ym3hcAoDkQ4rCUREob3TT4A2x2t7LZ3cqHWytj5/zliBEcnp7M4opaQpEIw5JdBz1Z22RPAWDzFw9Q+cMrNA9eTt7k2w7qmqJ77VMC8eabb7Jw4UK++OIL/vWvfxGJRBg6dCiTJk1i0qRJDB8+XMqXRJ+VecpJeLbvIP3YacQNGqh1OEJopqCggIcffljrMEQ3CIS8LF77AotW7+q97wu6Wdj2+vCSC2QkoofSKQrJ1miCl2gx8fCRI6MjAb4A63drK7uu0c2OVi+5cdEHxn9dU8YPdU3MGjuY4wdk8G5ZOWXNrRzdL7XdKt37Iy5jBFUr/4k9paTT3p/oHvuUQJSUlFBSUsKvf/1rPB4PS5Ys4YsvvuCtt95i3rx5xMXFcfjhh8cSioyMjK6OWwhNhTxe9BYzik6HzmCg6PprtQ5JCM0tXboUvV7PqFGj2LFjB3PmzKGiooLjjjuOq6++WuvwRCfS6wwsXf9yh/u+Wf8ykwZd0s0RiYOhKAopVjMpVjMTMnbNq3AHgjja5lMUuZx4QmFKEqKdkz7eXsVXlXVkOqwMSYrnlXVb+WhbZXRV7bbVtfPjHZj1P12alFRwDPbUwVjiMlEjIRq3LiEh54iufbOiU+z37BWbzcbUqVOZOnUqAFu2bOHLL7/kq6++4q677iIUCrF69epOD1SIniLY7GbV7Dtx5OeRf9UVXTr65vYEZNE40Su89dZb3HLLLVx00UWMGjWKWbNmsXz5ciZOnMiTTz6J0Wjk17/+tdZhik7iC7jxBd0d7wu6cXtr+PSHeZiNdtITikl3lZAaX4DRYOnmSMXBcO42/+HGUe1HCU7IySDLbmVE21oWP9Y1sqq+mVX1zbFj9IpCjtNGUYKTcWlJHDdgzwfM0eQhzLr/3EL9pv+Se9StpA85vWvekOg0Bzz9va6uju+++47ly5fzww8/sHr1aiKRCKNGjerM+IToUQKNTayaNQfP5i0EamvJPvMMzCnJez/xANQ0enlswXd8V1oT2zayOIVrzxxJiktKA0TP8vzzzzNz5kxuuukm6urqWLx4MTfccAOXXHIJzz33HK+++qokEIcQi8mJxejsMImwGJ3YzAmUVX+Fx98Y264oOnJTx/LLo54gHA6ytfY70l3FWM2dv7CY6HrHZKdzTHZ67PWVQwuYnJUaLYFqK4NqDATZ2NzKxuZWAI4bkMG3NQ3c9c0qJmQk8/tRJYQiEer9QSyuPBTdF5hsXfOZKjrXPicQGzdu5Ntvv4392bp1KwBFRUVMmDCBX//614wdOxarVW5sxKEpUN/Ayj/Oxrt9O8YEF0PunN1lyYPbE9gjeQD4rrSGxxZ8x+/PGy0jEaJH2bRpE7feeisACxcuRFXV2Ej10KFDZW7EISYcCTG28BwWrn6aJOcAUuLyqGneRJ17C2MLzyEY8nHcyJupbCylsmEtlY2lePwNGPTREYjq5g289PkVWE3x3HDyp4DKF2ueIzW+kPSEYuKsaTK3spfp57DRz2HjmP7RpEJVVWq8fkob3axvdFPoipY+lTY0U+HxUev1A7ChqYULP/6GHOcInj/7eMyuHD7cvI3sSDXFuaPQyc9Bj7RPCcTYsWNxu92oqkpWVhaHH3441113HePHjycxMbGrYxRCc/7aOlb+cRa+8gpMSYkMuWsO1qzMLvt6TS3+PZKHnb4rraGpxS8JhOhR4uLiaG2NPmX8/PPPyczMJCcnB4CtW7eSkJDwM2eL3sZksDJx4MWMy70As9lCpDWCzq7D7/disBgw6E0M7v8LBvf/BRC9mWzx1RIIeQHwBZpx2bNw2TNRFIV69zY+W/lE7PpWk6ut9KmY9IQS0l3FJDr6o5NWn72Goiik2iyk2iwckZkS2z4jN4silzM2N2JHixd928Rua0IuZfW13LG0FADrd59QmBAfaylb5HKSF+fApN9zETxvKIxBp9ASCOEwGQhFVKzdsK5FX7VPCcT48eOZOHEiEyZMIDs7u6tjEqJH8VVVs/KPs/BXVWNOTWHIXbOxpKfv/cSD0OoNHtR+Ibrb+PHjmTdvHuvXr+ejjz7i4osvBuA///kPjzzyCJMmTdI4QtHZ9GEj4Q+CeD6qB48KNgXDsRb0Jxrhf+7bFEXBad11E5mbNo5rT3iHcCT6u0wFhg04kcrGtdQ0l+ENNFJW9TVlVV/HzjHqLaS6CjlpzGyS43Jpaq3AbkmUlrG9jMNo4LDUXQ+fp2anMTEzmeZA9GfBG9GRp2tgW9iBN2Lkh7omfqhrih1vUBT+efwEMuxW/ru9mjSbmfx4By+VbmbB+m24gyGcRgNnFmbzq5Kcn53ELQ7cPiUQjz76aFfHIUSP5dm2DX9tHYaUVGxX/Y5anZ14T6BLRwDs1p9fuGdv+4Xobrfffjs33ngjjz/+OBMmTODyyy8H4N577yUzM5MbbrhB4whFZwr5IoTf8xJ6y7tro0cl9KYXVNBPNMP2MJgUFLOCYlPQZRtQIyq0qGAEzAp6XfR3WZKzPyePuzN67bCf6qaNsdKnqsZSqhrXEQz72FH3I1aTC4AFi2+gqnE95xzxCPnpE9hSvRwVlXRXMRaTs5u/I+JgWPR6LNbojf6g5EReOvkkWht3UGtKZXV1Nevq69no1bGu0U0oEiHNZkFVVe5dvoY/jB7EFxWbeW51Wex67mCIZ9ten1ecIyMRXUDWEBdiL8IFg/h25AwW1xpoeWk1sLrLJzPHO8yMLE7psIxpZHEK8Q554iZ6loSEBJ599tk9tv/jH/8gM7Pryv2ENvQGBf9Hvg73hT7yYTrRhueuJnCrACgZemz3JUCriuea+l0HG4gmGSYFTIBZwXyJg8y8QSQvyWHQpmkYJlvQHa2ncekOfKvcmEvtqMMiKHU6Ssonk7Q+h1BVgFWrP2S7+3tCej9Wh4vExP4kZmaTllBEuqsEpzVF5lX0EjqDGWdyHubWGpqW3kKav5krT30Ws3MkDf4gOkXBGwozMT2JMWmJ3LV0VYfXWbB+GxcOzO3m6PsGSSCE6IBn61Z0ZjNBZ0J0MnNjXLt/LV09mdlpM3HtmSM77ML0mzNHyvwH0WtI8nBoUlvVaNlSRzwqakuEHTlgrYQkvR5dSrRmXQ38zzkhIKSi7n6tcPQ/kXUhwt8F0A8xodOZcGxNwvyZjbApgGG4iQv6zyfw7xZYAX6aOYrL2l06QpiHTzgeFLj8w1cJRXyoxjB6ixGTzQ4mFcWsQzEpGI61YhhmIvSdn/DqIPqBRgyjzEQqQoTXBFGS9RiGmVC9ESJlITArKEYFzNHEJ5oAKWBAkpROpOjbRttVlUjIj6IoJFqin39Wg55Z44bQ4AvgDoY6PN8dDNESDJFwkCtniz1JAiHE/2gt28zKO+agt1hIvuFmzSYzp7is/P680bIOhBA9nKqqfe6mUbErYFM6TiJsCkq8jlkD3cw8KZ4zi6MT6G//ooL8eBOTH3aQazdDQI0mFAEVNQAEVPCr6LKi5SaGqRb0Q43oi6K3KroCA4ZjLOiLozeVujgd+iHGXdfw0/bfCJFABBWVIQOOp6pxHeaQHUPEFE1YvKA2RLMUlQgA+nFmvP4mvN83Yf7UhqJTYJSZ8IYQgb+2oh9mxDDMRKQijO9PzfwkBaz3J6BL1eN7tJlIRbgt0dgtyTApKCbQDzNhGGMmsj1E+IcASroewygzamuE8Logik2HvtiIGlFRqyLtr2HsG4mK0eJi0EmPEwn5scT3w++uRGe0YrTsav3rMBlwGg24gyFcJiNJFhN1vgCNgSBOoyG2EJ7oXPJdFWI37vUbWD37LkItLVhSU/BGfr5usqsnMzttJkkYhOjhFEWhoaGBUKjjp6CHIpc9AcOxluich/9hONZCKBDm5iEWnMYANTU17GgN89m2Fj7bBhWNLfy6xAbAhuYQAxx6jDYFbG0XaG37k9H2B6AG6Nf2B2/0dRpw0U9FqICqZ5JyPQDeHD+NzRtobtpBsrmAeEMmq8reoaLmB/JcR5CbfCTr1n/KxpZPyS4YRlNrDf6FXnLrx9O/eBj6NAvumgBKE5jTgCAoQSDQ9v+Rtpt5lej8DkCtCqPuCPMT4zQocToY05akvOJBP9IUHfUoD+P/ixslVYftgUTwq3hvbtjj7e2cR6KYFMyXO9GXGAm87SG8NhhLNBQTsXkomBR0WXoMY82ozRHCq4MoDgX9EBNqRCWyNYxiAiVdj6JToomZgWgypSGTPTr53tu4hdVvXYnJnsKgk/4PvSn6AxOKqPx6cB79zTZGpScSbAljdOj5trKerX4PoYiKcc+mTeIgSQIhRJvmtaWsnnM3YY8HZ3Exg2bdTmVr5GfPkcnMoi/z+/2YzTIfByAUCvWpBMLtaybhxOjIQuhDX7suTKYTbTS0NJBmVgGVUCiCQ69ydbGZpXUhxibpCIVC+MIqty5tRa+DU/ubOHVAFz4sidcTHz+A+OwBAAQIkZ9zLMmewRj0FoKWEL6KFioy17A5dWn0nEoo5SMojL6MW5JBkiOPpDMKGN7vdBSl7a5UBcJgUA0kxSWhxEVvuE2XO6ElghoE/LtGWwhES7n0BW0jK6k6DBPN6Aa03ZIZFXR5BpSEtusHAasSPbetvAuVaPISUFFRQY2mKZEtISIrf/rBln60CcPYtiTlCXfb3BQTtKr47mgEwPZ8EgDeG+pRm9To3BTjriRkZ1JiPN6KYYyZ0DI/4e8Cu5V1sdu8FgUlXsEw0owaUolsDEUTmRw9iqKguiNgbCsD20uiooaDhINeQn434WBrLIGwGvScPiCbwDseAh81gEclYFMYcayVsScmoRgO/ZEaLWieQDQ2NvLQQw/x2Wef0dLSQnFxMTfccAOjR48GYM2aNdxzzz2sXLkSl8vF+eefzyWXXKJx1OJQ07RqFavvnEvE5yNu8CAG/uE2DDYr8UpAJjOLbuH2BHpdudqUKVOYN28eI0eOZN68eZxxxhmkpaVpHZboBn6/nwYacB4Xh+0kG6pHRbEpBANhGloa8Pv97Y63GxSOzjBydMauhy4VnggOo0JDQMXedpOnqioPr/FT4NQxNtlAmrXrHh3rFD0J9v6x1wMzplOSfjwt/hrqWjZS17qRupZN1LVuotVfQ7OvgmZfBTXuDYzIPpNwJMiCZZeRYOvP1IG3YrA4aTHV4lSSAQV9/327xdIPMqEftOvfuj7HgHW2K/ZaidNhfyp6U6+GVAjuSkJ2Jia6jOhoufE4K4ZRpmhJmF9FDaq7JS+g6982qm5W0A0yomtLUtQwKAk61LCKom/7uwi0BbAzUWmNJik7R1TUlrakpSxEaFH7v+/dKf300QSiRcV3TxMoYPtr9P14bm6AFhXrvS6ULAO+h5qIbAm3n1fSlrTozKkMznsO09EudDVOPG9WYnA4MIyxEFri77AjmAIYT7BFkx/RqTRPIH73u99RV1fHQw89RGJiIv/4xz+45JJLeOONN0hMTOSiiy5i2rRpzJkzhxUrVjBnzhxcLhennXaa1qGLQ0Tz6jWsnn03kUCA+GFDGXj7Legt0dVSZTKz6A41jd4Of8a6stNXZ3C73VRXVwPw+OOPc+SRR0oC0Yf4/X78/hr0ej0Gg4FQfYhwOLz3E9vkOvU8dbiNTe4IqZbojezW1ghfVof4shrq/SoXFEQf0mxrjZBpU9B3cd2/oig4Lak4LankJB8e2+4LNsWSiZ0jDw2tm2n11xAK+zDoLKhqhP/74DRAIc1VSLqrpG0xvBJS4vLQ6w9+xFoxKGBQwBqtYvpf+kIjFO796+hzDVhv2TWPQOfSYXuk/cLAtocTIAiqv/08lZ0jKbq2BEk/1AQWpd3oyq7jQElqSwIjKkq6HnafM7RzUr2pLWlpUlEboiP/HZV+GZJSMVrtuH/cgO4/LoJOP8bj7AR/oiNY8EMfxpNsHe4TB0fTBGLLli18+eWXvPzyy4waNQqI9hJfuHAh7777LhaLBZPJxOzZszEYDOTn57NlyxaeeeYZSSBEp7FmZWJOS8WckkLJLb9H/z8lGTKZWXQltyewR/IAXd/pqzMMGzaMG264gfvuuw9VVbn66qsxmTqOVVEUPv74426OUHSHcDi8X4nD7nSKQkHcrrlmLpOOiwpMLK0NMzY5eovSGlK5cZkHh0HhsCQ9Y5INDEvQY9Z331NlizGerISRZCWM3BWrfQAnDX8Qb7ABRVFo8dcSjoQIRwJsq13BttoVsWP1OiMpcfmx1bWLsiYTb+vaBUkPlmLV/WSisjt9iRF9yd6TFl2iHtuf269Ib3s6KVqi1Xa6+RpndKQjNjFebZfE7ExaQs56mtM/x5ZdiLVl4s92BMOjQpyMQHQ2TROIhIQEnn76aYYMGRLbpigKqqrS1NTEypUrGTNmDAbDrjDHjx/PU089RV1dHUlJSVqELQ4xxvh4ht5zJ3qbDZ2x41+CMplZdJWmFr9mnb4O1oMPPshf//pXGhsbefPNNxk0aBCJiYl7P3EfSHlr3xRvUjihn4kT+u3atr01glkHTUGVTytDfFoZwqSDYQl6ri2xYDdqc3No0JlIjSuJvXZaUrn51EXUNm+mqrGUysZSKhtKqWxciz/YQmXjWiob1wKQFJdLvC2dz1b+H3XuzYwpOIv+KaOIRELodJoXh3QbRdc2Z6KNLkUPKT99/E4JR45FyQnhTB+BotP9bEcwbJI8dAVNf0rj4uI46qij2m17//332bp1K5MmTeIvf/kLRUVF7fanpqYCUF5eLgmEOGC1i5cQamkl/dhpQDSJEEILe+vk1dWdvg5GWloaN998MwBff/01v/3tbykpKdnLWftGylvFTsXxep6daGdNU5iltWGW1Yao8atsckewtd3FLK8LsaUlwphkPf1sOs1anOp1RtJchaS5ChnGiUB0Xkdja3k0gWgopapxLemu6L3N+vJFVDauZVD2sQB89P1fWLP9Y9JdxaTFSqCKcdmz+kTb1v3h6j8BgEBTA4ZjzITe2rOMyXisJTrxvO/kZN2mR31Lly9fzm233cbUqVOZMmUK99577x7D4Ts7fvzvBC0h9lXN54tY9/CjoKpYMzOIHzJY65BEH7a3Tl69pdPXp59+CkBzczMrVqzA7XaTkJDAsGHDcDgc+3UtKW8V/8uoUxiWYGBYgoGLC0xsaY1Q799VS/9ReZBldWEWbIbnJtqxGaKlTxYd6DVuQ6ooCgmOLBIcWQzsN7XdvmnDr6e8fjX9koYBUNW4Dre3Bre3hvUVX8SOMxsdpLuiyURaQglZiUNIjsvpzrfRYzVVLyV5xjRQIPShf7eOYGaMJ9qik7FFp+sxCcTHH3/MjTfeyPDhw3nooYcAsFgsBAKBdsftTBxsNpkUI/Zf1SefsuGxJ0BVSZ1yNHEDO+dpqRAHKt5hPmQ6fT399NM88cQT+P1+1La2kkajkSuuuIKrr756n68j5a3i5yiKQo5DT85ueem4ZAMRFXQK2No6Or22OcB/K4OMTDRwebEZSzfOmdhXuWljyU0bG3t99qRHqGpa11b6FC1/qmnaiD/Ywpaa5WypWQ7AwH7HcPqE+6hqXM/SDa/QP2UUwwacoNXb0FRizhFUrF6AddgA4k8cTaQlgM5homnLUurXbiV14MnojT23GUVv1SMSiJdeeol77rmHY445hgceeCA26pCenh7r8LHTztfS6UPsr8r/fMjGJ54CIO0Xx5J/xWXR2kkhNHSodPp6/fXXeeihhzj99NM56aSTSE5Opqamhrfeeot58+aRmZnJzJkz9+laUt4q9tfOFrE7E1eAtU1hWkJQ2hzG3ParfnldiGpfhNFJBlIsPe/3v8loIzt5BNnJI2LbwuEgNe4yqhrWtiUVpfRv27+tdgXfbfoXzZ4qhg04gZqmTfxz8Y2kJ5RERyza/mszuzR5P91B0RnYvvRJwn43BosLky2ZgKeWkK8RvdlJ2mAZlewKmicQ//jHP7jrrrs4//zzue2229DtdkM3ZswYXnnlFcLhMHp9tEvDkiVLyM3NlQ8IsV/K3/03Zc88C0DGidPJvfRiqScVPcah0Onrr3/9K+eccw6zZs2KbcvLy2PcuHFYLBb+9re/7XMC8b+kvFXsq91/r989ysq6pggtoV2lTu9vD7KiIcyLGwM8P9GOWa/gC6uYdfTYzwS93ki6q4h0VxHD/2dfVuJgJg68mCRnDgAVjWuoc2+mzr2ZVVs/iB0XZ0vfrQQq+t94W0aPfc/7IxxwE/a7AQj5Ggn5Gnft87sJB1rQWRN+4mxxoDRNIMrKypg7dy7HHHMMl19+OXV1dbF9FouF0047jfnz53P77bdz6aWX8sMPP/DCCy8wZ84cDaMWvU3zmrWx5CFr5skMuOD8Q+KXpji09PZOX1u2bOGWW27pcN/UqVN5/fXXD+i6Ut4qDpReURjo0rfbNjJJjy+i4jAosTawL20M8E1tiDHJekYnGRiSoMeo8byJfZWROIiMxEGx14UZR3DOEY/FOkBVNa6lvmUbzZ5Kmj2VrCv/PHbs70/5HIvJyeK1L+C0plCUeSRm4/7NV+oJ9CYnerMzlkS022d2ojf1vvfUG2iaQPznP/8hGAzy0Ucf8dFHH7XbN3PmTP70pz8xf/587rnnHmbOnElKSgo33XTTAT/FEn2Ts6SYrFNPQTEY6H/u2ZI8CNEF0tLS2L59e4f7tm3btt8TqUHKW0Xni7aINRHerdRpZWOI+oDKf8pD/Kc8hEUPIxL1TM8yMeh/EpCezmqKoyBjIgUZE2Pb/MEWqhrXxVrKVjauIxT2YzE5CYZ8fPrjPFQ1zHUnvo/Z6ODdZXehoIstgpcaX4DRYNHwXf08NRIiY9jZbF/6zB77MoadjRoJQScs4ifa0zSBuOKKK7jiiit+9phhw4bx6quvdlNE4lChqipqKITOaERRFAb86jxJHIToQlOmTOHRRx+luLiYESNGxLZ/9913PPbYY0yZMmW/riflraIr7b6i9Z8Ps7GyMczS2hDL6sI0BFS+qgkzPjkC6AmrKh/sCDIq0UCGrefNm9gbs9FB/5RR9E8ZFdumqtHVnoNhH4fln0Zjyw6c1lQikTArt35AMOSNHasoepKdOaQnFJPWVgaVnlCC1RTX7e+lI3qjlaxRFwFQ8cMrhP1u9GYnGcPOJmvURegMvacRRW+iqLvPOOoDpk6NtlD75JNPNI5EdBVVVdnywou0bCpj0B9uRfcTK+N2Brcn0Kvr1kXf1BW/B1taWjjnnHPYsGEDmZmZpKSkUFNTQ3l5Ofn5+bz00ku4XK59ulZZWRkzZsxg8uTJ7eZUwK7ypeOPP54pU6bEyltnz57NnDlzDniE+mC/JzU1NYRCoQM6V/QcETW6vsTS2hAzsk04jAprGsP8cUX0hvriAhPT+5lix+p2S0QMBgMpKfuwCloPFo4EWbv907YSqOikbY+/ocNjjx91K6MLzmBH3UpafHVkJQ3BYdEugQ8HvSg6A+FAC3qTAzUSku5L+2l/fg9qPolaiM6kqiplzz5PxTvvAdDw7XckjR/XJV+rptHbYeeca88cSYpLfmmJvsXhcPDaa6/x+uuvs3TpUpqamhg2bBiXXHIJp556KhbLvpdASHmr0IpOUSiI01MQ1750aViCnlWN4VhJU70/wk3LvIxK0vOLLCP5zt5V6vRT9Dojg/v/gsH9fwFEP1Pd3hoqG9e2lUFFk4rG1h0kOfsDsHzjP/l+8zscMegyJg+5ko2VSyir+jrWASrR0R+druu/PzuThdiEaSlb6lKSQIhDhhqJsOmpZ6j84EMA8q74dZclD25PYI/kAeC70hoeW/Advz9vtIxEiD7HbDZz7rnncu655x7UdaS8VfQkA1167nBZaQ2p2Nrug5fXhWkMqnxaGWJUkoF8J4QiKu9sbGJCpp0k66Fxe6UoCnG2VOJsqRRlHhnb7gu4MeijpUFxtnRS4wvIbJvMvb58EUs3vBI71qi3kOoqJN21q7Vsanx+7HzROx0aP+Giz1PDYTY8/iTVn3wKikLBNVeRNm3/aq73R1OLv8OFvyCaRDS1+CWBEEL0CTqdDp1ORyQSIRKJaB1Ol7EbdpUrTU43kGZVWFYbZnhCNKtY1RjiT982AXDbuDROyIvOEVBV9ZCbg2cxOWP/P3nIlUwecmXsdV76eCJquG3UYj3BsI8ddT+yo+7H2DE6xUBBxiTOmvQQwZCPHfUrSXcVt7uu6NkkgRC9nhoOs/6RedR8vhB0Ogqvu5bUyUfu/cSD0OoNHtR+IYTo7QwGA0abFZvFgjvoI95owePzEfR4D/n5IEadwrAEA8MSdt1GRVQYlGRmdZ2fwUnRkr2KliC/+e8OJmXZmZRlZ3iKFUMvaRF7oIoyj4yNVkQiYepatrS1lC2NtZf1Bhox6KMP2SobS3nxs1/jsCTz25M+JBwJsnjtC7HRCocleZ8SsEDIi15nwBdwYzE5CUdCmAxSTtxVJIEQXaI7Jxf7qqqoX7YcRa+n6IbrSZ44oUu+zu7s1p+vrdzbfiGE6M0MBgPOBBcvrP+aVzd9izvox2k0c1beKC4oHIe7ofGQTyL+18gkI8eWZFLnDZFoiY5KfFHeSnlLkAWljSwobcRp1DE+M5pMTO3vOORGJv6XTqcnJS6PlLg8hg44HoiOyDR7qwiHo2u5+IJu4m0ZpMTnA1DbvJnPVj4Ru4bdnEiaq6jd6tqJjmwUZVdHrFDYz+K1L7B0/cv4gm4sRidjC89h4sCLpFSqi0gCITpdd08utmZmMnj2Hwk0NJA0bmynX78j8Q4zI4tTOixjGlmcQrxDfmEJIQ5dRpuVF9Z/zfzSJbFt7qA/9vqs/iMINe+5sFdfsPv8hxPz4kizGVi0vYXF5R4a/WE+2uJmS3OAaQOi5Tpr6324zHrS7X3jwZOiKMTb0mOvCzMmUXjie4TD0ZF7naJnSP/jqWwspc69mVZ/PZuqvmJT1Vexc0wGG2muIk4aOweHJZnFa//KotW71oHwBd0sXP00AIeXXCAjEV1AEgjRqbprcnEkEMBbXoE9ZwAAzqLCg77m/nDaTFx75sgOE6XfnDlS5j+IPqu1tRW73Q7A+++/T2VlJUcffTQ5OTnaBiY6jU6nw2ax8Oqmbzvc/+qmb7moaDz+ltZDek7EvrAadBzZz8GR/RyEIyqr6nx8saOVDPuu268Hl9Wwus7HqFQrj03tp2G02tK3dU1Kic9j5vh7AAiGvFQ3bYiVPlU2llLdtJ5AyMO22hXYzQnodQaWrn+lw2t+s/5lJg26pNveQ18iCYToVN0xuTjs97Pmnj/Rsn4Dg++chbOw4KCud6BSXFZ+f95oWQdCCKJrN1x++eVMnz6d66+/nr/85S88/fTTqKrKww8/zHPPPcdhhx2mdZiiE+h0OtxBH+6gv8P97qCfOn8rT29ZhgWFPFsiBfZE0i19e4KsXqcwLMXKsJRdT8ODYRWTTkGnQPJuIxfzvqvBG1I5IsvO+Ey7FuH2CEaDlaykoWQlDY1ti0RC1Lq3UOfejNnooNVXjy/Y8WiXL+jGF2zBbk7orpD7DEkgRKfq6snFIY+XNXfPpXnVanQWCxF/xx9g3cVpM0nCIATwwAMPoNfrmTp1KsFgkJdffpnjjz+eO++8k1tuuYWHH36YF198UeswRSeIRCLEGy04jeYOkwin0YzLZOWL2jIaA9EF2EocKfx58C+oD3h4rXwVefZE8mwJZFvjMXbDGgE9lVGv8Pi0fjT6w3iD0dGaUETl3U3NuAMR1jX4YwnEhgY/yTYDLnPf/X4B6HQGUuPzSW2bM2ExObEYnR0mERajE4vR0d0h9gmSQIhO1ZWTi0Otrayecw/u0lL0NhuD7riduIElB3w9IUTnWbp0Kffccw9Dhw5lyZIluN1uzjrrLBwOB2effTbXXnut1iGKThKJRPD4fJyVN6rdHIidzsobRYOnhbMzh7KptZ5NnnpKHMkAbGit592q0tixBkVHf2s8ufZE8m0J5NkTKXYko99tgmxf4DLr2yUGsw5P54sdrRS4dj2gmvt1Fesb/RzZz8E9kzK0CLNHCkdCjC08h4Wrn8ZmduGwJNPiq8Xjb2Rs4TmEIyH0ur4xv6Q7SQIhOtWBTC7el45NQbeb1bPvomXDRvR2O4Pn3KFZ6ZIQYk/BYJD4+HgAPv/8c6xWa6xkKRwOYzDIx82hJOjxckFhdKHOn+rCND2taI/zUsx2Tk4vYVNrA5s89bSGg2zyNLDJ08AnRBOKV0efhV6BJ8q+JtVs57jUIhyGvjPSa9ApHJ5p5/DdSpd8oQhhVSWigsO4K7l6/LtaFAUmZdkZnGRBf4i3iO2IyWBl4sCLGZx1Ai5HOn5PALPNRKO7Eld8WqxdrOhc8htddKr9nVy8Lx2bgs1uVt0xh9ayMgxOJ4PvnIUjL7d73pAQYp8UFxfz4YcfkpOTw7///W8mTZqEwWAgGAzy97//naKiPW8mRe9mikQ4LymTiwvG4Pa14LQ48DfswPQzE6dzbQlcMmA0EG3nWeVvZZOnnjJPA5ta61FRMer0NAV9fFC9HgWYnlYMwB/WfIxJpyPPlhgrgUozH/qtUAEsBh0vHD+AytYgoYgKQCAc4c0NTXhCEf6+pgGXWc/ETDuT+tkZk27Daug7ozg61UjLijTWfRMg5AODJUDmuDQSJ8nIQ1eRBEJ0un2dXLyvHZv0FjPG+DiM8fEMvms29gH9u/PtCCH2wW9+8xuuuuoq/v73v2MymbjssssA+MUvfkFtbS1PPvmkxhGKzuS0GPF/8TzBz54hZHNhcCTjaalF9TTin3wZzjFn09Dy8+tAKIpCusVBusXBhMT2v9cVFC7IHkl9wINNb8QfDrGyuYoIKssay2PH2fVGcm0J5NoTOTalgAE2V1e83R7jf1u93jw2lS92tLKkvJVGf5j3ypp5r6yZR47OYnS6DX8ogjsYaTdB+1ATDqhs+zLI1oW75liGfLD18+jr7AlG9KZDP8ncXXesfn7o/kQJTe3L5OJ97dikM5koue1mAvUNWDPSOzxeCKGt8ePH88477/Djjz8yfPhwsrKyALjgggsYP348xcXFGkcoOotOp8NstdH8VbR1puppRPU0xvYHvnqFuEkXYH7nXsKWeNSEfkQSB6Cm7vsoVJzRzGmZg2OvDTod9w46NjanYlNrA1u8jbSGg6x0V7PSXc1oVxYDcPFk2Tesa60lz5ZIbtu8ilxbAhb9oXXLY9LrmDbAybQBTkIRle+rvSza0cryCjeZSgs1Na0sqQ5w3w8eCuP0XFFiJT+u934PVBVQQdFBOAA1XxuJ+PQMOcNB+TcdN2gp/zpI/0lGamtrUVW1ewPWiMFgICGh67tO9d6fJNHr/VxHprhgCzVvv03WJeeiKAp6s1mSByF6sNNPP52rr76a6dOnt9t+wQUXaBSR6Co6nY6I143q67h1pupzo7Y2YKjfjK56AwCRlEL8p/8lekDAg/HLp1Fd/Yi4slBd/VDj0kH/0+UmekXHQGcKA50psW3BSJjt3uZoQuFpIN+WCMDalho2eRrY0FofO1YBMi1x5NkTOCIxh/GJ2d3ylLa7GHQKh6XbOCzdRk0NhEIhQhEoa46OAq1vDuPURwiFQmxrjfBReZAxyXoGxusx9MB5E/46HYEGHYFGPYEmHcFGHYEmHamTfMQVBVGBhpVWbCk6gh6VkK/j64R8EPSpqKra51ZG72qSQAjN/FRHpvigm3N2fIRnSwvlSQ6yZp7czZEJIfbXtm3bcDikXWJfEIlE0FnjUSxOVJ8bxeZCcSSjtpUwKRYniiMJ/6DpKCkbURq3oyblxM5XGndgWPtxu2uqig41Lh3VlUXE1Y/QiNNgL+VIRp2eXHsCufYEpu62/ZbCI9nYGk0qdo5YNAR97PA1s8PXTLYlnvGJ2XxUs4GXt//ItNR8ftlvOJ5QAHc4QKrJfsgkFmfmmDgmw8DapgiJ5uiciK9rQvx7R5B/7why5wgrg1x6/GGVkAp2Q/e8bzUCoRaFQFM0SQg26TAlhHENjj5YrPrcir92z3a1gcboe1B0kDTGj8kJRrsVgyWaLBhtYHQoBFtUgh4wWMBoUYi09u0FDbuCJBBCMx11bEoINHNO+YfEhTyY0tNJnjRRwwiFEPvqhBNO4KmnniIzM5Ps7GytwxFdKBKJ4Pd6ME+5gmBCBpa8sbh9jTgtLnybvsHYUIHf5yU08FgY2MEFrPEEx/wSpWE7usbtKE3lKEEvSlM5NJWj37KU0MjTY4eb/j0HxdscG62I/jcLNT4TOujOlG5xkm5xMjFpQGxbQ8AbSyhGxEdboG5qbaAu6CEQiT6ZXt5Uzv0bvsChN5FrT4hO1m4rgepnjeu1rWUTzDoOT90V+0CXninpBtY1hymOi25fWhvisbV+Brv0XFxgpp+9c95rOADBRh1hv4I9OwxA3TIzDd+bUMPtkxVbv1AsgbCkhlEMKqb4CCZXBOPO/8btSgQSRwQA8Hv9DJhsxBoHrnwDEXcQndNI48YQ3mbwef19fkX0riAJhNDM/3ZsSgo0cs6Oj3CEvZgyMhk+905MibJ6pBC9webNm1m2bBnHHnssFouFxMTEdvsVReHjjz/+ibNFb+MJRrCPOZW/r3+NVz65BHewFafRztk5J3LRmNNpbWr5yXNVZyqh0efstkGF1jp0jTtQGneguKvAGh/brasqRfE1o6subX8dlLZrnUu4pG0MwtsEkRDYEmG3UYQEk5XDTFYOc2XGtv0qewRHJucQZ4i2F68PeDEoOlrCAX5sruLH5qrYsSZFT4kzhbsHTkNVVda31jHA6sLcC+dVDHbpGezStyvh2uCOEFbhx4YwzrbigC0tYZbUhBidZCDPqUO3l1GZ1m16Ag16dCaV+JJoIlD5sQ3PdgN6S4S8X0V/JnRGFTWsoOhUjPG7kgNLSjh2rdRJP1GT1AFf0EPmqHiCH5Thn78NvCGwGnBM7k/CcTk0uH/6Z1EcuN73ky8OKTs7NtWs3UDlg39GDXuxZGcz9O45mFzxe7+AEKJHyMjIYMaMGVqHIbqJzmLg+fWv8cz6V2Pb3MHW6GsFzs4+AfZcpLpjigKOZCKOZOg3vP0+VcV/0lyUxu1tCcb2XYlGoDWabOx2Y2tY/R+M3/wN1WhtK4dqm2PRVhqlxmeC0QKAzWBikDM1du7JGQOZnlbEVm9TW+lTA2WeespaG/BGQgQj0Rvc2oCHG1d9gA6FLGtcbJSi0JnC2HgH8SbrgX1Tu9nuZVoXFpj5RaaR9e4w8abo6MPimhCvbwny2pYgiSaF8fEGRhnN9A8biDRH5yZYM8IkDI2OBNQusRBo1GNODscSCJMrjL9ehyk+QiQEOgM4i4LYc4MYHSqdMahj05kJflBG+P2yXRu9IcLvbwJFxXZkJn7/vv4win0lCYTQnJUQtQ/fj9rixp6Xy+A5szDGObUOSwixH+69916tQxDdRKfTYTNbeWXzux3uf6XsXS7IP5VrS+8k25zOtTnnAfBlw7dU+GqIMzo4JmkCiqLgCXsJRsI4DLaOS4QUBTUpBzUph3ZFKKoK3kZ0jTuIuLJ2bfe3oCq6aElUzQZ0NRvaXS6SnI//jEeiL4I+9Gs/iiUYqiMZo05Pvj2RfPuuEbSIqlLpd+MNR0ud6gIe4g0WmkI+tnmb2OZt4vO6zdGDV/2HDGscTx9xDum2OFbWl5NscZBmdfb4eRUZNh3pZh2BxuifQU06diRHWFEf5shtcYxbGV3Yrm63cyIqJAyN/r8tO4QpMYI5MQSBMKiQPN5PygQ/ujo/SnkIJRjBGFQhGEEJRQgnmQn3j17X/GUNuuZgdF9QRQlGIBTBPyaJUHEcqCpxj66DYITgoHi80zOjHcFsFnyfbevwPYX/uw3LcbnovDopY+pkkkAIzRlsNgb86jyqPvyYQXfcjkEmYgrRK5SXl5OSkoLRaKS8vHyvx2dmZu71GNHz6XQ63MEW3MHWDve7g600BJoIEKQ20BDb/nndUhY3foddb+XY5Oj8tg9qvuC57a+joPB/Q2bRz5LOZs8OXtjxJvEGB06Dnbjd/nu4awSKohBUw+iscWBrX+YamnAxoXHnozRXthu1iM632EHE1S92rNK4A9MXT8VeqwYzanxmdKSibeQinDMWnclGpiUudlyJM4W/jTqN+qA3tgDeptYGNnkbqPS5qfO3kmyJ3hTf9M1b1PhamH/EuQxPyuKDbatRFIXi+FSyHQndO68iEr1xJ6ASaVEx6MOgQGvITsP3JvRVXoxNPnSRCDo1QqEpzLhBPvyJJn7UKdAESZ4txKstGO166g4fwP/VtDLhn7X8Ymsz+ZG2G/+wCu9BMN9B67k5AFg+qcJU2rxHSP7RiXjbEgjTyib01XuWLoUKnYQAFAXFG4p+DV90NEin06G2BKJlSx3xhlA9oWjnMEkgOpUkEEIzu9dfpk2bSurRk1H0u7ouuD2BvS5G1116UixC9BRTp07l1VdfZdiwYUyZMmWvT1jXrFnTTZGJrhSJRIg3OnAa7biDrbhMcSSbE6j1N9AYaMZptJNsTuCK7LPRqbt+JobHFWPRm9vdNHvD0RtGFRWH3gZAVaCWpU0/7vF1TYqRNw57DIAPahbx1LZXsett3F9yI/2tmezwVfFqxfvEGRxMThxLQe7hhFWVlS3riTM4iDPYcSrmXTc+ikI4Z3w0wWiuQAn5UerK0NXtKoXxXvAimKJxGb98BkKBaIKR0I9kVxZJcemMbhsBMRgMWF1xbGttwKDT4w0FcJmsNAa8FMZH28/OL13ClpZoe1mL3kihM5kiVxrFrlRKrCnkqg5MIQWCEdRAOPokH9APjZ4f3thIZF09BCIYjhmAYjMS2eYm+MGm6LGBCFaPHwJhVBXqziom5FFIKK/G9p+KPb6n4RQz7ukDad1qpH/zNtK8le0P+AJ0xU4GnpKAMq4Z5z+bMG5sIaIY2JHv59u6AEM8YayeDm7gg7tu2NV4I+FEE6pRB0YdqlEBg45wqiV2jH9kAjpPCNWoaztOQTXoCGfsKglruTAPVa+gti2MF4lEUBJMYDV0nERYDSg2AxGvJA+dTRIIoYmmH1ey+a9/Y+Afb8PkcgG0Sx5qGr17rFI9sjiFa88cSYqre+tLe1IsQvQkc+fOjXVcmjt3bo8v0RCdIxKJ4PF7uaL4XHJMGRyWOpRgiw+jw8Kyqh/YEqzE4/dSYsttd94JqZM54X+u9cusGZyVMR13uJU4w/+3d9/xUVX5/8df906f9EoooUuoQugtokF0XXf1x6J+RWEVbOu6YMUCuhTF8gVEUVkbLrqosArr166IuwrKUgIoVUQgSAmBkDLJ9Ln398ckE8YECJhkUj7PxyOPZO6duXPuMCT3Peeczwn2Pre3teYv7a7H4S+jxF+Kw19Gsb8Uhcr3V4k/ODG2LODEZghehOZ5jvNlwX8ByIhqT+eotrg0Nw/98HTYc9oNVmKN0Tzc+XbaX/Ywx7wneO/ISmIDOoO1ODqVulAKD5LrOobNqxHrdWGIsWL48SuUMheKvw2KbkYzHkY3O9Fj22BwD4eUbliSO9LJG8DrLcDg01hy+TUEki0Yjrtxv7aZZx1t0bytMfh0LAEFs6ZyxYVbWW4KMCIvjke/71Bl6ohuN2Kfe1Hwtd9TiP+DnwAwDGkFNiPe4x70zfmh+59c/PTAO1EoJojvHf5/M4CKYlXQbcE5CynDXFiOGPAeigarUn6RH7yID7SwYij/vMwzIAlf9zh0i0pLu8rfBtvZ3kph5/BEWscZcSkwabObMkVhYJqJO8ufr/SSNAyXtqz2/VTBOzDptPsBAmnhf3ODFcE8GC5KJ/Dxvir3N1yUjscpVZjqgqI3l6X5yo0cGazUsGrVqgi3pPkq2vIdO2c/ieb1kvbb39DptlvC9jucXuYs2VjtKtWZGSlMGde/3j79b0htEaK2yO/Bqn7ta3Ls2LFmtVCVxWIhISZY+Sbwn8rKN4aL0jFd2oFCR3HdTVz1aWheH2WaiyKji1bWVExHvRQWFfBD0U/4PB7Ot51HkhLLidQA9zleoMRfyrhtvYn1WbEEjFgCRvraMrBqJrYP9HCH81livVZWfDsRi98YHId/0lOW3NoJpehblu7/nlvWjwLgs7Yf83PcZuJ9Olfterjapvq7fImhgx3F1BE+TK32Pm/dYGKzP5+oH0t5ZEs7DFYjRZqHQt2N26BTagzw7LBCepva0l+z0jPXhN0fxSGtLaUlRlSXl3hPAZpiwBhjIGWEl4CqU/SjhaMHEzAl6LQe5cBo0PD5DQS8KqY4HbUOPkIucGss2etl0wk/4ztauLiVCV3XmbTeSSubSnZLE4NTaveJjUYjSXEJ+D/fT+Df4e9F4yXtKSgubFb/N41GIykpKWe+YzXO5veg9ECIenViYw67npyD7vOR0L8fHSZUXaW2uNRT7QU7wOYfjlFc6qm3i/aG1BYhGrq8vDw2bdqE1+sNbdM0DZfLxcaNG5k/f34EWydqk91kxfdZNZVvyj8FjhqShv9IEbpRQUsJ9hAYcsswnPBA+QTZ4CRZnUBrG75uwap7tk8OByfc+rTK+/k13CNa4O0TnO8Q9787UDSIPj+emCuDcxrsnxwh9qCTdlR8il0GlKGMSOXFC2YCEPvlTtSyylKh4AN8JPaO4aq0S3G5nFhdKoRP1yag6CgaBDJG8v2x9fwUXYDH4OerJPg6LdjufqVHaJPYAUtSCosP/Au36qOTy8nwY+sJnCjmpVbJJHZpT5LPzcjiUozte3HkorEUUcrQ2CguLYKEnlnYbkgj4IItX+dy9IgTtdhErDOGO9cHC4s81PNlvCl+rlOG0W9L8BLObzBxLC2FxBYWTCkqgYxi/H4/UW2ho1paXqRKRUfFCBipu8+Nk6wqd3a34td0tPKnyS3TyHPp5LkCtI1SQwHiP3k+OkSrtI1Sf1Xvpd/vp6C4kJgLWmH9TQd0px/FbsTj9DS78FCfJECIelOwdh0/zH0a3e8ncfAgMu67G9VUdTXqMpfvtMc50/7a1JDaIkRD9sknnzBlyhT8fn/oYuDkeU4dO3aMZPNELVJVFYvNgvvfp6l8M6oDMUsP4EsxUzYuOJTJsqEA885qJtL2SwwFCMNBJ8a8qhNpFddJF/4mFTwaBCovhLVEMwFv4KTx82owvCRbKp9nSAoE9OD4e5OKbgwO0YltZePGuNGg65Tc5gmN0ddNKpoBAqqOqfzj+v/X6wo2dy2mxO8l1X8el/pb4fCXceiq9qQkdEWPs7Lo8/UEdI3xLS7gEst9BI7n8m3plxxU8+hZ6mOEowRDzHm8UfIJ3+3ez0VHf0equyXvps/jiH03Y4rTGbjvf/nl/5iA2c/vk3qz3rcHYktJu9LFTj2fJ/Z9Qr8W6bww7BoOlBZy77qP6WCPp6M9kU5RCbS1xWNSq67qXJeMamUgaBelMn+AjQ3HA/RNCrbD4dNZuMuDBqRaFfonGemfbKB7nCHssTXl9/sp9BejutTghGmXJsOW6pgECFEvjq/5hh/mPQOaRvLwYZx392RUY/Vvvyhb1VBxNvtrU0NqixAN2UsvvUT37t2ZMWMGb775Jn6/n1tvvZWvvvqK+fPnM3Xq1Eg3UdQSVVXRnf7TV74p9aK3tKOf9Csy0MqGz6+HLs4rLvIDbeyh+7hHpKJ4tLCLeEwqWmzlgYrvygCjCiddaDqvrKyudCqeIcmnv4OioJ00qRdALf+q0C+ux2kPEWeOYd3lK3B4nXhLwFlmJb63gTsO98K51UL8wc5sN5jpN6iYZNdKOlh0+hYOASDN3Zrc6B/xuH8gVv8ai36Qt1sVsSM2n4RkE6+MepyLyGbmFx+xJP8jlgB/O/9RrlB64ddLuParOwkEEtnvsLPTcdKcCEWhrS2e1tZoWtui+F2LrsQao+ptzpKiKKRHGUiPqgwxZX6dvkkGvi8MkO/W+fiQj48P+bAb4OkBdpKt51adStMkONQXCRCizmk+H7lvLgVNI+XCEZw3+Y6wCdO/FBdtITMj5ZTzDmxWIwfzHVUqItVFpaQztSUu2lLNo4Rofvbt28fcuXPp3r07Q4YM4dVXX6VTp0506tSJgoICXnzxRYYNGxbpZopaoGkait14+so3cWYc17ULu5jzDE3BM/T0x/Z3iT39HQDM9ftpek34ShVO7DJxzOXGWaDhKgA9AChuhk+zc0mr4Rw95uOH/waH93m1FkzudgO+9jp5qX5siTA76QYcZYMxtskjsehnAsf30925hwR3GS073R56LvOhXSRHQ0BR6LJoAl2s0bzTNo3/i3ODfhiUOBSiQI/GQAwB3cA+ZyH7nIXo+Hn7yKu8k/ks03atwqQGKPAdIsFs5IbWo+gZ245Sv5P1xd8TY4gqr14VLKMbZbDVWuhIs6k82MuGJ6DzfWGAjcf9bCwIYFYhyRJ8jlVHfKw56qd/soHsNBM245mf22AwYDQa8fv9BAKBM95fnDsJEKLOqSYTPWY+wtHPVtL2umtPGx4AYuxmJl2TWW3loz+P6c3f3v2OdduPhrYP6tGCm6/sxcLl39V6paTTtWXyNZky/0GIcqqqEl9eUa19+/bs3bsXTdNQVZWsrCxWrFgR2QaKWtPcKt/oOvhKgour+YrDv7f6jRNriobmVSjYZAIqL1oVA9gSFXxOHUuMQkJnI31uUrElqZhswYthk00hfWhF70oiCSlDoH3lc19L+FBAXdO4K3Ukk47vRzu+H01RwF3Khfv20cFmpNioUti9O0WJrSnxl5JQmM+lh4vZE9eKt81e9uoODKoJh9/LnrIT5c8STb4bppasJt60gRYWG7vKvkdX8kCpnAhvUFTe6j2PKKONrY7dvH/0S2KN0YxOu5g21jTcAQ/fO34g1hhNS0sqcaYzr+lkMSgMSDYyINmIpusUeCrP9b/H/GwtCrCrOEB2WvA1OuLUKPPrdIxRUU8KMxaLBVNUNFajCYcvQILJgNvvw1dWKqtQ1xEJEKLOuI7kYWuZBoA1NZV246+v8WNT4m1MGdc/rEfBZjVWCQ8A7VrG8fw73/Hdj+G9BJt/OMZz/9z8qyslVdcWWQdCiHAdO3YkJyeHAQMG0L59e3w+Hzt37qRHjx6UlJSETawWjZ/DXUbSJe0BTln5pjEKeKF0nwlfkYqttZ+oNsFAcODdKPRA1U/AvUUq1hQNU5xGfHcfiW2isCcr2JJVrHEKyknDrMxRCuaos+89OflTf0VVsV12T+i27vOgnfgZ+/Fc0o/nEji+H2eb3+JPaAuAae8CjLnraQVcUP4YjwL6lpt5IqULOakd2BaTRJ7HRbHPT5HPTZHPjUI6bWxGPHohJ9wtCOhGdPUAdoOVEp+bXNcR1hZtAQgtCpjnOc6sPQsBuCX9aq5sEazoM/H7qZhVc2gtjooejUuSh9Ha2gKf5mOP8wCxxmgSTXFAcBjZxM4WNhb4KfHpod6Hjw/5+OSQjwSzwtz+duLMChaLhai4eN7YfoLlu4tx+DRiTCpjusTxxx6JUFwkIaIOSIAQdeLw+x+yf/EbZEy5h6Qhg8/pGDF2c9hF+sF8R5XwANC1XQL//GJ3tceorUpJv2yLECLctddey/Tp03E6ndxzzz0MGjSIqVOnctVVV7FkyRJ69Dj92HHRuDTWyjd6AFBBUcBTqFL0nRlvsUrSAA/2VgF0v0L+V8Fea12DqDYBFAUsSQG0gII5TsMcr2GKC2COD/4MoBqg5YU+UlLqd16cYrJgaNEZQ4vOoW3OY8eg/PX39/odWouM8lW5DwZX3y7JQ/GW0f/QZjItdryDJwLgOfYThz+exY8J6eyNbcGfvdEo8V0Zp5RRqPuZ1ulWFEXhxf0b+G/hQVJNFxFrUthaXEJAz0fBR2d7Wxz+svIgAB7NS773RNWGAwPietHa2oKj3gKm7JoDwI2tR3NVy0sBmJP7vxhQGRDfC7gMgCPufExqLCgBDnqOUBKIJsXWleXbT7B4e2Vodfi00O1rM2IlQNQBCRCi1h1c/i9y31gCQOlPe885QPzSqSoeeX2n7yaXSklC1L2rr74ar9fLwYMHAZg1axa33nors2fPpnXr1kybNi3CLRS1rSFXvtE1cB814C0KDjXyFqv4ilR8DpW2V5VhSdDQ/VCyO/jBkKfAj71VAINNJ6qdD2O0jq1V5XCkNlc6aYzrJOrJHQkk/6KeU8CHUpKHUnQILJXDjGylx+juLKC7swAOlT8eeMoaxx57AoO2vg/XvUSex4Ff1znu9XHcC3vLtgJbUYCW1h6MbtGVrMQMSv0ePFqA57s/Qom/tPyrDEf5z2mW4KR2n+YnzZJMib+U2PJhTwFd48eyXHR00m2VC9AVGt/CYTtAqZ7AAz8cp11Ua95q+zzLdx+u9vyX7y7mxh5JuAwGmRNRyyRAiFqj6zo/L3uHn99eBkD6tdeQfu01tXb8U1U8MptOX61BKiUJUT+uv75ymGLbtm355JNPKCwsJDExMYKtEnUtUpVvNP9JcxPKQ0JMZx9R6QHQ4eCHdtCrXvX7ilQsCcHeg8R+bszxGtYWwYtLRYFWl7qqPKYxhodTMpjQE9LRE9LDNmsdBuOa8BZq0SGU8t4Ktegg7YoO0b7oIKhG3LZ45vW4jKOeMn7+eCZ7VSN77PHssSdRYLJy2O1AO/QdqqeEVYqBRUd2kGiy0TEqkY72BDpGtadvbAItLNGhoVkd7G14tddjQPA6osLjGXfj8JeRZIoPbesZcx6xxujyMKLQPa4zpb4AjlN8kOjwaZT6AhiNRgkQtUwChAhzrpWMdF3nwJK3OPhucKJk23HXkX71mFpt26kqIu3KLaT3eSlV5kCAVEoSor599dVX5OTkUFxcTHJyMkOGDJEA0cSpankPRB0ECV2HgFPBW6QGJyWnBS8CD31qw3nACIRf2RujdaLSA8H7tgx+N8cFMMVroeFHBnvwIlU1QVI/mZsTxhqLlhYLad3Ct2sBFOcJUBQUIM1so63BQFbhzyhHtgJQaLTwky2BNh4HFm8ZRYMmoAAnfC5OFB1iY9Gh0OGiVCMDEtK5p/MwfFqAg+4S0q1xGNXgB4IGRaVXTJcqzZvQ5g9htw0GA3FmIzEmtdoQEWNSiTYZKGygQ+oaMwkQIuRYkavaakNnqmSk6zr7F7/B4ffeB6D9xBtofeUVtd6+U1VEyj1SzF+u7l1tFSaplCRE/SgqKuKWW25h69atGI1G4uPjKSoqYuHChWRlZfH8889jNsv/xabEaDQSZY3BarPgc2uYrCpul4cyt+Oc5kBofijLNYYmJke1DR4j991oNI+CPd1H68uCvQPBdd0UVLOOKU7DHB+ckxA27Oh3zto4TQGgGtCjU8Jue698MvizpxSl6BBRRYfoXd5roRUdZEKr7lyT0pl9zkJyt7zHvpIj/GRPYL81jjLAv28t5m3vszehLXdZ40g223kt8w8EdI3P8vfQ0Z5Ae3sCVsOpL1UDgQBuv4+rMuL4+7aqE/evyojD7fdJ70MdkAAhgGDPwy8vzKHmlYy08glKHW+9mZaXX1Zn7TxdRSSplCRE5Dz++OMcOHCA559/npEjR6IoCpqm8cUXX/DII48wf/58HnjggUg3U9QSo9FIQlwSB7/xcXi9E78bjFZoNdBIm2FJFBYXVAkRug7+UiU0H8FbZMBbrJJwvic07ChvVXBRudhuXqLa+lEUMCcECLgUjNGVw1uSB7lJGerGYNOb1vCixsgSjd4ig0CLjCq7bED3mFR6tuqBQfejFB0k8PMmftb8KOgYXMUUO4uwdxxGui048frw0d28mLsRAEXXaW0w0skaS/vYNDrGt6KTPYFYU+WCf7rHzfjuieg6Vaowje+eiKestF5ehuZGAoQAoLjUU+1iaXDmSkaKotDx1ptJHj6MuJ51X2nlVBWRpFKSEJHz1Vdfcd9993HxxReHtqmqyiWXXMKJEyd4/vnnJUA0IVHWGA5+4+PA15Uhwe8mdDu5dyyHdpZgtOnYK0qhLo/Ce6JqGVN7Kz9R6QFUE0S182Gw6NhbVh63ze+rTmA2xeiIxkPrNAytU+VCkq28TpTiw3gLDzKg6BDL4lvj6BhcZTBQkseA4sP8ZE/ghMnGQS3AQWchXzkLIW8nAC93G0VabAtW5u8hu0N3vtm/mfOSE/hX93YUeXzEW0z8N/8A7+3P5XdtukFZRE67SZMAIYAzVyr65X49ECDv089occkoVJMJRVXrJTwIIRqu5OTkare3bNkSp1OGkzQVqqpitVk4vL76f9PD6/2kD7NzYp0dSwsf9jbBYUemaA1vkYoptnI+gilew9aicnhJk5/ALILMdvSUzgRSKsvPVtSDat+uP9Pt8ShFBykqPMReZyF7fW5+MprZY0ugxGghNToZXddZnreTq3sM4qVd3+DweYg320iyRlHgLqPI6yLGZOEPHfpgkCpMtU4ChADOXKno5P2a38+Pzyzg+OpvKNm5i4z77jnNI4UQzcHo0aP529/+xsCBA4mKigpt9/v9LFmyhNGjR0ewdaI2qaqKz63hd1e/3+8Gn1Mnpr0GtsqLthYXulDNoJy+cJ5o7izRaG36QJs+xAC9y7/wuYO9FmUFqKoBd8DPRakdKPK6cPiCw6iLvC6KvJUh1OHzUOp1SxWmOiABQgCnrnAE4ZWMNJ+P3fPmU7B2HYrBQNLQIfXdVCFEA2S1Wtm/fz/Z2dlkZ2eTmppKYWEha9asIS8vj7i4OB566CEgOOzx8ccfj3CLxbnStOCEaaOVakOE0QrmaIWUC5xhVZkM1qr3FaLGTFb05I6Yyte1sBqMXNe6NwnWaGJMllCIOFmMyUK02UphiYxhqm0SIARw6gpHJ1cy0nw+dj01l8ING1GMRro+cB+JAwdEsNVCiIbi/fffJzo6OAhh3bp1YfvS0tLYtGlT6LYiY1IaNU3TcLs8tBpoDJsDUaHVQCNul6fBLConmq5AIIDX5+F/Ovbl1R/WVtn/Px374vV5pPehDkiAECGnq3AU8HjY9cT/UrR5C6rZTNeH7iehb2akmyyEaCC+/PLLSDdB1KMyt4M2w5NAUTi8zldZhWmQiTbDjBQWFUS6iaKZ8JS5ubHLIACW7d2Ew+chxmThfzr25cYugygrdkS4hU2TBAgR5lSVjHbPmx8MDxYL3R5+iPjze0WgdUIIIRoKDYjuaWDgMBMep47FrlBYrCH9DqI+BQIBTAEY37YPEzOGUOp1E2224nO5MAWQ3oc6IgFC1Eir3/8Ox67dZDxwH3E9uke6OUIIISLIaovh2y0+vs7xY7dCtF2h1KnjdMMF/Yxkdo3B5yiUCkqibmkasSYLgVXrUFd+SyAlgai0ZAJ5x1GPFRK4ZCgxw/pQWCq9ELWtQQWIhQsXsnbtWv7xj3+Etu3cuZPZs2ezbds24uPjGT9+PDfddFMEW9k8xfXqSb9X/obBYgnb7nB6ZfE2IYRoRlRVxWazsH5rsIyr0w1Od+W6DOu3+hnax86rK2LQNB2LWScpViO7vwt/ADb9YMFq1rGag/usoS8Ns0nKtjZbAQ3F7UFxeVBdwe/Bn90oLi+Ky43i9qA6PaH7KQYD5odvw7MmJ3iMY4VoxypXpA6szsFy8WBUZ5nMyallDSZALF68mAULFjBgQOWk3MLCQiZMmMDFF1/MzJkz2bJlCzNnziQ+Pp4xY8ZEsLVNn6/Ewe75z9L+xj8S1a4tQJXwcKzIVe2k60nXZJISb6vX9gohhKgfqqri9mi4vdXvd3uDgSLKppJ/QsflAWP5+nEuj8J3P1qqfyCgoPP7LCctEgNs3GnhRIn6i5BR/mXRibJpRNtkQbkGKRA4KQCcFAbcHpTyABDcXhkOVM/p16OqjtIyGb3UCa6qFZgAcHnQXR5UVZUAUcsiHiCOHj3KtGnTyMnJoUOHDmH7/vnPf2I2m5kxYwZGo5FOnTqRm5vLK6+8IgGiDnmLitk+fSbO/bnsLiigzzPzUNTwwt0Op7dKeIDgqtXP/XMzU8b1l54IIYRogjRNw2pRsZqpNkRYzcEhTZcMdOB0g9uroJb3KhhU6NnRg9ur4PaquL0KHq+C26vg8yvoKJhNwVCQV2DgSMGpL1My2nq5INPNgTwj/9lkpVVygIsHunC6FTbstFTTwxF+W3o6asgfqPzk3+VGcXtRnO5gGKguJLg8KN6zDwMAOqBbzeg2C7rNimYzo9us6DYLms0S3G61oNstaFYLRNlIjYsGm6X6EGGzoNgsaM7SX/caiCoiHiC2b99OXFwc77//Pi+88AKHDh0K7du4cSMDBgzAaKxs5uDBg3nppZcoKCggKSkpEk1u0rwnCtn21xm4fj6IKSGejPvuqRIeAIpLPdWuGQHBEFFc6pEAIYQQTZCmabhcHgb2MvJ1TtUyrgN7GXG5PNgsGrZfdDbYrTpDelX/aXFAA49XwWoOBog+XTx0cvpwexTcPgW3R8HjU3B7VDw+hWh78BNll0fB41PxBYK3y9wKuw+c/u/PVdmlJMRo/DvHSkmZesqwkRCjER+jEdBAAar5c9i4+PxVhgFVXvy7qw8Dvqr/xjWhQ3kQsFQJAJrNil4eDrST7qNbzWf9IntcbgxZ/Qh8/m2VfYasfnhcbul9qAMRDxAViw5VJy8vjy5duoRtS01NBeDw4cMSIGqZp6CAbQ/PwH34MOakRHo+OhNb61bV3rfMdfpPF860XwghROPldjkY2if4N3j9Vj9ub7DnYWAvI0P7mCg+hzKuBjUYMCq0SQ0AZ66g06GVj5SEQKiXw2bRGdDNXd7LoYT1cni8wfBREVIKig0UOgynPHbvzh4G9vCw52cTX2+xYTYGh09ZTDpWi1b+Xcdu0enXLXiRmnskgNWskByvYDDUTTeHruvg9YHTjV7mRC9zQ/l387HjmMpcYUOE1PL5A4rv3CoS6YoSfsFf3gsQDAO/DAnBHgPdYqqXxOVwu0jKDpZxDazOCfZE2CwYsvphyB5EUXFRnbehOYp4gDgdt9uN2Rz+KYKlfBy+x3OK8W7inLjz89n+yAzceUexpCTT49GZ2FqmnfL+UTbTaY93pv1CCCEaL7/fT3FRAf26xTC8rx23R8dqUXC5PBQXFeD3n9un1ufCbIJEU+UnzNE2nT5dTjFBA9C0yonaF2S6cLrVKkGj4ue46OBx3d7gA7x+Ba+/IhRUBg+bRaNfNzeBgM7r/xdcnvu+G+zYbbBohYsyl47dqmCzKtitYLMo2G0KNotCeppKi5gAzkIXPocLm9eF0R0MA3qZMxQSKHOjl7nQnS4odcEpXuNTzzAJ0lUlbBiQXt4boFVc+NvCQ0EwDJgb7Ox2v99PQXERMcP6YLl4MLrLg2Kz4HG5KSouqtf3YnPSoAOE1WrF6w3/JVARHOx2eySa1GTl/uNN3HlHsaa1oMejM7CW9/ScSly0hcyMlGqHMWVmpBAXfaZfYUIIIRozv99PqaMQZ5mKqqqUObRGMVTk5A/FUxM0qMHKFb06e8lo5ztpGJUS1sNhKD+mx6eTFAdut445P4+A00XRiWTK/AaKHNVP+L7w+LfEn9jMxoRMvkoZCsRg0nzYAm5sAVf5dzfxvmKyCtajofBDTKfgdt2LzaJjtxkw200QbcetQsBiOmkeQXgPQVMsdeX3+yksdaA6y4ITpp2ljeK92Jg16ACRlpZGfn5+2LaK2y1atIhEk5qsTrffhqKqtPvjOCw1GBoWYzcz6ZrMaqswTb4mU+Y/CCFEM6FpjSM4nDVdR/H4UFxuDC4PZpeH2Ir5AaFJxRUVhdwo//aiuj3cVD4Xw781eJixpjhcBisug638+0lfqpVkd3C4l99gQtE1dEXFp5rwqSZKTDGh5qTavGT/oS1lJjsffBFTpblGI9hVBbMxwPBeTlLiNXYfMOFwqrRL9JEcr+FwKniKldB8D6OhaWWJJvtebIAadIAYMGAAS5cuJRAIYDAEuwrXrl1Lhw4dZP5DLfCVlGCKjQXAaLfT5e47z+rxKfE2pozrL+tACCFEM6aqaqhMZoO9eNP1YPWgk9cWKK8mVFlStGpFIUU/xzKxJiNE2VCibCRH2cAe/LliW9jP9nYQZSPbbCIb8JSXwa34crmDC/RZzGYM3RKgVKNdS09wnyd4X00Ljmgq8euASkUm2HPQxKFjRmLsGsnxGjv2mfl+T+UIAYMaXp0qOT7AoB4evD7YfcCM1ayF9lV8NxkbbuhoFO/FJqJBB4gxY8bw6quvMm3aNG6++Wa+//57Xn/9dWbOnBnppjV6Zfv2s+2vM2n9/66gzZjR53ycGLtZAoMQok7JIqMNk9FoJMZqw2KzojvdKHYrHpcbh9tVt+PONS148V++pkConGjo56rVhBS395zDgG4yVpkTEJpIbLOilU8mVmPsJLZuFQwG5nOfB2i1gNWikBhX/f7YaJUbrqxca0nX9fL51MHAceRoUWjuRvuWPmKjNBJig5OngxPVNdxeBU1TCGgKZW6FMnf5scqPWepSWbvNWu3zq6qO1aRz+TAn8TEam34wU+ZWqymXq2E169gswdBRl4xGI1ZbDDabBbcnWGbY5fLgdjlkDkQdadABIikpiVdffZXZs2czevRoUlJSuP/++xk9+twveAWU7vmJ7dNn4S8t5fg3a2l1xe9QTTLpWQjR8Mgiow2T0WgkKS6ewKp1wVWAT6p8k5Q9iIKaTl7VtNAFvlp+4X/KtQUqqgq5vZzrB+C62RS+tsApqgkFQ0LwPqFV8GrwmigJsefYsnOnKAoWM1jMCgmxYFa00Pzq7h18QGVVxP7dPPTv5kHXwR+gyuRxS/mlgEGFjq2rzvkIaMHg4fQoGI3BuJF7xMTx4lO/Rj07ehjSy8O+w0ZydllonepnSE8PpS6FnfvMYeVzTw4eNZ2qYTQaiYtP4tstPtZvdf6iIlhSvU/qby4aVIB48sknq2w7//zzWbZsWQRa0zQ5ftjN9pmPEihzEpPRhe5/fVjCgxCiwZFFRhu2GKuNwKp1BFaeVHvf5QnW4tchtn83HLsPBFcYPrkn4JfBwHPqaklnollMvygnWn0VoZPXH6hpGGjqFCU4yspk1ImxV+2ZiYvWGNnfVWW7309oTQ67Jfi4Hh29OJxqlSpWFbcrSuaWuVQKHQYSYoK9IyVlKltOtyq5ojN2VClRNp2vN1txexWslmDvR2UpXZ3zOsTy7RZf2Jokbi+h2/26xVDqKDz3F0tUq0EFCFG3irfvYMes2WhuN7Hdu9HtkWkY7bYzP1AIIeqZLDLacKmqisVmDfY8VCOwJgdL9kCiv1gPZVUvQqujWc3lQ4Kq6wmoXFsgLAwYGvuqbo2P0QjRRp1oW2Xo6NL29Os+VYwc69jaR0JsILTSuM2i06Ojpzx0qKEFAz0eBV9AQdeDk70BDh0zUuqq+u9tt0LfnmbWb3VW+9zrt/oZ3teOs0yVORG1TAJEM1H0/VZ2PvYEmsdD3Pm96DbtQQzW6sc3CiFEpMkiow2XqqroTndw2FJ1XB70Mhdap3T8Zc5Tri0Q6jE4h9WHReNRMQzJbtWxWysXskuI0Rh6qlXJy4dXVXQYDenlxun+5aKACnabEadbx32Kjiy3F9wePTSxWtQeCRDNhOvnn9E8HuL7ZtL1wSkYLLJOgxCicZJFRiNL0zSU4Gpo1YcImwUlNgrHZUPkok2cE4MBok7q5Wjfsvo5DKqqEm2Pwmqm2hBhNQcnpJc55H1Y2yRANBMtL/8t5sREEvr3kzkPQohGTRYZjSxN0/C43Biy+gXnPPyCIasfHpdbwoOoc5qm4XJ5GNjLGDYHosLAXkZcLo+8F+uA9Bk2YYWbNuMvLQ3dThoyWMKDEKLRk0VGI8/hdmHIHoThkqHBnggIVmG6ZCiG7EE43DWb+yDEr+V2ORjax8QF/YxYyzsmrWa4oJ+RoX1MuF2OyDawiZIeiCbq2Ner2T1/AdGdOtJj1gyZLC2EaDJkkdHI8/v9FBQXETOsD5aLB6O7PCg2Cx6Xm6KalnAVohb4/X6Kiwro1y2G4X3tuD06VouCy+WREq51SHogmqD8L//N7vkLQNOwp6djsMhCb0KIpmPMmDGUlpYybdo09uzZw4oVK3j99de57bbbIt20ZsXv91NY6iC/oIACt5P8ggIKS2XhLlH//H4/pY5CCo7n4yoroOB4PqWOQnkv1iHpgWhi8j5fyU8LXwJdp8Wlo+j0p1tRpLqFEKIJkUVGGxZN02SMuWgQ5L1YfyRANCFHPvqYvS8vAoKTpjvcMhGlJss4CiFEAyaLjAohRMMiAaKJyPv081B4aPX/rqD9jX+U8CCEEEIIIWqdBIgmIq5XT0wJ8bS4eCRtrx8r4UEIIYQQQtQJCRBNhK11KzKffRpTXFykmyKEEEIIIZowmV3bSOm6Tu6bb1O4aXNom4QHIYQQQghR16QHohHSdZ19ixZz5IMPUc1m+r74PBapfS6EEEIIIeqBBIhGRtc09r78KnmffAZA+4k3SHgQQgghhBD1RgJEI6IHAuxZ+BL5X6wCRaHzX26nxcUjI90sIYQQQgjRjEiAaCT0QIAfF7zAsf98BarKeXf+hdQLR0S6WUIIIYQQopmRANFI5H2+MhQeMu69i+ThwyLdJCGEEEII0QxJgGgk0i4ZhWPXbpKGDCJp8KBIN0cIIYQQQjRTEiAaMM3nQ1FVFIMBxWCgy92TI90kIYQQQgjRzMk6EA1UwONh52NPsOf5v6FrWqSbI4QQQgghBCABokEKuFzsmDWboi3fcfzbtbgOHY50k4QQQgghhABkCFOD43c62TFrNo6duzDYbHSf/jD29DaRbpYQQgghhBCABIgGxV9ayvYZj1H6448YoqLoMeMRYrqcF+lmCSGEEEIIESIBooHwlZSwffosyvbuwxgTTY+Z04nu1DHSzRJCCCGEECKMBIgGomx/Ls4DP2OKi6XHrBlEtW8X6SYJIYQQQghRhQSIBiL+/F50feA+rC1bypwHIYQQQgjRYEmAiCDPseNofj+2lmkAJA4cEOEWCSGEEEIIcXpSxjVC3EePsnXqI2x7eDruo/mRbo4QQgghhBA1IgEiAlxHjrB16l/x5OejmowoqvwzCCGEEEKIxkGGMNUz58GDbHt4Br7CQmxtWtNj1gwsSYmRbpYQQgghhBA1IgGiHpXlHmD7IzPwFRdjb9eWHrOmY46Pj3SzhBBCCCGEqDEJEPWkbH8u2x6ejt/hIKpDB3rM+ium2NhIN0sIIYQQQoizIgGinpgTEzEnJWJNa0GPGY9gjI6OdJOEEEIIIYQ4axIg6okpNoaes6ajGI0Yo6Ii3RwhhBBCCCHOiZT/qUPF27aT9+nnodumuDgJD0IIIYQQolGTHog6UrTlO3bOfhLN68WclEjigP6RbpIQQgghhBC/mgSIOnBiYw67npyD7vOR0K8v8b3Pj3SThBBCCCGEqBUSIGpZwbr1/PC/89D9fhIHDSBjyr2oJlOkmyWEEEIIIUStkABRi45/8y275z2DHgiQNGwIXe65C9UoL7EQQgghhGg65Oq2lpTlHuCHufNB00i58ALOm/wXFIMh0s0SQgghhBCiVkmAqCVR7drS5qo/4C04Qec7/iThQQghhBBCNEkSIH4lPRAIhYW2110Luo6iSnVcIYQQQgjRNMmV7q9w+IMP2fbXmQQ8HgAURZHwIIQQQgghmrRGcbWraRoLFiwgKyuL3r17M3HiRHJzcyPapoMr3mPfq3+nZNt2jn+9OqJtEUIIIYQQor40igCxcOFCli5dymOPPcayZctQFIVbbrkFr9cbkfb8vOwdcl//BwDp/3M1qRePjEg7hBBCCCGEqG8NPkB4vV5ee+01Jk2axIgRI+jatSvz58/n6NGjrFy5sl7bous6uUve4sBbSwFoe/1Y2l53LYqi1Gs7hBBCCCGEiJQGHyB27dpFWVkZgwcPDm2LjY2le/fubNiwod7aoes6+xe/wcF3lgPQ/sY/kn7NVfX2/EIIIYQQQjQEDb4KU15eHgAtW7YM256amsqRI0fqrR3eghPkr/oSgA633ESr3/223p5bCCGEEEKIhqLBBwiXywWA2WwO226xWCguLq63dliSk+gxczple/fSYtTF9fa8QgghhBBCNCQNPkBYrVYgOBei4mcAj8eDzWar17ZEd+pIdKeO9fqcQgghhBBCNCQNfg5ExdCl/Pz8sO35+fmkpaVFoklCCCGEEEI0Ww0+QHTt2pXo6GjWrVsX2lZSUsKOHTvo379/BFsmhBBCCCFE89PghzCZzWbGjRvH3LlzSUxMpHXr1syZM4e0tDRGjRoV6eYJIYQQQgjRrDT4AAEwefJk/H4/Dz/8MG63mwEDBrBo0aIqE6uFEEIIIYQQdatRBAiDwcCUKVOYMmVKpJsihBBCCCFEs9bg50AIIYQQQgghGg4JEEIIIYQQQogakwAhhBBCCCGEqDEJEEIIIRolTdNYsGABWVlZ9O7dm4kTJ5KbmxvpZgkhRJMnAUIIIUSjtHDhQpYuXcpjjz3GsmXLUBSFW265Ba/XG+mmCSFEkyYBQgghRKPj9Xp57bXXmDRpEiNGjKBr167Mnz+fo0ePsnLlykg3TwghmrRGUca1NuXn5xMIBBg5cmSkmyKEEBFx5MgRDAZDpJvxq+zatYuysjIGDx4c2hYbG0v37t3ZsGEDl19++Vkd79f+bdA0DV3Xz+mxomlQFAVVjfznsvJebN5+zfvwbP42NLsAYbFYpHtbCNGsGY3GRr8QZ15eHgAtW7YM256amsqRI0fO+ni/9m9DQ7hwFALkvSjO3dn8bWh2AWLjxo2RboIQQohfyeVyAVT5Y2exWCguLj7r48nfBiGEqDmJqUIIIRodq9UKUKXXwOPxYLPZItEkIYRoNiRACCGEaHQqhi7l5+eHbc/PzyctLS0STRJCiGZDAoQQQohGp2vXrkRHR7Nu3brQtpKSEnbs2EH//v0j2DIhhGj6mt0cCCGEEI2f2Wxm3LhxzJ07l8TERFq3bs2cOXNIS0tj1KhRkW6eEEI0aRIghBBCNEqTJ0/G7/fz8MMP43a7GTBgAIsWLWr0FaaEEKKhU3QpFiyEEEIIIYSoIZkDIYQQQgghhKgxCRBCCCGEEEKIGpMAIYQQQgghhKgxCRBCCCGEEEKIGpMAIYQQQgghhKgxCRBCCCGEEEKIGpMAUQOaprFgwQKysrLo3bs3EydOJDc3N9LNqlULFy5k/PjxYdt27tzJuHHj6NOnDxdeeCGLFi2KUOt+vaKiIv76179ywQUX0LdvX8aOHcvGjRtD+5vKuRYUFDBlyhQGDx5MZmYmt956K3v27Antbyrn+Uv79u0jMzOTFStWhLY1pXM9dOgQGRkZVb7eeecdoGmdq/j1xo8fz4MPPljtvgcffLDK7/pzVZvHOhsHDx4kIyMjbBVy0bTccccdXHPNNVW2jx07loyMDNavXx+2/dNPPyUjI4Mbb7wx7D353nvvkZWVRa9evXjjjTfqvN3NiQSIGli4cCFLly7lscceY9myZSiKwi233ILX641002rF4sWLWbBgQdi2wsJCJkyYQPv27Vm+fDmTJk3i2WefZfny5RFq5a9zzz338N133/H000/z7rvv0qNHD2666SZ++umnJnWut99+Oz///DOvvPIK7777LlarlRtvvBGXy9WkzvNkPp+P++67D6fTGdrW1M71hx9+wGKxsHr1atasWRP6+v3vf9/kzlUIIYYOHcqOHTtwu92hbQ6Hg++//56WLVvy9ddfh91/48aNdOzYkeeee47nnnsutH327NkMHTqUTz/9lNGjR9db+5sDWYn6DLxeL6+99hpTpkxhxIgRAMyfP5+srCxWrlzJ5ZdfHuEWnrujR48ybdo0cnJy6NChQ9i+f/7zn5jNZmbMmIHRaKRTp07k5ubyyiuvMGbMmAi1+Nzk5ubyzTff8Pbbb9O3b18Apk2bxtdff82HH36I1WptEudaWFhImzZtuP322znvvPMA+POf/8yVV17Jjz/+yNq1a5vEef7Sc889R1RUVNi2pvT+Bdi9ezcdOnQgNTW1yr7XX3+9SZ2rEEIMGTIEn8/H1q1bGTBgAADffvstsbGxXH311Xz++efcd999oftv2LCBYcOGERMTE3ackpISBg4cSOvWreu1/c2B9ECcwa5duygrK2Pw4MGhbbGxsXTv3p0NGzZEsGW/3vbt24mLi+P999+nd+/eYfs2btzIgAEDMBorM+bgwYPZt28fBQUF9d3UXyUhIYGXX36Znj17hrYpioKu6xQXFzeZc01ISODpp58OhYfjx4+zaNEi0tLS6Ny5c5M5z5Nt2LCBZcuW8dRTT4Vtb2rn+sMPP9C5c+dq9zW1cxX1JyMjI2zYH0B2dnboE9xAIMCcOXMYMWIEPXv25De/+Q1vv/122P39fj+PPfYY/fr1Y/DgwTz99NP4/f7Q/pycHCZMmEC/fv3o2bMnv/vd7/jwww9D+x988EGmTJnCU089xZAhQ+jduzd//vOfOXbsWOg+u3fv5o9//CN9+vTh0ksv5b///W9dvByiAenYsSNpaWls2rQptG316tUMHTqUrKwsdu3aRX5+PhAMCbt372b48OGhYXUVw9wApk6dGvrZ6/UyZ84csrKyyMzM5JprrmHNmjX1f4JNgASIM8jLywOgZcuWYdtTU1M5cuRIJJpUa7Kzs5k3bx7p6elV9uXl5ZGWlha2reLTz8OHD9dL+2pLbGwsI0aMwGw2h7Z98sknHDhwgOHDhzepc63wyCOPMGzYMD799FNmz56N3W5vcudZUlLC/fffz8MPP1zl/2dTO9fdu3dTUFDAddddx9ChQxk7diyrV68Gmt65iobjrbfe4tNPP2X+/Pl89tlnjBs3jhkzZoTNH9u0aRPHjx9n6dKlPPHEEyxfvpwnn3wSCPZyT5w4ka5du7JixQr+7//+j169evHQQw9x/Pjx0DE++eQTioqKWLJkCc8//zw5OTnMnz8fCA5bufHGG4mOjuadd97hr3/9KwsXLqzfF0JExJAhQ9i8eXPo9po1a8jKyqJnz57Ex8eHfgfm5ORgMBgYOHBg6L4tW7YMBYOpU6eGfn7ooYdYvXo1c+bM4V//+heXXXYZf/rTn/jPf/5TfyfWREiAOAOXywUQdvEJYLFY8Hg8kWhSvXC73dWeM9DozzsnJ4epU6cycuRIsrOzm+S53nDDDSxfvpwrrriCO+64g+3btze585wxYwZ9+vTh97//fZV9TelcvV4v+/fvp7S0lLvuuouXX36ZXr16ccstt7B27domda6i9nzwwQdkZmZW+frggw9qfIwDBw5gt9tJT0+ndevWjBs3jr///e9hQ15TUlJ46qmnOO+887jooou48847Wbp0KS6XC6/Xy1/+8hemTJlCu3bt6NSpE7fddlvoPV0hOjqaWbNm0alTJ7KysrjyyivJyckB4KOPPsLlcoWeY9iwYUydOrXWXifRcFUECF3X2bNnD3l5eQwbNgxVVRkyZEgoQGzYsIHMzEzsdnvosQaDgZSUFABiYmJISUkhNzeXDz/8kNmzZzN48GDat2/PhAkTuPzyy6XwxDmQORBnYLVageAf8YqfIfiH2WazRapZdc5qtVaZJF5xMXLyf9LG5osvvuC+++6jd+/ePP3000DTPNeK4S6PPvooW7ZsYcmSJU3qPN977z02btx4youhpnSuZrOZDRs2YDQaQ0GhZ8+e/PTTTyxatKhJnauoPdnZ2WFjxCvMnTuXoqKiGh3j+uuv54svvuCCCy6gZ8+eDBs2jMsuu4ykpKTQfXr27BkKrADnn38+Pp+P/fv3061bN8aMGcOSJUvYs2cP+/fvZ+fOnUBweFSFdu3aYTKZQrdjYmLw+XxAsPetffv2YWPbMzMza/YiiEZtyJAhFBUVsXfvXtasWUPXrl1DoWD48OGhv+EbN24kOzv7jMfbsWMHAH/84x/Dtvt8PmJjY2u59U2fBIgzqBgakZ+fT9u2bUPb8/Pz6dq1a6SaVefS0tJC4wsrVNxu0aJFJJr0qy1ZsoTZs2czatQo5s6dG7oYayrnWlBQwNq1a7nsssswGAwAqKpKp06dyM/PbzLnCbB8+XIKCgq48MILw7ZPnz6dRYsW0apVqyZzrlB9EOjSpQtr1qxpUv+uovZERUXRrl27arefHCB0XQ/bX3HhDtC+fXs+//xz1q9fzzfffMOqVat48cUXeeKJJ0IVbSp+11TQNA0IBt+ffvqJsWPH0r17d4YNG8bIkSNJSEjg6quvDnvML3vQfumXbTx5vo9oulJTU+ncuTObN29mzZo1DB8+PLRv+PDhTJs2jW3btrFjxw4eeeSRMx6v4n305ptvVim8oaoyIOdsySt2Bl27diU6Ojqs3nRJSQk7duygf//+EWxZ3RowYAA5OTlhnxKtXbuWDh06hH361Fi89dZbPProo1x//fU888wzYX+wmsq55ufnc++994bVx/b5fOzYsYNOnTo1mfOE4KeoH3/8Me+9917oC2Dy5Mm8/PLLTepcd+3aRWZmZti4c4Bt27bRuXPnJnWuon6ZTCYcDkfodmlpKSdOnAjdfuONN/j8888ZNmwY999/Px988AFDhgzh448/Dt1n586dodAAwSGiVquV9PR03n77bZKSkli8eDG33HILI0aMCM19+GUoOJVu3bqxb9++sHZt3br1nM9ZNC4Vw5hycnLCAkRFcZClS5cSFRVFjx49znisigIj+fn5tGvXLvS1YsUKKXt9DiRAnIHZbGbcuHHMnTuXVatWsWvXLu6++27S0tIYNWpUpJtXZ8aMGUNpaSnTpk1jz549rFixgtdff53bbrst0k07a/v27ePxxx9n1KhR3HbbbRQUFHDs2DGOHTuGw+FoMufatWtXhg8fzsyZM9m4cSO7d+/mgQceoKSkhBtvvLHJnCcEP1k/+Q9AxSetSUlJtG7dukmda5cuXTjvvPNC/64//fQTTzzxBFu2bOFPf/pTkzpXUb8yMzNZtmwZ27dvZ/fu3dx///1hn+4XFBQwa9YsVq1axaFDh/j666/ZsWNH2BCiI0eOMHXqVH788Uc+++wznnvuOW6++WbMZjNpaWnk5eXx1VdfcejQIT7//HNmzJgBUON1lC6//HKSkpK499572bVrF+vXr+fxxx+v1ddBNFwnB9aKMuwVhg8fzkcffcTQoUNr1INQMU9n+vTprFq1ip9//plFixbx0ksvVVtMRpye9APWwOTJk/H7/Tz88MO43W4GDBjAokWLztjt2pglJSXx6quvMnv2bEaPHk1KSgr3339/o1yI5bPPPsPn87Fy5UpWrlwZtm/06NE8+eSTTeJcFUXhmWeeYd68edx11104HA769+/Pm2++SatWrQCaxHnWRFN6/6qqyosvvsjcuXO56667KCkpoXv37vz9738PlSZsKucq6teMGTOYOXMm1157LYmJiUyYMCFsQca//OUv+P1+Hn30UY4fP05KSgrXXXddWDgdOXIkBoOBa665BpvNxtixY/nzn/8MBMea7927l/vvvx+v10v79u255557WLBgAd9//z0XXHDBGdtot9t54403mDVrFmPHjiUuLo4777zzlCtti6Zl0KBBeL1esrKyqlxzDR8+nMWLFzNs2LAaH2/+/PnMnz+f6dOnU1xcTHp6Oo8++qismXMOFL2m/YhCCCGEEEKIZk+GMAkhhBBCCCFqTAKEEEIIIYQQosYkQAghhBBCCCFqTAKEEEIIIYQQosYkQAghhBBCCCFqTAKEEEIIIYQQosYkQAghhBBCCCFqTAKEEEIIIYQQosYkQIhm6cEHHyQ7O/uU+8ePH8/48eNr7fkOHjxIRkYGK1asqLVjnovs7GxZwVUIIYQQv4oECCGEEEIIIUSNSYAQQgghhBBC1JgECCHOoLrhTOvWrSMjI4N169YBoGkazz77LNnZ2fTs2ZPs7GyefvppfD5f2OOOHj3Kbbfdxvnnn8+IESNYsGABgUAgtN/tdjNv3jwuueQSevbsSd++fZkwYQI7d+4M3efBBx/kxhtvZPny5Vx66aX07NmTK664gq+++irsuXbt2sWECRPIzMzkoosu4v3336/tl0YIIYQQzZAx0g0QIpL8fn+123VdR1GUGh/nlVde4c033+SBBx4gPT2d7777jvnz52MymZg0aVLofs899xxXXnklL7zwAps3b+bFF18kEAhw9913A3D//fezYcMG7r33Xtq2bcv+/ft59tlnufvuu/nkk09Cbdq2bRv5+flMnjyZ6Ohonn32WSZPnszXX39NXFwcR48eZdy4cbRt25Y5c+ZQWlrK3LlzKSgo+BWvlhBCCCGEBAjRjB06dIgePXqccv/AgQNrfKz169fTo0cPxowZE3qszWYjOjo67H5DhgzhiSeeACArK4vS0lLeeOMNJk6ciM1mo6ysjEceeYTf/va3oeOUlZXx5JNPcuzYMVJTUwFwOBysWLGCtm3bAmC32xk3bhz//e9/ufTSS1m8eDF+v59XXnmFpKQkADp06MA111xT43MSQgghhKiOBAjRbKWkpPC3v/2t2n3Tp08/q2MNGjSIefPmcd111zFq1CguuOACxo0bV+V+FcGgwiWXXMLrr7/Oli1bGDFiBIsWLQIgPz+f3Nxc9u7dy7///W+AsOFQiYmJofAAkJaWBoDL5QIgJyeHPn36hMIDQO/evWnVqtVZnZcQQgghxC9JgBDNltlsplevXtXui4qKOqtj3XzzzURFRbF8+XKeeuopnnzySbp06cLUqVMZMmRI6H7Jyclhj0tMTASguLgYgNWrV/P444+zd+9eoqKiyMjICLVF1/XQ42w2W9hxKoY2aZoWOl6bNm2qtDMlJeWszksIIYQQ4pdkErUQNXDyRGcAp9MZdltVVa6//npWrFjBN998wxNPPIHH42HSpEl4vd7Q/UpKSsIed/z4cQCSkpI4cOAAd9xxB127dmXlypVs2rSJt99+m4suuuis25uQkBA69smKiorO+lhCCCGEECeTACHEGURHR5OXlxe2bdOmTWG3r732Wh577DEgGAb+8Ic/cP311+NwOCgtLQ3db/Xq1WGP++ijj7DZbPTu3Ztt27bh8Xi47bbbwoYnVTzm5B6IMxk8eDCbN2/m6NGjoW179uzh559/rvExhBBCCCGqI0OYhDiDiy66iC+//JLZs2dz8cUXk5OTw3vvvRd2nwEDBvDaa6+RnJxMZmYmR48e5e9//zsDBw4kMTEx1GPx+eef06JFC4YOHcqaNWtYtmwZd955J9HR0fTo0QOj0cicOXOYOHEiXq+XFStW8J///Aeo2utxOjfccAPvvvsuN910E5MmTSIQCPDMM89gMplq62URQgghRDMlAUKIMxgzZgwHDhzgX//6F8uWLWPgwIE8++yzjB07NnSfO++8E7PZzPLly3nhhReIiYkhOzube++9N+xYDz74IJ9++imLFy8mJSWFhx56iBtuuAGAdu3aMW/ePJ5//nluv/124uLi6NOnD//4xz8YP348GzduJCMjo0ZtTkhI4O2332b27Nk8+OCDREVFcfPNN/Pxxx/X3gsjhBBCiGZJ0c9mXIQQQgghhBCiWZM5EEIIIYQQQogakwAhhBBCCCGEqDEJEEIIIYQQQogakwAhhBBCCCGEqDEJEEIIIYQQQogakwAhhBBCCCGEqDEJEEIIIYQQQogakwAhhBBCCCGEqDEJEEIIIYQQQogakwAhhBBCCCGEqDEJEEIIIYQQQoga+/8NwIL3EtQ1KwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plotting code - don't get sidetracked by this, it's not that important\n", "plt.figure(figsize=(8,4))\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, errorbar=None, color=[0.9,0.9,0.9])\n", "sns.lineplot(data=socks.T, legend=False, marker='o') # ci=None switches off errorbars\n", "plt.ylabel('pairs of socks owned')\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "233aae63", "metadata": {}, "source": [ "**Note-** \n", "\n", "* In most couples (11/14) the husband owns more socks that the wife\n", "* There is one couple in which the husband has an extreme number of socks.\n", "\n", "#### Test statistic\n", "\n", "We need to sumamrize the sex difference as a number. \n", "\n", "Here, we are interested in the within-pairs difference in sock ownership, so we need the mean difference in pairs of socks owned, for [husband - wife].\n", "\n", "This mean difference, mean(Husband-Wife), is calculated as follows:" ] }, { "cell_type": "code", "execution_count": 4, "id": "0c6ac6e6", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "4.571428571428571" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "socks['difference'] = socks.Husband-socks.Wife # create a new column for the within-couple difference\n", "socks.difference.mean()" ] }, { "cell_type": "markdown", "id": "24f48262", "metadata": {}, "source": [ "On average, the husbands own 4.6 more pairs of socks than their wives.\n", "\n", "### Is the result statistically significant? \n", "\n", "In our particular sample, husbands have (on average) more socks than their wives - our test statistic (mean(Husband-Wife)) is 4.6. \n", "\n", "However, even if there were no true sex differences in sock ownership, if we picked 10 random couples we wouldn't expect the average difference in number of socks owned to be *exactly zero* - just by chance the wives would sometimes have more socks than the husbands, or vice versa. \n", "\n", "Is it plausible that the sex difference we observed (husbands have 4.6 more pairs of socks than their wives) arose due to chance from a situation in which there is no true sex difference in sock ownership (ie, the *null hypothesis is true*?)\n", "\n", "In that case, it would be just chance whether the person with more socks in each couple was the man or the woman.\n", "\n", "\n", "### Permutation\n", "\n", "To work out how likely it would be to get a test statistic (mean difference in socks) as large as 4.6, we are going to do is shuffle the data around to create many new (re)samples in which, **within each couple**, we we randomly decide who gets which label, *husband* or *wife*. \n", "\n", "Within each resample, some couples will retain their true labels and some will be relabelled; which couples are which will differ between resamples.\n", "\n", "For each resample, we then calculate the mean difference, mean(Husband-Wife).\n", "\n", "After repeating the process thousands of times, we can ask on what proportion of random resamples we get a values of mean(Husband-Wife) at least as large as our observed difference, 4.6.\n", "\n", "An important consideration in the shuffling process is that, whilst we want to shuffle the data as if sex was irrelevant (randomly allocate the labels 'husband' and 'wife'), we want to **preserving the non-sex-related features of the dataset**:\n", "\n", "**Non-sex-related (to be retained in shuffled datasets):**\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 is reasonable\n", " * ie the distriution of differences has positive skew\n", " \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", "**Sex effects (will be 'shuffled out'):**\n", "\n", "* in most couples the husband has more socks\n", "* the outlier individuals with loads of socks tend to be male\n", "\n", "To achieve shuffles that **preserve the non-sex-related information** whilst **ignoring the sex of the sock owner**, we only move people within pairs (exchange the labels 'husband' and 'wife'). We *never* move people between pairs (swap someone from one couple into another).\n", "\n", "\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", "* We assume there are no sex differences, so we can swap men for women\n", "* We assume that it *does* matters which couple you are in (high-sock men tend to have high-sock wives), therefore we cannot swap people between couples.\n", "\n", "Therefore, the only shuffling that we are allowed is to swap the labels 'Husband' and 'Wife' within couples. \n", "\n", "Conceptually, to generate a new simulated dataset (one *resample*), we go through each couple in turn and 'flip a virtual coin' to decide whether the labels 'Husband' and 'Wife' are flipped.\n", "\n", "### Visualizing random shuffles in the paired data\n", "\n", "Here are the original data and four random permutations:\n", "\n", "\n", "\n", "Look closely. Each coloured line represents one couple. Note that some of the lines are left-right flipped, but the lines never move up or down overall (because we are only flipping people within couples).\n", "\n", "The grey bars show the group means in each permutation. The values of the mean difference (Husband-Wife) for each permutation are shown in red text.\n", "* Sometimes the mean number of socks is higher for the husbands (positive differences), sometimes for the wives (negative differences). \n", " * These two situations should happen equally often in the permuted data as the 'high sock' partner in each couple is assigned equally often the labels 'husband' and wife'.\n", "* The group with the higher mean does tend to be the group with the extreme 48-sock individual in it.\n", "\n" ] }, { "cell_type": "markdown", "id": "ae1e1fe3", "metadata": {}, "source": [ "## Run 10,000 shuffles\n", "\n", "To get a fine-grained estimate of how likely it is to get a test statistic of 4.6 due to change, we run a *lot* of shuffles - say 10,000.\n", "\n", "We can do this using a function from the library `scipy.stats`, called `scipy.stats.permutation_test()`\n", "\n", "First of all we will run this (using the code block below) and learn about the output. Afterwards we will have a closer look at the Python syntax including defining a function." ] }, { "cell_type": "code", "execution_count": 5, "id": "200125b4", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PermutationTestResult(statistic=4.571428571428571, pvalue=0.0072992700729927005, null_distribution=array([-0.71428571, -1.28571429, 3.28571429, ..., 3. ,\n", " 2.42857143, 2. ]))\n" ] } ], "source": [ "def mDiff(x, y):\n", " return np.mean(x-y)\n", "\n", "results = stats.permutation_test((socks.Husband, socks.Wife), mDiff, permutation_type='samples', alternative='greater', n_resamples=10000)\n", "print(results)" ] }, { "cell_type": "markdown", "id": "a13b89d1", "metadata": {}, "source": [ "The results are displayed below the code box; there are number of outputs from the test (contained in the parentheses)\n", "\n", "* **statistic** - the test statistic (mean difference) from our actual data, ie 4.6\n", "* **pValue** - the proportion of shuffles in which a result greater than 4.6 was obtained (due to chance)\n", "* **null distribution** - this contains all 10,000 mean-difference values, for the 10,000 shuffles\n", "\n", "### The Null Distribution\n", "\n", "The **null distribution** is the set of 10,000 values we get for mean difference (Husband-Wife), from our 10,000 shuffles (the equivalent of the red numbers on the 4-permutation example above)\n", "\n", "Let's plot the null distribution in a histogram:" ] }, { "cell_type": "code", "execution_count": 6, "id": "ea965aee", "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAGgCAYAAABFdAY8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAs2klEQVR4nO3df1CWdb7/8dcNePMjJJGAm9xSw4yoxIq7pVOEo4vnHD11Is5+nQrW1Fj7sdCPDZqEEifJNtlI83AqfpRHayHDsbLTbEQ7U9twCCjrtEhuHmIrQJRINn7dInz/cLk7BK2AwHXfXM/HDDPcnx/4vhgnX12fz/W5LAMDAwMCAAAwCQ+jCwAAAJhKhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqXkYX4Gqio6PlcDgUHBxsdCkAAGCUjh49KqvVqpqamtOOJfz8QG9vr06ePGl0GQAAYAz6+vo02nObCT8/EBISIkmqqKgwuBIAADBay5YtG/VY9vwAAABTIfwAAABTIfwAAABTIfwAAABTIfwAAABTIfwAAABTIfwAAABTIfwAAABTIfwAAABTIfwAAABTIfwAAABTIfwAAABTIfwAAABTIfwAAABTIfwAAABTIfwAAABTIfwALqK/f8CQuQBgNl5GFwDgFA8Pi0rLD+loe9eY5gUH+mlV/MJJqgoAph/CD+BCjrZ3qelYp9FlAMC0xrIXAAAwFcIPAAAwFcIPAAAwFUPDT1VVlS666KIRv5YtWyZJOnjwoJKSkrR48WItWbJERUVFQ35Gf3+/tm/frtjYWEVFRWnt2rVqbGw04nIAAIAbMDT8XH755frjH/845Ku4uFheXl6644471N7erjVr1mjevHkqKytTamqqtm3bprKyMufPyM/PV0lJiTZv3qzS0lJZLBalpKTI4XAYeGUAAMBVGfq0l9VqVXBwsPPziRMntGXLFi1fvlw///nP9eyzz8pqtSo7O1teXl4KDw9XY2OjCgoKlJiYKIfDoeLiYqWnpysuLk6SlJeXp9jYWJWXl2vlypVGXRoAAHBRLrXn58UXX1Rzc7MeeughSVJNTY3sdru8vL7PaDExMWpoaFBbW5vq6+vV2dmpmJgYZ39AQIAiIyNVXV095fUDAADX5zLhp7e3V88884xWr16tkJAQSVJLS4tsNtuQcYN9TU1NamlpkSSFhYUNG9Pc3DwFVQMAAHfjMuHn1VdfVW9vr5KTk51tPT09slqtQ8Z5e3tLOhWWuru7JWnEMb29vZNcMQAAcEcuE3727dun5cuXKzAw0Nnm4+MzbOPyYKjx8/OTj4+PJI04xtfXd5IrBgAA7sglws8333yjjz76SCtWrBjSbrPZ1NraOqRt8HNoaKhzuWukMT9cLgMAAJBcJPx8+OGHslgsuuqqq4a02+121dbW6uTJk862yspKzZ8/X0FBQYqIiJC/v7+qqqqc/R0dHaqrq1N0dPSU1Q8AANyHS4Sf+vp6nXfeecOWqhITE/Xdd98pMzNTn3/+ufbu3audO3dq/fr1kk7t9UlKSlJubq4qKipUX1+v++67TzabTfHx8UZcCgAAcHEu8Vb3Y8eOadasWcPag4KCVFhYqJycHCUkJCg4OFgZGRlKSEhwjklLS1NfX5+ysrLU09Mju92uoqKiYZugAQAAJMkyMDAwYHQRrmTwtRoVFRUGVwIz2vHyATUd6xzTnHPPOUu/+n+LJ6cgAHATY/n32yWWvQAAAKYK4QcAAJgK4QcAAJgK4QcAAJgK4QcAAJgK4QcAMKL+/vE/DHwmc4HJ5hLn/AAAXI+Hh0Wl5Yd0tL1rTPOCA/20Kn7hJFUFnDnCDwDgRx1t7xrz2VOAq2PZCwAAmArhBwAAmArhBwAAmArhBzA5nugBYDZseAZMjid6AJgN4QcAT/QAMBWWvQAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAAgKkQfgAALqW/f8CQuTAPL6MLAADg//LwsKi0/JCOtneNaV5woJ9WxS+cpKownRB+AAAu52h7l5qOdRpdBqYplr0AAICpEH4AAICpuET42bdvn1asWKHLLrtMK1eu1JtvvunsO3jwoJKSkrR48WItWbJERUVFQ+b29/dr+/btio2NVVRUlNauXavGxsapvgQAAOAmDA8/r776qjZs2KBVq1Zp//79WrFihe6//3599NFHam9v15o1azRv3jyVlZUpNTVV27ZtU1lZmXN+fn6+SkpKtHnzZpWWlspisSglJUUOh8PAqwIAAK7K0A3PAwMD2rZtm1avXq3Vq1dLku6++259+OGH+uCDD/TBBx/IarUqOztbXl5eCg8PV2NjowoKCpSYmCiHw6Hi4mKlp6crLi5OkpSXl6fY2FiVl5dr5cqVRl4eAABwQYbe+fnf//1fff3117r++uuHtBcVFWn9+vWqqamR3W6Xl9f3GS0mJkYNDQ1qa2tTfX29Ojs7FRMT4+wPCAhQZGSkqqurp+w6AACA+zA0/HzxxReSpK6uLq1bt05XX321fv7zn+udd96RJLW0tMhmsw2ZExISIklqampSS0uLJCksLGzYmObm5kmuHgAAuCNDw893330nSXrwwQf1L//yLyouLtY111yju+66S5WVlerp6ZHVah0yx9vbW5LU29ur7u5uSRpxTG9v7xRcAQAAcDeG7vmZMWOGJGndunVKSEiQJF188cWqq6vT888/Lx8fn2EblwdDjZ+fn3x8fCRJDofD+f3gGF9f36m4BAAA4GYMvfMzuKS1cOHQ48gXLFigr776SjabTa2trUP6Bj+HhoY6l7tGGvPD5TIAAADJ4PATGRmps846Sx9//PGQ9kOHDun888+X3W5XbW2tTp486eyrrKzU/PnzFRQUpIiICPn7+6uqqsrZ39HRobq6OkVHR0/ZdQAAAPdh6LKXj4+Pbr/9dv37v/+7QkNDtWjRIr3xxht6//339cILL2jBggUqLCxUZmambr/9dn3yySfauXOnNm3aJOnUXp+kpCTl5uZq9uzZmjNnjrZu3Sqbzab4+HgjLw0AALgow19setddd8nX11d5eXk6cuSIwsPD9fTTT+unP/2pJKmwsFA5OTlKSEhQcHCwMjIynPuDJCktLU19fX3KyspST0+P7Ha7ioqKhm2CBgAAkFwg/EjSmjVrtGbNmhH7Fi1apNLS0h+d6+npqfT0dKWnp09WeQAAYBox/PUWAAAAU4nwAwAATIXwAwAurr9/wJC5wHTlEnt+AAA/zsPDotLyQzra3jWmecGBfloVv/D0AwGTIfwAgBs42t6lpmOdRpcBTAssewEAAFMh/AAAAFMh/AAAAFMh/AAAAFMh/AAAAFMh/AAAAFMh/AAAAFMh/AA/wGm6ADC9ccgh8AOcpgsA0xvhBxgBp+kCwPTFshcAADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVw8PP119/rYsuumjY1549eyRJBw8eVFJSkhYvXqwlS5aoqKhoyPz+/n5t375dsbGxioqK0tq1a9XY2GjEpQAAADfgZXQBn332mby9vfX222/LYrE422fOnKn29natWbNGP/vZz7Rp0yYdOHBAmzZt0qxZs5SYmChJys/PV0lJibZs2aLQ0FBt3bpVKSkp2r9/v6xWq1GXBQAAXJTh4efQoUOaP3++QkJChvXt3LlTVqtV2dnZ8vLyUnh4uBobG1VQUKDExEQ5HA4VFxcrPT1dcXFxkqS8vDzFxsaqvLxcK1eunOrLAQAALs7wZa/PPvtMCxYsGLGvpqZGdrtdXl7fZ7SYmBg1NDSora1N9fX16uzsVExMjLM/ICBAkZGRqq6unvTaAQCA+zE8/Bw6dEhtbW265ZZb9A//8A+6+eab9d5770mSWlpaZLPZhowfvEPU1NSklpYWSVJYWNiwMc3NzVNQPQAAcDeGLns5HA598cUX8vX1VUZGhvz8/PTaa68pJSVFzz//vHp6eobt2/H29pYk9fb2qru7W5JGHHP8+PGpuQgAAOBWDA0/VqtV1dXV8vLycgaYSy+9VIcPH1ZRUZF8fHzkcDiGzOnt7ZUk+fn5ycfHR9KpEDX4/eAYX1/fKboKAADgTgxf9vLz8xt252bhwoU6cuSIbDabWltbh/QNfg4NDXUud4005ofLZQAAAJLB4ae+vl6XX365ampqhrR/+umnWrBggex2u2pra3Xy5ElnX2VlpebPn6+goCBFRETI399fVVVVzv6Ojg7V1dUpOjp6yq4DAAC4D0PDz8KFC3XhhRdq06ZNqqmp0eHDh7VlyxYdOHBAd9xxhxITE/Xdd98pMzNTn3/+ufbu3audO3dq/fr1kk4tmyUlJSk3N1cVFRWqr6/XfffdJ5vNpvj4eCMvDQAAuChD9/x4eHjomWeeUW5uru699151dHQoMjJSzz//vC666CJJUmFhoXJycpSQkKDg4GBlZGQoISHB+TPS0tLU19enrKws9fT0yG63q6ioiAMOYRr+fjPU3z8gDw/L6QcDAIw/5HD27Nl67LHHfrR/0aJFKi0t/dF+T09PpaenKz09fTLKA1yer9VLHh4WlZYf0tH2rjHNXXh+oJbHzJ2kygDANRkefgBMjKPtXWo61jmmOcGzeCoSgPkY/rQXAADAVCL8AAAAUyH8AAAAUyH8AAAAUyH8AAAAUyH8AAAAUyH8AAAAUyH8AAAAUyH8AAAm1OArVwBXxQnPAIAJxStX4OoIPwCAScErV+CqWPYCAACmQvgBAACmQvgBAACmQvgBAACmQvgBgGmKR86BkfG0FwBMUzxyDoyM8AMA0xyPnANDsewFAABMhfADAABMhfADAABMhfADAABMhfADAABMhfADAABMhfADAABMhfADAABMhfADTBBeJQAA7oETnoEJwqsEAMA9EH6ACcarBAD31d8/IA8Py5TPxdQi/AAA8DfjvXsbHOinVfELJ6kqTLRxhZ/q6mpFRkbqrLPOGtbX0dGh9957TytXrjzj4gAAmGrjuXsL9zKuDc+/+MUvdPjw4RH76urq9NBDD51RUQAAAJNl1Hd+HnzwQTU3N0uSBgYGlJ2dLX9//2HjvvjiC51zzjkTVyEAAMAEGvWdn3/8x3/UwMCABga+f5R38PPgl4eHhxYvXqwtW7ZMSrEAAABnatR3fpYuXaqlS5dKkpKTk5Wdna3w8PBJKwzA9MZTNQCMMq4Nz7t27ZroOiRJDQ0Nuummm/Twww/rpptukiQdPHhQOTk5+vTTTzVr1iwlJydr3bp1zjn9/f3asWOH9uzZo46ODl155ZXauHGj5s7lzBTAlfFUDQCjjCv8dHd365lnntEf/vAHdXd3q7+/f0i/xWLR22+/PaafeeLECT3wwAPq6vr+P4Tt7e1as2aNfvazn2nTpk06cOCANm3apFmzZikxMVGSlJ+fr5KSEm3ZskWhoaHaunWrUlJStH//flmt1vFcHoApwlM1AIwwrvCTk5OjsrIyXXXVVbr44ovl4XHmb8l4+umnhz06//LLL8tqtSo7O1teXl4KDw9XY2OjCgoKlJiYKIfDoeLiYqWnpysuLk6SlJeXp9jYWJWXl/O4PQAAGGZc4eett97Sfffdp1/+8pcTUkR1dbVKS0u1b98+LVmyxNleU1Mju90uL6/vy4yJidGzzz6rtrY2ff311+rs7FRMTIyzPyAgQJGRkaqurib8AACAYcZ1y6avr0+LFi2akAI6OjqUkZGhrKwshYWFDelraWmRzWYb0hYSEiJJampqUktLiyQNmxcSEuJ8LB8AAOD/Glf4ufbaa/Xuu+9OSAHZ2dlavHixrr/++mF9PT09w/bteHt7S5J6e3vV3d0tSSOO6e3tnZD6AADA9DKuZa8VK1Zo48aN+uabbxQVFSVf3+EvZbzxxhtP+3P27dunmpoavf766yP2+/j4yOFwDGkbDDV+fn7y8fGRJDkcDuf3g2NGqgkAAGBc4efee++VdCq87Nu3b1i/xWIZVfgpKytTW1vbkH0+krRx40YVFRXp3HPPVWtr65C+wc+hoaHq6+tztp1//vlDxkRERIz+ggAAgGmMK/xUVFRMyB+em5urnp6eIW3Lly9XWlqaVqxYoTfeeEMlJSU6efKkPD09JUmVlZWaP3++goKCNHPmTPn7+6uqqsoZfjo6OlRXV6ekpKQJqREAAEwv4wo/c+bMmZA/PDQ0dMT2oKAgzZkzR4mJiSosLFRmZqZuv/12ffLJJ9q5c6c2bdok6dRen6SkJOXm5mr27NmaM2eOtm7dKpvNpvj4+AmpEQAGcSo1MD2MK/zs2LHjtGN+9atfjedHDxEUFKTCwkLl5OQoISFBwcHBysjIUEJCgnNMWlqa+vr6lJWVpZ6eHtntdhUVFXHAIYAJx6nUwPQw4eHH399fISEh4w4/n3322ZDPixYtUmlp6Y+O9/T0VHp6utLT08f15wHAWHAqNeD+xhV+6uvrh7V1dXWptrZW2dnZevjhh8+4MAAAgMlw5u+l+Bs/Pz/Fxsbq7rvv1hNPPDFRPxYAAGBCTVj4GRQWFqbDhw9P9I8FAEmSv98M9fcPjHv+mcwFMD2Ma9lrJAMDA2publZBQcGEPQ0GAD/ka/Vi4zGAMzKu8BMRESGLZeRHNgcGBlj2AjDp2HgMYLzGFX7uvvvuEcOPv7+/lixZonnz5p1pXQAAAJNiXOEnNTV1ousAAACYEuPe8+NwOLR3715VVVWpo6NDgYGBio6OVkJCgvPN6wAAAK5mXOGno6NDv/jFL1RfX69zzz1XwcHBamho0P79+/Xiiy/qpZde0syZMye6VgAAgDM2rkfdf/vb36qlpUW7d+/WO++8o9LSUr3zzjvavXu32tratG3btomuEwAAYEKMK/xUVFTo3nvvVXR09JD26OhopaWl6a233pqQ4gAAACbauMJPZ2enzjvvvBH7zjvvPH377bdnUhMAAGN2pgdgwjzGtefnggsu0B/+8Addc801w/oqKio0d+7cMy4MAICxOJMDMCVp4fmBWh7Dv19mMK7ws27dOt1///1yOBy6/vrrdc455+jYsWN6/fXXtWfPHmVnZ09wmQAAjM54D8AMnuU7CdXAFY0r/KxYsUJffPGFnnnmGe3Zs8fZPmPGDN19991atWrVhBUIAAAwkcYVfrq6unTXXXcpKSlJBw4c0PHjx9Xc3KxVq1bp7LPPnugaAQAAJsyYNjwfPHhQN954o1544QVJUkBAgK677jpdd911euqpp3TLLbfwRncAAODSRh1+vvzyS9122206fvy4FixYMKTParVqw4YN6uzs1C233KKWlpYJLxSAa+HJGgDuatTLXs8995wCAwNVUlKiWbNmDenz9fVVUlKS/vmf/1n/9m//pmeeeYZNz8A0dyZP1vBUDQAjjTr8VFZW6o477hgWfP6voKAgrVmzRi+++OJE1AbADYznyRqeqgFgpFEvex09enRU5/csXLiQZS8AAOCyRh1+Zs+erdbW1tOO++abb/7u3SEAAAAjjTr82O127d2797Tj9u3bp4svvviMigIAAJgsow4/ycnJqqqq0uOPP67e3t5h/Q6HQ7/5zW/03nvv6dZbb53QIgEAACbKqDc8X3bZZXrooYf02GOP6dVXX9XVV1+tn/zkJzp58qSamppUVVWl9vZ23XPPPYqNjZ3MmgEAAMZtTCc833rrrYqIiFBRUZEqKiqcd4DOOussXXvttVq7dq2ioqImpVAAAICJMObXW1x55ZW68sorJUnt7e3y8PDglRYAAMBtjOvdXoMCAwMnqg5gQvX3D8jDw2J0GQAAF3RG4QdwVZw8DAD4MYQfTFucPAwAGMmY3uoOAADg7gg/AADAVAg/ADDJ/P1mqL9/wOgyAPwNe34AYJL5Wr3YhA+4EMIPAEwRNuEDroFlLwAAYCqEH7gs9kgAACaD4ctebW1tevzxx/Xee++pt7dXdrtdGRkZWrBggSTp4MGDysnJ0aeffqpZs2YpOTlZ69atc87v7+/Xjh07tGfPHnV0dOjKK6/Uxo0bNXcua+Tujj0SAIDJYPidnzvvvFNffvmlCgoK9Morr8jHx0e33Xaburu71d7erjVr1mjevHkqKytTamqqtm3bprKyMuf8/Px8lZSUaPPmzSotLZXFYlFKSoocDoeBV4WJMrhHYixf7R09RpcNAHBhht75aW9v109+8hPdeeeduvDCCyVJd911l/71X/9Vf/7zn1VZWSmr1ars7Gx5eXkpPDxcjY2NKigoUGJiohwOh4qLi5Wenq64uDhJUl5enmJjY1VeXq6VK1caeXkAAMAFGXrnJzAwUE8++aQz+Bw7dkxFRUWy2WxasGCBampqZLfb5eX1fUaLiYlRQ0OD2traVF9fr87OTsXExDj7AwICFBkZqerq6im/HgAA4PoM3/Mz6OGHH9bLL78sq9Wq//iP/5Cfn59aWlq0cOHCIeNCQkIkSU1NTWppaZEkhYWFDRvT3Nw8NYUDAAC3Yvien0GrV69WWVmZbrjhBt19993605/+pJ6eHlmt1iHjvL29JUm9vb3q7u6WpBHH9Pb2Tk3hAADArbjMnZ/Bp7seffRRHThwQLt375aPj8+wjcuDocbPz08+Pj6SJIfD4fx+cIyvLweDAQCA4Qy989PW1qb9+/fr5MmTzjYPDw+Fh4ertbVVNptNra2tQ+YMfg4NDXUud400xmazTXL1AADAHRkaflpbW/XrX/9aH3zwgbPtxIkTqqurU3h4uOx2u2pra4eEo8rKSs2fP19BQUGKiIiQv7+/qqqqnP0dHR2qq6tTdHT0lF4LAABwD4aGn4iICF177bXatGmTampqdOjQIT344IPq6OjQbbfdpsTERH333XfKzMzU559/rr1792rnzp1av369pFN7fZKSkpSbm6uKigrV19frvvvuk81mU3x8vJGXBgAAXJShe34sFoueeuop/fa3v9W9996rv/71r4qOjtaLL76oc889V5JUWFionJwcJSQkKDg4WBkZGUpISHD+jLS0NPX19SkrK0s9PT2y2+0qKioatgkaAABAcoENzzNnzlR2drays7NH7F+0aJFKS0t/dL6np6fS09OVnp4+SRUCAIDpxGUedQcAAJgKhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8ApuHvN0P9/QNGlwHAYF5GFwAAU8XX6iUPD4tKyw/paHvXmOYuPD9Qy2PmTlJlAKYS4QeA6Rxt71LTsc4xzQme5TtJ1QCYaix7AQAAUyH8AAAAUyH8AAAAUyH8AAAAUzE8/Hz77bd65JFHdN111+mKK67QzTffrJqaGmf/wYMHlZSUpMWLF2vJkiUqKioaMr+/v1/bt29XbGysoqKitHbtWjU2Nk71ZQAAADdhePi5//779fHHH+vJJ5/UK6+8oksuuUTr1q3T4cOH1d7erjVr1mjevHkqKytTamqqtm3bprKyMuf8/Px8lZSUaPPmzSotLZXFYlFKSoocDoeBVwUAAFyVoY+6NzY26v3339fvfvc7XXHFFZKkzMxMvfvuu9q/f798fHxktVqVnZ0tLy8vhYeHq7GxUQUFBUpMTJTD4VBxcbHS09MVFxcnScrLy1NsbKzKy8u1cuVKIy8PAAC4IEPv/AQGBuq5557TpZde6myzWCwaGBjQ8ePHVVNTI7vdLi+v7zNaTEyMGhoa1NbWpvr6enV2diomJsbZHxAQoMjISFVXV0/ptQAAAPdgaPgJCAhQXFycrFars+3NN9/UX/7yF1177bVqaWmRzWYbMickJESS1NTUpJaWFklSWFjYsDHNzc2TXD0AAKec6atTeO3K1HKpE55ra2u1YcMGLVu2TEuXLtWWLVuGBCNJ8vb2liT19vaqu7tbkkYcc/z48akpGgBgemfy6pTgQD+til84SZVhJC4Tft5++2098MADioqK0pNPPilJ8vHxGbZxube3V5Lk5+cnHx8fSZLD4XB+PzjG15ej6AEAU2s8r07B1DP8aS9J2r17t1JTU3XdddepoKDAGWRsNptaW1uHjB38HBoa6lzuGmnMD5fLAAAAJBcIPy+99JIeffRR3XrrrXrqqaeGLGHZ7XbV1tbq5MmTzrbKykrNnz9fQUFBioiIkL+/v6qqqpz9HR0dqqurU3R09JReBwAA48F+oaln6LJXQ0ODHnvsMcXHx2v9+vVqa2tz9vn4+CgxMVGFhYXKzMzU7bffrk8++UQ7d+7Upk2bJJ3a65OUlKTc3FzNnj1bc+bM0datW2Wz2RQfH2/UZQEAMGrsF5p6hoaf3//+9zpx4oTKy8tVXl4+pC8hIUGPP/64CgsLlZOTo4SEBAUHBysjI0MJCQnOcWlpaerr61NWVpZ6enpkt9tVVFQ0bBM0AACujP1CU8fQ8HPHHXfojjvu+LtjFi1apNLS0h/t9/T0VHp6utLT0ye6PAAAMA0ZvucHAABgKhF+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+AACAqRB+MKn6+weMLgEAgCG8jC4A05uHh0Wl5Yd0tL1rTPMWnh+o5TFzJ6kqAICZEX4w6Y62d6npWOeY5gTP8p2kagAAZseyFwAAMBXCDwAAMBXCDwAAMBXCDwAAMBXCDwAAMBXCD06Ls3oAANOJSz3qnp+fr8rKSu3atcvZdvDgQeXk5OjTTz/VrFmzlJycrHXr1jn7+/v7tWPHDu3Zs0cdHR268sortXHjRs2dyxkxE4WzegAA04nL3Pl54YUXtH379iFt7e3tWrNmjebNm6eysjKlpqZq27ZtKisrc47Jz89XSUmJNm/erNLSUlksFqWkpMjhcEz1JUxrg2f1jOWrvaPH6LIBABjG8Ds/R44cUWZmpmprazV//vwhfS+//LKsVquys7Pl5eWl8PBwNTY2qqCgQImJiXI4HCouLlZ6erri4uIkSXl5eYqNjVV5eblWrlxpxCUBAAAXZvidnz/96U86++yz9dprrykqKmpIX01Njex2u7y8vs9oMTExamhoUFtbm+rr69XZ2amYmBhnf0BAgCIjI1VdXT1l1wAAANyH4Xd+li5dqqVLl47Y19LSooULFw5pCwkJkSQ1NTWppaVFkhQWFjZsTHNz8yRUCwAA3J3hd37+np6eHlmt1iFt3t7ekqTe3l51d3dL0ohjent7p6ZIAADgVlw6/Pj4+AzbuDwYavz8/OTj4yNJI47x9eXFmAAAYDiXDj82m02tra1D2gY/h4aGOpe7Rhpjs9mmpkgAAOBWXDr82O121dbW6uTJk862yspKzZ8/X0FBQYqIiJC/v7+qqqqc/R0dHaqrq1N0dLQRJQMAABfn0uEnMTFR3333nTIzM/X5559r79692rlzp9avXy/p1F6fpKQk5ebmqqKiQvX19brvvvtks9kUHx9vcPUAAMAVGf60198TFBSkwsJC5eTkKCEhQcHBwcrIyFBCQoJzTFpamvr6+pSVlaWenh7Z7XYVFRUN2wQNAAAguVj4efzxx4e1LVq0SKWlpT86x9PTU+np6UpPT5/M0gAAwDTh0steAAAAE43wAwAATIXwAwAATIXwAwAATIXwAwAATIXwAwAATIXwAwAATIXwAwAATIXwAwCAm/L3m6H+/oFxzz+Tue7MpU54BgAAo+dr9ZKHh0Wl5Yd0tL1rTHODA/20Kn7hJFXm2gg/AAC4uaPtXWo61ml0GW6DZS8AAGAqhB8AAGAqhB8AAGAqhB8AAGAqhB8AAEzIzI/J87QXAAAmZObH5Ak/AACYmBkfk2fZCwAAmArhBwAAmArhBwAAmArhBwAAmArhBwAAmArhBwAAmArhxyTc+TAqAAAmEuf8mMR4D7JaeH6glsfMnaSqAACYeoQfExnPQVbBs3wnqRoAAIzBshcAADAVwg8AADAVwg8AADAVwg8AADAVwg8AADAVwo8b4aweAADOHI+6uxHO6gEA4MwRftwMZ/UAAIzm7zdD/f0D8vCwjGv+mcydCIQfAAAwJr5Wr3GvRgQH+mlV/MJJqmx0CD8AAGBcxrMa4QqmxYbn/v5+bd++XbGxsYqKitLatWvV2NhodFkAAMAFTYvwk5+fr5KSEm3evFmlpaWyWCxKSUmRw+EwujQAAOBi3D78OBwOFRcXKzU1VXFxcYqIiFBeXp6OHDmi8vJyo8sDAAAuxu3DT319vTo7OxUTE+NsCwgIUGRkpKqrqw2sDAAAuCLLwMCAW5+c99Zbbyk1NVUff/yxfHx8nO333HOPenp69Oyzz47p51122WU6efKkwsLCJrrUCdHZfUInx3jY4QwvD/l6ezGXucxlLnNd8M8221xPD4vO8p0xpjmj0dzcLE9PT/3P//zPace6/dNe3d3dkiSr1Tqk3dvbW8ePHx/zz/P29nbpvUJn8heGucxlLnOZ67p/ttnmTjQvL69hWeBHx05yLZNu8G6Pw+EYcuent7dXvr5jP9yvpqZmwmoDAACux+33/AwuT7W2tg5pb21tlc1mM6IkAADgwtw+/ERERMjf319VVVXOto6ODtXV1Sk6OtrAygAAgCty+2Uvq9WqpKQk5ebmavbs2ZozZ462bt0qm82m+Ph4o8sDAAAuxu3DjySlpaWpr69PWVlZ6unpkd1uV1FR0ag3PgEAAPNw+0fdAQAAxsLt9/wAAACMBeEHAACYCuEHAACYCuEHAACYCuEHAACYCuEHAACYCuEHAACYCuFnjIqKirRs2TItWrRIN910k/77v//b6JLcXk1NjS6++OIhryjB6DU3N+v+++/XNddcI7vdrnXr1unPf/6z0WW5jf7+fm3fvl2xsbGKiorS2rVr1djYaHRZbufbb7/VI488ouuuu05XXHGFbr75Zl4UPQEaGhp0+eWXa+/evUaXMq0QfsYgPz9fTz/9tO655x699tprWrRoke688059+eWXRpfmtv76178qIyND/f39RpfilhwOh375y1+qra1Nzz77rF566SXNnDlTq1ev1jfffGN0eW4hPz9fJSUl2rx5s0pLS2WxWJSSkiKHw2F0aW7l/vvv18cff6wnn3xSr7zyii655BKtW7dOhw8fNro0t3XixAk98MAD6urqMrqUaYfwM0pdXV0qKChQenq6brjhBs2bN08PP/ywzjvvPNXW1hpdntvKzs7WeeedZ3QZbqumpkaHDh3SE088oUsvvVQXXnihnnjiCXV1demdd94xujyX53A4VFxcrNTUVMXFxSkiIkJ5eXk6cuSIysvLjS7PbTQ2Nur999/Xxo0bFR0drQsuuECZmZkKDQ3V/v37jS7PbT399NM666yzjC5jWiL8jFJNTY26u7u1cuVKZ5unp6dee+013XjjjcYV5sZeffVVffTRR9qwYYPRpbitCy+8UM8995xCQ0OHtA8MDOj48eMGVeU+6uvr1dnZqZiYGGdbQECAIiMjVV1dbWBl7iUwMFDPPfecLr30UmebxWLh7+EZqK6uVmlpqX7zm98YXcq0RPgZpS+++EJnn322PvvsM9188826+uqrlZycrA8//NDo0tzSV199pZycHD3xxBP8n80ZCA4OVlxc3JC2//zP/1Rvb6+uueYag6pyHy0tLZKksLCwIe0hISFqbm42oiS3FBAQoLi4uCEvk37zzTf1l7/8Rddee62Blbmnjo4OZWRkKCsra9jfTUyMafFW94nw1VdfadmyZT/af88996inp0ePPPKIfv3rX+vcc89VaWmpVq9erX379ik8PHwKq3Vtp/tdvvvuu8rIyNCqVasUHR2tr776agqrcy+n+13+8Y9/VHBwsPPzW2+9pby8PCUnJysiImIqSnRr3d3dkjTkH21J8vb25o7FGaitrdWGDRu0bNkyLV261Ohy3E52drYWL16s66+/3uhSpi3Cz9+Ehobqv/7rv360v6KiQj09PdqwYYPz/7QvueQSffTRR9q9e7c2btw4VaW6vNP9Lvfs2aOuri6lpqZOYVXu6XS/y9mzZzu//93vfqdHH31UK1as0EMPPTQV5bk9Hx8fSaf2/gx+L0m9vb3y9fU1qiy39vbbb+uBBx5QVFSUnnzySaPLcTv79u1TTU2NXn/9daNLmdYIP38zY8aMv3v3pq6uTpJ00UUXOdssFovCw8O5c/EDp/td7t27V62trfrpT38q6dT+FElKSUnRVVddpcLCwimp0x2c7nc5KDc3VwUFBUpOTlZmZqYsFssUVOf+BpcUWltbdf755zvbW1tbuXM2Drt371ZOTo7i4+OVm5s77I4aTq+srExtbW1asmTJkPaNGzeqqKhIb7zxhjGFTTOEn1GKjo6WxWLRgQMH9E//9E+STv2j/fnnn+vqq682uDr3smvXLvX19Tk/HzlyRMnJydq8ebMzEGH0tm7dqsLCQmVkZGjdunVGl+NWIiIi5O/vr6qqKmf46ejoUF1dnZKSkgyuzr289NJLevTRR5WcnKwNGzbIw4MtpeORm5urnp6eIW3Lly9XWlqaVqxYYVBV0w/hZ5TCwsKUmJiozZs3y9fXV3PnztWuXbv01Vdf6ZZbbjG6PLcyZ86cIZ89PT0lnVri+eFTS/j7qqqqVFhYqOTkZN1www06evSos8/Pz4/N5KdhtVqVlJSk3NxczZ49W3PmzNHWrVtls9kUHx9vdHluo6GhQY899pji4+O1fv16tbW1Oft8fHw0c+ZMA6tzLz/238CgoKBh/+3E+BF+xiA7O1s7duxQVlaWjh8/rsjISBUXF+uCCy4wujSY1OAZKrt27dKuXbuG9P3qV79iX9UopKWlqa+vT1lZWerp6ZHdbldRURFLNmPw+9//XidOnFB5efmw85ESEhL0+OOPG1QZMDLLwOCGCwAAABNgURYAAJgK4QcAAJgK4QcAAJgK4QcAAJgK4QcAAJgK4QcAAJgK4QcAAJgK4QcAAJgK4QcAAJgK4QcAAJgK4QcAAJjK/weQtEWgnwBTiQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.histplot(results.null_distribution)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "b28f3a42-1a88-403e-8363-1921349f45a1", "metadata": {}, "source": [ "We can see that indeed, only in a few cases did the mean difference exceed 4.6.\n", "\n", "The distribution of differences of means from shuffled data looks pretty odd (with two peaks). This probably reflects the fact that whether the difference (Husbands-Wives) is positive or negative is strongly influenced by whether the hich-sock person (with 48 pairs of socks) happens to be labelled as a husband or a wife; one peak in the null distribution will be shuffles where that person was labelled as male, and the other rpeak will be shuffles where they were labelled as female.\n", "\n", "### The $p$ value\n", "\n", "We can calculate the exact proportion of cases in which the mean difference in socks owned for [Husband-Wife] exceeds the value we observed in our original sample, 4.6, as follows:" ] }, { "cell_type": "code", "execution_count": 12, "id": "6308bbce-3a50-42d4-a63a-dfabdb057788", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "0.0072" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check which values is the null distribution are greater than 4.6 and code as 1/0\n", "# then use the numpy function np.mean() to get the avergae of those ones and zeros (ie he proportion of 1s)\n", "np.mean(results.null_distribution >= 4.57)" ] }, { "cell_type": "markdown", "id": "b972dcbc-065d-4b83-a5b5-badcc22fb372", "metadata": {}, "source": [ "The proportion of shuffles in which the mean difference is greater than 4.6 is 0.004, 0.4%, or 40 of the 10,000 shuffles; we can check the last one a different way (just as a sanity check!)" ] }, { "cell_type": "code", "execution_count": 13, "id": "32e8dd58-cc0d-480b-a067-70a8a564b559", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "72" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sum(results.null_distribution >= 4.57)" ] }, { "cell_type": "markdown", "id": "be3cfbb8-b609-427e-b52a-cb34ac6c1bb0", "metadata": {}, "source": [ "#### Interpretation\n", "\n", "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.4% 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.4% or p=0.0004.\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": "code", "execution_count": 9, "id": "966fc15b-1cb0-4b49-ab71-861257d7cde8", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "10000\n" ] } ], "source": [ "print(results.null_distribution.size)" ] }, { "cell_type": "markdown", "id": "1fdec81e", "metadata": {}, "source": [ "## Syntax of `stats.permutation_test()`\n", "\n", "Let's have a look at the Python code to run the permutation test.\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* that calculates that test statistic. We created that function, `mDiff()`, using a function definition, before running the permutation test\n", "\n", "### Arguments\n", "\n", "When running the permutation test, we have to specify some arguments\n", "\n", "* **n_resamples** the number of shuffles to do.\n", "* **alternative** - whether to do a one-tailed test (and in whcih direction: 'greater' counts how many shuffles result in test statistic values greater than the one obtained from unshuffled data) or a two tailed test\n", "* **permutation_type** - pairings, samples, or independent. These are for test of correlation, paired data and unpaired data respectively.\n", "\n", "Here we are testing for a mean difference in paired data, so we choose `permutation_type='samples'`; the samples in this case are the men and women (husbands and wives) - we are permuting the data by moving people between samples (flipping husband/wife labels) but not moving people between couples.\n" ] }, { "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": "deb12a97-6048-4e09-a9f7-3fefc682419b", "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.18" } }, "nbformat": 4, "nbformat_minor": 5 }