{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Categorical Ordinary Least-Squares \n",
    "====\n",
    "\n",
    "## Unit 12, Lecture 3\n",
    "\n",
    "*Numerical Methods and Statistics*\n",
    "\n",
    "----\n",
    "\n",
    "#### Prof. Andrew White, April 19 2018"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Goals:\n",
    "---\n",
    "\n",
    "1. Regress on categorical data\n",
    "2. Learn that unordered categories have to be split into dummy variables\n",
    "3. Know when to test interaction variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import random\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from math import sqrt, pi, erf\n",
    "import seaborn\n",
    "seaborn.set_context(\"talk\")\n",
    "#seaborn.set_style(\"white\")\n",
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Regression with discrete domains\n",
    "===\n",
    "\n",
    "Let's say we have a potential drug molecule that improves your test-taking abilities. We give 10 people the drug and 7 are given a plcebo. Here's their results on an exam:\n",
    "\n",
    "Drug Group | Control Group\n",
    "---------- | ---------:|\n",
    "35 | 27\n",
    "38 | 29\n",
    "25 | 34\n",
    "34 | 32\n",
    "42 | 35\n",
    "41 | 22\n",
    "27 | 19\n",
    "32 |\n",
    "43 |\n",
    "36 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Did the drug make a difference?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We know how to solve this from hypothesis testing - Wilcoxon Sum of Ranks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RanksumsResult(statistic=2.04939015319192, pvalue=0.04042397933690852)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import scipy.stats as boop\n",
    "\n",
    "drug = np.array([35, 38, 25 , 34, 42, 41 , 27, 32, 43, 36])\n",
    "control = np.array([27, 29, 34 , 32, 35 , 22, 19])\n",
    "\n",
    "boop.ranksums(drug, control)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's try regressing it! To regress, we need two dimensions. Right now we only have one. Let's create a *category* or *dummy* variable indicating if the exam score came from the drug or control group:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1. 35.]\n",
      " [ 1. 38.]\n",
      " [ 1. 25.]\n",
      " [ 1. 34.]\n",
      " [ 1. 42.]\n",
      " [ 1. 41.]\n",
      " [ 1. 27.]\n",
      " [ 1. 32.]\n",
      " [ 1. 43.]\n",
      " [ 1. 36.]\n",
      " [ 0. 27.]\n",
      " [ 0. 29.]\n",
      " [ 0. 34.]\n",
      " [ 0. 32.]\n",
      " [ 0. 35.]\n",
      " [ 0. 22.]\n",
      " [ 0. 19.]]\n"
     ]
    }
   ],
   "source": [
    "drugged = np.concatenate( (np.ones(10), np.zeros(7)) )\n",
    "exam = np.concatenate( (drug, control) )\n",
    "\n",
    "print(np.column_stack( (drugged, exam) ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEiCAYAAAAiQw8CAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XuYHGWZ9/HvbyaTk4QJh4ABloWgYXdNXhU0HhNRCC9kAUF2WVAXYYnxvEpA9HVliairIgZBZdcYFUVXhEURJIJBhAQRwQNKUCFLAigGiUCGQCbJZOZ+/3hqoOnp6Zmema7umvl9rquvzlQ9VXV3TafueaqegyICMzOzPLQ0OgAzMxs7nHTMzCw3TjpmZpYbJx0zM8uNk46ZmeXGScfMzHLjpGNmZrlx0jEzs9w46ZiZWW6cdMzMLDdOOmZmlhsnHTMzy42TjpmZ5cZJx8zMcjOu0QE0iqQdpKT7RKNjMTMrkJ2BnogYUv7QWJ1PR1IPoPb29kaHYmZWGB0dHQAREUO6UzZmazrAE+3t7e2bNm1qdBxmZoUxdepUOjo6hnyHyM90zMwsN046ZmaWGycdMzPLzVh+pmNmVmgdW7pYsWYDGzdvY9qUCSyYNZ32yW2NDqsqJx0zs4KJCJauvJdlq9bR0iK2bu9m4vhWllx9N4vmzWDx/JlIanSYFTnpmJkVzNKV97J89Xq27eh5elnn9m4Alq9eD8AZhx/YkNgG4mc6ZmYF0rGli2Wr1tHZ1V1xfWdXN8tWraOjsyvnyAbHScfMrEBWrNlAS0v1W2ctLWLFXRtyiqg2TjpmZgWycfM2tm6vXMvptXV7Nxs3b8spoto46ZiZFci0KROYOL61apmJ41uZNmVCThHVxknHzKxAFsyaTk9P9TEze3qCBbOn5xRRbZx0zMwKpH1yG4vmzWBSW+XazqS2VhbNm0H7pObsr+OkY2ZWMKcf9nxm7b1zxXWz9t6Z0w97fs4RDZ6TjplZwVxww1rWPFR5oOc1Dz3BBTeszTmiwXPSMTMrEPfTMTOz3LifjpmZ5cb9dMzMLDfup2NmZrlxPx0zM8uN++mYmVmu3E/HzMxy4346ZmaWC/fTMTOz3LifjpmZ5cb9dMzMLDfup2NmZrlxP50RJGmJpJB0Z4V18yXdJqlT0iOSvihpaiPiNDNrFPfTGSGSXgB8APhzhXWHACuAPwBHA2cCxwDXSmqaz2BmlofF82eycO7+TBjXwqTxrQiYNL6VCeNaWDh3fxbPn9noEPuliOrVtFyCSInjVuAOYDYwNSJeVLL+dqANODgierJl84EfAidGxLeHcMxN7e3t7Zs2bRqJj2BmlruOLV2sWLOBjZu3MW3KBBbMnl73Gs7UqVPp6OjoiIgh3WkaN9IBDdHpwD7A4cDVpSsk7Q28FDijN+EARMRKSQ8BxwM1Jx0zs6Jrn9zGSXP2bXQYNWl40pE0AzgXeFNEPCH1aX8+K3tfU2Hzu0rWl+93oCpMey1xmpnZ8DX0eYhShvkScH1EXNVPsd2y98cqrHusZL2ZmTW5Rtd03gq8BPi7QZTt7+FTxeUD3W/MakKu7ZiZ5ahhSUfS7sB5wCeAp0qaP48DWrOftwKPZssr1Wh2pXINyMzMmlAjb6/tQ6ppfAJ4vOT1KtJzmseBJcDdWflKz25mU/lZj5mZNaFG3l77X+C1FZZ/FtgJWAg8GBF/lPRz4E2SPlvSZPpQYG/gO3kFbGZmw9OwpBMRTwI3lS/vbXUWEaXrPkDqk/MtScuAvYBPAT8Drqh3rGZmNjIK0Zs/Im4EjgL2A64FlmbvR0ZE9eFWzcysaTS69VofEXFIP8uvA67LNxozMxtJhajpmJnZ6OCkY2ZmuXHSMTOz3DjpmJlZbpquIYGZmQ1On6kNZk2nfXJzTt7Wy0nHzKxgIoKlK+9l2ap1tLSIrdu7mTi+lSVX382ieTNYPH8mFUbsbwpOOmZmBbN05b0sX72ebTuenmKMzu2py+Ly1esBOOPwAxsS20D8TMfMrEA6tnSxbNU6Orsq94vv7Opm2ap1dHR25RzZ4DjpmJkVyIo1G2hpqX7rrKVFrLhrQ04R1cZJx8ysQDZu3sbW7dVH/9q6vZuNm7flFFFtnHTMzApk2pQJTBzfWrXMxPGtTJsyIaeIauOkY2ZWIAtmTaenp7+JlJOenmDB7Ok5RVQbJx0zswJpn9zGonkzmNRWubYzqa2VRfNm0D6pOfvruMm0mVnBLJ4/E6BPP52enmDh3P2fXt+MFFG9mlZxI2kCsDuwMSK2j3hUOZC0qb29vX3Tpk2NDsXMbEj6jEgwe3rdazhTp06lo6OjIyKmDmX7mmo6kg4CzgdeDbQC84EbJe0BfAv4RETcMJRAzMysNu2T2zhpzr6NDqMmg36mI+lFwGrgAODrpesi4hFgEvCWEY3OzMxGlVoaEpwL/Al4AfBBoLx30o+AOSMUl5mZjUK1JJ25wJci4kmg0oOgB4G9RiQqMzMblWpJOhOBjirrdx5mLGZmNsrV0pDgPuDgKutfB/x2eOGYmdlgjfb5dP4bOFvS5cCvsmUBIOkM4AjgvSMbnpmZlRsr8+mcT2oifT3we1LCuUDSNOC5wErg4hGP0MzMnmVMzKeTdQKdD5wJdAJbgZnAX4CzgKMioqf/PZiZ2XCNqfl0ImJHRFwQES+JiOdExOSIeGFEfCYidtQrSDMzS8bEfDqSdpJ0n6T31TsgMzPr35iYTyfrm7Mb8GR9wzEzs2rG0nw6twEvqVcgZmY2sLE0n84HgRMknaoRaIsn6ZWSrpf0kKStkjZKulHSkWXlbpIUFV6XDTcGM7OiGUvz6SwFHgeWA+dJug/YUlYmIuLQQe5vF+Ae4KvAw9nPi4AVkk6KiNKkshY4uWz7v9QQe0MVsQOXmTWv0vl0eiLo6g7aWkWL1PTz6dSSdGaQ+uY8mP2853AOHBHXAteWLpN0DbCelHxKk86WiLhtOMdrhCJ34DKz5lc+H9pQ5kfL26CTTkTsV8c4eo+xQ1IH0JwNzGtU5A5cZta8eq8t27ufSTJd2b+b/dpSUz+depDUImmcpL0kfYTU4fSCsmIHSnpc0g5JayV9WFJT358qegcuM2tORb+21DRzKICknYHDSLfbANYBKyNi8xBjuBw4Pvv3E8AJEXFdyfrVpFttvwd2Ao4lze1zMHBclTgHmoe6fYjxDkotHbiKNvOfmTVO0a8ttU5XvRD4DOni3/upA3hS0uKI+PIQYjgL+BRp/LY3ApdLektEfAsgIs4uK/99SX8GPiTp1RFxyxCOWXdF78BlZs2p6NeWWqarPgZYBmwEFpPGYZsPnA48AiyTdHStAUTEuoi4IyKuiYiTSAOKfkFStdi+lr2/osp+p1Z7UX1uoGEregcuM2tORb+21PJM5yzgd8CLIuLCiPhR9roIOIh0++sDIxDT7aTm09OqlOmNu2kHGC16By4za05Fv7bUknReCFySDYnzLNnznK9lZYYs63R6CLAJeLRK0d4+O03bjLroHbjMrDkV/dpSa0OCak+vamogLumbwAPAL0gdPacDbyHNQPqerPn0XNJICFdmZZ8DvB44FbgiIn5SY/y5Ov2w53Pbuke54/7H+6ybtffOnH7Y8xsQlZkVXZGvLbXUdH4NvEXSc8pXSNoJOCUrM1g/BV5Lek70I+ALWTzHRMTnszK9Y3OfS+pIehnpVt5i4KQajtUQF9ywljUPPVFx3ZqHnuCCG9bmHJGZjQZFvrZosD1YJR0LfIc0JM1FwG+zVS8A3gM8D3hDRHyvDnGOOEmb2tvb2zdtGqhl9dB0bOlizn/c8KyOoeUmjGvh9n87rGmrwWbWfBp9bZk6dSodHR0dWYOsmtUyc+hVwLuBvYDPATdkr4uyZe8uSsLJQ9EnWjKz5lT0a0tNz3Qi4mJJ/01qKr0/6RnPfaTOoXVtglw0RW9Lb2bNqejXlppHJIiITcAVdYhlVOltS99Z5cvRzG3pzaw5Ff3aUkvn0BdLeleV9e+S9KKRCav4it6W3syaU9GvLbW0XjsH+Psq648E/n144YweRW9Lb2bNqejXllqSzkuBm6usvxmYM7xwRpfF82eycO7+TBjXwqTxrQiYNL6VCeNamn6iJTNrXkW+ttTSZHobqYXal/pZ/1bgcxExcQTjq5t6N5ku1Wfm0NnTm/avEDMrjkZcW4bbZLqWhgSPkPrk9GcW8NhQghjt2ie3NeUQ42ZWbEW8ttRye+0GYKGkPolH0t8Bp2VlzMzMKqqlpvMx4A3AHZK+AtxJGm/txcC/ANuBj454hGZmNmoMOulExH2SDgUuAd5Ztvpu4NSIaN4Bf8zMrOFqHZHg58CsrD/O80kjEtwTEbUM9GlmZmNUzSMSAETEnaTbazYIfVqYzJpO+2S3XjOz4SnitWXQTab7bCjNAE4E9iaNOP2ViOgcwdjqKo8m0xHB0pX3smzVOlpaxNbt3Uwc30pPT7Bo3gwWz59JmrfOzGzwGnltqWuTaUmnAf8KHBkRfypZPp80zcFk0i22AN4m6ZWVZhYdq5auvJflq9c/awjy3vGSlq9eD8AZhx/YkNjMrLiKfG0ZqMn0UcCOsoQj4IukhPMJ4BhS44JZwOn1CbN4OrZ0sWzVOjq7Kg/K19nVzbJV6+jo7Mo5MjMrsqJfWwZKOi8EVpYteyWwH3BpRHw4Ir4fEacBPwaOHfkQi6noc16YWXMq+rVloKQzDVhXtuxVpNtpl5ctX0GaPdQo/pwXZtacin5tGSjp7ADGly17afb+07LljwLNOYFDA/TOeVFNM895YWbNqejXloGSzv2k22kASGoF5gJrI+LxsrK7AX8Z0egKrOhzXphZcyr6tWWgpHMl8A+S3p2Nr/ZJ0i2371QoOwdYP8LxFVbR57wws+ZU9GvLQJ1DLwJOBi7MfhbwB+AzpYUktZMmeFs60gEWWe+cFpXa0jf7nBdm1ryKfG0ZsHOopCnAIlIjgfuA5RGxqazMy4E3Af8ZEb+tU6wjyvPpmFnRFXE+nSGPSFB0eSYdM7PRYrhJp5b5dMzMzIbFScfMzHLjpGNmZrlx0jEzs9w0LOlIeqWk6yU9JGmrpI2SbpR0ZIWy8yXdJqlT0iOSvihpSA+xzMyscRpZ09kFuAc4AziC1Cx7G7BC0om9hSQdQhrX7Q/A0cCZpJGtr5XkmpqZWYE0VZNpSeNIoxqsjYjXZctuB9qAgyOiJ1s2H/ghcGJEfHuIx3KTaTOzGuXaZFrSGyX9JLvF1V3htWMoQfSKiB1AB9CVHW9v0gCjl/YmnKzcSuAh4PjhHM/MzPI10DA4T5P0YeAjwJ+BW4HyAT+HJLtF1gLsAbwNmEm6hQZpYjiANRU2vatkvZmZFcCgkw7wTuAm4IiIGMkp6S7nmRrLE8AJEXFd9vNu2ftjFbZ7DDiov51KGui+WXstQZqZ2fDVcnttZ+DyEU44AGeRRqg+htRg4HJJJ5WV6e/BU/M8kDIzswHVUtP5FfBXIx1ARKzjmdlJr5F0DfAFSd8mTQwHz9R4Su1K5RpQ736rPuTKakKu7ZiZ5aiWms6HgbdL6veW1gi5ndScehpwd7as0rOb2VR+1mNmZk1q0DWdiLhZ0mnAbZJ+SppVtHyi7oiI04YajCQBhwCbgEcjYoeknwNvkvTZkibThwJ7U3kyuabTZ/jxWdNpn+ypDcxs7Bl0Px1JLwOuo/otqYiI6pN3P7O/bwIPAL8gTXM9HXgLqaPoeyLi81m515H65FwJLAP2Aj4FPAi8KiLKE99gP0/d++lEBEtX3ltxoqVF82aweP5MUp41MyuG4fbTqeWZzoWk/jOvB1aXT+Q2BD8lTfz2NlIi6wB+DhwTEdf0FoqIGyUdRWqufS2wGbgKOGuoCScvS1fey/LV69m24+kuRnRuTyEvX51m9j7j8AMbEpuZWSPUUtPZAiyJiPPqG1I+6l3T6djSxZz/uOFZCafchHEt3P5vh3kWUTMrjDxHJHgE2D6Ug4xFK9ZsoKWl+q2zlhax4q4NOUVkZtZ4tSSdrwBvzsZHswFs3LyNrdur3/3bur2bjZu35RSRmVnj1ZJAbgGOIrVeu5g0MGefq2pErBqh2Apt2pQJTBzf+vQznEomjm9l2pQJOUZlZtZYtSSdG0r+vZy+owEoWzao1muj3YJZ01ly9d1Vy/T0BAtmT88pIjOzxqsl6ZxatyhGofbJbSyaN4Plq9fT2dW3tjOprZWFc/d3IwIzG1Nq6Rz6tXoGMhotnj8ToGI/nYVz9396vZnZWNFUk7jlKc9J3PqMSDB7ums4ZlZIeXYOBUDSnsBLSOOj9Wn9FhFfH0ogo1n75DZOmrNvo8MwM2u4WiZxawG+ACykelNrJx0zM6uoln46Z5KGrPkWaYw0AR8E3gWsJQ1hM3+kAzQzs9GjlqTzFuD6iDgZ+EG27BcR8V/AwcDu2buZmVlFtSSdGTyTbHoHFGsDiIingK+Sbr2ZmZlVVEvS6SSNMg3wJKkj6B4l6x+mDjOLmpnZ6FFL67UHgAMAIqJL0v+S5r65NFt/GPDnkQ1vdPAkbmZWD0W8ttQytcFngGMj4oDs5w8D5wI3kxoVzAXOj4gP1CnWEeVJ3MysqBp5bcmzn875wA8lTYiIbcAnSLfX3kwa+HMZcM5QghitPImbmdVDka8tg36mExEbIuL6LOEQEd0R8a8RsWtETIuId0TE1vqFWiwdW7pYtmpdxXHXADq7ulm2ah0dnV0V15uZVVL0a0stDQkGJKl9JPdXZJ7EzczqoejXlkEnHUk/kvTcKutfBdw5IlGNAp7EzczqoejXllpqOq8Efi1pQelCJWcDP65xf6Na7yRu1XgSNzOrVdGvLbUkiZcBjwLXSFoqqU3S3sCNwEeA7wMvqkOMhbRg1nR6eqq3DPQkbmZWq6JfW2ppSPAb0jA3lwDvI421dicpGb07It4QEY/XI8gi6p3EbVJb5b9IJrW1smjeDE9xYGY1Kfq1paapDSKiU9LbgZnAq0ijErwnIi6uR3BF50nczKweinxtqWkSN0kHAJcBBwH/Dbwa2IfUSfRjUaAZ4TyJm5kVXSOuLcPtHFrLiARvAi4mDfa5KCKuyJpIfxl4A2lkgjdGRHO20yuTZ9IxMxsthpt0amlIcCnwW+DFEXEFQER0RMQ/AO8E5gC/HkoQZmY2NtSSdD4NzI2I+8tXZHPqvAwP+GlmZlUMuiHBQAN5RsQaSS8dfkhmZjZajVhnTkmTgb1qKH+opEsk3SNpi6Q/SvqOpNll5W6SFBVel41U7GZmlo+qNR1J24GTI+Ky7OcpwDeBf4uIu8qKHwd8HajeVfYZbwd2Ay4AfgfsCZwF3CHpkIi4raTsWuDksu3/MsjjmJmNSkWcT2eg22vjeHZtaDxwFPDZETj2uyLikdIFkn4IrAfeDxxfsmpLWRIyMxuz+ptPZ8nVdzf9XF01dQ4dSeUJJ1u2SdJaUt8fMzOrYEzMp5MHSdOAWcCaslUHSnpc0g5JayV9WFJz1yHNzOqg6PPpNKymU06pLriMlAjPL1m1mjQKwu+BnYBjSSMgHEx6jtTf/gbq9em5f8yscGqZT+ekOfvmFNXgNU3SIfUDOhY4NSJ+17swIs4uK/d9SX8GPiTp1RFxS55Bmpk1UtHn0xlM0llQMnnbZNIgn/8oqXwag4OHGoSkjwNnAO+NiEsGscnXgA8BrwAqJp2BhmjIakKu7ZhZofTOp9NZJfE083w6g0k6b8xepd7WT9maB/yUdC4pgZwVERcNcrPeZ1E9VUuZmY0yC2ZNZ8nVd1ct08zz6QyUdF5bz4NLOgc4Gzg7Ij5dw6a9fXbcjNrMxpTe+XSWr15fsTHBpLZWFs7dv2lHsq+adCLi5nodWNIZwBLSjKM3SHp5yeptEfErSXOBDwJXAg8AzwFeD5wKXBERP6lXfGZmzWrMzKczogeWbgJe08/qByJiP0nPAy4EXgjsTrqddg/pmc7nIqL607Tqx/fUBmZWaKN6Pp3RxknHzKx2ec6nY2ZmNixOOmZmlhsnHTMzy42TjpmZ5cZJx8zMcuOkY2ZmuXHSMTOz3DjpmJlZbpx0zMwsN046ZmaWGycdMzPLjZOOmZnlxknHzMxyM5iZQ22Y+gw/Pms67ZObc4IlM7N6ctKpo4hg6cp7+0y0tOTqu1k0bwaL589EUqPDNDPLjZNOHS1deS/LV69n246ep5d1bk/zzi1fvR6AMw4/sCGxmZk1gp/p1EnHli6WrVpXcQ5zgM6ubpatWkdHZ1fOkZmZNY6TTp2sWLOBlpbqt85aWsSKuzbkFJGZWeM56dTJxs3b2Lq9ci2n19bt3WzcvC2niMzMGs9Jp06mTZnAxPGtVctMHN/KtCkTcorIzKzxnHTqZMGs6fT0RNUyPT3BgtnTc4rIzKzxnHTqpH1yG4vmzWBSW+XazqS2VhbNm0H7JPfXMbOxw02m62jx/JkAffrp9PQEC+fu//R6M7OxQhHVbwGNVpI2tbe3t2/atKnux+ozIsHs6a7hmFkhTZ06lY6Ojo6ImDqU7V3TyUH75DZOmrNvo8MwM2s4P9MxM7PcOOmYmVlunHTMzCw3TjpmZpabhiUdSYdKukTSPZK2SPqjpO9Iml2h7HxJt0nqlPSIpC9KGlLLCTMza5xG1nTeDuwLXAAcCSzOfr5D0st7C0k6BFgB/AE4GjgTOAa4VpJramZmBdKwfjqS9oiIR8qWTQXWAzdGxPHZstuBNuDgiOjJls0HfgicGBHfHuLxc+unY2Y2Wgy3n07DagrlCSdbtglYC+wDIGlv4KXApb0JJyu3EngIOD6faM3MbCQ0VedQSdOAWcC3skWzsvc1FYrfVbK+0r4GqsK01xygmZkNS9M8E5EkYBkppvOzxbtl749V2OSxkvVmZlYAzVTT+TRwLHBqRPyubF1/D576fSA10P3GrCbk2o6ZWY6aoqYj6ePAGcB7I+KSklWPZu+VajS7UrkGZGZmTarhNR1J5wIfAs6KiIvKVt+dvc8itVYrNRu4tc7hjYg+o0zPmk77ZI8ybWZjT0OnNpB0DrAEODsiPtZPmTtINbKXljSZPhS4ATgpIi4b4rHr3mQ6Ili68t6K8+ksmjeDxfNnkh5lmZkVQ2GnNpB0BinhfB+4obRDKLAtIn6V/fsDpFrOtyQtA/YCPgX8DLgiv4hrt3TlvSxfvZ5tO55u7U3n9m4Alq9eD8AZhx/YkNjMzBqhkc90js7ejwJ+Wvb6bm+hiLgxK7MfcC2wNHs/MiK6c4y3Jh1buli2ah2dXZVD7OzqZtmqdXR0duUcmZlZ4zSsphMRh9RQ9jrguvpFM/JWrNlAS0v1W2ctLWLFXRs8wZuZjRlN0XptNNq4eRtbt1eviG3d3s3GzdtyisjMrPGcdOpk2pQJTBzfWrXMxPGtTJsyIaeIzMwaz0mnThbMmk5PT/WWgT09wYLZ03OKyMys8Zx06qR9chuL5s1gUlvl2s6ktlYWzZtB+yT31zGzsaPhnUNHs8XzZwJU7KezcO7+T683MxsrGto5tJHynE+nz4gEs6e7hmNmhVTYzqFjSfvkNjeLNjPDz3TMzCxHTjpmZpabsfxMpwdQe7un1DEzG6yOjg6AiIghVVrGctLZQarpPZHTIXuzW0dOxysCn5PKfF4q83mpLO/zsjPQExFDahMwZpNO3rKZSgec0XQs8TmpzOelMp+Xyop2XvxMx8zMcuOkY2ZmuXHSMTOz3DjpmJlZbpx0zMwsN046ZmaWGycdMzPLjfvpmJlZblzTMTOz3DjpmJlZbpx0zMwsN046gyRpJ0kXSdogqVPSzyUdM8htD5B0laQOSZslrZD0d/2U/VdJ90raJuk+SWdJatrf01DPi6SFkq6W9EC23dpsP9MqlI1+Xm+vz6cavmGclyX9fNaH+yk/Vr4v91f5Hvy+rGyhvi+S9pF0oaRbJD2ZxXpIDdsX6/oSEX4N4gWsBB4FTgNeB3wd6AYWDLDdHsCfgDuBY4GjgNuAjcA+ZWU/nO3zXOAQ4ENAF/DJRn/+OpyXh4BvAG8EXgO8A3gYuB+YWlY2gMuAl5e99mj056/DeVmSfd7Dyj7rQRXKjqXvy4sr/P7fmp2rT5aVLdT3JfvdPQJcB3wvi/+QQW5buOtLw094EV7AguyLcFzJMgG3AL8bYNvzgE5gr5Jlu5GmVPjPsmWdwIVl2388+2LsM9zP0WTnpc8FIEs+AbynbHkAn230583pvPQmnakDlBtT35d+9ndRtr+ZBf++tJT8+9gak07hri9NWw1vMseR5qr4Xu+CSL+xrwF/019VtmTblRHxp5JtHwWuAd5QUu4IYGK2z1KXAOOAQd3Ky9mQz0tEPFJh8R3Z+z4jGWQDDOf7Mlhj6vtSTtJ4Ui35loi4d6QDzVNE9Axj88JdX5x0BmcW8NsKX47flKzvQ9Ik4ABgTYXVvwH2kLRHyT4CuLu0UESsJf2FUvEYDTak81LF67L3Sufr5OwZwFZJP5N0Qo37ztNInJffSerOnn18qeR7UnqMsfx9OZb01/tX+llfpO/LkBT1+uKkMzi7AY9VWP5YyfpKdiHdPhjMtrsBWyJiW4Wyj1c5RiMN9bz0IWlX0u2StcDlZau/CbwbOBw4mfSf5NuS3ltrwDkZznm5j3Sv/VRgPnAxcCJwm6Rdyo4xZr8vwL8AT9L3uwLF+74MVSGvL0OabnSMqjZ0w0DDOgx22+Eco1GGHbOkycBVwK7AvPL/GBHx5rLy/wPcBHxM0rKI6Kwp4nwM6bxExKVli26UdBvwQ+BdwMeGe4wGG4nvyz6khPzViHiqz06K+X0ZjkJdX1zTGZxHqfyXwK7Ze6W/NCD9BRGD3PZR4DmSJlQou0uVYzTSUM/L07JbBFeTWictiIjfDLBJ7z3wbwA70Zy3kYZ9XkpFxEpgA/CKsmOMue9L5hTStau/W2vPUoDvy1AV8vripDM4dwN/W6E9++zsvdI9VbK/qNZR+Ys+G9hY8kD9blJV+QWlhSQ9D5jU3zEabEjnpZekiaSHyq8AjoqIW2s4du8xh/MQtl6GdV760cKzP+uY+74ASBIp6fx+FH1fhqSo1xcnncH5LjAVOLps+cnAPRHx2wG2nS/pub0LsucXRwPfKSn3A2Ab8M9crIWtAAAJN0lEQVRl278F2EFqjdJshnxesr+4rgLmAq+PiJsHe9DsovUmYDNlD0abxHC+L31IOhzYk9T/oteY+r6UeA3p4fmgajlQiO/LcBTv+tLoNupFeJH+QrgR+AvpAeZrSU0Ne4CjS8rdRNYKtGTZnqROj78EXg/8PfBTUnV337Ky55C+AEtI/7k+CGwHPt3oc1CH83IN6dbAR+jbie+AknJnAl8CTiJ1aDsRuDnb9p2NPgd1OC+/Ak4n9WmZn30nNpMaWJR3mh0z35eSdV8n9SvZs5/1hfu+ZHH/Q/b6VBbrOdnPRw7wfSnc9aXhJ7soL2Bn4PPZL3hr9ks+tqxMxf8swPNJt5GeILW4+QHwggrlBLwvu8BsA9YD/4+SzmPN9hrqecn+Y/X3uqSk3NHAalIP6y5gE/Cj0otUM76GcV6+lf3+n8ouCPcBFwC7juXvS7Z8SnZevldl/0X9vvT3f+H+QZyXQl1fPJ+OmZnlxs90zMwsN046ZmaWGycdMzPLjZOOmZnlxknHzMxy46RjZma5cdIxM7PcOOlY05F0SNnc9t2SHpe0RtLXJB2RjcFVCJJ2kfTvku6QtEnSdkl/lHSlpDcM9bNIep+kU0Y4XLO6cudQazqSDgF+TOqdv4LUk3oKcCBp8q59gRuAf4yITQ0Kc1AkzSH1Ft+DNJr2zaSe43uThrp5OfCuiLh4CPu+n9Rj/ZCRites3jyfjjWzX0bEN0oXSFpMmhd+MSkpHVltB5LagNaI2Fq3KPs/9nNJY8xNBF4TEbeUFfmopP9LGlp+TJDUCkyIiC2NjsUaw7fXrFAiojsizgBuAY6Q9OredZKWZLfjXiBpqaQ/ksb3ermk/bJ1S8r3WbLdfmXLXyPpp9m0xw9LujDbd8X9VPB+Ug3nAxUSTu/nuT4iLis55j9JulrSg5K2SfqLpKsk/Z+y2AL4a+A1Zbci9ysp8xJJ3832sU3SPZL+TVKfPzYlHS/p19n0zg9KOkfSYdk+Tykru7ukL0j6Q3ar8A/Zz7uVlTsl2/4wSWdLuo/0+zghO9aDFaY5QNIJ2XblIyLbKOCajhXVl4FXk0bVLb+gf5M0RfFnSIMmbqh151ky+yFpoqxPkgaOPAF4VQ27OZ40aOfXatjm3aQJtZaRBsU8AFgE/ETSQZHmtIc0RP0FpBGbP16y/cYs/gWkYe//l3QeHiPNW3Qu8CLgH3s3kPRPpFrjfaRRv3eQhrwvn4IASe3ArcDzSNML/JI0Ad87gNdJmhMRm8s2Ox9oI43+/ARwT/bvz5FG0r6+rPy/AB3A/1Q7UVZQjR5d1S+/yl+kIekDOLNKmYOyMleWLFuSLbsJGFdWfr9s3ZIK++rdbr+SZbeT/iqfUbKsDfhJf/sp2+eUrNxvavzsz6mw7G9JowJfXLb8fuCmCuUnkhLWqgrn4fQsrkOyn8cBDwF/BnYpKbcTaYKwAE4pWf5xKkwTQJpKO4CPliw7JVt2DzC5rHw7acToy8uW/xXQXf5Z/Ro9L99es6J6InvfucK6z0bEjqHuWNKewEtJQ+iv610eEV3AhYPcTW9cT1QtVSYinspikKSdJe1Oqr3cA7xskLuZT5pn5avA1Ox22O7ZvlZkZQ7P3g8G9iJNJ/F4SRxPAv9VYd/HZfEsK1v+RVKt67gK2/xnlD3DiYgO4Arg9VlcvU4l3fb/8oCf0grJSceKqtpF/d5h7nv/7P2eCusqLaukN64ptRxY0oslfZ80cVsH6QK/kTT98GAbHPxt9v6Vku17X7/P1u2Zvdf6WfcnzfL5rKSe/XwPMKPCNv39PpYB44E3w9NTUZ8K3BkRv+hnGys4P9Oxoup9sF7pwlipZVS1vgHl/w+G3QcoIjZLegD4G0mTIs1nX5WkfUm3xJ4APkr6bE+RYv8s6ZbXYPTG/37gzn7K/KmsbD1VbKkWEbdKWgOcRvp8h5Jug747h5isQZx0rKhOy96vHWT5x7L3XSusK//rvPeW2oEVylZa1p/vkJ6h/DN9b0dVchwpsRwTET8uXZG1DNtWVr6/RNrb2OCpiLhhgGOuz94H+1nXAQdKGlda28laxM3kmXM3WF8CLsz6M51Geo72zRr3YQXi22tWKJJaJZ1Parm2IiJ+MpjtIrWoepjUwurpv+4lzSB1OC0t+2fg56TnDTNKyrYB760h3PNIt7TOk/SKfj7P4ZJOzH7s7l1cVuatwHMrbP4klZPo9cAjwAcl9VkvaZKk3tt+Pye17jtF0i4lZXYC3l5h31cB04CFZcvfmi3/boVtqrmUlGjeT0q6V0aTd/i14XFNx5rZQZLenP27dESCvyY1Z35jjfv7PPAx4AeSriI9QH87sIbUcKDUmcBK4FZJF5Oer5xAegYB1W/XpQIRD0s6ijQiwS3ZMXtvn+0FHEFKnu/INvkB6VbUpZI+T2qu/SrSyAX30ff/623AaZI+CvwO6AGuiYinJJ1MShD3SPoKqen0VOBvgDeQLvA3RcQOSWeSahe3S/oyqcn0KcCjpGc4pZ/1PFJz6y9IOgj4FanJ9Gmk24HnDXReys7R45L+h+y5DrC8lu2tgBrdfM4vv8pfPNNkuvfVTbro303q83JEP9stoazpc9n6caSL4gbSX9e/JPVFqbgd8DrShX0rqUnxhaQWZAGcVcPn2RU4h1Sr6CD13fkjqR/KMWVl55H6HW0m9Q26FphFagZ+f1nZPYArSbcOe8o/Q7bdN0hNordnn+FW4Gxg17J9nQD8hnQL78Es3uOyfZ5QVnYacHH2Gbqy9y8Au5eVO4WS5tlVzs/crNxasqG5/Bq9L4+9ZlYDSceTksVJUTKSwGgk6QxSx85XRMRtdTzOHOBnwIci4hP1Oo41Bycdswqy5z4TomTMtuyZzk3AHOCvIuLhBoU3oiSNB7ojortk2U6kms/OwF4Rsb2Ox/86cCKw72g5p9Y/P9Mxq2wC8ICkb5KeVewG/BOpqfanRtnFcQbpOddlpNZs00nD4OwPvKMeCUfSc0i3Nl9Aep6zbJSdU+uHazpmFSiNhvwl4DWki7BIyWdZDGEagmaWNcf+PKnRwh6khgR3ARdExOV1OuZ+pAT3JKkBxcKIqGn0BismJx0zM8uN++mYmVlunHTMzCw3TjpmZpYbJx0zM8uNk46ZmeXGScfMzHLz/wGzakWesTZHtwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(drugged, exam, 'o')\n",
    "plt.xlim(-0.1, 1.1)\n",
    "plt.xlabel('Drug Category')\n",
    "plt.ylabel('Exam Score')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "It looks like we can indeed regress this data! The equation we're modeling is this:\n",
    "\n",
    "$$y = \\beta_0 + \\beta_1 \\delta_d + \\epsilon $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.014285714285715 28.28571428571429\n"
     ]
    }
   ],
   "source": [
    "dof = len(drugged) - 2\n",
    "cov = np.cov(drugged, exam, ddof=2)\n",
    "slope = cov[0,1] / cov[0,0]\n",
    "intercept = np.mean(exam - slope * drugged)\n",
    "\n",
    "print(slope, intercept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEiCAYAAAAiQw8CAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucHFWZ//HPM5PJjSSTBBIyCbIQJChMRILiNYEVwg+yiCC7LKg/LkuM91UDoj8VyaKuK2IQFXYNEcHLgrAggkQwyGKCiIhcEyCJJNxCIIFkJrfJzGTm+f1xqklPT3VP98x0ddfM9/169aszVaeqTtd06plT9ZxzzN0RERFJQk2lKyAiIoOHgo6IiCRGQUdERBKjoCMiIolR0BERkcQo6IiISGIUdEREJDEKOiIikhgFHRERSYyCjoiIJEZBR0REEqOgIyIiiVHQERGRxCjoiIhIYoZUugKVYma7CUF3a6XrIiKSImOATnfvVfywwTqfjpl1AlZfX1/pqoiIpEZzczOAu3uv7pQN2pYOsLW+vr6+qamp0vUQEUmNsWPH0tzc3Os7RHqmIyIiiVHQERGRxCjoiIhIYgbzMx0RkVRr3tnOkhUb2LStlQmjhzGnsYH6kXWVrlZBCjoiIinj7ixcuppFy9ZSU2Psautg+NBaFty2knmzpjJ/9jTMrNLVjKWgIyKSMguXrmbx8nW07u58fVlLWwcAi5evA+D84w+pSN16omc6IiIp0ryznUXL1tLS3hG7vqW9g0XL1tLc0p5wzYqjoCMikiJLVmygpqbwrbOaGmPJExsSqlFpFHRERFJk07ZWdrXFt3IydrV1sGlba0I1Ko2CjohIikwYPYzhQ2sLlhk+tJYJo4clVKPSKOiIiKTInMYGOjsLj5nZ2enMmd6QUI1Ko6AjIpIi9SPrmDdrKiPq4ls7I+pqmTdrKvUjqrO/joKOiEjKfP64g2mcMiZ2XeOUMXz+uIMTrlHxFHRERFLm8rvXsGJ9/EDPK9Zv5fK71yRco+Ip6IiIpIj66YiISGLUT0dERBKjfjoiIpIY9dMREZHEqJ+OiIgkRv10REQkUeqnIyIiiVE/HRERSYT66YiISGLUT0dERBKjfjoiIpIY9dMREZHEqJ9OPzKzBWbmZvZozLrZZvaAmbWY2UYz+5GZja1EPUVEKkX9dPqJmR0GfBF4JWbdMcAS4AXg/cAFwMnAHWZWNZ9BRCQJ82dPY+7MAxk2pIYRQ2sxYMTQWoYNqWHuzAOZP3tapauYl7kXbqYlUokQOO4H/gJMB8a6+1uz1j8I1AFHuntntGw28DvgDHf/ZS+O2VRfX1/f1NTUHx9BRCRxzTvbWbJiA5u2tTJh9DDmTG8oewtn7NixNDc3N7t7r+40DenvCvXS54H9gOOB27JXmNkU4O3A+ZmAA+DuS81sPXAaUHLQERFJu/qRdZx51P6VrkZJKh50zGwqcAnwYXffatYt/7wxel8Rs/kTWetz99tTE6a+lHqKiEjfVfR5iIUIczVwl7vfmqfY3tH75ph1m7PWi4hIlat0S+ejwNuAQ4som+/hU+zynu43Ri0htXZERBJUsaBjZvsAlwLfAnZkpT8PAWqjn3cBr0XL41o044lvAYmISBWq5O21/QgtjW8BW7Je7yE8p9kCLABWRuXjnt1MJ/5Zj4iIVKFK3l77G/D3Mcu/B4wC5gLPu/uLZvYQ8GEz+15WyvSxwBTglqQqLCIifVOxoOPu24F7c5dnss7cPXvdFwl9cq43s0XAZODbwJ+Bm8pdVxER6R+p6M3v7vcAJwEHAHcAC6P3E9298HCrIiJSNSqdvdaNux+TZ/mdwJ3J1kZERPpTKlo6IiIyMCjoiIhIYhR0REQkMQo6IiKSmKpLJBARkeJ0m9qgsYH6kdU5eVuGgo6ISMq4OwuXrmbRsrXU1Bi72joYPrSWBbetZN6sqcyfPY2YEfurgoKOiEjKLFy6msXL19G6+/UpxmhpC10WFy9fB8D5xx9Skbr1RM90RERSpHlnO4uWraWlPb5ffEt7B4uWraW5pT3hmhVHQUdEJEWWrNhATU3hW2c1NcaSJzYkVKPSKOiIiKTIpm2t7GorPPrXrrYONm1rTahGpVHQERFJkQmjhzF8aG3BMsOH1jJh9LCEalQaBR0RkRSZ09hAZ2e+iZSDzk5nzvSGhGpUGgUdEZEUqR9Zx7xZUxlRF9/aGVFXy7xZU6kfUZ39dZQyLSKSMvNnTwPo1k+ns9OZO/PA19dXI3Mv3EyL3chsGLAPsMnd2/q9Vgkws6b6+vr6pqamSldFRKRXuo1IML2h7C2csWPH0tzc3OzuY3uzfUktHTObAVwGvBeoBWYD95jZROB64FvufndvKiIiIqWpH1nHmUftX+lqlKToZzpm9lZgOXAQ8NPsde6+ERgBnN2vtRMRkQGllESCS4CXgMOALwG5vZN+DxzVT/USEZEBqJSgMxO42t23A3EPgp4HJvdLrUREZEAqJegMB5oLrB/Tx7qIiMgAV0oiwTPAkQXWvw94sm/VERGRYg30+XT+G7jIzG4EHomWOYCZnQ+cAHy2f6snIiK5Bst8OpcRUqTvAp4mBJzLzWwCMAlYClzV7zUUEZEuBsV8OlEn0NnABUALsAuYBrwKXAic5O6d+fcgIiJ9Najm03H33e5+ubu/zd33cveR7n64u3/X3XeXq5IiIhIMivl0zGyUmT1jZp8rd4VERCS/QTGfTtQ3Z29ge3mrIyIihQym+XQeAN5WroqIiEjPBtN8Ol8CTjezc60fcvHM7N1mdpeZrTezXWa2yczuMbMTc8rda2Ye87qhr3UQEUmbwTSfzkJgC7AYuNTMngF25pRxdz+2yP2NA1YBPwFejn6eBywxszPdPTuorAHOytn+1RLqXlFp7MAlItUrez4d890M7WihpXYUNWZVP59OKUFnKqFvzvPRz/v25cDufgdwR/YyM7sdWEcIPtlBZ6e7P9CX41VCmjtwiUgVcofNa2H9wxy97vccXfMwh9qz3Gbv4iL/GL2ZHy1pRQcddz+gjPXIHGO3mTUD1ZlgXqI0d+ASkSqw9SVY/zC89HD0/gjsasKIHrBHD0gOr1lLe1sIONV+ban4dNVmVkM4dROBjxE6nF6QU+wQM9sCjCa0hK4Dvu3uVRucMh24sgNOtkwHrrkzq/feq4gkaOfmKLg8sifIbH85tugOH84TfiCPdU7l8c6DeMynvr6u2q8tJQcdMxsDHEe43QawFljq7tt6WYcbgdOif28FTnf3O7PWLyfcansaGAWcQpjb50jg1AL17Gke6vpe1rcopXTgStvMfyLSR63bYcNjWS2Yh2HLs/Fla4fCpOkweQZMmcEdrzVw4R9a2NGW/1ZaNV9bSp2uei7wXcLFP3NFdWC7mc139x/3og4XAt8mjN/2IeBGMzvb3a8HcPeLcsr/xsxeAb5sZu919/t6ccyyS3sHLhHpJ7tb4ZUVe26PrX8YXl0FcaOGWQ1MeDNMOeL1IMPEw2DI0NeLPPP7NexsW13wkNV8bSk66JjZycAiQsvma8CKaNVhwGeARWa20d1vL6UC7r422ifA7VEywZVm9ssCY7ldB3wZeBcQG3TcfWwPn6eJMrZ2Mh24WgoEnmruwCUivdDZAZtWdW3BvLwCOvM8CRg/dU9wmTwDGt4CQ/cqeIi0X1tKaelcCDwFvCMaoSDj92b2E0Ln0S8CJQWdGA8CJwETgFfylMn0L6raAUbnNDaw4LaVBctUcwcuEemBO2xZ17UFs+ExaN8RX3705Ci4HLHnfcS4kg+b9mtLKUHncOCSnIADgLtvM7PrgNxbYSWJOp0eAzQBrxUomumzU7Vp1JkOXIuXr4sdDXZEXS1zZx5YlQ/6RCTG1g1dWzAvPQItW+LLjhjXtQUzZQaMntQv1Uj7taXURIJCT8ZLShA3s18AzwF/JXT0bADOJsxA+pkofXomYSSEm6OyewEfAM4FbnL3P5ZY/0R9/riDeWDta/zl2e5fzMYpY/j8cQdXoFYi0qOdm0NQyc4m25Zn1Oa6vWDyW7NaMDNg3AFQxj54ab62lBJ0HgPONrMr3b1L+9HMRgHnRGWK9Sfgw4Q06XqgGXgIODnruVDmt3wJsA/hdtoqYD7wgxKOVRGX372GFeu3xq5bsX4rl9+9pmpz6UUGjbYd4bZYdn+YLeviy9YOhX0bu7Zg9pkGNYUH4Oxvab62WLE9WM3sFOAWwpA03weejFZlEgneCHzQ3X9dhnr2OzNrqq+vr29q6imzunead7Zz1L/fnbefDsCwITU8+JXjqrYZLDLg7G4LmWTZLZhNTxfIJHtTFFyibLJ9D4MhlX1AX+lry9ixY2lubm7uKVkrn1JGJLjVzD5NSG/+AXtupxmwA/h0WgJOEtRPR6TCOjvg1dVdWzCvrICOtvjy4w7s2oKZ9BYYNirZOhch7deWkp7puPtVZvbfhGmrDyQEnGcInUOby1C/1FI/HZEEuYfOldnDxWx4DNryTAE2uqFrC2byETByfKJV7q20X1tKHpHA3ZuAm8pQlwEl7bn0IlVt28vdxyRr2RxfdvjYri2YyTNgTHWmExcj7deWUjqHHgG8292vzLP+U8Af3f3R/qpcmqU9l16karRs2dMPJvO+7aX4snV7QcPhXfvDjDuwrJlkSUv7taWUls7FwFAgNugAJwLHAh/sa6UGgrTn0otURNvO7mOSbV4bX7amDiY1dm3BTDgk8UyypKX92lJK0Hk7IWstnz8An+1bdQaW7ImWsufT6ez0qp9oSaTsdrfBxpVZt8kegU1PxWeSYSGTLLsFs29jxTPJKiXN15ZSUqZbCRlqV+dZ/1HgB+4+vB/rVzblTpnO1m3m0OkNVftXiEhZdHaGTLLcMck68jzsHndAzphkh1dlJlmlVeLakljKNLCR0Ccnn0Ygz5O8wa1+ZF1Vpi6KlIU7ND0fBZi/hhbMhkfzZ5KNmpT1oD/KJktJJlmlpfHaUkrQuRuYa2ZXu3uXp1hmdihwHqHzqIgMJts35mSSPQw78wydOLy++5hkYyYnW1+pqFKCzjcISQJ/MbNrgEcJHUSPAP4FaAO+3u81FJHq0dIUWi3Zz2G2vhhftm5kuC32epA5IgzlP4AyyaR0pYxI8IyZHQtcC3wyZ/VK4Fx3X9OPdRORSmrbCS8/0bUF89rf4svW1IUhYrqMSXYI1JbcFVAGuFJHJHgIaDSztwIHE0YkWOXupQz0KSLVpqMdNj7ZtQWz8UnwuA6IFga5nDIDphy5Z0yyulTkEEmF9erPkKgDqDqBFqlbhkljA/Ujlb0mFdLZGVosXTLJnoDdu+LLj/27PbfHJs8Iw/gPG51snSVWGq8tRadMd9vQbCpwBjCFMOL0Ne7e0o91K6skUqbdnYVLV8fm0s+bNZX5s6dhur8t5eQOzS90fdC/4TFojR8Wn70mZt0iOzIEmr32TrbO0qNKXlvKmjJtZucB/wqc6O4vZS2fTchUG0m4xebAx8zs3XEziw5WC5euZvHydV2GIM+Ml7R4eZivo1rnvJCU2r6pawtm/cOw89X4ssPqQ6ulSybZFD3oT4E0X1t6ur12ErA7J+AY8CNCwPkWYcroUwmzeX4eZbABodm7aNnavHNetLR3sGjZWubOnKqOotI7u5rhpUe7DnrZ/EJ82SEjoOEtXdOVx0+Fmppk6yx9lvZrS09B53Dgxpxl7wYOAH7q7l+Nlv3GzA4ATkFBB0j/nBdSZdpbwnOX7BbMa3mSRWuGwMRDu7ZgJrxZmWQDRNqvLT19CycAuaPtvYdwOy03GC0BvtZP9Uq9tM95IRXU0Q4bn+p6m2zjU9C5O6awwT4Hd23BTJquTLIBLO3Xlp6Czm7CyNLZ3h69/yln+WvA4Bx9L0ba57yQhHR2wuZnurZgXn48fyZZ/f57hoqZMgMa3grDxyRbZ6motF9bego6zxJup/0QwMxqgZnAGnffklN2byDPE8vBJ+1zXkgZuEPzi11bMC89Bq15Jt3da0L3IWP22ifZOkvVSfu1paegczPwNTO7H7iHkCwwAbgmpuxRwLr+rV56pX3OC+kHO17tPibZjk3xZYeNCZlk2UGmfj9lkkk3ab+29BR0vg+cBVwR/WzAC8B3swuZWT3wD8DC/q5gmqV5zgsp0a6t3ccka34+vuyQ4TDpLV1bMOMPUiaZFC3N15YeO4ea2WhgHvBG4Blgsbs35ZR5J/Bh4D/d/cky1bVfaT4d6bX2Xd3HJHt1DSG/JofVwr6Hdm3BTHwz1Or3L32Xxvl0ej0iQdolGXQkxTp2h9kss2+TbXwyTyYZsPfBXVswk6ZD3Yhk6yxSRklO4iYysHV2wua1XVswGx6H3XlGd6p/w56pkzNjkg2vT7bOIimjoCODkztsXZ/zoP/R/JlkI/fp2oKZPANGTUi2ziIDgIKODA47XgvDxGS3Yra/El926OjuY5LVv0GZZCL9QEFHBp7WbWEk5exWTNNz8WVrh3Ufk2zvNyqTTKRMKhZ0zOzdwMVAI6Fj6TbgCeA77v7bnLKzCWO6HR6V+xXwxdwsOhmEdrfCyyu6tmA2rSJvJtnEQ7v26J94qDLJRBJUyZbOOGAV8BPg5ejnecASMzvT3W8AMLNjCOO63Qp8FZgMfJswg+lMd48falUGno7d8Oqqri2YV1ZCZ3t8+b3f2H1MsqEjk62ziHRRVSnTZjaEMKrBGnd/X7TsQaAOODITYKKWz++AM9z9l708llKmq5l7lEn2SFYm2WPQvjO+/Jj9uo9JNqJXGZ0iUkCiKdNm9iHgU8DBhFtiudzde916cvfdZtYMtEfHm0IYYPT87BaNuy81s/XAaUCvgo5Uma0v5WSSPQK78vxBMHLv7mOSjZqYbH1FpFeKDhBm9lXg34BXgPuB3AE/e8XMaoAaYCLwMWAacEG0ujF6XxGz6RNZ6yVNdm7eM1RMJshsfzm+bCaTLLs/zNj9lUkmklKltEo+CdwLnODueW6i98qNhBYLwFbgdHe/M/o505raHLPdZmBGvp2aWU/3zdSLLwmt28NtsewH/VuejS9bOyw8d8luwex9sDLJRAaQUoLOGODGfg44ABcSEgMmAR8CbjSzs939+qwy+R48Vc8DKQmZZK+s2HN7bP3D4cF/XK6H1YTMsewWzMRDYUju9E0iMpCUEnQeAd7Q3xVw97XsmZ30djO7HbjSzH5JmBgO4p8fjSe+BZTZb8GHXFFLSK2d3ursCKnJ2S2YV1ZCR1t8+fEH5YxJ9hZlkokMQqUEna8CN5vZLe7+cLkqBDwInESYtyczU1EjIVst23TCsyUpN3fYsq5rC2bDY9C+I778mCndxyQbMS7ZOotIVSo66Lj7H8zsPOABM/sTYVbR3BmE3N3P621lzMyAY4Am4LUom+0h4MNm9r2slOljgSnALb09VpK6DT/e2ED9yCrukLh1Q87slo9AS568kRHjYMqRWdlkR8DoScnWV0RSo+h+Omb2DuBOCt+ScnevLXJ/vwCeA/5KmOa6ATgbOAH4jLtnpsh+H6GVczOwiD2dQ58H3uPu+ScKL3z8svfTcXcWLl0dO9HSvFlTmT97GlbpLKydm7PGJIvet22ILzt0VOj/kukPM/kIGHeAMslEBpEk++lcQeg/8wFgeT8MQfMnwsRvHyMEsmbgIeBkd789U8jd7zGzkwjp2ncQhsG5FbiwtwEnKQuXrmbx8nW07t7zIL2lLVR58fIws/f5xx+SXIXadnQfk2xLnhnGa4fCvo2hFZO5TbbPwVBT1N8UIiKxSmnp7AQWuPul5a1SMsrd0mne2c5R/353l4CTa9iQGh78ynHlmelvd1vIJMtuwWx6On8m2YQ3RbfIolbMvo3KJBORbpJs6WwE8qQmSa4lKzZQU1P4tlNNjbHkiQ2cedT+fTtYZwe8ujpnTLIV+TPJxh3YNZOs4XAYulff6iAiUoRSgs41wEfM7IfunmeuXsnYtK2VXW2F7/7tautg07bW0nbsHjpXZg8Xs+ExaNseX350Q9cWzOQjYOT40o4pItJPSgk69xFSmR8ws6sIA3N2u6q6+7J+qluqTRg9jOFDa19/hhNn+NBaJoweVnhH217pmkm2/mFoydM9afjY7rNbjmnow6cQEelfpQSdu7P+vZjuowFYtExPmoE5jQ0suG1lwTKdnc6c6VlBoaUpZ3bLR8KUynHq9gq3xTJpylNmhNtmyiQTkSpWStA5t2y1GIDqR9Yxb9ZUFi9fR0t799bOuLrdXHh4K/WPXr0nyGx+Jn5nNXUwqbFrC2bCIcokE5HUKaVz6HXlrMhANH/2NACuWbaaN9W8wCEdf2PGkLU08gzTbD01Kzpixs+2kEmW3YLZtxGG9HAbTkQkBapqErckJTmJW9tNH2XoyhvjV447oGsLpuFwGDaq7HUSEemNRCdxAzCzfYG3EaaX7jbmvLv/tDcVGciG7ncErLwRRk3KetAfZZMpk0xEBpFSOofWAFcCc4kJNhnFDoNTaYlOV71zM+zeBWMml/9YIiJl1NeWTimzY11AGLLmesIYaQZ8iTB99RrCEDaze1OJAW/keAUcERFKCzpnA3e5+1nAb6Nlf3X3/wKOBPaJ3kVERGKVEnSmsifYZAbwqgNw9x3ATwi33kRERGKVEnRaCKNMA2wndASdmLX+Zcows6iIiAwcpWSvPQccBODu7Wb2N8LcNz+L1h8HvNK/1RsYUjeJm4ikQhqvLaVkr30XOMXdD4p+/ipwCfAHQlLBTOAyd/9imerarzSJm4ikVSWvLUn207kM+J2ZDXP3VuBbhNtrHyEM/LkIuLg3lRioqm4SNxEZENJ8bSn6mY67b3D3u6KAg7t3uPu/uvt4d5/g7p9w913lq2q6NO9sZ9GytbHjrgG0tHewaNlamlvaY9eLiMRJ+7WllESCHplZfX/uL81KmcRNRKRYab+2FB10zOz3ZjapwPr3AI/2S60GgLJN4iYig1rary2ltHTeDTxmZnOyF1pwEfC/Je5vQMtM4lZIUZO4iYhkSfu1pZQg8Q7gNeB2M1toZnVmNgW4B/g34DfAW8tQx1Sa09hAZ2fhzMBuk7iJiPQg7deWUhIJHicMc3Mt8DnCWGuPEoLRp939g+6+pRyVTKPMJG4j6uL/IhlRV8u8WVOpH1HdOfUiUl3Sfm0paWoDd28xs48D04D3EEYl+Iy7X1WOyqVdZhK3uFz6uTMPfH29iEgp0nxtKWkSNzM7CLgBmAH8N/BeYD9CJ9FveIpmhEtyaoNuvYanN1TtXyEikh6VuLb0tXNoKSMSfBi4ijDY5zx3vylKkf4x8EHCyAQfcvfqzNPLkeh8OiIiA0SS8+n8DHgSOMLdbwJw92Z3/0fgk8BRwGO9qYSIiAwOpQSd7wAz3f3Z3BXRnDrvQAN+iohIAUUnEvQ0kKe7rzCzt/e9SiIiMlD1W2dOMxsJFD0ns5kda2bXmtkqM9tpZi+a2S1mNj2n3L1m5jGvG/qr7iIikoyCLR0zawPOcvcbop9HA78AvuLuT+QUPxX4KVC4q+weHwf2Bi4HngL2BS4E/mJmx7j7A1ll1wBn5Wz/apHHEREZkNI4n05Pt9eG0LU1NBQ4CfhePxz7U+6+MXuBmf0OWAd8ATgta9XOnCAkIjJo5ZtPZ8FtK6t+rq6SOof2p9yAEy1rMrM1hL4/IiISY1DMp5MEM5sANAIrclYdYmZbzGy3ma0xs6+aWXW3IUVEyiDt8+lUrKWTy0JbcBEhEF6WtWo5YRSEp4FRwCmEERCOJDxHyre/nnp9au4fEUmdUubTOfOo/ROqVfGqJugQ+gGdApzr7k9lFrr7RTnlfmNmrwBfNrP3uvt9SVZSRKSS0j6fTjFBZ07W5G0jCYN8/pOZ5U5jcGRvK2Fm3wTOBz7r7tcWscl1wJeBdwGxQaenIRqilpBaOyKSKpn5dFoKBJ5qnk+nmKDzoeiV7WN5ypY84KeZXUIIIBe6+/eL3CzzLKqzYCkRkQFmTmMDC25bWbBMNc+n01PQ+ftyHtzMLgYuAi5y9++UsGmmz47SqEVkUMnMp7N4+brYZIIRdbXMnXlg1Y5kXzDouPsfynVgMzsfWECYcfRuM3tn1upWd3/EzGYCXwJuBp4D9gI+AJwL3OTufyxX/UREqtWgmU+nXw9sdi9wdJ7Vz7n7AWb2RuAK4HBgH8LttFWEZzo/cPfCT9MKH19TG4hIqg3o+XQGGgUdEZHSJTmfjoiISJ8o6IiISGIUdEREJDEKOiIikhgFHRERSYyCjoiIJEZBR0REEqOgIyIiiVHQERGRxCjoiIhIYhR0REQkMQo6IiKSGAUdERFJTDEzh0ofdRt+vLGB+pHVOcGSiEg5KeiUkbuzcOnqbhMtLbhtJfNmTWX+7GmYWaWrKSKSGAWdMlq4dDWLl6+jdXfn68ta2sK8c4uXrwPg/OMPqUjdREQqQc90yqR5ZzuLlq2NncMcoKW9g0XL1tLc0p5wzUREKkdBp0yWrNhATU3hW2c1NcaSJzYkVCMRkcpT0CmTTdta2dUW38rJ2NXWwaZtrQnVSESk8hR0ymTC6GEMH1pbsMzwobVMGD0soRqJiFSegk6ZzGlsoLPTC5bp7HTmTG9IqEYiIpWnoFMm9SPrmDdrKiPq4ls7I+pqmTdrKvUj1F9HRAYPpUyX0fzZ0wC69dPp7HTmzjzw9fUiIoOFuRe+BTRQmVlTfX19fVNTU9mP1W1EgukNauGISCqNHTuW5ubmZncf25vt1dJJQP3IOs48av9KV0NEpOL0TEdERBKjoCMiIolR0BERkcQo6IiISGIqFnTM7Fgzu9bMVpnZTjN70cxuMbPpMWVnm9kDZtZiZhvN7Edm1qvMCRERqZxKtnQ+DuwPXA6cCMyPfv6Lmb0zU8jMjgGWAC8A7wcuAE4G7jAztdRERFKkYv10zGyiu2/MWTYWWAfc4+6nRcseBOqAI929M1o2G/gdcIa7/7KXx0+sn46IyEDR1346FWsp5AacaFkTsAbYD8DMpgBvB36WCThRuaXAeuC0ZGorIiL9oao6h5rZBKARuD5a1Bi9r4gp/kTW+rh99dSEqS+5giIi0idV80zEzAxYRKjTZdHivaP3zTGbbM5aLyIiKVBNLZ3vAKcA57r9pFZcAAAPvElEQVT7Uznr8j14yvtAqqf7jVFLSK0dEZEEVUVLx8y+CZwPfNbdr81a9Vr0HteiGU98C0hERKpUxVs6ZnYJ8GXgQnf/fs7qldF7IyFbLdt04P4yV69fdBtlurGB+pEaZVpEBp+KTm1gZhcDC4CL3P0becr8hdAie3tWyvSxwN3Ame5+Qy+PXfaUaXdn4dLVsfPpzJs1lfmzpxEeZYmIpENqpzYws/MJAec3wN3ZHUKBVnd/JPr3FwmtnOvNbBEwGfg28GfgpuRqXLqFS1ezePk6Wne/nu1NS1sHAIuXrwPg/OMPqUjdREQqoZLPdN4fvZ8E/Cnn9atMIXe/JypzAHAHsDB6P9HdOxKsb0mad7azaNlaWtrjq9jS3sGiZWtpbmlPuGYiIpVTsZaOux9TQtk7gTvLV5v+t2TFBmpqCt86q6kxljyxQRO8icigURXZawPRpm2t7Gor3BDb1dbBpm2tCdVIRKTyFHTKZMLoYQwfWluwzPChtUwYPSyhGomIVJ6CTpnMaWygs7NwZmBnpzNnekNCNRIRqTwFnTKpH1nHvFlTGVEX39oZUVfLvFlTqR+h/joiMnhUvHPoQDZ/9jSA2H46c2ce+Pp6EZHBoqKdQyspyfl0uo1IML1BLRwRSaXUdg4dTOpH1iktWkQEPdMREZEEKeiIiEhiBvMznU7A6us1pY6ISLGam5sB3N171WgZzEFnN6GltzWhQ2aiW3NCx0sDnZN4Oi/xdF7iJX1exgCd7t6rnIBBG3SSFs1U2uOMpoOJzkk8nZd4Oi/x0nZe9ExHREQSo6AjIiKJUdAREZHEKOiIiEhiFHRERCQxCjoiIpIYBR0REUmM+umIiEhi1NIREZHEKOiIiEhiFHRERCQxCjpFMrNRZvZ9M9tgZi1m9pCZnVzktgeZ2a1m1mxm28xsiZkdmqfsv5rZajNrNbNnzOxCM6va31Nvz4uZzTWz28zsuWi7NdF+JsSU9Tyvj5fnU/VdH87Lgjyf9eU85QfL9+XZAt+Dp3PKpur7Ymb7mdkVZnafmW2P6npMCdun6/ri7noV8QKWAq8B5wHvA34KdABzethuIvAS8ChwCnAS8ACwCdgvp+xXo31eAhwDfBloB/6j0p+/DOdlPfBz4EPA0cAngJeBZ4GxOWUduAF4Z85rYqU/fxnOy4Lo8x6X81lnxJQdTN+XI2J+/x+NztV/5JRN1fcl+t1tBO4Efh3V/5git03d9aXiJzwNL2BO9EU4NWuZAfcBT/Ww7aVACzA5a9nehCkV/jNnWQtwRc7234y+GPv19XNU2XnpdgGIgo8Dn8lZ7sD3Kv15EzovmaAztodyg+r7kmd/34/2Ny3l35earH+fUmLQSd31pWqb4VXmVMJcFb/OLPDwG7sOeFO+pmzWtkvd/aWsbV8Dbgc+mFXuBGB4tM9s1wJDgKJu5SWs1+fF3TfGLP5L9L5ff1ayAvryfSnWoPq+5DKzoYRW8n3uvrq/K5okd+/sw+apu74o6BSnEXgy5svxeNb6bsxsBHAQsCJm9ePARDObmLUPB1ZmF3L3NYS/UGKPUWG9Oi8FvC96jztfZ0XPAHaZ2Z/N7PQS952k/jgvT5lZR/Ts4+qs70n2MQbz9+UUwl/v1+RZn6bvS6+k9fqioFOcvYHNMcs3Z62PM45w+6CYbfcGdrp7a0zZLQWOUUm9PS/dmNl4wu2SNcCNOat/AXwaOB44i/Cf5Jdm9tlSK5yQvpyXZwj32s8FZgNXAWcAD5jZuJxjDNrvC/AvwHa6f1cgfd+X3krl9aVX040OUoWGbuhpWIdit+3LMSqlz3U2s5HArcB4YFbufwx3/0hO+f8B7gW+YWaL3L2lpBono1fnxd1/lrPoHjN7APgd8CngG309RoX1x/dlP0JA/om77+i2k3R+X/oiVdcXtXSK8xrxfwmMj97j/tKA8BeEF7nta8BeZjYspuy4AseopN6el9dFtwhuI2QnzXH3x3vYJHMP/OfAKKrzNlKfz0s2d18KbADelXOMQfd9iZxDuHblu7XWRQq+L72VyuuLgk5xVgJvjslnnx69x91TJfqLai3xX/TpwKasB+orCU3lw7ILmdkbgRH5jlFhvTovGWY2nPBQ+V3ASe5+fwnHzhyzLw9hy6VP5yWPGrp+1kH3fQEwMyMEnacH0PelV9J6fVHQKc6vgLHA+3OWnwWscvcne9h2tplNyiyInl+8H7glq9xvgVbg/+Zsfzawm5CNUm16fV6iv7huBWYCH3D3PxR70Oii9WFgGzkPRqtEX74v3ZjZ8cC+hP4XGYPq+5LlaMLD86JaOZCK70tfpO/6Uukc9TS8CH8h3AO8SniA+feEVMNO4P1Z5e4lygLNWrYvodPjw8AHgH8A/kRo7u6fU/ZiwhdgAeE/15eANuA7lT4HZTgvtxNuDfwb3TvxHZRV7gLgauBMQoe2M4A/RNt+stLnoAzn5RHg84Q+LbOj78Q2QoJFbqfZQfN9yVr3U0K/kn3zrE/d9yWq9z9Gr29Hdb04+vnEHr4vqbu+VPxkp+UFjAF+GP2Cd0W/5FNyysT+ZwEOJtxG2krIuPktcFhMOQM+F11gWoF1wP8jq/NYtb16e16i/1j5XtdmlXs/sJzQw7odaAJ+n32RqsZXH87L9dHvf0d0QXgGuBwYP5i/L9Hy0dF5+XWB/af1+5Lv/8KzRZyXVF1fNJ+OiIgkRs90REQkMQo6IiKSGAUdERFJjIKOiIgkRkFHREQSo6AjIiKJUdAREZHEKOhI1TGzY3Lmtu8wsy1mtsLMrjOzE6IxuFLBzMaZ2dfM7C9m1mRmbWb2opndbGYf7O1nMbPPmdk5/VxdkbJS51CpOmZ2DPC/hN75Swg9qUcDhxAm79ofuBv4J3dvqlA1i2JmRxF6i08kjKb9B0LP8SmEoW7eCXzK3a/qxb6fJfRYP6a/6itSbppPR6rZw+7+8+wFZjafMC/8fEJQOrHQDsysDqh1911lq2X+Y08ijDE3HDja3e/LKfJ1M/s/hKHlBwUzqwWGufvOStdFKkO31yRV3L3D3c8H7gNOMLP3ZtaZ2YLodtxhZrbQzF4kjO/1TjM7IFq3IHefWdsdkLP8aDP7UzTt8ctmdkW079j9xPgCoYXzxZiAk/k8d7n7DVnH/Gczu83MnjezVjN71cxuNbO35NTNgb8Djs65FXlAVpm3mdmvon20mtkqM/uKmXX7Y9PMTjOzx6LpnZ83s4vN7Lhon+fklN3HzK40sxeiW4UvRD/vnVPunGj748zsIjN7hvD7OD061vMx0xxgZqdH2+WOiCwDgFo6klY/Bt5LGFU394L+C8IUxd8lDJq4odSdR8Hsd4SJsv6DMHDk6cB7StjNaYRBO68rYZtPEybUWkQYFPMgYB7wRzOb4WFOewhD1F9OGLH5m1nbb4rqP4cw7P3fCOdhM2HeokuAtwL/lNnAzP6Z0Gp8hjDq927CkPe5UxBgZvXA/cAbCdMLPEyYgO8TwPvM7Ch335az2WVAHWH0563AqujfPyCMpH1XTvl/AZqB/yl0oiSlKj26ql565b4IQ9I7cEGBMjOiMjdnLVsQLbsXGJJT/oBo3YKYfWW2OyBr2YOEv8qnZi2rA/6Ybz85+xwdlXu8xM++V8yyNxNGBb4qZ/mzwL0x5YcTAtaymPPw+ahex0Q/DwHWA68A47LKjSJMEObAOVnLv0nMNAGEqbQd+HrWsnOiZauAkTnl6wkjRt+Ys/wNQEfuZ9Vr4Lx0e03Samv0PiZm3ffcfXdvd2xm+wJvJwyhvzaz3N3bgSuK3E2mXlsLlsrh7juiOpiZjTGzfQitl1XAO4rczWzCPCs/AcZGt8P2ifa1JCpzfPR+JDCZMJ3Elqx6bAf+K2bfp0b1WZSz/EeEVtepMdv8p+c8w3H3ZuAm4ANRvTLOJdz2/3GPn1JSSUFH0qrQRX11H/d9YPS+KmZd3LI4mXqNLuXAZnaEmf2GMHFbM+ECv4kw/XCxCQdvjt6vydo+83o6Wrdv9F7qZz2QMMtnl6Ae/bwKmBqzTb7fxyJgKPAReH0q6nOBR939r3m2kZTTMx1Jq8yD9bgLY1xmVKG+Abn/D/rcB8jdt5nZc8CbzGyEh/nsCzKz/Qm3xLYCXyd8th2Eun+PcMurGJn6fwF4NE+Zl3LKllNsppq7329mK4DzCJ/vWMJt0E8nUCepEAUdSavzovc7iiy/OXofH7Mu96/zzC21Q2LKxi3L5xbCM5T/S/fbUXFOJQSWk939f7NXRJlhrTnl8wXSTLLBDne/u4djrovei/2sa4FDzGxIdmsnyoibxp5zV6yrgSui/kznEZ6j/aLEfUiK6PaapIqZ1ZrZZYTMtSXu/sditvOQUfUyIcPq9b/uzWwqocNpdtlXgIcIzxumZpWtAz5bQnUvJdzSutTM3pXn8xxvZmdEP3ZkFueU+SgwKWbz7cQH0buAjcCXzKzbejMbYWaZ234PEbL7zjGzcVllRgEfj9n3rcAEYG7O8o9Gy38Vs00hPyMEmi8Qgu7NXuUdfqVv1NKRajbDzD4S/Tt7RIK/I6Qzf6jE/f0Q+AbwWzO7lfAA/ePACkLiQLYLgKXA/WZ2FeH5yumEZxBQ+HZdKOD+spmdRBiR4L7omJnbZ5OBEwjB8xPRJr8l3Ir6mZn9kJCu/R7CyAXP0P3/6wPAeWb2deApoBO43d13mNlZhACxysyuIaROjwXeBHyQcIG/1913m9kFhNbFg2b2Y0LK9DnAa4RnONmf9VJCuvWVZjYDeISQMn0e4XbgpT2dl5xztMXM/ofouQ6wuJTtJYUqnT6nl165L/akTGdeHYSL/kpCn5cT8my3gJzU55z1QwgXxQ2Ev64fJvRFid0OeB/hwr6LkFJ8BSGDzIELS/g844GLCa2KZkLfnRcJ/VBOzik7i9DvaBuhb9AdQCMhDfzZnLITgZsJtw47cz9DtN3PCSnRbdFnuB+4CBifs6/TgccJt/Cej+p7arTP03PKTgCuij5De/R+JbBPTrlzyErPLnB+Zkbl1hANzaXXwH1p7DWREpjZaYRgcaZnjSQwEJnZ+YSOne9y9wfKeJyjgD8DX3b3b5XrOFIdFHREYkTPfYZ51pht0TOde4GjgDe4+8sVql6/MrOhQIe7d2QtG0Vo+YwBJrt7WxmP/1PgDGD/gXJOJT890xGJNwx4zsx+QXhWsTfwz4RU7W8PsIvjVMJzrhsI2WwNhGFwDgQ+UY6AY2Z7EW5tHkZ4nrNogJ1TyUMtHZEYFkZDvho4mnARNkLwWeS9mIagmkXp2D8kJC1MJCQSPAFc7u43lumYBxAC3HZCAsVcdy9p9AZJJwUdERFJjPrpiIhIYhR0REQkMQo6IiKSGAUdERFJjIKOiIgkRkFHREQS8/8B/Rxw/p4WYR8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(drugged, exam, 'o')\n",
    "plt.plot(drugged, slope * drugged + intercept, '-')\n",
    "plt.xlim(-0.1, 1.1)\n",
    "plt.xlabel('Drug Category')\n",
    "plt.ylabel('Exam Score')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Well that's interesting. But can we get a $p$-value out of this? We can, by seeing if the slope is necessary. Let's take the null hypothesis to be that $\\beta_1 = 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.032\n"
     ]
    }
   ],
   "source": [
    "s2_e = np.sum( (exam - slope * drugged - intercept)**2 ) / dof\n",
    "s2_b = s2_e / np.sum( (drugged - np.mean(drugged)) ** 2 )\n",
    "\n",
    "slope_se = np.sqrt(s2_b)\n",
    "T = slope / slope_se\n",
    "\n",
    "#The null hypothesis is that you have no slope, so DOF = N - 1\n",
    "print('{:.2}'.format(2 * (1 - boop.t.cdf(T, len(drugged) - 1))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We get a $p$-value of $0.032$, which is quite close to what the sum of ranks test gave. However, we can now do more dimensions than the sum of ranks test."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Multiple Categories\n",
    "---\n",
    "\n",
    "Now let's say that I'm studying plant growth and I have two categories: I fertilized the plant and I put the plant in direct sunlight. Let's turn that into two variables: $\\delta_s$ for sunlight and $\\delta_f$ for fertilizer. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now you might believe that there is an *interaction* between these two *factors*. That is, having sunlight and fertilizer has an effect beyond the sum of the two individually. I can write this as $\\delta_{fs}$. This category is 1 when both fertilizer and direct sunlight is provided. It turns out:\n",
    "\n",
    "$$\\delta_{fs} = \\delta_s \\times \\delta_f $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Our model equation will then be:\n",
    "    \n",
    "$$ g = \\beta_0 + \\beta_s \\delta_s + \\beta_f \\delta_f + \\beta_{sf} \\delta_s \\delta_f + \\epsilon$$\n",
    "\n",
    "and here's some data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "growth = [13.4,11.4, 11.9, 13.1, 14.8, 12.0, 14.3, 13.2, 13.4, 12.2, 1.9, 3.1, 4.4, 3.9, 4.2, 3.2, 2.1, 2.1, 3.4, 4.2, 5.6, 4.5, 4.7, 6.9, 5.1, 3.2, 2.7, 5.2, 5.4, 4.9, 3.5, 2.3, 3.4, 5.4, 4.1, 4.0, 3.2, 4.1, 3.3, 3.4]\n",
    "fertilizer = [1.0,1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
    "sunlight = [1.0,1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.   1.  13.4]\n",
      " [ 1.   1.  11.4]\n",
      " [ 1.   1.  11.9]\n",
      " [ 1.   1.  13.1]\n",
      " [ 1.   1.  14.8]\n",
      " [ 1.   1.  12. ]\n",
      " [ 1.   1.  14.3]\n",
      " [ 1.   1.  13.2]\n",
      " [ 1.   1.  13.4]\n",
      " [ 1.   1.  12.2]\n",
      " [ 1.   0.   1.9]\n",
      " [ 1.   0.   3.1]\n",
      " [ 1.   0.   4.4]\n",
      " [ 1.   0.   3.9]\n",
      " [ 1.   0.   4.2]\n",
      " [ 1.   0.   3.2]\n",
      " [ 1.   0.   2.1]\n",
      " [ 1.   0.   2.1]\n",
      " [ 1.   0.   3.4]\n",
      " [ 1.   0.   4.2]\n",
      " [ 0.   1.   5.6]\n",
      " [ 0.   1.   4.5]\n",
      " [ 0.   1.   4.7]\n",
      " [ 0.   1.   6.9]\n",
      " [ 0.   1.   5.1]\n",
      " [ 0.   1.   3.2]\n",
      " [ 0.   1.   2.7]\n",
      " [ 0.   1.   5.2]\n",
      " [ 0.   1.   5.4]\n",
      " [ 0.   1.   4.9]\n",
      " [ 0.   0.   3.5]\n",
      " [ 0.   0.   2.3]\n",
      " [ 0.   0.   3.4]\n",
      " [ 0.   0.   5.4]\n",
      " [ 0.   0.   4.1]\n",
      " [ 0.   0.   4. ]\n",
      " [ 0.   0.   3.2]\n",
      " [ 0.   0.   4.1]\n",
      " [ 0.   0.   3.3]\n",
      " [ 0.   0.   3.4]]\n"
     ]
    }
   ],
   "source": [
    "print(np.column_stack((fertilizer, sunlight, growth)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "You can see we have 10 examples of each possible combination of categories. Let's now regress it using multidimensional least squares. Our columns will be $[1, \\delta_s, \\delta_f, \\delta_s \\delta_f] $, so we'll be doing 4-dimensional regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 1. 1. 1.]\n",
      " [1. 1. 1. 1.]\n",
      " [1. 1. 1. 1.]\n",
      " [1. 1. 1. 1.]\n",
      " [1. 1. 1. 1.]\n",
      " [1. 1. 1. 1.]\n",
      " [1. 1. 1. 1.]\n",
      " [1. 1. 1. 1.]\n",
      " [1. 1. 1. 1.]\n",
      " [1. 1. 1. 1.]\n",
      " [1. 0. 1. 0.]\n",
      " [1. 0. 1. 0.]\n",
      " [1. 0. 1. 0.]\n",
      " [1. 0. 1. 0.]\n",
      " [1. 0. 1. 0.]\n",
      " [1. 0. 1. 0.]\n",
      " [1. 0. 1. 0.]\n",
      " [1. 0. 1. 0.]\n",
      " [1. 0. 1. 0.]\n",
      " [1. 0. 1. 0.]\n",
      " [1. 1. 0. 0.]\n",
      " [1. 1. 0. 0.]\n",
      " [1. 1. 0. 0.]\n",
      " [1. 1. 0. 0.]\n",
      " [1. 1. 0. 0.]\n",
      " [1. 1. 0. 0.]\n",
      " [1. 1. 0. 0.]\n",
      " [1. 1. 0. 0.]\n",
      " [1. 1. 0. 0.]\n",
      " [1. 1. 0. 0.]\n",
      " [1. 0. 0. 0.]\n",
      " [1. 0. 0. 0.]\n",
      " [1. 0. 0. 0.]\n",
      " [1. 0. 0. 0.]\n",
      " [1. 0. 0. 0.]\n",
      " [1. 0. 0. 0.]\n",
      " [1. 0. 0. 0.]\n",
      " [1. 0. 0. 0.]\n",
      " [1. 0. 0. 0.]\n",
      " [1. 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "N = len(growth)\n",
    "dof = N - 4\n",
    "\n",
    "x_mat = np.column_stack( (np.ones(N), sunlight, fertilizer, np.array(sunlight) * np.array(fertilizer)) )\n",
    "print(x_mat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 3.67  1.15 -0.42  8.57]\n"
     ]
    }
   ],
   "source": [
    "import scipy.linalg as linalg\n",
    "\n",
    "#I'll use the lstsq convienence function\n",
    "#It gives back a bunch of stuff I don't want though\n",
    "#I'll throw it away by assigning it all the an underscore\n",
    "beta,*_ = linalg.lstsq(x_mat, growth)\n",
    "\n",
    "print(beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "That tells us about what you'd expect: the combination is quite important! Let's get individual confidence intervals on the factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta_0: 3.67 +/- 0.65\n",
      "beta_s: 1.15 +/- 0.92\n",
      "beta_f: -0.42 +/- 0.92\n",
      "beta_sf: 8.57 +/- 1.3\n"
     ]
    }
   ],
   "source": [
    "s2_e = np.sum( (growth - x_mat.dot(beta))**2 ) / dof\n",
    "s2_b = linalg.inv(x_mat.transpose().dot(x_mat)) * s2_e\n",
    "s_b = np.sqrt(np.diag(s2_b))\n",
    "\n",
    "names = ['beta_0', 'beta_s', 'beta_f', 'beta_sf']\n",
    "for bi, si, ni in zip(beta, s_b, names):\n",
    "    print('{}: {:.3} +/- {:.2}'.format(ni, bi, si * boop.t.ppf(0.975, dof)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We could go farther and justify the existence of the coefficients. As you can see from the confidence intervals, the $\\beta_s$ and/or $\\beta_f$ may be unnecessary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Some notes on our regression\n",
    "---\n",
    "\n",
    "What we've learned about so far is part of a large field: the design and analysis of experiments with statistics. We've used simple analysis, hypothesis tests on existence of coefficeints, but that is not the most robust measure. For example, if two coefficients are unneeded, do you remove both or one at a time? A better way is to use an *Analysis of Variance* (ANOVA) and *F-test*, a hypothesis test for comparing models. You can find more information on this in your Langley statistics book. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Unordered categorical data\n",
    "===\n",
    "\n",
    "When working with categorical data that has more than one category, but they are unordered, you must create a series of dummy variables. For example, if your category is color from three possibilities: `red`, `blue`, `orange`, it's *incorrect* to assign a number to each color. This would induce an ordering and imply that `red` is greater than `blue` or vice-versa. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "You should instead create dummy variables for each value in the category:\n",
    "\n",
    "$$\n",
    "[\\delta_{\\textrm{red}}, \\delta_{\\textrm{blue}}, \\delta_{\\textrm{orange}}]\n",
    "$$"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}