robust_regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# math \n",
    "import numpy as np\n",
    "from scipy.linalg import norm\n",
    "\n",
    "# random\n",
    "from numpy.random import rand, randn\n",
    "np.random.seed(6)\n",
    "\n",
    "# plotting\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# robust regression \n",
    "from sklearn.linear_model import LinearRegression, QuantileRegressor, Lasso, Ridge, HuberRegressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot line y = w*x+b\n",
    "def plotline(w, xmin=-100,xmax=100,nsamples=100,label=''):\n",
    "    x = np.linspace(xmin, xmax, nsamples)\n",
    "    y = w[0] * x + w[1]\n",
    "    plt.plot(x, y, color='k', label=label)\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot function y = f(x)\n",
    "def plotfunc(f, xmin=-100,xmax=100,nsamples=100,label=\"\"):\n",
    "    x = np.linspace(xmin, xmax, nsamples)\n",
    "    y = f(x)\n",
    "    plt.plot(x, y, color='k', label=label)\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVUElEQVR4nO3df4xddZ3G8eexzKaDdjNqu2s7BdslzQDaQNkBlImEKKTAkrYWNxazKmSTRgOrbsyYYoiQuG6bLTG7Li4EBZUsS2Owztalbtc4JrJLNJ0yQIHaWNhV5hZlqFso2yG0+Nk/7h2Y3s5M59e558f3/Uom3HvOmXs/147nuef7/ZxzHBECAKTrLXkXAADIF0EAAIkjCAAgcQQBACSOIACAxJ2WdwHTtXDhwli2bFneZQBAqezZs+fFiFg03rrSBcGyZcs0MDCQdxkAUCq2fzXROoaGACBxBAEAJI4gAIDElW6OYDzHjh3T0NCQXn311bxLydT8+fO1dOlStbW15V0KgAqpRBAMDQ1pwYIFWrZsmWznXU4mIkKHDh3S0NCQli9fnnc5ACqkEkHw6quvVjoEJMm23vnOd2p4eDjvUoDk9Q3WtHXXfh08PKIlHe3qXd2ldas68y5rxioRBJIqHQKjUviMQNH1DdZ08/a9Gjn2uiSpdnhEN2/fK0mlDQMmiwFgGrbu2v9GCIwaOfa6tu7an1NFs0cQZOS2227T7bffPuH6vr4+Pf300y2sCMBcOHh4ZFrLy4AgyAlBAMyNvsGaerb0a/mmh9SzpV99g7VM329JR/u0lpdBkkGQ1R/OV77yFXV1denyyy/X/v31w8RvfOMbuvDCC3Xeeefp2muv1dGjR/XII49ox44d6u3t1fnnn69nnnlm3O0ATG50vL52eEShN8frswyD3tVdam+bd8Ky9rZ56l3dldl7Zi25IMjqD2fPnj3atm2bBgcHtX37du3evVuStH79eu3evVuPP/64zjnnHN1zzz265JJLtGbNGm3dulWPPfaYzjrrrHG3AzC5PMbr163q1Ob1K9XZ0S5L6uxo1+b1K0s7USxVqGtoqib7w5nNP+TDDz+sD3/4wzr99NMlSWvWrJEkPfnkk7rlllt0+PBhvfLKK1q9evW4vz/V7QC8Ka/x+nWrOku942+W3BFBln8447V3Xn/99brjjju0d+9e3XrrrROe/TzV7QC8qYrj9XlILgiy+sO59NJL9f3vf18jIyM6cuSIfvCDH0iSjhw5osWLF+vYsWO6//7739h+wYIFOnLkyBvPJ9oOwMSqOF6fh+SCIKs/nAsuuEAf/ehHdf755+vaa6/VBz7wAUnSl7/8ZV188cW64oordPbZZ7+x/YYNG7R161atWrVKzzzzzITbAZhYFcfr8+CIyLuGaenu7o7mG9Ps27dP55xzzpRfo8ynh0/3swKAJNneExHd461LbrJYqt5EDwDMRnJDQwCAE1UmCMo2xDUTKXxGAK1XiSCYP3++Dh06VOkd5ej9CObPn593KQAqphJzBEuXLtXQ0FDlr9U/eocyAJhLlQiCtrY27toFADOU2dCQ7TNs/8T2PttP2f7sONvY9tdsH7D9hO0LsqoHADC+LI8Ijkv6fEQ8anuBpD22fxQRY6+9fJWkFY2fiyXd2fgvAKBFMjsiiIjnI+LRxuMjkvZJam7eXyvpvqj7maQO24uzqgkAcLKWzBHYXiZplaSfN63qlPTcmOdDjWXPt6IuACiaPK58kHkQ2H6bpO9J+lxEvNy8epxfOakH1PZGSRsl6cwzz5zzGgGgCEbvlzJ6qfzR+6VIyjQMMj2PwHab6iFwf0RsH2eTIUlnjHm+VNLB5o0i4u6I6I6I7kWLFmVTLADkLI8b7UjZdg1Z0j2S9kXEVyfYbIekTzS6h94n6aWIYFgIQJLyutFOlkNDPZI+Lmmv7ccay74o6UxJioi7JO2UdLWkA5KOSrohw3oAoNCWdLSrNs5OP+sb7WQWBBHxnxp/DmDsNiHpxqxqAIAy6V3ddcIcgdSaG+1U4sxiACiq6XQBjS6vXNcQgOyV+WZLVTaTLqA87pdSiauPAikb3dnUDo8o9ObOpm+wlndpycurC2i6CAKg5Mqys0lRXl1A00UQACVXlp1Niibq9sm6C2i6CAKg5Mqys0lR7+outbfNO2FZK7qAposgAEquLDubFK1b1anN61eqs6NdltTZ0a7N61dq3apO9Q3W1LOlX8s3PaSeLf25zunQNQSUXF4th5ia8bqA8rqm0EQIAqAC8mg5xMxNNsGfx78jQ0MA0GJFm+AnCACgxYo2wU8QAECLFW2CnzkCAGixok3wEwQAkIMiTfAzNAQAiSMIACBxBAEAJI4gAIDEEQQAkDiCAAASRxAAQOIIAgBIHEEAAIkjCAAgcQQBACSOIACAxHHROQBJ6BusFeZqn0VDEACovKLdI7hoGBoCUHmT3SMYBAGABBTtHsFFQxAAqLyi3SO4aAgCIEF9gzX1bOnX8k0PqWdLv/oGa3mXlKmi3SO4aDILAtv32n7B9pMTrL/M9ku2H2v8fCmrWgC8aXTitHZ4RKE3J06rHAbrVnVq8/qV6uxolyV1drRr8/qVTBQ3ZNk19G1Jd0i6b5JtHo6IazKsAUCTySZOq7xjLNI9gosmsyOCiPippN9l9foAZoaJUzTLe47g/bYft/1D2++ZaCPbG20P2B4YHh5uZX1A5TBximZ5BsGjkt4dEedJ+kdJfRNtGBF3R0R3RHQvWrSoVfUBlcTEKZrldmZxRLw85vFO2/9ke2FEvJhXTUAKRsfJudxCeWR9eYzcgsD2uyT9NiLC9kWqH50cyqseICVlnDhN9VpBrbg8RmZBYPsBSZdJWmh7SNKtktokKSLukvQRSZ+2fVzSiKQNERFZ1QOgvFK+VlArurwyC4KIuO4U6+9Qvb0UACaVasur1Jour7y7hgDglFJueW1FlxdBAKDwUm55bUWXF0EAoPBSbnltxeUxuDENkLNUu2GmI/WW16y7vAgCIEcpd8NMVxlbXsuCIECSivItPOVuGBQHQYDkFOlbeMrdMLNRlCCvCiaLkZwi3b825W6YmWrV/RRSunkPQYDkFOlbeMrdMDPViiBP7eY9BAGSU6Rv4dw5a/paEeRFOmpsBeYIkJze1V0nzBFI+X4LpxtmepZ0tKs2zk5/LoO8SEeNrcARAZLDt/Bya8VwWpGOGluBIwIkiW/h5dLcJXTtn3bqJ78YzqxrqGhHjVkjCICElaENc7x23+/tqWV6FJfamcwEAZCoIp1PMZm8TrpL6aiRIAASleUOdi6PNFKbuM0Dk8VAorLawc51D35qE7d5IAiARGW1g53rHnxOusseQQAkKqsd7FwfadDumz3mCIBEZdUZk8UJXylN3OaBIAASlsUONrUe/CogCADMqdR68KuAIAAw5xjKKRcmiwEgcQQBACSOIACAxBEEAJA4ggAAEkcQAEDiCAIASBznEQAFVYabxqAaCAKggMpy0xhUQ2ZDQ7bvtf2C7ScnWG/bX7N9wPYTti/IqhagbOb6Us7AZE4ZBLZvsv32Gbz2tyVdOcn6qyStaPxslHTnDN4DqCTuyoVWmsoRwbsk7bb9XdtX2vZUXjgifirpd5NsslbSfVH3M0kdthdP5bWBquOuXGilUwZBRNyi+rf2eyRdL+mXtv/W9lmzfO9OSc+NeT7UWHYS2xttD9geGB4enuXbAsXHXbnQSlOaI4iIkPSbxs9xSW+X9KDtv5vFe493ZBETvP/dEdEdEd2LFi2axVsC5cBdudBKp+wasv0ZSZ+U9KKkb0rqjYhjtt8i6ZeSvjDD9x6SdMaY50slHZzhawEnKXv7JZdyRqtMpX10oaT1EfGrsQsj4ve2r5nFe++QdJPtbZIulvRSRDw/i9cD3kD7JTB1pwyCiPjSJOv2TbTO9gOSLpO00PaQpFsltTV+7y5JOyVdLemApKOSbphO4cBkJmu/JAiAE2V2QllEXHeK9SHpxqzeH2mj/RKYOq41hEqi/RKYOoIApdM3WFPPln4t3/SQerb0q2+wdtI2tF8CU8e1hlAqU50EHn1c5q4hoFVcH6ovj+7u7hgYGMi7DOSkZ0u/auOM83d2tOu/Nn1wwt8reyspMFu290RE93jrOCJAqcxkEphWUmByzBGgVGYyCcyVPIHJEQQolZlMAtNKCkyOIECpzOQaPLSSApNjjgClM91r8PSu7jphjkCilRQYiyBA5dFKCkyOIEASpnsUMZftprSuougIAqDJXLab0rqKMmCyGGgyl+2mtK6iDAgCoMlctpvSuooyIAiAJnPZbkrrKsqAIACazOWVS7kKKsqAyWKgyVy2m9K6ijLg6qMAkIDJrj7K0BAAJI6hISBHnGyGIiAIgJxwshmKgqEhICecbIaiIAiAnHCyGYqCIABywslmKAqCAMgJJ5uhKJgsBnLCyWYoCoIAyNF075MAZIGhIQBIHEEAAIljaGgWOCsUQBUQBDPEWaEAqoIgmKLmb/9HXzs+4VmhBAGAMsl0jsD2lbb32z5ge9M46y+z/ZLtxxo/X8qynpka/fZfOzyiUP3b//8ePTbutpwVCqBsMjsisD1P0tclXSFpSNJu2zsi4ummTR+OiGuyqmMujHdNmIlwViiAssnyiOAiSQci4tmIeE3SNklrM3y/zEz1Wz5nhQIooyyDoFPSc2OeDzWWNXu/7cdt/9D2e8Z7IdsbbQ/YHhgeHs6i1klN9C2/o71NnR3tsqTOjnZtXr+S+QFMqG+wpp4t/Vq+6SH1bOlX32At75IASdlOFnucZc33xXxU0rsj4hXbV0vqk7TipF+KuFvS3VL9VpVzXOcp9a7uOqFDSKp/+79tzXvY8WNK6DJDkWV5RDAk6Ywxz5dKOjh2g4h4OSJeaTzeKanN9sIMa5qRdas6tXn9Sr79Y8a49wCKLMsjgt2SVtheLqkmaYOkj43dwPa7JP02IsL2RaoH06EMa5oxrgmD2eDeAyiyzIIgIo7bvknSLknzJN0bEU/Z/lRj/V2SPiLp07aPSxqRtCEiWj70A2RtSUe7auPs9OkyQxFkekJZY7hnZ9Oyu8Y8vkPSHVnWABTBRPNMdJmhCDizGGgB7j2AIiMIgBZhnglFxWWoASBxBAEAJI4gAIDEEQQAkDiCAAASRxAAQOIIAgBIHEEAAIkjCAAgcQQBACSOIACAxBEEAJA4ggAAEkcQAEDiCAIASBxBAACJIwgAIHEEAQAkjiAAgMQRBACQOIIAABJHEABA4ggCAEgcQQAAiSMIACBxBAEAJI4gAIDEEQQAkDiCAAASRxAAQOIyDQLbV9reb/uA7U3jrLftrzXWP2H7gizrAQCcLLMgsD1P0tclXSXpXEnX2T63abOrJK1o/GyUdGdW9QAAxpflEcFFkg5ExLMR8ZqkbZLWNm2zVtJ9UfczSR22F2dYEwCgSZZB0CnpuTHPhxrLpruNbG+0PWB7YHh4eM4LBYCUZRkEHmdZzGAbRcTdEdEdEd2LFi2ak+IAAHVZBsGQpDPGPF8q6eAMtgEAZOi0DF97t6QVtpdLqknaIOljTdvskHST7W2SLpb0UkQ8P9eF9A3WtHXXfh08PKIlHe3qXd2ldatOGoECgCRlFgQRcdz2TZJ2SZon6d6IeMr2pxrr75K0U9LVkg5IOirphrmuo2+wppu379XIsdclSbXDI7p5+15JIgwAQNkeESgidqq+sx+77K4xj0PSjVnWsHXX/jdCYNTIsde1ddd+ggAAlMCZxQcPj0xrOQCkpvJBsKSjfVrLASA1lQ+C3tVdam+bd8Ky9rZ56l3d1bIa+gZr6tnSr+WbHlLPln71DdZa9t4AcCqZzhEUweg8QF5dQ0xWAyi6ygeBVN/h5rXTZbIaQNFVfmgob0xWAyg6giBjTFYDKDqCIGNTmaxmMhlAnpKYI8jTqSarmUwGkDeCoAUmm6xmMhlA3hgayhmTyQDyRhDkjMlkAHkjCHJWhDOfAaSNOYKc5X3mMwAQBAWQ55nPAMDQEAAkjiAAgMQRBACQOIIAABJHEABA4ly/f3x52B6W9KtTbLZQ0ostKKeoUv78KX92Ke3Pn/Jnl079+d8dEYvGW1G6IJgK2wMR0Z13HXlJ+fOn/NmltD9/yp9dmt3nZ2gIABJHEABA4qoaBHfnXUDOUv78KX92Ke3Pn/Jnl2bx+Ss5RwAAmLqqHhEAAKaIIACAxFUqCGxfaXu/7QO2N+VdTyvZPsP2T2zvs/2U7c/mXVOr2Z5ne9D2v+VdS6vZ7rD9oO1fNP4G3p93Ta1i+68bf/NP2n7A9vy8a8qS7Xttv2D7yTHL3mH7R7Z/2fjv26fzmpUJAtvzJH1d0lWSzpV0ne1z862qpY5L+nxEnCPpfZJuTOzzS9JnJe3Lu4ic/IOkf4+IsyWdp0T+d7DdKekzkroj4r2S5knakG9Vmfu2pCublm2S9OOIWCHpx43nU1aZIJB0kaQDEfFsRLwmaZuktTnX1DIR8XxEPNp4fET1HUEyNzmwvVTSn0n6Zt61tJrtP5R0qaR7JCkiXouIw7kW1VqnSWq3fZqk0yUdzLmeTEXETyX9rmnxWknfaTz+jqR103nNKgVBp6TnxjwfUkI7wrFsL5O0StLPcy6llf5e0hck/T7nOvLwJ5KGJX2rMTT2TdtvzbuoVoiImqTbJf1a0vOSXoqI/8i3qlz8cUQ8L9W/FEr6o+n8cpWCwOMsS6431vbbJH1P0uci4uW862kF29dIeiEi9uRdS05Ok3SBpDsjYpWk/9M0hwbKqjEWvlbScklLJL3V9l/kW1X5VCkIhiSdMeb5UlX8ELGZ7TbVQ+D+iNiedz0t1CNpje3/UX1I8IO2/znfklpqSNJQRIweAT6oejCk4HJJ/x0RwxFxTNJ2SZfkXFMefmt7sSQ1/vvCdH65SkGwW9IK28tt/4HqE0Y7cq6pZWxb9THifRHx1bzraaWIuDkilkbEMtX/3fsjIplvhRHxG0nP2e5qLPqQpKdzLKmVfi3pfbZPb/x/4ENKZKK8yQ5Jn2w8/qSkf53OL1fm5vURcdz2TZJ2qd45cG9EPJVzWa3UI+njkvbafqyx7IsRsTO/ktBCfyXp/saXoGcl3ZBzPS0RET+3/aCkR1XvnBtUxS81YfsBSZdJWmh7SNKtkrZI+q7tv1Q9HP98Wq/JJSYAIG1VGhoCAMwAQQAAiSMIACBxBAEAJI4gAIDEEQQAkDiCAAASRxAAs2T7QttP2J5v+62Na+O/N++6gKnihDJgDtj+G0nzJbWrft2fzTmXBEwZQQDMgcalHXZLelXSJRHxes4lAVPG0BAwN94h6W2SFqh+ZACUBkcEwBywvUP1S2Avl7Q4Im7KuSRgyipz9VEgL7Y/Iel4RPxL497Zj9j+YET0510bMBUcEQBA4pgjAIDEEQQAkDiCAAASRxAAQOIIAgBIHEEAAIkjCAAgcf8PCMT7CsmPdNQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate and plot data\n",
    "np.random.seed(0)\n",
    "n = 30\n",
    "\n",
    "def sample_data(num_points):\n",
    "    x = rand(num_points)*10\n",
    "    y = .2*x + .3*randn(num_points)\n",
    "    return x,y\n",
    "\n",
    "x,y = sample_data(n)\n",
    "\n",
    "\"\"\"function to plot the above data\"\"\"\n",
    "def plotdata(x=x,y=y, margin=.05):\n",
    "    plt.scatter(x,y, label=\"data\")\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.xlabel(\"x\") \n",
    "    plt.ylabel(\"y\")     \n",
    "    range_y = np.max(y) - np.min(y)\n",
    "    range_x = np.max(x) - np.min(x)\n",
    "    plt.ylim([np.min(y)-margin*range_y,np.max(y)+margin*range_y])\n",
    "    plt.xlim([np.min(x)-margin*range_x,np.max(x)+margin*range_x])\n",
    "    \n",
    "plotdata()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUz0lEQVR4nO3df2xd5X3H8c9nxl0MzZS2eC1xYM4Yc4JgSegFWqKykaYLXVkIiSqo1K50laJNZaVTlS5Z0dapY0Qzmro/qklpYXRqWlZB4qalbUobUH8wWOw6NIGQlbD+sEMXN5PXMJyRhO/+8DUEc21fJ/fc5977vF9SFPv4+Jzvzb05n3Oe5znncUQIAJCfX0ldAAAgDQIAADJFAABApggAAMgUAQAAmTordQGzce6550Z3d3fqMgCgqQwMDPwiIjonL2+qAOju7lZ/f3/qMgCgqdj+SaXlNAEBQKYIAADIFAEAAJlqqj6ASo4fP66hoSEdO3YsdSmFmTNnjhYsWKD29vbUpQBoIU0fAENDQ5o7d666u7tlO3U5NRcROnLkiIaGhrRw4cLU5QBoIUkDwPaPJR2VdFLSiYgozXYbx44da9mDvyTZ1hve8AaNjIykLgVAAn2Dw+rdeUCHRsc0f16HNqzq0ZplXTXZdiNcAVwTEb84kw206sF/Qqu/PgCV9Q0Oa9O2vRo7flKSNDw6pk3b9kpSTUKATmAAaFC9Ow+8dPCfMHb8pHp3HqjJ9lMHQEj6pu0B2+srrWB7ve1+2/3N0AzyiU98QnfeeeeUP+/r69OTTz5Zx4oANKtDo2OzWj5bqQNgeURcJumdkj5k++rJK0TElogoRUSps/NVdzI3HQIAQLXmz+uY1fLZShoAEXGo/PdhSdslXVH0PvsGh7V88y4t3PiAlm/epb7B4TPe5u23366enh6tXLlSBw6MX5p95jOf0eWXX64lS5Zo3bp1ev755/XII49ox44d2rBhg5YuXaqDBw9WXA8AJGnDqh51tLe9YllHe5s2rOqpyfaTBYDtc2zPnfha0u9L2lfkPic6VIZHxxR6uUPlTEJgYGBA9957rwYHB7Vt2zbt3r1bkrR27Vrt3r1bjz/+uBYvXqy77rpLV111lVavXq3e3l7t2bNHF154YcX1AEAa7+i9Y+2l6prXIUvqmtehO9Ze2hKjgN4oaXt5hMtZkr4QEd8ocofTdaic7j/od7/7Xd1www06++yzJUmrV6+WJO3bt0+33XabRkdH9dxzz2nVqlUVf7/a9QDkac2yrpod8CdLFgAR8YykJfXcZ1EdKpWGad58883q6+vTkiVLdM899+jhhx+u+LvVrgcAtZa6E7iuiuhQufrqq7V9+3aNjY3p6NGj+spXviJJOnr0qM477zwdP35cW7dufWn9uXPn6ujRoy99P9V6AFC0rAKgiA6Vyy67TDfeeKOWLl2qdevW6W1ve5sk6ZOf/KSuvPJKveMd79CiRYteWv+mm25Sb2+vli1bpoMHD065HgAUzRGRuoaqlUqlmDwhzP79+7V48eKqt1HkbdVFmu3rBIAJtgcqPWqnER4FUVdFdqgAQDPJqgkIAPCylgiAZmrGOh2t/voApNH0ATBnzhwdOXKkZQ+SE/MBzJkzJ3UpAFpM0/cBLFiwQENDQy39vPyJGcEAoJaaPgDa29uZKQsATkPTNwEBAE4PAQAAmSIAACBTBAAAZIoAAIBMEQAAkCkCAAAyRQAAQKYIAADIFAEAAJlKHgC222wP2v5q6loAICfJA0DSrZL2py4CAHKTNABsL5D0LkmfTVkHAOQo9RXApyR9TNKLU61ge73tftv9rfzIZwCot2QBYPs6SYcjYmC69SJiS0SUIqLU2dlZp+oAoPWlvAJYLmm17R9LulfSCtufT1gPAGQlWQBExKaIWBAR3ZJukrQrIt6bqh4AyE3qPgAAQCINMSVkRDws6eHEZQBAVrgCAIBMEQAAkCkCAAAyRQAAQKYIAADIFAEAAJkiAAAgUwQAAGSKAACATBEAAJApAgAAMkUAAECmCAAAyBQBAACZIgAAIFMEAABkigAAgEwRAACQKQIAADJFAABAppIFgO05tv/d9uO2n7D9N6lqAYAcnZVw3/8naUVEPGe7XdL3bH89Ih5NWBMAZCNZAERESHqu/G17+U+kqgcAcpO0D8B2m+09kg5LejAiHquwznrb/bb7R0ZG6l4jALSqpAEQEScjYqmkBZKusH1JhXW2REQpIkqdnZ11rxEAWlVDjAKKiFFJD0u6Nm0lAJCPlKOAOm3PK3/dIWmlpKdS1QMAuUk5Cug8SZ+z3abxIPpSRHw1YT0AkJWUo4B+KGlZqv0DQO4aog8AAFB/BAAAZIoAAIBMEQAAkCkCAAAyRQAAQKYIAADIFAEAAJkiAAAgUwQAAGSKAACATBEAAJApAgAAMpXycdAA0HT6BofVu/OADo2Oaf68Dm1Y1aM1y7pSl3VaCAAATaveB+O+wWFt2rZXY8dPSpKGR8e0adteSWrKEKAJCEBTmjgYD4+OKfTywbhvcLiwffbuPPDSwX/C2PGT6t15oLB9FokAANCUUhyMD42OzWp5oyMAADSlFAfj+fM6ZrW80REAAJpSioPxhlU96mhve8WyjvY2bVjVU9g+i5QsAGyfb/sh2/ttP2H71lS1AGg+KQ7Ga5Z16Y61l6prXocsqWteh+5Ye2lTdgBLaUcBnZD00Yj4ge25kgZsPxgRTyasCUCTmDjo1ntI5pplXYXto96jmpIFQEQ8K+nZ8tdHbe+X1CWJAABQlSIPxvWWYohpQ9wHYLtb0jJJj1X42XpJ6yXpggsuqG9hQBNopRuTWs1s3pvpRjUV9X4m7wS2/VpJ90v6SET8cvLPI2JLRJQiotTZ2Vn/AoEGlmIsPKoz2/cmxaimpFcAtts1fvDfGhHbUtYCNKMUZ42ozkz3KUy+Mpg/r0PDFQ72RY5qSjkKyJLukrQ/Iv4hVR1AM2u1G5NayVTvwcSVwOQrg2sWddZ9VFPKJqDlkt4naYXtPeU/f5CwHqDptNqNSa1kqvegza54ZfDQUyN1H2KachTQ9yQ51f6BVrBhVc8rRo5IzX1jUiuZ6r2ZfPCfcGh0rO6jmpJ3AgM4fa12Y1Irmeq96Wqgq7aGGAYK4PS10lj4VjPVe9MoV20EAADUUao7mCshAAC0vEa7Wa5RrtoIAAAtrdVm8aolAgDITKOdDReNm+WmRgAAGcnxbJib5aY24zBQ27fYfl09igFQrGad07ZvcFjLN+/Swo0PaPnmXbN61hE3y02tmvsA3iRpt+0v2b62/AgHAE2oGc+Gz/SBd602i1ctzRgAEXGbpIs0/tyemyX9yPbf2b6w4NqAlncmZ7anoxnPhs/0qoWb5aZWVR9ARITtn0v6ucZn8nqdpPvKM3h9rMgCgVprlE7QFO3xzfjoiOkeqrZ8866q3r/ZDrtslM9I0arpA/iw7QFJfy/p+5IujYg/lfRmSesKrg+oqUZ6fn6K9vhmPBue7uqkiPevkT4jRavmCuBcSWsj4ienLoyIF21fV0xZQDEaaUhgqvb4RrkJqVrXLOrU1kd/qpji57V+/xrpM1K0GQMgIv5qmp/tr205QLEaqRM0xQQgpyqqmaOW2+0bHNb9A8NTHvwn1PL9a6TPSNF4Giiy0kidoClHpxTVzFHr7VY6G6+klu9fI31GikYAICuNNCQwZXt8Uf0Ptd5uNWfdtX7/GukzUjTuBEZWGulJjBP1pNh3Uc0ctd7uVM1kbbZejCjk/Wu0z0iRCABkp9k6QYtQVP9Drbc71bDVoq+UcvmM0AQENKCibxArqpmj1tttxmGrzSTpFYDtuyVdJ+lwRFySshagUdTjBrGimjmK2G4uZ+MpOGKmAVYF7ty+WtJzkv6lmgAolUrR399ffGFAQss376rYjNI1r0Pf37giQUVodrYHIqI0eXnSK4CI+I7t7pQ1oPU0+238OY1DR1oN3wdge73tftv9IyMjqctBg2uF2/hzGoeOtBo+ACJiS0SUIqLU2dmZuhwkNlPn6FTj0D/yr3vq8rTNWshpHDrSYhgomkY1naPTNZM0y+xXOY1DR1oEAJpGNQ/pmmoc+lTrNypGvqAekjYB2f6ipH+T1GN7yPYHU9aDxlZN52il5pNqtyPVdvx9vSd7AWYr9Sig96TcP5pLNXeZntp8MtWVwFSdqbUcf5/j5OtoPg3fCQxMqLZzdM2yLn1/4wp96sals+pMreWDzJp18nXkhT4ANI3Zdo7Odv1ajr+vZlvNfr8Cmh8BgKYy287R2axfyweZzbQtmojQCGgCAspqOf5+pm3RRIRGwBUAUFbL8fczbYvHPaAREAAzoJ02L7Ucfz/dtlLPBwxIBMC0aKdFrUw+kbhmUafuHxh+1UQnPO4B9UQfwDRop0UtVHpA3f0Dw1r35i4mOkFSXAFMg3Za1MJUJxIPPTXC8/2RFFcA0+CxvKgFTiTQqAiAafBYXtQCJxJoVATANJiQGrXAiQQaFX0AM+CxvDhTPN8fjYoAAOqAEwk0IpqAACBTBAAAZIoAAIBMEQAAkCkCAAAylXpS+GttH7D9tO2NKWsBgNwkCwDbbZI+Lemdki6W9B7bF6eqBwByk/IK4ApJT0fEMxHxgqR7JV2fsB4AyErKG8G6JP3slO+HJF05eSXb6yWtl6QLLrhg1jthQhcAqCxlALjCsnjVgogtkrZIUqlUetXPp5N6QhfCB0AjS9kENCTp/FO+XyDpUC13kHJCl0qTgGzatld9g8OF7xsAqpEyAHZLusj2QtuvkXSTpB213EHK57BXEz59g8NavnmXFm58QMs37yIcANRVsiagiDhh+xZJOyW1Sbo7Ip6o5T5STrw9U/ikbp4CgKT3AUTE1yLityPiwoi4vdbbT/kc9pkmAWG+YQCptfSdwCkndJkpfJgmEEBqLT8fQKrnsM80CUjK5ikAkDIIgJSmC58Nq3pe0QcgMU0ggPoiABJhmkAAqREACTFNIICUWroTGAAwNQIAADJFAABApggAAMgUAQAAmSIAACBTBAAAZIoAAIBMEQAAkCkCAAAyRQAAQKYIAADIFAEAAJkiAAAgUwQAAGQqSQDYfrftJ2y/aLuUogYAyF2qK4B9ktZK+k6i/QNA9pLMCBYR+yXJdordAwDUBH0Attfb7rfdPzIykrocAGgZhV0B2P6WpDdV+NHHI+LL1W4nIrZI2iJJpVIpalQeAGSvsACIiJVFbRsAcOYavgkIAFCMVMNAb7A9JOmtkh6wvTNFHQCQs1SjgLZL2p5i3wCAcTQBAUCmCAAAyBQBAACZIgAAIFMEAABkigAAgEwRAACQKQIAADJFAABApggAAMgUAQAAmSIAACBTBAAAZIoAAIBMEQAAkCkCAAAyRQAAQKYIAADIFAEAAJlKMiew7V5JfyjpBUkHJX0gIkZT1HI6+gaH1bvzgA6Njmn+vA5tWNWjNcu6UpcFALOS6grgQUmXRMTvSPoPSZsS1TFrfYPD2rRtr4ZHxxSShkfHtGnbXvUNDqcuDQBmJUkARMQ3I+JE+dtHJS1IUcfp6N15QGPHT75i2djxk+rdeSBRRQBwehqhD+CPJX19qh/aXm+733b/yMhIHcuq7NDo2KyWA0CjKiwAbH/L9r4Kf64/ZZ2PSzohaetU24mILRFRiohSZ2dnUeVWbf68jlktB4BGVVgncESsnO7ntt8v6TpJb4+IKKqOWtuwqkebtu19RTNQR3ubNqzqSVgVAMxeqlFA10r6C0m/GxHPp6jhdE2M9mEUEIBm5xQn37aflvSrko6UFz0aEX8y0++VSqXo7+8vtDYAaDW2ByKiNHl5kiuAiPitFPsFALysEUYBAQASIAAAIFMEAABkigAAgEwlGQV0umyPSPrJDKudK+kXdSinEeX82qW8X3/Or13K+/VX89p/IyJedSdtUwVANWz3VxrulIOcX7uU9+vP+bVLeb/+M3ntNAEBQKYIAADIVCsGwJbUBSSU82uX8n79Ob92Ke/Xf9qvveX6AAAA1WnFKwAAQBUIAADIVMsEgO1rbR+w/bTtjanrqSfb59t+yPZ+20/YvjV1TfVmu832oO2vpq6l3mzPs32f7afKn4G3pq6pXmz/efkzv8/2F23PSV1TkWzfbfuw7X2nLHu97Qdt/6j89+uq3V5LBIDtNkmflvROSRdLeo/ti9NWVVcnJH00IhZLeoukD2X2+iXpVkn7UxeRyD9K+kZELJK0RJn8O9jukvRhSaWIuERSm6Sb0lZVuHskXTtp2UZJ346IiyR9u/x9VVoiACRdIenpiHgmIl6QdK+k62f4nZYREc9GxA/KXx/V+AEgmxlqbC+Q9C5Jn01dS73Z/jVJV0u6S5Ii4oWIGE1aVH2dJanD9lmSzpZ0KHE9hYqI70j670mLr5f0ufLXn5O0ptrttUoAdEn62SnfDymjA+CpbHdLWibpscSl1NOnJH1M0ouJ60jhNyWNSPrnchPYZ22fk7qoeoiIYUl3SvqppGcl/U9EfDNtVUm8MSKelcZPBiX9erW/2CoB4ArLshvfavu1ku6X9JGI+GXqeurB9nWSDkfEQOpaEjlL0mWS/ikilkn6X82iCaCZldu6r5e0UNJ8SefYfm/aqppLqwTAkKTzT/l+gVr8UnAy2+0aP/hvjYhtqeupo+WSVtv+scab/lbY/nzakupqSNJQRExc8d2n8UDIwUpJ/xkRIxFxXNI2SVclrimF/7J9niSV/z5c7S+2SgDslnSR7YW2X6PxjqAdiWuqG9vWeBvw/oj4h9T11FNEbIqIBRHRrfH3fVdEZHMWGBE/l/Qz2z3lRW+X9GTCkurpp5LeYvvs8v+BtyuTDvBJdkh6f/nr90v6crW/mGRO4FqLiBO2b5G0U+MjAe6OiCcSl1VPyyW9T9Je23vKy/4yIr6WriTU0Z9J2lo++XlG0gcS11MXEfGY7fsk/UDjI+EG1eKPhLD9RUm/J+lc20OS/lrSZklfsv1BjYfiu6veHo+CAIA8tUoTEABglggAAMgUAQAAmSIAACBTBAAAZIoAAIBMEQAAkCkCADgDti+3/UPbc2yfU342/SWp6wKqwY1gwBmy/beS5kjq0Phzee5IXBJQFQIAOEPlRzDslnRM0lURcTJxSUBVaAICztzrJb1W0lyNXwkATYErAOAM2d6h8UdRL5R0XkTckrgkoCot8TRQIBXbfyTpRER8oTw39SO2V0TErtS1ATPhCgAAMkUfAABkigAAgEwRAACQKQIAADJFAABApggAAMgUAQAAmfp/uLe8h06qLCUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# mess up the data\n",
    "outliers = True\n",
    "if outliers:\n",
    "    y[20] = 5\n",
    "    y[15] = -2\n",
    "\n",
    "plotdata()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Least squares will fit inliers _______ than Huber regression \n",
    "* A) better \n",
    "* B) worse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Least squares will fit outliers _______ than Huber regression \n",
    "* A) better \n",
    "* B) worse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Least squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.30576127701199385, -0.6102904304198062]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# X = vcat(x ones(length(x))]\n",
    "X = x.reshape(-1,1)\n",
    "lin = LinearRegression(fit_intercept=True).fit(X,y)\n",
    "w_quad = [lin.coef_[0], lin.intercept_]\n",
    "w_quad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlqklEQVR4nO3de3RV5Z3/8feTEEggSBAwSIKCxEEsCIHDpaRCAtmbOuNQxc7ItJ1ROy5aa3V0KhWK/fUyVWmZX9d0dWbN+mntz85PZ1KnKrW2U5IQjpGbkhBqkIsSQEzwwi2QQAK5PL8/Eo6AuZyQc84+yf681mJJdnae/f16wvN99rOfvbex1iIiIv6T4HUAIiLiDRUAERGfUgEQEfEpFQAREZ9SARAR8akBXgfQEyNHjrTjxo3rcp/Tp08zZMiQ2AQUZ/ycO/g7fz/nDv7OP5zcy8vLj1prR126vU8VgHHjxlFWVtblPsFgkNzc3NgEFGf8nDv4O38/5w7+zj+c3I0x73W0XVNAIiI+pQIgIuJTKgAiIj7Vp64BdKSpqYnq6moaGxsBGDZsGLt37/Y4qshKTk4mMzOTpKQkr0MRkX6kzxeA6upqhg4dyrhx4zDGUFdXx9ChQ70OK2KstRw7dozq6mrGjx/vdTgi0o94WgCMMQeBOqAFaLbWBnraRmNjY6jz74+MMYwYMYIjR454HYqIeGBtRQ1r1u3lcG0DY9JSWL5oIrdlZ0Sk7Xg4A8iz1h7tTQP9tfM/r7/nJyIdW1tRw8qXKmloagGgpraBlS9VAkSkCOgisIhInFqzbm+o8z+voamFNev2RqR94+X7AIwxB4ATgAX+j7X2qQ72WQYsA0hPT59RUFBw0feHDRtGVlZW6OuWlhYSExOjGXaXnnjiCVJTU3nwwQc7/P6rr75KVlYWN9xwQ4/a3bdvHydPnuxyn/r6elJTU3vUbn/i5/z9nDv03/wrazr/Nz8lYxgQXu55eXnlHU2xez0FlGOtPWyMuQooMsbssdaWXrhDe1F4CiAQCNhL73jbvXv3RRd9vb4IPGjQIAYNGtRpDOvWrSMpKYmZM2f2qN3k5GSys7O73MfPd0OCv/P3c+7Qf/NftbqEmtqGT23PSEvhgS/nAr3L3dMpIGvt4fb/fgy8DMyK9jHXVtSQs7qE8St+T87qEtZW1PS6zccff5yJEyeSn5/P3r1tp2ZPP/00M2fOZOrUqdxxxx2cOXOGzZs388orr7B8+XKmTZtGVVVVh/uJiAAsXzSRlKSLZzRSkhJZvmhiRNr3rAAYY4YYY4ae/zvgAjujeczzF1RqahuwfHJBpTdFoLy8nIKCAioqKnjppZfYtm0bAEuWLGHbtm386U9/YtKkSTzzzDPMnTuXxYsXs2bNGnbs2MGECRM63E9EBNou9D65ZAoZaSkY2kb+Ty6Z0i9WAaUDL7evcBkA/Ke19o/RPGBXF1Qu93/o66+/zu23387gwYMBWLx4MQA7d+7kscceo7a2lvr6ehYtWtThz4e7n4j4023ZGRHr8C/lWQGw1u4HpsbymIc7mEvranu4Olqmeffdd7N27VqmTp3Ks88+SzAY7PBnw91PRCTSfLUMdExaSo+2h2PevHm8/PLLNDQ0UFdXx+9+9zug7WL01VdfTVNTE88//3xo/6FDh1JXVxf6urP9RESizVcFIBoXVKZPn86dd97JtGnTuOOOO7j55psB+Kd/+idmz56N4zgXLflcunQpa9asITs7m6qqqk73ExGJNq+XgcbU+Xm0SN9WvWrVKlatWvWp7ffdd9+ntuXk5LBr166L9uloPxGRaPNVAYDoXlAREelLfDUFJCIin+gXBcDLx1nEQn/PT0S80ecLQHJyMseOHeu3neT59wEkJyd7HYqI9DN9/hpAZmYm1dXVoeflNzY29rvO8vwbwUREIqnPF4CkpKSL3pQVDAa7fWiaiIj0gykgERG5PCoAIiI+pQIgIuJTKgAiIj6lAiAi4lMqACIiPqUCICLiUyoAIiI+pQIgIuJTKgAiIj7leQEwxiQaYyqMMa96HYuIiJ94XgCAfwB2ex2EiIjfeFoAjDGZwF8Av/AyDhERPzJePkffGPMb4ElgKPCItfbWDvZZBiwDSE9Pn1FQUNBlm/X19aSmpkYh2vjn59zB3/n7OXfwd/7h5J6Xl1durQ1cut2zx0EbY24FPrbWlhtjcjvbz1r7FPAUQCAQsLm5ne4KtD0Ourt9+is/5w7+zt/PuYO/8+9N7l5OAeUAi40xB4ECYIEx5jkP4xER8RXPCoC1dqW1NtNaOw5YCpRYa7/iVTwiIn4TD6uARETEA3HxSkhrbRAIehyGiIiv6AxARMSnVABERHxKBUBExKdUAEREfEoFQETEp1QARER8SgVARMSnVABERHxKBUBExKdUAEREfEoFQETEp1QARER8SgVARMSnVABERHxKBUBExKdUAEREfEoFQETEp1QARER8SgVARMSnVABERHzKswJgjEk2xrxpjPmTMeZtY8wPvIpFRMSPBnh47LPAAmttvTEmCdhojPkfa+1WD2MSEfENzwqAtdYC9e1fJrX/sV7FIyLiN6atH/bo4MYkAuVAFvBv1tpHO9hnGbAMID09fUZBQUGXbdbX15OamhqFaOOfn3MHf+fv59zB3/mHk3teXl65tTbwqW9Yaz3/A6QBG4DJXe03Y8YM250NGzZ0u09/5efcrfV3/n7O3Vp/5x9O7kCZ7aBPjYtVQNbaWiAIfN7bSERE/MPLVUCjjDFp7X9PAfKBPV7FIyLiN16uAroa+FX7dYAE4AVr7asexiMi4itergJ6C8j26vgiIn4XF9cAREQk9lQARER8SgVARMSnVABERHxKBUBExKdUAEREfMrL+wBERCRK9u/fT1FREYWFhZ3uowIgItIP1NbWsmHDBgoLCykqKqKqqgqAsWPHdvozmgISEemDmpub2bRpE88++yxz585lxIgRLFmyhOeee44bb7yRn//85+zZs4f33nuv0zZ0BiAi0gdYa6mqqgqN8EtKSjh16hQJCQkEAgG+853v4Lous2fPZuDAgWG1qQIgIhKnTpw4wfr160Nz+QcPHgTg2muvZenSpTiOw8CBA1m8ePFlta8CICISJ5qamti6dWtolL9t2zZaW1sZOnQoCxYsYPny5TiOQ1ZWFsYYAILB4GUfTwVARMQj1lreeeed0Ah/w4YN1NfXk5CQwOzZs3nsscdwXZdZs2aRlJQU8eOrAIiI9MDaihrWrNvL4doGxqSlsHzRRG7Lzgj7548dO8b69etDo/xDhw4BcN111/GVr3wFx3FYsGABaWlpUcrgEyoAItJnne+Ml46tY9Xqkh53xpdzvJUvVdLQ1AJATW0DK1+qBOj0uOfOnWPz5s2hUX55eTnWWoYNG8aCBQtYuXIljuMwYcKEqMXdGRUAEemTLuqMx4bXGffWmnV7Q53/eQ1NLaxZtzd0TGste/bsCY3wg8Egp0+fJjExkTlz5vD9738fx3GYOXMmAwZ42wWrAIhInxROZxxph2sbOtz+/uEPKSgooLCwkMLCQmpqagC4/vrrueuuu3Bdl9zcXIYNGxaVuC6XCoCI9EmddcadbY+EMWkp1NQ2YJubaKzZReOBChoPVnDuoyr+5ueQlpZGfn4+juPgOA7jx4+PWiyRoAIgIn3S+c64o+2RZq1l165dTDxSyp9++3vOvFeJbT4LCYmkZEziS99YzoN/dweBQIDExMSIHz9aPCsAxpixwH8Ao4FW4Clr7c+8ikdE+pbliyZedEEWICUpkeWLJkak/Y8++oji4mKKioooKiri8OHDAGSMm8DQmX9O65jJjJ8yixWLs6N64TmavDwDaAa+Za3dbowZCpQbY4qstbs8jElE+ojzne6adXuBOjIuY0nmhRobG9m4cWPo4u2OHTsAuPLKK8nPz8d1XRzH4ZprrolQBp/W2yWmPeVZAbDWfgB80P73OmPMbiADUAEQkbDclp3BbdkZBINBHvhybo9+1lpLZWVlaHlmaWkpjY2NJCUlkZOTw+OPP47rumRnZ8dkWudylpj2lrHWRqXhHgVhzDigFJhsrT11yfeWAcsA0tPTZxQUFHTZVn19PampqVGKNL75OXfwb/61DU00N57hgzMwMDGB9GHJpKVE/q7ReBbuZ3/8+HHKysooKyujvLyc48ePA23P1gkEAgQCAaZOnUpKSmSuI9Q2NPHRyUbOtbR2+9ns/bCOcy2tn9o+MDGBiaOHdnqMcHLPy8srt9YGLt3ueQEwxqQCrwGPW2tf6mrfQCBgy8rKumwvGAySm5sbuQD7ED/nDv7M//yo8Rs3nOV/V7ad0KckJfLkkil9dl76cnT22Z85c4bXX389NMqvrGwbUY8cOTK0UsdxHDIzMyMe06Ujeuj6sxm/4vd01Bsb4MDqv+j0OOH83htjOiwAnq4CMsYkAS8Cz3fX+YvIp3mxFj6etba28tZbb4Xm8V9//XXOnj3LwIED+dznPsfq1atxHIdp06aRkBDd16F09dmc//6Fc/2xXNV0npergAzwDLDbWvtTr+IQ6cu8WAsfbw4fPswf//hHnn76aYqLi/n4448BmDx5Mt/4xjdwXZd58+YxePDg2MbVyWdwfm7/0rn+O2Zk8GJ5TdRWNXXEyzOAHOBvgUpjzI72bd+x1v7Bu5BE+hYvRo1eO336NKWlpaFpnbfffhuAq666CsdxcF2X/Px8xowZ42mcnX02icZ0eGawYc8RnlwyxTergDbSNr0lIpfp/Fr4tlXVbaI9aoy11tZWduzYEZrW2bhxI+fOnWPQoEHMmzePu+66i+HDh/PVr3416tM6PdHZfQqXdv7nHa5tCK1qihXdCSzSh53vLD7aux0DMRk1xkJ1dXVohF9cXMzRo0cBuOmmm3jwwQdxHIebb745tFonGAzGVecPF9+ncOGIfs26vXFz1qYCINLH3ZadQfDkuxxYnet1KJetvr6e1157LTTK3717NwCjR4/mlltuCU3rjB492uNIe6azEX0072DuCRUAEYm5lpYWtm/fHhrlb968maamJpKTk5k/fz733nsvjuMwefLk0KsP+4vOzgy8OGtTARCRmDh06FBohF9cXBy6CSs7O5uHH34Y13XJyckhOTk54seO9SMWuhPruf7OqACISFTU1dURDAZDnf7evW3r38eMGcPixYtxXZeFCxdy1VVXRTUOLx6x0FeoAIj4TLRGwy0tLZSVlYWmdbZs2UJzczODBw9m/vz5fP3rX8d1XSZNmhTTaR3dLNc5FQARH4n0aPjAgQOhDr+kpIQTJ05gjGH69OksX74cx3GYO3cugwYNimgePaGb5TrXbQEwxnyTtkc1nIhBPCISRb0dDZ88eZINGzaEOv19+/YBkJmZye23347jOOTn5zNy5MiIxt2bsxY/3iwXrnDOAEYD24wx24FfAuus10+QE5HL0tPRcHNzM9u2bQvN42/dupWWlhaGDBlCXl4eDzzwAK7rMnHixKhN6/T2rCXaL47py7otANbax4wx3wVc4B7gX40xLwDPWGuroh2gSH8W69Up4YyGq6qqQh1+SUkJJ0+exBhDIBBgxYoVuK7LnDlzGDhwYNTivFBvz1riadllvAnrGoC11hpjPgQ+pO2e8+HAb9rf4PXtaAYoEmnxsiTQi9UpHY2GB7Y08LkBh/n611+gqKiI/fv3A23PyP/rv/5rHMdhwYIFjBgxIioxdaerh6rlrC5h+aKJpHXTRk+XXcbL70i0hXMN4EHgLuAo8AtgubW2yRiTALwLqABInxFPSwK9WJ1yW3YGzU1N/OCXv+X9yjdoeX8HZ2r28pPWVoYOHUpeXh7/+I//iOu6ZGVlxcVNWJ2dtcAnn9+TcyP3xq54+h2JtnDOAEYCS6y171240Vrbaoy5NTphiURHPC0JjNXqFGst+/bto7CwkMLCQjZs2EBdXR0JCQnMmjUL5+5VuK7L7NmzSUqKvzeJ5d0wiue3HurwZSnQ9vl9dLIpYseLp9+RaAvnGsD/6uJ7uyMbjkh0xdOSwGiuTjl+/DglJSWhTv+999rGb+PHj+dLX/oSrutyZsRE/n3LR/y/2gbWbzzH8iEfR6SDi+T0ydqKGl4sr+m08z+vo1cpXq54+h2JNt0HIL4ST0sCI7k6pampidLS0tDF223btmGt5YorrmDhwoU8+uijuK7LhAkTgLaO9fEoTHNEevqko9F4RwYmRu5JoPH0OxJtKgDiK/G0JLA3q1Ostezdu/eiRyY3NjaSmJjI7Nmz+d73vofjOMyaNYsBAz79zzxa0xyRbjecUXdKUiLpwyK3IimefkeiTQVAfCXelgT2ZHXK0aNHWb9+fWiU//777wOQlZWF67rcc8895OXlMWzYsG7bitY0R6Tb7eqtWq3Whj6/tJPvXlb7HYm335FoUgEQ34mXJzF25+zZs2zevDk0yt++fTvWWtLS0li4cCGrVq3CcRyuu+46gsEgubm5YbcdrWmOSLfb2Wj8ySVTLvoMg8HIFQDoO78jvaUCIBInrLXs3r2bwsJCnnvxd1S8sZnWpkZMQiI33DSDH/zgB7iuy4wZMzqc1umJaE1zRLpdP43GveBpATDG/BK4FfjYWjvZy1hEvHDkyBGKi4tDq3UOHz4MwMArMxgyZSHJ46aTfM0UbOpQptw6hdkR6vii1bFGo12/jMa94PUZwLPAvwL/4XEcIjHR2NjIpk2bQvP4FRUVAAwfPpz8/Hxc1+Xf3x3CMXPFRT8XjXXo0epY1WH3HZ4WAGttqTFmnJcxSP8TT7fxW2vZuXNnaB6/tLSUhoYGBgwYQE5ODj/60Y9wXZfp06eTmNh2N+vjK37fYVv9cR26eMt4/WDP9gLwamdTQMaYZcAygPT09BkFBQVdtldfX09qamqkw+wT/Jw7tOXfnDiImhMNtF7we51gDBnDU0hLic1drsePH6e8vJyysjLKy8s5duwYANdccw2BQIBAIMDUqVMZPHhwhz+/98O6Dm9sGpiYwMTRQzv8GX32/s0/nNzz8vLKrbWBS7d7PQXULWvtU8BTAIFAwHa30qGnqyH6Ez/k3tXoPhgMsmprKzW1HT0X5hwZaYlRORtoaGhg48aNoXn8t956C4ARI0aEpnUcx2Hs2LFhtVd7yc1U8MnKl9xOYvfDZ98VP+ffm9zjvgCInBfOXaZdTZNE6m5Xay1vvfVWaB6/tLSUs2fPkpSURE5ODk888QSu65KdnU1CQs/vUNXKF4kVFQDpM8K5y7SrJ0d2tH+4Pvjgg4vuuv3oo48AuPHGG7nvvvtwXZd58+YxZMiQHmbVMV1IlVjwehnofwG5wEhjTDXwPWvtM17GJPErnLtMO1qHHm478MkUU/WREww5/i7XNx9g/47N7Ny5E4BRo0bhOE7oT0ZG5510PF2MFumI16uA/sbL40vfEs5dphdOn3R2JtDRXamtra38y6+LWP2LX1NXVU5j9S5oaWZPYhJTZszmxz/+MY7jMHXq1LCmdfz0THnpuzQFJH1GuHeZnp8+ubQTvnT/mpqai6Z1jhw5AkDSqHFcMf0vSR43jUFjP0PqqOF8+9sLehSrn54pL32XCoD0GT29OHrp/umDwRn+McFfrWHVV4rYtWsXAOnp6SxatIg/nriKQeOmMSD1yovauZz19+FMV2mKSLymAiB9Sk8ujra2tpLZ8gG32jcofKOQTZs28UZTE8nJycybN4977rkH13WZMmUKxhhyVpdE7EFm3U1XaYpI4oEKgPQrhw4doqioiKKiIoqLi0M3YU2dOpWHHnoI13XJyckhJeXTnXokH2TWXVuaIpJ4oAIgfVpdXR2vvfYahYWFrF27NvSM/Kuvvppbb70Vx3HIz88nPT2927Yiuf6+u7b89NpBiV8qAN3QPG18aWlpoby8PHQT1ubNm2lubiYlJYUpU6bw8MMP4zgOn/nMZzDG9Lj9SK6/76otP712UOKXCkAXNE8bHw4ePBharbN+/XpOnDgBwPTp0/nWt76F67rMnTuXrVu3xu3jAC4dSOTdMIoXy2t88dpBiV8qAF3QPK03Tp06RTAYDD1b59132972lJGRwW233YbruixcuJBRo0Z5HGl4OhpIvFhewx0zMtiw54jOLsUzKgBd0DxtbDQ3N1NWVhaa1tmyZQstLS0MHjyYvLw87r//flzX5YYbbrisaR2vdTaQ2LDnCJtW9Oz+ApFIUgHoguZpo2f//v0XTeucPHkSYwwzZszg0UcfxXEcPvvZzzJo0CCvQ+01DSQkXqkAdCFa7031o9raWjZs2BDq9KuqqgAYO3YsX/ziF0PTOiNGjPA40sjTQELilQpAF/RY3svX3NzMm2++GZrHf/PNN2lpaSE1NZW8vDweeughHMfhz/7sz/rktE5PaCAh8UoFoBt6LG94rLVUVVWFRvglJSWcOnWKhIQEZs6cycqVK3Fdlzlz5pCUFJs3c8ULDSQkXqkAyGU7ceIEJSUloYu3Bw4cAGDcuHEsXboU13VZsGABw4cP9zhS72kgIfFIBUDC1tTUxNatW0Oj/G3bttHa2soVV1xBXl4ejzzyCI7jkJWV1e+ndUT6AxUA6ZS1lnfffTc0jx8MBqmrqyMhIYHZs2fz3e9+F8dxmDVrlu+mdUT6AxUAucixY8dYv359aJR/6NAhACZMmMCXv/xlXNclLy+PtLQ0bwMVkV5TAfC5c+fOsWXLltA8fllZGdZahg0bxsKFC1m5ciWO4zBhwgSvQxWRCFMB8BlrLXv27Al1+MFgkNOnT5OYmMicOXP4/ve/j+M4zJw5kwED9Osh0p95/VL4zwM/AxKBX1hrV3sZT3919OhRiouLQ51+dXU1ANdffz133303juOQm5vLsGHDPI5URGLJswJgjEkE/g1wgGpgmzHmFWvtLq9i6i/Onj3Lpk2bQvP4FRUVWGsZPnw4CxcuxHVdHMdh3LhxXocqIh7y8gxgFrDPWrsfwBhTAHwBUAHoIWstu3bt4r//+7/58Y9/TGlpKWfOnGHAgAHMnTuXH/7wh7iuy4wZM0hMTPQ6XBGJE8Za682Bjfki8Hlr7b3tX/8tMNta+81L9lsGLANIT0+fUVBQ0GW79fX1pKamhr6ubWjio5ONnGtpZWBiAunDkklL6ftLFk+cOEF5eTllZWWUl5dz9OhRoO3ZOoFAgEAgwLRp0xg8eLDHkcbOpZ+9n/g5d/B3/uHknpeXV26tDVy63cszgI7uFPpUNbLWPgU8BRAIBGx3L/wIBoOhl4Ksrahh5fpKGpoSgAQAUpJaeHLJjTG5KzOSbxNrbGxk48aNoXn8HTt2AHDllVeSn5+P67qkpqZy5513RjCDvuXCz95v/Jw7+Dv/3uTuZQGoBsZe8HUmcDiSB/DyhS69fZuYtZbKysrQPH5paSmNjY0kJSWRk5PDE088geM4ZGdnh6Z1gsFg1PIRkf7HywKwDbjeGDMeqAGWAl+K5AG8fA57OMXn0jOEe2cMZ8CHOyksLKS4uJgPP/wQgEmTJvG1r30N13WZN2+eb091RSSyPCsA1tpmY8w3gXW0LQP9pbX27Ugew8vnsHdXfNZW1PDor7dRu/8tGg/uoOZgBZuPHARg5MiROI4T+pOZmRn1eEXEfzy9D8Ba+wfgD9Fq38vnsHdUfKxtJa3hA37yk5/wxNO/5uSBSmhpgsQBJGfeyJD5d3PNlNlU/GwZCQkJUY9RRPytX9/q6eVz2M8Xn7rjH9N4cAcNB7dz9uAODp05yZ+ApJHXMjT7z0kZP51BYz9DQlIyAKdAnb+IxES/LgAQ++ewnz59mtLSUl4rLOT07/6Hmqq9ACSlppEzfwFfvfMLOI7DX/3HHr0mUEQ81e8LQLS1trZSUVERWq2zadMmzp07x6BBg7j55pt54Gt/j+u6TJky5aKR/fJFVq8JFBFPqQBchvfff5+ioiKKioooLi4O3YR100038cADD+C6LjfffDMpKZ2P5vWaQBHxmgpAGOrr63nttddCL0bZs2cPAKNHj+aWW27BdV3y8/MZPXp0j9rVawJFxEsqAB1oaWlh+/btobtuN2/eTFNTE8nJycyfP597770X13WZPHmyXn0oIn2WCkC79957LzSPv379eo4fPw7AtGnTePjhh3Fdl5ycHJKTkz2OVEQkMnxbAE6dOkUwGAx1+u+88w4AY8aMYfHixTiOQ35+PldddZXHkYqIRIdvCkBLSwtlZWWhaZ0tW7bQ3NzM4MGDmT9/Pvfddx+u6zJp0iRN64iIL/TrAnDgwIGLpnVqa2sxxjB9+nQeeeQRXNdl7ty5DBo0yOtQRURirl8VgJMnT7Jx40ZeeOEFioqK2LdvHwCZmZksWbIE13VZuHAhI0eO9DhSERHv9ekC0NzczJtvvhka5b/xxhu0tLQwZMgQ8vLyQmvyJ06cqGkdEZFL9LkCUFVVFVqPX1JSwqlTpzDGEAgEWLFiBaNGjeK+++5j4MCBXocqIhLX+lQBqKysJCsrC4Brr72WO++8E8dxWLBgASNGjADaXoqizl9EpHt9qgCkpKTw05/+FNd1ycrK0rSOiEgv9KkCkJWVxf333+91GCIi/YIePC8i4lMqACIiPqUCICLiUyoAIiI+5UkBMMb8lTHmbWNMqzEm4EUMIiJ+59UZwE5gCVDq0fFFRHzPk2Wg1trdgNbxi4h4yFhrvTu4MUHgEWttWRf7LAOWAaSnp88oKCjoss36+npSU1MjGWaf4efcwd/5+zl38Hf+4eSel5dXbq391HR71M4AjDHFQEcvyV1lrf1tuO1Ya58CngIIBAI2Nze3y/2DwSDd7dNf+Tl38Hf+fs4d/J1/b3KPWgGw1uZHq20REek9LQMVEfEpr5aB3m6MqQY+C/zeGLPOizhERPzMq1VALwMve3FsERFpoykgERGfUgEQEfEpFQAREZ9SARAR8SkVABERn1IBEBHxKRUAERGfUgEQEfEpFQAREZ9SARAR8SkVABERn1IBEBHxKRUAERGfUgEQEfEpFQAREZ9SARAR8SkVABERn1IBEBHxKRUAERGf8uSdwMaYNcBfAueAKuAea22tF7FcjrUVNaxZt5fDtQ2MSUth+aKJ3Jad4XVYIiI94tUZQBEw2Vp7E/AOsNKjOHpsbUUNK1+qpKa2AQvU1Daw8qVK1lbUeB2aiEiPeFIArLWF1trm9i+3AplexHE51qzbS0NTy0XbGppaWLNur0cRiYhcHmOt9TYAY34H/Npa+1wn318GLANIT0+fUVBQ0GV79fX1pKamRjzO8yprTnb6vSkZw6J23HBEO/d45+f8/Zw7+Dv/cHLPy8srt9YGLt0etQJgjCkGRnfwrVXW2t+277MKCABLbBiBBAIBW1ZW1uU+wWCQ3NzcngccppzVJdTUNnxqe0ZaCptWLIjaccMR7dzjnZ/z93Pu4O/8w8ndGNNhAYjaRWBrbX43Ad0F3AosDKfzjxfLF01k5UuVF00DpSQlsnzRRA+jEhHpOa9WAX0eeBSYb60940UMl+v8ah+tAhKRvs6TAgD8KzAIKDLGAGy11n7do1h67LbsDHX4ItLneVIArLVZXhxXREQ+oTuBRUR8SgVARMSnVABERHxKBUBExKc8vxO4J4wxR4D3utltJHA0BuHEIz/nDv7O38+5g7/zDyf3a621oy7d2KcKQDiMMWUd3fHmB37OHfydv59zB3/n35vcNQUkIuJTKgAiIj7VHwvAU14H4CE/5w7+zt/PuYO/87/s3PvdNQAREQlPfzwDEBGRMKgAiIj4VL8pAMaYzxtj9hpj9hljVngdTywZY8YaYzYYY3YbY942xvyD1zHFmjEm0RhTYYx51etYYs0Yk2aM+Y0xZk/778BnvY4pVowxD7f/zu80xvyXMSbZ65iiyRjzS2PMx8aYnRdsu9IYU2SMebf9v8PDba9fFABjTCLwb8AtwI3A3xhjbvQ2qphqBr5lrZ0EzAHu91n+AP8A7PY6CI/8DPijtfYGYCo++f9gjMkAHgQC1trJQCKw1Nuoou5Z4POXbFsBrLfWXg+sb/86LP2iAACzgH3W2v3W2nNAAfAFj2OKGWvtB9ba7e1/r6OtA/DNCwuMMZnAXwC/8DqWWDPGXAHMA54BsNaes9bWehpUbA0AUowxA4DBwGGP44kqa20pcPySzV8AftX+918Bt4XbXn8pABnA+xd8XY2POsALGWPGAdnAGx6HEkv/AnwbaPU4Di9cBxwB/m/7FNgvjDFDvA4qFqy1NcA/A4eAD4CT1tpCb6PyRLq19gNoGwwCV4X7g/2lAJgOtvlufasxJhV4EXjIWnvK63hiwRhzK/Cxtbbc61g8MgCYDvy7tTYbOE0PpgD6sva57i8A44ExwBBjzFe8japv6S8FoBoYe8HXmfTzU8FLGWOSaOv8n7fWvuR1PDGUAyw2xhykbepvgTHmOW9DiqlqoNpae/6M7ze0FQQ/yAcOWGuPWGubgJeAuR7H5IWPjDFXA7T/9+Nwf7C/FIBtwPXGmPHGmIG0XQh6xeOYYsa0vVj5GWC3tfanXscTS9baldbaTGvtONo+9xJrrW9GgdbaD4H3jTET2zctBHZ5GFIsHQLmGGMGt/8bWIhPLoBf4hXgrva/3wX8Ntwf9Oql8BFlrW02xnwTWEfbSoBfWmvf9jisWMoB/haoNMbsaN/2HWvtH7wLSWLoAeD59sHPfuAej+OJCWvtG8aY3wDbaVsJV0E/fySEMea/gFxgpDGmGvgesBp4wRjz97QVxb8Kuz09CkJExJ/6yxSQiIj0kAqAiIhPqQCIiPiUCoCIiE+pAIiI+JQKgIiIT6kAiIj4lAqASC8YY2YaY94yxiQbY4a0P5t+stdxiYRDN4KJ9JIx5kdAMpBC23N5nvQ4JJGwqACI9FL7Ixi2AY3AXGtti8chiYRFU0AivXclkAoMpe1MQKRP0BmASC8ZY16h7VHU44GrrbXf9DgkkbD0i6eBinjFGPN3QLO19j/b30292RizwFpb4nVsIt3RGYCIiE/pGoCIiE+pAIiI+JQKgIiIT6kAiIj4lAqAiIhPqQCIiPiUCoCIiE/9f2qY6PYpmksYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the fit\n",
    "plotdata()\n",
    "plotline(w_quad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.535009239560774e-16"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(y - lin.predict(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.06276304347967676"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.median(y - lin.predict(X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ℓ_1 regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhmUlEQVR4nO3df3iV5Z3n8feXEEggCCo0QFBBQEBRiElOOrJ1E6uDs/6oYnsVO9P6oy1uW2ptZ5iRsd25Zrut7LLTa3pd7e6srW47U9tstyJVpyPjVDN06kp+EAQk4VdFSQBFJJRAICG594+THJNwkpyEnNznnPvzui4uyHOePM/3NvH+Ps/93M/9NeccIiISnjG+AxARET+UAEREAqUEICISKCUAEZFAKQGIiARqrO8AhmLq1Klu9uzZA+5z6tQpJk6cODoBpZiQ2w5htz/ktkPY7U+k7bW1te8556b13Z5WCWD27NnU1NQMuE9lZSVlZWWjE1CKCbntEHb7Q247hN3+RNpuZm/F264hIBGRQCkBiIgESglARCRQafUMIJ729nYaGxs5c+YMAJMnT6a+vt5zVCMrJyeHWbNmkZ2d7TsUEckgaZ8AGhsbmTRpErNnz8bMOHnyJJMmTfId1ohxznHs2DEaGxuZM2eO73BEJIN4TQBmdgA4CXQA55xzxUM9xpkzZ2KdfyYyMy699FKOHj3qOxQR8WBjXRPrN+3mUHMrM6fksmb5Au4qLBiRY6fCHUC5c+69CzlApnb+3TK9fSIS38a6JtZu2EFrewcATc2trN2wA2BEkoAeAouIpKj1m3bHOv9ure0drN+0e0SObz7rAZjZm8BxwAH/yzn3RJx9VgGrAPLz84sqKip6fT558mTmzZsX+7qjo4OsrKxkhj2gb3/72+Tl5fHwww/H/fyFF15g3rx5LFy4cEjH3bdvHydOnBhwn5aWFvLy8oZ03EwScvtDbjtkbvt3NPX///y1BZOBxNpeXl5eG2+I3fcQ0DLn3CEz+xDwkpk1OOc299yhKyk8AVBcXOz6vvFWX1/f66Gv74fA48ePZ/z48f3GsGnTJrKzsykpKRnScXNycigsLBxwn5DfhoSw2x9y2yFz2//Yupdpam49b3vBlFy+/MdlwIW13esQkHPuUNff7wLPApFkn3NjXRPL1r3MnEf/kWXrXmZjXdMFH/Nb3/oWCxYs4Oabb2b37uit2Q9+8ANKSkpYsmQJ99xzD6dPn+bVV1/lueeeY82aNSxdupT9+/fH3U9EBGDN8gXkZvce0cjNzmLN8gUjcnxvCcDMJprZpO5/A38I7EzmObsfqDQ1t+L44IHKhSSB2tpaKioqqKurY8OGDVRXVwOwYsUKqquref3111m0aBFPPvkkN9xwA3feeSfr169n27ZtzJ07N+5+IiIQfdD7+IprKZiSixG98n98xbUZMQsoH3i2a4bLWOCnzrkXk3nCgR6oDPc/6G9+8xvuvvtuJkyYAMCdd94JwM6dO/n6179Oc3MzLS0tLF++PO73J7qfiITprsKCEevw+/KWAJxzvwOWjOY5D8UZSxtoe6LiTdO8//772bhxI0uWLOFHP/oRlZWVcb830f1EREZaUNNAZ07JHdL2RNx44408++yztLa2cvLkSZ5//nkg+jB6xowZtLe38/TTT8f2nzRpEidPnox93d9+IiLJFlQCSMYDleuvv55PfvKTLF26lHvuuYePfOQjAHzzm9+ktLSUW265pdeUz5UrV7J+/XoKCwvZv39/v/uJiCSb72mgo6p7HG2kX6t+7LHHeOyxx87b/oUvfOG8bcuWLWPXrl299om3n4hIsgWVACC5D1RERNJJUENAIiLygYxIAD6XsxgNmd4+EfEj7RNATk4Ox44dy9hOsrseQE5Oju9QRCTDpP0zgFmzZtHY2BhbL//MmTMZ11l2VwQTERlJaZ8AsrOze1XKqqysHHTRNBERyYAhIBERGR4lABGRQCkBiIgESglARCRQSgAiIoFSAhARCZQSgIhIoJQAREQCpQQgIhIoJQARkUB5TwBmlmVmdWb2gu9YRERC4j0BAF8B6n0HISISGq8JwMxmAbcBP/QZh4hIiMznOvpm9gvgcWAS8GfOudvj7LMKWAWQn59fVFFRMeAxW1payMvLS0K0qS/ktkPY7Q+57RB2+xNpe3l5ea1zrrjvdm/LQZvZ7cC7zrlaMyvrbz/n3BPAEwDFxcWurKzfXYHoctCD7ZOpQm47hN3+kNsOYbf/QtrucwhoGXCnmR0AKoCbzOwnHuMREQmKtwTgnFvrnJvlnJsNrAReds79ia94RERCkwqzgERExIOUKAnpnKsEKj2HISISFN0BiIgESglARCRQSgAiIoFSAhARCZQSgIhIoJQAREQCpQQgIhIoJQARkUApAYiIBEoJQEQkUEoAIiKBUgIQEQmUEoCISKCUAEREAqUEICISKCUAEZFAKQGIiARKCUBEJFBKACIigVICEBEJlLcEYGY5ZlZlZq+b2Rtm9te+YhERCdFYj+c+C9zknGsxs2zg38zsn5xzr3mMSUQkGN4SgHPOAS1dX2Z3/XG+4hERCY1F+2FPJzfLAmqBecD3nXN/EWefVcAqgPz8/KKKiooBj9nS0kJeXl4Sok19Ibcdwm5/yG2HsNufSNvLy8trnXPF533gnPP+B5gCvAIsHmi/oqIiN5hXXnll0H0yVchtdy7s9ofcdufCbn8ibQdqXJw+NSVmATnnmoFK4Fa/kYiIhMPnLKBpZjal69+5wM1Ag694RERC43MW0Azgx13PAcYAP3fOveAxHhGRoPicBbQdKPR1fhGR0KXEMwARERl9SgAiIoFSAhARCZQSgIhIoJQAREQCpQQgIhIoJQARkUApAYiIBEoJQEQkUEoAIiKBUgIQEQmUEoCISKCUAEREAuVzOWgRkbSzsa6J9Zt2c6i5lZlTclmzfAF3FRb4DmtYdAcgImlrY10Ty9a9zI6mEyxb9zIb65qSfr61G3bQ1NyKA5qaW1m7YUfSz5ssSgAikpZ6dsYwOp3x+k27aW3v6LWttb2D9Zt2J+2cyaQEICJpyUdnfKgr2SS6PdUpAYhIWvLRGc+ckjuk7b61tLSwefPmfj/XQ2ARSUszp+TGhn/6bk+WNcsXsHbDjl53HrnZWaxZviBp50zUuXPn2LlzJ1VVVbE/b7zxBp2dnf1+j7cEYGaXAX8PTAc6gSecc9/1FY+IpBcfnXH3bB/fs4Ccc7z11lts2bKFZ555hm984xvU1tbS2hpNiJdccgmRSIS7776bSCTC7bffHvc4Pu8AzgF/6pzbamaTgFoze8k5t8tjTCKSJnp2xnCSglHqjO8qLEjaOfqbYvr+++/3urKvqqri6NGjAGRnZ1NcXMyqVauIRCKUlpZy5ZVXYmaDns9bAnDOHQYOd/37pJnVAwWAEoCIJKS7M66srOTLf1zmO5wL0j2r6XRrK23v7Ke+Zg/3/WQvE39/gMMHDwBgZixatIjbbrst1tkfO3aMW265ZVjnNOfcCDZheMxsNrAZWOyc+32fz1YBqwDy8/OLKioqBjxWS0sLeXl5SYo0tYXcdgi3/c2t7Zw7c5rDp2Fc1hjyJ+cwJTfbd1ijKlV/9s2t7bxz4gxtHZ1xfzadnZ0cPHiQ+vp6/l/tdt7cv4emtw/Q2REd1pp88aXMmTefkqXXsmjRIq666iomTpzY6xyJtL28vLzWOVfcd7v3BGBmecC/At9yzm0YaN/i4mJXU1Mz4PEqKyspKysbuQDTSMhthzDb333V+MWFZ/mbHdEb+tzsLB5fcW3avp06HKn4s+/+2fR8RjH2TDMfv6wVe28/VVVVVFdXc/LkSQBsXC7jZ8xn3IyrGD9jAeNmzGfspKkY8Oa62/o9TyJtN7O4CcDrLCAzywaeAZ4erPMXkfMNNBc+pASQih7/5Vbe37eDtsN7OHt4D22H9tDRcoz/CmSNHUtO/pWMmfcR5s29hodX/gd+tvsch35/9rzjJHNWk89ZQAY8CdQ7577jKw6RdJZpLyalq/b2dnbs2MGWLVtiD2l37aoHoiMsYy+ewfjLFzO+6+r+olnzOeuyot8L/I9tZ7inqIBnaptGdVaTzzuAZcCngR1mtq1r2186537lLySR9OJjLnzonHPs37+/14ycrVu3cvZs9Op92rRpRCIRTs4s4ezFVzJuxnyyci+KfX+WGWf7DL23tnfwSsNRHl9x7ahOMfU5C+jfgMHnKYlIv7rnwkdnVUelyotJmeLo0aPnTcF8//33AcjNzaWoqIgvfelLlJaWUlJSwuzZszGzuM8AcrOzzhuy63aouTWpU0zj0ZvAImmsu7N4Z/dWDNJ+eWLfTp8+TV1dXa+hnDfffBOAMWPGcM0113D33XdTWlpKJBLhmmuuYezY+N1ofy+Nrd+0O2Xu2pQARNLcXYUFVJ7Yy5vrynyHklY6Ojqor6/v1dnv2LGDjq4pmJdffjmRSIQvfvGLRCIRrr/++iFPNe3vij5VlpNQAhCRjOeco7GxsVdnX1NTw6lTpwCYMmUKJSUlrF27lkgkQklJCdOnT09KLKmynAQoAYhIBmpubqa6ujrW2f/m1dc4/t67AFhWNvMWXcODDz5IJBIhEokwb948xowZvcWRR3usvz9KACKS1s6ePUt9fX2vlTB37/6gJkDB7Ll0zLiWSwqjL1mNmzaHrNwcbgrsZbl4lABEApPONW07OzvZu3dvrxk527Zto62tDYDp06dTWlrKZz7zGSKRCMXFxdz2d1vPe+iql+WilABEAtJ3amJ3GUUgJTvDI0eO9Orsq6uraW5uBiAvL4/i4mIeeeQRJkyYwIMPPsisWbPOWwVTL8v1b9AEYGariS7VcHwU4hGRJErlpSNaWlqora3t1eG//fbbAIzJyiI3fw42+8PMnbuYL3/yVlavKCMrK/o2bWVlJZdddlnc4+pluf4lcgcwHag2s63AU8Am53sFOREZllS5Gh6setWVV17JDTfcwCOPPMKZKXP40W44S3QVzXPA321v54p5RxJKWqlcxcu3QROAc+7rZvYN4A+BB4DvmdnPgSedc/uTHaBIJhvt8XgfV8M9q1d1d/Y9q1ddeumlRCIRVqxYEZuVM3Xq1Nj3L1v3MmcZ/hh+Kk27TDUJPQNwzjkzOwIcIZqALwZ+0VXB68+TGaDISEuVh6A+xuNH42p4oOpVOTk5XH/99Tz00EOxgiZz5swZsHpVf3cnTc2tLFv3MmuWL2DKIDENddplqvyOJFsizwAeBu4D3gN+CKxxzrWb2RhgL6AEIGkjlR6C+hiPH+mr4TNnzrBt27ZeV/f79u0DotWrrr76am6//fZYZ7948WKys4dWrKa/uxb44Of3+A1Zw4o/nlT6HUm2RO4ApgIrnHNv9dzonOs0s/iVhkVSVCo9BPU1Hj/cl5A6OztpaGjodWX/+uuvc+5cdCG6WbNmUVJSwuc+9zlKS0spKipi0qRJFxxv+cJpPP3a2/T34LG1vYN3TrRf8Hm6pdLvSLIl8gzgPw3wWf3IhiOSXKnyEBT8z04ZbJijqanpvCmY3dWrLrroIkpKSlizZk1sFcyZM2fGjvvYpt0cenHzBd9hbKxr4pnapn47/25tHZ3DOn48qfQ7kmx6D0CC4rvT7cnn7JS+wxwH3znGw3/zDzxzyQlONe6mqqqKpqYmALKzs1myZAmf/vSnY6tgXnXVVXGXThjp4ZN4V+PxjMsauWUcUul3JNmUACQoqTQl0NfslLa2Nv7qqed5d88Ozh7aQ9vhPbQfOwg4fgLMnz+fsrKyWGe/ZMkScnJyEjr2SA+fJHLVnZudRf7kcUM+dn9S6Xck2ZQAJCipNiUw2YuCDVa9asyEyYyfcRUTFn2E8TMXMH76fPZ8d+WwzzfSwyf9XY1nmdHpXOznN+XE3mEdP55U+x1JJiUACU6qrMSYDD2rV23ZsoXq6upY9aoJEyZQVFTE6tWref7wRE5dNJusi6b1moJZcIHDHCM9fNLf1fjjfRZyq6wcuQQAmf070pMSgEgKSmQe+unTp9m6dStVVVU8//zzPPDAAxw4cACIVq9avHgxK1asiA3lXH311bHqVf+un3KFFzrMMdLDJyFdjfvgNQGY2VPA7cC7zrnFPmMRSRXxHqQ++ottvLW3nryTb8Wu7nfu3BmrXpWfn8+NN97I6tWrY9WrJk6c2O85ktWxJuO4oVyN++D7DuBHwPeAv/cch0jK+G8vNnDy2OHYA9qzh/fQdmQfj7SfAaLVqyKRCHfccUdsCmZ9fT1lZWVDOk+yOlZ12OnDawJwzm02s9k+Y5DMk26v8R8/fpyamprY27RbXv4Nnaeaox9mjWXch+aSd90tjJ9xFa9+5yHmzZt33tIJ9fV6JUeGznwv7NmVAF7obwjIzFYBqwDy8/OLKioqBjxeS0vLkAs3Z4qQ2w7R9p/LGk/T8VY6e/xejzGj4OJcpuQObQmCZGhra2P//v3U19fT0NBAQ0MDBw8ejH1++eWXM/OKucy6cj6z517FzMuvYOzYaNzjssawYHr8N2v1sw+3/Ym0vby8vNY5V9x3e8ongJ6Ki4tdTU3NgPtUVlYO+VY4U4TQ9oGu7isrK3nstc5+140pGOW7ge7qVT3Xydm2bRvt7dFlC7qrV3Wvk1NcXMzkyZPPewYA8We+9BTCz34gIbc/kbabWdwE4PsZgEjCEnnLdKD55sle1OvIkSO9Ovvq6mpOnDgBfFC96qtf/Wqs0y8oKIi7CqZmvshoUQKQtJHIW6YDrRwZb//hamlpoaamptcLVt1DOVlZWVx33XWsXLky1tkvXLgwVr0qEXqQKqPB9zTQnwFlwFQzawT+yjn3pM+YJHUl8pZpvHnoiR4H4g8x3bb4Q+dVr9q1a1ev6lXLli2LFTMpLCxkwoQJsWP95384wMwp7+gqXlKO71lA9/o8v6SXRN4y7Tl80t+dQH9vpW6sa+LRZ7Zz8r1DtB3azY7De/jU9/fS8d7vaDsTnYLZXb3q4x//OJFIhJKSkl7Vq3oeK5Q15SV9aQhI0kaib5l2D5/09zC15/7Hjh2LXdX/7c/+id+/XU9n6+8BsLHjGJc/l0uLbuM7q6Md/mDVq7qFtKa8pC8lAEkbQ3042nf//IljuOuysxz41//Lp9ZHO/39+6Nlrc2MsZdeRu68SHRRtBlXkT31CixrLAasXHnbkGJNZLgq3d5XkMyjBCBpJdGHox0dHTQ0NHB8WxXXvVXFmaoqardvZ0uP6lWRSITPf/7zsepVt36/esQWMhtsuEpDRJIKlAAkIzQ1NbF582ZefPFFqqqqqKmpSah6VU8juZDZYMfSEJGkAiUASTsnTpw4bwrmoUOHgKFVr+prJOffD3askMoOSupSAhiExmn9amtrY/v27b06+4aGBrrfYJ8/fz7l5eWUlpaSlZXFgw8+mHD1qnhGcv79QMcKqeygpC4lgAFonHZ0OefYt29fr86+rq4uVr1q2rRplJaWcu+998aWTrjkkkti319ZWXlBnX8y9b2QKF84jWdqm4IoOyipSwlgABqnTa533323V2dfVVXF8ePHgd7Vq7qHci6//PKEpmCmmngXEs/UNnFPUQGvNBzV3aV4owQwAI3TjpzTp09TW1vbq7PvW73qnnvuiVu9Kt31dyHxSsNRfvvoTZ6iElECGJDGaYeno6ODXbt29VoYrWf1qiuuuIJIJJJw9ap0pwsJSVVKAAMY6fqmmcg5x8GDB3t19rW1tZw6dQr4oHrVnXfeGVs6IT8/33PUo0sXEpKqlAAGoGV5z3f8+HGqq6t7DeW88847AIwbN46lS5fy4IMPxoZy4lWvCo0uJCRVKQEMIuRlec+ePcu2bdt6dfZ79uyJfb5w4UKWL18e6+yvu+46xo0b5zHi1KQLCUlVSgACRKtX7dmzp1dnH6961X333derepUkJuQLCUldSgCBOnz4cK/OfrjVq0QkfSkBBCDZ1atEJD0pAWSY9vb2hKpXdXf2S5cuZcKECZ6jFhEflADSmHOOAwcOxKZgvvTSS+zfv5/W1uiUw0SrV4lImJQA0kjP6lXdf9577z0AcnJymDt3Lg899FDs6j7R6lUiEibfReFvBb4LZAE/dM6t8xlPKmltbaWurq5XZ9+zetXVV1/NHXfcEevsFy9ezG9/+1vKysr8Bi4iacNbAjCzLOD7wC1AI1BtZs8553b5ismX7upVPTv77du3c65P9apVq1YRiUQoKipi0qRJnqMWkXTn8w4gAuxzzv0OwMwqgI8BGZ8Ampqaei2dMJzqVSIiF8q6C2uM+onNPg7c6pz7XNfXnwZKnXOr++y3ClgFkJ+fX1RRUTHgcVtaWsjLy4t93dzazjsnztDW0cm4rDHkT85hSm72CLdm4Hh2795NQ0ND7E/3uP3YsWOZO3cuCxcuZNGiRSxatIhZs2YlVL2qv3P1bHtoQm5/yG2HsNufSNvLy8trnXPFfbf7vAOI93TyvGzknHsCeAKguLjYDTbGXVlZGRsH31jXxNpf76C1fQwQ7VRzszt4fMXVSXkrs2/1ql9vfpXGA/ugK8nOuHxOr6UTlixZMqIFTHq2PUQhtz/ktkPY7b+QtvtMAI3AZT2+ngUcGskTJLOgy2DVqy66+FI6ps5l8rJPMX7mAsZNn0/eRVP4+IprtSSAiKQEnwmgGphvZnOAJmAl8KmRPMFIrsM+1OpV9/50H4dOnOl1jL7JR/WGRcQnbwnAOXfOzFYDm4hOA33KOffGSJ5juOuwj0T1qsMndsY9dnfyUb1hEfHN63sAzrlfAb9K1vETWYc9kepVpaWlQ65eNVjyUb1hEfEto98E7rsO+4zJOTywZBLt+15lzU/Pr1518cUXj1j1qsGSj8oEiohvGZ0Ajh8/zoSjb/BH56qo2hXt8Fd1Va8aP348hYWFfPaznyUSiYx49arBioCoTKCI+JYxCaC7etWGDRt48skn41avuvXWW2Od/WhUrxqoCIjKBIqIb2mZADo7O9m7d2+vcfue1atmzJhBaWkp999/P5FIJCWrV6lMoIj4llYJoKmpiVtuueW86lUlJSV87WtfIxKJ0NHRwSc+8QnPkSZGZQJFxKe0SgBHjhxh5syZ3HvvvbGhnL7VqyorK/0FKCKSRtIqARQWFlJbW+s7DBGRjDC8Vcc8Ge4iaSIicj71qCIigVICEBEJlBKAiEiglABERAKlBCAiEiglABGRQCkBiIgESglARCRQSgAiIoFSAhARCZQSgIhIoJQAREQC5SUBmNknzOwNM+s0s2IfMYiIhM7XHcBOYAWw2dP5RUSC56UegHOuHhixAuwiIjJ05pzzd3KzSuDPnHM1A+yzClgFkJ+fX1RRUTHgMVtaWsjLyxvJMNNGyG2HsNsfctsh7PYn0vby8vJa59x5w+1JuwMws38Bpsf56DHn3C8TPY5z7gngCYDi4mJXVlY24P6VlZUMtk+mCrntEHb7Q247hN3+C2l70hKAc+7mZB1bREQunKaBiogEytc00LvNrBH4A+AfzWyTjzhERELmaxbQs8CzPs4tIiJRGgISEQmUEoCISKCUAEREAqUEICISKCUAEZFAKQGIiARKCUBEJFBKACIigVICEBEJlBKAiEiglABERAKlBCAiEiglABGRQCkBiIgESglARCRQSgAiIoFSAhARCZQSgIhIoJQAREQC5aUmsJmtB+4A2oD9wAPOuWYfsQzHxrom1m/azaHmVmZOyWXN8gXcVVjgOywRkSHxdQfwErDYOXcdsAdY6ymOIdtY18TaDTtoam7FAU3NrazdsIONdU2+QxMRGRIvCcA598/OuXNdX74GzPIRx3Cs37Sb1vaOXtta2ztYv2m3p4hERIbHnHN+AzB7Hvg/zrmf9PP5KmAVQH5+flFFRcWAx2tpaSEvL2/E4+y2o+lEv59dWzA5aedNRLLbnupCbn/IbYew259I28vLy2udc8V9tyctAZjZvwDT43z0mHPul137PAYUAytcAoEUFxe7mpqaAfeprKykrKxs6AEnaNm6l2lqbj1ve8GUXH776E1JO28ikt32VBdy+0NuO4Td/kTabmZxE0DSHgI7524eJKD7gNuBjybS+aeKNcsXsHbDjl7DQLnZWaxZvsBjVCIiQ+drFtCtwF8A/945d9pHDMPVPdtHs4BEJN15SQDA94DxwEtmBvCac+4/eoplyO4qLFCHLyJpz0sCcM7N83FeERH5gN4EFhEJlBKAiEiglABERAKlBCAiEijvbwIPhZkdBd4aZLepwHujEE4qCrntEHb7Q247hN3+RNp+hXNuWt+NaZUAEmFmNfHeeAtByG2HsNsfctsh7PZfSNs1BCQiEiglABGRQGViAnjCdwAehdx2CLv9Ibcdwm7/sNuecc8AREQkMZl4ByAiIglQAhARCVTGJAAzu9XMdpvZPjN71Hc8o8nMLjOzV8ys3szeMLOv+I5ptJlZlpnVmdkLvmMZbWY2xcx+YWYNXb8Df+A7ptFiZl/t+p3faWY/M7Mc3zElk5k9ZWbvmtnOHtsuMbOXzGxv198XJ3q8jEgAZpYFfB/4I+Bq4F4zu9pvVKPqHPCnzrlFwIeBLwXWfoCvAPW+g/Dku8CLzrmFwBIC+e9gZgXAw0Cxc24xkAWs9BtV0v0IuLXPtkeBXzvn5gO/7vo6IRmRAIAIsM859zvnXBtQAXzMc0yjxjl32Dm3tevfJ4l2AMEULDCzWcBtwA99xzLazOwi4EbgSQDnXJtzrtlrUKNrLJBrZmOBCcAhz/EklXNuM/B+n80fA37c9e8fA3clerxMSQAFwMEeXzcSUAfYk5nNBgqBLZ5DGU1/C/w50Ok5Dh+uBI4C/7trCOyHZjbRd1CjwTnXBPx34G3gMHDCOffPfqPyIt85dxiiF4PAhxL9xkxJABZnW3DzW80sD3gGeMQ593vf8YwGM7sdeNc5V+s7Fk/GAtcD/9M5VwicYghDAOmsa6z7Y8AcYCYw0cz+xG9U6SVTEkAjcFmPr2eR4beCfZlZNtHO/2nn3Abf8YyiZcCdZnaA6NDfTWb2E78hjapGoNE5133H9wuiCSEENwNvOueOOufagQ3ADZ5j8uEdM5sB0PX3u4l+Y6YkgGpgvpnNMbNxRB8EPec5plFj0cLKTwL1zrnv+I5nNDnn1jrnZjnnZhP9ub/snAvmKtA5dwQ4aGYLujZ9FNjlMaTR9DbwYTOb0PX/wEcJ5AF4H88B93X9+z7gl4l+o6+i8CPKOXfOzFYDm4jOBHjKOfeG57BG0zLg08AOM9vWte0vnXO/8heSjKIvA093Xfz8DnjAczyjwjm3xcx+AWwlOhOujgxfEsLMfgaUAVPNrBH4K2Ad8HMz+yzRpPiJhI+npSBERMKUKUNAIiIyREoAIiKBUgIQEQmUEoCISKCUAEREAqUEICISKCUAEZFAKQGIXAAzKzGz7WaWY2YTu9amX+w7LpFE6EUwkQtkZv8FyAFyia7L87jnkEQSogQgcoG6lmCoBs4ANzjnOjyHJJIQDQGJXLhLgDxgEtE7AZG0oDsAkQtkZs8RXYp6DjDDObfac0giCcmI1UBFfDGzzwDnnHM/7apN/aqZ3eSce9l3bCKD0R2AiEig9AxARCRQSgAiIoFSAhARCZQSgIhIoJQAREQCpQQgIhIoJQARkUD9fzWRQniGtheAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "l1 = QuantileRegressor(quantile=.5, fit_intercept=True, alpha=0).fit(X,y)\n",
    "w_l1 = [l1.coef_[0], l1.intercept_]\n",
    "\n",
    "# plot the fit\n",
    "plotdata()\n",
    "plotline(w_l1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1036974185259609"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(y - l1.predict(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.845146423486767e-12"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.median(y - l1.predict(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Quantile regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjNklEQVR4nO3df3TV9Z3n8ecnIYRAIFHAAAHlN0GM/Aok97LWRG1xZlx/YH/Y2pZ66nK2U9RuWxwZ3JmzO9vRc5ids3O2M3vWVredM06zbUVrbWeoU83alvAj4Yf8iCkgAgmoQA0kEEIS3vtHkisX8uMGcu/33vt5Pc7xSL753u/3/fHi5/35fr6fH87MEBER/2QEHYCIiARDCUBExFNKACIinlICEBHxlBKAiIinhgUdwGCMGzfOpk6d2u85Z8+eZdSoUYkJKMn4XHbwu/w+lx38Ln8sZa+trT1pZuMvP55SCWDq1KnU1NT0e05VVRXl5eWJCSjJ+Fx28Lv8Ppcd/C5/LGV3zh3u7bi6gEREPKUEICLiKSUAERFPpdQ7gN60t7fT0NDA+fPnAcjLy6Ouri7gqIbWiBEjmDx5MllZWUGHIiJpJOUTQENDA6NHj2bq1Kk452hubmb06NFBhzVkzIxTp07R0NDAtGnTgg5HRNJIoAnAOfce0Ax0Ah1mVjLYa5w/fz5S+acj5xxjx47lxIkTQYciIgF4ZUcj6zfWc6yplUn5OaxZPof7FxYOybWT4QmgwsxOXssF0rXy75Hu5ROR3r2yo5G1G3bT2t4JQGNTK2s37AYYkiSgl8AiIklq/cb6SOXfo7W9k/Ub64fk+i7I/QCcc4eAjwAD/reZPdfLOauAVQAFBQWLKysro36fl5fHzJkzIz93dnaSmZkZz7D79dd//dfk5uby+OOP9/r71157jZkzZ1JUVDSo6x44cIDTp0/3e05LSwu5ubmDum468bn8Ppcd0rf8uxv7/n++uDAPiK3sFRUVtb11sQfdBbTMzI45524AXnfOvWNmb116QndSeA6gpKTELp/xVldXF/XSN+iXwNnZ2WRnZ/cZw8aNG8nKymLJkiWDuu6IESNYuHBhv+f4PBsS/C6/z2WH9C3/umffoLGp9Yrjhfk5PPZwOXBtZQ+0C8jMjnX/+0PgZWBpvO/5yo5Glj37BtOe+gXLnn2DV3Y0XvM1v/Od7zBnzhzuuusu6uu7Hs2+973vsWTJEubPn8+DDz7IuXPn2LRpE6+++ipr1qxhwYIFHDx4sNfzREQA1iyfQ05WdI9GTlYma5bPGZLrB5YAnHOjnHOje/4MfArYE8979rxQaWxqxfj4hcq1JIHa2loqKyvZsWMHGzZsYNu2bQCsWLGCbdu2sWvXLubOncvzzz9POBzm3nvvZf369ezcuZMZM2b0ep6ICHS96H1mRTGF+Tk4ulr+z6woTotRQAXAy90jXIYB/2xm/xrPG/b3QuVq/4P+5je/4YEHHmDkyJEA3HvvvQDs2bOHp59+mqamJlpaWli+fHmvn4/1PBHx0/0LC4eswr9cYAnAzN4F5ifynsd66Uvr73isehum+ZWvfIVXXnmF+fPn84Mf/ICqqqpePxvreSIiQ82rYaCT8nMGdTwWn/jEJ3j55ZdpbW2lubmZn//850DXy+iJEyfS3t7Oiy++GDl/9OjRNDc3R37u6zwRkXjzKgHE44XKokWL+NznPseCBQt48MEHue222wD4q7/6K0pLS/nkJz8ZNeTzoYceYv369SxcuJCDBw/2eZ6ISLwFPQw0oXr60YZ6WvW6detYt27dFce/9rWvXXFs2bJl7Nu3L+qc3s4TEYk3rxIAxPeFiohIKvGqC0hERD6WFgkgyOUsEiHdyyciwUj5BDBixAhOnTqVtpVkz34AI0aMCDoUEUkzKf8OYPLkyTQ0NETWyz9//nzaVZY9O4KJiAyllE8AWVlZUTtlVVVVDbhomoiIpEEXkIiIXB0lABERTykBiIh4SglARMRTSgAiIp5SAhAR8ZQSgIiIp5QAREQ8pQQgIuIpJQAREU8FngCcc5nOuR3OudeCjkVExCeBJwDgCaAu6CBERHwTaAJwzk0G/gT4fpBxiIj4yAW5jr5z7qfAM8Bo4Ntmdk8v56wCVgEUFBQsrqys7PeaLS0t5ObmxiHa5Odz2cHv8vtcdvC7/LGUvaKiotbMSi4/Hthy0M65e4APzazWOVfe13lm9hzwHEBJSYmVl/d5KtC1HPRA56Qrn8sOfpff57KD3+W/lrIH2QW0DLjXOfceUAnc4Zz7pwDjERHxSmAJwMzWmtlkM5sKPAS8YWZfDCoeERHfJMMoIBERCUBSbAlpZlVAVcBhiIh4RU8AIiKeUgIQEfGUEoCIiKeUAEREPKUEICLiKSUAERFPKQGIiHhKCUBExFNKACIinlICEBHxlBKAiIinlABERDylBCAi4iklABERTykBiIh4SglARMRTSgAiIp5SAhAR8ZQSgIiIp5QAREQ8FVgCcM6NcM5tdc7tcs7tdc79l6BiERHx0bAA790G3GFmLc65LOC3zrl/MbPNAcYkIuKNwBKAmRnQ0v1jVvc/FlQ8IiK+cV31cEA3dy4TqAVmAn9vZn/WyzmrgFUABQUFiysrK/u9ZktLC7m5uXGINvn5XHbwu/w+lx38Ln8sZa+oqKg1s5IrfmFmgf8D5ANvArf0d97ixYttIG+++eaA56Qrn8tu5nf5fS67md/lj6XsQI31UqcmxSggM2sCqoC7g41ERMQfQY4CGu+cy+/+cw5wF/BOUPGIiPgmyFFAE4Efdr8HyAB+bGavBRiPiIhXghwF9DawMKj7i4j4LineAYiISOIpAYiIeEoJQETEU0oAIiKeUgIQEfFUkMNARUQkDs6cOcPWrVuprq5m06ZNfZ6nBCAiksLMjAMHDrBp06ZIhb9nzx4uXryIc4558+b1+VklABGRFHLu3Dlqamoilf1bb71FU1MTAGPGjKGsrIwVK1YQCoUoLS0lLy8P51yv11ICEBFJUmbG0aNHo1r3O3fupKOjA4DZs2dTVlbGAw88QDgcZu7cuWRmZsZ8fSUAEZEk0dbWxo4dOyKVfXV1NY2NjQCMHDmSpUuXsmbNGsLhMGVlZYwbN46qqirKy8uv6n5KACIiAXn//fcjlf2mTZuora2lra0NgKlTp3L77bcTCoUIh8PceuutDBs2tFW2EoCISAJ0dHTw9ttvR7XuDx06BMDw4cMpKSlh9erVhMNhQqEQEydOjHtMSgAiInFw6tQpNm/eHGndb926lXPnzgEwceJEwuEwq1evJhQKsWjRIrKzsxMeoxKAiMggvLKjkfUb6znW1Mqk/BzWLJ/DvfMnsm/fvqjWfX19PQCZmZksWLCAr371q5HunBtvvLHPkTmJpAQgIimrpzJ+aEoz6559gzXL53D/wsK43m/tht2cbWmm7Vg9exvr+OL36rEP93Ou5QwAY8eOJRwOs3LlSsLhMCUlJYwaNSpuMV0LJQARSUk9lXFreydMgcamVtZu2A0wpEng0olWf/YPP+EPh/bSfuIwYIAja/xNXH/zbfzDn36GUCjErFmzkqJ1HwslABFJSes31ndV/pdobe9k/cb6a0oA586dY9u2bZGunOrqak6ePAmAGz6S7MIiRs4Ok104l+xJs8nIHoUDVq78k2spTiCUAEQkJR1rah3U8d6YGUeOHImaaLVr167IRKs5c+Zwzz33EA6HCYfD/IefHefYmbYrrjMpP+fqChEwJQARSUmT8nNo7KWy768ybmtrY/v27VEV/vHjx4GuiValpaU8+eSTkYlWY8eOjfr8kxfyP+526paTlcma5XOGqFSJFVgCcM5NAf4RmABcBJ4zs78LKh4RSS1rls8ZsDI+fvx4VGVfW1vLhQsXAJg2bRoVFRWR1n1xcfGAE616upYuHwUUzxfP8RTkE0AH8C0z2+6cGw3UOudeN7N9AcYkIini0soYmpk4OotPT+2k4Xcv84X1XZX+e++9B0B2djYlJSU8/vjjkYlWEyZMuOr7xqvC722IaTyTS2AJwMyOA8e7/9zsnKsDCgElABEZ0MmTJxnWuIPbz2zihb/5F37/+9+zuXui1aRJkwiHw5EKf8GCBYFMtBqMqFFNxG9U06WcmcXlwoMKwrmpwFvALWZ25rLfrQJWARQUFCyurKzs91otLS3k5ubGKdLk5nPZwd/yN7W203H+HMfPwfDMDAryRpCfkxV0WEOqs7OTw4cPs3fvXvbu3cu+ffs4evQo0DXRavr06RQXFzNv3jzmzZvHDTfckBRDMZta2/ng9HkudF4c8Lupf7+ZC50Xrzg+PDODORNG93mPWP7eV1RU1JpZyeXHA08Azrlc4P8B3zGzDf2dW1JSYjU1Nf1e71pWxkt1Ppcd/Cx/T6vxT4va+O+7ux7oc7IyeWZFccr2SwOcPn2aLVu2RPrvN2/ezJkzXW3D8ePHR2bUhkIhSkpK2Lp1a9J995e36KH/72baU7+gt9rYAYee7XuIaSx/751zvSaAQEcBOeeygJeAFweq/EXkSvEaC59IZsb+/fujXtbu3bsXMyMjI4NbbrmFL3zhC5EKf8aMGUnRuh9If99Nz+8v7eu/mlFN1yrIUUAOeB6oM7O/DSoOkVQ2FGPhE+3s2bNXTLQ6deoUAPn5+ZSVlfHZz36WUCjE0qVLGTNmTMARX52+voOevv3L+/ofXFzIS7WNCR1iGuQTwDLgS8Bu59zO7mN/bma/DC4kkdQSRKtxMMyMw4cPXzHRqrOzq5IrKirivvvui3TpFBUVkZGREXDUQ6Ov7ybTuV6fDN585wTPrCj2ZhTQb+nq3hKRq9QzFr5rVHWXICcmnT9/nu3bt0etitkz0WrUqFGUlpaydu1aQqEQZWVlXH/99YHEmQh9zVO4vPLvcaypNa5DTHujmcAiKaynsvigfjsOEj4x6dixY1E7Wm3fvj0y0Wr69Onceeedkdb9LbfcMuQ7WiWzviaNrd9YnzRPbf58GyJp6v6FhVSd3s+hZ8vjep/29nZ27doV1bo/fPgw0DXRasmSJTzxxBORl7UFBQVxjScV9NWiT5blJJQARKRXJ0+ejGrdb9u2jdbWrpbr5MmTCYfDfOMb34hMtBo+fHjAEaeGZFpOQglAROjs7GTv3r1Rrfv9+/cDMGzYMBYtWsSqVasirfspU6YEHPHgJHqJhYEkuq+/L0oAIh5qamqKTLTatGkTW7Zsobm5GYAbbriBcDjMo48+SjgcZvHixeTkJMeooqsRxBILqUIJQCTNmRn19fWR1v3GN3/D0Xd/D2bgMpg2ey5f/OIXI6376dOnp8REq1ilw2S5eFECEEkzLS0tV0y0+sMf/gBA7pg8bPws8v7dw2RPKiJ74myyc0fzqRRfOqI/qThZLlEGTADOudV0LdXwUQLiEZFBMDMOHTrE66+/zk9+8hOqq6vZtWsXFy92LSp2880388ADD0Ra9199ufGKHa1SoTV8LX34yT5ZLkixPAFMALY557YDLwAbLegV5EQ8df78eWpra6Nm1n7wwQcA5ObmUlpayrp16yITra677rqozx//4bu9XjeZW8PX2ocfy8YxvhowAZjZ0865/wx8CngE+K5z7sfA82Z2MN4BiqSzgVq2DQ0NUSNztm/fTnt7OwAzZszgU5/6FKFQiKysLB555BEyMzP7vV8qtoavtQ8/mYZdJpuY3gGYmTnn3gfep2vO+XXAT7t38HoyngGKDLVkGRJ4ecu24VQz3/juT/n59U2cPVrHpk2bImvejxgxgiVLlvDNb36TUChEKBTihhtuiFyrqqpqwMofUrM13N+iasuefYM1y+eQP8A1BjvsMln+jsRbLO8AHgdWAieB7wNrzKzdOZcB7AeUACRlJNOQwO+8tJlTdTtpa3yHtsY6Lry/H+u4wAvAlClTIv324XCY+fPnD8lEq1RsDff11AIff3/PhAdOfrFKpr8j8RbLE8A4YIWZHb70oJlddM7dE5+wROIjqCGBnZ2d7NmzJ6rv/uDB7h7UjGEML5hB7oI/IrtwLiMmFXHkH1bGLZZkmYQUq4qi8by4+Uivm6VA1/f3wen2IbufT8NGY3kH8Bf9/K5uaMMRia9EDQn86KOP2Lx5c6Sy37JlCy0tLQAUFBQQDofpmH0n56+fyfCCGWRkfbxfbWGC+uPj1c0xlNd9ZUcjL9U29ln59+htK8Wr5dOwUc0DEK/E4yXoxYsXqa+vj2rd19V1tY0yMjKYP38+X/7ylyNdOtOmTcM51+eWgYnoj49XN8dQX7e31nhvhmcO3R4Cqfii/GopAYhXhuIlaHNzM1u3bo1U9ps3b+ajj7qmyVx33XWEQiEefvjhyI5WfW3YHWR/fLy6OYb6urG0unOyMinIG7qF6FLxRfnVUgIQrwy20jUz3n333ahVMXfv3h010erBBx+MvKydPXv2oHa0Cqo/Pl7dHEN93f521bpoFvn+8k/vv6rr9yYVX5RfLSUA8U5/lW5rays1NTVRY+8//PBDAEaPHk1paSlPP/00oVCI0tLSKyZapYp4dXMM9XX7ao0/c9nSFVVVQ5cAIPVelF8tJQDx2tGjR6+YaNXR0bW94syZM7n77rsjrft58+bFNNZ+KMR7HHq8ujmG+ro+tcaDEGgCcM69ANwDfGhmtwQZi6S/CxcusHPnzqiXtQ0NDUDXRKulS5fyrW99i3A4TFlZWdREq0RKxDj0eFWs8biuL63xIAT9BPAD4LvAPwYch6ShDz74IKp1X1NTw/nz5wG48cYbWbZsGeFwODLRKisrK+CIuyRqHHq8KlZV2Kkj0ARgZm8556YGGYOkh46ODg4cOMC+ffv48S9+zebN1bT94TgAw4ZlUVKymK997WuRoZiFhclbQfk0Dl2C5YJe2LM7AbzWVxeQc24VsAqgoKBgcWVlZb/Xa2lp6XPYXbrzqexnzpxh37597Nu3jz179lBXVxdp3Y/Oy2farKLuf+Zw07SZTC3IIz8nOVr4A6l/v7nXiU3DMzOYM2F0r5/x6bvvjc/lj6XsFRUVtWZWcvnxpE8AlyopKbGampp+z6mqqqK8vHxogksx6Vr2ixcv8s4771BdXc2Pfv46mzZV03riCAAZmZnMv/VWwuEw+fn5/MuFuZx0eb3uaFWYIi8Q+5ogdvnIl0ul63cfK5/LH0vZnXO9JoCg3wGIXKG5uZktW7ZETbRqamoCICNnNNmTisifW072pCLybiziLx5ayv0LC6mqquLFfz1LX5sZpsqiXhr5IomiBCCBMjMOHjwY9bK2Z6KVc46bb76Zz3zmM4RCIf7nngxODRsb1bpvg6iXo/2tHAmps6iXXqRKIgQ9DPRHQDkwzjnXAPylmT0fZEwSXz0TrXpm1VZXV3PixAmga6JVWVkZTz/9NOFwmNLSUvLz8yOf/a9P/aLX1v2lL0d7G4fe3/mXG+qFzNSKl2QW9Cigzwd5f4m/o0ePRlX2O3bsiEy0mjVrFn/8x38cmWh188039zvRKpZZppd2n/T1JNDXrNShHH/v05rykrrUBSRD5sKFC+zYsSNqolVjYyMAOTk5LF26lG9/+9uRiVbjx48f1PVjnWXa030y2NU2h3L8vU9rykvqUgKQq/b+++9TXV0dqexrampoa2sD4KabbuK2226LTLS69dZbr3mi1WBfjg72/KEcfx/LtdRFJEFTApCYdHR0sHv37qjW/aFDhwAYPnw4ixYt4utf/3qkO2fSpElxiWOwL0cHc/5QLmQ20LXURSTJQAlAenXq1KmoHa22bt3K2bNnAZgwYQLhcDhS4S9atIgRI0YEHPG1G8qFzAa6lrqIJBkoAQgXL16krq4uas37+vp6ADIzM1mwYAGPPPJIZBmFm266qdeJVqluKMffD3QtLfcgyUAJYADp2E975swZtm7dGqnsN2/ezOnTpwEYO3YsoVCIlStXEgqFWLJkCaNGjQo44sQZyvH3/V3Lp20HJXkpAfQjHfppeyZaXToUc/fu3ZgZzjnmzZvH5z73uUjrftasWWnZug/a5Q2JiqLxvFTb6MW2g5K8lAD6kYr9tOfOnYtMtOrp0jl58iQAY8aMoaysjBUrVkR2tMrLyws44vTXW0PipdpGHlxcyJvvnEirp0tJLUoA/Uj2flozi5po9atf/YqDBw9GJlrNnj2be+65JzIyZ+7cuQnb0Uo+1ldD4s13TvC7p+4IKCoRJYB+JVs/bVtb2xUTrY4dOwbAyJEjmT17Nk8++SShUIiysjLGjRsXSJwSLdkbEuIvJYB+xGvf1FgdP348aqJVbW1tZKLV1KlTKS8vj/Td33rrrfz2t7/1dkncZJZsDQmRHkoA/UjksrwdHR28/fbbUUMx33vvPQCys7NZvHgxjz32GKFQiFAoxMSJE4c8BomPoBsSIn1RAhhAvJblPXXqVFTrfuvWrZw7dw6ASZMmEQ6HeeyxxwiHwyxcuJDs7Owhj0ESQ+v7S7JSAkiAixcvsm/fvqg173smWg0bNowFCxbw6KOPRl7WTpkyRUMx04zW95dkpAQQB2fOnGHLli2Ryv7SiVbjxo0jHA7zla98hXA4TElJCSNHjgw4YhHxkRLANTIzDhw4EDUyZ8+ePZGJVsXFxXz+85+PtO5nzJih1r2IJAUlgEE6d+4c27Zti1T41dXVkYlWeXl5hEIhPv3pTxMOh1m6dCljxowJOGIRkd4pAfTDzDhy5EhU637nzp10dnaN5igqKuLee++NtO6LiorIyMgIOGoRkdgoAVyira2N7du3R1X4x48fB2DUqFGUlpby1FNPRXa0uv766wOOWETk6gW9KfzdwN8BmcD3zezZRN7/+PHjUZV9bW0tFy5cAGD69OnccccdkYlWxcXFDBumfCki6SOwGs05lwn8PfBJoAHY5px71cz2xeN+7e3tUROtqquroyZaLVmyhCeeeCLSup8wYUI8whARSRpBNmmXAgfM7F0A51wlcB8wJAng5MmTUROttm3bFploVVhYSDgc5oknniAUCrFw4UKGDx8+FLcVEUkZzsyCubFznwbuNrNHu3/+ElBqZqsvO28VsAqgoKBgcWVl5RXX6uzs5PDhw+zdu5edO3eyf/9+jh49CnTtaDV56nSmzpzDrDlzKV08n5k3pueEnJaWFnJzc4MOIzA+l9/nsoPf5Y+l7BUVFbVmVnL58SCfAHobDH9FNjKz54DnAEpKSqy8vJzTp09fMdHqzJkzAOTn53P77bezevVqOsbO5If7M2ljGO8B7wG/3Z/JM8WzEjIrM9G7iVVVVXm9GJzP5fe57OB3+a+l7EEmgAZgyiU/TwaO9feBw4cPU1xczN69ezEzMjIyKC4u5uGHH44MxTxy5AgVFRUALHv2DdqIXoUxURu6pMNuYiKS3oJMANuAWc65aUAj8BDwhf4+8NFHH7FkyRI++9nPRiZajR49Ouqcnq4fCHYd9lh2E0vH/YZFJHUElgDMrMM5txrYSNcw0BfMbG9/n1mwYAG//OUvY75HkOuwD5R89IQgIkELdNqqmf3SzGab2Qwz+85QX3/N8jnkZEVvgZioddj7SjI9x/t7QhARSYS0Xrfg/oWFPLOimML8HBxQmJ/DMyuKE9LCHij5aJtAEQla2k9tDWod9oE2AdE2gSIStLRPAEHqL/lom0ARCZoSQEC0TaCIBE0JIEDaJlBEgpTWL4FFRKRvSgAiIp5SAhAR8ZQSgIiIp5QAREQ8pQQgIuIpJQAREU8pAYiIeEoJQETEU0oAIiKeUgIQEfGUEoCIiKeUAEREPKUEICLiKSUAERFPBZIAnHOfcc7tdc5ddM6VBBGDiIjvgnoC2AOsAN4K6P4iIt4LZEcwM6sDcM4FcXsREQGcmQV3c+eqgG+bWU0/56wCVgEUFBQsrqys7PeaLS0t5ObmDmWYKcPnsoPf5fe57OB3+WMpe0VFRa2ZXdHdHrcnAOfcvwETevnVOjP7WazXMbPngOcASkpKrLy8vN/zq6qqGOicdOVz2cHv8vtcdvC7/NdS9rglADO7K17XFhGRa6dhoCIingpqGOgDzrkGIAT8wjm3MYg4RER8FtQooJeBl4O4t4iIdFEXkIiIp5QAREQ8pQQgIuIpJQAREU8pAYiIeEoJQETEU0oAIiKeUgIQEfGUEoCIiKeUAEREPKUEICLiKSUAERFPKQGIiHhKCUBExFNKACIinlICEBHxlBKAiIinlABERDylBCAi4qlA9gR2zq0H/j1wATgIPGJmTUHEcjVe2dHI+o31HGtqZVJ+DmuWz+H+hYVBhyUiMihBPQG8DtxiZrcCvwfWBhTHoL2yo5G1G3bT2NSKAY1NrazdsJtXdjQGHZqIyKAEkgDM7Fdm1tH942ZgchBxXI31G+tpbe+MOtba3sn6jfUBRSQicnWcmQUbgHM/B/6vmf1TH79fBawCKCgoWFxZWdnv9VpaWsjNzR3yOHvsbjzd5++KC/Pidt9YxLvsyc7n8vtcdvC7/LGUvaKiotbMSi4/HrcE4Jz7N2BCL79aZ2Y/6z5nHVACrLAYAikpKbGampp+z6mqqqK8vHzwAcdo2bNv0NjUesXxwvwcfvfUHXG7byziXfZk53P5fS47+F3+WMrunOs1AcTtJbCZ3TVAQCuBe4A7Y6n8k8Wa5XNYu2F3VDdQTlYma5bPCTAqEZHBC2oU0N3AnwG3m9m5IGK4Wj2jfTQKSERSXSAJAPgukA287pwD2Gxm/zGgWAbt/oWFqvBFJOUFkgDMbGYQ9xURkY9pJrCIiKeUAEREPKUEICLiKSUAERFPBT4TeDCccyeAwwOcNg44mYBwkpHPZQe/y+9z2cHv8sdS9pvMbPzlB1MqAcTCOVfT24w3H/hcdvC7/D6XHfwu/7WUXV1AIiKeUgIQEfFUOiaA54IOIEA+lx38Lr/PZQe/y3/VZU+7dwAiIhKbdHwCEBGRGCgBiIh4Km0SgHPubudcvXPugHPuqaDjSSTn3BTn3JvOuTrn3F7n3BNBx5RozrlM59wO59xrQceSaM65fOfcT51z73T/HQgFHVOiOOf+U/ff+T3OuR8550YEHVM8OedecM596Jzbc8mx651zrzvn9nf/+7pYr5cWCcA5lwn8PfBHwM3A551zNwcbVUJ1AN8ys7lAGfB1z8oP8ARQF3QQAfk74F/NrAiYjyf/HZxzhcDjQImZ3QJkAg8FG1Xc/QC4+7JjTwG/NrNZwK+7f45JWiQAYClwwMzeNbMLQCVwX8AxJYyZHTez7d1/bqarAvBmwwLn3GTgT4DvBx1LojnnxgCfAJ4HMLMLZtYUaFCJNQzIcc4NA0YCxwKOJ67M7C3gD5cdvg/4YfeffwjcH+v10iUBFAJHL/m5AY8qwEs556YCC4EtAYeSSP8DeBK4GHAcQZgOnAD+T3cX2Pedc6OCDioRzKwR+BvgCHAcOG1mvwo2qkAUmNlx6GoMAjfE+sF0SQCul2PejW91zuUCLwHfMLMzQceTCM65e4APzaw26FgCMgxYBPwvM1sInGUQXQCprLuv+z5gGjAJGOWc+2KwUaWWdEkADcCUS36eTJo/Cl7OOZdFV+X/opltCDqeBFoG3Ouce4+urr87nHP/FGxICdUANJhZzxPfT+lKCD64CzhkZifMrB3YAIQDjikIHzjnJgJ0//vDWD+YLglgGzDLOTfNOTecrhdBrwYcU8K4ro2VnwfqzOxvg44nkcxsrZlNNrOpdH3vb5iZN61AM3sfOOqcm9N96E5gX4AhJdIRoMw5N7L7/4E78eQF+GVeBVZ2/3kl8LNYPxjUpvBDysw6nHOrgY10jQR4wcz2BhxWIi0DvgTsds7t7D7252b2y+BCkgR6DHixu/HzLvBIwPEkhJltcc79FNhO10i4HaT5khDOuR8B5cA451wD8JfAs8CPnXNfpSspfibm62kpCBERP6VLF5CIiAySEoCIiKeUAEREPKUEICLiKSUAERFPKQGIiHhKCUBExFNKACLXwDm3xDn3tnNuhHNuVPfa9LcEHZdILDQRTOQaOef+GzACyKFrXZ5nAg5JJCZKACLXqHsJhm3AeSBsZp0BhyQSE3UBiVy764FcYDRdTwIiKUFPACLXyDn3Kl1LUU8DJprZ6oBDEolJWqwGKhIU59yXgQ4z++fuvak3OefuMLM3go5NZCB6AhAR8ZTeAYiIeEoJQETEU0oAIiKeUgIQEfGUEoCIiKeUAEREPKUEICLiqf8PpDZ/xDNCHSMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "q = .8\n",
    "quantile = QuantileRegressor(quantile=q, fit_intercept=True, alpha=0).fit(X,y)\n",
    "w_quantile = [quantile.coef_[0], quantile.intercept_]\n",
    "\n",
    "# plot the fit\n",
    "plotdata()\n",
    "plotline(w_quantile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.27211613802894624"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(y - quantile.predict(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.3906153070161097"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.median(y - quantile.predict(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.288787692705806e-09"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.quantile(y - quantile.predict(X), q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Huber regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjJklEQVR4nO3dfXRV13nn8e8jISGBANmABQhs3mSBeRW6N2kCsVFiB8+YIbbsjF+S2I3HC49Tx3Ha0oY6bddM65pZdLoms9Jpx409SVfcqFkJpo6TlLjFctqkaaQrAcKAeLENljAYbCRLIKG3PX9c6ZoLerlCujr3av8+a7FA5x6d82xL3s85+zxnb3POISIi/skIOgAREQmGEoCIiKeUAEREPKUEICLiKSUAERFPTQg6gOGYMWOGmz9//qD7nD9/nsmTJ49NQCnG57aD3+33ue3gd/sTaXskEjnrnJt5+fa0SgDz58+nurp60H0qKytZv3792ASUYnxuO/jdfp/bDn63P5G2m9nx/rZrCEhExFNKACIinlICEBHxVFo9A+hPZ2cnDQ0NtLe3AzBt2jQOHjwYcFSjKycnh7lz55KVlRV0KCIyjqR9AmhoaGDKlCnMnz8fM6OlpYUpU6YEHdaocc7x3nvv0dDQwIIFC4IOR0TGkUATgJm9BbQA3UCXcy403GO0t7fHOv/xyMyYPn06Z86cCToUEQnAztpGtu+q52RTG3Pyc9myoZg7SwpH5dipcAdQ5pw7O5IDjNfOv894b5+I9G9nbSNbd9TR1tkNQGNTG1t31AGMShLQQ2ARkRS1fVd9rPPv09bZzfZd9aNyfAtyPQAzexM4Bzjg/zrnnu1nn83AZoCCgoLSioqKuM+nTZvG4sWLY193d3eTmZmZzLAH9Wd/9mfk5eXxxBNP9Pv5yy+/zOLFi1myZMmwjnv06FGam5sH3ae1tZW8vLxhHXc88bn9Prcdxm/76xoH/n9+ReE0ILG2l5WVRfobYg96CGitc+6kmV0HvGJmh5xzP790h96k8CxAKBRyl7/xdvDgwbiHvkE/BJ44cSITJ04cMIZdu3aRlZVFOBwe1nFzcnIoKSkZdB+f34YEv9vvc9th/Lb/qW27aWxqu2J7YX4uX/7cemBkbQ90CMg5d7L373eBF4GPJPucO2sbWbttNwu+9mPWbtvNztrGER/z6aefpri4mFtvvZX6+uit2d/8zd8QDodZtWoVd999NxcuXOCXv/wlL730Elu2bGH16tUcO3as3/1ERAC2bCgmNyt+RCM3K5MtG4pH5fiBJQAzm2xmU/r+DXwa2J/Mc/Y9UGlsasPx4QOVkSSBSCRCRUUFtbW17Nixg6qqKgDKy8upqqpi7969LF26lOeee46Pf/zjbNq0ie3bt7Nnzx4WLVrU734iIhB90PtM+QoK83Mxolf+z5SvGBdVQAXAi70VLhOAv3PO/WMyTzjYA5Wr/Q/6L//yL9x1111MmjQJgE2bNgGwf/9+vv71r9PU1ERraysbNmzo9/sT3U9E/HRnSeGodfiXCywBOOfeAFaN5TlP9jOWNtj2RPVXpvmbv/mb7Ny5k1WrVvHtb3+bysrKfr830f1EREabV2Wgc/Jzh7U9ETfffDMvvvgibW1ttLS08KMf/QiIPoyePXs2nZ2dvPDCC7H9p0yZQktLS+zrgfYTEUk2rxJAMh6orFmzhnvvvZfVq1dz991384lPfAKAP/mTP+GjH/0ot912W1zJ53333cf27dspKSnh2LFjA+4nIpJsQZeBjqm+cbTRfq36qaee4qmnnrpi+2OPPXbFtrVr13LgwIG4ffrbT0Qk2bxKAJDcByoiIunEqyEgERH50LhIAEFOZzEWxnv7RCQYaZ8AcnJyeO+998ZtJ9m3HkBOTk7QoYjIOJP2zwDmzp1LQ0NDbL789vb2cddZ9q0IJiIymtI+AWRlZcWtlFVZWTnkpGkiIjIOhoBEROTqKAGIiHhKCUBExFNKACIinlICEBHxlBKAiIinlABERDylBCAi4iklABERTykBiIh4KvAEYGaZZlZrZi8HHYuIiE8CTwDAV4CDQQchIuKbQBOAmc0F7gC+FWQcIiI+siDn0TezHwDPAFOA33XObexnn83AZoCCgoLSioqKQY/Z2tpKXl5eEqJNfT63Hfxuv89tB7/bn0jby8rKIs650OXbA5sO2sw2Au865yJmtn6g/ZxzzwLPAoRCIbd+/YC7AtHpoIfaZ7zyue3gd/t9bjv43f6RtD3IIaC1wCYzewuoAD5pZt8NMB4REa8ElgCcc1udc3Odc/OB+4DdzrnPBxWPiIhvUqEKSEREApASS0I65yqByoDDEBHxiu4AREQ8pQQgIuIpJQAREU8pAYiIeEoJQETEU0oAIiKeUgIQEfGUEoCIiKeUAEREPKUEICLiKSUAERFPKQGIiHhKCUBExFNKACIinlICEBHxlBKAiIinlABERDylBCAi4iklABERTykBiIh4KrAEYGY5ZvZrM9trZq+b2X8LKhYRER9NCPDcF4FPOudazSwL+Fcz+6lz7lcBxiQi4o3AEoBzzgGtvV9m9f5xQcUjIuIbi/bDAZ3cLBOIAIuBv3TO/X4/+2wGNgMUFBSUVlRUDHrM1tZW8vLykhBt6vO57eB3+31uO/jd/kTaXlZWFnHOha74wDkX+B8gH3gVWD7YfqWlpW4or7766pD7jFc+t905v9vvc9ud87v9ibQdqHb99KkpUQXknGsCKoHbg41ERMQfQVYBzTSz/N5/5wK3AoeCikdExDdBVgHNBr7T+xwgA/i+c+7lAOMREfFKkFVA+4CSoM4vIuK7lHgGICIiY08JQETEU0oAIiKeUgIQEfGUEoCIiKeUAEREPKUEICLiKSUAERFPKQGIiHhKCUBExFNKACIinlICEBHxlBKAiIingpwOWkQk7eysbWT7rnpONrUxJz+XLRuKubOkMOiwroruAEQkbe2sbWTttt3UNTazdttudtY2Jv18W3fU0djUhgMam9rYuqMu6edNFiUAEUlLl3bGMDad8fZd9bR1dsdta+vsZvuu+qSdM5mUAEQkLQXRGZ/sTTaJbk91egYgImkpiM54Tn5u7I7j8u2ppK2tjT179lBVVUV1dfWA+ykBiEhaCqIz3rKhmK076uLuPHKzMtmyoThp5xxKV1cXtbW1VFVVxf7s37+f7u5ojLNmzRrwewNLAGY2D/hbYBbQAzzrnPtGUPGISHoJojPuq/YJqgqou7ubw4cPx3X2NTU1dHZ2AnDNNdcQDofZuHEjoVCIcDhMYWEhZtbv8YK8A+gCfsc5V2NmU4CImb3inDsQYEwikiYu7YyhhcIx6ozvLClM2jkuLTGdPS2HB5fncm1bQ6yzj0QitLa2AjB58mRKS0u56667KC8vJxQKsXDhwgE7+/4ElgCcc+8A7/T+u8XMDgKFgBKAiCSkrzOurKzky59bH3Q4I/L8zyL80bdeorWhnounjnDi1FH+re0DALKzs1m9ejUPPfQQ4XCYcDhMcXExmZmZVFZWsn79+qs6pznnRrEJV8fM5gM/B5Y75z647LPNwGaAgoKC0oqKikGP1draSl5eXpIiTW0+tx38bX9TWydd7Rd45wJkZ2ZQMC2H/NysoMMaU6n6s29q6+R0czsd3T1xP5vm5mbq6+vj/pw9exaAjIwMZs+9nnkLFnPDwiIWFd3I+tAysrL6/5km0vaysrKIcy50+fbAE4CZ5QGvAU8753YMtm8oFHKDPdEGRpQN053PbQc/299XC/+lJRf5n3XRG/rcrEyeKV+Rtm+nXo1U/Nn3/WzOt7bQcfoYF985Qve7R5n0wXFONxyP7VdcXEwoFOKnpyaRNetGsgsWkJGVE/vcgDe33THgeRJpu5n1mwACrQIysyzgh8ALQ3X+InKlwWrhfUoAqaK9vZ29e/dSVVXFf3/+JZpOHKLzvQYgeqGdOfU6sq9fwhfKH6DmwnRa8uYxvWAG92wo5s1d9WNe1RRkFZABzwEHnXN/EVQcIulsvL2YlE46Ozs5cOBAXEVOXV0dXV1dAGRMzmfirCImLb2ZibMWkz2riMzJ+QBUZWXSltuN8eEbzHeXFvLDSOOYVjUFeQewFvgCUGdme3q3/YFz7ifBhSSSXtLlxaR019PTEyu/rK6upqqqitraWtrb2wHIz88nFAqxZcsWwuEwoVCIe797mJPN7VccK9Os37u2Vw+d4ZnyFWNaYhpkFdC/Eh3eEpGr1FcLH62qjgr6xaR055zj+PHjsav66upqIpEIH3wQrU+ZNGkSa9as4bHHHotV5CxatOiK8svfuz2j3/cULu/8+5xsaktqiWl/9CawSBrr6yxO19dgkPbTEwfh1KlTccM41dXVsYqc7OxsVq1axec+97lYZ7906VIyMzOHPO5AL41tD2CsfyBKACJp7s6SQiqbj/DmtvVBh5Ly3n//fSKRSFyH39gYnT00IyODZcuWsWnTpthbtCtWrGDixIlXfb6BruhTZToJJQARGZdaW1upqamJjdlXVVVx7Nix2OdFRUXccsstsTH7kpISJk+enPS4gp5O4lJKACKS9jo6Ovj1r38d95D24MGD9PT0ADBj1hx6pi8i/5Z1zFm8nD948D/yuZtvCizesR7rH4gSgIikla6urrjyy+rqavbu3Rsrv7zuuusIh8Pcc889hMNhzkws5H+8dio25HIe+NNXTjB5yrSU6ISDpAQg4pl0WtO2p6eHo0ePxo3Z19bW0tYWfYg6bdo0QqEQn/3sZ7n77rsJh8PMmzcvriJn7bbdelluAEoAIh7pm56gr0PsewkJCLwzdM5x4sSJuDH7SCRCc3MzALm5uaxZs4ZHH300Nm6/ePFiMjIyBp0OQS/LDWzIBGBmjxOdquHcGMQjIkmUSlNHnD59Om7MvqqqijNnzgCQlZXFypUruf/++wmHw7TkXc/fH3E0tnTg8nO5ZWkxN96YWLx6WW5gidwBzAKqzKwGeB7Y5YKeQU5ErkpQV8NNTU2xjr7v77fffhuIll/edNNN3HHHHbFa+5UrV8bKL3fWNrJ9BHctqbiKV6oYMgE4575uZn8IfBr4IvBNM/s+8Jxz7tjg3y0igxnr8fixuBo+f/78FUsUHj16NPb54sWLWbduXazWvqSkZNDpjEd615JKZZepJqFnAM45Z2angFNE3zm/BvhB7wpev5fMAEVGW6o8BA1iPH60r4YvXrzIvn374oZxDhw4ECu/nDt3LuFwmIcffphQKEQoFOKaa64Z1jkGujtpbGpj7bbdbNlQTP4Qxxhu2WWq/I4kWyLPAJ4AHgLOAt8CtjjnOs0sAzgCKAFI2kilh6BBjMeP5Gq4u7s7Vn7Z1+Hv27ePjo4OAGbMmEE4HKa8vDz2kHawBckTNdBdC3z483vm40NPzZCoVPodSbZE7gBmAOXOueOXbnTO9ZjZxuSEJZIcqfQQNKjx+ESuhp1zsfLLvs6+pqaGCxcuADB16lRKS0t58sknY+P2119//bDWo01U2ZKZvPCrEwz04LGts5vTzZ2jdr5U+h1JtkSeAfzRIJ8dHN1wRJIrlUoCg65O6RvmaDx3gRnWyqdmtJDx3htUV1dTXV1NU1MTEC2/LCkp4ZFHHol19kVFRWRkZAx63NEYPtlZ28gPI40Ddv59Orp7rur4/Uml35Fk03sA4pWgO91LBVWdcubMGf53xU/56x+8wvnG6ALkx883EQEyJ0xg1cqV3HvvvbGHtMuWLWPChMS6itEePunvarw/2Zn9J6OrkUq/I8mmBCBeSaWSwLGoTmlubo5d0fc9pD1x4kTvp0bW9HnkLihl4uwismcVMb/oJv7tD2+/6vON9vBJIlfduVmZFEzLHvaxB5JKvyPJpgQgXkm1ksDRnBTswoULsfLLvg7/8OHDsc8XLlzIxz72MZ544gn+vKaLrOsWkjFxUtwxTp0f+mp7MKM9fDLQ1XimGT3OxX5++c1Hrur4/Um135FkUgIQ76TKTIwj0dHRQV1dXeyqvrKykuPHj9PdHe3ACwsLCYVCPPjgg7GKnGuvvTb2/Tu27U7KMMdoD58MdDX+TPmKuJ9hZeXoJQAYH78jiVACEElBlz5InT01m/tuzGRq64lYh793795Y+eX06dNZtGgRDzzwQOwh7ezZswc9frKGOUb7uD5djQch0ARgZs8DG4F3nXPLg4xFJBU45/jrl37Bn37nZVoa6ul45wjHTx/j3zqji4tPmTKF0tJSvvKVr8Qe0s6fP5/XXnttwMnQ+pOsjjUZx/XlajwIQd8BfBv4JvC3AcchMuacczQ2NsaN2VdXV3PuXHTeRZuQTdZ1C8hbeRvZs4qYd+MKIn/+4IDll8OVrI5VHXb6CDQBOOd+bmbzg4xBxp9UfY3/7NmzcYuYVFVVcerUKQAmTJjAihUruOeee9hxYiLZs28ka8b1WOaH/4s2wah1/iIAFvTEnr0J4OWBhoDMbDOwGaCgoKC0oqJi0OO1trYOOrHUeOZz2yHa/q7MiTSea6Pnkt/rDDMKr8klPzdrzGI5f/48hw8f5tChQ9TX11NfXx/r7M2M66+/nuLiYoqLi1myZAmLFi2KzX5Zf6ql3xebsjMzKJ41pd/z6Wfvb/sTaXtZWVnEORe6fHvQQ0BDcs49CzwLEAqF3FDjnIMtDDHe+dD2wa7uKysreepXPTQ29TcvTAeF+ZlJuRtoa2tjz549cbNf1tfXxz5fsGABN998c6waZ82aNUydOnXA4zVd9jIVfFj5sn6A2H342Q/G5/aPpO0pnwBE+iTylulg9eajMalXZ2cndXV1cS9W7d+/P1Z+OXv2bMLhMJ///OcJh8OUlpYyY8aMYZ1DlS8yVpQAJG0k8pbpYDNH9rf/YLq7u6mvr48bs9+zZw8XL14E4NprryUUCrFx48ZY+eWcOXNG0MIP6UGqjIWgy0C/B6wHZphZA/DHzrnngoxJUlcib5n2V4eeyHGcc7z55pv8n+/v4ns/eZX33jpI5+lj9HRE983Ly6O0tJTHH3881tkvWLBg0NkvU/VhtEifoKuA7g/y/JJeEnnL9NLhk4HuBObk53Ly5Mm4Mfvq6mref//96A6ZWWRft5DJyz9FXuGNbH1oI1+682YyMxOfc96nOeUlfWkISNJGom+Z9g2f9HXCrR+co+OdI1w8dYTu00dpOfcmhVtPA5CZmcny5cspLy9n99k8LkybT/bMG7DMDyuGKo44vjyMzh/8mlNe0pcSgKSNRB6OtrS0UFNTE7uyP/eLX3G68UTs88L5iyjbcFusImf16tVMmhSdEG3B137MxH7OezUTmSUyXKUhIgmaEoCklUsfjra3t7Nnzx6++c0Xqaqq4rXXXuPEiRP0vdtyww038ImPfYRw+LcIh8OsWbOGadOmDXjs0ZzIbKhjaYhIUoESgKSFzs5OXn/99bjyy7q6Orq6ugAoKChg4cKFPPzww7Gr+5kzZw7rHKM5kdlQx9IQkaQCJQBJOT09PRw+fDiu/LK2tpb29uiEaPn5+YRCIbZs2RKryCksLBz2hGiXG836+6GO5dOyg5K6lACGoHHa5HLOcfz48biKnEgkQktLCwCTJk2itLSUxx57LNbZL1q0KCmLj8Po1t8Pdiyflh2U1KUEMAiN046+U6dOXVF+efbsWQCys7NZtWpV7C3acDjM0qVLh1V+maouv5AoWzKTH0YavVh2UFKXEsAgNE47Mu+//z6RSCSuw29sbASis1ouW7aMTZs2xea1X7FiRWxCtPGkvwuJH0Yaubu0kFcPndHdpQRGCWAQGqdNXGtra2w92r4/x44di31eVFTELbfcEntAW1JSwuTJkwOMeOwMdCHx6qEz/OJrnwwoKhElgEFpnLZ/Fy9eZO/evXEPaQ8ePEhPT3QK43nz5hEOh3nkkUdiE6Ll5+cHG3SAdCEhqUoJYBDJWjc1nXR1dXHgwIG4Mft9+/bR2dkJwMyZMwmHw9xzzz2xq/uCgoKAo04tupCQVKUEMAjfpuXt6enh6NGjcZ19TU0NbW3RzmvatGmEQiF++7d/O/aQdt68eUmryBkvdCEhqUoJYAjjdVpe5xwnTpyIe7EqEonQ3NwMQG5uLmvWrOHRRx+NPaRdvHixliS8Cr5dSEj6UALwxOnTp+M6+6qqKs6cOQNAVlYWK1eu5P77748N49x0001MmKBfj9EyXi8kJL3p//BxqKmpKdbZ9/399ttvA9Hyy6VLl3LHHXfEhnFWrlw5LssvRWRwSgBp7vz587Hyy5dffpnNmzdz5MiR2OeLFy9m3bp1sWGckpISbxfPFpF4SgBp5OLFi9TV1cUN4xw4cCBWfjlz5kzWrVvHF7/4RUKhEKFQiGuuuSbgqEUkVSkBpKju7m4OHDgQN26/b98+Ojo6AJgxYwbhcJjy8vLY1f2hQ4dGNBmaiPhFCSAFOOdi5Zd9HX5NTQ0XLlwAYOrUqZSWlvLkk0/GHtLecMMNV5RfHjp0KIjwRSRNBb0o/O3AN4BM4FvOuW1BxjMWnHM0NDTE1dpXV1fT1NQEQE5ODiUlJbG3aMPhMEVFRSq/FJFRF1gCMLNM4C+B24AGoMrMXnLOHQgqpmQ4c+bMFbNfnj4dXY92woQJrFy5knvvvTc2jLNs2TKVX4rImAiyp/kIcNQ59waAmVUAnwHSNgE0NzdfMfvliRPR9WjNjKVLl3L77bfHhnFWrVpFTk5OwFGLiK+sb/3UMT+x2T3A7c65R3q//gLwUefc45fttxnYDFBQUFBaUVEx6HFbW1vjyhyb2jo53dxOR3cP2ZkZFEzLIT83a8Txt7e3c+TIEerr62N/+mrtAebMmUNxcTFLliyhuLiYoqKi2OLjyXJ5233jc/t9bjv43f5E2l5WVhZxzoUu3x7kHUB/E8hckY2cc88CzwKEQiE3VJVLZWVlrBJmZ20jW/+5jrbODCA6hp6b1c0z5TcN663Mjo6OWPll30Pa119/ne7u6NwuhYWFhEIhHn300djV/bXXXhtbBOSn77Qxpy2HLRuKkvo26KVt95HP7fe57eB3+0fS9iATQAMw75Kv5wInR/MEV7OgS3d3N4cOHYobs9+zZ0+s/HL69OmEQiE2bdoU6+znzJlzxXG0mpiIpLogE0AVUGRmC4BG4D7ggdE8wVDzsDvneOONN+LG7Gtqajh//jwAU6ZMobS0lCeeeCJWkTN//vyEZr9MJPlovWERCVJgCcA512VmjwO7iJaBPu+ce300z3HpPOzOObpb36PjnSNkn3uTT3/6G1RXV3Pu3DkAJk6cSElJCQ8//HCsIqe4uPiqyy+HSj66QxCRoAVab+ic+wnwk2Qc++zZs3wq7yR/9Y8/43xjPR3vHKH7fLSzz8jMZOaKFXGLmCxfvpysrJE/HO4z1CIgWm9YRII2LgrOP/jgAyKRCNXV1fz4xz/m+PHjvPXWW0C0/DJnxjxy5q/muoXLePTu23jiP99Gbm5yV2MaahEQLRMoIkFLuwTQ1tbGnj174ipy6uvr6StnnT17NuvWreNLX/oS4XCYNWvWMHXq1DGPc6hFQLRMoIgELa0SwIEDB5g6dSpdXV0AzJo1i3A4zAMPPBAbytm/f3/KlIMNtgiIlgkUkaClVQKYMGECX/3qV2MVOYWF6TtWrmUCRSRoaZUAbrzxRp5++umgwxg1WiZQRIKkKSZFRDylBCAi4iklABERTykBiIh4SglARMRTSgAiIp5SAhAR8ZQSgIiIp5QAREQ8pQQgIuIpJQAREU8pAYiIeEoJQETEU0oAIiKeUgIQEfFUIAnAzD5rZq+bWY+ZhYKIQUTEd0HdAewHyoGfB3R+ERHvBbIimHPuIICZBXF6EREBzDkX3MnNKoHfdc5VD7LPZmAzQEFBQWlFRcWgx2xtbSUvL280w0wbPrcd/G6/z20Hv9ufSNvLysoizrkrhtuTdgdgZv8EzOrno6ecc/+Q6HGcc88CzwKEQiG3fv36QfevrKxkqH3GK5/bDn633+e2g9/tH0nbk5YAnHO3JuvYIiIycioDFRHxVFBloHeZWQPwMeDHZrYriDhERHwWVBXQi8CLQZxbRESiNAQkIuIpJQAREU8pAYiIeEoJQETEU0oAIiKeUgIQEfGUEoCIiKeUAEREPKUEICLiKSUAERFPKQGIiHhKCUBExFNKACIinlICEBHxlBKAiIinlABERDylBCAi4iklABERTykBiIh4KpA1gc1sO/CfgA7gGPBF51xTELFcjZ21jWzfVc/Jpjbm5OeyZUMxd5YUBh2WiMiwBHUH8Aqw3Dm3EjgMbA0ojmHbWdvI1h11NDa14YDGpja27qhjZ21j0KGJiAxLIAnAOfcz51xX75e/AuYGEcfV2L6rnrbO7rhtbZ3dbN9VH1BEIiJXx5xzwQZg9iPg751z3x3g883AZoCCgoLSioqKQY/X2tpKXl7eqMfZp66xecDPVhROS9p5E5Hstqc6n9vvc9vB7/Yn0vaysrKIcy50+fakJQAz+ydgVj8fPeWc+4fefZ4CQkC5SyCQUCjkqqurB92nsrKS9evXDz/gBK3dtpvGprYrthfm5/KLr30yaedNRLLbnup8br/PbQe/259I282s3wSQtIfAzrlbhwjoIWAj8KlEOv9UsWVDMVt31MUNA+VmZbJlQ3GAUYmIDF9QVUC3A78P3OKcuxBEDFerr9pHVUAiku4CSQDAN4GJwCtmBvAr59x/DSiWYbuzpFAdvoikvUASgHNucRDnFRGRD+lNYBERTykBiIh4SglARMRTSgAiIp4K/E3g4TCzM8DxIXabAZwdg3BSkc9tB7/b73Pbwe/2J9L2G5xzMy/fmFYJIBFmVt3fG28+8Lnt4Hf7fW47+N3+kbRdQ0AiIp5SAhAR8dR4TADPBh1AgHxuO/jdfp/bDn63/6rbPu6eAYiISGLG4x2AiIgkQAlARMRT4yYBmNntZlZvZkfN7GtBxzOWzGyemb1qZgfN7HUz+0rQMY01M8s0s1ozeznoWMaameWb2Q/M7FDv78DHgo5prJjZV3t/5/eb2ffMLCfomJLJzJ43s3fNbP8l2641s1fM7Ejv39ckerxxkQDMLBP4S+A/ADcB95vZTcFGNaa6gN9xzi0FfgP4Lc/aD/AV4GDQQQTkG8A/OueWAKvw5L+DmRUCTwAh59xyIBO4L9ioku7bwO2Xbfsa8M/OuSLgn3u/Tsi4SADAR4Cjzrk3nHMdQAXwmYBjGjPOuXecczW9/24h2gF4s2CBmc0F7gC+FXQsY83MpgI3A88BOOc6nHNNgQY1tiYAuWY2AZgEnAw4nqRyzv0ceP+yzZ8BvtP77+8AdyZ6vPGSAAqBty/5ugGPOsBLmdl8oAT494BDGUv/C/g9oCfgOIKwEDgD/L/eIbBvmdnkoIMaC865RuDPgRPAO0Czc+5nwUYViALn3DsQvRgErkv0G8dLArB+tnlX32pmecAPgSedcx8EHc9YMLONwLvOuUjQsQRkArAG+CvnXAlwnmEMAaSz3rHuzwALgDnAZDP7fLBRpZfxkgAagHmXfD2XcX4reDkzyyLa+b/gnNsRdDxjaC2wyczeIjr090kz+26wIY2pBqDBOdd3x/cDognBB7cCbzrnzjjnOoEdwMcDjikIp81sNkDv3+8m+o3jJQFUAUVmtsDMsok+CHop4JjGjEUXVn4OOOic+4ug4xlLzrmtzrm5zrn5RH/uu51z3lwFOudOAW+bWXHvpk8BBwIMaSydAH7DzCb1/j/wKTx5AH6Zl4CHev/9EPAPiX5jUIvCjyrnXJeZPQ7sIloJ8Lxz7vWAwxpLa4EvAHVmtqd32x84534SXEgyhr4MvNB78fMG8MWA4xkTzrl/N7MfADVEK+FqGedTQpjZ94D1wAwzawD+GNgGfN/M/gvRpPjZhI+nqSBERPw0XoaARERkmJQAREQ8pQQgIuIpJQAREU8pAYiIeEoJQETEU0oAIiKeUgIQGQEzC5vZPjPLMbPJvXPTLw86LpFE6EUwkREysz8FcoBcovPyPBNwSCIJUQIQGaHeKRiqgHbg48657oBDEkmIhoBERu5aIA+YQvROQCQt6A5AZITM7CWiU1EvAGY75x4POCSRhIyL2UBFgmJmDwJdzrm/612b+pdm9knn3O6gYxMZiu4AREQ8pWcAIiKeUgIQEfGUEoCIiKeUAEREPKUEICLiKSUAERFPKQGIiHjq/wPRQz4DA6Pw4QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "huber = HuberRegressor(fit_intercept=True, alpha=0).fit(X,y)\n",
    "w_huber = [huber.coef_[0], huber.intercept_]\n",
    "\n",
    "# plot the fit\n",
    "plotdata()\n",
    "plotline(w_huber)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.06140586728499734"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(y - huber.predict(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.03137143891216554"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.median(y - huber.predict(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'hcat' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/85/fz3xj3ln08n9hgtdxwtq9q7r0000gn/T/ipykernel_41736/143207705.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mhcat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mw_quad\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mw_l1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mw_huber\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mw_quantile\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'hcat' is not defined"
     ]
    }
   ],
   "source": [
    "hcat(w_quad, w_l1, w_huber, w_quantile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot quantiles of distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = .1\n",
    "l = QuantileLoss(quantile=q)\n",
    "w_nocov = proxgrad(1/n*l, ZeroReg(), ones(n,1), y, stepsize=.9, maxiters=50)\n",
    "\n",
    "scatter(y, zeros(n), axis = nothing, shape=:x, label=\"data\")\n",
    "scatter!([w_nocov], [0], color = :red, shape=:square, label=\"$(q)th quantile\") \n",
    "ylims!((-1,1))\n",
    "m,M = minimum(y), maximum(y)\n",
    "xlims!((minimum(y)-(M-m)/10, maximum(y)+(M-m)/10))\n",
    "# arrow(0, 0, maximum(y)+(M-m)/20, 0, head_width=0.05, head_length=0.1, fc=\"k\", ec=\"k\")\n",
    "# arrow(0, 0, minimum(y)-(M-m)/20, 0, head_width=0.05, head_length=0.1, fc=\"k\", ec=\"k\")\n",
    "# savefig(\"quantile_$(q).pdf\")"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}