1. Simulations/GLM pt. 1 - Linear Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Linear Regression via `sklearn` and `statsmodels`**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **0.** Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model    import LinearRegression\n",
    "from sklearn.metrics         import mean_squared_error, r2_score\n",
    "from numpy.random            import uniform, normal\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **1.** Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case of linear regression $Y$ is a real number and it is modeled as:\n",
    "\n",
    "$\\qquad\n",
    "\\begin{cases}\n",
    "        \\;Y\\sim N(\\mu,\\sigma^2)\\\\\n",
    "        \\;\\mu = \\beta_0 + \\beta_1 \\cdot X\n",
    "\\end{cases}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "\n",
    "# Parameters\n",
    "n_sample = 500\n",
    "beta_0   = 15\n",
    "beta_1   = 2\n",
    "std_dev  = 0.7\n",
    "liminf_x = 0\n",
    "limsup_x = 10\n",
    "\n",
    "# Data simulation\n",
    "x  = uniform(liminf_x, limsup_x, size=n_sample)\n",
    "mu = beta_0 + beta_1*x\n",
    "y  = normal(mu, std_dev)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **2.** Linear Regression with `sklearn`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True model: y ~ N(μ, σ^2)\n",
      "True μ = 15 + 2 * x\n",
      "Regression μ = 14.89 + 2.03 * x\n",
      "Mean squared error: 0.4848\n",
      "Coefficient of determination: 0.9855\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAFNCAYAAADRi2EuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/kElEQVR4nO3de5zWc/7/8ceraWTKYUJSQ6pFqFQaZHPMUgij5BSLddwcFv1SYddxt5LzslYOOR9jB4VE7fJtkakpicKqaEKhQTXVNL1/f3yuq6655vOZua6ZuU5zPe+3WzfN5/2+Pp/3mJt6er/fn9fbnHOIiIiISOI1S/UARERERLKFgpeIiIhIkih4iYiIiCSJgpeIiIhIkih4iYiIiCSJgpeIiIhIkih4iUiDmNmhZrYo1eNoCsxsgZkdkepxiEjiKHiJSEzMbImZ/S76unPuPedcl1SMKZqZ3WhmlWa22szKzey/ZnZwqscVK+dcV+fcv1M9DhFJHAUvEclIZtY8oOl559w2wE7ADODFBDzbzEx/fopI3PQHh4g0iJkdYWbLIr5eYmb/z8w+NrOfzex5M9s6on2gmc2NmJHaL6JtlJn9z8x+NbNPzezkiLZzzWymmd1lZj8CN9Y2LufcRuBpoMDM2oTusb2ZPWJm35pZmZndamY5obYcM7vDzH4ws8VmdpmZuXDAM7N/m9lfzWwmsBbobGZ7m9k0M/vJzBaZ2akR4z0u9D38GnrW/wtd38nMJoe+/5/M7L1wiIucVTSzFmZ2t5ktD/2628xaRP47N7PhZrYi9P2cV7+foIgkk4KXiCTCqcAAoBOwH3AugJn1Ah4FLgZ2BB4EXg0HCuB/wKHA9sBNwFNm1i7ivgcBXwFtgb/WNgAz2wr4PfAjsCp0+TFgI7AH0As4Brgg1HYhcCzQE9gfKPK57dnARcC2wEpgGvAMsDNwOvAPM9s31PcR4GLn3LZAN2B66PpwYBnQJvR9XAv4nd12HdAnNJ4ewIHA9RHtu+D9eyoAzgfuN7PWwf9GRCQdKHiJSCLc65xb7pz7CXgNLzyAF1oedM596Jyrcs49DqzHCxg4514MfW6Tc+554Au8wBG23Dn3d+fcRudcRcCzTzWzcqACL0yd4pzbaGZtgeOAK51za5xzK4C78AITeGHxHufcMufcKmCsz70fc84tCM2mDQCWOOcmhsZTCrwEDAn1rQT2NbPtnHOrnHNzIq63A3Z3zlWG9sj5Ba+hwM3OuRXOuZV4QfTsiPbKUHulc+51YDWQFnvtRCSYgpeIJMJ3Eb9fC2wT+v3uwPDQMlt5KCDtBrQHMLPfRyxDluPNFO0Uca9vYnj2C865fLzZpE+A3hHPzgW+jbj/g3izVYTGEHl/v2dFXtsdOCjqexmKNxMFMBgv6C01s/9EbPIfD3wJvGVmX5nZqIDvoz2wNOLrpaFrYT+GAmBY5L9nEUlTQZtTRUQS4Rvgr865GsuEZrY78BBwFPC+c67KzOYCFtHNb2bIl3PuBzO7CCgxs2dCz14P7BQVWMK+BXaN+Ho3v9tGfS//cc4dHfD8j4CTzCwXuAx4AdjNOfcr3nLjcDPrBkw3s4+cc+9E3WI5XrhbEPq6Q+iaiGQwzXiJSDxyzWzriF/x/s/bQ8AlZnZQ6M3AVmZ2vJltC7TCCzYrAUKbxbs1ZLDOuUXAVOAa59y3wFvAHWa2nZk1M7PfmNnhoe4vAH8yswIzywdG1nH7ycBeZna2meWGfh1gZvuY2VZmNtTMtnfOVQK/AJtC39dAM9vDzAz4GagKt0V5FrjezNqY2U7AX4CnGvLvQ0RST8FLROLxOt7eqfCvG+P5sHOuBG/f1X14G96/JLTx3jn3KXAH8D7wPdAdmNkIYx4PXGRmO+Nttt8K+DT0/El4+63AC4VvAR8DpXjf60a8YOT3vfyKtzn/dLyZqO+AcUD4RYGzgSVm9gtwCd4yJMCewNt4e7LeB/7hnJvh84hbgZLQeOYDc0LXRCSDmf+eThGR7GZmxwL/dM7tnuqxiEjToRkvERHAzPJCtbeam1kBcAPwr1SPS0SaFs14iYgAZtYS+A+wN94y6hTgT865X1I6MBFpUhS8RERERJJES40iIiIiSaLgJSIiIpIkGVFAdaeddnIdO3ZM9TBERERE6jR79uwfnHNt/NoyInh17NiRkpKSVA9DREREpE5mtjSoTUuNIiIiIkmi4CUiIiKSJApeIiIiIkmi4CUiIiKSJApeIiIiIkmi4CUiIiKSJApeIiIiIkmSEXW8RERERCIVl5YxfuoilpdX0D4/jxH9u1DUqyDVw6qTgpeIiIhklOLSMka/PJ+KyioAysorGP3yfIDN4Stdg5mWGkVERCSjjJ+6aHPoCquorGL81EXAlmBWVl6BY0swKy4tS8Foq1PwEhERkYyyvLyi1ut1BbNUSljwMrOtzWyWmc0zswVmdlPo+mNmttjM5oZ+9UzUGERERKTpaZ+fV+v1uoJZKiVyxms90M851wPoCQwwsz6hthHOuZ6hX3MTOAYRERFpYkb070Jebk61a3m5OYzo3wWoO5ilUsKCl/OsDn2ZG/rlEvU8ERERyQ5FvQoYM6g7Bfl5GFCQn8eYQd03b56PDmYXfziJJeMGMnP0UTBsWIpG7THnEpeFzCwHmA3sAdzvnBtpZo8BB+PNiL0DjHLOrff57EXARQAdOnTovXTp0oSNU0RERJqW4tnfsOnSyxj04avVrm9qnkuzyg0JfbaZzXbOFfq1JXRzvXOuyjnXE9gVONDMugGjgb2BA4AdgJEBn53gnCt0zhW2adMmkcMUERGRpmLdOhg4kKLCDjVCF8ATvQem9O3GpLzV6JwrB2YAA5xz34aWIdcDE4EDkzEGERERacJWrYL99oO8PJgyxbfLhYOu58Yjzk/p240JK6BqZm2ASudcuZnlAUcD48ysnXPuWzMzoAj4JFFjEBERkSbu66+9wPXzz4FdBp01njkF+2z+OpVvNyaycn074PHQPq9mwAvOuclmNj0UygyYC1ySwDGIiIhIUzRvHvTsGdzerBmnXTWRD5vvWKMplW83Jix4Oec+Bnr5XO+XqGeKiIhIEzdtGhxzTHB7p07wwQew886cUVrGxxFHC0H1shOpoLMaRUREJO1En7V4b+V8ev/lyuAPHH64t7erVavNl8LlJdLpzEYFLxEREUkrmw/B3rCRYR+8yDXvPhHc+eyz4dFHoXn1SBMd3O46rWdaHJKt4CUiIiJp5Y43PuX6yfcwdO6bwZ2uvx5uvhnMajRtDm6hJcbwIdlAysOXgpeIiIikhVff/5KdzzmT9774KLBP6XVjuWybA7yZrHEzfJcOazskW8FLREREsttPP/FLYR9OXPxFYJc/DP4Lc7r3Zf3GTVSEykEEzWSl8yHZCl4iIiKSGkuWQLdusGYN2wV0OensO5jX3jt7cWsjppms9vl5lPmErCZ9SLaIiIiIrzlzvL1ZnTrBmjU1mtfnNOewix6i48jJzGvfhYL8PAb3LmDV2krf20XPZEUfkg2Qm2OsWb+RTqOm0Hfs9JQdG6QZLxEREUmOqVNhwIDA5sWt2zH4rNv5qeX2m68VhEpAhJcU/UTPZEWXkchvmcvqdRspr/CCWyo322vGS0RERBJr4kRvhisodPXrx2szv+C4Sx+tFrrCxU79NstH94lW1KuAmaP6sXjs8bTcqjmVm1y19vASZbIpeImIiEjjc25LuYc//MG/z3nnQWUlvPMOJ/x2D8YM6k5Bfh6GN9M1ZlB3inoV1LopPtynNum02V5LjSIiItJ4qqpYcspZdCx+LrjPDTd4v6JqcBX1KvANUUGb5Qvy82JaKkynzfaa8RIREZGGq6jwzlBs3jw4dD38sDcTduONvoVPg/htlo/nzMWGfr4xacZLRERE6u+HH6BvX/j888Au555yI1/0PpSZ5/er1yMaeuZiOp3ZqOAlIiIi8fvqK9h3X1i/PrDLwHPu5pNd9gDAGrifKmgZMlmfbywKXiIiIhK7khI44IDA5nW5LTj6D/fzTf4u1a43dD9V9KHXqZqxaigFLxEREanb66/D8ccHNn+1464sfOF1NrTegR9eng8R5R8aup8qnQ+9jpeCl4iIiAR75BG44ILA5nc79uLCQdezPrcFBbNWMnNUd6Bx91Ol86HX8VLwEhERkerCbx7efHNgl2d69Of6Y4axqdmWtwXDdbEaez9VOtXhaigFLxEREfFs3OjNbj3+eHCfW26hb87BlP28rkZToupipVMdroZSHS8REZFst3Yt9OsHubnBoWviRG8m7PrrGTFg76TWxUqnOlwNpRkvERGRbLVyJfTp45WGCHDVeWM5/PKzqi0dJrsuVjrV4Wooc87V3SvFCgsLXUlJSaqHISIi0jR8+SXss4+3tBjg+HPvYUHb3wDe7FIsZyKKx8xmO+cKfdsUvERERLLErFlw0EHB7dtsw+BhDzLbtq/RlJ+XS6sWzTN+xikZagte2uMlIiLS1L32mnc2YlDo6toVfvwRfv2VOT6hC6C8opKy8gocXh2tq56fS8dRU+g7djrFpWWJG3sTo+AlIiLSVD34oBe4TjzRv/3YY73DrT/5BHbYAYj9TcHwelm4mKnCV2wUvERERJoS5+C667zAdckl/n0uuQSqqrxq9FtvXa3J7w3CuoSLmUYrLi2j79jpdNLM2GZ6q1FERKQp2LgRzjsPnnoquM+YMTBqVK23iX6DMNad4NHFTJvSMT+NScFLREQkk61eDccdB++9F9zniSfg7LNjvmVk5fm+Y6f7Fi+NFr1E2ZSO+WlMWmoUERHJRN9/Dx07wrbbBoeuadO8pcc4Qle0WJYe/YqZNqVjfhqTgpeIiEgm+fxzaNYMdtkFli717zN3rhe4fve7Bj+uqFcBg3sXYFHXw18X5Of51vgK2qSficf8NCYFLxERkUzwwQfehvkuXbxQFS0/3wtizkGPHo366BkLV9bY6+XwQtfMUf18lw6b0jE/jUnBS0REJJ298ooXuA4+2L+9Rw/46SdYtQo6dEjIEOqzbFjUq4Axg7pTkJ+HETwzlm20uV5ERCRGxaVlyTsv8B//gEsvDW4/4QR44YUa5SASoX1+nu8G+7qWDSM36YtHM14iIiIxCJdHiKze3uiFQ53zyj2YBYeuSy/1anC9+mpSQhdo2bAxKXiJiIjEoLbyCA1WWQlnnOFtmh83zr/PuHFeMLvvPq9fEmnZsPFoqVFERCQGCSmP8Ouv0L8/vP9+cJ9nnvFCWYpp2bBxKHiJiIjEoL77nHx99x0UFkJZLcuU06fDkUfGf29Ja1pqFBERiUGj7HNauNDbv9WuXWDomv58qOhpPUKXzkZMf5rxEhERiUH0GYZxvdU4cyYcckhg809523Hcuffy3XY7kTd/I2P2LIt7WU9nI2YGc35F2NJMYWGhKykpSfUwRERE4vPyyzB4cGDz/La/4cwz/savLVpVux4uTBqPoDMV63MvaRgzm+2cK/Rr04yXiIhIY7v3XvjTnwKb39zrYK444Ro2NM/1ba/Phv2gg6yz/WzEdKPgJSIiEtKgAqnOwYgRcMcdgV1e+O0gRh5yLs5q32Id74b94tIyDGoc61Ofe0liaXO9iIgIDSiQumEDnHqqV1srKHTdfjs4x1b33cvWW/nPcoXVpzDp+KmLfEOXgYqcphkFLxEREepRIPWXX6BPH2jRAl580b/Pc895M2HDhwP+hUjP6tOhwYVJg5YTHdpYn2601CgiIk1arMuHMRdIXb4cevf2anEF+fe/4fDDfZsSUYg0qMZYgZYZ046Cl4iIZIx492DFU2KhzgKpn34KXbvWPsBPPqm7Tz2/l9qM6N+l2vcJOksxXamchIiIZIToEAVeuKhtaS6oxELrlrm03Kp5tdAD+N5/Qse1HHpBcEkI2raF2bOhIPbQVJ/vJZZ7NlaQk4aprZyEgpeIiGSEWOtURQaQWP+GC4ce2FIgdeg3s7j1mZuDP3TggTBtGmy3XTzfBqCaW02d6niJiEjGi2UPlt9MUiwqKqsY/sI8NjnHVZ9M4YopDwT2ndKlL7NuuZebhuwf1zOCxhzLdWk6FLxERCTlYlkmi+WQar83E2PiHNe+8xDnl7wS2GXCASfztyP/AGbkzPmOm4bE/5jw9xk0E6eaW01fwoKXmW0NvAu0CD1nknPuBjPrBDwH7AjMBs52zm1I1DhERCS9xboBPpYN5HXNGBXk57Fm/UbKKyoByK2q5O+v3saAz98P/MxNR13IxMKTql2rqsc2nbpm47QZPjskcsZrPdDPObfazHKB/zOzN4Crgbucc8+Z2T+B84HgOV0REWnSaqufFRm8YjmkOmhWDLbsnyouLeOvz3zAI0+NYr/vvgwc17CTRjF1n0N9Q1aOWVzfY3jcQaGrQJvhs0bCgpfzdu2vDn2ZG/rlgH7AmaHrjwM3ouAlIpK14tnvVFcNrBH9u3Dl83OD71dWRtExPSn64YfAeww5cywf7daNgvw8zti7DU998HWNPmcctBsQ35uEQd+ngTbUZ5GEVq43sxwzmwusAKYB/wPKnXMbQ12WAYr3IiJZLGhfU332OxX1KqB1y5pH8uy1cgmLxw2EXXeFgNB11AUP0HHkZD7ardvmZb/C3Xeo8RdlM6Bw9x3iPmKoMb9PyVwJDV7OuSrnXE9gV+BAYO9YP2tmF5lZiZmVrFy5MlFDFBGRFBvRvwt5uTnVrjVkv9MNJ3TdfL8+X3/MknEDeevRy/w7t28Py5dTPGcZ636zV41je8ZPXcSmqI9swls2jPeIocb+PiUzJeWtRudcuZnNAA4G8s2seWjWa1fA938NnHMTgAng1fFKxjhFRMRfIotzxrJ3K977FUx9hQNGXxrc6eCDYepU2HZb7zPt/M80rE/Zh6C2xv4+JTMl8q3GNkBlKHTlAUcD44AZwCl4bzaeAwS/uysiIikXz7E79dUY5xcWl5ax7LpbuOyNBzkgqNNpp8GTT0JuzeVIP3WVsKirvEW0RJzTKJklkUuN7YAZZvYx8BEwzTk3GRgJXG1mX+KVlHgkgWMQEZEGindJLek2beKr08+jaP9dueyNB/37XHMNbNoEzz0Xc+iC2pcHtXQo9ZHItxo/Bnr5XP8Kb7+XiIhkgERWWQ9awoxpaXP9ejj1VHj1VToH3P+OEy9n+Cv31nt8sSwPaulQ4qGzGkVEpFZB5wrm5+ViBqvWVm7++sYTu8YcPIIOih7cu4CXZpcFHyBdXg5HHAHz5gXe++KTr2XqXr/FgMVjj49pPCKNpbazGhP6VqOIiGQ+vyW13GbGr+s3bg5dAOUVlYx4cV5gOYVoQUuYz374je/1x59/D3bYAVq3Dgxdg4aOp+PIyUzd67eASjVI+tFZjSIiUiu/5ba1G6qHrrDKTa5GxfkgQUuV0ZXi916xmDcnXl7rvY656EE+b139mbk5pv1WknYUvEREpE7Rb+N1GjUlsG+se7+C3hjMMaPKOQ5eOo9nn7su+AYdOsCsWfSduMD3Pq22aq79VpJ2FLxERLJULBvYg/rUdiZirMt7QYde37x6LkPuHh38wUMPhddfh222AWB5uf8e4J8ras7IiaSagpeISBaKpTZXbX1G9O/CiEnzqKyqviyY2yz25b3oJcxr5hbzx6kPB39g6FCYOLFGOYi6am2JpBNtrhcRyUKx1OaqrU9RrwLGn9Kj2rmI+Xm5jB/SI67lvaIe7Zj59SQWjxsYGLru7zOEjte8Rs89f0/xJytqtKuelmQSzXiJiGShWGpz1dWnQVXY162DwYO9JcMAfz5mGE/2Om7z1+G3JsPPDtNRPJJJFLxERLJQLMtzCVnCW7XK26O1YEFglwsG/Zm39zzIty3orUkdxSOZQkuNIiJZKJbluUZdwlu6FLbbzqvDFRC6is6+g44jJweGrrDGqJgvkiqa8RIRyULRy3P5LXNxDq56fi7jpy6qtlTXoCW8uXOhV43T47bIyYHPPqPTI58T6zkqQTNuMR0zJJJiCl4iIlkqvDxX1xuO9Qov06bBMccEt3fqBB98ADvvDED7/G8Cy1NECnprMpa3NEXSgZYaRUSy3I2vLqjzDceYPfEEmAWHriOOgNWr4auvNocu8F/WjFbbW5OxvKUpkg404yUi0sTEs+RWXFpGeUCh0Zj3UjkHf/sbXH99cJ9zzoGHH4bm/n/thMc3/IV5NY4MAijIz2PmqH6Bt4/lLU2RdKDgJSLShMS75FbbjFCdby9WVcGwYTBhQnCfP/8ZbrrJmwULGG9kSDzjoN14aXZZjWr2dW3oVxFVyRRaahQRaULiXXKrbUYoMOxUVMCAAd7sVVDomjDBmwm7+eZaQ9fol+dTVl6BwwuJL80uY3DvAgry8zC8ma4xg7rXuU9LRVQlU2jGS0SkCYl3yS1opqh1y9yaYefHH+GQQ2DhwuABTJ4Mxx8f01iDQuKMhStrXVb0oyKqkikUvEREmpB4l9yCDqq+4YSuWzotXgxdu3ozXUFmzYIDDohrrI29L0tFVCUTaKlRRKQJiXfJrahXAWMGdfdf2pszx1sm7NzZP3S1aAFffuktKcYZuiA4DGpfljRlCl4iIk1IrUEqRq2+XgwXXwy9e/t32HNPWLHCO2/xN7+p91i1L0uykZYaRUSamHiW3CLfguz6/f/44yuT6LdoJrhNNfp+tMf+/P7Ea9mhTWtGLNtAUZuGjxO0L0uyizmfeinpprCw0JWUlKR6GCIiGc2vvtf4Nxey2/xZ/PGDSRy+eI7v55aedBrHdT2bNRH74HObGdts3ZzytZUKTCJRzGy2c67Qr00zXiIiWSC6vtfyVWuYNmYCf5/5PPsvr1lq4j+d9ufwiXfCYYdx5rgZrIna8F65ybFqrVd4VcfziMROwUtEJAuESzc0r9rIiZ/9h0s+eIm9fvy6Wp8qa8brXfryzz6nUN6lGzMPPxyI7S3DcK0wBS+R2il4iYgkWDxH+CTKTytXcc7H07hw1svs+svKam0bcpozqdvvePCgQSxt3Z683BzG9O+yedyxbkjR8TwidVPwEhFJoHiP8Gl0q1bBfffx3wfvoPWan6s1rd4qj1f6nMCO113D/XN+Znl5BQWhYAjUqO9VF5WBEKmbgpeISALVdoRPIoJXeJaq6ptvuGL+FIbMfp3ctWtoHdHnh5bbM7H3iUw68ARGD/0tA3oVMOCY6vfpO3Z6YOjKz8tlzYaNVFZtmQtTGQiR2Ch4iYgkUGNXZ69NcWkZDz70BpfNfJFBC6bTompjtfa17Xbln4Un89Ceh7NDm9aMrmXJM2h8Bsy94Zi0WD4VyUQKXiIiCRR0hE8zMzqNmtJ4oWX2bLY792qmfPIezaJ2ZX3VtiOdb7+FlqedxtW5uVxdy23q2tcVXk7U8Twi9aPK9SIiCeRXnR2gyjkcW/Z8FZeWxX9z52DGDDjmGCgspN8n71YLXbPb7835g//M7865F846C3Jza71deD+aX1AELSeKNAbNeImIJFB0dfZmZlRFFa6Oe8/Xpk3wyiswdqx3OHWUGZ1780CfIczatSuYURDjpne//WhhBVpOFGkUCl4iIgkWuSzXadQU3z4x7fnasAGeeQbGjYOFC6u3NWvGsqNP4LIOxzB3h903X45nlqq2fV0zR/XzbdNeL5H4aKlRRCSJgkou1FqKYc0auOce2GMPOO+86qGrRQu45BL4/HN2fbOYc/94Ur0PyI53bJFLkw1eNhXJEprxEhFJohH9u9SojxU4K/XTT3DffXDvvfDjj9Xbtt0Whg2DP/0J2rXbfLkhm97jGhvJL5Uh0hQoeImIJFH0ni/f5blly+DOO2HCBG+2K9LOO8OVV8If/wj5+ckfW4RklsoQaSoUvEREkixwVmrRIrjtNnjySaisrN7WsSOMGOEtNeYlrkJ8PDNmQaUyVMFeJJj2eImIpFpJCZxyCuyzDzz6aPXQ1b07PP00fPGFt7SYwNAVL79SGSo5IVI7zXiJiKSCczB9OowZA++8U7O9b18YPRqOOw7Mkj++GMS7NCkiCl4iIsm1aRMUF3s1uD76qGb78cfDqFFwyCFJH1p9qIK9SHwUvEQkY2R0zagNG+Cpp7w9XIsWVW9r1gxOPx1GjoT99kvN+EQkKRS8RCQjhGtGhcsXhGtGAekdvlavhocegjvugLKo+lYtWsD558Pw4dC5c2rGJyJJpeAlIhkh42pG/fgj/P3v3q+ffqrett12cOmlXg2utm1TMz4RSQkFLxHJCBlTM+qbb7bU4Fq7tnpb27Zw1VVepfntt0/N+EQkpRS8RCQjpH3NqIULvTMUn3oKNm6s3ta5s1eD65xz0qochIgkn+p4iUhGSNuaUbNmwaBBsO++8Nhj1UPXfvt5h1ovWuTNcil0iWQ9zXiJSEZIRM2oeN6SrNZ3+625LX8FfSc97NXiinbooV4NrgED0rYGl4ikhjnnUj2GOhUWFrqSkpJUD0NEmpDotyTBm0Eb3LuAGQtXVgtjAKNfns/69Rvo//n7/PHDSez33Zc1bzpwoFeDq2/fwGdmbDkMEYmZmc12zhX6til4iUg26jt2uu+eMQMi/1TMy81hW9vIER9N4+JZL/Gbn6qXhNhozXi7Rz+2Gj2Sfqf+LvB5QUFvzKDuCl8iTUxtwUtLjSKSFaJnm/xCF1QPXa3Wr+WMWW9ywUfF7LK6ekmIdc234vn9juahA05mWf4u5M3fyJg9ywJDVMaVwxCRhFDwEpEmz6/4avTMVqQd1v7MObNf45w5k8lft7pa2y8tWvHE/sczsfeJ/Ngqf/P1ukJUxpTDEJGEUvASkSbPb7bJUXNZseCXFVww61+cPu8t8jaur9Z/RavWPHxAEc/0PJbVLVr6Pqe2EJX25TBEJCkSFrzMbDfgCaAt3p9tE5xz95jZjcCFwMpQ12udc68nahwiIkGByAEF+Xm0/HIRV5UW03/edHKqqtfgWpLfjgcPGszL3fqxvvlWtT4nOkRFLm/mt8wlt5lRuWlL1EuLchgiklSBwcvMXgeGOeeW1PPeG4Hhzrk5ZrYtMNvMpoXa7nLO3V7P+4qIxCVotunoX5bw0MJp8MorNdoW7NyZ5486k6d2PYBNzXJqtEeLDlHRy5ur1laSm2Pk5+Xyc0Wl3moUyVK1zXhNBN4ys8eB25xzlfHc2Dn3LfBt6Pe/mtlngP6EEZGkG9G/y5YQ5ByHLZ7DpbNe4qClH9fsfPjhMHo0XY85hpvNeCfg7cfWLXNpuVXzwNIQfsublVWOVi2aM/eGYxr9exSRzBAYvJxzL5rZG8CfgRIzexLYFNF+Z6wPMbOOQC/gQ6AvcJmZ/R4owZsVW1W/4YtIJktWXauiXgVQVcXsOx/mtOnP0O37/9XsdOKJXg2ugw+udrlaaAvJy83hhhO61jpWbaYXET917fHaAKwBWgDbEhG8YmVm2wAvAVc6534xsweAW/C2V9wC3AH8wedzFwEXAXTo0CHex4pImisuLWPEpHlUVnl7nsrKKxgxaR5AXOGrzvC2fj088QRFt91G0ZdRRU9zcmDoULjmGuja1ff+9a2Yr830IuInsICqmQ0A7gReBW52zq2N++ZmucBkYKrfDFloJmyyc65bbfdRAVWRpqfXzW+xam3NHQytW+ZS+hdvKa6uUFVrUdI9toMHH4Q774Rvv63+kLw8uOACGD4cdt89Id+fCqaKZK/6FlC9DhjinFtQz4ca8AjwWWToMrN2of1fACcDn9Tn/iKS2fxCV+R1v9pbo1+eD1SfhYreR5X380+UDx8JpVOgvLz6zfPz4bLL4IoroE2bxvtmfCTibEkRyXy17fE6tIH37gucDcw3s7mha9cCZ5hZT7ylxiXAxQ18jog0QbFUeo/cL7Xrz99zwax/cdrH02rU4KJdO7j6arj4Yth224SPPayoV4GClohUk7A6Xs65/8OrTxhNNbtEhPy8XMoras565eflArFtTm+fn0erLz7j4g9f4qRP/0NzF7UNdY89KD3jIq7aqjtLf6ii/f0fadZJRFKqWaoHICLZ6cYTu5LbrPr/m+U2M2480dvkHrQJffP1999n0htjeevRyxi8YEa10FW+dzd44QWKn32HM113lqyuwuEtV171/Fw6jppC37HTKS4t832GiEiiKHiJSEoU9Spg/JAeFOTnYXgV5McP6bF5NmpE/y7k5VYvXJrXvBnjt1nu1dr67W9p9+60au2zO/dk5v1Pk//pxzBkCOPf/tL3qCDYsmdM4UtEkklnNYpIytS2Bypyc/r3P63mzG9mMXzOv9j+8099OhfByJH07tOn2uW6ambVdbC1iEhjU/ASkbRVtM+OFM36Ah4eD/+LKnravLlXg2vkSNhnH9/PB9XSiqSCpiKSTApeIk1QsirCJ8wvv8A//wl33QXffVe9rWVLuPBC7y3FOoor+1Wdj6aCpiKSTApeIk1MLPWv0taKFXDPPXD//fDzz9XbWreGyy/3fu20U0y3i1yuLCuvwNiyxwtqHmwtIpJoCl4iTUws9a/SzowZ8Pe/wxtvwLp11dvat/cqzF94Yb1qcEXuI8v4mUARyXgKXiJNTEYdzvzss3Dmmf5te+1F6ekXctGmvVm5wuCv79K6ZW6dh1PXRgVNRSTVFLxEmpiMOJz5ttu8TfE+FhbsxT0HDGZWj0NZtW4TkSVRV62trNdB2iIi6UJ1vESaGN/6V+mwl2nTJu+cRLPA0HXLgGEMGHoHb3Tpy49RoSusssoxfuqixI5VRCRBNOMl0sSk3eHM69fDkCHw2mvBfe67j76/7lNn6YewtFw2FRGJgYKXSBOUFnuZysu9CvMffxzcp7gYTjoJgOWjpsR867RaNhURiYOCl4g0rq+/hv32q1kOIsKgs8bzfdf9GdGhC0Wha7EUOwXIzbHUL5uKiNSTgpeINI5586Bnz8BmZ8Zxl0zgs+3aeRei6ov5FTvNzTGaNzMqKr3dXg19q1FEJNUUvESkYd5+G44+Orh9993hww85ZOKCGjNakfXF0m5vmohIAih4iUj9PPUUnH12cPuhh8Lrr8M22wCwvLzEt1vkRvm02JsmIpJAKichIjUUl5bRd+x0Oo2aQt+x0ykuLdvSOGaMVxIiIHS93PVIul07meJ7nt0cuiB4Q7w2yotINlHwEpFqwmc9lpVX4PDOerz2pXksPuVsL3Bde63v5/5+8Gl0vOY1rh44nNVV1Ki1lbb1xUREkkhLjSIZIlnnDEae9dhi4wb++a+/cuRXswP7X9v/Up7peWyN69G1trSHS0REwUskI4RnocKBqCzqjcBYPh8ZeI7cuw0zFq70DUDLyyvYvuJXXnx6JHv9+HXwTV99leJd9+f5F+aBczWa/ZYQtYdLRLKdgpdIBoichQqLfCOwNn6h7akPtgSqaiFuh40suGsILTfUUk/rww/hwAM337fKJ3RpCVFExJ/2eIlkgKAjcmI5OscvtEXrvOxzivbfFTp29A1dlc1ymFb8njezdeCBtd43x4wxg7prZktExIdmvEQyQFBV91jeCKwtnB26eA5PvvCXwPYl+e249NL7uHDwQTWCVNB9Nzmn0CUiEkDBSyQD+FV1j3U5zy+0DZ7/Dne8flfwh448El57jY6tWhF0gmJDwqCISLbSUqNIBijqVcCYQd0pyM/DgIL8vJiX8zaXcXCOy/77HEvGDQwOXeeeC5WVMH06tGoV230jaG+XiEjtzPlsjE03hYWFrqTEv+q1iNShqoolQ86m47+eDeyy8KKr2Pufd3h1uuKQrBIXIiKZxMxmO+cK/dq01CjSBBWXlnHP5I+56dHrOGxJKR0D+o0acDnvHV7EiP5d2DvO0AUqDyEiEi8FL5EMEssM05QZn7DfkGOZ8eOywPtcfPpNTN29t/dFLTXBNKMlItK4tNQokiGi63EBGDC0TwduLeoOixfDvvvCunXBN/noI/q+/YvvpviC/DxmjupX6/PycnNUKkJEpA61LTVqc71IhvCrm+WAecXTvb1ZnTv7hq71ObkcevHDdBo5GQoLY64JVlvRVhERqR8tNYrEKVXLb9HB6Ij/lfDYpBsD+/9vhwJOGXobq1puD3gzWhB7GYiGFG0VERF/Cl4icWjomYkNEQ5Mp857i9vevDew37sde3HxoOuoyN1687XIMg+x1gRTnS4Rkcan4CUSh4acmdggzvHI4snsPSG46Omknv255uhhbGrm1dYyvKXIgqhZufA/65q1a0jRVhER8afgJRKHpC+/VVXBBRfAY4+xd0CXew8/m4lHDmVVxcZq18OhK3LDfFgsZSBiDWgiIhI7BS+ROCRt+W3tWjjhBK+CfIBrjr2CmYedxIj+XSh/fq5vn4YGQtXpEhFpXHqrUSQOiTomp7i0jL5jp7P/Fc/wzU67esf1BIWuN98E57jt9XuYOaofRb0KAoNfMzOKS8saNDYREWk8mvESiUNty2/1fduxuLSMfzw8lRn/vJitNm0M7HfeZQ9w0h9O8L2n334sgCrnkrb5X0RE6qYCqiKNoN7FRj/6CA48MLB5Te7W9D//fpZt37bOexaXljH8hXlU+fw3HbTXS0REGp8KqIokWNzFRqdM8YqeBoSuz3fsQI8rnqXr1ZM2h6667lnUq4BNAf8jpdpbIiLpQUuNIo0g5rcdJ0yAiy8OvM+/O/Xm4pOvZX1ui7ifBaq9JSKS7jTjJdIIgoJNfstccA6uv96b4QoIXU/3HEDnEa9w7qk31Rq6ansWJG7zv4iINA7NeIk0ghH9uzBi0jwqq7Ys9eVsquIvL94JN/QP/Nxth/2efxx8qtffDHOO7fNyWbNhY7V7hdUVolR7S0QkvWlzvUgj6XnTW5RXVJK3YR0TJ91In28+Cex79fFX8XK3o6pdM2Dx2OOBLedBlpVXkGNGlXM1KtCLiEh6qm1zvWa8JOsk6pDr5j+s5L0nh7Pbz98H9vnT+bfxyk77+rZFLiGqcKmISNOk4CVZJSGHXH/xBXTpQklts8elpdCzJ0eWlvH6i/Oo3FS9b26OaR+WiEgW0OZ6ySpxl32ozQcfeBvm99rL20Af5ZcWrTjqssconrMMevYEvHA3fkgP8vNyN/dr3TKX8af00AyXiEgW0IyXZJVGOeT6lVegqCiw+ctdOjH4tDFss0sb32XMyGXE8LLnVc/PZfzURdrDJSLSxCl4SVZpUJ2rBx6AYcOC248/HiZNYo+tt2ZeDGOp77JnovaoiYhI4mmpUbJK3HWunINRo7wlxaDQNWwYVFXB5Mmw9dYxj6U+y57hsFZWXoFjS1jTQdgiIplBwUuySlGvAsYM6k5Bfh6Gd4ah79mHlZVw5pnQrBmMG+d/s7FjvWB2//1evzjVZ9mzUfeoiYhI0iVsqdHMdgOeANoCDpjgnLvHzHYAngc6AkuAU51zqxI1DpFotZZqWL0aBgyAmTODb/DUUzB0aIPHUZ9lz0bZoyYiIimTyBmvjcBw59y+QB/gUjPbFxgFvOOc2xN4J/S1SGp99x106ADbbhscut5+25vhaoTQBfU73icolOksRhGRzJCw4OWc+9Y5Nyf0+1+Bz4AC4CTg8VC3x4GiRI1BpE6LFnn7t9q1g2++8e/z8cde4DrqKP/2eop52TOCzmIUEclsSXmr0cw6Ar2AD4G2zrlvQ03f4S1FiiTXzJlwyCHB7Tvu6BU93W23hA4j3gr1OotRRCSzJTx4mdk2wEvAlc65X8xsc5tzzpmZb7lvM7sIuAigQ4cOiR6mNAExlVn4179g0KDgm/TqBdOnQ35+QsfaEDpOSEQkcyX0rUYzy8ULXU87514OXf7ezNqF2tsBK/w+65yb4JwrdM4VtmnTJpHDlCagzjIL993nLSkGha6TToJ162DOnLQOXSIiktkSFrzMm9p6BPjMOXdnRNOrwDmh358DvJKoMUj28C2zsGEja6682gtcl1/u/8ErrvBqcBUXQ4sWiR+oiIhktUQuNfYFzgbmm9nc0LVrgbHAC2Z2PrAUODWBY5AsEVlOoXnVRu6efAcDF74X/IHbb4fhw5MwMhERkS0SFrycc/8HWEBz474eJlmvfX4e5d//yFPP/5le39ZSTPS55+C005I3MBERkQg6q1Ey37ff8vYdZ5L3w/fBfWbMgCOOSNqQRERE/Ch4Seb67DPYd18AAsuHzp8P3bolbUgiIiK10VmNknnee8/bMB8KXdFWtszn8CuepHjOMoUuERFJK5rxkpSIqeZWtBdfhFOD38WYt8uenHX6rfzaohXgvemoelciIpJOFLwk6cI1t8LlH8I1twD/oHT33XDVVYH3e2Ov33LFiSOozMmtdj3yTcd6BT0REZFGpqVGSTrfmluVVYyfGvE2onNwdagGV1DouuoqqKri1vNuqRG6YMvB0XUWVxUREUkSBS9JusiZqEhl5RXsNaKYGd0Pg2bN4K67/G9w111eMLvzTmjWrM6Do2MKeiIiIkmgpUZJuvb5eZRFha9t1q/l6eeuo8d3XwR/8MUX4ZRTql0KLyFWVFaRY0aVcxRELSUGBb2g6yIiIomi4CVJN6J/l817vNr++gOTH7uSNmvLgz/w7rtw6KE1LkfvFatybvNMV+T+Lb+gF74uIiKSTFpqlKQr6lXAfT22Ysm4gXz4j3MDQ9fR5//DW1L0CV0Q+xJiXUuRIiIiyaIZL2l0tb5B+J//wBFHBJ4Z9f02O3DC7+9ixbY7UlDHjFSsS4jhZ+utRhERSTUFL2lUQaUiCqa+ygGjhwV+bm7B3gwdcjNrWrQEYpuRimcJsahXgYKWiIiknJYapVFFL/9dMOtlPrv12ODQddppsGEDS157m/y2O2JAQX4eYwZ1rzMoaQlRREQyjWa8pFEtL6/A3Cb+8s5DnDf7tcB+Dx44iCeKhjFiwN4U5ebWa0ZKS4giIpJpzDmX6jHUqbCw0JWUlKR6GFKX9et5d/9+HPbpfwO73Nr/Eh7uOXDz13m5OTHNbomIiGQKM5vtnCv0a9OMlzRceTn06welpRwW0OXD8RO4euNvauzJCr+FqOAlIiLZQHu8pP6WLYMdd4TWraG01LfLxZfcS/GcZRz0/y5UIVMREcl6mvGS+H38MfToUXufhQuhSxcejLikQqYiIpLtNOOVAYpLy+g7djqdRk2h79jpqTvcefp079DqgND13fY788Zbc7yip11qvlmotxBFRCTbacYrzQXVxQKq7YuqtWhpHH18PfMMDB0a2PxRwb6cO+RG1rRoSd6737N+pzLf++otRBERyXZ6qzHN9R073Xd5riA/j5mj+gE1wxnUfFswlj6RikvLWH7tjQx78+HAsb3Vox/Djr6CjTnV83vk2ERERLKN3mrMYEEbz8vKK+g0agrt8/NYu2Fj4JmFkbNMdfUBYNMmvjr9DxS9+HjwoEaNgr/9jYtHv45fbNdmeREREX8KXmkuaEM6gIPANqgegOp8o3DdOjjlFJgyhc4B9xt/0p8YUXx3nWMLb5av99KmiIhIE6XN9WnOb0N6rCLfFtw+L9e3z66sg/32g7w8mDLFt8+Fg66n48jJ/GPvo+scW3izfHhps6y8YnNAHP3y/NS9GCAiIpIGNOOV5qI3pMe6Iy/6bUGz6u3tf1nBm49eznbr1wTeY9BZ45lTsM+Wz0SVfahts3zfsdNjW9oUERHJIgpeaSpomS5os31+Xi6tWjT3XdYrLi1j1dpKAPZZ8RVvTLwi+MHNmjHt5X9zRcmaGhvx/co+BJ2xqGKpIiIiNSl4paHaSkiM6N/F9+3EG0/sGvh24ogX53HI4lKeeuHPwQ/t1Ak++AB23pmjgTEdGrY/S8VSRUREatIerzRU1xuIYwZ1pyA/D8Mr3VDbIdMf3XIPX4w5PjB0/dC7D6xeDV99BTvv3Gjfg4qlioiI1KQZrzRU1zJd0PLeZs7BmDFw3XX8NaDLS12PpPnEiZx0wO412mIt2lobFUsVERGpScErDdV7ma6qCi69FB58MLDLvQefxp2HngVmLPEJXRBHza861BkQRUREsoyCVxoK2scVuExXUQEnnwxTpwbec3T/y3i254DNX7du6V9eArQxXkREJFEUvNJQzMt0P/0EhxwCn30WeK/3736M369oQ2XVlkIUuTnGDSd0DfyMNsaLiIgkhoJXmokuI3HXaT1rBq4lS6BbN1gTXIOLDz+EAw/kYGB8nBXk455xExERkZgoeKWROje1z5kDvXsH3yA3Fz79FPbYo9rlePdaaWO8iIhIYih4pZGgTe3v3vskRY+NDv7gHnvAf/8Lbdo02li0MV5ERKTxKXjFyK+SPDTurFD05vUhH09j/Bv3BH+gXz947TVo2bLezxQREZHkUfCKgd8S4IgX54GxedN6fWpdRWufn0fZqrVc/t/nGP5/Twd3PO88mDABmuvHJyIikkn0N3cM/JYAKzfVPK66QYdAV1Xx9IcP0bH4ueA+N9zg/Yo+8VpEREQygoJXDOKpX7W8vCLwgGtfFRVw0kkwbRodg2768MNw/vnxDltERETSjM5qjEE89avyW+Yy+uX5lJVX4NiyBFlcWla94w8/QJcu3v6sadN87zX83L9RPGeZQpeIiEgToeAVgxH9uxDL4l5ebg7OEXjcDuAdRr311t4biJ9/7nufgefcTceRk3mp7X7+oU1EREQykoJXDIp6FVBzR9cWBhTk5zFmUHd+rqj07bPjZx97e7N+8xtYv75Ge0XzFhx68cN0HDmZT3bZUoerWmgTERGRjJb1e7xi3Y9VEHCMTkF+HjNH9dv89fipi6r1O+J/H/HYpJuCB9ClC8cNGcOnlVsFdtEZiSIiIk1DVs94hctE1LkfC2+5MS83p9o1v2N0wv1OnfcWS8YNDA5dxxwDa9fCwoV8VkvoAp2RKCIi0lRk9YxXUKV4v5IQMR2j4xxFxRMouvXm4IdeeCE88ADkbAlxQYdSQ/AZiXG9OSkiIiJpIauDV9ASXtD1wGN0Nm6ECy6Axx8PfNYdh57Fy8eey4gBe1OUU33mzO9QaoDWLXO54YSuNZ5Z55mOIiIikpayOngFzTSFl/bqnFVauxYGDoQZMwKfMfL4K3m+2++8L35e5xuQ4j2UOp6ZOhEREUkf5lxt7+ulh8LCQldSUtLo942eOQJvaW/MoO4ANdoMcEDX5ut54fGrabVsaeC9fz/kJt7t3Nu3LXpDfrw6jZri+5alAYvHHl/v+4qIiEjDmdls51yhX1tWz3jVNtPUd+z0GrNKHVYt5+2H/0jupiq/2wFw/Dl3syCiHISf+r6lGJ6BC4rK2oQvIiKS3rI6eEHwvq3IcNRz+SKKnxwefJNttoH58+n0zwW11vsKq09A8pudixS0CV9ERETSR9YHryDt8/PYu+TfPPLSLYF9Fu3UgS6LSmGHHUKf+Srw7cSw+gYkv31dYQV6q1FERCQjJCx4mdmjwEBghXOuW+jajcCFwMpQt2udc68nagz1Nfe6ccz826jA9hmde3PJydeRt20r5oZCF/i/nZibY7Taqjk/V1Q2qOxD0PKkQYP2i4mIiEjyJHLG6zHgPuCJqOt3OeduT+Bz68c5uP56+Nvf6BnQ5amex/LnY/6IM6/u7KYNGykuLdscpOJ9OzEedb2BKSIiIukvYcHLOfeumXVM1P0b1bPPwplnBjaPO/wcHjx4CJuiNnBVVrkaJRwCa301kN9smvZ1iYiIZJZU7PG6zMx+D5QAw51zq1Iwhi1++QXOO8+36arjr+Zf3ULLeAG75svKK6rNeiVKImfTREREJDkSWscrNOM1OWKPV1vgB7wYcwvQzjn3h4DPXgRcBNChQ4feS5cG18xqkLVr2dCuPVv98vPmS0NPu5WZHXtW65ZjRlXAv6tw7a9wCNJxPiIiItkrbep4Oee+D//ezB4CJtfSdwIwAbwCqokaU/GiVTw49Ha6Ll3Agrad+WznzjX65OXmMLh3AS/NLvN9szBcNR7gptcWsGpt5eY2HecjIiIiYc2S+TAzaxfx5cnAJ8l8vp/xUxfx2XbtmNT9d76hy4DBvQu4taj75or2fsIBKzJ0hUUGMxEREcleCQteZvYs8D7QxcyWmdn5wG1mNt/MPgaOBK5K1PNjVVcVeQfMWOhVvyjqVUBBwFuEOWaBdbZieY6IiIg0fYl8q/EMn8uPJOp59RVUpiFSZGgKeruwttAVfo6IiIhkt6QuNaajEf27kJebU2ufyNBU1KuAMYO6U5Cfh+FVjQ9/HURlH0RERAR0ZFC1Mg1l5RUY1StH+IWmoFpdfmcp5uflcuOJXbWxXkRERBS8oHqQqm8pCNXZEhERkboktI5XYyksLHQlJSWpHoaIiIhInWqr45X1e7xEREREkkXBS0RERCRJFLxEREREkkTBS0RERCRJFLxEREREkkTBS0RERCRJFLxEREREkkTBS0RERCRJMqKAqpmtBJYm6PY7AT8k6N5Sf/q5pC/9bNKXfjbpSz+b9JSon8vuzrk2fg0ZEbwSycxKgqrLSuro55K+9LNJX/rZpC/9bNJTKn4uWmoUERERSRIFLxEREZEkUfCCCakegPjSzyV96WeTvvSzSV/62aSnpP9csn6Pl4iIiEiyaMZLREREJEmyNniZ2QAzW2RmX5rZqFSPRzxmtpuZzTCzT81sgZn9KdVjki3MLMfMSs1scqrHIluYWb6ZTTKzhWb2mZkdnOoxicfMrgr9WfaJmT1rZlunekzZysweNbMVZvZJxLUdzGyamX0R+mfrRI8jK4OXmeUA9wPHAvsCZ5jZvqkdlYRsBIY75/YF+gCX6meTVv4EfJbqQUgN9wBvOuf2Bnqgn1FaMLMC4Aqg0DnXDcgBTk/tqLLaY8CAqGujgHecc3sC74S+TqisDF7AgcCXzrmvnHMbgOeAk1I8JgGcc9865+aEfv8r3l8gBakdlQCY2a7A8cDDqR6LbGFm2wOHAY8AOOc2OOfKUzooidQcyDOz5kBLYHmKx5O1nHPvAj9FXT4JeDz0+8eBokSPI1uDVwHwTcTXy9Bf7mnHzDoCvYAPUzwU8dwNXANsSvE4pLpOwEpgYmgZ+GEza5XqQQk458qA24GvgW+Bn51zb6V2VBKlrXPu29DvvwPaJvqB2Rq8JM2Z2TbAS8CVzrlfUj2ebGdmA4EVzrnZqR6L1NAc2B94wDnXC1hDEpZLpG6h/UIn4YXj9kArMzsrtaOSIM4r85DwUg/ZGrzKgN0ivt41dE3SgJnl4oWup51zL6d6PAJAX+BEM1uCtzTfz8yeSu2QJGQZsMw5F54ZnoQXxCT1fgcsds6tdM5VAi8Dv03xmKS6782sHUDonysS/cBsDV4fAXuaWScz2wpvs+OrKR6TAGZmeHtVPnPO3Znq8YjHOTfaOberc64j3n8v051z+j/3NOCc+w74xsy6hC4dBXyawiHJFl8DfcysZejPtqPQiw/p5lXgnNDvzwFeSfQDmyf6AenIObfRzC4DpuK9ZfKoc25Bioclnr7A2cB8M5sbunatc+711A1JJO1dDjwd+h/Jr4DzUjweAZxzH5rZJGAO3hvbpaiCfcqY2bPAEcBOZrYMuAEYC7xgZucDS4FTEz4OVa4XERERSY5sXWoUERERSToFLxEREZEkUfASERERSRIFLxEREZEkUfASERERSRIFLxHJKma2m5ktNrMdQl+3Dn3dMcVDE5EsoOAlIlnFOfcN8ABe/R5C/5zgnFuSskGJSNZQHS8RyTqhY6lmA48CFwI9Q0e6iIgkVFZWrheR7OacqzSzEcCbwDEKXSKSLFpqFJFsdSzwLdAt1QMRkeyh4CUiWcfMegJHA32Aq8ysXWpHJCLZQsFLRLKKmRne5vornXNfA+OB21M7KhHJFgpeIpJtLgS+ds5NC339D2AfMzs8hWMSkSyhtxpFREREkkQzXiIiIiJJouAlIiIikiQKXiIiIiJJouAlIiIikiQKXiIiIiJJouAlIiIikiQKXiIiIiJJouAlIiIikiT/HxVZRpJfWYwEAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Train and test\n",
    "x_train, x_test, y_train, y_test = train_test_split(x, y, test_size = 0.3, random_state = 42)\n",
    "\n",
    "# Linear Regression\n",
    "#   Note:   It is needed to give both the fit and predict methods 2D arrays. \n",
    "#           Right now x_train and x_test are currently only 1 dimensional.\n",
    "#           What we shoudl do to fix it is: \n",
    "x_train = x_train.reshape(-1, 1)\n",
    "x_test  = x_test.reshape(-1, 1)\n",
    "\n",
    "# Create linear regression object\n",
    "model = LinearRegression()\n",
    "\n",
    "# Train the model using the training sets\n",
    "model.fit(x_train, y_train);\n",
    "\n",
    "# Make predictions using the testing set\n",
    "y_pred = model.predict(x_test)\n",
    "\n",
    "# The coefficients\n",
    "print('True model: y ~ N(μ, σ^2)')\n",
    "print('True μ =', beta_0, '+', beta_1, '* x')\n",
    "print('Regression μ = %.2f' % model.intercept_, '+ %.2f' % model.coef_[0], '* x')\n",
    "\n",
    "# The mean squared error\n",
    "print('Mean squared error: %.4f' % mean_squared_error(y_train, model.predict(x_train)))\n",
    "\n",
    "# The coefficient of determination: 1 is perfect prediction\n",
    "print('Coefficient of determination: %.4f' % r2_score(y_test, y_pred))\n",
    "\n",
    "# Plot outputs\n",
    "plt.figure(figsize=(10,5));\n",
    "plt.scatter(x_test, y_test)\n",
    "plt.plot(x_test, y_pred, color='red', linewidth=3);\n",
    "plt.title('Linear Regression');\n",
    "plt.xlabel('X');\n",
    "plt.ylabel('Y');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **3.** Linear Regression with `statsmodels`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.987\n",
      "Model:                            OLS   Adj. R-squared:                  0.987\n",
      "Method:                 Least Squares   F-statistic:                 2.620e+04\n",
      "Date:                Wed, 15 Jun 2022   Prob (F-statistic):               0.00\n",
      "Time:                        16:05:37   Log-Likelihood:                -369.94\n",
      "No. Observations:                 350   AIC:                             743.9\n",
      "Df Residuals:                     348   BIC:                             751.6\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         14.8912      0.074    202.316      0.000      14.746      15.036\n",
      "x1             2.0254      0.013    161.863      0.000       2.001       2.050\n",
      "==============================================================================\n",
      "Omnibus:                        0.403   Durbin-Watson:                   2.041\n",
      "Prob(Omnibus):                  0.818   Jarque-Bera (JB):                0.505\n",
      "Skew:                           0.071   Prob(JB):                        0.777\n",
      "Kurtosis:                       2.879   Cond. No.                         11.8\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "# Train and test\n",
    "x_train, x_test, y_train, y_test = train_test_split(x, y, test_size = 0.3, random_state = 42)\n",
    "\n",
    "# Add (explicitly) the constant to the model\n",
    "x_train = sm.add_constant(x_train)\n",
    "x_test  = sm.add_constant(x_test)\n",
    "\n",
    "# Fit and summarize OLS model\n",
    "model = sm.OLS(y_train, x_train)\n",
    "model = model.fit()\n",
    "print(model.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some statistical tests.\n",
    "\n",
    "* **Overall F-test for significance**<br>\n",
    "$\\quad$ `p-value = Prob (F-statistic) = 0.00` <br>\n",
    "The F-statistic test is used to test the significance of regression coefficients in linear regression models. If the p-value is less than a certain significance level (*e.g.* 0.05, 0.01), the sample data provide sufficient evidence to conclude that the regression model fits the data better than the model with no independent variables.\n",
    "\n",
    "* **T-test for the intercept and the coefficients**<br>\n",
    "$\\quad$ **`Intercept`** $\\;$ `p-value = P>|t| = 0.00`<br>\n",
    "$\\quad$ **`Coefficient`** $\\;$ `p-value = P>|t| = 0.00`<br>\n",
    "The t-test helps to determine the correlation between the response and the predictor variables. The test is used in linear regression to test the null hypothesis $H_0$ that the slope or the coefficient is equal to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True model: y ~ N(μ, σ^2)\n",
      "True μ = 15 + 2 * x\n",
      "Regression μ = 14.89 + 2.03 * x\n",
      "Mean squared error: 0.52\n",
      "Coefficient of determination: 0.99\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAFNCAYAAADRi2EuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/jklEQVR4nO3de5yXc/7/8ceraWTKYUKoIZNFqFQaZHPMUgij5GyxjpvDol8q7DquSuu4rJVDzsfYQSFRu3xbZGpKDoVV0YRCg2qqmen9++P6fOozn7mumc9nZj6n+Tzvt1s3zXVdn+t6j7mpp/f7db3e5pxDRERERBKvVaoHICIiIpItFLxEREREkkTBS0RERCRJFLxEREREkkTBS0RERCRJFLxEREREkkTBS0SaxMwONrOFqR5HS2Bmn5jZYakeh4gkjoKXiMTEzBab2e+ijzvn3nXOdU3FmKKZ2Q1mVmVmq8yswsz+a2YHpnpcsXLOdXPO/TvV4xCRxFHwEpGMZGatA04955zbAtgOmAG8kIBnm5npz08RiZv+4BCRJjGzw8xsacTXi83s/5nZR2b2s5k9Z2abR5wfZGZzI2ak9ok4N8rM/mdmv5rZp2Z2YsS5c8xsppndaWY/AjfUNy7nXDXwFFBgZh1C99jazB42s2/NrNzMbjGznNC5HDO73cx+MLNFZnapmblwwDOzf5vZX81sJrAG2NXM9jSzaWb2k5ktNLOTI8Z7TOh7+DX0rP8XOr6dmU0Off8/mdm74RAXOatoZm3M7C4zWxb6dZeZtYn8d25mw81seej7ObdxP0ERSSYFLxFJhJOBgUAXYB/gHAAz6w08AlwEbAs8ALwSDhTA/4CDga2BG4EnzaxjxH0PAL4CdgD+Wt8AzGwz4PfAj8DK0OFHgWpgN6A3cBRwfujcBcDRQC9gX6DY57ZnARcCWwIrgGnA08D2wKnAP8xs79C1DwMXOee2BLoD00PHhwNLgQ6h7+MawG/vtmuBvqHx9AT2B66LOL8j3r+nAuA84D4zax/8b0RE0oGCl4gkwj3OuWXOuZ+AV/HCA3ih5QHn3AfOuRrn3GPAOryAgXPuhdDnNjjnngO+wAscYcucc393zlU75yoDnn2ymVUAlXhh6iTnXLWZ7QAcA1zhnFvtnFsO3IkXmMALi3c755Y651YCY33u/ahz7pPQbNpAYLFzbmJoPGXAi8DQ0LVVwN5mtpVzbqVzbk7E8Y7ALs65qlCNnF/wOgO4yTm33Dm3Ai+InhVxvip0vso59xqwCkiLWjsRCabgJSKJ8F3E79cAW4R+vwswPLTMVhEKSDsDnQDM7PcRy5AVeDNF20Xc65sYnv28cy4fbzbpY6BPxLNzgW8j7v8A3mwVoTFE3t/vWZHHdgEOiPpezsCbiQIYghf0lpjZfyKK/McDXwJvmtlXZjYq4PvoBCyJ+HpJ6FjYj6EAGBb571lE0lRQcaqISCJ8A/zVOVdnmdDMdgEeBI4A3nPO1ZjZXMAiLvObGfLlnPvBzC4ESs3s6dCz1wHbRQWWsG+BnSK+3tnvtlHfy3+cc0cGPP9D4AQzywUuBZ4HdnbO/Yq33DjczLoD083sQ+fc21G3WIYX7j4Jfd05dExEMphmvEQkHrlmtnnEr3j/5+1B4GIzOyD0ZmA7MzvWzLYE2uEFmxUAoWLx7k0ZrHNuITAVuNo59y3wJnC7mW1lZq3M7Ddmdmjo8ueBP5lZgZnlAyMbuP1kYA8zO8vMckO/9jOzvcxsMzM7w8y2ds5VAb8AG0Lf1yAz283MDPgZqAmfi/IMcJ2ZdTCz7YC/AE825d+HiKSegpeIxOM1vNqp8K8b4vmwc64Ur+7qXryC9y8JFd475z4FbgfeA74HegAzm2HM44ELzWx7vGL7zYBPQ8+fhFdvBV4ofBP4CCjD+16r8YKR3/fyK15x/ql4M1HfAeOA8IsCZwGLzewX4GK8ZUiA3YG38Gqy3gP+4Zyb4fOIW4DS0HjmA3NCx0Qkg5l/TaeISHYzs6OBfzrndkn1WESk5dCMl4gIYGZ5od5brc2sALge+FeqxyUiLYtmvEREADNrC/wH2BNvGXUK8Cfn3C8pHZiItCgKXiIiIiJJoqVGERERkSRR8BIRERFJkoxooLrddtu5wsLCVA9DREREpEGzZ8/+wTnXwe9cRgSvwsJCSktLUz0MERERkQaZ2ZKgc1pqFBEREUkSBS8RERGRJFHwEhEREUkSBS8RERGRJFHwEhEREUkSBS8RERGRJFHwEhEREUmSjOjjJSIiIhKppKyc8VMXsqyikk75eYwY0JXi3gWpHlaDFLxEREQko5SUlTP6pflUVtUAUF5RyeiX5gNsDF/pGsy01CgiIiIZZfzUhRtDV1hlVQ3jpy4ENgWz8opKHJuCWUlZeQpGW5uCl4iIiGSUZRWV9R5vKJilUsKCl5ltbmazzGyemX1iZjeGjj9qZovMbG7oV69EjUFERERank75efUebyiYpVIiZ7zWAf2dcz2BXsBAM+sbOjfCOdcr9GtuAscgIiIiLcyIAV3Jy82pdSwvN4cRA7oCDQezVEpY8HKeVaEvc0O/XKKeJyIiItmhuHcBYwb3oCA/DwMK8vMYM7jHxuL56GB20QeTWDxuEDNHHwHDhqVo1B5zLnFZyMxygNnAbsB9zrmRZvYocCDejNjbwCjn3Dqfz14IXAjQuXPnPkuWLEnYOEVERKRlKZn9DRsuuZTBH7xS6/iG1rm0qlqf0Geb2WznXJHfuYQW1zvnapxzvYCdgP3NrDswGtgT2A/YBhgZ8NkJzrki51xRhw4dEjlMERERaSnWroVBgygu6lwndAE83mdQSt9uTMpbjc65CmAGMNA5921oGXIdMBHYPxljEBERkRZs5UrYZx/Iy4MpU3wvuWDwddxw2HkpfbsxYQ1UzawDUOWcqzCzPOBIYJyZdXTOfWtmBhQDHydqDCIiItLCff21F7h+/jnwksFnjmdOwV4bv07l242J7FzfEXgsVOfVCnjeOTfZzKaHQpkBc4GLEzgGERERaYnmzYNevYLPt2rFKVdO5IPW29Y5lcq3GxMWvJxzHwG9fY73T9QzRUREpIWbNg2OOir4fJcu8P77sP32nFZWzkcRWwtB7bYTqaC9GkVERCTtRO+1eE/VfPr85YrgDxx6qFfb1a7dxkPh9hLptGejgpeIiIiklY2bYK+vZtj7L3D1O48HX3zWWfDII9C6dqSJDm53ntIrLTbJVvASERGRtHL7659y3eS7OWPuG8EXXXcd3HQTmNU5tTG4hZYYw5tkAykPXwpeIiIikhZeee9Ltj/7dN794sPAa8quHculW+znzWSNm+G7dFjfJtkKXiIiIpLdfvqJX4r6cvyiLwIv+cOQvzCnRz/WVW+gMtQOImgmK503yVbwEhERkdRYvBi6d4fVq9kq4JITzrqdeZ28vRc3N2KayeqUn0e5T8hq0Ztki4iIiPiaM8erzerSBVavrnN6XU5rDrnwQQpHTmZep64U5OcxpE8BK9dU+d4ueiYrepNsgNwcY/W6arqMmkK/sdNTtm2QZrxEREQkOaZOhYEDA08vat+RIWf+jZ/abr3xWEGoBUR4SdFP9ExWdBuJ/La5rFpbTUWlF9xSWWyvGS8RERFJrIkTvRmuoNDVvz+vzvyCYy55pFboCjc79SuWj74mWnHvAmaO6s+iscfSdrPWVG1wtc6HlyiTTcFLREREmp9zm9o9/OEP/tecey5UVcHbb3Pcb3djzOAeFOTnYXgzXWMG96C4d0G9RfHha+qTTsX2WmoUERGR5lNTw+KTzqSw5Nnga66/3vsV1YOruHeBb4gKKpYvyM+LaakwnYrtNeMlIiIiTVdZ6e2h2Lp1cOh66CFvJuyGG3wbnwbxK5aPZ8/Fpn6+OWnGS0RERBrvhx+gXz/4/PPAS8456Qa+6HMwM8/r36hHNHXPxXTas1HBS0REROL31Vew996wbl3gJYPOvouPd9wNAGtiPVXQMmSyPt9cFLxEREQkdqWlsN9+gafX5rbhyD/cxzf5O9Y63tR6quhNr1M1Y9VUCl4iIiLSsNdeg2OPDTz91bY7seD511jffht+eGk+RLR/aGo9VTpveh0vBS8REREJ9vDDcP75gaffKezNBYOvY11uGwpmrWDmqB5A89ZTpfOm1/FS8BIREZHawm8e3nRT4CVP9xzAdUcNY0OrTW8LhvtiNXc9VTr14WoqBS8RERHxVFd7s1uPPRZ8zc030y/nQMp/XlvnVKL6YqVTH66mUh8vERGRbLdmDfTvD7m5waFr4kRvJuy66xgxcM+k9sVKpz5cTaUZLxERkWy1YgX07eu1hghw5bljOfSyM2stHSa7L1Y69eFqKnPONXxVihUVFbnS0tJUD0NERKRl+PJL2Gsvb2kxwLHn3M0nO/wG8GaXYtkTUTxmNts5V+R7TsFLREQkS8yaBQccEHx+iy0YMuwBZtvWdU7l5+XSrk3rjJ9xSob6gpdqvERERFq6V1/19kYMCl3dusGPP8KvvzLHJ3QBVFRWUV5RicPro3Xlc3MpHDWFfmOnU1JWnrixtzAKXiIiIi3VAw94gev44/3PH320t7n1xx/DNtsAsb8pGF4vCzczVfiKjYKXiIhIS+IcXHutF7guvtj/mosvhpoarxv95pvXOuX3BmFDws1Mo5WUldNv7HS6aGZsI73VKCIi0hJUV8O558KTTwZfM2YMjBpV722i3yCMtRI8uplpS9rmpzkpeImIiGSyVavgmGPg3XeDr3n8cTjrrJhvGdl5vt/Y6b7NS6NFL1G2pG1+mpOWGkVERDLR999DYSFsuWVw6Jo2zVt6jCN0RYtl6dGvmWlL2uanOSl4iYiIZJLPP4dWrWDHHWHJEv9r5s71AtfvftfkxxX3LmBInwIs6nj464L8PN8eX0FF+pm4zU9zUvASERHJBO+/7xXMd+3qhapo+fleEHMOevZs1kfPWLCiTq2XwwtdM0f19106bEnb/DQnBS8REZF09vLLXuA68ED/8z17wk8/wcqV0LlzQobQmGXD4t4FjBncg4L8PIzgmbFso+J6ERGRGJWUlSdvv8B//AMuuST4/HHHwfPP12kHkQid8vN8C+wbWjaMLNIXj2a8REREYhBujxDZvb3ZG4c657V7MAsOXZdc4vXgeuWVpIQu0LJhc1LwEhERiUF97RGarKoKTjvNK5ofN87/mnHjvGB2773edUmkZcPmo6VGERGRGCSkPcKvv8KAAfDee8HXPP20F8pSTMuGzUPBS0REJAaNrXPy9d13UFQE5fUsU06fDocfHv+9Ja1pqVFERCQGzVLntGCBV7/VsWNg6Jr+XKjpaSNCl/ZGTH+a8RIREYlB9B6Gcb3VOHMmHHRQ4Omf8rbimHPu4buttiNvfjVjdi+Pe1lPeyNmBnN+TdjSTFFRkSstLU31MEREROLz0kswZEjg6fk7/IbTT7uVX9u0q3U83Jg0HkF7KjbmXtI0ZjbbOVfkd04zXiIiIs3tnnvgT38KPP3GHgdy+XFXs751ru/5xhTsB21kne17I6YbBS8REZGQJjVIdQ5GjIDbbw+85PnfDmbkQefgrP4S63gL9kvKyjGos61PY+4liaXiehEREZrQIHX9ejj5ZK+3VlDo+tvfwDk2u/ceNt/Mf5YrrDGNScdPXegbugzU5DTNKHiJiIjQiAapv/wCfftCmzbwwgv+1zz7rDcTNnw44N+I9My+nZvcmDRoOdGhwvp0o6VGERFp0WJdPoy5QeqyZdCnj9eLK8i//w2HHup7KhGNSIN6jBVomTHtKHiJiEjGiLcGK54WCw02SP30U+jWrf4Bfvxxw9c08nupz4gBXWt9n6C9FNOV2kmIiEhGiA5R4IWL+pbmglostG+bS9vNWtcKPYDv/ScUruHg84NbQrDDDjB7NhTEHpoa873Ecs/mCnLSNPW1k1DwEhGRjBBrn6rIABLr33Dh0AObGqSe8c0sbnn6puAP7b8/TJsGW20Vz7cBqOdWS6c+XiIikvFiqcHym0mKRWVVDcOfn8cG57jy4ylcPuX+wGundO3HrJvv4cah+8b1jKAxx3JcWg4FLxERSblYlsli2aTa783EmDjHNW8/yHmlLwdeMmG/E7n18D+AGTlzvuPGofE/Jvx9Bs3EqedWy5ew4GVmmwPvAG1Cz5nknLvezLoAzwLbArOBs5xz6xM1DhERSW+xFsDHUkDe0IxRQX4eq9dVU1FZBUBuTRV/f+U2Bn7+XuBnbjziAiYWnVDrWE0jynQamo1TMXx2SOSM1zqgv3NulZnlAv9nZq8DVwF3OueeNbN/AucBwXO6IiLSotXXPysyeMWySXXQrBhsqp8qKSvnr0+/z8NPjmKf774MHNewE0Yxda+DfUNWjllc32N43EGhq0DF8FkjYcHLeVX7q0Jf5oZ+OaA/cHro+GPADSh4iYhkrXjqnRrqgTViQFeueG5u8P3Kyyk+qhfFP/wQeI+hp4/lw527U5Cfx2l7duDJ97+uc81pB+wMxPcmYdD3aaCC+iyS0M71ZpZjZnOB5cA04H9AhXOuOnTJUkDxXkQkiwXVNTWm3qm4dwHt29bdkmePFYtZNG4Q7LQTBISuI86/n8KRk/lw5+4bl/2Kdtmmzl+UrYCiXbaJe4uh5vw+JXMlNHg552qcc72AnYD9gT1j/ayZXWhmpWZWumLFikQNUUREUmzEgK7k5ebUOtaUeqfrj+u28X59v/6IxeMG8eYjl/pf3KkTLFtGyZylrP3NHnW27Rk/dSEboj6yAW/ZMN4thpr7+5TMlJS3Gp1zFWY2AzgQyDez1qFZr50A3/81cM5NACaA18crGeMUERF/iWzOGUvtVrz3K5j6MvuNviT4ogMPhKlTYcstvc909N/TsDFtH4LONff3KZkpkW81dgCqQqErDzgSGAfMAE7Ce7PxbCD43V0REUm5eLbdaazm2L+wpKycpdfezKWvP8B+QRedcgo88QTk1l2O9NNQC4uG2ltES8Q+jZJZErnU2BGYYWYfAR8C05xzk4GRwFVm9iVeS4mHEzgGERFponiX1JJuwwa+OvVcivfdiUtff8D/mquvhg0b4NlnYw5dUP/yoJYOpTES+VbjR0Bvn+Nf4dV7iYhIBkhkl/WgJcyYljbXrYOTT4ZXXmHXgPvffvxlDH/5nkaPL5blQS0dSjy0V6OIiNQraF/B/LxczGDlmqqNX99wfLeYg0fQRtFD+hTw4uzy4A2kKyrgsMNg3rzAe1904jVM3eO3GLBo7LExjUekudS3V2NC32oUEZHM57ekltvK+HVd9cbQBVBRWcWIF+YFtlOIFrSE+cwH3/gef+y5d2GbbaB9+8DQNfiM8RSOnMzUPX4LqFWDpB/t1SgiIvXyW25bs7526Aqr2uDqdJwPErRUGd0pfs/li3hj4mX13uuoCx/g8/a1n5mbY6q3krSj4CUiIg2Kfhuvy6gpgdfGWvsV9MZgjhk1znHgknk88+y1wTfo3BlmzaLfxE9879Nus9aqt5K0o+AlIpKlYilgD7qmvj0RY13eC9r0+qZVcxl61+jgDx58MLz2GmyxBQDLKvxrgH+urDsjJ5JqCl4iIlkolt5c9V0zYkBXRkyaR1VN7WXB3FaxL+9FL2FePbeEP059KPgDZ5wBEyfWaQfRUK8tkXSi4noRkSwUS2+u+q4p7l3A+JN61toXMT8vl/FDe8a1vFfcsyMzv57EonGDAkPXfX2HUnj1q/Ta/feUfLy8znn105JMohkvEZEsFEtvroauaVIX9rVrYcgQb8kwwJ+PGsYTvY/Z+HX4rcnws8O0FY9kEgUvEZEsFMvyXEKW8Fau9Gq0Pvkk8JLzB/+Zt3Y/wPdc0FuT2opHMoWWGkVEslAsy3PNuoS3ZAlstZXXhysgdBWfdTuFIycHhq6w5uiYL5IqmvESEclC0ctz+W1zcQ6ufG4u46curLVU16QlvLlzoXed3eM2ycmBzz6jy8OfE+s+KkEzbjFtMySSYgpeIiJZKrw819Abjo0KL9OmwVFHBZ/v0gXefx+23x6ATvnfBLaniBT01mQsb2mKpAMtNYqIZLkbXvmkwTccY/b442AWHLoOOwxWrYKvvtoYusB/WTNafW9NxvKWpkg60IyXiEgLE8+SW0lZORUBjUZjrqVyDm69Fa67Lvias8+Ghx6C1v5/7YTHN/z5eXW2DAIoyM9j5qj+gbeP5S1NkXSg4CUi0oLEu+RW34xQg28v1tTAsGEwYULwNX/+M9x4ozcLFjDeyJB42gE78+Ls8jrd7Bsq6FcTVckUWmoUEWlB4l1yq29GKDDsVFbCwIHe7FVQ6JowwZsJu+mmekPX6JfmU15RicMLiS/OLmdInwIK8vMwvJmuMYN7NFinpSaqkik04yUi0oLEu+QWNFPUvm1u3bDz449w0EGwYEHwACZPhmOPjWmsQSFxxoIV9S4r+lETVckUCl4iIi1IvEtuQRtVX39ct00XLVoE3bp5M11BZs2C/faLa6zNXZelJqqSCbTUKCLSgsS75Fbcu4Axg3v4L+3NmeMtE+66q3/oatMGvvzSW1KMM3RBcBhUXZa0ZApeIiItSL1BKkbtvl4EF10Effr4X7D77rB8ubff4m9+0+ixqi5LspGWGkVEWph4ltwi34Ls9v3/+OPLk+i/cCa4DXWu/XC3ffn98dewTYf2jFi6nuIOTR8nqC5Lsos5n34p6aaoqMiVlpamehgiIhnNr7/X+DcWsPP8Wfzx/UkcumiO7+eWnHAKx3Q7i9URdfC5rYwtNm9NxZoqBSaRKGY22zlX5HdOM14iIlkgur/XspWrmTZmAn+f+Rz7LqvbauI/Xfbl0Il3wCGHcPq4GayOKniv2uBYucZrvKrteURip+AlIpIFwq0bWtdUc/xn/+Hi919kjx+/rnVNjbXita79+Gffk6jo2p2Zhx4KxPaWYbhXmIKXSP0UvEREEiyeLXwS5acVKzn7o2lcMOsldvplRa1z63NaM6n773jggMEsad+JvNwcxgzounHcsRakaHsekYYpeImIJFC8W/g0u5Ur4d57+e8Dt9N+9c+1Tq3aLI+X+x7HttdezX1zfmZZRSUFoWAI1Onv1RC1gRBpmIKXiEgC1beFTyKCV3iWquabb7h8/hSGzn6N3DWraR9xzQ9tt2Zin+OZtP9xjD7jtwzsXcDAo2rfp9/Y6YGhKz8vl9Xrq6mq2TQXpjYQIrFR8BIRSaDm7s5en5Kych548HUunfkCgz+ZTpua6lrn13TciX8WnciDux/KNh3aM7qeJc+g8Rkw9/qj0mL5VCQTKXiJiCRQ0BY+rczoMmpK84WW2bPZ6pyrmPLxu7SKqsr6aodCdv3bzbQ95RSuys3lqnpu01BdV3g5UdvziDSOOteLiCSQX3d2gBrncGyq+SopK4//5s7BjBlw1FFQVET/j9+pFbpmd9qT84b8md+dfQ+ceSbk5tZ7u3A9ml9QBC0nijQHzXiJiCRQdHf2VmbURDWujrvma8MGePllGDvW25w6yoxd+3B/36HM2qkbmFEQY9G7Xz1aWIGWE0WahYKXiEiCRS7LdRk1xfeamGq+1q+Hp5+GceNgwYLa51q1YumRx3Fp56OYu80uGw/HM0tVX13XzFH9fc+p1kskPlpqFBFJoqCWC/W2Yli9Gu6+G3bbDc49t3boatMGLr4YPv+cnd4o4Zw/ntDoDbLjHVvk0mSTl01FsoRmvEREkmjEgK51+mMFzkr99BPcey/ccw/8+GPtc1tuCcOGwZ/+BB07bjzclKL3uMZG8ltliLQECl4iIkkUXfPluzy3dCnccQdMmODNdkXafnu44gr44x8hPz/5Y4uQzFYZIi2FgpeISJIFzkotXAi33QZPPAFVVbXPFRbCiBHeUmNe4jrExzNjFtQqQx3sRYKpxktEJNVKS+Gkk2CvveCRR2qHrh494Kmn4IsvvKXFBIauePm1ylDLCZH6acZLRCQVnIPp02HMGHj77brn+/WD0aPhmGPALPnji0G8S5MiouAlIpJcGzZASYnXg+vDD+ueP/ZYGDUKDjoo6UNrDHWwF4mPgpeIZIyM7hm1fj08+aRXw7VwYe1zrVrBqafCyJGwzz6pGZ+IJIWCl4hkhHDPqHD7gnDPKCC9w9eqVfDgg3D77VAe1d+qTRs47zwYPhx23TU14xORpFLwEpGMkHE9o378Ef7+d+/XTz/VPrfVVnDJJV4Prh12SM34RCQlFLxEJCNkTM+ob77Z1INrzZra53bYAa680us0v/XWqRmfiKSUgpeIZIS07xm1YIG3h+KTT0J1de1zu+7q9eA6++y0agchIsmnPl4ikhHStmfUrFkweDDsvTc8+mjt0LXPPt6m1gsXerNcCl0iWU8zXiKSERLRMyqetyRrXbv15tyWv5x+kx7yenFFO/hgrwfXwIFp24NLRFLDnHOpHkODioqKXGlpaaqHISItSPRbkuDNoA3pU8CMBStqhTGA0S/NZ9269Qz4/D3++MEk9vnuy7o3HTTI68HVr1/gMzO2HYaIxMzMZjvninzPKXiJSDbqN3a6b82YAZF/Kubl5rClVXPYh9O4aNaL/Oan2i0hqq0Vb/Xsz2ajR9L/5N8FPi8o6I0Z3EPhS6SFqS94aalRRLJC9GyTX+iC2qGr3bo1nDbrDc7/sIQdV9VuCbG29WY8t8+RPLjfiSzN35G8+dWM2b08MERlXDsMEUkIBS8RafH8mq9Gz2xF2mbNz5w9+1XOnjOZ/LWrap37pU07Ht/3WCb2OZ4f2+VvPN5QiMqYdhgiklAKXiLS4vnNNjnqLisW/LKc82f9i1PnvUle9bpa1y9v156H9ivm6V5Hs6pNW9/n1Bei0r4dhogkRcKCl5ntDDwO7ID3Z9sE59zdZnYDcAGwInTpNc651xI1DhGRoEDkgIL8PNp+uZAry0oYMG86OTW1e3Atzu/IAwcM4aXu/VnXerN6nxMdoiKXN/Pb5pLbyqjasCnqpUU7DBFJqsDgZWavAcOcc4sbee9qYLhzbo6ZbQnMNrNpoXN3Ouf+1sj7iojEJWi26chfFvPggmnw8st1zn2y/a48d8TpPLnTfmxolVPnfLToEBW9vLlyTRW5OUZ+Xi4/V1bprUaRLFXfjNdE4E0zewy4zTlXFc+NnXPfAt+Gfv+rmX0G6E8YEUm6EQO6bgpBznHIojlcMutFDljyUd2LDz0URo+m21FHcZMZbwe8/di+bS5tN2sd2BrCb3mzqsbRrk1r5l5/VLN/jyKSGQKDl3PuBTN7HfgzUGpmTwAbIs7fEetDzKwQ6A18APQDLjWz3wOleLNiKxs3fBHJZMnqa1XcuwBqaph9x0OcMv1pun//v7oXHX+814PrwANrHa4V2kLycnO4/rhu9Y5VxfQi4qehGq/1wGqgDbAlEcErVma2BfAicIVz7hczux+4Ga+84mbgduAPPp+7ELgQoHPnzvE+VkTSXElZOSMmzaOqxqt5Kq+oZMSkeQBxha8Gw9u6dfD44xTfdhvFX0Y1Pc3JgTPOgKuvhm7dfO/f2I75KqYXET+BDVTNbCBwB/AKcJNzbk3cNzfLBSYDU/1myEIzYZOdc93ru48aqIq0PL1vepOVa+pWMLRvm0vZX7yluIZCVb1NSXfbCh54AO64A779tvZD8vLg/PNh+HDYZZeEfH9qmCqSvRrbQPVaYKhz7pNGPtSAh4HPIkOXmXUM1X8BnAh83Jj7i0hm8wtdkcf9em+Nfmk+UHsWKrqOKu/nn6gYPhLKpkBFRe2b5+fDpZfC5ZdDhw7N9834SMTekiKS+eqr8Tq4iffuB5wFzDezuaFj1wCnmVkvvKXGxcBFTXyOiLRAsXR6j6yX2unn7zl/1r845aNpdXpw0bEjXHUVXHQRbLllwsceVty7QEFLRGpJWB8v59z/4fUnjKaeXSJCfl4uFZV1Z73y83KB2IrTO+Xn0e6Lz7jogxc54dP/0NpFlaHuthtlp13IlZv1YMkPNXS670PNOolISrVK9QBEJDvdcHw3clvV/n+z3FbGDcd7Re5BRegbj7/3HpNeH8ubj1zKkE9m1ApdFXt2h+efp+SZtznd9WDxqhoc3nLllc/NpXDUFPqNnU5JWbnvM0REEkXBS0RSorh3AeOH9qQgPw/D6yA/fmjPjbNRIwZ0JS+3duPSvNatGL/FMq/X1m9/S8d3ptU6P3vXXsy87ynyP/0Ihg5l/Ftf+m4VBJtqxhS+RCSZtFejiKRMfTVQkcXp3/+0itO/mcXwOf9i688/9bm4GEaOpE/fvrUON9Qzq6GNrUVEmpuCl4ikreK9tqV41hfw0Hj4X1TT09atvR5cI0fCXnv5fj6ol1YkNTQVkWRS8BJpgZLVET5hfvkF/vlPuPNO+O672ufatoULLvDeUmygubJf1/loamgqIsmk4CXSwsTS/yptLV8Od98N990HP/9c+1z79nDZZd6v7baL6XaRy5XlFZUYm2q8oO7G1iIiiabgJdLCxNL/Ku3MmAF//zu8/jqsXVv7XKdOXof5Cy5oVA+uyDqyjJ8JFJGMp+Al0sJk1ObMzzwDp5/uf26PPSg79QIu3LAnK5Yb/PUd2rfNbXBz6vqooamIpJqCl0gLkxGbM992m1cU72NBwR7cvd8QZvU8mJVrNxDZEnXlmqpGbaQtIpIu1MdLpIXx7X+VDrVMGzZ4+ySaBYaumwcOY+AZt/N61378GBW6wqpqHOOnLkzsWEVEEkQzXiItTNptzrxuHQwdCq++GnzNvffS79e9Gmz9EJaWy6YiIjFQ8BJpgdKilqmiwusw/9FHwdeUlMAJJwCwbNSUmG+dVsumIiJxUPASkeb19dewzz5120FEGHzmeL7vti8jOnelOHQslmanALk5lvplUxGRRlLwEpHmMW8e9OoVeNqZcczFE/hsq47egaj+Yn7NTnNzjNatjMoqr9qrqW81ioikmoKXiDTNW2/BkUcGn99lF/jgAw6a+EmdGa3I/mJpV5smIpIACl4i0jhPPglnnRV8/uCD4bXXYIstAFhWUep7WWShfFrUpomIJJDaSYhIHSVl5fQbO50uo6bQb+x0SsrKN50cM8ZrCREQul7qdjjdr5lMyd3PbAxdEFwQr0J5EckmCl4iUkt4r8fyikoc3l6P17w4j0UnneUFrmuu8f3c3w88hcKrX+WqQcNZVUOdXltp219MRCSJtNQokiGStc9g5F6PbarX889//ZXDv5odeP01Ay7h6V5H1zke3WtLNVwiIgpeIhkhPAsVDkTlUW8ExvL5yMBz+J4dmLFghW8AWlZRydaVv/LCUyPZ48evg2/6yiuU7LQvzz0/D5yrc9pvCVE1XCKS7RS8RDJA5CxUWOQbgfXxC21Pvr8pUNUKcdtU88mdQ2m7vp5+Wh98APvvv/G+NT6hS0uIIiL+VOMlkgGCtsiJZescv9AWbdeln1O8705QWOgbuqpa5TCt5F1vZmv//eu9b44ZYwb30MyWiIgPzXiJZICgru6xvBFYXzg7eNEcnnj+L4HnF+d35JJL7uWCIQfUCVJB993gnEKXiEgABS+RDODX1T3W5Ty/0DZk/tvc/tqdwR86/HB49VUK27UjaAfFpoRBEZFspaVGkQxQ3LuAMYN7UJCfhwEF+XkxL+dtbOPgHJf+91kWjxsUHLrOOQeqqmD6dGjXLrb7RlBtl4hI/cz5FMamm6KiIlda6t/1WkQaUFPD4qFnUfivZwIvWXDhlez5z9u9Pl1xSFaLCxGRTGJms51zRX7ntNQo0gKVlJVz9+SPuPGRazlkcRmFAdeNGngZ7x5azIgBXdkzztAFag8hIhIvBS+RDBLLDNOUGR+zz9CjmfHj0sD7XHTqjUzdpY/3RT09wTSjJSLSvLTUKJIhovtxARhwRt/O3FLcAxYtgr33hrVrg2/y4Yf0e+sX36L4gvw8Zo7qX+/z8nJz1CpCRKQB9S01qrheJEP49c1ywLyS6V5t1q67+oaudTm5HHzRQ3QZORmKimLuCVZf01YREWkcLTWKxClVy2/Rweiw/5Xy6KQbAq//3zYFnHTGbaxsuzXgzWhB7G0gmtK0VURE/Cl4icShqXsmNkU4MJ08701ue+OewOveKezNRYOvpTJ3843HIts8xNoTTH26RESan4KXSByasmdikzjHw4sms+eE4Kank3oN4Oojh7Ghlddby/CWIguiZuXC/2xo1q4pTVtFRMSfgpdIHJK+/FZTA+efD48+yp4Bl9xz6FlMPPwMVlZW1zoeDl2RBfNhsbSBiDWgiYhI7BS8ROKQtOW3NWvguOO8DvIBrj76cmYecgIjBnSl4rm5vtc0NRCqT5eISPPSW40icUjUNjklZeX0GzudfS9/mm+228nbricodL3xBjjHba/dzcxR/SnuXRAY/FqZUVJW3qSxiYhI89GMl0gc6lt+a+zbjiVl5fzjoanM+OdFbLahOvC6cy+9nxP+cJzvPf3qsQBqnEta8b+IiDRMDVRFmkGjm41++CHsv3/g6dW5mzPgvPtYuvUODd6zpKyc4c/Po8bnv+mgWi8REWl+aqAqkmBxNxudMsVrehoQuj7ftjM9L3+GbldN2hi6Grpnce8CNgT8j5R6b4mIpActNYo0g5jfdpwwAS66KPA+/+7Sh4tOvIZ1uW3ifhao95aISLrTjJdIMwgKNvltc8E5uO46b4YrIHQ91Wsgu454mXNOvrHe0FXfsyBxxf8iItI8NOMl0gxGDOjKiEnzqKrZtNSXs6GGv7xwB1w/IPBztx3ye/5x4Mne9WaYc2ydl8vq9dW17hXWUIhS7y0RkfSm4nqRZtLrxjepqKwib/1aJk66gb7ffBx47VXHXslL3Y+odcyARWOPBTbtB1leUUmOGTXO1elALyIi6am+4nrNeEnWSdQm161/WMG7Twxn55+/D7zmT+fdxsvb7e17LnIJUY1LRURaJgUvySoJ2eT6iy+ga1dK65s9LiuDXr04vKyc116YR9WG2tfm5pjqsEREsoCK6yWrxN32oT7vv+8VzO+xh1dAH+WXNu044tJHKZmzFHr1ArxwN35oT/Lzcjde175tLuNP6qkZLhGRLKAZL8kqzbLJ9csvQ3Fx4Okvd+zCkFPGsMWOHXyXMSOXEcPLnlc+N5fxUxeqhktEpIVT8JKs0qQ+V/ffD8OGBZ8/9liYNIndNt+ceTGMpbHLnomqURMRkcTTUqNklbj7XDkHo0Z5S4pBoWvYMKipgcmTYfPNYx5LY5Y9w2GtvKISx6awpo2wRUQyg4KXZJXi3gWMGdyDgvw8DG8PQ9+9D6uq4PTToVUrGDfO/2Zjx3rB7L77vOvi1Jhlz2atURMRkaRL2FKjme0MPA7sADhggnPubjPbBngOKAQWAyc751Ymahwi0ept1bBqFQwcCDNnBt/gySfhjDOaPI7GLHs2S42aiIikTCJnvKqB4c65vYG+wCVmtjcwCnjbObc78Hboa5HU+u476NwZttwyOHS99ZY3w9UMoQsat71PUCjTXowiIpkhYcHLOfetc25O6Pe/Ap8BBcAJwGOhyx4DihM1BpEGLVzo1W917AjffON/zUcfeYHriCP8zzdSzMueEbQXo4hIZkvKW41mVgj0Bj4AdnDOfRs69R3eUqRIcs2cCQcdFHx+2229pqc775zQYcTboV57MYqIZLaEBy8z2wJ4EbjCOfeLmW0855xzZubb7tvMLgQuBOjcuXOihyktQExtFv71Lxg8OPgmvXvD9OmQn5/QsTaFthMSEclcCX2r0cxy8ULXU865l0KHvzezjqHzHYHlfp91zk1wzhU554o6dOiQyGFKC9Bgm4V77/WWFINC1wknwNq1MGdOWocuERHJbAkLXuZNbT0MfOacuyPi1CvA2aHfnw28nKgxSPbwbbOwvprVV1zlBa7LLvP/4OWXez24SkqgTZvED1RERLJaIpca+wFnAfPNbG7o2DXAWOB5MzsPWAKcnMAxSJaIbKfQuqaauybfzqAF7wZ/4G9/g+HDkzAyERGRTRIWvJxz/wdYwOnmfT1Msl6n/Dwqvv+RJ5/7M72/raeZ6LPPwimnJG9gIiIiEbRXo2S+b7/lrdtPJ++H74OvmTEDDjssaUMSERHxo+Almeuzz2DvvQEIbB86fz507560IYmIiNRHezVK5nn3Xa9gPhS6oq1om8+hlz9ByZylCl0iIpJWNOMlKRFTz61oL7wAJwe/izFvx90589Rb+LVNO8B701H9rkREJJ0oeEnShXtuhds/hHtuAf5B6a674MorA+/3+h6/5fLjR1CVk1vreOSbjo0KeiIiIs1MS42SdL49t6pqGD814m1E5+CqUA+uoNB15ZVQU8Mt595cJ3TBpo2jG2yuKiIikiQKXpJ0kTNRkcorKtljRAkzehwCrVrBnXf63+DOO71gdscd0KpVgxtHxxT0REREkkBLjZJ0nfLzKI8KX1usW8NTz15Lz+++CP7gCy/ASSfVOhReQqysqiHHjBrnKIhaSgwKekHHRUREEkXBS5JuxICuG2u8dvj1ByY/egUd1lQEf+Cdd+Dgg+scjq4Vq3Fu40xXZP2WX9ALHxcREUkmLTVK0hX3LuDenpuxeNwgPvjHOYGh68jz/uEtKfqELoh9CbGhpUgREZFk0YyXNLt63yD8z3/gsMMC94z6fottOO73d7J8y20paGBGKtYlxPCz9VajiIikmoKXNKugVhEFU19hv9HDAj83t2BPzhh6E6vbtAVim5GKZwmxuHeBgpaIiKSclhqlWUUv/50/6yU+u+Xo4NB1yimwfj2LX32L/B22xYCC/DzGDO7RYFDSEqKIiGQazXhJs1pWUYm5Dfzl7Qc5d/argdc9sP9gHi8exoiBe1Kcm9uoGSktIYqISKYx51yqx9CgoqIiV1pamuphSEPWreOdfftzyKf/DbzklgEX81CvQRu/zsvNiWl2S0REJFOY2WznXJHfOc14SdNVVED//lBWxiEBl3wwfgJXVf+mTk1W+C1EBS8REckGqvGSxlu6FLbdFtq3h7Iy30suuvgeSuYs5YD/d4EamYqISNbTjJfE76OPoGfP+q9ZsAC6duWBiENqZCoiItlOM14ZoKSsnH5jp9Nl1BT6jZ2eus2dp0/3Nq0OCF3fbb09r785x2t62rXum4V6C1FERLKdZrzSXFBfLKBWXVS9TUvjuMbX00/DGWcEnv6wYG/OGXoDq9u0Je+d71m3XbnvffUWooiIZDu91Zjm+o2d7rs8V5Cfx8xR/YG64Qzqvi0YyzWRSsrKWXbNDQx746HAsb3Zsz/Djryc6pza+T1ybCIiItlGbzVmsKDC8/KKSrqMmkKn/DzWrK8O3LMwcpapoWsA2LCBr079A8UvPBY8qFGj4NZbuWj0a/jFdhXLi4iI+FPwSnNBBekADgLPQe0A1OAbhWvXwkknwZQp7Bpwv/En/IkRJXc1OLZwsXyjlzZFRERaKBXXpzm/gvRYRb4tuHVeru81O7EW9tkH8vJgyhTfay4YfB2FIyfzjz2PbHBs4WL58NJmeUXlxoA4+qX5qXsxQEREJA1oxivNRRekx1qRF/22oFnt851+Wc4bj1zGVutWB95j8JnjmVOw16bPRLV9qK9Yvt/Y6bEtbYqIiGQRBa80FbRMF1Rsn5+XS7s2rX2X9UrKylm5pgqAvZZ/xesTLw9+cKtWTHvp31xeurpOIb5f24egPRbVLFVERKQuBa80VF8LiREDuvq+nXjD8d0C304c8cI8DlpUxpPP/zn4oV26wPvvw/bbcyQwpnPT6rPULFVERKQu1XiloYbeQBwzuAcF+XkYXuuG+jaZ/vDmu/lizLGBoeuHPn1h1Sr46ivYfvtm+x7ULFVERKQuzXiloYaW6YKW9zZyDsaMgWuv5a8Bl7zY7XBaT5zICfvtUudcrE1b66NmqSIiInUpeKWhRi/T1dTAJZfAAw8EXnLPgadwx8FnghmLfUIXxNHzqwENBkQREZEso+CVhoLquAKX6Sor4cQTYerUwHuOHnApz/QauPHr9m3920uACuNFREQSRcErDcW8TPfTT3DQQfDZZ4H3eu+uR/n98g5U1WxqRJGbY1x/XLfAz6gwXkREJDEUvNJMdBuJO0/pVTdwLV4M3bvD6uAeXHzwAey/PwcC4+PsIB/3jJuIiIjERMErjTRY1D5nDvTpE3yD3Fz49FPYbbdah+OttVJhvIiISGIoeKWRoKL2d+55guJHRwd/cLfd4L//hQ4dmm0sKowXERFpfgpeMfLrJA/NOysUXbw+9KNpjH/97uAP9O8Pr74Kbds2+pkiIiKSPApeMfBbAhzxwjwwNhatN6bXVbRO+XmUr1zDZf99luH/91TwheeeCxMmQGv9+ERERDKJ/uaOgd8SYNWGuttVN2kT6JoanvrgQQpLng2+5vrrvV/RO16LiIhIRlDwikE8/auWVVQGbnDtq7ISTjgBpk2jMOimDz0E550X77BFREQkzWivxhjE078qv20uo1+aT3lFJY5NS5AlZeW1L/zhB+ja1avPmjbN917Dz7mVkjlLFbpERERaCAWvGIwY0JVYFvfycnNwjsDtdgBvM+rNN/feQPz8c9/7DDr7LgpHTubFHfbxD20iIiKSkRS8YlDcu4C6FV2bGFCQn8eYwT34ubLK95ptP/vIq836zW9g3bo65ytbt+Hgix6icORkPt5xUx+uWqFNREREMlrW13jFWo9VELCNTkF+HjNH9d/49fipC2tdd9j/PuTRSTcGD6BrV44ZOoZPqzYLvER7JIqIiLQMWT3jFW4T0WA9Ft5yY15uTq1jftvohK87ed6bLB43KDh0HXUUrFkDCxbwWT2hC7RHooiISEuR1TNeQZ3i/VpCxLSNjnMUl0yg+Jabgh96wQVw//2QsynEBW1KDcF7JMb15qSIiIikhawOXkFLeEHHA7fRqa6G88+Hxx4LfNbtB5/JS0efw4iBe1KcU3vmzG9TaoD2bXO5/rhudZ7Z4J6OIiIikpayOngFzTSFl/YanFVaswYGDYIZMwKfMfLYK3iu+++8L35e6xuQ4t2UOp6ZOhEREUkf5lx97+ulh6KiIldaWtrs942eOQJvaW/M4B4Adc4Z4IBurdfx/GNX0W7pksB7/37ojbyzax/fc9EF+fHqMmqK71uWBiwae2yj7ysiIiJNZ2aznXNFfueyesarvpmmfmOn15lV6rxyGW899EdyN9T43Q6AY8++i08i2kH4aexbiuEZuKCorCJ8ERGR9JbVwQuC67Yiw1GvZQspeWJ48E222ALmz6fLPz+pt99XWGMCkt/sXKSgInwRERFJH1kfvIJ0ys9jz9J/8/CLNwdes3C7znRdWAbbbBP6zFeBbyeGNTYg+dV1hRXorUYREZGMkLDgZWaPAIOA5c657qFjNwAXACtCl13jnHstUWNorLnXjmPmraMCz8/YtQ8Xn3gteVu2Y24odIH/24m5OUa7zVrzc2VVk9o+BC1PGjSpXkxERESSJ5EzXo8C9wKPRx2/0zn3twQ+t3Gcg+uug1tvpVfAJU/2Opo/H/VHnHl9Zzesr6akrHxjkIr37cR4NPQGpoiIiKS/hAUv59w7ZlaYqPs3q/Hj4eqrA0+PO/RsHjhwKBuiCriqalydFg6Bvb6ayG82TXVdIiIimSUVNV6XmtnvgVJguHNuZQrGsEl5eWDouvLYq/hX99AyXkDVfHlFZa1Zr0RJ5GyaiIiIJEdC+3iFZrwmR9R47QD8gBdjbgY6Ouf+EPDZC4ELATp37txnyZLgnllN8uOPsN12tQ6dccotzCzsVetYjhk1Af+uwr2/wiFI2/mIiIhkr7Tp4+Wc+z78ezN7EJhcz7UTgAngNVBN1JhKvl7L1FNu4PiyN7mn36l8tv2uda7Jy81hSJ8CXpxd7vtmYbhrPMCNr37CyjVVG89pOx8REREJS2rwMrOOzrlvQ1+eCHyczOf7GT91IeWFRbxe6BtMMWBInwJuKe5B0S7bcMVzc32vCwes+oKZgpeIiEh2a5WoG5vZM8B7QFczW2pm5wG3mdl8M/sIOBy4MlHPj1VDXeQdMGOB1/2iuHcBBQFvEeaYBfbZiuU5IiIi0vIl8q3G03wOP5yo5zVWUJuGSJGhKejtwvpCV/g5IiIikt0SNuOVKUYM6Epebk6910SGpuLeBYwZ3IOC/DwMr2t8+OsgavsgIiIioC2DarVpKK+oxKjdOcIvNAX16vKr8crPy+WG47upvktEREQUvKB2kGpsKwj12RIREZGGJLSPV3MpKipypaWlqR6GiIiISIPq6+OV9TVeIiIiIsmi4CUiIiKSJApeIiIiIkmi4CUiIiKSJApeIiIiIkmi4CUiIiKSJApeIiIiIkmi4CUiIiKSJBnRQNXMVgBLEnT77YAfEnRvaTz9XNKXfjbpSz+b9KWfTXpK1M9lF+dcB78TGRG8EsnMSoO6y0rq6OeSvvSzSV/62aQv/WzSUyp+LlpqFBEREUkSBS8RERGRJFHwggmpHoD40s8lfelnk770s0lf+tmkp6T/XLK+xktEREQkWTTjJSIiIpIkWRu8zGygmS00sy/NbFSqxyMeM9vZzGaY2adm9omZ/SnVY5JNzCzHzMrMbHKqxyKbmFm+mU0yswVm9pmZHZjqMYnHzK4M/Vn2sZk9Y2abp3pM2crMHjGz5Wb2ccSxbcxsmpl9Efpn+0SPIyuDl5nlAPcBRwN7A6eZ2d6pHZWEVAPDnXN7A32BS/SzSSt/Aj5L9SCkjruBN5xzewI90c8oLZhZAXA5UOSc6w7kAKemdlRZ7VFgYNSxUcDbzrndgbdDXydUVgYvYH/gS+fcV8659cCzwAkpHpMAzrlvnXNzQr//Fe8vkILUjkoAzGwn4FjgoVSPRTYxs62BQ4CHAZxz651zFSkdlERqDeSZWWugLbAsxePJWs65d4Cfog6fADwW+v1jQHGix5GtwasA+Cbi66XoL/e0Y2aFQG/ggxQPRTx3AVcDG1I8DqmtC7ACmBhaBn7IzNqlelACzrly4G/A18C3wM/OuTdTOyqJsoNz7tvQ778Ddkj0A7M1eEmaM7MtgBeBK5xzv6R6PNnOzAYBy51zs1M9FqmjNbAvcL9zrjewmiQsl0jDQvVCJ+CF405AOzM7M7WjkiDOa/OQ8FYP2Rq8yoGdI77eKXRM0oCZ5eKFrqeccy+lejwCQD/geDNbjLc039/MnkztkCRkKbDUOReeGZ6EF8Qk9X4HLHLOrXDOVQEvAb9N8Ziktu/NrCNA6J/LE/3AbA1eHwK7m1kXM9sMr9jxlRSPSQAzM7xalc+cc3ekejzicc6Nds7t5JwrxPvvZbpzTv/nngacc98B35hZ19ChI4BPUzgk2eRroK+ZtQ392XYEevEh3bwCnB36/dnAy4l+YOtEPyAdOeeqzexSYCreWyaPOOc+SfGwxNMPOAuYb2ZzQ8eucc69lrohiaS9y4CnQv8j+RVwborHI4Bz7gMzmwTMwXtjuwx1sE8ZM3sGOAzYzsyWAtcDY4Hnzew8YAlwcsLHoc71IiIiIsmRrUuNIiIiIkmn4CUiIiKSJApeIiIiIkmi4CUiIiKSJApeIiIiIkmi4CUiWcXMdjazRWa2Tejr9qGvC1M8NBHJAgpeIpJVnHPfAPfj9e8h9M8JzrnFKRuUiGQN9fESkawT2pZqNvAIcAHQK7Sli4hIQmVl53oRyW7OuSozGwG8ARyl0CUiyaKlRhHJVkcD3wLdUz0QEckeCl4iknXMrBdwJNAXuNLMOqZ2RCKSLRS8RCSrmJnhFddf4Zz7GhgP/C21oxKRbKHgJSLZ5gLga+fctNDX/wD2MrNDUzgmEckSeqtRREREJEk04yUiIiKSJApeIiIiIkmi4CUiIiKSJApeIiIiIkmi4CUiIiKSJApeIiIiIkmi4CUiIiKSJApeIiIiIkny/wG4xEjE34PjewAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make predictions using the testing set\n",
    "y_pred = model.predict(x_test)\n",
    "\n",
    "# The coefficients\n",
    "print('True model: y ~ N(μ, σ^2)')\n",
    "print('True μ =', beta_0, '+', beta_1, '* x')\n",
    "print('Regression μ = %.2f' % model.params[0], '+ %.2f' % model.params[1], '* x')\n",
    "\n",
    "# The mean squared error\n",
    "print('Mean squared error: %.2f' % mean_squared_error(y_test, y_pred))\n",
    "\n",
    "# The coefficient of determination: 1 is perfect prediction\n",
    "print('Coefficient of determination: %.2f' % model.rsquared)\n",
    "\n",
    "# Plot outputs\n",
    "plt.figure(figsize=(10,5));\n",
    "plt.scatter(pd.DataFrame(x_test)[1], y_test)\n",
    "plt.plot(pd.DataFrame(x_test)[1], y_pred, color='red', linewidth=3);\n",
    "plt.title('Linear Regression');\n",
    "plt.xlabel('X');\n",
    "plt.ylabel('Y');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.10 64-bit (system)",
   "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.9.10"
  },
  "vscode": {
   "interpreter": {
    "hash": "9805bccb659ffcb6952b8b29440036fac091f8170d71a5a402ad9a8f2988209f"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}