04_05_Regression_TP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 04. / 05. TP Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### enseignant: Anastasios Giovanidis 2022 - 2023\n",
    "#### date: 11 Octobre 2022 and 18 Octobre 2022"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### student name or binome:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the TP related to regression. We will need to import the following libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import math\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will work with linear models of the form \n",
    "\n",
    "$y = \\beta_0+\\beta_1 x + \\epsilon$. \n",
    "\n",
    "In the following, the random error $\\epsilon$ gets a Normal distribution $\\mathcal{N}(0,\\sigma^2)$, where $\\sigma$ is the standard deviation. It is independent of $x$ and every realisation is i.i.d. We remind that $Var(\\epsilon)=\\sigma^2$.\n",
    "\n",
    "As a means of example let us consider the linear model $y=b_0+b_1x+\\epsilon$, with $b_0=2$ and $b_1=3$. Also let the noise standard deviation be $\\sigma=30$. (Note that: at the end of the example, you can change these values to your own choice)\n",
    "\n",
    "$y(x)$ is what we measure, when the single feature takes value $x$. \n",
    "\n",
    "The sample size is $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 20 # sample size\n",
    "b0 = 2 # real intercept (unobserved)\n",
    "b1 = np.array([3]) # real slope (unobserved)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation: Generating \"synthetic\" data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For the Error Vector we need to define:\n",
    "mue, sigmae = 0, 30 # error mean and standard deviation\n",
    "# If you want to fix a 'random' vector and not have different vectors per call, then use:\n",
    "np.random.seed(199)\n",
    "Er = np.random.normal(mue, sigmae, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 33.24214448  -9.66833454  68.66231565 -55.3634741  -46.02241369\n",
      "  32.69834474  36.28414203  -8.77906033 -28.89190806   6.07129635\n",
      "  48.70306245 -69.33409712  31.09989668  70.21805685  -6.69218164\n",
      " -19.60951679  24.44625946 -12.92376476  -0.37686216  39.06455214]\n"
     ]
    }
   ],
   "source": [
    "print(Er)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We suppose that the feature input $x$ ranges within an interval $\\left[x_l,x_h\\right]$. The $n$ values of $x$ are chosen uniformly randomly within the interval as a means of simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 78.57500336 117.84471417  47.98856377 101.29831222  55.66108205\n",
      "  65.38116798  74.337341    39.82041465 117.71944956  31.32134616\n",
      " 112.80089673   5.39719021  50.41908488  99.59823395  69.68557585\n",
      "  13.90856407  13.01791564 107.00400554  53.31687507  64.69449614]\n"
     ]
    }
   ],
   "source": [
    "xl, xh = 0, 120\n",
    "np.random.seed(199)\n",
    "x = np.random.uniform(xl,xh,n)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having prepared the n-sized vectors Er, and x, we can now find the 'real' values of the response $y$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[270.96715455 345.86580798 214.62800697 250.53146256 122.96083245\n",
      " 230.84184868 261.29616502 112.6821836  326.26644061 102.03533484\n",
      " 389.10575263 -51.1425265  184.35715132 371.01275871 204.36454591\n",
      "  24.11617543  65.50000638 310.08825187 161.57376303 235.14804055]\n"
     ]
    }
   ],
   "source": [
    "y = b0+b1[0]*x+Er\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us plot the line $2+3*x$ as well as the generated points $y$ that include the error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fcdf73cf280>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZrklEQVR4nO3de5Bc5Xnn8e8zQggGjJFAEKHRzIhExAgCxjVFsHEcL4TlYozwxsCI8a5ScUp/hFTw1u4a5IkrcVKzxSZZV7xVdrLjqzYedAELS4sBAyKOE5e5SMbYCFnWBF1BIMniqtF9nvxxTkNPn9MzPT19+lz696lS9fQzPd3va8Y/vXrf95zX3B0RESmWtrQbICIijadwFxEpIIW7iEgBKdxFRApI4S4iUkAnpd0AgLPPPtu7u7vTboaISK5s3Lhxv7vPjvteJsK9u7ubDRs2pN0MEZFcMbMd1b5X87SMmU0zs2fN7MHw+Swze8zMtoaPM8teu8zMhs1si5ldO7Xmi4jIZE1mzv1OYHPZ87uB9e6+AFgfPsfMFgK9wEXAdcBXzGxaY5orIiK1qCnczawD+BjwtbLyImB5+PVy4Oay+kp3P+Lu24Bh4PKGtFZERGpS68j974DPAqNltXPdfQ9A+HhOWJ8L7Cp73e6wNoaZLTWzDWa2Yd++fZNtt4iIjGPCcDezG4G97r6xxve0mFrkBjbuPujuPe7eM3t27GKviIjUqZaR+5XATWa2HVgJXGVm3wZeNbM5AOHj3vD1u4F5ZT/fAbzcsBaLiGTI0BB0d0NbW/A4NJR2iwIThru7L3P3DnfvJlgofcLdPwWsA5aEL1sCrA2/Xgf0mtkMM5sPLACebnjLRURSNjQES5fCjh3gHjwuXZqNgJ/KFar3ANeY2VbgmvA57r4JWA28ADwC3OHuJ6baUBGRrOnvh5GRsbWRkaCeNsvC/dx7enpcFzGJSN60tQUj9kpmMDoarTeamW1095647+neMiIidersnFy9mRTuIiJ1GhiA9vaxtfb2oJ42hbuISJ36+mBwELq6gqmYrq7geV9f2i3LyI3DRETyqq8vG2FeSSN3EcmdrO4tzxKN3EUkV0p7y0tbEEt7yyGbI+i0aOQuIrmS5b3lWaJwF5Fc2blzcvVWpXAXkVzJ8t7yLFG4i0iuZHlveZYo3EUkV7K8tzxLtFtGRHInq3vLs0QjdxGRAlK4i4ikYfQYPHENvPJEIm+vaRkRkWY6cRRWzXj3+SuPw+2Nv/W6wl1EpBlOHIFVp0Trtx1J5OM0LSMikqQTh+FeiwT7Sn+b7s85bdNPTuT+OBq5i4gk4fghWN0erd96kKFV7YnfH0fH7ImINFLVUB+Bk04FgjtZ7tgRfUlXF2zfXvtHjXfMnkbuIiKNcHwEVp8Wrd92CKaNnZJpxv1xFO4iIlNx9DW4f1a0HhPqJZ2d8SP3Rt4fRwuqIiL1OHIgWCitDPbbDgdbG6sEOzTn/jgauYuITMbh/bBmdrReNqc+kdKiaX9/MBXT2RkEeyNvqaBwFxGpxeG9sObcaH2c6ZfxJH1/HE3LiIiM59CrwfRLZbDXMP0ynqTPgdXIXUQkzqE98MB50fptR2DayVN662acA6t97iIi5UZegu92ROu9R6FtekM+QvvcRUSa5eAuWBuzF7H3GLQ1Niq1z11EJGlvbIbvLYzWEwj1Eu1zFxFJyuvPBwullcHeezxYKE0o2EH73EVEGu+1n8HDl0brvcehbVpTmqB97iIijXLgWXjkA9H64hNgzZ/E0D53EZGp+NUzwfRLZbAvPhFMv1hb4nvO06CRu4gU0/4n4dEPRuuLR8HsnafN2HOeBu1zF5Fi2fcjeOzD0XpFqJc0as95GrTPXUSK75XH4YlrovUqoV7SjD3naZhwzt3MTjGzp83sOTPbZGZfCOuzzOwxM9saPs4s+5llZjZsZlvM7NokOyAiLW7XmmBOvTLYF4+Gc+rVgx2q7y1v5J7zNNSyoHoEuMrdLwXeD1xnZlcAdwPr3X0BsD58jpktBHqBi4DrgK+YWXP2F4lI69ixKgj1f/n9sfXbvaZQL2nGnvM0TBjuHng7fDo9/OPAImB5WF8O3Bx+vQhY6e5H3H0bMAxc3shGi0gL2/J/glD/Ue/YeinUJ6mvDwYHgzl2s+BxcDDfi6lQ45x7OPLeCPwG8GV3f8rMznX3PQDuvsfMzglfPhd4suzHd4e1yvdcCiwF6Mz7v39EJHmb/zc8+9+j9ToCvVLSe87TUFO4u/sJ4P1mdibwgJldPM7L4/4tFPlf390HgUEIdsvU0g4RaUGb7oHnlkXrDQj1IpvUbhl3f93MfkAwl/6qmc0JR+1zgL3hy3YD88p+rAN4uRGNFZEW8vMvwM//IlpXqNeklt0ys8MRO2Z2KvB7wC+AdcCS8GVLgLXh1+uAXjObYWbzgQXA0w1ut4gU1XP9wZx6ZbDXOafeqmoZuc8Blofz7m3Aand/0Mx+DKw2s08DO4FbANx9k5mtBl4AjgN3hNM6IiLVPfs/YPPfRusK9LroClURSdfGz8CWL0XrCvUJ6QpVEcmeZ/4Ytv59tK5QbwiFu4g015N/CC9+M1pXqDeUwl1EmuNHfbDj3mhdoZ4IhbuIJOtfboFd90frCvVE6bAOkYLJzMET/3xTsKWxMti1pbEpNHIXKZBMHDzxxLXwyqPRugK9qbQVUqRAUj144vHfhb0/jNYV6onRVkiRFpHKwRPfvwJ+9VS0rlBPlcJdpEA6O+NH7oncePXhy+C1n0brCvVM0IKqSMZNZoG0KQdPPLgwWCitDHYtlGaKRu4iGTbZBdJSrb8/mIrp7AyCvSGLqWvPh4PbonUFeiZpQVUkw1JdIC1ZMwcOvxKtK9RTpwVVkZxKZYG05L6ZcOz1aF2hnguacxfJsGoLoYmeTLmqPZhTrwz2DM+pZ+bCrQzRyF0kwwYGxs65QwILpCUrpoGPRusZDfSSTFy4lUEauYtkWF8fDA4Gc+xmwePgYIND614L/lQGe4ZH6uX6+8f+5QfB8/7+dNqTFVpQFWlV98adZU8uAr1cWxvExZgZjMb8Q6RItKAqIu8qSKiXNPXCrRzRtIxIqyhNv1TKyfRLNU25cCuHNHIXKbqCjdQrJXrhVo5p5C5Sp8xvvyvoSD1OX19wUdfoaPDY6sEOGrmL1CXT2+8KPlKX2ijcRepQbfvdkiXB16kEvEJdyijcRepQ7fL/EydSGMEr1CWG5txF6jDeNrumXUDTQnPqMnkKd5E6xG2/K5fojb0U6lIDTcuI1KE05bJkSTAVUymRC2g0/SKToHAXqVMp4BO/sZdCXeqgcBeZgkQvoFGoyxQo3EWmqK8vgbs0xlGoyyRoQVUyIfNXezaDFkqlgTRyl9Rl+mrPZtBIXRKg+7lL6jJxCHQaFOoyRbqfu2RaqodAp0GhLk2gOXdJXSqHQKehypz6kGtOXRpPI3dJXVMPgU5DlZG69QWBXrrStSXWF6RpJhy5m9k8M/snM9tsZpvM7M6wPsvMHjOzreHjzLKfWWZmw2a2xcyuTbIDkn9NOQQ6DVVG6tbn7wQ76DBnScaEC6pmNgeY4+4/MbP3ABuBm4E/AA64+z1mdjcw093vMrOFwArgcuA84HHgAnePuUg7oAVVKZRx5tRb+TBnabwpLai6+x5gT/j1W2a2GZgLLAI+Gr5sOfAD4K6wvtLdjwDbzGyYIOh/PLVuiGRcDQulOsxZmmVSC6pm1g1cBjwFnBsGf+kvgHPCl80FdpX92O6wVvleS81sg5lt2LdvXx1NF8mISVx8pMOcpVlqXlA1s9OB7wCfcfc3zaqMUiDuG5F/iLr7IDAIwbRMre0QyYw6tjTqMGdplprC3cymEwT7kLuvCcuvmtkcd98TzsvvDeu7gXllP94BvNyoBoukLi7Up50Ctx2q6ccbfi8akRi17JYx4OvAZnf/Ytm31gHhiZEsAdaW1XvNbIaZzQcWAE83rskiKRlv+qXGYBdpllrm3K8E/jNwlZn9NPxzA3APcI2ZbQWuCZ/j7puA1cALwCPAHePtlBEZTyZuKBYX6m0zdEMvybRadsv8K/Hz6ABXV/mZAUBLRDIlqd9QLG6UfvJM+OSBJny4yNToxmGSWandUCwu1E+dA5/Q0pFki24cJrnU9BuKxYX6afNh0YsJfaBIchTukllNu+AnLtTPuBBufKHBHyTSPLorpGRW4hf8xC2UzrwsWCRVsEvOaeQumZXYBT9xI/WzPwT/8UdTfGOR7FC4S6Y19IKfuFA/9yq4en2DPkAkOxTuUnxxoT7nevgPDzW/LSJNonCX4ooL9Y5PwEfWROsiBaMFVcm12CtY4xZKu3qDhVIFu7QIjdwltyqvYN3+P2NG6hd+Fi77X81tmEgGKNwlt/r7g2D3oZhQv+jP4NK/an6jRDJC4S65FTdS//x9f8nA2s/ryDppeQp3yZ+YhdI/+upX+foP/ggI7j0j0uoU7pIfMaG+9Bvf4qvrl7zzXEfWiQQU7pJ9cVsaPzQE3bfzuw6PDuvIOpFKCnfJJndYEbNT98pV0HXrO091ZJ1IPIW7ZEu1UP+dNTDvE81vj0hOKdwlG6qF+kfWQcfHm98ekZxTuEu6qoX6Rx+C865vfntECkLhLunwUVgxLVq/aj382lXNb49IwSjcpbmqhfq1z8BZsUdBikgdFO7SHKMnYGXMr9t1P4FZlzW/PSIFp3CXZFUL9eufg5mXNL89Ii1C4S7JGD0OK6dH6zf8HM68uPntEWkxCndprNFjsPLkaP1jL8B7L2x+e0RalMJdGuPEUVg1I1q/cQuccUHz2yPS4nQSk0zNiSPBvV8qg/3jw8HJRykGe+wpTSItQiN3qc+Jw7Dq1Gj9phfh9PnNb0+FylOaduwInoPuRSOtwdw97TbQ09PjGzZsSLsZUosTR2DVKdH6ou1wWnZupN7dHQR6pa4u2L692a0RSYaZbXT32AtENHKX2lQbqd+8C9o7mt+eCezcObm6SNEo3GV8xw/B6vZIee6fvMT0M85jwLM5zdHZGT9y7+xsfltE0qBwl3jHR2D1aZFyx52/4qX9s4Inr2V3HntgYOycO+iUJmkt2i0jYx0/GOx+qQz2Tx6g+3P+brCHRkagv7+J7atRXx8MDgZz7GbB4+Bg9v4SEkmKFlQlcOwtuO+MaP2Tr8PJ7wWCLYVxvy5mMDqabPNEJEoLqlLdsTfhvvdG67e8AdPHhr3msUXyQ9MyBVTTxTtH3wimXyqD/Za3gouPpkdH8QMDwbx1Oc1ji2TThOFuZt8ws71m9nxZbZaZPWZmW8PHmWXfW2Zmw2a2xcyuTarhEq908c6OHcEUSuninXcC/uhrQajff+bYH7z17TDUT6/63prHFsmPCefczewjwNvA/3P3i8PaXwMH3P0eM7sbmOnud5nZQmAFcDlwHvA4cIG7nxjvMzTn3jjVLt655H0HeO7zZ0W/cetBOCm61VFEsm9Kc+7u/kMz664oLwI+Gn69HPgBcFdYX+nuR4BtZjZMEPQ/rqvlMmmVF+mcdfp+9v/f2dEX3joCJ8VclCQihVDvguq57r4HwN33mNk5YX0u8GTZ63aHtQgzWwosBejUilzDlBY9z37PPvb9wznRF9x2CKbF3D5ARAql0QuqFlOLnfdx90F373H3ntmzY0aWUpe/HngdH7JIsK/ww8GcuoJdpCXUG+6vmtkcgPBxb1jfDcwre10H8HL9zZOahQult767tg3Ags8fYcidxX0x91oXkcKqd1pmHbAEuCd8XFtWv9fMvkiwoLoAeHqqjZRxHDkA34lZKO09Bm0nsfX25jdJRNI3Ybib2QqCxdOzzWw38OcEob7azD4N7ARuAXD3TWa2GngBOA7cMdFOGanT4f2wJmY6q/c4tE1rfntEJFNq2S2zuMq3rq7y+gFAl7Uk5fA+WBOzUKpQF5Eyuv1AXhzeC2vOjdYXnwDThcYiMpbCPesOvQIPzInWFeoiMg6Fe1Yd2gMPnBetK9RFpAYK96wZeQm+G3Ns3eLR4IYuIiI1ULhnxchu+O68aF2hLiJ1ULin7eBOWNsVrSvURWQKFO5peXs7rJsfrSvURaQBtDJXpqZDLqbq7W3B/dQrg33xaHDvFwW7iDSARu6h0iEXIyPB89IhF9Cgwyje+jf4/78RrWukLiIJ0AHZoWqHXHR1wfbtU3jjN7fCgxdUFA1u14nSIjI1OiC7BpWHXExUn9Cbv4QHf3NsrW069B6t8w1FRGqncA+VDrmIq0/KG7+A7104tjbtVLhtpO62iYhMlhZUQwMD0F5xlGh7e1CvyRsvBAul5cE+/YxgkVTBLiJNppF7qLRo2t8fTMV0dgbBPuFi6uvPw0O/NbY24yz4/f2JtFNEpBYK9zJ9fZPYGfPaz+DhS8fWTvk1+E97Gt4uEZHJUrhP1ms/hYcvG1tr74Cbd6XSHBGROAr3Wh3YCI9U7Dg6rRsWbUulOSIi41G4T+RXz8D3Lx9be88C+Pgv02mPiEgNFO7V7H8KHr1ibO29C+Fjm9Jpj4jIJCjcK73+c3jokrG1My+BG55Lpz0iInVQuJcc3Anrfh38+Lu1mR+A6zem1yYRkTop3A/ugHXng5fd6+WCP4WeL6XXJhGRKWrdcI+7n/pvfw1+/dOpNEdEpJFaL9zffjGYfil3xTfh/D9IpTkiIklonXCPu5/6Fcvh/P+STntERBJU/BuHvTUc3NCrPNg/+I/BDb3GCfamnMokIpKQ4o7c4w7J+OC3Yf7EN49J/FQmEZGEFe8kpje3wIPvG1v70L3Qvbjmt0jsVCYRkQZqjZOY4g7JuHIVdN066bdq+KlMIiJNlv9wf2MzfG/h2NqH74POT9b9lg07lUlEJCX5XlA9cmBssH/4/mChdArBDg04lUlEJGX5HrmffCZc9DmY1QPzPtGwt637VCYRkYwo3oKqiEiLGG9BNd/TMiIiEkvhLiJSQAp3EZECSizczew6M9tiZsNmdndSnyMiIlGJhLuZTQO+DFwPLAQWm9nC8X9KREQaJamR++XAsLu/6O5HgZXAooQ+S0REKiQV7nOBXWXPd4e1d5jZUjPbYGYb9u3bl1AzRERaU1LhbjG1MRvq3X3Q3XvcvWf27NkJNUNEpDUlFe67gXllzzuAlxP6LBERqZBUuD8DLDCz+WZ2MtALrEvos0REpEIi95Zx9+Nm9ifA94FpwDfcfVMSnyUiIlGJ3TjM3R8CHkrq/UVEpDpdoSoiUkAKdxGRAlK4i4gUkMJdRKSAFO4iIgWkcBcRKSCFu4hIASncRUQKSOEuIlJACncRkQJSuIuIFJDCXUSkgBTuIiIFpHAXESmgXIf70BB0d0NbW/A4NJR2i0REsiGx+7knbWgIli6FkZHg+Y4dwXOAvr702iUikgW5Hbn3978b7CUjI0FdRKTV5Tbcd+6cXF1EpJXkNtw7OydXFxFpJbkN94EBaG8fW2tvD+oiIq0ut+He1weDg9DVBWbB4+CgFlNFRCDHu2UgCHKFuYhIVG5H7iIiUp3CXUSkgBTuIiIFpHAXESkghbuISAGZu6fdBsxsH3AQ2J92WxJwNupXnqhf+VPUvtXSry53nx33jUyEO4CZbXD3nrTb0WjqV76oX/lT1L5NtV+alhERKSCFu4hIAWUp3AfTbkBC1K98Ub/yp6h9m1K/MjPnLiIijZOlkbuIiDSIwl1EpIBSD3czu87MtpjZsJndnXZ76mVm88zsn8xss5ltMrM7w/osM3vMzLaGjzPTbms9zGyamT1rZg+Gz4vSrzPN7H4z+0X43+6DReibmf3X8PfweTNbYWan5LFfZvYNM9trZs+X1ar2w8yWhVmyxcyuTafVE6vSr78Jfw9/ZmYPmNmZZd+bdL9SDXczmwZ8GbgeWAgsNrOFabZpCo4D/83dLwSuAO4I+3I3sN7dFwDrw+d5dCewuex5Ufr1JeARd38fcClBH3PdNzObC/wp0OPuFwPTgF7y2a9vAddV1GL7Ef7/rRe4KPyZr4QZk0XfItqvx4CL3f0S4JfAMqi/X2mP3C8Hht39RXc/CqwEFqXcprq4+x53/0n49VsEITGXoD/Lw5ctB25OpYFTYGYdwMeAr5WVi9CvM4CPAF8HcPej7v46BegbwVkNp5rZSUA78DI57Je7/xA4UFGu1o9FwEp3P+Lu24BhgozJnLh+ufuj7n48fPok0BF+XVe/0g73ucCusue7w1qumVk3cBnwFHCuu++B4C8A4JwUm1avvwM+C4yW1YrQr/OBfcA3wymnr5nZaeS8b+7+EvC3wE5gD/CGuz9KzvtVplo/ipQnfwg8HH5dV7/SDneLqeV6b6aZnQ58B/iMu7+ZdnumysxuBPa6+8a025KAk4APAH/v7pcR3N8oD1MV4wrnoBcB84HzgNPM7FPptqopCpEnZtZPMM07VCrFvGzCfqUd7ruBeWXPOwj++ZhLZjadINiH3H1NWH7VzOaE358D7E2rfXW6ErjJzLYTTJtdZWbfJv/9guD3b7e7PxU+v58g7PPet98Dtrn7Pnc/BqwBPkT++1VSrR+5zxMzWwLcCPT5uxch1dWvtMP9GWCBmc03s5MJFg3WpdymupiZEczdbnb3L5Z9ax2wJPx6CbC22W2bCndf5u4d7t5N8N/nCXf/FDnvF4C7vwLsMrPfDEtXAy+Q/77tBK4ws/bw9/JqgjWgvPerpFo/1gG9ZjbDzOYDC4CnU2hfXczsOuAu4CZ3Hyn7Vn39cvdU/wA3EKwM/xvQn3Z7ptCPDxP8U+lnwE/DPzcAZxGs6G8NH2el3dYp9PGjwIPh14XoF/B+YEP43+27wMwi9A34AvAL4HngH4EZeewXsIJg3eAYwQj20+P1A+gPs2QLcH3a7Z9kv4YJ5tZL+fEPU+mXbj8gIlJAaU/LiIhIAhTuIiIFpHAXESkghbuISAEp3EVECkjhLiJSQAp3EZEC+ncnvXo4Vu1NFwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,y,'o',color='blue')\n",
    "plt.plot(x,b0+b1*x,color='orange')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Sample generator is written in what follows as a function, to be called at later stages. When the seed = seedit takes an integer value, the random vectors do not change on each call. By letting seed = () the random vectors $x,y,Er$ change on each call. The parameter size determines the size of vector $x$ i.e. the number of different features included (if $size>1$ then the model is multi-linear)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def genSample(n,b0,b1,sigmae,xLow,xHigh,seedit=199,size=1):\n",
    "    if type(seedit)==int:\n",
    "        np.random.seed(seedit)\n",
    "        Er = np.random.normal(mue, sigmae, n)\n",
    "        x = []\n",
    "        for k in range(size):\n",
    "            np.random.seed(seedit+k)\n",
    "            x.append(np.random.uniform(xl,xh,n))\n",
    "    else:\n",
    "        np.random.seed()\n",
    "        Er = np.random.normal(mue, sigmae, n)\n",
    "        x = []\n",
    "        for k in range(size):\n",
    "            np.random.seed()\n",
    "            x.append(np.random.uniform(xl,xh,n))\n",
    "    y = b0+Er\n",
    "    for k in range(size):\n",
    "        y +=b1[k]*x[k]\n",
    "    #\n",
    "    # Output\n",
    "    if size==1:\n",
    "        return (x[0], y, Er)\n",
    "    else:\n",
    "        return (x, y, Er)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x= [ 78.57500336 117.84471417  47.98856377 101.29831222  55.66108205\n",
      "  65.38116798  74.337341    39.82041465 117.71944956  31.32134616\n",
      " 112.80089673   5.39719021  50.41908488  99.59823395  69.68557585\n",
      "  13.90856407  13.01791564 107.00400554  53.31687507  64.69449614]\n",
      "y= [270.96715455 345.86580798 214.62800697 250.53146256 122.96083245\n",
      " 230.84184868 261.29616502 112.6821836  326.26644061 102.03533484\n",
      " 389.10575263 -51.1425265  184.35715132 371.01275871 204.36454591\n",
      "  24.11617543  65.50000638 310.08825187 161.57376303 235.14804055]\n"
     ]
    }
   ],
   "source": [
    "n, b0, b1 = 20, 2, np.array([3]) \n",
    "(x,y,Er) = genSample(n,b0,b1,sigmae,xl,xh,seedit=199,size=1)\n",
    "print(\"x=\",x)\n",
    "print(\"y=\",y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can also put the $x,Er,y$ vectors in a panda table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>78.575003</td>\n",
       "      <td>270.967155</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>117.844714</td>\n",
       "      <td>345.865808</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>47.988564</td>\n",
       "      <td>214.628007</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>101.298312</td>\n",
       "      <td>250.531463</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>55.661082</td>\n",
       "      <td>122.960832</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>65.381168</td>\n",
       "      <td>230.841849</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>74.337341</td>\n",
       "      <td>261.296165</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>39.820415</td>\n",
       "      <td>112.682184</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>117.719450</td>\n",
       "      <td>326.266441</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>31.321346</td>\n",
       "      <td>102.035335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>112.800897</td>\n",
       "      <td>389.105753</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>5.397190</td>\n",
       "      <td>-51.142527</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>50.419085</td>\n",
       "      <td>184.357151</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>99.598234</td>\n",
       "      <td>371.012759</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>69.685576</td>\n",
       "      <td>204.364546</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>13.908564</td>\n",
       "      <td>24.116175</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>13.017916</td>\n",
       "      <td>65.500006</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>107.004006</td>\n",
       "      <td>310.088252</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>53.316875</td>\n",
       "      <td>161.573763</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>64.694496</td>\n",
       "      <td>235.148041</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             x           y\n",
       "0    78.575003  270.967155\n",
       "1   117.844714  345.865808\n",
       "2    47.988564  214.628007\n",
       "3   101.298312  250.531463\n",
       "4    55.661082  122.960832\n",
       "5    65.381168  230.841849\n",
       "6    74.337341  261.296165\n",
       "7    39.820415  112.682184\n",
       "8   117.719450  326.266441\n",
       "9    31.321346  102.035335\n",
       "10  112.800897  389.105753\n",
       "11    5.397190  -51.142527\n",
       "12   50.419085  184.357151\n",
       "13   99.598234  371.012759\n",
       "14   69.685576  204.364546\n",
       "15   13.908564   24.116175\n",
       "16   13.017916   65.500006\n",
       "17  107.004006  310.088252\n",
       "18   53.316875  161.573763\n",
       "19   64.694496  235.148041"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lucy = {'x': x, 'y': y}\n",
    "df_slr = pd.DataFrame(data=lucy)\n",
    "df_slr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Part A.** \n",
    "With the above synthetic data: \n",
    "* a) Compute the sample means $\\bar{x}$, $\\bar{y}$, and the sample variance $s_{xx}$ and covariance $s_{xy}$ .\n",
    "\n",
    "* b) Find the coefficient estimates for the least-squares line $\\hat{y}=\\hat{\\beta}_0+\\hat{\\beta}_1x$, using the formulas presented in the course. Compare these with the unobserved real values of $(\\beta_0,\\beta_1)=(2,3)$. What do you observe?\n",
    "\n",
    "* c) Plot the regression line, least-squares line and the data samples.\n",
    "\n",
    "You can work with either numpy or pandas.\n",
    "\n",
    "You can compare your results with the ones from the **Simple Linear Regression using scikit-learn** that follows."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note** For numpy commands:\n",
    "\n",
    "- Command x.var() in numpy has by default 'ddof=0', so that the denominator is 'n' (Population Statistic). To change this to sample variance, one needs to set the parameter 'ddof=1', i.e. Sxx = x.var(ddof=1).\n",
    "\n",
    "- Command cov(x,y) in numpy has by default 'ddof=1', so that the denominator is 'n-1' (Sample Statistic). To change this to population covariance one needs to set the parameter 'ddof=0', i.e. Sxy = np.cov(x,y,ddof=0)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer 1.A** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11f43ff28>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt4VNW9//H3N0C5KILlYikICXKXu0FRKOUggoo/UBQQoqK1Ymtb23NaqpbzU7zQarXa01OljVVBTQW1iIioCEIBWy+AqMj9EiCA3JSLRJTLOn/snZjJTMgkmZk9M/m8nifPzF6zs2ftTPhksdbae5lzDhERSV8ZQVdARETiS0EvIpLmFPQiImlOQS8ikuYU9CIiaU5BLyKS5hT0IiJpTkEvIpLmFPQiImmuZtAVAGjcuLHLzMwMuhoiIill2bJle51zTcrbLymCPjMzk6VLlwZdDRGRlGJmW6LZL+quGzOrYWYfmNlsfzvLzN41s/VmNt3MvuWX1/a3N/ivZ1bmBEREJDYq0kf/c2B1ie0HgEecc22Bz4Eb/fIbgc+dc22AR/z9REQkIFEFvZm1AIYAf/O3DRgAvOjvMhW43H8+zN/Gf/1Cf38REQlAtH30fwR+DdT3txsB+51zx/ztAqC5/7w5sA3AOXfMzA74++8teUAzGweMA2jZsmXYGx49epSCggKOHDkS9cmIJEKdOnVo0aIFtWrVCroqIlEpN+jN7DJgt3NumZn1LyqOsKuL4rVvCpzLBXIBsrOzw14vKCigfv36ZGZmov8QSLJwzrFv3z4KCgrIysoKujoiUYmm66YPMNTM8oFpeF02fwQamlnRH4oWwA7/eQFwJoD/egPgs4pW7MiRIzRq1EghL0nFzGjUqJH+pynF8vIgMxMyMrzHvLygaxSu3KB3zt3hnGvhnMsErgbecs7lAAuAq/zdxgIv+89n+dv4r7/lKrmMlUJekpF+L6VIXh6MGwdbtoBz3uO4cckX9lW5MvY24L/MbANeH/wTfvkTQCO//L+A26tWRRGR5DRhAhQWhpYVFnrlyaRCF0w55xYCC/3nm4BzI+xzBBgRg7olnYkTJ3Lqqafyq1/9KuLrM2fOpF27dnTq1CnBNRORIGzdWrHyoKTNvW6SoZ9s5syZrFq1KvFvLCKBiDBh8KTlQUmLoI9nP9mkSZNo3749AwcOZO3atQA8/vjj9OrVi27dunHllVdSWFjIv/71L2bNmsX48ePp3r07GzdujLifiKSPSZOgXr3Qsnr1vPJkkhZBH69+smXLljFt2jQ++OADZsyYwfvvvw/A8OHDef/99/nwww/p2LEjTzzxBBdccAFDhw7lwQcfZMWKFZx11lkR9xOR9JGTA7m50KoVmHmPubleeTJJipuaVVW8+skWL17MFVdcQT3/T/bQoUMBWLlyJf/93//N/v37+eKLLxg8eHDE7492PxFJXTk5yRfspaVFiz6e/WSRptJdf/31/PnPf+bjjz/mrrvuKnNOdbT7iUj0kmE8LtWkRdDHq5+sX79+vPTSS3z55ZccOnSIV155BYBDhw7RrFkzjh49Sl6J37L69etz6NCh4u2y9hORykmVeevJJi2CPl79ZD179mTUqFF0796dK6+8ku9973sA3HvvvZx33nlcdNFFdOjQoXj/q6++mgcffJAePXqwcePGMvcTkcpJlXnrycYqedFqTGVnZ7vSC4+sXr2ajh07BlQjkZPT72cwMjK8lnxpZnDiROLrEzQzW+acyy5vv7Ro0YtI9ZAq89aTjYJeRFJGqsxbTzYKehFJGakybz3ZpMU8ehGpPlJh3nqyUYteRCTNKehFRIJw7HDC3kpBn2LuvPNO5s2bV+XjLFy4kMsuuwyAWbNmcf/991f5mOAttTdgwAAOHjwYk+MVmTJlCj/96U9jeszSSv5MZs+ezV133RXX95Nq6sBq+LvB86fC5sRc6aWgj4JzjhNVnKR77Nix8neKwj333MPAgQNjcqwiQ4cO5fbbY7M+zJw5c+jWrRunnXZaTI4XC8ePH6/w9wwZMoRZs2bpjqMSO/tXegH/aon1KiwxEaygL0N+fj4dO3bklltuoWfPnmzbto25c+dy/vnn07NnT0aMGMEXX3wBeOHWoUMH+vbty6233lrcKpw4cSLjxo1j0KBBXHfddRw/fpzx48fTq1cvunbtyl//+lcAdu7cSb9+/ejevTudO3dm8eLFHD9+nOuvv57OnTvTpUsXHnnkEcC7f86LL74IwPz58+nRowddunThBz/4AV999RUAmZmZ3HXXXfTs2ZMuXbqwZs2ak55rydby9ddfz6233soFF1xA69ati98L4MEHHyyue1mt3by8PIYNG1b8M+zcuXPxaw899BATJ04EoH///tx2222ce+65tGvXjsWLFwPeWsE33HADXbp0oUePHixYsKD4+7dt28bFF19M+/btufvuuwE4fPgwQ4YMoVu3bnTu3Jnp06cX/wzuuece+vbty+9//3vOPfebNXLy8/Pp2rUrAK+//nrxZzdjxozifcyM/v37M3v27JP+7ETKdWQ3PF8f5nT5puyUVuQ5R2b/0Qm5Z0/qzLqZ1z+8rOVIaHcLHCuEhZeGv976eu/ryF5YclXoawMXlvuWa9eu5amnnuKxxx5j79693HfffcybN49TTjmFBx54gIcffphf//rX3HzzzSxatIisrCxGjx4dcoxly5axZMkS6tatS25uLg0aNOD999/nq6++ok+fPgwaNIgZM2YwePBgJkyYwPHjxyksLGTFihVs376dlStXArB///6Q4x45coTrr7+e+fPn065dO6677jomT57ML37xCwAaN27M8uXLeeyxx3jooYf429/+Vu75Ftm5cydLlixhzZo1DB06lKuuuoq5c+eyfv163nvvPZxzDB06lEWLFtGvX7+Q73377beL/4CV59ixY7z33nvMmTOHu+++m3nz5vHoo48C8PHHH7NmzRoGDRrEunXrAHjvvfdYuXIl9erVo1evXgwZMoQtW7bw3e9+l1dffRWAAwcOFB+/Tp06LFmyBIDp06ezadMmWrduzfTp0xk5ciRHjhzhpptu4q233qJNmzaMGjUqpH7Z2dksXryYkSNHRv2zEyn25aewZAQcWg/HvEYhdZvDFQXF9+wp+g9j0T17ID4zitSiP4lWrVrRu3dvAN555x1WrVpFnz596N69O1OnTmXLli2sWbOG1q1bk5WVBRAW9EOHDqVu3boAzJ07l6effpru3btz3nnnsW/fPtavX0+vXr146qmnmDhxIh9//DH169endevWbNq0iZ/97Ge8/vrrYV0ha9euJSsri3bt2gEwduxYFi1aVPz68OHDATjnnHPIz8+v0HlffvnlZGRk0KlTJ3bt2lVc97lz59KjRw969uzJmjVrWL9+fdj3fvbZZ9SvXz+q94lUxyVLlnDttdcC0KFDB1q1alUc9BdddBGNGjWibt26DB8+nCVLltClSxfmzZvHbbfdxuLFi2nQoEHx8UsG98iRI3n++ecBL/RHjRrFmjVryMrKom3btpgZ11xzTUj9mjZtyo4dO6I6F5Finy33umheagZ7lsDp3aHvCzDGwRUFQOLv2ZM6LfqTtcBr1jv563UaR9WCL+2UU04pfu6c46KLLuK5554L2eeDDz6o0DH+93//N+J96RctWsSrr77Ktddey/jx47nuuuv48MMPeeONN3j00Ud5/vnnefLJJ0OOdTK1a9cGoEaNGhUeHyj63pLv45zjjjvu4Oabbz7p99asWZMTJ06QkZFR/LxI6ds0R6rjyc6r9C2jzYx27dqxbNky5syZwx133MGgQYO48847gdCf/ahRoxgxYgTDhw/HzGjbti0rVqyIeBvqkvUt+iMtUq79K0O7ZwC63gedw9M70WvNqkUfpd69e/P222+zYcMGAAoLC1m3bh0dOnRg06ZNxS3Soj7iSAYPHszkyZM5evQoAOvWrePw4cNs2bKFpk2bctNNN3HjjTeyfPly9u7dy4kTJ7jyyiu59957Wb58ecixOnToQH5+fnF9nnnmGb7//e/H4cy/qfuTTz5ZPC6xfft2du/eHbZf+/bt2bRpEwBnnHEGu3fvZt++fXz11VdR9Xf369ev+JbO69atY+vWrbRv3x6AN998k88++4wvv/ySmTNn0qdPH3bs2EG9evW45ppr+NWvfhX2cypy1llnUaNGDe69997iln6HDh3YvHkzGzduBAj7I75u3bqQMQaRiI4Vwj+HhYb8ae29FnyEkIfE37MndVr0AWvSpAlTpkxh9OjRxYOe9913H+3ateOxxx7j4osvpnHjxiGDfqX98Ic/JD8/n549e+Kco0mTJsycOZOFCxfy4IMPUqtWLU499VSefvpptm/fzg033FDcIv7d734Xcqw6derw1FNPMWLECI4dO0avXr340Y9+FLfzHzRoEKtXr+b8888H4NRTT+XZZ5+ladOmIfsNGTKEhQsX0qZNG2rVqsWdd97JeeedR1ZWVlS3ar7lllv40Y9+RJcuXahZsyZTpkwpbvn37duXa6+9lg0bNjBmzBiys7N54403GD9+PBkZGdSqVYvJkyeXeexRo0Yxfvx4Nm/eDHg/w9zcXIYMGULjxo3p27dv8ZgIwIIFC8J+7iLFdrz2zdhgRi3vscN/Qc8/lPutkyaF9tFDnO/Z45wL/Oucc85xpa1atSqsLFkdOnTIOefciRMn3I9//GP38MMPB1yj4OzYscMNHDgw6GpU2aeffuoGDBhQ5uup9PspMZY/zbk8Qr82PV3hwzz7rHOtWjln5j0++2zFqwIsdVFkrFr0MfD4448zdepUvv76a3r06FFuP3Y6a9asGTfddBMHDx5Mqrn0FbV161b+8IfyW2ZSjRzeAi9nhpePqdyaHom8Z48WHhGpBP1+ViO7/gnz+4eXVzLgwZszP2GCN/jasqXXZVOZ0I924ZGkbtE75046K0IkCMnQOJIE2PgEvPvDb7bPvBKaD4XW11XpsImeQw9JPOumTp067Nu3T/+oJKk459i3bx916tQJuioSLxuf8ObBlwz5MwbA916scshDMOveJm2LvkWLFhQUFLBnz56gqyISok6dOrRo0SLoakisfb0fXusJhzeHllehiyaSRM+hhyQO+lq1ahVfbSoiEjef3A8f3hFadmprGLoxLm/XsqXXXROpPF6SNuhFROJqxR2wqtTtuS9bB6e1jevbJnwOPQp6Ealu9rwNb/YNL49xF01ZigZcYzHrJloKehGpHj75HXz4m/DyBAV8SYle9zZpZ92IiMTE0lu9WTQlQ77rvV7Aj3Hk5Xn3g0/EfeGDoha9iKSn934MG/4SWlazPoz8ZpnLIOa0ByFpr4wVEamUXQtg/oDw8ghdNJmZkWfAtGoFFVzGIRBpcWWsiEjU3rkRNj0ZXn6SPvgg5rQHodygN7M6wCKgtr//i865u8wsC5gGfBtYDlzrnPvazGoDTwPnAPuAUc65/DjVX0Squ1c7w4FPQsvO+RO0/1m53xrEnPYgRDMY+xUwwDnXDegOXGxmvYEHgEecc22Bz4Eb/f1vBD53zrUBHvH3ExGJrZWTvEHWkiFfq6HXgo8i5MGb1livXmhZvOe0B6HcFr1/z2N/ZVtq+V8OGACM8cunAhOBycAw/znAi8CfzcxcMgwGiEjqm9UGvohw1WolpkkGMac9CFH10ZtZDWAZ0AZ4FNgI7HfOFS1GWgA09583B7YBOOeOmdkBoBGwt9QxxwHjAFqm2/+TRCT2XmoOX5ZarL3F5dDvpSodNtFz2oMQVdA7544D3c2sIfASEOlG3EV/TiPdVzjsT61zLhfIBW/WTVS1FZHq5+9l3Ko8gAudUlWFZt045/ab2UKgN9DQzGr6rfoWQNGf2gLgTKDAzGoCDYDPYldlEakWFlwKO18LL1fAV1g0s26aAEf9kK8LDMQbYF0AXIU382Ys8LL/LbP87X/7r7+l/nkRidor7eHQuvByBXylRdOibwZM9fvpM4DnnXOzzWwVMM3M7gM+AJ7w938CeMbMNuC15K+OQ71FJN1E6qI5rSNctirxdUkz0cy6+QjoEaF8E3BuhPIjwIiY1E5E0t/Lmd7C2yXVagAj9gdSnXSkK2NFJBgvNISjB8LL1UUTcwp6EUmsSF00p/eAS5Ynvi7VhIJeRBIjUsCriyYhFPQiEl8zz4TCgtCy2o3gyr2R95eY08IjImkkqRbR+Lt5X6VDfoxTyCeYWvQiaSJpFtGI1EXTciT0nZ7ASkhJWnhEJE0EvohGpIA/4z/gwrcS8ObVkxYeEalmAltEI1LAn5IJwzbH+Y0lWgp6kTSR8EU0IgX8dwbCgDfj9IZSWRqMFUly0Q6wJmwRjaJB1pLOvMobZFXIJyW16EWSWEUGWOO+iEakFvyZw+F7/4jRG0i8aDBWJIkFPsAKkQO+dhO4cneCKiBl0WCsSBoIbIAVIgd8+5/DOX9MwJtLLKmPXiSJlTWQGtfVNyP1wWeN9frgkzTkk+pCsSSkFr1IEps0KbSPHuI0wAqRW/AtroB+M+LwZrGTNBeKJTG16EWSWE4O5OZ6ffJm3mNubowDLFILvsXlXgs+yUMevMHnkn8IwdueMCGY+iQjDcaKVFeRWvDdfgdn3574ulRBRgZEijEzOHEi8fVJJA3GikhkkQK+y0ToclfCqxILCb9QLAUp6EWqi0gB/90h0H924usSQwkdx0hRCnqRdBcp4Hs+DB3+M/F1iYO4XyiWBjQYK1JJST+lL9Iga5d7vEHWNAn5Ijk53gVkJ054jwr5UGrRi1RCUk/pi9SCP/s30E19GdWVgl6kEsqa0jd2rPc8kLCPFPDtfgbZf0p8XSSpKOhFKqGsWxAcPx5Ayz5SwJ/3JJx1Q4IqIMlOffQilXCyqXsJu1gnUh98n+leH7xCXkpQ0ItUQqR7v5cU15uORQr4Hg95Ad9qZBzfWFKVum5EKqGoW2bsWK+7prS4XKwTqYvm+69C80vj8GaSThT0IpVUFPZxv1gnUsCfmwttborhm0g6U9CLVEFcL9aJFPAX5EHmmBgcXKoTBb1IFeXkxOFukqV1fwA6/TqGbyLViQZjJSkk/VWmiRBpkPV7M7xBVoW8VIFa9BK4pL7KNBEiteD7vQwthia+LpKWdD96CVxSLIAdhIiDrH+FNuMSXxdJSbofvaSMQBfATjTn4LkIPaYD/wlN+yW+PlItqI9eAhfIAtiJ5pzXgi8V8ovcP7w+eIW8xFG5QW9mZ5rZAjNbbWafmNnP/fJvm9mbZrbefzzdLzcz+5OZbTCzj8ysZ7xPQlJbpKtM02bhiDICfugfXsZyHJeMG149B54loaJp0R8Dfumc6wj0Bn5iZp2A24H5zrm2wHx/G+ASoK3/NQ6YHPNaS1pJyALYiVZGwF/1Py9gOY5XlnsDrVrEWhKh3D5659xOYKf//JCZrQaaA8OA/v5uU4GFwG1++dPOG+V9x8wamlkz/zgiEcV8LnpQ3Al4rkZ4+eD3yWiSHXER67Qci5CkUqHBWDPLBHoA7wJnFIW3c26nmTX1d2sObCvxbQV+mYJe0ldZAT/oHWh8HqBFrCU4UQ/GmtmpwD+AXzjnDp5s1whlYe0YMxtnZkvNbOmePXuirYZIcjlx3O+iKRXyF77lDbL6IQ9pPhYhSS2qFr2Z1cIL+Tzn3Ay/eFdRl4yZNQN2++UFwJklvr0FsKP0MZ1zuUAuePPoK1l/kWCcOAbTaoWXX/oRNOwS8Vu0iLUEpdygNzMDngBWO+ceLvHSLGAscL//+HKJ8p+a2TTgPOCA+uclbZw4CtO+FV5+yYdwetdyvz1txiIkpUTTou8DXAt8bGYr/LLf4AX882Z2I7AVGOG/Nge4FNgAFAJa6kZSX5kBvwJO75b4+ohUQDSzbpYQud8d4MII+zvgJ1Wslwjg3Qcn0K6OsgL+wgVwRv8EVkSk8nQLBElagd7s7FghPH9KePn/Ww/128T5zUViSzc1k6QVyM3Ojh+B6XXDyy9bA6e1j9ObilSObmomKS+hNzs7+gW8UD+8/NKPoWHnOLyhSOIo6CVpJeQCo7ICfshqaNAhhm8kEhzdvVKSVlwvMPr6c+9Cp9Ihf/l270InhbykEbXoJWnF5QKjowfhhQbh5VfsgLrNqnBgkeSloJekFrMLjMoK+MvWwWltY/AGIslLQS/p7ctd8NJ3wsuH74Y6TRJfH5EAKOglPX29H148Pbx8+C6o0zS8XCSNaTBWUlpenjffPiPDe3wxb683yFo65K/Y6Q2yKuSlGlKLXlJWyStnm5y2m/zfnhG+0/A9UKdx4isnkkQU9JKyJkyAhrW3c/jxFuEvjjgItSLMjxephhT0kpqO7CX/t+GDqfVuOMyRo/U4MSaAOokkKQW9pJYje2BGeD/7KT/4gsKvvJuQtWqV6EqJJDcFvaSGw1vh5fAEP/3HR9h/sHbxtpbmEwmnoJfkdngbvBzh5jajjkCN2vzZaWk+kfIo6CU5ldFFUxTwRbQ0n0j5FPSSXA5vgZczw8tLBbyIRE9BL8nh0EZ4JcLKTVcfhQz9mopUhf4FSbAKt8PMCPPgrz4GGTUSXx+RNKSgl2CUNciqFrxIzOlflCTWwfUwu114+ejjYLr1kkg8KOglMT7/CF7rFl6ugBeJOwW9xFdZffCjT4BZ4usjUg0p6CU+ypomqYAXSTgFvcTW5yvgtR7h5Qp4kcAo6CU2DqyBVzuGlyvgRQKnUTCpmv0rvRWdSof86BPeik4BhXzplafy8gKphkhSUIteKmffUnijV3j5GJf4upRScuUpgC1bvG3QfXGkejLngv+HmZ2d7ZYuXRp0NSQan38Ir3UPL0+CgC+SmemFe2mtWkF+fqJrIxI/ZrbMOZdd3n5q0Ut0yrgffDIFfJGtWytWLpLu1EcvJ/fZcq8PvlTIZ/7GkZcE/xuMpGWEOyucrFwk3SnoJbL9n3gB//o5xUVvrhyM5TgsxxX3eyfjIOekSd5KUyVp5SmpztR1I6G+2ASzzgorzvyNC+v3Liz0VndKtgHOovpo5SkRjwZjxbPzTVgwKLQs41tw9Vfe0wyI9KtiBidOJKB+IhJGg7ESnc8+gNd7hpa1HAl9p4cWtYw8k0X93iLJr9w+ejN70sx2m9nKEmXfNrM3zWy9/3i6X25m9icz22BmH5lZz7KPLPEQ9YVCuxb4ffClPqIxLizkQf3eIqksmsHYKcDFpcpuB+Y759oC8/1tgEuAtv7XOGBybKop0Si6UGjLFq+bJeKA6YHVXsDPH/BNWYPOXsCfZKpkTg7k5npz0c28x9xc9XuLpIKo+ujNLBOY7Zzr7G+vBfo753aaWTNgoXOuvZn91X/+XOn9TnZ89dHHxkkvFPr3XFgwOPzFJJwHLyLRibaPvrLTK88oCm//salf3hzYVmK/Ar9MEiDSBUEXtHub/N9aaMif9cNyW/Aikj5iPRgb6Q5WEdPEzMbhde/QUiN6MVFywLRryw/58HelblVQ8xQY+UXiKyYigapsi36X32WD/7jbLy8AziyxXwtgR6QDOOdynXPZzrnsJk2aVLIaUtKkSdDv7HdxeRYS8ntcb6/1rpAXqZYqG/SzgLH+87HAyyXKr/Nn3/QGDpTXPy8x8sVmcsz45296FxflvXszec7RJOffAVZMRIJWbteNmT0H9Acam1kBcBdwP/C8md0IbAVG+LvPAS4FNgCFwA1xqLOU9NkyeL3UWEyLK6DfDHLGBFMlEUku5Qa9c250GS9dGGFfB/ykqpWSKOx7H944N7Ts/Kch69pg6iMiSUtXxqaave/C3N6hZZ1ug+73B1MfEUl6CvpUcXAdzG4fWtbpDuj+22DqIyIpQ0Gf7Ha8BgsvDS3r/P+h6z3B1EdEUo6CPlntXgTzvh9a1vcFaHlVMPURkZSloE82B1bBq2eHlnUcDz1+H0x9RCTlKeiTxa4FoTcaAw2yikhMKOiD9ul8eGtgaFm/mdBiWDD1EZG0o6APSqQVnXo8BB1/GUx9RCRtaXHwEqJetKMqPv/Qux98yZDv8QfvXjQKeRGJA7XofUWLdhQWettFi3ZAjBbXyP87/KvUgXr+ETr8PAYHFxEpmxYH95100Y78Khx4+6vwz8tCy/q/Bt8tvWiXiEjFaHHwCoq0aMfJysu19z2Ye15oWec7oevdlTygiEjlKOh9JRftKF1eIZFm0Zz7V2gzrtJ1ExGpCgW9b9Kk0D56gHr1vPKobJsBi68MLRswD74TdpNPEZGEUtD7igZcJ0zwumtatvRCvtyB2M158O9rQsvOfwayrom8v4hIginoS8jJqcAMm73vwNzzQ8t6T4HWYyPuLiISFAV9RW18Ct79QWhZ76nQ+rpg6iMiUg4FfbQ2PwP/LhXmF70NTS4Ipj4iIlFS0Jdn6z9gSalbA/f6C7S9OZj6iIhUkIK+LDvnwoLBoWV9n4eWIyLvLyKSpBT0peVPg3+VWg994D+hab9g6iMiUkUK+iKfzoNNUyC/xJ3MLnwLzviPwKokIhILCvpdC2F+iTBvNRra/QSa9AmsSiIisVR9g/7gOvjkt7B56jdlFy+Fb58TXJ1EROKg+gX9wXWw8j7IfwZq1IH2v4A2N0ODDkHXTEQkLqpP0Be8AouGes9r1IVmg6H3U1C3WbD1EhGJs/RfYerAGvh7xjch3+BsGLoZ/uP1MkM+IStNiYgkSPq26PevhLevhgOrAH9xlUs/hoadT/ptcV9pSkQkwdKvRb9rIbw9BuZ0hQOfQMfxMHyXtyZrOSEP3t0rS96qGLztCRPiU10RkXhLnxZ96VsVdLoNOvwS6jSp2GFivdKUiEjAUj/o938CK++FrdO/KbtsDZzWvlKHi9lKUyIiSSK1u24+nQ9zOsOOV71pkpet8bpoKhny4C02Uq9eaFmFVpoSEUkyqd2ib3C2t+B2+1uhdqOYHLLSK02JiCQpc84FXQeys7Pd0qVLg66GiEhKMbNlzrns8vZL7a4bEREpl4JeRCTNxSXozexiM1trZhvM7PZ4vIeIiEQn5kFvZjWAR4FLgE7AaDPrFOv3ERGR6MSjRX8usME5t8k59zUwDRgWh/cREZEoxCPomwPbSmwX+GUiIhKAeAS9RSgLm8PgLFm+AAAGA0lEQVRpZuPMbKmZLd2zZ08cqiEiIhCfoC8Aziyx3QLYUXon51yucy7bOZfdpEnF7kcjIiLRi0fQvw+0NbMsM/sWcDUwKw7vIyIiUYj5LRCcc8fM7KfAG0AN4Enn3Cexfh8REYlOXO5145ybA8yJx7FFRKRidGWsiEiaU9CLiKQ5Bb2ISJpT0IuIpDkFvYhImlPQi4ikOQW9iEiaU9CLiKQ5Bb2ISJpT0IuIpDkFvYhImlPQi4ikOQW9iEiaU9CLiKQ5Bb2ISJpL2aDPy4PMTMjI8B7z8oKukYhIcorLwiPxlpcH48ZBYaG3vWWLtw2QkxNcvUREklFKtugnTPgm5IsUFnrlIiISKiWDfuvWipWLiFRnKRn0LVtWrFxEpDpLyaCfNAnq1Qstq1fPKxcRkVApGfQ5OZCbC61agZn3mJurgVgRkUhSctYNeKGuYBcRKV9KtuhFRCR6CnoRkTSnoBcRSXMKehGRNKegFxFJc+acC7oOmNke4DCwN+i6xEFjdF6pROeVWqr7ebVyzjUpb6ekCHoAM1vqnMsOuh6xpvNKLTqv1KLzio66bkRE0pyCXkQkzSVT0OcGXYE40XmlFp1XatF5RSFp+uhFRCQ+kqlFLyIicRB40JvZxWa21sw2mNntQdensszsTDNbYGarzewTM/u5X/5tM3vTzNb7j6cHXdfKMLMaZvaBmc32t7PM7F3/vKab2beCrmNlmFlDM3vRzNb4n9356fCZmdl/+r+HK83sOTOrk4qfmZk9aWa7zWxlibKIn495/uRnyUdm1jO4mp9cGef1oP97+JGZvWRmDUu8dod/XmvNbHBF3y/QoDezGsCjwCVAJ2C0mXUKsk5VcAz4pXOuI9Ab+Il/LrcD851zbYH5/nYq+jmwusT2A8Aj/nl9DtwYSK2q7n+A151zHYBueOeY0p+ZmTUHbgWynXOdgRrA1aTmZzYFuLhUWVmfzyVAW/9rHDA5QXWsjCmEn9ebQGfnXFdgHXAHgJ8jVwNn+9/zmJ+dUQu6RX8usME5t8k59zUwDRgWcJ0qxTm30zm33H9+CC8wmuOdz1R/t6nA5cHUsPLMrAUwBPibv23AAOBFf5dUPa/TgH7AEwDOua+dc/tJg88M7xbkdc2sJlAP2EkKfmbOuUXAZ6WKy/p8hgFPO887QEMza5aYmlZMpPNyzs11zh3zN98BWvjPhwHTnHNfOec2AxvwsjNqQQd9c2Bbie0CvyylmVkm0AN4FzjDObcTvD8GQNPgalZpfwR+DZzwtxsB+0v8Uqbq59Ya2AM85XdL/c3MTiHFPzPn3HbgIWArXsAfAJaRHp8ZlP35pFOe/AB4zX9e5fMKOugtQllKTwMys1OBfwC/cM4dDLo+VWVmlwG7nXPLShZH2DUVP7eaQE9gsnOuB95tOFKqmyYSv896GJAFfBc4Ba9bo7RU/MxOJi1+L81sAl5XcF5RUYTdKnReQQd9AXBmie0WwI6A6lJlZlYLL+TznHMz/OJdRf999B93B1W/SuoDDDWzfLyutQF4LfyGfrcApO7nVgAUOOfe9bdfxAv+VP/MBgKbnXN7nHNHgRnABaTHZwZlfz4pnydmNha4DMhx38x9r/J5BR307wNt/dkA38IbcJgVcJ0qxe+3fgJY7Zx7uMRLs4Cx/vOxwMuJrltVOOfucM61cM5l4n0+bznncoAFwFX+bil3XgDOuU+BbWbW3i+6EFhFin9meF02vc2snv97WXReKf+Z+cr6fGYB1/mzb3oDB4q6eFKBmV0M3AYMdc4VlnhpFnC1mdU2syy8web3KnRw51ygX8CleCPMG4EJQdenCufRF++/Ux8BK/yvS/H6s+cD6/3Hbwdd1yqcY39gtv+8tf/LtgF4AagddP0qeU7dgaX+5zYTOD0dPjPgbmANsBJ4Bqidip8Z8BzeOMNRvJbtjWV9PnhdHI/6WfIx3qyjwM+hAue1Aa8vvig//lJi/wn+ea0FLqno++nKWBGRNBd0142IiMSZgl5EJM0p6EVE0pyCXkQkzSnoRUTSnIJeRCTNKehFRNKcgl5EJM39H2lRH9bBBm9xAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Answer 1.c\n",
    "plt.plot(x,y,'o',color='blue',label='data')\n",
    "plt.plot(x,b0+b1[0]*x,'--',color='orange',label='regression line (unobsrvd)')\n",
    "# UNCOMMENT THE LINE BELOW!\n",
    "#plt.plot(x,b0_est+b1_est*x,color='green',label='least-squares')\n",
    "plt.legend(loc=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Part B. (Accuracy)** \n",
    "You can now evaluate the accuracy of your prediction, using the following metrics:\n",
    "\n",
    "* d) Evaluate the Residual Standard Error (RSE) and the Mean Squared Error (MSE). \n",
    "What do they tell about the accuracy of the model?\n",
    "\n",
    "**Note:** The Mean Squared Error (MSE) can be calculated as $MSE := \\frac{1}{n}RSS = \\frac{n-2}{n}RSE^2$.\n",
    "\n",
    "* e) Evaluate the $R^2$-statistic. \n",
    "What does it tell about the accuracy of the model?\n",
    "* f) Determine the $95\\%$ confidence intervals for the coefficients. \n",
    "Do the real coefficients $\\beta_0$ and $\\beta_1$ fall in the calculated intervals?\n",
    "\n",
    "(again) You can compare your results with the ones from the **Simple Linear Regression using scikit-learn** that follows."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer 1.B** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple Linear Regression using scikit-learn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to import appropriate libraries from sklearn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets, linear_model\n",
    "from sklearn.metrics import mean_squared_error, r2_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = x\n",
    "#print(x)\n",
    "# if x only 1 dimension then we need to add axis\n",
    "if len(np.shape(x1))==1:\n",
    "    x_sc = x1[:,np.newaxis]\n",
    "#print(x_sc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create linear regression object\n",
    "regr = linear_model.LinearRegression()\n",
    "# Train the model using the training sets\n",
    "regr.fit(x_sc,y)\n",
    "# Make predictions using the same x set\n",
    "y_pred = regr.predict(x_sc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Slope: 3.15\n",
      "Intercept: -1.05\n"
     ]
    }
   ],
   "source": [
    "# The coefficient $\\beta_1$\n",
    "print('Slope: %.2f' % regr.coef_)\n",
    "# The intercept $\\beta_0$\n",
    "print('Intercept: %.2f' % regr.intercept_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean_squared_error: 1441.63\n",
      "R2_score: 0.89\n"
     ]
    }
   ],
   "source": [
    "# The mean squared error\n",
    "print(\"Mean_squared_error: %.2f\"\n",
    "      % mean_squared_error(y, y_pred))\n",
    "# Explained variance score: 1 is perfect prediction\n",
    "r2n = r2_score(y, y_pred)\n",
    "print('R2_score: %.2f' % r2n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAveElEQVR4nO3deXxU1f3/8ddJCIQQkF2BEELZIUCAsKksiiyKLAURMCBoC1VU3AVBRdD8tGJLtVIttYLVACIiYuUrggVxATRYtSyyyRZANgEJJJDl8/vjTiYzk8k+ez7PxyOPyT1zZ+acEN5czjn3HCMiKKWUCi1h/q6AUkopz9NwV0qpEKThrpRSIUjDXSmlQpCGu1JKhaBK/q4AQN26dSUuLs7f1VBKqaCydevWUyJSz91zARHucXFxpKam+rsaSikVVIwxBwt7rsTdMsaYcGPMf40x/7Yd1zbGrDXG7LE91nI493FjzF5jzC5jzMDyVV8ppVRplabP/X5gp8PxdOBTEWkBfGo7xhjTFhgDtAMGAX8zxoR7prpKKaVKokThboyJAQYDrzsUDwPetH3/JjDcoXypiFwSkf3AXqCbR2qrlFKqREra5/4X4DGgukPZlSJyDEBEjhlj6tvKGwGbHc5Ls5U5McZMBiYDxMbGFvjArKws0tLSyMzMLGEVlfK+yMhIYmJiiIiI8HdVlCpSseFujLkZOCEiW40xfUvwnsZNWYEFbERkAbAAIDExscDzaWlpVK9enbi4OIxx95ZK+ZaIcPr0adLS0mjatKm/q6NUkUrSLXMNMNQYcwBYClxvjHkbOG6MaQBgezxhOz8NaOzw+hjgaGkrlpmZSZ06dTTYVcAwxlCnTh3936RykpICcXEQFmY9pqT4u0aWYsNdRB4XkRgRicMaKP2PiIwDVgETbKdNAD6wfb8KGGOMqWKMaQq0AL4uS+U02FWg0d9J5SglBSZPhoMHQcR6nDw5MAK+PHeoPg/0N8bsAfrbjhGR7cAyYAfwMXCPiOSUt6JKKRVoZs6Eixedyy5etMr9rVQ3MYnIBmCD7fvTQL9CzksGkstZt4Dz9NNPEx0dzSOPPOL2+ZUrV9KyZUvatm3r45oppfzh0KHSlfuSri3jQStXrmTHjh3+roZSykfcTPQrstyXQibcvTWokZycTKtWrbjhhhvYtWsXAP/4xz/o2rUrHTt2ZOTIkVy8eJGvvvqKVatW8eijj5KQkMC+ffvcnqeUCh3JyRAV5VwWFWWV+52I+P2rS5cu4mrHjh0Fygrz9tsiUVEi1pCG9RUVZZWXR2pqqsTHx8uFCxfk3Llz0qxZM5k7d66cOnXKfs7MmTPl5ZdfFhGRCRMmyLvvvmt/rrDzVHArze+mCn1vvy3SpImIMdZjeXOnNIBUKSRXA2LhsPIqalAjKans7/v555/z29/+lijbP81Dhw4FYNu2bTzxxBOcPXuW9PR0Bg50v3xOSc9TSgWvpKTy5Yy3hES3jDcHNdxNfZs4cSKvvPIK//vf/5g1a1ah855Lep5SqnQCdW55IAmJcPfWoEbv3r15//33ycjI4Pz583z44YcAnD9/ngYNGpCVlUWKw29V9erVOX/+vP24sPOUUmUXyHPLA0lIhLu3BjU6d+7M6NGjSUhIYOTIkfTq1QuAZ555hu7du9O/f39at25tP3/MmDHMnTuXTp06sW/fvkLPU0qVXSDPLQ8kxuqT96/ExERx3axj586dtGnTpsTvkZJi/eEeOmRdsScnB2Y/mAp+pf3dVJ4VFmZdsbsyBnJzfV+fsjp49iBxL8Wx856dtK5btos/Y8xWEUl091xIDKhC4A5qKKU8KzbW6opxVx4Mfsn4hTov1LEfVwmv4pXPCYluGaVUxRHQc8uLkJmdiZlt7MHerFYzZJbQtJZ3VhgNmSt3pVTFkPc/9GDphs2VXMLnOG9GJ7O83x2u4a6UCjrB0A0rIoTNce4cyXkqhzDjmw4TDXellPKw2HmxHP71sP04c2YmVSp5p2+9MBruSinlIQPeGsDan9baj89MO0PNyJqFvyDzFETUgPDKHq+LDqgGkaNHj3LLLbd45L369u1L3vTTm266ibNnz3rkfVeuXMmcOXM88l6O4uLiOHXqlMff15Hjz+SGG27gzJkzXv08FTpi58ViZht7sB964BAySwoP9vSfYLGBFfXg5BdeqZOGewmICLnlnECbnZ1d7no0bNiQ5cuXl/t9XK1evZqaNWt65L1eeOEFpkyZ4pH38oSy/tzHjx/P3/72Nw/XRoWaYUuHYWYbexfM93d9j8wSGl/R2P0Ljqy2Qv3LsdZxl5fgquu9UjcN90IcOHCANm3aMGXKFDp37szhw4eZO3cuXbt2pUOHDsyaNct+7jPPPEPr1q3p378/Y8eO5cUXXwSsK8EZM2bQp08fXnrpJbZu3UqfPn3o0qULAwcO5NixYwC8/PLLtG3blg4dOjBmzBgAPvvsMxISEkhISKBTp06cP3+eAwcOEB8fD1h7zN5xxx20b9+eTp06sX79egAWLVrEiBEjGDRoEC1atOCxxx4rtq15V8V5bZ40aRLt2rVjwIABZGRkALBv3z4GDRpEly5d6NWrFz/++GOB99m9ezdVqlShbt26gLW2juM/RtHR0QBs2LCBvn37csstt9C6dWuSkpLIu5nu008/pVOnTrRv354777yTS5cu2V8/d+5cunXrRrdu3di7dy8A7777LvHx8XTs2JHevXvbfwajRo1iyJAhDBgwgNGjR7N69Wr7+0ycOJH33nuPjIwMxowZQ4cOHRg9erS9rWAtErdkyZJif3aqYpq+bjpmtmHVrlUAvD7kdWSW0OHKDu5fsD/FCvXPBlvHVy8mRYS4gVO9tz5OYctF+vKrvEv+esP+/fvFGCObNm0SEZE1a9bIpEmTJDc3V3JycmTw4MHy2WefyTfffCMdO3aUixcvyq+//irNmzeXuXPniohInz595O677xYRkcuXL0vPnj3lxIkTIiKydOlSueOOO0REpEGDBpKZmSkiImfOnBERkZtvvlm++OILERE5f/68ZGVlyf79+6Vdu3YiIvLiiy/KxIkTRURk586d0rhxY8nIyJCFCxdK06ZN5ezZs5KRkSGxsbFy6NChAu3r06ePfPPNNyIi0qRJEzl58qTs379fwsPD5b///a+IiIwaNUreeustERG5/vrrZffu3SIisnnzZrnuuusKvOcbb7whDz30kP3YdQnkatWqiYjI+vXrpUaNGnL48GHJycmRHj16yOeffy4ZGRkSExMju3btEhGR8ePHy7x58+x1fPbZZ0VE5M0335TBgweLiEh8fLykpaU5/ewWLlwojRo1ktOnT4uIyIoVK+T2228XEZFLly5JTEyMXLx4Uf70pz/Z/wy+//57CQ8Pt/9MRESaN2/utGxzHn//bir/eX3r68LT2L+mr51e9AtOfS3y/SyRFKyv3a+JiOeWKScklvz94WnYNjv/eJBtuYKPHe68jZ8FHZ6G9xtChnVVTK3OcONW2DIZ9v0j/9zhRyCqYZEf2aRJE3r06AHAJ598wieffEKnTp0ASE9PZ8+ePZw/f55hw4ZRtWpVAIYMGeL0HqNHjwZg165dbNu2jf79+wOQk5NDgwYNAOjQoQNJSUkMHz6c4cOHA3DNNdfw0EMPkZSUxIgRI4iJiXF63y+++IL77rsPgNatW9OkSRN2794NQL9+/bjiiisAaNu2LQcPHqRx40L+m+iiadOmJCQkANClSxcOHDhAeno6X331FaNGjbKf53hFnefYsWPUq1evRJ/TrVs3e5sSEhI4cOAA1atXp2nTprRs2RKACRMmMH/+fB544AEAxo4da3988MEHAevnNHHiRG699VZGjBhhf//+/ftTu3ZtAG688UamTp3KpUuX+Pjjj+nduzdVq1Zl48aNTJ06FbD+DDp0cL7qql+/PkePHqVOnTqoim3dT+vo/1Z/+/GQlkNYNXZV4S/YuwC+/oP1/VUDYHgaRDWyP+2tZcodBU+4d3ja+nJ1m5ubAX57tGBZ9wXWVylUq1bN/r2I8Pjjj/OHP/zB6Zx58+aV6D1EhHbt2rFp06YC53z00Uds3LiRVatW8cwzz7B9+3amT5/O4MGDWb16NT169GDdunVERkY61acwVarkT7kKDw8vVb+z62szMjLIzc2lZs2afPfdd0W+tmrVqpw7d85+XKlSJftYhYhw+fLlIutYVJvAefnlvO9fe+01tmzZwkcffURCQoK9jo5/dpGRkfTt25c1a9bwzjvv2P+RcH1PV5mZmfZ/tFXFtO3ENtq/2t5+HHtFLAcfcLP2QZ7ME7Diyvzjnm9B03EFTvPF3qva515CAwcO5I033iA9PR2AI0eOcOLECa699lo+/PBDMjMzSU9P56OPPnL7+latWnHy5El7uGdlZbF9+3Zyc3M5fPgw1113HS+88IJ9Y499+/bRvn17pk2bRmJiYoE+7t69e9uXEd69ezeHDh2iVatWXml7jRo1aNq0Ke+++y5gBfX3339f4Lw2bdrY+8LB6svfunUrAB988AFZWVlFfk7r1q05cOCA/T3eeust+vTpY3/+nXfesT/27NkTsMYCunfvzpw5c6hbty6HDx8u+MZYK3YuXLiQzz//3L5piuPPcNu2bfzwww/280WEn3/+mbi4uCLrrELTz+k/Y2Ybp2CXWVJ4sG+aaJv9ciU0GgLXfWJdeLoJdvDN3qvBc+XuZwMGDGDnzp32UImOjubtt9+ma9euDB06lI4dO9KkSRMSExPtXSKOKleuzPLly5k6dSrnzp0jOzubBx54gJYtWzJu3DjOnTuHiPDggw9Ss2ZNnnzySdavX094eDht27blxhtvtA/AAkyZMoW77rqL9u3bU6lSJRYtWuR0NexpKSkp3H333Tz77LNkZWUxZswYOnbs6HRO7969efjhhxERjDFMmjSJYcOG0a1bN/r16+d0Ne1OZGQkCxcuZNSoUWRnZ9O1a1fuuusu+/OXLl2ie/fu5Obm2gc7H330Ufbs2YOI0K9fPzp27Oj2fxgDBgzg9ttvZ+jQoVSubM0pvvvuu7njjjvo0KEDCQkJdOvWzX7+1q1b6dGjB5Uq6V+RiuTC5QtEPxftVFbkUgEnv4S11+Yf3/AZ1O9d7OckJ1tr0Dt2zXh8fZzCOuN9+RWIA6qlcf78eRERuXDhgnTp0kW2bt3q5xr5z9SpU2Xt2rX+rka5TZ06VdatW+f2uWD63VQlk52T7TRQytMU/YLUB/MHSVMQ2f5CqT/TE3uvEhIDqgFs8uTJ7Nixg8zMTCZMmEDnzp39XSW/mTFjBlu2bPF3NcotPj6efv36+bsaygfMbOdxl6wns6gUVkg0pk6F3X/NPx52EKqVrS/F2+vjaLh7wOLFi/1dhYBx5ZVX2jcSD2aTJk3ydxWUl7mGepFLBRxYDF85JHGHZyD+iXJ9vrc3GNJwV0pVKK6hvve+vTSr3cz9yetvhGMf5x8P2QPVm5e7Dnn7wOb1ueftAwsVcSqkUkqVg2uofzbxM3o3KWTw88RGWJc/U4tuf4fmkz1WF53nrpRS5dToz404ej7/3peFwxYyMWGi+5PX9bGCPU+PN+E3t3u8Tr6Y567hrpQKSSOXjWTFzhX240d6PsLcAXPdn/zZUDjyYf5x3/+DhoO8Vjdf7AOrNzEVIW+hK09ZtGgRR4+6uXtWKeUxcz6bg5lt7MHer2k/ZJa4D/b1N1k3H+UFe7e/WzcfeTHYwTf7wOqVuw8tWrSI+Ph4GjYsek0bb8ubBxsWpv+2q9Dx6jevMmV1/nLTURFRXJhxwf3JSyMh12F9pEFbobbvpjD7Yh9Y/dtdQoUt9zt8+HC6dOlCu3btWLDAWrsmJyeHiRMnEh8fT/v27Zk3bx7Lly8nNTWVpKQkEhISnJaXBWvRrd69e5OQkEB8fDyff/45AAsXLqRly5b06dOHSZMmce+99wKFL6ebnp5Ov3796Ny5M+3bt+eDDz4ASr6E8YULFxg8eDAdO3YkPj7efsu/UoHq058+xcw2TsEus8R9sO993bpSzwv2Ds9YV+o+DPY8SUlw4ADk5lqPHp/zXtjdTb78CtQ7VPOWqC1suV8RsS8re/HiRWnXrp2cOnVKUlNT5YYbbrC/T95StI7L7Lp68cUX7UvaZmdny6+//ipHjx6Vxo0by4kTJ+TSpUty9dVXyz333CMihS+nm5WVJefOnRMRkZMnT0qzZs0kNze3xEsYL1++XH7/+9/b3/fs2bPl+yGGoED43VQie07vKfldpY53k6Yg8vN/nJ72xN2i/kAo3KHqOo3JE4pcM8JBYcv99u7dm5dffpn3338fgMOHD7Nnzx5atWrFTz/9xH333cfgwYMZMGBAsZ/RtWtX7rzzTrKyshg+fDgJCQl8+umn9O3b176M7ujRo+3L+hbaJhFmzJjBxo0bCQsL48iRIxw/fhwo2RLGvXr14pFHHmHatGncfPPN9OrVq0Q/I6V85UzGGWq/UNuprNC/yz/Mgm0O2z5evRjixjqd4os55/4QNOFe0iD2ymcXstzvhg0bWLduHZs2bSIqKoq+ffuSmZlJrVq1+P7771mzZg3z589n2bJlvPHGG06v3bJli/395syZw9ChQ9m4cSMfffQR48eP59FHH6VGjRqFLklb2HK6KSkpnDx5kq1btxIREUFcXByZmZlAyZYwBmvRrNWrV/P4448zYMAAnnrqqTL+5JTynOzcbCKeiXAqy30q1/3fkaWVIddhFdJuC6C5+7uOfTHn3B+CJtz9aeDAgTz55JMkJSURHR3NkSNHiIiI4Ny5c9SqVYuoqCh+/PFHNm/eDMCpU6eoXLkyI0eOpFmzZkycOBGA6tWrc/78eQC6d+/utHrhwYMHadSoEZMmTeLChQt8++23TJs2jfvvv5/Tp09To0YN3n33XftKjHnL6d56661Oy+meO3eO+vXrExERwfr16znobr5VEW3Kzs6mdu3ajBs3jujoaBYtWuSdH6pSpeD6P/eMmRlEVooseOLiMMDhQnDg11Cna5Hv7Ys55/5QbLgbYyKBjUAV2/nLRWSWMaY28A4QBxwAbhWRM7bXPA78DsgBporIGq/U3kcKW+530KBBvPbaa3To0IFWrVrZuzyOHDnCHXfcYb+yfu655wBrEPSuu+6iatWqbNq0yWkjiA0bNjB37lwiIiKIjo7mX//6Fw0aNODpp5+mZ8+eNGjQgM6dO5OTkwNQ6HK6SUlJDBkyhMTERBISEmjdunWp2rR3714effRRwsLCiIiI4NVXX/XCT1SpknEN9cMPHiamRkzBE1fGwkWHtfwTX4GW95ToM3wx59wfjBS/+40BqolIujEmAvgCuB8YAfwiIs8bY6YDtURkmjGmLbAE6AY0BNYBLUUkp7DPSExMlNTUVKeynTt30qZNm3I0LfQsWrSI1NRUXnnlFX9XpULT303vcw31L+/8kqsbX13wxMUuXTL9v4J6PUv1Wa597mDNOV+wIPC7ZYwxW0Uk0d1zxV6520Zk022HEbYvAYYBfW3lbwIbgGm28qUicgnYb4zZixX0BfeXU0opB66h/o8h/+D3nX9f8ETXeepdXoJWU8v0mb6Yc+4PJepzN8aEA1uB5sB8EdlijLlSRI4BiMgxY0x92+mNgM0OL0+zlbm+52RgMkBssP//x0cmTpxo779XKpS4hvrkzpP5+5C/FzzR9UrdQ8sEeHttdX8oUbjbulQSjDE1gfeNMfFFnO5uekeBvh8RWQAsAKtbppDPLXIDY6V8rbhuTFU6kc9Gcikn/wq8Wa1m7J26t+CJrqF+7TKIHeXl2gW3Us2WEZGzxpgNwCDguDGmge2qvQFwwnZaGtDY4WUxQKkXVImMjOT06dPUqVNHA14FBBHh9OnTREa6maWhSmXokqF8uPtDpzK3051dQ93DS++GspLMlqkHZNmCvSpwA/BHYBUwAXje9viB7SWrgMXGmD9jDai2AL4ubcViYmJIS0vj5MmTpX2pUl4TGRlJTIyb2RqqRJ77/Dlm/GeGU1mJQr3/l1DPzYCqKlRJrtwbAG/a+t3DgGUi8m9jzCZgmTHmd8AhYBSAiGw3xiwDdgDZwD1FzZQpTEREBE2bNi3ty5RSAWj1ntUMXjzYqaxEoZ44H1pOKXieKlaxUyF9wd1USKVU8Nt1ahet5zvfa1GiUL9+LVx1gxdrFhrKNRVSKaVK61zmOWr+saZTWYlCvdd70HiE9ypWgWi4K6U8JldyCZ8T7lzmbv2XArNflkPsSC/XrmLRcFdKeYTrXPX0x9OpVrma80muoT4oFWp38XLNKiYNd6VUubiG+oH7D9CkZhPnk7x085EqnO7EpFSISUmBuDgIC7MeU1K88zlmtnEK9vUT1iOzxDnYFxvnYB+80yd7lCq9clcqpPhi4wnXK/W/3vhX7u12r/NJeqXudzoVUqkQEhfnfvnaJk2sfTrLwzXUx3UYx1u/fcv5JNdQv+kHqNm+fB+sCqVTIZWqILyx8USdF+rwS8Yv9uOG1Rty5KEjzie5hvqALVC3W9k/VJWbhrtSIcSTG0+MXj6aZduXOZUVmKteINQ3Qd0epf8w5XE6oKpUgCvNAGlysrXRhKOoKKu8pOZtmoeZbZyCXWaJc7C7DpQO3WcNlGqwBwy9clcqgJV2gLQ8G0+s+2kd/d/q71RW/JX6ZqjbvQQtUb6mA6pKBTBvDpDm2ffLPpr/tblTWbGhPvwwROnqmP6mA6pKBSlvDJDmSb+cTvXnqjuVFRvqN/8INVqV/8OV12mfu1IBrLCB0PLsTJkruZjZxinYc57KyQ92kYJ96sMOWX3qARrsvrpxK5jolbtSASw52bnPHUo/QOrIda76r9N/pXoVW8iLwBKX671hB6FaYO9x7Isbt4KR9rkrFeBSUso2QOrINdT33LeH5rVt/ezuQn3oPoj+TTlq7Tu+GJcIVEX1uWu4KxXCXEN9zbg1DGg2wDpwF+ojTkBkPR/VzjPCwqymuDIGcnN9Xx9f0gFVpSoY11B/sf+LPHz1w9aB5MIS5zXXGbofouN8UzkP8+SNW6FEw12pEOIa6iPajOC9W9+zDtyF+shTUKWOj2rnHZ4elwgVGu5KhYBOf+/Edz9/Zz9uWrMpP93/k3XgLtRHHIfI+r6roBeV58atUKZTIZUqo0CYfnff6vsws41TsMsssYI9N9uazugY7CNPWVMaQyTY8yQlWYOnubnWY0UPdtArd6XKxN/T7xb+dyF3rrrTqcw+Tz03B5a6/NUe9StEON+wpEKbhrtSZTBzpnMfL1jHEyZY33sr4L86/BXXvHGNU1l+qGfB0srOL7jlDFSu6Z3KqICm4a5UGRR2+39Ojneu4A+ePUjcS3FOZUWG+ugMCI/0XAVU0NFwV6oMCpt+B9YV/MyZngn3C5cvEP1ctFOZPdRzLsM7VZxfMOocRNQo/weroKfhrlQZuJt+56i8C3uJCGFznOc7ZD+ZTXhYOGRnwDKXRdv1Sl250HBXqgzyrsonTLC6YlyV5wYa17nqvzz2C7Wq1oKcS7DY5a+shroqhIa7UmWUF/CeuoHGNdS33b2NdvXbweVzBZfe1VBXxdBwV6ocPHEDjWuov3fre4xoMwKyLxYM9TGXISyinLVWFYGGu1LllJRUtsFT11Cfce0Mkvslw+WzeqWuyk3DXQUETyxrGyxcQ/3a2Gv5/I7P3Xe/jMmCMP1rqkpPf2uU3/n7bk9fuerFqzh+4bhTmcwS91fq2v2iyknXc1d+F+qbLYxbMY6U/zkvPCOzBDJPwIornU8emwNGl3xSJaPruauA5s1NoP3pb9/8jXtW3+NUZl2pn3FzpZ4NYS4rNypVDhruyu9CbbOFjQc30mdRH6cymSWQvr9AqDedmcOzz4aRpBfrysM03JXfhcpmC4fPHSb2L87/IlndL6cKhLpJygWsslAcX1D+V+z1gjGmsTFmvTFmpzFmuzHmflt5bWPMWmPMHttjLYfXPG6M2WuM2WWMGejNBqjgl5QECxZYfezGWI8LFgRP2GVmZ2JmG6dgl1mCPGy7Ul+Rvydp3IxcTJKQF+yQvxaNUp5U7ICqMaYB0EBEvjXGVAe2AsOBicAvIvK8MWY6UEtEphlj2gJLgG5AQ2Ad0FJE3NykbdEBVRWM3K3/kvVkFpUyjsAHcc4nj80FYyr0Zs7K88o1oCoix4Bjtu/PG2N2Ao2AYUBf22lvAhuAabbypSJyCdhvjNmLFfSbytcMpQKH61z1k4+epK5chKUu0xdtoZ4n1MYXVOAqVZ+7MSYO6ARsAa60BT8icswYk7dvVyNgs8PL0mxlru81GZgMEKu/2SpIuIZ66qRUulQNh5X1nE90CfU8oTK+oAJficPdGBMNvAc8ICK/Gje/uHmnuikr8B9REVkALACrW6ak9VDKH1xDfdGwRUxo2h0+auN84m1F/yrrZs7KV0oU7saYCKxgTxGRFbbi48aYBrar9gbACVt5GtDY4eUxwFFPVVgpX3IN9cmdJ/P3a+6CjzvDDocnigl1R2Vdi0ap0ig23I11if5PYKeI/NnhqVXABOB52+MHDuWLjTF/xhpQbQF87clKK+VtrqEeVzOO/Xeugw+bw8cL8p8oRagr5UsluXXiGmA8cL0x5jvb101Yod7fGLMH6G87RkS2A8uwrms+Bu4paqaMUkVJSbGWJwgLsx5TUop7RfnUfL5mgWCX+1LZX++AFex5bhMNdhXQdG0ZFbBcFxQDa/DRG3Pgx743lqXbljqVyZSvYO3VzidqoKsAUtRUSA13FbB8saCY2/Vf7t4I63o7n6ihrgKQLhymgpI3FxT74tAX9FrYy6lM7lgOX9ziHOwa6ipIabirgOWNG36O/HqEmHkxTmUy5UtYe40V7Hk01FWQ07XoVMBKTrb62B2V9YafS9mXMLONU7DLpNVIC6xgz6MDpSpE6JW7ClieuuHHdfbLpUlrqbyhP2y4Kb9QA12FGB1QVSHLNdSPjltCgy1jnU/SUFdBTAdUVYXiGupfDnqKq/fNAcdg11BXIU7DXYUM11Bf0m8mYw4lw745+YUa6qqC0HBXQS0lBcbtdQ71R9sM4IXsT+CQw8irhrqqYDTcVdCqlxzHqez8uZI1fmnFue67IPuT/JM01FUFpeGugs7498fz9g9v248jM2uR0f4MsCv/JA11VcFpuKug8ZfNf+HBNQ86lUkLgDP2Y5Mk1pZ1t/m2bkoFGg13FfDW7F3DoJRBTmVWqOezNp226MZeSmm4qwC269QuWs9v7VTmGOpZEk3Nyed1yzql3NBwVwHnTMYZar9Q26nM9Uqd24QIrH0adcs6pQrScFcBIzs3m4hnIpzKcps77DNdoxXc/KPT87plnVLuabirgFBg/ZfmUNmxSGe/KFUqGu7Kr1xD/eRvoG64Q4GGulJlouGu/MI11Pc0geaVbQeR9WHEcd9XSqkQouGufMo11Dc0gj6Oa7brlbpSHqHhrnwidl4sh389bD/+Z3248wrbQbNJ0H2BfyqmVIjScFdeNerdUSzfsdx+/HBNeLGe7SC6OQzd45d6KRXqNNyVVzy78VmeXP+k/bhvVVift8NdnR4wcJN/KqZUBaHhrjxq+Y7ljHp3lP24ioHM5raDmGHQe6Vf6qVURaPhrjwi9WgqXf/R1anMfldpjdZw807fV0qpCkzDXZVL2q9pNJ7X2KnMHuot74PEl31fKaUUYf6ugApOFy5fwMw2TsEuLWzBXrOjNaXRz8GekgJxcRAWZj2mpPi1Okr5lF65q1LJlVzC54Q7l+Wt/xL/FHSY7Z+KuUhJgcmTsa8YefCgdQy6Fo2qGIyI/28aSUxMlNTUVH9XQxWj0PVfGg6Gvv/2T6UKERdnBbqrJk3gwAFf10Yp7zDGbBWRRHfP6ZW7KpZrqJ/9DVwRDsT8Fnqv8E+linHoUOnKlQo12ueuCmVmG6dgPxhn9am/smYuYeOEuNtXBGw/dmG7MekuTaqi0HBXBXR8raNTqKc2tkL9fKVpVJskPJHyCCL5/diBGPDJydauTI50lyZVkWi4K7vx74/HzDb8cPwHAD5oYIV6l/4r4DZh8Mznnba0A2vAcuZMP1S2GElJsGCB1cdujPW4YIEOpqqKQwdUFckbk3li/RP247/UhftrAYnzoeUUe3lYGLj7dTEGcnN9UFGllBMdUFVuLd22lLHvjbUfT7kC5tcHrlkKTUYXOD821v0MFO3HVirwaLdMCCru5p0vDn2BmW3swd6nqtX9Mn/AC9bNR26CHbQfW6lgUuyVuzHmDeBm4ISIxNvKagPvAHHAAeBWETlje+5x4HdADjBVRNZ4pebKraJu3uk2aA8tX2lpP7d6GPzaDOj8Z2j9YLHvnddfPXOmNaUwNtYKdu3HVirwFNvnbozpDaQD/3II9xeAX0TkeWPMdKCWiEwzxrQFlgDdgIbAOqCliOQU9Rna5+45bm/eqXoaptV1KpIWFNr9opQKDuXqcxeRjcaYOJfiYUBf2/dvAhuAabbypSJyCdhvjNmLFfS6eLePON2kE34Jnox0el5aAJ3/Aq3v92W1lFI+VtYB1StF5BiAiBwzxtS3lTcCNjucl2YrK8AYMxmYDBCrI3IeYw16CjztPJyS2xxMv3VwVT8/1Uwp5UueHlA1bsrc9vuIyAIRSRSRxHr16rk7RZXB0TsjnYL9UnOYsGExixENdqUqkLKG+3FjTAMA2+MJW3ka4Li4dwxwtOzVUyXV6e+dMLMNWXIJgPRmMGbNElrOFPr/fqwOeipVwZS1W2YVMAF43vb4gUP5YmPMn7EGVFsAX5e3kqpwo5ePZtn2Zfbj402h/tAfoGZ7lozzY8WUUn5VkqmQS7AGT+saY9KAWVihvswY8zvgEDAKQES2G2OWATuAbOCe4mbKqLKZ8ekMnvviOfvx7ibQ4qaNUL+XH2ullAoUJZktM7aQp9x24IpIMqC3tXjJq9+8ypTV+UsCbG4M3UftgerNi3iVUqqi0eUHgsQHP37A8HeG249XNYAht2yF2p39VymlVMDScA9wm9M20/OfPe3Hr9WHP4z+Fmp38mOtlFKBTsM9QO057bxUwIxakDxOQ10pVTIa7gHmePpxrvrTVfbjMdGwZPxnUL+3H2ullAo2Gu4BIv1yOtWfq24/7lIFUu/dD9Fx/quUUipoabj7WXZuNhHPRNiPqxq4OPV/UDPej7VSSgU7DXc/ERHC5ris//LwQUy0rrOjlCo/3azDQXGbXHhKzeeqOQX75Qf3ILNEg10p5TF65W5T1CYXnlqXpcdfG7PllzT78fkpqUTX6+KZN1dKKQe6QbaN200ugCZN4MCB8r33+EVdeftgfvuOTf6SqxpcXb43VUpVeLpBdgk4bXJRgvKSeOrdATyzY639+MfxK2n1m2Flf0OllCohDXcba5ML9+Wl9frSTkza9Z39+MvbVnJ1Cw11pZTv6ICqTXIyREU5l0VFWeUl9dGGKZjZxh7sK4a+hswSDXallM/plbtN3qDpzJlWV0xsrBXsJRlM/WbdaLp9mb+m+iv9n+Oeq6d7qaZKKVU8DXcHSUmlmxmz7/u5NF/5mP34se738cdBL3uhZkopVToa7mVw6tuZ1Pvw/9mPR7YYyPLbPvZjjZRSypmGeylc3HI31T5+zX7cvl4bfpiyw481Ukop9zTcSyDnm3uptHq+/TjchJP9VLYfa6SUUkXTcC+CrO5E2DffOZXlPpWLMcY/FVJKqRLScHdn55+o994jnHLY2vvSE5eoHF7Zf3VSSqlS0HB3tOsV/rj2Pqafzi86N/0cNarU8F+dlFKqDDTcAU5tZuWKnvz2WH7RkYeO0LB6Q//VSSmlyqFih/uxtXy2aTp9//ctALWr1uLHe3ZRr1o9P1dMKaXKp2KG+8FlfPef0XRyWBTswP0HaFKzif/qpJRSHlSxwv3kl+z9agotvvvBXrR9ynba1mvrx0oppZTnVYyFw7Y/z7F/GSq/eq092Df9bhMySwoNdl/tyqSUUt4Q2lfuP3/Kme+eosPWr0iz3XP0cdLHDGw+sMiX+WJXJqWU8qbQ3IlpWzIXD62k1/++5dvMXACWjFzCmPgxJXq5N3dlUkopT6k4OzEd/4ysdX0ZdhT+z3bV/dcb/8q93e4t1dt4Y1cmpZTypeAPdxH4eR253z7C7/btZdEZq/ip3k8x+7rZZXpLT+7KpJRS/hD04S6Lw5h+Gl6whfrdiXcz/6b55Vr/JTnZuc8dSr8rk1JK+VPQh3vYXutxRJsRvHPLO1QKK3+TyrMrk1JKBYKgD/clI5cwvPVwIitFevR9S7srk1JKBZKgD/eSzoBRSqmKpGLcxKSUUhWMhrtSSoUgr4W7MWaQMWaXMWavMWa6tz5HKaVUQV4Jd2NMODAfuBFoC4w1xujqXEop5SPeunLvBuwVkZ9E5DKwFBjmpc9SSinlwlvh3gg47HCcZiuzM8ZMNsakGmNST5486aVqKKVUxeStcHd3e6jTCmUiskBEEkUksV493flIKaU8yVvhngY0djiOAY566bOUUkq58Fa4fwO0MMY0NcZUBsYAq7z0WUoppVx45Q5VEck2xtwLrAHCgTdEZLs3PksppVRBXlt+QERWA6u99f5KKaUKp3eoKqVUCNJwV0qpEKThrpRSIUjDXSmlQpCGu1JKhSANd6WUCkEa7kopFYI03JVSKgRpuCulVAjScFdKqRCk4a6UUiFIw10ppUKQhrtSSoUgDXellApBQR3uKSkQFwdhYdZjSoq/a6SUUoHBa+u5e1tKCkyeDBcvWscHD1rHAElJ/quXUkoFgqC9cp85Mz/Y81y8aJUrpVRFF7ThfuhQ6cqVUqoiCdpwj40tXblSSlUkQRvuyckQFeVcFhVllSulVEUXtOGelAQLFkCTJmCM9bhggQ6mKqUUBPFsGbCCXMNcKaUKCtord6WUUoXTcFdKqRCk4a6UUiFIw10ppUKQhrtSSoUgIyL+rgPGmJPABeCUv+viBXXRdgUTbVfwCdW2laRdTUSknrsnAiLcAYwxqSKS6O96eJq2K7hou4JPqLatvO3SbhmllApBGu5KKRWCAincF/i7Al6i7Qou2q7gE6ptK1e7AqbPXSmllOcE0pW7UkopD9FwV0qpEOT3cDfGDDLG7DLG7DXGTPd3fcrKGNPYGLPeGLPTGLPdGHO/rby2MWatMWaP7bGWv+taFsaYcGPMf40x/7Ydh0q7ahpjlhtjfrT92fUMhbYZYx60/R5uM8YsMcZEBmO7jDFvGGNOGGO2OZQV2g5jzOO2LNlljBnon1oXr5B2zbX9Hv5gjHnfGFPT4blSt8uv4W6MCQfmAzcCbYGxxpi2/qxTOWQDD4tIG6AHcI+tLdOBT0WkBfCp7TgY3Q/sdDgOlXa9BHwsIq2BjlhtDOq2GWMaAVOBRBGJB8KBMQRnuxYBg1zK3LbD9vdtDNDO9pq/2TImEC2iYLvWAvEi0gHYDTwOZW+Xv6/cuwF7ReQnEbkMLAWG+blOZSIix0TkW9v357FCohFWe960nfYmMNwvFSwHY0wMMBh43aE4FNpVA+gN/BNARC6LyFlCoG1YezVUNcZUAqKAowRhu0RkI/CLS3Fh7RgGLBWRSyKyH9iLlTEBx127ROQTEcm2HW4GYmzfl6ld/g73RsBhh+M0W1lQM8bEAZ2ALcCVInIMrH8AgPp+rFpZ/QV4DMh1KAuFdv0GOAkstHU5vW6MqUaQt01EjgAvAoeAY8A5EfmEIG+Xg8LaEUp5cifwf7bvy9Quf4e7cVMW1HMzjTHRwHvAAyLyq7/rU17GmJuBEyKy1d918YJKQGfgVRHphLW+UTB0VRTJ1gc9DGgKNASqGWPG+bdWPhESeWKMmYnVzZuSV+TmtGLb5e9wTwMaOxzHYP33MSgZYyKwgj1FRFbYio8bYxrYnm8AnPBX/croGmCoMeYAVrfZ9caYtwn+doH1+5cmIltsx8uxwj7Y23YDsF9ETopIFrACuJrgb1eewtoR9HlijJkA3AwkSf5NSGVql7/D/RughTGmqTGmMtagwSo/16lMjDEGq+92p4j82eGpVcAE2/cTgA98XbfyEJHHRSRGROKw/nz+IyLjCPJ2AYjIz8BhY0wrW1E/YAfB37ZDQA9jTJTt97If1hhQsLcrT2HtWAWMMcZUMcY0BVoAX/uhfmVijBkETAOGishFh6fK1i4R8esXcBPWyPA+YKa/61OOdlyL9V+lH4DvbF83AXWwRvT32B5r+7uu5WhjX+Dftu9Dol1AApBq+3NbCdQKhbYBs4EfgW3AW0CVYGwXsARr3CAL6wr2d0W1A5hpy5JdwI3+rn8p27UXq289Lz9eK0+7dPkBpZQKQf7ullFKKeUFGu5KKRWCNNyVUioEabgrpVQI0nBXSqkQpOGulFIhSMNdKaVC0P8HRR0Nt++s6d0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot outputs\n",
    "f = plt.figure()\n",
    "plt.scatter(x, y,  color='blue',label='data')\n",
    "plt.plot(x, b0+b1[0]*x,'--', color='orange', linewidth=1, label='regression line (unobsrvd)')\n",
    "plt.plot(x, y_pred, color='green', linewidth=1, label='least-squares')\n",
    "\n",
    "plt.legend(loc=2)\n",
    "\n",
    "plt.show()\n",
    "# Export as pdf\n",
    "#f.savefig(\"LinearRegression.eps\", format='eps',bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple Linear Regression using statsmodels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** Statsmodels does not has as default the calculation of a constant bias.\n",
    "Uses OLS = Ordinary Least Squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   0.891</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.885</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   146.9</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Mon, 10 Oct 2022</td> <th>  Prob (F-statistic):</th> <td>4.29e-10</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:35:49</td>     <th>  Log-Likelihood:    </th> <td> -101.11</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    20</td>      <th>  AIC:               </th> <td>   206.2</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    18</td>      <th>  BIC:               </th> <td>   208.2</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>   -1.0524</td> <td>   19.328</td> <td>   -0.054</td> <td> 0.957</td> <td>  -41.658</td> <td>   39.553</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x</th>     <td>    3.1469</td> <td>    0.260</td> <td>   12.122</td> <td> 0.000</td> <td>    2.602</td> <td>    3.692</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 1.073</td> <th>  Durbin-Watson:     </th> <td>   2.484</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.585</td> <th>  Jarque-Bera (JB):  </th> <td>   0.845</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.199</td> <th>  Prob(JB):          </th> <td>   0.656</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.075</td> <th>  Cond. No.          </th> <td>    161.</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   R-squared:                       0.891\n",
       "Model:                            OLS   Adj. R-squared:                  0.885\n",
       "Method:                 Least Squares   F-statistic:                     146.9\n",
       "Date:                Mon, 10 Oct 2022   Prob (F-statistic):           4.29e-10\n",
       "Time:                        15:35:49   Log-Likelihood:                -101.11\n",
       "No. Observations:                  20   AIC:                             206.2\n",
       "Df Residuals:                      18   BIC:                             208.2\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const         -1.0524     19.328     -0.054      0.957     -41.658      39.553\n",
       "x              3.1469      0.260     12.122      0.000       2.602       3.692\n",
       "==============================================================================\n",
       "Omnibus:                        1.073   Durbin-Watson:                   2.484\n",
       "Prob(Omnibus):                  0.585   Jarque-Bera (JB):                0.845\n",
       "Skew:                          -0.199   Prob(JB):                        0.656\n",
       "Kurtosis:                       2.075   Cond. No.                         161.\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# using pandas\n",
    "X = df_slr[\"x\"]\n",
    "y = df_slr[\"y\"]\n",
    "\n",
    "X = sm.add_constant(X) ## let's add an intercept (beta_0) to our model\n",
    "\n",
    "# Note the difference in argument order\n",
    "model = sm.OLS(y, X).fit()\n",
    "y_pred_OLS = model.predict(X) # make the predictions using SAME train data\n",
    "\n",
    "# Print out the statistics\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** We observe the same estimates for intercept $\\hat{\\beta}_0$ named 'const coef', slope $\\hat{\\beta}_1$ named 'x coef', $R^2$-statistic named 'R-squared'."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following Dataset related to advertisement is given in pandas and numpy format:\n",
    "\n",
    "**Download the Advertising.csv dataset from the link:**\n",
    "\n",
    "https://www.statlearning.com/resources-first-edition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "directory = #YOUR DIRECTORY HERE\n",
    "prefix = \"Advertising.csv\"\n",
    "filename1 = directory+prefix\n",
    "dataAd = np.loadtxt(filename1, delimiter=\",\",skiprows=1,usecols=[1,2,3,4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pandas\n",
    "pdAd = pd.DataFrame(dataAd, columns=[\"TV\",\"radio\",\"newspaper\",\"sales\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "# numpy\n",
    "TV = pdAd.iloc[:,0].values\n",
    "radio = pdAd.iloc[:,1].values\n",
    "news = pdAd.iloc[:,2].values\n",
    "sales = pdAd.iloc[:,3].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "i) Perform simple linear regression between all pairs x-y, where y=sales. For this you need to use (a) scikit-learn tools, and (b) statsmodels. Compare the results and plot the regression lines and data points. Especially, check how much is the $R^2$ measure per feature, to understand how much of the output $y$ is explained per feature.\n",
    "\n",
    "ii) Compare the values of the slopes from the 3 simple linear regressions, with the slopes from a multi-linear regression that includes all 3 features. Find the $R^2$ from the multi-dimensional model, and compare it with the above values. What do you observe?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3 (Sample size & Error Variance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we will study the influence of sample size $n$ and the influence of error variance on the quality of the regression.\n",
    "\n",
    "i) If the initial value for the sample size was $n$, choose $n_{new} = 2n$, so that the new sample includes the $n$ initial samples and $n$ new ones. Set $n=20$ as before. What do you observe about the quality of the prediction? ($R^2$-statistic)\n",
    "\n",
    "ii) What happens when the Error standard deviation is doubled? tripled?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use the given function for synthetic data generator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** Normally, for small n_plus we need to average out many random realisations to see the improvement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4 (Unbiased Estimators)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use $L=15$ different samples of size $n=20$, with error $sigmae = 50$. Repeat the regression process and plot all least-squares curves, on the same plot. What do you observe? (for more emphasis, choose larger $sigmae$, e.g. double or triple). \n",
    "\n",
    "Take care to save the coefficients of each iteration for further use.\n",
    "\n",
    "With this, average over the $L$ values of each coefficient and draw the average curve. Compare this curve with the orginal regression line $\\beta_0+\\beta_1x$. What can you see? If you increase $L$ what happens to the average curve?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5 (Multiple Linear Regression)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now work on problems of multi-linear regression. The regression plane (>2 dimensions) is given by the following expression, for $k=4$ features:\n",
    "\n",
    "$y = \\beta_0 + \\beta_1x_1 + \\beta_2x_2 + \\beta_3x_3 + \\beta_4x_4 + \\epsilon$.\n",
    "\n",
    "First create synthetic data for $(x_1,x_2,x_3,y,Er)$. Let $x_4=np.ones(n)+Er$, hence a variable that does not influence $y$.\n",
    "\n",
    "* (a) For the above multi-regression model, estimate the unknown coefficients $(\\beta_0,...,\\beta_4)$. Determine the F-statistic, $R^2$ and the p-values for the four features.\n",
    "\n",
    "* (b) Use a multi-regression model only for $(x_2,x_4)$. Estimate the coefficients, and determine the F-statistic, $R^2$ and the p-values for the two features. How does $p-$ and $F-$ compare with (a)? What does this imply?\n",
    "\n",
    "* (c) Use $4$ single-regression models, one for each variable $x_i$. How do the estimated slopes $\\beta_1$, $\\ldots$, $\\beta_4$ compare with the estimated slopes you get using multi-regression in (a)?\n",
    "\n",
    "Compare the answers for cases (a), (b), (c). What do you observe? Explain.\n",
    "\n",
    "\n",
    "**Note 1** Use scikit-learn and create a linear regression object.\n",
    "\n",
    "**Note 2** For the p-values and the F-statistic, you can use the scikit-learn command\n",
    "sklearn.feature_selection.f_regression\n",
    "\n",
    "https://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.f_regression.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Synthetic data creation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100 # sample size\n",
    "b0 = 2 # real intercept (unobserved)\n",
    "b1 = np.array([3, 7, 20]) # real slope (unobserved)\n",
    "xl, xh = 0, 120\n",
    "sigmae = 30 # error mean and standard deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "(xM,yM,ErM) = genSample(n,b0,b1,sigmae,xl,xh,seedit=199,size=len(b1))\n",
    "x_extra = 1+ErM\n",
    "xM.append(x_extra)\n",
    "yM += xM[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data Visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataM = {'x1': xM[0], 'x2' : xM[1],'x3' : xM[2] ,'x4': xM[3], 'error': ErM, 'y_mult': yM}\n",
    "df_slrM = pd.DataFrame(data=dataM)\n",
    "#df_slrM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer**\n",
    "\n",
    "(Here I will just show the question (a) using statsmodels, you need to do it with scikit-learn.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>y_mult</td>      <th>  R-squared:         </th> <td>   1.000</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   1.000</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>4.563e+31</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 13 Oct 2021</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>12:19:59</td>     <th>  Log-Likelihood:    </th> <td>  2680.5</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   100</td>      <th>  AIC:               </th> <td>  -5351.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    95</td>      <th>  BIC:               </th> <td>  -5338.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     4</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    1.0000</td> <td> 1.78e-13</td> <td> 5.63e+12</td> <td> 0.000</td> <td>    1.000</td> <td>    1.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    3.0000</td> <td> 1.65e-15</td> <td> 1.82e+15</td> <td> 0.000</td> <td>    3.000</td> <td>    3.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>    <td>    7.0000</td> <td> 1.57e-15</td> <td> 4.47e+15</td> <td> 0.000</td> <td>    7.000</td> <td>    7.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x3</th>    <td>   20.0000</td> <td> 1.61e-15</td> <td> 1.24e+16</td> <td> 0.000</td> <td>   20.000</td> <td>   20.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x4</th>    <td>    2.0000</td> <td> 1.91e-15</td> <td> 1.05e+15</td> <td> 0.000</td> <td>    2.000</td> <td>    2.000</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 5.179</td> <th>  Durbin-Watson:     </th> <td>   0.472</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.075</td> <th>  Jarque-Bera (JB):  </th> <td>   2.835</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.172</td> <th>  Prob(JB):          </th> <td>   0.242</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.250</td> <th>  Cond. No.          </th> <td>    350.</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                 y_mult   R-squared:                       1.000\n",
       "Model:                            OLS   Adj. R-squared:                  1.000\n",
       "Method:                 Least Squares   F-statistic:                 4.563e+31\n",
       "Date:                Wed, 13 Oct 2021   Prob (F-statistic):               0.00\n",
       "Time:                        12:19:59   Log-Likelihood:                 2680.5\n",
       "No. Observations:                 100   AIC:                            -5351.\n",
       "Df Residuals:                      95   BIC:                            -5338.\n",
       "Df Model:                           4                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          1.0000   1.78e-13   5.63e+12      0.000       1.000       1.000\n",
       "x1             3.0000   1.65e-15   1.82e+15      0.000       3.000       3.000\n",
       "x2             7.0000   1.57e-15   4.47e+15      0.000       7.000       7.000\n",
       "x3            20.0000   1.61e-15   1.24e+16      0.000      20.000      20.000\n",
       "x4             2.0000   1.91e-15   1.05e+15      0.000       2.000       2.000\n",
       "==============================================================================\n",
       "Omnibus:                        5.179   Durbin-Watson:                   0.472\n",
       "Prob(Omnibus):                  0.075   Jarque-Bera (JB):                2.835\n",
       "Skew:                           0.172   Prob(JB):                        0.242\n",
       "Kurtosis:                       2.250   Cond. No.                         350.\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# using pandas\n",
    "Xm = df_slrM[[\"x1\", \"x2\", \"x3\", \"x4\"]]\n",
    "Ym = df_slrM[\"y_mult\"]\n",
    "Xm = sm.add_constant(Xm) ## let's add an intercept (beta_0) to our model\n",
    "\n",
    "# Note the difference in argument order\n",
    "model = sm.OLS(Ym, Xm).fit()\n",
    "y_pred_OLS = model.predict(Xm) # make the predictions using SAME train data\n",
    "\n",
    "# Print out the statistics\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do the same using scikit-learn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The END."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}