{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Simple logistic regression\n", "\n", "This notebook follows John H McDonald's [Handbook of Biological Statistics](http://www.biostathandbook.com/simplelogistic.html) chapter on simple logistic regression.\n", "\n", "This notebook is provided with a CC-BY-SA license." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import io\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import numpy as np\n", "import scipy.stats\n", "import scipy.special\n", "import seaborn as sns\n", "sns.set_style('white')\n", "sns.set_context('notebook')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Spider data from *Suzuki et al. (2006)*:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Grain size (mm)Spiders
00.245False
10.247False
20.285True
30.299True
40.327True
\n", "
" ], "text/plain": [ " Grain size (mm) Spiders\n", "0 0.245 False\n", "1 0.247 False\n", "2 0.285 True\n", "3 0.299 True\n", "4 0.327 True" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = \"\"\"Grain size (mm)\tSpiders\n", "0.245\tabsent\n", "0.247\tabsent\n", "0.285\tpresent\n", "0.299\tpresent\n", "0.327\tpresent\n", "0.347\tpresent\n", "0.356\tabsent\n", "0.36\tpresent\n", "0.363\tabsent\n", "0.364\tpresent\n", "0.398\tabsent\n", "0.4\tpresent\n", "0.409\tabsent\n", "0.421\tpresent\n", "0.432\tabsent\n", "0.473\tpresent\n", "0.509\tpresent\n", "0.529\tpresent\n", "0.561\tabsent\n", "0.569\tabsent\n", "0.594\tpresent\n", "0.638\tpresent\n", "0.656\tpresent\n", "0.816\tpresent\n", "0.853\tpresent\n", "0.938\tpresent\n", "1.036\tpresent\n", "1.045\tpresent\n", "\"\"\"\n", "df = pd.read_table(io.StringIO(data))\n", "df.Spiders = df.Spiders == 'present'\n", "df.head()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAroAAAHxCAYAAACCih8rAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XtYVHXix/HPCMJggiJe1vK2dnEwCxBNab2UoKWWoL+V\nVHKtRM0Sy3LTrA21UtxM87duXtJ8FFpLzZK8rCZuPV1XXTEviKldlDRyUjMUGIX5/dHPWRFRBubG\n8f16Hp/H+Z7vzPnMOaEfT985Y7Lb7XYBAAAABlPL2wEAAAAAd6DoAgAAwJAougAAADAkii4AAAAM\niaILAAAAQ6LoAgAAwJAougAAADAkii4AAAAMiaILAAAAQ6pRRddms+n+++/Xtm3bKpzz0UcfKSEh\nQVFRUYqPj9eWLVs8mBAAAAC+osYUXZvNpqeeekoHDx6scM7+/fuVkpKigQMHKjMzU4mJiRo7dqz2\n79/vwaQAAADwBTWi6B46dEiJiYnKy8u74ry1a9cqJiZGSUlJat68uZKSktSpUydt2LDBQ0kBAADg\nK/y9HaAytm7dqpiYGD355JOKiIiocF7//v117ty5cuMFBQXujAcAAAAfVCOK7uDBgys1r3Xr1mUe\nHzhwQF9++aWGDBnijlgAAADwYTVi6UJVnDhxQikpKYqOjlZsbGylnvPggw/qwQcfdHMyAAAAeEKN\nuKLrLKvVqocfflgmk0lz5syp9POOHTvmxlQAAADwJMNd0c3Pz1dSUpJKSkqUnp6u0NBQb0cCAACA\nFxiq6BYWFio5OVm1a9dWRkaGGjZs6O1IAAAA8JIav3TBarUqODhYgYGBmj9/vvLy8rRs2TKVlpbK\narVKksxms+rWrevlpAAAAPCkGndF12QylXncpUsXx31yN23apKKiIiUmJqpr166OXy+//LI3ogIA\nAMCLTHa73e7tEL7iwt0ZsrKyvJwEAAAA1VXjrugCAAAAlUHRBQAAgCFRdAEAAGBIFF0AAAAYEkUX\nAAAAhkTRBQAAgCFRdAEAAGBIFF0AAAAYEkUXAAAAhkTRBQAAgCFRdAEAAGBIFF0AAAAYEkUXAAAA\nhkTRBQAAgCFRdAEAAGBIFF0AAAAYEkUXAAAAhkTRBQAAgCFRdAEAAGBIFF0AAAAYEkUXAAAAhkTR\nBQAAgCFRdAEAAGBIFF0AAAAYEkUXAAAAhkTRBQAAgCFRdAEAAGBIFF0AAAAYEkUXAAAAhkTRBQAA\ngCFRdAEAAGBIFF0AAAAYEkUXAAAAhkTRBQAAgCFRdAEAAGBIFF0AAAAYEkUXAAAAhkTRBQAAgCFR\ndAEAAGBIFF0AAAAYEkUXAAAAhkTRBQAAgCFRdAEAAGBIFF0AAAAYEkUXAAAAhkTRBQAAgCFRdAEA\nAGBIFF0AAAAYEkUXAAAAhlSjiq7NZtP999+vbdu2VTgnJydHiYmJioyM1MCBA7V3714PJgQAAICv\nqDFF12az6amnntLBgwcrnFNYWKiRI0eqY8eOWr16tSIjIzVq1CgVFRV5MCkAAAB8QY0ouocOHVJi\nYqLy8vKuOG/dunUKCgrSn//8Z7Vu3VrPPfecrrvuOv3zn//0UFIAAAD4Cn9vB6iMrVu3KiYmRk8+\n+aQiIiIqnLdr1y5FR0eXGWvfvr2ys7OVkJDg7pjwETt2fqV5i1fodGGpQoJqafTwRLWPjHB6jiv2\n48z8i7efP3tSJj9/+QUGVzmfKzJX9Tmu5M79u/ocepun8zmzP18/dtcSX/qZqgmM+J6upKL3W1OP\ng9/kyZMnezvE1dx2223q0qWL/P39NXfuXA0YMEA33HBDuXnvvPOOWrRooZiYGMfYzp07deTIEd13\n331X3c+yZcskScOGDXNdeHjUjp1f6Zlp6SoKjdF5c3P9Wut6bVi3Vu1ubqKmv/tdpee4Yj/OzL94\n+8kCu37I/1nX/b5nlfO5IrOrjlV1uHP/rj6H3ubpfM7sz9eP3bXEl36magIjvqcrqej9+pee1l8X\nrq+Rx6FGLF2orKKiIgUEBJQZCwgIkM1m81IieNq8xStkbtZFJpNJkmQymWRu1kXzFq9wao4r9uPM\n/Iu3//TdDrW4rVe18rkic1Wf40ru3L+rz6G3eTqfM/vz9WN3LfGln6mawIjv6Uoqer8z56bX2ONg\nqKIbGBhYrtTabDaZzWYvJYKnnS4sdfwgXmAymXS6sNSpOa7YjzPzL97u5x9Q7XyuyFzV57iSO/fv\n6nPobZ7O58z+fP3YXUt86WeqJjDie7qSit5vSa06NfY4GKroNmnSRMePHy8zZrVa1ahRIy8lgqeF\nBNWS3W4vM2a32xUSVMupOa7YjzPzL95ect5W7XyuyFzV57iSO/fv6nPobZ7O58z+fP3YXUt86Weq\nJjDie7qSit6vX+nZGnscfD+hEyIiIpSdnV1mLDs7W5GRkV5KBE8bPTxRRXmfOn4g7Xa7ivI+1ejh\niU7NccV+nJl/8fbGrdrr8O5N1crnisxVfY4ruXP/rj6H3ubpfM7sz9eP3bXEl36magIjvqcrqej9\njh8ztMYeB5P90oru4ywWi9LT09WxY0dJv12xDQ4OVmBgoAoKCnTPPfeob9++euCBB7R8+XJt3LhR\nH374YaWWL8TGxkqSsrKy3Poe4F47dn6l+W+u1C9nS65414WrzXHFfpyZf/H282dPqpa/v2oFuP6u\nC86+b1ccq+pw5/5dfQ69zdP5nNmfrx+7a4kv/UzVBEZ8T1dS0futqcehxhXd8PBwLVu2zFF0LRaL\n0tLSHLcP2717t1JTU/XNN9+oTZs2mjJliiwWS6Vem6ILAABgHDWu6LoTRRcAAMA4DLVGFwAAALiA\nogsAAABDougCAADAkCi6AAAAMCSKLgAAAAyJogsAAABDougCAADAkCi6AAAAMCSKLgAAAAyJogsA\nAABDougCAADAkCi6AAAAMCSKLgAAAAyJogsAAABDougCAADAkCi6AAAAMCSKLgAAAAyJogsAAABD\nougCAADAkCi6AAAAMCSKLgAAAAyJogsAAABDougCAADAkCi6AAAAMCSKLgAAAAyJogsAAABDougC\nAADAkCi6AAAAMCSKLgAAAAyJogsAAABDougCAADAkCi6AAAAMCSKLgAAAAyJogsAAABDougCAADA\nkCi6AAAAMCSKLgAAAAyJogsAAABDougCAADAkCi6AAAAMCSKLgAAAAyJogsAAABDougCAADAkCi6\nAAAAMCSKLgAAAAyJogsAAABDougCAADAkCi6AAAAMCSKLgAAAAypRhRdm82mSZMmqWPHjuratauW\nLFlS4dwPP/xQffv2VVRUlJKSkpSTk+PBpAAAAPAVNaLozpgxQzk5OUpPT1dqaqrmzp2rTZs2lZt3\n8OBBjR8/XqNGjVJmZqYsFotGjhyp4uJiL6QGAACAN/l80S0sLNSqVav0/PPPy2KxKC4uTsnJycrI\nyCg399NPP9XNN9+sfv36qXnz5nrqqadktVp18OBBLyQHAACAN/l80c3NzVVJSYkiIyMdY9HR0dq1\na1e5ufXr19fBgwe1Y8cO2e12vfvuuwoODlaLFi08GRkAAAA+wN/bAa7m+PHjql+/vvz9/xs1LCxM\nxcXFOnnypEJDQx3jffr00ZYtWzRkyBD5+fmpVq1aWrhwoYKDg70RHQAAAF7k81d0CwsLFRAQUGbs\nwmObzVZm/NSpU7JarUpNTdXKlSuVkJCgiRMn6sSJEx7LCwAAAN/g80U3MDCwXKG98DgoKKjM+MyZ\nM9WmTRsNHjxYbdu21dSpUxUUFKTVq1d7LC8AAAB8g88X3SZNmujUqVMqLS11jFmtVpnNZoWEhJSZ\nu3fvXlksFsdjk8kki8Wio0ePeiwvAAAAfIPPF93w8HD5+/tr586djrHt27erXbt25eY2bty43B0W\nvv32WzVr1sztOQEAAOBbfL7oms1mxcfHKzU1Vbt379bmzZu1ZMkSDRs2TNJvV3cv3Cd34MCBWrly\npdasWaPDhw9r5syZOnbsmBISErz5FgAAAOAFPn/XBUl69tlnNWXKFA0bNkzBwcF64oknFBcXJ0nq\n0qWL0tLSlJCQoD59+qiwsFALFixQfn6+wsPDtWzZMjVo0MDL7wAAAACeZrLb7XZvh/AVsbGxkqSs\nrCwvJwEAAEB1+fzSBQAAAKAqKLoAAAAwJIouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAA\nDImiCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImi\nCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAA\nAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwJP/KTszOztaKFSv0yy+/\n6L777lOfPn3cmQsAAAColkpd0f3888/14IMP6tChQyouLtb48eM1ffp0d2cDAAAAqqxSV3TnzJmj\noUOHauLEiZKkLVu2aOzYsTp37pyeeeYZmc1mt4YEAAAAnFWpK7oHDhzQ4MGDHY979Oih2bNn6/33\n31fnzp313HPPyWq1Kjw83G1BAQAAAGdU6opu/fr19f3336tly5aOsZ49eyoyMlKfffaZAgICVKdO\nHT3++ONuCwoAAAA4w2S32+1Xm5SWlqZ169bp6aefVlxcnOrWreuJbB4XGxsrScrKyvJyEgAAAFRX\npa7oPvHEE/r55581adIk1a1bV3Fxce7OBQAAAFRLpa7oXnDy5EkFBATouuuuu+x2u90uk8nksnCe\nxhVdAAAA43DqCyNCQ0PVr18/nTp1qty2/Px8de7c2WXBAAAAgOqo1NKF9evX65NPPpEk/fDDD5o6\ndaoCAwPLzPnhhx9q9NVcAAAAGEulim5UVJTefvttXVjlcPToUdWuXdux3WQyqU6dOpoxY4Z7UgIA\nAABOqlTRbdq0qZYtWyZJGjp0qObOnat69eq5NRgAAABQHZUquhdLT093Rw4AAADApZwuut98842m\nTp2qHTt26Ny5c+W279u3zyXBAAAAgOpwuuimpqbq559/1tNPP62QkBB3ZAIAAACqzemi+9VXX2n5\n8uW69dZb3ZEHAAAAcAmn7qMr/XYv3YvvuAAAAAD4IqeL7oMPPqhZs2apoKDAHXkAAAAAl3B66cLn\nn3+u7du364477lBYWJgCAgLKbOfrcwEAAOALnC660dHRio6OdkeWCtlsNk2ePFkffvihzGazHnnk\nET388MOXnbt//35NmTJFe/fuVcuWLfXcc8+pU6dOHs0LAAAA73O66I4ZM8YdOa5oxowZysnJUXp6\nuvLy8jRhwgTdcMMN6tWrV5l5BQUFGj58uGJjYzVjxgy9//77GjNmjDZu3KgGDRp4PDcAAAC8x+k1\nupKUm5urZ599VoMGDVJ+fr7eeustbd261dXZJEmFhYVatWqVnn/+eVksFsXFxSk5OVkZGRnl5q5e\nvVrXXXedpkyZoubNmyslJUWtWrXSnj173JINAAAAvsvportnzx4NHDhQeXl52rNnj2w2m/bt26dH\nHnlEH3/8scsD5ubmqqSkRJGRkY6x6Oho7dq1q9zcbdu2qUePHmXGVq5cqW7durk8FwAAAHyb00V3\n5syZeuSRR5Senu64zdhLL72kpKQk/e1vf3N5wOPHj6t+/fry9//vKouwsDAVFxfr5MmTZeYeOXJE\noaGheuGFF9SlSxcNGjRIO3bscHkmAAAA+L4qXdFNSEgoN56UlKRDhw65JNTFCgsLy93Z4cJjm81W\nZvzs2bNatGiRGjdurEWLFqlDhw4aPny48vPzXZ4LAAAAvs3polu7du3L3kP32LFjCgoKckmoiwUG\nBpYrtBceX7o/Pz8/hYeHa8yYMbJYLBo/frxatWqlNWvWuDwXAAAAfJvTRTcuLk6vvfaaTp8+7Rg7\ndOiQXn75Zd11112uzCZJatKkiU6dOqXS0lLHmNVqldlsVkhISJm5jRo1UuvWrcuMtWrVSseOHXN5\nLgAAAPg2p4vuhAkTdObMGXXu3FmFhYUaMGCA7rvvPvn5+emZZ55xecDw8HD5+/tr586djrHt27er\nXbt25eZGRkYqNze3zNg333yjG264weW5AAAA4Nucvo9u3bp19fbbb+uLL75QTk6OSktLdcstt6hr\n166qVatKdyu7IrPZrPj4eKWmpmratGnKz8/XkiVLlJaWJum3q7vBwcEKDAzUoEGDlJGRoblz56pf\nv3567733lJeXp379+rk8FwAAAHybyW6326vzAidOnNDWrVvVrl07NWvWzFW5yigqKtKUKVO0ceNG\nBQcHKzk5WUOHDpUkWSwWpaWlOT4gl52drRdffFGHDh3SjTfeqOeff17t27ev1H5iY2Ml8TXGAAAA\nRuB00f3666+VkpKil156SW3atFHfvn11/PhxBQQEaOHChercubO7srodRRcAAMA4nF5rMGPGDLVs\n2VKtW7fW2rVrde7cOX388ccaPny4XnvtNXdkBAAAAJzmdNHNzs7WhAkTFBYWpk8++UTdu3dXkyZN\nNGDAgHIfBAMAAAC8xemiW6tWLQUEBOj8+fPaunWrYmJiJElnzpyR2Wx2eUAAAACgKpy+60JkZKQW\nLFigBg0aqLi4WN26dVN+fr5mzZqlyMhId2QEAAAAnOb0Fd2//OUvysnJ0fLlyzVp0iQ1aNBACxcu\n1KFDh9xyH10AAACgKqp9ezHpt1uM1atXT35+fq7I5DXcdQEAAMA4qvQND0VFRXr//ff16quv6tSp\nUzp48GCZrwQGAAAAvM3pNbpWq1WDBg2S1WqVzWZTYmKi3nzzTe3Zs0dLly7VjTfe6I6cAAAAgFOc\nvqKblpamm266SV988YUCAwMl/XZv3ZtvvlmvvPKKywMCAAAAVeF00f3yyy81duxYBQUFOcbq1aun\nCRMmaMeOHS4NBwAAAFSV00X3zJkzqlOnzmW3nT9/vtqBAAAAAFdwuuh27NhRy5cvLzN27tw5zZs3\nT+3bt3dZMAAAAKA6nP4w2oQJE5SUlKStW7fq3Llzmjx5sr755hv9+uuvysjIcEdGAAAAwGlOF90b\nb7xRmZmZWr58uRo3bqzS0lL17t1bQ4YMUbNmzdyREQAAAHCa00V3zJgxGjdunJ544gl35AEAAABc\nokp3XbhwWzEAAADAVzlddPv376+ZM2fqwIEDstls7sgEAAAAVJvTSxc+/vhjHT58WBs3brzs9n37\n9lU7FAAAAFBdThfd0aNHuyMHAAAA4FJOF93+/fu7IwcAAADgUk4XXUnasGGDli5dqq+//lp+fn5q\n27atRowYoS5durg6HwAAAFAlTn8YbdWqVXr66ad1/fXXa9y4cXr88cdVr149jRo1Sps3b3ZHRgAA\nAMBpJrvdbnfmCb169dKQIUP00EMPlRlftGiRMjMzlZmZ6cp8HhUbGytJysrK8nISAAAAVJfTV3Tz\n8/N11113lRvv2bOnvv/+e1dkAgAAAKrN6aLboUMHrV+/vtz4p59+qujoaJeEAgAAAKrL6Q+jdejQ\nQfPmzdOePXt0xx13qHbt2tq9e7fWrl2rAQMGaO7cuY65Y8aMcWlYAAAAoLKcXqPbo0ePyr2wyVTj\n1rqyRhcAAMA4nL6iu2XLFnfkAAAAAFzK6TW6AAAAQE1A0QUAAIAhUXQBAABgSBRdAAAAGFK1i+65\nc+e0e/dunTlzxhV5AAAAAJdwuugeO3ZMjzzyiHbt2qWioiL1799fAwcOVI8ePbRv3z53ZAQAAACc\n5nTRnT59un799Vc1aNBAGzZs0NGjR/WPf/xDPXv21CuvvOKOjAAAAIDTnL6P7pdffqmlS5eqWbNm\nmjlzprp27ar27dsrNDRUAwYMcEdGAAAAwGlOX9E9d+6c6tWrJ7vdri+++EJ33nmnJKm0tFT+/k73\nZgAAAMAtnG6mbdu21apVq9SoUSOdPn1a3bt3l81m0xtvvCGLxeKOjAAAAIDTnC66EyZM0KOPPqqT\nJ09qxIgR+t3vfqfJkycrKytLixYtckdGAAAAwGkmu91ud+YJhw8fVrNmzVRQUKCQkBBJ0rfffqvQ\n0FDVr1/fLSE9JTY2VpKUlZXl5SQAAACoLqfX6CYlJWnPnj2OkitJv//972t8yQUAAICxOF10a9eu\nzYfOAAAA4POcbqz9+/dXcnKy4uPj1bJlS5nN5jLbExISXBYOAAAAqCqn1+he6c4KJpOpRn87Gmt0\nAQAAjMPpK7q5ubnuyAEAAAC4lNNrdC84evSoPvnkExUVFennn392ZSYAAACg2py+omuz2TRhwgRt\n2LBBtWrV0saNGzVjxgwVFBRo7ty5qlu3rjtyAgAAAE5x+oruvHnzlJubq6VLlyowMFCSNHToUB0+\nfFgzZ850eUAAAACgKpwuuuvWrdNf/vIXderUyTHWqVMnvfzyy3yICwAAAD7D6aKbn5+vFi1alBtv\n2rSpfvnlF5eEupTNZtOkSZPUsWNHde3aVUuWLLnqc/Ly8hQVFaVt27a5JRMAAAB8m9NF98Ybb9QX\nX3xRbnzdunW66aabXBLqUjNmzFBOTo7S09OVmpqquXPnatOmTVd8zuTJk1VUVOSWPAAAAPB9Tn8Y\nLSUlRePGjdPBgwdVUlKi9957T99++602btyo2bNnuzxgYWGhVq1apcWLF8tischisSg5OVkZGRnq\n1avXZZ+TmZmps2fPujwLAAAAag6nr+jefffd+t///V/t2bNHfn5+Wrx4sY4cOaLZs2frnnvucXnA\n3NxclZSUKDIy0jEWHR2tXbt2XXb+yZMn9eqrr2rq1Kly8rswAAAAYCBOX9GVpG7duqlbt26uznJZ\nx48fV/369eXv/9+oYWFhKi4u1smTJxUaGlpmflpamvr37++2ZRQAAACoGSpVdN9///1Kv2BCQkKV\nw1xOYWGhAgICyoxdeGyz2cqMf/7558rOztaLL77o0gwAAACoeSpVdCdOnFjmsclkkt1uV1BQkPz9\n/fXrr7/Kz89PoaGhLi+6gYGB5QrthcdBQUGOseLiYk2ePFmpqanlijEAAACuPZUqurm5uY7fr127\nVosXL9b06dNlsVgkSd99950mTpyovn37ujxgkyZNdOrUKZWWlqpWrd+WFFutVpnNZoWEhDjm7dq1\nS0eOHFFKSkqZtbkjRoxQQkKCJk+e7PJsAAAA8F1Or9GdOXOm5syZ4yi5ktSqVSs999xzGj16tIYO\nHerSgOHh4fL399fOnTvVvn17SdL27dvVrl27MvMiIiLK3XKsZ8+eevnllxUTE+PSTAAAAPB9Thfd\n06dPO77692KlpaVuuW+t2WxWfHy8UlNTNW3aNOXn52vJkiVKS0uT9NvV3eDgYAUGBqp58+blnt+4\ncWM1aNDA5bkAAADg25y+vVinTp00depU5eXlOcYOHTqkKVOm6K677nJlNodnn31W7dq107Bhw/Ti\niy/qiSeeUFxcnCSpS5cu2rBhw2WfZzKZ3JIHAAAAvs9kd/Jms/n5+Ro+fLgOHTrkWCN7+vRp3X77\n7Vq4cKHq1avnlqCeEBsbK0nKysrychIAAABUl9NLF5o0aaI1a9bo888/14EDByT9to62c+fOXEEF\nAACAz6jSF0b4+fmpa9eu6tq1q6vzAAAAAC5RqaIbHh6uTz/9VGFhYbJYLFe8crtv3z6XhQMAAACq\nqlJFd9q0aQoODnb8niUKAAAA8HVOfxjNyPgwGgAAgHFUaY3uhg0btHTpUn399dfy8/NT27ZtNWLE\nCHXp0sXV+QAAAIAqcfo+uqtWrdLTTz+t66+/XuPGjdPjjz+uevXqadSoUdq8ebM7MgIAAABOc3rp\nQq9evTRkyBA99NBDZcYXLVqkzMxMZWZmujKfR7F0AQAAwDicvqKbn59/2W9A69mzp77//ntXZAIA\nAACqzemi26FDB61fv77c+Keffqro6GiXhAIAAACqy+kPo3Xo0EHz5s3Tnj17dMcdd6h27dravXu3\n1q5dqwEDBmju3LmOuWPGjHFpWAAAAKCynF6j26NHj8q9sMlU49a6skYXAADAOJy+ortlyxZ35AAA\nAABcyuk1uhc7ceKENm3apB07drgqDwAAAOASlS66f//739WpUyfHnRV27NihXr16aezYsRoyZIge\nfvhhFRUVuS0oAAAA4IxKFd133nlH8+fPV2JiosLCwiRJkyZNktls1tq1a/Xxxx/rzJkzWrhwoVvD\nAgAAAJVVqaK7cuVKTZw4UU8//bTq1q2r3bt367vvvtPQoUN10003qUmTJho9erTWrVvn7rwAAABA\npVSq6B46dEh/+MMfHI+//PJLmUwmde/e3TF200036ejRo65PCAAAAFRBpdfomkwmx++3b9+uevXq\nyWKxOMbOnDmjoKAg16YDAAAAqqhSRfeWW25x3Fnh9OnT+ve//13mCq8kbdiwQbfccovrEwIAAABV\nUKn76CYlJSk1NVX79u1Tdna2bDabhg0bJknKz8/XBx98oMWLF+vll192a1gAAACgsipVdPv16yeb\nzably5erVq1amj17tm6//XZJ0oIFC7RixQqNGDFC8fHxbg0LAAAAVJbTXwF8qfz8fAUEBCg0NNRV\nmbyGrwAGAAAwDqe/AvhSTZo0cUUOAAAAwKWq9RXAAAAAgK+i6AIAAMCQKLoAAAAwJIouAAAADImi\nCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAA\nAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAAAEOi\n6AIAAMCQKLoAAAAwJIouAAAADKlGFF2bzaZJkyapY8eO6tq1q5YsWVLh3I8++kgJCQmKiopSfHy8\ntmzZ4sGkAAAA8BU1oujOmDFDOTk5Sk9PV2pqqubOnatNmzaVm7d//36lpKRo4MCByszMVGJiosaO\nHav9+/d7ITUAAAC8yeeLbmFhoVatWqXnn39eFotFcXFxSk5OVkZGRrm5a9euVUxMjJKSktS8eXMl\nJSWpU6fQ2BLGAAAa2klEQVRO2rBhgxeSAwAAwJv8vR3ganJzc1VSUqLIyEjHWHR0tBYsWFBubv/+\n/XXu3Lly4wUFBW7NCAAAAN/j81d0jx8/rvr168vf/7+dPCwsTMXFxTp58mSZua1bt1abNm0cjw8c\nOKAvv/xSMTExHssLAAAA3+DzRbewsFABAQFlxi48ttlsFT7vxIkTSklJUXR0tGJjY92aEQAAAL7H\n54tuYGBguUJ74XFQUNBln2O1WjVs2DCZTCbNmTPH7RkBAADge3y+6DZp0kSnTp1SaWmpY8xqtcps\nNiskJKTc/Pz8fCUlJamkpETp6ekKDQ31ZFwAAAD4CJ8vuuHh4fL399fOnTsdY9u3b1e7du3KzS0s\nLFRycrJq166tjIwMNWzY0JNRAQAA4EN8vuiazWbFx8crNTVVu3fv1ubNm7VkyRINGzZM0m9Xd4uL\niyVJ8+fPV15enqZPn67S0lJZrVZZrVbuugAAAHANMtntdru3Q1xNUVGRpkyZoo0bNyo4OFjJycka\nOnSoJMlisSgtLU0JCQnq3bu3vvvuu3LPT0hI0PTp06+6nwsfWsvKynJpfgAAAHhejSi6nkLRBQAA\nMA6fX7oAAAAAVAVFFwAAAIZE0QUAAIAhUXQBAABgSBRdAAAAGBJFFwAAAIZE0QUAAIAhUXQBAABg\nSBRdAAAAGBJFFwAAAIZE0QUAAIAhUXQBAABgSBRdAAAAGBJFFwAAAIZE0QUAAIAhUXQBAABgSBRd\nAAAAGBJFFwAAAIZE0QUAAIAhUXQBAABgSBRdAAAAGBJFFwAAAIZE0QUAAIAhUXQBAABgSBRdAAAA\nGBJFFwAAAIZE0QUAAIAhUXQBAABgSBRdAAAAGBJFFwAAAIZE0QUAAIAhUXQBAABgSBRdAAAAGBJF\nFwAAAIZE0QUAAIAhUXQBAABgSBRdAAAAGBJFFwAAAIZE0QUAAIAhUXQBAABgSBRdAAAAGBJFFwAA\nAIZE0QUAAIAhUXQBAABgSBRdAAAAGBJFFwAAAIZE0QUAAIAhUXQBAABgSBRdAAAAGFKNKLo2m02T\nJk1Sx44d1bVrVy1ZsqTCuTk5OUpMTFRkZKQGDhyovXv3ejApAAAAfEWNKLozZsxQTk6O0tPTlZqa\nqrlz52rTpk3l5hUWFmrkyJHq2LGjVq9ercjISI0aNUpFRUVeSA0AAABv8vmiW1hYqFWrVun555+X\nxWJRXFyckpOTlZGRUW7uunXrFBQUpD//+c9q3bq1nnvuOV133XX65z//6YXkAAAA8CZ/bwe4mtzc\nXJWUlCgyMtIxFh0drQULFpSbu2vXLkVHR5cZa9++vbKzs5WQkOD2rLi6HTu/0st/naevvzum2oF1\n1LpZmCY9PULtIyO8Hc1rduz8SvMWr9DpwlKFBNXS3V2i9K9Psx2PRw9PdMnxuXQ/F7/ulbZV9TVx\ndZc7fpI4pgDgIj5fdI8fP6769evL3/+/UcPCwlRcXKyTJ08qNDTUMf7TTz/plltuKfP8sLAwHTx4\n0GN5UbEdO7/SYxPn6GyJWS07PyyTyaQiu11PTn5Tr01+5Jr8y3zHzq/0zLR0mZt1kSnApB+Pf68v\nF2xQy8j7ZQow6Yzdrmempeuvk1St43Ppfi5+XUkVbrvSPq/0mtfiuXTW5Y7fYxPnKMAcrJDWd3NM\nAcAFasTShYCAgDJjFx7bbLYy40VFRZede+k8eMdvV6nsanFbL5lMJkmSyWRSSOu7NW/xCi+n8455\ni1f8VnT+/3gc/z77t5J70fExN+tS7eNz6X4uft0rbavqa+LqLnf8Thfafyu5HFMAcAmfL7qBgYHl\niuqFx0FBQZWaazab3RsSlXK6sFR+/gGOv8Qv+O0v+FIvpfKu04WlZY6Hu47Ppfu5+HWvtK2qr4mr\nu9zx4+cDAFzL54tukyZNdOrUKZWW/vcPeqvVKrPZrJCQkHJzjx8/XmbMarWqUaNGHsmKKwsJqqWS\n8zbZ7fYy43a7XSFBPv+foluEBNUqczzcdXwu3c/Fr3ulbVV9TVzd5Y4fPx8A4Fo+/6dneHi4/P39\ntXPnTsfY9u3b1a5du3JzIyIilJ2dXWYsOzu7zAfZ4D2jhycqJMikw7s3Of4yt9vtOv3NvxwfwrnW\njB6eqKK8Tx3Ho1HLKH2/84Myx6co79NqH59L93Px615pW1VfE1d3ueMXEmTS6W/+xTEFABcx2S+9\nfOCDUlNTtWPHDk2bNk35+fmaOHGi0tLSFBcXJ6vVquDgYAUGBqqgoED33HOP+vbtqwceeEDLly/X\nxo0b9eGHH1Zq+UJsbKwkKSsry91v6Zq1Y+dXmvbKPH397TH5B9bRjc3C9Cx3XdD8N1fql7Mljrsu\nfPTZTsdjV9514eL9XHrXhYq2VfU1cXWXO36SOKYA4CI1ougWFRVpypQp2rhxo4KDg5WcnKyhQ4dK\nkiwWi9LS0hy3D9u9e7dSU1P1zTffqE2bNpoyZYosFkul9kPRBQAAMI4aUXQ9haILAABgHD6/RhcA\nAACoCoouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAw\nJIouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwJIou\nAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAA\nDImiCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImi\nCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAA\nAEOi6AIAAMCQKLoAAAAwJIouAAAADImiCwAAAEOi6AIAAMCQKLoAAAAwpBpRdGfOnKmYmBh16tRJ\nr7zyyhXn7ty5U4MGDVJUVJR69+6tlStXeiglAAAAfIm/twNczZtvvql169bp9ddf17lz5zR+/Hg1\nbNhQDz/8cLm5VqtVI0eO1JAhQ/TXv/5Ve/bs0bPPPqvGjRure/fuXkgPAAAAb/H5K7rp6el64okn\nFBUVpTvuuEPjx49XRkbGZedu3rxZjRo10pNPPqkWLVqoT58+io+P19q1az2cGgAAAN7m01d0f/rp\nJx07dkwdOnRwjEVHR+vo0aOyWq1q2LBhmfndunVT27Zty73Or7/+6vasAAAA8C0+fUX3+PHjMplM\naty4sWOsYcOGstvt+vHHH8vNv/7663X77bc7Hv/8889av3697rzzTo/kBQAAgO/w+hXd4uJi5efn\nX3bb2bNnJUkBAQGOsQu/t9lsV33dlJQUNW7cWA888EClshw/flznz59XbGxspeYDAADA85o2bVrh\nUtaLeb3ofvXVV/rTn/4kk8lUbtv48eMl/VZqLy24QUFBFb7m2bNnNXr0aB0+fFjLly9XYGBgpbIE\nBATIbrc7+xYAAADgg0x2H252P/30k7p3766srCxdf/31kqS8vDz17NlTn3zySbk1upJUUFCg5ORk\n5eXlaenSpbrxxhs9HRsAAAA+wKfX6DZu3FhNmzbVf/7zH8fY9u3b1bRp08uWXLvdrjFjxuiHH35Q\nRkYGJRcAAOAa5vWlC1czaNAgzZw5U02aNJHdbtesWbM0fPhwx/YTJ07IbDarTp06WrlypbZu3ap5\n8+apbt26slqtkqTatWurXr163noLAAAA8AKfXrogSaWlpXrllVe0evVq1apVS4mJiRo3bpxje48e\nPTRgwACNGTNGycnJ+uyzz8q9RseOHbVs2TJPxgYAAICX+XzRBQAAAKrCp9foAgAAAFVF0QUAAIAh\nUXQBAABgSBRdAAAAGBJFFwAAAIZE0b1G2Gw2TZo0SR07dlTXrl21ZMmSCud+9NFHSkhIUFRUlOLj\n47VlyxYPJoUrOHO+L8jLy1NUVJS2bdvmgYRwJWfO9/79+zVkyBBFRESoX79++ve//+3BpHAFZ873\nhx9+qL59+yoqKkpJSUnKycnxYFK4ks1m0/3333/FP6NzcnKUmJioyMhIDRw4UHv37vVgQt9E0b1G\nzJgxQzk5OUpPT1dqaqrmzp2rTZs2lZu3f/9+paSkaODAgcrMzFRiYqLGjh2r/fv3eyE1qqqy5/ti\nkydPVlFRkYcSwpUqe74LCgo0fPhw3XzzzVq7dq169uypMWPG6MSJE15Ijaqq7Pk+ePCgxo8fr1Gj\nRikzM1MWi0UjR45UcXGxF1KjOmw2m5566ikdPHiwwjmFhYUaOXKkOnbsqNWrVysyMlKjRo3iz3U7\nDO/s2bP222+/3b5t2zbH2Ouvv24fOnRoubkzZ860jxgxoszYI488Yp89e7bbc8I1nDnfF6xZs8Y+\nePBgu8VisW/dutUTMeEizpzvpUuX2nv16lVm7I9//KP9448/dntOuIYz53vJkiX2//mf/3E8Ligo\nsLdp08a+Z88ej2SFaxw8eNAeHx9vj4+Pv+Kf0StXrrTHxcWVGevVq5f9vffe80RMn8UV3WtAbm6u\nSkpKFBkZ6RiLjo7Wrl27ys3t37+/nn766XLjBQUFbs0I13HmfEvSyZMn9eqrr2rq1Kmy8/0xNY4z\n53vbtm3q0aNHmbGVK1eqW7dubs8J13DmfNevX18HDx7Ujh07ZLfb9e677yo4OFgtWrTwZGRU09at\nWxUTE6N33nnnin9G79q1S9HR0WXG2rdvr+zsbHdH9Gn+3g4A9zt+/Ljq168vf///nu6wsDAVFxfr\n5MmTCg0NdYy3bt26zHMPHDigL7/8UkOGDPFYXlSPM+dbktLS0tS/f3/ddNNNno4KF3DmfB85ckS3\n3XabXnjhBW3ZskXNmjXTM888o/bt23sjOqrAmfPdp08fbdmyRUOGDJGfn59q1aqlhQsXKjg42BvR\nUUWDBw+u1LyffvpJt9xyS5mxsLCwKy53uBZwRfcaUFhYqICAgDJjFx7bbLYKn3fixAmlpKQoOjpa\nsbGxbs0I13HmfH/++efKzs7WY4895rF8cC1nzvfZs2e1aNEiNW7cWIsWLVKHDh00fPhw5efneywv\nqseZ833q1ClZrValpqZq5cqVSkhI0MSJE1mTbVBFRUWX/W/jSn/PXwsouteAwMDAcv+hX3gcFBR0\n2edYrVYNGzZMJpNJc+bMcXtGuE5lz3dxcbEmT56s1NTUcn84ouZw5ufbz89P4eHhGjNmjCwWi8aP\nH69WrVppzZo1HsuL6nHmfM+cOVNt2rTR4MGD1bZtW02dOlVBQUFavXq1x/LCcyr6b8NsNnspkW+g\n6F4DmjRpolOnTqm0tNQxZrVaZTabFRISUm5+fn6+kpKSVFJSovT09HL/qxu+rbLne9euXTpy5IhS\nUlIUFRWlqKgoSdKIESM0efJkT8dGFTnz892oUaNyy5NatWqlY8eOeSQrqs+Z8713715ZLBbHY5PJ\nJIvFoqNHj3osLzynSZMmOn78eJkxq9WqRo0aeSmRb6DoXgPCw8Pl7++vnTt3Osa2b9+udu3alZtb\nWFio5ORk1a5dWxkZGWrYsKEno8IFKnu+IyIitGnTJq1Zs0aZmZnKzMyUJL388ssaO3asRzOj6pz5\n+Y6MjFRubm6ZsW+++UY33HCD23PCNZw5340bNy63PvPbb79Vs2bN3J4TnhcREVHug2fZ2dllPrh4\nLaLoXgPMZrPi4+OVmpqq3bt3a/PmzVqyZImGDRsm6bd/8V24r+L8+fOVl5en6dOnq7S0VFarVVar\nlbsu1CCVPd8BAQFq3rx5mV/Sb385NmjQwJtvAU5w5ud70KBB2r9/v+bOnavDhw9rzpw5ysvLU79+\n/bz5FuAEZ873wIEDtXLlSq1Zs0aHDx/WzJkzdezYMSUkJHjzLcCFLj7f99xzj3799VdNmzZNhw4d\n0ksvvaSzZ8+qd+/eXk7pZV6+vRk8pLCw0D5x4kR7VFSUvVu3bvZly5Y5trVp08Zxn717773XbrFY\nyv2aOHGit6KjCip7vi/FfXRrJmfO944dO+z9+/e333777fb+/fvb//Of/3gjMqrBmfO9atUqe+/e\nve3t27e3JyUl2fft2+eNyHCRS/+MvvR879q1y96/f397RESEPTExkfNtt9tNdjs3zgQAAIDxsHQB\nAAAAhkTRBQAAgCFRdAEAAGBIFF0AAAAYEkUXAAAAhkTRBQAAgCFRdAEAAGBIFF0AAAAYEkUXAAAA\nhkTRBYDLWL16tYYOHaqYmBjddttt6tWrl6ZNmyar1Vrt1547d65iY2NdkPLynn32Wf3pT39y2+tf\n8Pjjj2vz5s1u388Fy5Yt07Rp0zy2PwA1H18BDAAXsdvtevzxx/Wf//xHo0ePVteuXXXdddfpwIED\nev3113X06FG99957atCgQZX3UVhYqKKiIoWGhrow+X8VFBSotLRUISEhbnl9SVq7dq1WrFihZcuW\nuW0flyopKVHfvn01bdo0tW/f3mP7BVBzUXQB4CJvvvmmZs2apVWrVslisZTZVlxcrPvuu0/33HOP\nxo8f76WE3ldaWqp77rlHL7zwgrp27erRfb/11lvasGGDMjIyPLpfADUTSxcA4CIZGRlKSEgoV3Il\nKTAwUMuWLdMTTzwhSfrhhx9ksVi0cOFC/eEPf1DPnj115swZff3113r00Ud1xx13qF27doqLi9OS\nJUscr/O3v/1NPXr0KPMamzZtUmJiom6//Xb16NFDK1asqDBjaWmpXnnlFd1111267bbb1Lt3b739\n9tuO7RMnTnQsXRg6dKgsFku5X8OGDXPMf/fdd9WnTx9FRESob9++WrZsma50DWTjxo06ffq07rzz\nTsdYjx499MYbb2jUqFGKjIxUjx49tHnzZmVlZenee+9VVFSUkpOTdeLECUnS1q1bdeutt2rz5s26\n9957FRERoYceekg//vijXnrpJXXs2FF33nmn5s+fX2bf9957r7Kzs7Vnz54K8wHABRRdAPh/R44c\n0dGjRxUTE1PhnKZNm6p27dplxt5//30tW7ZMr732mvz8/PTII48oNDRU77zzjtavX6/evXtrxowZ\nys3NlSSZTCaZTKYyr5GWlqbHHntM69ev1913360pU6bohx9+uGyGt956S5s2bdKcOXO0adMmPfjg\ng5oyZYp27NjheP0L/v73v+uzzz5z/Hruuefk7++v0aNHS5LeeecdvfLKK0pJSdG6dev05JNP6o03\n3tCrr75a4THIysrSnXfeKT8/vzLjr7/+uvr27asPPvhA4eHhmjBhghYsWKBXX31VCxYs0O7du/XG\nG2845peUlGj+/PmaNWuWli1bpn379ik+Pl6BgYFatWqVBg0apNdee00HDhxwPCcsLEzt2rVTVlZW\nhfkA4AKKLgD8v59//lmSyq2/ffTRRxUVFeX4df/995fZnpSUpBtvvFG33nqrzp49q4ceekgvvPCC\nfv/736tFixYaM2aMJOnrr7+ucN8PP/yw7rrrLjVr1kxPPvmkSkpK9NVXX1127pEjRxQUFKTrr79e\nTZs2VVJSkt588021atWq3NyQkBCFhYUpLCxMR44c0cyZM/XCCy+oc+fOkqR58+bpscceU+/evdWs\nWTP17NlT48aNU3p6umw222X3/9VXX+mWW24pN3733XerX79+at68uRITE3X27FmNGzdOt956q+64\n4w7deeedZUqrJD355JNq27atIiIi1LlzZ9WpU0d//vOf1bJlS40aNUqSyj3n5ptv1s6dOys8lgBw\ngb+3AwCAr7jw4bBTp06VGX/xxRdVVFQkSVq6dKn+9a9/ldneokULx+8bNGigwYMH64MPPlBOTo4O\nHz6s3NxcmUwmlZaWVrjv1q1bO34fHBwsSRUWzaSkJG3evFndu3dXeHi4/vCHP6hPnz5X/IBcXl6e\nHn/8cQ0ePFiJiYmSpBMnTujHH3/UrFmzNHv2bMdcu92uc+fOKS8vr0yuC6xWq8LCwsqNX3wcgoKC\nJEnNmzd3jJnNZsc/JqTfrjxf/Jw6deqoWbNmjseBgYGXPQ4NGjSo8B8BAHAxii4A/L/mzZurUaNG\n2rp1q3r37u0Yb9SokeP39evXL/c8s9ns+L3ValViYqIaNmyoHj16qEuXLrrtttvUvXv3K+47ICCg\n0jlbtmypDz/8UFu3btVnn32mjz76SG+88YamT5+uhISEcvMLCgo0atQo3XbbbZowYYJj/MI63EmT\nJl12uUbTpk0vu/+KSru/f/m/UmrVuvL/OLz0OZcu6bic8+fPX/V1AUBi6QIAONSqVUtDhw7V+++/\nr/379192ztGjR6/4GmvXrtXp06f19ttv69FHH1VcXJzjCrGrbnKTnp6ujRs3KiYmRuPHj1dmZqZi\nYmK0YcOGcnNLSko0duxY+fn5adasWWWK5IUlDYcPH1bz5s0dv3bv3q3Zs2dXmLdRo0aOD5V5w8mT\nJ8v84wMAKsIVXQC4yIgRI5Sbm6ukpCSNGDFC3bt3V3BwsPbv36+33npLn3/+uf74xz9W+Pzf/e53\nKiws1Pr16xUdHa1Dhw4pLS1NJpOpwqUIzjpx4oRef/11mc1mWSwWHTp0SPv27dNDDz1Ubu6UKVO0\nf/9+vfnmmyosLNTZs2cd2xo2bKjk5GS99tpratq0qbp166bc3FxNmTJFcXFx5T50d0FERIT27t17\n1ZxXK/ZVLf579+5Vz549q/RcANcWii4AXMRkMmnWrFnauHGj3n33XaWnp+uXX35Ro0aN1KFDB2Vk\nZCg6OrrM/Ivde++9ysnJUVpams6cOaMbbrhBf/zjH5WVlaXdu3frgQceuOw+KzN2QUpKis6fP6+X\nXnpJVqtVDRs2VFJSkkaOHFnu+StWrJDJZCqzpMFut8tkMmnfvn16+OGHZTablZ6errS0NDVq1EiD\nBg1yfIDucuLi4vSXv/xFJSUljjsvOPseKrP9cnNOnDihgwcPasaMGVd9LgDwhREAAKecP39e9957\nr5555hn16tXLo/tevHixtmzZorfeesuj+wVQM7FGFwDgFH9/f40ZM6bMl2B4gs1m09tvv61x48Z5\ndL8Aai6KLgDAaQkJCapXr542bdrksX3+4x//UPfu3dWhQweP7RNAzcbSBQAAABgSV3QBAABgSBRd\nAAAAGBJFFwAAAIZE0QUAAIAhUXQBAABgSBRdAAAAGBJFFwAAAIZE0QUAAIAh/R89zm6wtMqlogAA\nAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "df.plot.scatter('Grain size (mm)', 'Spiders')\n", "plt.ylabel('Spiders present?')\n", "sns.despine()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I will analyse this with the *scikit-learn* package." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import sklearn.linear_model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "scikit-learn has a logisitic regression classifier which uses regularization. To eliminate regularization, we set the regularization parameter `C` to $10^{12}$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-1.64748882] [[ 5.12118846]]\n" ] } ], "source": [ "# C=1e12 is effectively no regularization - see https://github.com/scikit-learn/scikit-learn/issues/6738\n", "clf = sklearn.linear_model.LogisticRegression(C=1e12, random_state=0)\n", "clf.fit(df['Grain size (mm)'].reshape(-1, 1), df['Spiders'])\n", "print(clf.intercept_, clf.coef_)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is in agreement with the equation John reports:\n", "$$\n", "probability of spider presence = \\frac{e^{-1.6476+5.1215(grain \\; size)}}{(1+e^{-1.6476+5.1215(grain \\; size)}}\n", "$$" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def plot_log_reg(x, y, data, clf, xmin=None, xmax=None, alpha=1, ax=None):\n", " if ax is None:\n", " fig, ax = plt.subplots()\n", " else:\n", " fig = ax.figure\n", " ax.scatter(data[x], data[y], color='black', zorder=20, alpha=alpha)\n", " if xmin is None:\n", " xmin = x.min()\n", " if xmax is None:\n", " xmax = x.max()\n", " X_test = np.linspace(xmin, xmax, 300)\n", "\n", " loss = scipy.special.expit(X_test * clf.coef_ + clf.intercept_).ravel()\n", " ax.plot(X_test, loss, linewidth=3)\n", "\n", " ax.set_xlabel(x)\n", " ax.set_ylabel(y)\n", " fig.tight_layout()\n", " sns.despine()\n", " return fig, ax" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxAAAAIYCAYAAAD93QokAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xd4lHXe/v1zSjIJ6Qkk9F4SOgSkCBaaq64UC6uiYkEs\nK+vqug+iqyCuLvxWb2933XtFQFRQV8GGsAgK6or0XgOEHgghIaGkTjJzPX9Q1mECDCGZayZ5v44j\nxySf7zWTM6AwJ1ezGIZhCAAAAAB8YDU7AAAAAIDgQYEAAAAA4DMKBAAAAACfUSAAAAAA+IwCAQAA\nAMBnFAgAAAAAPqNAAAAAAPAZBQIAAACAzygQAAAAAHwWVAXC6XTqlltu0erVqy+4zQ8//KChQ4eq\nS5cuGjJkiJYsWeLHhAAAAED1FjQFwul06umnn1Z6evoFt9mxY4fGjBmjO+64Q3PnztXw4cP1u9/9\nTjt27PBjUgAAAKD6CooCsXv3bg0fPlwZGRkX3W7evHnq1auXRowYoUaNGmnEiBHq0aOHFixY4Kek\nAAAAQPVmNzuAL1atWqVevXrp97//vTp16nTB7YYNG6bS0lKveX5+flXGAwAAAGqMoCgQd911l0/b\nNW/e3OPrXbt2acWKFbr77rurIhYAAABQ4wTFIUwVkZubqzFjxig1NVX9+/c3Ow4AAABQLVTLApGT\nk6ORI0fKYrHozTff9Pl599xzj+65554qTAYAAAAEt6A4hOlyZGVl6b777pPNZtPMmTMVFxfn83Mz\nMzOrMBkAAAAQ/KrVHoiioiKNGjVKISEhmjVrlmrXrm12JAAAAKBaCfo9EDk5OYqKipLD4dDbb7+t\njIwMffDBB3K73crJyZEkhYWFKTIy0uSkAAAAQPALuj0QFovF4+s+ffqcu8/DokWLVFxcrOHDh6tv\n377nPl555RUzogIAAADVjsUwDMPsEIHi7NWaFi9ebHISAAAAIDAF3R4IAAAAAOahQAAAAADwGQUC\nAAAAgM8oEAAAAAB8RoEAAAAA4DMKBAAAAACfUSAAAAAA+IwCAQAAAMBnFAgAAAAAPqNAAAAAAPAZ\nBQIAAACAzygQAAAAAHxGgQAAAADgMwoEAAAAAJ9RIAAAAAD4jAIBAAAAwGcUCAAAAAA+o0AAAAAA\n8BkFAgAAAIDPKBAAAAAAfEaBAAAAAOAzCgQAAAAAn1EgAAAAAPiMAgEAAADAZxQIAAAAAD6jQAAA\nAADwGQUCAAAAgM8oEAAAAAB8RoEAAAAA4DMKBAAAAACfUSAAAAAA+IwCAQAAAMBnFAgAAAAAPqNA\nAAAAAPAZBQIAAACAzygQAAAAAHxGgQAAAADgMwoEAAAAAJ9RIAAAAAD4jAIBAAAAwGcUCAAAAAA+\no0AAAAAA8BkFAgAAAIDPKBAAAAAAfEaBAAAAAOAzCgQAAAAAn1EgAAAAAPiMAgEAAADAZxQIAAAA\nAD6jQAAAAADwGQUCAAAAgM8oEAAAAAB8RoEAAAAA4DMKBAAAAACfBVWBcDqduuWWW7R69eoLbrNt\n2zYNHz5cnTt31h133KGtW7f6MSEAAABQvQVNgXA6nXr66aeVnp5+wW2Kioo0evRode/eXZ9//rk6\nd+6sRx55RMXFxX5MCgAAAFRfQVEgdu/ereHDhysjI+Oi282fP1/h4eH64x//qObNm+v5559XRESE\nvvnmGz8lBQAAAKq3oCgQq1atUq9evfTJJ5/IMIwLbrdp0yalpqZ6zLp27ar169dXdUQAAACgRgiK\nAnHXXXdp7NixcjgcF93u6NGjSkxM9JglJCQoKyurKuMBQWfNmjW6/fbb1bdvX91+++1as2ZNubPL\neb4vr9+vXz/Vrl1btWvXVv/+/S/6Pa70+12pirymGTkrI/flvIY/foZAVNGfOxB+vQIhQ3nOzzVz\n5syAzAmgHEaQadOmjbFq1apy10aOHGn8/e9/95i9+eabxgMPPODTa/fr18/o16/fFWcEAtnq1auN\nxo0bG5LOfSQlJRl169b1mDVu3NhYvXq1T8//5bYXev2EhASPmSSjbt265X6PK/1+F8p+Jb9Gl3pN\nM3JWRu7LeQ1//AyBqKI/dyD8egVCBl9z2Wy2gMuJ6svlchulZS6j2FlmFBaXGvmFTuNkQYlx/FSx\nkXuiyMg5XmhkHSswMnPyjUNHTxkHjpw09mWeMPYcOm6kH8wzdh7INbbvO2Zs3ZNjbE7PNjbuOmps\n3HXUyC90mv2j+YW9EruI6RwOh5xOp8fM6XQqLCzMpERA4Jk0aZIOHDjgMStvL92BAwc0efJkzZ49\n+5LP/+W2vr6+JB05cqTc73Gl3+9C2X1Vkdc0I2dl5L6c1zAMo8p/hkBU0V9Xf/yeX0ogZChPeblc\nLpfH14GQE97cbkOlLrfKytwqLXOrzHX64+znpWfnZW6Vnjf/5aPL5ZbLbajMZcjldsvlMuRyG7+Y\nu+U+83h6fnq7Mpdx3vwC6+de0/OxzH16varYbRa98GBPdU1OvPTGQaxaFYikpCRlZ2d7zHJyclSn\nTh2TEgGB53IO6Stv2ws9/+z8cg8ZvNT2Ff1+V3LoYkVe04ycl5vhSl/DuMA5aNX9MNGK/rr64/f8\nUgIhw5V8f7NzBiqX21CJs0wlpS45S91ylrrOfJz+vKTMpdJS95n1Mx9l/92upNSl0rJfrrs9trvQ\nm/7SMneVvvmuDspchuYs2UWBCCadOnXS1KlTPWbr16/Xo48+alIiIPAkJSVd0bYXev7Z+eW8vi/b\nV/T7XW6Oy/meFXlOVeS83AxX+hoXKhCV+TMEoor+uvrj9/xSAiHDlXx/s3NeKcMwVFrmVkFxqYqK\ny1RYXKYiZ5lKnC4VlZSpxFmmohKXip1lKnaeeSw57/HsWsl/H51lbrN/NFxEXNTFz9mtDizGhf5G\nCFDJycmaOXOmunfvLun0HoaoqCg5HA7l5+frhhtu0M0336zf/OY3+vjjj7Vw4UJ9++23Ph3G1L9/\nf0nS4sWLq/RnAMy0Zs0a3XbbbR6HDyQlJclisejIkSPnZo0bN9Znn32mbt26XfL5v9z2Qq9fVlam\nY8eOebxW3bp19fXXX3t9jyv9fhfK7quKvKYZOSsj9+W8hqQq/xkCUUV/Xf3xe34pgZDB11w2m83j\nMKZAyFlS6lJ+oVOnCktVUFSqguJSFRaXqfC8x7MFobz1MldQvc3CJVgtktVqkdVikc12+tFqPf1h\ns1rUtlmCRg/toLjo6n34fNAViJSUFH3wwQfnCkRycrImTZqkoUOHSpI2b96s8ePHa8+ePWrTpo1e\neuklJScn+/TaFAjUFGvWrNHkyZOVlZWlpKQkjR07VpK8Zhd7s3yxbS/0+mPHjtXGjRslnd5jOHny\nZJ/eHFTk+13pm46KvKYZOSsj9+W8hj9+hkBU0Z87EH69AiGDL7kGDx6suXPnVknO4pIynShw6lSh\n81whOPt4enbmsajUY51/6b88VqtFIXar7DZruY8hNqvsZx5tNots1rOPFtltVlnPPNqslnPrdpvl\nvPl/18/NzryO3WqV1WaR/bztPNbPvNH/5Zv+X5aA8wuBx8wiWSwWs3+ZA0LQFYiqRIEAACDwlZS6\ndCK/RCfznTqeX6IT+SU6ke88/Vjwi8/zS3Q83ylnqevSLxrkbFaLQkNscoTYFBpiVYj9v5+HhtjO\nfFjPbRNit55Z/8Wa/cxzbbbTb/R/+ab/IsXAfuZzm5U31zVFtToHAgAABK/ikjLlniw+76NEuSc8\nZ0UlZWZHrRC7zapwh02OUPt/H0PtcoTaFO6wKyzUprAzj45Q25m109uGhdoVdvYx9PSjI/TMm3+7\nVTZbUNzaC9UEBQIAAFQpwzBUUFym7LxCZR8vUnZekbLzCpVzvFh5p4p17ERgFwOLRarlsKtWeIgi\nwkIU7rArIjzk3Oz0o121HCGKCLcr/MxjLUeIxzzEbjP7RwEqBQUCAABcEbfbUO7JYh05VvDfgnC8\nyKMwBEI5sFktiqoVqshaIeU+RoWHKLJW6LnZ2XmtsBAOzwF+gQIBAAAuyVnqUlZuoY4cK1DmsQId\nOVaozJwCZeUWKOtYoSknHNusFsVEhiom0nH6I8KhmKjQ04+RDsX+ci0yVOEOOyfBApWAAgEAACRJ\nLpdbWXmFyjiar4ysfGUcPXW6LOQU6NjJYvnrsivhDrvio8OUEBOmuKgwxceEKT7aofjoMMVHhyk2\nyqHYSIciwkMoBIAJKBAAANQwhcWlOpSdf7ooHD1dFDKO5utwdoHKXFW3J8FikeKiwlQnLlx1YsNV\nOzZcCTFh54pBfHSY4qLDFO7g7QkQyPg/FACAaqq0zKWDWfnal3lC+zJPad/hEzqQdUrHThRXyfcL\nDbEpKT5cdWJrnSsJpx9Pf50QE64QO1cLAoIdBQIAgCBnGIayjxdpX+ZJ7c88qX2HT2pv5kkdys6X\n2125xx3FRjqUlFBL9RIiVDchQnUTaqluQoTq1Y5QXJSDQ4qAGoACAQBAEDEMQ9l5RUrPOH764+Bx\npWec0KlCZ6V9j+iIUDVMjFTDxCg1qBNxriAkxddSrbCQSvs+AIITBQIAgACWc7xIuw7mKT3jxJmy\ncFwnC668LFgtUlJ8hBokRqphYqQaJUWpYWKkGtSJVEykoxKSA6iuKBAAAASI0jK39hw6ru378pS2\nP1dp+3Ir5XyFpPhaalovWk3rRatJvWg1TopSvdoRCg3hxmYALh8FAgAAk+SdLFba/tzThWFfrtIz\njqv0Cu6nEBFmV9P6MWpSN0pN68eoWb1oNa4bxWFHACoVBQIAAD85dqJIm9NztCk9R1t2H1PmsYIK\nv1ZEeIhaNYxVi4YxatUoTi0axigpvhYnMQOochQIAACqSO7JYm1Oz9Hm3TnanJ6jwzkVKwwR4SFq\n2TBGLRvGqmWjWLVsGEtZAGAaCgQAAJUkv6hUG3dla+PObG1Kz9Gh7PwKvU6jpCilNI1XcpM4JTeN\nV4M6kbJaKQsAAgMFAgCACnK5De3OOK51O45qXdpR7TiQd9n3XQh32NSmcbzaNI1TStN4tWkcp8ha\noVWUGACuHAUCAIDLcPxUidamZWld2lGt35l92fdfiAizq32L2mrforY6tEhQ0/oxsrF3AUAQoUAA\nAHARhmEo42i+Vm49olVbjyhtf66My9jJEO6wq13zBHVsWVsdWtZWMwoDgCBHgQAA4Dwul1vb9uZq\n1bYjWrn1iDIv4+TnELtV7ZonqHOrOurQsrZaNIiRzWatwrQA4F8UCAAAdPombht3ZevnjYe1cmum\nThWW+vzchomR6tomUV2TE9WueYLCQvnrFUD1xZ9wAIAa62xpWLrxkFZsOaKCIt9KQ7jDps6tE0+X\nhjaJSoyvVcVJASBwUCAAADVKmcutDTuz9dOGQ1q5JVMFxWU+Pa92TJiualdXPdrVU4eWCQqx26o4\nKQAEJgoEAKDaMwxDuw4e1/drD+qnDYd0It+3Kye1aBijHm3r6qp2ddW8QQw3bgMAUSAAANXYkWMF\n+nFdhr5fe1CHsn07ETq5SZyu7tRAvTvWU2IchyYBwPkoEACAaqWwuFQ/bTisJWsOaNveXJ+ek9I0\nXld3qq/eHeqrTlx4FScEgOBGgQAABD3DMLTzQJ4WrtivnzYcUrHTdcnntGkcp2u6NFDvjvVVO5bS\nAAC+okAAAILWqUKnvl97UItW7Nf+I6cuuX3dhFq6rmsjXZ/aUPXrRPohIQBUPxQIAEBQMQxD2/bm\nasGyfVq2+bBKy9wX3T6qVoj6dG6g67s2UnLTOE6EBoArRIEAAASFklKXflyXoflL92rP4RMX3dZq\nkbql1NWAqxqpW0pdhdi5EzQAVBYKBAAgoB3NK9S/f96rRSv3X/Lu0EnxtTSwR2MN6N5YCTGc1wAA\nVYECAQAIOIZhaMvuY/p66R6t3JIpt3Hhbe02q3p1qKdBPRqrY8s6slo5RAkAqhIFAgAQMFwut37e\ndFiffZ+uPYcufphS3YRauvnqZro+tZFiIh1+SggAoEAAAExX7CzTd6sO6Isfd+tobuFFt+3aJlG/\n7tNMqclJ7G0AABNQIAAApjmRX6J//7xXXy/dq1OFzgtuF+6wqX+3xrq5TzM1TIzyY0IAwPkoEAAA\nv8vOK9Ln3+/SolUH5Cy98E3f6tWO0C19mqt/90aqFRbix4QAgAuhQAAA/CY7r0izl+zUtysPqMx1\n4fs3tGoUq9v6tVLP9vVk4zAlAAgoFAgAQJU7mleoOYt36dtV+1XmuvAllbqlJOnW61uqffMEbvgG\nAAGKAgEAqDJHcws1e8kufXeR4mCzWnRt14Yadl1LNa0X7eeEAIDLRYEAAFS6YyeK9Mm3Oy+6x8Fu\ns+pXPZvo1utbqU4cN30DgGBBgQAAVJr8olJ9/v0uffWfPRc8OTrEbtUNPZvo9n6tuFs0AAQhCgQA\n4Io5S12a//NezV68U6cKS8vdhuIAANUDBQIAUGEut6Hv1xzUhwvTlHO8qNxtQuxW/apXU912fUuK\nAwBUAxQIAMBlMwxDa7Zn6b3523TgyKlyt7FaLRp4VWPdNagNxQEAqhEKBADgshzMOqVpX23Ruh1H\nL7hN7471dO+NKdw1GgCqIQoEAMAn+YVOfbxoh+b9vFdud/lXVurQorZG3pyiNk3i/ZwOAOAvFAgA\nwEW5XG4tXLlfsxak6VShs9xtmtWP1sib26prm0RuAAcA1RwFAgBwQRt3ZWvaV1u0L/Nkuevx0Q6N\nvLmtruvaSFYrxQEAagIKBADAy7ETRZr61Rb9vPFwueshdquGXddSt/drpXAHf5UAQE3Cn/oAgHNc\nbkP//nmvZi7YrqKSsnK36d2xnh74dTvVTYjwczoAQCCgQAAAJEm7Dubp/+ZsVHrGiXLXm9aL1uih\nHdShZW0/JwMABBIKBADUcAVFpZq1YLvmL9sro5yLK0XVCtG9N7XVoB5NZOM8BwCo8SgQAFBDGYah\npRsPa9pXm5V7sqTcbQZ0b6z7f91WMZEOP6cDAAQqCgQA1EDHThTp/+Zs0qptR8pdb5QUqcdv66T2\nLThcCQDgiQIBADWIYRhasuagpn61RQVFpV7roXar7hzURkOvbakQu9WEhACAQEeBAIAaIud4kd6a\nvUFr046Wu56anKhHb+3I1ZUAABcVFP+85HQ69dxzz6l79+7q27evZsyYccFtv/32W918883q0qWL\nRowYoW3btvkxKQAEHsMwtHDFfv32r0vKLQ+xkQ6Nva+bxo/qSXkAAFxSUOyBmDx5srZt26aZM2cq\nIyNDY8eOVYMGDTRo0CCP7dLT0/XMM8/o5ZdfVpcuXfTee+9p9OjRWrx4sRwOTgAEUPMczS3U32dv\n0Iad2eWuX9e1oR4e2kHREaF+TgYACFYBXyCKioo0Z84cTZ8+XcnJyUpOTtaoUaM0a9YsrwKxdOlS\ntWrVSoMHD5YkPf300/rwww+Vnp6udu3amREfAExhGIYWrz6od77cpKISl9d6fLRDj9/WST3a1zMh\nHQAgmAV8gUhLS5PL5VLnzp3PzVJTUzVlyhSvbWNjY5Wenq5169apS5cu+uyzzxQVFaXGjRv7MzIA\nmOpkgVP/mLNByzZllrver1sjPTykvSJrsdcBAHD5Ar5AZGdnKzY2Vnb7f6MmJCSopKREeXl5iouL\nOze/6aabtGTJEt19992y2WyyWq165513FBUVZUZ0APC79TuO6n//ta7c+zokxITpiTs6q1tKkgnJ\nAADVRcAXiKKiIoWGev4r2dmvnU6nx/z48ePKycnR+PHj1alTJ3388cd69tln9cUXXyg+Pt5vmQHA\n35ylLr3/722a+5895a4P6N5Yo4a0V0R4iJ+TAQCqm4C/CpPD4fAqCme/Dg8P95i/9tpratOmje66\n6y61bdtWEydOVHh4uD7//HO/5QUAf9t7+ISe+t8fyy0PUbVC9dz93fXknV0oDwCAShHweyCSkpJ0\n/Phxud1uWa2n+05OTo7CwsIUHR3tse3WrVt13333nfvaYrEoOTlZhw8f9mtmAPAHt9vQ3J926/35\n21Xmcnutd22TqCfv7KL46DAT0gEAqquALxApKSmy2+3asGGDunbtKklas2aN2rdv77VtYmKi0tPT\nPWZ79+5Vx44d/ZIVAPzlZIFT//uvdVq9LctrLdRu1f2/bqdf92kmi8ViQjoAQHUW8AUiLCxMQ4YM\n0fjx4/Xqq68qKytLM2bM0KRJkySd3hsRFRUlh8OhO+64Q88995zat2+vLl266NNPP1VmZqaGDh1q\n8k8BAJVn+95c/b9Za5RzvMhrrVn9aP1hRKqa1I0u55kAAFy5gC8QkjRu3Di99NJLGjlypKKiovTk\nk09qwIABkqQ+ffpo0qRJGjp0qG666SYVFRVpypQpysrKUkpKij744ANOoAZQLbjdhj7/IV0zF2yX\n2214rFks0rBrW+qeG5MVYreZlBAAUBNYDMMwLr1ZzdC/f39J0uLFi01OAgCeTuSX6I2P12lt2lGv\ntdhIh/4woqs6t040IRkAoKYJij0QAFCTbd1zTH+dtUbHThR7rXVoUVvP3JPKidIAAL+hQABAgDIM\nQ599f+FDlu4c2Ea/GdhGNisnSgMA/IcCAQABqLC4VP/7r/VavjnTay02yqFn7k5Vp9Z1TEgGAKjp\nKBAAEGAyjp7SKzNWKeNovtdax5a19cyIVMVxyBIAwCQUCAAIIMs3Z+qNj9epqKTMY26xSHcNbKPh\nHLIEADAZBQIAAoDLbeijhWn69LudXmuR4SH64z3d1DWZqywBAMxHgQAAk50qdOq1D9dqXTmXaG1W\nP1rP3X+V6iZEmJAMAABvFAgAMNHewyf06nurdORYodfadV0b6rd3dFJYKH9UAwACB38rAYBJVm7J\n1GsfrlWx0+Uxt1ktemhwe/26TzNZLJzvAAAILBQIAPAzwzD0+ffpev/f22R43t5BsVEOjb23m9q3\nqG1OOAAALoECAQB+VFrm0luzN2rJmoNea22axGncyO5KiAk3IRkAAL6hQACAnxw/VaJX31ul7fty\nvdb6dWukJ+7opBC7zYRkAAD4jgIBAH6w9/AJvfzuSmXnFXnMLRZp5E1tdev1LTnfAQAQFCgQAFDF\nLnSydFioTc+MSFWP9vVMSgYAwOWjQABAFTEMQ1/9Z4/e/XqL18nSdeLC9cKDPdSsfow54QAAqCAK\nBABUAZfb0LSvNmve0r1eaylN4zXu/u6KiwozIRkAAFeGAgEAlazYWabXP1yrFVuOeK1xsjQAINhR\nIACgEh0/VaI/v7tSOw7kea3de2OK7ujfipOlAQBBjQIBAJXkUHa+JkxdriPHCj3mdptFT97ZVdd1\nbWhSMgAAKg8FAgAqwba9x/Tnd1fqVGGpxzwiPETP33+VOrTkztIAgOqBAgEAV+jnjYf1+kdrVVrm\n9pgnxoVr/Kiealw32qRkAABUPgoEAFyBeUv36J0vN3tdprVFwxiNf6in4qK50hIAoHqhQABABRiG\noY8W7tC/vt3htdYtJUn/373dFO7gj1gAQPXD324AcJlcbkNTPt+kBcv3ea3d2KupHhnWQTab1e+5\nAADwBwoEAFyG0jKXXv9wnX7edNhrbcSvkvWbAa25TCsAoFqjQACAjwqLS/XKjFXalJ7jMbdYpMdu\n66QbezU1JxgAAH5EgQAAHxw/VaKXpi1XesYJj7ndZtUzI1J1daf6JiUDAMC/KBAAcAlZuYV6ccoy\nHc4p8JiHO2x6/oEe6tSqjknJAADwPwoEAFzEwaxT+tPby5R7sthjHhMZqgmjeqllo1iTkgEAYA4K\nBABcwJ5DJ/TClGU6WeD0mCfGhWviI73VoE6kSckAADAPBQIAypG2P1cTpq5QQVGpx7xJ3Si9NLqX\nEmLCTUoGAIC5KBAAcJ5N6dl6efpKFTtdHvPWjWM14eFeiqoValIyAADMR4EAgF9Ysz1Lf3lvlZxl\nbo95+xYJeuHBHqoVFmJSMgAAAgMFAgDO+HnjYb324RqVuQyPeWpyop4d2V1hofyRCQAAfxsCgKTF\nqw/ob5+sl9uzO6h3x3p6ZkQ3hdit5gQDACDAUCAA1HgLlu3V/322yWver1sj/W54Z9lslAcAAM6i\nQACo0eYt3aMpX2z2mt/Yu6keHdZRVqvFhFQAAAQuCgSAGmvuf3Zr6ldbvOa3XtdS9/+6rSwWygMA\nAOejQACokb78cbemz/UuD3cPaqM7B7WhPAAAcAEUCAA1zuffp2vGvK1e83tuTNZvBrQxIREAAMGD\nAgGgRvlsyS69N3+b1/y+m1J0R//WJiQCACC4UCAA1BizF+/UB//e7jW//+a2uq1fKxMSAQAQfCgQ\nAGqET77doVnfpHnNH7ylnYZd19KERAAABCcKBIBq7+OFafpo0Q6v+UOD22votS1MSAQAQPCiQACo\ntgzD0EcLd+hf33qXh4eHttfgvpQHAAAuFwUCQLX1r0Xll4dHh3XQzX2am5AIAIDgR4EAUC3NXryz\n3MOWHr+to27s3cyERAAAVA8UCADVzpc/ppd7taXf3t5Jv+rV1P+BAACoRqxmBwCAyjRv6R5Nn+t9\nk7jHKQ8AAFQKCgSAauOb5fs05YvNXvPRQzvoRsoDAACVggIBoFr4btV+/WPORq/5g7e00y19OWEa\nAIDKQoEAEPR+WHtQf/t0g9f83htTuEkcAACVjAIBIKgt3XhIb3y8TobhOb9rUBsNH9DanFAAAFRj\nFAgAQWv55ky9Nmut3OeVh9v7tdJdg9qYEwoAgGouKAqE0+nUc889p+7du6tv376aMWPGBbfdsWOH\n7r77bnXq1EmDBw/WypUr/ZgUgL+s3nZE/2/marnOaw9Drmmh+25KkcViMSkZAADVW1AUiMmTJ2vb\ntm2aOXOmxo8fr7feekuLFi3y2i4/P18PPfSQWrVqpXnz5mngwIF64oknlJuba0JqAFVl485s/eX9\n1SpzeZaHm69upocGt6M8AABQhQK+QBQVFWnOnDn605/+pOTkZA0YMECjRo3SrFmzvLb9/PPPFRER\noZdeekmNGjXSmDFj1LRpU23ZssWE5ACqQtq+XP15xkqVlrk95jf0bKLRQztQHgAAqGIBfyfqtLQ0\nuVwude61APh5AAAgAElEQVTc+dwsNTVVU6ZM8dp29erV6tevn8ds9uzZVZ4RgH/sPXxCE6atULHT\n5THv162RHr+tk6xWygMAAFUt4PdAZGdnKzY2Vnb7f7tOQkKCSkpKlJeX57HtwYMHFRcXpxdffFF9\n+vTRnXfeqXXr1vk7MoAqkHH0lF6cslwFRaUe86s71dfvftOF8gAAgJ8EfIEoKipSaGiox+zs106n\n02NeWFioadOmKTExUdOmTVO3bt300EMPKSsry295AVS+o7mFeuHtZTqeX+Ix75aSpD/cnSob5QEA\nAL8J+ALhcDi8isLZr8PDwz3mNptNKSkpeuKJJ5ScnKxnnnlGTZs21VdffeW3vAAqV97JYv1pyjLl\nnCj2mLdvkaBnR3ZXiD3g/xgDAKBaCfi/eZOSknT8+HG53f89YTInJ0dhYWGKjo722LZOnTpq3ry5\nx6xp06bKzMz0S1YAletUoVMvTFmmzJwCj3mrRrF64cEecoTYTEoGAEDNFfAFIiUlRXa7XRs2bDg3\nW7Nmjdq3b++1befOnZWWluYx27Nnjxo0aFDlOQFUrsLiUo1/Z7n2HznlMW9SN0oTHu6lWmEhJiUD\nAKBmC/gCERYWpiFDhmj8+PHavHmzvvvuO82YMUMjR46UdHpvREnJ6eOi77zzTu3YsUNvvfWWDhw4\noDfffFMZGRkaPHiwmT8CgMtUUurSy++u1K6Dxz3m9WpHaOIjvRUdEXqBZwIAgKoW8AVCksaNG6f2\n7dtr5MiRevnll/Xkk09qwIABkqQ+ffpowYIFkqT69etr+vTpWrJkiW655Rb9+OOPmjp1qhITE82M\nD+AylJa5Nen91dqy+5jHvHZMmP78SG/FR4eZlAwAAEiSxTAM49Kb1Qz9+/eXJC1evNjkJEDN5HIb\nem3WGi3deNhjHhMZqkm/7aOGiVEmJQMAAGcFxR4IANWf223oH7M3eJWHiPAQvfxIb8oDAAABggIB\nwHSGYejdr7fq21UHPOZhoTZNGNVTzerHmJQMAACcjwIBwHSzF+/SV//Z7TELsVv1pwd6KLlpvEmp\nAABAeSgQAEy1YPk+zVyw3WNmtVr07H3d1al1HXNCAQCAC6qUAlFaWqrNmzeroKDg0hsDwBk/bzys\nf3620Wv+1J1ddFW7uiYkAgAAl1KhApGZmakHH3xQmzZtUnFxsYYNG6Y77rhD/fr10/bt2y/9AgBq\nvA07j+q1D9fq/OvAjR7aQdelNjInFAAAuKQKFYi//OUvOnXqlOLj47VgwQIdPnxYH330kQYOHKi/\n/vWvlZ0RQDWz80CeXpmxSmUut8f8NwNb65a+zU1KBQAAfGGvyJNWrFih999/Xw0bNtRrr72mvn37\nqmvXroqLi9Ott95a2RkBVCMHs05pwtQVKna6POY39m6qETckm5QKAAD4qkJ7IEpLSxUTEyPDMLR8\n+XL17t1bkuR2u2W3V6iTAKgBsvOK9OKUZTpV6PSY9+3cQI8M6yiLxWJSMgAA4KsKvdtv27at5syZ\nozp16ujkyZO69tpr5XQ6NXXqVCUn8y+IALydyC/Ri+8sU86JYo9559Z19NRdXWWzUh4AAAgGFSoQ\nY8eO1aOPPqq8vDw9/PDDqlu3riZMmKDFixdr2rRplZ0RQJArLC7VS9NWKONovse8deNYPXf/VQqx\nc0VpAACChcUwzr8GyqUdOHBADRs2VH5+vqKjoyVJe/fuVVxcnGJjYys9pL/0799fkrR48WKTkwDV\nR2mZSxOnrdSGXdke80ZJkZr0276Kjgg1KRkAAKiICv2z34gRI7Rly5Zz5UGSmjVrFtTlAUDlc7kN\nvf7ROq/yUCcuXBNH96Y8AAAQhCpUIEJCQjhZGsBFGYahtz/fpJ83HvaYR0eEauLoXqodG25SMgAA\ncCUq1AKGDRumUaNGaciQIWrSpInCwsI81ocOHVop4QAErw+/SdM3y/d5zMIdNk14uKcaJkaZkgkA\nAFy5Cp0DcbErLVkslqC9GzXnQACVY+5/dmvqV1s8ZnabVRNG9VSn1nVMSgUAACpDhfZApKWlVXYO\nANXE92sPepUHq0V65p5UygMAANXAFV078fDhw/rpp59UXFysY8eOVVYmAEFq9bYj+t9/rfeaP357\nJ13dsb4JiQAAQGWr0B4Ip9OpsWPHasGCBbJarVq4cKEmT56s/Px8vfXWW4qMjKzsnAAC3La9xzTp\ngzVyuz2PirzvphTd0LOpOaEAAEClq9AeiH/+859KS0vT+++/L4fDIUm69957deDAAb322muVGhBA\n4Nt7+IQmTl8pZ6nLYz7kmha6vV8rk1IBAICqUKECMX/+fL3wwgvq0aPHuVmPHj30yiuvcAIyUMMc\nOVag8e8sV0FRqcf8+tSGevCWdrJYLCYlAwAAVaFCBSIrK0uNGzf2mterV08nTpy44lAAgkPeqWK9\n+M5y5Z0q8Zh3b5uk3/2mi6xWygMAANVNhQpEixYttHz5cq/5/Pnz1bJlyysOBSDwFRaXasLUFcrM\nKfCYt20Wr7H3dZfddkXXaAAAAAGqQidRjxkzRk899ZTS09Plcrn0xRdfaO/evVq4cKHeeOONys4I\nIMA4S116ZcYq7Tnkucexab1ovfBQTzlCbCYlAwAAVa1C/0R4/fXX629/+5u2bNkim82m6dOn6+DB\ng3rjjTd0ww03VHZGAAHE5Tb0+kdrtSk9x2OeGF9LL43upcjwEJOSAQAAf6jQHghJuuaaa3TNNddU\nZhYAAc4wDL39+SYt25TpMY+NdOjl0b0UHx1mUjIAAOAvPheIL7/80ucXHTp0aIXCAAhsHy3coW+W\n7/OYhTvsGv9wT9Wvw/1fAACoCXwuEM8++6zH1xaLRYZhKDw8XHa7XadOnZLNZlNcXBwFAqiG5i3d\no399u8NjZrdZ9fwDV6llw1iTUgEAAH/zuUCkpaWd+3zevHmaPn26/vKXvyg5OVmStG/fPj377LO6\n+eabKz8lAFP9tP6Q3vlys8fMYpGeuSdVnVrVMSkVAAAwQ4VOon7ttdc0YcKEc+VBkpo2barnn39e\nU6ZMqbRwAMy3fsdR/c/Ha2UYnvPHbuukqzvWNycUAAAwTYUKxMmTJ+VwOLzmbrdbxcXFVxwKQGDY\neSBPr763SmUuz/Yw4lfJurFXU3NCAQAAU1WoQPTo0UMTJ05URkbGudnu3bv10ksv6brrrqusbABM\nlHH0lCZMXaFip8tj/uurm+k3A1qblAoAAJitQpdxnTBhgh566CENHDhQ0dHRkk7vlejYsaNeeOGF\nSg0IwP+OnSjSi+8s16lCp8f8ms4N9PDQDrJYLCYlAwAAZqtQgUhKStJXX32lZcuWadeuXZKklJQU\n9ezZkzcWQJA7VejUi+8sV3Zekce8c+s6+v1dXWW18v84AAA1WYVvJGez2dS3b1/17du3MvMAMFGx\ns0wvT1+pA0dOecxbNYrVuJHdFWKv0FGPAACgGvG5QKSkpGjp0qVKSEhQcnLyRfc0bN++vVLCAfCf\nMpdbkz9Yo+37cj3mDepEavyonqoVFmJSMgAAEEh8LhCvvvqqoqKizn3OoUpA9eF2G/r7pxu0ZnuW\nxzwhJkwTR/dSTKT3VdcAAEDNZDGM86/uXnP1799fkrR48WKTkwD+9e7XW/XFD+kes8jwEE16oo+a\n1I02KRUAAAhEFT4HYsGCBXr//fe1c+dO2Ww2tW3bVg8//LD69OlTmfkAVLHPv9/lVR5CQ2x68aGe\nlAcAAOClQmdEzpkzR3/4wx9Uv359PfXUU/rtb3+rmJgYPfLII/ruu+8qOyOAKvLdqgOaMW+bx8xq\ntWjcyO5KaRZvUioAABDIKnQI06BBg3T33Xfr/vvv95hPmzZNc+fO1dy5cysrn19xCBNqklVbj+iV\n91bJ7fb8I+Cpu7qqX7dGJqUCAACBrkJ7ILKyssq94/TAgQO1f//+K80EoIpt3XNMkz9Y7VUeHhrc\njvIAAAAuqkIFolu3bvr3v//tNV+6dKlSU1OvOBSAqrMv86Renr5CzjK3x/y261tq6LUtTUoFAACC\nRYVOou7WrZv++c9/asuWLbrqqqsUEhKizZs3a968ebr11lv11ltvndv2iSeeqLSwAK5MVm6hxr+z\nTAXFZR7zgVc11sib25qUCgAABJMKnQPRr18/317cYgmq8wk4BwLVWd6pYo19a6kycwo85j3a1dW4\nkd1ls3GXaQAAcGkV2gOxZMmSc5/n5uZq9erVql27NocvAQEqv6hU499Z7lUe2jaL1x/v7UZ5AAAA\nPrusdw3/+Mc/1KNHj3MnSq9fv16DBg3S73//e91zzz164IEHVFxcXCVBAVRMsbNME6et0N7DJz3m\nTetF64WHesoRYjMpGQAACEY+F4hPPvlEb7/9toYPH66EhARJ0rhx4xQWFqavv/5aP/zwgwoKCvTO\nO+9UWVgAl6e0zK1J76/W9n25HvN6CRF6aXQvRYaHmJQMAAAEK58LxOzZs/Xss8/qD3/4gyIjI7V5\n82bt27dP9957r1q2bKmkpCQ99thjmj9/flXmBeAjl9vQGx+v09q0ox7z+GiHJj7SS/HRYSYlAwAA\nwcznArF7925dffXV575esWKFLBaLrr322nOzli1b6vDhw5WbEMBlMwxDUz7fpJ82HPKYR9UK0cRH\neqtuQoRJyQAAQLC7rHMgLBbLuc/XrFmjmJgYJScnn5sVFBQoPDy88tIBqJCZC7ZrwfJ9HrOwUJvG\nj+qpJnWjTckEAACqB58LROvWrbVu3TpJ0smTJ7Vy5UqPPRKStGDBArVu3bpyEwK4LJ9/n67Zi3d5\nzOw2q55/4Cq1aRJvUioAAFBd+HwZ1xEjRmj8+PHavn271q9fL6fTqZEjR0qSsrKy9PXXX2v69Ol6\n5ZVXqiwsgItbtHK/Zszb6jGzWqQ/3pOqzq0TTUoFAACqE58LxODBg+V0OvXxxx/LarXqjTfeUMeO\nHSVJU6ZM0aeffqqHH35YQ4YMqbKwAC7s502H9Y/ZG7zmT9zRWb071jchEQAAqI4qdCfq82VlZSk0\nNFRxcXGVkck03IkawWr9jqOaOH2lylxuj/mDt7TTsOtampQKAABUR5Vy+9mkpKQqLQ9Op1PPPfec\nunfvrr59+2rGjBmXfE5GRoa6dOmi1atXV1kuIBCk7c/Vq++t8ioPd/RvRXkAAACVzudDmMw0efJk\nbdu2TTNnzlRGRobGjh2rBg0aaNCgQRd8zoQJE7grNqq9fZkn9dLUFSp2ujzmN/ZqqntvTDEpFQAA\nqM4qZQ9EVSoqKtKcOXP0pz/9ScnJyRowYIBGjRqlWbNmXfA5c+fOVWFhoR9TAv53KDtfL7y9TPlF\npR7zvp0b6JFbO3pcdhkAAKCyBHyBSEtLk8vlUufOnc/NUlNTtWnTpnK3z8vL0+uvv66JEyeqEk7v\nAAJSVm6h/vTPn3U8v8Rj3jU5UU/d1VU2K+UBAABUjYAvENnZ2YqNjZXd/t+jrRISElRSUqK8vDyv\n7SdNmqRhw4apZUuO/Ub1dOxEkV54e5lyTngeoteueYLG3dddIfaA/98aAAAEsYB/p1FUVKTQ0FCP\n2dmvnU6nx3zZsmVav369Hn/8cb/lA/zpRH6JXpiyXJnHCjzmLRvF6sWHeijMERSnNQEAgCAW8AXC\n4XB4FYWzX4eHh5+blZSUaMKECRo/frxX4QCqg/yiUr34znIdzDrlMW9aL1oTR/dSrbAQk5IBAICa\nJOD/uTIpKUnHjx+X2+2W1Xq67+Tk5CgsLEzR0dHnttu0aZMOHjyoMWPGeJz78PDDD2vo0KGaMGGC\nv6MDlaaopEwvTV2uPYdOeMwb1InQxEd6KaoWpRkAAPhHwBeIlJQU2e12bdiwQV27dpUkrVmzRu3b\nt/fYrlOnTlq0aJHHbODAgXrllVfUq1cvv+UFKltJqUt/fnel0vZ7nvOTGBeulx+5WnFRYSYlAwAA\nNVHAF4iwsDANGTJE48eP16uvvqqsrCzNmDFDkyZNknR6b0RUVJQcDocaNWrk9fzExETFx8f7OzZQ\nKUrL3Jr0/mptSs/xmMdHO/TnR69WnbjwCzwTAACgagT8ORCSNG7cOLVv314jR47Uyy+/rCeffFID\nBgyQJPXp00cLFiwo93lcBx/BzOVy6/WP1mrN9iyPeXREqF5+pLfq1Y4wKRkAAKjJLAY3Szinf//+\nkqTFixebnAQ1ndtt6M1P1mvJmoMe84gwu1557Gq1aBhrUjIAAFDTBcUeCKAmcbsN/WPORq/yEBZq\n04SHe1EeAACAqSgQQAAxDEP//HyTFq3c7zEPtVv1wkM9lNyU83kAAIC5KBBAgDAMQ1O+2Kxvlu/z\nmNttFo27/yp1bFnHlFwAAAC/RIEAAoBhGJr21RbN/3mvx9xmtWjcyKvULSXJpGQAAACeKBCAyQzD\n0Ltfb9Xcn/Z4zK1Wi8be101XtatrUjIAAABvFAjARIZh6IN/b9eXP+72mFutFv3xnlT16lDfpGQA\nAADlo0AAJvpwYZrmLNnlMbNapKfv6qo+nRqYlAoAAODCKBCAST5etEOffLvTY2axSL+/q6uu7drQ\npFQAAAAXR4EATPDpdzv10cI0j5nFIv1ueBddn9rIpFQAAACXRoEA/GzOkl2auWC71/y3t3fWgKsa\nm5AIAADAd3azAwA1ySff7tCsb9K85o/d1lE39GxiQiIAAIDLQ4EA/MAwDH28aIc+XrTDa2300A66\nqXczE1IBAABcPgoEUMUMw9CH36Tpk+92eq2NGtJet/RtbkIqAACAiqFAAFXIMAy9P3+bPvs+3Wvt\nkWEd9Os+lAcAABBcKBBAFTl7h+nzbxInSY/f1lE3ctgSAAAIQhQIoAoYhqFpX23R3J/2eMwtltNX\nW+KEaQAAEKwoEEAlMwxDU77YrPk/7/WYn73PA5dqBQAAwYwCAVQil9vQP2Zv0LerDnjMrWfuMM1N\n4gAAQLCjQACVpMzl1v98tE4/bTjkMbdapKfvTtW1XRualAwAAKDyUCCASuAsdWnyB2u0atsRj7nV\natEzI1LVt3MDk5IBAABULgoEcIWKSsr0yoyV2rgrx2Nut1k19r5u6tm+nknJAAAAKh8FArgC+UWl\nmjhthbbvy/WYO0Jtev7+q9SlTaJJyQAAAKoGBQKooBP5JXrxneXac+iEx7xWmF3jR/VU22YJJiUD\nAACoOhQIoAKOnSjSC1OW6WBWvsc8qlaoJo7upZaNYk1KBgAAULUoEMBlyswp0IvvLNORY4Ue8/ho\nhyY+0ltN6kablAwAAKDqUSCAy7A747gmTF2h4/klHvPEuHD9+dGrVa92hEnJAAAA/IMCAfhoc3qO\nXn53pYpKyjzmDepE6s+P9lbt2HCTkgEAAPgPBQLwwbJNh/XXWWtV5nJ7zJvXj9GE0T0VFxVmUjIA\nAAD/okAAl/DN8n3652cb5TY85x1a1NbzD1yliPAQU3IBAACYgQIBXIBhGPrku5368Js0r7VeHerp\nmRGpCg2xmZAMAADAPBQIoBxut6F3vtys+T/v9Vq7oWcTPXZbJ9msFhOSAQAAmIsCAZzHWerSGx+v\n09KNh73WfjOwtUbckCyLhfIAAABqJgoE8AsnC5x6ZcZKbdub6zG3WKTRQzvo132am5QMAAAgMFAg\ngDOOHCvQhKnLdSi7wGNut1n01F1ddU2XhiYlAwAACBwUCEDSjv25evndlTqR7/SYhztsGjfyKnVp\nk2hSMgAAgMBCgUCNt3xzpl77cK2cpS6PeUJMmMaP6qlm9WNMSgYAABB4KBCo0eb+tFvTvtoi47x7\nPDStF63xo3pyd2kAAIDzUCBQI7ncht79eovm/meP11rn1nU0bmR31QrjBnEAAADno0CgxiksLtX/\nfLROK7ce8Vob0L2xfntHJ9ltVhOSAQAABD4KBGqUo7mFevndldqXedJr7Z5fJWv4gNbc4wEAAOAi\nKBCoMbbtPaZX31vldaUlu82iMcO7qF+3RiYlAwAACB4UCNQI3606oH/M2aAyl+fZ0lG1QvTsyO7q\n2LKOSckAAACCCwUC1ZrLbei9eVv15Y+7vdYaJUXqTw/2UP3akSYkAwAACE4UCFRbhcWl+uustVqz\nPctrLTU5UX+8p5siwrnSEgAAwOWgQKBaOpyTr1dmrNKBI6e81oZc00IP3NJONisnSwMAAFwuCgSq\nndXbjuj1D9eqoLjMY263WfTYbZ00qEcTk5IBAAAEPwoEqg2329An3+7QR4t2eK1F1QrVc/d3V/sW\ntU1IBgAAUH1QIFAt5BeV6vUPyz/foXHdKL3wYA/VTYgwIRkAAED1QoFA0NuXeVKvzlilzGMFXmt9\nOzfQmOGdFe7gP3UAAIDKwLsqBLUf12Xo77M3qMTp8phbrRY98Ou2GnJNC+4sDQAAUIkoEAhKpWVu\nvTd/q+b+Z4/XWkxkqMbe210dWnK+AwAAQGWjQCDoZOUW6q8z12jHgTyvtdaNYzVu5FWqHRtuQjIA\nAIDqjwKBoLJs02H97dMNKigq9Vr7Va+mGj20vULsNhOSAQAA1AwUCASF0jKX3p27VfN+3uu1FmK3\n6tFbO3J/BwAAAD+wmh3AF06nU88995y6d++uvn37asaMGRfc9ocfftDQoUPVpUsXDRkyREuWLPFj\nUlSFwzn5+uPffyq3PNSvHaHXfncN5QEAAMBPgmIPxOTJk7Vt2zbNnDlTGRkZGjt2rBo0aKBBgwZ5\nbLdjxw6NGTNGzz77rK655hr95z//0e9+9zt99tlnatOmjUnpcSV+Wn9If5+9QUUlZV5r13ZpqMdv\n76haYSEmJAMAAKiZAr5AFBUVac6cOZo+fbqSk5OVnJysUaNGadasWV4FYt68eerVq5dGjBghSRox\nYoSWLFmiBQsWUCCCTFFJmabP3aKFK/Z7rYXarRo9rKMG9WjMJVoBAAD8LOALRFpamlwulzp37nxu\nlpqaqilTpnhtO2zYMJWWep9cm5+fX6UZUbl27M/V6x+tU2aO943hGiZGaux93dW0XrQJyQAAABDw\nBSI7O1uxsbGy2/8bNSEhQSUlJcrLy1NcXNy5efPmzT2eu2vXLq1YsUJ333233/Ki4lwutz79bqf+\n9d1Oud2G13q/bo302K0dFcZdpQEAAEwT8O/EioqKFBoa6jE7+7XT6bzg83JzczVmzBilpqaqf//+\nVZoRV+5wTr7+58N15d7bISzUpkdv7aj+3RubkAwAAAC/FPAFwuFweBWFs1+Hh5d/s7CcnBw98MAD\nslgsevPNN6s8IyrOMAwtWnlA077arGKny2u9TeM4PT2iq+rXjjQhHQAAAM4X8AUiKSlJx48fl9vt\nltV6+qqzOTk5CgsLU3S093HwWVlZuu+++2Sz2TRz5kyPQ5wQWI6fKtFbszdo5dYjXmtWq0V3Dmit\n4QNay2YLiqsNAwAA1AgBXyBSUlJkt9u1YcMGde3aVZK0Zs0atW/f3mvboqIijRo1SiEhIfrggw8U\nHx/v77jwgWEY+mnDIU35YrNOFngfhlavdoT+cHdXtWnC7x8AAECgCfgCERYWpiFDhmj8+PF69dVX\nlZWVpRkzZmjSpEmSTu+NiIqKksPh0Ntvv62MjAx98MEHcrvdysnJOfcakZEcAhMIck8W6//mbCx3\nr4Mk3dCziR4a3F7hnCgNAAAQkCyGYXhf7ibAFBcX66WXXtLChQsVFRWlUaNG6d5775UkJScna9Kk\nSRo6dKhuvPFG7du3z+v5Q4cO1V/+8pdLfp+zJ1svXry4UvPj9F6HJWsOaupXW1RQ5H2p3ZjIUI25\no7N6tK9nQjoAAAD4KigKhL9QIKpGdl6R/jFng9amHS13vXfHenr01o6KiwrzczIAAABcLo4TQZVx\nuw19u2q/ps/dqqKSMq/1mMhQPXZrJ13dqb4J6QAAAFARFAhUif2ZJ/V/n23Utr255a5f26WhHh7a\nXjGRDj8nAwAAwJWgQKBSFZeU6V/f7tCXP+6Wq5y7ScdHO/T4bZ041wEAACBIUSBQaVZuydSULzcr\nO6+o3PUB3RvrocHtFFkrtNx1AAAABD4KBK7Y0bxCvfPF5gtemrVuQi09dmsndU1O9HMyAAAAVDYK\nBK7IguX7NH3uFpU4XV5rdptFt13fSncMaC1HiM3/4QAAAFDpKBCosHU7jur/5mwsd61Di9p67LaO\napQU5edUAAAAqEoUCFTY3kMnvGYxkaF6aHB7Xde1oSwWiwmpAAAAUJUoEKiwqzvV10eLdshZ6pLF\nIv2qZ1Pdd1MKJ0kDAABUYxQIVFjdhAhNfW6ANuzMVpsmcWpQJ9LsSAAAAKhiFAhckfjoMPXr1sjs\nGAAAAPATq9kBAAAAAAQPCgQAAAAAn1EgAAAAAPiMAgEAAADAZxQIAAAAAD6jQAAAAADwGQUCAAAA\ngM8oEAAAAAB8RoEAAAAA4DMKBAAAAACfUSAAAAAA+IwCAQAAAMBnFAgAAAAAPqNAAAAAAPAZBQIA\nAACAzygQAAAAAHxGgQAAAADgMwoEAAAAAJ9RIAAAAAD4jAIBAPj/27v3qCjrPI7jn1EUrDSTi2Fq\nbe0m5gUQtNBFDa9oBW5GKKtloWZq6VlKq03TdQ2TsLbWLpYeAbcLallpq6Hbnl07G+I9FVPbs4hX\nRrQ0bgnP/tFhcpoBnjGYGeH9OodzmN/ze57nO98z+OPjPM8AAIBpBAgAAAAAphEgAAAAAJhGgAAA\nAABgGgECAAAAgGkECAAAAACmESAAAAAAmEaAAAAAAGAaAQIAAACAaQQIAAAAAKYRIAAAAACYRoAA\nAAAAYBoBAgAAAIBpBAgAAAAAphEgAAAAAJhGgAAAAABgGgECAAAAgGkECAAAAACmESAAAAAAmEaA\nAAAAAGAaAQIAAACAaQQIAAAAAKYRIAAAAACYRoAAAAAAYNoVESAqKir09NNPq3fv3oqOjtaKFStq\nnLt//34lJCQoLCxM9913n/bt2+fGSgEAAIDG7YoIEIsWLdL+/fuVmZmpuXPn6tVXX9WmTZsc5pWW\nlsltcTIAABoUSURBVGrSpEnq3bu31q5dq7CwME2ePFllZWUeqBoAAABofLw+QJSWlmr16tX64x//\nqJCQEA0ePFjJycnKyspymLt+/Xq1atVKTzzxhG6++WY988wzuvrqq/X3v//dA5UDAAAAjY/XB4j8\n/HxVVlYqLCzMNhYREaE9e/Y4zN2zZ48iIiLsxnr16qWdO3c2eJ0AAABAU+D1AaKoqEht27aVj4+P\nbczf31/l5eU6e/as3dzTp08rKCjIbszf31+nTp1yS62Au+Tl5SkmJkYBAQEKCAjQoEGDlJeX5+my\nnMrLy9Po0aMVHR2tmJgYDRo0SNHR0Ro9erSpmi/dv3ofZ2Nm9rsS/dL+AQBQ33zqnuJZpaWlatmy\npd1Y9eOKigq78bKyMqdzfz4PuJLl5eXprrvusgvGW7Zs0d13362PP/5YkZGRHqzOXl5enu69914V\nFBQ43b5t2zatWbOmxpqd7f/vf/9bFotFJ0+erPE4zvar61ze6Jf2DwCAhuD170D4+vo6BIDqx61a\ntTI118/Pr2GLBNwoNTXV6btqJ0+e1KJFizxQUc1SU1Nr/OVXkgoKCmqt2dn+p06dsgsPzo7jbL+6\nzuWNfmn/AABoCF4fINq3b69z586pqqrKNma1WuXn56c2bdo4zC0qKrIbs1qtCgwMdEutgDvUdkme\nt12uZ6ae+no+l86taT9v609dfmn/AABoCF4fILp27SofHx/t2rXLNpaXl6fu3bs7zA0NDXW4YXrn\nzp12N2ADV7r27dtf1jZPMFNPfT2fS+fWtJ+39acuv7R/AAA0BK8PEH5+foqLi9PcuXO1d+9e5eTk\naMWKFXrggQck/fgOQ3l5uSRp2LBhOn/+vBYuXKgjR45owYIFKikpUWxsrCefAlCvZs+e7fSXxuuv\nv16zZs3yQEU1mz17tjp37lzj9s6dO9das7P927dvr+uvv77W4zjbr65zeaNf2j8AABqCxTAMw9NF\n1KWsrEzz5s3Txo0b1bp1ayUnJ2vcuHGSpJCQEKWmpio+Pl6StHfvXs2dO1fffPONunTponnz5ikk\nJMTUeQYNGiRJ2rx5c8M8EaCe5OXladasWdq9e7ekH999W7RokVfeTJuXl6dFixbp1KlTatGihSwW\niyoqKtS+fXvNmjWrzpov3b96H0kOYz8/jrP9vLE/dfml/QMAoL5dEQHCXQgQAAAAQO28/hImAAAA\nAN6DAAEAAADANAIEAAAAANMIEAAAAABMI0AAAAAAMI0AAQAAAMA0AgQAAAAA0wgQAAAAAEwjQAAA\nAAAwjQABAAAAwDQCBAAAAADTCBAAAAAATCNAAAAAADCNAAEAAADANAIEAAAAANMIEAAAAABMI0AA\nAAAAMI0AAQAAAMA0AgQAAAAA0wgQAAAAAEwjQAAAAAAwjQABAAAAwDQCBAAAAADTCBAAAAAATCNA\nAAAAADCNAAEAAADANAIEAAAAANMIEAAAAABMI0AAAAAAMI0AAQAAAMA0AgQAAAAA0wgQAAAAAEwj\nQAAAAAAwjQABAAAAwDQCBAAAAADTCBAAAAAATCNAAAAAADCNAAEAAADANAIEAAAAANMIEAAAAABM\nI0AAAAAAMI0AAQAAAMA0AgQAAAAA0wgQAAAAAEwjQAAAAAAwjQABAAAAwDQCBAAAAADTCBAAAAAA\nTCNAAAAAADCNAAEAAADANAIEAAAAANMIEAAAAABMI0AAAAAAMI0AAQAAAMA0AgQAAAAA066IAJGW\nlqaoqCjdfvvtWrx4ca1zd+3apcTERIWHhys2NlbZ2dluqhIAAABo/Hw8XUBdli9frvXr12vp0qX6\n4YcflJKSooCAAE2YMMFhrtVq1aRJkzR27Fi98MIL+uqrr/TUU08pKChIAwYM8ED1AAAAQOPi9e9A\nZGZm6vHHH1d4eLj69OmjlJQUZWVlOZ2bk5OjwMBAzZgxQ507d9aIESMUFxenTz75xM1VAwAAAI2T\nV78Dcfr0aZ04cUKRkZG2sYiICB0/flxWq1UBAQF28/v376/bbrvN4Tjnz59v8FoBAACApsCr34Eo\nKiqSxWJRUFCQbSwgIECGYejkyZMO8zt06KCePXvaHp85c0YbNmxQ37593VIvAAAA0Nh5/B2I8vJy\nnTp1yum2kpISSVLLli1tY9XfV1RU1Hnc6dOnKygoSPfff7+pWoqKinTx4kUNGjTI1HwAAADA2wUH\nB9d4C8Dl8HiA2L17t8aPHy+LxeKwLSUlRdKPYeHnwaFVq1Y1HrOkpERTpkxRQUGB3nnnHfn6+pqq\npWXLljIMw9WnAAAAADQZFsOLf2M+ffq0BgwYoM2bN6tDhw6SpMLCQg0ZMkT/+te/HO6BkKQLFy4o\nOTlZhYWFWrlypW655RZ3lw0AAAA0Wl59D0RQUJCCg4O1fft221heXp6Cg4OdhgfDMDRt2jQdO3ZM\nWVlZhAcAAACgnnn8Eqa6JCYmKi0tTe3bt5dhGEpPT9fDDz9s215cXCw/Pz9dddVVys7OVm5url57\n7TVdc801slqtkqQWLVro2muv9dRTAAAAABoNr76ESZKqqqq0ePFirV27Vs2aNVNCQoJmzpxp2x4T\nE6Pf/e53mjZtmpKTk7V161aHY/Tu3VsZGRnuLBsAAABolLw+QAAAAADwHl59DwQAAAAA70KAAAAA\nAGAaAQIAAACAaQQIAAAAAKY1+QCRlpamqKgo3X777Vq8eHGtc3ft2qXExESFh4crNjZW2dnZbqrS\nsyoqKvT000+rd+/eio6O1ooVK2qcu3//fiUkJCgsLEz33Xef9u3b58ZKvYsrffv8888VHx+v8PBw\nxcXFacuWLW6s1Hu40rNqhYWFCg8P17Zt29xQoXdypW8HDx7U2LFjFRoaqnvuuUdffvmlGyv1Hq70\n7LPPPtPIkSMVHh6upKQk7d+/342Vep+Kigrdfffdtf7MsRY4MtM31gJ7ZnpWjbXgR2Z6Vi/rgNGE\nvf3228bAgQONHTt2GF9++aURHR1tLF++3OncoqIio3fv3saSJUuM//3vf8b69euNnj17Gp9//rmb\nq3a/+fPnG3FxccaBAweMzz77zOjVq5exceNGh3klJSVGv379jBdeeME4cuSIsWDBAqNfv35GaWmp\nB6r2PLN9y8/PN7p3725kZWUZBQUFRlZWltGtWzcjPz/fA1V7ltmeXerhhx82QkJCjNzcXDdV6X3M\n9u38+fNGv379jDlz5hgFBQXGX/7yFyMyMtI4c+aMB6r2LLM9O3TokNGzZ09j3bp1RkFBgTF//nyj\nX79+RllZmQeq9rzy8nJj6tSptf7MsRY4MtM31gJ7Znp2KdYCcz2rr3WgSQeIgQMHGh988IHt8bp1\n64yYmBinc9955x1jxIgRdmPPPvuskZKS0qA1elpJSYnRs2dPY9u2bbaxpUuXGuPGjXOYm52dbQwe\nPNhubOjQoXY9bipc6VtaWpoxceJEu7GHHnrIWLJkSYPX6U1c6Vm1devWGWPGjGnSi4YrfVu5cqUx\ndOhQu7HRo0cb//znPxu8Tm/iSs9WrFhh3HvvvbbHFy5cMLp06WJ89dVXbqnVmxw+fNiIi4sz4uLi\nav2ZYy2wZ7ZvrAU/MduzaqwF5ntWX+tAk72E6fTp0zpx4oQiIyNtYxERETp+/LjtL1hfqn///nr+\n+ecdxs+fP9+gdXpafn6+KisrFRYWZhuLiIjQnj17HObu2bNHERERdmO9evXSzp07G7xOb+NK30aN\nGqU//OEPDuMXLlxo0Bq9jSs9k6SzZ8/qxRdf1Pz582U04T9n40rftm3bppiYGLux7Oxs9e/fv8Hr\n9Cau9Kxt27Y6fPiwduzYIcMwtGbNGrVu3VqdO3d2Z8leITc3V1FRUXrvvfdq/ZljLbBntm+sBT8x\n2zOJtaCa2Z7V1zrgc1lVNgJFRUWyWCwKCgqyjQUEBMgwDJ08eVIBAQF28zt06KAOHTrYHp85c0Yb\nNmzQY4895raaPaGoqEht27aVj89PLxV/f3+Vl5fr7Nmzuu6662zjp0+f1q233mq3v7+/vw4fPuy2\ner2FK327+eab7fY9dOiQ/vOf/2js2LFuq9cbuNIzSUpNTdWoUaP061//2t2lehVX+nb06FH16NFD\nc+bM0ZYtW9SxY0c9+eST6tWrlydK9xhXejZixAht2bJFY8eOVfPmzdWsWTO9+eabat26tSdK96gx\nY8aYmsdaYM9s31gLfmK2ZxJrQTWzPauvdaBRvwNRXl6ugoICp18lJSWSpJYtW9rmV39fUVFR53Gn\nT5+uoKAg3X///Q33BLxAaWmpXY+kmvtUVlbmdG5d/WyMXOnbpYqLizV9+nRFRERo0KBBDVqjt3Gl\nZ1988YV27typRx991G31eStX+lZSUqK33npLQUFBeuuttxQZGamHH35Yp06dclu93sCVnp07d05W\nq1Vz585Vdna24uPjNXv2bBUXF7ut3isNa8Ev15TXAlewFriuvtaBRv0OxO7duzV+/HhZLBaHbSkp\nKZJ+XCx+vnC0atWqxmOWlJRoypQpKigo0DvvvCNfX98GqNx7+Pr6OvyjX1Ofaprr5+fXsEV6IVf6\nVs1qtWrChAmyWCx6+eWXG7xGb2O2Z+Xl5Xruuec0d+5ch19SmiJXXmvNmzdX165dNW3aNElSSEiI\ntm7dqnXr1mnSpEnuKdgLuNKztLQ0denSxfa/e/Pnz1dsbKzWrl2r5ORk9xR8hWEt+GWa+lpgFmvB\n5amvdaBRB4g+ffooPz/f6bbTp08rLS1NVqvVdmlS9WVNgYGBTve5cOGCkpOTVVhYqJUrV6pTp04N\nVru3aN++vc6dO6eqqio1a/bjG1ZWq1V+fn5q06aNw9yioiK7MavVWmM/GzNX+iZJp06d0vjx49W8\neXNlZmY6XK7TFJjt2Z49e3T06FFNnz7d7jrPiRMnKj4+Xs8995y7S/coV15rgYGBDpdJ3HTTTTpx\n4oTb6vUGrvRs3759Gj9+vO2xxWJRSEiIjh8/7taarySsBZePtcA81oLLU1/rQKO+hKk2QUFBCg4O\n1vbt221jeXl5Cg4Odrj/QZIMw9C0adN07NgxZWVl6ZZbbnFnuR7TtWtX+fj4aNeuXbaxvLw8de/e\n3WFuaGiow01yO3futLtRsalwpW+lpaVKTk5WixYtlJWV5fT11xSY7VloaKg2bdqkdevW6aOPPtJH\nH30kSfrzn//c6O9JcsaV11pYWJjDf6p88803uuGGGxq8Tm/iSs+CgoIcrt3/73//q44dOzZ4nVcq\n1oLLw1rgGtaCy1Nf60CTDRCSlJiYqLS0NOXm5urLL79Uenq6HnjgAdv24uJi270S2dnZys3N1YIF\nC3TNNdfIarXKarXq22+/9VT5buHn56e4uDjNnTtXe/fuVU5OjlasWGHrk9VqVXl5uSRp2LBhOn/+\nvBYuXKgjR45owYIFKikpUWxsrCefgke40rfXX39dhYWFev7551VVVWV7bTW1T94w27OWLVuqU6dO\ndl/Sj7/otWvXzpNPwSNcea0lJibq4MGDevXVV1VQUKCXX35ZhYWFuueeezz5FNzOlZ7dd999ys7O\n1rp161RQUKC0tDSdOHFC8fHxnnwKXoe14PKwFriOtcB1DbIOuPYps41LZWWlkZqaavTp08e44447\njPT0dLvtd955p/HKK68YhvHTHyj5+Vdtn1HfWJSWlhqzZ882wsPDjf79+xsZGRm2bV26dLH7bO89\ne/YYo0aNMkJDQ42EhATjwIEDnijZK5jt2/Dhw52+tmbPnu2p0j3GldfapZryZ38bhmt927FjhzFq\n1CijZ8+exqhRo4zt27d7omSPc6Vnq1evNmJjY41evXoZSUlJTfrftWo//5ljLTCntr6xFjhX12ut\ntrlNVV09q491wGIYTfhDcwEAAAC4pElfwgQAAADANQQIAAAAAKYRIAAAAACYRoAAAAAAYBoBAgAA\nAIBpBAgAAAAAphEgAAAAAJhGgAAAAABgGgECAAAAgGkECAC4gqxdu1bjxo1TVFSUevTooaFDh2rh\nwoWyWq2/+NivvvqqBg0aVA9VOvfUU09p/PjxDXb8alOnTlVOTk6Dn6daRkaGFi5c6LbzAYCnWQzD\nMDxdBACgdoZhaOrUqdq+fbumTJmi6OhoXX311Tp06JCWLl2q48eP64MPPlC7du0u+xylpaUqKyvT\nddddV4+V/+TChQuqqqpSmzZtGuT4kvTJJ5/o/fffV0ZGRoOd4+cqKys1cuRILVy4UL169XLbeQHA\nUwgQAHAFWL58udLT07V69WqFhITYbSsvL9ddd92lYcOGKSUlxUMVel5VVZWGDRumOXPmKDo62q3n\nXrVqlT799FNlZWW59bwA4AlcwgQAV4CsrCzFx8c7hAdJ8vX1VUZGhh5//HFJ0rFjxxQSEqI333xT\n/fr105AhQ/T999/r66+/1iOPPKI+ffqoe/fuGjx4sFasWGE7ziuvvKKYmBi7Y2zatEkJCQnq2bOn\nYmJi9P7779dYY1VVlRYvXqyBAweqR48eio2N1bvvvmvbPnv2bNslTOPGjVNISIjD1wMPPGCbv2bN\nGo0YMUKhoaEaOXKkMjIyVNv/eW3cuFHfffed+vbtaxuLiYnRsmXLNHnyZIWFhSkmJkY5OTnavHmz\nhg8frvDwcCUnJ6u4uFiSlJubq27duiknJ0fDhw9XaGioHnzwQZ08eVILFixQ79691bdvX73++ut2\n5x4+fLh27typr776qsb6AKCxIEAAgJc7evSojh8/rqioqBrnBAcHq0WLFnZjH374oTIyMvTSSy+p\nefPmeuihh3Tdddfpvffe04YNGxQbG6tFixYpPz9fkmSxWGSxWOyOkZqaqkcffVQbNmzQnXfeqXnz\n5unYsWNOa1i1apU2bdqkl19+WZs2bdLvf/97zZs3Tzt27LAdv9pf//pXbd261fb1zDPPyMfHR1Om\nTJEkvffee1q8eLGmT5+u9evXa8aMGVq2bJlefPHFGnuwefNm9e3bV82bN7cbX7p0qUaOHKmPP/5Y\nXbt21axZs/TGG2/oxRdf1BtvvKG9e/dq2bJltvmVlZV6/fXXlZ6eroyMDB04cEBxcXHy9fXV6tWr\nlZiYqJdeekmHDh2y7ePv76/u3btr8+bNNdYHAI0FAQIAvNyZM2ckyeH+hkceeUTh4eG2r7vvvttu\ne1JSkm655RZ169ZNJSUlevDBBzVnzhz96le/UufOnTVt2jRJ0tdff13juSdMmKCBAweqY8eOmjFj\nhiorK7V7926nc48ePapWrVqpQ4cOCg4OVlJSkpYvX66bbrrJYW6bNm3k7+8vf39/HT16VGlpaZoz\nZ47uuOMOSdJrr72mRx99VLGxserYsaOGDBmimTNnKjMzUxUVFU7Pv3v3bt16660O43feeafuuece\nderUSQkJCSopKdHMmTPVrVs39enTR3379rULA5I0Y8YM3XbbbQoNDdUdd9yhq666Sk888YRuvPFG\nTZ48WZIc9vnNb36jXbt21dhLAGgsfDxdAACgdtU3NZ87d85u/E9/+pPKysokSStXrtQ//vEPu+2d\nO3e2fd+uXTuNGTNGH3/8sfbv36+CggLl5+fLYrGoqqqqxnPffPPNtu9bt24tSTX+Ap+UlKScnBwN\nGDBAXbt2Vb9+/TRixIhab+wuLCzU1KlTNWbMGCUkJEiSiouLdfLkSaWnp2vJkiW2uYZh6IcfflBh\nYaFdXdWsVqv8/f0dxi/tQ6tWrSRJnTp1so35+fnZQpr04zsll+5z1VVXqWPHjrbHvr6+TvvQrl27\nGsMVADQmBAgA8HKdOnVSYGCgcnNzFRsbaxsPDAy0fd+2bVuH/fz8/GzfW61WJSQkKCAgQDExMfrt\nb3+rHj16aMCAAbWeu2XLlqbrvPHGG/XZZ58pNzdXW7du1eeff65ly5bp+eefV3x8vMP8CxcuaPLk\nyerRo4dmzZplG6++z+Hpp592etlWcHCw0/PXFIZ8fByXumbNan8D/uf7/PzSLmcuXrxY53EBoDHg\nXzoA8HLNmjXTuHHj9OGHH+rgwYNO5xw/frzWY3zyySf67rvv9O677+qRRx7R4MGDbe9o1NeH8WVm\nZmrjxo2KiopSSkqKPvroI0VFRenTTz91mFtZWanHHntMzZs3V3p6ut0v6NWXNhUUFKhTp062r717\n92rJkiU11hsYGGi7GdoTzp49axfqAKCx4h0IALgCTJw4Ufn5+UpKStLEiRM1YMAAtW7dWgcPHtSq\nVav0xRdfaPTo0TXuf/3116u0tFQbNmxQRESEjhw5otTUVFkslhovSXJVcXGxli5dKj8/P4WEhOjI\nkSM6cOCAHnzwQYe58+bN08GDB7V8+XKVlpaqpKTEti0gIEDJycl66aWXFBwcrP79+ys/P1/z5s3T\n4MGDHW4WrxYaGqp9+/bVWWddgelyA9W+ffs0ZMiQy9oXAK4kBAgAuAJYLBalp6dr48aNWrNmjTIz\nM/Xtt98qMDBQkZGRysrKUkREhN38Sw0fPlz79+9Xamqqvv/+e91www0aPXq0Nm/erL179+r+++93\nek4zY9WmT5+uixcvasGCBbJarQoICFBSUpImTZrksP/7778vi8Vid2mTYRiyWCw6cOCAJkyYID8/\nP2VmZio1NVWBgYFKTEy03fjtzODBg/Xss8+qsrLS9klMrj4HM9udzSkuLtbhw4e1aNGiOvcFgCsd\nf0gOANAoXLx4UcOHD9eTTz6poUOHuvXcb7/9trZs2aJVq1a59bwA4AncAwEAaBR8fHw0bdo0uz+O\n5w4VFRV69913NXPmTLeeFwA8hQABAGg04uPjde2112rTpk1uO+ff/vY3DRgwQJGRkW47JwB4Epcw\nAQAAADCNdyAAAAAAmEaAAAAAAGAaAQIAAACAaQQIAAAAAKYRIAAAAACYRoAAAAAAYBoBAgAAAIBp\nBAgAAAAAphEgAAAAAJj2f7Nv8cLOjrlOAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_log_reg(x='Grain size (mm)', y='Spiders', data=df, clf=clf, xmin=0, xmax=1.5);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Hypothesis testing\n", "\n", "To test if *Grain size* is a significant factor, we use the [**likelihood ratio test**](https://en.wikipedia.org/wiki/Logistic_regression#Evaluating_goodness_of_fit).\n", "\n", "We calculate the likelihood of the model with the grain size (the alternative model):" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def log_reg_null_model(X, y):\n", " X_ = np.zeros_like(X)\n", " if X_.ndim == 1:\n", " X_ = X_.reshape(-1, 1)\n", " clf = sklearn.linear_model.LogisticRegression(C=1e12)\n", " clf.fit(X_, y)\n", " return clf\n", "\n", "clf0 = log_reg_null_model(df['Grain size (mm)'], df['Spiders'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The likelihood ratio test operates by calculating the test statistic $D$ from the likelihoods of the null and alternative models:\n", "$$\n", "D = -2 \\log{ \\frac{L(H_0)}{L(H_1)} }\n", "$$\n", "The test statistic is then approximately chisquare distributed.\n", "\n", "*scikit-learn* has a log-loss function that can help us do that. \n", "The log-loss is defined as the negative log-likelihood, so we can rewrite:\n", "$$\n", "D = 2 (-\\log{L(H_0)} + \\log{L(H_1)}) \\Rightarrow \\\\\n", "D = 2 (logloss(H_0) - logloss(H_1))\n", "$$" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import sklearn.metrics" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def log_reg_lik_ratio_test(X, Y, clf0, clf1, df=1):\n", " if X.ndim == 1:\n", " X = X.reshape(-1, 1)\n", " y_prob0 = clf0.predict_proba(X)\n", " loss0 = sklearn.metrics.log_loss(Y, y_prob0, normalize=False)\n", " y_prob1 = clf1.predict_proba(X)\n", " loss1 = sklearn.metrics.log_loss(Y, y_prob1, normalize=False)\n", " D = 2 * (loss0 - loss1)\n", " return scipy.stats.distributions.chi2.sf(D, df=df)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.033243767135570736" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "log_reg_lik_ratio_test(df['Grain size (mm)'], df['Spiders'].astype(np.float64), clf0, clf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "John indeed reports 0.033.\n", "\n", "Note that the log-loss calculation in equivalent to:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "15.3157739197 15.3157739197\n" ] } ], "source": [ "_ = clf.predict_proba(df['Grain size (mm)'].reshape(-1, 1))\n", "df['prob_absent'], df['prob_present'] = _[:,0], _[:,1]\n", "lik = df.loc[df['Spiders'], 'prob_present'].prod() * df.loc[~df['Spiders'], 'prob_absent'].prod()\n", "print(\n", " -np.log(lik), \n", " sklearn.metrics.log_loss(\n", " df['Spiders'], \n", " clf.predict_proba(df['Grain size (mm)'].reshape(-1, 1)), \n", " normalize=False\n", " )\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Second example" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
LocationLatitudeMpi90Mpi100p, Mpi100
0Port Townsend, WA48.1471390.748
1Neskowin, OR45.21772410.577
2Siuslaw R., OR44.0108711830.521
3Umpqua R., OR43.71871750.483
4Coos Bay, OR43.53976710.628
\n", "
" ], "text/plain": [ " Location Latitude Mpi90 Mpi100 p, Mpi100\n", "0 Port Townsend, WA 48.1 47 139 0.748\n", "1 Neskowin, OR 45.2 177 241 0.577\n", "2 Siuslaw R., OR 44.0 1087 1183 0.521\n", "3 Umpqua R., OR 43.7 187 175 0.483\n", "4 Coos Bay, OR 43.5 397 671 0.628" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = \"\"\"Location\tLatitude\tMpi90\tMpi100\tp, Mpi100\n", "Port Townsend, WA\t48.1\t47\t139\t0.748\n", "Neskowin, OR\t45.2\t177\t241\t0.577\n", "Siuslaw R., OR\t44\t1087\t1183\t0.521\n", "Umpqua R., OR\t43.7\t187\t175\t0.483\n", "Coos Bay, OR\t43.5\t397\t671\t0.628\n", "San Francisco, CA\t37.8\t40\t14\t0.259\n", "Carmel, CA\t36.6\t39\t17\t0.304\n", "Santa Barbara, CA\t34.3\t30\t0\t0\n", "\"\"\"\n", "df = pd.read_table(io.StringIO(data))\n", "df.head()" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArkAAAHqCAYAAAAArFZwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X1Y1GWi//HPEMLgCCqgLqSefCpAVJRI2bXtBK5t+YS7\naltZLGl62nTbdrV86ISWofZ0bSe31PKwpm1XZW4+bWbqbpaxrY8rMVoKVnIkElAoAsZkfn/0c3Yn\nkGZ0YIab9+u6vGLuuYf54PVt/Pj1/t5fi9PpdAoAAAAwSJC/AwAAAAC+RskFAACAcSi5AAAAMA4l\nFwAAAMah5AIAAMA4lFwAAAAYh5ILAAAA47Spkjt58mRNnjzZ3zEAAADQzIL9HaAllZSU+DsCAAAA\nWkCbOpMLAACAtoGSCwAAAONQcgEAAGAcSi4AAACMQ8kFAACAcSi5AAAAMA4lFwAAAMah5AIAAMA4\nlFwAAAAYh5ILAAAA41ByAQAAYBxKLgAAAIxDyQUAAIBxKLkAAAAwDiUXAAAAxqHkAgAAwDiUXAAA\nABiHkgsAAADjUHIBAABgHEouAAAAjEPJBQAAgHEouQAAADAOJRcAAADGoeQCAADAOJRcAAAAGIeS\nCwAAAONQcgEAAGAcSi4AAACMQ8kFAACAcQKq5DocDo0ZM0Z79uy54By73a5JkyYpKSlJEydOVEFB\nQQsmBAAAQGsQMCXX4XDot7/9rY4dO3bBOTU1NZo2bZpSUlK0fv16JSUlafr06aqtrW3BpAAAAAh0\nAVFyCwsLNWnSJBUXFzc5b8uWLQoLC9Ps2bPVu3dvzZ8/XzabTVu3bm2hpAAAAGgNAqLk/uMf/1Bq\naqpeeeUVOZ3OC847dOiQkpOT3caGDBmiAwcONHdEAAAAtCLB/g4gSbfccotH87744gtdeeWVbmNR\nUVFNLnEAAABA2xMQZ3I9VVtbq5CQELexkJAQORwOPyUCAABAIGpVJTc0NLRBoXU4HLJarX5KBAAA\ngEDUqkput27ddOrUKbexsrIydenSxU+JAAAAEIhaVckdNGhQg4vMDhw4oKSkJD8lAgAAQCAK+JJb\nVlamuro6SdINN9ygL7/8Ujk5OSosLNSiRYv09ddf68Ybb/RzSgAAAASSgCu5FovF7fHw4cP15ptv\nSpI6dOig5cuXa+/evfr5z3+u/Px8Pf/886zJBQAAgBuLs6mNaQ2Tnp4uSdqxY4efkwAAAKA5BdyZ\nXAAAAOBSUXIBAABgHEouAAAAjEPJBQAAgHEouQAAADAOJRcAAADGoeQCAADAOJRcAAAAGIeSCwAA\nAONQcgEAAGAcSi4AAACMQ8kFAACAcSi5AAAAMA4lFwAAAMah5AIAAMA4lFwAAAAYh5ILAAAA41By\nAQAAYBxKLgAAAIxDyQUAAIBxKLkAAAAwDiUXAAAAxqHkAgAAwDiUXAAAABiHkgsAAADjUHIBAABg\nHEouAAAAjEPJBQAAgHEouQAAADAOJRcAAADGoeQCAADAOJRcAAAAGIeSCwAAAONQcgEAAGAcSi4A\nAACMQ8kFAACAcSi5AAAAMA4lFwAAAMah5AIAAMA4lFwAAAAYh5ILAAAA41ByAQAAYBxKLgAAAIxD\nyQUAAIBxKLkAAAAwDiUXAAAAxqHkAgAAwDiUXAAAABiHkgsAAADjUHIBAABgHEouAAAAjEPJBQAA\ngHEouQAAADAOJRcAAADGoeQCAADAOJRcAAAAGIeSCwAAAONQcgEAAGAcSi4AAACMExAl1+FwaN68\neUpJSdG1116r3NzcC859++23NWrUKA0ePFi33Xab7HZ7CyYFAABAaxAQJXfp0qWy2+1as2aNsrOz\ntWzZMm3btq3BvGPHjmnWrFmaPn26Nm7cqLi4OE2bNk11dXV+SA0AAIBA5feSW1NTo3Xr1unBBx9U\nXFycRowYoalTp2rt2rUN5r733nvq16+fxo4dqx49eui3v/2tysrKdOzYMT8kBwAAQKDye8k9cuSI\nzp07p6SkJNdYcnKyDh061GBup06ddOzYMe3fv19Op1Ovv/66wsPD1bNnz5aMDAAAgAAX7O8Ap06d\nUqdOnRQc/K8oUVFRqqur0+nTp9W5c2fX+E033aSdO3fq1ltv1WWXXaagoCCtXLlS4eHh/ogOAACA\nAOX3M7k1NTUKCQlxGzv/2OFwuI2fOXNGZWVlys7O1muvvaaMjAzNmTNHFRUVLZYXAAAAgc/vJTc0\nNLRBmT3/OCwszG38iSee0FVXXaVbbrlFCQkJevjhhxUWFqb169e3WF4AAAAEPr+X3G7duunMmTOq\nr693jZWVlclqtSoiIsJtbkFBgeLi4lyPLRaL4uLidPLkyRbLCwAAgMDn95IbHx+v4OBgHTx40DW2\nd+9eJSYmNpjbtWvXBjspHD9+XN27d2/2nAAAAGg9/F5yrVarxo0bp+zsbOXn52v79u3Kzc1VZmam\npG/P6p7fB3fixIl67bXXtGHDBn322Wd64oknVFJSooyMDH/+CAAAAAgwft9dQZLmzp2rhQsXKjMz\nU+Hh4br33ns1YsQISdLw4cO1ZMkSZWRk6KabblJNTY1WrFih0tJSxcfH68UXX1RkZKSffwIAAAAE\nEovT6XT6O0RLSU9PlyTt2LHDz0kAAADQnPy+XAEAAADwNUouAAAAjEPJBQAAgHEouQAAADAOJRcA\nAADGoeQCAADAOJRcAAAAGCcgbgYBAACAxjmdThUUlauiqlaREVb17x0li8Xi71gBj5ILAAAQoPLy\nTyp3k10l5dWusZgom7LGJCh1QKwfkwU+lisAAAAEoLz8k1qyeo9bwZWkkvJqLVm9R3n5J/2UrHWg\n5AIAAAQYp9Op3E121Tsbf77eKeVutsvpvMAEUHIBAAACTUFReYMzuN9VUlYt+/GKFkrU+lByAQAA\nAkxFVa1n8yo9m9cWUXIBAAACTGSE1bN5HT2b1xZRcgEAAAJM/95RiomyNTknJtqmhF6RLZSo9aHk\nAgAABBiLxaKsMQkKusB2uEEWKWt0AvvlNoGSCwAAEIBSB8RqTmaKYqLdz+jGRNs0JzOFfXK/BzeD\nAAAACFCpA2I1LDFGBUXlOl1Vp8iOViX0iuQMrgcouQAAAAHMYrEosU+0v2O0OixXAAAAgHEouQAA\nADAOJRcAAADGoeQCAADAOJRcAAAAGIeSCwAAAONQcgEAAGAcSi4AAACMQ8kFAACAcSi5AAAAMA4l\nFwAAAMah5AIAAMA4lFwAAAAYh5ILAAAA41ByAQAAYBxKLgAAAIxDyQUAAIBxKLkAAAAwDiUXAAAA\nxqHkAgAAwDiUXAAAABiHkgsAAADjUHIBAABgnGBvX1BUVKSHH35Y+/fv19mzZxs8f/jwYZ8EAwAA\nAC6W1yU3Oztb5eXl+t3vfqeIiIjmyAQAAABcEq9L7j//+U+9/PLL6t+/f3PkAQAAAC6Z12tyO3fu\nrHbt2jVHFgAAAMAnvC65kydP1lNPPaWvvvqqOfIAAAAAl8zr5Qrvv/++9u7dq2uuuUZRUVEKCQlx\ne37Hjh0+CwcAAABcDK9LbnJyspKTk5sjCwAAAOATXpfcGTNmNEcOAAAAwGe8LrmS9OGHH2rVqlX6\n+OOPFRwcrL59+yozM1MDBw70dT4AAIzhdDpVUFSuiqpaRUZY1b93lCwWi79jAUbyuuT+4x//0J13\n3qkrr7xSP/rRj1RfX6/9+/fr1ltv1erVq1nKAABAI/LyTyp3k10l5dWusZgom7LGJCh1QKwfkwFm\nsjidTqc3L7jlllt05ZVXauHChW7jCxcu1LFjx7RmzRqfBvSl9PR0SVwcBwBoWXn5J7Vk9R7VN/In\nbpBFmpOZQtEFfMzrLcTsdrvuuOOOBuOTJ0/Whx9+6JNQAACYwul0KneTvdGCK0n1Til3s11ennMC\n8D0u6mYQp0+fbjBeUVHRYDsxAADauoKicrclCo0pKauW/XhFCyUC2gavS+7111+vRx55RIWFha6x\nY8eOadGiRUpLS/NpOAAAWruKqlrP5lV6Ng+AZ7y+8Ow3v/mNsrKyNHr0aIWHh8tisaiqqkpxcXG6\n//77myMjAACtVmSE1bN5HT2bB8AzXpfcjh07at26dXr33Xd19OhROZ1OXXXVVRo+fLiCgrw+MQwA\ngNH6945STJStySULMdE2JfSKbMFUgPkuap/coKAgXXfddbruuut8nQcAAKNYLBZljUlocneFrNEJ\n7JcL+JhHp17j4+NVXl4uSYqLi1N8fPwFf10Mh8OhefPmKSUlRddee61yc3MvOPejjz7SrbfeqkGD\nBmns2LH64IMPLuo9AQBoKakDYjUnM0Ux0Ta38ZhoG9uHAc3EozO5OTk5Cg8PlyQtXrzY5yGWLl0q\nu92uNWvWqLi4WA888IAuv/xyjRw50m3eV199pSlTpig9PV1Lly7VG2+8oRkzZuitt95SZCT/zAMA\nCFypA2I1LDFGBUXlOl1Vp8iOViX0iuQMLtBMPCq548ePd31tsVh00003Ndgu7Ouvv9arr77qdYCa\nmhqtW7dOq1atUlxcnOLi4jR16lStXbu2Qcldv369bDab60YUM2fO1K5du/Thhx/qxz/+sdfvDQBA\nS7JYLErsE+3vGECb4FHJraioUG3tt1ubzJ07V/369VPnzp3d5hw+fFhPPfWUfvnLX3oV4MiRIzp3\n7pySkpJcY8nJyVqxYkWDuXv27GmwTdlrr73m1fsBAADAfB6V3F27dmnOnDmyWCxyOp2aMGFCgzlO\np/OiLkQ7deqUOnXqpODgf0WJiopSXV2dTp8+7VamT5w4oQEDBuihhx7Szp071b17d91///0aMmSI\n1+8LAAAAc3lUcjMyMnT55Zervr5emZmZ+p//+R917NjR9bzFYlH79u115ZVXeh2gpqamwdKH848d\nDofb+Ndff60XXnhBd9xxh1544QVt3rxZU6ZM0datW9WtWzev3xsAAABm8ngLsZSUFEnSrFmzlJiY\nqJiYGJ8ECA0NbVBmzz8OCwtzG7/ssssUHx+vGTNmSPp2p4fdu3drw4YNmjZtmk/yAAAAoPXz+u4N\ny5cvb1BKL0W3bt105swZ1dfXu8bKyspktVoVERHhNrdLly7q3bu329gVV1yhkpISn+UBAABA6+d1\nyR00aJB27tzpswDx8fEKDg7WwYMHXWN79+5VYmJig7lJSUk6cuSI21hRUZEuv/xyn+UBAABA6+f1\nHc86dOigxx57TMuXL9cVV1yh0NBQt+dffPFFr76f1WrVuHHjlJ2drZycHJWWlio3N1dLliyR9O1Z\n3fDwcIWGhuoXv/iF1q5dq2XLlmns2LH685//rOLiYo0dO9bbHwMAAFwip9OpgqJyVVTVKjLCqv69\no9j3FwHD65Lbvn17ZWRk+DTE3LlztXDhQmVmZio8PFz33nuvRowYIUkaPny4lixZooyMDMXGxmrV\nqlV65JFH9Pzzz6tPnz56/vnn1bVrV5/mAQAATcvLP6ncTXaVlFe7xmKibMoak8Ad3BAQLE6ns5E7\naZspPT1dkrRjxw4/JwEAoPXKyz+pJav3qL6RBhFkEbcqRkDw+kyuJJWUlOill17Sxx9/rODgYPXr\n108333yzYmM5oAEAMJnT6VTuJnujBVeS6p1S7ma7hiXGsHQBfuX1hWcfffSRxo4dqw0bNqhdu3Zy\nOp1av369xo4dq6NHjzZHRgAAECAKisrdlig0pqSsWvbjFS2UCGic12dyH3vsMQ0dOlRPPvmk66Kz\nuro6zZo1S0888USjt+MFAABmqKiq9WxepWfzgObi9Znc/fv3a+bMmW67KoSGhuqee+7Rvn37fBoO\nAAAElsgIq2fzOno2D2guXpdcm82ms2fPNhhvbAwAAJilf+8oxUTZmpwTE21TQq/IFkoENM7rkjts\n2DA99thjOnPmjGusoqJCjz/+uFJTU30aDgAABBaLxaKsMQkKusA1ZUEWKWt0Ahedwe+83kLs888/\n1y9+8QtVVlbqiiuukMVi0fHjx9WpUyetWbNG3bt3b66sl4wtxAAA8I28/JPK3WxXSdm/7ZMbbVPW\naPbJRWC4qH1yq6urtWHDBh09elROp1NXXXWVxowZow4dOjRHRp+h5AIA4Dvn73h2uqpOkR2tSugV\nyRlcBIyL2ifXZrMpIyNDx48f12WXXaZevXo1uL0vAAAwm8ViUWKfaH/HABrldck9e/ascnJy9Prr\nr7suNrNarbrjjjt03333+TwgAAAA4C2vS+5TTz2lrVu3at68eRo8eLDq6+u1f/9+PfPMMwoLC9N/\n/dd/NUdOAAAAwGNel9wNGzYoJydH119/vWssPj5eXbp0UU5ODiUXAAAAfuf1FmJ1dXXq2bNng/G+\nffuqsrLSJ6EAAACAS+F1yc3IyNDTTz8th8PhGnM6nVq9erXGjx/v03AAAADAxfB6ucKZM2f017/+\nVWlpaRo4cKCCg4Nlt9v1f//3fxo0aJDuuOMO19wXX3zRp2EBAAAAT3hdckNCQjR69Gi3sZSUFKWk\npPgsFAAAAHApvC65ixcvbo4cAAAAgM9c1M0gSkpK9NJLL+njjz9WcHCw+vXrp5tvvlmxsdzGDwCA\n5nD+7mIVVbWKjLCqf+8o7i4GNMHrkvvRRx9p8uTJslqtGjhwoOrr67V+/Xq99NJLevnll9WvX7/m\nyAkAQJuVl39SuZvsKimvdo3FRNmUNSZBqQM4wQQ0xuJ0Op3evGDKlCkKCwvTk08+6bqVb11dnWbN\nmiWHw6EVK1Y0S1BfSE9PlyTt2LHDz0kAAPBMXv5JLVm9R/WN/GkdZJHmZKZQdIFGeL2F2P79+zVz\n5kxXwZWk0NBQ3XPPPdq3b59PwwEA0JY5nU7lbrI3WnAlqd4p5W62y8vzVUCb4HXJtdlsOnv2bIPx\nxsYAAMDFKygqd1ui0JiSsmrZj1e0UCKg9fC65A4bNkyPPfaYzpw54xqrqKjQ448/rtTUVJ+GAwCg\nLauoqvVsXqVn84C2xOsLz373u9/plltu0fXXX68rrrhCFotFx48fV6dOnZSTk9McGQEAaJMiI6ye\nzevo2TygLfG65MbExGjLli3asGGDjh49KqfTqYkTJ2rMmDHq0KFDc2QEAKBN6t87SjFRtiaXLMRE\n25TQK7IFUwGtg9cld8KECVq0aJFuvfXW5sgDAAD+P4vFoqwxCU3urpA1OoH9coFGeL0m98SJE2rf\nvn1zZAEAAN+ROiBWczJTFBNtcxuPibaxfRjQBK/3yX3++ee1a9cuTZkyRT179pTV6r4OKJDvesY+\nuQCA1ur8Hc9OV9UpsqNVCb0iOYMLNMHrktu/f3+dO3fu2xf/2/9cTqdTFotFhw8f9m1CH6LkAgAA\ntA1er8nNzc1tjhwAAACAz3hdcq+55hrX12fOnNFll12m8PBwn4YCAKAtOb8UoaKqVpERVvXvHcVS\nBOASeV1yJemFF17Qiy++qFOnTkmSunfvrrvuukuTJk3yaTgAAEyXl39SuZvsbtuExUTZlDUmgYvK\ngEvgdclduXKlnn32Wd1+++0aPHiw6uvrtW/fPuXk5MjpdOrmm29ujpwAABgnL/9ko9uDlZRXa8nq\nPeyeAFwCr0vuSy+9pAULFigjI8M1NmLECPXp00crV66k5AIA4AGn06ncTfZG97+VpHqnlLvZrmGJ\nMSxdAC6C1/vkVlZWatCgQQ3GU1JSVFpa6pNQAACYrqCovMk7mUlSSVm17McrWigRYBavS256errW\nrFnTYHzTpk1KS0vzSSgAAExXUVXr2bxKz+YBcOf1coWoqCi9/PLL2rdvn6655hoFBwfrww8/1N69\ne5Wenq65c+e65i5evNinYQEAMEVkhPX7J0mK7OjZPADuvC65hw8fVlJSkiTpyJEjrvGrr75alZWV\nqqys9F06AAAM1b93lGKibE0uWYiJtimhV2QLpgLM4XXJbWypAgAA8I7FYlHWmIRGd1eQpCCLlDU6\ngYvOgIvk9ZpcAADgG6kDYjUnM0Ux0Ta38ZhoG9uHAZfoom4GAQAAfCN1QKyGJcaooKhcp6vqFNnR\nqoRekZzBBS4RJRcAAD+zWCxK7BPt7xiAUViuAAAAAON4fSa3oqJCJSUlqqurU/v27dW1a1dFRnLl\nJwAAAAKHxyX3z3/+s1auXKlPPvlE0re3I5S+/SeWXr16afr06Ro3blyzhAQAAAC84VHJXbNmjZ54\n4gn98pe/1LBhw9S1a1eFhITI4XDoiy++0Pvvv6/s7GxVV1fr1ltvbe7MAAAAQJM8Krl//OMftWDB\nAo0fP77Bc3369FFqaqp69eqlZ599lpILAAAAv/PowrOKigoNGjSoyTmDBg3SqVOnfBIKAAAAuBQe\nldzExETl5uaqvr6+0eedTqdeeOEFxcfH+zQcAAAAcDE8Wq4wb9483Xnnndq1a5dSUlIUExPjtiZ3\n7969+uqrr7Rq1armzgsAAAB8L4vz/DYJ3+PMmTN65ZVXtG/fPn3++eeqra1VaGioYmJidPXVV2vC\nhAkBv5VYenq6JGnHjh1+TgIAAIDm5PEWYp06ddL06dObMwsAAADgEx6X3G+++Ubbtm3Tnj17VFJS\nIofDobCwMHXr1k0pKSn6yU9+ouBg7hIMAAAA//PowrPi4mKNGjVK8+bN00cffSSr1aouXbqoXbt2\nOnLkiObOnauxY8fq5MmTzZ0XAAAA+F4ercmdNm2azp07p9///vcKDw9v8HxVVZXuu+8+tWvXTsuX\nL2+WoL7AmlwAAIC2waMzuXv27NH999/faMGVpIiICM2ePVt79uzxaTgAAADgYnhUcsPDw1VaWtrk\nnJMnT8pqtfokFAAAAHApPCq5EyZM0Jw5c/TKK6/o008/lcPhkCQ5HA6dOHFCr7/+uubPn6+f/exn\nzRoWAAAA8IRH2yHMnDlTQUFBeuyxx/T11183eN5ms+m2227Tvffe6/OAAAAAgLc8vhmEJJ09e1aH\nDx9WaWmpampqZLVa9YMf/EBxcXEKCQm56BAOh0MLFizQ22+/LavVqjvvvFNZWVlNvqa4uFhjxozR\nypUrlZKS4tH7cOEZAABA2+DVxrbt2rXTwIEDfR5i6dKlstvtWrNmjYqLi/XAAw/o8ssv18iRIy/4\nmgULFqi2ttbnWQAAAND6ebQmtznV1NRo3bp1evDBBxUXF6cRI0Zo6tSpWrt27QVfs3HjxkaXTQAA\nAACSh2dyvdkazNOlA+cdOXJE586dU1JSkmssOTlZK1asaHT+6dOn9eSTT2rVqlUaPXq0V+8FAACA\ntsGjkvvb3/5WZWVlkqSmlvBaLBYdPnzYqwCnTp1Sp06d3G4JHBUVpbq6Op0+fVqdO3d2m79kyRKN\nHz9effv29ep9AAAA0HZ4VHI3btyoKVOmKCgoSL///e9lsVh8FqCmpqbBRWvnH5/fquy8999/XwcO\nHNAjjzzis/cHAACAeTwquZ07d9Zzzz2ncePGKS8vTxMnTvRZgNDQ0AZl9vzjsLAw11hdXZ0WLFig\n7OzsS9rJAQAAAObz+MKzbt26adasWT7ffqtbt246c+aM6uvrXWNlZWWyWq2KiIhwjR06dEgnTpzQ\nzJkzNXjwYA0ePFiSdNddd2nBggU+zQQAAIDWzastxCZMmKAJEyb4NEB8fLyCg4N18OBBDRkyRJK0\nd+9eJSYmus0bNGiQtm3b5jb2k5/8RI8++qhSU1N9mgkAAACtm1cltzlYrVaNGzdO2dnZysnJUWlp\nqXJzc7VkyRJJ357VDQ8PV2hoqHr06NHg9V27dlVkZGRLxwYAAEAA86jkLlu2TFOmTFFYWJiWLVvW\n5NwZM2Z4HWLu3LlauHChMjMzFR4ernvvvVcjRoyQJA0fPlxLlixRRkZGg9f58gI4AAAAmMOj2/qm\npaXp9ddfV+fOnZWWlnbhb2axBPQtc7mtLwAAQNvg0ZncnTt3Nvo1AAAAEIguak1ufX298vLy9PHH\nHysoKEj9+/fX1Vdf7etsAAAAwEXxuuR+8cUXmjJlio4ePaqOHTvq3Llz+uqrrzRkyBCtWLFC4eHh\nzZETAAAA8JjH++Se9/DDDyskJER/+ctf9MEHH2jv3r3atGmTamtrlZOT0xwZAQAAAK94XXLff/99\nLVy4UL1793aN9evXTw899BAXdAEAACAgeF1ybTabzp4922C8Xbt2ateunU9CAQAAAJfC65L761//\nWg899JAOHz7sGjtx4oQWLVp0UXvkAgAAAL7m0T65/+66665TeXm5zp07J5vNpuDgYFVWVsrpdDa4\nOcO/F+FAwD65AAAAbYPXuyv85je/aY4cAAAAgM94XXLHjx/fHDkAAAAAn/Go5C5btsyjb2axWHTP\nPfdcUiAAAADgUnlccoOCgvSDH/ygyXmUXAAAAAQCj0rupEmT9Pbbb0uSRo0apVGjRikuLq5ZgwEA\nAAAXy+PdFc6dO6e///3v+stf/qLt27crMjJSo0eP1qhRo3TFFVc0c0zfYHcFAACAtsHrLcQk6ezZ\ns3rvvff05ptvaseOHerZs6duuukmjRo1SrGxsc2R0ycouQAAAG3DRZXcf+dwOPT666/rySefVHV1\ndcDtjfvvKLkAAABtg9dbiJ33xRdfaNu2bdq6dav27dun//iP/9Dtt9/uy2wAAADARfGq5JaWluqt\nt97S1q1bdeDAAfXo0UM33nijHnzwQS5EAwAAQMDwqOSuXr1aW7du1T//+U/Fxsbqxhtv1Pz589W/\nf//mzgcAAAB4zaM1uXFxcWrXrp1++MMfasCAAU3OnTFjhs/C+RprcgEAANoGj87knt8x4ejRozp6\n9OgF51ksloAuuQAAAGgbPCq5O3fubO4cAAAAgM8E+TsAAAAA4GsXvYUYALM5nU4VFJWroqpWkRFW\n9e8dJYvF4u9YAAB4hJILoIG8/JPK3WRXSXm1aywmyqasMQlKHRC4dzUEAOA8lisAcJOXf1JLVu9x\nK7iSVFJerSWr9ygv/6SfkgEA4DlKLgAXp9Op3E121V9gY8F6p5S72a5LvBs4AADNjpILwKWgqLzB\nGdzvKimrlv14RQslAgDg4lByAbhUVNV6Nq/Ss3kAAPgLJReAS2SE1bN5HT2bBwCAv1ByAbj07x2l\nmChbk3POjIawAAAWAklEQVRiom1K6BXZQokAALg4lFwALhaLRVljEhR0ge1wgyxS1ugE9ssFAAQ8\nSi4AN6kDYjUnM0Ux0e5ndGOibZqTmcI+uQCAVoGbQQBoIHVArIYlxqigqFynq+oU2dGqhF6RnMEF\nALQalFwAjbJYLErsE+3vGAAAXBSWKwAAAMA4lFwAAAAYh5ILAAAA41ByAQAAYBwuPAPQajidThUU\nlauiqlaREVb17x3Fjg8AgEZRcgG0Cnn5J5W7ya6S8mrXWEyUTVljEti7FwDQAMsVAAS8vPyTWrJ6\nj1vBlaSS8motWb1Hefkn/ZQMABCoKLkAAprT6VTuJrvqnY0/X++Ucjfb5XReYAIAoE2i5AIIaAVF\n5Q3O4H5XSVm17McrWigRAKA1oOQCCGgVVbWezav0bB4AoG2g5AIIaJERVs/mdfRsHgCgbaDkAgho\n/XtHKSbK1uScmGibEnpFtlAiAEBrQMkFENAsFouyxiQo6ALb4QZZpKzRCeyXCwBwQ8kFEPBSB8Rq\nTmaKYqLdz+jGRNs0JzOFfXIBAA1wMwgArULqgFgNS4xRQVG5TlfVKbKjVQm9IjmDCwBoFCUXQKth\nsViU2Cfa3zEAAK0AyxUAAABgHEouAAAAjEPJBQAAgHEouQAAADAOJRcAAADGoeQCAADAOJRcAAAA\nGIeSCwAAAOMERMl1OByaN2+eUlJSdO211yo3N/eCc//2t78pIyNDgwcP1rhx47Rz584WTAoAAIDW\nICBK7tKlS2W327VmzRplZ2dr2bJl2rZtW4N5H330kWbOnKmJEydq48aNmjRpkn7961/ro48+8kNq\nAAAABCq/l9yamhqtW7dODz74oOLi4jRixAhNnTpVa9eubTB38+bNSk1N1W233aYePXrotttu09Ch\nQ/Xmm2/6ITkAAAACVbC/Axw5ckTnzp1TUlKSayw5OVkrVqxoMHf8+PE6e/Zsg/GvvvqqWTMCAACg\ndfH7mdxTp06pU6dOCg7+V9+OiopSXV2dTp8+7Ta3d+/euuqqq1yPjx49qr///e9KTU1tsbwAAAAI\nfH4vuTU1NQoJCXEbO//Y4XBc8HUVFRWaOXOmkpOTlZ6e3qwZAQAA0Lr4veSGhoY2KLPnH4eFhTX6\nmrKyMmVmZspisejpp59u9owAAABoXfxecrt166YzZ86ovr7eNVZWViar1aqIiIgG80tLS3Xbbbfp\n3LlzWrNmjTp37tyScQEAANAK+L3kxsfHKzg4WAcPHnSN7d27V4mJiQ3m1tTUaOrUqWrXrp3Wrl2r\n6OjolowKAACAVsLvJddqtWrcuHHKzs5Wfn6+tm/frtzcXGVmZkr69qxuXV2dJGn58uUqLi7W4sWL\nVV9fr7KyMpWVlbG7AgAAANxYnE6n098hamtrtXDhQr311lsKDw/X1KlTdfvtt0uS4uLitGTJEmVk\nZOjGG2/UJ5980uD1GRkZWrx48fe+z/kL1Hbs2OHT/AAAAAgsAVFyWwolFwAAoG3w+3IFAAAAwNco\nuQAAADAOJRcAAADGoeQCAADAOJRcAAAAGIeSCwAAAONQcgEAAGAcSi4AAACMQ8kFAACAcSi5AAAA\nMA4lFwAAAMah5AIAAMA4lFwAAAAYh5ILAAAA41ByAQAAYBxKLgAAAIxDyQUAAIBxKLkAAAAwDiUX\nAAAAxqHkAgAAwDiUXAAAABiHkgsAAADjUHIBAABgHEouAAAAjEPJBQAAgHEouQAAADAOJRcAAADG\noeQCAADAOJRcAAAAGIeSCwAAAONQcgEAAGAcSi4AAACMQ8kFAACAcSi5AAAAMA4lFwAAAMah5AIA\nAMA4lFwAAAAYh5ILAAAA41ByAQAAYBxKLgAAAIxDyQUAAIBxKLkAAAAwDiUXAAAAxqHkAgAAwDiU\nXAAAABiHkgsAAADjUHIBAABgHEouAAAAjEPJBQAAgHEouQAAADAOJRcAAADGoeQCAADAOJRcAAAA\nGIeSCwAAAONQcgEAAGAcSi4AAACMQ8kFAACAcSi5AAAAMA4lFwAAAMYJiJLrcDg0b948paSk6Npr\nr1Vubu4F59rtdk2aNElJSUmaOHGiCgoKWjApAAAAWoOAKLlLly6V3W7XmjVrlJ2drWXLlmnbtm0N\n5tXU1GjatGlKSUnR+vXrlZSUpOnTp6u2ttYPqdHaOZ1OfVhYpl0HivVhYZmcTqe/IwEAAB8J9neA\nmpoarVu3TqtWrVJcXJzi4uI0depUrV27ViNHjnSbu2XLFoWFhWn27NmSpPnz52vXrl3aunWrMjIy\n/BEfrVRe/knlbrKrpLzaNRYTZVPWmASlDoj1YzIAAOALfj+Te+TIEZ07d05JSUmuseTkZB06dKjB\n3EOHDik5OdltbMiQITpw4ECz54Q58vJPasnqPW4FV5JKyqu1ZPUe5eWf9FMyAADgK34vuadOnVKn\nTp0UHPyvk8pRUVGqq6vT6dOn3eZ+8cUX6tq1q9tYVFSUSktLWyQrWj+n06ncTXbVX2BlQr1Tyt1s\nZ+kCAACtnN9Lbk1NjUJCQtzGzj92OBxu47W1tY3O/e484EIKisobnMH9rpKyatmPV7RQIgAA0Bz8\nXnJDQ0MblNTzj8PCwjyaa7VamzckjFFR5dlFihWVXMwIAEBr5veS261bN505c0b19fWusbKyMlmt\nVkVERDSYe+rUKbexsrIydenSpUWyovWLjPDsL0SRHfmLEwAArZnfS258fLyCg4N18OBB19jevXuV\nmJjYYO6gQYMaXGR24MABt4vWgKb07x2lmChbk3Niom1K6BXZQokAAEBz8HvJtVqtGjdunLKzs5Wf\nn6/t27crNzdXmZmZkr49U1tXVydJuuGGG/Tll18qJydHhYWFWrRokb7++mvdeOON/vwR0IpYLBZl\njUlQkKXx54MsUtboBFksF5gAAABaBb+XXEmaO3euEhMTlZmZqUceeUT33nuvRowYIUkaPny43nzz\nTUlShw4dtHz5cu3du1c///nPlZ+fr+eff541ufBK6oBYzclMUUy0+xndmGib5mSmsE8uAAAGsDjb\n0F5J6enpkqQdO3b4OQkCgdPpVEFRuU5X1Smyo1UJvSI5gwsAgCH8fsczwF8sFosS+0T7OwYAAGgG\nAbFcAQAAAPAlSi4AAACMQ8kFAACAcSi5AAAAMA4lFwAAAMah5AIAAMA4lFwAAAAYh5ILAAAA41By\nAQAAYBxKLgAAAIxDyQUAAIBxKLkAAAAwDiUXAAAAxqHkAgAAwDiUXAAAABiHkgsAAADjUHIBAABg\nHEouAAAAjEPJBQAAgHEouQAAADAOJRcAAADGoeQCAADAOJRcAAAAGIeSCwAAAONQcgEAAGAcSi4A\nAACMQ8kFAACAcSxOp9Pp7xAtZeDAgfrmm28UExPj7ygAAAC4gJiYGK1du/aSvkewj7K0CiEhIWpD\nnR4AAKDNalNncgEAANA2sCYXAAAAxqHkAgAAwDiUXAAAABiHkgsAAADjUHIBAABgHEouAAAAjEPJ\nBQAAgHEouQAAADAOJRcAAADGMbbkfvbZZ5oyZYoGDx6stLQ0rVq1yvVccXGxsrKyNHjwYI0ePVq7\nd+/2Y1K0dk0da4sWLVJcXJzi4+Nd/33ppZf8mBYmmDZtmubOnet6bLfbNWnSJCUlJWnixIkqKCjw\nYzqY5rvH2913393gc+2dd97xY0K0Ztu3b29wPN17772SLv2zzciS63Q6NW3aNEVHR2vDhg1asGCB\nnnvuOW3ZskWS9Ktf/Updu3bV66+/rrFjx2rGjBn6/PPP/ZwardH3HWtFRUWaNWuW3nvvPe3evVvv\nvfeeJkyY4OfUaM22bNmiXbt2uR7X1NRo2rRpSklJ0fr165WUlKTp06ertrbWjylhiu8eb9K3n2tP\nPvmk2+faD3/4Qz8lRGt37NgxpaWlaffu3a7j6dFHH/XJZ5uRJbesrEwJCQnKzs5Wz5499eMf/1ip\nqanat2+f/v73v6u4uFgPP/ywevfurWnTpikpKUnr1q3zd2y0Qk0da5JUWFiohIQERUVFuX6Fhob6\nOTVaq8rKSj3++OMaOHCga2zLli0KCwvT7Nmz1bt3b82fP182m01bt271Y1KYoLHjzeFwqLi4WImJ\niW6fa+3atfNjUrRmhYWF6tevnyIjI13HU4cOHXzy2WZkye3SpYueeuoptW/fXpK0b98+7d27V9dc\nc43++c9/qn///m5FIzk5WQcPHvRXXLRijR1re/bs0dChQ/XVV1+ptLRUV1xxhX9DwhhLly7VuHHj\n1KdPH9fYoUOHlJyc7DZvyJAhOnDgQEvHg2EaO96KiopksVjUvXt3PyaDSQoLC9WrV68G4774bDOy\n5P67tLQ0TZ48WUlJSRo5cqROnTqlrl27us2JiopSaWmpnxLCFOePtcGDB2vkyJEqLCyUxWLRc889\np+uuu07jxo3TG2+84e+YaKXy8vK0b98+3XPPPW7jX3zxBZ9p8LkLHW9FRUXq0KGD7r//fg0fPlwT\nJ05ssJwB8Mbx48f17rvv6oYbbtBPfvITPfXUUzp79qxPPtuCfR020DzzzDMqKyvTggULlJOTo5qa\nGoWEhLjNCQkJkcPh8FNCmOL8sZadna1HH31UiYmJCgoKUp8+fXT77bfrH//4h/77v/9bHTp00IgR\nI/wdF62Iw+HQggULlJ2d3eDzq7a2ls80+FRTx1tRUZHq6up07bXXatq0aXr77bd1991369VXX1X/\n/v39lBit1cmTJ1VbW6vQ0FA9/fTTKi4udq3H9cVnm/El9/z/dHPmzNGsWbM0YcIEVVVVuc1xOByy\nWq3+iAeDnD/W5s6dq9mzZ+uBBx5QWlqaIiIiJElXXnmlPvnkE7388suUXHjlmWeeUWJiYqMX94SG\nhjb40OczDZeiqeNtxowZyszMVHh4uCTpqquu0ocffqhXXnlFDz/8cEtHRSsXGxurDz74wPXnZFxc\nnOrr6zV79mwNHTr0kj/bjCy55eXlOnDggFuR6Nu3r86ePasuXbqosLDQbX5ZWZm6dOnS0jFhgKaO\nterqanXq1Mltfu/evfXBBx+0dEy0cn/5y19UXl6uwYMHS5LOnj0rSXrrrbc0evRonTp1ym0+n2m4\nFE0db/v373cV3PP69OnT4M9VwFPnC+55ffr0UV1dnaKjoy/5s83INbnFxcWaOXOm229Ofn6+oqKi\nlJycrIKCAre/Hezbt09JSUn+iIpW7kLHWmRkpF588UVlZWW5zT98+HCjC+yBpqxdu1abNm3Sxo0b\ntXHjRqWlpSktLU0bNmzQoEGDGlyIceDAAT7TcNGaOt7mzp2r+fPnu80/cuQIn2u4KO+9956GDh2q\nuro615jdblfnzp119dVXa//+/W7zvf1sM7LkDhgwQImJiZo7d64KCwv1zjvv6IknntDdd9+tlJQU\nxcTEaM6cOTp27JhWrlyp/Px89i7FRWnqWLv++uu1Z88e5ebm6sSJE/rTn/6kjRs3aurUqf6OjVYm\nJiZGPXr0cP2y2Wyy2Wzq0aOHbrjhBn355ZfKyclRYWGhFi1apK+//lo33nijv2OjlWrqeEtPT9fG\njRv1xhtv6LPPPtOyZcu0f/9+3X777f6OjVZo8ODBCgsL0/z583X8+HG98847evzxx3XXXXdp5MiR\nl/zZZnE6nc5mzO83p06d0iOPPKK8vDyFhYVp8uTJmjZtmiTpxIkTmjdvng4dOqSePXtq/vz5GjZs\nmJ8To7Vq6ljbuXOnnn76aX366ae6/PLLdd9997EeF5fs/N2nFi9eLOnbfz3Izs5WUVGRrrrqKi1c\nuFBxcXH+jAiDfPd4W7dunZ5//nl9/vnn6tu3r+bNm9dgqyfAU4WFhcrJydHBgwdls9n0i1/8Qr/6\n1a8kXfpnm7ElFwAAAG2XkcsVAAAA0LZRcgEAAGAcSi4AAACMQ8kFAACAcSi5AAAAMA4lFwAAAMah\n5AIAAMA4lFwAAAAYh5ILAAAA41ByAaCZpKWladmyZRf9+mPHjumdd95xPY6Li9Mbb7whSfrmm2/0\nxz/+8VIjun1PADAJJRcAAtT06dOVn5/verx7927ddNNNkqTNmzdr6dKl/ooGAAEv2N8BAACNczqd\nbo+joqJcX9fX17d0HABoVTiTCwB+4HA4tHTpUqWnpysxMVFDhw7Vb37zG50+fVrSt0sdSkpKtGzZ\nMt1xxx2S/rW04M9//rPmzZsnp9Op+Ph47dmzR8uWLVNaWprbezzzzDNuY6Wlpbr77rs1ZMgQ/ed/\n/qc2b97cINdf//pX/exnP9OgQYM0cuRIPf3003I4HM34OwEAzYOSCwB+8Pjjj2v79u1aunSp3n77\nbS1dulR5eXlavny5JGndunXq1q2b7rzzTv3hD39we+2oUaM0b948WSwW7d69W0lJSZIki8XiNs9i\nsbjGzp07pylTpqiyslJ/+tOf9PTTT2vVqlVur9m1a5fuu+8+3XLLLdqyZYsWLFigrVu36oEHHmjO\n3woAaBYsVwAAPxg4cKB++tOfKjk5WZIUExOjH/3oR/r4448lSZGRkQoKClL79u0VHh7u9tqQkBDX\nWGRkpEfv9/7776uwsFBvv/22unfvLklavHixMjIyXHNWrFihm2++WRMnTpQkde/eXQsWLFBmZqZm\nz56t2NjYS/uhAaAFUXIBwA/GjBmjvLw8Pfnkk/rkk09UVFSk48eP6+qrr26W9zt69KgiIiJcBVf6\ndvlDWFiY67Hdbld+fr5effVVt9cGBQWpsLCQkgugVaHkAoAfPPTQQ9q2bZvGjx+v9PR03XPPPVq1\napVKS0t99h7ffPPN984JDv7XHwP19fWaOnWqxo8f32Bely5dfJYLAFoCJRcAWtiZM2f06quv6ve/\n/71++tOfusYLCwtls9lcj7+7xvbfffe5du3aqbq62m3sk08+cX0dHx+vqqoqFRYWqk+fPq7nv/zy\nS9ecfv366fjx4+rRo4dr7IMPPtCaNWu0cOFCWa1W735QAPAjSi4ANKNPP/1U7777rttYaGioIiIi\ntH37diUkJKimpkZr166V3W53XUQmSe3bt9enn36q8vJyt+3Dzj8nSQUFBerbt6+SkpJUWVmp//3f\n/9UNN9ygd999V++++646deokSRo2bJgGDhyo2bNnKzs7W0FBQXr00Ud12WWXub7nXXfdpfvuu09/\n+MMfNGrUKJWUlGj+/Pnq2bNng/cHgEBncX53I0YAgE+c3wbsu2JjY7Vo0SItXrxYn332mTp27Kih\nQ4eqb9++WrlypXbv3q3Q0FC99tpreuyxx3T55ZfrjTfeUHx8vOtisaqqKt111106fPiwHn/8cd1w\nww169tln9ac//UnV1dW69tprlZycrBdffFE7duyQJFVWVuqRRx7R3/72N1mtVk2fPl0rVqzQrFmz\nXBegvfXWW1qxYoWOHTumjh07Kj09XbNmzVKHDh1a9PcOAC4VJRcAAADGYZ9cAAAAGIeSCwAAAONQ\ncgEAAGAcSi4AAACMQ8kFAACAcSi5AAAAMA4lFwAAAMah5AIAAMA4lFwAAAAYh5ILAAAA41ByAQAA\nYJz/Bw4QhS8UCMlvAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "df.sort_values('Latitude').plot('Latitude', 'p, Mpi100', ls='', marker='o')\n", "plt.ylabel('Mpi100 proportion')\n", "plt.legend().set_visible(False)\n", "plt.xlim(30, 50)\n", "plt.ylim(-0.1, 1.1)\n", "sns.despine()" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
AlleleLatitudeLocation
0048.1Port Townsend, WA
1048.1Port Townsend, WA
2048.1Port Townsend, WA
3048.1Port Townsend, WA
4048.1Port Townsend, WA
\n", "
" ], "text/plain": [ " Allele Latitude Location\n", "0 0 48.1 Port Townsend, WA\n", "1 0 48.1 Port Townsend, WA\n", "2 0 48.1 Port Townsend, WA\n", "3 0 48.1 Port Townsend, WA\n", "4 0 48.1 Port Townsend, WA" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rows = []\n", "for i, row in df.iterrows():\n", " for _ in range(row['Mpi90']):\n", " rows.append({'Location':row['Location'], 'Latitude': row['Latitude'], 'Allele': 0})\n", " for _ in range(row['Mpi100']):\n", " rows.append({'Location':row['Location'], 'Latitude': row['Latitude'], 'Allele': 1})\n", "raw_df = pd.DataFrame(rows)\n", "raw_df.head()" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-7.59881547] [[ 0.17754715]]\n" ] } ], "source": [ "clf = sklearn.linear_model.LogisticRegression(C=1e12, random_state=0)\n", "clf.fit(raw_df['Latitude'].reshape(-1, 1), raw_df['Allele'])\n", "print(clf.intercept_, clf.coef_)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which is very close to McDonald's intercept of -7.6469 and slope of 0.1786." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxAAAAIYCAYAAAD93QokAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xt0VeWd//HPSUIuJCEJud8DgiSAJBCSQABRAihQBemq\nqyrW2jK0f2jVolXHzk8dR5fXWWN1uqqj45rCrKmjjbUVlUsQlAQIt4RLEiQYQkLukAuB3M/+/cGA\npIGTHUhOzjl5v9ZyrZ69v/vsb6yQ8znPfp7HYhiGIQAAAAAwwW24GwAAAADgPAgQAAAAAEwjQAAA\nAAAwjQABAAAAwDQCBAAAAADTCBAAAAAATCNAAAAAADCNAAEAAADANAIEAAAAANOcKkB0dnbqjjvu\n0J49e65as23bNq1YsULTp0/X8uXLtXXrVjt2CAAAALg2pwkQnZ2d+vWvf63S0tKr1hw9elQPP/yw\nfvSjH+mvf/2r7r77bv3qV7/S0aNH7dgpAAAA4LqcIkAcP35cd999tyorK23WffbZZ5o9e7buu+8+\nxcbG6r777lNGRoa++OILO3UKAAAAuDaP4W7AjPz8fM2ePVuPPvqokpOTr1p31113qaurq8/x1tbW\noWwPAAAAGDGcIkDcc889purGjx/f6/WxY8e0a9cu3XvvvUPRFgAAADDiOMUjTNfizJkzevjhh5Wa\nmqqsrKzhbgcAAABwCS4ZIBoaGvTAAw/IYrHozTffNH3dqlWrtGrVqiHsDAAAAHBuTvEI00DU1tbq\nJz/5idzd3bVu3ToFBQWZvra6unoIOwMAAACcn0uNQLS1tWn16tUaNWqU1q9fr5CQkOFuCQAAAHAp\nTj8C0dDQIH9/f3l5eekPf/iDKisr9cc//lFWq1UNDQ2SJG9vb/n5+Q1zpwAAAIDzc7oRCIvF0uv1\n3LlzL+3zsGnTJrW3t+vuu+/WvHnzLv3z4osvDkerAAAAgMuxGIZhDHcTjuLiak05OTnD3AkAAADg\nmJxuBAIAAADA8CFAAAAAADCNAAEAAADANAIEAAAAANMIEAAAAABMI0AAAAAAMI0AAQAAAMA0AgQA\nAABGvKeeekqJiYn6/e9/3+dca2urpk6demnPsP48/fTTevrpp/sc37t3rxYuXNjn+GeffaZFixZp\n+vTpeuihh9TY2Njr/Ouvv67Zs2crIyNDr732msmfaOgQIAAAADAkDMNQUd0x5Z7co6K6Y3Lk/Yst\nFotGjRqlrVu39jm3fft29fT0mH6vZ555Rs8880yvY0ePHtWjjz7a59/BwYMH9dvf/lYPP/ywPvzw\nQzU3N/cKH//5n/+pDRs26Pe//73eeust/e1vf9MHH3wwwJ9ucHkM690BAADgkvIrC7SuMFu1rfWX\njoX7her+5JVKj0kZxs6uLjU1Vfn5+aqrq1NYWNil41u2bFFKSorq6upMvY+fn1+v13/605/06quv\nKi4uTmfPnu117r//+7+1ZMkS3XnnnZKk1157TbfeeqtOnTql6OhorVu3To888oimT58uSXr88cf1\n5ptv6sEHH7yeH/W6MAIBAACAQZVfWaA38t7tFR4kqba1Xm/kvav8yoJBv+epU6eUmJiozz77TDff\nfLPS09P14osvymq1mn6PyMhIJSUl6auvvrp0rKurSzt27NCCBQsuHcvPz9f8+fO1bt06ZWRkaO7c\nufrDH/5w6fzfP8K0Y8cOvfrqq3rggQf63LOgoEBpaWmXXkdERCgyMlKFhYWqq6tTdXW1Zs6ceel8\namqqqqqq1NDQYPrnGmwECAAAAAwawzC0rjD7qo8rGYah9TbOX69///d/15tvvqm3335bmzZt0u9+\n97sBXZ+VlaWcnJxLr/Py8jRx4kSFhIT0qjt9+rQ+/fRT/dd//Zeef/55vffee/roo4+u+J5vv/32\nFec+SFJ9fX2v0Q5JCgkJUU1Njerr62WxWHqdDwkJkWEYqqmpGdDPNZgIEAAAABg0xfWlfUYe/l5N\na71KGkqH5P6/+c1vNH36dKWnp+uRRx656of6q1mwYIF27dqltrY2SVJOTo4WLVrUp66np0cvvfSS\nEhMTlZWVpQceeEAffvjhgPttb2+Xp6dnr2Oenp7q7Oy81MPl5y/+787OzgHfa7AQIAAAADBoGtub\nzNW1NQ/6vS0Wy6W5ApI0depUnTlzps+qRrYkJiYqNDRUubm5MgxDW7duveLogY+Pj2688cZe9zp+\n/PiAe/by8uoTBjo7O+Xt7S0vL69Lry8/d/H+w4UAAQAAgEET5B1ors4nYEju7+Hx/RpBF+c/uLkN\n7CPvggULlJOTo4KCAgUHBys2NrZPzahRo3q97unpGfB9JCksLKzPfIaGhgaFhYUpPDxchmH0On/x\nsabQ0NAB32uwECAAAAAwaJJCJyjcz/aH2wi/UCWGTBj0exuGoZKSkkuvDx06pLCwMAUEDCysZGVl\nafv27dq8efMVH1+SpJaWFlVVVfW616RJkwbcc0pKivbt23fpdXV1tWpqapSSkqKwsDBFRUX1Or93\n715FRkb2mZNhTwQIAAAADBqLxaL7k1fKYrFc9fwqG+ev14svvqjDhw8rLy9Pv/vd73TfffdJujAa\n0dDQoK6urn7fIy0tTT09Pfrf//3fq05+NgxD//RP/6Rjx45p48aNWr9+vVatWjXgfu+55x59+umn\n+vjjj1VSUqInn3xSt956q6KioiRJP/7xj/X6668rPz9fu3fv1r/+679ecTUne2IfCAAAAAyq9JgU\nrc1co/WF2aq5bEJ1hF+oVg3xPhBLlizRL37xCxmGoXvvvVdr1qyRdOGb/aysLK1bt67XsqlX4u7u\nrptvvlkHDhxQYmLiFWssFovmzZune++9V76+vlq7dq2WLl064H5TUlL0z//8z3rzzTfV3NysuXPn\n6oUXXrh0fvXq1WpsbNTDDz8sNzc33X333cMeICyGI28JaGcXtye/fOkuAAAAXBvDMFRcX6qm9mYF\n+QQoMWTCkI08nDp1SgsXLlROTs6lb+//3ltvvaX58+dr2rRp13Wv/Px8PfDAAyouLr6u93FWjEAA\nAABgSFgsFk0Om2i3+9n6Xry1tVW7du3SL3/5S7v146oIEAAAAHAJtkY3/Pz89Mc//lHu7u527Mg1\n8QjTZXiECQAAALCNVZgAAAAAmEaAAAAAAGAaAQIAAACAaQQIAAAAAKYRIAAAAACYRoAAAAAAYBoB\nAgAAAIBpBAgAAAAAphEgAAAAAJhGgAAAAABgGgECAAAAgGkECAAAAACmESAAAAAAmEaAAAAAAGAa\nAQIAAACAaQQIAAAAAKYRIAAAAACYRoAAAAAAYBoBAgAAAIBpBAgAAAAAphEgAAAAAJhGgAAAAABg\nGgECAAAAgGkECAAAAACmESAAAAAAmEaAAAAAAGAaAQIAAACAaQQIAAAAAKYRIAAAAACYRoAAAAAA\nYBoBAgAAAIBpBAgAAAAApjlVgOjs7NQdd9yhPXv2XLWmqKhId999t1JSUvSjH/1IR44csWOHAAAA\ngGtzmgDR2dmpX//61yotLb1qTVtbm9asWaO0tDRlZ2crJSVFv/jFL9Te3m7HTgEAAADX5RQB4vjx\n47r77rtVWVlps27Dhg3y8fHRE088ofHjx+uZZ56Rr6+vvvzySzt1CgAAALg2pwgQ+fn5mj17tj78\n8EMZhnHVuoMHDyo1NbXXsRkzZujAgQND3SIAAAAwIjhFgLjnnnv05JNPysvLy2ZdXV2dwsLCeh0L\nDg5WbW3tULYHYAg1NjZqy5Yt+uijj7RlyxY1NjYOd0tX1NLSol27dmnTpk3atWuXWlpahrslwCFY\nLJY+/1x011139Tp+1113Dfg9HElAQECvHgMCAoa7pRHviSee6PX/yRNPPDHcLbkEpwgQZrW3t8vT\n07PXMU9PT3V2dg5TRwCuR2Njoz799FOdO3dOnp6eOnfunD799FOHCxEtLS3KycnRuXPn5O7urnPn\nziknJ4cQgRHvah/0L4aFv/zlL72O/+Uvf+kTImy9hyMJCAjo82e+paWFEDGMnnjiCb3++uu9jr3+\n+uuEiEHgUgHCy8urT1jo7OyUt7f3MHUE4Hrs27dPQUFBvY4FBQVp3759w9TRlRUVFWnMmDG9jo0Z\nM0ZFRUXD1BHg+P4+PPR33NFd7QsDvkgYPn8fHvo7DvNcKkCEh4ervr6+17GGhgaFhoYOU0cArsfV\nRhoccQRiIMcBAHBmLhUgkpOT+0yYPnDggFJSUoapIwDX4+9HH/o7Plz+fvShv+MAADgzpw8QDQ0N\n6ujokCTddtttOnv2rF566SUdP35c//Iv/6Lz589ryZIlw9wlgGuRmpraZ7ShsbGxz2prw23y5MlX\nfPZ58uTJw9QR4PhWrFgxoOOOji8SHM/jjz8+oOMwz+kCxN9Pmpo7d66++OILSZKfn5/+8Ic/aO/e\nvfrhD3+oQ4cO6T/+4z+YAwE4qaCgIC1fvly+vr7q7OyUr6+vli9f7pAjEFlZWfL19VVPT498fX2V\nlZXFBweMeFdbet0wDH3yySd9wsKKFSv0ySefmH4PR9Lc3HzFuVDNzc3D1BFee+21PmHh8ccf12uv\nvTZMHbkOi+FofwKHUVZWliQpJydnmDsBAAAAHJPTjUAAAAAAGD4ECAAAAACmESAAAAAAmEaAAAAA\nAGAaAQIAAACAaQQIAAAAAKYRIAAAAACYRoAAAAAAYBoBAgAAAIBpBAgAAAAAphEgAAAAAJhGgAAA\nAABgGgECAAAAgGkECAAAAACmESAAAAAAmEaAAAAAAGAaAQIAAACAaQQIAAAAAKYRIAAAAACYRoAA\nAAAAYBoBAgAAAIBpBAgAAAAAphEgAAAAAJhGgAAAAABgGgECAAAAgGkECAAAAACmESAAAAAAmEaA\nAAAAAGAaAQIAAACAaQQIAAAAAKYRIAAAAACYRoAAAAAAYBoBAgAAAIBpBAgAAAAAphEgAAAAAJhG\ngAAAAABgGgECAAAAgGkECAAAAACmESAAAAAAmEaAAAAAAGAaAQIAAACAaQQIAAAAAKYRIAAAAACY\nRoAAAAAAYBoBAgAAAIBpBAgAAAAAphEgAAAAAJhGgAAAAABgGgECAAAAgGkECAAAAACmESAAAAAA\nmEaAAAAAAGAaAQIAAACAaQQIAAAAAKYRIAAAAACY5hQBorOzU//4j/+otLQ0zZs3Tx988MFVazdv\n3qxly5Zp+vTpuu+++1RUVGTHTgEAAADX5hQB4pVXXlFRUZHWrVunZ599Vm+//bY2bdrUp660tFSP\nP/64fvGLX+ivf/2rEhMTtWbNGnV0dAxD1wAAAIDrcfgA0dbWpo8//li//e1vlZiYqIULF2r16tVa\nv359n9odO3Zo4sSJuvPOOxUbG6tf//rXamhoUGlp6TB0DgAAALgehw8QJSUl6unpUUpKyqVjqamp\nOnjwYJ/awMBAlZaWav/+/TIMQ3/+85/l7++vuLg4e7YMAAAAuCyP4W6gP/X19QoMDJSHx/etBgcH\nq6OjQ42NjQoKCrp0fOnSpdq6davuvfdeubu7y83NTe+++678/f2Ho3UAAADA5Tj8CERbW5s8PT17\nHbv4urOzs9fxpqYmNTQ06Nlnn9VHH32kFStW6KmnntKZM2fs1i8AAADgyhw+QHh5efUJChdf+/j4\n9Dr++uuva9KkSbrnnns0efJk/fM//7N8fHyUnZ1tt34BAAAAV+bwASI8PFxNTU2yWq2XjjU0NMjb\n21tjxozpVXvkyBElJiZeem2xWJSYmKiqqiq79QsAAAC4MocPEElJSfLw8FBBQcGlY3v37tXUqVP7\n1IaFhfVZcamsrEwxMTFD3icAAAAwEjh8gPD29tby5cv17LPP6tChQ9qyZYs++OADPfDAA5IujEZc\n3OfhRz/6kT766CN9+umnOnnypF5//XVVV1drxYoVw/kjAAAAAC7D4VdhkqSnn35azz//vB544AH5\n+/vrkUce0cKFCyVJc+fO1csvv6wVK1Zo6dKlamtr0zvvvKPa2lolJSXpj3/8o8aOHTvMPwEAAABg\nf13dPTry3WntK6lTfWObnnog7brf02IYhjEIvbmErKwsSVJOTs4wdwIAAABcm9PNbdpTVKu9xbUq\nPFav9s6eS+f+9sby635/pxiBAAAAAHBlVquh0som7SmqVX5Rjb471Tyk9yNAAAAAAE6mraNbBd/W\na09RjfYW16rxbIfd7k2AAAAAAJxA3Znz2lNUo/ziWh0qbVBXt7X/iyRZLNLE2EDNTAwflD4IEAAA\nAIAD6rEaOnayUflFNdpTVKsT1S2mr/X19tCMxHDNTArXjElhCvT3GrS+CBAAAACAg2jr6Nb+o3XK\nP3Lh0aSWc52mr40O9VXa5AilT45Q0rix8nAfmh0bCBAAAADAMGo82678I7XafaRaBd/Wm340yd3N\noinjg5U2OUJpk8MVHeo3xJ1eQIAAAAAA7KyqvlW7Dldr1+EalZSfkdmNFfxHe2pmUpjSJkdo+qQw\n+fmMGtpGr4AAAQAAAJdkGIaK60vV2N6kIO9AJYVOkMViGZZeLi61ejE0VNSeNX1tXIS/0v9vlGFS\n/Fi5uw3Pz3ARAQIAAAAuJ7+yQOsKs1XbWn/pWLhfqO5PXqn0mBS79NDVbdWh0gbtOlKt3YdrdKal\n3dR17m4W3XRDiNKnXAgNEcG+Q9zpwBAgAAAA4FLyKwv0Rt67Mv7uuaDa1nq9kfeu1mauGbIQcb69\nS/uK67TrcLX2ltTqfHu3qet8vNw1IzFcs6ZGamZS+LA8mmQWAQIAAAAuwzAMrSvM7hMeLj+/vjBb\nadHJg/Y4U+v5TuUX1Si3sFoHvq0zPQk60N9LGVMiNGtqpKZNCJHnKPdB6WeoESAAAADgMorrS3s9\ntnQlNa31KmkoVVLoxGu+T3Nrh3YdrlbewWoVHqtXj9XcLOjoUF/NmhqpWVMjdWNckNyGeT7DtSBA\nAAAAwGU0tjeZq2trHvB7n25u085DF0LDke8aZDIzaFJ80KWRhthw/wHf19EQIAAAAOAygrwDzdX5\nBJiqqz1zXnkHq5R3sEol5Y2mrnF3syh5Yqhm3RSpjCkRGjvG29R1zoIAAQAAAJeRFDpB4X6hNh9j\nivALVWLIhKueP1Xfeik0lFaaG6kY5eGmGZPClDktUumTI+Q32nPAvTsLAgQAAABchsVi0f3JK6+4\nCtPF86uSV/aZQF3V0KodBVXaUXhKZVUtpu7l5emumUnhmnNTlFKTwjTa23FXThpMBAgAAAC4lPSY\nFK3NXKP1hdmquWwkIsIvVKsu2wei5vQ55RZW6ZvCUzpucqRhtLeH0qdEKPOmKM1IDJOXk6ycNJgI\nEAAAAHA56TEpSotOVnF9qZramxXkE6DEkAlqaGrXJ9tKtaPwlL49aW7Ctf9oT82aGqHMaVFKnhii\nUR4jLzRcjgABAAAAl2SxWDQ5bKJON7cpt7BK/1nwjemJ0EH+Xpp1U6Tm3BSlqTcEy93dbYi7dR4E\nCAAAALicxpZ25R2s0jeFVSoqO62r7CvXS4CfpzKnRWleSrQmjwuWuxPu0WAPBAgAAAC4hObWDuUd\nrNKOwiodPm5unwb/0Z7KnBapecnRjDSYRIAAAACA0zrX1qWdh6q1/UClDpY2yGoiNfj6jFLmTZGa\nmxytaRND5EFoGBACBAAAAJxKZ1eP9hbXavuBSu0pqlVXt7Xfa0Z7e2jW1EjNTY5Syo1hGuVBaLhW\nBAgAAAA4vB6rocOlDdq2v1J5h6p0vr2732t8vNyVPjlS81KiNH1SmDxH4JKrQ4EAAQAAAIdkGIaO\nVTRp+4FKfXPglBrPdvR7jZenu9KSwjUvJVqpSeEjcp+GoUaAAAAAgEOprDur7ftPafuBSlU3nOu3\n3sPdotTEcN08PVrpkyPk7cVH3KHEv10AAAAMu9PNbfqm4JS2769UqYldoS0Waer4EM2fEa3MaVHy\nH+1phy4hESAAAAAwTFrbupRbWKWvD1Tq0PEGU3s1jI8O0PzpMZqXEq3QIJ+hbxJ9ECAAAABgN13d\nVu0vqdVX+yq1+0iNunv6X0EpMthXN8+I1vzpMYoN97dDl7CFAAEAAIAhdXEy9Fd7K/R1wSm1nOvs\n95pAfy/dnBKt+TNiNDE2UBYLu0I7CgIEAAAAhkRd43lt21epr/ZVqLKutd96Hy8PZU6L1PzpMZo2\nIYRdoR0UAQIAAACD5nz7hXkNX+27MK+hPx7ubkqbHK75M2I0k2VXnQIBAgAAANelp8eqA9/W66t9\nFdp1qFqdJnaGTkoYq1tnxmpecpT8WEHJqRAgAAAAMGCGYaisqkVb91Zo+4FKNZnY5C0ieLRuTY3V\nLakxigrxs0OXGAoECAAAAJh2urlN2/dXauveCpXXnO233tdnlOalROvW1BglJYxlMrQLIEAAAADA\npq7uHu0+UqOcPRXaX1Iraz/7Nbi7WTQzKVy3zoxVWlK4PJnX4FIIEAAAAOjDMAwdr2xWzp6T2n6g\nUmfPd/V7zcTYQC2YGat5KdEK8POyQ5cYDgQIAAAAXNJ0tkPb9lcoZ0+FTlS39FsfGuSjW2bE6NbU\nWDZ5GyEIEAAAACNcd49Ve4pqlbPnpPYW16qnn2eUfLzcNWdatBbMjNWU8cFyc2New0hCgAAAABih\nyqqatWXPSW3fX6nm1v53h77phhBlpcUqc1qUfLz4GDlS8f88AADACNJyrlPb91cqZ+9JHa9s7rc+\nLMhHC2bGKSstVhHBvnboEI6OAAEAAODienqs2n+0Tlv2nFT+kRp199h+RMlzlLvmTItUVlqcbroh\nhEeU0AsBAgAAwEXVnD6nzfkntSX/pM60tPdbn5QwVllpcZqXEqXR3qPs0CGcEQECAADAhXR29Wjn\noWpt2l2ug6UN/dYHB3hrwcxYZaXFKTqU3aHRPwIEAACACyiratam3eXatq9SrW2292wY5eGm2VMj\nlZUep+SJoXLnESUMAAECAADASZ1v79L2A6e0aXe5Siua+q2fEBuoRelxujklWn6jPe3QIVwRAQIA\nAAwZwzBUXF+qxvYmBXkHKil0giwWvu2+HoZhqKjsjDbtLlfuwSp1dPbYrPfzGaVbZ8ZqUXqcxkUF\n2KlLuDICBAAAGBL5lQVaV5it2tb6S8fC/UJ1f/JKpcekDGNnzqnxbLu+2luhTbtP6lR9a7/1yRND\ntDgjXrOmRspzlLsdOsRIQYAAAACDLr+yQG/kvSvD6L1caG1rvd7Ie1drM9cQIkzosRo6cLROm3aX\nK/9ITb87RI8d461F6XFamB7Hng0YMgQIAAAwqAzD0LrC7D7h4fLz6wuzlRadzONMV1F75rw255cr\nJ/+kGpptL7/q7mZR+pQILUqP04xJYXJ3d7NTlxipCBAAAGBQFdeX9nps6UpqWutV0lCqpNCJdurK\n8XX3WJV/pEZf7jyhgmP1ukr+uiQ61FeL0uO1YGasgsZ426VHQCJAAACAQdbY3v9qQJLU2NY8xJ04\nh5rT57Rpd7m25J9U49kOm7Weo9w1NzlKizPiNXncWEZwMCwIEAAAYFAFeQeaq/MZuSsCXRxt2Lir\nXAe+ret3tGFCTIAWZ8Tr5ukx8vVhh2gMLwIEAAAYVEmhExTuF2rzMaYIv1AlhkywY1eOofbMeW3c\ndcLUaIOvzyjdOiNGizLiNT565IYtOB6nCBCdnZ167rnntHnzZnl7e+tnP/uZHnzwwSvWHj16VM8/\n/7yOHDmi+Ph4PfPMM8rIyLBzxwAAjFwWi0X3J6+84ipMF8+vSl45Yh6/6e6xak9Rjb7caW60ISlh\nrG6fHa85ydHyYvlVOCCnCBCvvPKKioqKtG7dOlVWVurJJ59UdHS0Fi9e3KuutbVVP//5z5WVlaVX\nXnlFf/nLX/TQQw9p48aNGjt27DB1DwDAyJMek6K1mWu0vjBbNZeNRET4hWrVCNkHYqCjDQtmxuq2\njHjFR46xU4fAtXH4ANHW1qaPP/5Y77//vhITE5WYmKjVq1dr/fr1fQJEdna2fH199fzzz0uSHn74\nYX399dc6fPiwbr755uFoHwCAESs9JkVp0ckqri9VU3uzgnwClBji2jtRX+toQ+a0KHl7OvzHMkCS\nEwSIkpIS9fT0KCXl+28qUlNT9c477/Sp3bNnjxYsWNDr2EcffTTkPQIAgCuzWCyaHOb6S7XWnjn/\nfyspletMSz+jDd4eunVmrG6flcBoA5ySwweI+vp6BQYGysPj+1aDg4PV0dGhxsZGBQUFXTpeUVGh\nm266Sf/v//0/bd26VTExMfrNb36jGTNmDEfrAADAhV0abdhVrgNHzY023DYrXnOSGW2Ac3P4/3rb\n2trk6enZ69jF152dnb2Onz9/Xu+9955+8pOf6L333tNnn32mn//85/ryyy8VHh5ut54BAIDrqm9s\n08ZdJ7R5AKMNt81KUAKjDXARDh8gvLy8+gSFi699fHx6HXd3d1dSUpIeeughSVJiYqJyc3P16aef\nas2aNfZpGAAAuByr1dCBb+v0Rd4J7SmqkbWf0YbE+CDdNitBc1MYbYDrcfj/osPDw9XU1CSr1So3\nNzdJUkNDg7y9vTVmTO8kHxoaqvHjx/c6lpCQoOrqarv1CwAAXEdza4e25J/Ul7tOqOb0eZu1vt4e\nujU1VrfNZrQBrs3hA0RSUpI8PDxUUFBwaS7D3r17NXXq1D61KSkp2rNnT69j3333ne644w679AoA\nAJyfYRgqOdGoz3eWKbewSl3dVpv1jDZgpHH4/8q9vb21fPlyPfvss3rppZdUW1urDz74QC+//LKk\nC6MR/v7+8vLy0o9//GOtX79eb7/9tu6880598sknqqys1J133jnMPwUAAHB059u7tH1/pT7PO6ET\n1S02a3283HVLaqyWzE7QuKgLu0QbhqGiumNqbG9SkHegkkJde8lajFwW40pbRDqY9vZ2Pf/889q4\ncaP8/f21evVq3X///ZIuzHN4+eWXtWLFCknSgQMH9MILL+j48eO64YYb9Nvf/tb0KkxZWVmSpJyc\nnKH5QQAAgMM5Ud2iz/PKtG1fpdo6um3WJkSO0dLMBM2fEaPR3qMuHc+vLNC6wmzVXrZpXrhfqO4f\nIZvmYWRyiZSgAAAgAElEQVRxigBhLwQIAABGhq7uHuUerNbnuWUqPnHGZq2Hu5vmJkdpSWaCkhLG\n9hlVyK8s0Bt57+pKH6ksFovWZq4hRMClOPwjTAAAAIOl5vQ5fbnzhLbsOanm1k6btRHBo3X7rAQt\nTI9TgJ/XFWsMw9C6wuwrhoeL59cXZistOpnHmeAyCBAAAMCl9VgN7Suu1ed5Zdrfz4ZvbhYpbXKE\nlmQmaPqNYXJzs/2hv7i+tNdjS1dS01qvkoZSJYW6/o7cGBkIEAAAwCU1trRrU365Nu4qV31jm83a\nIH8vLc6I1+JZ8QoLGm3+Hu1N5uramk2/J+DoCBAAAMBlGIahw9+d1ue5Zdp5qFo9/ez4Nm1CiJZk\nJmjW1Eh5uLsN+H5B3oHm6nwCBvzegKMiQAAAAKd3rq1LW/dW6IudJ1RRe9Zmra+3hxakxWnJ7ATF\nhvtf132TQico3C/U5mNMEX6hSgyZcF33ARwJAQIAADitE9Ut2pBbpm37KtTe2WOzdkJMgJZkjtPN\nKdHy9hqcj0AWi0X3J6+0uQrTquSVTKCGSyFAAAAAp9LdY9XOQ9XakFumI9+dtlnr6eGmm6fHaElm\ngm6MCxqSftJjUrQ2c43WF2ar5rKRiAi/UK1iHwi4IAIEAABwCmda2rVx5wl9uatcZ1rabdZGh/pq\nSeY4Zc2Mld9ozyHvLT0mRWnRySquL1VTe7OCfAKUGMJO1HBNBAgAAOCwDMNQUdkZbcgtU97BKpuT\not3cLJo1NUJLM8dp2oQQu394t1gsmhzGUq1wfQQIAADgcNo7urVtf6U25JbpRHWLzdpAfy/dNite\nt89KUEigj506BEYuAgQAAHAYVfWt2pBXppz8kzrX3m2zNilhrJbNGafMaVEa5THwJVgBXJtrChDb\nt2/Xe++9p7KyMn344YfKzs5WXFycli9fPtj9AQAAF3dxp+gNuRd2irbFc5S7bpkRo2Vzxml8NHsr\nAMNhwAEiNzdXDz30kJYtW6bCwkJZrVZ1d3fr6aeflmEYWrFixVD0CQAAXEzLuU5t3l2uz3eeUN2Z\n8zZrI4N9tXROghamxdllUjSAqxtwgHjrrbe0du1a/fSnP9XGjRslSY899pj8/Pz0/vvvEyAAAIBN\nxyoatSG3TN8cOKXObutV6ywWKTUxXMvmjNOMSWFyc2NFI8ARDDhAHD16VK+++mqf47fffrvefvvt\nQWkKAAC4lq7uHn1TUKXPc8t09GSjzVo/n1FalBGvpZkJigj2tVOHAMwacIDw9/dXXV2d4uLieh0v\nLS1VQADPIgIAgO/VNZ7XlztPaNPucjW3dtqsvSEmQMsyx2ne9Gh5e7LOC+CoBvyn84477tBLL72k\nl156SRaLRefOndPXX3+tF154QUuXLh2KHgEAgBMxDEOFx+q1IbdM+UdqZGPrBnm4WzQ3OVrL5o7T\npLigK+7dYBiGiutL1djepCDvQCWFskEbMJwGHCAeffRR1dTUXJrrcNddd8kwDN1yyy167LHHBr1B\nAADgHM63dylnT4U+zytTZV2rzdqQAG/dnpmgxRnxCvL3vmpdfmWB1hVmq7a1/tKxcL9Q3Z+8Uukx\nKYPWOwDzLIZh2Phe4OrKy8tVXFwsq9WqG2+8URMmTBjs3uwuKytLkpSTkzPMnQAA4DzKa1q0IbdM\n2/ZVqK2jx2bttAkhWjZnnDKmRMjd3fbeDfmVBXoj711d6aOKxWLR2sw1hAhgGFzzA4bx8fGKj48f\nzF4AAICT6O6xavfhGm3ILdOh4w02a3283LVgZpyWZiYoLmKMqfc3DEPrCrOvGB4unl9fmK206GQe\nZwLszFSASExMNP2Hs7i4+LoaAgAAjquxpV0bd5fry50ndLq53WZtbLifls0Zr1tTYzTae9SA7lNc\nX9rrsaUrqWmtV0lDqZJCJw7ovQFcH1MB4uKEaQAAMPIYhqHiE2e0IbdMeQer1N1z9aef3dwsypgS\noR/MHaebbgi55s8Pje1N5uramq/p/QFcO1MBYuXKlUPdBwAAcDDtnd3avv+UPs8t03dVtj+oB/p5\nafGseN0+K0GhQT7Xfe8g70BzdT4sIQ/Y2zXNgdi+fbvef/99fffdd/rwww+VnZ2tuLg4LV++fLD7\nAwAAdlbV0Kov8k5oc/5JnWvrslmbGB+kZXPHa860SI3ycB+0HpJCJyjcL9TmY0wRfqFKDHH+RVwA\nZzPgAJGbm6uHHnpIy5YtU0FBgaxWq7q7u/X000/LMIxLy7sCAADnYbUa2ldSqw25ZdpXUmez1tPD\nTfNnxGjpnHGaEGNupGCgLBaL7k9eaXMVplXJK3nEGhgGAw4Qb731ltauXauf/vSn2rhxoyTpscce\nk5+fn95//30CBAAATuTs+U5t3n1SX+wsU83p8zZrI4JHa8nscVqUESf/0Z5D3lt6TIrWZq7R+sJs\n1Vw2EhHhF6pV7AMBDJsBB4ijR4/q1Vdf7XP89ttv19tvvz0oTQEAgKF1vLJJG3LLtH1/pTq7rVet\ns1ikGZPC9IO54zVjUpjc3Oz7jX96TIrSopNVXF+qpvZmBfkEKDGEnaiB4TTgAOHv76+6ujrFxcX1\nOl5aWqqAACYyAQDgqLq6rco9WKXPc8tUfOKMzVpfn1FalB6npZnjFBnia6cOr8xisWhyGEu1Ao5i\nwAHijjvu0EsvvXRpaddz587p66+/1gsvvKClS5cORY8AAOA6nG5u0xc7T2jjrnI1ne2wWTs+KkDL\n5o7TzdOj5e15zfvNAnBhA/6b4dFHH1VNTc2luQ533XWXDMPQLbfcoscee2zQGwQAAANnGIYOf3da\nG3aUaefhalmtV9+7wcPdosxpUVo2Z5ySEsba7fEgwzBUXF+qxvYmBXkHKimUR5MAZ2AxrrZHfD/K\ny8tVXFwsq9WqG2+8URMmOP8yallZWZKknJycYe4EAIBr09bRrW37KvRZbplO1py1WTt2jLeWZCbo\ntox4BY3xtlOHF+RXFmhdYXavZVrD/UJ1P5OjAYd3zQHCFREgAADOqrLurD7PO6GcPSd1vr3bZu3U\nG4L1gznjlTE1Qh7ubnbq8Hv5lQU2l2ddm7mGEAE4MFOPMCUmJpoeUiwuLr6uhgAAgDk9VkN7i2r0\nWW6ZCr69+oZrkuTl6a4FqbFaNmec4iPH2KnDvgzD0LrC7CuGh4vn1xdmKy06mceZAAdlKkBcnDAN\nAACGX3Nrhzbnn9QXeWWqa2yzWRsV4qtlc8ZpQVqc/HxG2anDqyuuL7W5u7Qk1bTWq6ShVEmhrLwE\nOCJTAWLlypVD3QcAAOjHsYpGfbajTN8UnFJXP3s3pCVFaNnccUqZGGr3vRtsaWxvMlfX1jzEnQC4\nVqYCxEA2iHvooYeuuRkAANBbV3ePdhRWacOOMh092Wiz1n/0KC3OiNftsxMUETy8ezdcTZB3oLk6\nH/aWAhyVqQCRnZ1t6s3a29sJEAAADIL6xjZ9sbNMm3aXq7m102btDTEB+sGc8Zo3PVpeo9zt1OG1\nSQqdoHC/UJuPMUX4hSoxxPlXdwRclakAsXXrVpvnjx07pj/96U/629/+NihNAQAwEhmGoYOlDdqQ\nW6bdh6tlY+sGebi7aW7Khb0bJsUFOc1cRYvFovuTV9pchWlV8kqn+XmAkeiat5js7OzUl19+qT/9\n6U86cOCALBaLFi5cOJi9AQAwIpxv79JXeyu0Ia9MFbWtNmtDAry1JHOcFmfEK9Dfy04dDq70mBSt\nzVyj9YXZqrlsJCLCL1Sr2AcCcHgDDhAnTpzQhx9+qE8++URNTU2yWCxauXKlfvnLXyo2NnYoegQA\nwCVV1J7Vhtwybd1bobYO23s3TJsQomVzxiljSoTch2HvhsGWHpOitOhkFdeXqqm9WUE+AUoMYSdq\nwBmYChDd3d3avHmzPvzwQ+3evVvu7u6aO3euli1bpqeffloPPvgg4QEAABN6eqzKL6rRZzvKdLC0\nwWatj5e7bk2N1dI54xQfMXx7NwwVi8WiyWEs1Qo4G1MB4pZbblFra6syMjL0wgsvaNGiRQoIuLA6\nwlNPPTWkDQIA4AqaznZo0+5yfbHzhBqabO/dEB3qpx/MHacFM2M12nv4924AgMuZChBnz55VcHCw\noqKiFBgYKB8fn6HuCwAAl/DtyUZ9tuM7fVNQpe6eq+/d4GaR0qdE6AdzxmvaxBAe5QHgsEwFiNzc\nXH3++ef685//rP/5n/+Rr6+vsrKytHTpUv6CAwDg77R3dmtHwSltyDuh0grbG6eN8fXU4ox4LZmd\noLCxo+3UIQBcO4txpTXUbDh+/Lg+/vhj/e1vf1NDQ4MsFot++MMfavXq1UpISBiiNu0jKytLkpST\nkzPMnQAAnNGp+lZ9kXdCOXtOqrWty2btjXGBWjZnnOYmR8vTwfduAIDLDThAXNTT06Nt27bpk08+\n0bZt22S1WpWZman33ntvsHu0GwIEAGCgLkyKrtXneWUq+Pbqm6NJF/ZuuHl6tJbNGacb44Ls1CEA\nDK5r3gfC3d1dWVlZysrK0pkzZ/Tpp5+a3rEaAABn19jSrk27y/XlzhNqaG63WRsa5KMlsxO0OCNe\nAX7OuXcDAFx0zSMQrogRCACALYZh6Mh3p/V53gnlHaxSj62toiXNmBSmpZkJmjk5Qu5uzBkE4Bqu\neQQCAICR4nx7l77aV6nP88p0suaszVo/n1FamB6nJZkJigrxs1OHAGA/BAgAAK7iRHWLPs8r07Z9\nFWrr6LFZe2NcoJbMHqd506PlxaRoAC6MAAEAwGW6uq3aeahKn+ed0JHvTtus9fRw0/wZMVqSmaCJ\nsUyKBjAyECAAAJBU13heG3eVa9OucjW1dtisjQrx1ZLMccpKi5X/aE87dQgAjoEAAQAYsaxWQwXf\n1uvzvDLtKaqRrTnRF3eKXpo5TskTQ+XGpGgAIxQBAgAw4pw936kt+Sf1Rd4JVZ8+Z7M20N9Lt82K\n120ZCQoN8rFThwDguAgQAIARwTAMHato0ud5ZfrmwCl1dltt1k+9IVhLM8dp1tRIjfJws1OXAOD4\nnCJAdHZ26rnnntPmzZvl7e2tn/3sZ3rwwQdtXlNZWak77rhD7777rtLS0uzUKQDA0Zxv79L2A6f0\n5c4T+u5Us81aHy8PLZgZqyWZCYqPGGOfBgHAyThFgHjllVdUVFSkdevWqbKyUk8++aSio6O1ePHi\nq17z3HPPqb3d9s6gAADX9d2pZn2584S27e9/CdaEyDFampmg+TNiNNp7lH0aBAAn5fABoq2tTR9/\n/LHef/99JSYmKjExUatXr9b69euvGiD++te/6vz583buFAAw3No7u7Wj4JS+3FmuoycbbdZ6uFs0\nZ1q0lmQmaPK4sbJYmBQNAGY4fIAoKSlRT0+PUlJSLh1LTU3VO++8c8X6xsZGvfHGG3r//ff1gx/8\nwF5tAgCGUXlNi77ceUJf7a3QufZum7VhQT66bVaCFmXEKcjf2z4NAoALcfgAUV9fr8DAQHl4fN9q\ncHCwOjo61NjYqKCg3hv3vPzyy7rrrrs0YcIEe7cKALCjzq4e5R2s0pe7yvvd8M3NIqVNjtDtsxM0\nfVKY3FmCFQCumcMHiLa2Nnl69t6k5+Lrzs7OXsfz8vJ04MABvfDCC3brDwBgX1X1rfpyV7m25J/U\n2fOdNmuDA7y1OCNei9LjWYIVAAaJwwcILy+vPkHh4msfn+9/GXR0dOi5557Ts88+2ydwAACcW3eP\nVbsP1+iLnWUqPNZgs9ZikaZPCtOS2QlKSwqXuztLsALAYHL4ABEeHq6mpiZZrVa5uV34JdDQ0CBv\nb2+NGfP9EnsHDx5URUWFHn74YRnG91uJ/sM//INWrFih5557zt6tAwCuU+2Z89q464Q2559U09kO\nm7WBfl5alBGnxRnxigj2tVOHADDyOHyASEpKkoeHhwoKCjRjxgxJ0t69ezV16tRedcnJydq0aVOv\nY4sWLdKLL76o2bNn261fAMD16emxam9xrb7YeUL7j9bpsu+ErmjahBAtyUxQxhQ2fAMAe3D4AOHt\n7a3ly5fr2Wef1UsvvaTa2lp98MEHevnllyVdGI3w9/eXl5eXYmNj+1wfFhamsWPH2rttAMAA1Z05\nr035F+Y2nG62vY+P/+hRykqL0+2zExQd6menDgEAkhMECEl6+umn9fzzz+uBBx6Qv7+/HnnkES1c\nuFCSNHfuXL388stasWJFn+tY0xsAHFt3j1W7j9Ro065yHfi2/9GGyePGasnsBGVOi5LnKHf7NAkA\n6MViGP39dT1yZGVlSZJycnKGuRMAcG1V9a3atLtcOXsq1NRqe26Dr7eHbp0Zq9tnJyg+YozNWgDA\n0HOKEQgAgPPr7OpR3qFqbdpVrkPHba+kJEk3xgVqyewEzU2Jlrcnv64AwFHwNzIAYEiV17Ro065y\nfbWvQmfPd9ms9fX20PwZMbptVoLGRwfYqUMAwEAQIAAAg669o1s7Ck9p465ylZQ39luflDBWt82K\n15zkKEYbAMDB8bc0AGDQlFY2adOucm0/UKnz7d02a/1He2rBzFgtzohTHHMbAMBpECAAANflfHuX\ntu+v1Mbd5Tpe2dxvffLEEN2WkaBZN0VolAcrKQGAsyFAAAAGzDAMFZWd0Zb8k/qm8JQ6Onts1gf5\ne2lhepwWpccrMoRdogHAmREgAACmnW5u09a9FdqSf1JVDeds1rpZpBmJ4VqcEa+0yeHycGeXaABw\nBQQIAIBNXd1W7S2u0abdJ7W/pFbWfnYPCgn00eL0OC1Mj1dokI99mgQA2A0BAgBwReXVLdqy56S+\n2leh5tZOm7XubhalT4nQ4ox4TZ8UJnc3i526BADYGwECAHDJubYufV1wSlvyy/XtyaZ+66ND/bQw\nPU5ZM2MVNMbbDh0CAIYbAQIARjir1dDh7xq0Of+k8gqr1NlttVnv4+WuucnRWpQer8SEIFksjDYA\nwEhCgACAEaqu8fylCdG1Z873Wz9lfLAWpsVpTnKUfLz49QEAIxW/AQBgBOns6tHuwzXanF+ugmP1\nMvqZED12jJey0uKUlRan6FA/+zQJAHBoBAgAcHGGYai0sklb91Ro2/5KtbZ12ay/OCF6UXqcZkwK\nkzvLrwIALkOAAAAXdbq5Tdv2VSpnb4Uqas/2Wx8X4a9F6fG6NTVGAX5edugQAOCMCBAA4ELaO7u1\n63CNtu45qcJj9f3u2TDa20Pzp8doYXqcJsYGMiEaANAvAgQAODnDMFRUdkY5e05qR2GV2jq6+71m\n2oQQLUqP06ybIuXtya8CAIB5/NYAACdVc/qctu6t0Na9FaZWUQobO1oLUmOVlRariGBfO3QIAHBF\nBAgAcCLn27u0o7BKW/dW6Mh3p/ut9/Hy0NzkKC2YGavJ44Llxg7RAIDrRIAAAAfXYzVU+G29tu6t\n0M7D1ers6rFZb7FIyRNDlTUzlkeUAACDjt8qAOCgyqqatW1fpbbtr9SZlvZ+62PD/bRgZpxuTY1R\ncICPHToEAIxEBAgADs0wDBXXl6qxvUlB3oFKCp3g0isF1Te2afuBSm3bV6Hymv6XXvUfPUrzp8fo\n1pmxrKIEALALAgQAh5VfWaB1hdmqba2/dCzcL1T3J69UekzKMHY2uFrbupRbWKVt+yt0+Hj/8xrc\n3SyamRSurLRYzUyK0CgPNnoDANgPAQKAQ8qvLNAbee/KMHpvZFDbWq838t7V2sw1Th0iurp7tKeo\nVtv2V2pPUa26e6z9XjMhJkC3zozV/Ols9AYAGD4ECAAOxzAMrSvM7hMeLj+/vjBbadHJTvXIjtVq\n6EjZaW3fX6kdhVU619bV7zUhgT6aPz1at6bGKj5yjB26BADANgIEAIdTXF/a67GlK6lprVdJQ6mS\nQifaqatrV17doq/2VWj7gVNqaGrrt97X20NzkqN1S2qMprD0KgDAwRAgADicxvYmc3VtzUPcybVr\naGrT1wcurKBUVtXSb72Hu5vSJofrlhkxmpkULs9R7nboEgCAgSNAAHA4Qd6B5up8Aoa4k4Fpbu1Q\n3sEqfV1wSke+O62rPIHVy9QbgnXLjBjNmRYlv9GeQ98kAADXiQABwOEkhU5QuF+ozceYIvxClRgy\nwY5dXdm5ti7tOlytrw+cUsGxelmt/aeGuAh/3TIjRvNnxCgsaLQdugQAYPAQIAA4HIvFovuTV15x\nFaaL51clrxy2CdTtHd3aU1Srrwsqtbe4ztQKSsEB3po/PUa3pMYoIXKMU03+BgDgcgQIAA4pPSZF\nazPXaH1htmouG4mI8AvVqmHYB6Kru0f7Sur0zYFT2l1Uo47Onn6vGe3toTnTojR/Roym3hAidyZD\nAwBcAAECgMNKj0lRWnSyiutL1dTerCCfACWG2G8n6p4eqwpLG/TNgVPaeahK59q7+73Gc5S70ieH\n6+bp0UpNZDI0AMD1ECAAODSLxaLJYfZbqtVqNVRUdlpfF5xSbmGVWs519nuNh7tFqYnhmpcSrfQp\nEfLx4q9WAIDr4rccgBHPajX07clG7Sis0o7CUzrd3N7vNW5uFiVPCNHN06M1a2okKygBAEYMAgSA\nEclqNVRSfka5hVXKO1ilBhOhQZKmjA/WvJRozZkWpUB/ryHuEgAAx0OAADBi9FgNlZw4ox2Fp5R3\nsFpnWsyFhomxgbp5erTmJkcrJNBniLsEAMCxESAAuLQeq6Gi704r9+CFkYbGsx2mrkuIHKN5KdGa\nlxKtyBDffusNw1Bxfaka25sU5B2opFD7TfYGAMCeCBAAXE5Pj1WH/y807DxUrSaToSEmzE9zkqM0\nLyVa8RFjTN8vv7JA6wqze218F+4XqvuHYblZAACGGgECgEvo6bHq0PEG5R6s1s5DVWpu7X/1JOnC\nrtBzpkVpTnLUgELDRfmVBVfc8K62tV5v5L2rtZlrCBEAAJdCgADgtDq6elRwtE55h6q1p6hGZ893\nmbouIXKM5iRHac60KMWG+1/z/Q3D0LrC7Cvuln3x/PrCbKVFJ/M4EwDAZRAgADiVc21d2lNcq52H\nqrS/pE7tJnaElqRxUd+Hhpiwaw8NlyuuL+312NKV1LTWq6ShVEmh9tvLAgCAoUSAAODwGs+2a/fh\nGu08VK2DpfXq7rnyN/5/74aYgAuPJ02LUlSo3+D31d5krq6tedDvDQDAcCFAAHBINafPaeehau08\nVK2S8jO6ylNCfUyICdCc5Av7NJhZPel6BHkHmqvzCRjSPgAAsCcCBACHYBiGymvO/l9oqFJZVYup\n69ws0uTxwZp9U6RmTY1UWNDoIe70e0mhExTuF2rzMaYIv1AlhkywW08AAAw1AgSAYdPdY1Vx2Rnt\nPlKj/CM1qj59ztR1Hu5uSrkxVJk3RSp9SoQC/IZnR2iLxaL7k1decRWmi+dXJa9kAjUAwKUQIADY\n1bm2Lu0vqdPuIzXaW1Krc23mVk7y8fLQzKRwzb4pUqmJYRrtPWqIOzUnPSZFazPXaH1htmouG4mI\n8AvVKvaBAAC4IAIEgCFXe+a88v9vlOHQ8Qb1WM1NaAjw81TGlEjNvilSyRNDNMrDfYg7vTbpMSlK\ni05WcX2pmtqbFeQToMQQdqIGALgmAgSAQWe1GiqtbFL+kRrtPlKjE9Xm5jNIUmiQj2ZPvRAaksYF\ny93NOT6EWywWTQ5jqVYAgOsjQAAYFB1dPSo8Vq/8IzXaU1SjMy0dpq+dEBOg9CmRypgSoXFRY/jm\nHgAAB0aAAHDN6hrPa29xrfYW16rwWIM6u8xt6ubh7qbkiSHKmBKhtMkRCgn0GeJOAQDAYCFAADCt\nu8eq4hNntK+4VnuKa3Wy5qzpa/1HeyptcrgypkRo+qQw+Xjx1w8AAM6I3+AAbGo82679JXXaU1yr\ngqN1OtfebframDC/S6MMiQljnWY+AwAAuDoCBIBeLk6Avvho0rGKJtPXulmkpHHBypgSofQpEYoO\n9RvCTgEAwHAgQABQ6/lOHfi2XnuLa7W/pE5NreYnQPuP9lRqYphmJoVr+qQwjfH1HMJOAQDAcHOK\nANHZ2annnntOmzdvlre3t372s5/pwQcfvGLttm3b9G//9m8qLy9XXFycHnnkES1YsMDOHQOOrcdq\n6FhFo/aX1Gn/0TodO9kok1szSJLGRwcoLSlcM5PCNTEuiEeTAAAYQZwiQLzyyisqKirSunXrVFlZ\nqSeffFLR0dFavHhxr7qjR4/q4Ycf1lNPPaWbb75ZX3/9tX71q1/pz3/+syZNmjRM3QOOoaGpTfuP\nXggMhd/Wq9XkDtCS5OPlrpQbL4wypCaGKTiAVZMAABipHD5AtLW16eOPP9b777+vxMREJSYmavXq\n1Vq/fn2fAPHZZ59p9uzZuu+++yRJ9913n7Zu3aovvviCAIERp6OrR0eOn74UGipqza+YJEnRoX5K\nm3xhlGHyuGCN8nAbok4BAIAzcfgAUVJSop6eHqWkpFw6lpqaqnfeeadP7V133aWurr7fqra2tg5p\nj4AjMAxDJ2vP6sDROu0vqdP/b+/uo6Ks8/+PvwaR2wG5RxDwBlQ0EpTMamvbxbusNq3NatvSY5mm\n2VE7Wqm7X63M2+qsJzOtU+6qW78t11U320rdXVPXSk3JNO9AQFRuBrm/F+b3BzI1gnahwAzyfJzD\nca7PfOaaN5wPl/Pic32u63Bqnqou1Bp+vZuri+Kig3TTxVOTwoK8W7BaAADQVjl9gMjNzZWfn59c\nXX8sNTAwUJWVlcrPz5e/v7+tvUePHnavPXHihL766is9+uijrVYv0JryiyuUfMKi707k6sCxHFkK\nK5r0+qjOPhrQO0T9e4fohh6Bcu/YoYUqBQAA1wunDxDl5eVyc7O/qkv9dlVV1WVfd/78eT377LNK\nTEzU4MGDW7RGoLWUVVTr+9Q8JZ/IVfLxXKU34UZukmT27KiEXsG20MAdoAEAQFM5fYBwd3dvEBTq\ntz09G//wY7FYNG7cOJlMJi1btqzFawRaSvWFWh3PyNfB47lKPpGr4xn5qmnC5ZJcTFLvrgHq3ztE\nA13X3f8AAB5FSURBVHoHKyaSKyYBAIBr4/QBIjQ0VAUFBaqtrZWLS90iTovFIg8PD/n6+jbon52d\nrTFjxqhDhw5au3at3SlOgLOrrbUqPavIFhgOp+apoqqmSfsI9ve0zTDE9wyW2bNjC1ULAADaI6cP\nEH369JGrq6sOHjyoAQMGSJL27dunuLi4Bn3Ly8s1fvx4dezYUWvWrFFAQEBrlws0idVq1bm8Uh06\naVHyCYuST+SqqPTyp+Y1xtPdVTdGBym+Z5D69w5RRIhZJhOzDAAAoGU4fYDw8PDQyJEjNXfuXC1Y\nsEDZ2dlavXq1Fi1aJKluNsLHx0fu7u5auXKlMjMztWbNGtXW1spisdj2YTabHfltAJLqAkNWXpkO\npVh06KRFh1IsymviwmfXDibFdgtQQs9gxfcMVs9IP3XowCVWAQBA6zBZrdYm3H/WMSoqKvTSSy/p\n888/l4+Pj8aPH6/HH39ckhQbG6tFixZp1KhRGjFihNLS0hq8ftSoUVq4cOHPvk/9Yuvt27c3a/1o\nv6xWq7LPl+m7i2Hh+5OWJl8pSaq783N8z2Al9AxW3+4B8nB3+uwPAACuU20iQLQWAgSuVX1gqJ9d\nOJSSJ0tBeZP30znQS/EXZxj6xQSpk9m9BaoFAABoOv6MCVyD+lOSDqfWhYVDKRbl5jc9MPj5uF9c\nxxCs+J5B6hzITdwAAIBzIkAATVBTa1X6uSIdTs3T4VN5+uFUns4XVTZ5P/WB4cboQMVFB7HwGQAA\ntBkECOAKqqprdOJ0gS0wHE07r7KKC03ej5/ZXXHRgboxJkg3EhgAAEAbRoAAfqK0vFo/pJ3XkVN5\nOpyap+MZBbpQU9vk/RAYAADA9YoAgXbLarUqt6Bcx9Ly6wLDqTylnSvS1VxWwM/srht61AeGQEWG\n+hAYAADAdYkAgXaj+kKNUs4U6mhavo6mndfR9PNNvgdDvbBAb/XtEaAbugfqhh6BCgvyJjAAAIB2\ngQCB69b5ogodTTuvH9LO61h6vk5mFqj6QtNPRzKZpO5hneoCQ49A9e0eqABfjxaoGAAAwPkRIHBd\nuFBTq7SzRTqaXhcYjqbnK+d82VXtq6Ori3pF+atv9wD17R6oPt0C5O3ZsZkrBgAAaJsIEGhzrFar\nLAUVOn46Xycy8nUsI18nTheosqrmqvZn9uyo2G4BtsDQM9JPbh07NHPVAAAA1wcCBJxeaXm1Tp4u\n0LGMfB3PyNeJ0/lXde+FelGdfRTbNUB9uvmrd9cAdQk2y8WF9QsAAABGECDgVKov1CrtXKGOZxTo\n+MXAkJlTctX783R3Ve8of8V2C1CfbgHq1dVfZk5HAgAAuGoECDhMba1VZy0lOplZaAsLqWcKr2qh\nc73wIG/Fdguo++rqr6jOvurA7AIAAECzIUCgVdTWWnUmt0QpmQU6mVmok5kFSj1TqPLKpt/VuZ63\nh6t6RvmrV5S/ekX6KbZbgDqZ3ZuxagAAAFyKAIFmV1Nr1ZmcYqWcqQsKKZmFSj1ToPLKq1vkLEmu\nHVzUo4uvekX6XwwNfgoPYu0CAABAayNA4JrU1FqVmVP848zC6QKdOluoiqu8IlK9LsFm9Yryq5td\niPJX93BfdXTlykgAAACORoCAYVXVNcrIKlbq2UKdOltYN7NwtvCqL59aL8DXXTER/urV1a9uhiHS\nT2Yvt2aqGgAAAM2JAIFGFRRXKvVsodLOFir1TJFOnStUZk6Jamut17TfAF8PxUT4KSaik6Ij/RQT\n4cddnQEAANoQAkQ7V1Nr1dncEp06W6jUM4U6da5Ip84UKr/46u+zUC+ok4eiI/wUE+mn6C6dFBPh\nJ3/CAgAAQJtGgGhHyiqqlXYxIJw6V6TUM4VKzypWVfW1nYIkSUF+noqJqAsJ0RF+io7oJH8fwgIA\nAMD1hgBxHaqqrlFmTonSs4qUfq5I6VnFSs8qUm5++TXv22SSwgK91T28k7p38b14OpIfl08FAABo\nJwgQbVhNTa3O5ZXWBYRzRRcDQ7HOWUp0jUsVJEnubh3UrbOvunfppO7hvuoR3kldw3zl6c6wAQAA\naK/4JNgGWK1W5RaUKyOrWGkXg0LGuWKdzim+prs2/1SAr4d6XAwK3cPr/g0LMnMXZwAAANghQDiR\nmlqrss+XKjO7RKez6wJCZnaJTucUq6zi6u/Y/FMuLiZFhfqo28UZhfrAwClIAAAAMIIA4QBV1TU6\nk1tiCwens4uVmVOiM7klzTajIEmhAV7q2tlXXcN8FNXZV107+ygixMwN2QAAAHDVCBAtqLS8+uIs\nQrFOXwwLmdklyj5f2ixrFOr5+7ira2dfRYX5qGtnX3UL81VkqA9rFQAAANDs+IR5jaov1Cr7fKnO\n5pbqTG7dLELd42KdL7r2eyn8lLeHa91MQpivunX2UVSYr6JCfTj9CAAAAK2GAGFAba1VeYUVOptb\nojOWn4aEEmWfL7vmuzNfysfLTZGhZkWG+igixEeRoWZ17eyrwE4eMplY1AwAAADHIUBcovpCrbbv\nzbALCWctpc1ys7VLBfl5KjLkYlAI9bE9ZkYBAAAAzooAcYnzRRX60/870Gz7c3ExKSzQu8GMQkQI\naxQAAADQ9vAJtpn4md0VHuytLsFmhQeb1SXYW+HBZoUHmdXR1cXR5QEAAADNggDRBJ7uHerCQdAl\nISHYLLNnR0eXBwAAALQ4AkQjIkLMDWYSugSb5e/jziJmAAAAtGsEiEuEBnjp7RcGO7oMAAAAwClx\ncj4AAAAAwwgQAAAAAAwjQAAAAAAwjAABAAAAwDACBAAAAADDCBAAAAAADCNAAAAAADCMAAEAAADA\nMAIEAAAAAMMIEAAAAAAMI0AAAAAAMIwAAQAAAMAwAgQAAAAAwwgQAAAAAAwjQAAAAAAwjAABAAAA\nwDACBAAAAADDCBAAAAAADCNAAAAAADCMAAEAAADAMAIEAAAAAMPaRICoqqrS7NmzNXDgQN1xxx1a\nvXr1ZfseOXJEDz30kBISEjR69GgdPny4FSsFAAAArm9tIkAsXrxYR44c0dq1azV37lwtX75cX3zx\nRYN+5eXlmjBhggYOHKgNGzYoISFBEydOVEVFhQOqBgAAAK4/Th8gysvLtX79ev3hD39QbGyshgwZ\novHjx2vdunUN+m7ZskWenp6aOXOmevTooTlz5sjb21ufffaZAyoHJKvVqiM5J7Q7Y6+O5JyQ1Wp1\ndEkAAADXxNXRBfyco0ePqqamRgkJCba2xMRErVq1qkHf7777TomJiXZtAwYM0IEDBzRq1KgWrxX4\nqW8yD2pt8gZll+Ta2kLNwXo8/gHdHJFwhVcCAAA4L6efgcjNzZWfn59cXX/MOoGBgaqsrFR+fr5d\n35ycHIWEhNi1BQYGKjs7u1VqBep9k3lQr//vHe09sl+7/rFLu/6yXbv+sUt7j+zX6/97R99kHnR0\niW3GsWPH9Nxzz2nMmDF67rnndOzYMUeX1KijR49q2rRp+v3vf69p06bp6NGjji4JcAomk6nBV72h\nQ4fatQ8dOrTJ+3AmI0eOtKtx5MiRji6p3UtOTtYTTzyhe+65R0888YSSk5MdXdJ1wekDRHl5udzc\n3Oza6rerqqrs2isqKhrte2k/oCVZrVatTd6g9PR0Hf/n9/K50EE+Xl7yudBBx//5vdLT07UueQOn\nMxlw7NgxvfjiiyovL5e3t7fKy8v14osvOl2IOHr0qJ5//nmVlZXJ29tbZWVlev755wkRaPcu90G/\nPixs27bNrn3btm0NQsSV9uFMRo4cqc2bN9u1bd68mRDhQMnJyZo8ebLKysrk4+OjsrIyTZ48mRDR\nDJw+QLi7uzcIAPXbnp6ehvp6eHi0bJHAT/yQe1LZJbnK2J+h4BB/u+eCQ/yVsT9DWSW5Omo56aAK\n245Vq1apc+fOdm2dO3du9BRGR1q5cmWjda5cudJBFQHO79Lw8HPtzu7S8PBz7Wh5y5YtU2RkpF1b\nZGSkli1b5qCKrh9OHyBCQ0NVUFCg2tpaW5vFYpGHh4d8fX0b9M3NzbVrs1gsCg4ObpVaAUnKryio\ne1BU2XiHi+355YWtVFHbZbFYmtTuKJced36uHQDQ8i53Cjuntl87pw8Qffr0kaurqw4e/PGc8X37\n9ikuLq5B3/j4eB04cMCu7cCBA3YLsIGW5u/hV/fA173xDhfb/T07tVJFbVdQUFCT2h3lcn+k4I8X\nAOA4oaGhTWqHcU4fIDw8PDRy5EjNnTtXhw4d0rZt27R69WqNHTtWUt1fIisr6/6iO3z4cBUXF2vB\nggVKSUnR/PnzVVZWphEjRjjyW0A70yc4RqHmYEUlRik3x36hf25OvqISo9TZHKzYoBgHVdh2TJw4\nUVlZWXZtWVlZmjhxooMqatzTTz/daJ1PP/20gyoCnN+QIUOa1O7s7rvvvia1o+VNnTpVp0+ftms7\nffq0pk6d6qCKrh9OHyAkadasWYqLi9PYsWP1yiuvaOrUqbYDzO23365//etfkiSz2ayVK1dq3759\n+u1vf6tDhw7p3XffZQ0EWpXJZNLj8Q+oa9eu6vWbOBW71qi4rEzFrjXq9Zs4de3aVY/FP+B0CwCd\nUe/evbVo0SJ5enqqtLRUnp6eWrRokXr37u3o0uzExsZqyZIl8vLyUmlpqby8vLRkyRLFxsY6ujTA\noS53sQir1aqtW7c2CAtDhgzR1q1bDe/DmWzatKlBWLjvvvu0adMmB1WE+Ph4rVixQl5eXiouLpaX\nl5dWrFih+Ph4R5fW5pmszvYb6ECDBw+WJG3fvt3BleB68E3mQa1L3qCsn9wHorM5WI9xHwgAANCG\nOf2N5IC26uaIBA3sEq8fck+qoKJQ/p6dFBsUw8wDAABo0wgQQAsymUzqG9LT0WUAAAA0mzaxBgIA\nAACAcyBAAAAAADCMAAEAAADAMAIEAAAAAMMIEAAAAAAMI0AAAAAAMIwAAQAAAMAwAgQAAAAAwwgQ\nAAAAAAwjQAAAAAAwjAABAAAAwDACBAAAAADDCBAAAAAADCNAAAAAADCMAAEAAADAMAIEAAAAAMMI\nEAAAAAAMI0AAAAAAMIwAAQAAAMAwAgQAAAAAwwgQAAAAAAwjQAAAAAAwjAABAAAAwDACBAAAAADD\nCBAAAAAADCNAAAAAADCMAAEAAADAMAIEAAAAAMMIEAAAAAAMI0AAAAAAMIwAAQAAAMAwAgQAAAAA\nwwgQAAAAAAwjQAAAAAAwjAABAAAAwDACBAAAAADDCBAAAAAADCNAAAAAADCMAAEAAADAMAIEAAAA\nAMMIEAAAAAAMI0AAAAAAMIwAAQAAAMAwAgQAAAAAwwgQAAAAAAwjQAAAAAAwjAABAAAAwDACBAAA\nAADDCBAAAAAADCNAAAAAADCMAAEAAADAMAIEAAAAAMMIEAAAAAAMaxMB4rXXXtOtt96qQYMGaenS\npVfse/DgQT3yyCPq37+/RowYoY8//riVqgQAAACuf66OLuDnvP/++9qyZYtWrFih6upqzZgxQ0FB\nQRo3blyDvhaLRRMmTNCjjz6qJUuW6Pvvv9esWbMUEhKiO++80wHVAwAAANcXp5+BWLt2raZOnar+\n/fvr5ptv1owZM7Ru3bpG+27btk3BwcGaNm2aoqKidPfdd2vkyJH65JNPWrlqAAAA4Prk1DMQOTk5\nOnfunG666SZbW2Jios6ePSuLxaKgoCC7/r/85S/Vt2/fBvspLi5u8VoBAACA9sCpZyByc3NlMpkU\nEhJiawsKCpLValVWVlaD/uHh4erXr59tOy8vT59++qluu+22VqkXAAAAuN45fAaisrJS2dnZjT5X\nVlYmSXJzc7O11T+uqqr62f0+++yzCgkJ0cMPP2yoltzcXF24cEGDBw821B8AAABoS8LCwi67HMAo\nhweI5ORkjRkzRiaTqcFzM2bMkFQXFi4NDp6enpfdZ1lZmSZNmqSMjAx9+OGHcnd3N1SLm5ubrFZr\nU78FAAAAoN0wWZ34E3NOTo7uvPNObd++XeHh4ZKkzMxMDR06VDt37mywBkKSSkpKNH78eGVmZuov\nf/mLoqOjW7tsAAAA4Lrl1GsgQkJCFBYWpv3799va9u3bp7CwsEbDg9Vq1ZQpU3TmzBmtW7eO8AAA\nAAA0M4efwvRzHnnkEb322msKDQ2V1WrVG2+8oSeffNL2/Pnz5+Xh4SEvLy99/PHH+uabb/T222/L\nbDbLYrFIkjp27KhOnTo56lsAAAAArhtOfQqTJNXW1mrp0qXasGGDXFxc9NBDD2n69Om255OSkvTA\nAw9oypQpGj9+vHbv3t1gHwMHDtSaNWtas2wAAADguuT0AQIAAACA83DqNRAAAAAAnAsBAgAAAIBh\nBAgAAAAAhhEgAAAAABjWbgNERkaGnnzySfXv319JSUl67733bM9lZmZq3Lhx6t+/v+69995Gr+wE\nOMqVxu78+fMVGxurPn362P7961//6sBqgYYmTJigWbNm2baPHDmihx56SAkJCRo9erQOHz7swOqA\nxl06bidNmtTgeLtjxw4HVgj8aNu2bQ3G59SpUyU1zzG3XQYIq9WqCRMmKCgoSJs2bdK8efP09ttv\na8uWLZKkyZMnKyQkRH//+9913333acqUKcrKynJw1cDPj93U1FTNmDFDu3bt0u7du7Vr1y49+OCD\nDq4a+NGWLVv05Zdf2rbLy8s1YcIEDRw4UBs2bFBCQoImTpyoiooKB1YJ2Lt03Ep1x9vXX3/d7nh7\n2223OahCwN7JkyeVlJSk3bt328bnq6++2mzH3HYZICwWi/r27au5c+cqKipKv/zlL3Xrrbdq//79\n+uqrr5SZmamXX35ZPXr00IQJE5SQkKD169c7umzgimNXklJSUtS3b18FBgbavtzd3R1cNVCnsLBQ\nS5cuVb9+/WxtW7Zskaenp2bOnKkePXpozpw58vb21meffebASoEfNTZuq6qqlJmZqbi4OLvjbceO\nHR1YKfCjlJQU9ezZUwEBAbbxaTabm+2Y2y4DRHBwsN544w15eXlJkvbv3699+/bp5ptvVnJysm64\n4Qa7D12JiYk6ePCgo8oFbBobu3v37tWgQYNUUlKi7OxsdevWzbFFApexePFijRw5UtHR0ba27777\nTomJiXb9BgwYoAMHDrR2eUCjGhu3qampMplMioiIcGBlwOWlpKSoe/fuDdqb65jbLgPETyUlJemx\nxx5TQkKChg0bptzcXIWEhNj1CQwMVHZ2toMqBBpXP3b79++vYcOGKSUlRSaTSW+//bbuvPNOjRw5\nUhs3bnR0mYAkac+ePdq/f7+eeeYZu/acnByOuXBalxu3qampMpvNev7553X77bdr9OjRDU5xAhzp\n1KlT2rlzp4YPH66hQ4fqjTfeUHV1dbMdc12bs9i26M0335TFYtG8efO0YMEClZeXy83Nza6Pm5ub\nqqqqHFQh0Lj6sTt37ly9+uqriouLk4uLi6Kjo/X444/rm2++0R//+EeZzWYNGTLE0eWiHauqqtK8\nefM0d+7cBsfXiooKjrlwSlcat6mpqaqsrNQdd9yhCRMmaOvWrZo0aZI++ugj3XDDDQ6qGKhz9uxZ\nVVRUyN3dXcuWLVNmZqZt/UNzHXPbfYCo/0V/8cUXNWPGDD344IMqKiqy61NVVSUPDw9HlAdcVv3Y\nnTVrlmbOnKkXXnhBSUlJ8vX1lST16tVLaWlp+vDDDwkQcKg333xTcXFxjS4wdXd3b/AfF8dcOIMr\njdspU6Zo7Nix8vHxkST17t1b33//vf72t7/p5Zdfbu1SATvh4eH6+uuvbZ8HYmNjVVtbq5kzZ2rQ\noEHNcsxtlwEiLy9PBw4csPtQFRMTo+rqagUHByslJcWuv8ViUXBwcGuXCTRwpbFbWloqPz8/u/49\nevTQ119/3dplAnY+/fRT5eXlqX///pKk6upqSdLnn3+ue++9V7m5uXb9OebCGVxp3H777be28FAv\nOjq6wecHwFHqw0O96OhoVVZWKigoqFmOue1yDURmZqaeffZZux/goUOHFBgYqMTERB0+fNgune3f\nv18JCQmOKBWwc7mxGxAQoDVr1mjcuHF2/X/44YdGF1EBrWndunX65z//qc2bN2vz5s1KSkpSUlKS\nNm3apPj4+AaL9w4cOMAxFw53pXE7a9YszZkzx67/0aNHOd7CKezatUuDBg1SZWWlre3IkSPy9/fX\nTTfdpG+//dau/9Ucc9tlgLjxxhsVFxenWbNmKSUlRTt27NBrr72mSZMmaeDAgQoLC9OLL76okydP\n6p133tGhQ4e4lj6cwpXG7q9//Wvt3btXq1ev1unTp/XBBx9o8+bNGj9+vKPLRjsXFhamyMhI25e3\nt7e8vb0VGRmp4cOHq7i4WAsWLFBKSormz5+vsrIyjRgxwtFlo5270rgdPHiwNm/erI0bNyojI0PL\nly/Xt99+q8cff9zRZQPq37+/PD09NWfOHJ06dUo7duzQ0qVL9dRTT2nYsGHNcsw1Wa1WawvV79Ry\nc3P1yiuvaM+ePfL09NRjjz2mCRMmSJJOnz6t2bNn67vvvlNUVJTmzJmjW265xcEVA3WuNHb//e9/\na9myZUpPT1eXLl00ffp01j/A6dTfzXfhwoWS6mbR5s6dq9TUVPXu3VsvvfSSYmNjHVki0MCl43b9\n+vV69913lZWVpZiYGM2ePbvB5TEBR0lJSdGCBQt08OBBeXt765FHHtHkyZMlNc8xt90GCAAAAABN\n1y5PYQIAAABwdQgQAAAAAAwjQAAAAAAwjAABAAAAwDACBAAAAADDCBAAAAAADCNAAAAAADCMAAEA\nAADAMAIEAAAAAMMIEACABpKSkrR8+fKrfv3Jkye1Y8cO23ZsbKw2btwoSbpw4YL+/Oc/X2uJdvsE\nALQeAgQAoNlNnDhRhw4dsm3v3r1bd999tyTpk08+0eLFix1VGgDgGrk6ugAAwPXHarXabQcGBtoe\n19bWtnY5AIBmxAwEAKBJqqqqtHjxYg0ePFhxcXEaNGiQpk2bpvz8fEl1pz+dO3dOy5cv15gxYyT9\neLrRP/7xD82ePVtWq1V9+vTR3r17tXz5ciUlJdm9x5tvvmnXlp2drUmTJmnAgAH61a9+pU8++aRB\nXf/5z3/0wAMPKD4+XsOGDdOyZctUVVXVgj8JAGifCBAAgCZZunSptm3bpsWLF2vr1q1avHix9uzZ\no5UrV0qS1q9fr9DQUD3xxBN666237F57zz33aPbs2TKZTNq9e7cSEhIkSSaTya6fyWSytdXU1OjJ\nJ59UYWGhPvjgAy1btkzvvfee3Wu+/PJLTZ8+Xb/73e+0ZcsWzZs3T5999pleeOGFlvxRAEC7xClM\nAIAm6devn+666y4lJiZKksLCwvSLX/xCx48flyQFBATIxcVFXl5e8vHxsXutm5ubrS0gIMDQ+/3v\nf/9TSkqKtm7dqoiICEnSwoULNWrUKFufVatW6eGHH9bo0aMlSREREZo3b57Gjh2rmTNnKjw8/Nq+\naQCADQECANAkv/nNb7Rnzx69/vrrSktLU2pqqk6dOqWbbrqpRd7vxIkT8vX1tYUHqe6UKE9PT9v2\nkSNHdOjQIX300Ud2r3VxcVFKSgoBAgCaEQECANAk//d//6cvvvhC999/vwYPHqxnnnlG7733nrKz\ns5vtPS5cuPCzfVxdf/wvrLa2VuPHj9f999/foF9wcHCz1QUAIEAAAJqgoKBAH330kf70pz/prrvu\nsrWnpKTI29vbtn3pmoafuvS5jh07qrS01K4tLS3N9rhPnz4qKipSSkqKoqOjbc8XFxfb+vTs2VOn\nTp1SZGSkre3rr7/W2rVr9dJLL8nDw6Np3ygA4LIIEACARqWnp2vnzp12be7u7vL19dW2bdvUt29f\nlZeXa926dTpy5IhtQbQkeXl5KT09XXl5eXaXcK1/TpIOHz6smJgYJSQkqLCwUO+//76GDx+unTt3\naufOnfLz85Mk3XLLLerXr59mzpypuXPnysXFRa+++qo6dOhg2+dTTz2l6dOn66233tI999yjc+fO\nac6cOYqKimrw/gCAa2OyXnqxbgBAu1d/KdZLhYeHa/78+Vq4cKEyMjLUqVMnDRo0SDExMXrnnXe0\ne/duubu76+OPP9aSJUvUpUsXbdy4UX369LEtfC4qKtJTTz2lH374QUuXLtXw4cO1YsUKffDBByot\nLdUdd9yhxMRErVmzRtu3b5ckFRYW6pVXXtF///tfeXh4aOLEiVq1apVmzJhhW0z9+eefa9WqVTp5\n8qQ6deqkwYMHa8aMGTKbza36swOA6x0BAgAAAIBh3AcCAAAAgGEECAAAAACGESAAAAAAGEaAAAAA\nAGAYAQIAAACAYQQIAAAAAIYRIAAAAAAYRoAAAAAAYBgBAgAAAIBhBAgAAAAAhhEgAAAAABhGgAAA\nAABg2P8HGXELsQbiTR4AAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plot_log_reg(x='Latitude', y='Allele', data=raw_df, clf=clf, xmin=30, xmax=50, alpha=0.02)\n", "df.sort_values('Latitude').plot('Latitude', 'p, Mpi100', ls='', marker='o', ax=ax)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0.19685284] [[ 0.]]\n" ] } ], "source": [ "clf0 = log_reg_null_model(raw_df['Latitude'], raw_df['Allele'])\n", "print(clf0.intercept_, clf0.coef_)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "7.0572187317602719e-20" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "log_reg_lik_ratio_test(raw_df['Latitude'], raw_df['Allele'], clf0, clf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Multiple logistic regression\n", "\n", "Now following the [chapter of the same name](http://www.biostathandbook.com/multiplelogistic.html) in the same book.\n", "\n", "Here is an example using the data on bird introductions to New Zealand." ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
speciesstatuslengthmassrangemigrinsectdietclutchbroodswooduplandwaterreleaseindiv
0Cyg_olor115209600.01.21112.026.01.0001629.0
1Cyg_atra112505000.00.5610.016.01.00011085.0
2Cer_nova18703360.00.0710.014.01.000138.0
3Ans_caer07202517.01.10312.023.81.0001110.0
4Ans_anse08203170.03.4530.015.91.000127.0
\n", "
" ], "text/plain": [ " species status length mass range migr insect diet clutch \\\n", "0 Cyg_olor 1 1520 9600.0 1.21 1 12.0 2 6.0 \n", "1 Cyg_atra 1 1250 5000.0 0.56 1 0.0 1 6.0 \n", "2 Cer_nova 1 870 3360.0 0.07 1 0.0 1 4.0 \n", "3 Ans_caer 0 720 2517.0 1.10 3 12.0 2 3.8 \n", "4 Ans_anse 0 820 3170.0 3.45 3 0.0 1 5.9 \n", "\n", " broods wood upland water release indiv \n", "0 1.0 0 0 1 6 29.0 \n", "1 1.0 0 0 1 10 85.0 \n", "2 1.0 0 0 1 3 8.0 \n", "3 1.0 0 0 1 1 10.0 \n", "4 1.0 0 0 1 2 7.0 " ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data=\"\"\"\n", "Cyg_olor 1 1520 9600 1.21 1 12 2 6 1 0 0 1 6 29\n", "Cyg_atra 1 1250 5000 0.56 1 0 1 6 1 0 0 1 10 85\n", "Cer_nova 1 870 3360 0.07 1 0 1 4 1 0 0 1 3 8\n", "Ans_caer 0 720 2517 1.1 3 12 2 3.8 1 0 0 1 1 10\n", "Ans_anse 0 820 3170 3.45 3 0 1 5.9 1 0 0 1 2 7\n", "Bra_cana 1 770 4390 2.96 2 0 1 5.9 1 0 0 1 10 60\n", "Bra_sand 0 50 1930 0.01 1 0 1 4 2 0 0 0 1 2\n", "Alo_aegy 0 680 2040 2.71 1 . 2 8.5 1 0 0 1 1 8\n", "Ana_plat 1 570 1020 9.01 2 6 2 12.6 1 0 0 1 17 1539\n", "Ana_acut 0 580 910 7.9 3 6 2 8.3 1 0 0 1 3 102\n", "Ana_pene 0 480 590 4.33 3 0 1 8.7 1 0 0 1 5 32\n", "Aix_spon 0 470 539 1.04 3 12 2 13.5 2 1 0 1 5 10\n", "Ayt_feri 0 450 940 2.17 3 12 2 9.5 1 0 0 1 3 9\n", "Ayt_fuli 0 435 684 4.81 3 12 2 10.1 1 0 0 1 2 5\n", "Ore_pict 0 275 230 0.31 1 3 1 9.5 1 1 1 0 9 398\n", "Lop_cali 1 256 162 0.24 1 3 1 14.2 2 0 0 0 15 1420\n", "Col_virg 1 230 170 0.77 1 3 1 13.7 1 0 0 0 17 1156\n", "Ale_grae 1 330 501 2.23 1 3 1 15.5 1 0 1 0 15 362\n", "Ale_rufa 0 330 439 0.22 1 3 2 11.2 2 0 0 0 2 20\n", "Per_perd 0 300 386 2.4 1 3 1 14.6 1 0 1 0 24 676\n", "Cot_pect 0 182 95 0.33 3 . 2 7.5 1 0 0 0 3 .\n", "Cot_aust 1 180 95 0.69 2 12 2 11 1 0 0 1 11 601\n", "Lop_nyct 0 800 1150 0.28 1 12 2 5 1 1 1 0 4 6\n", "Pha_colc 1 710 850 1.25 1 12 2 11.8 1 1 0 0 27 244\n", "Syr_reev 0 750 949 0.2 1 12 2 9.5 1 1 1 0 2 9\n", "Tet_tetr 0 470 900 4.17 1 3 1 7.9 1 1 1 0 2 13\n", "Lag_lago 0 390 517 7.29 1 0 1 7.5 1 1 1 0 2 4\n", "Ped_phas 0 440 815 1.83 1 3 1 12.3 1 1 0 0 1 22\n", "Tym_cupi 0 435 770 0.26 1 4 1 12 1 0 0 0 3 57\n", "Van_vane 0 300 226 3.93 2 12 3 3.8 1 0 0 0 8 124\n", "Plu_squa 0 285 318 1.67 3 12 3 4 1 0 0 1 2 3\n", "Pte_alch 0 350 225 1.21 2 0 1 2.5 2 0 0 0 1 8\n", "Pha_chal 0 320 350 0.6 1 12 2 2 2 1 0 0 8 42\n", "Ocy_loph 0 330 205 0.76 1 0 1 2 7 1 0 1 4 23\n", "Leu_mela 0 372 . 0.07 1 12 2 2 1 1 0 0 6 34\n", "Ath_noct 1 220 176 4.84 1 12 3 3.6 1 1 0 0 7 221\n", "Tyt_alba 0 340 298 8.9 2 0 3 5.7 2 1 0 0 1 7\n", "Dac_nova 1 460 382 0.34 1 12 3 2 1 1 0 0 7 21\n", "Lul_arbo 0 150 32.1 1.78 2 4 2 3.9 2 1 0 0 1 5\n", "Ala_arve 1 185 38.9 5.19 2 12 2 3.7 3 0 0 0 11 391\n", "Pru_modu 1 145 20.5 1.95 2 12 2 3.4 2 1 0 0 14 245\n", "Eri_rebe 0 140 15.8 2.31 2 12 2 5 2 1 0 0 11 123\n", "Lus_mega 0 161 19.4 1.88 3 12 2 4.7 2 1 0 0 4 7\n", "Tur_meru 1 255 82.6 3.3 2 12 2 3.8 3 1 0 0 16 596\n", "Tur_phil 1 230 67.3 4.84 2 12 2 4.7 2 1 0 0 12 343\n", "Syl_comm 0 140 12.8 3.39 3 12 2 4.6 2 1 0 0 1 2\n", "Syl_atri 0 142 17.5 2.43 2 5 2 4.6 1 1 0 0 1 5\n", "Man_mela 0 180 . 0.04 1 12 3 1.9 5 1 0 0 1 2\n", "Man_mela 0 265 59 0.25 1 12 2 2.6 . 1 0 0 1 80\n", "Gra_cyan 0 275 128 0.83 1 12 3 3 2 1 0 1 1 .\n", "Gym_tibi 1 400 380 0.82 1 12 3 4 1 1 0 0 15 448\n", "Cor_mone 0 335 203 3.4 2 12 2 4.5 1 1 0 0 2 3\n", "Cor_frug 1 400 425 3.73 1 12 2 3.6 1 1 0 0 10 182\n", "Stu_vulg 1 222 79.8 3.33 2 6 2 4.8 2 1 0 0 14 653\n", "Acr_tris 1 230 111.3 0.56 1 12 2 3.7 1 1 0 0 5 88\n", "Pas_dome 1 149 28.8 6.5 1 6 2 3.9 3 1 0 0 12 416\n", "Pas_mont 0 133 22 6.8 1 6 2 4.7 3 1 0 0 3 14\n", "Aeg_temp 0 120 . 0.17 1 6 2 4.7 3 1 0 0 3 14\n", "Emb_gutt 0 120 19 0.15 1 4 1 5 3 0 0 0 4 112\n", "Poe_gutt 0 100 12.4 0.75 1 4 1 4.7 3 0 0 0 1 12\n", "Lon_punc 0 110 13.5 1.06 1 0 1 5 3 0 0 0 1 8\n", "Lon_cast 0 100 . 0.13 1 4 1 5 . 0 0 1 4 45\n", "Pad_oryz 0 160 . 0.09 1 0 1 5 . 0 0 0 2 6\n", "Fri_coel 1 160 23.5 2.61 2 12 2 4.9 2 1 0 0 17 449\n", "Fri_mont 0 146 21.4 3.09 3 10 2 6 . 1 0 0 7 121\n", "Car_chlo 1 147 29 2.09 2 7 2 4.8 2 1 0 0 6 65\n", "Car_spin 0 117 12 2.09 3 3 1 4 2 1 0 0 3 54\n", "Car_card 1 120 15.5 2.85 2 4 1 4.4 3 1 0 0 14 626\n", "Aca_flam 1 115 11.5 5.54 2 6 1 5 2 1 0 0 10 607\n", "Aca_flavi 0 133 17 1.67 2 0 1 5 3 0 1 0 3 61\n", "Aca_cann 0 136 18.5 2.52 2 6 1 4.7 2 1 0 0 12 209\n", "Pyr_pyrr 0 142 23.5 3.57 1 4 1 4 3 1 0 0 2 .\n", "Emb_citr 1 160 28.2 4.11 2 8 2 3.3 3 1 0 0 14 656\n", "Emb_hort 0 163 21.6 2.75 3 12 2 5 1 0 0 0 1 6\n", "Emb_cirl 1 160 23.6 0.62 1 12 2 3.5 2 1 0 0 3 29\n", "Emb_scho 0 150 20.7 5.42 1 12 2 5.1 2 0 0 1 2 9\n", "Pir_rubr 0 170 31 0.55 3 12 2 4 . 1 0 0 1 2\n", "Age_phoe 0 210 36.9 2 2 8 2 3.7 1 0 0 1 1 2\n", "Stu_negl 0 225 106.5 1.2 2 12 2 4.8 2 0 0 0 1 2\n", "\"\"\"\n", "cols = 'species status length mass range migr insect diet clutch broods wood upland water release indiv'.split()\n", "ind_vars = cols[2:]\n", "df = pd.read_table(io.StringIO(data), names=cols, header=None, delim_whitespace=True, na_values=['.'])\n", "df.dropna(inplace=True)\n", "y_true = df['status']\n", "X = df[ind_vars]\n", "df.head()" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def log_reg_fit(X, y):\n", " clf = sklearn.linear_model.LogisticRegression(C=1e12)\n", " if X.ndim == 1:\n", " X = X.reshape(-1, 1)\n", " clf.fit(X, y)\n", " return clf" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Adding release with p-value 4.9e-09\n", "Adding upland with p-value 0.0044\n", "Adding migr with p-value 0.011\n" ] }, { "data": { "text/plain": [ "['release', 'upland', 'migr']" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def forward_selection(X, y_true, pvals_threshold = 0.05):\n", " ind_vars = X.columns\n", " clf0 = log_reg_null_model(y_true, y_true)\n", " included_vars = []\n", " excluded_vars = list(ind_vars)\n", "\n", " while excluded_vars:\n", " selected_var, pval = None, 1\n", " for var in excluded_vars:\n", " X_ = X[included_vars + [var]]\n", " clf = log_reg_fit(X_, y_true)\n", " pval_ = log_reg_lik_ratio_test(X_, y_true, clf0, clf)\n", " if pval_ < pval:\n", " selected_var, pval = var, pval_\n", " if pval < pvals_threshold:\n", " print(\"Adding {} with p-value {:.2g}\".format(selected_var, pval))\n", " included_vars.append(selected_var)\n", " excluded_vars.remove(selected_var)\n", " X_ = X[included_vars].copy()\n", " X_['placeholder'] = 0\n", " clf0 = log_reg_fit(X_, y_true)\n", " else:\n", " break\n", " return included_vars\n", "\n", "forward_selection(X, y_true)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I get a similar result to the one presented by John, but I need a smaller P value threshold (I used 0.05 instead of 0.15). With 0.15 I get several more variables. Not sure why... Maybe the P value calculation in the book is using Wald test (used by SAS in the example given in the book), whereas I'm using the likelihood ratio test (which I believe is a better choice here)." ] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:Py4Eng]", "language": "python", "name": "conda-env-Py4Eng-py" }, "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.5.2" } }, "nbformat": 4, "nbformat_minor": 1 }