{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Models\n",
    "\n",
    "There's a few linear models out there that we felt were generally useful. This\n",
    "document will highlight some of them. \n",
    "\n",
    "## LOWESS\n",
    "\n",
    "Lowess stands for LOcally WEighted regreSSion and has historically \n",
    "been used for smoothing but you can also use it\n",
    "for machine learning where you're interested in interpolating. Here's \n",
    "a demo; "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsMAAAEICAYAAAC6S/moAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeXxU5fXH8c/JBmFLQNaEJSiLsqgg4oIrLrhVqdqqrVZbra2tdrO02NpqrVaqrdZW+7O2drG21VYtdccFlcUVBARBENkDSFjCmpDt+f1xZ2ASZiYzyczczMz3/XrNK7PcufckN5Azzz3Pecw5h4iIiIhINsrxOwAREREREb8oGRYRERGRrKVkWERERESylpJhEREREclaSoZFREREJGspGRYRERGRrKVkWEQygpldZ2afmtkuMzuomW2vMrNZrTzeiWa2NMrrZWbmzCyvNcdp68zsFDNb14r3P2hmP0lkTIH93mpmjyZ6vyKSeZQMi0gjZrbKzKoCSeU2M3vOzPr5HVc0ZpYP3AOc6Zzr5JzbkuxjOudmOueGhsSwysxOT/Zxm0pEYp8q4WJ1zn3dOfdzv2ISEVEyLCLhfMY51wnoA3wK/M7PYGIYXe0FtAc+TEE4acfMcv2OIZ1l+ui+SLZTMiwiETnnqoEngGHB58ysnZn9yszWBMoSHjSzwkj7MLOvmtkSM9tpZovNbHTg+clm9knI858Nec9VZjbbzO41sy3ArZGOa2ZDgGC5QqWZTQ9XomBmr5vZNc19z2b2NzO7MXC/NLCfbwYeH2JmW80sJ7Q8wMz+DvQHngmMqP8gZJdfDMS82cx+HOW4EX+uZva8mf06ZNvHzOzPZnYY8CBwXOC4lYHX/2pm/xd4327gVDM718zmmdkOM1trZreG7C/487rWzNab2QYz+36T2H4TeG194H67CN9H2PPaTKy3h7z/q2a2PPBzftrMSkJec2b2dTP72MwqzewBM7Mop7PAzB4JxPKhmY0J2VeJmT1pZhVmttLMvhXy2q1m9oSZPWpmO4CrzKzIzB4O/GzKzex2fcgQyQxKhkUkIjPrAFwCvB3y9BRgCHAkMAgoBX4a4f2fA24FvgR0Ac4HgiUMnwAnAkXAz4BHzaxPyNuPAVbgjfreEem4zrllwPDAe4qdc+Nb/A173gBOCdw/ORDDSSGPZzrnGkLf4Jy7AlhDYETdOXdXyMsnAEOB04CfBpLCcKL9XL8CXGFm483si8BY4NvOuSXA14G3AsctDtnfF/B+bp2BWcBuvPNQDJwLXGdmE5vEcCowGDgT+KHtL/v4MXBsILYjAse/OcL3Efa8NhMrAGY2HrgT+DzeVYnVwGNNNjsPOBo4PLDdhAhxgPf79ljge34auD9wnBzgGWAB3s/5NOA7Zha6rwvwPggWA/8A/grU4Z2bUYGfUbMfrkQkDTjndNNNN9323YBVwC6gEqgF1gMjA68ZXlJ1SMj2xwErI+xrGl7SFstx5wMXBO5fBawJeS3qcYEywAF54R4HnnsduCZk/7MixHEIsA1vsOBB4GvAusBrfwO+F7h/SvD5kJ/b6SGPgzH0DXnuXeDSMMds9ucKXASsBTYDJ4Q8f8D3gpe4PdLMz/s3wL1NYj005PW7gIcD9z8Bzgl5bQKwKtzPIYbzGi7W2wP3HwbuCnmtU+B3sCzw2DX53v8NTI5w3FuBV0IeDwOqAvePCf39Cjx3E/CXkPfOCHmtF7AXKAx57jLgtWT/e9RNN92Sf1MdlIiEM9E590rgMvAFwBtmNgxoADoAc0OuThsQ6XJxP7xE6gBm9iXge3iJGHiJT/eQTdaG3O8R53FbzDn3SaC04Ei8Ec6fA1eb2VC8keHfxrnLjSH39+B9n03F8v09g1e7vdQ5F8uEudCfH2Z2DN7o8wigAGgH/CfKe1YDIwP3SwKPQ18rIYwYzms0JcD7wQfOuV3mlcmU4n3YgNh+nkFNt20fKJ0ZAJQESzUCcoGZIY9DfxYDgHxgQ8j5yWmyjYikKZVJiEhEzrl659xTQD3e5f7NQBUw3DlXHLgVOW+yXThr8UZaGzGzAcAfgeuBg5x3yXwRXgK47/Ah9+M97u7A1w4hz/WO+s029gZwMVDgnCsPPL4S6Io30hmOi/B8LGL5/u4AlgB9zOyyGI7b9Pl/4pUK9HPOFeGNejettw3tGtIf76oAga8DIry2TwzntbmfUaPjmFlH4CCgvJn3xWst3qh7ccits3PunJBtXJPt9wLdQ7bv4pwbjoikPSXDIhKReS7ASwKXOK9W9o/AvWbWM7BNaZNay1B/Ar5vZkcF9jUokDB1xEs2KgL7+DLeiGVY8R7XOVeBl0Bdbma5ZvYVwiTlUbyBl9DNCDx+PfB4lnOuPsJ7PgUOjuMYofFG/f7M7CTgy3g1v1cCvzOz0pDj9jWzgmYO0xnY6pyrNrOxeDXFTf3EzDqY2fDA8R4PPP8v4GYz62Fm3fFqmcP18G3uvDYX67+AL5vZkYEJer8A3nHOrWrme4vXu8BOM/uheZMwc81shJkdHW5j59wG4CXg12bWxbwJlIeY2ckJjktEfKBkWETCecbMdgE78EYkr3TOBduW/RBYDrwdmGn/Ct4EsQM45/4TeP8/gZ3AVKCbc24x8GvgLbwEaSQwu5mYYj5uwFeBSXgT9oYDbzaz/1Bv4CWPwWR4Ft4o84yI7/Amft0c6HLw/SjbRRL2+zOzLsAjwPXOuXLn3Ey82tq/BDopTMdrKbfRzDZH2f83gNvMbCdeMvvvMNu8EYjhVeBXzrmXAs/fDswBPgAW4pUy3N70zTGc16ixOudeAX4CPAlswPsAc2mU76lFAh9ozsMrhVmJNzL/J7xJf5F8Ca+8ZDFeTfkTeJP8gguw7Ep0nCKSGuZca67siYhIujOzMrykMN85V+dvNCIiqaWRYRERERHJWkqGRURERCRrqUxCRERERLKWRoZFREREJGv5tuhG9+7dXVlZmV+HFxEREZEsMXfu3M3OuR7hXvMtGS4rK2POnDl+HV5EREREsoSZrY70msokRERERCRrKRkWERERkaylZFhEREREspaSYRERERHJWkqGRURERCRr+dZNQkREJBNMnVfO3dOWsr6yipLiQiZNGMrEUaV+hyUiMVIyLCIi0kJT55Vz01MLqaqtB6C8soqbnloIoIRYJE2oTEJERKSF7p62dF8iHFRVW8/d05b6FJGIxEvJsIiISAutr6yK63kRaXtUJiHSBrWmBtGv+kXVTUo2KikupDxM4ltSXOhDNCLSEhoZFmljgjWI5ZVVOLwaxO8+Pp+yyc8xbsp0ps4rj+u9Nz21MOp7khVzKo4r4rdJE4ZSmJ/b6LnC/FwmTRjqU0QiEi8lwyIpMHVeOeOmTGdgDAltuBpEF/jaXJLpV/2i6iYlW00cVcqdF46ktLgQA0qLC7nzwpG6KiKSRlQmIZJk8c42b67WMJhkxvPeZNcvqm5SstnEUaVKfkXSmJJhkSSLNmoa7g9opBrEUJGSTL/qF5s7bjz1xKo9FhGRVFIyLJJk8Y6aTpowtNFIcjiRkttw701W/WJo0lpUmE9+rlFb7/a9HjxuPCPj6tkq6SJRH9r04U/Ef0qGRZIs3tHa4B/Cu6ctpbyyCmN/zTAcmNw2/WN60VGlvPZRRcL/uDZNfnfX1O1LfiurasnPMbp2yKdyT22j446bMj3mkfF4R9FFUiXa739LP7Tpw59I26BkWCTJWjJaG1qDGG3kKNwf0yfnlid8Ak/T41RW1R6wTW2Do0NBHvN+emaj5+MZGY+0bXllFQMnP6eRM/FFLL//LfnQpg9/Im1Ds8mwmfUDHgF64Q1QPeScu6/JNqcA/wNWBp56yjl3W2JDFUlPoSO9LRmtjTY5J9F/TCMl3uGOE04wmQ3dT44Z9c4dsG24kfFo9dKhLdtAI2eSOvH+/jcn+O8j0u+6Jp6KpFYsI8N1wI3OuffNrDMw18xeds4tbrLdTOfceYkPUST9NU1og63WWlvKEG89cryjzMHEM9Y/ziXFhQfsJ1wiHGlkPJZ6aY2cSSqE/ls58Dc4vFgmqjb999HS/YhI4jSbDDvnNgAbAvd3mtkSoBRomgyLSAwSWScYTz1yc8eNNsocS4eLYIIbaRQt14wG56Im/01H0SMlIRo5k2SKJWFtKtaJqs2NMmvBDpHUi6tm2MzKgFHAO2FePs7MFgDrge875z4M8/5rgWsB+vfvH2+sIimXjJneiSxtaK4eublyhdDjRhtlvveSIw84Tn6O0al93gET5r77+Pyw+2lwjpVTzt0XV6SR8dBR9HFTpmupW0m5WMoiIv3+hxPrKHOpauJFfBFzMmxmnYAnge8453Y0efl9YIBzbpeZnQNMBQY33Ydz7iHgIYAxY8bEeuVJxBfJmumdyAUqotUjN42/ne2he34lB+Vtp73tpX1ODe2sBhrgul+8wllFOeyuL2BPQ3v2NBSypa4LW+uL6FnUJa6651h6Djf9uU76zwJ+9syHByQWqWwVJxIU7d+iQVwfjGMdZS4tLmT25PHxhioiCRBTMmxm+XiJ8D+cc081fT00OXbOPW9mvzez7s65zYkLVSS1kjXTO9ELY+wbSa3dAds/gh0vw7wlFM19i0f6b6FH3jZ65G2jY251i/ZfWdeJZX/tzdiDBjN7wnAoGg7FI6FLt7DbN5fAhvu51jY4tu3xZuiH+9ChPqySSpH+jbYkYY1llDnS1Rz9voukRizdJAx4GFjinLsnwja9gU+dc87MxgI5wJaERiqSYslaYjgho51Vn8Lmt7zbtnmwYwnsWbf/9ZwC+uX05NPabizYM4SKumIq6rpSUduVrfVFVDW0o6qhHTUuH4cBkEsD7XP20il3L4W2h4PydtAtr5IeedvoW7CJss2LqN81nVy31zuG5UDxEdBjHPQ4AXqdCu17NpvAxvLzC/3QoaVuJdUSeUUinlFm9R0W8UcsI8PjgCuAhWYWLAb8EdAfwDn3IHAxcJ2Z1QFVwKXOhZlCLpJGkrW0cdyjnc7B9kWwaYaX/Fa8CbsDXQxz8qFoJPQ8FYqGQdFh0GUYdBrIlXfNCBt/cCJbpH+gwT/Q5VsPfG+/4gJmfmMAVH4A2+Z78XzyZ1h2v7dBtzFQcjYTB0xk4g9PBbMD9hHLZDzQJDnxT2uvSMTSWjDcKLP6Dov4w/zKWceMGePmzJnjy7FFYhGu1q8wPzfhC1qEVb0ZNr4MG6bBxpegakMggD7Q/Tjofrz3tdtoyG3fovgjTU4rLS6MONHHYN8kuH0aar3EeMM0WP8CbHkbXAN0GgQDLoGyy6Ho0KhxhaMaSklHsfx+R/p/ZODk52L/dycicTGzuc65MeFe0wp0IhGktF7VOS+hXPdfL6HcOhdwUNANep8OfSZAr/HQcQBT56/n7v8tZX3lNkqK34y5TVnT+KNdCo60IEDYUfGcfDjoaO824mbYuwXW/hfWPA6L74QP7/DKKA75KvS/+IC4mi5tGxpHOKqplLasNa0Fk3U1SkSi08iwiF9cg1dmsPYp77Z7lVeH2/046D3BS4C7HQU5ufvekujR6kiJZcKOU/UprHwEPvkT7FwG7Q6Cwd+EId+E9j2bjSNcvL6N1ovEoDWju/r9FkmeaCPDSoZFUsk52PwmrPqHN3pavdEbWe19BvS7EErPh/Y9Ir49WmlDoksKEjoC6xxsegM+uhfKn/ZKOw7+Cgz/EXSIfZ+p/P4l8yXjKkNrf0d15UMkOVQmIeK3ncth5d9h1aOwawXkFkLJOdDvIu9rQVFMu0lWh4twEtrFwQx6neLdtn8EH/0alj8EnzwMg6+DYZOhsFezu0nl9y+ZLVmdG1rbiULdU0RST8mwSLLs3QKrH/eS4C1vA+a1HxvxEy8Jzu+8b9NYR4Myoqaw6FA45o/eqPCin8Oy33lJ8Ygfw9BvR5wQCBny/UubkKzODeqNLZJ+lAyLxCimhLV+L6x/zkuA1z/ndVooGg5HToGyL0KHvmH3G+sIVUatyNZpIBz7Z29UeN4kmD8ZPv4DjL4H+k0M+5aM+v7FV8m8yqDRXZH0omRYJAZRE9YjS7w64JV/hzX/hppt0L4XDL4eBl4BXY8M2283KJ4RqowcdeoyBE7+H2x8BeZ+B2Z+FvpOhDG/O+DDQ0Z+/+ILXWUQkSAlwyIxCJew9rI1bH37n7B29v464L6f9RLg3qdDTmz/vOIdocrYUafep8PZ8+Cj38DCW+DZYd6I+uCve102AuL5/jUZSSLRVQYRCVIyLBKDYGJanLuD84pncmHxdEZ3XEqDM+g4Hkb81OsGEVIHHCuNUIXIyYdhk6D/RfDu12HON73ey8f8GTr2i2tXWtpWotFVBhEJUjIs0pz6vXyhz3ucXDCNUzrPoSCnjo+qBnDnhqt4p+Espn7xkrBvi3VUUiNUYXQ6GE6dBp/8Ed7/Hjw/EsbcDwMvj3kXWtpWmpMOV1l0dUMk+ZQMi4TjGmDTTK8f8Jr/cEePSirquvK3Lefx323jWVw9kML8PO68cGTYt8czKqkRqgjMYNC10Os0ePsqeOsK+PRVLynO69js29WGTdKdrm6IpIaSYZEg56ByAaz6J6z+F+xZ5yVdfSdC2RW8ufFQ/vrSJ6yvrqK0mYQ13lHJdBih8k3nQ+C012HRbV4rti3vwgn/gaJhUd+m8hNJd7q6IZIaSoZFdq2C1f/0RoG3L6bO5fLGztHMqLmSMeOu5DNHDQbgghK4YPSAmHapUckEy8mFw38GPU+EN78ILx4Nx/0V+n8u4ltUfiLpTv+PiKSGkmHJTnvWecshr37Ma4sGbCk8mvs3fJOpW49nW723Ity/p35CfU6HsKMw0Wr5NCqZJL1Ph7Pnw8yLYdbnvYU7Rt7mJctNqPxE0p3+HxFJDXPO+XLgMWPGuDlz5vhybMlSOz+BtU/B2idhyzvec0UjoOwLMOAyxt2/IuwfntLiQmZPHt/ouaa1fAD5OUan9nlU7qmlqDCf3TV11Nbv//dVmJ/LnReOVDKWCPV7Yc4N3gS7knNh3D8hv0vcu9HkJAlqi78L4f6f0f8jIi1jZnOdc2PCvaaRYclczsH2xfsT4MoF3vPdjoIj7vCWRO6y/5L5+soPw+4m3CXJcLV8tQ2ObXtqAaisqiU/x+jaIZ/KPbVt5o9rxshtB2P/AN1GeUnxyyfCyc/G1X5Nk5MkqK3+LujqhkhqKBmWzNJQB5vf8pZCXjcVdiz1nu9+PIz6tdcLuFNZ2LfGc0kylpq92gZHh4I85v30zLi+BYmRGQy+DjoNglkXw0vHeglxt1ExvV2TkySoLf8uaHKtSPIpGZb0V70ZNrzoJcAbpnnLIVse9DwJhnzL6wbRoaTZ3cQz4SpS4tyUJrqkQJ8z4IxZ8Pq58MqJcOJT0Kf5DyCanCRB+l0QyW5KhiX91NfAlrdhw8uw8ZVA/a+D9j2h9HwoPRd6nwkFRXHtNp5LkuES53A00SVFikfCmW/D62fDG+fBcX+HAeEXQwlq7kpAW6whleTQRDWR7KZkWNo+1wCVi7wFFza8DJvegPo9YDnQbay3FHLpuV4tsOW06lCxXpJsmjhHmjCnNl4p1KEETn8D3jgfZl8Ge7fAkG9E3DzalYC2WkMqyaE2fCLZTd0kpO2p3wtb53grwFXMgorZUFvpvdZ5CPQ+w2ux1esUKCg+4O1+jehpJLFteHrucrrN+xIndHiLB7Z9ldLjb4l4HiKds3FTpsfcWUQyQ7r8+02XOEXammjdJJQMi/9qtnu9foPJ75Z3oWGv91qXQ6HHCd6t16nQsX/UXTXXikh/SDJb8PzX1NZwb/9fc37xDH696csccvLP4zrPAyc/R7j/GQ1YOeXchMUraay+Bio/gOqNUF0BtTu8qxMdB3qrJhZ0Tfgh1WpNpOXUWk3ajrrdsHWeN/IbvAU7PlgedBsNQ64PJMDjoH2PuHYfbVY4oEvfGW7/+c/lu2tupMEZN/b8C394y2DUwzHvRzWkElZNJax+HDa8ABtfhbpdkbftfhwMuMxbJbGwd0IO35a7XoikMyXDkjy1u2D7ItgSmvgu8WqAAQpL4aAxMOCLXuLb/RjI69iqQ0abFa4/JJkv9PzXk8v31n6PBnL4Wtc/w8IyGPmTmPYTroY0P8fYU1PHwMnP6apCmmrxlaGqT2HpvbDs91C3EzoOgLLLvXKtDv28D+15naGqHHat9P7fW/MfmPsteP87cPBX4PCftzopVtcLkeRQMiyt11ALO5ZB5ULvj0Dw664V+7dp3xO6HQ39LvYS4G5HQWGfhIcSbURPf0gyX9Pz30Au31/7HQoL8jh74U+BBhh5S7P7iTRBMrioiq4qpJ8WTYqs3wsf/gKW3AUNNdDvczDsB9B1lNfnuqn23aHrEdBvIoy42Vv05+M/wMe/95Z+H/4jOPR73qIxLaArFiLJoZphiV3dbtj5sZf47lwG25fA9oWw4yMvIQawXG9Vt6IRXrurohFe8ltYGv6PR4JFq6m7e9rSsH9Ics1ocE6jfRkg4vn/7DAm1twCK/7qdR8ZeWtcv4+aUJf+4j6Hm9+Gd672EtoBX/B+Z7oMbtnBdyyD+T+Adf/zBgVOfDKu1RKDVDMs0nKqGZbY1VXB7tWwa/n+pHfnMu9+VXnjbTv09xLeknOgaCQUj/AmvLVw1CMRmusVHK43cH3gA6FG+9Jf1PPvHgZyYNFt3u/o8B/FvF9dVUh/MZ9D1wALb4VFt0OHvnDK81BydusO3mUInDQV1v4X3roSXhwN4x6D3qfFtRstzyySHEqGs019dSDZXQW7A7ddK/ffr/608fYFXaHzUO8/7S5DvdZmnYdA50GQ1yHl4cciUq/gpn9Icsz2JcJBqiFOfxF7RVsOHPNH73L3gh9DwUEw+Gsx7VOXp9NfTOewdge8eTmUPwMDr4Qxv4X8LokLot9noWgYzLwQXjsTxvw+5t/BIC3PLJJ4SoYzSV0VVK33RnD3BL5WrYc95bBnjZf0Vm9s/J6cfG+Et2MZlJ7nfe04EDoN9JLfdgf58Z0kTegfkoGTnwu7jUb70ktck6IsB479s7dk93vXeR/2Bny+2WNoUYb01+w53LncW7Bl5zI46ncw5JvJKe3qMhTOfAdmXwrvfR1ogMHXtWhXahUpkhjNJsNm1g94BOgFOOAh59x9TbYx4D7gHGAPcJVz7v3Eh5uFXIPXzqd6E+zd5H2t3hQ+6a3ZduD7cztAYYk3+7n03ECyG7h1KoP2fSAnN7XfUxuh0b70F++kqGDysHX7l3ls8GpGzr6cnIJi6HNm1OPo8nT6i3oOKxfC9DPA1cH4l72e5smU38mrG571OXjvG97/80O+GdcutEqiSOI0O4HOzPoAfZxz75tZZ2AuMNE5tzhkm3OAG/CS4WOA+5xzx0Tbb9ZOoHMN3qW4vVtgb8X+5HZfolvROOndWwGu/sD9WC607w0dSr1kt7B0//3Q5/K7pGTiWlsWafREk1HSXzyTopqe7y45u/j3oJsYVLiRvDOmQ/djUxKztDFb5sBrEyC3PYx/FYoOTd2x62u8hLj8aRj7Rxh0Tcxv1aROkfi0agKdc24DsCFwf6eZLQFKgcUhm10APOK8zPptMys2sz6B92am+r3eSGzNNqjZCnu37r8f6XFtYPtgn92m8jp7Lcja9/RGbg8a691v13P/8+167P+apSO68Yhl9ESjfekrnoltTftM72joxBUrbuO/g39I39fPgdNnQvHwpMUqbVDFbHjtbGjXHU57BTodnNrj5xbACf+BGRd4ZTudDobesSWymtQpkjhx1QybWRkwCninyUulwNqQx+sCzzVKhs3sWuBagP79oy+rmzSuAWp3Qu12b4Q2+LVm+4HP1W6P/HpDTZSDGBQUQ0G3wK2r959c6ON23aB9r8bJbZ4uzydacwttaDJKeoun1CVcklBR15XLlt/GzNE/8SY0nTHbKx+SzLd1Hrx+jtfv/LRXvc4Rfsgt8DpLvHw8zLrYqyeOoYWbyrxEEifmZNjMOgFPAt9xzu1oycGccw8BD4FXJtGSfbTKrlXw9MDmt7McyOvilRgUFHlfC3t7Ex/yu0B+4LnQhLddSKKbX6RR2zZCoyeZLZ6JbZGSh4aOA+HUafDySV5ydOab3r9tyVw7lnmlEfnFMP4V/xLhoIIiOPkZmHYMvHEeTHjb+1sShSZ1iiROTMmwmeXjJcL/cM49FWaTciC0g3jfwHNtS7vuXsP90IQ2v2h/wht8Lq9T1tfZZgqNnmS2eEpdoiYPxaVw0lNegjTzIjjlBW/ETjLP7rXeZDnwJsu1YPGLpOh0MJz4FEw/Dd78Epz8dNS/QyrzEkmcWCbQGfA3YKtz7jsRtjkXuJ79E+h+65wbG22/WTuBTlJKk+QkVLOtqFb8Dd6+Cg7+MhzzsD4UZ5raHfDSOK/V5GmvQbfRfkd0oKW/hbnf9tq7Db3e72hEMkZrV6AbB1wBLDSz+YHnfgT0B3DOPQg8j5cIL8drrfbl1gYtkggaPZFQzdaIH3wl7FrhrVLX6RAY8ePUBSfJ1VAHsy7xlo8/9cW2mQgDDLkBNkyDed+HnidB18P9jkgk4zU7MpwsGhkWkTbJOXjrClj1Dzj+n1B2md8RSWs5B3Ouh49/D2MfgkFf9Tui6Ko3wfOHe4seTXivza72KZJOoo0M56Q6GBGRNs3MK5HoeZJXMrFplt8RSWstu99LhA/7fttPhMHrMHTcI7B9Mcz/od/RiGQ8JcOS9qbOK2fclOkMnPwc46ZMZ+q8tjd3U9JMbjs48b9ev+8ZF8COj/2OSFrq09fg/e9C6flwxBS/o4ldnzNhyLdg2QNeP2QRSRolw5LWghPkyiurcOxfVEMJsbRau25wyvNeq8UZn/GWRZf0snuNVyfceTAc//f0a3l5xB3QsT+8cw3UV/sdjUjGUjIsaS3aohoirdb5EDjxSdj5Ccy+DBrCLI0ubUboVaJTp7zItmnne0nkSVO9tpnpJr8THP2gN+lv0R1+RyOSseJagU6krdGiGpJI4VuvnQRH/x7evRbm/wBG/9rvMCVE8JyVV1ZhQHBK+Nc73UPX6gW83fdhju2SxgtRlJwFZZfD4ikw4PNQPDLips22DhSRsJQMS1rTohqSKE17UiKMgjMAACAASURBVAdLbgAmjvoqVC6Ej+7xkpGDr1Li0QY0PWfBRHhi8Wtc0u1lfvfpJTy2pozZJ/kXY0KMvhc2vAjvfg3OmOWV7jQR/fdXv5ci0ahMQtLapAlDKcxvXAeoJUmlJZotuRl9D/Q+Hd79GjNm/0+16m1AuHM2sKCcO0of4J1dw/nNp1/IjKtE7bvDkXfB5re8ln9hqGRMpOWUDEtamziqlDsvHElpcSEGlBYXanU5aZFmS25y8mDc49ChPyNWXEVXNjTarqq2nhv/vUBdTVKo6TkrsFruH/BLalw+3147iXpyM+cq0cFXQrejYd4PoHbnAS+rZEyk5ZQMS9qbOKqU2ZPHs3LKucyePF6JsLRIpKSp0fPtusHJz5BHDX8su51CazzDv945jRSnUNNzdlOfPzO8cAXfX/sdNtZ2Jz/H2FNTlxkfUCwHxvwOqjfCotsPeDmm318RCUvJsIgIcZTcFB3Kz7bczGHtV/Krfveyv1K1MV2iTr7Qc3Zml7f4cvdneLjiAl7deQzFhflgsG1PbeZ8QOl+DBx8FSy994De1yoZE2k5JcMiIsRXcnPiyV/iV5u+wrnFs/lmz39H3KcuUSdX8JyNPmgHd/W9j4/2DqHHCfewasq5dGyXR2194w8qGfEB5Yg7Iac9vP+9Rk+rZEyk5dRNQkQkYOKo0piSh4mjSpnqbmXa3FXc2OtRllUP5OUdYw/YLtIlanWiSJyJR/RkYsX9UGkUn/0ch3Y+GMjgGtrC3jDiZm+Z5k9fh16n7Hsp1t9fEWlMI8MiIi0wcXRfJnz5GXIOGs3vD76H4R0bX36PdIlaqyYm2MJbYfObMPYP0HnQvqczuoZ2yA3QoS/MmwSuwe9oRNKekmERkZbKK4QT/0t+QQceP2wKh3ata/YStVpgJdDGV+DDO+GQq6HsskYvZXQNbV4hHH47bJ0Da/7jdzQiaU9lEiIirdGxH5zwJJ2mj+fFsf8HJz8HObkRN8/Yy/epVvUpvHk5FB0GR/32gJeDH0Qythyl7HL46New4EfQ97OQW+B3RCJpS8mwiEhr9TwBxjzgLdm84CYYdVfETbVqYgK4BnjrCqjdDuNfgbwOYTfL6BranFxvIY7Xz4blD8LQb/kdkUjaUpmEiEgiDPoqDP4GLLkbVv0z4mYZffk+BabOK+fB/7sWNr7MlE3fYOrKrn6H5J8+E6DXeK/vcO2uA16eOq+ccVOmZ0afZZEkUjIsIpIoR/0Gep4E71wNW+eG3UQtsFpu6rxyHn/hca4p+gvPVJ7Ig+WnZvfkQzOvdnhvBXz8QKOXNFFTJHbmXPiG8ck2ZswYN2fOHF+OLSKSNNUV8OIYoAEmzIHCXn5HlDHOumsqf+p5LfUul/M+vo+dDR0B7wPF7MnjfY7OR6+dDVvfg/NXQn5nAMZNmR62HCfrf1aStcxsrnNuTLjXNDIsIpJI7XvASVNh7xaYdRHU1/gdUWZwju92/iU987Zxw5of7EuEQZMPGfkz7/dt2e/2PaWJmiKxUzIsIpJo3UbBsX+BitkwVxObEmLZA0woeptfbrySD6qGNHop6ycfdh8LJefBkl9BzXYgw/ssiySYkmERkWQYcAkMuwmW/wE+ftDvaNLb1nkw70Y2djqNf26/qNFLmnwYcPitULMNlnpt5jRRUyR2SoZFRBIsOIv/kMeO5c2qY2h47wbYNMPvsNJT7U6YfQm060HvMx/jzgsP1+TDcLodBX0v8HoP11RqoqZIHDSBTkQkgYKz+IOrzHXO2c3/Bt9IaYcq2p03FzoO8DnCNPPml2D1P+C017xOHRLZtvnwwigY8VM4/Gd+RyPSpmgCnYhIijRdbnlnQ0euWXUztbXVMOOzULfHx+jSzIq/waq/w4hblAjHouuR0O8iWPob2LvV72hE0oaSYRGRBAo3W3/F3r7csPr73sjdO1eDT1fk0sr2xfDeN6DnKTD8x35Hkz5G3gK1O+Cje/yORCRtKBkWEUmgSLP1l+WeBEfeCasfgyWRl2sWvGRu5oVez9zj/+EtPSyxKR4J/T8PS++D6s1+RyOSFpQMi4gkULhZ/Pk5xp6aOgY+MpyXd5+Km38TlD/vU4RtnHPw9ldg53IY9zh0KPE7ovQz8hao2w0f/crvSETSgpJhEZEEajqLv7gwHwy27anFYdyw4pssqT6Y2pmXQuVCv8Ntez66B9Y+CUf+Enqd7Hc06aloGAy4FJbdr9FhkRg0mwyb2Z/NbJOZLYrw+ilmtt3M5gduP018mCIi6WPiqFJmTx7Pyinn0rFdHrX1+2uEq117rl55M5U1BfD6ObCn3MdI25hNM2iY90Om7zmJgX8byrgp05k6Tz+fFhlxszc6vPQ+vyMRafNiGRn+K3BWM9vMdM4dGbjd1vqwREQyQ7gJdRtqe3DlilugphJeP9frpZvt9qyn+rWLWbW3DzesuB6HUV5ZxU1PLVRC3BJFw7zOEst+6/2eiUhEzSbDzrkZgHq0iIi0QKQJddvbD4cTnoDti2DW56ChNsWRtSENtTD78zTU7uTaVT9id0OHfS9V1dZz97SlPgaXxkbc7E1GXHa/35GItGmJqhk+zswWmNkLZjY80kZmdq2ZzTGzORUVFQk6tIhI2xV1WdySCXD0g7BhmtdGLFtbrs37AVTM5gdrv8Xyvf0PeDnc6LrEoOuRUHIefHQv1O7yOxqRNisRyfD7wADn3BHA74CpkTZ0zj3knBvjnBvTo0ePBBxaRKRta3ZZ3EHXwPCb4ZM/wYe/8DVWX6z8u7dIxNBvM88mhN0kx4yBk59TDXFLjPgx1GyF5Q/6HYlImxXTcsxmVgY865wbEcO2q4AxzrmoU1i1HLOISIBz8NaXYNWjcNzfYeDlfkeUGptmwfTToMcJcOqLTF2wqdFS1uEU5uc2/jAhzZt+hte55PyVkBe+bEck0yV1OWYz621mFrg/NrDPLa3dr4hI1jCDYx6GXqfCO1+B9S/6HVHy7VoBMz8LHcvgxCcgJ/+AUfRc709LI6ohboHhN0P1p97VBxE5QCyt1f4FvAUMNbN1Zna1mX3dzL4e2ORiYJGZLQB+C1zqYhluFhGR/XIL4MSnoGi4t/rappl+R5Q8NZXw+nng6uHkZ6Gg676XQtvSNUT4U6Ia4jj1Ohl6nOitfFi/1+9oRNqcWLpJXOac6+Ocy3fO9XXOPeyce9A592Dg9fudc8Odc0c45451zr2Z/LBFRDJQQTGcOg069vdarm2d63dEiVdfzebnz6J2+8dctmQS436/NmIdcKROHJGelyhG3Ax71sHKR/yORKTN0Qp0IiJtSfueMP4VaNcNXpsA2z7wO6LEaahn/fOfpfued/jemu/y1u7Do/YSjtqJQ+LT+wy2tT+C8jdvYdDkpzUZUSSEkmERkbamQ18vIc5pD9PHw7YFfkfUes7BnG9QsvNFbi2/lme2719qOVIdcLOdOCRmU+ev58cfT6Q0bwOfKX5DC5qIhIipm0QyqJuEiEgzdi6HV0+Fuj1w2qte39h05BzMnwxL7uKBTZ/j7o1XHrCJASunnJv62LLEuCnTKa/cwwuDb6DA6jhj2QM0kEtpcSGzJ4/3OzyRpEtqNwkREUmSzoPgtNchryO8Oh42v+N3RPELSYQZ/A3+Wf218JuBLt0nkTfp0Pjdpks5pP06zi56M+R5keymZFhEpC3rfAic/obXceHV8bD+Bb8jip1zsOCmfYkwY+5n0oRDD6gDDtKl++QJTjp8cftxLK/uy/U9H8do0GREEZQMi4i0fZ0GwhmzocsQeON8b9W2tq6hHuZ+Gxb/EgZfB2PuB7NGdcDhqI9wcgQnIzaQywObPs9hhas4u+tcTUYUQcmwiIjvps4rZ9yU6dGXHC7s7Y0Q9zzJW63ug1vBNaQ81pjUV8PsS2HZ72Dod/clwkHBXsIHLqnh0aX7xAv9EPJM5cmU1/XhjsFTmXhkid+hifguz+8ARESy2dR55Y2WIA6WCgAHdk3I7wKnPA/vfg0W/QwqF8Bxj0B+51SHHdnerTBjIlTMhNH3wKHfjbhpSXEh5WESX126T46Jo0r3/04t/xTevRY2vgx9zvQ3MBGfaWRYRMRHd09bui8RDopaKpDbDo79C4z+DZQ/Ay8dB9s/SkGkMdj6Prx4FGx5B8Y9FjURBvUR9tXAK70Wfot+7tV2i2QxJcMiIj6KVBIQtVTADA79trdaXfVGeHE0fPx//iY1nzwMLx0Prg5OnwEDLmn2Leoj7KPcAjjsh1AxCzbN8DsaEV+pz7CIiI+8/q/hE9/S4kImTRgaPTncsx7e/jJsfAlKzoOxf4AOKawDrd4Ec74Fax6H3mfA8f+A9j1Sd3xpuboqeHogFI+E8S/7HY1IUqnPsIhIGxWuVCCovLKKSf9ZwKjbXoo8ua5DCZz6Ahx1n1f/+exQWHw31NckN3DnYMUj8OxhsO6/cPjtcMoLSoTTSV4hHPZ92PgKbH7b72hEfKNkWETER821GqttcGzbU4sjSh9ey4Gh34JzP4Rep8L8H8DzI2HNf7wWZ4n26Wvwyknw9pVQdBicPR9G/Bhywif10oYN+jq0OwgW3eF3JCK+UTIsIuKz5lqNhYo6ua7zIXDy017HCcuBWZ+HZw+Fj//gLencGg11sP5FeOVUb/GPXSvg6Ae9+uCiw1q3b/FPfiev/d36Z2HrPL+jEfGFkmERkTYi1pZizfbhLTkbzlkEJzwBBcXw3tfhyR4w61JY+xTUbI8toPq9UPEWzJsE/+sPr58NO5Z4nSzO/wQGf81LuiW9Dbke8ovgQ40OS3ZSn2ERkTZi0oShjXoORxJT0pyTC/0vgn4Xej1/V/0L1j7pTXQD6DwYuh0FHQdCQZGXDDXUehPi9lbA9kWwZQ407AXLg5JzYOCXoPQ8r71bHKbOK+fuaUtZX1lFSSyTAiW1CopgyA3w4e1Q+SEUD/c7IpGUUjcJEZE2JDRxLCrMZ3dNHbX1+/+fLszPbXn7sYY6r43W5jdh61yvL3DVeq8dWpDlQMFB0OkQ6DEOehwPPU6C9t1b/P00TfBb9T1IclRvhqfLoO9EOP5Rv6MRSbho3SSUDIuItGFJH1V1DuqroKYScvKhoFurJ8KFxpxjRn2YvzOlxYXMnjy+VceRBJs3CT66B85bCp0H+R2NSEJFS4ZVJiEi0oY1WkK3GS1KnM0gr4N3S4CmI8HhEmGIoe5ZUu/QG2HZ/fDhnXDsw35HI5IymvkgIpIBgkloeWVV9DZsSRZueelwYp0sKClU2BsO+SqsfAR2r/Y7GpGUUTIsIpIBwiWhUduwJUksI76F+blMmjA0BdFI3A6b5F0tWHyX35GIpIySYRGRDBApCU11OUKkEd9cMwyvVliT59qwjv1g4FXwycPeUt8iWUDJsIhIBoiUhKa6HCHc8tKF+bn8+vNHsHLKucyePF6JcFs3/CZw9eo7LFlDE+hERDJAuB7F8ZQjJKprRfA96iucXpqe/0cPv4yByx+Cw26ETgf7HZ5IUqm1mohIGmmatJx6aA9e+6hiX19iM6jcUxtXEqpewNkt3PkfULiN6UOuIbfsEjjubz5GJ5IY0VqrqUxCRCRNhOsY8ejba/Y9rqyqpbq2gXsvOTKucoS2MvlO/BHu/K+u6srj28+HVY/C9sU+RSaSGkqGRUTSRCxty1qSxLaVyXfij0jn+VdrJ0JuB/jgpymOSCS1VDMsIpImYk1Og9vFWgdcUlxIeZh9xzr5Lumr5ElSRTr/hZ17ewtxLPqZt3x3t6N8iE4k+TQyLCKSJmJNTkuKC+NahCNSB4hYJt+1lcU+pOWinv/Dvuct0b3gZp+iE0m+ZpNhM/uzmW0ys0URXjcz+62ZLTezD8xsdOLDFBGRcElLU8EkJp464ImjSrnzwpGUFhfG3QtY9cbpL+r5z+8CwybDhhdh00y/QxVJima7SZjZScAu4BHn3Igwr58D3ACcAxwD3OecO6a5A6ubhIhI/KJ1kwgtURg4+Tki/e9u0OpyhmAc4S6vB4+xcsq5Ldq3tDF1e+DpQ6DzIDh9hrdCnUiaidZNotmaYefcDDMri7LJBXiJsgPeNrNiM+vjnNvQomhFRCSiiaNKY0pgI9WBAo3KGYL7jEe4Vlzhji8ZIq8DjPgJzPkmbJgGJWf5HZFIQiWiZrgUWBvyeF3guQOY2bVmNsfM5lRUVCTg0CIiEk4sJRUtLWdorqtFPIt9iD+mzitn3JTpDJz8HOOmTG++xvuQa6BjGSz4Mfi0PoFIsqR0Ap1z7iHn3Bjn3JgePXqk8tAiIlmlaR1oJC1pnxbtPfHUG4s/WjTpMbcARt4K296HtU+kKlSRlEhEMlwO9At53DfwnIiI+GjiqFJmTx7PyinnUhqhbKEl5QyR3lNaXBjXYh/ijxZPeiy7HIpGwPzJUL83iRGKpFYikuGngS8FukocC2xXvbCISNvSmvZpydyXpF6LF1nJyYVRv4JdK2DZA0mITMQfzU6gM7N/AacA3c1sHXALkA/gnHsQeB6vk8RyYA/w5WQFKyIiLRMcrY11cYxoC2nEuy9pW1q1yErJBOgzARb9HA6+EtodlIQIRVKr2dZqyaLWaiIibVO4bhGF+bmqBc4QrT6/lYvghSNgyA1w1G+SGKlI4kRrraYV6EREpBEtpJHZWrPICgDFI+Dgq71SiR0fJzVWkVRotkxCRESyS4trSiVtxNqvOqLDb4PVj8H734NTnklcYCI+0MiwiIg0Eql2NMcs9r60ktaa7UNc2BtG/hTWPwvlz/sTpEiCKBkWEZFGIi3YUe9c7H1pJW3F3Id4yLegy1CY+221WpO0pmRYREQaaVpTmmsHLtuhGuLMFXPNeG4BjL4Pdi2Hj+5NYYQiiaVkWEREDhC6YEdDhK5DqiHOTHHVjJdMgL4XwIe3wx5dKZD0pGRYRCRLNVsXGhCphrglq9dJ2xf3+R59LzTUwfvfTWJUIsmjZFhEJAvFXBeKVpzLNnGf704DYcRPYM1/oPy5FEQoklhKhkVEslA8vYRb3ZdW0kqLzvdhk6BoGLz3DajbnbJYRRJBfYZFRLJQvL2EW92XVtJK3Oc7twCO/gO8ciJ8cAuM/lXyghNJMI0Mi4hkIdUBS8L1PAEGXQtLfwNb5/kdjUjMlAyLiGQh1QFLUhw5Bdp1h3e+AvU1fkcjEhMlwyIiWUh1wJIUBV1h7B9g23z48Bd+RyMSE9UMi4hkKdUBS1L0vQDKLocP7/Dudxvld0QiUWlkWERERBLrqPu8com3r1S5hLR5SoZFREQksdp1g7EPQeVCWHir39GIRKVkWERERBKv72fgkKth8RT49HW/oxGJSMmwiIiIJMfo30DnwfDm5bB3i9/RiISlZFhERESSI78TjPsX7N0E71wDzvkdkcgBlAyLiIhI8nQbDUfcCeumwvIH/Y5G5ABKhkVERCS5Dv0u9DkL5n4HtrzndzQijSgZFhERkeSyHDj+USjsAzMvgurNfkckso+SYREREUm+dgfBCU9A9SZ48zJoqPc7IhFAybCIiIikykFj4OgHYOMr8MHNfkcjAmg5ZhEREUmlQ66GLe96/Ye7HAYHf8nviCTLaWRYREREUmvM/dBrPLx7DWya6Xc0kuWUDIuIiEhq5eTDiU9Ax4Ew87Owc7nfEUkWUzIsIiIiqVfQFU55zluI4/VzvIl1Ij6IKRk2s7PMbKmZLTezyWFev8rMKsxsfuB2TeJDFRERkYzSeRCc/AzsWQevnQW1O/yOSLJQs8mwmeUCDwBnA8OAy8xsWJhNH3fOHRm4/SnBcYqIiEgm6nE8nPgkVC6EN86H+mq/I5IsE8vI8FhguXNuhXOuBngMuCC5YYmIiEjWKDkbjnsENs2AWZ+H+hq/I5IsEksyXAqsDXm8LvBcUxeZ2Qdm9oSZ9UtIdCIiIpIdyi7zehCXPwOzLob6vX5HJFkiURPongHKnHOHAy8Dfwu3kZlda2ZzzGxORUVFgg4tIiIiGWHwdXD0772EeOZFKpmQlIglGS4HQkd6+wae28c5t8U5F/wI9yfgqHA7cs495Jwb45wb06NHj5bEKyIiIpls8HVw9IOw/jmYMRHqdvsdkWS4WJLh94DBZjbQzAqAS4GnQzcwsz4hD88HliQuRBEREckqg78Gx/wJNr4Mr54G1Zv9jkgyWLPJsHOuDrgemIaX5P7bOfehmd1mZucHNvuWmX1oZguAbwFXJStgERERyQKHXA0nPAmVC+CVE2D3ar8jkgxlzjlfDjxmzBg3Z84cX44tIiIiaWLTLHjjM5DbHk6aCt2P8TsiSUNmNtc5Nybca1qBTkRERGI2dV4546ZMZ+Dk5xg3ZTpT55U3/6bW6HkCnDELcgvhlZNh5aPJPZ5kHSXDIiIiEpOp88q56amFlFdW4YDyyipuemph8hPi4uFw1nvQ/Th46wqYNwka6pJ7TMkaSoZFREQkJndPW0pVbX2j56pq67l72tLkH7zdQTD+JRj8DVjyK3j1FNi9ttm3iTRHybCIiIhEFSyNKK+sCvv6+gjPJ1xOvrcwx/H/gG0L4IUjofzZ1BxbMpaSYREREYkotDQikpLiwhRGBJR9Ac6aCx37e5Pr3rkWanemNgbJGEqGRUREJKJwpRGhCvNzmTRhaAojCugyBM58Cw6bBJ/8CZ4fCRunpz4OSXtKhkVERCSiaCUQpcWF3HnhSCaOKk1hRCFy28Oou7xuEzkFMP00eOsqqK7wJx5JS0qGRUREJKJIJRClxYXMnjzev0Q4VI/j4ez5MOwmWP1PeHYofPwgNEQe0RYJUjIsIiIiEU2aMJTC/NxGz4WWRqS873AkeR3gyF/A2Qug65Hw3nWBCXbPg08LjEl6UDIsIiIiEU0cVcqdF46ktLgQo3FphG99h6MpOgzGvwonPAH11fDGuTD9dKh407+YpE3TcswiIiLSIpHarQVLKHxXXwPLH4JFt8HeCuh9Ooy4xVvVTrKKlmMWERGRhIs0uS5lfYebk1sAQ6+HC1bCqF9B5Qfwyonw0jhY86RqigVQMiwiIiItFGlyXcr7DjcnryMcdiOcvxKO+i1Ub4RZF8Mzg2HxXVC9ye8IxUdKhkVERKRFmptc1+bkdYChN8B5y+DEJ6FDX5j/Q5jaF2Zd4k22q6/xO0pJsTy/AxAREZH0FGyrdve0payvrKKkuJBJE4a2jXZr0eTkQr8Lvdv2xbD8j7Dyb7Dm31DQFfpdBAMuhZ4nQ45SpUynCXQiIiIi9TWw8SVY/Tismwp1u6B9T+j3Oeh7AfQ80VvkQ9JStAl0+rgjIiIiklsAped5t7oqWP88rH4MVjwMHz/gJcI9T4Y+E7xbl8PAzO+oJQGUDIuIiIiEyiuE/hd5t7rd8OkbsGEabJwG73/P26ZDX+h5KvQ4DrofB0UjVFKRpnTWRERERCLJ6wil53g3gN2rYcNLgeT4JVj19/3bHTTWS4y7HQ1dD4eOZWDqVdDWKRkWERGRpJk6rzz9JthF03EADPqqd3MOdq+CzW95t4o3YfEvwQX6F+d1guKRUHy4dysaDp0HQ2EflVi0IUqGRUREJCmCyzVX1XrJYXC5ZiBsQpx2ibMZdBro3cq+4D1XtxsqF3kLfARvqx+H5X/Y/768jtBpkJcYdx4MnQdBx/7QoT906OeVaUjKKBkWERGRpLh72tJ9iXBQVW09d09bekCSG2/i3GbldYTux3i3IOegqtxr47bz4/23ygVe5wpX13gf7Xp4SXHHftC+l3dr1xMKm9zPL9YIcwIoGRYREZGkiGe55ngS57Rj5k2469AX+pzZ+LWGWtizFnavhT1rYPeawOM1sPMTr/yiugII0wo3J99LjAu6Bm7FXoIcvF8QuJ9fDPmdILejl6zvu3XyFiLJ8rpmJcMiIiKSFCXFhZSHSXzDLdccT+KcUXLyodPB3i2ShnrYuxn2boLqT73lo6s/3X+/phJqtnlJdM1C737t9thjyC1snCTndvRKNXIKIKed9zW3XePHoc+Fey0nHyzP67Bhefsf9zypzZWBKBkWERGRhAmt+y0qzCc/16it3z+qaXglEOOmTG9UExxP4px1cnK9sojCXsDI2N7TUA91O/cnynW7vHrm4K1+d+PHTW/1VVC3Bxq2QUMN1O/1vjbsbfI4zuWrJ66FvL5x/wiSScmwiIiIJETTut/Kqlryc4yuHfLZtqcWY//F/qY1wZMmDG30XoDC/FwmTRia4u8iQ+Tk7i+VoCx5x3HOK/VomiS7OmioC3yt3f+4XY/kxdJCWo5ZREREEmLclOlhR3dLA6O74V7LNaPBOUqKCzn10B689lFF+nSTkLSh5ZhFREQk6VpS91sfGJQrr6ziybnl3HnhSCXAklLZPX1QREREEiZSfW9JcWFMtb/B7hEiqRRTMmxmZ5nZUjNbbmaTw7zezsweD7z+jpmVJTpQERERadsmTRhKYX5uo+eCdb/hXgsn47tHSJvTbJmEmeUCDwBnAOuA98zsaefc4pDNrga2OecGmdmlwC+BS5IRsIiIiLRNwfKGaKvIBV/LMdtXIhFK3SMk1WKpGR4LLHfOrQAws8eAC4DQZPgC4NbA/SeA+83MnF+z80RERMQXE0eVRqz5DX2taecJUPcI8UcsZRKlwNqQx+sCz4XdxjlXB2wHDkpEgCIiIpJ5Jo4q5c4LR1JaXIjhdZzQ5DnxQ0q7SZjZtcC1AP3790/loUVERKSNiTaKLJIqsYwMlwP9Qh73DTwXdhszywOKgC1Nd+Sce8g5N8Y5N6ZHj7bXdFlEREREskssyfB7wGAzG2hmBcClwNNNtnkauDJw/2JguuqFRURERKSta7ZMwjlXZ2bXA9OAXODPzrkPzew2YI5z7mngYeDvZrYc2IqXMIuIiIiItGkx1Qw7554Hnm/y3E9D7lcDn0tsaCIiIiIiyaUV6EREREQka5lfpb1mVgGs9uXg0B3Y7NOxJXl0XjOTzmvm0rnNTDqvmSudYBSQVQAAAwtJREFUz+0A51zY7g2+JcN+MrM5zrkxfschiaXzmpl0XjOXzm1m0nnNXJl6blUmISIiIiJZS8mwiIiIiGStbE2GH/I7AEkKndfMpPOauXRuM5POa+bKyHOblTXDIiIiIiKQvSPDIiIiIiJKhkVEREQke2VlMmxmnzOzD82swcwyrkVINjKzs8xsqZktN7PJfscjrWdmfzazTWa2yO9YJHHMrJ+ZvWZmiwP/D3/b75gkMcysvZm9a2YLAuf2Z37HJIljZrlmNs/MnvU7lkTLymQYWARcCMzwOxBpPTPLBR4AzgaGAZeZ2TB/o5IE+Ctwlt9BSMLVATc654YBxwLf1L/XjLEXGO+cOwI4EjjLzI71OSZJnG8DS/wOIhmyMhl2zi1xzi31Ow5JmLHAcufcCudcDfAYcIHPMUkrOedmAFv9jkMSyzm3wTn3fuD+Trw/rqX+RiWJ4Dy7Ag/zAzfN0s8AZtYXOBf4k9+xJENWJsOScUqBtSGP16E/riJtnpmVAaOAd/yNRBIlcCl9PrAJeNk5p3ObGX4D/ABo8DuQZMjYZNjMXjGzRWFuGjEUEfGZmXUCngS+45zb4Xc8khjOuXrn3JFAX2CsmY3wOyZpHTM7D9jknJvrdyzJkud3AMninDvd7xgkZcqBfiGP+waeE5E2yMzy8RLhfzjnnvI7Hkk851ylmb2GV/evSbDpbRxwvpmdA7QHupjZo865y32OK2EydmRYssp7wGAzG2hmBcClwNM+xyQiYZiZAQ8DS5xz9/gdjySOmfUws+LA/ULgDOAjf6OS1nLO3eSc6+ucK8P7+zo9kxJhyNJk2Mw+a2brgOOA58xsmt8xScs55+qA64FpeJNx/u2c+9DfqKS1zOxfwFvAUDNbZ2ZX+x2TJMQ44ApgvJnND9zO8TsoSYg+wGtm9gHeIMXLzrmMa8MlmUfLMYuIiIhI1srKkWEREREREVAyLCIiIiJZTMmwiIiIiGQtJcMiIiIikrWUDIuIiIhI1lIyLPL/7daBAAAAAIAgf+tBLooAgC0ZBgBgK33PRjFsvnRYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pylab as plt\n",
    "from sklego.linear_model import LowessRegression\n",
    "\n",
    "n = 100\n",
    "xs = np.linspace(0, np.pi, n)\n",
    "ys = 1 + np.sin(xs) + np.cos(xs**2) + np.random.normal(0, 0.1, n)\n",
    "\n",
    "mod = LowessRegression(sigma=0.1, span=0.5).fit(xs.reshape(-1, 1), ys)\n",
    "\n",
    "xs_new = np.linspace(-1, np.pi + 1, n*2)\n",
    "preds = mod.predict(xs_new.reshape(-1, 1))\n",
    "\n",
    "plt.figure(figsize=(12, 4))\n",
    "plt.scatter(xs, ys)\n",
    "plt.plot(xs_new, preds, color='orange')\n",
    "plt.title(\"Be careful with extrapolation here.\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The line does not look linear but that's because internally, during prediction,\n",
    "many weighted linear regressions are happening. The gif below demonstrates how \n",
    "the data is being weighted when we would make a prediction. \n",
    "\n",
    "![](_static/lowess-rolling.gif)\n",
    "\n",
    "### Details on `sigma`.\n",
    "\n",
    "We'll also show two different prediction outcomes depending on the hyperparameter\n",
    "`sigma`. \n",
    "\n",
    "![](_static/lowess-pred-1.gif)\n",
    "![](_static/lowess-pred-2.gif)\n",
    "\n",
    "You may be tempted now to think that a lower sigma always has a better fit, but you\n",
    "need to be careful here. The data might have gaps and larger sigma values will be\n",
    "able to properly regularize.\n",
    "\n",
    "![](_static/lowess-two-predictions.gif)\n",
    "\n",
    "Note that this regression also works in higher dimensions but the main downside of \n",
    "this approach is that it is really slow when making predictions. \n",
    "\n",
    "If you really want to get advanced there's also a hyperparameter `span` but you'll\n",
    "really need to know what you're doing. It was added for completeness but the authors \n",
    "of this package have yet to find a proper usecase for it. \n",
    "\n",
    "### Details on `span`.\n",
    "\n",
    "The `span` parameter can be used to force that you'll only include a certain\n",
    "  percentage of your dataset. Technically without a `span` you are still using\n",
    "  all the data in your dataset, albeit with small weights if they are far away.\n",
    "  \n",
    "The effect of the `span` parameter on the weights can be seen below.\n",
    "\n",
    "![](_static/lowess-span-1.png)\n",
    "  \n",
    "This will also effect the predictions.</p>\n",
    "\n",
    "![](_static/lowess-span-2.png)\n",
    "\n",
    "You may need to squint your eyes a bit to see it, but lower spans cause more \n",
    "  jiggles and less smooth curves."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ProbWeightRegression\n",
    "\n",
    "> Note, this is a somewhat experimental feature. We found this feature to be plausibly useful but we've not seen it cause a big \"win\" yet. \n",
    "\n",
    "Let's say that you're interested in combining a few models for a regression task. \n",
    "You could put them together in an ensemble. Say we've got predictions $y_1, y_2, y_3$, each of which come from respectable models, then you may want to combine these together in another linear regression. This way the new, hopefully best, prediction $y_*$ is defined via;\n",
    "\n",
    "$$y^* = w_1 y_1 + w_2 y_2 + w_3 y_3$$ \n",
    "\n",
    "This can be a valid way of reweighing. But there's one issue:\n",
    "technically the weights $w_1, w_2, w_3$ can sum to a number that isn't one. Since that's \n",
    "numerically totally possible we need to be aware that we can end up in a strange situation. \n",
    "\n",
    "The `ProbWeightRegression` adresses this. It assumes that every input it receives is\n",
    "the output of a model and it will ensure that they are reweighed with a constraint in mind. For this usecase, it would optimise;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation} \\label{eq1}\n",
    "\\begin{split}\n",
    "\\text{min} & \\sum_i (y_i - (w_1 y_{1,i} + w_2 y_{2,i} + w_3 y_{3,i}))^2 \\\\\n",
    "\\text{subject to} & \\\\\n",
    "& \\text{ } w_1 + w_2 + w_3 = 1 \\\\ \n",
    "& \\text{ } w_1 \\geq 0, w_2 \\geq 0, w_3 \\geq 0\n",
    "\\end{split}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final positivity constraint is optional in our model. \n",
    "\n",
    "Here's an example usage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_regression\n",
    "import pandas as pd\n",
    "\n",
    "X, y = make_regression(n_samples=1000, n_features=10, random_state=42)\n",
    "df = pd.DataFrame(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've turned the array into a dataframe so that we can apply the `ColumnSelector`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.pipeline import Pipeline, FeatureUnion\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "\n",
    "from sklego.meta import EstimatorTransformer\n",
    "from sklego.linear_model import ProbWeightRegression\n",
    "from sklego.preprocessing import ColumnSelector\n",
    "\n",
    "pipe = Pipeline([\n",
    "    ('models', FeatureUnion([\n",
    "        ('path1', Pipeline([\n",
    "            ('select1', ColumnSelector([0,1,2,3,4])),\n",
    "            ('pca', PCA(n_components=3)),\n",
    "            ('linear', EstimatorTransformer(LinearRegression()))\n",
    "        ])),\n",
    "        ('path2', Pipeline([\n",
    "            ('select2', ColumnSelector([5,6,7,8,9])),\n",
    "            ('pca', PCA(n_components=2)),\n",
    "            ('linear', EstimatorTransformer(LinearRegression()))\n",
    "        ]))\n",
    "    ])),\n",
    "    ('prob_weight', ProbWeightRegression())\n",
    "])\n",
    "\n",
    "grid = GridSearchCV(estimator=pipe, param_grid={}, cv=3).fit(df, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>div.sk-top-container {color: black;background-color: white;}div.sk-toggleable {background-color: white;}label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.2em 0.3em;box-sizing: border-box;text-align: center;}div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}div.sk-estimator {font-family: monospace;background-color: #f0f8ff;margin: 0.25em 0.25em;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;}div.sk-estimator:hover {background-color: #d4ebff;}div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 2em;bottom: 0;left: 50%;}div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;}div.sk-item {z-index: 1;}div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;}div.sk-parallel-item {display: flex;flex-direction: column;position: relative;background-color: white;}div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}div.sk-parallel-item:only-child::after {width: 0;}div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0.2em;box-sizing: border-box;padding-bottom: 0.1em;background-color: white;position: relative;}div.sk-label label {font-family: monospace;font-weight: bold;background-color: white;display: inline-block;line-height: 1.2em;}div.sk-label-container {position: relative;z-index: 2;text-align: center;}div.sk-container {display: inline-block;position: relative;}</style><div class=\"sk-top-container\"><div class=\"sk-container\"><div class=\"sk-item sk-dashed-wrapped\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"1f1cc5db-07d1-470d-aa5a-253d7d316602\" type=\"checkbox\" ><label class=\"sk-toggleable__label\" for=\"1f1cc5db-07d1-470d-aa5a-253d7d316602\">GridSearchCV</label><div class=\"sk-toggleable__content\"><pre>GridSearchCV(cv=3,\n",
       "             estimator=Pipeline(steps=[('models',\n",
       "                                        FeatureUnion(transformer_list=[('path1',\n",
       "                                                                        Pipeline(steps=[('select1',\n",
       "                                                                                         ColumnSelector(columns=[0,\n",
       "                                                                                                                 1,\n",
       "                                                                                                                 2,\n",
       "                                                                                                                 3,\n",
       "                                                                                                                 4])),\n",
       "                                                                                        ('pca',\n",
       "                                                                                         PCA(n_components=3)),\n",
       "                                                                                        ('linear',\n",
       "                                                                                         EstimatorTransformer(estimator=LinearRegression()))])),\n",
       "                                                                       ('path2',\n",
       "                                                                        Pipeline(steps=[('select2',\n",
       "                                                                                         ColumnSelector(columns=[5,\n",
       "                                                                                                                 6,\n",
       "                                                                                                                 7,\n",
       "                                                                                                                 8,\n",
       "                                                                                                                 9])),\n",
       "                                                                                        ('pca',\n",
       "                                                                                         PCA(n_components=2)),\n",
       "                                                                                        ('linear',\n",
       "                                                                                         EstimatorTransformer(estimator=LinearRegression()))]))])),\n",
       "                                       ('prob_weight',\n",
       "                                        ProbWeightRegression())]),\n",
       "             param_grid={})</pre></div></div></div><div class=\"sk-parallel\"><div class=\"sk-parallel-item\"><div class=\"sk-item\"><div class=\"sk-serial\"><div class=\"sk-item\"><div class=\"sk-serial\"><div class=\"sk-item sk-dashed-wrapped\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"722adf51-b4eb-4aea-acce-ff4fd0835955\" type=\"checkbox\" ><label class=\"sk-toggleable__label\" for=\"722adf51-b4eb-4aea-acce-ff4fd0835955\">models: FeatureUnion</label><div class=\"sk-toggleable__content\"><pre>FeatureUnion(transformer_list=[('path1',\n",
       "                                Pipeline(steps=[('select1',\n",
       "                                                 ColumnSelector(columns=[0, 1,\n",
       "                                                                         2, 3,\n",
       "                                                                         4])),\n",
       "                                                ('pca', PCA(n_components=3)),\n",
       "                                                ('linear',\n",
       "                                                 EstimatorTransformer(estimator=LinearRegression()))])),\n",
       "                               ('path2',\n",
       "                                Pipeline(steps=[('select2',\n",
       "                                                 ColumnSelector(columns=[5, 6,\n",
       "                                                                         7, 8,\n",
       "                                                                         9])),\n",
       "                                                ('pca', PCA(n_components=2)),\n",
       "                                                ('linear',\n",
       "                                                 EstimatorTransformer(estimator=LinearRegression()))]))])</pre></div></div></div><div class=\"sk-parallel\"><div class=\"sk-parallel-item\"><div class=\"sk-item\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><label>path1</label></div></div><div class=\"sk-serial\"><div class=\"sk-item\"><div class=\"sk-serial\"><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"ee9ceae5-59b7-4c88-9fbf-9505349fa887\" type=\"checkbox\" ><label class=\"sk-toggleable__label\" for=\"ee9ceae5-59b7-4c88-9fbf-9505349fa887\">ColumnSelector</label><div class=\"sk-toggleable__content\"><pre>ColumnSelector(columns=[0, 1, 2, 3, 4])</pre></div></div></div><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"e0129048-6eb8-47bc-a8c0-ae9d96d63cd1\" type=\"checkbox\" ><label class=\"sk-toggleable__label\" for=\"e0129048-6eb8-47bc-a8c0-ae9d96d63cd1\">PCA</label><div class=\"sk-toggleable__content\"><pre>PCA(n_components=3)</pre></div></div></div><div class=\"sk-item sk-dashed-wrapped\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"fe8293b5-d914-4a59-8c39-e2b4e18851c8\" type=\"checkbox\" ><label class=\"sk-toggleable__label\" for=\"fe8293b5-d914-4a59-8c39-e2b4e18851c8\">linear: EstimatorTransformer</label><div class=\"sk-toggleable__content\"><pre>EstimatorTransformer(estimator=LinearRegression())</pre></div></div></div><div class=\"sk-parallel\"><div class=\"sk-parallel-item\"><div class=\"sk-item\"><div class=\"sk-serial\"><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"41c53c2a-d2c2-4241-b549-8f7653a0029f\" type=\"checkbox\" ><label class=\"sk-toggleable__label\" for=\"41c53c2a-d2c2-4241-b549-8f7653a0029f\">LinearRegression</label><div class=\"sk-toggleable__content\"><pre>LinearRegression()</pre></div></div></div></div></div></div></div></div></div></div></div></div></div><div class=\"sk-parallel-item\"><div class=\"sk-item\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><label>path2</label></div></div><div class=\"sk-serial\"><div class=\"sk-item\"><div class=\"sk-serial\"><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"86b9c93e-12be-4fa6-865f-bb874ff99129\" type=\"checkbox\" ><label class=\"sk-toggleable__label\" for=\"86b9c93e-12be-4fa6-865f-bb874ff99129\">ColumnSelector</label><div class=\"sk-toggleable__content\"><pre>ColumnSelector(columns=[5, 6, 7, 8, 9])</pre></div></div></div><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"ec46eb32-4807-466f-a47b-dba193ff5315\" type=\"checkbox\" ><label class=\"sk-toggleable__label\" for=\"ec46eb32-4807-466f-a47b-dba193ff5315\">PCA</label><div class=\"sk-toggleable__content\"><pre>PCA(n_components=2)</pre></div></div></div><div class=\"sk-item sk-dashed-wrapped\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"5efe2360-3f3b-4c67-bc8f-e2619f06d229\" type=\"checkbox\" ><label class=\"sk-toggleable__label\" for=\"5efe2360-3f3b-4c67-bc8f-e2619f06d229\">linear: EstimatorTransformer</label><div class=\"sk-toggleable__content\"><pre>EstimatorTransformer(estimator=LinearRegression())</pre></div></div></div><div class=\"sk-parallel\"><div class=\"sk-parallel-item\"><div class=\"sk-item\"><div class=\"sk-serial\"><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"6ceb127c-c870-4944-b6c3-5f47c0b5af1b\" type=\"checkbox\" ><label class=\"sk-toggleable__label\" for=\"6ceb127c-c870-4944-b6c3-5f47c0b5af1b\">LinearRegression</label><div class=\"sk-toggleable__content\"><pre>LinearRegression()</pre></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"8d1aaabb-dda3-4fe0-99dc-c4b60597b332\" type=\"checkbox\" ><label class=\"sk-toggleable__label\" for=\"8d1aaabb-dda3-4fe0-99dc-c4b60597b332\">ProbWeightRegression</label><div class=\"sk-toggleable__content\"><pre>ProbWeightRegression()</pre></div></div></div></div></div></div></div></div></div></div></div></div>"
      ],
      "text/plain": [
       "GridSearchCV(cv=3,\n",
       "             estimator=Pipeline(steps=[('models',\n",
       "                                        FeatureUnion(transformer_list=[('path1',\n",
       "                                                                        Pipeline(steps=[('select1',\n",
       "                                                                                         ColumnSelector(columns=[0,\n",
       "                                                                                                                 1,\n",
       "                                                                                                                 2,\n",
       "                                                                                                                 3,\n",
       "                                                                                                                 4])),\n",
       "                                                                                        ('pca',\n",
       "                                                                                         PCA(n_components=3)),\n",
       "                                                                                        ('linear',\n",
       "                                                                                         EstimatorTransformer(estimator=LinearRegression()))])),\n",
       "                                                                       ('path2',\n",
       "                                                                        Pipeline(steps=[('select2',\n",
       "                                                                                         ColumnSelector(columns=[5,\n",
       "                                                                                                                 6,\n",
       "                                                                                                                 7,\n",
       "                                                                                                                 8,\n",
       "                                                                                                                 9])),\n",
       "                                                                                        ('pca',\n",
       "                                                                                         PCA(n_components=2)),\n",
       "                                                                                        ('linear',\n",
       "                                                                                         EstimatorTransformer(estimator=LinearRegression()))]))])),\n",
       "                                       ('prob_weight',\n",
       "                                        ProbWeightRegression())]),\n",
       "             param_grid={})"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn import set_config\n",
    "set_config(display='diagram')\n",
    "grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the `ProbWeightRegression` indeeds sums to one. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.03102466, 0.96897535])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid.best_estimator_[1].coefs_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Least Absolute Deviation Regression\n",
    "\n",
    "Imagine that you have a dataset with some outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x26699696348>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA54AAAD4CAYAAACEyjk9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAT+0lEQVR4nO3df4xd6XkX8O+DvZVm28AEdvrDs1m8VMFQuolcBigs0LBp5E2IErMqUkN/hBBkIaBUiLpZq1LzR4V2kfnRotJGVrpsK6Lkj9S4QaV1V12lC20SmI03cZqt2yhptzsbWC/BBaWDsuu8/OFx1p6MPddz73tnzr2fj2R57rln5jy6eu/M/Z73Pc+p1loAAACglz+y2wUAAAAw2wRPAAAAuhI8AQAA6ErwBAAAoCvBEwAAgK72T/Ngd9xxRzt48OA0DwkAAMCUPPnkky+01pY2b59q8Dx48GBWV1eneUgAAACmpKp+b6vtltoCAADQleAJAABAV4InAAAAXQmeAAAAdCV4AgAA0NVUu9oCAMBecebcWk6evZDnLq3nwOJCjh85lKOHl3e7LJhJgicAAHPnzLm1nDh9PusvXk6SrF1az4nT55NE+IQOLLUFAGDunDx74Suh86r1Fy/n5NkLu1QRzDbBEwCAufPcpfVb2g6MR/AEAGDuHFhcuKXtwHgETwAA5s7xI4eycNu+67Yt3LYvx48c2qWKYLZpLgQAwNy52kBIV1uYDsETAIC5dPTwsqAJU2KpLQAAAF0JngAAAHQleAIAANCV4AkAAEBXgicAAABdCZ4AAAB0JXgCAADQ1bbBs6oeqarnq+pTWzz3Q1XVquqOPuUBAAAwdKPMeD6a5P7NG6vqVUnekOSZCdcEAADADNk2eLbWnkjyhS2e+jdJfjhJm3RRAAAAzI4dXeNZVW9JstZa+8QI+x6rqtWqWr148eJODgcAAMCA3XLwrKrbk/xIkh8dZf/W2qnW2kprbWVpaelWDwcAAMDA7WTG85uT3J3kE1X1u0nuTPLxqvrGSRYGAADAbNh/q9/QWjuf5OuvPt4InyuttRcmWBcAAAAzYpTbqbw/yUeSHKqqZ6vqnf3LAgAAYFZsO+PZWnvbNs8fnFg1AAAAzJwddbUFAACAUQmeAAAAdCV4AgAA0JXgCQAAQFeCJwAAAF0JngAAAHQleAIAANCV4AkAAEBXgicAAABdCZ4AAAB0JXgCAADQleAJAABAV4InAAAAXQmeAAAAdCV4AgAA0JXgCQAAQFeCJwAAAF0JngAAAHQleAIAANDVtsGzqh6pquer6lPXbDtZVb9VVZ+sqv9YVYt9ywQAAGCoRpnxfDTJ/Zu2PZbkW1trr0ny20lOTLguAAAAZsS2wbO19kSSL2za9iuttZc2Hn40yZ0dagMAAGAGTOIaz7+X5Jcm8HMAAACYQWMFz6r6kSQvJXnfTfY5VlWrVbV68eLFcQ4HAADAAO04eFbV25O8Ocn3tNbajfZrrZ1qra201laWlpZ2ejgAAAAGav9Ovqmq7k/yriTf0Vr7w8mWBAAAwCwZ5XYq70/ykSSHqurZqnpnkp9M8ookj1XVU1X1ns51AgAAMFDbzni21t62xeaf6VALAAAAM2gSXW0BAADghgRPAAAAuhI8AQAA6ErwBAAAoCvBEwAAgK4ETwAAALoSPAEAAOhK8AQAAKArwRMAAICuBE8AAAC6EjwBAADoSvAEAACgK8ETAACArgRPAAAAuhI8AQAA6ErwBAAAoCvBEwAAgK4ETwAAALoSPAEAAOhK8AQAAKCrbYNnVT1SVc9X1aeu2fbHq+qxqvqdjf9f2bdMAAAAhmqUGc9Hk9y/aduDSX61tfbqJL+68RgAAAC+yrbBs7X2RJIvbNr81iQ/u/H1zyY5OuG6AAAAmBE7vcbzG1prn0+Sjf+//kY7VtWxqlqtqtWLFy/u8HAAAAAMVffmQq21U621ldbaytLSUu/DAQAAsMfsNHj+z6r6piTZ+P/5yZUEAADALNlp8PxQkrdvfP32JL8wmXIAAACYNaPcTuX9ST6S5FBVPVtV70zycJI3VNXvJHnDxmMAAAD4Kvu326G19rYbPPX6CdcCAADADOreXAgAAID5JngCAADQleAJAABAV4InAAAAXQmeAAAAdCV4AgAA0JXgCQAAQFeCJwAAAF0JngAAAHQleAIAANCV4AkAAEBXgicAAABdCZ4AAAB0JXgCAADQleAJAABAV4InAAAAXQmeAAAAdCV4AgAA0JXgCQAAQFeCJwAAAF3tH+ebq+qfJvn7SVqS80ne0Vr7f5MoDAC43plzazl59kKeu7SeA4sLOX7kUI4eXt7tsgBgWzue8ayq5ST/JMlKa+1bk+xL8t2TKgwAeNmZc2s5cfp81i6tpyVZu7SeE6fP58y5td0uDQC2Ne5S2/1JFqpqf5Lbkzw3fkkAwGYnz17I+ouXr9u2/uLlnDx7YZcqAoDR7Th4ttbWkvzLJM8k+XySP2it/crm/arqWFWtVtXqxYsXd14pAMyx5y6t39J2ANhLxllq+8okb01yd5IDSb62qr53836ttVOttZXW2srS0tLOKwWAOXZgceGWtgPAXjLOUtvvTPK51trF1tqLSU4n+SuTKQsAuNbxI4eycNu+67Yt3LYvx48c2qWKAGB043S1fSbJt1fV7UnWk7w+yepEqgIArnO1e62utgAM0Y6DZ2vtY1X1wSQfT/JSknNJTk2qMADgekcPLwuaAAzSWPfxbK29O8m7J1QLAAAAM2jc26kAAADATQmeAAAAdCV4AgAA0JXgCQAAQFeCJwAAAF0JngAAAHQleAIAANCV4AkAAEBXgicAAABdCZ4AAAB0JXgCAADQleAJAABAV4InAAAAXQmeAAAAdCV4AgAA0JXgCQAAQFeCJwAAAF0JngAAAHQleAIAANCV4AkAAEBXYwXPqlqsqg9W1W9V1dNV9ZcnVRgAAACzYf+Y3/8TSX65tfZdVfU1SW6fQE0AAADMkB0Hz6r6o0n+epK/mySttS8l+dJkygIAAGBWjLPU9k8luZjk31fVuap6b1V97eadqupYVa1W1erFixfHOBwAAABDNE7w3J/k25L8dGvtcJIvJnlw806ttVOttZXW2srS0tIYhwMAAGCIxgmezyZ5trX2sY3HH8yVIAoAAABfsePg2Vr7H0l+v6oObWx6fZJPT6QqAAAAZsa4XW1/IMn7NjrafjbJO8YvCQAAgFkyVvBsrT2VZGVCtQAAADCDxrnGEwAAALYleAIAANCV4AkAAEBXgicAAABdCZ4AAAB0JXgCAADQleAJAABAV4InAAAAXQmeAAAAdCV4AgAA0JXgCQAAQFeCJwAAAF0JngAAAHQleAIAANCV4AkAAEBXgicAAABd7d/tAgAAgJs7c24tJ89eyHOX1nNgcSHHjxzK0cPLu10WjEzwBACAPezMubWcOH0+6y9eTpKsXVrPidPnk0T4ZDAstQUAgD3s5NkLXwmdV62/eDknz17YpYrg1gmeAACwhz13af2WtsNeNPZS26ral2Q1yVpr7c3jlwQAAFx1YHEha1uEzAOLC7tQzexw3ex0TWLG8weTPD2BnwMAAGxy/MihLNy277ptC7fty/Ejh3apouG7et3s2qX1tLx83eyZc2u7XdrMGit4VtWdSf5mkvdOphwAAOBaRw8v56EH7sny4kIqyfLiQh564B6zc2Nw3ez0jbvU9seT/HCSV9xoh6o6luRYktx1111jHg4AAObP0cPLguYEuW52+nY841lVb07yfGvtyZvt11o71Vpbaa2tLC0t7fRwAAAAE3Gj62NdN9vPOEtt703ylqr63SQfSHJfVf2HiVQFAADQietmp2/HS21bayeSnEiSqnpdkh9qrX3vhOoCAAAGaAjdYq/Ws9frvGoIr+l2xr6dCgAAQPJyt9irjXuudotNsueC0lCumx3Sa3ozEwmerbUPJ/nwJH4WADBss3BmnvEYA/PrZt1ijYGdmZXX1IwnADAxs3Jmnp0zBuabbrGTNyuv6Vj38QQAuNbQ7o135txa7n348dz94C/m3ocfd/P4CRjaGGCydIudvFl5TQVPAGBihnRm/urM3Nql9bS8PDMnfI5nSGOAydMtdvJm5TUVPAGAiRnSmXkzc30MaQwweUcPL+ehB+7J8uJCKsny4kIeeuAey6zHMCuvqWs8AYCJOX7k0HXX9yV798z80GbmhtKwZ0hjgD6G0i12SGbhNRU8AYCJGdK98Q4sLmRti5C5F2fmhtSwZ0hjAJieaq1N7WArKyttdXV1ascDALiRzWEuuTIztxeXsN378ONbhuTlxYX8+oP37UJFAFurqidbayubt5vxBADm0pBm5oa2LBhgM8ETAJhbQ7luakjLggG2oqstAMAeNyu3UwDmlxlPAIA9bkjLggG2InhChtOiHpg873+GYijLggG2Ingy94bUot4HZJisIb3/AWDIXOPJ3Dt59sJ1rfSTZP3Fyzl59sIuVbS1qx+Q1y6tp+XlD8hnzq3tdmkwWEN5/wPA0JnxZO4NpUX9zT4gm5kZj5nk+TWU9z8ADJ0ZT+bejVrR77UW9T4g92Emeb4N5f0PAEMneDL3htKi3gfkPiy1nG9Def8DwNBZasvcG0qL+uNHDl3XBCXxAXkSzCT3MZTly0N5/wPA0AmekGG0qPcBuY8DiwtZ2yJkmkneuaF1ih3C+x8Ahk7whAHxAXnyhjaTPISZRI2wAIDNdhw8q+pVSX4uyTcm+XKSU621n5hUYQDTMKSZ5KHMJFq+DABsNs6M50tJ/llr7eNV9YokT1bVY621T0+oNoCpGMpM8lBmEi1fBgA223FX29ba51trH9/4+v8meTrJ3vnkc4vOnFvLvQ8/nrsf/MXc+/DjbqUA7DlDmUnUKRYA2Gwi13hW1cEkh5N8bIvnjiU5liR33XXXJA43cUNZvjY0Q7gWDYZkKDOJQ1q+DABMR7XWxvsBVV+X5NeS/PPW2umb7buystJWV1fHOl4P9z78+JYf5pYXF/LrD963CxUN3+Ywn1yZ8XjogXt8+IQd8r4CAPa6qnqytbayefuOl9pu/NDbkvx8kvdtFzr3sqEsXxuSm12LBuzM0cPLeeiBe7K8uJDKlZNjQicAMATjdLWtJD+T5OnW2r+eXEnTN5Tla0MizEMfQ2mEBABwrXFmPO9N8n1J7quqpzb+vWlCdU2VRhiTd6PQLswDAMD82fGMZ2vtvyapCdayazTCmLzjRw5teS2aMA8AAPNnIl1tZ8FQlq8NpVOsMA8AAFwleA7I0G77MpQwTx9DOUkCAEB/Y3W1Zbp0imUorp4kWbu0npaXT5KcObe226UBALALBM8B0SmWoXCSBACAawmeA6JTLEPhJAkAANcSPAfEbV8YCidJAAC4luA5IEcPL+ehB+7J8uJCKsny4kIeeuAeDVvYc5wkAQDgWrraDoxOsQyB2+kAAHAtwRPowkkSAACustQWAACArgRPAAAAuhI8AQAA6ErwBAAAoCvBEwAAgK6qtTa9g1VdTPJ7UzvgztyR5IXdLoK5ZxyyFxiH7BXGInuBccheMIRx+Cdba0ubN041eA5BVa221lZ2uw7mm3HIXmAcslcYi+wFxiF7wZDHoaW2AAAAdCV4AgAA0JXg+dVO7XYBEOOQvcE4ZK8wFtkLjEP2gsGOQ9d4AgAA0JUZTwAAALoSPAEAAOhqLoNnVd1fVReq6jNV9eAWz1dV/duN5z9ZVd+2G3Uy20YYh9+zMf4+WVW/UVWv3Y06mX3bjcVr9vsLVXW5qr5rmvUxH0YZh1X1uqp6qqp+s6p+bdo1MvtG+Nv8x6rqP1XVJzbG4Tt2o05mW1U9UlXPV9WnbvD8ILPK3AXPqtqX5N8leWOSb0nytqr6lk27vTHJqzf+HUvy01Mtkpk34jj8XJLvaK29JsmPZcAXk7N3jTgWr+73L5KcnW6FzINRxmFVLSb5qSRvaa39uSR/e+qFMtNG/H34j5J8urX22iSvS/Kvquprploo8+DRJPff5PlBZpW5C55J/mKSz7TWPtta+1KSDyR566Z93prk59oVH02yWFXfNO1CmWnbjsPW2m+01v73xsOPJrlzyjUyH0b5nZgkP5Dk55M8P83imBujjMO/k+R0a+2ZJGmtGYtM2ijjsCV5RVVVkq9L8oUkL023TGZda+2JXBlbNzLIrDKPwXM5ye9f8/jZjW23ug+M41bH2DuT/FLXiphX247FqlpO8reSvGeKdTFfRvmd+KeTvLKqPlxVT1bV90+tOubFKOPwJ5P82STPJTmf5Adba1+eTnnwFYPMKvt3u4BdUFts23xPmVH2gXGMPMaq6m/kSvD8q10rYl6NMhZ/PMm7WmuXr5zkh4kbZRzuT/Lnk7w+yUKSj1TVR1trv927OObGKOPwSJKnktyX5JuTPFZV/6W19n96FwfXGGRWmcfg+WySV13z+M5cOWt1q/vAOEYaY1X1miTvTfLG1tr/mlJtzJdRxuJKkg9shM47krypql5qrZ2ZTonMgVH/Nr/QWvtiki9W1RNJXptE8GRSRhmH70jycGutJflMVX0uyZ9J8t+mUyIkGWhWmceltv89yaur6u6Ni8G/O8mHNu3zoSTfv9Ex6tuT/EFr7fPTLpSZtu04rKq7kpxO8n3O6NPRtmOxtXZ3a+1ga+1gkg8m+YdCJxM2yt/mX0jy16pqf1XdnuQvJXl6ynUy20YZh8/kyqx7quobkhxK8tmpVgkDzSpzN+PZWnupqv5xrnRm3Jfkkdbab1bVP9h4/j1J/nOSNyX5TJI/zJWzWzAxI47DH03yJ5L81MZM00uttZXdqpnZNOJYhK5GGYettaer6peTfDLJl5O8t7W25a0GYCdG/H34Y0kerarzubLc8V2ttRd2rWhmUlW9P1e6Jt9RVc8meXeS25JhZ5W6slIAAAAA+pjHpbYAAABMkeAJAABAV4InAAAAXQmeAAAAdCV4AgAA0JXgCQAAQFeCJwAAAF39fy1OUk9E2khEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "np.random.seed(0)\n",
    "X = np.linspace(0, 1, 20)\n",
    "y = 3*X + 1 + 0.5*np.random.randn(20)\n",
    "X = X.reshape(-1, 1)\n",
    "\n",
    "y[10] = 8\n",
    "y[15] = 15\n",
    "\n",
    "plt.figure(figsize=(16, 4))\n",
    "plt.scatter(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple linear regression will not do a good job since it is distracted by the outliers. That is because it optimizes the mean squared error\n",
    "$$ \\sum_i \\left(y_i-\\textrm{model}(x_i)\\right)^2 $$\n",
    "\n",
    "which penalizes a few large errors more than many tiny errors. For example, if y-model(x) = 4 for some single observation, the MSE here is 16. If there are two observations with y_1 - model(x_1) = 2 and y_2 - model(x_2) = 2, the MSE is 8 in total, which is less than for one larger error. Note that the sum of the errors is the same in both cases.\n",
    "\n",
    "Hence, linear regression does the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x266ab252948>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA54AAAD4CAYAAACEyjk9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAbWElEQVR4nO3de5Cd50Hf8d+j1cWr6+qWi2XLMo4QjqM4Sraum7QlJWUcKCUupTOk3Jqm4+m0hUxbDHGZgT+YDumoF+hQYFxIA1Mm/BFcQ2cAkYGBtMGhlXGIA0aQhhIsp1i2dbGt1W319I9dIR15z+45u+c9u+85n8+MRrvvHj3nOTuPpP2e9z3PKbXWAAAAQFPWrfYEAAAAGG3CEwAAgEYJTwAAABolPAEAAGiU8AQAAKBR64d5Z3v27KkHDhwY5l0CAAAwJE8++eQLtda9Nx8fangeOHAgx48fH+ZdAgAAMCSllD9d6LhLbQEAAGiU8AQAAKBRwhMAAIBGCU8AAAAaJTwBAABo1FB3tQUAgLXi8adO5uixE3nuzExunZrMww8cyoNH9q32tGAkCU8AAMbO40+dzCOPPZ2Zy7NJkpNnZvLIY08nifiEBrjUFgCAsXP02Im/iM5rZi7P5uixE6s0IxhtwhMAgLHz3JmZvo4DKyM8AQAYO7dOTfZ1HFgZ4QkAwNh5+IFDmdww0XFscsNEHn7g0CrNCEabzYUAABg71zYQsqstDIfwBABgLD14ZJ/QhCFxqS0AAACNEp4AAAA0SngCAADQKOEJAABAo4QnAAAAjRKeAAAANEp4AgAA0Kglw7OU8tFSyvOllM8v8LXvKaXUUsqeZqYHAABA2/VyxvNjSd5788FSyu1JvjbJlwY8JwAAAEbIkuFZa/1UkpcW+NJ/SPK9SeqgJwUAAMDoWNZrPEsp35jkZK3193q47UOllOOllOOnTp1azt0BAADQYn2HZyllc5LvT/IDvdy+1vporXW61jq9d+/efu8OAACAllvOGc+7ktyZ5PdKKf83yW1JfreU8oZBTgwAAIDRsL7fP1BrfTrJ6659Ph+f07XWFwY4LwAAAEZEL2+n8vEkTyQ5VEp5tpTyweanBQAAwKhY8oxnrfX9S3z9wMBmAwAAwMhZ1q62AAAA0CvhCQAAQKOEJwAAAI0SngAAADRKeAIAANAo4QkAAECjhCcAAACNEp4AAAA0SngCAADQKOEJAABAo4QnAAAAjRKeAAAANEp4AgAA0CjhCQAAQKOEJwAAAI0SngAAADRKeAIAANAo4QkAAECjhCcAAACNWjI8SykfLaU8X0r5/A3HjpZS/rCU8rlSyn8rpUw1O00AAADaqpcznh9L8t6bjn0yyVtqrW9N8kdJHhnwvAAAABgRS4ZnrfVTSV666div1VqvzH/6mSS3NTA3AAAARsAgXuP5D5P8ygDGAQAAYAStKDxLKd+f5EqSn1vkNg+VUo6XUo6fOnVqJXcHAABACy07PEsp35nkG5J8a621drtdrfXRWut0rXV67969y707AAAAWmr9cv5QKeW9Sb4vyVfXWs8PdkoAAACMkl7eTuXjSZ5IcqiU8mwp5YNJfizJtiSfLKV8tpTykw3PEwAAgJZa8oxnrfX9Cxz+6QbmAgAAwAgaxK62AAAA0JXwBAAAoFHCEwAAgEYJTwAAABolPAEAAGiU8AQAAKBRwhMAAIBGCU8AAAAaJTwBAABolPAEAACgUcITAACARglPAAAAGiU8AQAAaJTwBAAAoFHCEwAAgEYJTwAAABolPAEAAGiU8AQAAKBRwhMAAIBGCU8AAAAatWR4llI+Wkp5vpTy+RuO7SqlfLKU8sfzv+9sdpoAAAC0VS9nPD+W5L03Hftwkl+vtR5M8uvznwMAAMBrLBmetdZPJXnppsPvS/Iz8x//TJIHBzwvAAAARsRyX+P5+lrrl5Nk/vfXdbthKeWhUsrxUsrxU6dOLfPuAAAAaKvGNxeqtT5aa52utU7v3bu36bsDAABgjVlueP55KeWNSTL/+/ODmxIAAACjZLnh+UtJvnP+4+9M8ouDmQ4AAACjppe3U/l4kieSHCqlPFtK+WCSjyT52lLKHyf52vnPAQAA4DXWL3WDWuv7u3zpPQOeCwAAACOo8c2FAAAAGG/CEwAAgEYJTwAAABolPAEAAGiU8AQAAKBRwhMAAIBGCU8AAAAaJTwBAABolPAEAACgUcITAACARglPAAAAGiU8AQAAaJTwBAAAoFHCEwAAgEYJTwAAABolPAEAAGiU8AQAAKBRwhMAAIBGCU8AAAAaJTwBAABo1PqV/OFSyj9P8o+S1CRPJ/lArfXCICYGAHR6/KmTOXrsRJ47M5Nbpybz8AOH8uCRfas9LQBY0rLPeJZS9iX57iTTtda3JJlI8i2DmhgAcN3jT53MI489nZNnZlKTnDwzk0ceezqPP3VytacGAEta6aW265NMllLWJ9mc5LmVTwkAuNnRYycyc3m249jM5dkcPXZilWYEAL1bdnjWWk8m+bdJvpTky0nO1lp/7ebblVIeKqUcL6UcP3Xq1PJnCgBj7LkzM30dB4C1ZCWX2u5M8r4kdya5NcmWUsq33Xy7WuujtdbpWuv03r17lz9TABhjt05N9nUcANaSlVxq+zeT/Emt9VSt9XKSx5K8czDTAgBu9PADhzK5YaLj2OSGiTz8wKFVmhEA9G4lu9p+Kcn9pZTNSWaSvCfJ8YHMCgDocG33WrvaAtBGyw7PWuvvlFI+keR3k1xJ8lSSRwc1MQCg04NH9glNAFppRe/jWWv9wSQ/OKC5AAAAMIJW+nYqAAAAsCjhCQAAQKOEJwAAAI0SngAAADRKeAIAANAo4QkAAECjhCcAAACNEp4AAAA0SngCAADQKOEJAABAo4QnAAAAjRKeAAAANEp4AgAA0CjhCQAAQKOEJwAAAI0SngAAADRKeAIAANAo4QkAAECjhCcAAACNEp4AAAA0akXhWUqZKqV8opTyh6WUZ0opf2VQEwMAAGA0rF/hn//RJL9aa/3mUsrGJJsHMCcAAABGyLLDs5SyPclfT/IPkqTWeinJpcFMCwAAgFGxkkttvyLJqST/pZTyVCnlp0opW26+USnloVLK8VLK8VOnTq3g7gAAAGijlYTn+iRvT/ITtdYjSV5N8uGbb1RrfbTWOl1rnd67d+8K7g4AAIA2WslrPJ9N8myt9XfmP/9EFghPAAAA+nD+fPLCC8mLL3b+PjubfOhDqz27ZVl2eNZa/18p5c9KKYdqrSeSvCfJHwxuagAAAC1Wa/eI7Pbxiy8mMzMLjzc1NX7hOe+7kvzc/I62X0zygZVPCQAAYI2pNXn11f4j8sKF7mPu2pXs3p3s2ZPcfnvytrfNfXzt2I0f7949d/uWWlF41lo/m2R6QHMBAABoXq3JK6/0H5EXLy48XimdEXnHHcnb3754RO7cmaxf6XnA9hifRwoAAIyeWpOXX+4/Ii91eSfIUubC8Fok3nlnMj3dPSL37Jm7BHZiYriPu2WEJwAAsDbUmpw7t3gwLnTs8uWFx1u3rjMi77orue++pSNy3Ure/IOFCE8AAGDwak3Onu0/Iq9cWXi8iYnOiDx4MLn//oUvY712bMcOEblGCE8AAGBxV6/2H5EvvbR4RN4YjIcOJe985+IRuX27iGwx4QkAAOPk6tXkzJn+I3J2duHx1q/vDMa77+6+oc6NEVnKcB83q0p4AgBAW129mpw+3X9EXr268HgbNnQG4z33LB2R27aJSJYkPAEAYC2YnV1eRNa68HgbN3YG4+HDS0fk1q0ikkYITwAAGLTZ2bko7CciT5/uHpGbNnUG4733dt+V9drHW7aISNYM4QkAAIu5cmV5EdnNLbd0RuL+/UtH5ObNIpJWE54AAIyPy5f7j8gzZ7qPNznZGYkHDnS/jPXGiIQxIzwBAGiny5cXD8aFgvLs2e7jbd7cGYlf8RWLR+Tu3SISeiQ8AQBYfZcu9R+R5851H2/Lls5IfNOblo7IycnhPV4YM8ITAIDBunix/4h8+eXu423b1hmMX/mVS0fkLbcM7/ECSxKeAAB0d+FC/xH5yivdx9u+vTMUv+qrlo7ITZuG93iBRghPAIBxMTPT+4Y6135/9dXu4+3YcT0UX/e65M1v7r6hzrWI3LhxeI8XWDOEJwBAG50/339Enj/ffbypqeuR+IY3JG95y+IRuWuXiAR6JjwBAFZTrcuLyJmZ7mPu3Hk9Em+9NXnrW5eOyA0bhveY6dvjT53M0WMn8tyZmdw6NZmHHziUB4/sW+1pQc+EJwDAoNQ6d2lqvxF54UL3MXftuh6Jt9+evO1tS0fkej/ijZLHnzqZRx57OjOXZ5MkJ8/M5JHHnk4S8Ulr+FcJAGAhtc5tktNvRF68uPB4pXRG5B13JG9/++IRuXOniCRHj534i+i8ZubybI4eOyE8aQ3/kgEAo6/Wubfr6DciL11aeLxSrm+Ws2dPcuedyfT00hE5MTHcx81IeO7MwpdVdzsOa9GKw7OUMpHkeJKTtdZvWPmUAAAWUWty7lz/EXn58sLjrVvXGZF33ZXcd1/3iNyzZ24jnnXrhvu4GVu3Tk3m5AKReevU5CrMZnR43exwDeKM54eSPJNk+wDGAgDGSa3J2bP9R+SVKwuPNzHRGZEHDyb33794RO7YISJZ0x5+4FDHazyTZHLDRB5+4NAqzqrdvG52+FYUnqWU25L8rST/Osm/GMiMAIB2unq1/4h88cXFI/LGSDx0qPtlrNeObd8uIhk510LI2bnB8brZ4VvpGc8fSfK9SbZ1u0Ep5aEkDyXJ/v37V3h3AMBQXL2anDnTe0S+8ELy0kvJ7OzC461f3xmMd9/dW0SWMtzHDWvUg0f2CaIB8rrZ4Vt2eJZSviHJ87XWJ0sp7+52u1rro0keTZLp6em63PsDAJbp6tXk9On+I/Lq1YXH27ChMxjvuWfpiNy2TUQCa4bXzQ7fSs54vivJN5ZSvj7JLUm2l1L+a6312wYzNQDgNWZn+4/I06e7R+TGjZ3BePjw0hG5dauIBFrN62aHb9nhWWt9JMkjSTJ/xvN7RCcA9GF2du7MYr8RWbtcQLRpU2cw3ntv9w11rn28ZYuIBAaqDbvFtu11s234ni7F+3gCwCBcudI9IrsF5Zkz3SPylls6I/HIkaUjcvNmEQmsqjbtFtuW18226Xu6mIGEZ631N5P85iDGAoBVd/ny8iKym8nJzki8447eIrKlRuGZeVbGGhhfdosdvFH5njrjCcBou3z5+tt29BqRZ892H2/z5s5IvPPOxSNy9+5WR2S/RuWZeZbPGhhvdosdvFH5ngpPANrj0qXXhuJSEXnuXPfxtmzpjMS77uq+oc61iJy04+Fi2vbMvDNzg9e2NcBg2S128Ebleyo8AVgdFy/2H5Evv9x9vK1bOyPx4MGlI/KWW4b3eMdEm56Zd2auGW1aAwye3WIHb1S+p8ITgJW7cKH/iHzlle7jbd/eGYyHDi0dkZs2De/x0lWbnpl3Zq4ZbVoDDF7bdottg1H5ngpPADpduNDbW3vceOzVV7uPt2PH9VDcuze5++6lI3LjxuE9XgaqTc/Mt+3MXFsuC27TGqAZbdkttk1G4XsqPAFG2cxM/xF5/nz38aamrofi61+f3HPP4hG5a5eIHDNtema+TWfm2nRZcJvWADA8pXZ7/7AGTE9P1+PHjw/t/gBGyvnzvUXkjR/PLHLmZufO7sF488fXInLDhuE9XmjYzTGXzJ2Z++FvOrzmIuldH/mNBSN539RkPv3hr1mFGQEsrJTyZK11+ubjzngCDFutvUfkjccuXOg+5q5d1yPxttuSt71t6Yhc778Axlubzsy17bJggJv5qQNgJWqde31jvxF58eLC45XSGZH79ydHjiwekTt3ikhYpra8bqpNlwUDLMRPKgDX1Dq302q/l7NeurTweKVc3yxnz57kwIHkHe9YOiInJob6sIG1z4Y9QNsJT2A01Tr3no/9RuTlywuPt25dZ0TedVdy332LR+TUlIgEBqJNlwUDLER4QtqzRf3YqjU5d66/iHzxxe4ROTHRGZEHDyb337/4Jjs7dszFJyPH33/aoi2XBQMsRHgy9tq0Rf1I/IBca3L2bP8ReeXKwuNNTHRG4qFDS+/Uun27iCRJu/7+A0CbCU/G3tFjJzpeM5MkM5dnc/TYiTX1g+ea/AH56tWFI3KxoHzxxWR2duHx1q/vjMS77+5+GeuNEVnKcB83I6Mtf/8BoO2EJ2OvLVvUN/4D8tWryZkz/Ufk1asLj7dhQ2ck3nPP0hG5bduqRORInElmWdry9x8A2k54MvbaskV9Xz8gz872H5EvvdQ9Ijdu7AzGw4eXjsitW1txJnJNnklmaNry9x8A2k54MvbW/Bb1s7PJ6dO5/9KpXH7++eyceTk7z5/LrplzmZo5l9tmzyfve/S1EVnrwuNt2tQZjPfeu3REbtnSiohcDpdajrc1//cfAEaE8GTsDXWL+itX5qJwsQ11bj52+nRSaz6+wHAX1m/M1d27k9nXz4XikSPdN9S59vHmzSMbkcvhUstmtOXyZW9RAQDDITwhy9yi/sqV669z7Cciu5mc7IzEO+54TTh++kzyn585lz++sjGb3vD6fPffvtcPyCvkUsvBa9vly96iAgCaJzwhmXu/x34j8syZ7uNt3twZkXfeufiZyN275/7MEt41/4vBadullm04k+jyZQDgZssOz1LK7Ul+NskbklxN8mit9UcHNTFYtkuX+o/Is2e7j7dlS2ck3nXX0hE56WxZW7TpUsu2nEl0+TIAcLOVnPG8kuRf1lp/t5SyLcmTpZRP1lr/YEBzg+Tixf4j8ty57uNt3doZiQcPdt9Q51pE3nLL8B4vq6Itl1q25Uyiy5cBgJstOzxrrV9O8uX5j18upTyTZF+SVoZnGy5fa72LFxcPxoWC8uWXu4+3bVtnJB46tHREbto0vMcLA9aWM4ltu3wZAGjeQF7jWUo5kORIkt9Z4GsPJXkoSfbv3z+Iuxu4tly+tqZcuLBkRP75F0/m9J99OVtfPpNdF17O5kuL/HC8Y8f1UNy7N7n77qUjcuPG4T1eWAPaciaxTZcvAwDDUWq39/rrdYBStib5rST/utb62GK3nZ6ersePH1/R/TXhXR/5jQV/mNs3NZlPf/hrVmFGQzYz0/+ZyFdf7T7e1FRe2boj/6fekhc3bcvpzdvz0uT2vLJlR979rrtz5MibOiNy1y4RCT24+UmyZO5M4g9/02FRBwCsCaWUJ2ut0zcfX9EZz1LKhiS/kOTnlorOtawtl6/15Pz5/iPy/Pnu4+3cef1s4xvfmBw+vPiZyF27kg0b8kCXmP/E1sl8+u+OQcxDA5xJBADaaiW72pYkP53kmVrrvx/clIZvTV6+VuvyInJmkVjetet6KO7bl9x778LxeGNErl/eEhmpmIc1pC0bIQEA3GglZzzfleTbkzxdSvns/LF/VWv95ZVPa7ga3wij1rlLU/uNyAsXFh6vlM6I3L8/OXJk8YjcuXPZEbkcazLmAQCAVbGSXW3/Z5IywLmsmr4uX6s1eeWV/iPy4sWF77yU65vl7NmTHDiQvOMdS0fkxERz35ABsKslAABwzYo3F+rHWt1cKEly5kzyxBO9ReSlSwuPsW5dZ0Qu9lrIa79PTfUVkW1625c2zRUAAFi5bpsLCc9rnngieec7r38+MbG8iFy3rrEp2tGSNvHEAwDA+GlkV9uRcvjwXHxeC8rt2xuNyOU4euxER3Qmyczl2Rw9dsIP9Kwp3hsXAIAbra2yWk1btyb335+86U2Nn7lcLjvF0haLPUkCAMD4WXt1RVfddoS1UyxrjSdJAAC4kfBskYcfOJTJDZ0bEdkplrXIkyQAANxIeLbIg0f25Ye/6XD2TU2mJNk3NWljIdYkT5IAAHAjmwu1zINH9glN1ry+3hsXAICRJzyBRniSBACAa1xqCwAAQKOEJwAAAI0SngAAADRKeAIAANAo4QkAAECjSq11eHdWyqkkfzq0O1yePUleWO1JMPasQ9YC65C1wlpkLbAOWQvasA7vqLXuvfngUMOzDUopx2ut06s9D8abdchaYB2yVliLrAXWIWtBm9ehS20BAABolPAEAACgUcLztR5d7QlArEPWBuuQtcJaZC2wDlkLWrsOvcYTAACARjnjCQAAQKOEJwAAAI0ay/Aspby3lHKilPKFUsqHF/h6KaX8x/mvf66U8vbVmCejrYd1+K3z6+9zpZTfLqXcuxrzZPQttRZvuN1fKqXMllK+eZjzYzz0sg5LKe8upXy2lPL7pZTfGvYcGX09/N+8o5Ty30spvze/Dj+wGvNktJVSPlpKeb6U8vkuX29lq4xdeJZSJpL8pyRfl+TNSd5fSnnzTTf7uiQH5389lOQnhjpJRl6P6/BPknx1rfWtSX4oLX4xOWtXj2vx2u3+TZJjw50h46CXdVhKmUry40m+sdZ6T5K/N/SJMtJ6/Pfwnyb5g1rrvUneneTflVI2DnWijIOPJXnvIl9vZauMXXgmuS/JF2qtX6y1Xkry80ned9Nt3pfkZ+uczySZKqW8cdgTZaQtuQ5rrb9daz09/+lnktw25DkyHnr5NzFJvivJLyR5fpiTY2z0sg7/fpLHaq1fSpJaq7XIoPWyDmuSbaWUkmRrkpeSXBnuNBl1tdZPZW5tddPKVhnH8NyX5M9u+PzZ+WP93gZWot819sEkv9LojBhXS67FUsq+JH8nyU8OcV6Ml17+TfzKJDtLKb9ZSnmylPIdQ5sd46KXdfhjSe5O8lySp5N8qNZ6dTjTg7/QylZZv9oTWAVlgWM3v6dML7eBleh5jZVS/kbmwvOvNjojxlUva/FHknxfrXV27kl+GLhe1uH6JO9I8p4kk0meKKV8ptb6R01PjrHRyzp8IMlnk3xNkruSfLKU8j9qreeanhzcoJWtMo7h+WyS22/4/LbMPWvV721gJXpaY6WUtyb5qSRfV2t9cUhzY7z0shank/z8fHTuSfL1pZQrtdbHhzNFxkCv/ze/UGt9NcmrpZRPJbk3ifBkUHpZhx9I8pFaa03yhVLKnyT5qiT/azhThCQtbZVxvNT2fyc5WEq5c/7F4N+S5Jduus0vJfmO+R2j7k9yttb65WFPlJG25DospexP8liSb/eMPg1aci3WWu+stR6otR5I8okk/0R0MmC9/N/8i0n+WillfSllc5K/nOSZIc+T0dbLOvxS5s66p5Ty+iSHknxxqLOElrbK2J3xrLVeKaX8s8ztzDiR5KO11t8vpfzj+a//ZJJfTvL1Sb6Q5Hzmnt2CgelxHf5Akt1Jfnz+TNOVWuv0as2Z0dTjWoRG9bIOa63PlFJ+NcnnklxN8lO11gXfagCWo8d/D38oycdKKU9n7nLH76u1vrBqk2YklVI+nrldk/eUUp5N8oNJNiTtbpUyd6UAAAAANGMcL7UFAABgiIQnAAAAjRKeAAAANEp4AgAA0CjhCQAAQKOEJwAAAI0SngAAADTq/wMQd8/bwSZXkAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "x = np.array([0, 1]).reshape(-1, 1)\n",
    "plt.figure(figsize=(16, 4))\n",
    "plt.scatter(X, y)\n",
    "plt.plot(x, LinearRegression().fit(X, y).predict(x), 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By changing the loss function to the mean absolute deviation \n",
    "$$ \\sum_i \\left|y_i-\\textrm{model}(x_i)\\right|\\enspace, $$\n",
    "\n",
    "we can let the model put the same focus on each error. This yields the least absolute deviation (LAD) regression that tries to agree with the majority of the points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x266ac8b6b88>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA54AAAD4CAYAAACEyjk9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df3Dkd33f8ddb+0Na7erH6YfvfLqfPhthjH0+o7ouTgqFMmcowS6lM6H5QV06TqctYZrigpOZkBmageT6i05CqEscwpRx2hDXIR3gQmEIDcGkZ2w4gjlwDTa+c3y6H7rTSitpf7z7x0p3Wt3+3v2u9rv7fMzc3Oqrve9+pPla1nO/n+/na+4uAAAAAACCMrDdAwAAAAAA9DbCEwAAAAAQKMITAAAAABAowhMAAAAAECjCEwAAAAAQqGgnX2xqasoPHDjQyZcEAAAAAHTIk08+ec7dp7du72h4HjhwQCdOnOjkSwIAAAAAOsTMni+3nam2AAAAAIBAEZ4AAAAAgEARngAAAACAQBGeAAAAAIBAEZ4AAAAAgEB1dFVbAAAAoFs8/tRpHTt+SmcWMto9ntCDR2d135GZ7R4W0JMITwAAAPSdx586rYceO6lMNi9JOr2Q0UOPnZQk4hMIAFNtAQAA0HeOHT91JTo3ZLJ5HTt+aptGBPQ2whMAAAB958xCpqHtAFpDeAIAAKDv7B5PNLQdQGsITwAAAPSdB4/OKhGLlGxLxCJ68OjsNo0I6G0sLgQAAIC+s7GAEKvaAp1BeAIAAKAv3XdkhtAEOoSptgAAAACAQBGeAAAAAIBAEZ4AAAAAgEARngAAAACAQBGeAAAAAIBAEZ4AAAAAgEARngAAAACAQNUMTzN7xMzOmtl3ynzufWbmZjYVzPAAAAAAAGFXzxnPT0q6Z+tGM9sr6U2SXmjzmAAAAAAAPaRmeLr7VyVdKPOp/yjp30jydg8KAAAAANA7mrrG08zeJum0u3+rjuc+YGYnzOzE/Px8My8HAAAAAAixhsPTzIYl/YqkX63n+e7+sLvPufvc9PR0oy8HAAAAAAi5Zs54HpJ0UNK3zOxHkvZI+qaZ7WrnwAAAAAAAvSHa6D9w95OSrtv4eD0+59z9XBvHBQAAAADoEfXcTuVRSV+XNGtmL5rZu4MfFgAAAACgV9Q84+nu76zx+QNtGw0AAAAAoOc0taotAAAAAAD1IjwBAAAAAIEiPAEAAAAAgSI8AQAAAACBIjwBAAAAAIEiPAEAAAAAgSI8AQAAAACBIjwBAAAAAIEiPAEAAAAAgSI8AQAAAACBIjwBAAAAAIEiPAEAAAAAgSI8AQAAAACBIjwBAAAAAIEiPAEAAAAAgSI8AQAAAACBIjwBAAAAAIEiPAEAAAAAgSI8AQAAAACBqhmeZvaImZ01s+9s2nbMzL5nZt82s/9pZuPBDhMAAAAAEFb1nPH8pKR7tmz7oqRXu/ttkr4v6aE2jwsAAAAA0CNqhqe7f1XShS3b/tTdc+sfPiFpTwBjAwAAAAD0gHZc4/lPJH2+DfsBAAAAAPSglsLTzH5FUk7Sp6s85wEzO2FmJ+bn51t5OQAAAABACDUdnmb2LklvlfQz7u6VnufuD7v7nLvPTU9PN/tyAAAAAICQijbzj8zsHknvl/Q6d19u75AAAAAAAL2kntupPCrp65JmzexFM3u3pN+SNCLpi2b2tJl9POBxAgAAAABCquYZT3d/Z5nNvxvAWAAAAAAAPagdq9oCAAAAAFAR4QkAAAAACBThCQAAAAAIFOEJAAAAAAgU4QkAAAAACBThCQAAAAAIFOEJAAAAAAgU4QkAAAAACBThCQAAAAAIFOEJAAAAAAgU4QkAAAAACBThCQAAAAAIFOEJAAAAAAgU4QkAAAAACBThCQAAAAAIFOEJAAAAAAgU4QkAAAAACBThCQAAAAAIFOEJAAAAAAgU4QkAAAAACFTN8DSzR8zsrJl9Z9O2CTP7opn9YP3vHcEOEwAAAAAQVvWc8fykpHu2bPuApC+5+02SvrT+MQAAAAAA16gZnu7+VUkXtmy+V9Lvrz/+fUn3tXlcAAAAAIAe0ew1njvd/SVJWv/7ukpPNLMHzOyEmZ2Yn59v8uUAAAAAAGEV+OJC7v6wu8+5+9z09HTQLwcAAAAA6DLNhufLZna9JK3/fbZ9QwIAAAAA9JJmw/Ozkt61/vhdkv64PcMBAAAAAPSaem6n8qikr0uaNbMXzezdkj4i6U1m9gNJb1r/GAAAAACAa0RrPcHd31nhU29s81gAAAAAAD0o8MWFAAAAAAD9jfAEAAAAAASK8AQAAAAABIrwBAAAAAAEivAEAAAAAASK8AQAAAAABIrwBAAAAAAEivAEAAAAAASK8AQAAAAABIrwBAAAAAAEivAEAAAAAASK8AQAAAAABIrwBAAAAAAEivAEAAAAAASK8AQAAAAABIrwBAAAAAAEivAEAAAAAASK8AQAAAAABIrwBAAAAAAEivAEAAAAAAQq2so/NrN/JemfSnJJJyXd7+4r7RgYAAAo9fhTp3Xs+CmdWcho93hCDx6d1X1HZrZ7WAAA1NT0GU8zm5H0i5Lm3P3VkiKSfrpdAwMAAFc9/tRpPfTYSZ1eyMglnV7I6KHHTurxp05v99AAAKip1am2UUkJM4tKGpZ0pvUhAQCArY4dP6VMNl+yLZPN69jxU9s0IgAA6td0eLr7aUn/TtILkl6SdMnd/3Tr88zsATM7YWYn5ufnmx8pAAB97MxCpqHtAAB0k1am2u6QdK+kg5J2S0qa2c9ufZ67P+zuc+4+Nz093fxIAQDoY7vHEw1tBwCgm7Qy1fbvSvqhu8+7e1bSY5Je255hAQCAzR48OqtELFKyLRGL6MGjs9s0IgAA6tfKqrYvSLrLzIYlZSS9UdKJtowKAACU2Fi9llVtAQBh1HR4uvs3zOwzkr4pKSfpKUkPt2tgAACg1H1HZghNAEAotXQfT3f/oKQPtmksAAAAAIAe1OrtVAAAAAAAqIrwBAAAAAAEivAEAAAAAASK8AQAAAAABIrwBAAAAAAEivAEAAAAAASK8AQAAAAABIrwBAAAAAAEivAEAAAAAASK8AQAAAAABIrwBAAAAAAEivAEAAAAAASK8AQAAAAABIrwBAAAAAAEivAEAAAAAASK8AQAAAAABIrwBAAAAAAEivAEAAAAAASK8AQAAAAABIrwBAAAAAAEqqXwNLNxM/uMmX3PzJ4xs7/VroEBAAAAAHpDtMV//1FJX3D3d5hZXNJwG8YEAAAAAOghTYenmY1K+tuS/rEkufuapLX2DAsAAAAA0CtamWp7g6R5Sb9nZk+Z2SfMLLn1SWb2gJmdMLMT8/PzLbwcAAAAACCMWgnPqKQ7JP2Oux+RtCTpA1uf5O4Pu/ucu89NT0+38HIAAAAAgDBq5RrPFyW96O7fWP/4MyoTngAAAACABiwvS+fOSefPl/6dz0vvfe92j64pTYenu/+1mf3YzGbd/ZSkN0r6bvuGBgAAAAAh5l45Iis9Pn9eymTK7298vP/Cc917JH16fUXb5yTd3/qQAAAAAKDLuEtLS41H5MpK5X1OTEiTk9LUlLR3r3T77cXHG9s2P56cLD4/pFoKT3d/WtJcm8YCAAAAAMFzl9LpxiNydbX8/sxKI3L/fumOO6pH5I4dUrTV84Dh0T9fKQAAAIDe4y4tLjYekWsV7gRpVgzDjUg8eFCam6sckVNTxSmwkUhnv+6QITwBAAAAdAd36fLl6sFYbls2W35/AwOlEXnokHTnnbUjcqCVm3+gHMITAAAAQPu5S5cuNR6RuVz5/UUipRF5003SXXeVn8a6sW1sjIjsEoQnAAAAgOoKhcYj8sKF6hG5ORhnZ6XXvrZ6RI6OEpEhRngCAAAA/aRQkBYWGo/IfL78/qLR0mC8+ebKC+psjkizzn7d2FaEJwAAABBWhYJ08WLjEVkolN9fLFYajLfcUjsiR0aISNREeAIAAADdIJ9vLiLdy+8vHi8NxltvrR2RqRQRiUAQngAAAEC75fPFKGwkIi9erByRg4OlwXj4cOVVWTceJ5NEJLoG4QkAAABUk8s1F5GVDA2VRuK+fbUjcniYiESoEZ4AAADoH9ls4xG5sFB5f4lEaSQeOFB5GuvmiAT6DOEJAACAcMpmqwdjuaC8dKny/oaHSyPxhhuqR+TkJBEJ1InwBAAAwPZbW2s8Ii9frry/ZLI0Em+8sXZEJhKd+3qBPkN4AgAAoL1WVxuPyMXFyvsbGSkNxle8onZEDg117usFUBPhCQAAgMpWVhqPyHS68v5GR0tD8ZWvrB2Rg4Od+3oBBILwBAAA6BeZTP0L6mz8vbRUeX9jY1dD8brrpFe9qvKCOhsRGY937usF0DUITwAAgDBaXm48IpeXK+9vfPxqJO7aJb361dUjcmKCiARQN8ITAABgO7k3F5GZTOV97thxNRJ375Zuu612RMZinfua0bDHnzqtY8dP6cxCRrvHE3rw6KzuOzKz3cMC6kZ4AgAAtIt7cWpqoxG5slJ5nxMTVyNx717p9ttrR2SUX/F6yeNPndZDj51UJpuXJJ1eyOihx05KEvGJ0OCnEgAAQDnuxUVyGo3I1dXy+zMrjcj9+6U77qgekTt2EJHQseOnrkTnhkw2r2PHTxGeCA1+kgEAgN7nXrxdR6MRubZWfn9mVxfLmZqSDh6U5uZqR2Qk0tmvGz3hzEL5adWVtgPdqOXwNLOIpBOSTrv7W1sfEgAAQBXu0uXLjUdkNlt+fwMDpRF56JB0552VI3JqqrgQz8BAZ79u9K3d4wmdLhOZu8cT2zCa3tFt181m81ktZZeUXktraa34d3otfWVbei0tk+n+I/dv2xhb0Y4znu+V9Iyk0TbsCwAA9BN36dKlxiMylyu/v0ikNCJvukm6667qETk2RkSiqz14dLbkGk9JSsQievDo7DaOKtxauW42V8hVDMOq27PVn7OWrzDDYpOJxER/hqeZ7ZH09yT9uqRfasuIAABAOBUKjUfk+fPVI3JzJM7OVp7GurFtdJSIRM/ZCKFuOjsXBhuBWC4M3/+5J3S+cFmFSEZuK3Kt6IJW9Av/K6v/8dxY1ahczVe4jruMARtQKp668icZSyoVT2lqeEr7x/cXt8dSSsaT1zwnFS+/PazM3Zv/x2afkfRhSSOS3lduqq2ZPSDpAUnat2/fa55//vmmXw8AAHRIoSAtLNQfkefOSRcuSPl8+f1Fo9WDsVJEmnX26wbQcflCvuYU04rbqzxnJVdlteitfECmIQ1oSDdMTlaMvnrCcPP2wcigrM9+jpnZk+4+t3V702c8zeytks66+5Nm9vpKz3P3hyU9LElzc3PNVy4AAGhOoSBdvNh4RBYK5fcXi5UG4y231I7IkREiEgi5jUAsF4GtTDFtJBAHbKBs4E0kJrR3dG/dcfgLn/qOzi2azIdkGpIpLpNpZjyhr73nDQF+F/tXK1Nt75b0NjN7i6QhSaNm9t/c/WfbMzQAAHCNfL7xiLx4sXJExuOlwXjrrbUjMpUiIoEuVvDClcir+/rDtbTS2erPyeTqX0XXZGWjb3xoXHtG9xS3x+o7a7h5+1B0qC1nEH/tzXu4brbDmg5Pd39I0kOStH7G831EJwAADcjni2cWG43ISpfJDA6WBuPhw5UX1Nl4nEwSkcA2KXhBy9nl+sOwzmmmjQbiRuBtjr6xoTHNjM40Nc30y9+9pI/+7+f10qWV4vWob+i+61HDdt1st63A2wzu4wkAQDvkcpUjslJQLixUjsihodJIPHKkdkQODxORQAA2ArGd1x+m19Jazi43NI5y0Tc6OKrrU9c3df1hKp5SIppo6zWIjz91Wr/22WebWi220+47MtN1YyqnlRV4u0lbwtPdvyLpK+3YFwAA2y6bbS4iK0kkSiNx//76IjKkeuGdebRmu44Bd79yBrGdU0yXsksNjSMZS5aNvl2pXU1PMU3EEhqw7l+x+djxUyXTVyUpk83r2PFT/BxoUq98TznjCQDobdns1dt21BuRly5V3t/wcGkkHjxYPSInJ0MdkY3qlXfm0bx6joGNQGwoDiucTdz8+eXsslz1r2U5HBsuG307UzubXsl0ODYcikAMypmF8tN8K21Hbb3yPSU8AQDhsbZ2bSjWisjLlyvvL5ksjcRDh6rf7mNysnj2EhWF7Z15zs5W5+7K5DINheGnnjiltBbl8VUVlJFrRQVb0Ts/u6qxL11dFbWRQExEE2Wjb3p4uukppv0eiEHZPZ7Q6TJBtHucn53N6pXvKeEJANgeq6uNR+TiYuX9pVKlkXjTTbUjcmioc19vnwjTO/O9dHbW3bWSW2nsjGEdt7loJhBX8zHZQEIDKt6mYsCHFPUxDeQGde/sbMNTTIdjw4oMRAL87qGdHjw6y2qxbdYr31PCEwDQupWVxiMyna68v9HR0mCcna0dkYODnft6UVGY3pnfjrOzG4HYjimmm5+zlF1SwSvcMqeMoehQ2fCbHJ4sPo7VjsOt25KxpCIDEd39kS+XPQZmxhP6Lz/F/RF7XdhWiw2DXvmeEp4AgFIrK/Xd2mPztqUqC2+MjV0Nxelp6eaba0dkPN65rxdtFaZ35qudnXV3reZX23qLi43tjQbiRuAVCkO6mDblcnElYim9+vr9euWB6xqeZroRiEEJ0zGAYIRltdgw6YXvKeEJAL0sk2k8IperLO8/Pn41FHfulG65pXpETkwQkX2m0+/MbwRiM3G4mHxOS9kluVbktqKCVuSWkWxVsQ+tKO/52gNYNxgZLBt9e0b3NHX94caqqNGB4q9qG9OCBzZibk2afz6iX3rNrV33y2ivnJ0B0F7mle4fFoC5uTk/ceJEx14PAHrK8nJ9Ebn5cabKdXU7dlQOxq2PNyIyFuvc14ue4u5ay6+1/TYX6bV0Q4EYj8SvxF0hP6hziyYvDGlAgzJPKDYwrJ+8cUaHZ3Y1FIwbgRiUatNXv/YBpq8C6B5m9qS7z23dzhlPAOg09/ojcvO2lZXK+5yYuBqJe/ZIt99eOyKj/C8A5W0EYrunmeYKubrHEI/Ey0bf7pHdTd3iYmNbLFL65klYVrUN06JNAFAOv3UAQCvci9c3NhqRq6vl92dWGpH79klHjlSPyB07iMg+tZZfazwM61jJtJFAjA5ENRIfuSb6rk9d39IU03ikM1O0w3LdVJgWbQKAcvhNBQA2uBdXWm10OuvaWvn9mV1dLGdqSjpwQHrNa2pHZITbBvSabD5b9ab3zZ5JzBaydY8hOhAtG307Uzt1KH6orlVMy23vVCD2OxbsARB2hCeA3uRevOdjoxGZrfCL/MBAaUQeOiTdeWf1iBwfJyJDJpvPtv02F40GYsQiGhkcuSbwrktepxt23HBle6NnEgnEcGPBHgBhx+JCgMJzjU/fcpcuX24sIs+frxyRkUhpRFZblXXj8dhYMT7RFXKFXNvCcH7pkhYyi8p5RrL6p5hGLFI7/JqYZhqPxGVmAX73AAAIDosLARVsLFG/MX3p9EJGDz12UpK6Lj57IpDdpUuXGo/IXIUgiERKI3F2tnZEjo4SkR2SK+RKbnDfrpVMV/MVrpEtY8AGroTd5uibGp5S1K/TmaU1JQqDMh+SaUiDA8P6B3cc0k8e2ls1GAlEAADqR3ii7x07fqrkmhlJymTzOnb8VFdFXVcGcqFQPiKrBeX581K+wq0PotHSSLz55tpnJUdHi9dSoiX5Qj6QKabNBOLWs4GTw5PaH9/f9JnEwchgxUC8+yNf1vjqtQu2fOt7CX3sPm5RAQBAuxCe6HthWaI+8EAuFKSFhcYjslAov79YrDQSb7mldkSOjGxLRIbpTPJGILbzFhfptbRWclVu1bKFycpG346hHdo7urfpKaZD0aGOn0EMy3//AACEHeGJvheWJeob+gU5n288Ii9cqByR8XhpMN56a+2ITKVCcSYyqDPJ+UJey9nltt/motFALBd940PjmhmdKW5vYiXT7QjEoITlv38AAMKO8ETf6/ol6vN56eJF3bU2r+zZs9qRWdSO5cuayFzWeOay9uSXpXsfvjYiKy0cNjhYGoyHD9eOyGQyFBHZjN/8wjNayi7KtaqCZeRa0Up+Rb/8+W8pP3hj02cTM7nGzpiVi76xoTHNjM40PcU0EU30TCAGpev/+wcAoEcQnuh7HV2iPpcrRmG1BXW2brt4UXLXo2V2txKNqzA5KeV3FkPxyJHKC+psPB4eDmVEFrxw5QxiO6eZLmeXpTInt17OSu/4w2u3X7mVxaboG4mP6PrU9U1PMU3EEhqw3lrsKCzTl7lFBQAAncHtVIBm5XJXr3NsJCIrSSRq3trjawvSf33msn6Qi2tw10794k8d7rpfkAteUCabafsU0+XsckPjGI4NV4++9Smm//0v55VeiWpgfUXTjb+vS43r0+9+Xcm/7cVADMLW6ctS8Szih99+a9cdrwAAoL24nQpQTTbbeEQuLFTe3/BwaUQePFj9TOTkZPHf1HD3+p92cHctZ5fbvpLpUnapoXFsBOLWONyZ2tn0PRGHY8N1B+JPTK9HUq40kj705lt1287ui6QwnEkMy0rRAACgc5oOTzPbK+lTknZJKkh62N0/2q6BAU1bW2s8Ii9dqry/ZLI0Eg8dqh2RifYtTOLuyuQybY3DjTOIrvpnPCSiibLRNz083dT1h6l4qqFADEqYplp25S11ymClWAAAsFUrZzxzkv61u3/TzEYkPWlmX3T377ZpbIC0utp4RF6+XHl/qVRpJN50U9WprZqclIaG6hrqRiAuraWVvvjXbbvNxdLaUsOBWC76poeni9tjjV1/uBGIkYFI3WMIm/uOzHRVuFUSljOJrBQLAAC2ajo83f0lSS+tP140s2ckzUgKZXiGYfpa6K2uVg/GckG5uFh5fyMjpZE4O1s7IgcH5e5aya1Uib6XtJR9VukX0ko/u2l7tvbZxkYCcSg6VDb6Jocnm55imowlezoQ+11YziSyUiwAANiqLdd4mtkBSUckfaPM5x6Q9IAk7du3rx0v13Zhmb7WVVZWakbky8+d1sUfv6TU4oImVhY1vFbll+OxsauhOD0t3XyzfGpSqxNjSk+OKD2W0NJYQumRQaWHo1pKRJQurFY4S/is0mtPF7f/KK309699TsEr3K+yjMHIYNnom0hMtDTFNDrAJdZoTFjOJIZp+jIAAOiMlle1NbOUpD+T9Ovu/li153brqrZ3f+TLZX+ZmxlP6GsfeMM2jKjDMpnGz0QuXV1AxiWtRqV0XFqKSenJEZ0dS+rZaFQvJ4Z0LhnXhURcl4ejmjm4QyM7k1efG8kprTUt5cqvgtpoIFaNviammCbjSQIRXYPVYgEAQLcLZFVbM4tJ+iNJn64Vnd0sLNPX6rK8XDEi/fw5rZ0/q/TFl5W+fF5Li+eVXrqotK8WIzC+Ho/xTY9Tg0qPxJU+FNPSLRGlh0zpWFJLkYTSA+vR6KvK++brzhbX/5QXvxgvG30zozMtTTGNRWKBf3uB7cSZRAAAEFatrGprkn5X0jPu/h/aN6TO67bpa+6utdyqli7NK/3yj5U+d0ZL518qBuPCWS1dXg/GpYtaWrm8fv3hkpbyK0pH8lfPJm4NyQkpP1X/OOKRuJKb74W4Hn4zdcThez79XZmGZJ7QgIZkPqQBFf/86CP3BvfNA3pcWBZCAgAA2KyVM553S/o5SSfN7On1bb/s7p9rfVid1cpCGGv5tTpucbGodPqili6fU3r9LONS5pLSK4vFBWtyGaULK8Wzh5ZTOppXrtodJgYkjaz/kRR100ghpqQNKxUZUiqaVDKe1PVDo0oNjyuZ3KHUyIRSg6N1TzNNxpOKR+JNf08/Olp5+jIAAACA/tLKqrZ/LsnaOJZtc9+RGV1YeVn/9kuP6OLKopJDWR3al9CfvPCHevT/lYnKlUWl1xa1lF1W1nN1v040L6XWin+S2auPdyqmQzakVHSsGIODI0olxpRMjiuVmlRqbErJsWmlJnYpNbVbycldSg2NXQnGVgIxKKxqCQAAAGADq6asu2XHBf0w/9tSTFrMD+jPn4spVYgqlRtQcs2VWnHtXM4puZxVaqVQjMe1LSGZNaUGR5QcHlMqNaHUyJSSY1NK7dip1MQuxad3XXu7j/FxKVL/7S8ef+q0jv3eKZ1ZeL6rr+/iWjQAAAAAGwjPdXecH9S53yhGZDxfkEVy0uTYtfeCPFjm/pCbI3Kg2hzZ1oTtti9ci9bfuDcuAAAANhCe62KHj2jyy1+/GpSjo4FGZDOOHT9VMnVVkjLZvI4dP8Uv9OgqYXuTBAAAAMHqrrLaTqmUdNdd0o03Bn7mslk9ddsX9LRqb5IAAACg/3RfXaGiSrd32a7bvgCV8CYJAAAANiM8Q+TBo7NKxEoXImKlWHQj3iQBAADAZoRniNx3ZEYffvutmhlPyFS8J+aH334r18yh6/AmCQAAADZjcaGQYaVYhAG30wEAAMBmhCeAQPAmCQAAADYw1RYAAAAAECjCEwAAAAAQKMITAAAAABAowhMAAAAAECjCEwAAAAAQKHP3zr2Y2byk5zv2gs2ZknRuuweBvsdxiG7AcYhuwbGIbsBxiG4QhuNwv7tPb93Y0fAMAzM74e5z2z0O9DeOQ3QDjkN0C45FdAOOQ3SDMB+HTLUFAAAAAASK8AQAAAAABIrwvNbD2z0AQByH6A4ch+gWHIvoBhyH6AahPQ65xhMAAAAAECjOeAIAAAAAAkV4AgAAAAAC1ZfhaWb3mNkpM3vWzD5Q5vNmZv95/fPfNrM7tmOc6G11HIc/s378fdvM/sLMDm/HONH7ah2Lm573N8wsb2bv6OT40B/qOQ7N7PVm9rSZ/ZWZ/Vmnx4jeV8f/m8fM7E/M7Fvrx+H92zFO9DYze8TMzprZdyp8PpSt0nfhaWYRSb8t6c2SXiXpnWb2qi1Pe7Okm9b/PCDpdzo6SPS8Oo/DH0p6nbvfJulDCvHF5OhedR6LG8/7DUnHOztC9IN6jkMzG5f0MUlvc/dbJP3Djg8UPa3On4f/QtJ33f2wpNdL+vdmFu/oQNEPPinpniqfD2Wr9F14SrpT0rPu/py7r0n6A0n3bnnOvZI+5UVPSBo3s+s7PVD0tJrHobv/hX7f4nMAAAKsSURBVLtfXP/wCUl7OjxG9Id6fiZK0nsk/ZGks50cHPpGPcfhP5L0mLu/IEnuzrGIdqvnOHRJI2ZmklKSLkjKdXaY6HXu/lUVj61KQtkq/RieM5J+vOnjF9e3NfocoBWNHmPvlvT5QEeEflXzWDSzGUl/X9LHOzgu9Jd6fia+QtIOM/uKmT1pZj/fsdGhX9RzHP6WpJslnZF0UtJ73b3QmeEBV4SyVaLbPYBtYGW2bb2nTD3PAVpR9zFmZn9HxfD8iUBHhH5Vz7H4nyS9393zxTf5gbar5ziMSnqNpDdKSkj6upk94e7fD3pw6Bv1HIdHJT0t6Q2SDkn6opn9H3e/HPTggE1C2Sr9GJ4vStq76eM9Kr5r1ehzgFbUdYyZ2W2SPiHpze5+vkNjQ3+p51ick/QH69E5JektZpZz98c7M0T0gXr/33zO3ZckLZnZVyUdlkR4ol3qOQ7vl/QRd3dJz5rZDyW9UtJfdmaIgKSQtko/TrX9v5JuMrOD6xeD/7Skz255zmcl/fz6ilF3Sbrk7i91eqDoaTWPQzPbJ+kxST/HO/oIUM1j0d0PuvsBdz8g6TOS/jnRiTar5//NfyzpJ80sambDkv6mpGc6PE70tnqOwxdUPOsuM9spaVbScx0dJRDSVum7M57unjOzf6niyowRSY+4+1+Z2T9b//zHJX1O0lskPStpWcV3t4C2qfM4/FVJk5I+tn6mKefuc9s1ZvSmOo9FIFD1HIfu/oyZfUHStyUVJH3C3cveagBoRp0/Dz8k6ZNmdlLF6Y7vd/dz2zZo9CQze1TFVZOnzOxFSR+UFJPC3SpWnCkAAAAAAEAw+nGqLQAAAACggwhPAAAAAECgCE8AAAAAQKAITwAAAABAoAhPAAAAAECgCE8AAAAAQKAITwAAAABAoP4/srXOgQKqlNEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "from sklego.linear_model import LADRegression\n",
    "\n",
    "x = np.array([0, 1]).reshape(-1, 1)\n",
    "plt.figure(figsize=(16, 4))\n",
    "plt.scatter(X, y)\n",
    "plt.plot(x, LinearRegression().fit(X, y).predict(x), 'r')\n",
    "plt.plot(x, LADRegression().fit(X, y).predict(x), 'g')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### See also\n",
    "scikit-learn tackles this problem by offering a variety of [robust regressors](https://scikit-learn.org/stable/auto_examples/linear_model/plot_robust_fit.html). Many of them use an indirect approach to reduce the effect of outliers. [RANSAC](https://en.wikipedia.org/wiki/Random_sample_consensus), for example, samples random points from the dataset until it consists of all inliers.\n",
    "The closest thing to LADRegression that scikit-learn offers is the [HuberRegressor](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.HuberRegressor.html) with a loss function that is partly a squared and partly an absolute error. However, it is more complicated and requires hyperparameter tuning to unleash its full potential.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QuantileRegression\n",
    "\n",
    "This is an extension of the LADRegression (see above). While the LADRegression outputs a line that over- and underestimates around 50% of the data, the Quantileregression yields lines with different over- and underestimation shares. This can be used for creating simple confidence intervals around predictions. As an example, consider the following:\n",
    "1. Create a QuantileRegression with quantile=0.1,\n",
    "2. create a QuantileRegression with quantile=0.9,\n",
    "\n",
    "and around 80% of the data is between these two lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PolyCollection at 0x255934fab88>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABINUlEQVR4nO2dd5wdV3n3v2dmbt2uXXVZsmXJvcrCNjbNNBeKaaGEYgiJ4aUECElewkteWpI3IZAAgZA4xtSA6UGmGGyDsTFucpctyZLVV217uf3OnPePM3Pb3m3au0Wa5/v5XO1q7t07M7f8zjPP+Z3nUVprBEEQhHBhzfcBCIIgCHOPiL8gCEIIEfEXBEEIISL+giAIIUTEXxAEIYQ4830AU6Grq0uffPLJ830YgiAIxxUPPfRQr9Z6cb37jgvxP/nkk9m8efN8H4YgCMJxhVJq73j3SdpHEAQhhIj4C4IghBARf0EQhBAi4i8IghBCRPwFQRBCiIi/IAhCCBHxFwRBCCEi/oIghIrH9g/yxIGh+T6MeUfEXxCEUPH3v9jKP966db4PY945Llb4CoIgNIpM3kVHJO6VV0AQhFBRcD2KnnQwFPEXBCFU5Iseroi/iL8gCOEiV/QouiL+Iv6CIISKvOvhaRF/EX9BEEKF5PwNIv6CIIQKyfkbRPwFQQgV+aJH0fPm+zDmHRF/QRBCg+dpip7GlQlfEX9BEMJD3jURv+T8RfwFQQgRgfhLzl/EXxCEEJEvSuQfIOIvCEJoCMRfIn8Rf0EQQkShlPMXt4+IvyAIoUEi/zIi/oIghIac5PxLiPgLghAaAreP1sbzH2ZE/AVBCA2FYjnXH/boX8RfEITQEET+IHl/EX9BEEJDviryD7fjR8RfEITQUCn+EvkLgiCEhMq0j+T8BUEQQoJE/mVE/AVBCA0S+ZcR8RcEITRURf4hr+kv4i8IQmgQt08ZEX9BEEJDQXz+JUT8BUEIDXlZ4VtCxF8QhNCQk8i/hIi/IAihQSL/MlMWf6XUTUqpo0qpLRXbFimlblNK7fB/dvjblVLqi0qpnUqpx5VSGyr+5jr/8TuUUtc19nQEQRDGpzrnLxO+U+XrwFU12z4C3KG1Xg/c4f8f4GpgvX+7HvgKmMEC+DhwCXAx8PFgwBAEQZhtqiJ/sXpODa31XUB/zeZrgW/4v38DeFXF9m9qw31Au1JqOXAlcJvWul9rPQDcxtgBRRAEYVaQFb5lZprzX6q1PuT/fhhY6v++Ethf8bgD/rbxtguCIMw6ssK3TMMmfLXWGmjYq6mUul4ptVkptbmnp6dRTysIQojJF8sSJZH/zDjip3Pwfx71t3cDJ1U8bpW/bbztY9Ba36C13qi13rh48eIZHqYgCIKJ/G1LARL5z1T8NwGBY+c64KcV29/mu34uBYb89NCvgJcqpTr8id6X+tsEQRBmnXzRJRmxAXH7OFN9oFLqu8ALgC6l1AGMa+cfge8rpd4J7AVe7z/8F8A1wE4gDbwDQGvdr5T6NPCg/7hPaa1rJ5EFQRBmhXzRIxG1GckVQx/5T1n8tdZvGueuF9V5rAbeO87z3ATcNNX9CoIgNIq865GMBpF/uMVfVvgKghAaCkVN3E/7iM9fEAQhJEjkX0bEXxCE0JAveiSjJtsd9py/iL8gCKEh50/4grh9RPwFQQgNhYq0j0T+giAIISFf9EhEJOcPIv6CIISIvFtO+0jkLwiCEAJcT+N6Wtw+PiL+giCEgqCRS8ntIz5/QRCEE5+cX8s/5hjZE7ePIAhCCMhXiL9jKcn5z/cBCIIgzAVBI5eoY2FbSnL+830AgiAIc0GhWBZ/ifxF/AVBCAlB5B+xJfIHEX9BEEJCkPOP2haObYn4z/cBCIIgzAW5YnXOX9I+giAIIaDgVuf8xeopCIIQAirTPhL5i/gLghAS8sXayF/EXxAE4YSn1ucvkb8gCEIIqHL7WBau1PYRBEE48an1+UvkLwiCEAKqavvY4vYR8RcEIRTkxedfhYi/IAihID/G5y/iLwiCcMJTEJ9/FSL+giCEgrzroRTYljJuHxF/QQgnR4ezDKTy830YwhyRL3pEbQullET+iPgLIeY9//0wn7jlyfk+DGGOyBU9on4LR6ntA858H4AgzBeHhrI4tprvwxDmiILrlfr32paSBu7zfQCCMF8MZwuM5orzfRjCHJEvekRsP/K3xe0j4i+EEs/TjOaKjGZF/MNC3i2nfWyZ8BXxF8JJKl9EayTyDxHBhC8gPXwR8RdCyrAf8Yv4h4dCVeQvaR8RfyGUjGQLAGQLXqnDk3Bik6vM+VuKYsjdPiL+QigZqcj1pyT6DwX5okT+lYj4C6FkOFMo/S6pn3CQr7B6Ss5fxF8IKZWRv4h/OCi45QlfW5q5iPgL4STI+QNi9wwJtT5/ifwbgFJqj1LqCaXUo0qpzf62RUqp25RSO/yfHf52pZT6olJqp1LqcaXUhkYcgyBMh+EKwR+RyD8USM6/mkZG/ldorS/QWm/0//8R4A6t9XrgDv//AFcD6/3b9cBXGngMgjAlhisif5nwDQf5mto+4vaZPa4FvuH//g3gVRXbv6kN9wHtSqnls3gcgjCGqpy/pH1CQb7G5+9ps9I7rDRK/DXwa6XUQ0qp6/1tS7XWh/zfDwNL/d9XAvsr/vaAv00Q5ozhTIGu5iggE75hoXaFL4Crwyv+jarq+RytdbdSaglwm1JqW+WdWmutlJrWq+wPItcDrF69ukGHKQiGkWyRZW1xekfzVVcBwolLbW0fANfTROz5PKr5oyGRv9a62/95FPgJcDFwJEjn+D+P+g/vBk6q+PNV/rba57xBa71Ra71x8eLFjThMQSgxki3QnoiSjNqS8w8J9SL/MDt+Ziz+SqkmpVRL8DvwUmALsAm4zn/YdcBP/d83AW/zXT+XAkMV6SFBmBOGs0VaEw7NMUfSPiGg6Hp4mqqcPxBqr38j0j5LgZ8opYLn+47W+lal1IPA95VS7wT2Aq/3H/8L4BpgJ5AG3tGAYxCEaTGSLdASi9Acc8TqGQIKvsiX3D52EPmH1/EzY/HXWu8Czq+zvQ94UZ3tGnjvTPcrCDNhJFukJe7QHHck7RMC8kUj8hG7JvKXtI8ghIeC65HOu7QmTOQvVs8Tn5zrAlT5/GFh5vz/3y+28vnbn571/UgPXyF0BGLfEndoijn0p9LzfETCbBNE/jF7rNtnobHpsYMMpPO845JltEXyoIuQWNbw/Yj4C6FjpCT+EVpijlg9G8B//O4Z2hMR3njxwrRlj8n5L4TI3ytAMVW+5QfJpXs5POShUWy640e89YwhwIbT3gNWYz2pIv5C6AhKO7QGOf+8iP9M+dFDB+hqji1Y8R8/5z+LE77ag2LaCLubgvww5Psg1w/5fnAz1Y9XNgdHm9D+etjv7uziLRe2oDL76zz5zBHxF0JHIP4t8XLOX2uN71gTjoF03mWwokfCQiMQ/4ZH/m62HLkXRo2o5/rMz8IwaA0K/H/AjoEVg0grxDrHPN2+HnN8r1xfZNMOhyd6FOc1z+wQx0PEXwgdIzU5/6KnyRU94mFd6tkA0gv86ilfM+FrBeI/mc/fK1anZgqDZXHPDYAuUBJ2NFgRsOJG5BPLQU3PU7N/2DzX+y8qcttum+8+5XDexdN6iikj4i+EjqCLV1siQkvcfAVGc0UR/xmQzrul6Hohki/6Of/a2j6u5wt7GoqjUBgpi3u+39yHwpQvw4i5FQM7DvEuUI39zOwfUUQtzakdmpetc9m0w+ZjFyqaGroXg4i/EDoqI//mmC/+2SJdzbH5PKwJ2d+fxrYUK9oT830oY3D9K6cc1WWTFwxujnx2CIBodj/0bMfu7wMsinu+D+kg9+6LfCDudhIi7TCH6cD9wxarWjWWgjed6fLDbQ4/293CG85u/L5E/IXQEYh/c6xC/BfwQq+nDg7zxhvu5YLVHXzzT2YpBzADKlM+Q5kCi1vmeBD1XDOhWnLNDFVH716OfHczsILowB/A8XCKLUAXbqQTknN7uBOxf1ixqsVcZWxY5rG+w+M7T7fxhlnYl4i/EDqGswWaojaObZXEf6HaPXf1jPK2m+5nOFtkf//CXI+Qybul34cy+caLv9bliVU3mFjtg1wv5AdMqqb0WEykbsfNLdoBVoS8fzUSbVoKCY0dNf8vYgELJ121f1hx3hJzPErBRy8rECv2GkNCg/cl4i+EjpFsgZZ4BIBmP+e/EEs8dA9meMuN96M1vPjMpfx+Z8+CdCWlKsR/MH2Mjp8qz/so5AfL4p7rp0qgNWBHfedMHBKtk6Zm8q65P+qn6IPMlLtwdJ/hHAzmFKtby5PQV6zxTFpqFt5zEX8hdAR1fYAFnfb5i+89ykiuyM3XX8p9u/q5fesRhjNF2pKR+T60KirTPuOK/7ied1/g3Qwl14zWYDkmcrdikFg644nVgi/ywZy+7bcXWUhz1PtHzPmf1Do3ByXiL4SO4WyB1oQf+QdpnwUo/lsPDfOaC1dy9oo29vSalM+h4cwCFP+KyH9kAFK5as97YcDk4UuedwBV4Xlvg1jXrB5jzj/EqG1EvxT5L6DqDoHNszLyn01E/IXQMZItsqjJtHAM0j4LrbhbKldkOFtkue/uWdYWB+DQUJYzlrXO+v4PD2VZ0hIr+eGBsZ73/CDk+0jvHwLMSqTBA/dB26B5fJXnfcWcumZqCcanmC/6vuNzYUX+w+agTmoR8ReEWWEkW2RNp3FOJyI2llp4Of9DQ1kAlvuiH/w87G+fFbQGN03/0CDP+5fH+MerOnjN+myN5730YFAOWDHSuZbS1iGvHZKztCR1BozJ+fvj0ELK+e8fVrRENf5bPeuI+AuhYzhToNWP+JVSC7KbVyDyy1qNEixuiWGp8qBwzLi5mui9ohxBfhCAA71R8u5qHtq5i9csHTa5d6fJOGfqkPbK+fjB3MKajA4Icv6B+Jcif71wjnf/sJqzlA+I+AshxEz4lvPmzQuwsuehIbPwaHmbSftEbIvFLTEOD2Um+rMpeN7z/gODBU0Rf1FTuRzBEdco49bBFohPbttM+y9dc0QzOIsXJjMh74KldEn0g5z/QmrktX9Esa5DxF8QZoVswSXveiW3D7Agu3kFkf/StrL4LmtLmMg/cM0EtshcX7laZHGUkrBrwKooRxBdZFw0k3AkZaLh7f0WngZrkuA4XTAPWN6sGVqgkX/eLUf9AH4XR4oLZMJXaxP5X7HanfzBDULEXwgVpXLOierIf0GlfbwCBweG6GpyiKV3mgJi+T6W23l29QA778NfzWR+WlF/UVPCOGdmOLF6NG3+PlVQ7B9WrGmbWCFTvrtz2YIWf1Ul/o5lzmm2cv6uV04tTYWeNORcxepJXutGIuIvhIogvdNaFflHGJrLcsS6TjGxfK8fwRvP++HDK1kWt+HgL83Eqh1jWXMH9xyMQ3Ll9HanjUA3R6f2+CDyB9jaN7n4ZwqQcDSL4pp9QwtV/CFaIcaz6fb51S6LD98R5cZr8ly6cmo7CJw+q+bI6QMi/kLIqCzqFtAcs+keaHDphFKd99Fx6rwz1vNux0ue90OZmBGC5KrSUy5rcRjJK0bzUxdygF/vtvjg7VHuekuWxVOoY3MkpVjX4fHMgGJrr8VVaycWsFRB0RSB9vj0JnynGx3PhLxX9vhDhdunwVr7hwMW7/91lLynuG2PPWXx3zccLPAS8ReEWSEo59xaM+Gbyk0z1+oVKnLvo2ZiNd9bUee9Jo1kRSpWrK6cNDVzeFTxrOXVwrG8yQjD4ZRiXXTqIrF7yCJTVDxyxOKlp0wuRkdSipPbNJ6GbX2Tq3OmCImIpi2mGc5NTdT/+0mbj98d4ZzFmktXuLxwjcfFK2Zv9nVMzn8WIv/Hjij+7JdRTm7XRG3NgwenPrIFq3sl8heECegZyR1z8bDK/r0BzbHI2Jy/9kzJgVITjxFTiiCI3ovpCgEve95NnffFMypHkCmYCHpZU7UQLGv2xX90eq6QwIHz+NGpif/RtGLDMo+oDVt6Jo/k0wVF0oG2mEajGMmbq4CJuK/bIhkxZRa++pjDfz4Cj74zS9ssFQQdm/M3P12vfH7ZIqQLsOgYqmb3puHtP4/REdd86xU5vrXF4SsPO6QK0DSFBdn7hxVLmzTxOVRkEX/huOLR/YO86sv38P13PZuLT1k07b8fKbVwdEqe92YrxWiuiHf0Hqx8YIscotTAIyAQ90hL3RZ8jeJwquyeqSSI/A+lppdXDyZhHztaHYlqDfcdtLh0hVcax/Iu9GWMEC1r0vzimcikaaZ0AZIRaPeFezCnaI9PPDjtGLC4aKnH116e54fbbP7yN1EGs4q22OxEvoUpuH0+e7/Dph0Od745y3QraDzRYzGQVXz5pXmWNsHG5R6uVjxy2OI5J00+4O4bVqyeo5o+AQus64IwV2TyLr/ZdmS+D2PaPLZ/EIBNj3VP/mDPNfn1zCEY2Ql9DzF89EkAWru/Dju+Aru/SXPabEv3PAbZoyZqTyw3+fbKW3yxEX5rGgn3Y+DQqFGmZTXiv6SpHPlPh+Gc+fnEUQtd8ZR37LV4009j3NtdloEe3+mzNKk5o9OI0fZJUj/pIiQjuiT4k3n9ix7sGlCctsg8f6sv+MFxTpe9Q2pS107ehUjFadSr6nk0rTiaVnz3qfGv2u7eb7Hxa/Exx9qXqU7bXLTMw1KaBw9NTWJ3D1pzmu8HEf/Q8ssth/iTr2/miQND830o02Ln0VEAfvXkETzXrxSZ7YHRPTD4BBy5E/b9iJ7Hb8Lb/iV45ibY+z3o/hn03sPI6CCW0jQlO0qi3tzcBsCoWgLRdmOZnGbv1UYSiP+KGvGPO9CZ0KUrg6kSTMIO5lRpYhHgt3uNyG3vL2874ldwWNKkObPL7H9r38T7SxXUmMg/oD8DB0aq/37fsCLvKdYtMs/f6o+lI/npnZen4XP3Ozz/v+N8cfPESYxczYRvKfKvEP9s0Wy84dEI4635e+iwRW9GlXL0AYH4dybNPlqicEbn1MT/0KgZeM5bXFm2WkMhB+khyAxP+hzHgqR9Qkow8Xn3zh7OXdU2z0czAW6+asXqzu5uLGXy/g///kY2LvVXvAYVI60YKS/OC36whD87v4MPXlw9kTviRmiOgrIr3D6++IzmK3q1ziOBuNfm/INt0438h3ImjXMkpXjsqMWaNhet4a79Rph2DVqAeZ2O+vte2qRZ2axpiWq29pXvr0emAElH0+ZH/kPZ8vF98vcRHjtqceeby6Hy0/5gc1rHsUf+wzn40O1R7thr0x7TfPtJh/+1oThuzjzvmro5AUqZ+YZKt0+2aCyrR1KKH26zecs5Y885GMh609Wflb60psPJkMwPQToLhSx/1FbkgX05ijtHcFyzjULO/1m+RYdz/DSa5/StaXiq4nHB83eugvM/NvUXZ4qI+IeUdMF8sO/Z2ct7XrBu/g6k5HkPJlaDOu+B5706h7Cz9xSuXFPgjn1xftm9mI2njP2CPnnQIlVQfO2JCNdf6Fblb4dz0FrjlGmOmP+PzKHVfyIOjSo64vUn/5Y16aqc/4OHLDYfsnjLOUVaxslGDeXgWctdbttt89hRi1eud9kzpEre8t2DFZF/uiz+SsGZnR7bphz5+2mfChHf1mexZ8jiSAqW+l3Id/ab/QaRf3Dcw1OM/HMuvPbHMXYPKT713Dxr2zVvuSXGLTtt/uiM+oNUwYVYZTZHe7RYWWK5URgahUKWU7IuZ7dncdwsBx7K4Xoj2DWi/cbuPK+N5Djj3jQ8WBbwjxRyfNTR8L3yLt4BvMMG7qo5GDsKkThEYhCJk84nGKSFyOJOiMb9++LgxPyfs7PqV8Q/pGT9GrcP7hkgW3CJR2bWLGNc6rbgq/G8V5YjKLXgG1vnfSgHPRmHC5Zrctrj1t0OH7vcHeOafMJ3qAzlFN/favP288pfnuG8GiOSQeSfmiTy19oIz2w7Mg6nxjp9ApY1ax45Uk4lfPL3Ebb0WHz1MYcPX1Lg9We4Y2yWQzlFVwLO7vJ4/Kh5bX63zzxow1KXXZXin1JELE2H79Y5o1Pzo+32hGUeMkUz4Rs4dQb9yN/TsNtf9PXwYYurTzWR/tMDipUtXskFU4r884yP9qCYg0KOw705YoMeXz5vhCvb0uhclg+1FBnanEOnR1GFrP/Ysjh/LpWnI5uB72bMtmKeRyMadmJuwCdr97nZ/+lEwTGC3FpIMKjijNitdCxaXBLqTbubGCzGeftFTknYB9w4b/9VG2+9wOF150bKgm9Vf9c+/JMohaji+S8c59InvX+CF+bYEfEPKUEDjnzR46G9A1y+bgbNNKpa8JXrvBuBH6zvebf8HqtT8LwH7BzwI8YOj4645jd7o2zpUZy7pFoonzhqsaxJs7LF48bHHN5yjlua4BvJl8UmoNm/EhidSHwwbpDvbXW46y3Td4NMh4OjaozTJ2B5k6Y/q8gWjT1wS4/FH59V5Ol+xd/cGeWXz7h88xXlE3E9GM4pWmOa85bA97baFD24a7/NyW0eL1jj8S8PREqOnSMpxZKkLgn9mV0eqS0OB4brlx7Q2rh9miKaiG2uooZyCrTHwf48HW6GJpXl0J5RSKagkGXlUZcN8Qw8MQqFHK2FLP/PKXDhzjQczdRPjxTLwrgG+FkM2G5uCvhAcOejGLGOlCPojIpzsNhGpH0xSxfHSts/93Azpy+L8fIzjGD/9e9b6GiJ85Hn2rztV60MuXE2vUmZ+kiY+YGr/zOOqxV/trzA/7m8/Lm+cV+Mxc2at59Zfu07gIGWGL8e1Lyuuf6Hq+DC4z0Wf3zW3NX0CRDxDynpgktzzCFbcPn9zt6JxT9oweemyytWcz2mkFil5z1Ytars8oIm3/P+30/aDGYV773o2GvoPDNgFGldh2bDMg/7Ts0vd9mcu6T6OR/vsTh3icfrTi/yrltj/PIZm1esN1+u4ZxiZc1CmuBKYKQw/iC0a1Bxw6MOBU/xq102rz69/pc178KjRyw2LvcmLYg2HodHFRcsqW9fCRxAR1KK/3naxlaaD11coCsB//fuCN/a4lBwy+0KR3zNaYtpTmnTfP0Jh6d6Ffd2W/zRGS5r283z7RnUnNWWpTA8yoXxHPSOQCHLJcUsr7Q8Rp4YgbbMmHy1l89xczTP2mcysC/LnXaWlh1ZeDrHKuD+wO+/178B/xsgBzxo/q+cGC+2E1jpWDnVkWwDZ0mViAe3naMJ/nFzCx+8zOacFeZvMirOi37YxoUro3z5qurPw39tdvjXAw73viQb9JwB4NuPxnlFq8vL15l83313x9iQ8FAdBZ57usM/3BthMJ8prVk4nFK4fgno3kzNhG/aXCXV8qzlHr/da5spqTqfh+39imxRccHSuS8vKuIfUrJ5l/ZkhOVtLdyzs9dPzfjiXkyZqL2yBV9tOiQQ90jrpJ73vgz83T3GQfGik926X5KpsKNfEbU1q1pMad5LV3rcusvmry4plr5YI3mTw37VaS4vOcVjbbvHDY86XLnW5SfbbfYOKc7qqv6iNUUmj/z//p4IMRsWJz1+uL2++OddeO+vo9y22+aTz81z3bnTj+ayRejPThD5+9sPjmhu257jFSsGWZwfgVSW5zoFDllF0tuHabN8gR7J8jGnyBX7UyyOZLkpkid5W4YfWDnWdqeJ7M/xZCxL0y0msv5isKNN5sda4ItR4OmKg3DK0bO240CCXLQdFsXYnG5CReNceZrDg71JfrK7mTOXRbjncJIvvszmaCHO225t44OX2bzyrIh5LmXxqm/GuHSFx+deNPnEy/69Frd7Md6zJAu+IykBvOIshxsfdegecUsDvNawaYfNs1Z4LKvpMWOrsW6fIKW3zp+MfmbA4iJ/pfWBCqdUb7r8u9bQl1V0JuqL/4+2OzwzWH9h3qN+Cu9CEX9h1vCKfuRuUjPp0V6SKsflXQN8YXOMwS030B7z/OjdrxRZU+e9HjkXvnS/w+WrvHHrmNz4qEPWzwt/9v4IN14zSX5lHHYOWKxtL9dkv2qty9/eFWXHgOI0f/LwyR4LjeLcxSby/tPzi3z0d1Eu/1acHt9O964LqyPDpqqcP/xku80TPRZ/vrFAexzu3Gdxx16bjz67QKYIn3/Q4cCIqlqKX3Dhfb7wr271+Kd7I7zoZK/0mO29mo//xmNJJMOGzjTndqRZHs3QamdJksXyUxvZkSx/6xS54nAKflsbaee4JJ/lqViO5K9z/BqgD/ixOYaXAi+NAveVz60deKMdJzIYIxqPs8xK0JtJkKaDM5YuBSfGd7Y1c96KKJesifDpB1o4c1m0Ikcd5zWb2ti4OspHn2+VxDrg0LDiDd+O889n5/mjM1y+tSlKtghXbshzy10RNlk2zzsnz63dMbY4WXryil06xuplWYiUX7+W2NQnfDP+21ebenvrOS7/9ajDfzzi8OnnmUFke79i54DFp88d+5lzrOraPtmK+ZxTfaF+ZlBx0XKzLXD6rOvwqiL/0YJxE3XVWRm80R84Hj5ssa5jbDDwyBGLroSe07IOASL+JwpaV5cjKI6Uo/f8gF/nvUw6vZKE7fCcFVk+T5x7+1eVJuSmylAO3n1rlHu7bX620+P2N+XGTDYOZuGbWxxets7lzE7NP98f4aFD5WhqOuwcrE6HXLnW5eN3a3620+YvLjaKsKXHHMC5vmf6Nae73PCoR3scPnNFnhes9sZcfsds4wEfKZj0zv++M0LR1fzmaZe/edYI33usyNUtWd6xeJihkRy7LY8dfxhl1fI0FHO4uSz3783zypEc/3dJmsXRDLtyeSI/yqCjWXQhy+lunpuDHQ6Mf45tKF5vx4kOxSBbjrBpWlSKtL+9tZmcipEhzgcuc4jFzWOeGk7w13e38IkrbDaujkIkxl0HE7zt5wl+cE2OZy33+LufRvlDt81lK11e+EIjiF/bE+OSuMf5Zxb46t0J/mpZAVaXB8hiS4xtGQ2RsQIalHNOOka82mLlRWi7BxVr202KDuDhIxbZQjl1V0lrVE/Z6pnxnyNZo16rWjRvOdvl20/avOHMIucs1tyyw6TGrl47VnhtS9dE/hDzz2NVi6nPs3Og/GE5MKJQaM5b4nH3/vKkbeDxX1RnVfNK/0qtNk0U8OgRiwuWjv1MzgUi/scBQ+kC8ahFTFV0aSqMGlHP9Zqf+UHA/yQHufcgcq9T5z3jRUlE4PzlMZoimt8fsLlyrUmRfPUxh5uuyY2ZSK3k0Ci8/Wcxdg0q/uiMIj/Y5nDrLouXrasW9Zsed0gVFO+7qMjqVs3XHnf4zP0ON1+bn9YHPls0l92vO718TEuScNlKj58+bfOhi3KoYpZ9h1wua8rTNToMg1nihSx3XuJHzkNZeLhmMtF3hdwSydG1I4u9LctjTpZEIHQPwpXBDn8NXcDno0C3uWkUBRVjnZcg3hSjPRIDJ05Tewf39q/ktI4ojw8kGdBxXntuhCVtMbIqzu5Ugt5igr5inFv3NbFzNMFP3qi440CCD94R547XZEvRZyUR4N+2xxnJK159WpHY2eU0SXRAsUXH6bbybEwasRvKG5EKyiacv8TjD902z1tdfp/Wtmt2DVpVHv9KVrRUi2Al6UD8/Si8LaYZygXrB0xxuiVNsKrF45HDFo4FK5u9MeUiWmPQXbNwytNwb7fFZSurxTHoHJaIjH19PnxJgV88Y/Oxu6L86NU5btlpc9kqj6461UwdVV7h62kTvccr2jyubdM8M1COZg6MKJY2mcV3/RlKDqi+dPUCr0riDkQtTb3um0NZeGbQ4jWnz4/HWMR/IeG5VakZCsMc6O3lld/J8fJTRvnUJYco1wHWoCJG3K04JJZOq5hYpmhcHREbLl3hcedeizdvMlE8wLefdPinJWM/lOkCfOdJc2mdc+HrL89z6QqPhw5b/PvDEa45NVf6og7n4GuPO1x5SjnP//6NBT5+d5S79mqevzxdd9FLWZjLv2dGc/xbpMDFB9JwpJwOuTFr/l59wxzrp4ID/dl4Z66qPNbBrUd1saOQYNCNs/GkCGcsi+E5ce7rSdBfTPCyMyIo/7G/2p/kY/e28u8vt/j5ngRf3xLhry4pVE1mr/TgQ/8T5aFum6Sj+c61OZYsNa9BHDiz4ohOPVXxih/GufHJAjH/G1lb2qGS5c2akX4zr1FJl59z7qvo9Bh47gMP/mWrPP7rMc2L15T/9pR2zf88bZUWl40R/2bNXfusupOWGX9VbCD+7THNYM4s/OoetXh9u9nPhmUeDx60aI+X/f2VtEY122rSPvccsHjrLTFueV22KhApRf51HFdtMfjoZQX+4o4oH7srwr5hi/ddVD/NaFvl2j5BXb9KG++pHV7pShJM8LGqxaMzoXG1YiALnQno9V/vrjo5f6XMwFavyU1Qa2k+JntBxH9uqet599vv5ftM5cgKz3vWtXjPrSfRn43yu+4mSJ7UsENJ+w04AC5b5XLH3igDWc1nrsjzwCGLn++0+cRzCgQNr7Tr8s1Hinzv8QKFXI6rF6f4s7NHWV3MwjNZPrs8x21PF9n321HWJI04Hz6S50s6x8ZsGn5ixPpthSxvjOeI/WaK0Y4yYh0nzhkqQYv2r2biXRCJY1lxvrGtifVLYly6JsLH72/lilMjXHV6ZKzIO3FjA6xzyfEP34/xVK/Fs1e6vPmqPChT++SyOof03DZIb47z57/VHEpZ/On5Bd6zoXoewbbgn64o8Ne/gb+8pMgFS8cX83OXmLTEjY85XLHGpTWmJ6wEuapFM5DVXL6qWjRaY6ZDVV9FimHYF51W34P/3JM8HnpHdfXMte0eI3njAoJ6kb9HuugwlBtbrTNI+wST5u1xKHqKJ3st/7nN9g1LPTbtcDiS1jxn1dgUTGudnH9Q5mIgW73+ohT5j6Nerz7N5eatLt99yiFiaa6ssxAQ/Jy/X9Uz6z8k7pT3c2qH5pe7jK027pjIf+Myr5Tb702bSd5SaYc64m9ek/odzh45YpXSSPOBiP8E3LH1COm8yyvOXzGlxxddj4HRNIsTxXE87wPsHbL4uwe7eNc5g2xckjGr/UrRe2uVMH3yzgiP9zo8f7XL7/bZHBqF5c3j7z/ghkccfrlD84lLhzl/UX3f9DU5lzNTadic5s3ZLBuX5FjXkqFpX46rUlneqXJ438+CMn+v3CLXAdcBxIBh4N7yPjcAGyLg7VHoSJxhnUDlE6yJx2iKRyFixFpF4jx6NMm9R5O8/9kOdqy8kvF3h5L8+xPNvGR9hNef69DaEjevj1J85QGHLz3ksPXaLFRc4ESBzaNRbjhi8fmL8nzPjXH1aTmYQiXFSlqimpit+YfnFyZNRyUj8LJTXb6/zeE1pxf56GXFun+zrkPz49dObXL7Ly4u8KvdMX620+H0RRMf+8cuL5AulIuTBVgKOuJUif9QDmJ29Wrh2rLJgUAHV31La9IXQd66e3Rstc4gCg+EOLjCCBaindJuziXI+3tasb5O5N8S1YzkqFpM1u8vFkvXxAmZghHp8ay0SsGnn1vgZT+weP5qj7Zxykvbqhz5Z+sMKKd2aDyt2DNknDqHU8Ym3JUs5/FPpyz+45WCNqmwsdsfPWKxfpEed2X2bCPiPw5buod497cfwvU0qxclOf+kdnNH4HkPovf8MDrXy23bh/nMfRY7ByO8+KRRPnRBH2d35qs874eLS3nzbQkOjFjcfaiJr1yZ54o1/hfdLTIwlGJoNEchl2XroQIHthX4/LpRNnaluaHbpf/+FMvb0qVcdXd/DruYZVmswnddyPL2bJbrVRFuH//8/hJgEBiyiEfinB+JQcpExy2JGD12G33Eec7JETwnzje3NzPqxXn3JQ5OtDZtYn7/+rYmPnFvM6fGNc8MWrz93CJ/c1mhSqwBurfbfP5glJevylZN/P306QgPFWzu36L4wtOad19Y5H9tKGJhnD5rWnX1En2fa08r8qvdMW541Hycz108/Ujqg88qki2aFMhU+PAlBc7o9HjrOe4x+/krWb9I8+rTXH603Zkw5QNlsa5HZ0JXTS4OZtWktfWD57v/oEXM1qWrhIAVJXup4uyu6n3XRv5BfZ+HD1dH/md2auKOJltUrK8zuLXGzPxJqlBed9GfCfZR/QKni2rcqD/g9E7Nt1+RZ9UElTJN5G9+D9JXlYNk2e5pVoUXPcWqVs1iP8IP7J59GeUHD/X30xajbjG+J3osrlgzfTtwo5g38VdKXQV8ASMNN2qt/3E+jkNrzd07elncEuPM5a0AjKZHef93NtOZtFFo/vLme7nldXni2njeRwuKzYccBlNFRlJFHjlss6fPY0NzinevzfBod57v/yLHhYvSbOhMsyqeoZDLsmNvgX9yM5yzJMPAaA7n9iz5SAbHzWHpIh2YVYEA64FXRoED5vapCGaRjLKM6yMSJ5tKMOLFaVsWI9HSApE4W4eS3J1Kcs0ZDo/1J7n7SJLl7TE+dJlTleO+5LttvO4ch7+6TI1Jgyhg8wMO/7bZ4Z5zc9zTbfGJkShffmkOZ934wvr6c+CLj0BPRvEfV+XGbf+33v9S7eyv9j4/3W/x7JUef/PsAv/yQIR/vj9CX0bxt5cX2NGv6k6AArxwjUdzRPO7fTYrmz06j6EZx7On2G4vYGkT/Mn5jf3ifvBZRTbtsGdk++tK6JJogsk1T1Yjf0WzcbYM5kxN+dqrmBX+8RysU1CuPPlqfgaVPR85YrG8ySvl5SM2nLfE44GDdl2/e1BvaThXLsAWRNTp6oyaWY3sTP4aTdZC0bbKPv8g8o9XVP5c26ZRmECmK2keuKoi8u/xX+e+TP18f0BbTFdVTgWz396MmheLZ8C8iL9Syga+DLwEI28PKqU2aa2fmpUdam2qQ+ZGIT8C2SHI9NJz9DDff2APu48M0UKWF6wY4fLOfh7ZD+8ddnnh0gEo5tjTrxn5bpp4JEsxnyNWyPKCiiqHbwWTCikAB+F1CmPLGIHCsM2IipMlzlKdYHl7lJamOMm2Nv5wJMFv0klSmHTIaUujdLTEsaIx7Gic05ZFifqFnj7wuxa2DiX49Zs9UIo791q84+fmm3a14/KVF+cpuPCn345x8lLNu5+XZzXw+D0OX3gswvtXZqpqmB9xE0SjBVD1V9y+9nSXL26OcPNWmx9uszl/icc1k1hBkxH4n9fmiEc0SyboFRuI+I4Bi6t8h5LrmUVcbz3H5awuzX9dnedT90S46XGHzoRmz5DixSfXF9u4Y2yfP9rucM7i+fsyzZSTWjU3vyrPyuZjzwEvSuhSwTYwaZ/JxN+24OQ2zdP9aky+H8ykZtTWHBypI/41tssgLXQ4pbisRnyvXe/SGqVumqOlorLnyhazLRD/2sg/U1SlwWYmOBVVPQPxj1UoYiICK32n08pmcwyrWjRtMYhYuiryXzSh+JfnXgKCngcT/d1sM1+R/8XATq31LgCl1M3AtUBDxT8/eBjr3y/GKoxi6bHCsRh4LxihBuiBQo/NuTrO+fE4rUWTzuhLJng41c6SaIwnskmsSJxnnxyhoyVGUzJOLF7tHAluWRXnNwdi/Hi7wwOHLD73wjyn+W30IsDFBXj8MYdzFns87yRvwr6n55zk8NP9EY6kMyxtglt22rRENdedW+RLD0V44KBF94jiUMriH15QTjAu9S/ZU3lKuc/JJswA1rRpLl7u8qWHHDyt+OyLclOyZtar/1JLMmKsfzsq7IP7hhU5t7xYSyn428sL9Kbhn+83b9BErQuvXW/E/1hSPguJi5bN7Pg7E7U5/7HlLOqxtl3zdD8sqWNXtJS5OuiuF/kXjJUxKCfRXjHQ1Kan3ny2y5vPrj+AV9f0N3/XH0T+NTl/E/lPekqTUhX5u2PTPmAClWcGLE71z2VFi6l2alw+ZfFfM8Hnvi2mGckril55nmbAn89oj437Z7POfIn/SqCyVN0B4JJG72RYx/h5+mJSJBjVcdLESRFnVCfIqRhnLdFcd16WxS02RJu4p6+Tj/xhMWva4Bsvz5da3SzNwdtujnGo1+La9UU+/bzCmLxoPeLANad6XHNq/Um/RATev3FqtW4uWeECEe47aHPlKS6/3mVz9aku791Q5IfbHP7unggFz6RUXlDh4W7xB7bRgirlYzMT+KQree3pLg8csnn+andMFDdT1ndodvSXR7vgsrhystNS8NkXFRjIKn5/wC51fqrH5as8PvLsAq857dhrB50IdCU0o4WyQ2UoB2d1TUX8PcAulV2uZUXz+JF/peWycn4hmOydCvUqe/b50fHYyH/yz+5UcKyyxbNs9ax+3nUdmvsPWuzze+wGef2uZHXkv2GCQbut4qommBQOJrPDGPlPilLqeuB6gNWrVx/Tc7S3trH6Lf9Oc9yhJe7QHDO3ZMQmqvIVlssU5Pu5vKOPu1b3o/PDWNnym9mq4ZsvbWJ/KsEVJ4Oyo8x1E7Szukxjjfu7zaTcaEHxyvUuiQj81SUFPvwbEzp95orqxVNB6YLKujXjrZCs5eXrXe4/WJxRMbbxWNfhcU+3g+uZCOxpfyCodYLEbPjPq/Pctc+aMKVjW/DuC8Mt/FC2G/ZnFCtatJ/zn/zvgonuemkfMBHv7/ePndEMWjgGxB3jLsq5asKJ6VpKNf2ruoD5aZ+a2ClTqF9HZ7oYt49v9Szl/Ksfc2q7R7ZortxXtZQ1ocufWPc09GcnzvkHA+JQrpweGvDnCzom6XU8m8yX+HcDlab1Vf62ElrrG4AbADZu3HhMr5BjW1xxxpJx7k2Ak8Cs2SyjADWmVEKK9Z39rM/1QWEQMn3BUVaspg0sm/6tweu1HQs2LvO476DFUM58yYNJylef7vL1JzyOphXX1iz8KZUrLpQvp4PL6Mmip6YI/MuLZ2f14fpFmryrODBiLpm395nJxnoLd5oiTLv0RFhZVLHQa3HSRM2TNVMHI3JQv3sYmMj/SIqqiqFgIv/a9GFbDI6mJ3Yl1RJE/kEV0kyh7MBJ1aZ9inBSgyL/wO2TreP2gfL81P5hiw1Ly8FFV9J8Zgezxr460WAURP6Vds9S5B9C8X8QWK+UOgUj+m8E/niejmUsSoGTNDcWj71fe/7gEKzGDRqU9JcblQSLUkqDQ8x41q24GSiOYXC4ZKXHb/dFODCieMOZ5Rr1loJvvDxHpqjG2M2ag7RPZeQ/hZz/bBPY6HYMGPF/ut8q5fuFYydwOvVmFEM53345yYQvwIVLzQK/K+vUwAHj9dcoDqdUVaNxU8u/+rHtcc1glmk5WWoj/76KVpDp2rRPgYZM+Fbn/M3PxJi0TznoqDyfLn9xV29pgdf4+yk1ucmVA7BSzn8SG+5sMi9ff611USn1PuBXGKvnTVrrJ+fjWI4JZYHTZG4TDg6p8uAQDAz5AVMLPxgUzBPWXDnUHxwu9fP+edekfCoxucSxX7bAf10ZPdUuyZ8PgsnbHf0Wzz3JY/eQ4iXjrMQUpk65xIMqrSqdStpHKXj9meO//pV2z1rxr72CbI9rVrcxoYGhlqhthDco7lZpVx0b+U/u858KjmKM26c28l8UD0pWqKo1A11JTd5TpS5oU4v8y9/pgaz5bs52V7iJmLdda61/AfxivvY/q1QNDnWoWiiWNmUdgraGweBAOUoAwIpxTnucpkiM1qiesiukujm5oVSMawpe6dmiNWZSDDsGFLsHFUVPcbpE/jOmMucfpBmmEvlPRuVCr0rSfo2oSj74rCKFYxjHW6LlEg+BY6ktpsdE/lP1+U+GqeppnrveIi8wg+K6Do/Nh6vXXyz2I/3tfWaEm1D84+UJ34D+rJrXfD8s4AnfExplQaTZ3OpRKvCWLtcAyvXi5Af48EXDdEbTWJnKEs2qfNVgx0zBN//KoV6LwnSxekn+fLG+w2Nnv8V2f7J3IjePMDWaIsaT35ehIvJvoPiPjhXi2lTJsTrDWn1LJJQne09q1VWRv+tBzm2Uz78y529+1lule2qHZvPhmrSPP+AFn92JI3/zc6gilTWQUSyax5QPiPgvTCwbrBaItIy5609WU6f650i5QFxh0GzzrxyaigDrSWUypr6PHSdbMJ/wRnyBZsK6RZrvPWWxvU9hK83aCXz8wtRQquxEmU7aZzISETM5OVb81YRF6KZDZU3/IPJf3ap55HB5n6VGLo3y+fsfuVzRuJTqTcVdvNzjjj26lPqCcnptW5+p8d8xgZDHbGMhHayN/OfR5gki/scnEwwOQEXXrjTRYoqovY1RtQicEcj3k05ngSUkC92Q9jBXDvGKOYe5GRXWd5hqkXfuszmlffzaKML0CBZ6BWmfqbh9psKKlrFe/3SxcXNHLTETEYPx+EctzZKkJlUsTx5MVMt/ulRX9VTj5t9fc7rLq0+vruEURP57hhSLEpPPb7TFxub8p+OGmg1E/E9ELAesVtNfF2iO72LUWQZrzgUgvf9pYAeJta/1q3YOl11K+UEzWV2JssoppQYODsGk71O9Fi87VTz6jaIzoenPKAb9NENrnXIKx8KKZs3eoTr59wYIMZjIf99QOe2zKGEmRdMFSr0EslNcozIVKnv4Bovi6qFUhTfDpyMOtjJ1/TsTk6e52mO6RvwnL7g324j4h4CmmE0qVxbXbNFc3sZaVtS3nHqF6gnp/FB5Qrow6A8Owd9pU7nUivnVS6NTHhwqqzuKzbNxLEpodvSbtE9TRFf58mfCimbNvd3lEDfvmkqXjRBiqK7p3+/Xy0lGzD5yfn/dIPJvxIBT2cM3W6wu6jYZljIOu570xPn+gMqyznnXlLGYT48/iPiHguZYhNFc2X6RzrskIzZqvLUGVgSibeZWj8rBoZiCwlCFlbUfvBxVsVJpcAiuHMzHriNezk9PVsNemDrlnH91rZ2ZsrLFTMgO54xQl1s4NmYfLX7OX2vj8++M65JVOV3wxb+0QHHm+6uu6jl+2mc8uhKanvTURLw1ZjqBQbmom7h9hFmnOWYzmitbJtJ5l0R0Bm/9ZIODm6+ekM4PmQnpOoPDuraV9GaSnNYyCEW/98E02lEKY+lMGEfMwVE1pRpUU2WFX23UPG/ZgtmonH9rDAp+lN+fgTWtuvTcqYK5Esg00KnmKF12+7jjp33GoyupoY+6/YFraY9rnuzxr2oWQF0fEPEPBc0xh94Kr2e24JKIzmJtIjtqbtH2+ve7udLgcOZJe3mib4g1y5ZBcRCyvaBrVvVUNMQxVlYZHCYiEJVdgxantDXuiqq00GtEcUanHtO8faaUa/pX5/yhHPE3co1KpdvH5Pyn95xBU5eppX3KE77B6t6JHEJzgYh/CGiKOezpS5f+n84XSUbm8a23fRGPdvCBq5fwhsuz2Mtay/dXDA7lVpj+VUO2B3Qwf+EvhFNOzTqHcA8OgQ3xSEpx4dLGPW9lO0dobP4djNsHTDOg0YKiK1EZ+ePv04/8G+Hzr6ntM11XVOD4mWrOP11U5F3j9IH5resDIv6hoDnmMFox4WvSPgtDINuTUdqTNXaUisGhLm62es4hP1BeHZ05An6DmHJdpYg/ER30BF4Y5z5bVIpRIxZ4BXQljf0yyF2nGui8gXLkHziKFiXKjezTfmHCTCnyn/n+bAWuVmgdTPhO7++DRu5TFX8wxd2CBWzi8xdmneaYU+X2yeRdEo2ygMwHdtzcYovG3qe1mVMo1VVKV9tYKweHACtSrsZ6AgwOlUXGGin+ljK9cR/vMSnDTIMnfIP5id2D5vkXxTVN0eraVI2cZC51ttPHlvYJIv+uCYq6BQRNW4Zyqpz2mcdGLiDiHwqaYg7pvIvraWxLkSm4tM338t7ZQqmKwaFz7P1amyuHyrRSpVOpanDw00rB4BCkltTc9nKYLtWRf2Ofe8Myjx9stSl6FZF/gz5KbdHywikwg1h15F+uwdMotw8Yx89Ei7zG4/JVLteuL3JW1+TzKq0Vxd36s6ZP8XzHXyL+IaDZb0yayhdpjUdM5L9A0j5zjlKmj4OTmGBwyFTUVUqVrxwKA5A5TNWVg8Yv1R0MDtF5HxzijpkoTRUmb94+XS5c6vGNJxy296tyqYWGrfCtFn/j8/drU/kRf6ZoFlc1wq/g+G7kcuQ/vb9f2gRfeMnU+l2U0j5Z08hlvvP9IOIfCpr9T3UqZ8Q/nXdJhlX8J6Oql0PX2PvrNPopDw6DfqOf2l4Os9vopx6diUD8G/u8G5aage+RwxYFfwxsVHXYYCVykPbpjGsc/2MaRP5B28hGvIS2ZY676E1/kdd0aatI+5i6PrO2qykj4h8CmvzIfzRbhDbIFI7znP98MqNGPwPVjX5KaaWZN/qppTMB+4YbV9cn4KRWTVdC8/CRclPzRkX+cQcillmg5lia1lj5lUpVRP6NqoHvVKR9cseQ9pkOwfswlDNWz4naPs4VIv4hoNmvmBY4fjIzXeQljM8xN/rxJ6SPsdFPLUHev9FpH6XggqUejxy2WNHsYqnGFeRTytT0788aD3xQUyfp6Cqff6OuNGwVPGf9Wv6NJLiqMZH/2F7V84EoQAhojpnQLJVzKboeedeTtM98MUuNfsp1lczgUBb/xp/ChmUet++J0D3i0dSgFExAa0zT75d2CEhGypPL6WJjavlDOfIP1j9O1+0zHSK2mYcZzCm/lr+IvzAHNJUi/wJpv8WSpH0WKFNu9BMMDqbRj1nrEAwO0Ol0Aoto9w5CrsKtpGau1kHe/54DdsMbAgURcmXpg+ZoOfLPNDLy98U/NQeRP5g6S0dTZgCb79W9IOIfClr8yH8055LN++Ivkf/xyaS9HMzg8CLrCD3qMK2rngOF/ppGPz4a3xobq3YrTcJ5SzxspTmaVg0tHwFlx0+l+FdF/oXGzTEEaZ/ASTTdRV7TpTUGe4f9NQyS8xfmgiDyT+WKpH3xl7TPCYo/OFx0WgsXnbZu7P0VjX5MRdbhignpfrNdWVRNStc0+klG4IxOzZO9jUvBBASRf+XCqaCmP5iouTPZ4Mh/DtI+YOZfnvAXyM13RU8Q8Q8FJbePiL9Q0+hnDKXBIVUzOFQ3+tnQtZgne9tJWjlzX4Ma/bSOE/kHjWmyxcalfRwVrB6em7RPW6y8r/nu3wsi/qEg5lhEbMVorkjGz/nHJecv1GPSwcH0ctgwuJ9vbdtLMtkGTa0Na/QT1PepnBBtikD3iPk98Pk3Artmwjc22zn/inOSCV9hTlBK0RRzGM0WyZQif3nrhWPA7+Vw4ToH2EuyaRGsuKh8v1eorqs0zUY/Lb4Cd1ZF/rqivEPjff6lyH8WF3lBte12vou6gYh/aGiKOn7O33j9Je0jzIQ1nUmWtMRY1FwzQWxFTB+HcXs5TNzop1UlgCUsUkchlQHLoclaQqqQRHvurPj8gwnfRjuXaglstwo9Kxbc6SLiHxJa4o6kfYSGoZTie+96Nq3TDcMnafTTld4H9z/BsnUvgVYX8oMkEz2kC5p8uhdXN5NkCNID/oEce6OfUuSfn5ucfzCf0Ror73s+EfEPCU1+Tf+MTPgKDeKUrnEWqs2AK89dxc2tTaxZVS6619y1k4K3naGT3gncQWLJRljdVt3oJ9c37UY/ZZ+/+Tnbbp+grPNCyPeDiH9oaIo5DKXz4vYRFjQR2+LStdXVVoPPau+ImZlNJlsguar+E0y10Q8KJx8DTmIkWwAc4rYLzN73Isj5LwSbJ4j4h4aWmEP3QLqU9pFFXsLxQpNvTugdzQGTfHan0ejHdvuBblKu8V3GC0dBu9V/08BGP20lG+sxP0VDEfEPCU0xu5T2sRRE7QWQdBSEKZD0FykG4n/MTrWaRj9OayvQTUq1A0PEznwv6NwkjX405cp7QaOfeLmnwwS9HNp9b79E/sKc0hyLkMq5fi1/BzUHNeUFoRHURv6NSlnalvkOjOaKRB0Ly7aAmTT6OUTVymitqxr9tEWNM0py/sKc0hyzSeWN1VNSPsLxRCnn76/GapRTzfG9nqlckfhU7DfH3OinHwoDtHp9XLRkJRctGoBUat4a/QSI+IeEppiD1uYLJBU9heOJoDxJ70hjI3/HCsTfLdW/mhGTNPqxtMeP1mfKC+CKo37ZjIGJG/3o2blSEPEPCUErx57RnDh9hOOKQPx7Gp72MdF+Kl+ks3ax2mww1V4OpTkHv9GPdmelL7SIf0horoielrQugOWFgjBFmmrSPo1KWwaRv9YQdxZAQDRZL4cGI5aPkBBMmvWMSOQvHF8kY7UTvo2JWYMJX4B4JHxSGL4zDilB2ifvepLzF44rgs9rv194v1GfX6dC/GMh/E6I+IeE5op6tdK8XTiesC1FImLjepqYY1VF7DN93oAw1roS8Q8JTRXinwzhB104vgncOI1MWTpWWf6mZPU8wQjfGYeU6shfxF84vgjy/I1MWdq2RP7HjFLqE0qpbqXUo/7tmor7/kYptVMptV0pdWXF9qv8bTuVUh+Zyf6FqSPiLxzPBBF/Iz+7TsgnfBuR/P1XrfVnKzcopc4C3gicDawAbldKnebf/WXgJcAB4EGl1Cat9VMNOA5hAuIRkyt1PS1pH+G4IwheGtmBrjLnH0YTxGwNd9cCN2utc1rr3cBO4GL/tlNrvUtrnQdu9h8rzDJKqZJfWiJ/4XgjsHs28rNrK0n7zJT3KaUeV0rdpJTq8LetBPZXPOaAv2287cIc0DwLXyBBmAuCwKWRE76WpQiCf7F61kEpdbtSakud27XAV4BTgQuAQ8DnGnVgSqnrlVKblVKbe3p6GvW0oSbw+ssiL+F4I0j3NPqzGzh+JOdfB631i6fyREqp/wJ+5v+3Gzip4u5V/jYm2F673xuAGwA2bty4MGqgHucEds9ERHz+wvFFYPVs9GfXthS4C6S8wxwzU7fP8or/vhrY4v++CXijUiqmlDoFWA88ADwIrFdKnaKUimImhTfN5BiEqSNpH+F4ZfYif795ewjTPjMdRj+jlLoAU4d0D/AuAK31k0qp7wNPAUXgvVqb/mhKqfcBv8I0y7xJa/3kDI9BmCJlx0T4PujC8c1smRUCr7+kfaaJ1vqtE9z398Df19n+C+AXM9mvcGyU0z4i/sLxRXKWPrthjvzDN9yFGIn8heOV5lko7wBlr38YAyIR/xAhOX/heGW23T6xEKZ9wnfGIaY1Yb5AlUXeBOF4oOT2aXBFWjvEaR9RgRDx2g2rWNWRpDUeme9DEYRpMetuH7F6Cicync0xrjl3+eQPFIQFxupFSZpjDmsXj9P/9hgpR/7hk0KJ/AVBWPCsaE+w5ZNXTv7AaRLmtE/4hjtBEAQfxxbxFwRBCB12iGv7hO+MBUEQfGTCVxAEIYTYliLqWFgNagp/PCHiLwhCaHEsFcrm7SDiLwhCiLEtFcrJXhDxFwQhxDgi/oIgCOHDtqxQOn1AFnkJghBi3vbsNQxmCvN9GPOCiL8gCKHleactnu9DmDfCeb0jCIIQckT8BUEQQoiIvyAIQggR8RcEQQghIv6CIAghRMRfEAQhhIj4C4IghBARf0EQhBCitNbzfQyTopTqAfbO4Cm6gN4GHc7xQhjPGcJ53mE8ZwjneU/3nNdoreuuZDsuxH+mKKU2a603zvdxzCVhPGcI53mH8ZwhnOfdyHOWtI8gCEIIEfEXBEEIIWER/xvm+wDmgTCeM4TzvMN4zhDO827YOYci5y8IgiBUE5bIXxAEQahAxF8QBCGEnNDir5S6Sim1XSm1Uyn1kfk+ntlCKXWSUuq3SqmnlFJPKqU+4G9fpJS6TSm1w//ZMd/H2miUUrZS6hGl1M/8/5+ilLrff8+/p5SKzvcxNhqlVLtS6odKqW1Kqa1KqWef6O+1UupD/md7i1Lqu0qp+In4XiulblJKHVVKbanYVve9VYYv+uf/uFJqw3T2dcKKv1LKBr4MXA2cBbxJKXXW/B7VrFEEPqy1Pgu4FHivf64fAe7QWq8H7vD/f6LxAWBrxf//CfhXrfU6YAB457wc1ezyBeBWrfUZwPmY8z9h32ul1Ergz4GNWutzABt4Iyfme/114KqabeO9t1cD6/3b9cBXprOjE1b8gYuBnVrrXVrrPHAzcO08H9OsoLU+pLV+2P99BCMGKzHn+w3/Yd8AXjUvBzhLKKVWAS8DbvT/r4AXAj/0H3IinnMb8DzgqwBa67zWepAT/L3GtJxNKKUcIAkc4gR8r7XWdwH9NZvHe2+vBb6pDfcB7Uqp5VPd14ks/iuB/RX/P+BvO6FRSp0MXAjcDyzVWh/y7zoMLJ2v45olPg/8NeD5/+8EBrXWRf//J+J7fgrQA3zNT3fdqJRq4gR+r7XW3cBngX0Y0R8CHuLEf68DxntvZ6RxJ7L4hw6lVDPwI+CDWuvhyvu08fSeML5epdTLgaNa64fm+1jmGAfYAHxFa30hkKImxXMCvtcdmCj3FGAF0MTY1EgoaOR7eyKLfzdwUsX/V/nbTkiUUhGM8P+31vrH/uYjwWWg//PofB3fLHA58Eql1B5MSu+FmFx4u58agBPzPT8AHNBa3+///4eYweBEfq9fDOzWWvdorQvAjzHv/4n+XgeM997OSONOZPF/EFjvOwKimAmiTfN8TLOCn+v+KrBVa/0vFXdtAq7zf78O+OlcH9tsobX+G631Kq31yZj39jda6zcDvwVe5z/shDpnAK31YWC/Uup0f9OLgKc4gd9rTLrnUqVU0v+sB+d8Qr/XFYz33m4C3ua7fi4FhirSQ5OjtT5hb8A1wNPAM8D/me/jmcXzfA7mUvBx4FH/dg0mB34HsAO4HVg038c6S+f/AuBn/u9rgQeAncAPgNh8H98snO8FwGb//f4foONEf6+BTwLbgC3At4DYifheA9/FzGsUMFd57xzvvQUUxtH4DPAExg015X1JeQdBEIQQciKnfQRBEIRxEPEXBEEIISL+giAIIUTEXxAEIYSI+AuCIIQQEX9BEIQQIuIvCIIQQv4/ohQrH5Uc7awAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from skbonus.linear_model import QuantileRegression\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "np.random.seed(123)\n",
    "X = np.arange(100).reshape(-1, 1)\n",
    "y = 2*X.ravel() + X.ravel()*np.random.standard_cauchy(100)\n",
    "\n",
    "q_10 = QuantileRegression(quantile=0.1).fit(X, y)\n",
    "q_90 = QuantileRegression(quantile=0.9).fit(X, y)\n",
    "lad = QuantileRegression().fit(X, y)\n",
    "\n",
    "plt.plot(X, y)\n",
    "plt.plot(X, lad.predict(X))\n",
    "plt.fill_between(X.ravel(), q_10.predict(X), q_90.predict(X), alpha=0.33, color='orange')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}