linear_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "\n",
    "\n",
    "<a href=\"https://colab.research.google.com/github/saskeli/x/blob/master/linear_regression.ipynb\"><img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open and Execute in Google Colaboratory\"></a>\n",
    "\n",
    "|                                                -                                                |                                                -                                                |                                                -                                                |\n",
    "|-------------------------------------------------------------------------------------------------|-------------------------------------------------------------------------------------------------|-------------------------------------------------------------------------------------------------|\n",
    "|            [Exercise 10 (linear regression)](<#Exercise-10-(linear-regression&#41;>)            |                 [Exercise 11 (mystery data)](<#Exercise-11-(mystery-data&#41;>)                 | [Exercise 12 (coefficient of determination)](<#Exercise-12-(coefficient-of-determination&#41;>) |\n",
    "|    [Exercise 13 (cycling weather continues)](<#Exercise-13-(cycling-weather-continues&#41;>)    |                                                                                                 |                                                                                                 |\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Machine learning: linear regression\n",
    "\n",
    "## Linear regression\n",
    "Regression analysis tries to explain relationships between variables. One of these variables, called dependend variable, is what we want to \"explain\" using one or more *explanatory variables*. In linear regression we assume that the dependent variable can be, approximately, expressed as a linear combination of the explanatory variables. As a simple example, we might have dependent variable height and an explanatory variable age. The age of a person can quite well explain the height of a person, and this relationship is approximately linear for kids (ages between 1 and 16). Another way of thinking about regression is fitting a curve to the observed data points. If we have only one explanatory variable, then this is easy to visualize, as we shall see below.\n",
    "\n",
    "We can apply the linear regression easily with the [scikit-learn](https://scikit-learn.org/stable/) package. Let's go through some examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we make the usual standard imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-24T19:31:02.144483Z",
     "start_time": "2020-06-24T19:31:01.933556Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sklearn   # This imports the scikit-learn library"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we create some data with approximately the relationship $y=2x+1$, with normally distributed errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-24T19:31:02.150290Z",
     "start_time": "2020-06-24T19:31:02.145821Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  0.           0.52631579   1.05263158   1.57894737   2.10526316\n",
      "   2.63157895   3.15789474   3.68421053   4.21052632   4.73684211\n",
      "   5.26315789   5.78947368   6.31578947   6.84210526   7.36842105\n",
      "   7.89473684   8.42105263   8.94736842   9.47368421  10.        ]\n",
      "[  2.76405235   2.45278879   4.08400114   6.39878794   7.07808431\n",
      "   5.28588001   8.26587789   8.21706384   9.31783378  10.88428271\n",
      "  11.67035936  14.03322088  14.39261667  14.80588554  16.18070534\n",
      "  17.12314801  19.33618434  18.68957858  20.26043612  20.14590426]\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(0)\n",
    "n=20   # Number of data points\n",
    "x=np.linspace(0, 10, n)\n",
    "y=x*2 + 1 + 1*np.random.randn(n) # Standard deviation 1\n",
    "print(x)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we import the `LinearRegression` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-24T19:31:02.253488Z",
     "start_time": "2020-06-24T19:31:02.152025Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can fit a line through the data points (x, y):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-24T19:31:02.350861Z",
     "start_time": "2020-06-24T19:31:02.254992Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xl4lNX9/vH3h0WNVomigkQRUETCLlFBCypYwKWSovWLtcqvBAKKVnahrWxVQcBdBMIiYikqS5AKGAkoWFmzsIQooohIQEAwgBAhy/n9kQkNmJCQyWSSmft1Xbkyc+ZZPnOJ95ycOc95zDmHiIgEj0r+LkBERMqWgl9EJMgo+EVEgoyCX0QkyCj4RUSCjIJfRCTIKPhFRIKMgl9EJMgo+EVEgkwVfxdQkEsvvdTVqVPH32WIiFQYiYmJPzrnLivOtuUy+OvUqUNCQoK/yxARqTDM7LvibquhHhGRIKPgFxEJMgp+EZEgo+AXEQkyCn4RkSCj4BcRCTIKfhGRIKPgFxEpB9IefJCE3/62TM6l4BcR8aMjR47wxBNPsG3OHLISEzl69KjPz1kur9wVEQkGH330Eb169eL777+nb1gYdevWpfIFF/j8vOrxi4iUsQMHDtCtWzfuuusuLrjgAj7//HOuvfZaKleuXCbnV/CLiJQR5xxz584lPDycf//73zzzzDMkJyfTunXrMq1DQz0iImVgz5499OnTh9jYWFq2bMnHH39Ms2bN/FKLgl9ExEsLktMYF7eV3ekZ1AoNYVDHBkS2CANye/lvvfUW/fv35/jx44wdO5Z+/fpRpYr/4lfBLyLihQXJaQydv5mMzGwA0tIzGDp/MwDNQk8QHR1NfHw8bdu2ZerUqdSvX7/AY9Temc6JrGwGjFl+ygeHLxQ5xm9mV5nZJ2aWamZbzOwpT/slZrbUzLZ5fl9cyP7dPNtsM7Nupf0GRET8aVzcVjIysxkWH8Ow+BgAMjKzeWbOeho3bszatWt58803+eSTTwoN/aHzN3Mi69QPjgXJaT6ruThf7mYBA5xz4UAroI+ZhQNDgGXOufrAMs/zU5jZJcBw4GbgJmB4YR8QIiIV0e70DADC920nfN/2k+2Hsytz2223sWXLFh577DEqVSo4bvM+OPLLyMxmXNxWn9Vc5FCPc24PsMfz+IiZfQGEAZ2B2z2bvQ18Cjx92u4dgaXOuYMAZrYU6ATMLoXaRUT8rlZoCGme8M/v4nMcixYtwszOuH/eB0fq5fUKbPeFsxrjN7M6QAtgLVDD86EA8ANQo4BdwoDv8z3f5WkTEQkIgzo2YPDcDae0nVelEiO6NC8y9OF/Hxyj7oz+VbuvFHsev5n9BpgH9HXOHc7/mnPOAc6bQsws2swSzCxh//793hxKRKRMZGRksGr2q+xZMJ4tF9ci9fJ6hIWGMOb+psX+cnZQxwaEVD31wq2QqpUZ1LGBL0oGitnjN7Oq5Ib+LOfcfE/zXjO7wjm3x8yuAPYVsGsa/xsOAriS3CGhX3HOxQAxABEREV59iIiI+NqKFSvo0aMHX3/9NT169OD+ceMIDQ2l+1keJ+8DorDpoL5QZPBb7t8q04AvnHMv5XtpIdANGOP5/UEBu8cBz+f7QrcDMNSrikVE/Ojw4cM8/fTTTJo0iXr16rFs2TLatWvn1TEjW4T5NOhPV5yhnluBR4B2ZrbB83M3uYH/OzPbBtzpeY6ZRZjZVADPl7r/BNZ7fkblfdErIlLRLFq0iEaNGhETE0P//v3ZtGmT16HvD5Y7PF++REREuISEBH+XISIVxJmunC0NP/74I3379mXWrFk0atSIqVOn0qpVq1I7fmkws0TnXERxttUibSJSoeVdAJWWnoGjdC+Acs7x7rvv0rBhQ95//32GDx9OYmJiuQv9s6XgF5EKrbArZ729ACotLY3IyEgeeugh6tatS2JiIiNGjODcc88tjbL9Smv1iEiFlv/K2YLaz5ZzjqlTpzJw4EAyMzMZP348ffv2LbO18suCevwiUqEVdqFTSS6A+uabb2jfvj3R0dHccMMNbNq0iQEDBgRU6IOCX0QquNK4ACo7O5uktm1Z0qABiYmJxMTEsGzZMq699trSLrdcUPCLSIUW2SKM0V2acE6V3PAPCw1hdJcmxZ7Vk5KSQovInhz/eh8RtRpy/eA5XBZxd6GLqgUCjfGLSIUX2SIMaocC8PmQ4s2rP3HiBKNHj+aluSsJ7dAH27gegL1HMk+up1+WF1WVpcD9SBMRKcS6deto2bIlI0aMoGbHXljVU2fq+HpZZH9T8ItI0Dh27BgDBw6kdevW/PTTT/znP//hRJULgNxlkfMvjezLZZH9TUM9IhIUPvnkE3r06MH27dvp3bs3Y8aMoVq1aoxOWV7myyL7m3r8IhLQDh06RHR0NO3ataNSpUp88sknTJw4kWrVqgH+WRbZ39TjF5GAtXDhQh577DF++OEHBg4cyMiRIzn//PNP2cYfyyL7m4JfRALOvn37+Otf/8p7771HkyZNWLBgATfeeGOh25f1ssj+pqEeEQkYDpg1axbh4eHMnz+fkSNHkpCQcMbQD0bq8YtIQDh+/DhfffUVf/7zn7n55puZNm0ajRo18ndZ5ZKCX0QqtJycHKZMmUJ2YiLOOV5++WWefPLJgFtfpzQp+EWkwtq2bRs9e/ZkxYoVtGvXjilTplCvXr2idwxyRY7xm9l0M9tnZin52t7LdxvGHWa2oZB9d5jZZs92uqWWiJSKrKwsxo8fT9OmTdmwYQNTp04lPj5eoV9MxenxzwDeAGbmNTjn/i/vsZm9CBw6w/53OOd+LGmBIiL5bdq0iaioKBISEujcuTNvvvkmtWrV8ndZFUqRPX7n3EqgwBukm5kBDwKzS7kuEZFTHD9+nGHDhtGyZUt27tzJ+++/T2xsrEK/BLwd428D7HXObSvkdQd8bGYOmOyci/HyfCIShNasWUNUVBSpqan8+c9/5pVXXqF69er+LqvC8nYe/0Ocubf/W+fcDcBdQB8za1vYhmYWbWYJZpawf/9+L8sSkUBw9OhR+vXrxy233MKRI0dYtGgR77zzjkLfSyXu8ZtZFaAL0LKwbZxzaZ7f+8wsFrgJWFnItjFADEBERIQraV0iUvEsSE771ZIJvznwBT179mTHjh08/vjjjB49mosuusjfpQYEb4Z67gS+dM7tKuhFM7sAqOScO+J53AEY5cX5RCQALUhOY+j8zWRkZgOQlp5B/3cT2PuflwmrWpUVK1bQtm2hgwVSAsWZzjkbWA00MLNdZhbleakrpw3zmFktM1vseVoD+K+ZbQTWAYuccx+VXukiEgjGxW0lIzObYfExDIvP/Rowx6pw9X1PsnHjRoW+DxTZ43fOPVRI+/8roG03cLfn8XagmZf1iUiAy7vhSfi+7ae0Z3AeISGBuya+P+nKXRHxG+ccF1XN5lBm5VPufgWBfSMUf1Pwi4hf7Ny5k969e/P1jmNcfs9Tp9wBK9BvhOJvWpZZRMpUTk4Ob775Jo0aNWLFihWM7t2FF7tGEBYaggFhoSGM7tIkqNbHL2vq8YtImfnqq6/o0aMHn332Gb/73e+IiYmhTp06AHRpeZV/iwsi6vGLiM9lZWXxwgsv0LRpUzZv3sxbb71FXFzcydCXsqUev4j41MaNG+nevTtJSUl06dKFCRMmULNmTX+XFdQU/CLitYKuvO3UsDrPPvssL7zwAtWrV2fu3Lncf//9/i5VUPCLiJcKuvJ28NwNDFo7i6+XzaZbt2689NJLXHLJJX6uVPIo+EXEK/mvvAUYdWc0J7KBhp34aFA3Onbs6N8C5VcU/CLilcKuvLULqiv0yykFv4h4peZF57Ln8HFdeVuBKPhFpMRiY2PZ+eHbVG79iK68rUA0j19EztoPP/zAAw88QJcuXbg84zv6t6mpK28rEPX4RaTYnHPMnDmTfv36cezYMZ5//nkGDhxI1apV6evv4qTYFPwiUiw7duygV69efPzxx9x6661MnTqV66+/3t9lSQloqEdEzignJ4fXX3+dxo0bs2rVKt544w1Wrlyp0K/A1OMXkUJ98cUX9OjRg1WrVtGpUycmT55M7dq1/V2WeKk4t16cbmb7zCwlX9sIM0szsw2en7sL2beTmW01s6/NbEhpFi4ivpOZmclzzz1H8+bN+fLLL5k5cyaLFy9W6AeI4gz1zAA6FdD+snOuuedn8ekvmlllYAJwFxAOPGRm4d4UKyK+l5SUxI033sg//vEPOnfuTGpqKo888ghm5u/SpJQUGfzOuZXAwRIc+ybga+fcdufcCeBdoHMJjiMiZSAjI4MhQ4Zw0003sXfvXmJjY3n//fepUaOGv0uTUubNl7tPmNkmz1DQxQW8HgZ8n+/5Lk+biJQzY2bHM7t1JPXjU7m2378ZP2cFkZGR/i5LfKSkwT8RuAZoDuwBXvS2EDOLNrMEM0vYv3+/t4cTkWI4cuQIkU89y5sJhwg/uJtG+77ll8oX8GzcdhYkp/m7PPGREgW/c26vcy7bOZcDTCF3WOd0aUD+e6ld6Wkr7JgxzrkI51zEZZddVpKyROQsLFmyhEaNGpFwIoxKVc875bWMzGzGxW31U2XiayUKfjO7It/TPwApBWy2HqhvZnXN7BygK7CwJOcTkdJz4MABHn30Ue6++24uvPBCqla7HIDUy+udstBa3qqbEniKnMdvZrOB24FLzWwXMBy43cyaAw7YAfTybFsLmOqcu9s5l2VmTwBxQGVgunNui0/ehUiQK+gOWKevleOcY86cOTzxxBP89NNPPPPMM/z973+n3cufk5aeccoia6DVNQNZkcHvnHuogOZphWy7G7g73/PFwK+meopI6SnoDlhD528GOBn+e/bsoU+fPsTGxtKyZUvi4+Np2rQpAIM6Njhlf9DqmoFOSzaIVHD574CVdxesvDF65xzTp0+nYcOGLFmyhLFjx7JmzZqToQ+5Hw6juzTR6ppBREs2iFRwhd0Ba3d6Bh06dCA+Pp62bdsydepU6tevX+AxIluEKeiDiIJfpBwozhh9YWqFhpBWwBexWYf3s3btWiZOnEh0dDSVKukPfMmlfwkifpY3Rp+WnoHjf2P0xZ1HP6hjA0KqVj6lLSfzF2qnJ7NlyxZ69+6t0JdT6F+DiJ+daYy+OCJbhPHP+xpSCYcDco78yP/Vzea/s17hqquuKnJ/CT4a6hHxszON0RdHQkIC/+wRxV++30KD664jccTd6CJIORP1+EX8rLD58kXNoz927BiDBw/m5ptv5scff6TeBx9w19atCn0pkoJfxM8KGqMvah79ihUraNasGePGjSMqKorU1FTuu+8+X5cqAUJDPSJ+ljd755xZlTmRlU3YGWb1HDp0iKeffprJkydTr149li1bRrt27cq6ZKngFPwi5UBkizCoHQrA50MKDvJFixbRq1cv9uzZw4ABAxg1ahTnn39+WZYpAUJDPSLl3P79+3n44Ye59957CQ0NZfXq1YwfP16hLyWm4Bcpp5xzvPvuu4SHhzNnzhyGDx9OUlISN91U0CroIsWnoR6RcigtLY3HHnuM//znP9x4441MmzaNJk2a+LssCRDq8YuUIw6YMmUK4eHhxMfHM378eFavXq3Ql1KlHr9IOZGRkcHWr74iesUKbr/9dqZMmcK1117r77IkACn4RfwsOzubV199lXMSE6lUuTIxMTH06NEDM/N3aRKgFPwifpSSkkJUVBTr1q3j97//PRMnTiQsTMsji29pjF/ED06cOMHIkSO54YYb2L59O7Nnz+aDDz5Q6EuZKM49d6cD9wL7nHONPW3jgN8DJ4BvgL8459IL2HcHcATIBrKccxGlV7pIxbRu3TqioqJISUnhT3/6E6+++iqXXnqpv8uSIFKcHv8MoNNpbUuBxs65psBXwNAz7H+Hc665Ql+C3bFjxxgwYACtW7cmPT2dDz/8kFmzZin0pcwVGfzOuZXAwdPaPnbOZXmergGu9EFtIgHjk08+oUmTJrz00kv07NmTLVu2cM899/i7LAlSpTHG3x1YUshrDvjYzBLNLPpMBzGzaDNLMLOE/fv3l0JZIv536NAhoqOjadeuHZUqVeLTTz9l0qRJXHTRRf4uTYKYV8FvZn8HsoBZhWzyW+fcDcBdQB8za1vYsZxzMc65COdchNYTl0CwcOFCwsPDmTZtGoMGDWLTpk3cdttt/i5LpOTBb2b/j9wvfR92zrmCtnHOpXl+7wNiAS0yIgFv3759dO3alc6dO1O9enXWrl3L2LFjCQk5841VRMpKiYLfzDoBg4H7nHPHCtnmAjO7MO8x0AFIKWmhIuWdc45Zs2YRHh7O/PnzGTVqFAkJCUREaF6DlC/Fmc45G7gduNTMdgHDyZ3Fcy6w1HN14RrnXG8zqwVMdc7dDdQAYj2vVwH+7Zz7yCfvQsTPdu3aRWLbtuz/9lvqt2rFtGnTCA8P93dZIgWyQkZp/CoiIsIlJCT4uwyRIuXk5DBlyhQGDRrEh0ePUq9uXa7YupXKlSsXvbNIKTKzxOJOm9eVuyIltG3bNtq1a0fv3r25vtMjXFC7EbtOhNB23AoWJKf5uzyRQin4Rc5SVlYW48aNo2nTpmzYsIGnXvoXRxrcS2Z27l/PaekZDJ2/WeEv5ZaCX+QsbNq0idatWzN48GA6duxIamoq649fQUZmzinbZWRmMy5uq5+qFDkzBb9IMRw/fpxhw4bRsmVLvvvuO9577z1iY2OpVasWu9MzAEi9vB6pl9c7uU9eu0h5o2WZRYqwZs0aoqKiSE1N5ZFHHuHll1+mevXqJ1+vFRpCWnoGo+489eL0WqGaty/lk3r8IoX4+eef6du3L7fccgtHjhxh8eLFzJw585TQBxjUsQEhVU+dxRNStTKDOjYoy3JFik09fpECxMfH07NnT3bs2EGfPn0YPXo0F154YYHbRrbIXUN/XNxWdqdnUCs0hEEdG5xsFylvFPwi+fz0008MHDiQ6dOnc91117Fy5UratGlT5H6RLcIU9FJhKPhFgAXJaRzt/Ti/ZEPcjZF0fTqct0b04bzzzvN3aSKlTmP8EvRmfLqF/rMTuGbPdhrt206VapezoWojPvrigL9LE/EJBb8ELeccM2fOZNicBHIqnfrHr+bhSyBT8EtQ2rlzJ3fffTfdunWj0m9yZ+loHr4EC43xS1DJyclh0qRJPP300zjneO2115h9NITd6b9oHr4EDfX4JWhs3bqV2267jT59+tC6dWtSUlJ48sknGdzxes3Dl6Ci4JeAl5WVxZgxY2jWrBlbtmxhxowZxMXFUadOHSB3KuboLk0ICw3BgLDQEEZ3aaLpmRKwNNQjAW3Dhg1ERUWRlJREly5dmDBhAjVr1vzVdpqHL8FEPX4JSL/88gt///vfiYiIIC0tjblz5zJv3rwCQ18k2BQr+M1supntM7OUfG2XmNlSM9vm+X1xIft282yzzcy6lVbhIoX5/PPPadGiBc8//zyPPvooqamp3H///f4uS6TcKG6PfwbQ6bS2IcAy51x9YJnn+SnM7BJy79F7M3ATMLywDwgRb/3888/89a9/pU2bNmRkZBAXF8f06dO55JJL/F2aSLlSrOB3zq0EDp7W3Bl42/P4bSCygF07Akudcwedcz8BS/n1B4iI1+Li4nivVi2uef11nnzySVJSUujQoYO/yxIpl7z5creGc26P5/EPQI0CtgkDvs/3fJenTaRUHDx4kP79+/P222+zNiSEBs2bU+3VV/1dlki5ViqzepxzzsycN8cws2ggGqB27dqlUZYEuHnz5tGnTx8OHDjA3/72NyI++4xKlTRfQaQo3vxfstfMrgDw/N5XwDZpwFX5nl/pafsV51yMcy7CORdx2WWXeVGWBLo9e/Zw//3388ADDxAWFsb69et57rnnFPoixeTN/ykLgbxZOt2ADwrYJg7oYGYXe77U7eBpEzlrzjlmzJhBeHg4ixYtYsyYMaxdu5bmzZv7uzSRCqW40zlnA6uBBma2y8yigDHA78xsG3Cn5zlmFmFmUwGccweBfwLrPT+jPG0iZ2XHjh107NiRv/zlLzRu3JiNGzfy9NNPU6WKrkEUOVvF+r/GOfdQIS+1L2DbBKBHvufTgeklqk6CXk5ODhMmTGDo0KGYGRMmTKB3796/GtZZkJxG7Z3pnMjKZsCY5br1ocgZaFBUyq0vv/ySNm3anJybn5KSwuOPP15g6A+dv5kTWdkApKVnMHT+ZhYkF/h1kkjQU/BLuZOZmcnzzz9Ps2bN+PLLL3n77bdZvHgxV199dYHbj4vbSkZm9iltupGKSOE0QCrlSlJSEt27d2fjxo08+OCDvPbaa9SoUdAlIv+Td8OU/DdRyd8uIqdS8Eu5kJGRQWKbNiQmJbG3Rg1iY2OJjCzoYvBfqxUaQlp6hm6kIlJMGuoRv/vss89o3rw5WYmJdKxRg9TU1GKHPsCgjg10IxWRs6DgF785fPgwffr0oW3btmRmZtKsaVOub9CAiy8+u3X8dCMVkbOjoR7xiyVLltCrVy927dpF3759efbZZ7ngnntKfDzdSEWk+BT8UqYOHDhAv379eOedd2jYsCGrVq2iVatW/i5LJKhoqEfKhHOO999/n4YNGzJ79myeeeYZkpOTFfoifqAev/jc7t27efzxx/nggw9o2bIl8fHxNG3a1N9liQQt9filVCxITuPWMcupO2QRt45ZzoLkNJxzTJs2jfDwcOLi4hg7dixr1qwpMPQXJKeRtDOdNdsPnNxfRHxDwS9ey1syIS09A0fukglPz9vIzf/3BD169KB58+Zs2rSJQYMGFbiompZcEClbCn7xWt6SCcPiYxgWHwPA8SzHnktvZOLEiSxfvpz69esXuX9+WnJBxHc0xi9ey1saIXzf9lPaq1x0Gb17Fz1FU0suiJQtBb947Ypq57H70C+/Cu7iLpmgJRdEypaGesQr69ev5+CnM8jJ/IVRd0afDO+zWTJBSy6IlC0Fv5TIsWPHGDx4MK1ateJIynK6hZ9T4iUTtOSCSNky51zJdjRrALyXr6keMMw590q+bW4n916833qa5jvnRhV17IiICJeQkFCiusT3Pv30U3r27MnXX39Nz549GTduHNWqVfN3WSJBzcwSnXMRxdm2xGP8zrmtQHPPCSsDaUBsAZt+5py7t6TnkfLj8OHDDB48mMmTJ3PNNdewfPly7rjjDn+XJSJnqbSGetoD3zjnviul40k5s2jRIsLDw5kyZQr9+/dn06ZNCn2RCqq0gr8rMLuQ11qb2UYzW2JmjUrpfFJG9u/fz8MPP8y9995LaGgoq1at4sUXX+T888/3d2kiUkJeB7+ZnQPcB8wp4OUk4GrnXDPgdWDBGY4TbWYJZpawf/9+b8sSLznnmD17NuHh4cyZM4fhw4eTlJTEzTff7O/SRMRLpdHjvwtIcs7tPf0F59xh59zPnseLgapmdmlBB3HOxTjnIpxzEZdddlkplBVcClorp6TS0tLo3Lkzf/rTn6hXrx5JSUmMGDGCc845pxQrFhF/KY3gf4hChnnMrKaZmefxTZ7zHSiFc0o+Ba2VU5K1bpxzTJkyhfDwcOLj43nxxRdZtWoVjRs39k3hIuIXXl25a2YXAL8DeuVr6w3gnJsEPAA8ZmZZQAbQ1ZV0/qgUKv9aOQCj7ow+udZNcefCf/PNN2y4/XaO7tpFyzvuYMqUKVxzzTW+LFtE/MSr4HfOHQWqn9Y2Kd/jN4A3vDmHFK2wtXKKs9ZNdnY2r776Kv/4xz+IO3GCW667jqeWLcPzh5qIBCCt1RMA8ta6Kaj9TFJSUoiKimLdunX8/ve/56b9+zn33HNBoS8S0LRkQwA427VuTpw4wYgRI7jhhhvYvn07s2fP5oMPPsgNfREJeOrxB4C8cfxzZlXmRFY2YaEhDOrYoMDx/XXr1tG9e3e2bNnCww8/zCuvvMKllxY40UpEApSCv5xYkJzGuLit7E7PoNYZgrswkS3CoHYoAJ8Pafer148ePcqwYcN45ZVXqFWrFh9++CH33FP0WvkiEngU/OVA3nTMvLtQ5U3HBEplhcrly5fTs2dPtm/fTu/evXnhhRe46KKLvD6uiFRMGuMvBwq6dWFp3HowPT2dnj170r59eypXrsynn37KxIkTFfoiQU49/nLAm+mYhVm4cCGPPfYYP/zwA4MGDWLkyJGEhOiOViKiHn+5UNi0y5LcevBEZiZdu3alc+fOVK9enbVr1zJ27FiFvoicpOAvB0rj1oPOOfbu3cu6deuIjY3ln//8J4mJiUREFOu+DCISRDTUUw6czXTMguzcuZPevXvT4csvuaJmTZKXLSM8PNyXJYtIBabgLyeKmo5ZkJycHCZPnszgwYPJycmh4yuv8MATT1C5cuWidxaRoKXgr6C2bdtGjx49WLlyJe3bt2fKlCnUrVvX32WJSAWgMf4KJisri7Fjx9K0aVM2btzItGnTWLp0qUJfRIpNPf4KZOPGjXTv3p2kpCQiIyOZMGECtWrV8ndZIlLBqMdfARw/fpxnnnmGiIgIdu3axZw5c5g/f75CX0RKRMFfzq1evZoWLVrw7LPP8tBDD5GamsoDDzxQ6uvlL0hOI2lnOmu2H/D61o0iUr4p+Mupn3/+mb59+3Lrrbdy9OhRlixZwsyZM6levXrRO5+lvLWCTmSdulaQwl8kMHkd/Ga2w8w2m9kGM0so4HUzs9fM7Gsz22RmN3h7zkC3dOlSmjRpwquvvsrjjz9OSkoKnTp18tn58tYKyq801goSkfKptL7cvcM592Mhr90F1Pf83AxM9PyW02RmZfHNN9/QoUMHrrvuOlauXEmbNm18ft68NYFSL69XYLuIBJaymNXTGZjpucn6GjMLNbMrnHN7yuDcFUZsbCw11q0jMzOTIUOGMHz4cM4777wyOXferRtH3Rn9q3YRCTylMcbvgI/NLNHMogt4PQz4Pt/zXZ42Afbu3csf//hHunTpwrfVqtHwoYcYPXp0mYU+lM5aQSJScZRGj/+3zrk0M7scWGpmXzrnVp7tQTwfGtEAtWvXLoWyyjfnHO+88w59+/bl2LFjPPfcczw4aBBVq1Yt81ry1gTy5g5gIlJxeB38zrk0z+99ZhYL3ATkD/404KpmmgnqAAAJXklEQVR8z6/0tJ1+nBggBiAiIsJ5W1d59t1339GrVy/i4uK45ZZbmDZtGtdff71fa4psEaagFwkSXg31mNkFZnZh3mOgA5By2mYLgUc9s3taAYeCdXw/JyeHCRMm0LhxY/773//y+uuv89lnn/k99EUkuHjb468BxHouJqoC/Ns595GZ9QZwzk0CFgN3A18Dx4C/eHnOCmnr1q1ERUXx+eef06FDByZPnkydOnX8XZaIBCGvgt85tx1oVkD7pHyPHdDHm/NUZJmZmYwfP56RI0dy/vnnM2PGDB599NFSv/JWRKS4tEibDyUnJxMVFUVycjJdunRhwoQJ1KxZ099liUiQ05INPvDLL7/wt7/9jRtvvJHdu3czd+5c5s2bp9AXkXJBPf5SsiA5jXFxW0lLz4BjB9kfv4pHH32UF198kYsvvtjf5YmInKTgLwULktMYMm8Tgz/K/Wpj1J3RhEUO4r4/Nlfoi0i5o6GeUjBqwQZ+ycohfN92wvdtB+BEDlrkTETKJQW/Fw4ePEi3bt048EsOkLvIWf6FzrTImYiURxrqKaF58+bRp08fDhw4wHUD7uMo52mRMxGpENTjP0t79uzh/vvv54EHHiAsLIz169fz3P+10iJnIlJhqMdfTM45ZsyYQf/+/cnIyGDMmDEMGDCAKlWq0NyzjRY5E5GKQMFfDN9++y29evVi6dKltGnThqlTp3Ldddedso0WORORikJDPWeQnZ3Na6+9RuPGjVm9ejVvvvkmn3766a9CX0SkIlGPvxBffPEFUVFRrF69mk6dOjF58uSguE+AiAQ+9fhPk5mZyXPPPUfz5s3ZunUrM2fOZPHixQp9EQkY6vHnk5iYSPfu3dm0aRMPPvggr732GjVq1PB3WSIipUrBD2RkZDBy5EiuGDuWJ0JCuCw2lsjISH+XJSLiEwET/HmLpJ3tdMqVK1fSo0cPtm3bxpc1a3LNNddQRaEvIgEsIMb4FySnMXT+ZtLSM3BAWnoGQ+dvZkHyr27te9Lhw4fp06cPt912G1lZWcTHx9OgQQOqVAmYz0IRkQKVOPjN7Coz+8TMUs1si5k9VcA2t5vZITPb4PkZ5l25BRsXt5WMzGyGxccwLD4GgIzM7EIXSVu8eDGNGzdm4sSJ9O3bl82bN9O+fXtflCYiUu54073NAgY455I8N1xPNLOlzrnU07b7zDl3rxfnKVLeYmh5K2Oe3p7nxx9/pF+/fvzrX/8iPDycVatW0apVK1+WJiJS7pS4x++c2+OcS/I8PgJ8Afjl0tXCFkPLa3fO8d577xEeHs67777LsGHDSEpKOiX0FySnkbQznTXbD3DrmOVnHCYSEanISmWM38zqAC2AtQW83NrMNprZEjNrVBrnO92gjg0KXSRt9+7d/OEPf6Br165cffXVJCYmMnLkSM4999yT2+Z9R3AiKxso3ncEIiIVldffZJrZb4B5QF/n3OHTXk4CrnbO/WxmdwMLgPqFHCcaiAbO+mKpvNk758yqzImsbMJCQxjY4Tp+TFxC+MCBHD9+nLFjx9KvX78Cv7zN+44gv7zvCLT+jogEGq+C38yqkhv6s5xz809/Pf8HgXNusZm9aWaXOud+LGDbGCAGICIiwp1tLZEtwqB2KADvPFiHnj27sXz5cm677TamTp3KtddeW+i+ed8F5L+JSv52EZFAUuLgNzMDpgFfOOdeKmSbmsBe55wzs5vIHVo6UNJzFsU5x660NBo3bkyVKlWYNGkSPXv2pFKlM49o1QoNIS09QzdSEZGg4E2P/1bgEWCzmW3wtP0NqA3gnJsEPAA8ZmZZQAbQ1Tl31r354vjpp5/4NjmZw0eO0P7ee5k4cSJXXnllsfYd1LEBQ+dvPmW4RzdSEZFAVeLgd879F7AitnkDeKOk5zgboaGh7K1Vi7p167Jw4UJy/yApnrxxfN1IRUSCgfmoA+6ViIgIl5CQ4O8yREQqDDNLdM5FFGfbgFiyQUREik/BLyISZBT8IiJBRsEvIhJkFPwiIkFGwS8iEmQU/CIiQUbBLyISZMrlBVxmth/4roS7Xwr8ahG4AKf3HPiC7f2C3vPZuto5d1lxNiyXwe8NM0so7tVrgULvOfAF2/sFvWdf0lCPiEiQUfCLiASZQAz+GH8X4Ad6z4Ev2N4v6D37TMCN8YuIyJkFYo9fRETOIGCC38w6mdlWM/vazIb4ux5fM7OrzOwTM0s1sy1m9pS/ayorZlbZzJLN7EN/11IWzCzUzOaa2Zdm9oWZtfZ3Tb5mZv08/65TzGy2mZ3n75pKm5lNN7N9ZpaSr+0SM1tqZts8vy/2xbkDIvjNrDIwAbgLCAceMrNw/1blc1nAAOdcONAK6BME7znPU8AX/i6iDL0KfOScux5oRoC/dzMLA/4KRDjnGgOVga7+rconZgCdTmsbAixzztUHlnmel7qACH7gJuBr59x259wJ4F2gs59r8inn3B7nXJLn8RFywyDg7xVpZlcC9wBT/V1LWTCzakBbYBqAc+6Ecy7dv1WViSpAiJlVAc4Hdvu5nlLnnFsJHDytuTPwtufx20CkL84dKMEfBnyf7/kugiAE85hZHaAFsNa/lZSJV4DBQI6/CykjdYH9wFue4a2pZnaBv4vyJedcGjAe2AnsAQ455z72b1VlpoZzbo/n8Q9ADV+cJFCCP2iZ2W+AeUBf59xhf9fjS2Z2L7DPOZfo71rKUBXgBmCic64FcBQf/flfXnjGtTuT+6FXC7jAzP7s36rKnsudcumTaZeBEvxpwFX5nl/paQtoZlaV3NCf5Zyb7+96ysCtwH1mtoPc4bx2ZvYv/5bkc7uAXc65vL/m5pL7QRDI7gS+dc7td85lAvOBW/xcU1nZa2ZXAHh+7/PFSQIl+NcD9c2srpmdQ+4XQQv9XJNPmZmRO+77hXPuJX/XUxacc0Odc1c65+qQ+994uXMuoHuCzrkfgO/NrIGnqT2Q6seSysJOoJWZne/5d96eAP9CO5+FQDfP427AB744SRVfHLSsOeeyzOwJII7cGQDTnXNb/FyWr90KPAJsNrMNnra/OecW+7Em8Y0ngVmeTs124C9+rsennHNrzWwukETu7LVkAvAqXjObDdwOXGpmu4DhwBjgfTOLIneF4gd9cm5duSsiElwCZahHRESKScEvIhJkFPwiIkFGwS8iEmQU/CIiQUbBLyISZBT8IiJBRsEvIhJk/j/xbT4xhdYAcAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4cc748c550>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "model=LinearRegression(fit_intercept=True)\n",
    "model.fit(x[:,np.newaxis], y)\n",
    "xfit=np.linspace(0,10,100)\n",
    "yfit=model.predict(xfit[:, np.newaxis])\n",
    "plt.plot(xfit,yfit, color=\"black\")\n",
    "plt.plot(x,y, 'o')\n",
    "# The following will draw as many line segments as there are columns in matrices x and y\n",
    "plt.plot(np.vstack([x,x]), np.vstack([y, model.predict(x[:, np.newaxis])]), color=\"red\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear regression tries to minimize the sum of squared errors $\\sum_i (y[i] - \\hat{y}[i])^2$; this is the sum of the squared lengths of the red line segments in the above plot. The estimated values $\\hat{y}[i]$ are denoted by `yfit[i]` in the above code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-24T19:31:02.355766Z",
     "start_time": "2020-06-24T19:31:02.352189Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters: [ 1.88627741] 2.13794752053\n",
      "Coefficient: 1.88627741448\n",
      "Intercept: 2.13794752053\n"
     ]
    }
   ],
   "source": [
    "print(\"Parameters:\", model.coef_, model.intercept_)\n",
    "print(\"Coefficient:\", model.coef_[0])\n",
    "print(\"Intercept:\", model.intercept_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the coefficient is the slope of the fitted line, and the intercept is the point where the fitted line intersects with the y-axis.\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "Note that in scikit-learn the attributes of the model that store the learned parameters have always an underscore at the end of the name. This applies to all algorithms in sklearn, not only the linear regression. This naming style allows one to easily spot the learned model parameters from other attributes.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameters estimated by the regression algorithm were quite close to the parameters that generated the data: coefficient 2 and intercept 1. Try experimenting with the number of data points and/or the standard deviation, to see if you can improve the estimated parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiple features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous example had only one explanatory variable. Sometimes this is called a *simple linear regression*. The next example illustrates a more complex regression with multiple explanatory variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-24T19:31:02.360803Z",
     "start_time": "2020-06-24T19:31:02.357229Z"
    }
   },
   "outputs": [],
   "source": [
    "sample1=np.array([1,2,3])   # The three explanatory variables have values 1, 2, and 3, respectively\n",
    "sample2=np.array([4,5,6])   # Another example of values of explanatory variables\n",
    "sample3=np.array([7,8,10])   # ...\n",
    "y=np.array([15,39,66]) + np.random.randn(3)   # For values 1,2, and 3 of explanatory variables, the value y=15 was observed, and so on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try to fit a linear model to these points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-24T19:31:02.366671Z",
     "start_time": "2020-06-24T19:31:02.361957Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  5.69493795e+00,   3.36972233e+00,   4.20919214e-03]), 0.0)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model2=LinearRegression(fit_intercept=False)\n",
    "x=np.vstack([sample1,sample2,sample3])\n",
    "model2.fit(x, y)\n",
    "model2.coef_, model2.intercept_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's print the various components involved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-24T19:31:02.373120Z",
     "start_time": "2020-06-24T19:31:02.368636Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x:\n",
      " [[ 1  2  3]\n",
      " [ 4  5  6]\n",
      " [ 7  8 10]]\n",
      "b:\n",
      " [[  5.69493795e+00]\n",
      " [  3.36972233e+00]\n",
      " [  4.20919214e-03]]\n",
      "y:\n",
      " [[ 12.44701018]\n",
      " [ 39.6536186 ]\n",
      " [ 66.8644362 ]]\n",
      "product:\n",
      " [[ 12.44701018]\n",
      " [ 39.6536186 ]\n",
      " [ 66.8644362 ]]\n"
     ]
    }
   ],
   "source": [
    "b=model2.coef_[:, np.newaxis]\n",
    "print(\"x:\\n\", x)\n",
    "print(\"b:\\n\", b)\n",
    "print(\"y:\\n\", y[:, np.newaxis])\n",
    "print(\"product:\\n\", np.matmul(x, b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Polynomial regression\n",
    "It may perhaps come as a surprise that one can fit a polynomial curve to data points\n",
    "using linear regression. The trick is to add new explanatory variables to the model.\n",
    "Below we have a single feature x with associated y values given by third degree polynomial,\n",
    "with some (gaussian) noise added. It is clear from the below plot that we cannot explain the data well with a linear function. We add two new features: $x^2$ and $x^3$. Now the model has three explanatory variables, $x, x^2$ and $x^3$. The linear regression will find the coefficients for these variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-24T19:31:02.501116Z",
     "start_time": "2020-06-24T19:31:02.375078Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficients: [-36.65414588 -20.17228669   0.15359003]\n",
      "Intercept: -167.160466064\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD8CAYAAACPWyg8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzsnXd8VFX6/99n0icJBBIgQEihSwIo3VWwIIJYsGCNiq6Ka/erP2tULMTVtay6tmVdpWUVRXdFDSAIdjHSpJcQkhAgQhLSe+b5/XFvwiQkpE0yM8l5v173lbnnnnvuM5OZ+7nneZ5zjhIRNBqNRqNpDhZnG6DRaDQa90OLh0aj0WiajRYPjUaj0TQbLR4ajUajaTZaPDQajUbTbLR4aDQajabZaPHQaDQaTbPR4qHRaDSaZqPFQ6PRaDTNxtPZBrQVISEhEhkZ6WwzNBqNxq3YsGFDloj0aKxehxWPyMhI1q9f72wzNBqNxq1QSqU1pZ52W2k0Go2m2Wjx0Gg0Gk2z0eKh0Wg0mmbTYWMe9VFRUUFGRgalpaXONsUt8PX1JSwsDC8vL2ebotFoXIxOJR4ZGRkEBgYSGRmJUsrZ5rg0IkJ2djYZGRlERUU52xyNRuNidCq3VWlpKcHBwVo4moBSiuDgYN1L02jciISEBCIjI7FYLERGRpKQkNBm1+pUPQ9AC0cz0J+VRuM+JCQkMHv2bIqLiwFIS0tj9uzZAMTGxjr8ep2q56HRaDQdlbi4uBrhqKa4uJi4uLg2uZ4Wj3YmICAAgEOHDjFz5kwnW6PRaDoK6enpzSpvLQ4RD6VUkFJqqVJql1Jqp1LqdKVUd6XUKqXUXvNvN7OuUkq9oZRKVkptUUqNsmtnlll/r1Jqll35aKXUVvOcN1QH8Kf06dOHpUuXtuk1Kisr27R9jUbjOoSHhzervLU4qufxOrBCRIYCI4GdwKPANyIyCPjG3Ae4ABhkbrOBdwCUUt2BOcB4YBwwp1pwzDq32Z03zUF2O43U1FRiYmIAmD9/PpdffjnTpk1j0KBBPPzwwzX1vv76a04//XRGjRrFlVdeSWFhIQDPPvssY8eOJSYmhtmzZyMiAJx99tncf//9jBkzhtdff73935hGo3EK8fHxWK3WWmVWq5X4+Pg2uV6rA+ZKqa7AJOAmABEpB8qVUjOAs81qC4BvgUeAGcBCMe5268xeS2+z7ioRyTHbXQVMU0p9C3QRkXVm+ULgUmB5a+x+5ovt7DiU35omTmBYny7MuTi6Redu3ryZTZs24ePjw5AhQ7jnnnvw8/Nj7ty5rF69Gn9/f1588UVeffVVnnrqKe6++26eeuopAG644Qa+/PJLLr74YgDKy8v1vF4aTSejOigeFxdHeno64eHhxMfHt0mwHByTbRUFHAU+UEqNBDYA9wG9ROSwWScT6GW+7gscsDs/wyw7WXlGPeUnoJSajdGbabOuWlsxefJkunbtCsCwYcNIS0sjNzeXHTt2cMYZZwCGKJx++ukArF27lr/97W8UFxeTk5NDdHR0jXhcffXVznkTGo3GqVx33XVcd9117ZIp6Qjx8ARGAfeIyK9Kqdc57qICQEREKSUOuNZJEZF5wDyAMWPGnPR6Le0htBU+Pj41rz08PKisrEREmDJlCh9++GGtuqWlpdx5552sX7+efv368fTTT9caj+Hv799udms0Gtdhf95+blt1Gy9MfIGxoWPb9FqOiHlkABki8qu5vxRDTP4w3VGYf4+Yxw8C/ezODzPLTlYeVk95h2fChAn89NNPJCcnA1BUVMSePXtqhCIkJITCwsI2D7xrNBr3ICkziSPFRwi1hrb5tVotHiKSCRxQSg0xiyYDO4BlQHXG1Czgc/P1MuBGM+tqApBnurdWAucrpbqZgfLzgZXmsXyl1AQzy+pGu7Y6ND169GD+/Plce+21jBgxgtNPP51du3YRFBTEbbfdRkxMDFOnTmXs2LZ9wtBoNO5BUmYSvf17ExYY1njlVqKqs3Ra1YhSpwLvAd5ACnAzhjB9DIQDacBVIpJjCsCbGBlTxcDNIrLebOfPwONms/Ei8oFZPgaYD/hhBMrvkUYMHzNmjNQNGu/cuZNTTjml1e+3M6E/M43GPbCJjbOWnMWksEnEn9nyDCul1AYRGdNYPYdMTyIim4H6Lja5nroC3NVAO+8D79dTvh6IaaWZGo1G02HZe2wvuWW5jAsd1y7X0yPMNRqNpgOQlJkEoMVDo9FoNE0nKTOJfoH96B3Qu12up8VDo9Fo3JwqWxUbMje0W68DtHhoNBqN27MrZxcFFQVaPDQajUbTdKrjHW09MNAeLR6aWpM0ajQa9+PXzF/p37U/Paw92u2aWjw6MFVVVc42QaPRtDEVtgo2/rGxXXsdoMWj3SkqKuLCCy9k5MiRxMTEsGTJElasWMHQoUMZNWoU9957LxdddBEATz/9NC+//HLNuTExMaSmpgJw6aWXMnr0aKKjo5k3b15NnYCAAB588EFGjhzJL7/8woYNGzjrrLMYPXo0U6dO5fBhY67KDRs2MHLkSEaOHMlbb73Vfh+ARqNxKNuztlNSWdKu8Q7ohGuY17D8Ucjc6tg2Q4fDBS+ctMqKFSvo06cPX331FQB5eXnExMSwZs0aBg4c2OQZcd9//326d+9OSUkJY8eO5YorriA4OJiioiLGjx/PK6+8QkVFBWeddRaff/45PXr0YMmSJcTFxfH+++9z88038+abbzJp0iQeeuihVr91jUbjHJwR7wDd82h3hg8fzqpVq3jkkUf44Ycf2L9/P1FRUQwaNAilFNdff32T2nnjjTcYOXIkEyZM4MCBA+zduxcwZuS94oorANi9ezfbtm1jypQpnHrqqcydO5eMjAxyc3PJzc1l0qRJgLEeiEajcU+SMpMY3G0w3Xy7NV7ZgXTenkcjPYS2YvDgwWzcuJHExESeeOIJJk8+YQaXGjw9PbHZbDX71bPpfvvtt6xevZpffvkFq9XK2WefXXPM19cXDw8PAESE6Ohofvnll1rt5ubmOvptaTQaJ1BeVc7mI5u5cvCV7X5t3fNoZw4dOoTVauX666/noYce4ueffyY1NZV9+/YB1Fq7IzIyko0bNwKwceNG9u/fDxiurm7dumG1Wtm1axfr1q2r91pDhgzh6NGjNeJRUVHB9u3bCQoKIigoiB9//BGAhISENnu/Go2m7fj96O+UVZW1e7wDOnPPw0ls3bqVhx56CIvFgpeXF++88w5ZWVlceOGFWK1WJk6cSEFBAQBXXHEFCxcuJDo6mvHjxzN48GAApk2bxrvvvsspp5zCkCFDmDBhQr3X8vb2ZunSpdx7773k5eVRWVnJ/fffT3R0NB988AF//vOfUUpx/vnnt9v712g0juO3zN+wKAujQ0e3+7UdMiW7K+KuU7J/++23vPzyy3z55ZfONgVwj89Mo+mszFo+i7KqMj666COHtdnUKdm120qj0WjckJLKErZkbXGKywq028rlOPvsszn77LOdbYZGo3FxNh3ZRKWtst1TdKvRPQ+NRqNxQ37L/A1P5cmoXqOccn0tHhqNRuOGJGUmER0Sjb+Xv1Our8VDo9Fo3IiEhAQiB0fye+bvJC1NclqqvRYPjUajcRMSEhKYPXs22X7ZKA/FoXWHmD17tlMExGHioZTyUEptUkp9ae5HKaV+VUolK6WWKKW8zXIfcz/ZPB5p18ZjZvlupdRUu/JpZlmyUupRR9nsynz77bc1EyTWZfr06XqUuEbTCYmLi6O4uJiA6ABsFTaKk4spLi4mLi6u3W1xZM/jPmCn3f6LwN9FZCBwDLjFLL8FOGaW/92sh1JqGHANEA1MA942BckDeAu4ABgGXGvW7bQkJiYSFBTkbDM0Gk07k56eDkDA8ACKdhUhFVKrvD1xiHgopcKAC4H3zH0FnAssNassAC41X88w9zGPTzbrzwA+EpEyEdkPJAPjzC1ZRFJEpBz4yKzrlixcuJARI0YwcuRIbrjhBm666SaWLl1aczwgIKDmdX5+PhdeeCFDhgzhL3/5S808V5GRkWRlZdXbnkaj6biEh4fjFeKFbx9fCrcW1ipvbxw1zuM14GEg0NwPBnJFpNLczwD6mq/7AgcARKRSKZVn1u8L2E/SZH/OgTrl4+szQik1G5gNjX+YLya9yK6cXY29r2YxtPtQHhn3SIPHt2/fzty5c/n5558JCQkhJyeHBx54oMH6SUlJ7Nixg4iICKZNm8Znn33GzJkzT9qeRqPpuMTHx/PgwgcBKNhqTGNktVqJj49vd1ta3fNQSl0EHBGRDQ6wp1WIyDwRGSMiY3r0aL/lGJvKmjVruPLKKwkJCQGge/fuJ60/btw4+vfvj4eHB9dee23NRIYtbU+j0bg3sbGx/Cn2T9iO2ajIrCAiIoJ58+YRGxvb7rY4oudxBnCJUmo64At0AV4HgpRSnmbvIww4aNY/CPQDMpRSnkBXINuuvBr7cxoqbzEn6yG0J/bTrttsNsrLy2uOGd48GtzXaDSdi4qqCg55HuLq8VfzlO0pp9rS6p6HiDwmImEiEokR8F4jIrHAWqDaxzIL+Nx8vczcxzy+RozZGZcB15jZWFHAICAJ+A0YZGZveZvXWNZau53BueeeyyeffEJ2djYAOTk5REZGsmGD0WlbtmwZFRUVNfWTkpLYv38/NpuNJUuWcOaZZzbankaj6bhsPrqZoooizuh7hrNNadO5rR4BPlJKzQU2Af82y/8NLFJKJQM5GGKAiGxXSn0M7AAqgbtEpApAKXU3sBLwAN4Xke1taHebER0dTVxcHGeddRYeHh6cdtppvPjii8yYMYORI0cybdo0/P2PjxYdO3Ysd999N8nJyZxzzjlcdtlljbY3f/78dn5XGo2mvfjp4E94Kk/Gh9Yb9m1X9JTsmpOiPzONxnWYuWwmgd6BfDDtgza7hp6SXaPRaDoQR4qPsPvYbs7se2bjldsBLR4ajUbjBvx86GcALR7OoqO66doC/VlpNK7DTwd/IsQvhMHdBjvbFKCTiYevry/Z2dn6ptgERITs7Gx8fX2dbYpG0+mpslXx86GfOaPPGS6Tst+pVhIMCwsjIyODo0ePOtsUt8DX15ewsDBnm6HRdHq2ZW8jvzzfZVxW0MnEw8vLi6ioKGebodFoNM3ip4M/YVEWJvSe4GxTauhUbiuNRqNxFxISEoiMjMRisfDml28SKqEE+brObNpaPDQajcbFqF70KS0tDYu/BUtfC9u+2Oa0VQPrQ4uHRqPRuBjViz4BBEQHoCyKYxuPOWXRp4boVDEPjUajcQfsF3cKGB5AZUElJaklpNP+iz41hO55aDQajYtRsx6RMsSjcFshiHMWfWoILR4ajUbjYsTHx2O1WvHt54tXVy8KtxY6bdGnhtBuK41Go3Exqhd3embFMwAE5QYxd95cpyz61BBaPDQajcYFiY2NZUW3FZRVlfHx9o+dbc4JaLeVRqPRuCBZJVn8fvR3zg0/19mm1IsWD41Go3FB1h5YiyBaPDQajUbTdNakryEsIIxBQYOcbUq9aPHQaDQaF6OwvJBfD//K5PDJLjOLbl20eGg0Go2L8eOhH6mwVbisywq0eGg0Go3LsSZ9Dd19uzOyx0hnm9IgrRYPpVQ/pdRapdQOpdR2pdR9Znl3pdQqpdRe8283s1wppd5QSiUrpbYopUbZtTXLrL9XKTXLrny0Umqrec4bylX7cRqNRtNKKqoq+CHjB87udzYeFg9nm9Mgjuh5VAIPisgwYAJwl1JqGPAo8I2IDAK+MfcBLgAGmdts4B0wxAaYA4wHxgFzqgXHrHOb3XnTHGC3RqPRuBxJmUkUVhRybj/XdVmBA8RDRA6LyEbzdQGwE+gLzAAWmNUWAJear2cAC8VgHRCklOoNTAVWiUiOiBwDVgHTzGNdRGSdGOvHLrRrS6PRaDoUa9LX4Ofpx/je451tyklxaMxDKRUJnAb8CvQSkcPmoUygl/m6L3DA7rQMs+xk5Rn1lGs0Gk2HwiY21h5Yy5l9z8TX09fZ5pwUh4mHUioA+BS4X0Ty7Y+ZPQZx1LVOYsNspdR6pdR6vU65RqNxN7ZlbeNoyVGXzrKqxiHioZTywhCOBBH5zCz+w3Q5Yf49YpYfBPrZnR5mlp2sPKye8hMQkXkiMkZExvTo0aN1b0qj0WjamW/Sv8FTeTKx70Rnm9Iojsi2UsC/gZ0i8qrdoWVAdcbULOBzu/IbzayrCUCe6d5aCZyvlOpmBsrPB1aax/KVUhPMa91o15ZGo9F0GNakr2FM6Bi6+nR1timN4ohZdc8AbgC2KqU2m2WPAy8AHyulbgHSgKvMY4nAdCAZKAZuBhCRHKXUc8BvZr1nRSTHfH0nMB/wA5abm0aj0XQYUvJSSM1P5bpTrnO2KU2i1eIhIj8CDY27mFxPfQHuaqCt94H36ylfD8S0wkyNRqNxadakrwHgnH7nONmSpqFHmGs0Go0LsDZ9LdHB0YT6hzrblCahxUOj0WicSEJCAlHDo9iStYWNSzaSkJDgbJOahBYPjUajcRIJCQnMnj2b/B7G6Ib0NenMnj3bLQREi4dGo9E4ibi4OIqLiwkcFUhZZhllh8ooLi4mLi7O2aY1ihYPjUajcRLp6el4BHoQMCyA/A35tcpdHUek6mo0Go2mBYSHh1PQvwDlochbl1er3NXRPQ+NRqNxEvHx8XQ7oxulB0spPVAKgNVqJT4+3smWNY4WD41Go3ES5844F7+BfqhdCqUUERERzJs3j9jYWGeb1ijabaXRaDROYnmqMVnGN29+Q7+F/Rqp7VronodGo9E4icSUREaEjKBfF/cSDtDiodFoNE5hX+4+dh/bzfT+051tSovQ4qHRaDROIHF/IhZlYWrkVGeb0iK0eGgcSkJCApGRkVgsFiIjI91ipKxG096ICIkpiYwLHUeIX4izzWkRWjw0DqN6qoW0tDREhLS0tJqpFrSoaDTH2Za1jYzCDKZHuafLCnS2lcaBVE+1YE9xcTH33XcfJSUlNceqRQVwi5REjcbRJO5PxNvizXkR5znblBajex4ah9HQlArZ2dn1ioo7zN+j0TiaKlsVK1JXMDFsIoHegc42p8Vo8dA4jOZOqeAO8/doNI7mtz9+I6sky61dVqDFQ+NA4uPjsVqttcqsVivBwcH11neH+Xs0GkeTmJKIv5c/k8ImOduUVqHFQ+MwYmNjmTdvHhEREbWmWnj99dfrFRV3mL9Ho3Ek5VXlrE5bzeTwyfh6+jrbnFahA+YahxIbG9tgEDwuLo709HTCw8OJj4/XwXJNp+OHgz9QUFHg9i4rcKOeh1JqmlJqt1IqWSn1qLPt6ew0N/U2NjaW1NRUbDYbqampWjg0nZLElES6+3ZnfO/xzjal1biFeCilPIC3gAuAYcC1SqlhzrWq83Ky8RwajaZ+CssL+S7jO86POB9PSxs5fbL3wdJboDinbdq3wy3EAxgHJItIioiUAx8BM5xsU6elofEcOvVWo2mYxP2JlFWVccmASxzfeHkRfPMsvD0B9qyEw787/hp1cJeYR1/ggN1+BuD+/T43oLyqnN05uzlScoSc0hyyS7Ipn1hOvwv74RHggVQIVUVVVBZWUlJUwoe7PiTIJ4jwLuEM7jYYL4uXs9+CRuMSLN2zlCHdhhATEuO4RkVgx/9gZRzkH4QR18CUZyAw1HHXaAB3EY8moZSaDcwGnQbaUmxiY8+xPaw7tI51h9eRdCiJCipq1el+RnfKc8upLKjEw88Dnz4+ePh74GH14Plfn6+p5+PhQ3RwNCN7jGREjxGM6DGCntae7f2WNBqnsz17OztzdvL4+MdRSjmm0SO7YPlDsP976DUcrvg3RJzumLabgLuIx0HAfsL7MLOsFiIyD5gHMGbMGGkf09wfm9j48eCPLNu3jKTDSRwrOwZAMMFkrc0ib1se5VnlVOZX4l3pzazrZ7FgwYJariur1co7/3yH6VdMJ7c0l+S8ZH4/8jtbsraweOdiKrYbAhTZJZJpUdOYHjWdqK5RTnm/Gk178+meT/HEk8cufYzYPbGtyzgszoHv/ga//Qu8A2D6yzDmz2DxcLzhJ0GJuP49VinlCewBJmOIxm/AdSKyvaFzxowZI+vXr28nC92T8qpyvkr5igXbF7Avbx/BvsGc0fcMJvSewLjQcYyPHk9aWtoJ50VERBAfH9/k1NvyqnJ25uzk9yO/8/3B70k6nIQgDAsexvSo6VwQdYHukWg6LMUVxUz8z0Sy12WT9u7x35PVam3ekrNVFbD+A/j2eSjNg1E3wrlPgX/9g3BbilJqg4iMabSeO4gHgFJqOvAa4AG8LyInHWGmxaNh8sry+GTPJyTsTCCrJIsh3YZwU8xNTI2cWitGYbFYqO/7oZTCZrO1+PpHio+wYv8KEvcnsj17O2ITinYWYdlo4ZmbnuH6669vcdsajavx373/5amfnyJlbgrFybUTTSIiIkhNTW28keTVsOJxyNoNUZNg6l8h1IGxEzuaKh7ukm2FiCSKyGARGdCYcHR2GhqDUV5Vzlub32LK0im8vvF1BgUN4p9T/sknF3/CRf0vOiG43VDcqLXxpJ7WntwYfSMX513MgacPcHTZUXxCfbDeYOXZ1GeJWxRHpa2yVdfQaFyFpXuWUnqw9AThgCbM73Z0NyRcCYuvAFsFXPMh3LiszYSjWYhIh9xGjx4tnZHFixeL1WoVoGazWq3y0qKX5NL/XSox82PkwW8flJ3ZO1vc1uLFix1ia0RERE27ykNJ0JlBMvD5gRIzP0amLZ0mDy58UCIGRIhSSiIiIhx23WoWL14sERFt175Gsztnt8TMj5GB1wys9Tuq3iIiIuo/MT9TZNm9Ik8HiTwfJvLTGyIVZe1iM7BemnCPdfpNvq22zioe9jdkQJSXkl5X9pLo96Pl3I/Ple8OfNes9tryBquUOvEHpZAuo7vIlPlTJGZ+jAx5bYh0O7uboBwrXG0tjJrOSd3fy6wFs+S0hafJvEXzmvZ9Ky0QWfO8yNzeIs90F0l8WKTwaLu+h6aKh9vEPJpLZ4152Mcp/Ab40feWvvj28eXY98fY8c4Ounh3cbKFx4mMjGwwIA+Q5ZdFz8t64j/En5L0EjI/zCSkOKRpPuJWXNsR7Ws6H9UzL1RnISovxdDXhxLjH8PHN35MQkJCw0kmVZWwaSGs/SsUHYFhl8LkpyB4QLu/jw4XMG8unVU8IiMjSUtPo9flvQi5MISKYxUc+uAQwQXBLndTrPtjg+MZKDfccEONCHYZ04XQq0Px7uFN/oZ8fpz7I/269Guo2SbRVskAms5L3QeSrqd3pd/t/Sj+oJh93+6r/ySbDXYug7XxkLUHwk+H8+dCWKP37jajwwXMOyvNnYDw6blPE3lPJD0u7sGxH46RHJeMLcVGfHy8y60j3tAU7rGxsbWC8vnr89n7+F4yP8kkMCaQGZ/P4NX1r1JYXtjia7dVMoCm81I3+N397O6UZZax//v9J1YWMTKo/nU2fDILlAWu+Q/cvNypwtEsmuLbcsetI8Q8muuXP1ZyTK7/6vqaAJ19nMLdfPwN2fvuoncl7oc4iZkfI+d9cp78lPGTQ9t31c9D4/rYxxu9e3tLzPwYCZkecmJQPG2dyPsXiMzpIvL3GJFN/xGpqnSKzfWBDpi7v3jUDX5Xb9VfRvvgXOTISJk0f5KMWjhKVuxf0ey2XJGTBes3H9ksl/z3EomZHyNP/vik5JXlObR9jaa52D+QhF4dKtHvRUtgr8Dj36uDG0UWX2mIxt8Givw6r90yqJpDU8VDxzxcmJP55RctWlQTL/CN9CXy/yJRnoobu9zIo9efuNxJR/Txl1WV8c7md/hg+weE+IUw5/Q5br+0p8a9SUhIIO6pOPz+zw85IDw54kliz42Bb1+A3YngGwRn3Avj/wLe/s42t150wLyl4lFyDCxe4BPgeKOaSWPZSGlpaQSMCCD8rnAqCypJeyWNUO/QegPjHTm7aHvWdp746QmSc5O5ZMAlPDz2Ybr6dHW2WZpOypcpX/LYD4/xzqhHOHP7ctj1Jfh2hdPvhvG3G69dGB0wbynfvQSvxRh/S3Kdakp8fHyDa3+np6cTMDKAiPsiKDtcRspzKZQdLmtwxOrJ2nK1QHpj1LV388rNLLloCbePuJ3ElESuWHYF6zPdu9epcU9EhIWb3yVS+fCnT++C/T/A2Y/BfVvgrIddXjiaRVN8W+64tTjmcWC9SMLVhl/y+TCR1c+KFGa1rC0H0JBfPupPUTLsn8NkwNMDxOJraVIMo7623C1w3Ji927K2yfRPp8uIBSPk7c1vS6ULBSI1HZy0dfLrwgskZn6MfPJaf5E18SLFx5xtVbNBxzxaGfM4vAV+eBl2LAMvP2PK4z/d0y6LrDRGal4qV/3vKvKz8tn33D4q84x5oJo9Syfu585qir1FFUXMXTeXL1O+ZEyvMbww8QV6+fdqZ0s1nQIRSPkWfngFUn/grt592GYN4OvLvsQnwD2/c9pt1Vp6j4CrFsKd6+CUi2Hd2/DaCPjiPmOdYAfTVNdRVkkWf1n9F3x9fLmr1130Dep7whiJ5tCQm6vRCducRFPs9ffy568T/0r8mfFsz97OzC9m8n3G9+1loqYzYLPBzi/hvcmw6FLITibl7If43teTa4b/2W2Fo1k0pXvijpvDU3Wzko2Jyp7tITKnq8hH1xsuLgfQVNdRYXmhXLnsShm7eKxsObLFIdd2txTe5tqbkpsiM5fNlJj5MfJi0otSXlXevgZrOhblJSLrPxB5Y9TxcRq//VukolTm/DRHRi8aLdkl2c62slWgx3m00TiP/EyR1c+I/LWf8eV5f7rI7pUiNluLm2zKDbG8slxuW3mbjFwwstmTG56MjhbzqI/SylKJXxcvMfNjZNbyWXK0uH0nmtN0AIpzRL57yRifMaeLyLsTRbYuFamsEBGRo8VHZdTCUfLMz8842dDWo8WjrQcJluaL/PymyCunGF+mf4wV+e19kbKiZjdV7+yyIEopERGpslXJo98/KjHzY+SzPZ85+p243WC5huxt7H18ue9LGbNojJy75FzZfGSzM0zXuBtZycbMtnN7G7/zRZeL7Pv2hIfFNze9KTHzYyQlN8VJhjoOLR5HtcawAAAgAElEQVTtNcK8okxk84ci75xpfLleiDR6JnkHm9xEYz2Pf2/9t8TMj5F//v7PNnoT7k9TeyS7snfJ1KVT5bSFp8nHuz92krUal8ZmMwQi4WrDRf1MsMins0UOb623enFFsZz54Zly9zd3t7OhbUNTxUNnWzkKEUj72Qis7/rKWIw++nJjJGnY6JOeerLZZYefP5xZy2dxTvg5vHLWKyil2vqduCXNyRrLK8vj4e8f5udDP3PFoCt4fPzjeHt4t5OlGpelogS2fgLr3oUj28EaDGNugbG3nDTL8uPdH/Pcuuf4YOoHjAl1k0kNT4IeYe7M6Uly9kPSPNi4CMoLoM9pMPZWiLnCSPuth/rm+r/kyku48osrERE+ueQTl1qLw9Vo7vQrVbYq3tz8Ju9tfY/hIcN57ZzX6Gnt2R6malyNnP2w/n3YtBhKcqBXDEy4A2JmgpfvSU+1iY0Z/5uBv5c/H174YYd4uNPi4QpzW5Xmw5Yl8Nt7cHQX+HWD0643xox073/SU0WEh79/mFVpq5g/bT6n9jy1nYx2T1o6XmV12moe//FxAr0D+ce5/2BY8LA2tFLjMtiqYO8q47eZvNqYEn3odBg3GyInQhNFYG36Wu5dey9/m/Q3Loi6oI2NPhGbTSiuqKKwtJLCsgoKy4zXMX27EGRtWW+6XcRDKfUScDFQDuwDbhaRXPPYY8AtQBVwr4isNMunAa8DHsB7IvKCWR4FfAQEAxuAG0SkXCnlAywERgPZwNUiktqYbS4hHtWIQOqP8Nu/jNxwscGAc2H0TTDkAvDwOuGUz/Z+xpyf53DfqPu4dfit7W+zm3Ey119jY1925+zm7jV3k1uay/MTn2dKxJS2NlfjLAoyYXMCrJ8PeekQEGr8DkfPgi59mt3crOWzOFx0mMTLE/G0eDb5vMoqG4VllRSUVlJYZm6llRSUVVJk9/q4KNjVtf9bXkl9t/BFt4xj4qAezX4/0H7icT6wRkQqlVIvAojII0qpYcCHwDigD7AaGGyetgeYAmQAvwHXisgOpdTHwGci8pFS6l3gdxF5Ryl1JzBCRP6ilLoGuExErm7MNpcSD3vyD8GGBbBxIRQcAv+ecOp1MOrGmiUnU3JTuPrLqxnZcyTzpszDovRYzqZw0mU+GyGrJIv71t7HlqNbuOe0e7ht+G0dwgWhwVjidd83xu9uzwqQKqN3MfZWGHphvQ9v9SEilFXaam7ivx/ZypO/3cYlYX9hdLcZ9d/g67wuMMWgtKLxmayVggBvT/x9PAnw9STAx5NA82+Aj1Hexbe+414MCQ2kq1/T3teJ121nt5VS6jJgpojEmr0OROSv5rGVwNNm1adFZKpZ/phZ9gJwFAg1hej06nrV54rIL0opTyAT6CGNGO6y4lFNVaXRXd64APasrPlCl556HdelfUJ2aQ5LL15KD2vLnh40x2mqqJRVlTHn5zl8lfIV0Sqadc+uIz2l+UKkcRGOpcKmxcimBFTBIaqsIeQNvpLDA2aS7RNR62n/ZE/4BaWVFJUb+5W247cdv7AFeFj3U5j8KNiOx0Y8LarmZl7rhu/rRYCPB4G+XjXHAnw9CTSFIND3+M0/wNcTq5cHFkv7P8A0VTya3s9qnD8DS8zXfYF1dscyzDKAA3XKx2O4qnJFpLKe+n2rzzGFJc+sn+VA29sfD08YMs3Y8g/D5sWwcSEv//A4e7sE8nbACHpk7YN+IU32v2pOpK47Ky0tjdmzZwOcIAY+Hj789cy/UpBawPe27/G41gPLG5aTnqNpO8orbXY3+Ip6nuANF4/9zb6yJJ/hBd8xqXg1p1VtxSaK72wj+KjqKtaUjqJinSes+wP444Tr+XpZCPT1ItDuST7c33rCU32Arye5lXuZt28nMyJu4Zqp59UIQYCPJz6elk7Ra22056GUWg3Ul6cWJyKfm3XigDHA5SIiSqk3gXUistg8/m9guXneNBG51Sy/AUM8njbrDzTL+wHLRSRGKbXNPCfDPLYPGC8iJ4iHUmo2MBsgPDx8dH0BVFdmbdoa7v32PmZ59+H/pWyBimIjsD7yOhh5NQTp9bWbS0sC6ZGRkRzrcYyw28KozK8k7dU0yg6Xuexkka6EiFBcXmXc1Os8vRs3ePPpvo7vvsbXb7dfXtm4a8eiINDHwjleO7mI7ziz4hd8KeOIV182dZvGrtCLsQX2tXv6P+7uqRaCAB9P/L098PRounv41q9vZe+xvSy/fDlWL2vjJ7gRDut5iMh5jVzoJuAiYLKdK+kg0M+uWphZRgPl2UCQUsrT7H3Y169uK8N0W3U169dn6zxgHhhuq8bemytRVFHE3KR4hnQbwn0XfgiVpcaMvr9/CGvnGlv4n2D4FTDsMvAPdrbJbkFLJn5MT09H0oT92fuJuD+C/k/0J/0f6aTvds3JIh1BZZWNorIq4wm/jjunqKyuAFTaCUDFCQJga8Ivz9vDYtzQ7dw7vbv61rhvAny8CPQ1buoBdm6e6nMCvT0IzNuB767/obYuNeKHvl1hxHVw6nX0DBvLVKWY2gafVdLhJH49/CsPjXmowwlHc2iV28rMnHoYOEtEiu0OLQP+o5R6FSNgPghIAhQwyMysOghcA1xn9lbWAjMxMq5mAZ/btTUL+MU8vqaxeIc78vbmtzlafJRXz34VLw8vI4h3WqyxHUuDLR8bA5i+ehCWPwL9z4HhM42An0+gs813WcLDw+vteYSHN9yLqz6nZH8J+57bR8T/RRDx/yIo+7ysLU1tNtUB3Lr+ecO1U1HHn18niFtHIEoqqpp0zePBWuOmHujjSc9A3zo+e886Pv/aPn5/Hw98PD1a9qaP7oFtS2Hbp5CdDBZPGHgeTHseBl/Q6LiM1iIivLHpDXpae3L10Ebzdjo0rc22SgZ8ON4TWCcifzGPxWHEQSqB+0VkuVk+HXgNI1X3fRGJN8v7YwhHd2ATcL2IlCmlfIFFwGlADnCNiKQ0ZpvLB8zt2J2zm6u/vJrLBl3GnNPnNFxRBP7YBlvNH0/eAfD0hUFTYNilMOh88NUDCe1pSQpv3XMsVgtR90XhN8SPMy1nsvzJ5S3K6KrGZhMjAHtCwLb+AK6Ru18nmGs+5VdUNf77bSiA619zY/eoceHYB2+r94+7djydEsAlKxl2LoPtn0HmVkBB1ERj0O0pl4C1e7uZ8n3G99z1zV08OeFJrhpyVbtdtz3RgwTdRDxsYuPG5TeSnp/OF5d90fS1t202yEgyhGTnF1CYCR7eMGAyDJthjB/xC2pb492ElqTw1j3n2bnP8lXVcnawjbykUjL/W4XFw4pflyBuu/NeRo4ZX49//0T3TlFZFYVllSe9djV+Xh4n3sB9au/bH/f3NsSg7tO/2wVwReDIDsNtu3OZ8Rqg7xijtx19mVMWZbOJjau/vJrC8kKWXbYML0vLUmFdHS0ebiIen+75lKd/eZrnzniOSwde2rJGqoVkx+fGln8QLF7G09mQ6TB4GgT1a7ydDoiIUFIzAreBp327/PsT3T/2AdwqvEO+wafHaiqLBlCScUOtFE0wc/PtbvClBbmk79tDcV4OVi/FGeNGM+bU4bWe/usKRKCPF/4+zQvguj22Ksj4DXYvNwQjJwVQEH46DLvEWJCta5hTTVyZupL/993/4/kzn+fiARc71Za2RIuHG4jHsdJjXPy/ixnQdQDzp813zNOhzQaHNsKO/8GuRMgxVz0MHQ5DLjR6JL1Hunz6b5VN6gyyqjjpAKz6fPsF5tN+UwK4Xh7qBN987Rx9T1587hlsZUVYB+YQMjWf8iw4+EE5FUdLOJCyhwAfT6zeHjX/x9aMeu8UlBVA8jfGwL29X0NxthHDiDzTcEcNvQgCXWNFvipbFZctuwwLFj695FM++vCjFg9IdXW0eLiBeDz505N8ue9LPrn4EwZ2G9g2F8naC7sTjSe6A78aU6ME9oaBk41AY/+zjTm3HERZZVUj7psTg7nVA7DsBaC4vGkBXCMbx7N2QLYB901dMah+wvf38cTXq/EArn3ar3+0P+F3h2MrsVH+UTkpv54YhnO39eHbHBHI2mMMjk1ebUzZU1UOvkFG3G7wNOM76YLu1s+TP+eJn57g1bNf5ciPRzr0Q4EWDxcXj41/bGTWilncHHMzD4x+oH0uWpRljGZPXgX71kBpHigPJGws5VHnUND3LHK7nkJBuTTi3qkz/46du6cpAVwPizLTMOsGb2vv18rWqdMTqD7u0Y4B3Lo9Cd9+vkQ+GIl/kD9vnv8mE3pPqFW/uTP9dkhKcmH/d0YPI/kbyM8wyoMHweCpRk+43wRj0KwLkpCQQNxTcfj8xQePSg8e7/c4TzzxRId+KNDi4cLiUWGr4KovrqKwopDPZ3ze4lzxiipb7RG29WXr1HnCr65bXFJGROlORpVv4HTZxEiL8eScJ1bW2Ybxky2an23RJEtfjAxrYwRudf59dbpm9Yhc/0ae8O2P+Xq5WQDXjrqB9EfmPsLXgV+TmpfKs2c8W8sX3il7HuXFcGAd7P8B9n8PhzYZU+/4dIGoSUbPYuBktxjwWv2w4DPeh76z+pL6Siq2fbZaPQ57OspDgRYPFxUPEeG9LR/wxua/89ioF4jpdkaDT/i1B2dV1PHnV1LWhBG4dQO4/nVv7uaNvYfKp3/BesJyf6Nn1jr8iownxCprTyojJuIZdQYekX+CkCFg6USB3CaQX57P/639P5Iyk2pNqtgpYh7lRXBwA6T+BKk/wIEksFUYsYu+o40JCAdOhrCxTZ6A0FWIjIzkwNEDDH5hMGWHy9j/1/0AeHh4UFV1olu1ozwUaPFwsHjUF8Ctnju/vmCuvRjY9w4KynKxDnyJqpIISg7cRPVTfV28PJSdv77OhGrV2TkNTKhm//Tv19LJ1Y6lGk+O1VuhOReQXzcjAyZ8gjHivfdI8NSr8FVUVfDkz0/yVcpXzBw8k7jxcXhaPFs1069LUnjU6Fmkr4P0X+Dw72CrNNbD6D3SEIuos4zvh0+As61tFRaLhdDrQ+l+Tnf2zdlH6YHSmmNWq7XDPhRo8WiheCxal8aKbYdbHMC1envUn63j48WhtH2s++MdAoZncnBxL8oPVeFFFQ/eeyeXXnSBEQcwz2nxCNy2QMRInUz/BdJ+Mf5WZ3F5eEPoCAgbY+Thh42GblEun83VFlSPPn5v63tMCpvES5Necu/pKypKjUF5B9cbvYuDG8wUWsDDx/if9xtvPEz0G+eSge7W0P9P/bHeZiV7dTaZ/8msKY+IiCA+Pr5jPRTYocWjheLxr+9TWLk9s8bFU/2EX9/MmtUjc6t9/41NrhYVE4X1Xiv5G/PJ+GdGTblbdncL/jBE5OB6yFgPhzZDZYlxzBpsLL3be6QhLL1HdCpB+Xj3x8T/Gs/Q7kN5a/JbhPiFONukxqkoMQbjHd4CmVvg4EZjNgObOaAxsA/0HWW4nyKqe5w+zrW5DamyVTFt0TQycjPY8+gebKWGi7gj9TAaQouHC8Y8+lzfh+7ndmfv43sp/6O8prxDBNqqKoybT4b5lHpos7H0rpg9Np+uxliT0OHQ8xRj6zG0w0ynUtc9dfNzN/OV5Su6+3bnnfPeIaprlLNNNBCBvAw4utv4f2VuNbasPbX/V31ONWIWfUcbotGCVfbcibr/vyueu4KvbV9zieUSEp5M6JA9jIbQ4uFi4nGo8BDnLzmfYz8e49D8Q7WOuWXPoylUlJpPs78bT7OHtxj7FXbZKl3CoOdQQ0hCBkHwQGML6OU2PZWGAuPP/PMZvvD6giqp4h/n/oPTep7WfkZVlhtxq5x9xgSCR3bB0Z2GaJQXHq/Xpa/ROwwdbvQQQ4dDUITbfPaOoO7/zyPQg8EvDqa/f3++uOELt80MbClaPFxMPJ786Um+SP6ClLgUCg4W1JR3hm5wLWw2yE2DIzuNm9mRncaNLWsPVNnNWusdYKxlEjwQukUYqZ1B4caNrWu/Np89tTmcLCX3h60/cOfqOzlUeIj4ifFMi5zmmIuKQHGOsQ53bjrkHrATi33GpJli15v172mK9CnQY8jxnl87TiroqtT9//W9tS9dJ3Sl5J0SUtY3Ogdrh0OLhwuJR0peCpd9fhmxp8TSd0/fDhtoaxW2KsOdUn3zy042t31Gua2idv2AXoYrJbCPMUleYG/o0tt4HdDLiLtYg9vFL9/YYMDc0lzuW3sfG49s5P5R9/PnmD83/DQrYvQMirONuFJhpvG34LCR8VaQaXweeQdq9+DAcDcF9zdEt/sACB5w/K8WiQax//9ZB1vp/3h/jn5xlCOfHXF/d3ILcMYytJoGeGvTW/h6+HLr8FvpPra7Fov6sHgYPYxuETDg3NrHbFXGTTO3+ik7zdjyDxtP2+m/QElO/e16BxoLZ1mDwa+7EWPx6WKsgeLbxbjh+gSCl5+xefrW/mvxMMYsKA+712ZShNhAbIwZEsahgxlYFHh7gK+nwtcTIsN6w741BFWWMa/vBTxZVsxrG1/jwO4viPMfild5IZTlGz2I4hxDMEpyjCk7Tvh8PA1RDOgFPQYbg+2C+hm9sKBw47VvUKdyNzmKmjVfPKDPDX0ozyrnyBdHTrrmi0aLR5uzI3sHX6d9ze0jbqe7r376axEWD+ja19giTq+/TkWp8ZSefxiKjho34uIsKMo+/ro4y0g1LcuH0vzabrJWkHQ1QH0LchXAossAY9GbF4Cwbl35F/vIzN7Fy0UQ4B1oCFu3SCMwbQ02egl+3Y/3ogJ7G+V6cGabEB8fz+zZs/Gb5IdvP1/SXk/Dz9OP+Ph4Z5vm0mjxaGP+sekfdPHuwqzoWc42pWPj5WvcgLtFNv2cyjJDRMryjWV/K0qNdGP7v7ZKIwvJVmn0gMRmpq8qowdisYCysC5pPZ9+9l+OZGXj3zWEq667gbPPm2b0YDx9wMsPi08X7vUJJCzta577NZ4be0fx9uS3CfVv/7UpNMeJjY0lW7KZVzKPgs0FdM/pTvw87U5uDB3zaCMSEhJ46p9PYb3FStnXZcyZPkd/GTU1/HLoFx749gH8PP34x+R/EB0c7WyTOi3lVeXEJsaSWZTJp5d8Sk9rT2eb5FSaGvPQ/eA2oDr1T84QKnIrSF6azOzZs0lISHC2aRoX4fQ+p7PwgoV4Wby4aflNrExdWXMsISGByMhILBYLkZGR+nvTxvx9w9/ZlbOL5854rtMLR3PQPY82IDIykiNeRxjwxAAOLT5EzmojmNthx3NoWkx2STb3r72fzUc3c+epdxK4NZDbb7+9w86b5GpUr0l+7dBreXz84842xyXQqbpOFA+LxULYXWH4D/Vn9wO7kXLjM+4QI8k1Dqe8qpxnfnmGZfuWUbG1gj1v7EEqav8u9YOH4zlafJSZX8wk2C+YDy/8EB+PjjvdSnNoV7eVUupBpZQopULMfaWUekMplayU2qKUGmVXd5ZSaq+5zbIrH62U2mqe84YyE+GVUt2VUqvM+quUUo5b9q6NiDg1gi6jupCzJqdGOACd+qepF28Pb+aeMZf7R92PZ7QnUY9F4RlUO5clPT3dSdZ1TGxiI+7HOIorinlp0ktaOFpAq8VDKdUPOB+w/3ZfAAwyt9nAO2bd7sAcYDwwDphjJwbvALfZnVc9FPdR4BsRGQR8Y+67NJPunwRVkPPN8bEHVqtVp/5pGkQpxS3Db6F0SSk+fXwY8NQAfCOPj6LXDx6OZcH2Bfxy+BceHvcwA4IGONsct8QRPY+/Aw8D9v3sGcBCMVgHBCmlegNTgVUikiMix4BVwDTzWBcRWSeGH20hcKldWwvM1wvsyl2SY6XH2GHZwUjvkfQN6otSioiICO2z1tSioaD4Mzc8w+FXDiM2oX9cf4LODGr0wUMH2JvHtqxtvLHxDaZETGHmoJnONsdtadU4D6XUDOCgiPxeZ7qFvsABu/0Ms+xk5Rn1lAP0EpHD5utMoFdrbG5rluxeQmlVKc/OeJYBN+onGs2J1J2ILy0tjdmzZwPUPGA8Ef8E5ReVE3ZrGKep07jqmqta3JbmOPnl+Tzy/SOEWEOYc/qcTjfpoSNptOehlFqtlNpWzzYDeBx4qu3NNDB7JQ1G+JVSs5VS65VS648ePdpeZtVQVlXGh7s+ZFLYJN0V1jRIXFzcCetgFxcXExcXBxg3/f079rP3hb3cHHMzm2QTN628iT+K/mh2W5rjlFeVc9+a+zhUdIgXJ75IV5+uzjbJrWlUPETkPBGJqbsBKUAU8LtSKhUIAzYqpUKBg0A/u2bCzLKTlYfVUw7wh+nWwvx75CS2zhORMSIypkePHo29NYfzxb4vyCnN4abom9r92hr3oaHgd91yT4snD4x+gJfPepm9x/Zy9ZdXs+GPDS1qq7NjExuP//g46/9YT/wZ8YzqNarxkzQnpcUxDxHZKiI9RSRSRCIxXE2jRCQTWAbcaGZdTQDyTNfTSuB8pVQ3M1B+PrDSPJavlJpgZlndCHxuXmoZUJ2VNcuu3KWwiY0F2xcwLHgYY3o1muWm6cQ0FPxuqHxq5FT+M/0/BHgHcOvKW1m4fWHNLLDNbauz8sr6V1iZupIHRj/A9P7TnW1Oh6CtRpgnYvRMkoF/AXcCiEgO8Bzwm7k9a5Zh1nnPPGcfsNwsfwGYopTaC5xn7rsc3x34jtT8VG6Kvkn7UTUnJT4+Hqu19trmjQXFB3YbyGVFl1Gyo4SX1r/EsLhh/Gvxv1rUVmdj0Y5FLNyxkOuGXqe9Ao5ERDrkNnr0aGlPbky8UaZ8MkUqqira9boa92Tx4sUSEREhSimJiIiQxYsXN1rfarUKIN3P6y7D/jVMhr42VOIXxje7rc7Eiv0rZPj84XL/mvulsqrS2ea4BcB6acI9Vs9t1QQaS4XccnQLG49s5IZhN+Bp0RMVaxonNjaW1NRUbDYbqampjWZG2QfGc1bnkPJcClUlVfyn8j/kDMshOSW5yW11Fjb8sYHHf3icU3ueyl8n/hUPi4ezTepQaPFohOpUyLS0NESkJhXSXkAWbF9AoFcglw+63ImWajoydQPgpeml7Ht6H7k/5TJvyzz+vPLPHC483MDZnY9dObu4Z8099AnowxvnvIGvp+ssW9xR0OLRCI2lQh4qPMTq9NVcOeRK/L38nWGiphNQXwDcVmbDc7UnL058kT3H9nD5ssv5dM+n9S6J25mIXxTPzE9mkpOZw7ZntvHVp1/VHNMDKh1IU3xb7rg5KuahlKoeW1JrU0qJiMgbG9+QEQtGyKGCQw65nkZTH/Yxj+rNarXWxDfS89Pl5hU3S8z8GLl15a2SUZDRKWMhjy18TKLfi5aBzw0Uz26etT6nxj5DjQFNjHk4/SbfVpujxCMiIqJe8YiIiJCKqgo5d8m5cseqOxxyLY3mZDQmBlW2Klmya4mMWzxOTpt/moReECqojnujrPt53LPwHomZHyNRj0aJxWo54fd6st+y5jhaPBwkHid7WlmTtkZi5sfIN2nfOORaGo0jOFRwSIbGDTVupI9FiXcv71bfKF2tF1Prd6mQ0GtDJWZ+jPS7s58orxO9Beas3yf1ImgMtHg4MFW3oR/OnavvlHOWnKPTczUuh1JKgs4MklPePkWG/WuY9Lqil1h8LY3eKOv7rruiu6e6F6E8lfS7o5/EzI+R0OtCxcPTo8Hehe55NI2miodeDKqFHC48zLTPpnHr8Fu557R72uw6Gk1LiIyMJC0tDc8gT0KvCiXoT0FU5FVQ9W0Vu5buqjdtte4ki2AMOPTz8yM7O/uE+s5coMpiseDd15uwW8Pwi/Qjc0kmWcuzAMPm+lZiBOp9f3rG69roNczbmP8m/xcR0em5GpekeuR5ZW4lGfMy2PfsPqqyq/Cd4cs1X13Db5m/nXBOQ5mF9QkHnHz+rLbMaqq0VTLw+oEMeHoAXt28SHs9rUY4qpc/iIiIOGE5hNjY2AaPaVpAU7on7ri15QjzyqpKmfzxZLn969vb7BoaTWup64JatGiRJKYkypRPpkjM/Bi595t7ZVf2rpr6DcUEGtqq3Vrt6ebam7NXrv7iaomZHyORd0eKR4CHy7jSOgromEfbicd3B76TmPkxsip1VZtdQ6NpK0oqSuSfv/9Txi0eJzHzY+T2VbdL0uGkBmMCwcHB9YrBHXfcUW95cHCwQ2MLixcvloioCOlxYQ+Jfi9axs4fKyv2r3C5IH5HoanioWMeLeCeNfew9ehWVl25Ci+LV5tcQ6Npa/LK8liyewkJOxPIKc2hN73ZNG8TWb9k1ayaYx8viIuLIz09nfDwcOLj44mLiyMtLa3J11NKYbPZGjyekJBwwjVsYuPBfz1I1yld8e3nS976PI59fIx3X31Xu5vaiKbGPJzeQ2irra16HpmFmTJiwQh5bcNrIuJ6KYwaTVOo9b0dECEPLHxApi2dJjHzY2TYy8MkZFqIRA6PPOn3uSVurpPZY9+LUd5KQi8IlaGvGCnHA+MHSpexXXSGVDtAE3seeha/ZvK/5P9hExuXD7pcLwGqcUtO+N7uS+Pdv7zLu/98lx5n9WDh9oVsDdmKQrEqeBVeu704L+I8uvt2r9VOeHh4vT2P4OBgSkpKTshqOtk08dXBeg9/D7pP7k7wecF4dvGkaG8RBxcdpOD3glpriOrFrpyPdls1gypbFdM/m05ElwjmnT+vJh2yLs5MYdRoGqMp39vUvFSWpy5n+f7l7M/bj4fyYEKfCZwVdhYjQkYwuNtgPv7o4wZTX+FEN1f1A1Ut91REOHc9dxcvffISAdEBWAdYUZ6K/M35ZH2VRfHe4hPsrGurxrE01W2lxaMZ/HjwR+5YfQcvn/UyUyOnYrFYqO/za8y3q9E4k+Z8b0WEPcf2kLg/kZWpKzlYaKwO7W3x5pTgU/DK8uKHj3/g4JaDhHYNZc4jc5gVO+uEBdFsYiOnNIeFny7k+defxxZoI+CUAPyH+eMZYDhASlJLKNxeSO7PuZQdLAMa7sXoFNu2Q4tHG4jH/WvvZ/uTZwsAAAn1SURBVNORTayeuRovDy/d89C4JS393ooIh4sOszVrK1uPbmVr1lZ2ZO+gtKq0Vj0vixddvLvQxacLPh4+5JTkkF2aTZVU1apXcayCwm2FFG4rxPuwN0VZRc3uxWgcT1PFQ8c8msjR4qN8e+Bbboy+ES8PI8MqPj6+3m67XgJU48q09HurlKJPQB/6BPRhauRUACpsFSQfSyatII38snzyy42toLyA/LJ8yqrKOKX7KYT4hRDiF8IdN95BZW4lFbkVVGRV1Gp70aJFDYqEFgsXpClRdXfcHJ1t9f7W9yVmfoyk5qXWKtfZVhp3xFnfWz2/lOuDHufhWLfVlV9cibfFm4QL9eIxGk1LaWj+LB3DcB303FYOZF/uPnbl7GJ6/+nONkWjcWv0/FIdh1aLh1LqHqXULqXUdqXU3+zKH1NKJSuldiulptqVTzPLkpVSj9qVRymlfjXLlyilvM1yH3M/2Twe2Vqbm0vi/kQsylLj59VoNC0nNjaW1NRUbDYbqampWjjclFaJh1LqHGAGMFJEooGXzfJhwDVANDANeFsp5aGU8gDeAi4AhgHXmnUBXgT+LiIDgWPALWb5LcAxs/zvZr12Q0RITElkQu8JhPiFtOelNRqNxmVpbc/jDuAFESkDEJEjZvkM4CMRKROR/UAyMM7ckkUkRUTKgY+AGcpICj8XWGqevwC41K6tBebrpcBkVTeJvA3ZkrWFjMIMpkdpl5VGo9FU01rxGAxMNN1J3ymlxprlfYEDdvUyzLKGyoOBXBGprFNeqy3zeJ5Z/wSUUrOVUuuVUuuPHj3ayrdmkJiSiLfFm8nhkx3Snkaj0XQEGh3noZRaDYTWcyjOPL87MAEYC3yslOrvUAubgYjMA+aBkW3V2vYqbZWsSF3BWf3OIsA7oNX2aTQaTUehUfEQkfMaOqaUugP4zMwNTlJK2YAQ4CDQz65qmFlGA+XZQJBSytPsXdjXr24rQynlCXQ167c5SYeTyCnN4cKoC9vjchqNRuM2tNZt9T/gHACl1GDAG8gClgHXmJlSUcAgIAn4DRhkZlZ5YwTVl5nisxaYabY7C/jcfL3M3Mc8vkbaaXDKV/u/ItArkIlhE9vjchqNRuM2tHZ6kveB95VS24ByYJZ5Y9+ulPoY2AFUAneJGBPbKKXuBlYCHsD7IrLdbOsR4COl1FxgE/Bvs/zfwCKlVDKQgyE4bU5pZSnfpH/D1MipeHt4t8clNRqNxm1olXiYGVPXN3AsHjhhshwRSQQS6ylPwcjGqlteClzZGjtbwncZ31FUUaSzrDQajaYe9AjzBkhMSaSHXw/G9Pr/7d1NiBxlHsfx7y+jBsbX9QUJmkxPFlfISYOIB/WgoplhNStesgRWUZjLCoqIRPriZQ4qehBEGVFQaRVEg6KsuIZl9aJGszEvxmwmMYmGSeILRHEkRvx7qGdCOelOUibzVJn5faCYmqcr9I//01X/6erq1JHvxmhmNtu4eXSxb/8+3tv1HksGl9A3p6/uOGZmjePm0cWqnas48MsBX2VlZtaDm0cXb257k9YZLRads+jIG5uZzUJuHtPs+WEPq3evZnhw+JBbaZqZWcHNY5q3tr9FEAwNDtUdxcyssdw8ppnbN5dr5l9D68xW3VHMzBrLdxI0M7ODfCdBMzObMW4eZmZWmZtHSafTodVqMWfOHFqtFp1Op+5IZmaNdKz/MeIJo9PpMDIywuTkJAA7duxgZGQEwPdYNjObxu88kna7fbBxTJmcnKTdbteUyMysudw8kp07d1YaNzObzdw8kgULFlQaNzObzdw8ktHRUfr7+38z1t/fz+joIbckMTOb9dw8kuXLlzM2NsbAwACSGBgYYGxszB+Wm5l14W+Ym5nZQf6GuZmZzRg3DzMzq8zNw8zMKnPzMDOzytw8zMysshP2aitJXwE7fuc/Pxf4+jjGOV6cqxrnqsa5qmlqLji2bAMRcd6RNjphm8exkPTR0VyqlptzVeNc1ThXNU3NBXmy+bSVmZlV5uZhZmaVuXl0N1Z3gB6cqxrnqsa5qmlqLsiQzZ95mJlZZX7nYWZmlbl5lEh6QNIuSWvTMlx67H5J45I2S7ohc66HJX0maZ2klZLOSuMtST+W8j6ZM1fKsCTVZFzSitzPX8oxX9J/JH0qaaOku9J4zznNmG27pPXp+T9KY2dL+rekLennnzJnurhUk7WSvpN0dx31kvSMpL2SNpTGutZHhcfS622dpMWZc9W+L/bIlf/YFRFe0gI8ANzbZXwR8AkwFxgEtgJ9GXNdD5yU1h8EHkzrLWBDjfXqS7VYCJySarSopizzgMVp/XTg/2neus5p5mzbgXOnjT0ErEjrK6bmtMZ53A0M1FEv4Gpgcfm13Ks+wDDwL0DAFcAHmXPVvi/2yJX92OV3HkdnKfBSROyPiM+BceDyXE8eEW9HxM/p1/eBC3M99xFcDoxHxLaI+Al4iaJW2UXERESsSevfA5uAC+rIcpSWAs+m9WeBv9WY5Vpga0T83i/VHpOIeBf4dtpwr/osBZ6LwvvAWZLm5crVhH2xR716mbFjl5vHoe5Mb0mfKZ1KuAD4orTNl9R3YLqd4i+vKYOS/ifpv5KuypylSXU5SFILuBT4IA11m9OcAnhb0seSRtLY+RExkdZ3A+fXkGvKMuDF0u911wt616dJr7km7YuQ+dg165qHpHckbeiyLAWeAP4MXAJMAI80JNfUNm3gZ6CThiaABRFxKXAP8IKkM3JlbiJJpwGvAHdHxHfUOKclV0bEYmAI+Kekq8sPRnF+oZbLHiWdAtwEvJyGmlCv36izPr00cF/MPm8nzfQTNE1EXHc020l6Cngj/boLmF96+MI0li2XpNuAvwLXpp2JiNgP7E/rH0vaCvwFyHULxRmvSxWSTqZoHJ2IeBUgIvaUHi/PaTYRsSv93CtpJcVpgz2S5kXERDrtsjd3rmQIWDNVpybUK+lVn9pfc03cFw8zbzNWr1n3zuNwpp07vRmYuprhdWCZpLmSBoGLgA8z5loC3AfcFBGTpfHzJPWl9YUp17ZcuYDVwEWSBtNfsMsoapWdJAFPA5si4tHSeK85zZXrVEmnT61TfOC6gaJOt6bNbgVey5mr5O+UTlnVXa+SXvV5HfhHuurqCmBf6fTWjGvqvljLsSvH1QF/lAV4HlgPrEtFn1d6rE1xpcJmYChzrnGK85Zr0/JkGr8F2JjG1gA31lCzYYorm7YC7Rrn7kqKUxvrSnUaPtycZsq1kOJql0/SXLXT+DnAKmAL8A5wdg01OxX4BjizNJa9XhTNawI4QHFO/o5e9aG4yurx9HpbD1yWOVft+2KPXNmPXf6GuZmZVebTVmZmVpmbh5mZVebmYWZmlbl5mJlZZW4eZmZWmZuHmZlV5uZhZmaVuXmYmVllvwJ6KZucu5VHjQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4cc539c320>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x=np.linspace(-50,150,50)\n",
    "y=0.15*x**3 - 20*x**2 + 5*x - 4 + 5000*np.random.randn(50)\n",
    "plt.scatter(x, y, color=\"black\")\n",
    "model_linear=LinearRegression(fit_intercept=True)\n",
    "model_squared=LinearRegression(fit_intercept=True)\n",
    "model_cubic=LinearRegression(fit_intercept=True)\n",
    "x2=x**2\n",
    "x3=x**3\n",
    "model_linear.fit(np.vstack([x]).T, y)\n",
    "model_squared.fit(np.vstack([x,x2]).T, y)\n",
    "model_cubic.fit(np.vstack([x,x2,x3]).T, y)\n",
    "xf=np.linspace(-50,150, 50)\n",
    "yf_linear=model_linear.predict(np.vstack([x]).T)\n",
    "yf_squared=model_squared.predict(np.vstack([x,x2]).T)\n",
    "yf_cubic=model_cubic.predict(np.vstack([x,x2,x3]).T)\n",
    "plt.plot(xf,yf_linear, label=\"linear\")\n",
    "plt.plot(xf,yf_squared, label=\"squared\")\n",
    "plt.plot(xf,yf_cubic, label=\"cubic\")\n",
    "plt.legend()\n",
    "print(\"Coefficients:\", model_cubic.coef_)\n",
    "print(\"Intercept:\", model_cubic.intercept_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Linear and squared are not enough to explain the data, but the linear regression manages to fit quite well a polynomial curve to the data points, when cubic variables are included!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <div class=\"alert alert-info\">Exercise 10 (linear regression)</div>\n",
    "\n",
    "This exercise can give two points at maximum!\n",
    "\n",
    "Part 1.\n",
    "\n",
    "Write a function `fit_line` that gets one dimensional arrays `x` and `y` as parameters. The function should return the tuple `(slope, intercept)` of the fitted line. Write a main program that tests the `fit_line` function with some example arrays. The main function should produce output in the following form:\n",
    "\n",
    "```\n",
    "Slope: 1.0\n",
    "Intercept: 1.16666666667\n",
    "```\n",
    "\n",
    "Part 2.\n",
    "\n",
    "Modify your `main` function to plot the fitted line using matplotlib, in addition to the textual output. Plot also the original data points.\n",
    "\n",
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <div class=\"alert alert-info\">Exercise 11 (mystery data)</div>\n",
    "Read the tab separated file `mystery_data.tsv`. Its first five columns define the features, and the last column is the response. Use scikit-learn's `LinearRegression` to fit this data. Implement function `mystery_data` that reads this file and learns and returns the regression coefficients for the five features. You don't have to fit the intercept. The `main` method should print output in the following form:\n",
    "\n",
    "```\n",
    "Coefficient of X1 is ...\n",
    "Coefficient of X2 is ...\n",
    "Coefficient of X3 is ...\n",
    "Coefficient of X4 is ...\n",
    "Coefficient of X5 is ...\n",
    "```\n",
    "\n",
    "Which features you think are needed to explain the response Y?\n",
    "\n",
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <div class=\"alert alert-info\">Exercise 12 (coefficient of determination)</div>\n",
    "\n",
    "This exercise can give two points at maximum!\n",
    "\n",
    "Using the same data as in the previous exercise, instead of printing the regression coefficients, print the *coefficient of determination*. The coefficient of determination, denoted by R2, tells how well the linear regression fits the data. The maximum value of the coefficient of determination is 1. That means the best possible fit.\n",
    "\n",
    "Part 1.\n",
    "\n",
    "Using all the features (X1 to X5), fit the data using a linear regression (include the intercept). Get the coefficient of determination using the `score` method of the `LinearRegression` class. Write a function `coefficient_of_determination` to do all this. It should return a list containing the R2-score as the only value.\n",
    "\n",
    "Part 2.\n",
    "\n",
    "Extend your function so that it also returns R2-scores related to linear regression with each single feature in turn. The `coefficient_of_determination` (https://en.wikipedia.org/wiki/Coefficient_of_determination) function should therefore return a list with six R2-scores (the first score is for five features, like in Part 1). To achieve this, your function should call both the `fit` method and the `score` method six times. \n",
    "\n",
    "The output from the main method should look like this:\n",
    "```\n",
    "R2-score with feature(s) X: ...\n",
    "R2-score with feature(s) X1: ...\n",
    "R2-score with feature(s) X2: ...\n",
    "R2-score with feature(s) X3: ...\n",
    "R2-score with feature(s) X4: ...\n",
    "R2-score with feature(s) X5: ...\n",
    "```\n",
    "How small can the R2-score be? Experiment both with fitting the intercept and without fitting the intercept.\n",
    "\n",
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <div class=\"alert alert-info\">Exercise 13 (cycling weather continues)</div>\n",
    "\n",
    "Write function `cycling_weather_continues` that tries to explain with linear regression the variable of a cycling measuring station's counts using the weather data from corresponding day. The function should take the name of a (cycling) measuring station as a parameter and return the regression coefficients and the score. In more detail:\n",
    "\n",
    "Read the weather data set from the `src` folder.  Read the cycling data set from folder `src` and restrict it to year 2017. Further, get the sums of cycling counts for each day. Merge the two datasets by the year, month, and day. \n",
    "Note that for the above you need only small additions to the solution of exercise `cycling_weather`.\n",
    "After this, use forward fill to fill the missing values.\n",
    "\n",
    "In the linear regression use as explanatory variables the following columns `'Precipitation amount (mm)'`, `'Snow depth (cm)'`, and `'Air temperature (degC)'`. Explain the variable (measuring station), whose name is given as a parameter to the function `cycling_weather_continues`. Fit also the intercept. The function should return a pair, whose first element is the regression coefficients and the second element is the score. Above, you may need to use the method `reset_index` (its counterpart is the method `set_index`).\n",
    "\n",
    "The output from the `main` function should be in the following form:\n",
    "\n",
    "```\n",
    "Measuring station: x\n",
    "Regression coefficient for variable 'precipitation': x.x\n",
    "Regression coefficient for variable 'snow depth': x.x\n",
    "Regression coefficient for variable 'temperature': x.x\n",
    "Score: x.xx\n",
    "```\n",
    "\n",
    "Use precision of one decimal for regression coefficients, and precision of two decimals for the score.\n",
    "In the `main` function test you solution using some measuring station, for example `Baana`.\n",
    "\n",
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional information\n",
    "\n",
    "* The [scikit-learn](https://scikit-learn.org/stable/) library concentrates on machine learning. Check out library [statsmodels](http://www.statsmodels.org/stable/index.html) for a more statistical viewpoint to regression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary (week 5)\n",
    "\n",
    "* `pd.concat` and `pd.merge` can both combine two DataFrames, but the way the combining is done differs. The function `pd.concat` concatenates based on *indices* of DataFrames, whereas `pd.merge` combines based on the *content* of common variable(s).\n",
    "* The option `join=\"outer` to `pd.concat` can create missing values, but `join=inner` cannot. The former gives the union of indices and the latter gives the intersection of indices.\n",
    "* With `pd.concat` overlapping indices can:\n",
    "\n",
    "     * cause an error\n",
    "     * cause renumbering of indices\n",
    "     * create hierarchical indices\n",
    "* Merging can join elements\n",
    "\n",
    "     * one-to-one\n",
    "     * one-to-many\n",
    "     * many-to-many\n",
    "* In grouping a DataFrame can be thought to be split into smaller DataFrames. The major classes of operations on these groups are:\n",
    "\n",
    "     * aggregate\n",
    "     * filter\n",
    "     * transform (retains shape)\n",
    "     * apply\n",
    "* Series which are indexed by time are called time series\n",
    "* Linear regression can be used to find out linear relationships between variables\n",
    "     * can have more than one feature (explanatory variable)\n",
    "     * fitting polynomials is still linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "\n",
    "\n",
    "<a href=\"https://colab.research.google.com/github/saskeli/x/blob/master/linear_regression.ipynb\"><img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open and Execute in Google Colaboratory\"></a>\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}