{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Introduction to Linear Regression\n",
"\n",
"Downloaded and adapted from [NbViewer](http://nbviewer.ipython.org/github/justmarkham/DAT4/blob/master/notebooks/08_linear_regression.ipynb)\n",
"\n",
"*Adapted from Chapter 3 of [An Introduction to Statistical Learning](http://www-bcf.usc.edu/~gareth/ISL/)*\n",
"\n",
"## Libraries\n",
"\n",
"[statsmodels](http://statsmodels.sourceforge.net/) : basic tools in statistics\n",
"\n",
"[pandas](http://pandas.pydata.org) : tools for real database analysis\n",
"\n",
"[matplotlib](http://matplotlib.org) : best visualization library\n",
"\n",
"[scikit-learn](http://scikit-learn.org/stable/) : most used library in machine learning"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# imports\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# this allows plots to appear directly in the notebook\n",
"%matplotlib inline\n",
"#figsize = (16,8)\n",
"figsize = (8,5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example: Advertising Data\n",
"\n",
"Let's take a look at some data, ask some questions about that data, and then use linear regression to answer those questions!"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"
\n",
" \n",
" \n",
" | \n",
" TV | \n",
" Radio | \n",
" Newspaper | \n",
" Sales | \n",
"
\n",
" \n",
" \n",
" \n",
" 1 | \n",
" 230.1 | \n",
" 37.8 | \n",
" 69.2 | \n",
" 22.1 | \n",
"
\n",
" \n",
" 2 | \n",
" 44.5 | \n",
" 39.3 | \n",
" 45.1 | \n",
" 10.4 | \n",
"
\n",
" \n",
" 3 | \n",
" 17.2 | \n",
" 45.9 | \n",
" 69.3 | \n",
" 9.3 | \n",
"
\n",
" \n",
" 4 | \n",
" 151.5 | \n",
" 41.3 | \n",
" 58.5 | \n",
" 18.5 | \n",
"
\n",
" \n",
" 5 | \n",
" 180.8 | \n",
" 10.8 | \n",
" 58.4 | \n",
" 12.9 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" TV Radio Newspaper Sales\n",
"1 230.1 37.8 69.2 22.1\n",
"2 44.5 39.3 45.1 10.4\n",
"3 17.2 45.9 69.3 9.3\n",
"4 151.5 41.3 58.5 18.5\n",
"5 180.8 10.8 58.4 12.9"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# read data into a DataFrame\n",
"#data = pd.read_csv('http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv', index_col=0)\n",
"data = pd.read_csv('Advertising.csv', index_col=0)\n",
"data.head()\n",
"#data[199:200]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What are the **features**?\n",
"- **TV**: advertising dollars spent on TV for a single product in a given market (in thousands of dollars)\n",
"- **Radio**: advertising dollars spent on Radio\n",
"- **Newspaper**: advertising dollars spent on Newspaper\n",
"\n",
"What is the **response**?\n",
"- **Sales**: sales of a single product in a given market (in thousands of widgets)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(200, 4)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# print the shape of the DataFrame\n",
"data.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are 200 **observations**, and thus 200 markets in the dataset."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAAFHCAYAAAC1VKUoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvX+UHNV17/s9zTDjQRqQWiMLuBiJDHYAw5Oa4cny1cUz\nedGoncRRIrSuYhyyRlauFb+VQJAaNIgxWMuauWJkD9azbxIix44UEuwQE3kpP263Br+ZuVFyY11j\nScYgYnCAZ8cGI4/t2M/zGOPZ74+q6q6uPlVdP05Vneren7VqSdPddepUnW+dX/ucvQURgWEYhmGY\n7JJLOwMMwzAMw0SDG3OGYRiGyTjcmDMMwzBMxuHGnGEYhmEyDjfmDMMwDJNxuDFnGIZhmIwTW2Mu\nhHiTEOJLQoizQohnhRAHzc/zQogpIcTXhRAnhRDL4soDwzAMw7QDIs595kKIS4joJ0KIDgCnANwD\nYAuAC0R0SAgxAmA5Ed0XWyYYhmEYpsWJdZqdiH5i/rcTwEUAvg+jMT9mfn4MwK/HmQeGYRiGaXVi\nbcyFEDkhxFkArwKYJqJnAKwiolfNn7wKYFWceWAYhmGYVqcjzsSJaBHAOiHEZQAqQohfcHxPQgj2\nJ8swDMMwEYi1Mbcgoh8KIf4OQD+AV4UQlxPRK0KIKwB81/l7buAZJ0Qk4r4G646xk4TmANYdU09Y\n3cW5mr3XWqkuhOgGMATgDIATAIbNnw0D+ILsfCKKdHz4wx9O9XzOg7p7SJKoeVV536rT0jFPuqaV\nNDrdezukpWOeououzpH5FQCOCSFyMDoNjxLRF4UQZwA8LoT4bQAvAdgex8VfeumlVM/nPKg5P6uo\nvG9VaemYJ53Tyhq6Pkcd09IxT1GJrTEnoqcB3Cz5fA7ApriuyzAMwzDtRst6gNuxY0eq53Me1Jyf\nVVTet6q0dMyTzmllDV2fo45p6ZinqMTqNCYsQgjSMV9MOgghQAktgGPdMUBymjOvxbpjAETTXcuO\nzGdmZlI9n/Og5vysovK+VaWlY550Titr6PocdUxLxzxFpWUbc4ZhGIZpF3iandEenmZnkoan2Zk0\n4Gl2hmEYhmljWrYx18HWy3nQx56UNDra5HTMk85pZQ1dn6OOaemYp6i0bGPOMAzDMO0C28wZ7WGb\nOZM0bDNn0oBt5gzDMAzTxrRsY66DrZfzoI89KWl0tMnpmCed08oauj5HHdPSMU9RadnGnGEYhmHa\nBbaZM9rDNnMmadhmzqQB28wZhmEYpo1p2cZcB1sv50Efe1LS6GiT0zFPOqeVNXR9jjqmpWOeotKy\njTnDMAzDtAtsM2e0h23mTNKwzZxJA7aZMwzDMEwb07KNuQ62Xs6DPvakpNHRJqdjnnROK2vo+hx1\nTEvHPEWlZRtzhmEYhmkX2GbOaA/bzJmkYZs5kwZsM2dSoVKpYPPmbdi8eRsqlUra2WHaDNZfNuBy\nSggi0u4wshWN6enpVM9v9TyUy2Xq7l5FwFECjlJ39yoql8vKr09EZOohE7qzUHHfqtPSMU9h03LT\nn6p8JaU5Uqi7tMtERrlcps7O5U3riaTzpeOzIoqmOx6ZM6GYnDyC+fkJAMMAhjE/P4HJySNpZ4tp\nE1h/2WBy8ggWFnaByyl+2GbOhGLz5m2YmtoC4yUFgGMYGjqBkyefUH4ttpkzTuLWH9vM1ZBkPdEK\nRNFdh+rMMO1BqbQLp04NY37e+Lu7ewSl0rF0M8W0Day/bMDllBwtO82uw/7oVs5DsVjE8eNGL3to\n6ASOHz+GYrGo/PpZRcd9rDrmKWxabvprV70B6ZeJjGKxiP379zStJ5LOl47PKio8MmdCUywWI72Y\nDBMF1l82WL9+Pfbu3Zt2NloetpkzSqhUKtWFLaXSLqWVLNvM24c4dRQEtplnF100FIYouuPGnIlM\npVLB1q3D5upiwy4WdTrNDjfm7UHcOgoCN+bZRCcNhYGdxkhoZXu1bnnw2iakiz0paXS0yemYJ3ta\nKrabtaveAP3LN4m0/GpI1/uLQss25gzDMAzTLvA0OxMZnmZnVKDTFClPs2cTnTQUBraZM6nDC+AY\nFeiyeIkb8+yii4bCEEl3Yf3AxnmAfbO3TB7YN7s+aemYJ13TSkpz1OK+2XVNS8c8EbFvdiZlOCoS\nowLWERME1ouDsL2AOA8oHCExwSiXyzQ0dBsNDd3WEN1I9p0zelUut5wKhY2RIiM5QQZH5kwwjOha\nK6s66uhYQYXCRqkOnee56TUKSWmOYtJdXM8lCm55CpNXv1Ebs0YU3aXecEszxZVqKni9IG7fDQ3d\nZn5G5nGUgA1KXy5uzFufQmFAoqMbPSvqOCv0LDfmOjZ0bnkKm1dZvTM0dFsCdxIvWjbmAN4CYBrA\nMwC+BuAu8/P9AL4F4Ix5vFtybuSHooOtN2t5kL0g/f23un5n9aYbK+Hbqt+zzVyftHTMk5VWPt8n\n0VGfZ0XtpdeoZLExt8pERUOnunyD1B9eeVV1j7q+C1F0F6dv9p8C2E1EZ4UQSwE8JYSYAkAAHiai\nh2O8NpMQzqhIwAiAYwBeSS9TTOZYvfpyzM3dY/vkHgA/n1Z2GM3haGwSwvYCgh4AvgBgE4APAyg1\n+a2yng7jnzDT7NZ3hcIA5XIrCChldsqTdZcehs18GQEbzOOSplriaXY57TDNbqWp27qAqETRXSL7\nzIUQawDMAng7gBKA9wP4IYAvmw37Dxy/pyTyxTTitUez2f7NuPZ38j7z9sCun4GBmzE7+xUA3lrK\nuubMaynXnY57rd3ypGNe0yKK7mIPgWpOsX8ewO8T0Y+FEH8E4CPm1wcATAL4bed5O3bswJo1awAA\ny5Ytw7p16zA4OAig5gvX6++zZ8/i7rvv9v171edbDA4Opna+/Vzr70qlgtHRcQDA+PhoXQzowcFB\nFItFdHV1Sc+3Qk46fRFbf8u+d17fT/4PHz6Ms2fPVss/SaLqzvo7zH3HqQPAeK5h7yeJ++vq6sL9\n999Z/X7jxubnnzt3rprGuXPncP78+VDv7czMDI4ePQoAmdUdALz++uvS93tmZibx+s+pO7f6w1nu\nSegujvsL+x4p013YIb2fA8DFACoA7nb5fg2ApyWfR56uyNrisyTyEGZKS4d7QAan2XVcYKNjnqKk\nJdPzxMSEkjwlpTlSqLuJiQll0+s6lG+caemYJyJNp9mFEALGSqjvEdFu2+dXENF3zP/vBvC/E9H7\nHOdSXPlqZzZv3oapqS0wIgoBwDEMDZ3AyZNPpJmtpvA0OyMjTj1ncZo9q+83U0PXEKgbAdwB4BeE\nEGfM45cATAghviqEOAdgAMBuz1QYrWEvTEwYWDeMF6yPEIQd0sd5gKfZY8mD6ml2P+nxNLs+aemS\nJ6duOjuXh5oO5mn2elplmt2rXmm1d8FJFN2xb/Y2olgs4vhxY+ptaOhEoNCAsp7y5OQRM9TgMAAj\n7KC1KpVh3HDqZmFhV1U3QUZkMj2vX78+/hvQlPXr14d+v3Vi374DmJ+/BsAJAJdzveKXsL2AOA/w\nfl+tCOLKNQ6XisjgyJxxx003Ou2PTkpzxLqro1wuUy63vKoBYBUBpZZw1eqHKLqLfWsak33qR1LA\n/LzxGXthYsLgphs3nWVxdMmEY3LyCBYXP47aIj4glyuhVPqL9DKVEVp2mt25lzHp83XKQ5TFJF55\ncJu2t1/v0KFDEe8gm6goO9VppZEnmfacutm/f0+kBpv1ZqCL5pxlHjVfa9feWNWHbu9CpVLBLbe8\nS4+FemGH9HEe4AVwkc63uzncuXNnpKnL6enpQNOfqhY32UEGp9l1XGCTdJ786sZKK8w0+9jYWN20\nrAq9EWVzml0HzUVdlOjl9nVo6Dbq779VSfmqeFa1vI4oMwtF0V3qDbc0U2xDCk1jfPEVSuzafv0g\nx2FHz2JjzoTTQhB/24Z9VY2+nWSxMdcBFe+/UwM6raWwo1tdxzbzFsNpd1xcfERJupYrRoaJkyA6\nM+yrb405R0zSODWwefM2XkvhA7aZx3S+LnkArkYutxuGM75j5mKjXbHloVTahe5uKwzqMXR27g50\nvVZBF/tlHOn4TcupBTftRcvXRtTC7h6DEL/XlnoD9NCcrMw3bXqHsnwB4fLVkIqCZ1W71/sQpm5V\nDY/MWwznSuHOzik8+GAJs7MnzO/j3XtqLW6y9oVu2nQf96DbFKcWVGuvpvU7ADyCXO557NixnfWW\nIrIytwI3haW+TjuP7u6jWuyase51dHQc+fzzsdetzUgkBGpQ2Ed2NMKGFNQ1FCH7Zm99dNNsFn2z\nR0XX9x/QO28qiaI7bswZAMbLsnXrsGmbMvb+6uJBihvz1kZH7bVbY65jGbQjugZaSRUd7NVZyoOb\na1Yd7iGL6GC/jCsd1WmNjo4rcwvcrnoDot174/u/Q5kLVR11p2OeotKyjTnDMAzDtAs8zc4A0Hua\njafZWxsdtcfT7OmXQTvCNnNGCbouMuHGvPXRTXvt1pgD+pVBO8I2cwk62HrTzoPTb3AzH+3FYhEn\nTz6BkyefUOYLWRd7UtLoaJPTMU/N0goaV6Bd9QZEv/disVjdJz06Oh7J13hc/vJb/V2IAu8zbzGs\n3vWFC6/imWe+joWFnQCux+zsbwH4KRYWDgMATp0a5mk0RgtOnz6N/fsfrk7xWtoEUDf1a/+cR5Dq\nqZ9qP4+tW8PVEc4p+9nZ3Vi7dm3TdHhmICJh/cDGeaCFfBUnSb0P4w0NfoOBGwm4zTyyEyMY7Ju9\npZH5uC4UNlI+3yf5fCARP91JaY400l0zX+NxxmfQ1f960kTRXctOs7cj9dtLrpT84tsAtpjHMVy4\n8KprWlHCpjKMRTgdPY1z557F3NzKhm9efvlbyraxMf6xRttTU1swNbUFW7cOK60X3LbGhs1rW9Zd\nYXsBcR7gEKihzq/vEZcJ6K2G5xNiGQGlut5yX986aU+7vpc8EqmXrOI5IoMjcx3CUcaVjt+0GsPh\nrqRCYWOD3iYmJlwi/ZUJqB+tFQobPUd9qu4xKc2RQt1FvXev9z7IaDtMGGS/6bvdozVrUChspM7O\nlU1H+Dq+n0TRdMc28xai0S/7G7j66jKuueZ5XLjwdpw5c1Pd71988Zv4xjfuBlBvQ6/vJc9gfv56\njlLEBMYZwW9hAThz5hEAW+r0tn79+jp/3hcu3IgzZwCgCCNgx37k86/hscfsdnTjGkZwi/T9dLcC\ndr/qc3OvYXw83JqaMPEZnHVXkHKtt9E/AuCjaMsIa2F7AXEe0MSGlEXc7FqNcc6XN4zUrZ5wHHF6\no4AMjswZuY6M9RrBRnbO0VWQmOdhSUpzlBHdJWHTDluu9TrTq+4KShTd8T7zNsK+WvTChVdx5swH\nYPVggWMoFD6F3t5VuHDhe3jmmXPVle9pO5DgfebZxLmqGbgLQAHACgDXYGjoRZw8+YTruWmubG7H\nfeZOnGUA6LmLYPPmbZia2gKjLqsAuAPAxwCkX3cFJZLuwvYC4jzANnPleXD2ehvtWsvqbE2WfbO/\n/9ZIPXC2meuTVhp5snTX13cTAZdW9QVcSmNjYw1pRR11s81cDc51DF4j8WZlFrfuZHVZoTDgqSEd\n30+iaLpjm3kb4BwhWfbKejvlWpw5837Y7Zu9vSdw//13YnBwMKWcM1mnWCyiWCxi8+Zt+MY3SqjN\nBAGzsycwOlr7rZtOszKqaiUef/xv69Y7uNmedSizxhjqn2tPzYTtBcR5IAM2pCzhZgO396gLhQFt\nbU3I4MicqcePvnRaq5GU5khT3fkpi3K5LPUFoEu9kUWi6I5H5m3KhQvfq+tRd3bejc7Oe7GwYHwf\ndpVw2rZORj1Ry7RSqeCZZ84BuKf6WWfnvSiVHlWZTUYhzVaX10bk1/hKj+uFBAjbC4jzANvMleZB\nthJVtl/X2gNstzUFyYPsOhMTE5HvARkcmetokwuTjtsq5iBp1UZ5ZXO18QYqFDY25EvFimm2mavB\nqjfcbOH1ZepdZkHs737ypQId30+iaLpjD3AZx4+3I8umNDR0AkNDJ3D8+DH09q5q+F1v76qGQCtB\nkHlxevzxvw2cDqMPKj1zGfvGnwDwQan+ZDq1dNi2Xr1SRBZ4SfIrGL4AHkE+f0BqL6+3vwfTEJd7\nAML2AuI8oKENSUeijGTi2Dcal80TGRyZtwoqyjSq1tLw252U5iijugtSJmE11I7+2qPoLvWGW5qp\nDIo7DaJWtKqdb8T18nFjnh6qyjSK1tJYGMeNeXP8lmlYDem0IDIpouiuZafZdYjDrXse/EyjuZ0v\nm/6STZN2dXVFvYVMomO85DDpWGVaKHwK+fwBXHfddaHS8tKajs8qiyT9HP1NwwNdXV2u5hN/VABs\nA/AIXnzx6wHOc6cVNdeyjXk7UCrtQnf3CIxVwu9ELlfCwMDNodIKYpuKO4ISow+VSgX79h0wo5g9\ngDNn3o+tW4dx+vRp19+rsHHa0xkYuNnU+TEAx8yV1btCp92qPProo1ix4lqsWHEtduzYoZWt2W/D\nb6dU2oXOzrtheHTbAuCDeOmlb2pxP1oSdkgf54GMTjulwdjYmOlnPfw0aNBpMK996zzN3jrUynOD\nr+lOlVPyznTGxsZi98duJynNkSLdjY2NkdPDnhF7Idu2Zp39X8RBFN3xPvOMMzv7FSwufhxRogQ5\no1uFjTSkKh1GD2rleSLg79XraHb2hKsfdwZ4+OE/BfAJ2D3sGeX2sUy/h729K9LOQmZo2Wl23e3V\nWc5DbXq/+bTn3Nxrka6fVXS0yYVPZxeA+vLetOkdSvIULV/xppU9zitLSZcycdY1nZ27lZhYdLk/\nlfDIPONEiQMcNo1GX8i1BS3OdLZv3xPshhhtqOliAsAdyOVKWLv2Rhw8KF/YqEKLKtNpJ/bseT8+\n9KH/CuB685O7AHwAtc52Np9fmNjo7UpsIVCFEG8B8GcA3gyAABwhok8IIfIA/hLAagAvAdhORD9w\nnEtx5asVUeEq0Z7GwMDNmJ39Sl16fq8Rh9tGDoGaHlZ5XrjwKoy+/xsAOtDbu0JavqrKP233n1kM\ngTo+Pm5OtwO/+qv/Cd/+9o8ABHtfAT3DnLYLWoZABXA5gHXm/5cC+BcY3cZDAPaan48AeEhyrsIl\nBYydZntD3RYfpem8AbwALlVqmigR0CvVgWqfBWmTlOYoRd01C4Ochk+BdieK7hIRq5FHfAHAJgDP\nAVhFtQb/OclvIz8UXfyi65QHP6uNG1eqj0SKjMTxzPVJK2w6NU3YtTGtbBeDjs8qi4150HtvfNc3\nSMs3CHHFZ7BI+12IO60ouktkAZwQYg2AAoAvmQ35q+ZXrwJodNLMxIJaP9sMY8C6Yiw4PkOKhO0F\n+D1gTLE/BeDXzb+/7/h+TnKOsp4OU8NvjOL6qbeV1Nd3E+VyK1Lbt4oMjsyzjn2qtGZmkU+zt6Lb\nzaQ0RzHpzs9UdxzT7K2ohSSJortYV7MLIS6GESbpUSL6gvnxq0KIy4noFSHEFQC+Kzt3x44dWLNm\nDQBg2bJlWLduHQYHBwHUtgLw38H+rq0SNrawdHZ+BhcuvA233PIubN/+HuzduxfFYhH79+/B449/\nGouLOTzzzE/xjW/8MgAglzuCtWtvwHvfu6duNbPq/B4+fBhnz56tln+SsO6A119/3YxVvQMAcOrU\nJzE6eieOH6/ghz+8HD09f4re3hXYtGkPzp07hwsXXkUutxuLi4auuruPolQ6ps39+Pl7ZmYGR48e\nBYDM6+7QoUN44IGHsLDwcQDA7OztOHDgPuzdu7fu99ZK8dHRcQDA+PjnAMD2t7FLJcj1S6VdmJ29\nHQsL5wFcb25j3IOZmZlUdDw5eQRzc69V67ckr5+47sL2ApodAASM1ewfd3x+CMCI+f/7ENMCON3s\n1SrOD7OwxJnG2NgY5fN91NNzNXV0LPFhPx+J1Mtmm7k+aflNx8/oyop3XRvZlUiI5dTTczUVChtD\n2cxVLJximzlRobDRtH/fRka88WTfW2c5OtftRCljv/lqto5Dx/eTSN+R+UYYTnW/KoQ4Y362D8BD\nAB4XQvw2zK1pMeahZbD8oRv2KGM/d9CABePj43jwwUnTYxxg+HS/HEAx016imHSot49WQHQMP/rR\nR3DmDLBly2/hxIlHfetJhb4Z4zmeO/csAOsdH4ZRDSdHsVisKzdrJBq2jO3b5zZtekd1hOtFW3qj\nDNsLiPMA2y4bUBHu1LB716dh9N7l6ekSTxgZHJlnGb/lXq/JRn0WCgO+r6mbrTUpzZFi3cmeYy63\nQostYmHKuN3Cp0bRXcu6c80ScUSacqYzOXkEi4tvlZz1bbi5ZJWFNG3pni0DwH+517va/Lbj26fx\n9NPnlUbtUvWetBtr196Y2ffWGGHfAcPP/AnMz9/ha6eEyoiSmSFsLyDOA21kM/fqeQbZI+72vZWG\n0VMtEVD7jRDLqFDY6Gm/0uE5IoMjcx1tcnHkybKBLl16BQH5qv3cHsGr2Wiq0f7uX99x3WNSmiOF\nuvPzHIOkpQq7VoLmzbD/13ZQAJdSobDR13W9Ikrq+H4SRdMd+2ZPGb+2HWu0tG/fQbz88rewevW1\nTdPZt+8gcrlF5PMrMTBwM06d+qTZy30Eudzz+MhH7sHo6GgCd8m0KpZO9+17FWfPngPRIwC+BXsE\nL5mmnXbQvXv3uvr7B9rUBhoCr7gJSeF0EWvtfCkWixgdvRMPP3wAALBnz50+8tYB4GOoRYM7D+Cf\nfeVDRUTJTBG2FxDngTayXQax7Xj1bN1sZX5jQjtXmerkkhEZHJm3AtbOh3y+j8bGxqS/ca5oz+VW\nUE/P1Z6aDjNCS9oGmpTmKGHdBXmvw9QBXmWbdLln0W4eRXepN9zSTLVRpRpE4F7idKZjTC+VQnUS\nOjtXUmfnsshTdargxjx5xsbG6qbKgUulDbpMk4XCRk9NJ7kQKiyt2JgHeYZxLDxLutx1WcAbhCi6\na9kFcNZ2iLTO95uG12KjIHlwprN27Q0AbgLQPA2nC8aFhY9iYeE68+/VkdxzqniOWUTlfatKK0g6\nRvQta6p8GMAnqhG5mqXV27sq4MLJ5nG4/S7Ka1e9Ac3vPYjb3dHRcWUueufmXgt1HtBY7vv37/E9\nTa6qbm2GLppjm7kGOPdlutEszrM9ndqezh0AXs50TGNGX9w06aVp5zmdnUdQKn226bX8vidMesj0\nsH37Htfv/NRJ9nIP2nC2lWbCDunjPMDTna72qjhsXjzNzrpz4neanSi8bdXu+12X9RkWSWmOWmya\n3TrX79qcKGm1IlF0l3rDLc1Um1eqqm09foMu8AK49tadE+cCuDg0oatds9Uac6vsCoWNVCgMxLYA\nTiW6aiNOuDGXoMP+6Ojxo6dDrcK0v4Q7d+5MNdY07zPXJ62o/rbr/bBfSoXCQOTKVYX/fzvtvs9c\nRphGMepztNdBYeOZyxbM9fffGilfFjq+n0TRdMc28xbD6f9YiN8D0e+gbfZaMrFQWzx1OYAREP2f\nOHPmemzdyj7UdSfpPfrOOmh2djfWrl3LGombsL2AOA+0yXSnl13c2ZP2a1eU9WaNCEq1v/P5Pi2m\nz/2CDI7MWwVLo/l8n6mlcNuLvGyozpjafqeB4yQpzVHMuiuXy2bZbSAjgpq/MvNKr1k9FHZ/t8zU\nx9PsAXQU9sQ4j3aoVP24r7QvEIqyF93uPMZY1FTK1MvBjXk6OKfWDe1sCFRR+6mQ7fbczs6VWlTe\nrdCYO5+94Ra1FPq5+m1cVe4nT9tunzTcmEvQwdbrlYbfmNF+f2vRONJZXh3VGz10f45k/NyDH9hm\nrk9aQdNp1F2Jli69wnRINOKrwY1L56ru0Y0sNubOe5c9z3y+z1ejKHuOfstHVgep8vSmw/qRONOK\noruWdRrTrjgdJRw4cB9GR0dx8uQT6O9fC8ORDMOE4Sa8853vxN///WfR3/9PHEkvg/T3x2+7ltVB\nrJEECNsLiPNAG0x3JrXnM0o6ukxxIYMjc10J6qcgqu7S0LkKktIcxai7OLa4GumVCNhAudwKV98D\naec1q0TRXeoNtzRTbVCpEsUf9CBKOjq9XNyYqyFMmarQXRo6j0orNOZE6p+nV1jRqOhS9mnCjbkE\nHWy9OuUh6Isi2+turTIO8rKxzTz9tKyy7+m5UmqXDFOJ6nR/caSVtca8XC5Tf/+tyjzquT3HMOsa\ndCxfHfNEFE13bDPXkEqlgs2bt+Geex5EpVIJfN7mzdvqzrP2fU5NbcHU1BZs3TocKF2Lc+e+FjkN\nJlnsZf+jH13W8P2FC9/zpQ03bTW7dtBzmOBYZfzUU/8RU1Nb8KEPHcLU1DWYmroGv/zLv4mbbx5s\n+vy5rFqAsL2AOA+0+HSnHVV7K73OU7FVJEhIVdUggyNzXagv+zIZ25NqGikUNjbVRtjp+Syuz7BI\nSnOkQHdy3xIbCVC/ViGq+U23ctaNKLpLveGWZqoFK1UZshejUBjwvf3D/lJ4Ndi178pkOP3YQIXC\nRl/5s67hp9KPC27MwyPbXmZ3GmSU6wZTF2VpufrpDAbRoxOd1mdYZL8x7/P9/IN29v00yLLfZL2D\nlwTcmEvQyV7thts+0Npn075HSm6dgOnpaSqXy9TRsYSAWiS0zs6Vgfabhq1w2Waeblr15TZSV27l\ncrnOSQuwjIToaVih3MxHdpROKRFRf/+tyjqK7Wgzrz3/EfM5XkrAjZ7PtL6j3lhWUXygu+0zV9HB\nY5u5+8G+2TVj9erLMT8/Ysb8PY/u7qMYGLgTmzdvA2DEBJb5WgY+he7uEddYwT/7WReAh6vnLCwE\n889s7R2dnDxi5oP3F2cBe7nNzb2G8fFauU1OHsHCwkdhaQIAiB7B+Pgnccstt1R/5xajulKpYHLy\nCJ566hzm5+9AED0y6rDKeHR0HPn88xgY2IsnnvjvOHduNxYXjd/Yn7/Td3pn593o7LwXCwu131ox\nyMPgrJ8WFs5X6w0vLD0BxloOmT/5+++/M3S+Wp6wvYA4D7TQdKfXVJEfF4YyV65uU95u1zJ6xMHc\ncOoEMjgyzwLy6dnbfE2jy12Fln3p0YmO069JaY5i1J13fVBf7oXCRmXP1G0E7lXOzi1vNRfU2auv\nohBFd6nXkMLSAAAgAElEQVQ33NJMtUil6qeSalYxyV+8gaZ+3QuFAcrn+6rxi43Fa/UL2rJih+LG\nXD3lcpn6+m4g4DJbg7yK3OzmTuQdgQ2BTDDO/Ojk/6AVGnM3wq6hceJWZl5l5GZLb2y8S7HtZ9cZ\nbswl6GAzj2oLnJ6e9uzlur1I9XbQS6mjY4n5WXDPTTo8xyw25jra5OzrHzo7l5EVeMNohJcRsM13\nxSnTdthofEHur5ndtR1t5hZ+772+/OVraJql1axTZa+fmsUzd5s5lM0WsM3c/WCbuebI7JWWvVpm\ns260g57HG2/MolC4CL29LwK4EqXSfrZ3tzGGRq4D8EHUdHIM+fwB9PeTr/UQ27e/B88+W28Tf+wx\nXkeRBYrFIt7+9rU4c+b9CLuGplmMdHv9NDMz4yPFjQBGqn/lcrtx8OBnWU9BCNsLiPNAi0x3qpoS\nDGIjdJv+zLK9CRkcmeuMqjUUaWwd4ml2NYTxPaHyfDu1MlXv8z1rRNFd6g23NFMtVKkmXeE1TrP3\nUmfnskzbm7gxV4ufaVad4QVw0VHh/EV1EJd221MugxtzCTrYetPKg30B3LXX3hjp5fBjO/N6Cdlm\nHk9aYSo/ezqWwxhjkeRArHbuLKaVxcY86L17achPWm6L2eKyc6tMS8c8EUXTXSCbuRDiIgBLiOjf\no0zttxP2vZOl0q5EbEBOe9Xg4GAs13HuVz11ali7+NatqFkVz91tzUXSVCoVc3/0ysTejyTIgu6i\nasB5vpsuu7q6IufVThp1aiZo1toD+CyASwEsAfAsgH8DsDds78HPgRaZ7vSzjzypqXe/236C5Kve\nblYmYAPl833K7wkBe6thNZsV3RnPvWRuKzL+32xfd5rTmGG2MKVNUM2RRrrzet5BNBCtPqBIdnSv\nPLWyb/gwurMOPwI9Z/77mwAmAVwM4OmwF/SVqYxUqs1wc86QZAXmV/xhXpL6/arx3VOIxjyUZrOi\nO8NpUM3eDfRSX99Nnnt702o0va6dROUflpCNeeq68xpABHXME74+iK88o1xD586jRdyN+TOmKP8K\nwKD52VfDXtBXplrEZu62F9evGJPc6+4WcMMrD7WXw31ldBo287CaVVmpxmmTk/nSvuiiFWY5lBvK\noFYBTiurZP3en1flm2a+mhGyMU9Vd15+KYI2grJ6o5mXOLfGUuW7EMV3R/0zmFbW2dDFZu4nnvkf\nA3gJwFIA/0MIsQbAD32c1/Zs3/4edHePADgG4Bi6u0ewevVVsVwrSjziSqWCc+eehbHveAuMvaNP\nNz3P8gmdz78WJstx0tKa7e1d0fDZz372VtT2jbuV/2kA2wA8ggsXXo0tf34plXaZ70cZ1vtRKu1K\nO1tR0FZ3Rnk/AqP8w8Qrfxrnzj3rGffeqg+Ghk5gaOhELOtnZHVqxjWjjqCtPwABoCNs78HnNZT1\ndNLGaaOJw47u9Gvsd4rVuqZstiCXW+E7H3FPXyHiymK/ms2K7pzP2+kXHdhQt9VMthWto8O7fFXZ\nFoN4CtNpyjOq5igl3cme99jYWODtqkadYs32lEiIvLIp9KhlHvZ8nmYHLgfwaQBl8+8bAPx22Av6\nylRGKtWw+G3g/abVLCiB2xaS2jVlrhQHIt2TSoIKPKxms6S7+o5Yqa7sgKsaKuu+vnW+y5j3EIer\nVHXRnfN5y9fuDHieby9/IZYRsFRJY552g6q7FuNuzMsAfgOm7QeGTehrYS/oK1MZtZk7heJ3f7Cb\nPctPHtz9Gg943kfjSvRe6Qumw9qDEI15KM2qrFST2sfaOEq3gqWU6nylG43+SJ1G8vk+aZp+7KvN\n8hSkwnSmFaXCTdlmnqru/L3rtbJ3e77G7+u1YsRHDxeoyZ4vt0XBfsvb6SshbZ2oTitKY+5nn3kv\nEf2lEOI+U3U/FUK84eM8CCE+A+BXAHyXiG4yP9sP4L8AsAyt+4io7Cc9nZHtsdy/f0/dHm+3fZjR\nqfdrDNyNZ54xrudvj2cRwLDpm3ttK8QqD63ZrGHZKd/3vt/F3NxKGLZEADiGubmPYWrK0NmVV74Z\nc3N/BOB68/t7sHr1zyvPT9Q98FnwXeCBlrpzxncA7sLc3AcwNXVTgOe7FMDHABwB8G2sXXuDojIx\nbPGLix8A4L+8M66TeGjW2gOYAbACwBnz7w0AZv30FADcCqAA2/YMAB8GsKfJecp6OknhZzTj9hvn\n6CqXW+7bN3HtXCv61QoCxppOg6mc7op76grBR+ahNKuj7oL4CPA2m2w07eYbyLCpu9tMo6zriLo9\nSZftakE1Rxrqzl5eY2NjrmYZ5/N1ln9n50pTO2rqivq6Llzccl10opowurMOP6vZSwD+BsDPCSH+\nCcCjAO7y2VH4BwDfl3wl/JzfLhSLRYyO3olcrgTgESwu7sT4+CddV6XbV64DwPHjx7B06WMAvgXg\nRgC3+LqmipWnVg/Za5VrCoTWrE4Eebb28jR2FzwNY+XyNgBPo7d3FU6c+ByGhq7E0NCVePDBezA5\neUS6+0GmDQA6lrNuaKM7p3bGxz+JUmmXuZvmprrfXrjwvbq/neX/4IO/j7e/fS3y+QMoFD5VrSvC\n7KBxpr127Y2qbpnx0+LDsP3caB4XB+ktAFiDxpH5SwDOwVgsskxyTuQeTtK2XtloxhnHN6gDjf7+\nW31dx2u1qtt9+B3xNXsOsj3P9h5yWr7Zw2hWhe4sVNx37dlOe44+nGU5NjZGwKU2PVxKY2NjdfHM\nozsEGXHNS9C0nXbQKDNGaftmT1N3zWzTQ0O3SR0OFQob69Ipl8vU339rdUQf1QmNV5l47cTxSksX\nnahOK6zuiDxs5kKIbQAIxija+hcA3iaEABH9dcj+wx8B+Ij5/wMwPCX9tvNHO3bswJo1awAAy5Yt\nw7p166r2Zys+rtffZ8+eDfR7v+dbvqQBYHx8FMViETMzM+jq6sLx48cwOXkEc3OvYfv2PVi/fn3d\n+VavtHb+ser5c3P2vdozAM7X/rJd34gjvAPAagCDmJ8HJibux8LCTthjmF99dbmatvP+Dh06hAce\neAgLCx8HAMzO3o4DB+7D3r17Az2v119/HefOfc3M6wwA4/u5udfqfMIHff6HDx/G2bNnq+XvFxWa\njao7VX9XKhWcPXsGxkztajN35+t0MjMzg9OnT2P//odN2+F5TE9vQ3f3CgCfsJ33CczOnsCSJUsA\n2GNRG9/Pz09gcvJIdX2FW/78lLPsPXDGtXamb7+foOefPn0aTz75JQDApk3vQGdnZ6jnPTMzg6NH\njwJAZnVn4VWf9PauAvAfYIyjVgIYRi53GocOHcKTT34JFy58D08//WW88cYWAEV88YslLC7+F9jr\nG6v+qsUzn8H8/I5qPHNn/s6ePVv9215/bt1axPj4J7G4OARgArncdzE6uhtdXV1N648wOlPZPrjd\nX9Dzo+iuAbdWHsBRAH/qdvjtLcAxMvfzHTS0XRI5e4MlyuVWhIo41Txtd7u5rMfd0/OWhs+8tp6o\nsjcZ6ZQo7CpXv8BnbzWqZnXSnd9nK9+V4B2rPOjqZqL0txTZ82G3A8eVJ7+aI01157XuQTaz5+63\nQK6lRg2VAmsorL28lQmiO+cR6qRAF2icZr/C9v/dAB6TnBPDY4pOTcDx+CJ3OmqQOXgpFAYaFqP0\n9d1EzabO5PehojG3noexRc7rumGJIvAgh0668/ts68vSW59uGjKm5EtNtZz2Hl1Zh7fZYq6wJKU5\nilF3buXlZx+6oSWjkZZNgzsHNnazjpuGGjueV5GbC+J2JfbGHMB7AOwF8KB1+DzvswC+DWABwDcB\n7ATwZwC+CsNm/gUAqyTnRX4ocdjMa2L01xgGzYObXVLmxGHp0iuoUNhoexkbo2i55UGVvctPOina\nzANrVmWlGvW+65/tiM8AOfZR1BgBV1FHx5tpbGyMyuUydXbWKuXOzpXVeOZhGsQ0bI7yRmdD3d+y\ndSZhCFuppqm7sGUif67XV99payW80+lUoTBA+XwfLV16hWd9aOXLbTBkDERKvgZFHM/c/Wi6z1wI\n8ccAugH8HwA+BeA/A/hSs/NMhd4u+fgzfs7Vkdp+zWsSvW7NxjkMwHhdfvzjR/Dccy848mXsuTT8\nFbvvX7ds97WYwOFWsqtKRzVRNKsL9mc7N/dadX2F1+8uXLgIzzxzLxYWnoax3/xjeOMNYHx8BNdd\ndx0WFnbB0tDCAtDbewK9vaswNXVTQ7pZIZd7HouLhta7u0ewffue1PKSVd2VSrvwxS/ejsVF65N7\nAbwJ+fwBPPbYMXO3Te33zj3eudxu+InlUF9/1uozAHXXYkLSrLWHOUWOmlejpQBOhe09+Dmg0XSn\nE6NHujHQCkz7uV7TlG4jXfdpsFoPOO0p0DhB8H3moTSrs+78UvP21mgTl42e/Myu6KItN3tvHHkL\nqjnKoO7s5To8POxq4nMiq49q9u/mppogkSPbjTC6sw4/Aj1t/vvPMJZBvgnAC2Ev6CtTGahUg1Zw\nfqe2Zem6u+w0KulWa7ydhGjMQ2k2C7qT4ccGWihsbBpwR6ZlXRa+2fOTRMciZGOeuu7CORkK1jGS\n62vAd7nopimdiLsxfwDAchjeJ14B8B0AB8Je0FemNLWZB8G+V9Otgm3WG3Xuvy0UBswecMm3rSnt\n55CSb/ZQmlVZqablm122OtnSx8TEROCG0G3lu9OHQhR0tF+GbMxT1d3ExIRvz33N6qOoa2Xs+I1R\n4Qe2mbsfXvvM1wP4JhEdMP9eCsMw8hyAw27nMXab0g4A1+PUqWFcd911kdIsFotVr0uGH+5/BPDn\nAIqYn7+puseznWlHzTrXU8zPA7OzJ6RrGbq6uqq+BKIwN7cSDzzwENauXdv2mgP00d3jj/9tgxb2\n7TuI5557Tmk8CBVrZaz6jFGIWysP4AyAvPn/d8HoZW4DMAbg82F7D34OZHS60yLoNKeK9FvZ5gT/\n+8wjaTaLuotbC14mHtacXrpzm0UJu1aCSZ4gunMeXr7Zc0Q0Z/7/NwD8MRE9QUQfAvDWaF2I7BHG\nD7Gd3t5VSnyhA8aq0O7uERgrlo+Zq9d3hUqrxWg7zcatBWsUls8fAPCIeZ1kRlRR37kE0UJ3Mi2s\nXn259LeqYjPYyVB5tSZurTyAr8H0LQzgXwAM2L57Jmzvwc8BzWzmQXuxtd+PROr1qogZnbbNO0mb\neVTNqtCdRZI2uSS04HwHOjvVefqL6guhWVph8Ks50kh309PTDVpI6jl6XUdH+7SOeSKKNjL32mf+\nWQCzQogLAH4C4B8AQAjxVgA/iN6NyA4yu6SXjdrugz2ff96XTalSqdhsULua/j6qzcm63oULrwLo\nQG/vCl/X1ZxMa7ZSqWDfvgN4+eVXsHr1VTh4cJ9neQTVTBScdtJNm+6L9XpB37mU0UZ3snphdPRO\nPPzwAfz0pwt485vfXC1Dlc9SdXklqe2WwaulB/BOAFsBLLF99jYAN4ftPfg5oJnt0ssuqWKrTNL2\nq9r1SmR3A6ur3QzBRkmhNZum7gwPbcvqyqOzc6X2vtK9iPJupL0uJIjmSGPdNa558O9tLQhhy8vP\nVlwdtR0XQXVnP2JrkKMcujXmbuIKM/0uq9yahRFVjZdb2mbBEtIgisCDHGnpruZI4yrfOvBbedrd\nblruf5MgaoWcdoWelOYoZt15OZyygjGlNSAJ4iSrlRdb2uHGXIJqW2/UvZpeHQK36EEqbDFBfcxb\n/r2t/GXJZh71UFmp+r1vd9/q3jrwU+EZo317nPtLq3Hu476/IBWyW1phGpm045mHOVTpTnbvsoEC\nYHyWy61w9UkQ5jm6lZdbWm4a8dIO28zdj6a+2RmDqDZqN5sSACwu7gAwUv1tLrcbpdJnI+TWm5qP\n5DsA3GP75h7U9q4b+bv//jtjywfj1MXlAO6oftfZeS9KpUel59XK0Phb5o9/cvIIFhY+Cnuc+4WF\nWZ1tz3XwXmQVvIHGd7wXwAgWF3fg4Yf/VFovhXnvVZWXH20zEsL2AuI8oNk0uxvO0bZbDHIi973n\nxvTqBjKiXAULI+pn5OL2G+vzQmFjdRo2rnCSUUEGR+Z+adRFiYAV1NPzlmq0My9Xq17f1bRVCzEJ\nbPAsV1XuUtOeJo9KUpqjmHUni6gI9JHdHbRfc02zuBJBdeOlkVr9NECFwsbAekzK7a9qougu9YZb\nmqmMNOZE3jHI7TRu7VnmmAINtjDFz3WDVKgqKt+4XqBWacz9LPYx4kJvI2ADCZGnjo4lgcvEbdET\n0Os5za66Ac5qhUrUGo25tV5CiDzVXEDXx653m2Z3puP1myi6adYhDZNuljuS3JhLSMrWG8S+Yxeu\nzJblDJriZUu0R20zvHKVGnrURt5GfI+2ZS+W0+4f9MVjm7ncr75sFGKMkrZRfaznZQRMBJotkWmy\no+PNdO21N3pWairs3GHQMa0sNuZea3QMHV1HHR2XNYx0m733zXThZ+1QmI6dLF0/8er96FhHzRFF\n0x3bzENg3wNp7NP2h92mtHnztobv+/vdfV3XX/N7WFz8OOzxgA3vXFf6zkuz/Mmub49hbPl4tn7v\ntiag3W3uzudmrI04hvn5iart2jo2b96Gqal/gzPWM3A/gLWR8vELv/CfcP/9d2JwcDBSOn7gPcLp\n43wfDU7gjTfuQ2/vCZw8+UT1U/s7bJVbV1eXsrycPn0a+/c/7Fp3WDh1wwQkbC8gzgMaTrPbbcz2\n6fHOzpXm/mD3KSg3m3XYKXDZ6vdcbkXkaXavnnPU3nkUkMGRuYXX1iDnTIzbzgbDlNI4Ra5CW3ZU\nmVrqzUkrQ9k80yYpzVEKupPtevCKtuas8/x+75UXP3nwYwKQofM0e7N6NoruUm+4pZnSrDFvtn3I\nLZavH1uTn6kn2SIp+zS718I7v4vkak5kNlS3rHjnof6FjPMFar3G/CqybJfO52SshXCaUMpVnVmo\n0paTqHZu+f1uUK6JuMl6Y+62bkJWBvLFuQMNa3yses5qyGXfyxas+WnM3X6Tlo7jwE8dyY25BJU2\n83qh+R+B9vff6vu3XnYrt5XwfsTqjGEtu05txWv9yny7zduPEJvZ3sKSxcbc7bnVFqNdQsDVBFxF\nfX031J1bLjtXok83aEe2SrnZTEgSdkL3EaG7/nW0X2axMXdbo9NsRbiszHK5ZVTb3VImYAMtXXpF\ntUEvFDZK6yRZHeEWZ71ZHuLUShqa83OPUXTHNvPA7IJ9L7CKPZBu9mjLbiXbd3nwoD9/7w888BAW\nFj4OAJidfS+Ai829x87Yxv8Iu612cbHe5u0nhjHvC27Eem5GDPqVMPbxfxlAB4CPAAC+8Y27MD4+\njtHR0eo5jz32B6YmXgFwHt3dR+t0ZqzV+B8APmZ+cg8uXPj5xO7LDadWa74LmKTx+z7KymxxcTOM\n6GvLAXwSwB348Y9fwJkz7wdg+MIwQrbXePnlV1zXzTSrO3hvuQLC9gLiPKD1NPtR6uhYQj09V1M+\n3+c6vS07z22a0a9NKfpqULmHMTdbrXMaPa1pK2RwZO6kfvdC497efL6v4RyvZy7bDWGfhldB1CnO\nQmHAcz2JziSlOUq5vrPbvXt6rqZGvwSWi2HZ6H1FXdlGdUut49S4SniaXRP8LvaQnVcoDNDSpVdQ\nT89bqFAYCD3FFJQgjXlf3zoClpPb1ql28JMdl+7K5bK5XzxvlsHlvhpzL9ym2VVViKrKPKsVdNYb\n82BrZYwyNvaj1+vyoousuk5uV3ea8JxOtJKMBxCWJDXKC+BCoMJm7m5f9r8Xt1xuHg3LreKMeg/G\ntZfbrrusoSPiXDHqfAGnp6cjLU5hm3mZenreYjbkluOOLjIWwNUcxXjN8Mie4djYWEMaw8PDng1w\nkLJoNsry63/AD2wzV6O7IOtbiNy8D1ozKSPm/7vMTv51poa907QGL3afCp2dy5U1kqpt5io6rbrs\nM0+94ZZmSoPG3G3RRpDpzVpDKB8R2/GzeMy4/kYzAlbjCF+W1s6dOxt6z/YoWn4qbbcFeH5egnZu\nzBsXv9VWpgtxSdVUMzw87FmufgOtNHPNGST4SzOzi66VITfm7qvTLaz6QebC2Wi0B6jmibA2CAGW\nUF/fOmmnzV7nNC6OGwm08NfPPUbF61kFnRXlxjyhSjUoXkKvCdUu8F5PX+pujblzj7GffPmJd+1n\ny5KRzgYyXIYubypmWZpJhm3NYmNO5L3XF9hQnRlxlqsQy5pOT4ZpzIPlW767wU8+4tIBUXJTolls\nzC3kustToTDQMBNnNNrXEbCBOjous9UN8ql1WdnKptf9xHnw6wo7TpLWbzO4MVeEfBtRua6Qg2wJ\nskbBQvSQfYoqqB92IvdOgfPazcRpdEYurTbmwCUkRPNFSs5KNMmXIKuNuTz85FVkTF/aA100lqvV\n2Lvpw7kfPYqTDSe1si1Ts+A/SeogybUbWW7MjZkVu6+Cmj/2ml28VrZGY24MDsbGxmydQn9lK9/a\ntsKznBrzKHdJHTdprwdywo25hDBTH/WinK5WqvZC9ip8e4NnjLhqYhViqbmYpH61qJd4G/e6h2nM\n66e4li69gpwzC296U951tOPlH56n2d2Znp6WzuIYlepKH415zVuXzNxSc/JzIwF56uu7qarPoGXp\nxE/ZpjHN7qfjwNPsBob2NpCxc8I+St5AzlkXwyZerj5PY8HmZWQF5/GaCSSSl0tf3zrK5/son++j\nnTt3+jqnWUQ/5z1GQdc1H1F0x/vMm5DPv4b+/hPVPY+Tk0dw5ZUr8N3vPoiLL74Ye/bciWKx2LBX\n/Itf3I3FxSFYey6JgMsuO4C5uQ8CCL4Xu1TahdnZ92JhoRabWBbvemDgZkxN3WX75JMYGLi/+pcQ\nFwM4CLvP5osvfrDOV7Mf/Ow7b3d6e1cB2ADgAABrj3kRwE0A9qO7+0Xs2XMnPvKRj9WVK3AvgEcB\nvCJNt+Z3+3IzzYfxjW8AW7cafgNkZTk+Po6JiT/ExRd3Y8+e91f3tMsIUrasAz05ePABsz66Bobe\nLDYC+AwAZ2yHIwC24KmnzgEAxsfvx+c+9/d4/vmLIcQDuPban8PBg49K/an/678+D8Be59yFl18W\neOON/wsA8Od/vhvbt29vqotc7nmUSvvD3G4kWsY/RtheQJwHNJlml2/Pqu+tWr9x62k2WzQ2Njbm\nu1foZwFcMzNAEvuTVYMMjsyJ7JqRjbxXVO3iVrn29FxtrmFo3B5op6Y1f9OgspXvXqvndYWn2YNh\n6ap+OruXgHUudVUw8199eZQIWEHG4rnrmupSZmfPoiZVE0V3qTfc0kxpsADO2bg2q0D92o2cU/Gq\nK6dmU5HGArj6QDG67wPNamNOZDxvw7RRW5fgjCntXKDof39wc7MLEUkXxgXd164LvAAuOLWFZjfa\ntGc3/1xmNsD+zH8W7gs8/ekyyf3dWYEbcwlB7IMyQTnPrwnXskXdRk47k7Nhtm8Lc47Ay+Wysq1E\nzvupD4DQuAo5yEuk0sd9WLLYmDttcrVVwleRc6WvNcsSxNYtG3W5zfTUdDatrDFXUa7WPff336rd\nPuQsNub+/OWPkWEv76Hu7iuqftud9ZBVJl71hHyB54YGfxY67zPXLS1uzCX4ecC1UaoRLUyIfHWq\nR7boyKiQ61el28NSOsXvtkioFjbVuwcbViT2fExMTIRKI2oeVJ1PlP3GnMi55bG+zIVY7ulV0OsZ\n+pnpqU2zj1RnBqJOaUbVZr0nxe2Uy63w9J0Qd76ctEJj3rjNtkzG4sv6+kumm507d/rc4irfUqmy\nDnK7xygje27MMybuZhg9y/qVnUIscxWGX3uzn21ctVXt+myL0JUsNuZuyLfkNNoY3cLqeuFlYrG2\nHDWLJxAn9Z1aqyOr5zuQxcbcQm4rv1SqM2sFedCtp0G2L8aBblvKVMGNeUjctgW52YtkArdi+Foe\n1WS9XLfpKOP/RnjBoE5kvGg1W1TWG3NnedS2Dd1WLf/G0XptBsja/9usTHVzgOGkPn/BFvElTVYb\nc6/1FB0db3ZtzJ1pNM4glSif7/Ns7MN0QMOiu9bDwo25BD9TH0al2hhcQLa3l6jRUUdn50pzP2Zt\n2spwEFO/V925ir02ze7dq3TGIveDH5t5EHiaPRx2E4vd455sWtNpYzSc+Njt6tubOvZxs6H7na4P\ne39BqK+Ay2RM0W4gY/pfTaXc7tPsXgt1e3qubqiv7GZCInv9sZ1qU+glsu+GkDkocqvTyuUy9fff\nqqyBV+WGNe13wQ1uzCX4tZkb0axqIhTCqGy9HXVsoFxuBfX13STp6VohA8sE3ErW9JNzdNZs9Gw0\nAu4VsxvNnMYEhRvzcFj37eb+1zntbdeDsfrdXobXe1Zczi1CbvbntCswZ0dTiKV0xRU/R0IsDaxz\nlfmSkf3GvN58YW076+xcSX19N1VnEp3Punb+NFmzRrIRvXNqvq/vBnIuDC4UBszyHqnrBKhw0BJ1\nmj3td8ENbswjUC6XzUY5T8bWDfk+S1lPULaYCbjRrJyae06y46zcw/Y8W3H6KYuNOZFlH7+soZLr\n6XmLZ0XU2AGQzx5ZZKnMZf64o1bwcZDFxpxItvfbCrtb23bm5VfCvZ5z910hXwdSktaPzdy8Br1X\n3XQTFW0bcxiuhl4F8LTtszyAKQBfB3ASwDLJefE8KRf8VIZyG1F9bHNr2sqIDe6/cpU59ZCN+v1U\n0K24MCSLjXltet2+v9yo5Hp6rvYsW+fUfEfHEs/V7llqzLOS16w25kTOCGYD0ga1uQ+DmtaGh4cb\n6if7IkpZmRqzQ41b3pyOtHQs+zTRuTG/FUDB0ZgfArDX/P8IgIck50V+KEGmPgwxlsjYQ95HwI10\n7bU31v3GrZE0bJUDddNW9VNVzUUr68FedNGK0NOPvDUt/Ur12mtvpJrJpVauQuSllZxsEZK9DL1G\nIX47cH5NT1YjUChsdB31hC1XWcXf339rqLRktPs0uxO3UXOzwUVPz5XV8LzN/GG4LYar6XLEbOD9\nRatzRw4AACAASURBVFMLeo9h4Gn2cEJd42jMnwOwyvz/5QCek5wT+aEEecDGyPgSsk9t5nKXhna2\n4hRxs4ZYPl2/gVTsv027MW7HxrxcLpsdscYVxfWVnL+OWpBG2EubzdJx5svLvWeUfebOe49rH3IU\nWqUxJyLJ7gn3RrS+7rIWvnnv+PFyF2xfAKfC4yU35no15t+3/V/Y/7Z9ruzh+MHoWfrfouan4gxi\nz2l8GRpDr7YzWWvM3RYh2ddOBOkYprndxx69TRVZsHVmsTF3Yp9l8bN7hsht+6C3L4DazKa/sNC6\nl32aRNFdqlHTiIiEEJRmHoLijI526pQRqcoZdSdIJB4rgtXDDx/Av//7j/DGG8MIE1mN0Y0igGMA\n9qOn59v4q7+qRZ3yow+/WssaLROlSmOc2unsvBeFwqfQ27sqRGS7mo7z+dfw2GOy828C8DHz/8cA\nvChPics+NtJozF8VQlxORK8IIa4A8F3Zj3bs2IE1a9YAAJYtW4Z169ZhcHAQADAzMwMAnn+fPXsW\nd999t6/fb9r0DkxPfwRvvGGFoTyPXO6/YWBgH26+eRAvvPA8enouxqpVa/Dyy9/C/PwOAKsBDGJ+\nHhgdHUdXVxdef/117Nt3EC+88DxWrcrjAx/4Lezdu9dXfjdu3IjR0VFUKhVs2XI7FhaMounuPopN\nm/ZgZmYm0P1bf1v/9/t7Hc4/fPgwzp49Wy3/JImqO8AKV3s7FhbOAwC6u1/EX/3VMXR1dWFmZgav\nv/46JiePYG7uNWzf/h7s3bsXlUoFo6PjAIDx8VEUi0WMjo47tHYe7373dvT1XYU/+IOH0dXV5Ss/\n9ufqdT+bNr0Ds7O7sbBgPY27ALwb3d0jKJWOKdWF/W9nmlHSC/LeO+/n6NGjAJBZ3VmfzczMmNqZ\ngBHm9BAWFnrxwgvfxMGDD1R1aD//9OnTePLJL+HChVfR0fH7eOONCoBrYGigpuPHHjvWcH69bp6G\nEJ/Ciy9eg0qlgmKxiMOHD2NhYQFPPvml6u/Xr1+fSn1k/R1WJ7K/m71Xieku7JDe74HGafZDAEbM\n/9+HhBfAuS3wGRsbo6VLr6COjjdTX9862rlzp2OluhVtSD4db+xZX2H7/WWUy70p1FSSlUcVASjS\ntnm3o82ciGjnzp2mk443U1/fDXXT67IwuDJbomwhpaG/S6mjY4lSe6N1vb6+G6in52rq6bma+vpu\nUL4Azn4tK20d7ZdJaY4U6s7uZKq2ir1+itzyo2HHqcmOjiXU3f3m6gI457S4bKrcy2nRxMSEsh02\nbDP30FHYE30lDnwWwLcBLAD4JoD3w9ia9iRS2JrmtsDH6YHLcMEq21ZhLSCpd7VprWhv/P2yRP0V\ntypZa8zHxsZMd6wbyAo3aXna8uuvwKoo6/W6iqwFTMBVyuzY9Q6RaotA49jWmJWtk1lrzBs9Py4j\nISwfB/KtaZYea/vIiWqe+eTl41V+btsOs7IdUQe0bcxDZyqGxrxcloccdYu/K//tABnb15aQ3T2n\ne9pXZTZutE5kqTF32wpk+cD215jX/GCPjY2ZHrjqHX+obMy9XICqrnTDVuxJL5zKWmMu309+Gcm2\nRxpeKQckA5tyUw14lV+ajXmrLKyLortctEl6fbHbVqzFIHNzK32fn88vRWfnvTAWcxyDYTs6A+Bn\nAP4AwP8E8D+xsHAYk5NHsHr15QDusf3+HgACq1dfpeQe0koj7fOzxuTkESwufhyGnXsYwASAf6x+\nXyrtQnf3CCyddHePYM+e99s+uwfApzA39wCmprZgfPyT+MVfLAB4FsArsLTY0fF9lEq7AuVNZVk0\nS6tSqWDz5m3YvHkbKpWK52/n5l5rmtbWrcOYmtqCqakt2Lp12DXNdtMbYDyfc+e+Bsu2bbFkyWUA\nfgPAbtTqpREAG821P5ZNfRjG4rX9MCZSAWAmcD5k2i6VdmHTpndIPw+Ds17fvHkbbr55EFu2vNeX\nPmTpREUbzYXtBcR5QLHNvNYzlG0Bk0+zW446LIcw3d3W1JO8l2nYzK1pLcOumcsFt2u63UNaaaR9\nPlFyoyQVupPbufOece/tn8lmePr7bzUdeTTa4IPgtY4k6DR7M/u721RsmH3mQUZ27Wgzr20Nq80I\n5XLLaWxszKzXtlHNrWvJ1YxoOb4yzpH7yGhmJpFpe3p6WtnI2c03e9DtvGwzz4i4ndRPI1r7IQcI\nuK4aerSZ2Nz2DjsrKqc3OCY6WWrMGyuZywI1vmnZF2sLQzdWQ1mG1W+zewhasafxTLLXmFt1U31s\ncbdybdbh8iofHaa03X0iJPfOxAE35k2oCde/Yxj3NI6SV2QqRj1ZasyJolV2WVkg5oXqxjeNZ5Kl\nxjzs89GhUQ6LvDHfkNl3xoIbcwmyEKa1rRO1MKbObRpu51tpWOL3E+lJhynqtPPQbtPsKrYVxrV1\nK86pRee7kYa7Wr9p+SFLjTlRvdvUsKNp+290dLFrT6tx9f5Kz1gCSeRJBdyYS7A/YLtIh4eHpXsh\nndj3bDq/d9srLLMVqbqHtNJI+3yi7DTm5XLZtvZixDP0bZBRke6Nud/3IUi+wo4a27UxJ/KzJsK9\nzmtsHJe7lmHQsolTd2FnFvz4XEjj/rgx98ApUj+Re5pVTm5hBbM8NaozWWnMZbqQxY6WadJthigL\ntMK0upMsNuZu+Cmfxt+U6gY9lkZ1KJs4Sfv+ouguVd/sSTA5ecS2BQNYXASARwKdMz8PPPhgCYuL\nkwCAXK7UcM7i4lvrfj85eYR9ELcZL7/8LV+fyTT54IMl3HLLLalrplKpYHLyCABjq1Ea+ZG9f/w+\nJc0/mlss6zW6du2NocpGB135Icvaa4t95o18Hdaex87Oe+v2PFYqFTz11Dk492zWGuthLC7uQC5X\n27tp/H9jwDxEvYdk0kj7/CxR72vgPgD3mJ81Z3HxrdXKzsnMzEygfdtu+NkbHmY/t9v+YlX5CkI7\n6c2J2737KZ9SaRc6O+8G8E7zONeQzuLiW6Wd02YcOnTIt66aoap8W1FzLT8yL5V24dSpYczPW5/c\nA2APgBMAvo23v/1t1V5XLdLQHQD+CMD1AIBcbjcWF3faUr0Ja9fegN7eEwCAgYESxsc/ifn5mwCg\nGpiCaS8OHnwAW7a8FwsLjwD4ITo738DBgw80/K5U2oUvfvF2c5YIMBx53AG3SFOnT5/G/v0Pxx49\nLeyopFgs4vjxY7aRV7S8Od9Zfp+i4b98LgbwQQBALncXiH4PhhUAsDS6evU/Y35+JFDZPP7432Zm\ntJtp7YWdn4/zQAxb0xp9EDfajurtRmUCNlBPz9W+VudmeZuH7iAjNnOiRh246WJsbMxcZ1Fz5OGm\nm6T2WevkQzvt9ykpzVECNnM/yMq+u/sKEmJ5g0az4CcgCmlqL4ruUm+4pZmKQdyWQxejAi1JG2WZ\n6JxBCbixTp4sNeZ2wnjLkpFUZRjn4p+svT/cmBv+/4XIe0bP80Pai8qyBDfmEry2zuRyy6Ue2spl\ne5CMEbKCZIStOHXY1pV2Htppa5qd6elpZY2wqhCSae3nblaZ67jnN4uNedSwtPWeCy+pDnpyueWR\nXbCq6szpvk0zKlF01/I2c0C+eri390SDzaZYLGLt2htw5swjMOxHx2AEt5DbMhkmCdavX6/UJu1F\nsVjUxhbPJIdlV3/f+34Xc3OvA3gPjOArRn0Ztbzi0BXjIGwvIM4Dsflmbz5C4ikh/UAGR+ZErCWL\nrNlMibI5MleBMTu5InPl1SpE0Z0wztcLIQSpzFelUjFXGV8HAOjsfA4nTnwOxWJRuv8xK3si2wUh\nBIhIJHAdpboD6vfXDgzcjNnZrwBIVldp67m2S8RYjd/dPRLLanyVJKU581rKddcML11++ctfxoMP\nTpr7zLNRXq1CJN2F7QXEeSAGm7k9xKnlZtNr5NQK9mod8tDONnM7UUbp6myhI8pmB4LmyctmqqP9\nMinNkULd+b33Rvv4pQ2Lgr3cWceVryTT0jFPRNF01zY284WFj8Ky2S0soNorZVsekwRp2Y3rrzuD\n+fnrU9E420z1walFgxMAPlbV5f3334m9e/emlEMmDC3rAW5wcBCA3aNbPU89dU76ufN8FXlI63wd\n8qDiHrLI4OBgnde2Cxe+FyktRblSlE6wPDXzXqdSI+2qNyCd5+jHM6E9raieDFXdY0tqLuyQPs4D\niqadDMccVmCVXse00nUEbDP/394LlHQHGZxmb4xCtazO1JNUzOm0F+Flde96Upojxbrzg59p9mbn\nWwGnOjuX+S7bNLTYTv4NUm+4pZlSIO6JiQnHqkzDoxuwvCpcYx/5Nsrn+xoKuxXs1TrkoV1t5v39\ntzasCLbiLfutWGqV30ikys+q0KLEWHfit1z9rGQPoxG3hoFt5v6wN3JBwjc3dgR6zbpVXrb2tFTs\nalDl34Bt5hni8cf/1gyOYlGEsWf8QwCGzL8B4BH096/FyZNPJJ1Fps3o7V0VSGc12+ZqAIOh7eyW\nvXpmZgavv/46Nm/eBiDbOzXc1iDcf/+d6WYsIzjXMIyO+jtPbm8/glp9Gpy4dlu0nX+DsL2AOA8o\n6KkavcCSOfq29yKtz8oEHCUh8pmYfmlnkMGRudsOiiC0QpzwuK4Z9971pDRHinWnAq+pabnb1w2h\np9n9xL0IS7v5N0i94ZZmSlHAC0MkJQKuMm3kAwRYjfwGApZRR8cSbsw1J7uN+TJTZxuos3NZYJ3V\nV34lyuVWUKEwEFqvaVVucdgt4+6YtGtj3kxzjWtBVlbNR2HWgMSpybTXi4SBG3MJExMTVCgMUD7f\nR296U57qF8D1EnBFdXTuZeeJgg725rTzwDbzaJVUuVyma6+90RYzIHzUPlV5Ioq+992e17Bpye6Z\nbebRqDWuO8lYX9SouaCdM6981WZQb6sOtLw0qcq/AdvMM0KlUsEDDzyEhYWPm5/sgeFn2G7nOQzD\nzpORWLVMW1IsFnHZZeNYXLwHbrY/p4c1t3jn27e/B88+GywWtWpked2/f0+o7T28dz0ungbweQCf\ngExzKp/7wMDNmJo6ZF4LAO7CwIC6/e1tpZGwvYA4D0TsqTZO3WyQ2HmuyszUS7uDDI7MVU7xNZuK\nDBp7IM2tOlmxYyalOVKsu6jUfLM31plxlFNW9JAUUXTXkiPzRjYCuNv29z1YuvRivPOdJ2KNQMW0\nL1YUKhWRzkqlXTh1aljJiLqtRipMYIzIkTfizJlbAIxUP8/ldqNU+mx6GWOaE7YXEOeBiD1VY/FR\nzd7T2bmSOjqWBFqM1Ar2ah3y0K42c9U2Oa8Rtd9ZAB3shLK8TkxMpJ4vJ0lpjhTqTtW91+pPY6Fw\nLreCxsbGYslX0Bks9s3eZiPzYrGIAwfuw5NPngAAlEqPAoBtlLSfRydMpvAaUaucBYgbWV67urpS\nzhVjp1Z/fgnAlbHWl1nSru60RQhUJttkOQQqk01aPQQqoydRdNeygVacRHXwzzBRYQ0yrQTrWTPC\nzs/HeSCGeOZBVxa3gr1ahzywzdwgrXjmcaTTDmklpTnS0GbuJ60getaxfHXME1E03bX0yNzqOb7v\nfb9r89Fr7HG1bDQMkwT1fqKzoUEeeWWL06dPJ1ZeWdRzyxO2FxDnAaXuXGu+g3kvYzZBBkfmTuLe\nT6t6/3gWXWGqJCnNUSz1XfzlpUrPafs90I0ouku94ZZmSoG4C4UBm9jKZHfn2m4VU9ZphcY87rje\nqtOWVdaFwoCS/GaBrDXmSTtf8as5FVsq2wluzB2Uy2US4lKHuI2AK/l8X2R/wn7Rwd6cdh7YZl4j\n7Cik2TP0W5EHKQtZmrncCu3jQberzdworxFljbmf59hMz7XGekTaWIfpgLDN3P1oyX3mk5NHQLQO\nwG7bp38CYBP6+4n3MTKxUqlUMDo6jnx+ZV185ix5XyuVduGLX7wdi4vWJyNYXNyBffsOYHLyCObm\nXsP4+Ghm7qfVKZV2YXb2diwsXA+g0UtglJjhbuc203PNrr4awGDrxxNPm7C9gKgHgJcAfBXAGQCn\nHd9F6t0UChvNafVtBORNm3mJgEsjeTJi0gEZGpm3Usxw4z3aQEY0qzIZITG9o7e1CklpjhTpjsh9\npBxFH1HObTby5mn2RqLoLs3G/EUAeZfvIj2Qmr2cnfi3AllqzFs5ZrgRgKM93qcsNuZuRNFklHP9\nNNa8AK6eKLpLe2taLB6WentXADgfKY2ZmZnI+YiaRivkQcU9ZJMZdSn5eIbFYhEnTz6BkyefcJ3G\nDFoWlqvNoaETGBo6gbVrb7SnFigtL1RqpH31pt9ztPTT3/9pDA2dkIbl9aNb1flSmY7qtKKQps2c\nADwphPgZgD8mok+pSrhmP9oJ4J7q52nEb2bai1qEsx0AXkZn5724cOFt2Lx5W2BbpQ7Y7aK1WOQA\ncB7d3Uf5fdIEt3UaQLSoe1Ej9hWLRXR1dYWKV88EI83GfCMRfUcIsRLAlBDiOSL6B+vLHTt2YM2a\nNQCAZcuWYd26dVVBWD0ht7+7urpwxx1bcfToUSwurgSwH0LMYXR0L4rFYtPznT0tv7/X8e/BwcHM\nnX/48GGcPXu2Wv5JEkV3MzMz6OrqqgaOePHFMl566Sc4c+YDAIDZ2dtx4MB92Lt3r+/0VP5tfRb2\n/K6uLuzfv6cawGjTpj11QVJ00X2Y/MzMzODo0aMAkDndHTp0CA888BAWFj4OoFFnYcsNMBrj/fv3\n4PHHP212FIzAOEF0ZH2mQ33mvL+o+Yl6fyp1p0WgFSHEhwH8mIgmzb8par42b96GqaktMDwUAYAx\nXXjy5BPRMsskTlYDrbAGs0uWAq2wzlqHzAVaEUJcIoToMf+/BMBmAE+rvMbc3GuRznf23tJIoxXy\noOIesojK+1aVlo550jmt7DGjLiVNy6TV34UopDXNvgrAcSGElYe/IKKTKi+wfft78OyzI6FtPQwT\nlaj2Robxg3OdBuusPdFimt2JqunOKI4SGH3I6jQ7wBrMKlmaZgdYZ61CFN21dGPOtAZZbsyZbJK1\nxpxpDTJnM08CHWy9nAd97ElJo6NNTsc86ZxW1tD1OeqYlo55ikrLNuYMwzAM0y7wNDujPTzNziQN\nT7MzacDT7AzDMAzTxrRsY66DrZfzoI89KWl0tMnpmCed08oauj5HHdPSMU9RadnGnGEYhmHaBbaZ\nM9rDNnMmadhmzqQB28wZhmEYpo1pyca8Uqngllvehc2bt6FSqYRKoxXs1TrkQRd7UtKotslVKhVs\n3rwtdU23U1pZw3nvUTSja5mwzdydNEOgxkIt5vIOANfj1KlhHD9+jN0bMpnl9OnT2L//YczPTwAA\na5ppSq0eZM20Cy1nM+dwgK1Hu9vMWdPJk3WbOWsmm7DNnGEYhmHamJZrzEulXejuHgFwH4BjZjjA\nXYHTaQV7tQ550MWelDQq73vTpneYmj6GtDXdTmllDfu91+rBcJrRtUzYZu5OyzXmxWIRx48fQ3//\nP2Fo6ATbiZjMs379ehw/bkyTsqYZP1j1IGumfWg5mznTerS7zZxJnqzbzJlswjZzhmEYhmljWrYx\n18HWy3nQx56UNDra5HTMk85pZQ1dn6OOaemYp6i0bGPOMAzDMO0C28wZ7WGbOZM0bDNn0oBt5gzD\nMAzTxrRsY66DrZfzoI89KWl0tMnpmCed08oauj5HHdPSMU9RadnGnGEYhmHaBbaZM9rDNnMmadhm\nzqQB28xtqAgVyTBMvPB7mhz8rNsEItLuMLIVnHK5TN3dqwg4SsAIdXevonK5HCqt6enpUOepTKMV\n8qDiHkw9aKs7GSruW3VauuSp/j09Sp2dy0O/pyrzZScpzZFC3cnu3fms/daJumglrrR0zBNRNN21\n1Mh8cvKIGb93GMC7MT8/gcnJI2lni2EYG/Xv6TAWFnbxexoTzmfNdWLr0lI2c47h25qwzby1yMJ7\n2io28yw8a6ZGFN11qM5MmpRKu3Dq1DDm542/jbB/x9LNFMMwdfB7mhz8rNuHlppmt4f96+//dKSw\nf62wx1uHPOiyBzNpdNzHqkuenOE59+/foyw8Z7vqDZDfe9hQqLpoJa60dMxTVFpqZA4Y4i0Wi5iZ\nmcHg4GDa2WEYRoL1ngL6VIativ1ZM61LS9nMmdaEbeZM0rSKzZzJFrzPnGEYhmHamJZtzHWw9XIe\n2ncKVUebnI550jmtrKHrc9QxLR3zFJWWbcwZhmEYpl1gmzmjPWwzZ5KGbeZMGrDNnGEYhmHamFQa\ncyHEu4UQzwkhnhdCjMRxDR1svZwHfexJSaOjTU7HPOmcVtbQ9TnqmJaOeYpK4o25EOIiAP8NwLsB\n3ADgdiHE9aqvc/bs2VTP5zyoOT+rqLxvVWnpmCed08oauj5HHdPSMU9RSWNkvh7AC0T0EhH9FMDn\nAPya6ov84Ac/SPV8zoOa87OKyvtWlZaOedI5rayh63PUMS0d8xSVNBrz/wDgm7a/v2V+xjAMwzBM\nCNJozBNZtvnSSy+lej7nQc35WUXlfatKS8c86ZxW1tD1OeqYlo55ikriW9OEEBsA7Ceid5t/7wOw\nSEQTtt/wPg2mjqS2psV9DSY7JLk1LYnrMNkgrO7SaMw7APwLgF8E8G0ApwHcTkTnE80IwzAMw7QI\niUdNI6I3hBC/B6AC4CIAn+aGnGEYhmHCo6UHOIZhGIZh/KOVBzghxH4hxLeEEGfM45ds3+0zncw8\nJ4TY7JFGKIc0QoiXhBBfNa972vwsL4SYEkJ8XQhxUgixzPb7zwghXhVCPG37zOv3dfl3OT/Q/Qsh\n3iKEmBZCPCOE+JoQ4q6A+fhNl/N950MI8SYhxJeEEGeFEM8KIQ4GzMOvuJwfWQs+y/0/m/f/MyHE\nzY7vAl8nikOkoJpqklZgbbikE7h8feTtIrNM/ybiPQZ6Z5uktUwI8XkhxHnzPt8R5R59XO+j5rXO\nCSH+Wghxme27xHSnSieONFWVr7IyMZ/pM0KIp4UQjwkhuvykFfSd9Co7l7TU6YCItDkAfBjAHsnn\nNwA4C+BiAGsAvAAgJ/ndReZ3a8zfngVwvc9rvwgg7/jsEIC95v9HADxk++5WAAUATzf7vUv+3yU5\nP9D9A7gcwDrzN0thrEW4PkA+XgRQkJwfNB+XmN93APhnAP8p4LNYIjk/khYCaO46AG8DMA3g5ijX\niaK/oJrykVYgbTRJy3f5+szbHgB/AeBExHsM9M42SesYgJ22+7wsyj36uN6QpScADzV5P2LTnUqd\nxFC+SsrEfC7/CqDL/PsvAQz7SQvR6/lck7SU6UCrkbmJbCXfrwH4LBH9lIhegnFj6yW/i+qQxnnt\nLTAEBfPfX7e+IKJ/APB9n7+X5X9Bcr4sD27nryeiV4jorJmfHwM4D2PPvt98/AuALsn5QfPxE/P7\nThgVy/cDPoubJOcHyoPkd74goueI6OuKrhNJfwE11SytoNrwSitI+XoihLgKwC8D+BPUyjdUWlaS\njr8Dp2WOhm4los8AxroeIvphxHx5QkRTRLRo/vklAFeZ/09Udyp1AqgrX8Vl8u8AfgrgEmEswL4E\nxuLrpmkpqOerZSdLS6UOdGzM7zSnHD5tm764EoZzGQs3RzNRHNIQgCeFEF8WQnzA/GwVEb1q/v9V\nAKuapOH2e7/5B0LevxBiDYxe35fC5MN2/j8HzYcQIieEOGtea5qIngmYh6sk54d+FooIc504HCIF\n1WADPrXhdX6Q8m3GxwHcC2DR9lnYtFS8swBwDYDXhBB/KoT4ihDiU0KIJRHyFZSdAP7e/H9quouq\nExNV5ausTIhoDsAkgP8HRiP+AyKaCpkvr/uJWjdF0kEavtmnTLuF89gC4I9gFOI6AN+BUQBuyFbu\nRVnNt5GICgB+CcDvCiFurUvYmPvwnb6P38u+C3X/QoilAJ4A8PtE9KOg+TDP/7x5/o+D5oOIFolo\nHYxe5buEEL8QMA/O8weD5sHjOy/N/arXeUGv4+P7SATVIBBZG9bvopavlZf3APguEZ2BfNYl6D2q\nemc7ANwM4A+J6GYA/y+A+yLkC4A/3QkhRgEsENFjHknFrjsVOlFcvsrKRAjRB+BuGNPVVwJYKoS4\nI2S+6ghZz8vyGFkHaWxNG/LzOyHEnwD4G/PPfwPwFtvXV5mfOXH+7i2o79145es75r+vCSGOw5jS\neFUIcTkRvSKEuALAd5sk4/Z7X/knomr6fu9fCHExjJfwUSL6Qoh8vGKe/+fW+WHyYZ73QyHE3wHo\nD/MsbOffQkQzYfIgw6/mHAS+juQc3/rzIKgGqwTURlN8lq8X/xHAFiHELwN4E4BLhRCPhs2ToncW\nMMroW0T0v8y/Pw9gH4BXwj4rM1+euhNC7IAxJf2Lto8T151CnagsX5VlcguAfyKi75n3+9cA3hky\nLXjcT5iyU6YDrabZzQdjsRWAtervBID3CiE6hRDXAHgrDGczTr4M4K1CiDVCiE4Av2Ge2+y6lwgh\nesz/LwGw2bz2CRgLJWD++wV5ClXcfu8r/0HvXwghAHwawLNEdDhkPj7oPD9gPr5hTYELIbphLOg4\nEyAPPw/Ddl93vhDi8iDPwvksQ2IfTYS5Tij9NSGoBgEAIbThlk5vwPJ1hYjuJ6K3ENE1AN4L4P8m\not8Kk5bCdxZE9AqAbwoh3mZ+tAnAMzA6kIGfvR+EEO+GMR39a0T0/9m+SlR3qnQCqC1fxWXyHIAN\nQohu8343AXg2ZFpAxHrejlIdkKLVmSoOAH8G4KsAzpkPaJXtu/thLAJ4DkDRI41fgtE4vABgn8/r\nXgNj5eBZAF+zzgOQB/AkgK8DOAlgme2cz8KwvyzAsFe9v8nv6/IvOX9n0PuHsap40cz3GfN4d4B8\n7JGc/0tB8gFj8dpXzDS+CuBeH8/OnsbvuJwfWQs+y36r+fznYcxS/Pco1wmjv7CaapJWYG24pBO4\nfH3mbwC11c6B00KId7ZJemsB/C9Tb38NY+V0pHtscr3nAbxsK5s/TEN3qnSiunxVlwmAvTA6VWPI\nkgAABTVJREFUA0/DWLR2sZ+0gr6TXmUnSWunSh2w0xiGYRiGyThaTbMzDMMwDBMcbswZhmEYJuNw\nY84wDMMwGYcbc4ZhGIbJONyYMwzDMEzG4cacYRiGYTION+YaIIRYIWqhPr8jaqE/F0VjGL27hRB/\nmFZemewgjLCuZ4QRJvSvTZedQc6fEWZYWCHE3wkhLo0np4yOmPXPx2x/3yOE+HCaeWLc4cZcA4jo\ne0RUIMPP9CMAHjb//zswPCnZ+Q0AXv57GcbiJ6au/jcYkaN+J+D5VScURPQrRPTvSnPH6M7C/9/e\n2YVYVUZh+HmRsHFKo5+bgcwwunAaaZoaJbIL8y4qy1IwAiEKMjCyiCgiIiooKiVzKMFCBqYfCzKI\nTERz/Bm9GH9mMm9Coe5CtHAaoYvVxVqbsz00MCnTOfuwHhhmn7X3+c63D4u99lr7O+sFHpR0Tbxu\niaYkcuW0liODeXNStBX9Cri3cD65qlGHme1t0LyS6nIAmAsgqVfSfrka1b6iZWa0u/xM0vHoX91W\nvFnSKUlXx/Za1URDnmnEyST/C38DHwPP1u+QdJ2krZIOxd+dYT8maaac05IeC/sWSUskdUo6GBWj\no5LmRhvaE5L6w/e+jNbBSHolxh+R9FHp83dLWhfjjEi6I+ztkjbHZwzLBbyQtErSNkk7gR1T/s01\ngAzmTYy5dN8hvAk/eJb+eeNmlFQRSdPw3uWjYfoZ14q+DXgVeDPsTwHnzGxe2HtKw1iM1QOswkVN\nFgJPSLp1qs8haRgbgUf/5RHLeuB9M+sFHsb1ywH24S1iO4FfYhvcV/bj1aH1UXnsoSYecjPwYfje\nn8DqsG8ws14z6wLa5Mps4P7YFuOsBjaH/WVgp5ktABYD70iaEfu6gWVmdoHqX6uQwbz5GaBWal8R\nr5NkMrRJOoxLyF6PP8IBuArYKmkEeA+YF/ZFQD+AmY3gvdjLCL84f21m42Y2hvfMXkTSkpjLoW4B\n1tTtWgJsCP/6BrgyBG8Ggbtxn+gD5kvqAM6Y2V94heglSS8Ac6wmLvKrmR2I7X5qNwGLJQ1JOoYH\n58JXIa6FZjaIK7TNwm9aX4x57QKmA7Px4L/DzM5e+rfSnGQwb362AfdI6gZmmGsFJ8lkGI/M5Qbg\nPPBA2F/Hs5cu4H5K5XQm0KEuYXXHiBZ5lppMyDrgcaC9ZBOwoFjrY66WNgbsoRbMdwO/45n7IICZ\nDQD34cJG30kqsmSrG9skTccrA8ti3ccmXFp1IooxHirNa46ZnQj72EWce2XIYN7kmNk5/A7zE3Lh\nW3IRmNk4nlm9ERKQM3H1JvCSecEeYCWApFuA+fVD4RflpfF8vR1YGrakRTGzM8AXeEAvAuYPlLL1\n4lGLmf0GXAvcZGYngb3A88CPcdyNZnbSzD7AM/quGGK2pIWxvRL3qcvj807HLzEeKU1LeKUSSXcB\nZ2OB5va6eXWXjm9pMpg3J/WZzgDu9FliT/4L5dXoR3A5xeXA28BbkoaBaaXj+oArJB0HXsN1si8c\n0CtDn+JrOYaATWZ2dArPIWkc5evQu3iQLlgD3B6L2H4CniztG8LlQcGDeUf8B1guaTTK4J14CR9c\nvvXp8L1ZQJ+Z/YFn46PA98DBurmdDx/eiN9ogFedLouFeKO4HxfHt3QFKSVQkyRJkoYRv9L5Nh77\nTPY9u4DnzGx4quZVNTIzT5IkSRpNZpWXSGbmSZIkSVJxMjNPkiRJkoqTwTxJkiRJKk4G8yRJkiSp\nOBnMkyRJkqTiZDBPkiRJkoqTwTxJkiRJKs4/TDISZIdFdRgAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# visualize the relationship between the features and the response using scatterplots\n",
"fig, axs = plt.subplots(1, 3, sharey=True)\n",
"data.plot(kind='scatter', x='TV', y='Sales', ax=axs[0], figsize= figsize)\n",
"data.plot(kind='scatter', x='Radio', y='Sales', ax=axs[1])\n",
"data.plot(kind='scatter', x='Newspaper', y='Sales', ax=axs[2])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Questions About the Advertising Data\n",
"\n",
"Let's pretend you work for the company that manufactures and markets this widget. The company might ask you the following: On the basis of this data, how should we spend our advertising money in the future?\n",
"\n",
"This general question might lead you to more specific **questions**:\n",
"1. Is there a relationship between ads and sales?\n",
"2. How strong is that relationship?\n",
"3. Which ad types contribute to sales?\n",
"4. What is the effect of each ad type of sales?\n",
"5. Given ad spending in a particular market, can sales be predicted?\n",
"\n",
"We will explore these questions below!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simple Linear Regression \"Sales\" régréssé contre \"TV\"\n",
"\n",
"Simple linear regression is an approach for predicting a **quantitative response** using a **single feature** (or \"predictor\" or \"input variable\"). The model takes the following form:\n",
"\n",
"$$\\boxed{Y = \\beta_0 + \\beta_1X + \\xi}$$\n",
"\n",
"What does each term represent?\n",
"- $Y$ is the response\n",
"- $X$ is the feature\n",
"- $\\beta_0$ is the intercept\n",
"- $\\beta_1$ is the coefficient for $X$\n",
"\n",
"Together, $\\beta_0$ and $\\beta_1$ are called the **model coefficients**. To answer the questions, we must \"estimate\" (or \"learn\") the values of these coefficients. And once we've learned these coefficients, we can use them to predict Sales!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Estimating (\"Learning\") Model Coefficients\n",
"\n",
"Generally speaking, coefficients are estimated using the **least squares criterion**, which means we are to find the line (mathematically) which minimizes the **sum of squared residuals** (or \"sum of squared errors\"):"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What elements are present in the diagram?\n",
"- The black dots are the **observed values** of x and y.\n",
"- The blue line is our **least squares line**.\n",
"- The red lines are the **residuals**, which are the distances between the observed values and the least squares line.\n",
"\n",
"How do the model coefficients relate to the least squares line?\n",
"- $\\hat\\beta_0$ is an estimator of the **intercept** $\\beta_0$\n",
"- $\\hat\\beta_1$ is an estimator of the **slope** (coefficient of $x$) $\\beta_1$\n",
"\n",
"Predictions at the points of the design are\n",
"$$\\hat y_i = \\hat\\beta_0 + \\hat\\beta_1 X_i$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's use **Statsmodels** to estimate the model coefficients for the advertising data:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Intercept 7.032594\n",
"TV 0.047537\n",
"dtype: float64"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# this is the standard import if you're using \"formula notation\" (similar to R)\n",
"import statsmodels.formula.api as smf\n",
"\n",
"# create a fitted model in one line\n",
"lm = smf.ols(formula='Sales ~ TV', data=data).fit()\n",
"\n",
"# print the coefficients\n",
"lm.params"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Interpreting Model Coefficients\n",
"\n",
"How do we interpret the TV coefficient $\\beta_1$?\n",
"- A \"unit\" increase in TV ad spending is **associated with** a 0.047537 \"unit\" increase in Sales.\n",
"- Or more clearly: An additional $1,000 spent on TV ads is **associated with** an increase in sales of 47 widgets.\n",
"\n",
"(Note that if an increase in TV ad spending was associated with a **decrease** in sales, $\\beta_1$ would be **negative**)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Using the Model for Prediction\n",
"\n",
"Let's say that there was a new market where the TV advertising spend was **$50,000**. What would we predict for the Sales in that market?\n",
"\n",
"$$\\hat y = \\hat\\beta_0 + \\hat\\beta_1x$$\n",
"$$\\hat y = 7.032594 + 0.047537 \\times 50$$"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"expected Sales = 9409.444 widgets\n"
]
}
],
"source": [
"# manually calculate the prediction\n",
"sales = 1000*(7.032594 + 0.047537*50)\n",
"print(\"expected Sales = {} widgets\".format(sales))"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([ 9.40942557])"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# use the model to make predictions on a new value\n",
"X_new = pd.DataFrame({'TV' : [50]})\n",
"lm.predict(X_new)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plotting the Least Squares Line\n",
"\n",
"Let's make predictions for the **smallest and largest observed values of x**, and then use the predicted values to plot the least squares line:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"[]"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEPCAYAAABCyrPIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXt4VOW1/z9vwMTIxRiigDfAW1VEGWkVy7HEKkRrDwo8\n2urRH9gqeurdIFGsFmtQQaOeaqsHDwq2Xh8pHmxrLh6JLXrUygFFxFsL1EulXJSWGg2Q9ftj78lM\nkpnJntnvntk7rM/z7Cez98zs/Z03M+/a71rrXa8RERRFURSlqNACFEVRlHCgBkFRFEUB1CAoiqIo\nLmoQFEVRFEANgqIoiuKiBkFRFEUBAjQIxpjdjTGvGmNWGmPeNsbc5h4vN8Y0GWPeM8Y0GmPKgtKg\nKIqieMcEOQ/BGLOHiHxhjOkNLAOmAxOATSIy1xhTA+wlItcFJkJRFEXxRKAuIxH5wn1YDPQCPsMx\nCAvd4wuBM4PUoCiKongjUINgjCkyxqwENgBLRWQ1MFBENrgv2QAMDFKDoiiK4o3eQZ5cRNqAkcaY\nPYEGY8xJnZ4XY4zWzlAURQkBgRqEOCKy1RjzW2AUsMEYM0hEPjXGDAb+1vn1aiQURVFyQ0RMru8N\nMsuoIp5BZIwpBcYBK4AlwBT3ZVOAZ1K9X0RCv/3kJz8puAbVqTpVp2qMb34JcoQwGFhojCnCMTy/\nFJH/McasAJ4yxvwQWAecHaCGQFm3bl2hJXhCddpFddolCjqjoNEGgRkEEVkFHJvi+BbglKCuqyiK\nouSGzlT2wdSpUwstwROq0y6q0y5R0BkFjTYIdGJarhhjJIy6FEVRwowxBgljUHlXoLm5udASPKE6\n7aI67RIFnVHQaAM1CIqiKAqgLiNFUZQeg7qMFEVRFCuoQfBBVPyKqtMuqtMuUdAZBY02UIOgKIqi\nABpDUBRF6TFoDEFRFEWxghoEH0TFr6g67aI67RIFnVHQaAM1CIqiKAqgMQRFUZQeg8YQFEVRFCuo\nQfBBVPyKqtMuqtMuUdAZBY02UIOgKEqPoqGhgfHjJzN+/GQaGhoKLSdSaAxBUZQeQ0NDAxMnTqGl\nZQ4ApaU1LF68kKqqqgIryw9+YwhqEBRF6TGMHz+ZpqYJJJZtX8i4cUtobFxUSFl5Q4PKBSQqfkXV\naRfVaZco6IyCRhsEtqayoihKvqmunsayZVNoaXH2S0trqK5eWFhREUJdRoqi9CgaGhqoq5sHOAZi\nV4kfgMYQFEVRsqanGg2NIRSQqPgVVaddVKdd8q0znonU1DSBpqYJTJw4pdv01Ki0pV80hqAoyi5F\nXd08Ny3VyURqaXGO9ZRRgh/UZaQoyi5FT05N9esy0hGCoii7FJqJlB6NIfggKn5F1WkX1WmXfOmM\nl7Soq5vHDTdczrhxSxg3bomnmcxRaUu/qEFQFCVQbNcWyuV8nQPJs2ffS3X1NKqrp1FXN0/rHsUR\nkUA24ABgKbAaeAu4wj0+C/gIWOFup6Z4ryiKEn3q6+ultHSgwAKBBVJaOlDq6+vzfr5x4ya57xF3\nWyCx2Fir2sKA23fm3G8HGUPYDlwtIiuNMX2B5caYJkCAu0TkrgCvrShKCLCd0WPzfOvXf6TZRp0I\nzGUkIp+KyEr38TZgDbCf+3TOUfAwERW/ouq0i+q0Sz50VldPo7S0BlgILKS0tIYhQwZ5fn9U2tIv\neckyMsYMBWLAK8AY4HJjzP8DXgeqReTzfOhQFCW/2M7oyfV8VVVVLF68MGl2svMep1S2HW09gcDn\nIbjuomagVkSeMcbsA2x0n74FGCwiP+z0Hglal6Io6bFZ2sF2mYgways0oa5lZIzZDfgN8JyI3JPi\n+aHAsyIyotNxmTJlCkOHDgWgrKyMkSNHUllZCSSGb7qv+7pvf3/u3LnceOPttLbeDUBx8dXccst1\nzJgxIxT6dD+x39zczIIFCwAYOnQoN998sy+DEGSWkQEeAe7udHxw0uOrgcdSvNdCvD14li5dWmgJ\nnlCddunpOlNl5IwbN8muuCSCbs/6+noZN26SjBs3Kecsoqj8zwlxltEY4DzgTWPMCvfYTOAcY8xI\nnGyjtcDFAWpQFGUXpvOSmsuWTdmlltTMFq1lpChKB8K0LrFfH39PrluUCq1lpCiKVVJl5BTKGOjd\nfZ7x428KakNjCFZRnXZRnXZJp9NGLMPWTOmotCU+Ywhay0hRehDJdX5ee+01T6/ryTV84qOdbArZ\n7SptkxI/1iSojYiMEBQlTHi9G7ZdXygoCqEzKm2TDnyOEAre+acUpQZBUbLGq4sl32mlfrCRMpoN\nUWqbVPg1COoy8kF8gkjYUZ12iYpOp3xY+MnUnlVVVTQ2LqKxcVFBg8nR+Z/7Q7OMFKWH0LnOT3Hx\nPKqrH+/2dUHU8IlqSYhdfTU1nYegKD0Irx1xkB12mOYx5EJUjRmEvJZRrqhBUJTosqtNBgsTfg2C\nxhB8EBW/ouq0i+q0SxR0RkGjDTSGoCiKVXZ1P3yUUZeRovRwvPjEw7xmgeIdjSEoyi5Ops7XS4A3\n6kFgJYFfg1DwSWipNiIyMS0q9U1Up13CpDPTzNqlS5d6mmiVeE29wCSB0RKLjcnbZwhTe6YjChpF\ndGKaouzS1NXNc+/spwDOXX58tJAdq9xzTAAu4Y033s5Yx8dGvZ/4OaZPv2nXqxkUVvxYk6A2IjJC\nUJRC090IwEttnvr6eikqGuC5ZIONej9RrxkUVtBaRoqy6+K1w++uHlAsNtazQbBR7yfqNYPCil+D\noC4jH0QlN1l12iVMOjOVd47rzFQPKO62ASguvgpYCCx0U0Wn+dbnzbXU7Ps6QROm/3mQ6DwERYk4\nVVVVOWUEdc4uKi6+lljsQSoqBmZcJc3rPINMK551PMcaSksX6FyFEKBpp4qyi5JLiYl4iuumTZuB\nHa7xSD3PoLvz61wF++iayoqi5IXU8xVuzLkjz3VkEwRqnFz8BCCC2ohIUDkqucmq0y750GljYZju\ndGab6ZNtINjr+Qv9f/eis9AavYLPoLKOEBQlZGTyvdskHpBO3BnbvUbQ57dFx7kc0NLiHAuj1qDR\nGIKihIywlo/uqSUuUrX3WWOf4qmqf4FLL4X+/QspLyu0/LWiKHkhU4profEzc7q6ehqlpTXAQr7B\nT/hVr0t44uUmmDkTHnkkGMFhxY+/KagNjSFYRXXaJWidHX3a1VJUNEBisbFZxxIK1Z7Zxj/86vQ9\n67mlRd6cPl1W7dE3HhyRNmNEJkwQeeklKxrzBToxTVF6FvE78VjsQYqKHqKtrY4VKy5g4sQpoa/5\nE3crNTVNoKlpQl4051zP6cMPnVHAAQcw4s47OeqLbWymD3P4DkcUD6DhRz+Cb34zUO2hw481CWoj\nIiMERQmSKJZ3KITmrK7Z1ibywgsikyaJFBW1jwje67enTOWHsjtfRKatU4GOEBRF2ZVJjgGkLbux\nbRvcfz8cdRR8+9vw619DURGccw68/DKXHv9tFnAiX1JaiI8QHvxYk0wbcACwFFgNvAVc4R4vB5qA\n94BGoCzFewOyn3aJil9RddolXzr9+sb96Mx1HkQuxfa8zJfoTkva17z7rsgVV4j0798+GpDBg0Vu\nvlnkk088647Kd5OwVjsFBgEj3cd9gXeBI4C5wAz3eA1we4r3BtNalonKl0R12iWfOv1MUMtVp19D\nlElzqnPPmTPHrpYdO0SefVakqkqS/EgiJ54o8uSTIq2tWeuOyncztAahy4XgGeAU4B1goCSMxjsp\nXhtAUymK4oUg4wCZzp2qQ85Ky+bNInfeKXLQQQkjUFoqcuGFIitXWtEfdvwahLzMVDbGDAViwKuu\nMdjgPrUBGJgPDYqi2CGIuj/pZmd74o034L774NFHaS/BOmyYM6nsggugvNy3vl2FwA2CMaYvsAi4\nUkT+YUxiEp2IiDEm5ZTkqVOnMnToUADKysoYOXIklZWVQKI2eaH348fCoifd/j333BPK9ttV2nPu\n3Lk89dRvKC/fm+rqaZSUlORFb/xYtu8/4IAyjLkMcX+ZxcVXc8op1wFOxz1hwjm0tk4DjmDZsinM\nmnUNxx13nKfzV1dP48UXz6G1dQ1wBKWlNQwbdgY33DA7KXW0mZaWqdTVzUv5+lNOuYbm55+ncvNm\nuO8+mpctc84PUFVF89ixcNxxVJ58srX2XLlyJVdddZW189nab25uZsGCBQDt/aUv/AwvutuA3YAG\n4KqkY+8Ag9zHg4mwyygqfkXVaZdsdBZyqchc2rO2tlaKivYSqBYYLUVFA6S2trb9eRvupFRBZa+u\npBcee8wJCA8enHAL9e8vcuWVTgA5IKLy3SSsMQTAAI8Ad3c6PheocR9fR4SDyorSHVGaS+BlbeWg\nPk9Gw9nW5swYPucckd12SxiC4cNF7r9f5B//8H39noJfgxDkPIQxwHnAScaYFe52KnA7MM4Y8x7w\nbXdfUQqGnzo4hcS27rq6ebS1HZrxNZ5y/nMgZZ2kb30LHnoIRo2CMWPg8cdh506YNAleeAFWrYJL\nLoG+fX1fX3HxY02C2ojICCEqw0jVmZ5cXDp+XUa1tbW+1zrwWsM/m7RV5+6/WiBx3qKivbqdR+CX\nLu25dq3IjBki5eWJ0UBFhcjMmSLr1/u+Xi5E5TdEWF1GvkSpQbCK6kxPLi6QbHUmd6C1tbVWYgpe\ndM+ZMyerayWMTOr4QbaftbtrxV83Z84cxy3U2OgUlDMmYQi+/nWRhQtFWlqy0mGbqPyG1CAoig/y\n7eO3dT0v58nlWn7u/r2OtpJf14/75erd+sk/9t8/YQSKi0XOP1/k1Vezur7i3yDoimnKLk119TSW\nLZvSnr7u+MQ95r8XkKB0+1nn2OvKY3V18xjSciWX8kemsJB+27fBR/+A/faDf/93uOgi2Gcfvx9F\nyQU/1iSojYiMEKIyjFSdmcnmrri+vl5GjToxZ/+5zTTU7nRn6zLyS7cjkh07RBYvluXle0vSi+Qu\nDpCbjz5eZPv2wLT5JSq/IdRlVDii8iVRnXZIdOY1vjpY20HZdGQKKgehIa2x27hR5PbbRQ48sN0I\nbAP5BSfJcGqluLhr4DpshP27GcevQdA1lRXFI2Fd6zgZL2UlglwbOfn6N084iROWL3fSRb/6ynnB\nIYfApZfyPwccwJz/fCyjTiV7/K6prDEERekhpKsHlMqH78XXnwtVJ51ElVtSgssvdw4aA6efDpdd\nBuPHQ1ERJwMnT57s+3qKXXSBHB8k14wJM6rTwe9ErsSkrOvoblJWISa7dV1Kcqq3pSSTyFn3xx/D\nTTfBgQfCv/0b/O//QlkZXHMNvPce/OY3cOqpzqI0nYjC9zMKGq3gx98U1IbGEKyiOu0Fc70ElQtV\nv6hrULcmZZppOn1Z625rE/n970XOOkukV69EoHjECJF580S2bfOkOwrfzyhoFNGgsqJ4Ip/zDQpV\nvyibDt3X2gPbtjkd/tFHJ4xAr14iZ5/tGIi2tiA/ppIBvwZBYwiK0kOI1wNKBJXTB4pzmm/wpz/B\nL37h1Bf6/HPn2D77wMUXO9t+++WsPYg1FpQc8GNNgtqIyAghKsNI1WnXjeNlDeBClbxOxmt7xkcL\nsdgYKS7eu6Pu3/1O5LnnRL7znY4lJU44QeTRR0W+/NK3zmzmS+QrZbczUfkNoS6jwhGVL4nqdPDb\nmcTfP2rUiZ4msAXVcXk9t5f27Gy8iovLJBYbK2dWflfWXHyxyCGHJIxASYnI1Kkir79u8dOIjBp1\noidXVdTWligEahAUJQey7bDDctdvW0fnuMFwamXJ/sNE+vRJGIIDD3Qmlm3caPGTpNeQziBEaW2J\nQqEGQVGyJJdONcjOKPsS1fZ0xGJjpBfHyyS+Ji+wpySdWOTkk0WeecYpOREgXv8fahC6x69B0HkI\nPohKbrLq7EjXfP05WebrN3f7Cq/5/PHJZE1NE2hqmsDEiVOszVvorj2XPvkk331zBWtZzSLe5SS2\n8g925/6i3Vk2bx48/zyccQb06mVFTzpKSkq6Lo6TIqgc1OI8XojKb8g3fqxJUBsRGSFExa+oOjuS\na1noxF1sTbeBT68jkGy1ZHPutO356qsi558vX5mi9tHAGgbJpdwr/dia9zvvbBcc0qByelCXkdLT\nsd0J5OqH96ojm07e5poFnY8n7zcuWeIsNPONb7QbgZ0gzxCTk/kXgYfVFdMDUIOg9GiCCuYGeaeZ\nTSdvcwZ1x2yhvaW4uEz2p05m8135G0kpo+XlIjNmSPOCBUkrpFUUPGCu+EcNQgGJyjAyyjrDGEi0\nPQ/BhnHq2E5tUkmNzKKf7CDhFnqv354iDz0k8sUXXa4di42RWGxs3l0xItH4fkZBo4h/g6AzlRXF\nMtnMGI6/3sbM3D58yXk8wGXcx1GsphloozdP8j3u4yD6Hv82jRdcEMi1lZ5BVushGGN6AX1E5O/B\nSdL1EJQEQdbu7zG89x7ra2ooe+a/2RPnd/MJRTzYq5gHdt7Op5Ttsu22q5XECHw9BGPM48DFwE7g\nj8Cexpj/EJG5uV5UUbwSv9u+/vrbWL/+I4YMOaTQksLBzp3w3HPOugMNDQxxD68qG8CSAw7mG7fd\nxOjevXmpbh4j6H6U0hPxuj6EkkR3PiXgDffvvwF1wG7AKj9+Kg/XtOlWC4yo+BWjrjPorCBbOvPC\n5s0id94pctBBiSBxaanIhReKrFjR4aWZ2jPo1M1srhFUe6aLP+Xy+aPyGyIPMYTexpjdgDOBn4vI\ndmOM+nOUvJHLCl897u7wjTec0cCjjzoNADBsGFx6KVxwAZSXezpNPtqloaGBCRO+T2vr4QC8+OL3\nWbLkiVC0/aZNG3rW98I23VkM4ArgY+A5nBXWhgJ/8GOFPFzTvulUIksumUa5ZicVauJTSlpbZeX1\n18ubZQMk6YPIawP2keU335xTSQm/WVte2icWG9MhjRUqJBYbk7VWv6QaWcZiY0OXtWYTgh4hiMjP\ngJ/F940x64GTrFsmRUlDdfU0li2b0n5j7JQsWGj9OqEZVXz6Kcybx5f/8R8cs2ULAFsp5mHgF/yU\n9zcPovT2GhYff3xetXltn/XrPwXuJD6ic47dkjedcVJle2W7pOguR3cWAxgEzAfq3f0jgR/6sUIe\nrmndcgZBVPyKYdKZ6Q4zk858VCf1evccSHu2tYm89JLIOeeI7LZb+2jgLfaVS/iF9OVfs76zTaXT\nz0Q4r+2T6i48Fhublc6gyPXzh+k3lAnyEENYADwM3ODuvw885RqJjBhjHgJOB/4mIiPcY7OAC4GN\n7suuF5F6DzqUiOPnDtxrvnxymuENN1zOiy8uAUKcZdPSAk884cQH/u//nGNFRTBpEteu38idy38A\nTAWe7/LWTZs2M378ZMB7SmW2cyRy4bbbrmfChPNpbXX2i4uv5bbbfmn1GrmSj88fabqzGMDr7t8V\nScdWerE2wIlAjKSsJOAnwDXdvC8A26kUmqBnHfstA9H5/UVFe0ltba01fR1Yu1ZkxgynjES8QSoq\nRGbOFFm/PoWeaoH+XUpT5LPchN81mwtNGDXZhqBLV+DU+h0QNwjAaOBFzxdwgtCdDUJ1N+8JprWU\nghK0QbBx/traWikqGiAwWqDabkfb1ibS2CgyYUKH5Sjf6V8mc4aPcgrQdSK5E5syZYqUlx8s5eUH\ny8EHH1mQ4GhUO9WwLHAUNPkwCKOAl4Gt7t/3gWM8XyC1QVgHvIHjdipL8Z7gWswiUfEr2tTpp0Po\n7kfpV6cNg+DlHFnr3LpV5N57Rb72tcRooLhYPv72t+XEknLPd9zJ6x0bU+aOGizq7IYozuuIay4v\nP9jXdyMqv3W/BsFLltFyY8xY4GvuoXdFZLtHj1Qq7gd+6j6+BWey2w87v2jq1KkMHToUgLKyMkaO\nHEllZSWQWKyi0PtxwqIn3f7KlSutnO+rr75yYwBTASe/fPjwYygqauPss7/LjBkzMr4/7r+94YbZ\nAMye7fhvbbVnIhtpDQClpQuorl6Y9flgDc7A2Nlfu/ZPfP3r36K8fG+qq6exZs0ab+cbOBB+/nOa\n58+HL790zrbffjRXVcHpp3PrA4/yh6/uAneecXyhnpKSki7nu+yyGlpb78DJ3GlG5GKMmYcTmltD\ncfE8qqsf73B9v+2ZvP/aa68xa9ZdbvxnDS++eA5Lljze4f9X6O9n5u/rmuQW6bDv5XwrV64Mze85\neb+5uZkFCxYAtPeXvkhnKYDJwKROf+OPJ3m1OHQaIXh5joiMEHY1Ot4910sYSyb7vYtNV0ba8+fc\nscNZdvKUUyTpdlQ2H320yNNPi2zf3v7SbEY0qe5w+/U7MG/umzBWne2OKHxfbUOAI4R/BbdSVmp+\nnYsBMsYMFpG/ursTgVW5nEcpNPNIzjX3Mns4H/it3tk5C2XTpsNYseIiuv2cmzbB/Plw//2wfj0A\n/wR+SSU/52T+9P59LO7bl6rezk+uoaGBTZs2U1RUTVvbKmBExvkVQ4YMYsuW6UlHpnPIIV+jsXFR\nzp9116IKmEJ5+S2MGnWMZhelw4816W4DHgc+AVqBD4EfAI8Ab+LEEJ4BBqZ4XzDm0zJR8Sva0tnx\n7nm09TtGPzqD8m+nujMeNerExAuWLxeZOlWkpCQxIjjkEPnFYUfLnvw8ZfukymaKxcZ0u2aCM1IZ\nLTBaiovL8lojKMigbFC/I5uao/JbJx8L5ADfBWYAN8U3Pxf1cD37LRUAUfmSBBFUjsXGdAhy2ugg\nctUZZGeV6txzb71V5NFHRU44IWEEjBE5/XSR554T2bkzo4sl1XPl5Qd3qzlbo6dBZXuao/JbD9wg\nAP/p3tV/hJMh9BYw389FPVwzgKZSbBPEWse5nC8f8xvGjZsk3z/xNPng3HNFBg5MGIKyMpFrrhF5\n//0u70lnpFLphdE91q+t5I98GIRV7t833b99gWV+LurhmtYbSgk3+SipkBNtbSK//73IWWeJ9OqV\nMAQjRojMmyeybVvGz5TKwHX+rDDQDXqGP1CrhJt8GITX3L+vAPsBuwMf+Lmoh2vab6kAiMowMgo6\nnU69JqdOPRCX0bZtTod/9NEJI9Crl8hZZ8nSe+5xDIUP6uvr3cyh0a4xsD+yicL/XSQaOqOgUcS/\nQSjyEHd+1hizF3AHsBxY6waLlR5OQ0MD48dPZvz4yTQ0NBRaTlrimUHjxi1h3Lgl/iqU/ulPMH06\n7L8/TJsGb74J++wDN97oZA899RQccwyYnFcpbNf82GM/p7R0LfApsNDNMprm67yK4ot0lgI4Dhic\ntD8FaALuBQb4sULdbURkhNCTyfdU/4KWFti5U16vrZX/rRgkO5PmDsjo0U7w+MsvA7t0VEpBZHJ/\nRUH/rgJBuYyAFUC5+/hbwF9xJqbVAk/7uWi3otQgFBw/fvlcO4m8dy6ffSZy992ybd99241AC73l\nkV67y8s/+1nw148I6Yz1rlIfKEoEaRDeSHr8c2BWqueC2KJiEKLiV8xFp58Vx3LtJPI2D2HVKpFL\nLpHtu+/ebgjWMUBmcLsMYGO3n7Un/99Tke67kOk7EoY1lW0SBY0i/g1CppnKvYwxu4lTt+gUINm5\n6WUdBSXC5LpKWS7rH/vF0zoLO3bAf/+3s+6AWwumN/A8R3IfO3iW62jjgsA07kqEZuU5JXvSWQqc\nBXFeBpbguI+K3OOHAi/5sULdbURkhNDTycWFk+ukKz9kHM1s2CBSWyuy//7to4F/9uotD5b0lyOY\nKmGscxOE6yyIKrXpjkex7lFPgSDTToETcOoN9Uk6dhhwrJ+LditKDUJk6ZpjXyHW1xXoRKoO6LLj\nKkXOP1+kuLjdEGzbbz+5erd+0o/7k7TFUz6rpbz84IIER5M769raWut+eRu+/myCykEZhKBiTD0p\nMB6oQSjUFhWDEBW/Yr7XQ8g1x95v6YoS5sn5XCR/NL3bjYAY4yxI09go40+Z2KWjitdkCirW0V17\npapr1N06B16J68z3HXu2BshLewYVwPZ63qj81v0aBI0FKClJXps4vl6vV99wVVUVo0YdQ1PTBJwq\nk8FSdeSRrDrjJAYsupyy7V85NXrLy+HCC+GSS2DYMADkjge6vLe8fCOjRi0JpPqll/bqHHNpawPo\nqjMqxL83hx9+CPAwFRUDrLRtULGpQsS8Qo0faxLURkRGCD0VG75hm3d0Ke+y29pEXnhBZNIkkaKi\nxIhg5EiR+fNFvvgiUE1e8NJeqV7jLOEZLpdR9teplqKiARKLjbVyraBGOT0t3oG6jBTb5JJmmAob\nvtnOnVnF7nvLW5dfLjJ8eMII9O4tcs45Ii+91G1JiXxOsPLSXqk669raWutaamtr29djrq2tTfs6\nP+2Q+Lz14tRnCr9R62lzKdQgFJCo+BWz1ZmuIwv6x5NKZ1zLobwrd3OlfE5pwhAMGiQya5bIJ5/4\num4QPu9szht0WelsdPj5/ya+N9ndOGTTnoUKKkflt64GoYBE5UuSrc5MHUN9fb3EYmOlvPzgbhd1\n8a1zxw758cgT5DmOkqTeRV7uvbvUjjhObp01y0oHke3Ix2ZQOUiyDSr7nZ0ei4113V2HW23PMGQB\nReW3rgZByRqvmULpXCvJC+MUF+9t/0e6ebPInXeKHHRQuxH4J8Uyj2/JMfQRJwunWqC/ldFKT/Mj\ndyZog9D5BsKYMjGmr5X/TS6jljAYkEKhBkFJSaYO3Y9bIBYb26XTiMXG2hG9cqXIhReKlCa5hYYN\nk3cuvFAmVX7XTWWNp2R27bxisbE511AK2o8c/3/EYmPadQZ1jVz/551fV1y8t6s3s+ZUhiTX/4WX\nc3cXt+pJMYFsUYNQQMI6jOz6w97L2qQhp1PuOhM5Z1pbRZ58UuTEE2VpkltIqqpk+c03S9UpE9s7\nlY7a7WbnZHNXmbsLrlqCmhGdqiOcM2dOh+e9fL6E4RorvXv3EThKoNzVnlqz3+9UpvbM9txBjfbC\n+lvvjBqEAhLWL0nXH0VN+4/C7w8mFhvToVODConFxmSt8YXHHpMFBx0hm4oTBeaWlpaKXHGFyLvv\nps2+SU5rTHYZ2ZzQ1R25B+mDc02l+r/267dvzgan6/85/apufu/KM7VntudWg6AGQelEph9FNj+w\nVHeVtbUysm0cAAAV+UlEQVS1Ykw/gf0FjpLi4jLvP/62NpGXXpJPKivlq6TRwGrTS1ZfdpnI3//e\n7WdI1pScnul0YPlNifVKoQyCn3WaU40E4/pTaQ6yPbM5t7qM1CAonejuR+E1qJz5Dt25K8+U097O\nF1+IPPSQyLHHthuBHRhZxESp5AWBhz1N2LLlO853p1EIl5HfdZpTxYr8GJh8okFlNQgFIczDyOQf\nRbIv2SupOuRUd40ZO5u1a0VmzBApL283BFJRIY8NPUwOoK7LeZLbMxejlqkjSH4uVWdny+edjnwF\nlTvWkFqas0HonE1mTJn1NOM4Yf4dxYmCRhE1CAWlEF8SJ997jPTtO1j69TvQU2kAWwvk9O072O1s\nJqW/+2xrE2lqEjnjjI4lJb7+dZGFC0VaWtJ29p112sqU6vz6RPA5fwYhX3T8rDX2S4YEQJjbM04U\nNIqoQdilcO7aysQJpiZcD0HMBeiaqVQmvXsnsnigomP8YOtWkXvvFTn88IQRKC4WOe88kVdfTXn+\nVLEAL5/Df+ZJtRuE7pl+5l3ZZbKr49cgaLXTCFFXN4/W1sPdvUuIV2hsbbVfobGqqorFixe2Vzzd\ntOkYVqy4oP2aAMOHP0zVkCFw2WWwcCFs2+Y8sd9+vH/KKdywdiMfrP4QLrmGioqB7VVT4+fPpoKq\nXUZwzDFHUlGxBCCQSqeKEkn8WJOgNiIyQsj3MNK50x0t8Rr+mTJwku8MbaxVnBw/KGKHTOAKWV6+\njySJEBk7VuTpp6Xht7/1vEhOx7v3pQILJBYbk/EO16/LyO+IoHOsI0x34zZdRvkiCu6YKGgUUZdR\nQcn3lySTyyjTSlt+F56J5/0PoK/M4CxZy4CEEdhjD5GLLxZ5883296VOgUydstjVIHhz52TbEXf3\neic2461GU7ZF4/JJKgMb9jIcUehso6BRJOQGAXgI2ACsSjpWDjQB7wGNQFmK9wXSWD2BdEHlICbk\nxM8ZY7nM5wJpoVe7Idi2774id98t8tlnad/nxSDYDvjmQueMmi7xkTQENQnKD2HUpOQPvwYh6BjC\nw8C9wCNJx64DmkRkrjGmxt2/LmAdPYa47z2ZhoYGli9/A3gbp8kHAMP8Xai1lZP++iGzqOWbfABA\nG4ZXKwZy/CMP06eqCoqKUr61unoay5ZNoaUlfmQ6MIXS0hqqqxd2+TwdYxVHsWKFP+nZ4sRm7iA5\nPtLa+kCXuEznVeTCSOe2T9XmipIWP9bEywYMpeMI4R1goPt4EPBOivfYN50BEIZhZLpJT9C/fdJY\nVjo//ljkpptEBg5sHw1sYQ+5kyoZXjJA6uvrPU9si+f8x2MCXtJjC+GGSTfLN/nOOpVvPpObrpDE\n237UqBNDoac7wvA76o4oaBQJuctIUhuEz5Iem+T9pOP2WyoAwvAl8VIWoVudbW0iv/+9yNlnO6uP\nxU8yYoS8deWV8q/fntDeoQfZYSf75vMZqPXiMkrnmw9bUDmZMHw/vRAFnVHQKBJ+l1FGRESMMZLq\nualTpzJ06FAAysrKGDlyJJWVlQA0NzcD6L67D78DXgEmuPvNwBqSaW5u7vL+HVu3sur6m6n4YA0H\nbP+SSqCtqIhnRozgia8Mnw88mOrTTuOaM0varzd+/GRaWqYCQ4BKWlrghhtmU1JS4klvQ0MDN9ww\nG4CjjjqIZ59dxvbtLZx11qnMnz8fgJKSEmbOvDxv7VdSUsItt0zniSceZv36jygvH8RFF53f7i5q\nbm5my5aNdMRp36qqKkpKSjqcr9Dfh6jtx4+FRU+6/WStYdBTWVlJc3MzCxYsAGjvL33hx5p42Ujt\nMhrkPh5MhF1Gtsm2HIOIU2zOyTrKok7OBx/InydPli2Y9tHApxj5KaNl2G79O9wtdz6Pn6Bl13o7\n/QUmd3FxhZEwZhQpSmeIoMtoLlDjPr4OuD3FewJoKvvYHEZ6KSbXfR5/vcBoKS8/uOs8hJ07RZ57\nTuT000VMwhC8zGg5l19JMQ+6bqfUcxySdea6YlpqX/0AiZfB6NdvXyttGRTqmw+GKOiMgkaRkLuM\njDGPA2OBCmPMh8BNwO3AU8aYHwLrgLOD1BA2OmeqxGfrnnvupe5sXSfTpaUF7rrrli7HMs9IrgI+\nZdSoJYnXbN0KTz8NF10EHzjZQpSU0DBgIDd8cgHLmeW+N1UmSgPwAMuXb6ShoSHputuBB5Ie++FQ\nYB4Jd1d4iWd4Jbs3FKVH4ceaBLURkRFCtmQeBXS9M/dSXTStK2PVKpFLLhHp06d9NCAHHihy220i\nGzemcN84s4mLi8vcEUBqF5Rfl1HypDOnRHO1+9nD7TJSlChA2F1GOYnqoQYhc0npereDzM5lJJJw\nZVSdfKas+PGPRSorJekiIiefLLJ4scj27Snf1zk1NFFGuWvH73fiU21trTv5bLRrDPaUvn0HqzFQ\nFAuoQSgguS+lmG4UUC9wlPTuvU97+QRPaY0bNojMni2y//4JI9Cnj8iPfiSyerW18tdxDZ2roGZb\n3z+Imkv5RHXaJQo6o6BRRA1CQcn2S9L57ri7dYK7zWR59VWR8893ykzHe+7DDhP52c9EPv88pU6v\nefOZsmqSF3vJlJGU7TWj8qNTnXaJgs4oaBRRgxAZutbsSSw/maqiaFp3TEuLs9DMN76RMALGiEyY\nINLY6GQTedTQXQfe3epeXtxHmq6pKPlDDUJE6K7zzOS3FxGRv/xFZOZMkb33ThiCvfYSufZakT//\nOeU1O9+Z2+7AvZxPi60pSv7waxBSVydTPNF5BmOuxBeJ2bLlTJxCcAuBhZTuPoPZpxwPkyfD0KFw\n662wcSOMHAnz58NHH8HcuTCsayG7+DmbmibQ1HQoEydOYdOmzd1qqaubl5Tq6ixcE0+T7Ux19TRK\nS2sSektrfBV9s9WeQaM67RIFnVHQaANdMS1PZKpC2bETHkcfbuTiPdYya+8+9Kupcd7Quzd873vO\n6mQnnADGZLxex3M209JyBPAgpaU11iphdq5UmmrlMa2+qSgRws/wIqiNHugyEkmsZeAswjK23RUT\nd6scyrtyN1fK55Qm3EKDBonMmiXyySdZXStdRlN3cYEgfP5hLgCnKD0JfLqMjHOOcGGMkTDq8kvn\n9YNLS2tY/PRDVLz+Optvns34ttb213525JHsddNNMHEiFBf7vhZcAVwEjHCum2Hd4lSzqRVFCT/G\nGEQks/sgE36sSVAbERkh+JmHsBeb5Rq+Jx+XJmYSf1lUJL/bd4i8dN99VvTF78z79dvXTWkNd2A3\nKql9qtMuUdAZBY0iGlSOHEfzF+ZxER+xL3U8yb4t/+SLgQPhjjso2biR0z5exzcvvTTjORoaGhg/\nfjLjx0+moaEh7euqqqpobFzEYYcdDIyw/EkURelx+LEmQW1EZITgmdZWkSeflM1HHSVJt+lSz1Fy\nOldJn933yWqWb7Y+fp0LoCi7Bug8hPDywmOPyYKDjpCNJbu3G4Hte+whD5TsKYdxW8Y5CemCsLnm\n9WtgV1F6Pn4NgrqMfJAyN1mEV+6+m9+V782Yc89lyp/XUPHVl7xtevH2ZZfR+9NPWfStk3mPwR3e\ntmnTZsaPn8yxx1YyYcL33fkDE5g4cUpGt5BXnXH3UWPjotAGiaOS66067RIFnVHQaAOdh2CLlhZ4\n4gm23noro911B3Zi+DVnch+XsVT+wrh3n6WxX78uufnFxdeyevV2WlvvcU82HRgEVHVZA0Hz+hVF\nCQpNO/XLunVw//3wX/8FW7YAsJF+PMhgHuBiPuQa94ULGTduCY2Ni2hoaOD6629h/fpPGTJkf2AH\nK1ZcRHwhHGfm7xJgUYf3xdG0UEVRUuE37VRHCLnywQcwfTo8+yy0tTnHRo3iji/hxtUX8xUH4nTw\nA4DEnXzn+QEtLTUcfvghKS7wCbCQ4uJrqa7+ZYdn4it3KYqi2ERjCLnSpw/Nzz4LvXrBeefBK6/A\nH//I0XWzKSq9EfgUOI+iompisYfbJ4I5JSXOwxkBLHEf9+5QEwiuBXoBDzB8+GG+O/+o+D9Vp11U\npz2ioNEGOkLIlcGDYdYsmDYNBg7s8NThhx/O+vW3MGTIIG677dEOHfqmTRuA3wN3ukem8/e/D2bx\n4oWce+6lbNmyN/BLnPWRF1JRsSQvH0dRFEVjCBZJWZqiU4mIY4+tZMWKC0iOFxhzNc899zhAt+9X\nFEVJh8YQQkTHCqN0yRACqKgY0OV9Il+jrm4ejY2Luq0eqiiKEhQaQ/BBLn7F6uppFBVdTSJeUAOM\naX8+iPkCUfF/qk67qE57REGjDdQgWMTLgjFVVVX89KfVFBVVAw8A51Fa+itfC8soiqLYQGMIlulu\njkD8eWf1sh1UVAzUuQSKolhBYwgRInXQ+UY1BoqihAJ1Gfmgs1+x4zrGTh2i2bNnt5eqvv762zyv\nVxykzrCiOu2iOu0RBY020BGCRTp2+E6W0U03VdPWVgfgxg0URVHCicYQLDF79mx+/OO7gLvoWJPo\nTuAwd99QVPQCbW13AzrPQFEUu/iNIRTMIBhj1gF/B3YC20XkuKTnImUQGhoa+M53zqGt7QckjADA\nlYAAP3P3p3Pwwftw0EGHA1qYTlEUu/g1CIWMIQhQKSKxZGMQJeJ+xbq6ee5d/53Ar3DSSX9MaWkx\njjGY4m530r//3nlflyAq/k/VaRfVaY8oaLRBoYPKOVuy8FIFXEJRUQuHH35kl2dTzVRWFEUJA4V0\nGf0Z2IrjMvpPEXkw6bnIuYyS00mNuYqDDtqX/v33YvXq92htvQPQmIGiKMES5RjCYBH5qzFmb6AJ\nuFxE/uA+FymDAB0nnK1e/Ub76mfFxVcxfPgxVFQM0JiBoiiBEtmJaSLyV/fvRmPMYuA44A/x56dO\nncrQoUMBKCsrY+TIkVRWVgIJf16h9+PHmpubKSkpobFxEePHT6a19WJgCFBJaysUFc1n5szLC6b3\nnnvuCWX7ZWrPMOjR9szPfhTac+XKlVx11VWh0RPfb25uZsGCBQDt/aUvRCTvG7AH0M993Ad4CRif\n9LxEgaVLl3Y5Nm7cJIEFAuJuC2TcuEn5F5dEKp1hRHXaRXXaIwoaRUTcvjPnvrkgLiNjzDBgsbvb\nG3hURG5Lel4KocsGXtZEUBRFCYLIxhAyEWWDAN0XuFMURQmCKM9DiDzJvtpkgljTwA/pdIYN1WkX\n1WmPKGi0gRoERVEUBVCXkaIoSo9BXUaKoiiKFdQg+CAqfkXVaRfVaZco6IyCRhuoQVAURVEAjSFY\nQdNMFUUJAzoPocDoRDRFUcKCBpULSHNzM3V18wqyTnI2RMX/qTrtojrtEQWNNlCDoCiKogDqMvKN\nuowURQkLGkMIARpUVhQlDGgMoYDE/Yphq13Umaj4P1WnXVSnPaKg0QZqEBRFURRAXUaKoig9BnUZ\nKYqiKFZQg+CDqPgVVaddVKddoqAzChptoAZBURRFATSGkDOaaqooStjQeQgFQCejKYoSRjSoXAAS\n9YuGENb6RclExf+pOu2iOu0RBY02UIOgKIqiAOoyygl1GSmKEkY0hlAgNKisKErY0BhCgaiqqmLm\nzMtDW78omaj4P1WnXVSnPaKg0QZqEBRFURRAXUaKoig9BnUZKYqiKFYoiEEwxpxqjHnHGPO+Maam\nEBpsEBW/ouq0i+q0SxR0RkGjDfJuEIwxvYD7gFOBI4FzjDFH5FuHDVauXFloCZ5QnXZRnXaJgs4o\naLRBIUYIxwEfiMg6EdkOPAGcUQAdvvn8888LLcETqtMuqtMuUdAZBY02KIRB2A/4MGn/I/eYoiiK\nUkAKYRB6TPrQunXrCi3BE6rTLqrTLlHQGQWNNsh72qkxZjQwS0ROdfevB9pEZE7Sa3qM0VAURckn\nkSpdYYzpDbwLnAx8ArwGnCMia/IqRFEURelA73xfUER2GGMuAxqAXsB8NQaKoiiFJ5QzlRVFUZT8\nE6qZysaYWcaYj4wxK9zttKTnrncnsr1jjBlfSJ2untBOrjPGrDPGvOm24WvusXJjTJMx5j1jTKMx\npizPmh4yxmwwxqxKOpZWU6H+32l0hu57aYw5wBiz1Biz2hjzljHmCvd4qNo0g85QtakxZndjzKvG\nmJXGmLeNMbe5x0PTnhk02mtLEQnNBvwEuCbF8SOBlcBuwFDgA6CogDp7uRqGuppWAkcUuv2S9K0F\nyjsdmwvMcB/XALfnWdOJQAxY1Z2mQv6/0+gM3fcSGASMdB/3xYnLHRG2Ns2gM4xtuof7tzfwCvAv\nIWzPVBqttWWoRgguqSLkZwCPi8h2EVmH88GOy6uqjkRhcl3ndpwALHQfLwTOzKcYEfkD8JlHTQX7\nf6fRCSH7XorIpyKy0n28DViDM58nVG2aQSeEr02/cB8W49z0fUb42jOVRrDUlmE0CJcbY94wxsxP\nGp7tizOBLU6hJ7OFfXKdAM8bY143xlzkHhsoIhvcxxuAgYWR1oF0msL2/4YQfy+NMUNxRjWvEuI2\nTdL5insoVG1qjCkyxqzEabelIrKakLVnGo1gqS0LUcuoyRizKsU2AbgfGAaMBP4K1GU4VSGj4WGP\nxI8RkRhwGnCpMebE5CfFGU+G6jN40FRIvaH9Xhpj+gKLgCtF5B8dhISoTV2dT+Po3EYI21RE2kRk\nJLA/8C1jzEmdni94e6bQWInFtixE2uk4L68zxvwX8Ky7+zFwQNLT+7vHCkVnPQfQ0RIXFBH5q/t3\nozFmMc4wcYMxZpCIfGqMGQz8raAiHdJpCtX/W0Ta2ypM30tjzG44xuCXIvKMezh0bZqk81dxnWFt\nUwAR2WqM+S0wihC2ZyeNXxeR5vhxv20ZKpeR2+BxJgLxTI8lwPeNMcXGmGHAoTgT2grF68Chxpih\nxphi4HuuxoJjjNnDGNPPfdwHGI/TjkuAKe7LpgDPpD5DXkmnKVT/7zB+L40xBpgPvC0i9yQ9Fao2\nTaczbG1qjKmIu1qMMaXAOGAFIWrPdBqNMYOSXuavLYOOimcZQX8EeBN4A6fhByY9NxMnKPIOUBUC\nrafhZEx8AFxfaD1JuobhZBasBN6KawPKgeeB94BGoCzPuh7HmZneihN/uSCTpkL9v1Po/EEYv5c4\n2SVt7v95hbudGrY2TaPztLC1KTAC+D9X55vAte7x0LRnBo3W2lInpimKoihAyFxGiqIoSuFQg6Ao\niqIAahAURVEUFzUIiqIoCqAGQVEURXFRg6AoiqIAahAUJSXGmAFJ5YT/mlReuK1zGWFjzFXGmF8U\nSqui2EINgqKkQEQ2i0hMnJpQDwB3uY8vBr7f6eXfAx7Lt0ZFsY0aBEXxRry88CLgdOOsDR6v4Lmv\niCwrkC5FsYYaBEXJAhHZglMP5jvuoe8DTxZOkaLYQw2ComTP4yTcRt9z9xUl8qhBUJTsWQKcbIyJ\n4SxpuKLQghTFBmoQFCVLxFngZSnwMBpMVnoQahAUxRudywI/jlOOWN1FSo9By18riqIogI4QFEVR\nFBc1CIqiKAqgBkFRFEVxUYOgKIqiAGoQFEVRFBc1CIqiKAqgBkFRFEVxUYOgKIqiAPD/AUd7480b\nVee2AAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# create a DataFrame with the minimum and maximum values of TV\n",
"X_new = pd.DataFrame({'TV': [data.TV.min(), data.TV.max()]})\n",
"#X_new.head()\n",
"# make predictions for those x values and store them\n",
"preds = lm.predict(X_new)\n",
"#preds\n",
"# first, plot the observed data\n",
"data.plot(kind='scatter', x='TV', y='Sales')\n",
"\n",
"# then, plot the least squares line\n",
"plt.plot(X_new, preds, c='red', linewidth=2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Confidence intervals\n",
"\n",
"Statsmodels calculates 95% confidence intervals on the model coefficients. It means that if we assume that the model $Y=\\beta_0+\\beta_1 X +\\xi$ holds then with probability tending to 95%, the \"true\" coefficients $\\beta_0$ and $\\beta_1$ belong to these intervals."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"
\n",
" \n",
" \n",
" | \n",
" 0 | \n",
" 1 | \n",
"
\n",
" \n",
" \n",
" \n",
" Intercept | \n",
" 6.129719 | \n",
" 7.935468 | \n",
"
\n",
" \n",
" TV | \n",
" 0.042231 | \n",
" 0.052843 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" 0 1\n",
"Intercept 6.129719 7.935468\n",
"TV 0.042231 0.052843"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# print the confidence intervals for the model coefficients\n",
"lm.conf_int()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Hypothesis Testing and p-values\n",
"\n",
"We want to study the following test problem:\n",
"- **null hypothesis:** There is no relationship between TV ads and Sales $$H_0 : \\beta_1 = 0$$\n",
"- **alternative hypothesis:** There is a relationship between TV ads and Sales $$H_1 : \\beta_1\\neq0$$\n",
"\n",
"How do we test this hypothesis? Intuitively, we reject the null (and thus believe the alternative) if the 95% confidence interval **does not include zero**. Conversely, the **p-value** represents the probability that the coefficient is actually zero:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Intercept 1.406300e-35\n",
"TV 1.467390e-42\n",
"dtype: float64"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# print the p-values for the model coefficients\n",
"lm.pvalues"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If the 95% confidence interval **includes zero**, the p-value for that coefficient will be **greater than 0.05**. If the 95% confidence interval **does not include zero**, the p-value will be **less than 0.05**. Thus, a p-value less than 0.05 is one way to decide whether there is likely a relationship between the feature and the response. (Again, using 0.05 as the cutoff is just a convention.)\n",
"\n",
"In this case, the p-value for TV is far less than 0.05, and so we have a high confidence in rejecting ($H_0$). As a consequence, we highly believe that there is a strong connection between TV ads and Sales.\n",
"\n",
"(Note that we generally ignore the p-value for the intercept)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## How Well Does the Model Fit the data?\n",
"\n",
"The most common way to evaluate the overall fit of a linear model is by the **R-squared** value. R-squared is the **proportion of variance explained**, meaning the proportion of variance in the observed data that is explained by the model, or the reduction in error over the **null model**. (The null model just predicts the mean of the observed response, and thus it has an intercept and no slope.)\n",
"\n",
"$$R^2=\\frac{\\|\\hat y - \\bar{Y}_n\\|_2^2}{\\|Y-\\bar{Y}_n\\|_2^2}$$\n",
"where $\\hat y$ is the prediction vector via least squares.\n",
"\n",
"R-squared is between 0 and 1, and higher is better because it means that more variance is explained by the model. Here's an example of what R-squared \"looks like\":"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can see that the **blue line** explains some of the variance in the data (R-squared=0.54), the **green line** explains more of the variance (R-squared=0.64), and the **red line** fits the training data even further (R-squared=0.66). (Does the red line look like it's overfitting?)\n",
"\n",
"Let's calculate the R-squared value for our simple linear model:"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0.61187505085007099"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# print the R-squared value for the model\n",
"lm.rsquared"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Is that a \"good\" R-squared value? It's hard to say. The threshold for a good R-squared value depends widely on the domain. Therefore, it's most useful as a tool for **comparing different models**."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Multiple Linear Regression\n",
"\n",
"Simple linear regression can easily be extended to include multiple features. This is called **multiple linear regression**:\n",
"\n",
"$$Y = \\beta_0 + \\beta_1X^{(1)} + ... + \\beta_n X^{(k)} + \\xi$$\n",
"\n",
"Each $X^{(j)}$ represents a different feature ( = coordinates of the design $X=(X^{(1)}, \\ldots, X^{(k)})$, and each feature has its own coefficient. In this case:\n",
"\n",
"$$\\boxed{Y = \\beta_0 + \\beta_1 \\times TV + \\beta_2 \\times Radio + \\beta_3 \\times Newspaper +\\xi}$$\n",
"\n",
"Let's use Statsmodels to estimate these coefficients:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Intercept 2.938889\n",
"TV 0.045765\n",
"Radio 0.188530\n",
"Newspaper -0.001037\n",
"dtype: float64"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# create a fitted model with all three features\n",
"lm = smf.ols(formula='Sales ~ TV + Radio + Newspaper', data=data).fit()\n",
"\n",
"# print the coefficients\n",
"lm.params"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"How do we interpret these coefficients? For a given amount of Radio and Newspaper ad spending, an **increase of $1000 in TV ad spending** is associated with an **increase in Sales of 45.765 widgets**.\n",
"\n",
"A lot of the information we have been reviewing piece-by-piece is available in the model summary output:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"OLS Regression Results\n",
"\n",
" Dep. Variable: | Sales | R-squared: | 0.897 | \n",
"
\n",
"\n",
" Model: | OLS | Adj. R-squared: | 0.896 | \n",
"
\n",
"\n",
" Method: | Least Squares | F-statistic: | 570.3 | \n",
"
\n",
"\n",
" Date: | Sun, 27 Sep 2015 | Prob (F-statistic): | 1.58e-96 | \n",
"
\n",
"\n",
" Time: | 20:32:01 | Log-Likelihood: | -386.18 | \n",
"
\n",
"\n",
" No. Observations: | 200 | AIC: | 780.4 | \n",
"
\n",
"\n",
" Df Residuals: | 196 | BIC: | 793.6 | \n",
"
\n",
"\n",
" Df Model: | 3 | | | \n",
"
\n",
"\n",
" Covariance Type: | nonrobust | | | \n",
"
\n",
"
\n",
"\n",
"\n",
" | coef | std err | t | P>|t| | [95.0% Conf. Int.] | \n",
"
\n",
"\n",
" Intercept | 2.9389 | 0.312 | 9.422 | 0.000 | 2.324 3.554 | \n",
"
\n",
"\n",
" TV | 0.0458 | 0.001 | 32.809 | 0.000 | 0.043 0.049 | \n",
"
\n",
"\n",
" Radio | 0.1885 | 0.009 | 21.893 | 0.000 | 0.172 0.206 | \n",
"
\n",
"\n",
" Newspaper | -0.0010 | 0.006 | -0.177 | 0.860 | -0.013 0.011 | \n",
"
\n",
"
\n",
"\n",
"\n",
" Omnibus: | 60.414 | Durbin-Watson: | 2.084 | \n",
"
\n",
"\n",
" Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 151.241 | \n",
"
\n",
"\n",
" Skew: | -1.327 | Prob(JB): | 1.44e-33 | \n",
"
\n",
"\n",
" Kurtosis: | 6.332 | Cond. No. | 454. | \n",
"
\n",
"
"
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Sales R-squared: 0.897\n",
"Model: OLS Adj. R-squared: 0.896\n",
"Method: Least Squares F-statistic: 570.3\n",
"Date: Sun, 27 Sep 2015 Prob (F-statistic): 1.58e-96\n",
"Time: 20:32:01 Log-Likelihood: -386.18\n",
"No. Observations: 200 AIC: 780.4\n",
"Df Residuals: 196 BIC: 793.6\n",
"Df Model: 3 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"------------------------------------------------------------------------------\n",
"Intercept 2.9389 0.312 9.422 0.000 2.324 3.554\n",
"TV 0.0458 0.001 32.809 0.000 0.043 0.049\n",
"Radio 0.1885 0.009 21.893 0.000 0.172 0.206\n",
"Newspaper -0.0010 0.006 -0.177 0.860 -0.013 0.011\n",
"==============================================================================\n",
"Omnibus: 60.414 Durbin-Watson: 2.084\n",
"Prob(Omnibus): 0.000 Jarque-Bera (JB): 151.241\n",
"Skew: -1.327 Prob(JB): 1.44e-33\n",
"Kurtosis: 6.332 Cond. No. 454.\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# print a summary of the fitted model\n",
"lm.summary()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What are a few key things we learn from this output?\n",
"\n",
"- TV and Radio have significant **p-values**, whereas Newspaper does not. Thus we reject the null hypothesis for TV and Radio (that there is no association between those features and Sales), and fail to reject the null hypothesis for Newspaper.\n",
"- TV and Radio ad spending are both **positively associated** with Sales, whereas Newspaper ad spending is **slightly negatively associated** with Sales. (However, this is irrelevant since we have failed to reject the null hypothesis for Newspaper.)\n",
"- This model has a higher **R-squared** (0.897) than the previous model, which means that this model provides a better fit to the data than a model that only includes TV."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Feature Selection\n",
"\n",
"How do I decide **which features to include** in a linear model? Here's one idea:\n",
"- Try different models, and only keep predictors in the model if they have small p-values.\n",
"- Check whether the R-squared value goes up when you add new predictors.\n",
"\n",
"What are the **drawbacks** to this approach?\n",
"- Linear models rely upon a lot of **assumptions** (such as the features being independent), and if those assumptions are violated (which they usually are), R-squared and p-values are less reliable.\n",
"- Using a p-value cutoff of 0.05 means that if you add 100 predictors to a model that are **pure noise**, 5 of them (on average) will still be counted as significant.\n",
"- R-squared is susceptible to **overfitting**, and thus there is no guarantee that a model with a high R-squared value will generalize. Below is an example:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0.89719426108289568"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# only include TV and Radio in the model\n",
"lm = smf.ols(formula='Sales ~ TV + Radio', data=data).fit()\n",
"lm.rsquared"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0.89721063817895219"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# add Newspaper to the model (which we believe has no association with Sales)\n",
"lm = smf.ols(formula='Sales ~ TV + Radio + Newspaper', data=data).fit()\n",
"lm.rsquared"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**R-squared will always increase as you add more features to the model**, even if they are unrelated to the response. Thus, selecting the model with the highest R-squared is not a reliable approach for choosing the best linear model.\n",
"\n",
"There is alternative to R-squared called **adjusted R-squared** that penalizes model complexity (to control for overfitting), but it generally [under-penalizes complexity](http://scott.fortmann-roe.com/docs/MeasuringError.html).\n",
"\n",
"So is there a better approach to feature selection? **Cross-validation.** It provides a more reliable estimate of out-of-sample error, and thus is a better way to choose which of your models will best **generalize** to out-of-sample data. There is extensive functionality for cross-validation in scikit-learn, including automated methods for searching different sets of parameters and different models. Importantly, cross-validation can be applied to any model, whereas the methods described above only apply to linear models."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"OLS Regression Results\n",
"\n",
" Dep. Variable: | Sales | R-squared: | 0.897 | \n",
"
\n",
"\n",
" Model: | OLS | Adj. R-squared: | 0.896 | \n",
"
\n",
"\n",
" Method: | Least Squares | F-statistic: | 859.6 | \n",
"
\n",
"\n",
" Date: | Fri, 25 Sep 2015 | Prob (F-statistic): | 4.83e-98 | \n",
"
\n",
"\n",
" Time: | 15:39:26 | Log-Likelihood: | -386.20 | \n",
"
\n",
"\n",
" No. Observations: | 200 | AIC: | 778.4 | \n",
"
\n",
"\n",
" Df Residuals: | 197 | BIC: | 788.3 | \n",
"
\n",
"\n",
" Df Model: | 2 | | | \n",
"
\n",
"\n",
" Covariance Type: | nonrobust | | | \n",
"
\n",
"
\n",
"\n",
"\n",
" | coef | std err | t | P>|t| | [95.0% Conf. Int.] | \n",
"
\n",
"\n",
" Intercept | 2.9211 | 0.294 | 9.919 | 0.000 | 2.340 3.502 | \n",
"
\n",
"\n",
" TV | 0.0458 | 0.001 | 32.909 | 0.000 | 0.043 0.048 | \n",
"
\n",
"\n",
" Radio | 0.1880 | 0.008 | 23.382 | 0.000 | 0.172 0.204 | \n",
"
\n",
"
\n",
"\n",
"\n",
" Omnibus: | 60.022 | Durbin-Watson: | 2.081 | \n",
"
\n",
"\n",
" Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 148.679 | \n",
"
\n",
"\n",
" Skew: | -1.323 | Prob(JB): | 5.19e-33 | \n",
"
\n",
"\n",
" Kurtosis: | 6.292 | Cond. No. | 425. | \n",
"
\n",
"
"
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Sales R-squared: 0.897\n",
"Model: OLS Adj. R-squared: 0.896\n",
"Method: Least Squares F-statistic: 859.6\n",
"Date: Fri, 25 Sep 2015 Prob (F-statistic): 4.83e-98\n",
"Time: 15:39:26 Log-Likelihood: -386.20\n",
"No. Observations: 200 AIC: 778.4\n",
"Df Residuals: 197 BIC: 788.3\n",
"Df Model: 2 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"------------------------------------------------------------------------------\n",
"Intercept 2.9211 0.294 9.919 0.000 2.340 3.502\n",
"TV 0.0458 0.001 32.909 0.000 0.043 0.048\n",
"Radio 0.1880 0.008 23.382 0.000 0.172 0.204\n",
"==============================================================================\n",
"Omnibus: 60.022 Durbin-Watson: 2.081\n",
"Prob(Omnibus): 0.000 Jarque-Bera (JB): 148.679\n",
"Skew: -1.323 Prob(JB): 5.19e-33\n",
"Kurtosis: 6.292 Cond. No. 425.\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# print a summary of the fitted model that includes only TV and Radio in the model\n",
"lm = smf.ols(formula='Sales ~ TV + Radio', data=data).fit()\n",
"lm.rsquared\n",
"lm.summary()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Linear Regression in scikit-learn\n",
"\n",
"Let's redo some of the Statsmodels code above in scikit-learn:"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2.93888936946\n",
"[ 0.04576465 0.18853002 -0.00103749]\n"
]
}
],
"source": [
"# create X and y\n",
"feature_cols = ['TV', 'Radio', 'Newspaper']\n",
"X = data[feature_cols]\n",
"y = data.Sales\n",
"\n",
"# follow the usual sklearn pattern: import, instantiate, fit\n",
"from sklearn.linear_model import LinearRegression\n",
"lm = LinearRegression()\n",
"lm.fit(X, y)\n",
"\n",
"# print intercept and coefficients\n",
"print lm.intercept_\n",
"print lm.coef_"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"[('TV', 0.045764645455397615),\n",
" ('Radio', 0.18853001691820442),\n",
" ('Newspaper', -0.001037493042476266)]"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# pair the feature names with the coefficients\n",
"zip(feature_cols, lm.coef_)"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"12.202667011892375"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# predict for a new observation\n",
"lm.predict([100, 25, 25])"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0.89721063817895219"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# calculate the R-squared\n",
"lm.score(X, y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that **p-values** and **confidence intervals** are not (easily) accessible through scikit-learn."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.10"
}
},
"nbformat": 4,
"nbformat_minor": 0
}