A Geometric Interpretation of Ordinary Least Squares Regression
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plan for Today:\n",
"\n",
"- Derive the OLS regression from scratch"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import numpy.linalg as la # Easily access the linear algebra api\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"import statsmodels.formula.api as smf"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generating Fake Data \n",
"\n",
"Let's generate some fake data where we know the true coefficient values. It's useful to know the population values when playing around with estimators because we can easily assess if our answers are correct or not."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"np.random.seed(123) # set seed to repoduce results\n",
"N = 100 # set sample size\n",
"x = np.random.normal(0,1,N) # simulate a random variable\n",
"e = np.random.normal(0,1,N) # simulate error\n",
"y = 1 + 3*x + e # generate some y that is a function of x and random error"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Combine $X$ in into a design matrix with a columns of ones, which we'll use to estimate the constant."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 1. , 0.99734545],\n",
" [ 1. , 0.2829785 ],\n",
" [ 1. , -1.50629471],\n",
" [ 1. , -0.57860025]])"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X = np.vstack([np.ones(N),x]).T\n",
"X[1:5,:]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot the data to get a feel for what is going on."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3YAAAEyCAYAAAC2+0LeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X+MnPd9H/j3JzSTbtxct6jVulzblQ41aDhRbNYLXwumh0hxQ8cJYoZOG6XF9YIcoPNdAjSBwYKqgbS4IhDviPRn0qZqYvTauomNWmaMSq3sgAKcGOc2lClFTiz2dI7TeBQktBs6yXmDSPS3f3DXWq1nd2d3np15npnXCyA488zD5/nucEDNW5/v9/Ot1loAAAAYrq+a9wAAAACYjmAHAAAwcIIdAADAwAl2AAAAAyfYAQAADJxgBwAAMHCCHQAAwMAJdgAAAAMn2AEAAAzcy+Y9gL284hWvaHfeeee8hwEAADAXTzzxxOdaa3fsd16vg92dd96Zq1evznsYAAAAc1FVvz7JeaZiAgAADJxgBwAAMHCCHQAAwMAJdgAAAAMn2AEAAAycYAcAADBwgh0AAMDACXYAAAAD1+sNygEAAI7K5WujXHrsep67uZETqys5f+Zkzp5am/ewDkWwAwAAls7la6M88PDT2Xj+VpJkdHMjDzz8dJIMMtyZigkAACydS49d/3Ko27Lx/K1ceuz6nEY0HcEOAABYOs/d3DjQ8b4T7AAAgKVzYnXlQMf7TrADAACWzvkzJ7Ny/NhLjq0cP5bzZ07OaUTT0TwFAABYOlsNUnTFBAAAGLCzp9YGG+R2MhUTAABg4AQ7AACAgRPsAAAABs4aOwAAoLcuXxstTIOToyTYAQAAvXT52igPPPx0Np6/lSQZ3dzIAw8/nSTC3Q6mYgIAAL106bHrXw51Wzaev5VLj12f04j6S7ADAAB66bmbGwc6vswEOwAAoJdOrK4c6PgymzrYVdXJqnpy26/fraof2nHON1fVF7ad8yPT3hcAAFhs58+czMrxYy85tnL8WM6fOTmnEfXX1M1TWmvXk7wxSarqWJJRkg+OOfUXWmvfMe39AACA5bDVIEVXzP113RXzW5L8f621X+/4ugAAMAja83fr7Kk1798Eul5jd1+Sn9nltb9QVU9V1b+vqq/f7QJVdX9VXa2qqzdu3Oh4eAAAcHS22vOPbm6k5cX2/JevjeY9NBZctda6uVDVVyd5LsnXt9Z+a8dr/12SL7XWfr+q3pbkH7bWXrvfNdfX19vVq1c7GR8AABy10xevZDSmY+Pa6ko+duHeOYxoGFQ5d1dVT7TW1vc7r8uK3bcl+cTOUJckrbXfba39/ubjR5Mcr6pXdHhvAACYO+35D06VsxtdBrvvzS7TMKvqlVVVm4/fvHnfz3d4bwAAmDvt+Q/OJuTd6CTYVdXLk/ylJA9vO/bOqnrn5tPvTvLJqnoqyT9Kcl/rag4oAAD0hPb8B6fK2Y1OumK21v7/JH9ix7Gf3Pb4x5P8eBf3AgCAvlrm9vyHXSd3YnVl7LpEVc6D6Xq7AwAA6LWjbtSxjO35t9bJbU2p3Fonl2Tf9+L8mZMv+bOJKudhdL3dAQAA9JZGHUdjmnVyZ0+t5cFzd2dtdSWV2x1EHzx399KF42mp2AEAsDT2CiCCxOFNu05uGaucXVOxAwBgaWjUcTR0A50/wQ4AgKUhgBwN3UDnT7ADAGBpCCBHwzq5+bPGDgCApbHM2xEcNevk5kuwAwBgqQggLCJTMQEAAAZOxQ4AAOiNo95AflEJdgAAQC9sbSC/tdfg1gbySYS7fZiKCQAA9MJeG8izN8EOAADoBRvIH56pmAAAC8Yapfnx3k/nxOpKRmNCnA3k96diBwCwQLbWKI1ubqTlxTVKl6+N5j20I3P52iinL17JXRceyemLV+b2sy7je981G8gfnmAHALBAlm2NUp/C1LK990fh7Km1PHju7qytrqSSrK2u5MFzd6t6TsBUTACABbJsa5T2ClOzDgPTvvemcd5mA/nDUbEDAFggu61FWtQ1Sn0KstO8932qPDJMgh0AwALZbY3SPa+7oxfr0LrWpyA7zfow0ziZlmAHALBAxq1Reseb1vKBJ0YLWQ3qU7ONadaH9anyyDBZYwcAsGB2rlE6ffFKb9ahdW1r/H1Zm3bY9WHa/DMtwQ4AYMEtejVoEZptnD9zMg88/PRLArg2/xyEqZgAAAuuT+vQGE+bf6alYgcAsOCOohqkNX/3FqHyyPwIdgAAC67rdWhbrfm3guJWM5bt99rtzwmDcDQEOwCAJdBlNegwm4IfNgwCk7HGDgCAAzlMMxb7tMHR6izYVdVnqurpqnqyqq6Oeb2q6h9V1bNV9ctV9ee6ujcAALNzmGYsi96ZE+at64rdPa21N7bW1se89m1JXrv56/4k/7TjewMAMAOH2RS8L505L18b5fTFK7nrwiM5ffHKQmzSDsls19i9Pcm/bK21JB+vqtWq+tOttd+c4RgAACam2cd4h2nGcs/r7sh7P/5f0rYdm/U+bX1d5+dzRhe6DHYtyYerqiX5Z621h3a8vpbkN7Y9/+zmsZcEu6q6P7crennNa17T4fAAACbX1xDQFwdpxnL52igfeGL0klBXSd7xptm29z9M05ej5nNGV7qcivlNrbU/l9tTLn+gqv7Hw1yktfZQa229tbZ+xx13dDg8AIDJafbRnXHvZUvy+DM3ZjqOPq7z8zmjK51V7Fpro83ff7uqPpjkzUk+uu2UUZJXb3v+qs1jAAC908cQME/TTBfsy3u5+rXH8ztffH7s8Xnpy3vD8HVSsauql1fV1209TvKtST6547QPJfnrm90x/3ySL1hfBwD0VV+afWyZZ9OPremCo5sbaXlxuuCkY+jLe9nawY7PQl/eG4avq6mYfyrJL1bVU0n+U5JHWmv/oareWVXv3Dzn0SSfTvJskn+e5H/v6N4AAJ07TOfHrm2FuTsvPJIfft+Thw5W05p2umAf3ssk+cLGV1br9jo+C315bxi+TqZittY+neQNY47/5LbHLckPdHE/AICjdpjOj13a2VRjZ1Fplk0/dpsWOLq5kbsuPLLvezPv93LLidWVjMb8LPOsjvXlvWH4qs2z9ryP9fX1dvXqV+x1DgCw8E5fvDI2hGxXSX7t4rcf6voHWTM3yVhWjh/Lg+fu7nUg2RmWk2GMm+VWVU/ssk/4S3S9QTkAAB2YpHnGYStNB10zN2664E5D6OR49tRaHjx3d9ZWV1JJ1lZXhDoWxiw3KAcAYB9blbT95lRNsw7roPu57ZwuuNvYhtDJ8SD778GQCHYAAD0xbqrgdpXba+3WplyHdZgW+9sD0W5TM3VyhPkR7AAAemJcJW3LtGFuu2mbiJw/c3LsWjWdHGF+BDsAgJ7YrWJWST524d7O7nPP6+7Iv/74fxl7fBI6OUL/CHYAAD0xq3b8jz9z40DHx7FWDfpFV0wAgJ6Y1WbVh1ljB/Sbih0AMBcH2UdtWcxqimMfN+oGpiPYAQAzt7P749Y+akmEuxlMcZym+YlADv0k2AEAM3fQfdTY20HD1mErgwI59JdgBwDMnDVe3Tls2DpMZVAgh/7SPAUAmLnd1nJZ43Vwe4Wtrgnk0F+CHQAwc7Pq/rgMZhm2BHLoL8EOAJi5s6fW8uC5u7O2upJKsra6kgfP3W063yHMMmwJ5NBf1tgBAHNhg+tuTNPh8qBmtR0DcHCCHQDAgM06bAnk0E+CHQDAEZjlfm/CFiDYAQB0zH5vwKwJdgDAQpplxWyn3bYgeNf7n0oi3AHdE+wAgIUz74rZblsN3GpN5Q44ErY7AAAWziw27b58bZTTF6/krguP5PTFK7l8bfTl1/baauCoNg8HlptgBwAsnKPetHurIji6uZGWFyuCW+Fu3H5vRzEOgC2CHQCwcI560+79KoJbG7AfqzrScQBsEewAgIUzrmLW5abdk1QEz55ay4/9lTcc6TgAtmieAgA9N8/ujkN11Jt2n1hdyWhMuNtZiZv15uHA8qrW2rzHsKv19fV29erVeQ8DAOZmZ3fH5HbF58Fzd88lHOwXMpclhPbt7wVYXFX1RGttfb/zpp6KWVWvrqrHq+pXq+pXqupvjDnnm6vqC1X15OavH5n2vgCwDGbR3XFS+zUM2e/1RbK1hm5tdSWVZG11ZddQt1f3TICudDEV84Uk72qtfaKqvi7JE1X1kdbar+447xdaa9/Rwf0AYGkcdXfHg9grZJ49tbbv64vm7Km1fX+uee+nByyPqSt2rbXfbK19YvPx7yX5VBL/UgFAB466u+NB7Bcy+xRC+6JPFVdgsXXaFbOq7kxyKsl/HPPyX6iqp6rq31fV1+9xjfur6mpVXb1x40aXwwOAwTnq7o4HsV/I7FMI7QthF5iVzoJdVf3RJB9I8kOttd/d8fInkvyZ1tobkvzjJJd3u05r7aHW2nprbf2OO+7oangAMEgHWct11PYLmX0KoX0h7AKz0sl2B1V1PLdD3Xtbaw/vfH170GutPVpV/6SqXtFa+1wX9weARTbJWq5ZjSPZvXW/1v5f6fyZk2O7Zy5z2AWOxtTBrqoqyU8n+VRr7e/tcs4rk/xWa61V1Ztzu1L4+WnvDQDM1riQuSxbHByGsAvMShcVu9NJ/qckT1fVk5vH/laS1yRJa+0nk3x3kv+tql5IspHkvtbnDfQAgIno+ri/WVdcBW1YTlMHu9baLyapfc758SQ/Pu29AGAIlumL9bJtcdB3gjYsr07W2AEAt/X5i/VRBE5dH/tF0Ibl1el2BwCw7Pq6b9lW4Bzd3EjLi4Hz8rXRVNedtOvj5WujnL54JXddeCSnL16Z+r6MJ2jD8hLsAKBDff1ivVvgfNf7n5oqZE2yxcFRhUq+ku0VYHkJdgDQob5+sd4tWN5qbaqQNck+e32tYi4iewnC8rLGDgA61Nd9y06srmS0S7ibdg3Wfl0f+1rFXES2V4DlJdgBQIf6+sV6XODcbtKQdZgGLLuFynlXMRdVXza0B2ZLsAOAjvXxi/XWeN71/qdya8xWspOErMN2/OxrFRNgkQh2ALAktsLXYUPWYVvpH7aKuUz7AQJMS7ADgCWyM2Stfu3xtJb88PuezKXHru8ZnqZZK3fQKmaf9wME6CNdMQFgyZw9tZaPXbg3f/973pg/eP5Lubnx/ETbEMyy46dOmgAHI9gBwJI6aHiaZSt9nTQBDsZUTAAGwXqr7h00PM2y46dOmgAHI9gB0Ht9XG+1CEHzMOFpVh0/ddIEOBhTMQHovb6tt9oKmqObGxOtTeurWU6tPKizp9by4Lm7s7a6kkqytrqSB8/dPbjwDDArKnYA9F7f1lsdtu1/3/R1M/UtfdwPcKgWocIM7E2wA6D35r3eaueX4nFjSYbZ2EN4Wnx9nMoMdM9UTAB6b55TBsdNu6xdzj1s0Lx8bZTTF6/krguP5PTFK4Ob0km/9W0qM3A0VOwA6IW9porNc8rguC/FLUlt/r7lsEFTNYWj1repzMDREOwAmLtJws28pgzu9uW35XZDj2mD5qKs16O/5j2VGZgNwQ6AuetzuNntS/Ha6ko+duHeqa+vmsJRs3UELAdr7ACYuz6Hm6Ne37db1eSrqqy1oxO2joDloGIHsKT61P68z1PFjnp937hqSpLcas1aOzqj+yksPsEOYAn1rWFH36eKHeWX4q3rvuv9T+VWay95rS/TUQHoP1MxAZZQ39qfL/tUsbOn1vKlHaFuSx+mowLQfyp2AEuoj2vapq2K9Wlq6WH0eToqAP0n2AEsoSGFiEkCW9+mlh5G36ejdmHo4Rugz0zFBFhCR93psStbgW10cyMtLwa2nd0i+za19DDGTUd9x5vWcumx67nrwiM5ffHKoLtkTvp3CcDhdFKxq6q3JvmHSY4l+anW2sUdr39Nkn+Z5E1JPp/ke1prn+ni3gAc3FF3ejysnRWdL/7hCxPtbzfPqaVdVqG2T0ddhCrkdn3eqxBgEUwd7KrqWJKfSPKXknw2yS9V1Ydaa7+67bT/JcnvtNb+bFXdl+T/TPI9094bgMPrW/vzcUFmNzsD27ymlh5l+Fq0INTHdZ0Ai6SLqZhvTvJsa+3TrbU/TPKzSd6+45y3J/m/Nx//2yTfUlXVwb0BFsbla6OcvnhlIabdHca4ILObnYFtXlNLj3IK6KIFod1Cdh/XdQIMURfBbi3Jb2x7/tnNY2PPaa29kOQLSf7EuItV1f1VdbWqrt64caOD4QH0n/VHkweWcYFtXtslHGX4WrQgNJR1nQBD1buumK21h5I8lCTr6+vjN/UBWDCLNu3uMHabTrm6cjwv/5qX7buGbR5TS49yCuiidcns67pOgEXRRbAbJXn1tuev2jw27pzPVtXLkvyx3G6iAkAWb9rdYewWZP7Od379SxqKXHrsen74fU/2IhgcZfhaxCDUt3WdAIuki2D3S0leW1V35XaAuy/JX91xzoeS/M9J/p8k353kSmtNNQ5gU9/2lZvHfmP7BZk+dok86vAlCAEwqeoiX1XV25L8g9ze7uA9rbUfrar/I8nV1tqHquqPJPlXSU4l+a9J7mutfXq/666vr7erV69OPT6AvtsZWpLblZ9ZrBPr81i2O33xytjwu7a6ko9duPclx2yEDcCiqKonWmvr+53XyRq71tqjSR7dcexHtj3+gyR/uYt7ASyiPk276+t6v0mnq/axsgcAR613zVMAllVfpt31db3fpNNV+xpMAeAodbHdAQALpK9t9idtl9/XYAoAR0nFDmDJ7Lf+bNJOj7NexzbpdNVJKns7x37P6+7I48/cmPs0WAA4rE6apxwVzVMAujVpY5T9QltfG6wk+49t3Os79eVnAYBJm6cIdsCudBZcPAfpLDmL6xyVvT67u419p778LAAst5l2xQQWj86Ci6mr9Wd9X8e2VyOaScfYl58FACaheQow1l6dBRmurhqj9LXByiQmHeMQfhYA2CLYAWP1vSLD4UzaWXJW15mHcWPfaSg/CwBsMRUTGGvSPcNmzbq/6UzaWXK/97lPG6of1Lix64oJwNBpngKM1ceuh30c0yLyPgNAf0zaPMVUTGCss6fW8uC5u7O2upLK7Q6B8/5ib93f3i5fG+X0xSu568IjOX3xSi5fGx3qOt5nABgeUzGBXe3VWXAerPvbXZddTL3PADA8gh0wGH1d99cHe1XZ9gp249bSeZ8BYHhMxQQGY8idGKe13zTLw1TZtqp8o5sbaXmxynfP6+5Y2vcZAIZKsAMGY1br/rpaq9bleMYFsO3jOsy+crtV+R5/5kbv1lcCAHszFRMYlKNe99flWrVx1z7M9gCTTLM8f+bk2E6We1XZ9qry9W19JQCwN8EOYJvDrlXbzzSBcZJplnvtK7dboLSWDgAWh2AHsM1+Ieooq267mTSAjauy7RUoD1PlAwD6yRo7gG32Wqs2yVq33UyzhcA9r7sjtePYpAFsv0BpLR0ALAYVO4Bt9qpizaLqttPla6N84IlR2rZjleQdb5psDdx+gdJaOgBYDCp2ANvsVcWapup22K0axoXJluTxZ27se8/kcN0yAYDhUbED2GG3KtY0zUb2am6yl2nCZGIdHQAsC8EOYELThqSDTHvcatLSdnl90orbYQPlPB22QQ0ALLNqbbevDfO3vr7erl69Ou9hAEtikkAxi9Cxs5PlTivHj+Udb1rL48/cWLjwM+5nXzl+TFMXAJZWVT3RWlvf7zwVO6AX5l2lmXSfuVk0Gxm3rm7L2upK7nndHfnAE6Mj2UR93o5qH0EAWHSapwBzN802Al3ZK1DM2m7r5yrJxy7cm8efudGbsXZt2jWFALCsBDtg7voQqvoUKPbrZNmnsXZNF08AOJypgl1VXaqqZ6rql6vqg1W1ust5n6mqp6vqyaqyaA54iS6DyuVro5y+eCV3XXgkpy9embjqt1tw+KqqmVYOk/23Rljk8HPYbSEAYNlNW7H7SJJvaK19Y5L/nOSBPc69p7X2xkkW/gHLpaugMs2UznGBIklutTbzaaF77aW321gXJfzs97MDAON11hWzqr4ryXe31v7amNc+k2S9tfa5g1xTV0xYDuM6IR4/Vnn5V78sX9h4fuJmKqcvXhm7z9za6ko+duHeicbxrvc/lVtj/l2c9BqzMu9mMwDAbMyjK+b3J3nfLq+1JB+uqpbkn7XWHtrtIlV1f5L7k+Q1r3lNh8MD+mrnXmurX3s8v/8HL+TmxvNJJu/6OO2UzrOn1vLD73tyqmvMyiy6cwIAw7FvsKuqn0/yyjEvvbu19nOb57w7yQtJ3rvLZb6ptTaqqj+Z5CNV9Uxr7aPjTtwMfQ8ltyt2E/wMQE9MU0XaHlROX7yS3/ni8y95fZKW9ydWV8ZW7A4ypbOLawAAzNq+a+xaa29prX3DmF9boe77knxHkr/WdpnX2Vobbf7+20k+mOTNnf0EQC90uWXBbtWx0c2NPZuidLH2bJHXrwEAi2varphvTfI3k3xna+2Lu5zz8qr6uq3HSb41ySenuS/QP11uWbBXdWyv0NhF4w3NOwCAIZqqeUpVPZvka5J8fvPQx1tr76yqE0l+qrX2tqr673O7Spfcnvr5b1prPzrJ9TVPgeG468IjGfevSSX5tYvffqBrjWumMk7fGpoAAHRtJs1TWmt/dpfjzyV52+bjTyd5wzT3Afqvy7VpO5up7Pa/n/rW0AQAYF6m3ccOIEn3a9POnlrLxy7cm1+7+O1ZW+ANuQEAutDldgfAEttZZdutK+ZhOmeeP3PyK6ZmLnJDE3vUAQAHJdgBndlvb7Wda+cm3Z9u0tC4CA7yHgmAAMCWqZqnHDXNU2CxnL54Zew6PE1QXjTpezSuwczK8WM6eALAgpm0eYo1dsDM7NbsRBOUF036HnW5vQQAMHyCHTAzuzU70QTlRZO+R0IyALCdYAfMTNedMxfRpO+RkAwAbKd5CjAzR9kEZVEaiUz6Hi1bp1AAYG+ap8ASWJTQs5tlbSSy6H+vAMDkzVNU7GDBHXaLgSHZq5HIovyM4+y3vQQAsDyssYMFtwzdEzUSAQCWnWAHC24ZQo9GIgDAshPsYMEtQ+jRbRMAWHaCHfTY5WujnL54JXddeCSnL17J5WujA19jGULP2VNrefDc3VlbXUklWVtdWfjGKQAA22meAj3VVdOTo9xioE80EgEAlplgBz3VZadHoQcAYLGZigk9tQxNTwAA6IZgBz21DE1PAADohmAHPbUMTU8AAOiGNXbQU8vS9AQAgOkJdtBjmp4AADAJUzEBAAAGTrADAAAYOMEOAABg4KyxA9jh8rWRpjUAwKAIdgDbXL42ygMPP52N528lSUY3N/LAw08niXAHAPSWYMcgqahwVC49dv3LoW7LxvO3cumx6z5jAEBvTbXGrqr+TlWNqurJzV9v2+W8t1bV9ap6tqouTHNP2KqojG5upOXFisrla6N5D40F8NzNjQMdBwDogy6ap/z91tobN389uvPFqjqW5CeSfFuS1yf53qp6fQf3ZUntVVEZosvXRjl98UruuvBITl+8IqDO2YnVlQMdBwDog1l0xXxzkmdba59urf1hkp9N8vYZ3JcFtUgVFdXH/jl/5mRWjh97ybGV48dy/szJOY0IAGB/XQS7H6yqX66q91TVHx/z+lqS39j2/LObx8aqqvur6mpVXb1x40YHw2PRLFJFZdGqj4vg7Km1PHju7qytrqSSrK2u5MFzd1tfBwD02r7NU6rq55O8csxL707yT5P83SRt8/cfS/L90wyotfZQkoeSZH19vU1zLRbT+TMnX9K1MBluRWWRqo+L5OypNUEOABiUfYNda+0tk1yoqv55kn835qVRkldve/6qzWNwKFtfuBehK+aJ1ZWMxoS4IVYfAQCYn6m2O6iqP91a+83Np9+V5JNjTvulJK+tqrtyO9Ddl+SvTnNfWJSKyiJVHwEAmJ9p97H7v6rqjbk9FfMzSf7XJKmqE0l+qrX2ttbaC1X1g0keS3IsyXtaa78y5X1hISxS9REAgPmp1vq7jG19fb1dvXp13sMAAACYi6p6orW2vt95s9juAAAAgCMk2AEAAAycYAcAADBwgh0AAMDACXYAAAADJ9gBAAAMnGAHAAAwcNNuUM6CuHxtZJNsAAAYKMGOXL42ygMPP52N528lSUY3N/LAw08niXAHAAADYComufTY9S+Hui0bz9/Kpceuz2lEAADAQQh25LmbGwc6DgAA9ItgR06srhzoOAAA0C+CHTl/5mRWjh97ybGV48dy/szJOY0IAAA4CM1T+HKDFF0xAQBgmAQ7ktwOd4IcAAAMk6mYAAAAAyfYAQAADJxgBwAAMHCCHQAAwMAJdgAAAAMn2AEAAAycYAcAADBwgh0AAMDACXYAAAADJ9gBAAAMnGAHAAAwcIIdAADAwL1smj9cVe9LcnLz6WqSm621N4457zNJfi/JrSQvtNbWp7kvAAAAL5oq2LXWvmfrcVX9WJIv7HH6Pa21z01zPwAAAL7SVMFuS1VVkr+S5N4urgcAAMDkulpj9xeT/FZr7f/d5fWW5MNV9URV3b/Xharq/qq6WlVXb9y40dHwAAAAFte+Fbuq+vkkrxzz0rtbaz+3+fh7k/zMHpf5ptbaqKr+ZJKPVNUzrbWPjjuxtfZQkoeSZH19ve03PgAAgGW3b7Brrb1lr9er6mVJziV50x7XGG3+/ttV9cEkb04yNtgBAABwMF1MxXxLkmdaa58d92JVvbyqvm7rcZJvTfLJDu4LAABAugl292XHNMyqOlFVj24+/VNJfrGqnkryn5I80lr7Dx3cFwAAgHTQFbO19n1jjj2X5G2bjz+d5A3T3gcAAIDxuuqKCQAAwJwIdgAAAAMn2AEAAAycYAcAADBwgh0AAMDATd0Vc5lcvjbKpceu57mbGzmxupLzZ07m7Km1eQ8LAABYcoLdhC5fG+WBh5/OxvO3kiSjmxt54OGnk0S4AwAA5spUzAldeuz6l0Pdlo3nb+XSY9fnNCIAAIDbBLsJPXdz40DHAQAAZkWwm9CJ1ZUDHQcAAJgVwW5C58+czMrxYy85tnL8WM6fOTmnEQEAANymecqEthqk6IoJAAD0jWB3AGdPrQlyAABA75iKCQAAMHBMD0uRAAADlElEQVSCHQAAwMAJdgAAAAMn2AEAAAycYAcAADBwgh0AAMDACXYAAAADJ9gBAAAMXLXW5j2GXVXVjSS/vs9pr0jyuRkMh8Xmc0QXfI7ois8SXfA5ois+S/P1Z1prd+x3Uq+D3SSq6mprbX3e42DYfI7ogs8RXfFZogs+R3TFZ2kYTMUEAAAYOMEOAABg4BYh2D007wGwEHyO6ILPEV3xWaILPkd0xWdpAAa/xg4AAGDZLULFDgAAYKkJdgAAAAM3+GBXVX+3qn65qp6sqg9X1Yl5j4lhqqpLVfXM5ufpg1W1Ou8xMTxV9Zer6leq6ktVpTU0B1JVb62q61X1bFVdmPd4GKaqek9V/XZVfXLeY2G4qurVVfV4Vf3q5n/X/sa8x8TeBh/sklxqrX1ja+2NSf5dkh+Z94AYrI8k+YbW2jcm+c9JHpjzeBimTyY5l+Sj8x4Iw1JVx5L8RJJvS/L6JN9bVa+f76gYqH+R5K3zHgSD90KSd7XWXp/kzyf5Af8m9dvgg11r7Xe3PX15Et1gOJTW2odbay9sPv14klfNczwMU2vtU6216/MeB4P05iTPttY+3Vr7wyQ/m+Ttcx4TA9Ra+2iS/zrvcTBsrbXfbK19YvPx7yX5VJK1+Y6Kvbxs3gPoQlX9aJK/nuQLSe6Z83BYDN+f5H3zHgSwVNaS/Ma2559N8j/MaSwAX1ZVdyY5leQ/znck7GUQwa6qfj7JK8e89O7W2s+11t6d5N1V9UCSH0zyt2c6QAZjv8/S5jnvzu3pB++d5dgYjkk+RwCwCKrqjyb5QJIf2jFTjp4ZRLBrrb1lwlPfm+TRCHbsYr/PUlV9X5LvSPItzSaP7OIA/ybBQYySvHrb81dtHgOYi6o6ntuh7r2ttYfnPR72Nvg1dlX12m1P357kmXmNhWGrqrcm+ZtJvrO19sV5jwdYOr+U5LVVdVdVfXWS+5J8aM5jApZUVVWSn07yqdba35v3eNhfDb0oUVUfSHIyyZeS/HqSd7bW/B9ODqyqnk3yNUk+v3no4621d85xSAxQVX1Xkn+c5I4kN5M82Vo7M99RMRRV9bYk/yDJsSTvaa396JyHxABV1c8k+eYkr0jyW0n+dmvtp+c6KAanqr4pyS8keTq3v2cnyd9qrT06v1Gxl8EHOwAAgGU3+KmYAAAAy06wAwAAGDjBDgAAYOAEOwAAgIET7AAAAAZOsAMAABg4wQ4AAGDg/hs73SfQREs5qwAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(15,5))\n",
"plt.scatter(X[:,1],y)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Geometric Interpretation of Linear Regression\n",
"\n",
"The following discussion will tie together what we've learned up to this point about linear algebra and apply it one of the most commonly used statistical learning model in our toolkit: linear regression."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Design Matrix: Column Space of a Matrix of Explanatory Variables"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Say we had a some variables in matrix $\\textbf{X}$ that we wanted to use to predict some outcome $\\textbf{y}$. That is, we have some variables that we want use to help train a model to predict future values of some outcome, like housing prices or the number of crimes.\n",
"\n",
"$\\textbf{X}_{n\\times p}$ is a matrix with dimensions $n$ (the number of observations) and $p$ (the number of variables). We call this matrix a \"design matrix\". It holds all of the variables we intend to use in our prediction. \n",
"\n",
"First, as we've learned, if the column space of $\\textbf{X}$ — $C(\\textbf{X})$ — is of full rank (meaning each column vector is linearly independent) then we can think of $C(\\textbf{X})$ as spanning the space in $\\Re^p$ (which we'll represent for the sake of convenience as a two dimensional plane in $\\Re^2$)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Outcome as a Vector "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we've seen in previous lectures, our outcome variable can be represented as a column vector $\\vec{\\textbf{y}}$ with some magnitude and direction.\n",
"\n",
"The dimensions of this column vector are $\\textbf{y}_{n\\times 1}$, where $n$ corresponds with the number of observations. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Using a linear combination of $\\textbf{X}_{n\\times p}$ to approximate $\\textbf{y}_{n\\times 1}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**_Our goal_** is to find a linear combination of the vectors in matrix $\\textbf{X}_{n\\times p}$ that yields $\\textbf{y}_{n\\times 1}$.\n",
"\n",
"However, there might not be enough information in $\\textbf{X}_{n\\times p}$ to do this. \n",
"\n",
"Put differently, $\\textbf{y}_{n\\times 1}$ could be thought of as shooting off the plane of $\\textbf{X}_{n\\times p}$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given that we can't find a linear combination of $\\textbf{X}_{n\\times p}$ to get to $\\textbf{y}_{n\\times 1}$. We can approximate the relationship by taking the **projection** of $\\textbf{y}_{n\\times 1}$ ( $\\hat{\\textbf{y}}$) on $C(\\textbf{X}_{n\\times p})$.\n",
"\n",
"This projection is going to be some **linear combination** of the column vectors of $\\textbf{X}$. Recall that we can get anywhere in the span of $\\textbf{X}$ by arbitrarily scaling and adding the linearly independent vectors.\n",
"\n",
"$$ c_1\\vec{x}_1 + c_2\\vec{x}_2 + \\dots + c_p\\vec{x}_p $$\n",
"\n",
"We want to find the constants that generate the projection of $\\hat{\\textbf{y}}$. We'll call this vector $\\beta$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Finding $\\beta$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Recall that to find the projection of $\\textbf{y}_{n\\times 1}$ onto $C(\\textbf{X}_{n\\times p})$, we need to shoot an **orthogonal vector** from our plane to the vector $\\textbf{y}_{n\\times 1}$. We'll call this vector $\\epsilon$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Another way to think about $\\epsilon$ is that it is the distance between vectors $\\textbf{y}$ (the outcome) and $\\hat{\\textbf{y}}$ (our prediction).\n",
"\n",
" \n",
" \n",
"\n",
"\\begin{equation}\n",
"e_{n \\times 1} = \\textbf{y}_{n \\times 1} - \\hat{\\textbf{y}}_{n \\times 1}\n",
"\\end{equation}\n",
"\n",
" \n",
" \n",
"\n",
"Which we can re-express as\n",
"\n",
" \n",
" \n",
"\n",
"\\begin{equation}\n",
"e_{n \\times 1} = \\textbf{y} - \\textbf{X}\\beta\n",
"\\end{equation}\n",
"\n",
" \n",
" \n",
"\n",
"Recall that two vectors are **perpendicular (orthogonal) to one another when their dot products equal 0**. \n",
"\n",
"> Do you recall why?\n",
"\n",
"We want $\\epsilon$ epsilon to be orthogonal to the span of $\\textbf{X}_{n\\times p}$. Note that we must transpose $\\textbf{X}_{n\\times p}$ to multiply it by $\\epsilon_{n\\times 1}$\n",
"\n",
" \n",
" \n",
"\n",
"\\begin{equation}\n",
"\\textbf{X}_{p \\times n}^T \\cdot e_{n \\times 1} = 0\n",
"\\end{equation}\n",
"\n",
" \n",
" \n",
"\n",
"Now, let's rearrange so that we can isolate our vector $\\beta$, which will yield the equation for $\\beta$.\n",
"\n",
" \n",
" \n",
"\n",
"\\begin{equation}\n",
"\\textbf{X}_{p \\times n}^T \\cdot ( \\textbf{y}_{n \\times 1} - \\textbf{X}_{n \\times p} \\beta_{p \\times 1} ) = 0\n",
"\\end{equation}\n",
"\n",
" \n",
" \n",
"\n",
"\\begin{equation}\n",
"\\textbf{X}_{p \\times n}^T \\textbf{y}_{n \\times 1} - \\textbf{X}_{p \\times n}^T \\textbf{X}_{n \\times p} \\beta_{p \\times 1} = 0\n",
"\\end{equation}\n",
"\n",
" \n",
" \n",
"\n",
"\\begin{equation}\n",
"\\textbf{X}_{p \\times n}^T \\textbf{y}_{n \\times 1} = \\textbf{X}_{p \\times n}^T \\textbf{X}_{n \\times p} \\beta_{p \\times 1}\n",
"\\end{equation}\n",
"\n",
" \n",
" \n",
"\n",
"\\begin{equation}\n",
"(\\textbf{X}_{p \\times n}^T \\textbf{X}_{n \\times p})^{-1} \\textbf{X}_{p \\times n}^T \\textbf{y}_{n \\times 1} = \\beta_{p \\times 1}\n",
"\\end{equation}\n",
"\n",
"\n",
" \n",
" \n",
"\n",
"This provides a formula to find our constants for the linear combination of $\\textbf{X}$ that generates our projection $\\hat{\\textbf{y}}$\n",
"\n",
" \n",
" \n",
"\n",
"$$ \\beta = (\\textbf{X}^T \\textbf{X})^{-1} \\textbf{X}^T \\textbf{y}$$\n",
"\n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.98091425, 2.98340745])"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"B = la.inv(X.T.dot(X)).dot(X.T.dot(y))\n",
"B"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generating a prediction ( $\\hat{\\textbf{y}}$ )\n",
"\n",
"As we've seen earlier, we can now generate a projection of $\\textbf{y}$ by dotting the vector $\\beta$ by $\\textbf{X}$.\n",
"\n",
"$$ \\hat{\\textbf{y}}_{n \\times p} = \\textbf{X}_{n \\times p} \\cdot \\beta_{p \\times 1}$$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"y_hat = X.dot(B)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's plot this line through the data. As we see, we produced a linear fit. In fact, we produced the best linear fit. One that minimizes the sum of the squared error."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3YAAAEyCAYAAAC2+0LeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XuUXFWZ9/Hfk6aBDqjtQBTTSRNUjEETibSobxQhIIGYBS3oEDJLh3Fm4g0dFaMJ+DqOjhLNjDdEHUad0VEIvBAahgDhEh0xA0qH3IQmCCRAOlzCJeGSBpPO8/7R3Ul1d9WpU1Wnzq2+n7V6ddU5u87ZXZRYP/bezzZ3FwAAAAAgu8Yk3QEAAAAAQG0IdgAAAACQcQQ7AAAAAMg4gh0AAAAAZBzBDgAAAAAyjmAHAAAAABlHsAMAAACAjCPYAQAAAEDGEewAAAAAIOP2S7oDQQ499FCfNGlS0t0AAAAAgESsXr36SXcfV65dqoPdpEmT1N3dnXQ3AAAAACARZvZQmHZMxQQAAACAjCPYAQAAAEDGEewAAAAAIOMIdgAAAACQcQQ7AAAAAMg4gh0AAAAAZBzBDgAAAAAyjmAHAAAAABmX6g3KAQAAAKBeutb0asmKjdq6vU/jW1u0YNZkdU5vS7pbVSHYAQAAAGg4XWt6tWjZBvXt6pck9W7v06JlGyQpk+GOqZgAAAAAGs6SFRv3hrohfbv6tWTFxoR6VBuCHQAAAICGs3V7X0XH045gBwAAAKDhjG9tqeh42hHsAAAAADScBbMmq6W5adixluYmLZg1OaEe1YbiKQAAAAAazlCBFKpiAgAAAECGdU5vy2yQG4mpmAAAAACQcQQ7AAAAAI3nid9J179FutSkP/0o6d7UjKmYAAAAAPKv/yVp43eltQtHn9uzK/7+RIxgBwAAACC1utb0Vl/g5Nk/SXd9Rtp6/ehzLeOlt/1QajtNMou20wkg2AEAAABIpa41vVq0bIP6dvVLknq392nRsg2SVDzcuUubL5W6PyHtenb0+fazpOlLpIMm1rPbiSDYAQAAAEilJSs27g11Q/p29WvJio37gt2LT0rrL5Duv6T4RTp+IL3+o9KYfEeffP91AAAAADJr6/a+oscn7vqDtPxT0o57Rp889P9IHRdJf/HWOvcuXQh2AAAAAFJpfGuLerf36eVjntcXX/Of+qtDbtx3ckdBw6MWSm86X2p+Wex9TIuag52ZTZZ0ecGh10r6srt/t6DN8ZKukbRp8NAyd/9qrfcGAAAAkFOPdGlV+/ul9tGndu7XprEz/k1qe1/8/UqpmoOdu2+UdLQkmVmTpF5JVxdpepu7z6n1fgAAAAByaM9uaWlzYJPTe5fqb056V/iqmA0k6qmYJ0p6wN0fivi6AAAAQCbUVJ6/0Ty2Ulp5Yunzr5wunXyH1LS/pIEpgCgu6mA3V9JlJc6908zWSdoq6fPufnexRmY2X9J8SWpvLzLuCgAAAKRUxeX5G9GKt0tP/aH0+SmfH9iSABUxd4/mQmb7ayC0vcndHx9x7uWS9rj782Y2W9L33P3Ictfs6Ojw7u7uSPoHAAAA1NuMxSvVW6SSY1tri1YtnJlAj1Jg5xapK3jfuJte97/6p1/3McpZhJmtdveOcu2iHLE7VdJdI0OdJLn7swWPrzezH5rZoe7+ZIT3BwAAABJVqjx/qeO5tf4fpT8G1Ep8xVHS+wYm8DHKGY0og93ZKjEN08wOk/S4u7uZHStpjKSnIrw3AAAAkLih8vzFjuda/4vS5WX+xuNvkMafMupwqE3IUVYkwc7MDpL0XkkfLTj2MUly9x9L+oCkj5vZbkl9kuZ6VHNAAQAAgJRYMGvysNEnSWppbtKCWZMT7FWdPPAf0u8/EtzmrJf2Fj4phVHOaEQS7Nz9BUmHjDj244LHP5D0gyjuBQAAAKTV0AhTbqtiXmqBp2994T36uwcW7Pu7y4Q6qYFHOSMWdVVMAAAAINXqvR1B5/S2/AS57X+Urp8a3GbOfep6YGzV6+QaapSzjgh2AAAAaBgU6ghh5SzpsZuC28wbvqpqyYqVVa+Ty/0oZ0wIdgAAAGgYFOooYvcL0hUHB7eZsVQ6/KySp2tdJ5erUc6EEOwAAADQMCjUMei+i6Xuc4PbzP2zNKY51OVYJ5c8gh0AAAAaRkMHkDKFT/Tav5He8bOqLs06ueQR7AAAANAwGiqAPH2XdOMxwW1O3ywddHjNt2KdXPIIdgAAAGgYuQ8gK94pPXVH6fNjmgemWNYB6+SSRbADAABAQ8lVAPnzDunK1uA2x10jTTgtnv4gMQQ7AAAAIEvuWSKt/UJwm7m7pTFN8fQHqUCwAwAAANLMXbpsTHCbN3xK6vh+PP2ps3pvIJ9XBDsAAAAgbbb9r3TzjOA2nb3S2PHx9CcmbCBfPYIdAAAAkAbltiM44FDpzG3x9CUhbCBfPYIdAAAAkISdW6SuicFtTlghvebkePqTAmwgXz2CHQAAQM6wRik5Zd/7m4+Ttt0WfJGz+yUrs6Yupxp6A/kaEewAAABypBHXKKUlyBZ7789ftk6dPROCX3jIsdKs38fQw/RrqA3kI0awAwAAyJFGW6OUpiA79N7/7aFX6/+O/2lw49M3SwcdHku/siT3G8jXEcEOAAAgRxptjVJqguylplXtktoD2szzwEukZeQxabnaQD5GBDsAAIAcabQ1SokF2WfWSTccHdhk4ZZztfTpU9TW2qJVC2cGtk3TyCOyiWAHAACQI6XWKJ3wxnGasXhl7kaDYg2y5bYjkHRUz3Lt3LVvZC7s+rDUjDwiswh2AAAAOVJsjdIJbxynq1b35nI0qK7FNvbskpbuH9ym6UDprH3B8htVTqdstCm0iB7BDgAAIGdGrlGasXhlbkeDIi+20f0p6b4fBLc57UHp4CNK9qeaezfaFFpEj2AHAACQc3kfDaq52EaIKZblCp/UijL/qBXBDgAAIOcYDRrh8d9It54Q3OaYi6TJ58bSHYky/6gdwQ4AACDn6jEalLnS/GFG5c7eI1mIdnVCmX/UgmAHAACQc1GPBlVbmj/WMLj7BemKg4PbjG2XOh+qz/2BmBHsAAAAGkCUo0HVlOaPZZ+2G6ZLz6wNbnPaJungSdHcD0gRgh0AAAAqUk0xlrrt05aCwidAGkQW7Mxss6TnJPVL2u3uHSPOm6TvSZotaaekc9z9rqjuDwAAgHhUU4wlssqcm34p3f6h4DbTl0hTPl/ZdYGMi3rE7gR3f7LEuVMlHTn483ZJPxr8DQAAgAypphhLTZU5IxyVy1zRFyCkOKdini7pF+7uku4ws1Yze427PxpjHwAAAEIjBBRXTTGWE944Tr+642EVxq+SYfClp6WrDinfkQqnWMayzq8KfM4QhSiDnUu6ycxc0r+5+yUjzrdJeqTg+ZbBY8OCnZnNlzRfktrb2yPsHgAAQHhpDQFpUUkxlq41vbpqde+wUGeSzjym4BqX7Sd5f7GXF9y0Vxo7vqr+SnVc51cDPmeIypgIr/Uud3+rBqZcftLMjqvmIu5+ibt3uHvHuHHjIuweAABAeEEhAJUp9l66pH/eOW1gmuWlVjrUzfN9PzWEOinCdX4R4nOGqEQ2YufuvYO/nzCzqyUdK+m3BU16JU0seD5h8BgAAEDqpDEEJKmW6YJD79nHx/0/ffE1Pw9ufMz3pMmfrrW7RbWObdYzO3cVPZ4UPmeISiTBzswOkjTG3Z8bfHyypK+OaHatpHPNbKkGiqbsYH0dAABIq5qKfdRBkuuwapoueKlp07QyN4hpOwIvcZtSx+OQts8ZsiuqqZivlvQ7M1sn6Q+Slrv7jWb2MTP72GCb6yU9KOl+Sf8u6RMR3RsAACByC2ZNVktz07Bj5So/Rq1rTa9mLF6pSQuX67OXr1Xv9j659gWrrjXxTH6qaLrgCw/vm14ZUM1ySs8N6pqyJdY95nb0jR6tCzoehzR8zpAPkYzYufuDkt5S5PiPCx67pE9GcT8AAIB6q6byY5RGjpKNjD9xFv0oNS2wd3ufjli4XJumzSl7jeuPXKuv3/rU3vfywjPir/yYxtGxpD9nyI84tzsAAADIlEoqP0at2CjZSLWsw6pkamexQLQ5RJgrHI2bLWn226rubiSq2X8vDkl+zpAfBDsAAIAUChPaqh1pqnTN3IJZk3X46vdq+tiewOt+95lP6DOfvLiqPsWB0THkGcEOAAAgRYZG0sqtPKtlpCn0fm6Da+Q6JWls8WtNWn/d3scm6TNV9Sg+jI4hrwh2AAAAKTFyJG0k08Bau7YaR5pKjQaOe2mtdOmJZV8/4+FbU7dWDWh0BDsAAICUCFpXV2uYK1S4Zi7UWrkzn5QOOGTv0wVFAmga1qoBjYxgBwAAkBKlRtJM0qqFM6O5ibtWtZ8otZdpF7ANAWvVgPQh2AEAAKRE3crxB+wnN+Rbj35Y1+z669ABkrVqQLoQ7AAAAFIi0nL8IcJcYeETSTJVv30CgGQR7AAAQCIq2UetUdQ0xfGhy6VVc8u3m+easXglxU+AnCHYAQCA2FW6j1ojqWiKY4hROX3wWan5ZcMO1TIySCAH0olgBwAAYhd6HzUMt2e3tLS5bLOuKVsC38dqRwYJ5EB6EewAAEDsSlV/LHW8oYUYlfveEx/Wdx77y73PW3rKh61qip8QyIH0ItgBAIDY1a36Y16EmWI5uB1BsfVy9QpbBHIgvcYk3QEAANB4FsyarJbmpmHHGnqD655vD4S5oZ9S5vm+n0Fxhq1SwZtADiSPETsAABA7NrhWuFG5s16Umg4IbBLn6Gek2zEAiBTBDgAAJKLhNrje3SddMbZ8u4LRuDDiDFsEciC9CHYAAAD1EmZU7m0/lI78eNW3iDtsNVwgBzKCYAcAABClCgqfRIWwBYDiKQAAALXo/lSowidTem5Q15QtkYc6AJAYsQMAADnVtaa3ftMTQ4zKHbH+Wvmw/4ber/OuWCeJzbwBRI9gBwAAcqdrTe+wgiK92/u0aFn5TbtLeukp6apDy7cbHI07YuFyFRuX63evrR8AUAJTMQEAQO4sWbFxWJVIad+m3aEVTq8sEeq+/OSXdMT66zTj4VsHplkOCtpqoOJ+AEAIBDsAAJA7VW/aHXKT8K4pWzSl5wb9Yus75No3Iti1pldS8Q3YK+oHAFSIqZgAACB3Qm/avfJk6bGby19wRMGToBHBwgqV512xTv0+elJmPTYPB9DYCHYAACB3AjftjmA7gjAjgkPhLq7NwwE0NoIdAAApV9fqjjlVuGm3vfCQfjflIwMnegJeVME2BGFHBOPePBxA4zIvMj0gLTo6Ory7uzvpbgAAkJiR1R2lgRGfC8+Ymkg4KBcyUxNCw4zKzbxFOuzEqi6ftn8uAPLLzFa7e0e5djWP2JnZREm/kPRqSS7pEnf/3og2x0u6RtKmwUPL3P2rtd4bAIC8K7eWK07lthCIfIuBSkUwxTKsSkbiUhN2AeRaFFMxd0s6z93vMrOXSVptZje7+z0j2t3m7nMiuB8AAA2j6uqOdVAuZMYeQq9uk/q2lm8XUZgbqbBISimJh10ADaPmYOfuj0p6dPDxc2bWI6lN0shgBwAAKhS6umMMyoXMWEJojKNyUUjTiCuAfIt0HzszmyRpuqTfFzn9TjNbZ2Y3mNmbAq4x38y6zax727ZtUXYPAIDMKbYfWlJVFUuFyaHj5c5X5em7Qu8tt/cnRdI04gog3yILdmZ2sKSrJH3G3Z8dcfouSYe7+1skXSSpq9R13P0Sd+9w945x48ZF1T0AADKpc3qbLjxjqtpaW2SS2lpbEivQUS5kRhZCC4PcjccUb3NKd2rDXKG6hF0AKCKS7Q7MrFkDoe5X7r5s5PnCoOfu15vZD83sUHd/Mor7AwCQZ2HWcsXVD6l0wZCaSvtnbIplWIH76QFAhKKoimmSfiqpx92/XaLNYZIed3c3s2M1MFL4VK33BgAA8SoWMquq+hgmyDWNlc56oYbeJo997ADEJYoRuxmSPiRpg5mtHTx2vqR2SXL3H0v6gKSPm9luSX2S5nqaN9ADAAChVFT1MaejcuXEPeLK9gpAY2KDcgAAItZIX6xnLF5ZtGpnW2uLVp3j0sqTyl8kh2EuKWycDuRPbBuUAwCAfdK8b1k9AufI6o6bpxVsWbuyxItOe0A6+LU13RfFsb0C0LgIdgAARCitX6zrFTjHtx6oVe3lR+W6pmzZFyof3qwFsw4gaNQB2ysAjYtgBwBAhNL6xbpU4DzvinWSKgx3BWvlVrUXb/JUS4cOef+dktI9ipk3adrQHkC8It2gHACARpfWfctKBct+dy1atkFda3qDLxBik/Aj1l+nGQ/fqq4pW/aGOil4FBPRStOG9gDixYgdAAARSuu+ZaVGcqQSU0U3/VK6/UPlL1xQ+GTTvOJN0jqKmUdsrwA0LoIdAAARSusX62KBs9DW7X2htiO44cg1+udbn973t63pLfu3MT0wXmnZ0B5AvAh2AABELI1frIf6c94V69Q/uNWRaY82TTut/IsHR+WqXSuX1lFMAMgTgh0AAA2ic3qbOnsmlG848Uzp3VeOOlxtxc9qRzEbaT9AAKgVwQ4AgLwLMcXyiPXXlQ1PtayVq3QUk0qaAFAZqmICAJA3Pd8OVcWya8oWTem5QZPWXyfXvvBUqkJmnBU/qaQJAJVhxA4AgDwIMSqnDz4nNR+89+mSxSsrmloZ51o5KmkCQGUIdgCATGC91Qj9L0mXH1i+XcF2BCNVGp7irPhJJU0AqAzBDgCQemlcb5VI0AwzKnfUQunoC0NdrprwFFfFTyppAkBlCHYAgNSrthpjvcQaNMOEuYBRuSBpDk9p3Q8QANKKYAcASL20rbeqa9Bc/Tlp43fKt6syzBVKe3hK436AWcVUZiD/CHYAgNRLer3VyC/Fxfoi1RA0w4zKzf2zNKa5uusHIDzlXxqnMgOIHtsdAABSb8GsyWppbhp2LK4pg0Nfinu39+3dEqBUDAsdNHc9G2o7As3zfT91CHVoDGwdATQGRuwAAKkQNFUsySmDxb4UuyQb/D2kbNAMMSp3/pZP6tKnT1VLc5MuPGOqOqvqMTBc2qYyA6gPgh0AIHFhpoolNWWw1Jdfl9TW2hIcNEOEuRkP3zpqameShWGQP0lPZQYQD4IdACBxaat6WajUl+K21hatWjhz+MFVZ0sPLS1/0YLCJ1sXLi/ahNEURCXN1U8BRIdgBwBIXJqnipX9Uhym8MnZeyQr3q5UcBxjpq41vYkHW2Rf2qufAogGwQ4AGlSayp+nearYyC/F0w55Qde0nSX1aOCnlJDbERQLjpLU707lQkSG6qdA/hHsAKABpa38edqninX2TFBnu6T2gEbvulJqP7Pyaw++3+ddsU79PjwMpmU6KgAg/Qh2ANCA0ramLZVTxcJMsYxgk3Bp4O//7OVri55Lw3RUAED6EewAoAGlcU1brVPFap5aesN06Zni4WqYiMLcSGmejgoASD+CHQA0oCyFiDCBreqppTGOypWT9umoUUjTuk4AyBuCHQA0oKyEiLCBLfTU0h090vKjyt84pjBXqNh01BPeOE5LVmzUZy9fm/kglLZ1nQCQN5EEOzM7RdL3JDVJ+om7Lx5x/gBJv5B0jKSnJJ3l7pujuDcAoHKpXNOm0SM6O/+8O1RgC5xaGmZU7qT/kV51XCR9ruV9LJyOmrcglLZ1nQCQNzUHOzNrknSxpPdK2iLpTjO71t3vKWj2t5KecffXm9lcSd+UdFat9wYAVC9t5c+LBZlSRga5kVNLN0+bU/6GEYzK1TN85S0IpXFdJwDkSRQjdsdKut/dH5QkM1sq6XRJhcHudElfGXx8paQfmJm5e/xzXQAgpRp9/VGxIFPKyLWAvz18tpraXyr/woinWNYzfOUtCGVpXScAZFEUwa5N0iMFz7dIenupNu6+28x2SDpE0pMjL2Zm8yXNl6T29qANgwAgP/I27a4aYQPL3rWABVMsm0o1rvNauXqGr7wFoays6wSArBqTdAdGcvdL3L3D3TvGjRuXdHcAIBZBIz+NolRgaW1pVltri6a3bNTmaXPUM+VUdfZMKH2heb7vp85K9TmK8LVg1mS1NA+PrFkOQp3T23ThGVPV1toik9TW2qILz5jaMP/hAgDqLYoRu15JEwueTxg8VqzNFjPbT9IrNFBEBQCg/E27q0axEZ1Qa+Vmb5Ba31zHnpVWz1GotBa4qUXa1nUCQJ5EEezulHSkmR2hgQA3V9K8EW2ulfTXkm6X9AFJK1lfBwD7pG3aXRLr/YauHzgaN2jS+uskDYSoC6e8Up3T69q1kuodvghCAICwag52g2vmzpW0QgPLHH7m7neb2Vcldbv7tZJ+Kum/zOx+SU9rIPwBAAalaf1R7Ov9CtbKdZZq8xfHaMbab40Kv6UKlcQZTAlfAIA0iGQfO3e/XtL1I459ueDxi5I+GMW9ACCP0jTtLpYy+2H2lhuxRm7rb5YXbTZyuiqFaAAAjSiSYAcAqF1aRn7qst5v6wrpN6eUbxdQ8CTsdNW87f8GAEAYBDsAwDCRrfcLMyrXuUUaGy5shZ2uSiEaAEAjItgBQIMpt/4sbIAadZ2T36DOeyeqrCq3IQg7XTVMMB3Z9xPeOE6/vndb4tNgAQColqW5OGVHR4d3d3cn3Q0AyI2R68+kwcqSI/YTKxf+hq7TM+XU8jedeKb07isj/TuClPsbi50fqdh7AgBAEsxstbt3lG1HsANQShIl71FfMxavLDqa1dbaolULZ4a7SBWFT+IW9Nkt9R6MVNF7AgBAnYQNdkzFBFAUlQXzqar1Zw/+XLrjnLLXnrT+OpmkTYvfV13nIhRUiCbsWjvW5AEAsoRgB6AoKgvmU+jCKCFG5ab9came3XNw8HVSqNR7UKwdAABZMSbpDgBIJyoL5tOCWZPV0tw07FhLc5MWnPz6gTA39FPKPJfmubqmbNGupleMvk4CG6pXqth7MFJW/hYAAIYwYgegqMhK3keMdX+1Kawsuar9xH0n7i3e/kdPfEC/fPHjo97nNG2oXqlifacqJgAg6yieAqCosNUTG71PmRNiimXXlC28zwAApATFUwDUJI0jMqz7C1Z0NPPApdKaz5d/cUEVyyWLV/I+AwCQMQQ7ACUFVRZMAuv+Sisczdw8bc7AwZ6AF5z1otR0QNFTvM8AAGQPwQ5AZqR13V/i+l9UZ88EdU4p067I3nLFRvl4nwEAyB6CHYDMWDBrctG1X41QvXBkABtW+KSEr/TO18+fOq3kvnKl9io885g2XbW6tyHfZwAAsopgByAz4lr3l7bKm0MBrGfKqVJ7cNtJ668b9rwtYJSt1JrFX9+7TReeMTVV7wEAAAhGsAOQKfVe91dqFGvo3rVeu6KwdPc3pHUXqFMKnmY5z/f2Wwo/yha0li5t6ysBAEAwgh0AFKhX5c3QgTHEdgSvXX+N9qhJJu2dZhk0mlkqULKWDgCA/CDYAUCBchUhq52mWSowXnzTGnX2TCj7+pFTLKXRAazYKFtQoGzkNYsAAOQNwQ4ACgSNYtUyTbMwMO7djiDIzJulw06SJH2pa4NMD6uwpmXYABY0Arlq4cy9bVhLBwBAthHsAKBA0ChWLdM0N4UJcyW2I7hqde+wUGeSzjwm3Bq4ciOQrKUDACAfxiTdAQBIk87pbbrwjKlqa22RaaCq5IVnTFXn9LbKNu6+Z8nAermhnyLuf2miuqZsGQh0RUKdVHzEzSX9+t5tof6eUuvlWEcHAEC+MGIHACOUGsUqW2wkROGTGQ/foq3bXww97bGiMFkE6+gAAGgMBDsACGlkSHpl0w6tedNfDZy8NOCFBaNxq0Lea6hIS/FxvPAjbnHt/ReltO0jCABAFph7qa8Nyevo6PDu7u6kuwGgQYQKFCFG5TR7vdQ6taZ+jBxlK9TS3KQzj2nTr+/dlrvwU+xvb2lu2jsdFgCARmNmq929o1w7RuwApELSozSBFS9DbEdQao1cNYqtqxvS1tqiE944Tlet7q3LJupJq9c+ggAA5B3BDkDiatlGICqFgeJTr7pM5x32q4ETPSVeMPEM6d1X1aUvpdbPmaRVC2dqxuKVuQ0/ta4pBACgURHsACQuDaM0q9pPlNrLNIpwVC5IuSIteQ4/ZQvUAACAomra7sDMlpjZvWa23syuNrPWEu02m9kGM1trZiyaAzBMlEGla02vZixeqSMWLteMxSvVtaa3eMOdW8tuRyBJr9uwfN+WBDFZMGuyWpqbhh0rrGSZ5y0Myv3tAACguFpH7G6WtMjdd5vZNyUtkvTFEm1PcPcna7wfgByKapSm7JTOEIVP3tnzH3p017iCIx77tNBylSzzvIVBFqt4AgCQBpFVxTSz90v6gLv/VZFzmyV1VBrsqIoJNIZilRCbm0wH7b+fdvTtCv3lfsbilaMC4uZpc8p3YHA0rmtNr867Yp36i/x7sa21RasWzgzx18Qj6WIzAAAgHklUxfyIpMtLnHNJN5mZS/o3d7+k1EXMbL6k+ZLU3l5uwQuAPBg5StM6tlnPv7hb2/t2SQpfTGXr9j594bD/1CdedWXwDY9aJB39jaL9+Ozla0teO01KbaIOAAAaU9lgZ2a3SDqsyKkL3P2awTYXSNot6VclLvMud+81s1dJutnM7nX33xZrOBj6LpEGRuxC/A0AUqKWUaTCoDJj8Uo9s3PXsPOBxVQGp1humhZwg5Br5CjeAQAAsqhssHP3k4LOm9k5kuZIOtFLzOt0997B30+Y2dWSjpVUNNgByKYotywoNTrWu71PRyxcrrcc8ry62uaWvc6k9dft29w65L3zvH4NAADkV01TMc3sFElfkPQed99Zos1Bksa4+3ODj0+W9NVa7gsgfaLcsqDYqFmotXIfeFpdd+/cO2rYVsXaM4p3AACALKp1jd0PJB2ggemVknSHu3/MzMZL+om7z5b0aklXD57fT9Kl7n5jjfcFkDJRblkwNGrWM+XU8o1HTLHsnP7KmkMY69cAAEDW1BTs3P31JY5vlTR78PGDkt5Sy30ApF8ka9M2XiSt/rQ6JXVOKd5k/uYLdNOz75QkmaRN8yrvKwAAQN5EWRUTQAOrem1aiL3lZjx8KwVNAAAAAhDsAEQi7Nq0G+5WAqBvAAAPG0lEQVRco1P/9Nbgi41/n3T8dXufLiiyz12eC5qwRx0AAKgUwQ5AZEquTbvxbdLT3ZKkkqvmzuqTmg4seV2pMQqaVFJdlAAIAACGEOwARM9dumxM2WaT1g+MyrW1tmhViVA3pFEKmoStLhrl9hIAACD7CHYAovHwldLvPhjYZPZ939c9L7521PFqKmfmVdjqolFuLwEAALKPYAegeiEKnxRuR7Bj8UrpRYqgBAlbXTTK7SUAAED2lZ8rBQBDdvYOhLmhn2I6Lh4Ic0M/BRbMmqyW5qZhx/JcBKUaYd+jUmGYkAwAQGNixA5AsNvPkTb9PLjN3N3SmKbgNqpvEZS8FBIJ+x5Vvb0EAADIJXP38q0S0tHR4d3d3Ul3A8i8ikKP75EuKxPSXvlW6dTV0Xe0SiMLiUgDIefCM6ZmMtyFlZcwCwAASjOz1e7eUa4dI3ZAzoWqnvhIl3Tb+4MvdNqD0sFH1LOrVWvUQiKNUikUAACUR7ADcq5U6OnsmSD1lHnxvPSO6BeikAgAAGh0BDsg54bCzWHNT+qOKecENz6uS5pwev07FbGwlSQBAADyimAH5NmaL2jTtCXBbc7ulyzbBXIpJAIAABodwQ5IsYqLY+zpl5YG/8/68qffq688/rlcFRapZ7VNAACALCDYASkVquiJJG29UfrNqYHXWv6GDfrGLU/sDT0XnpG/0EMhEQAA0MgIdkBKBVZ6vP9N0q4dpV/cOk2avW7v0/dJel/ZIrkAAADIKoIdkFKFFR1HFT7ZVeQFJ98uHfqOuvcLAAAA6UOwA1LqsxOX69Ov/FFwo7P3SGbxdAgAAACpRbAD0mJE4ZNPv3J0k28+/neafPxXWEsGAACAYQh2QJKeuE265bjAJjMfWa5NzziVHgEAAFASwQ6I2y3vkZ74benzbzhX6rho79OVMXQJAAAA2UawA+qt71Hp6vHBbebcJ738yHj6AwAAgNwh2AH1cPc3pHUXlD5/0CTp9E2xdQcAAAD5RrADorBnl7R0/+A2x10jTTgtnv4AAACgoRDsgGo9dqu08qTgNmf1SU0HxtMfRKZrTa+WrNiordv7KFoDAAAygWAHVOKGY6Rn7ip9fsoXpOnfjK8/iFzXml4tWrZBfbv6JUm92/u0aNkGSSLcAQCA1CLYIZNiG1F54RHpmvbgNqdvlg46PPp7IxFLVmzcG+qG9O3q15IVGwl2AAAgtWoKdmb2FUl/L2nb4KHz3f36Iu1OkfQ9SU2SfuLui2u5Lxpb3UdU1n1Juvvrpc+3TpVmr6/9Pkilrdv7KjoOAACQBlGM2H3H3f+l1Ekza5J0saT3Stoi6U4zu9bd74ng3mhAkY+o9L8kXV5mHdzxN0rjZ1V+7RBYz5Uu41tb1FskxI1vbUmgNwAAAOHEMRXzWEn3u/uDkmRmSyWdLolgh6pEMqLSu1z6nznBbc56SWoqU+myRqznSp8FsyYP+2ciSS3NTVowa3KCvQIAAAgWRbA718w+LKlb0nnu/syI822SHil4vkXS20tdzMzmS5ovSe3tZdY2oSFVPaLy35Ol5+4rfX7qP0lTv1xj7yrDeq70GXrfGUUFAABZUjbYmdktkg4rcuoCST+S9DVJPvj7XyV9pJYOufslki6RpI6ODq/lWsin0CMqzz8oXfu64It1bpHGJveFnfVc6dQ5vY0gBwAAMqVssHP3Mht1DTCzf5d0XZFTvZImFjyfMHgMqErgiMrqz0kbv1P6xYe8Q5p1e0w9LY/1XAAAAIhCrVUxX+Pujw4+fb+kPxZpdqekI83sCA0EurmS5tVyX2DviMruPumKsVKPBn6KOfE30qvfE2PvwmM9FwAAAKJQ6xq7b5nZ0RqYirlZ0kclyczGa2Bbg9nuvtvMzpW0QgPbHfzM3e+u8b5oZE/eId30zuA2c3dJY9K/TSPruQAAABAFc0/vMraOjg7v7u5OuhtImrv0h/nSAz8p3ebob0lHLYivTwAAAEAMzGy1u3eUa5f+IQ00pp29UteE0uf3O0jq7JX2f0V8fQIAAABSimCH9Lj/EukPHy19ftrXpDd/Kb7+AAAAABlBsENy+l+Ulr9Zev6B0m3mbJRe/ob4+gQAAABkEMEO8XriNumW40qfHz9Hes81ko2Jr08AAABAxhHsUF/u0u0fljb/snSb42+Qxp8SX58AAACAnCHYIXovPCRdM6n0+QMOkU7bJDW/LLYuAQAAAHlGsEM0Nl4krf506fNsRwAAAADUDcEOkqSuNb2VbZK9+wXpv4+U+h4t3ea0B6SDXxt9ZwEAAAAMQ7CDutb0atGyDerb1S9J6t3ep0XLNkjS8HD32Epp5YmlL9T+l9KMpZJZPbsLAAAAYASCHbRkxca9oW5I365+LVlxrzqf/4z0yJWlXzzzVumwmXXuIQAAAIAgBDto6/a+vY8n7v+Ybnvj3+07+ciIxmMnDOwtt9/YeDoHAAAAoCyCHTTrVQ9oQesSve7ALcUbHPN9afKn4u0UAAAAgNAIdo2o/yVp43eltQslST8+bHSTmX/6uT592onBBVQAAAAApALBrlE8d7+0+jPS1uWjzx14mH5/yNf0uduP0NbtLw5UxTytTFVMAAAAAKlBsMsrd2nzpVL3J6Rdz44+P/ED0lv/RTrocEnS2yWtek+8XQQAAAAQDYJdnrz0lLTuS9L9Py5+/pjvS0d+XBrDP3YAAAAgT/iGn3WP/8/AqNyOe0afO+QdUsdF0iEd8fcLAAAAQGwIdlnT/6LU86/S+i8VPz9lgfTmL0nNL4+3XwAAAAASQ7DLgh33Sqv/QXrsptHnxk6QOn4otc2RzOLvGwAAAIDEEezSyPdIm/5LuvMTUv/O0ecPnysd/S3poInx9w0AAABA6hDs0uLFJ6R150sP/LT4+bf9UHrd31P4BAAAAMAopIQkPXbrwKjcc/eNPjfuXQNVLP9ievz9AgAAAJApBLs47e6Ter4lbfhK8fNvOl86apHUfHCs3QIAAACQbQS7ett+t7T609LjK0efO2iS1HGx1DY79m4BAAAAyA+CXdR8j/TgfwxMsdzz59HnJ31IOnqxNHZ8/H0DAAAAkEsEuyj0PSatXSht+vnoc7bfQOGT135EGtMUf98AAAAA5B7BrlpbV0jdn5Cef3D0uVcdL3V8X2qdGnu3AAAAADSemoKdmV0uafLg01ZJ29396CLtNkt6TlK/pN3u3lHLfRP12K3SypNGH3/zl6WjvijtNzb+PgEAAABoaDUFO3c/a+ixmf2rpB0BzU9w9ydruV8q9L808Pvg10tvu1h6zcnJ9gcAAABAw4tkKqaZmaS/lDQziuulWttsaZ4n3QsAAAAA2GtMRNd5t6TH3f1PJc67pJvMbLWZzQ+6kJnNN7NuM+vetm1bRN0DAAAAgPwqO2JnZrdIOqzIqQvc/ZrBx2dLuizgMu9y914ze5Wkm83sXnf/bbGG7n6JpEskqaOjg6ExAAAAACijbLBz9yKVQvYxs/0knSHpmIBr9A7+fsLMrpZ0rKSiwQ4AAAAAUJkopmKeJOled99S7KSZHWRmLxt6LOlkSX+M4L4AAAAAAEUT7OZqxDRMMxtvZtcPPn21pN+Z2TpJf5C03N1vjOC+AAAAAABFUBXT3c8pcmyrpNmDjx+U9JZa7wMAAAAAKC6qqpgAAAAAgIQQ7AAAAAAg4wh2AAAAAJBxBDsAAAAAyDiCHQAAAABkXM1VMRtJ15peLVmxUVu392l8a4sWzJqszultSXcLAAAAQIMj2IXUtaZXi5ZtUN+ufklS7/Y+LVq2QZIIdwAAAAASxVTMkJas2Lg31A3p29WvJSs2JtQjAAAAABhAsAtp6/a+io4DAAAAQFwIdiGNb22p6DgAAAAAxIVgF9KCWZPV0tw07FhLc5MWzJqcUI8AAAAAYADFU0IaKpBCVUwAAAAAaUOwq0Dn9DaCHAAAAIDUYSomAAAAAGQcwQ4AAAAAMo5gBwAAAAAZR7ADAAAAgIwj2AEAAABAxhHsAAAAACDjCHYAAAAAkHEEOwAAAADIOHP3pPtQkpltk/RQmWaHSnoyhu4g3/gcIQp8jhAVPkuIAp8jRIXPUrIOd/dx5RqlOtiFYWbd7t6RdD+QbXyOEAU+R4gKnyVEgc8RosJnKRuYigkAAAAAGUewAwAAAICMy0OwuyTpDiAX+BwhCnyOEBU+S4gCnyNEhc9SBmR+jR0AAAAANLo8jNgBAAAAQEMj2AEAAABAxmU+2JnZ18xsvZmtNbObzGx80n1CNpnZEjO7d/DzdLWZtSbdJ2SPmX3QzO42sz1mRmloVMTMTjGzjWZ2v5ktTLo/yCYz+5mZPWFmf0y6L8guM5toZr82s3sG/3/tH5LuE4JlPthJWuLu09z9aEnXSfpy0h1CZt0s6c3uPk3SfZIWJdwfZNMfJZ0h6bdJdwTZYmZNki6WdKqkoySdbWZHJdsrZNR/Sjol6U4g83ZLOs/dj5L0Dkmf5N9J6Zb5YOfuzxY8PUgS1WBQFXe/yd13Dz69Q9KEJPuDbHL3HnffmHQ/kEnHSrrf3R909z9LWirp9IT7hAxy999KejrpfiDb3P1Rd79r8PFzknoktSXbKwTZL+kORMHMvi7pw5J2SDoh4e4gHz4i6fKkOwGgobRJeqTg+RZJb0+oLwCwl5lNkjRd0u+T7QmCZCLYmdktkg4rcuoCd7/G3S+QdIGZLZJ0rqR/jLWDyIxyn6XBNhdoYPrBr+LsG7IjzOcIAIA8MLODJV0l6TMjZsohZTIR7Nz9pJBNfyXpehHsUEK5z5KZnSNpjqQTnU0eUUIF/04CKtEraWLB8wmDxwAgEWbWrIFQ9yt3X5Z0fxAs82vszOzIgqenS7o3qb4g28zsFElfkHSau+9Muj8AGs6dko40syPMbH9JcyVdm3CfADQoMzNJP5XU4+7fTro/KM+yPihhZldJmixpj6SHJH3M3fkvnKiYmd0v6QBJTw0eusPdP5Zgl5BBZvZ+SRdJGidpu6S17j4r2V4hK8xstqTvSmqS9DN3/3rCXUIGmdllko6XdKikxyX9o7v/NNFOIXPM7F2SbpO0QQPfsyXpfHe/PrleIUjmgx0AAAAANLrMT8UEAAAAgEZHsAMAAACAjCPYAQAAAEDGEewAAAAAIOMIdgAAAACQcQQ7AAAAAMg4gh0AAAAAZNz/B0Uj+CtnDdY6AAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(15,5))\n",
"plt.scatter(X[:,1],y)\n",
"plt.plot(X[:,1],y_hat,color=\"orange\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## What is $e$? \n",
"\n",
"It's our error! We can think of it as the distance between our approximation of $\\textbf{y}$ and the actual observed $\\textbf{y}$. Can you see why $\\epsilon$ minimizes the squared distance? \n",
"\n",
"Recall that we can get a vectors squared length by dotting the vector by itself.\n",
"\n",
" \n",
" \n",
"\n",
"$$ e_{1 \\times n}^T \\cdot e_{n \\times 1} = ||e||^2$$\n",
"\n",
" \n",
" \n",
"\n",
"When $\\epsilon$ is orthogonal to $\\textbf{X}$ we're minimizing the length of $\\epsilon$. Any other (non-orthogonal) vector would have a greater length than our current vector $\\epsilon$. Consider the many grey lines below. Each one has an angle greater or less than 90 degrees, and each has a length longer than the orange line."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"e = y - y_hat"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that when we plot the values of $e$ we can see that they are normally distributed with mean of 0. This is a key assumption of the linear model. "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAAEyCAYAAACmpOSfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAADdBJREFUeJzt3X+IZXd5x/HP02xsiwpGMqRpjF2xwRIa3ZQhtSjFn22U0mhoS/OHTWlgFQwoCMU1UC2txGLVQinSlQTzR+oPmpWIBGsaAqlQ007sNrvJak1FacKandRaI4WWTZ7+Mdd0dt3pTGbunbvfva8XDHPPuWfmPIcN+86599yz1d0BgLPdj817AADYCsECYAiCBcAQBAuAIQgWAEMQLACGIFgADEGwABiCYAEwhD27ubMLL7yw9+7du5u7BOAs98ADDzzR3Uubbberwdq7d29WVlZ2c5cAnOWq6ttb2c5LggAMQbAAGIJgATAEwQJgCIIFwBAEC4AhCBYAQxAsAIYgWAAMQbAAGIJgATCEXb2XILDmwKEj8x7hFDdfe8W8R4BNOcMCYAiCBcAQBAuAIQgWAEMQLACGIFgADEGwABiCYAEwBMECYAiCBcAQBAuAIQgWAEPYNFhV9RNV9Q9V9c9V9VBV/eFk/Uuq6v6qeqSqPlNVz5n9uAAsqq2cYf13ktd19yuS7EtydVW9MsmfJPlYd/9skv9IcsPsxgRg0W0arF7zg8ni+ZOvTvK6JH89WX9bkrfMZEIAyBbfw6qq86rqcJITSe5O8q9JvtfdJyebPJrkkg1+dn9VrVTVyurq6jRmBmABbSlY3f1Ud+9L8qIkVyX5ua3uoLsPdvdydy8vLS1tc0wAFt2zukqwu7+X5N4kv5TkBVX1w3+x+EVJHpvybADwjK1cJbhUVS+YPP7JJG9Mcixr4fqNyWbXJ7lzVkMCwJ7NN8nFSW6rqvOyFrjPdvcXqurhJJ+uqj9O8k9JbpnhnAAsuE2D1d0PJrnyDOu/mbX3swBg5tzpAoAhCBYAQxAsAIYgWAAMQbAAGIJgATAEwQJgCIIFwBAEC4AhCBYAQxAsAIYgWAAMQbAAGIJgATAEwQJgCIIFwBAEC4AhCBYAQxAsAIYgWAAMQbAAGIJgATAEwQJgCIIFwBAEC4AhCBYAQxAsAIYgWAAMQbAAGIJgATAEwQJgCJsGq6ourap7q+rhqnqoqt41Wf+Bqnqsqg5Pvt48+3EBWFR7trDNySTv6e6vVtXzkzxQVXdPnvtYd//p7MYDgDWbBqu7jyc5Pnn8ZFUdS3LJrAcDgPWe1XtYVbU3yZVJ7p+surGqHqyqW6vqginPBgDP2HKwqup5Se5I8u7u/n6Sjyd5aZJ9WTsD+8gGP7e/qlaqamV1dXUKIwOwiLYUrKo6P2uxur27DyVJdz/e3U9199NJPpHkqjP9bHcf7O7l7l5eWlqa1twALJitXCVYSW5Jcqy7P7pu/cXrNntrkqPTHw8A1mzlKsFXJXlbkiNVdXiy7n1JrquqfUk6ybeSvH0mEwJAtnaV4JeT1Bmeumv64wDAmbnTBQBDECwAhiBYAAxBsAAYgmABMATBAmAIggXAEAQLgCEIFgBDECwAhiBYAAxBsAAYgmABMATBAmAIggXAEAQLgCEIFgBDECwAhiBYAAxhz7wHAObvwKEj8x7hGTdfe8W8R+As5QwLgCEIFgBDECwAhiBYAAxBsAAYgmABMATBAmAIggXAEAQLgCEIFgBDECwAhiBYAAxh02BV1aVVdW9VPVxVD1XVuybrX1hVd1fVNybfL5j9uAAsqq2cYZ1M8p7uvjzJK5O8s6ouT/LeJPd092VJ7pksA8BMbBqs7j7e3V+dPH4yybEklyS5Jsltk81uS/KWWQ0JAM/qPayq2pvkyiT3J7mou49PnvpOkos2+Jn9VbVSVSurq6s7GBWARbblYFXV85LckeTd3f399c91dyfpM/1cdx/s7uXuXl5aWtrRsAAsri0Fq6rOz1qsbu/uQ5PVj1fVxZPnL05yYjYjAsDWrhKsJLckOdbdH1331OeTXD95fH2SO6c/HgCs2bOFbV6V5G1JjlTV4cm69yX5UJLPVtUNSb6d5LdmMyIAbCFY3f3lJLXB06+f7jgAcGbudAHAEAQLgCEIFgBDECwAhiBYAAxBsAAYgmABMATBAmAIggXAEAQLgCEIFgBDECwAhiBYAAxBsAAYgmABMATBAmAIggXAEAQLgCHsmfcAsFsOHDoy7xGAHXCGBcAQBAuAIQgWAEMQLACGIFgADEGwABiCYAEwBMECYAiCBcAQBAuAIQgWAEMQLACGsGmwqurWqjpRVUfXrftAVT1WVYcnX2+e7ZgALLqtnGF9MsnVZ1j/se7eN/m6a7pjAcCpNg1Wd9+X5Lu7MAsAbGgn72HdWFUPTl4yvGCjjapqf1WtVNXK6urqDnYHwCLbbrA+nuSlSfYlOZ7kIxtt2N0Hu3u5u5eXlpa2uTsAFt22gtXdj3f3U939dJJPJLlqumMBwKm2Fayqunjd4luTHN1oWwCYhj2bbVBVn0rymiQXVtWjSd6f5DVVtS9JJ/lWkrfPcEYA2DxY3X3dGVbfMoNZAGBD7nQBwBAEC4AhCBYAQxAsAIYgWAAMQbAAGIJgATAEwQJgCIIFwBAEC4AhCBYAQxAsAIYgWAAMQbAAGIJgATAEwQJgCIIFwBAEC4AhCBYAQxAsAIYgWAAMQbAAGIJgATAEwQJgCIIFwBAEC4AhCBYAQxAsAIYgWAAMQbAAGIJgATAEwQJgCJsGq6puraoTVXV03boXVtXdVfWNyfcLZjsmAItuK2dYn0xy9Wnr3pvknu6+LMk9k2UAmJlNg9Xd9yX57mmrr0ly2+TxbUneMuW5AOAU230P66LuPj55/J0kF220YVXtr6qVqlpZXV3d5u4AWHQ7vuiiuztJ/z/PH+zu5e5eXlpa2unuAFhQ2w3W41V1cZJMvp+Y3kgA8KO2G6zPJ7l+8vj6JHdOZxwAOLOtXNb+qSR/n+RlVfVoVd2Q5ENJ3lhV30jyhskyAMzMns026O7rNnjq9VOeBQA25E4XAAxBsAAYwqYvCQIssgOHjsx7hGfcfO0V8x5hrpxhATAEwQJgCIIFwBAEC4AhCBYAQxAsAIYgWAAMweewmJmz6fMrwPicYQEwBMECYAiCBcAQBAuAIQgWAEMQLACGIFgADEGwABiCYAEwBMECYAiCBcAQBAuAIQgWAEMQLACGIFgADEGwABiCYAEwBMECYAiCBcAQBAuAIezZyQ9X1beSPJnkqSQnu3t5GkMBwOl2FKyJ13b3E1P4PQCwIS8JAjCEnZ5hdZIvVVUn+cvuPnj6BlW1P8n+JHnxi1+8w92xmQOHjsx7BICZ2OkZ1qu7+xeSvCnJO6vql0/foLsPdvdydy8vLS3tcHcALKodBau7H5t8P5Hkc0mumsZQAHC6bQerqp5bVc//4eMkv5Lk6LQGA4D1dvIe1kVJPldVP/w9f9XdX5zKVABwmm0Hq7u/meQVU5wFADbksnYAhiBYAAxhGne6AJganyVkI86wABiCYAEwBMECYAiCBcAQBAuAIQgWAEMQLACGIFgADEGwABiCYAEwBMECYAiCBcAQBAuAIQgWAEMQLACGIFgADEGwABiCYAEwBMECYAiCBcAQBAuAIQgWAEMQLACGIFgADGHPvAfYjgOHjsx7BIBdd7b93XfztVfs6v6cYQEwBMECYAiCBcAQdhSsqrq6qr5eVY9U1XunNRQAnG7bwaqq85L8RZI3Jbk8yXVVdfm0BgOA9XZyhnVVkke6+5vd/T9JPp3kmumMBQCn2kmwLknyb+uWH52sA4Cpm/nnsKpqf5L9k8UfVNXXZ73PM7gwyRNz2O88OebFsYjH7ZjPAh+a3q/6ma1stJNgPZbk0nXLL5qsO0V3H0xycAf72bGqWunu5XnOsNsc8+JYxON2zItpJy8J/mOSy6rqJVX1nCS/neTz0xkLAE617TOs7j5ZVTcm+Zsk5yW5tbsfmtpkALDOjt7D6u67ktw1pVlmaa4vSc6JY14ci3jcjnkBVXfPewYA2JRbMwEwBMECYAgLEayq+qOqerCqDlfVl6rqp+c9026oqg9X1dcmx/65qnrBvGeatar6zap6qKqerqpz+hLgRbyXZ1XdWlUnqurovGfZLVV1aVXdW1UPT/7bfte8Z5qXhQhWkg9398u7e1+SLyT5g3kPtEvuTvLz3f3yJP+S5MCc59kNR5Ncm+S+eQ8ySwt8L89PJrl63kPsspNJ3tPdlyd5ZZJ3Lsif9Y9YiGB19/fXLT43yUJcadLdX+ruk5PFr2Ttw93ntO4+1t3zuJvKblvIe3l2931JvjvvOXZTdx/v7q9OHj+Z5FgW9DZ4M78109miqj6Y5HeS/GeS1855nHn4vSSfmfcQTM2Z7uX5i3OahV1SVXuTXJnk/vlOMh/nTLCq6m+T/NQZnrqpu+/s7puS3FRVB5LcmOT9uzrgjGx23JNtbsraywq37+Zss7KVY4ZzTVU9L8kdSd592qtGC+OcCVZ3v2GLm96etQ87nxPB2uy4q+p3k/xaktf3OfKhu2fxZ30u29K9PDk3VNX5WYvV7d19aN7zzMtCvIdVVZetW7wmydfmNctuqqqrk/x+kl/v7v+a9zxMlXt5LoiqqiS3JDnW3R+d9zzztBB3uqiqO5K8LMnTSb6d5B3dfc7/32hVPZLkx5P8+2TVV7r7HXMcaeaq6q1J/jzJUpLvJTnc3b8636lmo6renOTP8n/38vzgnEeauar6VJLXZO2f2ng8yfu7+5a5DjVjVfXqJH+X5EjW/g5LkvdNbo23UBYiWACMbyFeEgRgfIIFwBAEC4AhCBYAQxAsAIYgWAAMQbAAGML/AkPt7e3z0vb4AAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(7,5))\n",
"plt.hist(e,alpha=.6,bins=10)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"e.mean().round(2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Assessing Fit: $R^2$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Recall a projection is really a scaled down version of the vector being projected (e.g. $y$, the outcome vector). Keeping this in mind, there is an intuitive way to think of the model fit. \n",
"\n",
"\n",
"\n",
"We'll first need to center our outcome and prediction variables. \n",
"\n",
" \n",
" \n",
"\n",
"$$\\textbf{y}^* = \\textbf{y}-\\bar{y}$$ \n",
"\n",
"$$\\hat{\\textbf{y}}^* = \\hat{\\textbf{y}} - \\bar{\\hat{y}}$$\n",
"\n",
" \n",
" \n",
"\n",
"We can take the length of the projection and divided it by the length of the observed vector. \n",
"\n",
" \n",
" \n",
"\n",
"$$ \\frac{||\\hat{\\textbf{y}}^*||^2}{||\\textbf{y}^* ||^2} $$\n",
"\n",
" \n",
" \n",
"\n",
"If the lengths are the same, then this ratio will equal 1 — that is, if $\\textbf{y}$ lands on the plane created by $C(\\textbf{X})$, then our projection and $\\textbf{y}$ are equivalent (i.e. $\\textbf{y}$ is a linear dependent on $C(\\textbf{X})$). \n",
"\n",
"If the $\\textbf{y}$ doesn't project onto $C(\\textbf{X})$) at all then this ratio will equal 0 (meaning that all $C(\\textbf{X})$ is completely orthogonal to our outcome vector $\\textbf{y}$).\n",
"\n",
"This ratio is our **R-Squared**, which we use to assess model fit."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.9233409734326696"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"y_hat_center = y_hat - y_hat.mean()\n",
"y_center = y - y.mean()\n",
"\n",
"# R Squared\n",
"r_squared = y_hat_center.dot(y_hat_center)/y_center.dot(y_center)\n",
"r_squared"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An alternative specification of $R^2$ uses the length of the error vector instead.\n",
"\n",
" \n",
" \n",
"\n",
"$$ 1- \\frac{||\\hat{ e }||^2}{||\\textbf{y}^* ||^2}$$\n",
"\n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.9233409734326684"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"1 - (e.dot(e)/y_center.dot(y_center))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Calculating Uncertainty\n",
"\n",
"Calculating the uncertainty in our estimates of $\\beta$ relies on the Gauss-Markov assumptions and some manipulations of the resulting equations. (We won't not cover these in class, but see [here](https://web.stanford.edu/~mrosenfe/soc_meth_proj3/matrix_OLS_NYU_notes.pdf) for more details on how these estimates are derived and why they are unbiased given the assumptions). That said, let's think through how to calculate the uncertainty. \n",
"\n",
"Recall that our orthogonal vector $\\textbf{e}$ is the error in the model. \n",
"\n",
"We can calculate our estimate of uncertainty by taking the length of $\\textbf{e}$ and dividing it by our sample size $N$ adjusted for our degrees of freedom (i.e. the number of parameters we had $p$). \n",
"\n",
" \n",
" \n",
"\n",
"$$\\hat{\\sigma}^2 = \\frac{||e||^2}{n-p} = \\frac{e^T\\cdot e}{n-p} $$\n",
"\n",
" \n",
" \n",
"\n",
"This will provide use with a scalar that we'll then apply to our matrix of explanatory variables to calculate the **standard error** of our coefficients $\\beta$.\n",
"\n",
" \n",
" \n",
"\n",
"$$ cov(\\beta) = \\hat{\\sigma}^2 (\\textbf{X}^T \\textbf{X})^{-1} $$\n",
"\n",
"\n",
" \n",
" \n",
"\n",
"$$ D = diag(cov(\\beta)) $$\n",
"\n",
" \n",
" \n",
"\n",
"$$ se(\\beta) = \\sqrt{D} $$\n",
"\n",
"\n",
" \n",
" \n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.9598505223222045"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Let's calculate the uncertainty\n",
"n = X.shape[0]\n",
"p = X.shape[1]\n",
"sigma2 = e.T.dot(e)/(n-p)\n",
"sigma2"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0.00960405, -0.00020442],\n",
" [-0.00020442, 0.0075405 ]])"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cov_B = sigma2*la.inv(X.T.dot(X))\n",
"cov_B"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.00960405, 0.0075405 ])"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"D = cov_B.diagonal()\n",
"D"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.09800024, 0.08683606])"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"se = np.sqrt(D)\n",
"se"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### How did we do?\n",
"\n",
"Let's compare our calculations those yielded by the `statsmodels` module, which can estimate a range of statistical models for us in Python."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Our Manually Generated Values:\n",
"\n",
" Intercept: 0.9809 (0.098)\n",
" x: 2.9834 (0.087)\n",
"\n",
" R-Squared: 0.923\n",
"\n"
]
}
],
"source": [
"# Use \n",
"\n",
"print(f'''\n",
"Our Manually Generated Values:\n",
"\n",
" Intercept: {B[0].round(4)} ({se[0].round(3)})\n",
" x: {B[1].round(4)} ({se[1].round(3)})\n",
"\n",
" R-Squared: {r_squared.round(3)}\n",
"''')"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"
OLS Regression Results
\n",
"
\n",
"
Dep. Variable:
y
R-squared:
0.923
\n",
"
\n",
"
\n",
"
Model:
OLS
Adj. R-squared:
0.923
\n",
"
\n",
"
\n",
"
Method:
Least Squares
F-statistic:
1180.
\n",
"
\n",
"
\n",
"
Date:
Mon, 28 Oct 2019
Prob (F-statistic):
1.84e-56
\n",
"
\n",
"
\n",
"
Time:
13:22:17
Log-Likelihood:
-138.83
\n",
"
\n",
"
\n",
"
No. Observations:
100
AIC:
281.7
\n",
"
\n",
"
\n",
"
Df Residuals:
98
BIC:
286.9
\n",
"
\n",
"
\n",
"
Df Model:
1
\n",
"
\n",
"
\n",
"
Covariance Type:
nonrobust
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
coef
std err
t
P>|t|
[0.025
0.975]
\n",
"
\n",
"
\n",
"
Intercept
0.9809
0.098
10.009
0.000
0.786
1.175
\n",
"
\n",
"
\n",
"
x
2.9834
0.087
34.357
0.000
2.811
3.156
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
Omnibus:
5.027
Durbin-Watson:
1.860
\n",
"
\n",
"
\n",
"
Prob(Omnibus):
0.081
Jarque-Bera (JB):
5.131
\n",
"
\n",
"
\n",
"
Skew:
-0.308
Prob(JB):
0.0769
\n",
"
\n",
"
\n",
"
Kurtosis:
3.924
Cond. No.
1.13
\n",
"
\n",
"
Warnings: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y R-squared: 0.923\n",
"Model: OLS Adj. R-squared: 0.923\n",
"Method: Least Squares F-statistic: 1180.\n",
"Date: Mon, 28 Oct 2019 Prob (F-statistic): 1.84e-56\n",
"Time: 13:22:17 Log-Likelihood: -138.83\n",
"No. Observations: 100 AIC: 281.7\n",
"Df Residuals: 98 BIC: 286.9\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept 0.9809 0.098 10.009 0.000 0.786 1.175\n",
"x 2.9834 0.087 34.357 0.000 2.811 3.156\n",
"==============================================================================\n",
"Omnibus: 5.027 Durbin-Watson: 1.860\n",
"Prob(Omnibus): 0.081 Jarque-Bera (JB): 5.131\n",
"Skew: -0.308 Prob(JB): 0.0769\n",
"Kurtosis: 3.924 Cond. No. 1.13\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model = smf.ols('y ~ x', data=pd.DataFrame(dict(y=y,x=X[:,1]))).fit()\n",
"model.summary()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# OLS as transformation matrices\n",
"\n",
" \n",
" \n",
"\n",
"We've seen that we can generate a projection of $\\textbf{y}$ onto $C(\\textbf{X})$ using the following formula.\n",
"\n",
" \n",
" \n",
"\n",
"$$\\hat{\\textbf{y}} = \\textbf{X}\\beta$$\n",
"\n",
" \n",
" \n",
"\n",
"Where\n",
"\n",
" \n",
" \n",
"\n",
"$$ \\beta = (\\textbf{X}^T \\textbf{X})^{-1} \\textbf{X}^T \\textbf{y}$$\n",
"\n",
" \n",
" \n",
"\n",
"So we can re-write this as \n",
"\n",
" \n",
" \n",
"\n",
"$$\\hat{\\textbf{y}} = \\textbf{X}(\\textbf{X}^T \\textbf{X})^{-1} \\textbf{X}^T \\textbf{y}$$\n",
"\n",
" \n",
" \n",
"\n",
"$$\\hat{\\textbf{y}} = \\textbf{H} \\textbf{y}$$\n",
"\n",
" \n",
" \n",
"\n",
"$\\textbf{H}$ is known as the \"Hat Matrix\". $\\textbf{H}$ transforms the vector $\\textbf{y}$ into $\\hat{\\textbf{y}}$."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"H = X.dot(la.inv(X.T.dot(X)).dot(X.T))\n",
"y_hat2 = H.dot(y)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3oAAAEKCAYAAABE7ieQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X+c3Hdd4PHXO5stnRTbLbbKZZtcKifhKIFGFlqsKKXVgNISq/yo4IH30Jx68ksuQKDS1qtXJAiiKHcRC3qttVDi2iISQIHTagObbts0lJzyq+0GzgDdKmRLN5v3/TEz6exkZnZ2d3a/szOv5+Mxj8x85jvf73v322T73vfn8/5EZiJJkiRJ6h2rig5AkiRJktRZJnqSJEmS1GNM9CRJkiSpx5joSZIkSVKPMdGTJEmSpB5joidJkiRJPcZET5IkSZJ6jImeJEmSJPUYEz1JkiRJ6jGriw5gPs4444zcsGFD0WFIkiRJUiH27dv3jcw8c67jVlSit2HDBsbGxooOQ5IkSZIKERFfbec4p25KkiRJUo8x0ZMkSZKkHmOiJ0mSJEk9xkRPkiRJknqMiZ4kSZIk9ZgV1XVTkiRJkpbS6PgEO/cc5NDkFGuHSmzfspGtm4eLDmveCq3oRcTrIuJARNwTETdGxMlFxiNJkiSpf42OT7Bj934mJqdIYGJyih279zM6PlF0aPNWWKIXEcPAq4GRzHwKMAC8tKh4JEmSJPW3nXsOMjU9M2tsanqGnXsOFhTRwhW9Rm81UIqI1cAa4FDB8UiSJEnqU4cmp+Y13s0KS/QycwJ4B3Af8DXgocz8eFHxSJIkSepva4dK8xrvZkVO3TwdeCFwNrAWOCUiXt7guG0RMRYRY4cPH17uMCVJkiT1ie1bNlIaHJg1VhocYPuWjQVFtHBFdt28GPhyZh4GiIjdwA8D19celJm7gF0AIyMjudxBSpIkSeouS9UZs3qOXui6WWSidx9wfkSsAaaAi4CxAuORJEmS1OWqnTGrTVOqnTGBjiV7KzGxq1dYopeZeyPiZuAO4CgwTqVyJ0mSJKm/Vat2E5NTDEQwk8nwUIkjjxxt2hmzFxK0Til0w/TMvBK4ssgYJEmSJHWX+qrdTJZXcE206H65EjtjLqWit1eQJEmSpFka7Wc3l5XYGXMpmehJkiRJ6irzrc6t1M6YS6nQqZuSJEmSetdCu2OuHSq1nKY5VBrklMesXvGdMZeSiZ4kSZKkjltMd8ztWzbO+myt0uAAV116jondHJy6KUmSJKnjGq2zq3bHnMvWzcNce9kmhivr7gYiABgeKnHtZZtM8tpgRU+SJElSxzVbZ9fu+rte2c+uKFb0JEmSJHVcsy6YdsdcHiZ6kiRJkjpu+5aNlAYHZo3ZHXP5OHVTkiRJUsdVp10upOumFs9ET5IkSdKScJ1dcZy6KUmSJEk9xkRPkiRJknqMiZ4kSZIk9RgTPUmSJEnqMSZ6kiRJktRj7LopSZIkdZnR8QmuuuUAk1PTAJy+ZpArLznHDpZqm4meJEmS1EVGxyfY/qG7mD6Wx8cePDLN9pvvAjDZU1sKnboZEUMRcXNEfCEi7o2IZxUZjyRJklS0nXsOzkryqqZnkp17DhYQkVaioit67wY+lpk/GxEnAWsKjkeSJEnquNHxCXbuOcihySnWDpXYvmVj08rcocmppudp9Z5Uq7BELyJOA34UeCVAZj4CPFJUPJIkSVKn1a+1A5iYnGLH7v1A42mYa4dKTDRJ6NYOlZYmUPWcIqdung0cBt4fEeMR8b6IOKX+oIjYFhFjETF2+PDh5Y9SkiRJWoDR8Ql27N4/K8mrmpqeaToNc/uWjQyuihPGBweC7Vs2djxO9aYiE73VwA8B783MzcB3gDfVH5SZuzJzJDNHzjzzzOWOUZIkSVqQnXsOMjU90/T9ZtMwt24eZueLnsZQafD42OlrBtn5s0+zEYvaVuQavQeABzJzb+X1zTRI9CRJkqSVaK71dK2mYW7dPGxSp0UprKKXmV8H7o+Iav35IuDzRcUjSZIkdVKrRK40OOA0TC2pQrdXAF4F3BARdwPnAv+j4HgkSZKkjti+ZSOlwYETxk9fM8i1l22yYqclVej2Cpl5JzBSZAySJEnqby/7o3/kti9+6/jrC57wOG74pcVv71xN5NrdVkHqpMg8cTPGbjUyMpJjY2NFhyFJkqQeUZ/kVXUq2ZM6LSL2ZeacxbKiN0yXJEmSOq7dDcobJXmtxqWVwkRPkiRJPaW6f111a4O5NiiXelHRzVgkSZKkjmq0f12rDcqlXmSiJ0mSpJ7SbP+6RuMXPOFxDY9tNi6tFCZ6kiRJ6inN9q9rNH7DLz3rhKTORizqBa7RkyRJUk/ZvmXjrDV60HqDcpM69SITPUmSJBWm3e6Y8+H+dZKJniRJkgowOj7B1bce4MEj08fHOtkdc+vmYRM79TXX6EmSJGlZVbc/qE3yquyOKXWGiZ4kSZKWVaPtD2o165opqX1O3ZQkSdKCVdfYTUxOEUBWxk9fM8iVl5zTcPrkXIlcs66ZktpnoidJkqQFuWJ0Pzfcft/x5C5r3nvwyDTbb74LOHG93dqhEhNNkr1W3TEltc+pm5IkSZq30fGJWUleI9Mz2XC93fYtGykNDpwwPlQa5NrLNtlEReoAK3qSJEmat517DrZM8qoaTdN0+wNp6ZnoSZIkCShPxbxx7/3MZDIQweXnreOarZsaHttuw5Rm6+3c/kBaWoUnehExAIwBE5n5gqLjkSRJ6nW1DVQGIpjJ5JSTBvjOI492wpzJ5Prb7wNomOy1WmdXNTgQrreTClJ4oge8BrgXOLXoQCRJknrN6PgEV91ygMmp8p51awZXMX0smZ4pT7ycyfKftUlerRv33t8w0du+ZSM7du9vuk1Cq66bkpZeoYleRJwF/BTwW8CvFxmLJElSLxkdn+DqWw+csCn5kelj8zpPNRGs5zo7qbsVXdH7XeANwPcUHIckSVLPGB2faFltm4+BiKbvuc5O6l6FJXoR8QLgXzJzX0Q8p8Vx24BtAOvXr1+m6CRJkrpbdZ1do2razj0HO5LkAVx+3rqOnEfS8iqyoncBcGlE/CRwMnBqRFyfmS+vPSgzdwG7AEZGRtrp4itJktSTapuoBI9uUD4xOcWO3fuBcpWt3Y6Y9QZWBceOJQlzdt2U1N0KS/QycwewA6BS0ftv9UmeJEmSGq+3q//t99T0DDv3HGTr5uE5O2KuCjitNMiDR6aPd90cdo2d1FOKXqMnSZKkFuaz3q5ayWvVEdNumFJ/6IpELzM/DXy64DAkSZK6znzW21U3J7cjpqSuSPQkSZL6Qbk6dzdTlS0OVgX83HnrW66Da3e9XWlwYNbm5HbElPrbqqIDkCRJ6gej4xP8+k13Hk/yAI4lXH/7fVwxur/p56pVukaqGx8MD5W49rJNJnaSjrOiJ0mS1AGttjuA8jTKZluV37j3/qZVvWbr7YZKg1x1qWvtJDVmoidJkrRI9Q1T6rc7gNZTMGey+Q5SrreTtBAmepIkSW1qVrVr1DCldrsDoOWWBwMRDcerXG8nab5coydJktSGatVuYnKK5NGq3ej4RNNqXe349i0bm/6P1+Xnret8wJL6momeJElSG1pV7Zo1TKkd37p5mHe+5FxKg4/+79eqgJef37rrpiQthFM3JUmS2tCqaveul5x7QsOU+u0OwCmYkpaPFT1JkqQ2tKrabd08zLWXbWJ4qETgdgeSimdFT5Ik9YW5tj+YS6NtDmqrdlbrJHUTEz1JktTz2tn+YC5ucyBpJTHRkyRJPa+d7Q/aYdVO0krRdI1eRJwcEa+IiEuj7I0R8ZGIeHdEnLGcQUqSJC1GO9sfSFIvaVXR+1NgGjgFeD1wD/Ae4EeADwAvWOrgJEmSai10nV2zzcqbNViRpJWuVaL35Mx8SkSsBh7IzB+rjH8sIu5ahtgkSZKOW8w6u7kaqUhSr2mV6D0CkJlHI+JQ3XszDY6XJElalFYVu8Wss7ORiqR+0yrROysifg+ImudUXvuvoiRJWpTR8QmuvvUADx6ZBqA0uIqjx5LpmQROrNgtdp2djVQk9ZNWid72mudjde/Vv563iFhHeR3g9wMJ7MrMdy/2vJIkqbs0qtIBbL/5ruNJHcDU9LETPltbsXOdnSS1r1Wi978z88R/cYGIGOrAtY8Cr8/MOyLie4B9EfGJzPx8B84tSZIKUpvYDa0Z5NsPH2X62Owq3cmDq2Ylea1UK3aus5Ok9jXdXgEYi4jz6gcj4heBOxZ74cz8WmbeUXn+b8C9OCVUkqQVrdowZWJyigQePDJ9PMmrmpqeOT5dsx3Vit3WzcNce9kmhodK5XUkQyWuvWyT0zElqYFWFb1XA7si4rPAG4F/D/wh8ADwo50MIiI2AJuBvQ3e2wZsA1i/fn0nLytJkjrgitH93Lj3fmayvQrdfNRX7FxnJ0ntaVrRy8y/B54O/D/gi8AtwJWZ+aLMfKBTAUTEY4EPA6/NzH9tEMeuzBzJzJEzzzyzU5eVJEkdcMXofq6//b55J3lDpUEGB+KE8VXA6WsGrdhJ0iK1qugB/CxwOfBe4MeBl0TEWGZ+qxMXj4hBykneDZm5uxPnlCRJy+fGvffP+zOlwQGuuvQcgFldN4dKg1x16TkmdpLUAU0TvYj4JPAwcHFmfjkirgD+K/C5iPjtzNy1mAtHRAB/DNybme9czLkkSdLiVBuoTExOMRDBTCbDbew1104lb3AgOOWk1Tw0NX3C/nUmdZK0NFpV9P4gM/+i+qLSgfP3I+JDwO8Ai0r0gAuAnwf2R8SdlbE3Z+ZHF3leSZLUhtrkLijvdQSPJm/1+9g1Uk0KGwlwY3JJKkjTRK82yasb/zrwssVeuLIG8MTJ+ZIkaclVu2NWtypoVper3ceukcvPW8f1t993wvjLz1/PNVs3dSpcSdI8zbVGT5Ik9aCdew7O2o+ulUMNNimvqiZz1a6bAxFcft46kzxJKpiJniRJfahV8lavuo9dM9ds3WRiJ0ldptWG6ZIkqUfNlbxV1e9jJ0laGeas6EXEmZQ3TH8ycHJ1PDOfu4RxSZKkBhbaHbPe9i0bZ63RA443ZFnMeSVJ3aGdqZs3ADcBPwX8MvAK4PBSBiVJksqqid2hySlOKw3ynUeOMj1Tbp0yn+6Y9arHVc9td0xJ6i2Rc+x/ExH7MvPpEXF3Zj61Mva5zHzGskRYY2RkJMfGxpb7spIkLatm2x7MZXioxG1vcsKNJPWySn42Mtdx7VT0pit/fi0ifgo4BDxuMcFJkqTG2t32oJH5NFiRJPW2dhK9ayLiNOD1wO8DpwKvW9KoJEnqU/PZ9qBeuw1WJEm9b85ELzM/Unn6EHDh0oYjSVLvqV1nN9dauIVW5eyOKUmq1TTRi4g3ZObbI+L3aTBzJDNfvaSRSZK0QtUmdkNrBvn2w0eZPtZe45S1QyUmWiR7g6uCx568mgePTNsdU5LUVKuK3r2VP+1+IklSG0bHJ7jqlgNMTk0fH3vwyPQJx01Nz7Bzz8GGiVmrbQ9M6CRJ7Wqa6GXmrZWn+zPzjmWKR5KkFaF+OuaFTzqTD++baHt9XbMpmm57IEnqhHaasfxORDweuBm4KTPvWeKYJEnqavWdMScmp7jh9vvm1SGzVeOUrZuHTewkSYvSTjOWCyuJ3ouB/xURp1JO+K5Z8ugkSSrAXM1TGnXGnE+SZ+MUSdJSa6eiR2Z+Hfi9iPgU8AbgrYCJniSpZzTbpLxR85T5dsYcHAhOOWk1D01NOxVTkrQs5kz0IuI/Ai8Bfgb4JnAT5T31JEnqCXNtUl7fPKVZZ8xo8NnT1wxy5SXnmNhJkpZVOxW964A/B7Zk5qFOXjwinge8GxgA3peZb+vk+SVJgoVNxaxXW8Vr1BmzNDjAzzx9mE994bBNVCRJhWtnjd6zluLCETEA/AHw48ADwOci4pbM/PxSXE+S1F86PRWztnmKnTElSd2urTV6S+SZwD9n5pcAIuLPgRcCJnqSpEXp1FTMqkbNU+yMKUnqZkUmesPA/TWvHwDOKygWSdIKU1uxG4hgJvP4huKdmIrpJuWSpJWsnWYsL8rMD801tlQiYhuwDWD9+vXLcUlJUhcbHZ/g6lsP8OCR6eNjM1mu2VWnZbazablTMSVJvaydit4OoD6pazQ2XxPAuprXZ1XGZsnMXcAugJGRkflsUyRJWsFqK3ZVAaxaFcwca/7jYGp65niFrxmnYkqSel3TRC8ing/8JDAcEb9X89apwNEOXPtzwA9GxNmUE7yXAj/XgfNKkla4+jV2VQktk7yqmUxKgwNOxZQk9a1WFb1DwBhwKbCvZvzfgNct9sKZeTQifg3YQ3l7hesy88BizytJWvnaWWPXSu1aPadiSpL6UdNELzPvAu6KiD/LzOlmxy1GZn4U+OhSnFuStHK1s91BM9VpmU7FlCT1s1VtHLMhIm6OiM9HxJeqjyWPTJLUt2obpbQjovzn8FCJay/bZIInSep77TRjeT9wJfAu4ELgF2gvQZQk9ZFq85ROTJVstN1BrVNOGuDIIzNOyZQkqYl2Er1SZv5NRERmfhW4KiL2AW9d4tgkSStEffOU6jYHwIKSsNrtDmq7bg5EcPl567hm66YORC1JUu9qJ9H7bkSsAv6p0jxlAnjs0oYlSVpJGjVPmZqeYeeegwuutrnGTpKkhWsn0XsNsAZ4NfDfgecCr1jKoCRJ3eeK0f3cuPd+ZjJPqKw1a56ymKYqkiRp4eZM9DLzc5Wn36a8Pk+S1MMarbUb++q3uP72+44fM5N5/PU1Wzexdqg0a4pl1XybqkiSpM5otWH6rZT3lm0oMy9dkogkSctudHyCq245wOTU7N10qmvtvnu0cVOUG/fezzVbNzVsnlLd5kCSJC2/VhW9d1T+DOCPgF9c+nAkSUutvmJ34ZPO5KbP3s/0sca/22u1cflMlj9T2zzFDcolSSpeqw3TP1N9HhHfrn0tSVqZGnXHvOH2+5pP35jDQHUDO2yeIklSN2mnGQu0mMIpSepOjdbaNeqO2c4/8KXBVUxNHzth/PLz1nUoWkmS1Emt1ug9rublQEScTnkaJwCZ+a2lDEyS1L76NXZrBlcxfSyZnimncdW1dq2mYTZTGhzg2ss2MfbVbzXtuilJkrpLZDb+XW5EfJnyL3qjwduZmT+wlIE1MjIykmNjY8t9WUnqKqPjE1x96wEePPJoUvfw9DFOrLedaCDi+Lq6dgyVBrnq0nOckilJUpeIiH2ZOTLXca3W6J3d2ZAkSQvRrCNm1ZEGUyqbmcmkNDhwQnfMn3n6MB+562vHr3H6mkGuvMQET5KklardNXqSpGVWX7nrhOGatXr13TGdhilJUu8w0ZOkLlTfHbMTqvva2R1TkqTe16oZy0eBX83MryxfOJLUPxp1xazdj26xSd6qgNNKg0wemXZfO0mS+kyrit77gY9HxJ8Ab8/Mzs0dkqQ+12g/ux279wPl/egOTU7N63yrAh6z+tEtEFxjJ0lSf2vVjOVDEfHXwG8AYxHxv+HRpm6Z+c6FXjQidgKXAI8AXwR+ITMnF3o+SVppGlXspqZn2LnnIFs3D7N2qMREm8menTElSVK9udboPQJ8B3gM8D3QVvfudnwC2JGZRyPit4EdwBs7dG5J6nrNKnbV8e1bNjZdo2e1TpIkzaXVGr3nAe8EbgF+KDOPdOqimfnxmpe3Az/bqXNL0nJqtc6ulWYVu7VDJYBZa/Xme25JkqRWFb23AC/KzANLHMN/Bm5q9mZEbAO2Aaxfv36JQ5Gk9s21zq6VRhW7alfMKrtjSpKkhYrMXJoTR3wSeHyDt96SmX9ZOeYtwAhwWbYRyMjISI6NjXU2UElqolqtq6+8XfCEx3HDLz2LC972tw2rcsNDJW5703PbPr8VO0mS1K6I2JeZI3Mdt2T76GXmxa3ej4hXAi8ALmonyZOk5dRqH7vbvvgtXvZH/zjnOru5WLGTJElLZVURF62s/3sDcGkn1/5JUqfMtY/dbV/81vH1dPWajUuSJC2XQhI94D2Uu3h+IiLujIj/WVAcktRQO1W57Vs2UhocmDVWv85OkiSpCEs2dbOVzPwPRVxXUn+pXWM3EMFMJsNtroVrZx87O2NKkqRuVUiiJ0lLrX6N3UxlKXC7nTFb7WMH5YYs1XOY2EmSpG5joidpRWvWubLVGrup6Rl27jnYMkGrrdY167opSZLUrUz0JK1YV4zu54bb76Patre2WjfXGrt21uBZrZMkSSuViZ6kFaO2endaaZDJqekTjqlW6+ZaY2dnTEmS1MtM9CR1pfopmRc+6Uw+vG/i+HTMRkle1aHJKd71knObrrGzM6YkSep1JnqSCtVojR0wK0mbmJyaNUVzLmuHSiessZtv101JkqSVzERPUmHqO2NW19g9ZvWqEypx7SZ5AceTRdfYSZKkfmWiJ6kwjTpjTk3PNO2WOZcAXnb+epM7SZLU90z0JC2ZZlsfVLXT+bJWMLuyN7gqeOzJq5k8Mu1m5ZIkSTVM9CQtiWbTMuHRPeqadcY8fc0gD08fm1XZKw0O8DNPH+ZTXzjcNHGUJElSmYmepCXRbFpm7Ubl27dsPKEzZmlwgCsvOef4OUzqJEmS5s9ET9K8zDUds6rZtMza8drOmI3OZ2InSZK0MCZ6ktrWznTMqmbTMus3KrczpiRJUuetKjoASStHq+mY9bZv2UhpcGDWmBuVS5IkLQ8relIfanf6Zb12pmNWzTUtU5IkSUvHRE/qM/OZflmv3emYVU7LlCRJKkahUzcj4vURkRFxRpFxSP1kPtMv6zkdU5IkaWUorKIXEeuAnwDuKyoGqR/NZ/plPadjSpIkrQxFTt18F/AG4C8LjEHqO/OdflnP6ZiSJEndr5CpmxHxQmAiM+8q4vrSSjU6PsEFb/tbzn7TX3HB2/6W0fGJeZ/D6ZeSJEm9b8kqehHxSeDxDd56C/BmytM22znPNmAbwPr16zsWn7TSLKaJSi2nX0qSJPW+yMzlvWDEJuBvgCOVobOAQ8AzM/PrrT47MjKSY2NjSxyhVLxG2x/s3HOw4ZTL4aESt73puQVEKUmSpOUWEfsyc2Su45Z9jV5m7ge+r/o6Ir4CjGTmN5Y7FqkbNavc1XfKrGqniYokSZL6S6HbK0g6UbPtDwYiGh7fbhMVSZIk9Y/CN0zPzA1FxyB1k2YVuplMSoMDs5JAm6hIkiSpESt6UpdpVqEbHipx7WWbGB4qETWvbaIiSZKkeoVX9CTNtn3LxhPW5FUrd+5hJ0mSpHaY6Eldxu0PJEmStFgmelIXsnInSZKkxXCNniRJkiT1GBM9SZIkSeoxJnqSJEmS1GNM9CRJkiSpx5joSZIkSVKPseumesLo+ITbEUiSJEkVJnpa8UbHJ2ZtMD4xOcWO3fsBTPYkSZLUl5y6qRVv556Dx5O8qqnpGXbuOVhQRJIkSVKxTPS04h2anJrXuCRJktTrTPS04q0dKs1rXJIkSep1Jnpa8bZv2UhpcGDWWGlwgO1bNhYUkSRJklQsm7Foxas2XLHrpiRJklRmoqeesHXzsImdJEmSVFHY1M2IeFVEfCEiDkTE24uKQ5IkSZJ6TSEVvYi4EHgh8LTM/G5EfF8RcUiSJElSLyqqovcrwNsy87sAmfkvBcUhSZIkST2nqETvicCzI2JvRHwmIp7R7MCI2BYRYxExdvjw4WUMUZIkSZJWpiWbuhkRnwQe3+Ctt1Su+zjgfOAZwAcj4gcyM+sPzsxdwC6AkZGRE97X0hgdn7CLpSRJkrRCLVmil5kXN3svIn4F2F1J7D4bEceAMwBLdl1gdHyCHbv3MzU9A8DE5BQ7du8HMNmTJEmSVoCipm6OAhcCRMQTgZOAbxQUi+rs3HPweJJXNTU9w849BwuKSJIkSdJ8FLWP3nXAdRFxD/AI8IpG0zZVjEOTU/MalyRJktRdCkn0MvMR4OVFXFtzWztUYqJBUrd2qFRANJIkSZLmq7AN09W9tm/ZSGlwYNZYaXCA7Vs2FhSRJEmSpPkoauqmuli14YpdNyVJkqSVyURPDW3dPGxiJ0mSJK1QTt2UJEmSpB5joidJkiRJPcZET5IkSZJ6jImeJEmSJPUYEz1JkiRJ6jEmepIkSZLUY9xeYRFGxyfca06SJElS1zHRW6DR8Ql27N7P1PQMABOTU+zYvR/AZE+SJElSoZy6uUA79xw8nuRVTU3PsHPPwYIikiRJkqQyE70FOjQ5Na9xSZIkSVouJnoLtHaoNK9xSZIkSVouJnoLtH3LRkqDA7PGSoMDbN+ysaCIJEmSJKnMZiwLVG24YtdNSZIkSd2mkEQvIs4F/idwMnAU+NXM/GwRsSzG1s3DJnaSJEmSuk5RUzffDlydmecCb628liRJkiR1QFGJXgKnVp6fBhwqKA5JkiRJ6jlFrdF7LbAnIt5BOdn84YLikCRJkqSes2SJXkR8Enh8g7feAlwEvC4zPxwRLwb+GLi4yXm2AdsA1q9fv0TRSpIkSVLviMxc/otGPAQMZWZGRAAPZeapc31uZGQkx8bGlj5ASZIkSepCEbEvM0fmOq6oqZuHgB8DPg08F/indj60b9++b0TEV5u8fQbwjY5Ep8XwPnQP70V38D50B+9Dd/A+dAfvQ/fwXnSHlXYf/n07BxVV0fsR4N2UE82HKW+vsG+R5xxrJ7PV0vI+dA/vRXfwPnQH70N38D50B+9D9/BedIdevQ+FVPQy8++BpxdxbUmSJEnqdUVtryBJkiRJWiK9lOjtKjoAAd6HbuK96A7eh+7gfegO3ofu4H3oHt6L7tCT96GQNXqSJEmSpKXTSxU9SZIkSRI9luhFxLkRcXtE3BkRYxHxzKJj6lcR8aqI+EJEHIiItxcdTz+LiNdHREbEGUXH0q8iYmfl78PdEfEXETFUdEz9JCKeFxEHI+KfI+JNRcfTjyJiXUR8KiI+X/m58JqiY+pnETEQEeMR8ZGiY+lXETEUETdXfjbcGxHPKjqmfhQRr6v8m3RPRNwYEScXHVMn9VSiB7wduDozzwXeWnmtZRYRFwIvBJ6WmecA7yg4pL4VEeuAnwDuKzqWPvcJ4CmZ+VTg/wI7Co6nb0TEAPAHwPOBJwOXR8TQ/LG1AAAJ00lEQVSTi42qLx0FXp+ZTwbOB/6r96FQrwHuLTqIPvdu4GOZ+STgaXg/ll1EDAOvBkYy8ynAAPDSYqPqrF5L9BI4tfL8NMobs2v5/Qrwtsz8LkBm/kvB8fSzdwFvoPx3QwXJzI9n5tHKy9uBs4qMp888E/jnzPxSZj4C/DnlX0RpGWXm1zLzjsrzf6P8P7XDxUbVnyLiLOCngPcVHUu/iojTgB8F/hggMx/JzMlio+pbq4FSRKwG1tBjuUOvJXqvBXZGxP2Uq0j+1rwYTwSeHRF7I+IzEfGMogPqRxHxQmAiM+8qOhbN8p+Bvy46iD4yDNxf8/oBTDAKFREbgM3A3mIj6Vu/S/kXgMeKDqSPnQ0cBt5fmUL7vog4peig+k1mTlDOF+4DvgY8lJkfLzaqzipkw/TFiIhPAo9v8NZbgIuA12XmhyPixZR/U3LxcsbXL+a4D6uBx1GenvMM4IMR8QNpi9eOm+M+vJnytE0tg1b3IjP/snLMWyhPYbthOWOTukVEPBb4MPDazPzXouPpNxHxAuBfMnNfRDyn6Hj62Grgh4BXZebeiHg38CbgN4oNq79ExOmUZ3icDUwCH4qIl2fm9cVG1jkrLtHLzKaJW0T8KeV55wAfwmkJS2aO+/ArwO5KYvfZiDgGnEH5t1fqoGb3ISI2Uf6H666IgPJUwTsi4pmZ+fVlDLFvtPo7ARARrwReAFzkLz2W1QSwrub1WZUxLbOIGKSc5N2QmbuLjqdPXQBcGhE/CZwMnBoR12fmywuOq988ADyQmdWq9s2UEz0tr4uBL2fmYYCI2A38MNAziV6vTd08BPxY5flzgX8qMJZ+NgpcCBARTwROAr5RaER9JjP3Z+b3ZeaGzNxA+YfKD5nkFSMinkd5qtSlmXmk6Hj6zOeAH4yIsyPiJMoL7W8pOKa+E+XfOP0xcG9mvrPoePpVZu7IzLMqPxdeCvytSd7yq/wsvj8iNlaGLgI+X2BI/eo+4PyIWFP5N+oieqwpzoqr6M3hl4B3VxZUPgxsKziefnUdcF1E3AM8ArzCCob63HuAxwCfqFRYb8/MXy42pP6QmUcj4teAPZQ7ql2XmQcKDqsfXQD8PLA/Iu6sjL05Mz9aYExSkV4F3FD5BdSXgF8oOJ6+U5k2ezNwB+VlFePArmKj6qzw/78lSZIkqbf02tRNSZIkSep7JnqSJEmS1GNM9CRJkiSpx5joSZIkSVKPMdGTJEmSpB5joidJPSTK/j4inl8z9qKI+FjdcV+JiDNqXj8nIj4yx7nPrWy0PJ94Lo2IFbsRcES8eY73PxoRQ22eaygifrUzkc15rRsj4u6IeN1yXK9yzddGxJqa121/byRJnef2CpLUYyLiKcCHgM2U90sdB56XmV+sOeYrwEhmfqPy+jnAf8vMF7Q47ysrn/m1JQt+kSJiIDNnOni+b2fmYxuMB+Wfocfmca4NwEcy8ykN3ludmUcXE2vNuR4P/H1m/od5fGbR16//b0qSVCwrepLUYzLzHuBW4I3AW4E/rU3y5hIRz4yIf4yI8Yj4h4jYWNnU9zeBl0TEnRHxkrrP3B4R59S8/nREjETEKyPiPZWxSyJib+W8n4yI729w7ZMj4v0Rsb9y3IWV8ePnqbz+SCU5JSK+HRG/ExF3Ac+qO9+nI+JdETEWEfdGxDMiYndE/FNEXFNz3GhE7IuIAxGxrTL2NqBU+XpviIgNEXEwIv4UuAdYV62MVs57dyX+UyrnqU/o3gY8oXK+nZUq6t9FxC3A55vFUfM1/lZE3FX5Xn9/ZfxFEXFPZfz/VA7/ODBcuc6zK5XY2yvx/UVEnF7zvfndiBgDXhMRH4iI91aO/VIlvusq37cP1MTy3sr380BEXF0ZezWwFvhURHyqMna8ahwRv16J856IeG1lbEPl3H9UOdfHI6JU/9+EJGmBMtOHDx8+fPTYAzgFOAjsBx7T4P2vVN67s/L4Z8rVJoBTgdWV5xcDH648fyXwnibXex1wdeX5vwMO1n8GOJ1HZ5L8IvA7Dc7zeuC6yvMnAfcBJ9dfG/gI8JzK8wRe3CSuTwO/XXn+GuBQJb7HAA8A31t573GVP0uUk7jq+LdrzrUBOAacX/d9PKPy/BrgHcAfADsaxLIBuKfm9XOA7wBn14w1iyOBSyrP3w5cUXm+HxiuPB9qcp27gR+rPP9N4Hdrvjd/WHPcB4A/BwJ4IfCvwCbKvxTeB5xbF+NA5RxPrf9e1L4Gnl6J8xTgscABytXmDcDRmvN+EHh50X93fPjw4aNXHquRJPWczPxORNxEOVH5bpPDLsy6qZuV8dOAP4mIH6ScYAy2cckPUq4kXQm8GLi5wTFnATdFxL8DTgK+3OCYHwF+v/I1fCEivgo8cY5rzwAfbvH+LZU/9wMHMvNrABHxJWAd8E3g1RHx05Xj1gE/WBmv99XMvL3JdX4T+BzwMPDqOWKu+mxm1n4fmsXxCOXkFspJ149Xnt8GfCAiPgjsrj95RJxGOQH8TGXoTyhP6626qe4jt2ZmRsR+4P9l5v7KeQ5QTszuBF5cqTauppw0P5lyMtnMjwB/kZnfqZxrN/Bsyvfly5l5Z83XtaHFeSRJ8+DUTUnqXccqj/n678CnsryW7BLKFbWWMnMC+GZEPBV4CScmEFBO4N6TmZuA/9LOeWscZfbPrNrPPpyt1+VVE91jNc+rr1dXktyLgWdl5tMor2lsFtt3WlzneylXrL6nxeebnm+OOKYzs7qofoZykkVm/jJwBeWkcF9EfG+b1z3h+hVzfa/OpvwLgYsy86nAXzG/+1iv9hrHvy5J0uKZ6EmS6p0GTFSev7Jm/N8oJzHN3AS8ATgtMxtVeGrP+4om5/g74GUAEfFEYD3lKahfAc6NiFURsQ545pxfRftOAx7MzCMR8STg/Jr3piOinYomwP8CfgO4AfjtBu/P9f1rFUdDEfGEzNybmW8FDlNO+I7LzIeAByPi2ZWhnwc+w8KdSjk5fKiyTvD5Ne81+/r+DtgaEWsi4hTgpytjkqQlZKInSar3duDaiBhndoXlU8CTo0EzloqbgZdSnsbZyFXAhyJiH9CsM+MfAqsqUwdvAl5ZmXp6G+Wpnp8Hfg+4Y35fUksfo1ytupdyw5TaqZm7gLsj4oZWJ4iI/0S56vZnlXM8IyKeW3tMZn4TuK3SkGTnPONoZmeUG9fcA/wDcFeDY15ROe5u4FzKU0wXJDPvolxp/ALwZ5TvS9Uu4GPVZiw1n7mD8vq/zwJ7gfdl5vhCY5AktcftFSRJkiSpx1jRkyRJkqQeY6InSZIkST3GRE+SJEmSeoyJniRJkiT1GBM9SZIkSeoxJnqSJEmS1GNM9CRJkiSpx5joSZIkSVKP+f8TafWXqYbV0QAAAABJRU5ErkJggg==\n",
"text/plain": [
"