week05b-wcoa-cruise-regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Implementing linear regression in Python\n",
    "\n",
    "## Example: 2007 West Coast Ocean Acidification Cruise "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy import stats\n",
    "import statsmodels.formula.api as smf\n",
    "import pingouin as pg\n",
    "\n",
    "import PyCO2SYS as pyco2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data\n",
    "\n",
    "##### Load 2007 data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename07 = 'wcoa_cruise_2007/32WC20070511.exc.csv'\n",
    "df07 = pd.read_csv(filename07,header=29,na_values=-999,parse_dates=[[6,7]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['DATE_TIME', 'EXPOCODE', 'SECT_ID', 'SAMPNO', 'LINE', 'STNNBR',\n",
       "       'CASTNO', 'LATITUDE', 'LONGITUDE', 'BOT_DEPTH', 'BTLNBR',\n",
       "       'BTLNBR_FLAG_W', 'CTDPRS', 'CTDTMP', 'CTDSAL', 'CTDSAL_FLAG_W',\n",
       "       'SALNTY', 'SALNTY_FLAG_W', 'CTDOXY', 'CTDOXY_FLAG_W', 'OXYGEN',\n",
       "       'OXYGEN_FLAG_W', 'SILCAT', 'SILCAT_FLAG_W', 'NITRAT', 'NITRAT_FLAG_W',\n",
       "       'NITRIT', 'NITRIT_FLAG_W', 'PHSPHT', 'PHSPHT_FLAG_W', 'TCARBN',\n",
       "       'TCARBN_FLAG_W', 'ALKALI', 'ALKALI_FLAG_W'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df07.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear regression: four methods in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create $x$ and $y$ variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = df07['PHSPHT']\n",
    "y = df07['NITRAT']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'nitrate')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEJCAYAAACdePCvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABGzklEQVR4nO29e3wU9b3//5rZTUgCCRBzBRIWAhiFANWAXNQIKheBIqDn+Ci1aH+1oBYf5fR4qXpsq7VwkB4rpYrl0WOtB885IqhIjwEkyreahBAUEpCYmLDJEnLZkJDb5rK78/n9MTuTmb3OJrvJJvt+Ph627szuzCfz8PF+f+Z9eb05xhgDQRAEEXbwQ70AgiAIYmggB0AQBBGmkAMgCIIIU8gBEARBhCnkAAiCIMIU/VAvQCvd3d04f/48EhMTodPphno5BEEQwwK73Q6z2YxZs2YhKipKdW7YOIDz589j48aNQ70MgiCIYcn+/fuRnZ2tOjZsHEBiYiIA8Y9ISUkZ4tUQBEEMD+rr67Fx40bZhioZNg5ACvukpKRg0qRJQ7wagiCI4YW70DklgQmCIMIUcgAEQRBhCjkAgiCIMIUcAEEQRJhCDoAgCCJMIQdAEAShgTJjMw6cKEeZsXmolxIwhk0ZKEEQxFBQZmxGXrEJx4uqIQgMej2Pl7csRqYhXtNvSyubkJWRoOn7gw05AIIgCA+UGZvx3N4v0WsV5GM2m4DSyiafBl36rc0m+OU0BhMKAREEQXigtLIJNluf8ecA6PU8sjISfP42r9iEXqsAgfU5jVCD3gAIgiA8kJWRAL2eh80mQKfjcde8NCzNTte0+z9eVC1/1ul4xMZE4sCJ8pAKB5EDIAiC8ECmIR4vb1msOY4vxfzLa1pgt/eNWx8Trcfr758DAxAZweORtVlot/QiNiYS7ZbeIXMK5AAIgiC8kGmI15zwfW7vl7BaBTCncy3tvfK/91oF/On9c/JnjgMi9H1OYTCdATkAgiAIDfiq6Dn4WYUqWawVxgCrTcDeQyVgzL8qo4FCDoAgiLDHl3FX7u45Dpg/MwU3ZybLO/bC83UoPF/f/wUwwM7E9watVUaBgBwAQRBhjZZyzdLKJjm0wxhQeL4ehefr5aogZaVQf1CGjDgOmqqMAgGVgRIEEdZIpZ7eyjWzMhLA85zLcQbxN84xf0+4XsEVQevFAgA5AIIgwhqp1JPnPNf4ZxrisWX9bLjxAW4dgye02HbGxB6CwYBCQARBhDXuSj2dcwJlxma0W3rx6IY5OH6qGuWma/Lvx0RHoLWz1/MN+sXgvAaQAyAIIqxwl/BVlno65wQeWZuFNz8ogd3OoNNxWHt7hsoBBNr463UclmanB/SaHu81KHchCIIIAbQmfKWcgNUq4MCJctgcTV02O8Pfv7wU8HWtXGjA0uy0QReOIwdAEETY4Gzc84pNLsY2KyMBOh0PwZHcbWzpUp3v7rUHdE0Reh5Ls9M0N5wFEkoCEwQxYvCm2V9mbIa5xSInbRmAY6eMLt/NNMTj5sykwVguoiJ1+N2jQ6cSSm8ABEGMCMqMzfjl61/Isfrtj93qNq7PFPlVuwDs+7AUj9ybJX83t8A4sKYuAKkJMahrsvj8XvQonVvjn1tgRH7JFSyaPQErFhoGtBZvkAMgCGJEkFdco4rV5xXXABDDPuYWixz6cabcdA3P7f0SL29ZDAB4XaHT0x9mTo3Ht9Utmr57rb0XuQVGlQZQboFR1gr6utwMAEFzAuQACIIYIajr8Vvae+Rdv07Hi6EfgckhIJtCrVPMB9QgcXyMzwLMSD2P6enj0NbZi8sNHS7fb+vsVV3bGwxw0QDKL7mi+k5+yRVyAARBEN5Ymp2GT0/XyBU+42Oj+nb9dgHLFxiQOD5abvTKKzbhaKERAhMN8aenTfjpvVnQ8WJoyBO9NgEXqtzPBdbxwMTEMTA1dGhaM88BgsDkjuLSyiYsmj1B3vkDwKLZEzQ+Af8JugP493//d7S0tGDHjh3Iz8/H9u3b0dPTg5UrV2Lbtm3Bvj1BEGFCpiEeP703C/klVzB14lhYuq3QOay5XlFpI2Gsa1PlA2w2AWfKGpAxcRzqrnag3WLzew3LbhHLOU+dr/f6JnHd2ChMTxuHmzOTse+jUtlpKUtAh30OoKCgAB988AHuuOMOdHd349lnn8U777yD1NRUbN68GSdPnkROTk4wl0AQRJhQZmzGvo9K0WsV8HW5GRwAnY7D8gUGF+NfZmzGG4dKVEaaAQNK/kZGiE4GjvtKYSAO6r7eCD2PZ340T16PITXOpf5/xUJDUA2/RNAcwLVr1/Dqq69iy5YtKCsrQ0lJCSZPnoy0NPEBrVmzBrm5uW4dQFtbG9ra2lTH6usHlpUnCGJkInX2mlsssCr0+BnEOH/91U4Y69qQV1yDlvYeAEBtYwcEH6prY8dEorVDW5cvB+CRtWIl0YET5apr3zg1HunJsWhp78H42CjZGSk7ku+/c4bff3cgCJoDeOGFF7Bt2zbU1dUBABobG5GYmCifT0pKQkNDg9vfvv3229izZ0+wlkYQxDDCm1a/srxTp+PBcVCFdQCxkkYZU9dKu0W7xAMDUFXbigMnyhEbEwme5yA43gAuVDXjm6pmRET0dR5r6UgeDILiAA4cOIDU1FQsXLgQhw4dAgAIggCO68vSM8ZUn5Vs2rQJ69atUx2rr6/Hxo0bg7FcgiBCFF+GUtnZC7uAG6bEe0zQ+sOUCbGormvX/H29jsPxomoIgljNc/3k8ap1KJO8mYZ41XwBq3XwBsC4rDsYF/2///s/mM1mrF27Fq2trbBYLKitrYVOp5O/YzabkZTkvtsuLi4OcXFxwVgaQRDDCE9a/dIbgSTlLL0BaK2/98WlK+3wsD9VERujx61zJgFgOFpYLa/T6mZAjFJqOjYmUs4LMMfnoSAoDuCtt96S//3QoUMoKirCb37zGyxbtgzV1dWYNGkSjhw5gg0bNgTj9gRBjBCUBl6v5xEbE+nyRiBJOZtbLDhaWB2wezuHktyRet0YTJ04FlW111QVR3fPn4xLV9ocjonD3fPTsTQ7Xd7lt1t65XAVx/kXbgokg9YHMGrUKOzYsQNbt25FT08PcnJysGLFisG6PUEQwxBnrX7nNwKpeUvaWZ8oNvVrMLsWOE78R1Bcvtx0DeWma+A4QMerK47cVfdIZGUkIELPw2oVwHPcyHoDULJ+/XqsX78eALBw4UIcPnw42LckCGIE4aySKb0RcBxw9FQNmMAQESHq9t+ZnYbC8/VytU8guWVmCjYsmY6Dn1Xg1IV61RsCY2JDV+L4aLczBtz9TY+szcLeQyUQGMO+j0phSI0b9DwAdQITBDFskN4I8oprcPRUjVxu2WsVZP2cYFJa2YTxsVEez3vayZcZmx1jHpkcCmq39EJgDIwBVtvQJILJARAEMeyov2rxWcc/UJQlpXodh+KLDSi6UO8IBXFgTkkCgcHtTr7M2Ixn3/hSTgwfL6rB9sduFRPBjkswNjSJYHIABEGEDO52yupzNfj0tAk2b2I9gYKJk7oSx0fLCWaBwdHW69752Nzs5Esrm1RVQXY7k6uZpC5hDkOTCCYHQBBESOBpp6xsnFImeJ0lFgKNUj8ot8AIjqvxWRqkLPWUkBO+jr+L48TdviE1DhERvEoHaLChiWAEQYQE3nbKUuOUEo7nkJoQM+D7SvLQzghMnCmQW2DEvo9KfYacOK5PDkKJJFKXljxGdlr7PioFALy8ZTF+uPKGkdUJTBAE4S/OO2WdjpN3xVkZCS4yD4LANE3d8kZa8hh8/7YMtFt6ERsTifySKzhbYQZjogPKLagGz9fIks1eYe7DOEqROul7Uqjo/jtnDNk4SIAcAEEQQ4Szxk+mIR6/e3SxnAOYOnEc8oprcPCzCkfFTP/uM2VCLK6YLeixug5zv9zYgX0flco7cENqHC5cuirLNDCIjsZTvCk1IQaNzV0QHKWo7sI4Uu+CBAf3oaKhgBwAQRCDjieNH+mf3AIj3jh4rt9GX4LngczJ1+Gm65Nx8LPv5ONjx0SiraMXjKkTt31lpiYcL6qG3S5OEJs97Tp8Xd7kcv2505MwdeJYWbvf3W7eWa7irnlpLgnuoYIcAEEQg447jR/lAHexQSoAN2LA0UIj9HoeG5ZMQ2llEyprW0XjDzFu77wblxzB1Ilj5UatcxWuxj9Cz2PqxLHyQJcLl666beZy7mYOBcMvQQ6AIIhBx1njR2mASyubIGgR4tGA5ERsNgGjoyOwYFYqvjNdk0sv505PxKLZE+Rks9I4t1t6wRyNWsrVcBxw45R4pCfHoaq21aMjU+KtK3goIQdAEMSg425XLOUEYmMioeM5zYPVlbibBwCod/l6SYOH51Q7eGe5aclJOWsLMSZq/F+oaoZOx6lE4PoT1/c270DL+YFADoAgiCFBuSt2zgkkjo/2u8Jn5tR43HFTGt44VKIq2ZyRPk5VnqnU4Png8+9UbwkHP6tAT69dnsX78pbF+OvfL3icMWC3M8yblYwZ6eP7ZaB9zTsI9uAYcgAEQQwZzuMcGURdn/6Ud35zqRkPrZqJdTkZqoTv3fMnqxzN8VPVsDusvmomMOubCSxNEDOkxqHCdE3+jo7n5N9KjI8d1e+Rjt5yIVrODxRyAARBDAnK3S3HcwPu6mUMsjy0UmKhqvaaPKrxzQ9KPIaWnI/ml1xBu6VXLuHkOWDZLekAOBw7ZYQgiL0KS7PT+71mb7kQLecHCjkAgiCGBOXulutHvN8dLe09aGnvUU3bOl4kNnJxnOvu3RuLZk+AITVOZYCl8s2l2Wlu4/L+xut9VQgFu4KIHABBEAHBX+Mn7W6dJR6UpCbEuISDvjcjARcuNfe9OQgMAhN342fKGl3GMco7fsbAc/BYXhoZwWPNrVNRVdsq5wAAuDXA7qp6+huv91UhFMwKInIABEH4REulilbjl1tgRH7JFUydOBbT08bhm6pmj+Gf2JhINHAWldEurbyK7BuSAQDjY6MwdeJYtFt6YW7pQm6B0evfMTFxNEyNnS7HJR0fyegr8WWA+/IYXUGN1wcDcgAEQXhFi3H3layUpJxrGtrlihop0eoJjhObrZydg83O5GQtIGr1b143G+aWLo9loBLujD8Ajzo+vsgtMMoVRTp+4CWhnghWKSg5AIIgvKKlEsXd8PYDJ8plI/jM61/A3o84v6fySyU2O8MbB8VpYBzPgROYpoHugJgkBuBRx8cbUseylFewCwy3zEyWy0gDZaiDWQpKDoAgCK9oqUSRZtxKoR1lc9VN1yf5NP4cAMOEWFy60i4f82TE3emyySEiO+vT+dHwt61wDHzpz866tLJJ1W/AAThT1gi73bMshD8MRmiJHABBEF7RUokiSR7bbAJKvmuSZ93abAJqzR0+7zEpaQyS40erHIAnJJMr7d51OnE8ozQkrLVDDOXEROlh6baJ3+VE2YepE8fiw5OVsnqnNPClP2RlJCAioq+reN6NySi6UB8QQ63c9fNBDC2RAyAIwie+EqGllU0KvXtRQZMDg17PY2LiGJgavDuBy40dqG+2+DXliwFYMCsFUZE6fP5Vrct5yfgDYgPXD5ZnItMQjwWzUgMST3d2jADw1beNAanZV4bdIDAsX5COxPExlAMgCCL0qKlvU30WBIaxoyMxKXkMoiJ1Pn/PANjtAq4bG4Wm1m7VuZlT41FmbHFbw3/qQr2meL+yG9hTCWd/nILztQJVs+8cdguWfDQ5AIIg+oU0wL2lvRunLtS7nG/t7EWrhiSuhMCAWRnXuezm77gpDQ+tmom84hqUVTf7nSeIcMz29fZ3BCrJGqia/cGSkCYHQBAEAP92wc4D3APFP85dwR03TcTX5Wa0dvSC48T5uY+szcKJYpNjqArnNanMc8DyBQa0tHejua1btft3R7D1dvrLYEhIkwMgCELzLlgl3hZg4w+I6pr/OHulT6zNYZTzS67IRpoXGNKS1XmFpPHRaGnvgd2RKBUrkWpgswkw1rV5rcgJtt5OKEMOgCAITbtgpZPQ6Xi/EraAqJoZoefR2NLl9XvOsX6OE+8kHdXpxMTyFXOn/N2rrV2Yd2MKxseOwtLsdL929aE8sSvYkAMgCELTLlhpVIV+7P5b2nvcHv/ejAScLW9ycSaSg7ELTDWP124XUOSUc7ALopRzZISYMPV3Vx+qE7uCDTkAgiB87oLLjM0wt1jE3XiAxjVKdHbb3Eo4eLqL9ILAQSzvFIS+twNpt3//nTPCdlfvD+QACGKEozW5K52T5uMCor5+S3sPTl9s6JeUgxbi46LA8xwEjdfX8aKz0Ot5PLI2C1W11/DpaZMc/5d2++G6q/cHcgAEMYLxp8RRWdnDO9ps/ZDP7zc3Z4rKnkqBN09Iwm/tll6VQ5Pi/rTb9w9yAAQxgvEnGZpXbJIre4Jp+JXJYw6iCuf42ChNvxUEhnZLr8sIRmm3X2ZsVonQSUPmnR0GIUIOgCBGMP4lQ4Nj9TmIej2AqNzpfJfYmEhkZSTgWFG1KsykV/xGPublb3DWz+E4Dja7AMbENUREBH6o+nCHHABBjGCk5G5ecQ365NNEnHMDUyeOA1AdkPvyHDB/ZgoAsfwT4HC00OjyPQbgzQ9KsP2xW7HjsVvlnMP42CgszU6Dsa5N1NsXRH2hR9ZmeTTgyrcdZmdQuhqG0GryChWC6gBee+01HD16FBzH4b777sPDDz+M/Px8bN++HT09PVi5ciW2bdsWzCUQBAHIXbQnimvw8pbFAKDKDTyyNguH/1EZsPslXxeDr75tlJUy783JgE7Huy0ftdkZ8opNeOy+OQCgckqS5DJDX/jHE8q3HXdvAFqbvII1fCUUCZoDKCoqQmFhIQ4fPgybzYZ77rkHCxcuxLPPPot33nkHqamp2Lx5M06ePImcnJxgLYMgwh53eQAA8jGrTcDr758LaACotb1XVge1CwwfnKwE5/UOzG3COjYmUjXgPTYm0uMVlDMJpIHu/uYAgjl8JRQJmgOYP38+/va3v0Gv16OhoQF2ux1tbW2YPHky0tJEYaY1a9YgNzeXHABBBBHnnXF5TQsAyBrzANAfUYfxsZFoaXe/I7f02FSfBYE5BaD60Os4j927AOQeAY7zPrZROZPgwqWreHnLYpdksfQ9Tzv8UNUFChZBDQFFRERg9+7d+M///E+sWLECjY2NSExMlM8nJSWhoaHB5XdtbW1oa1PLy9bX+y4RIwjC1cBJeYC//v0CLlQ1y+WWHAesWGDAmbIGn/IM7uiw2Hx/yQHPATwvirgp3wN4XizrlIysu4R1hMJ5mVu6UGZsdmuUSyubYLWJIR+rBjkLdzv8cNMFCnoS+IknnsAjjzyCLVu2wGg0OnQ9RBhjqs8Sb7/9Nvbs2RPspRHEiMOTgTPWtbnM12UMOFvRiKtt3R6u5h2rXft7ww1T4vHQqpkorWxCeU2L7ISYwFBVew2A525kKYn96WkTjhYa5TyGs3GPjYmUu4kZcx8u8rXDDzddoKA5gMrKSvT29uKGG25AdHQ0li1bhtzcXOh0fcMhzGYzkpKSXH67adMmrFu3TnWsvr4eGzduDNZyCWLYU2ZsxrtHy2C1CnLVS15xDUorm1BYWuf2N3VNlkFZW3pynKpWv/hig1wS+ulpkzzwxF33rpQMttu9h2baLb1yj4HUX+CM1vnGI93wSwTNAVy+fBm7d+/Gf//3fwMATpw4gQceeAA7d+5EdXU1Jk2ahCNHjmDDhg0uv42Li0NcXFywlkYQIw5p52+1CX0hFg44eqrGa/xdKzPSxqLc1Kr5+zwHcA55B0meWWrQyjTE4+75k/FJgRGAKO7mK9auxXBLM3oltVJzi8UlXKTc4cfGRMp5hnAx+M4EzQHk5OSgpKQE9957L3Q6HZYtW4ZVq1YhPj4eW7duRU9PD3JycrBixYpgLYEgRgRayhKl0IZSUE0QAKm5q78VPrExetw6ZxLOlLnm6rzBACyfPxmJ46MRGxMpJ2elsNTS7DScKK7xS63TV2imr+fBhONF1ThaWI0TxSY5XKR8jlkZCWFV7eOJoOYAtm7diq1bt6qOLVy4EIcPHw7mbQlixKC1LFHaIUvhn0DR2WVDboFR0zWnTIhFTX0HBMagc4gJZWUkIK/YJJeESslZX2qd7pyeltCMsnfAuZpIOctgSmqcU6jMFDZxfyXUCUwQIYzWskTn3a8tQMqd/mgCZU6+Do9tmKtK2H56ugZ2oS9ZrEzOejLouQVGsfuXMUT0Y3fuLlzkPMug3HQNgCRTweN4UTUEgYXd2wA5AIIIYTzFvp13yNLnpdlpmDpxLA59/h3qmjqDvj6pRp/nOUydOFbegUsduO46f33V8u89VCJP+vJUzukNT+Ei5zckjgPmTk9EynUxOFpYHTa1/0rIARBECOPOmDmPZrw5MwlnyhphswvgOA4cGPyo0Ow3UybE4p5FU+Xd+r6PSmFIjVOVYwIAz0v5CLHpy1u8XwrfyL/lvH/fU37E+e3C+Q1J2u3/YHkmgD6pjHCo/VdCDoAgQgRfxkySOja3WFThDKWOPgvwtC5vXG7sRFXtNTDG5OHtcvcu+soxl99ikFYnl3u6Q5w61gWdTmwY43kOW9b3NYkpnw8AF2PuK3QjPcel2Wluew0oB0AQxJDgK9nrvOvnec6heDm4pCWPgamhA4Bo8Gsa2sW3DsbU3bsRfWGrpdlpfmnw6HQ8VixMVzkLZZmru7ccLaEbpQPxNE8g3CAHQBAhgK9kr/I87AKWLzAAgIuGfjDhOeD7t2XgzQ9K5CauC1XN4AAXqWZ/a+2d/77E8TEuf78k8+D8lqNF6TPcRN60Qg6AIEIAX41OWRkJ8txcgQGfnTEhbnQkxo2JxNXWnqCvT+cIx6xYaEBVbavcxAWIoR7G1FLNknHVanTlMlbHDt9ZxiErIwE8x8GuMP6S4b9rXprP0NK7R8v6SlGt4ZXo9QY5AIIIAbw1OpUZm5FXbFLtfLt77eju9V/ArT/MSB+n2t1LTVxSRQ3H9e3AlWEWf5Q1JSlnafiLlFBW9gBsWT9bPq/F8AOK0JFVUYoK77LS4QQ5AIIIUcqMzTj4WQWKvmlQVcYMNlW1agkIZzkFSWtfOb0rIkIcMuPtrcY56d1u6ZWHv7jbpa9YaJA1/rUma+UOacUxX7LS4QQ5AIIIAZxj1OJu+NyglHP6QhCYW9VM57cUZf1+r1VAu6XX61tNf4a/+JusVYaW5NnAYVbq6Q1yAAQRAjiHSw7/o3LIjD/HATdOice31S1yuMWXwXSu3wdEA+7JYPsa/gIA+SVXVGEgb0hhMudSU09vKxT/FyEHQBAhgPPUrtrGjiFZhzSeo8J0DZvXzdZsMKUktV3hBLTM77Va+5K+htQ4RCi6dc9WmOXJXs5vD86Ncc++IZaIAsDxohpsf+xWvzSEwhV+qBdAEOGG1NBVZuwb0CLtVG+YEo9Ivc4vDZ5AIRl/qamr3dKL+++c4TPJeuBEOQBgy/rZ0PEcOACREWrZCnd/7yNrs8TKJkcXMSCWj86dkSi/CVit4kwD5f2e2/sl/uuTi3hu75eyM5CMPwDY7Ux+oyC8Q28ABDGIKHerEXoev3u0T6r49YNncelKe9DurezOVfqXGenjkDFxLKZOHKeSbfYV9pFF2xxJ35e3LMaOx2/1KFvhXArabul16SK+/84Z+MHyTIeekOvAGHeho6yMBPHNweEEdD7kJog+NDmAzs5O7Nq1C1VVVXjttdfwH//xH3j66acxevToYK+PIEYUecUm2VBZHTLEAPDMn75QhU+Cwfol0zA6OsJFn19Z4qmssgGgGuLiLMWgEm2z9hlwTw1szqWgnnofvA2McfebTEM8fvfoYrc5AMI7mhzAb3/7WyQlJeHq1asYNWoUOjo68MILL+D3v/99sNdHECMMZyPPcPCziqAbfwCwdFsxOjoChtQ4j9U5UrxcKcmsVwjO2e2i4b0zOw2Coi+B593vurMyEqDT8WAOiQfld7z1PngaGKMUdVM+S4rz9w9NDuDixYvYvn07Tp48iejoaOzatQurV68O9toIYtjha3rX0ux0HC+qgd3OoNNxWJqdjlf+qzjo6+I5MZQiGfCXtyx20cORcCfJrBScs1oFtLT3yGEXDsC8G5M93ltqYJP+P7fAiPySK1g0ewJWLDR4nW/g6VlKzkGa+AUgLMXcBoomB8Dz6lyx3W53OUYQ4Y4WvZlMQzy2P3arKpTSdC34Hb03TInHxUvNmrpy3ZV0KmEAii82YPO62aiqvYZPT5tQdKEeX33biEfWZqkqh6RrMYj9BAc/q5CdydflZgBig5c7tJaQ5hWbxM5kmwCe65OsUKJlrGY4oskBzJs3D6+88gq6u7vxj3/8A/v378f8+fODvTaCGFb4M71LFk1769SgVPykJ8ehvOaa21CMs3GUhquLJZricHfmWKS0VkFgqKptRf1Vizz8xWoTsPdQCRjrk2d2jtk3t3ar1pVfcsWjA3CHJBnNO9bEcRxa2rvlRi87Y9h7qETVP0BCcJ7R5AD+9V//FX/+858RGxuLV199Fbfddhsef/zxYK+NIIYVvgTdlM1KUyeOw5myBlVopb9wAG6ZlYKbM5Px+VcmXKhqdvlOTJRerLiBmFQ11rV5rNIBgDuz0yElVAHIjVRS8lgao2h3VOpwnDi8RdrtK6t6lKEcY12bPI4RABbNnqD571SulXeUmwoCw+lv1M/QuXPZH02icEOTAzh58iQef/xxldH/8MMPce+99wZrXQQx7PAl6PbL179QzOqtDth9dToON2cmo6q2FRcvuRp/DqKejyQbLTDIu2TXcEqNajqWVFHjXCVkbrHgaGG1XFY6d3oiFs2e4LaMVPl76f+VOQCtKNfKHH+LQz1aRYRTD4L0xgCNXc3hhFcHkJeXB5vNhp07dzrqdcWHbrPZ8Mc//pEcAEE44S5uXWZsFg1jkHT77XamSto6o9NxWDR7Akq+a5K/IzAmOyrlWwvAud0tOw9TKTM2qxzFD5ZnItMQr0msbcVCg1+GX8K5W1oQmGr0pDTjV1qL85CZ5QvSqUTUCa8O4OLFiygsLMTVq1fxt7/9re9Hej0eeuihYK+NIIYF3hKMkhHqtfZP2Gd8rCiI1trRKw5DQZ8OPmMMgiMOLrgZBenY9ILjOBhS4/rklBmTBdH6yiprAIiD3Z3DWJ5i6O7edoJZjuky19ch7sbzYtew0hEBvofMED4cgBT22b9/PzZu3DhYayKIYYNyVKFzBYo0iMTqwfgrhc88nV9zWwYA4L8+uShXvackjMb6O6bJu+3YmEh5SpcS6YVAaqS6/84Zbnfoxro2HDtVIzsG50qeAyfK3b4VDEXtvVRZZFeEgBg4LF+QjqkTx6mmj/nKyRAacwD3338/jh8/js7OTgBiGWhNTQ22bdsW1MURRCijNPAMYgXKn94/B0CMlSsFytzh1fhD3MEXnq9TCacxAPVXO7Hvo1K5lr/M2OxStslzgE7Hy3X/7uLx0t/gXPMvaQBJhJohdRaeE0PTnCr/oHxLcW4aI/rQ5AC2bdsGk8kEs9mMG2+8EefOnaMyUCKsUU6acjYtf3r/HKZMiPVq/H2ROD4ajS1dKK+5hvKaa9iwZBqqalvxdblZFkmTdrvvHi1ziYU/umGOV1kHCeeaf55z7ej11ZQ12KimgzneWgDmsdLn2Ckj7AKQW1CNx+6b06/8w0hFcyfwsWPH8Otf/xoPP/wwBEHAr3/96yAvjSBCF3eTppQMVNTN3KJuDiv9rgl33zJZbp5iAGrq27A/t0zeCUvD2ZVhKF+CbMqaf+m3vnoXQgHn6WAAVElp6ZgosyH+hgF4/f1zmmcMhAOaHEBSUhL0ej0MBgPKy8uxcuVKtLcHT7WQIEId5ZD2YOB81craVmTUtsp5Aw7Aya9rVTv/uTMSVUlQCW918KG2u/cHZ6fk7u9oblM3njGA+gAUaHIAMTEx+Pjjj5GZmYn33nsPU6dOhcViCfbaCCKkSU0YDVNDcAa3SBU8EnaBoaW9GxGOWDzHqYev6HjOrfEHfMfwQ21331/c/R13z5+M8ppr8mcdjyHPYYQSmhzACy+8gPfeew9PPvkk3n//fTz44IOUACbCCmcp5Gde/0KuRBkoUk27RFryGExMHOPSJXymrBE/vVes0JG6cpXVR0rpA+VOeDjv8geKFAo7XlSN+LgobFgyPaz+fl9ocgAHDx7EU089BQD4wx/+EMz1EETI4RxDn542bsDGn+eAdXeI+vydXVYc/Ow7+dz3b8uAITUOZ8oanSZd9VXolBmbVXINvnRvRsouvz/0t/EsHNAk6fn5558HeRkEEbo4x9AbrvY//JmWPAYLZqVg+QIDFsxKRVZGAqpqW+XzHCdOyso0xOOn92aB5/p+K4m4SUb+aKERJxwDZTytlUYjEt7Q9AYwadIk/PjHP8ZNN92kmgL28MMPB21hBBEqKGPoHAfYhf6Xd5oaOmBq6ADHAZ+ergFjTFXFE6GI0bdbelXJ4CkT4gD0b8oWQbhDkwMYN24cAKC2tjaYayGIkESKoUta9i3tvX79ngMQNyYSrR19v5Pm4AKQ1TSVOjYAEBsTqaryqai5huf2folH1mZ5NPLhHO8n/EeTA7jzzjtx1113qY59+OGHwVgPQYQcUlLVuaRQKwxAW4er09DpxPiO4FCpdK7iabf0qga4SzLL7ZZer0Y+nOP9hH8EVQ10z549+OSTTwAAOTk5eOqpp5Cfn4/t27ejp6cHK1eupGoiIqRRav14k27whbuf3j1fVKf0ZMiVTVrSW4JyEDoZeWKgBE0NND8/H1988QU++OADcByHn/zkJzhy5Ah27dqFd955B6mpqdi8eTNOnjyJnJycgPwxBBFoSiubBmz8dbw4VlWq6OEgatY7a+07owznxMZEqgTaCCIQBE0NNDExEc888wwiI0U524yMDBiNRkyePBlpaWkAgDVr1iA3N5ccABGSSMNElGEYf0hLHoNZU69zmarljyEfbjt9mr07vPDqAD766COsXbsWPT09eOutt1zOe6sCmj59uvzvRqMRn3zyCX74wx8iMTFRPp6UlISGhgaX37a1taGtrU11rL5+4KPzCEJCMlTOBjm3wIj8kisAGM5WNIExscu2P68A0aP0qhr9kW4Qafbu8MOrA6iuFsfWVVRU9PsGFRUV2Lx5M5566inodDoYjUb5HGPiMAtn3n77bezZs6ff9yQIbzjH9aWQzJpbp6oasiQ8TdryhVS142wIg7VLHurdN83eHX54dQBPPPEEAHEewKuvvopr1675dfEzZ87giSeewLPPPotVq1ahqKgIZrNZPm82m5GUlOTyu02bNmHdunWqY/X19TSUhggIznF9qbrm06IaTb/3FRIaOyYSbR1iDb/VKs7Z9dWpO1BCYfdNPQjDD01loL/61a+wfv163HjjjW537O6oq6vD448/jldffRULFy4EAMyZMweXLl1CdXU1Jk2ahCNHjmDDhg0uv42Li0NcXJwffwZBaMe5vh4AwAGtndrq+zkOiI7UwdJjdznHc8APV9wgT+hiAD49bZJDQf7ukrXu6kNh9009CMMPTQ5Ar9f73fX7l7/8BT09PdixY4d87IEHHsCOHTuwdetW9PT0ICcnBytWrPBvxQQxQJzr69OSx+Byo3ZVT4HBxfjrdBwyJ4/HQ6tmItMQj6raVnxSYATQN5LR3zGF/uzqQ2X3PdyS1uGOJgcwffp0fPvtt7j++us1X/j555/H888/7/bc4cOHNV+HIAKNPGLRUZZ5ubFjQGWeALDjsVtVhm9pdho+PV0Dm02QNXwA/3bJ/uzqafdN9AdNDsBkMmHDhg2YMGECRo0aJR//+OOPg7YwgggWmYZ43DUvDbkF1eJQcQ/GPzUhBnVNvoXfdLw6LFpmbEZesQmCQzOIOd1A6y7Z31097b4Jf9E8E5ggRhJTJ44Dz9d4rfDxZPw5ADdOjceFqmYAooGXdudlxmaXWQGCwPoVk6ddPRFsNDkAGgBPjCRyC4x44+A59LO6EwxiIjlCz8NuV+/O//r3Cy6zAgYSk6ddPRFMNDkAghgp5BYY8af3z/n1G3dln4Xn66HjgWlp43D3/Mny7l96K5BIGBuFp380j4w4EZKQAyBGLM4llGXGZuw9VOLXNXgOuGFKPMqMLS7hIrsAlNdcg7GuDYbUOOQVu/YR/PPd15PxJ0IWcgDEiEIp8bDvo1I5gbrm1qk4+dVlTV29aclj5GHvAgO+qWqGTsdhwawUAOJsXptNkN8KpGavlvYe1XWmTIijUYRESEMOgBgxKOvmOU4ctM4A9FoFtxIPnkgYG4WGZossw8wgJnJnpI+X5/HmFdfgWFEN7Ipmr5sz1V3tyfExOHCinBK4RMhCDoAY9ki7fnOLRa6bB2Pg+inj+XV5EzYsmQZLtw3Hi6rlgS3KWn7RoHOqZq/xsVHQ6zjY7Qy8jkPxxQYUXagnYTQiZCEHQAwr3MX1pV2/TseD4znAUYUzkOauqtpWvLh5EZZmp3ksw1yanYYTxTVymGlpdpr8fXOLBUcLq0kYjQhpyAEQwwZ30gjKblnYBUxLG4fymmser3Hd2FG42trj8bzEotkTAHgvw/RUpy85phPFpiGXZiAIb5ADIEIa5Y7fnTSCsltWp+MRoefBeZDv5wBcUwx05wDwPAfGmJworqptxaLZEzQnbz05iGA3cQ219DMxMiAHQIQszjv+R9ZmuZVGuDM7DS3tPTh9scGlDl8JAyAwJhp9R1z/p/dmBW3UYrCauEJB+pkYGZADIEKW0somuRLHahXQbulV7aoBqKp+fJV4cgD0Oh6C4KjuYQyG1LhhZzxDQfqZGBnwQ70AgvBEbEykXMQjyS9kGuJx/50zkGmIR16xCb1W0RAKgqPqxws8z+HmzCQw1lfaWVrZFOS/IvBIYS+eG5jMBEHQGwAxJLir5nGOabdbeuV4PscBZ8oacLyoGvFxUYiK1OHzr2rl6zEAqdd5V++0C0ws1QwB3fyBQCJxRKAgB0AMOu5i+/s+KoXVKoDjgPkzU7BhyXRxp6sTjTUHDoXn671et67J4jEBDIhORFmq6Ww8h1NilUTiiEBADoAYVMqMzXj3aJkqtn/8VHVf1y0ThdaKLzZg87rZYIzJ8XotePva+jumqUo1ndf1S4eMs07HYbvTgBeCGImQAyAGDWnnb1Xo6DAAlbXXXL5rszMc+vw7F2llrTg3AS+YlYKHVs/0+P284hrYHPey2ZlqkDtBjFTIARCDhlS94rxLtwvuv1/X1OnX9ZVGnwGqcs8NS6a7/Y0U9nEWchOvRhAjG3IAxKChbNrieQ4cx6lUNTkAk5LHwNJtQ3Nrt18yPrExeiy7xYBDn38HxgC9jpOF3JThI2WcH+grI+V5TtbxkWQdCGKkQw6ACBrOSVXn6hVADL18etoEu13s5K1r6pQVNv3Rcps2aRw+/qIKYKKGf2rCaHHYO9Tlnsrk853ZaX0yEgLD8gWTkTg+ZlgkgQkiEJADIIKCp25V5+qVTEM8lmano7SyCeU1LTh1vl42+ikJowEwTYPZvy7vq+dnDLKeP6eolXduoAI4VUno0ux0MvxEWEEOgAgK/nSrSsf/+9i3qh1/fVOnyxtAzCgdLD12t9eR5J+VIaW50xPxg+WZAABziwU6HQ/YXdU7addPhCPkAIigoIz3KxuucguMyC+54iK4lldsgtWmzga7C/9E6HWABwew/o5psHRb5ZCSXs/Lxl8Z61++wICl2WkeS0IJIlwgB0AEhUxDPB5ZmyUb+0xDvGog+9flZgDAioUG5BYYcbTQKP+Wc/yPc7VQakIMAKC1s1d1fHS0HlkZCVgwK1UVUpJ29QdOlKti/Ynjo8noEwTIARAa8dUl607aYd9HpbDaBJR8J8bn80uuqH5zvKgahtQ47D1UAqWOG5P/R427XIBOx6HXKqDoQj2++rZRlWuQ8PQ2QhDhDjkAQoU7Q+9LftjToBaro+bfzhj2HirB7GnXqe5VXnMNu9/7WtOgdnd8b0YiUq4bjaOFRq+5BtLOIQj3kAMII7Ts4p0NOQCVdIPNJiCv2KS6jtLYWxWDWniOg90Rx7ELDGcrXJU3pWodf4lQxPeVYxml3b27ElQy/AShhhxAmKBliIhz5U5ecQ1OOJKzUl2+TserBqW/vGWxKNvs2MQz1ifbvGX9bLxx8BwE5qjpH8CMXomk8dG4OTNZlcR13t3TwBSC0AY5gDBBS1lmbEwkOIecJsdxaGnvkaUbOE4sqUy5LsZl2DnQ17TFQZRxBgBDapx8vQDYfvA8hyd/mO02xKM8RgNTCEIbNBAmTPA1RERK2gpCX8jmTFkjeJ4Dz/WFXJZmp7tcJzYmUpR2ABAR0Xftg59V9Du+7wzHAY+un63JkAdrYEqZsRkHTpSjzOh57CRBDCfoDSBM8JUIlYXaFMfsdgHLFxiQOD5a9RvpOrExkbKUgzRr95G1WXIY5tQF7/r9Et40/CVumZni16D2QCd9KaxEjETIAYQR3hKh0q5ZSuZygNwt6y7kAojNVVJyGBAHrlc5pJ3zik2aYv463vcs3wgvap6eCHTSl8JKxEiEHAABQL1rjo2JRLul1+vu2d0bAwNwvKgGS7PT0dLe7fOeqQneRzgCYtJXmgE8lFAvATESIQcQJmgZd+i8a5Zi3u7m9ioNojRkHRCHqRz8rAKnv2nwuSYtIm+NLV3484elMKTGDakToF4CYiQSVAfQ0dGBBx54AHv37sWkSZOQn5+P7du3o6enBytXrsS2bduCeXvCgZZGrrxiEwAmK2I6/2bNrVPx4clKCII4MjH7hmTcdH0S2i29uFClTorWmjsClvwFQifkQr0ExEgjaA7g3LlzeP7552E0GgEA3d3dePbZZ/HOO+8gNTUVmzdvxsmTJ5GTkxOsJRAOvMWvy4zNePaNL2UhtqOF1Xh0wxy0W3rl31itgjxoBRB3+d4GtOt1gZ2mxes4CrkQRBAIWhnoe++9h1/96ldISkoCAJSUlGDy5MlIS0uDXq/HmjVrkJub6/a3bW1tuHz5suqf+nptFSWEiLJk0VtZpNTFKyEwYO+hEsTGRMq/4XlO+2QWANX1/evu5QDoeNdhjMvmT6adN0EEgaC9Abz88suqz42NjUhMTJQ/JyUloaHBfZz47bffxp49e4K1tBGPu5CPp/h1VkYCIhzVPxJSNc+d2ekAGGKiIvDhyUrYBaapZFMYQPhn2S0GAEwl6UzjGQkiOAxaElgQBLEr1AFzdJu6Y9OmTVi3bp3qWH19PTZu3BjUNY4USiub5PJMq1XU7nGu5ZfINMTjd48uxsHPKlD0TQMYY+A54OipGggCk98ABI3GXwnPAeAAwcPQdwCIGx2JDksvGBObyKSyU2dJZ4IgAs+gOYCUlBSYzWb5s9lslsNDzsTFxSEuLm6wljascVfdExsTKUdsGIBjp4zioHQPDUyZhng89/AtyC0w4vD/q3TM0hWvIDBAsLO+i/mBwIAIHY8bDONcEsUS7Z290Ok43D1/ssuQFjL8BBFcBs0BzJkzB5cuXUJ1dTUmTZqEI0eOYMOGDYN1+xGHVLnjLMwGuOruCwJkJU9P1TRlxma8+UEJbPbAVe8AgM0u4ObMZMTGRLpNHEtD22lIC0EMPoPmAEaNGoUdO3Zg69at6OnpQU5ODlasWDFYtx9RSDH+XmtfbMVZvRPo6+ZljMlOwrmaRhrROCpSB/sAjD/HifdzDv8zBpwpa8AdN6Wh6EK96ryO58AY8yrjTBBE8Ai6A8jLy5P/feHChTh8+HCwbzlikYyjucUCmyJpKxl6gHNR75Q0890Z1d/vL8bnX9WqrtNfDKlxqK5rc3vuQlUzLl5qVkWQIvQ8fnpvlqrjOFB6O+RECEIb1Ak8TFAaR52OF0szBQadjsdd89KwNDsdgHo4yg+WZ3ocfJ5bYFQZf8DvEL+KS1fcG38J5c6fA3DXvDQXcbdA6O2QaBtBaIccwDBBaRzhQaUTcB2O4onjRdX9XktqQgwarlpcwj1akCSjJYelJBB6OyTaRhDaIQcwTHA2ju5UOgHv1TNlxmYc/KwCtY0d6Ojq7fda6pos/QoX6d1U+ygJhN4OibYRhHbIAYQwzrHsgRjHMmMzfvn6FwGr8vHnKqLhT5d1hrwx0PJPEm0jCO2QAwhRPMWy+2vQ8opNmoz/HTdNxJmyRrRbrP26j5KZU+ORnhzncccfLKiHgCC0QSMhQxR3sez+UmZs1hzz/39f1wbE+ANAenIsHrtvDgDQKEWCCEHoDSBECWQsu7SySbM+TwBVnAFwVJVDECEMOYAhwleteiBj2ZIz6bUK4ACMi41ES3v/k8Ce0PGA3dGeoNdxWJqdRlU5BBHCkAMYAnILjNh7qASCwBAR4XlXHKhYdqYhHmtunYqDn30HBqClvRc8DzBBbCBbPDvVpSfAXxbMSsGGJdORV1wDgFPF/fvzJkPNXAQRfMgBDDJlxmbsPVQiT8yyWgdnV1xV26r6PDFxDJbcnIasjAT89e8XBnRtDpCHtieOj1EZ7f68yVDYiCAGB3IAg0xpZRMEhaYyzwdv2pVyF71o9gR8Xd6nxnq5oQOdXVaUVjah8nKrl6v4JnqU+J+RJ6Pt75sMhY0IYnAgBzDIKAew8ByHLetnD8i4iaqg6rCLdOx4UQ1sdoYIPY/fPboYC2alyIqcDMDBz74LyN80f2ZyQI02NXMRxOBADmCQCWRy17m5K7fQiJzvTUR+aZ1KKdRqE7D7va8xMXEMHBJCA4LnOYyPjURXtx3zZybjFxuzUWZsDpjRpmYughgcyAEMAYFK7pZWNqkknBmDx2SuqaEDpob+zeoFgBnp45CVkYDR0REeJ4sF0mhTMxdBBB9yAMOYrIwE6HScS4evv6MbfbFhyTQ8tHomgL68AuCqMOrOaEvfj42JVEk/EwQx9JADGEY4l0ZmGuKxed1s/On9c6rv5XxvIrp77S4DWPrDHTdNVBl/f6pzpO9bpRkFgNeyV4IgBhdyAMMET8a33eLa0FVd3w5jXduA3gJGRfCYljYO3b12vP7+WXlIuz+JXun70jp8jaUkCGJwIQcwDCgzNuPdo2WwWgUwiL0DecU1yDTEIzYm0uX7voaz+CI1IQaNLV2qQe7Hi2qwed1svxK9UjWP8g2AqnoIInQgBxDiuJNxZgByC6txvvIqJiaN8Tvm700OguNEvX9nbHaG/JIreGRtluZYvjIxTDkAggg9yAGEAN5kD/KKa9zKODMGmBo7YGrsgN5NItgbERE8nn3oFhSer0N+aR30PHCtoxcJ46JQXdfuUev/bIUZFy5d9SuGT9U8BBG6kAMYYjzF9iWnUNPQ7vMa3pQ+eY7DnOnXIWtaossuPNMQjwWzUuX7d/faVbOGb85MAgA0t3WjwnQNjDpzCWJEQQ5giPGk+//sG2L1jBa8VfoIjOHCpWb8YPkNbks0lbkFT7OGnZ0UxfAJYmRADmCIUcoe6HQ8zC0WHPysQrPx14LSsUihJgB9JZroS9C6m95FnbkEMTIhBzBIeIrzS8b14GcVKLpQj9yCanD9mbjuAcmwx8ZEqnbxd2anyyWaHAfMnZ6IHyzP9GjcKZZPECMPcgCDgLs4P6DejZ/+pkEO5Xiq6Bk7OhKtndoGucxIG4e7b5ksx/ydQ00AU5V0ejP+BEGMTMgBDALOxjevuAYnik2K3XiaPB/AGzdMicdX3zaqhN4klKWgOh2HjEljYUiNUxl1pcFfmp0uN3dRWIcgwhNyAINAbEwkOI4D59h1A5zcHGW1CThbYfZ5DYnk8dEwNXaqjs2cGo+HVs1EXrEJLe3dKL7YgKOF1ThRbJKrijzF8cnwE0T4Qg4gQEgx/s4uK6pqW7Fo9gQYUuOQV2zC8aJqCAIDz3NYc+tU1Jo7+uQRGNDc1qPpHpKWvzPpyXGykT9wolzWAHIu2aQ4PkEQSsgBBABZ9EwqpwTwdblZNSQdEOv1PzxZ6RLuGROtR0+vvV/35jlgaXaa/JmGqRAEoRVyAAFAFj1zOm53CtVznPumraut2t4AnImN0eOF/2+h26oiiu0TBOELcgABQIrxO5fv8I5D0lGO5yD4IdngDZ2OczH+EhTqIQhCC+QABkiZsRn7PiqFIDDoeA4ZE+NQbuobsp44PhqNLV0AEBDjf8dNE5GeEke7e4IgBgw5AI0ok7yf5F9CV68d0yeNxYKsCXL4RxCYyvgLDLLxB+BRZE0LSeOjcf+dM7BioWEAVyEIguiDHIAHlJ27gCib4Fx/X25qlQ0+h8CPYuR5YNqkcbh7/mQy/ARBBJywcgDK+bRVtdcAcHIFjdLY//XvF/DNpWZZJuHGKfE+tXkYAENqLC5d8a3eqZXltxjw2H1zAnY9giAIJUPiAD7++GO88cYbsNls2LRpEzZu3Bj0e7or1QSAY6eMYqKWiZ2ydrugUtdkDKrJWN4IpPGPjOBV5Z0EQRCBZtAdQENDA1599VUcOnQIkZGReOCBB3DLLbdg2rRpQbunJHvsTkJBWaoZSAXOgfC9Gd6F2QiCIALBoDuA/Px8LFiwAOPGjQMALF++HLm5ufjZz34mf6etrQ1tbeq5tvX17rtgfSHt/N0Z/1BEr+PI+BMEMSgMugNobGxEYmKi/DkpKQklJSWq77z99tvYs2dPQO4nNWmFKjwHPLphDqpqWwEwLM1OJ+NPEMSgMOgOQBAEsWnKAWNM9RkANm3ahHXr1qmO1dfX9ytXIEkjhOIbwIJZKdiwZDoZfIIghoRBdwApKSkoLi6WP5vNZiQlJam+ExcXh7i4uIDcTymN0NllRel3TYgfG4WbM5PRbukVj1U2IcIxNGV8bBSmThyLdksvaurbUFp5FZERPFrbe2Dp6Z9ej8SGJdOQct1o5JdcwaLZE6i0kyCIIWXQHcCiRYvwxz/+Ec3NzYiOjsaxY8fw0ksvBfWeoSaNQIafIIhQYNAdQHJyMrZt24Yf/ehHsFqtuO+++zB79uzBXgZBEETYMyR9AGvWrMGaNWuG4tYEQRCEA36oF0AQBEEMDeQACIIgwhRyAARBEGHKsBGDs9vFEsz+dgQTBEGEI5LNlGyokmHjAMxmMwAMinAcQRDESMNsNmPy5MmqYxxjgVSwDx7d3d04f/48EhMTodPp/Pqt1EW8f/9+pKSkBGmFwWG4rp3WPfgM17UP13UDw2PtdrsdZrMZs2bNQlRUlOrcsHkDiIqKQnZ29oCukZKSgkmTJgVoRYPLcF07rXvwGa5rH67rBkJ/7c47fwlKAhMEQYQp5AAIgiDCFHIABEEQYUpYOIC4uDj87Gc/C5jC6GAyXNdO6x58huvah+u6geG9dmAYVQERBEEQgSUs3gAIgiAIV8gBEARBhCkjzgF8/PHHuOeee7Bs2TLs37/f5fzFixexfv16LF++HM899xxsNtsQrNIVX+ves2cPlixZgrVr12Lt2rVuvzNUdHR0YPXq1bh8+bLLuVB93hLe1h6qz3zPnj1YtWoVVq1ahZ07d7qcD+Vn7mvtofrMX3vtNdxzzz1YtWoV3nrrLZfzofzMvcJGEPX19WzJkiWspaWFdXZ2sjVr1rCKigrVd1atWsW+/vprxhhjv/zlL9n+/fuHYKVqtKx78+bN7KuvvhqiFXrm7NmzbPXq1WzmzJnMZDK5nA/F5y3ha+2h+My//PJL9s///M+sp6eH9fb2sh/96Efs2LFjqu+E6jPXsvZQfOanTp1iDzzwALNarayrq4stWbKEVVZWqr4Tqs/cFyPqDSA/Px8LFizAuHHjEBMTg+XLlyM3N1c+X1tbi+7ubsydOxcAsH79etX5ocLXugHg/PnzePPNN7FmzRq8+OKL6OnpGaLVqnnvvffwq1/9ymWuMxC6z1vC29qB0HzmiYmJeOaZZxAZGYmIiAhkZGTgypUr8vlQfua+1g6E5jOfP38+/va3v0Gv1+Pq1auw2+2IiYmRz4fyM/fFiHIAjY2NSExMlD8nJSWhoaHB4/nExETV+aHC17o7Oztxww034Mknn8QHH3yAtrY2vP7660OxVBdefvlljxIdofq8JbytPVSf+fTp02VDYzQa8cknnyAnJ0c+H8rP3NfaQ/WZA0BERAR2796NVatWYeHChUhOTpbPhfIz98WIcgCCIIDjOPkzY0z12df5ocLXukaPHo19+/YhIyMDer0eP/7xj3Hy5MmhWKpfhOrz1kKoP/OKigr8+Mc/xlNPPQWDwSAfHw7P3NPaQ/2ZP/HEEygoKEBdXR3ee+89+fhweOaeGFEOICUlRZaNBkT5U+XrvfP5pqYmj6//g4mvdV+5cgXvv/++/JkxBr0+9HX8QvV5ayGUn/mZM2fw0EMP4Re/+AXWrVunOhfqz9zb2kP1mVdWVuLixYsAgOjoaCxbtgzffvutfD7Un7k3RpQDWLRoEQoKCtDc3Iyuri4cO3YMt99+u3x+4sSJGDVqFM6cOQMA+Oijj1Tnhwpf646KisIrr7wCk8kExhj279+Pu+++ewhXrI1Qfd5aCNVnXldXh8cffxy7du3CqlWrXM6H8jP3tfZQfeaXL1/G888/j97eXvT29uLEiRO4+eab5fOh/Mx9MjS55+Bx+PBhtmrVKrZs2TL25z//mTHG2E9+8hNWUlLCGGPs4sWLbMOGDWz58uXsX/7lX1hPT89QLlfG17pzc3Pl888880zIrFtiyZIlciXNcHjeSjytPRSf+UsvvcTmzp3Lvv/978v/vPvuu8PimWtZeyg+c8YY2717N1u5ciVbvXo12717N2Ns+P137g6SgiAIgghTRlQIiCAIgtAOOQCCIIgwhRwAQRBEmEIOgCAIIkwhB0AQBBGmkAMgwo5nnnkGf/nLX4J+nwcffNBvTRiTyYStW7cGaUUEoYYcAEGEEFeuXMGlS5eGehlEmEB9AMSI5dSpU9i1axcmTJiAqqoqREVFYceOHdi3bx86OjrQ1NSEpqYmTJ8+Hb///e8RExOD4uJi7Ny5E11dXYiIiMDPf/5z3H777TCbzXj66afR0tICAMjJycHPf/5zHDp0CLm5uRAEAVeuXEFycjJ27NiB5ORkPPjgg0hJSYHRaERTUxMWLlyI3/72t+B5Hnv37sWJEyfQ3d2Nrq4uPP3001i6dClWrFiBhoYGzJs3D3/5y1/w1VdfYdeuXejq6gLP8/jZz36GJUuWDPGTJUYMQ9uHRhDBo7CwkGVmZrLTp08zxhh799132bp169jTTz/N7rvvPmaxWJjNZmPr1q1jH3zwAWtubmYLFy5kZ8+eZYwxVl5ezubPn89qamrYnj172L/9278xxhjr7OxkP//5z1lbWxs7ePAgmzt3LquqqmKMMfbKK6+wrVu3MsYY++EPf8geffRRZrPZmMViYYsXL2anT59mly9fZg8++CDr6upijDF25MgRtnr1annNq1atYowxdu3aNbZs2TK5S7m+vp7dfvvtrLa2dpCeIDHSGXqlJYIIIpmZmbLk84YNG/Diiy8iKSkJd911F6KjowGIMsXNzc0oKSlBeno65syZIx+/6aabUFRUhNtuuw0//elPUVdXh0WLFuEXv/gFYmNjAQCLFy/GlClTAAD/9E//hLVr18r3v+eee6DT6RAdHQ2DwYCrV68iOzsbO3fuxMcff4zq6mqcO3cOnZ2dLms/e/YszGYzHn/8cfkYx3H49ttvMWHChOA8MCKsIAdAjGh0Op3LMZ7nVSqTHMeBMQa73e4i48sYg81mw+zZs3HixAkUFBSgsLAQ999/P/bt2+dyD0EQVJ/d3efChQt47LHH8NBDD2Hx4sWYN28efvOb37is0263IyMjAwcOHJCPNTQ0ID4+vh9PgiBcoSQwMaIpKytDWVkZAOB///d/8b3vfQ9xcXFuvzt37lxUVVWhpKQEgKhbf/r0acyfPx+7du3C66+/jrvuugvPPfccpk2bhoqKCgBAYWGhPADkf/7nf3zG6E+fPo1Zs2bh4Ycfxvz583HixAnY7XYAojOxWq3yeqqrq3H69GkA4tzZ5cuXD5thI0ToQ28AxIgmISEBf/jDH1BbW4v4+Hjs3LkTe/bscfvd+Ph4vPbaa3jppZfQ3d0NjuOwfft2TJkyBZs2bcIzzzyD1atXIzIyEtdffz1WrVqFI0eOIDk5GU8++STMZjOmTZuGF1980euaVq9ejWPHjmHlypUQBAFLlixBa2srOjo6MG3aNIwaNQr33XcfDhw4gN27d2Pnzp3o6ekBYww7d+7EpEmTgvGoiDCEqoCIEcupU6fw0ksv4ciRI0G7x6FDh3D06FG8+eabQbsHQQQLCgERBEGEKfQGQBAEEabQGwBBEESYQg6AIAgiTCEHQBAEEaaQAyAIgghTyAEQBEGEKeQACIIgwpT/HxjBol0PwT3DAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(x,y,'.')\n",
    "plt.xlabel('phosphate')\n",
    "plt.ylabel('nitrate')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Create a subset where both variables have finite values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "ii = np.isfinite(x+y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0       True\n",
       "1       True\n",
       "2       True\n",
       "3       True\n",
       "4       True\n",
       "        ... \n",
       "2343    True\n",
       "2344    True\n",
       "2345    True\n",
       "2346    True\n",
       "2347    True\n",
       "Length: 2348, dtype: bool"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ii"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Method 1: Scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = stats.linregress(x[ii],y[ii])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinregressResult(slope=14.740034517902107, intercept=-3.9325720551997954, rvalue=0.9860645445968036, pvalue=0.0, stderr=0.05292378356970245, intercept_stderr=0.10258209230911518)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.740034517902107"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.slope"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exercise: Draw the regression line with the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f8be85239d0>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD7CAYAAAB37B+tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABAW0lEQVR4nO2dfXhTdZr3vycnaUE0Y9RisUWcRdgKRVGzSp1La1WgFJGBtlP6MsOMIogO7so+Kray3Q5bBHWdR8taZUe9kNIaG2VFn1pAq4wLsRi0YIAqIi+2UIjasaNC2iTn+aP+fp5zcvLatE3a+3NdXJLz+stRv+fO/SpIkiSBIAiCiFt0g70AgiAIom+QkBMEQcQ5JOQEQRBxDgk5QRBEnENCThAEEefoB/qGZ8+ehcPhQFJSEkRRHOjbEwRBxCUejwdOpxPp6ekYMWKEYt+AC7nD4UBxcfFA35YgCGJIsGnTJpjNZsW2ARfypKQkvpjk5OSBvj1BEERc0tHRgeLiYq6hcgZcyJk7JTk5GampqQN9e4IgiLhGyyVNwU6CIIg4h4ScIAgiziEhJwiCiHNIyAmCIOIcEnKCIIg4h4ScIIhhhd1uR1VVFex2+2AvJWoMePohQRDEYGC322G1WmGxWODxeGAwGGCxWHyKa/yda7PZkJGREdLxAw0JOUEQQx673Y6CggK4XC7IZ+nYbLagwszO7enpCUv8BxJyrRAEMeSx2Wzo6enhIi4IAgwGAzIyMoKea7Va4XK54PF40NPTA5vNFvE6PC5XxOcGgoScIIghT0ZGBgwGA0RRREJCAkpKSkKyrO12OywWi+IF8Nprr+Guu+4Ky8d+tqMDO+fm4sPfFOFvLXv79F20INcKQRBDHrPZDIvFErKfm/nEW1pa0NPTw7e73W4cOnQIhw4dwjvvvIPKykp0dnbCZDKhs7PT59qSJOGjhXei57suvu0XU9Kj/v1IyAmCGBaYzeaQfNs1NTUoKyuDx+NBoNn0brcbDz/8MP8sCAISExNRUVGBzs5OmE854f10P99/+R+X4uLpt/XtS/iBhJwgiCFDKNklLHvF6XQiKSkJeXl5/NiamhqsWLEioID7Q5IknD17Fv9e9ig2TL0WXtm+62tfhn7UqEi+UkiQkBMEMSQIJbvEbrcjPz8f3d3dfNsrr7yCBQsWID09HaWlpRGJOONV8zTF5x/PPRfTN22I+HqhQkJOEMSQgGWmeDwe/lkt5OwYOT09PaipqYEoivzccJlqPB+lE9MU24o+3o3XNm+O6HrhQkJOEMSQgGWmAPCbWsiOkVvkQK9bJFIRV1vhAPAb+4cAgNbW1gHJOSchJwhiSKCVmaLlMy8oKIDT6cThw4dx6NAhfn64LpXnr7wGpoQExTYm4IyGhgaUlJRE+I1Ch4ScIIi4REuk5Zkpap95RUUFVq5ciZ6eHoiiiGuvvTbie6ut8F3ffoP/++Uhn+NycnIivkc4kJATBBF3hBLYlPvMvV4vqqqquEvF7Xajubk57PsGcqMAwNSpU3HDDTdg//79yMnJGRBrHCAhJwgiDlGLtNVq9RHyjIwMHsCUJAltbW0R32+kTsSGa/5Jse2Nkyewqf04/5yQkICKiopB6cNCJfoEQcQcNTU1KCoqQk1Njc8+u92O9vZ26HS98iVJEurq6jRL5kePHt3ntbxqnuYj4r+xf6gQ8VGjRqG+vn7QmmmRRU4QRExRU1PDKyZ37NgBANxFIXepyIOTbrcb5eXl3CK22+3Izc2F2+2OeB2ZFybhvl+OV2x75MCnOPzjDz7HGgwGTRGvqalBQ0NDv7tZSMgJgogpGhoafD6npaXBZrOhvb2du1QEQVAc19LSgtzcXN7/pC8iHswXrua7775DTU2Not9KoBdStCEhJwgipsjJyeHCBwCTJ0/mVrgoihBFEQAgiiIkSfJpalVWVobKysqI7q0l4AX2DxEsMVGSJJSVlUGSJB581XohDbqQr127Fp2dnVizZg127dqFxx57DC6XC7NmzcIDDzzQL4sjCGL4UVJSgmPHjuHtt9/GrFmzYDQaFRWbRUVFSElJ4QU/1dXVaGxs5Od7PB689957SE1NDSvAqSXihZ98FFTEgd6Xitfrhdfb22HFZrP5vJD6MxUxJCG32WzYvHkzbr75Zpw9exalpaXYuHEjxowZgyVLlmDHjh3IzMzst0USBDF8sNvtePHFF9Hd3Y3nn38eZrOZW+EGg0HR5EoLSZIUwh4MLQE/kj8PJSUlmH7XXX6vZTAYcOuttyIpKQnp6ekoLy/n6ZDy3PaY8JH/7W9/w5///Gfcc889aG1txb59+zBu3DiMHTsWADBnzhw0NjZqCnlXVxe6uroU2zo6OqK0dIIghiI2mw3d3d3cum1uboYoipgyZQpuuOEGWK1WVFdX8+O3bt0a8b0CibjdbkdTUxPfLggCUlJSkJKSgokTJ/IXCitMYu1r5SJeUlISG5Wd//Zv/4YHHngAJ0+eBACcPn0aSUlJfP/o0aNx6tQpzXM3bNiAdevWRWmpBEHEM8FazLL9JpMJOp2OCznQ6y5paWlBS0tL2PdVXwsIHMz87RUTUFVVhfb2dkXAVJIknDhxAl9//TVKS0u5iMfCPM+AQl5fX48xY8YgIyMDr7/+OgDA6/UqosWSJPlEjxkLFy7EvHnzFNs6OjpQXFzc13UTBBFHBBM89f477riDa05fmDRpEj777DP+edzIc/DE5CsVx7x0/CjePt3rKdDr9bBYLPB4PNydI8fr9fK5nWazWfHrQZKkkIY59wcBhbyhoQFOpxNz587Fd999hx9//BHt7e2KL+h0Ov0m3RuNRhiNxuiumCCIuEOrxSz7Z0ZGhmK/JEnYsmULBEHoU29wADhw4AAvHApkhZ933nn41a9+haSkJNTW1vJ1jh8/XtFYSz202WQycWvf6/XCZDL1ab2RElDIX3rpJf73119/Hbt370ZFRQVmzJiBY8eOITU1FW+99RZyc3P7faEEQcQv6hazJpPJp6EV2y8IArdwo8FD/zAB15yvFNi7Wuz4u8xtctFFFyErKwsOh0MRWF20aBFvtKXX67FgwQJFsLWzs5O7bgRBgMPhiMqawyXsPPLExESsWbMGy5Ytg8vlQmZmJrKzs/tjbQRBDBHULWbVFrrD4UB+fj6AXut4/fr1Pn7tSNCywhd8vNtn25EjR/Dwww9Dp9NBr9ejqKiICzYrRtLy7WdkZECv16O7uxuSJMFisQTNqukPQhby+fPnY/78+QB6F79ly5Z+WxRBEEMP9fBjuQVeW1sLr9fL3SCRDnlgBHKjZGdnY+nSpVi9erVPB0Sv1wuPx4OUlBTN1riAb9C2oKAAGzduBNBbkDQYfnJqmkUQxIDDLPSioiI+nYf9sz9FnGGz2TBx4kT+4lDjz9ddU1OD3NxcPP744ygoKIDdbkd6ejrfP1h+cirRJwhi0Dh+/HifhZvhT8DlYm0wGNDU1ITt27f7zbbzeDwoLy9HWlqajyVeVlbGUxK7u7t50Jb5yXU6HTo7O6PyfcKBhJwgiKhjt9thtVoBwMdnzPZZLJY+NbZijBRFbLha2Wb2wN+78O+fHeCff/vb3yIlJQUtLS3YunVr0ECqPMWQYbPZFH57nU7Hs1cSEhIUVZ0DDQk5QRBRxW63Iz8/n0/jsVgsvFc3yxd3uVxcTPuSZhiKG4VloZhMJjQ1NYV0L1EUfQSZBTZ7enogCAJuu+02ANqzQgcaEnKCIKIKy0hhqAtozp4963NOuGI+Peli3D3ul4pt/36oFa1/7/K5Xk9PDzZt2gSdTheSG0cQBBQUFPgVZEmSIEkStm3bhvfff58XNw3WUAmAhJwgiD6izuJgOePMIhdFES0tLbjrrrvw1Vdf+ZwfrjWu2Wp2TzMMBgOKi4uRnp6OhoYGfPDBB4piHcD/CyM5ORlff/01vF4vEhISkJeX53OMzWZTvAjUVZ6DCQk5QRAR46/0vr6+HlarFZ9//jk++uijsLoRaiEIAizXXu+zfdFn+/H377/nGS8pKSkoKSlBWloampubefm8IAjQ6XT41a9+pWgty5g+fTp/AeTk5GgKM3tBSZLEA5uD5RNXQ0JOEETEaJXey90Mubm5USns0RLx39g/5Ba2WlTlfuuuri5eYPTBBx/4XCchIUHRhra5udknY0V9TZPJ5NPpcDAhIScIImLUpfdy61Sd5REJWm6UN6+YgOPHjytS/m688Ubk5OTwdED2MjGbzaiqquJWtBxBEHDddddh4sSJcDgcmi8kNYPtC/cHCTlBEBGjlbEhb0fLytfDRafT4ZVrrvPZ/jtHCywrHwEA7jrR6XS48MILUVZWxn3c8t7g7GWjDnRKkoTm5mY0NzdDr9creqxE4i4JtU1vf1jxghStzjQh0tbWhltvvRXvvvsuUlNTB/LWBEH0M2qf+ZgxY3DkyJGwruEvpXDq1KmoqKjgIlhTU4OysjJeFcpg/nCv1wuDwYD6+noA0CzJl5OdnY2pU6dGJLThtumNpG95IO0ki5wgiD7DrM329nYeYPR4PGGJ+C/POQdrJyl7hb/81TG8dap3qE1hYSEXP7vdjrq6Os2CIhb4BHqrL61WK/Ly8rB3715ezSmKos8LICkpCcuWLQvvi/+Ev1hBqPv7Cgk5QRB9Qm5tsha04aJlhf/XqETs+EnEWYvYqqoqmEwmrFy50q/LRmsiEBNSSZIgiiIKCwsBAHV1dfB4PHwWqPw7heMGCRQrCGV/XyEhJwhCQbgiJrc2/fUv8UfphDRM/cX5im13tdjxo8yqBnqt7FdeeYWnEgYr7GFT7eUCLRdS1jYgLy/P57tG4gYJVt3Z39WfJOQEMYwIJSAXqojV1NSgoaEBkydP5sLpL+RmNBp9BrH7K+zRsqgBKKpFdTqd5r0EQfAJdrL1awmpPAtF7h6KxA0SLKOlPzNeSMgJYpgQikgH8+Wyhleff/45DxxqFdioGTt2LA4ePAiv1xuwMrOkpATp6emKLoNa+HthCIKAiooKzcn1gYSUBU69Xm9UMlj80V+ZKyTkBDFMCCXgpjWSraqqiouZvBlWqAiCgP379wMI3OSKCXdnZycWL16M5557zq+/PVCyXbhtZNXtad1uN2bMmIEzZ874rfKMhGhkrviDhJwghgmhBNzMZjMqKiq4y4RVOxoMhpBFXN3PRJIkv1a4/Div14u6ujpIkgRBEHDJJZfgxIkTvElVsHsyt0q4FrS6cEkQBDQ1NcHj8fit8gyHvrpsQoGEnCCGCaEE3Ox2OxfvnTt3wuv1cpH7/PPPg97DaDTiiiuu4G6XkToRG64J3CucIQiCwp3S1tbGr6n2r+v1eu6TNxgMWLVqVcQl8xkZGUhISODFRbfddhu2b98eFcGVW+GiKPaby4aEnCCGEcECbjabTdErnImPwWCAy+UKev3vv/8en3zyid/KTHWvcAYr4tHKRpGLuF6vR2FhIc9EiYa/Wf2CA4D3338fQN8FV+3OKioqQkpKCvnICYKIPuznv81mU7gxPB4PEhMTYTQaQ7qO1+vFbedfgN9feplie+lBBxLGXQq9Xh+0iCcQhYWFWLNmDf+sFsNIg4nqF1y0UgXV7iz1tKRoQUJOEMMcf6XuDJfLhdOnT+P06dNBrxUomJnwxRdYvHgx1q9fH/KIN3kqol6v1+wTzohmMDFaqYIDNT2IhJwghhjhWKXqjI1I0RLw/9RJOHDwZ1+4x+PB/v37Qx7xJooipk+fzv3V/ibeM/q7DD5SBqJjIgk5QQwhQrVK5ZkUfZ1iryXiCz7erXCjsKwSADwrJTExETk5OXj99df5eYIg8ECmwWBAUlIS3+fxeAKKc3+XwccyJOQEMYQIxSplYt/d3d0vg48FQcCkSZNw4MABLtqsrF5ePJSTk4OJEyfi+uuv51kugiBg8uTJmDJlCnejsO6FwcQ5FoYgDxYk5AQxhAjFKrXZbLxDYaRoulJEiVvhrACIoXWv//mf/4EgCBBFEQkJCXC73fB6vdi3bx9aW1t5YDAccY7VwQ/9DQk5QcQJofi+mfBZrVbNc00mE9rb2yNeQ6BgZnZ2tqZg+7P45VN7ioqKcPz4cT4wWT7UeLiKcziQkBNEHBBuRkZ9fT26u7tRV1eHO+64A2+88UaffOETR52L/7giXbHttRNtsJxo45+TkpJC6n44ZswYZGRkoKGhgX8f5kZpbm7m24aTj7uvkJATRBwQTkaGvKjH6/UqgomREMgKl5Oeno7p06ejsbEx4PU6OjrQ0NAQcodC+a+JWBp4HEuQkBNEHBBORobJZIo4gCln7aQp+OU5oxTb7nG04NuzZzWPdzgcWLp0Kd555x2fdEZ5UFWSJPT09KCzs9NnIg9zo9jtdj5Eory8nPv0dTodEhISotpwaigQODGTIIiYgPm+i4qKkJ+fr9jHRM9utwPoFdS+8qp5mo+Iv3nFBMzOz/ebz11XVwcAeO2115CdnQ1RFKHT6TBixAgsXbqU9xkBenPE/b2MmBvpiSeeQFlZmSIwK/efEz8TkkX+9NNPY+vWrRAEAXl5efjDH/6AXbt24bHHHoPL5cKsWbPwwAMP9PdaCWLYU19fj56eHtTX18NisQCAwndeUVGB9957L+Lra7lRfudo6W0o1WLH4sWLodfr+dg0OW63G1arFWvWrMHSpUt5DjibwiM/vqCgIKBriLmRdDodf3EwizxU/3l/Tq2PNYIK+e7du/Hhhx9iy5YtcLvdyMnJQUZGBkpLS7Fx40aMGTMGS5YswY4dO5CZmTkQayaIYYlc4CRJwlNPPYVLL71UsW3FihVRzQu/67P9OPuTK8Xr9eL5558PGtBUB2bz8vJgMpkUGS3p6el+z8/IyOATh/R6Pe9sGI6PvD97f8ciQYX8uuuuw8svvwy9Xo9Tp07B4/Ggq6sL48aNw9ixYwEAc+bMQWNjIwk5QfQjzE/OWsv+9a9/hU6n47nYofTt1iLUYCYAPpdT6z4JCQnc+lYHZoGf+6bodLqwhj/46wceyOKO1XL9/iIkH7nBYMAzzzyD2bNnIyMjA6dPn1aUzo4ePRqnTp3yOa+rqwttbW2KPx0dHdFbPUEMYdS+bzb0gcE6Bno8HhQVFeGSSy4J6/qjRNFHxE+5zvoVcQDctaH2k4uiiFWrVsFsNvMXDmt/m5GRwXt+i6IIvV6P9vZ2/r3U2Gw2uN1uSJIEt9ut6Q+X+9ELCgp8rqW1hqFMyFkr999/P+6++27cc889OHr0qOLnFSvDVbNhwwasW7cuOisliGGEP9eAw+HwsYYlScIHH3yAkydPhnz9cKxwOePHj8eTTz4Jm82GlpYWbN26lac5siCrv2pMVqhksVhQW1vL/fxqS1nuhvF6vTCZTD7rCGZxD7dy/aBCfvjwYXR3d+OKK67AyJEjMWPGDDQ2Nioi0E6nE6NHj/Y5d+HChZg3b55iW0dHB4qLi6OwdIIYmtjtdjz11FOKbA2r1Qqr1eo3kHn06NGQrj03+RIUp16q2Fbx2QHs/3uXnzOULFq0SJEi2NTUhO7ubkiSBIvFwsvqtaoxzWYzbDYb/xUBaLs8Ojs7g7ph5OmYoihyC18t5kNdwBlBXSttbW149NFH0d3dje7ubrz77rtYsGABjhw5gmPHjsHj8eCtt97CTTfd5HOu0WhEamqq4k9ycnK/fBGCiHXUrhJ/xxQUFPBSdcamTZuwceNGPv4sEl41T/MR8d/YPwwo4oIgIDs7G5mZmVi7di3S0tL4dzCbzSgoKOC/xll3wkCE4vII5oZhvvGKigoUFRUBAGprazVdLMOFoBZ5ZmYm9u3bh1//+tcQRREzZszA7NmzccEFF2DZsmVwuVzIzMxEdnb2QKyXIOKSULMomMtALuL90Wa2wP4hQgmLGgwGLF26FGazGTU1NcjNzYXX6+VFOXl5eTwlUkuY1QHJUFwe8n4xajcMAMUMzEmTJvFmW0DvL5fh4k6RE5KPfNmyZT4VWBkZGdiyZUu/LIoghhqhZlHIM1OiUZ0ZiS98woQJmDat97z09HTYbDa0trYqBlC4XC7YbDYsW7bMrzCzyUNy4Q+1CZY/NwwA/hw9Hg/27t0LSZKg0+kgiiIsFgs8Hs+wSDmUQyX6BDEA+Cux10qhy8/Ph9Pp1Cx1D5VIg5kAMH36dJSVlcFutyM/Px89PT0+g5ElSeJBSC1hVk8e6u7uDjsF0N8zk7/omIjfeOONuPTSS1FbWztsUg7lkJATxACg5VKQu1tEUcQtt9yCpqYmuN3ukLoI+iMSEZ80aRIOHjwIAHjxxRcxc+ZMWK1WdHd3A/DNHxcEIWAuuM1mU7iHdDpdwBRArRdaKNkvzPpevnw5gNCHUAw1SMgJIsr4K1RRN4Rqb29XuAmCdQ0MRqRWuF6vx7XXXovPPvsMHo/Hby+T6667Dnv27OGukkC9Utrb2/mQCZ1Oh8rKSp9uhux8tSjLXSL+sl/MZjMvPgrWPXE4QEJOEFEkWFBTbYWzUvS++MP/8dzzsCptsmLbGx0nsKntuN9zUlNTeQaMx+PB559/DkEQfHqZWCwWvtaJEydi/vz5Acvk1d+vuLiYpyTK97Mxc2wEnPz7B3OJyF8E/ronDjdIyAkiigQLaqr7pYwfPx5Ab71GJKPXIrHCRVHEsmXLsHLlSp4D3tzczEv9Kyoq+Jrr6+sV2SPBgojq75+SkuLz/f2NmRMEIahLZLj1UAkVEnKCiCLB+oazhlAejwderxeHDh0CgLB94s9deQ0uSEhQbLu7ZQ++c/cEPE+v16OyshIlJSVwOByoqalR9AmXJAmdnZ0KqzclJSVoEY/6+7Nqb3VVZkZGBi/2YTABLygoUFjvalihFBuaIUnSsApoBoKEnCCiiL+gJptw43A4NPPCw3GtROoLZ77qkpISAOA54PKhDQaDASaTyac1bqhDLVg/mNLSUrjdbqxcuVLR9MpsNqOyspKnJer1+qACDvxsiTMRB/yX7w9HSMgJoh/REqBI6UtKIUOeaSJ/6chbxFqtVsV6Ozs7AwYR1cFd+cuqu7sbVqtVcU5JSQnS0tLCCkoyl438GYbbRXEoQ0JOEFFE7cO9+eabeT/vvtBXERcEQTPTRB0ctNvtsFgsijRDeVWmGi2fdSgECkpqZf3IXTbykW/DKcUwECTkBBFF1MHMbdu29el6fRFwFrxcvHgxjEZjSNYvayHLCPYrQiu4m5eXp8h2AeDT0IptUwt2oGpQrV8P5B/vhYScICIkmOUIIKJMFAA4T6/HC1OVIvWD240/tITWFIoV70iSBKPR6JOmp0bux5cHI+UBxWDflwU3zWazIttl06ZNqKurU/jn/VnygapBh2tqYSiQkBNEBMjL1w0GA+rr67nIXHXVVWhubo742pFa4Swv+9e//jUaGhr8NrJSo7aCFy9ejPXr1ysKf/yl/bHgJju/vLycBzeZdc8mGpWVlSn2aU0RCqcalPgZEnKCiAB5+ToL6AFAbm5uxP1RFqSMxfwxKYptqz9vRUvX3xTb9Hq9wmKW9xtZvnw5zGYzFi5cqKierKqq8smiYfvUVrDRaMRrr72msL6rqqr85sd3dnZy3zWrCmWTguTWvdfrVezTyoRJSEjoHfSsqgYlAkNCThBRorq6ekCaXN122204c+YMJk+ejBdffJFbyUzEgZ/dEHJrW6/X834urBw+Pz9f0wpWuzHkA5FFUVRYyv5EWZ1qKA9OylvVMobbVJ9oQkJOEBoEGuwLQBHQY5Pi77vvvrDvE26vcJ1Ox4W4ubkZFRUVfgN/Wh0I5f1cvF4vnE4nt4IFQcBtt90W8tpramrQ0NCAnJwcvwIcLNWQ9TKX9xsnwoeEnCBUhFIGzgJ6chfFiRMnwrpPJL7w8ePH48svv+Qujs7OTr+BTHUHQjWSJKGpqQmrVq2Cw+GAxWLB9u3b8f777/u8IFhvcDbwubq6mr8UduzYgbVr1/pdh78gpdpPbrVaeYGSuniJEewFO1whIScIFaEMgVALyvz580POUOlLZeaiRYuwcuVKTReHek1sZJq8QZU8m4aJssPhwPHjx3lgUpIklJWVQZIk/iJTu086OjoUa2toaPAR3UCwDomiKPKMF6fTyatM1cFRdg71WdGGhJwgVATrl1JZWYnnnnuOC93VV18dcpZKIBEXBAEzZ85EVlYWNmzYgAMHDiiO0+l+HrHLRLi1tdWnt7k8nS8/Px9ArysIAE8xLC8v53neFouFi7hOp+MBSvZi0poE1NraipaWFr6enJyckL4/oBRknU4HQRDg9XrxzjvvKPLW5cFRto5QpiwNR0jICUKFv6Cb3W5XuBSAXr9zKCIeihVuMBiQlZUFh8OB1tZWn+O9Xi8aGhq4z9vj8XCr1Z+bQu7Dlwsic520t7ejtraWi/iNN96InJwcLvTq4KU8mAqA+8jDscbla5Vn36gbacmDo3ILnj0rSk38GRJygtBAq3S9oKAg7HL7yecZUf6PkxTbtp7uwAvHjyrcHZIkwe12K4KTarxeLyZPnoydO3f6pPSpf0UA0MzT1rLa5VN1WOZLKL1QSkpKwhJwhnytbIScuoeKPJVS3eO8qKgoaJOt4QYJOUH8hN1u5+lwaqFgVmQ4BHOjXHbZZZg1axZPIWRDFvyh0+lgNBo1U/rU6Xzp6ek+7iEt14S/4cn9WUUpX2tdXZ1P+qNer8ell17KtwXrcU6QkBMEF/BXXnmFi3VtbS1Wr16NkpIS/rOe+Y6D9R/R6hX+tPssrrrhBoxwtPCA3rFjx/Diiy9yN4fJZEJpaalCzJmwsbxvJrZaFnNraysXxoSEBM3URH/53gMtjPJMGIYgCLjyyitx4MAB1NbW8pTEYD3OCRJyYpjjr82sx+PBI488gmPHjuEvf/kLr+IMRiArfGdLC+69917s378ff/3rX+H1etHd3a1oHytHFEWsXr1aIdrAz1Wa8nQ/rZxxdWpirBXcMIFmz9ZgMGDKlCn49NNPfX41VFRU4JFHHuHuJ3k2C0FCTgxztPpcM7xeL6qrq0PqIx5qSuGuXbtQWFiIHTt28Ht0dXX5vEwEQUBRURH3QfvLTJEHMEPpUxJLjafkzbWAnzNr5D579h3ee+89/v3cbjdWr16N119/fRBWHZuQkBPDGmYV+vNNR1PEAcDhcMDhcHA3jU6nw65du3xEPDExkQsbI1D6nTxnPJ76lGi9WLR+Najz1nfv3q3ZGne4QkJODHsi7VYYioCr51O63W5eFs+yMBwOBxdxvV6PwsJCzayMQPntseY26Qta4l5YWKjIWwcoj1wOCTkxrFB3/oukW+G5oh4vXq0UEJfHg4V79yi2TZgwAePHj1fknQPgZfHyPG6g1xIvLCzEmjVrfNaqHq6gJdax5DaJNiUlJTh27BgvxEpMTKQ8chkk5MSwQe1jvuqqq8IWcS0r/D9FCROvSMM9Gf+EZ599lm9ftGgR0tLS0NTUpAiWejweHoisqamBIAh8dBlzpwTq/z1UxToYZWVlmDlz5pD41RFtSMiJYYO6ovDjjz8O+VytXuFbjaPQnXIJSn8S36eeeooX+LDBwGazGatWrVKkFbIeKXa7HeXl5dxXXlFRQeXoQRjOL7JAkJATwwbmY2a54KEW+AQr7GEVkm63m4u4vLycDV5gTJrUW+nJxNrr9UIQBMVE+GD9XghCji74IQQxdMjPz0dycnJIx75qnuYj4v/a0aYIaLIXAhNkVl4uTw00mUyKgOfevXtRUFAAk8kEg8EAURT9Bi8ffPBB6vJHBCUki3zdunV4++23AQCZmZl46KGHsGvXLjz22GNwuVyYNWsWHnjggX5dKEH0BbnPWRCEoMdrWeELPt6teaxOp4Moirz6Uj6pB+i1yNUDjXt6etDZ2Tlsg5dEdAkq5Lt27cL//u//YvPmzRAEAYsWLcJbb72FJ598Ehs3bsSYMWOwZMkS7NixA5mZmQOxZoIIG7XP2R/h9grXqr7UShtkOd7MapeX25NYE30lqJAnJSVhxYoVSPipd8T48eNx9OhRjBs3DmPHjgUAzJkzB42NjSTkRMzCZk4GEnJ/Is66FOr1vf+7uN1uCIKA6dOnY+nSpT6tXdXI0wZNJpPf0WwEESlBhXzChAn870ePHsXbb7+NkpISJCUl8e2jR4/GqVOnfM7t6upCV1eXYpu6Qosg+gLLtVYLJNve1dWFXbt2ITk5GbfccotPTjcQ3Aq//PLLMW3aNMVwhnCFON4sbxqpFl+EnLVy6NAhLFmyBA899BBEUcTRo0f5PtaVTM2GDRuwbt26qCyUINQwv7fcZcG6/q1cuTJoo6srzj0PFWmTFdv+36mT2PDVMcW2w4cP46uvvuLVlvJBE/0hdoMtojRSLf4IScj37NmD+++/H6WlpZg9ezZ2794Np9PJ9zudTowePdrnvIULF2LevHmKbR0dHSguLu7jsgmi1zJmIg70NqDq6elBVVVVUBEPxRd+4YUXorOzE16vFy6XC1artd/nR8aCiFIOe/wRNP3w5MmTuO+++/Dkk09i9uzZAHp7Uxw5cgTHjh2Dx+PBW2+9hZtuusnnXKPRiNTUVMWfUFO/CCIY6rQ+9quwra3N7znrplztI+J3t+zR7JHy0EMPcb+4JEmwWCyw2+0AlGLX09PDJ/D4w263o6qqip/vj3Cv2x+wHHattEgiNglqkb/wwgtwuVy8/wMALFiwAGvWrMGyZcvgcrmQmZmJ7Ozsfl0oQahRp/Vdfvnl+OKLL/weH2pGymWXXYann34aZrMZDocDNTU1fNgxs07DKdgJx8qOhUKgodSAa7gQVMgfffRRPProo5r7tmzZEvUFEUSoZGRkQK/XczfK4cOHNdvOhpNSuHbtWsUcyry8PFgsFt6pUD5VJ1SxC8dVESsiGm/B2eEOlegTcYvZbEZBQQG3mOVuFkY4Iq7X65GWlsY/sxFwWtdl9w9F7MK1sklEiXAhISfiFrvdDqfTySexywlVwFNTU7lPXZIkbi3X1NTgkUceUYi43LUSDrFiZRNDFxJyIi6x2+3Iz8/3yU4ZJYp46ep/Umz7rqcHd6t6hTPa29sB9AZK5V0J1SIO9M1nTVY20Z+QkBMxj1ZetdVq9RHxQFZ4amqqZoEa86nLfevq+ZcAMGbMGDz33HMkxkRMQkJOxDRaGR+tra18qg4A5I1JwW9SxirOq/jsAPb//WfRbmtr0yxaY/3DgZ9dJ2qxB4B/+Zd/IREnYhYSciImYVZ4e3u7IuNj9erV2L17NxffQFb4pEmTcPDgQYXVrdPpcOWVV2LKlClIT0+Hw+GAxWLh/VNMJpNP7nZqaio6Oztp2C8Rs5CQEzGDvG9KeXk5uru7ecMqURQBgA9J1hLwAvuHkCcfHjhwAHq9Hl6v16eEXy7I6enpKCsrg9frRXl5Oe68807s2LGD7+/o6MATTzxB5epEzEJCTgwKar+3ul+4x+PhlrQgCEhJSeGByXBSCtPT05GdnR2w6yCb4MNK/I1GI9auXYuGhgaMHDkS27dvp3J1IqYhIScGHC2/t7xoRqfTKXzXkiShra0t7F7hAFBYWKgo8NFCK8/bbDajpKQEdrsd77//vmIfQcQaJOTEgCC3wLUqHeViKooiLr/8chw4cICfH2xupiiKkCQJoiji6quvhsvlCknEgcB53v2dAz7YnQ6JoQEJOdHvqC3wiooKzUrH/Px8OJ1OvPPOO1zEAw17YEiShMWLF8NoNEYsiIHyvPsrBzwWOh0SQwMScqLfkbeblSTJZ1YlAIV/3O12I+3c8/AnVa/wNzpOYFPbcW6BMz+6TqeD0WjEsmXLBuPrRQy1iyWiRdA2tgQRCHV7Vq12rfJ2s16vF11dXQp3gtVqxdmzZ+HxeOB2u/GqeZqPiP/G/iE2tR0H0OtKWbx4MRITEyGKIhISEuLSd03tYoloQRY5ETFaLpPy8nIfV4G83awgCKiuruZTpVJSUnivk6fTr8KYESMV97i7ZQ++c/cotnm9XhiNxrjvX0I9WIhoQUJORITdbsdTTz3FXSZerxd1dXXcVcBysisqKni72Z6eHkUpPMtGAcLvUhhoAn1NTQ0aGhqQk5MTUrBzMKEeLEQ0ICEnwkY9KxPoFeVPP/1UkTLY0tKC3Nxc3HbbbYq8cDnBBFyehgj0Vms+9thjfsWvpqYGDz/8MADwop5YF3OC6CvkIyfChgXpmKuEwSxxOW63G42NjT5tZoHQrHC5iAuCgDvuuCOgBdvQ0BDwM0EMRcgiJ8JGnfMNwMdtEohAAs5mZLrdbn599hKQJAkmkwmAdmWozWbD5MmTFeX1OTk5kXxFgogrSMiJoKhFUx2kA3rbylosFng8HgiCwP3mwM/ukXNFPV68WmlNq3uFp6amYunSpairq0NiYiJcLhf27t3L0wxZ86pAQdZ7770X+/fvjwsfOUFEAxJyIiD+ilbUQTqz2Yy8vDzYbDa0tLRg69atAHpF/LLLLsPaCy/2ubZWMPPo0aO8gRVz3TARZ2mG6vzrhoYGxWej0ahoc0sQQx3ykRMBkYtmT08Pb/GqlS/Opss3NTVxN8u85Et8RHzVZwcC9khxu92KIKpOp8ONN94Ii8UCoHeqjyiKPP86JyeH8rGJYQ1Z5ERAtBpKsTFrzEqvr69XTO7p6enN+w7kCx8xYgTOnj2reU/WN0Xeenb58uUAfq4AFUURRUVFyMvLg9lsRlpaGuVjE8MWEnJCgZY/vKKigudlm81mrFixgo9Z6+7uhtVq5QOLa2trYbn2ep/rqnuF+xNxVn5fUFCA9PR0RevZqqoqhQslJSWFizblYxPDGRLyIYq/rnqBuu1p+cMB8CEPO3fuBAA4nU7FeZ9//jnsdjvKyspQpxp8DARvNQsorXCPx4OUlBSfQKXWrwOCIEjIhyT+BFmeWaLVbU/d3Ir5w+XVm2VlZbj22msV92tuboZr1WOonap8MYQq4EVFRUhPT1dknjCRVr94qKSdIHwhIY9DgvWwVmd1WK1W1NfXw+VyKXK9rVar4jrq5lYmkwlpaWm8TwrQG4iUBziB8Mrr5ej1elRWVnLLW+3nDpQxQxDEz5CQxxmh9LBmLgjWmMrpdCoKdpgfWm2dy5tbsZxts9mMxYsX49lnn+XXZy+ISAUcALKzs7F06VKfFEb1LwRq80oQwaH0wzjDXzqgmptvvpkX5jQ1NfF0vYSEBJSUlKCgoAAej0dxnYyMDCQkJPi0hv3yyy8V1/7Hc8/1EXHribaQRFwQBNx777144YUXgooytXkliNAgizzOCBbwYxa73I3i8XhQVFSElJQUhduivr4eQK+fur29Ha2trcjPzwcAntZnt9uxbds2fv2+WOEAsHTpUpSVlYV0bH/5xGm8GjHUICGPM4KJG7PY5W4Ug8HAhVl9HRYA3bRpE6+kFEUR6enpfOiD1+vFU5OvROrIcxT3WtRiR9dPPVEAKHzpWrBJPuF+32iKLY1XI4YiJORxSCBxUze0Kigo8BFx+XVsNpuia6EkSXC73SgtLUVaWhqcTmdQK9xkMmHUqFG8t7g/WB/xwYT87sRQhIQ8zgjmFtCy2Fk5vbpTYEZGhiIwKremPR4PXKsew52q62u5UTo7O9HZ2Rl07aF2R+xPKBedGIqEJOTff/89FixYgOeeew6pqanYtWsXHnvsMbhcLsyaNQsPPPBAf6+TQHC3gN1uh9VqBaD0ccvPufPOO7F+/Xp4vV7o9XrccsstuPnmm3H8+HE+uR7ouy9cC7fbPegWMOWiE0ORoEK+d+9ePProozh69CiA3tLq0tJSbNy4EWPGjMGSJUuwY8cOZGZm9vdahz2B3AKs/wkrna+trcXq1avR2dnJz5EkCc899xy3vLu7u9HY2Ki4R38IOEMUxZiwgCkXnRhqBE0/fPXVV1FeXo7Ro0cDAPbt24dx48Zh7Nix0Ov1mDNnjo8YMLq6utDW1qb409HREd1vMMSRdxkMlI7HRJ7h8XhQVlYGk8nEz9HpdH7dG6NE0UfE28+e6ZOIy6cHCYKAwsJCElCC6AeCWuSVlZWKz6dPn0ZSUhL/PHr0aJw6dUrz3A0bNmDdunV9XOLwRWuAgjo9kMFEnlnkQG91psPhwM0334yOjg78wz/8A7Zs2aKYuAP0jxWenZ2NpKQkRdFRXl5en65JEIQ2YQc71XMaWfWgFgsXLsS8efMU2zo6OlBcXBzubYcl8t4nrM+JJEmaomg2m1FfX4/q6mps376du09qamq4YLe0tCj+Xd1+8Rj8buw4xXX+rXU/Wr//e8RrFgQBiYmJvGqTDZsgfzRB9B9hC3lycrKi+53T6eRuFzVGozHsvOHhitYMypaWFkVaoHwSvVbQ0Gw244UXXkBlZSWqq6s1c7ojtcJ1Oh3Gjx+PL774wq97hg2AWL58ObWXJYgBJGwhv+qqq3DkyBEcO3YMqampeOutt5Cbm9sfaxsWsEwTuQuCzaB0uVz8OFaowyxyf0FDu92O9evX+xVbLQFX9wr3x6JFi+BwOLBx40bFdr1ez9fFRJyqJwli4AhbyBMTE7FmzRosW7YMLpcLmZmZyM7O7o+1DVmYyJlMJi7YcuFlMyjl1ZmJiYmoqKhQDFpg1NTU4Mknn0RXVxeMRiP3gavpiy/c6/WitLQUc+fOVWxnQUx1+X80qifpZUAQoRGykDc1NfG/Z2RkYMuWLf2yoKGOXORYUyt1OX1OTg6am5sBBK/OrKmpwcMPP8w/q4c+ANELZno8HmzevFmxTRRFn7VFo3qSSukJInSosnOAkYucTqeDTqdTjDcLdwZlXV1dwPsFE3E2pT5U1K1wKysr/bbRBSKvnqRSeoIIHRLyAUYtcv7cJYGChHa7HdXV1Th8+DC+++47zWM0feF7mnH55Zf3+TsIgoCbbrpJEdSUE43qSSqlJ4jQISEfAKI5rsxutyMvL09R/CPn0pHn4MnJVyq2PX/0S7z79WkAwKFDhxT7wu1/wvz1/kSc0ddsFSqlJ4jQISHvZ6I9rsxqtfoV8f4qr09NTUVWVpbPVPv+hlIXCSI0aEJQPxPqRJ9QsNvtfJCynLvH/dJHxP/wyUdR65GSnp6ONWvWIC0tLSrXIwgiupBF3s9E09fLeofL6c8mV4ykpCTKIiGIGIaEvA9otY1VE01fb0ZGBkRRhMfj6VcB1+v1cP80+Ye1A4g0i4RywQmi/yEhjxB121iLxYL6+nq/Yh4NETObzbjllltw59d/89lXsu9jJCcn97m7ZHZ2NrKyslBXV4fk5GTFpPtwf1mQFU8QAwP5yCNE3Ta2r/7vUNg5N9dHxH9j/xAvXnQ+6uvr/QZBwyErKwvl5eX49NNP8f777/Pt7JfFgw8+GLIgRzM+QBCEf8gijxB129i++r+1XBBs2w+dnbhxz17F8Tu+duK/jh4GAGzfvh1JSUl9FvLMzEzFIApA6UIJ95cF5YITxMBAQh4hrG1sMB95KDAXRHd3NwRBwPTp07ll/HL6VJ/j/1OUsPvYl/yzx+NRtKuNhMzMTNTW1sJut0dNfCkXnCAGBhLyPhAt37e87zgANDY24vuP7D4i/uD+fTh25kfNa4Qj4oIgYObMmcjKyvLJC4+2+FIuOEH0PyTkMUBGRgZ0Oh0X8v7KSJk/fz4mTpzIrWx/Yq0lvvKOjQNZFEQQRHBIyAcBtT/cbDZj8eLFuGZnM4w/uTUYL174C2zdtq3P95w/fz6qqqr4/cPJJpG7frxeL3Q6HRISEigLhSBiBMpaGWCYKD7++OPIzc1FTU0NAODm3R/7iPjD35zCtu3b++T7BoBJkybhxx9/xIoVK/hLJJxsEnY8+8Xg9XopC4UgYgiyyAcQu92Op556ig+S8Hq9+GX9ZuysV/b4jlZhz7nnnoszZ87gwIEDOHDgAIDefPdVq1aFFdBk2SdszTqdjrJQCCKGICGPIoGqGFkBkXzyj9oXftbjwe8++Sjk+wmCgIsuukhzmIQgCPjhhx98rPnu7m40NDT4bZ+rhTwASj5ygog9SMijRDC/c3V1Nc859xfM1OlC93SxdrJ/+ctfsHXrVmzevBkdHR2QJIlfR2v4MgB88MEHaG5uDsvHTdknBBG7kJBHCX+9SFg/lm3btiF1xEg8lX6V4rx1R77AX7/5GoB/4WXodDoUFxf7tJM1m80wGo144okn+P3lk4duueUWAEBHRwf27dun8HGTOBNE/ENCHiW0qhjl/ViikVIoCAJSUlJQUlKi2G6329He3g5RFPn9tVwn6l8N5OMmiKEBCXmUUPuRbTYbWlpasHBMCqYnXaw49veffIQf/Uy69wcbzMxeEMwXD4CLsyiKKCoqGpBOjARBxA4k5GGiDmiqP7e2tqKsrAxerxevXHMdoBLxSDJSEhIS+GBmAAqrmgVQmUslJSWlX0ewEQQRe5CQh4HaNVFRUYHy8nLF59LSUtRd/U8+58oFPDU1FSdOnND0icsrPHU6Ha688koUFhZyd0pVVZVCuIHw28sSBDG0ICEPA6vVynPAAaChoYFXO0qShLq6uqAiDgAXXXQRRo4c6TMIOTs7G0uXLoXVaoXT6URTUxM+/fRTtLa2Ii0tDWaz2ccXn5eXxwc/kLuEIIYnJOQqmKukq6sL+/fvR05ODtLS0mC1WlFXV8dFXBAEnD17llvPr1xznc+1/LlRWlpaNLdnZWVx10dVVRW2b9/ukwXjz89NAk4QwxcSchnMdSK3unfs2AG9Xg+Px6MorpEkCc3NzRih0+FllYi/6zyF548dCfv+nZ2d/O+BenmTn5sgCDkk5DJYLri6GpLNr2SIogiv1xvVLoUjRozwEWvKMCEIIhRIyGWYTCYIggBBEHzEXL5tmukC/PNl4xX7/8/+fTjup1d4IAKlDJLlTRBEKJCQ/4Tdbkd5eTm8Xi9EUcTYsWNx5MjP7pELLrgA33zzTdSscHlvcBJrgiD6wrATcnkw86WXXsKZM2cwcuRIZGZmKlq1Hj16VHHeqpRLcdEvJyi2RSLgqampWLZsmU91JkEQRKQMaSFnfU6cTieSkpKQnp6O8vJyRTATAM6cOYPGxkYA4A2n5Pv7aoWLoogpU6Yo8sEJgiCiRZ+E/M0330R1dTXcbjcWLlyI4uLiaK0rIPKxYw6HAwB41SMLDra2tuKRRx5RFN1o+b7VeL1eCIIAIDoj1wRBQFFREdasWRPWeQRBEKESsZCfOnUKf/7zn/H6668jISEBCxYswPXXX4/LL788mutTwCxsi8UCt9utEGn5FHl/gh3qpB1JknxE/Nvubtyz7+Ow15yYmMhfMgRBEP1BxEK+a9cuTJs2Deeffz4AYObMmWhsbMQf//jHaK1NgVaOtxx1jnekRDOlMDMzE8uXL6dgJkEQ/UrEQn769GkkJSXxz6NHj8a+ffsUx3R1daGrq0uxraOjI6L7+cvxjhZJCYn4ryuvVmz7j88PYl/XdyFf4/rrr8eePXvg9XqRkJBAIk4QxIAQsZDLfclArxUs/wwAGzZswLp16yJfnQxW6egJs/1rKETDCk9ISEBpaSkAUBEPQRADSsRCnpycDLvdzj87nU6MHj1acczChQsxb948xbaOjo6IgqKs0rG6uppnmACAXq9HYWEhfvjhBx4AvfTSS5GUlITzzjsPW7ZsQVtbm+Y1ZyZdjLvG/VKx7bcf74YryKQeNsDB4/FoFvSQgBMEMZBELOQ33HADqqqq8O2332LkyJHYtm0bVq1apTjGaDTCaDT2eZEMs9mMF154ATU1Nairq0NycjKWLl0aUDjLyspQWVmJzZs3Y8SIEfj+++/hcrnwl4mTfI5lVvhll12Gc845B52dnbj44ovhcDjgdrthNBoxd+5cnwwZEm6CIAYTQeqD0/nNN9/E888/j56eHuTl5eHuu+8Oek5bWxtuvfVWvPvuu0hNTY301hHzyT8vx49Hjym2/eqN1wZ8HQRBEOEQSDv7lEc+Z84czJkzp0+LG0h2zs1VfBZHnYNptRsHaTUEQRDRYUhXdjLUAg6QFU4QxNBhSAu5t6cHtrwFim3/+ND/wUW/onFoBEEMHYaskJMVThDEcGHICbnL+TXsi5Yotl338osw/OIXg7QigiCI/mVICTlZ4QRBDEeGhJD/rWUv9pf/SbHthv+x+lSaEgRBDEXiXsjVVvjFM6fj8nvvGaTVEARBDDxxK+Sn3nkXX1Q9q9hGbhSCIIYjcSfkkiRh16+V/b3TV/8Jv5g8eZBWRBAEMbjElZCfOdmBj++5T7GNrHCCIIY7cSXk3+7+CABgMJlg/u9q6AyGQV4RQRDE4BNXQj5m9iwkZ8+AmJg42EshCIKIGeJKyHV6PaCPqyUTBEH0O7rBXgBBEATRN0jICYIg4hwScoIgiDiHhJwgCCLOISEnCIKIc0jICYIg4pwBz+XzeDwAgI6OjoG+NUEQRNzCNJNpqJwBF3Kn0wkAKC4uHuhbEwRBxD1OpxPjxo1TbBMkSZIGchFnz56Fw+FAUlISRFEM69yOjg4UFxdj06ZNSE5O7qcV9g/xunZa98ATr2uP13UD8bF2j8cDp9OJ9PR0jBgxQrFvwC3yESNGwGw29+kaycnJSE1NjdKKBpZ4XTute+CJ17XH67qB2F+72hJnULCTIAgiziEhJwiCiHNIyAmCIOKcuBJyo9GIP/7xjzAajYO9lLCJ17XTugeeeF17vK4biO+1A4OQtUIQBEFEl7iyyAmCIAhfSMgJgiDinJgV8jfffBM5OTmYMWMGNm3a5LP/4MGDmD9/PmbOnImysjK43e5BWKUvwda9bt06ZGVlYe7cuZg7d67mMYPF999/j9tvvx1tbW0++2L1eTMCrT1Wn/m6deswe/ZszJ49G48//rjP/lh+5sHWHqvP/Omnn0ZOTg5mz56Nl156yWd/LD/zgEgxSEdHh5SVlSV1dnZKP/zwgzRnzhzp0KFDimNmz54tffLJJ5IkSdIjjzwibdq0aRBWqiSUdS9ZskT6+OOPB2mF/mlpaZFuv/12afLkydJXX33lsz8Wnzcj2Npj8Znv3LlTKigokFwul9Td3S397ne/k7Zt26Y4JlafeShrj8Vn3tzcLC1YsEDq6emRzpw5I2VlZUmHDx9WHBOrzzwYMWmR79q1C9OmTcP555+Pc845BzNnzkRjYyPf397ejrNnz2Lq1KkAgPnz5yv2DxbB1g0ADocDzz//PObMmYM//elPcLlcg7RaJa+++irKy8sxevRon32x+rwZgdYOxOYzT0pKwooVK5CQkACDwYDx48fjxIkTfH8sP/Ngawdi85lfd911ePnll6HX6/HNN9/A4/HgnHPO4ftj+ZkHIyaF/PTp00hKSuKfR48ejVOnTvndn5SUpNg/WARb9w8//IArrrgCDz74IDZv3oyuri48++yzg7FUHyorK/22TojV580ItPZYfeYTJkzggnH06FG8/fbbyMzM5Ptj+ZkHW3usPnMAMBgMeOaZZzB79mxkZGTg4osv5vti+ZkHIyaF3Ov1QhAE/lmSJMXnYPsHi2DrGjVqFP77v/8b48ePh16vx5133okdO3YMxlLDIlafdyjE+jM/dOgQ7rzzTjz00EO47LLL+PZ4eOb+1h7rz/z++++HzWbDyZMn8eqrr/Lt8fDM/RGTQp6cnMzb3QK9bRvlP5vV+7/++mu/P6sHkmDrPnHiBKxWK/8sSRL0+gHvWxY2sfq8QyGWn/mePXvw+9//Hv/6r/+KefPmKfbF+jMPtPZYfeaHDx/GwYMHAQAjR47EjBkz8Nlnn/H9sf7MAxGTQn7DDTfAZrPh22+/xZkzZ7Bt2zbcdNNNfH9KSgoSExOxZ88eAMAbb7yh2D9YBFv3iBEj8MQTT+Crr76CJEnYtGkTpk+fPogrDo1Yfd6hEKvP/OTJk7jvvvvw5JNPYvbs2T77Y/mZB1t7rD7ztrY2PProo+ju7kZ3dzfeffddXHvttXx/LD/zoAxOjDU4W7ZskWbPni3NmDFDWr9+vSRJkrRo0SJp3759kiRJ0sGDB6Xc3Fxp5syZ0vLlyyWXyzWYy+UEW3djYyPfv2LFiphZNyMrK4tnfsTD85bjb+2x+MxXrVolTZ06Vbrjjjv4n9ra2rh45qGsPRafuSRJ0jPPPCPNmjVLuv3226VnnnlGkqT4++9cCyrRJwiCiHNi0rVCEARBhA4JOUEQRJxDQk4QBBHnkJATBEHEOSTkBEEQcQ4JOUEQRJxDQk4QBBHnkJATBEHEOf8fSZjuEW/zY6sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(x,y,'k.')\n",
    "plt.plot(x,result.slope*x+result.intercept,'r-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Method 2: statsmodels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ordinary least squares fit using [statsmodels](https://www.statsmodels.org/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "smres = smf.ols('NITRAT ~ PHSPHT',df07).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>NITRAT</td>      <th>  R-squared:         </th> <td>   0.972</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.972</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>7.757e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 22 Feb 2022</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>21:22:25</td>     <th>  Log-Likelihood:    </th> <td> -4993.9</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  2210</td>      <th>  AIC:               </th> <td>   9992.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  2208</td>      <th>  BIC:               </th> <td>1.000e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -3.9326</td> <td>    0.103</td> <td>  -38.336</td> <td> 0.000</td> <td>   -4.134</td> <td>   -3.731</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>PHSPHT</th>    <td>   14.7400</td> <td>    0.053</td> <td>  278.514</td> <td> 0.000</td> <td>   14.636</td> <td>   14.844</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>874.728</td> <th>  Durbin-Watson:     </th> <td>   0.269</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>5147.310</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-1.766</td>  <th>  Prob(JB):          </th> <td>    0.00</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 9.589</td>  <th>  Cond. No.          </th> <td>    4.90</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                 NITRAT   R-squared:                       0.972\n",
       "Model:                            OLS   Adj. R-squared:                  0.972\n",
       "Method:                 Least Squares   F-statistic:                 7.757e+04\n",
       "Date:                Tue, 22 Feb 2022   Prob (F-statistic):               0.00\n",
       "Time:                        21:22:25   Log-Likelihood:                -4993.9\n",
       "No. Observations:                2210   AIC:                             9992.\n",
       "Df Residuals:                    2208   BIC:                         1.000e+04\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -3.9326      0.103    -38.336      0.000      -4.134      -3.731\n",
       "PHSPHT        14.7400      0.053    278.514      0.000      14.636      14.844\n",
       "==============================================================================\n",
       "Omnibus:                      874.728   Durbin-Watson:                   0.269\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             5147.310\n",
       "Skew:                          -1.766   Prob(JB):                         0.00\n",
       "Kurtosis:                       9.589   Cond. No.                         4.90\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smres.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Method 3: pingouin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>names</th>\n",
       "      <th>coef</th>\n",
       "      <th>se</th>\n",
       "      <th>T</th>\n",
       "      <th>pval</th>\n",
       "      <th>r2</th>\n",
       "      <th>adj_r2</th>\n",
       "      <th>CI[2.5%]</th>\n",
       "      <th>CI[97.5%]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Intercept</td>\n",
       "      <td>-3.932572</td>\n",
       "      <td>0.102582</td>\n",
       "      <td>-38.335853</td>\n",
       "      <td>6.540950e-247</td>\n",
       "      <td>0.972323</td>\n",
       "      <td>0.972311</td>\n",
       "      <td>-4.133740</td>\n",
       "      <td>-3.731405</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>PHSPHT</td>\n",
       "      <td>14.740035</td>\n",
       "      <td>0.052924</td>\n",
       "      <td>278.514375</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.972323</td>\n",
       "      <td>0.972311</td>\n",
       "      <td>14.636249</td>\n",
       "      <td>14.843820</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       names       coef        se           T           pval        r2  \\\n",
       "0  Intercept  -3.932572  0.102582  -38.335853  6.540950e-247  0.972323   \n",
       "1     PHSPHT  14.740035  0.052924  278.514375   0.000000e+00  0.972323   \n",
       "\n",
       "     adj_r2   CI[2.5%]  CI[97.5%]  \n",
       "0  0.972311  -4.133740  -3.731405  \n",
       "1  0.972311  14.636249  14.843820  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pg.linear_regression(x[ii],y[ii])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}