unit03_mult_lin_reg/demo1_glucose.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demo:  Predicting Glucose Levels using Mulitple Linear Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this demo, you will learn how to:\n",
    "* Fit multiple linear regression models using python's `sklearn` pachage.  \n",
    "* Split data into training and test.\n",
    "* Manipulate and visualize multivariable arrays.\n",
    "\n",
    "We first load the packages as usual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diabetes Data Example\n",
    "To illustrate the concepts, we load the well-known diabetes data set.  This dataset is included in the `sklearn.datasets` module and can be loaded as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets, linear_model\n",
    "\n",
    "# Load the diabetes dataset\n",
    "diabetes = datasets.load_diabetes()\n",
    "X = diabetes.data\n",
    "y = diabetes.target"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can print a description of the data as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ".. _diabetes_dataset:\n",
      "\n",
      "Diabetes dataset\n",
      "----------------\n",
      "\n",
      "Ten baseline variables, age, sex, body mass index, average blood\n",
      "pressure, and six blood serum measurements were obtained for each of n =\n",
      "442 diabetes patients, as well as the response of interest, a\n",
      "quantitative measure of disease progression one year after baseline.\n",
      "\n",
      "**Data Set Characteristics:**\n",
      "\n",
      "  :Number of Instances: 442\n",
      "\n",
      "  :Number of Attributes: First 10 columns are numeric predictive values\n",
      "\n",
      "  :Target: Column 11 is a quantitative measure of disease progression one year after baseline\n",
      "\n",
      "  :Attribute Information:\n",
      "      - Age\n",
      "      - Sex\n",
      "      - Body mass index\n",
      "      - Average blood pressure\n",
      "      - S1\n",
      "      - S2\n",
      "      - S3\n",
      "      - S4\n",
      "      - S5\n",
      "      - S6\n",
      "\n",
      "Note: Each of these 10 feature variables have been mean centered and scaled by the standard deviation times `n_samples` (i.e. the sum of squares of each column totals 1).\n",
      "\n",
      "Source URL:\n",
      "https://www4.stat.ncsu.edu/~boos/var.select/diabetes.html\n",
      "\n",
      "For more information see:\n",
      "Bradley Efron, Trevor Hastie, Iain Johnstone and Robert Tibshirani (2004) \"Least Angle Regression,\" Annals of Statistics (with discussion), 407-499.\n",
      "(https://web.stanford.edu/~hastie/Papers/LARS/LeastAngle_2002.pdf)\n"
     ]
    }
   ],
   "source": [
    "print(diabetes.DESCR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The target values are stored in the vector `y`.  The attributes for the diabetes data are stored in a data matrix, `X`.  The size is is number of samples (442) x number of attributes (10)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "num samples=442  num attributes=10\n"
     ]
    }
   ],
   "source": [
    "nsamp, natt = X.shape\n",
    "print(\"num samples={0:d}  num attributes={1:d}\".format(nsamp,natt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In class exercise: Print the (normalized) ages of the first 5 subjects?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In-class exercise:   Print the attributes S1-S3 for subjects 10-15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In class exercise:  Create a scatter plot of the target variable, `y` vs. the BMI.  Does there seem to be a relation?  What about `y` vs. the age?  Which is a better predictor?  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using Simple Linear Regression for Each Feature Individually\n",
    "\n",
    "As a first attempt to predict the glucouse level, we could try *one attribute at a time*.  That is, for each attribute $x_k$, we could attempt to fit a [simple linear regression](../simp_lin_reg/readme.md) model:\n",
    "$$ \\hat{y} = \\beta_{0,k} + \\beta_{1,k}x_k$$\n",
    "where $\\beta_{0,k}$ and $\\beta_{1,k}$ are the coefficients in the simple linear regression model using only the attribute $x_k$.\n",
    "\n",
    "Now, we saw in class, that the goodness of fit in a simple linear regression model is given by the coefficient of determination.  Let $R_k^2$ be the coefficient for predicting $y$ from the $k$-th predictor $x_k$:\n",
    "$$R^2_k = \\frac{|s_{x_k,y}|^2}{s_{x_k}^2 s_y^2}.$$\n",
    "The following code computes $R^2_k$ for each variable $k$ as well as the coefficients in the linear model, $\\beta_{0,k}$ and $\\beta_{1,k}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 0 Rsq=0.035302 beta0=152.133484 beta1=304.183075\n",
      " 1 Rsq=0.001854 beta0=152.133484 beta1=69.715356\n",
      " 2 Rsq=0.343924 beta0=152.133484 beta1=949.435260\n",
      " 3 Rsq=0.194908 beta0=152.133484 beta1=714.741644\n",
      " 4 Rsq=0.044954 beta0=152.133484 beta1=343.254452\n",
      " 5 Rsq=0.030295 beta0=152.133484 beta1=281.784593\n",
      " 6 Rsq=0.155859 beta0=152.133484 beta1=-639.145279\n",
      " 7 Rsq=0.185290 beta0=152.133484 beta1=696.883030\n",
      " 8 Rsq=0.320224 beta0=152.133484 beta1=916.138723\n",
      " 9 Rsq=0.146294 beta0=152.133484 beta1=619.222821\n"
     ]
    }
   ],
   "source": [
    "ym = np.mean(y)\n",
    "syy = np.mean((y-ym)**2)\n",
    "Rsq = np.zeros(natt)\n",
    "beta0 = np.zeros(natt)\n",
    "beta1 = np.zeros(natt)\n",
    "for k in range(natt):\n",
    "    xm = np.mean(X[:,k])\n",
    "    sxy = np.mean((X[:,k]-xm)*(y-ym))\n",
    "    sxx = np.mean((X[:,k]-xm)**2)\n",
    "    beta1[k] = sxy/sxx\n",
    "    beta0[k] = ym - beta1[k]*xm\n",
    "    Rsq[k] = (sxy)**2/sxx/syy\n",
    "    \n",
    "    print(\"{0:2d} Rsq={1:f} beta0={2:f} beta1={3:f}\".format(k,Rsq[k],beta0[k],beta1[k]))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the best $R^2_k = 0.34$.  That is, only about 34% of the variance is explained by a linear model with any one variable.  We can see this somewhat poor fit in the scatter plot as well where there is a significat variation from the regression line.\n",
    "Also notice that beta0 are all the same no matter which attribute is used to predict the target!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2de3hU1bn/PyvDBCYEmYAYIaJoa6FSjqCpWu2vB+iF3sSIVrS2Xmq1PV4qqGi8VMFL4Ui91Nra2tbWVitYLxGlFq2S0yMttWBAxMoRRYVR8QJBQoZkMlm/P2Z2svfMXvsys+eSZH2eJ0+SfV1r75nvfve73vW+QkqJRqPRaPoXFaVugEaj0WiCR4u7RqPR9EO0uGs0Gk0/RIu7RqPR9EO0uGs0Gk0/ZFCpGwCw7777ynHjxtmu27NnD0OHDi1ug4rMQOgjDIx+6j72D/pKH9euXfuBlHKU3bqyEPdx48axZs0a23XNzc1MnTq1uA0qMgOhjzAw+qn72D/oK30UQrypWqfdMhqNRtMP0eKu0Wg0/RAt7hqNRtMP0eKu0Wg0/RAt7hqNRtMPcRV3IcQQIcTzQoj1QoiNQogF6eW/E0JsEUKsS/9MTi8XQog7hBCbhRAvCiGOKHQnNBqNRmPFSyhkBzBdStkmhAgDzwkhnkyvmyelfChj+68Ah6Z/jgbuSv/WaDQaTZFwtdxlirb0v+H0j1Oe4BOA36f3Ww1EhRCj82+qRqPRaLwivORzF0KEgLXAx4GfSSmvEEL8DvgMKcv+GaBRStkhhHgCWCSlfC697zPAFVLKNRnHPA84D6C2tvbIJUuW2J67ra2N6urqHLvXNxgIfYSB0U/dx/6Bnz62xhNs37WXzmQ3laEKaocPIRoJF7iFKaZNm7ZWSllvt87TDFUpZRKYLISIAo8KIT4FXAm8C1QCdwNXANcDwu4QNse8O70f9fX1UjUbrK/MFMuHgdBHGBj91H3sH3jtY1NLjCuf2UA8UYHhCImEkyycdRgNU+oK20gXfEXLSClbgWbgy1LKd9Kulw7gt8BR6c22AWNNux0AvB1AWzUajaasWLxiE/FE0rIsnkiyeMWmErWoFy/RMqPSFjtCiAjwBeAVw48uhBBAA/BSepdlwBnpqJljgF1SyncK0nqNRtOvaWqJcdyiZzm4cTnHLXqWppZYqZtk4e3WuK/lxcSLW2Y0cG/a714BPCilfEII8awQYhQpN8w64Pvp7f8MfBXYDLQDZwffbI1G099paolx5SMbeizjWGucKx/ZAFByl4fBmGiEmI2Qj4lGStAaK67iLqV8EZhis3y6YnsJXJB/0zQazUDGyeVRLuI+b8Z4ywMIIBIOMW/G+BK2KkVZpPzVaDT9k6aWGItXbOLt1jhjohHmzRjvWZjL2eVhYPQl1z4WEi3uGo2mIOTrVilnl4eZhil1ZSHmmejcMhpNH6dcBx3zjSSZN2M8kXDIsqxcXB59AW25azR9mHIedMzXrVLOLo++gBZ3jaYPU86DjkG4VcrV5dEX0G4ZjaYPU86DjtqtUlq0uGs0fRiVFVwOg44NU+pYOGsSddEIAqiLRlg4a5K2xIuEdstoNH2Yco6zBu1WKSVa3DWaPowedNSo0OKu0fRxSm0d5zNRSVM4tLhrNJqcKedQzIGOFneNJkfsLNZoqRtVZMo5FHOgo6NlNJocMCzWWGscSa/F2hpPlLppRaWcQzEHOlrcNZocUFms23ftLVGLSkM5h2KWLVLCtm0FP40Wd40mB1SWaWeyu8gtKS16opIPpITly+Goo+DYY6Gjo6Cn0z53Tcko1ygLL+1STa2vDGXbS+XYz6DapEMxPSAl/PnPMH8+rFnTu/y3v4Xvf1+5W75ocdeUhHKNsvDaLtXkodrhlTkdr5gE3aZihWKW40PSEZWoAwweDDt2FPT02i2jKQnlWljYa7tUU+ujkXBOxysm5dgmN1QD2OWS3tiClPDEEyn3y9e/bhX2wYPhBz+A11+Hq64qaDO05a4pCeUaZeGnXXYWa3Pzqzkfr5CYrV6p2KbU196JPhFyafjU58+HtWut6wYPTrlgLr8cxowpSnO05a4pCeUaZRF0u8qhn5lWr4pSX3snyuUhaYthqX/603D88VZhHzwYLr44ZanffnvRhB20uGtKRDGjLPxUKgq6XdMmjPK1vBDYWb2ZlHuESzk8JLOQEh5/3F7UhwwpmagbuIq7EGKIEOJ5IcR6IcRGIcSC9PKDhRD/FEK8KoRYKoSoTC8fnP5/c3r9uMJ2QdMXKVY6WL++2qDbtfKV930tLwRO1m1fScVbViGXhqjX18PMmdmiPmdOSUXdwIvPvQOYLqVsE0KEgeeEEE8ClwC3SSmXCCF+AZwD3JX+vVNK+XEhxKnAfwOzC9R+TR+mGFEWufhqg2xXObgTVGGbddEIqxqnF60d+VAWIZeGqC9YAC+8YF03ZEivT3306OK1yQFXcZdSSqAt/W84/SOB6cA308vvBeaTEvcT0n8DPATcKYQQ6eNoNL7JJwSu1OIaRKm5fCn3nO9eKVn2S0PU58+HlhbruiFD4L/+C+bNKxtRNxBeNFcIEQLWAh8HfgYsBlZLKT+eXj8WeFJK+SkhxEvAl6WU29LrXgOOllJ+kHHM84DzAGpra49csmSJ7bnb2tqorq7OsXt9g4HQR8itn63xBLGdcbpNn9MKIairiWSFHdqx6d3dtrNGBQKJpDJUQe3wIZ6O5YXMPubbfvNxtu/aS2eyO6c257u/mYHweW1ra6N66FBG/v3vjLv3Xoa9ao2CSlZW8vbMmWw97TQ6R4woUSth2rRpa6WU9XbrPIVCSimTwGQhRBR4FPik3Wbp38JhnfmYdwN3A9TX18upU6fanru5uRnVuv7CQOgj5NbP4xY9S6w1lLW8LhpiVaP7sVozJuzYEQknWTjrsECsQnMfm1piLF69iVhriJAQJKWkLgd3QlNLjCuf2UA8UYExTBZkm/3S7z+vUrLhppuY9MgjSks9dPnljN1/f8aWpoWe8BXnLqVsFUI0A8cAUSHEICllF3AA8HZ6s23AWGCbEGIQMBwo7FQsTcHJ1TVi3q9xcjetLTFfgpSvWyXTV1uRFlkzhYiXzpwFmpSyxxXi9zx9Isa7PyAlLFsGCxYwKVPUI5Fe98v++5emfT7xEi0zKm2xI4SIAF8A/g2sBE5Ob3Ym8Fj672Xp/0mvf1b72/s2uc4OzNyvM9nte1ZhECFwDVPqWNU4nS2LvmZxj5gJ2gcf5CzQUo8b9HukhKYmOOIIaGiwWuuRCFxySSr65ZZb+oywg7c499HASiHEi8C/gKellE8AVwCXCCE2AyOB36S3/w0wMr38EqAx+GZrikmuQhWEwAUdAleseOkgBTlaZe8bL+dJR30Cs6ifeCKsW9ezKjl4MFx6KWzZ0udE3cBLtMyLwBSb5a8DR9ks3wt8I5DWacqCXIUqCIELOgSuWJEjQUXJNLXEaNvblbU8HBJ9LtqlbJASHnssFdJoEnQgZamffz6rjz2W42bNKk37AkLnltG4kqtQBSVwQYbAFSteOqiHyOIVm0h0Z7uShlYO0v52v3R394r6+vXWdWlRZ948qK0l0dxckiYGiRZ3jStehMpuwLXU8dWqQeBixEsH9RBRveXsGmDl/PLCTdQvuAAuuwxqa0vTvgKhxV3jiptQqfKDL5w1iYWzJvXsVxmqKNo093LIox7EQ6QcJkH1Wbq7Uz71BQvgxRet6/qxqBtocR/geA1xdBIqp4HTVY3Te/Zrbm5magFyx9i1v9jhg0Y7Th27m6sXPevJSvdy7Uv99lMq8irM4STqVVW9or7ffsE3vIzQ4j6ACcq6LVWonlP7i9kmSzvGeruOXq99ocYIyrmqUc6fSy3qFnTK3wFMULHY+YQX+knHm4lT+4uZIjaX6+h1n0KIcLlXNfJ9Pbu74eGHYfJkOOkkq7BXVaUGSbdsgZtvHjDCDtpyH9AEZd3m6jpQWWhr3tzBylfedxU0p/bfNnty0dwZfq6jIdZ2fvTMfZyuz/IX32Fne2pQNRoJM3/mRM+i78dlZfdwiXo6S+54vp7d3fDooylLfcMG67qqKrjwwlSs+gASdDNa3AcwQYYqgn/XgUpk7l/9Vk8yItUreVNLzDaVgNH+YqaI9XodM8VadSwD1fW5b/VblmWt8QTz/pSKAvHSP6/ieU3TBtt7sfDY7Fw/QeJ6Pbu74ZFHUqL+0kvWjQxRv+wyGFW8gijliBb3AUyQg3W5RIaoRCZTrjOtSkMk7YTd3P5ChTxmWrPTJozi4bUx5XV0s9bt9gF/b1CJbul5sNjLw6ipJWYRdoN4Isn2XYUNw1R+Lr94KDz0kL2oDx3aa6kPcFE30OI+gCl1AQSVyNhhFjpV2biQEL5CLXPxZ9u5Sh5eG+OkI+vS1ZV2WzI/erHWAdtskX6uD3h/GHh5qC9esUlZbzUzhXLQ4wKZn8u6fQZza/g1jjprnhZ1H2hxH+CUrAAC9iIjsMkPTUro3Czgbil9CXsuERkqV8nKV95nVeN0mpubuej0qY7bZxKNhG0rIvm5PuDdnebloe70oKgM9cZhFGo+QcOUOhoOH50aKL2+UYt6Dmhx15QMO5FRuTimTRjly1/tRq5x8H4Hob1Y03s6u2iySYWsuj5L/7WVRNIq8eEKf7lm3B7qqrcGAdQOH9Lzf0HmExjRLwsWwMaN1nVDh8JFF6VEfd99czv+AEGLu6ak2IlM/UEjsqxKNws4HBJMmzCKyQueojU9Nb+mKsx1x9tHkeQqxn4Hob24VhJJtb9cdX0WPL4x52gZL6jeGk4/5kCikQ97lgU6n0CLeqBocdeUHXaCNnfpOsXWKZLdkj+ufguzN3hne4J5D9lHkeQaKeR3ENpuezv8Zsos1CQm8/JoVZjBgyrYFU9Ytmk2JdUKJOKquzs1UHr99dmiXl2dEvVLLtGi7hMt7ppAKPSMRzcL2CZxIqC2inONFPI7CJ25Pdj7zIcHVMPVC07x82aX2M72BJFwiNtmT1b2L6+IKy3qBUWLuyZvipGky6sFbIedVZxPpJBfy9m8/ZTrn+pxp5gRdpWHA8T88FWVGnzgn1t9lyDM6Tomk72i/vLL1nVa1ANDi7smb9wG1ZpaYsxftpFzPh7nrMbljr5wFWYR8RMeCGoXQb7ujVxmb7baCLtqeVBvQ3b1XO1QLXdzGamuY1b7v/BxGl77h1rUf/CDlKiPHOnYl3LNiVNuaHHX5IT5S6YKzXu7NU5TS4x5f1pvKTjh5At3whCRppYYc1x88AaFqlikTHPsMnsz19msXt+G7MTPSzgmpOYJqGb8+sXc/oruJFP+sYJP3XI2fGCdXetV1DOPCaVJ49yX0InDNL7JTDylYkw0oqwkZPjCc6FhSh1RhY/a7N6oqQqz+OTDC/LFV72tbN+113E/rzVhc0lGpkoI5uVNJxIOcdrRYwOrV7t4xSY6Ojo5/uX/YcU9F3Lnspv5uFnYq6vh6qvhjTfgpptchd04ZlBFxwcC2nLvg5T61dSLJWiIglOUSz7pd+fPnGg7kFesYiCqtmfO3szEq486lxBDlfipLPKQEHRLaWmDXRiq7+uZTHLE3//CD/6+hEM/3GpZ1VYZoXreJTB3ridBN1Oq1NJ9FS3uBaJQApzPq2lQbXL6MgnIKpqhshzzSb9brqkTzLM3VXjx9auOL4HjFMVAVPclKSWRcMjTgzCvcYhkEh58EK6/np++8opl1e7KCL87ciZ//vxsnrzhxJwOr6tS+cP1kyiEGCuEWCmE+LcQYqMQ4uL08vlCiJgQYl3656umfa4UQmwWQmwSQswoZAfKkULmy8711TTINqm+THXRCFsWfc1SfWnejPGEK7JDQYLwhTdMqWNV4/SscxYDlXvFPHsz6OMbqO6d031ZOGsSddEIwvR/YNcrmYQHHoBPfQq++U0wCfvuygh3fGY2n/3+Pfz882fxvVlH5Xwary4tTQovlnsXcKmU8gUhxDBgrRDi6fS626SUPzZvLIQ4DDgVmAiMAf4qhPiElNJ/DFsfpZAl3nJ9NQ2yTaqwxHabafTG3/OXbST1UeqdOQopKzQXyzvoNyO/x1O9OUR3vZpzG1THt7NW7e6dU8x5QXIImSx1Mix1hg1j0+zvcPHoqWzqDBckoZiOlnHGVdyllO8A76T/3i2E+DfgdDVPAJZIKTuALUKIzcBRwD8CaG+foJC+wVxfTYNsk1mwjan+kIqCUZWKM2Y2vpFOqpWveynIqIlcj2cnmM3NwYi7+fgHNy63HbjOvHdFE79kEv74R7jhBltRZ84cmDOH8SNG8Jdgz1zSRHd9DSEVsa22GwsxDvgb8CngEuAs4CNgDSnrfqcQ4k5gtZTyvvQ+vwGelFI+lHGs84DzAGpra49csmSJ7Tnb2tqorq721alSs+nd3bYDa5WhCsbvPyxruZ8+tsYTxHbG6TbdtwohqKuJKCNIcmmTF/wes/Wj3WxvTw06CgTSRrK8tEd1XoHggBHO18HP8Yz21A4f4vmYXu9lazzB9l176Ux2u56jEPcuJ5JJ9lu5kgN/9zuqY1aXUFdVFdtOOoltJ59M1z77FK9NBaKv6M60adPWSinr7dZ5HlAVQlQDDwNzpJQfCSHuAm4gNcZzA3AL8B1SY2qZZH2LpZR3A3cD1NfXy6lTp9qet7m5GdW6cqXVJoe3MYA11cbq8NvHXFwSftvkhbMblyNthm0EsGXR1Kw2x95ey+IXQzgN9djt6/W8AJFwkoWzDvM2qSZ93ZyO53bMTLzcy6aWGPOeWk+iuwLjWoQrEiz+hv05CnHv/NC05i3WL/4F3/rrH/jYjm3WlfvsAxdfzKA5cxg3YgTjCt6a4tAXdScTT+IuhAiTEvb7pZSPAEgpt5vW/wp4Iv3vNmCsafcDgLcDaW0fodCvx7m8mhaiTX5cRItXbOLUse5viSr3ktv0eQOnWqAq14tb3pqgxksM5i/bmBX7n+iWXPnIi8rMkFACX3MyyZpFP+c/bruZhg+top6oHkb4krkpF0xNTWHb0Y8oZhizq7gLIQTwG+DfUspbTctHp/3xACcCRjb9ZcAfhRC3khpQPRR4PtBW9wHK0Tfo1KZcPnR+kka93Rq3PvJtUO3rdfq85VwZOA0oe8lbE2QstXmcwtqebtu87lDkz1MyCUuWwA03UL/JGoX1UWUVr3zteK6ddBp/WXBCcdrTTyj2DFsvlvtxwLeBDUIIY0bKVcBpQojJpFwubwDfA5BSbhRCPAi8TCo84oKBFCnTF8lnQBG8WZQpi3x31nK7iTSZeJ0+bz2XFacBZS95a4oVSx3kG4Jvurp6RJ3/+z/Lqo8qq7jn0ydwT/0JnPvpIWzaoKfI+KWQUXR2eImWeQ57P/qfHfa5Cbgpj3Zpikg+HzqvFuW8GeOJ/XutZZnXGaVuk6bMdrzK+ndzIZnz1gRVNFxFTVXYNjMklGi2pYOotw0eyq/rZ3JP/Ql8NMQYYOzK+WFX6tnVpaTYM2z141dTlA9dw5Q6mt59mbpoyPGLnVksQkp1zdBMYRfASUfaP2xULqRpE0ZlxdovnDWpoAJ03fETlYnPijrbsqsrNfnohhvg1YwQzuHDYc4c/mfaN/jlX7darluFyG0C2kBP/FXsGbZa3DVF+9ClCkFPVa7P/PKrrFuwLxQtgZWvvG+7vZ0LKbNea09mx1mTbAtWB0XDlDrWvLmD+1e/5emtwwu+LGI3UZ87Fy6+GKJRvgYk9olajl1Xk8xJjIvplijHN4S8CpvkgBZ3DdMmjOK+1W/ZLi8mXn3rdQ7RLbHWOAc3Lrf9Qme6kI5b9GxOYhOEcNzYMCmYJF34sIh9iLqZzOtmLrPnh2K5Jcr1DaHYUU9a3DVKa1e1XEW+ouflSy6AVY3TOW7Rs0qBN+fOAfUXOhexcROO1njCc0qFoCJgXC3irq7eGaWbN1t3dhD1oCnWG2KxBy79UMyoJ53PXROIReWWmKypJcamd3dzcONyjlv0rG3CMi9fcmMbp8RaBm4J1VTnc2qHk3A0tcSI7YwXJGGcE6r7tH1HG/z+9/DJT8KZZ1qFffhwWLAglU/9uusKLuxQvMRfOjVwCi3umpxELhM30bvykQ10JrsdRc9NsM1C0DClzpLpUIXTFzoXsXESjsUrNlnSQkBxiklk3qdQd5KTNjzDynvOzxb1aLRX1K+9tiiibpB5zwLPTpkmiM9zf0C7ZTSBDPS4iZ6X1+RMn6QRLbMrnugZAF28YhNzl67rcXkYA58qN020Kqx0k+TiA3VyLagmahXaYjTuX2dHJw0bm7no70sY1/qOdaNoNFXK7gc/SFntJaIYboliD1yWK1rcNYEM9LiKng12y1Vffjdft90XOhwStO3t6om6sfPD+xUbJ+FIWejZE7UKbTE2TKpl7ON/ovbnizlgR0amjzIR9WKiUwOn0OKuAfK3qNxEL9+BNDfr3+4LvaejK2uqfzyRZM7SdT1pB4LO0WM3UatgFmNXF9x3H9x4I0e+9pp13QAUdTPlmP6j2Ghx1wSCm+ilLOaunu39iF5TS0wZGWO2/jO/0Ac3Llce0ymaxi3qx044jH1OHSt7apbWBVhEJNYa7znugcMqub3rJY74w8/ATtQvvRQuumhAirqmFy3umh7yCWV02tf4vX3TC1k1Vr0c1xBhO5ys/1yyPeYSI23ZZ2xvzdIghN3yNpTs4uSNK7nw70uzfeo1NSlLXYu6Jo0Wdw0QfGWkuUvXMWfpOov12rzrVddc7Zk4TWxys/5zyfaYS4x0EHHVdg9H47iDkl2cuHElF/5jKQe1vmvdsaam11LvB0UyNMGhxV0DuBfedrLo7fY1AgLND4lcgu6cIk3cwuhyyfaYS4x0rpOhzFFBbXu7enK8G9cssbeDbyhEvXVINb/+9Ilc9sSdWtQ1tmhxz5FyzF2RT5tUQmQIjZNF7xbqZzwkbjrG/7QKlWulLhrx1De/2R5zmUXpdx+3HDqDkl18/cWnuOjvSzlw13bLutYh1fzq0ydy75HHM7x2JJdpYdco0JOYcsBtNmZfbFO0yr5+Z4XA0aIHb1EvucZ6qyY2tXd2+breXifQ5DKxye8+KlfToGQX33jxKZ791fdY/OQdFmHfOWQYN3/uDD77/Xv42bGzSVYPG3Bx2xp/aMs9B8oxd0W+bVIVN+pWLDeLtRffttewR7u3j4WzJjF/2UZLWOPO9oTvZFBewuNyiZE27wO7XaNkMh90KZ/6s7aWeufwKL89ehY/PezLxIcMDSwKR9P/0eKeA+WYu8KpTV7cNbsUpd9UmMU607etLKCxKyMTYQaqQd2FsyYxdPAg25j1QjxQc4mRNvZpbm7motOnOm5ruHEGJbuY9dKzXPgPG/dLZBixcy5g4o+u4nvDhqXKnGk0PtDingPFTrrvBVWbhkfCnqJgVPtHI2E6urpdfdVmQVQ9TJqbncXd6e0j6AdqKcdMLp9+CP+6/id877kHGJsh6juHDGPpf57CAdfM4+uf1W4XTe5occ+BUueusBMmVZuEg8/cLGaq/efPnAj0WuUhISw+dztBzHV2oJOAB/lALVm+70QC7r2XE266iRPeeMOyqnN4lMorLqfmwgv5/rBhhWuDZsCgB1RzoFjZ7exQDZwCtm1q9Vir06lPRu6WSDhEUlrD9YIcRHbK5hdkuli3sM/A6eyEX/8aPvEJOPfcVEZGg5EjYeFCKre+BVdeCVrYNQGhLfccKVXuCidhWtU4PatNfvK6OPVp/rKNnt4A8nF3OL0RBZEMyjyV3w6vLh7PfezsTOVTv+kmq6BDStQvuwwuuEALuqYguIq7EGIs8Htgf6AbuFtK+RMhxAhgKTAOeAM4RUq5UwghgJ8AXwXagbOklC8Upvn9Cy+i4RSP3tQSsw3ts8uWuKejS1mOzq5dmYOZdu3x4+5w6qvd8nx95HZx7pl4cfF46mNnJ9x7b0rU33zTeoCRI2HePDj//LIT9XKcu6HJHS+WexdwqZTyBSHEMGCtEOJp4CzgGSnlIiFEI9AIXAF8BTg0/XM0cFf6t8YBr8LolC/FbvvMSJYKAYmk7BFrL/5mJ3fF8EjYsp1X696pr3YPgnx95G71Wb26eBz7OHEUo594As46Sy3qF1wA1dWe2lxMyrXuqCZ3XH3uUsp3DMtbSrkb+DdQB5wA3Jve7F6gIf33CcDvZYrVQFQIMTrwlvczvPqBnaoVqfzGZp+5Xdy6m7/ZyV1h1Axtaol5jmjx6/N2q/J03KJnHcv3ufUBYEjY/quQeXy7B2s4meA/mx+FQw9l/C23WIV95EhYtCjllrniirIUdijBOISm4Aipmr1it7EQ44C/AZ8C3pJSRk3rdkopa4QQTwCLpJTPpZc/A1whpVyTcazzgPMAamtrj1yyZIntOdva2qgu0y9EULS1tbFll9qqnFRnzfLXGk+wdUe75+0BNr27m85kt2M77Pbzum+FSBW7yywzBxCqEIyJRuja28476mYr27AhtsvxvOZzVghBXU2EaCRMazzB9l176Ux2IxBInD/r5n0hdZ1jO+O2fQKo6Erwyb89S/1jD7HPh9Zi4p3Dh7N19mzebmggGSn/8m5O1zjzngyU72Rf6OO0adPWSinr7dZ5HlAVQlQDDwNzpJQfCaGsXGm3IuvbIaW8G7gboL6+Xk6dOtX2YM3NzajW9Ream5tZ8lK3MoeK3aQYlRWp2v7sxuVIhxc11X4ArS0x5i5d5yKNqVQFdm8GQytDdMsk50+AWzaoP3KqNlyt6KuR3zz7OCHmzTiUK5/ZQDxRgZ+gsLpoiFWNqTakrnH2W1I4meDkDc9wwT+WcsBH2aJeedVVVJ5/Ph+rruZjns9cWlTX2O6eDJTvZF/vo6dPvRAiTErY75dSPpJevN1wt6R/v5devg1rJckDgIzaX5pM/Ib6+d3eabDQzd/cMKWO04850LEQNahTFezpTDr6u93aoOqrnbCDum4rpB4IXgtqZ7pywskE31z3JCvvPo+FK+60Cvu++8LNN/PPBx6Ayy+H6mrPLqNyIMhQU0154Cru6eiX3wD/llLealq1DDgz/feZwGOm5WeIFMcAu6SUGZUFNJn4jZ03bw9YJhfZiYjKVx+NhD3F6N/YMPV/AWIAACAASURBVInTjzmQkPqNzXGdEyEhOOlI59DSwYN6P6o1VWFL3zNxqtvaLSVbFn3Ncd/Mv82i/qMVP7MVdbZsgXnzelww5ZhczolSzt3QFAYvbpnjgG8DG4QQ69LLrgIWAQ8KIc4B3gK+kV73Z1JhkJtJhUKeHWiL+xG9pdl2c/WiZ5k3YzyrGqd73t9cws4tyiHfOPGmlhgPr40preVIOMRJR9bx8NpYVpz64EEVylBKSFUuenhtjPqDRrhGygDsTaT8//nUbfUyy/iKaeN4YcFtnPvcUup2W90vH0b2YcnUU7ngT7fC0KFZ5ynH5HJu6Lqj/QtXcU8PjKpMss/bbC+BC/JsV78nszRbrqFnfkQkny+vUyihOUth/UEjsh4gkF1DNRNVm90mbRnbqOq2qsTb8WHX0QG//S0zf/QjZm7dajn3B1XD+eVRJ/HwUcdz7amfthV2KM/kcpqBhZ6hWiL8iLLT5JJ8Z1t6RXU8AZa3DXNxjMUrNjF36TrGRCOcdGQdlXvfyMoY6XYON5FUPbC8vKlk7dvRAb/4BfzoR5Ah6h01I/nVMSfz8wlfpGa/Gq51eeuJVoWzinAYy0FPGNIUHi3uJcKrZec0uQRQimXQGSr9JO6ya/PDa2MsPHYIWxZ9URnpY3esfBKGeX5T6ehg3YJbGf3z26jdZXW/sN9+cPnlDP7+97lw6FAudD8aoM6PL6WeMKQpDjpxWIFRRUw4Jcky42ThL16xyVbYBQQe5aCKppg2YVRW/1Rt3r5rr/JYAHs6sqsr5RLF4TlKpaMD7rqL9oMOZvLCqyzC/sHQKC/N/SG8/nqqALXC/aJClR9/VzyhJwxpioK23AuIykJb8+YO9nRk+5/tRCsX360keAuwYUoda97cwQP/3EpSSkJCcMSBwy0DqJn1VjMxJkIZbVvw+EaL66I1nl1dye9AsCeruKMDfvMbWLgQtm2jyrT/+1VRfnn0LO6f/FUGDx/G0J/+MyfXidMbR6H98drlowEt7gVFZaHdv/qtLIu7pirMdcdPzPoSurklVBNPDIL6omdGyySlZNVrO7K2iyeSyslFlaHeF8WGKXUsXrEpyy9tN+7gZyDYcSzjsH0tom7m/aoovzj6JO6f8hX2hoek9osnfOXgMZNPJE8+aJePxkCLewFRWWJ2rpSqykG22Q+nTRjF0n9tJZHs3SscEpYoFFVESJBfdLfEW2aSUhIJh7LaVTu80rJdvhas+VpFq8JIiW3IZWVXgunPLoefnZkl6uy3H3fUz+LnE77QI+oq/IQyOr1xrHlzR9YDPqgJQ30xBFNTGLS4FxCnDI6ZGLVOM8V46fNbs63g9L9uLosgv+h+XAZGaGRmu6IZNVTzGSzNvFZ2kSmDuzo55cWnOP8ff2J024fWlbW1qURe3/seB27aiXhkA3h4eHmtSQvqDJcPr41ZhF2A6yQur+gQTI2BFvcCYvdq7hTdYifGCZs5/Ylu2SPQTi6LIL/oXh9UxmCuXbsya6iqrk+sNc5x6UldTmmIVW8SXkWdqpS3vWFKVc8xDcFu7+yyfWDY1aSdu3Qda97cwY0Nk5TXxandElj5yvv2O/ikHOv7akqDFvcCYmdZT5swKmsWZ4VIuVnmLl2nOlQWTgJtWJaqePJcvuh2QmyHn8HczFzz5gefW6EPOwFzE/UNp3+fi4cfxZbtkjF3rLY8PDIfRnYzY1U1aSVw/+q3qD9oBFGsZFr5hZ6XUOr6vpryQYt7Hnh5PbezYDNncdbVJHsGGL26cVQC7VZxyM0yVvVJFeGSiSpniwrj2Hax706FPswM7upk9voVnL/6T+zfljHIW1vLhm/9F+cOnsK7yRC0W2vAGm2waxdku7xUD2CZ3vamY3oHje3cbIWelxBEOUJN/0CLe47kM1iZKfjNzc2AoiRehYB09SQDJ0vMyV3hZhl7qZBkDPrOX7Yxa/AyHwsxl0IfTqL+fnUN737vYt6YdTqXP7nZ9pq4jT/YPZjda7AOtWxrZ+VnCnzQlrXOEaMBPYkpZwoxEcUuM9/ibxzO4pMP95StT+WuMMi0GDPb66dPQwen7AIjE2S+WQS9Tup6uzXO4K5Ozlz7OP/zy+9y/V9/aRH294bWcPvXz+cfTz3PpB9fx6K/veXoSrKbEew0AWrejPHKREt2bbVDgs6+qCk42nLPkUJFJbjlSlFh567wgrm9qgeDeXmmdW+EPeb76u/JV7x3Lz94+S+ctvKPWZb6B8NGsO8N17Lfeecxx1T5yO1+mAXZadLZylfe73FzHPuxEfz9tR321rcpIkjlY6+LRnxl/9RockGLe46UW1SCkzvGKe2uub2qyUfmPO0q637B4xvz8vM6+or37oVf/QoWLmTuO9bSAO8NreHXx53Cp+ZfxszPfNy2f6qHVubDw8uks1hrnB17Ojn9mAMtgm+0tenJlzlu0bO83RpneCRMOCQ8u9Q0miDR4p4jbpZmsaeAO1moC2elQvTsBlrbO1P5XBqm1ClztZuXq86zsz3RM9Ca62SprLeWeBzuuCNVYDpD1D8YNpKfHXUSzf/ZwMVf/w9mKs6jivKxmxHsddJZPJFk5SvvZ1nfTS0xYjvjPaX5jIfp0MoQ7Z1JPbipKSpa3HPEydIsxRRwJwt18YpNzJsxnoWzJmUNhO5s783nUufgRnALr8wkr1mR8XjKUrcRdUaPhsZG9j33XK6LRLjO5VB+okf8TjrLJFV4JfsKtXcmuW32ZC3qmqKixT0PVP7xUkwBd4pDNx4uC2dNYujgQVnuGaNtqreRaRNGeYpxz8T3+EM8DnffnRL1d9+1rhs9Gq68Er77XYjkFm7pht9JZ5m83Rq3Vg9OY4RJanHXFBMt7h7w62IpxRTwzAlBmRgC7tQ2lZXrJ6+MGc/jD15E/dxzYYhz7pd88TrpTOU3T/V3t+2x9fR/TbHR4u5CLi6W4ZGw7eDl8Ei4cA2l10I9uHG5rbXpNEuyQoge33tmv/zMnDXwNHBYJqJuxtx/c256Y7C5zuHhPm/GeLa+vMb2uHr6v6bY6Dh3F3KJZxeKQGjV8qBRPUQMS9yuUEZSSq58ZINtYQu/wlRTFXaO3Y7H4fbb4ZBDYM4cq7CPGQM//WmqSMZFFxVV2M0YD3XjQegl5LNhSh0jh1ZmxcHrCBlNKdDi7oLqddqYwm832aVVMT1ftdwLXqsLNbXE2NOZXQikglRkzNyl6xgSrrB90KgeWnYPhHBI/aQy0hdntaGjo1fU5861F/XXXoMLLyyZqBvkOkltTDTCbbMn60lKmpLj6pYRQtwDfB14T0r5qfSy+cC5gJHK7iop5Z/T664EzgGSwA+klCsK0O6ioXJjGDlaINtVE3QMvB/X0OIVmyxx1Qbd9KbFdcoNYzy0Yq1xiyvipCPrsuK65yjcNVkPxPZ2+OUvOeaGG2DnTuu6urqU++Wcc2DIEK5p2mCp9nTa0WM9ZVvMBaexlHzGTbwO4OqKSZpC4sVy/x3wZZvlt0kpJ6d/DGE/DDgVmJje5+dCiGwfQB/Czmq1i6AwW3W51P10wo8Vme/AnfmhZcS3GwWu580Yz5ZFX2NV43QaptQpk4T1PMTa2+G221KW+iWXUGkW9ro6uPNO2LwZLrigR9jvW/2WpdrTfavf4pom/zNv3TC7XSS9D0y/NW4LdX6NJl9cxV1K+Tcgu56aPScAS6SUHVLKLcBm4Kg82ldy7PK9qGK9DWG12yefV3M/VmS+4qPqm93DRPUQa/zcgRZRZ/v23g1sRN3ggX9utT23ank+uD0w580Yn0raZiJcIQLznesi2ZpCI6RiVqJlIyHGAU9kuGXOAj4C1gCXSil3CiHuBFZLKe9Lb/cb4Ekp5UM2xzwPOA+gtrb2yCVLltieu62tjerqar/9Kiib3t3dU+zZTGWogvH7D/N9PLc+qs43qEJQIQSdyW4qQxXUDk8JZWxnnG4P9zUXJtUNt/zfGk+wfddeOpPdVHUlOOq5pxn/8J+sVjrQse++/N/JJ7PzxBPprrSW2zPYENvl+bx+MLfRuE5bd7Q7nqs1nmDbzjjm74cQggNqIkQdop68fl4L1ddiUI7fyaDpK32cNm3aWillvd26XEMh7wJuIGXo3QDcAnwHbBPm2aqMlPJu4G6A+vp6OXXqVNsTNTc3o1pXKloVhRwWzprE1Bysc7c+2p0vHBIgjUpNFek2JFNvDPu7x2rnQl00wkWn27SzvR1+8Qv47/+G997L2KkOrrqKweecw4f/+EdPP+38zbevWK9MgVC3tTsnn3RTS4wrn9lAPFGB+ToNHhSxDVc1+pgad8j2KNZFQ6xqnGp7ntQM1SRLXnJv69U2+evN5y9nyvE7GTT9oY85ibuUsuc9WwjxK+CJ9L/bsM7ROwB4O+fWlQi3ga5cCiLkM3hmd749HV3KmaaGT9yMuUDIkHAF8UT2m4ATtmMGe/akRP3mm7NF/YAD4Kqr4DvfgcGDaWqJsf3d3ZzduJxoVZi2vV09JQQNf/Mxh9Sw6jV7D2CsNc6cpetY8PjGrJwwTqjcH0PCFbZFvI0++nGFWQa8x6qzSZrvua6YpCk0OYm7EGK0lNJI+nEi8FL672XAH4UQtwJjgEOB5/NuZRHxGpnipyBCELlmMs93cONy2+1i6QLOThWh7KoeOZE1ccdF1Nd9+3zmVE3hzTeTjLltVc+bw/kTupFU2EbrxBNJ3vgwzreOObAnWsYOcy4cL9dOJdKt7Qlumz1Z+cBVRTxFq8I9WR+dZvDaZZPMLHwC3g2EYkTW6Oid/oWXUMgHgKnAvkKIbcB1wFQhxGRSLpc3gO8BSCk3CiEeBF4GuoALpJT5+QKKTCHywhTimE5JrtzEz0tEjeFmshzDg6W+7IgZXPrYKyTaeh9k961+y0OPUu26sWESNzZMYpzi4QWpa3fpg+uZu3Sdqwg5haU6PaBtq2KFBG17u7KyX6rcXaqIKlXdVhXFSERXimR3msLiKu5SytNsFv/GYfubgJvyaVQxKUYBY6+v+H4sJ6dEYW4PDlU/Q0LQLWX2uffsgbvuSon6++9bdxo7NuV+OftsGDyYaxc81eNu8Yuf9AzmME0nEcrV/ZGZqyckhO38AXNqAi/k8jkqRiK6UiS70xSWAZ1bplgFjFW5ZrxUAQJ1Aec1b+5QWsVOIqISPFtL3aOoG9j10yt7TLnla6rCjpOtzDiJUD4Fo41t3DJiGqkJvGSTHB7Jduu4taUYiehKkexOU1gGtLirChhnks9AlyodQGbMtF/LqaklxsNr1RNenB5GroLnJupXXw1nnWURdeOtIx8SSdnT3+uOn8i8h9bbWst2OIlQPgWjvWTErDP53mE3dYoIpXCFYE9n70C4V9dHMap+lVtlMU3+DGhxdxKEaCTMrngi74ElVTqA6iHW/Ct+LSe3snpe3A5ZfdqzB37+c1i8OEvU2/cfQ9X8a7NEHbLfOuwYWhmiMtSNwHm8wDwRDKxuESMlgZ0LpFAi5Ga5mpOJNUypo7m5uSeU0RyhNCYaob2zK+ttxBg/ALXAFyOyRkfv9D8GtLg7iczQwYNYd92X8j6HU7SGl7aoRMutrJ6fSJ6fPr6OL6x8mO//61Fq9rRa1seGjeJnx57CE0fM4PqjjqDBxlp3i7wJhwQ3nTiJ6K5X2bJoKqCO2DH73e0eQHYPkkKKkNNnxCn9L3iPcDIychr72B0HrA8682zWIHzi+bivNOXJgM4K6SQIQfkaveYo8ZuPRjX4WJeOAvHC46v+j82XXcuDi77Jlc2/tQj7tn1GceWMC5n6vbv54+Sv8FF3hcXtkpkSV0VdNMLikw+3tKmpJUa7jasKev3uKoJO7eCGKrdQLji9XbilHmiYUtfTlszB5KDy0TRMqWNV43RL/iBN32VAi7sxcGdHUK/5bqLd1BJj07u7e1LxRiNhV9G6pmmD7cCl59wnbW1w8818dsZRXPbsPYyMf9SzyhD1aefdzQOTv0wi1Ht9zA88r77oTJEwHgqqwdJEUjJn6TrHtMbFxPwwAetAqV9xVeXSN3AzKAqdj8ZrWmlN32BAu2UArjt+YkFf851edw2hM0/uiYRDjsWUm1piygiZLrcQxLa2Xp/6Bx9QY1q1bZ/9uPMzp/DwpM9bBN2M+YHnJkQCmDZhVNZyryX7VIONpYjHNtwrdq4kP+GCxjaXPmifZsHNoChkRIuOc+9/DHhxL4avURWtkUuEjDH4ZodEMYGprQ1+9jP48Y/hgw8s+9iJejQSpqOrO6ttezp6QxWdfNFGWx5eG6P+oBGeBo7tsLsWpYzHDkJcVeGVXgyKQka06Dj3/seAF3fIL1QuH3LJX+I2WcbyhTREffFi+PBD64YHHUTLt8/nzOQn+ai71zsXCYeYP3MiAAse32hxn7TGe6f+O02ism1Lev8KHxN+IPtaeLlmqslg+U6vD0pcczUoChnRouPc+x9a3EuIH7Hw6s4AaN3+Yaro9I9/bCvqXHMNnHEGUyorud5B8Bav2GQbumckJzO2eTtdcMIOQxyaWmLEdsZJymyfs2rCD2RfC7drdk3TBtucLmve3GGJO8/F7RCkuHoxKOweRgtnTSrIW6aOc+9/aHEvIYZYpNLwpFCJhRcLamhHO2e0LOfc5x8F0yApAOPGpSYfnXEGmPKpO4mMmzXnJRmZIQ6pdLjZEh4SgltOORzw5qpwEtimlphF2A3iiaRtMjK/bodihguqfOALZ03qebAGiY5z739ocS8hhihs3/RCz+QelVg4xuSbRH1EhqhvHV7LfZ//FoddcSEyHGbxrc95FiY/1pybOLzdGrcmg07TLWWWT93NnWJXz9UY8FS9AahcQX7dDsVy4RXbB67j3PsfA17cS53mtGFKHc2myT0q7MRzZHcHs//5GOc+/yg1e3dbtt86vJaffmY2j3xqOl2hQYQfe9lU3MObW8KPNWfkujEXtz7pyF4hTD0QdmftZ35QeJm0ZNRztQsTdRLqYs9szZdS+MBLNfakKQwDWtyDCv8K8gGhOpbZstq1/UMu/PcKvrP6YSp3WWeUZoq6gSqjoZMl6MeaM3LdmItbG7lvVr7yfuoNIMNy9/La78eCVb1pCOC0o8dm5XopZ7eD9oFr8mVAi3sQr75Bxge7HavhY8No6F4Nv78FdlgrFr01vJafHjubRydaRd0NN0vQqzWnupaqmHy3qftu7Yu1xrOyK9q9aQjg9GMO5MaGSVm5XsrZ7aB94Jp8GdDiHsSrbxAPiNZ4QjkgGU8k+dljL9Dw5L1wi1rUmyZOJxEaRJ2iBJ8KwxLM9+3DzzWLRsKOg4LmtqhCJwX0XC/zYKNTNElfcjtoH7gmXwa0uAfx6ptP3LWxLrYzbluMubqjnTNeeCIV/ZLhU38zuj93fmY2j06c1mOpG9P97ZJrWQtq99Le2cU1TRvyDhN0m9Rkxvzgybw2malyVcKuqnLUl3Ki2PU9c6C4EJExmoHBgBb3IF593R4Qbq4WuxDB6o52zlz7OOf+61Gie9usBz7kEOaNP94i6gbmEMXMwc3Znx5L/UEjmL9so0Vcd7YnbF0nft8+vExqMtPUEstqS6w1bhvKCNYqUYWollVs7D4X5vugp/9r8mXAJw7LN8OgW2Iwt2RPZkGq7mjngr8v5blffId5//sHq7Afcgjccw+88gp//9xMW7+6+YGiGtwcOjg4f7yZzARbTgytDHHlI/bJz1ShjN1S9mQrVJ2jLw02epmUFmRSMM3Ao19Z7rn4jfP1w7r5Rt3cNmOiEcLt27nw73/mu/9qyrLU9xxwEENvmA+nnw7hVO4Xt4k8dompDKHwI9gSGNe4nJAQnHb0WG5smOR6LRqm1HFw43KlSIdDgnCownc5PrNw94fBRq/3oS+9jWjKi34j7qXMauf0gHB023z0Eb/a+hcOufEuhuyxivpbI0ZzxzGzef7Yr/C56GhW3vK/WdPQzW6NIeGKnin2ThN2/PjGDZJS9rgM3ATeqc8CweKTD2fu0nWO+2f61DOFuxSDjUHPh/B6H/rS24imvHAVdyHEPcDXgfeklJ9KLxsBLAXGAW8Ap0gpdwohBPAT4KtAO3CWlPKFwjTdipv7I/OLabesEOJgZ2WOSu7lV1v/B8adzGE7d1q2f6tmNHcedyoPf3IqyYoQ7E7Y+mJPOrKOjq7unuU72xNKf7WB0U8/vnEzD/xzqydxV1nWB4yo7BlnUAlbJBxSzkA1U8zIl0IYDl7uQ197G9GUF14s998BdwK/Ny1rBJ6RUi4SQjSm/78C+ApwaPrnaOCu9O+C4xQPnfnFnPen9SB6J/YEZeU7WXeLV2xi9/YPuOilJznzn49SuXuXdeePfQx++EO+tW0Mb+3udDyPKleKk7ALyGpPrDXumLQrk8zzeZlwZV4X3fUqoBa2mqow1x0/sewGEFWGw6UPrmfu0nU5GQd218guWqbcroWm7+Aq7lLKvwkhxmUsPgGYmv77XqCZlLifAPxeSimB1UKIqBBitJTynaAarEL1mmvUmzSTGQ4I+eftcLTuDqmmIbEKfncrtFpnlPKxj/Hvk0/mkzfeCIMGsVVRZzMTP2lzoVfAzRN/bp89GcDRks7EqM5jF+limXBlY1k3N7/as944b18QMpXhkFnuDiDq47h9Ke5e0/fI1edeawi2lPIdIcR+6eV1wFbTdtvSywou7ipXgB/3g5/Bq6aWmCXfuZ0FPKjtIz64/GpY81i2qH/84/DDH8I3v8mmp1fy3R//zXHSTr5EI2FllkE/rpp5D623jZcH/1WJSpXiwS9e/ONG3286ZkAHoGnKCCE9CEnacn/C5HNvlVJGTet3SilrhBDLgYVSyufSy58BLpdSrrU55nnAeQC1tbVHLlmyxPbcbW1tVFdXO7avNZ5g+669dCa7EQgkkspQBbXDh/Qs90JlqILx+w9z3a41nmDbzjiqa1fZvofDVzzB5CeXMaR9j2Vde10db37727z3hS8gQyFa4wkS8XbeDSAoQiAQIhU2aKZCCCqEfRm+ylBKjLxeIy9Mqhtuu9zLvbSjNZ4gtjNu6VeFENTVRIgqCoUHid35VRw8PJRTH43zGJ9X4/NbjP75Jdf72JfoK32cNm3aWillvd26XC337Ya7RQgxGngvvXwb1vRQBwBv2x1ASnk3cDdAfX29nDp1qu2JmpubUa2DtDvkmQ3EExUYYfuRcKgnXt12tmaFsPjczftM9WANplIFZM8o3WdvG2etfZwz/tXE8A6rqHPoofDDH1J12ml8ctAgPmk61qlj4ZYN9rdClc0wc300EmZPZ1dWgrBoJMz8mROZu3Sdi289QItzwx7b3DFu91KF6nrXRUOsavR/vFzwkhKhLhrhprqKnPpo/zlOsnDWYWXnusn1PvYl+kMfcxX3ZcCZwKL078dMyy8UQiwhNZC6q9D+drfcLir/rt0yr1+iTPfNPnvbOHvNMs5Z8xj7ZIh624EHU33jAjjtNBiUfblVec4NklI6DnompSQSDiGEfebH3XtThUByCYHMhyBDUcuhBJz5s2RnMPREtqQHjZ2wczHpGqaaoPESCvkAqcHTfYUQ24DrSIn6g0KIc4C3gG+kN/8zqTDIzaRCIc8uQJstuH3xnXy1uX5pDKF0EvUtNWNYe8YFnPzjy21F3XwsuzznZtycAfFEUukvT0rZEzqZmfLWCeOB4ieaxq5dQYhTtCqcVe4Pgo8B9+rXdxoQNgaNnc5hN/ahui96EpMmV7xEy5ymWPV5m20lcEG+jfKD0yShQk1suurY/dl8zY84+/mmLFF/vWYMdxx3Go9/8nN0V4S47JoVyvS2TS0x9nR0UWiM1Ls1VWFAEk84+9fNgp7v0G6+bwtNLTHa9mZfo3BIBBoD7vezkmuki8pC72vFRDTlT58f2nfK7eI2sck3ra0wfz5fm3ksF//v/RZhf31EHXO+filf/O5dNE2cRrIilFWk2QgjhF4x8TsNX0U0Ek5lfnRgZ3uCvQ7CXheNUBeN+BJ0t3MKsPTbL4tXbLKNzBlaOShQd0XgnxUFTmGVTjmKNBq/9Hlxd0r+5VbowbPopEWdceNgwQLYZZqA9IlPwH338cVzft4j6nZkCoWXxFFeiYRDzJ85MTVQ7IKTcE+bMMqXG6CmKszikw93TBYmoaffRt76gxuXe77+qvbsCuih2NQSU+bSh9RnJZ+HUyYqS9z43OaTxE6jMdMvcsuoXpGdBhE9uWh27oTbb0/9fGQtPM0nPsGaMy5krpjAtg2dVIQEuITKmasH+XV3RCNhOrq6sx4IRjQMQLuLu8WNh9fGGB4Je3qbiIRDPbNJjes3TjEB6+20QJrz1nt1kRWy3JzdwKgdQeYomjZhlG2K5WkTRulJTZpA6fOWuxN2LhszytfunTvhuutSlvr111uFPW2pN93/NN/eeyhbP+pE4m3GqFE9yGnLaCRs+3o+f+bELMvu9tmTWXfdl3ryteRLPJFECByvmdGPIw4czvxlGxnXuJxxjcuZcv1TypjsMdEIi1dsyooT9+L2cEupnA9e356CdM+sfOV9X8s1mlzpF5a7isw8KnZYXvudLPXx4+Haa2n6xGdZ/NfNxDZs8N0eN/k3RNxos58In6CiKlrbE5x+zIG2uWsMJLDqNWu5v53tCUIVgnCFsPjIDSGeu3SdbcinlxquUJhUBX6uWVDXtxzCOjUDg34t7tDrslH5VcdEIylRv+02+MlPskT9tREH8IcvnsHky74HoRDz/rTedoDPTGbaAy/hhHUe4u+dQvWcXFBGJEZdNEJ7Z5dtWKFBVWXINbukimS3ZJ+qMFWVg7LamLJ8s0M+vbhXCuWu8BP7H1TUihc3UylTLWj6D/1e3A3s8qfs39XOb15/GsadmC3qI8fyk2NP5YkJn6W7IkTksZcB6SrshkjbTVKx+1LXRSNMqqvgotOnOobjAY6hevNmjLd98IRDqRzqThNwzNvu6cxvkLe1PUHLtV/KWj5vxnhikthpIgAAEFVJREFU/7ZmoSh1NIjdZ0I1ezmodroVGillXQJN/2LAiLv59X7PO+8xZ8NyTn++iXBGkQwmTOC6w0/iDwceRbcp8sWLb9b4kqoszcwvteGD3/RuN61pa80uHG/B4xv5KN6lrK5knM+cyMwgkZQ92zSZzmFY82arfk9Hl6fBVKc3EZWF2zCljqZ3X6YuGiobi7QQs5dzPad5uZ6pqgmCASPuAA0HRfjka8s44A+/YmhHu3XlhAlw7bVwyin8/uq/+HZLmN0q5rS6djnOM/Oodya7HaM2nNwoZl9tq2I7I1rFfA4jrtocbnewQ7phc4HqaRNGsfT5rbZvCU4WbjQSZlXj1J6HzNyl61i8YlNJRV71IC5ke5zcTNonrwmKgSHuO3bAbbeRuO12xmdY6pv3HcvOSxv59LyUTx3852GpqQqzqnG66yu1k/8/15h3s6Xs5M9VWYRzlq7j0gfXc9rRYx3K48EtpxxuEaT6g0ZYcrp7LbSh3Q7OFDL0UzOw6N/ivmMH3Hor3HEH7N6NOVDv1ZFjuePYU1k+4bOMpppVoV4XjJ/85uGQ4LrjeyNcvLxSB2mFmS1lVbvdBlGNGqnHfWwEO/Z0ZrmOTj/mwMDK3AXhdijUgGM5DGT2h+LfmvKgf4r7hx+mol/Som7m1ZFj+clxp/Hn8cf1+NQzxTbTL+rkojEPVjpNmDq4cXmPYASZodEsPsbfmVWSnITdzOrXd3LLKYfnJXBuAuk0a3jygqeYP9PZ+i+U5V/I4/q5nn2tSpWmfOlf4v7hh72WelvGQOlhh3Ht4bO474BPWwZKwf6V12yZqsIo66IRy5fOKfe6xFrcOjNDYz7ZFzPbPX/Zxpz2TUqZV9ihF4F0erC1xhOp+raoBbVQA46FOG6uDww9U1UTBP1jhuqHH8LVV6dmlP7oR1ZhP+wwWLIEXnyRI+Z9n8GDKy27ennldZolaeQmObhxuadZqvFEkpWvvG+ZbVoZquD0Yw702tseqsLW29fUEmPygqdyTkYWEu65aZzwknzLbdZwols6zgYt1IBjIY5brGRkGo0d/UPc77/fXtSXLoUNG2D2bAiFHJOMuTHEJKTRSJiFsyYBqfBGt5QCmcRa4zRMqWNV43S2LPoa4/cfxo0Nk9Ipeb2TSMqepFZesky6Sfcxh9T4On8mTi4Xo53GPcjlOKAeWMx3wLEQx9WRL5pS0j/E/dxzYf/9U3+bRf2UU6DC2kWzqK5qnO45usPst+7oSiXoyjWzoyoN7nXHT3TN62LGbOV6aYvb28EbH3oTHfPbijm7o5MQmh88DVPqHDNJOh2nULlmCnHcQj2INBov9A9xj0RSvnYHUc8VpxDCXAdFJdj6xTPfLLxgWIFu1mBNVZgbG3K3mA2Mh53xthJrjTNn6TpO/9U/mDZhlHK/eCLJ9l17e/5XbRuucI6Vz+fty4lCHLeQSc80Gjf6z4DqaaqCUf4xRzgEMchpR2s8QVNLzDHE0CnPuIFhBToNVBrpeSEl8qroGS8WpeoNYdVrO2h5q9Vx385k6o2nqSXG0ue32m4z+6ixroJaqAHHoI+rI180paR/WO4BkmmZ+sVLwQyDXNLdmjFbgapta6rCFgv0a/8x2vZYIReL2cDJunfLJ18ZSn3cVNWVoP+lvvXrBtRogkKLu4mmlhiXPrg+rwpJi79xON865kBPbhUv6W7NroKaqjDRSNjWbdAwpY6TjqzLinipqrS+nKnEc9hgb2XrcvUXR8IhaocPAZz7rQcb7VGNc2g0KvqPWyZPDIvdSzijE0aulPqDRvS8jgsBdoZqEOlur2na4Jh7PTO22mvZOtXkm3kzxjNn6TplezLzuUNvtajorlcBZxdSkIONQcw4LYdZqzplgyYX8rLchRBvCCE2CCHWCSHWpJeNEEI8LYR4Nf07v/i6AHGyftyiTbzGgJu/eMbr+K2nTC7IwNo1TRu4b/Vbrg8kc2y1SjwrhMgKqzQPmhoFvhum1HHcx0aoTyawvF0Y1aIANr27m4Mbl9Pe2WX7wXNLPOYHpz4U8xhBoOPlNbkQhFtmmpRyspSyPv1/I/CMlPJQ4Jn0/yXH7Yvq5A4QeCujZ5D5xStUhMf9NrU4VRj9U/nmk1Iy76H1lrTAZsx9uv/cz/AtRVhlIikZOniQxcfc1BJj3p/W05nsRpJOhyAgYpo7YBTbDsoSDUIQy0VUdby8JhcK4ZY5AZia/vteoBm4ogDn8YXb9HInV0Eujhq7fDVehMurG6CpJearXYbFbhzrkgfXZbmKEknp6HIx9+nGhknKik2ZfZ+/bGOWq6Zbwt5ENyLdtqDdHUEIYrmIqs4UqcmFfC13CTwlhFgrhDgvvaxWSvkOQPr3fnmeIxDcvqjzZoz3FekCKdeKalZpLl88P24AP9ajwJo9smFKne0YgBuZffI6SUc1a1ZCwdwdQUwgKpdJSDpeXpMLQuYxgCiEGCOlfFsIsR/wNHARsExKGTVts1NKmeV3Tz8MzgOora09csmSJbbnaGtro7q6Ouc2Gmx6d3dPnLWZylAF4/cfBsC/3/mILo+qJxAcMCL1Jd+2I4402bDGumjEWzoBo49e2miwIbbL07EBRg6tzBIkP/tDyidfV2PtU2s8QWxnnG7TZ8jYDmD7rr2W/tRGYLuD0WvXz1xxapvX+5LLMYL6vNq1xbielaEKaocP8dyPoClUH8uJvtLHadOmrTW5xC3k5ZaRUr6d/v2eEOJR4ChguxBitJTyHSHEaOA9xb53A3cD1NfXy6lTp9qeo7m5GdU6P7Ta1A41KhFNTbsDzm5c7tnVIYAti75IU0uMnzy33lJzM1W39LCe44Kzu8XoY+r82S9TArht9qGW/ds7I7aTkQYPqiAkemPOhYDTjz6Ab061zk6d4zHBmJvbxK5fYJQUrMD8cnjppC5u2aD+yKWu6VTXNnmlFNEyQX1eyxndx75BzuIuhBgKVEgpd6f//hJwPbAMOBNYlP79WBANzRcvswX95Fk3LOHFKzZZhB2sdUvBeyib6vzRqnDW/uEKQTgksgo5H3HgcFa9tqNnmZRwX3rg1Zx+YP7MibYFtc3URSOsapzueB3sxhKOW/RsTnMFgnZ3BDHjVKff1fRV8vG51wLPCSHWA88Dy6WUfyEl6l8UQrwKfDH9f1ngNlvQzrcZsvHDm/2dXgbdFjy+0VPUhcq3KmV2Gb5EtySRlD0hmkYEzurXd9q254F/Wqf7N0ypY/E3Du9J4JXZSy8+XVVoqVvkUU1VOGt8Q/uQNZpgydlyl1K+Dhxus/xD4PP5NKpUqKx7u2XGtlFFrhbDCm1qiSlzubhVgDLONdchgsUodG20SRXtYhfKabZK/bofnN5GVG8glaEKtiz6Wk7n02g0/tAzVDNQvYar/M1te7uylpsn4zhFtbhVgDJYvGKTo7vIHNKpqgblNgnLr/vBKbRUVQe0dnhvoRTt7tBoCovOLZMHqgRYQyt787Q4uSi8uiHcEoiZz3Pa0WNt16uW54qTO0o1aatU0R0azUBEW+554CVPi3KQNBL2bLma3TVuOVmMQVMj30xICE47eqxrLne/uE2ssbPMm5tfDbQNGo1GjRb3PPAyc1Dlopg/c6KvcxlimenrNo5nfgu4sWFS4GKeybwZ47OibdwKbWg0muKh3TJ54GXmYNB5ZQqVpyYnMt34+dXX1mg0AaIt9zzwWmmnEBV+Sj0Y6SW+X6PRlA4t7nlSDkJbCsolqZZGo7FHu2U0vmlqiVGhCK3UmQo1mvJAi7vGF04Vq/QsU42mfNDirvGFqmJVSIjSDexqNJostLhrfKHyqXdLqYVdoykjtLhrfFEuBSw0Go0zWtw1vtBVgTSavoEOhdT4wmtsv0ajKS1a3DW+Gaix/RpNX0K7ZTQajaYfosVdo9Fo+iFa3DUajaYfosVdo9Fo+iFa3DUajaYfosVdo9Fo+iFa3DUajaYfIqRNdr+iN0KI94E3Fav3BT4oYnNKwUDoIwyMfuo+9g/6Sh8PklKOsltRFuLuhBBijZSyvtTtKCQDoY8wMPqp+9g/6A991G4ZjUaj6YdocddoNJp+SF8Q97tL3YAiMBD6CAOjn7qP/YM+38ey97lrNBqNxj99wXLXaDQajU+0uGs0Gk0/pCzEXQgxQgjxtBDi1fTvGsV2fxFCtAohnshYfrAQ4p/p/ZcKISqL03Lv+OjjmeltXhVCnGla3iyE2CSEWJf+2a94rXdGCPHldNs2CyEabdYPTt+Xzen7NM607sr08k1CiBnFbLcfcu2jEGKcECJuum+/KHbb/eChn58TQrwghOgSQpycsc72s1tu5NnHpOleLiteq3NASlnyH+BmoDH9dyPw34rtPg8cDzyRsfxB4NT0378A/qvUfcqlj8AI4PX075r03zXpdc1Afan7YdPmEPAacAhQCawHDsvY5nzgF+m/TwWWpv8+LL39YODg9HFCpe5TwH0cB7xU6j4E2M9xwH8AvwdONi1XfnbL6SefPqbXtZW6D15/ysJyB04A7k3/fS/QYLeRlPIZYLd5mRBCANOBh9z2LzFe+jgDeFpKuUNKuRN4GvhykdqXK0cBm6WUr0spO4ElpPpqxtz3h4DPp+/bCcASKWWHlHILsDl9vHIjnz72JVz7KaV8Q0r5ItCdsW9f+ezm08c+RbmIe62U8h2A9G8/LoeRQKuUsiv9/zagHGvAeeljHbDV9H9mX36bfh38YRkJh1ubLduk79MuUvfNy77lQD59BDhYCNEihPgfIcT/K3Rj8yCf+9Gf7qUTQ4QQa4QQq4UQ5WhE9lC0GqpCiL8C+9usujrfQ9ssK0l8ZwB9dOrL6VLKmBBiGPAw8G1Sr42lxsv1V21TNvfOhXz6+A5woJTyQyHEkUCTEGKilPKjoBsZAPncj/50L504UEr5thDiEOBZIcQGKeVrAbUtUIom7lLKL6jWCSG2CyFGSynfEUKMBt7zcegPgKgQYlDaYjoAeDvP5uZEAH3cBkw1/X8AKV87UspY+vduIcQfSb1eloO4bwPGmv63u/7GNtuEEIOA4cAOj/uWAzn3UaYctR0AUsq1QojXgE8Aawreav/kcz+Un90yI6/PnJTy7fTv14UQzcAUUj78sqNc3DLLAGN0/UzgMa87pr88KwFjVNvX/kXESx9XAF8SQtSko2m+BKwQQgwSQuwLIIQIA18HXipCm73wL+DQdMRSJanBxMwoAnPfTwaeTd+3ZcCp6UiTg4FDgeeL1G4/5NxHIcQoIUQIIG3tHUpqsLEc8dJPFbaf3QK1Mx9y7mO6b4PTf+8LHAe8XLCW5kupR3TTI9AjgWeAV9O/R6SX1wO/Nm33v8D7QJzUE3hGevkhpERhM/AnYHCp+5RHH7+T7sdm4Oz0sqHAWuBFYCPwE8ooqgT4KvB/pCyYq9PLrgdmpv8ekr4vm9P36RDTvlen99sEfKXUfQm6j8BJ6Xu2HngBOL7Ufcmzn59Of/f2AB8CG50+u+X4k2sfgWOBDel7uQE4p9R9cfrR6Qc0Go2mH1IubhmNRqPRBIgWd41Go+mHaHHXaDSafogWd41Go+mHaHHXaDSafogWd41Go+mHaHHXaDSafsj/B3GVZ4k5+esUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Find the index of the single variable with the best R^2\n",
    "imax = np.argmax(Rsq)\n",
    "\n",
    "# Regression line over the range of x values\n",
    "xmin = np.min(X[:,imax])\n",
    "xmax = np.max(X[:,imax])\n",
    "ymin = beta0[imax] + beta1[imax]*xmin\n",
    "ymax = beta0[imax] + beta1[imax]*xmax\n",
    "plt.plot([xmin,xmax], [ymin,ymax], 'r-', linewidth=3)\n",
    "\n",
    "# Scatter plot of points\n",
    "plt.scatter(X[:,imax],y)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the calculations above could have been done without a for-loop using Python broadcasting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove the means\n",
    "ym = np.mean(y) \n",
    "y1 = y-ym  # a column vecotor each minus mean\n",
    "Xm = np.mean(X,axis=0) # averaging over column, resulting a row vector of dimension natt\n",
    "X1 = X - Xm[None,:] # minus the same mean in each column\n",
    "\n",
    "# Compute the correlations per features\n",
    "syy = np.mean(y1**2)\n",
    "Sxx = np.mean(X1**2,axis=0) #a row vector with each element indicating the variance of one attribute\n",
    "Sxy = np.mean(X1*y1[:,None],axis=0) #a row vector with each element indicating the covarance on one attribute to the targer\n",
    "\n",
    "# Compute the coefficients and R^2 value per feature\n",
    "beta1 = Sxy/Sxx # element wise division, resulting a row vector containing  beta1 for each attribute\n",
    "beta0 = ym - beta1*Xm # element wise multiplication, resulting a row vector containing beta0 for each attribute\n",
    "Rsq = Sxy**2/Sxx/syy #a row vector containing Rsq for each attribute"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us see whether we get the same result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 0 Rsq=0.035302 beta0=152.133484 beta1=304.183075\n",
      " 1 Rsq=0.001854 beta0=152.133484 beta1=69.715356\n",
      " 2 Rsq=0.343924 beta0=152.133484 beta1=949.435260\n",
      " 3 Rsq=0.194908 beta0=152.133484 beta1=714.741644\n",
      " 4 Rsq=0.044954 beta0=152.133484 beta1=343.254452\n",
      " 5 Rsq=0.030295 beta0=152.133484 beta1=281.784593\n",
      " 6 Rsq=0.155859 beta0=152.133484 beta1=-639.145279\n",
      " 7 Rsq=0.185290 beta0=152.133484 beta1=696.883030\n",
      " 8 Rsq=0.320224 beta0=152.133484 beta1=916.138723\n",
      " 9 Rsq=0.146294 beta0=152.133484 beta1=619.222821\n"
     ]
    }
   ],
   "source": [
    "for k in range(natt):\n",
    "    print(\"{0:2d} Rsq={1:f} beta0={2:f} beta1={3:f}\".format(k,Rsq[k],beta0[k],beta1[k]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Improvements with a Multiple Variable Linear Model\n",
    "\n",
    "One possible way to try to improve the fit is to use multiple variables at the same time.  We can fit the multipe linear model using the `sklearn` package.  \n",
    "\n",
    "For reasons that we will explain in the next demo, we need to split the data into two parts:  one part for training the model and a second part for testing the fit.  In this example, we will use `ns_test=300` samples for training and the remaining `ns_test=442-300=142` for testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "ns_train = 300\n",
    "ns_test = nsamp - ns_train\n",
    "X_tr = X[:ns_train,:]     # Gets the first ns_train rows of X\n",
    "y_tr = y[:ns_train]       # Gets the correspoinding rows of y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To fit the linear model, we first create a regression object and then fit the data with regression object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regr = linear_model.LinearRegression()\n",
    "regr.fit(X_tr,y_tr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see the coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "152.34786451820108"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regr.intercept_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -16.57607993, -254.66532396,  560.98630022,  278.91811152,\n",
       "       -393.41357305,   97.05460405,  -19.0023093 ,  169.46450327,\n",
       "        632.95050374,  114.21638941])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regr.coef_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next compute the RSS (per sample) and the R^2 on the training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RSS per sample = 0.485281\n",
      "R^2 =            0.514719\n"
     ]
    }
   ],
   "source": [
    "y_tr_pred = regr.predict(X_tr)\n",
    "RSS_tr = np.mean((y_tr_pred-y_tr)**2)/(np.std(y_tr)**2)\n",
    "Rsq_tr = 1-RSS_tr\n",
    "print(\"RSS per sample = {0:f}\".format(RSS_tr))\n",
    "print(\"R^2 =            {0:f}\".format(Rsq_tr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see R^2 is higher than the best single variable model -- more than 51% of the target variance is explained by the model instead of just 34% for the best single variable model.  \n",
    "\n",
    "We also create a scatter plot of predicted vs. actual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deXxU1fXAvydhgACWQEXEgGKVoiACQtGKtYBaUKumuLf+3Itrq5SiQVuBqiWKStVarRX3BdyKKCAiEBcstiIgIiLUBQmIWgiyDDBJ7u+P92aYzLz35s1kVnK+n08+mXnreXdm7rn3nHPPEWMMiqIoigJQlGsBFEVRlPxBlYKiKIoSQZWCoiiKEkGVgqIoihJBlYKiKIoSoVmuBWgMe++9t+natWtK527bto3WrVunV6AMUkjyFpKsUFjyFpKsUFjyFpKs0Dh5Fy1a9K0xpoPjTmNMwf7169fPpMr8+fNTPjcXFJK8hSSrMYUlbyHJakxhyVtIshrTOHmB94xLv6rmI0VRFCWCKgVFURQlgioFRVEUJYIqBUVRFCWCKgVFURQlQsaUgoi0FJF/i8hSEVkuIuPt7Y+KyGcissT+62NvFxG5R0RWi8gHInJEpmRTFEVRnMnkOoWdwBBjzFYRCQBvi8gse99oY8zzMcefCHSz/44E7rf/K4qiKFkiYzMFOxx2q/02YP955ek+DXjcPm8hUCoinTIln6IoSkESCkFlJXt9/HFGLi8mg/UURKQYWAQcDNxnjLleRB4Ffow1k5gLVBhjdorIK0ClMeZt+9y5wPXGmPdirjkCGAHQsWPHflOmTElJtq1bt9KmTZvUHiwHFJK8hSQrFJa8hSQrFJa8hSBrm1Wr6D5xInutWsXq009n7dVXp3SdwYMHLzLG9Hfc6baqLZ1/QCkwHzgM6AQI0AJ4DLjJPmYGcEzUOXOBfl7X1RXN+UkhyWpMYclbSLIaU1jy5rWswaAxN9xgTHGxMR07GvP884W9otkYUwNUAcOMMettuXYCjwAD7MPWAl2iTusMrMuGfIqiKHnLggXQpw/8+c9w/vmwYgWcfnrGbpfJ6KMOIlJqvy4Bjgc+DvsJRESAcuBD+5TpwPl2FNJRwGZjzPpMyacoipLXbNkCv/kN/OQnsGMHzJ4NDz8M7dpl9LaZjD7qBDxm+xWKgGeNMa+IyDwR6YBlQloCXG4fPxM4CVgNbAcuyqBsiqIo+cvs2TBiBHz5paUYbr0VsuTvyJhSMMZ8APR12D7E5XgDXJUpeRRFUfKejRvhd7+Dxx6DQw6Bt96CgQOzKoKuaFYURckHXngBevSAJ5+EG2+ExYuzrhCgwIvsKIqiFDzr18PVV8OLL8IRR8Crr1qO5RyhMwVFUZRcYAw88og1O5gxAyor4d13c6oQQGcKiqIo2efzzy1H8pw5VnTRQw/BD3+Ya6kAnSkoiqJkj7o6uOceOOww+Ne/4L77oKoqbxQC6ExBURQlO6xYAZdeCu+8A8OGwd//Dvvvn2up4tCZgqIoSiYJhax1Bn36wMcfwxNPwMyZeakQQGcKiqIomWPRIrjkEli6FM46C+69F/bZJ9dSeaIzBUVRlHQTDEJFBRx5JHz9NfzznzB1at4rBNCZgqIoSnp5803Ld7BqlTVLuOMOKC3NtVS+0ZmCoihKOvjuO7jqKvjpT6G2Fl5/3Qo1LSCFAKoUFEVRGs+sWVaY6f33w7XXwrJlcNxxuZYqJdR8pCiKkir/+x+MHGlFFPXoYYWbHnVUrqVqFDpTUBRFSRZj4Nln4dBD4Zln4Kab4P33C14hgM4UFEVRkmPdOrjySnjpJejf3/IdHH54rqVKGzpTUBRF8YMxMHmyZSaaPRsmTrRSVexBCgF0pqAoipKYTz+FX/8a5s2zooseeggOPjjXUmUEVQqKomSVaYurmTh7JetqguxXWsLood0p71uWa7GcqauzViHfeCMUF8MDD1jKoWjPNbKoUlAUJWtMW1zNmBeXEQzVAVBdE2TMi8sA8k8xLF9uLT579104+WRLIXTunGupMk7G1J2ItBSRf4vIUhFZLiLj7e0Hisi7IrJKRKaKSHN7ewv7/Wp7f9dMyaYoSm6YOHtlRCGECYbqmDh7ZY4kcmDXLvjTn6BvX/jvf+Hpp+Hll5uEQoDMOpp3AkOMMb2BPsAwETkKuA2YZIzpBmwCLrGPvwTYZIw5GJhkH6coyh7EuppgUtuzzn/+Y0UUjR0LZ5wBH30E554LIrmWLGtkTCkYi63224D9Z4AhwPP29seAcvv1afZ77P3HiTShT0JRmgD7lZYktT1rbN/ODx54wFpnsHEjTJ9uzRA6dMitXDlAjDGZu7hIMbAIOBi4D5gILLRnA4hIF2CWMeYwEfkQGGaMWWvv+y9wpDHm25hrjgBGAHTs2LHflClTUpJt69attGnTJrUHywGFJG8hyQqFJW8hyQrx8tYEQ1RvClIf1e8UiVDWroTSkkAuRKR0yRJ+eMcdtKquZt3Pf85/L7uMugJo48Z8FwYPHrzIGNPfcacxJuN/QCkwH/gJsDpqexdgmf16OdA5at9/ge97Xbdfv34mVebPn5/yubmgkOQtJFmNKSx5C0lWY5zl/ef7a83RE+aarte/Yo6eMNf88/212RfMGGNqaoy57DJjwJiDDjKL77orN3KkSGO+C8B7xqVfzUr0kTGmRkSqgKOAUhFpZoypBToD6+zD1tpKYq2INAPaAhuzIZ+iKNmjvG9Z7iONXnkFLr8c1q+H3/8exo+n5t//zq1MeUImo486iEip/boEOB5YgTVjOMM+7ALgJfv1dPs99v55tkZTFEVJD998A7/8JZxyCrRrZ61InjgRWrXKtWR5QyZnCp2Ax2y/QhHwrDHmFRH5CJgiIrcAi4HJ9vGTgSdEZDXWDOGcDMqmKEpTwhiYMgV++1vYvBnGj7cqozVvnmvJ8o6MKQVjzAdAX4ftnwIDHLbvAM7MlDyKojRR1q6FK66wTEYDBlj5iw47LNdS5S177lptRVGaNvX18OCD0LMnzJ0Ld91l1TtQheCJprlQFGXPY/VqK0dRVRUMGWIph4MOyrVUBYHOFBRF2XOorYU77oBevayiN//4h1XvQBWCb3SmoCjKnsGyZVYCu//8B049Ff72NyjLsyR7BYAqBUVRUmba4mrGv7ycTdtDAJSWBBh3as/srkPYuRP+/Gfrr107K8rorLOaVL6idKJKQVGUlJi2uJrRzy8lVLd7OVFNMMTo55YCWUqF/e671uxg+XI47zyYNAn23jvz980i2a4/oT4FRVFSYuLslQ0UQphQvcl8Kuxt2+B3v4Mf/9had/DKK/DEE3ukQhjz4jKqa4IYdtefmLa4OmP3VKWgKE2IaYurGVg5jwMrZjCwcl6jOhevdNcZTYU9b55VF3nSJCtVxfLlVhGcPZBc1J9Q85GiNBHSXfVsv9ISql06/4ykwq6pgdGjrfrI3brBG2/Ascem/z4O5KqEqHf9idYZuacqBUVpIniNOlPp4EYP7R7nUwAIFAmjh3ZvlKxxvPSStSp5wwa47joYNw5K0q94nDp/IGclRN0UbybrT6j5SFGaCOmuelbet4yJZ/SmXavddRBKSwJMPLN3+jrLr7+Gc86B8nKr4M2778Jtt2VEIdQEQ472+/EvL89ZCdHRQ7tTEihusK0kUJx+pRuFzhQUpYmQiVFnOtJgO5pm+uwHTz0F11wDW7fCzTfD9ddDIHOFeDZs3kEw1HCcHAzVxSmEMNkoIRpuWyfTVVXVqozcU5WCojQRRg/t3sAMApkfdSbCyc9x9yPzOGrxo+z79jyrPObkydCjR8Zl2VVXTzLGk2yVEM12/QlVCorSRPAadeaKaD+HmHp+teRVKqoeodgYuPtuuOoqKC5OcJX00LzYWSGIWJm3o8m1Ms0kqhQUpQmRF1XPogibYA7cWE3lrHs4cu1y3jqgDzcMu5q3fntJVmXp2LYlJYF4c1GsQsj0qu1cRTqFUaWgKErO6LJXc058/WlGvv00O4sDjD7xGp7rdTxl7bJfCa20JMCE4T0iHXKRCHUOxR9bt2iW1k46WgmUtgqwdUctoXrrvtmMdAqjSkFRlNywdCnTnxxF6cfLePWHP+aPJ1zBN23a59Q0Ez2TOrBihuMx6XQwx/pUwjmkomlM2HAqqFJQFCW77NwJt9wClZWUtm/Pv29/gJtrD+bbzTsoS8FckilzSzbWCDitHXEiG5FOYVQpKIqSPd55By69FFasgAsugDvvZMD3v8+CFC+X7lXa0WQjWstvZ5+tSCfQxWuKomSDrVutNQfHHGMls3v1VXj0Ufj+9xt12UzmBirvW8aE4b0oKy1BgLLSEiYM75VWM46fzj7b5rSMzRREpAvwOLAvUA88aIy5W0TGAb8GvrEPvcEYM9M+ZwxwCVAH/NYYMztT8imKkiXmzIERI+Dzz+Hqq626B3vtlZZLp3uVdiyZjtZymo0EioXWzZuxORja46KPaoFRxpj3RWQvYJGIzLH3TTLG3BF9sIj0AM4BegL7Aa+LyA+NMYkNboqixJHr0EY2baL7bbdZs4Lu3eGtt6yZQhrJRW6gdJKPa0cyphSMMeuB9fbrLSKyAvB60tOAKcaYncBnIrIaGAD8K1MyKsqeSiZt7b745z/hyivZ9+uvYcwYuOkmaNkypUuFlVt1TZBiO0w07JDOx1XayZJva0fEOMThpv0mIl2BN4HDgN8BFwLfAe9hzSY2ichfgYXGmCftcyYDs4wxz8dcawQwAqBjx479pkyZkpJMW7dupU2bNimdmwsKSd5CkhUKS16/sq78aoudtqEhzYuL6L5vekw3TjTfuJGD77mHfd54gy0HH8ziq66ivk+flK9XEwxRvSlIvUM/VSRCWTtrRrBh8w521dXTvLiIjm1bUlqSfI6kQvoeQOPkHTx48CJjTH+nfRmPPhKRNsALwLXGmO9E5H7gZsDY/+8ELgacCqrGfROMMQ8CDwL079/fDBo0KCW5qqqqSPXcXFBI8haSrFBY8vqV9aKKGRiHOBIBPqtMfH7SGAOPPw4jR8L27fDnP7PX739P/YIFvuR1M3UNrJxHdY17mouy0mIWVAxJyyMU0vcAMidvRpWCiASwFMJTxpgXAYwxG6L2/wN4xX67FugSdXpnYF0m5VOUPZWs2tq/+AIuuwxmz4aBA60iOIccAlgj/YGV8+I6+2gl0LYkwLZdtZG6DNGmrkQO42zG7zcVMhaSKiICTAZWGGPuitreKeqwXwAf2q+nA+eISAsRORDoBvw7U/Ipyp5MVvLw19fDX/8KPXvC22/DvffCm29GFMK0xdVUbwrG1Sf4w7RlDeoW1ARDcYV6wmGliZTYfqUlaS0xqmR2pjAQ+D9gmYgssbfdAJwrIn2wTEOfA5cBGGOWi8izwEdYkUtXaeSRkg3cRrOFTMajWlauhEsugQULYOhQZv9mPH9ato11N8yK3Gvi7JWc0yW+s3/m3S8dcwrFsq4myKSz+8Q5ksOUBIoZfEiHnDrUcx7hlQEyGX30Ns5+gpke59wK3JopmRQllt2jWWtUnYsEZJkiI1EtoRDccQeMHw+tWsGjjzKt13GM+eeHcR1zMFTX0CBs40chgDULiFZuTtFH6S4xmgw5j/DKELqiWUkb4Wn8surNGZvGp9tUMHH2yrjIlmyVWswWaWuzxYthwAC44QY45RT46CO44AImvvaJY8dcLE5jQueRYizhWcDAynmMnGoZGv5ydh/uPKs3ZaUlrKsJRhSFE9nwNWRyNXUu0dxHSlpoMGrqkplRUyZGZutqgo6j2Vw7MNNllki1zaLv37V1MX//7BV++PgDVp3kF16A4cMjx7q1VZ0xFMUohkCRUA/U1TdUxEVA21YBaraHKG0VYGeojicXronsr64JMvq5pSA0cEi7kY3Fa5leTZ0rdKagpIVsjJoycQ+3ziOXK2LDHXmsg9ZthO81E0ilzaLv32/tch76y6X88JH7+OLkM6zZQZRCAPe2KistoaxdSYPcQW1aNotTCGAphMU3/YxJZ/dhR6ie7aH4NRahehPnkHYiUCSeDvV0zZzy8buTDnSmoKSFbIyaMnGP0UO7U71iUYNtuV4R69WR33pUw3FcoplAojZzmpFMnL2Soq1bGP/mY1zw/gy+bNuR8866mc/6Hs2Cdu3iruW1qrh08yoWVAyKbHerUVBj1xHwm0raizYt3YvgeLVXaZL32RNWUzuhMwUlLWRj1JSJe5T3LYsbzaY7E2ayJKP8Es0EvNrMbUZy0PtvM/vhq/i/92fySL9TGHrxX3n7wL6uciWTTTTRZ5iOQURYwTiRztlmNrKo5gKdKShpIRujpkzdo7Qk0GA0m2uSWXiWSIF4tVlsB9k2uIWbZvyD0z+cx+r2nTnjV7fzfudDPe8fxm+kU6LP0O3ZwTILRfsU3PCS07u9WieQPp58y1uUDnSmoKSF6FETZGbUlOmRWb4sgkpm4VmikbdXm0U6SGM48eO3ef2hKzj1oze498dnc/qI+xoohHQp+ESfodOzg6W4J57Zm4ln9I6c265VwFIUUSSSc0/1A6QTnSkoaSM8aqqqquI3vxqU0Xukm3yKOfdaeFZVtarBsX5mT25ttl9pCbvWVnPznPsZ9sm/WNbxIC44609s7t6T8fZMIhOLsrw+Qz+L7qJfJxul5dlem1e5npfORWr5vuBNlYKi4G1rzsUP1q/yS3nlsjHcu2MJB00eS4vaXUwYdCEP/egXNG/RnAn2+amGwDaQpXfyTmO/906lcw3vHzd9OTVBy/fQMuBtMEnngCGfBh9uqFJQFJJz7ibbGaVzZOiWkiOp6332GYwYwRGvv863fY/kssFX8G5g77SMgGM7vOpNdUxbXB13zca2SWM71521u0NeN20PMebFZUw42jkbazoHDPk2+HBClYKi4N+5m2xnlO5RZqNSctTVwX33WUVviorgb39j78suY0pRelyLTh1evTFxHV6iNvGjMBrTubqdu2Gzc9RSOkOhC2HBmzqaFQX/zt1kQxrTGQLZqJQcH30EP/kJXHMN/PSnsHw5XHGFpRzShN8Oz6tN/C7ca0zn6naMU1EiSK9zuhAc3aoUFAX/kU3JdkY5H2WGQnDLLdC3L3zyCTzxBMyYAfvvn/T9E+G3w/N6Dr9KtDGdq9sxzYudu8N0piHPSkrzRqLmowJhT0zvnG/4sc0nW7wmncVurHO2+L/WokVw8cXwwQdw9tlwzz2wzz5J39cvTpE9RRKfcsKrTfyuwK6uCSI0LM3ot3N1i0Dq2La54/HpTEOe8ZTmaUCVQgGwJ6d3LjSSXUCXzgV3vlNyBIMwbpyV4rpjR5g2DU47Len7JYtTh1fWri7uO+rWJoMP6eBaayF6BXb4PAMRxVDm0Lm6+SbcOuZSj5DUdIZC5/uCNzUfFQBNIb1zLkgl1XeyC+jSueDOV0qON9+E3r3h9tutWcJHH2VFIUTLuKBiCJ9VnsyCiiGUlgQcj4ltk9P7lfHCompHheC2Aht2K4QFFUMcndluvonyvmWMHto9MjuZOHtlJES1qeM5UxCR33ntjy6zqWSOfE3vXMhMW1zN6OeWEqo3kVTf105dwntfbOSW8l6e5yY70kvnyNA1Jcd330FFBdx/Pxx4ILz+Ohx3XFrumQli22Rg5TzHRHjFIhHFF66rEEuyOaHCEU5+w2ebGolmCnvZf/2BK4Ay++9yoEdmRVPCFELEQqExbvpySyHE8OTCNYVX43fmTKtO8gMPwMiRsGxZXisEJ9zyHdUZE+mkk/kdJPJNeIXPNnU8lYIxZrwxZjywN3CEMWaUMWYU0A/onA0BFcsGG1usJN8iFgoNL1NBwXQM334L550HJ58M3/sevPMO3HUXtE4+sVuucavSFr09nTmhCmG9QK7w61PYH9gV9X4X0NXrBBHpIiLzRWSFiCwXkWvs7e1FZI6IrLL/t7O3i4jcIyKrReQDETkihefZI8nH9M57MnnfMRgDU6dCjx7W/7Fj4f334aij0nqbVBIExp7j107vVrc5ensy/plECkRn3+74jT56Avi3iPwTy7fzC+DxBOfUAqOMMe+LyF7AIhGZA1wIzDXGVIpIBVABXA+cCHSz/44E7rf/K+RfeudCp12rAJtc8u6n2jEkk7oh1TQPzb/9FsrLYfp06N8f5s6FXt4+kFRIZdUxkLKdvswlTLUs5rNIV04ov+GzTRFfSsEYc6uIzAJ+Ym+6yBizOME564H19ustIrICyx9xGjDIPuwxoApLKZwGPG6MMcBCESkVkU72dRSlAY3NnTP2lJ6Mem5pXGnIRKUcveTxm84ipdQXxsDkyQy49lorXcUdd1irk5tlJqo80SIyJ/lbBop8pblwIhO1MpLNxuoUPtsUEeMybYs7UOQYoJsx5hER6QC0McZ85vPcrsCbwGHAGmNMadS+TcaYdiLyClBpjHnb3j4XuN4Y817MtUYAIwA6duzYb8qUKb7kj2Xr1q20adMmpXNzQSHJm2lZa4IhqjcFG4TpFolQ1s4aVW7YvINddfU0Ly6iY9uWjmGR4eusrwny/RaGDUEoLhL2Ky1xPd6LlV9tcUyT0Ly4iO777pXysQAtq6vpfuedtFu8mP8ddhirKyoIlmW281pWvdl1X/PiIteUELF0LIENQehV1jbhsTXBkO/PLln8XLuQfmPQOHkHDx68yBjT32mfr2GGiIzFikDqDjwCBIAngYE+zm0DvABca4z5TlwcSljrUGKJ01jGmAeBBwH69+9vBg0a5OMJ4qmqqiLVc3NBqvLmInd7ptt2YOW8yEK+aEpLDDtr6wmGigi7y0oCdUwY3sPzmauqqjjLRV6/7XdRxQyMg4tOgM8qB6V2bF0d3H03/OEP1ozg739n2cEHM2jIENdn8YOfZ7qxcp6rOWddTdBRfidG9aplypd7Zay+hh+mLa5m9GtLCdXv/l4EikJMPLPh96Kp9AmJ8Oto/gVwKrANwBizDitU1RMRCWAphKeMMS/amzeISCd7fyfga3v7WhpG43cG1vmUb48jVYdd7DX8JBcrNNwcwTXBUNqSz4Fz+42cuoSuDo7XZByXvo798EM4+mgYNcoKL/3oIxgxotEJ7BJ9J8Lfu3AaiWjC5hw3+UtLAnHO3Xyw0zuFH4fqDeOmL8+RRPmN32/YLtvWbwBEJGHMm1hTgsnAiphFbtOBC+zXFwAvRW0/345COgrY3FT9CU4/3OpNwaQ783Rm6MwnknUEpxpN5LaCFnZ3pn+Ytsy1Ew0UCdt31cZF73hGxuzaBePHwxFHwKefwtNPW07lzlYEeDgHVqolQ/1mKA0/a/iZoiN93OQfd2rPuOigsnYlWbPTu0VLuQ2odAWzM369VM+KyN+BUhH5NXAx8FCCcwYC/wcsE5HwUsQbgEr7epcAa4Az7X0zgZOA1cB24CLfT7GH4TcvfSL21FhsN6dky0CRY0RRqtFEidopGKrjqYVrIooiOhdPaUmAbbtqI/I4OZPjTDi166DfMGuW8Mtfwl/+Ah06RO6XjhxYyWYojU4jESZRZE+sScaLdJk3G6xQx2qb0c8tTfo6iv/ooztE5ATgOyy/wk3GmDkJznkbZz8BQNxyS3smcpUfefZ00tWZl7qEXfrpJFP9sU5bXM2Gr7ZwUcWMjPkw3DolIE5ZBIqFbTut0Xqy8rhl84wm1ukV7kQhfiQanWahQWTM9u1w000waRJ06mTNDE45Je5eE2ev5JwuzjmwGvtMfjKURpOO1B3pLEDkZSJyCz9u1yo9Tuw9DV/mIxG5zRgzxxgz2hjze2PMHBG5LdPCNVXSsbBm2uJqtu6ojdseKE5s43WzO4dNJW6mi/B5u+rqM+rD8Mp+GW2+aNcqAMbqnFORx8lM4od1NUH/Hez8+XD44XDnnfDrX1vFbxwUguO5CbY74WW6yvaCrnSaN71MRGNP6UmguOH4NFAsjD2lZ9L3aQr49Smc4LDtxHQKouzG6YebrMNu4uyVjrl9WjdvlnK5wqcWrvF0WmfDh+En+2U4S2er5s3i2iAYqmPUs0uTyoia7Ihyv9KSxB3s5s1w2WUQjiSaN8/KXdTWOXRz2uLquFQncdf0gdeq4GwUgIm2+7vNwtJt3izvW8bEM3o3eOaJZ/TWNQkuJMqSegVwJXCQiHwQtWsv4J1MCtaUScfCGrcf1mYfzjW3c2NVTDBUx7jpyyNyZcOHkUxtXrf71hnj20xR3reMibNXOpofSgJFgLguuHJdjPXyy3D55fDVV/D731uO5VatXGUIK0KvtNLJ4Gb6yXQBmFhzkRupzEwSmYjyvYZBPpHIp/A0MAuYgJWOIswWY8zGjEmlxH2JEznsYmlMxS8/tvQwNcFQJI1BOquMxRJdccsJJwXg9RzJ2OLdlMuOUD2Tzu7j2YlG77txwN6cNHE0PPOMlZpi2jT40Y8S3t9JEULDtNLpIpOdp9tzRJPqzGTsKT0Z/fxSQnW7FaeaiFLDUykYYzYDm0XkbmCjMWYLgIjsJSJHGmPezYaQSvI0Jm2A07mxpQ+jCXeu4fOstFfJ3dMLPyNMJ8Xj9BzR+J3BeCm7RKkUyvuWWSkqnnkGzjrLqnswfrxV+6C5c/lHv3LWR6WVLgS82lvA98zEKwgin8tcFgp+Q1LvB6Kzlm5z2KbkEY35kTidO/iQDjy5cI3j8eEfe/i8DSvfb/AjBxpVXzrRCNNL8bRoFp+PJ4zfGUwyCja2w/pjn70Ydt+f4JVX4MgjYfJkq/ZBEmRyBpZN3J4jNuTVi0QRS6oEGo9fpSAmKkmSMaZeRLS+cx4T3TmVtgqwbWctI6cuYeLslb46Zacf2IwP1icMcS3vW0bV5lWRdA3pCDv0GmE61eZ1uq8Tgw/p4LovGr8KNvqeYur5adWLHDP2YWqLDM3uugt++1soTj6aKRPJ4nJBOp4jGZ+Skhp+O/ZPReS3WLMDsJzPn2ZGJKWxxHaI0R15Y2LBx57SM+kfdWN/xOGoGycnq9cI04/9ev7H3yS8fxg/o9DwPbturKby1Xs56ssPWXDA4fzlrNE8N/I83/dyunf4+rDFVRHmO+kw8eypCzLzCb9K4XLgHuAPWKbludiZSpX8I1GHmOrIKpUfdao/4mmLqxk3fblnKh90q7EAACAASURBVIJtO2tdc/X76STS3ZFs2LiVX//nJUa9/SS7igNcN+y3PHv4CXgkgfRNWClVVVU1KrlcLhIkRtNYE8+eYkrLZ/yuaP4aOCfDsihporEdolfHkeyPOpUfsd/QxZpgyHXW4yeCKq0dyQcf8PLT13Fo9Upe63YUfzjhCr7e6/vpv08jSNWUl2tFEs2eYkrLZzwXr4nIdfb/e+1SmQ3+siOikix+w06dSHdm1VQWRPkx/YRxWxw3emh3AkXeI/TRQ7szbXE1K7/aknKCOXbutMph9uvHD7Z/y8jhYxjxixsjCiH2WVMpcelFMtdLZXFhvmXaTaYkp5IaiWYKK+z/73kepeQViUIxwx2V0wgw3Y68dJqckj7eQyeEFzWNeXEZVx5Sj6EoeX/LwoVwySXw0Ue82ud4xhx7Mab992knULM9FPes6cz1k8r13NqpuibomhsqHx27GmWUWRKtU3jZ/v9YdsRR0kFsR1zaKoAx1mpmr3q66Yjpd5MnHSYnN0qj0lAkWuQGllIce0rP1Du8bduswjd33832ffZl5DnjmX1Av/AFKAkUM+nsPpT3LYuM5NfVBB0d5sFQHeNfXu54PyelXRq1P1n5vdo1ehYAu79DjXXshlN9p6tudT6ZsvZUEqW5eBn3NUsYY05Nu0RKWkjUEQ+snOfYoRS7RPqkyy7u50c9emh3rp26xOUK8WzdURsxZ/jxRYTNDW738FRIc+daies++wyuvJLTOgxj1Y6GVli3WsZO7QpWdFisw9xtFjDh6N2muGQ77EQzyGjZw7I0xrGbTKpvP7OedM+0FGcSJcS7A7gT+AwIAv+w/7YCH2ZWNCWTeOUF8uMDSMU27tc+Xd63LKkkdKF6q9aEH19EWenuoi/FLlFBjttrauDSS+H4463SmG+8Affdx+odzj8ht/oEbsTa9d1mARs274i8Tzaraaw93o3o70ZjkuRNnL2yQR1tcPdh+PF37KlFo/INT6VgjHnDGPMG0NcYc7Yx5mX775fAMdkRUckEbh1H2HHn5chL1fmYzI967Ck9k0pb7ZWuOozYsoaVmNvIvc4Y+ox/bffzvPQS9OgBjz4K118PS5fCsccC3h1zMia32GPdzt1VVx95nUqHHZ1FtsyHUmmMYzeZmYyfY3WNQnbwu06hg4j8wBjzKYCIHAj4Ww6q5CVeoX2xPolwpx293alzj7WNx5qKkklmFyuDqw3TJtyRud0jOndTWIm5ZdYEyxZe+WgV/a+fSuc5L0Pv3lZ20379Ghzn1Y6JfBtO8odpWxJwXaMxsHJeWvL9+A3vTNWxaz3TFsftsd8NPwWhdI1CdvCrFEYCVSISXsXcFbgsIxIpWcGrQ3Gy3Y6cuoT3vtjILeW9XEdmYdt4Kc72X7ekel7mjrCc4TrIToRrIW/aHoq7R0mgmBbNihyroLVoVmSPtGOKERnDL5bP56a5/6B1aAfccgtcdx0E4k1aiTpmp0pwdXWG+qhrBIria2V4rXerrgky+vmlkftH29wnzl7JyKlLfCmIZJWKlz/Iad/ood2pXrGowTVKAsUMPqRD3HcjUCQEiqVBltNYBZXpNQrJOMX3ZPwuXntVRLoBh9ibPjbG7MycWEqm8OPodavV+9TCNfQ/oL3nqH/Us0u59rAQf5m9NM48E13DOEz0j9pLNjcnaUmgiNp6ExllRt8jnA5ipItDeXMwxKSz+/DlR7sjrvf77mtunX0fgz9dxKL9DuG6E69h7o2XO54fJpn6BIMP6cDU/3xJfVTn52Tgr3GZwYQJ1ZkGM7NUnbB+ZwFe1wfnSLYJw3tR1q6EstLihGHPoXpDaUmA1i2auX43M5kJNR31r/cUfCkFEWkF/A44wBjzaxHpJiLdjTGveJzzMPBz4GtjzGH2tnHAr4Fw0pkbjDEz7X1jgEuAOuC3xpjZKT6T4oLfjsOryE44oZ5b5E5YEbjZ68OdtZ/ZiVP2S6d1FbEKKlohTJy90tX0FE59/feV7yOmnl8tnkXFG49SZOoZd9wIHj/iZExRsWsqDT/EdroDK+c1GA2D1cHHhpH6CcuNNrdkej1BIn+Q277RveOv5VUAasnYn3nKkak1Cumof72n4Nd89AiwCPix/X4t8BzgqhSAR4G/Ao/HbJ9kjLkjeoOI9MBKo9ET2A94XUR+aIzxF7qh+MJvx+Fm3wXrB13etyxhXiI33JLY+ZHNqUNwmwUkWncRPUP5wZZvefaZO/nRl8t5s2tfbhh2NWvbdmwgW6KOwW/8vF9nqZ/w0VSumyqpXL+6Jkj1prq40bcf/0G2WVcThC4u25sYfms0H2SMuR0IARhjgniuFwVjzJuA3+pspwFTjDE7jTGfAauBAT7PVXzi94ftMsi39gFdK2YgQsI0ErF42X9T7dTcOpJiEdcONRJB06sj3HYbx111Gb1r1vL7k67l/LP+1EAh+JHBbzRWMnWWy/uWcXo/b0Uk9jWdzg9TJJKWlBpeUVZen4FTSKoxZLwWdLIkG967JyPGqwcIHyTyDnAcsMAYc4SIHAQ8Y4zx7LhFpCvwSoz56ELgO6zUGaOMMZtE5K/AQmPMk/Zxk4FZxpjnHa45AjtDa8eOHftNmTLF35PGsHXrVtq0aZPSuU7UBENs2LyDXXX1NC8uomPblpSWJFfw3Yt0yLvyqy0NQhrDNC8uovu+e0XeL6ve7Ot6IoJA3A+/YwlsiOlHi4vEijJxaRO/ssVSEwxRvSnYQIYih84omi7tW9G5+gsOvu02SlevZu2RRzH/gsvZWtqO2vr48xLJ4Ed2JzmjaVYkdIppH6frxrZtkQhl7ezIK4/rRx+byvfSrZ3d7h3+DJy+C2B9Bpn8vSRLTTBEKLidrxzaNpdyedGYPmHw4MGLjDH9nfb5NR+NBV4FuojIU8BArM49We4HbsYacN6MtTDuYpxnHY7fbmPMg8CDAP379zeDBg1KQQyr5nGq58YybXE1Y+YuIxgqIjz5KgnUMWF4j7TZI9Mhb41D9tGSQDEThvdikG3Xt2z0/usnOa2AHtWrljuXNbyGdZ8eDHJpj0SyeeGWw8nJJt+idhc3LHqUU9+ayv9afo/Ly8fQ84wjuXNZMwJfCghxETCJZLioYgbGYdItECk2ZEVPea+7iP3OOF3XqW3LSotZUDGkQTu41aAolhD1ZldKTtpko48sO/0WB3lLWPArf5XWssm0WXOY8klxwUQfpbMPiybhr1+sZPAfA8OBo7C+69cYY75N9mbGmA1R1/0Hu30Sa2lo0esMrEv2+rki10nD/NqzkwlD9YubQzmWRO3RmMgSN+djrDP8iLUruH3W3Ry8cS0zj/gZY465iM0le9HTDkn1EwHjhJ/4eT+26WCojmujquP5zQEVPia6HQ6smOF4bPjzSiW6xsvJ67bPKSQ1X9Ncl5YEWFAxKNdi5JyESsEYY0RkmjGmH+D8TfOJiHQyxqy33/6C3akypgNPi8hdWI7mbsC/G3OvbJLLlZbJhiK6/XiTSckQjVuuJCec2mPa4mrGv7w84ngsLQlEksk1hvK+ZZHrttoVZPSbj3PBoldY9729GXlRJdP2OcxxKuonAiYWP/HzXs77WMKf4en9ynhhUXXCz8UpLYcfhZLpgUt53zKmffVRXEhqPo++Ff+O5oUi8qNkLiwizwD/ArqLyFoRuQS4XUSWicgHwGCsRXEYY5YDzwIfYZmpriqkyKNcOqnSlQ/GS4F5pZto0Ux8O5xj22Pa4mpGP7+0QWdZEwwx+rmlacnXP/aUnhy3ZgmvTb6KCxa9wuNHnMxplz3AT39zXlo/Mz+pIHzqzQjBUB3zP/7GV64iJ6XslALDiUwPXKzRt5VWY0HFEFUIBYBf4/Fg4HIR+RzYhr0+yBhzuNsJxphzHTZP9jj+VuBWn/LkFZlcaRk2DZ3TZQs3xqQ3gPTNUtxGluF4/+jRfDTbQ/UEioXSkkAkNff3W++kJCAJ22Pi7JVxMfuwO8FdozqQTZsov+cPlD/zKF906MLZv6pk3WH9+aPHauPGfGaJ4uc3pxC+Gw7/jV7V7ZQ2wimHUaw5zs3H4KSoNTV108avUjgxo1IUOJlaadnANNTF2TSUrnwwiXIhTZy90tX8EaoztG7RLGJ2qaqqYsLwbgnbw0txNWoE++KLcNVV8M03MGYMB9x0E8+1bNngkOjPDLZElF+mOsBka0SEz4nGLW2EmyKLTYGRSAlqamoFEtdTaAlcDhwMLAMmG2Nqvc5pqmRipaUfB/bgQzrw1MI1rqkj/JJIsSXqpGP3+2kPr44yJdPbV1/B1VfDCy9Anz4wcyb07et6eFjGqqoqfvOrQcnfD/8jazelO2F4L8DfrKUxNno/A5dcB0wo+UGimcJjWAvW3sKaLfQArsm0UIpFItPQtMXVvLCouoFCEOD0ft4dsltHFv6LTqzmNxImlU589NDujH5+aZwJySlBnCfGwGOPwe9+B9u3w4QJMGqUYwK7dJLMyNpvp5yos29MhEwiRa2pqRVIrBR6GGN6QWRBWcFEBO0JJDINuSWum//xN3HnhEnUkbnt94qEcSvCk6iTC7+PjT4ad2pP32kl5IvPuXPu/Ry56j045hh46CHonp2Qx2RH1qmEdDaGZP0DmppagcRKIWJENsbUilc+XyXtJHJgpzKyS9SRue0PR8KEF4WFQ1GdbPE1wZC9mM/fCDrZznDa4mpueGEpZ747neveeAwjws3DrqTXLRWUd3dIYONxHS8nfiL8zORy5bRNxT+Q6dTUSmGQSCn0FpHv7NcClNjvw9FH38uodE0EL3MOuDtDUxnZJerIvPb77cA3bN5hr+7eTTpt01OfmMNjz07kR9Uf8caBR3DD0KupbrsPZXNWUd7Pn1Lw48SPPtbp8/Fq/1w7bVPxD2QqYEIpLDyVgjHGfz1EJSX8pIx2c4amMrJLpEjSYUKw8vXEL4FptG06FIKJE3n0nrEEAy353ckjebHnkEhFmuqaIAdWzPDVmfntNL0+Hy8nfyqdcjpnFqn6BzKVmlopHPwnuVEyQmMiPsr7lvHeFxt55t0vqTOGYpE4J3NsRzP4kA5xvoFoRTJ6aHdGP7eUUFRiuESO39h7nHeAv0ygSfH++3DJJbBkCQsOO5brBv2ab1u3izssOksp4NoWbmvJYjtNr9KjO0L1rk5+t3oTbs76dM8s1D+gpIrfFc1KhmhMxEc4+ii6sM1TC9fwh2nLIvtjUzq/sKia0/uVRVbJtmsVoEWzIkZOXcLAynm898VGYvN9xuf/bChD7D3q6q2yk9GkbJsOBmHMGBgwgB1r1zHmvPFcfPJ1/M9BITQ4LWZVd6ycbsR2ml6lR72c/E6pJ7y2p2tlehinFc3qH1D8oDOFHONnRDdtcTUbvtrCRTGmEbfooycXruHJhWsc8xKFncbhrJqxo9MnF66Jk6Wufnfpx9hZwbadtQ4yGFo3b+aaWM63meTtt63ZwSef8MVpZ3PWwaezoVmryHO61XwOE92h+8nt5NRpti0JJFVMKHxPt3xQbtvTHQ6q/gElVVQp5JhEfoFwx33lIfUYiqiuCXLt1CW+Kp8l6oCSSYK3aXuIPuNfY9uu2si6Aq91C26J5XyZSbZssWYH990HXbvCnDn88r0iNsQWA2J3iofGZikNF+UJj8zDCnDbrvi1moEioXWLZo7tH3Yyuyksp5QU4fPSbe5R/4CSCmo+yjGJkqm5ddyplMIMY4A/TFuW9Ci0JhhyzFXkhFtnltBMMmsW9OwJf/sbXHMNLFsGxx/vOZL2Yypxk6ddq0CDvEDRVdPccjO1admMcaf2dL2nW11oAVfzTaJnmLa4moGV81hWvbnRVdQUxQudKeQBXiO6TK0mfXLhGlo3L2bbrvQnoy0Sd8e02/NsX78Bzj8fnngCDj0UFiyAH/84st9rJO3HVOI2IzMmvnJcWEm5yVqzPeR5T7e60V7q1Ot6yYTPKkpjUaWQ56SSSM0v23fVURIoTqmOghth57VbZxX3PMZw0soF3DL3AdixFf7wB+uvRYsG5yUysyUylbh1um4d+LqaoGsNhNJWAc97en1m6ap1oTmJlEyhSiHPCXeG4J6HMFAsvs060RisEMrYWPtAsdC6ubPN3M81a4IhBlbOaxAGO//jbyIdbRFWRFOHrRu55bW/MXTVQjYdejg88zj07u143XQ4Tp063fDCwFjCTnTHZ0zQ1E4KLEwqnbnmJFKyiSqFPCfceaz96D3H/e1aBRh7Sk9GPbvUdwW0MMUizP/4mzizRjgV9rhTe6ZWorPeREbKsRFNm7aHwBjOXDaHP86bTPO6ELcNuZhDKm/itN4HeF43E45Tr3TUbrOIRLURwjK6rVVIV60LXXOgZAJ1NBcA5X3L6Ny+VVzsf6BYGHuKlTwu1i7uh3OP7JIwrcWE4b1cY+t9FlxrQOear3hi6h+ZOOseVuxzICdedC/3/2g4t8/9r+9rhJ2uB1bMSOh0TXRsed8yytqVODr6G1OdrbxvmWekUTLomgMlm+hMoUAoLQkw8YweruYTt9FkuBB99L5iEc49sgu3lPdi/sffeI5Cw9d3SnGNgeIioa4+sUIqqq/jgvdfYfSbj1MnRdz4syt5us8wjFjjEr+j52RW/vo91i0dtZ80Il5rLtKVYC7bBYGUpo0qhQLCy3ziZscWwbMD8dNxlfctc1wXUQ+ID4Vw8LdruH3W3RyxbiXzf9CPG4ZezfrvdWhwjN/RczJO18Y6aBP5MfzkrfI6PxnSURBIUfyQMaUgIg8DPwe+NsYcZm9rD0wFugKfA2cZYzaJlZP7buAkYDtwoTHm/UzJli3SnTrZ63rh/7Gd96btoYQRL0CDlNixi7jA3Y7upRICdSEuX/g8V/9rKtuat+Kan4/ipR6DIgnswoSVkJ/2Ssbpmg4HrZci9qN0dAGZUmhkcqbwKPBX4PGobRXAXGNMpYhU2O+vx6rq1s3+OxK43/5fsKQ7wZnfGgVbdsRHzER38l4pur3kdTNPOaXSAOi1fhW3z7qbQ7/5nOmHHsv440bwv9aljs/mVJIy+v7RcicqQB+tWPwWq08VjQpS9kQyphSMMW+KSNeYzacBg+zXjwFVWErhNOBxY4wBFopIqYh0Msasz5R8mSYdseXRHdyoXnUEQw2djdHXCyshtwikcCfr1um7yTvq2aWAs5lJsFJpRKd0aBHaydHPPMGVM6fzTetSLh3+R17v5q7fy+zFZwMr5znef9z05eysrY/sc3q+6JlGtIxex6YDjQpS9kTEpBC14vvillJ4Jcp8VGOMKY3av8kY005EXgEqjTFv29vnAtcbY+LiMEVkBDACoGPHjv2mTJmSkmxbt26lTZs2KZ3rh2XVm1339Spr63luTTDE+pogtVH2+o4lsMFlANqrrC0rv9pi1zFwRhCMg7GneXER3ffdy1PeIhHK2lkd3YbNO1zvs9+KDznuofso3bCeDwefwIJzLmBXa+827tK+FaUlAc/7e9GsSOhUWkJpScC1DcLP3ry4iI5tW1Ja0rB2s9N3oSYYijyr23k1wRDVm4INIr/CbRV7bLrI9Pc23RSSvIUkKzRO3sGDBy8yxvR32pcvjman4EZHbWWMeRB4EKB///5m0KBBKd2wqqqKVM/1w42V8xxHkWWlJZ6OwmmLq20zUcNZwahetdy5LP7jCl/voooZGJcIY69VywJ8VjnIVd7d9ylmQcUQAAbGHNtm53Yqqh7h9CWzWFO6L/PH/ImR9UfAp9b1jz6oPe+v2Rw3y/jVUfvzmxMt01Gi+7tREihmwvAeDOpb5toG1jOe7HqN2O/C7s+giHDUdkmgjgnDe/iuypYpMv29TTeFJG8hyQqZkzfb6xQ2iEgnAPv/1/b2tUB0HcXOwLosy5ZWUo0tTyZzafT13EwWRQItmrl/zOHznOSNpromGInzj7aZD/7vf3ht8pWcu3Q2//hROUMv+itrex4e2W+A5eu2xCX9m3R2H/of0D6yhmD7rloCRfEO6HatvEfc0f6SxqwriCaZ2gblfctYUDGEzypPZkHFEHUqKwVPtpXCdOAC+/UFwEtR288Xi6OAzYXsT4DE2U/d8JPi2el6Tp16oFgoFvFMV7FtZy3TFlcnXKgGu/0Qpa0CtNu+mUkv38Ejz49nS4tWnH7eRG4dcinB5i3jzgvfP7rzBBoUvdm0PQRirRmIfr6xp8RnI40l3GbpWuSlDmSlKZPJkNRnsJzKe4vIWmAsUAk8KyKXAGuAM+3DZ2KFo67GCkm9KFNyZZNUwhG9kqlZppJ4xRI2YQRDdZFooDI7d0+i/EU1wfiQVa/UFsFdtRy/9E3Gvf539tq5nb8MPJf7fnwWzVq2REJ1riGqsQ52p9F4eHGck6knHDLrROxCu8aac9SBrDRlMhl9dK7LruMcjjXAVZmSJd+Jtku3LQm4JriLrb8cPjc24kbwLoATS3QUU+y6hWg6bvmWW167nxNWv8uSTt24/sRrWNmhK+1aBTj58E6OVdvCxI6yXdNSB0ORmUuYsFyxzwrOC+0aa8JJ10pkRSlE8sXR3GSJ7ehqgiGrspdDrYMXFlXT/4D2CUfcqcSTRXfS4Y414lA2hnOWzuaG+Q8TqK/jlsEX83D/06gvskw1xliyeRE7yvaaEbmF7WarxGQm7pNth7SipIoqhRyQaIFVqN5QH4oPrXRa55AuO7eTaWT00O787aHZjHvlbo5e8wH/2r8XFcN+wxft9mtwXCITVaAovujO6KHdU8oimq0Vwum8TzoWMoa/M+d02cKNlfNUqSgZQ7OkZplwBxF2sCZb4L26Jtgg42e67Nzbd9U2zCBaV0f5/KnMfPgqDt/wX24YejW/POfWOIWQiNKSABPP7B3XgZX3LXONLNrTbPfJRDM5Ef2dgYYlQxUl3ahSyDJ+Q069ooAMVscwcuoSWjUvclzkkSzhHEnTFlfDhx/C0UfDqFE0O/542qz+mD+/em8ko6lfykpLWDL2Z64jWqfIoj3Rdt/YaKbGKhVFSQY1H2UZPx1BSaCYI/Zvyzv/3ejpHzDAqq+3eV5LgFY+azHX7tjBt6NvhDefhrZt4Zln4OyzIwnsypIoDeqnc0/Fdl+ItvnGRjNpiKySTVQpZBmvxHL1xkTKV76wqDolh3EsBthVW5+wZGfvdSu5bdY9HPLtF7zW+zhq75rESUN6RfZPW1zN9l3uJUGjEcTXmgxIznaf7iSD2aKx0UwaIqtkEzUfZRm3BVZ3ntU7srBr/sffJF0C04tQvaF182YNFtKVBKyPvmVoBzfOe4gXnxxN2x1bufj0mxgxbCSj5ldHbNbhzji2iH2rQJHjKuTO7Usy0kkXqhkl1YWMYbTympJNdKaQZfyYTDJhFtgcDLFk7M8i7/v+6TX6rF5C5av3cEDNVzzVZxiVgy5iS4vWQMNIJzc/SLvWLRg9tHvcs5RuXpV2+cGfGSVfzUuNiWbSymtKNlGlkAMSdRBeMfyp0sDUsHkzo1+8i18unc3npZ0459w/s3D/w+POCXe2ieo4xz5LVVVmlIJbu7QtCUTWVESn8S4U85IftPKaki3UfJQjvArKJ0pOl4jYaKQGpoaXXybYrTtnfzCHBwYMZ9jF9zoqBLA614GV8yjNk9BRx/xORcK2XbURZRHrNUnVvOT1+SjKnozOFHLAtMXVjH5uKSG7XkJ1TZDRz1nFbNxq+3ZpX8fnlSc0MI+4uY0Nlt26gQmlc3PWDiun8+yX+LxDV677vztZ1qlbQlmra4IEiiTOUZ0Lm7ZTu2zfVRvn64glnOE1mQinQnRoK0o6UKWQA8ZNXx5RCGFC9YZx05e71vadNmtOg45t0tl9XJPElZWWRDKRYgw8/TQ7h/yGfb7bwp3H/IoHjjqDULH/IjChekNpSYDWLZrl3FYf2y4HVsxIeE50Lig/HXw6quYpSqGiSiEHuKWFcNs+bXE1azcGqa6xTCfhmcXZA7rwwqJq91DHL7+EK66AGTNY3eVQrh1+C6s6HOAqVzgs1mkGEuuojpYtztHseof0k8j/Eu1jCJOog9d1AUpTRn0KBcC46cvjSmmG6g2vLF3vHOrYuxM88AD07Anz58OkSZx6TqWnQgArtYabn6CtQ3nJ2JQd4VF4olxI6cTJzxD2qZSVlria2Lw6+HQV61GUQkSVQg5wy/njtt1rZhFX+avNdhgyxJohDBgAy5bBtdeyb/vEtVzDoY6xaw8AtsXmRsLdzPLlxu1Zc846rQGYdHYfPrfboyyFDl7XBShNGVUKOWDsKT0JFDfseAPFwthTeqZ+0dpamDgRDj8cliyByZNhzhz4wQ+AxBFN4U6vvG8ZbVrGWxVDdSYuisdrtJ3NpG1eJTFT6eAbu9hMUQoZ9SnkgGRz/lgziPgUE5GZxdKlcMklsGgRnHYa/O1vsN/ubKZOldnatQpgjOUriL1/jUs0T6wSSGTPzwfnbKq1EbKVoltR8g1VCjkimU5n7Ck9WbtiUYNtgWJh/NCD4Y9/hMpKaN8enn0WzjiDaUvWMfHxeZFKbtt21UbCSeuMoSRQzNhTerre32+uHaecPrHkg3NWO3hF8Y+ajwqA8r5ldG5X0sCc8VC3Wk696Odwyy1w7rnw0Udw5plMW7KugfO3JhiKS4QXvaDLaZGWX5NLtJnFDXXOKkphoTOFPMQxzLMkwIKKQbBtG9x4I9xwD3TuDDNnwoknRs71W69hXU3QdZHWhOG9mDC8ly+TSzL1kxVFyX9yohRE5HNgC1AH1Bpj+otIe2Aq0BX4HDjLGLMpF/LlEteO+uhieP11+PWv4fPP4corYcIE+N73Gpzv11yzX2mJ5yKtWIdtIjRpm6LsGeTSfDTYGNPHGNPffl8BzDXGdAPm2u+bHE4ddWDLZg65/XY44QQIBODNN+G+++IUAvgz14RH8OlepBWOAupV1jZppaIoSn6QTz6F04DH7NePAeU5lCVnxHbIP/vkX8yZfCU/fHMeVFRYkUY/+Ynrwf1lZAAADHhJREFU+W5J49q1CsSFV+bLIi1NPqco+YMYlwLxGb2pyGfAJqwMBH83xjwoIjXGmNKoYzYZY9o5nDsCGAHQsWPHflOmTElJhq1bt9KmTeIFXdlm5Vdb2FVXT8nmGn762IN0+/c7fLN/V5b+5mqa9e/n6xo1wRAbNu9gV109zYuL6Ni2JaUOK5JrgiGqNwWpj/oOFIlQ1q7E8Xi/hNvWjxyZkiEVeQuBQpIVCkveQpIVGifv4MGDF0VZaRqQK0fzQGPMOhHZB5gjIh/7PdEY8yDwIED//v3NoEGDUhKgqqqKVM/NJDXvr2Xh+L9Q8drfKQnt4PZjz+eJgWdyc88WlGdA3kwUpamqqqKmbTfGzF1GMFREeEJaEqhjwvAeDa5v1UGIX1RXVlpsOdazQL5+F5woJFmhsOQtJFkhc/LmRCkYY9bZ/78WkX8CA4ANItLJGLNeRDoBX+dCtpyyZg3lN15G+auvsmz/How84TcED/4hN0dVM0t3J56pGH6/mUY1+Zyi5BdZ9ymISGsR2Sv8GvgZ8CEwHbjAPuwC4KVsy5Yz6ustx3HPntS+8SaTTrma8rMnEDz4hw06fbcEdPlog/fb2eeLX0NRFItczBQ6Av8UkfD9nzbGvCoi/wGeFZFLgDXAmTmQLfusXAmXXgpvv83XRx3Luf0v5r+t9wYa5v4vpbDy/DdmVbSub1CU3JF1pWCM+RTo7bD9f8Bx2ZYnZ4RCcOedMG4clJTAI4/wi/VdqN68o8Fh4U7/1qOKCsrU4rezTzU3kaIomUFXNOeCxYutBHaLF8Pw4ZbpaN99WedSRczq9Fv7Hn3nA8l09pqbSFHyB1UK2WTHDrj5ZrjtNth7b3j+eTj99MjuRJ1+oZlatLNXlMIjnxav7dksWAB9+sCf/wznnWclsItSCJA497/m+VcUJdPoTCHTbN0KN9wAf/0r7L8/vPoqDB3qeKiXyaWqalXkGFUCiqJkClUKmWT2bLjsMlizBq6+2polJFiBqJ2+oii5RM1HmWDjRrjwQhg2DFq2hLfegnvuSagQFEVRco0qhXTzwgvQowc8+aRlNlqyBAYOzLVUiqIovlDzUbpYv94yEb34IvTta/kO+vTJtVSKoihJoTOFxmIMPPqoNTuYMcOql/zuu6oQFEUpSHSm0Bg+/xxGjIA5c+CYY+Chh6B7fq4ZUBRF8YPOFFKhrg7uvRcOOwz+9S9rRfIbb6hCUBSl4NGZQrKsWGElsHvnHSu66IEH4IADci2VoihKWtCZgl9CIbj1VstX8PHH8PjjMHOmKgRFUfYodKbgh/ffh4svtuojn3mmZTrq2DHXUimKoqQdnSl4EQxCRQUMGAAbNljhps8+qwpBUZQ9Fp0puPHWW5bv4JNPrDTXEydCu3a5lkpRFCWj6Ewhlu++g6uugmOPhV27rHDThx5ShaAoSpNAlUI0s2ZZYab33w/XXgsffgjHH59rqRRFUbKGKgWA//0Pzj8fTjrJSlq3YAFMmgStW+daMkVRlKySd0pBRIaJyEoRWS0iFRm9mTGW4/jQQ+GZZ+CPf7RKZP74xxm9raIoSr6SV45mESkG7gNOANYC/xGR6caYj9J9r+bffmvVR542Dfr1g9dfh8MPT/dtFEVRCoq8UgrAAGC1MeZTABGZApwGpFcpzJzJgAsvtNJV3H47jBwJzfKtKRRFUbKPGGNyLUMEETkDGGaMudR+/3/AkcaYq6OOGQGMAOjYsWO/KVOmJH2fkupqDrzrLj4bOZJg587pET7DbN26lTYFUqSnkGSFwpK3kGSFwpK3kGSFxsk7ePDgRcaY/k778m14LA7bGmgtY8yDwIMA/fv3N4MGDUrpRlVlZaR6bi6oqqoqGHkLSVYoLHkLSVYoLHkLSVbInLz55mheC3SJet8ZWJcjWRRFUZoc+aYU/gN0E5EDRaQ5cA4wPccyKYqiNBnyynxkjKkVkauB2UAx8LAxZnmOxVIURWky5JVSADDGzARm5loORVGUpki+mY8URVGUHKJKQVEURYmgSkFRFEWJoEpBURRFiZBXK5qTRUS+Ab5I8fS9gW/TKE6mKSR5C0lWKCx5C0lWKCx5C0lWaJy8BxhjOjjtKGil0BhE5D23Zd75SCHJW0iyQmHJW0iyQmHJW0iyQubkVfORoiiKEkGVgqIoihKhKSuFB3MtQJIUkryFJCsUlryFJCsUlryFJCtkSN4m61NQFEVR4mnKMwVFURQlBlUKiqIoSoQmqRREZJiIrBSR1SJSkWt5YhGRz0VkmYgsEZH37G3tRWSOiKyy/7fLoXwPi8jXIvJh1DZH+cTiHrutPxCRI/JA1nEiUm237xIROSlq3xhb1pUiMjSbstr37yIi80VkhYgsF5Fr7O15174esuZl+4pISxH5t4gsteUdb28/UETetdt2qp22HxFpYb9fbe/vmgeyPioin0W1bR97e/q+B8aYJvWHlZL7v8APgObAUqBHruWKkfFzYO+YbbcDFfbrCuC2HMp3LHAE8GEi+YCTgFlYVfWOAt7NA1nHAb93OLaH/X1oARxof0+KsyxvJ+AI+/VewCe2XHnXvh6y5mX72m3Uxn4dAN612+xZ4Bx7+wPAFfbrK4EH7NfnAFPzQNZHgTMcjk/b96ApzhQGAKuNMZ8aY3YBU4DTciyTH04DHrNfPwaU50oQY8ybwMaYzW7ynQY8biwWAqUi0ik7krrK6sZpwBRjzE5jzGfAaqzvS9Ywxqw3xrxvv94CrADKyMP29ZDVjZy2r91GW+23AfvPAEOA5+3tsW0bbvPngeNExKlkcDZldSNt34OmqBTKgC+j3q/F+4ucCwzwmogsEpER9raOxpj1YP0YgX1yJp0zbvLla3tfbU+zH44yxeWVrLa5oi/WKDGv2zdGVsjT9hWRYhFZAnwNzMGardQYY2odZIrIa+/fDHw/V7IaY8Jte6vdtpNEpEWsrDYpt21TVApOmj7f4nIHGmOOAE4ErhKRY3MtUCPIx/a+HzgI6AOsB+60t+eNrCLSBngBuNYY853XoQ7bsiqzg6x5277GmDpjTB+s+u8DgEM9ZMqpvLGyishhwBjgEOBHQHvgevvwtMnaFJXCWqBL1PvOwLocyeKIMWad/f9r4J9YX94N4emg/f/r3EnoiJt8edfexpgN9g+uHvgHu00YeSGriASwOtmnjDEv2pvzsn2dZM339gUwxtQAVVj291IRCVehjJYpIq+9vy3+TZFpI0rWYbbJzhhjdgKPkIG2bYpK4T9ANzvioDmWA2l6jmWKICKtRWSv8GvgZ8CHWDJeYB92AfBSbiR0xU2+6cD5dnTEUcDmsBkkV8TYWn+B1b5gyXqOHXVyINAN+HeWZRNgMrDCGHNX1K68a183WfO1fUWkg4iU2q9LgOOx/CDzgTPsw2LbNtzmZwDzjO3VzZGsH0cNDATL9xHdtun5HmTLm55Pf1ie+k+w7Ik35lqeGNl+gBWhsRRYHpYPy5Y5F1hl/2+fQxmfwTILhLBGKJe4yYc1rb3PbutlQP88kPUJW5YP7B9Tp6jjb7RlXQmcmIO2PQZr2v8BsMT+Oykf29dD1rxsX+BwYLEt14fATfb2H2App9XAc0ALe3tL+/1qe/8P8kDWeXbbfgg8ye4IpbR9DzTNhaIoihKhKZqPFEVRFBdUKSiKoigRVCkoiqIoEVQpKIqiKBFUKSiKoigRVCkoSgJE5BciYkTkkATHXSgi+zXiPoNE5JVUz1eUdKBKQVEScy7wNtZCRy8uBFJWCoqSD6hSUBQP7Lw+A7EWvZ0Ttf06sWpeLBWRShE5A+gPPGXnuS8Rqy7G3vbx/UWkyn49QETeEZHF9v/u2X8yRXGmWeJDFKVJUw68aoz5REQ22sVLOtrbjzTGbBeR9saYjSJyNVYdgXBhJLdrfgwca4ypFZHjgT8Dp2f+URQlMaoUFMWbc4G/2K+n2O+LgEeMMdsBjDHJJklrCzwmIt2w0kQE0iSrojQaVQqK4oKIfB+rAMthImKwqvYZrKygfvLD1LLbRNsyavvNwHxjzC/sOgRVaRJZURqN+hQUxZ0zsKpZHWCM6WqM6QJ8hpU++WIRaQVW/WT7+C1YZSnDfA70s19Hm4faAtX26wszI7qipIYqBUVx51ysehbRvIAVYTQdeM+ujPV7e9+jwANhRzMwHrhbRN4C6qKucTswQUQWYM0+FCVv0CypiqIoSgSdKSiKoigRVCkoiqIoEVQpKIqiKBFUKSiKoigRVCkoiqIoEVQpKIqiKBFUKSiKoigR/h/KZUGP56IbqQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(y_tr,y_tr_pred)\n",
    "plt.plot([0,350],[0,350],'r')\n",
    "plt.xlabel('Actual')\n",
    "plt.ylabel('Predicted')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Validate on Test Data\n",
    "\n",
    "As we will explain in the next lab, one should not evaluate the performance of a model on the data on which it is trained.  It is important to see how the model works on *new* data that is independent of the training data set.  For this reason, we will evaluate the model on the *test* samples that were not used in training."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RSS per sample = 0.492801\n",
      "R^2 =            0.507199\n"
     ]
    }
   ],
   "source": [
    "X_test = X[ns_train:,:]\n",
    "y_test = y[ns_train:]\n",
    "y_test_pred = regr.predict(X_test)\n",
    "RSS_test = np.mean((y_test_pred-y_test)**2)/(np.std(y_test)**2)\n",
    "Rsq_test = 1-RSS_test\n",
    "print(\"RSS per sample = {0:f}\".format(RSS_test))\n",
    "print(\"R^2 =            {0:f}\".format(Rsq_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the model predicts new samples almost as well as it did the training samples.  We can also plot the fit on the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEGCAYAAACevtWaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deXhU5fXHPycxQCBoUDFioO6ioBUkxQW1BHetQt1tf3UpFitQlSolLq1QtyguFaUiivsCVDSioKhAXFBQMeyCorgQEFQIEgkQkvf3x70TJ8md/c7MvZPzeZ48mbnr974zc+655z3vecUYg6IoipJZZKVbgKIoiuI+atwVRVEyEDXuiqIoGYgad0VRlAxEjbuiKEoGslO6BQDsvvvuZp999olr359//pl27dq5KyiJ+Emvn7SCv/T6SSv4S6+ftEJieufPn/+DMaaj40pjTNr/evXqZeJl9uzZce+bDvyk109ajfGXXj9pNcZfev2k1ZjE9AIfmxB2VcMyiqIoGYgad0VRlAxEjbuiKEoGosZdURQlA1HjriiKkoFENO4i0kZEPhSRhSKyVERG2cufEJFVIrLA/uthLxcRGSMiK0VkkYgckeyLUBRFURoTTZ77NqCfMaZaRHKA90TkNXvdcGPMC022Pw040P47EnjI/q8oiqKkiIieu51OWW2/zbH/wtUJ7g88Ze83F8gXkU6JS1UURckgamuhtJT2y5cn5fBioqjnLiLZwHzgAGCsMWaEiDwBHI3l2c8ESowx20TkVaDUGPOeve9MYIQx5uMmxxwEDAIoKCjoNXHixLguoLq6mry8vLj2TQd+0usnreAvvX7SCv7S6weteZ9/TtfRo2n/+eesPOccVg8dGtdxiouL5xtjihxXhhrd5PQH5AOzgUOBToAArYEngX/Z20wDjg3aZybQK9xxdYSqN/GTVmP8pddPWo3xl15Pa62pMeaGG4zJzjamoMCYF17wxghVY0wVUA6caoxZax9/G/A40NvebDXQJWi3zsCaWM6jKIqSccyZAz16wO23w8UXw6efwjnnJO100WTLdBSRfPt1LnAisDwQRxcRAQYAS+xdpgIX21kzRwGbjDFrk6JeURTF62zeDH/7Gxx3HGzdCjNmwGOPQYcOST1tNNkynYAn7bh7FjDZGPOqiMwSkY5YoZkFwF/t7acDpwMrgS3AZe7LVhRF8QEzZsCgQfDtt5aBv+02SFF/QETjboxZBPR0WN4vxPYGGJK4NEVRFJ+yYQP8/e/w5JNw8MHw7rvQp09KJegIVUVRFDeZMgW6dYNnnoEbb4SKipQbdvDIZB2Koii+Z+1aGDoUXnwRjjgCXn/d6kBNE+q5K4qiJIIx8Pjjlrc+bRqUlsK8eWk17KCeu6IoSvx89ZXVYfrmm1Y2zKOPwkEHpVsVoJ67oihK7NTVwZgxcOih8MEHMHYslJd7xrCDeu6Koiix8emncPnl8P77cOqp8PDD8KtfpVtVM9RzVxRFiYbaWitPvUcPWL4cnn4apk/3pGEH9dwVRVEiM38+DBwICxfC+efDAw/AHnukW1VY1HNXFEUJRU0NlJTAkUfC+vXw0kswaZLnDTuo564oiuLMO+9YsfXPP7e89rvvhvz8dKuKGvXcFUVRgvnpJxgyBH77W9ixA956y0px9JFhBzXuiqIov/Daa1Z640MPwTXXwOLFcMIJ6VYVFxqWURRF+fFHGDbMyoDp1s1KczzqqHSrSgj13BVFabkYA5MnwyGHwPPPw7/+BZ984nvDDuq5K4rSUlmzBgYPhpdfhqIiK7b+61+nW5VrqOeuKErLwhiYMMEKv8yYAaNHWyUEMsiwg3ruiqK0JL78Ev7yF5g1y8qGefRROOCAdKtKCuq5K4qS+dTVwX/+A4cdBh99BOPGWQY+Qw07qOeuKEqms3SpNQhp3jw44wzLsHfunG5VSSei5y4ibUTkQxFZKCJLRWSUvXxfEZknIp+LyCQRaWUvb22/X2mv3ye5l6AoiuLA9u3w739Dz57wxRfw3HPwyistwrBDdGGZbUA/Y8zhQA/gVBE5CrgTuM8YcyCwERhobz8Q2GiMOQC4z95OURSPUlZRSZ/SWexbMo0+pbMoq6hMt6TE+egjKwPm5pvh3HNh2TK46CIQSbeylBHRuBuLavttjv1ngH7AC/byJ4EB9uv+9nvs9SeItKAWVRQfUVZRyfUvLqayqgYDVFbVcP2Li/1r4LdsYb9x46w89Q0bYOpUy2Pv2DHdylKOGGMibySSDcwHDgDGAqOBubZ3joh0AV4zxhwqIkuAU40xq+11XwBHGmN+aHLMQcAggIKCgl4TJ06M6wKqq6vJy8uLa9904Ce9ftIK/tLrFa0rvtvM9rr6ZstbZWfRdc/2De+9ojcc+QsWcNDdd9O2spI1v/sdX1xxBXUe1wyJtW1xcfF8Y0yR07qoOlSNMXVADxHJB14CDnHazP7v5KU3u4MYY8YD4wGKiopM3759o5HSjPLycuLdNx34Sa+ftIK/9HpF62Ul0zAOD/ACrCrt2/A+Eb1lFZWMnrGCNVU17JWfy/BTujKgZ2F8gp3YtAlGjLBmRNp/fxbcey89hg1jL/fOkFSS9V2IKRXSGFMFlANHAfkiErg5dAbW2K9XA10A7PW7ABvcEKsoirvslZ8b0/JYSXrY59VXoXt3eOQRuO46WLSIqp493Tm2z4kmW6aj7bEjIrnAicCnwGzgXHuzS4CX7ddT7ffY62eZaGI/iqK4RrSdpMNP6UpuTnajZbk52Qw/pasrOkbPWEFNbV2jZTW1dYyesSKxA3//PfzhD3DmmdChgzXCdPRoaNs2seNmENGEZToBT9px9yxgsjHmVRFZBkwUkVuBCmCCvf0E4GkRWYnlsV+YBN2KooQg4C0HjGrAWwaahUMC75MVNllTVRPT8ogYAxMnwlVXWeGYUaOsmZJatUpAZWYS0bgbYxYBzZ5zjDFfAr0dlm8FznNFnaIoMRPOW3Yy2gN6FrobAw9ir/xcKh0MeVxhn9Wr4corrVBM795WfZhDD3VBZWai5QcUJcNw3VtOAFfCPvX1MH68FVufORPuvdeqt66GPSxafkBRMgxXveUEiSXs45hV077GKvRVXg79+llGfv/9U3wV/kSNu6JkGMNP6doo5g7udpLGSjRhn6b9BN9tqOaz4TdT994zZLdubWXDDBzYokaYJooad0XJMJLdSZoMgvsJun7/FXe+dj891n7Ou4ccw3FvToZC72r3KmrcFSUDSWYnaTJYU1VDqx21DPlgMoPnTmZTmzyGnvUPph18HKvUsMeFGndFUVJGqNGqJ/70Fdf97y66/vANL3Yv5pZ+l7Ox7S4UpqGfIFNQ464oSkpwyr//98SPOPSeUYx/7lG+y9uNy869mdn7/wZIbz9BJqDGXVGUlNA0//7orxdS+voD7F31HVx5JfMvuorP5qxBfNJP4HXUuCuKkhICefY7b63m+tmPcdGiN/iyw15c8IdSJv13BL8DfnfcwekVmUGocVcUJSXslZ9Lt49mc+sb/2X3n6sYd+Q53NfnD+zeMT/d0jISNe6KoiSf9euZNPt+Or8xlU877sPlZ/+TxZ0O1Lh6EtHyA4qiJA9j4Jln4JBD6Fz+OssGD+eKqx9mSacD6dA2h9Y7ZTFs0oLMmd7PQ6jnrqSdpE/m0MJJW/t++y389a8wfbo17d2ECXTr1o13iK1ypZu0pO+aGnclraTrR95SiKZ9XTd49fXWrEj/+If1+v77YcgQyP6lgFislSvdoKV91zQso6SVpE3moACR29f1mZI++wz69oXBgy1vfckSq/Z6duPKkE6FzSC5lStb2ndNjbuSVrxUnjYTidS+rhm8HTvgrrvg8MNh8WJ47DF44w3Yd99mm5ZVVDpOtAzJrVzZ0r5ratyVtJLsOTxbOpHaN5RhC+VZO7JwIRx5pDVJ9WmnwbJlcNllISs4jp6xAqd5NwWSkjkTmHIw1FyfmfpdU+OupJVkz+HZ0onUvqEMm0Dk0My2bfDPf0JRkTVL0v/+B1OmQKdOYXcLdUMxuB/7Dg47OZHJ3zU17kpaGdCzkDvOPozC/FwEKMzP5Y6zD8vIDq50EKl9h5/S1TFEYiB8aOb996FnT7j1VvjjHy1v/dxzo6q3HuqGkowiYU5hp+DzZfJ3TbNllLTjt/K0fiNc+w7oWcg1kxY4rnP0sKur4cYb4YEHoEsXeP11OOWUmPSkcjKRUE8JAswp6ef6+bxERM9dRLqIyGwR+VRElorI1fbykSJSKSIL7L/Tg/a5XkRWisgKEYntk1cUJaWE8pibedhvvgmHHQZjxlipjUuWxGzYIbVPay25Tycaz30HcK0x5hMRaQ/MF5E37XX3GWPuDt5YRLoBFwLdgb2At0TkIGOM87ORoihpJaInvXEjXe+80/LSu3aFd9+FY49N6Jypelrz2pSDqSSi526MWWuM+cR+vRn4FAj3qfQHJhpjthljVgErgd5uiFUUxX3CetIvvQTdurHnG2/A9dfDggUJG/ZU0pL7dMSYUAlCDhuL7AO8AxwK/B24FPgJ+BjLu98oIg8Cc40xz9j7TABeM8a80ORYg4BBAAUFBb0mTpwY1wVUV1eTl5cX177pwE96/aQV/KXX61pbbdjAAWPGsMfbb7P5gAOoGDKE+h490i0rKrzetk1JRG9xcfF8Y0yR40pjTFR/QB4wHzjbfl8AZGN5/7cBj9nLxwL/F7TfBOCccMfu1auXiZfZs2fHvW868JNeP2k1xl96Pau1vt6YJ54wpkMHY1q3Nub2243Zvt27eh3wk1ZjEtMLfGxC2NWosmVEJAeYAjxrjHnRvimsC1r/CPCq/XY10CVo987AmmjOoyhKGvn6a7jiCpgxA/r0gUcfhYN18gy/Ek22jGB5358aY+4NWh48UuH3wBL79VTgQhFpLSL7AgcCH7onWVEUV6mvhwcfhO7d4b33rDTHd95Rw+5zovHc+wB/AhaLSCAh9gbgIhHpgTXe4SvgCgBjzFIRmQwsw8q0GWI0U0ZRvMmKFTBwIMyZY6U1Pvww7L13TIdoSWV0/URE426MeQ8cB7FND7PPbVhxeEVRvEhtLdx9N4waBW3bwhNPwMUXRzXCNJhMLaObCTcsLT+gKC2Nigro3RtuuAHOPNMqHXDJJTEbdsjMMrqul0FOE1p+QFFaClu3Wp766NHQsaNV5Ovssx03bea5Hu4cWc3EMrrpmEgkGajnrigtgffes2qtl5byao+TOPyC/9Dns3xHb9TJc63cWOO4bSYO78+UG5Yad0XJZDZvhqFD4bjj+Lm6hoF/uI2hJw5lU5u8kOEGJ8+13hjHUEsmlmzOlBuWhmUUpQmJdqZ5pjNuxgwYNMiaqPqqq+i/68msbOJ81tTWMXLq0kb6ovFcg68xv20OrXfKYlNNbdTX65k2ciBT6tGocVeUIBLN/vBE9siGDTBsGDz1lJWr/t57cMwxfFEyzXHzqppayioqG/TtlZ/rOLlFwHNteo0bt9SSm5PNfRf08E8bhSGgwas3n2jRsIyiBJFo9kdas0eMgRdegEMOgeees+quV1TAMccA4cMKwfqcQi1ZIg2eq6/bKEoG9CxkTkk/VpWewZySfr4z7KDGXVEakWhnWto649auhXPOgfPOg86d4aOPrFmS2rRp2CRcWCFYn1MlxcIOuQ0Gzrdt1MJQ464oQSTamZbyzjhj4PHHoVs3eO01uPNOmDcPHCo4DuhZSIe2OVHpa+q55ufmhNw20vJot/Nbh6XXUeOuKEEkmv2R0uyRVavg5JPhz3+2ZkhauBD+8Q/YKXRX2s1ndk9Yn6/aqAWjHaqKEkSinWkp6Yyrq4OxY63JM7Ky4L//tao5ZkX21dzQ54s2UtS4K0pTEp0CLqlTyC1bBpdfDh98AKedBuPGwa9+lXJ9nm4jBVDjrrQwvJxfHZbaWiuefsst0L49H99yP1dnd2fNfxezV/5K/1yHkjLUuCtJxUvG1Ov51SGZP9+Kqy9aBBdcwPS/XM+15Wupqd0K+Og6lJSiHapK0vBadT0/5Fc3oqYGRoywKjh+/z2UlcHEidz20Y8puY6yikr6lM5i35JprPhus++qIrZ01LgrScNrxjTd+dUxGct33rEKfd11l+W1L1sG/fuH1evmdTS9MW+vq/dl2duWjBp3JWaCjVSf0lkhf/ChjI3T0PZUkM786pvKFjNs0oLIxvKnn2DwYPjtb2HHDnjrLXjkEcjPj6jXzevw2o1ZiR017kpMxBJqCWVsxD5OqklXfnVZRSXPzv0G02R5M2M5fbo1j+m4cVZtmMWL4YQTmh1v+CldycluPLFGTra4eh3pfspREkeNuxITsXh0w0/p6jg/o4G0eIBOw+rvOPuwpHdCjp6xoplhD7CmqgZ++AH+7//gjDNg553h/ffh3nuhXbvQB21ywNo6w8dfb3BNs44i9T8Rs2VEpAvwFLAnUA+MN8bcLyK7ApOAfbAmyD7fGLNRRAS4Hzgd2AJcaoz5JDnylVQTi0c3oGch10xa4LB1+jzAZOZXh8oMCnmtxvCnb+dBt8tg40a4+WZrYFLr1mHPM3rGCmrrm98unp37DUV77+rK9Xm57K2XMrC8TDSpkDuAa40xn4hIe2C+iLwJXArMNMaUikgJUAKMAE4DDrT/jgQesv8rGUCkcrBNKYxxe7/ilGZ5zaQFjJy6lPy2OWzcUtto+3YbN/DIiw9z0sp5UFQEM2daJQSiINTNIvBE1NTQxWMMm44ibZWdlfSnnGh0+jadNQ1EDMsYY9YGPG9jzGbgU6AQ6A88aW/2JDDAft0feMpYzAXyRaST68qVtBBr3Lql1BFxCleBVSu9euuOX2LkxnDBwhn8ccRQ+n6zAO6+2xptGqVhh/A3xqaGP5F01ODiYV33bJ+w8QzXER+tTu3ojR4xJlQ00GFjkX2Ad4BDgW+MMflB6zYaYzqIyKtAqTHmPXv5TGCEMebjJscaBAwCKCgo6DVx4sS4LqC6upq8vLy49k0HydJbVVPLuk1b2V5XT6vsLAp2adOokl88hNIa67mSoS0WvalgceWmsOuzs4QO67/j2PEP0mXZYn449FC+KCmhpjB2g1lVU8u3G7Y4rmuVnUXXPds3vF/x3Wa219VH3C4SibZtVU0tlRtrqA+yN1kiFHbIJT83J2qd4dr5sMJdXNGaahLRW1xcPN8YU+S0LuoRqiKSB0wBrjHG/GSF1p03dVjW7A5ijBkPjAcoKioyffv2jVZKI8rLy4l333SQDL1lFZVcP3MxNbVZBB7GcnPquOPsbgl5W9q20XNj6ayQKZ5Z9XX8+eOpDJ77nFWx8eGHWXLAAfTt1y/u891UtrhZBk5uTjZ3nH0YfYM+88tKpmEcHtAFWFXaN+rzJdq2fUpnUVmV3Wx5YX42c0r6Rq0zVDsX5ufytz/2dUVrqkmW3qiyZUQkB8uwP2uMedFevC4QbrH/r7eXrwa6BO3eGVjjjlzFCX1UTT9O4SeAg77/ihefGc5NsydYaY3LllnzmkZRwTEctw44jPsu6BEx88crWS+ROuKj1dlSwnxuEE22jAATgE+NMfcGrZoKXAKU2v9fDlo+VEQmYnWkbjLGrHVVtdIIr+ck+z27IRr9gfejXlnKxi215NTVMviD/zHkg8lUt2nLR7c/yG9KBkPoJ96YzxtN5o9Xsl4idcRHq1PLBUdPNGGZPsCfgMUiEshruwHLqE8WkYHAN8B59rrpWGmQK7FSIS9zVbHSjFgzWFKJ17MbwhnQsopKRk5dSlXNL5ku4fQHjG35M9Poct3f2X/dKt44/AR23Hsvp/f7dcy63Gq31jtlNRynQ9scbj6ze8rbPpLxjsVoa7ng6Iho3O2O0VDuRrPhc8bqoR2SoC4lBrzinTkRLmSU7h9oOAMKNGvTACH1b9kC//oXfe+7Dzp1gqlTOfnMMxudL2C8SnrUU1VRGbIN3Gi3ptcHsLW2eadlKojGeKvRdhct+ZsBePlR1csho0h9FU6GPUAz/bNnw1/+Al98Yc2KdOedsMsuDaubGtpAbRlw9sTdaDev3VjVeKcWNe4Zgld/OF4OGSViQBv0b9pkzVs6fjzsvz/MmgXFxc22j9XQutFuXr6xKslHa8soScXL2Q3hMjTCGdEG/a+8At26waOPwnXXWZNpOBh2iN3QutFuXsmUUdKDGnclqaSrWFc4AiMlK6tqmnUmBQxoqNTGDm1zuKd4LwaMHg5nnQW77QZz58Lo0dC2bchzxmpo3Wi3VN9Yoy0FraQGDcsoMRNramMiISO30yibxr4NVraAwTKgTY/f6NwnH8SA5e/A+X+w6q6PGgUlJdCqVcTzxtPp7cYk1M2uIUl9MV7PimqJqHFXYiKVP+JI2SxNjVZ+qAMF4RT7Dhj2OSWNR4w2Mq6rV8OVV8Krr8KRR8KECVbt9ShxoxBXvAXAUmFcvdZ5q6hxTzlVNbX0KZ3luayWaEnlj3jUK0sdzzXqlaVsra1vZvTvOKZ5GKUpMXcy1tdbMyENH86O2lr+e8aV3N/tVPZ8ZR3Dt+fHdM3Bhra8vLxRmYBIpOKmmshTkp86b/0+qC5a1LinkLKKSio31jTU2PDjo2u4qfPKwuRtx0pZRWWzMrkBnJbX1NaxbpPz9sFEm4VSVlHJ88/OZNjk0Rz17RK+OKw3fz3+r3yetweQ+s8u2TfVeG4ewUYyS4Q6hyKEXuu8bUnhI+1QTSGjZ6xoVBUP/FcDJtyP1c0JlONpE6eqgk2JppPx5Y++5rPhN/Pk/X+h2/pV/OPUqzjxtH82GPYATp9dsjoVk+kZl1VUcu3khTHVJ2paotfJsHslKyqYllSHSY17CvHTo2soQmWRgLs/knBtEqpccKvsyF/niFkoixZx0O9P4R8zJ/DOvkdw4sD/MvnwkzEhasIE60ykdnokkpXWGNDsZJwh9OcQqn59tohnsqKcyITfYLRoWCaFWD/EzSGWp5Z4446BbZI9fV6o8El+bg4jz+rumHlSsEvkrBUI0cm4bRvcfjvcfjsdW7VjyFkjmHbwsRELfe0SdKNJJHQS6fNIVomJUEY6QKjvZqjPud4YVpWekZCmZOLlQXVuo557Chl+SleymhiLdDy6uuFhZocwem79SEKFT0ae1T2k9x33BCBz58IRR8C//w0XXsif/v440w45LqoKjj9v39HQbvF6hdF8HskaLxBOW7jvpl8HSHl5UJ3bqHFPIQN6FlLYITftA3oSiTuGe4x380cSyZgFTwE3p6RffG34888wbBgcc4yVtz5tGjz9NFec3dvRALTNaf5zqa0zDe0Wr8GL5vNIVoZHKG3ZImG/m05GUrBuTF4ewOTFQXXJQsMyKSY/N4c5JX0bLUt1alYiccdwsdZ4fiThrj2pOdozZ1qFvlatgsGD4Y47YOedG84LzfPoh0UIRcUbOon0ebiV4eHU1qE0R/osg9soMNI3cLv3egaKV+swuY167mkmmZ1woUjkkTpcrNVppvpwmSPpuHaqquDyy+HEE60p795+G8aObTDsAYKfDIaf0pXRM1Y0nyvSJtBu8XqFkT4PNzI8QrU10KAZrJt04NiRPodAGxXm5zZrm0zNQPET6rmniYAX5dS5k+yRfYl0zsWSJx7J2wxltEa9sjQ5TzIvv2yNMl2/HkaMgJtvhtzwNzSnmujBNG23eLzCSJ9Hssv/Bkbmxvt00JIyUPyEeu5pINiLCkUyfxiJxB1DpUJu/HlbIw89Gm8z1DVu3FLrrje/bh1ccAEMGAB77AHz5kFpaUTDDuGzSdyK10b6PNzovIxkgBN5OvBr52qmo8Y9DURKP4Pk/zDi7ZAc0LOQc3oVNqumuKW2vpExDnXjCjYy+W2jy26J+xHfGHj6aassb1kZ3HorfPQR9OoV9SFCGUWB+DtyHQh8Hvdd0AOAYZMWNNwoU1H+NxHvuyVloPgJNe5pINIPxus/jNnLvw8ZfwbLGEeTKhli3IwjMT/JfPMNnHEGXHwxdO0KCxbAjTdCTmzpkqn0SqOJiyer/G8i19mSMlD8RMSYu4g8BvwOWG+MOdReNhL4C/C9vdkNxpjp9rrrgYFAHXCVMWZGEnT7mlBxa3AuO+s1ojG0dcaQm5MdNq6/KWji6UhEbUzr62HcOCumXl8P998PQ4ZAduSiYk6kcn7aSHHxZJb/TfQ6W0oGip+IpkP1CeBB4Kkmy+8zxtwdvEBEugEXAt2BvYC3ROQgY0z4GEQLI970s2QRaypmuJtTgMBNKtxxQx0nOK0Oojcyud9+C337wrvvwkknWVPf7bNPxP3Ckcqa6MnumAxngL08D68SHxGNuzHmHRHZJ8rj9QcmGmO2AatEZCXQG/ggboUZiJd+SOGyWkLVR3e6OQUTMMaRvLlQN7lzehUye/n30bfNjh1wzz385p//hHbt4PHH4ZJLohph2pRQN7pUfDbpHhqv3ndmISaKwKdt3F9tEpa5FPgJ+Bi41hizUUQeBOYaY56xt5sAvGaMecHhmIOAQQAFBQW9Jk6cGNcFVFdXk5eXF9e+6cBreld8t9mxmmKr7CwK8ySk1qqaWtZt2sr2unp2yhKrMmC9oVV2FgW7tIm6FEDwcWLdFyBv5Uq63nUX7T//nLVHH82qa69l+267Rb0/WJ7xhp9rMQ49CVkiFHbIjb+0QQicvgdVNbVUbqxpVDk0WeePFa99b8PhJ62QmN7i4uL5xpgip3Xx5rk/BNyC9fR8C3AP8GdolkQBOPe9GWPGA+MBioqKTN++feMSUl5eTrz7pgOv6b2sZBrGoV9dgMdPbecprY3YuhVuuQXuvBN23x1eeIEVu+0Ws96byhbzzNxvgNAx+cL87GajihMl1PfAqxNJeO17Gw4/aYXk6Y3LuBtj1gVei8gjwKv229VAl6BNOwNr4lanJI2AEYk06tKTvP8+DBwIy5db4Zd774Vdd4Xy8pgP9fy8byNu4+aYg0C7X9hlMzeWzmpmvDU0orhFXKmQItIp6O3vgSX266nAhSLSWkT2BQ4EPkxMovfx26zvkQZReTEVs6yikhNHTeOJorOoP/ZYtlRthtdfhyeesAx7nISqYx6MWze6pu2eknILSoslonEXkeexOkS7ishqERkI3CUii0VkEVAMDA8TyMIAABubSURBVAMwxiwFJgPLgNeBIZmeKZOW+igJkopRl25SVlHJK6Of4Il7/8zF81/lqZ5n0PvC+9hn9o6Eb6ah8vEDuHmja0mzACnpJ5psmYscFk8Is/1twG2JiPITfpz1PdKoS0+xcSM5fxnIhPkz+GLXzpz/x1I+7ty9YXU02T3huOjILnbMvTlujzlwK9XRq3F5xVto4bAE8WPRpHSn3EXNiy/CkCGcsm49Y486jzF9LmLbTs1nW6qprePayQu55/jYM0puHXAYYMXe64whW4SLjuzSsNxN3Gj3ljTBs5IYatwTJJofbLCnVdKjnqqKSscfYiiPzG1PLZWjLuPiu+9g6FCYMgV69ODy80ZS3rZz2F3qjKFyYw1lIdo2HLcOOCwpxrwpbrS7H58UlfSgtWUSJFLNjqYx+e119Y4x+VCx+5vKFrse0/dsLRBjrA7Sbt3g1VetCTQ+/JABl/0u5KTcwdQb4+n4dXC7Q3zt7scnRSU9qOeeIJFGm0braYXaLhAuiLR/PLrTbsyD+eoruOIKeOMNOPZYePRR6Nq1UfngbBHqjGlWniAYrxu5QLuXl5fztz/2jXl/34TUlLSjxt0FwhnKaD2tUNuFStXzuhFrSsjQUn29NRPS9ddb5QIefNCaUCMrq1l8OVCM7JxehY43Pch8I+f5kJriGTQsk2SiLaUabqLiWI6bDuKdTu+tF9+G44+Hq66yvPUlS6wKjlnW1zLU08zs5d9zz/mHNwvVZIlkvJHzbEhN8RzquSeBYC91l9wccrKF2rpfvEwnTytcEa0p8ys966nFM53eTnU7uOyDyRx/x/Owcx48+ST86U/NCn2Fe+pxCocVdqhrEUbOcyE1xZOocXeZpsauqqaWnCyhQ9scqrbU0io7y9HTcjJWxQd3ZPby7xvFm71W7z2aPoVgI939u5Xc9doYuq//kuld+3D621OgoMDx2Pltc9i4pXnN9ywR9i2Z1qx/ozyO8gPxoHnmih9Q4+4yTsautt7QttVOVPzrZKtIUJia2gEj4RRvFqD44I6eMiTR9CnslZ/LD99XcfX7zzNo3otsaLsLV/z+Bpb0PoHTQxj2sopKqrfucFwXiLWnI8db88wVv6DG3WXcSlUbOXVps5uEAZ6d+w1Fe+/qGUMSTfbGHXtsosvoq9h3QyWTDjuJ2/oNpLb9LtwRlC7a1BMePWMFtfWR676kOsc7HXnm+qSgxIN2qLqMG3NullVUUhViCjoDnsrlDpvnv3kzDB3K8QPPZo82WVw98C5KTr+a9nt2bAhNhepsjTTTUzCBG2dVTW3SC7ilOs/cj7WLFG+gnrvLuDUKMRxeSoMMmef/3SLofwWsXg1XX027W2/l/rw87m+yfyhPONDHEA175edSVlFJ5cYaKqusG02ywiWpzjPXEalKvKhxdxk3ptCLZLzjNSTJerxvlL3x448wbBg8/TQccgjMmQNHHx1y33D5/U0n2M7JFjA0CtcEbpxWjXT3B3s1JdV55joiVYkXNe5JINFUtXATUOfmZFN8cEf6lM6KyUhH0xHYzPgfHkO1ZmPghResmjAbNsBNN1l/rVvHda2hJtgG5xvnsEkLGk8TY+O2EQy01cipSxtCZ21ykhfd1BGpSryocfcgoSagzs/N4XeHd2qU9x5t+CHS472T8a/cWBddIa61a2HwYCgrg169rBIChx8e97VGmmDbaZll7DaHWO4+23b8Mu/sxi21UX0G0RaQC8YrI1ID/RnaqesftEPVgziNQvzPBT1YcPPJTFu0Nq4JHyI93jsZ/4iFuIyBxx6zwi+vvw533QVz50Zt2MG9EZfDT+lKVpNBUMkygvFMuhFtAbmmeGFE6i/9Gdqp6yfUc/coTl5rWUWl46AeiC5OH+7xPubY7qpVMGgQvPWWVULgkUfgoIPCagiFGyMuB/QspOy7ZRTmZyfdu4wnDp5Ix2i6R6Smqj9DcRc17j6hrKKSaycvDLk+Uvgh0uN91LHdujqruNcNN0B2Njz0kGXks9L/EJifm8Ockr5JP088cXA/d4yuqapJSX+G4i7p/0UqEQk80odLDYwUfhjQs5BzehU2FCLLFuGcXr94hE756s0KcS1bZhX4uuYa+O1vYelS+OtfXTXsfphsPFINfyfcGP+QLvysvSUTzQTZj4nIehFZErRsVxF5U0Q+t/93sJeLiIwRkZUiskhEjkim+JZCuAmtwfJYo8mWmTK/suEGUWcMU+ZXNhhPp9huYYdcBvQs5OUPV/HISZex7deHU7VoGR/fOgamTYMuDu5cAvhlwE48cfB4bgheIZX9GYp7RBOWeQJ4EHgqaFkJMNMYUyoiJfb7EcBpwIH235HAQ/Z/JQHCjdbMzclm5FndHdcFZ2dkOQwKqqmtY+TUpSFz38vLyyl/djpd/z6U/utXMfWQ4xl1wiC2bN2NOxasCWnM4s2n99OAnVjj4E3HP4QqIOdFUtmfobhHRONujHlHRPZpsrg/0Nd+/SRQjmXc+wNPGWMMMFdE8kWkkzFmrVuCWyLhRmuGMhBOhcecqKqpbcjXbpRWefCu7DduHIWT/8f37fK5/Ox/8taB9n06jMFNpLCWn+PS0RB8QwhXQM6LpKo/Q3EPMVEM8baN+6vGmEPt91XGmPyg9RuNMR1E5FWg1Bjznr18JjDCGPOxwzEHAYMACgoKek2cODGuC6iuriYvLy+ufdNBPHoXV24Kue6wwl0cl6/4bjPb6+od10Vi7xXLOPGxsbStrGRJ8UnMufAStrdrrtnp3KHO2yo7i657tg973kT2Bf98F6pqatmxdQtrt1jXVrBLG/Jzc9ItKyx+aVvwl1ZITG9xcfF8Y0yR0zq3s2Wcpg1yvHsYY8YD4wGKiopM37594zpheXk58e6bDuLRe2PprJCjOEPNw3lZyTRMjP3ledu2UFL+OGcteA32248F99zDkO09qPwy+nOHOq8Aq0qdtQaoauL1gxV2uuPsw6Lycv3wXSirqOT6mYsZfDDcs9j6+eXm1HHH2d08HebwQ9sG8JNWSJ7eeNMc1olIJwD7/3p7+WoaJ011BtbEL08Bd7MzskUQaOYpFn/xEW9MGMxFC2fw/LHnwqJFVB1xRMznTiSzwgsDdpJNPAOgFCUe4jXuU4FL7NeXAC8HLb/Yzpo5Ctik8fbEcTM7457zD2dV6Rm0a215jR22bOK+V+7m8RdGsbl1W879v9HkjvkPtGsX17kTzQoZ0LOQOSX9WFV6BnNK+mWUYYfM71dQvEPEsIyIPI/Vebq7iKwGbgZKgckiMhD4BjjP3nw6cDqwEtgCXJYEzb4nnmySRLMzAtP2jZ6xgmGTFmCM4cxP32HkWw/TftsW/tPnIsYefT612Tm85DAFYDSplsHzxrbJyaJqS61mVjRBC4EpqSKabJmLQqw6wWFbAwxJVFQmU1VTy/Uz45+mLZYbQ6hp+wo2/8CtbzzESSvnsaDTgYw47WpWdNwHsDzzWHGaNzY3J5v7LuihRr0JgZHC8MsUgpozriQDLT/gItEY3nWbtlJT2zgaFm0udyJphqNnrKBm+w4uXDiDG2Y/Rk59HbcW/5nHivpTn2WFUcQ+Zp/SWQw/pSv5YY/Y5Ng+yU9PN4H2WLfiEwT0yUZJGmrcXSJaw2ul+jXv6ogm5pqIEd1p1Zc89/oDHPPNIj741WGUnPo3vu6wV8N64Ze0poD2O47JdjxWtNo1juzMgJ6FlG/6PGL2kKIkghp3l4jW8LbKdu7DjibmGpcRrauD++9nxuM3sF2yuf6UoUw8/GSMWDoCYZimceCa2jrWbXKuQOmkXePIiuIttHCYS0RreAt2aRN3NknMaYZLlsAxx8C111J11HGc9ddxPN/j1AbDHjhvKO3RDoLyc90URclU1Li7RLSGNz83J+5cbicjKkDxwR0bb7h9O4waBUccAV9+Cc8/z55vv8E1l/ZzPG8o7aGeMprSEvLTFcVvaFjGJaKZDq2sopJ1322m9PUF7JWfG3M2yYCehXz89QaenftNQ3zcAFPmV1K0967WsT78EAYOtLz2P/wB7r8fdt+9YX+n84XSXrBLq5i0qTFXFO+gnrtLRPJeAx2u2+vqEypnO3v5983qOdTU1jHmlYVw7bVw9NGwcSO88go8+2yDYY9Hu9frnSiKEhr13F0knPfqVrqgU3z86K8XUfr6GKj6Dq64Au68E3ZxLigWi/by8s9jOoaiKN5BjXuKcCtdMDgzpf22n7l+9mP8YeEMVu+6F8yeDUkoQBRvfXZFUdKHhmVShFtTlQU6VU9YOY83H72SCxa9yYSjz6XilbeTZtj9MDuSoiiNUeOeItxKFxzQuRVvfjSOCVNuoSp3ZwYNfpDdxv6HM485wE25DWgVQ0XxJxqWcZlQIYyEh50bA889B1dfTeeffoJ//5uDR4xgQqvoM1riQUefKoo/UePuIpFKEMQ97Pzbb+HKK61JqY88kpnX3c6/VsKaf73ZUO1x9vLvkxIT19GniuJPNCzjIq6HMOrrYdw46N7d6iy97z7KHpzM0IXbG8XAn5n7TdJi4jr6VFH8iRp3F3E1hPH559Cvn+Wx9+4NixfDNdcw+q2VzW4gTXEzJq6jTxXFn2hYxkXiCWE0jdH/44T92f/ZRzhw7Gi2Zefw4DnX0e2GqxmwX2cg+huFmzFxHX2qKP5DjbuLRFOCIJimMfqdP1vKfvdfwaHfreSNA4/ippOuZH373ch5YRGINNSBcbqBNEVj4orSstGwjIs4hTDO6VXI6Bkr2LdkGn1KZ1FV80sZ3UCMvtWOWv7+ztNMfXIYe/70A4P7lzDo9zeyvv1uANTWG0ZOXQo4x8CbojFxRVHUc3eZUFPbgdXZWbmxjrKKSgb0LGRNVQ1HVH7Kna+N4cAfv2VK92JuOeEvVOXu3Oy4gZtCqLlRY8mW0RGnipL5qHFPIk7ZM/XGWPVkDsrnzncf49z3X2Jt+9259NyRlO9fFNVxE4mBJzJVn6Io/iEh4y4iXwGbgTpghzGmSER2BSYB+wBfAecbYzYmJtOfhOrU3HfBB3DoQM7/6iue6/U7bj/uYqpbtw17rA5t3anQqPOdKkrLwI2Ye7ExpocxJuB2lgAzjTEHAjPt9y2Spp2aO2+tpt8jD/DMpJsgJwfeeYe2j4xjl4LdGmL0/3fUr8jJlkb75WQLN5/ZveF9WUUlfUpnNcTxY8lp1xGnitIySEZYpj/Q1379JFAOjEjCeTxPcPbMyZ99wC1vPkTHLVV8dtkQDho7GnJzGUDzcEjR3ruGjIknGlbREaeK0jIQY5pO/RDDziKrgI1YEwI9bIwZLyJVxpj8oG02GmM6OOw7CBgEUFBQ0GvixIlxaaiuriYvLy+ufVNBzZr1HDRmDPvPm8MPe+/L4qv+hhzRM+7jrfhus+Pcpq2ys+i6Z/uI+1fV1FK5sYb6oM89S4TCDrnNJufwets2xU96/aQV/KXXT1ohMb3FxcXzg6ImjUjUc+9jjFkjInsAb4rI8mh3NMaMB8YDFBUVmb5xlqstLy8n3n2TijHw9NNwzTXw889w223sPnw4MmdOQnovK5mGcYimCURdsybabBnPtm0I/KTXT1rBX3r9pBWSpzch426MWWP/Xy8iLwG9gXUi0skYs1ZEOgHrXdDpL775xpoR6fXXrWnvJkyAQw5x5dBuhFV0xKmiZD5xd6iKSDsRaR94DZwMLAGmApfYm10CvJyoSN9QXw9jx1qFvt59F8aMsf67ZNhBC3kpihIdiXjuBcBLIhI4znPGmNdF5CNgsogMBL4Bzktcpg9YsQIuvxzeew9OOgnGj4d99nH9NE6DmHQQkqIoTYnbuBtjvgQOd1j+I3BCIqJ8RW0t3HMPjBwJubnw+ONwySUgEnHXeNGwiqIokdARqolQUQEDB1r/zz7bCsnsuWe6VSmKomjhsLjYuhVuvBF+8xtYswZeeAGmTFHDriiKZ1DPPVbmzLG89RUrrPDLvffCrrumW5WiKEoj1HOPlupquOoqOO44y3N//XV44gk17IqieBI17tEwYwYceig8+CAMHQpLlsApp6RblaIoSkjUuIdjwwa49FI49VRo0+aX3HUfDW1WFKVlosY9FFOmQLdu8MwzcMMNsGAB9OmTblWKoihRoR2qTVm71gq9vPgi9OxpxdZ79Ei3KkVRlJhQzz2AMVYHabduMG0alJbCvHlq2BVF8SXquQN89RUMGgRvvgnHHguPPgpdtVaLoij+pWV77nV18MADVibMBx9YI0zfflsNu6Iovqfleu6ffmoV+nr/fSsbZtw42HvvdKtSFEVxhZbnudfWwm23WbH05cvhqadg+nQ17IqiZBQty3P/5BP4859h4UI47zwrJFNQkG5ViqIortMyPPeaGigpgd69Yd06K81x8mQ17IqiZCyZ77m/+64VW//sM6vg1+jR0KHZfN2KoigZReZ67j/9BEOGwPHHw/btVprjo4+qYVcUpUWQmcb9tdes9MaHHoJrrrEKfZ14YrpVKYqipIzMMu4//ggXXwynn24V95ozB+67D9q1S7cyRVGUlJI04y4ip4rIChFZKSIlyToPYJUOmDwZDjkEnn8e/vlPa+q7o49O6mkVRVG8SlI6VEUkGxgLnASsBj4SkanGmGVun6vVDz9Y85eWlUGvXvDWW/DrX7t9GkVRFF+RrGyZ3sBKY8yXACIyEegPuGvcp0+n96WXWmUE7roLhg2DnTI/AUhRFCUSYoxx/6Ai5wKnGmMut9//CTjSGDM0aJtBwCCAgoKCXhMnToz5PLmVlex7772sGjaMms6d3RGfZKqrq8nzyWQfftIK/tLrJ63gL71+0gqJ6S0uLp5vjClyWpcsN1ccljW6ixhjxgPjAYqKikzfvn3jOlF5YSHx7psOysvLfaPXT1rBX3r9pBX8pddPWiF5epPVoboa6BL0vjOwJknnUhRFUZqQLOP+EXCgiOwrIq2AC4GpSTqXoiiK0oSkhGWMMTtEZCgwA8gGHjPGLE3GuRRFUZTmJC21xBgzHZierOMriqIoocmsEaqKoigKoMZdURQlI1HjriiKkoGocVcURclAkjJCNWYRIt8DX8e5++7ADy7KSTZ+0usnreAvvX7SCv7S6yetkJjevY0xHZ1WeMK4J4KIfBxq+K0X8ZNeP2kFf+n1k1bwl14/aYXk6dWwjKIoSgaixl1RFCUDyQTjPj7dAmLET3r9pBX8pddPWsFfev2kFZKk1/cxd0VRFKU5meC5K4qiKE1Q464oipKB+Nq4p3QS7jgQka9EZLGILBCRj+1lu4rImyLyuf2/Qxr1PSYi60VkSdAyR31iMcZu60UicoQHtI4UkUq7fReIyOlB6663ta4QkVNSqdU+fxcRmS0in4rIUhG52l7uufYNo9WT7SsibUTkQxFZaOsdZS/fV0Tm2W07yS43joi0tt+vtNfv4wGtT4jIqqC27WEvd+97YIzx5R9WKeEvgP2AVsBCoFu6dTXR+BWwe5NldwEl9usS4M406jseOAJYEkkfcDrwGtYsW0cB8zygdSRwncO23ezvQ2tgX/t7kp1ivZ2AI+zX7YHPbF2ea98wWj3ZvnYb5dmvc4B5dptNBi60l48DrrRfDwbG2a8vBCZ5QOsTwLkO27v2PfCz594wCbcxZjsQmITb6/QHnrRfPwkMSJcQY8w7wIYmi0Pp6w88ZSzmAvki0ik1SkNqDUV/YKIxZpsxZhWwEuv7kjKMMWuNMZ/YrzcDnwKFeLB9w2gNRVrb126javttjv1ngH7AC/bypm0baPMXgBNExGkq0FRqDYVr3wM/G/dC4Nug96sJ/4VMBwZ4Q0Tm2xOCAxQYY9aC9aMC9kibOmdC6fNqew+1H18fCwpxeUqrHQboieW1ebp9m2gFj7aviGSLyAJgPfAm1tNDlTFmh4OmBr32+k3AbunSaowJtO1tdtveJyKtm2q1ibtt/WzcI07C7QH6GGOOAE4DhojI8ekWlABebO+HgP2BHsBa4B57uWe0ikgeMAW4xhjzU7hNHZalVLODVs+2rzGmzhjTA2t+5t7AIWE0pVVvU60icihwPXAw8BtgV2CEvblrWv1s3D0/CbcxZo39fz3wEtaXcF3gMcv+vz59Ch0Jpc9z7W2MWWf/cOqBR/glNOAJrSKSg2UsnzXGvGgv9mT7Omn1evsCGGOqgHKs+HS+iARmlwvW1KDXXr8L0Yf4XCNI66l2KMwYY7YBj5OEtvWzcff0JNwi0k5E2gdeAycDS7A0XmJvdgnwcnoUhiSUvqnAxXZv/lHApkB4IV00iUX+Hqt9wdJ6oZ0lsS9wIPBhirUJMAH41Bhzb9Aqz7VvKK1ebV8R6Sgi+fbrXOBErH6C2cC59mZN2zbQ5ucCs4zde5kmrcuDbvCC1TcQ3LbufA9S1WucjD+snuXPsOJtN6ZbTxNt+2FlFCwElgb0YcX6ZgKf2/93TaPG57Eet2uxPIaBofRhPS6Otdt6MVDkAa1P21oW2T+KTkHb32hrXQGcloa2PRbrcXoRsMD+O92L7RtGqyfbF/g1UGHrWgL8y16+H9ZNZiXwP6C1vbyN/X6lvX4/D2idZbftEuAZfsmoce17oOUHFEVRMhA/h2UURVGUEKhxVxRFyUDUuCuKomQgatwVRVEyEDXuiqIoGYgadyUjEJE6u7reEhH5n4i0TeBYfUXkVfv1WRKm4qiI5IvI4DjOMVJErotXo6JEQo27kinUGGN6GGMOBbYDfw1eaQ8Kifn7boyZaowpDbNJPlbVQUXxFGrclUzkXeAAEdlHrBrl/wU+AbqIyMki8oGIfGJ7+HnQMDfAchF5Dzg7cCARuVREHrRfF4jIS3Zt7oUicgxQCuxvPzWMtrcbLiIf2UWhRgUd60ax6p+/BXRNWWsoLRI17kpGYdcOOQ1rdB9YRvQpY0xP4GfgJuBEYxV0+xj4u4i0waqdciZwHLBniMOPAd42xhyOVVt+KVZN9i/sp4bhInIy1nD83lgFt3qJyPEi0gurREZPrJvHb1y+dEVpxE6RN1EUX5Brl1UFy3OfAOwFfG2suthgFZfqBsyxy3m3Aj7Aqs63yhjzOYCIPAMMojn9gIvBqvQHbJLmM2mdbP9V2O/zsIx9e+AlY8wW+xyeqYOkZCZq3JVMocZYZVUbsA34z8GLsOppX9Rkux64VwJWgDuMMQ83Occ1Lp5DUSKiYRmlJTEX6CMiBwCISFsROQhYDuwrIvvb210UYv+ZwJX2vtkisjOwGcsrDzAD+HNQLL9QRPYA3gF+LyK5drXQM12+NkVphBp3pcVgjPkeuBR4XkQWYRn7g40xW7HCMNPsDtWvQxziaqBYRBYD84HuxpgfscI8S0RktDHmDeA54AN7uxeA9saaxm4SVsXFKVihI0VJGloVUlEUJQNRz11RFCUDUeOuKIqSgahxVxRFyUDUuCuKomQgatwVRVEyEDXuiqIoGYgad0VRlAzk/wFBulfGSpm93wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(y_test,y_test_pred)\n",
    "plt.plot([0,350],[0,350],'r')\n",
    "plt.xlabel('Actual')\n",
    "plt.xlabel('Predicted')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Manually Computing the Coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can manually compute the regression coefficients using the least-squares matrix formula in class.  Of course, using the `sklearn` package is easier.  I illustrate the manual computation only to show what is happening \"under the hood\" of the linear regression routine. \n",
    "\n",
    "To compute the coefficients manually, we first construct the matrix $\\mathbf{A} = [ \\mathbf{1} ~ \\mathbf{X} ]$.  This is performed with the `ones` command and `hstack`.  Note that after we do this, `A` has 11 columns -- one more column that the data matrix `X_tr`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(300, 11)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ones = np.ones((ns_train,1))\n",
    "A = np.hstack((ones,X_tr))\n",
    "A.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next use the `lstsq` method to solve $\\mathbf{y} \\approx \\mathbf{A\\beta}$.  This will find the desired least-squares fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = np.linalg.lstsq(A,y_tr,rcond=None)\n",
    "beta = out[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 152.34786452,  -16.57607993, -254.66532396,  560.98630022,\n",
       "        278.91811152, -393.41357305,   97.05460405,  -19.0023093 ,\n",
       "        169.46450327,  632.95050374,  114.21638941])"
      ]
     },
     "execution_count": 152,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that `beta[0]` mathches the intercept in `regr.intercept_` from the regression fit and `beta[1:]` matches the coefficients in `regr.coef_`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -16.57607993, -254.66532396,  560.98630022,  278.91811152,\n",
       "       -393.41357305,   97.05460405,  -19.0023093 ,  169.46450327,\n",
       "        632.95050374,  114.21638941])"
      ]
     },
     "execution_count": 153,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regr.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "152.34786451820108"
      ]
     },
     "execution_count": 154,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regr.intercept_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us verify this is the same as direct matrix calculation.  For a small problem like this, it makes no difference.  But, in general, using a matrix inverse like this is *much* slower computationally than using either the `lstsq` method or the `LinearRegression` class.  So, **do not ever** solve a least squares problem like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 152.34786452  -16.57607993 -254.66532396  560.98630022  278.91811152\n",
      " -393.41357305   97.05460405  -19.0023093   169.46450327  632.95050374\n",
      "  114.21638941]\n"
     ]
    }
   ],
   "source": [
    "beta = np.linalg.inv(A.T.dot(A)).dot(A.T.dot(y_tr))\n",
    "print(beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## In-Class Simple Exercise\n",
    "\n",
    "You are given target values `y` and features `x1` and `x2` below.  Fit the model on the first 4 data points and test the model on the fifth data point.  You may want to use the following steps\n",
    "\n",
    "*  Construct the training training data `X_tr,y_tr`\n",
    "*  Create a regression object `regr = linear_model.LinearRegression()`\n",
    "*  Fit the model with the `regr.fit()` method\n",
    "*  Predict the value on the test value with the `regr.predict()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = np.array([0,1,3,5,4])\n",
    "x2 = np.array([0,0.7, 4.3, 15.1, 13.2])\n",
    "y = np.array([-2, -0.9, 1.5, 18, 13])\n",
    "\n",
    "# TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}