{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# The Current Controversy of the 2013 World Championships\n", "> Some swimmers said that they felt it was easier to swim in one direction versus another in the 2013 World Championships. Some analysts have posited that there was a swirling current in the pool. In this chapter, you'll investigate this claim! References - Quartz Media, Washington Post, SwimSwam, and Cornett, et al. This is the Summary of lecture \"Case Studies in Statistical Thinking\", via datacamp.\n", "\n", "- toc: true \n", "- badges: true\n", "- comments: true\n", "- author: Chanseok Kang\n", "- categories: [Python, Datacamp, Statistics]\n", "- image: images/linreg_swim_lane.png" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "import dc_stat_think as dcst\n", "\n", "plt.rcParams['figure.figsize'] = (10, 5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction to the current controversy\n", "- Task\n", " - Investigate improvement of individual swimmers moving from low- to high-numbered lanes in 50m events\n", " - Compute the size of the effect\n", " - Test the hypothesis that on average there is no difference between low- and high-numbered lanes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### ECDF of improvement from low to high lanes\n", "Now that you have a metric for improvement going from low- to high-numbered lanes, plot an ECDF of this metric. I have put together the swim times of all swimmers who swam a 50 m semifinal in a high numbered lane and the final in a low numbered lane, and vice versa. The swim times are stored in the Numpy arrays `swimtime_high_lanes` and `swimtime_low_lanes`. Entry i in the respective arrays are for the same swimmer in the same event." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "swimtime_high_lanes = np.array([24.62, 22.9 , 27.05, 24.76, 30.31, 24.54, 26.12, 27.71, 23.15,\n", " 23.11, 21.62, 28.02, 24.73, 24.95, 25.83, 30.61, 27.04, 21.67,\n", " 27.16, 30.23, 21.51, 22.97, 28.05, 21.65, 24.54, 26.06])\n", "\n", "swimtime_low_lanes = np.array([24.66, 23.28, 27.2 , 24.95, 32.34, 24.66, 26.17, 27.93, 23.35,\n", " 22.93, 21.93, 28.33, 25.14, 25.19, 26.11, 31.31, 27.44, 21.85,\n", " 27.48, 30.66, 21.74, 23.22, 27.93, 21.42, 24.79, 26.46])" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAE9CAYAAABDUbVaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAXN0lEQVR4nO3df7Cld30X8PeHDSFAK4RkUUyy2cSmLYGh0lxDHOyAlkBgxgSn6ASsYoVmOho7Y4s1DIgQtTPEcehoM8W0ZcAONqVV2x1NjUgDaofQ7C0QSGLKsmTNEpSQrnQChbDuxz/uWTh7f+zd7N3nPOfc+3rN3Lnn+bH3fs537t597/f7PJ+nujsAAMzWU8YuAABgJxLCAABGIIQBAIxACAMAGIEQBgAwAiEMAGAEZ41dwJN1/vnn9969e8cuAwBgU8vLy1/p7t3rHVu4ELZ3797s379/7DIAADZVVYc2OmY5EgBgBEIYAMAIhDAAgBEIYQAAIxDCAABGIIQBAIxACAMAGMFgIayq3ldVX66qz25wvKrqX1bVgaq6t6p+cKhaAADmzZAzYe9Pcs1Jjr86yWWTjxuS/MKAtQAAfNvyoSO59a4DWT50ZLQaBuuY393/rar2nuSU65L8m+7uJHdX1bOr6nnd/aWhagIAWD50JH/9l+7OE0eP5eyznpIPvvmqXHHxuTOvY8xrwi5I8vDU9uHJvjWq6oaq2l9V+x999NGZFAcAbE93H3wsTxw9lmOdfOvosdx98LFR6hgzhNU6+3q9E7v7tu5e6u6l3bvXfQYmAMApuerS83L2WU/JrkqeetZTctWl541Sx5gP8D6c5KKp7QuTPDJSLQDADnHFxefmg2++KncffCxXXXreKEuRybghbF+SG6vq9iQvSfJV14MBALNwxcXnjha+jhsshFXVryZ5eZLzq+pwkn+c5KlJ0t3vTXJHktckOZDk60l+bKhaAADmzZB3R75+k+Od5O8O9f0BAOaZjvkAACMQwgAARiCEAXNtHrpaAwxhzLsjAU5qXrpaAwzBTBgwt+alqzXAEIQwYG7NS1drgCFYjgTm1rx0tQYYghAGzLV56GoNMATLkQAAIxDCAABGIIQBAIxACAMAGIEQBjuUTvQA43J3JOxAOtEDjM9MGOxAOtEDjE8Igx1IJ3qA8VmOhB1IJ3qA8QlhsEPpRA8wLsuRAAAjEMIAAEYghAEAjEAIAwAYgRAGO5SO+QDjcnck7EA65gOMz0wY7EA65gOMTwiDHUjHfIDxWY6EHUjHfIDxCWGwQ+mYDzAuy5EAACMQwgAARiCEAQCMQAgDABiBEAbblI74APPN3ZGwDemIDzD/zITBNqQjPsD8E8JgG9IRH2D+WY6EbUhHfID5J4TBNqUjPsB8sxwJADACIQwAYARCGADACIQwAIARDBrCquqaqnqwqg5U1U3rHN9TVXdV1Ser6t6qes2Q9cB2pkM+wGIZ7O7IqtqV5NYkVyc5nOSeqtrX3fdPnfb2JB/q7l+oqsuT3JFk71A1wXalQz7A4hlyJuzKJAe6+2B3P5Hk9iTXrTqnk/yJyetnJXlkwHpg29IhH2DxDNkn7IIkD09tH07yklXnvDPJf6mqv5fkmUleMWA9sG0d75D/raPHdMgHWBBDhrBaZ1+v2n59kvd397+oqj+f5Feq6oXdfeyEL1R1Q5IbkmTPnj2DFAuLTId8gMUzZAg7nOSiqe0Ls3a58U1JrkmS7v54VZ2T5PwkX54+qbtvS3JbkiwtLa0OckB0yAdYNENeE3ZPksuq6pKqOjvJ9Un2rTrnfyX54SSpqucnOSfJowPWBAAwFwYLYd19NMmNSe5M8kBW7oK8r6purqprJ6f9dJIfr6pPJ/nVJH+ru810AQDb3qAP8O7uO7LSdmJ63zumXt+f5KVD1gAAMI90zAcAGIEQBgtMl3yAxTXociQwHF3yARabmTBYULrkAyw2IQwW1PEu+bsquuQDLCDLkbCgdMkHWGxCGCwwXfIBFpflSACAEQhhAAAjEMIAAEYghAEAjEAIgwWlWz7AYnN3JCwg3fIBFp+ZMFhAuuUDLD4hDBaQbvkAi89yJCwg3fIBFp8QBgtKt3yAxWY5EgBgBEIYAMAIhDAAgBEIYTASzVYBdjYX5sMINFsFwEwYjECzVQCEMBiBZqsAWI6EEWi2CoAQBiPRbBVgZ7McCQAwAiEMAGAEQhgAwAiEMACAEQhhMAO64wOwmrsjYWC64wOwHjNhMDDd8QFYjxAGA9MdH4D1WI6EgemOD8B6hDCYAd3xAVjNciQAwAiEMACAEQhhAAAjEMIAAEYghMHAdMsHYD3ujoQB6ZYPwEYGnQmrqmuq6sGqOlBVN21wzl+rqvur6r6q+rdD1gOzpls+ABsZbCasqnYluTXJ1UkOJ7mnqvZ19/1T51yW5K1JXtrdR6rquUPVA2M43i3/W0eP6ZYPwAmGXI68MsmB7j6YJFV1e5Lrktw/dc6PJ7m1u48kSXd/ecB6YOZ0ywdgI0OGsAuSPDy1fTjJS1ad871JUlW/m2RXknd2938esCaYOd3yAVjPkCGs1tnX63z/y5K8PMmFSf57Vb2wu//vCV+o6oYkNyTJnj17znylAAAzNuSF+YeTXDS1fWGSR9Y557e6+1vd/YUkD2YllJ2gu2/r7qXuXtq9e/dgBQMAzMqQIeyeJJdV1SVVdXaS65PsW3XObyb5i0lSVednZXny4IA1AQDMhcFCWHcfTXJjkjuTPJDkQ919X1XdXFXXTk67M8ljVXV/kruS/IPudg8/ALDtVffqy7Tm29LSUu/fv3/sMuCULR864u5IgB2qqpa7e2m9Yzrmw4B0zAdgI54dCQPSMR+AjQhhMKDjHfN3VXTMB+AEliNhQDrmA7CRk4awqjprcpcjcJp0zAdgPZstR/7e8RdV9a8GrgUAYMfYLIRNP3ropUMWAgCwk2wWwhariRgAwILY7ML876+qe7MyI/ZnJq8z2e7uftGg1QEAbFObhbDnz6QK2KZ0ywdgIycNYd19KEmq6tlJLpvs/oPu/urQhcGi0y0fgJM56TVhVXV2Vb0/yUNJbkvyi0keqqr3VdXZw5cHi0u3fABOZrML89+e5KlJLuruF3f3n02yJyszaP9o6OJgkemWD8DJVPfGN0BW1WeTXNndX1+1/7uS3N3dLxy4vjWWlpZ6//79s/62cFpcEwaws1XVcncvrXdsswvzj60OYEnS3Y9XlfYVsAnd8gHYyGYhrKvq3JzYtPW4YwPUAwCwI2wWwp6VZDnrhzAzYQAAp2mzFhV7Z1QHAMCOslmLildV1evW2f+Gqrp6uLIAALa3zVpUvCvJx9bZ/ztJbj7z5cBiWz50JLfedSDLh46MXQoAc26za8Ke0d2Prt7Z3f+7qp45UE2wkHTIB+DJ2Gwm7JyqWhPUquqpSZ4+TEmwmHTIB+DJ2CyE/fskvzg96zV5/d7JMWBCh3wAnozNliPfnuSfJjlUVYcm+/Yk+eV4bBGc4IqLz80H33yVDvkAnJKTPrbo2ydVPT3J90w2D3T3Hw9a1Ul4bBEAsChO9tiizVpU/EySTELX93f3Z44HsKr62TNeKQDADrHZNWHXT71+66pj15zhWgAAdozNQlht8Hq9bQAATtFmIaw3eL3eNgAAp2izuyN/oKr+KCuzXk+fvM5k+5xBK4MFsHzoiLshATgtmz3Ae9esCoFFo0M+AFux2XIksAEd8gHYCiEMTpMO+QBsxWbXhAEb0CEfgK0QwmALrrj4XOELgNNiORIAYARCGADACIQwAIARCGEAACMQwuAULB86klvvOpDlQ0fGLgWAbcLdkbAJnfEBGIKZMNiEzvgADEEIg03ojA/AEAYNYVV1TVU9WFUHquqmk5z3uqrqqloash44Hcc74//UK7/PUiQAZ8xg14RV1a4ktya5OsnhJPdU1b7uvn/Ved+d5CeTfGKoWmCrdMYH4EwbcibsyiQHuvtgdz+R5PYk161z3j9JckuSbwxYCwDAXBkyhF2Q5OGp7cOTfd9WVS9OclF3/8cB6wAAmDtDhrBaZ19/+2DVU5K8J8lPb/qFqm6oqv1Vtf/RRx89gyUCAIxjyBB2OMlFU9sXJnlkavu7k7wwyUer6qEkVyXZt97F+d19W3cvdffS7t27BywZvkODVgCGNGSz1nuSXFZVlyT5YpLrk7zh+MHu/mqS849vV9VHk7ylu/cPWBOcEg1aARjaYDNh3X00yY1J7kzyQJIPdfd9VXVzVV071PeFM0GDVgCGNuhji7r7jiR3rNr3jg3OffmQtcCTcbxB67eOHtOgFYBBeHYkrON4g9a7Dz6Wqy49z1IkAGecEAYb0KAVgCF5diQAwAiEMACAEQhhAAAjEMIAAEYghMEqOuUDMAvujoQpOuUDMCtmwmCKTvkAzIoQBlOOd8rfVdEpH4BBWY6EKTrlAzArQhisolM+ALNgORIAYARCGADACIQwAIARCGEAACMQwmBCp3wAZsndkRCd8gGYPTNhEJ3yAZg9IQyiUz4As2c5EqJTPgCzJ4TBhE75AMyS5UgAgBEIYQAAIxDCAABGIIQBAIxACGPb0fkegEXg7ki2FZ3vAVgUZsLYVnS+B2BRCGFsKzrfA7AoLEeyreh8D8CiEMLYdnS+B2ARWI4EABiBEAYAMAIhDABgBEIYAMAIhDAWmu74ACwqd0eysHTHB2CRmQljYemOD8AiE8JYWLrjA7DILEeysHTHB2CRCWEsNN3xAVhUgy5HVtU1VfVgVR2oqpvWOf5TVXV/Vd1bVR+pqouHrAcAYF4MFsKqaleSW5O8OsnlSV5fVZevOu2TSZa6+0VJfiPJLUPVAwAwT4acCbsyyYHuPtjdTyS5Pcl10yd0913d/fXJ5t1JLhywHgCAuTFkCLsgycNT24cn+zbypiS/PWA9AABzY8gL82udfb3uiVU/mmQpycs2OH5DkhuSZM+ePWeqPkawfOiIuxkBIMOGsMNJLpravjDJI6tPqqpXJHlbkpd19zfX+0LdfVuS25JkaWlp3SDH/NPhHgC+Y8jlyHuSXFZVl1TV2UmuT7Jv+oSqenGSf53k2u7+8oC1MAd0uAeA7xgshHX30SQ3JrkzyQNJPtTd91XVzVV17eS0f57ku5L8elV9qqr2bfDl2AZ0uAeA76juxVrdW1pa6v37949dBqfJNWEA7CRVtdzdS+sd0zGfmdLhHgBWeIA3AMAIhDAAgBEIYQAAIxDCAABGIIQxE8uHjuTWuw5k+dCRsUsBgLng7kgGp1M+AKxlJozB6ZQPAGsJYQxOp3wAWMtyJIO74uJz88E3X6VTPgBMEcKYCZ3yAeBEliMBAEYghAEAjEAIAwAYgRAGADACIWwb0I0eABaPuyMXnG70ALCYzIQtON3oAWAxCWELTjd6AFhMliMXnG70ALCYhLBtQDd6AFg8liMBAEYghAEAjEAIAwAYgRAGADACIWxgutkDAOtxd+SAdLMHADZiJmxAutkDABsRwgakmz0AsBHLkQPSzR4A2IgQNjDd7AGA9ViOBAAYgRAGADACIQwAYARC2CqaqwIAs+DC/CmaqwIAs2ImbIrmqgDArAhhUzRXBQBmxXLkFM1VAYBZEcJW0VwVAJgFy5EAACMQwgAARiCEAQCMYNAQVlXXVNWDVXWgqm5a5/jTqurXJsc/UVV7h6wHAGBeDBbCqmpXkluTvDrJ5UleX1WXrzrtTUmOdPf3JHlPkncPVQ8AwDwZcibsyiQHuvtgdz+R5PYk160657okH5i8/o0kP1xVNWBNAABzYcgQdkGSh6e2D0/2rXtOdx9N8tUkazqkVtUNVbW/qvY/+uijA5ULADA7Q4aw9Wa0+jTOSXff1t1L3b20e/fuM1IcAMCYhgxhh5NcNLV9YZJHNjqnqs5K8qwkfzhgTQAAc2HIjvn3JLmsqi5J8sUk1yd5w6pz9iV5Y5KPJ3ldkt/p7jUzYdOWl5e/UlWHBqj3VJyf5Csjfe95ZUzWMiZrGZO1jMlaxmQtY7LWoo3JxRsdGCyEdffRqroxyZ1JdiV5X3ffV1U3J9nf3fuS/HKSX6mqA1mZAbv+FL7uaOuRVbW/u5fG+v7zyJisZUzWMiZrGZO1jMlaxmSt7TQmgz47srvvSHLHqn3vmHr9jSR/dcgaAADmkY75AAAjEMKenNvGLmAOGZO1jMlaxmQtY7KWMVnLmKy1bcakNrkOHgCAAZgJAwAYgRA2paqeU1UfrqrPTT6fu8F5b5yc87mqeuPU/n9WVQ9X1eOzq3oYW3n4elW9dbL/wap61SzrHtLpjklVnVdVd1XV41X187Oue0hbGJOrq2q5qj4z+fyXZl37ULYwJldW1acmH5+uqr8y69qHspXfJ5PjeyZ/f94yq5qHtoWfk71V9cdTPyvvnXXtQ9nivzsvqqqPV9V9k98r58yy9tPW3T4mH0luSXLT5PVNSd69zjnPSXJw8vncyetzJ8euSvK8JI+P/V62OA67knw+yaVJzk7y6SSXrzrn7yR57+T19Ul+bfL68sn5T0tyyeTr7Br7PY08Js9M8heS/ESSnx/7vczJmLw4yZ+evH5hki+O/X7mYEyekeSsyevnJfny8e1F/tjKmEwd/3dJfj3JW8Z+P2OPSZK9ST479nuYszE5K8m9SX5gsn3eovy7YybsRNMPFP9Akteuc86rkny4u/+wu48k+XCSa5Kku+/u7i/NpNJhbeXh69club27v9ndX0hyYPL1Ft1pj0l3f627/0eSb8yu3JnYyph8sruPP0HjviTnVNXTZlL1sLYyJl/vlWfoJsk5WecRbgtqK79PUlWvzcp/du+bUb2zsKUx2aa2MiavTHJvd386Sbr7se7+fzOqe0uEsBP9yeMhavL5ueuccyoPJl90W3n4+nYdnzP2QPpt5EyNyY8k+WR3f3OgOmdpS2NSVS+pqvuSfCbJT0yFskV22mNSVc9M8g+TvGsGdc7SVv/uXFJVn6yqj1XVDw1d7IxsZUy+N0lX1Z1V9ftV9TMzqPeMGLRZ6zyqqv+a5E+tc+htp/ol1tm3Xf7HetxWHr6+XcfnjD2QfhvZ8phU1QuSvDsr/5PdDrY0Jt39iSQvqKrnJ/lAVf12rzS1XmRbGZN3JXlPdz++zSaBtjImX0qyp7sfq6orkvxmVb2gu//oTBc5Y1sZk7OycsnHn0vy9SQfqarl7v7ImS3xzNtxIay7X7HRsar6P1X1vO7+UlUdvyZjtcNJXj61fWGSj57RIsf3ZB6+frhOfPj6qfzZRbSVMdmutjQmVXVhkv+Q5G929+eHL3cmzsjPSXc/UFVfy8r1cvuHK3cmtjImL0nyuqq6Jcmzkxyrqm9096Lf4HLaY9IrFz19M0m6e7mqPp+VmaCd/HNyOMnHuvsrSVJVdyT5wSRzH8IsR57o+APFM/n8W+ucc2eSV1bVubVy9+QrJ/u2k28/fL2qzs7KBZD7Vp0zPVbTD1/fl+T6yV0slyS5LMnvzajuIW1lTLar0x6Tqnp2kv+U5K3d/bszq3h4WxmTSyb/sKSqLk7yfUkemk3ZgzrtMenuH+ruvd29N8nPJfnZbRDAkq39nOyuql1JUlWXZuV37MEZ1T2krfyOvTPJi6rqGZO/Qy9Lcv+M6t6ase8MmKePrKwtfyTJ5yafnzPZv5Tkl6bO+9tZueD8QJIfm9p/S1YS+bHJ53eO/Z62MBavSfIHWblb5W2TfTcnuXby+pys3K10ICsh69KpP/u2yZ97MMmrx34vczImD2Xlf2yPT342Lp91/fM0JknenuRrST419fHcsd/PyGPyN7Jy8fmnkvx+kteO/V7GHpNVX+Od2SZ3R27x5+RHJj8nn578nPzlsd/L2GMyOfajk3H5bJJbxn4vp/qhYz4AwAgsRwIAjEAIAwAYgRAGADACIQwAYARCGADACIQwgCRV9ZNV9UBVfXDsWoCdQYsKgCRV9T+z0tfuC2PXAuwMZsKAHa+q3pvk0iT7qurvj10PsDOYCQNIUlUPJVnqyfPnAIZmJgwAYARCGADACIQwAIARCGEAACNwYT4AwAjMhAEAjEAIAwAYgRAGADACIQwAYARCGADACIQwAIARCGEAACMQwgAARvD/AVt6jn4isBGzAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Compute the fractional improvement of being in high lane: f\n", "f = (swimtime_low_lanes - swimtime_high_lanes) / swimtime_low_lanes\n", "\n", "# Make x and y values for ECDF: x, y\n", "x, y = dcst.ecdf(f)\n", "\n", "# Plot the ECDFs as dots\n", "_ = plt.plot(x, y, marker='.', linestyle='none')\n", "\n", "# Label the axes and show the plot\n", "_ = plt.xlabel('f')\n", "_ = plt.ylabel('ECDF')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Oooo, this is starting to paint a picture of lane bias. The ECDF demonstrates that all but three of the 26 swimmers swam faster in the high numbered lanes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Estimation of mean improvement\n", "You will now estimate how big this current effect is. Compute the mean fractional improvement for being in a high-numbered lane versus a low-numbered lane, along with a 95% confidence interval of the mean.\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "mean frac. diff.: 0.01051\n", "95% conf int of mean frac. diff.: [0.00630, 0.01596]\n" ] } ], "source": [ "# Compute the mean difference: f_mean\n", "f_mean = np.mean(f)\n", "\n", "# Draw 10,000 bootstrap replicates: bs_reps\n", "bs_reps = dcst.draw_bs_reps(f, np.mean, size=10000)\n", "\n", "# Compute 95% confidence interval: conf_int\n", "conf_int = np.percentile(bs_reps, [2.5, 97.5])\n", "\n", "# Print the result\n", "print(\"\"\"\n", "mean frac. diff.: {0:.5f}\n", "95% conf int of mean frac. diff.: [{1:.5f}, {2:.5f}]\"\"\".format(f_mean, *conf_int))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "it sure looks like swimmers are faster in lanes 6-8." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### How should we test the hypothesis?\n", "Q: You are interested in the presence of lane bias toward higher lanes, presumably due to a slight current in the pool. A natural null hypothesis to test, then, is that the mean fractional improvement going from low to high lane numbers is zero. Which of the following is a good way to simulate this null hypothesis?\n", "\n", "A: Subtract the mean of `f` from `f` to generate `f_shift`. Then, take bootstrap replicate of the mean from this `f_shift`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Hypothesis test: Does lane assignment affect performance?\n", "Perform a bootstrap hypothesis test of the null hypothesis that the mean fractional improvement going from low-numbered lanes to high-numbered lanes is zero. Take the fractional improvement as your test statistic, and \"at least as extreme as\" to mean that the test statistic under the null hypothesis is greater than or equal to what was observed." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p = 0.00028\n" ] } ], "source": [ "# Shift f: f_shift\n", "f_shift = f - f_mean\n", "\n", "# Draw 100,000 bootstrap replicates of the mean: bs_reps\n", "bs_reps = dcst.draw_bs_reps(f_shift, np.mean, 100000)\n", "\n", "# Compute and report the p-value\n", "p_val = np.sum(bs_reps >= f_mean) / 100000\n", "print('p =', p_val)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A p-value of 0.0003 is quite small and suggests that the mean fractional improvment is greater than zero. For fun, I tested the more restrictive hypothesis that lane number has no bearing at all on performance (item (1) in the previous MCQ), and I got an even smaller p-value of about 0.00001. You can perform that test, too, for practice if you like." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Did the 2015 event have this problem?\n", "You would like to know if this is a typical problem with pools in competitive swimming. To address this question, perform a similar analysis for the results of the 2015 FINA World Championships. That is, compute the mean fractional improvement for going from lanes 1-3 to lanes 6-8 for the 2015 competition, along with a 95% confidence interval on the mean. Also test the hypothesis that the mean fractional improvement is zero." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "swimtime_high_lanes_15 = np.array([27.7 , 24.64, 23.21, 23.09, 26.87, 30.74, 21.88, 24.5 , 21.86,\n", " 25.9 , 26.2 , 24.73, 30.13, 26.92, 24.31, 30.25, 26.76])\n", "swimtime_low_lanes_15 = np.array([27.66, 24.69, 23.29, 23.05, 26.87, 31.03, 22.04, 24.51, 21.86,\n", " 25.64, 25.91, 24.77, 30.14, 27.23, 24.31, 30.2 , 26.86])" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "mean frac. diff.: 0.00079\n", "95% conf int of mean frac. diff.: [-0.00191, 0.00336]\n", "p-value: 0.28205\n" ] } ], "source": [ "# Compute f and its mean\n", "f = (swimtime_low_lanes_15 - swimtime_high_lanes_15) / swimtime_low_lanes_15\n", "f_mean = np.mean(f)\n", "\n", "# Draw 10,000 bootstrap replicates\n", "bs_reps = dcst.draw_bs_reps(f, np.mean, 10000)\n", "\n", "# Compute 95% confidence interval\n", "conf_int = np.percentile(bs_reps, [2.5, 97.5])\n", "\n", "# Shift f\n", "f_shift = f - f_mean\n", "\n", "# Draw 100,000 bootstrap replicates of the mean\n", "bs_reps = dcst.draw_bs_reps(f_shift, np.mean, 100000)\n", "\n", "# Compute the p-value\n", "p_val = np.sum(bs_reps >= f_mean) / 100000\n", "\n", "# Print the results\n", "print(\"\"\"\n", "mean frac. diff.: {0:.5f}\n", "95% conf int of mean frac. diff.: [{1:.5f}, {2:.5f}]\n", "p-value: {3:.5f}\"\"\".format(f_mean, *conf_int, p_val))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Both the confidence interval an the p-value suggest that there was no lane bias in 2015." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The zigzag effect\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Which splits should we consider?\n", "Q: As you proceed to quantitatively analyze the zigzag effect in the 1500 m, which splits should you include in our analysis?\n", "\n", "A: \n", "You should include all splits except the first two and the last two. You should neglect the last two because swimmers stop pacing themselves and \"kick\" for the final stretch. The first two are different because they involve jumping off the starting blocks and more underwater swimming than others." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### EDA: mean differences between odd and even splits\n", "To investigate the differences between odd and even splits, you first need to define a difference metric. In previous exercises, you investigated the improvement of moving from a low-numbered lane to a high-numbered lane, defining f = (ta - tb) / ta. There, the ta in the denominator served as our reference time for improvement. Here, you are considering both improvement and decline in performance depending on the direction of swimming, so you want the reference to be an average. So, we will define the fractional difference as $f = 2 \\frac{(t_a - t_b)}{(t_a + t_b)}$.\n", "\n", "Your task here is to plot the mean fractional difference between odd and even splits versus lane number." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "f_13 = np.array([-0.01562214, -0.0146381 , -0.00977673, -0.00525713, 0.00204104,\n", " 0.00381014, 0.0075664 , 0.01525869])\n", "f_15 = np.array([-0.00516018, -0.00392952, -0.00099284, 0.00059953, -0.002424 ,\n", " -0.00451099, 0.00047467, 0.00081962])\n", "lanes = np.array([1, 2, 3, 4, 5, 6, 7, 8])" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAE9CAYAAACGIy/LAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df5xddX3n8deHJCajEiKYYH6AiTWNEaEERtBifYhpCHaRRFBKRTfdssu2W1e7rtlH8uhWXNxiNFKRLWs3CjbFVgWahlhbIgTsrlSEiQkEAnmQrbSZSSSREH7toCF+9o97BifhzsydH3fOzZnX8/G4jzvne7733M+9anzf8z3f843MRJIkSUe/Y8ouQJIkSSPDYCdJklQRBjtJkqSKMNhJkiRVhMFOkiSpIgx2kiRJFTG+7AJawWtf+9qcPXt22WVIkiQNaPPmzT/JzKn19hnsgNmzZ9PR0VF2GZIkSQOKiH/ua59DsZIkSRVhsJMkSaoIg50kSVJFeI1dHw4ePEhnZycvvPBC2aWMikmTJjFr1iwmTJhQdimSJGmIDHZ96Ozs5Nhjj2X27NlERNnlNFVm8uSTT9LZ2cmcOXPKLkeSJA2RQ7F9eOGFFzjhhBMqH+oAIoITTjhhzJydlCSpqgx2/RgLoa7HWPqskiRVlcGuhe3atYtzzz2X+fPnc8opp/DFL34RgP3797No0SLmzp3LokWLeOqppwB49NFHefvb387EiRP5/Oc//9JxXnjhBc466yx+5Vd+hVNOOYUrr7yylM8jSVJVrd/SxTmr7mLOim9zzqq7WL+lq5Q6DHYtbPz48VxzzTU88sgj3HvvvVx//fVs376dVatWsXDhQh577DEWLlzIqlWrADj++OO57rrr+MQnPnHYcSZOnMhdd93FAw88wNatW7n99tu59957y/hIkiRVzvotXaxct42uA90k0HWgm5XrtpUS7koNdhFxfkTsiIidEbGizv6JEfHNYv8PImJ20X5CRNwdEc9FxJ8e8ZrvFsfcWjymjcZnaUZSnz59OmeccQYAxx57LPPnz6erq4vbbruNZcuWAbBs2TLWr18PwLRp03jrW9/6spmtEcGrX/1qoDbb9+DBgw69SpI0QlZv3EH3wUOHtXUfPMTqjTtGvZbSgl1EjAOuB94DvBn4rYh48xHdLgeeysw3Al8APlu0vwD8EfAJ6rssM08vHntHvvrDjUZSf/zxx9myZQtnn302TzzxBNOnTwdq4W/v3oE/4qFDhzj99NOZNm0aixYt4uyzzx6x2iRJGst2H+geVHszlXnG7ixgZ2b+U2b+DPgGsOSIPkuAtcXftwILIyIy8/nM/B61gFe6Zif15557josvvphrr72WyZMnD+kY48aNY+vWrXR2dnLffffx0EMPjUhtkiSNdTOmtA2qvZnKDHYzgV29tjuLtrp9MvNF4GnghAaO/dViGPaPoo8xx4i4IiI6IqJj3759g6++l2Ym9YMHD3LxxRdz2WWXcdFFFwFw4oknsmfPHgD27NnDtGmNjzZPmTKFd73rXdx+++3Drk2SJMHyxfNomzDusLa2CeNYvnjeqNdSZrCrF7hyCH2OdFlmngr8WvH4cL1OmbkmM9szs33q1KkDFtufZiX1zOTyyy9n/vz5fPzjH3+p/cILL2Tt2tqJzLVr17JkyZEnOg+3b98+Dhw4AEB3dzd33nknb3rTm4ZVmyRJqlm6YCafuehUZk5pI4CZU9r4zEWnsnTBkeermq/MlSc6gZN6bc8CdvfRpzMixgPHAfv7O2hmdhXPz0bEX1Eb8v2LkSq6nuWL57Fy3bbDhmNHIqnfc8893HTTTZx66qmcfvrpAFx99dWsWLGCSy65hBtuuIGTTz6ZW265BYAf//jHtLe388wzz3DMMcdw7bXXsn37dvbs2cOyZcs4dOgQP//5z7nkkku44IILhlWbJEn6haULZpYS5I5UZrC7H5gbEXOALuBS4INH9NkALAO+D7wfuCsz+zxjV4S/KZn5k4iYAFwA3NmM4nvr+Q9y9cYd7D7QzYwpbSxfPG/Y/wG/4x3voK+Pu2nTppe1ve51r6Ozs/Nl7aeddhpbtmwZVi2SJKn1lRbsMvPFiPgIsBEYB9yYmQ9HxFVAR2ZuAG4AboqIndTO1F3a8/qIeByYDLwiIpYC5wH/DGwsQt04aqHuy6PxeVolqUuSpLGrzDN2ZObfAX93RNsne/39AvCBPl47u4/DnjlS9UmSJB1NXHlCkiSpIgx2kiRJFWGwkyRJqgiDnSRJUkUY7FrYrl27OPfcc5k/fz6nnHIKX/ziFwHYv38/ixYtYu7cuSxatIinnnoKgEcffZS3v/3tTJw4kc9//vOHHWv27Nkv3Q+vvb191D+LJElqPoNdCxs/fjzXXHMNjzzyCPfeey/XX38927dvZ9WqVSxcuJDHHnuMhQsXsmrVKgCOP/54rrvuOj7xiU/UPd7dd9/N1q1b6ejoGM2PIUmSRonBbqQ8eDN84S3wqSm15wdvHvYhp0+fzhlnnAHAsccey/z58+nq6uK2225j2bJlACxbtoz169cDMG3aNN761rcyYcKEYb+3JEk6+hjsRsKDN8O3PgpP7wKy9vytj45IuOvx+OOPs2XLFs4++2yeeOIJpk+fDtTC3969ewd8fURw3nnnceaZZ7JmzZoRq0uSJLWOUm9QXBmbroKD3Ye3HeyutZ92ybAP/9xzz3HxxRdz7bXXMnny5CEd45577mHGjBns3buXRYsW8aY3vYl3vvOdw65NkiS1Ds/YjYSnX74+a7/tg3Dw4EEuvvhiLrvsMi666CIATjzxRPbs2QPAnj17mDZt2oDHmTFjBlAbrn3f+97HfffdN+zaJElSazHYjYTjZg2uvUGZyeWXX878+fP5+Mc//lL7hRdeyNq1awFYu3YtS5Ys6fc4zz//PM8+++xLf3/nO9/hLW95y7BqkyRJrceh2JGw8JO1a+p6D8dOaKu1D8M999zDTTfd9NJtSgCuvvpqVqxYwSWXXMINN9zAySefzC233ALAj3/8Y9rb23nmmWc45phjuPbaa9m+fTs/+clPeN/73gfAiy++yAc/+EHOP//8YdUmSZJaj8FuJPRcR7fpqtrw63GzaqFumNfXveMd7yAz6+7btGnTy9pe97rX0dn58uHfyZMn88ADDwyrFkmS1PoMdiPltEtGZKKEJEnSUHmNnSRJUkUY7CRJkirCYNePvq5vq6Kx9FklSaoqg10fJk2axJNPPjkmAk9m8uSTTzJp0qSyS5EkScPg5Ik+zJo1i87OTvbt21d2KaNi0qRJzJo1vPvuSZKkchns+jBhwgTmzJlTdhmSJEkNcyhWkiSpIgx2kiRJFWGwkyRJqgiDnSRJUkUY7CRJkirCYCdJklQRBjtJkqSKMNhJkiRVhMFOkiSpIgx2kiRJFWGwkyRJqohSg11EnB8ROyJiZ0SsqLN/YkR8s9j/g4iYXbSfEBF3R8RzEfGnR7zmzIjYVrzmuoiI0fk0kiRJ5Sot2EXEOOB64D3Am4Hfiog3H9HtcuCpzHwj8AXgs0X7C8AfAZ+oc+gvAVcAc4vH+SNfvSRJUusp84zdWcDOzPynzPwZ8A1gyRF9lgBri79vBRZGRGTm85n5PWoB7yURMR2YnJnfz8wE/gJY2tRPIUmS1CLKDHYzgV29tjuLtrp9MvNF4GnghAGO2TnAMQGIiCsioiMiOvbt2zfI0iVJklpPmcGu3rVvOYQ+Q+qfmWsysz0z26dOndrPISVJko4OZQa7TuCkXtuzgN199YmI8cBxwP4BjjlrgGNKkiRVUpnB7n5gbkTMiYhXAJcCG47oswFYVvz9fuCu4tq5ujJzD/BsRLytmA37r4HbRr50SZKk1jO+rDfOzBcj4iPARmAccGNmPhwRVwEdmbkBuAG4KSJ2UjtTd2nP6yPicWAy8IqIWAqcl5nbgd8D/hxoA/6+eEiSJFVe9HMCbMxob2/Pjo6OssuQJEkaUERszsz2evtceUKSJKkiDHaSJEkVYbCTJEmqCIOdJElSRRjsJEmSKsJgJ0mSVBEGO0mSpIow2EmSJFWEwU6SJKkiDHaSJEkVYbCTJEmqiPFlFyBJklrb+i1drN64g90HupkxpY3li+exdMHMsstSHQY7SZLUp/Vbuli5bhvdBw8B0HWgm5XrtgEY7lqQQ7GSJKlPqzfueCnU9eg+eIjVG3eUVJH6Y7CTJEl92n2ge1DtKpfBTpIk9WnGlLZBtatcBjtJktSn5Yvn0TZh3GFtbRPGsXzxvJIqUn+cPCFJkvrUM0HCWbFHB4OdJEnq19IFMw1yRwmHYiVJkirCYCdJklQRBjtJkqSKMNhJkiRVxICTJyJiEnAB8GvADKAbeAj4dmY+3NzyJEmS1Kh+g11EfAp4L/Bd4AfAXmAS8MvAqiL0/efMfLC5ZUqSJGkgA52xuz8zP9XHvj+JiGnAySNbkiRJkoai32CXmd8eYP9eamfxJEmSVLKGblAcEb8MLAde3/s1mfnuJtUlSZKkQWp05YlbgD8Dvgwcal45kiRJGqpGg92LmfmlplYiSZKkYWn0Pnbfioj/EBHTI+L4nkdTK5MkSdKgNBrsllG7xu4fgc3Fo2O4bx4R50fEjojYGREr6uyfGBHfLPb/ICJm99q3smjfERGLe7U/HhHbImJrRAy7RkmSpKNFQ0OxmTlnpN84IsYB1wOLgE7g/ojYkJnbe3W7HHgqM98YEZcCnwV+MyLeDFwKnELtpsl3RsQvZ2bP9X/nZuZPRrpmSZKkVtbQGbuIeGVE/NeIWFNsz42IC4b53mcBOzPznzLzZ8A3gCVH9FkCrC3+vhVYGBFRtH8jM3+amT8CdhbHkyRJGrMaHYr9KvAz4FeL7U7gvw/zvWcCu3ptdxZtdftk5ovA08AJA7w2ge9ExOaIuKKvN4+IKyKiIyI69u3bN6wPIkmS1AoaDXa/lJmfAw4CZGY3EMN873qvzwb79PfaczLzDOA9wO9HxDvrvXlmrsnM9sxsnzp1aqM1S5IktaxGg93PIqKNIjxFxC8BPx3me3cCJ/XangXs7qtPRIwHjgP29/fazOx53gv8DQ7RSpKkMaLRYPcp4HbgpIj4S2AT8F+G+d73A3MjYk5EvILaZIgNR/TZQG1GLsD7gbsyM4v2S4tZs3OAucB9EfGqiDgWICJeBZwHPDTMOiVJko4Kjc6K/U5EbAbeRm0Y9GPDnXWamS9GxEeAjcA44MbMfDgirgI6MnMDcANwU0TspHam7tLitQ9HxM3AduBF4Pcz81BEnAj8TW1+BeOBv8rM24dTpyRJ0tEiaifABugUsQH4OrAhM59velWjrL29PTs6vOWdJElqfRGxOTPb6+1rdCj2GuDXgO0RcUtEvD8iJo1YhZIkSRq2Rodi/wH4h+Kmwu8G/h1wIzC5ibVJkjQq1m/pYvXGHew+0M2MKW0sXzyPpQuOvAOX1PoaCnYAxazY9wK/CZzBL24cLEnSUWv9li5WrttG98Ha4kVdB7pZuW4bgOFOR51GV574JvAItbN111O7r91/bGZhkiSNhtUbd7wU6np0HzzE6o07SqpIGrpGz9h9Ffhgr7VYJUmqhN0HugfVLrWyRidP/G9g5QivFStJUulmTGkbVLvUyspcK1aSpNItXzyPtgnjDmtrmzCO5YvnlVSRNHRlrhUrSVLpli6YyWcuOpWZU9oIYOaUNj5z0alOnNBRqdFr7JqxVqwkSS1h6YKZBjlVQqPB7koOXyv2HOC3m1WUJEmSBq/RGxTfERE/ZATXipUkSdLIavgGxZn5JPDtJtYiSZKkYWh08oQkSZJanMFOkiSpIgYd7CLiimYUIkmSpOFp+Bq7Xn4XWDPShUiSmmP9li5Wb9zB7gPdzJjSxvLF87y1h1RRQwl23phYko4S67d0sXLdtpcWue860M3KddsADHdSBQ3lGrv3jngVkqSmWL1xx0uhrkf3wUOs3rijpIokNdOgg11mdjajEEnSyNt9oHtQ7ZKObs6KlaQKmzGlbVDtko5uBjtJqrDli+fRNmHcYW1tE8axfPG8kiqS1Ez9Tp6IiIv625+Z60a2HEnSSOqZIOGsWGlsGGhWbM9EiWnArwJ3FdvnAt8FDHaS1OKWLphpkJPGiH6DXWb+G4CI+FvgzZm5p9ieDlzf/PIkSZLUqEavsZvdE+oKTwC/3IR6JEmSNESN3qD4uxGxEfg6kMClwN1Nq0qSJOlo8uDNsOkqeLoTjpsFCz8Jp10y6mU0FOwy8yPFRIpfK5rWZObfNK8sSZKko8SDN8O3PgoHi/tDPr2rtg2jHu4aXlKsmAHrZAlJradFfilLGqM2XfWLUNfjYHetvZWCXUQ8S23ota7MnDziFUnSYLTQL2Wpsvzx1L+n+1iUq6/2JhpoVuyxABFxFfBj4CYggMuAY5tenSQNpIV+KUuV5I+ngR03q/a91GsfZY3Oil2cmf8zM5/NzGcy80vAxc0sTJIa0kK/lKVK6u/Hk2oWfhImHLFM34S2WvsoazTYHYqIyyJiXEQcExGXAYeG++YRcX5E7IiInRGxos7+iRHxzWL/DyJidq99K4v2HRGxuNFjSqqYvn4Rl/BLWaokfzwN7LRL4L3XwXEnAVF7fu91rTsrFvgg8MXiAfC9om3IImIctZscLwI6gfsjYkNmbu/V7XLgqcx8Y0RcCnwW+M2IeDO1W66cAswA7oyInvvqDXRMSVWy8JOHDxNBab+UpUpqoWHGlnbaJS0xNN3QGbvMfDwzl2Tma4vH0sx8fJjvfRawMzP/KTN/BnwDWHJEnyXA2uLvW4GFERFF+zcy86eZ+SNgZ3G8Ro4pqUpa6JeyVEktNMyogTV0xi4iZgH/AziH2izZ7wEfy8zhnIedCfT+CdAJnN1Xn8x8MSKeBk4o2u894rU9CyEOdEwAIuIK4AqAk08+eWifQFJraJFfylIl9fxvy1mxR4VGh2K/CvwV8IFi+0NF26JhvHfUaTvy1ip99emrvd4ZyLq3a8nMNcAagPb29j5v6SJJ0pjnj6ejRqOTJ6Zm5lcz88Xi8efA1GG+dydwUq/tWcDuvvpExHjgOGB/P69t5JiSJEmV1Giw+0lEfKiYFTsuIj4EPDnM974fmBsRcyLiFdQmQ2w4os8GYFnx9/uBuzIzi/ZLi1mzc4C5wH0NHlOSJKmSGh2K/R3gT4EvUBva/MeibciKa+Y+AmwExgE3ZubDxc2QOzJzA3ADcFNE7KR2pu7S4rUPR8TNwHbgReD3M/MQQL1jDqdOSZKko0XUToCNbe3t7dnR0VF2GVJ9LuUjSeolIjZnZnu9ff0OxUbEf42I4/vZ/+6IuGC4BUrqQ89SPk/vAvIXS/k8eHPZlUmSWtBAQ7HbgG9FxAvAD4F9wCRq17SdDtwJXN3UCqWxzHVQJUmD0G+wy8zbgNsiYi61e9hNB54BvgZckZnd/b1e0jC5lI8kaRAamjyRmY8BjzW5FklHcikfSdIgNHq7E0llcCkfSdIgGOykVuY6qJKkQWj0PnaSyuJSPpKkBg35jJ23OZEkSWotwxmKfeuIVSFJkqRhG+gGxR8onuccuS8zr2xWUZIkSRq8gc7YrSye/7rZhUiSJGl4Bpo8sT8i7gbmRMSGI3dm5oXNKUuSNGJcb1gaMwYKdr8BnAHcBFzT/HIkSSOqZ73hnqXpetYbBsOdVEEDBbsbMvPDEfHlzPyHUalIkjRyXG9YGlMGusbuzIh4PXBZRLwmIo7v/RiNAiVJw+B6w9KYMtAZuz8DbgfeAGwGote+LNolSa3K9YalMaXfM3aZeV1mzgduzMw3ZOacXg9DnSS1OtcblsaUfs/YRcTkzHwG+MN6Q6+Zub9plUmShq/nOjpnxUpjwkBDsX8FXEBtGDZxKFaSjj6uNyyNGf0Gu8y8oHh+2coTkiRJai0DDcWe0d/+zPzhyJYjSZKkoRpoKLbnpsSTgHbgAWrDsacBPwDe0bzSJEmSNBgDzYo9NzPPBf4ZOCMz2zPzTGABsHM0CpQkSVJjBrpBcY83Zea2no3MfAg4vTklacx58Gb4wlvgU1Nqzw/eXHZFkiQdlQYaiu3xSER8BfgatdmwHwIeaVpVGjtcx1KSpBHT6Bm7fwM8DHwM+ANge9EmDU9/61hKkqRBaeiMXWa+AHyheEgjx3UsJUkaMY2esZOao6/1Kl3HUpKkQTPYqVyuYympFTiJSxXR6OQJqTlcx1JS2ZzEpQoZcrCLiCsyc81IFqMxynUsJZWpv0lc/tuko8xwhmJjyC+MOD4i7oiIx4rn1/TRb1nR57GIWNar/cyI2BYROyPiuoiIov1TEdEVEVuLx28MtUZJ0hjhJC5VyJCDXWb+r2G87wpgU2bOBTYV24eJiOOBK4GzgbOAK3sFwC8BVwBzi8f5vV76hcw8vXj83TBqlCSNBU7iUoU0FOwi4uqImNJr+zUR8d+H8b5LgLXF32uBpXX6LAbuyMz9mfkUcAdwfkRMByZn5vczM4G/6OP1rcELciWptTmJSxXS6Bm792TmgZ6NImgNZ5jzxMzcUxxrDzCtTp+ZwK5e251F28zi7yPbe3wkIh6MiBv7GuKF2jWCEdERER379u0b6ufoX88FuU/vAvIXF+Qa7qQRtX5LF+esuos5K77NOavuYv2WrrJL0tHktEvgvdfBcScBUXt+73VeX6ejUqOTJ8ZFxMTM/ClARLQBE/t7QUTcCbyuzq4/bPA9613Dl/20Q22I9tPF9qeBa4DfqXfwYuLHGoD29vas12fYvCBXarr1W7pYuW4b3QcPAdB1oJuV62pLWy9dMLO/l0q/4CQuVUSjwe5rwKaI+Cq10PQ7/GIota7M/PW+9kXEExExPTP3FEOre+t06wTe1Wt7FvDdon3WEe27i/d8otd7fBn42/5qbDovyJWabvXGHS+Fuh7dBw+xeuMOg52kMaehodjM/Bzwx8B84BTg00XbUG0Aema5LgNuq9NnI3BecT3fa4DzgI3F0O2zEfG2Yjbsv+55fRESe7wPeGgYNQ6fF+RKTbf7QPeg2iWpyhq+j11m/j3w9yP0vquAmyPicuBfgA8AREQ78LuZ+W8zc39EfBq4v3jNVZm5v/j794A/B9qKmnrq+lxEnE7trOLjwL8foXqHZuEnD7/pJXhBrjTCZkxpo6tOiJsxpa1Ob0mqtqhNLB2gU8TbgP9B7YzdK4BxwPOZObm55Y2O9vb27OjoaM7BH7zZVRWkJjryGjuAtgnj+MxFpzoUK6mSImJzZrbX29foGbs/BS4FbgHaqQ1/vnFkyqs4L8iVmqonvK3euIPdB7qZMaWN5YvnGeokjUmDGYrdGRHjMvMQ8NWI+Mcm1iVJDVu6YKZBTpJoPNj9v4h4BbA1Ij4H7AFe1byyJEmSNFiN3qD4w0XfjwDPAycBFzerKEmSJA3egGfsImIc8MeZ+SHgBeC/Nb0qSZIkDdqAZ+yKa+qmFkOxkiRJalGNXmP3OHBPRGygNhQLQGb+STOKkiRJ0uA1Gux2F49jgGObV44kSZKGqt9gFxE3ZeaHgQOZ+cVRqkmSJElDMNA1dmdGxOuB3ynWbD2+92M0CpQkSVJjBhqK/TPgduANwGYgeu3Lol2SJEktoN8zdpl5XWbOB27MzDdk5pxeD0OdJElSC2noBsWZ+XvNLkSSJEnD0+jKE5IkSWpxBjtJkqSKMNhJkiRVhMFOkiSpIgx2kiRJFWGwkyRJqgiDnSRJUkUY7CRJkirCYCdJklQRBjtJkqSKMNhJkiRVhMFOkiSpIgx2kiRJFWGwkyRJqgiDnSRJUkUY7CRJkirCYCdJklQR48suQFL/1m/pYvXGHew+0M2MKW0sXzyPpQtmll2WJKkFlXLGLiKOj4g7IuKx4vk1ffRbVvR5LCKW9Wr/44jYFRHPHdF/YkR8MyJ2RsQPImJ2cz+J1Fzrt3Sxct02ug50k0DXgW5WrtvG+i1dZZcmSWpBZQ3FrgA2ZeZcYFOxfZiIOB64EjgbOAu4slcA/FbRdqTLgacy843AF4DPNqF2adSs3riD7oOHDmvrPniI1Rt3lFSRJKmVlRXslgBri7/XAkvr9FkM3JGZ+zPzKeAO4HyAzLw3M/cMcNxbgYURESNauTSKdh/oHlS7JGlsKyvYndgTzIrnaXX6zAR29druLNr689JrMvNF4GnghHodI+KKiOiIiI59+/YNsnxpdMyY0jaodknS2Na0YBcRd0bEQ3UeSxo9RJ22HKnXZOaazGzPzPapU6c2WJI0upYvnkfbhHGHtbVNGMfyxfNKqkiS1MqaNis2M3+9r30R8URETM/MPRExHdhbp1sn8K5e27OA7w7wtp3ASUBnRIwHjgP2D6ZuqZX0zH51VqwkqRFl3e5kA7AMWFU831anz0bg6l4TJs4DVjZ43O8D7wfuysyBzvJJLW3pgpkGOUlSQ8q6xm4VsCgiHgMWFdtERHtEfAUgM/cDnwbuLx5XFW1ExOciohN4ZUR0RsSniuPeAJwQETuBj1Nntq0kSVJVhSe0oL29PTs6OsouQ5IkaUARsTkz2+vtc0kxSZKkijDYSZIkVYTBTpIkqSIMdpIkSRVhsJMkSaoIg50kSVJFGOwkSZIqwmAnSZJUEQY7SZKkijDYSZIkVYTBTpIkqSIMdpIkSRVhsJMkSaoIg50kSVJFGOwkSZIqwmAnSZJUEQY7SZKkijDYSZIkVYTBTpIkqSIMdpIkSRVhsJMkSaoIg50kSVJFGOwkSZIqwmAnSZJUEQY7SZKkijDYSZIkVYTBTpIkqSIMdpIkSRVhsJMkSaqI8WUXIK3f0sXqjTvYfaCbGVPaWL54HksXzCy7LEmSjjoGO5Vq/ZYuVq7bRvfBQwB0Hehm5bptAIY7SZIGqZSh2Ig4PiLuiIjHiufX9NFvWdHnsYhY1qv9jyNiV0Q8d0T/346IfRGxtXj822Z/Fg3P6o07Xgp1PboPHmL1xh0lVSRJ0tGrrGvsVgCbMnMusKnYPkxEHA9cCZwNnAVc2SsAfqtoq+ebmXl68fjKyJeukbT7QPeg2iVJUt/KCnZLgLXF32uBpXX6LAbuyMz9mfkUcAdwPkBm3puZe0alUjXVjCltg2qXJEl9KyvYndgTzIrnaXX6zA/L3iEAAAcaSURBVAR29druLNoGcnFEPBgRt0bESX11iogrIqIjIjr27ds3mNo1gpYvnkfbhHGHtbVNGMfyxfNKqkiSpKNX0yZPRMSdwOvq7PrDRg9Rpy0HeM23gK9n5k8j4nepnQ18d72OmbkGWAPQ3t4+0HHVJD0TJJwVK0nS8DUt2GXmr/e1LyKeiIjpmbknIqYDe+t06wTe1Wt7FvDdAd7zyV6bXwY+23DBTeKtPAa2dMFMvxNJkkZAWUOxG4CeWa7LgNvq9NkInBcRrykmTZxXtPWpCIk9LgQeGYFah6znVh5dB7pJfnErj/VbusosS5IkVVRZwW4VsCgiHgMWFdtERHtEfAUgM/cDnwbuLx5XFW1ExOciohN4ZUR0RsSniuN+NCIejogHgI8Cvz2Kn+llvJWHJEkaTZHp5WXt7e3Z0dEx4seds+LbdS8KDOBHq/7ViL+fJEmqvojYnJnt9fa5VmwTeSsPSZI0mgx2TeStPCRJ0mhyrdgm8lYekiRpNBnsmsxbeUiSpNHiUKwkSVJFGOwkSZIqwmAnSZJUEQY7SZKkijDYSZIkVYTBTpIkqSIMdpIkSRVhsJMkSaqIyKy3TP3YEhH7gH9u8tu8FvhJk9/jaOd31D+/n4H5HfXP72dgfkf98/sZ2Gh8R6/PzKn1dhjsRklEdGRme9l1tDK/o/75/QzM76h/fj8D8zvqn9/PwMr+jhyKlSRJqgiDnSRJUkUY7EbPmrILOAr4HfXP72dgfkf98/sZmN9R//x+Blbqd+Q1dpIkSRXhGTtJkqSKMNg1WUTcGBF7I+KhsmtpRRFxUkTcHRGPRMTDEfGxsmtqNRExKSLui4gHiu/ov5VdUyuKiHERsSUi/rbsWlpRRDweEdsiYmtEdJRdT6uJiCkRcWtEPFr8e/T2smtqJRExr/jvTs/jmYj4g7LraiUR8Z+Kf6MfioivR8SkUupwKLa5IuKdwHPAX2TmW8qup9VExHRgemb+MCKOBTYDSzNze8mltYyICOBVmflcREwAvgd8LDPvLbm0lhIRHwfagcmZeUHZ9bSaiHgcaM9M70FWR0SsBf5PZn4lIl4BvDIzD5RdVyuKiHFAF3B2Zjb7HrBHhYiYSe3f5jdnZndE3Az8XWb++WjX4hm7JsvM/w3sL7uOVpWZezLzh8XfzwKPADPLraq1ZM1zxeaE4uEvsl4iYhbwr4CvlF2Ljj4RMRl4J3ADQGb+zFDXr4XA/zXUvcx4oC0ixgOvBHaXUYTBTi0jImYDC4AflFtJ6ymGGbcCe4E7MtPv6HDXAv8F+HnZhbSwBL4TEZsj4oqyi2kxbwD2AV8thvO/EhGvKruoFnYp8PWyi2glmdkFfB74F2AP8HRmfqeMWgx2agkR8Wrgr4E/yMxnyq6n1WTmocw8HZgFnBURDusXIuICYG9mbi67lhZ3TmaeAbwH+P3iMhHVjAfOAL6UmQuA54EV5ZbUmoph6guBW8qupZVExGuAJcAcYAbwqoj4UBm1GOxUuuK6sb8G/jIz15VdTysrhoe+C5xfcimt5BzgwuIasm8A746Ir5VbUuvJzN3F817gb4Czyq2opXQCnb3OhN9KLejp5d4D/DAznyi7kBbz68CPMnNfZh4E1gG/WkYhBjuVqpgYcAPwSGb+Sdn1tKKImBoRU4q/26j9A/JouVW1jsxcmZmzMnM2tSGiuzKzlF/KrSoiXlVMTqIYYjwPcKZ+ITN/DOyKiHlF00LACVz1/RYOw9bzL8DbIuKVxf+vLaR2zfioM9g1WUR8Hfg+MC8iOiPi8rJrajHnAB+mdpalZxr9b5RdVIuZDtwdEQ8C91O7xs5bemgwTgS+FxEPAPcB387M20uuqdX8R+Avi/+dnQ5cXXI9LSciXgksonY2Sr0UZ3tvBX4IbKOWr0pZgcLbnUiSJFWEZ+wkSZIqwmAnSZJUEQY7SZKkijDYSZIkVYTBTpIkqSIMdpI0SBHx3MC9JGn0GewkSZIqwmAnSUMUEa+OiE0R8cOI2BYRS4r22RHxSER8OSIejojvFKuGEBG/FBG3R8TmiPg/EfGmcj+FpCrxBsWSNEgR8VxmvjoixgOvzMxnIuK1wL3AXOD1wE6gPTO3RsTNwIbM/FpEbAJ+NzMfi4izgc9k5rtL+zCSKmV82QVI0lEsgKsj4p3Az4GZ1JbvgtqC4FuLvzcDsyPi1dQWBr+ltpwkABNHsV5JFWewk6ShuwyYCpyZmQcj4nFgUrHvp736HQLaqF3+ciAzTx/VKiWNGV5jJ0lDdxywtwh151Ibgu1TZj4D/CgiPgAQNb8yCnVKGiMMdpI0dH8JtEdEB7Wzd4828JrLgMsj4gHgYWBJE+uTNMY4eUKSJKkiPGMnSZJUEQY7SZKkijDYSZIkVYTBTpIkqSIMdpIkSRVhsJMkSaoIg50kSVJFGOwkSZIq4v8DP8JIq3MOv0gAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot the fractional difference for 2013 and 2015\n", "_ = plt.plot(lanes, f_13, marker='.', markersize=12, linestyle='none')\n", "_ = plt.plot(lanes, f_15, marker='.', markersize=12, linestyle='none')\n", "\n", "# Add a legend\n", "_ = plt.legend((2013, 2015))\n", "\n", "# Label axes\n", "_ = plt.xlabel('lane')\n", "_ = plt.ylabel('frac. diff. (odd - even)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "EDA has exposed a strong slope in 2013 compared to 2015!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### How does the current effect depend on lane position?\n", "To quantify the effect of lane number on performance, perform a linear regression on the f_13 versus lanes data. Do a pairs bootstrap calculation to get a 95% confidence interval. Finally, make a plot of the regression.\n", "\n", "Note that we could compute error bars on the mean fractional differences and use them in the regression, but that is beyond the scope of this course." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "slope: 0.00447 per lane\n", "95% conf int: [0.00393, 0.00502] per lane\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmgAAAEvCAYAAADxWj0AAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9aXCd13nn+T9333BXXOwrQXARRS0UKVpWbHlRvKTTNtPtbruT6rhqnHimq12dDzOqbs10plJOauK0ejKV6aRdcSeeJO6knZTtqFVq24ptWV5kLSRFiaRIAgQJgAAu9ntxd9z1zIc/jg9FkRIlAsTC51d1C/e+97znPYAs+VfP85znKK01BEEQBEEQhK2DY7MXIAiCIAiCILwRETRBEARBEIQthgiaIAiCIAjCFkMETRAEQRAEYYshgiYIgiAIgrDFEEETBEEQBEHYYrg2ewHrSWtrqx4YGNjsZQiCIAiCILwtJ0+eXNJaJ6/33Y4StIGBAZw4cWKzlyEIgiAIgvC2KKUmb/SdpDgFQRAEQRC2GCJogiAIgiAIWwwRNEEQBEEQhC2GCJogCIIgCMIWQwRNEARBEARhiyGCJgiCIAiCsMUQQRMEQRAEQdhiiKAJgiAIgiBsMUTQBEEQBEEQthgiaIIgCIIgCFsMETRBEARBEAQA0BrIZICZmc1eyc46i1MQBEEQBOEdY8SsVAJiMb42GRE0QRAEQRDuTLQGlpeB1VUgHudriyCCJgiCIAjCnUWzSTGrVillra2bvaI3IYImCIIgCMKdQbMJLC0BtRqQSAA+32av6IaIoAmCIAiCsLNpNChm9TqjZV7vZq/obRFBEwRBEARhZ9JoAIuLjJy1tgIez2av6KYRQRMEQRAEYWdRrzNipjXFzO3e7BW9Y0TQBEEQBEHYGdRqFDMASCYB1/bVnO27ckEQBEEQBIC7MZeWAIdj24uZYfv/BoIgCIIg3JlUKmyX4XQC7e38uUMQQRMEQRAEYXtRqTBi5nYDHR2MnO0wRNAEQRAEQdgelMtAOs3dmJ2dO1LMDCJogiAIgiBsbUolipnPB3R1AUpt2KOePDWDJ54ZQWqljK6oH499dC+O3d+9Yc+7EeuinkqpjymlRpRSY0qpf3ed771Kqb9d+/4lpdTA2vWEUuqHSqmCUuqPr7nnubU5X117ta3HWgVBEARB2CYUi8D0NCNn3d1smbHBcvb4t05jZqUMDWBmpYzHv3UGT56a2bBn3ohbFjSllBPAnwD4OIC7APwLpdRd1wz7HICM1no3gP8HwB+sXV8F8NsA/rcbTP9rWuv71l4Lt7pWQRAEQRC2AYUCxaxSoZglEhsqZgCASgVPfPt1lGvNN1wu1xp44pmRjX32dViPCNqDAMa01pe11lUAXwfwyWvGfBLAX669/waADyullNa6qLX+KShqgiAIgiDcyeTzFLN6Hejp4UHmGy1mJkqXzSKVr113SGqlvLFruA7rIWjdAKau+jy9du26Y7TWdQBZAImbmPv/W0tv/rZS1/8npJT6vFLqhFLqxOLi4jtfvSAIgiAIm0s2S0lqNilm0ejGPk9rIJPhM6tVRujqdXSFrl+a3xX1b+x6rsN6CNr1xEm/izHX8mta64MA3rf2+pfXG6S1/orW+rDW+nAymXzbxQqCIAiCsEVYWaEkKUUxi0Q29nmNBrCwAKRSbNERCjGdOj4OzM7isQNB+K9ppeZ3O/HYR/du7Lquw3oI2jSA3qs+9wBI3WiMUsoFIAIg/VaTaq1n1n7mAfwNmEoVBEEQBGE7c3X0yuWimIXDG/vMahWYnaWcBYOUs6UlYGICmJ+nuNVqOBap4vePRNHtd0AB6I768fv/5OCm7OJcjzYbxwEMK6UGAcwA+AyAX71mzFMAPgvgBQCfAvCs1vqGEbQ1iYtqrZeUUm4Avwzg++uwVkEQBEEQNgOt2SqjXAZiMb42mmKRUToTLcvlgJkZplTdbvZTKxQ4rqUFqNdxLNHEsUfDwMMPb/z63oJbFjStdV0p9QUAzwBwAviq1vp1pdQXAZzQWj8F4M8BfE0pNQZGzj5j7ldKTQAIA/AopY4B+AiASQDPrMmZE5Sz/3KraxUEQRAE4TajNY9jWl1l0X/iZkrQb/F52SzFKxAA/H6+T6UYSQuF+JqbA7xevkolblBYXATOnOEmhe0uaACgtf42gG9fc+3/vOr9KoB/doN7B24w7QPrsTZBEARBEDaBZpNiVq1SzFpbb8/zKhWmMV0uRuzMWZ3JpE2txmKMnhWLHDc+DoyOshFuRwcjbZuMnCQgCIIgCML60WyyvqtWY7TM59vY51WrlDCt+SylgKkpm7bs67MCFo9TyBYXuc4LF1iDlkgAbW3AlSvAxYsbu96bRARNEARBEIRbp9GgmDUaFB6vd2OfVyoxIuZwMEJWLFLMlKJsxWKUrcVFpjrrdW4UKJWAc+eYcm1r4waFixeZ4mw0GFlrb9/Ytd8EImiCIAiCILx7Gg0bkWptpeBsJCsrrClzuShj6TSvBQLA7t1cy4ULjJ45nUxXZrOUuQsXWJMWDFLIjh+nqDmdbPHR3s7fZwscwi6CJgiCIAjCO6depwwBFDO3e+OedXV9mdNJiUqleD2RAO65Bzh/HnjpJdtLbX6eAnf5so2ihULA5CRlrdmkrPX12Y0EpRJr0Db69IKbQARNEARBEISbp1ZjKhNg4b1rA1XCPKvR4OdCgaLm8wG7djHS9dprrDHzeimNMzOsRzt/nvd7vZSxs2cZLfN4KHWJBIWsVmOac3CQfdLGx29PC5C3QQRNEARBEIS3p1qlLDkcGy9mpr6s0bCbDqpVHgF16BBbZBw/zmia1kxXmp+XL9sI2Py8lcmWFh68rhRFL59n9KxW43yLi4wEtrdv7O92k2z+CgRBEARB2LpUKrZVRXs7f24U2SzFqVbja3GRqdOuLhb0nzsH/OAHFCit+b3Xy92Xy8ucY3WVn0slpjW7uhghS6dZj9bVxWjZzAxr0mIxCmejYRvYdnVt3O94k4igCYIgCILwZlZXKT1uN+uyNqpwvtmkPBWLlLJCga9oFDh4kEJ48iTwyitcS7VKkQKYjiwWmdpMp5miBHhvXx+vZ7OUuQMHeIrBlSuMmHV1cVPBygrvi8eBI0cYZTMp1U1EBE0QBEEQBEu5TNnxeIDOzo0Ts1qNAlgsUrpWVph+bG+nTC0tAc8/z/UEAhTGfJ6fp6d5X6XCzQLlMgv9e3uZysxkmN4cHAT27qWUvfoqxa2/n/fNzzP61t8P7NvHqFyzSTnbyCjhTSKCJgiCIAgCU4LpNAvwu7o2bidjuUwxKxT4vlCgVA0PM4o1MgI8/TRFyeu145aWGOnK5eyxTA4HpWtw0KZiazUKXqXCerSZGWBggBK2uMgD0js6gEceYc1ZvU4hazb5HqAQbjIiaIIgCIJwJ1MsMuLk99si+o0gl6MAZrOMhjWbrP3au5eCdOIEW2CEQpSzfJ6vpSWmJFdW+CoW2cest5dj0+k3Stj4OFOi7e3A/v28Z3aW4jk8zN2fSlHutOY6ymXKYDDI68HgxvwN3gEiaIIgCIJwJ1Io2AavGyVmpr4snaYE1usUpd5eFv3ncsCzz/K7WIzjTdPbuTmmLzMZrrXZtLVlpjluo8E6tXIZuHSJactduxghm55mNK6rC/jQh3iv1pTBep2iB/C5wSDX5ffzs5zFKQiCIAjCbSWX4ysUAnp6NuYZpont/DzlSimmE9vaKEpjY8Df/R2jVYEAa9Dm55mWvHLFilm1ynV2dXHcygrFbXAQuOsu9jo7eZLf33cfo23T07xn/37+fk6nbcdRqTCVG4uxD5rHY/uiNZuMxM3O8v5NRgRNEARBEO4EslmKWTi8cWK2ukrBMc1ifT4W4ZsjoH72M9aFhUL8nE5TmLJZnoc5N8c5HA6eCGDEaXmZP++5h+NNtGzvXkbUpqbYgqO3l9GycNimMSsVyp7bzbRnWxuFraWFsphKUfScTrshYiNPRbhJRNAEQRAEYSdjzq4MhykwG0Eux/ox07csHmd9WSLBOrLvfpcSFo9zvGl4m0oxmpbJsLg/GGS61eu1mwn6+lhbNjICnDpFuTx0iN9NTPCeAwe449Tt5m5Mk8IsFjlfPM7rXi/Tn6USRXF0lBKptU31xuMUtk1GBE0QBEEQdhpaU3pKJUaJNiJiZhrFXrzIqFcoxNRjPM7o14ULwPe+x7HBIIUrlWLd2Ogoi/lLJUpVLMaIllKUPYeD0bFSieOWl5mybDQYZRsZoZA98gif63bznlKJEbx4nN87nZwzmeS4ixcZxXO7+V00yohadzeF7aWXmFaVo54EQRAEQVg3tGbasFymZJiI1XpSrzNadvkyRSeRsGIGUHImJlhwD9gO/vk8z82cm+P1lhZGx9xuRs9KJYrSrl1MYZ45w/To4cOMuI2PM8q1dy9TlR4P763X2X5Da+7k7Ozk/MEgxXRuDjh9mhIZCPC5ySTnXllhA9yXXuL4jg5G2qQPmiAIgiAIt4zWjDKtrlKUEon1f8bqKvD664yamWhZJEIRXFwE/uEfmMb0eln3VShQdKanmTIsFChVySRFye1mCtLjAYaGbDuMTAbYs4fRsqUlimBrK3D0KCNeHg/Xk8txbHc3o2tK8dXTw5+nT1Py3G4rX/39/O7kSX7vcjGC1tHxxh2eRi43EaW13uw1rBuHDx/WJ06c2OxlCIIgCMLtwRTQV6sUs40Qi6UlykytRplJJilmHg9TjefPU6aU4linkxG8s2cZaWs2GbUKh1nv1WxyrGmGazYG7NrFefN5ypfbTaFqb6fQOZ0UuqUlCtfAAGWw0bCHnF+4wPRpvc65kklG1GIxytrEBH+n1lZeMw5UrXIut5vrbGvjmjcYpdRJrfXh630nETRBEARB2G40mxSVWo3RMp9vfefXmqJz6RJFrK+PwhKPU6BOnWI9mdb8vLpKgZqd5W7KxUXKYjzO+7RmjZjfT/lZXOROz0CA0a9qlfPMzFCs7r2X97rdtudZpUJh6+/nNbNDNJdjijKXs812k0lbV/bii7w3GmWkzvR7q1a5Jp/PHpgeifB3yeVui6C9FRJBEwRBEITtgkn7NRoUM693fecvlylf6TSjTJ2dTGcGApSvsTFG7EwKU2uKzoULtug/GuXLtK3w+ThXucw0ZrNJsQqFKEOmYWx3N+UtEuHnXI7Pam+3qUmAa4pEKGWTk1yDaXzb28u05+nTXJ/Px/Sl+TtVq4yuBQKMwiWTtkGuaaQbCHD+jTqD9CreKoImgiYIgiAIWx0TRWo2bU+x9cQU0jcarOGKx5nqq1YZKZuYsA1uazXK0sICxWxujjKTSFB6Gg3KVzDIn2bnZjhMuVKK8xaL/L6jgy+Ph+I3N0ehGhpiFKtatYecj41RIPN5Cll/vz0F4ZVXKK+m8Ww0yuu1mj1MPRjkd62tFMiVFUqiaQ1Sq/Ga1kzBbjAiaIIgCIKwHTEd+QFKxXo2UG00WENmdkf29VFggkFGk+bmGDVbWqLgNJuMkJkdnNksBSuRYLG9Obg8Huc9hQLX29Nja+NKJUpQayulzDSUXV6msPX28lWrUbT6+/n+Jz+hEHq93MXZ2Ul5e/112xQ3GmVEzOXiPYUC34fD/K69nRE1c8anw8HomVL8fRsNrjMavS3RM+A21KAppT4G4I8AOAH8mdb6S9d87wXwVwAeALAM4NNa6wmlVALANwAcAfAXWusvXHXPAwD+AoAfwLcB/JbeSTYpCIIgCDeiVqPkAFY61ot8ntGylRWK0j33WIHKZChfMzPs1N9oUH5mZ7kbc3qaa4tGgd27GX2KxWxK88oVzh+JUKSM6GSzlLWODttQtlDg+FiMjWg9Hlvwn0gAL78MHD/OZ+zZw1q19nau7+WXuQ6zm9Q0m83lGHGLROwGg1DojelVswkgk6EYer185nr+jdeBW16NUsoJ4E8A/CKAaQDHlVJPaa3PXTXscwAyWuvdSqnPAPgDAJ8GsArgtwHcvfa6mi8D+DyAF0FB+xiA79zqegVBEARhy1Kt2i775kii9aDRoGSdP2/Td2YXpEkrptOUn+VlykqxSIGanWX0yuWiyHi9XF8ySdlJpTguHOacfj+FqVq1R0vt30+Jq1Y5l9PJsXfdRQFraeHnyUnguee4ho4O9kDr7ubfZGSEaUyfz9bGORwUvfl5zmFaZkSjjD5mMpRDv5/SZ9K0pk/cFjjS6Uashy4+CGBMa30ZAJRSXwfwSQBXC9onAfzO2vtvAPhjpZTSWhcB/FQptfvqCZVSnQDCWusX1j7/FYBjEEETBEEQdiKVCiXE5WLUZ73EzBxpNDNDSRkYYArT4aCkZDI8x3JsjGswOzEXFihe+TzTn/39FB6TmiwUeF8mQ9Hp7bUHjJu6rmSSkS+Hg9empnjvAw9wnNPJ9GejwZ2WP/gBhenuu4EPfpAyZ9p4OJ2MivX0cP25HKXMbEAYHqaAKcVnpVK2ZUapxN/VRP5aW9fnb7vBrIegdQOYuurzNICjNxqjta4rpbIAEgCW3mLO6Wvm7F6HtQqCIAjC1mF1ldEit5tRofWofTK90S5domCFw0whulyMnpVKfO7Zs4x8uVwcZ4RrdtbWdPX18fvubs5z5Qo3BkQiTC3GYhQ483sEg7afmTk4PRymqPl8NoUZj7N+7JvfZKRueBj4pV9imtMcx9Rscu6BAc5nJHZ5mXMcPMiImdvN9ZvNCuEwfxaL/H1MXdw2Yz0ETV3n2rW1Yjcz5l2NV0p9HkyFoq+v7y2mFARBEIQtQrnMlKLHs35iVi5TUiYnGUmKROyxR0pZwXnlFUqOz2c3ASwvU2YcDka5TAF9d7fdGGCOdertpQR5vYy05XIUoHvu4bOyWW5sGBhgnRpA0ert5fNefJHrTCR4OkAyyV2ir73GqJnPxyiikarFRf6tYjHu7Ozs5JjVVa692WT0LhRiZC+b5e/evb3jOushaNMAeq/63AMgdYMx00opF4AIgPTbzHn1ya7XmxMAoLX+CoCvANzF+Y5WLgiCIAi3k1KJsuHz2U76t4Lp3zU3R5ExKVKv1xa91+tsOnvmDCNYtRqFyNRulcuMfHV1URi7umy07OxZitHAAKNWLS0Uo8VFjjWpzXyestTZySazgI28NRpMU/7wh5S8ffuA97+fz56YYETO4+HcplXHwoI9z7Onx+7arNf59zOCGQgwUlYs8vv1+JtuEdZD0I4DGFZKDQKYAfAZAL96zZinAHwWwAsAPgXg2bfakam1nlVK5ZVS7wHwEoBfB/Cf1mGtgiAIgnD7Mek20+n+ViXCRI/m5ihYLpeNxAUCtsP/8eNMXTocjF6ZdhaLixS1ZJJSFIlQwgoFG4Fra6McJZMUqOVlzhWJsOi/Xqdwulz8bFptJBKUtslJ4LvfZYSttxf4+McZxbtyhU1mAYpkfz8lMJdjrZzTyb/RXXcxPak1o2IzM3yGOd5pdZXPMb3Vdhi3LGhrNWVfAPAM2Gbjq1rr15VSXwRwQmv9FIA/B/A1pdQYGDn7jLlfKTUBIAzAo5Q6BuAjaztA/xVsm43vQDYICIIgCNuNQoFF68HgrYuZ1pS8dJov0yvMnCgQCDCiduUKDwM3uxWnp60grqxwbDzONZlzKqengVdfpXz19dlu/aUSI10OByNz3d0Uv2KRNWjhMNfV0sIIWzrNVOX0tD2yKRTims6d4/qcTq65o4O/w9QUn9HeDhw5wnnM7szZWf7uPp89yNzj4djb1Ktss5BGtYIgCIKw3uRyjPq0tDAKdCtUKhSflRVKi6nT8nopSIEAr585w1RiqcQUYSbD90tLnCMcpoyFQoyWra4yjai1PdbJnFKwtMTnRSKUuUaDr64uRtVWVxkNNDsiJya4E7TZZN2ZqTcrFilV5hzOzk7OPzXFv5FJn5reaJWKbRpr+qKZczdNr7UdhJwkIAiCIAi3g2zW9v4yZ0q+G7SmIOVyjFiVy5QXn4+CFY9TVmZmGPm6coXjFxcpayZiBlCCIhFGnRIJpkXn5rhGc6h4NMpnLCzYnZYuF9+bprSNBq9FIlzD/DylLJvlPMPDb3xus2k3FiQSHD83Z9t2mD5spl9ZtWoFzAhaNLqle5XdKht+koAgCIIg3NGY6FYkwujRu6VaZbTMRMpyOQqM10uRisf53fnzPA1gYYHiYyTOrCMQYA1ZNEoRqlYpc6kUrx0+bFtULC5yE0EoxPVrTREcHmYatFzm82MxCtjly5zH52NLjcFB256j2eTvYXZR5vOMrs3N8fOjj/I5pq7MnJZgUr9K8TnrfdboNkQETRAEQRDeDaYmrFSi9PT0vP09N5onm6XMrK6yLiuTsY1VOzspSlNTwLPPUqZSKfvsfJ73mwPJ+/uZioxGKXhnz1KmurtZX2aiZTMzFDeThvX7+X1PD+fz+ex3s7NsGlsus3bswAGKYKnE36HZtLs6HQ72MjP9yh5+2J55WSxy7UbklGLULBrl84SfI4ImCIIgCO8ErRnlMscFvdsmqCZaVqnYHZFLS5SWri4KVbnMhrMXLjByNTnJa+UyJapctmnEaJSCpNQbo2VHjjC96XZz/vPnrQzFYkxP7t/PddRqvN7fz8jaxYtcY0sL52g0+H2txvcOh91QMDHBjQDhMHuimV2l5ngnc4/TaU8GCARu+R/HTkUETRAEQRBuBq0ZFVpdpZQlEu9uDlNXZmRoYYGpyVCIPcJaWylizz/PaNnICGWpXn/jveEwJai7m6JVKlHmjOAdPcoxxSKFzRzlZI5HOnCAz0ynKVGtrbZ9x8gI5SoYtI1rnU4KVr3OeXftYnr18mWmQIeGgF/4BVu7lk7biKDTyevRqD0SSnhLRNAEQRAE4a0wRydVqxSzd3OWY61mW2PU6xSYK1cYAevsBN73Pl4bGbFiduECpcvcW61ScsxuzO5uRsUWFjg2FAIOHbJ9wZaXeR3gvckkJWrPHtaEVasc191to2WFAlOVwSDnNvViJpU5NMT7Ll7kc3t62N/M57PyWShQBgGKWyLBCNwO7FW2kYigCYIgCML1aDYZTarVKBnvpkbKRLxMzVU2a8+/3LuXgpNK8ezJ6Wng1ClGu4zs5HJ8HwxSlvr6KGiNBsc3mxSsw4eZMjSHmJfLthXGrl383hyDtLxMiTPvMxnKk9Zv3CGayVC0enspbWNjbOWRSDBSZlK7pRJ/h3LZnp9pdo7u8F5lG4kImiAIgiBcTaNBMWs0bBPYd0K9TvGp123UaHqaUapolNEyj4dRsqefZorw1CnKmEkN1uuMYLW02GiZx8PvLl+mBB08yB2aZr0zMzal2NMDPPggo2aTk3y+OcqpWGQErNGgXPl8/D3DYbuDNJnk/BMTjOoFg0yJmma7tRrTm8Ui3wcC/N12YK+yzUIETRAEQRAASpE5fNs0bH0n5POULKX4Wl2lTFUqFKmjRylAr7zCKNfZsxSlWo2ik8/bhq6mviwc5lzz80wttrezVUUkwvHj46xfKxYpckeOMM05Pc1nmSL+QgGYm8OT01U8MamQqrvR5fXgsXvacCy+JmpeL1trpNOUvcVFttB4+GFKV7PJ78xJAsGgFUiX6MR6I39RQRAE4c6mXqeMABSzd9IY1RzeXa1S6JRibdbUFIXnvvsYnRobA/7H/2B685VX+DxT29Zo2FMBolF7lmW5zPYWLhc3D+zeTZlbXqZALS/zeXv3Ah/8INdw5Qrrznp6GBmrVChrhQKenKnh8eU4yppRvZkK8PjJPPBwG44NtlMm02kK3Uc/auvK8nnbNNfn4zpNNE7YMOQkAUEQBOHOpFazbS1M5/ybxdRzmd2J5pBx01X//vspTOfOcWfl2BgjZiZStrpK+YlE7HmaLS2UQ9M+Ix7nWZbhMOdfWOB683lG0h55hNI2MsLvzaaBWo0RsUKBkbe1OrSHnytiptR806/S7dV4/tEIDyc3O1PLZcpaJmNPD4jH33m6V3hL5CQBQRAEQTBUqxQdp5Pd9G+2ZsrUh1UqlB6Xy9Z+NZu2zcTMDPDDHzIidfEi67iMMFUqtl4rFmOaMBCwzWoBphmHhnhtcZFRsbk5PvOee4APf5ifp6dZx5ZM2jYbKyu2o39fH9ezdvxTqnT9gEyqolgXZ+rKlpd5fyzGdfj96/N3F94RImiCIAjCnUGlQqFyuRiBulkxKxYpPg4HZapUYkQsk2HK7/77KVojI8Bf/zXFaWTEHjheqVC2wmFG6kykzO+3KdJIhDVq0ajtZzY3x0hbby/wm7/JMa+/Dpw8yRTm4CDnzmY51kTvDh/mHLkc19jdDeRy6PIuYKby5l+vK+xhWnR1lc/v7+fvI2wqImiCIAjCzmZ1lVEht9t2t387Gg3KzeoqpczvZ2RqdJRSlEwyOlUssqbs/HlK0qVLFC6zI9Pj4dholPOa/mLFIscMD1O0TD3aq69S7AIBFuc/9BCjcJcuMcXY10epq9XY2mJ+nsK3Zw/nzWQ4TzLJOjazprY2PPaxvXj822Mo1xo//zX9ToXHDsVZUya9yrYUImiCIAjCzsTUUXk8Ny9mJk2oFIWlVGKq0uyw7OnhXBMTwLe+xfTj7Cx/Li1RnppNRsvicQpZo8GfSnGeaJQ7LaNRitrYmD2fcmgI+LVf4z0XL1LYzKHnpmZuZobfDw2xBs30Wouvidb8POeLRLhJIZkEHA4cW10Fcnk88dIcUqUmusJePPaxfTh26F2eISpsKLJJQBAEQdhZlEoUM9Pf6+2iQs0mI0/lMiNXbjeFK5ulFPn9dlfka68xYmZ2Uk5O2ka0LhdlKJnkXF4v56vXWfc2MMD0oVKsC5uaolzF46wB27OHslYqcd3xOO9bXWXaNJ9nurO7m9dLJYpfOGyjfR4Pn2N2WdbrjKItLnJ9HR2cWxrIbglkk4AgCIKw8ykWKSqBgG2o+laUyxwPsCC+UmG0rFqlwLS0UGgKBRb9m7qzqSkKXL3Oe0Mh20i2XOa1ZJLy5XIxjRmJcP7z5+0GhbvvZhozm2XUa3SUzzM91C5dYjQvEmGkzOXi71guU8ocDn6uVBjV6+nhms3mgtlZvm9rY9NZ6VW2rZAImiAIgrC9yecpOcEgU4FvJWZXR8v8fsrc+DhlyhwIHg5Tiu7Xjo8AACAASURBVMbGgBdeoOysrFCglpc5jzljsrubz9eaz240+L67my+tKV/msPLOTtaVxeMUqFqNMtfSwu9NDZlSjLaZCFyjwXXV6xRGI149PbZ7vzlGypwZ2tPzznq6CbcdiaAJgiAIO49cjlLS0kIZeStWV5n2BBgtazZZR1arMXXp9Vq5e/FFRrqyWXuIeKlkDw/v6+M9hYJNR+bzvNbTY5vMnjnDMX4/a86GhriOlRUKVmenbQQ7OmpPCjDRskqFP51OPiebpaSZfmVeL6+PjPBnOMxo3bs5M1TYcoigCYIgCNuLbJZS09LCmqwboTUjUlefN2nOpXS7KTgeD+eZn2fR/9wc5erSJRbaNxoUrL4+ClW1yu99PkbZALbOOHCA0pdKsRWGw8HdmR/4AO/P5Shm8ThfuRyjXZUK5+nq4nqaTf50OimP5TLXODREsQyHeX1ykvP5fIy0tbTclj+9cPsQQRMEQRC2B5kMa64ikbeOmFUqjJaZtKNSTGM2GoxuBQI2MnXqFLv9ZzKUpgsXbO2YOQ+zpYWRsEKB97a08LvubspULseNA40G5eu97+VPh8NuPOjp4fulJf4MBjnG7aaAORx8XyhQwLSmfMZifCllG9O63fxu9+7b97cXbjsiaIIgCMLW5eooWDR6YzHT2h4aburDZmYYzXK7KXWlEsfmckxjzs1RiGZnGTGrVilT/f283+Ph2FqN730+SlkiwbnMLky/nzsw29oob1dvMmhpYcRveZlzh0L8PTweSqDbzbmU4j1dXVbKfD4e7/Taa1x3ZyfwwAPSq+wOQQRNEARB2HpozShYuUxZicevP65atUcTRaOUnMlJvlpa7O7MYpHRsfPnWVdm+pstLPCeeJzF+pEIhazZ5H0mNWpSqakUo3Gm6e3+/RQvn48RuUCA0bV8nusyombq3MyrVOL8Dgefa87iNEJnatKSSZ7HebOnHgg7hnURNKXUxwD8EQAngD/TWn/pmu+9AP4KwAMAlgF8Wms9sfbd4wA+B6AB4N9orZ9Zuz4BIL92vX6jXQ6CIAjCDkJris3qKqXJHN597ZhsltEvr5fj5ucpX14vJadWowAtLwOnTzPatbJC6Rsft6lHc6yR6TlWqdhUaH8/ZWtpiQedN5scd/CgjYb5/YyGud2cM5fjGhMJfna5KGeBAH+nYpFjEwkKZSBgJXJykuIWiQD79nFe4Y7llgVNKeUE8CcAfhHANIDjSqmntNbnrhr2OQAZrfVupdRnAPwBgE8rpe4C8BkABwB0Afi+UmqP1tqcQ/FBrfXSra5REARB2OKYo45Mi4jW1jePqVYpWI0GJaalhRGt6WnKUjJJAVpaYspydJTpy3KZ4+bm2KYimWRELBTiy5zR6fUyKjYwwHkmJzlHKMSCfxMdC4c51uHgvaurfIZpUFuvUyLjcf5e2SylzJxz6fPZaN/UFJ/j9/O5cgamsMZ6RNAeBDCmtb4MAEqprwP4JICrBe2TAH5n7f03APyxUkqtXf+61roCYFwpNbY23wvrsC5BEARhq9NsUo5qNUqZ1/vG77W2Rxl5PBSzTIYC5nbbczLzeUbGLlxg7dn8PCNsU1O2BYYp6jfPyWYpcIkEcOQI55mYAF56yR4PNTjIsZEIhUspSlixyLSjEcls1u7IdDrtoeyxGJ/r89mI29wco30ulxT7CzdkPQStG8DUVZ+nARy90RitdV0plQWQWLv+4jX3dq+91wD+QSmlAfyp1vor67BWQRAEYSvQaFBiGg1KzrXpvFqN0bJ63dZmLSxQqEIhilCjYTvwT0wwklYosMZsYYERt9ZWClgoxB5jWlPajBw99BDHnj9PWUwm2WfMyF9Xly3kT6et4Pl8tpbN76fMZTJ8flsbTwnweu0GhXSaa9SapwXcd58ctyS8JeshaNfbTnLt8QQ3GvNW9z6stU4ppdoAfE8pdUFr/eM3PVypzwP4PAD09fXd/KoFQRCE20+9TjFrNq8vZiZa5nJRqnI5K1QeDyNSxSIL/E0Kc26O10zUzOGwEpVMUsxmZ9lwtrWVUub1sr3Gj39MgRoctIX+7e18zuoq065mnr4+ipY5e7OtzUbTwmHWjfl8jLJFIrx/dpYRvUTCip8g3ATrIWjTAK7uFNgDIHWDMdNKKReACID0W92rtTY/F5RSfw+mPt8kaGuRta8APOppHX4fQRAEYb2p1xlxAihJV4tKvU4Rqtdtn7KVFY4PBvnZ6aTsnD/PBq/m8PClJZtebGnh3KbxayhECZubYxuMo0eZBj11is/v76d0OZ0ca1p4pNOUvbY2RrqyWa6n0eB6fD5Gy/x+il0wyMhYIEA5m5/nMyMRNpj1+2//31vY9qyHoB0HMKyUGgQwAxb9/+o1Y54C8FmwtuxTAJ7VWmul1FMA/kYp9YfgJoFhAC8rpYIAHFrr/Nr7jwD44jqsVRAEQbid1GqUKKUYhbr6wO58nhEyp5OSY9KTWlOgzG7MkRG+UilKmmlhYdpU+HyMeHV1Ubrm51mLFg6zaSzAY5eefZbSddddfKbLRSkLBN5Yj3bvvZx7cZHRr1CI86+s8J7+fp4cUK9znW43hW1hgb+HaW4r/cqEW+CWBW2tpuwLAJ4B22x8VWv9ulLqiwBOaK2fAvDnAL62tgkgDUoc1sb9HbihoA7gX2utG0qpdgB/z30EcAH4G631d291rYIgCMJtolqlmDmdlCLTx6teZ4SqWrXRqGKRomaatzqdlJ2zZxnxSqUoS9ks76tWbad/c6pAMkkJ+8lPGLX65CeZAn3pJT7H7J50OrmDsr3d7viMRIC9eylaS0vcVWl2iZoTBHbtssc5AZynWmWkzNSlDQxIvzJh3VBa75ys4OHDh/WJEyc2exmCIAh3LqZlhctFaTHCUihQsJxOpvyKRQpOs8mxzSZTiKOjtuh/dpb3FIsUO8DWonV1UZryeXba93qZwlSKpwdUKkxfJhIUL6eT9yjF6JvTSSlrb2fELZejwJkNAdUqxW9wkM9tNHivaZuhFKNq5sQBQXgXKKVO3qjPqwiaIAiCcOuYgnq3m2LmcFBq0mnKUiDAz9Uq05ZOp5We+XkW8I+NUczM7kgjZc2mLbwfHqZUnT3LOrSeHtaJjY4CCwt4su0AnogfQgpedKkaHkvkcazTyRRktUppO3DA1q6ZprT5PAv+29spbj4fP5tnl8u8PxKhlAUCm/rnFnYGImiCIAjCxlAuU8zM+ZcOByNeKyt87/VyDMDaMq2tqI2NMYU5OmpbZJTLHON08qc5n3LvXoreiRO8fs89jFxdusTr/f14smUIj9cHUIZNM/rRxO+353DswwcphouLfLbfTwk0OzD376dYFov2qCetuR5zfmY4LHVlwrryVoImZ3EKgiAI7xzTF8w0gG02+Xl11aYKTb0YwGhYs8lo2ZUrlLKREUayKhV+73RSuioVCtRdd1H6RkeBH/6Q7w8dohCOjLBx7NDQz5/3RKH/DXIGAGU48EQ2imOTk5y72aQ81mqMxg0M2JMAFhaYQjXHQLW0sHZN6sqETUAETRAEQbh5ikWmCwMBilm5zFoxpShlSlnZAvh9ucy+ZXNzLOQfH2e0TGsKk89nI1atrYyW1Wocu7pKSbr3XkbZRkb4+dAhCpffz/sLBaQa1xep1CpsRG9w0O7szOe5AcHj4fPcbs7X1fXmEw0E4TYjgiYIgiC8Pfk8a7ZMG4mVFbag8HgoPvW63YVpDhyfnqYAXb7MfmSzs7yutT1E3IjavfdyJ+bEBHditrRQlOp1Rt2aTW4KCId5r99vI3TJJPCe96Dr6XnMlJpvWnqXD8Cjj/LeTIZrMrtF3W6uo71dzsEUthQiaIIgCMKNyeUoZi0tTDFmMoyIuVx2V6PZlVkoUOTGxljrdfYs5Syfp8RpzbHNJiUuGAQOH2aU7PXX+WptZepxaYkp064uRur8fj5Ta9t/7P77WRu2tAQsLuKxQYXHz2uUm7ZOzO924LEP9vN3WF5m5M/nY4QsGGTRv9SVCVsQETRBEAThzWSzFCvT2b9QoEgZKfN6GTErFChIU1OMTE1PA6dPM7oGUKb8fo43Ozj37WNvtEuXgJ/9jLIVi1G2ikVG2QYHec3Ik9k4sHcvz7KsVrnGWo3CmErhWCAAPDiIJ86XkcpW0BV04bGDIRzr8drGtF4v53XJ//0JWxv5X6ggCIJgyWQoSV6v3ZFp6smUokRls4ysLS8zQpZOM4V58SK/M9GycJj3lkqUvPvu4/uLF9np3+9nVM7smkwkKGZut53D4WDNWX8/r6+s8Bn1OmvaajUe4/TpTwNeL47NzOBYX5bPCwR4v8vFdft8m/u3FYR3gAiaIAjCnY7WVsyMFNVq9vtgkKKVyTBCZs7CTKVYtD8+bhvNtrTYNGa1SrFKJDjmhRc4tzk6CbBF+ZEI73c67WHnw8N8dqXCKJnptZZOc96PfISRuNlZRuM8Hj6/rY0yGQ7zsyBsQ0TQBEEQ7lS0puysrPCziZqZIv5YjGlOc3D55cuUocuXKVzLyxQhrdnywuOhRHm9LOgvFtki4/x5Rr9MRM3jsQ1fAc7hcHCO7m5ed7kYbVtdZQp1YYH3HzoE7N7NKNr0NDcQxGJMe5oDyyMRzicI2xgRNEEQhDsNrSk98/OUKb+f1x0Oyo3HQ/mamOBrZoZRsytX+LlaZYoxHGY0zOzi7OhgxGtigtEys9syFKKEBYPccel2M1rmcDD12NbG6+a0AdOrbGaGY/bvBz7+cV6fmODRTtEon9do2Bo2t3vT/qSCsN6IoAmCIGwTnjw1gyeeGUFqpYyuqB+PfXQvjt3fffMTNJuMgM3PM/VndjD6fBSccpkRtbk5Hhh+5Qp/plIUutVVRsDicd5frVLwenooVKOjnMPn4xiHw0bLwmFKnMPB9+3tHBMMUqxM1/65OUbOenp44Hk0Sik7d46yl0zaeSIRK5eCsMMQQRMEQdgGPHlqBo9/6wzKtQYAYGaljMe/dQYA3l7SGg2mJVdWKDzJJEXLFM6n02yNYaRsfJwSNztr22dEozzH0uOhICWTjJBNTgLPP2/F69pomWlcGwiw1qylhWOMmFUqfNbyMiXxwQd5OkAqxbXMzjLC1tJiW3q0tEhrDGHHI4ImCIKwDXjimZGfy5mhXGvgiWdGbixo5TLFa3WV8tPVZaNlpgFsKsUI1eXLfL+8zJfpvN/aSjFrNilZ8Tjr0V591fYy6+5mNM3l4vehkP3c12fr04JBftdocP75ea5n927gH/0jiuDMDNt0JJNAb69tgBuLSV2ZcEchgiYIgrANSK2Ub+66KfyfnOROzNZWvmIxylAux52Xk5OUssuXuTtzefmNTWl37bLpw9ZWStWFC8CpUxSv7m7+XF2lOPX2MqrVaDB92dtrD0s3cpXL8XmNBmXxH/9jPmNqiinMWIzz1mpW9qSuTLhDEUETBEHYBnRF/Zi5jqR1RdckqlbjTsdUitGutjamJaPRNx5SPj5OQVtctP3MzFFKySTFDKCkBQK858c/ZkQsHgcOHrQNYsNhilapxOjY3r02QtbSwu8LBZsqTSSAI0d4z/w8d2G2tLyxrsw0xhWEOxwRNEEQhG3AYx/d+4YaNADwu5147P19TFHOzbGurKODcub1MsV57hwjVOfP8+fKCoVqaYnvfT42hzUtL+JxStvZsxQ+j4ed/10uSl06zQJ+015jYICf63VeM7JlUpihENOcu3dzPcvLXIfZIKA11xAOS12ZIFyFCJogCMI2wNSZ/XwXZ4sHjx2K4RgWgJUAm7rG4xy8sMDas0uXGC1bXqYcmRqvep1Cdv/9TDGa8ykvXGDBf6NB0XvkEY6fn6f07dvHvmixGHD33Uw/lsu2a38mw4ib18u13HMPr6fTnCMS4bzNpm2NYU4pEAThDYigCYIgbBOO7Q7jWNteRqfyeSDkA4YOULKqVUbJJiaAM2eY6jTnZ5rCf7eb0azeXopRIMC5XniBEbhgELj3Xka9xsZ4JNPgIOvCtObuyr4+PtvsAnW7KWb5PKNgvb2sWTPd/91uymCzaZvfejyb/acUhC2PCJogCMJWpl5nBKpY5KtcZnRq3z5K0vIycPIko1+jo7YtRqHASNrqKuu8HnjAdugHOP7sWUbLBgaAT3zCttlIJBghK5cZObvnHtaHFQq8NxJhGnRhge/7+niP12sPOzeROXPkUjC4aX9CQdiOiKAJgiBsRQoFWy/WaFC02tttinB0lFGuV16xB4hXq6wty2Y5prMTOHCA4tRsMtL18sscEwoB73kP5e3sWc6zaxcL+JXiAeQDA5RD08fM72ekLBLh3IEAI2Jm16XLZY9ZCgZtI1xBEN4xImiCIAhbhUaDQpTPU3ocDrsjMx6nYD3/PIVqbIxpxEKB0Swjcx4PpWzXLspRLsfWGCMjnHNoCHjvexkpGxlhUf/ddzOFmUgwWgZQCAsFzlet2p5oTielzOvl/B4P05jmjM1YzEbpBEF418i/RYIgCJtNsUj5KhQoO6Z/mOnMPz4OPPssG7jm8xSxYpHpzVKJopZMAg89xChbNssNAq+/zjRkNAocPcqo1ugoha2nhy0vnE4bLcvnOa8RLI+HETatee/Vxzc1GrzX1JV5vZv6JxSEnYYImiAIwmbQaFDKTGqypYUiVatRhppNm3qcnuYYEynL5ShlTicjYvfeS3G6dInjJyeZlhwYoJilUmwQGwpRxkzz2Lvu4vtczrbpCAa5DoARsrY2pjLNgegOhz3kPBTa1D+hIOxk1kXQlFIfA/BHAJwA/kxr/aVrvvcC+CsADwBYBvBprfXE2nePA/gcgAaAf6O1fuZm5hQEQdiWlEpMY+ZyjERFoxSjYpFCVC4DP/oRo1/VKsfncrynUuHLnFk5PEx5+9nPmLJcXqZgHTrEerHpac7T2sprgQB3Wfb2cm5znFNLCyNv5mzMcNimM2s1yqTLxbWaGjNBEDaUWxY0pZQTwJ8A+EUA0wCOK6We0lqfu2rY5wBktNa7lVKfAfAHAD6tlLoLwGcAHADQBeD7Sqk9a/e83ZyCIAjbA1Ogn05TjCIRRreqVVvQPz5uC/i1tudhFgq2Hm1gwPYWO36cr/l5ft/ZCXzoQ3zGwoI9fikaZRSsu5vSNTdnO/j39jKCVqnwZ1eXbdlRrfJaMCh1ZYKwCazHv3EPAhjTWl8GAKXU1wF8EsDVMvVJAL+z9v4bAP5YKaXWrn9da10BMK6UGlubDzcxpyAIwtbGdM43HfvjccpONssGsNks68rGxhilKhZth38jSa2tTEUODzNN+YMfsKO/icCZBrWLi5S8cJgtOJJJXu/v5zpMjdvAAMWr0WDELh7n5oB6nRsDymWuta2NPwVB2BTWQ9C6AUxd9XkawNEbjdFa15VSWQCJtesvXnNv99r7t5tTEARh63F1tMykI3fvphQtLbFAf2aG9WUrKzZatrTEdGa9zrG7d1O0fD7gpz9lGjOT4fzBIHD4MMeasy7b2hhFM4ejJ5OMlo2PM3q2Zw+f1Wiwdqyri59LJc7h9zPN2dKy2X9BQRCwPoJ2vSY3+ibH3Oj69Qocrp2TEyv1eQCfB4C+vr4br1IQBGEjWV1lFGtlhbKTSFDOABbtX7zIHZQzM7yWzdoDyysVXmtr44HjnZ2855lnrOg5HLze1mYjbcGgPRkgGKR01WoUPtOA1uD18n6fj2vNZu3OzGhU6soEYYuxHoI2DaD3qs89AFI3GDOtlHIBiABIv829bzcnAEBr/RUAXwGAw4cPX1fiBEEQNgStGdVaXGQ6Mh63uySrVUbLTp/m8UvVKuUpk2HdWKnEz6EQ+5YNDXHMq68Czz1nu/Y7nZQ2r9dKXSLB6FprK+9vbbWHk/f28nvTm6yzkyJmom2mZUZ3N8cIgrAlWQ9BOw5gWCk1CGAGLPr/1WvGPAXgswBeAPApAM9qrbVS6ikAf6OU+kNwk8AwgJfByNrbzSkIgrA5VCosxE+nKTvt7Sz8B3jtuecoZ9ksrxWLTDfmcoxeOZ2Mdg0Ps45sbg74znc4vl7nPaZ+zJy72WjwXEwjXOYcTHMMVEcHP7vdFDbTxd80sg2HGVHz+zfjLyYIwjvklgVtrabsCwCeAVtifFVr/bpS6osATmitnwLw5wC+trYJIA0KF9bG/R1Y/F8H8K+11g0AuN6ct7pWQRCEd43WNkpVrTJKtX8/BatSAc6d46HjFy/anmFmvNmpGYtxF2YkwmL8s2eZ8iyX+QynkycABAJMU6ZSFLK777aHjMfjHG8OHx8a4vhwmN8BdndoSwvvb2mRI5cEYZuhtN45WcHDhw/rEydObPYyBEHYSVSrjHCZaFlbG2VIa0bRzp1je4zlZY5vNChdmQwjW34/pcvsplxepsTl85zb7bYHjptoWEsLo2WJhO3mHwhQ+pxOpjXjcd5nmspqzXsDAXvd6dy8v5sgCG+LUuqk1vrw9b6TxjaCIAjXYqJls7OUotZW1om53ZSg115jA9gLF/jZ6WQaMZWyac32djaT9fkoZa+9RnEz471e1pY5nRS9mRlGw5JJ3uN2U8I8HtaqaW13acbj/Nxs2vliMdaVeTyb+7cTBGFdEEETBEEwmGiZ6cjf28voVaPBXZWjo3zNzVGMHA6mMBcWWOsViQAHDzIatrxMwZuasn3NPB67E3N5mQX/ra08Q9PlYhoyHKakORx8hs/HWjXT2b9aZR1bvc4oWVcXo2aCIOwoRNAEQbiz0do2jq3VKE8HD1KY0mngxReBkREKlRlTrVK+FhcpUj09wHvfy6jYzAzw0ksUt0KB84TDlKxKhXPOz7P5bGsrr7ndfG4iwdRotcqUaH8/v1td5bhq1fYwC4elrkwQdjAiaIIg3JnUapQs00+sv5/yU6kwSnb5MttjlEpMXTablDSzGzOZBI4cYSozk2GEbWrKNpNtaeGuyZYWilwqxRTkhz/MZ2htTwJYXeVzlAIefpjyZa41GhQ/k9qUujJBuCMQQRME4c4im+XB4o0Go1b33sso2MwMW2OMjrKYf2GBhfz5PN/Pz1Oo+vuB972P98zOMsI2P0/pcrlsF39zvFK9Dtx3H2vECgWO6+ujhOVylLB9+5hObTR4zUTmEgnO5/Vu9l9NEITbjAiaIAg7n3qdAra4yIiWaU2RywFnzvAszLk5RtWWllj7lUpRwMplStJ730thKhbtiQC5HEUqEGCa0+mklM3OUuSOHqWQra4yYrZ/PyVsdZVRuwcfpHzlcrxHKXuQejC42X81QRA2ERE0QRB2Ltks047VKuu2Dh2yBf+Tk8ClS5SybJYiNj/P6NrSEovz+/t5JqZSlLvnn2eas1ajWCUSjJYVCrwnEuEzEglKVy5HcYvFGIlzOFjf1tXFz0tLNh26axeL/qWuTBAEiKAJgrDTqNeB6WnKTyhkG78uLADHj7O2LJ2mZJVKNiI2Pc3IVlsb8Au/QMkqldhKY37eCpbfz1oyh4PXFhfZs+zIEaZAMxnOMzjIz80mJe6RR6zojY7a45aSSaZGBUEQrkL+qyAIws7A1JbV62xlcegQRWl8nPJlzsNcXWUUbX6e515mMpSuoSG+AKY3L1zgnI0GU5fRKF/VKiNmwSA7/Le32zSm388IWaNBOdyzhxK2tEQxNLs19+xhhE4QBOEGiKAJgrB9aTSYwjS1ZcPDlKCpKba6GBtjmrFW49hKhccrjY7yfU8P8MEP2pTk2bMUtkKBUma6/AcCFLBymfLX1UVZW1lhBK6nh3LmdPJMzAMH+MypKda4mc0ILS2b/RcTBGGbIIImCML2I5tlRKxepxw98ABlaWSEUbTpaUa6lLJHMpneZH4/d00ODlKoJiaA8+cpWs0mX8EgI2BK2Qja0BDlyzSQrdcZCQOsHLa1Mfp25gwFbmiIGwwcjs38awmCsA0RQRMEYXtwbbRs3z7K1JUrjIiNjVGcGg17gLk5jimXY6rxAx+gZC0tUaLyeUa66nVKlMdDMXM4WBcWCjFilkgwglatUsK6u/l9Tw8lrVDgpoPlZbbQuOsuRt8EQRDeJSJogiBsbUy0rNGgGN1/PyNiZ89S2Ex3f7+fgjUxAZw8ySiaz8co1vAw55qa4kHljQZlq15nlCwUYhrT7bb1Zp2dvN80le3utt/t28d7Ll8GTp/m2IcekiOXBEFYN0TQBEHYelwbLdu7l6I0PQ2cO0cxyucZxYpEGBF77jm2zVhd5c7Io0cZ4cpmeU+5zHlNKtPcGwxSzEzbDHPmJUDpi8cpXkNDbLsxN8fnRyKMlMXj0hpDEIR1RwRNEIStQzbLXZeNBlOFHR2s6Tp9mqnM2VlGveJxStiZM8A3v8nUotfLHZWDg+zSPzvLSJrDwRRnucxIWDhM4TI9x4JBylYsRlHzeDjO5+NmgH37+MwrV5gyHRqypw8IgiBsECJogiBsLqZx7PKyjZZls4xUnTnDlGShYOvBUing6acZxWo2KVbDwxQnE3kz3+Vytv1FNMqxPh9fHg8lLxzmc/N53hMOs31GLMZ1jY5S/N7/fjlySRCE24YImiAIm0M2a0Wqt5fpxYUFRssuX6agKcXar2gUeOEF4JVXKFKBAFORbW28N5PhTkyl2Fx2ZYVRsJYWilZrK6XM6aSstbXxeQB3drpcFMN9+/h5dpbPOXjQjhMEQbiNiKAJgnD7aDZZxG+iZUNDFLVUinJ28SIjXrEYZenSJeBrX+NGAJfLRrv6+/lzZobNZpXiPKUSI2OxGKUukWDUq9m09WXJJO9bWGAK9Zd/mTI3M8Mdn8PDrF+TujJBEDYRETRBEDYeU1vWbLKuKxTicUvnzrF32cICI1v9/dyR+b3vcZdmqcQasWiUMtXTw92XFy9y3mqV82jNqFpXFyUsFuOzvF5+7u5mHdnEBAXwwAEeATUxwdqyZBL48IclhSkIwpZBBE0QhI2h0aAApdOMdvX0MG04N8do1cgIx3V0APfdx9Tml7/MDl/G+AAAIABJREFUnZsuF9ORnZ0UqUiEc50+zXsKBe7GdLkob+Ew5zE9zEIhpjHDYT5vbIzy9olPcA3pNHeEHjpEORMEQdhiiKAJgrC+rKzYvmWmqWsmw2L706dZuN/Swp2QS0vA979PWatUuIvS76eUDQ3ZMyxNzzJTyO/zsXA/meTLtMoIh/lSihKYTgP33MN1mBMGhoZ4GLrswhQEYQsjgiYIwq1Tr3PHYzrNVGMyyVRiJsN05Pg4I2JGvF5+GfjDP2QtmtYUqtZW4MgR1pCdOQP86Eecu1KhoJkdloGAPQ8zEKDQAZS0fJ4S1toKvO99lMFSiev7yEcYWRMEQdgGiKAJgvDuyWQYLTPtLjo6KEWXLrF4v1plcf6RI4yEPfkkf1YqjLB5vdwMcOQII2wvvWS/q9X4DJeLEbdEglLW30/RMlE1l4tiVixS/vbtYxSvXGaUrqdHCv4FQdh23JKgKaXiAP4WwACACQD/XGuduc64zwL492sff09r/Zdr1x8A8BcA/AC+DeC3tNZaKfU7AH4TwOLaPf+71vrbt7JWQRDWiVqN0bJMhtGueJzXcjlGvmZnmYLcs4cCdeIE8I1vcHy1ykhaIgF86EOMgD3/PPDXf81IGkCZqtcpbx4PU6TDw8DAAL+fmWENmlmL2w3s3k2pU4oS94EPSMG/IAjbGqXNfxTfzc1K/QcAaa31l5RS/w5ATGv9b68ZEwdwAsBhABrASQAPaK0zSqmXAfwWgBdBQft/tdbfWRO0gtb6P76T9Rw+fFifOHHiXf8+giC8Bek067hqNUa0/H6mDycn2WEfoBwlk4ySvfgif5ZKlKdAgBGto0dZtP/CC4x8AfZgcaUoaqb57OHDfH/lCov9fT47NhajwAFMed59N8VPomWCIGwTlFIntdaHr/fdraY4PwngA2vv/xLAcwD+7TVjPgrge1rr9NpivgfgY0qp5wCEtdYvrF3/KwDHAHznFtckCMJ6YVpTpNOMfMXjlK1sFvjxj5lKDIV4JmW5DJw6xdTm4iJFzulkivGRRzjupz8F/vN/5ryBAAXL7baRtYEB7qw8fJjPOHOGDWiDQdsw1uu152geOMCUpxE1QRCEHcKtClq71noWALTWs0qptuuM6QYwddXn6bVr3Wvvr71u+IJS6tfB6Nv/er3UqSAIG4SJlpVKLMxva2Nt2Guvsejf7WbUqrubEbRvfpPXs1mKVjTKaNm993KDwDPPMALmclGmQiFGysplzv/AA8Cjj+LJUghPnFhC6r8vowtVPBZN4FhfnM8GKGW7dzN9as7SFARB2IG8raAppb4PoOM6X/0fN/mM6/0XVL/FdQD4MoDfXfv8uwD+bwD/0w3W93kAnweAvr6+m1ySIAhvolqllC0t8XN7O4UolQJ++EPuyoxGeSRSpcJI2WuvUbwaDaY89+4FHnyQ70+eBP70Tylhbjejbx4Pi/lXVxkt++AHgYceAopFPPm91/D4ShJl5QQAzMCLx1eSgDONYwc6GaXr6pJomSAIdwRvK2ha60dv9J1Sal4p1bkWPesEsHCdYdOwaVAA6AFTodNr76++nlp75vxVz/gvAJ5+i/V9BcBXANagvc2vIwjC1WjNVhfT0yzyj8e5E7NcZjpybo5C1NnJn6kU8Pd/zxqyYpHi1dZGedq7l9GyH/2I8ynFdGR3N9Od5TLneP/7gV/8RQrahQs8yqlYxBPBD/5czgxlOPBErRvHPv5hiZYJgnBHcaspzqcAfBbAl9Z+/vfrjHkGwP+llIqtff4IgMe11mmlVF4p9R4ALwH4dQD/CQCM9K2N/xUAZ29xnYIgXE2lAkxNccel00kpi8XYMPbsWcpQOMxC/WqVh5S/+qqNroXDlLJ77uH9Z88Cf/mXTIkC7EPm8TBSViyyDu3RR4H3vIf1ZN/9LvC3f8tI2913Aw4HUueu/5+jVK4iciYIwh3HrQralwD8nVLqcwCuAPhnAKCUOgzgf9Fa/8aaiP0ugONr93zRbBgA8K9g22x8B3aDwH9QSt0HpjgnAPzPt7hOQRBMtGxigtGytjZGsdJp4Ac/YP2Yz0eZ8njYzuLrX2fas1KhWHV1MVLW3s7zM597jhsCmk3KVns768xKJUrV0aOMlu3bx/H/7b9x3r4+NpJNp3m9tRVdASdmSs03Lbsr6r/dfylBEIRN55babGw1pM2GIFyH1VVGy6anKV5dXRSq06fZUNblYtF+LEap+slPWFuWzXJ8NMoUZ28vv7940R46bo5mCgR4dFKlwrEPP8xXVxcl7vvfZ53avfdyvvl5fh4eZgStrw9PnpnH4986g3Kt8fOl+91O/P4/OYhj93ff8NcTBEHYrmxkmw1BELYizSbTkePj7DXW0cEo1uXLwLe/zXown4/S5fNRuJ5+muKkNfucDQ3x+1CITWZfeokROJeLYpZIsMZMKaY5h4aAD3+YElavM/p24QKF7MgRjlta4vxHjjBFGov9fMlGwp54ZgSplTK6on489tG9ImeCINyRSARNEHYS5TJla3qaKcmeHtaAvfoqrzkcFKZ4nBL33e/ag8qDQXukUlsbJWtujpJnOvb7fIyqBQIUrWCQLTI+9CE+6/XXeZzT8jLP3RweZqQtnWZk7f772SbDNJsVBEG4g3mrCJoImiBsdxoN1oGNjVGGurpYxH/uHJu8VquMdHV0ULBefZVd/E2TWZ+PotXaSnEqFCh52ayNlgWDVsxqNT7jfe9jjZnLBXznO5wTAO67j4JnzsMcHqbEtV2vTaIgCMKdi6Q4BWEnUiyyhmxujpGvwUGmKF98kYX3LhePXUokmFp8+mmmOD0eCldvL8UtvtYIdmGBtWpK8d62Nv4MBPg8p5Np0kceYTpzaQn46lc5ZyLB9hnRKOdYWaGo3Xuv9C0TBEF4F4igCcJ2otFgL7KxMaYo+/ooS2fPAj/7GUXLnGOpNdtjHD/OyFosxp5kTicjZwAjXGfO8HuHg9eDQUbgmk3OEYkwAvbww/z+1Cng936PO0H37QP+6T9lOnR6mtG2T3yC6U5BEAThXSOCJgjbgVwOGB1lLVciQQEbH+fuyHSakbDOTqYrUyn2GDORtWSSsuV0MhpWrXLMygqjZR4P50wmKWnlMp85OMhTAe6+m2nNp5+mzDmdLPLfu5ctMyYmKGq/8Rs22iYIgiDcEiJogrBVqdfZg2x8nKlGk5J8/XXge9+zBf979rBu7ORJts4A2I+st5cy5lr717xYpFCZgv9g0J4ckMvx+3iczWcfeIB1ZpOTwJe/zIa2HR2MlgUCbLUxOgq89708sFwayQqCIKwrsklAELYaS0usLcvlKFrBID+fPctUZCjEgv56nZGw48fZBiMW4yuX4zw+H6NhmQylzBy95Pcz1RkK2SOZBgcpegcOML35ox8BL7/MlOn99zNilk5TFtvb2U4jkdjcv5MgCMI2RzYJCMJWp1JhXdn8vN1xCQAnTjB6ZdKQpkv/T39KaXM4gP5+ylahwJfWlLSlJaY2fT6+YjGeHJDLcUNArcZC/qEhYP9+pjy//W1GxuJxe17m6Chr2Q4eBD7/ea5PEARB2FAkgiYIm4XWTDmao5Q6OhgVu3SJDV6VYoG+18vNAdPTjKLlchzb0kKhazY5JptlhE1rCpnpW9bVReGamuL37e3sRTY4yE0Gr77KKFw6zVqyo0dZZ3b+POd43/s4XtKYgiAI64pE0ARhK5HPszVFOk3JiscpX888w+hYKMRoWbNJcTt1it/7fIxoLS1RxopFpjCLRcqTz0ehczo578AA51tY4Dx9fdxIcPAghe/FF4FvfIM1ZYcOcTPA9DT7p3V3A5/6lKQxBUEQNgmJoAnC7aBe527HhQW+N41cR0ZYR+b1Mo3pcPA1Pc1C/FLJnoE5P0/R0pqSp7U9B9Pc29nJCNnEBJ+TSHAjQV8fa8ymppg2XV7mtfvuY+pzbIyid999rDmTNKYgCMKGIxE0QdgsFheZwiwWGd0KBFho//LLjJB5PIyYuVyMhl28SGELhylmi4tMgwJMbVarlLJYjBEz07tsYIDpy5kZylxnJyXLFP2fOwf81//K5+zdC/zSL/F54+PcRPDQQ0x5Op2b+ucSBEEQiAiaIKw35jzMbJYCFQhQtEz9mN/P6JfXS2GanWUES2tGv7q6OP7cOUbMKhWOi0YpZeaA8mSSInbpEiNj3d08/9Lno5g1Gizuz2SYRn3/+zm/mbu7G/iVX2GUTerLBEEQthQiaIKwHpgO/2Z3ZCjEaNf4OGvGlGLEzOFglEopbgRYXmZkzUjZxASjbcUiC/TjcYqYaQDr81GyXC5KnVIs4M/nOc++fYyi/eQnFL5du9hs1u1minRsjKnOj3xEmsoKgiBsYUTQBOHdojWjZFNTjJq5XLw2OWkL8wGblgwEKHGnT3NsOMx+ZktLvCeb5f2mHYbpWdZsUr7a21mblkpRvO6/nxsNwmGK2fnz7F8WifC7lhbeu7TE5z3wANOYbvem/tkEQRCEt0cETRDeKabWa2WFETGleKzS3Bz7kAEUNqeTsqQUI1e5nO3en80yWjY/z+hbKMR+Zl4vxQrg3K2tvHbxIsfddZft+t/Xx+japUtcT28vBczvpxwuL1P2Hn3UHuMkCIIgbAtE0AThZmg0mIJcWOB7p5OCNjlJKavXKWVaU7BiMXuoOcDUpN/PerCZGe7O9Hp5qLjfz1coRLHy+9nnbGGBY3ft4kHlqRTFrrubKc2xMdalHThgd3SWy1xXf7/dICAIgiBsO0TQBOFGmHYWc3O2Jiyf567MTIapy2qVshUIsHWGOQ6pUKDEeb38fnSUgudwMJo1OGgjZFpT8EIhitmlS5TA/fttirJaZcQsnWaas6eH8wBMX1YqlLP9+yl00iZDEARhWyOCJgjXsrpKmUqnbfuL6WmmI02krFDgd62tjI4tLzOaVqlQ5BwOGwGrVhlRu/tuezpAJELZazYpZek0I2S9vcB73sMoWC7HCFgyaXeE7t5td3W6XJzb62Wj2Y4OaZMhCIKwQxBBEwSA4pXJULTqdUa1UilGz8plSluhwB2ara1MM2YyjHYZ0XI6OW5sjELl9TLq1dLC7zo6OC6bpah1d/NEgclJ9ibbtYvPqlTYH63RYORuYIBji0VKms9nj2zas4c1bdImQxAEYUchgibcuZhDxTMZSpbbzchZKmVFLZejNAUCVpIyGYqVESZzpqY5F7OjgzVgTidrxOJxzlcqMQ1qNhXU6+zcX6/bEwQaDY6Lxfi8UonCZnZ0Nho83HxggBsOBEEQhB2JCJpw51EuU7JM2nBlxe7CrNUoXoWC7cjvdPLzxYsUttVVClM6zTRmoUAR272bdWReL/uaVav2APPOTqZJJycZVRsc5Pc+HyUsm6WADQ3xWj5vd306HJzzrrsobR7PZv8FBUEQhA3mlgRNKRUH8LcABgBMAPjnWuvMdcZ9FsC/X/v4e/9/e+ce3NZ5nvnnJQmCBG8AQYJ3ipR1v1iWzDhOfKnrOFGSjRN3N02TTrbeaTNuupdp2hk1cbMzze42jRt3drs7nek04zRJdzveZtPEqZ20rl3HYye+yrJl3a+UxPsNIAmQIAkQ3/7x4OTQMi1KoWRA0vObwRDn4DvnfEQU6ef3fb/3c859J3/+qwB+A0DEOVe9ZHwQwN8AuBnABIBfc86dWc1cxXWOl8JMpShouRz7l8XjFLS5OQqRl8Ls6uI1Q0M876Ul02mKnCd3sRjTk6WljJRFIpS2eJyfAYyupdN4rH0XHp6OYHCqFK2zOeypHMV9U0Ps5N/Tw2dPT/u7D+Ryfk+0pia1yRBCiOuIVW2WbmZfBxB3zj1kZl8CReuL542pB7AXQA8AB+A1ADc75xJmdiuAswBOnCdo/x7Ajc65z5vZpwH8inPu11aajzZLF2/BS2EmkxSzXI5ylUj4jWRTKT8iFosxOjUxQcFKJChkZkx9Tk7yGq9pbE0Nr2tuZupxfJxiFQ5T/ubmuNqysxOPJQJ4cCyMtPNrxSrN4Wtby3Ff7RzHhkJ8VmkpI3CdnbyXEEKIa5ILbZa+WkE7BuAu59yQmbUAeNY5t/G8MZ/Jj/nt/PFf5cc9umRM6jxBexLAV5xzL5pZGYBhAI1uhclK0AQAP4WZTvuRs/FxvrxomRdFa2ykcM3NUd6mpyljZWUUO68+raKCstTSws8iEUa+hoeZqoxGKW8DA5S2zZs5PpMBAgHcdqwOA9m3r7BsK8/hZ7cYxay6mvfv6GC6UwghxDXNhQRttTVoTc65IQDIS1psmTFtAPqWHPfnz12In1/jnMua2RSAKIDxVc5XXKtks4x6eZuLT08zSub9nJmhhC0ssHVFeztFaniYn4+NccziIiVudpYCV1HBbZSqq/1eZ9ksr8lmKWoDA8CxY4x63XUX55PJMBKWX1QwmI0sO+3BhRKgsYHXevInhBDiumfFfw3M7GkAzct89OWLfMZy6/9XCttd9DVm9gCABwCgs7PzIqckrgm8vTBnZvz6LU+0RkcpWl4UrbSUUa5IhNf09/Pn+DilbGGB18/P+w1km5oocXV1PO7vZ01aQwOFrb+f6dBNm1j0PzHBSJgZnzkyQqHbsAGtz6UxMPf2X6G1ppybmatNhhBCiCWsKGjOuXve6TMzGzGzliUpztFlhvUDuGvJcTuAZ1d4bD+ADgD9+RRnHYD4O8zvGwC+ATDFucJ9xbXA7CzlK5ejjI2NMQ05NfXzgnzkcnxVVTGNacb2GYODlCxvWyZvcYC3VdKaNbwmFKKILS5yvBct6+8HDh3iasrdu3nfZJKfe6svzbhgYPNm3j+RwJ7WUjx4JoB0zv81KgOl2PPRLZIzIYQQb2O1+ZR/AHA/gIfyP3+4zJgnAfyJmXk5ng8BePAi7/sigE8CeGal+jNxjZPJsB4sk2G0yyva9zYuj8cpU4uLFKVolKnM8XFuSj42xvfeH6NkkpJWXs4i/2jUrzPz9tEcHOT5UIhF/2Vl3A1g926/JUcoxPPxOO8VizFdmUxyXpWVQCSC+zqqgTVleHhfAoNTc2gNV2LP7o24b+dK2X4hhBDXI6tdJBAF8F0AnQDOAfhV51zczHoAfN4597n8uN8E8If5y77qnPtW/vzXAfw6gFYAgwAecc59xcwqAPxvADvByNmnnXOnV5qPFglcYzj31nqweNzfl3J8nNGsbJYRqFyOQtbQ4G/VlEhwjFcPlskwygZwBaaXqqyq4vtsluNranwpm5xk0f7OnRSxRMIfOz3N+5aWUuoaG/nc6WmKXm2tv/G5tmESQghxHldsFWexIUG7Rjg/hTky4m+zdOYMPzOjwAUCrBWrqPDrzvr6KElm/Hxqyt/QvL6e8lRZ6b8fGOBzYjFePzDA+23bxsax4+N+DVs8TgGsruZ82tooeGfPMnrX2OgLXns7768UphBCiGW4kqs4hbg8ZDKUn2zWT2GO5xftDg8zsgX4m4TX1jJq5Qnc4CAFbXHR3wLJayYbibCnWGUlxaqpidcND/NneTmfd+AAo2Uf/Siv81Z1RqOcy9gYxauigvK1sEBhrKriCkxvQUFHB8cJIYQQvyASNFE4vI3DZ2cpaKkUpWku37j11Cl+VlFBYSovZ3rRjCsm33yTYuZtLh4K8XwqRRHr7OS5ykpGx0IhFvyfOkXpMvMFa8sWjllY4NyiUUboJib4uRcVa27mHI8c4fhNmxhda8i3yggECvqVCiGEuDaQoIl3n5kZitnCAsVsdJQiVFLC9OTYmC9k5eWUrbo6itepU/7+l85RnObnWRvmbbdUVUUpC4eZckwmKXLBIO+bSjEi1t0NbN3qb0JeU8N7xuMcHwpR1OrrKV69vbyus5ObnHt7dTY2ahsmIYQQlxUJmnh38FKYXiPZmRmmLc1YL3bmjL/dUVkZhai+nuI0NgYcPcpx2SylraKC90ul3tp41qv9CgT8BrK1tbxuYoLvt2+n9HlRt9pa3isep2jV1vptNiYmuEl6VRWFzpPGjg5KoxBCCHEFkKCJK0cux6L7mRm/N1l/P9OWAHDuHCNSwSClzBOsmhped/Cg31qjtJTSNT9PkQoEmGIMBChMTU18TUywYL+8nFI2O0ux6+4GNmzwV21GozyfTjPt6a3EjMU4l95e3qe5mdEygELW1kY5FEIIIa4gEjRx+fFSmKmUXy82MsLolLfKcnGRMhYIcEwkwnNnzvgrJUtK+AoEKFO5HNOWDQ0UqpoaphtLSih+Bw9SzObn+fxIhFLm9UXzUpUjI3yZ8VpvY/KREeDwYc5j3Tpg/Xp/ZWZzs7ZhEkII8a6hf3HE5WFhwd9YPJvl69w5Rr8Avz1GZSVFJ5v1Vz0ODgJvvMFolrdNUnm5v39mMEhBMuP5lhYej40BJ0/6uwYsLPDalhaKWCTCqFpjIwv7vRSmGQVv/XpK25tvAk8/TeHbsoXPM6O4eYsJhBBCiHcRCZr4xfFSmN5+mMEgI1leD7LJSUpaSYkvPdkshcmrO5uZYWG+t3rSOb9GrbaWAldWRnm64QaK1alTlDqAES4vmuZF5NragK4uXxLn5vy+adXVwI038twrrzC619kJvO99/uboHR0cJ4QQQhQICZq4dFIpik0ySfGanGQd1+IiZezMGdZ+lZdTjGZnKVDBIMVqYICylM1SxrzUoZfWrK2laJWVUZ5iMa7cPHLE30MzGGQ0rrKSz4hGmZZsamL92IkTlKzKSj6jq4srNk+fBn70I0b2Nm8Gduzg57W1FLvy8kJ+s0IIIQQACZq4WLztlRIJPxo1OEhJKy/nZ6Oj/nZGqRQjWuEwx54+TSEDKGLO+QKXTvsRsJISRs26uyloJ0/613oNar22GBUVHOdtSn78OMWxpoYLAZxjtKy1FXjpJeCb36S07dzJZ2UyTHHGYtqGSQghRFEhQRPvjLf/5cQEo0xetGx8nJ8BLPjPZBjx8jrve9sdTU2xTgzwt2qqqnrr2KoqSldJCaNcDQ185v79foSuro5jFhZ4TSTCVhnt7Wy/8dpr/uKBqSle89738jlPPw384AeMxN1zD6XNjNfW1am+TAghRFEiQRNvx+von0pRdtJpHs/OMno1Pk4RKiujRE1NUXYiEV4zNESBS6UoVSUlFLbZWYpZZaWf8qytpZh50bKDB3ltXR2L+L39OLNZFvBv38577t/PVKrXXiOV4n3uvJOC98QTXJW5fTtw773+YoOODj5fCCGEKGIkaILMz1No4nG/JmxqimnL8nJK0eioXzc2Nua3x2hspCAlkxSh2VlGuurqmIacnvZXcIbDvF97O1OLiQQjYNPTHNvWxjHxOJ/f0gLccQfl69gx4Pnned/OTt4zk6GEtbQwxfntb3N+t9zCKFo6zSjdunXahkkIIcRVgwTteiaXo2iNjPB9IEDh8dpjBAIULy8dOTRE+QqHGbWam+Nn8/OUJa+3WVMTjxMJf4PyQMDfSLy8nKnJffuYcozFgI0bee94nNdt28b6sdJSYO9eyllLC/e+9NKmt9/OqNi+fUxjVlcD738/n7W4yJRnNKptmIQQQlx1SNCuR6am2A7DS/uZMTo2PU25yWYpSuk005kTE36z12iU101MUMIWFnhNSwvPj4/7qdBw2N96KRzm+Jdf5vPDYcpWOMy5nD1Lebv3XkbLensZLfMWAniRuepq4OabOdcXXqC4tbUBH/6wL5UdHUyhCiGEEFcpErTrhbk5itDkJFOY5eUUsHPn+Hkw6AtWPM5WGNkso1sbNlB+kklG22ZmGNnyiveHh3mfigoel5VR5Jqb+ZzDh4EXX2Qkq7mZrS1mZ7kKNJXi8bZtlLlXXwUOHGAUbtcuiuPsLHugtbbyd/jhDxnNW7cO+NjH+HsEgxwTDBb2exZCCCEuA+acK/QcLhs9PT1u7969hZ5G8bC4SNEaH6dsVVdTwkZGKEbV1TyfSFDKenspYd4WSsEgj8fHOd6rK4tGeext3+T1IqutpYDV1fF+Bw8y0hWJMApWX89I2cwMo2o330ypGhkBXn+dadYbbuC8xsf5c9MmzuP0aaYy5+Z4rrWVqVVvkYDaZAghhLjKMLPXnHM9y32mCNq1yMQE218sLFCeKiooQf39FJmKCorVuXNcCTk4yKhXdzdXSs7NMXKVSPC9t0l4VZXfBLaiwl+JGY1yoYAZa8teeIH3a2sDenoYAevrY/TuppsoWPX1lLLHHmOac+dOztdbHXrHHRS5I0d4z2CQc/NSl83NFD+1yRBCCHENogjatcLMDOVpdtZva5FIsKB+fp4SlMlQ3s6epWQtLDAStXUrx8TjFDlvF4D6ehbaj40xjeltlVRa6kfLKir4nGPH/Nqy7m6mRnt7GUFrb2e6cu1aRuT27mUErrvbX405N8f3zc2c45EjFEdvBWZZGd+3t/OnEEIIcZWjCNq1SiZD2Uok/AaupaU8Pn2aglVTw9ThoUOUnokJSlRPDyUrHmeN2PQ0BayxEVizhpJ3/DgjX6GQv/1Sc7N/z95eplBLSih6O3b4uwZMTlL8tm7lPQ8fBh5/nEK3bZufxozHmdasrGRd2ZNPUuIaGxltKylhhK61VdswCSGEuG5QBO1qI5djZGloiFGn+npKUSJB4SkvZ0rSazZ74AAlrrSU/cI2bKA89fVR1tJpylBDA687c4bpTYDH3s+WFt5jcpLRt2TS70fW3Mz5xOOMnO3axZWYuRx7nI2OUrA2beJch4f97Zy8th6nTnF8ayvlrbycz4zF1CZDCCHENcmFImgStKuFeJwik0wyolVTw3PxOFOVTU2Moo2Osjj/yBGmKru72bTVjNcPDjIdWlND+amv5/HJk/6m5uXlTCm2t7P2a2GB9Wujo75EtbRQpE6c4L03bGC0rKWFkvfmm5SxLVvY9mJykq+mJgrdzAylbGCAUTVvP8xIxG9WK4QQQlzDSNCuVmZnKTuJBKNIsRjrvBIJilp9vV8DdvIkVzmOjDAleOed/HxU1sIvAAAXTUlEQVR42E85entQxmJMVx4/zvoyr2WGc35tmXN8xrlzlKnycoqZJ1sjIxy7cyclMBjEY0+/iYd7HQZzZWgNlWHPe5twX/Uspc5bZDA15derRaN+urS1lWO0DZMQQojrBAna1UQmw/Tj6CgL91taGMGKxylGVVVMScbjjIa9/DIlzoxbG23dSrE6fpxj5uY4vr2dkbexMUauFhYYAauspPy1tPBnLkep8wr3o1E/6nX6NOfU1eVvrzQ6Cuzfj8fGS/DgdBPSOf9XqSwFvnZ3B+7bFGU69fhx/n7t7f6m5d3dvHeZyiGFEEJcX1wxQTOzegB/B6ALwBkAn3LOJZYZdz+A/5w//GPn3Hfy578K4DcARJxz1UvG/zsADwMYyJ/6C+fcIyvN56oVtFzOb4ORSlGKvE3JEwlGuJqbKVUjI1wFefAgI2wbNwK33caI2NGjrAXzUpidnf4+madOUZJKSxlZy+UoaNEo55BKUbampjimsZHylMtRGMvLufXS2rW87vBhRu1qa4Ht23HbP01gIJl526/WVlmCn23I78PZ2uqnaNetY92b2mQIIYS4TrmSgvZ1AHHn3ENm9iVQtL543ph6AHsB9ABwAF4DcLNzLmFmtwI4C+DEMoLW45z7j5cyn6tK0LzC/nPnKEWhEGVpcpIilcvxuKKCEa0jRyhmo6N+CtPbEqm31xe5zk6+slnWd/X1MRJWW8vo2+Kin+LM5fisRIKLBUIhRsW6u/1dB7zVmS0tXOl54AAjc17PtIoKYHAQ3d8fw3J/kgxA7yeb+JyGBgpldfUyI4UQQojriyvZZuMTAO7Kv/8OgGcBfPG8MbsBPOWci+cn8xSADwN41Dn3Uv7cKqdxFTEz46cwva2PysoYGfNaYHR38/jwYW595KUw3/Me4FOf4j0OH2YhfiZD4br1VgpWPM4GsPG4v+WSc35fs9JSXu+1uMjl+Mxt2yhwZ88yFbllC7B+PSNxJ08C//iPvN/27cDdd1PoxsZ4LhJBa/UkBlJvj6C1BsHo2V13aRsmIYQQ4iJZraA1OeeGAMA5N2RmsWXGtAHoW3Lcnz+3Ev/GzO4EcBzA7znn+la6oGiZn/dbY8zPU8paWnjspQk7OxnpOnECeOIJSlIqxVTgZz/LMYcOAc8847e42LSJhfXJJKWvr4+Rs3CY5xcWeF0oxHlMTzNalkpR+FpbGdFKJhkxC4cpeh0dnMvhwxTF9nbgnnv8LZy8lGcsxgjaxAT2NM7gwZkA0s6X7coyw56Pbwe2dhToixdCCCGuTlYUNDN7GkDzMh99+SKfsVx4bKW86uNghG3ezD4PRufufof5PQDgAQDo7Oy8yCm9C2SzjDD191OMolG/qP7MGUarmpsZRTt5ktsjHTzIyFY4DPzSL1GezpwB9u+nWAWDvMf73sfrBgZ4XTzOzxobmcI0Y+QrFPI3QE8kKIehEKNg7e0s+j92jI1i772XAjYwQAkEgM2bueXS4qK/SCEY9Pe+HB7m3MNh3Hf3NuDVs3j4+DwG54DWcAX27N6E+3ZejIsLIYQQYimrrUE7BuCufPSsBcCzzrmN5435TH7Mb+eP/yo/7tElY1JLa9DOu74UrHOrW2k+V7IG7bHXB/Dwk8cwOJlGa7gSe3ZvfLt85HIUoYEByllNDaUpmWS0zGu+WlnJMQcOMCp27hyFZ/NmrsRMp1nwPzHBFGZTE9OekQjvOzTEa7y2GDU1/ubn4TAlanKSc5ma4rxiMa7w9HYICIV43NnJ+3gLDBobKXD19bzn9DSv9xrgmnHBwdAQI21NTYyoBYO8XyRyRb5/IYQQ4lrjSi4SeBjAxJJFAvXOuT84b0w9uDBgV/7UPnCRQHzJmNR5iwRavNSpmf0KgC86525daT5XStAee30AD37/ANKZxZ+fqwyU4mv/ejvuuym/MnF4mNIVCFBanGP0bHGRqcTaWsrSG2+w4P/0aYpYWxulrKmJwjYy4vc4a2/nQoD5edaGDQxQukIhfp7J+DVkS/fanJzkCs9AgCnSTZsYiRsb41y2beP4sTGKWSbDcRs3UhSTSdapOccUZiTipzxnZni/igqKWTjMejXtjymEEEJcEldS0KIAvgugE8A5AL/qnIubWQ+AzzvnPpcf95sA/jB/2Vedc9/Kn/86gF8H0ApgEMAjzrmvmNnXAHwcQBZAHMDvOOeOrjSfKyVotz30DAYm028731YdwM/uClGS6usZwerroxzFYkxhzs4ydXn4MNOJySQjaLt2UWwGBrgKM5lkwX1TE2WpspLS19dH0Ssp8aNliQSfFYkwquWtxJyaovQ1NXHlZU0NBWxxkfK1dq0fARsY4LU33sj0aybDOWQyTM/W1vLlbVweDHK3gIUFpmmbmniswn8hhBDiF0KNaldJ95d+9M4tJH5nGyNRY2P+CszFRUrQ/v1MJw4NMZrV1cXO+wAXAyQSFJ5YjOnClhZKkhct82QuFmPkKpPhM5qb/W2dvFSmFy3bsYNRuL4+Cti2bWxvMT1NWZufZ1pz82YK4fw8pW5+nhGzcJjyd/YsxTEW4/ipKaY8W1v5O6qxrBBCCLEqrmSbjeuC1nDlshG01iCYOly7lvLV3892FIcOUcq8XmYf+ADl6+xZtsCYmaEIrV3LF0Ch+tnPfJnzGs2OjTG92NDAqNXIiN/3LJ3m/X/5l5kqPXKE2z2tWQPs3s10ZX8/z4dCjJY1NvJ5MzOUvNlZylZjI38ePsyo2Zo1wM03c3GAt9n5TTfxnkIIIYS4okjQLoI9uze+vQat1LBn9wagPssGsm++SVnLZFiftXUrpS2d5vmTJ5mm9PazDIcZlTp0iNEyr8i/q4vnMxlGtHbt8ncaeP11illZGcf19Pibjo+P+5uVJ5O879wcn3fPPSzyX1zk+GyWP0MhilgqxXsvLjIK5+23OTXFOrjmZomZEEII8S6iFOdF8pZVnNUB7OkG7hs9xLoybzPxWIzCU1vLVZaJBOWqsZGvNWuYSuzro5QNDVHmqqpYy+XVlnktObzo1dgYha26mlGstWspZcPDHOcV7Q8O8rlVVTwXy7elW1igfHmpzEiENXODg5x/TQ3TlgsLlLqKCl7b2EipFEIIIcRlRzVol5PvfpcblI+OUmRqaxlh8va8HB6m6Hibmnd1cVw8ztTk4CAlqLKSkjQ1xft6zWrLyvwVockkBam9nSs9S0uZrpyfp4C1t1O6Tp+mJDY383wgwOvSaf6cnOQzYjFGzY4f5/29hrkzM/zcm1M0KjETQgghrjCqQbuc7N/vpwK9vTIHByk9ZWX+6s3GRtZ3nTlD4RoZobRVVDBiNT3Nz1taKGbT036qcnGRIvWe97Dof2CAz62r43FVFQXx1VcpVTfcwFRmLsc5LiwwPRqP8/M1ayiFR49SIjs7GYlLpd66YjMa1eblQgghRBEgQbtUvFWaU1MUL+eYMuzo8Ju+jo9TngYHOTYYpLSl0xSiYJCiVVvLIv7nnqNUlZVxMcAtt1CWDh/mZ93drCNLp5keTaW4YOD22yl8ZkyBAoyGzc5SAjdv5hxfeokRsc5ORtfm5/1nVVUxaiYxE0IIIYoGCdqlMjfH2rHqaq6c7Ozk+2SSEarhYdaMeXVlwaC/YjIWo2xls2xKOz5OUQoGGZG75Ra/wL+3l33SduzguDfe4LjubkbLFheZ8lxY4M94nHLY2MjUZ28vt2yqqWG/skyGklZdzXNVVVyoIDETQgghig4J2qVSX08JisUoSSMjjHR50bJAgJ8vLFCISkq4urKpyZemXI6yFY2yL9oNN1DunnuO0vf+91PiBgZY9N/YyP03Kyv5jFyOr8VFtsQoK6MoLu2/1tDAHmgzM/w8EvEXGoTDhf4WhRBCCHEBJGiXyoYNlLFDhyhniQTTjCUlFJ90mpGs+npg/Xpes3cv8OKLfuSqoYHpyZIS9i3r7+d9162jcB09Sqnq7mZq1FvIkc3ymtlZpjm91ZoTE4ywLSywpm3zZo7xiv7Tac6xqalw35sQQgghLhoJ2qWydy93AcjlmB6sr2fkyjkK0rp1XLl58iTw9NNMP1ZWUsq8zdCPHQNeeIHX3ngj7zswwFRmNMrFAVVVrBUrKWHkq7SUqzGzWd6ruZnp1J/+lDLX1saxi4u8b0MDxSwY5HshhBBCXDVI0C6VyUl/7835eQpRVRXTmIEA8PzzwCuvcMVlTQ1TinfcwfGvvgo8/jgja7ffzrqxs2cpel1djH4BlLBczq8xm5rimLY2ft7Xx/09q6sphF7Rf1sbr52dZTQvGi3Y1ySEEEKIXxwJ2qVSXc304uwsRau7mynJJ5+kJNXVsd5r7VpuwXTuHGvCKiooUxs2cJHBqVNcxblzJ0Vubs5PZQYCfEY6zehbVxefd+IEBbGhgfdJp9kgt6uLtWbT034TWiGEEEJctUjQLpXKSuC22yhTzz/PprXRKMWsvJwF/mvXMor2+OMs3u/poVgNDzMq1tHBaJeZ33qjtJT39HqThcNciDA97deXxWKUs/l5yt26dbzv+LjfYFYIIYQQVz0StEulrAx44glGuSIRvlpbgQ9+kMX6r7/OSNf69ZSw4WGmMWtr2Ry2tpayNT/P+rCKCr6fnaWwedsrTUyw+a23f6dHUxNlMB7nIgWv3kwIIYQQ1wwStEsllaKUlZQwPbljBxcO/PjHjGBt2cJ04/g4BaytjX3JvLTl/DyjZeXlPPY2V29uZj3b4CD7qFVUMBLnrcBsbWX0bnyc0uftYiCEEEKIaw4J2qVSV8ei/1yO6c3eXorUjTdSrEZGWFO2Y4ffe2xmxt/mKZOhmOVyHFdezpRmby9bdoTDbxWz7m4K3fg4Fws0NPAaIYQQQlyzSNAulZYW4NlnKVzr1vnF+ZkMo1xtbYx0eeerqvjyNiQHmOYsKWE07dQppjdjMRb7Z7N+2jKXo5jlchIzIYQQ4jpCgnaplJZyBeXkpF8DtmULBcs5Rrmco0yVlFDSSkp4HAjw+mSSfc+co/BFIn6bjLo6StrICJ/X0MDrhBBCCHHdIEG7VJJJClRbG+vGamqYjkwk/H0uZ2b83mRVVSz+LynhmMFBRtja29laIxTy68syGbbg8BYLlOl/HiGEEOJ6RAZwqWzcSKkqLWW0LJHgsbdhuhcpKyvzdxsYGWGqMhIB1qzx22SsX89xCwsUt5ISRuJKSwv9WwohhBCigEjQLpXaWi4GqKxkdCydZhozFGL0bHGRopXLseN/Ok3pWrOGKc2mJqZFzShqIyOUtKYmiZkQQgghAEjQLh0z1olNTTElGQqxncbcHCNn8/PcPSCXY+pycZGrMdvbGWUDOHZiguNbWih0QgghhBB5JGiXSiLhN6nNZChbXv3Y6dMUtrY2nq+upqQFg7w2naaYBYMSMyGEEEK8I6syBDOrN7OnzOxE/mfkHcbdnx9zwszuz58LmdmPzOyomR0ys4eWjA+a2d+Z2Ukze9nMulYzz8tKNMqoWDpNUZua4sblySS3dfL24tyxgz3MgkG20ejv5+KBtjZ/twAhhBBCiGVYrSV8CcC/OOfWA/iX/PFbMLN6AH8E4L0AbgHwR0tE7s+cc5sA7ARwm5l9JH/+twAknHPrAPwPAH+6ynlePhYXGTEbGgIOHGD9WEcHa9NaWylmLS2sJ5uZoZjNzVHMGhqYIhVCCCGEuACrTXF+AsBd+fffAfAsgC+eN2Y3gKecc3EAMLOnAHzYOfcogJ8AgHNuwcz2AWhfct+v5N9/D8BfmJk559wq57t6vMayzc2UMi+NGQr5Y5JJRtaqqvxN0YUQQgghLpLVClqTc24IAJxzQ2YWW2ZMG4C+Jcf9+XM/x8zCAO4F8D/Pv8Y5lzWzKQBRAOOrnO/qaW2lfIXDlLSlTWSnp/mqruaiACGEEEKIX4AVBc3MngbQvMxHX77IZywXPvp5JMzMygA8CuB/OedOX8w1583vAQAPAEBnZ+dFTmkVRKNsmbG0hmxqimJWWysxE0IIIcSqWVHQnHP3vNNnZjZiZi356FkLgNFlhvXDT4MCTGM+u+T4GwBOOOf+/LxrOgD05wWuDkD8Heb3jfw90NPTc+VToLW1/vtEgnVmdXWsQxNCCCGEuAysdpHAPwC4P//+fgA/XGbMkwA+ZGaR/OKAD+XPwcz+GJSvL1zgvp8E8ExR1J8BbDYbj7P4PxBgxKymptCzEkIIIcQ1xGoF7SEAHzSzEwA+mD+GmfWY2SMAkF8c8N8AvJp//VfnXNzM2sE06RYA+8zsDTP7XP6+3wQQNbOTAH4fy6wOLRhjY2ydsbTxrBBCCCHEZcSKJTB1Oejp6XF79+4t9DSEEEIIIVbEzF5zzvUs95m6pQohhBBCFBkSNCGEEEKIIkOCJoQQQghRZEjQhBBCCCGKDAmaEEIIIUSRIUETQgghhCgyJGhCCCGEEEWGBE0IIYQQosiQoAkhhBBCFBkSNCGEEEKIIkOCJoQQQghRZFxTe3Ga2RiAs1f4MQ0Axq/wM6529B1dGH0/K6Pv6MLo+1kZfUcXRt/Pyrwb39Ea51zjch9cU4L2bmBme99pY1NB9B1dGH0/K6Pv6MLo+1kZfUcXRt/PyhT6O1KKUwghhBCiyJCgCSGEEEIUGRK0S+cbhZ7AVYC+owuj72dl9B1dGH0/K6Pv6MLo+1mZgn5HqkETQgghhCgyFEETQgghhCgyJGgXiZn9tZmNmtnBQs+lGDGzDjP7iZkdMbNDZva7hZ5TsWFmFWb2ipntz39H/6XQcypGzKzUzF43sycKPZdixMzOmNkBM3vDzPYWej7FhpmFzex7ZnY0//fR+wo9p2LCzDbm/+x4r2kz+0Kh51VMmNnv5f+OPmhmj5pZRUHmoRTnxWFmdwJIAfgb59y2Qs+n2DCzFgAtzrl9ZlYD4DUA9znnDhd4akWDmRmAKudcyswCAH4K4Hedcy8VeGpFhZn9PoAeALXOuY8Vej7FhpmdAdDjnFMPq2Uws+8AeN4594iZlQMIOecmCz2vYsTMSgEMAHivc+5K9xC9KjCzNvDv5i3OubSZfRfAj51z336356II2kXinHsOQLzQ8yhWnHNDzrl9+fdJAEcAtBV2VsWFI6n8YSD/0n8hLcHM2gH8KwCPFHou4urDzGoB3AngmwDgnFuQnF2QDwA4JTl7G2UAKs2sDEAIwGAhJiFBE5cdM+sCsBPAy4WdSfGRT9+9AWAUwFPOOX1Hb+XPAfwBgFyhJ1LEOAD/bGavmdkDhZ5MkbEWwBiAb+XT5I+YWVWhJ1XEfBrAo4WeRDHhnBsA8GcAzgEYAjDlnPvnQsxFgiYuK2ZWDeDvAXzBOTdd6PkUG865RefcTQDaAdxiZkqX5zGzjwEYdc69Vui5FDm3Oed2AfgIgP+QL78QpAzALgB/6ZzbCWAGwJcKO6XiJJ/+/TiA/1fouRQTZhYB8AkA3QBaAVSZ2WcLMRcJmrhs5Ouq/h7A3zrnvl/o+RQz+bTLswA+XOCpFBO3Afh4vsbq/wK428z+T2GnVHw45wbzP0cB/ADALYWdUVHRD6B/SWT6e6CwibfzEQD7nHMjhZ5IkXEPgF7n3JhzLgPg+wDeX4iJSNDEZSFfAP9NAEecc/+90PMpRsys0czC+feV4F8ERws7q+LBOfegc67dOdcFpl6ecc4V5L9cixUzq8ovwkE+dfchAFpZnsc5Nwygz8w25k99AIAWKi3PZ6D05nKcA3CrmYXy/659AKypfteRoF0kZvYogBcBbDSzfjP7rULPqci4DcC/BaMe3vLtjxZ6UkVGC4CfmNmbAF4Fa9DUSkJcCk0Afmpm+wG8AuBHzrl/KvCcio3/BOBv8/8/uwnAnxR4PkWHmYUAfBCMDokl5KOv3wOwD8AB0JMKsqOA2mwIIYQQQhQZiqAJIYQQQhQZEjQhhBBCiCJDgiaEEEIIUWRI0IQQQgghigwJmhBCCCFEkSFBE0IIIYQoMiRoQgghhBBFhgRNCCGEEKLI+P+Cf8XXOIFmHQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Compute the slope and intercept of the fra diff/lane curve\n", "slope, intercept = np.polyfit(lanes, f_13, 1)\n", "\n", "# Compute bootstrap replicates\n", "bs_reps_slope, bs_reps_int = dcst.draw_bs_pairs_linreg(lanes, f_13, 10000)\n", "\n", "# Compute 95% confidence interval of slope\n", "conf_int = np.percentile(bs_reps_slope, [2.5, 97.5])\n", "\n", "# Print slope and confidence interval\n", "print(\"\"\"\n", "slope: {0:.5f} per lane\n", "95% conf int: [{1:.5f}, {2:.5f}] per lane\"\"\".format(slope, *conf_int))\n", "\n", "# x-values for plotting regression lines\n", "x = np.array([1, 8])\n", "\n", "# Plot 100 bootstrap replicate lines\n", "for i in range(100):\n", " _ = plt.plot(x, bs_reps_slope[i] * x + bs_reps_int[i], \n", " color='red', alpha=0.2, linewidth=0.5)\n", " \n", "# Update the plot\n", "_ = plt.plot(lanes, f_13, marker='.', markersize=12, linestyle='none')\n", "plt.draw()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The slope is a fractional difference of about 0.4% per lane. This is quite a substantial difference at this elite level of swimming where races can be decided by tiny differences." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Hypothesis test: can this be by chance?\n", "The EDA and linear regression analysis is pretty conclusive. Nonetheless, you will top off the analysis of the zigzag effect by testing the hypothesis that lane assignment has nothing to do with the mean fractional difference between even and odd lanes using a permutation test. You will use the Pearson correlation coefficient, which you can compute with `dcst.pearson_r()` as the test statistic. " ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p = 0.0\n" ] } ], "source": [ "# Compute observed correlation: rho\n", "rho = dcst.pearson_r(lanes, f_13)\n", "\n", "# Initialize permutation reps: perm_reps_rho\n", "perm_reps_rho = np.empty(10000)\n", "\n", "# Make permutation reps\n", "for i in range(10000):\n", " # Scramble the lanes array: scrambled_lanes\n", " scrambled_lanes = np.random.permutation(lanes)\n", " \n", " # Compute the Pearson correlation coefficient\n", " perm_reps_rho[i] = dcst.pearson_r(scrambled_lanes, f_13)\n", " \n", "# Compute and print p-value\n", "p_val = np.sum(perm_reps_rho >= rho) / 10000\n", "print('p =', p_val)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The p-value is very small, as you would expect from the confidence interval of the last exercise." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }