{ "metadata": { "name": "", "signature": "sha256:7c9ecc18673a4d35de67d0c3e73cb58f642e7bf5d3a33d1d9f32c351473a681b" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "These are some notes on simple regression, multiple regression, and the general linear model.\n", "\n", "For a more background, please see [The general linear model and fMRI: Does love last forever?](http://matthew.dynevor.org/_downloads/does_glm_love.pdf).\n", "\n", "This page starts by setting up a simple regression. Then we show how the simple regression gets expressed in a *design matrix*. Once we have that, it's easy to extend simple regression to multiple regression. By adding some specially formed regressors, we can also express group membership, and therefore do analysis of variance. This last step is where multiple regression becomes the general linear model." ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "The example regression problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's imagine that we have measured scores for a \"psychopathy\" personality trait in 12 students. We also have some other information about these students. For example, we measured how much sweat each student had on their palms, and we call this a \"clammy\" score. We first try and work out whether the \"clammy\" score predicts the \"psychopathy\" score. We'll do this with simple linear regression." ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Simple linear regression" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We first need to get our environment set up to run the code and plots we need." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Set up notebook to show data plots in the notebook\n", "%matplotlib inline" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "code", "collapsed": false, "input": [ "# Import numerical and plotting libraries\n", "import numpy as np\n", "import numpy.linalg as npl\n", "import matplotlib.pyplot as plt" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "code", "collapsed": false, "input": [ "# Get t and gamma distribution code from scipy library\n", "from scipy.stats import t as t_dist, gamma" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "code", "collapsed": false, "input": [ "# Get symbolic mathematics routines\n", "from sympy import (Matrix, Symbol, Eq, MatAdd, MatMul, ones,\n", " init_printing)\n", "# Set up notebook to output symbolic math expressions nicely\n", "init_printing(use_latex=True)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here are our scores of \"psychopathy\" from the 12 students:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "psychopathy = [11.416, 4.514, 12.204, 14.835,\n", " 8.416, 6.563, 17.343, 13.02,\n", " 15.19 , 11.902, 22.721, 22.324]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "These are the skin-conductance scores to get a measure of clamminess for the handshakes of each student:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "clammy = [0.389, 0.2 , 0.241, 0.463,\n", " 4.585, 1.097, 1.642, 4.972,\n", " 7.957, 5.585, 5.527, 6.964]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We happen to believe that there is some relationship between `clammy` and `psychopathy`. Plotting them together we get:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "plt.plot(clammy, psychopathy, '+')\n", "plt.xlabel('Clamminess of handshake')\n", "plt.ylabel('Psychopathy score')" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 7, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEPCAYAAACukxSbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtUlHXiP/D3I0iZoKghg1HBAS/AzDAj47DmbSAVTc1L\nhrcSFahVq9XYyE5nv6K1pttpW6zdcivtrmFuppaeshxEzTWVS2qZIbMiikcDRQRDhs/vD36OcZ0Z\nmZkHfN6vc+YcZxieeSPy5uNnnufzkYQQAkREpAid5A5ARESew9InIlIQlj4RkYKw9ImIFISlT0Sk\nICx9IiIFcVvpFxcXIy4uDlFRUVCr1Vi9ejUAICMjA8HBwdDr9dDr9dixY4e7IhARUSOSu87TLy0t\nRWlpKXQ6HSorKxETE4PNmzcjKysLfn5+ePrpp93xskRE1Apvdx1YpVJBpVIBAHx9fREREYGSkhIA\nAK8HIyKSh0fm9C0WC3Jzc/GHP/wBAPDaa68hOjoaycnJuHjxoiciEBERPFD6lZWVmDp1KjIzM+Hr\n64v58+ejqKgIeXl5CAoKQlpamrsjEBHRdcKNampqxOjRo8Wrr77a7MeLioqEWq1u8nhYWJgAwBtv\nvPHGmxO3sLAwu73stpG+EALJycmIjIzEokWLbI+fPXvW9ufPPvsMGo2myecWFhZCCNHub0uXLpU9\nA3MyJ3My4/VbYWGh3W522xu5e/fuxYcffgitVgu9Xg8AWLFiBdavX4+8vDxIkoTQ0FCsWbPGXRGI\niKgRt5X+0KFDUVdX1+TxsWPHuusliYjIDl6R2wYmk0nuCA5hTtdiTtfqCDk7QkZHue3irLaQJAnt\nMBYRUbvmSHdypE9EpCAsfSIiBWHpExEpCEufiEhBWPpERArC0iciUhCWPhGRgrD0iYgUhKVPRKQg\nLH0iIgVh6RMRKQhLn4icYjbLnYDagqVPRE5h6XdsLH0iIgVx2yYqRHTrMJtvjPCXLbvxuMlUf7uV\nmc231tfI0iciuxqXe0aGTEFkcKuVPqd3iIgUhCN9InLKrTTqbcmtPJ3F7RKJiFqRkdFxprO4XSIR\nETXA0iciakVHn85pjNM7RES3CE7vEBFRAyx9IiIFYekTESkIS5+ISEFY+kRECsLSJyJSEJY+EZGC\nsPSJiBSEpU9EpCAsfSIiBWHpExEpCEufiEhB3Fb6xcXFiIuLQ1RUFNRqNVavXg0AKCsrw6hRo9Cv\nXz+MHj0aFy9edFcEIiJqxG2rbJaWlqK0tBQ6nQ6VlZWIiYnB5s2bsW7dOtx5551IT0/HqlWrUF5e\njpUrVzYMxVU2iYicJusqmyqVCjqdDgDg6+uLiIgIlJSUYMuWLUhKSgIAJCUlYfPmze6KQEREjXhk\nPX2LxYIRI0bgyJEjuOeee1BeXg4AEEKgZ8+etvu2UBzpExE5zZHudPvG6JWVlXjooYeQmZkJPz+/\nBh+TJAmSJDX7eRm/25TSZDLBdKttX0NE1EZmsxnm6zu4O8itI/1r165h/PjxGDt2LBYtWgQAGDBg\nAMxmM1QqFc6ePYu4uDj89NNPDUNxpE9E5DRZ5/SFEEhOTkZkZKSt8AHgwQcfxHvvvQcAeO+99zBp\n0iR3RSAiokbcNtLfs2cPhg8fDq1Wa5vCeemll2A0GpGYmIhTp04hJCQEWVlZ8Pf3bxiKI30iIqc5\n0p3cGJ2I6BbBjdGJiKgBlj4RkYKw9ImIFISlT0SkICx9IiIFYekTESkIS5+ISEFY+kRECsLSJyJS\nEJY+EZGCsPSJiBSEpU9EpCAsfSIiBWHpExEpCEufiEhBWPpERArC0u8gnNz7mIioWXZL/8qVK3jh\nhReQmpoKADhx4gS2bdvm9mDUEEufiFzBbunPnTsXPj4+2LdvHwCgT58+eP75590ejIiIXM/b3hMK\nCwuRlZWFDRs2AAC6du3q9lBUz2y+McJftuzG4yZT/Y2IyFl2S/+2225DdXW17X5hYSFuu+02t4ai\neo3LPSNDpiBEdMuwW/oZGRkYM2YMTp8+jZkzZ2Lv3r149913PRCNiIhcrdXSr6urQ3l5OTZt2oT9\n+/cDADIzMxEQEOCRcHQDp3OIyBUkIYRo7QkxMTE4dOiQp/IAACRJgp1YRETUiCPdabf0lyxZgjvv\nvBPTpk1r8CZuz549XZOyuVAsfSIip7mk9ENCQiBJUpMDnzx5su0JWwrF0icicppLSl8OLH0iIuc5\n0p12z96pqanBG2+8gd27d0OSJIwYMQJ//OMf0blzZ5cFJSIiz7A70k9OTkZtbS2SkpIghMAHH3wA\nb29vvP322+4LxZE+EZHTXDK9o9VqUVBQYPcxV2LpExE5z5HutLv2jre3N3755Rfb/cLCQnh7250V\nIiKidshue7/88suIj49HaGgoAMBisWDdunVuD0ZERK7n0Nk7V69exfHjxyFJEvr164fbb7/dvaE4\nvUNE5DSXTO+8/vrrqK6uRnR0NLRaLaqrq/Gvf/3LZSGJiMhz7I70o6OjkZ+f3+AxnU6HvLw894Xi\nSJ+IyGkuGenX1dWhrq7Odt9qteLatWsOBZg3bx4CAwOh0Whsj2VkZCA4OBh6vR56vR47duxw6FhE\nRNQ8Z3bWs1v6CQkJmD59Or755hvs3LkT06dPx5gxYxw6+Ny5c5uUuiRJePrpp5Gbm4vc3FyHj0VE\nRM1zpvTtnr2zatUq/Pvf/8Ybb7wBABg1ahRSUlIcOviwYcNgsViaPM6pGyIiedgtfS8vL8yfPx/z\n589HWVkZiouL4eXl1aYXfe211/D+++/DYDDglVdegb+/f5uOZ4/ZzPXoiejW0tJ2qvbYnd4ZMWIE\nKioqUFZWhpiYGKSmpmLx4sU3lxLA/PnzUVRUhLy8PAQFBSEtLe2mj+UoZ/7rQ0TUEZhM9VuoZmQA\nS5c6vp2q3ZH+pUuX0K1bN7z99tuYPXs2li1b1uCNWWf17t3b9ueUlBRMmDCh2edl/O4rMJlMMHGo\nTkTUgNls/v83F5a+1WrF2bNnkZWVhRdffBEAmqyv74yzZ88iKCgIAPDZZ5+1+Asko427gLf0X5/G\nm40TEXVU1wfE16ewlzkwz2O39P/v//4PCQkJGDJkCIxGIwoLC9G3b1+HAs2YMQPZ2dm4cOEC7r77\nbixbtgxmsxl5eXmQJAmhoaFYs2aNQ8dyVuNyb+PvECKidsuZgawiNlG5Pu9FRHQrc8nFWbcCTucQ\nEdVTxEifiEgJXDLSt1qtLgtERETyslv6ffv2xTPPPINjx455Ig8REbmR3dLPy8tD3759kZKSgtjY\nWKxZswYVFRWeyEZERC7m1Jy+2WzGrFmzUF5ejocffhh/+ctfEB4e7vpQnNMnInKaS+b0a2tr8fnn\nn2PSpElYtGgR0tLScPLkSUyYMAEPPPCAy8ISEZH72b04q1+/fjCZTEhPT8d9991ne3zq1KnIzs52\nazgiInItu9M7ly9fhp+fn6fyAOD0DhHRzXCkO+2O9Kurq7F69WpYLBbU1tbaDrx27VrXpCQiIo+x\nW/oTJ07E8OHDMWrUKHTqVP8WQFsWXCMiIvnYnd5x9ybozeH0DlHbcOMgZXLJ2Tvjx4/HF1984bJQ\nROR+3DiIWtLi9I6vr69tGmfFihXw8fFB586dAdT/NuEFWkREHU+LpV9ZWenJHETURtw4iBxh943c\n+++/H998843dx4hIXtw4iBzRYulXV1ejqqoK58+fR1lZme3xiooKlJSUeCQcERG5Voulv2bNGmRm\nZuLMmTOIiYmxPe7n54cnnnjCI+GI6OZwOodaYveUzdWrV+Opp57yVB4Azp2yyVPTiIjqOdKdDq2y\neeTIERw7dgxXr161PTZ79uy2J2wplBOlz/1viYjquWQZhoyMDGRnZ+Po0aMYN24ctm/fjqFDh7q1\n9ImIyD3slv6nn36K/Px8DBw4EOvWrcO5c+cwa9YsT2RrEU9NI5Ifp1Y7Jrul36VLF3h5ecHb2xuX\nLl1C7969UVxc7IlsLeKpaUTyY+l3THZLf9CgQSgvL0dqaioMBgO6du3aYF19IiLqOJzaLtFisaCi\nogJardadmXj2DlE71XhqdenS+j9zarV9cMnZO0II/Oc//8GePXsgSRKGDRuGyZMnuzRok1AuWGWT\nvwyI3ItnzrU/Llllc8GCBVizZg20Wi3UajXWrFmDBQsWuCyku3CVQSKipuzO6e/atQvHjh2zbaAy\nZ84cREZGuj0YEbVv/J90x2S39MPDw3Hq1CmEhIQAAE6dOoXw8HB357opPJWTyHP4M9Ux2S39iooK\nREREwGg0QpIkHDhwAIMGDcKECRMgSRK2bNniiZwO4amcRESts1v6y5cvB3BjX9zfv0nAvXKJiDoW\nh07ZLC0txffffw9JkmA0GtG7d2/3huLZO0RETnPJ2TtZWVmIjY3Fxo0bkZWVBaPRiI0bN7ospLuw\n8ImImrI70tdqtdi5c6dtdH/+/Hncf//9KCgocF8oF4z0iYiUxiUjfSEEAgICbPd79erFQiYi6qDs\nvpE7ZswYJCQkYObMmRBC4JNPPsHYsWM9kY2IiFzMrcswzJs3D1988QV69+6NH374AQBQVlaGadOm\n4X//+x9CQkKQlZUFf3//hqE4vUNE5DSX7Zx1s3JycuDr64vZs2fbSj89PR133nkn0tPTsWrVKpSX\nl2PlypVOByciooZcMqe/adMm9O3bF926dYOfnx/8/PzQrVs3hwIMGzYMPXr0aPDYli1bkJSUBABI\nSkrC5s2bHToWERG1nd05/fT0dGzbtg0REREuecFz584hMDAQABAYGIhz58655LhERGSf3ZG+SqVy\nWeE3JkkSr+olIvKgFkf6mzZtAgAYDAZMmzYNkyZNgo+PD4D6sp4yZcpNvWBgYCBKS0uhUqlw9uzZ\nFq/uzfjdwjkmkwkmXm1FRNSA2WyG2cl15Ft8I3fOnDkN1ttpPCJft26dQy9gsVgwYcKEBm/k9urV\nC88++yxWrlyJixcv8o1cIiIXkP3snRkzZiA7OxsXLlxAYGAgli9fjokTJyIxMdG2XDNP2SQicg2X\nlH5SUhIyMzNtxVxeXo60tDSsXbvWdUkbh2LpExE5zSWnbObn5zcYiffo0QOHDx9uezoiIvI4h9be\nKSsrs90vKyuD1Wp1aygiInIPu+fpp6WlYfDgwUhMTIQQAhs3bsTzzz/viWxERORiDr2Re/ToUXz7\n7beQJAnx8fFu3xidc/pERM5zyRu5Tz/9NJKTkxEVFeXScK1h6RMROc8lb+RGRETgscceg9FoxJtv\nvolLly65LCAREXmWw+fp//TTT3j33Xfx8ccfY+jQoUhNTUVcXJx7QnGkT0TkNJeM9AHAarXip59+\nwo8//oiAgABER0fj73//O6ZNm+aSoER0g5NX1RM5xe5If/Hixdi6dSvi4+ORkpICo9Fo+1j//v1x\n/Phx14fiSJ8ULCOj/kbkLEe60+4pm1qtFi+++CK6du3a5GP//e9/bz4dERF5nN3S79u3r+03xwcf\nfIDDhw9j0aJFuPfee5usmUNEN8dsvjGts2zZjcdNpvobkavYnd7RaDTIz8/HDz/8gDlz5iAlJQVZ\nWVnIzs52XyhO75CCcXqHbpZL3sj19vZGp06dsHnzZixcuBALFy7E5cuXXRaSiIg8x+70jp+fH1as\nWIEPP/wQOTk5sFqtuHbtmieyESkSp3PIneyO9LOysnD77bdj7dq1UKlUKCkpwTPPPOOJbESKxNIn\nd2pxTr+6uhpvvvkmfvnlF2i1WiQnJ8Pb2+5/DFwTinP6REROa9PaO4mJifDx8cHQoUOxY8cO3Hvv\nvcjMzHRL0CahWPpERE5rU+lrNBrbvra1tbUYNGgQcnNzXZ+yuVAsfSIip7Xp7J3fT+V4alqHiIjc\nq8WRvpeXF+644w7b/erqanTp0qX+kyQJFRUV7gvFkT4RkdPatAwDt0QkIrr1OLTKJhER3RpY+kRE\nCsLSdyGug05E7R1L34VY+kTU3rH0iYgUhCfgtxHXQSeijoSl30aNy53roBNRe8bpHSIiBWHpuxCn\nc4iovbO7XaIcuAwDEZHzXLJdIhER3TpY+kRECsLSJyJSEJY+EZGCsPSJiBREtouzQkJC0K1bN3h5\neaFz5844cOCAXFGIiBRDttKXJAlmsxk9e/aUKwIRkeLIOr3Dc/GJiDxLttKXJAkjR46EwWDAW2+9\nJVcMIiJFkW16Z+/evQgKCsL58+cxatQoDBgwAMOGDbN9PON3K5eZTCaYuMYBEVEDZrMZZic38mgX\nyzAsW7YMvr6+SEtLA+DYpcRmM9e6ISL6vXa7DENVVRUuX74MALhy5Qq++uoraDQap47BXaqIiJwn\ny/TOuXPnMHnyZABAbW0tZs2ahdGjR8sRhYhIUWQp/dDQUOTl5Tn9edylioiobTrUzlncpYqIqG24\nDAMRkYJ02NLndA4RkfPaxSmbjXHnLCIi57XbUzaJiEgeLH0iIgVh6RMRKQhLn4hIQVj6REQKwtIn\nIlIQlj4RkYKw9ImIFISlT0SkICx9IiIFYekTESkIS5+ISEFY+kRECsLSJyJSEJY+EZGCsPSJiBSE\npU9EpCAsfSIiBWHpExEpCEufiEhBWPpERArC0iciUhCWPhGRgrD0iYgUhKVPRKQgLH0iIgVh6RMR\nKQhLn4hIQVj6REQKwtInIlIQlj4RkYLIUvo7duzAgAED0LdvX6xatUqOCEREiuTx0rdarXjiiSew\nY8cOHDt2DOvXr8ePP/7o6RguYTab5Y7gEOZ0LeZ0rY6QsyNkdJTHS//AgQMIDw9HSEgIOnfujOnT\np+Pzzz/3dAyX6Cj/EJjTtZjTtTpCzo6Q0VEeL/2SkhLcfffdtvvBwcEoKSnxdAwiIkXyeOlLkuTp\nlyQiouuEh3333XciISHBdn/FihVi5cqVDZ4TFhYmAPDGG2+88ebELSwszG4HS0IIAQ+qra1F//79\n8c0336BPnz4wGo1Yv349IiIiPBmDiEiRvD3+gt7eeP3115GQkACr1Yrk5GQWPhGRh3h8pE9ERPJp\nd1fkdoQLt+bNm4fAwEBoNBq5o7SquLgYcXFxiIqKglqtxurVq+WO1KyrV68iNjYWOp0OkZGReO65\n5+SO1CKr1Qq9Xo8JEybIHaVFISEh0Gq10Ov1MBqNcsdp0cWLFzF16lREREQgMjIS+/fvlztSE8eP\nH4der7fdunfv3m5/jl566SVERUVBo9Fg5syZ+O2335p/okvfpW2j2tpaERYWJoqKikRNTY2Ijo4W\nx44dkztWE7t37xaHDx8WarVa7iitOnv2rMjNzRVCCHH58mXRr1+/dvn3KYQQV65cEUIIce3aNREb\nGytycnJkTtS8V155RcycOVNMmDBB7igtCgkJEb/++qvcMeyaPXu2eOedd4QQ9d/3ixcvypyodVar\nVahUKnHq1Cm5ozRRVFQkQkNDxdWrV4UQQiQmJop333232ee2q5F+R7lwa9iwYejRo4fcMexSqVTQ\n6XQAAF9fX0RERODMmTMyp2reHXfcAQCoqamB1WpFz549ZU7U1OnTp/Hll18iJSUFop3Pirb3fJcu\nXUJOTg7mzZsHoP69vu7du8ucqnU7d+5EWFhYg+uM2otu3bqhc+fOqKqqQm1tLaqqqnDXXXc1+9x2\nVfq8cMt9LBYLcnNzERsbK3eUZtXV1UGn0yEwMBBxcXGIjIyUO1ITixcvxssvv4xOndrVj00TkiRh\n5MiRMBgMeOutt+SO06yioiIEBARg7ty5GDhwIFJTU1FVVSV3rFZt2LABM2fOlDtGs3r27Im0tDTc\nc8896NOnD/z9/TFy5Mhmn9uu/vXywi33qKysxNSpU5GZmQlfX1+54zSrU6dOyMvLw+nTp7F79+52\nd9n7tm3b0Lt3b+j1+nY/it67dy9yc3Oxfft2/POf/0ROTo7ckZqora3F4cOHsWDBAhw+fBhdu3bF\nypUr5Y7VopqaGmzduhUPP/yw3FGaVVhYiH/84x+wWCw4c+YMKisr8dFHHzX73HZV+nfddReKi4tt\n94uLixEcHCxjoo7v2rVreOihh/DII49g0qRJcsexq3v37hg3bhwOHjwod5QG9u3bhy1btiA0NBQz\nZszAt99+i9mzZ8sdq1lBQUEAgICAAEyePBkHDhyQOVFTwcHBCA4OxqBBgwAAU6dOxeHDh2VO1bLt\n27cjJiYGAQEBckdp1sGDB3HfffehV69e8Pb2xpQpU7Bv375mn9uuSt9gMODEiROwWCyoqanBJ598\nggcffFDuWB2WEALJycmIjIzEokWL5I7TogsXLuDixYsAgOrqanz99dfQ6/Uyp2poxYoVKC4uRlFR\nETZs2ID4+Hi8//77csdqoqqqCpcvXwYAXLlyBV999VW7PMtMpVLh7rvvxs8//wygfr48KipK5lQt\nW79+PWbMmCF3jBYNGDAA+/fvR3V1NYQQ2LlzZ8tTpB56c9lhX375pejXr58ICwsTK1askDtOs6ZP\nny6CgoKEj4+PCA4OFmvXrpU7UrNycnKEJEkiOjpa6HQ6odPpxPbt2+WO1URBQYHQ6/UiOjpaaDQa\n8be//U3uSK0ym83t9uydkydPiujoaBEdHS2ioqLa7c+QEELk5eUJg8EgtFqtmDx5crs9e6eyslL0\n6tVLVFRUyB2lVatWrRKRkZFCrVaL2bNni5qammafx4uziIgUpF1N7xARkXux9ImIFISlT0SkICx9\nIiIFYekTESkIS5+ISEFY+uSU0tJSTJ8+HeHh4TAYDBg3bpztgjpPXwQ0btw4VFRUePQ1W7Jx40ZE\nRkbi/vvvb/C42Wx22zLMJpMJhw4duumPN+bOrNR+eHznLOq4hBCYPHky5s6diw0bNgAACgoKcO7c\nOVmWy/jiiy88/poteeedd/D222/jvvvu89hrSpLU6npVXMuKmsORPjls165d8PHxwWOPPWZ7TKvV\nYujQoQ2eZ7FYMHz4cMTExCAmJgbfffcdgPqR5IgRIzBp0iSEhYVhyZIl+OCDD2A0GqHVanHy5EkA\nwJw5c7BgwQIMHjwYYWFhMJvNSEpKQmRkJObOnWt7nZCQEJSVlcFisSAiIgKPPfYY1Go1EhIScPXq\nVQD1C1GNHTsWBoMBw4cPx/HjxwHUj8w1Gg10Oh1GjBgBADh69ChiY2Oh1+sRHR2NX375pcnfwfr1\n66HVaqHRaLBkyRIAwPLly7F3717MmzcP6enpDZ4vSRIqKyvx8MMPIyIiAo888ojtYy+88AKMRiM0\nGg0ef/xx2+MmkwlLlixBbGws+vfvjz179gCoX6Ji+vTpiIyMxJQpU1BdXQ2gfoXSOXPmQKPRQKvV\nIjMz03asjRs3NjlOS9+f3/v+++8xcOBAFBUV4dChQzCZTDAYDBgzZgxKS0ubPJ86EE9eJkwdW2Zm\npli8eHGzHysqKrJtKlNVVWXbzOHnn38WBoNBCCHErl27hL+/vygtLRW//fab6NOnj1i6dKnt2IsW\nLRJCCJGUlCRmzJghhBDi888/F35+fuLIkSOirq5OxMTEiPz8fCHEjc1CioqKhLe3t+3xxMRE8eGH\nHwohhIiPjxcnTpwQQgixf/9+ER8fL4QQQqPRiDNnzgghhLh06ZIQQognn3xSfPTRR0KI+k09qqur\nG3yNJSUl4p577hEXLlwQtbW1Ij4+XmzevFkIIYTJZBKHDh1q8veya9cu0b17d1FSUiLq6urE4MGD\nxZ49e4QQQpSVldme9+ijj4qtW7fajvXnP/9ZCFG/LMnIkSOFEPUbuCQnJwsh6peu8Pb2FocOHRIH\nDx4Uo0aNsh3r+tfT0nFa+/6MHz9e7N27V8TExIji4mJRU1MjBg8eLC5cuCCEEGLDhg1i3rx5Tb5O\n6jg4vUMOc3S6oKamBk888QTy8/Ph5eWFEydO2D42aNAgBAYGAgDCw8ORkJAAAFCr1di1a5ftda7P\nLavVaqhUKttiXFFRUbBYLNBqtQ1eMzQ01PZYTEwMLBYLrly5gn379jVYDrempgYAMGTIECQlJSEx\nMRFTpkwBAAwePBh//etfcfr0aUyZMgXh4eENXuP7779HXFwcevXqBQCYNWsWdu/ejYkTJwJoeeMS\no9GIPn36AAB0Oh0sFguGDBmCb7/9Fi+//DKqqqpQVlYGtVqN8ePHA4At08CBA2GxWAAAOTk5+NOf\n/gQAtlE9AISFheHkyZN46qmnMG7cOIwePdr22s0dp7Xvz48//ojHH38cX3/9NVQqFY4cOYKjR4/a\n1ma3Wq22r4U6Jk7vkMOioqIcemPw1VdfRVBQEAoKCnDw4MEGe3Xedttttj936tTJdr9Tp06ora21\nfczHx6fJc5p7XnPH9fLygtVqRV1dHXr06IHc3Fzb7ejRowCAN954Ay+++CKKi4sRExODsrIyzJgx\nA1u3bkWXLl3wwAMP2H4JXSdJUoNiF0I0+EXY0i/F5rJdvXoVCxcuxKZNm1BQUIDU1FTblNTvP8fL\ny6vB19vcLxZ/f3/k5+fDZDLhzTffREpKSqvHae37ExQUhC5dutiWORZCICoqyvb3V1BQgB07djT7\ndVLHwNInh8XHx+O3335rsBtTQUGBba74uoqKCqhUKgDA+++/D6vV6tGcQH1Z+fn5ITQ0FJ9++qnt\nsYKCAgD1c/1GoxHLli1DQEAATp8+jaKiIoSEhODJJ5/ExIkT8cMPPzQ45qBBg5CdnY1ff/0VVqsV\nGzZssL0f4KzrBd+rVy9UVlZi48aNdj9n+PDh+PjjjwEAR44csX0t1/NMmTIFL7zwAnJzc1s9Tmvf\nH39/f2zbtg3PPfccsrOz0b9/f5w/f962afm1a9dw7Ngx579gajdY+uSUzz77DDt37kR4eDjUajWe\nf/5526Yd10e6CxYswHvvvQedTofjx4832K2rpdFw4zNRHBlBN/785u5/9NFHeOedd6DT6aBWq7Fl\nyxYAQHp6uu0N2SFDhkCr1SIrKwsajQZ6vR5Hjx5tsklKUFAQVq5cibi4OOh0OhgMBrunOLZ0ho2/\nvz9SU1OhVqsxZsyYVrexvP758+fPR2VlJSIjI7F06VIYDAYA9duMxsXFQa/X49FHH8VLL73U6nHs\nfX969+7vvLA0AAAAV0lEQVSNbdu2YeHChcjPz8enn36KZ599FjqdDnq9vtk3fqnj4NLKREQKwpE+\nEZGCsPSJiBSEpU9EpCAsfSIiBWHpExEpCEufiEhBWPpERArC0iciUpD/B08EDZeqEri2AAAAAElF\nTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "It looks like there may be some sort of straight line relationship. For example, if we draw some sort of straight line through the points, it looks kind of OK. The line we chose was just a guess by eye. It has intercept $10$ and slope $0.9$:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "plt.plot(clammy, psychopathy, '+')\n", "def my_line(x):\n", " # My prediction for psychopathy given clamminess\n", " return 10 + 0.9 * x\n", "x_vals = [0, max(clammy)]\n", "y_vals = [my_line(0), my_line(max(clammy))]\n", "plt.plot(x_vals, y_vals)\n", "plt.xlabel('Clamminess of handshake')\n", "plt.ylabel('Psychopathy score')\n", "plt.title('Clammy vs psychopathy with guessed line')" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 8, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEZCAYAAAB7HPUdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVHX+P/DXcFNU5A7DRQUhk9swCOJ6BdxMTd3M+/2u\nldrmZTPdLmq2Xratfra7lW2lW6ZmW5lpum0mapalCaLgFZm8AYLc5CYwvH9/8PXUMCCgAwPO6/l4\nzOPBzJw5532Y4TWHz+ecz0clIgIiIrIIVuYugIiImg5Dn4jIgjD0iYgsCEOfiMiCMPSJiCwIQ5+I\nyIIw9BvBihUrMHnyZHOXcd/T6XSwsrJCZWWluUvBpk2b0Ldv3ybdZmhoKA4ePFjr87GxsXjvvfea\nsCLziI+PR4cOHWp93srKChcvXgQAPPnkk3j55ZebqrRmiaF/l7Zs2YKoqCg4ODjA29sbjzzyCA4f\nPgwAUKlUZq6OGlNz+bI5deoU+vXrB6DmAw2VSsXPYjVvvfUWnn/+eXOXYVYM/bvw2muvYeHChXj+\n+edx/fp1XL58GfPmzcOXX34JAOD1bpaB7zO1RAz9BsrPz8fy5cvx5ptvYvjw4bC3t4e1tTWGDBmC\ntWvX1via0aNHw8vLC05OToiJiUFKSory3LRp0zB37lw88sgjcHBwQN++fZGRkYGnn34azs7OCAoK\nQmJiorK8n58f/va3v0Gj0cDBwQEzZ85EZmYmBg8eDEdHRwwYMAB5eXkAgCFDhuAf//iHQS0ajQZf\nfPGFUY2DBw/GP//5T4PHwsPDsWPHDgDAwoUL4enpCUdHR2g0GiQnJ9e4r7GxsVi2bBl69OgBR0dH\nDB8+HLm5uQCA0tJSTJo0CW5ubnB2dkZ0dDSuX7+OTz75BFFRUQbree211zB8+HAAQElJCRYvXgw/\nPz84OTmhb9++uHXrlrLs5s2b0alTJ7i7u2P16tXK47du3cKCBQvg4+MDHx8fLFy4EGVlZQCqmgR8\nfX2xZs0auLu7w9/fH1u2bFFeu3v3bkRERMDR0REdO3bEypUrleduH107OTmhffv2OHLkiHJE/cwz\nz8DFxQWdO3fG3r17AaDO/fut/fv3Q6PRKPcHDBiA6Oho5X7fvn2xc+dOAFWfhX379mHv3r1Ys2YN\nPv74Yzg4OCAiIkJZXqfToU+fPmjfvj0GDhyIGzdu1Pi+AcBf//pXeHt7w9fXF++++65Bs0j1pqLq\nzVlnzpzBgAED4Orqiq5du+KTTz5Rnvvqq68QEhKC9u3bw9fXF6+++ioAIDs7G0OHDoWzszNcXV3R\nr18/5Yv02rVrGDlyJDw8PNC5c2f8/e9/V9ZXUlKCadOmwcXFBSEhITh69Git+1TdtGnT8MILLwD4\n9TPw2muvwdPTE97e3ti0aZOy7K1bt/CnP/0JnTp1glqtxpNPPonS0tJ6b6vZEmqQPXv2iI2Njej1\n+lqXWb58uUyaNEm5v3HjRiksLJSysjJZsGCBaLVa5bmpU6eKm5ubHD9+XEpLS6V///7SqVMn+fDD\nD6WyslKef/55iYuLU5b38/OTnj17yvXr1+Xq1avi4eEhERERkpiYqLx+5cqVIiKyfft26dGjh/La\nxMREcXV1lfLycqOaP/jgA+ndu7dyPzk5WZycnKSsrEz27t0rkZGRkp+fLyIiZ86ckfT09Br3PSYm\nRnx8fCQ5OVmKiopk5MiRyu/i7bfflmHDhklJSYlUVlbK8ePHpaCgQG7duiUuLi5y+vRpZT1arVY+\n++wzERGZO3euxMXFybVr10Sv18sPP/wgt27dkrS0NFGpVDJnzhwpLS2VEydOSKtWreTMmTMiIvLC\nCy9Iz549JSsrS7KysqRXr17ywgsviIjI/v37xcbGRhYvXixlZWVy4MABadu2rZw9e1ZEROLj4+XU\nqVMiIpKUlCSenp6yY8cOERHR6XSiUqkMPgMbN24UW1tbeffdd6WyslLeeust8fb2FhGR0tLSO+7f\nbxUXF0vr1q3lxo0bUlZWJh4eHuLr6yuFhYVSXFws9vb2kpOTo3wW9u3bJyIiK1askMmTJxu9FwEB\nAXL+/HkpKSmR2NhYWbp0aY3v2549e0StVktKSooUFxfLxIkTRaVSSWpqqoiIxMbGynvvvWewv336\n9BERkcLCQvH19ZVNmzaJXq+XhIQEcXNzU/ZXrVbLd999JyIieXl5cvz4cRERWbp0qTzxxBNSUVEh\nFRUVyjJ6vV66desmq1atkvLycrl48aJ07txZ/vvf/4qIyLPPPiv9+vWT3NxcuXz5soSEhEiHDh1q\n3C8RMdiPadOmGX0Gli9fLhUVFfLVV19JmzZtJC8vT0REFixYII8++qjk5ubKzZs3ZdiwYbJs2bJa\nt9NSMPQbaPPmzaJWq++4TPXQ/63c3FxRqVRSUFAgIlUfwjlz5ijP//3vf5fg4GDlflJSkjg5OSn3\n/fz8ZMuWLcr9kSNHyty5cw1eP3z4cBERKSkpEWdnZ7lw4YKIiCxevFjmzZtXY10FBQXStm1buXTp\nkoiI/PnPf5aZM2eKiMi+ffukS5cucuTIkTt+2YlUhcNv/zBSUlLEzs5O9Hq9vP/++9KrVy9JSkoy\net0TTzwhzz33nIiInDp1SpydnaWsrEz0er3Y29vX+JrboX/16lXlsejoaPn4449FRCQgIED27Nmj\nPPff//5X/Pz8ROTXP/ji4mLl+TFjxsiqVatq3K+nn35aFi5caLDd6qEfGBio3C8qKhKVSiWZmZl3\n3L+a9O3bVz777DP54Ycf5OGHH5axY8fK3r175dtvvxWNRqMs99vQr+kzFxsbK3/5y1+U+2+++aYM\nGjSoxm1Onz5d/vznPyv3L1y4UO/Q37Ztm/Tt29dgfXPmzFEOPjp27CgbNmxQDhpue/HFF+XRRx9V\nPp+3HTlyRDp27Gjw2OrVq2X69OkiIgZfACIi77zzjvj6+ta4XyLGof/888+LSNVnwN7e3uB99PDw\nkB9//FEqKyulbdu2yutERL7//nvx9/evdTstBZt3GsjV1RXZ2dn17sTT6/VYunQpAgMD4ejoCH9/\nfwBV/9re5uHhofzcunVrg/v29vYoLCw0WKenp6fB87+937p1a2X51q1bY8yYMfjwww8hIti2bVut\nZxU5ODhgyJAh2Lp1KwBg27ZtmDhxIgCgf//+mD9/PubNmwdPT088/vjjuHnzZq37/NszKTp27Ijy\n8nLcuHEDkydPxsCBAzFu3Dj4+Pjg2WefRUVFBQBg6tSpSvPKhx9+iLFjx8LW1hbZ2dkoLS1FQEBA\nrdtTq9XKz23atFH2/9q1a+jUqZNBLdeuXVPuOzs7w97eXrnfqVMn5fkff/wRcXFx8PDwgJOTEzZs\n2HDHppGa6gCg1FLb/tUkJiYG8fHxOHToEGJiYhATE4MDBw7g4MGDiI2NvWMNd6qpps/Sbenp6Qbv\nm6+vb7238csvv+DHH3+Es7OzctuyZQsyMzMBAJ9++im++uor+Pn5ITY2FkeOHAFQ1RQWGBiIhx9+\nGAEBAVi3bp2yvmvXrhmsb82aNbh+/TqAqve1+mfsbrm6usLK6tcYvP35ycrKQnFxMSIjI5UaBg8e\nbPB321Ix9BuoZ8+eaNWqFT7//PN6Lb9lyxbs3LkT+/btQ35+PtLS0gCYthPwTuuaOnUqPvroI3zz\nzTdo06YNevToUeuy48ePx9atW/HDDz+gtLQUcXFxynNPPfUUjh07hpSUFJw7dw6vvPJKreu5dOmS\nwc+2trZwc3ODjY0NXnzxRSQnJ+P777/Hrl278MEHHwAAfve738HOzg4HDx7E1q1blS8nNzc3tG7d\nGhcuXKj37+M2b29v6HQ6g1q8vb2V+7m5uSguLlbu//LLL/Dx8QEATJgwAcOHD8eVK1eQl5eHJ554\nQvmiv5szYmrbv5rExMRg//79Ssjf/hI4cOAAYmJianzNvZ6l4+XlhcuXLyv3f/szALRt2xZFRUXK\n/YyMDOXnjh07IiYmBrm5ucrt5s2bSh9RVFQUduzYgaysLAwfPhxjxowBALRr1w5/+9vfkJqaip07\nd+K1117Dt99+i44dO8Lf399gfQUFBdi1a5dSa/XPWEPU53fl5uYGe3t7pKSkKDXk5eWhoKCgQdtq\njhj6DeTo6IiXXnoJ8+bNwxdffIHi4mKUl5djz549ePbZZ42WLywsRKtWreDi4oKioiL8+c9/Nnje\nlOFfk549e0KlUuFPf/oTpkyZcsdlH3nkEfzyyy9Yvnw5xo0bpzx+7Ngx/PjjjygvL0ebNm3QunVr\nWFtb17gOEcHmzZtx+vRpFBcX48UXX8To0aOhUqkQHx+PkydPQq/Xw8HBAba2tgbrmTx5MubPnw87\nOzv06tULQNU51jNmzMCiRYuQnp4OvV6PH374QemQvZPx48fj5ZdfRnZ2NrKzs/HSSy8Zhe3y5ctR\nXl6OQ4cOYffu3Rg9ejSAqvfN2dkZdnZ2+Omnn7BlyxYlLNzd3WFlZYXU1NQ6a/itmvavJr169cLZ\ns2dx9OhRREdHIzg4WDmavt2JXJ1arYZOpzP6PNX38zVmzBhs3LgRZ86cQXFxMVatWmXwvFarxWef\nfYaSkhJcuHDBoFN3yJAhOHfuHDZv3ozy8nKUl5fj6NGjOHPmDMrLy/HRRx8hPz8f1tbWcHBwUN7z\nXbt24cKFCxARtG/fHtbW1rC2tkZ0dDQcHBzw17/+FSUlJdDr9Th16hSOHTum1LpmzRrk5eXhypUr\nBp28dZGqJu06l7OyssLs2bOxYMECZGVlAQCuXr2Kr7/+ut7baq4Y+ndh0aJFeO211/Dyyy/Dw8MD\nHTt2xJtvvonHHnsMgOH50VOmTEGnTp3g4+OD0NBQJYRvq34udU3nVtd1ZFLX66dMmYKTJ09i0qRJ\nd1yPnZ0dRowYgX379mHChAnK4wUFBZgzZw5cXFzg5+cHNzc3PPPMM7XWMnnyZEybNg1eXl4oKyvD\nG2+8AaDq6HD06NFwdHREcHAwYmNjDUJ48uTJSE5ONqrzb3/7G8LCwtC9e3e4urpi2bJlyh/unX43\nzz//PKKioqDRaKDRaBAVFWVwjrZarYazszO8vb0xefJkbNiwAV26dAEAvPnmm3jxxRfRvn17rFq1\nCmPHjlVe16ZNGzz33HPo3bs3XFxc8OOPP9brfatt/6pr06YNIiMjERISAhsbGwBVXwS3f/c1uf1l\n5erqanCmUF2fjdsGDRqEP/7xj4iLi0OXLl3Qs2dPAECrVq0AVJ29ZWdnB09PT0yfPh2TJk1S1uXg\n4ICvv/4a27Ztg4+PD7y8vLBs2TLli3nz5s3w9/eHo6Mj3nnnHXz00UcAgAsXLmDAgAFwcHBAr169\nMG/ePMTExMDKygq7du1CYmIiOnfuDHd3d8yZM0c5yl6+fDk6deoEf39/DBo0CFOmTLnj56Cuv7fa\nrFu3DoGBgfjd736nnBl37ty5WpdvMRqrs+DSpUsSGxsrwcHBEhISIuvXrxeRqg4nHx8f0Wq1otVq\nDTraqHF88MEHRh1tjaV6h19DFBcXi4ODg1HHXmPYv3//HTv/GkNT7t+9SklJEWtr6zo77qnlsWms\nLxNbW1u8/vrr0Gq1KCwsRGRkJAYMGACVSoVFixZh0aJFjbVp+o3i4mL885//xPz585tsm3KXTVZv\nvfUWoqOj79hp25I19/37/PPP8cgjj6C4uBjPPvss/vCHPxh0ctL9odFCX61WK2cOtGvXDkFBQbh6\n9SoAXsnYVP773/9i5MiRGDBggEFzTWO7m05FPz8/qFQq5WKwptCUQxSYY/8a6p133sH06dNhbW2N\n2NhYvPnmm+YuiRqBSpoggXU6HWJiYpCcnIxXX30VGzduhKOjI6KiovDqq6/CycmpsUsgIiI0QUdu\nYWEhRo0ahfXr16Ndu3Z48sknkZaWhsTERHh5eWHx4sWNXQIREd3WmB0GZWVl8vDDD8vrr79e4/Np\naWkSGhpq9HhAQIAA4I033njjrQG3gICAOnO50Y70RQQzZ85EcHAwFixYoDyenp6u/Pz5558jLCzM\n6LWpqanK+bTN+bZ8+XKz18A6WSfrZI23b/W5dqTROnIPHz6MzZs3Q6PRKKP+rV69Glu3bkViYiJU\nKhX8/f2xYcOGxiqBiIiqabTQ79OnT43j0wwePLixNklERHXgSbj3oKGDX5kL6zQt1mlaLaHOllBj\nfTXJKZsNpVKp0AzLIiJq1uqTnTzSJyKyIAx9IiILwtAnIrIgDH0iIgvC0CcisiAMfSIiC8LQJyKy\nIAx9IiILwtAnIrIgDH0iIgvC0CcisiAMfSJqkPh4c1dA94KhT0QNwtBv2Rj6REQWpNEmUSGi+0d8\n/K9H+CtX/vp4bGzV7X4WH39/7SNDn4jqVD3cV6wwUyFmcL+FPpt3iIgsCI/0iahB7qej3trcz81Z\nnC6RiOgOVqxoOc1ZnC6RiIgMMPSJiO6gpTfnVMfmHSKi+wSbd4iIyABDn4jIgjD0iYgsCEOfiMiC\nMPSJiCwIQ5+IyIIw9ImILAhDn4jIgjD0iYgsCEOfiMiCMPSJiCwIQ5+IyII0WuhfvnwZcXFxCAkJ\nQWhoKN544w0AQE5ODgYMGIAuXbrg4YcfRl5eXmOVQERE1TTaKJsZGRnIyMiAVqtFYWEhIiMjsWPH\nDmzcuBFubm5YsmQJ1q1bh9zcXKxdu9awKI6ySUTUYGYdZVOtVkOr1QIA2rVrh6CgIFy9ehU7d+7E\n1KlTAQBTp07Fjh07GqsEIiKqpknG09fpdIiJicGpU6fQsWNH5ObmAgBEBC4uLsp9pSge6RMRNVh9\nsrPRJ0YvLCzEyJEjsX79ejg4OBg8p1KpoFKpanzdit9MShkbG4vY+236GiKiexQfH4/42zO411Oj\nHumXl5dj6NChGDx4MBYsWAAA6Nq1K+Lj46FWq5Geno64uDicOXPGsCge6RMRNZhZ2/RFBDNnzkRw\ncLAS+ADwhz/8Af/+978BAP/+978xfPjwxiqBiIiqabQj/e+++w79+vWDRqNRmnDWrFmD6OhojBkz\nBpcuXYKfnx+2b98OJycnw6J4pE9E1GD1yU5OjE5EdJ/gxOhERGSAoU9EZEEY+kREFoShT0RkQRj6\nREQWhKFPRGRBGPpERBaEoU9EZEEY+kREFoShT0RkQRj6REQWhKFPRGRBGPpERBaEoU9EZEEY+kRE\nFoShT0RkQRj6LUQD5z4mIqpRnaFfVFSEVatWYfbs2QCA8+fPY9euXY1eGBli6BORKdQZ+tOnT4ed\nnR2+//57AIC3tzeee+65Ri+MiIhMz6auBVJTU7F9+3Zs27YNANC2bdtGL4qqxMf/eoS/cuWvj8fG\nVt2IiBqqztBv1aoVSkpKlPupqalo1apVoxZFVaqH+4oVZiqEiJqNisoKnM0+i4SMBBxPP46EjASk\n5aYh7em0er2+ztBfsWIFBg0ahCtXrmDChAk4fPgwNm3adK91ExFRHUrKS3Dy+kkkpCcgIaPqdur6\nKfg4+CDCKwLd1N2wrM8yRKgjoFKp6rXOO4Z+ZWUlcnNz8emnn+LIkSMAgPXr18Pd3f3e94YahM05\nRPe3vNI8JGYkGgR8ak4qurh2QTevbohQR2CSZhLCPcPh0MrhrrejEhG50wKRkZH4+eef73oDd0Ol\nUqGOsoiIWqz0m+lVwf5/AX88/TiuF11HuDocEeqIqptXBELcQ9DKpv7N6fXJzjpDf+nSpXBzc8PY\nsWMNOnFdXFzqXUhDMfSJ6H4gIriYe1EJ+OMZx5GQnoCKygpEeEUYBPwDLg/A2sr6nrZnktD38/Mz\naitSqVS4ePHiPRV3x6IY+kTUwpTry3Em+4zSuZqQkYATGSfQvlV7JeBvN9P4tvetdxt8Q5gk9M2B\noU9EzVlxeTGSMpMM2t9TslLQoX0HpYM1wisCWrUWbm3cmqwuk4R+WVkZ3nrrLRw8eBAqlQoxMTF4\n4oknYGtra9JiDYpi6BNRM5FbkmvQ/n77FMkg9yCD5hmNpwbt7NqZtVaThP7MmTNRUVGBqVOnQkTw\n4YcfwsbGBu+++65JizUoiqFPRE1MRHDt5jWjDtackhyjDtZg92DYWduZu2QjJgl9jUaDpKSkOh8z\nJYY+ETWmSqlEak6qwQVOCekJAGDUwRroEggrVcsYm7I+2VnnxVk2Nja4cOECAgMDAVRdkWtjU+fL\niIiahXJ9OVKyUow6WF3sXZSAn999Prp5dYO3g3ejdLA2J3Wm9yuvvIL+/fvD398fAKDT6bBx48ZG\nL4yIqKGKyopwIvOEQfv76azT8HPyUzpYh3cdDq1aCxf7xjvtvDmr19k7paWlOHv2LFQqFbp06YLW\nrVs3blFs3iGiOtwovmHUwfpL3i8I8Qgx6mBtY9vG3OU2CZO06f/jH//AxIkT4ezsDADIzc3F1q1b\nMXfuXNNVWr0ohj4R/R8RwZWCK0YdrPm38qFVaw0CPsgtCLbWjXdmYXNnktAPDw/HiRMnDB7TarVI\nTEy89wprK4qhT2SRKqUS52+cN+pgtbGyMehg7ebVDf7O/i2mg7WpmKQjt7KyEpWVlbCyqvrl6vV6\nlJeX16uAGTNmYPfu3fDw8MDJkycBVI3a+e677yqDtq1ZswaDBg2q1/qI6P5xq+IWkrOSDZpnkjKT\n4N7GXQn4hb9biAh1BLwcvMxdbrMWH1//QRnrDP2BAwdi3LhxePzxxyEi2LBhQ71Devr06Xjqqacw\nZcoU5TGVSoVFixZh0aJF9auQiFq8m7duGnWwns0+i87OnZWhCUYFj4JWrYVTaydzl9vimDT0161b\nh3feeQdvvfUWAGDAgAGYNWtWvVbet29f6HQ6o8fZdEN0/8oqyjJqf7968ypCPUIRoY5AD58eeCLq\nCYR5hMHe1t7c5VqcOkPf2toaTz75JJ588knk5OTg8uXLsLa+t5Hg/v73v+ODDz5AVFQUXn31VTg5\nNe43e0O+BYmofkQEl/IvGY0gWVhWqDTPDHlgCJ7v9zy6unWFjRWv7zGl2qZTrUudHbkxMTH48ssv\nUVFRgcjISLi7u6N37954/fXX67UBnU6HYcOGKW36169fV9rzX3jhBaSnp+O9994zLMrEHbkrVnCq\nQaJ7oa/U49yNcwYXOCVmJKKVdSujK1j9nfzv+wucmpvbGWeSjtz8/Hy0b98e7777LqZMmYKVK1ci\nLCzsrovz8PBQfp41axaGDRtW43IrfpPSsbGxiOWhOlGTKK0oxanrpwza309mnoS6nVq5wOmZXs8g\nQh0Bz3ae5i7XosXHx//frf4HtnWGvl6vR3p6OrZv346XX34ZAO7pWzw9PR1eXlU98Z9//nmtXyAr\n7vHQvLZ/fapPNk5kyQpuFRhN0Xf+xnkEugQqHazjQsch3DMcjq0dzV0uVXP7gPh2E/bKerTz1Bn6\nL774IgYOHIjevXsjOjoaqampeOCBB+pV0Pjx43HgwAFkZ2ejQ4cOWLlyJeLj45GYmAiVSgV/f39s\n2LChXutqqOrhzuYdsnSZhZlGHawZhRkI8wxDhDoCvTv0xvzo+Qj1CEVrm8a96p5MqyEHshYxiQrb\n9MmSiAh0eTqjC5xKK0qN2t8fdH3wnqfoo+aDM2f9H569Q/erisoKnM0+a9TB2ta2rdEVrB0dO7KD\n9T7H0Ce6j5SUl+Dk9ZMG7e+nrp+Cj4OPwRysWrUWHm096l4h3XdMEvp6vf6ez8tvKIY+Wbq80jyj\nDtbUnFR0ce2idLBGeEUg3DMcDq0czF0uNRMmCf3OnTtj5MiRmD59OoKDg01aYK1FMfTJgqTfTDfq\nYL1edN1oir4Q9xC0smll7nKpGTNJ6BcUFGDbtm3YtGkT9Ho9ZsyYgfHjx6N9+/YmLdagKIY+3YdE\nBBdzLxp1sFZUVhh1sD7g8gA7WKnBTN6mHx8fj4kTJyI3NxejR4/GCy+8oEyjaEoMfWrpyvXlOJN9\nxqiD1bGVo1EHq297X3awkkmY5IrciooK7N69Gxs3boROp8PixYsxYcIEfPfdd3jkkUdw7tw5kxVM\n1BIVlxcjKTPJoP09JSsFHdp3UK5gHdplKLRqLdzauJm7XLJwdYZ+ly5dEBsbiyVLlqBXr17K46NG\njcKBAwcatTii5ia3JNeo/V2Xp0NXt65KB+s07TRoPDVoZ9fO3OUSGamzeefmzZtwcGjaswPYvEPm\nJiK4dvOa0QiSOSU5Rh2swe7BsLO2M3fJRKZp079+/Tr+9a9/QafToaKiQlnx+++/b7pKqxfF0Kcm\nVCmVSM1JNWh/T0hPAACjDtZAl0BO0UfNlklCv2fPnujXrx8iIyOVKRNVKhVGjhxpukqrF8XQp0ZS\npi9DSlaKQfv7iYwTcLF3Mepg9XbwZgcrtSgmCf3GngS9Jgx9MoWisiKjKfpOZ52Gn5Of0sEa4RUB\nrVoLF3sXc5drUhx6xDKZ5OydoUOHYvfu3RgyZIjJCiMytRvFNww6WBMyEvBL3i8I8QhBhDoCkV6R\nmNVtFjSeGrSxbWPuchsdQ59qU+uRfrt27ZR/bYuKimBnZwdbW9uqF6lUKCgoaLyieKRPtRARXCm4\nYtTBmn8rH1q11qD9PcgtCLbWtuYu2Sw4sqxl4oBr1KJVSiXO3zhvdAWrjZWNUQdrZ+fOFt/BWn3i\noOXLq37mxEGWwySh//vf/x779u2r8zFTYuhbnlsVt5CclWzQPJOUmQT3Nu4GI0hGqCPg5eBl7nKb\nPR7pW6Z7atMvKSlBcXExsrKykJOTozxeUFCAq1evmq5Ksjg3b9006mA9m30WnZ07Kx2so4JHQavW\nwqm1k7nLJbqv1Br6GzZswPr163Ht2jVERkYqjzs4OGD+/PlNUhy1fFlFWUZXsF69eRWhHqGIUEeg\nh08PPBH1BMI8wmBva2/ucu8bbM6h2tTZvPPGG2/gj3/8Y1PVA6BhzTs8S6F5EBFcyr9k1P5eWFZo\n1P7e1a0rbKzqPHGMiBrIZB25p06dQkpKCkpLS5XHpkyZcu8V1lZUA0KfbZdNT1+px7kb54xGkGxl\n3coo4P2d/HmBE1ETMcl5+itWrMCBAweQnJyMIUOGYM+ePejTp0+jhj41H6UVpTh1/ZRB+/vJzJNQ\nt1MrAf+jCh3XAAAXrklEQVRMr2cQoY6AZztPc5dLRHWoM/T/85//4MSJE+jWrRs2btyIzMxMTJw4\nsSlqq1X1U9Nu46lp96bgVoHRFH3nbpzDAy4PKGfOjAsdh3DPcDi2djR3uWRmbFptmeoMfXt7e1hb\nW8PGxgb5+fnw8PDA5cuXm6K2WlUPdzbvNFxmYaZRB2t6YTo0nhpEqCPQu0NvzI+ej1CPULS2aW3u\ncqkZYui3THWGfvfu3ZGbm4vZs2cjKioKbdu2NRhXn5o3EYEuT2fUwVpaUao0zzz64KNYEbsCD7o+\nyCn6iO5zDboiV6fToaCgABqNpjFr4tk7d6misgJnss8YNM8kZiSirW1bow7WTo6d2MFKDcarfps3\nk5y9IyL47LPP8N1330GlUqFv37547LHHTFqoUVEmuCL3fv8yKCkvwcnrJw0C/tT1U/Bx8DG4glWr\n1sKjrYe5y6X7EM+ca35McvbO3LlzkZqaivHjx0NEsGHDBvzvf//Dm2++abJCG8P9FPp5pXlGHayp\nOano4tpFuYJ1kmYSwj3D4dCqaWc5I6KWpc7Q379/P1JSUpQJVKZNm4bg4OBGL8xSpd9MNxpB8nrR\ndWWKvphOMVjwuwUIcQ9BK5tW5i6XLNj9clBlaeoM/cDAQFy6dAl+fn4AgEuXLiEwMLCx67orLelU\nThHBxdyLRh2sFZUVSvPMqKBR+Ev/v+ABlwfYwUrNTnP7m6L6qbNNv1+/fjh69Ciio6OhUqnw008/\noXv37mjfvj1UKhV27txp+qJM0KbfnNoby/XlOJ192qiD1bGVo9EUfb7tfdnBSkR3xSRt+i+99JKy\nMgAGK2Q4GSsuL0ZSZpJBwCdfT0ZHx45KwA/tMhRatRZubdzMXS4RWZh6nbKZkZGBo0ePQqVSITo6\nGh4ejXs2SEs5eye3JNfoAiddng5d3boqV7BGeEVA46lBO7t2jVsMEVk8k5yyuX37djzzzDOIiYkB\nABw8eBCvvPIKRo8ebbpKqxfVzCZRERFcu3nNqIP1RskNoyn6gt2DYWdtZ+6SicgCmST0NRoNvvnm\nG+XoPisrC7///e+RlJRkukqrF2XG0K+USqTmpBqMIJmQngAARhc4BboEWvwUfUTUfJikTV9E4O7u\nrtx3dXVtVkfh96JMX4aUrBSD9vcTGSfgYu+iBPz87vPRzasbvB282YdBRC1enaE/aNAgDBw4EBMm\nTICI4OOPP8bgwYObojaTKiorMpqi73TWafg5+SkBP7zrcGjVWrjYu5i7XCKiRtGowzDMmDEDu3fv\nhoeHB06ePAkAyMnJwdixY/HLL7/Az88P27dvh5OT4Tyo99q8c6P4hlEH66X8Swh2DzbqYG1j2+au\nt0NE1JyYbOasu3Xo0CG0a9cOU6ZMUUJ/yZIlcHNzw5IlS7Bu3Trk5uZi7dq1DS4cqPpCulJwxegC\np/xb+UYdrEFuQbC1tm2U/SQiag5MEvqffvopli5diszMTGVlKpUKBQUF9SpCp9Nh2LBhSuh37doV\nBw4cgKenJzIyMhAbG4szZ87UWbi+Uo/zOecNmmcS0hNgY2Vj1MHa2bkzO1iJyOKYpCN3yZIl2LVr\nF4KCgkxSVGZmJjw9q6bV8/T0RGZmZo3LHU8/bhDwSZlJcG/jrgT8gh4L0M2rG7wcvExSFxGRJagz\n9NVqtckCvzqVSlXrGTFTPp+ijCA5KngUtGotnFo71bgsERHVT62h/+mnnwIAoqKiMHbsWAwfPhx2\ndlUXHalUKowYMeKuNni7WUetViM9Pb3Wq3tHXR8FXAfyT+QDsYCTHwOfiOi34uPjEX97lMl6qrVN\nf9q0aQbj7VQ/It+4cWO9NlC9TX/JkiVwdXXFs88+i7Vr1yIvL++uO3KJiOhXZj97Z/z48Thw4ACy\ns7Ph6emJl156CY8++ijGjBmjDNfcGKdsEhFZIpOE/tSpU7F+/XolmHNzc7F48WK8//77pqu0elEM\nfSKiBqtPdtZ5XuOJEycMjsSdnZ1x/Pjxe6+OiIiaXJ2hLyLIyclR7ufk5ECv1zdqUURE1DjqPGVz\n8eLF6NmzJ8aMGQMRwSeffILnnnuuKWojIiITq1dHbnJyMr799luoVCr079+/0SdGZ5s+EVHDmaQj\nd9GiRZg5cyZCQkJMWtydMPSJiBrOJB25QUFBmDNnDqKjo/H2228jPz/fZAUSEVHTqvd5+mfOnMGm\nTZuwZcsW9OnTB7Nnz0ZcXFzjFMUjfSKiBjPJkT4A6PV6nDlzBqdPn4a7uzvCw8Px2muvYezYsSYp\nlIh+1cCr6okapM4j/YULF+LLL79E//79MWvWLERHRyvPPfjggzh79qzpi+KRPlmwFSuqbkQNZZKh\nlTUaDV5++WW0bdvW6Lkff/zx7qsjIqImV2foP/DAA8o3x4cffojjx49jwYIF6NSpk9GYOUR0d+Lj\nf23WWbny18djY6tuRKZSZ/NOWFgYTpw4gZMnT2LatGmYNWsWtm/fjgMHDjReUWzeIQvG5h26Wybp\nyLWxsYGVlRV27NiBefPmYd68ebh586bJiiQioqZTZ/OOg4MDVq9ejc2bN+PQoUPQ6/UoLy9vitqI\nLBKbc6gx1Xmkv337drRu3Rrvv/8+1Go1rl69imeeeaYpaiOySAx9aky1tumXlJTg7bffxoULF6DR\naDBz5kzY2NT5j4FpimKbPhFRg93T2DtjxoyBnZ0d+vTpg71796JTp05Yv359oxRqVBRDn4iowe4p\n9MPCwpR5bSsqKtC9e3ckJCSYvsqaimLoExE12D2dvfPbppymatYhIqLGVeuRvrW1Ndq0aaPcLykp\ngb29fdWLVCoUFBQ0XlE80iciarB7GoaBUyISEd1/6jXKJhER3R8Y+kREFoShb0IcB52ImjuGvgkx\n9ImouWPoExFZEJ6Af484DjoRtSQM/XtUPdw5DjoRNWds3iEisiAMfRNicw4RNXd1TpdoDhyGgYio\n4UwyXSIREd0/GPpERBaEoU9EZEEY+kREFoShT0RkQcx2cZafnx/at28Pa2tr2Nra4qeffjJXKURE\nFsNsoa9SqRAfHw8XFxdzlUBEZHHM2rzDc/GJiJqW2UJfpVLhoYceQlRUFP71r3+ZqwwiIotituad\nw4cPw8vLC1lZWRgwYAC6du2Kvn37Ks+v+M3IZbGxsYjlGAdERAbi4+MR38CJPJrFMAwrV65Eu3bt\nsHjxYgD1u5Q4Pp5j3RAR/VazHYahuLgYN2/eBAAUFRXh66+/RlhYWIPWwVmqiIgazizNO5mZmXjs\nsccAABUVFZg4cSIefvhhc5RCRGRRzBL6/v7+SExMbPDrOEsVEdG9aVEzZ3GWKiKie8NhGIiILEiL\nDX025xARNVyzOGWzOs6cRUTUcM32lE0iIjIPhj4RkQVh6BMRWRCGPhGRBWHoExFZEIY+EZEFYegT\nEVkQhj4RkQVh6BMRWRCGPhGRBWHoExFZEIY+EZEFYegTEVkQhj4RkQVh6BMRWRCGPhGRBWHoExFZ\nEIY+EZEFYegTEVkQhj4RkQVh6BMRWRCGPhGRBWHoExFZEIY+EZEFYegTEVkQhj4RkQVh6BMRWRCG\nPhGRBWHoExFZEIY+EZEFYegTEVkQs4T+3r170bVrVzzwwANYt26dOUogIrJITR76er0e8+fPx969\ne5GSkoKtW7fi9OnTTV2GScTHx5u7hHphnabFOk2rJdTZEmqsryYP/Z9++gmBgYHw8/ODra0txo0b\nhy+++KKpyzCJlvJBYJ2mxTpNqyXU2RJqrK8mD/2rV6+iQ4cOyn1fX19cvXq1qcsgIrJITR76KpWq\nqTdJRES3SRP74YcfZODAgcr91atXy9q1aw2WCQgIEAC88cYbb7w14BYQEFBnBqtERNCEKioq8OCD\nD2Lfvn3w9vZGdHQ0tm7diqCgoKYsg4jIItk0+QZtbPCPf/wDAwcOhF6vx8yZMxn4RERNpMmP9ImI\nyHya3RW5LeHCrRkzZsDT0xNhYWHmLuWOLl++jLi4OISEhCA0NBRvvPGGuUuqUWlpKXr06AGtVovg\n4GAsW7bM3CXVSq/XIyIiAsOGDTN3KbXy8/ODRqNBREQEoqOjzV1OrfLy8jBq1CgEBQUhODgYR44c\nMXdJRs6ePYuIiAjl5ujo2Gz/jtasWYOQkBCEhYVhwoQJuHXrVs0LmrSX9h5VVFRIQECApKWlSVlZ\nmYSHh0tKSoq5yzJy8OBBOX78uISGhpq7lDtKT0+XhIQEERG5efOmdOnSpVn+PkVEioqKRESkvLxc\nevToIYcOHTJzRTV79dVXZcKECTJs2DBzl1IrPz8/uXHjhrnLqNOUKVPkvffeE5Gq9z0vL8/MFd2Z\nXq8XtVotly5dMncpRtLS0sTf319KS0tFRGTMmDGyadOmGpdtVkf6LeXCrb59+8LZ2dncZdRJrVZD\nq9UCANq1a4egoCBcu3bNzFXVrE2bNgCAsrIy6PV6uLi4mLkiY1euXMFXX32FWbNmQZp5q2hzry8/\nPx+HDh3CjBkzAFT19Tk6Opq5qjv75ptvEBAQYHCdUXPRvn172Nraori4GBUVFSguLoaPj0+Nyzar\n0OeFW41Hp9MhISEBPXr0MHcpNaqsrIRWq4Wnpyfi4uIQHBxs7pKMLFy4EK+88gqsrJrVn40RlUqF\nhx56CFFRUfjXv/5l7nJqlJaWBnd3d0yfPh3dunXD7NmzUVxcbO6y7mjbtm2YMGGCucuokYuLCxYv\nXoyOHTvC29sbTk5OeOihh2pctll9ennhVuMoLCzEqFGjsH79erRr187c5dTIysoKiYmJuHLlCg4e\nPNjsLnvftWsXPDw8EBER0eyPog8fPoyEhATs2bMH//znP3Ho0CFzl2SkoqICx48fx9y5c3H8+HG0\nbdsWa9euNXdZtSorK8OXX36J0aNHm7uUGqWmpuL//b//B51Oh2vXrqGwsBAfffRRjcs2q9D38fHB\n5cuXlfuXL1+Gr6+vGStq+crLyzFy5EhMmjQJw4cPN3c5dXJ0dMSQIUNw7Ngxc5di4Pvvv8fOnTvh\n7++P8ePH49tvv8WUKVPMXVaNvLy8AADu7u547LHH8NNPP5m5ImO+vr7w9fVF9+7dAQCjRo3C8ePH\nzVxV7fbs2YPIyEi4u7ubu5QaHTt2DL169YKrqytsbGwwYsQIfP/99zUu26xCPyoqCufPn4dOp0NZ\nWRk+/vhj/OEPfzB3WS2WiGDmzJkIDg7GggULzF1OrbKzs5GXlwcAKCkpwf/+9z9ERESYuSpDq1ev\nxuXLl5GWloZt27ahf//++OCDD8xdlpHi4mLcvHkTAFBUVISvv/66WZ5lplar0aFDB5w7dw5AVXt5\nSEiImauq3datWzF+/Hhzl1Grrl274siRIygpKYGI4Jtvvqm9ibSJOpfr7auvvpIuXbpIQECArF69\n2tzl1GjcuHHi5eUldnZ24uvrK++//765S6rRoUOHRKVSSXh4uGi1WtFqtbJnzx5zl2UkKSlJIiIi\nJDw8XMLCwuSvf/2ruUu6o/j4+GZ79s7FixclPDxcwsPDJSQkpNn+DYmIJCYmSlRUlGg0Gnnsscea\n7dk7hYWF4urqKgUFBeYu5Y7WrVsnwcHBEhoaKlOmTJGysrIal+PFWUREFqRZNe8QEVHjYugTEVkQ\nhj4RkQVh6BMRWRCGPhGRBWHoExFZEIY+NUhGRgbGjRuHwMBAREVFYciQIcoFdU19EdCQIUNQUFDQ\npNuszSeffILg4GD8/ve/N3g8Pj6+0YZhjo2Nxc8//3zXz1fXmLVS89HkM2dRyyUieOyxxzB9+nRs\n27YNAJCUlITMzEyzDJexe/fuJt9mbd577z28++676NWrV5NtU6VS3XG8Ko5lRTXhkT7V2/79+2Fn\nZ4c5c+Yoj2k0GvTp08dgOZ1Oh379+iEyMhKRkZH44YcfAFQdScbExGD48OEICAjA0qVL8eGHHyI6\nOhoajQYXL14EAEybNg1z585Fz549ERAQgPj4eEydOhXBwcGYPn26sh0/Pz/k5ORAp9MhKCgIc+bM\nQWhoKAYOHIjS0lIAVQNRDR48GFFRUejXrx/Onj0LoOrIPCwsDFqtFjExMQCA5ORk9OjRAxEREQgP\nD8eFCxeMfgdbt26FRqNBWFgYli5dCgB46aWXcPjwYcyYMQNLliwxWF6lUqGwsBCjR49GUFAQJk2a\npDy3atUqREdHIywsDI8//rjyeGxsLJYuXYoePXrgwQcfxHfffQegaoiKcePGITg4GCNGjEBJSQmA\nqhFKp02bhrCwMGg0Gqxfv15Z1yeffGK0ntren986evQounXrhrS0NPz888+IjY1FVFQUBg0ahIyM\nDKPlqQVpysuEqWVbv369LFy4sMbn0tLSlElliouLlckczp07J1FRUSIisn//fnFycpKMjAy5deuW\neHt7y/Lly5V1L1iwQEREpk6dKuPHjxcRkS+++EIcHBzk1KlTUllZKZGRkXLixAkR+XWykLS0NLGx\nsVEeHzNmjGzevFlERPr37y/nz58XEZEjR45I//79RUQkLCxMrl27JiIi+fn5IiLy1FNPyUcffSQi\nVZN6lJSUGOzj1atXpWPHjpKdnS0VFRXSv39/2bFjh4iIxMbGys8//2z0e9m/f784OjrK1atXpbKy\nUnr27CnfffediIjk5OQoy02ePFm+/PJLZV1/+tOfRKRqWJKHHnpIRKomcJk5c6aIVA1dYWNjIz//\n/LMcO3ZMBgwYoKzr9v7Utp47vT9Dhw6Vw4cPS2RkpFy+fFnKysqkZ8+ekp2dLSIi27ZtkxkzZhjt\nJ7UcbN6heqtvc0FZWRnmz5+PEydOwNraGufPn1ee6969Ozw9PQEAgYGBGDhwIAAgNDQU+/fvV7Zz\nu205NDQUarVaGYwrJCQEOp0OGo3GYJv+/v7KY5GRkdDpdCgqKsL3339vMBxuWVkZAKB3796YOnUq\nxowZgxEjRgAAevbsib/85S+4cuUKRowYgcDAQINtHD16FHFxcXB1dQUATJw4EQcPHsSjjz4KoPaJ\nS6Kjo+Ht7Q0A0Gq10Ol06N27N7799lu88sorKC4uRk5ODkJDQzF06FAAUGrq1q0bdDodAODQoUN4\n+umnAUA5qgeAgIAAXLx4EX/84x8xZMgQPPzww8q2a1rPnd6f06dP4/HHH8f//vc/qNVqnDp1CsnJ\nycrY7Hq9XtkXapnYvEP1FhISUq+Owddffx1eXl5ISkrCsWPHDObqbNWqlfKzlZWVct/KygoVFRXK\nc3Z2dkbL1LRcTeu1traGXq9HZWUlnJ2dkZCQoNySk5MBAG+99RZefvllXL58GZGRkcjJycH48ePx\n5Zdfwt7eHo888ojyJXSbSqUyCHYRMfgirO1LsabaSktLMW/ePHz66adISkrC7NmzlSap377G2tra\nYH9r+mJxcnLCiRMnEBsbi7fffhuzZs2643ru9P54eXnB3t5eGeZYRBASEqL8/pKSkrB3794a95Na\nBoY+1Vv//v1x69Ytg9mYkpKSlLbi2woKCqBWqwEAH3zwAfR6fZPWCVSFlYODA/z9/fGf//xHeSwp\nKQlAVVt/dHQ0Vq5cCXd3d1y5cgVpaWnw8/PDU089hUcffRQnT540WGf37t1x4MAB3LhxA3q9Htu2\nbVP6AxrqdsC7urqisLAQn3zySZ2v6devH7Zs2QIAOHXqlLIvt+sZMWIEVq1ahYSEhDuu507vj5OT\nE3bt2oVly5bhwIEDePDBB5GVlaVMWl5eXo6UlJSG7zA1Gwx9apDPP/8c33zzDQIDAxEaGornnntO\nmbTj9pHu3Llz8e9//xtarRZnz541mK2rtqPh6mei1OcIuvrra7r/0Ucf4b333oNWq0VoaCh27twJ\nAFiyZInSIdu7d29oNBps374dYWFhiIiIQHJystEkKV5eXli7di3i4uKg1WoRFRVV5ymOtZ1h4+Tk\nhNmzZyM0NBSDBg264zSWt1//5JNPorCwEMHBwVi+fDmioqIAVE0zGhcXh4iICEyePBlr1qy543rq\nen88PDywa9cuzJs3DydOnMB//vMfPPvss9BqtYiIiKix45daDg6tTERkQXikT0RkQRj6REQWhKFP\nRGRBGPpERBaEoU9EZEEY+kREFoShT0RkQRj6REQW5P8DWD9d2y5RIIoAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "What does our straight line relationship mean?\n", "\n", "We are saying that the values of `psychopathy` can be partly predicted by a straight line of formula `10 + clammy * 0.9`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To make this more general, let's call our `psychopathy` data $y$ - a vector with 12 values, one for each student. $y_1$ is the value for the first student (= 11.416) and $y_i$ is the value for student $i$ where $i \\in 1 .. 12$.\n", "\n", "Our `clammy` score is a predictor. Lets call the clammy scores $x$ - another vector with 12 values. $x_1$ is the value for the first student (= 0.389) and $x_i$ is the value for student $i$ where $i \\in 1 .. 12$.\n", "\n", "Our straight line model says:\n", "\n", "$y_i \\approx c + bx_i$\n", "\n", "where $c$ is the intercept and $b$ is the slope. For the guessed line above:\n", "\n", "$y_i \\approx 10 + 0.9 x_i$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the $\\approx$ above, we are accepting that we will not succeed in explaining all the variation in our psychopathy data. We can rephrase this by saying that each observation is equal to the predicted value (from the formula above) plus some error for each observation:\n", "\n", "$y_i = c + bx_i + e_i$" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Simple regression in matrix form" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We're going to rephrase the simple regression model in matrix form. Let's make the data and predictor and errors into vectors.\n", "\n", "$y$ is the vector of values $y_1 .. y_{12}$:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "y = Matrix(psychopathy)\n", "y" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left[\\begin{matrix}11.416\\\\4.514\\\\12.204\\\\14.835\\\\8.416\\\\6.563\\\\17.343\\\\13.02\\\\15.19\\\\11.902\\\\22.721\\\\22.324\\end{matrix}\\right]$$" ], "metadata": {}, "output_type": "pyout", "png": "iVBORw0KGgoAAAANSUhEUgAAAEoAAAErCAMAAABO0eUIAAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\nTlMAMquZdlQQQO0wRM0iu+/dZol8bIlZUvoAAAtcSURBVHgB7Zxrk7M4DoUJEN6Zzg12+f+/dXWx\nfSQZ0yTVU9U1Sz6kDRIntgA/yDbdXVb59N3Hn5cqdN1lHUb6XD9W6mY+vl9Z6vK5Co6co9RtykaU\n0p7LkArFsizjUqoRpO6PNUmhlKW710OKxTLdvrpuuGW7k5oeQ69SKGXHrhufLGUst4W2H6pPJSdF\n22OulSnRbvrc51s6KvnM610t6fu41NIFqcfTKR2v1XyPUuvtsowDhevNWk1jF6Sm9TaTyLNcjkcb\nyBH2DZxWCevwzFfPQakrXz1eqltf3LJ5zU08JjVxpaLUU66oeeVm8ueY1LXnz/rsKWLlgtEz+G6t\nWIAi7K+ra4qVGg/Uik5U9l29VPegZk+tMzj0z/XVc+NReqV+bHmsK5tg6ZZh6MvdXMUq1+CDvzHs\nH0jkQ06pHInv/56x+j5G2eOfj1WhJ/3kbZy6S+7DC00J6pmwqVqNWiV6itOTny2ko+oKTdlgfXh7\nW0rpyXb63JZ+FOIZmtJu58N+m1KgJ7vYhyUQ1/s0pUDPtpT3aUkZeorUPI9D6uFKrYJPQ8rSk10e\nBLp7QlSWij4NKUtPdpHPQ6BXnkm2fOqwO3pmpW7Qx5ZUq00flvqz/imHdJ6etL8X4izaQpWqfOTw\n/8QHSE9P8lEID/rcpVKVj0jVDZTdmZ60obfHS2mYw85OxkeOaUgJPZWmd9aan/qEZ6UyYUWHfWID\n2ZDpqTS9D0s/yJOPoWnxyUINKZjfKW3W6h0B+J5SiMV3pR+O1UQc+IHPnfr/fyof/Lx6Pxwr18DE\nzPswPB54ZO2wPQ2UJBWL4WqoVWbmxMRK3R032WzfKD+aXjkVNFx1UmBmL+f0WTI+bM+yb0kP85ar\nTop+PnVIT6FCX1J9bD+EsJpJmMyVDm1I3aT/hRS2V5H6WiUNdFxtSHF86AFDO3Td0O1plQ76S55H\nPFf3pL5Wyd2KEm/fdd+F/wSu7kmZsyNyvG2lAld3pIaUaOVayfakteIGRq6ylEUqkDLahyGSS9sa\nq+t6rbgakVqkrhzee74Suy5va6vndaq42mqg3hD6tMctLNujXKJlJMZwtSF1fy3LMvR0oSaulu3u\nxTfOM9+EhqteqjDzpsNkfM0LV832tNCoV1LK7OV6V1e77Pzsy9fqM4101Cl1PHw/HKsTqcdCX4W9\nkYZO47hopnqnG/0hkHDpa3UPZqRKPWyHLHDoSUGeva8yKmfSV/Z3tQJS2WRxeZUBvTuBdBRsyZBh\nSV/ZO0jRdulFPS4X6ecn6nS+njTiqJ1Y6LFdrZyUw+XXylnhkgcbtYFHpTwuu35d56LU3TR0Jn3d\naWDAJWVN65pO7mV5iZJLX3ekAi67ebjcVg6TfJbSt3cpfW1LVbjkuMxrGS9+llJKX5NUjdQKl3rk\nhS4EReODSi59JaltpEZcTsRw/tC4qI5nP+gadekrGdsXg0tDU2b5uKfc9EVZucY+pa/fSJlU9UvO\n3kwP+TNfW1d+vHLpayVVkEqWjEtNVS/9kG7nr4EenuVUIn0l97qBsvOzrxirz1TkqFPqePB+OFb2\nHjxeicoz3oOVw/EdP9xAl1oer0X03KgV6AnnkowW4rppVHHckCr0hFIe5C3EDdOoDSnQE1JKV0Pc\nMI3akAI9ixQGeRNx4zRqQyrQk71A1yQVp1EbUoGe5GXomqTiNGpLytBTXCxdVYoeGbhbxiRcS8rT\nk5sXJyyradQi5R9r50BPR9fcQBk7wISlSNUjRaCnOHi6Jqk4jVpq5W4cQ09x8HRNUnEadVuqAz3F\nQb9yMpqk4jRqQwr01Gc7ccvJaJKK06gNqQ70PDyN2pKS/R98bfQMH6ictXozaD8c9hOpx+K/GXaP\nS5uMZgtSVfzMhlTEJZLRYjGp6q5UxCWS0WIxqeqeVIXLkkHCYlLVPakKl0UqWDRV3ZOqcNnnZDRY\nNFXdkapxmedSnaWkqrtSaSWNn0WhZDSC1KSqoscXg0dqXHWkP8vJaLQgVRWfDaTKAD9WHSEZBUiR\nqqJ9VWpJ08wyHAuIIxmFBanqrlTEJZJRWJSwnKqaz8aNU1YdKVJNMlosSFX3pbDqKM6lwoJUFVob\ntYLxvdIpdTxevzdWJ1KPncXqDGoaennRyBB/sopnqCatfl/sr3IaOus8js4rklxgqAwJh31eCmno\noD1R4ZZnqCatfl89fpWe8bVhWK7sGJqSVrePah5jldMFjtGXLrrionwyQ5G08lBdHgzelQpTjWm4\n1yatlChqf80/tVMrDBmzY2GoSVrLPnbYkcrZnLjJlzLUJq2UwmIImGtlkYpYjWnCG1IdM9QlrWQD\nVyNSIWWnEzDca5LWyNVmrKZyeXK1CkNN0lr2pXo3pWhVf3LhP4Ghshn2tcN+TVLK1cBQSVrDviBl\nBn7veUBCuWoZmoeE7T6qeWygadS7xVPqeMR+b6xOpB47i9tnMA3vBmZ2WJmE1Ur4mS2pzNXIzK6s\nTDKrlfakwNWKmWVlElYrQWm7Z0i9cmQmViZhtdJBKXEzzMTKJKxWekcKzHQrk1girF7aCjtmjB0z\n7XIiVoqrl/al6AAwM0o5upEnS1mk8s+hVrxRmGlWJvH+uHqpi0hlpywVmVlWJolPefaSQ+hrr4GR\nmdoiWplEx+XVSlnnG6nITLMyqaxWOigFZipXsTIJq5f2pcBVMFO5WlYmmdVK0NqMFczvlE6p49H6\nvbE6kXrsLG6fwYRUO+Sb5JJF4OpGr7b7q4zUDkO+SShbpielbvc0D5XrXNcKSMXrM+oNSy9w7X2W\nV0vRYTnRqfrcbNEXARe/9PYjqUmnwMbymqXUel8qD/nmcORa6ajkgmSXHXal8pBvUcpSiyzBuaU3\nmpJ5V0p8sP6INnMUb5ShXzdqtY9UrD8yUvze9WXRJdKpUrtIxZBv9i61kh3pheZs3Gsghnyzt5e6\n+YtlTwpDvlFKXoS5+HHf3TOIId+E1BL2kW7FSWbEy49sXgxAKtYfKVJhoTdSdYnSvpSxvlXcjNVb\nCsX5lCqh+Lbww7H6/0JqAh0FGSXaAPia+WUIewYdUa78pwY+dQZ8zfzSSQF0KLEQraMs4MObLzG/\ndFJ0TO68bYl2A3zt/PKYlAFfO788JpXWhhJidvLLg1IFfDEpNPnlQSl6vVjBF6VMfnlUKoNvJ79k\nKXsPNs4gXxB0SdCjSzu/jEsEd6UYfO388mADAb52fvmtlCLQgK+8+RLzSy8F0KGkCAT4mvmll9Lo\nfvh9Sh0P3O+NVfVYe7xR1nNzpMg6HC//3lhh4SIGYVGSFpoJUWex2HUNxCAsShorTIhai8eukwIk\nUVIpTIjCErHrpDAIi5JKYULUWdDnkpuTwiAsSiol3zIh6ixtKT0Mg7AoicWkysnyjRQgiZIomQnR\nbPlGCpBEiaXshGi27EthEBYlqZSZEC2WSurPX3+LM3/hFVKU1Jhr4nzs8MB//3LvhmMQFiVVwoSo\nsVS1wo2DQViUVIn/zZmWrKUtBUiilLPKPCEKCym3pTAIi5K+TMqDOfLyWmcswK7U113tqS0f/jml\njgfujNW/K1bAJUraQsyFeov2E+LjLgbgEqWkxD2MTBUGCzpEDy/gEiWVwlyot9gXTl2tgEuUVApz\noc6CV2LIzUkBlyiplHxLA53Fri7yUnoYQIqSWDAXmuY+zSsx5OBqJQdkXIYZTjcXqj5mdREfWUvh\nnKAkv2HmQtXiVxfVUgWX9QxnmQtVn7i6KNYKIEWJ6uTmQtViVhdppYMUcIkSO9q50GQxq4u2pIBL\nlMTPzIV6ixrEx9UKuERJkYq5UFjk+PxKDG04KeASpYTUMhdqLDS/nf6zYF0r2fPxl6vVxypnrd4M\nnYZdlnT6mYa3dMq/wKU3IvmT3yF/S0Sd5V/gjmP3P/kG0DHexS6MAAAAAElFTkSuQmCC\n", "prompt_number": 9, "text": [ "\u23a111.416\u23a4\n", "\u23a2 \u23a5\n", "\u23a24.514 \u23a5\n", "\u23a2 \u23a5\n", "\u23a212.204\u23a5\n", "\u23a2 \u23a5\n", "\u23a214.835\u23a5\n", "\u23a2 \u23a5\n", "\u23a28.416 \u23a5\n", "\u23a2 \u23a5\n", "\u23a26.563 \u23a5\n", "\u23a2 \u23a5\n", "\u23a217.343\u23a5\n", "\u23a2 \u23a5\n", "\u23a213.02 \u23a5\n", "\u23a2 \u23a5\n", "\u23a215.19 \u23a5\n", "\u23a2 \u23a5\n", "\u23a211.902\u23a5\n", "\u23a2 \u23a5\n", "\u23a222.721\u23a5\n", "\u23a2 \u23a5\n", "\u23a322.324\u23a6" ] } ], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "$x$ is the vector of values $x_1 .. x_{12}$:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = Matrix(clammy)\n", "x" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left[\\begin{matrix}0.389\\\\0.2\\\\0.241\\\\0.463\\\\4.585\\\\1.097\\\\1.642\\\\4.972\\\\7.957\\\\5.585\\\\5.527\\\\6.964\\end{matrix}\\right]$$" ], "metadata": {}, "output_type": "pyout", "png": "iVBORw0KGgoAAAANSUhEUgAAAD8AAAErCAMAAABesavUAAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\nTlMAMquZdlQQQO0wRInN3SJm77t8bMVussMAAAoNSURBVHgB7VzpkrM2EJQBk8RgMAnv/6yZQ9Ic\nSMIkW5WvUvDDC2haF2KamYYNj522Llzb3gwL4bH3A2zPa/CwIqjbEf+4CBXz1eDHfhp6VdkyTdNG\n/RqHYepfiCObJdVg8TPYju9cuPRg9txX+MW90GHxBw6XT7Ix+PWDVtOMv7gNO7a9v6ESrCMsUNKR\nTbfhMWwGv9EleO4jl4XXB3o87oCayH6EmvAIjpONwe+Ef1GrsYrY/9e+QaXTCrVRTcNOc2HbH3ca\n5WsfMhh2Zh76vq8AD+ET26cBWfyyTwh88B/cDY/pTXCYwH2fcVwTXDCoNLWh+1/AI+CDF3TtH/OO\n0xFm6N1zoouC59X1H7lh1//w+YxhxZlZd9iDP8PwmNIcaXzg8T9TWVjoKm9wTEAYGa0lqKMrzv+b\n5nZNZWGnvW1fRxo0wOKshTndbKb9gdZGTzMMjYQP1ffelxDX27aEFXcfcIo3gw9vWqAwX7RoworN\nPXGqXzT3K+wNsDfOqR9m/gA2wR1J98+bl1Lf9xutlEfXx/tnSqeoA7b92KkLf278T/q/CxMfTe/5\nL8//gUlg1W+KWdJM1+bPMcmIt9OUnG4Cw98K3jNJR5TwIf+g0FW8Z5IPMAd6rcMIKu17JpnJk32N\nLzNJeGfPmMdQbr/MBC9FLKmCK3h2zwnJf8v4IpP05I0tnK7fb/tv7uyBSaB8SC7f2P6p+EsVeCYB\nN448yoSk7GrrxzNJeBAND4k2chXl8QfPJMsbnqT6jlZRxuJODe+YZOZnxe/xppHGQa39BsQU3fiy\n/zOT1DjA+RsH8k4Nq3rRAo8kP/f8X2+nVvKz17/NGnwPQk/oYYx7ZNtvs0b0X8um3KjBt1ljoMfJ\ncevzwyt0weCbrLGsc3SAQ639JmtM4QzfZI11OcW3WGOEdXbWfhnPs46RzRm+wRpPJN4aPvOHjz8A\nw6wxUmBVwhv+qLLGs8Nt/3QwCxiWyv1qrv8Ja3A40sCfsAaHni18izWmbd8xfuq7z/7OgZTpPw7u\n4nbj/3v/e/NHeiz3/IGLObHGNA0T2klOC0vt+nX8gQaBncI4QxiPmQWV08JSg/f8gQbMGmFGD7QB\nAUhOC0st3vMHlEfWWHPCQnJaR7znD7CIrLHZyIdzcliB7n+BPxJr7PMDMoMx5wWeOCalLP7o/xNr\nQDoFMyYfSv9ITsu1f8Qn1hg5E9RzGinltI79Jzcv+SthDcwhYmfTCDCnRRuOv8YfijU4a7dC1kxy\nWoRv8YdiDZ5/bD/ntKT9/Px24A+0IdbgnGYPV1FyWlSBvn6HqIMsmDU2mJoR519yWgW84w+wSKwR\nIG3Fma1XzmlhBbZ9qvLSz43/7/kjr/9LV46Nzfr/B/hf7/onvsDBzBCZPcjtif4h59Di2P8YZWBh\n+GDYTs+M5DBR/1Dn0OKAj3yBZdD+1HHKQfQPOUcWHi9RBhbnlInoH3KuiJcow+Cz/qHqLOETX1AZ\n5uvXqCd1Sf9Q59DI9j/xRYQHFF0WdtpJ/9DnDvjEFwlPfzf0/aJ/yDncw/bz85PwBRnFnx64V+sf\neBrP0abjb8UXXNhFreel9A85RzZ6/IovGM+s0++j0j/yOTbR+AiSNBuz9Bt6IfqHnDu0z3jiC9Y/\n6FmHBBPRP+RcEZ/4gvWPpZ+6nphW9A85hxUc+8+9+Pb3xt/8kZ+fvl00yu7H14/KTYFnyXo67aHP\neLxBCcSNO+HaN7kp8P+oR6KeLqr5yplgftR394/LTYHWiY2gni6qec+OKzGDa9/kNoLEI6Kac78p\nksLKm/gcj3jV/EWceIZX8YhTzcXFtdpX8YRVzUmIx9bb/Vd4o5onb5zwmT/ghM4NmXyWUs0HogVE\nA7vBTOj1q/EHFSQmHvUTRmv8MXQVPZ1V8zGtHWy/iZd4RKvm8LoKInmr4Jk/JJ+lVfNnHS+5KeYP\niUeUar5wGFluP/Xr27++/9/ikt2Nt+s3zcu3f3H+7vgj5a9g1gx/gPOOWSucT45MFqOku/Xn+CNn\nrehykN9wSrrBe/6QrBXiOTJxSrrBo5HKjYNziU/JCI+RiVPSm3iTtYqRiVPSm3idtbKRSVbSW3id\ntbKRiSjpbTxNBmWtbGQiHhzxNf6gN9/QRb+CjUyUkt7kj5y1spGJVtJb/Q85a2UiE6OkN/GStcIF\nEFNXvIqTkl7BR/7IWSvCU2TilHSL9/whWaucyXJKusVTM5d+bvy/9//6+efS5IMxzv/NH4o/st6B\nM1mKOpr8ARjxTcWoo8Uf2KTOX5WijhP+MPmrUtRxwh86f1WMOtr8YVmiGnXIm1j2/rcs4d7VVVFH\njT8sS2D+Ud7V1VGHXCNsP/OHZQm8GuWoo8IfliUQTtsh6qjxh2GJhIZBUKwpUUeTPxJLELwYdbT5\nA17wxuCQ+aMUdTT5I7NE4PijEHXc/CFXFvbs+jdFXx38CvibPyJ/uNyU4g9RQpyNuX4uN6X4Q1QP\nZ2PwLjelslaiejgbg3e5Kf2tR8p3gSCOWzV/Ba/659yU5o+MJ7jYmPapTHJT+lsPUULQSGwOeMUS\nmj+UEgK9z58fHO4/zRKWP7BZfn9U22D7mT/AROem4FDxBxxF1UPb+PhDmAnseUP+MKqHtnHjF5ZI\naOYPrXoYG4c3uSn1rYdWPYyNw6fclOcPrXokG+qhw+fclOcPpXpkG6zA4WXUX+7d+Dv++Lfxk+Ct\nto3+45r+LXo3L185phvwVP/Oendc/fn4S/3bemmlVX+pf1fxX+rf1stf17+Nl4c5UMdf6d80bdHL\nxymMXr+lf8v1J0zWtmMNTf0b/U9+fjJeHtByLN//yTmq3rw/or08lubjxFjw2mo+R3DrP7WXx2I5\n/k7/Fi/P/l+Ov9S/xcvf+jdfn7Pf/wN/6eens/H6cv/85MvPjn+t+RO+4H7Lcc46teIHiRLisDN/\nSNapFT9IlBDxmT8k69SKH9y7UYo/JOvUjB+q/OGyTrX4oc4fNKCsWtTiB8UXBLDHOevUiB8AVuQP\nrC49tydvjOcK66/EH2ias046fiD++P0PLNd8UTiWrFPqBxr99bt6/nPcYLhCsk71+EH4gttXxyrr\nVI8fhC88f+isUyN+qPKHzjrd8QNfHfotrD9Ver5749X9cz5dB4ufnj+td+v8EyQ9Dko4dsa2b/Vu\nnX+yJeJBLD5/pcHjFP2ipISTjcEbvRuKRb8wJVrjMHijd2MSm/w8/q8RU6I1DoPXejf2jv0h/q8P\nXWI0Do0HpyVfaSA+vzWlS6zGYfGidyNc8k/6+w2rcWi80rsZLvkn+X7DaRwGz+PFrzT0hvmnXOI1\nDoPPereGk36RS7zGYfBO71b5J1eSvgT06zfr3cwfSr/IJdQ1Xhi4a9qHK5a+0vD5JykBo/glYAFP\ntV/5ce1fgZLtjUf/S5t/eDyby/z/o+DbBNzSt01nuFQOuV3cwt/aAr3b4fkTPwAAAABJRU5ErkJg\ngg==\n", "prompt_number": 10, "text": [ "\u23a10.389\u23a4\n", "\u23a2 \u23a5\n", "\u23a2 0.2 \u23a5\n", "\u23a2 \u23a5\n", "\u23a20.241\u23a5\n", "\u23a2 \u23a5\n", "\u23a20.463\u23a5\n", "\u23a2 \u23a5\n", "\u23a24.585\u23a5\n", "\u23a2 \u23a5\n", "\u23a21.097\u23a5\n", "\u23a2 \u23a5\n", "\u23a21.642\u23a5\n", "\u23a2 \u23a5\n", "\u23a24.972\u23a5\n", "\u23a2 \u23a5\n", "\u23a27.957\u23a5\n", "\u23a2 \u23a5\n", "\u23a25.585\u23a5\n", "\u23a2 \u23a5\n", "\u23a25.527\u23a5\n", "\u23a2 \u23a5\n", "\u23a36.964\u23a6" ] } ], "prompt_number": 10 }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\epsilon$ is the vector of errors $e_1 .. e_{12}$:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "epsilons = [Symbol('e_{0}'.format(i)) for i in range(1, 13)]\n", "epsilon = Matrix(epsilons)\n", "epsilon" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left[\\begin{matrix}e_{1}\\\\e_{2}\\\\e_{3}\\\\e_{4}\\\\e_{5}\\\\e_{6}\\\\e_{7}\\\\e_{8}\\\\e_{9}\\\\e_{10}\\\\e_{11}\\\\e_{12}\\end{matrix}\\right]$$" ], "metadata": {}, "output_type": "pyout", "png": "iVBORw0KGgoAAAANSUhEUgAAACsAAAErCAMAAABw9IpGAAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\nTlMAMquZdlQQQO0wRIki3e9mzbt8bDjT79EAAAWaSURBVGgF7ZvtlpwgDIZRkW3Fj5nW+7/WJigk\nMUPcMz32tFv5seL4TowBHyDDumZNpXX1Mm0S55q181D6utQtKGhX1DaGjC4t/5c2RO+7kR6/1HQc\nhgl0Y0RF8yw6rCjtOEObhA4i3nTtZGvbGdpoCUnUn2hX1i1OtcnT7d5n2nnTelSfaf0DRGGL2ZnW\n+c4vy+d82FRvaGf+Td0WdHXsprVNz7h/ptqNtKp2a7eQ3HGoxuEtlsTu2W0vf8UusSQCf1p8+3JR\nbUEsQVY1KwMbaj/Wj/xV54gl8wBcW+FPLj+OY8DKWAKEWJnDygfOEucm833jLHGxy/fHo7LLWbII\n6QstWRrA9YENfdpu0TaPvu+f1rMVqZvTGEznpg9MlqqGD0fppXaDZw+r7kwfjNCI9/zhVZ+kGB1r\nv9nG77HEx5a3qPKBWNLBqNzxWYzSEks6eIXwcil4UmEJajrOKJslYWavpu7rjCXBPzhKtFaw5OTZ\nOEuwa7MZtIpDeerk67ju0xn8uK51yN9lZo1haJsI7cZQbdkt3uSKYTdLyvFK7c0Ssz+UNsiVK9vi\nwOq3WAJ+cjzoPkksAamYEmgtsQSWDpFPH5K2yhIfhNZiSdNLrYovY0l0J1piyTKeaQtLAvDpxG5u\nfNd3UFaEcC7K33whHUcRB1O7PNf2UzwTN0gnpt2D/Eqt6A+H+/JT1R/4xUP9Sn8P79vhznT6BXx4\nhyVN24w5C5CCoeJALBlgXs2SHBtTxdhCLBn6hg2wYFnNU8nUwJY3r32gNc4AWTV7TUYsCTBu9xyU\n6tmIJXjfcFyTyb5T8iVPsBv4mkzZTU+R/qS8FAeEoY0QMbHQMrQuRpz6UbG0pNpqf4v2ZskXm5fE\n6NmwqecaxBL3GNzIs4PYJyss8QAdHMBLMViCyV1RVF8nlqz9p1kyrpDVGHgGWNktLAlpyWLzobAk\nyT65bnlgbG275dEH4F7P13rK3yKFPthFm6lMe6xadv+k9mbJl2JJ7AMU1oFUPyOWtCnfxwZ61FZY\n0g5QeB7RYAmipOdzE+UDsQSkwR6PaV4CWpEK0pwsLAHp2e9ZZY0D2ihc0HZBUopM2djamVP9pJ9N\nPGIn2uLMVlHxPVznp1dqb5b8te32zhrHLfFkzUAsWSDnGRhK0ntRYQmuAxwfYw2WdFPjRr54Uv2X\nWBKm9WRuz1iyPFZcY5Si7BJL4uKayc575nlJSI/1YJMpZbfccUzh4jsd6lr3QF85/Qxt6CD5aj5b\n8UFXDLtKfKX2ZsnXYknfyUFW9R1iCa5bPP81FLUVluC0XiQS6ywZ024hcx1QWBI2rfUeE0tw7dSY\naxFiCeZVvOkDm5c03qvfu2VupbySOOXDtUspKr7linvCmvBp/kZF2gHyJRwl9nyHvpdqhg8H5bV9\n/WbJtfFlfbI6z8AGl/2B7T1TveGgJTaUeb2HMS5/TdilXGrep7aAsM9zCNSW/kC5VFCkLAVMNyBX\nsC+gRC6+sAEESbulPXPiRPhAbNi1TTI57csBoSU2SO2euBFaxobNhw0P80ttjg75m55tX2ZIu0et\nS6uRHIe6NmWikNUivqw/ZMNlnxqkX0vOpGo3f4sdr9SKsYXd81hVY8tRwM6v9PdFfNmdqfqP+/Ae\nS/K6Px8hHCIOxJLMh3xMgRNaYsn+boKE7Q9F7ce371vANUuE9uc3nCTso65iibZbtIolhlaxxNAq\nllhauJZLfv58hM9FzLIOj/2eqc7HujbzIR+TlarddFX+ubVbPO44vIzDb7GE72cV8dUsEftZhVaz\nROxnFVrNErGfVWrLNIjzpmxeEdpXLGH7WYX2FUvYz6xC+4olbD+r1LL3dmcI389a1yaWiP2sFW1m\niNjPWtEyZ6h6rTZl3fN0lG5KtfK/MAH/0cV7mRgmHdbS/8J4734BM9tr9lzGzCYAAAAASUVORK5C\nYII=\n", "prompt_number": 11, "text": [ "\u23a1e\u2081 \u23a4\n", "\u23a2 \u23a5\n", "\u23a2e\u2082 \u23a5\n", "\u23a2 \u23a5\n", "\u23a2e\u2083 \u23a5\n", "\u23a2 \u23a5\n", "\u23a2e\u2084 \u23a5\n", "\u23a2 \u23a5\n", "\u23a2e\u2085 \u23a5\n", "\u23a2 \u23a5\n", "\u23a2e\u2086 \u23a5\n", "\u23a2 \u23a5\n", "\u23a2e\u2087 \u23a5\n", "\u23a2 \u23a5\n", "\u23a2e\u2088 \u23a5\n", "\u23a2 \u23a5\n", "\u23a2e\u2089 \u23a5\n", "\u23a2 \u23a5\n", "\u23a2e\u2081\u2080\u23a5\n", "\u23a2 \u23a5\n", "\u23a2e\u2081\u2081\u23a5\n", "\u23a2 \u23a5\n", "\u23a3e\u2081\u2082\u23a6" ] } ], "prompt_number": 11 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can rephrase our model as:\n", "\n", "$$y = c + b x + \\epsilon$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bear with with us for a little trick. If $o$ is a vector of ones, then we can rewrite the formula as:\n", "\n", "$y = co + bx + \\epsilon$\n", "\n", "because $o_i = 1$ and so $co_i = c$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "c, b = Symbol('c'), Symbol('b')\n", "o = ones((12, 1))\n", "# Compile expression without any simplification or reordering\n", "rhs = MatAdd(MatMul(c, o, evaluate=False),\n", " MatMul(b, x, evaluate=False),\n", " epsilon)\n", "Eq(y, rhs)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left[\\begin{matrix}11.416\\\\4.514\\\\12.204\\\\14.835\\\\8.416\\\\6.563\\\\17.343\\\\13.02\\\\15.19\\\\11.902\\\\22.721\\\\22.324\\end{matrix}\\right] = c \\left[\\begin{matrix}1\\\\1\\\\1\\\\1\\\\1\\\\1\\\\1\\\\1\\\\1\\\\1\\\\1\\\\1\\end{matrix}\\right] + b \\left[\\begin{matrix}0.389\\\\0.2\\\\0.241\\\\0.463\\\\4.585\\\\1.097\\\\1.642\\\\4.972\\\\7.957\\\\5.585\\\\5.527\\\\6.964\\end{matrix}\\right] + \\left[\\begin{matrix}e_{1}\\\\e_{2}\\\\e_{3}\\\\e_{4}\\\\e_{5}\\\\e_{6}\\\\e_{7}\\\\e_{8}\\\\e_{9}\\\\e_{10}\\\\e_{11}\\\\e_{12}\\end{matrix}\\right]$$" ], "metadata": {}, "output_type": "pyout", "png": "iVBORw0KGgoAAAANSUhEUgAAAVkAAAErCAMAAABDxgAjAAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\nTlMAMquZdlQQQO0wRM0iu+/dZol8bIlZUvoAACAASURBVHgB7V2Jkuq6kjTrfcPS4Bn+/1unFtUm\nS5xC0OfyIkxEt2VbzkylhTF2Ik+bB7220/DrxghD2y/Zx9GWWClJLxEmOBRv89jt4XVIqWhWuuP2\n20dz3Z8WLtnH0ZZYf2Kn9S8RJjgUb/PYpAT8odJ91Nkm+xDaOy3JEmY5EC/UPR/FQCuVJZtdKdia\ned7PYkxWmuALqt+vhjyEZi057ub9TpQB02We5yu9K4/7/bw7ITnVuRQZU5bQOGTL9rRy9nJ9FGet\npBverlS0NcczKNydS4WsNMXjglNqyFO6oQHNsM7g4vGmtl2wUxwed/hP3WOLq39g9vIjdbLyjSNQ\nL2aCs8frbsvOWsm22P+gs37NeYYFV/Z70Av3jvHIg2ja6vsP6p5lp0/7B/bXxw3sRXenC6zZUp3t\nq/KVA3GevIKzUG8vfdaVePPL/VxESJ37Q/Y31cju9EqNVyrIUGUITbGudKJz0MacfuDddXyAnzM1\n4gge4xzMS50soXJU7ahn887OU+3slXa6Imal6QZc8Eo/5uyDnD1RT1U+OhqcHlc43s138Jk83j/o\nmJvflV6vQjcKaWfvl4Wzj/MGPiaKsLy0SoVX+ilnjw86mp4ee0925kPs43EHY6fpp/RZOjzk5Xu9\ngHKc9/tdeO8WyqyzR9BY9Vl4c5FA+sQFuC/qs5cHfgLASQ9NqKmb+UbGwkfY40GnIDOdJZ7F/az8\n6OwJPyUvzCZHbOIjO3xd6zVWwpq47cJZOkTtfsrZRFYaE+v/NvvYfhKshrPYhh88DbvvNucHHnan\nM/SWw0wnDLg8eTouHLAJnlrgSQae3m12Wzh2uxfi+brmp5Wg+gE1Vc7SRy0qevFA5dix2GZPNzSg\nCdaRO2t1NJh+oBPc8Qh8f1B3gO9Km1mOxUPObn8A485d6zDi7JG6e+0sH6jur+704MXvODvxcfYg\nrk0XOhJeYZ7fYRtdU04zB/ssf1Jyi4acPWzx9fjZ0ieC9GY+N/jGPjvxl5q7nFFNDypdH/djec9v\ny+fWdJYLUUN91h3IpyFneafQNwUoirN8uriTc6+stL/RZ/ckSr8eTiz9Bifg5TvX9TLdsbjRc/Ks\nfDniUDN++COSOtyrztLpNZvBJ4Dm7HQF3CMew+mVlVaqy8QrlX0G64bQDOuGnyz4icX679hFD3ga\ncKLzgjuU9lA60ukNSckSGgdstscztyOfdT11drf9edzojWKlW3m7zNfHA1fZGvi03e22+NlGr6w0\nqV+mptQjv+ssnmayNNZ/2u12V/qo3Wx35YoMqOdFL8k3vbjZfrenvQbFp84Sxei/t50NxENosdUB\n748zWcIux+psx+P3nZUPHCZAvO5e6GjoLM5KqzbvsA+hdbAqxvZslrDDcdndHnzqVOBXZ8XnN50V\nGJ2uzooVq7PiRGvaeae2qi6Wrc4uLHELVmedGb3Pz2wX8lA9rFCnN5MlzO49xMvW7Wkqy7PSKpgO\n+xBaB6tibM9mCbMcq7Pi8+qsONGaZvtTa9u/5KyGNkDDeX+EK+bl8rZFLWBZFJiVFreKxyKDH0Jr\nObvIdMDVgqte7DAxWcIWh6FYCfEadUtog+r9YByOrpNNPmpRroFGKJtLlxy7h882NPA4LF1eZTqO\neIFplvsgWit/ca3F4WC02HaWQxul0nne7umSfIxahDpYdcgLt18j/BBao9V1pmNLt1Z+4nf8V+Q3\nOIpRcdJ01kIbWFkuumPZLqDGOrhuyAvnbIAfQ2u0us50/NDdqq0Pk6H4vPwGBwNU/5vOWmgDa7ed\njXWw3nc6W2c6znTf5l9y1oU2yNn7XWN92merOl/rbDvTMd30Dhkqp1e2Y1R9Np3kABIf2kBOTOtc\nyiFfnK3rYL2sNKzrXl6pwMPqITSPxRTt5MHJ3xksWrKEkcOSHPPuuivJC0JEvFg3hDYKK9jLN9Ol\n6T7YIXWy0qR+mXp2gYdVQ2gei+HbzvpTH5GTJQwcluSY4UN+S3mmArh0NoQ2hBYiO3R6UJrerJOV\npphc8Eo/72wz07GLOSHWkZXv9UJSVJIcmEKyG8EAiXj/PP6x1sbQBiznHVHOALnpizq0+f8C1MDL\nsztnh9A8VpFSZzpg8d5/JqviLGHgsCQH5pnkmEmYiBf2QgxtQB1Owuz4kM9NX9QhqOxO17ZwwbM7\nZ4fQPFahqTMdcI8VvztyaMZLyRIGjni8trwuACNeqMtkEtqAOf4Oe+NDiGt6SUeYuKw024JKnt3B\nD6F5rEJTZzrKl3L+7uOlZAkDh09yTNONv6kyasdZCm1wCIIi/pQngS1c0ycJdoi8rDSpX6ZeqYMf\nQvNYQlNlOi43+C1InRvEulnCwOGTHJC+EE7BC3VxoYQ2OARx2c1bPp/wUQupY2BZabYFlYzdw6cb\nGtAMyxZXmY4z/ygw5gaxdlZ+5HBJjnswtnc0MGH5UlZahRiV6sohtA6Wgj4rZAl7HCc4Ez25a2iI\n16v7TEdjXVZatWmHfQitg1UxtmezhB2OzflwOOCXKnmtznonpPxs2nGWrrX6887VWXHxzT4rMDpd\nnRUrVmfFida0805tVV0sW51dWOIWfJuzR7iF+IHXxX8DyeN12IfQOlgpMVnCLAfivbOnnejs28lt\ngsUO+xBaB6tibM9mCbMciJet21akS7PSdAMudNiH0DpYFWN7NkuY5VidFZ9/3dmSpbhUeQebr5IR\nlufISpOmlGnoAx9LclQqJ9OPtCra+AavG1SNcbNVn5UsRZ13cPNVMsJufLzvrLCjviE020uVSqcf\nwYtoz5cmNA5E6r+Cs5alqPMONl8lI1yeY8gLd5Q3dpQ7hKatrlROph+xl6OLvECoHIj05BWchXrl\nCmmdd7D5mIzweY4hL5yzxo56h9C01VEl3Bjx+Q0T7a4HpwmVA1U+eXWcrfMONh+TET7PMeTFLzkb\nVULqz+c3TPS/4CzvizrvAPMxGRHyHF/kbFSpHYvb40R/wNl0ksNz1XkHnA/372Oe44ucDSrVWG6P\nF+1bm5UfjwYuybGf+Whe+BDP1/Vc9rFfujDcpg+aY54jK01byoUO+xCaYAWVysft8aJ9a7OEwkGw\nluTY4U+SfY7hibOhHuDQvE9GVHmOrDRtKRe8UtfSITTB8iqVjtsTRDu+sU8wS3LsYI8FzTjjswnG\nVecdyrwlI+o8RzYKoU3lQpt9GkJTLFOpbKw/irbWTmlC5UBkS3LgnI2nADN1kkO56ryDzPNbCke7\nqPMcYY8hUe4l/QxrK3u1/3NI7rhmKmXToj+Kdnyxw8lWjanXO4UkB42noFugHb6ucPH3P8s76Hyd\njLDMxzc5W6us8hsiWlqLdmTle7cml+Q47nEsJXt1nLW8Q4lzWP5BkxEMYnmOrDQjp5JX6lo6hGZY\nqnKhH0lFtOMbczYkOZ59gmmWwuUdKM7h5i0ZARJ9nmPIC/+OUXZs/BCaOWsqa/0mOvClCY0DZbok\nBzZFRn4pDYh1sf7Qa8gL72xgHUJ7pyVZwjYHHWMvMsAWtgTx2nVDOzMzWWkVVod9CK2DVTG2Z7OE\nHQ7Mz95luD1kWJ0Vn990FgbfnLeU3i6Aq7MfclZgdLo6K1a82WcFRqers2LF6qw40Zp2Pl1aVRfL\nfsPZbDZhISYuyEYh4lbw87NmjmQIrYNVMbZns4RZjjXJIT7/Rp995z0kumCaleY2wWKHfQitg1Ux\ntmezhFkOxMvWbSvSpVlpugEXOuxDaB2sirE9myXMciycddkGzTuAFHsmiT2nJI7W8ZE+a+zZhgab\nrNV1kgOrSXvkkTCuJbA2S2gcgXoxUzkbsg3+bg39bgSfSeKeUzK50TpekFZpcEo9e7ahAc2wqiQH\n1eL26CNhfEtekG8cgXoxE5yNWQoX0nDPJLHnlMC9Zhmtg2CHvHDHosg+hKatrpMcKLC0Rx8J41sC\nq7OEykGN7v8LziK9/vDf8g6w2J5JYs8piWNK5KVVYrxSY083NKApVp3kgFqlPfZIGN8SWP/3nLW8\nA9DqM0m4IfzUojBaR14aQ+h/dQOWfMzZOskB2KU91SNhSktg/V9z1uUd0IStPJMEZyZ+Tsm0daN1\n5KURgv37DWcbSQ5pT/VImNISkDPo7MtJDp93IBfkmSRwYNTnlPjROvLSzFQq/Yazy7yBtAdu3uBg\n6jzsu7UEFo05a0kOgPBxA8Jrtc3nHdAAeyYJzvFzSrCko3XkpfFm+r/FDiuzDVUcLAjW0llpz5Gf\nrKCPhLGWZAmFg4gtyQGz4Uf4PWdD3gE28s8kQUh8Tkl5ldE6Br0wNxDuU8fZRZLD2oNPX0OtZcw8\na8mQs5bkgFN+fdQYMpCzPpvAbYt5B6hmzySx55TE0TryUQjitX9Ldlr30SSHa489Esa1BBmzhF5v\nSHLsj3peJXi+f7OzMe8AuwO+A+MLnqegzykBu2k0iTJax1f12frpKq499kgY1xJs21CfdUmOzSE6\ni3hLZ5EJnLMjsj2ThBfic0riaB15aQyu/zvs2YYqDhYUa5HkwLUk3R4J41qCa7OEykGQePyGoxj8\nzdD/9BAJ88+cdQNz2DNJ7Dkl/EVXRuvIS0Mp7uWVfuo4C79EiE9XIT7Ob+gjYVxLcPWQs5bkgLO6\nZ876bIOENHhgDnsmiT2nBKyV0TpekUattH/mrGdPN9SAoGRYVZID1kl77JEwriWwfshZTXLgWd0z\nZ4PMV2ey0ipccyOsGELrYAXg3kyWsM1xgOfi7B4YopUX4rXrSo30NCutAuywD6F1sCrG9myWsM9x\neXqcbbOmlmalVWAdpUNoHayKsT2bJexy3K8PvMwqL8Tr1pVKuWlWWoXWYR9C62BVjO3ZLGGWY3VW\nfF6dFSda02x/am27OttyRZZ9m7Phm7CIfH2a/eJdIXfYh9A6WBVjezZLmOVAvHf2tFOZfTu5TbDY\nYR9C62BVjO3ZLGGWA/GydduKdGlWmm7AhQ77EFoHq2Jsz2YJsxyrs+LzX3DWQhtCCt1axqnUqIXE\nIaxOVpptQaXQBxQ+/TU+oAUsEw117BEx1jxbhiBZ+ZEj0IeZRp/V0IarWEIdGrXQOISrk5XmNsGi\nU6rwsHgIzWEBhE+iWOjEmmfLUEeWMHLglu3X0tkDXVW4+GHWl4NYaBzCoWaluU2wqEo/meRAYJ9E\nsdCJa95YEEX1IseT19JZC23oZhbqKBdQLQ6hdfI73W2CRa/0c9dnNblR2HTQdNc8XUZ1sh3D6y3g\nzcnS2Sq0gVtZqKM0vYpDMHJWWqXDK/2ksyYaCdVF1zxdRoqy8r3eqilhdulsFdqA2hKCgGJpehWH\nYMSstMD/a33WiUZCC51YJsWWYY2s/MrZdJIDOSy0gXPhaSvsbIhDUB38l5WmG3DBK/1cn5XkhpC5\n0Ik2zy17Qb7XC7Es/IXdBe6Gbbabizy7nUgbfTaGNvBYACct5WajOEvXeDUOwfK/ylkvmuXBf3pE\nTNW88tiYfMcIzlqS4wQjs4fDC9oRf9NQhzYsBAHa5GhQxSFIe/YnFNpQLnh212eH0BQriFZCDJ3U\nzdMgSpZQORDWkhynw8bfuYWOvBjL00IbpMiFIGC+NN3iEFSH/31Rn42iQZ49IsaaZ8tIf1Z+6LPW\nTU8SuhFHEC/UdaENquNCEDBfnLU4hODANCvNbYJFz+767BCaYEXRQKKhE9c8XcaCsoTCQVtZkuO0\n2+//8BQrC20wI/2XUEdpusUhXJ2sNLcJFr3STznLFCIa5vjLOT4ixppny6h6Vr7X68bkOMIn2cGi\nL9zRQt3JQhs8nAXR1oNYaByC2/CSNLcJFj37Z511SRR7RIw1z5a9JN/rnSzJgRhHftQXwdFbONSF\npm53Mz/elkMcFoJwUYt5t9u6B14gWHanM7H+N3YHP4hmWCaaG2GhE2ueLXtFvueAw+Nuz+daV+iz\nR004FjtiXW3yq4W3nQ2EQ2jvtCRL2OGg81ofOEC8Tt3Q0MRMVloF1WEfQutgVYzt2Sxhh2OGc67w\nsKXVWfH5TWcneDYWBxML4Orsp5wVHJmuznonpPxs2jkaLDZZnRVL3j0aCI5MV2e9E1J+Nn2lz2az\nCc/4YF02ClHBdNiH0DpYFWN7NkuY5ViTHOLzejQQJ1rT7Du1te3qbMsVWfb1zsaUhg88yJo4nAU3\nLLvTxYYyDW58MslhWMBkY3RQCb7kT5sbXCDBFwvJyg96SxtaE8Sr6tYpDQs86JpqOIvXpFUyHPsn\nkxweCxh1jI7jD+QpLngt8c5PvpWL13/B2TqlYYEHXVMNZ/EZZz+Z5IhY4OIPapzh+uyWSlu4krrD\njgsL6H/+4prrCWXL9mTZZxcpDaGGncxa4I4lPvndXb4l7OxOr4R4pZ+8Puuw4N4iNQIv2D8o/IM/\nPuajwCzXQrPyvd6qKWF26ewipaHOVmtsOAtGzEoL/PFY5NwYQvOtdlgyhvzpAb+Npev+e7mSetKx\nubOEnqNqSphdOrtIaWjgoVpjw1l8vbNujA6++zXL6MZ2f2XQ2XySY5nSkMBDWBOGs/h6Z91IEjN8\nZMMH2oP76mw/OhxzVpMccBSf94bGx+3Qv+tBK9g1DEHUa2w4i/8mZyFJC7cCS5+VO7nYgCFnLckx\nnU9wylHcEryQTZjqQSu4MgYe6jVuYA6slI1CMKD+9+zu2DiE1sGCTkEfVyfqqff9fjPzKPJ0g7Ao\nyRJ6DpfkwEeMYORJX40kB3143uVA5EMQuqYazoLRsjtdubng3zHO2SG0DhZ8gtEt8IOOyr/lO1Y+\nu5wl9BxuTA4e88e1DfFC3alOaVjgwdZUw1l8v7NsIT4Thl9nOuE5yrcEXDjmrN6geRz+mOSoUxoW\neLA11XAWrDYrjWvrf79ff63P2hgdNNLFhs/MN9yVX5Pv9VqS44JQJ4q7lZYt++ykKQ3+KuACD7qm\nGs7iNWnqKRe80tecPV758YoG2MAq32d0jI49XE440phdOAwUdxsCyHYMz2FJjiMdaOQsGQEbztqg\nFXUIwtbE4Sy4aVlpXFv/m9LXkxz8plYod1wzLG6EjdEBGZRrSbfRJ7tsnZVvenFLTXKQqfrhBGta\nzgrXi9OstAo2KtWVTTQMT/mXfiKVhR0sv0m33CKs+XDjDge9Df7UZ7vkz1e0pD3fgtZ2lDbRDv60\nBsfBlE+kwtPBSqhof4JVfM/0nuDb3IFOngobNuAdPU500wu3vlPssDfRqpbO/hMD4TtYHea4uEVY\n8dEGPY47pOH8jv4XnIVvgXu5wNR3o9XQqWrpeQtnOm6AkX/V2biffvU4O5/dSz+EL/gpfXZvm04f\nSDgLH/vQSfTSZn8v1Y1uzrcIqz1J23X0LjD/ep+l5+v5kWw6SlsNrfosn19b7PrfPRrU1v5tZ3UU\nRRUy7iw9qrNcEmC4DpZyPSu0duW7fTabTXimC9alohAu4V7gOuwVGh9abj90gCmHFrpADFezTVcH\nyyo8KUXCFh9vnOX4xSTHfHOvcmMi9DBW2ulnrS5UHQ3ovDy8CTpYT/y0VS3Cd/vsO3pMWe6SRvmW\n4vpZh73V0NpZvBmHt2L11cHS9c8KLcL/Jmd32HmPu4+cdcFouNNVTzrQta92lvMPVcoBf1o6z1ft\naWWIjriwtdOXnWSeYUB7vzi4YemLJlrsQ8fdHs11rw4W1ZAQCs6w/gtcQbjqPm4RRj4mChy8qPkf\n8VxdyT9UKYfwTBXA4cudVZ6jJa3JGRc22LFCE63VUo/WwYIqGkKh6qT/iG+fWb/qtwhbfI7DUy/K\nwVnLP1QpB/jtIvZXvl0DF3j4Ez4ubHuxIFwsUKXGjnVaDZ1OsYtmsaCehlBwG9a/pe+iP3LvqkXY\n4lO9C/a4IDiLpHyRo0o5xOiGDNFR5Tla0iJZc84rfe367BKug4X7yV0lK/r54u72UY4HWfmeY6nA\nlnScpQqWcuD6Et0Io13Iwk4vM6ZeySv9NWdDCKXoP1Mn+lectZQDu1KiG3G0C8tzZHd6ZfFfcdaH\nUKL+m1yJzMr3eqEp6SSH6zUu5QAIGt3wo13oQnQrKw3rupdX6tiH0DpYGEHDIzTdX/X68RKv3CPM\nEnoOvPUFhxkck2M+HOHl2oV4vq61zacceAOObqASGaIDipbnyEpz7Fhss4/tpw5WCKFE/XZfPCvf\nc0B2Az7Z6fQchqeBl3ubI57PJpizPuVQzMDoxmK0C81zZKMQlbNt9rFcSAdLzmrwmSpR/86syMr3\nHC7JsT3BC4Na8qqTHOas7U2oq9ENP9qFLmSw7E4X6jL1fcDYP9tny8gREE/x+uFMqFzOQClZ+V6v\nS3LgWenBj8uBeL6uti2kHOyZKn60iyrPkZVWHJVJkz3fUIGhaQcLvgTTWSv0Wa9/OuBXY+4feULP\nMelRGul5r5EQ+Nd1NqQcyoNW6JkquCV/U6jyHN/srIVQqOGsnL/j7suZblZ+cJau40PfJ1T+6kFF\n+Nd1VlIOnISooht8pK4WZqUJdZl6pfqOyXehgNbAKkkODaFQfdJ/ucG1kN1W7lJm5XsOS3IA7MaO\n2UgSnbX8A3/qYQ1OQvjohj6nxC8c9MIfixz7IJq12rBYv4VQ4HTmCkNr3SFDS6+3nLUkB8DK+Ru6\nVjvLywb/Z3d6BW9uhBVDaB2sANybyRJ2OTCY616xz7oVrxez0irkjtIhtA5WxdiezRJ2OcI1+LXP\nmstvO3vz51yrsx901qCohHuq27+run+Yze70CqbDPoTWwaoY27NZwizH6qz4vDorTrSm2f7U2vY3\nnM1mE1p63LIYhXArnhc77ENoHaznAsraLGGW4xeTHKn2dI/y2S4UWL6tz76jxzVsyIvVWedgr7g6\nG5xBO5Z9tgQqYlIDL57Pe463xBAEI37C2T8kOYL0xkxsiT4SBmq2RuKIjcjKjxwNEWVRy1mJc1RJ\nDbhJg7cm8MZPFYJgrKy0SoxTKsRYYwjNYQGEu3bfGomjakSWMHJUjXGzS2ctUFElNfzAFjEEwXhZ\naY4di6rUiHHxEJpiIYJ/ukprJI5EkgNh6lfgqFe6+aWzsLJcJq2SGm5giyoEwXhDXjhnjRjxhtB8\nqyVvQuJaI3FUjcgSeg7C7vx75ixtYkmNctMHo8BVCIKxs9IqJV7pJ698+7xJcySOqhFZ+V5v1ZQw\n+0dnLanhBrZgCA1B8GxWWuD/tT4b8xrdkTgmbURWfuVsOslBzdauE5IabmALqmUhiC90NuY1IHKI\nGhsjcVgjxpzVJMd0h/Mmd1O8ulvDFslxluYsqVE76z56qWZWWiGRie8DukvfPs7GvAaYChf8WyNx\nWCOy8r1el+S4Q/Du6O+EIV7IJmCDXQPhPi0/NQOq+YEtYLgrD4NbZaMQWNe9PLsjHkJTrJjXQLb2\nSByuEVlC5UBUe7oKjk0/uZ+5NZ6uAhVKA6ukRhzYwocgkGOwl/3OcTbmNVge/F+MxOEbMdRnbeiJ\n3W0zXVww5OnRoEpqlDNvHtgihCBYelaaNpQL/t3l+uwQmmCFvIbx1SNxhEZkCYWDYC3Jcbw9EiOo\nlwZWSQ049UY0+n1bDEGw9qw0rq3/vdJPOcvgrJ/KzZE4YiOy8r1eG5MD3rLnB3491Rfihbq4pjTQ\nkholDqEDW1QhCEbLSlNuLnj2zzpLeQ2W3hqJo2pEVr7X65IckInd3OQjCVvWcNZCEJbUqAa2qEIQ\nbFFWGtfW/6bUiIswrZMtGJbmNUoSpTESR9WIrHzPAZ2wPF3lSB9eeGFFXg1nZdWr06y0Cjcq1ZVD\naB0sBX1WyBK2OS50prRfnW1Y/J6z0xmPsT4zh3jtvdAgf74oK61C6bAPoXWwKsb2bJaww4EXgOc/\nfoK1qf+0NCutwukoHULrYFWM7dksYZYD8bJ124p0aVaabsCFDvsQWgerYmzPZgmzHKuz4vPqrDjR\nmmb7U2vb33A2m01o6XHLslEItwkWO+xDaB2sirE9myXMcqxJDvH5N/rsO+8h0QXTrDS3CRY77ENo\nHayKsT2bJcxyIF62bluRLs1K0w240GEfQutgVYzt2SxhlqPtbAlU+OeqFDlljQYjnMqsNLcJFoPS\nTyU5FmOKmGBtVFUnKz/ordriZ1vOaqDCnqtSNpE1FoxwWFlpbhMsOqUCj4uH0BSrHlPECdZGVXWy\nhMqBKp+8ls66QIU9V4UQbI0FIxx0VprbBIuq1OBx8RCaYtVjijjB2qiqTpZQOVDlk9fSWagsl0n9\nzQfGKGssGOGgs9LcJlj0SoUYFg+hKRb/NEufnOKeqaLPqfn7T1eBVkkDe87WwQg0aNCLX3KWBNmY\nIl5wbJTVye5K3XvE0f/3vM/CAFB+aC31vApGMHxWWiXGK5VdClWG0DwW/IzZmJxgfVgMrbU6WcLI\nkR+Tg+ikgfJcFVNY1sRgRFmdlWZoVPJKhRhWDKF5rMmPKeIEh0a5OlnCwGFjcsAPz+PPQxEv1oVW\nuQbCnpcfpqIPssYHI3A5vrLSuLb+9+wCP4rmseKYIpXg0ihfJyvfc7gkB95g3/sB7xAvZBOwxa6B\ncK/WD8ika1wwQjzKRiGkfpl6doUfzIUELB+rgP1uz1QB4tIoP+5IVr7ncEkOHFHtIqP9YNMQL+wF\nXFgauCVtNgqbrcFKGozgma/rs5YrEoGU5AiN8nWG+qwmOS407JcOWAeUiNdz1p6rotpcp4KsTfyY\nzUpTMC54dgc/hOaw4pgiTIWCfaNCnSyh44DxaeSX90d29k93GEsD+ahx8++qssYHI9SorDTd4Fed\njWOKmGDfqFAnKz84a2Ny4CBImz8+qaL4Z89VKUkOPU74YIQalZWmG/yqs3FMERNsjfro01WmA5y/\n6TMIsWVoR9gLcHTf/jxuNPzoZTdvORPKSQ5b44IRatTbzhp8EabQ2YJriT45haWbYGuUjTuC+Fn5\njgO2kiQHuAgfkf7DvuFsthl1vay0aruoVFcOoXWwFPRZIUvY4cBQMj4HRF+rs2LFm85e4STLhgkG\n0NXZDzl7gnGUfJBjdVaMHTzO0ZEWMgAAClVJREFU6uaLwtpnxZI3jwYCo9PVWbFidVacaE07n9ut\nqotlv+FsNpuwEBMXZKMQcas1yVH50ZjN7vRq004/G0LrYFWM7dksYZYD8bJ124p0aVaabsCFDvsQ\nWgerYmzPZgmzHAtn7Y6/lUCK3a+nEo+UPzb0QtWwoNQ4sw0NaIalqYKy3uaP+/28o+HgbBlWyhIa\nR6BezFTO2h1/K+E27n69Dh9Rjx+RlVZpcEo95xCaYWmqoLDZPF3soke/2TKslCU0jqoh1Wxw1u74\nW4nr2/36O/0obMZLi4NDL1QKVGnkzDY0oCkWPPNjK0P2Ug2dP9BjRC6oX5dRjSyhcQTqxUxwFtba\ntWcrwWILGFzpcjeNQTw49EKlwSt1nNmGBjTDihfloRdIvZmumhzx5p4uo3VZQuMQyPY056y7X8/3\nJ+hJcoNDL1RCvNK/4OzpgTdv6Z7tNzhb7nHMMLY7P4r79NBhbF8eeuHvOFtFJSxlAA8+uPPNcFuG\nkgb6bHfUCMFr9xrXf2Afw8kZHJge+3qQg4GHaCCSe7XZ0w11SP78MaQKoI6b3z3oaa5hGdR43Vkb\nNSKIKDO5owGYWkZeWDhrtz+z0ioZv+MskYSoBCyh+ftuc3786Ci8WicrX/Xa8z9gl/I17z38Jkya\nl3VW7tdXw0f48SOy0oS6TFUpzLv3yRCaxwK4EJUo83c8ut4f+ttjrZMlVA4bNWJTRri/g60HOXoj\nnr9uYG2zEkjhF468wMfZQ3kWqh964QPXDRznEJq2JKQKQLzNs6Ub0G/LqHVZQuXQrAFsfqAgEQwe\nAadRZTAZxNO9AMutbVYiYvyH9+v53c/DR4w9REPhuNBmTx/2Appi+VQB1tB5yRvBHVRdxhAv91nN\nGsD25OyRbjDSA51gEeKpHpg3P62E1X7gTgTdmnTDR5QHRcoZeVYawPlXm/1NZ32qAMlsHhsCLxj2\nxZbRkqx81WtZg+LshrrrrTzR8I/OLkZe0OEjpsGhF6gd9k+VwiK3N7MNNSAoKZalCli/zZ/o0gQ+\nvtSWEUaWUDkoFHbkm1/UZ4uz5Xw0Omt3/K1U36/Hk7gtXZEZHHoheOHccEEHqJJtaEDTVoNtMSph\n85vtrlyRsWUIkiU0Dssa8NGA4wY/TWeDzFdnstIqXFMaVgyhdbACcG8mS9jkKMdZ+gQrJ3SI16zb\nE9BfnpVWIXTYh9A6WBVjezZL2OQo5wZoqpwbrM6Kz+85S1cAMZkczmebe0EI89OstAqxwz6E1sGq\nGNuzWcIGx2V3e2zx8DrDIVzGRkS8Rt02+fOlWWkVSod9CK2DVTG2Z7OEWY7VWfF5dVacaE2z/am1\n7W84G37T0CLNLcv+hKJC67APoXWwKsb2bJYwy9H8BUib+k9Lszu9wun0syG0DlbF2J7NEmY5EC9b\nt61Il2al6QZc6LAPoXWwKsb2bJYwy7E6Kz7/rrMWzbASMbtRFuIai16kv3hLU8o09AGDyzY0oBlW\nTGngPZEjXJrD70iq37UJQbKExhGoFzOhz9pTMazEW9goC2GNj16kpVUanFIPl21oQDOsmNKAC7T4\nBFY8mTf91ibCyBIaR6BezARnLZphJd7CRllwa2L04m1nI1y2oaFF1uqY0oA+K8kO029tIowsoXEE\n6sVMcNaiGVbiLWwkhrjGXVB921lgcnDZhoYWWavlbpSs1nnTb22iSllC4xDs9jQ4a9EMK7nNaJSF\nuMZZ8V/ibNQ//dWRIzSaYc/HYHvdr52kzvc6201yUFtEvxtdYqDPvpTkgJ+L6Z0zK5EaN8qCrvla\nZ11yg8THedXv2vS6sy7JIU9rlylwhqMB707tmhbSwBVy1xPLuuZrnUWVcdALPy/6fZtedtaSHJI3\nkCmRL5y1p2JYiWq6URZszXc7qykNaoAlO1S/a1P6Y0I/wSzJIXkDmwIhOvvPf/6nUMOns36IWolX\nyn72dfyH+fR//1GUVwqBXW51TGNoilWlNGJqw1pmbcoTKscyyRGcRTt0L8AavAlPT1WxEttkoyz4\nNd/aZ6uURkhtmH5rE7Tx5aOBfh7BxnwfzKYFz5y1p2JYiY0F+zn5EOMb3+psldLwqQ3XMm0TtvFl\nZxdJDgARhwueOmvRDCvJoBEyEoOtQTnf6qylNOokh9cvbcKWvO7sIskBID1nLZphpfIMDR1lwa2J\n0Yu0NGqH/dP9GuGyXciAoGRYltLgJIrOe/06ugSCZAmNo05yAEjPWSQYf2WlVQymNKwYQutgBeDe\nTJawySGOyhRIEK9ZtyegvzwrrULosA+hdbAqxvZslrDJcaC8AfTZMgWK1Vnx+Q1nJW8gU4Jcnf2A\nswIRpquzYscbfVYgwnR1VuxYnRUnWtPmp0urYmPZ6mzDFF20OqtWQKHjRrYLeageVqjTm8kSdvQu\nYBEvW3excVyQlRa36rEPob3TkixhlmN1Vnb07zqrKQfLOxTiyzzPVxpd1erQKr50hMWstIIok9AH\nPpbkgJQw/FKThw5BJhtTJK55XX7QK41oTEOftZSDlXgbunZ0wOFV6zV28fh9Zz+Y5DieIQ9DDz2m\nBrgxReKa1+UPOWspByuxs3v6NegDfgNZrdnbg5DfdfbDSQ4QftVbejamCGQ63JoB+UPOWsrBSuzs\nCX9cTZc645rL/azq33UWmNzl3iE0bfXdjwQLwDamSFgzIl852Jnuf2yA1rWUg5XclnQ0iGvm6Tud\nvdo1J2yAG1MkrBmRr245Z1rF4CxX0JRDneQ4l/s1UK3UuV++1NnHeTPveewoahTfF4MxRaD32poh\n+c7ZzyQ5NvPNjC1JiCOE+76yz8JxC39XLD/bdmOK+DVj8s3ZZZJj3l2f/GrJPi2txJ15/pGzmLIG\nPwu+1Fm6vb7T8SF0TJHjw9aMyVdnl0kOfEYz348HYxZHA005+OE22NnppygtdQ5o9Fc6O+FZDDYO\no8j8kmeA2JpB+ersMslxk6EKkLJ21lIOVoJqFEGYrjwSR1lzpJjkdzrLx9W7f74BtALHFNE1o/LV\n2WWSA0+gLrI3K2ct5WAl3AH8HrqSUllz2OLr8UO/i4Q6CDXwUqWw7cfOuvgMwPdZVIZjiuiaUfmq\nt5XkgDthEvOJzlrKwUpkFn8fuOFZYnMNVvomZ7l9Oz33sjFFqjWvf9FRZ1tJDjhrKsMbxKOBpRys\nxEmIO37UHjDob2vI8unxld8Upiscqo54bsD67Rkgtob0vy5fnW0lOeA3zWxLdZy1lIOVSpLjtNvt\nrvh54NbA3Hx9POhxIVB+u8/aOCCjaNpqULbb8cgh9ZgitmZMvnE0khx3NTb2WWjQG6+3nQ3cQ2jW\n6oCVmskSNjlKguME75RTOTVFvGbdlJpQKSstbNRlH0J7pyVZwiYHJzg258PhgOlyfK3Osg/5g1nD\nWUlw0I/O5AxpdfYDzgpEmK7Oih1vHQ0ExE1XZ8WM33EWf5Tqv6wJW3Z6IwA5wmS34nob3lZ/H4GX\nKOn1GgzVXmKlQF4iTHAoHoyDjy+7epGS4yvB5Q58+UXp8pJ9HG2JlZLxEmGCQ/D+H+0TYvUUV6rk\nAAAAAElFTkSuQmCC\n", "prompt_number": 12, "text": [ "\u23a111.416\u23a4 = c\u22c5\u23a11\u23a4 + b\u22c5\u23a10.389\u23a4 + \u23a1e\u2081 \u23a4\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a24.514 \u23a5 \u23a21\u23a5 \u23a2 0.2 \u23a5 \u23a2e\u2082 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a212.204\u23a5 \u23a21\u23a5 \u23a20.241\u23a5 \u23a2e\u2083 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a214.835\u23a5 \u23a21\u23a5 \u23a20.463\u23a5 \u23a2e\u2084 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a28.416 \u23a5 \u23a21\u23a5 \u23a24.585\u23a5 \u23a2e\u2085 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a26.563 \u23a5 \u23a21\u23a5 \u23a21.097\u23a5 \u23a2e\u2086 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a217.343\u23a5 \u23a21\u23a5 \u23a21.642\u23a5 \u23a2e\u2087 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a213.02 \u23a5 \u23a21\u23a5 \u23a24.972\u23a5 \u23a2e\u2088 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a215.19 \u23a5 \u23a21\u23a5 \u23a27.957\u23a5 \u23a2e\u2089 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a211.902\u23a5 \u23a21\u23a5 \u23a25.585\u23a5 \u23a2e\u2081\u2080\u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a222.721\u23a5 \u23a21\u23a5 \u23a25.527\u23a5 \u23a2e\u2081\u2081\u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a322.324\u23a6 \u23a31\u23a6 \u23a36.964\u23a6 \u23a3e\u2081\u2082\u23a6" ] } ], "prompt_number": 12 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can do the same calculation with matrix multiplication.\n", "\n", "Call $X$ the matrix of two columns, where the first column is the column of ones ($o$ above) and the second column is $x$. Call $B$ the column vector:\n", "\n", "$$\n", "\\left[\n", "\\begin{array}{B}\n", "c \\\\\n", "b \\\\\n", "\\end{array}\n", "\\right]\n", "$$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "X = o.row_join(x)\n", "B = Matrix([c, b])\n", "Eq(y, MatAdd(MatMul(X, B), epsilon))" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left[\\begin{matrix}11.416\\\\4.514\\\\12.204\\\\14.835\\\\8.416\\\\6.563\\\\17.343\\\\13.02\\\\15.19\\\\11.902\\\\22.721\\\\22.324\\end{matrix}\\right] = \\left[\\begin{matrix}1 & 0.389\\\\1 & 0.2\\\\1 & 0.241\\\\1 & 0.463\\\\1 & 4.585\\\\1 & 1.097\\\\1 & 1.642\\\\1 & 4.972\\\\1 & 7.957\\\\1 & 5.585\\\\1 & 5.527\\\\1 & 6.964\\end{matrix}\\right] \\left[\\begin{matrix}c\\\\b\\end{matrix}\\right] + \\left[\\begin{matrix}e_{1}\\\\e_{2}\\\\e_{3}\\\\e_{4}\\\\e_{5}\\\\e_{6}\\\\e_{7}\\\\e_{8}\\\\e_{9}\\\\e_{10}\\\\e_{11}\\\\e_{12}\\end{matrix}\\right]$$" ], "metadata": {}, "output_type": "pyout", "png": "iVBORw0KGgoAAAANSUhEUgAAAUIAAAErCAMAAACciHo8AAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\nTlMAMquZdlQQQO0wRM0iu+/dZol8bIlZUvoAACAASURBVHgB7V3porM6jiRr92Q5J5nJ+z/rSLIt\nlbwQxJfcvt0NPxIvslQUBAiFzLR7ybKf1i6/ycHi4X8csBspiqLrBBoXoNSQu9fhSMsJhseKDx6+\nfy0e9McBu5GiKLpOoHEBSg25e+1g5NriI0LhJwJ2gQZQdMdD41JaOKSzvZ6LFyvllt0hF7Tn+Tw+\nCxsB8IOAJbB9nw/P46EEoObb8/m8y6/lfDw+Dxe2FJtbGRNAUYaMvh3KkRG1VxTe7q9MoZV09O9d\nitpzvtI6HK65PwAewKkzjYKFK9F1/lV+brwNT68HfcrW3HP3D1VvP8UmgAID9cqAstetbY7C8/2w\nTxRaSS2n4w9TCD3XJ9Xvide0Mcx4tqTgwFlvwOOHW59lI03HF++Br1/ikWmcbtSzF5v9ChTsYW5R\nlHNG1OcopPqx7IVQSi5uj2vGmW0er7LpxSCw/RGcBUxh8PMuFwknhXT5od3+/CLingLlTGRyjerF\nJoACI/XKiLLXX9qWU/icKgrvsv2LI9kYWpkvILg5Cl9C4UX2PfUoP+TL604HnOeDCBUyjy85LoZ+\nC+pxUECUAxNpXkzh41ZT+Lru6GCfoYfAI7gZCs8vOeJdXkdcg2s6DL5eD2Jwmn7yXii/7BAKdNor\nI0rqPz+Px4P73eVBSyk802r4vZB+ULIKcoIkb4GfEIKbofD24oMtXTDIlwDePX+FQTqfvF5yYfCU\na6lroTmAQhzOfCDKabrwWe2WAJUDbxq8lEIeW1Mox5/DTz6FB8AjuBiFhOL5w1c5j8Pu+uJD43Sl\njXt6ymma25dfnZL17IIo6ZzPZ3++wNod9nQIhmUhhSeG7SmUMyNjjh+FENwMhee0+1U/5OmHNtqD\nj5KPl2w++ouwe5bj5bco3P9QmEfaW05rKDzLDlxRmI5CjxXbfyGFUzoWngo9000ORXeqp11/pz35\nWux7e2E6s6V9bxWFpz0vr5+9HNjzjpPOyF/cC6d0Lf8oFyzTS0r31+Ocf677fBKZruUeybf2Qjge\nT6soTOzLpTUVM4Xpeu1QLm0C4JfuhUfxrf9/pgTgl65H87+R+216cHGnl6gBFGmlxp+IcvpJpzTZ\nhaIUypVsipMuwZTC6U5uz3yYlSUAHsHNHAun6ZeP4Xz6SCgevNOd+OR7kbPxg0pHKp3l2iCKIsEe\nfyJKCkOG53RRM0vhYf/z+pVfh5V+82/keX+9uMt6pufhsOfzjCyrKARnxQ9+87VYipBQXA6Hw11O\nX7v9Id9mIBCpKYoCA/XKjsLpeDjKFiTLWQp7npa2raJwqfPFdgEU73x6CsF6oxDImC2OKSyH/zSc\nt9rQdjZC3RnY/p8JWAOQegBFdzw0DlDeDr+vdGWSbTcKgTRfHFDojai2UdhQUho2CgsTq783CldT\nVwZuFBYmVn9vFK6mrgzcKCxMrP7eKFxNXRn4hxSq6k7+rscz3arNN1ZVhefbtyWYfAcuah048Oj8\n+UojytM/47v+PzfbAAob1C85lH0TaR1cF2bVXUx++KEjucszOeEcbcgwAB7AOY8zOCtR/sz3Pp7l\nhjmMC6CAUd0ioOz2l8Y+hUl1zzbX5/4ot4u9cO5syDQAXsF5jwVT57sW5fdyC/7H/1nlcQEUnTCu\nSVG61rbSpdBUdx5Qbghz2e7ueRvqCoBHcOaR/Y+WWpT/Efli3z5PFUAxilXaEWVp6313KTTVnYf0\nKfQ2ZBYAj+CWUViL8leRAP6+FILqLhQ+Hvp4la5wZfNlCvui/PSrqgrjlCWwIcuQ0TduaLJZLMWz\nLZ06iljH7vnRi1s+chcKaxsyC4BHcMUjRxouHVGebC+oCeWxARTDaLkDUaIU/zzcD1k7F0sO6W0n\nVN01zD1pp2WFOzYB8BiweNRIvUKfwuqSQNen52FFG6IEKf5Jp9a9PIWSfbYUOtVdIx+SRpZXuGfz\nTQq7ovzBP5eRoAZQ6LoNCo5Ck+L5wRBTDGksh/zH6x/mxavu1J4Yz9dgicLGhof/b9Z2zdWwhAEX\n7YWNKE+uj3ia01ABFDpmUECUk0nx/BRKOa7JSA7p6PaqO9mkhxYO6cidVrixYVeB7Y8Bl1FYi/Ik\no/F/o/R4A0cvSwBFGTL6RpSTP+zaI480mEM62+SwqO5US3/jftOvH1cYbHhQADwGRI8pdu+zFuXz\nv8t0yY8DAihwWK+MKJ0UT/p2+rOWRg0oFNU96d/ygLM8M0AjcIWLMp/DB8AjOPTYW5HcVonyt196\ncL1+xIpNAyhmokkXonRSPOnnOLZLYVHdk/59Ozz36SyOwnmxUWcB8AYOPaqnXqES5a8pWcg/YsXj\nAih6YbDNUHIrSPEPx+Doh4yuFpYD4D24hf6XmQVQvHM4QnmhK7oL3CTikCPbdzF8fwD8ZwL68LkW\nQNEdD40DlLvr6XTivxtl2SgsTDTfAwrl3h9ewG0UNtSVhgGFpVu/NwqVirqwUVgzEq5vFIYpqwdE\nKDyTvvTnyw2v2OfdfSZgN0YARXc8NC5FySGX0g3uO8XA5cRnAnYwfPHSuhtNGrfTyZCbpRt6o/Dz\nFGZl/FZJ3Vb3orjp8Wt/yHNSvI9FN7ccKo0NLgIohtzljpV7YVHGa6kb6l4Ut5vvAfAArgTsr4+P\nNQEKts+xnYsAin5IawWU1tgpcUi1NWW8lrqt7kVx0OMD4DsBO9Do3CBSu2XFGwo2b7L0uTGAgs3n\nFkU5Z0R9jkKq59t3tdRtdSeKox4fAI/g5u4XulgEzlBQxWKjiwAK8jG7IMo5wwGFtdRtdSeKox4f\nAI/gcP1roC4WdRoKqlhsdBFAUUer64iy7sP6gMJkUkvdVHeiuNPjA+ARHK4/AqOyi2V9CRXERhcB\nFOaxX0KUjCaeFd9K3Sx9o6Lr9fgAeASH61+tCsayriTAY2x0EUBhLvslROmk+OMzHZTzMA6JtgjH\nzrbJmOu4Wl6PD4AfBaxWBWNZV0KFsRFzAIW57JcQJUjxB85DRAl7hkJnR1GkDqJ4pccHwCM4XP9q\nVSCW9SRULja6CKAwn/0SoqTpcEpW/IG2novCFdScDU4tdee6ZqrXenxABO8HbFdEY1lXQuFjG+bQ\nAwHmtF9ClCDFs7GlSFOlluIVTi11l3r6HVGmeq3Huy3TB1VacftqwNIJ3xpL2zIKHxtdBFCo00EB\nUXopXlKkdRSHRNsCJ/15Mqlb65Uobnp8AHwvoCKCQhWLkIr8WFCV2AUzjwyggEDdIqJEKf585ElJ\nbBlQaFJ31uNN+lZRXJyYHh8Aj+Bw/Q1WLmmsBgUblNjoIoCiiVY1IEovxc+dTlQZB6lb9HiomyhO\nl7c5U56DB8AbOA1Ywc9Vi1WjsNjeRQBFP6S1GkpuAymef7dlOgXq4ZDe1nzESgHwnwnYhRdA0R0P\njX2Uchy8lZll2HyjEEjzxT6FEz9f+ChTRG0Ues6q2oBCmuvtuZckkmy/7YUVcVYdUGgGG4UNF1XD\nRmFFSLy6URjnrBoRoXCp5lyFqKoBEfwzAav4qRpA0R0PjUtRblI8kOaLkb1wqa2PUNcCF7WfCVgD\nkHoARXc8NC5FuV3UAGm+uJpCkLVV6ibXNiG8TRLv0uXX/UcmzxDQrwJHrSaoZ4OCqkyOb3i491+/\nFzpZO92uY2A5+4QnhIdJ4iFdniwC4GH7uoASCT8qKV66EiqdHB/wcHcABQbqlQFlr1vb3A/ZpHju\nB5kdJoS3SeInTZcXbwHwCs4HVFClUEvx3J5R6eT4gIe7AyjYfG5RlHNG1OcopLrdezOpm5ptQnib\nJN7leofAIzgL2CKtpXiyyKhscnzAww7+VhSa1E3IdEJ4hsmTgfKUoD6PMAB+KYW1FE8xM6pqcvyE\nh5EFULD53IIo5+yGeyFI3Tx+XyaEF2d5knhIlw+BR3Aze2FHii+oqsnxEx6G9j0Kw1I8St3CWpkQ\nnk6JZZJ4TJcPgV9IYasjF1SkA/DPIE3Fq3i+SiFMUE8TBQgl+WO0F6LULdB0QniupUnixUVOl/9r\nKCyozmmGa50c3/B8ay+0CepprV0u7eh04qRuGoQTwjNzPEl8WnK6/DcobKR4Q8WvjeGQeQ4nw/Mt\nCk2K56dr3LQaHBI153Ro8lI3gbUJ4W2SeJcuHxLB24Bpc9SflRQPqGxyfMMjowMPBNTR6jqidFL8\n8ewo7EvxXuom3oloXmhea50k3qfLf2MvrCeoB1Q2Ob7hEYTf2gshK3538hRyyNHRvUjdhM0mhE+N\nPEm8S5f/CoWNFM8sCQCbHN/wfJVCmKD+SXtUOYxxzDkKITXeJoS3SeJduvxXKKwnqBeSkgCvk+Mb\nnq9SaBPU03XVHIUoaxeZPaXG24TwNkm8pcuXjSFr8f7DdnsM2BlXSfFkUVDZ5PiGhx1864esUjxf\nV81R2FmLxU0B8EbhYu9LDQMo3rnsozzRawYOL37IsCz1D7m0x78D4Pvg4iE7IwIoOqNd0xjlbfZY\n6JyEKgHwY3ChiD3jAIrecGwbonzcX3zbrywccmhbjBZ9B8B/JmAXVQBFdzw0LkW5UQik+eJGoedj\nRW2jcAVpfkiEQvdn0LsJ1AL/Tj8TsIstgKI7HhqXoqz/I4OLYDFwIF+6fYMI2DyA4p33pSi308mQ\nyY3CITVLO/6AQlPdLVjRv004L0q42gR+Qg7cnBTPzjU0lU38N5DWxtYBFGw+tziUM4adH7LcxMLL\nbxqeVXkVzlUJN9cB8ABOPZqjqoQPBJj4byCtjQcGUFRxmiqgbPqwoaXQXsNudk0CuirhZhMAr+De\nSPHsHB8IMPEfQK5+IMCg90uKst+trS2FprqrkanyWbI0JVxtItsfwc2IoOzbQnNNlWsAqW1sEEEh\n9jMfiHLGTEJ620p158GmyucVrpRwCbBmL6SBbyi00BxE6QKQ2hZFIfYzH56WsWG7F1aqOw0t+jcV\n8wpXSri4/wqFEJqD7FX8t0cDrI0tAijYfG6pKFwsxbNPU90lQtG/uZIodEq42NBHADyCm90LMTTH\nAfFfQUJbDAU7nFsQJWTF7/a7W3nBqgzv7IX2GvYUoOjfXCsUyj1HVcLVUxrw9hPBzVKIodWriP8V\nyDUPBKjDQQFRQlb8hWbidQcPptA/l12r7qZ/U6i8wrUSziACD4pjwDkKXWhdURb/a5D6QEAAhToc\nFBAlZMVfTjvU79KKO7pBdRfXoH9TPa+wKeEW/vM/ZB+aIpn4b48GWJtACaAw6P2So8V2vEt5gqKM\nan7IoLqLDejfVM8UmhJe/HzjWOhDUySdKx9AaltC8jUK+Zgiy+VwPL55WYep7mUQfRdVPlNoSrjZ\nBMDj9p37IYvzEpoqJv4bSGsT6wAKg94vIUrIij9TCuMJH+1q9kJ4DXtKRRf/dQK6KuEWPQAewb2l\nEB4IMPHfHg2wNoESQGHQ+yVEWWXFn9OrS9K4lsLJVPekwpv+DcK5e018FLyBA4/99SjSe4Ji4r+B\ntDb28C0KVYqf7rQXnvWRshzS1qi/GstaA+A/E7ALK4CiOx4aByg5pdsJyZ29ELxEigHwA3CRaCPb\nAIqRi9I+QPmkSxr3qomNwsJY8z2gcKK3hOjpmQdtFDbUlYYRhaW/fG8UFiaa743ChpJow0ZhlLHG\nPkLhUs25CeIaAiL4ZwK66KXSoBBpxQ7/T37jUftC2zIcv5ei/I+X4n9Oux3cWjlT7bSMwsheuNQW\nN1FbDlyRfSZgC4FaGhQ/8JamPGK3UdjlLjf+XSj0Mjtq3aXHp6IL/Ab8eE3dXvhGinfdliUvJZ72\nafe7P9A78Z75cNeg+JfshbXMblq39lSp6H9A4RspvurWLPnzDz0vfuNbXo/0RrxyS/TvQWEts5vW\nrT1VKvpqCt9I8XW3Zcnv5XC2p9t2B94VTR79W1DYyOyqtVhPlYq+mkIa+OZ+oeu2LPmXvKKU0wnT\nD/hZThp/CwobmV0prHosFf2volCz5M/pFvCx3La76OSqfwsKG5ldte6qx1LR/yIKIUs+KSbPMpml\n3YefoZBOhUfZWddd1CyX4luZvWjdrselov9FFEKW/FNSU6+vtPc9LRdpSOGNp9a9yu9/FYWQFU/b\nwgJ2bnbVCedCDz1H8Et3u70Ab6nofz2F9KQhiUB5Lyx6HsMYUiipnOmRvzUUQlb89ULXApmWEtJp\nzlNPZueHRG5Nj6XGs6eACI4B3fkCgOUidrsseZovfvdME/tKlma2b1Dk68L8lnGxWkghogQpnud3\n5wdVdOGQ7ko3i7WPcphB/Vt+BtxTpaKLt2b7a4ymgAGRo8awOmFXWfL0AJU8V4APcTYoMoUgpU4L\nKUSUkBWf5tQAqBzS2U61zG5at/VUqejirwEPUaoiBoxQmLiiyfGzv6tcLJzLZTW3NigShXkHTuNW\nUag3e16nt1J8LbOb1m09VSq6IGvAJ7y9z7UUWpa8vHh9l9TcXdo5U6AGRd4L849Ksg/XUKhZ8TeO\ndsHXg3NIXCMCojJ7kuJB69aeKhVd0Dfge+SlNgy4bC/MTwXohPVH+ut8lrlqeMqTtJX7KMqxkHfY\n82H1RY1mxZ/l+FsuSDlmh0JLOK/1b+vxqeh98NLa/TAK30jx1p2gWJY8PQpwz88HycmyxGk2ZKaQ\ndDd6vDNZrdkLTYoX9vRMQR57FBY4se8G/Hi4UTi2WdnToCgUgr8OhfygTL0MUMq+/24vrH0tqzfg\nx8MG4MYDlvc0KJZReMKrlBxtgPJCf4VOcm2S7TjkwHY5brFswI/HfyZg13+D4uMUTo/D883DcV1k\nCxob8OMx/94U1uv1vb3weYUFzpmEYKOw3gypvu2FfV4Crf/VFC7VnOf5bETwsflnAnb9NyjenU7S\n8eb3R4467nizFOX3pPjnLyx651vWezsWdjd/+we/b8atG4V9bv6rj4Vup0jSdyVw0w1VErvvOstS\nylT3bWspdFp7u3Wq7vIsgOzMcujy749vUORj4fme3k8t4/Cec4oY+HfSQOSQQGGRviuBG6elZxdy\n467S4xvwTSxt6ATUPl8oeHKrPgsgdUFRvT++QaGnk3RzkQd2/iN/jEKTviuBm9Re3gOTJkC3kuV+\noW9bdSy0gMJJ/dF067MAbJlQ+PfHtyiUQnv7X4fCC+pJGQZs6BqYq1d7YVHGK4F7ctp7zlR3bZ37\nxS6OqyC4ZfcL03B7FoDqGYV7f3wHRaHwove6e3uhg5criLLXX9oGFEq3CdzJumjvmKle2jrgS4j2\nG8FFKHTPAmQU7v3xHRSFwqfdaO7shS1Gd3zrdWvbHIUmcCfzrL27THXT45ujkMZoCmspxGcBHIpJ\n32rfoCgUXvckeaTz4ToKF0vxsFOAwE0kqPYOmeraxiQ14BvmtGElhfgsAKAgt+n98ey/QZEppLEk\nWmWxpT0jKzQoIErIin+ezrSAHYdEW6MQBe40IGnvfJi86g5qenwDHqJUxX7AyihVDQ9nvtmzAB6F\nyaANikxhks7S7NRr9kKT4mlWCFqUgLTVUHM2yChw57Vj7b3OVFc9vhHBu5RIYz9g197wULc9JeBR\nwFvtGxSZQnk1eZZCF1KIKEGK319o4YdkylJL8QbZtizZqvYOmeralnw127+EaL9X7oX2lACgIO/4\nVvsGRaZQjtmXdGGzkEJECVI8H09PmBnPIdFWKXQCt01LD5nqlR7fgG+pKy3dgKXTfysebtZnAQAF\nrQ//SUnbc3wsFMktPxWyikKV4ilY2oKMiJchhU7gzlnxPC29LHJpXenxfwGF9iyAoajeH9+gyHsh\nS27ydDGNXEOhSvEcOV3OCwb6GFJYBO4kglfauxxNq7YGfAnRfsf3wizF67MA4lNQ2Fvtpa1BkSmk\nFzxM93w7cA2FKsVTlB2cS6jqKTTpO03MwqiSCI7ae8lUx7bOT4gH9xejEAL2TK07obBnATRVH95i\nzx5GFJ4PR2ZRljUUmhRPkfFHXVOYY6z6asCPvRiFY5uVPQ2KvBeiu1UUggN+thEWvxdCR7jYgB97\n+DenkBM2YNkoBDJ8cbihf/GSZvshe9ZcbUihs9oorOjA6kZhZmNLpsXdYkm5OamtT+mO7IVLNef5\nVWhE8LH5ZwJ2/QdQdMdD41KU35PiAUxdXLp963EL6s1euGDMwGQpyv/4i5oBPwuaNwoXkDRv8mcU\nZv3ba+38xPzzmB6a9/q3QAn8hBy4SmtvVgsnqO/lwnsoARRNpKrBoaz6sNr7IRf9u9La6X4/3W48\n88xplf4tDgPgAVyJhZh8GW799nLhKygBFD5MWwOUbSe0tBSa/l1p7RMkpct97x8n4wTAKziLBZB8\nESeo7+XCv5PivbdITVG+GdRSSAPyreJKa58sKb3SvyXGGgot1gimm6C+lwtfQQmgGIUs7Z+gUHyZ\n1p71AxYgKv1bDAPgEZy7s1/Q2zfK/t1c+ApKAIUF6ZcQZd8itXLIxhZXy7R2SEpPY1X/lmoAPAbE\nWC1QL7gPc+FnpPjW59IWREljFkvx4l9Xy2ntkJQuVqZ/S/UbFHrBfRrmwhuUAApBPfPhKbSs+Add\nloAGWt34zw6VQqqb1l5TCKdKHhcAj+AwVg5vX15w57fukGDXyYU3KAEUFqZfQpQwQf2DUiHPqJ5w\nSKc5szu3Wqq1u6R0ypNHNzSoEcH7wLgVA7pY1RAvuHNnPxceoARQVMGaKqIEKZ5nFJ4wgayW4sVR\nXq1Ka59cUjrq3zIosP1x+85Q6AV3XcUmFx6hBFCow0EBUYIUf/jdTTd8/J5DOlv2l1er0trdq9ud\n/i0YAuAx4AyFTnC39axz4R2UAArz2C8hysneFX/+fS2YDjevVqW102OlHEueTkn/uo6YghoAj+Bm\nKExrllBIuZsL76EEUCT3409ECRPUT4/ri/+h6cIhnS335NUyrT0r4ZqUXunf4i0AHgO+pVAE9wSg\nlwtfQQmgUA4GBUQJE9TTM4O73/SmizSwQ6Hp36a1V0nplf69nkKLNViPIvsnAJ1c+ArKtyhUKf4s\nZxK+W1CWDoWlK/gdAO+2bzDMG/MAijee2h+nDLjJhchxo/Adfdw/2NBXPg7ig0nbXjhkc0Ah37N8\nvj2dDL3OdQR+QgNwc96X9gVQvHO5FOW2Fw6Z3CgcUrO0Y6NwKVNDuwiFSzXnYTDpCIjgnwnYhRNA\n0R0PjUtRblI8kOaLkb1wqa2PUNcC58LPBKwBSD2AojseGpei3M7IQJov/hmFWR7Hqemz+9yjmrhF\nDWx/B25Gim9y8y2sQqtsAigMer/kUPZNpLW3F6o8blPTZw+lxzRxcx0AD+CKR/MDpTo3H8IqtMom\ngAICdYuAsttfGlsKQR63qenF3HpMEy9+1mkn5tH8QKnOzYewCq2y+VtQSKtQ7uHh7e20ZrnHNHFb\n4QB43L4lljmyUsrv0MnnJwir0CqbAAqL0y8hyr5Fam33Qmovq6U41UPqqTVx6Q6AR3AlloaoC5ab\nj2E9NLMJoKgD1XVEWfdhfZ5CfTV7GZJXuNLEpTcAHsG9pRCEQgirs+ZLbLMJoCjrNPpGlGSzToov\nU9NbkLzCXhNP3QHwCO4dhZibD2EdNLAJoLC16pcQJWTFU+apzyDjkN6W/LnV0tewc6DSg5p4BhAA\njwGLx+yl/vK5+VXYDA1tAijqUHUdUYIUz5LrEedH45BOc2ZHbrX0NeyuBzTxHDoggmNAFyu7gq8q\nN9+HzdDQJoAConSLiBKkeE7LdW+mnZPiq9ewcxy3wlkTz/ED2x+3r/PYroo96qF9HNZBQ5sACnU4\nKCBKk+JvL54b++2Mwnm1bGp6DeJW2KZw4v4AeATnPGqcUvC5+amVwyI0ZxNAUWKMvhGlSfHnROE7\n+SmvVvrB/+LzI7kHNfGCIAAewc1T6HPzLSxCczYBFAX46BtRghTPr0nYvZ0aPK+WTU2fpXj9IaMm\nXhAEwCO4eQp9br6FNWj0qNfsBPUFXvwbUYIUf6IrKH1jEnvlFXe29MTH/uf1Kynk9hr2pIRbD2ji\nBdsqCs1jceO/dfL5Woo3aJa/z0MDKHyktuZpOR6O8nQH0UVvWimTVJSQ3rb1tawlAP4zAbuwAii6\n46FxgJKfzuQZ1nXhkANbtVlWCID/TMAurACK7nhoHKC802Oid3ljRbbdKATSfHFA4YXeFY9KfO9Y\n6B0trgW2/wDc4lAzhgEUM16kaynKbS8cMrlROKRmacdG4VKmhnYRCpdqzsNg0hEQwT8TsAsngKI7\nHhqXotykeCDNFyN74VJbH6GuBc6FnwlYA5B6AEV3PDQuRdmckU3WtRL5NQVXSnJpuTqZ2oFzYWAN\nuKhqcW63+vl4fB5k1iJrY6N/PYUm61qJkYGCO8rrDoAHCn0YDoWLqsW50epyq0bekmptbBRAgYF6\nZUDZ69Y2txearGulZGkK7jCvOwBewdVhFFYuqFpc108y89iN78N5mwCKOlhdV5R1R1V3FFKf3Xuy\nEjWbgjvM6w6AR3AuTAVu8lKnvRF+esr//DNP9+1tAijqYHUdUdZ9WF9GISi4+rLxOq87AB7BraPw\n8mIJT5S7fxMK8412SmMd53V/hcJKyDb1mGaXfiTt09p41wigwD2pV7YNPZSQZdiyvdCSqeuk5FXJ\n1AYOjxyd1XBqMfVD/fCS+dJdG1l8g0LL5u5ATCH7a+R/YUXBrSk09SwAvh+wi495sxcciInUH4fd\n9fWjczGqTQDFKF5pLyhtYnW6tZrutR4p9aSYyVYrttxoxFlJjLOCO87rDoDvB1RMVcEJ2dTH9Qcf\nAR8vzSZUmwCKKkxTLSj3P7TuD75ZvTvsZXM+iL+THoM5JP4ZNOKspM5ZwR3mdQf+nfYDahgtOLWY\nWq2euNvRZCXWJsMCKDTMoFBQ4kvAT0IhJXXTRUqZn6H+j2zEWUlDsIKbfrnpZeNrk6nL9mXHnTAa\nD9VibtR6eQSERDJtS6O+sBdaNjfte0zhWcQnfWc8h+yvkVs3U3CHed0B8P2AiQP8RLWY263+k269\n0zwJ1iYjAygwUq9cUOLE6kLh7Cs4WwAABvRJREFUTnbA3zKv8FsKm2TqUV53AHwBx8DdlqrWxNTi\nhMLqF/lnza/dtjYZHEBRBWuqBSVOrI4UltmtPYUm61qpVnD5ImnPB4P1ydQFHGrWzRpwg6nFCYXV\nd3t6S7GckK2NR3yBQs3mJveZQv4N/PQpZBBrlwB4o3BtsOG4AIqhj9zRQ5mPhXI6KVdUfi9853Wu\nPwC+B27OdaAvgOKd1x5KoXCSucH1jLxROGSyS6FMksJPaPrrwp7t0PGwI7D9PxOwiySAojseGluU\nt8Pva8+HwCe+Ln7bC4E0X2wp9P2ltlFYmGi+NwobSqINEQrdc9nRQGofeFD8MwE1MhYCKHBYr7wU\nZfdx9Z7Dt22BA/nS7fs2ZmsQQNEO9i1LUW7HQs8b1DYKgYx1xXUUmrhuJYkPidOuB3T0wE/IgQMX\nzap6mR2leUUByHh4AEUTrWpwKKs+rLofss1UbqVkbInT2ON09AB4AOdcIDApe5md/trzG0352tZQ\nGDIZEUDRRKsaAGXV46uOQpup3ErJ3BKnrafS0QPgFVzlwkOjmpfZoW4oDJmMDqBoolUNirJqr6uO\nQpup3EppgCVOux682RcAj+DQRQ2uktlBdjcUhkxGB1A00aoGRFl1uaqj0GYqtxJYS+K068H1D4BH\ncOgCYqWiSjy5R+sOBSUxlHt3//JjYQJqk6ZbSXog2SL34Pp/hcKhFC+AFJ8hC6DIm2X4ZRs6JsXT\nFlW5xUoSBhKnS8+3KQTpXSD4ekGRHgxJVHyDQpDiyxtVyzcH5ZBGN7eYuG4lbi+qGdp8m0KO1Zfi\nEQUi+wKFJsUXHbl8C7yWQpup3EpiConT2vOXUKgye0IsUjwXDQXkqn6BQpPis3ZCsdPd6wSIQ/7j\nn/+TKvRpM5VbKXXaPmk9SOH//VOdvCu4gHoDvR1VyexedjcUhmyaAijaeL6loGyleE8hh4Qfsonr\nVkqOLXEaepDCwPaHgLMiaCWzO9ndUBgyQhpA4QlrawWlnhvIpOx95ZtHcchiSyURtnnSdCtl35o4\njT1fprCS2VF2BxSKrKxPRvynX4WWRoonx0MKbaZyK5Vs7pI4bT3k6csUmsxeS/GIoiATyr6wFzZS\n/ByFJq5bKb/oXBOnocfkekYfAF+273opHlAosm9R2EjxcxQKiLUfqyhcG2w4LoBi6CN32IY2y/ID\nLt/cwyF7tjZqaSkA/jMBu8ACKLrjobGH8pTfW1e+2XyjEEjzxZbCoiOX72S/Ueh5g1pLIXRCcaMQ\nyPDFjULPx4raRuEK0vyQjULPx4raRuEK0vyQjULPx4raRuEK0vyQdRSqwE0Pih/uOK2SvfPc96Q7\nJhw78L/AgZuT4ulhSErVkhR8WT3Lzfc9q1B4wtqaQ9l2a4u7LjSB20rJUu6ZyKuSqx673bmOwnkp\n/nylZ8Ll3YUCA3Lzfc8qFMrBoLCKQhO4rZT82zvPfQ+8hnwNhe+leAp/V33OcvNJpIeedSgGzGnz\nKgpN4LZS8mjvPHc9+BryNRSSc7zlqOhz4YHTBFKb5ea7npUo6mh1fRWFJnBbCRzLD9n14GvIv0Dh\nPd8XyRggN9/1rEQBa9YtrqIweVKB214fnjrsnee5x72G/AsUvq47eiN2SZEpaYz8imnsWYuiyxs0\nGoUfkuLdO8+TCO5fQ/55CuluP6cKasql5eZjz2oUwFa3qBS2UvzzcLd3nfOKq614srOblVIIe+d5\n6uEj+lWP9d+gUATSgyZu64vOz+nVzdKzGkWXN2gstLRSPL+MMQm0bN5QqAK3Sd3qtrzzPNlUryH/\nPIXTK+WgwzzcZXZ161mPQldrUCgUtlI8vx7ZZsStKTSB20oUwr3zPPXUryH/AoVJR37gJNKEhXPz\ntecPUAyY0+ZCYSvF85wQN92wFYUmcFuJfeI7z3NP/RryL1CYzrsPBZvWjnPztecPUChXg4JSmB5h\nFCuTnU76FIan0ARuK8lQeOe577HXkH+BwgTzoJc2lptf9axCMWBOmwuFPSment5yjzQW28kEbisl\nEdzeeW49EkteQy6lL1A43WkHOPMZOaGw2dWtZz0K5WpQKLT0pHhKZNRRbi80gdtKWYrXd55DD/3V\nv79eMg87+VtFoVfzFZUWaBZ3ycDPKGBSd+1ZjUKDDAqFwp4U/zAG2zPywN/75lUUvncbtAigeOdZ\nKQTDfCy80K/jUm4gub0QbOPFAPgeuHjA7ogAiu54aOyhTBL87no6nfiJW1k2CoE0X2wpLBK8pL/Q\nf5K0bBQWJprvlsLGRBo2Cvu8UOtG4ZCapR0xCjm1Df/HLI2S7X5lvB4c3o7eJXtNxHk7YJFBFMU7\npwtQakia4JgXuy/3znndT3/+eambh/U/Dtj1HEXRdQKNC1CWkP8PkQYkFXK2tDQAAAAASUVORK5C\nYII=\n", "prompt_number": 13, "text": [ "\u23a111.416\u23a4 = \u23a11 0.389\u23a4\u22c5\u23a1c\u23a4 + \u23a1e\u2081 \u23a4\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a24.514 \u23a5 \u23a21 0.2 \u23a5 \u23a3b\u23a6 \u23a2e\u2082 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a212.204\u23a5 \u23a21 0.241\u23a5 \u23a2e\u2083 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a214.835\u23a5 \u23a21 0.463\u23a5 \u23a2e\u2084 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a28.416 \u23a5 \u23a21 4.585\u23a5 \u23a2e\u2085 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a26.563 \u23a5 \u23a21 1.097\u23a5 \u23a2e\u2086 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a217.343\u23a5 \u23a21 1.642\u23a5 \u23a2e\u2087 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a213.02 \u23a5 \u23a21 4.972\u23a5 \u23a2e\u2088 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a215.19 \u23a5 \u23a21 7.957\u23a5 \u23a2e\u2089 \u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a211.902\u23a5 \u23a21 5.585\u23a5 \u23a2e\u2081\u2080\u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a222.721\u23a5 \u23a21 5.527\u23a5 \u23a2e\u2081\u2081\u23a5\n", "\u23a2 \u23a5 \u23a2 \u23a5 \u23a2 \u23a5\n", "\u23a322.324\u23a6 \u23a31 6.964\u23a6 \u23a3e\u2081\u2082\u23a6" ] } ], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": {}, "source": [ "In symbols:\n", "\n", "$$y = X B + \\epsilon$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We still haven't found our best fitting line. But before we go further, it might be obvious that we can easily add a new predictor here." ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Multiple regression" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's say we think that psychopathy increases with age. We add the student's age as another predictor:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "age = Matrix([22.5, 25.3, 24.6, 21.4,\n", " 20.7, 23.3, 23.8, 21.7,\n", " 21.3, 25.2, 24.6, 21.8])\n", "age" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left[\\begin{matrix}22.5\\\\25.3\\\\24.6\\\\21.4\\\\20.7\\\\23.3\\\\23.8\\\\21.7\\\\21.3\\\\25.2\\\\24.6\\\\21.8\\end{matrix}\\right]$$" ], "metadata": {}, "output_type": "pyout", "png": "iVBORw0KGgoAAAANSUhEUgAAADUAAAErCAMAAABJkzsdAAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\nTlMAMquZdlQQQO0wRCK7ie/dzWZ8bFsxDvgAAAfZSURBVGgF7Vtrl6QoDKXUcmartNRd/v9v3TyA\nPMAuunfOzpnT+qEbhUtCwGsSqHCLdA2h59q4cQi3OE5w3XtA4cCmQ0TUrQsgjQ6Deo7jslAXUuK2\n+zSH2/hIQIOacWhrhDoppXYvHM6UboJBDTM+fr1CkFJquK/D9Mwgi3pt+HyAYUoptbQGNrL2iMIQ\nJaX3KG6xERbLUgrDcUxjMbWRRahHXBkdpBTCAlo80U501ahtSVVBSunJQuOGmwo1FpCUcjdjTGb0\nqKkYS0oAGnZE0lRiwaHuIzx7Yo9SwmYvQo3ZThZ1Q1DA6ZQSPgn0PGyEhVuDem7ruo4DjFlKc8Sm\nT4Qdr7w6DGrndwdQUgobjfQ5rsNICw7gVhY+6LqMrC4ENrpQ2lTfwxrCnZY1YfnCul4KpxtrKO60\nrMlr/h6PZEiDUtxpWTNMEeXEJm8o7ixEwL0/Xkjj9KrhAyNLcadDEfZEQ+6XuNOyJlfszANeFlUy\nd1rWhIrbuhWQ1ZBQijsLa1JFWF+Zss24sFJzZ2FNRsE3KjEHon7EH+kxsBqbwbEmU2RYyP7Q9m/z\nhS3c6VgzROLPJU+z1bBwp2PN8CL235o8L9wprMlTe+BaupcPs5GluFNYk1n0MYJ7kD9fDcsXw3xQ\nMLI+aGerLpS2B1pjBn/sc9cT/LD/6h/2SvwV8yUsilL5Dcjy13Vam++yYlFsq7ggzDss3bHpbygW\nBdDELxXL2tGRW7KPZcalWBSY/dhzI6TN7KBQJwalWBRcraBQC/jC6jIofs4e6PHUqLjf1qk45o23\nkll0htkXWUADSAGv/CmqZbHlcPAaRRw1aj40K4pZ9I4zI6j0DTqy44uyahadyWFWKObHI/Nhm0Xv\nA17xNeQYg21oZImGhUXJmmqW72lc9Nx99YRFqTbiLKcP5AI6z20bKhaFWV5iHMDezKJhHcchL97G\nfCUlPvxXz9eHzVPlhdJWQmt8Mxbd/TdDWNWsDc2izyUHdsV4wqoGJSw6L+PgUYpVDcqw6ORQmlUN\nyrCoR2lWNSgeQYrjHcqwao3KcbxFWVatUdlSFmVZtUIVX9SgHKsiqmZRHKBGeVZts6hHeVa1GmoW\n1bLYuMKqBmVYNKEkzAjEqtSBQSkWHYdX3JBEM4sWVkWYQVE/PX8ulLbS162h16Hu8bzs1uF5Q1Pz\ndQ3FBzAdfnBjZYkvOo/gbZWvd7htwwgx/Urei1uHikV3cLTmrTiFB2eqYgqjjSxh0YP8wTV7rGFk\nfI69DUpYdKF69mhwfKxZ21sWFmVVHjlmZcM8slfVer+ARedI4eijBJIE026wjziQRZ+cG72VFClp\niW82X2Zc9AhZtIWagQPzhSjzNScWnVmI0XAqBoVeXUyUInoe111bIzM5ynMa5mwoNznUJ2zOM1yj\nCotONMslLoGWkPrH9nwZWYpFN1xRmM/ILHo/RSkWnVfYJ6DVm3zRZ3FfQZyRlRV4+/9CaRN93Rrf\niUVtHrQzL0pZLMkyduZFXR60My/q8qDB3ZtZFhalhaA0dPcGxYsm+aIQL6t3HqvkvkYlX9TkQQFi\n7muUMJjkQVkJuUdUzaLcSvKg/v6ERWmrSPKgPi/qNMws6vKgnXlRlwftzItKHpRZVO5piEZDxaKS\nB73yomktnP4zNjxt5Su+jvrOLDpP03qSCRRf1LMob8QOzUyg+KKeRe/krz2zK2XmS3xRx5oQsuFE\nz809qXMWfUTcAlqzh2hk8cJpsugQ41FADd/mhEXHGCVNUstqs+gx3vaIW2d0IaqDRQ/01I+Ys6md\nLMrNb9k3dRqesGj2lCmAJqnaMy++aMWiHHEsKXAxssQXFdZkFn2Q/Y4cPBjUByx6g0jqZPWyWTv+\nGlkd7bnJhdKmQmt8JxaFscsLD69FX0QfbF60L6L3edG+iB4UNHm2vojeo2iK30X02MjIItS7iL6J\nEuptcG9mUS8r8xT2iOuwzaIe1RHRQ38elXXIslSOKLNojeqJ6GtUV0RfNPxURF/lRa+IHie0ef3f\nPP+NWbTzpJNjUeebmvkSX9SzqPNNDUp8UVgRhgGcb2pQxhc1KFpXEuEbFK+59u4S1HVE9FbD7oi+\n1rAnoq9REuGf+KLehi7Cd9Y4Y1EX4VtU8UW9LOebGpT4ooL6VETvWVQifOiw8f3Cp+8uo+G7xqX+\nQhVTQAGt8SexqOwTVSc/y8l6HJ+ZZYnNgzv5qWoqlPBhcCc/VU2FUnyY95GwDVyqBm+NhooPHUrV\nVCh8kE7Nt05+ynl6I4tAHJtXJz+hLkXtUKpRyqNzJz+lpkKV3Xbo0p78VDWI0uvw5OQndKHP07td\nzsyH/uSnO09vNSx86E9+uvP0BiV86E9+Sg3o6myoYnN38lPVVCjqp+eP0bAHQG0ulDYVWsPERLry\ntOy/5qcNTcXvs7ywKCrEb0BW7Tr5mS0R/syTnz//KgNgrvRnlNLvicrJz39+6lxlPj/vzyiF6+Rn\nsWtd+H3vcq1L+8mlobbLr7CGYVGd4YQvSMc+EWzK25Of5IC82yfysXnfPhEM3USj3ftEBtW9T2RQ\n8LPQrn0iqyH6il37RFZWvU+kso7qV0gG5faJ/NoQj9Og+vaJnOVz/vXNPpFDhfR7zDf7RIL61D6R\nj82vfSIw5Mnl18ZJM/eYUXQm0YV3rmG+Lb9SBybBK2985vr2f/qV+jSFfwG5JaSsUYBipgAAAABJ\nRU5ErkJggg==\n", "prompt_number": 14, "text": [ "\u23a122.5\u23a4\n", "\u23a2 \u23a5\n", "\u23a225.3\u23a5\n", "\u23a2 \u23a5\n", "\u23a224.6\u23a5\n", "\u23a2 \u23a5\n", "\u23a221.4\u23a5\n", "\u23a2 \u23a5\n", "\u23a220.7\u23a5\n", "\u23a2 \u23a5\n", "\u23a223.3\u23a5\n", "\u23a2 \u23a5\n", "\u23a223.8\u23a5\n", "\u23a2 \u23a5\n", "\u23a221.7\u23a5\n", "\u23a2 \u23a5\n", "\u23a221.3\u23a5\n", "\u23a2 \u23a5\n", "\u23a225.2\u23a5\n", "\u23a2 \u23a5\n", "\u23a224.6\u23a5\n", "\u23a2 \u23a5\n", "\u23a321.8\u23a6" ] } ], "prompt_number": 14 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now rename the `clammy` predictor from $x$ to $x_1$. Call the `age` predictor $x_2$. Call the slope for `clammy` $b_1$ (slope for $x_1$). Call the slope for age $b_2$ (slope for $x_2$). Our model is:\n", "\n", "$$y_i = c + b_1 x_{1, i} + b_2 x_{2, i} + e_i$$\n", "\n", "In this model $X$ has three columns (ones, $x_1$, and $x_2$), and the $B$ vector has three values $c, b_1, b_2$. This gives the same matrix formulation, with our new $X$ and $B$: $y = X B + \\epsilon$.\n", "\n", "This is a *linear* model because our model says that the data $y_i$ comes from the *sum* of some components ($c, b_1 x_{1, i}, b_2 x_{2, i}, e_i$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can keep doing this by adding more and more regressors. In general, a linear model with $p$ predictors looks like this:\n", "\n", "$$y_i = b_1 x_{1, i} + b_2 x_{2, i} + ... b_p x_{p, i} + e_i$$\n", "\n", "In the case of the models above, the first predictor $x_1$ would be a column of ones, to express the intercept in the model.\n", "\n", "Any model of the form above can still be phrased in the matrix form:\n", "\n", "$$y = X B + \\epsilon$$" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Solving the model with matrix algebra" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The reason to formulate our problem this way is so we can use some basic matrix algebra to find the \"best\" line." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's assume that we want to find the line (intercept and slope) that gives the smallest \"distance\" between the estimated values (from the line), and the actual values (the data).\n", "\n", "We'll define 'distance' as the squared difference of the predicted value from the actual value. These are the squared error terms $e_1^2, e_2^2 ... e_{n}^2$, where $n$ is the number of observations - 12 in our case." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Revising: our model is:\n", "\n", "$$\n", "y = \\bf{X} B + \\epsilon\n", "$$\n", "\n", "Where $y$ is the data vector $y_1, y_2, ... y_n$, $\\bf{X}$ is the design matrix of shape $n, p$, $B$ is the parameter vector, $\\beta_1, \\beta_2 ... \\beta_p$, and $\\epsilon$ is the error vector giving errors for each observation $\\epsilon_1, \\epsilon_2 ... \\epsilon_n$.\n", "\n", "Each column of $\\bf{X}$ is a regressor vector, so $\\bf{X}$ can be thought of as the column concatenation of $p$ vectors $x_1, x_2 ... x_p$, where $x_1$ is the first regressor *vector*, and so on." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In our case, we want the vector $B$ such that the errors $\\epsilon = y - \\bf{X} B$ have the smallest sum of squares $\\sum_{i=1}^n{e_i^2}$. $\\sum_{i=1}^n{e_i^2}$ is called the *residual sum of squares*.\n", "\n", "It might or might not be obvious that this also means that the design matrix $\\bf{X}$ should be orthogonal to the errors - meaning that $\\bf{X}^T \\epsilon$ should be a vector of zeros.\n", "\n", "If that is the case then we can multiply $y = {\\bf X} \\boldsymbol{\\beta} + \\epsilon$ through by $\\bf{X}^T$:\n", "\n", "$$\n", "\\bf{X}^T y = \\bf{X}^T X B + \\bf{X}^T \\epsilon\n", "$$\n", "\n", "The last term now disappears because it is zero and:\n", "\n", "$$\n", "\\bf{X}^T y = \\bf{X}^T \\bf{X} B\n", "$$\n", "\n", "If $\\bf{X}^T \\bf{X}$ is invertible then there is a unique solution:\n", "\n", "$$\n", "B = (\\bf{X}^T \\bf{X})^{-1} \\bf{X} y\n", "$$\n", "\n", "It turns out that, if $\\bf{X}^T \\bf{X}$ is not invertible, then are an infinite number of solutions, and we have to choose one solution, taking into account that the parameters $B$ will depend on which solution we chose. The *pseudoinverse* operator gives us one particular solution. If $\\bf{A}^-$ is the pseudoinverse of matrix $\\bf{A}$ then the general solution for $B$, even when $\\bf{X}^T \\bf{X}$ is not invertible, is:\n", "\n", "$$\n", "B = \\bf{X}^-y\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So what line do we estimate for `psychopathy` and `clammy`?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "X = np.column_stack((np.ones(12), clammy))\n", "X" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 15, "text": [ "array([[ 1. , 0.389],\n", " [ 1. , 0.2 ],\n", " [ 1. , 0.241],\n", " [ 1. , 0.463],\n", " [ 1. , 4.585],\n", " [ 1. , 1.097],\n", " [ 1. , 1.642],\n", " [ 1. , 4.972],\n", " [ 1. , 7.957],\n", " [ 1. , 5.585],\n", " [ 1. , 5.527],\n", " [ 1. , 6.964]])" ] } ], "prompt_number": 15 }, { "cell_type": "code", "collapsed": false, "input": [ "B = npl.pinv(X).dot(psychopathy)\n", "B" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 16, "text": [ "array([ 10.07128585, 0.99925723])" ] } ], "prompt_number": 16 }, { "cell_type": "code", "collapsed": false, "input": [ "plt.plot(clammy, psychopathy, '+')\n", "def my_best_line(x):\n", " # Best prediction for psychopathy given clamminess\n", " return B[0] + B[1] * x\n", "x_vals = [0, max(clammy)]\n", "y_vals = [my_best_line(0), my_best_line(max(clammy))]\n", "plt.plot(x_vals, y_vals)\n", "plt.xlabel('Clamminess of handshake')\n", "plt.ylabel('Psychopathy score')" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 17, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEPCAYAAACukxSbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYVPe5B/DvCOIGiiA7MRBwYQdnGK5xA2PUVE3coqIt\nLqipJuk1sbG2aRs1uUZvbptq0iY2JmatFpPGxtj4NIvjlkZRNsVdmYhsioCIgMDwu38QTjwsDsjM\nnIHz/TzPPI8zDGdeRL/8eOec96cRQggQEZEqdFO6ACIish2GPhGRijD0iYhUhKFPRKQiDH0iIhVh\n6BMRqYjVQj83NxcJCQkICwtDeHg4Nm/eDABYs2YN/P39ERMTg5iYGOzdu9daJRARURMaa52nX1hY\niMLCQkRHR6OiogJarRa7du1CSkoKXFxc8Oyzz1rjZYmI6C4crXVgb29veHt7AwCcnZ0REhKCvLw8\nAACvByMiUoZNevpGoxHp6en4r//6LwDAa6+9hqioKCQnJ6OsrMwWJRAREWwQ+hUVFZg5cyY2bdoE\nZ2dnLFu2DDk5OcjIyICPjw9Wrlxp7RKIiKiRsKKamhoxfvx48eqrr7b48ZycHBEeHt7s8aCgIAGA\nN9544423dtyCgoLM5rLVVvpCCCQnJyM0NBQrVqyQHi8oKJD+/OmnnyIiIqLZ5168eBFCCLu/vfDC\nC4rXwDpZJ+tkjY23ixcvms1mq72Re/jwYXz44YeIjIxETEwMAGD9+vXYvn07MjIyoNFoEBgYiC1b\ntlirBCIiasJqoT9y5EjU19c3e/yRRx6x1ksSEZEZvCK3A+Lj45UuoU1Yp2WxTsvqDHV2hhrbymoX\nZ3WERqOBHZZFRGTX2pKdXOkTEakIQ5+ISEUY+kREKsLQJyJSEYY+EZGKMPSJiFSEoU9EpCIMfSIi\nFWHoExGpCEOfiEhFGPpERCrC0CeidjEYlK6AOoKhT0TtwtDv3Bj6REQqYrVNVIio6zAYflzhr137\n4+Px8Q23rsxg6FpfI0OfiMxqGu5r1ihUiAK6WuizvUNEpCJc6RNRu3SlVW9runI7i9slEhHdxZo1\nnaedxe0SiYhIhqFPRHQXnb2d0xTbO0REXQTbO0REJMPQJyJSEYY+EZGKMPSJiFSEoU9EpCIMfSIi\nFWHoExGpCEOfiEhFGPpERCrC0CciUhGGPhGRijD0iYhUxGqhn5ubi4SEBISFhSE8PBybN28GAJSU\nlODhhx/G4MGDMX78eJSVlVmrBCIiasJqUzYLCwtRWFiI6OhoVFRUQKvVYteuXdi2bRsGDBiAVatW\nYePGjSgtLcWGDRvkRXHKJhFRuyk6ZdPb2xvR0dEAAGdnZ4SEhCAvLw+fffYZ5s+fDwCYP38+du3a\nZa0SiIioCZvM0zcajRgzZgxOnjyJgQMHorS0FAAghICbm5t0XyqKK30ionZrS3ZafWP0iooKzJgx\nA5s2bYKLi4vsYxqNBhqNpsXPW3PHppTx8fGI72rb1xARdZDBYIChcQf3NrLqSr+2thaTJ0/GI488\nghUrVgAAhg4dCoPBAG9vbxQUFCAhIQFnzpyRF8WVPhFRuyna0xdCIDk5GaGhoVLgA8Cjjz6K9957\nDwDw3nvvYerUqdYqgYiImrDaSv/QoUMYPXo0IiMjpRbOyy+/DL1ej1mzZuHy5csICAhASkoKXF1d\n5UVxpU9E1G5tyU5ujE5E1EVwY3QiIpJh6BMRqQhDn4hIRRj6REQqwtAnIlIRhj4RkYow9ImIVISh\nT0SkIgx9IiIVYegTEakIQ5+ISEUY+kREKmL1TVSIiMjyakw1OFF0AkfyjuD89fN4deKrbfo8hj4R\nkZ0TQuBS6SUczTuKI3lHcDTvKDKLMhHUPwh6Pz3i/OJQL+rbdCyOViYisjPXK68jNT8VR64cwdH8\noziadxQ9HXtKAa/300Pro4VLj+Zb0HKePhGRHbtddxsZhRnSCv5I3hEUVRRB56tDnF8c4vwbQt7X\nxdfssRj6XYjBAHBveKLOrV7U4/z187I2Tfa1bAx2H9wQ8D+s4ocOGAqHbg7tPn5bstNsT//WrVv4\n4x//iMuXL+Ott97C+fPncfbsWUyePLndBdG9Y+gTdT5Xb11tCPgf2jSpeano17Of1KaZEz4Hw3yG\noXf33jaryWzoL1y4EFqtFt9++y0AwNfXFzNnzmToExHdobK2EukF6TiSd0RaxZdVlyHWNxZxfnF4\nKvYp6Kfq4eXspWidZkP/4sWLSElJwY4dOwAAffr0sXpR1MBgaLgBwNq1Pz4eH89VP5GS6kU9Tl87\nLWvTnL1+FmEeYdD76TF50GSsi1+HQe6D0E1jX5dDmQ39Hj16oKqqSrp/8eJF9OjRw6pFUYOm4b5m\njUKFEKlc/s18WZvmWP4xePbxlNo0C6IXINo7Gj0deypdqllmQ3/NmjWYOHEirly5grlz5+Lw4cN4\n9913bVAaEZHtVdRU4Hj+cVmbpqq2Sgr45x58DrG+sXDv7a50qffkrmfv1NfXY+fOnXjooYfw3Xff\nAQDi4uLg4eFh3aJ49k4zfCOXyPLq6uuQfTVb1qa5WHoRUV5RsnPiH+j/ADQajdLlmmWRUza1Wi2O\nHz9u0cLMYegTkaUJIZBbnitr06QVpMG/r78s4CO9IuHk4KR0uffEIqG/evVqDBgwALNnz5a9ievm\n5maZKlsqiqFPRB10o/oGjuUfk7Vp6kW97Hz4WL9YuPZ0VbpUi7FI6AcEBDT7tUaj0eDSpUsdr7C1\nohj6RNQOtaZanLh6Akeu/Bjwl29cRoxPjCzkB/Yb2CnaNPeKV+QSUZcjhEBOWY6sTZNZmInA/oHQ\n++qlsQXhnuFw7KaumZIWCf2amhq88cYbOHDgADQaDcaMGYOf//zn6N69u0WLlRXF0CeiH5RUlSA1\nL1VawR/NOwonBydZH17nq2s2fEyNLBL6ycnJqKurw/z58yGEwAcffABHR0ds3brVosXKimLoE6nS\n7brbyCzKlLVpCisKofXVyto0fn39lC7VLlkk9CMjI5GVlWX2MUti6BN1fUIInC85L2vTnLx6EoPd\nB8vaNCEDQu5p+JgaWWTgmqOjIy5cuIDg4GAADVfkOjqqq09GRB137dY16Xz4I3lHkJqXir49+kpt\nmllhszDMZxj6OHHUizWZXel//fXXWLhwIQIDAwEARqMR27Ztw9ixY61XFFf6RJ1aVW0V0grSZBc9\nlVSVINYvVmrTxPrFwtvZW+lSuxSLnb1TXV2Ns2fPQqPRYPDgwejZ07rzJRj6RJ1HvajHmeIzsjbN\nmeIzCPUIlbVpBrsPtrvhY12NRUL/9ddfx7x589C/f38AQGlpKbZv347ly5dbrtKmRTH0iexWwc0C\nWZvmWP4xePT2kNo0cf5xnWb4WFdjkdCPiopCZmam7LHo6GhkZGR0vMLWimLoE9mFxuFjd7ZpbtXe\n+jHgf2jTDOg9QOlSCRZ6I7e+vh719fXo1q3h1zKTyYTa2to2FbBo0SLs2bMHnp6eOHHiBICGqZ1b\nt26Vhra9/PLLmDhxYpuOR0TWY6o3IftatqxNc6HkAiI8IxDnF4fpIdOxYdwGBPUP6tJXtXZG7RnI\naDb0J0yYgDlz5uCJJ56AEAJbtmxpc0gvXLgQTz/9NJKSkqTHNBoNnn32WTz77LNtq5CILE4IgSvl\nV2RtmrSCNPi5+Emr+KXapYjyjuq0w8fUxKKhv3HjRvz1r3/FG2+8AQB4+OGHsXjx4jYdfNSoUTAa\njc0eZ+uGyLbKb5c3DB+746Knuvo6xPk3tGh+O+q30Pnq0L9Xf6VLJSszG/oODg5YtmwZli1bhpKS\nEuTm5sLBoWMXSrz22mt4//33odPp8Ic//AGurtadcsdZ9KQmjcPH7uzDf1/2PWJ8YqD31WNuxFz8\naeKfcH+/+9mm6cRa207VHLNv5I4ZMwa7d+9GXV0dtFotPDw8MGLECLz66qttegGj0YgpU6ZIPf2r\nV69K/fzf/e53KCgowNtvvy0vysJv5K5Zw60GqWsSQsBYZpQFfEZhBgJcA2SzacI9w9HdwXrzskhZ\njRlnkTdyb9y4gb59+2Lr1q1ISkrC2rVrERERcc/FeXp6Sn9evHgxpkyZ0uLz1tyR0vHx8YjnUp0I\npVWlSM1PlbVpHLs5Sm2aFxNehNZXi749+ipdKtmAwWD44db2ha3Z0DeZTCgoKEBKSgpeeuklAOjQ\nr4QFBQXw8fEBAHz66aet/gBZ08GleWu/+jTdbJzIXjUOH7tzFV9wswBaXy30vnosilmENye/Cf++\n/kqXSgppXBA3trDXtqHPYzb0f//732PChAkYMWIE9Ho9Ll68iEGDBrWpoMTEROzfvx/FxcW47777\nsHbtWhgMBmRkZECj0SAwMBBbtmxp07Haq2m4s71D9kwIgQslF2Rn05y8ehKD3AZB76dH/P3xWPXg\nKoR6hHL4GDXTnoWsKjZRYU+f7E3j8LHGkE/NT4Wzk7PsoicOH6P2skhPvytgO4eUVFVbhfTCdFmb\npriyGLG+DcPHlscuR6xvLHxcfJQulVRAFSt9IlupF/U4W3xW1qY5U3wGIQNCZGfTDBkwhMPHyOIs\nMnvHZDJ1+Lz89mLoU2dRWFHYMLLgh5A/ln8M7r3dZQEf4x2DXt17KV0qqYBFQv+BBx7AjBkzsHDh\nQoSGhlq0wFaLYuiTHbpVcwvHC+TDx27evimbLhnrGwuPPh5Kl0oqZZHQLy8vx44dO/Duu+/CZDJh\n0aJFSExMRN++1jsPmKFPSjPVm3Dq2ilZm6Zx+Nidq/hgt2Be1Up2w2KbqDQyGAyYN28eSktL8fjj\nj+N3v/udtI2iJTH0ydak4WNXfhw+5uPiIzubJtIrEj0ceyhdKlGrLBL6dXV12LNnD7Zt2waj0Yik\npCTMnTsXhw4dwm9+8xucO3fOokUDDH2yrsbhY3e2aWpNtQ07PP2w05POVwe3Xm5Kl0rULhbr6cfH\nx2Px4sV48MEHZR97+umn8dprr3W80qZFMfTJQurq63CiSD58zFhmRLR3tKxNE+AawDYNdXoWCf2b\nN2/CxcXFooWZw9CneyGEwPc3vpe1aTIKM3C/6/2ygI/wjODwMeqSLBL6V69exVtvvQWj0Yi6ujrp\nwO+8847lKm1aFEOf2qCsukx2VevRvKPopukm9eD1fnrofHXo17Of0qUS2YRFQn/48OEYPXo0tFqt\ntGWiRqPBjBkzLFdp06IY+tREjakGmYXy4WN5N/Og9dHKVvH+ff3ZpiHVskjoW3sT9JYw9NVNCIGL\npRdlFz2duHoCwW7B0hutej89Qj1C4dhNFZNE2o0bB6mTRWbvTJ48GXv27MGkSZMsVhjRnYori5u1\nafp07yOt4GeEzsAwn2FwdnJWutROg6FPrWl1pe/s7Cz9mnzr1i04OTmhe/eGN780Gg3Ky8utVxRX\n+l1WdV010gvSZRc9NQ4fu7NNw+FjHcPJsurUoZV+RUWFxQsidakX9Th3/ZysTXO6+DSGDhgKva8e\n44PG47ejf4uhA4Zy+JgFcOMgaguzPf2HHnoIX3/9tdnHLFoUV/qdUlFFkdSeOZJ3BKl5qXDr5Sab\nTcPhY7bBlb46dWilX1VVhcrKSly7dg0lJSXS4+Xl5cjLy7NcldQpVdZW4nj+cVmbpnH4mN5PjxVx\nKxDrFwvPPp7mD0ZENtNq6G/ZsgWbNm1Cfn4+tFqt9LiLiwueeuopmxRH9sFUb8Lp4tOyi57Ol5xH\nuGc49L56PDrkUbw09iUMchvE0yXtBNs51Bqz7Z3NmzfjF7/4ha3qAdC+9g7PUrC8vPI8WZvmeP5x\neDt7y9o0UV5RHD5GZGcsNmXz5MmTOHXqFKqrq6XHkpKSOl5ha0W1I/TZu+yYm7dvNhs+dtt0WzZd\nMtYvlsPHiDoBi5ynv2bNGuzfvx/Z2dmYNGkSvvjiC4wcOdKqoU/WUVdfh5NXT8raNDllOQ3Dx3z1\nmBU2C/83/v8Q6BrINg1RF2U29D/++GNkZmZi2LBh2LZtG4qKijBv3jxb1NYqnppmnhACl29clrVp\n0gvSMbDfQGkVvzx2OSK9Ijl8jO4JW6udk9nQ79WrFxwcHODo6IgbN27A09MTubm5tqitVU3Dne2d\nhuFjqXmpsjaNRqORLnZaM2YNh4+RRTH0OyezoR8bG4vS0lIsWbIEOp0Offr0aTZXn2yrxlSDrKKs\nhoue8htaNXk38zDMZxj0vnokRSXh9Z+8jvv63sc2DRHJtGu7RKPRiPLyckRGRlqzJp69cwchBC6V\nXpK1abKKshDUP0h2Ng2Hj5EtNG2tvvBCw5/ZWrUPFjl7RwiBf/zjHzh06BA0Gg1GjRqFadOmWbTQ\nZkVZ4IrczvrD4Hrl9WbDx3p17yW1aeL84qD11XL4GCmOZ87ZH4ucvbN8+XJcvHgRiYmJEEJgy5Yt\n+PLLL/GXv/zFYoVaQ2cI/eq6amQUZsjaNNcqr0Hnq4PeV4+l2qXY+uhW+Lr4Kl0qEXURZkN/3759\nOHXqlLSByoIFCxAaGmr1wrqaelGP89fPy9o0p66dwhD3IdD76TEucByeH/U8hrgPgUM3B6XLJTLL\n3hdV1DKzoR8cHIzLly8jICAAAHD58mUEBwdbu657Yk+ncl69dVU2XTI1PxWuPV2lFk1ieCJifGLQ\nu3tv2xZGZCEM/c7JbE9/9OjRSE1NhV6vh0ajwdGjRxEbG4u+fftCo9Hgs88+s3xRFujp27LfWFlb\nibSCtB+Hj105ghu3bzQMH/thp6dY31h4OXvZpiAiUiWL9PTXrVsnHQyA7IBqPB3QVG/CmeIzsjbN\nuevnEOYRBr2fHpMHTca6+HUY5D6IM+KJyO606ZTNwsJCpKamQqPRQK/Xw9PTuuNy7ensnfyb+bI2\nzfGC4/Ds4ymbTRPlHYWejj07/mJERB1gkVM2U1JS8Nxzz2HMmDEAgAMHDuCVV17B448/brlKmxal\n0CYqFTUVsuFjR64cQXVddcNG3He0adx7u9u8NiIicywS+pGRkfjqq6+k1f21a9fw0EMPISsry3KV\nNi3KBqFfV1+H7KvZsjbNpdJLiPKKku3V+kD/B1TZxiKizsciPX0hBDw8PKT77u7unW4rQyEEcstz\nZW2a9MJ0+Pf1lwJ+mW4ZIrwi4OTgpHS5RERWYzb0J06ciAkTJmDu3LkQQuDvf/87HnnkEVvUds9u\nVN9Aan6qrE0DQGrT/H7M76Hz1cG1p6vClRIR2ZZVxzAsWrQIe/bsgaenJ06cOAEAKCkpwezZs/H9\n998jICAAKSkpcHWVh2972js1phqcKDoha9Pk3shtGD52x2waDh8joq7OYjtn3auDBw/C2dkZSUlJ\nUuivWrUKAwYMwKpVq7Bx40aUlpZiw4YNbSpcCIGcshxZmyazKBMP9H9AeqM1zi8OYZ5hHD5GRKpj\nkdD/5JNPsHr1ahQVFUkH02g0KC8vb1MRRqMRU6ZMkUJ/6NCh2L9/P7y8vFBYWIj4+HicOXOmxcJL\nqkqkXZ6O5jcMIevh0EO2gtf6aOHSw6VNtRARdWUWCf2goCB8/vnnCAkJuacimoZ+//79UVpaCqBh\n5e7m5ibdv7Pw4M3BKKooahg+dsfZNH59/e6pDiKirs4iZ+94e3vfc+Cbo9FoWu2zfzr7U4QMCOHw\nMSIiC2o19D/55BMAgE6nw+zZszF16lQ4OTWczqjRaDB9+vR7esHGto63tzcKCgpavbr34798LP05\nPj4e8ZzuREQkYzAYYGicMtlGrbZ3FixYIJu303RFvm3btja9QNP2zqpVq+Du7o5f/epX2LBhA8rK\nytr8Ri4REbVO8bN3EhMTsX//fhQXF8PLywvr1q3DY489hlmzZknjmjt6yiYRETWwSOjPnz8fmzZt\nkoK5tLQUK1euxDvvvGO5SpsWxdAnImq3tmSn2dm/mZmZspV4//79kZaW1vHqiIjI5syGvhACJSUl\n0v2SkhKYTCarFkVERNZh9pTNlStXYvjw4Zg1axaEENi5cyeef/55W9RGREQW1qY3crOzs/HNN99A\no9Fg7NixVt8YnT19IqL2s8gbuc8++yySk5MRFhZm0eLuhqFPRNR+FnkjNyQkBEuXLoVer8ebb76J\nGzduWKxAIiKyrTafp3/mzBm8++67+Nvf/oaRI0diyZIlSEhIsE5RXOkTEbWbRVb6AGAymXDmzBmc\nPn0aHh4eiIqKwh//+EfMnj3bIoUS0Y/aeVU9UbuYXek/88wz2L17N8aOHYvFixdDr9dLHxsyZAjO\nnj1r+aK40icVW7Om4UbUXhaZshkZGYmXXnoJffr0afaxI0eO3Ht1RERkc2ZDf9CgQdJPjg8++ABp\naWlYsWIF7r///mYzc4jo3hgMP7Z11q798fH4+IYbkaWYbe9EREQgMzMTJ06cwIIFC7B48WKkpKRg\n//791iuK7R1SMbZ36F5Z5I1cR0dHdOvWDbt27cKTTz6JJ598Ejdv3rRYkUREZDtm2zsuLi5Yv349\nPvzwQxw8eBAmkwm1tbW2qI1IldjOIWsyu9JPSUlBz5498c4778Db2xt5eXl47rnnbFEbkSox9Mma\nWu3pV1VV4c0338SFCxcQGRmJ5ORkODqa/cXAMkWxp09E1G4dmr0za9YsODk5YeTIkdi7dy/uv/9+\nbNq0ySqFNiuKoU9E1G4dCv2IiAhpX9u6ujrExsYiPT3d8lW2VBRDn4io3Tp09s6drRxbtXWIiMi6\nWl3pOzg4oHfv3tL9qqoq9OrVq+GTNBqUl5dbryiu9ImI2q1DYxi4JSIRUdfTpimbRETUNTD0iYhU\nhKFvQZyDTkT2jqFvQQx9IrJ3DH0iIhXhCfgdxDnoRNSZMPQ7qGm4cw46EdkztneIiFSEoW9BbOcQ\nkb0zu12iEjiGgYio/SyyXSIREXUdDH0iIhVh6BMRqQhDn4hIRRj6REQqotjFWQEBAejbty8cHBzQ\nvXt3HD16VKlSiIhUQ7HQ12g0MBgMcHNzU6oEIiLVUbS9w3PxiYhsS7HQ12g0GDduHHQ6Hd566y2l\nyiAiUhXF2juHDx+Gj48Prl27hocffhhDhw7FqFGjpI+vuWNyWXx8POI544CISMZgMMDQzo087GIM\nw9q1a+Hs7IyVK1cCaNulxAYDZ90QEd3JbscwVFZW4ubNmwCAW7du4d///jciIiLadQzuUkVE1H6K\ntHeKioowbdo0AEBdXR3mzZuH8ePHK1EKEZGqKBL6gYGByMjIaPfncZcqIqKO6VQ7Z3GXKiKijuEY\nBiIiFem0oc92DhFR+9nFKZtNcecsIqL2s9tTNomISBkMfSIiFWHoExGpCEOfiEhFGPpERCrC0Cci\nUhGGPhGRijD0iYhUhKFPRKQiDH0iIhVh6BMRqQhDn4hIRRj6REQqwtAnIlIRhj4RkYow9ImIVISh\nT0SkIgx9IiIVYegTEakIQ5+ISEUY+kREKsLQJyJSEYY+EZGKMPSJiFSEoU9EpCIMfSIiFWHoExGp\nCEOfiEhFGPpERCrC0CciUhGGPhGRiigS+nv37sXQoUMxaNAgbNy4UYkSiIhUyeahbzKZ8NRTT2Hv\n3r04deoUtm/fjtOnT9u6DIswGAxKl9AmrNOyWKdldYY6O0ONbWXz0D969CiCg4MREBCA7t27Y86c\nOfjnP/9p6zIsorP8Q2CdlsU6Lasz1NkZamwrm4d+Xl4e7rvvPum+v78/8vLybF0GEZEq2Tz0NRqN\nrV+SiIgaCRv7z3/+IyZMmCDdX79+vdiwYYPsOUFBQQIAb7zxxhtv7bgFBQWZzWCNEELAhurq6jBk\nyBB8/fXX8PX1hV6vx/bt2xESEmLLMoiIVMnR5i/o6IjXX38dEyZMgMlkQnJyMgOfiMhGbL7SJyIi\n5djdFbmd4cKtRYsWwcvLCxEREUqXcle5ublISEhAWFgYwsPDsXnzZqVLalF1dTXi4uIQHR2N0NBQ\n/PrXv1a6pFaZTCbExMRgypQpSpfSqoCAAERGRiImJgZ6vV7pclpVVlaGmTNnIiQkBKGhofjuu++U\nLqmZs2fPIiYmRrr169fPbv8fvfzyywgLC0NERATmzp2L27dvt/xEi75L20F1dXUiKChI5OTkiJqa\nGhEVFSVOnTqldFnNHDhwQKSlpYnw8HClS7mrgoICkZ6eLoQQ4ubNm2Lw4MF2+fcphBC3bt0SQghR\nW1sr4uLixMGDBxWuqGV/+MMfxNy5c8WUKVOULqVVAQEB4vr160qXYVZSUpJ4++23hRAN3/eysjKF\nK7o7k8kkvL29xeXLl5UupZmcnBwRGBgoqqurhRBCzJo1S7z77rstPteuVvqd5cKtUaNGoX///kqX\nYZa3tzeio6MBAM7OzggJCUF+fr7CVbWsd+/eAICamhqYTCa4ubkpXFFzV65cwb/+9S8sXrwYws67\novZe340bN3Dw4EEsWrQIQMN7ff369VO4qrv76quvEBQUJLvOyF707dsX3bt3R2VlJerq6lBZWQk/\nP78Wn2tXoc8Lt6zHaDQiPT0dcXFxSpfSovr6ekRHR8PLywsJCQkIDQ1VuqRmnnnmGbzyyivo1s2u\n/ts0o9FoMG7cOOh0Orz11ltKl9OinJwceHh4YOHChRg2bBiWLFmCyspKpcu6qx07dmDu3LlKl9Ei\nNzc3rFy5EgMHDoSvry9cXV0xbty4Fp9rV/96eeGWdVRUVGDmzJnYtGkTnJ2dlS6nRd26dUNGRgau\nXLmCAwcO2N1l759//jk8PT0RExNj96vow4cPIz09HV988QX+/Oc/4+DBg0qX1ExdXR3S0tKwfPly\npKWloU+fPtiwYYPSZbWqpqYGu3fvxuOPP650KS26ePEi/vSnP8FoNCI/Px8VFRX46KOPWnyuXYW+\nn58fcnNzpfu5ubnw9/dXsKLOr7a2FjNmzMBPf/pTTJ06VelyzOrXrx8mTZqEY8eOKV2KzLfffovP\nPvsMgYGBSExMxDfffIOkpCSly2qRj48PAMDDwwPTpk3D0aNHFa6oOX9/f/j7+yM2NhYAMHPmTKSl\npSlcVeu++OILaLVaeHh4KF1Ki44dO4YHH3wQ7u7ucHR0xPTp0/Htt9+2+Fy7Cn2dTofz58/DaDSi\npqYGf/9L6PIRAAAIqElEQVT73/Hoo48qXVanJYRAcnIyQkNDsWLFCqXLaVVxcTHKysoAAFVVVfjy\nyy8RExOjcFVy69evR25uLnJycrBjxw6MHTsW77//vtJlNVNZWYmbN28CAG7duoV///vfdnmWmbe3\nN+677z6cO3cOQEO/PCwsTOGqWrd9+3YkJiYqXUarhg4diu+++w5VVVUQQuCrr75qvUVqozeX2+xf\n//qXGDx4sAgKChLr169XupwWzZkzR/j4+AgnJyfh7+8v3nnnHaVLatHBgweFRqMRUVFRIjo6WkRH\nR4svvvhC6bKaycrKEjExMSIqKkpERESI//3f/1W6pLsyGAx2e/bOpUuXRFRUlIiKihJhYWF2+39I\nCCEyMjKETqcTkZGRYtq0aXZ79k5FRYVwd3cX5eXlSpdyVxs3bhShoaEiPDxcJCUliZqamhafx4uz\niIhUxK7aO0REZF0MfSIiFWHoExGpCEOfiEhFGPpERCrC0CciUhGGPrVLYWEh5syZg+DgYOh0Okya\nNEm6oM7WFwFNmjQJ5eXlNn3N1uzcuROhoaF46KGHZI8bDAarjWGOj4/H8ePH7/njTVmzVrIfNt85\nizovIQSmTZuGhQsXYseOHQCArKwsFBUVKTIuY8+ePTZ/zda8/fbb2Lp1Kx588EGbvaZGo7nrvCrO\nsqKWcKVPbbZv3z44OTlh6dKl0mORkZEYOXKk7HlGoxGjR4+GVquFVqvFf/7zHwANK8kxY8Zg6tSp\nCAoKwurVq/HBBx9Ar9cjMjISly5dAgAsWLAAy5cvx/DhwxEUFASDwYD58+cjNDQUCxculF4nICAA\nJSUlMBqNCAkJwdKlSxEeHo4JEyaguroaQMMgqkceeQQ6nQ6jR4/G2bNnATSszCMiIhAdHY0xY8YA\nALKzsxEXF4eYmBhERUXhwoULzf4Otm/fjsjISERERGD16tUAgHXr1uHw4cNYtGgRVq1aJXu+RqNB\nRUUFHn/8cYSEhOCnP/2p9LEXX3wRer0eEREReOKJJ6TH4+PjsXr1asTFxWHIkCE4dOgQgIYRFXPm\nzEFoaCimT5+OqqoqAA0TShcsWICIiAhERkZi06ZN0rF27tzZ7DitfX/ulJqaimHDhiEnJwfHjx9H\nfHw8dDodJk6ciMLCwmbPp07ElpcJU+e2adMm8cwzz7T4sZycHGlTmcrKSmkzh3PnzgmdTieEEGLf\nvn3C1dVVFBYWitu3bwtfX1/xwgsvSMdesWKFEEKI+fPni8TERCGEEP/85z+Fi4uLOHnypKivrxda\nrVZkZmYKIX7cLCQnJ0c4OjpKj8+aNUt8+OGHQgghxo4dK86fPy+EEOK7774TY8eOFUIIERERIfLz\n84UQQty4cUMIIcTTTz8tPvroIyFEw6YeVVVVsq8xLy9PDBw4UBQXF4u6ujoxduxYsWvXLiGEEPHx\n8eL48ePN/l727dsn+vXrJ/Ly8kR9fb0YPny4OHTokBBCiJKSEul5P/vZz8Tu3bulY/3yl78UQjSM\nJRk3bpwQomEDl+TkZCFEw+gKR0dHcfz4cXHs2DHx8MMPS8dq/HpaO87dvj+TJ08Whw8fFlqtVuTm\n5oqamhoxfPhwUVxcLIQQYseOHWLRokXNvk7qPNjeoTZra7ugpqYGTz31FDIzM+Hg4IDz589LH4uN\njYWXlxcAIDg4GBMmTAAAhIeHY9++fdLrNPaWw8PD4e3tLQ3jCgsLg9FoRGRkpOw1AwMDpce0Wi2M\nRiNu3bqFb7/9VjYOt6amBgAwYsQIzJ8/H7NmzcL06dMBAMOHD8f//M//4MqVK5g+fTqCg4Nlr5Ga\nmoqEhAS4u7sDAObNm4cDBw7gscceA9D6xiV6vR6+vr4AgOjoaBiNRowYMQLffPMNXnnlFVRWVqKk\npATh4eGYPHkyAEg1DRs2DEajEQBw8OBB/Pd//zcASKt6AAgKCsKlS5fwi1/8ApMmTcL48eOl127p\nOHf7/pw+fRpPPPEEvvzyS3h7e+PkyZPIzs6WZrObTCbpa6HOie0darOwsLA2vTH46quvwsfHB1lZ\nWTh27Jhsr84ePXpIf+7WrZt0v1u3bqirq5M+5uTk1Ow5LT2vpeM6ODjAZDKhvr4e/fv3R3p6unTL\nzs4GALzxxht46aWXkJubC61Wi5KSEiQmJmL37t3o1asXfvKTn0g/hBppNBpZsAshZD8IW/uh2FJt\n1dXVePLJJ/HJJ58gKysLS5YskVpSd36Og4OD7Ott6QeLq6srMjMzER8fjzfffBOLFy++63Hu9v3x\n8fFBr169pDHHQgiEhYVJf39ZWVnYu3dvi18ndQ4MfWqzsWPH4vbt27LdmLKysqRecaPy8nJ4e3sD\nAN5//32YTCab1gk0hJWLiwsCAwPx8ccfS49lZWUBaOj16/V6rF27Fh4eHrhy5QpycnIQEBCAp59+\nGo899hhOnDghO2ZsbCz279+P69evw2QyYceOHdL7Ae3VGPDu7u6oqKjAzp07zX7O6NGj8be//Q0A\ncPLkSelraaxn+vTpePHFF5Genn7X49zt++Pq6orPP/8cv/71r7F//34MGTIE165dkzYtr62txalT\np9r/BZPdYOhTu3z66af46quvEBwcjPDwcDz//PPSph2NK93ly5fjvffeQ3R0NM6ePSvbrau11XDT\nM1HasoJu+vkt3f/oo4/w9ttvIzo6GuHh4fjss88AAKtWrZLekB0xYgQiIyORkpKCiIgIxMTEIDs7\nu9kmKT4+PtiwYQMSEhIQHR0NnU5n9hTH1s6wcXV1xZIlSxAeHo6JEyfedRvLxs9ftmwZKioqEBoa\nihdeeAE6nQ5AwzajCQkJiImJwc9+9jO8/PLLdz2Oue+Pp6cnPv/8czz55JPIzMzExx9/jF/96leI\njo5GTExMi2/8UufB0cpERCrClT4RkYow9ImIVIShT0SkIgx9IiIVYegTEakIQ5+ISEUY+kREKsLQ\nJyJSkf8HBD7t5b2s2TMAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 17 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our claim was that this slope and intercept minimize the sum of squared error:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "fitted = X.dot(B)\n", "errors = psychopathy - fitted\n", "np.sum(errors **2)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$252.92560645$$" ], "metadata": {}, "output_type": "pyout", "png": "iVBORw0KGgoAAAANSUhEUgAAAHUAAAAPBAMAAAAok50oAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpm7MhCriUTv3c12\nVGZoascqAAACCUlEQVQoFZ2SvYsTQRiHn83HbXJJNouljeEU78QiC1f40bid7WKhxCsS/DhQi4RD\nvCaQ/QeEXKGgcmAvYgpRG00MioiCQeyNhYUohkPihYtnfGcmuaR2YGd/72/eZ3fmnRdrYdkn9Tq2\nhFZAfLHtmkDbsP7G50z1kCzZ1esyx3IG4Cz2gORo1DJKlgpda2PGZiWMNHjOeZf4M45KwrxnAC7D\nHVI3XmKULF2BizM2T4l6Th27QarCsiS89QzAXSj7c+IYJe8h1IKp7WyJFy2R6dNzRWJ999AAzWDC\naiXH/QUfwqkdrUt+tkRilxMKxU5NWAkeBXM/274IUTLdk/92daDt7OPqEsUciT8MT7dD+KLYMZDY\nJhpYA8kWJaMGzZYOtF18QbJzoEJsOz7sskm8JewYUGWTcVIeo1KtjPrvxC72idwy7CjgSGAj7Bjg\noFLkXcaKH69qnUmQd7P3cQbFitrzDvTCT2NWAekSXJOdhlrp71AOpnaygfM3myOzqyrR+1pRrAH4\niOWKmQ+0MuwTprbU2RnIAZy+qnPv5upq+XbLANJgaTcHD1SrpWUffHNjWybQttxrpJ6uI831Wc4r\nCUkPvcL+9bWr7MPZMKrY4F2w0p2xOUahw0MKPtFKfFPYrGcAmqPRDtbihcCo+RL2wnG5pT1b4ktw\n7r1MrB32pTdO/e5qQIz/Hv8ASJq1UDNm3SwAAAAASUVORK5CYII=\n", "prompt_number": 18, "text": [ "252.92560645" ] } ], "prompt_number": 18 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Is this sum of squared errors smaller than our earlier guess of an intercept of 10 and a slope of 0.9?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "fitted = X.dot([10, 0.9])\n", "errors = psychopathy - fitted\n", "np.sum(errors **2)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$255.75076072$$" ], "metadata": {}, "output_type": "pyout", "png": "iVBORw0KGgoAAAANSUhEUgAAAHUAAAAOBAMAAADjz06NAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpm7MhCriUTv3c12\nVGZoascqAAAB0klEQVQoFV2Qv2vUYByHn9xdLr0fpofgoFOtwxUcPBAsdMrs0kyCOCQUXXSwiJjB\nov0DBCvUQYugCC4OvUFws2dBRVyC/8C5ulgVbA8H4+f9vj0HM3zyPt9Pnrx5QzB/NqHzrrHgAzi/\nUhSrNgn7uwmNF/2igOJDgrHVZnGBaEKrqkY+5O5VVZXbpFOGz4iF+1xarw0xttosrsAjOrfe4kPu\nSWh6+Aq3qfUg5zX1AcZWm8VjuJ40pUjw1yrc9HAOxmkEURL/cJ2x1Waxk/7vwszIu/twv5Rzj/qG\ncw9ZtVlutJ02v+0mWDimro/QJPwpdyQeMPtqbYEpq9a1nSpmDqinwcSHje9I1yQ4gGV9YvSF7A2t\ncsqqzVK2B2699C8Iv7s1S9NnM8h+UXt4yL721il78kxPNwu6Qz85pn3dN3+C2afEk9Czr83q5nAN\n7q5bOMu90kDnHZfwHFpD4j94th2dBZ8JelvaMrVwkywHA0l7KaF+r/5zPHEvEbvaLBpzdHtz8BIL\njVnWyuAjnIaGtjui827g2dVmcaK4cZWjxA98ZDrrWKVN2mW4KU8ui1ws8exqs9ipqt8E/ZXURzuH\n4zqjTcK194mWT+RG85fBs6vN+gt1HLUZeQ63dAAAAABJRU5ErkJggg==\n", "prompt_number": 19, "text": [ "255.75076072" ] } ], "prompt_number": 19 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Contrasts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can combine the values in the $B$ vector in different ways by using a *contrast* vector. A contrast vector $C$ is a vector of weights $c_1, c_2 ... c_p$ for each value in the $B$ vector. Assume that all vectors we've defined up until now are column vectors. Then a scalar value that is a linear combination of the $B$ values can be written:\n", "\n", "$$\n", "C^T B\n", "$$\n", "\n", "We now return to our original question, whether clamminess of handshake predicts psychopathy score.\n", "\n", "If clamminess does predict psychopathy, then we would expect the slope of the best fit line between `clammy` and `psychopathy` would be different from zero.\n", "\n", "In our model, we have two predictors, the column of ones and `clammy`. $p = 2$ and $B$ is length 2. We could choose just the second of the values in $B$ (therefore $b_1$) with a contrast:\n", "\n", "$$\n", "\\left[\n", "\\begin{array}{C}\n", "1 \\\\\n", "0 \\\\\n", "\\end{array}\n", "\\right]\n", "$$\n", "\n", "This is the slope relating `clammy` to `psychopathy`. Now we might be interested if our *estimate* of this slope is different from zero." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To test whether the estimate is different from zero, we can divide the estimate by the variability of the estimate. This gives us an idea of how far the estimate is from zero, in terms of the variability of the estimate. We won't go into the estimate of the variability here though, we'll just assume it (the formula is in the code below). The estimate divided by the variability of the estimate gives us a t statistic." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With that introduction, here's how to do the estimation and a t statistic given the data $y$, the design $\\bf{X}$, and a contrast vector $C$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def t_stat(Y, X, C):\n", " \"\"\" betas, t statistic and significance test given data, design matrix, contrast\n", " \n", " This is OLS estimation; we assume the errors to have independent\n", " and identical normal distributions around zero for each $i$ in \n", " $\\epsilon_i$ (i.i.d).\n", " \"\"\"\n", " # Make sure X, Y, C are all arrays\n", " Y = np.asarray(Y)\n", " X = np.asarray(X)\n", " C = np.atleast_2d(C)\n", " # Calculate the parameters\n", " B = npl.pinv(X).dot(Y)\n", " # The fitted values\n", " fitted = X.dot(B)\n", " # Residual sum of squares\n", " RSS = ((Y - fitted)**2).sum(axis=0)\n", " # Degrees for freedom is the number of observations n minus the number\n", " # of independent regressors we have used. If all the regressor columns\n", " # in X are independent then the (matrix rank of X) == p\n", " # (where p the number of columns in X). If there is one column that can\n", " # expressed as a linear sum of of the other columns then (matrix rank of X)\n", " # will be p - 1 - and so on.\n", " df = X.shape[0] - npl.matrix_rank(X)\n", " # Mean residual sum of squares\n", " MRSS = RSS / df\n", " # calculate bottom half of t statistic\n", " SE = np.sqrt(MRSS * C.dot(npl.pinv(X.T.dot(X)).dot(C.T)))\n", " t = C.dot(B)/SE\n", " # Get p value for t value using t distribution function\n", " ltp = t_dist(df).cdf(t) # lower tail p\n", " p = 1 - ltp # upper tail p\n", " return B, t, df, p" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 20 }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, does `clammy` predict `psychopathy`? If it does not, then our estimate of the slope will not be convincingly different from 0. The t test divides our estimate of the slope by the error in the estimate; large values mean that the slope is large compared to the error in the estimate." ] }, { "cell_type": "code", "collapsed": false, "input": [ "X = np.column_stack((np.ones(12), clammy))\n", "Y = np.asarray(psychopathy)\n", "B, t, df, p = t_stat(Y, X, [0, 1])\n", "t, p" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 21, "text": [ "(array([[ 1.91438925]]), array([[ 0.04229476]]))" ] } ], "prompt_number": 21 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Dummy coding and the general linear model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So far we have been doing *multiple regression*. That is, all the columns (except the column of ones) are continuous vectors of numbers predicting our outcome data `psychopathy`. These type of predictors are often called *covariates*.\n", "\n", "It turns out we can use this same framework to express the fact that different observations come from different groups.\n", "\n", "That's the same as saying that we can express *analysis of variance* designs using this same notation.\n", "\n", "To do this, we use columns of *dummy variables*." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's say we get some new and interesting information. The first 4 students come from Berkeley, the second set of 4 come from Stanford, and the last set of 4 come from MIT. Maybe the student's college predicts if they are a psychopath?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How do we express this information? Let's forget about the clamminess score for now and just use the school information. Our model might be that we can can best predict the psychopathy scores by taking a representative psychopathy score for the relevant school:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "y_i = \\mu_{Berkeley} + e_i \\space\\mbox{if}\\space 1 \\le i \\le 4\n", "$$\n", "\n", "$$\n", "y_i = \\mu_{Stanford} + e_i \\space\\mbox{if}\\space 5 \\le i \\le 8\n", "$$\n", "\n", "$$\n", "y_i = \\mu_{MIT} + e_i \\space\\mbox{if}\\space 9 \\le i \\le 12\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can code this with predictors in our design using *indicator variables*. The \"Berkeley\" indicator variable vector is 1 when the student is from Berkeley and zero otherwise. Similarly for the other two schools:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "berkeley_indicator = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0]\n", "stanford_indicator = [0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0]\n", "mit_indicator = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1]\n", "X = np.column_stack((berkeley_indicator,\n", " stanford_indicator,\n", " mit_indicator))\n", "# Show the matrix nicely with LaTeX markup\n", "Matrix(X)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left[\\begin{matrix}1 & 0 & 0\\\\1 & 0 & 0\\\\1 & 0 & 0\\\\1 & 0 & 0\\\\0 & 1 & 0\\\\0 & 1 & 0\\\\0 & 1 & 0\\\\0 & 1 & 0\\\\0 & 0 & 1\\\\0 & 0 & 1\\\\0 & 0 & 1\\\\0 & 0 & 1\\end{matrix}\\right]$$" ], "metadata": {}, "output_type": "pyout", "png": "iVBORw0KGgoAAAANSUhEUgAAAFgAAAErCAMAAABtirTdAAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\nTlMAMquZdlQQQO0wRM0iu+/dZol8bIlZUvoAAATdSURBVHgB7V3hepswDHQTSrs2adKN93/WAYmE\nTmDfvABr9yk/ZpXDF+cszMXRkvTUjY9DWulxufGl9NQdm/7xvBJvug5sh24gflqLc+K5OuL3dsJm\nUXv8bI75QQCMxKdzVyJ+79VqL6fZE94PAGyJ2/PxUCK+fgwMn+8ZYoQtcd+hKRGfx7x5zp2CcA1x\nNxK/dZn0QbiCuO2OgwhvXbOohYMriE/d58D4dGtm5A7+CsTtbah5KcYXJHDFiNNN4+fs5I1TIHAN\n8eU8CHvNpRvCNcTNeIEccxcIwjXE6TJc0h/ZxQJgID4ePrrL4TrLJDnQfvbLYZY3AQzEQrBGG8Sq\nYkgRUqgCGuyWFWFYRPMwLHpXKEnhHImoJ62DXR7/E+J6w/LSvcjrKY24zrD8xLcKRWJ0JDIWbRGu\n0DihI1FCCRCuIQ7DIhrWtU7jus6ls4NY1dlUirYpvTPXMdQFp/5t5i77FXXDKpy9qcZWijAsMgvF\nOwhuoUgXbRF2k1ckxi0UJZQA4Qpit4UifNI6uILYORIhlNbBX4E4DItMTr+NGDsskxgVkcvjip7k\n1CBWgb6nFGFYxgncbfLCsMgFU/QV6Eiki7YIu8krEqMjUUIJEK4gdo5E+KR1cAWxcyRCKK2DvwJx\nGBaZnDAskxJ1kcvjus6ls4NY1fmeUpjdWH0lDwduN/ZhPiX4nhrrDgsUoeiLgqBkaKA7SAFFKMB3\n/6Nc4wLdLTEWocyJSY0LdrfE6DjmxP2Rku/A7pYYHUc1MXY3xM5x1BK77obYOY5aYtd9U+L7mxti\nSMZXUJg8193uxpIalYG5QOw+MTJSJPxMp1Zj190Skx0UNmLsbonZRz5ECuwOxFCEsiAFqXGB7kC8\nwPXXh4JYpQspQIowLKMc+2QFOA6dBgjCsIgcpVtTGBZRKQyLKJHCsExShGExWthwn4XePuOjcYxY\nFRykCMMyyrFPVoRhkc0mdByakTYIw6JqFKQIw6IqhWGZpAjDYrSw4T4LvX3GR+MYsSoYhgWkkJup\nHlwjgHQLwyIah2GRqnTnOBZzLgyLylKQIgyLqpTCsBgtbAjrsQUejYNYFdxUithhGXXeVGO5Sydm\nWBie0vSJEYwYilA0caaA4bbExRJjEcrEJxHBscTFEjPDwnAoZrDEWIQi45xahueIXRHKRHiPGN6f\nZu61ZsTMsDB8R+I/rGFxRSozqXDENTUsvMZlWWNXhDIfEa9xyRBjEcqcmOG5rHCOY06MjmQBz6Rb\n/+Ut5e9ZYTiUuJg8XhrC3x8LYtUupAApwrCMckBWMN/A8PAVmmK5ZZP5BobnFnrmGxjeD3z51sR8\nA8P3J2a+geHzEctawXwDw0FjW6DPfAPDgdhe0sw3MDxLTH0DfIHbdFlM0XK6ha+IL8a454jN4ylt\nVog2JY7/fTzO0KYax0bITGO20cHw7ELPDAnDc8TMkDC8V2D5DsIMCcP3J2aGhOHzEb+8/hgzhBkS\nhoPGv16nWxMzJAwHYrtWMEPC8CxxGBb5aHPMjNX+sZO3GulAFMQqZ0jxv0jBdlAYHjssmgmZ239i\nhoThuTsIMyQM7wcehkVnb1kKVyuiZ2tQZVjssskMCcNzk8e+v43iWeIwLGFYbnlv81ivhDWCjYk3\n+gGddvjBm6Z5W0OBgWP8AZ2mSb8B5j+mG59c1oYAAAAASUVORK5CYII=\n", "prompt_number": 22, "text": [ "\u23a11 0 0\u23a4\n", "\u23a2 \u23a5\n", "\u23a21 0 0\u23a5\n", "\u23a2 \u23a5\n", "\u23a21 0 0\u23a5\n", "\u23a2 \u23a5\n", "\u23a21 0 0\u23a5\n", "\u23a2 \u23a5\n", "\u23a20 1 0\u23a5\n", "\u23a2 \u23a5\n", "\u23a20 1 0\u23a5\n", "\u23a2 \u23a5\n", "\u23a20 1 0\u23a5\n", "\u23a2 \u23a5\n", "\u23a20 1 0\u23a5\n", "\u23a2 \u23a5\n", "\u23a20 0 1\u23a5\n", "\u23a2 \u23a5\n", "\u23a20 0 1\u23a5\n", "\u23a2 \u23a5\n", "\u23a20 0 1\u23a5\n", "\u23a2 \u23a5\n", "\u23a30 0 1\u23a6" ] } ], "prompt_number": 22 }, { "cell_type": "markdown", "metadata": {}, "source": [ "These indicator columns are *dummy variables* where the values code for the group membership." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the $B$ vector will be:\n", "\n", "$$\n", "\\left[\n", "\\begin{array}{B}\n", "\\mu_{Berkeley} \\\\\n", "\\mu_{Stanford} \\\\\n", "\\mu_{MIT} \\\\\n", "\\end{array}\n", "\\right]\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When we estimate these using the least squares method, what estimates will we get for $B$?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "B = npl.pinv(X).dot(psychopathy)\n", "B" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 23, "text": [ "array([ 10.74225, 11.3355 , 18.03425])" ] } ], "prompt_number": 23 }, { "cell_type": "code", "collapsed": false, "input": [ "print(np.mean(psychopathy[:4]))\n", "print(np.mean(psychopathy[4:8]))\n", "print(np.mean(psychopathy[8:]))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "10.74225\n", "11.3355\n", "18.03425\n" ] } ], "prompt_number": 24 }, { "cell_type": "markdown", "metadata": {}, "source": [ "It looks like the MIT students are a bit more psychopathic. Are they more psychopathic than Berkeley and Stanford?\n", "\n", "We can test whether the mean for the MIT students is greater than the mean of (mean for Berkeley, mean for Stanford):" ] }, { "cell_type": "code", "collapsed": false, "input": [ "B, t, df, p = t_stat(psychopathy, X, [-0.5, -0.5, 1])\n", "t, p" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 25, "text": [ "(array([[ 2.34035606]]), array([[ 0.02199731]]))" ] } ], "prompt_number": 25 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ah - yes - just as we suspected." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The model above expresses the effect of group membersip. It is the expression of an analysis of variance (ANOVA) model using $y = XB + \\epsilon$." ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "ANCOVA in the General Linear Model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our formulation $y = XB + \\epsilon$ makes it very easy to add extra regressors to models with group membership. For example, here is a simple ANCOVA model (analysis of covariance)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ANCOVA is a specific term for the case where we have a model with both group membership (ANOVA model) and one or more continuous covariate. \n", "\n", "For example, we can add back our clamminess score to the mix. Does it explain anything once we know which school the student is at?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "X = np.column_stack((berkeley_indicator,\n", " stanford_indicator,\n", " mit_indicator,\n", " clammy))\n", "X" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 26, "text": [ "array([[ 1. , 0. , 0. , 0.389],\n", " [ 1. , 0. , 0. , 0.2 ],\n", " [ 1. , 0. , 0. , 0.241],\n", " [ 1. , 0. , 0. , 0.463],\n", " [ 0. , 1. , 0. , 4.585],\n", " [ 0. , 1. , 0. , 1.097],\n", " [ 0. , 1. , 0. , 1.642],\n", " [ 0. , 1. , 0. , 4.972],\n", " [ 0. , 0. , 1. , 7.957],\n", " [ 0. , 0. , 1. , 5.585],\n", " [ 0. , 0. , 1. , 5.527],\n", " [ 0. , 0. , 1. , 6.964]])" ] } ], "prompt_number": 26 }, { "cell_type": "markdown", "metadata": {}, "source": [ "No - it looks like the MIT students also have clammy hands, and once we know that the student is from MIT, the clammy score is not as useful." ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Displaying the design matrix as an image" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can show the design as an image, by scaling the values with columns.\n", "\n", "We scale within columns because we care more about the seeing variation within the regressor than between regressors. For example, if we have a regressor varying between 0 and 1, and another between 0 and 1000, without scaling, the column with the larger numbers will swamp the variation in the column with the smaller numbers." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def scale_design_mtx(X):\n", " \"\"\"utility to scale the design matrix for display\n", "\n", " This scales the columns to their own range so we can see the variations \n", " across the column for all the columns, regardless of the scaling of the \n", " column.\n", " \"\"\"\n", " mi, ma = X.min(axis=0), X.max(axis=0)\n", " col_neq = (ma - mi) > 1.e-8\n", " Xs = np.ones_like(X)\n", " mi = mi[col_neq]\n", " ma = ma[col_neq]\n", " Xs[:,col_neq] = (X[:,col_neq] - mi)/(ma - mi)\n", " return Xs" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 27 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we can display this scaled design with a title and some default image display parameters:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def show_design(X, design_title):\n", " \"\"\" Show the design matrix nicely \"\"\"\n", " plt.imshow(scale_design_mtx(X),\n", " interpolation='nearest',\n", " cmap='gray') # Gray colormap\n", " plt.title(design_title)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 28 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can then see our ANCOVA design above at a glance:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "show_design(X, 'ANCOVA')" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAG4AAAEKCAYAAADzSWADAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEpNJREFUeJztnXtMU+f/x99FCm0pRQpClYtGKCgXuUztplNRh4pR1InJ\nlAleplvcH+qc0SXqNHET4jeZOpa5GGLwAl7YMtlEvGyiMGbYBJ3GRXDCLCiY4RVRufTz+2M/O2sL\nwqEoj3xeSRP6nOd8zufpi+ecp+ecp0dGRARGOBxedgKMNFicoLA4QWFxgsLiBIXFCQqLE5RuKS4m\nJgZarRaNjY3msnnz5sHBwQG//fabuezKlStwcLBswtGjRzF69GhoNBp4eXkhJiYGP/zwg3l5VVUV\nEhMT4enpCbVaDYPBgMOHDwMAqqurIZfLcfXqVaucZsyYgZUrV1rkI5fLUVNTY7d2dwjqZlRUVJBS\nqaTg4GA6ePCguTw5OZk8PDxowoQJ5rLy8nKSyWTm9wcPHiSNRkPp6el07949IiI6deoULVq0iIiI\n6urqqH///rRgwQKqra2lR48eUVZWFmk0GsrOziYiookTJ9L69estcqqrqyNnZ2e6ePEiERHV19eT\nWq2mIUOG0ObNm7vmg3gO3U7chg0baOrUqbRx40aaMmWKuXzevHn00UcfkU6no1OnThGRpTiTyUR+\nfn70v//9r9XYa9asofDwcKvy1NRU6t+/PxERZWZmUkBAgMXyr776iqKjo83vMzIyKDw8nPbs2UNh\nYWGS29oZup24gIAA2rNnD5WVlZFcLqebN28S0b/i1qxZQ9u2baM333yTiCzF/fnnnySTyaiysrLV\n2AaDwao3ERFdvXqVZDIZlZWVUUNDA7m5uVFhYaF5+euvv05bt241vx83bhxt3LiR7t27RwqFgs6e\nPWuXtneEbnWMKywsRHV1NeLj46HX6xESEoK9e/eal8tkMrz//vu4du0a8vLyLNatq6sDAPTt27fV\n+HV1dTaXPyn7559/oFQqMWvWLOzatQsAUF5ejpKSEsyZMwcAcO3aNeTn52PWrFlwdXXFxIkTzXVf\nJN1KXEZGBiZMmABXV1cAwKxZs5CRkQEAoP8/F+7k5IS1a9di7dq1kMlk5nU9PDwAADdu3Gg1vqen\nJ65fv25V/mQdT09PAEBycjIOHjyIx48fY/fu3Zg0aZJ52e7duxEWFoagoCBzjpmZmWhubu5U2zvM\nC+/jrdDQ0EAajYbUajXpdDrS6XTk7u5ODg4OdP78efOukoioqamJAgMDKTU11eIY5+/v365jnMlk\nsihPSUkxH+OexAoICKD9+/fTwIED6dtvvzUvCwoKIqVSac7R09OTZDIZHTp0yI6fxvPpNuIyMzNJ\nq9WS0Wik2tpaqq2tpZqaGho9ejStWLHCQhwR0d69e0mr1VqMKrOzs8nNzY127txJd+/epZaWFioo\nKKDFixcT0b+jQ39/f5o/fz7V1NTQw4cPKTMzkzQaDR04cMAinw0bNlD//v3Jw8ODGhsbiYioqKiI\nHB0d6eLFixY5JiYm0syZM1/Ap/Qf3UbcpEmT6OOPP7YqP3DgAOl0Opo7dy6tXbvWXG4ymSgsLIwc\nHBws6ufl5dGoUaNIrVZTnz59aOzYsZSbm2tefu3aNZo9ezZptVpycXGh4cOHU05OjtV2KyoqyMHB\ngZYsWWIu++CDDyghIcGqbnFxMSkUCrp9+7aktktBRsQXUkWkWw1OmPbD4gSFxQkKixMUx64KHBkZ\nifPnz3dV+FcKKePDLhtVymSy5ya0fv16rF+/vl2xuhsODg5wdGz7/765ufm5dRobGyWJ412loLA4\nQZEsLi8vD4MGDYJer0dqaqqkGDExMVI3/9J59sq71DpSkXSMa2lpQXBwME6cOAEfHx8MGzYMWVlZ\nGDx48H+B23GMa3eS3fAY5+TkZJc4L/QYV1xcjMDAQAwYMAByuRzvvPMODh06JCUUIxFJ4qqrq+Hn\n52d+7+vri+rqarslxTwfSeK6466rpyHpC7iPjw+MRqP5vdFohK+vr1W9p7+jxcTECD0YsRcmkwkm\nk6nTcSQNTpqbmxEcHIyffvoJ/fr1w/Dhw3lwIhGpgxNJPc7R0RFpaWmYOHEiWlpasHDhQgtpTNfz\nUk95dSRWd+Nl9zg+cyIoLE5QWJygsDhBYXGCwqPKbgCPKnsQLE5QWJygsDhBYXGCwuIEhcUJCosT\nFBYnKCxOUFicoLA4QWFxgsLiBIXFCQqLExQWJyhdNgf8VcdeV+Wl3iXAPU5QWJygsDhBYXGCwuIE\nRZI4o9GIsWPHIjQ0FGFhYdi2bZu982Keh5Qfubxx4waVlpYSEdH9+/cpKCiILl26ZFFHYmibAOh2\nL5lMZpeX1M9JUo/T6XSIjIwEAKjVagwePNjmj1QzXUenj3GVlZUoLS2FwWCwRz5MO+nUmZP6+nok\nJCRg69atUKvVVst58r41ZK/5FCQxUlNTE6ZMmYK4uDgsW7bMOvArPunDnqe8pHxOksQREZKTk+Hh\n4YEvvvjCdmAW1y5eqLjCwkKMHj0aQ4YMMTdg06ZNmDRp0n+BWVy7eKHi2hWYxbULqeL4zImgsDhB\nYXGCwuIEhW9dkIizs7Nd4jx69EjSetzjBIXFCQqLExQWJygsTlBYnKCwOEFhcYLC4gSFxQkKixMU\nFicoLE5QWJygsDhBYXGCwuIEhcUJSpfeV2kvuuP9mVu2bLFLnGXLlvF9lT0JFicoLE5QWJygdEpc\nS0sLoqKiMHXqVHvlw7STTonbunUrQkJCuuVsmlcdyeKqqqqQm5uL9957z27Ddab9SBa3fPlybN68\nGQ4OfJh8GUiaO/Djjz/Cy8sLUVFRyM/Pt3NKrzbl5eW4cuVKp+NIEldUVIScnBzk5ubi0aNHuHfv\nHpKSkrBr165OJ/Sqo9frodfrze+PHj0qKY6k/dznn38Oo9GIiooK7Nu3D+PGjWNpLxi7HKB4VPni\n6fT8uDFjxmDMmDH2yIXpADwkFBQWJygsTlBYnKDwFXCJ+Pn52SWO0WjkK+A9CRYnKCxOUFicoLA4\nQWFxgsLiBIXFCQqLExQWJygsTlBYnKCwOEFhcYLC4gSFxQkKixMUFico/NwBiWi1WrvEMRqNktbj\nHicoLE5QWJygsDhBkSzuzp07SEhIwODBgxESEoIzZ87YMy/mOUgeVS5duhSTJ09GdnY2mpub8eDB\nA3vmxTwHSeLu3r2LgoICZGRk/BvE0RFubm52TYxpG0m7yoqKCvTp0wfz589HdHQ0Fi1ahIaGBnvn\nxrSBJHHNzc0oKSnBkiVLUFJSAhcXF6SkpNg7N6YNJInz9fWFr68vhg0bBgBISEhASUmJXRN7Vamv\nr0dNTY35JRVJ4nQ6Hfz8/FBWVgYAOHHiBEJDQyUn0ZNQq9XQ6XTml1Qkjyq//PJLJCYmorGxEQEB\nAdi5c6fkJJiOw/PjJBIREWGXOOfPn+f5cT0JFicoLE5QWJygsDhB4VGlRNLT0+0SZ+HChTyq7Emw\nOEFhcYLC4gSFxQkKixMUFicoLE5QWJygsDhBYXGCwuIEhcUJCosTFBYnKCxOUFicoAhxBdxedMdH\ngspkMr4C3pNgcYLC4gSFxQmKZHGbNm1CaGgowsPDMWfOHDx+/NieeTHPQZK4yspK7NixAyUlJbhw\n4QJaWlqwb98+e+fGtIGk+XEajQZyuRwNDQ3o1asXGhoa4OPjY+/cmDaQ1OO0Wi1WrFgBf39/9OvX\nD71798Zbb71l79yYNpAk7q+//sKWLVtQWVmJ69evo76+Hnv37rV3bkwbSNpV/v777xgxYgQ8PDwA\nAG+//TaKioqQmJho1+ReRfLz85Gfn9/5QCSBc+fOUWhoKDU0NJDJZKKkpCRKS0uzqAOg2726I1Lz\nkrSrjIiIQFJSEoYOHYohQ4YAABYvXiwlFCMRPsn8kuGTzD0MFicoLE5QWJygsDhB4ecOSCQnJ+el\nbp97nKCwOEFhcYLC4gSFxQkKixMUFicoLE5QWJygsDhBYXGCwuIEhcUJCosTFBYnKCxOUFicoLA4\nQeFbFyRiMple6va5xwkKixMUFicobYpbsGABvL29ER4ebi67desWYmNjERQUhAkTJuDOnTtdniRj\nTZvi5s+fj7y8PIuylJQUxMbGoqysDOPHj+fHSL8k2hQ3atQouLu7W5Tl5OQgOTkZAJCcnIzvv/++\n67JjWqXDx7ja2lp4e3sDALy9vVFbW2v3pJjn06nBiUwm65YTGHsCHf4C7u3tjZqaGuh0Oty4cQNe\nXl5dkdcry8WLF3Hx4sVOx+mwuPj4eGRkZGDVqlXIyMjA9OnTO51ETyIsLAxhYWHm9/v375cUp81d\n5ezZszFixAhcvnwZfn5+2LlzJ1avXo3jx48jKCgIP//8M1avXi1pw0znaLPHZWVl2Sw/ceJElyTD\ntB8+cyIoLE5QWJygsDhBYXGCwlfAJRIZGflSt889TlBYnKCwOEFhcYLC4gSFxQkKixMUFicoLE5Q\nWJygsDhBYXGCwuIEhcUJSo8S196HEbWn3pkzZ+xSRyosTmI9FsdIgsWJih0fhWZBRETES39enAiv\niIgISZ9vlz2GjOlaeFcpKCxOUOwmrq1J/Xl5eRg0aBD0ej3c3d0xZMgQREVFYfjw4eY606ZNg5OT\nE5ydnbF06VKr+Pn5+XBycoJcLodSqcTGjRtt5hESEgJHR0colUqUlpbajOPq6gpXV1colUp4e3tj\n27ZtNmPNnz8fKpUKCoUCgYGBNusdO3YMvXr1glKphEKhwOjRo63qPHr0CDqdDs7OzlAoFFiwYIHN\nvNzc3BAVFYWoqKhW22fGXoORlStXUmpqKhERpaSk0KpVq4iIqLm5mQICAqiiooIaGxtJLpdTUVGR\nxbo5OTmkUqmooqKCCgsLSaVS0aVLlyzqnDx5kkaMGEElJSUUFhZmM4fDhw+TwWCgkpISCggIIIPB\nYFXn5MmTFBsbS6WlpUREdP/+fQoKCrLa3uHDh2ncuHFUWlpKZ86coaFDh9qsd/LkSYqLiyMioqam\nJjIYDFRQUGAVa8KECUREVFhYSC4uLlZ1Tp48SVOnTrXZLlvYrce1Nqm/uLgYgYGBGDBgAORyOVxc\nXHDkyBGLddPT06HX6zFgwACMHDkSKpXK5nPFtVqt1Y8JPJvD8uXL4e7uDqVSiTt37tico+7s7Gy+\noVWtVmPw4MG4fv26VazFixcjMjISBoMB9+/fx8CBA63qAUCvXr0AAI2NjWhpaYFWq7WK9aSXRUVF\noaWlxeZPSlEHxol2E9fapP7q6mr4+fmZ68nlcnz99dcYOnQoduzYAQAwGo3w9/c31/Hy8kJ5eblF\nfJlMhqKiIsTFxaGyshKXLl2yyuHZbfn6+qKqqspmnIiICEyePBnHjx9HaWkpDAZDm7E8PDxw7tw5\nq3oymQy//PILlEolNBoNwsPDERISYhXLx8cHkZGR8Pb2hpeXF1xcXNrMy1b7nqZDt6DHxsaipqbG\nqvyzzz6zSuLJpP5nJ/evW7cOly9fxrp16xAbG4tBgwbZ/AGAZ8uio6NhNBpx8+ZNxMTEYPr06Sgr\nK7Na79n/2tbiqFQqfPfdd5g6dSr27dsHtVrdaqz6+npcuHABn376qVW96OhoVFVVQaVSITs7G3Pn\nzkVSUhJiYmKs8jh37hzu3r0Lf39/nD17Fq+99prNvI4cOdJq+57QoR53/PhxXLhwweoVHx9vntQP\nwGJSv4+PD4xGoznG/fv34evriz59+mDGjBkoLi6Gn58f/v77b3OdmzdvIjAw0GLbrq6uUKlU5r+b\nmppw69YtizrPbquqqgo+Pj424zQ1NeGbb76BSqWyOaB4EqupqQkzZ86EQqHAu+++a1Xv6bwSEhKg\nUChw+vTpVvNyc3ODXC5HdXV1q3Hi4uJstu9p7LarfDKpH4DFpP6hQ4eivLwclZWVuHPnDrKyshAf\nH48HDx7g2LFjCA8Px4IFC8x1CgoK0NDQgMTERIv4tbW15h7w8OFDEJHVsSQ+Ph67du0CADQ0NKB3\n797m3ffTcUwmExYuXAitVguNRmMV5+n2LFy4EO7u7hg4cKBVLAC4dOkSbt++DQA4ffo0Hj58iJEj\nR1rUGTNmDNLT0wH8O3psbGy0+md5un3FxcU222dBu4cxz6Guro7Gjx9Per2eYmNj6fbt20REVF1d\nbR6R+fv7k06no4iICOrbty9NmzbNvP6UKVNILpeTk5MTffjhh0REtH37dtq+fTsREaWlpZFGoyFH\nR0eSyWTk5eVF6enpFnWIiPR6PfXq1YtkMhl5e3tb1UlLS6MBAwYQAFKpVKTX6ykyMpJyc3OtYs2Y\nMYMAkEKhoODgYJv1PvnkE1IoFKRQKEilUtnM/Y8//iBPT09ycnIihUJBS5cutdm+0NBQioiIoDfe\neIN+/fXXNj9vPuUlKHzmRFBYnKCwOEFhcYLC4gSFxQkKixMUFico/weGlEy0qR3uxgAAAABJRU5E\nrkJggg==\n", "text": [ "" ] } ], "prompt_number": 29 } ], "metadata": {} } ] }