content/11_statistics/rank_based_methods.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rank based methods\n",
    "\n",
    "Consider a case where you hypothesized that your lizard would grow over time, but you had no idea if that growth \n",
    "would go up at a constant rate (linear) or whether it would slow down over time (perhaps logarithmically) or if lizards might have growth spurts at various points.\n",
    "\n",
    "In such a case, we might want a method that could account for any of these possibilities: we just suppose that the lizards keep getting bigger, and never get smaller, but beyond that we don't have any special feelings about whether the increase is linear, logarithmic, or goes up in fits and spurts.\n",
    "\n",
    "In such a case we might choose to use a __rank based__ method. Instead of predicting the exact length of a lizard based on it's age, we would instead merely assert that the youngest lizard is smallest, the oldest lizard is largest, and the ones in the middle sort out according to age. In that case we could replace the y value for the smallest lizard with a 1, the y value for the second largest lizard with a 2, etc. These would then be the *ranks*  of the lizards lengths.\n",
    "\n",
    "If we were doing a normal linear regression of length on age, we would be doing Pearson regression. If we instead regress the *rank* of age on the *rank* of length then we are performing Spearman regression.\n",
    "\n",
    "It turns out that by transforming interval or ratio values into __ranks__ we can produce a range of new methods that make fewer assumptions about our distribution.\n",
    "\n",
    "For each common *parameteric* statistical test from statistics there is likely an equivalent non-parametric or rank-based alternative.\n",
    "\n",
    "| Parametric Test  | Non-parametric Test          |\n",
    "|------------------|------------------------------|\n",
    "|Pearson regression|Spearman regression           |\n",
    "|T-test            |Mann-Whitney U                |\n",
    "|Paired T-test     |Mann-Whitney Signed Rank test |\n",
    "|ANOVA             |Kruskall-Wallace              |\n",
    "|Principal Coordinates Analysis (PCoA) | Non-metric multidimensional scaling (NMDS)|\n",
    "\n",
    "\n",
    "If you can understand the difference between Pearson and Spearman regression, it will be much easier to understand all the methods in the table above. \n",
    "\n",
    "Let's consider an example dataset in which each x step increases y exponentially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1a25ea15c0>"
      ]
     },
     "execution_count": 206,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucXHV9//HXZ3ZndzbZSyDXJSRsVAgJJCQhASyCaMBiS0UQiy0IlGrUWiot9fLjUapia2u12lZrMQXFS6pVQMRWKTTlYopQkoCSEO4kIWQ2N5LZ3ezM7lw+vz/O2clms5fZsLMzs/N+PshjZ875zpzPLMn3M+d7NXdHREQEIFLqAEREpHwoKYiISJ6SgoiI5CkpiIhInpKCiIjkKSmIiEiekoKIiOQpKYiISJ6SgoiI5NWWOoDRmjZtmre1tZU6DBGRirJhw4a97j59pHIVlxTa2tpYv359qcMQEakoZratkHJqPhIRkTwlBRERyVNSEBGRPCUFERHJK1pSMLOYmf2fmf3KzDab2WcHKVNvZv9uZi+Y2WNm1laseEREZGTFHH3UA7zd3bvMLAqsM7Ofu/uj/cr8IbDf3d9kZu8DvgBcXsSYRESGtCWe4N5Nu3j1QJLZUxq48NSZLGhtKXVY4xpX0ZKCB1u6dYVPo+Gfgdu8XQx8Jnx8B/A1MzPXdnAiY0KV3OhiWv3wy7Q0RGltiZFIpln98MusOndeSWMb77iKOk/BzGqADcCbgH9298cGFJkNvALg7hkzSwBTgb3FjEtkrE3kSs7dcQ++0bk7OQcnOBacP/TcgVxYngHH3R0Hnm3vYM2j22mK1dIUqyWeSPKPa5/n98+Yy5tmNOXL9r33oTjCn/2v3b/sYeUPL9M/zvzjAe+35rFtZHI5etJZetJZANLZHN9+ZBu/d8bcQ3EM+N0ceezIkod9jsGODfM+P9rwCt29GcyCEi0NdQDcu2lX5SUFd88CS8xsCvBjMzvV3Tf1K2KDvWzgATNbBawCmDt37hEvkOpRLpWvu5PJOZmss3nnAb71v1tprK+lKVbDzgPd/MP9z3P5ijnMm94YVpJBZZrNeb7SzLmHzznsfN/jnDvZvtfmgue5fudyDrncgGO5Q4/XbtnFwd4M+7pqgpiBVDrL5//zGd5y4rR8JXp4Bd8X2+HHx9IvX9pHTzpLLFqTP5ZKZ/naAy/y5jdMHfsLFuhXryRorK9htx2qltydl3qyzGzeXbK4nol30lhfw4HuND6lgZaGOppitbx6IFmU643LjGZ3P2BmDwIXAv2Twg5gDrDDzGqBFuC1QV6/GlgNsHz5cjUtValCv/m6O+msk8nlgp/ZHJmck87myPQ7ns4Of/7Q64PHvZlcmAiC432GquRuefilo6rkzCBiRk3E8o8j+Z/h44gd/tyMSMSorQ2OmUEynWXq5DoiZkBwrCVWy/7uNG+a0YgZGEb4X/51Rt9PwnOH4ug7buFJG+K1fZ9jsNc+097BtMZGaiKHV757Onu5ZOnsQ3Hlr9Pv/cID1u9c/7L54jb4cTPr97jf6w16Mlk6kmlaGqL5a3akMjQ3RFl17hsOi6P/dQfGOdixQT9Pv5gGHhtYNjEgrs5UhtlTGo686BgoWlIws+lAOkwIDcD5BB3J/d0DXA38ErgM+B/1J0ifTDbHwZ4snT1punoyfHPdVjqTaVLpLDsPBN+QD/Zk+MLPn+UtJ047rIIfrYgZtTVGtMaojUSCnzURojURJtdHqI1EDzsfPI5QGzE270wwffpkaiMRIhb8Y3aHPV09vO+MOfkKscYOVeQWOfS8fxKI2OGVxOvx0p6DR1QmiWSaN85oYuWCmWNyjaNx8qzmQeOaP6uJtmmTSxbXxUuOY/XDL1NXm6MpVktnKkN3b5Yrz5rL5PrSrQh04akzWf3wywD5uBLJNJevOL4o1yvmJ20Fvh32K0SAH7r7f5jZzcB6d78HuA34rpm9QHCH8L4ixiNlJJtzulKZfIXflcrQmcrQGT7u6klzsCd72Gs27+ygpaGWeq+hJqxEm2O1JFIZ5h47KaikB1TqtRHLH49GIkRrBz/f/1vraJ1yXMugldyJM5pobSnOt7lCjHdlUulxLWhtYdW58w5rnrx8xfEl7xsa77is0r6YL1++3LUgXvG9nrb7bM6Dir4nQ2cqHVb+hyr+rp403b3ZI9qr66MRmupraYzV0lQfpTFWG7bT19IUi/LNdS/R1ZM9ovJtaYjypxecNJYff1T6N2v1r+RKPWqlL7Zy6IOplLgmMjPb4O7LRyynpCADDVfJnTSzOV/hB5V8ut+3++D5YBV+XW2E5lhtWNFH+1X2QcXfGKulvrZm8IAKiKvUFYoqOSl3Sgpy1L5y/3Ps6+qhN5ujJxN0sCa609REjNPbjhm0wj9UwR9dhV8oVb4iR6fQpFBx+ylIcb12sJfHt75GNpcDgnb7utoIk+trONib5cx5Uw+r7JvGsMIvxILWFiUBkSJSUhDcnZ2JFOu3vsZLew7iOaexvpY3zWjKD7Psa7t/8xtLN45cRIpPSaGK5XLOS3u7WL91P/FEioa6Gs58w7GcN38a3/nldnoyOepqI2UzOkREik9JoQqlszm2xDvYuG0/+7uDO4C3nTyDha3N1NUGC+fGojVlNzRPRIpPSaGKJHuz/GrHAX71ygG6e7PMbI7x24un8abpjUQGjNNX271IdVJSqAKJ7jQbt+9n884E6awzb9pkTj/hGI4/pmHMZs+KyMSgpDCB7epIsWHbfp7b1UnEjPmzmjj9hGOY1lhf6tBEpEwpKUww7s7Wfd1s2LafV17rpq42wuknHMOSOVNoikVHfgMRqWpKChNENuc80x50Hu/t6qUpVsu5J03jlONaDlu9U0RkOEoKFa4nk2XTqwme2H6AzlSGaY11vOOUmZw8q/l1LfImItVJSaFCdabSPPnKAX69I0FvJsecYyexcsFM2qZOUuexiBw1JYUKs7erhw3b9vNseyc5d06aGXQez2yOlTo0EZkAlBQqgLuzY3+SDdv28/Leg0RrjEWzW1g29xhaJqnzWETGjpJCGcvlnBf2BMtQ7OpIMamuhje/cSqnHT+Fhjp1HovI2FNSKEO9mRxPh8tQJJJppkyKsnLBDBa0NhOtiZQ6PBGZwJQUSmjg3gDnzZ9GMp3j1zsSJHuztLbEOPekabxh2pHLUIiIFIOSQon030XsmElRNr16gPs2t7N07hTOfMNUlrcdy3EtMY0kEpFxpaRQIvdu2kVLQ5TaiPHUjgQYHDclRmOslouXzC51eCJSpZQUSuTVA0laW2LEDyRxYOmcKURrIsQTqVKHJiJVTL2WJTJ7SgOd4Wb3DdEa6mtr6ExlmD2lodShiUgVU1IokQtPnUkimWZPVw+T62tIJNMkkmkuPHVmqUMTkSqmpFAiC1pb+L0Vc6iJREilc7Q0RFl17jxtbCMiJaU+hRJqaojy5jdM5cqzTmB6k/Y4EJHS051CCbUnUtTVRpg6ua7UoYiIAEoKJRVPpJjZHNPENBEpG0oKJZLO5tjT2UNri1Y3FZHyoaRQIrs6UuTcmaWkICJlREmhRNrDSWq6UxCRcqKkUCLxRIqWhiiT6jQATETKR9GSgpnNMbMHzGyLmW02s48NUqbFzH5qZr8Ky/xBseIpJ+5OeyKluwQRKTvF/JqaAW5w941m1gRsMLP73f3pfmU+Cjzt7r9jZtOBZ81sjbv3FjGukuvsCZa3aNWSFiJSZop2p+DucXffGD7uBLYAA5f/dKDJgvWhG4HXCJLJhKb+BBEpV+PSoG1mbcBS4LEBp74G3APsBJqAy909Nx4xlVI8kaI2Ykxr1CxmESkvRe9oNrNG4E7genfvGHD6N4EngeOAJcDXzKx5kPdYZWbrzWz9nj17ih1y0bUnksxsjlGjSWsiUmaKmhTMLEqQENa4+12DFPkD4C4PvAC8DJw8sJC7r3b35e6+fPr06cUMuegy2Ry7Ono0P0FEylIxRx8ZcBuwxd2/PESx7cDKsPxMYD7wUrFiKgd7unrI5lz9CSJSlorZp3A28H7gKTN7Mjx2IzAXwN1vAT4H3G5mTwEGfNLd9xYxppLr21lNdwoiUo6KlhTcfR1BRT9cmZ3AO4oVQzlqT6RoitXSFIuWOhQRkSNoRvM4iydStLZofoKIlCclhXF0sCdDRzKtpiMRKVtKCuMorklrIlLmlBTGUXsiRU3EmKGtN0WkTCkpjKN4Isn0pnpqa/RrF5HypNppnORyzq6OlPoTRKSsKSmMk71dPaSzmrQmIuVNSWGc5DuZmzUcVUTKl5LCOIknUkyur6G5QTutiUj5GrGGMrMZBEtWHAckgU3A+mpY4nostSeSzGppIFgSSkSkPA2ZFMzsbcCngGOBJ4DdQAx4N/BGM7sD+PtBlsOWAZK9WfZ3pzlldkupQxERGdZwdwq/BXzQ3bcPPGFmtcBFwAUES2PLMNo7wkXwmtXJLCLlbcik4O4fH+ZcBri7KBFNQPFEEjOYqaQgImVuuOYjA95LsI/yHcDbgYuBZ4Bb1KdQuPiBFNMa66mrVb++iJS34ZqP/hmYAdQRJIN64KcEzUrzgY8VPboJIJdz2jtSnDyrqdShiIiMaLikcI67Lwq31GwHWt2918z+jaDjWQrwWncvvZmcZjKLSEUYrj0jA+DuaeBxd+8Nn2eA7DjENiG051dG1aQ1ESl/wyWFdjNrBHD3C/sOmtksoLfYgU0U8USKWLSGYyZppzURKX/DjT565xCnOgmGo0oB2hNJWltimrQmIhVh2OEwZvax8Od1fcfc/aC77y52YBNBKp1l38Fe9SeISMUYaYxkl5l9HDg4HsFMNLs7enDXTmsiUjmGTApm9mmCoaefA+ab2V+OW1QThCatiUilGTIpuPtnCTqULwDS7n7zuEU1QcQTKY6dXEcsWlPqUERECjLSKqkPufsvzEybCo+SuxNPpHjj9MmlDkVEpGAj9SksDH8uKHYgE82B7jSpdFbzE0SkoqijuUjyO61NUX+CiFQOdTQXSXtHkrraCMdOqit1KCIiBVNHc5HEEylmNceIRDRpTUQqx0jNRw+7+y+AB8chlgmjN5Njb2ev5ieISMUZNim4+31m9ofu/t99x8ysJmxakiHs7kyRc9dMZhGpOIXs+rLSzH5mZq1mdirwKDDi5gBmNsfMHjCzLWa2uW/JjEHKnWdmT4ZlHhpl/GVJK6OKSKUaaZ4C7v77ZnY58BTQDfyeu/9vAe+dAW5w941m1gRsMLP73f3pvgJmNgX4OnChu283sxlH9zHKy85EiimTojTUadKaiFSWEe8UzOxEgl3W7gS2Au83s0kjvc7d4+6+MXzcCWwBZg8o9vvAXe6+PSxX8QvtuXt+ZVQRkUpTSPPRT4Gb3P1DwFuB54HHR3MRM2sDlgKPDTh1EnCMmT1oZhvM7KrRvG856khlONiTZZaajkSkAg3ZfGRmre4eB85w9w4Ad3fg783snkIvEG7Ucydwfd/7DLj+6cBKoAH4pZk96u7PDXiPVcAqgLlz5xZ66ZLo6084TncKIlKBhrtT+KaZPQrcGHYG5xOIuz9fyJuH+zvfCaxx97sGKbIDuDfco2Ev8DBw2sBC7r7a3Ze7+/Lp06cXcumSiSeSRGuMqY1aLkpEKs9wk9feCZxHMEfhEuBRM7vLzFaZ2Yhf1y3Yauw2YIu7f3mIYj8BzjGz2rCf4kyCvoeK1Z5IMaM5Ro0mrYlIBRp29JG7p4B7wz+Y2TzgncDXzGyWu58xzMvPBt4PPGVmT4bHbgTmhu99i7tvMbN7gV8DOeBWd9/0ej5QKWWyOXZ39rB07pRShyIiclRGHJLan7u/TDCE9OtmNuyiPu6+Dhjx67K7fxH44mjiKFd7unrI5lwjj0SkYhUy+ugIZvaUu/eOdTCVbueBoJNZI49EpFINN/ro0qFOAbOKE05la0+kaIrV0lg/qhswEZGyMVzt9e/AGsAHOaf2kUHEE0ktbSEiFW24pPBr4EuDdfya2fnFC6kydfVk6ExlWHaC8qWIVK7h+hSuBwZONutzSRFiqWjtiSSAOplFpKINeacQ7qMw1Ln1xQmncsUTKWoixnRNWhORCjaq0UdmtrFYgVS6eCLFjKZ6amuOakCXiEhZGG0Npmm6g8jmnN0dKW2qIyIVb7RJ4T+LEkWF29fVQzrrGnkkIhVvtElhqDWMqtrORN+kNd0piEhlGzIpmNnZ/bbSPNPM7gfWm9krZvbmcYyx7LUnkkyur6E5pklrIlLZhqvFvgL8LtBI0Gz0bndfZ2bLgK8SLHgnBJ3MrS0NBAvDiohUruGSQtTdnwIwsz3hAneEey6r8TyU7M1yoDvNotktpQ5FROR1G65Pof+5/zfg3LArpFaTeDhpTf0JIjIRDJcUbgo3vsHd7+47aGZvBL5T7MAqRXsiRcSMGU1KCiJS+Yab0TzoPszu/iLwd0WLqMLEEymmNdVRV6tJayJS+YYbfbTazBYNcW6ymV1rZlcUL7Tyl8s57R0prXckIhPGcB3NXydoQloEbAL2ECyZfSLQDHyTYGntqrXvYC+9mRyzmtXvLiITw3DNR08Cv2tmjcByoBVIAlvc/dlxiq+stYeT1nSnICITxYizrdy9C3iw+KFUnngiSUNdDVMmRUsdiojImFDv6OvQ15+gSWsiMlEoKRylVDrLvq5eZjWr6UhEJo7R7qcQMbPmYgVTSXZ19PUnqJNZRCaOEZOCmf2bmTWb2WTgaeBZM/t48UMrb/FECjOY2aKd1kRk4ijkTmGhu3cA7wZ+BswF3l/UqCpAeyLF1Ml11NfWlDoUEZExU0hSiJpZlCAp/MTd04AXN6zy5u7EEylmqelIRCaYQpLCN4CtwGTgYTM7AegoZlDlbn93mlQ6q/kJIjLhFDJP4Z+Af+p3aJuZva14IZU/rYwqIhNVIR3NM83sNjP7efh8IXB10SMrY+2JFHW1EaZO1griIjKxFNJ8dDvwX8Bx4fPngOuLFVAlCHZa06Q1EZl4CkkK09z9h0AOwN0zQLaoUZWx3kyOvV09ajoSkQmpkKRw0MymEo44MrOzgMRILzKzOWb2gJltMbPNZvaxYcquMLOsmV1WcOQlsqsjhbsmrYnIxDRiRzPwZ8A9wBvN7H+B6UAhlXcGuCHc07kJ2GBm97v70/0LmVkN8AWCJqqy1x7OZNbyFiIyERUy+mijmb0VmA8Y8Gw4V2Gk18WBePi408y2ALMJZkX3dx1wJ7BilLGXRDyR4phJURrqNGlNRCaeEZOCmV014NAyM8PdC96n2czagKXAYwOOzwYuAd5OBSQFdyd+IMkJUyeXOhQRkaIopPmof2UdA1YCG4GCkkK4Sc+dwPXhchn9/QPwSXfPDjeSx8xWAasA5s6dW8hli6IjmaG7V5PWRGTiKqT56Lr+z82sBfhuIW8eLo9xJ7DG3e8apMhy4AdhQpgG/JaZZdz97gExrAZWAyxfvrxkS2zEO4JJa0oKIjJRFXKnMFA3wT7Nw7Kgpr+NYPvOLw9Wxt3n9St/O/AfAxNCOYknUkRrjGmNWhlVRCamQvoUfsqhBfAiwELghwW899kEq6k+ZWZPhsduJFhlFXe/ZdTRllh7IsXM5hiRiCaticjEVMidwpf6Pc4A29x9x0gvcvd1BKOVCuLu1xRathQy2Rx7OntYNveYUociIlI0hfQpPDQegZS73Z09ZHOumcwiMqENmRTMrJPB900wwN29qrbljCf6tt9UUhCRiWvIpODuTeMZSLlrT6Robogyuf5o+uZFRCpDwTWcmc0gmKcAgLtvL0pEZSqeSHLcFK13JCITWyH7KbzLzJ4HXgYeItiF7edFjqusdKbSdKYy6k8QkQmvkFVSPwecBTwXzitYCfxvUaMqM+1hf8JxWhlVRCa4QpJC2t33AREzi7j7A8CSIsdVVuKJFLURY3qTJq2JyMRWSJ/CgXD9ooeBNWa2m2C+QtVoT6SY0VxPjSaticgEV8idwsUES1v8KXAv8CLwO8UMqpxkc86ujhSz1HQkIlWgkDuFVcCPwlnM3y5yPGVnb1cPmZxrfoKIVIVC7hSagf8ys1+Y2UfNbGaxgyonfZPWNPJIRKrBiEnB3T/r7qcAHwWOAx4ys/8uemRlIn4gSWN9LU2atCYiVaCQO4U+u4F2YB8wozjhlJ94IsWslhjDbQIkIjJRFDJ57SNm9iCwlmAjnA+6++JiB1YOunszJJJpjpuipiMRqQ6FtImcQLCV5pMjlpxgDvUnaOSRiFSHQpbO/tR4BFKO2hMpImbM0KQ1EakSo+lTqDrxRIrpTfVEa/RrEpHqoNpuCLlw0prmJ4hINSmko/mPzazq9qDcd7CX3kxO8xNEpKoUcqcwC3jczH5oZhdalYzNbNdOayJShQqZvPYXwInAbcA1wPNm9nkze2ORYyupnYkkDXU1tDRESx2KiMi4KahPwd2dYOJaO8EKqccAd5jZ3xUxtpJqTwT9CVVyYyQiAhTWp/AnZrYB+DuCzXUWuftHgNOB9xQ5vpJIpbO8drCXVs1PEJEqU8jktWnApe6+rf9Bd8+Z2UXFCau01J8gItWqkMlrfznMuS1jG055iCdSmMGMZk1aE5HqonkKg2jvSDK1sZ762ppShyIiMq6UFAZwd+KJFK3NajoSkeqjpDDA/u40PWlNWhOR6qSkMMDOA0lAncwiUp2UFAZoT6Soj0Y4dnJdqUMRERl3SgoDxDs0aU1EqlfRkoKZzTGzB8xsi5ltNrOPDVLmCjP7dfjnETM7rVjxFKInk2VfVw+zmjVpTUSqUzF3o88AN7j7RjNrAjaY2f3u/nS/Mi8Db3X3/Wb2TmA1cGYRYxrW7o4e3NWfICLVq2hJwd3jQDx83GlmW4DZwNP9yjzS7yWPAscXK55CHNp+U0lBRKrTuPQpmFkbsBR4bJhifwj8fIjXrzKz9Wa2fs+ePWMfYCieSHLs5DpiUU1aE5HqVPSkYGaNwJ3A9e7eMUSZtxEkhU8Odt7dV7v7cndfPn369KLE6e60J1K6SxCRqlbMPgXMLEqQENa4+11DlFkM3Aq80933FTOe4XQkM3T3ZtWfICJVrZijj4xgY54t7v7lIcrMBe4C3u/uzxUrlkLsTAST1nSnICLVrJh3CmcD7weeMrMnw2M3AnMB3P0W4C+BqcDXw3kBGXdfXsSYhtSeSFFXG2HaZK2MKiLVq5ijj9YBw84Ac/cPAB8oVgyjEU+kmNkcIxLRpDURqV6a0Qykszn2dPaoP0FEqp6SArC7s4ecu/oTRKTqKSkA7QmtjCoiAkoKQNCf0NIQZVJdUUfoioiUPSUFgpFHuksQEVFSoDOVpjOVUX+CiAhKCvlF8FpbtFy2iIiSQiJFbcSY3qRJayIiVZ8U2hNJZjbHqNGkNRGR6k4K2Zyzu6NH/QkiIqGqTgp7OnvI5Fwjj0REQlWdFOJaGVVE5DBVnRTaEymaYrU0xaKlDkVEpCxUdVKIa6c1EZHDVG1SONiTIZFMqz9BRKSfqk0KmrQmInKkqk0K7YkUNZq0JiJymKpNCvFEkulN9URrqvZXICJyhKqsEXM5Z1eHOplFRAaqyqSw92AP6awmrYmIDFSVSaG9r5O5WZ3MIiL9VeVWY/FEikl1NTQ3VOXHFyGdTrNjxw5SqVSpQ5ExFovFOP7444lGj25SblXWiu3hpDUzrYwq1WnHjh00NTXR1tamfwcTiLuzb98+duzYwbx5847qPaqu+SiVzvLawV7NT5CqlkqlmDp1qhLCBGNmTJ069XXdAVZdUjg0aU2dzFLdlBAmptf7/7UKk0ISM5jZrKQgIoX5/Oc/f9jz3/iN3xjxNY2NjSOW+cxnPsOXvvSlYcvcfffdPP300yO+11ipuqTQnkgxrbGeutqq++giR21LPMFX7n+OP//Rr/jK/c+xJZ4odUjjamBSeOSRR8bt2koKReTutHek1HQkMgpb4glWP/xyfgHJRDLN6odfft2J4Xvf+x5nnHEGS5Ys4UMf+hDZbJZt27Zx4oknsnfvXnK5HOeccw733XcfW7du5eSTT+bqq69m8eLFXHbZZXR3dwOwdu1ali5dyqJFi7j22mvp6ekBoK2tjU9/+tMsW7aMRYsW8cwzzwBw8OBBrr32WlasWMHSpUv5yU9+AsDtt9/OpZdeyoUXXsiJJ57IJz7xCQA+9alPkUwmWbJkCVdccQVw6C6gq6uLlStX5q/R917D+eu//mvmz5/P+eefz7PPPps//q//+q+sWLGC0047jfe85z10d3fzyCOPcM899/Dxj3+cJUuW8OKLLw5abiyZu4/pGxbb8uXLff369Uf12n1dPXznl9t4xykzOeW4ljGOTKRybNmyhQULFgDw4LO72dPZM2TZ+za3k+zN0lBXkz/W9/wdp8wa9DXTm+o5b/6MYa//iU98grvuuotoNMof/dEfcdZZZ3HVVVdx6623cu+993LmmWfywgsv8I1vfIOtW7cyb9481q1bx9lnn821117LwoUL+eM//mNOPPFE1q5dy0knncRVV13FsmXLuP7662lra+OGG27guuuu4+tf/zobN27k1ltv5cYbb2ThwoVceeWVHDhwgDPOOIMnnniCH/3oR9x888088cQT1NfXM3/+fNatW8ecOXNobGykq6srH3/f80wmQ3d3N83Nzezdu5ezzjqL559/HjM74jUAGzZs4JprruGxxx4jk8mwbNkyPvzhD/Pnf/7n7Nu3j6lTpwLwF3/xF8ycOZPrrruOa665hosuuojLLrsMYMhyQ/3/7WNmG9x9+ZD/U0JVdaeglVFFRm9/d5pY9PCqIhaNsL87fdTvuXbtWjZs2MCKFStYsmQJa9eu5aWXXgLgAx/4AJ2dndxyyy2HtbfPmTOHs88+G4Arr7ySdevW8eyzzzJv3jxOOukkAK6++moefvjh/GsuvfRSAE4//XS2bt0KwH333cff/u3fsmTJEs477zxSqRTbt28HYOXKlbS0tBCLxVi4cCHbtm0b9nO4OzfeeCOLFy/m/PPP59VXX2XXrl1Dlv/FL37BJZdcwqRJk2hubuZd73pX/tymTZs455xzWLRoEWvWrGHz5s2Dvkeh5Y5W0eYpmNkc4DvALCAHrHb3fxxQxoBHtzazAAAKlklEQVR/BH4L6AaucfeNxYqpPZEiFq3hmEnaaU2kz3Df6AF27E+SSKZpaTj07yaRTHNKQ5T3Lp9zVNd0d66++mr+5m/+5ohz3d3d7NixAwiaZ5qamoAjR9WYGSO1dNTXB6sg19TUkMlk8te+8847mT9//mFlH3vssXz5ga8Zypo1a9izZw8bNmwgGo3S1tY24nDQoUYHXXPNNdx9992cdtpp3H777Tz44IOvq9zRKuadQga4wd0XAGcBHzWzhQPKvBM4MfyzCviXYgTS10n2zw+8wJOv7OeZ9o5iXEZkQrrw1JkkkmkSyTQ59/zjC0+dedTvuXLlSu644w52794NwGuvvZb/Vv7JT36SK664gptvvpkPfvCD+dds376dX/7ylwB8//vf5y1veQsnn3wyW7du5YUXXgDgu9/9Lm9961uHvfZv/uZv8tWvfjWfUJ544okR441Go6TTR94ZJRIJZsyYQTQa5YEHHhjxzuLcc8/lxz/+Mclkks7OTn7605/mz3V2dtLa2ko6nWbNmjX5401NTXR2do5YbqwULSm4e7zvW7+7dwJbgNkDil0MfMcDjwJTzKx1LOPo6yTbf7CHaI1hMCadZCLVYkFrC6vOnUdLQ5R4IkVLQ5RV585jQevR98stXLiQv/qrv+Id73gHixcv5oILLiAej/PQQw/x+OOP5xNDXV0d3/rWt4I4Fizg29/+NosXL+a1117jIx/5CLFYjG9961u8973vZdGiRUQiET784Q8Pe+2bbrqJdDrN4sWLOfXUU7nppptGjHfVqlUsXrw439Hc54orrmD9+vUsX76cNWvWcPLJJw/7PsuWLePyyy9nyZIlvOc97+Gcc87Jn/vc5z7HmWeeyQUXXHDY+7zvfe/ji1/8IkuXLuXFF18cstxYGZeOZjNrAx4GTnX3jn7H/wP4W3dfFz5fC3zS3YfsSR5tR/NX7n+ORDLI8FviHSxobQKMloYof3rBSUfxaUQq32AdkeVs69atXHTRRWzatKnUoVSEsu5oNrNG4E7g+v4Joe/0IC85IkuZ2SozW29m6/fs2TOq6796IElTrJaIwTGTokyur6UpVsurB5Kjeh8RkWpQ1KRgZlGChLDG3e8apMgOoH9P1fHAzoGF3H21uy939+XTp08fVQyzpzTQmcrQFIsyf1YztZEInakMs6doBJJIpWhra9NdwjgpWlIIRxbdBmxx9y8PUewe4CoLnAUk3D0+lnEUo5NMRGSiKubS2WcD7weeMrMnw2M3AnMB3P0W4GcEw1FfIBiS+gdjHURfJ9m9m3bx6oEks6c0cPmK419XJ5nIRODuWhRvAnq9/cRFSwph5/Gwf+M8iP6jxYqhz4LWFiUBkX5isVh+ZqwSw8TRt59CLHb0S/lU5SY7ItXu+OOPZ8eOHYx24IaUv76d146WkoJIFYpGo0e9M5dMbFW19pGIiAxPSUFERPKUFEREJK/i9lMwsz3A8KtODW0asHcMwxkr5RoXlG9simt0FNfoTMS4TnD3EWf/VlxSeD3MbH0ha3+Mt3KNC8o3NsU1OoprdKo5LjUfiYhInpKCiIjkVVtSWF3qAIZQrnFB+camuEZHcY1O1cZVVX0KIiIyvGq7UxARkWFURVIws2+a2W4zK6sF2c1sjpk9YGZbzGyzmX2s1DEBmFnMzP7PzH4VxvXZUsfUn5nVmNkT4c59ZcHMtprZU2b2pJkVvjVgkZnZFDO7w8yeCf+evbkMYpof/p76/nSY2fWljgvAzP40/Du/ycy+b2ZHv7LcGDKzj4UxbS7276oqmo/M7Fygi2A/6FNLHU+fcD/qVnffaGZNwAbg3e7+dInjMmCyu3eFGyWtAz4W7qNdcmb2Z8ByoNndLyp1PBAkBWC5u5fV2HYz+zbwC3e/1czqgEnufqDUcfUxsxrgVeBMdz/a+UdjFctsgr/rC909aWY/BH7m7reXOK5TgR8AZwC9wL3AR9z9+WJcryruFNz9YeC1UscxkLvH3X1j+LgT2ALMLm1UwZLm7t4VPo2Gf8ri24OZHQ/8NnBrqWMpd2bWDJxLsNkV7t5bTgkhtBJ4sdQJoZ9aoMHMaoFJDLITZAksAB519253zwAPAZcU62JVkRQqgZm1AUuBx0obSSBsonkS2A3c7+5lERfwD8AngFypAxnAgfvMbIOZrSp1MKE3AHuAb4XNbbea2eRSBzXA+4DvlzoIAHd/FfgSsB2IE+wEeV9powJgE3CumU01s0kEG5PNGeE1R01JoQyYWSPBXtbXu3tHqeMBcPesuy8h2Df7jPAWtqTM7CJgt7tvKHUsgzjb3ZcB7wQ+GjZZllotsAz4F3dfChwEPlXakA4Jm7PeBfyo1LEAmNkxwMXAPOA4YLKZXVnaqMDdtwBfAO4naDr6FZAp1vWUFEosbLO/E1jj7neVOp6BwuaGB4ELSxwKBFu8vitsv/8B8HYz+15pQwq4+87w527gxwTtv6W2A9jR7y7vDoIkUS7eCWx0912lDiR0PvCyu+9x9zRwF/AbJY4JAHe/zd2Xufu5BE3hRelPACWFkgo7dG8Dtrj7l0sdTx8zm25mU8LHDQT/WJ4pbVTg7v/P3Y939zaCZof/cfeSf5Mzs8nhQAHC5pl3ENzyl5S7twOvmNn88NBKoKSDGAb4Pcqk6Si0HTjLzCaF/zZXEvTzlZyZzQh/zgUupYi/t6rYec3Mvg+cB0wzsx3Ap939ttJGBQTffN8PPBW23wPc6O4/K2FMAK3At8ORIRHgh+5eNsM/y9BM4MfhXse1wL+5+72lDSnvOmBN2FTzEvAHJY4HgLBt/ALgQ6WOpY+7P2ZmdwAbCZpnnqB8ZjbfaWZTgTTwUXffX6wLVcWQVBERKYyaj0REJE9JQURE8pQUREQkT0lBRETylBRERCRPSUFkjJlZ18ilRMqTkoKIiOQpKUjVMrMVZvbrcP+IyeFa9acOKPMFM/ujfs8/Y2Y3mFmjma01s43hPgoXD/L+5/Xf88HMvmZm14SPTzezh8IF9P4rXEYdM/sTM3s6jOsHRfvwIkOoihnNIoNx98fN7B7gr4AG4HvuPnB5ih8QrMz69fD57xKsA5UCLnH3DjObBjxqZvd4AbNBw/Wuvgpc7O57zOxy4K+BawkWrJvn7j19S42IjCclBal2NwOPE1TyfzLwpLs/YWYzzOw4YDqw3923hxX758PVUHME+2DMBNoLuOZ84FTg/nBpjBqCpZoBfk2wLMXdwN2v65OJHAUlBal2xwKNBBsJxQiWlx7oDuAyYBbBnQPAFQRJ4nR3T4crtw7cujHD4U20fecN2Ozug22N+dsEG+O8C7jJzE4JN1YRGRfqU5Bqtxq4CVhDsGb9YH5AsCrrZQQJAqCFYG+HtJm9DThhkNdtAxaaWb2ZtRCsugnwLDC9b79kM4ua2SlmFgHmuPsDBBsJTSFIWCLjRncKUrXM7Cog4+7/Fq4I+4iZvd3d/6d/OXffHC6N/aq79zXzrAF+ambrgScZZGlxd38l3Of31wTr3z8RHu81s8uAfwqTRS1Bv8VzwPfCYwZ8pQy3z5QJTqukiohInpqPREQkT0lBRETylBRERCRPSUFERPKUFEREJE9JQURE8pQUREQkT0lBRETy/j8ddITMfXatPwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib.pyplot import plot,xlabel,ylabel,legend\n",
    "from random import randint\n",
    "from math import log\n",
    "\n",
    "def simulate_y(x):\n",
    "    return (3 - 1/x**3)\n",
    "\n",
    "xs = list(range(1,10))\n",
    "ys = []        \n",
    "for x in xs:\n",
    "    ys.append(simulate_y(x))\n",
    "\n",
    "#Plot the result\n",
    "plot(xs,ys,'o-',alpha=0.50,label=\"exponential data\")\n",
    "xlabel(\"x values\")\n",
    "ylabel(\"y values (3-1/x**3)\")\n",
    "legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pearson r =  0.6191269478109577\n",
      "Pearson p= 0.07542274289051659\n",
      "Spearman r =  1.0\n",
      "Spearman p= 0.0\n"
     ]
    }
   ],
   "source": [
    "#Now let's fit Pearson and Spearman regressions\n",
    "from scipy.stats import spearmanr,pearsonr,linregress\n",
    "\n",
    "pearson_r,p = pearsonr(xs,ys)\n",
    "print(\"Pearson r = \",pearson_r)\n",
    "print(\"Pearson p=\",p)\n",
    "\n",
    "spearman_rho,p = spearmanr(xs,ys)\n",
    "print(\"Spearman r = \",spearman_rho)\n",
    "print(\"Spearman p=\",p)\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although we did not specify that the data were exponential specifically, we see a perfect (r = 1.0) Spearman regression, with a p value of 0.0. This is because the y variable is always monotonically increasing as x increases - which is exactly what Spearman regression is testing for. In contrast, Pearson regression assumes that we have a linear relationship between x and y. Data that are, in truth, exponential (and we know this because we simulated them) are fairly poorly explained by a linear regression. Indeed the linear regression is not statistically significant when applied to these data.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1a27cfcba8>"
      ]
     },
     "execution_count": 208,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VOX1+PHPmWSyb+xbgkHLviMBFAUVFbTU3aq1LrWK1dblW0vbrz+1ltp+v7Wt/aqttVSrbaVqQYuodasb4kKVRWVREAkkgBCyJzNJZjm/P+5kCCEbmMkkM+f9euVF5t5n7j0JcM/c59zneURVMcYYYwBc0Q7AGGNM92FJwRhjTJglBWOMMWGWFIwxxoRZUjDGGBNmScEYY0yYJQVjjDFhlhSMMcaEWVIwxhgTlhipA4tICrASSA6dZ5mq/qRZm+8DVwN+oAS4SlV3tHXcvn37an5+fkRiNsaYWLVmzZr9qtqvvXYRSwpAPXCKqtaIiBtYJSIvqOp7TdqsA6aqqkdErgPuBi5q66D5+fl88MEHkYvaGGNikIi0+YG7UcS6j9RRE3rpDn1pszavq6on9PI9IDdS8RhjjGlfRGsKIpIgIuuBfcArqrq6jebfBl6IZDzGGGPaFtGkoKoBVZ2EcwcwTUTGtdRORL4JTAV+1cr+BSLygYh8UFJSErmAjTEmzkWyphCmqhUi8gYwD9jQdJ+InAr8P2C2qta38v7FwGKAqVOnHjLXt8/no7i4mLq6us4O3XQDKSkp5Obm4na7ox2KMTEvkk8f9QN8oYSQCpwK/LJZm8nAH4F5qrrvSM9VXFxMZmYm+fn5iMiXitt0L6pKaWkpxcXFDBs2LNrhGBPzItl9NAh4XUQ+At7HqSk8JyKLROSsUJtfARnAUhFZLyIrjuREdXV19OnTxxJCDBIR+vTpY3eBxnSRiN0pqOpHwOQWtt/R5PtTO+t8lhBil/3dGtN1bESzMcZ0d+WFULiqS05lSSEO/eIXvzjo9fHHH9/uezIyMtptc+edd/LrX/+6zTbLly9n06ZN7R7LGBOy9Eq4dyI8813QQ56z6XRxmRQ276nkt69s4QdLP+S3r2xh857KaIfUpZonhXfeeafLzm1JwZg2lG6Dt34Dj3wV/A3OtmGz4bRFcPkz0AVdqXGXFDbvqWTxyu1Uen0Myk6h0utj8crtXzoxPPbYY0ybNo1JkyZx7bXXEggE2LFjB8OHD2f//v0Eg0FOPPFEXn75ZQoLCxk1ahRXXHEFEyZM4IILLsDjcQZ2v/rqq0yePJnx48dz1VVXUV/vPKWbn5/PT37yE6ZMmcL48eP55JNPAKitreWqq66ioKCAyZMn88wzzwDw6KOPct555zFv3jyGDx/OD3/4QwB+/OMf4/V6mTRpEpdeeilw4C6gpqaGOXPmhM/ReKy2/PznP2fkyJGceuqpfPrpp+Htf/rTnygoKGDixImcf/75eDwe3nnnHVasWMHChQuZNGkS27Zta7GdMXGlajes/BX84QS4fwq8uggC9VDzhbN/6rdg5k3QK79LwhHtgtuRzjR16lRtPvfR5s2bGT16NABvfLqPkuoWhzsA8PLGL/A2BEhNSghva3x9+tiBLb6nX2YyJ43s3+oxN2/ezA9/+EOefvpp3G43119/PTNmzODyyy/noYce4sUXX2T69Ol89tln/PGPf6SwsJBhw4axatUqZs6cyVVXXcWYMWP43ve+x/Dhw3n11VcZMWIEl19+OVOmTOHmm28mPz+fW265hRtuuIEHHniAtWvX8tBDD3HrrbcyZswYvvnNb1JRUcG0adNYt24dS5cuZdGiRaxbt47k5GRGjhzJqlWryMvLIyMjg5qamnD8ja/9fj8ej4esrCz279/PjBkz2Lp1KyJyyHsA1qxZw5VXXsnq1avx+/1MmTKF73znO/zgBz+gtLSUPn36AHDbbbcxYMAAbrjhBq688krmz5/PBRdcANBqu5Z+x41/x8b0eCVbIDHJudDveBcemQd502HM2TD6LMjJ6/RTisgaVZ3aXrsuGbzWnZR7fOSkHvxjp7hdlHt8R3zMV199lTVr1lBQUACA1+ulf38niVx99dUsXbqUBx98kPXr14ffk5eXx8yZMwH45je/yX333cdpp53GsGHDGDFiBABXXHEFv//977n55psBOO+88wA49thjefrppwF4+eWXWbFiRbgvv66ujp07dwIwZ84csrOzARgzZgw7duwgL6/1f2yqyq233srKlStxuVzs2rWLvXv3MnBgy8nyrbfe4txzzyUtLQ2As846K7xvw4YN3HbbbVRUVFBTU8PcuXNbPEZH2xnT4+37BDYth43LoWQzTLsWzrzbSQb/tQmyh0Q7QiAGk0Jbn+gBisu9VHp9ZKceGB1b6fUxNtXNhVOPLDurKldccQX/8z//c8g+j8dDcXEx4HTPZGZmAoc+ZikitHfXlpycDEBCQgJ+vz987qeeeoqRI0ce1Hb16tXh9s3f05olS5ZQUlLCmjVrcLvd5Ofntzs+oLXHRa+88kqWL1/OxIkTefTRR3njjTe+VDtjerQ/nwE73wEEjjoezrjbuSMAcLm6TUKAOKwpzBs3gEqvj0qvj6Bq+Pt54wYc8THnzJnDsmXL2LfPGZRdVlbGjh3OLLU/+tGPuPTSS1m0aBHXXHNN+D07d+7k3XffBeDxxx/nhBNOYNSoURQWFvLZZ58B8Le//Y3Zs2e3ee65c+dy//33hxPKunXr2o3X7Xbj8x16Z1RZWUn//v1xu928/vrr4Z+hNbNmzeKf//wnXq+X6upqnn322fC+6upqBg0ahM/nY8mSJeHtmZmZVFdXt9vOmB5JFb74GF67Cx47/8DTQqPOhDN/Dbd8At/6F0y/FrIGRTfWVsTcnUJ7Rg/KZsGsYby4YS+7KrwMyUnlooJcRg/KPuJjjhkzhrvuuovTTz+dYDCI2+3m97//PYWFhbz//vu8/fbbJCQk8NRTT/HII49w8sknM3r0aP7yl79w7bXXMnz4cK677jpSUlJ45JFHuPDCC/H7/RQUFPCd73ynzXPffvvt3HzzzUyYMAFVJT8/n+eee67N9yxYsIAJEyYwZcqUgy7El156KV/72teYOnUqkyZNYtSoUW0eZ8qUKVx00UVMmjSJo446ihNPPDG872c/+xnTp0/nqKOOYvz48eFEcPHFF3PNNddw3333sWzZslbbGdOjlG6D9UucrqGybSAuyD8B6iohNQeOP7RO1l3FXKG5JygsLGT+/Pls2LCh/cYG6Hl/xybGqcKeDyFzoPO1cTksuwqGnegUi0d9DTLaXeSsS1mh2RhjOpMq7F4Hm55xCsblhXDKbTBrIYyYBz/YAul9ox3ll2ZJIQry8/PtLsGYniTghwdmQOlWcCXC0SfBibfAqPnOfneK8xUDLCkYY0xTwSDsWuPcDVTvgQv+DAmJMP4CyM6FkWdCWu9oRxkxlhSMMQbgiw1OsXjTM1C1CxKS4Jg5EPBBghtO+nG0I+wSlhSMMfEpGISi1dB/tPOE0I534P2HnEQw5w4YeQakHPlTiT2VJQVjTPwIBmDne07X0KYVzvxCZ/0OplwGky6BiRdDSla0o4yquBu8FimNk8rt3r07PK9PvOvIlNzGdBlvBdwzBh49E9b+FfIK4PyHYew5zv7kzLhPCGB3Cp1u8ODBLFu2LKLn8Pv9JCa2/FfX1r6OCgQCJCQktN+wHV05JbcxBwn4Ycfbzh0BwPzfOl1EEy6EwZNh+FxIbn+NkHhkdwqdrLCwkHHjxgGtT18NzkR2xx13HFOmTOHCCy8Mz0C6aNEiCgoKGDduHAsWLAhPX3HSSSdx6623Mnv2bO69996DznnnnXeyYMECTj/9dC6//HICgQALFy6koKCACRMm8Mc//hGAYDDI9ddfz9ixY5k/fz5nnnlmOIHl5+ezaNEiTjjhBJYuXcq2bduYN28exx57LCeeeGJ4qu6lS5cybtw4Jk6cyKxZswDYuHFjeNrwCRMmsHXrVuDA3ZOqsnDhQsaNG8f48eN58sknAXjjjTc46aSTuOCCCxg1ahSXXnppu/M/GdOmov/AszfBb0bCX8+CD58Af/2B6SZOvwvGnW8JoQ2xeafwyFcP3Tb2HJh2DTR4YMmFh+6f9A2YfCnUlsI/Lj9437eeP+JQ1q9ff9D01TfccAOpqancdddd/Pvf/yY9PZ1f/vKX3HPPPdxxxx1873vf4447nGWsL7vsMp577jm+9rWvAVBRUcGbb77Z4nnWrFnDqlWrSE1NZfHixWRnZ/P+++9TX1/PzJkzOf3001mzZg2FhYV8/PHH7Nu3j9GjR3PVVVeFj5GSksKqVc6Sf3PmzOHBBx9k+PDhrF69muuvv57XXnuNRYsW8dJLLzFkyBAqKioAePDBB7npppu49NJLaWhoIBAIHBTb008/zfr16/nwww/Zv38/BQUF4YSybt06Nm7cyODBg5k5cyZvv/02J5xwwhH/vk2cCfhg+0pnSonEZNjyEny0FEbOc0YWf+U0SEqLdpQ9SmwmhW6kpemrKyoq2LRpU3jq7IaGBo477jgAXn/9de6++248Hg9lZWWMHTs2nBQuuuiiVs9z1llnkZqaCjh3IR999FH4LqCyspKtW7eyatUqLrzwQlwuFwMHDuTkk08+6BiNx6+pqeGdd97hwgsPJM/GxX5mzpzJlVdeyde//vXwVN7HHXccP//5zykuLua8885j+PDhBx131apVXHLJJSQkJDBgwABmz57N+++/T1ZWFtOmTSM3NxeASZMmUVhYaEnBtM3f4CSCTf+ET54Hbzl84x8wYq4zx9CsH4A7NdpR9lixmRTa+mSflNb2/vQ+X+rOoLmWpq9WVU477TQef/zxg9rW1dVx/fXX88EHH5CXl8edd9550NTV6enprYfdZJ+qcv/99x+yNsHzz7f9czUeIxgMkpOTc9D6D40efPBBVq9ezfPPP8+kSZNYv3493/jGN5g+fTrPP/88c+fO5aGHHuKUU045KJ7WHO703ibOlW2HxSdBXQUkZTqPjY49x1myEpy6gflSrKYQBTNmzODtt98OT5Ht8XjYsmVLOAH07duXmpqaIy5Yz507lz/84Q/h6bG3bNlCbW0tJ5xwAk899RTBYJC9e/e2unZBVlYWw4YNY+nSpYBzUf/www8B2LZtG9OnT2fRokX07duXoqIiPv/8c44++mhuvPFGzjrrLD766KODjjdr1iyefPJJAoEAJSUlrFy5kmnTph3Rz2biiL8ePn0R/nmds0QlQM5RMP5CuPhxWPgZnP8nGPXVmJliojuIzTuFbq5fv348+uijXHLJJeFumbvuuosRI0ZwzTXXMH78ePLz88MruR2uq6++msLCQqZMmYKq0q9fP5YvX87555/Pq6++yrhx4xgxYgTTp08Pd201t2TJEq677jruuusufD4fF198MRMnTmThwoVs3boVVWXOnDlMnDiR//3f/+Wxxx7D7XYzcODAcE2k0bnnnsu7777LxIkTERHuvvtuBg4cGC5eG3OQba85BeJPX4D6KmcA2ZQrnH0uF3z119GNL8bZ1NlxpqamhoyMDEpLS5k2bRpvv/12q8ttdif2dxzDfF7Y/hYMPw1EYMUNzsCyUfMPdA0lJkU7yh4v6lNni0gKsBJIDp1nmar+pFmbZOCvwLFAKXCRqhZGKiYD8+fPp6KigoaGBm6//fYekRBMDPJ5YesrzjxDW16EhhpY8IYzhuDUn8JX73HmGzJdLpLdR/XAKapaIyJuYJWIvKCq7zVp822gXFW/IiIXA78EWn/ExnxptgayibriNQQfnY/L76EmIZvP+55GztTzGTrAGd8TzRlIN++pPGhVxnnjBnypVRl7YlwRSwrq9EvVhF66Q1/N+6rOBu4Mfb8M+J2IiB5Bn5aqtrqIvOnZekIXp11MWtFQ64wd2LQccgvg+Bv4RHPZn3k62/qcQmm/qVTVQ+XnPhbkeqL6O9u8p5LFK7eTnepmUHYKlV4fi1duZ8GsYXEVV0QLzSKSAKwBvgL8XlVXN2syBCgCUFW/iFQCfYD9h3OelJQUSktL6dOnjyWGGKOqlJaWkpLSfZ8u6WkXk2tOGMbIQVmoKooz2FfR8KDfxtcHvuegtjRp33yfOjtJ2vIcKZ88TXLha4jfSyCtHzW9xlFV5uHJdfupGvRfZKYkonXOeQTlif8UccXxzvQq4XOFXx143bi5eZwHtzl4e1v7Gps8vnoH9f4AngbB0+A8Gl3vD/Dntwu5uCDvkOM0fe+h59L229DygZq3f2ptMbX1fucdquSkOfWVFzfs7XlJQVUDwCQRyQH+KSLjVLXpkmMtXcEP+VgoIguABQBDhw495A25ubkUFxdTUlLSOYGbbiUlJSU8wK0rP/mqKoGg4g99BQKKLxgkEFR8gWB435L3iggEg9T5Anh9AVTVuZisKuSCY3MJqnOsoEJQlaA6F9Smf7bVJhhu23Q/0Oy1ht7XuG3V1v14G/wkJyaE/1PV+QLc+dwmjju6T6f+rpL8NQyu/pDCXs6AzLM3PcKAms182Hc+W/rOYXfWJJQEWFPM6u1lZCQnHPQBTlWpqa8JX/Ci4eNdVWQkJ1Ba03BQXIX7PazMPqzPqYel6edYaXJJbNy+ZW81mckJVNcJgpCTlkRmSiK7KrwRiadLHklV1QoReQOYBzRNCsVAHlAsIolANlDWwvsXA4vBefqo+X63282wYcMiELnpTjbvqeSPK7eT7nbRK83NF1Ve7n11KxdNzWNY3wz8waBzAQ+ELuLBIL5Akwt7IBjafvBF3R9ovW1HrN1Z3uJFrrDew+CcUsD5D+4SwSUgIuHvXSIH7XO5BBHn0tC8TYLLhcvV+PrAvubHlNDrd7aVclSfNFwSGo4kgCqltT6OP8a5qxZxNjuhN30t4e0S2k6z1wkNVWTu+DeZnz9PWtGbuAL1FF72HwJZQ3CN/AP7UnLIdiUwrdn7qrw+qut9ZKUcKCTX1PvITHFzybShTc7Dgbg5cMEMxxKKs0mTJvsONGranjbe4w8EqarzkZV6IK4qr4/sVDfXn3zMwcdt8r6mx2p6/EO3t9ymPYGgUhmKo1F1nZ8hOZEZtR3Jp4/6Ab5QQkgFTsUpJDe1ArgCeBe4AHjtSOoJpmdzPlkH8TQEqK33U9vgp7Y+gCf0Z229H0+Dnxc3fkFNnZ8U94EZXOt8Af7w5udtfvJ1iZCYICS6hASX4E5wkeByXicmuEhLcpHgSsSdICS4XKHtEmrjavG9bpeLhNB2r89PjddPdnqSMxpUoNrrJzvNzY1zhocv2l1t/c6KQy4mlV4feb3Tmf5l7xS2vARPfhMCDZA5GKZeBWPPIT93uDOWgNYnnPt6QS6LV24nqJCZkkh1nR9fAC44NpeB2dHrJpw/cRCLV24nweUKx1VTH+Ab04eSnPjlZw0+UvPGDWDxyu3Agd9XpdfHRQW5ETlfJO8UBgF/CdUVXMA/VPU5EVkEfKCqK4CHgb+JyGc4dwgXRzAecxg6o5smGFQ8vtCFvt7f5kW/pU/liS4hLTmR9KQEskPdCsf0SycpMYEEl/OpGGB/TQMXT8tr+SLucuFyRfaCfO7kISxeuZ2khkD4P211vZ9LpjsxRUunXUy85fDJv5xi8dhznckjB02CgmuccQRDpoYSQceMHpTNglnDDvr3dVFBbtQL8xaXIyYGr5nO1bRA2fRismDWMEYNzKIhEMRTHwhf3Gsb/E1e+6ltCOCp94f61w89foo7gfTkBNKTEklPTiAt9Gd6ciLpSYmkJTnfJye6DvqE/dtXtrT4yTc71c1/nTaiK341rYr6Uz6dHZeqs17xxn/C529A0A/ZQ2H2QphyebtvN91PRwevWVIwh/jtK1vYV1WHP6g0+IM0BIJUeX0kJggF+b3xBQ79N5PgkvDFPC2p8YJ/6EU/zZ1AYsKRTbnVVrLqDhfgHq+2FPash6/McV4vPgk8pTDmHOeOYPCUgzvGTY8S9RHNpufxB4J8VlLD29v248LpB09wCUkJLlLdLmobAozPzSE9ybnQZyQnkhb6xJ/idkW837y73t73aDUl8Mmzzsji7W85o4gXbnMWobn0KWcgmSWCuGJJwbC/pp6Pd1XyyZ5q6nwB0twJpCclkN/X6b+HA900s0f0i2qsowdlWxLoLGv/Bs/eCBqE3kfDzJucO4Kk0DTs6Z372KrpGSwpxKkGf5Ate6vZsKuSPZV1JLiEY/plMH5INnPHDuBPb23H6wuSmOCK+NMOpgtU74XNK5w7guNvcBakyZsOJ3zfSQQDxtkdgQEsKcQVVWVvVT0bdlXy6d5qGvxB+mQkMWtEP8YMyiI16cBjd9ZNEwP89bDmL85TQzveART6jnAeIwXoNwLm3B7VEE33Y0khDtT5AnzyRTUf76pkf3U97gRh+IBMxg3JZnB2Sou1AOum6aEqd0HpZzBsFrjc8Pa9Tn1g9o+cO4L+Nv24aZslhRilquyq8LJhVyVb99bgDyr9s5I5ZVR/Rg7MPGgAmOnhKoud9Qc2LYei1ZDeH275BFwJcO1Kqw2Yw2JJIcZ4Gvxs2l3Fhl2VlHt8JCW6GDski3GDs+mf1X0nlTNH6K174NWfOt8PGAcn3+bcEbhCSd8SgjlMlhRiQDCo7CzzsGF3Jdv21RJUZUhOKgXDejO8fyZJibYUd0wo3+EUizcuhzPvhiHHQv6JcMrtzliCvl+JdoQmBlhS6MGq6nzhu4LqOj+pSQlMGprDuMFZ9MlIjnZ4pjPUV8MHf3YSwe61zrZBk5x1CgDyCpwvYzqJJYUeJhBUtu+vZcOuSgpLa1GFob3TmDWiH0f3TT/i0cKmGyn7HGr2wdAZ4EqEN38FfYc7y1SOORt624zAJnIsKfQQFZ4GNuyqYtOeSmrrA2QkJzItvzdjB2eTnWZr2fZ4pduceYY2PQNffAT9x8D174I7FW7+KKpLVJr4YkmhG2ucdmLDriqKyjyIwLC+6Ywbks2wPukRn/3TdJEXfgSrH3S+z50Gp//cuSNoZAnBdCFLClHU2gyW+2ucAWabQ9NOZKe6Of6YPowZnEVmit0V9Gglnzr1gU3PwDeehJw8GH465BwFY86CbBs1bqLLkkKUNF8/t7y2gV+9tIXxg7MgNBFd47QTeb1Tbe3pnqx2P/znT844gpJPAIGhx4G3zEkKX5lzYGZSY6LMkkKUvLhhb3gK6B2lteyvaaC23s+HxZXcMGc4owdlkpZkfz09kirs2+RMMzFkivP6rd84cw2d8SsY/TXIGhTtKI1pkV11omRXhZdB2Snsr6lnb1U9fTOSGDEgg5r6AMce1Sva4ZnDpQp7NxzoGirdCkefDJcvh4x+sHArpNrfq+n+LClEyZCcVCq9Pqq8PtwJwlf6Z1Dpjdxi3CbCll7pdA+JC/JPgBnXOXcEjSwhmB7CkkKUNK6fu6+qnr6ZSVR6bXrqHkHVWZ1s43Jn8fpvvwwpWTDuPDj6JCcRpPeNdpTGHDFLClEyelA2Fxfk8ZtXttDgV7JT3TY9dXdWUQT/Wex0DVXscAaVDZsNnv1OUmj6CKkxPZglhShKTUrguKP78K2Z+eSkJUU7HNOUKuxaA+40GDAGGmrgvQecu4FZC2HUV238gIlJ7SYFEekPzAQGA15gA/CBqgYjHFvMKyrzkpmSSHaqjT3oFoJB2PXBgWJxVTFMuBjO+yP0G+WsXZyaE+0ojYmoVpOCiJwM/BjoDawD9gEpwDnAMSKyDPiNqlZ1RaCxRlUpLvdwVJ90G4PQXTxyBhS9BwlJcMwpcMptMPIMZ5+IJQQTF9q6UzgTuEZVdzbfISKJwHzgNOCpCMUW00prG/A0BMjrbU8bdblgwFmMZuNy2PkOXPMGJCTClMth6lUwch6kWG3HxKdWk4KqLmxjnx9YHpGI4kRRmQeA3F5pUY4kjpRscYrFm1dAzV5ISIbhp0FdhfPE0ORLox2hMVHXVveRABcCCiwDTgHOBj4BHmyvpiAiecBfgYFAEFisqvc2a5MNPAYMDcXya1V95Ih/mh6kqNxLdqrb6gmRFPA7dwJZQ6DPMVC9G9b9zUkEY86BEXMhOTPaURrTrbTVffR7oD+QhJMMkoFncbqVRgI3tXNsP3CLqq4VkUxgjYi8oqqbmrT5LrBJVb8mIv2AT0Vkiao2HOHP0yMEg049YXh/uyB1uoAfdqxyuoY2P+s8Mnr8jXD6z5xVyhZucxayN8a0qK2kcKKqjhcRN/AFMEhVG0Tk7ziF5zap6h5gT+j7ahHZDAwBmiYFBTJDdyUZQBlOMolpJTX11PuCVk/obMEg/O5YKC90HiUdMde5Ixh+mrPflWAJwZh2tJUU/ACq6hOR9xs/vauqX0QCh3MSEckHJgOrm+36HbAC2A1kAhfFw6OuxeVWT/jSAj7Y/qZzR1C6Da56AVwuOO57kDEAvnIqJNnv15jD1VZS+EJEMlS1RlXnNW4UkYFAh7t3RCQD5wmlm1t4fHUusB6nXnEM8IqIvNW8nYgsABYADB06tKOn7raKyrz0Tk8iI9nGDh623evg/Yfgk+fBWw5Jmc7TQj6vs0rZtGuiHaExPVpbTx+d0cquapzHUdsV6np6Cliiqk+30ORbwP+qqgKfich2YBTwn2axLAYWA0ydOlU7cu7uKhBUdlV4GT3I6gkd4q+Hz9+AgROc6ab3fwabVjjjB8ac44wncKdEO0pjYkabH1VF5CZVvVdEblDV+wFUtRaobe/AoTrBw8BmVb2nlWY7gTnAWyIyAKeA/fnh/AA9zd6qOhr8QfKs66h1vjrY9pozqvjTF6C+0lmi8vjvOauTjTkLEpOjHaUxMam9/osaEVkIlB7BsWcClwEfi8j60LZbcR4/RVUfBH4GPCoiHwMC/EhV9x/BuXqM4nIvYPWEVjV44J7RztiBlGwYPd+5Izh6trPfkoExEdXWOIWfAOnAjcC9InKHqi7q6IFVdRXOhb6tNruB0zt6zFhQVOahb2YyqUkJ0Q4l+nxe+OzfTrHY54VL/u4Uh0/8PvQfC8NmQaJNFGhMV2qrpvBTEbkLZyqLuYeTEEzL/IEguyu8jM+N8ykUtr8FH/zZWY/AVwupvWHsOc7MpCIws70hMMaYSGkyBenFAAAarUlEQVSv++hNVX1LROyevRPsqazDH1TyesdZ11FDLWx92SkKp2Q7i9RsfxMmXOh0DeWfAAk2stuY7qC9pDAGeAUYDfw78uHEtuJyLyLEx5Kb9TWw9SWna2jrK+D3wvkPw/gLoOBqmH6dMwmdMaZbiWSh2TRTVO6hf2YKKe4YrSc0dv9U7YH7JjuJIL2/M9HcmHPgqOOddu44SIrG9FARKzSbg/kCQb6orGPy0Bibk7+uCra86NwRpPWCs3/vjCc48ftw1EwYOsOZXsIY0yNYobmL7K7wEghq7IxP+PQFWPtX+OxVCNRD5mCYdMmB/bN/GL3YjDFHrL3uo5VWaO4cxeVeXCIM7qn1BG+F87TQuPOdWsCOd2DPh1DwbadrKLfAmXvIGNOjtZkUVPVlEfm2qj7cuE1EEoDbVPWnEY8uhhSVeRiYnUxSYg+6cHrK4NN/OSOLt70OQZ/TNTRsFpz033DqTy0RGBNjOvL4xxwROR/4NtAHeAR4M6JRxZh6f4C9VfUU5PeKdijtaywWf/ExLD4Jgn7IGQozvgNjzoUhU5x2NgOpMTGp3aSgqt8QkYuAjwEPcImqvh3xyGLIrnIvQe3G4xNqS+GTZ507goHj4bRF0H8MnPB9Z+K5wZOdRGGMiXntJgURGY6zytpTOOMVLhORdarqiXRwsaK43EuiSxiY3c1m8/zwSfjw784IYw1Ar2FwzBxnnysBTvl/0Y3PGNPlOtJ99CzwXVV9NTTz6feB94GxEY0shhSVexiYnYI7Icr97zX7nNlHJ1zkfPLf9hpUFMEJNzvF4oHj7Y7AmDjX1jiFQaElNac1LnoTWvfgNyKyoqsC7OnqfAFKquuZcXSf6ARQvRc2r3C6hna8DRqEwVOg3wiY/1tnIJklAmNMSFt3Cn8WkV7AGyLyIrBKVRuX6NzaJdHFgOJyD6p0bT2hsVi87TX423mAQr9RMOuHMOZs6DvcaWfFYmNMM22uvCYiKcBJwLnAr0VkJ/Ai8KKq7uyaEHu2onIv7gRhYFaE6wmVuw7cEYw8E2be6IwdOOm/nUTQf1Rkz2+MiQntjVOoI5QEAERkGHAG8DsRGaiq0yIfYs9WXOZhcE4qCa4IddGsXgwblkHRauf1gHGQFuqqSs6Ek34UmfMaY2LSYU1TqarbgQeAB0TEVj9pR229n/01DYwalNV5B63YCUX/cWYbBWdwmc8Dp9zuFIv7fqXzzmWMiTtHNHexiHysquM7O5hY07j05pee76i80OkW2rgcdq8FcTlrE6T1hkset1lHjTGdpq2nj85rbRcwMDLhxJbicg9JiS76Zx7B1FGNxeIPn4B/XutsGzQJTr3TqRGk9Xa2WUIwxnSitu4UngSWANrCvm42Cqt7KirzkNsrFVdH6wml22DTcueOYMb1zqyj+Sc6cwyNORt6D4tswMaYuNdWUvgI+LWqbmi+Q0ROjVxIsaG6zke5x8f43HbWTwgGYNU9sPEZ2Puxsy23AFJCdYjsIc7gMmOM6QJtJYWbgapW9p0bgVhiSlFZqJ7Qu4XunZItUPIJjDnLmU5i43Jwp8HcX8DosyAnr4ujNcYYR1vjFN5qY98HkQkndhSXe0hxJ9AvI1RP2PfJga6hks2QlAkj5kFiElz9b6sNGGO6hcN6+khE1qrqlEgFE0uKyjzk5qQgIrDq/+DfPwEEhh4HZ9wNo7/mJASwhGCM6TYO95FUmySnLaqwdyN1Hz7NueuXUTb758AQ544gKR1GzXcWqTHGmG7qcJPC8x1tKCJ5wF9xHl8NAotV9d4W2p0E/B/gBvar6uzDjCn6Gmrhrd84YwlKPyNZXOzLmkK/rHRnf/9RNs2EMaZHONykcM9htPUDt6jqWhHJBNaIyCuquqmxgYjk4IyQnqeqO0Wk/2HGEx2qzvrE1V/AyHmQmALrH3cmmjvuu7xOAZ950rhm9NHRjtQYYw5LW4PXZgIP4XzKvwq4CzhGRNzA11X13bYOHJp2e0/o+2oR2QwMATY1afYN4OnGyfVUdd+X+FkiSxV2r3OKxZuecUYZ5wyFEXOdJ4huXAfuFFSVbW9tJ7dXqlNPMMaYHqStO4XfAl8HMnC6jc5R1VUiMgW4H5jZ0ZOISD4wGVjdbNcIwC0ibwCZwL2q+teOHjfiNDRuTwReuhXeewBciTBsNpx4C4z86oG1CNzOeL5yj4+aev+Xn9rCGGOioK2k4FbVjwFEpERVVwGEuoM6/LiMiGTgLOV5c+NiPc3OfywwB0gF3hWR91R1S7NjLAAWAAwdOrSjpz4ywSDsWnPgjuCSJ2DgOBh3PgwY60xL3TjFRAuKypxVSlscn2CMMd1cW0mh6dqR/91sX4dmSA11NT0FLFHVp1toUoxTXK4FakVkJTAROCgpqOpiYDHA1KlTW5p248ur3X+gWFy1C1xuZ9K5oN/ZnzvV+WpHcbmXzJREslPdEQnTGGMiqa2kcLuIpKmqR1WXN24UkWNwnipqU2g954eBzaraWoH6GZy1GRJxEs10nG6ryAsGnTUI/HVwzMmhYvESGHo8zLnDeYw0tZ0pKppRVYrKPeT3Sbd6gjGmR2prRHOL6zCr6jbg7g4ceyZwGfCxiKwPbbsVGBo6zoOqujm01OdHOAXth1qaa6nTBAOw871Q19AKqPkCcqc5SSE5A36wFRKPYEbTkP01DXgbAtZ1ZIzpsdp6+mgxcH9jXaHZvnTgIqBeVZe09P5QDaLdj8uq+ivgVx2O+Mt46mrY+LRzV/CVU2Hsuc7TQ42+REIAKCpvrCdYkdkY0zO11X30AE4X0nhgA1CCM2X2cCAL+DPO1No9x7FXwOj5MHyuc2fQyYrLveSkuclKsXqCMaZnaqv7aD3w9dDTQ1OBQYAXp0bwaRfF17mOPilihw4GleJyD8P7Z0bsHMYYE2ntjmhW1RrgjciH0rOV1NRT7wtaPcEY06O52m9iOqJxfEKuDVozxvRglhQ6SXG5lz4ZSWQkH+50UsYY030cVlIQEZeIZEUqmJ4qEFR2VXjJ7WVdR8aYnq3dpCAifxeRrNBjqJuAT0VkYeRD6zn2VtXR4A/afEfGmB6vI3cKY0JzFp0D/Atn8NllEY2qh7F6gjEmVnQkKbhDcxidAzyjqj4gMvMP9VDF5V76ZSaTmpQQ7VCMMeZL6UhS+CNQCKQDK0XkKKD5bKdxyx8IstvqCcaYGNGRcQr3Afc12bRDRE6OXEg9y57KOvxBtaktjDExoSOF5gEi8rCIvBB6PQa4IuKR9RBF5R5EYEiO3SkYY3q+jnQfPQq8BAwOvd4C3BypgHqa4nIvA7JSSHFbPcEY0/N1JCn0VdV/4Extjar6gUBEo+ohfIEgX1TWWT3BGBMzOpIUakWkD6EnjkRkBlAZ0ah6iN0VXgJBtfEJxpiY0ZE5Gb4PrACOEZG3gX7ABRGNqocoKvPiEmGw1ROMMTGiI08frRWR2cBInEVzPg2NVYh7xeUeBmWnkJRoU0gZY2JDu0lBRC5vtmmKiKCq7a7THMvq/QH2VtVTkN8r2qEYY0yn6Uj3UUGT71OAOcBaIK6Twq5yL0G18QnGmNjSke6jG5q+FpFs4G8Ri6iHKCr3kugSBmWnRDsUY4zpNEfSGe7BWac5rhWXexiUk0pigtUTjDGxoyM1hWc5MAGeCxgD/COSQXV3db4AJdX1zDi6T7RDMcaYTtWRmsKvm3zvB3aoanGE4ukRiss9qGL1BGNMzOlITeHNrgikJykq8+JOEAZmWT3BGBNbWk0KIlJNy+smCKCqGrfLchaVexjSK5UEl0Q7FGOM6VStVklVNVNVs1r4yuxIQhCRPBF5XUQ2i8hGEbmpjbYFIhIQkW4/Urq23k9pTYOtsmaMiUkdqSkAICL9ccYpAKCqO9t5ix+4JTQiOhNYIyKvqOqmZsdNAH6JMxNrt1dc7gWw+Y6MMTGpI+spnCUiW4HtwJs4q7C90N77VHWPqq4NfV8NbAaGtND0BuApYF/Hw46eojIPSYku+mcmRzsUY4zpdB15yP5nwAxgi6oOwxnR/PbhnERE8oHJwOpm24cA5wIPHs7xoqmo3ENur1RcVk8wxsSgjiQFn6qWAi4Rcanq68Ckjp5ARDJw7gRuVtXmazv/H/AjVW1zfQYRWSAiH4jIByUlJR09daerqvNR4fFZPcEYE7M6UlOoCF3YVwJLRGQfTr2gXSLixkkIS1T16RaaTAWeEBGAvsCZIuJX1eVNG6nqYmAxwNSpU1t6IqpLFJeF6gm9bapsY0xs6khSOBvwAv8FXApkA4vae5M4V/qHgc2qek9LbULdUY3tHwWea54QupOicg+pSQn0y7B6gjEmNnUkKSwAloZGMf/lMI49E7gM+FhE1oe23QoMBVDVHlNHAFBVisqcekLozsYYY2JOR5JCFvCSiJQBTwDLVHVve29S1VU4A906RFWv7GjbaKjy+qmu8zM13+oJxpjY1W6hWVV/qqpjge8Cg4E3ReTfEY+smykq9wCQ18vqCcaY2HU48z7vA74ASoH+kQmn+yoq85CenEDv9KRoh2KMMRHTkcFr14nIG8CrOE8IXaOqEyIdWHeiqhSVe8jrlWb1BGNMTOtITeEonDEG69ttGaPKPT5q6wM2PsEYE/M6MnX2j7sikO6sqCxUT7DxCcaYGGdrSXZAUbmHzJREslPd0Q7FGGMiypJCO1SV4nIveb2tnmCMiX0dKTR/T0R6dUUw3dH+mga8DQGbKtsYExc6cqcwEHhfRP4hIvMkzj4uN45PyLV6gjEmDnRk8NptwHCceYyuBLaKyC9E5JgIx9YtFJV5yElzk5Vi9QRjTOzrUE1BVRVn4NoXODOk9gKWicjdEYwt6oJBZVeF17qOjDFxo91HUkXkRuAKYD/wELBQVX0i4gK2Aj+MbIjRU1JTT70vSF5vSwrGmPjQkcFrfYHzVHVH042qGhSR+ZEJq3toHJ+Qa/MdGWPiREcGr93Rxr7NnRtO91JU7qFPRhLpyR3JncYY0/PZOIVWBILK7oo6qycYY+KKJYVW7K2qo8EftKktjDFxxZJCK4rKPIjAkBy7UzDGxA9LCq0oKvfSNyOZ1KSEaIdijDFdxpJCC/yBIHsqvPYoqjEm7lhSaMGeyjr8QbWlN40xcceSQguKyp16wuAcSwrGmPhiSaEFxWVeBmSlkOK2eoIxJr5YUmimwR/kiyobn2CMiU+WFJrZU+klEFQbn2CMiUuWFJopKvOS4BIGZVtSMMbEn4glBRHJE5HXRWSziGwUkZtaaHOpiHwU+npHRCZGKp6OKir3MDArhaREy5fGmPgTySufH7hFVUcDM4DvisiYZm22A7NVdQLwM2BxBONpV50vwN6qOltlzRgTtyI2/aeq7gH2hL6vFpHNwBBgU5M27zR5y3tAbqTi6YjdFV5UsSKzMSZudUkfiYjkA5OB1W00+zbwQlfE05qici+JLmFQdko0wzDGmKiJ+EIBIpIBPAXcrKpVrbQ5GScpnNDK/gXAAoChQ4dGKFJnErxBOakkJlg9wRgTnyJ69RMRN05CWKKqT7fSZgLOMp9nq2ppS21UdbGqTlXVqf369YtIrN6GACXV9Ta1hTEmrkXy6SMBHgY2q+o9rbQZCjwNXKaqWyIVS0fsqnCW3rRJ8Iwx8SyS3UczgcuAj0VkfWjbrcBQAFV9ELgD6AM84OQQ/Ko6NYIxtaqozEtSoosBWVZPMMbEr0g+fbQKkHbaXA1cHakYDkdRuYfBOSkkuNoM2RhjYppVVIHaej+lNQ32KKoxJu5ZUgCKy72A1ROMMcaSAs6jqMluF/0ykqMdijHGRJUlBZx6wpCcVFxWTzDGxLm4TwpVdT4qPD7rOjLGGCwpUFwWqidYkdkYYywpFJV7SE1KoG9GUrRDMcaYqIvrpKCqFJV5yO2VSmjwnDHGxLW4TgqVXh/VdX7rOjLGmJC4Tgo2PsEYYw4W10mhqMxDenICvdLc0Q7FGGO6hbhNCqpKUbmHvF5pVk8wxpiQuE0KZbUN1NYHrOvIGGOaiNukEK4nWJHZGGPC4jYpFJV7yExJJCs14iuSGmNMjxGXSUFVKS73ktfb6gnGGNNUXCaFkpp6vA0B6zoyxphm4jIpNNYTcnunRjkSY4zpXuIyKRSVechJc5OVYuMTjDGmqbhLCsFgqJ5gXUfGGHOIuEsK+6rrafAHbXyCMca0IO6SQlG5B4DcXlZPMMaY5uIuKRSXe+iTkUR6so1PMMaY5uIqKQSCyi6rJxhjTKviKil8UVWHL6Dk2aOoxhjTooglBRHJE5HXRWSziGwUkZtaaCMicp+IfCYiH4nIlEjFA86jqCIwJMfuFIwxpiWR7Fj3A7eo6loRyQTWiMgrqrqpSZszgOGhr+nAH0J/dqrNeyp5ccNe3t62n8zkRApLaxg9KLuzT2OMMT1exO4UVHWPqq4NfV8NbAaGNGt2NvBXdbwH5IjIoM6MY/OeShav3E65p4FEgUSXsHjldjbvqezM0xhjTEzokpqCiOQDk4HVzXYNAYqavC7m0MTxpby4YS/ZqW4SRFCEQTmpZKe6eXHD3s48jTHGxISIJwURyQCeAm5W1armu1t4i7ZwjAUi8oGIfFBSUnJY599V4SUzJRERyElzk5mSSGZKIrsqvId1HGOMiQcRTQoi4sZJCEtU9ekWmhQDeU1e5wK7mzdS1cWqOlVVp/br1++wYhiSk0p1nZ/MFDejBmaR6HJRXednSI49gWSMMc1F8ukjAR4GNqvqPa00WwFcHnoKaQZQqap7OjOOeeMGUOn1Uen1EVQNfz9v3IDOPI0xxsSESD59NBO4DPhYRNaHtt0KDAVQ1QeBfwFnAp8BHuBbnR3E6EHZLJg1jBc37GVXhZchOalcVJBrTx8ZY0wLIpYUVHUVLdcMmrZR4LuRiqHR6EHZlgSMMaYD4mpEszHGmLZZUjDGGBNmScEYY0yYJQVjjDFhlhSMMcaEifMAUM8hIiXAjiN8e19gfyeG01m6a1zQfWOzuA6PxXV4YjGuo1S13dG/PS4pfBki8oGqTo12HM1117ig+8ZmcR0ei+vwxHNc1n1kjDEmzJKCMcaYsHhLCoujHUArumtc0H1js7gOj8V1eOI2rriqKRhjjGlbvN0pGGOMaUNcJAUR+bOI7BORDdGOpSkRyROR10Vks4hsFJGboh0TgIikiMh/ROTDUFw/jXZMTYlIgoisE5Hnoh1LIxEpFJGPRWS9iHwQ7XgaiUiOiCwTkU9C/86O6wYxjQz9nhq/qkTk5mjHBSAi/xX6N79BRB4XkZRoxwQgIjeFYtoY6d9VXHQficgsoAZnPehx0Y6nUWg96kGqulZEMoE1wDmquinKcQmQrqo1oYWSVgE3hdbRjjoR+T4wFchS1fnRjgecpABMVdVu9Wy7iPwFeEtVHxKRJCBNVSuiHVcjEUkAdgHTVfVIxx91VixDcP6tj1FVr4j8A/iXqj4a5bjGAU8A04AG4EXgOlXdGonzxcWdgqquBMqiHUdzqrpHVdeGvq8GNtPJa1QfCXXUhF66Q1/d4tODiOQCXwUeinYs3Z2IZAGzcBa7QlUbulNCCJkDbIt2QmgiEUgVkUQgjRZWgoyC0cB7qupRVT/wJnBupE4WF0mhJxCRfGAysDq6kThCXTTrgX3AK6raLeIC/g/4IRCMdiDNKPCyiKwRkQXRDibkaKAEeCTU3faQiKRHO6hmLgYej3YQAKq6C/g1sBPYg7MS5MvRjQqADcAsEekjImk4C5PltfOeI2ZJoRsQkQyctaxvVtWqaMcDoKoBVZ2Es272tNAtbFSJyHxgn6quiXYsLZipqlOAM4Dvhrosoy0RmAL8QVUnA7XAj6Mb0gGh7qyzgKXRjgVARHoBZwPDgMFAuoh8M7pRgapuBn4JvILTdfQh4I/U+SwpRFmoz/4pYImqPh3teJoLdTe8AcyLcijgLPF6Vqj//gngFBF5LLohOVR1d+jPfcA/cfp/o60YKG5yl7cMJ0l0F2cAa1V1b7QDCTkV2K6qJarqA54Gjo9yTACo6sOqOkVVZ+F0hUekngCWFKIqVNB9GNisqvdEO55GItJPRHJC36fi/Gf5JLpRgar+t6rmqmo+TrfDa6oa9U9yIpIeelCAUPfM6Ti3/FGlql8ARSIyMrRpDhDVhxiauYRu0nUUshOYISJpof+bc3DqfFEnIv1Dfw4FziOCv7eIrdHcnYjI48BJQF8RKQZ+oqoPRzcqwPnkexnwcaj/HuBWVf1XFGMCGAT8JfRkiAv4h6p2m8c/u6EBwD+d6wiJwN9V9cXohhR2A7Ak1FXzOfCtKMcDQKhv/DTg2mjH0khVV4vIMmAtTvfMOrrPyOanRKQP4AO+q6rlkTpRXDySaowxpmOs+8gYY0yYJQVjjDFhlhSMMcaEWVIwxhgTZknBGGNMmCUFYzqZiNS038qY7smSgjHGmDBLCiZuiUiBiHwUWj8iPTRX/bhmbX4pItc3eX2niNwiIhki8qqIrA2to3B2C8c/qemaDyLyOxG5MvT9sSLyZmgCvZdC06gjIjeKyKZQXE9E7Ic3phVxMaLZmJao6vsisgK4C0gFHlPV5tNTPIEzM+sDoddfx5kHqg44V1WrRKQv8J6IrNAOjAYNzXd1P3C2qpaIyEXAz4GrcCasG6aq9Y1TjRjTlSwpmHi3CHgf5yJ/Y/OdqrpORPqLyGCgH1CuqjtDF/ZfhGZDDeKsgzEA+KID5xwJjANeCU2NkYAzVTPARzjTUiwHln+pn8yYI2BJwcS73kAGzkJCKTjTSze3DLgAGIhz5wBwKU6SOFZVfaGZW5sv3ejn4C7axv0CbFTVlpbG/CrOwjhnAbeLyNjQwirGdAmrKZh4txi4HViCM2d9S57AmZX1ApwEAZCNs7aDT0ROBo5q4X07gDEikiwi2TizbgJ8CvRrXC9ZRNwiMlZEXECeqr6Os5BQDk7CMqbL2J2CiVsicjngV9W/h2aEfUdETlHV15q2U9WNoamxd6lqYzfPEuBZEfkAWE8LU4uralFond+PcOa/Xxfa3iAiFwD3hZJFIk7dYgvwWGibAL/thstnmhhns6QaY4wJs+4jY4wxYZYUjDHGhFlSMMYYE2ZJwRhjTJglBWOMMWGWFIwxxoRZUjDGGBNmScEYY0zY/wcVFTAmuRB1CQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Plot the result\n",
    "plot(xs,ys,'o-',alpha=0.50,label=\"exponential data\")\n",
    "results = linregress(xs,ys)\n",
    "plot(xs,[results.slope*x+results.intercept for x in xs],'--',label=\"linear regression\")\n",
    "xlabel(\"x values\")\n",
    "ylabel(\"y values (3-1/x**3)\")\n",
    "legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The flexibility of rank-based methods has allowed us to detect an exponential trend, even though our model (Spearman regression) didn't explicitly make any assumptions about the data being exponential. In contrast, Pearson regression got poor results because it tried to explain an exponential distribution with a straight line.\n",
    "\n",
    "#### What exactly is the difference between Spearman and Pearson regression?\n",
    "\n",
    "Having seen this difference in performance, we might be wondering how exactly these types of regression differ. It turns out that the *only* difference of significance is that Spearman correlation coverts all the x and y values into ranks (smallest y, next smallest y, etc) and then just does Pearson regression on the ranks. \n",
    "\n",
    "We can prove to ourselves (and any other skeptical friends) that these two things are really the same. Here's how: we'll calculate a Spearman regression, then convert the data to ranks and calculate a Pearson regression. If we're right, these should give exactly the same result. We can then plot the rank data to show what Spearman regression is 'seeing' when it finds a perfect correlation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pearson r =  0.6191269478109577\n",
      "Pearson p= 0.07542274289051659\n",
      "Spearman r =  1.0\n",
      "Spearman p= 0.0\n",
      "Pearson r (calculated on ranks) =  1.0\n",
      "Pearson p (calculated on ranks)= 0.0\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import rankdata\n",
    "\n",
    "pearson_r,p = pearsonr(xs,ys)\n",
    "print(\"Pearson r = \",pearson_r)\n",
    "print(\"Pearson p=\",p)\n",
    "\n",
    "spearman_rho,p = spearmanr(xs,ys)\n",
    "print(\"Spearman r = \",spearman_rho)\n",
    "print(\"Spearman p=\",p)\n",
    "\n",
    "rank_ys = rankdata(ys)\n",
    "rank_xs = rankdata(xs)\n",
    "pearson_r,p = pearsonr(rank_xs,rank_ys)\n",
    "print(\"Pearson r (calculated on ranks) = \",pearson_r)\n",
    "print(\"Pearson p (calculated on ranks)=\",p)\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above code demonstrates that simply converting our x and y data into ranks (with rank 1 being the smallest number) is sufficient to turn Pearson regression into Spearman regression, and therefore to let it to detect any increasing trend in the dataa, not just a linear one. \n",
    "\n",
    "> Spearman regression is just normal Pearson linear regression, but with the input x and y data converted into ranks. \n",
    "    \n",
    "This was for me a very handy thing to realize - it mean that I didn't really have to learn two totally different approaches to regression, just an extremely simple trick that represents a different view into the data. The same is true of the statistical tests listed at the beginning of the chapter. In a single stroke, this halves the number of statistical tests we have to understand.\n",
    "\n",
    "> We don't have to memorize dozens of tests, we just have to realize that for every 'standard' parametric statistical test, there is generally some sort of rank-based alternative that can accomodate non-linear trends, just like Spearman regression is the rank based version of Pearson regression.\n",
    "\n",
    "#### Visualizing rank transformed data to understand what Spearman regression sees\n",
    "\n",
    "As a final way to make sure we understand how this rank transform is working, let's try replotting our data once it has been converted into ranks. Once we do, it will likely be immediately obvious why Spearman regression sees a perfect linear correlation in our exponential data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1a2a027e48>"
      ]
     },
     "execution_count": 210,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8lOW5+P/PlX2f7CEkQEBZAhgWA+IGiILWXaulVluPXdS239bWrz2tPefr1tNfz7E9p9rTc+qXY116tNWq1VO19QvIVhCURUAg7IsEYvZM9mSW6/fHPMSAASaQySSZ6/165ZWZJ8/z3NcMwzX33HM/1y2qijHGmKEvKtwBGGOM6R+W8I0xJkJYwjfGmAhhCd8YYyKEJXxjjIkQlvCNMSZCWMI3xpgIYQnfGGMihCV8Y4yJEDHhDqC77OxsLSoqCncYxhgzaGzcuLFGVXOC2XdAJfyioiI2bNgQ7jCMMWbQEJFDwe5rQzrGGBMhLOEbY0yEsIRvjDERYkCN4ffE4/FQXl5Oe3t7uEMx5owkJCRQWFhIbGxsuEMxEW7AJ/zy8nJSU1MpKipCRMIdjjG9oqrU1tZSXl7O6NGjwx2OiXAhTfgich/wDUCA/1LVJ3p7jvb2dkv2ZtASEbKysqiurg53KGYAKqtw8862So40tFGQnshVk/MozneFrL2QjeGLyGQCyX4mMAW4VkTGnuG5+jI0Y/qVvX5NT8oq3CxadQB3m4d8VwLuNg+LVh2grMIdsjZD+aVtMbBOVVtV1QusBG4KYXvGGDNovLOtEldiLFECn7jbcSXG4kqM5Z1tlSFrM5QJfxswW0SyRCQJuBoYceJOInK3iGwQkQ0D9WNvdHQ0U6dOZfLkydx66620traGO6RB5amnnuJ3v/tduMM4awcPHmTy5MnhDsMMEYfrWqlr7mD70UYqG9vx+ZXUhBiONLSFrM2QjeGrapmI/AuwBGgGtgDeHvZbBCwCKC0tPesV1UMxJpaYmMjmzZsBuP3223nqqae4//77z+qcPp+P6OjoszpHqNvwer3ExJz9S+Tee+8963P0tb56bMaciYM1LdQ0d+Bu81CUncyIjCSiowR3m4eC9MSQtRvSefiq+ltVna6qs4E6YE8o2+uPMbFLL72UvXv3AvDCCy8wc+ZMpk6dyj333IPP5wPgm9/8JqWlpUyaNImHH36469iioiIee+wxLrnkEl555RV+9atfMXHiREpKSvjiF78IQF1dHTfeeCMlJSXMmjWLrVu3AvDII4/w1a9+lblz5zJmzBh+9atf9RhfSkoKDz30EBdccAFr165l48aNzJkzh/PPP58rr7ySiooKANavX09JSQkXXnghP/jBD7p6rs899xy33nor1113HQsWLADg5z//OTNmzKCkpKTr8bS0tHDNNdcwZcoUJk+ezMsvvwzAj370o67H9MADD3TF/otf/AKAzZs3M2vWLEpKSrjpppuor68HYO7cufzwhz9k5syZjBs3jr/97W+feWyq2hXreeed19XmihUrmDt3LrfccgsTJkzg9ttvR/WzfYe5c+fy4x//mDlz5vDkk0/y5ptvcsEFFzBt2jSuuOIKKisrg36u9+/fz7Rp01i/fj3bt2/veh2UlJSwZ09IX+ZmEGvr9PHOtk94/cMjTC5IIz89gYykOETA3ebB3ebhqsl5IWs/1LN0clW1SkRGAjcDF57N+VbsqqK6qeOkf1+8/RPaOn00t3u6trV1+vi3xbtZMGlYj8fkpMYzd3xuUO17vV7++te/ctVVV1FWVsbLL7/MmjVriI2N5Vvf+hYvvvgiX/nKV/jpT39KZmYmPp+Pyy+/nK1bt1JSUgIE5mSvXr0agOHDh3PgwAHi4+NpaGgA4OGHH2batGm88cYbLFu2jK985Stdny527tzJ8uXLaWpqYvz48Xzzm9/8zNzulpYWJk+ezGOPPYbH42HOnDn8z//8Dzk5Obz88sv8wz/8A8888wx33XUXixYt4qKLLuJHP/rRcedYu3YtW7duJTMzk8WLF7Nnzx4++OADVJXrr7+eVatWUV1dzfDhw3n77bcBcLvd1NXV8frrr7Nz505EpOsxdfeVr3yFf//3f2fOnDk89NBDPProozzxxBNdz+8HH3zAX/7yFx599FGWLl163LF/+tOf2Lx5M1u2bKGmpoYZM2Ywe/ZsAD788EO2b9/O8OHDufjii1mzZg2XXHLJZ9pvaGhg5cqVANTX17Nu3TpEhKeffprHH3+cf/3Xfz3pc33Mrl27+OIXv8izzz7L1KlT+c53vsN9993H7bffTmdnZ9cbvzHHqCp7q5pZtrOKdo+fC0ZnMnN0Jnuqmo4bkVg4ozCks3RC/Zn2NRHJAjzAt1W1PpSN1bd6SE88/iElxEZR3+o5yRHBaWtrY+rUqUCgh/+1r32NRYsWsXHjRmbMmNG1T25u4I3jj3/8I4sWLcLr9VJRUcGOHTu6Ev7ChQu7zltSUsLtt9/OjTfeyI033gjA6tWree211wCYN28etbW1uN2BTyjXXHMN8fHxxMfHk5ubS2VlJYWFhcfFGh0dzec//3kgkJi2bdvG/PnzgcAQT35+Pg0NDTQ1NXHRRRcB8KUvfYm33nqr6xzz588nMzMTgMWLF7N48WKmTZsGQHNzM3v27OHSSy/lgQce4Ic//CHXXnstl156KV6vl4SEBL7+9a9zzTXXcO211x4Xm9vtpqGhgTlz5gBw5513cuutt3b9/eabbwbg/PPP5+DBg5/5d1i9ejW33XYb0dHR5OXlMWfOHNavX09aWhozZ87sei6mTp3KwYMHe0z43Z//8vJyFi5cSEVFBZ2dncfNk+/puQaorq7mhhtu4LXXXmPSpEkAXHjhhfz0pz+lvLycm2++mbFjz2gymhmimju8LN9Zxd6qZvLSErhpei65qQkAFOe7QprgTxTShK+ql/bl+U7XEy+vb8Pd5sGV+Gmv193mYVJiLLeWfub74qB1H8M/RlW58847+dnPfnbc9gMHDvCLX/yC9evXk5GRwd/93d8dd5VwcnJy1+23336bVatW8ec//5mf/OQnbN++vcehiGPT+uLj47u2RUdH4/V+5isREhISusbtVZVJkyaxdu3a4/Y5NoxyMt1jVFUefPBB7rnnns/st3HjRv7yl7/w4IMPsmDBAh566CE++OAD3n33XV566SV+/etfs2zZslO21d2xx3eyx9bTc3Pisac6/sTH9p3vfIf777+f66+/nhUrVvDII4+c9nwul4sRI0awZs2aroT/pS99iQsuuIC3336bK6+8kqeffpp58+YF8YjNUKaqbD/ayKo91fh8yqVjs5k+MoOoqPBN0x1StXSumpzXNQ7mVw3pmNjll1/Oq6++SlVVFRAYez906BCNjY0kJyfjcrmorKzkr3/9a4/H+/1+Dh8+zGWXXcbjjz9OQ0MDzc3NzJ49mxdffBEIjE1nZ2eTlpZ2RjGOHz+e6urqroTv8XjYvn07GRkZpKamsm7dOgBeeumlk57jyiuv5JlnnqG5uRmAI0eOUFVVxdGjR0lKSuKOO+7ggQceYNOmTTQ3N+N2u7n66qt54oknPvMm6XK5yMjI6Bqf/+///u+u3n4wZs+ezcsvv4zP56O6uppVq1Yxc+bMXj0n3bndbgoKCgB4/vnngzomLi6ON954g9/97nf8/ve/BwLj+WPGjOG73/0u119/fdf3LiZyuVs9/GnTEZbsqCQ7JZ47Zo2itCgzrMkeBkFphd4ozndx9+zR/TImNnHiRP7pn/6JBQsW4Pf7iY2N5T/+4z+YNWsW06ZNY9KkSYwZM4aLL764x+N9Ph933HEHbrcbVeX73/8+6enpPPLII9x1112UlJSQlJQUdCLqSVxcHK+++irf/e53cbvdeL1evve97zFp0iR++9vf8o1vfIPk5GTmzp2Ly9Xzc7RgwQLKysq48MLA1y8pKSm88MIL7N27lx/84AdERUURGxvLb37zG5qamrjhhhtob29HVfnlL3/5mfM9//zz3HvvvbS2tjJmzBieffbZoB/PTTfdxNq1a5kyZQoiwuOPP86wYcPYuXPnGT0/jzzyCLfeeisFBQXMmjWLAwcOBHVccnIyb731FvPnzyc5OZkdO3bwwgsvEBsby7Bhw3jooYfOKB4z+Pn9yubyBt7bW4OIcHlxLucVuAbMxXdyqo/J/a20tFRPXAClrKyM4uLiMEU0dDU3N5OSkgLAP//zP1NRUcGTTz4Z5qiGLnsdD321zR0sLavkaEM7o7OTmVecS1pC6AvmichGVS0NZt8h1cM3wXv77bf52c9+htfrZdSoUTz33HPhDsmYQcnnVzYcrOP9A3XExURx1eRhTBiWOmB69d1Zwo9QCxcuPG7GijGm9yob21m8o5Kapg7GD0tl7vgckuIGbloduJF1o6oD8t3SmGAMpGFT0zc8Pj/r9tey8VA9yXExXDdlOOfmpoQ7rNMa8Ak/ISGB2tpasrKyLOmbQedYPfyEhIRwh2L6yOG6Vt4tq6S+1cN5BS4uGZtNQmxoy6T0lQGf8AsLCykvL7d64mbQOrbilRncOrw+Vu+pYWu5G1diLLecX8iIzKRwh9UrAz7hx8bG2kpBxpiwOlDTwrtllTR3eJk+KoMLx2QRFzP4LmMa8AnfGGPCpa3Tx8rdVZRVNJGVEsfCkhHku0JXzTLULOEbY8wJVJXdlc2s2FVFh9fPrDFZzCjKICZ68PXqu7OEb4wx3TS1e1i2s4r91S0McyVwRXEeOanxpz9wELCEb4wxBHr1244Eip2pKrPH5TBtRHrY69/0JUv4xpiI19DaydKyKg7XtTIiM4krinNJT4oLd1h9zhK+MSZi+f3Kh4cbWLsvUOzsiuI8JhekDdlrfkK94tX3ga8DCnwE3KWq7ac+yhhjQq+muYMlOyr5xN3OmJxk5k3IJbUfip2FU8gSvogUAN8FJqpqm4j8Efgi8Fyo2jTGmJ6UVbi7yqbnuxIoSE+kqqmDuJgorj4vn3F5KUO2V99dqId0YoBEEfEAScDRELdnjDHHKatws2jVAVyJsaQmxLD+QB1L2jx8fnoBd8wqIjFucJRF6Ashm1SqqkeAXwAfAxWAW1UXh6o9Y4zpyTvbKklNiKGhtZMdRxuJjYmiOD+V5g5fRCV7CGHCF5EM4AZgNDAcSBaRO3rY724R2SAiG6xejjGmr+2ubOJAdTMV7nby0uIpKXRRmJHEkYa2cIfW70J52dgVwAFVrVZVD/An4KITd1LVRapaqqqlOTk5IQzHGBNJ2j0+lu6opLalkw6fn4n5aYzOTiEmKoqmdi8F6YO3RMKZCuUY/sfALBFJAtqAy4ENpz7EGGPO3r7qZpaVVdHS6eWGKfmsP1iPAn5Vmtq9uNs8LJwReRVMQ5bwVfV9EXkV2AR4gQ+BRaFqzxhjWju9rNhVza5PmshOjef6qcPJS0tg+qiMrlk6BemJLJxRSHG+K9zh9rsBv4i5Mcacjqqyq7KJFbuq6fT6mTk6kxlFmUQPobIIJ2OLmBtjIkZju4flTrGzfFcCV0zMIztlaBQ762uW8I0xg5Kq8tERN3/bU4OqMmd8DlMLh1axs75mCd8YM+jUt3SytKyS8vo2RmYmcUVxHq6koV0WoS9YwjfGDBp+v7Lp43rW7qslOlqYPzGPScOHbrGzvmYJ3xgzKFQ3BYqdVTa2c05uCvMm5JISbymsN0574ZWI3CciaRLwWxHZJCIL+iM4Y4zx+vy8t6+G37//MU3tHq4pyee6knxL9mcgmGfsq6r6pIhcCeQAdwHPAlYXxxgTUkcb2lhaVkltcyfF+WnMGZcTcfVv+lIwCf/Y4NjVwLOqukVswMwYE0Kd3kCvfvPhBlLiY7hxWgGjs5PDHdagF0zC3ygiiwkUQXtQRFIBf2jDMsZEqo9rW1lSVkljm4epI9K56Nws4mOsV98Xgkn4XwOmAvtVtVVEsggM6xhjTJ9p9/hYtbua7UcbyUiK5dbSQgozksId1pASTMJXYCJwLfAYkAwkhDIoY0xk2VvVzPKdVbR2+phRlMkFYzKJjQ5lMd/IFEzC/08CQzjzCCT8JuA1YEYI4zLGRICWjkCxs92VTeSkxnPD1OHkpll/MlSCSfgXqOp0EfkQQFXrRSQuxHEZY4YwVaWsoomVu6vx+PxcfG4254/KiIhiZ+EUTML3iEg0gaEdRCQH+9LWGHOGGts9vFtWycGaVoanJzB/4jAyk60P2R+CSfi/Al4HckXkp8AtwD+GNCpjzJCjqmwpd7Nmbw0Al03IZUqhy8oi9KPTJnxVfVFENhJYsUqAG1W1LOSRGWOGjLqWTpbuCCxAMioricuL83AlWrGz/nbahC8iI4FW4M3u21T149McNx54udumMcBDqvrEGcZqjBlkfE6xs3X7aomJjmLBpDwm5luxs3AJZkjnbQLj90JgOuZoYBcw6VQHqeouAvP3cb4DOEJgaMgYM0SVVbi7lhJMT4wlITaK6KgoxualcNn4XJKt/k1YBTOkc173+yIyHbinl+1cDuxT1UO9PM4YM0iUVbhZtOoAqQkx+Hx+Nhysw6vKdy47lysn54c7PEMQ1TJPpKqb6P0c/C8Cf+htW8aYweOdbZXERAVKIxx1t1OYmcS0EensqGgKd2jGEcwY/v3d7kYB04HqYBtw5uxfDzx4kr/fDdwNMHLkyGBPa4wZQDq9fjYeqqfT6yMhNpri/FRciXH4VTnS0Bbu8IwjmAG11G63vQTG9F/rRRufAzapamVPf1TVRcAigNLSUu3FeY0xA8DBmhbe3VmFx+cnLTGWCcPSui6gamr3UpCeGOYIzTHBjOE/epZt3IYN5xgz5LR7fKzcXc2Oo41kJsfxrcvO4Y0Pj9Lc4SU1IYamdi/uNg8LZxSGO1TjOGnCF5E3ca6u7YmqXn+6k4tIEjCf3n/Ja4wZwPZUNrF8VxVtnX4uGJ3JzNGZxERHkZsa3zVLpyA9kYUzCinOd4U7XOM4VQ//F2d7clVtBbLO9jzGmIGhucPL8p1V7K1qJjctnhun5ZGb+mmxs+J8lyX4AeykCV9VV/ZnIMaYgUtV2VHRyMrd1fh8yqVjs5k+MoMoK3Y2qAQzS2cs8DMCNfG73spVdUwI4zLGDBDutkCxs0O1rRSkJ3LFxDwrdjZIBTNL51ngYeCXwGUEVruyt3VjhjhVZfPhBt7bVwvAvAm5lFixs0EtmISfqKrviog4V8o+IiJ/I/AmYIwZgmqbO1haVsnRhnaKsgPFztISrNjZYBdMwm8XkShgj4j8LwI1cXJDG5YxJhx8fmXjoXrW7a8lNjqKKycNozg/1Xr1Q0QwCf97QBLwXeAnBIZ17gxlUMaY/lfV2M7iHZVUN3UwLi+VueNzrNjZEBPMv6ZXVZuBZgLj98aYIcTj8/P+/jo2HqonKS6a66YM59zclHCHZUIgmIT/byKSD7wCvKSq20MckzGmn5TXt7J0RyX1rR4mF7i4dGw2CbHR4Q7LhEgwpRUuE5FhwBeARSKSBrysqv8U8uiMMSHR4fWxZm8NWw67SUuM5fPTCxmZlRTusEyIBTVAp6qfAL8SkeXA3wMPAZbwjRmEDtS08G5ZJc0dXqaNTOeic7KJi+l1pXQzCAVz4VUxsJDA4uW1wEvA/w5xXMaYPtbWGSh2VlbRSFZKHF84bwTDrZJlRAn2wqs/AAtU9WiI4zHG9DFVZU9VM8t3VtHu8XPBmExmFgWKnZnIEswY/qz+CMQY0/eaO7ws21nFvqpm8tISuHl6Hjmp8eEOy4SJTbI1ZghSVbYfbWTVHit2Zj5lCd+YIcbd6mFpWSUf17VSkJHI/OI8MqzYmSG4L20nq+q2/gjGGHPm/H5lc3kD7+2tQUS4vDiX8wqs2Jn5VDA9/KechcifA36vqg3BnlxE0oGngckEVs/6qqquPZNAjTGfKqtwH7ey1IVjMjlY20qFu53R2cnMK861YmfmM4L50vYSpyb+V4ENIvIB8KyqLgni/E8C76jqLc6bhl3ZYcxZKqtws2jVAVyJsQxLi2fXJ40s3v4JF56Txa2lI5gwzIqdmZ4Fe+HVHhH5R2AD8CtgmgReUT9W1T/1dIxzRe5s4O+cc3QCnX0RtDGR7J1tlbgSY4mJErYfbaS100eeK56E2GiK89PCHZ4ZwIIZwy8hUDTtGmAJcJ2qbhKR4cBaoMeED4wBqoFnRWQKsBG4T1VbTjj/3cDdACNHjjzTx2FMxDhc34rf76fC3UFcdBTjh6XiSoylwt0e7tDMABfMlRe/BjYBU1T126q6CcC5COsfT3FcDDAd+I2qTgNagB+duJOqLlLVUlUtzcnJ6fUDMCaSHK5rpba5gwM1reSmxlMywkVGUhxN7V4K7KpZcxrBjOHPPsXf/vsUh5YD5ar6vnP/VXpI+MaY0+vw+li9p4at5W7G56VyqLaVrJR4okRwt3lwt3lYOKMw3GGaAS5k8/BV9RMROSwi41V1F3A5sCNU7RkzVO2vbmbZziqaO7xMH5XBRedksbeq6bhZOgtnFFKc7wp3qGaAC/WFV98BXnRm6OzHFlAxJmiBYmdVlFU0kZ0SxzUlI8h3BYZtivNdluBNr4U04avqZqA0lG0YM9SoKrsrm1m+q4pOr59ZY7KYOTqTaCuLYM7SSRO+iLxJ4GKpHqnq9SGJyJgI1tTuYdnOKvZXtzDMlcAVxVbszPSdU/Xwf+H8vhkYBrzg3L8NOBjCmIyJOKrKtiOBYmeqyuxxOUwbkW7FzkyfOmnCV9WVACLykxNm6rwpIqtCHpkxEaKhtZOlZVUcrmtlRGYSVxTnkp5kxc5M3wtmDD9HRMao6n4AERkN2IR5Y86S3698eLiBtfsCxc6uKM5jckGalUUwIRNMwv8+sEJE9jv3i4B7QhaRMRGgprmDJTsq+cTdzpicZOZNyCXVip2ZEAvmwqt3nOJpE5xNO1W1I7RhGTM0+fzKBwfqWH+wjriYKK4+L59xeSnWqzf9IphaOknA/cAoVf2GiIx1LqZ6K/ThGTN0VLjbWLqjkprmTorzU5kzLpfEuOhwh2UiSLCLmG8ELnTulwOvAJbwjQmCx+fnvX21fPhxPSnxMdwwdThjclLCHZaJQMEk/HNUdaGI3Aagqm1inz+NCcrhulaW7KjE3eahpNDFJWOziY+xXr0Jj2ASfqeIJOJchCUi5wA2hm/MKbR7AsXOPjriJj0pllvOL2REpq3/Y8IrmIT/MPAOMEJEXgQuxlnUxBjzWfuqm1lWVkVLp5fSogxmjckiNjqYSuTGhFYws3SWiMgmYBYgBBYxqQl5ZMYMMq2dXlbsqmbXJ01kp8Zz/dTh5KUlhDssY7oEWzytAIh29p8tIpxsaUNjIo2qsvOTJlburqbT6+fCc7KYUWTFzszAE8y0zGeAEmA74Hc2Kydf2tCYiNHY7mFZWRUHalrIdyUwf2IeWSlW7MwMTMH08Gep6sSQR2LMIKKqbC13s3pvDarKnPE5TC20YmdmYAsm4a8VkYmq2uvVqkTkINAE+ACvqlptfDPo1bd0sqSskiP1bYzMTOKK4jxcSVYWwQx8wST85wkk/U8ITMcUQFW1JMg2LrMvec1gVVbh7lpKcLgrgZGZSVS424mOFuZPzGPScCt2ZgaPYBL+M8CXgY/4dAzfmCGvrMLNolUHcCXGkpYYw/qD9SzeUckNU4fz5QuLSIkP9QqhxvStYF6xH6vqn8/w/AosFhEF/q+qLjrD8xjT797ZVklqQgyNbR6ONrQREy1MGJZKu8dnyd4MSsG8aneKyO+BN+l2hW2Q0zIvVtWjIpILLBGRnap63OIpInI3cDfAyJEjg4/cmBDbU9VES7uXdq+f7JR4irKSiIoSjjS0hzs0Y85IMAk/kUCiX9BtW1DTMlX1qPO7SkReB2YCq07YZxGwCKC0tPSka+ga0186vX7e21dDbXMnflUm5qd1rUDlbvNQkJ4Y5giNOTOnTPgiEg1sVdVf9vbEIpIMRKlqk3N7AfDYmYVpTP84VNvC0rIqGts8XHPeMDZ+3ICI4Felqd2Lu83DwhmF4Q7TmDNyyoSvqj4RuR7odcIH8oDXnRkMMcDvVfWdMziPMSHX7vGxanc12482kpEUy62lhRRmJDFj9KezdArSE1k4o5DifFe4wzXmjAQzpPOeiPwaeBloObZRVTed6iBnDdwpZxeeMaG3t6qJZTuraOv0M6MokwvGZHYVOyvOd1mCN0NGMAn/Iud39+EYBeb1fTjG9J+WjkCxs92VTeSkxnPj1DxyrdiZGcKCqZZ5WX8EYkx/UVXKKgLFzjw+Pxefm835ozKs2JkZ8oIpnuYiUBN/trNpJfCYqrpDGZgxoeBu87BsZyUHa1oZnp7A/InDyEyOC3dYxvSLYK+03QZ8wbn/ZQLr3N4cqqCM6WuqypZyN2v2Bqp8zB2fw9QR6VYWwUSUYNe0/Xy3+4+KyOZQBWRMX6tr6WTpjsBMm1FZSVxenIcr0YqdmcgTTMJvE5FLVHU1gIhcDLSFNixjzp7Pr2z6uJ51+2qJjhYWTMpjYr4VOzORK5iEfy/wO2csH6AeuDN0IRlz9qoa21lSVklVYwdj81K4bHwuyVb/xkS4k/4PEJH7VPVJIEVVp4hIGoCqNvZbdMb0ktfn5/0DdWw4WE9iXBTXluQzNi813GEZMyCcqstzF/Ak8O/AdEv0ZqA70tDG0h2V1LV0MnF4GnPG5ZAQGx3usIwZME6V8MucFatyRGRrt+29XQDFmJDq8Pp4b28tW8obSE2I5ebpBYzKSg53WMYMOCdN+Kp6m4gMA/4fcH3/hWRM8A7WtLC0rJLmDi9TRqRz8TnZxMVEhTssYwak0xVP+4Ru9XBEJAMYoapbT36UMaHX7vGxcnc1O442kpkcx62lI6xssTGnEcyVtisI9PBjgM1AtYisVNX7QxybMT3aU9nE8l2BYmczR2dywehMYqKtV2/M6QQzT82lqo0i8nXgWVV9+IQxfWP6RXOHl+U7q9hb1UxuWjw3TssjN9WKnRkTrGASfoyI5BMorfAPIY7HmM9QVbYfbWTVnmp8PuWSsdmcPzKDKCt2ZkyvBJPwHyPwxe1qVV0vImOAPaENy5gAd5uHd8sqOVTbSkF6IldMzLNiZ8bvNOl1AAAXcElEQVScoWDKI78CvNLt/n7g8yc/4njOMokbgCOqeu2ZBGkij9+vbClv4L19tQDMm5BLSaHLyiIYcxaC+dI2B/gGUNR9f1X9apBt3AeUAWlnEJ+JEGUVny4lmJkcS1JsNIpQlJ3EvAlW7MyYvhDMkM7/AH8DlgK+3pxcRAqBa4CfAjarx/SorMLNolUHSEuIQf1+PjhQh8enfGvOOVxdkm+9emP6SDAJP0lVf3iG538C+HvAipmYk3pnWyVx0cLHda20dvoY7kokMzmO3VXNXGPJ3pg+E8zk5bdE5OrenlhErgWqVHXjafa7W0Q2iMiG6urq3jZjBjmPz8+HH9dzsKYFj8/PuLxUxualkpEcx5EGq8JtTF8KJuHfRyDpt4lIo4g0iUgwhdQuBq536vG8BMwTkRdO3ElVF6lqqaqW5uTk9Cp4M7iV17fy4rpDtHv9pCTEMGVEetcMnKZ2r105a0wfO23CV9VUVY1S1URVTXPun/YLWFV9UFULVbUI+CKwTFXv6IOYzSDX4fWxbGclr2wox6dw7+wxpCXG0dLhw6+Ku82Du83DVZPzwh2qMUNKUCtCODV0xgJdlzWq6qpQBWWGrgM1LbzrFDubNjKdi5xiZ/npCV2zdArSE1k4o5DifNfpT2iMCVow0zK/TmBYp5BALZ1ZwFpgXrCNqOoKYMUZRWiGhLbOQLGzsopGslLi+MJ5IxjebcimON9lCd6YEAumh38fMANYp6qXicgE4NHQhmWGClVlT1Uzy3dW0e7xc8GYTGYWWbEzY8IhmITfrqrtIoKIxKvqThEZH/LIzKDX3OFl2c4q9lU1k5eWwM3T88hJjQ93WMZErGASfrmIpANvAEtEpB44GtqwzGB2YrGzS8dmM92KnRkTdsHU0rnJufmIiCwHXMA7IY3KDFruVg9Lyio5XNdKQUYi84vzyLBiZ8YMCKdM+CISBWxV1ckAqrqyX6Iyg47fr2wub+C9vTWICJcX53JegRU7M2YgOd0Sh34R2SIiI1X14/4KygwuNc0dLN1RSYW7ndHZycwrziUtwYqdGTPQBDOGnw9sF5EPgJZjG1XVFjaPcD6/sv5gHR8cqCMuJorPnTeM8Xmp1qs3ZoAKJuHbFEzzGZ+421lSVklNUwfjh6Uyd3wOSXFBXcdnjAmTYP6HXn1itUwR+RfAxvMjkMfnZ+2+WjZ9XE9yXAzXTx3OOTkp4Q7LGBOEYK5+md/Dts/1dSBm4Dtc18oL6w6x8VA9k4e7+PKFoyzZGzOInLSHLyLfBL4FjBGRrd3+lAqsCXVgZuBo9/hYs7eGreVuXImx3HJ+ISMyk8IdljGml041pPN74K/Az4AfddvepKp1IY3KDBj7q5tZtrOK5g4v00dlcNE5WcRaWQRjBqWTJnxVdQNu4LZj20TkblVd1B+BmfBq7fSyclc1Oz9pIjsljmtKRpDvsvr0xgxmvZ1WcS9gCX8IU1V2VTaxYlc1nV4/s8ZkMXN0JtFWFsGYQa+3Cd/+1w9hTe0elu2sYn91C8NcCcyfmEd2ihU7M2aoCKYefryqdjh3r3O2Zdo4/tChqmw7Eih2pqrMHpfDtBHpVuzMmCEmmB7+n0TkRlX1qGq5iOQDbwHnn+ogEUkAVgHxTjuvqurDZx2xOStlFe7jVpa6+JwsDtW1Ul7fxojMJK4oziU9yYqdGTMUBZPw3wBeEZHPAyOAPwMPBHFcBzBPVZtFJBZYLSJ/VdV1Zx6uORtlFW4WrTqAKzGWYWnx7KlqYsmOT5hRlMkt549gckGalUUwZggLpjzyf4lIHIHEXwTco6rvBXGcAs3O3VjnR888VHO23tlWiSsxlthooayiieYOL9kp8STFRXNeoS0vaMxQd6oLr+7vfpdA734zMEtEZqnqv53u5CISDWwEzgX+Q1Xf72Gfu4G7AUaOHNm76E2vlNe34leloqGdmGhhbG4KGcmxVLg7Tn+wMWbQO1UPP/WE+6+fZPtJqaoPmOqsmPW6iExW1W0n7LMIZ6pnaWmpfQIIkQp3G3UtndS1dFKYkUhRVjIx0VG42zwUpNv8emMiwakuvOqzKpmq2iAiK4CrgG2n2d30oU6vn7X7a/nw43rOzUnmYFQUOakJREUJ7jYP7jYPC2cUhjtMY0w/CGZa5jgCX9IWdd9fVeed5rgcwOMk+0TgCuBfzipa0yuH61pZsqMSd5uHkkIXl4zNZn9183GzdBbOKKQ438bvjYkEwczSeQV4Cnga8PXi3PnA8844fhTwR1V9q/chmt5q9/j4254ath1xk550fLGz4nyXJXhjIlQwCd+rqr/p7YlVdSswrfchmbOxr7qZZWVVtHR6KS3KYNYYK3ZmjAkIJuG/KSLfIvClbdd0DrvSdmBp7fSyYlc1uz5pIjs1nuunDicvLSHcYRljBpBgEv6dzu8fdNumwJi+D8f0lqqy85NAsTOPz89F52RRWmTFzowxnxXMhVej+yMQ03uN7R6WlVVxoKaFfKfYWZYVOzPGnERQ1TJFZDIwEegaI1DV34UqKHNqqsrWcjer99agqswZn8PUQit2Zow5tWCmZT4MzCWQ8P9CYD3b1YAl/DCob+lkSVklR+rbGJmZxBXFebiSYsMdljFmEAimh38LMAX4UFXvEpE8AlM0TT/y+5VNH9ezdl8t0dHC/Il5TBpuxc6MMcELJuG3qapfRLwikgZUYV/Y9quqpnaW7KikqrGDc3JTmDchl5T43q5dY4yJdMFkjQ1OLZz/IlAIrRn4IKRRGQC8Pj8fHKhj/cF6EmKjuLYkn3NzU6xXb4w5I6dM+BLILD9T1QbgKRF5B0hzLqoyIXS0oY2lZZXUNndSnJ/GnHE5JMZFhzssY8wgdsqEr6oqIm/grG6lqgf7I6hI1un1896+GjYfbiAlPoabphVQlJ0c7rCMMUNAMEM660RkhqquD3k0Ee5QbQtLy6pobPMwdUQ6F52bRXyM9eqNMX0jmIR/GXCPiBwCWggshqKqWhLSyCJIu8fHqt3VbD/aSEZSLLeWFlKYkRTusIwxQ0wwCf9zIY8igu2tamLZziraOv3MHJ3JBaMzibFiZ8aYEAimtMKh/ggk0rR0eFm+q4o9lc3kpMZz49Q8cq3YmTEmhGwydz9TVcoqmli5uxqvz8/F52Zz/qgMK3ZmjAk5S/j9yN3mYdnOSg7WtDI8PYH5E4eRmRwX7rCMMREiZAlfREYQqLczDPADi1T1yVC1N9CUVbi7LSWYwOjsZI40tANw2YRcphS67AIqY0y/CuW3g17gf6tqMTAL+LaITAxhewNGWYWbRasO4G7zkJ4Uy4aD9Ty1cj+qfu6YNYqpI9It2Rtj+l3IeviqWgFUOLebRKQMKAB2hKrNgeKdbZWkJcTQ3OHlSH0rUSKMy0vB41NciVbZ0hgTHv0yhi8iRQTWt32/h7/dDdwNMHLkyP4IJ+T2VTfT0uGhtdNPZnIco7OTiY6SriEdY4wJh5BP+BaRFOA14Huq2nji31V1kaqWqmppTk5OqMMJKa/Pz5q9NdQ0ddDU7mVcXgrj8lKJjY6iqd1LQXpiuEM0xkSwkPbwRSSWQLJ/UVX/FMq2wu1IQxtLd1RS19LJgol5bDniJjoqCr8qTe1e3G0eFs4oDHeYxpgIFspZOgL8FihT1X8LVTvh1uH18d7eWraUN5CaEMvN0wsYlZV8wiydRBbOKKQ43xXucI0xESyUPfyLgS8DH4nIZmfbj1X1LyFss18drGlhaVklzR1epoxI5+JzsomLCYySFee7LMEbYwaUUM7SWU2g0NqQ0+7xsWJXNWUVjWQmx/GF0hEMt/F5Y8wAZ1fa9tKeykCxs3aPnwtGZzLTip0ZYwYJS/hBau7wsnxnFXurmslNi+em6XnkplqxM2PM4GEJ/zRUle1HG1m1pxqfT7lkbDbnj8wgyoqdGWMGGUv4p+Bu9bC0rJKP61opyEhkfnEeGVbszBgzSFnC74Hfr2wpb2DN3hpEhHkTcimxYmfGmEHOEv4Japs7WFpWydGGdkZnJzOvOJe0BKt/Y4wZ/CzhO3x+ZcPBOt4/UEdsdBRXTR7GhGGp1qs3xgwZlvCBysZ2Fu+opKapg3F5qVw2IYekOHtqjDFDS0RnNY/Pz7r9tWw61EBSXDTXTRnOubkp4Q7LGGNCImITfnl9K0t3VFLf6mFygYtLx2aTEBsd7rCMMSZkIi7hd3h9rNlbw5bDblyJsXx+eiEjs5LCHZYxxoRcRCX8AzUtvOsUO5s+KoMLx2R1FTszxpihLiISflunj5W7qyiraCIrJY6FJSPId1mxM2NMZBnSCV9V2V3ZzIpdTrGzMZnMLLJiZ8aYyDRkE35zh5d3yyrZX91CXloCN0/PIyc1PtxhGWNM2IRyxatngGuBKlWdHKp2TnRisbPZ47KZNsKKnRljTCh7+M8BvwZ+F8I2jltKMDs5jqT4aHx+KMxIZP7EPNKTrNiZMcZAaFe8WiUiRaE6PwSS/aJVB0hLiAGU9w/U0uFT7p09huumDLeyCMYY082gHsN/Z1slKfHRlNe30dzhJS8tgazkOPZVt1iyN8aYE4Q94YvI3cDdACNHjuzVsUca2sh3xVPX0skwVwrZKXH4NbDdGGPM8cI+P1FVF6lqqaqW5uTk9OrYgvREmtp9nJubSnZKPCA0tXspsAXFjTHmM8Ke8M/GVZPzcLd5cLd58Kt23b5qcl64QzPGmAEnZAlfRP4ArAXGi0i5iHytr9sozndx9+zRuBJjqXC340qM5e7ZoynOd/V1U8YYM+iFcpbObaE6d3fF+S5L8MYYE4RBPaRjjDEmeJbwjTEmQljCN8aYCGEJ3xhjIoQlfGOMiRCiquGOoYuIVAOHzvDwbKCmD8PpKxZX71hcvWNx9c5QjGuUqgZ11eqASvhnQ0Q2qGppuOM4kcXVOxZX71hcvRPpcdmQjjHGRAhL+MYYEyGGUsJfFO4ATsLi6h2Lq3csrt6J6LiGzBi+McaYUxtKPXxjjDGnMOgTvog8IyJVIrIt3LEcIyIjRGS5iJSJyHYRuS/cMQGISIKIfCAiW5y4Hg13TN2JSLSIfCgib4U7lu5E5KCIfCQim0VkQ7jjOUZE0kXkVRHZ6bzWLhwAMY13nqdjP40i8r1wxwUgIt93XvfbROQPIpIQ7pgAROQ+J6btoX6uBv2QjojMBpqB36nq5HDHAyAi+UC+qm4SkVRgI3Cjqu4Ic1wCJKtqs4jEAquB+1R1XTjjOkZE7gdKgTRVvTbc8RwjIgeBUlUdUPO3ReR54G+q+rSIxAFJqtoQ7riOEZFo4Ahwgaqe6fU1fRVLAYHX+0RVbRORPwJ/UdXnwhzXZOAlYCbQCbwDfFNV94SivUHfw1fVVUBduOPoTlUrVHWTc7sJKAMKwhsVaECzczfW+RkQ7/giUghcAzwd7lgGAxFJA2YDvwVQ1c6BlOwdlwP7wp3su4kBEkUkBkgCjoY5HoBiYJ2qtqqqF1gJ3BSqxgZ9wh/oRKQImAa8H95IApxhk81AFbBEVQdEXMATwN8D/nAH0gMFFovIRmcN5oFgDFANPOsMgz0tIsnhDuoEXwT+EO4gAFT1CPAL4GOgAnCr6uLwRgXANmC2iGSJSBJwNTAiVI1Zwg8hEUkBXgO+p6qN4Y4HQFV9qjoVKARmOh8pw0pErgWqVHVjuGM5iYtVdTrwOeDbzjBiuMUA04HfqOo0oAX4UXhD+pQzxHQ98Eq4YwEQkQzgBmA0MBxIFpE7whsVqGoZ8C/AEgLDOVsAb6jas4QfIs4Y+WvAi6r6p3DHcyLn4/8K4KowhwJwMXC9M1b+EjBPRF4Ib0ifUtWjzu8q4HUC463hVg6Ud/uE9iqBN4CB4nPAJlWtDHcgjiuAA6paraoe4E/ARWGOCQBV/a2qTlfV2QSGp0Myfg+W8EPC+XL0t0CZqv5buOM5RkRyRCTduZ1I4D/BzvBGBar6oKoWqmoRgWGAZaoa9t4XgIgkO1+84wyZLCDwMTysVPUT4LCIjHc2XQ6EdVLACW5jgAznOD4GZolIkvP/83IC362FnYjkOr9HAjcTwuctZGva9hdnsfS5QLaIlAMPq+pvwxsVFwNfBj5yxssBfqyqfwljTAD5wPPO7Iko4I+qOqCmQA5AecDrgRxBDPB7VX0nvCF1+Q7wojN8sh+4K8zxAOCMRc8H7gl3LMeo6vsi8iqwicCQyYcMnKtuXxORLMADfFtV60PV0KCflmmMMSY4NqRjjDERwhK+McZECEv4xhgTISzhG2NMhLCEb4wxEcISvgk5EXlERB44zT7pIvKtfornDyKyVUS+3x/t9dD+3P6qCNqfbZmBb9DPwzf9x7lgRVQ1FPVu0oFvAf/ZQ7vRqurri0ZEZBhwkaqO6sUxMU5hK2MGNevhm1MSkSKn1vp/ErhoZYSI/EZENpxYU9+pG/+oiGxy6sdP6OF83xCRvzpX+nb3z8A5Tg31nzs90+Ui8nvgI+fYN5wCZtu7FzETkWYR+akE6vyvE5E8Z/utTp3xLSKyytl9MZDrtHOpiEx1jtkqIq87NVcQkRUi8v+JyErgPhF5znncy0Vkv4jMkcBaDGUi8ly3WBaIyFrnOXjFqaeEiFwlgbr1qwlcTdnTc32/iDzj3D7PiT3phH3eF5FJ3e6vEJHzRWSmiLznFFJ7r9sVuN2PPe6TlnP+Iuf2HRJYK2GziPxfCRTZi3Ye9zbn3zMsn4hMH1JV+7Gfk/4ARQQqWM7qti3T+R1NoB5PiXP/IPAd5/a3gKed248ADwD/C/gzEH+SdrZ1uz+XQEGw0T20m0igvEGWc1+B65zbjwP/6Nz+CChwbqefpJ2twBzn9mPAE87tFcB/dtvvOQJ1foRAEa5G4DwCnaaNwFQgG1hFYM0BgB8CDwEJwGFgrHP8H4G3engOopzjbwI2ECjaduI+3wcedW7nA7ud22lAjHP7CuC1bs/jW93/Hbqda5vzfBQDbwKxzvb/BL4CnE+goirdn0P7Gbw/1sM3wTikxy+S8gUR2UTg8vRJwMRufztWKG4jgWRyzJcJFNT6vKp2BNnuB6p6oNv974rIFmAdgRKyY53tncCxceru7a4BnhORbxB4czqOiLgIJLGVzqbnCdSYP+blEw55UwOZ7yOgUlU/0sDw1nanzVkEnos1TkmNO4FRwAQChbv2OMf3WBjOOdffAf8NrFTVNT3s9kfgVuf2F/i0GqULeEUCK7/9ksC/S7AuJ5Dc1ztxX06g/PJ+YIyI/LuIXEXgTc4MYjaGb4LRcuyGiIwm0Fufoar1znBG96XijiVzH8e/vrYR6AUXAgdEZASBXiXAUwRKw56q3bkEeq4XqmqriKzo1q7HSaTHtauq94rIBQQWVtksIlN78ZiPa/+Ex+bvdvvY/Rin7SWqelv3g5x2g61hMpbACm7De/qjqh4RkVoRKQEW8mm9mp8Ay1X1JmeYZkUPh3s5fhj32PMnwPOq+uCJB4jIFOBK4NsE3mC+GuTjMAOQ9fBNb6URSIRuZ6z8c0Ee9yGB5PRnERmuqodVdarz8xTQBKSe4ngXUO8k+wkEetOnJCLnqOr7qvoQUMMJC0uoqhuoF5FLnU1fJrDi0JlaB1wsIuc67SeJyDgCFUlHi8g5zn639XSw84njSQKfMrJE5JaTtPMSgcViXKr6kbPNRWA5QQh8SujJQZwSyiIynUBteIB3gVvk06qNmSIySkSygShVfQ34Pwys8svmDFjCN72iqlsIJO/twDMEhk2CPXY1gU8HbzvJpPvfagkMhWwTkZ/3cPg7QIyIbCXQmw1mHd6fO182biMwNr6lh33udPbbSuATyGPBPp4TqWo1gWT7B+d864AJqtoO3E3gca8GTrbk3y8JfG+wG/ga8M/HkvAJXiVQRvqP3bY9DvxMRNbQw/CV4zUg0xm2+Saw24l7B/CPBFb12kpgMY58AstyrnD2fw74zCcAM7hYtUxjjIkQ1sM3xpgIYQnfGGMihCV8Y4yJEJbwjTEmQljCN8aYCGEJ3xhjIoQlfGOMiRCW8I0xJkL8/+xYkDQ/zKVjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Plot the result\n",
    "plot(rank_xs,rank_ys,'o-',alpha=0.50,label=\"Pearson regression on ranks\")\n",
    "results = linregress(rank_xs,rank_ys)\n",
    "xlabel(\"rank-transformed x values\")\n",
    "ylabel(\"rank-transformed y values\")\n",
    "legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If that's what Spearman regression sees, it's no wonder it finds a perfect trend!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The flexibility of rank-based methods trades off with statistical power\n",
    "\n",
    "So if rank-based methods are so flexible, why would we ever use a method like Pearson regression that makes strong assumptions about the data?\n",
    "\n",
    "The answer lies in the notion of [statistical power](https://en.wikipedia.org/wiki/Power_of_a_test). Statistical power refers to one's chances of rejecting an untrue hypothesis given some set of data and some particular statistical test. So if you have more statistical power you are more likely to get statistically significant results if your null hypothesis is in fact incorrect. Thus, statistical power is the same idea as *statistical sensitivity* or the *true positive rate* when doing classification tasks, and is the opposite of the False Negative Rate or FNR (indeed in that context statistical power is 1 - FNR).\n",
    "\n",
    "Informally, if we have more statistical power, we're less likely to miss a real pattern in our data. \n",
    "\n",
    "One major reason that Pearson regression is still widely used - even though Spearman regression is more flexible - is that *if* our data is actually linear, Pearson regression will have greater *statistical power* to detect the significance of that linear trend. \n",
    "\n",
    "Let's test this out by doing a new simulation. In this case, the data are definitely linear, but have a lot of noise. This causes each observed data point to vary from the y-value predicted from our linear regression by a random amount. We can simulate this random variance by adding a number drawn from a normal distribution to each y-value.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1a2a0cfb38>"
      ]
     },
     "execution_count": 211,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEKCAYAAAA8QgPpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl0ldW9//H3lxBIJJjIHOaIiECUiAFEEFEckFqoiBXrrXq1clet1fbXZa/a1lqqva3t7b21g5ZWr7aiqKgUtQUnlDohIAEDKIMgBMLcBJAEMuzfH/tJOMSEHE9ypuTzWuus85x9npznCwfz9dnDd5tzDhERkS+qTbwDEBGR5KQEIiIiEVECERGRiCiBiIhIRJRAREQkIkogIiISESUQERGJiBKIiIhERAlEREQi0jbeAURTly5dXP/+/eMdhohIUlm+fPke51zXxs5r0Qmkf//+LFu2LN5hiIgkFTP7NJzz1IUlIiIRUQIREZGIKIGIiEhEWvQYSH0qKiooKiqivLw83qFIFKSlpdG7d29SU1PjHYpIi9fqEkhRUREdO3akf//+mFm8w5Fm5Jxj7969FBUVkZOTE+9wRFq8VteFVV5eTufOnZU8WiAzo3Pnzrq7FImRVpdAACWPFkzfrUjstMoEIiIiTacEkgDuuecefvWrXzX4/rx581izZk0MIxKRpOMcfPoOlO+P2SVb3SD6F7W2uJQFhTvZVlJGr6x0JuZ2Z3B2ZkxjmDdvHpdddhlDhgyJ6XVFJAk4B1vfh4LZsO8TOOt6OH1aTC6tO5DjWFtcyqzFmygtqyA7M43SsgpmLd7E2uLSJn/2fffdx6BBg7jwwgv5+OOPAfjTn/7EiBEjGDZsGFdccQWHDh3inXfeYf78+dx+++3k5eWxcePGes8TkVaoaBm8+F14/adw5DMY+10YennMLq8EchwLCneSmZ5KZnoqbcxqjxcU7mzS5y5fvpw5c+awYsUKnnvuOZYuXQrA1KlTWbp0KStXrmTw4ME8/PDDnHPOOUyePJlf/vKXFBQUMGDAgHrPE5FW6OO/w+H9MOY2uPyPcMoEaJMSs8urC+s4tpWUkZ2Zdkxbx7S2bCspa9Ln/vOf/+Tyyy/nhBNOAGDy5MkAFBYW8sMf/pCSkhIOHjzIJZdcUu/Ph3ueiLQwxatg5ZMw+hbI7AXn3ArtMiAlPr/KlUCOo1dWOqVlFWSmH13VfKC8kl5Z6U3+7Pqmm15//fXMmzePYcOG8eijj/LGG2/U+7PhniciLcSOQih4AnasgvST4LNdPoGkZ8U1LHVhHcfE3O6UllVQWlZBtXO1xxNzuzfpc8eNG8fzzz9PWVkZBw4c4IUXXgDgwIEDZGdnU1FRwezZs2vP79ixIwcOHKh93dB5ItLCOAev/gQW3AGlW2DkDLjiYeh5ZrwjA3QHclyDszOZMS7nmFlYV43o3eRZWMOHD+eqq64iLy+Pfv36ce655wLw05/+lFGjRtGvXz9OP/302qQxffp0brrpJh544AHmzp3b4Hki0kKUbIGsvmAGnU+B7DNg0CRo2z7ekR3DnHPxjiFq8vPzXd0NpdauXcvgwYPjFJHEgr5jSVp7N/rpuFvfh4n/BT1Oj0sYZrbcOZff2Hm6AxERibd9m/wYx5Z3oV0HOPPr0GlAvKNqlBKIiEg8VVXAyz+EqiOQ9zUYMsUnkSSgBCIiEmslW2HdAsi/EVJSYfydcFI/aN8x3pF9IUogIiKxUroNVs2BjYv8gPiAC6DzAOiRG+/IIqIEIiISbUcOwfuzYONr0CYVcqfC0KlxX8fRVFFfB2JmfcxskZmtNbPVZnZb0N7JzF4xs/XB80lBu5nZA2a2wcxWmdnwkM+6Ljh/vZldF+3YRUSapPKwf26b5gsdDp4M0x6G/BuSPnlAbBYSVgLfc84NBs4GvmVmQ4A7gNeccwOB14LXAJcCA4PHDOBB8AkH+DEwChgJ/Lgm6SSbjIwMALZv3860abGpmpnozjnnnHiHINJ8PtsD7/4Bnv2Gv/to0wYu+18YeZNfSd5CRL0LyzlXDBQHxwfMbC3QC5gCjA9Oewx4A/jPoP0vzi9Qec/MsswsOzj3FefcPgAzewWYCDwZ7T9DtPTs2ZO5c+dG9RqVlZW0bVv/13y898JVVVVFSkrTi7e98847Tf4Mkbg7tA8+fMYXOQQYeDFUV/rjNi2v8EdM/0Rm1h84E1gCdA+SS02S6Rac1gvYGvJjRUFbQ+1Ja/PmzeTm+sGzRx99lKlTpzJx4kQGDhzI97///drzXn75ZUaPHs3w4cO58sorOXjwIAAzZ85kxIgR5ObmMmPGDGoWhY4fP5677rqL8847j9/85jfHXPOee+5hxowZXHzxxVx77bVUVVVx++23M2LECM444wz++Mc/AlBdXc3NN9/M0KFDueyyy5g0aVJtsuvfvz8zZ85k7NixPPPMM2zcuJGJEydy1llnce655/LRRx8B8Mwzz5Cbm8uwYcMYN24cAKtXr2bkyJHk5eVxxhlnsH79euDoXZlzjttvv53c3FxOP/10nnrqKQDeeOMNxo8fz7Rp0zjttNO45ppraMmLYCUJHdgBz94IH70IAybA5bNg9Lcg7cR4RxY1MRtEN7MM4FngO865/cfZu7q+N9xx2uteZwa+64u+ffs2Htg/7vh8W865cNqXfP/lKz/+/PunXAgDL4TyUlj0X8e+d+nPG79mAwoKClixYgXt27dn0KBBfPvb3yY9PZ17772XV199lQ4dOvCLX/yCX//619x9993ccsst3H333QB8/etf58UXX+TLX/4yACUlJbz55pv1Xmf58uW89dZbpKenM2vWLDIzM1m6dCmHDx9mzJgxXHzxxSxfvpzNmzfz4YcfsmvXLgYPHswNN9xQ+xlpaWm89dZbAEyYMIGHHnqIgQMHsmTJEm6++WZef/11Zs6cycKFC+nVqxclJSUAPPTQQ9x2221cc801HDlyhKqqqmNie+655ygoKGDlypXs2bOHESNG1CafFStWsHr1anr27MmYMWN4++23GTt2bMR/3yJNVlYCuz+CvmdDxx6Qdw30H+uPW4GYJBAzS8Unj9nOueeC5p1mlu2cKw66qHYF7UVAn5Af7w1sD9rH12l/o+61nHOzgFngS5k04x8j6iZMmEBmpq+zNWTIED799FNKSkpYs2YNY8aMAeDIkSOMHj0agEWLFnH//fdz6NAh9u3bx9ChQ2sTyFVXXdXgdSZPnkx6uq8o/PLLL7Nq1arau4vS0lLWr1/PW2+9xZVXXkmbNm3o0aMH559//jGfUfP5Bw8e5J133uHKK6+sfe/wYT9wOGbMGK6//nq++tWvMnXqVABGjx7NfffdR1FREVOnTmXgwIHHfO5bb73F1VdfTUpKCt27d+e8885j6dKlnHjiiYwcOZLevXsDkJeXx+bNm5VAJD7K98Pq52HtfF/w8Kt/gfYZMdsJMFFEPYGYv9V4GFjrnPt1yFvzgeuAnwfPfwtpv8XM5uAHzEuDJLMQ+FnIwPnFwJ1NDvB4dwxt2x///bTMJt1x1NW+/dFCaSkpKVRWVuKc46KLLuLJJ48d6ikvL+fmm29m2bJl9OnTh3vuuYfy8vLa9zt0aHgla+h7zjl++9vffm5PkZdeeum4sdZ8RnV1NVlZWRQUFHzunIceeoglS5bw0ksvkZeXR0FBAV/72tcYNWoUL730Epdccgl//vOfueCCC46JpyH1/f2IxNThg7BmHqz5G1SUQc44GHa1Tx6tUCzGQMYAXwcuMLOC4DEJnzguMrP1wEXBa4C/A58AG4A/ATcDBIPnPwWWBo+ZNQPqLdnZZ5/N22+/zYYNGwA4dOgQ69atq00WXbp04eDBgxEPxl9yySU8+OCDVFRUALBu3To+++wzxo4dy7PPPkt1dTU7d+5scM+RE088kZycHJ555hnAJ4CVK1cCsHHjRkaNGsXMmTPp0qULW7du5ZNPPuHkk0/m1ltvZfLkyaxateqYzxs3bhxPPfUUVVVV7N69m8WLFzNy5MiI/mwiza68FFY9DT2Hw5TfwXnfh6w+jf9cCxWLWVhvUf/4BcCEes53wLca+KxHgEeaL7rE17VrVx599FGuvvrq2q6he++9l1NPPZWbbrqJ008/nf79+zNixIiIPv8b3/gGmzdvZvjw4Tjn6Nq1K/PmzeOKK67gtddeIzc3l1NPPZVRo0bVdq/VNXv2bL75zW9y7733UlFRwfTp0xk2bBi3334769evxznHhAkTGDZsGD//+c95/PHHSU1NpUePHrVjODUuv/xy3n33XYYNG4aZcf/999OjR4/agXmRmDpyyA+K79/m9xvP7AXTHoEOXeIdWUJQOXdp0MGDB8nIyGDv3r2MHDmSt99+mx49En9wUN+xNFlFuU8chc/C4QPQZySMvytuW8fGmsq5S5NddtlllJSUcOTIEX70ox8lRfIQabKdq2HRz3x3Va+z/MyqrqfGO6qEpAQiDdJe69JqVB6Bsn1++m1WX+g6CE6/ErrpTvZ4WmUCcc5xnHUoksRacpesREFVBax/GVbO8SVGvvwbX1J9wt2N/6y0vgSSlpbG3r176dy5s5JIC+OcY+/evaSlpcU7FEl0VZWw4RVY9ZSvW9VtCJx5jd+DXMLW6hJI7969KSoqYvfu3fEORaIgLS2tdrGhSIM+WQTv/t53VZ1zK/Q8U8kjAq0ugaSmppKTkxPvMEQklqqr4JM3IKWdL1V08vlwQie/nkOJI2KtLoGISCtSXQ2bF0PBE7B/O/Qe4RNISls/w0qaRAlERFqm7StgySwo3er3Gz//Lug7Ot5RtShKICLScjjnZ1a1becHygHG3wH9xqirKgqUQEQk+TkHW9+Hgsd919RZ10PvfH/cAjdyShRKICKSvJyDbcthxV9h70a/ELDTyf49M911RJkSiIgkr+WP+npVGd1gzG0w4AJo0/QtliU8SiAiklyKV0JGd3+3cfJ46JjtdwltJYUOE4n+xkUkOewohBWPw85Cv+X02d+ETjn+IXGhBCIiiW3XWp84ildCehaMnAGnTox3VIISiIgkuk2L4V+bYcSNMGiS32paEoISiIgklr0b/R3H0Msh+wzI+xoMvw5SVSQz0SiBiEhi2PeJLzmy5T1olwEnn+fb23eMb1zSICUQEYm/d38PH/8DUk/wOwAOmQLtToh3VNIIJRARiY/SIujY068UP6k/DJsOQ74C7TPiHZmESQlERGKrdBusfAI+eRPGfhdOmeCn5UrSUQIRkdjYXwwrn4SNr/t9OXKnqqR6klMCEZHocw4W3ev35BgyBU6f5vcgl6SmBCIi0XFwN6x+Hs78Nz8gPua7fhfAEzrFOzJpJkogItK8Du2DVU/BugX+da+zoPdZ0OWU+MYlzS7qhfLN7BEz22VmhSFt95jZNjMrCB6TQt6708w2mNnHZnZJSPvEoG2Dmd0R7bhF5AuqroL3/wRzb/BTck+5EKb+yScPaZFicQfyKPA74C912v/HOfer0AYzGwJMB4YCPYFXzezU4O3fAxcBRcBSM5vvnFsTzcBFJAyVh315kTYpULLFV8gdNt1Xy5UWLeoJxDm32Mz6h3n6FGCOc+4wsMnMNgAjg/c2OOc+ATCzOcG5SiAi8VJeCoXPwbqFMPm3kNEVLvyJdgBsReI5BnKLmV0LLAO+55z7F9ALeC/knKKgDWBrnfZR9X2omc0AZgD07du3uWMWkcMH/OD4mvlQWQ455x19T8mjVYnXt/0gMADIA4qB/w7a69t/0h2n/fONzs1yzuU75/K7du3aHLGKSI0jh+C5GbDqab/n+JTfw3m3+7sPaXXicgfinNtZc2xmfwJeDF4WAX1CTu0NbA+OG2oXkWg6cgiKlvrihu1O8NNyuw3RRk4SnwRiZtnOueLg5eVAzQyt+cATZvZr/CD6QOB9/B3IQDPLAbbhB9q/FtuoRVqZijJY+6Lfc/zIQeh8CmT2UtkRqRX1BGJmTwLjgS5mVgT8GBhvZnn4bqjNwH8AOOdWm9nT+MHxSuBbzrmq4HNuARYCKcAjzrnV0Y5dpFWqPAwfveQTR3mpX8eRd41PHiIhzLl6hxJahPz8fLds2bJ4hyGSXA4f8Gs5ug6CvH+DbqfFOyKJMTNb7pzLb+w8rUQXae2qKvxU3G3LYcLdfgOnrzwEHTrHOzJJcEogIq1VVSWsf9mXHTm0F7oP9XcfaScqeUhYlEBEWqOSrfDK3fDZbuh6mt+XI3sYWH0z5kXqpwQi0lpUV8GBYsjsDR2z/ayq0bdAr+FKHBIRJRCRlq66Gja96TdzqiiDKx6Gtu3ggh/EOzJJckogIi2Vc7D5n1DwhN9//KR+cNa/Q0pqvCOTFkIJRKSlKl4Jb94PWX1h/J3Q7xx1VUmzUgIRaSmcg61L/Iyq077kB8UvvMcvBFTikChQAhFJds5B0TJY8VfY9wl0OhkGTfJJo3eja8FEIqYEIpLM9qyH9/7gnzO6w9jvwMnn645DYkIJRCTZOAdVR4JdANv6elXn3AoDLoAU/SctsaN/bSLJZMeHsOJxOKGL34ejUw5M/bM2cpK4UAIRSQY71/jEsWMVpJ8EOeOOvqfkIXGiBCKS6Na+CEsegrRMGPENGHSp774SiTMlEJFEtGc9WBvoPMCv36g6DKddpsQhCUUJRCSR7N0IBbNh6/vQewRc+GM4oRPkXhHvyEQ+RwlEJBH8azOsmA1b3oV2Hfy+44MnxzsqkeNSAhFJBNuWQ3EBDLsahkyB9hnxjkikUUogIvFQug1WPuHLjAy4wI9vDLzY7wYokiSUQERiaf92X1Z94yJIaQddTvXtbdtrgFySjhKISKwUPOnvOtqkwtDL/cB4ela8oxKJ2BdKIGbWBshwzu2PUjwiLcvB3X5QvN0Jfj+O074Mp0/zM6tEklyjS1jN7AkzO9HMOgBrgI/N7PbohyaSxD7bC+/+AZ77Bqx9wbf1OwdGzVDykBYjnDuQIc65/WZ2DfB34D+B5cAvoxqZSDI6tA8+fAY+/ge4aj8wPuCCeEclEhXhJJBUM0sFvgL8zjlXYWYuynGJJKd3HvBTck+5EM6YDh27xzsikagJpwrbH4HNQAdgsZn1A8IeAzGzR8xsl5kVhrR1MrNXzGx98HxS0G5m9oCZbTCzVWY2PORnrgvOX29m14V7fZGoKi+FZf/nxzoA8m+Ay2fBmNuUPKTFazSBOOcecM71cs5Nct6nwPlf4BqPAhPrtN0BvOacGwi8FrwGuBQYGDxmAA+CTzjAj4FRwEjgxzVJRyQuDh+A5Y/B3Buh8Fm/CBD8/uMnZsc3NpEYabQLy8y6Az8DejrnLjWzIcBo4OFwLuCcW2xm/es0TwHGB8ePAW/gx1amAH9xzjngPTPLMrPs4NxXnHP7gphewSelJ8OJQaRZrZzjk0ZFGfQ/F/Ku9olDpJUJZwzkUeD/gB8Er9cBTxFmAmlAd+dcMYBzrtjMugXtvYCtIecVBW0NtYvERuXhowv9DuyAnmf6siOdcuIbl0gchTMG0sU59zRQDeCcqwSqohRPfRs5u+O0f/4DzGaY2TIzW7Z79+5mDU5aoYoyWPU0PH0d7F7n2865Fc6/S8lDWr1wEshnZtaZ4Be2mZ0NlDbxujuDrimC511BexHQJ+S83sD247R/jnNulnMu3zmX37Vr1yaGKa1W5WH4cK4f4/jgL9B10NE7EO0AKAKE14X1/4D5wAAzexvoCkxr4nXnA9cBPw+e/xbSfouZzcEPmJcGXVwLgZ+FDJxfDNzZxBhE6lddDfO/7etW9TzTl1bvOijeUYkknEYTiHPuAzM7DxiE70r62DlXEe4FzOxJ/CB4FzMrws+m+jnwtJndCGwBrgxO/zswCdgAHAL+PYhhn5n9FFganDezZkBdpFlUHoHN//SL/tq0gWHTIaMHdB8S78hEEpb5CU/HOcHs2vranXN/iUpEzSg/P98tW7Ys3mFIIquqgPUvw6qn/CryS+6D7GHxjkokrsxsuXMuv7HzwunCGhFynAZMAD4AEj6BiDSoqhI2vAqr5sBne6DraXDu96DHGfGOTCRphNOF9e3Q12aWCfw1ahGJxMqqpyC9k59V1fNMsPom+4lIQyLZD+QQfqW4SPKoroZPFvnuqot+Cm3bwZf+G9JPUuIQiVA4K9Ff4OiaizbAEODpaAYl0myqq/3geMETsH8bnNQfDu315UZUVl2kScK5A/lVyHEl8KlzrihK8Yg0n7ISWHgXlGyBrH5+8V/f0brjkBZtbXEpCwp3sq2kjF5Z6UzM7c7g7MyoXCucMZA3o3JlkWhwDkq3+tpUaZnQeYAvOdJ/rBJHFMTyl5U0bm1xKbMWbyIzPZXszDRKyyqYtXgTM8blROV7aXBJrZkdMLP99TwOmJm2tJXE4hxsWQIv3AYvfMfffZj5mVU55yp5REHNL6vSsopjflmtLW5qoQqJ1ILCnWSmp5KZnkobs9rjBYU7o3K9Bu9AnHMdo3JFkebknN/AacXjsHcDZHSH0d+C9vrnG22hv6yA2ucFhTt1FxIn20rKyM5MO6atY1pbtpWUReV6Yc/CCirm1kbmnNsSlYhEvoj92+DVeyCjm9/E6eTzISWSyYXyRcX6l5U0rldWOqVlFbXJHOBAeSW9stKjcr1Gq8KZ2WQzWw9sAt7E7074j6hEIxKO4pW+Qi5AZm+4aKbfBXDgRUoeMdQrK50D5ZXHtEXzl5U0bmJud0rLKigtq6DaudrjibnR2R0znLKiPwXOBtY553LwK9Hfjko0IsezoxAW3AkLfwAf/92XWgfoNVyJIw5i/ctKGjc4O5MZ43LITE+luLSczPTUqA2gQ3hdWBXOub1m1sbM2jjnFpnZL6ISjUh9SrfBe3/wdx5pmTByBpw60S8GlLip+WUVOgvrqhG9Nf4RZ4OzM2P2HYSTQErMLANYDMw2s1349SAi0VWzC2Bqui+tPuJGGDTp6L4cEnex/GUliSecBDIFKAO+C1wDZAIzoxmUtHJ7N/pZVUc+g0t/4VeMX/GwNnISSTDhJJAZwDPB6vPHohxPi6UFV2HY94kvObLlPWiXAblTwVWDpSh5iCSgcBLIicBCM9sHzAHmOueisyqlhYr16tCktOU9eP1e312Vdw0MmQztOsQ7KhE5jnBKmfwE+ImZnQFcBbxpZkXOuQujHl0LoQVXDSjZCuWl0CP36Naxp31JiwBFksQXmfu4C9gB7AW6RSeclkkLruoo3eY3ctq4yFfHnfxbPzA+bHq8IxORLyCccu7fxN95dAXmAjc559ZEO7CWJNarQxPWgR2wcg5sfA3apPoxjqFTVadKJEmFcwfSD/iOc64g2sG0VBNzuzNr8SbA33kcKK+ktKyCq0b0jnNkMbZnPWx6EwZPhtwrtB+HSJIz51zjZyWp/Px8t2zZsiZ9RnPNnmqVs7A+2+NLjmR0g9On+cKHZf9S4hBJcGa23DmX39h5qv9wHM05e6pVLbg6tA8+fMaXGwEY8hX/bKbkIdKCKIEch2ZPRWDdQljyEFRX+eKGZ1zl70BEpMUJZxD9FmC2c+5fMYgnoWj2VJjKSvxzepbfOjZnnN8FsGOP+MYlIlEVzh1ID2CpmX0APAIsdC154CSEZk81onw/rH4e1s6Hk8fDOd+Gbqf5h4i0eI3Wh3DO/RAYCDwMXA+sN7OfmdmApl7czDab2YdmVmBmy4K2Tmb2ipmtD55PCtrNzB4wsw1mtsrMhjf1+o1RueoGHD7ga1U9e6Mf6+gzCoZMiXdUIhJjYRUYCu44dgSPSuAkYK6Z3d8MMZzvnMsLGfG/A3jNOTcQeC14DXApPpENxNfnerAZrn1csa6tnzQ++Ktfz9FzOEz5HZz3fcjqG++oRCTGwhkDuRW4DtgD/Bm43TlXYWZtgPXA95s5pinA+OD4MeAN4D+D9r8Eyew9M8sys2znXHEzX/8YiTp7KqbTgo8cgo9e9OVGugyE06+EQROh08nRuZ6IJIVwxkC6AFOdc5+GNjrnqs3ssiZe3wEvm5kD/uicmwV0r0kKzrniYC92gF7A1pCfLQraoppAElHMijNWlPvEUfis77aqrvIJJKMrvjCBiLRm4RRTvPs4761t4vXHOOe2B0niFTP76Djn1lfv4nOD+WY2A9/FRd++LbNbJSbTiz/+hx/nKC/1W8bm/Rt0PbV5PltEWoS4rgNxzm0PnneZ2fPASGBnTdeUmWXjiziCv+PoE/LjvYHt9XzmLGAW+JXo0Yw/XqI2vbjyCKSk+gV/h/ZBpxxfWr3b4KZ9roi0SHHbpcfMOphZx5pj4GKgEJiPH3MheP5bcDwfuDaYjXU2UBrt8Y9E1SsrnQPlx+4q3KTpxVUV8NFL8Nw3YMu7vm3Y1XDxvUoeItKgeN6BdAeeN1+JtS3whHNugZktBZ42sxuBLcCVwfl/ByYBG4BDwL/HPuTE0GzFGasqYcOrvrT6Z3ug2xBID0qNaAdAEWmEiikmqWaZhbXgLtixCrqc6jdz6nmmSquLiIoptnQRTS+uroJNi6HfOX4Dp6Ff8Xty9DpLiUNEvjAlkNaguho2L4aCJ2D/dnBVcMqF0GdkvCMTkSSmBNKSOQeb3/KJo3QrnNQPzr8L+o6Od2Qi0gIogbR0a+b55/P+E/qPVVeViDQbJZCWxDnY+j4UzoXzf+DLq5//A0jL0qwqEWl2SiAtgXOwbblfOb53A3TMhs92+wSiHQBFJEqUQJJdVQUsuBN2f+R3/htzGwy4ANqkxDsyEWnhlECSVckWX0I9JRW6DvKzqk65EFL0lYpIbOi3TbLZUQgFs2HHhzDl935m1cib4h2ViLRCSiDJYtdaP8ZRvNKPbYyc4cc6RETiRAkkGRw+AAt/AKnpMOJGGDTJryQXEYkjJZBEtWcDbP4nnHU9tO8IF97ja1alpjXygyIisaEE0oyapcDh3o1+5fjWJdAuA067zO8AmH1GdIJOEjHdwldEwqLVZc2kZpvZ0rKKY7aZXVtcGt4HHNoHi34GL9wGOwt9ddxpjwTbx7ZuTf67FZGo0B1IM4l4m9nKw348o10H+NdmGDYdhnwF2mfEIOrkEJMtfEXkC1MCaSZfeJvZ0m2w8knYsw6m/MEnka88pJIj9YjaFr4i0iQesNQtAAAM9UlEQVRKIM2kV1Y6pWUVtf93DA1sM7u/2CeOja9DSjsYfBlUV/oFgEoe9Qr771ZEYkq/sZrJxNzulJZVUFpWQbVztccTc7sfPWnPenh+hp9dNWQKTHsY8m/QzKpGhPV3KyIxpy1tm1G9M4U6HoGST6F3vi96+OEzvuSIihx+IZqFJRI74W5pqwQSLYf2waqnYN0CPx33ysdUp0pEkoL2RI+Xsn/Bh3Pho5fAVcPAi+CMq5Q8RKTF0W+15nZgB6ydDwMm+Cm5HXvEOyIRkahQAmmq8v2w+jm/L8fIm6DbYN9dpTEOEWnhlEAidfgArJ4Ha/4GleV+Eyfn/J7jSh4i0googURiyxL456+gogz6j4W8r/nNnUREWhElkHAdOQRHPvO1qTrlQM/hfoyjU068IxMRiYukW0hoZhPN7GMz22Bmd0T9ghVlsOoZmHsDvPOAb8voBuffqeQhIq1aUt2BmFkK8HvgIqAIWGpm851za5r9YpWH4aMX4cNn4fB+6HUW5F3T7JcREUlWSZVAgJHABufcJwBmNgeYAjR/AlkzHz54DHrm+cTRbXCzX0JEJJklWwLpBWwNeV0EjAo9wcxmADMA+vZtwsD2aZOg+xDWVvdmwcqdbCtZqRIaIiIhkm0MxOppO6YWi3NulnMu3zmX37VrEzZjateBtdW9tZGRiEgDki2BFAF9Ql73BrZH62KhGxm1Mas9XlC4M1qXFBFJGsmWQJYCA80sx8zaAdOB+dG62LaSMjqmHdvLp42MRES8pEogzrlK4BZgIbAWeNo5tzpa1+uVlc6B8spj2rSRkYiIl1QJBMA593fn3KnOuQHOufuieS1tZCQi0rCkSyCxNDg7kxnjcshMT6W4tJzM9FRmjMvRLCwREZJvGm/MDc7OVMIQEamHEkg9arZPXb29lP3llWSmt2VIdqbWgIiIhFAXVh1ri0uZtXgTm3YfpGhfGfvLKtiy9xCb9xzUGhARkRBKIHXUrP3YceAw7VPbkJmeSlpqCjv2H9YaEBGREEogddSs/ThYXkn7tv6vp33bNuwvr9AaEBGREEogddSs/chIa8vhymoADldWc2JaqtaAiIiEUAKpo2btR4+O7TlcUU1pWQXlFVX0OLG91oCIiITQLKw6atZ+LCjcyaGKqtpZWP27ZGgWlohICCWQemjth4hI49SFJSIiEVECERGRiCiBiIhIRJRAREQkIkogIiISESUQERGJiKbx1qOmGu+2kjJ6ZaVr/YeISD10B1JHTTXe0rIKsjPTKC2rUBVeEZF6KIHUUVONNzM9lTZmtceqwisiciwlkDpqqvGGUhVeEZHPUwKpo6YabyhV4RUR+TwlkDpqqvGWllVQ7VztsarwiogcSwmkjppqvJnpqRSXlpOZnsqMcTmahSUiUoem8dZD1XhFRBqnOxAREYlIXBKImd1jZtvMrCB4TAp5704z22BmH5vZJSHtE4O2DWZ2RzziFhGRo+LZhfU/zrlfhTaY2RBgOjAU6Am8amanBm//HrgIKAKWmtl859yaWAYsIiJHJdoYyBRgjnPuMLDJzDYAI4P3NjjnPgEwsznBuUogIiJxEs8xkFvMbJWZPWJmJwVtvYCtIecUBW0NtYuISJxELYGY2atmVljPYwrwIDAAyAOKgf+u+bF6Psodp72+684ws2Vmtmz37t3N8CcREZH6RK0Lyzl3YTjnmdmfgBeDl0VAn5C3ewPbg+OG2utedxYwCyA/P7/eJCMiIk0Xr1lY2SEvLwcKg+P5wHQza29mOcBA4H1gKTDQzHLMrB1+oH1+LGMWEZFjxWsQ/X4zy8N3Q20G/gPAObfazJ7GD45XAt9yzlUBmNktwEIgBXjEObc6HoFHQvuLiEhLZM613F6e/Px8t2zZsrjGULO/SGZ6Kh3T2nKgvJLSsgqVRxGRhGVmy51z+Y2dl2jTeFuc0P1FgNrnBYU7I04guqMRkUSgUiZR1tz7i2jHRBFJFEogUdbc+4tox0QRSRRKIFHW3PuLaMdEEUkUSiBR1tz7i2jHRBFJFBpEj4Hm3F9kYm53Zi3eBHDMrK6rRvRuls8XEQmX7kCSjHZMFJFEoTuQJKQdE0UkEegOREREIqIEIiIiEVECERGRiCiBiIhIRJRAREQkIi26Gq+Z7QY+bcJHdAH2NFM40ZDo8UHix5jo8YFibA6JHh8kVoz9nHNdGzupRSeQpjKzZeGUNI6XRI8PEj/GRI8PFGNzSPT4IDlirEtdWCIiEhElEBERiYgSyPHNincAjUj0+CDxY0z0+EAxNodEjw+SI8ZjaAxEREQiojsQERGJiBJIPcxsopl9bGYbzOyOOMbxiJntMrPCkLZOZvaKma0Pnk8K2s3MHghiXmVmw2MQXx8zW2Rma81stZndloAxppnZ+2a2MojxJ0F7jpktCWJ8yszaBe3tg9cbgvf7RzvG4LopZrbCzF5M0Pg2m9mHZlZgZsuCtkT6nrPMbK6ZfRT8exydYPENCv7uah77zew7iRRjRJxzeoQ8gBRgI3Ay0A5YCQyJUyzjgOFAYUjb/cAdwfEdwC+C40nAPwADzgaWxCC+bGB4cNwRWAcMSbAYDcgIjlOBJcG1nwamB+0PAd8Mjm8GHgqOpwNPxei7/n/AE8CLwetEi28z0KVOWyJ9z48B3wiO2wFZiRRfnVhTgB1Av0SNMew/S7wDSLQHMBpYGPL6TuDOOMbTv04C+RjIDo6zgY+D4z8CV9d3Xgxj/RtwUaLGCJwAfACMwi/Yalv3OwcWAqOD47bBeRbluHoDrwEXAC8GvzQSJr7gWvUlkIT4noETgU11/x4SJb564r0YeDuRYwz3oS6sz+sFbA15XRS0JYruzrligOC5W9Ae17iDrpQz8f+Hn1AxBt1DBcAu4BX8HWaJc65mb+DQOGpjDN4vBTpHOcT/Bb4PVAevOydYfAAOeNnMlpvZjKAtUb7nk4HdwP8F3YB/NrMOCRRfXdOBJ4PjRI0xLEogn2f1tCXDVLW4xW1mGcCzwHecc/uPd2o9bVGP0TlX5ZzLw/+f/khg8HHiiGmMZnYZsMs5tzy0+TgxxOt7HuOcGw5cCnzLzMYd59xYx9gW39X7oHPuTOAzfHdQQ+L530o7YDLwTGOn1tOWcL+HlEA+rwjoE/K6N7A9TrHUZ6eZZQMEz7uC9rjEbWap+OQx2zn3XCLGWMM5VwK8ge9TzjKzmh05Q+OojTF4PxPYF8WwxgCTzWwzMAffjfW/CRQfAM657cHzLuB5fCJOlO+5CChyzi0JXs/FJ5REiS/UpcAHzrmdwetEjDFsSiCftxQYGMyCaYe/3Zwf55hCzQeuC46vw4871LRfG8zeOBsorbk1jhYzM+BhYK1z7tcJGmNXM8sKjtOBC4G1wCJgWgMx1sQ+DXjdBZ3Q0eCcu9M519s51x//b+1159w1iRIfgJl1MLOONcf4PvxCEuR7ds7tALaa2aCgaQKwJlHiq+NqjnZf1cSSaDGGL96DMIn4wM+AWIfvK/9BHON4EigGKvD/R3Ijvr/7NWB98NwpONeA3wcxfwjkxyC+sfjb6lVAQfCYlGAxngGsCGIsBO4O2k8G3gc24LsT2gftacHrDcH7J8fw+x7P0VlYCRNfEMvK4LG65r+JBPue84Blwfc8DzgpkeILrnsCsBfIDGlLqBi/6EMr0UVEJCLqwhIRkYgogYiISESUQEREJCJKICIiEhElEBERiYgSiEgcmdnBeMcgEiklEBERiYgSiEgYzGxEsC9DWrAye7WZ5dY55xdmdnPI63vM7HtmlmFmr5nZB+b31JhSz+ePt2AvkOD178zs+uD4LDN7MyhkuDCk9MWtZrYmiGtO1P7wIg1o2/gpIuKcW2pm84F7gXTgcedcYZ3T5uDrWP0heP1VYCJQDlzunNtvZl2A98xsvgtjFW9Qa+y3wBTn3G4zuwq4D7gBXzAwxzl3uKZci0gsKYGIhG8mvlZaOXBr3TedcyvMrJuZ9QS6Av9yzm0JksDPggq21fiy3N3xmwo1ZhCQC7ziS4+Rgi9vA75sx2wzm4cv3yESU0ogIuHrBGTgdzZMw5cNr2suvshhD/wdCcA1+IRylnOuIqi8m1bn5yo5tku55n0DVjvnRtdzrS/hd62cDPzIzIa6o3uIiESdxkBEwjcL+BEwG/hFA+fMwVfVnYZPJuBLru8Kksf5+K1M6/oUGGJ+z/NMfEVZ8DvRdTWz0eC7tMxsqJm1Afo45xbhN6PKwic3kZjRHYhIGMzsWqDSOfeEmaUA75jZBc6510PPc86tDkqfb3NHy2/PBl4ws2X4isUf1f1859xWM3sa3y21Hl9BGOfcETObBjwQJJa2+HGWdcDjQZsB/+P8ficiMaNqvCIiEhF1YYmISESUQEREJCJKICIiEhElEBERiYgSiIiIREQJREREIqIEIiIiEVECERGRiPx/+9r2ZwK/ut0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "from numpy import array\n",
    "\n",
    "def simulate_noisy_linear_data(xs,slope=1,intercept=0,variance=18):\n",
    "    \"\"\"Return a list of x and y values for a noisy linear regression\n",
    "    xs -- x values to calculate\n",
    "    slope -- the slope (m) of the line in the equation y=mx+b\n",
    "    intercept -- the intercept (b) of the line in the equation y=mx+b\n",
    "    variance -- the variance of random deviations from the line    \n",
    "    \"\"\"  \n",
    "    ys = []\n",
    "    n_xs = len(xs)\n",
    "    xs = array([float(x) for x in xs])\n",
    "    ys = xs*slope + intercept\n",
    "        \n",
    "    noise_mean = 0.0 #the error distribution centers on 0 error\n",
    "    random_distribution = norm(loc=noise_mean,scale=variance)\n",
    "        \n",
    "    noise = random_distribution.rvs(size = n_xs)\n",
    "    ys += noise     \n",
    "    return xs,ys\n",
    "\n",
    "x_data = [0,1,5,7,8,35,51,53,100,101,250,500,750]\n",
    "xs,ys = simulate_noisy_linear_data(xs= x_data,slope=3,intercept=7,variance=600)\n",
    "plot(xs,ys,'o',alpha=0.5,label=\"data\")\n",
    "xs,true_ys = simulate_noisy_linear_data(xs= x_data,slope=3,intercept=7,variance=0)\n",
    "plot(xs,true_ys,'--',alpha=0.75,label=\"linear regression\")\n",
    "xlabel(\"x values\")\n",
    "ylabel(\"y values\")\n",
    "legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Spearman r =  0.5549450549450549\n",
      "Spearman p= 0.04900449559833098\n",
      "Pearson r =  0.574394575075116\n",
      "Pearson p= 0.04006154239129748\n"
     ]
    }
   ],
   "source": [
    "#Now let's fit Pearson and Spearman regressions\n",
    "from scipy.stats import spearmanr,pearsonr\n",
    "spearman_rho,p = spearmanr(xs,ys)\n",
    "print(\"Spearman r = \",spearman_rho)\n",
    "print(\"Spearman p=\",p)\n",
    "\n",
    "pearson_r,p = pearsonr(xs,ys)\n",
    "print(\"Pearson r = \",pearson_r)\n",
    "print(\"Pearson p=\",p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exact result may vary each time you run the simulation, but in general the Pearson correlation will tend to be stronger (on average across many simulations) than the Spearman correlation. What should we take away from these results?\n",
    "\n",
    "We showed earlier in this section that if, in truth, the trend is [monotonic](https://en.wikipedia.org/wiki/Monotonic_function) (always either flat or headed in one direction) but **not** linear, then a Spearman regression has more power than a Pearson regression.\n",
    "\n",
    "However, our second simulation shows that if a trend actually is linear, then Pearson correlation will have more power on average than Spearman correlation.\n",
    "\n",
    "> Even though Spearman regression is more flexible than Pearson regression, Pearson regression has slightly more statistical power if the data are really linear\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantifying the difference in power between Spearman and Pearson regression using Monte Carlo\n",
    "\n",
    "One major advantage of studying statistics after learning some coding techniques is that you can using simulations to explore the behavior of different statistical tests under circumstances you care about. Let's try this with Spearman and Pearson regression.\n",
    "\n",
    "There are a variety of questions we could address by simulation: \n",
    "\n",
    "- How much do outliers impact the statistical power of Pearson and Spearman regression?\n",
    "- How many samples do I need to detect a true significant trend with Spearman or Pearson regression, given either a linear or non-linear trend and at different levels of variance?\n",
    "- If errors are **not** normally distributed, how often is the result of Pearson regression wrong?\n",
    "\n",
    "We'll explore some of these more in the **Exercises**. For now let's focus on exploring further the statistical power of each of these methods. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running simulations for variance 20\n",
      "Running simulations for variance 40\n",
      "Running simulations for variance 80\n",
      "Running simulations for variance 160\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Test</th>\n",
       "      <th>Variance</th>\n",
       "      <th>Sample Size</th>\n",
       "      <th>Statistical Power</th>\n",
       "      <th>Monte Carlo Trials</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>20</td>\n",
       "      <td>40</td>\n",
       "      <td>0.956667</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>20</td>\n",
       "      <td>40</td>\n",
       "      <td>0.933333</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>20</td>\n",
       "      <td>60</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>20</td>\n",
       "      <td>60</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>20</td>\n",
       "      <td>80</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>20</td>\n",
       "      <td>80</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>20</td>\n",
       "      <td>100</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>20</td>\n",
       "      <td>100</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>40</td>\n",
       "      <td>40</td>\n",
       "      <td>0.456667</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>40</td>\n",
       "      <td>40</td>\n",
       "      <td>0.436667</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>40</td>\n",
       "      <td>60</td>\n",
       "      <td>0.913333</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>40</td>\n",
       "      <td>60</td>\n",
       "      <td>0.913333</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>40</td>\n",
       "      <td>80</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>40</td>\n",
       "      <td>80</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>40</td>\n",
       "      <td>100</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>40</td>\n",
       "      <td>100</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>80</td>\n",
       "      <td>40</td>\n",
       "      <td>0.113333</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>80</td>\n",
       "      <td>40</td>\n",
       "      <td>0.106667</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>80</td>\n",
       "      <td>60</td>\n",
       "      <td>0.370000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>80</td>\n",
       "      <td>60</td>\n",
       "      <td>0.356667</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>80</td>\n",
       "      <td>80</td>\n",
       "      <td>0.720000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>80</td>\n",
       "      <td>80</td>\n",
       "      <td>0.680000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>80</td>\n",
       "      <td>100</td>\n",
       "      <td>0.933333</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>80</td>\n",
       "      <td>100</td>\n",
       "      <td>0.916667</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>160</td>\n",
       "      <td>40</td>\n",
       "      <td>0.080000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>160</td>\n",
       "      <td>40</td>\n",
       "      <td>0.076667</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>160</td>\n",
       "      <td>60</td>\n",
       "      <td>0.126667</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>160</td>\n",
       "      <td>60</td>\n",
       "      <td>0.123333</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>160</td>\n",
       "      <td>80</td>\n",
       "      <td>0.250000</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>160</td>\n",
       "      <td>80</td>\n",
       "      <td>0.243333</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>Pearson</td>\n",
       "      <td>160</td>\n",
       "      <td>100</td>\n",
       "      <td>0.453333</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31</th>\n",
       "      <td>Spearman</td>\n",
       "      <td>160</td>\n",
       "      <td>100</td>\n",
       "      <td>0.426667</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        Test  Variance  Sample Size  Statistical Power  Monte Carlo Trials\n",
       "0    Pearson        20           40           0.956667                 300\n",
       "1   Spearman        20           40           0.933333                 300\n",
       "2    Pearson        20           60           1.000000                 300\n",
       "3   Spearman        20           60           1.000000                 300\n",
       "4    Pearson        20           80           1.000000                 300\n",
       "5   Spearman        20           80           1.000000                 300\n",
       "6    Pearson        20          100           1.000000                 300\n",
       "7   Spearman        20          100           1.000000                 300\n",
       "8    Pearson        40           40           0.456667                 300\n",
       "9   Spearman        40           40           0.436667                 300\n",
       "10   Pearson        40           60           0.913333                 300\n",
       "11  Spearman        40           60           0.913333                 300\n",
       "12   Pearson        40           80           1.000000                 300\n",
       "13  Spearman        40           80           1.000000                 300\n",
       "14   Pearson        40          100           1.000000                 300\n",
       "15  Spearman        40          100           1.000000                 300\n",
       "16   Pearson        80           40           0.113333                 300\n",
       "17  Spearman        80           40           0.106667                 300\n",
       "18   Pearson        80           60           0.370000                 300\n",
       "19  Spearman        80           60           0.356667                 300\n",
       "20   Pearson        80           80           0.720000                 300\n",
       "21  Spearman        80           80           0.680000                 300\n",
       "22   Pearson        80          100           0.933333                 300\n",
       "23  Spearman        80          100           0.916667                 300\n",
       "24   Pearson       160           40           0.080000                 300\n",
       "25  Spearman       160           40           0.076667                 300\n",
       "26   Pearson       160           60           0.126667                 300\n",
       "27  Spearman       160           60           0.123333                 300\n",
       "28   Pearson       160           80           0.250000                 300\n",
       "29  Spearman       160           80           0.243333                 300\n",
       "30   Pearson       160          100           0.453333                 300\n",
       "31  Spearman       160          100           0.426667                 300"
      ]
     },
     "execution_count": 213,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsvXmY1OWVt3+fqq7qfd9ZmgZEEREFW1HAFRfc992o0YxJJpmZJO87v9/MvJmMZpKZSSaZN+vEcWYSiUtc4zJGjXEBUVEBQRRwQbppGuhuet+XqjrvH6e66b0LpLp6ee7r4qKr61vVB61+Pt/nPOd8jqgqDofD4XAAeGIdgMPhcDjGD04UHA6Hw9GLEwWHw+Fw9OJEweFwOBy9OFFwOBwORy9OFBwOh8PRixMFh8PhcPTiRMHhcDgcvThRcDgcDkcvcbEO4FDJycnR4uLiWIfhcDgcE4pNmzbVqGruaNdNOFEoLi5m48aNsQ7D4XA4JhQisjuS61z6yOFwOBy9OFFwOBwORy9OFBwOh8PRixMFh8PhcPQy6kGziJQApwPTgHbgQ+BlVa2LcmwOh8PhGGOG3SmIyO0i8h7wt0Ai8DFQDawA/iQiq0WkaITX/1pEqkXkw2GeFxH5mYjsFJGtIrLk8/1THA6Hw/F5GWmnkAwsV9X2oZ4UkROBeUD5MK+/H/gF8Nthnr8w/Pp5wFLgV+G/HQ6HwxEjhhUFVf3lSC9U1S2jPP+6iBSPcMnlwG/V5oG+LSIZIlKoqvtHel/HBCbUDYFWOgIt1AU6qNNO6oPdBHAjYccDwVCAlkArzcF22jUQ63AcQ1CSXMySzOOj+jMiOVM4GruLz1fVhSKyCLhMVb/3OX/2dGBPn8cV4e8NEgURuQu4C6CoaNiMlSMWhAIQaIVAS+8fDbTSHmijLtRNnQap0xB1KHUeH3XeRNo8vn5vITEKfaqiqhDqJBjsoDvUTRdKOx46PHH95Nn9fxl/ZLRXxl4UgP8E/hr4DwBV3SoiDwOfVxSG+swNecuoqvcB9wGUlJS428poEwqGF/j+i33v425b+FtDndThoc6baH/iEnu/7vClgPjB60fER7p4yBIPCySOLI+PLE8CWXGJZHqT8Hm8sf4XT1oCoQDVnTVUdTVQFWihMthBFQHavQl4vQl4gUKEAo+PfG8S+b5UCuKzyfClI+JkYSoSiSgkqeq7Az4gR2JvWQHM7PN4BrDvCLyvYyhCQQi2QvfAxb7P190tdk3AjpEUaPbEH1zwfWnU+TOoS5xOnS+VLm8CePzg8eMRHxmeOLI8fmZ44snyeMkCsoAMJqCfygRDVWkJtFDZWUtVdxNVgVYqQ13UooTC1/iBPPFyXFwS+d4U8v3p5MfnEO+Nj2XojnFGJL+rNSIyl/BdvIhcwxApnsPgWeDrIvIIdsDc6M4TDhEN9bt7P7joD7Xwtw35FiFvPE3+LOr8WeHFPp06Xyp1ccnUexPp9sb3Lvxe8ZCBLfSzwn/3Xfjd/f7YEAgFqOmso7K7nqruZqqCHVRqgL7/hzMQ8j0+FngTyPelke/PIsuf6e7+HaMSiSh8DUvdzBeRvUApcPNoLxKR3wFnATkiUgH8A+ADUNV7geeBi4CdQBvwxcOIf/KhIVvAh0jZDFrsg+2gQ2TTvH6IS7E/8TmEkotp8GXYnX54wa+LS6LOk0C9x0uwz0vjgExsoZ9L/4U/HdftONa0dLdQ1VlLZXcjVcE2qkJdHNBQ792/DyVPvMz3JpIfl0yBP4P8+BwSvAkxjdsxcYlEFFRVzxWRZMCjqs0iMjuCF9042ptigjP56V3oR0nbdLdAsG3ohd7js0XelwLx2ZBcdHDhj0smGJdCvS+VOm8SdV4/ddD7pwF6FxEwZc4CcoFj6L/wp+EOGGNBMBSkprOOqu56KrubqQq1UxUK0NLnmjSgQHwc7UshPy6VAn8WWfGZeMRJtePIEYkoPAksUdXWPt97AjgpOiFNEFRtAR81bRP+e9iFPjm82GdC0syDj8OLfe/XXj8BoB5b6Guh38LfSP9T+nhskS8EjqP/wp+CW/hjSVugjcrOGqq6GqkMtvbe/ffs2OKAXDzM8ySQ70sh35dOfnw2SXFJsQzbMUUYVhREZD62nqSLyFV9nkoDJufeVNVSMqOlbQKt4YU+NPg9PHEHF3RfOiRO77PwD1jsPX4YkOPt4uDCP/BPE/0X/kRskZ8JnED/hT8Jt/DHmpCGqO2so7KrjqpAC1WBdiq1m+Y+16QC+eJjblwyBb5U8v2Z5Piz8LiKLEeMGGmncAxwCXaGeGmf7zcDfxbNoKJCdzN01gxe3HsX/ZbhF3rxHlzQfWmQOK3/4t53sffED1roB9LB8At/84Brk7FFvpj+i34WJgqO8UF7oD2c+2+gKtBKVaiTag31lul5gVwR5ngTyPcmU+C3u//kuORYhu1wDGKkjuZngGdE5DRVXT+GMUWH+veh8uWDj8UbvpsPL+iJBUOnbeKSwZsw6kI/kHaGXvTrgNYB16Zii/xRDF74XbHg+CKkIeo666nqqqcy0ERVsJ2qUDeNfa5JBgokjlPiksn3pVDgyyInPguvu/t3TAAiOVPYLCJfw1JJvWkjVb0jalFFg4zjIGl6nxz9oS/0fVGsZGq4hX+gYVQ6tsjPp/+in4nVjzvGHx3BDqp6cv+Blt67/+7w8x4gB2GWJ578uGTyfekUxGeT4kuJZdgOx+ciElF4APgIuAD4LlaOuiOaQUUFf6b9OQQUaGH4hb+zz7UCvTX8Cxlcw9/f2MExnlBV6rsaqOyqpaq7mcrw3X9DnxOcJJR88XFSXJLl/n0Z5MZnE+dxbXmOyUUkn+ijVPVaEblcVVeHLS7+GO3AxgrFDnCHW/i7+1zr4WANfxGueWsi0hnspLqzlsquBqqCLVQGO6nWIF3h5wW7+5/h8YcFIJ38+CxS41Jd45djShCJKPSsiw0ishCoxM49JxQtQBWDF/16+nt2eDm40M+h/8LvmrcmDqpKQ3cjVZ11VHY3URVsozLUTX2fu/8EoEC8LI5LJj8uhQJ/Jrnx2fg8bl/nGGeoQls5xOdBXHRLTCIRhftEJBP4NmZNkQJ8J6pRRYEtQM8xc0/zVg5wNP0X/lTcwj/R6Ap2Ud1VS1VXA5XdLVSFOqjSYG96T4AshGkenwmAL5UCfzZpvjR39+8Y32gIGrdDzXpo2wuF50Hu8qj+yFFFQVX/K/zl69jN84TkOKye3zVvTVxUlabupoO5/0ArVaFu6tDe+/94IF+8LPImU+BLId+fQZ4/G7/XHec7JhDBTqh7D2rfga4GczGYfglknhD1Hz2iKIiIF8hU1ZrwYz9wG/AtVT026tEdQTLDfxwTg+5QNwc6a6nsqj9o+axBOvpckwXke/wc702iwJdGfnyWs3x2TGy6m6DmHajbaMKQPAsKV0HaMZ+rWvJQGKmj+QZshkKriHwK3I1VIm0gAkM8hyMSVJXmQHM4999IVaCt1/K55+7fD+TjYWHY8rnAn0FefLazfHZMHtr3W4qoITzSPn0B5JxmZfTYQe5aYAU2iSyajLRT+DZwkqruFJElwHrgBlV9KsoxOaYI9Z11PNjwAbV9Dn8zEAo8cRzXM/DFn0Wms3x2TEZUoXkn1LwFLaXmbpx9CuScCv4M4KAY7MBSowuIrSh0qepOi13fE5FSJwiOI0V9Rw33N35AF8KF8TnO8tkxdQgFoGGr7Qw6Dph1TuF5kHWSNdXSXwwSsBkESxkba5uRRCFPRL7V53FK38eq+m/RC8sxmalvr+L+xm10eeK4LW0BBQk5sQ7J4Yg+gTao3QC175rPWmIBzLwK0o+DsAVKJbAG6xbuEYNTGVsH0pFE4T+xCs3hHjsch0xD2z5WN26nyxvPrenHURCfFeuQHI7o0llru4L69yHUDanzIHcZJBf3Hh6PBzHoYSRDvHvGMhDH5KehdQ/3N31MpzeBWzOOp9CfHuuQHI7o0NNsdmA9NH8MeKycNOc0SMjtvWw/liYaD2LQgzNucYwJDS1l3N/8CR1xKdyWcTyFPrfpdExCNASNO+zwuG0vxCVB3hmQdbI5MocZj2LQgxMFR9RpbPqM1a076fClc2v68RT63AwBxyQj2An1m6Hm7T7NZhdD5ok2YTHMQDE4GztAHg9i0IMTBUdUaWz8mPtby2j3Z3BrxiKmRdm3xeEYU3qbzTZBsMNmpxeugrSjoc/s7IkgBj2MKgoikg/8EzBNVS8UkQXAaar631GPzjGhaWrYxv1te2iLz+bWjOOZ5spNHZOF9kpLETV8CGi42WxZb7NZDxNJDHqIZKdwP/Ab4P+EH38CPAo4UXAMS1Pd+9zfsY+2hDxuTV/IdNd97JjoDNtstnTQrJb9WDXRx0wcMeghElHIUdXHRORvAVQ1ICLBKMflmMA01W3i/vZqWhML+ELGQqY7K2rHRCaCZrMeJrIY9BCJKLSKSDY2jwYRORX6jaR1OHppqnmX1Z01tCZN5wvpxzHDTSZzTFQCbWZMV/sudLeEm82uhPSFvc1mPUwGMeghkt/Y/4XNUZgrIm8CucA1UY3KMSFpPvAmq7saaU6eyRfSFjDDDap3TEQ6a62KqH5LuNnsKJi5DJJnD3Iq3YedGfSIwTnAKUxMMeghknkKm0TkTOAYbAzBx6raPcrLHFMJVZqr13J/oI3m5Fl8If1YZoobVeSYQKhC2x448FafZrNF4WazvEGXT0Yx6CGS6qP3sYPlR1X1s+iH5JhQqNJc+QqrQ100p8zmlrT5zHSOpo6JwqBms0TIPd0OkPs0m/XQVwwSmVxi0EMk6aPLgOuBx0QkhAnEY6paHtXIHOMfVVr2v8hqVZpS53JL6tEUOUFwTAQGNZtlWbNZxglWVTSAocRgKWZnPdmIJH20G/gh8EMRmQf8PfADbMa9Y6qiIVr2/oHV4qEp7RhuTpnrBMEx/ulusoPj2o0jNpv1sA87QP6EyS8GPURUGiIixcB12I4hCPx/0QvJMe4JBWnZ+yyrPX4a0o/lluTZzHKC4BjPtFeGJ5t9wMFms9MgacaQl09FMeghkjOFdwAf8DhwraruinpUjvFLKEDrnqdYHZdEQ8ZCbk6exaxYx+RwDIUqtHxmh8ctu8LNZieHJ5sNPbF9oBisxM4MpoIY9BDJTuE2Vf0o6pE4xj+hblrLn2S1L5WGjEXcnDyT4ljH5HAMJBSwHUHNeuioBl8qFJxrzWbDeG/txc4MprIY9DCsKIjILar6IHCRiFw08Hk3eW2KEeyitfxxVvszqM88kZuSpjtBcIwvBjWb5Q/bbNaDE4PBjLRT6PE3Hsr4Xof43iBEZBXwU+xQ+r9U9V8GPF8ErAYywtf8jao+H8l7O8aQYAdtux/lt/E51GWdxE2JBcyOdUwORw+ddeHJZqM3m/XgxGB4Rpq89h/hL19W1Tf7Piciy0d7YxHxAr8EzgMqgA0i8qyqbu9z2bex8tZfhd1Xnwd3AzquCLTTtvt3rE7Ipzb7ZG5KyGNOrGNyOHqazWrWQ9NHjNZs1sNe7MzgU5wYDEckZwo/B5ZE8L2BnALs7DmYFpFHgMuBvqKgQFr463TsnMcxXgi00lb2ML9NnEZt9incmJDrBMERW3qbzdZDW8WozWY9ODGInJHOFE4DlgG5IvKtPk+lEVmPwnRgT5/HFVhVV1/uBl4Skb/A0lXnRvC+jrGgu9kEIamImpxTuDE+m7mxjskxdQl29Wk2q7dms2kX2WSzIZrNenBicOiMtFPwAynha/qeKzQRmSHeUMm8gWcRNwL3q+qPwyL0gIgsVNVQvzcSuQu4C6CoqCiCH+34XHQ10l72IA+kzKEm+1RuiM90guCIDd3NUPtOn2azmVB4PqQdM2SzWQ99xSAJu9s8GScGkTDSmcJaYK2I3B/uaj5UKoCZfR7PYHB66E5gVfjnrReRBCAHqB4Qy33AfQAlJSURHXI7DpPOOtrLHuK3qfOozjmNG/3pHBXrmBxTj/aqPpPNQpB2LOQuG7bZrIcK7ADZicHhE8mZQpuI/CtwHH18n1T1nFFetwGYJyKzMeG+AbhpwDXl2I7ufhE5Nvz+ByKM3XGk6aihvewhHkibT3XuMm7wpTpBcIwdPc1mNeuh+bNws1nJiM1mPTgxOHJEIgoPYSZ4lwBfAW4jgoU7PKHt68AfsTOIX6vqNhH5LrBRVZ/FZjX8p4h8E0st3a6qbicQC9qr6Ch7iAfSj6Mqdzk3+FKYF+uYHFODw2g268GJwZFHRluDRWSTqp4kIltVdVH4e2tV9cwxiXAAJSUlunHjxlj86MlL2z46dj/MAxknUJmzjOt9yRwd65gck59Ae7jZ7J2DzWY5p0H68cM2m/VQgZ0Z7MTEYBl2gDz8kbMjvJaXjHZdJDuFnoE6+0XkYuxcYOTEnmPi0LqHjt2/44HMJVTmLue6uEQnCI7o0lkXnmy2+WCz2YzTIGXOsM1mPQwUg3NxYnCkiUQUvici6Viq5+dYSeo3oxqVY2xoKaOj/FEezCqhMmc518UlcEysY3JMXlr32OFxT7NZxvG2M0jMH/WlA8XgPCxN5MTgyBPJPIXnwl82YvOoHZOB5p10lj/Og9lL2Ze7nOu88U4QHEceDZkIHHirT7PZinCz2VAOOv1xYjD2RGKd/bMhvt2IHRY/c+RDckSdpo/pLH+SB3OXsS9nGdd5/cyPdUyOycXAZjN/ZkTNZj3swQ6QnRiMPZGkjxKA+dg8BYCrgW3AnSJytqp+I1rBOaJAwzY6K57mwdzT2ZtzGtd6fU4QHEeOnmazuk12kBxhs1kPTgxiTySicBRwjqoGAETkV8BL2P+vD6IYm+NIU/8+nXv/h4fyzmRvzmlc44nj2FjH5JgctFf1mWwWbjbLOc1EIQKcGIwfIhGF6ZgvUWP4cTIwTVWDItIZtcgcR5a6TXTue4GH8s+mInsp13jiWBDrmBwTG1WbaFbzljWbeXzWW5BzqnkTRcAe7MzgM5wYjBciEYUfAltEZA3mZ3QG8E8ikgy8HMXYHEeKmnfo2v8SDxecS0X2KVzt8TpBcBw+oQA0fmg7g/aqcLPZSsgqGbXZrIe+YpAMnA+U4MRgPBBJ9dF/i8jzWDmwAH+nqj0eRn8dzeAcR4DqN+iqeo2HCi9gT9ZJXO3xclysY3JMTHqbzd61s4OEPJh5RXiyWST3l04MJgKRVB8J5k80R1W/KyJFInKKqr4b/fAch40qVK+lq3odD09bRXnmSVzt8ThBcBw6g5rN5sKMKyJqNuvBicHEIRJ5/3cgBJwDfBdoBp7EUn+O8YgqVL5MV83bPDztYnZnnchV4mFhrONyTCxaeyab7eBQm816cGIw8YhEFJaq6hIR2QygqvUi4v6fjldUYd8LdNdu5HczL2N3xglcJcLxsY7LMXFoLYfKP5koeBMOqdmsh3KsmsiJwcQjIu+j8LxlBRCRXGzn4BhvaAj2Pkd33RYennklZRkLudIJgiNSNATV66B6DfjSYNqFkLk4omazHpwYTHwiEYWfAU8BeSLyfWzq2rejGpXj0NEQVDxNd/2H/K7oGsrSj+VKERbFOi7HxKC7Gfb8HlpKIXMRTLsYvJEbUJdjaaJdmBhcgImBLxqxOqJKJNVHD4nIJuywWYArVHVH1CNzRE4oCHuepLvxI35XfB2lafO5ApwgOCKjeSfseQpCXVZNlHFCxAfITgwmHxHVkanqR8BHUY7FcTiEAlD+GN1NO3mk+EZK0+ZxOXBCrONyjH9CQah6FQ68aYfHM6+BhNyIXurEYPIyrCiISDPhcwRsh9DzdRzgV9XICpMd0SPUDbsfIdBcyqOzb2FX6hwuB06MdVyO8U9XA5Q/Yc6l2SVQeIF1JI9CXzFIwYnBZGTYhV1V+5UaiEgq8OfAl7EzBkcsCXZC2cME2vbwyNzb+Cx5FpfhBMERAY07oOIZQKHoWsgYvXvFicHUIZLmtQzgG8CtwMPAyapaG+3AHCMQ7IDSBwm07+fROV9kZ/JMLgcWxzoux/gmFID9f4TaDZA0HYquMUvrEdiNVRM5MZg6jJQ+ysGmrV0P/BpYrKqNw13vGCMCbVD6AIGOAzw69w4+TZrOZThBcIxCRw3seQLaKyF3GeSvHHEO8kAxWAWchBODqcBIO4XdwAHgN0AbNj+h90lV/bfohuYYRHcLlP6WQFc9j839Ep8mFXApsCTWcTnGN/Xvw74/gMRB8U2QNvwU7lbgOWAHTgymKiOJwr9y8HA58lZGR3ToboJdqwl0N/PYUX/GJwl5XIr9wjocQxLsMjGofx9SimHmVdaUNgw7MEHowOrPT8WJwVRkpIPmu8cwDsdIdDXArtUEg+08ftRdfJKQwyU4QXCMQHulVRd11UL+WZB3xrCTzzqAF4D3gULgNiBvzAJ1jDdcWel4p7PWBEEDPH7UXXwcn8XF2GGfwzEIVbO33v9H8CbC7NtslzAMnwHPAC3AWcDpwPAnDY6pgBOF8UzHAShdTVCVx4+6i4/8GVyEs6d1DEOwAyqehcbtkDrPupPjkoe8tAv4E7AByAVuAKaNXaSOcYwThfFKeyWU/pageHniqLv4yJ/GRdikI4djEG0Vli7qboLC883iehirinLgaaAeWIZ54ruFwNFDJH0K+cA/YXOZLxSRBcBpqvrfUY9uqtK2F0ofIOhN4Im5d7LDl8qFOEFwDIGqzUiufMUOkefeAUkzhrw0ALwGvAVkALcDs8YsUMdEIZIbhPuxstT/E378CfAo4EQhGrSWQ9lDBL3JPHnUneyIS2YVsDTWcTnGH4FW2PM0NH8K6QtgxmU2/2AI9mM2BNXYedT5ODtrx9BEIgo5qvqYiPwtgKoGRCQY5bimJi27oOx3BH0Z/H7uF9kel8QqrDTQ4ehHS6lZXQfbYfrFkFUyZLooCLyBNaIlA7cAR41tpI4JRiSi0Coi2RwcsnMq4DqbjzRNn0L5o4Tis/n97NvZFpfIBThBcAxAQ1C9FqpfB382FN8MiQVDXnoA2x3sw2zULwQSxy5SxwQlElH4FvAsMFdE3sSKFa6JalRTjcYdUP4EoYQ8fj/ndrZ54zkfOC3WcTnGF91NUP4ktO6GzBNh2kVDTkVT4G3gFSxFdB2wYGwjdUxgIhmy856InAkcg1lof6yq3VGPbKrQ8AHseYpQ4jR+P/tWPvT6OR+rCnE4emn6BCqeBg3AzCshc+iJGfVYZdFuYD5wCWZX4XBESiTVR9cCL6rqNhH5NrBERL6nqu9FP7xJTt1m2PssoeRZPDXrZj70+jgPJwiOPoSCUPky1Ky3NNHMayAhZ9BlCrwH/JHweERs0FJk89McjoNEkj76e1V9XERWYM65PwJ+hSuI+XzUboC9fyCUOpenZ93IB544zgWWxzoux/ihqz48CGcvZJ9i/Qeewb+yTVh+dycwB7gcSB/bSB2TiKHNUPrTU2l0MfArVX2GCKvZRGSViHwsIjtF5G+GueY6EdkuIttE5OHIwp7gHFhvgpB2DE/PuomtnjhWAitiHZdj/NCwDT6912xOZl0P0y8aJAgKfAD8O5Yuugj4Ak4QHJ+PSHYKe0XkP4BzgR+ISDwRiImIeIFfAucBFcAGEXlWVbf3uWYe8LfAclWtF5HJ78NV/TpUvkoo4ziemXEVWz1eVmKeMw4Hoe7wIJyN1oRWdA34MwZd1gr8AdgOzMTSRdljG6ljkhKJKFyH2ar/SFUbRKQQ+OsIXncKsFNVdwGIyCPYznZ7n2v+DPilqtYDqGr1oQQ/oVC1IenV6whlnsCzMy7nffFwDk4QHGE6DoQH4VRB3grIO3vIQTgfY+miDuxObRmRbfkdjkgYafJaVp+Ha/p8rxPYGMF7Twf29HlcweBziKPD7/smZs54t6q+GMF7TyxUYf9LULMezTqJZ6dfwhYRzgbOiHVsjtijCvVbYN/z4PHD7FsgdXCLWQfwIrAFKMDm4+aPbaSOKcBIO4VNWNpS+vzdg2JnWiMxVOGDDngcB8zDXHtnAOtEZKGqNvR7I5G7gLsAioqKRvmx4wxVG3RSuxHNWcqzhavYIsJZwJmxjs0Re4Kd4UE4WyFldngQzuCZVrswi+tm7EbiTJzFtSM6jDRkZ/bnfO8KLN3ZwwysuXLgNW+H+x5KReRjTCQ2DIjlPuA+gJKSkoHCMn7RkFkZ129B81bwP/kr2RwWhLNiHJpjHNC+PzwIpw7yz4a80wcNwukCXgbeBXKAO7EtuMMRLSJyzBWRTGyx7nXbUtXXR3nZBmCeiMwG9mKW7TcNuOZp4EbgfhHJwdJJuyILfZwTCkLFU9DwIZp/Fv+TdybviXAmbocw5VGF2nctpRiXDHNuh+TBfqV7MJuKOszuZCVuPKYj+kTSvPYl4K+wO/0t2OdzPWbDPixh47yvY/00XuDX4Qa47wIbVfXZ8HPni8h2rPT1r1W19vP8g8YFoYDdATZ9hBaex3O5y3kP2/afhWsomtIE2qHiGWj6CNKOhhlXQFxS/0uwQ7w3sfLS24HiMQ7TMXUR1ZGzMSLyATbs621VPVFE5gP3qOr1YxHgQEpKSnTjxkjOuWNEqBt2PwbNn6LTLuS5nKVswiqMzsEJwpSmdY9VF3W3QMG5kHPqIGfTSuD3mMX1EqxbNH7sI3VMQkRkk6qOOsk3kvRRh6p2iAgiEq+qH4nIMUcgxslHsAt2/w5ay9AZl/KHrJOcIDgsXXTgDah6DXzpMPdOSOo//DKEWVyvAZKwPOvRYx6owxGZKFSISAaW//+TiNQz+MDYEeyAsoehbQ8640qez1zERqxL2QnCFKa7xc6Wmj+DjIUw/ZJBg3BqsLODvcBCrDM5afA7ORxjQiQuqVeGv7xbRF7D0pyTr5fg8xBoh7IHoX0/WnQtz6cvYAPmY7QSJwhTlpZd4UE4HTDjUshc0i9dpMA7WHWRH7gWOC42kTocvYzUvHaVqv4+/HWmqtar6tqxC22CEGiF0geg4wA663peSDuGDViX6bk4QZiSaAiq1sCBdRCfA7NvhYT+Di4N2Na7DEsTXYazuHaMD0baKXwbO/MCm9exJPrhTDC6m6H0t9DVgBbfxIupc3kXE4TzcIIwJelqhD1P2qztrCVQuKrfIBwFNmNbbcF8X07EfVYc44eRREGymGdYAAAgAElEQVSG+doB9stfuhoCLWjxzfwxpZh3sGlpThCmKE0fhwfhBKHoasg4vt/TzcD/AJ8AszFBGGx153DElpFEIVFEFmNeWwnhr3vXuik9ZKezznYIwQ60+Av8MXkmb2MNHOfjBGHKEQqEB+G8DYmF5mwa39+z9EPM1bQbm5V8Cu5z4hifjCQK+4F/C39d2edrsF3wiM1rk5aOGtshaBCdcxsvJRb2CsIFuF/0KUdnHZQ/bpYVOUuh4Lx+cw/aMDHYhnV/XoHZVTgc45WRvI/OHstAJgTtVbZDAHTO7fwpIY/1mPWrE4QpSMMHsPc58yuadQOkz+/39CeYxXU7VoW2HGdx7Rj/ROR95ADa9kHZAyA+dPat/Ckhh7ewNMAqnCBMKULdsO8FqHsPkotg5tXgPzjvrBM7SN6MWVvfglldOxwTAScKkdC6B8oeAm8COuc2XvZn8hbm/XEhThCmFB3Vli7qrDFX0/yz+zmblmKlpk1YJ/uZuF8yx8TCfV5Ho6UMdj8McSno7Nt4xZ/Om5ggXIQThCmDKtS/B/teBG88FN8CqXN7n+7GmtDewcZi3omdITgcE41DEgURuVtV745SLOOP5s/My8ific6+lVd9qbwBlOAEYUoR7IS9/wMNH5oQzLgSfAdbzSowm4pa7HzpXJzFtWPicqg7hcuAu6MQx/ij6WNzO03IRWd/gVfjklkHnARcjBOEKUPbPksXdTdCwUrIXdFrVRHEDOzeANKA27D+A4djInOoojA11sKGbdaVmliIFt/Ca3GJvYJwCVPlP8IURxVq34H9f7JdwZzb7VA5TBXW7l8FLMaqzxKGfCOHY2JxqKJwUlSiGE/Uv29dqUkzofhm1njjeR3z+HCCMEUItNlnoOkTSJsPMy6HuETALK7fxHYIidjYQOcj75hMHJIoqGooWoGMC+o2Wd15cjHMupE1Xj9rsTvBS3GCMCVo3Q3lT5rR4bQLIfuU3nRRLXZ2UIG5mV6Ms7h2TD5c9VEPNe9Y7XnqPJh1HWs8PtZggnAZThAmPRo6OAjHnwlHfcksK7D2/Xex6qI44Gps7oH7TDgmI04UAKrfMO+a9GNh5tWs9cSxBnOvdIIwBehutrkHLaVmYjf9Eis7xSyun8H6D+Zhn4fU2EXqcESdwxIFEfmiqv7mSAcz5qhC9Vrzvs9YCDOu5HWPl9eAE3CCMCVo3gl7noJQl50dZJ4IIiiwBetMVuyz0M8R0uGYpBzuTuEeYGKLgqrtDg68aQvBjMtYJx5exQThcpxPzaQmFITq12yXmJAHRbdDQi4ALZjF9cdAMfZZyIxVnA7HGDPS5LWtwz2FWbpMXFRh/4t2jpB9Mky7iHUivAIswgnCpKerITwIZw9knQTTVoHH2s22Ac9hHcoXYO63bnfgmEqMtFPIx34v6gd8X4C3ohZRtNGQVRjVvQc5p0Hh+bwRFoTjMWtjJwiTmMYdUPEMoDb3IGMhYE6mf8DmHkzHPge5MQvS4YgdI4nCc0CKqm4Z+ISIrIlaRNFEQ1Z/Xr8V8s6A/LN5U4SXMUG4EicIk5ZQAPa/BLXvQtI0KLrWqoyATzGL61bgbMzIzn0OHFOVkeYp3DnCczdFJ5woEgpayqBxOxScA3ln8BbwJ6y80AnCJKazNjwIp9J2hwXngsdLJ/BH4D0gD7gJKIxpoA5H7BnpTCFFVVtGenEk14wbqteaIBReALmn8RbwEiYIV+EEYdJSvxX2PQcSB8U3QdrRAJRhFteNwArgLFx9tmMcE2iz5tr0BYNGvR5pRvo9eEZEtmBl2ptUtRVAROZgu+zrgP8EnohqhEeK3GWQkA8Zx7EeE4TjcIIwaQl2wb7noX4LJM+CoqvBl0Y38ArwNpAF3AHMjGmgDscItFeaB1fDB5YC9cTHThRUdaWIXAR8GVguIplAAKvU+wNwm6pWRjW6I4k3ATKO420sZbAAJwiTlvYqSxd11UL+mZB3JoiHvZhNRQ02D+M8wB/TQB2OIdAQNH1k1ZGtu60yLvNEyF7aWzYdTUbcMavq88DzUY9ijHgba0ZagFkVeGMbjuNIo2pb7P0v2k3A7FshZTZBYC1mcZ0CfAGYO+IbORwxINBmVZG170J3kxVCFF5gghA2ZBwLpkwadSMmCMfiBGFSEuyAimft3Cj1KJh5JcQlU4XtDiox25JVOItrxzhjYIooZQ5Mv9h82GTscxlTRhRmYDYFl+AEYdLRtjc8CKcJCs+DnGWERHgLeA0TgRuA+bGN0uE4yLApolOswz6GTBlRKMA6lR2TCFWoWW92Jb40mHsHJM2gDtsd7MF2hpcAyTEN1OEI05MiqtsAXY3gz4DC8yFz8bApoi5gO1Y6fSbRT32OKAoi4gG2qurCKMfhcBwagVbY8zQ0f2ruttMvQ+MS2YD1nnixQoLjcTYVjnFAe6WdFTRsDaeIZkPhhVYiPUSKSIG9wGasy74TyMYqfaLNaAfNIRF5X0SKVLV8DOJxOEanpcwaEYPtlnvNKqFRhGeAXcBRmKtpWkyDdEx5elJEte/aZzaCFFEb8D4mBtWADyudX4KVTo/FDU4k6aNCYJuIvIs5AQCgqpeN9kIRWQX8FLtx+y9V/ZdhrrsGeBw4WVU3RhK4YwqiIah+3RoR/VlQfDOaWMD7wAvY3dUl2MxYtztwxIxDTBGFsJuZ97B6/yB2Bnop1lwbP4ahQ2SicM/hvLGIeIFfYuXgFcAGEXlWVbcPuC4V+EvgncP5OY4pQndTeBBOGWSeANMupsXr5zngI6AIM7HLimmQjilNe1W4iiiyFFE9NrNjC9ZZnwScghXExPKoeVRRUNW1IjILmKeqL4tIEpEV8JwC7FTVXQAi8gh21rt9wHX/CPwQ+N+HFLlj6tD0KVQ8BRqwUtPME9iOOTZ2AedjFteuEdEx5mgImj42MehJEWWcYCmixMETBgLADiw9tAvb0c7FPsPHMD4qf0aNQUT+DLgLuwmbizkL3wusHOWl07ECkB4qgKUD3nsxMFNVnxMRJwqO/oSCUPUKHHjLfsFmXkt7Qg4vAFuxvOaVxPauyjFFCbRD/XtQu8Hmc/gzrBw6c8mQKaJKLD30AWbTnoF5BZ0IpI9l3BEQiTB9DbvrfwdAVT8VkUh+D4dK62rvk1bZ9H+B20d9I5G7MGGiqKgogh/tmPB01UP5E9aDkH0yFJ7PTo+PZ7CDrbMwi2vXc+IYU9qr+lQRdYdTRKuGTBF1YCKwGdiHLbbzsUPj2Yzfc69IRKFTVbtE7J8gInH0WdxHoIL+XmMzsP82PaRi5yhrwu9dADwrIpcNPGxW1fuA+wBKSkoi+dmOiUzjdutOBph1HV3pC3gJ60rPBW4EpsUuOsdUQ0PQ9Ek4RVQaThEtGjJFpMBubFewHUsXFQAXYlMdx86s4vCJRBTWisjfAYkich7w59gI29HYAMwTkdlYye0NmGU9AKraCOT0PA4P7vnfrvpoChPqhv1/hNqNkDQDiq5mtz+Tp4EGYBlwDuMj7+qYAgxKEaWHU0SLIS6p36VNHCwlrcO66BeH/xQyfncFQxHJ79ffAHdiO6EvYwZ5/zXai1Q1ICJfx0xJvcCvVXWbiHwX2Kiqzx5+2I5JR0cN7Hnctue5ywnkn8OrHi/rsfzr7cCs2EbomCp0VJv9RG+KqNiM6dKO6ZciCgKfYELwKbZLKMZSm8diPQYTEVEdORsjIucAb6tq29iENDIlJSW6caPbTEwaVKHhfdj7B9uWz7iSfWnzeAo4AJRglRnO4toRVQaliOLCKaKlg1JENVh66H3sfCsVOzBezPguiRaRTapaMtp1kewUbgfuFZFaYF34zxuqWv/5QnRMebqbbW5ywweQUkxw5tW87ktlHeZVdAvWnexwRI1AO9RvtsPjEVJEXcA2TAz2YOXPx2BCcBSTqxw6kj6FWwFEZBpwDdaQNi2S1zocQxLsgANvQs3bdoeWfzbVeafzlHjYjx3IXcjEOJRzTFA6qk0I6t8fNkWkWLVMj/9QF3YIej72GU2JTeRRJ5I+hVuw6r/jsZ3TL7DdgsNxaIQCULfRrCoCbZBxPKGCc3jbn8mrWIroOmwIksNxxIkwRdTKwUPjA9jnssd/aAYT69D4cIjkbv8nwGdYw9prqloW1Ygckw9VSxFVvWpb9JQ5UHgeexILeQnbjs/HfIsm692XI4YMlSIqOBeylvSmiELYItfjPxTCBOAyTBDG2n8olkSSPsoRkeOAM4Dvi8g84GNV/ULUo3NMfJo/g8o/mXVwYiHMvpSq1Lm8iv3yJWNdyYuY/HdgjjFmYIooeZYZ06XN700R1WM7gi1YWWkSZruwBOuJmYpEkj5Kw/zGZmEVV+mYkDocw9O+H/b/CVp2mQXAzKuozzie10T4ALvzWon9ArrKIscRQ0M2Y6PmHfvs9aaIToHEAgC6MRPF94BS7GbkKGxU6zG4LvlI0kdv9PnzC1WtiG5IjglNVz1UvmrporgkmLaK5qwSXvfEsQmr0lgGrMAdJDuOIMEOqOtJEdXbJL4BKaL92K5gK2ZBkYk1Q56Im73Rl0jSR4ug1+LaWUw4hibQCtXrrPtTPJB3Ou25y3nTm8A7WKPPEmycYGpsI3VMJjoO2MFxvxTReb0ponYO+g/txxa8Y7HPYjEuZTkUkaSPFgIPYH0ZIiIHgNtU9cNoB+eYAAS7oPZtKzENdUHmYrryz+IdXxpvYmMEF2KOkOO5sccxgRgyRXR8uIqoAAXKsPTQDsx/qBC4CCuhdDvUkYkkfXQf8C1VfQ1ARM4Kf29ZFONyjHc0ZBUdVWusCS1tPsGClWxKyOV1oAU4Gjs3GOwq73AcBkOmiFZC1kkQl0QTdmC8GTtATsB2BD3+Q47IiEQUknsEAUBV14hIchRjcoxnVG2oSOXL0FkDyTMJFV3LB8lFrMF+GWdh/QbO5NxxROg4ELarft92psmz7LwgbT5Bj5ePMSHYieW3Z2NnBfOZuP5DsSQSUdglIn+PpZDA3AdKoxeSY9zSWm7lpa17ID4HnXUDH6cdw6siVGN3Y7dgk5hcrtbxuVC1FFHtO1bW3JsiOgUSCzmApYe2Ys1maViH7WLsANlx+EQiCndgc5p/H378OvDFqEXkGH90HIDKV6DpI/ClwoxLKc1czCvioQLIBq7FOpGdGDg+Fz0poroN0FnXJ0W0hM645F7/oQqsdLTHf2guk8t/KJaMKAoikotlA76jqg1jE5Jj3NDdZGcG9ZvB44eClezLXsorXj+fYXdnl2Elfe4X0vG56KgJD73vSREVQf5KNG0+ezxeNmOGdF1YU9kFWMOjy2MfeYYVBRH5EvBPWPf3bBG5y81AmCIMNKzLXkpN3hm8GpfEdqzr8wLgZJwrouNzMDBFJF5LEeUspSWxsNd/qAZrcFyIHRxPx+1Io8lIv9PfAI5T1QMiMgd4CHCiMJkZwrCuseAc1vgz2YId2p0FnMbU8oJxHGGCHVC/xQ6Pe1NE5xDKOomdccm8hw2vCWHzfC/H/Idc5/vYMJIodKnqAQBV3SUibh2YrAw0rEudS2vBeaxLLGBD+JJTsS5kt113HDYdNeEqoi39UkR1afPZ7PGyBWjGPmOnYmcFU9V/KJaMJAozRORnwz1W1b+MXliOMWOAYV3H7MtYnzqH9ZhHzInY7iA9pkE6Jiyq0LwznCLa2Zsi6s4+hR1J03gPazQTYB7WYHY0zn8olowkCn894PGmaAbiGGMGGNZ1F13NhvSFrBOhHaskOgcbKuJwHDKDUkSpaME57M8qYXNcEh9g/kNZWIPjCTj/ofHCsKKgqqvHMhDHGDHAsC44bRVbskpY64mjCSvtW4mN1nM4DplBKaKZtOevZGs4RVSJLToLsEPjWbhD4/GGKx6ZKgRa7QC5diOIB807g225y3nNG08tNlDkKswkzOE4JIZIEWnGQkpzlvFeYj4fYf5D04CLMf+hhJgG7BgJJwqTnQGGdZq5hJ35Z/OKL4VKIA+4Ecvjujs2xyER7OyTIqoFXyqNheexJXMJm+MSacDM507CDo0LYhutI0KcKExWhjCsKy84j1cSstmNWQFchdV+u8YzxyHRWRueaLYZgl0Ekov4uHAVm1Pm8pnHgwJzsDTksbhFZqIxUvPazxlhfoKrPhqnqJodReUrvYZ1lbOu59WkGXyCzUC+GMvnugoPR8T0pojetYYz8VKdtYT3sk9la0I2bdhB8RlYxZrzH5q4jCTiG8csCseRYYBhXV3xTbyWOo8PREgAzgVOwTUBOQ6BASmiTn8GH06/lPfSj2NvXEKv/9ASbHfgdp0TH1d9NBnoOGBW1k0fgy+V5hlXsDZzEe+JBy/WdLYcN1zEcQj0poi2oMFOytPms3nGZWxLnEm3x0Mezn9oshLJ5LVc4P/Hqsh6iwZU9ZwoxuWIhAGGdW0F5/FmzlLe8cShQAlmJ+zGXzpGJdQNbXusb6WlFNr20uJNYEvOcjZnLqLWn0E8JgKLcf5Dk5lIzoAeAh7FUtFfAW4DDkQzKMcoDDCs68o+lbfzzuTNuAS6sJK/s3F5XccIhILQvg9aS00IWveABmnyJlCaNp8duSv4JGUOIW88RdjNxQJc6nEqEIkoZKvqf4vIX6nqWmCtiKyNdmCOIQgFoHYDHHgdAu0EMk9gU8FKXvel0Yrlds/Bjb90DIEqdFTZLqBlF7TthmAXzR4/ZalHUTr9EsqSplPnzwJPHCmY8eFiXFf7VCMSUegO/71fRC4G9mG9To6xYoBhXSh1LlsLV7EmIZcGrOHsBsxR0uEA7DPTVWci0FpqfwfaaPH4KUueRVnBeZQlzaAmPgc8PhKw7uJTsM9TPi49NFWJRBS+JyLpwP8Cfo5Vnn0zqlE5DtLHsE4TC/lo7tW8mjyTA1iH6KVY1Yf7BXbQ3RTeCYSFoKuRNvFRljSdstwVlCYXcSA+F7zxxGMisAQTgQJc5ZDDGFUUVPW58JeNWKraMRa07bOKorBh3a7iG3gl9Rj2ipADXIc1BjkxmMIE2qG17ODhcGcN7RLH7sQCSjNPoiy5iKqEPPAm4hehCDOem43N03Yi4BiKSKqPVgN/1TOOU0QygR+r6h3RDm5KMsCwbu+MS3kl40R2ebykYwNHTsD9Qk9Jgl3QVn7wXKCjkg687E7IoSx9AaXJs6hKyEPjUvCJMBPrKi7GdpWuWdERCZGkjxb1nc+sqvUisjiKMU1NBhjWHchfyau5S9nh8ZMErMJKTJ1lwBQiFIT2ij6Hw3vpBMr9WZSlHUVp/tnsT8hDfWnEiYeZ2OyL2ZgIuM+K43CI5HPjEZFMVa0HEJGsCF+HiKwCfordpPyXqv7LgOe/BXwJM1E8ANyhqrsPIf6JTz/Dum4ask5mTf5ZvB+XiB/L152KG385JdCQDTvqORhu3U2XhtjjT6c0ZS5lOaeyL6GAkD8dr3iZgdlKzMYqP5wIOI4EkXyOfgy8JSJPhB9fC3x/tBeJiBf4JXAeUAFsEJFnVXV7n8s2AyWq2iYiXwV+CFx/KP+ACcsAw7qW9IWsKzyPjf50BCsHXAEkxTZKRzRRtc7hll0mBK1ldAc62eNPpyx5FqXTL2VvYiEhfwYej4/p2GeiGKs088U0eMdkJZKD5t+KyEasBF6AqwYs7MNxCrBTVXcBiMgjWEq897Wq+lqf698GbjmE2CcmAwzrOlJm81bxF3g7MY8AZiZ2Jm785aSlq/Fgw1hLKYHuVip8aZQmz6Cs4HwqEgsJ+jPxeOOZBizDRKCII9c41t3dTUVFBR0dHUfoHR3jiYSEBGbMmIHPd3i3DSO5pKapalM4XVQJPNznuSxVrRvlvacDe/o8rgCWjnD9ncALo4c8geljWNcdn8e7c77IG8lFtIuwEEsVZcc6RseRJdAKLWW9QhDsrKfCl0ZZYgFleWewJ3EagfgsxJNAoQhLsXRQEdFLGVZUVJCamkpxcTEirn5tMqGq1NbWUlFRwezZsw/rPUbaKTwMXILNZu5roS3hx3NGee+hPm1DWnGLyC3YOeqZwzx/F3AXQFFR0Sg/dhzSx7Au6Etj86wbWJt2NM3i4SisQqQw1jE6jgzBTmjdHT4T2EWwvZp9vlTKEnIpzTqZPUnT6fZnIXHJ5ItwMrYTmMXYTSPr6OhwgjBJERGys7M5cODwnYhGckm9JPz34cmN7Qz6NtnOwLqh+yEi5wL/BzhTVTuHieU+4D6AkpKSYWc8jDv6GNapx8+H0y/ltcwTqfN4mQlcgy0GjglMKBA2krOGsVDbXvZ7kylNyKEsfSHlhTPpis+CuFTyRViC7QRmEVvXWicIk5fP+/82kj6FV1R15WjfG4INwDwRmQ3sxZwYbhrwPouB/wBWqWr1IUU+nuljWKca4tO8s3k151QqvX7ysf8I83CNZxMSDUH7/t4zgVBbOZWeRMriMylLnsPuvLPojM8CXzq54ultFpuFs5jui9fr5fjjjycQCHDssceyevVqkpJcWcV4YKQzhQSs+CUn3LDWs4alYWXQI6KqARH5OvBHrCT116q6TUS+C2xU1WeBf8WGgT0eVrdyVb3s8/yDYsoAw7rdOUt5Je8MyuOSyQKuxsZfOjGYQKhC54HeXgFtLaNKfJT5MyhNnsXunKV0xOeAL51sTxzHY+mgYuyD7RiaxMREtmzZAsDNN9/Mvffey7e+9a3P9Z7BYBCvd3K06AUCAeLiYlNkPNJP/TLwDUwANnFwLWvCSk1HRVWfB54f8L3v9Pn63EMJdtwywLBuf8ZCXi1Yyaf+TFKxg5nFuI7SCUNX/UERaCnlgEKpP4Oy5BmUZV5Oe3wO+DPI8vhZgO0EinFzKw6X008/na1btwLw4IMP8rOf/Yyuri6WLl3Kv//7v+P1evnqV7/Khg0baG9v55prruGee+4BoLi4mDvuuIOXXnqJr3/961RXV3PvvfcSFxfHggULeOSRR6irq+OOO+5g165dJCUlcd9997Fo0SLuvvtuysvL2bVrF+Xl5XzjG9/gL/9y8JThlJQUvvzlL/Paa6+RmZnJI488Qm5uLlu2bOErX/kKbW1tzJ07l1//+td0d3dz4YUXsmnTJt5//31OPPFEdu/eTVFREXPnzuWDDz6gtbWVr3zlK5SXlwPwk5/8hOXLl3P33Xezb98+ysrKyMnJ4eGHHx4Uy1gw0pnCT4GfishfqOrPxzCmiUUfw7ra5GJem3UjHybmk4g1aJyCqycf93S39DaMacsuaoJdthNILGT3tFW0xueAP5MMbwLzObgTmBRlw/tehI7KI/ueCQUwbVVElwYCAV544QVWrVrFjh07ePTRR3nzzTfx+Xz8+Z//OQ899BC33nor3//+98nKyiIYDLJy5Uq2bt3KokWL7MclJPDGG28AMG3aNEpLS4mPj6ehwYwY/uEf/oHFixfz9NNP8+qrr3Lrrbf27lI++ugjXnvtNZqbmznmmGP46le/OqiUs7W1lSVLlvDjH/+Y7373u9xzzz384he/4NZbb+XnP/85Z555Jt/5zne45557+MlPfkJHRwdNTU2sW7eOkpIS1q1bx4oVK8jLyyMpKYkvfelLfPOb32TFihWUl5dzwQUXsGPHDgA2bdrEG2+8QWJi7E6cItmfVIpIqqo2i8i3MWPF76nqe1GObXzTx7CuKSGftXNuZ3PyLLwinIHVl49VNYnjEAl29JaJakspdd3NJgIJeZTln0VLfC74M0jzJnGUCMXYbiAjtlFPKtrb2znxxBMB2ynceeed3HfffWzatImTTz6595q8vDwAHnvsMe677z4CgQD79+9n+/btvaJw/fUH+10XLVrEzTffzBVXXMEVV1wBwBtvvMGTTz4JwDnnnENtbS2NjY0AXHzxxcTHxxMfH09eXh5VVVXMmNF/MoDH4+n9GbfccgtXXXUVjY2NNDQ0cOaZVjB52223ce211wKwbNky3nzzTV5//XX+7u/+jhdffBFV5fTTTwfg5ZdfZvv2g61eTU1NNDc3A3DZZZfFVBAgMlH4e1V9XERWYGNZfwT8ipF7DiYvfQzr2nxpvDHret5NPRr1eDkZm1DlcsnjjD6jJrWllIbOOksHxWdTln0KTQl54M8kNS6Z2eLpTQf1PUibtER4R3+k6Xum0IOqctttt/HP//zP/b5fWlrKj370IzZs2EBmZia33357v8a75OSDR/h/+MMfeP3113n22Wf5x3/8R7Zt24bq4ILFngqd+PiD3SBer5dAIDBq7KNV95x++umsW7eO3bt3c/nll/ODH/wAEeGSSy4BIBQKsX79+iEX/77/llgRidlmMPz3xcCvVPUZpuJUvkAr7HsBPv4Fnc2fsnb6Zfz0mL9gffqxLPR4+QvgQpwgjAtCQRsvWbUWdt1Pw47/y5Z9L/FUZy0/yVjET2ddx7OzruezGVdSlF3CJclFfN2XyrfEw9XYVjiLKSAI44yVK1fyxBNPUF1thYh1dXXs3r2bpqYmkpOTSU9Pp6qqihdeGLrHNRQKsWfPHs4++2x++MMf0tDQQEtLC2eccQYPPfQQAGvWrCEnJ4e0tLSI4wqFQjzxhLn8PPzww6xYsYL09HQyMzNZt24dAA888EDvruGMM87gwQcfZN68eXg8HrKysnj++edZvnw5AOeffz6/+MUvet9/oDjGmkh2CntF5D+Ac4EfiEg8U8m5uY9hXSAUYGP+WazLPpnWcH75HCAv1jFOdQaMmmzq2E+pN5kyfyZlqfOozz0d/Bkk+dIp9sSxHEsH5eAW/vHEggUL+N73vsf5559PKBTC5/Pxy1/+klNPPZXFixdz3HHHMWfOnN7FdSDBYJBbbrmFxsZGVJVvfvObZGRkcPfdd/PFL36RRYsWkZSUxOrVqw8pruTkZLZt28ZJJ51Eeno6jz76KACrV6/uPWieM2cOv/nNbwA7/AYTB4AVK1ZQUVFBZqZNTf/Zz37G1772NRYtWkQgEOCMM87g3nvvPZz/ZFFBhtpa9btApMe5+c/pIbEAABE1SURBVANV/VRECoHjVfWlsQhwICUlJbpx48bo/yANQd17UL2GUHcr7+csZU3uChp9KczGupDdTNIYMWDUZHPbXso8fsp8GZQmTacusRD8mST60pnl9femg/JwIgCwY8cOjj322FiHMWFISUmhpaUl1mEcEkP9PxaRTapaMtprR/U+ws5L14S/lwV0AmOwKseIPoZ12lnDjoxFvJp/NjXxmUzHHP1G8/dwRIE+oyZb2/ZQhodSfyZliQXUpB8P/gzi/RkUexN6rSPymUpbWofjyHAo3kd9b7Ii8T6aeIQN67R1D7tS5/JK0TXsS8gnV4Trgfm4O80xo8+oybbWPZRpkDJ/BmUJuVTnnwP+TPz+DGZ5E1kswmzcnGFHdJhou4TPSzS9jyYOfQzrKhILeWXO7ZQmzyJdhCuARbjFJur0jprcRXtrObtDHZTFZVCakENV7jLwZ+LzZ1AUl8KicJloIa4h0OE40kTT+2j808ewrtqXwauzbuCjtHkki5cLgZNw06yOKKEABFqguxkCzdY0Fmiio20v5d0tlPrS/l975x8c1XXd8c95u9pdIQmDUMAYDOK3k6CGkBAgtrEgNm1aD5lx7co000Abx51JRjX2jD2uf1LXnnGd2IY4jWNCA3ZrExo3dmzcOP5RKNgFK4AB/yAYpygxNQmCGKwFoZ+nf9y7q0VIQiCJ3bc6n5mdfe/te7vn6Omd8+65930vtfFSflf6OTQ2lGhsKGOiJcwLAspxWuyWBAyjf+k37aOcJkOw7iOJsWHUAnad92likRjzcNNfDrwxt72gtckF+3Sgr/eBP0lLc5JkWyPJtibqtY1kEKM+iFMfiZEMYnwcSXDwvKlofCiRgiFcWDCYy4Io43BJwJKyYZxb+lX7KOfIEKxLtrWycUQl24ZOR6KFfBG4GJv+Mo0qtDW2B/nUHb4P9s0tSZKtTdRrE/WKC/aROPVBzC8nSMbLOF40FoIYBHGIxCCIIRKjOFJAcRDnvCDGRb4cNBqTBDGMbDNwtI8+3gMf/pwTzUleL5vJlrJZtBaUMB03+XnPH2UJOarQ2pBRxsko57QkaWw+5u7q25qolyB9Z+8CfYz6yCCShedzIjLIB/m4D/oxgiBGSeCCfWlQwBgRSnBCccUZ70VYH81A57777uOpp54iEokQBAGPPfYYM2cOTJGEXKMnczQ/IiJTgU+RIeejqk/0p2F9TbO28kbROF4bfgkn4sOowE1/WZptw/oKVffUdWaw92UcbUlyouW4C/atTdQH0fZgH4lRH8SojxaRLBxNU0lhOsin7u6j4oN9JM5wiTK+i2A/CBudZZyezZs3s27dOrZv3048HufQoUM0NTX12+9lU4b6XNDX/vWko/luoBKXFP4Tp+bwGhCqpLBp8CfZOPiTTBZhHm74YijQNhfsM+7mU2UcbUlyvLXBl3GaqZeCdBnHBf0Y9QXFJONltEQS7cE+4u7uY0EBJRKjOBLngiCaDvAdg30CC/ZG33HgwAHKysrSukNlZWXpz8rLy6mqqmL9+vWAk5WYOHEidXV1ncpN19TUsGTJEhoaGigsLGTVqlVMmTKF1atX88ILL3DixAmOHTvGXXfdxd13382IESPYsWMHV111FRUVFSxfvpyGhgaeffZZJkyYwPPPP8+9995LU1MTw4YN48knn2TEiBEDSma7J+nlauAzwJuq+tciMgJY2WcWnCNmijARNyF6TtDWetLdfOYdfltzPcdaT5DUJurbWtoDfGYZJ1pCsnAEbUEiXatPlXISQQHFQYySIM6YIHJSgM8M+v01MbwRHl4E+lg4m/NxEghdMX/+fO655x4mT57M5ZdfTlVVVVo3CGDw4MHU1NTwxBNPsGTJEtatW8cNN9zQqdz0RRddxMaNG4lGo7zyyivcdtttaUXUzZs3s2vXLkpLS9mwYQM7d+5k9+7dlJaWMn78eK677jpqampYvnw5jzzyCMuWLeOSSy5hy5YtiAgrV67kgQce4MEHHwQGjsx2T5JCg6q2iUiLiAwGDhLCB9eKOEfTIbY1n9IpmxqV09qSJNl6osNInIw7+0icZMFgkoUXoJH4yWWcIMagoIDiSIwSifOJIOgy2FtnrZHLFBcXs23bNjZt2sT69eupqqri/vvvZ/HixQAsXLgw/X7jjTcCXctNHz16lEWLFrF3715EhObm5vQ+V1xxBaWl7QXiGTNmMHLkSAAmTJjA/PnzAaioqEi3TPbv309VVRUHDhygqamJcePaH9MaKDLbPUkKW0VkCPBD3CikJFDT55bkOq2NnYyxbx92Wd/qh11C+918+s4+QTJeyrGiC0/qmCUSR6SAokjM39nHGCldB3sbo2/0NdkRznYy1ZWVlVRWVlJRUcHjjz+eTgqZ0tSp5a7kpqurq5k7dy7PPPMMtbW1VFZWpj/rKEOdKZMdBEF6PQiCtGR2dXU1N910EwsWLGDDhg0sXbq00+PzWWa7Jx3N3/SLPxCRF4HBqrqrX6w516i6Zxa6GGPf1HyM+rZGV8bRoL1TNhXsowmS8eE0FJW3B3s/IifwJZziIMaQIMaFIp0GexuJYww09uzZQxAETJo0CXDS0WPHjk1/vnbtWm699VbWrl3L7NmzgXa56Ztvvjl9zLRp0zh69CijRo0CYPXq1b22LfP7zlRNFdpltq+99tpOZbYvvfTSU2S277jjDubMmXOSzHZqTomu/O5PzuiJZlWt7bgtNCT3wdHdaEuSxtbjJFsbqW9rpl4ip5ZxooNIFl5AY3FhulM2VcaJBAWU+GBfFhQwrotgbyNxDKNzkskk1dXVHDlyhGg0ysSJE1mxYkX688bGRmbOnElbWxtr1qwBupabvuWWW1i0aBEPPfQQ8+bN67VtS5cu5ZprrmHUqFHMmjWLffv2ndHx+SCz3aV0dsYTzetxo48yn2j+uapmRXv3bKWzt3+0i00nDpIsKKY5kjiljFMgBZSkyzgFp4zASQV7G4ljhJ1cls4uLy9n69atJ41IChO5IrPdL9LZ5NkTzUVDKhidcVffMdjHsGBvGIYxYJ5oniLClGwbYRhGt9TW1mbbhF6RC62E3tJlH6eIzBCR81MJQUS+JiI/E5Hv+sl2DMMwjDyju4EvjwFNACIyB7gf9xTzUWBFN8cZhpHjnG4aXiO89PbcdpcUIqr6B79cBaxQ1f9Q1TuBib36VcMwskYikeDw4cOWGPIQVeXw4cMkEonT79wF3XU0R0QkqqotuHnqr+/hcYZh5DCjR49m//791NXVZdsUox9IJBKnPGl9JnQX3NcA/y0ih4AGYBOAiEzElZAMwwghBQUFJ8k3GEYm3Y0+uk9EXsVNhfuStrc1A6D6XBhnGIZhnFu6LQOp6pZOtr3Xf+YYhmEY2cRkdwzDMIw0Xcpc5CoiUgf85iwPLwMO9aE52cR8yT3yxQ8wX3KV3vgyVlU/cbqdQpcUeoOIbO2J9kcYMF9yj3zxA8yXXOVc+GLlI8MwDCONJQXDMAwjzUBLCvkkz2G+5B754geYL7lKv/syoPoUDMMwjO4ZaC0FwzAMoxvyOimISERE3hSRdX59nIi8ISJ7RWStiMSybWNPEZFaEXlLRHaIyFa/rVREXvb+vCwiQ7Nt5+kQkSEi8rSI/EpEdovI7JD6McWfi9TrYxFZEkZfAETkRhF5R0TeFpE1IpII4/UiIjd4H94RkSV+WyjOiYj8SEQOisjbGds6tV0c3xWR90Vkl4hM7ys78jopADcAuzPW/wl4WFUnAR8BX8+KVWfPXFWdljEk7VbgVe/Pq34911kOvKiqFwGfwZ2f0Pmhqnv8uZgGfA44DjxDCH0RkVHA3wGfV9WpQAS4lpBdLyIyFfgG8AXc/9aVIjKJ8JyT1cCfdNjWle1fBib51/XAo31mharm5QsY7f+I84B1uNk2DwFR//ls4BfZtvMM/KkFyjps2wOM9MsjgT3ZtvM0PgwG9uH7ssLqRyd+zQdeD6svwCjgA6AUJ32zDvjjsF0vwDXAyoz1O4FbwnROgHLg7Yz1Tm3HzXezsLP9evvK55bCMtw/RJtfHwYcUScFDrAfdzGEBQVeEpFtIpKSMR+hqgcA/PvwrFnXM8YDdcAqX9ZbKSJFhM+PjlyLUxWGEPqiqv8HfAf4LXAAp4K8jfBdL28Dc0RkmIgMAv4UuJAQnpMMurI9lchT9Nn5ycukICJXAgdVdVvm5k52DdPQq4tVdTqu2fgtPxte2IgC04FHVfWzwDFytynfI3ydfQHwk2zbcrb4OvVXgHHABUAR7v+sIzl9vajqblzJ62XgRWAn0NLtQeGl3+JZXiYF4GJggYjUAj/GlZCWAUNEJKUMOxr4MDvmnTmq+qF/P4irXX8B+L2IjATw7wezZ2GP2A/sV9U3/PrTuCQRNj8y+TKwXVV/79fD6MvlwD5VrVPVZuCnwBcJ4fWiqv+iqtNVdQ7wB2Av4TwnKbqyfT+uFZSiz85PXiYFVf17VR2tquW4pv1/qepXgfXA1X63RcDPsmTiGSEiRSJSklrG1bDfBp7D+QEh8EdVfwd8ICJT/KYvAe8SMj86sJD20hGE05ffArNEZJCICO3nJXTXi4gM9+9jgKtw5yaM5yRFV7Y/B3zNj0KaBRxNlZl6TbY7Vs5Bx00lsM4vjwdqgPdxzf14tu3roQ/jcU3hncA7wO1++zBcZ/pe/16abVt74Ms0YCuwC3gWGBpGP7wvg4DDwHkZ28Lqyz8Av8LdbPwrEA/j9YKbIfJdf618KUznBJfADgDNuJbA17uyHVc++mfg18BbuJFjfWKHPdFsGIZhpMnL8pFhGIZxdlhSMAzDMNJYUjAMwzDSWFIwDMMw0lhSMAzDMNJYUjDyEhG53Stl7vIqpjP7+fc2iEiP584VkVlegXSHV4td6rcvEJFQP+VthJvo6XcxjHAhIrOBK4HpqtooImVArsk+Pw78haruFJEIMAVAVZ/DPZhkGFnBWgpGPjISOKSqjQCqeki9TIiI3CUiv/Sa+yv8E7ypO/2HRWSjv3OfISI/9Tr29/p9yv08EI/7FsjTXnjtJERkvohsFpHtIvITESnuxMbhuAeVUNVWVX3XH7tYRL7nlzPna2gQkcv80+0/8j68KSJf6Ye/nzGAsaRg5CMvAReKyHsi8n0RuSzjs++p6gx18wYU4loUKZrUaeb8ACcn8C1gKrBYRIb5faYAK1T1j4CPgW9m/rBvldwBXK5OwHArcFMnNj4M7BGRZ0Tkb0Uk0XEHbZ+v4U7/Pf8D3I6TbZkBzAW+7aVPDKNPsKRg5B2qmsRNfHM9Tqp7rYgs9h/P9bX8t3BCiZ/OODRVtnkLeEdVD/jWxv/SLj72gaq+7pf/Dbikw8/PAj4FvC4iO3B6NWM7sfEe4PO4BPaXOFXPU/CTxHwbqFInVjcfuNV/9wYgAYzp9g9iGGeA9SkYeYmqtuKC5gafABaJyI+B7+N0Yj7wnbuZd+iN/r0tYzm1nrpWOurCdFwX4GVVXdgDG38NPCoiPwTqMloj7otcC+DfgW+kyl/++/9cVfec7vsN42ywloKRd4ibP3lSxqZpwG9oTwCHfJ3/6lMOPj1jfEc2OIXU1zp8vgW4WEQmelsGicjkTmz8s1R/Bm5KxVbgSIfdVgGrVHVTxrZfANUZfSGfPQsfDKNLrKVg5CPFwCMiMgQ3ycr7wPWqesTflb+Fm970l2fx3btxrY7HcMqVJ82Nq6p1vlS1RkTifvMdwHsdvuevgIdF5Li38auq2prKEyIyFpe0JovI3/hjrgP+ETc3yC6fGGo5uV/EMHqFqaQaRg8RkXKcDPvULJtiGP2GlY8MwzCMNNZSMAzDMNJYS8EwDMNIY0nBMAzDSGNJwTAMw0hjScEwDMNIY0nBMAzDSGNJwTAMw0jz/45PJno0k/c3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Let's see how Spearman or Pearson regression perform \n",
    "#at regression on a linear dataset with a given variance and sample size\n",
    "variances = [20,40,80,160]\n",
    "sample_sizes = [40,60,80,100]\n",
    "n_trials = 300\n",
    "slope = 1\n",
    "intercept = 3\n",
    "results_columns = [\"Test\",\"Variance\",\"Sample Size\",\"Statistical Power\",\"Monte Carlo Trials\"]\n",
    "results = []\n",
    "\n",
    "for variance in variances:\n",
    "    print(f\"Running simulations for variance {variance}\")\n",
    "    pearson_results = []\n",
    "    spearman_results = []\n",
    "    for sample_size in sample_sizes:\n",
    "        #Test how many of the trials \n",
    "        #reflect a positive result\n",
    "        \n",
    "        pearson_successes = 0\n",
    "        spearman_successes = 0\n",
    "        \n",
    "        \n",
    "        \n",
    "        for trial in range(n_trials):\n",
    "    \n",
    "            xs,ys =\\\n",
    "              simulate_noisy_linear_data(xs = range(sample_size),slope=slope,\\\n",
    "                intercept=intercept,variance=variance)\n",
    "            \n",
    "            pearson_r,pearson_p = pearsonr(xs,ys)\n",
    "            spearman_rho,spearman_p = spearmanr(xs,ys)\n",
    "    \n",
    "            if pearson_p < 0.05:\n",
    "                pearson_successes += 1\n",
    "                \n",
    "            if spearman_p < 0.05:\n",
    "                spearman_successes +=1\n",
    "         \n",
    "        \n",
    "        pearson_label = f\"{sample_size}_samples_variance_{variance}_Pearson\"\n",
    "        spearman_label = f\"{sample_size}_samples_variance_{variance}_Spearman\"\n",
    "        \n",
    "        pearson_power = pearson_successes/n_trials\n",
    "        pearson_results.append(pearson_power)\n",
    "        \n",
    "        spearman_power = spearman_successes/n_trials  \n",
    "        spearman_results.append(spearman_power)\n",
    "        \n",
    "        \n",
    "        results.append([\"Pearson\",variance,sample_size,\\\n",
    "          pearson_power,n_trials])\n",
    "        results.append([\"Spearman\",variance,sample_size,\\\n",
    "          spearman_power,n_trials])\n",
    "        \n",
    "    \n",
    "    plot(sample_sizes,pearson_results,\"-\",color=\"orange\",label=f\"Pearson power\",alpha=0.50)\n",
    "    plot(sample_sizes,spearman_results,\"-\",color=\"cyan\",label=f\"Spearman power\",alpha=0.50)\n",
    "    \n",
    "    #Only call legend on the first iteration of the loop\n",
    "    if variance == variances[0]:\n",
    "        legend()\n",
    "    \n",
    "xlabel(\"Sample Size\")\n",
    "ylabel(\"Statistical Power (1 - False Negative Rate)\")\n",
    "from pandas import DataFrame\n",
    "df = DataFrame(results,columns=results_columns)\n",
    "df\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above results, the statistical power of Pearson regression (it's ability to reject the null hypothesis of no trend in noisy linear data) are plotted in orange for different variances, while the cyan lines represent the same for Spearman regression. As you can see, Pearson regression has *slightly* more power when its assumptions are met. (This is indicated by the orange line being sometimes a little bit higher than the cyan one)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "1. The significance test for Pearson regression is testing the chances of getting an *apparent* slope as large or larger than the one you see under the null hypothesis that the slope was 0. Try adjusting the slope of the line. If you make it much larger, does statistical power increase or decrease?\n",
    "1. Use the simulation to test whether adjusting the intercept of the true trend (i.e. the `intercept` parameter) has an effect on statistical power?  Once you get a result, try to interpret it.\n",
    "1. Finally, try adjusting the variance in the data. Does greater variance increase or decrease statistical power?\n",
    "1. You are studying a correlation between days since last feeding and body fat in pythons. If every day that the python didn't feed resulted in a loss of 1% body fat (a `slope` of -0.01), and the `variance` in body fat measurements was about 0.05, how many pythons do you think you'd have to measure in order to be pretty sure (say, with `statistical power` of at least 0.95) that you'd detect a significant association between these two parameters? Use the simulation to test this out. What about if you used Spearman correlation?\n",
    "1. Modify the code for exponential regression to add an error term to each data point (like we did in the `simulate_noisy_linear_data` function. You might call it `simulate_noisy_exponential_data`.  Use your new function to test how Spearman and Pearson regression perform if the data are actually exponential instead of linear.\n",
    "1. Pearson regression assumes a normal distribution of errors. Try changing the distribution that the errors come from to a unifrom (flat) distribution. You can do this by replacing the call to `norm` in the `simulate_noisy_linear_data` with the function `uniform` in scipy.stats. (You'll of course have to call `from scipy.statis import uniform`) before you call it.). How does this change impact the statistical power of Pearson vs. Spearman regression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Citations and Further Reading\n",
    "\n",
    "[1] Derryberry, Shou, Connover [\"Teaching Rank-Based Tests by Emphasizing Structural\n",
    "Similarities to Corresponding Parametric Tests\"](http://jse.amstat.org/v18n1/derryberry.pdf). *Journal of Statistics Education* Volume 18, Number 1 (2010)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}