{ "cells": [ { "cell_type": "markdown", "id": "c33197dd", "metadata": {}, "source": [ "# The Wilcoxon Sign-Rank Test\n", "\n", "\n", "The t-test is valid only when the data within each group (for independent samples t-test) or the pairwise differences (for paired samples t-test) are Normally distributed\n", "\n", "As we have seen in the lecture, many real life data distributions are normal, but many others are not.\n", "\n", "For non-Normal data we can use non-parametric tests, which do not assume that the data are drawn from a Normal distribution.\n", "\n", "\n", "The Wilcoxon Sign-Rank Test is a test for paired samples. It tests whether the median difference between the members of each pair is greater than zero. As such it is often considered to be a non-parametric equivalent for the paired samples t-test.\n", "\n", "The Wilcoxon Sign-rank test is not the same as the Wilcoxon Rank Sum test (Mann Whitney U test) which is for independent samples\n", "\n", "We will us a Python function called wilcoxon from the scipy.stats package to run the test" ] }, { "cell_type": "markdown", "id": "b563240d", "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": 14, "id": "93d0feed", "metadata": {}, "outputs": [], "source": [ "# Set-up Python libraries - you need to run this but you don't need to change it\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.stats as stats\n", "import pandas \n", "import seaborn as sns" ] }, { "cell_type": "markdown", "id": "8b8bd7b3", "metadata": {}, "source": [ "## Example: the Sign-Rank Test\n", "\n", "It has been argued that birth order in families affects how independent individuals are as adults - either that first-born children tend to be more independent than later born children or vice versa.\n", "\n", "In a (fictional!) study, a researcher identified 20 sibling pairs, each comprising a first- and second- born child from a two-child family. The participants were young adults; each participant was interviewed at the age of 21. \n", "\n", "The researcher scored independence for each participant, using a 25 point scale where a higher score means the person is more independent, based on a structured interview.\n", "\n", "Carry out a statistical test for a difference in independence scores between the first- and second-born children.\n", "\n", "Note that this is a paired samples design - each member of one group (the first-borns) has a paired member of the other group (second-borns).\n", "\n", "\n", "### Inspect the data\n", "\n", "The data are provided in a text (.csv) file.\n", "\n", "Let's load the data as a Pandas dataframe, and plot them to get a sense for their distribution (is it normal?) and any outliers\n" ] }, { "cell_type": "code", "execution_count": 15, "id": "da622d51", "metadata": {}, "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", " \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", "
FirstBornSecondBorn
01210
11812
21315
31713
489
51512
61613
758
8810
9128
10138
1159
12148
132010
141914
151711
1627
1757
181513
191812
\n", "
" ], "text/plain": [ " FirstBorn SecondBorn\n", "0 12 10\n", "1 18 12\n", "2 13 15\n", "3 17 13\n", "4 8 9\n", "5 15 12\n", "6 16 13\n", "7 5 8\n", "8 8 10\n", "9 12 8\n", "10 13 8\n", "11 5 9\n", "12 14 8\n", "13 20 10\n", "14 19 14\n", "15 17 11\n", "16 2 7\n", "17 5 7\n", "18 15 13\n", "19 18 12" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# load the data and have a look\n", "pandas.read_csv('https://raw.githubusercontent.com/jillxoreilly/StatsCourseBook/main/data/BirthOrderIndependence.csv')" ] }, { "cell_type": "markdown", "id": "5481f8e9", "metadata": {}, "source": [ "### Scatterplot\n", "\n", "In the case of paired data, the most effective way to get a sense of the data is a scatterplot:" ] }, { "cell_type": "code", "execution_count": 16, "id": "ce6a5214", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABkHElEQVR4nO3dd3xT1f/H8VcKpQUsQRBoyyyIQgGZsmQpMhW3IChDnIgIoj8QF/AVBcSJIjgYIl/BASgIIqAMkQ0tG2RUi9paWS1DWtqe3x/320LooClJkzTv5+ORB703555+bm7SfDj3DJsxxiAiIiLiRwI8HYCIiIhIQVMCJCIiIn5HCZCIiIj4HSVAIiIi4neUAImIiIjfUQIkIiIifkcJkIiIiPidop4OwBulp6fz119/ERISgs1m83Q4IiIikgfGGE6ePEl4eDgBAbm38SgBysZff/1F5cqVPR2GiIiI5MPhw4epVKlSrmWUAGUjJCQEsF7AUqVKeTgaERERyYukpCQqV66c+T2eGyVA2ci47VWqVCklQCIiIj4mL91X1AlaRERE/I4SIBEREfE7SoBERETE7ygBEhEREb+jBEhERET8jhIgERER8TtKgERERMTvKAESERERv6MESERERPyOEiARERHxOx5NgMaOHcv1119PSEgI5cuX54477mDfvn0OZYwxjBo1ivDwcIoXL067du3YtWvXJeueO3cukZGRBAUFERkZyfz58911GiIiIuJjPJoArVq1ioEDB7J+/XqWLVtGamoqHTt25PTp05llXn/9dd566y3ef/99Nm3aRGhoKB06dODkyZM51rtu3Tp69OhB79692bZtG71796Z79+5s2LChIE5LREREvJzNGGM8HUSGf/75h/Lly7Nq1SratGmDMYbw8HCGDBnC8OHDAUhOTqZChQqMHz+exx57LNt6evToQVJSEt9//33mvs6dO3PllVcye/bsS8aRlJSE3W4nMTFRi6GKiIj4CGe+v72qD1BiYiIAZcqUASAmJob4+Hg6duyYWSYoKIi2bduydu3aHOtZt26dwzEAnTp1yvGY5ORkkpKSHB4iIiJSeHlNAmSMYejQobRq1Yq6desCEB8fD0CFChUcylaoUCHzuezEx8c7dczYsWOx2+2Zj8qVK1/OqYiIiMjFjhyBhARPR5HJaxKgJ598ku3bt2d7i8pmszlsG2Oy7LucY0aMGEFiYmLm4/Dhw05GLyIiIjlavRrq14devSAtzdPRAF6SAA0aNIgFCxawYsUKKlWqlLk/NDQUIEvLTUJCQpYWnguFhoY6dUxQUBClSpVyeIiIiMhlSk+HV1+FG2+Ev/6CP//0mlYgjyZAxhiefPJJ5s2bx08//URERITD8xEREYSGhrJs2bLMfSkpKaxatYqWLVvmWG+LFi0cjgFYunRprseIiIiIC/39N3TuDC++aCVCffrApk0QFubpyAAo6slfPnDgQD7//HO+/fZbQkJCMltt7HY7xYsXx2azMWTIEF577TVq1qxJzZo1ee211yhRogS9evXKrKdPnz5UrFiRsWPHAjB48GDatGnD+PHjuf322/n2229Zvnw5a9as8ch5ioiI+JWffoL774f4eChRAj74APr29XRUDjyaAE2ePBmAdu3aOeyfPn06/fr1A2DYsGH8+++/PPHEExw/fpxmzZqxdOlSQkJCMsvHxsYSEHC+Matly5bMmTOHF198kZdeeokaNWrwxRdf0KxZM7efk4iIiF9LTYUnn7SSnzp14MsvITLS01Fl4VXzAHkLzQMkIiJyGbZtgylT4M03rRagAuKz8wCJiIiID1q6FD7++Px2/foweXKBJj/OUgIkIiIi+ZOaCi+8YHV2HjgQtm71dER55tE+QCIiIuKj/vgDevaEjAFGDz3klX19cqIESERERJyzeLE1rP3oUQgJgU8+ge7dPR2VU3QLTERERPLuhRfgllus5KdRI4iK8rnkB5QAiYiIiDP+t2A5gwbB2rVQo4Zn48kn3QITERGR3J0+DSVLWj8PHQrNmkGrVp6N6TKpBUhERESyl5ICQ4ZAkyZw6pS1z2bz+eQHlACJiIhIdg4dghtugHffhb17YeFCT0fkUkqARERExNHcudCwIWzeDFdeCQsWWEPeCxElQCIiImI5e9Zax+ueeyApCVq2hOho6NbN05G5nBIgERERsfzf/8GkSdbPw4fDypVQpYpHQ3IXJUAiIiJieeEFqFsXvv8exo2DwEBPR+Q2SoBERET81b//wuefn98ODbVWcu/c2XMxFRDNAyQiIuKP9u61ZnDesQOKFj0/m3OAf7SN+MdZioiIyHkzZ0LjxlbyU778+dmd/YgSIBEREX9x+jT07w99+8KZM3DTTdYor5tv9nRkBU4JkIiIiD/YtQuaNoXp063bXKNHw9KlEBbm6cg8Qn2ARERE/MHBg7B7t5XwfP45tGvn6Yg8SgmQiIhIYWWMtXYXwG23wSefWJMali/v2bi8gG6BiYiIFEbbtlmLlh4+fH7fQw8p+fkfJUAiIiKFiTHw4YfQrBmsXQvPPOPpiLySboGJiIgUFklJ8Oij8MUX1vYtt8AHH3g2Ji+lFiAREZHCYOtWa26fL76wJjacMMFaxf2qqzwdmVdSC5CIiIivW7HCWr4iJcVavPSLL6B5c09H5dWUAImIiPi65s3h2muhenWYNs0vZ3Z2lhIgERERX7RrF9SqBUWKQPHiVitQmTLnh71LrtQHSERExJcYA2+/DQ0bwtix5/eXLavkxwlqARIREfEVx45Bv36wcKG1vXOn42SHkmdqARIREfEFa9dCgwZW8lOsGEyaBLNnK/nJJyVAIiIi3iw9HV5/Hdq0sWZ1vvpqWL8ennhCyc9l8GgCtHr1arp160Z4eDg2m41vvvnG4XmbzZbtY8KECTnWOWPGjGyPOXv2rJvPRkRExA0OHoSXX4a0NOjZ05rvp2FDT0fl8zzaB+j06dPUr1+fBx98kLvvvjvL83FxcQ7b33//PQ899FC2ZS9UqlQp9u3b57AvODj48gMWEREpaDVrwvvvW319Hn5YrT4u4tEEqEuXLnTp0iXH50NDQx22v/32W2688UaqV6+ea702my3LsSIiIj4hPR3GjYObb4amTa19Dz/s2ZgKIZ/pA/T333+zaNEiHnrooUuWPXXqFFWrVqVSpUrceuutREVF5Vo+OTmZpKQkh4eIiEiB+/tva0bnF16AHj3g9GlPR1Ro+UwC9OmnnxISEsJdd92Va7latWoxY8YMFixYwOzZswkODuaGG25g//79OR4zduxY7HZ75qNy5cquDl9ERCR3P/1kjfJatsya2HDkSChZ0tNRFVo2Y4zxdBBg3baaP38+d9xxR7bP16pViw4dOvDee+85VW96ejqNGjWiTZs2TJw4MdsyycnJJCcnZ24nJSVRuXJlEhMTKVWqlFO/T0RExClpafDKK/Cf/1j9fOrUgS+/hMhIT0fmc5KSkrDb7Xn6/vaJiRB//vln9u3bxxdffOH0sQEBAVx//fW5tgAFBQURFBR0OSGKiIg4LykJbr8dVq60tvv3h/fegxIlPBqWP/CJW2BTp06lcePG1K9f3+ljjTFER0cTFhbmhshEREQuwxVXWLe5SpaEzz6DqVOV/BQQj7YAnTp1igMHDmRux8TEEB0dTZkyZahSpQpgNWd99dVXvPnmm9nW0adPHypWrMjY/62HMnr0aJo3b07NmjVJSkpi4sSJREdHM2nSJPefkIiIyKWkpsK5c1Y/n4AA+PRTOHLEWs1dCoxHE6DNmzdz4403Zm4PHToUgL59+zJjxgwA5syZgzGGnj17ZltHbGwsAQHnG7JOnDjBo48+Snx8PHa7nYYNG7J69WqaZgwlFBER8ZQ//oBevSAiwkp8wFrEtGxZz8blh7ymE7Q3caYTlYiISJ4sXgx9+sDRoxASAtu3Q7Vqno6qUHHm+9sn+gCJiIj4rHPnYNgwuOUWK/lp1MhazkLJj0f5xCgwERERnxQbC/fdB+vWWduDBsGECaCRxx6nBEhERMQd0tOtWZ337AG7HaZNg0tM5isFR7fARERE3CEgAN59F5o3h6goJT9eRgmQiIiIqxw6ZC1lkaFDB/jlF2vUl3gVJUAiIiKuMHcuNGwI99wDBw+e3x+gr1pvpKsiIiJyOc6ehSeftBKfpCRrLa/AQE9HJZegBEhERCS/9u+Hli0hY7WBYcNg1Sr432oG4r00CkxERCQ/5syBRx+FkyetmZxnzoSuXT0dleSREiAREZH82LDBSn5at4bPP4dKlTwdkThBCZCIiEheGQM2m/Xz+PFw9dXw2GNQVF+nvkZ9gERERPJi1ixrOYvUVGu7WDEYOFDJj49SAiQiIpKb06ehf3/o3Ru+/x6mT/d0ROICSltFRERysmsXdO8Ou3dbt75GjrSSIfF5SoBEREQuZgzMmGHd4vr3XwgNtTo633ijpyMTF9EtMBERkYuNHm219Pz7r7WcxbZtSn4KGSVAIiIiF+vRA0qVgldfhSVLoHx5T0ckLqZbYCIiIsZYrTwNGljbtWtDTAyUKePRsMR91AIkIiL+LSkJevWCxo3h55/P71fyU6gpARIREf8VFWUlPnPmWKO89uzxdERSQJQAiYiI/zHGWsC0eXM4cMBavPTnn621vcQvqA+QiIj4lxMn4OGHYe5ca/u226zJDXXLy6+oBUhERPzLN99YyU9gILz9trWt5MfvqAVIRET8S9++sH079OwJ11/v6WjEQ9QCJCIihduxY/DEE5CYaG3bbPDWW0p+/Fy+WoB+/PFHfvzxRxISEkhPT3d4btq0aS4JTERE5LKtWwf33QexsVYC9N//ejoi8RJOtwCNHj2ajh078uOPP3LkyBGOHz/u8BAREfG49HSYMAHatLGSnxo14JlnPB2VeBGnW4CmTJnCjBkz6N27tzviERERuTxHjlj9fBYvtrZ79ICPPrKWthD5H6cToJSUFFq2bOmOWERERC5PdDTceiv8+ScEBcHEifDII1a/H5ELOH0L7OGHH+bzzz93RywiIiKXp1Il699rr4WNG62JDZX8SDacbgE6e/YsH330EcuXL+e6664jMDDQ4fm33nrLZcGJiIj7pKUbNsYcI+HkWcqHBNM0ogxFAnwwWUhKOn9766qr4IcfoGpVuOIKz8YlXs3pFqDt27fToEEDAgIC2LlzJ1FRUZmP6Ohop+pavXo13bp1Izw8HJvNxjfffOPwfL9+/bDZbA6P5s2bX7LeuXPnEhkZSVBQEJGRkcyfP9+puERECrslO+NoNf4nen68nsFzoun58Xpajf+JJTvjPB2ac1assFp7Pv30/L46dZT8yCU51QKUlpbGqFGjqFevHmVcMGvm6dOnqV+/Pg8++CB33313tmU6d+7M9OnTM7eLFSuWa53r1q2jR48evPLKK9x5553Mnz+f7t27s2bNGpo1a3bZMYuI+LolO+MYMGsr5qL98YlnGTBrK5MfaETnumEeiS3P0tJgzBj4z3+sEV+TJkHv3hCg6e0kb2zGmIs/A7kKDg5mz549REREuDYQm4358+dzxx13ZO7r168fJ06cyNIylJsePXqQlJTE999/n7mvc+fOXHnllcyePTtPdSQlJWG320lMTKSURg2ISCGSlm5oNf4n4hLPZvu8DQi1B7Nm+E3eezssLg4eeAB++snafvBBeO89KFnSs3GJxznz/e10qlyvXj0OHTqU7+CctXLlSsqXL88111zDI488QkJCQq7l161bR8eOHR32derUibVr1+Z4THJyMklJSQ4PEZHCaGPMsRyTHwADxCWeZWPMsYILyhnLlkGDBlbyU7IkzJwJ06Yp+RGnOZ0Avfrqqzz77LN89913xMXFuTVx6NKlC//973/56aefePPNN9m0aRM33XQTycnJOR4THx9PhQoVHPZVqFCB+Pj4HI8ZO3Ysdrs981G5cmWXnYOIiDdJOJlz8pOfcgXq0CHo0gUSEqBePdi82brtJZIPTo8C69y5MwC33XYbtguGFhpjsNlspKWluSy4Hj16ZP5ct25dmjRpQtWqVVm0aBF33XVXjsfZLhrymBFbTkaMGMHQoUMzt5OSkpQEiUihVD4k2KXlClT16jB8OBw9aq3iXry4pyMSH+Z0ArRixQp3xJEnYWFhVK1alf379+dYJjQ0NEtrT0JCQpZWoQsFBQURFBTksjhFRLxV04gyhNmDiU88m6UTNJzvA9Q04vIHurjE999bo7yqV7e2x4zRvD7iEk4nQG3btnVHHHly9OhRDh8+TFhYzqMTWrRowbJly3j66acz9y1dulSzV4uIAEUCbIzsFsmAWVuxgUMSlJFWjOwW6fkO0OfOwQsvWOt5XX89rFkDxYop+RGXyddq8CdOnGDq1Kns2bMHm81GZGQk/fv3x263O1XPqVOnOHDgQOZ2TEwM0dHRlClThjJlyjBq1CjuvvtuwsLC+O2333j++ee56qqruPPOOzOP6dOnDxUrVmTs2LEADB48mDZt2jB+/Hhuv/12vv32W5YvX86aNWvyc6oiIoVO57phTH6gEaMX7nboEB1qD2Zkt0jPD4GPjbVWcF+3ztpu2hScG7AscmnGSZs2bTJlypQxFStWNHfeeae54447TKVKlUzZsmXNli1bnKprxYoVBus/IA6Pvn37mjNnzpiOHTuacuXKmcDAQFOlShXTt29fExsb61BH27ZtTd++fR32ffXVV+baa681gYGBplatWmbu3LlOxZWYmGgAk5iY6NRxIiK+JDUt3aw9cMR8E/WHWXvgiElNS/d0SMZ8+60xV15pDBhjtxvz9deejkh8iDPf307PA9S6dWuuvvpqPv74Y4oWtRqQUlNTefjhhzl06BCrV692bYbmAZoHSESkgKWkwHPPWZ2bwbrtNWfO+b4/Inng1nmANm/ezPDhwzOTH4CiRYsybNgwNm/e7Hy0IiIixkDGf6CHDLH6/Cj5ETdyug9QqVKliI2NpVatWg77Dx8+TEhIiMsCExERP2CM1bE5KAi+/BJ27IDbb/d0VOIHnG4B6tGjBw899BBffPEFhw8f5o8//mDOnDk8/PDD9OzZ0x0xiohIYZOcDIMGwcsvn99XvbqSHykwTrcAvfHGG9hsNvr06UNqaioAgYGBDBgwgHHjxrk8QBERKWQOHIAePWDrVmvx0r594eqrPR2V+BmnO0FnOHPmDAcPHsQYw9VXX02JEiVcHZvHqBO0iIibfPklPPwwnDwJZcvCp5/CLbd4OiopJJz5/s7XPEAAJUqUoHTp0thstkKV/IiIiBv8+y88/TR8+KG13aoVzJ4NlSp5Ni7JVlq6YWPMMRJOnqV8iDUzuMcnx3QxpxOg1NRURo8ezcSJEzl16hQAV1xxBYMGDWLkyJEEBga6PEgREfFhxsDNN8PatVaH5xEjYPRoKJrv/4OLGy3ZGZdlkswwb5kk04Wcfvc9+eSTzJ8/n9dff50WLVoAsG7dOkaNGsWRI0eYMmWKy4MUEREfZrPBI4/A/v0waxZ07OjpiCQHS3bGMWDW1izrxMUnnmXArK1MfqBRoUmCnO4DZLfbmTNnDl26dHHY//3333PfffeRmJjo0gA9QX2AREQu05kz8PvvULv2+X3Hj8OVV3ouJslVWrqh1fifHFp+LpSxUO6a4Td57e0wt06EGBwcTLVq1bLsr1atGsWKFXO2OhERKWx277bW7+rYEY4ePb9fyY9X2xhzLMfkB6y1quISz7Ix5ljBBeVGTidAAwcO5JVXXiE5OTlzX3JyMq+++ipPPvmkS4MTEREfM2MGNGkCu3ZBair89punI5I8SjiZc/KTn3LeLk99gO666y6H7eXLl1OpUiXq168PwLZt20hJSaF9+/auj1BERLzfqVMwcCDMnGlt33yz1d+nQgXPxiV5Vj4k2KXlvF2eEiC73e6wfffddztsV65c2XURiYiIb9mxA7p3h717rYkN//Mfa6RXgNM3GcSDmkaUIcweTHzi2SydoOF8H6CmEWUKOjS3yFMCNH36dHfHISIivmr8eCv5CQ+35vZp08bTEUk+FAmwMbJbJANmbcUGDklQRpfnkd0ivbYDtLOUnouIyOWZNMma3Tk6WsmPj+tcN4zJDzQi1O54myvUHlyohsDDZSyFUZhpGLyISC6iouDzz+H11605fqTQ8dWZoAtkKQwREfEzxsDkydaSFikpEBkJDz7o6ajEDYoE2GhRo6ynw3ArJUAiInJpiYnWba6vv7a2u3WD22/3bEwil0F9gEREJHebNkHDhlbyExgIb70F334LZQrHaCDxT3lqAZo4cWKeK3zqqafyHYyIiHiZadPg8cfh3DmoVg2++MKa5VnEx+UpAXr77bcdtv/55x/OnDlD6dKlAThx4gQlSpSgfPnySoBERAqTq6+GtDS46y6YOhX+93dfxNfl6RZYTExM5uPVV1+lQYMG7Nmzh2PHjnHs2DH27NlDo0aNeOWVV9wdr4iIuNuJE+d/btMGNmywbn8p+ZFCxOlh8DVq1ODrr7+mYcOGDvu3bNnCPffcQ0xMjEsD9AQNgxcRv5SebvXvefVVWLcOatXydERZuHN4trvq9tUh5b7IrcPg4+LiOHfuXJb9aWlp/P33385WJyIi3uDIEejXDxYtsrY/+8xKhLzIkp1xjF6422HF8jB7MCO7RV72BH3uqtudMcvlcXoUWPv27XnkkUfYvHkzGY1Hmzdv5rHHHuPmm292eYAiIuJma9ZYo7wWLYKgIJgyBcaM8XRUDpbsjGPArK0OiQRAfOJZBszaypKdcV5XtztjlsvndAI0bdo0KlasSNOmTQkODiYoKIhmzZoRFhbGJ5984o4YRUTEHdLTYexYaNcO/vgDrrnG6u/z2GNeNcNzWrph9MLd2S7QmbFv9MLdpKU7v7CBu+p2Z8ziGk7fAitXrhyLFy/m119/Ze/evRhjqF27Ntdcc4074hMREXeZMQOef976+YEHrFmer7jCoyFlZ2PMsSytKBcyQFziWTbGHHN69mJ31e3OmMU18j0T9DXXXKOkR0TEl/XpA3PmwH33WUtaeFGrz4USTuacSOSnXEHU7c6YxTWcToDS0tKYMWMGP/74IwkJCaSnpzs8/9NPP7ksOBERcaG0NGsun379oFgxKFoUfvjBaxOfDOVDgi9dyIlyBVG3O2MW13A6ARo8eDAzZszglltuoW7duti8/IMjIiJAfDzcfz/89BPs3WsNdwevT34AmkaUIcweTHzi2Wz71NiAULs1vNxb6nZnzOIaTidAc+bM4csvv6Rr167uiEdERFxt+XKrj8/ff0OJEtaILx9SJMDGyG6RDJi1FRs4JBQZ6dvIbpH5mlvHXXW7M2ZxDadHgRUrVoyrr77aJb989erVdOvWjfDwcGw2G998803mc+fOnWP48OHUq1ePkiVLEh4eTp8+ffjrr79yrXPGjBnYbLYsj7NndZ9VRPxMaiq89BJ07GglP/XqwZYt0Lu3pyNzWue6YUx+oBGhdsdbRqH2YCY/0Oiy5tRxV93ujFkun9MtQM888wzvvvsu77///mXf/jp9+jT169fnwQcf5O6773Z47syZM2zdupWXXnqJ+vXrc/z4cYYMGcJtt93G5s2bc623VKlS7Nu3z2FfcLDus4qIH/nzT+jVC1avtrYfeQTefReKF/dsXJehc90wOkSGumVWZXfV7c6Y5fI4nQCtWbOGFStW8P3331OnTh0CAwMdnp83b16e6+rSpQtdunTJ9jm73c6yZcsc9r333ns0bdqU2NhYqlSpkmO9NpuN0NDQPMchIlLo/PsvREVZw9o/+gh69vR0RC5RJMDmtmHj7qrbnTFL/jmdAJUuXZo777zTHbFcUmJiIjabLXMV+pycOnWKqlWrkpaWRoMGDXjllVeyrF12oeTkZJKTkzO3k5KSXBWyiEjBMeZ8p+arr4Yvv4QaNaBmTc/GJeKFnE6Apk+f7o44Luns2bM899xz9OrVK9cFzmrVqsWMGTOoV68eSUlJvPvuu9xwww1s27aNmjn8ERg7diyjR492V+giIu53+LA1yuvllyFjWaLOnT0bk4gXc3o1+Az//PMP+/btw2azcc0111CuXLnLC8RmY/78+dxxxx1Znjt37hz33nsvsbGxrFy50qkV2tPT02nUqBFt2rRh4sSJ2ZbJrgWocuXKWg1eRHzDwoXW3D7HjlnLWezeDUWKeDoqkQLnzGrwTo8CO336NP379ycsLIw2bdrQunVrwsPDeeihhzhz5ky+g87JuXPn6N69OzExMSxbtszphCQgIIDrr7+e/fv351gmKCiIUqVKOTxERLxeSgo88wzcdpuV/DRpAt9/r+RHJA+cToCGDh3KqlWrWLhwISdOnODEiRN8++23rFq1imeeecalwWUkP/v372f58uWULet8JzJjDNHR0YSFabihiBQiv/0GrVufn9Bw8GBrVffq1T0alrulpRvWHTzKt9F/su7gUb9fTDQlNZ2pPx/i5W93MvXnQ6Skpl/6IAHy0Qdo7ty5fP3117Rr1y5zX9euXSlevDjdu3dn8uTJea7r1KlTHDhwIHM7JiaG6OhoypQpQ3h4OPfccw9bt27lu+++Iy0tjfj4eADKlClDsWLFAOjTpw8VK1Zk7NixAIwePZrmzZtTs2ZNkpKSmDhxItHR0UyaNMnZUxUR8U6HD1uTGZ44AaVLw/TpkE33gcJmyc44Ri/c7bDIaJg9mJHdIv1yTp2xi3fz8c8xXJgDvrp4D4+0jmBE10jPBeYjnE6Azpw5Q4UKFbLsL1++vNO3wDZv3syNN96YuT106FAA+vbty6hRo1iwYAEADRo0cDhuxYoVmQlYbGwsAQHnG7JOnDjBo48+Snx8PHa7nYYNG7J69WqaNm3qVGwiIl6rUiXo1g3277cWM61a1dMRud2SnXEMmLU1y7IS8YlnGTBrq99NLDh28W4+XB2TZX+6IXO/kqDcOd0Jun379pQtW5aZM2dmTi7477//0rdvX44dO8by5cvdEmhBcqYTlYhIgTh40GrtyegKcOYMBAZaj0IuLd3QavxPDi0/F8pYV2vN8Jv8YoLBlNR0ar30Pbnd/Quwwd5XulCsqNM9XXyaWztBv/vuu6xdu5ZKlSrRvn17br75ZipXrszatWt599138x20iIjk4MsvrVteDz5ozfUD1ppefpD8AGyMOZZj8gPWOltxiWfZGHOs4ILyoM/W/ZZr8gNWS9Bn634rkHh8ldO3wOrWrcv+/fuZNWsWe/fuxRjDfffdx/33309xH55iXUTE65w9C08/DVOmWNvHjkFSEtjtno2rgCWczNtajnkt5+t+P5a37iZ5LeevnE6AAIoXL84jjzzi6lhERCTDr79C9+6wbZu1PWIE/Oc/UDRff7Z9WvmQvK3lmNdyvq5qmRIuLeevnL4FNnbsWKZNm5Zl/7Rp0xg/frxLghIR8Wv//S80amQlP+XKwZIl8Nprfpn8ADSNKEOYPZicevfYsEaDNY0oU5BheUzvFtW4VFenAJtVTnLmdAL04YcfUqtWrSz769Spw5SMZloREcmfM2fgxRfh9Glo1w6io6FTJ09H5VFFAmyM7GaNaLr4ez9je2S3SL/oAA1QrGgAj7SOyLXMI60j/K4DtLOcfnXi4+OznVSwXLlyxMXFuSQoERG/VaIEfPEFjBwJy5dDeLinI/IKneuGMfmBRoTaHW9zhdqD/W4IPFhD3B9rE5GlJSjABo+10TxAeeF0e2rlypX55ZdfiIhwzD5/+eUXwvVBFRFx3qefQloa9O9vbTdtaj3EQee6YXSIDGVjzDESTp6lfIh128tfWn4uNqJrJM90rMVn637j92NnqFqmBL1bVFPLTx45nQA9/PDDDBkyhHPnznHTTTcB8OOPPzJs2DCXL4UhIlKonToFAwfCzJkQFAStWlmLmUqOigTYaFHD+WWRCqtiRQN4qHXhXv7EXZxOgIYNG8axY8d44oknSElJASA4OJjhw4czYsQIlwcoIlIo7dhhjfLauxcCAqx+PzVqeDoqEb/h9EzQGU6dOsWePXsoXrw4NWvWJCgoyNWxeYxmghYRtzEGpk6FQYOseX7Cw+Hzz6FtW09HJuLz3DoTdIb4+HiOHTtGjRo1CAoKIp95lIiI/zAG+vaFRx6xkp/Ona1RXkp+RAqc0wnQ0aNHad++Pddccw1du3bNHPn18MMPqw+QiEhubDaoWROKFIFx42DRImueHxEpcE4nQE8//TSBgYHExsZSosT5WSZ79OjBkiVLXBqciIjPMwaOHz+//fzzsGULDB9u9f0REY9wuhP00qVL+eGHH6hUqZLD/po1a/L777+7LDAREZ+XmGjd7tq3D9avh+LFrdaf+vU9HZkUoLR0o6H7XsjpBOj06dMOLT8Zjhw5Uqg6QouIXJbNm6FHDzh0yFrC4pdf4OabPR2VFLAlO+MYvXC3w2r2YfZgRnaL9LvJG72N0+2vbdq0YebMmZnbNpuN9PR0JkyYwI033ujS4EREfI4xMHEitGxpJT9Vq8KaNUp+/NCSnXEMmLXVIfkBiE88y4BZW1myU6sneJLTLUATJkygXbt2bN68mZSUFIYNG8auXbs4duwYv/zyiztiFBHxDcePW7M5f/ONtX3HHTBtGlx5pSejEg9ISzeMXrib7MZHG6w1zEYv3E2HyFDdDvMQp1uAIiMj2b59O02bNqVDhw6cPn2au+66i6ioKGpoEi8R8WdPPGElP8WKWa1A8+Yp+fFTG2OOZWn5uZAB4hLPsjHmWMEFJQ6cbgECCA0NZfTo0a6ORUTEt40fDwcPwuTJ0Lixp6MRD0o4mXPyk59y4npOtwAtWbKENWvWZG5PmjSJBg0a0KtXL45fONRTRKSwO3oUZsw4v12lCmzYoORHKB8SfOlCTpQT13M6Afq///s/kpKSANixYwdDhw6la9euHDp0iKFDh7o8QBERr/TLL9CgATz4ICxceH6/Tf05BJpGlCHMHkxO7wYb1miwphFlCjIsuYDTCVBMTAyRkZEAzJ07l27duvHaa6/xwQcf8P3337s8QBERr5Kebs3i3LYt/PGHNbNz5cqejkq8TJEAGyO7Wd+VFydBGdsju0WqA7QHOZ0AFStWjDNnzgCwfPlyOnbsCECZMmUyW4ZERAqlhATo2hVGjIC0NOjVy5rVuUEDT0cmXqhz3TAmP9CIULvjba5QezCTH2ikeYA8zOlO0K1atWLo0KHccMMNbNy4kS+++AKAX3/9Ncvs0CIihcaqVdCzJ8TFQXAwvP++NeRdt7wkF53rhtEhMlQzQXshpxOg999/nyeeeIKvv/6ayZMnU7FiRQC+//57Onfu7PIARUS8Qlyc9ahdG778EurW9XRE4iOKBNhoUaOsp8OQi9iMMdnN0+TXkpKSsNvtJCYmUqpUKU+HIyKeYoxjC8/MmXD33VCypOdiEpEcOfP9raWIRUSy8+OP0KgRxMef39enj5IfkUJCCZCIyIXS0uDll6FDB4iOBk36KlIo5WsmaBGRQumvv6yRXatWWdsPPwxvvunZmETELZQAiYgA/PADPPAAHDkCV1wBH35oJUMiUigpARKRfElLN4VnaO9XX0H37tbP9etbo7yuucazMf1PoXqdRfCe93S+EqCZM2dit9u5/fbbM/d9++23JCYm0qdPnzzXs3r1aiZMmMCWLVuIi4tj/vz53HHHHZnPG2MYPXo0H330EcePH6dZs2ZMmjSJOnXq5Frv3Llzeemllzh48CA1atTg1Vdf5c4773T6PEUke0t2xjF64W6H1a7D7MGM7Bbpm5O7de5sJTw332zd8gr2jvWZCt3rLH7Pm97T+eoE3a9fP0aMGOGwb/jw4Tz44INO1XP69Gnq16/P+++/n+3zr7/+Om+99Rbvv/8+mzZtIjQ0lA4dOnDy5Mkc61y3bh09evSgd+/ebNu2jd69e9O9e3c2bNjgVGwikr0lO+MYMGurwx8wgPjEswyYtZUlO+M8FJmT1q+3hrkDhITApk0waZJXJT+F4nUW+R9ve097zTxANpvNoQXIGEN4eDhDhgxh+PDhACQnJ1OhQgXGjx/PY489lm09PXr0ICkpyWFdss6dO3PllVcye/bsPMWieYBEspeWbmg1/qcsf8Ay2LCm+V8z/CbvvU2TkgLPP2+19Lz1Fjz9tKcjyqJQvM4iFyio93ShmAcoJiaG+Pj4zLXGAIKCgmjbti1r167N8bh169Y5HAPQqVOnXI9JTk4mKSnJ4SEiWW2MOZbjHzAAA8QlnmVjzLGCC8oZv/0GbdqcH9n1558eDScnPv86i1zEG9/T+UqADh48yIsvvkjPnj1JSEgAYMmSJezatctlgcX/b/KxChUqOOyvUKFC5nM5HefsMWPHjsVut2c+KmtlZ5FsJZzM+Q9YfsoVqG++gYYNYcMGKF0a5s+HN97wdFTZ8unXWSQb3viedjoBWrVqFfXq1WPDhg3MmzePU6dOAbB9+3ZGjhzp8gBtFy00aIzJsu9yjxkxYgSJiYmZj8OHD+c/YJFCrHxI3vrH5LVcgUhOhsGD4c474cQJaNYMoqLgggEX3sYnX2eRXHjje9rpBOi5555jzJgxLFu2jGLFimXuv/HGG1m3bp3LAgsNDQXI0nKTkJCQpYXn4uOcPSYoKIhSpUo5PEQkq6YRZQizB5PTfydsWCM6mkaUKciwcrd7N3zwgfXzM8/A6tVQrZpHQ7oUn3ydRXLhje9ppxOgHTt2ZDukvFy5chw9etQlQQFEREQQGhrKsmXLMvelpKSwatUqWrZsmeNxLVq0cDgGYOnSpbkeIyJ5UyTAxshukQBZ/pBlbI/sFuldHXMbNoT33oOFC61bXhf8x81b+eTrLJILb3xPO50AlS5dmri4rEPVoqKiqFixolN1nTp1iujoaKKjowGr43N0dDSxsbHYbDaGDBnCa6+9xvz589m5cyf9+vWjRIkS9LpgdtY+ffo4DMkfPHgwS5cuZfz48ezdu5fx48ezfPlyhgwZ4uypikg2OtcNY/IDjQi1OzZVh9qDmfxAI8/PT3P2rHXLa/v28/sefxxuvdVzMeWD17/OIk7yuve0cdL//d//mVatWpm4uDgTEhJi9u/fb9asWWOqV69uRo0a5VRdK1asMFidvx0effv2NcYYk56ebkaOHGlCQ0NNUFCQadOmjdmxY4dDHW3bts0sn+Grr74y1157rQkMDDS1atUyc+fOdSquxMREA5jExESnjhPxJ6lp6WbtgSPmm6g/zNoDR0xqWrqnQzJm3z5j6tc3BoypVcuYc+c8HdFl88rXWeQyuPM97cz3t9PzAJ07d45+/foxZ84cjDEULVqUtLQ0evXqxYwZMyhSpIjLk7SCpnmARHzQ55/DY4/BqVNQrhx89hl06uTpqESkADnz/Z3viRAPHTrE1q1bSU9Pp2HDhtSsWTNfwXojJUAiPuTMGeuW1yefWNtt21rJUHi4Z+MSkQLnzPd3vhdDrV69OtWrV8/v4SIily8+Hjp0gJ07wWaDl16yHkW1zrOI5M7pTtD33HMP48aNy7J/woQJ3HvvvS4JSkQkT8qVg/LloUIFWLYMRo9W8iMieZKviRBvueWWLPs7d+7M6tWrXRKUiEiOTp+2RnoBFCkC//0vREdD+/YeDUtEfIvTCdCpU6ccJkDMEBgYqDW0RMS9du6E6693XMA0NNR6iIg4wekEqG7dunzxxRdZ9s+ZM4fIyEiXBCXirdLSDesOHuXb6D9Zd/Aoaen5GkNQKBToa2EMTJ1qJT979sCCBeDCiVdFxP84fbP8pZde4u677+bgwYPcdNNNAPz444/Mnj2br776yuUBiniLJTvjGL1wt8OKxmH2YEZ2i/S7SekK9LU4eRIGDLBudYE1tP2zz6BsWdf+HhHxK/kaBr9o0SJee+01oqOjKV68ONdddx0jR46kbdu27oixwGkYvFxsyc44BszaysUfloxJ2/1pZt4CfS22bYPu3eHXX63+PmPGwLBhEOB047WI+IECmQeoMFMCJBdKSze0Gv+TQ2vHhWxYU7mvGX5ToV+bqUBfi+RkqF4d/voLKlWCOXPghhsur04RKdSc+f7O93+jUlJS+OOPP4iNjXV4iBQ2G2OO5fiFD9b6LXGJZ9kYc6zggvKQAn0tgoJg8mRrDa/oaCU/IuJSTvcB2r9/P/3792ft2rUO+40x2Gw20tLSXBaciDdIOJnzF35+yvkyt78WW7bA8eNw883W9m23Qbdu1iSHIiIu5HQC1K9fP4oWLcp3331HWFgYNv1hkkKufEjwpQs5Uc6Xue21MAbefx+efRauuMJq8alc2XpOf2NExA2cToCio6PZsmULtWrVckc8Il6naUQZwuzBxCeezdLxF873e2kaUaagQytwbnktjh+Hhx6C+fOt7TZtrCRIRMSNnO4DFBkZyZEjR9wRi4hXKhJgY2Q3a46ri9siMrZHdoss9B2gwQ2vxYYN0KiRlfwUKwYTJ8K8eXDllS6LWUQkO04nQOPHj2fYsGGsXLmSo0ePkpSU5PAQKYw61w1j8gONCLU73toJtQf71RB4cNFrYQy89Ra0agW//WaN9lq7FgYN0i0vESkQTg+DD/jf/BsX9/0pTJ2gNQxecpKWbtgYc4yEk2cpH2Ld6vGHlp/sXPZr8eij8PHHcO+91r92u/uCFRG/4Mz3t9N9gFasWJHvwER8XZEAGy1qaAZiyOdrkZ5+fhLDd9+Ftm2hVy+1+ohIgdNEiNlQC5CIi6Wnw4QJsGoVfPedZnIWEbdw+0SIP//8Mw888AAtW7bkzz//BOCzzz5jzZo1+alORAqzf/6BW26B556D77+Hb7/1dEQiIs4nQHPnzqVTp04UL16crVu3kpycDMDJkyd57bXXXB6giPiw1auhQQNYsgSCg+GTT+COOzwdlYiI8wnQmDFjmDJlCh9//DGBgYGZ+1u2bMnWrVtdGpyI+Ki0NGvh0htvtNbyql0bNm2y5vtRfx8R8QJOd4Let28fbdq0ybK/VKlSnDhxwhUxiYive+IJ+Ogj6+d+/axZnkuW9GhIIiIXcjoBCgsL48CBA1SrVs1h/5o1a6hevbqr4hLxShoGn0cDBsDXX8Pbb0OfPp6Oxqe58z3nrrr1ORFf4HQC9NhjjzF48GCmTZuGzWbjr7/+Yt26dTz77LO8/PLL7ohRxCss2RnH6IW7HVZDD7MHM7JbpF9NhJittDTYuBFatLC2GzSA33/XkhaXyZ3vOXfVrc+J+Ip8DYN/4YUXePvttzl71nqDBwUF8eyzz/LKK6+4PEBP0DB4udiSnXEMmLU1y/pXGf+n9bfZoB389Zc1l8/atfDLL3D99Z6OqFBw53vOXXXrcyKe5vZh8K+++ipHjhxh48aNrF+/nn/++afQJD8iF0tLN4xeuDvbxT8z9o1euJu0dD+cUuuHH6zWnlWrICjISobksrnzPeeuuvU5EV+T79nISpQoQZMmTWjatClXqJlbCrGNMcccmvMvZoC4xLNsjDlWcEF5WmoqjBgBnTtb8/zUrw9btsDtt3s6skLBne85d9Wtz4n4mjz1AbrrrrvyXOG8efPyHYyIN0o4mfMf9fyU83mHD0PPntbtLrBGfL35pjXPj7iEO99z7qpbnxPxNXlKgOwXLFJojGH+/PnY7XaaNGkCwJYtWzhx4oRTiZKIrygfkrcv9ryW83nz5lnJT6lS1sSG997r6YgKHXe+59xVtz4n4mvylABNnz498+fhw4fTvXt3pkyZQpEiRQBIS0vjiSeeUIdhKZSaRpQhzB5MfOLZbPs32IBQuzXU1y8MGmT19Xn0UahRw9PRFErufM+5q259TsTXON0HaNq0aTz77LOZyQ9AkSJFGDp0KNOmTXNpcCLeoEiAjZHdIoHzo1kyZGyP7BZZeOc5+f13ay6fU6es7YAAGD9eyY8bufM95666/f5zIj7H6QQoNTWVPXv2ZNm/Z88e0tPTXRLUhapVq4bNZsvyGDhwYLblV65cmW35vXv3ujw28R+d64Yx+YFGhNodm+9D7cGFe2jvt99ao7w++wyGDfN0NH7Fne85d9Xtt58T8UlOT4T44IMP0r9/fw4cOEDz5s0BWL9+PePGjePBBx90eYCbNm0iLS0tc3vnzp106NCBey/R72Dfvn0Ot+TKlSvn8tjEv3SuG0aHyFD/mOE2JcVKeN5919pu2lQJkAe48z3nrrr96nMiPs3pBOiNN94gNDSUt99+m7i4OMBaHmPYsGE888wzLg/w4sRl3Lhx1KhRg7Zt2+Z6XPny5SldurTL4xH/ViTARosaZT0dhnsdOgQ9esDmzdb2M8/Aa69BsWKejctPufM95666/eJzIj7P6QQoICCAYcOGMWzYMJKSkgAKrPNzSkoKs2bNYujQodgusaJ0w4YNOXv2LJGRkbz44ovceOONOZZNTk4mOTk5czvjvET8zsqV1lw+SUlQpgx8+incequnoxIRcbl8T4QIVuJTkCO/vvnmG06cOEG/fv1yLBMWFsZHH33E3LlzmTdvHtdeey3t27dn9erVOR4zduxY7HZ75qNy5cpuiF7EB1x7rTWfzw03QHS0kh8RKbScXgvs77//5tlnn+XHH38kISGBiw+/sL+Oq3Xq1IlixYqxcOFCp47r1q0bNpuNBQsWZPt8di1AlStX1lpg4h+OHIGrrjq/vXevNcIrMNBzMYmI5IMza4E5fQusX79+xMbG8tJLLxEWFnbJW1Gu8vvvv7N8+fJ8zTTdvHlzZs2alePzQUFBBAUFXU54Ir5p9mx47DGYNg3uucfaV6uWZ2MSESkATidAa9as4eeff6ZBgwZuCCdn06dPp3z58txyyy1OHxsVFUVYmIZf5kdautFojgsUmtfj339h8GD4+GNre+bM8wmQiIgfcDoBqly5cpbbXu6Wnp7O9OnT6du3L0WLOoY8YsQI/vzzT2bOnAnAO++8Q7Vq1ahTp05mp+m5c+cyd+7cAo25MFiyM47RC3c7LHAYZg9mZLdIv5zPo9C8Hnv3QvfusGMH2Gzw4ovw8suejkpEpEA53Qn6nXfe4bnnnuO3335zQzjZW758ObGxsfTv3z/Lc3FxccTGxmZup6Sk8Oyzz3LdddfRunVr1qxZw6JFi7ROmZOW7IxjwKytWVZ3jk88y4BZW1myM85DkXlGoXk9Zs6Exo2t5KdCBVi6FP7zHyjq9P+FRER8mtOdoK+88krOnDlDamoqJUqUIPCijpLHjh1zaYCe4EwnqsIoLd3QavxPWb7sM2Ss6bNm+E2+efvHSYXm9di61Up+AG66Cf77XwgN9WxMIiIu5NZO0O+8805+4xIfsTHmWI5f9gAGiEs8y8aYY34x2VmheT0aNbImNbTb4fnn4YL1/ERE/I3TCVDfvn3dEYd4kYSTOX/Z56ecr/PZ18MY65ZX+/ZQqZK17403PBuTiIiXyNdEiAcPHuTFF1+kZ8+eJCQkALBkyRJ27drl0uDEM8qHBF+6kBPlfJ1Pvh4nT0Lv3tCvH/TsCampno5IRMSrOJ0ArVq1inr16rFhwwbmzZvHqVOnANi+fTsjR450eYBS8JpGlCHMHkxOvVlsWKOfmkaUKciwPMbnXo9t26BJE6uPT5EicMstEHBZk76LiBQ6Tv9VfO655xgzZgzLli2j2AWLI954442sW7fOpcGJZxQJsDGyWyRAli/9jO2R3SK9u8OvC/nM62EMfPghNGsGv/5q3fZatQqee04JkIjIRZz+q7hjxw7uvPPOLPvLlSvH0aNHXRKUeF7numFMfqARoXbH2zqh9mAmP9DIt+a9cQGvfz1OnoT77oPHH4fkZGsNr+hoa00vERHJwulO0KVLlyYuLo6IiAiH/VFRUVSsWNFlgYnnda4bRofI0MIx87ELePXrUaQI7N5tzeczbhwMHWpNcigiItlyOgHq1asXw4cP56uvvsJms5Gens4vv/zCs88+S58+fdwRo3hQkQCbdw/tLmBe9XoYYz0CAqBECfjyS0hMhObNPR2ZiIjXc/oW2KuvvkqVKlWoWLEip06dIjIykjZt2tCyZUtefPFFd8QoIhc7ccJau2v8+PP7atdW8iMikkdOzwSd4eDBg0RFRZGenk7Dhg2pWbOmq2PzGH+fCVq83MaN0KMH/PYbFC8OMTHWshYiIn7OrTNBZ6hRowbVq1cHwKa+BiLuZwy88w4MHw7nzkH16vDFF0p+RETyIV9jY6dOnUrdunUJDg4mODiYunXr8sknn7g6NhHJcOwY3H671bn53Dnr9tfWrdZ8PyIi4jSnW4Beeukl3n77bQYNGkSLFi0AWLduHU8//TS//fYbY8aMcXmQIn4tJcXq27N/PwQFwdtvW8Pd1fIqIpJvTvcBuuqqq3jvvffo2bOnw/7Zs2czaNAgjhw54tIAPUF9gMTrfPCBdfvryy+hQQNPRyMi4pWc+f52+hZYWloaTbJpdm/cuDGpWm9IxDWOHLHm9ckwYIA1saGSHxERl3A6AXrggQeYPHlylv0fffQR999/v0uCEvFrP/8M9etDt27WvD5g3e4qUcKzcYmIFCL5GgU2depUli5dSvP/zTmyfv16Dh8+TJ8+fRg6dGhmubfeess1UYr4g/R0GDsWXn7Z+rlWLfjnH7DbPR2ZiEih43QCtHPnTho1agRYcwGBtQ5YuXLl2LlzZ2Y5DY0XccLff0Pv3rBsmbXdty9MmgQlS3o2LhGRQsrpBGjFihXuiEPEf/30E9x/P8THW7e5PvjASoBERMRt8jUPEMCBAwf44Ycf+PfffwHI54TSIvL221byU6cObNqk5EdEpAA4nQAdPXqU9u3bc80119C1a1fi4uIAePjhh3nmmWdcHqBIoTd9Ojz7rLXERWSkp6MREfELTidATz/9NIGBgcTGxlLiglEpPXr0YMmSJS4NTqRQWrrUSngyXHUVTJigUV4iIgXI6T5AS5cu5YcffqBSpUoO+2vWrMnvv//ussBECp3UVBg50hrpZQy0bAl33eXpqERE/JLTCdDp06cdWn4yHDlyhKCgIJcEJVLo/PEH9OplzfED1lIWXbp4NiYRET/m9C2wNm3aMHPmzMxtm81Geno6EyZM4MYbb3RpcCKFwuLF1gzOP/8MISHWCu6TJ0Px4p6OTETEbzndAjRhwgTatWvH5s2bSUlJYdiwYezatYtjx47xyy+/uCNGEd/12mvwwgvWz40bW8lPjRqejUlERJxvAYqMjGT79u00bdqUDh06cPr0ae666y6ioqKooT/sIo4aN7aWsRg0CH75RcmPiIiXcHo1eH+g1eDlsiQkQPny57f37IHatT0Xj4iIn3Dm+ztPt8C2b9+e519+3XXX5bmsSKGSkgLDh8OMGbBlC1Svbu1X8iMi4nXylAA1aNAAm82GMcZhja+MxqML96Wlpbk4RBEfEBMDPXpYMzkDfP89DBzo2ZhERCRHeeoDFBMTw6FDh4iJiWHu3LlERETwwQcfEB0dTXR0NB988AE1atRg7ty57o5XxPvMnQsNG1rJT5kysGCBkh8RES+XpwSoatWqmY/XXnuNiRMn8thjj3Hddddx3XXX8dhjj/HOO+/wyiuvuDS4UaNGYbPZHB6hoaG5HrNq1SoaN25McHAw1atXZ8qUKS6NSSTT2bPw5JNwzz2QmGhNbBgVBd26eToyERG5BKdHge3YsYOIiIgs+yMiIti9e7dLgrpQnTp1iIuLy3zs2LEjx7IxMTF07dqV1q1bExUVxfPPP89TTz2llilxj4kTYdIk6+fhw2HlSqhSxaMhiYhI3jg9D1Dt2rUZM2YMU6dOJTg4GIDk5GTGjBlDbTd09ixatOglW30yTJkyhSpVqvDOO+9kxrp582beeOMN7r77bpfHJn5u8GBYsQKeekqzOouI+BinE6ApU6bQrVs3KleuTP369QHYtm0bNpuN7777zuUB7t+/n/DwcIKCgmjWrBmvvfYa1TNG11xk3bp1dOzY0WFfp06dmDp1KufOnSMwMDDb45KTk0lOTs7cTkpKct0JSOHx779Wi8+QIVC0KAQFWZ2dRUTE5zidADVt2pSYmBhmzZrF3r17McbQo0cPevXqRcmSJV0aXLNmzZg5cybXXHMNf//9N2PGjKFly5bs2rWLsmXLZikfHx9PhQoVHPZVqFCB1NRUjhw5QlhYWLa/Z+zYsYwePdqlsUshs3cvdO8OO3bAiRMwZoynIxIRkcvgdAIEUKJECR599FFXx5JFlwtuK9SrV48WLVpQo0YNPv30U4YOHZrtMRcOyYfsh+pfbMSIEQ71JSUlUbly5csJXQqTzz6DAQPg9GmoUAHatfN0RCIicpnylQD9+uuvrFy5koSEBNLT0x2ee/nll10SWHZKlixJvXr12L9/f7bPh4aGEh8f77AvISGBokWLZttilCEoKEgr2UtWp09bS1hMn25t33QT/Pe/kMc+aSIi4r2cToA+/vhjBgwYwFVXXUVoaKhDy4rNZnNrApScnMyePXto3bp1ts+3aNGChQsXOuxbunQpTZo0ybH/j0i29uyxhrfv3g0BATBypLWoaZEino5MRERcwOkEaMyYMbz66qsMHz7cHfE4ePbZZ+nWrRtVqlQhISGBMWPGkJSURN++fQHr1tWff/7JzJkzAXj88cd5//33GTp0KI888gjr1q1j6tSpzJ492+2xSiGTnm7N7hwWBp9/rtteIiKFjNMJ0PHjx7n33nvdEUsWf/zxBz179uTIkSOUK1eO5s2bs379eqpWrQpAXFwcsbGxmeUjIiJYvHgxTz/9NJMmTSI8PJyJEydqCLzkTVra+RaeOnVg/nxrhucLFzYVEZFCwenV4B966CGuv/56Hn/8cXfF5HFaDd4PbdsGvXrBhx9Cq1aejkZERPLB5avBX+jqq6/mpZdeYv369dSrVy9L35qnnnrK2SpFPMcY+Ogja1LD5GT4v/+DtWshl1GDIiLi+5xuAcpuGYzMymw2Dh06dNlBeZpagPxEUhI8+ih88YW13bUrfPopXHWVZ+MSEZF8cWsLUExMTL4DE/EaW7dCjx5w4IA1q/PYsTB0qDXiS0RECr18zQMk4tN27oQWLSAlxVq8dM4ca1tERPxGnhKgoUOH8sorr1CyZMkcZ2DO8NZbb7kkMBG3qVMHbr0VUlOtSQ7LlPF0RCIiUsDylABFRUVx7ty5zJ9zkttyEyIetXkz1KwJdrvVwXnWLAgOVmdnERE/5XQnaH+gTtCFiDHwzjswfDjcead1u0tJj4hIoeTWTtAiPuPYMXjwQViwwNpOT7f6/WjdNxERv6chL1I4rVsHDRpYyU+xYjBpEnz5pZIfEREBlABJYZOeDq+/Dq1bw+HDcPXVsH49PPGEbn2JiEgmJUBSuJw4Ae++a63r1bOnNd9Pw4aejkpERLyM+gBJ4VKmDMyeDfv2wcMPq9VHRESypQRIfFt6ujWLc9Wq8MAD1r42bayHiIhIDpQAie/6+2/o3RuWLYMSJeDGG6FiRU9HJSIiPkAJkPimFSugVy+Ij4fixeH99yE83NNRiYiIj1AnaPEtaWkwejTcfLOV/NSpY83y/OCD6u8jIiJ5phYg8R2pqdC5M/z4o7X90EMwcaJ1+0tERMQJagES31G0KFx/PZQsaa3l9cknSn5ERCRftBZYNrQWmBdJTYXjx6FcOWv73DmIjYUaNTwbl4iIeB1nvr/VAiTe648/rJFdt9xireEFEBio5EdERC6bEiDxTosXW2t5rVkDe/fCzp2ejkhERAoRJUDiXc6dg2HDrFafo0ehUSNrOYtGjTwdmYiIFCIaBSbe4/ff4b77rMVLAQYNggkTtIK7iIi4nBIg8R4PP2wlP3Y7TJsGd93l6YhERKSQ0i0w8R6TJ1sTHEZFKfkRERG3UgIknhMTY83lk+Hqq611vSIiPBeTiIj4Bd0CE8+YO9eayTkpCapVs1p+RERECohagKRgnT0LTz4J99wDiYnQvDnUrOnpqERExM8oAZKCc+AAtGwJkyZZ28OGwapVULWqZ+MSERG/o1tgUjC++sq65XXyJJQtCzNnQteuno5KRET8lBIgKRinTlnJT+vW8PnnUKmSpyMSERE/5tW3wMaOHcv1119PSEgI5cuX54477mDfvn25HrNy5UpsNluWx969ewsoasmUmnr+53794Msv4aeflPyIiIjHeXUCtGrVKgYOHMj69etZtmwZqampdOzYkdOnT1/y2H379hEXF5f5qKmOtgXrs8/guuus5SwAbDa4914oqkZHERHxPK/+NlqyZInD9vTp0ylfvjxbtmyhTZs2uR5bvnx5Spcu7cboJFunT1tLWEyfbm1PnAijR3s2JhERkYt4dQvQxRITEwEoU6bMJcs2bNiQsLAw2rdvz4oVK3Itm5ycTFJSksND8mHXLmja1Ep+bDYYNQpeftnTUYmIiGThMwmQMYahQ4fSqlUr6tatm2O5sLAwPvroI+bOncu8efO49tprad++PatXr87xmLFjx2K32zMflStXdscpFF7GWEnP9dfD7t0QGgo//ggjR0KRIp6OTkREJAubMcZ4Ooi8GDhwIIsWLWLNmjVUcrITbbdu3bDZbCxYsCDb55OTk0lOTs7cTkpKonLlyiQmJlKqVKnLitsvTJpkTW4I0KGD1f+nQgXPxiQiIn4nKSkJu92ep+9vn2gBGjRoEAsWLGDFihVOJz8AzZs3Z//+/Tk+HxQURKlSpRwe4oT777fW8Xr1VViyRMmPiIh4Pa/uBG2MYdCgQcyfP5+VK1cSkc9FMqOioggLC3NxdH7MGFi+3Fq/y2aD0qVhxw4IDvZ0ZCIiInni1QnQwIED+fzzz/n2228JCQkhPj4eALvdTvHixQEYMWIEf/75JzNnzgTgnXfeoVq1atSpU4eUlBRmzZrF3LlzmTt3rsfOo1BJSoLHHoM5c+DDD+HRR639Sn5ERMSHeHUCNHnyZADatWvnsH/69On069cPgLi4OGJjYzOfS0lJ4dlnn+XPP/+kePHi1KlTh0WLFtFVyy5cvqgo6N7dWtOraFH4919PRyQiIpIvPtMJuiA504nKLxgDH3wAQ4dCSgpUqWK1ALVo4enIREREMjnz/e3VLUDiBU6cgIcfhoxbiLfdZg15z8NcTCIiIt7KJ0aBiQft2AHz50NgILz9NnzzjZIfERHxeWoBkty1bg3vvw9NmlgTHYqIiBQCagESR8eOQa9esG/f+X0DBij5ERGRQkUtQHLeunVw330QG2uN9NqwwZrnR0REpJBRC5BAejpMmABt2ljJT40aMGWKkh8RESm01ALk744cgb59YfFia7tHD/joI9DwfxERKcSUAPmzAwegXTv4809rJud334VHHlHLj4iIFHpKgPxZ1arW44or4Msv4brrPB2RiIhIgVAC5G/++QfsdihWzJrb5+uvISTESoJERET8hDpB+5MVK6xWnuefP78vLEzJj4iI+B0lQP4gLQ1Gj4abb4b4eFiyBM6c8XRUIiIiHqMEqLCLi4OOHWHUKGu4e//+sHEjlCjh6chEREQ8Rn2ACrNly+CBByAhAUqWhMmToXdvT0clIiLicUqACqsTJ+DeeyExEerVs0Z51arl6ahERES8ghKgwqp0aWs25xUr4J13oHhxT0ckIiLiNWzGGOPpILxNUlISdrudxMRESvnSjMjff29NaHjjjZ6OREREpMA58/2tTtCFwblzMHw4dO0KPXvC3397OiIRERGvpltgvi421lrBfd06a/uee6yJDkVERCRHSoB82YIF0K8fHD9uJT1Tp8Ldd3s6KhEREa+nW2C+KC0Nhg6F22+3kp/rr4etW5X8iIiI5JESIF8UEGDN7QMwZAisWQPVq3s0JBEREV+iW2C+JDUVihYFm82a1PD++6FLF09HJSIi4nPUAuQLkpNh0CDrFlfGrAUhIUp+RERE8kktQN7uwAHo0cPq4wPW7a7WrT0bk4iIiI9TC5A3++ILaNTISn7KloXvvlPyIyIi4gJKgLzRv//C449b8/ucPAmtWkF0NNxyi6cjExERKRSUAHmj++6DDz+0Ojs//7y1nlelSp6OSkREpNBQHyBv9PzzsGULTJsGHTt6OhoREZFCRwmQNzhzBjZtgrZtre1mzeDgQQgK8mxcIiIihZRugXna7t3QtCl07gzbt5/fr+RHRETEbXwiAfrggw+IiIggODiYxo0b8/PPP+daftWqVTRu3Jjg4GCqV6/OlClTCihSJxgD06dDkyawaxeULg1JSZ6OSkRExC94fQL0xRdfMGTIEF544QWioqJo3bo1Xbp0ITY2NtvyMTExdO3aldatWxMVFcXzzz/PU089xdy5cws48lycOgV9+0L//taIrw4drFFerVp5OjIRERG/YDMmY2ph79SsWTMaNWrE5MmTM/fVrl2bO+64g7Fjx2YpP3z4cBYsWMCePXsy9z3++ONs27aNdevW5el3JiUlYbfbSUxMpFSpUpd/Ehfavt2a2HDvXmtNr//8B0aMsH4WERGRfHPm+9urv3VTUlLYsmULHS8aCdWxY0fWrl2b7THr1q3LUr5Tp05s3ryZc+fOZXtMcnIySUlJDg+3+fZbK/kJD7eGt7/wgpIfERGRAubV37xHjhwhLS2NChUqOOyvUKEC8fHx2R4THx+fbfnU1FSOHDmS7TFjx47FbrdnPipXruyaE8jO88/Diy9at7zatHHf7xEREZEceXUClMFmszlsG2Oy7LtU+ez2ZxgxYgSJiYmZj8OHD19mxLkoUgReeQXKlXPf7xAREZFcefU8QFdddRVFihTJ0tqTkJCQpZUnQ2hoaLblixYtStmyZbM9JigoiCANOxcREfEbXt0CVKxYMRo3bsyyZcsc9i9btoyWLVtme0yLFi2ylF+6dClNmjQhMDDQbbGKiIiI7/DqBAhg6NChfPLJJ0ybNo09e/bw9NNPExsby+OPPw5Yt6/69OmTWf7xxx/n999/Z+jQoezZs4dp06YxdepUnn32WU+dgoiIiHgZr74FBtCjRw+OHj3Kf/7zH+Li4qhbty6LFy+matWqAMTFxTnMCRQREcHixYt5+umnmTRpEuHh4UycOJG7777bU6cgIiIiXsbr5wHyBLfOAyQiIiJuUWjmARIRERFxByVAIiIi4neUAImIiIjfUQIkIiIifkcJkIiIiPgdJUAiIiLid5QAiYiIiN9RAiQiIiJ+RwmQiIiI+B2vXwrDEzImx05KSvJwJCIiIpJXGd/beVnkQglQNk6ePAlA5cqVPRyJiIiIOOvkyZPY7fZcy2gtsGykp6fz119/ERISgs1mc2ndSUlJVK5cmcOHDxfKdcYK+/lB4T9HnZ/vK+znqPPzfe46R2MMJ0+eJDw8nICA3Hv5qAUoGwEBAVSqVMmtv6NUqVKF9o0Nhf/8oPCfo87P9xX2c9T5+T53nOOlWn4yqBO0iIiI+B0lQCIiIuJ3lAAVsKCgIEaOHElQUJCnQ3GLwn5+UPjPUefn+wr7Oer8fJ83nKM6QYuIiIjfUQuQiIiI+B0lQCIiIuJ3lACJiIiI31ECJCIiIn5HCZAbfPDBB0RERBAcHEzjxo35+eefcy2/atUqGjduTHBwMNWrV2fKlCkFFKlzxo4dy/XXX09ISAjly5fnjjvuYN++fbkes3LlSmw2W5bH3r17Cyhq54waNSpLrKGhobke4yvXD6BatWrZXo+BAwdmW97br9/q1avp1q0b4eHh2Gw2vvnmG4fnjTGMGjWK8PBwihcvTrt27di1a9cl6507dy6RkZEEBQURGRnJ/Pnz3XQGl5bbOZ47d47hw4dTr149SpYsSXh4OH369OGvv/7Ktc4ZM2Zke13Pnj3r5rPJ6lLXsF+/flnibN68+SXr9ZZreKnzy+462Gw2JkyYkGOd3nT98vK94K2fQyVALvbFF18wZMgQXnjhBaKiomjdujVdunQhNjY22/IxMTF07dqV1q1bExUVxfPPP89TTz3F3LlzCzjyS1u1ahUDBw5k/fr1LFu2jNTUVDp27Mjp06cveey+ffuIi4vLfNSsWbMAIs6fOnXqOMS6Y8eOHMv60vUD2LRpk8O5LVu2DIB777031+O89fqdPn2a+vXr8/7772f7/Ouvv85bb73F+++/z6ZNmwgNDaVDhw6Z6/1lZ926dfTo0YPevXuzbds2evfuTffu3dmwYYO7TiNXuZ3jmTNn2Lp1Ky+99BJbt25l3rx5/Prrr9x2222XrLdUqVIO1zQuLo7g4GB3nEKuLnUNATp37uwQ5+LFi3Ot05uu4aXO7+JrMG3aNGw2G3fffXeu9XrL9cvL94LXfg6NuFTTpk3N448/7rCvVq1a5rnnnsu2/LBhw0ytWrUc9j322GOmefPmbovRVRISEgxgVq1alWOZFStWGMAcP3684AK7DCNHjjT169fPc3lfvn7GGDN48GBTo0YNk56enu3zvnT9ADN//vzM7fT0dBMaGmrGjRuXue/s2bPGbrebKVOm5FhP9+7dTefOnR32derUydx3330uj9lZF59jdjZu3GgA8/vvv+dYZvr06cZut7s2OBfI7vz69u1rbr/9dqfq8dZrmJfrd/vtt5ubbrop1zLeev2Myfq94M2fQ7UAuVBKSgpbtmyhY8eODvs7duzI2rVrsz1m3bp1Wcp36tSJzZs3c+7cObfF6gqJiYkAlClT5pJlGzZsSFhYGO3bt2fFihXuDu2y7N+/n/DwcCIiIrjvvvs4dOhQjmV9+fqlpKQwa9Ys+vfvf8lFf33p+mWIiYkhPj7e4foEBQXRtm3bHD+PkPM1ze0Yb5KYmIjNZqN06dK5ljt16hRVq1alUqVK3HrrrURFRRVMgPmwcuVKypcvzzXXXMMjjzxCQkJCruV99Rr+/fffLFq0iIceeuiSZb31+l38veDNn0MlQC505MgR0tLSqFChgsP+ChUqEB8fn+0x8fHx2ZZPTU3lyJEjbov1chljGDp0KK1ataJu3bo5lgsLC+Ojjz5i7ty5zJs3j2uvvZb27duzevXqAow275o1a8bMmTP54Ycf+Pjjj4mPj6dly5YcPXo02/K+ev0AvvnmG06cOEG/fv1yLONr1+9CGZ85Zz6PGcc5e4y3OHv2LM899xy9evXKdYHJWrVqMWPGDBYsWMDs2bMJDg7mhhtuYP/+/QUYbd506dKF//73v/z000+8+eabbNq0iZtuuonk5OQcj/HVa/jpp58SEhLCXXfdlWs5b71+2X0vePPnUKvBu8HF/5s2xuT6P+zsyme335s8+eSTbN++nTVr1uRa7tprr+Xaa6/N3G7RogWHDx/mjTfeoE2bNu4O02ldunTJ/LlevXq0aNGCGjVq8OmnnzJ06NBsj/HF6wcwdepUunTpQnh4eI5lfO36ZcfZz2N+j/G0c+fOcd9995Gens4HH3yQa9nmzZs7dCS+4YYbaNSoEe+99x4TJ050d6hO6dGjR+bPdevWpUmTJlStWpVFixblmij44jWcNm0a999//yX78njr9cvte8EbP4dqAXKhq666iiJFimTJUBMSErJkshlCQ0OzLV+0aFHKli3rtlgvx6BBg1iwYAErVqygUqVKTh/fvHlzj/9PJa9KlixJvXr1cozXF68fwO+//87y5ct5+OGHnT7WV65fxug9Zz6PGcc5e4ynnTt3ju7duxMTE8OyZctybf3JTkBAANdff71PXNewsDCqVq2aa6y+eA1//vln9u3bl6/PpDdcv5y+F7z5c6gEyIWKFStG48aNM0fWZFi2bBktW7bM9pgWLVpkKb906VKaNGlCYGCg22LND2MMTz75JPPmzeOnn34iIiIiX/VERUURFhbm4ujcIzk5mT179uQYry9dvwtNnz6d8uXLc8sttzh9rK9cv4iICEJDQx2uT0pKCqtWrcrx8wg5X9PcjvGkjORn//79LF++PF+JtzGG6Ohon7iuR48e5fDhw7nG6mvXEKwW2caNG1O/fn2nj/Xk9bvU94JXfw5d1p1ajDHGzJkzxwQGBpqpU6ea3bt3myFDhpiSJUua3377zRhjzHPPPWd69+6dWf7QoUOmRIkS5umnnza7d+82U6dONYGBgebrr7/21CnkaMCAAcZut5uVK1eauLi4zMeZM2cyy1x8fm+//baZP3+++fXXX83OnTvNc889ZwAzd+5cT5zCJT3zzDNm5cqV5tChQ2b9+vXm1ltvNSEhIYXi+mVIS0szVapUMcOHD8/ynK9dv5MnT5qoqCgTFRVlAPPWW2+ZqKiozBFQ48aNM3a73cybN8/s2LHD9OzZ04SFhZmkpKTMOnr37u0wSvOXX34xRYoUMePGjTN79uwx48aNM0WLFjXr168v8PMzJvdzPHfunLnttttMpUqVTHR0tMPnMjk5ObOOi89x1KhRZsmSJebgwYMmKirKPPjgg6Zo0aJmw4YNXnV+J0+eNM8884xZu3atiYmJMStWrDAtWrQwFStW9JlreKn3qDHGJCYmmhIlSpjJkydnW4c3X7+8fC946+dQCZAbTJo0yVStWtUUK1bMNGrUyGGYeN++fU3btm0dyq9cudI0bNjQFCtWzFSrVi3HD4GnAdk+pk+fnlnm4vMbP368qVGjhgkODjZXXnmladWqlVm0aFHBB59HPXr0MGFhYSYwMNCEh4ebu+66y+zatSvzeV++fhl++OEHA5h9+/Zlec7Xrl/GMP2LH3379jXGWENwR44caUJDQ01QUJBp06aN2bFjh0Mdbdu2zSyf4auvvjLXXnutCQwMNLVq1fJowpfbOcbExOT4uVyxYkVmHRef45AhQ0yVKlVMsWLFTLly5UzHjh3N2rVrC/7kTO7nd+bMGdOxY0dTrlw5ExgYaKpUqWL69u1rYmNjHerw5mt4qfeoMcZ8+OGHpnjx4ubEiRPZ1uHN1y8v3wve+jm0/e8ERERERPyG+gCJiIiI31ECJCIiIn5HCZCIiIj4HSVAIiIi4neUAImIiIjfUQIkIiIifkcJkIiIiPgdJUAiIiLid5QAifiAdu3aMWTIkMuq47fffsNmsxEdHe2SmNytX79+3HHHHZ4OI1+++eYbrr76aooUKcKQIUOYMWMGpUuXLpDf7Yr3iog/KOrpAETk0ubNm+fVi6uKo8cee4wHH3yQp556ipCQEIoWLUrXrl0vq85+/fpx4sQJvvnmG9cEKeLnlACJ+IAyZcp4OgTJo1OnTpGQkECnTp0IDw/P3F+8ePEcjzl37pxXJ7hpaWnYbDYCAnTTQAoPvZtFfMDFtzWqVavGa6+9Rv/+/QkJCaFKlSp89NFHDsds3LiRhg0bEhwcTJMmTYiKispS7+7du+natStXXHEFFSpUoHfv3hw5csTh9z755JM8+eSTlC5dmrJly/Liiy9y4RKCKSkpDBs2jIoVK1KyZEmaNWvGypUrM5/PuP3zww8/ULt2ba644go6d+5MXFxcZpm0tDSGDh2a+TuGDRvGxcsUGmN4/fXXqV69OsWLF6d+/fp8/fXXmc+vXLkSm83Gjz/+SJMmTShRogQtW7Zk3759DvUsWLCAJk2aEBwczFVXXcVdd92V53O5lJUrVxISEgLATTfdhM1mY+XKlVlugY0aNYoGDRowbdo0qlevTlBQEMYYvv76a+rVq0fx4sUpW7YsN998M6dPn2bUqFF8+umnfPvtt9hstsx6c5KamprrNTt+/Dh9+vThyiuvpESJEnTp0oX9+/dnPp8R73fffUdkZCRBQUH8/vvveXrfifgMly6tKiJu0bZtWzN48ODM7apVq5oyZcqYSZMmmf3795uxY8eagIAAs2fPHmOMMadOnTLlypUzPXr0MDt37jQLFy401atXN4CJiooyxhjz119/mauuusqMGDHC7Nmzx2zdutV06NDB3HjjjQ6/94orrjCDBw82e/fuNbNmzTIlSpQwH330UWaZXr16mZYtW5rVq1ebAwcOmAkTJpigoCDz66+/GmOMmT59ugkMDDQ333yz2bRpk9myZYupXbu26dWrV2Yd48ePN3a73Xz99ddm9+7d5qGHHjIhISHm9ttvzyzz/PPPm1q1apklS5aYgwcPmunTp5ugoCCzcuVKY8z5VbebNWtmVq5caXbt2mVat25tWrZsmVnHd999Z4oUKWJefvlls3v3bhMdHW1effXVPJ+LMSbLStcXSk5ONvv27TOAmTt3romLizPJyclm+vTpxm63Z5YbOXKkKVmypOnUqZPZunWr2bZtm/nrr79M0aJFzVtvvWViYmLM9u3bzaRJk8zJkyfNyZMnTffu3U3nzp1NXFxcZr05vVcudc1uu+02U7t2bbN69WoTHR1tOnXqZK6++mqTkpLicM1atmxpfvnlF7N3715z6tSpS77vRHyJEiARH5BdAvTAAw9kbqenp5vy5cubyZMnG2OM+fDDD02ZMmXM6dOnM8tMnjzZIQF66aWXTMeOHR1+z+HDhw1g9u3bl/l7a9eubdLT0zPLDB8+3NSuXdsYY8yBAweMzWYzf/75p0M97du3NyNGjDDGWF+mgDlw4EDm85MmTTIVKlTI3A4LCzPjxo3L3D537pypVKlSZgJ06tQpExwcbNauXevwex566CHTs2dPY8z5BGj58uWZzy9atMgA5t9//zXGGNOiRQtz//33m+zk5VyMMebaa6818+bNy7YOY4w5fvy4AcyKFSsy92WXAAUGBpqEhITMfVu2bDGA+e2337Ktt2/fvg4JYU4udc1+/fVXA5hffvkl8/kjR46Y4sWLmy+//DIzXsBER0c71H2p952IL1EfIBEfdd1112X+bLPZCA0NJSEhAYA9e/ZQv359SpQokVmmRYsWDsdv2bKFFStWcMUVV2Sp++DBg1xzzTUANG/eHJvN5lDPm2++SVpaGlu3bsUYk1k2Q3JyMmXLls3cLlGiBDVq1MjcDgsLy4w1MTGRuLg4h/iKFi1KkyZNMm/b7N69m7Nnz9KhQweH35OSkkLDhg1zfF3CwsIASEhIoEqVKkRHR/PII49kOV8gz+eyd+/ebI93VtWqVSlXrlzmdv369Wnfvj316tWjU6dOdOzYkXvuuYcrr7zS6bpzu2Z79uyhaNGiNGvWLPP5smXLcu2117Jnz57MfcWKFXN4LTPk9r4T8SVKgER81MWdZm02G+np6QBZ+s9kJz09nW7dujF+/Pgsz2UkDnmpo0iRImzZsoUiRYo4PHdhYpVdrHmJ8cLfA7Bo0SIqVqzo8FxQUJDD9oW/KyMJyDg+t47IeT0XVylZsqTDdpEiRVi2bBlr165l6dKlvPfee7zwwgts2LCBiIgIl/3enF53Y4xD0lS8eHGH7Qy5ve9EfIk6QYsUQpGRkWzbto1///03c9/69esdyjRq1Ihdu3ZRrVo1rr76aofHhV/OFx+3fv16atasSZEiRWjYsCFpaWkkJCRkqSM0NDRPsdrtdsLCwhx+T2pqKlu2bHE4n6CgIGJjY7P8nsqVK+f5dbnuuuv48ccfs33OFedyuWw2GzfccAOjR48mKiqKYsWKMX/+fMBqkUlLS8tTPblds8jISFJTU9mwYUPm80ePHuXXX3+ldu3arjsZES+nBEikEOrVqxcBAQE89NBD7N69m8WLF/PGG284lBk4cCDHjh2jZ8+ebNy4kUOHDrF06VL69+/v8EV7+PBhhg4dyr59+5g9ezbvvfcegwcPBuCaa67h/vvvp0+fPsybN4+YmBg2bdrE+PHjWbx4cZ7jHTx4MOPGjWP+/Pns3buXJ554ghMnTmQ+HxISwrPPPsvTTz/Np59+ysGDB4mKimLSpEl8+umnef49I0eOZPbs2YwcOZI9e/awY8cOXn/9dafOpVatWplJiStt2LCB1157jc2bNxMbG8u8efP4559/MpOSatWqsX37dvbt28eRI0c4d+5cjnXlds1q1qzJ7bffziOPPMKaNWvYtm0bDzzwABUrVuT22293+XmJeCvdAhMphK644goWLlzI448/TsOGDYmMjGT8+PHcfffdmWXCw8P55ZdfGD58OJ06dSI5OZmqVavSuXNnh/le+vTpw7///kvTpk0pUqQIgwYN4tFHH818fvr06YwZM4ZnnnmGP//8k7Jly9KiRQunJv575plniIuLo1+/fgQEBNC/f3/uvPNOEhMTM8u88sorlC9fnrFjx3Lo0CFKly5No0aNeP755/P8e9q1a8dXX33FK6+8wrhx4yhVqhRt2rRx6lz27dvnEJerlCpVitWrV/POO++QlJRE1apVefPNN+nSpQsAjzzyCCtXrqRJkyacOnWKFStW0K5du2zryss1Gzx4MLfeeispKSm0adOGxYsXe/VcRCKuZjPO3IgXEb/Srl07GjRowDvvvOPpUEREXEq3wERERMTvKAESERERv6NbYCIiIuJ31AIkIiIifkcJkIiIiPgdJUAiIiLid5QAiYiIiN9RAiQiIiJ+RwmQiIiI+B0lQCIiIuJ3lACJiIiI3/l/oyNENTIb0C0AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "birthOrder = pandas.read_csv('https://raw.githubusercontent.com/jillxoreilly/StatsCourseBook/main/data/BirthOrderIndependence.csv')\n", "\n", "plt.scatter(data = birthOrder, x=\"FirstBorn\", y=\"SecondBorn\")\n", "plt.xlabel(\"independence: first born\")\n", "plt.ylabel(\"independence: second born\")\n", "\n", "# add the line x=y (ie a line from point(50,50) to (110,110)) for reference \n", "plt.plot([0,20],[0,20],'r--')" ] }, { "cell_type": "markdown", "id": "850462ee", "metadata": {}, "source": [ "Comments:\n", " \n", "" ] }, { "cell_type": "markdown", "id": "1e0d04e2", "metadata": {}, "source": [ "### Check assumption of normality\n", "\n", "In the case of paired data, the assumption we would need to meet to use a t-test is that \n", "the differences between conditions (for each participant) are normally distributed - let's add a column to our pandas data frame to contain the differences" ] }, { "cell_type": "code", "execution_count": 17, "id": "ce83eb84", "metadata": {}, "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", " \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", "
FirstBornSecondBornDiff
012102
118126
21315-2
317134
489-1
515123
616133
758-3
8810-2
91284
101385
1159-4
121486
13201010
1419145
1517116
1627-5
1757-2
1815132
1918126
\n", "
" ], "text/plain": [ " FirstBorn SecondBorn Diff\n", "0 12 10 2\n", "1 18 12 6\n", "2 13 15 -2\n", "3 17 13 4\n", "4 8 9 -1\n", "5 15 12 3\n", "6 16 13 3\n", "7 5 8 -3\n", "8 8 10 -2\n", "9 12 8 4\n", "10 13 8 5\n", "11 5 9 -4\n", "12 14 8 6\n", "13 20 10 10\n", "14 19 14 5\n", "15 17 11 6\n", "16 2 7 -5\n", "17 5 7 -2\n", "18 15 13 2\n", "19 18 12 6" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "birthOrder['Diff'] = birthOrder.FirstBorn - birthOrder.SecondBorn\n", "birthOrder" ] }, { "cell_type": "markdown", "id": "3953720c", "metadata": {}, "source": [ "Now let's plot the differences to get a sense of whether they are normally distributed." ] }, { "cell_type": "code", "execution_count": 18, "id": "8f5e728d", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABI50lEQVR4nO3deVhV1eL/8TeCgKmQU4CGitmgaQ4HU1QyS1HKKS2tbmaD3Wi4peS9OWY/TU0zM0s0p8wGs9SGW1RSmjnQ4FjfskFT8SpkaoJDgsL5/bGSInAAgXXO2Z/X85zH3O5z+Bwegg9rr72Wn9vtdiMiIiLiIBVsBxAREREpbypAIiIi4jgqQCIiIuI4KkAiIiLiOCpAIiIi4jgqQCIiIuI4KkAiIiLiOAG2A3iivLw89uzZQ9WqVfHz87MdR0RERM6C2+3m0KFD1K5dmwoVTj/GowJUhD179hAZGWk7hoiIiJTArl27uPDCC097jgpQEapWrQqYT2BISIjlNCIiInI2srKyiIyMzP85fjoqQEU4edkrJCREBUhERMTLnM30FU2CFhEREcdRARIRERHHUQESERERx1EBEhEREcdRARIRERHHUQESERERx1EBEhEREcdRARIRERHHUQESERERx1EBEhEREcdRARIRERHHUQESERERx1EBEhEREcdRARIRERHHUQESERERx1EBEhEREcdRARIRERHHUQESERERx1EBEhEREcdRARIRERHHUQESERERx1EBEhEREcdRARIRERHHUQESERERx1EBEhEREcdRARIRERHHUQESERERx7FegJKSkoiKiiI4OBiXy8WqVatOeW56ejq33norl156KRUqVGDQoEFFnrdkyRIaN25MUFAQjRs35q233iqj9CIiIuKNrBagRYsWMWjQIEaMGMHGjRuJjY0lPj6etLS0Is/Pzs6mVq1ajBgxgmbNmhV5TmpqKv369aN///5s3ryZ/v3707dvX7744ouyfCsiIiLiRfzcbrfb1gdv3bo1LVu2ZMaMGfnHGjVqRK9evZgwYcJpn3v11VfTvHlzpk6dWuB4v379yMrK4oMPPsg/1rVrV6pVq8bChQuLfK3s7Gyys7Pz/56VlUVkZCSZmZmEhISU4J2JiIhIecvKyiI0NPSsfn5bGwHKyclh/fr1xMXFFTgeFxfH2rVrS/y6qamphV6zS5cup33NCRMmEBoamv+IjIws8ccXERERz2etAO3bt4/c3FzCwsIKHA8LCyMjI6PEr5uRkVHs1xw2bBiZmZn5j127dpX444uIiIjnC7AdwM/Pr8Df3W53oWNl/ZpBQUEEBQWd08cUERER72FtBKhmzZr4+/sXGpnZu3dvoRGc4ggPDy/11xQRERHfYq0ABQYG4nK5SElJKXA8JSWFtm3blvh1Y2JiCr3msmXLzuk1RURExLdYvQSWmJhI//79iY6OJiYmhlmzZpGWlkZCQgJg5ubs3r2bBQsW5D9n06ZNABw+fJhff/2VTZs2ERgYSOPGjQF4+OGHueqqq5g4cSI9e/bknXfe4eOPP2b16tXl/v5ERETEM1ktQP369WP//v2MGTOG9PR0mjRpQnJyMvXq1QPMwod/XxOoRYsW+f+9fv16XnvtNerVq8eOHTsAaNu2La+//jojR45k1KhRXHTRRSxatIjWrVuX2/sSERERz2Z1HSBPVZx1BERERMQzeMU6QCIiIiK2qACJiIiI46gAiYiIiOOoAImIiIjjqACJiIiI46gAiYiIiOOoAImIiIjjqACJiIiI46gAiYiIiOOoAImIiIjjqACJiIiI46gAiYiIiOOoAImIiIjjqACJiIiI46gAiYiIiOOoAImIiIjjqACJiIiI46gAiYiIiOOoAImIiIjjqACJiIiI46gAiYiIiOOoAImIiIjjqACJiIiI46gAiYiIiOOoAImIiIjjqACJiIiI46gAiYiIiOOoAImIiIjjqACJiIiI46gAiYiIiOOoAImIiIjjqACJiIiI4wTYDiAi4uuOHYPduyE3FypVggsvBD8/26lEnE0jQCIiZSAnB+bMgfh4qF4dGjaESy+FunWhXj24+27YtMl2ShHn0giQiEgpcrvhjTdg2DDYsQNatoQBA0wB8veHw4dh82b46CN48UX4xz9gwgQzKiQi5cfP7Xa7bYfwNFlZWYSGhpKZmUlISIjtOCLiJY4cgXvvhVdfhbZtYeBAiIoq+tzcXHj/fViwwFwOe+staN++fPOK+Jri/PzWJTARkVKwYwe0aQNLlsCoUTBu3KnLD5jRoB49YO5cM/pzzTUwf355pRURFSARkXO0Ywd06AC//QZJSabMnK3QUJg0CeLi4M474fXXyyymiPyF5gCJiJyDtDTo2NFc0nrmGahVq/ivUbEiJCbCiRNw++1QowZ07lz6WUXkTxoBEhEpoUOHoGtXyM6Gp58uWfk5qUIFGDLETJru3Ru+/770copIYSpAIiIlkJcHt90Gu3bBk09CWNi5v2ZAAIwebUaAbr7ZFCsRKRsqQCIiJfD44/Df/8Lw4WZtn9JSqRKMHAnffQdDh5be64pIQSpAIiLFtHw5PPEE3HUXxMSU/us3bAj//CdMnQopKaX/+iKiAiQiUiy//WYmKjdvDrfeWnYfp08faNEC7rvPbKUhIqVLBUhEpBgeeAAyM83lqQpl+B3Uzw8efhh27oTJk8vu44g4lQqQiMhZWroUFi6EQYPgggvK/uPVqwc33mgWVdy+vew/noiTqACJiJyFrCz417/MFhfFWejwXN1+O4SEwL//XX4fU8QJVIBERM7CqFFm/s9DD5nLU+WlUiWzQvSSJbBuXfl9XBFfpwIkInIG69fD88/DHXeUzno/xdW5M9SvDyNGlP/HFvFVKkAiIqfhdptRn6goc2eWDf7+pnwtWwarVtnJIOJrVIBERE7jzTdh7VpzO7q/v70csbFwySVm4UUROXcqQCIip3DsmJl83K4duFx2s1SoAAMGwOrVsGaN3SwivkAFSETkFJ55BvbsgXvvtZ3EaNPGzAWaNMl2EhHvpwIkIlKEAwfMJqc9ekBkpO00RoUK0LcvvPuudosXOVcqQCIiRZg4EY4fNzu+e5Jrr4VatbQ6tMi5UgESEfmbPXtg2jSzCnO1arbTFBQYCL17w8svQ3q67TQi3st6AUpKSiIqKorg4GBcLherznCP58qVK3G5XAQHB9OgQQNmzpxZ6JypU6dy6aWXUqlSJSIjIxk8eDDHtJugiJylJ54wRaNvX9tJitatm7kcNmeO7SQi3stqAVq0aBGDBg1ixIgRbNy4kdjYWOLj40lLSyvy/O3bt3PdddcRGxvLxo0bGT58OA899BBLlizJP+fVV19l6NChjB49mi1btjB37lwWLVrEsGHDyuttiYgX27EDZs+Gm2+GKlVspylalSrQqRPMnAknTthOI+Kd/Nxut9vWB2/dujUtW7ZkxowZ+ccaNWpEr169mDBhQqHzH330Ud599122bNmSfywhIYHNmzeTmpoKwIMPPsiWLVv45JNP8s955JFH+PLLL884unRSVlYWoaGhZGZmEhISUtK3JyJe6N57zdo/r75qtqHwVFu3wj33mC0yeve2nUbEMxTn57e1EaCcnBzWr19PXFxcgeNxcXGsXbu2yOekpqYWOr9Lly6sW7eO48ePA9C+fXvWr1/Pl19+CcDPP/9McnIy119//SmzZGdnk5WVVeAhIs6TlgYvvgg33eTZ5QegYUNo0gSmT7edRMQ7WStA+/btIzc3l7C/bawTFhZGRkZGkc/JyMgo8vwTJ06wb98+AG6++WbGjh1L+/btqVixIhdddBEdO3Zk6NChp8wyYcIEQkND8x+RnnLPq4iUqwkToHJl6NXLdpKz06MHLF+uW+JFSsL6JGi/v22r7Ha7Cx070/l/Pf7pp58ybtw4kpKS2LBhA0uXLuW9995j7Nixp3zNYcOGkZmZmf/YtWtXSd+OiHipXbtg7lwz8dnTR39O6tABzj/fzFkSkeIJsPWBa9asib+/f6HRnr179xYa5TkpPDy8yPMDAgKoUaMGAKNGjaJ///4MHDgQgKZNm3LkyBH++c9/MmLECCpUKNz5goKCCAoKKo23JSJeavJkU3x69rSd5OwFBpp1gV5+2SzaWLGi7UQi3sPaCFBgYCAul4uUlJQCx1NSUmjbtm2Rz4mJiSl0/rJly4iOjqbiH//nHz16tFDJ8ff3x+12Y3G+t4h4sF9/NaMovXvDeefZTlM8Xbua/MnJtpOIeBerl8ASExOZM2cO8+bNY8uWLQwePJi0tDQSEhIAc2nq9ttvzz8/ISGBnTt3kpiYyJYtW5g3bx5z585lyJAh+ed0796dGTNm8Prrr7N9+3ZSUlIYNWoUPXr0wN/mVs4i4rGefdb8ecMNdnOURMOGZpf4F1+0nUTEu1i7BAbQr18/9u/fz5gxY0hPT6dJkyYkJydTr149ANLT0wusCRQVFUVycjKDBw9m+vTp1K5dm2nTptGnT5/8c0aOHImfnx8jR45k9+7d1KpVi+7duzNu3Lhyf38i4vkyM+H556F7d/DWVS+6doWkJNi7Fy64wHYaEe9gdR0gT6V1gEScY+JEGDUKXnsNata0naZksrLMth1PPgmJibbTiNjjFesAiYjYlp0NzzwDcXHeW37AjFy1awfz59tOIuI9VIBExLFeecVcNvLUPb+K49pr4Ztv4NtvbScR8Q4qQCLiSHl5MGmSGTmpW9d2mnN35ZVQtSosXGg7iYh3UAESEUf673/hxx/Npqe+IDAQYmPNXCbN7BQ5MxUgEXGkSZOgaVO4/HLbSUrPtdfC9u3wx1aIInIaKkAi4jipqbB2rW/M/fmrZs3MZO7XXrOdRMTzqQCJiOM8/TRERsIpFp33Wv7+cPXV8PrrkJtrO42IZ1MBEhFH2bYN3nrLrJtTxNaAXq9jR3Nn26pVtpOIeDYf/N9fROTUpk416+Z06WI7Sdlo1AjCwmDxYttJRDybCpCIOMaBAzBvHvToAUFBttOUDT8/czfY4sXmVn8RKZoKkIg4xgsvwIkT0KuX7SRlq0MH+OUXM9FbRIqmAiQijpCTA9OmQefOUK2a7TRlq3FjqFUL3nzTdhIRz6UCJCKO8PrrkJEBN91kO0nZq1AB2rfXZTCR01EBEhGf53bD5MnQpg3Uq2c7Tfno0AH27IEvvrCdRMQzqQCJiM9bvtxsFHrjjbaTlJ8mTcylvrfftp1ExDOpAImIz5s8GRo2hJYtbScpP/7+EBMDS5dqbzCRoqgAiYhP++47+PBDM/fHz892mvLVvj1s3Qrff287iYjnUQESEZ82ZYq5I6pjR9tJyl/LllCpki6DiRRFBUhEfNYvv8DLL5t1fypWtJ2m/AUFQatWZusPESlIBUhEfNb06WYuTPfutpPY0749fPWVuSNMRP6kAiQiPunoUVOAunaFqlVtp7GnTRtTAt95x3YSEc+iAiQiPmn+fDh40BkLH55O1apwxRXw7ru2k4h4FhUgEfE5ubnw9NNmMcCICNtp7IuJgRUr4MgR20lEPIcKkIj4nHfegZ9/hr59bSfxDDExkJ0Nn3xiO4mI51ABEhGf4nbDpEnQvDlcdpntNJ7hwguhbl147z3bSUQ8hwqQiPiU1avN/lca/SmodWv473+1OarISSpAIuJTnnwSGjQwdz/Jn9q2hYwM2LjRdhIRz6ACJCI+4//+D5KTzeiP07a9OJMmTaBKFV0GEzlJBUhEfMakSRAWBtdeazuJ5wkIMKtC63Z4EUMFSER8QloaLFwIffqYH/ZSWJs2sGGD2SJExOlUgETEJ0yaBJUrQ7dutpN4rlatzJ8ffWQ3h4gnUAESEa+XkQFz5kDv3mb3cylatWrQqJGZJyXidCpAIuL1nnnGXPa64QbbSTxfq1bw4Ydw4oTtJCJ2qQCJiFc7cACSkqBnT2dvenq2WreGzEyzVpKIk6kAiYhXe/ZZM5px4422k3iHSy+F88/XZTARFSAR8Vq//QZTp0L37mZ+i5yZvz9ER6sAiagAiYjXmjoVcnLg5pttJ/EurVvDpk2Qnm47iYg9KkAi4pX+OvpTvbrtNN4lOtqslL1sme0kIvaoAImIV3rmGY3+lNT555u5QB9+aDuJiD0qQCLidX791RSgnj01+lNSrVqZEaDcXNtJROxQARIRrzN+vPnzllvs5vBmrVqZJQQ2bLCdRMQOFSAR8SppaWbdn759ITTUdhrv1bix2R1el8HEqVSARMSrPP64+cGtdX/Ojb8/tGypAiTOpQIkIl7j66/hpZfgH//Qnl+loVUrsyL0wYO2k4iUPxUgEfEKbjckJsKFF0KPHrbT+IZWrcwk6E8+sZ1EpPypAImIV3j/ffOD+t57zcancu7CwqBePa0HJM6kAiQiHu/4cXjkETNnJSbGdhrf4nKZeUBut+0kIuVLBUhEPN60abB1K9x3n1nBWEpPdLS5s27rVttJRMqXCpCIeLT//Q9Gj4ZevaBhQ9tpfE/z5uaSoi6DidOoAImIRxs8GIKD4c47bSfxTZUqQZMmKkDiPCpAIuKxPvgAFi+G++83a/9I2XC5YPlyM9dKxClUgETEI2Vmwj33mFu1O3a0nca3tWoFhw/D55/bTiJSflSARMQjDR5sFuh75BFNfC5rDRuabUU++sh2EpHyowIkIh4nORlefNHc9RUWZjuN7zu5LYbmAYmTlKgAbd++vbRziIgAkJFhJjxfeSVcd53tNM7hcsH69fDbb7aTiJSPEhWghg0b0rFjR1555RWOHTtW2plExKFyc80+XydOwKOP6tJXeXK5IC/PTIYWcYISFaDNmzfTokULHnnkEcLDw7n33nv58ssvSxQgKSmJqKgogoODcblcrFq16rTnr1y5EpfLRXBwMA0aNGDmzJmFzjl48CAPPPAAERERBAcH06hRI5KTk0uUT0TKz4QJsGIFDB8O1avbTuMs4eEQGQkpKbaTiJSPEhWgJk2aMGXKFHbv3s2LL75IRkYG7du35/LLL2fKlCn8+uuvZ/U6ixYtYtCgQYwYMYKNGzcSGxtLfHw8aWlpRZ6/fft2rrvuOmJjY9m4cSPDhw/noYceYsmSJfnn5OTk0LlzZ3bs2MHixYv54YcfmD17NnXq1CnJWxWRcvLee/DYY3DbbWY0Qsqf5gGJk/i53ee+A0x2djZJSUkMGzaMnJwcKlasSL9+/Zg4cSIRERGnfF7r1q1p2bIlM2bMyD/WqFEjevXqxYQJEwqd/+ijj/Luu++yZcuW/GMJCQls3ryZ1NRUAGbOnMlTTz3F999/T8WKFUv0frKysggNDSUzM5OQkJASvYaInL2vv4Z27cyqxP/v/0EF3Z5hxerVMGoUbNsGDRrYTiNSfMX5+X1O32bWrVvH/fffT0REBFOmTGHIkCFs27aN5cuXs3v3bnr27HnK5+bk5LB+/Xri4uIKHI+Li2Pt2rVFPic1NbXQ+V26dGHdunUc/2MFr3fffZeYmBgeeOABwsLCaNKkCePHjyc3N/eUWbKzs8nKyirwEJHysWcPdOsGtWubS18qP/Y0b27uCPv4Y9tJRMpeib7VTJkyhaZNm9K2bVv27NnDggUL2LlzJ0888QRRUVG0a9eOF154gQ0bNpzyNfbt20dubi5hf7vHNSwsjIyMjCKfk5GRUeT5J06cYN++fQD8/PPPLF68mNzcXJKTkxk5ciRPP/0048aNO2WWCRMmEBoamv+IjIw820+FiJyDX34xixzm5MATT5htGcSeKlWgUSPNAxJnCCjJk2bMmMFdd93FnXfeSXh4eJHn1K1bl7lz557xtfz+dpuH2+0udOxM5//1eF5eHhdccAGzZs3C398fl8vFnj17eOqpp3jssceKfM1hw4aRmJiY//esrCyVIJEy9uuvcM015rbrZ56BWrVsJxIw84DefdfckefvbzuNSNkpUQFKSUmhbt26VPjbWLXb7WbXrl3UrVuXwMBABgwYcMrXqFmzJv7+/oVGe/bu3VtolOek8PDwIs8PCAigRo0aAERERFCxYkX8//J/bqNGjcjIyCAnJ4fAwMBCrxsUFERQUNDp37SIlJqffoL4+D/Lj+5R8BwuFyxYABs3QnS07TQiZadEl8Auuuii/EtOf3XgwAGioqLO6jUCAwNxuVyk/G2sNSUlhbZt2xb5nJiYmELnL1u2jOjo6PwJz+3atWPr1q3k5eXln/Pjjz8SERFRZPkRkfK1ejW0aWM23nzuOahb13Yi+avGjeG883QZTHxfiQrQqW4cO3z4MMHBwWf9OomJicyZM4d58+axZcsWBg8eTFpaGgkJCYC5NHX77bfnn5+QkMDOnTtJTExky5YtzJs3j7lz5zJkyJD8c+677z7279/Pww8/zI8//sj777/P+PHjeeCBB0ryVkWklJw4AWPHwtVXm/VmnnvOTHwWzxIQAM2aqQCJ7yvWJbCT82T8/Px47LHHOO+88/L/LTc3ly+++ILmzZuf9ev169eP/fv3M2bMGNLT02nSpAnJycnUq1cPgPT09AJrAkVFRZGcnMzgwYOZPn06tWvXZtq0afTp0yf/nMjISJYtW8bgwYO54oorqFOnDg8//DCPPvpocd6qiE/Jy4Njx8zKysHB5b/C8oYN8OCD8MUXZqXn2283P2jFM7VsCbNnw9GjZjRIxBcVax2gjh07AmY15piYmAKXlAIDA6lfvz5Dhgzh4osvLv2k5UjrAIm3crvh//7PbGewdi38+CP8/DP8dWWHChUgJAQuvNBcfrr0UnPnzxVXQNOmpfsDb8sWs7rzyy9D/fqQmGg+hni27dvhrrvM7vB/W3lExKMV5+d3sX4HW7FiBQB33nknzz77rMqBiIfYvRvmz4e5c80Pr4oVTamJjDSbilatCifn+f/+Oxw6BPv2mdvQ33zTrMWTl2dGhi65xEx+bdnSrAvTrBn8cY/BWWf56CNTej791GxpMXgwXH+97iryFvXrQ82aZj0gFSDxVaWyErSv0QiQeItffoHx42HmTDOy06EDdOpkRlmKc2NjTg7s2GHuztq61fy5bZu5bAbmh+Gll0JUlJm3U6OGWbPH39+MLv32mxlt+vZb81w/PzOi1L07xMaC7j/wPuPHm6UKNm60nUTk7BXn5/dZF6DevXszf/58QkJC6N2792nPXbp06dmn9UAqQOLp8vJM6fnPf0zZ6NcPbrgBKlcuvY+Rm2tGc7Ztg7Q02LXL/EDcv9+Unpwcc07lymYBvYgIc0mtcWMzghQaWnpZpPx99BE8+STs3as1msR7lMklsNDQ0PzFBkP1nU3Eml27oH9/WLkSevSAgQPNJa7S5u9vCo1uU3emli3Nn8uXm4It4mt0CawIGgEST7VmjRnpqVAB/v1v7ZouZeuuu+Daa80dYSLeoMw3Q/399985evRo/t937tzJ1KlTWbZsWUleTkTOwssvm32zIiJgxgyVHyl7LVrAsmXm7kIRX1OiAtSzZ08WLFgAwMGDB7nyyit5+umn6dmzJzNmzCjVgCICL7xg1s7p1AkmT4Zq1WwnEidwucz8r59/tp1EpPSVqABt2LCB2NhYABYvXkx4eDg7d+5kwYIFTJs2rVQDijjd889DQoK59PXvf5tb3EXKQ7NmZi6YVoUWX1SiAnT06FGq/jHrctmyZfTu3ZsKFSrQpk0bdu7cWaoBRZzstdfgX/+Cm24yf5b3Cs7ibJUrm/WkPv7YdhKR0leiAtSwYUPefvttdu3axUcffUTcHytl7d27V5OGRUrJ8uVwxx3QpQvcd5/Kj9jhcsEnn5glD0R8SYkK0GOPPcaQIUOoX78+rVu3JiYmBjCjQS1atCjVgCJO9P335pJX8+YwZIjKj9jjcsHBg1oQUXxPibYjvPHGG2nfvj3p6ek0a9Ys//i1117LDTfcUGrhRJzo8GFTfqpXh8cf16ahYlejRmZ/uJQUs8CliK8o8bfW8PBwwsPDCxy78sorzzmQiJO53XD33ebOmxkztBO32BcQYCZDp6TAsGG204iUnhIVoCNHjvDkk0/yySefsHfvXvLy8gr8+8+6Z1KkRJKS4I03zMiPVmAWT9GypVkM8ehRlXLxHSUqQAMHDmTlypX079+fiIiI/C0yRKTkvv/e3Obes6fZ1FTEU7hcMH06rF6t3eHFd5SoAH3wwQe8//77tGvXrrTziDjS8eNw221m08mEBNtpRAqqX998bX78sQqQ+I4S3QVWrVo1qlevXtpZRBxr3DjYtMnMsQgOtp1GpCA/vz+3xRDxFSUqQGPHjuWxxx4rsB+YiJTMt9/C+PFw661w2WW204gUzeWCzZvh119tJxEpHSW6BPb000+zbds2wsLCqF+/PhX/tjb/hg0bSiWciK/Ly4N77jEbnN52m+00Iqd2cvPdTz6Bm2+2m0WkNJSoAPXq1auUY4g40+zZkJoKzzwDgYG204icWo0aEBVlbodXARJfUKICNHr06NLOIeI4v/4Kjz4K8fFmxWcRT+dymXlAbrdWJxfvV6I5QAAHDx5kzpw5DBs2jAMHDgDm0tfu3btLLZyILxs1ylwC++c/bScROTsuF/zvf/Djj7aTiJy7Eo0Aff3113Tq1InQ0FB27NjBPffcQ/Xq1XnrrbfYuXMnCxYsKO2cIj5l0yaYNQvuvx/OP992GpGz06yZWRk6JQUuvdR2GpFzU6IRoMTERO644w5++ukngv9yz258fDyfffZZqYUT8UVuNzz8sFnpWdPpxJtUqgRNmpgCJOLtSlSAvvrqK+69995Cx+vUqUNGRsY5hxLxZe+8A599ZkZ/tNGpeBuXC1asMIt3inizEhWg4OBgsrKyCh3/4YcfqFWr1jmHEvFVJ07A0KFmV23tHSzeyOWCQ4fgyy9tJxE5NyUqQD179mTMmDEc/+NXAD8/P9LS0hg6dCh9+vQp1YAivuSll+CHH2DgQNtJRErmkksgJESXwcT7lagATZ48mV9//ZULLriA33//nQ4dOtCwYUOqVq3KuHHjSjujiE/4/Xd47DHo2FETSMV7+fubZRu0LYZ4uxLNQAgJCWH16tWsWLGC9evXk5eXR8uWLenUqVNp5xPxGUlJsHcvTJpkO4nIuYmOhmefhcxMCA21nUakZIpdgPLy8pg/fz5Lly5lx44d+Pn5ERUVRXh4OG63Gz+tjiVSyJEj8OST0KUL1KljO43IuYmOhtxcMxladzKKtyrWJTC3202PHj0YOHAgu3fvpmnTplx++eXs3LmTO+64gxtuuKGscop4taQkOHhQ+32Jb4iIgMhIXQYT71asEaD58+fz2Wef8cknn9CxY8cC/7Z8+XJ69erFggULuP3220s1pIg3O3IEJk6Erl0hPNx2GpHS4XLBRx/ZTiFScsUaAVq4cCHDhw8vVH4ArrnmGoYOHcqrr75aauFEfEFSkpkr8Y9/2E4iUnqio+Hnn2HbNttJREqmWAXo66+/pmvXrqf89/j4eDZv3nzOoUR8xdGjZtKzRn/E1zRvbhby1GUw8VbFKkAHDhwgLCzslP8eFhbGb7/9ds6hRHzFvHlw4ADccovtJCKlq3JluPxyFSDxXsUqQLm5uQScZu1+f39/Tpw4cc6hRHzB8eNm9Oeaa6B2bdtpREqfywXLl2tbDPFOxZoE7Xa7ueOOOwgKCiry37Ozs0sllIgveO012LULHn/cdhKRshEdbUY5v/gC2re3nUakeIpVgAYMGHDGc3QHmAjk5cGECdC2LTRoYDuNSNm45BKzEOJHH6kAifcpVgF68cUXyyqHiE957z2z59dzz9lOIlJ2/P2hZUv48EMYO9Z2GpHiKdFeYCJyepMnQ5Mm5iHiy1q1gvXrYd8+20lEikcFSKSUffklrFoFN91kO4lI2WvVCtxu7Q4v3kcFSKSUPf202e+rXTvbSUTKXs2acNFF5jKYiDdRARIpRTt2wOLFcOONZn6EiBNER5uJ0G637SQiZ08FSKQUPfccVKlidn0XcYpWreCXX+Drr20nETl7KkAipeTwYZgzB66/HipVsp1GpPw0bQrBwfDBB7aTiJw9FSCRUvLSS2bn9169bCcRKV+BgdCihQqQeBcVIJFSkJcHzz4LsbFwwQW204iUvyuvhLVrISvLdhKRs6MCJFIKPvoIfvoJeve2nUTEjtat4cQJ+Phj20lEzo4KkEgpmDYNLr1UCx+Kc0VEQL16ugwm3kMFSOQcbd1q1kDp1Qv8/GynEbGnVStITtbt8OIdVIBEzlFSktkQsmNH20lE7GrdGvbsgW++sZ1E5MxUgETOwZEjMG8exMdDUJDtNCJ2XXGFbocX76ECJHIOXnvN3PXSs6ftJCL2BQaa3eGTk20nETkzFSCREnK74fnnISYGwsNtpxHxDK1bw5o1cPCg7SQip6cCJFJCn39ulv7X6I/In9q0gdxcszSEiCdTARIpoZkzoXZtsxGkiBgXXAANG8L779tOInJ6KkAiJbB/PyxaBN26QQX9XyRSwJVXmnlAubm2k4icmvVv3UlJSURFRREcHIzL5WLVqlWnPX/lypW4XC6Cg4Np0KABM2fOPOW5r7/+On5+fvTS5kxSyubPN9/c4+NtJxHxPG3amF8SvvzSdhKRU7NagBYtWsSgQYMYMWIEGzduJDY2lvj4eNLS0oo8f/v27Vx33XXExsayceNGhg8fzkMPPcSSJUsKnbtz506GDBlCbGxsWb8NcZi8PHP5q0MHOP9822lEPE/jxmZtLF0GE0/m53bbW7OzdevWtGzZkhkzZuQfa9SoEb169WLChAmFzn/00Ud599132bJlS/6xhIQENm/eTGpqav6x3NxcOnTowJ133smqVas4ePAgb7/99ilzZGdnk52dnf/3rKwsIiMjyczMJCQk5Bzfpfia5cvh2mth6lRo1sx2GhHPNG4c/PorbNpkO4k4SVZWFqGhoWf189vaCFBOTg7r168nLi6uwPG4uDjWrl1b5HNSU1MLnd+lSxfWrVvH8ePH84+NGTOGWrVqcffdd59VlgkTJhAaGpr/iIyMLOa7ESd54QWoX98s+iYiRYuJgc2bYdcu20lEimatAO3bt4/c3FzCwsIKHA8LCyMjI6PI52RkZBR5/okTJ9i3bx8Aa9asYe7cucyePfusswwbNozMzMz8xy79HyunsHcvvPUWXH+99v0SOZ0rr4SAAHjvPdtJRIpmfRK0399+irjd7kLHznT+yeOHDh3itttuY/bs2dSsWfOsMwQFBRESElLgIVKUl14yxedvA5Ei8jdVqphR0nfesZ1EpGgBtj5wzZo18ff3LzTas3fv3kKjPCeFh4cXeX5AQAA1atTg22+/ZceOHXTv3j3/3/Py8gAICAjghx9+4KKLLirldyJOkZdnLn916ADqyCJnFhMDs2fDoUNQtartNCIFWRsBCgwMxOVykZKSUuB4SkoKbdu2LfI5MTExhc5ftmwZ0dHRVKxYkcsuu4xvvvmGTZs25T969OhBx44d2bRpk+b2yDn59FPYts2s/SMiZ9a2LeTkwLJltpOIFGZtBAggMTGR/v37Ex0dTUxMDLNmzSItLY2EhATAzM3ZvXs3CxYsAMwdX88//zyJiYncc889pKamMnfuXBYuXAhAcHAwTZo0KfAxzv/jPuW/HxcprlmzoF49aNrUdhIR71C7NjRoAO++C3362E4jUpDVAtSvXz/279/PmDFjSE9Pp0mTJiQnJ1OvXj0A0tPTC6wJFBUVRXJyMoMHD2b69OnUrl2badOm0Uf/Z0kZ27fPTH6++25NfhYpjjZtzETo3Fzw97edRuRPVtcB8lTFWUdAnGHKFBg6FN580yzwJiJn57vv4IEH4LPPQOvSSlnzinWARLyF220uf7Vvr/IjUlyXXQY1asBp1qIVsUIFSOQM1qyBH34wa/+ISPFUqGAmQy9dan6ZEPEUKkAiZzBnDtSpAy1a2E4i4p3at4cdO+Drr20nEfmTCpDIaWRmwhtvQNeu5jdZESm+Fi3Mwoi6DCaeRN/SRU5j4ULIzjYFSERKpmJFczfY0qW2k4j8SQVI5DRmzzbfuIuxs4qIFKFdO3MJbPt220lEDBUgkVPYtAk2bIDrrrOdRMT7tW4NgYFmPS0RT6ACJHIKc+eakZ82bWwnEfF+lSpBq1aweLHtJCKGCpBIEX7/HV5+2ez6rtVrRUrHVVdBairs3m07iYgKkEiRli41d4Dp8pdI6WnbFgICNBlaPIMKkEgR5swxt+7WqWM7iYjvqFIFoqPNljIitqkAifzNtm3w6acQH287iYjvueoqWL0aMjJsJxGnUwES+Zt588xvqlddZTuJiO9p29YsKqq7wcQ2FSCRvzhxwhSga6+FoCDbaUR8T2ioubz8xhu2k4jTqQCJ/MWHH5qheU1+Fik7V18NK1dCerrtJOJkKkAifzFnDlxyiXmISNmIjTXLS2hNILFJBUjkDxkZ8N572vdLpKyFhJhFERcutJ1EnEwFSOQPL71k1ijp1Ml2EhHf17GjWRRx507bScSpVIBEALfbXP6KjYWqVW2nEfF9bduaGw00GVpsUQESAVatgq1b4frrbScRcYbKlc0GqboMJraoAIlgNj6tUweaNbOdRMQ5rrkGNm6EH36wnUScSAVIHO/gQbM0/3XXgZ+f7TQizhETYxYdfeUV20nEiVSAxPEWLoScHOjSxXYSEWcJDDQrrr/8spmHJ1KeVIDE8WbPNr+J1qhhO4mI83TubO4EW7vWdhJxGhUgcbQNG8wcBK38LGLHFVfABRfoMpiUPxUgcbQ5c6BWLbjySttJRJypQgWz997rr5tL0SLlRQVIHOvIEXj1VTP3x9/fdhoR5+rc2dyM8P77tpOIk6gAiWO9+SZkZenyl4htUVFm/70XX7SdRJxEBUgca9YsiI6GiAjbSUSka1dIToZffrGdRJxCBUgc6dtvzT5EWvlZxDNce62ZD6TJ0FJeVIDEkebMgWrVoF0720lEBMwO8e3awbx5WhNIyocKkDjOsWNm5/e4OKhY0XYaETmpa1f47jtYt852EnECFSBxnCVL4LffdPlLxNNER5tlKebNs51EnEAFSBznhRegRQuIjLSdRET+yt/fLEvx6qtmmQqRsqQCJI6yZQusWgXdutlOIiJFuf56OHwYFi2ynUR8nQqQOMrs2XD++dC+ve0kIlKU8HBo1QpmzrSdRHydCpA4xrFjMH++GWIPDLSdRkROpVs3+Oor2LTJdhLxZSpA4hiLF2vys4g3iImBmjXNiK1IWVEBEseYMQNcLk1+FvF0AQEQHw8vvwyHDtlOI75KBUgc4ZtvYO1a6N7ddhIRORvdusHRo6YEiZQFFSBxhJkzoUYNrfws4i0uuMD8//rcc1oZWsqGCpD4vMOHzW+R111nhtZFxDvccAN8/z2sWGE7ifgiFSDxeScXVdPkZxHv0qwZREXB9Om2k4gvUgESn+Z2m2+eMTEQFmY7jYgUh58f9OoFb78NO3faTiO+RgVIfNqaNWYCdM+etpOISEl07gxVqsC0abaTiK9RARKfNn26ue3d5bKdRERKolIlc0fY7NmQmWk7jfgSFSDxWRkZZuf3Hj2ggr7SRbzWDTfA77/D3Lm2k4gv0Y8F8Vlz5vy5u7SIeK+aNeHaa2HqVDhxwnYa8RUqQOKTjh+HpCTzTbNqVdtpRORc3Xgj7NoFb75pO4n4ChUg8UlvvQXp6dC7t+0kIlIaGjY0u8Q/+aQWRpTSoQIkPunZZ6FFC2jQwHYSESktt94KX38NH3xgO4n4AhUg8Tnr15t9v264wXYSESlNzZrB5ZfD+PG2k4gvUAESn/PccxAeDm3b2k4iIqXJz8+MAq1ZA6tW2U4j3k4FSHzKL7/AwoXm1nd/f9tpRKS0tWkDF10EY8faTiLeTgVIfMqMGWbNn27dbCcRkbJQoQLcdhukpEBqqu004s1UgMRnHDtmVn7u2lW3vov4squuMjc4PPaY7STizawXoKSkJKKioggODsblcrHqDBd2V65cicvlIjg4mAYNGjBz5swC/z579mxiY2OpVq0a1apVo1OnTnz55Zdl+RbEQ7z6KuzfD3362E4iImWpQgW4/Xb4+GMzH0ikJKwWoEWLFjFo0CBGjBjBxo0biY2NJT4+nrS0tCLP3759O9dddx2xsbFs3LiR4cOH89BDD7FkyZL8cz799FNuueUWVqxYQWpqKnXr1iUuLo7du3eX19sSC9xumDLFTHy+8ELbaUSkrMXGmrWBNAokJeXndttbUqp169a0bNmSGTNm5B9r1KgRvXr1YsKECYXOf/TRR3n33XfZsmVL/rGEhAQ2b95M6ikuBufm5lKtWjWef/55br/99iLPyc7OJjs7O//vWVlZREZGkpmZSUhISEnfnpSjZcvMlhfPPAPNm9tOIyLlYfVqGDXKjARde63tNOIJsrKyCA0NPauf39ZGgHJycli/fj1xcXEFjsfFxbF27doin5Oamlro/C5durBu3TqOHz9e5HOOHj3K8ePHqV69+imzTJgwgdDQ0PxHZGRkMd+N2DZpElx6qVknREScoV07sy7Qo49qdWgpPmsFaN++feTm5hIWFlbgeFhYGBkZGUU+JyMjo8jzT5w4wb59+4p8ztChQ6lTpw6dOnU6ZZZhw4aRmZmZ/9i1a1cx343YtGEDfPIJ9O1r1gkREWfw84OBA83ip4sX204j3ibAdgC/v/3EcrvdhY6d6fyijgNMmjSJhQsX8umnnxIcHHzK1wwKCiIoKKg4scWDTJ4MERHQoYPtJCJS3po3h9atYfhw6NULKla0nUi8hbURoJo1a+Lv719otGfv3r2FRnlOCg8PL/L8gIAAatSoUeD45MmTGT9+PMuWLeOKK64o3fDiMXbsgDfeMDtFa+FDEWcaOBC2bYNZs2wnEW9irQAFBgbicrlISUkpcDwlJYW2p9jDICYmptD5y5YtIzo6mop/qf1PPfUUY8eO5cMPPyQ6Orr0w4vHmDIFqlSB+HjbSUTEloYNzfeAUaPgt99spxFvYfU2+MTERObMmcO8efPYsmULgwcPJi0tjYSEBMDMzfnrnVsJCQns3LmTxMREtmzZwrx585g7dy5DhgzJP2fSpEmMHDmSefPmUb9+fTIyMsjIyODw4cPl/v6kbO3dC3PmmE1PK1WynUZEbLrrLsjO1hYZcvasFqB+/foxdepUxowZQ/Pmzfnss89ITk6mXr16AKSnpxdYEygqKork5GQ+/fRTmjdvztixY5k2bRp9/rLyXVJSEjk5Odx4441ERETkPyZPnlzu70/K1rPPmj+167uI1KgBt9xiNkP+8UfbacQbWF0HyFMVZx0BsSMzE+rWNcPefwwYiojDZWfDHXeYidHJybor1Im8Yh0gkXORlAS//w433WQ7iYh4iqAgeOAB+PBDeOcd22nE06kAidc5fBieftqM/vzt5j8Rcbh27SAmBh56CI4csZ1GPJkKkHidGTPMJbBbbrGdREQ8jZ+fGQX65Rd44gnbacSTqQCJVzlyxGx70bUrhIfbTiMinqhOHfjHP8wiqV9/bTuNeCoVIPEqM2bAwYPmm5uIyKnccgtERsLdd0Nuru004olUgMRrnBz96dJFoz8icnoVK8KQIWafsGnTbKcRT6QCJF7juefMKq+33WY7iYh4g8aNzTphI0bATz/ZTiOeRgVIvMLBgzBxInTrptEfETl7AwdC9eowYIAuhUlBKkDiFZ5+Go4d0+iPiBRPpUrw6KPw+edm70CRk1SAxOP9+itMnQq9emndHxEpvqZNoW9fGDkSvvnGdhrxFCpA4vHGjgW3G26+2XYSEfFWd91l7gq7+WaziryICpB4tK1bza3vt94KoaG204iItwoMNJOht20zd4eJqACJRxs2zExg7NPHdhIR8XZRUXDffWYvwbfftp1GbFMBEo/1+eeweDHceafZ5FBE5Fz16AGxsWbX+J9/tp1GbFIBEo/kdsPgwXDRRdC5s+00IuIr/PzgP/+BKlXMyPKxY7YTiS0qQOKRXnvNjAA98AD4+9tOIyK+pEoVePxx+O47+Ne/bKcRW1SAxOMcPgz//jdcdRW0aGE7jYj4ooYN4eGHYc4ceOEF22nEBhUg8TgTJ8L+/ZCQYDuJiPiy664z64s9+CCsXm07jZQ3FSDxKD/9ZDY87dsXIiJspxERX/fAA3D55dC7N+zYYTuNlCcVIPEYbjfcf79Z7fkf/7CdRkScICDAzAeqWBGuvx4yM20nkvKiAiQe44034OOPzaTE4GDbaUTEKc4/H8aPh7Q0uOkmOHHCdiIpDypA4hEyM2HQILM+R0yM7TQi4jT16pmRoBUr4J//NCPS4ttUgMQj/PvfkJVlJiOKiNjgcpnvRS++CI89ZjuNlDUVILFu+XKYPdv81nXBBbbT/Gn/fpg/3/zpbc9ft85s+rhuXfl/bNu2bjWjiVu3Ou/5y5dDt27mTycqja/buDjzveiJJ+D550stmnggFSCx6sgRuPtuaNYMune3naag/fvhpZfOrcDYev6WLfDLL+bP8v7Ytu3YAZs3l/yOHm9+/ubN5v+pzZtL9rG9XWl93d58s5kL9K9/mdcT3xRgO4A426OPQno6jB0LFVTHRcQD+PmZTVOPHoW77oLKleHGG22nktKmAiTWfPghTJ8ODz0EF15oO42IyJ/8/Mx+hMeOwS23mNvle/WynUpKk37nFiv27ze7vLdqBT172k4jIlKYvz8MGwbt25tLYu+8YzuRlCYVICl3breZZHj0qNmVWZe+RMRT+fvDyJGmBN14I7z5pu1EUlr0o0fK3fTpsHQpDBkCNWvaTiMicnonS9DVV5sJ0gsW2E4kpUFzgKRcbdgAjzxi9t2JjbWdRkTk7Pj7w9ChEBQEAwaYxVv/9S/bqeRcqABJufntN3MdPSoK7r3XdhoRkeLx9ze/wFWubG7e+OUXcwern5/tZFISKkBSLnJzzdDx/v0wYwYEBtpOJCJSfCdvka9eHcaNg927YdYss5mqeBcVICkXI0eajU4nTYKICNtpRETOTb9+pgRNmmQ2UV2yxGyqKt5Dk6ClzL3yCjz5pLnzy+WynUZEpHR07gxPPQVffQVt28K2bbYTSXGoAEmZ+vRTs5JqfDz07Ws7jYhI6Wre3OwZduiQWdfMqfuweSMVICkz331nVk5t1gwSEzVRUER8U926kJQEDRqYzVSfecasdyaeTQVIysSOHWZ4uGZNGD3aLCMvIuKrqlaFiROhTx/zC98tt8Dhw7ZTyemoAEmpS0+Ha681Iz4TJ0KVKrYTiYiUPX9/c4fY44/Df/9r5jx+/bXtVHIqKkBSqjIyTPk5dAgmT4YaNWwnEhEpXx06wMyZkJcHV15plv7QJTHPowIkpWbPHvM//r59pvyEh9tOJCJiR2Sk2fanSxe4/37o3h327rWdSv5KBUhKxc8/w1VXmeXhn3nGTAoUEXGywEAYPNgsmLh2LVx+uVkvSDyDCpCcsw0bICYGsrNN+alTx3YiERHP0bYtzJkDjRqZHeVvuUWjQZ5ABUjOyfvvm8teNWrAtGla5VlEpCjVq8P/+38wYgR88AFcdhnMn6+5QTapAEmJuN3mDq/u3c06P08/DdWq2U4lIuK5/PygUydTfKKj4c47zS+Q33xjO5kzqQBJsWVmmo1Nhw6Ff/wDxoyBSpVspxIR8Q7nnw/Dh5ubRXbuhBYt4OGH4cAB28mcRQVIimXdOvM/6/vvm7Uu7r4bKuirSESk2FwuMzdo4ECYO9esJD1lChw7ZjuZM+hHl5yVnBxTeGJizJ0NL7xghm5FRKTkKlY0I+ovv2y+p/7nP3DxxaYQnThhO51vUwGSM/ryS7OY1xNPwK23wnPP6U4vEZHSVK2auWV+3jxTgAYOhEsuMUXo+HHb6XyTCpCc0r59kJAAbdrA0aNmUa877zS/sYiISOmrWxceewxmzzaLKQ4caC6NPfOMWWFfSo8KkBTy++/w5JNw0UXwyivw4INmKfdLL7WdTETEGRo2NLfNz5tnFlD8z3/gwgvNKNHWrbbT+QYVIMl35Ii5nb1+fRg1yuzp9cor0Lu32eRPRETKV1SUueP2tdegWzd48UVziaxTJ3jjDbMArZRMgO0AYt/u3eby1syZZog1Ls7M9dE8HxERz1CrFtxzD9x+O6xYAcnJ0K+fuaW+Xz+zJEm7drortzhUgBzq+HH46CNznfn99yE4GOLjoU8fbWIqIuKpgoKga1fz2LkTli2Dt982d+ZGRJitNnr2NHszar7m6akAOcjx4/DZZ7B4Mbz5Juzfb+4yeOABM+pTubLthCIicrbq1TOjQnffDd9+C59+CosWmTt1Q0LM9/UuXczlsvr1baf1PCpAPszthh074JNPICUFPvwQsrLMCE/nznDNNeZasoiIeK8KFaBpU/N48EH46SdITTUL1y5dCnl5pixdfbW5TBYTYzZmdfrcTutXC5OSkoiKiiI4OBiXy8WqVatOe/7KlStxuVwEBwfToEEDZs6cWeicJUuW0LhxY4KCgmjcuDFvvfVWWcX3KIcPw5o18OyzZrfhyEhz++S998L//Z+ZzDxrlplMd++9Kj8iIr7Gz8+M7A8YYEaC3n4bxo41e499/rlZ2qRpUzN36KqrYNAgM7F63TpzI4yTWB0BWrRoEYMGDSIpKYl27drxwgsvEB8fz3fffUfdunULnb99+3auu+467rnnHl555RXWrFnD/fffT61atejTpw8Aqamp9OvXj7Fjx3LDDTfw1ltv0bdvX1avXk3r1q3L+y2WuhMnYM8e2L4dtm2DH36A7783m+lt327OCQoy5aZ9e2jSBJo3hypVrMYWERELqlY1Pwvatzd/P3r0z58bP/xgRoimTftzV/q6dc3Pj0suMb9AR0WZY3XrmonYvjTJ2moBmjJlCnfffTcDBw4EYOrUqXz00UfMmDGDCRMmFDp/5syZ1K1bl6lTpwLQqFEj1q1bx+TJk/ML0NSpU+ncuTPDhg0DYNiwYaxcuZKpU6eycOHC8nljZykvz3wxHjpkLk1lZsLBg2ZDvP374ddfzSMjw5Se//0P0tMhN9c838/PXM6KjDTtvm9fs3ZE/foQoIubIiLyN+edZ/ZzbNHiz2O//24mVO/YAbt2mceyZebnzdGjf54XEGB+5kREQFgYXHCBKUU1akD16mY169BQ8wgJMeWrShXzMT3xcpu1H5M5OTmsX7+eoUOHFjgeFxfH2rVri3xOamoqcXFxBY516dKFuXPncvz4cSpWrEhqaiqDBw8udM7J0lSU7Oxssv+ymEJmZiYAWVlZxXlLZ+Wrr8yEtJLw84OaNc0X3ckvvsDAP/89Pd08znAVUc7S/v3mz3ffNf+De9Pzv/7a/LlhQ8n2EzrX7Db9/LP5c+VK843cSc///vs//3zxxeJ/bG/nzV+3niIgwIz6REWZUaFDh/78ZXz/frNsyv/+Vzofq29fcydyaTr5c9t9ckjrdNyW7N692w2416xZU+D4uHHj3JdcckmRz7n44ovd48aNK3BszZo1bsC9Z88et9vtdlesWNH96quvFjjn1VdfdQcGBp4yy+jRo92AHnrooYceeujhA49du3adsYdYv1Di5+dX4O9ut7vQsTOd//fjxX3NYcOGkZiYmP/3vLw8Dhw4QI0aNU77PG+QlZVFZGQku3btIiQkxHYcr6LPXcnpc1dy+tyVnD53JeNLnze3282hQ4eoXbv2Gc+1VoBq1qyJv78/GRkZBY7v3buXsLCwIp8THh5e5PkBAQHU+GO881TnnOo1AYKCgggKCipw7Pzzzz/bt+IVQkJCvP4L2xZ97kpOn7uS0+eu5PS5Kxlf+byFhoae1XnW5nMHBgbicrlISUkpcDwlJYW2bdsW+ZyYmJhC5y9btozo6Ggq/rHk5anOOdVrioiIiPNYvQSWmJhI//79iY6OJiYmhlmzZpGWlkZCQgJgLk3t3r2bBQsWAJCQkMDzzz9PYmIi99xzD6mpqcydO7fA3V0PP/wwV111FRMnTqRnz5688847fPzxx6xevdrKexQRERHPY7UA9evXj/379zNmzBjS09Np0qQJycnJ1KtXD4D09HTS0tLyz4+KiiI5OZnBgwczffp0ateuzbRp0/JvgQdo27Ytr7/+OiNHjmTUqFFcdNFFLFq0yCfWACqJoKAgRo8eXegSn5yZPnclp89dyelzV3L63JWMUz9vfm732dwrJiIiIuI7fGhNRxEREZGzowIkIiIijqMCJCIiIo6jAiQiIiKOowLkw8aNG0fbtm0577zzTrmwY1paGt27d6dy5crUrFmThx56iJycnPIN6gXq16+Pn59fgcff97ETIykpiaioKIKDg3G5XKzS5nRn9Pjjjxf6+goPD7cdyyN99tlndO/endq1a+Pn58fbb79d4N/dbjePP/44tWvXplKlSlx99dV8++23dsJ6mDN97u64445CX4dt2rSxE7YcqAD5sJycHG666Sbuu+++Iv89NzeX66+/niNHjrB69Wpef/11lixZwiOPPFLOSb3DyeUaTj5GjhxpO5LHWbRoEYMGDWLEiBFs3LiR2NhY4uPjCyxnIUW7/PLLC3x9ffPNN7YjeaQjR47QrFkznn/++SL/fdKkSUyZMoXnn3+er776ivDwcDp37syhQ4fKOannOdPnDqBr164Fvg6Tk5PLMWE5O+NuYeL1XnzxRXdoaGih48nJye4KFSq4d+/enX9s4cKF7qCgIHdmZmY5JvR89erVcz/zzDO2Y3i8K6+80p2QkFDg2GWXXeYeOnSopUTeYfTo0e5mzZrZjuF1APdbb72V//e8vDx3eHi4+8knn8w/duzYMXdoaKh75syZFhJ6rr9/7txut3vAgAHunj17Wsljg0aAHCw1NZUmTZoU2DSuS5cuZGdns379eovJPNPEiROpUaMGzZs3Z9y4cbpU+Dc5OTmsX7+euLi4Asfj4uJYu3atpVTe46effqJ27dpERUVx88038/PPP9uO5HW2b99ORkZGga/BoKAgOnTooK/Bs/Tpp59ywQUXcMkll3DPPfewd+9e25HKjPXd4MWejIyMQpvEVqtWjcDAwEIbyjrdww8/TMuWLalWrRpffvklw4YNY/v27cyZM8d2NI+xb98+cnNzC31NhYWF6evpDFq3bs2CBQu45JJL+OWXX3jiiSdo27Yt3377bf5Gz3JmJ7/Oivoa3Llzp41IXiU+Pp6bbrqJevXqsX37dkaNGsU111zD+vXrfXKVaI0AeZmiJkv+/bFu3bqzfj0/P79Cx9xud5HHfU1xPpeDBw+mQ4cOXHHFFQwcOJCZM2cyd+5c9u/fb/ldeJ6/f+045evpXMTHx9OnTx+aNm1Kp06deP/99wF46aWXLCfzTvoaLJl+/fpx/fXX06RJE7p3784HH3zAjz/+mP/16Gs0AuRlHnzwQW6++ebTnlO/fv2zeq3w8HC++OKLAsd+++03jh8/Xug3KF90Lp/Lk3dGbN26Vb+h/6FmzZr4+/sXGu3Zu3evI76eSlPlypVp2rQpP/30k+0oXuXknXMZGRlERETkH9fXYMlERERQr149n/06VAHyMjVr1qRmzZql8loxMTGMGzeO9PT0/G8Wy5YtIygoCJfLVSofw5Ody+dy48aNAAW+yTpdYGAgLpeLlJQUbrjhhvzjKSkp9OzZ02Iy75Odnc2WLVuIjY21HcWrREVFER4eTkpKCi1atADM3LSVK1cyceJEy+m8z/79+9m1a5fPfp9TAfJhaWlpHDhwgLS0NHJzc9m0aRMADRs2pEqVKsTFxdG4cWP69+/PU089xYEDBxgyZAj33HMPISEhdsN7kNTUVD7//HM6duxIaGgoX331FYMHD6ZHjx7UrVvXdjyPkpiYSP/+/YmOjiYmJoZZs2aRlpZGQkKC7WgebciQIXTv3p26deuyd+9ennjiCbKyshgwYIDtaB7n8OHDbN26Nf/v27dvZ9OmTVSvXp26desyaNAgxo8fz8UXX8zFF1/M+PHjOe+887j11lstpvYMp/vcVa9enccff5w+ffoQERHBjh07GD58ODVr1izwC41PsXwXmpShAQMGuIFCjxUrVuSfs3PnTvf111/vrlSpkrt69eruBx980H3s2DF7oT3Q+vXr3a1bt3aHhoa6g4OD3Zdeeql79OjR7iNHjtiO5pGmT5/urlevnjswMNDdsmVL98qVK21H8nj9+vVzR0REuCtWrOiuXbu2u3fv3u5vv/3WdiyPtGLFiiK/rw0YMMDtdptb4UePHu0ODw93BwUFua+66ir3N998Yze0hzjd5+7o0aPuuLg4d61atdwVK1Z0161b1z1gwAB3Wlqa7dhlxs/tdrst9C4RERERa3QXmIiIiDiOCpCIiIg4jgqQiIiIOI4KkIiIiDiOCpCIiIg4jgqQiIiIOI4KkIiIiDiOCpCIiIg4jgqQiDiKn58fb7/9dv7fv//+e9q0aUNwcDDNmzc/5TER8S3aC0xEfMIdd9zBSy+9BEBAQADVq1fniiuu4JZbbuGOO+6gQgXz+156ejrVqlXLf97o0aOpXLkyP/zwA1WqVDnlMRHxLRoBEhGf0bVrV9LT09mxYwcffPABHTt25OGHH6Zbt26cOHECgPDwcIKCgvKfs23bNtq3b0+9evWoUaPGKY+JiG9RARIRnxEUFER4eDh16tShZcuWDB8+nHfeeYcPPviA+fPnAwUvgfn5+bF+/XrGjBmDn58fjz/+eJHHRMT3qACJiE+75ppraNasGUuXLi30b+np6Vx++eU88sgjpKenM2TIkCKPiYjv0RwgEfF5l112GV9//XWh4+Hh4QQEBFClShXCw8MBqFKlSqFjIuJ7NAIkIj7P7Xbj5+dnO4aIeBAVIBHxeVu2bCEqKsp2DBHxICpAIuLTli9fzjfffEOfPn1sRxERD6I5QCLiM7Kzs8nIyCA3N5dffvmFDz/8kAkTJtCtWzduv/122/FExIOoAImIz/jwww+JiIggICCAatWq0axZM6ZNm8aAAQPyF0IUEQHwc7vdbtshRERERMqTfiUSERERx1EBEhEREcdRARIRERHHUQESERERx1EBEhEREcdRARIRERHHUQESERERx1EBEhEREcdRARIRERHHUQESERERx1EBEhEREcf5/+qHDCajm6zOAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.kdeplot(birthOrder[\"Diff\"], color='b', shade=True)\n", "sns.rugplot(birthOrder[\"Diff\"], height=0.1, color='b')" ] }, { "cell_type": "markdown", "id": "41d53375", "metadata": {}, "source": [ "The distribution does not look very Normal, with a hint of bimodality (two peaks).\n", "\n", "Since there is also no theoretical reason to think the differences should be normally distributed, we will use a non-parametric test" ] }, { "cell_type": "markdown", "id": "30e5dff0", "metadata": {}, "source": [ "### Hypotheses\n", "\n", "Ho: the median difference in independence between first- and second-born siblings is is zero\n", "\n", "Ha: the median difference in independence is not zero\n", " \n", "This is a two-tailed test as the researcher's hypothesis (described above) is not directional.\n", "\n", "We will test at the $\\alpha = 0.05$ significance level" ] }, { "cell_type": "markdown", "id": "0ed79570", "metadata": {}, "source": [ "### Descriptive statistics\n", "\n", "We obtain the relevant descriptive statistics. For the t-test we reported the summary statistics (mean, s.d. and $n$) that went into the formula for $t$. The same approach won't quite work here as the test is not based on summary statistics. However, common sense suggests that reporting the median for each group, a measure of spread based on ranks (namely the IQR), and the sample size would give the reader an idea of the data.\n" ] }, { "cell_type": "code", "execution_count": 19, "id": "11de87b4", "metadata": {}, "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", "
FirstBornSecondBornDiff
count20.00000020.00000020.000000
mean12.60000010.4500002.150000
std5.3646012.4381834.120232
min2.0000007.000000-5.000000
25%8.0000008.000000-2.000000
50%13.50000010.0000003.000000
75%17.00000012.2500005.250000
max20.00000015.00000010.000000
\n", "
" ], "text/plain": [ " FirstBorn SecondBorn Diff\n", "count 20.000000 20.000000 20.000000\n", "mean 12.600000 10.450000 2.150000\n", "std 5.364601 2.438183 4.120232\n", "min 2.000000 7.000000 -5.000000\n", "25% 8.000000 8.000000 -2.000000\n", "50% 13.500000 10.000000 3.000000\n", "75% 17.000000 12.250000 5.250000\n", "max 20.000000 15.000000 10.000000" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "birthOrder.describe()" ] }, { "cell_type": "markdown", "id": "b7ab9c7d", "metadata": {}, "source": [ "### Carry out the test\n", "\n", "We carry out the test using the function wilcoxon from scipy.stats, here loaded as stats" ] }, { "cell_type": "code", "execution_count": 20, "id": "85858db8", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "WilcoxonResult(statistic=46.0, pvalue=0.026641845703125)" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "stats.wilcoxon(birthOrder['FirstBorn'],birthOrder['SecondBorn'],alternative='two-sided')\n", "#help(stats.wilcoxon)" ] }, { "cell_type": "markdown", "id": "3dcd664c", "metadata": {}, "source": [ "The inputs to stats.wilcoxon are:\n", "
    \n", "
  • the two samples to be compared (the values of FirstBorn and SecondBorn from our Pandas data frame birthOrder)\n", "
  • the argument alternative='greater', which tells the computer to run a one tailed test \n", "that mean of the first input (FirstBorn) is greater than the second (SecondBorn).\n", "
\n", " \n", "The outputs are a value of the test statistic ($T=164$) and pvalue ($p=0.0133$) - if this is less than our $\\alpha$ value 0.5, there is a significant difference.\n", "\n", "More explanation of how T is calculated below." ] }, { "cell_type": "markdown", "id": "69dcbf20", "metadata": {}, "source": [ "### Draw conclusions\n", "\n", "As the p value of 0.0133 is less than our alpha value of 0.05, the test is significant. \n", "\n", "We can conclude that the median difference in idenpendence is positive, ie the first borns are more independent" ] }, { "cell_type": "markdown", "id": "1a80a337", "metadata": {}, "source": [ "# How the Wilcoxon Sign-Rank test works\n", "\n", "You have seen how to carry out the Sign-Rank test using scipy.stats but you may be none the wiser about how the computer arrived at the test statistic and p value.\n", "\n", "In this section we will build our own version of the test step by step to understand how it worked.\n", "\n", "### How to do the test (if you were doing it with pencil and paper)\n", "\n", "1. Obtain the difference (in independence score) for each pair\n", "\n", "1. Rank the differences regardless of sign (e.g. a difference of +4 is greater than a difference of -3, which is greater than a difference of +2). Remove pairs with zero difference\n", "\n", "1. Calculate the sum of ranks assigned to pairs with a positive difference (first-born more independent than second-born) - this is $R+$\n", "1. Calculate the sum of ranks assigned to pairs with a negative difference (first-born more independent than second-born) - this is $R-$\n", "\n", "1. The test statistic $T$ is either:\n", " * $R+$ if we expect positive differences to have the larger ranks (in this case, that equates to expecting first-borns to have higher scores)\n", " * $R-$ if we expect negative differences to have the larger ranks (in this case, that equates to expecting second-borns to have higher scores)\n", " * The smaller of $R+$ and $R-$ for a two tailed test (as in the example, we have no a-prior hypothesis about direction of effect)\n", "\n", "1. $T$ is compared with a null distribution (the expected distribubtion of $T$ obtained in samples drawn from a population in which there is no true difference between groups)\n", "\n", "\n", "### Step 1: Obtain the differences\n", "\n", "We already did this and added a column 'Rank' do our dataframe\n", "\n", "### Step 2: Rank the differences regardless of sign\n", "\n", "We can obtain the absolute differences and save them in a new column of the dataframe as follows\n" ] }, { "cell_type": "code", "execution_count": 21, "id": "fd1031e3", "metadata": {}, "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", " \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", "
FirstBornSecondBornDiffAbsDiff
0121022
1181266
21315-22
3171344
489-11
5151233
6161333
758-33
8810-22
912844
1013855
1159-44
1214866
1320101010
14191455
15171166
1627-55
1757-22
18151322
19181266
\n", "
" ], "text/plain": [ " FirstBorn SecondBorn Diff AbsDiff\n", "0 12 10 2 2\n", "1 18 12 6 6\n", "2 13 15 -2 2\n", "3 17 13 4 4\n", "4 8 9 -1 1\n", "5 15 12 3 3\n", "6 16 13 3 3\n", "7 5 8 -3 3\n", "8 8 10 -2 2\n", "9 12 8 4 4\n", "10 13 8 5 5\n", "11 5 9 -4 4\n", "12 14 8 6 6\n", "13 20 10 10 10\n", "14 19 14 5 5\n", "15 17 11 6 6\n", "16 2 7 -5 5\n", "17 5 7 -2 2\n", "18 15 13 2 2\n", "19 18 12 6 6" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "birthOrder['AbsDiff'] = birthOrder['Diff'].abs()\n", "birthOrder" ] }, { "cell_type": "markdown", "id": "651ea1b0", "metadata": {}, "source": [ "Then we rank the absolute differences" ] }, { "cell_type": "code", "execution_count": 22, "id": "ff365ea2", "metadata": {}, "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", " \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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
FirstBornSecondBornDiffAbsDiffRank
01210224.0
118126617.5
21315-224.0
317134411.0
489-111.0
51512338.0
61613338.0
758-338.0
8810-224.0
91284411.0
101385514.0
1159-4411.0
121486617.5
132010101020.0
1419145514.0
1517116617.5
1627-5514.0
1757-224.0
181513224.0
1918126617.5
\n", "
" ], "text/plain": [ " FirstBorn SecondBorn Diff AbsDiff Rank\n", "0 12 10 2 2 4.0\n", "1 18 12 6 6 17.5\n", "2 13 15 -2 2 4.0\n", "3 17 13 4 4 11.0\n", "4 8 9 -1 1 1.0\n", "5 15 12 3 3 8.0\n", "6 16 13 3 3 8.0\n", "7 5 8 -3 3 8.0\n", "8 8 10 -2 2 4.0\n", "9 12 8 4 4 11.0\n", "10 13 8 5 5 14.0\n", "11 5 9 -4 4 11.0\n", "12 14 8 6 6 17.5\n", "13 20 10 10 10 20.0\n", "14 19 14 5 5 14.0\n", "15 17 11 6 6 17.5\n", "16 2 7 -5 5 14.0\n", "17 5 7 -2 2 4.0\n", "18 15 13 2 2 4.0\n", "19 18 12 6 6 17.5" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "birthOrder['Rank'] = birthOrder['AbsDiff'].rank()\n", "birthOrder" ] }, { "cell_type": "markdown", "id": "6e5c0c58", "metadata": {}, "source": [ "... phew! Let's just check in on our understanding:\n", "
    \n", "
  • In the largest ranked pair (rank 20), the difference is +10 (the first-born's independence score was 10 points higher than the second-born's).\n", "\n", "
  • Several pairs with a difference of 4 points between siblings share rank 11. In some of these pairs the first-orn was more independent, in some pairs the second-born was more independent. \n", "
      \n", "
    • Find them in the table.\n", "
    \n", "
  • In 13/20 pairs the first-born was more independent; in 7/20 pairs the second-born was more independent.\n", "
" ] }, { "cell_type": "markdown", "id": "dbdc69a5", "metadata": {}, "source": [ "### Work out the test-statistic $T$\n", "\n", "We need to add up all the ranks assigned to pairs with positive differences (first-born is more independent)\n", " to get $R+$ and all the ranks assigned to pairs with negative differences (second-born is more independent)\n", " to get $R-$ \n" ] }, { "cell_type": "code", "execution_count": 23, "id": "ee326925", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "R+ = 164.0\n", "R- = 46.0\n" ] } ], "source": [ "Rplus = sum(birthOrder.loc[(birthOrder['Diff']>0), 'Rank'])\n", "Rminus = sum(birthOrder.loc[(birthOrder['Diff']<0), 'Rank'])\n", "\n", "print('R+ = ' + str(Rplus))\n", "print('R- = ' + str(Rminus))" ] }, { "cell_type": "markdown", "id": "e6e4d1f8", "metadata": {}, "source": [ "T is the smaller of these - you might like to check that the value obtained below matches the 'test statistic' from the function scipy.stats.wilcoxon that we used above" ] }, { "cell_type": "code", "execution_count": 24, "id": "f3c3a51f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "T = 46.0\n" ] } ], "source": [ "T = min(Rplus, Rminus)\n", "print('T = ' + str(T))" ] }, { "cell_type": "markdown", "id": "359d9744", "metadata": {}, "source": [ "### Establish the null distriubtion\n", "\n", "To convert our $T$ to a $p$-value, we need to know how the probability of obtaining that value of $T$ due to chance if the null hypothesis were true.\n", "\n", "### What do we mean by \"If the null were true\"?\n", "\n", "If the null were true, the first- and second-born siblings should be equally likely to be more independent. If this were true, high- and low-ranked differences would occur randomly in favour of the first-born and second-born sibling. So each rank (in this case, the ranks 1-20, as we have 20 sibling pairs) might equally likely feed into $R+$ or $R-$\n", "\n", "### Estalbishing the null by complete enumeration\n", "\n", "For a much smaller sample size, we could work out all the possible ways the ranks could be randomly assigned to contribute to $R+$ or $R-$, and work out in what proportion of cases this simulated $U$ was as large as, or larger than, the value of $U$ from our experiment.\n", "\n", "This is not practical for larger samples as there are too many possible ways to deal the ranks (too many combbinations). In the example above where $n_1 = 9$ and $n_2 = 11$, the number of possible deals or combinations is 167960 - for larger $n$ the number of combinations becomes astronomical.\n", "\n", "### Establishing the null by simulation\n", "\n", "We can work out the null approximately by doing a lot of random 'deals' of the ranks to $R+$ or $R-$, calculating $T$ for each case, and work out in what proportion of cases this simulated $T$ was as large as, or larger than, the value of $T$ from our experiment.\n", "\n", "\"A lot\" might be 10,000 or more (this is still much lower than the number of possible deals under complete enumeration)\n", "\n", "Let's try it!" ] }, { "cell_type": "code", "execution_count": 25, "id": "61629cee", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,\n", " 18, 19, 20])" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Possible ranks are the numbers 1 to 20 (remember Python uses the counterintuitive syntax [1,21))\n", "ranks = np.arange(1,21)\n", "ranks" ] }, { "cell_type": "code", "execution_count": 26, "id": "620f262f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ True, False, False, True, False, True, False, False, False,\n", " False, True, True, True, False, False, False, True, True,\n", " True, True])" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# flip 20 virtual coind to assign each of these 20 ranks to R+ or R- \n", "isPos = (np.random.uniform(0,1,20)>0.5)\n", "isPos" ] }, { "cell_type": "code", "execution_count": 27, "id": "93ce8f34", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "R+ = 121\n", "R- = 89\n" ] } ], "source": [ "# calculate R+ and R-\n", "Rplus=sum(ranks[isPos])\n", "Rminus=sum(ranks[~isPos])\n", "\n", "print('R+ = ' + str(Rplus))\n", "print('R- = ' + str(Rminus))" ] }, { "cell_type": "code", "execution_count": 28, "id": "cd104a8d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "T = 121\n" ] } ], "source": [ "# T is the smaller of R+,R-\n", "T = min(Rplus, Rminus)\n", "print('T = ' + str(T))" ] }, { "cell_type": "markdown", "id": "5a1b1e46", "metadata": {}, "source": [ "Phew! That was one simulated sample. Now we will make a loop to do it 10 000 times!" ] }, { "cell_type": "code", "execution_count": 34, "id": "31dde7fd", "metadata": {}, "outputs": [], "source": [ "n = 20\n", "ranks = np.arange(1,(n+1))\n", "maxT = sum(ranks)\n", "\n", "nReps = 10000\n", "simulated_T = np.empty(nReps)\n", "\n", "for i in range(nReps):\n", " isPos = (np.random.uniform(0,1,20)>0.5)\n", " Rplus=sum(ranks[isPos])\n", " Rminus=sum(ranks[~isPos])\n", " T = min(Rplus, Rminus)\n", " simulated_T[i]=T" ] }, { "cell_type": "markdown", "id": "100c10aa", "metadata": {}, "source": [ "Now we can plot a histogram of the simulated values of T - note that the distribution looks like it has been 'sliced in half', because of the way we calculate $T$ (it is always the smaller of the two rank-sums)\n", "\n", "We can also calculate the p-value of our observed value of T, $T=46$ by working out the proportion of simmulated samples in which T<=46" ] }, { "cell_type": "code", "execution_count": 35, "id": "7e8ac4b7", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGyCAYAAAACgQXWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABFD0lEQVR4nO3de1hVdd7//9eWk4CIgsqGQqXCTPEUzviVmjyCmofMJitttLRuHfOAhyjHSiwFtRGY9M4O4y2m49BR75qZSjSzyJoURcX6ZgfGU5soI5BEIPj8/ujX/rZFkQ0bge3zcV3rutqf9d5rv9earnjNZ50sxhgjAAAAN9WisRsAAABoSIQdAADg1gg7AADArRF2AACAWyPsAAAAt0bYAQAAbo2wAwAA3BphBwAAuDXCDgAAcGuejd1AU1BVVaWvv/5aAQEBslgsjd0OAACoBWOMTp8+rbCwMLVoUcP8jWlEFRUVZtGiRaZz586mZcuWJiIiwixZssRUVlbaa6qqqszixYtNaGioadmypRkwYIDJzc112M7Zs2fNzJkzTXBwsPHz8zOjR482x48fr3Ufx48fN5JYWFhYWFhYmuFysb/5jTqzs2LFCj3zzDPasGGDunfvrr179+ree+9VYGCg5syZI0lauXKlUlJSlJ6eri5dumjp0qWKjY3VZ599poCAAElSfHy83njjDWVkZCg4OFjz58/XqFGjlJ2dLQ8Pj4v28ct2jh8/rtatWzfcDgMAAJcpLi5WeHi4/e/4hViMabwXgY4aNUohISFat26dfey2226Tn5+fNm7cKGOMwsLCFB8fr4ceekiSVFZWppCQEK1YsULTpk1TUVGR2rdvr40bN+qOO+6QJH399dcKDw/Xv/71Lw0bNuyifRQXFyswMFBFRUWEHQAAmona/v1u1AuUb7zxRu3YsUNHjhyRJB04cEBZWVm6+eabJUl5eXnKz89XXFyc/Ts+Pj4aMGCAdu/eLUnKzs5WRUWFQ01YWJiioqLsNecqKytTcXGxwwIAANxTo57Geuihh1RUVKSuXbvKw8NDlZWVWrZsme666y5JUn5+viQpJCTE4XshISE6evSovcbb21tt27atVvPL98+VnJysJUuWuHp3AABAE9SoMzsvvviiNm3apM2bN2vfvn3asGGD/vznP2vDhg0OdefeIWWMuehdUzXVLFy4UEVFRfbl+PHj9dsRAADQZDXqzM6DDz6ohx9+WHfeeackqUePHjp69KiSk5M1efJkWa1WST/P3oSGhtq/V1BQYJ/tsVqtKi8vV2FhocPsTkFBgWJiYs77uz4+PvLx8Wmo3QIAAE1Io87snDlzptp98R4eHqqqqpIkRUREyGq1KjMz076+vLxcu3btsgeZ6OhoeXl5OdTYbDbl5uZeMOwAAIDLR6PO7IwePVrLli1Tx44d1b17d+3fv18pKSmaMmWKpJ9PX8XHxyspKUmRkZGKjIxUUlKS/Pz8NGHCBElSYGCgpk6dqvnz5ys4OFhBQUFasGCBevTooaFDhzbm7gEAgCagUcPO6tWr9eijj2rGjBkqKChQWFiYpk2bpscee8xek5CQoNLSUs2YMUOFhYXq16+ftm3b5nBPfWpqqjw9PTV+/HiVlpZqyJAhSk9Pr9UzdgAAgHtr1OfsNBU8ZwcAgOanWTxnBwAAoKERdgAAgFsj7AAAALdG2AEAAG6NsAMAANwaYQcAALg1wg4AAHBrjfpQQQAA0Lx179lbNputxprQ0FAdPphzaRo6D8IOAACoM5vNprhlW2us2bZo7CXp5UI4jQUAANwaYQcAALg1wg4AAHBrhB0AAODWCDsAAMCtEXYAAIBbI+wAAAC3RtgBAABujbADAADcGmEHAAC4NcIOAABwa4QdAADg1gg7AADArRF2AACAWyPsAAAAt0bYAQAAbo2wAwAA3BphBwAAuDXCDgAAcGuEHQAA4NYIOwAAwK0RdgAAgFsj7AAAALdG2AEAAG6NsAMAANwaYQcAALi1Rg07nTt3lsViqbY88MADkiRjjBITExUWFiZfX18NHDhQhw8fdthGWVmZZs2apXbt2snf319jxozRiRMnGmN3AABAE9SoYWfPnj2y2Wz2JTMzU5J0++23S5JWrlyplJQUrVmzRnv27JHValVsbKxOnz5t30Z8fLy2bNmijIwMZWVlqaSkRKNGjVJlZWWj7BMAAGhaGjXstG/fXlar1b784x//0NVXX60BAwbIGKO0tDQtWrRI48aNU1RUlDZs2KAzZ85o8+bNkqSioiKtW7dOq1at0tChQ9WnTx9t2rRJhw4d0vbt2xtz1wAAQBPRZK7ZKS8v16ZNmzRlyhRZLBbl5eUpPz9fcXFx9hofHx8NGDBAu3fvliRlZ2eroqLCoSYsLExRUVH2mvMpKytTcXGxwwIAANxTkwk7W7du1Q8//KB77rlHkpSfny9JCgkJcagLCQmxr8vPz5e3t7fatm17wZrzSU5OVmBgoH0JDw934Z4AAICmpMmEnXXr1mnEiBEKCwtzGLdYLA6fjTHVxs51sZqFCxeqqKjIvhw/frzujQMAgCatSYSdo0ePavv27brvvvvsY1arVZKqzdAUFBTYZ3usVqvKy8tVWFh4wZrz8fHxUevWrR0WAADgnppE2Fm/fr06dOigkSNH2sciIiJktVrtd2hJP1/Xs2vXLsXExEiSoqOj5eXl5VBjs9mUm5trrwEAAJc3z8ZuoKqqSuvXr9fkyZPl6fn/2rFYLIqPj1dSUpIiIyMVGRmppKQk+fn5acKECZKkwMBATZ06VfPnz1dwcLCCgoK0YMEC9ejRQ0OHDm2sXQIAAE1Io4ed7du369ixY5oyZUq1dQkJCSotLdWMGTNUWFiofv36adu2bQoICLDXpKamytPTU+PHj1dpaamGDBmi9PR0eXh4XMrdAAAATZTFGGMau4nGVlxcrMDAQBUVFXH9DgAATghqH6K4ZVtrrNm2aKy+//Ybl/92bf9+N4lrdgAAABoKYQcAALg1wg4AAHBrhB0AAODWCDsAAMCtEXYAAIBbI+wAAAC3RtgBAABujbADAADcGmEHAAC4NcIOAABwa4QdAADg1gg7AADArRF2AACAWyPsAAAAt0bYAQAAbo2wAwAA3BphBwAAuDXCDgAAcGuEHQAA4NYIOwAAwK0RdgAAgFsj7AAAALdG2AEAAG6NsAMAANwaYQcAALg1wg4AAHBrhB0AAODWCDsAAMCtEXYAAIBbI+wAAAC35nTY2bBhg/75z3/aPyckJKhNmzaKiYnR0aNHXdocAABAfTkddpKSkuTr6ytJ+vDDD7VmzRqtXLlS7dq109y5c13eIAAAQH14OvuF48eP65prrpEkbd26Vb///e/1X//1X7rhhhs0cOBAV/cHAABQL07P7LRq1UqnTp2SJG3btk1Dhw6VJLVs2VKlpaVON3Dy5EndfffdCg4Olp+fn3r37q3s7Gz7emOMEhMTFRYWJl9fXw0cOFCHDx922EZZWZlmzZqldu3ayd/fX2PGjNGJEyec7gUAALgfp8NObGys7rvvPt133306cuSIRo4cKUk6fPiwOnfu7NS2CgsLdcMNN8jLy0tvvvmmPvnkE61atUpt2rSx16xcuVIpKSlas2aN9uzZI6vVqtjYWJ0+fdpeEx8fry1btigjI0NZWVkqKSnRqFGjVFlZ6ezuAQAAN+P0aaz//u//1iOPPKLjx4/r1VdfVXBwsCQpOztbd911l1PbWrFihcLDw7V+/Xr72K8DkzFGaWlpWrRokcaNGyfp5wukQ0JCtHnzZk2bNk1FRUVat26dNm7caJ9l2rRpk8LDw7V9+3YNGzbM2V0EAABuxOmZnTZt2mjNmjX63//9Xw0fPtw+vmTJEi1atMipbb3++uvq27evbr/9dnXo0EF9+vTR888/b1+fl5en/Px8xcXF2cd8fHw0YMAA7d69W9LPIauiosKhJiwsTFFRUfaac5WVlam4uNhhAQAA7qlOz9l5//33dffddysmJkYnT56UJG3cuFFZWVlObeerr77S2rVrFRkZqbffflvTp0/X7Nmz9cILL0iS8vPzJUkhISEO3wsJCbGvy8/Pl7e3t9q2bXvBmnMlJycrMDDQvoSHhzvVNwAAaD6cDjuvvvqqhg0bJl9fX+3bt09lZWWSpNOnTyspKcmpbVVVVen6669XUlKS+vTpo2nTpun+++/X2rVrHeosFovDZ2NMtbFz1VSzcOFCFRUV2Zfjx4871TcAAGg+nA47S5cu1TPPPKPnn39eXl5e9vGYmBjt27fPqW2FhoaqW7duDmPXXXedjh07JkmyWq2SVG2GpqCgwD7bY7VaVV5ersLCwgvWnMvHx0etW7d2WAAAuJx079lbQe1Daly69+zd2G26hNMXKH/22We66aabqo23bt1aP/zwg1PbuuGGG/TZZ585jB05ckSdOnWSJEVERMhqtSozM1N9+vSRJJWXl2vXrl1asWKFJCk6OlpeXl7KzMzU+PHjJUk2m025ublauXKls7sHAMBlwWazKW7Z1hprti0ae0l6aWhOh53Q0FB98cUX1W4zz8rK0lVXXeXUtubOnauYmBglJSVp/Pjx+vjjj/Xcc8/pueeek/Tz6av4+HglJSUpMjJSkZGRSkpKkp+fnyZMmCBJCgwM1NSpUzV//nwFBwcrKChICxYsUI8ePex3ZwEAgMuX02Fn2rRpmjNnjv7nf/5HFotFX3/9tT788EMtWLBAjz32mFPb+s1vfqMtW7Zo4cKFevzxxxUREaG0tDRNnDjRXpOQkKDS0lLNmDFDhYWF6tevn7Zt26aAgAB7TWpqqjw9PTV+/HiVlpZqyJAhSk9Pl4eHh7O7BwAA3IzFGGOc/dKiRYuUmpqqs2fPSvr5GpgFCxboiSeecHmDl0JxcbECAwNVVFTE9TsAgMtCUPuQWp3G+v7bby7Jduqitn+/nZ7ZkaRly5Zp0aJF+uSTT1RVVaVu3bqpVatWdW4WAACgodQp7EiSn5+f+vbt68peAAAAXK5WYWfcuHFKT09X69at7a9tuJDXXnvNJY0BAAC4Qq3CTmBgoP0BfYGBgQ3aEAAAl6vuPXvLZrPVWBMaGqrDB3MuTUNuolZh55cXdRpjlJiYqPbt28vPz69BGwMA4HJzOT375lJy6gnKxhhFRkba34cFAADQ1DkVdlq0aKHIyEidOnWqofoBAABwKaffjbVy5Uo9+OCDys3NbYh+AAAAXMrpW8/vvvtunTlzRr169ZK3t7d8fX0d1n///fcuaw4AAKC+nA47qamp9juzAAAAmjqnw84999zTAG0AAAA0DKev2fHw8FBBQUG18VOnTvHiTQAA0OQ4HXYu9N7QsrIyeXt717shAAAAV6r1aaynnnpKkmSxWPTXv/7V4cWflZWVeu+999S1a1fXdwgAAFAPtQ47qampkn6e2XnmmWccTll5e3urc+fOeuaZZ1zfIQAAQD3UOuzk5eVJkgYNGqTXXntNbdu2bbCmAABA3dXmHVvFp09fdDvFp0sU1D6k3ttpbE7fjbVz586G6AMAALhIbd6x9fLMwRfdjqmqcsl2GpvTFygDAAA0J07P7AAAgMbjLqeWLiXCDgAAzYi7nFq6lDiNBQAA3JrTMzuff/65du/erfz8fFksFoWEhCgmJkaRkZEN0R8AAEC91DrsFBUVadKkSXrjjTcUGBioDh06yBijb7/9VsXFxRo9erReeOEFtW7duiH7BQAAcEqtT2PNmjVLeXl5+vDDD1VYWKjPPvtMR44cUWFhoXbv3q28vDzNmjWrIXsFAABwWq1ndl5//XW9/fbb6tevX7V1/fr107PPPqvhw4e7tDkAAID6cuoCZYvFUqd1AAAAjaXWYWf06NG6//77tXfv3mrr9u7dq+nTp2vMmDEubQ4AAKC+ah12Vq9erbCwMP32t79VUFCQunbtquuuu05BQUHq16+fQkND7W9GBwAAaCpqfc1OmzZt9Oabb+rTTz/VRx99pPz8fEmS1WpV//791bVr1wZrEgAAoK6cfs7Oddddp+uuu64hegEAAHA5p8KOMUbbt2+v9lDBG264QUOGDOEiZQAA0OTU+pqdkydP6vrrr9eIESO0ZcsWffXVV/riiy+0ZcsWDR8+XH379tXJkycbslcAAACn1XpmZ8aMGQoKCtLx48cVGhrqsM5ms+nuu+/WAw88oK1bt7q6RwAAgDqrddjZsWOHPvjgg2pBR5JCQ0P15z//Wb/73e9c2hwAAO6ie8/estlsNdYUnz59ibq5vNQ67Pj6+ur777+/4PrCwkL5+vq6pCkAANyNzWZT3LKtNda8PHPwpWnmMlPra3buvPNOTZ48Wa+88oqKiors40VFRXrllVd07733asKECU79eGJioiwWi8NitVrt640xSkxMVFhYmHx9fTVw4EAdPnzYYRtlZWWaNWuW2rVrJ39/f40ZM0YnTpxwqg8AAOC+ah12Vq1apZEjR2rixIkKCgqSr6+vfH19FRQUpIkTJ2rkyJF68sknnW6ge/fustls9uXQoUP2dStXrlRKSorWrFmjPXv2yGq1KjY2Vqd/Nc0XHx+vLVu2KCMjQ1lZWSopKdGoUaNUWVnpdC8AAMD91Po0lre3t9auXasVK1Zo7969+uabbyT9/FDB6OhotW7dum4NeHo6zOb8whijtLQ0LVq0SOPGjZMkbdiwQSEhIdq8ebOmTZumoqIirVu3Ths3btTQoUMlSZs2bVJ4eLi2b9+uYcOG1aknAADgPpx+qGDr1q01eLDrzil+/vnnCgsLk4+Pj/r166ekpCRdddVVysvLU35+vuLi4uy1Pj4+GjBggHbv3q1p06YpOztbFRUVDjVhYWGKiorS7t27Lxh2ysrKVFZWZv9cXFzssv0BALiX2lxYHBoaqsMHcy5NQ3CaU2Hnxx9/1ObNm8/7UMG77rpL/v7+Tv14v3799MILL6hLly765ptvtHTpUsXExOjw4cP211GEhIQ4fCckJERHjx6VJOXn58vb21tt27atVvPL988nOTlZS5YscapXAMDlqTYXFm9bNPaS9IK6qfU1O5988om6dOmihIQEFRYWqmPHjrryyitVWFioBx98UNdee60++eQTp358xIgRuu2229SjRw8NHTpU//znPyX9fLrqF+c+ldkYc9EnNV+sZuHChSoqKrIvx48fd6pvAADQfNR6ZueBBx7QTTfdpA0bNsjb29thXXl5ue655x498MAD2rlzZ52b8ff3V48ePfT5559r7Nixkn6evfn1s30KCgrssz1Wq1Xl5eUqLCx0mN0pKChQTEzMBX/Hx8dHPj4+de4TAAA0H7We2fn3v/+tRx99tFrQkX6+ePlPf/qT/v3vf9ermbKyMn366acKDQ1VRESErFarMjMz7evLy8u1a9cue5CJjo6Wl5eXQ43NZlNubm6NYQcAAFw+aj2z07ZtW33++efq1q3bedd/8cUX1a6duZgFCxZo9OjR6tixowoKCrR06VIVFxdr8uTJslgsio+PV1JSkiIjIxUZGamkpCT5+fnZn+cTGBioqVOnav78+QoODlZQUJAWLFhgPy0GAMClUHy6REHtQy5Sw9ORG0utw87999+vyZMn65FHHlFsbKxCQkJksViUn5+vzMxMJSUlKT4+3qkfP3HihO666y599913at++vf7P//k/+uijj9SpUydJUkJCgkpLSzVjxgwVFhaqX79+2rZtmwICAuzbSE1Nlaenp8aPH6/S0lINGTJE6enp8vDwcKoXAADqylRV8XTkJqzWYScxMVG+vr5KSUlRQkKC/QJgY4ysVqsefvhhJSQkOPXjGRkZNa63WCxKTExUYmLiBWtatmyp1atXa/Xq1U79NgAAuDw4dev5Qw89pIceesj+DBzp54uEIyIiGqQ5AACA+nL6oYKSFBERQcABAADNglNh58SJE1q7dm21hwrGxMRo+vTpCg8Pb6g+AQAA6qTWt55nZWXpuuuu05YtW9SrVy9NmjRJd999t3r16qWtW7eqe/fu+uCDDxqyVwAAAKfVemZn7ty5uu+++5SamnrB9fHx8dqzZ4/LmgMAAKivWs/s5Obmavr06RdcP23aNOXm5rqkKQAAAFepddgJDQ3V7t27L7j+ww8/dHitAwAAQFNQ69NYCxYs0PTp05WdnX3ehwr+9a9/VVpaWgO2CgAA4Lxah50ZM2YoODhYqampevbZZ1VZWSlJ8vDwUHR0tF544QWNHz++wRoFAMDVuvfsLZvNVmMNr3lo/py69fyOO+7QHXfcoYqKCn333XeSpHbt2snLy6tBmgMAoCHZbDZe83AZqNNDBb28vLg+BwAANAu1vkD5Yr788ksNHkz6BQAATYvLwk5JSYl27drlqs0BAAC4RK1PYz311FM1rj958mS9mwEAAHC1Woed+Ph4hYaGytvb+7zry8vLXdYUAACAq9Q67HTq1EkrVqy44O3lOTk5io6OdlljAAAArlDra3aio6OVnZ19wfUWi0XGGJc0BQAA4Cq1ntl5/PHHdebMmQuu79atm/Ly8lzSFAAAgKvUOux069atxvVeXl7q1KlTvRsCAABwpTo9VBAAgKaOV0HgF4QdAIBb4lUQ+IXLHioIAADQFBF2AACAW6tV2AkKCrK/5XzKlCk6zTlOAADQTNQq7JSXl6u4uFiStGHDBp09e7ZBmwIAAHCVWl2g3L9/f40dO1bR0dEyxmj27Nny9fU9b+3//M//uLRBAACA+qhV2Nm0aZNSU1P15ZdfymKxqKioiNkdAADQLNQq7ISEhGj58uWSpIiICG3cuFHBwcEN2hgAAIArOP2cHV4JAQAAmpM63Xq+a9cujR49Wtdcc40iIyM1ZswYvf/++67uDQAAoN6cDjubNm3S0KFD5efnp9mzZ2vmzJny9fXVkCFDtHnz5oboEQAAoM6cPo21bNkyrVy5UnPnzrWPzZkzRykpKXriiSc0YcIElzYIAABQH07P7Hz11VcaPXp0tfExY8ZwPQ8AAGhynA474eHh2rFjR7XxHTt2KDw83CVNAQAAuIrTp7Hmz5+v2bNnKycnRzExMbJYLMrKylJ6err+8pe/NESPAAAAdeb0zM4f//hHZWRk6NChQ4qPj9ecOXOUm5urF198UdOmTatzI8nJybJYLIqPj7ePGWOUmJiosLAw+fr6auDAgTp8+LDD98rKyjRr1iy1a9dO/v7+GjNmjE6cOFHnPgAAgHup063nt956q7KysnTq1CmdOnVKWVlZuuWWW+rcxJ49e/Tcc8+pZ8+eDuMrV65USkqK1qxZoz179shqtSo2NtbhRaTx8fHasmWLMjIylJWVpZKSEo0aNUqVlZV17gcAALiPOoUdVyopKdHEiRP1/PPPq23btvZxY4zS0tK0aNEijRs3TlFRUdqwYYPOnDljv8W9qKhI69at06pVqzR06FD16dNHmzZt0qFDh7R9+/bG2iUAANCENHrYeeCBBzRy5EgNHTrUYTwvL0/5+fmKi4uzj/n4+GjAgAHavXu3JCk7O1sVFRUONWFhYYqKirLXnE9ZWZmKi4sdFgAA4J6cvkDZlTIyMrRv3z7t2bOn2rr8/HxJP7+X69dCQkJ09OhRe423t7fDjNAvNb98/3ySk5O1ZMmS+rYPAGgk3Xv2ls1mq7Gm+FeXPODy1mhh5/jx45ozZ462bdumli1bXrDOYrE4fDbGVBs718VqFi5cqHnz5tk/FxcXc9s8ADQjNptNccu21ljz8szBl6YZNHlOn8Z6/PHHdebMmWrjpaWlevzxx2u9nezsbBUUFCg6Olqenp7y9PTUrl279NRTT8nT09M+o3PuDE1BQYF9ndVqVXl5uQoLCy9Ycz4+Pj5q3bq1wwIAANyT02FnyZIlKikpqTZ+5swZp04NDRkyRIcOHVJOTo596du3ryZOnKicnBxdddVVslqtyszMtH+nvLxcu3btUkxMjCQpOjpaXl5eDjU2m025ubn2GgAAcHlz+jTWhU4RHThwQEFBQbXeTkBAgKKiohzG/P39FRwcbB+Pj49XUlKSIiMjFRkZqaSkJPn5+dnfvxUYGKipU6dq/vz5Cg4OVlBQkBYsWKAePXpUu+AZAABcnmoddtq2bSuLxSKLxaIuXbo4BJ7KykqVlJRo+vTpLm0uISFBpaWlmjFjhgoLC9WvXz9t27ZNAQEB9prU1FR5enpq/PjxKi0t1ZAhQ5Seni4PDw+X9gIAAJqnWoedtLQ0GWM0ZcoULVmyRIGBgfZ13t7e6ty5s/r371+vZt59912HzxaLRYmJiUpMTLzgd1q2bKnVq1dr9erV9fptAADgnmoddiZPnixJioiIUExMjLy8vBqsKQAAAFdx+pqdiIiIGp9t0LFjx3o1BAAA4EpOh53OnTvX+Awb3kkFAACaEqfDzv79+x0+V1RUaP/+/UpJSdGyZctc1hgAAIArOB12evXqVW2sb9++CgsL05NPPqlx48a5pDEAAABXcNmLQLt06XLed1wBAAA0Jqdnds59Q7gxRjabTYmJiYqMjHRZYwAAAK7gdNhp06bNeV/OGR4eroyMDJc1BgAA4ApOh52dO3c6fG7RooXat2+va665Rp6ejfYSdQAAgPNyOp0MGDCgIfoAAABoEHWaivnyyy+VlpamTz/9VBaLRdddd53mzJmjq6++2tX9AQAA1IvTYeftt9/WmDFj1Lt3b91www0yxmj37t3q3r273njjDcXGxjZEnwAAN9C9Z+8an8IvSaGhoTp8MOfSNITLgtNh5+GHH9bcuXO1fPnyauMPPfQQYQcAcEE2m01xy7bWWLNt0dhL0gsuH06HnU8//VQvvfRStfEpU6YoLS3NFT0BAC5jxadLFNQ+5CI1py9RN3AHToed9u3bKycnp9ozdXJyctShQweXNQYAuDyZqqqLzv68PHPwpWkGbsHpsHP//ffrv/7rv/TVV18pJiZGFotFWVlZWrFihebPn98QPQIAANSZ02Hn0UcfVUBAgFatWqWFCxdKksLCwpSYmKjZs2e7vEEAAID6cDrsWCwWzZ07V3PnztXp//+caUBAgMsbAwAAcIV6PfKYkAMAAJo6l731HAAAoCki7AAAALfGmzsBAC5Rm6cj83wcNAanw05eXp4iIiIaohcAQDNWm6cj83wcNAanT2Ndc801GjRokDZt2qSzZ882RE8AAAAu43TYOXDggPr06aP58+fLarVq2rRp+vjjjxuiNwAAgHpzOuxERUUpJSVFJ0+e1Pr165Wfn68bb7xR3bt3V0pKir799tuG6BMAAKBO6nw3lqenp2699Va99NJLWrFihb788kstWLBAV155pSZNmnTRi9QAAAAuhTqHnb1792rGjBkKDQ1VSkqKFixYoC+//FLvvPOOTp48qVtuucWVfQIAANSJ03djpaSkaP369frss890880364UXXtDNN9+sFi1+zk0RERF69tln1bVrV5c3CwAA4Cynw87atWs1ZcoU3XvvvbJareet6dixo9atW1fv5gAAAOrL6bDz+eefX7TG29tbkydPrlNDAAAAruT0NTvr16/Xyy+/XG385Zdf1oYNG1zSFAAAgKs4HXaWL1+udu3aVRvv0KGDkpKSXNIUAACAqzh9Guvo0aPnfV1Ep06ddOzYMZc0BQBoWnjvFZozp8NOhw4ddPDgQXXu3Nlh/MCBAwoODnZVXwCAJoT3XqE5c/o01p133qnZs2dr586dqqysVGVlpd555x3NmTNHd955p1PbWrt2rXr27KnWrVurdevW6t+/v9588037emOMEhMTFRYWJl9fXw0cOFCHDx922EZZWZlmzZqldu3ayd/fX2PGjNGJEyec3S0AAOCmnA47S5cuVb9+/TRkyBD5+vrK19dXcXFxGjx4sNPX7Fx55ZVavny59u7dq71792rw4MG65ZZb7IFm5cqVSklJ0Zo1a7Rnzx5ZrVbFxsbq9K+mSuPj47VlyxZlZGQoKytLJSUlGjVqlCorK53dNQAA4IacPo3l7e2tF198UU888YQOHDggX19f9ejRQ506dXL6x0ePHu3wedmyZVq7dq0++ugjdevWTWlpaVq0aJHGjRsnSdqwYYNCQkK0efNmTZs2TUVFRVq3bp02btyooUOHSpI2bdqk8PBwbd++XcOGDXO6JwAA4F6cDju/6NKli7p06eKyRiorK/Xyyy/rxx9/VP/+/ZWXl6f8/HzFxcXZa3x8fDRgwADt3r1b06ZNU3Z2tioqKhxqwsLCFBUVpd27d18w7JSVlamsrMz+ubi42GX7AQAAmhanw05lZaXS09O1Y8cOFRQUqKqqymH9O++849T2Dh06pP79++vs2bNq1aqVtmzZom7dumn37t2SpJCQEIf6kJAQHT16VJKUn58vb29vtW3btlpNfn7+BX8zOTlZS5YscapPAADQPDkddubMmaP09HSNHDlSUVFRslgs9Wrg2muvVU5Ojn744Qe9+uqrmjx5snbt2mVff+72jTEX/c2L1SxcuFDz5s2zfy4uLlZ4eHgd9wAAADRlToedjIwMvfTSS7r55ptd0oC3t7euueYaSVLfvn21Z88e/eUvf9FDDz0k6efZm9DQUHt9QUGBfbbHarWqvLxchYWFDrM7BQUFiomJueBv+vj4yMfHxyX9AwCAps3pu7F+HU4agjFGZWVlioiIkNVqVWZmpn1deXm5du3aZQ8y0dHR8vLycqix2WzKzc2tMewAAIDLh9MzO/Pnz9df/vIXrVmzpt6nsP70pz9pxIgRCg8P1+nTp5WRkaF3331Xb731liwWi+Lj45WUlKTIyEhFRkYqKSlJfn5+mjBhgiQpMDBQU6dO1fz58xUcHKygoCAtWLBAPXr0sN+dBQAALm9Oh52srCzt3LlTb775prp37y4vLy+H9a+99lqtt/XNN9/oD3/4g2w2mwIDA9WzZ0+99dZbio2NlSQlJCSotLRUM2bMUGFhofr166dt27YpICDAvo3U1FR5enpq/PjxKi0t1ZAhQ5Seni4PDw9ndw0AALghp8NOmzZtdOutt7rkx9etW1fjeovFosTERCUmJl6wpmXLllq9erVWr17tkp4AAIB7cTrsrF+/viH6AAAAaBBOX6AsST/99JO2b9+uZ5991v7qhq+//lolJSUubQ4AAKC+nJ7ZOXr0qIYPH65jx46prKxMsbGxCggI0MqVK3X27Fk988wzDdEnAABAndTpoYJ9+/bVgQMHFBwcbB+/9dZbdd9997m0OQBA/XTv2Vs2m63GmtDQUB0+mHNpGgIaQZ3uxvrggw/k7e3tMN6pUyedPHnSZY0BAOrPZrMpbtnWGmu2LRp7SXoBGovT1+xUVVWpsrKy2viJEyccbgkHAABoCpwOO7GxsUpLS7N/tlgsKikp0eLFi132CgkAAABXcfo0VmpqqgYNGqRu3brp7NmzmjBhgj7//HO1a9dOf//73xuiRwAAgDpzOuyEhYUpJydHf//737Vv3z5VVVVp6tSpmjhxonx9fRuiRwAAgDpzOuxIkq+vr6ZMmaIpU6a4uh8AAACXcjrsvPDCCzWunzRpUp2bAQAAcLU6PWfn1yoqKnTmzBl5e3vLz8+PsAMAAJoUp+/GKiwsdFhKSkr02Wef6cYbb+QCZQAA0OTU6d1Y54qMjNTy5curzfoAAAA0NpeEHUny8PDQ119/7arNAQAAuITT1+y8/vrrDp+NMbLZbFqzZo1uuOEGlzUGAADgCk6HnbFjxzp8tlgsat++vQYPHqxVq1a5qi8AwCVSfLpEQe1DLlJz+hJ1A7ie02GnqqqqIfoAADQSU1V10ZeFvjxz8KVpBmgALrtmBwAAoClyemZn3rx5ta5NSUlxdvMAAAAu5XTY2b9/v/bt26effvpJ1157rSTpyJEj8vDw0PXXX2+vs1gsrusSAACgjpwOO6NHj1ZAQIA2bNigtm3bSvr5QYP33nuvfve732n+/PkubxIAAKCunL5mZ9WqVUpOTrYHHUlq27atli5dyt1YAACgyXE67BQXF+ubb76pNl5QUKDT3JoIAACaGKfDzq233qp7771Xr7zyik6cOKETJ07olVde0dSpUzVu3LiG6BEAAKDOnL5m55lnntGCBQt09913q6Ki4ueNeHpq6tSpevLJJ13eIAAAQH04HXb8/Pz09NNP68knn9SXX34pY4yuueYa+fv7N0R/AIAL6N6zt2w2W401PPkYqEPY+YXNZpPNZtNNN90kX19fGWO43RwALiGbzcaTj4FacPqanVOnTmnIkCHq0qWLbr75Zvv/q7jvvvu47RwAADQ5ToeduXPnysvLS8eOHZOfn599/I477tBbb73l0uYAAADqy+nTWNu2bdPbb7+tK6+80mE8MjJSR48edVljAAAAruD0zM6PP/7oMKPzi++++04+Pj4uaQoAAMBVnA47N910k1544QX7Z4vFoqqqKj355JMaNGiQS5sDAACoL6dPYz355JMaOHCg9u7dq/LyciUkJOjw4cP6/vvv9cEHHzREjwAAAHXm9MxOt27ddPDgQf32t79VbGysfvzxR40bN0779+/X1Vdf7dS2kpOT9Zvf/EYBAQHq0KGDxo4dq88++8yhxhijxMREhYWFydfXVwMHDtThw4cdasrKyjRr1iy1a9dO/v7+GjNmjE6cOOHsrgEAADfkVNipqKjQoEGDVFxcrCVLlugf//iH/vWvf2np0qUKDQ11+sd37dqlBx54QB999JEyMzP1008/KS4uTj/++KO9ZuXKlUpJSdGaNWu0Z88eWa1WxcbGOryHKz4+Xlu2bFFGRoaysrJUUlKiUaNGqbKy0umeAACAe3HqNJaXl5dyc3Nd9vDAc29VX79+vTp06KDs7GzddNNNMsYoLS1NixYtsr93a8OGDQoJCdHmzZs1bdo0FRUVad26ddq4caOGDh0qSdq0aZPCw8O1fft2DRs2zCW9AgCA5snp01iTJk3SunXrGqIXFRUVSZKCgoIkSXl5ecrPz1dcXJy9xsfHRwMGDNDu3bslSdnZ2aqoqHCoCQsLU1RUlL3mXGVlZSouLnZYAACAe3L6AuXy8nL99a9/VWZmpvr27VvtnVgpKSl1asQYo3nz5unGG29UVFSUJCk/P1+SFBIS4lAbEhJif6ZPfn6+vL291bZt22o1v3z/XMnJyVqyZEmd+gQAAM2L02EnNzdX119/vSTpyJEjDuvqc3pr5syZOnjwoLKysqqtO3e7tXkPV001Cxcu1Lx58+yfi4uLFR4eXoeuAQBAU1ersHPw4EFFRUWpRYsW2rlzp8ubmDVrll5//XW99957Dk9mtlqtkn6evfn1BdAFBQX22R6r1ary8nIVFhY6zO4UFBQoJibmvL/n4+PDAxABALhM1OqanT59+ui7776TJF111VU6deqUS37cGKOZM2fqtdde0zvvvKOIiAiH9REREbJarcrMzLSPlZeXa9euXfYgEx0dLS8vL4cam82m3NzcC4YdAGjquvfsraD2ITUuxb+6KxXAhdVqZqdNmzbKy8tThw4d9J///EdVVVUu+fEHHnhAmzdv1v/+7/8qICDAfo1NYGCgfH19ZbFYFB8fr6SkJEVGRioyMlJJSUny8/PThAkT7LVTp07V/PnzFRwcrKCgIC1YsEA9evSw350FAE1J9569ZbPZaqwpPn1av39qR401L88c7Mq2ALdVq7Bz2223acCAAQoNDZXFYlHfvn3l4eFx3tqvvvqq1j++du1aSdLAgQMdxtevX6977rlHkpSQkKDS0lLNmDFDhYWF6tevn7Zt26aAgAB7fWpqqjw9PTV+/HiVlpZqyJAhSk9Pv2CPANCYbDab4pZtrbGGIAO4Tq3CznPPPadx48bpiy++0OzZs3X//fc7hI26MsZctMZisSgxMVGJiYkXrGnZsqVWr16t1atX17snAADgXmp9N9bw4cMl/fxcmzlz5rgk7AAAADQ0p289X79+fUP0AQAA0CCcfoIyAABAc0LYAQAAbo2wAwAA3BphBwAAuDXCDgAAcGuEHQAA4NYIOwAAwK0RdgAAgFsj7AAAALdG2AEAAG6NsAMAANwaYQcAALg1wg4AAHBrhB0AAODWCDsAAMCteTZ2AwDgTrr37C2bzVZjTfHp05eoGwASYQcAXMpmsylu2dYaa16eOfjSNANAEqexAACAmyPsAAAAt8ZpLABQ7a61CQ0N1eGDOZemIQAuQ9gBANXuWptti8Zekl4AuBZhBwBqqfh0iYLah1ykhjutgKaGsAMAtWSqqrjTCmiGuEAZAAC4NcIOAABwa4QdAADg1gg7AADArRF2AACAWyPsAAAAt0bYAQAAbo2wAwAA3BphBwAAuDXCDgAAcGuNGnbee+89jR49WmFhYbJYLNq6davDemOMEhMTFRYWJl9fXw0cOFCHDx92qCkrK9OsWbPUrl07+fv7a8yYMTpx4sQl3AsAANCUNWrY+fHHH9WrVy+tWbPmvOtXrlyplJQUrVmzRnv27JHValVsbKxO/+pFe/Hx8dqyZYsyMjKUlZWlkpISjRo1SpWVlZdqNwA0cd179lZQ+5AaF17gCbivRn0R6IgRIzRixIjzrjPGKC0tTYsWLdK4ceMkSRs2bFBISIg2b96sadOmqaioSOvWrdPGjRs1dOhQSdKmTZsUHh6u7du3a9iwYZdsXwA0XTabjRd4ApexJnvNTl5envLz8xUXF2cf8/Hx0YABA7R7925JUnZ2tioqKhxqwsLCFBUVZa85n7KyMhUXFzssAADAPTXZsJOfny9JCgkJcRgPCQmxr8vPz5e3t7fatm17wZrzSU5OVmBgoH0JDw93cfcAAKCpaLJh5xcWi8XhszGm2ti5LlazcOFCFRUV2Zfjx4+7pFcAAND0NNmwY7VaJanaDE1BQYF9tsdqtaq8vFyFhYUXrDkfHx8ftW7d2mEBAADuqcmGnYiICFmtVmVmZtrHysvLtWvXLsXExEiSoqOj5eXl5VBjs9mUm5trrwEAAJe3Rr0bq6SkRF988YX9c15ennJychQUFKSOHTsqPj5eSUlJioyMVGRkpJKSkuTn56cJEyZIkgIDAzV16lTNnz9fwcHBCgoK0oIFC9SjRw/73VkAAODy1qhhZ+/evRo0aJD987x58yRJkydPVnp6uhISElRaWqoZM2aosLBQ/fr107Zt2xQQEGD/Tmpqqjw9PTV+/HiVlpZqyJAhSk9Pl4eHxyXfHwAA0PQ0atgZOHCgjDEXXG+xWJSYmKjExMQL1rRs2VKrV6/W6tWrG6BDAADQ3DXZa3YAAABcoVFndgCgJt179pbNZqux5kzpWfn5tqyxhldBAJc3wg6AJqu2r3mIS3nrojUALl+cxgIAAG6NsAMAANwaYQcAALg1wg4AAHBrhB0AAODWCDsAAMCtces5gEZRm2fo8HwcAK5A2AHQKGr7DB0AqC9OYwEAALdG2AEAAG6NsAMAANwa1+wAcDkuPgbQlBB2ALgcFx8DaEo4jQUAANwaYQcAALg1TmMBcArX4wBobgg7AJzC9TgAmhvCDgA7Zm0AuCPCDgA7Zm0AuCPCDnCZYNYGwOWKsANcJpi1AXC54tZzAADg1pjZAdwAp6gA4MIIO4Ab4BQVAFwYYQdo4pi1AYD6IewAjai2Qeb3T+2osYZZGwC4MMIO0Ig4/QQADY+7sQAAgFtjZgdoIFxrAwBNA2EHaCCcogKApoGwA9QBszYA0HwQdoA6YNYGAJoPt7lA+emnn1ZERIRatmyp6Ohovf/++43dEgAAaALcYmbnxRdfVHx8vJ5++mndcMMNevbZZzVixAh98skn6tixY2O3h0ukNqeWzpSelZ9vyxprQkNDdfhgjgs7AwA0JrcIOykpKZo6daruu+8+SVJaWprefvttrV27VsnJyY3cXfNUm+DgqlDgqpBS24fvxaW8VWPNK7OHKqh9yEV/CwDQPDT7sFNeXq7s7Gw9/PDDDuNxcXHavXv3eb9TVlamsrIy++eioiJJUnFxscv7+23/G/RNfn6NNWfOlsmvpU+NNSFWqz7+8INL1lNxSYnGrnyjxpqtCWPUNrh9jTW12bfa/NaWBaM0ctmrF62pKP2xxhpjzEVrqiorNeiRv12S36KGGmqouSxqqqoa5G/sL9s0xtRcaJq5kydPGknmgw8+cBhftmyZ6dKly3m/s3jxYiOJhYWFhYWFxQ2W48eP15gVmv3Mzi8sFovDZ2NMtbFfLFy4UPPmzbN/rqqq0vfff6/g4GD7d4qLixUeHq7jx4+rdevWDdf4ZYxj3PA4xg2PY3xpcJwbXnM8xsYYnT59WmFhYTXWNfuw065dO3l4eCj/nNMyBQUFCgk5/3UXPj4+8vFxPLXSpk2b89a2bt262fyP3lxxjBsex7jhcYwvDY5zw2tuxzgwMPCiNc3+1nNvb29FR0crMzPTYTwzM1MxMTGN1BUAAGgqmv3MjiTNmzdPf/jDH9S3b1/1799fzz33nI4dO6bp06c3dmsAAKCRuUXYueOOO3Tq1Ck9/vjjstlsioqK0r/+9S916tSpztv08fHR4sWLq53ugutwjBsex7jhcYwvDY5zw3PnY2wx5mL3awEAADRfzf6aHQAAgJoQdgAAgFsj7AAAALdG2AEAAG6NsHMeTz/9tCIiItSyZUtFR0fr/fffb+yWmq3k5GT95je/UUBAgDp06KCxY8fqs88+c6gxxigxMVFhYWHy9fXVwIEDdfjw4UbquPlLTk6WxWJRfHy8fYxj7BonT57U3XffreDgYPn5+al3797Kzs62r+c4189PP/2kRx55RBEREfL19dVVV12lxx9/XFVVVfYajrFz3nvvPY0ePVphYWGyWCzaunWrw/raHM+ysjLNmjVL7dq1k7+/v8aMGaMTJ05cwr1wgfq+m8rdZGRkGC8vL/P888+bTz75xMyZM8f4+/ubo0ePNnZrzdKwYcPM+vXrTW5ursnJyTEjR440HTt2NCUlJfaa5cuXm4CAAPPqq6+aQ4cOmTvuuMOEhoaa4uLiRuy8efr4449N586dTc+ePc2cOXPs4xzj+vv+++9Np06dzD333GP+/e9/m7y8PLN9+3bzxRdf2Gs4zvWzdOlSExwcbP7xj3+YvLw88/LLL5tWrVqZtLQ0ew3H2Dn/+te/zKJFi8yrr75qJJktW7Y4rK/N8Zw+fbq54oorTGZmptm3b58ZNGiQ6dWrl/npp58u8d7UHWHnHL/97W/N9OnTHca6du1qHn744UbqyL0UFBQYSWbXrl3GGGOqqqqM1Wo1y5cvt9ecPXvWBAYGmmeeeaax2myWTp8+bSIjI01mZqYZMGCAPexwjF3joYceMjfeeOMF13Oc62/kyJFmypQpDmPjxo0zd999tzGGY1xf54ad2hzPH374wXh5eZmMjAx7zcmTJ02LFi3MW2+9dcl6ry9OY/1KeXm5srOzFRcX5zAeFxen3bt3N1JX7qWoqEiSFBQUJEnKy8tTfn6+wzH38fHRgAEDOOZOeuCBBzRy5EgNHTrUYZxj7Bqvv/66+vbtq9tvv10dOnRQnz599Pzzz9vXc5zr78Ybb9SOHTt05MgRSdKBAweUlZWlm2++WRLH2NVqczyzs7NVUVHhUBMWFqaoqKhmdczd4gnKrvLdd9+psrKy2gtEQ0JCqr1oFM4zxmjevHm68cYbFRUVJUn243q+Y3706NFL3mNzlZGRoX379mnPnj3V1nGMXeOrr77S2rVrNW/ePP3pT3/Sxx9/rNmzZ8vHx0eTJk3iOLvAQw89pKKiInXt2lUeHh6qrKzUsmXLdNddd0ni32VXq83xzM/Pl7e3t9q2bVutpjn9XSTsnIfFYnH4bIypNgbnzZw5UwcPHlRWVla1dRzzujt+/LjmzJmjbdu2qWXLlhes4xjXT1VVlfr27aukpCRJUp8+fXT48GGtXbtWkyZNstdxnOvuxRdf1KZNm7R582Z1795dOTk5io+PV1hYmCZPnmyv4xi7Vl2OZ3M75pzG+pV27drJw8OjWlotKCiolnzhnFmzZun111/Xzp07deWVV9rHrVarJHHM6yE7O1sFBQWKjo6Wp6enPD09tWvXLj311FPy9PS0H0eOcf2EhoaqW7duDmPXXXedjh07Jol/l13hwQcf1MMPP6w777xTPXr00B/+8AfNnTtXycnJkjjGrlab42m1WlVeXq7CwsIL1jQHhJ1f8fb2VnR0tDIzMx3GMzMzFRMT00hdNW/GGM2cOVOvvfaa3nnnHUVERDisj4iIkNVqdTjm5eXl2rVrF8e8loYMGaJDhw4pJyfHvvTt21cTJ05UTk6OrrrqKo6xC9xwww3VHptw5MgR+wuH+Xe5/s6cOaMWLRz/LHl4eNhvPecYu1Ztjmd0dLS8vLwcamw2m3Jzc5vXMW+0S6ObqF9uPV+3bp355JNPTHx8vPH39zf/+c9/Gru1ZumPf/yjCQwMNO+++66x2Wz25cyZM/aa5cuXm8DAQPPaa6+ZQ4cOmbvuuotbSevp13djGcMxdoWPP/7YeHp6mmXLlpnPP//c/O1vfzN+fn5m06ZN9hqOc/1MnjzZXHHFFfZbz1977TXTrl07k5CQYK/hGDvn9OnTZv/+/Wb//v1GkklJSTH79++3P06lNsdz+vTp5sorrzTbt283+/btM4MHD+bWc3fw3//936ZTp07G29vbXH/99fbbpOE8Sedd1q9fb6+pqqoyixcvNlar1fj4+JibbrrJHDp0qPGadgPnhh2OsWu88cYbJioqyvj4+JiuXbua5557zmE9x7l+iouLzZw5c0zHjh1Ny5YtzVVXXWUWLVpkysrK7DUcY+fs3LnzvP8Nnjx5sjGmdseztLTUzJw50wQFBRlfX18zatQoc+zYsUbYm7qzGGNM48wpAQAANDyu2QEAAG6NsAMAANwaYQcAALg1wg4AAHBrhB0AAODWCDsAAMCtEXYAAIBbI+wAcGCxWLR169YG/53OnTsrLS2twX/nfNLT09WmTZtG+e2avPvuu7JYLPrhhx8auxXArRB2gMtIQUGBpk2bpo4dO8rHx0dWq1XDhg3Thx9+aK+x2WwaMWJEI3Z5fpcqoKSnp8tisdS4vPvuuw3eBwDX8WzsBgBcOrfddpsqKiq0YcMGXXXVVfrmm2+0Y8cOff/99/aaX96EfLm64447NHz4cPvncePGKSoqSo8//rh9LCgoqDFaA1BHzOwAl4kffvhBWVlZWrFihQYNGqROnTrpt7/9rRYuXKiRI0fa6359Gus///mPLBaLXnrpJf3ud7+Tr6+vfvOb3+jIkSPas2eP+vbtq1atWmn48OH69ttv7dsYOHCg4uPjHX5/7Nixuueeey7YX0pKinr06CF/f3+Fh4drxowZKikpkfTz6Z17771XRUVF9tmVxMREST+/pTkhIUFXXHGF/P391a9fv2ozL+np6erYsaP8/Px066236tSpUxfsw9fXV1ar1b54e3vLz8+v2ti5+vfvr4cffthh7Ntvv5WXl5d27twpSdq0aZP69u2rgIAAWa1WTZgwQQUFBRfsJTExUb1793YYS0tLU+fOnR3G1q9fr+uuu04tW7ZU165d9fTTT19wm8DliLADXCZatWqlVq1aaevWrSorK3Pqu4sXL9Yjjzyiffv2ydPTU3fddZcSEhL0l7/8Re+//76+/PJLPfbYY/Xqr0WLFnrqqaeUm5urDRs26J133lFCQoIkKSYmRmlpaWrdurVsNptsNpsWLFggSbr33nv1wQcfKCMjQwcPHtTtt9+u4cOH6/PPP5ck/fvf/9aUKVM0Y8YM5eTkaNCgQVq6dGm9ej2fiRMn6u9//7t+/brBF198USEhIRowYICkn4PZE088oQMHDmjr1q3Ky8urMQDWxvPPP69FixZp2bJl+vTTT5WUlKRHH31UGzZsqNd2AbfSyC8iBXAJvfLKK6Zt27amZcuWJiYmxixcuNAcOHDAoUaS2bJlizHGmLy8PCPJ/PWvf7Wv//vf/24kmR07dtjHkpOTzbXXXmv/fO5b140x5pZbbrG/adkYYzp16mRSU1Mv2OtLL71kgoOD7Z/Xr19vAgMDHWq++OILY7FYzMmTJx3GhwwZYhYuXGiMMeauu+4yw4cPd1h/xx13VNvWhZxvX86noKDAeHp6mvfee88+1r9/f/Pggw9e8Dsff/yxkWROnz5tjPl/b6guLCw0xhizePFi06tXL4fvpKammk6dOtk/h4eHm82bNzvUPPHEE6Z///4X7Rm4XDCzA1xGbrvtNn399dd6/fXXNWzYML377ru6/vrrlZ6eXuP3evbsaf/nkJAQSVKPHj0cxmo6HVMbO3fuVGxsrK644goFBARo0qRJOnXqlH788ccLfmffvn0yxqhLly72matWrVpp165d+vLLLyVJn376qfr37+/wvXM/u0L79u0VGxurv/3tb5KkvLw8ffjhh5o4caK9Zv/+/brlllvUqVMnBQQEaODAgZKkY8eO1ek3v/32Wx0/flxTp0512P+lS5fa9x8AFygDl52WLVsqNjZWsbGxeuyxx3Tfffdp8eLFNZ5O8fLysv+zxWI571hVVZX9c4sWLRxO50hSRUXFBbd/9OhR3XzzzZo+fbqeeOIJBQUFKSsrS1OnTq3xe1VVVfLw8FB2drY8PDwc1rVq1UqSqvXRkCZOnKg5c+Zo9erV2rx5s7p3765evXpJkn788UfFxcUpLi5OmzZtUvv27XXs2DENGzZM5eXl593exY7jL8f8+eefV79+/Rzqzj0ewOWMsANc5rp16+by5+q0b99eNpvN/rmyslK5ubkaNGjQeev37t2rn376SatWrVKLFj9POL/00ksONd7e3qqsrHQY69OnjyorK1VQUKDf/e535912t27d9NFHHzmMnfvZVcaOHatp06bprbfe0ubNm/WHP/zBvu7//t//q++++07Lly9XeHi4pJ/3uybt27dXfn6+jDH2kJmTk2NfHxISoiuuuEJfffWVwwwSAEeEHeAycerUKd1+++2aMmWKevbsqYCAAO3du1crV67ULbfc4tLfGjx4sObNm6d//vOfuvrqq5Wamlrjg/Kuvvpq/fTTT1q9erVGjx6tDz74QM8884xDTefOnVVSUqIdO3aoV69e8vPzU5cuXTRx4kRNmjRJq1atUp8+ffTdd9/pnXfeUY8ePXTzzTdr9uzZiomJ0cqVKzV27Fht27ZNb731lkv39xf+/v665ZZb9Oijj+rTTz/VhAkT7Os6duwob29vrV69WtOnT1dubq6eeOKJGrc3cOBAffvtt1q5cqV+//vf66233tKbb76p1q1b22sSExM1e/ZstW7dWiNGjFBZWZn27t2rwsJCzZs3r0H2E2h2GvWKIQCXzNmzZ83DDz9srr/+ehMYGGj8/PzMtddeax555BFz5swZe53Oc4Hy/v377evPvYjWmOoXD5eXl5s//vGPJigoyHTo0MEkJydf9ALllJQUExoaanx9fc2wYcPMCy+8UO13pk+fboKDg40ks3jxYvtvPfbYY6Zz587Gy8vLWK1Wc+utt5qDBw/av7du3Tpz5ZVXGl9fXzN69Gjz5z//2eUXKP/in//8p5FkbrrppmrrNm/ebDp37mx8fHxM//79zeuvv+5wfM93bNeuXWvCw8ONv7+/mTRpklm2bJnDBcrGGPO3v/3N9O7d23h7e5u2bduam266ybz22mu17hlwdxZjLuEJbQAAgEuMu7EAAIBbI+wAAAC3RtgBAABujbADAADcGmEHAAC4NcIOAABwa4QdAADg1gg7AADArRF2AACAWyPsAAAAt0bYAQAAbo2wAwAA3Nr/BzyyMsrwhmuIAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "proportion of simulated datasets in which T<=46 = 0.0264\n" ] } ], "source": [ "sns.histplot(simulated_T)\n", "plt.xlabel('Simulated T value')\n", "plt.ylabel('frequency out of ' + str(nReps) + ' tries')\n", "plt.show()\n", "\n", "print('proportion of simulated datasets in which T<=46 = ' + str(np.mean(simulated_T<=46)))" ] }, { "cell_type": "markdown", "id": "5748fa14", "metadata": {}, "source": [ "We obtain a value of $T$ as extreme as the one in our real dataset about 2.5% of the time. \n", "\n", "The p value of the test, which is the same as the proportion of simulated datasets in which T<=46, is about 0.025 or 2.5% (it will vary slightly each time you run the code as the null distribubtion is built from random samples). This is not a bad match for the value we got from the built-in function, which was 0.0266 or 2.66%\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d21ce142", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.13" } }, "nbformat": 4, "nbformat_minor": 5 }