Linear regression sample from EFavDB/linear-regression

linear-regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpreting the results of linear regression\n",
    "Perform linear regression on simulated data to predict the weight of adult women using height and shoe size.  The results are discussed in a [blog post](http://efavdb.com/interpret-linear-regression)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# import some libraries. We use statsmodels.api.OLS for the linear regression since it contains a much \n",
    "# more detailed report on the results of the fit than sklearn.linear_model.LinearRegression\n",
    "\n",
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "import random\n",
    "from scipy.stats import t, f\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# set the population parameters for the simulated data\n",
    "\n",
    "# height in inches\n",
    "mean_height = 65\n",
    "std_height = 2.25\n",
    "\n",
    "# shoe size in inches\n",
    "mean_shoe_size = 7.5\n",
    "std_shoe_size = 1.25\n",
    "\n",
    "# correlation b/w height and shoe size\n",
    "r_height_shoe = 0.98\n",
    "# covariance b/w height and shoe size\n",
    "var_height_shoe = r_height_shoe*std_height*std_shoe_size\n",
    "\n",
    "# matrix of means\n",
    "mu = (mean_height, mean_shoe_size)\n",
    "# covariance matrix for two-regressor model\n",
    "cov = [[np.square(std_height), var_height_shoe],\n",
    "        [var_height_shoe, np.square(std_shoe_size)]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Generate the simulated data\n",
    "\n",
    "# number of data points\n",
    "n = 20\n",
    "random.seed(85)\n",
    "\n",
    "# height and shoe size\n",
    "X1 = np.random.multivariate_normal(mu, cov, n)\n",
    "# height, alone\n",
    "X0 = X1[:, 0]\n",
    "\n",
    "weight = -220 + np.random.normal(X0*5.5, 10, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA38AAAFHCAYAAAARAfVmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2cZGV95/3PT5nRBjNGsbMmgF2SgAMrA0wAXUXTgwxG\no5hoAhlvEx86umRcw01MTDDJzpDNboiraxLvzEZjB2MW2tFEoxhWAZ02sj6k5cHBAGY2WgOLBsuI\nI8ZeGOB3/1Gnx5qmH6ur6pyq+rxfr35N1anqc34z3fWd6zrnOtcVmYkkSZIkabA9quwCJEmSJEnd\nZ+dPkiRJkoaAnT9JkiRJGgJ2/iRJkiRpCNj5kyRJkqQhYOdPkiRJkoZAVzt/ETEZEfdExN6WbZsi\n4jMR8cWI+HBEPK7ltUsjYl9E3B4R53WzNknDzXyS1G8Wya2fjYgvRcRDEbG5zPokVV+3r/xdATx/\n3rZ3A2/KzFOBDwFvAoiIk4ELgJOAFwC7IiK6XJ+k4WU+Seo3C+XWrcDPAJ/qfTmS+k1XO3+ZeQNw\n77zNJxTbAa4HXlY8Ph94X2Y+mJl1YB9wVjfrkzS8zCdJ/Wah3MrML2fmPsATUpKWVcY9f/8QEecX\njy8Aji0eHwPc1fK+u4ttktQr5pMkSRpYZXT+XgO8PiJmgKOAB0qoQZIWYj5JkqSBdUSvD5iZ/0gx\nXj0iTgB+qnjpbuC4lrceW2x7hIjIbtYoqTyZWdrQJfNJ0mLKzKZOMJukwbTabOrFlb+gZRx6RIwW\nfz4K+G3gT4uXPgL8fESsj4inAj8G/P1iO83MSn/t2LGj9Bqs0Rr7qcbMUtolA51PVfq5W0v1a6la\nPVWqpWIOy60FXltU2f+O/fQzt0Zr7Ica29HVK38RcRUwDhwdEXcCO4AfiIjXAwl8MDPfA5CZt0XE\n+4HbgIPA9mz3byVJyzCfJPWbRXLrXuAdwJOAj0bELZn5gvKqlFRlXe38ZebLF3npjxd5/+8Dv9+9\niiSpyXyS1G+WyK2/6WkhkvpWGRO+DIXx8fGyS1iWNXaGNaqqqvRzt5aFVakWqFY9VapFvdEPP3Nr\n7AxrLE/048iliHDElTSAIoIcgEkVzCdpsJhNkqqonWzyyp8kSZIkDQE7f5IkSZI0BOz8SZIkSdIQ\nsPMnSZIkSUPAzp8kSZIkDQE7f5IkSZI0BOz8SZIkSdIQsPMnSZIkSUPAzp8kSZIkDQE7f5IkSZI0\nBOz8SZIkSdIQsPMnSZIkSUPAzp8kSZIkDQE7f5IkSZI0BOz8SZIkSdIQsPMnSZIkSWvUaDSYmZmh\n0WiUXcqi7PxJkiRJ0hpMTe1mbGwjW7dexNjYRqamdpdd0oIiM8uuYdUiIvuxbklLiwgyM8quYy3M\nJ2nwmE2SltJoNBgb28js7B5gE7CXkZEt7N9/B6Ojo107bjvZ5JU/SZIkSWpTvV5n/foazY4fwCbW\nrRujXq+XV9Qi7PxJkiRJUptqtRoPPFAH9hZb9nLw4H5qtVp5RS2iq52/iJiMiHsiYm/LtlMj4rMR\ncXNE/H1EnNHy2qURsS8ibo+I87pZm6ThZj5JkqROGB0dZXJyFyMjW9iwYTMjI1uYnNzV1SGf7erq\nPX8RcTbwXeC9mbmp2PZx4G2ZeW1EvAB4U2ZuiYiTgSuBM4FjgeuBExYaoO64dWkw9fK+GvNJ0kp5\nz5+klWg0GtTrdWq1Wk86fu1k0xHdKgYgM2+IiLF5mx8GHl88/kHg7uLx+cD7MvNBoB4R+4CzgM93\ns0ZJw8l8kiRJnTQ6OlrJq32tutr5W8QlwMcj4m1AAM8qth8DfLblfXcX2ySpV8wnSZI0sMqY8OWX\ngYsz8yk0G1p/XkINkrQQ80mSJA2sMq78vTIzLwbIzL+KiHcX2+8Gjmt537F8f8jVI+zcufPQ4/Hx\nccbHxzteqKTump6eZnp6uuwyWplPkqqYTZLUEV1f5D0iasDVmXlK8fwfgO2Z+amIeB5weWae2TKh\nwjNoDqe6DidUkIZKrydVMJ8krYQTvkiqospN+BIRVwHjwNERcSewA3gt8McR8Wjg/wKvA8jM2yLi\n/cBtwEGaDTBTSlJXmE+SJGnYdP3KXzd49koaTJ5dl1RFZpOkKmonm8qY8EWSJEmrFBGTEXFPROxt\n2faEiLg2Ir4cER+PiMcvtQ9Jw83OnyRJUn+4Anj+vG2/CVyfmU8DPglc2vOqJPUNO3+SJEl9IDNv\nAO6dt/klwF8Uj/8C+OmeFiWpr9j5k7SoRqPBzMwMjUaj7FIkaVlDmlk/lJn3AGTmPwM/VHI9kirM\nzp+kBU1N7WZsbCNbt17E2NhGpqZ2l12SJC3KzDrEWV0kLcrZPiU9QqPRYGxsI7Oze4BNwF5GRraw\nf/8djI6Odu24zqgnqR3dzqwqZVNEjNFcn3RT8fx2YDwz74mIJwN7MvOkBb4vd+zYcej5+Pg44+Pj\nPapaUidMT08zPT196Plll11WrXX+JPWner3O+vU1Zmc3FVs2sW7dGPV6vaudP0lqx5BlVhRfcz4C\nvAr4A+CVwIcX+8adO3d2sy5JXTb/pM1ll1226n047FPSI9RqNR54oA7MzSa+l4MH91Or1corSpIW\nMSyZFRFXAZ8BToyIOyPi1cDlwNaI+DLwvOK5JC3IK3+SHmF0dJTJyV1MTGxh3boxDh7cz+TkrkE8\ngy5pAAxLZmXmyxd56dyeFiKpb3nPn6RFNRoN6vU6tVqtJ42oKt1X0y7zSSpPtzLLbJJURe1kk50/\nSZVhA0tSFZlNkqqonWzynj9JkiRJGgJ2/iRJkiT1pUajwczMDI1Go+xS+oKdP0mSJEl9Z2pqN2Nj\nG9m69SLGxjYyNbW77JIqz3v+JFWG99VIqiKzSaqeRqPB2NhGZmf3AJuAvYyMbGH//jsGbqbfxXjP\nnyRJkqSBV6/XWb++RrPjB7CJdevGqNfr5RXVB+z8SZIkSeortVqNBx6oA3uLLXs5eHA/tVqtvKL6\ngJ0/SZIkSR3Ti0lYRkdHmZzcxcjIFjZs2MzIyBYmJ3cNzZDPdnnPn6TK8L4aSVVkNkkrNzW1m4mJ\n7axf37wyNzm5i23bLuza8RqNBvV6nVqtNnQdPxd5l9TXbGBJqiKzSVoZJ2HpLSd8kSRJklQKJ2Gp\nPjt/kiRJktbMSViqr6udv4iYjIh7ImJvy7b3RcRNxddXI+KmltcujYh9EXF7RJzXzdokDTfzSZKk\nznISlurr6j1/EXE28F3gvZm5aYHX3wp8OzN/LyJOAq4CzgSOBa4HTlhogLrj1qXB1Mv7aswnSSvl\nPX/S6gzzJCy91E42HdGtYgAy84aIGFviLRcA48XjlwDvy8wHgXpE7APOAj7fzRolDSfzSZKk7hgd\nHbXTV1Gl3fMXEc8B/jkzv1JsOga4q+UtdxfbJKmnzCdJksrXi/UCh02ZE75sA6ZKPL4kLcZ8kiSp\nRFNTuxkb28jWrRcxNraRqandZZc0ELo67HMxEfFo4KXA5pbNdwPHtTw/tti2oJ07dx56PD4+zvj4\neEdrlNR909PTTE9Pl13GYcwnSVXMJmmYNBoNJia2Mzu7h9nZ5nqBExNbOPfccw4NJ/W+wvZ0fZH3\niKgBV2fmKS3bfhL4jczc0rLtZOBK4Bk0h1NdhxMqSEOl15MqmE+SVsIJX6TempmZYevWizhw4MZD\n2zZs2Mz117+TM888k6mp3UxMbGf9+ubSEpOTu9i27cISKy5H5RZ5j4irgM8AJ0bEnRHx6uKlC5k3\npCozbwPeD9wGXANsN6UkdYv5JElSNS21XmDrVcEDB25kdnYPExPbvS9whbp+5a8bPHslDSbPrkuq\nIrNJ6r25q3vr1o1x8OD+Q1f3lrsqOEzaySY7f5IqwwaWpCoym6RyLHRfX6PRYGxsI7Oze4Dm/YAj\nI1vYv/+Oobv3r3Lr/EmSJElSOxZaL3B0dJTJyV1MTGw57KrgsHX82uWVP0mV4dl1SVVkNknV42yf\nDvuU1OdsYEmqIrNJUhVVbrZPSZIkSVI12PmTJEmSpCFg50+SJEmShoCdP0mSJEkaAnb+JEmSJGkI\n2PmTJEnqYxFxcUTcWnz9Stn1SKouO3+SJEl9KiL+LTABnAGcBrwoIo4vtypJVWXnT5IkqX+dBHw+\nM+/PzIeAvwNeWnJNkirKzp8kSVL/+hLwnIh4QkQcCbwQOK7kmiRVlJ0/SZKkPpWZdwB/AFwHXAPc\nDDxUalEaCI1Gg5mZGRqNRiX2o844ouwCJEmS1L7MvAK4AiAi/jNw10Lv27lz56HH4+PjjI+P96A6\n9aOpqd1MTGxn/foaDzxQZ3JyF9u2XVjaftQ0PT3N9PT0mvYRmdmZanooIrIf65a0tIggM6PsOtbC\nfJIGT9WzKSJGM7MREU8BPgY8MzO/M+89ZpNWpNFoMDa2kdnZPcAmYC8jI1vYv/8ORkdHe74fLa6d\nbHLYpyRJUn/764j4EvBhYPv8jp+0GvV6nfXrazQ7bACbWLdujHq9Xsp+1FkO+5QkSepjmfncsmvQ\n4KjVmkM0YS9zV+wOHtxPrVYrZT/qLK/8SZIkSQJgdHSUycldjIxsYcOGzYyMbGFycteqh2p2aj/q\nLO/5k1QZVb+vZiXMJ2nwmE0aRo1Gg3q9Tq1WW1OHrVP70SO1k012/iRVhg0sSVVkNkmqIid8kSRJ\nkiQtqKudv4iYjIh7ImLvvO1viIjbI+LWiLi8ZfulEbGveO28btYmabiZT5Ikadh0e7bPK4B3AO+d\n2xAR48CLgVMy88GIeFKx/STgAuAk4Fjg+og4wTEKGjaOje8Z80kqgRknSeXp6pW/zLwBuHfe5l8G\nLs/MB4v3fLPY/hLgfZn5YGbWgX3AWd2sT8Oj0WgwMzNDo9Eou5QlTU3tZmxsI1u3XsTY2EampnaX\nXdLAMp+k3mejGSdJ5Srjnr8TgedGxOciYk9E/Hix/Rjgrpb33V1sk9akXxobjUaDiYntzM7u4cCB\nG5md3cPExPbKd1gHjPmkodHrbDTjJKl8ZXT+jgCekJnPBN4EfKCEGjQk+qmxUa/XWb++RnMhVIBN\nrFs3Rr1eL6+o4WM+aSiUkY1mnCSVr9v3/C3kLuCDAJk5ExEPRcTRNM+kP6XlfccW2xa0c+fOQ4/H\nx8cZHx/vRq3qc3ONjdnZRzY2qnavSa1W44EH6sBemo2jvRw8uJ9arVZqXd00PT3N9PR02WW0Mp80\nFMrIxn7KuApmkyR1RNfX+YuIGnB1Zp5SPH8dcExm7oiIE4HrMnMsIk4GrgSeQXM41XXAghMquFaN\nVqrRaDA2tpHZ2T3MNTZGRrawf/8dlev8QXMY1sTEdtatG+Pgwf1MTu5i27YLyy6rZ3q9lpb5pGFV\nVjb2a8a5zp+kKqrcIu8RcRUwDhwN3APsAP6S5ix7pwH3A2/MzE8V778UmAAOAhdn5rWL7NcA04r1\nW2NjmGfC62UDy3zSsCsrG/sx4+z8SaqiynX+usUA02r1Y2NjGNnAknrLbFwZs0lSFbWTTWXc8yf1\n3OjoaMcaNjaWJEmS1I/KmO1T6lv9smyEJC3HPJOk4eOwT2mF+m3ymH7k0CqpN8yz1TGbJFVRO9nk\nlT9phVyjStKgMM8kaTjZ+ZNW6PA1qqDKa1RJ0lLMM0kaTnb+pGU0Gg1mZmYAmJzcxcjIFjZs2MzI\nyBYmJ3c5REpS3xkdHe1Jns3lZ6PR6Oh+JUnt8Z4/aQlz62CtX988Sz45uYtzzz3H2T67xPtqpN7q\n5uzFC+VnlddYXYrZJKmKXOdP6iAnROg9G1jSYBi0/DSbJFWRE75IHeSECJLUHvNTkqrJzp+0CCdE\nkKT2mJ+SVE12/qRF9GpCBHBSBEn9abHs6mV+SlqYbQstxHv+pGV0c0IEGKxJEdbK+2qk/rGS7Op2\nfvZKJ7MpIs4AngP8CDALfAm4LjPv7cT+lziu2TREbFsMByd8kfrMoE2KsFZ2/qT+MGzZ1YlsiohX\nA28AvgrcCHwDeCxwIvBsmp3A38nMO9dY7mLHN5uGxLB9PodZO9l0RLeKkbS8uUkRZmcfOSmCAS2p\nqsyuthwJPDszZxd6MSJOA04AutL502Brvcru51NL8Z4/qUROiiCpH5ldq5eZf7JYx694/ZbM/EQv\na9JgmJrazdjYRrZuvYixsY3cdNMtfj61KK/8SSWamxRhYmIL69aNcfDgfidFkFR5Zlf7ImIUeC1Q\no6UdlpmvKasm9a9Go8HExHZmZ/cUV/r2csklW3j72y/nkkv8fOqRvOdPqoBBmRRhrbznT+ovw5Jd\nHZ7w5TPAp2ne9/fQ3PbM/OtO7H+J45pNA2hmZoatWy/iwIEbD23bsGEz11//zkNDQAf98znMnPBF\nUl+z8yepijrc+bslM0/rxL5WeVyzaQB9f3KXvwaOAv6VkZGXObnLkGgnm5a85y8i/l1E/ElE7I2I\nRkTcGRHXRMTrI+LxaytXktpnPknqUx+NiBd2cocRcUlEfKnIwysjYn0n96/qGh0dZWLiF4AXAq8A\nXsjExCvs+GlRi175i4j/CXwN+DDwBQ6fkngL8GLgv2XmR3pT6mG1efZKGkArPYNlPknqpQ5f+buP\n5iWa+4GDQACZmRva3N+PADcAGzPzgYjYDfxtZr533vvMpgHksg7DrdNLPfxCZn5z3rbvAjcVX2+L\niCetskZJ6gTzSVJfyswf6MJuHw0cFREP01xS4mtdOIYqyGUdtFqLdv5aG1YR8WTgLCCBmcz85/nv\nkaReMZ8k9ZuI2JiZd0TE5oVez8yb2tlvZn4tIt5Gc33A7wHXZub1ayhVfeTwZVeaV/5c1kFLWXap\nh4j4JeA/Ap+kOTThHRHxu5n55yv43kngRcA9mbmp2LaD5hTH3yje9ubM/Fjx2qXAa4AHgYsz89rV\n/5UkDQvzSVIf+VXgdcDbFngtgXPa2WlE/CDwEmAMOAD8VUS8PDOvmv/enTt3Hno8Pj7O+Ph4O4dU\nhbjsynCZnp5menp6TftYdrbPiPgy8KzM/Jfi+dHAZzLzacvuPOJsmkOx3juvcXVfZv63ee89CbgK\nOBM4FrgeOGGhAeqOW++tYZnKW+Vb7dh180mdYs5pKVWeiTgifhZ4fma+tnj+C8AzMvM/zHuf2TSg\nGo0GN998MwCnn366GTZEOj7bZ+FfgPtant9XbFtWZt4A3LvASwsV+RLgfZn5YGbWgX00h3KpRFNT\nuxkb28jWrRcxNraRqandZZcktTKftGbmnHqhOOG01OsbIuLpbez6TuCZEfHYiAjgecDt7dSo/jOX\nXxdccCk//dPbuP76T5Zdkipuqdk+f7V4eBpwCs1Z9ZJmI2hvZr5qRQeIGAOunndm/VU0hyZ8AXhj\nZh6IiHcAn50bphAR7wauycwPLrBPz171gDNIqddWMdun+aSOMOe0Ep248hcRbweeAXyM5gLvDZqz\nFP8YzVmKx2hmzkwb+94B/DzN2UNvBn4pMw/Oe4/ZNGDML3X6yt8PFF//BPwNzYYVNBtZX22rwqZd\nwPHFAqf/zMJj31UBczNINQMFWmeQ6rZGo8HMzAyNRqPrx1JfMp/UEZ3OObNLi8nMS2jeZ/x14OeA\n/0TzPsATgHdm5nPb6fgV+74sM0/KzE2Z+cr5HT/1l5XmSJntNPWvpWb7vKwbB8zM1t/kPwOuLh7f\nDRzX8tqxxbYFedNy95U1g9TU1G4mJrazfn3z+JOTu9i27cKuHlPlaPfGZfNJndLJnDO7BkcnJlVY\nSGZ+i2a2/FnHd66BsJoccaZPtWOpYZ9X8/2z6Y+Qmeev6AARNZrDqk4pnj95bir2iLgEODMzXx4R\nJwNX0hwScQxwHU6oULq5EGqdQaqbjRmHMAy3VQz7NJ/UMZ3IObNrsFV5wpeVMpuqr50c6XU7TdXS\n6UXe37rGeoiIq4Bx4OiIuBPYAWyJiNOAh4E68O8BMvO2iHg/cBvNMevbTanybdt2Ieeee07PZsFz\nsVKtkPmkjulEzpldktaqnRzpdTtN/W/ZpR6qyLNXg8uz58PNs+vqV2bXYDOb1AvmiFaroxO+RMTV\nEfHiiFi3wGvHR8TvRsRr2ilUWszcYqUjI1vYsGEzIyNbBmaxUieC6BzzSVUzyNm1HLNtdSLiyIj4\nnYj4s+L5CRHxorLrUvmGOUfUO0vd8/dkmrNQvQz4Ft+fkrhGc4a9/y8zP9ybMh9Rm2evBtygLbjs\nRBArs4p7/swnVdKgZddyhiXbOnnlLyJ201zq4Rcz8+kRcSTwmWKW4a4xm/pHuzkybPmj9rJpRcM+\ni0kRfhiYBf4xM7/XToGdYoCpnziMY+XaCjHzSSrFMGVbhzt/X8jMMyLi5sw8vdj2xcw8tRP7X+K4\nZtMAG5YTMTpcp9f5OyQz65n52cy8peyGldRvXIenu8wnqRxmW9seiIgRihmLI+JHgfvLLUn9rNFo\nMDGxndnZPRw4cCOzs3uYmNjuUGwtaEWdP0ntO3wdHnAdHkmDwGxr207gY8BxEXEl8AngTaVWpL7m\niRithp0/qcu8gVvSIDLb2pOZ1wIvBV4FTAFnZOZ0mTWpv3kiRqux7D1/EXFxZv7Rctt6yXHrw2cQ\nbmIehL9Dt6127Lr5pEHXD7nRDzWuVYfv+fsE8LbMvKZl27sy83Wd2P8SxzWbBpiLvQ+nrkz4EhE3\nZebmedsO3aRcBgNsuHgT8/Boo/NnPmlgmX3V0eHO31eAu4BPZuZlxbZHZFmnmU2DbxhOxOhwHe38\nRcQ24OXA2cCnW176AeDhzHxeu4WulQE2PIZpNjmtaqkH80kDzeyrlg53/m4CzgL+GDgOeAWwx86f\npNVqJ5uOWOK1zwBfB54EvK1l+318f1Cx1FVzNzHPzj7yJmYbQEPNfNJAM/sGWmTmg8D2iHgVcAPw\nhHJLkjQsFu38ZeZ+YD/w73pXjnS4w29ibp799iZmmU8adGbfQPvTuQeZ+Z6IuBV4fYn1SBoiy872\nGREvjYh9EXEgIr4TEfdFxHd6UZzkbHJaivmkQWX2DZ6I2FA8/EBEPHHuC/gq8GslliZpiKxkwpf/\nDbw4M2/vTUnLc9z68PEm5uHQxoQv5pMGmtlXDZ245y8iPpqZL4qIr9Jc4L11f5mZx6+pyOWPbzZJ\nA6Zbs33+r8x89poq6zADTBpMbXT+zCdJXdfJCV/KYjZJg6fTs32+tHj4E8CTgb8B7p97PTM/2Gad\na2aASYNpFbN9mk+SeqbDs30+G7glM/81Il4BbAb+MDPv7MT+lziu2SQNmE53/q5Y4vsyM1+zmgN1\nkgEmDaZVdP7MJ0k90+HO317gVJoz+bwHeDdwQWb+RCf2v8RxzSZpwHRl2GcVGWDSYHJolaQq6vQ6\nf5m5OSL+I3B3Zk66yLukdnR6nb+5nf7xApsPAF/IzA+v5mCS1Enmk6Q+dF9EXEpzcffnRsSjgHUl\n1yRpSCy71APwWOA0YF/xtQk4FpiIiD/sYm2StBzzSVK/uZDmPcoTmfnPNDPrv5ZbkqRhsZLZPj8H\nPDszHyqeHwF8GjgbuDUzT+56lY+syaEL0gBqY7ZP80lS1zkkXVIVtZNNK7ny9wTgcS3PjwKeWDS2\n7l/4WySpJ8wnSZKkFVr2nj/gLcAtETFNc0HS5wL/JSKOAq7vYm2StBzzSZJUOY1Gg3q9Tq1WY3R0\ntOxypEOWvfKXmZPAs2iuo/Uh4OzMfHdm/mtm/vpS3xsRkxFxTzGt8fzX3hgRD0fEE1u2XRoR+yLi\n9og4b/V/HUnDxHyS1I8iYiQinlZ2HeqOqandjI1tZOvWixgb28jU1O6yS5IOWWqdv42ZeUdELDj1\ncGbetOzOI84Gvgu8NzM3tWw/lua6Nk8DfjwzvxURJwFXAWfSvPn5euCEhQaoO25dGkyrWOfPfJLU\nMx1e6uHFwFuB9Zn51Ig4DfjdzDy/E/tf4rhmUw80Gg3GxjYyO7uH5hxkexkZ2cL+/Xd4BVAd1+ml\nHn4VeB3wtgVeS+Cc5XaemTdExNgCL70d+HXgIy3bXgK8LzMfBOoRsQ84C/j8csfR4HCYhFbIfFJf\nMuME7KSZH9MAmXlLRDy1zILUOfV6nfXra8zOzp1T3MS6dWPU63U/86qERTt/mfm64s8tnTxgRJwP\n3JWZt0Yc1lE9Bvhsy/O7i20aElNTu5mY2M769TUeeKDO5OQutm27sOyyVEHmk/qRGafCwcw8MC9j\nvCQ3IGq15ucb9jJ35e/gwf3UarVS65LmLHvPX0QcGRG/HRHvKp6fEBEvaudgETECvBnY0c73a3A1\nGg0mJrYzO7uHAwduZHZ2DxMT22k0GmWXpgozn9QvzDi1+IeIeDnw6CKz3gF8puyi1Bmjo6NMTu5i\nZGQLGzZsZmRkC5OTu7zqp8pYyWyfVwA30pxUAZpnvD8AfLSN4/0oUAO+GM1TXscCN0XEWcV+n9Ly\n3mOLbQvauXPnocfj4+OMj4+3UY6qwmESw2l6eprp6em17MJ8Ul8w4/pLB7JpKW8AfovmcjRTwMeB\n/7SWHUbEicBumlcQAzge+J3M/OO1lap2bNt2Ieeee45DvFVJK1nk/QuZeUZE3JyZpxfbvpiZp67o\nABE14OrMPGWB174KbM7MeyPiZOBK4Bk0h1NdhxMqDA1vkBa0tci7+aS+YMb1t24s8h4RjwPIzO92\neL+PAv4P8IzMvKtlu9kkDZhuLfL+QDEcKouD/CgrXDw5Iq6iOZThxIi4MyJePe8tc2eoyMzbgPcD\ntwHXANtNqeHROkzicY87hcc85jm8/e2X2yjScswn9YWlhoI1Gg1mZmYcAjokIuKUiLgZ+AeaQ0Bv\njIind/AQ5wL/1Nrxk6Q5K7nydx7N4QknA9cCzwZelZnTXa9u8Zpsd/WR1cxu9853/hkXX/wm1q9/\nKg8+uN8JEYZMG1f+zCetWBVm2pxfg5PA9IcOL/XwGeC3MnNP8Xwc+C+Z+awlv3Hl+58EbszMXfO2\nm019pAp5peprJ5uW7fwVOz4aeCbNs+Cfy8xvtldiZxhg/WM1DRuHRamtEDOftAJV7GSZef2jw52/\nRwxNX81w9WX2vQ74GnByZjbmvZY7dnx/PivvR66uKuaVqmH+/ciXXXZZ5zt/EfE/gE8Bn87MO9qo\ns+NsXPUweVoTAAAd3klEQVSH1TZsZmZm2Lr1Ig4cuPHQtg0bNnP99e/kzDPP7F3hKk0bV/7MJy2r\nqp0sM69/dLjz9yHgJuAvi02vAH48M3+mA/s+n+aw9J9c4DWzqQ9UNa9UTd26528S+GHgHRHxlYj4\n64i4uK0KNVTmZrdrhhe0zm63kMPXxgHXxtEKmE9a1mqzqFfMvKH1GmAU+GDxNVps64RtNGcQVZ+q\nal5pcCzb+SvGpP9n4HeAPwPOAH65y3VpAKy2YePaOFot80krUYVO1kKTuph5wykz783MX8nMzcXX\nxZl571r3GxFH0pzs5YNrr1JlqUJeabCtZNjnJ4CjgM8CnwZuyMxv9KC2pWpy6EIfmJrazStf+Usc\nPPgw8MOsX/9N3vOedy47bt2bnIdXG8M+zSetyNw9NOvWjXHwYG8nk1ru/h0zr/o6POzzRODXaK4r\nemi95cw8pxP7X+K4ZlOfKDOv1F+6MuFLRLwd+HGa06f/L+DvgM9m5my7ha6VAVZ9h49Z/2HgOh77\n2Ndz553/aONGi2qj82c+acXK6GR5/85g6PSEL8CfAjcCD81tz8wbF/2mzhzXbOojnhTSSrSTTUcs\n94bMvKTY+Q8ArwKuAJ4MPKaNGjUk5sasz87OjVl/OevXv5V6vW6IqWPMJ63G6Ohoz/PnkVn4/ft3\nzMKh9WBm/veyi1C1lZFXGg7Ldv4i4j8Az6F5dr0O/DnN4VXSog4fs9482+2YdXWa+aSqMws1JyKe\nWDy8OiK2Ax+iOWoBgMz8VimFSRoqy3b+gMcC/43mgqEPdrkeDYi5iQwmJrYcNmbds1jqMPNJlWYW\nqsWNQNJckxTg11teS+D4nlckaeisaJH3qnHcev9wzLpWo5P31ZTFfNJCzML+ZjZJqqKuTPhSRQaY\nNJhsYEmqog5P+PJzwMcy876I+G1gM/CfMvPmTux/ieOaTdKA6dYi75IkSeqM3yk6fmfTXJdvkubs\nn5LUdXb+JEmSemdueYefAt6VmX8LrC+xHklDxM6fJElS79wdEe8ELgSuiYjHYHtMUo8YNhXRaDSY\nmZmh0WiUXYoktc0sk5Z1AfBx4PmZ+W3giRw+86ckdY2dvwqYmtrN2NhGtm69iLGxjUxN7S67JEla\nNbNMWl5mfi8zP5iZ+4rnX8/Ma8uuS9JwcLbPkjUaDcbGNjI7u4e5BYBHRrawf/8dTgeuoeNsn/3L\nLNMgM5skVZGzffaher3O+vU1mo0lgE2sWzdGvV4vryhJWiWzTJKk6rPzV7JarcYDD9SBvcWWvRw8\nuJ9arVZeUZK0SmaZJEnVZ+evZKOjo0xO7mJkZAsbNmxmZGQLk5O7HCYlqa+YZZL6hRNTaZh5z19F\nNBoN6vU6tVrNxpKGlvfV9D+zTIPIbBocU1O7mZjYzvr1zdEKk5O72LbtwrLLktrSTjbZ+ZNUGTaw\nJFWR2TQYnJhKg6ZyE75ExGRE3BMRe1u2/W5EfDEibo6Ij0XEk1teuzQi9kXE7RFxXjdrkzTczCdJ\nGi5OTCV1/56/K4Dnz9v2lsw8NTNPB/4W2AEQESfTXPj0JOAFwK6I6OuzbJIqzXySpCHixFRSlzt/\nmXkDcO+8bd9teXoU8HDx+HzgfZn5YGbWgX3AWd2sT9LwMp8kabg4MZUER5Rx0Ij4PeAXgW8DW4rN\nxwCfbXnb3cU2SeoZ80mSBte2bRdy7rnnODGVhlYpSz1k5m9n5lOAK4E3lFGDJC3EfJKkwTY6OsqZ\nZ55px09DqZQrfy2uonlfzU6aZ9KPa3nt2GLbgnbu3Hno8fj4OOPj492oT1IXTU9PMz09XXYZizGf\npCFV8WySpLZ1famHiKgBV2fmKcXzH8vM/108fgPwnMy8oJhQ4UrgGTSHU10HnLDQvMTDMl2x62Vp\n2PR6OnXzqTzmm/qJSz1IqqIqLvVwFfAZ4MSIuDMiXg1cHhG3RsQtwLnAxQCZeRvwfuA24Bpg+zCn\n1NTUbsbGNrJ160WMjW1kamp32SVJA8V8Ko/5JklSOVzkvYJchFTDyrPrg898Uz8ymyRVUeWu/Kk9\nLkIqaVCZb5IklcfOXwW5CKmkQWW+SeqmRqPBzMwMjUaj7FKkSrLzV0EuQippUJlvkrrF+4ml5XnP\nX4U5G56GjffVDA/zTf3EbKo+7yfWMGonm8pe509LGB0dHcrAslEoDb5hzbeVMAPVjoh4PPBu4OnA\nw8BrMvPz5VbVO3P3E8/OPvJ+Yj9H0vc57FOV4pANScPMDNQa/BFwTWaeBJwK3F5yPT3l/cTSyjjs\nU5XhkA05tErDzAysrqpnU0RsAG7OzB9d4j0Dn01TU7uZmNjOunVjHDy4n8nJXWzbdmHZZUld47BP\n9TWHbEgaZmag1uCpwDcj4gqaV/2+AFycmbPlltVb27ZdyLnnnuOwaWkJDvtUZThkQ9IwMwO1BkcA\nm4E/yczNwPeA3yy3pHKMjo5y5pln2vGTFuGVP1XG3BTwExNbDhuyYYBLGgZmoNbg/wB3ZeYXiud/\nBfzG/Dft3Lnz0OPx8XHGx8d7UZukDpmenmZ6enpN+/CeP1WOM90Nr6rfV7MS5pPWygysnn7Ipoj4\nFPDazPzHiNgBHJmZv9HyutkkDZh2ssnOn6TK6IcG1nLMJ2nw9EM2RcSpNJd6WAd8BXh1Zh5oed1s\nkgaMnT9Jfa0fGljLMZ+kwWM2VY9XyKX2sskJXyRJktQ3XA9Tap9X/iRVhmfXJVWR2VQdrocpfZ9X\n/iRJkjSw5tbDbHb8oHU9TEnLs/MnSZKkjmk0GszMzNBoNDq+b9fDlNbGzp8kSZI6otv3482thzky\nsoUNGzYzMrLF9TClVfCeP0mV4X01kqrIbFqZXt6P52yfUnvZdES3ipEkSdLwmLsfb3b2kffjdbqD\nNjo6aqdPaoPDPiVJkrRm3o8nVZ+dP0mSJK2Z9+NJ1dfVe/4iYhJ4EXBPZm4qtr0FeDFwP/BPwKsz\n8zvFa5cCrwEeBC7OzGsX2a/31EgDqJf31ZhPklbKe/5Wx/vxpN5oJ5u63fk7G/gu8N6WxtW5wCcz\n8+GIuBzIzLw0Ik4GrgTOBI4FrgdOWCipbFxJg6nHnT/zSdKK2PmTVEWVW+Q9M28A7p237frMfLh4\n+jmaDSmA84H3ZeaDmVkH9gFndbM+ScPLfJIkScOm7Hv+XgNcUzw+Brir5bW7i22SVAbzSZIkDZTS\nOn8R8VvAwcycKqsGSVqI+SRJkgZRKev8RcSrgBcC57Rsvhs4ruX5scW2Be3cufPQ4/HxccbHxztZ\noqQemJ6eZnp6uuwyDmM+SapiNklSJ3R1wheAiKgBV2fmKcXznwTeBjw3M/+l5X1zEyo8g+Zwqutw\nQgVpqPR6UgXzSdJKOOGLpCpqJ5u6euUvIq4CxoGjI+JOYAfwZmA9cF1EAHwuM7dn5m0R8X7gNuAg\nsN2UktQt5pMkSRo2Xb/y1w2evZIGk2fXJVWR2SSpiiq31IMkSZIkqRrs/EmSJEnSELDzJ0mSJElD\nwM6fJEmSJA0BO3+SJEmSNATs/EmSJKkrGo0GMzMzNBqNskuRhJ0/SZIkdcHU1G7GxjaydetFjI1t\nZGpqd9klSUPPdf4kVYZraUmqIrNp9RqNBmNjG5md3QNsAvYyMrKF/fvvYHR0tGd1SIPMdf4kSZJU\nunq9zvr1NZodP4BNrFs3Rr1eL68oSXb+JEmS1Fm1Wo0HHqgDe4stezl4cD+1Wq28oiTZ+ZMkSVJn\njY6OMjm5i5GRLWzYsJmRkS1MTu5yyKdUMu/563ONRoN6vU6tVjNQ1fe8r0YLMedUNrOpfX5+pe7x\nnr8h4yxakgadOSf1t9HRUc4880w7flJFeOWvTzmLlgaRZ9fVypxTVfRDNkVEHTgAPAwczMyz5r1u\nNkkDxit/XValhUqdRUvSWlQpzxZjzkmr8jAwnpmnz+/4SdIcO38rVLWhR86iJaldVcuzxZhz0qoE\ntuskLcNhnytQ1aFHU1O7mZjYzrp1Yxw8uJ/JyV1s23ZhafVIa9UPQ6uWU/WhVVXNs8WYc6qCfsim\niPgK8G3gIeBdmfln816vdDZJWr12sumIbhUzSOaGHs3OPnLoUZmNpW3bLuTcc89xFi1JK1bVPFuM\nOSet2LMz8+sRMQpcFxG3Z+YNZRclqVrs/K3A4UOPmmfKqzL0aHR01MaQpBWrcp4txpyTlpeZXy/+\nbETEh4CzgMM6fzt37jz0eHx8nPHx8R5WKGmtpqenmZ6eXtM+HPa5Qg49krqvH4ZWLacfhlaZZ9Lq\nVD2bIuJI4FGZ+d2IOAq4FrgsM69teU/ls0nS6rSTTXb+VsGFSqXuqnoDayX6pYFlnkkrV/Vsioin\nAh8Ckuaorisz8/J57+mLbJK0cnb+JPW1qjewVsJ8kgaP2SSpiiq3zl9ETEbEPRGxt2Xbz0bElyLi\noYjYPO/9l0bEvoi4PSLO62Ztkoab+SRJkoZNt9eDuQJ4/rxttwI/A3yqdWNEnARcAJwEvADYFRF9\nfZZNUqWZT5Ikaah0tfNXTDF877xtX87MfTQXI231EuB9mflgZtaBfTRnqpKkjjOfJEnSsOn2lb/V\nOAa4q+X53cU2SSqb+SRJkvpelTp/kiRJkqQuqdIi73cDx7U8P7bYtiAXKpX6XycWK+0R80kaIn2U\nTZK0Kl1f6iEiasDVmXnKvO17gF/LzBuL5ycDVwLPoDmc6jrghIXmJXa6Ymkw9Xo6dfNJ0kq41IOk\nKmonm7p65S8irgLGgaMj4k5gB80JFt4BPAn4aETckpkvyMzbIuL9wG3AQWC7KSWpW8wnSZI0bFzk\nXVJleHZdUhWZTZKqqHKLvEuSJEmSqsHOnyRJkiQNATt/kiRJkjQE7Pwto9FoMDMzQ6PRKLsUSaoc\nM1KSpP5h528JU1O7GRvbyNatFzE2tpGpqd1llyRJlWFGSpLUX5ztcxGNRoOxsY3Mzu4BNgF7GRnZ\nwv79dzA6OtrVY0vDyhn1+ocZqWFiNkmqImf77KB6vc769TWajRqATaxbN0a9Xi+vKEmqCDNSkqT+\nY+dvEbVajQceqAN7iy17OXhwP7VarbyiJKkizEhJkvqPnb9FjI6OMjm5i5GRLWzYsJmRkS1MTu5y\nOJMkYUZKktSPvOdvGY1Gg3q9Tq1Ws1EjdZn31fQfM1LDwGySVEXtZJOdP0mVYQNLUhWZTZKqyAlf\nJEmSJEkLsvMnSZIkSUPAzp8kSZIkDQE7f5IkSZI0BOz8SZIkSdIQsPMnSZIkSUPAzp8kSZIkDQE7\nf5IkSZI0BOz8SZIkSdIQsPMnSZIkSUPAzp8kSZIkDYGudv4iYjIi7omIvS3bnhAR10bElyPi4xHx\n+JbXLo2IfRFxe0Sc183aJA0380nSIImIR0XETRHxkbJrkVRd3b7ydwXw/HnbfhO4PjOfBnwSuBQg\nIk4GLgBOAl4A7IqI6HJ9XTM9PV12Ccuyxs6wxr418PlUpZ+7tSysSrVAteqpUi194mLgtrKLWIt+\n+JlbY2dYY3m62vnLzBuAe+dtfgnwF8XjvwB+unh8PvC+zHwwM+vAPuCsbtbXTf3wC2ONnWGN/WkY\n8qlKP3drWViVaoFq1VOlWqouIo4FXgi8u+xa1qIffubW2BnWWJ4y7vn7ocy8ByAz/xn4oWL7McBd\nLe+7u9gmSb1iPknqR28Hfh3IsguRVG1VmPDFoJJUVeaTpEqLiJ8C7snMW4AoviRpQZHZ3bZNRIwB\nV2fmpuL57cB4Zt4TEU8G9mTmSRHxm0Bm5h8U7/sYsCMzP7/APm2QSQMqM3vWcDGfJK1UL7NpNSLi\nvwCvAB4ERoAfAD6Ymb84731mkzSAVptNR3SrkBbzz0J9BHgV8AfAK4EPt2y/MiLeTnM41Y8Bf7/Q\nDqsawJL6jvkkqa9l5puBNwNExE8Ab5zf8SveZzZJ6m7nLyKuAsaBoyPiTmAHcDnwgYh4DbCf5gx6\nZOZtEfF+mjNVHQS2Z7cvS0oaWuaTJEkaNl0f9ilJkiRJKl8VJnxZVkQ8PiI+UCyu/A8R8YylFmOu\nUI1vKZ7fEhF/HREbqlZjy2tvjIiHI+KJVawxIt5QbLs1Ii6vWo0RcWpEfDYibo6Iv4+IM0qs78Si\njpuKPw9ExK9U6TOzRI2V+sysVpUWWY6IekR8ce53suRaFs2eHtex4O9dGbUU9VwSEV+KiL0RcWVE\nrC+xlouLfL211/8mETEZEfdExN6WbaXl1SL1/Gzxs3ooIjb3qpa1sO3UvRpbXrPttIYabTt1rMbV\nfWYys/JfwHuAVxePjwAeT/OenDcV234DuLyCNZ4LPKrYdjnw+xWrcUPx+FjgY8BXgSdWrUaaQ/Ou\nBY4otj+pYjU+Hvg4cF6x7QU0JwoprcaWWh8FfA04rmqfmUVqrNRnpo2/yyXA/wA+UoFavgI8oew6\niloWzJ6Sazr0e1fS8X+k+BmtL57vBn6xpFr+LbAXeAzw6CJvj+/h8c8GTgP2tmwrLa8WqedpwAnA\nJ4HNZfyc2vh7LPR/VaX+H1ikxkr9P7BYfmHbqRM/a9tOnalxVZ+Zyl/5K3qvz8nMKwCyucjyARZf\njLnnFqsxM6/PzIeLt32OZlBUqcbvFC/PrQ9UqiVq/GWaH7YHi+3frFiNB4CHaQYZwA/SXAeuCs4F\n/ikz76JCn5l5DtVYpc/MakX1FlkOKjC6Y5nsKVPrZ6MsjwaOiogjgCNp/kdehpOAz2fm/Zn5EPB3\nwEt7dfDMvAG4d97m0vJqoXoy88uZuY8+WUbBtlNXa7Tt1JkabTutTdttp9IbBivwVOCbEXFFcZnz\nXRFxJPBvcuHFmKtS48i897wG+J8l1DZnwX/HiDgfuCszby2xtjmL/axPBJ4bEZ+LiD1lDgtYpMYR\nmld83hrNiUPeAlxaYo2tLgSuKh5X6TPT6kJgaoHtZX9mVqtqiywncF1EzETEa0usYyX5WIbFfu96\nIjO/BrwNuJNmg+fbmXl9SeV8CXhOMbzpSJonMY4rqZY5P1TRvOoXtp06w7ZT92q07bQ2bbed+qHz\ndwSwGfiTzNwM/CvwmzyygVVmg2t+jd+j5Rc4In4LOJiZVy3y/b2w0L/jTprTQ+9oeV+ZZzUX+1kf\nQXP42jOBNwHvL6/EBWu8lOYZtosz8yk0w+zPyyuxKSLWAecDHyg2VekzAyxY49z2KnxmViyqucjy\ns4vf0RcCr4+Is0uqY6F8/M2SagEW/73rcQ0/SPOM8hjNIaCPi4iXl1FLZt5Bc2jTdcA1wM3AQ2XU\nsoTS86rP2HbqDNtOnWHbqYPW2nbqh87f/6F5duULxfO/pvkLdE9E/BuAaC7G/I2S6oNH1vhXwOkA\nEfEqmo2vUv5Tb7HYv2MN+GJEfJXmZeIbI6KssxqL1XgX8EGAzJwBHo6Io8spcdEafzEz/6ao8a+A\ns0qqr9ULgBtbhnpU6TMzZ67GxtyGCn1mVuPZwPkR8RWaZ+K2RMR7yywoM79e/NkAPkR5v5ML5WPZ\nE2Y84veuBOcCX8nMbxVDLT8IPKusYjLzisw8IzPHgW8D/1hWLYUq5lU/se3UGbadOsO2U2etqe1U\n+c5fcan1rog4sdj0POAf+P5izHD4Ysw9t0iNt0XET9IcBnZ+Zt5fVn2waI03ZuaTM/P4zHwqzQ/n\n6ZlZyi/2Ej/rvwHOgeZMR8C6zPyXitX4tWgurktEPI/yG04A2zh8SEBlPjMtDquxSp+Z1cjMN2fm\nUzLzeODngU/mAoss90oxLOlxxeOjgPNoDu3rucXysYxaWsz/bJThTuCZEfHYiAia/y63l1VMRIwW\nfz4F+Bm+P+SpZyVw+NWTsvNqqSv4VbiyvyTbTp1h26nrNdp2as+a2k59sc5fRJxKcxKFdTRnR3s1\nzRvl30/zvoT9wAWZ+e2K1fgFYD0w92H7XGZuL6fChWssbride/0rwBmZ+a2SSlzs3/F7NIcCnAbc\nD7wxMz9VsRqfDvwRzd/L/0tzEfCbS6zxSJqfi+Mz875i2xOp1mdmoRr3UaHPTDuK/8jemJnnl1jD\nU2le7Uuaw22uzMzSpvleLnt6XMsjfu/KEhE7aJ4sOEhzqOUvZebBkmr5O+CJRS2XZOZ0D499Fc2Z\nCY8G7qE5nO5vaA5p6nleLVLPvcA7gCfRvDJ6S2a+oBf1tMu2U2fYdupqjbadOlPjqtpOfdH5kyRJ\nkiStTeWHfUqSJEmS1s7OnyRJkiQNATt/kiRJkjQE7PxJkiRJ0hCw8ydJkiRJQ8DOnyRJkiQNATt/\naktEjEXErav8nn8fEa9Y5j2vjIh3LPLapct87ydaFra+YQX1TEXEjy73PkmSpLWy7aQqsPOntVjV\nIpGZ+c7M/B9r2O+bF/uGiHghzUV3v1sc6+wVHOe/A7+xgvdJ6hMR8dViUd5eHOujEbGhA/t5e0Sc\nXTx+V0RsXOb9r4+IV6/1uJJKYdtJpbLzp7U4omiofCkiPhYRjwGIiOMj4n9GxExEfCoiTiy274iI\nXy0enxkRX4yImyLiLfPOhB1TfP+XI+Ly4v2/D4wU7//LBWr5f4APzz2JiPuKP38iIvZExAci4vZ5\n3/tp4NyI8HMgDY5VNazWdKDMF2Xmd9ayj6Kj+ozMvKHY5+sy845lvu3PgTes5biSSmPbSaXyB6e1\nOAF4R2Y+HTgAvKzY/i7gP2TmmcCv0zxLNN+fA6/NzM3AQxzeYDsV+DlgE/DzEXFMZl4KfC8zN2fm\nLyywv2cDN7Y8b93facCvACcDPxoRzwLIzAT2FceT1Eci4sjiytvNEbE3In5u7iXgVyLixqKRNNeA\nekJEfKjY9pmIOKVlP5MR8bnie168wLGeXDTGbiqO9exi+1cj4onFsKybi9e/EhGfKF4/rzjWFyJi\nd0QcucBf5WXAx1qOtSciNheP74uI34uIW4r9jAJk5izw1Yg4o2P/oJJ6xbaTSmXnT2vxlcycO+t0\nI1CLiKOAZwEfiIibgXcC/6b1myLi8cDjMvPvi01XzdvvJzLzu5l5P3AbMLaCWp6Qmf+6yGt/n5lf\nLwLrFqDW8loD+JEV7F9StfwkcHdmnp6Zm2jpQAHfyMwfB/4U+LVi22XATZl5KvBbwHuL7b9FM3Oe\nCZwDvDUiRuYd6+XAx4oG16k0cwSKhlIxLOt04CzgLuBtEXF0se/nZeYZNDPyjQv8PeY3vlodBXwm\nM0+jebb9tS2v3Qg8Z5Hvk1Rdtp1UqiPKLkB97f6Wxw8Bj6V5QuHeopG0lFjFfud+T5f6ngfb2B80\na55d4nslVdOtNDtqvw/87dywycKHij9vBH6meHw28FKAzNxTXLF7HHAe8OKI+PXifeuBpwBfbtnf\nDDAZEeuAD2fmF4vt8zPpj4FPZuY1EfFTNM+Y/6+ICGAd8NkF/h4/TLMhtZD7M/Oalr/LuS2vfQN4\n2iLfJ6m6bDupVF7501o8IlAy8z6aw5F+9tCbIjbNe88B4DsRcWax6edXeLwHIuLRi7z25Yg4fqna\nFnEi8KUVvldSRWTmPmAzzU7g70XEb7e8PNdomd9gWUgALyuuIJ6emU/NzNaOH5n5aeC5wN3Ae2KB\nmfci4lXAcZl5Wct+ry2GW52emU/PzNfO/z6aDajHLlLbwZbHNr6kwWDbSaWy86e1WGxihVcAE8V9\nKl8Czl/gPb8EvDsibgKOpDnufbljvAu4dZGblv8W2LKC2g5tj4gfojkW/huLvFdSRUXEDwOzmXkV\n8F9pdgSX8mma2UREjAPfLGa4+zjN+1rm9nvaAsd6Cs2hpJPAu+cfKyJ+nOaQztZO4eeAZ0cxJXpx\nb+EJC9R1O/Bji/01l/j72PiS+pNtJ5UqmkN5pd6KiKPmxplHxG8AT87MS9awvycDf5GZz1/F9/y/\nwIHMvKLd40oqR0ScR7PT9zDwAHBRZt4cEV8BzsjMbxWdsv+amedExBNoTpZwPPCvwOsy80sR8Vjg\nD2nebxPAVzPz/HnH+kWaEzAcBO4DfiEz75w7FvBWmsNH5xpDX8jM1xWdzLcAj6HZePrtzPzovH2f\nDfz7uckYIuKTwK9l5k0R8Z3M3FBsfxnwU5n5muL5jcC5mXlvJ/49JVWfbSd1gp0/lSIiLgAupTmM\nqQ68KjP/ZY37/FmakzJ8d4XvfyXwl5n58FqOK0lrERF/B6x42Yji6uQlmfnK7lYmqUpsO6kT7PxJ\nklSi4h6e2cxc0TDOiHgesC8z7+xuZZKkQWPnT5IkSZKGgBO+SJIkSdIQsPMnSZIkSUPAzp8kSZIk\nDQE7f5IkSZI0BOz8SZIkSdIQ+P8B85b17FU08q8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fd056b81518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot simulated data\n",
    "plt.rcParams['figure.figsize'] = (15.0, 5.0)\n",
    "plt.subplot(1,3,1)\n",
    "# weight v. height\n",
    "plt.scatter(X1[:,0], weight)\n",
    "plt.xlabel('height (in)')\n",
    "plt.ylabel('weight (lb)')\n",
    "\n",
    "plt.subplot(1,3,2)\n",
    "# weight v. shoe size\n",
    "plt.scatter(X1[:,1], weight)\n",
    "plt.xlabel('shoe size (in)')\n",
    "plt.ylabel('weight (lb)')\n",
    "\n",
    "plt.subplot(1,3,3)\n",
    "# height v. shoe size\n",
    "plt.scatter(X1[:,0], X1[:,1])\n",
    "plt.xlabel('height (in)')\n",
    "plt.ylabel('shoe size (in)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Fit the linear models\n",
    "# add column of ones for intercept\n",
    "X0 = sm.add_constant(X0)\n",
    "X1 = sm.add_constant(X1)\n",
    "\n",
    "sm0 = sm.OLS(weight, X0).fit()\n",
    "sm1 = sm.OLS(weight, X1).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "linear model: weight ~ height\n",
      "\n",
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.788\n",
      "Model:                            OLS   Adj. R-squared:                  0.776\n",
      "Method:                 Least Squares   F-statistic:                     66.87\n",
      "Date:                Wed, 29 Jun 2016   Prob (F-statistic):           1.79e-07\n",
      "Time:                        14:28:08   Log-Likelihood:                -70.020\n",
      "No. Observations:                  20   AIC:                             144.0\n",
      "Df Residuals:                      18   BIC:                             146.0\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
      "------------------------------------------------------------------------------\n",
      "const       -265.2764     49.801     -5.327      0.000      -369.905  -160.648\n",
      "x1             6.1857      0.756      8.178      0.000         4.596     7.775\n",
      "==============================================================================\n",
      "Omnibus:                        0.006   Durbin-Watson:                   2.351\n",
      "Prob(Omnibus):                  0.997   Jarque-Bera (JB):                0.126\n",
      "Skew:                           0.002   Prob(JB):                        0.939\n",
      "Kurtosis:                       2.610   Cond. No.                     1.73e+03\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[2] The condition number is large, 1.73e+03. This might indicate that there are\n",
      "strong multicollinearity or other numerical problems. \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# print summaries of fit\n",
    "print('linear model: weight ~ height\\n')\n",
    "print(sm0.summary(), '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "linear model: weight ~ height + shoe_size\n",
      "\n",
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.789\n",
      "Model:                            OLS   Adj. R-squared:                  0.765\n",
      "Method:                 Least Squares   F-statistic:                     31.86\n",
      "Date:                Wed, 29 Jun 2016   Prob (F-statistic):           1.78e-06\n",
      "Time:                        14:28:08   Log-Likelihood:                -69.951\n",
      "No. Observations:                  20   AIC:                             145.9\n",
      "Df Residuals:                      17   BIC:                             148.9\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
      "------------------------------------------------------------------------------\n",
      "const       -333.1599    204.601     -1.628      0.122      -764.829    98.510\n",
      "x1             7.4944      3.898      1.923      0.071        -0.729    15.718\n",
      "x2            -2.3090      6.739     -0.343      0.736       -16.527    11.909\n",
      "==============================================================================\n",
      "Omnibus:                        0.015   Durbin-Watson:                   2.342\n",
      "Prob(Omnibus):                  0.993   Jarque-Bera (JB):                0.147\n",
      "Skew:                           0.049   Prob(JB):                        0.929\n",
      "Kurtosis:                       2.592   Cond. No.                     7.00e+03\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[2] The condition number is large,  7e+03. This might indicate that there are\n",
      "strong multicollinearity or other numerical problems. \n",
      "\n"
     ]
    }
   ],
   "source": [
    "print('linear model: weight ~ height + shoe_size\\n')\n",
    "print(sm1.summary(), '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta_hat: [-333.15990097    7.49444671   -2.30898743] \n",
      "\n",
      "degrees of freedom of residuals: 17 \n",
      "\n",
      "sigma_hat: 8.66991550428 \n",
      "\n",
      "standard error of beta_hat: [ 204.60056111    3.89776076    6.73900599] \n",
      "\n",
      "confidence intervals:\n",
      " [[ -7.64829352e+02   9.85095501e+01]\n",
      " [ -7.29109662e-01   1.57180031e+01]\n",
      " [ -1.65270473e+01   1.19090724e+01]] \n",
      "\n",
      "t-statistics: [-1.62834305  1.92275698 -0.34263027] \n",
      "\n",
      "p-values of t-statistics: [ 0.1218417   0.07142839  0.73607656]\n",
      "f-statistic: 31.8556171105 \n",
      "\n",
      "p-value of f-statistic: 1.77777555162e-06 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Calculate values manually for the two-regressor model: weight ~ height + shoe_size\n",
    "# OLS solution, eqn of form ax=b => (X'X)*beta_hat = X'*y\n",
    "beta_hat = np.linalg.solve(np.dot(X1.T, X1), np.dot(X1.T, weight))\n",
    "print('beta_hat:', beta_hat, '\\n')\n",
    "\n",
    "# residuals\n",
    "epsilon = weight - np.dot(X1, beta_hat)\n",
    "\n",
    "# degrees of freedom of residuals\n",
    "dof = X1.shape[0] - X1.shape[1]\n",
    "print('degrees of freedom of residuals:', dof, '\\n')\n",
    "\n",
    "# best estimator of sigma\n",
    "sigma_hat = np.sqrt(np.dot(epsilon, epsilon) / dof)\n",
    "print('sigma_hat:', sigma_hat, '\\n')\n",
    "\n",
    "# standard error of beta_hat\n",
    "s = sigma_hat * np.sqrt(np.diag(np.linalg.inv(np.dot(X1.T, X1)), 0))\n",
    "print('standard error of beta_hat:', s, '\\n')\n",
    "\n",
    "# 95% confidence intervals\n",
    "# +/-t_{1-alpha/2, n-K} = t.interval(1-alpha, dof)\n",
    "conf_intervals = beta_hat.reshape(3,1) + s.reshape(3,1) * np.array(t.interval(0.95, dof))\n",
    "print('confidence intervals:\\n', conf_intervals, '\\n')\n",
    "\n",
    "# t-statistics under null hypothesis\n",
    "t_stat = beta_hat / s\n",
    "print('t-statistics:', t_stat, '\\n')\n",
    "\n",
    "# p-values\n",
    "# survival function sf=1-CDF\n",
    "p_values = t.sf(abs(t_stat), dof)*2\n",
    "print('p-values of t-statistics:', p_values)\n",
    "\n",
    "# SSR (regression sum of squares)\n",
    "y_hat = np.dot(X1, beta_hat)\n",
    "y_mu = np.mean(weight)\n",
    "mean_SSR = np.dot((y_hat - y_mu).T, (y_hat - y_mu))/(len(beta_hat)-1)\n",
    "\n",
    "# f-statistic\n",
    "f_stat = mean_SSR / np.square(sigma_hat)\n",
    "print('f-statistic:', f_stat, '\\n')\n",
    "\n",
    "# p-value of f-statistic\n",
    "p_values_f_stat = f.sf(abs(f_stat), dfn=(len(beta_hat)-1), dfd=dof)\n",
    "print('p-value of f-statistic:', p_values_f_stat, '\\n')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}