{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import cvxpy as cvx\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "\n", "from sklearn import linear_model as lm\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exploratory data analysis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read in the *Engel food expenditure* data from the R package [`quantreg`](https://cran.r-project.org/package=quantreg)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "engel = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/quantreg/engel.csv', index_col=0)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
incomefoodexp
1420.157651255.839425
2541.411707310.958667
3901.157457485.680014
4639.080229402.997356
5750.875606495.560775
\n", "
" ], "text/plain": [ " income foodexp\n", "1 420.157651 255.839425\n", "2 541.411707 310.958667\n", "3 901.157457 485.680014\n", "4 639.080229 402.997356\n", "5 750.875606 495.560775" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "engel.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Scatter plot of income vs food expenditure" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XuU3HV9//Hna2/ZzYUQYcmJQLJgoQZtTcJysVoLEi7SSiinUFIrEba/BCuK0dMq1B/2p8Kh/aGpQAuJv6DiwVWoCGgjMUHQ1grJ5kIAF0uAjURywggYI5tkb+/fH/Odzezu7O5MsrOzs/t6nDNnZz77mZnPfCHz3s/t/VFEYGZmVoiKUjfAzMzKj4OHmZkVzMHDzMwK5uBhZmYFc/AwM7OCOXiYmVnBHDzMzKxgDh5mZlYwBw8zMytYVakbUCxHH310NDQ0lLoZZmZlY9OmTb+OiPp86o7b4NHQ0EBLS0upm2FmVjYk7ci3roetzMysYA4eZmZWMAcPMzMrmIOHmZkVzMHDzMwK5uBhNgalUu1s3LiLVKq91E0xy8nBw2yMaW5uZc6cVZx77n3MmbOK5ubWUjfJbAAHD7MxJJVqp6lpLfv2dbFnTwf79nXR1LTWPRAbc4oWPCQdL+lRSa2SnpF0bVL+JknrJD2X/JyRlEvSrZK2S9omaUHWay1J6j8naUmx2mxWam1te6ip6fvPsrq6gra2PSVqkVluxex5dAGfjIi5wJnARySdAnwaeCQiTgIeSR4DvA84KbktBe6AdLABPgucAZwOfDYTcMzGm4aG6XR09PQp6+zsoaFheolaZJZb0YJHROyKiM3J/b1AK3AssAj4elLt68DFyf1FwN2R9jhwpKRZwPnAuoh4LSJeB9YBFxSr3WalVF8/mdWrz6euroojjqihrq6K1avPp75+cqmbZtbHqOS2ktQAzAeeAGZGxC5IBxhJxyTVjgVeynrazqRssHKzcWnx4rksXDiHtrY9NDRMd+CwManowUPSVOA7wMcj4reSBq2aoyyGKM/1XktJD3kxe/bswhtrNkbU10920LAxrairrSRVkw4c90TE/Unx7mQ4iuTnK0n5TuD4rKcfB7w8RPkAEbEqIhojorG+Pq+swmZmdgiKudpKwGqgNSK+lPWrh4DMiqklwINZ5Vckq67OBPYkw1trgfMkzUgmys9LyszMrESKOWz1LuCDwFOStiZl1wM3A/dKagJ+CVya/G4NcCGwHWgHrgSIiNckfR7YmNT7XES8VsR2m5nZMBSRc/qg7DU2NoYPgzIzy5+kTRHRmE9d7zA3M7OCOXiYmVnBHDzMzKxgDh5mZlYwBw8zMyuYg4eZmRXMwcPMzArm4GFmZgVz8DAzs4I5eJiZWcEcPMzMrGAOHmZm40Qq1c7GjbtIpdqL/l4OHmZm40Bzcytz5qzi3HPvY86cVTQ3txb1/Rw8zMzKXCrVTlPTWvbt62LPng727euiqWltUXsgDh5mZmWurW0PNTV9v86rqytoa9tTtPd08DAzK3MNDdPp6OjpU9bZ2UNDw/SivaeDh5lZmauvn8zq1edTV1fFEUfUUFdXxerV51NfP7lo71nMY2jNzGyULF48l4UL59DWtoeGhulFDRzg4GFmNm7U108uetDIKNqwlaS7JL0i6emssm9L2prc2iRtTcobJO3L+t2dWc85VdJTkrZLulWSitVmMzPLTzF7Hl8DbgfuzhRExF9m7kv6IpC9FOD5iJiX43XuAJYCjwNrgAuAHxShvWZmlqei9Twi4ifAa7l+l/QeLgOah3oNSbOAIyLiZxERpAPRxSPdVjMzK0ypVlv9MbA7Ip7LKjtB0hZJP5b0x0nZscDOrDo7k7KcJC2V1CKpJZVKjXyrzcwMKF3wWEzfXscuYHZEzAc+AXxT0hFArvmNGOxFI2JVRDRGRGN9ff2INtjMzA4a9dVWkqqAS4BTM2URcQA4kNzfJOl54GTSPY3jsp5+HPDy6LXWzMxyKUXPYyHwbET0DkdJqpdUmdw/ETgJeCEidgF7JZ2ZzJNcATxYgjabmVmWYi7VbQZ+Bvy+pJ2SmpJfXc7AifL3ANskPQn8O3B1RGQm2z8M/D9gO/A8XmllZlZySi9iGn8aGxujpaWl1M0wMysbkjZFRGM+dZ3byszMCubgYWZmBXPwMDOzgjl4mJlZwRw8zMysYA4eZmZWMAcPMzMrmIOHmZkVzMHDzMwK5uBhZmYFc/AwM7OCOXiYJVKpdjZu3EUq1V7qppiNeQ4eZkBzcytz5qzi3HPvY86cVTQ3t5a6SWZjmoOHTXipVDtNTWvZt6+LPXs62Levi6amte6BmA3BwcMmvLa2PdTU9P2nUF1dQVvbnhK1yGzsc/CwCa+hYTodHT19yjo7e2homF6iFpmNfQ4eNuHV109m9erzqaur4ogjaqirq2L16vOpr59c6qaZjVlVpW6A2ViwePFcFi6cQ1vbHhoapjtwmA2jmGeY3yXpFUlPZ5X9o6RfSdqa3C7M+t11krZL+oWk87PKL0jKtkv6dLHaa+VpJJfX1tdP5rTTZjlwmOWhmMNWXwMuyFG+IiLmJbc1AJJOAS4H3pY8598kVUqqBP4VeB9wCrA4qWtW1OW13vNhNrSiBY+I+AnwWp7VFwHfiogDEfEisB04Pbltj4gXIqID+FZS1ya4Yi6v9Z4Ps+GVYsL8GknbkmGtGUnZscBLWXV2JmWDldsEV6zltd7zYZaf0Q4edwBvAeYBu4AvJuXKUTeGKM9J0lJJLZJaUqnU4bbVxrBiLa/1ng+z/Ixq8IiI3RHRHRE9wFdID0tBukdxfFbV44CXhygf7PVXRURjRDTW19ePbONtTCnW8lrv+TDLz6gu1ZU0KyJ2JQ//HMisxHoI+KakLwFvBk4CNpDueZwk6QTgV6Qn1f9qNNtsY1cxltdmglJT01qqqyvo7Ozxng+zHIoWPCQ1A2cBR0vaCXwWOEvSPNJDT23AMoCIeEbSvcDPgS7gIxHRnbzONcBaoBK4KyKeKVabrfzU108e8S927/kwG54iBp1CKGuNjY3R0tJS6maYmZUNSZsiojGfuk5PYlYA7/8wS3PwMMuT93+YHeTgYZYH7/8w68vBwywP3v9h1peDh1kevP/DrC8HD7M8+MwPs758nodZnrz/w+wgBw+zAhRjU6JZOfKwlZmZFczBw8zMCubgYWZmBXPwMDOzgjl42ITkHFVmh8fBwyYc56gyO3wOHlZWDrfHUIwcVe7F2ETk4GFlYyR6DCOdo8q9GJuoHDysLIxUj2Ekc1Q5065NZA4eVhZGqscwkjmqnGnXJrJh05NIqgX+Fng36bPH/wu4IyL2F7ltZr1GsscwUjmqnGnXJrJ8eh53A28DbgNuB+YC3xjuSZLukvSKpKezyv6vpGclbZP0XUlHJuUNkvZJ2prc7sx6zqmSnpK0XdKtklToh7TyN9JZbevrJ3PaabMOK0+VM+3aRKaIGLqC9GREvGO4shzPew/wO+DuiHh7UnYe8KOI6JL0TwAR8SlJDcD3M/X6vc4G4FrgcWANcGtE/GC4D9bY2BgtLS3DVbMyk0q19+kx9H88FtpkVq4kbYqIxnzq5pNVd4ukMyPi8eTFzwB+OtyTIuInSVDILvth1sPHgb8Y6jUkzQKOiIifJY/vBi4Ghg0eNj5lZ7Vtbm6lqWktNTUVdHT0sHr1+SxePLekbTKbKPIZtjoD+G9JbZLagJ8Bf5IMJW07jPe+ir5B4ARJWyT9WNIfJ2XHAjuz6uxMymyC80ons9LKp+dxwUi/qaR/ALqAe5KiXcDsiHhV0qnAA5LeBuSa3xh0nE3SUmApwOzZs0e20TamZFY67dt3sCyz0sm9ALPiy6fncVJE7Mi+AWdl3S+IpCXAnwEfiGTCJSIORMSryf1NwPPAyaR7GsdlPf044OXBXjsiVkVEY0Q01tfXF9o0KyNe6WRWWvkEjxsk3SFpiqSZkr4HvP9Q3kzSBcCngIsioj2rvF5SZXL/ROAk4IWI2AXslXRmssrqCuDBQ3lvG1+80smstPIZtvoT4JPA1uTxDRHRPNyTJDUDZwFHS9oJfBa4DpgErEtW3D4eEVcD7wE+J6kL6AaujojXkpf6MPA1oI70HIknyw3wmeJmpZRP8JhBetL8edLDRnMkKTPkNJiIWJyjePUgdb8DfGeQ37UAA5bwmoFXOpmVSj7DVo8DP4iIC4DTgDeTx1JdMzMbv/LpeSyMiF8CRMQ+4GPJBkAzb5Azm6Dy6Xm8JOmvJd0AIGk24LxW5nTkZhNYPsHj34B3Apk5jL3AvxatRVYWRmqTng9SMitP+QxbnRERCyRtAYiI1yXVFLldNsYNtUkv8/vhhrLGSnoRMytcPj2PzmQPRkB6TwbQM/RTbLwbbJPe5s278xrKcnoRs/KWT/C4FfgucIykG0mf53FTUVtlY16uTXorVpzF8uWP5RUQfJCSWXkbdtgqIu6RtAk4h3SuqYsjwjOjNmCTXiH5ppxexKy8DRo8JL0p6+ErQHP277J2gNsE1n+TXr4BIdNzaWpaS3V1BZ2dPU4vYlZGhup5bCI9zyFgNvB6cv9I4JfACUVvnZWVQgOC04uYla9Bg0dEnACQHAn7UESsSR6/D1g4Os2zclNoQHB6EbPylM9S3dOS5IUARMQPJH2+iG2yMtN/l7kDgtn4l89qq19L+oykBklzkoOcXi12w6w8jKVd5t5waDZ68gkei4F60st1HwCO4eBuc5vAxtJejbEUxMwmgnyW6r4GXCvpCKAnIn5X/GbZWJU9RDVWjoLNDmKZtjQ1rWXhwjkePjMrkmF7HpL+IElN8hTwjKRNkny+xgTU/6/7zZt3j4m9Gt5waDb68hm2Wgl8IiLmRMQc0qcKripus2ysyTVEtXz5Y6xYcXbRj4Idbi7DGw7NRl8+q62mRMSjmQcR8ZikKUVsk5XQYOdzDDZEtWDBMezYsbRoezXySZ7oDYdmo0/DnCaLpO8Cm4FvJEV/DTRGxMVFbtthaWxsjJaWllI3o6wM9UWdSrUzZ84q9u3r6q1fV1fFjh1Li/YlXeh7+mAqs8MjaVNENOZTN59hq6tIr7a6n/SKq3rgyjwbcpekVyQ9nVX2JknrJD2X/JyRlEvSrZK2S9omaUHWc5Yk9Z+TtCSf97bCDLdyKlcixGL/dV/oXEZ9/WROO22WA4fZKMhntdXrpI+ePZTVVl8Dbgfuzir7NPBIRNws6dPJ408B7wNOSm5nAHcAZyQ5tj4LNJJOl7JJ0kNJu2yE5LNyqtjpRPr3HDyXYTZ2FXW1VUT8BOifQHER8PXk/teBi7PK7460x4EjJc0CzgfWRcRrScBYB1yQz/tb/vL9oi7WX/eZlVznnHMfxx+/kpUrnyxJb8fM8pPPhHlmtdWjAJLOIr3a6o8O8T1nRsQugIjYJemYpPxY4KWsejuTssHKB5C0FFgKMHv27ENs3sR0qJPOIzHPkD1klnH11euAYNmyeU6eaDYGjaXVVspRFkOUDyyMWEWyjLixsXHolQA2QKHDUiN1jGxb2x6qqgZ2gq+99lEuueTkPrmyPCluNjbkM2H+gqT/neS2apD0GeDFw3jP3clwFMnPV5LyncDxWfWOA14eotyKIDPX0Na2Z8g0IyOZmiQ9ZNY9oLz/5LhTkJiNHYMGD0mZpbn/Sd/VVkeT52qrQTwEZFZMLQEezCq/Ill1dSawJxneWgucJ2lGsjLrvKTMiiDfL+iR3NVdXz+ZL3/5vQPKu7ujd85lLOXRMrOhh61OlTSH9Bf82aSHjzJDQbmGkgaQ1AycBRwtaSfpVVM3A/dKaiJ9qNSlSfU1wIXAdqCdJEBFxGtJCviNSb3P+RTDkdF/CKiQHFEjvRJq2bJ3AMG11z5KdXUF3d3RZ85lrOTRMrO0oYLHncDDwIlA9m67TBA5cbgXj4jBsu+ek6NuAB8Z5HXuAu4a7v1seJmAsXnzbpYvf6zPfMXv/d6RA76gq6oqWLPmBS688MQ+X9LF2NW9bNk8Lrnk5JxzGl62aza25LPD/I6I+PAotWfEeIf5QJkJ7qqqCvbu7ejzu7q6KjZt+iCnnvqNPqueAKZNq6arK3JOiI/mBHam/dnB6lAm6M0st0J2mA8bPMqVg0dfuVJ9ZDviiBrWr7+U7dt/M2SAKWY6kv7tzRWUvNrKrHhGOj2JjQO5JrizZYaAFi+ey44dS7nttvcybVp1nzqjleZ8qEl7pyAxGxscPCaIXHMG2VasOLv3C7m+fjIXXngiXV19e6UdHd1Fn2Pwqiqz8uDgMUFkJrgnTaoc8LupU6tZsOCYnPWrqw8urOvpCdav3zHg+SN5drgPdjIrDw4eE8jixXPZsuWKAQEkez9FtoUL51BVdbBuR0fPgF5A/yGmlSufPKxA4lVVZuXBwWOCmTv3KL761QvySjY4XC8g1xDT1Vev45xz7mX27JV84Qs/KziIlHsyxJHshZmNZV5tNY4NtTIpn1VLwx3GtHHjLs499z727OnI+fxM/UNZUluOq6pGKteXWal4tZUNm2Ykn1VLw/UChpuEBw55wrvcVlV5ot8mGgePcWgkv8gyS3fXr7+UHTuWDvhL+vrrz6Curopp02oGfY2JMOHtiX6baBw8xqFcX2QVFbBly+5BnzPUWH2uXkCmZ3PLLRuJCP7+70/jzjsXUls7cDVXPhPe5T5X4Il+m2gcPMpMPl+yDQ3TOXCgb4rzN97o4uKLH8yZJbfQVOf9ezb793dz001PcMklJ/PLXy7j859/V0ET3uMh1Xq5T/SbFcoT5mUknwnZVKqdL36xhVtu2UD3wCMymDSpki1bruDoo+toa9vD1Kk1A/JZDZeGZOPGXZxzzn190pdk0pucdtqs3nbkM+E93KR8uSnHiX6zjEImzPM5SdDGgHzSpTc3t3LFFf9BV+70VQAcONDN29/+VSoqYMqUGvbv76aiX/9zuFTnmzfvHpD3qv8QTfbpf0MZb6nW8/3cZuXOw1ZlIp89F1dd9fCQgSOjpwe6umDPng4OHOhm376+XZQDB7qZOjX3BHgq1c7y5Y8NKM9Ob1IIzxWYlScHjzIx3JdsW9seKivzOqNrgLq6KiZNquyd7K6oEKee+o2ccw+5gliu9Cb58lyBWXly8CgTub5kV6w4q/es8YaG6XR3556/yieoPPLIZWSmv9JDY7mX906dmh7qyjZYepN8Dbcc2MzGHs95lJHFi+eycOGcQU8CvOuuC1iyZA2dnQeDSHV1BZdeejLf/Oazg75uU9PbqampoLa2ss8qrf5zD5kJ+4qKdDCqra1E0oj0FDxXYFZeRn21laTfB76dVXQicANwJPC/gFRSfn1ErEmecx3QBHQDH4uItcO9T7mttipklc5QK5QgvZ/jN785wJFHTmLbthR/93c/GfL1amsr2bz5iiFXXeV6z8zKrblzjzqUj2xmY8yYXm0VEb8A5gFIqgR+BXwXuBJYERG3ZNeXdApwOfA24M3AekknR0SOhajlqdCcSLlWKGU2AZ533gmcd94JQDrIXHTRA3m14aWXfjvkmeS53nPSpEp+97vB81qZ2fhV6jmPc4DnI2LgIREHLQK+FREHIuJFYDtw+qi0bhQcSiqRXJPnb7zRxaJFD/SZ5E5/4Q/c8d3f/v3dXHzxgwCDzj14VZSZZSt18LgcaM56fI2kbZLukjQjKTsWeCmrzs6kbFzIZwlu/x3l9fWTWbHi7AHP27+/mw996Ae9dRsaptPVNXTiwoxM0AJyJiTMd1VUuacZMbP8lCx4SKoBLgLuS4ruAN5CekhrF/DFTNUcT885USNpqaQWSS2pVCpXlTFnqL/oB0vb0dzcyvLlj/ZOXGfr6OjpzWGV+cLPlW8ql+ES+Q23Kmo8pBkxs/yUsufxPmBzROwGiIjdEdEdET3AVzg4NLUTOD7reccBL+d6wYhYFRGNEdFYX19fxKaPnMH+ogdyDme1tr7aW95/yWzGb35zoPf+4sVzefDBi5kypXpAvaqqvv/58xmGGixVulOSm00spVyqu5isIStJsyJiV/Lwz4Gnk/sPAd+U9CXSE+YnARtGs6HFlr0EN7PaauPGXTnTdqxf30ZPz9Ar5JYs+QHd3dHbM5g/f+aA51RXV5C90q6mpmLAMFQhK8DGW5oRMxtaSXoekiYD5wL3ZxX/s6SnJG0DzgaWA0TEM8C9wM+Bh4GPjKeVVhn9/6Lv6OjpsywWoL29k09+8tEBGXP727+/u/ev/kwAWLHi7D69Gwi6ug4Gj4oKsXDhnN7HhQxBpVLtvP76fk+om00gJel5REQ7cFS/sg8OUf9G4MZit2us+OhH13P77Vt7H1dWQnc3fb7sh1NdXcHKlU9y001P9C4BXrHiLBYsmMnrr+/nssu+1+f42Jqayt5eQj5JGDOylxl3dXUnmw2rBiz1NbPxxTvMx5if/vRXfQIHkDO1+nA6O3u48cbH2b+/uzcALF/+GDt2LM0rT1Y+Q1C5gkxtbQX33fd+5s+f6cBhNo6VeqmuZWlubuW97/328BX7qamp4Jpr5vUZllqy5BSqq/uussoOAMOdTd7e3tnnufv2dQ4Ygsq1zLimppIZM2odOMzGOfc8xojW1le58sqHB/QIhlJRAc3Nf8bZZ8+mvn4yN9zwR6xc+SQ33vg499zTyt69fQNAdu8i1yQ9pHsT6aW+IntFtDRwWbA3DppNXO55jAHNza3Mn3/3sBPh/dXWVnLkkZP6/JV/001PsH9/d5/AMW1a7k19/SfpM5Pkl1zyIJ2dfYNCbW3VgD0gTqduNnG551FimXmDwQJHZrI8l/b2dFqRTC6sXHMVU6ZUs3z5qVx++VuHTGCYPX+Ry2A9isF6MGY2vrnnUWK55g2yDTdZnr0ZL3fOq05WrGgZ9HCn4doxZUr1sD2KwTYOmtn45eBRYg0N0w87M+2+fV2sXPlkn2GkadMOHiO7d2/nsDu+cwWe2tpK7r//Ih/QZGYDOHiU2K9/ve+QluL2d+ONj5NKtffmn7rttvcyZUrfUcmhclflmr+4664LOO+8E9yjMLMBPOdRYuvXD5WNPn/Zm/zq6yezf38Xb7zRd/5iuJVQnr8ws3w5eJRYOlXI4evo6Ob11/f3DkstX/7YgDorVpw1bEDwcbBmlg8Hj1EwVILB55//zWG99tSp1Rw40EVPT3DZZd+jo6OHa69dQEW/Aclp02pYsGDmYb2XmVmG5zyKbKgEg6lUOytWHNo567W1ldx557l85zsXUVVVSUdHT28q9Jtv3jBgyKqry5v3zGzkOHgU0XBnXKxcuZUDB/LfUQ7p7LdXX/2HbN58BcuWvYMZM2qHXOqb0dT0dg9HmdmIcfAooi1bdg847S+z4imVaudzn3u84NesrhbNzc/27tvItcQ2l9Wrn/bBTGY2Yhw8iqS5uZVFix7gjTdy55fasmX3gBQgw6msFAcO9PTpxQA593b0N9wRs2ZmhXDwKILMcFX/Y2Jrayu5/vozDuk1q6sr6O7ue55HVVU6IGT2djzyyKXceefCnCu4nLDQzEaSg0cR5Er1MWlSJT09wS23bGTOnFVs2/brnM+trBSVlQMz2Obqpezd28HmzbuBgylCli2bx44dS/n8599FbW0lU6dWM2lSZV7LdM3M8qXsc6zHk8bGxmhpObSVTIcrlWpnzpxVgyYZhPTchaQ+8xWVlenU54WcGFhXV8WOHUtzBoaVK7dy7bWPUlNTSVdXT28CRTOzXCRtiojGfOq651EE/VN91NRUUNn3XCY6O2PARHehR81C+kyPNWteGDAZnkq1s3z5Yxw40M3evQNXepmZHY6SBQ9JbZKekrRVUktS9iZJ6yQ9l/yckZRL0q2StkvaJmlBqdqdr8WL57Jp0wf5q796Kz09kVf+qtrayuEr9fPGG1189KM/GrCHJNfQmSfNzWyklLrncXZEzMvqJn0aeCQiTgIeSR4DvA84KbktBe4Y9ZYWqLm5lQUL7ubOO7fl3ZvoPyGer1w9C5/yZ2bFVOrg0d8i4OvJ/a8DF2eV3x1pjwNHSppVigYOJ5Vq54c/fDHnaqvhRBS2dHeonoVP+TOzYiplbqsAfigpgJURsQqYGRG7ACJil6RjkrrHAi9lPXdnUrYr+wUlLSXdM2H27NlFbn5adt6q9et30NS0looKhpwsH8ykSVVUVwf79g0edKqrK6ioEF/4wru54YafAgcDTv+ehbPkmlmxlDJ4vCsiXk4CxDpJzw5Rd+Da1XTw6VuQDkCrIL3aamSaObjm5lauuuphKirSE909PT10FR4zenV3B1Kuj5oOGvffv4iZMyf3BoJjj51KU9Naqqsr6OzsydmzcJZcMyuGkgWPiHg5+fmKpO8CpwO7Jc1Keh2zgFeS6juB47Oefhzw8qg2uJ9Uqp0lS9bQ2TkyMaqmpoK77rqA3/62g6uvXjfg95MmVTJzZnovR4Z7FmZWKiWZ85A0RdK0zH3gPOBp4CFgSVJtCfBgcv8h4Ipk1dWZwJ7M8FappNOLjFznJpMDa8GCY5gypXrA7web7Pb54WZWCqWaMJ8J/JekJ4ENwH9ExMPAzcC5kp4Dzk0eA6wBXgC2A18B/nb0mzxyLrroRCZN6rssd//+bpqa1jJ1ag09PQOD0pe/fLYDhJmNGSUJHhHxQkS8I7m9LSJuTMpfjYhzIuKk5OdrSXlExEci4i0R8QcRUZqt41m2bUsd8nPXrHmR22/Pfcb4Sy/9luuvP4Pa2kqmTath0qT0uR3Lls073CabmY0Ypyc5BK2trzJv3tfzSoU+mMpKUVFBn6Gv6mpRVVVJTU0FHR09XH/9GSxb9g73OMxsVDg9SRE1N7cyf/7dhxU4IL2yKiK9q/yII2qora1EUp+Do2666YkRarWZ2chy8ChAKtXOVVc9zIEDhW3+G0xXV3DZZb/Pffe9nwcfvHhAKnWnEzGzscrBowArV24teNf4cO6+++csWvQAjz760oCg5HQiZjZWlXKTYNlIpdrZsmU3N9546MNIFRXQM8hI1/793dx88waqq0V1tairqx5005+Z2Vjg4DGM5ubW3pQjh9PrGCxwZOvsDOrqqrjvvvczf/5MBw4zG7McPIaQOU72UPJUHarq6gpmzKh14DCzMc0DeNzzAAAJgklEQVRzHkPIdSbGSKmpqcj52p7nMLNy4OAxhIaG6SO2sqq/ysoKtm5d0nvWuNOmm1k5cfAYQn39ZD7+8VMP+fmLF781Z3kmSMydexSf+cw7+eUvl7F+/aXs2LHUZ4ybWVlw8BjCypVb+eIXNx7Sc5csOYUPfehtTJvWN8nhlClVPPDAoj5BwskNzazceMJ8ECtXbuXqq9cf8vPvv3879977P3R19R326umB+fNnHm7zzMxKysEjh1SqnY9+9EeH9Rp793YA6Ynx2toKamoqvXfDzMYNB48c0md1FJ67qrJSVFdX9NkPUlub3rcxY0atD2wys3HDcx45/OhHLw1fqZ+amgp+/OPLBxwj29nZw/z5Mz2nYWbjioNHP6lUO7feurng51VUiJqaClavPp+6uiovvTWzcc3DVv1s2bI750l+w9m/v5upU2t8rriZTQgOHlmam1u58sofcOBA4fMddXVV/O536Uny+vrJDhpmNq6N+rCVpOMlPSqpVdIzkq5Nyv9R0q8kbU1uF2Y95zpJ2yX9QtL5xWhXJo/VoQSODKcVMbOJohQ9jy7gkxGxWdI0YJOkdcnvVkTELdmVJZ0CXA68DXgzsF7SyRExonlDMnms9u3Lr35lZTrFSG1tlZfgmtmEM+rBIyJ2AbuS+3sltQLHDvGURcC3IuIA8KKk7cDpwM9Gsl0NDdPzOlr2Yx+bz5/+6Ym9G/08t2FmE1FJV1tJagDmA5lTlq6RtE3SXZJmJGXHAtlrZ3cydLA5JPX1k3n3u988bL2VK7f1nrXhtCJmNlGVLHhImgp8B/h4RPwWuAN4CzCPdM/ki5mqOZ6eczmUpKWSWiS1pFKpgtrT2voq69b9cth6NTWVPlfczCa8kgQPSdWkA8c9EXE/QETsjojuiOgBvkJ6aArSPY3js55+HPByrteNiFUR0RgRjfX19QW1acOGXXnV6+ryeRtmZqVYbSVgNdAaEV/KKp+VVe3PgaeT+w8Bl0uaJOkE4CRgw0i36/TTZ+Us/5u/eTuTJlUybZo3/ZmZZZRitdW7gA8CT0nampRdDyyWNI/0kFQbsAwgIp6RdC/wc9IrtT4y0iutAObOPYprrpnH7bdv7S275pp53HbbQm666T2eGDczy6KIwndTl4PGxsZoaWkp+Hmtra+yYcMuTj99FnPnHlWElpmZjU2SNkVEYz51vcO8n7lzj3LQMDMbhhMjmplZwRw8zMysYA4eZmZWMAcPMzMrmIOHmZkVbNwu1ZWUAnaUuh1FdDTw61I3YozwtUjzdTjI1yKt0OswJyLySs8xboPHeCepJd/12OOdr0War8NBvhZpxbwOHrYyM7OCOXiYmVnBHDzK16pSN2AM8bVI83U4yNcirWjXwXMeZmZWMPc8zMysYA4eY0hy/O4rkp7OKnuTpHWSnkt+zkjKJelWSduTo3sXZD1nSVL/OUlLSvFZDoek4yU9KqlV0jOSrk3KJ9S1kFQraYOkJ5Pr8H+S8hMkPZF8pm9LqknKJyWPtye/b8h6reuS8l9IOr80n+jwSKqUtEXS95PHE/U6tEl6StJWSS1J2ej/24gI38bIDXgPsAB4Oqvsn4FPJ/c/DfxTcv9C4Aekj+k9E3giKX8T8ELyc0Zyf0apP1uB12EWsCC5Pw34H+CUiXYtks8zNblfDTyRfL57gcuT8juBDyf3/xa4M7l/OfDt5P4pwJPAJOAE4HmgstSf7xCuxyeAbwLfTx5P1OvQBhzdr2zU/22U/EL4NuB/jIZ+weMXwKzk/izgF8n9lcDi/vWAxcDKrPI+9crxBjwInDuRrwUwGdgMnEF601dVUv5OYG1yfy3wzuR+VVJPwHXAdVmv1VuvXG6kj59+BHgv8P3kc02465C0O1fwGPV/Gx62GvtmRsQugOTnMUn5scBLWfV2JmWDlZelZMhhPum/uifctUiGarYCrwDrSP+1/JuI6EqqZH+m3s+b/H4PcBTj4DoA/wL8PdCTPD6KiXkdIH3a6g8lbZK0NCkb9X8bPgyqfClHWQxRXnYkTQW+A3w8In4r5fpo6ao5ysbFtYj0kcvzJB0JfBeYm6ta8nNcXgdJfwa8EhGbJJ2VKc5RdVxfhyzvioiXJR0DrJP07BB1i3Yt3PMY+3ZLmgWQ/HwlKd8JHJ9V7zjg5SHKy4qkatKB456IuD8pnpDXAiAifgM8Rnrc+khJmT/8sj9T7+dNfj8deI3yvw7vAi6S1AZ8i/TQ1b8w8a4DABHxcvLzFdJ/UJxOCf5tOHiMfQ8BmZUQS0iP/2fKr0hWU5wJ7Em6q2uB8yTNSFZcnJeUlQ2luxirgdaI+FLWrybUtZBUn/Q4kFQHLARagUeBv0iq9b8OmevzF8CPIj2g/RBwebIK6QTgJGDD6HyKwxcR10XEcRHRQHoC/EcR8QEm2HUAkDRF0rTMfdL/Tz9NKf5tlHryx7c+k17NwC6gk/RfBk2kx2ofAZ5Lfr4pqSvgX0mPgT8FNGa9zlXA9uR2Zak/1yFch3eT7kJvA7Ymtwsn2rUA/hDYklyHp4EbkvITSX/pbQfuAyYl5bXJ4+3J70/Meq1/SK7PL4D3lfqzHcY1OYuDq60m3HVIPvOTye0Z4B+S8lH/t+Ed5mZmVjAPW5mZWcEcPMzMrGAOHmZmVjAHDzMzK5iDh5mZFczBw6wAkv671G0wGwu8VNfMzArmnodZAST9Lvl5lqTHJP27pGcl3ZPsjEfSaZL+OzmHY4OkaUqfzfHV5ByGLZLOTup+SNIDkr4n6UVJ10j6RFLncUlvSuq9RdLDSTK8/5T01tJdBTMnRjQ7HPOBt5HOCfRT4F2SNgDfBv4yIjZKOgLYB1wLEBF/kHzx/1DSycnrvD15rVrSu30/FRHzJa0AriCdx2kVcHVEPCfpDODfSOd4MisJBw+zQ7chInYCJGnTG0in/94VERsBIuK3ye/fDdyWlD0raQeQCR6PRsReYK+kPcD3kvKngD9Msgv/EXBfVmbhSUX+bGZDcvAwO3QHsu53k/73JHKnth40n3y/1+nJetyTvGYF6bMr5h16U81Gluc8zEbWs8CbJZ0GkMx3VAE/AT6QlJ0MzCadnG9YSe/lRUmXJs+XpHcUo/Fm+XLwMBtBEdEB/CVwm6QnSZ/+V0t6jqJS0lOk50Q+FBEHBn+lAT4ANCWv+QywaGRbblYYL9U1M7OCuedhZmYFc/AwM7OCOXiYmVnBHDzMzKxgDh5mZlYwBw8zMyuYg4eZmRXMwcPMzAr2/wEaM0XhK2FRogAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "engel.plot.scatter(x='income', y='foodexp', color='darkblue')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Modelling" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define model matrix and outcome." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "X = np.hstack((\n", " np.ones(engel.shape[0])[:, np.newaxis],\n", " engel['income'][:, np.newaxis]\n", "))\n", "y = engel['foodexp'].values" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 1. , 420.15765084],\n", " [ 1. , 541.41170672],\n", " [ 1. , 901.15745665],\n", " [ 1. , 639.08022869],\n", " [ 1. , 750.87560583]])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X[:5, :]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([255.83942459, 310.95866706, 485.68001417, 402.99735554,\n", " 495.56077493])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y[:5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define a variable `betas` for the regression coefficients and a parameter `tau` representing the quantile." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "betas = cvx.Variable(X.shape[1])\n", "tau = cvx.Parameter()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define variables `u` and `v` for the positive and negative deviations." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "u = cvx.Variable(X.shape[0], nonneg=True)\n", "v = cvx.Variable(X.shape[0], nonneg=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the objective function and constraints." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "objective = cvx.sum(tau * u) + cvx.sum((1-tau) * v)\n", "constraints = [cvx.matmul(X, betas) + u - v == y]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set a value for $\\tau$." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "tau.value = 0.5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define and solve the quantile regression problem." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "8781.383825942727" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "problem = cvx.Problem(cvx.Minimize(objective), constraints)\n", "problem.solve()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Extract the regression coefficients." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([77.09390334, 0.56467913])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "betas.value" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define a new range of values to predict over." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "n_pred = 100\n", "X_pred = np.hstack((\n", " np.ones(n_pred)[:, np.newaxis],\n", " np.linspace(min(engel['income']), max(engel['income']), n_pred)[:, np.newaxis]\n", "))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 1. , 377.05836885],\n", " [ 1. , 423.3286179 ],\n", " [ 1. , 469.59886694],\n", " [ 1. , 515.86911599],\n", " [ 1. , 562.13936504]])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X_pred[:5, :]" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "y_pred = X_pred @ betas.value" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([290.01089332, 316.13873709, 342.26658085, 368.39442462,\n", " 394.52226839])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y_pred[:5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the original data and quantile regression line." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmczfX3wPHXmc0MY28sESNRkkJjK/1SCflWpI02WcYSCW1KkYqvSlHKFipiQohUJkJE9qUwyXwZ2WLs2+z3/fvjfu64M3Nn5s6YO3eW83w87mPufd/3/dxzP8w983mvYoxBKaWUcpePtwNQSilVuGjiUEoplSOaOJRSSuWIJg6llFI5oolDKaVUjmjiUEoplSMeSxwiEigiG0Vkh4jsEpERVnktEdkgIntFZI6IBFjlJazH0dbzoU7Hes0q3yMibT0Vs1JKqex58oojAbjbGHML0BBoJyLNgfeAscaYOsBpoIdVvwdw2hhzHTDWqoeI3Ah0BuoD7YAJIuLrwbiVUkplwWOJw9hdsB76WzcD3A18a5V/BXS07newHmM9f4+IiFX+jTEmwRizH4gGmnoqbqWUUlnz8+TBrSuDLcB1wGfA/4Azxphkq8ohoJp1vxpwEMAYkywiZ4GKVvl6p8M6v8alq666yoSGhubRp1BKqeJhy5YtJ4wxIdnV82jiMMakAA1FpBywEKjnqpr1UzJ5LrPyNESkF9ALoEaNGmzevDlXMSulVHElIgfcqZcvo6qMMWeAVUBzoJyIOBJWdeCIdf8QcA2A9XxZ4JRzuYvXOL/HFGNMmDEmLCQk24SplFIqlzw5qirEutJARIKA1kAUsBJ4xKrWFVhk3V9sPcZ6foWxr8C4GOhsjbqqBdQBNnoqbqWUUlnzZFNVVeArq5/DB5hrjFkiIruBb0TkXWAbMM2qPw2YKSLR2K80OgMYY3aJyFxgN5AM9LOawJRSSnmBFMVl1cPCwoz2cSilVM6IyBZjTFh29XTmuFJKqRzRxKGUUipHNHEopZTKEU0cSilVRHy37TumrZmWfcUrpIlDKaUKuSNnjvDwxId5aMJDTFs7DZvN5tH38+jMcaWUUp5js9mY+ttUXvn2FRKSE/hvp//y4r0v4uPj2WsCTRxKKVUI7fl3D+Ezwlmzdw2trm/FlKenUKdynXx5b00cSilViCQmJ/L+0vd554d3KBVQimldp9Ht9m7YFxPPH5o4lFKqkFj/v/X0nNGTXUd28VjYY3zS5RMql6mc73Fo4lBKqQLufPx5Xl/wOp+t+oxq5arxff/vuf+W+70WjyYOpZQqwJbsWELfWX05fOYw/Vr1Y1SnUZQOLO3VmDRxKKVUAXTs3DEGRAxg7ua51L+6PnN7z6VF7RbeDgvQxKGUUgWKMYYv1n7BS/Ne4mLiRd7p8A6vtHuFAL8Ab4eWShOHUkoVEHuP7aX3zN6s3LOSO+rcwZSnp3BD1Ru8HVYGmjiUUsrLkpKTGPPzGN5e8jYl/Eow+enJ9GzZ0+MT+XJLE4dSSnnRpv2bCJ8Rzo5DO3i48cOM7zKequWqejusLGniUEopL7gQf4E3F73JJ798QpWyVVjQdwEPNX7I22G5RROHUkrls6U7l9Ln6z4cOHmAvq368t+H/kvZkmW9HZbbNHEopVQ+iT0fy8BvBjJ742xuqHIDa15ZQ8s6Lb0dVo5p4lBKKQ8zxjDz95kMnjeYc3HnGP7AcF677zVK+Jfwdmi5oolDKaU8aF/sPnrP7M3yqOXcVvs2Pn/mc268+kZvh3VFNHEopZQHJKckM3b5WIYvHo6fjx8TnpxA7//rXWCH2OaEJg6llMpjWw9speeMnmz7ZxsP3vIgnz3xGdUrVPd2WHlGE4dSSuWRiwkXGb54OGOXjaVSmUp82+dbOjXulK97ZeQHTRxKKZUHft71M32+7sP+E/sJvyOc9x95n3Ily3k7LI/QxKGUUlfgxPkTDJ47mJnrZ1K3cl1WvbSKO6+/09theZQmDqWUygVjDLM3zGbgnIGciTvD0PZDeeP+Nwj0D/R2aB6niUMppXIo5kQMfWf1ZenOpTSr1YzPn/mcBtUbeDusfKOJQyml3JRiS+GTXz7hje/eQET4uPPH9LurH74+vt4OLV95bECxiFwjIitFJEpEdonIC1b5WyJyWES2W7f2Tq95TUSiRWSPiLR1Km9nlUWLyBBPxayUUpnZcXAHzUc1Z/Dcwdx1/V3sHrGbAfcMKHZJAzx7xZEMvGiM2SoipYEtIrLMem6sMWaMc2URuRHoDNQHrgaWi0hd6+nPgHuBQ8AmEVlsjNntwdiVUgqAuMQ43l7yNh9EfkDFUhWJCI/g8SaPF7khtjnhscRhjDkKHLXunxeRKKBaFi/pAHxjjEkA9otINNDUei7aGLMPQES+sepq4lBKedSKqBX0/ro30cej6XZ7N8Y8OoYKpSp4Oyyvy5e57yISCjQCNlhF/UXkDxGZLiLlrbJqwEGnlx2yyjIrV0opjzh18RTdv+zOPR/dgzGG5YOXM/3Z6Zo0LB5PHCISDMwHBhpjzgETgdpAQ+xXJB86qrp4ucmiPP379BKRzSKyOTY2Nk9iV0oVL8YY5myaQ7036zHj9xm82u5V/nzrT+6pd4+3QytQPDqqSkT8sSeNWcaYBQDGmGNOz38OLLEeHgKucXp5deCIdT+z8lTGmCnAFICwsLAMiUUppbJy8NRBnpv1HEv+WMKtNW8lcmAkDWs09HZYBZLHEofYe46mAVHGmI+cyqta/R8ADwE7rfuLgdki8hH2zvE6wEbsVxx1RKQWcBh7B/oTnopbKVW8pNhSmLByAq8vfB2bsfHhox8y4J4B+PnqbIXMePLM3A48DfwpItutsteBLiLSEHtzUwzQG8AYs0tE5mLv9E4G+hljUgBEpD8QCfgC040xuzwYt1KqmNh5eCc9v+rJhv0baHNjGyY9NYlaIbW8HVaBJ8YUvVadsLAws3nzZm+HoZQqoOKT4hn5w0hGLx1N2aCyjHt8HE82e7JYD7EFEJEtxpiw7OrptZhSqlhZ/fdqes3sxZ5/9/B086f58LEPCSkd4u2wChVNHEqpYuHMpTO8Ov9VpqyeQq2rahE5MJI29dt4O6xCSROHUqpIM8awYOsCno94nmPnjvFSm5d468G3KFWilLdDK7Q0cSiliqzDpw/Tb3Y/Fm1fRKMajfj++e+5teat3g6r0NPEoZQqcmw2G5N+ncSQBUNItiXz/iPvM6j1IB1im0f0LCqlipTdR3YTPiOcdf9bR+t6rZn01CRqV6rt7bCKFE0cSqkiISEpgf/+9F9G/TiK0oGl+bLblzzT4pliP8TWEzRxKKUKvbXRawmfEU7U0SieaPoEYx8fS6UylbwdVpGliUMpVWidizvHkAVDmLhqIjUq1ODHAT9yX4P7vB1WkaeJQylVKC3avojnZj3Hv2f/ZWDrgbzT4R2CA4O9HVaxoIlDKVWoHD1zlOcjnmf+1vncXP1mFj63kKa1mmb/QpVnNHEopQoFm83G1N+m8sq3rxCfFM9/O/2XF+99EX8/f2+HVuxo4lBKFXh7/t1Dr5m9WP33alpd34opT0+hTuU63g6r2NLEoZQqsBKTE3l/6fu8+8O7BAUEMfWZqXRv2V2H2HqZJg6lVIG0Yd8Ges7oyc7DO3ks7DE+7vwxVcpW8XZYCk0cSqkC5nz8eYYuHMqnKz+lWrlqLO6/mAduecDbYSknmjiUUgXGD3/8QN9ZfTl0+hD9WvVj5EMjKRNUxtthqXQ0cSilvO7YuWO88M0LzNk0hxur3sjaV9fSonYLb4elMqGJQynlNcYYvlz3JS/OfZGLiRd5u8PbvNruVQL8ArwdmsqCJg6llFdEH4+m98zerPhrBS2va8nnz3zODVVv8HZYyg2aOJRS+SopOYmPln3EW9+/RYBfAJOemkT4HeH4+Ph4OzTlJk0cSql8szlmMz2/6smOQzvo1LgT47uM5+pyV3s7LJVDmjiUUh53MeEiwxYNY9zycVQpW4WFzy2kY6OO3g5L5ZImDqWURy3duZQ+X/fhwMkD9LmzD6M7jaZsybLeDktdAU0cSimPiD0fy8BvBjJ742xuqHIDa15ZQ8s6Lb0dlsoDmjiUUnnKGMPM32cyeN5gzsWdY9j9w3i9/euU8C/h7dBUHtHEoZTKM/ti99Hn6z4s272M22rfxpSnp1C/Wn1vh6XymCYOpQqg2NhLxMScJTS0LCEhJb0dTraSU5IZt3wcwxYPw8/Hj0+f+JS+d/bVIbZFlMf+VUXkGhFZKSJRIrJLRF6wyiuIyDIR2Wv9LG+Vi4h8IiLRIvKHiDR2OlZXq/5eEenqqZiVKggiIqKoWXMK9947j5o1pxAREeXtkLK09cBWmo1qxsvfvsy99e5l94jd9LurnyaNIsyT/7LJwIvGmHpAc6CfiNwIDAF+McbUAX6xHgPcB9Sxbr2AiWBPNMBwoBnQFBjuSDZKFTWxsZfo0SOSuLhkzp5NJC4umR49IomNveTt0DK4lHCJl+e9TNNRTTly9gjz+szju37fUb1CdW+HpjzMY01VxpijwFHr/nkRiQKqAR2AVla1r4BVwKtW+QxjjAHWi0g5Ealq1V1mjDkFICLLgHZAhKdiV8pbYmLOEhDgQ1zc5TJ/fx9iYs4WqCar5buX0/vr3uyL3Uf4HeG89/B7lC+lf88VF/nSxyEioUAjYANQ2UoqGGOOikglq1o14KDTyw5ZZZmVK1XkhIaWJTHRlqYsKclGaGjBmPdw8sJJBs8dzIzfZ1CnUh1WvbSKO6+/09thqXzm8UZIEQkG5gMDjTHnsqrqosxkUZ7+fXqJyGYR2RwbG5u7YJXyspCQkkyb1pagID/KlAkgKMiPadPaev1qwxjDrPWzqDesHrM3zmZo+6H88dYfmjSKKY9ecYiIP/akMcsYs8AqPiYiVa2rjarAcav8EHCN08urA0es8lbpylelfy9jzBRgCkBYWFiGxKJUYdGlSz1at65ZYEZVxZyIoe+svizduZSmtZryyzO/0KB6A6/GpLzLk6OqBJgGRBljPnJ6ajHgGBnVFVjkVP6MNbqqOXDWatKKBNqISHmrU7yNVaZUkRUSUpImTap6NWmk2FIYt3wcN711E2v2ruHjzh+zbsg6TRrKo1cctwNPA3+KyHar7HVgNDBXRHoA/wCPWs/9CLQHooFLQDcAY8wpEXkH2GTVe9vRUa6U8owdB3cQPiOcTTGbaN+gPROenEDNijW9HZYqIDw5quo3XPdPANzjor4B+mVyrOnA9LyLTinlSlxiHG8veZsPIj+gQqkKRIRH8HiTx7E3IChlpzPHlVIArPxrJb1m9iL6eDTP3vYsYx4dQ8Xgit4OSxVAmjiUKuZOXTzFy/NeZvra6dQOqc3ywcu5p16GRgGlUmniUKqYMsYwb/M8no94npMXT/Jqu1cZ/sBwggKCvB2aKuA0cShVDB08dZDnZj3Hkj+WEFYzjMiBkTSs0dDbYalCQhOHUsVIii2FCSsn8PrC17EZGx8++iED7hmAn69+FSj36f8WpYqJnYd3Ej4jnPX71tO2flsmPjmRWiG1vB2WykP5tRy/Jg6lirj4pHjeXfIu70W+R7mgcszsMZMnmz2pQ2yLmIiIKHr0iCQgwIfERBvTprWlS5d6HnkvTRxKFWGr/15N+Ixw/j72N083f5qPHvuIq0pf5e2wVB5zXo7fsbJyjx6RtG5d0yNXHrrTilJF0JlLZ+g9szd3fnAnicmJRA6MZEaPGZo0iijHcvzOHMvxe4JecShVhBhjWLB1Af0j+nP83HFebPMiIx4cQakSpbwdmvKg/F6OX684lCoiDp8+zEMTHuKRSY9QpUwVNg7dyJhHx2jSKAbyezn+bK84RCQQeA5oiX0fjN+AicaYeI9EpJTKEZvNxuTVkxmyYAiJyYm89/B7DGo9CH8/f2+HpvJRfi7H705T1QzgPDDeER8wk8ur2iqlvCTqaBThM8JZG72We+rdw+SnJlO7Um1vh6W8JCSkZL4sxe9O4rjeGHOL0+OVIrLDUwEppbKXkJTA6J9GM/LHkQSXCOaLZ7+g621ddYityhfuJI5tItLcGLMeQESaAWs9G5ZSKjProtfRc0ZPoo5G0aVpF8Y9Po5KZSp5OyxVjLiTOJph35nvH+txDSBKRP7Evo3GzR6LTimV6lzcOYYsGMLEVROpUaEGPwz4gfYN2ns7LFUMuZM42nk8CqVUlhZtX0S/Wf04cvYIL9zzAu92fJfgwGBvh6WKKXcSRx1jzHLnAhHpaoz5ykMxKaUsR88c5fmI55m/dT4NqjVgft/5NLu2mbfDUsWcO4ljmIg8DLwEBANTgQRAE4dSHmKz2Zj22zRe/vZl4pPiGdlxJC+3fVmH2KoCwZ3EcSfwIrDdejzMGBPhuZCUKt72/LuHXjN7sfrv1bS6vhWTn5pM3Sp1vR2WUqncSRzlsXeQ/w+oDtQUETHGGI9GplQxk5icyAeRH/DOkncICghi6jNT6d6yuw6xVQWOO4ljPTDaGDNdRIKA97APx73No5EpVYxs2LeBnjN6svPwTh4Le4yPO39MlbJVvB2WUi65kzhaG2P+ATDGxAEDROT/PBuWUsXD+fjzvPHdG4xfMZ5q5aqxuP9iHrjlAW+HpVSW3EkcB0XkKeBaY8zbIlID0HWqlLpCS3Ys4bnZz3Ho9CH6terHqE6jKB1Y2tthKZUtdxLHBMAG3A28jX3dqvlAEw/GpVS+y69tN4+dO8YL37zAnE1zqH91fda+upYWtVt47P2UymtuzRw3xjQWkW0AxpjTIhLg4biUylf5se2mMYYv1n7BS/Ne4mLiRd7p8A6vtHuFAD/9dVKFizuJI0lEfLEvqY6IhGC/AlGqSMiPbTejj0fTe2ZvVvy1gjvq3MGUp6dwQ9Ub8uTYSuU3dxLHJ8BCoJKIjAQeAd7waFRK5SPHtpuOpAGXt9280sSRlJzEh8s+ZMT3IwjwC2DSU5MIvyMcHx/dQ00VXtkmDmPMLBHZAtwDCNDRGBPl8ciUyiee2nZzc8xmen7Vkx2HdtCpcSfGdxnP1eWuvqJjKlUQZPpnj4hUcNyA40AEMBs4ZpVlSUSmi8hxEdnpVPaWiBwWke3Wrb3Tc6+JSLSI7BGRtk7l7ayyaBEZktsPqoqm2NhLbNp0lNjYS7k+hqttN8eObUVMzNlcHfdiwkUGzxlMs1HNOH7+OAv6LmB+3/maNFSRIZlNABeR/dj7NQT7UuqnrfvlgH+MMbWyPLB9rscFYIYx5iar7C3ggjFmTLq6N2JPTE2Bq4HlgGONhb+Be4FDwCagizFmd1bvHRYWZjZv3pxVFVUE5HWHtmNU1datxxg0aFWujhu5M5I+X/ch5mQMfe7sw+hOoylb8squXJTKLyKyxRgTll29TK84jDG1jDHXApHAA8aYq4wxFYH7gQXZHdgYsxo45Wa8HYBvjDEJxpj9QDT2JNIUiDbG7DPGJALfWHVVMefcoX32bCJxccn06BF5xVceoaFlGTRoVY6PG3s+lqemPkW7j9sR6B/ImlfWMPGpiZo0VJHkTg9dE2PMj44HxpifsC98mFv9ReQPqymrvFVWDTjoVOeQVZZZuSrmHB3azhwd2vl5XGMMM3+fSb1h9Zi7eS7D7h/G9mHbaVmn5RXFoVRB5k7iOCEib4hIqIjUFJGhwMlcvt9EoDbQEDgKfGiVu1rFzWRRnoGI9BKRzSKyOTY2NpfhqcLCUx3aOTnuvth9tB3XlmemP0PdynXZ9uY2RnQYQQn/ElcUg1IFnTuJowsQgn1I7ndAJassx4wxx4wxKcYYG/A59qYosF9JXONUtTpwJItyV8eeYowJM8aEhYSE5CY8VYi46tCeNq3tFQ+fdee4ySnJjIkcw01v3cT6fev57InP+O2V36hfrf6VfiylCoVMO8czVBQpA9iMMRfcPrhIKLDEqXO8qjHmqHV/EPZZ6Z1FpD72EVuOzvFfgDrYrzj+xj4U+DD2zvEnjDG7snpf7RwvPjy1TEhmx12xYx39vunDXyf+5IFbHmDCExOoXqF6nr2vUt7kbud4tvM4RKQBMAOoYD0+AXQ1xuzM5nURQCvgKhE5BAwHWolIQ+zNTTFAbwBjzC4RmQvsBpKBfsaYFOs4/bF30PsC07NLGqp4CQkp6ZF1pdIf91LCJR5//3mWxHyJJJYiYFNXOt/1iiYNVSxle8UhIuuAocaYldbjVsAoY0yB3Y9DrzhUXlq+ezk9v+rFgVP7YU9T2PAfSCxJUJAfBw708uiCiErlpysejuuklCNpABhjVgGlriA2pQqFkxdO8uz0Z7l37L3YUoRSq/rDmkch0Z4o8mIUl1KFkTuJY5+IvGmNqgoVkTeA/Z4OTClvMcYwe8Ns6g2rx6yNs3i9/ev8NngjtsO109TLi1FcShVG7iSO7thHVS3APrIqBOjmyaCU8pYDJw/wn0/+w5NTn6TWVbXY8sYWRj40khpXV/TIKC6lCiOPjqryFu3jUDmVYkvh0xWfMvS7oQCMemgU/e7qh6+Pb5p6+bXZk1Le4PVRVUoVFn8c+oPwGeFs3L+R9g3aM/HJidSoWMNlXU+N4lKqMHFnP47JwOB0o6qmAAV2VJVSWXFcNVSuFsDE9WMY8/MYKpSsQER4BI83eRwRVwsWKKUc3EkcGUZViYiOqlJecaVNRY4VdX2qRXOp8RxMmRN0u70bYx4dQ4VS2e4W4JGYlCps3Ekc+0TkTWCm9fgpdFSV8oIrXUY9NvYS3fsuID5sEVy/Cc5WJGB5X94bNYYKpXL3hZ8fe5UrVdBktZGTI1GsIe2oqqvQUVUqn13pMurGGD7/ZQYJD7wHdbbAjrtgwWACT92Q67kYnljaXanCIKsrjltFpCbQFbgL+7pRjiFY2gis8tWV7Av+z8l/eG7Wc/zw5w/IxerwYw84ZV+d/0rmYnhyr3KlCrKsEsckYClwLeA8ttWRQK71YFxKpZHZcuenT8cTG3vJ5Rd1ii2FCSsn8PrC17EZGx8++iGVYu+l19Jf8C/jQ1KS7YrmYnhqaXelCjp31qqaaIzpm0/x5Amdx1E0OfoT/P19iI9PxhhDyZL+LvsWdh7eSfiMcNbvW0+bG9sw6alJ1Aqx73acl53ZzjE5EpH2cajCyt15HG5PACxMNHEUXbGxl9i27RgdOnxHfHxKarljwcHS5XwY+cNI3lv6HmWCyjDu8XE82exJjw6x1VFVqqjIswmAShUkISElKV8+kBIlfNMkDn9/Hxb8HsnYTa+x5989PNX8KT567CNCSnt+Uy+dFKiKG00cqtDJ0LcQEMfFRgvo8/3vhFYMZekLS2l7U1vvBahUEefOIodKFSjO27sG1YuCR8Zgq7OBwfcOZueInZo0lPIwveJQhdL/tStDq5Er+Gn399xU9Wa+6D6NsNBsm2aVUnlAE4e6IvndMWyz2Zi8ejJDFgwhMTmR0Z1GM/jewfj7+Xv8vZVSdpo4VK7lxXIbOUk8UUej6DWjF79F/8bdN9zN5Kcnc12l667kIyilckETh8oV5+U2HDOne/SIpHXrmgBuJQN3E09CUgKjfxrNqJ9GUSqgFNOfnc6ztz2rq9gq5SWaOFSuZLbcxuTJ2xk1amO2ySCrxOOcbNZFryN8Rji7j+6mc5POjOs8jsplKnv64ymlsqCjqlSuuFpuIzExhVGjNrq16J8j8ThzrPMEcC7uHP1m9aPl+y25kHCBJc8vIaJXhCYNpQoATRwqV5yHxDr24B46tHmWycBZVus8Ld6+mBuH3cjEXycy4O4B7Bqxi//c/B+Pfh6llPt0yRF1RZw7twFq1pxCXFxy6vOOpUBc9XWkX+dpzITGrIwfz7dbvqVBtQZM7TqVprWa5ttnUaq40yVHVL5wJARH8pg2rW2GRf8y6yDv0qUerVvXZN/+06w98R2vL32Q+KR4RnYcycttX87xEFtdM0qp/KGJQ+VI+i9nVyOjDhzo5fYX+KmUg7zyay9W/72aVte3YvJTk6lbpW6O49Kd+JTKP9pUpdyW/st57NhWDBq0yu2mKWeJyYl8EPkB7yx5h6CAIMY8MobuLbvnaohtbOylHDWRKaVc06YqlSuZNfe4Gj77wgsrCQjwTfN6d3bA27BvAz1n9GTn4Z08euujfNLlE6qUrZLjmBx0Jz6l8pfHRlWJyHQROS4iO53KKojIMhHZa/0sb5WLiHwiItEi8oeINHZ6TVer/l4R6eqpeJX9iqJmzSnce+88atacQkREVOpzmQ2fTUxMSVOW1Q545+PPMyBiAC1Gt+DMpTMs7r+YuX3mZpk0sorJQXfiUyp/eXI47pdAu3RlQ4BfjDF1gF+sxwD3AXWsWy9gItgTDTAcaAY0BYY7ko3KW85XFK7mYLj6ck5JMXz88d1phuRm1hn+wx8/UH94fT5d+SnPtXqOXSN28cAtD1xRTA6uhgZfyZawSqmseaypyhizWkRC0xV3AFpZ978CVgGvWuUzjL3DZb2IlBORqlbdZcaYUwAisgx7MorwVNzFVXbNPY4vZ1fbpHbqVCfTpqRj547xwjcvMGfTHOpfXZ+1r66lRe0WLmNw7O4H0KhR5Rw1QTlGaOmoKqU8L7/7OCobY44CGGOOikglq7wacNCp3iGrLLNylcfcae5x9eWcWf+DMYYv133Ji3Nf5GLiRd7u8DavtnuVAL8Al+8fERHFs8/+lBqDv78wfnzrDDElJKQQHOz6GLoTn1L5o6DMHHc1lMZkUZ7xACK9RGSziGyOjY3N0+CKA+fmnuBgf0qU8GXs2FYZvohDQkrSpEnV1KG4rvofoo9H0/qj1nT/sjv1r67PjmE7ePP+NzNNGo4mKeckkZRkGDhwBWPHtrJv2BRk/xvHxwduvXWmy74OpVT+yO/EccxqgsL6edwqPwRc41SvOnAki/IMjDFTjDFhxpiwkBDP7zNdFHXpUo+xY1uRlGQjIMCXQYNWZfoF7ar/oXvPH3nz23do8FYDNh/YzKSnJvHry79yQ9UbsnzfmJiz+Phk/BvB11do3LgyW7Y8jc1m/3shLi4lyzWwlFKel9+JYzHgGBnVFVjkVP6MNbqqOXDWatKKBNqISHmrU7yNVabyQGzsJTZtOpr6BRwbe4lBg1b98fKPAAAWNElEQVSRkJDC+fOXO6N//nl/hi/pDKOsrjpI4n1jeTdyGPfddB9Rb0fR+87e+Phk/18sNLRsamJwlpJiCA0ty4ULiQQGuh72q5TKfx7r4xCRCOyd21eJyCHso6NGA3NFpAfwD/CoVf1HoD0QDVwCugEYY06JyDvAJqve246OcpU7jj6JrVuPMWjQqjQzra+7rlyGzui4uGQ6dVqEzUaa2dipfSJ+iXBrJNRfgy2+NF88GcGzrTrnKCZHM1n6Po7p09ulNpXpcFulCg6dOV6MOGZ++/oKFy4kpXkuKMiPLVue5tZbZ6aZgZ2+jvNs7CGffs77a4dggk/hu7cFk7qPpeczzXIdn/OoqmuuKcOFC4kZljZJP6JLKZV3dOa4SsO5T8IVf38fLlxITB1y6+MjXLyYNrn4+oq9eSjwIoPnDubrHV9Tp9b1vNR8Gg+1aHPFI5pCQkrSpk0tIiKi6NhxUYZ1p3S4rVIFgyaOYiIm5ix+fpn3Nziafpo0qUrr1jXZtu0YHTp8R3z85ZnhFy4kMmHZdL6fPZZzcecYdv8wwpsO5Oih+NQ6V7pCbXY7A2rCUMr7NHEUE/Y+iRSXz5Uo4ZtmprXjL/9x4+6iT5/l9krBp6DlfL7c/zdhNZrxZfdp/LHah7q1v0q9MujR4yamTdtJQIAPCQkpDB3anN69b8nRl31hX3dKl3ZXxUFBmcehPMB51FRISEk+/vjuDHVKlPBl27ZnXPYXNG5cmeAyvtDgV3h4DFQ6QODWR/i07XwqBdTKMBz300+3pz6Oj0/hzTfXZrq+VGYK87pT7qyrpVRRoImjiHL1Jda79y1MmtSaEiV8CQ72JyjIjy++aEe9ehVdHuN8QAwX7/kImi2BI3Vg/ktI1O1cW6t8tk1fDjmdc1FY151yd10tpYoCbaoqRNxtBnHVT9Ct21IaNqxE794N6dSpbobjOB+7VBkY8f0IPlz2IWWqliPul66UONqIZJtJ/RJfsOBvzp9PdCtud5qanN+/MHaEF/YmNqVyQhNHIeHODneOL9/9+88CGdd4atjwKz755B56974FIHUC3fLlB1KPHVf+L8rdv4Tj8QfpeUdP3n/4fZIvlciwPtWgQasyxNi9+03Mnh2VpkMdsm9qyuyzFaYv3MLcxKZUTuk8jkLAnR3uHF++YIiLc90J7tC2bU1+/fUQJUr4kphoIzk5hSSfC9Dse6i7BTkXwpc9pvJMmwddvv7dd3/nzTfXpikLDvZnxYrHCA0ty+TJOxg1aoNbcy6K0u59OtdEFXY6j6MIya4ZJLs5GulFRh4AsK4MDNTeDs0XQYk42H43Zltrei36H/5fRLm8qhk5cn2GYzqWBwkJKckbb7Sgd+9b3GpqKkpNPIWxiU2p3NDEUQi4agZJTEzh9On41Oap9F++bgk+DbcvgGv+guPXwE+PwKmrAUhISUkzf8Jh27Zj+PpmXJDw9debpann7pyLotbEo3NNVHGgo6oKiPQLDjpLP9LI31+w2QyPPfY9NWtOYevW4xm+fLMkNqi/xj7Etso+ZMOD+EcOSE0aDukXEnTM6L54Me2VTWCgb2q/SU4V1lFUShVn2sdRALjT8Q2X13Lq2HFRmmapwEBfHn/8er755i9EJEPndBoVjkDLb6HSQfjnBljXCb/4Cixc2JFHHllMQsLl1zr3Nbjqi3C89/Tp7a64LV8nzinlfe72cegVh5dlN/4//SS+8uUDM8yfiI9P4auvdpOQYOOBB2qze3c35sy5nxdfvJWgIGs5ct8kCPsJOn4MpU/Biifg5+5woTzJyYZHHllMeHiDTP/yz7CMOlCqlD+LFnXMkw5g5w2ilFIFm/ZxeFlWncPOw2QdVyLnziVkOX9i3ry/GTHidh577AbuuqsGEybsgKrR0HI+lD0Bf4fBhvshoVSa1yUkpDBt2k62bHk6zaq0Dq76Imw2Q6NGldOU6ZWDUkWfJg4vy6xzODg4IMMkvu7dlyLiajfdtJYvj6FevYr4lUyg2Uu/serofDhbEX7sZZ8BngnHCrlNmlTN8JyjLyL9cFPn5OBuk5tSqnDTPo4CwNX4f19f4dlnl2boyxAh23kaASWEPqOFOQff58SFE7Sq+CRrx9cn/kLWLZPp50+4unrIrMxV30thnY+hVHGl8zgKEefx/8HBAQwfvpZ58/7OUC/LTm+HUmdIvG0hn+zaTY2S9fBb8hQbz19DfLqNm1wZO/auNBMKu3dfiq+vkJJiUjvA0w83dSQ9Hx8ydJwX1vkYSqmsaeIoIEJCSrJ8+QG6d1/qXoLIwAb1focmP4GPDf+tD/Lv7jtIjIcELieN0qUDSExMsb7oL79PcLA/jRtXSr166Nr1R5KSLl+Ndu36Y4Y5HdlNPCzM8zGUUpnTxFFAREWdpFu3pWmGw7qt3L9wx7dQ+QAcqgu/PUzShQoZqgUH+zN+/N00bVqVW2+dmea5lBTD1q3HuPPOOYiYNEkDICnJsG3bMdq0qZValtnEw1Kl/LE5LYiolCpaNHEUABERUblLGj7J0PAXuGUlgb6lSFn3BEm7GwKuO9BTUgzt21/rsqN77NhWDBq0yu1lS8B1x35goC8LFjxIo0aVNWkoVURp4vAyR3NPjpNG5f32q4xyx2lUuh2RI2aw7fcLdOqUcWa3Q48eN6V+madfVym7ZUsCAnwyDL3NbKSV81WJUqro0cThZTExZ0lOzkHS8I+Dpj9CvfVwvjz81JNth69nwfWH6dSpLrYsVh6ZNm0nw4bdlmaLWOerAlfLlpQs6YcxZNrspAv7KVX8aOLwssREW4b+hEyF/gktvoOg8/Dn/8GWtpAcAMCAASuoVatsapNTSootQyLIapRTZs1XjRtXzjYh6MJ+ShUvmjg8LLuZ1L/+ejD7g5Q8C7cthNBdcLIqLHsWTlyTpkpioo1OnRZjsxleeaUpo0ZlXPo8MTEly1FOevWglHKHJg4Pcmcm9Zo1h7I4gs3eJNXkJ3tH+Mb29isN4+uy9sWL9mG3I0asc/n80KHNs00GevWglMqOJg4PcbXvd/r9LXr0WMrSpTGuD1DumH0V2yoxcPg6WPswnLsqQ7XAQF+35n34+0uulz5XSilnmjg8ZPLkHVnOpI6KOsn06TszvtAnGW5ZaR9mmxQAvz4Ge8NwNcR2+PAW3Hbb1RmW+nDFnTWulFLKHbqsugdktr2q8659GzcezfjCSjHw0Di49WeIaQDzX4a9TchsXsbo0Rs5eTI+dSOk0qUDMo0pIMA3zaZMSimVW3rF4QExMWcpUSJjE1JSko2OHb8jOdkQHt7g8hP+8dDEGmJ7oSxEdoeD2a8qm5Bg3971wIFeHDjQi5iYs2zdeszlRL4LF5LYuvW4y5VvlVIqJ7xyxSEiMSLyp4hsF5HNVlkFEVkmInutn+WtchGRT0QkWkT+EJHG3og5J1zNqAb7zO24uBSSkmz2fTIAauyCR8bYk8au22H+S24lDQcfH/jxx30ANGlSld69G3LgQC+GDGmSoe6gQStdbk2rlFI54c2mqruMMQ2dlvAdAvxijKkD/GI9BrgPqGPdegET8z3SHAoJKUmPHjdlXSnoHNw9E9p8CQlBsLg/rO8AySVy9F4XLybz/PMrqFlzChERUanv36lT3QxNV+n3EFdKqdwoSE1VHYBW1v2vgFXAq1b5DGPfOGS9iJQTkarGGBedBN7nWF126tQ/M6lh4PqN0PQH+3aum9rBH60yHWLrDseOgM6jtkJDy5KcnHGDKF2tVil1pbyVOAzws4gYYLIxZgpQ2ZEMjDFHRaSSVbca4DxL7pBVViASR1TUSZYvP0DlyiU5fTqeQYNW4eOTyd4ZZWLt60tV3QdHr4U1j8C5kFy/t2N+iIPzqC13duxTSqnc8FbiuN0Yc8RKDstE5K8s6roaUpRhjQ4R6YW9KYsaNWrkTZTZeP755Xz66fbsK/okw82/QsPlkOJvTxh7mpC+pdDX1z76Kasd/vz9ffDxEd59tyXDhq0FLieO9FcUOhNcKeUJXkkcxpgj1s/jIrIQaAocczRBiUhV4LhV/RDgvL5GdeCIi2NOAaaAfetYT8YPsHbtYfeSRsg/cMc8qPAv7LsZfu8AcWVcVvXxEZKTXYfeo0cDunW7iYAAn9QkUK1acLZXFDoTXCmV1/I9cYhIKcDHGHPeut8GeBtYDHQFRls/F1kvWQz0F5FvgGbAWW/3b0RERNG1609ZV/JLgLClUH8tXCoDPz8L/9TP8iVBQf489lhdpk7NODHwxRfDqFevYpoyvaJQSnmDN644KgMLrZnMfsBsY8xSEdkEzBWRHsA/wKNW/R+B9kA0cAnolv8hX+bYqS8pKYv1y6+JgtsXQKmzENUcNrWHpMBsjx0Xl8SsWX/h7y9pVszt379hhqThoFcUSqn8lu+JwxizD8iwaJIx5iRwj4tyA/TLh9Cy5bjSyDRpBJ2H5ouh9nY4XRm+fw6Oh2Z5TB8fKFnSn+RkGzabSTNxLyDAhxUrHuf226vl4adQSqkrU5CG4xZYjiG23btndqVhoM5maPY9+CfCljaw4y6wZX96bTb7UiSdO9/AwoV704ySCgz0IyBAV4VRShUsmjiy4VgaXcRkMsT2BNw+H6pFw7+17COmzlbKWC8LiYk2ZszYnaFc510opQoiTRxZcF4aPQNJgQarofHP9iuL3zrBX83Ii8n4wcH+pKQYnXehlCqQNHFkISbmLH5+LqaRVDxkH2J71RGIuQnWdYRLeXNlULp0AOPH30379tdq0lBKFUiaOLKQmGjjwoWkywV+iXBrJNRfA3GlYdkzcKBB5gfIgp+f4OsrJCSk7TNJTrZp0lBKFWiaODKRYVZ4tT3Qcj6UPn15iG1iUI6PW7KkH8bAtGltad26JpMn72DkyPUEBPjqsiBKqUJB7KNdi5awsDCzefPmXL8+KuokN974hf1B4AX7ENvrtsGZEHvn97Frc3XcwEBfFi3qSKNGldMkh9jYSzqJTynldSKyxWnF8kzpFYcLo0dvAAxct9WeNPwTYGtr2HG3fa2pHCpVyg+bzX6V0aZNrQzP6yQ+pVRhookjnaiok8z+fi20+xaq/w3HasBvj8LpKtm+tm3bmqxefTjNKKygID8WLOiQ4SpDKaUKK51d5uTrWX/S4InuJD/4AVQ6YB8ttaSfW0kDYPDgsNT9v8uUCSAoyC/1KkOThlKqqNArDkts7CV6DppDyn9+hMN1Yd1DcLFcjo+jCw8qpYo6TRyWmJizBCZWImHhIDgbguttQDIXEOBDo0aVAe2zUEoVbZo4LKGhZe1bsNrcWy7E1xd8fHwICPDFZtNZ3kqp4kMTh+XEiThsWayU7hAY6MP7799J5871ALRJSilV7GjisGzc6N7eUCI+dO5cLzVRaMJQShU3mjgsTZtWdVnu5wfJyfbJeyKiTVJKqWJPE4elXr2K9O/fMM0yI92730SfPrcQHBzAhQuJ2iSllFLokiMZREWdZOPGozRtWjXT7VqVUqoo0iVHcqlevYqaMJRSKgs6c1wppVSOaOJQSimVI5o4lFJK5YgmDqWUUjmiiUMppVSOFMnhuCISCxzwdhwedhVwwttBFBB6Luz0PNjpebgsp+eipjEmJLtKRTJxFAcistmd8dbFgZ4LOz0PdnoeLvPUudCmKqWUUjmiiUMppVSOaOIovKZ4O4ACRM+FnZ4HOz0Pl3nkXGgfh1JKqRzRKw6llFI5oomjABGR6SJyXER2OpVVEJFlIrLX+lneKhcR+UREokXkDxFp7PSarlb9vSLS1Ruf5UqIyDUislJEokRkl4i8YJUXq3MhIoEislFEdljnYYRVXktENlifaY6IBFjlJazH0dbzoU7Hes0q3yMibb3zia6MiPiKyDYRWWI9Lq7nIUZE/hSR7SKy2SrL398NY4zeCsgN+D+gMbDTqex9YIh1fwjwnnW/PfATIEBzYINVXgHYZ/0sb90v7+3PlsPzUBVobN0vDfwN3FjczoX1eYKt+/7ABuvzzQU6W+WTgL7W/eeASdb9zsAc6/6NwA6gBFAL+B/g6+3Pl4vzMRiYDSyxHhfX8xADXJWuLF9/N7x+EvSW4T9FaLrEsQeoat2vCuyx7k8GuqSvB3QBJjuVp6lXGG/AIuDe4nwugJLAVqAZ9gldflZ5CyDSuh8JtLDu+1n1BHgNeM3pWKn1CssNqA78AtwNLLE+V7E7D1bcrhJHvv5uaFNVwVfZGHMUwPpZySqvBhx0qnfIKsusvFCymhkaYf9ru9idC6t5ZjtwHFiG/a/kM8aYZKuK82dK/bzW82eBihSB8wCMA14BbNbjihTP8wBggJ9FZIuI9LLK8vV3QzdyKrzERZnJorzQEZFgYD4w0BhzTsTVR7NXdVFWJM6FMSYFaCgi5YCFQD1X1ayfRfI8iMj9wHFjzBYRaeUodlG1SJ8HJ7cbY46ISCVgmYj8lUVdj5wLveIo+I6JSFUA6+dxq/wQcI1TverAkSzKCxUR8ceeNGYZYxZYxcXyXAAYY84Aq7C3U5cTEccffc6fKfXzWs+XBU5R+M/D7cCDIhIDfIO9uWocxe88AGCMOWL9PI79j4mm5PPvhiaOgm8x4Bjx0BV7e7+j/Blr1ERz4Kx1iRoJtBGR8tbIijZWWaEh9kuLaUCUMeYjp6eK1bkQkRDrSgMRCQJaA1HASuARq1r68+A4P48AK4y9AXsx0NkabVQLqANszJ9PceWMMa8ZY6obY0Kxd3avMMY8STE7DwAiUkpESjvuY/8/vZP8/t3wdkeP3tJ0cEUAR4Ek7H8R9MDeNvsLsNf6WcGqK8Bn2Nu8/wTCnI7THYi2bt28/blycR5aYr9s/gPYbt3aF7dzAdwMbLPOw05gmFV+LfYvvGhgHlDCKg+0Hkdbz1/rdKyh1vnZA9zn7c92BeekFZdHVRW782B95h3WbRcw1CrP198NnTmulFIqR7SpSimlVI5o4lBKKZUjmjiUUkrliCYOpZRSOaKJQymlVI5o4lAqB0RknbdjUMrbdDiuUkqpHNErDqVyQEQuWD9bicgqEflWRP4SkVnWjHdEpImIrLP20dgoIqXFvrfGF9Y+CttE5C6r7rMi8p2IfC8i+0Wkv4gMtuqsF5EKVr3aIrLUWthujYjc4L2zoIo7XeRQqdxrBNTHvsbPWuB2EdkIzAEeN8ZsEpEyQBzwAoAxpoH1pf+ziNS1jnOTdaxA7LN4XzXGNBKRscAz2NdlmgL0McbsFZFmwATsazYple80cSiVexuNMYcArKXPQ7Ev4X3UGLMJwBhzznq+JTDeKvtLRA4AjsSx0hhzHjgvImeB763yP4GbrVWCbwPmOa0QXMLDn02pTGniUCr3Epzup2D/fRJcL0+d6Zrw6Y5jc3pss47pg33viYa5D1WpvKN9HErlrb+Aq0WkCYDVv+EHrAaetMrqAjWwL7SXLeuqZb+IPGq9XkTkFk8Er5Q7NHEolYeMMYnA48B4EdmBfde+QOx9Er4i8if2PpBnjTEJmR8pgyeBHtYxdwEd8jZypdynw3GVUkrliF5xKKWUyhFNHEoppXJEE4dSSqkc0cShlFIqRzRxKKWUyhFNHEoppXJEE4dSSqkc0cShlFIqR/4fjIMsliJJHCYAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "engel.plot.scatter(x='income', y='foodexp', color='darkblue')\n", "plt.plot(X_pred[:, 1], y_pred, color='darkgreen')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare with the linear regression (OLS) line." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "betas_lr = lm.LinearRegression(fit_intercept=False).fit(X, y).coef_" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmcjWUbwPHfNZsZOxkSMSrZS4wtelPWtCgteFtkX5OlkIoQScpWtlARkyWyZqIo8dp3hkyMbDH2bcx27veP5xnOmH3MmfX6fj7zmXPuc5/nuc9TzjXPvVy3GGNQSimlksstoxuglFIqa9HAoZRSKkU0cCillEoRDRxKKaVSRAOHUkqpFNHAoZRSKkVcFjhExFtENovILhHZJyJD7PIyIrJJRA6JyFwR8bLLc9nPg+3X/ZyO9Z5dflBEmriqzUoppZLmyjuOcOBJY8zDQFWgqYjUBj4FxhhjygIXgPZ2/fbABWPMA8AYux4iUhFoBVQCmgITRcTdhe1WSimVCJcFDmO5aj/1tH8M8CSwwC7/Dnjeftzcfo79egMREbv8B2NMuDHmCBAM1HRVu5VSSiXOw5UHt+8MtgEPAF8BfwMXjTFRdpXjQAn7cQngGIAxJkpELgF32eUbnQ7r/J54FSlSxPj5+aXRp1BKqZxh27ZtZ40xvknVc2ngMMZEA1VFpCCwCKgQXzX7tyTwWkLlsYhIJ6ATQKlSpdi6dWuq2qyUUjmViBxNTr10mVVljLkIrAVqAwVFJCZglQRO2o+PA/cC2K8XAM47l8fzHudzTDXG+Btj/H19kwyYSimlUsmVs6p87TsNRMQHaAgEAWuAl+xqbYDF9uMl9nPs138zVgbGJUAre9ZVGaAssNlV7VZKKZU4V3ZVFQe+s8c53IB5xphlIrIf+EFEPgZ2ANPt+tOBWSISjHWn0QrAGLNPROYB+4EooLvdBaaUUioDSHZMq+7v7290jEMppVJGRLYZY/yTqqcrx5VSSqWIBg6llFIpooFDKaVUimjgUEqpbOKnHT8xfd30pCveIQ0cSimVxZ28eJIXJ73ICxNfYPr66TgcDpeez6Urx5VSSrmOw+Fg2p/T6LegH+FR4XzS4hP6NuqLm5tr7wk0cCilVBZ08N+DdJzZkXWH1lG/XH2mvj6VssXKpsu5NXAopVQWEhEVwaiVoxi2fBh5vPIwvc102tZti5VMPH1o4FBKqSxi498b6TCzA/tO7uMV/1cY33o8xfIXS/d2aOBQSqlM7sqNKwxcOJCv1n5FiYIlWNpjKc88/EyGtUcDh1JKZWLLdi2j6+yunLh4gu71uzOixQjyeefL0DZp4FBKqUzo9OXT9Azoybyt86h0TyXmdZ5HnfvrZHSzAA0cSimVqRhj+Gb9N7wz/x2uRVxjWPNh9GvaDy8Pr4xu2k0aOJRSKpM4dPoQnWd1Zs3BNTxW9jGmvj6V8sXLZ3Sz4tDAoZRSGSwyKpLRv4xm6LKh5PLIxZTXp9ChXgeXL+RLLQ0cSimVgbYc2ULHmR3ZdXwXL1Z7kQmtJ1C8YPGMblaiNHAopVQGuHrjKh8u/pDxv47n7gJ3s7DrQl6o9kJGNytZNHAopVQ6W7l3JV2+78LRc0fpWr8rn7zwCQVyF8joZiWbBg6llEonoVdC6fVDL+ZsnkP5u8uzrt866pWtl9HNSjENHEop5WLGGGb9bxZ95vfhcthlBj87mPeeeo9cnrkyummpooFDKaVc6HDoYTrP6szqoNU8ev+jfP3G11S8p2JGN+uOaOBQSikXiIqOYszqMQxeMhgPNw8mvjqRzv/pnGmn2KaEBg6llEpj249up8PMDuz4ZwfPPfwcX/33K0oWLpnRzUozGjiUUiqNXAu/xuAlgxmzagxF8xdlQZcFtKjWIl33ykgPGjiUUioN/LLvF7p834UjZ4/Q8bGOjHppFAVzF8zoZrmEBg6llLoDZ6+cpc+8PszaOIsHiz3I2nfW8ni5xzO6WS6lgUMppVLBGMOcTXPoNbcXF8Mu8n6z9/ngmQ/w9vTO6Ka5nAYOpZRKoZCzIXSd3ZWVe1dSq0wtvn7ja6qUrJLRzYIzZ+D4cahWzaWn0cChlFLJFO2IZvyv4/ngpw8QEca1Gkf3J7rj7uaesQ27fh3GjIFPP4V774W9e8GFA/Ium1AsIveKyBoRCRKRfSLytl3+kYicEJGd9k8zp/e8JyLBInJQRJo4lTe1y4JFZICr2qyUUgnZdWwXtUfUps+8PjxR7gn2D9lPzwY9MzZoREfDjBlQtix88AE0aAA//ujSoAGuveOIAvoaY7aLSD5gm4issl8bY4wZ7VxZRCoCrYBKwD3AahF50H75K6ARcBzYIiJLjDH7Xdh2pZQCICwijKHLhvJZ4GfclecuAjoG0LJGy4ydYmsM/Pwz9O9v3V3Urg1z50K99Ml75bLAYYw5BZyyH18RkSCgRCJvaQ78YIwJB46ISDBQ034t2BhzGEBEfrDrauBQSrnUb0G/0fn7zgSfCaZt3baMfnk0hfMUzthGbdsG/frBb7/BAw/AggXQooXL7zKcpcvadxHxAx4BNtlFPURkt4jMEJFCdlkJ4JjT247bZQmVK6WUS5y/dp5237ajwRcNMMawus9qZrw5I2ODRkgIvPoq+PvD7t0wfjzs2wcvvpiuQQPSIXCISF7gR6CXMeYyMAm4H6iKdUfyeUzVeN5uEim//TydRGSriGwNDQ1Nk7YrpXIWYwxzt8ylwocVmPm/mfRv2p89H+2hQYUGGdeo8+ehb18oVw4WLYL334e//4a33gIvrwxpkktnVYmIJ1bQmG2MWQhgjDnt9PrXwDL76XHgXqe3lwRO2o8TKr/JGDMVmArg7+8fJ7AopVRijp0/RrfZ3Vi2exnVS1cnsFcgVUtVzbgG3bgBX34Jw4fDpUvQti0MGQIlMz7nlcsCh1gjR9OBIGPMF07lxe3xD4AXgL324yXAHBH5AmtwvCywGeuOo6yIlAFOYA2g/9dV7VZK5SzRjmgmrpnIwEUDcRgHn7/8OT0b9MTDPYNWKzgcMGeOdWfxzz/w1FPWNNsqmWCdiM2VV6Yu8DqwR0R22mUDgdYiUhWruykE6AxgjNknIvOwBr2jgO7GmGgAEekBBALuwAxjzD4XtlsplUPsPbGXDt91YNORTTSu2JjJr02mjG+ZjGvQr7/Cu+/Cjh3WIr5vvoEnn8y49iRAjMl+vTr+/v5m69atGd0MpVQmdSPyBsOXD2fkypEU8CnA2JZjebXWqxk3xXbPHmtq7c8/Q+nSVvdU69aQznt3iMg2Y4x/UvV05bhSKkf5468/6DSrEwf/PcjrtV/n81c+xzefb8Y05vhxGDQIvv0WChSA0aOhe3fwztz5rjRwKKVyhIvXL9L/x/5M/WMqZYqUIbBXII0rNc6Yxly6ZI1bjBljjWn07QvvvQeFM3iNSDJp4FBKZWvGGBZuX8hbAW9x+vJp3mn8Dh899xF5cuVJ/8ZERMDkyTB0KJw7Z63L+Phj8PNL/7bcAQ0cSqls68SFE3Sf053FOxfzSKlHWPrWUqqXrp7+DTHGWuH93nvWGownn4RRo6B6BrQlDWjgUEplOw6Hg8m/T2bAwgFEOaIY9dIoejfsnTFTbNetg3fegc2boXJlWLECmjZN99XeaUkDh1IqW9l/cj8dZ3Zkw98baFihIZNfm8z9Re9P/4YcOAADBsDixXDPPVYW2zfeAPcMTsGeBjRwKKWyhfDIcD75+RNGrBhBPu98fNv2W96o80b6T7H991/46COYNg1y57bGMHr3th5nExo4lFJZ3vrg9XSc2ZGgU0H8t+Z/GdNyDEXzF03fRly9Cp9/Dp99BuHh0K0bfPgh+GbQVF8X0sChlMqyLoddZsDCAUxaO4lShUuxoucKnqryVPo2IirK6oYaPNi623jpJRgxwtpcKZvSwKGUypIW71xMt9nd+PfSv/Rq2IthzYeR1ztv+jXAGFi61FrxfeAA1K0LCxdCnTrp14YMooFDKZWlnLp4ircC3uLH7T/yUMmHWNRtETXL1Ez6jWlp0yYrp9S6dfDgg1a68+bNs/RMqZTQwKGUyhIcDgfT/pxGvwX9uBF5g09afELfRn3x9PBMv0b8/be1FmP+fChaFCZOhA4dwDMd25AJaOBQSmV6B/89SKdZnfjjrz+oX64+U1+fStli6TiGcPYsDBsGkyZZQWLQIGttRr586deGTEQDh1Iq04qIimDUylF8vPxjfLx8mPbGNNrVa5d+U2zDwmDcOPjkE2vWVIcO1lTb4sXT5/yZlAYOpVSmtOnwJjrM7MDeE3t5xf8VxrUax90F7k6fk0dHw6xZ1nTa48fh2Wdh5EioWDF9zp/JaeBQSmUqV25c4f1F7/Plmi8pUbAES3os4dmHn02fkxsDv/wC/frB7t1QowZ8/z08/nj6nD+L0MChlMo0lu9eTtfZXTl+4Tjd63dn+AvDye+TP31OvmOHFTBWr4YyZSAgAF55Jd03U8oKNHAopTLc6cunefuHt5m7ZS4Vi1dkff/11Lk/ndZDHD1qdUl9/z0UKgRjx0KXLpArV/qcPwvSwKGUyjDGGL7d8C195/XlWsQ1hjYfSv+m/fHy8HL9yS9etFZ4jx9vPe/Xz0pKWLCg68+dxWngUEpliOAzwXSe1ZnfDvxGvQfq8fUbX1O+eHnXnzg83Fp/8fHHcOECvP66NdW2VCnXnzub0MChlEpXkVGRfLHqCz5a+hFeHl5Mfm0yHR/riJurxxIcDpg3DwYOhCNHoHFja/vWqlVde95sSAOHUirdbA3ZSofvOrDr+C5aVGvBhNYTuKfgPa4/8dq1VoqQrVvh4YchMNAKHCpVdLqAUsrlroVfo++8vtQaUYvQq6Es6raIH7v+6PqgsW+ftQbjiSfg9Gn47jvYtk2Dxh3SOw6llEut3LuSLt934ei5o3R5vAsjW4ykQO4Crj3pyZNWmvMZMyBvXmvl99tvg4+Pa8+bQ2jgUEq5ROiVUHr90Is5m+dQ/u7yrOu3jnpl67n2pJcvWxspff65tU9Gz57w/vtQpIhrz5vDaOBQSqUpYwyz/jeLPvP7cDnsMoOeGcTAZgPJ5enCdRGRkTB1KgwZAqGh0LKlNdX2vvtcd84cTAOHUirNHA49TJfvu7Bq/yoevf9Rpr4+lUolKrnuhMZYe2EMGACHDlmpQT77zEoVolxGA4dSmVBo6HVCQi7h51cAX9/cGd2cJEVFRzF29VgGLRmEh5sHX/73S7o+3tW1U2w3bLBmSm3YABUqwJIl8MwzOWYzpYzksv+qInKviKwRkSAR2Scib9vlhUVklYgcsn8XsstFRMaLSLCI7BaRak7HamPXPyQibVzVZqUyg4CAIEqXnkqjRvMpXXoqAQFBGd2kRG0/up1aI2rx7oJ3aVShEfuH7Kf7E91dFzQOHYIXX7S2aj182Oqi2r3bmj2lQSNduHI6bhTQ1xhTAagNdBeRisAA4FdjTFngV/s5wFNAWfunEzAJrEADDAZqATWBwTHBRqnsJjT0Ou3bBxIWFsWlSxGEhUXRvn0goaHXM7ppcVwPv86789+l5oianLx0kvld5vNT958oWbika0545gz06GGlNg8MhKFDITgYOnYED+08SU8uu9rGmFPAKfvxFREJAkoAzYH6drXvgLVAf7t8pjHGABtFpKCIFLfrrjLGnAcQkVVAUyDAVW1XKqOEhFzCy8uNsLBbZZ6eboSEXMpUXVar96+m8/edORx6mI6PdeTTFz+lUB4X/T137RqMGWOt8g4Lg06drKm2xYq55nwqSekSpkXED3gE2AQUs4MKxphTIlLUrlYCOOb0tuN2WULlSmU7fn4FiIhwxCqLjHTg5+fidQ/JdO7qOfrM68PM/82kbNGyrH1nLY+Xc9FeFVFR8O231jatp07B889bmymVK+ea86lkc/nKcRHJC/wI9DLGXE6sajxlJpHy28/TSUS2isjW0NDQ1DVWqQzm65ub6dOb4OPjQf78Xvj4eDB9epMMv9swxjB742wqDKrAnM1zeL/Z++z+aLdrgoYxsHy5lUOqY0coXRrWrbNmT2nQyBRceschIp5YQWO2MWahXXxaRIrbdxvFgTN2+XHgXqe3lwRO2uX1bytfe/u5jDFTgakA/v7+cQKLUllF69YVaNiwdKaZVRVyNoSus7uycu9Kapapya9v/EqVklVcc7KtW62ZUmvXwgMPwIIF0KKFDnpnMq6cVSXAdCDIGPOF00tLgJiZUW2AxU7lb9izq2oDl+wurUCgsYgUsgfFG9tlSmVbvr65qVGjeIYGjWhHNGNXj6XyR5VZd2gd41qNY8OADa4JGkeOwH//a62/2LcPvvwS9u+3Zk9p0Mh0XHnHURd4HdgjIjvtsoHASGCeiLQH/gFetl9bATQDgoHrQFsAY8x5ERkGbLHrDY0ZKFdKucauY7voOLMjW0K20KxKMya+OpHSd5VO+xOdP2/ti/HVV+DubqU8798f8qfTdrEqVVw5q+pP4h+fAGgQT30DdE/gWDOAGWnXOqVUfMIiwhi6bCifBX5G4TyFCegYQMsaLZG0/qv/xg2YMMFKC3LpErRta6ULKemiqbwqTenkZ6UUAGsOrKHTrE4EnwnmzUffZPTLo7kr711pexKHA+bMsRIP/vMPPPWUNc22iovGTJRL6H4cSuVw56+dp/237Xny8ycxxrC6z2q+aftN2geNX38Ff39rq9YiRaznK1Zo0MiC9I5DqRzKGMP8rfN5K+Atzl07R/+m/Rn87GB8vNJ4z4rdu61xi5Urram1s2dDq1bg6q1ilcto4FAqBzp2/hjdZndj2e5l+Jf2J7BXIFVLpfHe28ePw4cfWrvuFSxo7ZHRrRt4e6fteVS608ChVA4S7Yhm4pqJDFw0EIdx8PnLn9OzQU883NPwq+DSJWvcYswYa0yjTx9rTKOQppjLLjRwKJVD7D2xl44zO7Lx8EaaVGrCpFcnUca3TNqdICICpkyxkg+ePQuvvQbDhoGfX9qdQyUqvdLxa+BQKpu7EXmDj5d9zKeBn1LQpyCz2s/i1Vqvpt0UW2OsFd7vvQd//w1PPmltplStWtLvVWkmICCI9u0D8fJyIyLCwfTpTWjduoJLzqWBQ6ls7I+//qDjzI78dfovXq/9Ol+88gVF8qXh/tvr1lkpQjZtgsqV4eefoUkTXe2dzpzT8cdkVm7fPpCGDUu75M5DpzUolQ1dvH6RzrM68/hnjxMRFUFgr0Bmtp+ZdkEjKAiaN4f//AeOHYMZM2DnTmjaVINGBohJx+8sJh2/K+gdh1LZiDGGhdsX0iOgB2cun6Fv474MeW4IeXLlSZsT/PuvtRfGtGmQJw8MHw69ekHuzLNXSE6U3un4NXAolU2cuHCC7nO6s3jnYqreW5Vlby2jeunqaXPwq1dh9GjrJzwcune3ptr6+qbN8dUdiUnH3759IJ6ebkRGOlyajj/JwCEi3kA3oB7WPhh/ApOMMTdc0iKlVIo4HA6m/DGFAQsHEBEVwacvfkrvhr3x9PC884NHRcH06dZdxunT8NJL8MknVspzlamkZzr+5NxxzASuABNi2gfM4lZWW6VUBgk6FUTHmR1ZH7yeBhUaMOW1Kdxf9P47P7AxsHSpteL7wAGoVw9++glq177zYyuX8fXNnS6p+JMTOMoZYx52er5GRHa5qkFKqaSFR4Yz8ueRDF8xnLy58vLNm9/Q5tE2aTPFdtMma6bUunXWjnuLFlkD4TrorWzJCRw7RKS2MWYjgIjUAta7tllKqYRsCN5Ah5kdCDoVROuarRnbcixF8xe98wP//be1FmP+fChaFCZNgvbtwTMNurxUtpKcwFELa2e+f+znpYAgEdmDtY3GQy5rnVLqpsthlxmwcACT1k6iVOFSLO+5nGZVmt35gc+etVZ4T5pkBYlBg+CddyBfvjs/tsqWkhM4mrq8FUqpRC3euZjus7tz8tJJ3m7wNh8//zF5vfPe2UHDwmDcOGuw++pV6+5iyBAoXjxtGq2yreQEjrLGmNXOBSLSxhjznYvapJSynbp4ircC3uLH7T9SpUQVfuz6I7Xuq3VnB42Ohlmz4IMP4MQJeO45GDkSKrgmPYXKfpKzcnyQiEwSkTwiUkxElgLPurphSuVkDoeDr//4mgqDKrBs9zKGPz+cbR9su7OgYQwEBlo5pNq2hXvugbVrYfFiDRoqRZJzx/E40BfYaT8fZIwJcF2TlMrZDv57kE6zOvHHX39Qv1x9prw2hQfvfvDODrpjB/TrB6tXw333wdy58PLLOlNKpUpyAkchrAHyv4GSQGkREWOMcWnLlMphIqIi+CzwM4YtG4aPlw/T3phGu3rt7myK7dGj1grv77+39sMYOxa6dgUvr7RruMpxkhM4NgIjjTEzRMQH+BRrOu6jLm2ZUjnIpsOb6DCzA3tP7OUV/1cY12ocdxe4O/UHvHDBGvQeP956/u671lTbggXTpsEqR0tO4GhojPkHwBgTBvQUkf+4tllK5QxXblzhg58+YMJvEyhRsARLeizh2YfvYAgxPBy++go+/hguXoQ33rA2VipVKu0arXK85ASOYyLyGnCfMWaoiJQCNE+VUndo2a5ldJvTjeMXjtO9fndGtBhBPu9Urp1wOKxxi4EDISQEGje2tm+tmsb7iCtF8gLHRMABPAkMxcpb9SNQw4XtUirdpde2m6cvn+btH95m7pa5VLqnEuv7r6fO/XVSf8A1a6yuqG3b4OGH4ZdfoFGjtGuwUrdJ1spxY0w1EdkBYIy5ICI6sqaylfTYdtMYwzfrv+Gd+e9wLeIaw5oPo1/Tfnh5pPKf0759VhLC5cvh3nth5kx49VVw0/3ZlGslJ3BEiog7Vkp1RMQX6w5EqWwhPbbdDD4TTOdZnfntwG88VvYxpr4+lfLFy6fuYCdPWmnOZ8yw0oJ8+in07Ane3mnSVqWSkpzAMR5YBBQVkeHAS8AHLm2VUukoZtvNmKABt7bdvNPAERkVyeerPmfI0iF4eXgx+bXJdHysI26puSu4fBk++ww+/9zaJ6NnT2v191133VEblUqpJAOHMWa2iGwDGgACPG+MCXJ5y5RKJ67adnNryFY6fNeBXcd30aJaCya0nsA9Be9J+YEiI+Hrr+GjjyA0FFq1smZN3Z8G+24olQoJ/tkjIoVjfoAzQAAwBzhtlyVKRGaIyBkR2etU9pGInBCRnfZPM6fX3hORYBE5KCJNnMqb2mXBIjIgtR9UZU+hodfZsuUUoaHXU32MmG03fXw8yJ/fCx8fD8aMqU9IyKVUHfda+DX6zO1DrRG1OHPlDAu7LuTHrj+mPGgYAwsXQqVK1latFSpYe2UEBGjQUBnLGBPvD3AEOGz/jgbOAufsx0cSep/T+/8DVAP2OpV9BLwTT92KwC4gF1AGa5W6u/3zN3Af4GXXqZjUuatXr25U9jdnzn7j4zPGFCgwzvj4jDFz5uy/o+OdOXPNbN580kyevCPVx125Z6Xx6+9n6IDpMquLuXjtYuoas369MY8+agwYU7GiMUuXGuNwpO5YSiUTsNUk8f1qjEn4jsMYU8YYcx8QCDxrjClijLkLeAZYmIyA9AdwPpnxqznwgzEm3BhzBAgGato/wcaYw8aYCOAHu67K4ZwHtC9diiAsLIr27QPv+M7Dz68AvXuvTfFxQ6+E8tq012g6rinent6s67eOSa9NokDuFHZ3/fUXvPgi1K0LR45YXVS7dsEzz2heKZVpJGeEroYxZkXME2PMz1iJD1Orh4jstruyCtllJYBjTnWO22UJlascLmZA21nMgHZ6HtcYw6z/zaLCoArM2zqPQc8MYuegndQrWy9lJz5zxuqOqljRWocxbBgcOgQdOoBHcuawKJV+kvN/5FkR+QD4HmtK7mtYXVapMQkYZh9nGPA50A5r0P12hvgDW7zJFUWkE9AJoJSmV8j2XDWgnZLjHg49TJfvu7Bq/yrq3F+Hr1//mkolKqXshNeuwZgx1pTasDDo1Mmaalus2J18DKVcKjl3HK0BX6wpuT8BRe2yFDPGnDbGRBtjHMDXWF1RYN1J3OtUtSRwMpHy+I491Rjjb4zx9/X1TU3zVBYS34D29OlN7nj6bHKOGxUdxejA0VT+qDIbD2/kq/9+xZ/9/kxZ0IiOhunT4cEHrey1jRtbC/omTtSgoTK/5AyEWGMm5AfyJre+/R4/Yg+OF3d63BtrXAOgErEHxw9jDYx72I/LcGtwvFJS59XB8ZwjZkD7zJlr6XLcX3euN+UHVDF0wDw74Vlz7NyxlB3Y4TBm2TJjKlWyBr7r1DHmzz/TsOVKpR7JHBxPsqtKRKoAM4HC9vOzQBtjzN4k3hcA1AeKiMhxYDBQX0SqYnU3hQCd7eC1T0TmAfuBKKC7MSbaPk4PrAF6d2CGMWZfUm1WOYevb26X5JW6/bjXw6/TctRbLAv5FonIg9eWNrR6oh8lC5dM/kG3brU2U1qzBh54ABYsgBYtdNBbZTliktiPSUQ2AO8bY9bYz+sDI4wxmXY/Dn9/f7N169aMbobKJlbvX02H7zpx9PwROFgTNj0NEbnx8fHg6NFOSQeuI0fg/fet9Re+vtYYRqdO4OmZPh9AqWQSkW3GGP+k6iVnjCNPTNAAMMasBfLcQduUyhLOXT3HmzPepNGYRjiihTxre8C6lyHCChRJzuI6dw769IFy5eCnn6zgERxszZ7SoKGysOQEjsMi8qGI+Nk/H2AtClQqWzLGMGfTHCoMqsDszbMZ2Gwgf/bZjONE7NXaCc7iunHDyil1//0wbpy1mdKhQ1aakPz50+lTKOU6yQkc7bBmVS3EmlnlC7R1ZaOUyihHzx3l6fFP8+q0VylTpAzbPtjG8BeGU+qeu5KexeVwWHt7lytnjWXUrQs7d8K0aVBClx+p7CM5SQ4vYG0Xmx9wGGOuur5ZSqWvaEc0X/72Je//9D4A41qNo/sT3XF3c79Zp3XrCjRsWDr+zZ5Wr7aCxY4dUL06fPMNPPlken8MpdKFy2ZVKZVV7D6+m44zO7L5yGaaVWnGpFcnUequ+BeRxpnFtXu3FTACA8HPD+bMgZYtdTMlla0lZ+X4FKDPbbOqpgKZdlaVUomJ2SK2WAkvJm0czehfRlM4d2ECOgb22RFuAAAgAElEQVTQskZLJDnTY48ftxbuffcdFCwIo0dDjx6QK5frP4BSGSw5gSPOrCoR0VlVKkPc6b7gMVvEupUI5nq1uZj8Z2lbty2jXx5N4TxJ7hYAly5Z6UHGjLHSnr/zDmc7vs2Ri274XY5GkxaonEBnVaksIyAgiNKlp9Ko0XxKl55KQEDK9hMLDb1Ou64LCfOfw7X6X2GMwWt1Vz59+sukg0ZEBIwfb82U+uQTeOklOHiQgEfaUurh+aluk1JZUWIbOc2yH64j9qyqIuisKpXO7jSNujGGr3+dSfizn0LZbbDrCVjYB+/z5RNfi2EMzJ9vZa19+22oWhW2bYNZswjN7Zvmqd2VygoS66qqLiKlgTbAE1gZbGOWmWuOBJWu7mRf8H/O/UO32d1Yvmc5cq0krGgP563psYlm1F23Dt5919p1r3JlWLECmja9mSLElXuVK5WZJRY4JgMrsXbfc87fERNA7nNhu5SKJaF05xcu3CA09Hq8X9TRjmgmrpnIwEUDcRgHn7/8OUVDG9Fp5a945ncjMtIRf0bdAwegf39YssRafzFjhrWIz909VjVXpXZXKrNLTq6qScaYrunUnjShuaqyp5iBbU9PN27ciMIYQ+7cnkREWAGgdesKN+vuPbGXjjM7svHwRhpXbMzk1yZTxrcMkMgA+7//wkcfWQv2cueG996zuqdyJ3z34NymmEDk3A6lspLk5qpKMnBkRRo4sq/Q0Ovs2HGa5s1/4saN6JvlMQkH8xV0Y/jy4Xy68lPy++RnbMuxvFrr1cSn2F69Cp9/bqUJCQ+Hbt3ggw9I7hSpO53ppVRmkdzAoXtSqizF1zc3hQp5kyuXe6zA4enpxsL/BTJmy3sc/Pcgr9V+jS9e+QLffIl8+UdFWZspDR4Mp09bM6U++cRKeZ7CNmnAUDmJBg6V5cQZW/AK49ojC+my9H/43eXHyrdX0qRyk4QPYIw1ftG/Pxw8CPXqWdlra9d2feOVygY0L4LKcpy3d/WpEAQvjcZRdhN9GvVh75C9iQeNTZvg8cfh+eet54sWwR9/aNBQKgX0jkNlSf9pmp/6w3/j5/1LqVz8Ib5pNx1/v0S6ZoODYeBAa01GsWIweTK0bw8e+k9AqZTSfzXqjqT3wLDD4WDKH1MYsHAAEVERjGwxkj6N+uDpkcDGSGfPwrBhMGkSeHnBoEHwzjuQL5/L26pUdqWBQ6VazFRULy+3eKfEJkdKAk/QqSA6zezEn8F/8mT5J5ny+hQeKJrAQPb169YmSiNHwrVr0KGDNQhevHiK2qeUiksDh0oV5xQgMSun27cPpGHD0gDJCgbJDTzhkeGM/HkkI34eQR6vPMx4cwZvPvpm/FNso6Nh5kwrc+2JE/Dcc1bwqKBrK5RKKxo4VKoklG5jypSdjBixOclgkFjgcQ42G4I30HFmR/af2k+rGq0Y22osxfIXi9sgY6w9Mfr1gz17oEYNa2+M//wnrT+6UjmezqpSqRJfuo2IiGhGjNicrKR/MYHHWUyeJ4DLYZfpPrs79UbV42r4VZa9tYyATgHxB40dO6BRI3jqKatbau5ca/aUBg2lXEIDh0oV5ymxMXtwv/9+7USDgbPE8jwt2bmEioMqMun3SfR8sif7huzj6YeejtuIo0fh9dehWjVrb++xYyEoCF555WYiQqVU2tOUI+qOOA9uA5QuPZWwsKibr8ekAolvrOP2PE+jJ1ZjzY0JLNi2gColqjCtzTRqlqkZ96QXLlgrvMePtwJEr17WYr6CBV32OZXKCTTliEoXMQEhJnhMn94kTtK/hAbIW7euQMOGpTl85ALrz/7EwJXPcSPyBsOfH867Td6NO8U2PBy++go+/hguXrQy1g4bBvfeC2jOKKXSiwYOlSK3fznHNzPq6NFOyf4CPx99jH6/d+KPv/6gfrn6THltCg/e/WDsSg6HNW4xcCCEhECTJtb2rQ8/fLNKWkwNVkolj3ZVqWS7/ct5zJj69O69NtldU84ioiL4LPAzhi0bho+XD6NfGk27eu3iTrFdu9baTGnrVmv3vVGjrIFwJ6Gh11PURaaUip92ValUSai7J77ps2+/vQYvr9ibGyVnB7xNhzfRYWYH9p7Yy8vVX2Z86/HcXeDu2JX27bPGLZYvJ7pESY5+NJ58XdriWyxvnOPpTnxKpS+XzaoSkRkickZE9jqVFRaRVSJyyP5dyC4XERkvIsEisltEqjm9p41d/5CItHFVe5V1R1G69FQaNZpP6dJTCQgIuvlaQtNnIyKiY5UltgPelRtX6BnQkzoj63Dx+kWW9FjCvC7zYgeNEyesVd4PPQR//snO1n0pcq4n1cYYSpeZFqtNMXQnPqXSlyun434LNL2tbADwqzGmLPCr/RzgKaCs/dMJmARWoAEGA7WAmsDgmGCj0pbzHUV8azDi+3KOjjaMG/dkrCm5CQ2GL9+9nEqDK/Hlmi/pVr8b+4bs49mHn71V4fJla/OksmVh1izo1Yuzm/bw6E8luXhDEl0XEt/U4MQG5ZVSd8ZlXVXGmD9ExO+24uZAffvxd8BaoL9dPtNYAy4bRaSgiBS3664yxpwHEJFVWMEowFXtzqmS6u6J+XKOb5vUFi3KJjgYfvryad7+4W3mbplLpXsqsb7/eurcX+dWhchI+Ppra8vW0FBO1X+W4LbvUP4p/xR1QcXM0NJZVUq5XnqPcRQzxpwCMMacEpGidnkJ4JhTveN2WULlKo0lp7snvi/nhMZEjDF8u+Fb+s7ry7WIawxtPpT+Tfvj5eEVU8HaC2PAADh0iNMVa/L8hVfZuLYErN2Mp+cWJkxoGKdN4eHR5M3rFe9n0J34lEofmWXleHzLfE0i5XEPINJJRLaKyNbQ0NA0bVxO4NzdkzevJ7lyuTNmTP04X8S+vrmpUaP4zam48Y2JBJ8JpuEXDWn3bTsq3VOJXYN28eEzH94KGhs2WLvuvfgieHpy6fsFlDncio1Rt/4miIw09Or1G2PG1Lc2bPKx/sZxc4Pq1WfFO9ahlEof6R04TttdUNi/z9jlx4F7neqVBE4mUh6HMWaqMcbfGOPv65vIPtMqQa1bV2DMmPpERjrw8nKnd++1CX5Bxzcm0q7DCj5cMIwqH1Vh69GtTH5tMr+/+zvli5e33vTXX1awqFsXjhyxuqh27eKvBx/FzT3u/4ru7kK1asXYtu11HA7r74WwsOhEc2AppVwvvQPHEiBmZlQbYLFT+Rv27KrawCW7SysQaCwihexB8cZ2mUoDoaHX2bLl1M0v4NDQ6/TuvZbw8GiuXLk1GP3LL0fifEnHmWVV5BgRT43h48BBPFX5KYKGBtH58c64ubnB6dPQvTtUrAi//AJDh8KhQ9bsKQ8P/PwK3AwMzqKjDX5+Bbh6NQJv7/in/Sql0p/LxjhEJABrcLuIiBzHmh01EpgnIu2Bf4CX7eorgGZAMHAdaAtgjDkvIsOALXa9oTED5Sp1YsYktm8/Te/ea2OttH7ggYJxBqPDwqJo0WIxDgexVmPfHBPxiIDqgVBpHY4b+fjm1QDerN/KevO1a/DFF9aivbAw6NzZ2oGvWOwMtzHdZG+++fPNMQ1PT2HGjKY3u8p0uq1SmYeuHM9BYlZ+u7sLV69GxnrNx8eDbdtep3r1WbFWYN9ex3k19oAvv2bU+gGYvOdxP1SHye3G0OGNWhAVBd9+awWJU6fghRespITlyiXavtDQ6+zYcRqAe+/Nz9WrEXFSm9w+o0splXZ05biKxXlMIj6enm5cvRpxc8qtm5tw7Vrs4OLuLlb3kPc1+szrw/e7vqdsmXK8U3s6L9RpjG8RH1i+3FrxvW8f1KkD8+dbYxrJ4Oubm8aNyxAQEMTzzy+Ok3dKp9sqlTlo4MghQkIu4eGR8JBWTNdPjRrFadiwNDt2nKZ585+4cePWyvCrVyOYuGoGS+eM4XLYZQY9M4iONXtx6vgNPHZuhxEfwtq13Lj3PsJnzKHAm61SvC9GUjsDasBQKuNllum4ysWsMYnoeF/Llcs91krrmL/8x4594lalvOeh6TS+PfIBZQo/wI4Pd1D+eisalfmCo3WfoVDjx7iyaSe9PV+k5KUe3N3tXz4evjHFM5+S2hkws7t9woFS2ZEGjmzM+UvM1zc348Y9GadOrlzu7NjxRrzjBdWqFSNvfneo8ju8OBqKHsV7+0t82eRHil3NT+jrXdkZ/gnNIvfwMQ0oEdaXsZG1OXc5mhs3ovnww/Vxcl4lJSvnnUos15dS2YkGjmwqvi+xzp0fZvLkhuTK5U7evJ74+HjwzTdNqVDhrniPccUrhGsNvoBay+BkWfjxHbz316TS8m8p6F+Z7tF/MIvqlKU/H9KUK3jHOUZK11xk1bxTSeX6Uio70TGOLCS5O9zFN07Qtu1KqlYtSufOVWnR4sE4x3E+dp78MGTpED5f9Tn5ixck7Nc2eJ98mJcitjK+wBryDD3F0SqP8fSemuzj7gTbESM5Kc6dz58VB8I1tbvKSTRwZBHJ2eEu5sv3yJFLQNwcT1Wrfsf48Q3o3NnaOS9m3GD16qM3jx1W6AAFn1nGmRvH6PBYB0a9OAr5eT2e7w8gz1/7oHR1Lk6dToVXDxJG7Bla7dpVZs6coFgD6pB0V1NCny0rfeFm5S42pVJK13FkAcnZ4S7myxcMYWHxD4LHaNKkNL//fpxcudyJiHAQFRVNpNtVqLUUHtyGXPbl2/bTeONuP+jXDwIDwc8PRoyAli35eMQmPvxwfaxj5s3ryW+/vYKfXwGmTNnFiBGbkrXmIjvt3qdrTVRWp+s4spGkukGSWqNxu8DAowD2nYGB+3dC7cWQKwx2PkmJ7f6w4GOM2YoULAijR0OPHpArF6Gh1xk+fGOcY8akB/H1zc0HH9Shc+eHk9XVlJ26eLJiF5tSqaGBIwuIrxskIiKaCxdu3Oyeuv3LN1nyXoC6C+HeA3DmXvIvf5YBF4PoxRcIMMajPm9smkWRsrey1u7YcRp397hrMwYOrBXrizK5ay6yWxePrjVROYHOqsokEpv/f/tMI09PweEwvPLKUkqXnsr27WfifPkmShxQaZ01xfbuw3j97xl6La/G3xe/5T3WsICHKEc/huR+niMXb/0vErOi+9q12Hc23t7uN8dNUiqrzqJSKifTMY5MIDkD33Arl9Pzzy+O1S3l7e1Oy5bl+OGHA4hInMHpWAqfhHoLoOgxOFqOl/+owCfh67ifc/wmZXnHNGMHJYHYYw3xjUXEnHvGjKZ33Jef3BljSinXSe4Yh95xZLCk5v/fvoivUCHvOKlDbtyI5rvv9hMe7uDZZ+9n//62zJ37DH37VsfHx05H7h4J/j/D8+Mg33keW96QjauuMy/8J67hybOeHVnU7SsO+PjF+5d/fCu68+TxZPHi59NkANh5gyilVOamYxwZLLHBYedpsjF3Ipcvh3PlSkSCx5s//y+GDKnLK6+U54knSjFx4i4oHgz1foQCZym/sxIjt0fR3LGa4xTgTV5hFtVxRLrhM2Mf27a9HisrbYz4xiIcDsMjj8ROka53Dkplfxo4MlhCg8N583rFWcTXrt1KJBlJA1evDqFChbvwyB1OrXf+ZO2pHyl2uhBDlpenw/X9XMOL93iKcdQjjFv7d8dkyK1Ro3icY8aMRdw+3dQ5OCS3y00plbXpGEcmEN/8f3d34c03V8YZyxAhyXUaXrmELiOFucdGEXYhlAmHKtPizwPkIopJ1GEYDTlL3jjvu339RHx3DwmVxTf2klXXYyiVU+k6jizEef5/3rxeDB68nvnz/4pTL9FB7xh5LhLx6CK+2rOfd4/czdu/e3O3YzfzeYj3eIq/KZLgW8eMeSLWgsJ27Vbi7i5ER5ubA+C3TzeNCXpubsQZOM+q6zGUUonTwJFJ+PrmZvXqo7RrtzJ5ASIOB1T4H/iv4LkT0Xw6Oy/lw/9lHWVozjNsphQA+fJ5ERERbX/R3zpP3ryeVKtW9ObdQ5s2K4iMvHU32qbNipt7YsRIauFhVl6PoZRKmAaOTCIo6Bxt264kPDwVQaPgv/DYAmrJUT5b5sNjFyIIIjfP8RJLqQhY4yJ583oyYcKT1KxZnOrVZ8U6RHS0Yfv20zz++FxETKygARAZadix4zSNG5e5WZbQwsM8eTxxOIyux1Aqm9LAkQkEBASlLmi4RUHVX7m/zG+M2u5GiyPwL+50pgXTqUk07rGqR0cbmjW7L96B7jFj6tO799pkpy2B+Af2vb3dWbjwOR55pJgGDaWyKQ0cGSymuyfFQaPYEYrUnMeHf5+l6yLBw8eLv1/tSt1F93D6evz/Wdu3r3zzy/z2vEpJpS3x8nKLM/U2oZlWznclSqnsRwNHBgsJuURUVAqChmcYPtWX0itqCwN+htxRwjRqkWfQUJq2rcvlhVOB+O8apk/fy6BBj8baItb5riC+tCW5c3tgDAl2O2liP6VyHg0cGSwiwhFnPCEhbqV38cbdCxi2+wYlr8NPUoH3eJoDFMPrwz0srfrAzS6n6GhHnECQ2CynhLqvqlUrlmRA0MR+SuUsGjhcLKmV1L//fizpg/hcpEm5WYw6+g8PHYVNHsX4Ly1YZ+67WSUiwkGLFktwOAz9+tVkxIi4qc8jIqITneWkdw9KqeTQwOFCyVlJvW7d8USO4OARv5WMuvE7DXc6CPbKzSs8z/yoqsTMlHJ27VokAEOGbIj3aO+/XzvJYKB3D0qppGjgcJH49v1u3z4w1lqI9u1XsnJlSLzvL5X/L4bnm8NrIdc46+lGz1wNmBzekMjb/pN5e7sna92Hp6ekOvW5Uko508DhIlOm7Ep0JXVQ0DlmzNgb530F5QoDfb+j59mjOK7BJ4XK8+mF1lwi7l3A4MF1ePTRe+Kk+ohPcnJcKaVUcmhadRdIaHtV5137Nm8+Fes1L6Lok3cRf3sMo++ZowQULcSDHm8z8EL7eIMGwMiRmzl37sbNjZDy5fOKtx6Al5c7ISGX7uyDKaUUesfhEiEhl8iVK24XUmSkg+ef/4moKEPHjlUAEBy0dN/KCI8llLkazs93e9DfPMOef+smeZ7w8Gjatw/k6NFOHD3aiZCQS2zffjrehXxXr0ayffuZeDPfKqVUSmTIHYeIhIjIHhHZKSJb7bLCIrJKRA7ZvwvZ5SIi40UkWER2i0i1jGhzSsS3ohqsldthYdFERjqYOHEX9Qlms9coAqLnczFPOI3KV6LZ2Y/YczrpoBHDzQ1WrDgMQI0axencuSpHj3ZiwIAacer27r0m3q1plVIqJTKyq+oJY0xVpxS+A4BfjTFlgV/t5wBPAWXtn07ApHRvaQr5+uamffvKCb5eiX9Z5jaFNUzB1+scr/kXoLpbD1YfeBOicqXoXNeuRfHWW79RuvRUAgKCbp6/RYsH43RdxYyxKKXUnchMXVXNgfr24++AtUB/u3ymsTYO2SgiBUWkuDHmVLxHyWAx2WWnTdsT57V7uMQQAmnLFi57wDsPufFleAPCtzUA4x7P0ZInZkdA51lbfn4FiIqKu0GUZqtVSt2pjAocBvhFRAwwxRgzFSgWEwyMMadEpKhdtwTgvEruuF2WKQJHUNA5Vq8+SrFiublw4Qa9e6/FzS323hn5uEE/1tKH33F3i2JsRRhRvDTnN7WEy76pPnfM+pAYzrO2krNjn1JKpUZGBY66xpiTdnBYJSIHEqkb3zzSODk6RKQTVlcWpUqVSptWJuGtt1bz5Zc7E3zdkyg6sYnBrMKXa8wpI7z/cC5Cgp6FVTW4vafQ3d2a/ZTYDn+enm64uQkff1yPQYPWA7cCx+13FLoSXCnlChkSOIwxJ+3fZ0RkEVATOB3TBSUixYEzdvXjwL1Oby8JnIznmFOBqWBtHevK9gOsX38ikaBhaMEeRvIzZTnLmiJevFsXtl2pAoHNISx/vO9ycxOiouJvevv2VWjbtjJeXm43g0CJEnmTvKPQleBKqbSW7oFDRPIAbsaYK/bjxsBQYAnQBhhp/15sv2UJ0ENEfgBqAZcyenwjICCINm1+jve1RznCZyznUY6yzyc3T9eDFXd5w4b/wj+VEj2uj48nr7zyINOmxV0Y2LevPxUq3BWrTO8olFIZISPuOIoBi+yVzB7AHGPMShHZAswTkfbAP8DLdv0VQDMgGLgOtE3/Jt8Ss1NfZGTsgeeyhDKSFbRgLyfdfGhf3YfvKl8n+mAdWNsMIr2TPHZYWCSzZx/A01NiZczt0aNqnKARQ+8olFLpLd0DhzHmMBAnaZIx5hzQIJ5yA3RPh6YlKeZOwzloFOUKg1hFZzYRhgcflLmbMf/5l+tXi8HydnDGL9FjurlB7tyeREU5cDhMrIV7Xl5u/PZbS+rWLeGqj6SUUimWmabjZloxU2zbtbt1p5GbCPrwO/34HR8imVKgDEManSA0XyjsbAy7ngBH0pfX4bBSkbRqVZ5Fiw7FmiXl7e2Bl5dmhVFKZS4aOJIQkxpdxHDjRjTuRNOWrQzhF+7hMj96PMjAOjf4q9zf8G8ZWPUSXCqa9IGdREQ4mDlzf5xyXXehlMqMNHAkwjk1OhieJohPWUElTrOBUrx8f0U2PLbVurP4swUcqEVaLMbPm9eT6Gij6y6UUpmSBo5EhIRcwsNDqM4xPmM5T/A3f1GEFnmfZlGD7eC7EUIqw4bn4Xra3Bnky+fFhAlP0qzZfRo0lFKZkgaOREhICFOufEtrdnKGPHR3e5ap1S8QVWUFhOWDVW/A0SqpOraHh+DuLoSHx56dFRXl0KChlMrUNHDE59w5fmvQmXq7fqIibgyjAZ8VL8GV/yyFfBcgqDZsaQYRPik+dO7cHhgD06c3oWHD0kyZsovhwzfi5eWuaUGUUlmCWLNdsxd/f3+zdevWlL8xLAwmTCD64+Fw5QozqMHgXPU4VWcNPLADLvrCupfg9H2pape3tzuLFz/PI48UixUcQkOv6yI+pVSGE5FtThnLE6R3HM5CQ2HQIHb7PsxrV+qx/4GTUHsyeIbD9oaw60mI9kzxYfPk8cDhsO4yGjcuE+d1XcSnlMpKNHA4K1WK4GX/o+YrC4hqugBK/gWnS8GfL8OFu5N8e5MmpfnjjxOxFvH5+HiwcGHzOHcZSimVVWngcPL97D28Ofojop/7GYybNVsqqI71OBn69PGnTZvKcRIPxneXoZRSWZUGDlto6HU69J5L9NMr4MSDsOEFuFYwxcfRxINKqexOA4ctJOQS3hFFCV/UGy75Ev82IAnz8nLjkUeKATpmoZTK3jRw2Pz8ClhbsDqSly7E3R3c3Nzw8nLH4dBV3kqpnEMDh+3s2TAcjqTreXu7MWrU47RqVQFAu6SUUjmOBg7b5s3J2xtKxI1WrSrcDBQaMJRSOY0GDlvNmsXjLffwgKgoa/GeiGiXlFIqx9PAYatQ4S569Kgaax/xdu0q06XLw+TN68XVqxHaJaWUUmjKkTiCgs6xefMpatYsnuB2rUoplR1pypFUqlDhLg0YSimVCN2XVCmlVIpo4FBKKZUiGjiUUkqliAYOpZRSKaKBQymlVIpky+m4IhIKHM3odrhYEeBsRjcik9BrYdHrYNHrcEtKr0VpY4xvUpWyZeDICURka3LmW+cEei0seh0seh1ucdW10K4qpZRSKaKBQymlVIpo4Mi6pmZ0AzIRvRYWvQ4WvQ63uORa6BiHUkqpFNE7DqWUUimigSMTEZEZInJGRPY6lRUWkVUicsj+XcguFxEZLyLBIrJbRKo5vaeNXf+QiLTJiM9yJ0TkXhFZIyJBIrJPRN62y3PUtRARbxHZLCK77OswxC4vIyKb7M80V0S87PJc9vNg+3U/p2O9Z5cfFJEmGfOJ7oyIuIvIDhFZZj/PqdchRET2iMhOEdlql6Xvvw1jjP5kkh/gP0A1YK9T2ShggP14APCp/bgZ8DMgQG1gk11eGDhs/y5kPy6U0Z8thdehOFDNfpwP+AuomNOuhf158tqPPYFN9uebB7SyyycDXe3H3YDJ9uNWwFz7cUVgF5ALKAP8Dbhn9OdLxfXoA8wBltnPc+p1CAGK3FaWrv82Mvwi6E+c/yn8bgscB4Hi9uPiwEH78RSg9e31gNbAFKfyWPWy4g+wGGiUk68FkBvYDtTCWtDlYZfXAQLtx4FAHfuxh11PgPeA95yOdbNeVvkBSgK/Ak8Cy+zPleOug93u+AJHuv7b0K6qzK+YMeYUgP27qF1eAjjmVO+4XZZQeZZkdzM8gvXXdo67Fnb3zE7gDLAK66/ki8aYKLuK82e6+Xnt1y8Bd5ENrgMwFugHOOznd5EzrwOAAX4RkW0i0skuS9d/G7qRU9Yl8ZSZRMqzHBHJC/wI9DLGXBaJ76NZVeMpyxbXwhgTDVQVkYLAIqBCfNXs39nyOojIM8AZY8w2EakfUxxP1Wx9HZzUNcacFJGiwCoROZBIXZdcC73jyPxOi0hxAPv3Gbv8OHCvU72SwMlEyrMUEfHEChqzjTEL7eIceS0AjDEXgbVY/dQFRSTmjz7nz3Tz89qvFwDOk/WvQ13gOREJAX7A6q4aS867DgAYY07av89g/TFRk3T+t6GBI/NbAsTMeGiD1d8fU/6GPWuiNnDJvkUNBBqLSCF7ZkVjuyzLEOvWYjoQZIz5wumlHHUtRMTXvtNARHyAhkAQsAZ4ya52+3WIuT4vAb8ZqwN7CdDKnm1UBigLbE6fT3HnjDHvGWNKGmP8sAa7fzPGvEoOuw4AIpJHRPLFPMb6f3ov6f1vI6MHevQn1gBXAHAKiMT6i6A9Vt/sr8Ah+3dhu64AX2H1ee8B/J2O0w4Itn/aZvTnSsV1qId127wb2Gn/NMtp1wJ4CNhhX4e9wCC7/D6sL7xgYD6Qyy73tp8H26/f53Ss9+3rcxB4KqM/2x1ck/rcmlWV466D/Zl32T/7gPft8nT9t6Erx5VSSqWIdlUppZRKEQ0cSimlUkQDh1JKqRTRwKGUUipFNHAopZRKEQ0cSqWAiGzI6DYold9G7C8AAAHBSURBVNF0Oq5SSqkU0TsOpVJARK7av+uLyFoRWSAiB0Rktr3iHRGpISIb7H00NotIPrH21vjG3kdhh4g8Ydd9U0R+EpGlInJERHqISB+7zkYRKWzXu19EVtqJ7daJSPmMuwoqp9Mkh0ql3iNAJawcP+uBuiKyGZgLtDTGbBGR/EAY8DaAMaaK/aX/i4g8aB+nsn0sb6xVvP2NMY+IyBjgDay8TFOBLsaYQyJSC5iIlbNJqXSngUOp1NtsjDkOYKc+98NK4X3KGLMFwBhz2X69HjDBLjsgIkeBmMCxxhhzBbgiIpeApXb5HuAhO0vwo8B8pwzBuVz82ZRKkAYOpVIv3OlxNNa/JyH+9NQJ5oS/7TgOp+cO+5huWHtPVE19U5VKOzrGoVTaOgDcIyI1AOzxDQ/gD+BVu+xBoBRWor0k2XctR0TkZfv9IiIPu6LxSiWHBg6l0pAxJgJoCUwQkV1Yu/Z5Y41JuIvIHqwxkDeNMeEJHymOV4H29jH3Ac3TtuVKJZ9Ox1VKKZUiesehlFIqRTRwKKWUShENHEoppVJEA4dSSqkU0cChlFIqRTRwKKWUShENHEoppVJEA4dSSqkU+T+l+LTHAjxosQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "engel.plot.scatter(x='income', y='foodexp', color='darkblue')\n", "plt.plot(X_pred[:, 1], y_pred, color='darkgreen')\n", "plt.plot(X_pred[:, 1], X_pred @ betas_lr, color='red')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }