KG_Linear_Regression_HW.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 1\n",
    "Create any multiple linear model. For example, it may be $y = 3x_1 + 5x_2 + \\mu$ or $y = 4x_1 - 3x_2 + 9x_3 + \\mu$. Then:\n",
    "    - (a) generate random 30 observations depend on the model.\n",
    "    - (b) fit a linear model depends on the observations generated by step (a), and record the coefficients.\n",
    "    - (c) repeat (a) and (b) 100 times. Calculate the mean and variance of the coefficients.\n",
    "Note that each time the generated observations should come from the sample distribution. If you do not understand, just use the `random.rand` function. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta0 mean: 3.01364\n",
      "beta0 variance: 0.03500\n",
      "beta1 mean: 5.01981\n",
      "beta1 variance: 0.03289\n"
     ]
    }
   ],
   "source": [
    "### your solution\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn import linear_model\n",
    "import pandas as pd\n",
    "\n",
    "# setup\n",
    "np.random.seed(100)\n",
    "coefficients = []\n",
    "ols = linear_model.LinearRegression()\n",
    "\n",
    "for i in range(100):\n",
    "    x = np.random.rand(30,2) # 30 random numbers in 2 columns\n",
    "    y = 3*x[:,0] + 5*x[:,1] + np.random.rand(30) # formula 3X + 5X + u\n",
    "    ols.fit(x,y)\n",
    "    coefficients.append(ols.coef_)\n",
    "\n",
    "# compute mean and variance on coeffecients\n",
    "\n",
    "beta0_coef_results = pd.Series(np.array(coefficients)[:,0])\n",
    "beta1_coef_results = pd.Series(np.array(coefficients)[:,1])\n",
    "\n",
    "print(\"beta0 mean: %.5f\" % beta0_coef_results.mean())\n",
    "print(\"beta0 variance: %.5f\" % beta0_coef_results.var())\n",
    "print(\"beta1 mean: %.5f\" % beta1_coef_results.mean())\n",
    "print(\"beta1 variance: %.5f\" % beta1_coef_results.var())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 2\n",
    "This question involves the use of simple linear regression on the Auto data set(In the file \"Auto.csv\"). Perform a simple linear regression with *mpg* as the response and *horsepower* as the predictor. \n",
    "    - (a) What is the coefficient of mpg? How to explain it?\n",
    "    - (b) What is the coefficient of determination of the model? How to explain it?\n",
    "    - (c) Write down the formula of the model.\n",
    "    - (d) What is the predicted mpg associated with a horsepower of 98?\n",
    "    - (e) Plot the response and the predictor.  as well as the least squares regression line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A: ####################################\n",
      "Coef of mpg is: -0.1578447334 -- it is derived from the formula of coefficients. It is an optimization formula to minimize RSS. \n",
      "B: ####################################\n",
      "Coef determination of our model is: 0.6059482579 -- it is based the calculation of R squared https://en.wikipedia.org/wiki/Coefficient_of_determination\n",
      "C: ####################################\n",
      "The fomular is Y = -0.15784X + 39.93586\n",
      "D: ####################################\n",
      "The horsepower of 98 based on our model predicts mpg at: 24.47 \n",
      "E: ####################################\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAikAAAF3CAYAAABg/9sEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYFNXZNvD7mZ4ZyYCCDERBpBv3oKgR1LgGhRhlMBq3\naEaDaEKEmOCX1zijZNFEdNBoNL6CLxoMOpNdo4ZBRXHHLeC+IagN0SgKiLKowHC+P6qL6aX2rq2r\n79919TXTp6urTnX30A+nnuccUUqBiIiIKG5qou4AERERkREGKURERBRLDFKIiIgolhikEBERUSwx\nSCEiIqJYYpBCREREscQghYiIiGKJQQoRERHFEoMUIiIiiiUGKURERBRLtVF3wIl+/fqpTCYTdTeI\niIjIB4sWLVqplOpvt11FBCmZTAYLFy6MuhtERETkAxFZ5mQ7Xu4hIiKiWGKQQkRERLHEIIWIiIhi\niUEKERERxRKDFCIiIoolBilEREQUSwxSiIiIKJYYpBAREVEsMUghIiKiWGKQQkRERLHEIKUSdXQA\nmQxQU6P97OiIukdERES+q4i1eyhPRwcwYQKwYYN2f9ky7T4ANDdH1y8iIiKfcSSl0kyZ0h2g6DZs\n0NqJiIgShEFKpVm+3F07ERFRhWKQUmkGD3bXTkREVKEYpFSaqVOBhobCtoYGrZ2IiChBGKRUmuZm\nYOZMIJ0GRLSfM2cyaZaIiBKH1T2VqLmZQQkRESUeR1KIiIgolhikEBERUSwxSCEiIqJYYpBCRERE\nscQghYiIiGKJQUpScNFBIiJKGJYgJwEXHSQiogTiSEq+Sh2N4KKDRESUQBxJ0VXyaAQXHSQiogTi\nSIrOzWhE3EZcuOggERElEIMUndPRCH3EZdkyQKnuEZcoAxUuOkhERAnEIEXndDQijvkfXHSQiIgS\niEGKzuloRFzzP5qbgWwW2LJF+8kAhYiIKhyDFJ3T0QjmfxAREYWCQYpbzP8gIiIKBYMUndOEWOZ/\nEBERhUKUUlH3wdaIESPUwoULgz1IJqMFJsXSaS3Hg4iIiHwhIouUUiPstuNIis4oQLFqJyIiokAx\nSNGlUu7aiYiIKFAMUnRdXe7aiYiIKFAMUnTptLt2IiIiChSDFB1Li4mIiGKFQYqOpcVERESxUht1\nB2KluZlBCRERUUxwJIWIiIhiiUEKERERxRKDFCIiIoolBilkraNDWzKgpkb7WbyWERERUUCYOEvm\n9EUXN2zQ7uuLLgJMMCYiosBxJIXMTZnSHaDoNmzQ2omIiALGIIXMLV/urp2IiMhHDFLI3ODB7tqJ\niIh8xCClEoWVzMqlAoiIKEIMUiqNnsy6bBmgVHcyaxCBCpcKICKiCIlSKuo+2BoxYoRauHBh1N2I\nh0xGC0yKpdNANht2b4iIiFwTkUVKqRF223EkpdIwmZWIiKpE4EGKiKRE5HkRmZO731dEHhCRJbmf\n2wfdh0RhMisREVWJMEZSJgN4Pe9+K4D5SqndAczP3SenmMxKRERVItAgRUQGAWgCcEte8wkAZud+\nnw3gxCD7kDhMZiUioioR9LT41wG4CMC2eW07KKXez/3+AYAdAu5D8jQ3MyghIqLEC2wkRUTGAvhQ\nKbXIbBullRYZlheJyAQRWSgiCz/66KOguklEREQxFeTlnsMAfEtEsgD+AuBoEWkHsEJEBgBA7ueH\nRk9WSs1USo1QSo3o379/gN2sAlzJmIiIKlBgQYpS6mKl1CClVAbA6QAeUkqdCeAeAONym40DcHdQ\nfSCEO/kbERGRj6KYJ6UNwDdEZAmA0bn71SPsUQ2uZExERBUqlCBFKfWIUmps7vdVSqlRSqndlVKj\nlVKrw+hDLJQzquE1uOHkb0REVKE442yYvI5qlBPccPI3IiKqUAxSwuR1VKOcSzac/I2IiCoUg5Qw\neR3VKOeSDSd/IyKiCsUgJUxeRzXKvWTT3KytkLxli/aTAQoREVUABilh8jqqwUs2RERUhRikhM1o\nVMOucoeXbIiIqAoxSImaXeWOHsCcdZZ2//bbjS/ZTJoE1NZqQUxtrXY/7jgTLhERWQh6gUGyY1e5\nM2FC9+N6AAMUBimTJgEzZnTf7+rqvj99ejD9LpcenNmdGxERVS3R1viLtxEjRqiFCxdG3Y1g1NRo\nIyhupNPaaIqutlYLTIqlUsDmzWV1LzCZjBaYFCs+NyIiShwRWaSUGmG3HS/3RM3LpGrFpcdGAYpV\nexxwJlwiIrLBICVqRpU7dooDm1TKeDuz9jjgTLhERGSDQUrUiit37DQ0AGPGFCacjhxpvO3IkfFN\nTGVZNRER2WCQEgf5ZcnptPl26TQwbhwwe3ZhNdBTTwGjRnWPnKRS2v2nnvK23k8YWFZNREQ2GKTE\njdkIQ3u7FsjMnWtcDbR0qZYk294ODBoEzJ/vfb2fsPgxEy7LmImIEoslyHGjf1FPmaIlkQ4erAUu\nertRRYzeXlzWayRJiaksYyYiSjSWIFcaq3LjQYPMgxhdkkp8WcZMRFSRWIKcVFblxnajJElLTGUZ\nMxFRojFIqTRmibWNjVpehtXzkpaYyjJmIqJEY5BSaYwSa+vrgU8/NR5lyU+6TVKAArCMmYgo4Rik\nxIGbCpX80l1Ay0XZuBHYtKl021QqeaMn+VjGTESUaAxSwlYckEyaZL0KspHm5u5RBKup77ds0bZ1\nGgTlb9evn3bzu7TXS1+stvOjjJmIiOJJKRX72/Dhw1UitLcr1dCglBaOaDeRwvv6LZ223lc6bfy8\n4n0YHbOhQWu365vdc/w4f6d98eP4REQUCwAWKgff/yxBDpNZyawREW10wIzd6skNDdqljylTnJXp\nOumbVWlvR4f53C52x3DaF5YWExElAkuQ48hNaaxdhYrV4/m5GU7LdJ30zWwbfVI1u0tW5faFpcVE\nRFWFQUqYzAKL4oUFnVSo2E2fr49iOC3TdVK2a7bNlCnOpuAvty8sLSYiqioMUsJkFlicd577ChWn\nlS1Oy3SNtrN7js7pyEc5fWFpMRFR1WGQEiazwOKww7zvz66yxWkwU7xdY6N2cxI4OR358NoXlhYT\nEVUlBilhKw4sAPclyE7pZbxnnaXdv/126zLd/L6tXKndnJT2uhn5YMkwERE5xCAlak7zOdxymszq\nB79HPsLsOxERxRZLkKNmVkqslyA7Ke01UsllvJXcdyIissUS5Ephlc9RzoiCWTLrsmXxH5FgCTIR\nEYFBSvSs8jnMLgVNnmw/ZbxVuW65l07crDXkBUuQiYgIDFKiZ5XPYTZysGqV/eiKVUlxOTkvYeSL\nsASZiIjAnJR4czONvp6vkZ/D0revFtCY8fLeh5UvMmmSFqx1dWmrOU+YAEyf7t/+iYgoMsxJSQK7\nCdby6bkm+aMcVgFK8Sy3ToWRL9LRAcye3b3Cc1eXdj/uuTREROQrBilxZnQpqMbkLaupAcaNK81h\nMeN1BC2MfJGgyrKJiKiiMEiJu+LJz8xWRt6ypXvkwQ9mybFh5IuwuoeIiMAgpfKk0/7sp7HR/DGr\n5FivE7e5qQgyG5Xp2zfYqiIiIooVJs5WGj2AcHpZJ5UyHmGZONE8EdXv5FijPjc0mAc3RtvX12sB\n06ZNzvZBRESxxcTZpDIayTAbFUmlgD59jB+bO9f8GH5fbnGbY2J0jttuWxig2O2DiIgqHoOUckya\nBNTWal+ktbXa/TAU56lcf71xnsjs2cDq1cb7MAo49EsyZqNrXpNjvQQ9xefo5jyIiCgRGKR4NWkS\nMGNGYZnsjBnhBSr5rPJEnFbj5OehGCknOdaPiqCwZqENejZdIiJyjEGKFx0dWkBiZObMcPuiKx55\n0PM0nFbjGF2S0ZW7qrEfFUFhVBVx9eXyMMAjIr8ppWJ/Gz58uIqN9nalGhqU0r7GjG9x096uVDqt\nlIj2s729dBsR43MRCa8PYezDSjpt/Bqk0/4eJ4mM/i4aGvx/j4goEQAsVA6+/6u6umfVui9QV1uD\n7XrUOX+S3VT1qRSweXPZfQtc/vT5gwcD69aZz1CbTmsjFkmvoqmpMc7HETGfn4Y0YS2XQESJwOoe\nB4Zf/iD2vXQeMq2deHjxh86eZJeoOWFC+R3zws1Qu9FljU8/1cp8jVTLZQ+uvuwdJ+AjogBUdZAy\nsHePrb+Pv/XfyLR24qd/e8H6SVZfWFZzjwTJbS6FUf7Jpk1ama/ZZHHVUO7L1Ze9Y4BHRAGo6iDl\nyYtHIdvWhGOG7rC17c7n3kOmtROZ1k5s3GwwxG/2RdbeXn6A4nQ0pHi7yZOt5yEpLpU2u1y1erX1\n0LzTFZmtRJVc6eS4XmfTJQZ4RBQMJ4krUd/CSpy99+X3VbplTslt0bLVhRsGkcDpNPHQSeJu/m3i\nROfb6gmiqZTx46lUOOfoNyZ1hiPoxGYiSgwwcda77Mr1GPnbR0raTx0+CFeful8wB3WaeGiXuJsv\nldJ+Oll4MH+KeRHz7cr5vESVXMmkTiKiWHGaOMsgxYJSCkMuNp4+PtvW5O/BnFaWWAUQXqVS2uy0\ngHaJyCwIKvdLParqGVbtEBHFCqt7fCAiyLY1IdvWhK8OLlwDR89bWb1+oz8Hc5p4qI+OOJFKaV/Q\ndvQv6qBmnNVFlVzJpE4ioorEIMWhf046DNm2Jtx05vCC9gN+8wAyrZ246/n3yjuA08RDJ5du8rd1\nMvLSt6/1jLOpFDBuXPkJpFElVzKpk4ioIjFIcenYfXZEtq0JL116TEH7BX99AZnWToy65hFvO3Za\nWWJWItzYaDzK0tUF9OxpPwJjNZ9FV5d2OajcSpyoqmdYtUNEVJGYk+JF0WytmdNvNNzs7SvGoKbG\n5xySjg7gnHOAjXmXmerrgSOOAObPN36OnnthlZsxeLB9Qq7bnJTiWW2LZ621e5yIiBKJOSlBMZg4\nLXvDacgOW4OzD80UbLrLJXORae3EkhVr/e1DcaCxcaN5gAJol3MA69wMo0sixdzMHmo3wRwX8yMi\nIhscSXHLQTnrq//9BE2/f6Jkk4kjd0XLsXsFc3wrjY3AypXdgUF+7kl+6bE+smG2/5oa4LbbnI12\n2L1OLAsmIqpaLEEOioty1s1dW7DblHsNd+O5hNns+Fby++bkEku/fuaLDeYHNV766eTSE8uCiYgS\njZd7guKinLU2VbO1hHlA3jpBQHcJ8/ovbFZMLp7OXb9047XPzc3aSMWWLdpPo2Bj9WrzfW3YoE3D\nbzfFvN3rxLJgIiKywSDFLY/lrE/l1glqO2lYQfvev7rffBVmo7yNtWuBurrC7azmQvFSamsXKKxa\nZZ9LYvc6sSyYiIhsBBakiEgPEXlWRF4UkVdF5LJce18ReUBEluR+bh9UHwJRZjnr6QcNRratCc9O\nGVXQrq/C3HzL092NRnOXbNwIbLdd4fFvu01bgbm4zNhrqa2TJNp8Risk271OLAsmIiIbgeWkiIgA\n6KmUWicidQCeADAZwEkAViul2kSkFcD2SqkWq33FKiclAJnWTsP2d646HuI0b2PSJO1LPn+yt3S6\ne2TCbalvR4d2WccsN8VIOm1/DJYdExFVvVglzopIA7QgZSKA2wCMVEq9LyIDADyilNrT6vlJD1J0\n5//pOcx56f2S9idmjMegTz/qbmhsBHr16v6i32038xLk+nrtssymTd1tTpNfOzqA8eMLn2tGpDAR\n1ugYdtVFdhjgEBElQiyCFBFJAVgEYDcANyqlWkRkjVKqT+5xAfCxft9MtQQpAICODjx96e9w+smX\nlTz08/k34/svztUCgo1lrhnkpNTXrEy4OCApvm92jHLKjssNcIiIKDZiEaTkdaYPgH8C+DGAJ/KD\nEhH5WClVkpciIhMATACAwYMHD1/mdm6QSpX3Rf55qg57XfjPkk3qN2/Em9ecVN5xnJT6WpU751/a\nMXtvio9RTtkx51UhIkqMWJUgK6XWAHgYwLEAVuQu8yD306CsBVBKzVRKjVBKjejfv38Y3YyHvFld\ne3RtQnbaWGSnjS3YZGNtPTItc5BpmYMvUrXejqOUefmwzqzKRw8M9DJms/WEip9fTtmx2Wy3bmbB\nJSKiihJkdU//3AgKRORLAL4B4A0A9wAYl9tsHIC7g+pDRTL5ws7+5UfItjWh5YW7Ctr3vPAuZFrm\nYOFOXzHfp1mJst1U9E7LhP3ezgjnVSEiqjpBjqQMAPCwiLwE4N8AHlBKzQHQBuAbIrIEwOjcfdLZ\nfJFPPGsksjechsduOrdgk1POvBqZljn4adNPS/e5/fbmox1G5cM6p2XCfm9nhPOqEBFVH6VU7G/D\nhw9XoWpvVyqdVkpE+9ne7u7xEI+/JZ1W6ZY5hjelXdTRbqlU4f38m4i7Png5f/05+X1x+9oF/boT\nEVEoACxUDr7/uXZPMSeL8MW0yuS0c67Hs1/eraR90e+/i8bPPjV/YnHyqdU5Au7P32h/Tp9LLL0m\nosSJVXVPuUINUip59d5cMDBv4DBMOPkXJQ9fM+danPzqQ4WNRkGC1TkC7s/fbuXmOLx2cRXjoJiI\nyCsGKV5V+uq9ef/rXrvLHhh2yjUlm6Q//i8enTmhe0ba4i87q3ME3J+/3crNcXnt4ijOQTERkUex\nKkGuKJW+em/eKsfbLn0D2d+eUFLCvGz7gVoJ8+k3ouuM75buw+ocvZy/3Wvj9LUrXhHaqnw6KVh6\nTURVjEFKsaSt3jthAgBsnW/l+88WTg636yVzkWntxBsf5OWsWJ2jl/O3WrDQ6WtntCK0Vfl0UsQ9\nKCYiChCDlGJJW713+vSCFZJ//tgfkf20E/dOPqJgs2OvexyZ1k5cMfd163M0ewwwH+XIfw7QvVqz\nm9fOaEVoq/LppKi0oJiIyEfMSalyXVsUdr1kruFj2d+eoK2qnEppoxbTp5duZLZacn09MGuW94UD\ngcI2u8TbJFe+sLqHiBKGibPk2hFXPYT/rP6spP2V352KXhs/00Zk8gMVq9JiQFuteeVK64Ma7aOu\nzvsiimFXvjCAICJyjUEKefb3/Y7Bz46bXNJ+y52XY/SbT3U32JUWA9ZVPU734VZYlS8sDyYi8oRB\nCnmXKzVe2dAbI35cmph68JC++OsPD7EvLQbsH3eyD7f0kuagRzlYHkxE5AmDFPKutlbLRcmTaZlj\nuOk708ZCzPbj5HJPUCMpU6cGP8oR9zlziIhiivOkkDEnc43kypbzZaeNRfbTTpy4/8CC9iEtc5Bp\nmYP/9N6hdD+nnWZ//DFjSqtX6uq0xNt8Zis5S1GIpFe+mFUDTZ7s31wrLA8mIgoUg5Rq4nSukaKy\nZaRSW5Nmrzv9q8i2NeEf5x1S8JQjzvsDMi1z8H8HndTdOLeoasjo+LNnA+PGFZY033qrVhmU33bb\nbcCoUYX7GzUKuP1241Jps8nOVq3yb64VlgcTEQWKl3viIKwKkQByKL7Y3IU9f35fSbuoLXjn6hO6\nL3t0dGjBSNFlJADaZaFevQrP/4orgNde695m6FDg61/XghC7smjA3WWkcnJIWN1DROQac1IqRZgV\nIkHlUOQCgt0uvAubU7UlD7+x9xr0OM+iVNmoP04/l8Vl0Tq78uji48U5h4SBEBElDIOUShFmhUhQ\nxyoKCG4+8NuYevS5JZv964+TMWzFW96PYySVAjZvNu9X/pf7unWlk84B8a7GYZkzESUQE2crRZgL\nyJmtobNuXXkJpEVT5f/gw+eQHbYGj190VMFmx599PTItc3DFyPHej1XM6PJRfr9yiy0imwWuv97/\nHJKgFz2s1uUAiIjAkZTohT3Xhtk09gH+71xlMhhy+o2GjxWv0Oya1UiKET8vnYQxysEyZyJKIF7u\nqRRRDOdHERjlzvFnx03G3/f9Rskmz19/Brb/fK12x4+clDCE8TpywjgiSiAGKZUk7MTIKP533tEB\nnHnm1rtP77wPTv9uW8lmNz0+E8cuXmCcO9KzJ/D5586qe8IQxuvInBQiSiBfc1JE5CSD2ygR+XL5\nXY3IpEnazKoi2s9Jk6LrS3HuhJcvH6PcCL1NP0cR7X7fvsb7yJ+EzOr1yT9Wv37arfgYxbkZzc3a\n//5zvvafV5CdNhavXXtywWbnHTEBme/Pxpnf+U1p/9avBwYN0o4xaBBw2GHm5+h3boiRICZzK34f\ngYJ8n4K5YIiIkk4pZXsD0AlgNYA7crdVAOYBWALgLCf7KOc2fPhw5auJE5XS/g9ceJs40d/jhKW9\nXamGhsJzqa9Xqq7O+Dzr6rTH89saGrT9KGX9+hgdy+iWvz+rfuZt99Vfz1Ppljklt81So20rUnoe\ntbXOjx/G617Ocf3eHxFRTAFYqBx8/zu63CMi9wP4nlJqRe7+DgBuA3AGgMeUUvsEED9t5fvlHoO1\naQC4T8KMCy/r3xhNoKb/79zq9Rk0qLxJ0hxc2rrlun/g8g++VLK7+Tf/ELuufs/Zsc2O7zc/L9Ux\n/4SIqoSvOSki8ppSamjefQHwqlJqqIg8r5T6anndteZ7kFK83ku+CsjRKeFlJWGrvAmr18dNUms5\nuRkdHcheeR1Gjr205KGWR27FxGfucHf8SpgQjZU8RFQl/J4n5RERmSMi40RkHIB7cm09Aawpp6OR\n0Nekcdoed15yIKyeY/Y6iJgv9Of2GHaam5G5+AJkbzgNb087vuChaSPHI9MyB3tf8Ddnx3e6ZlHU\n7HJcgp6ThYgoZpx+4/wIwK0A9s/dZgP4kVJqvVLqKMtnxpHBKr+W7XFnNElbfb22mrARuwnMzF6H\nVMp68jQ3x3AiN5FZDZS2CvO0sTjqrX9vfXj9Ng3I5FZhXl/Xw/z4lTIhmtWChZUSaBER+clJ4kru\nktCOAE4AcDyAHZ0+z4+b74mzSmlJoKmUlpyYSlVu0qyuvV2pdFpLLk2ntft6m36OQPdjdopfn169\njBNUa2qUamz0dgy787FIzH3gkCbDJNsn0vtp/ck/fnHCrX4TKa+PQTB6H5Xqfh+Lb+l0dH0lIvII\nPifOfh/ALwE8BEAAfB3Ar5VSs4IKnvIlfp6UMDnNzSjezixZVs+XyN9eL3Fevdp9/ofZjLj58hJJ\nP5n0E+y33TdLNjlx/4G47rGbu1dNttlPqNzmxxTNMVOA+SpEVIGc5qQ4HUVZDKAx734jgMVOnuvH\nLZCRlGrktMTVaDuz0Qj9f/tWZclOy2idlDfn76toe6ORlXTLHLUlqhJlp+do1Re714QjKURUgeDz\nSMqTAEYqpTbm7tcDeEQpdajnMMoFjqT4xGmJq9l2xZU9+synU6bYlyU7GbVwUkrd3t496mCy/dVj\nz8eNex9b0v7k9LMxcO1KrS9RVfe4LTO2ek048ywRVSi/S5BvAzAMwN0AFLTclJdyNyilri2rtzYS\nH6SEVR7rtMTVqqQ5nS7tp5MSaCeXJez2U/xFbnU+SuHVLw9B0/gbSh6+8qRhOOOgMiqPyuG2zNjq\nNckP2IiIKojfJchvAbgLwBZoQcrdAN4GsG3uRl6FWbXhdBp3s+30IKF4+n4npcZ+bLPbboUluFbT\n+6dS2PvDd5CdNhZLrj6h4OGL73wZmdZOHH3NI/Z98pvbqfSt3gsGKESUcE6DlLnQSo+/DeBUAKcA\nOEUpdZlS6rKgOlcVwiyPtSpx9bKd1fZOn+tmPw89VBjMrV1bWmatHyuvjLpuS9fWEuZ9uj7Z2v72\nR+uRae1EprUTX2x2WFpdLj9eW7vXk/OpEFFSOElcgZY4ezyAIQDS+s3Jc/24JTpxNuzyWLMS12Kj\nRhX2Z9Qo5/vt2VMrTfZS3m1Telxy69nTvJR86NDCbfPW+blj76MMk2xfWP6x87565bb83el7pm/L\n9X+IKObgc+LsE0qpwwOMlSwlOicljuu1TJoEzJhR2j5xIjB9uvVz9ctX+aNDbhM8vaxFVHysBQuM\nz6HIhz23x0Hn317Sfu7hQ/CLsUMNnlEmP14fK3H8PBERFfE7cXYUtMUE5wP4Qm9XSt1ZTiedSnSQ\nEvSXlhflLMDo5kvSLGHY6DVxI50G3n3X+ey4IlBdXRhy8VzDh7NtTd76YSToIILr/xBRBfA7cXY8\ntJyUY6Fd9jkewFjv3aOtmpu1gCSd1r5I0unoy0rNvtydfOkvX+6s3Sph2Og1cWP5cucBCgAoBRFB\ntq0J2WFrcMbSJwoe1vNWVq/f6K4fZn1z0+6W2wRcIqIYczqSslgptWcI/TGU6JGUOApjJMXtiILV\nysxG+3AzkgJogVLRCM6/dxqKU8+8qmTTG797AJr2HeB83/mCHkmJ48gcEVERv0dSnhSRAC7QUyyV\nswCj02oUtyMKjY32x84/lpvFInv10n4WVVod+N5ryE4bi9f/fkHB5j/603PItHbiO//3lPNj6MaM\ncdfuVhxH5oiIPHIapHwNwAsislhEXhKRl0XkpSA7RhGaPl1Lkk2ltPuplLOkWcD5l6TbyxLXX6+t\n7Jyvvl7rl9Gxis9BBNhmm9L91tYCN92k/W4SIH3pnbe0S0FtTdhhu+59PPPO6q2XgrpStc7Kfeca\n572YtnvR3Gw8nw0RUYVxGqQcC2B3AMegOx/l+KA6VZGSNjfF9OnapR2ltJ/5AYrduRp9SRY/Z8wY\n4zlR1q3Tti3eHgBmzSoMSGbN0vrl5Au5pgY455zS4OsHP7CflC6v/ZlLRiPb1oRfHV84sLjrz+5G\n5vQbsbTlMuv3PuicFK+S9vklomRwUqcc9S3286RU09wUXs7V7DkTJyrV2Fg690ldnVL19eW9nhMn\nGs+ros/f4mDBQtvjptMq22dHw/lWbpj/pulzYrdQYDV9fokoFuDnPClRi33ibDXNTeHlXK2eAzif\nE8XN62mW/Gu3XzfrKOWV+26BYJeWf5Vssk1tDRZfflx3QxwTW6vp80tEseDrPClRi32QUk1zU3g5\nV6vnAPaLEzo5htm+/d5vPpMv9x80T8UDg/YraX/50mOwbY+68BaUdKqaPr9EFAt+V/eQlWqam8LL\nuVo9x81r1Lev87wJPe/ECa/vk0kl083HpZFta8KtZx9Y8NCwS+ch09qJx9fWejteUKrp80tEFYVB\nih+8LAJXqbycq9VzjB6rqysNMmpqtAUFna4WbVaCXFP0kS/nfbKpZDpqry8j29aEF391TMHTzsr2\nQub0G3HlJRBuAAAgAElEQVT+8T8LdtVrp6rp80tElcVJ4krUt9gnzirlbhG4SuflXK2eU/zYxIml\nibNmN6uEU6OF/CJ+n4ySbNMtc9SWKBNnlYr8dSGi6gKHibMcSfFLpc1N4bTk1Gg7L+fq5jl/+xuw\n0eEU9Ealu3qf8xcYHDQIOOywwu3WrQMmTy48N7eluC63z151PLLTxuInC/5c0D7k9BuRae3Eu312\niKYEuNI+v0RUFZg4W42cVpiEUYnix2KC+RUoVvurq9Muy5gFQEaPW52vl9enKNn2jX5pHHvujSWb\nXf7wzTjzJ6cyWCCiRGJ1D5kLan0dP/tiRKSwCsUoIHCzP6fMztfL62MU2Ihgk9Rg95/dXbL54L4N\neOyio7z0mogothikkDmnJadhlKaaHaNYQwMwbpw2fbxV6a7T/blhdr5eX5/iEuSiQOeEs67BiwNL\n1/N845qT0GOnAdGXLBMRlclpkBKzWkgKhcEX49Z2L9sF0ZfGRm3hP7dziZjtr9w+ujmW3euzYIG2\nSrNS2s9evbT8mJy7b/8f7eeRJ2PyIeO3tu/1P3cCAP558UX4KsBAhYgSj4mz1chpyWkYpalmx7j+\nem+JnEb709XVlS5SaPe41fl6eX0mTdISevXZcLu6tACltuj/Cw0NOGHCt5G9ZRz+fcOZBQ99+7tX\nIfNyH1x6z6vmxyEiSgInJUBR3yqiBLnSOC05DaM01e9j6PvLL1U2KkFubNRu+cd12xe32+sl0cU3\n/fnF+8nbZgtgWsJMRFRJwBJk8kWQpal6+e6ZZ3Zf/jAqCy7e3q7ct7m5dJSjqwuYPVv7XT+flSu1\nW/65uT1fJys+5/fTbD0hpWyPKwCy08YiO20smp+fW/BYprUTmdZOrFz3hXV/iYgqCBNnq1EcFrlz\nWnqs9wtw1+eoFs2ze21TKfMkXKMApl8/YNWq0vbGRmDlSixa9jFOnvFkycM3nPFVHL/fwDJOhIgo\nOKzuIXNxWPXWTamw1WrJZn2OatE8u9e2Vy9g/frSx3v2LEie3aqjAzjnnMK5W+rrgVmzCoKzzzd1\nYa9f3Ffy9BHp7fGPiYe6Pw8iogAxSCFzcVj11u0qxYC7PocViNmUE5f008tr73LV5MPaHsJ7az4r\naV869TjUpniFl4iix1WQyVwcVr11u0qx2z5PnVpaqVNfb1+Z5Gaae32UI3/RQzN6P/v2NX7crB1w\nnSezoPVoZNua8OsT9i5o323Kvci0duLNFWstn09EFBcMUqpRHFa9NUsgLWa1WrJdn4tHLOxGDfV8\nErOVlosDmB/+0NkaQxGtKPy9QzLItjXh8aIZa4/53WPItHbi+geXhN4nIiI3AgtSRGRnEXlYRF4T\nkVdFZHKuva+IPCAiS3I/tw+qD2SiuVlL5EyntcsM6XS4SbNAd55JsZoaLSkU0EZbNmzQLnUA2oyz\n+ghMKqXdN+vzlCnApk2FbZs2de/L7DnFibz68Y0CGKPckvzzM3ptV6823t6s3Qc7921Atq0Jb18x\npqD9dw++iUxrJ3a5uDOwY3vidpFHIkqswHJSRGQAgAFKqedEZFsAiwCcCOBsAKuVUm0i0gpge6VU\ni9W+mJOSQPqkZsUmTtRWKy5OFq2t1YKD4hGYxkZt4jen0+Nb5X5YTamfSjkf/QHM9xNl1VFeXsvE\nCb/DvZ+WTmz30qXHYLsedcH1w04cKs+IKHCR56Qopd5XSj2X+30tgNcB7ATgBAC5CSswG1rgQtVm\n7lzz9smTSy+jbN5sHCSsWlV4SUbnJfejZ0/zx9wEKPpIkJExY9y1+8FgFGjG1DORHbYGfxx/YMGm\n+146D5nWTjyy+MPg+mPFajSLiKpOKNU9IpIB8BiAfQAsV0r1ybULgI/1+2Y4kpJAbqp7nCgeibCZ\nX8SQH4sT1tUBt95q/r9+L/0ql4PRm08/34R9L51Xsslx++yIGWcOD6ZfRuJQeUZEgYtNCbKI9ALw\nKICpSqk7RWRNflAiIh8rpUryUkRkAoAJADB48ODhy/xeNI6iVVtrPDrh9rKKzo8VnN0GTk5XZnZ6\njKD+Fl2+FrteMhddW0q3f/uKMaip8Tm4LBaHOXyIKHCRX+7JdaIOwB0AOpRSd+aaV+TyVfS8FcNx\nZaXUTKXUCKXUiP79+wfZTYqCWSDS1WV9ucSM0QrOTrbL56YsOpXS8iSmTw9u2QC/uHwt3rpiDLJt\nTfh/o/coaN/lkrnItHbiP6ttZgkuh9kCkevWMYGWqAoFWd0jAP4A4HWl1LV5D90DYFzu93EA7g6q\nDxRjZtU96bSWCFtXlLxpNQLh1wrOEyYYtxusUIzZs70FJGYBmJfAzCmPJeeTR++ObFsT5v2/Iwva\nj7jqYWRaO3HbU1l/+wl0V54Vvx5muUdElGxOViH0cgNwOAAF4CUAL+RuYwA0ApgPYAmABwH0tdsX\nV0FOoPZ2pRoaClcCbmjoXv23eHXhxkbj1YNTKX9XcJ44sXulYqOVk4v342XV5Pr6wnOorw9mdeni\n41r108F5bNrcFd4qzMWrWOu3dNr8OWGs2E1EvoDDVZADC1L8vDFISSijgMDsMaMvLP1mxuhLy+6L\nbOjQwn0PHWq9f6tAy00fovyCdXIeRU6ZscAwWFn/xSZ/+iRi/F6L+HYORBQdp0EK1+6haFjNh7Fg\ngfEcKkbMVg822n99vfb1lT/JW/4cHHvvDbz2Wum+hg4FXn21tN0uydPpnB9Rzw1SRrJq+9PL8PO7\nXilpv/3cg3DE7mXkkrntExNuiSpKbKp7/MAgJYGsvlTefbf8idPcrrKczbqvvLGrmnH6xRn1F6wP\nZb8ffPI5vnbl/JL2o/f6MmadfaDBM2y4DdxYukxUUWJR3UNkavly83YvJchO929k2TItUHDLrmrG\nrA/LlmnJuCLaT7Ngys05lMOHBSd37N0D2bYmZNuaCtofeuNDZFo7kWnthKv/ELldusHrOXAKfqJY\nY5BC0bD6UnFTCmxWFeN2RWcv8/DYVc1Y9UEPxKwCsrBWpfZ5wUk9WDlm6A4F7UMu1kqY/7vmM2c7\ncrP6s5dzsFtQkogixyCFomH1pWJWClx8OaauTitXdrr/+vrS0mYnhg41brf7377ZnB9OhLlyckAL\nTs783ghk25rQ8f2DC9oPbXsImdZO/PH6v5e1/wJezoFT8BPFn5Ps2qhvrO5JKKuKFrelwE73r7dZ\nVQs5re7xco5Wx/KruieOpbjt7eqzbXsbVgTt86v7SrYNpf9uK4iIyDdwWN3DkRRyz6/r+G6G8/2y\nYIGWmGumR4/C+wMGWO9v9Gjtf+76bfTowseLz9HsUlaNT3+KYV3CcPsZmDIFPdZ+guy0schOG1vw\n0NrPN2/NW9l4e4iXYHzIxSGigDmJZKK+cSQlRsKYj2LiROP/4dbWOj+uUT/t5lupqTFuHzXK+Bij\nRrnb3uo5xf+r9/qaepkEzS0vnwGTUYvrDz3dcHTl+QF7BNd/ndnnLH++HiIKBDhPCgUijHJZs8UH\njbidN8NqP1bbG/2deFks0EtptBthlOJ6+QzYPOedletx1G8fKXn4Oy/ej2n33RBMKXHUpd9mOjq0\nvBini1YSVSCWIFMwrEqH3TK7ZOCmBNmv/oT1peSmX15e06AvYXR0eCuZtqm+GdKvJ7JtTXjnLz8q\n2OSv+30TmZY5yFz0L2DSJPu+ubkE5edn2S+sOCIqwCCF3PHrS9DqH2M3Jchmx+3b1/k+3ByvXG76\n5SWw8LmcuID+npmx6q/D6huZOhXZG05DdtpY7PffxQWPZbZrQqa1Ex9++rl539x8uccxJ4UVR0SF\nnFwTivrGnJQY8SsnxSp3wo+cFLMFCa1yENzmmHjJSXHar3LyfIKqjrGqTvIzL6m9vWDffxj+LcO8\nlekPL7Xvm92ChHFb74cVR1QlwAUGKTB+fAna/WNcbgmy2f71/eXvN19x4GEVcHjZ3qpfcSsbLmbV\nd7/7a3CMlV/aznwVZq9f7nEr1w4j8ZkoBpwGKbzcQ+75UTpsN9Q+fTqwebP2T/Tmzdp9N8c12386\nXbrffOPHF16SGD/e+jwefLDw6+TBB623t+pX2OXY+Zzkc5j1PYjLZQb7bPzsU2R/e0LJ1PsAkLno\nX8i0zEGXFP2T5velm6Cn0Q/ych1RJXISyUR940hKAgU91O5l/2EM/3u5RBQ0p+dttF3Er9PJ0xcY\njqwsGrin/+93WJeH4ja6QxQAsASZYi/oUku3+4+yvDqV0kZ2ouDmvDs6gHHjjM/Bz9fJ5Xvx5NKV\n+O4tz5S0j/7KDrhlnEmVo9v3O64ly0QVyGkJMoMUIl0Yc4x4mVvFLbfBmdvzDuJ1Ku6zWZmzzTE2\ndW3B7lPuNXys5DJRHM67UnDuFvIZ50khciuMklSz/A2/8jrCKMX1+3Uy6rNZMGdzjLpUzdZVmIvp\nU+9/8tkm632Fdd6VgnO3UIQYpBDpwkhaNJtnxGr+ETe8zLPh9rydbO8mwdSoz0qVBiou3ws9WJk0\ncteC9v0um4dMayf+fuHV/p93EnHuFoqSk8SVqG9MnKXQhJG0aFRe7ZewSnGttnebYBpSWfbyVesN\nk2z3+p87/DnvpOLcLRQAsASZyAN9lWSltJ8LFlhvP2mSlgwrov20m7odAA47DBg0SHvOoEHafb94\nvSThtqzcanu3//M261tjo3UfXNq5b4M2ujJsTUH7Z7XbaCXMp98Ideut1juJYuVuO0GXRVfrZa5q\nFvRnyg0nkUzUN46kUCjcrorrZRXdOJZe+83t/7yN+lxfr1Rdnf/nUXSsgyf90XB0ZcmKteUdJyxR\nrkrO1aKTKaR/Q8ASZCKX3JYHeyknDqOMNepKDC/nWNzndeuAVavc7aOMvnXueRh+dOLFJe1nfm0w\nLj9xmPfjBS2MzxNLr6tLSO83q3uI3DJbfdmvdiCeK+/6zUuCafFllNWrjbdbtsx8GNrJELXJ69y0\neAGy08bitWtPLmhvf3r51qqgWArj8xTWZzZOlxiqWcz+jaqN5KhEcZRKmY+M+LE9YD4HiF/X9/Vy\nUT0nRC8XBcIbTdGPU85ojtnrJNLdnn9ugLPztpqDBUDDpi+QnTYWUKokMNHvv/GbY9GjLsSVs630\n7Ws84uRmtW07QX9mgXh8bkkTxvvtAkdSiHRuy4O9lBMHXcYal3LRchNMjV4nkdLJ1PRzc3reRvst\nNmoUgO4S5pMO2Kng4b1+cR8yrZ2Y//oKp2dT2cIovY7L55biV2rvJHEl6hsTZyk0bsuDvZQTB1nG\nmqRy0eLXyaxMWcTdeefvt0ePwu0t1lB69b1PDJNsR13zSFCvgL2w3u+gS6+T9LlNghBK7cESZLLk\npXS2GrgtDzZarTlKYZWLhpE/UDwak04bbzd4sLvzzt9v8WWR11837c7Qgdsh29aEd64cU9C+9MN1\n/uatuHltw3q/gy69DuI8mOPiXZxK7Z1EMlHfOJLiM5YUGguj9C4JJchRlTlbHdfL6tIDBxo/Z+BA\nx10yGllJt8xR/12zwf9zNJKUv2W/zyMOpfhkCSxBJlNxXIk3DpJSzhl0CXKUJalm5+blM+3jYo+3\nP5XFL+5+taT9wmP2wPlH7+58R9W6MrPf55GU1yXBuAoymQtjJd5KFMYqt0lYSTeO5+DlM231nHTa\nU5C3ev1GHPCbBwwfM1r0sES1rszs93kk5XVJMM6TQuaCXom3UoVxfT8JU4zH8Rz8/kx7XPG3b896\n21WYu7ZY/EegWldm9vs8kvK6EIOUqhT0SryVKozSu7iV93kRx3Pw8pkeONDZvj2WwurByqG7Fq5B\ntOslc5Fp7cSiZR+XPimIFakrgd/nkZTXhZg4W7WCXIm3koWxym0SVtL16xzsVlN2cwy7z7TR/szK\nmv1O4kyn1VODhxkm2X7vD8+4Ow+T/bt6L+L49+/330Wl/J35uQK5mRi+33CYOBt5AOLkxiCFKIGs\nKjD8rs4w2p+Xm9t/3A2Ou6lnL9OqIC4YWGXcvt9ePh8xfb+dBilMnCWiaFhVYADhVHu45bYCzqbK\nxGxulRev+w56f7G+ZHtfsLovPsKo5orp+83qHiKKN6sKDCCcag8v3OzHYZXJtfMW4/cPLS3ZbNbf\nL8XRby/0tyqF1X3xEUY1V0zfb1b3EFG8WVVghFXtYcavaiGH5/HTY/ZEtq0JC+75eUH7OadeikzL\nHJwx7rfujmuF1X3xEUY1V4W/3wxSiCgaVhUYYVR7mOnTx78KOJfnsdOUC5G94TRtJeY8T+2w59YS\n5rJHv1ndFx9hVHNV+vvtJHEl6hsTZ4kSys/qHi/H6tOnMJmwT5/u7f2qiCizeuP8K+40TLJdvmq9\nt/4oFctqj6rF6h7LG0dSiMg5fdE2ke5bOQtU+rmQmd2imUbHKp4rJf++0eKRdovWGT1e5jneMHgD\nsm1N+NsPDyloP+Kqh5Fp7cTMC37rfrHQN9/sTqbs6tLuR40LAjrj5fMUt4VQ3XASyUR940gKUQzY\nlfH6+b+zMBbaGzrU+DlDh3rrkx/lww728dnGzYYjK7teeJez8/ayGGPQqnVBwGo9b+V8JIXVPUTk\njF0Zr58ljW5LLcNYYNCuT34saudyH0Na/gUlpQPii393Mrb54nPjY8Sx2qNaFwSs1vMGS5CJyG9O\nynj9+vfEbaml3wsMGj3Hrk9+LGrn8bz/76CTcOVR55Q8/I/zDsGITF/D5xiK6vugWhcErNbzBkuQ\nichvdmW8fpY0ui21DKPM0q5PfpRNezzvHz57J7LTxuLxm84tePiUm55CprUTk//yvPM+RKFaFwSs\n1vN2gUEKETljV8brZ0mj21JLL2WWQ4e6a7frkx9l02We986frEB22li882nhTLZ3v/DfrSXMGDXK\neF9m7WGo1gUBq/W83XCSuBL1jYmzRDGhlz/mJ/oFVdLo9wKDRoqTZ82SZp32yY+yaZ/P+7SbnjRM\ntP2ooXc8kmZ1lbIgoN+q9LzBxFkiohjo6ACmTAGWL9eG8adOdV9q7cM+HnxtBb5/W+m/o1efsi9O\nHbGz9ZMnTQJmztSSk1MpbQQn7mWsfrzuSRHD14KJs0REUevo0L7QN2zobmto0L7wnX5J+LGPPOu+\n2Ix9fnV/Sfvgvg147KKjSp8waRIwY0Zp+8SJ8Q1UfH7NKlpMXwsGKUREUYugLNkNs1WY37piDFI1\nuSqgmK6ia6mKS3tLxPS1YJBCRBS1KMqSPbhi7uuY+djbJe33Tj4CXxnY2/yJcf3+qOLS3hIxfS1Y\ngkxEFLUoypI9uGTMV5Bta8L9FxxZ0H7c9Y8j0zIHlx91bumT4ryKLkt7u1X4a8EghYgoKFGUJZdh\nzx23RbatCW9dMaag/ZaDvo1MyxxkWuZ0N8Z5FV2W9nar8NeCQQoRUVCam7UExXRaG15Pp90nLPqx\nD5dSNYJsWxOybU3INBZ+wenBytprrw/s+GWL4DWLrQp/LZiTQkTVK4almbHU0YF/3HwPLvza90oe\nuvl7I/CNoTtE0KkYMPr8APxMOcDEWSIiKzEtzYydotdpZUNvjPhxR8lmB2a2x9/POzTs3kXH6PNT\nX68lqW7a1N3Gz5QhBilERFZiWpoZOxavU+b0Gw2f8s6VYyBWCxkmgd2q4Pn4mSrBIIWIyEpMSzNj\nx8Hr9NO/voA7n3+vZJPHfnYUBhfltCSGk1XBdfxMlWAJMhGRlQovzQyNg9fp2u/sj2xbE+6YeEjB\nJkde/TAyrZ2Y8chbQfYwGn6UkZMtBilEVJ0qvDQzNC5ep+Hpvsi2NWHx5ccWtE+7743uVZiTwuh1\nqa8H6uoK2/iZKktgQYqIzBKRD0Xklby2viLygIgsyf3cPqjjExFZqvDSzNB4eJ22qU1tLWGury38\nmtGDlc83GUy1X0mMXpdZs4Bbb+VnykeB5aSIyJEA1gG4TSm1T67tKgCrlVJtItIKYHulVIvdvpiT\nQkRU2W55/G1c3vl6SftfJ3wNB+/SGE4n4lhyHsc+hSAWibMikgEwJy9IWQxgpFLqfREZAOARpdSe\ndvthkEJElAzvfrwBh097uKS9adgA3Nh8QHAHjmPJeRz7FJK4BilrlFJ9cr8LgI/1+1YYpBARJYtS\nCkMunmv4WLatyf8DxrHkPI59Cknsg5Tc/Y+VUoZ5KSIyAcAEABg8ePDwZU7r0YmIqKKc9Ydn8PiS\nlSXtz04ZhS9v28Ofg8Sx5DyOfQpJXEuQV+Qu8yD380OzDZVSM5VSI5RSI/r37x9aB4mIKFy3n3sw\nsm1NuPXsAwvaD5o6H5nWTvz138vLP0gcS87j2KeYCTtIuQfAuNzv4wDcHfLxiYgopo7a68vItjXh\n1cu+WdDecsfLyLR24tAr53vfeRxLzuPYp5gJsrrnzwBGAugHYAWAXwG4C8DfAAwGsAzAaUqp1Xb7\nYk4KEVF1MptbZenU41Cbcvn/7DhW0sSxTyGIRU6KXxikEBFVt2nX3IEZH5Xmp8z58eHYZ6feEfSI\nysEghYiIkiGvVHdJ4874xvdnlGxy9qEZXPqtvSPoHHnBIIWIiJLBoFR3CwS7tPzLcPNASpjJVwxS\niIgoGWxKdUdf+yiWfriu5OEXf3UMen+prvR5FLm4liATERG5Y1Oq++BPv45sWxOu+87+BQ/vd9k8\nZFo7cd8rHwTdQwoIR1KIiCjeXE4fv3r9RhzwmwdK2vffuQ/u+tFhQfaUHOLlHiIiSg6PpbpmJczv\nXDkG2uosFAUGKURERDkX/eNF/G3huyXtj1w4Epl+PSPoUXVjkEJERFTkueUf46TpT5a0X3jMHjj/\n6N0j6FF1YpBCRERkYuPmLdjj5/caPsYS5uAxSCEiInJg71/eh/Ubu0ra3/jNsehRl4qgR8nHIIWI\niMiFWxe8g8v+9VpJ+59+cDAO3bVfBD1KLgYpREREHvx3zWc4tO2hkvZj994RN501PIIeJQ+DFCIi\nojIopTDk4rmGjzFvpTyccZaIiKgMIoJsWxOybU0YuWf/gscyrZ3ItHZixaefR9S76sCRFCIiIoce\nffMjjJv1bEn7Fd8ehu8ebDJ9P5Xg5R4iIqKAbNi4GUN/eX9Je79e22Dhz0dH0KPKwiCFiIgoBGZT\n7y+dehxqU8yqMMIghYiIKETXzluM3z+0tKT9X+cfjmGDekfQo/hikEJERBSBpR+uw+hrHy1pP+tr\nafzmxH0i6FH8MEghIiKK0JYtCrtcwhJmIyxBJiIiilBNTXcJ8147blvwmF7C/MmGTRH1rjJwJIWI\niCgk97z4X/zkz8+XtE9vPgBjhg2IoEfR4OUeIiKimFqzYSP2//UDJe37DuqNe84/PIIehYtBChER\nUQUwK2F++4oxqKmRkHsTDuakEBERVQA9b6V4xtpdLpmLTGsn3v5oXUQ9ix5HUoiIiGLkxf+swQk3\nLihpv2D07rhg9B4R9Mh/vNxDRERUwTZ1bcHuU+41fKzSS5gZpBARESXEfpfNwyeflZYrv/7rY/Gl\n+lQEPSoPgxQiIqKEuf3pZfjFXa+UtLefezAO371fBD3yhkEKERFRQn3wyef42pXzS9pHf+XLuGXc\ngRH0yB0GKURERFXArIQ5znkrLEEmIiKqAnoJ8+iv7FDQrk+9/8Enn0fUs/JxJIWIiChBFixdieZb\nnilpn3nWcByz944R9KgUL/cQERFVsc82duErv7yvpD0OwQqDFCIiIgIANN/yNBYsXVXQNuHIXdB6\n7F6RTL3PIIWIiIgKLFmxFk2/fwIbu7Zsbdtv5z6YPf5A9GmoD60fDFKIiIjI0LovNmNi+yI8vmRl\nQfs95x+GfQf1Cfz4DFKIiIjIklIK//vQUlzzwJsF7Vd8e1jJgod+YpBCREREjj2xZCXO/ENhVdDb\nV4wJJGfFaZBS6/uRiYiIqOIcvns/ZNua8P4nn+GUGU/hwMz2kPBzagswSCEiIqKtBvT+Eha0Hh11\nNwBwxlkiIiKKKQYpREREFEsMUoiIiCiWGKQQERFRLDFIISIiolhikEJERESxxCCFiIiIYolBChER\nEcUSgxQiIiKKJQYpREREFEsMUoiIiCiWGKQQERFRLDFIISIiolgSpVTUfbAlIh8BWBZxN/oBWBlx\nH6LA864uPO/qwvOuLnE677RSqr/dRhURpMSBiCxUSo2Iuh9h43lXF553deF5V5dKPG9e7iEiIqJY\nYpBCREREscQgxbmZUXcgIjzv6sLzri487+pScefNnBQiIiKKJY6kEBERUSwxSDEgIlkReVlEXhCR\nhbm2viLygIgsyf3cPup++klE9sydr377VEQuEJFLReS9vPYxUfe1XCIyS0Q+FJFX8tpM318RuVhE\nlorIYhH5ZjS9Lp/JeV8tIm+IyEsi8k8R6ZNrz4jIZ3nv+03R9bw8Judt+rlO+Pv917xzzorIC7n2\nJL3fO4vIwyLymoi8KiKTc+2J/hu3OO/K/htXSvFWdAOQBdCvqO0qAK2531sBTIu6nwGefwrABwDS\nAC4FcGHUffL5/I4EcACAV+zeXwBDAbwIYBsAQwC8BSAV9Tn4eN7HAKjN/T4t77wz+dtV8s3kvA0/\n10l/v4sevwbALxP4fg8AcEDu920BvJl7XxP9N25x3hX9N86RFOdOADA79/tsACdG2JegjQLwllIq\n6gn0AqGUegzA6qJms/f3BAB/UUp9oZR6B8BSAAeF0lGfGZ23UmqeUmpz7u7TAAaF3rGAmbzfZhL9\nfutERACcBuDPoXYqBEqp95VSz+V+XwvgdQA7IeF/42bnXel/4wxSjCkAD4rIIhGZkGvbQSn1fu73\nDwDsEE3XQnE6Cv/x+nFuqHBW0i5z5TF7f3cC8J+87d7NtSXROQDuzbs/JDcM/KiIHBFVpwJk9Lmu\nlvf7CAArlFJL8toS936LSAbAVwE8gyr6Gy8673wV9zfOIMXY4Uqp/QEcB+BHInJk/oNKGytLZFmU\niNQD+BaAv+eaZgDYBcD+AN6HNkScaEl+f82IyBQAmwF05JreBzA493fwUwB/EpHtoupfAKruc13k\nDI6j1ogAAATvSURBVBT+RyRx77eI9AJwB4ALlFKf5j+W5L9xs/Ou1L9xBikGlFLv5X5+COCf0Ib+\nVojIAADI/fwwuh4G6jgAzymlVgCAUmqFUqpLKbUFwM2owGFQh8ze3/cA7Jy33aBcW2KIyNkAxgJo\nzv3jjdzQ96rc74ugXaffI7JO+szic10N73ctgJMA/FVvS9r7LSJ10L6oO5RSd+aaE/83bnLeFf03\nziCliIj0FJFt9d+hJR29AuAeAONym40DcHc0PQxcwf+w9D/qnG9Dey2SyOz9vQfA6SKyjYgMAbA7\ngGcj6F8gRORYABcB+JZSakNee38RSeV+3wXaeb8dTS/9Z/G5TvT7nTMawBtKqXf1hiS937l8mz8A\neF0pdW3eQ4n+Gzc774r/G486czduN2hDwC/mbq8CmJJrbwQwH8ASAA8C6Bt1XwM4954AVgHondd2\nO4CXAbwE7Y95QNT99OE8/wxtqHMTtOvP51q9vwCmQPtfxmIAx0Xdf5/Peym06/Ev5G435bY9Off5\nfwHAcwCOj7r/Pp+36ec6ye93rv2PAM4r2jZJ7/fh0C7lvJT3uR6T9L9xi/Ou6L9xzjhLREREscTL\nPURERBRLDFKIiIgolhikEBERUSwxSCEiIqJYYpBCREREscQghYhM5VZKTercOEQUcwxSiCgQuZlN\nY69S+klUjRikEJGdlIjcLCKvisg8EfmSiOwvIk/nFuj7p75An4g8IiLXichCAJNF5FQReUVEXhSR\nx3LbpETkahH5d+75P8y1jxSRx0SkU0QWi8hNIlKTe+wMEXk5t69pubZTReTa3O+TReTt3O+7iMiC\n3O/Dc4unLRKR+/OmRS/oZ7gvJxE5xf9BEJGd3QGcoZT6gYj8DdpMlRcB+LFS6lER+TWAXwG4ILd9\nvVJqBACIyMsAvqmUek9E+uQePxfAJ0qpA0VkGwALRGRe7rGDAAwFsAzAfQBOEpEnAUwDMBzAxwDm\niciJAB7P9QPQVvVdJSI75X5/LLeOyQ0ATlBKfSQi3wEwFdpKsAX9JKJ4YpBCRHbeUUq9kPt9EYBd\nAfRRSj2aa5uN7lWzgbyF6wAsAPDHXHCjL3h2DIB9ReSU3P3e0AKhjQCeVUrpIyJ/hjbV9yYAjyil\nPsq1dwA4Uil1l4j0yq21tTOAPwE4ElqQcieAPQHsA+ABbVkTpKBNE2/UTyKKIQYpRGTni7zfuwD0\nMdswZ73+i1LqPBE5GEATgEUiMhyAQBuFuT//SSIyEtraI/ns1u14EsB4aGuuPA5tlOQQAP8DYDCA\nV5VSh9j1k4jiiTkpROTWJwA+FpEjcvfPAvCo0YYisqtS6hml1C8BfARtxON+ABNzl2MgInvkVhwH\ngINEZEguF+U7AJ6AtiLt10WkX27V1jPyjvc4gAsBPAbgeQBHAfhCKfUJtMClv4gckjtOnYjs7d/L\nQERB40gKEXkxDsBNItIAbXn38SbbXS0iu0MbPZkPbXXxlwBkADyXW17+IwAn5rb/N4D/BbAbgIcB\n/FMptUVEWnP3BUCnUuru3PaPQwt8HlNKdYnIfwC8AQBKqY25S0q/F5He0P69uw7ayq9EVAG4CjIR\nxULucs+FSqmxUfeFiOKBl3uIiIgoljiSQkRERLHEkRQiIiKKJQYpREREFEsMUoiIiCiWGKQQERFR\nLDFIISIiolhikEJERESx9P8BTFYnHD7Vjv0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1224374e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### your solution\n",
    "auto_file = pd.read_csv(\"Auto.csv\")\n",
    "x = auto_file[[\"horsepower\"]]\n",
    "y = auto_file[[\"mpg\"]]\n",
    "\n",
    "ols.fit(x,y)\n",
    "\n",
    "def print_hash(letter):\n",
    "    print(\"%s: ####################################\"% letter)\n",
    "\n",
    "print_hash(\"A\")\n",
    "print(\"Coef of mpg is: %.10f -- it is derived from the formula \\\n",
    "of coefficients. It is an optimization formula to minimize \\\n",
    "RSS. \" % ols.coef_)\n",
    "\n",
    "print_hash(\"B\")\n",
    "print(\"Coef determination of our model is: %.10f -- it is based \\\n",
    "the calculation of R squared https://en.wikipedia.org/wiki/Coefficient_of_determination\" % ols.score(x,y))\n",
    "\n",
    "print_hash(\"C\")\n",
    "print(\"The fomular is Y = %.5fX + %.5f\" % (ols.coef_, ols.intercept_) )\n",
    "\n",
    "print_hash(\"D\")\n",
    "print(\"The horsepower of 98 based on our model predicts mpg at: %.2f \" % ols.predict(98))\n",
    "\n",
    "print_hash(\"E\")\n",
    "# plt.figure(figsize=(x.max().values[0],y.max().values[0]))\n",
    "plt.figure(figsize=(9, 6))\n",
    "plt.scatter(x,y, color=\"red\")\n",
    "plt.xlabel(\"horsepower\")\n",
    "plt.ylabel(\"mpg\")\n",
    "plt.plot(x, ols.coef_*x + ols.intercept_)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 3\n",
    "Without using the sklearn library, calculate the coefficients manually by employing the numpy library, as well as the coefficient of determination and prediction. Note that the results should be exactly the same with the results in question (2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coef: -0.1578447334, Intercept 39.9358610212\n"
     ]
    }
   ],
   "source": [
    "### your solution\n",
    "Xnp = np.array(x.horsepower.values)  \n",
    "Ynp = np.array(y.mpg.values)\n",
    "\n",
    "A = np.vstack([Xnp, np.ones(len(Xnp))]).T\n",
    "m, c = np.linalg.lstsq(A, y)[0]\n",
    "\n",
    "print(\"Coef: %.10f, Intercept %.10f\" % (m,c))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}